Observations of meteors in the Earth's atmosphere: Reducing data from dedicated double-station wide-angle cameras
A. Margonis A. Christou J. Oberst
 


Observations of meteors in the Earth's atmosphere



















Department of Geodesy and Geoinformation Science, Technische Universität Berlin, Strasse des 17. Juni 135, Berlin, Germany

anastasios.margonis@tu-berlin.de
Armagh Observatory, College Hill, Armagh BT61 9DG, UK

aac@arm.ac.uk
Germany Aerospace Center, Institute of Planetary Research, Rutherfordstr. 2, 12489 Berlin, Germany

juergen.oberst@dlr.de
          





 
  Meteoroids entering the Earth's atmosphere can be observed as meteors, thereby providing useful information on their formation and hence on their parent bodies. We developed a data reduction software package for double station meteor data from the SPOSH camera, which includes event detection, image geometric and radiometric calibration, radiant and speed estimates, trajectory and orbit determination, and meteor light curve recovery. The software package is designed to fully utilise the high photometric quality of SPOSH images. This will facilitate the detection of meteor streams and studies of their trajectories. We have run simulations to assess the performance of the software by estimating the radiants, speeds, and magnitudes of synthetic meteors and comparing them with the a priori values. The estimated uncertainties in radiant location had a zero mean with a median deviation between 0.03^∘ and 0.11^∘ for the right ascension and 0.02^∘ and 0.07^∘ for the declination. The estimated uncertainties for the speeds had a median deviation between 0.40 and 0.45 km s^-1. The brightness of synthetic meteors was estimated to within +0.01m. We have applied the software package to 177 real meteors acquired by the SPOSH camera. The median propagated uncertainties in geocentric right ascension and declination were found to be of 0.64^∘ and 0.29^∘, while the median propagated error in geocentric speed was 1.21 km s^-1.

   

   Observations of meteors in the Earth's atmosphere:
 
Reducing data from dedicated double-station wide-angle cameras
    A. Margonis1, A. Christou2 and J. Oberst1,3
    Received 2018; accepted 2018
===================================================================================================================





§ INTRODUCTION



Observations of meteors in the Earth's atmosphere shed light on the properties of the population of meteoroids intercepting the orbit of our planet. The study of the temporal and spatial distribution of meteors requires sensitive optical systems that are able to monitor the night sky. Double station observations (i.e. observations of two cameras from different positions) are required to determine the trajectories and orbit parameters of the meteors. 

While algorithms for meteor data reduction are well established in the literature <cit.>, every camera may require an analysis system to account for the specific capabilities of the camera.

For observations of meteors in recent years, our team has used the Smart Panoramic Optical Sensor Head (SPOSH) camera <cit.>. The instrument features a highly sensitive CCD chip that delivers images of high photometric quality. With the wide-angle lens, the camera easily captures several hundreds of stars in one image, which requires sophisticated geometric calibration procedures.

To process the data from SPOSH, we developed a comprehensive software package, which allows us to carry out camera calibration, meteor detection, meteor trajectory determination, meteor photometric modelling, and orbit reconstruction. In this paper, we describe this software and demonstrate and assess its performance on synthetic and actual meteor data. 



§ THE SPOSH CAMERA



The SPOSH camera was designed to image faint transient noctilucent phenomena, such as aurorae, electric discharges, meteors, or impact flashes on dark planetary hemispheres from an orbiting platform <cit.>. The camera is equipped with a highly sensitive back-illuminated 1024×1024 CCD chip and has a custom-made optical system of high light-gathering power with a wide field of view (FOV) of 120×120^∘. The SPOSH camera system is accompanied by a sophisticated digital processing unit (DPU) designed for real-time image processing and communication with a spacecraft. Owing to the all-sky coverage and excellent radiometric and geometric properties  of the camera, a large number of meteors can be obtained for reliable event statistics.  

For outdoor tests and meteor monitoring, the camera is typically mounted on a tripod pointed vertically up at the sky taking one image every 2 s. For the determination of the meteor velocity, a mechanically rotating shutter with a known frequency is mounted in front of the camera lens. The shutter consists of two blades and has a rotating frequency of 250 RPM resulting in an exposure time of 0.06 s for every shutter opening. Double-station observations have been carried out routinely providing a large dataset of meteor images.



§ DATA REDUCTION


The reduction of the meteor data is performed by different software modules. The calibration software SPOSHCalib is a stand-alone software for the geometric calibration of SPOSH images <cit.>. The trajectory determination module is based on the MOTS software <cit.>, which was initially modified to process SPOSH data <cit.>. All modules were developed anew within the scope of this study. The interaction between the different modules can be seen in Figure <ref>.





 §.§ Meteor detection



Unlike video cameras, where a meteor only spends a fraction of their trajectory in each frame, exposures longer than one second often capture the whole meteor (e.g. a Perseid) in one image. The meteor detection algorithm that we used is based on the Hough transform technique for extracting linear features within images. This method has been used to detect meteors in photographic image data by previous authors <cit.>. 

The algorithm that we developed first generates 8 bit difference images between three consecutive frames, removing the background and highlighting only short temporal variations. Possible non-relevant information depicted in the margins of the images (e.g. surrounding mountains and man-made structures), typical within large FOVs, are removed by applying a circular mask. Background noise and stellar scintillation are filtered out by first applying an empirical threshold and then a median filter to the image, thus reducing the overall computation time of the algorithm. Each line, represented by a combination of ρ and θ, passing through each of the remaining pixels contributes to the parameter space H(ρ, θ), known also as voting space, by adding the value A_xy=1, 


    H(θ,ρ) = ∑_x∑_yA_xyδ(ρ, [ρ'])
    ,


where ρ is the distance from the origin to the closest point on the straight line, and θ is the angle between the x axis and the line connecting the origin with that closest point. The square brackets [ ] indicate rounding to the nearest integer and the normal representation of a line is


    ρ'=xcosθ+ysinθ,
 
with

    δ(ρ, [ρ'])= 
     
       1,   ρ = ρ' 
    
       0,    otherwise.
 
The event detection algorithm is triggered each time a certain threshold value is exceeded. This value is compared against the maximum value found in the parameter space of each image and represents the number of pixels lying on a line in the image space. 

Several criteria are used to mitigate the effect of false detections. Slow moving objects, such as airplanes and satellites, appear with a characteristic negative-positive-negative pattern in the difference images (Fig. <ref>). This pattern is compared against a predefined signal by the user, simulating the path of an airplane projected on the image plane. In this way, events appearing in three consecutive images moving with low apparent angular velocities of 0.6^∘ > v_ang > 2.2^∘ are rejected as slow-moving objects. This condition also affects meteors appearing close to their radiant position and/or close to the horizon. For every event, its time of occurrence, the central position of the line, and its direction within the image are saved together with the object name (meteor, slow-moving object, or star) in a text file. 



A quality parameter q_md was introduced to determine the threshold value for the Hough transform. The value of the threshold should ideally detect all meteors in the image data when applied to a meteor detection algorithm. At the same time slow-moving objects and random noise patterns resembling lines should be filtered out. To select a suitable value for this parameter, we balance the number of false positives against the number of meteors the algorithm failed to detect (false negatives) and the processing time it takes for the algorithm to scan the images. The value is computed, after applying various weights to the observed quantities, as follows:


    q_md = p_1w_1-p_2w_2-p_3w_3/100w_1-8w_3
    ,

where p_1 is the percentage of the detected meteors, p_2 the percentage of false detections, p_3 the processing time in minutes, and w_1, w_2, and w_3 the respective weights. The quality parameter is scaled to values between 0 and 1 (Eq. <ref>). The maximum value (q_md=1) for a threshold is reached when all meteors are detected (p_1=100), with no false detections (p_2=0), within a user-defined processing time. 




We tested the performance of our algorithm using various threshold values and applying these to > 14,000 images corresponding to eight hours of data from two observing sites. The results were compared with meteors identified after visual inspection of the images. The highest value of the quality parameter for this dataset was found for a threshold of 23 with w_1:0.6, w_2:0.3, and w_3:0.1. The threshold value corresponds to the highest number of pixels lying on a line in a given image. Applying these parameters, 70% of the visually identified meteors were successfully detected by the algorithm (true positives), while 15% of the detected events were false detections (false positives). Figure <ref> shows the calculated quality parameter for our dataset. High values are computed from data with a relative high signal-to-noise ratio in terms of detected meteors and false detections. This performance of the algorithm can be achieved under favourable weather conditions.





 §.§ Astrometry




  §.§.§ Camera calibration



The geometric calibration of the camera is performed by the SposhCalib software in a semi-automatic process using standard stars in the SPOSH images <cit.>. Stars are ideal calibration targets owing to their high abundance in the images and the precise knowledge of their position at a given time. The SPOSH images may feature up to several thousand stars, which are on average equally distributed over the whole image except image corners. By comparing the actual stars in the image with their expected positions based on a priori information about pointing and interior camera parameters, these parameters can be updated in a least-squares fashion. This provides an accurate knowledge of the interior, i.e. focal length and geometric distortion, and the exterior orientation (pointing) parameters. The coordinates of the stars are taken from the Tycho-2 and Hipparcos star catalogues <cit.>.

The transformation equations between the image coordinate system (x,y) and the camera coordinate system (X_cam,Y_cam,Z_cam) are described applying an equidistant camera model <cit.>,


    X_cam = x_c/√(x_c^2+y_c^2)sin(√(x_c^2+y_c^2)), 
        Y_cam = y_c/√(x_c^2+y_c^2)sin(√(x_c^2+y_c^2)), 
        Z_cam = cos(√(x_c^2+y_c^2)) .


A high number of standard stars is achieved by performing initially a pre-calibration with the help of at least six reference stars selected by the user. The pre-calibration step provides approximate values for the unknown parameters. After this step, a global calibration is performed using all point sources identified as standard stars in the image. 

The SPOSH images show significant radial and non-symmetrical distortion, mathematically expressed as


    Δ r_rad = A_1r^2 + A_2r^4 + A_3r^6



    Δ x_tan = B_1(r^2+2x^2)+2B_2xy, 
       Δ y_tan = B_2(r^2+2x^2)+2B_1xy



    Δ x_aff = C_1x, 
       Δ y_aff = -C_1y



    Δ x_sh = C_2y, 
       Δ y_sh = C_2x


with 


    r = √(x^2+y^2) .


The equations above describe the radial (Eq. <ref>) and non-symmetrical distortions (Eq. <ref>) and the deviations of the image coordinate system from an orthogonal, uniformly scaled coordinate system (Eqs. <ref>, <ref>). The outer and inner orientation of the camera and the distortion parameters introduced by the lens are determined by fitting a 6^th-order polynomial function. These distortion terms are added directly to the pixel coordinates of the stars. 

The average residual error for the star positions after the calibration is usually less than 0.25 pixel or 1.68 and usually consistent over the whole image. The displacement Δx_ij and Δy_ij in image coordinates due to radial distortion are stored in two separate TIFF files. These files serve as look-up tables in the subsequent steps providing the undistorted position of each pixel.



  §.§.§ Meteor path on the image plane



The projection of a trajectory of a meteor on an image plane can be seen as a meteor trail. By extending the trajectory before and after the luminous path, a line can be defined on the image plane representing the projection of that extended path. Once the line is defined in both images, its radiant can be determined (Section <ref>).

In order to speed up the process of defining the meteor line, a threshold is applied to each raw meteor image. The threshold is defined at 2σ of the noise level. The use of a relative low threshold ensures that fainter pixels belonging to the meteor trail are considered in the computation of the line. The line parameters defined in the meteor detection procedure (Section <ref>) are used to remove unwanted features in each meteor image, considering the proximity of each pixel to the detected line, its distance to the pixel with the maximum votes in the Hough transform, and its intensity value (Fig. <ref>). 


The positions of the remaining pixels are corrected for radial distortion by retrieving pixel-offset values from the look-up tables generated in the calibration step (Section <ref>). Owing to the equidistant projection model used by the lens system of the camera, perspective distortions in the image are evident that deflect the path of objects moving along a great circle from a straight line to a curved line. To efficiently detect linear features in the image, pixel coordinates are converted from an equidistant to a gnomonic projection, where straight lines in space preserve their straightness when projected on the image plane. 



The line along which the meteor is moving is computed by applying a customised Hough transform. As an input, we use the corrected image coordinates (in sub-pixel accuracy) belonging to the meteor trail. Lines running diagonal to the meteor trail results into a higher value in voting space than those parallel to the meteor motion, since more pixels lie along the diagonal line  (Fig. <ref>). We handle this effect as follows: First we apply a Hough transform to determine the top 20 lines intersecting the highest number of pixels. Then we perform a weighted Hough transform considering the intensity values. Unlike the standard Hough transform method, which searches for the line with the maximum votes V in parameter space, we defined a ratio coefficient calculated as the sum of intensity values I with respect to the number of pixels that are √(2)/2 pixels apart from each line parameter combination. A distance of √(2)/2 pixel is needed to identify which pixels lie on the line since the line does not cross the pixel centre (defined at 0.5 pix),


    V_max(θ,ρ) = ∑_i=1^n I_i/n
    .


The best-fitting line is defined as the line with the highest ratio. Since the point spread function (PSF) of an imaging system spreads the light of point sources to neighbouring pixels, the light emitted by a meteor also spreads to pixels located perpendicular to its motion. In order to account for the signal within these pixels, we define a buffer zone of 10 pixels perpendicular to the best line computed.

Occasionally, residual features may be located along the buffer zone. As a result, these remaining pixels affect the determination of the meteor line. In order to remove these features, the consecutive pixel-to-pixel distances are determined revealing gaps between features. Distances higher than a threshold indicate different pixel entities, where entity is a feature consisting of at least two neighbouring pixels; for example, the meteor trail or meteor segment is such an entity. Assuming that the meteor entity has the maximum number of pixels, we remove all secondary features from the line zone. Finally, the meteor line is determined using weighted least squares. The line parameters (slope plus intercept) and the middle point of the meteor trail, defined as the median of the chosen pixel coordinates, are saved in a text file. 





  §.§.§ Transformation to the spatial trajectory of the meteor



From the estimated parameters of the meteor line, the underlying image points are generated in sub-pixel accuracy and transformed from the gnomonic back to an equatorial projection. The pixel coordinates x_c, y_c, are normalised using the parameters of the interior orientation of the camera, i.e. (Section <ref>)

    x_n    =  (x_c-x_p)p_x /f, 
    
     y_n    =  (y_c-y_p)p_y /f
    ,

where x_p, y_p are the intersection of the optical axis with the image plane (principal point), p_x, p_y the pixel size, and f the focal length of the camera system. The points are first projected to the camera coordinate system (|x⃗|=1) using the equidistant projection equations. The vectors are then transformed to the local (horizontal) coordinate system, 

    𝐱_hor = 
     [ x_hor; y_hor; z_hor ]
      = R_ωϕκ·[ x_cam; y_cam; z_cam ]
    ,


where R_ωϕκ is the 3d rotation matrix that relates the camera to the local coordinate system. Finally, the pointing vectors are transformed to the common Earth-centred, Earth-fixed (ECEF) coordinate system. The z-axis becomes parallel to the north pole by rotating the local system by an angle 90 - ϕ_geo around the x-axis, with ϕ_geo the geocentric latitude of the camera location. The x-axis aligns with the direction of the prime meridian after rotating the system around the z-axis by an angle λ_geo, where λ_geo the geocentric longitude of the camera location,

    𝐱_geo = 
     [ x_geo; y_geo; z_geo ]
      = R_λ_gdϕ_gd·[ x_hor; y_hor; z_hor ]
    .




 §.§ Trajectory determination





  §.§.§ Meteor geometry



The trajectory of the meteor is determined using the 3D unit vectors of the defined points on the meteor line. These vectors are generated for each camera from the known camera orientation. The vectors point to the meteor trail and are defined in the geocentric coordinate system. Since the meteor line is initially defined in the images using a gnomonic projection, the intersection points of the direction vectors with a unit sphere lie on a great circle. A plane is fitted through all the unit vectors from each station by solving the standard plane equation using least-squares 

    n⃗(x⃗_0-x⃗) = n_xx + n_yy + n_zz + d = 0 
    ,

where x⃗_0 = 0 is the origin of the geocentric coordinate system, x⃗ is the direction vector, ⟨ n_x,n_y,n_z ⟩  are the vector components of the normal vector n⃗,⃗ and d is the distance from the plane to the origin and in this equation is equal to zero. The apparent radiant RA_app, Dec_app of the meteor is calculated as the cross-product of the two normal vectors n⃗_1 × n⃗_2, determined in (<ref>) with the subscripts indicating the two camera stations. The mean altitude of the meteor is computed by intersecting the direction vector of the central point of the meteor from the shuttered meteor station with the plane generated from the direction vectors to the meteor trail from the second station.

To determine the speed and duration of a meteor, each shuttered meteor image is compared with a database of synthetic meteors (see Section <ref>). These meteors have a fixed geometry and orientation with respect to the camera, i.e. the meteor is moving parallel to the x-axis of the camera system and at 100 km above the camera. The projection of the meteor position at time interval t=dt/2 coincides with the principal point of the camera. The database is created by varying two parameters: the speed and duration. The step size of the database is 0.1 km s^-1 for the velocity and 0.02 s for the duration of the meteor. From the known geometric relation between the image and meteor plane, the meteor image is transformed so that the meteor plane becomes parallel to the image plane and the distance between principal point and plane is adjusted to 100 km (Fig. <ref>). This normalised image is then compared with synthetic meteor images in the database accounting for (n_s/0.1)×(n_d/0.02) different combinations for speed and duration, where n_s and n_d are the resolution of our partitioning in speed and duration, respectively. For each combination, the meteor trail is time-shifted by 0.06 s to account for various beginning points. The Pearson correlation coefficient is calculated between a synthetic meteor in the database and the normalised image. For the best match we follow a top-down searching approach: first a coarse search is made and then gradually the step size is decreased around the parameters showing a higher correlation. The speed (V_obs) and duration of the meteor are derived from the synthetic image with the highest correlation.



A meteoroid experiences a deceleration when it reaches the denser layers of the Earth's atmosphere. This effect, so-called atmospheric deceleration, depends on the initial speed of the meteoroid and is more prominent for slower meteoroids. In our studies, we are focussing on the fast-moving Perseid meteoroids and therefore, deceleration is ignored here. The Earth's rotational velocity contributes an extra 0.004 ^∘/s to the calculated right ascension angle of the radiant and is also neglected in this study.

The speed of a meteoroid slightly increases as soon as it experiences the Earth's gravitational attraction, a phenomenon known as zenithal attraction. This effect is computed by performing two integrations following <cit.>: one integration backwards in time including the gravitational effects of the Earth-Moon system until the meteoroid reaches the Earth's sphere of influence and a second integration forwards accounting only for the masses of the Sun and the planets. As input for both integrations the state vector of the meteoroid is used. The new state vector yields the position and velocity of a meteoroid at the time it was recorded in the absence of the Earth-Moon system. The velocity vector now points to the geocentric radiant (RA_geo, Dec_geo).



 §.§ Heliocentric orbit


The orbital path of a meteoroid around the Sun, requiring knowledge of Earth and Sun positions, is computed using standard solar system ephemerides (DE-421). We use the SPICE software library <cit.> to access ephemeris data and retrieve the following geometric transformations.
First the state vectors are transformed from an Earth-centred to a Sun-centred ecliptic coordinate system in J2000, i.e.


    r⃗_eclip = r⃗_geo· R_geo2eclip
    .


The heliocentric position and velocity vector are simply computed as the vector sums

    r⃗_hel = r⃗_met + r⃗_earth
    ,


where r_met and r_earth are the state vectors of the meteoroid and the Earth with respect to the Sun in the heliocentric ecliptic coordinate system. Finally, the osculating elements of the orbit are determined using the heliocentric state vector from SPICE routines.



 §.§ Photometric reduction





  §.§.§ Meteor photometry



Photometric information on meteors is extracted by deconvolving the emitted light of a meteor from the registered signal in equal time intervals <cit.>. We remove the effects of radial distortions in the raw image and resample it using inverse distance weighting interpolation. The displacement values Δx_ij and Δy_ij for each pixel are determined from the geometric camera calibration (Section <ref>). To speed up the interpolation process, we limit the interpolation to a rectangular area around the meteor trail, for which the position is defined (see Fig. <ref>). The change of the angular velocity of a meteor owing to perspective distortion is taken into account by projecting the previously determined 3D meteor path (Section <ref>) to the image. 

The number of time intervals n_t for which the brightness of the meteor is estimated is computed as the ratio of the length of the rectangular area to the spatial sampling resolution defined by the user. A constant spatial sampling size ensures a stable numerical solution for meteors with low angular velocities, but at the same a high-resolution photometric profile for meteor with high angular velocities. From the estimated meteor velocity, the time the meteor needs to cross the rectangular area is calculated following an iterative process. The photometric model can now be applied to the meteor line in the interpolated image.



  §.§.§ Photometric calibration


For photometric calibration we use stars depicted in the image. Their positions in the image (pixel coordinates) are computed using the DAOPHOT routines <cit.> and transformed to the equatorial coordinate system at J2000. The stars are identified by querying the VIZIER database <cit.> and matching them to the brightest stars (m<8) found within a radius of 30 arcminutes (∼4 pixel) from their position. The flux of each star is measured by defining three circular areas around the light source: an inner circular area measuring the light coming from the star and an outer ring determined by two circular areas defining the sky background. We set the star aperture to a radius of 3×FWHM, which encloses nearly 100% of the stellar flux <cit.>. The instrumental magnitude m_inst is then defined as

    m_inst=A-2.5log_10((∑_i=1^nC_i)-nC_sky/t)
    ,

where A is an arbitrary constant, C_i is the DN value in the i^th pixel, C_sky is the mean sky background value, n is the number of pixels in each aperture, and t is the integration time of the frames.

To transform the computed instrumental magnitudes to a standard photometric system, we first convert the Hipparchos H_p magnitudes from the Vizier database into Johnson V magnitudes using the following expression <cit.>:

    V=H_p+a_1(B-V)+a_2(B-V)^2+a_3(B-V)^3+a_4
    ,

where B-V is the colour index of each star from the VIZIER database and α_i the transformation coefficients. The light emitted by each star is partially absorbed by the Earth's atmosphere. Therefore, the amount of the absorption for each star is proportional to the amount of atmosphere the light has to traverse to reach an observer on the Earth's surface. This means that light of stars appearing close to the horizon experiences a greater absorption than stars close to the zenith. The amount of atmosphere, called airmass, is calculated as

    X=sec(z)-0.00186(sec(z)-1) 
       -0.002875(sec(z)-1)^2 - 0.0008083(sec(z)-1)^3 ,

where z is the zenith angle <cit.>. The equation is taking into account the curvature of the Earth. Moreover, the attenuation of light is computed as a function of the wavelength due to Rayleigh scattering and therefore, the amount of attenuation for each star depends on its colour. To account for the colour difference of our star field, we apply a colour correction for each star using the colour indices from the star catalogue. The instrumental magnitudes are corrected for atmospheric effects and converted to absolute magnitudes,

    m_calib=m_inst+T_cCI-Xk-Z_pI
    ,

where X is the airmass, T_c the transformation coefficient, CI the colour index, k is the extinction coefficient given in magnitudes per unit airmass, and Z_p is a scaling factor. We compute three sets of correction parameters for U,V, and I colour corrections using least squares. 

We tested our photometric calibration module with a typical SPOSH image. We detected 393 stars in the image, where the faintest is of +6.3 magnitude. A high correlation between calibrated and standard stellar magnitudes using the V-I colour index was found, matching the spectral response of the system. Stars with an airmass greater than 3 were excluded from the procedure. Figure <ref> shows the calibrated magnitudes of the stars with respect to their catalogue magnitudes from a single frame. The standard deviation between catalogue and measured magnitudes was of 0.22 magnitudes.





§ ERROR PROPAGATION



The errors of the unknown parameters are calculated by applying error propagation. In the sections to follow we refer to these as the propagated uncertainty (or propagated errors) to distinguish  this uncertainty from the statistics of differences between estimated and a priori known parameters (estimated uncertainty) as well as the uncertainty in estimating a common property of the meteors, for example the radiant and speed of a shower, by taking the average over a number of meteors (observed uncertainty). We encounter the estimated uncertainty principally in tests with our synthetic data (Section <ref>). The unknown parameters are the apparent and geocentric radiant positions, observed, geocentric, and heliocentric speed of the meteoroid, and orbital elements of its orbit around the Sun. The observed quantities are the parameters of the meteor line ρ and θ. The general law of error propagation is of the form

    C_yy=∂ y/∂ xC_xx∂ y/∂ x^T
    ,

where C_xx is the stochastic model of the measurements and y are the parameters to be estimated.  The parameter C_yy is the variance-covariance matrix of the unknown parameters. The uncertainties of the direction and location of the meteor line on the image affect the uncertainties of the parameters and need to be carefully estimated. We used our synthetic meteor dataset (see Section <ref>) to estimate the line uncertainties. The distribution of the estimated uncertainties for ρ and θ are well approximated by a normal distribution with a standard deviation of 0.07^∘ for θ and 1.35 pixel for the distance of the projection to the line. These values depend highly on the length of the meteor trail, PSF, and resolution of the CCD. 



§ SOFTWARE VALIDATION





 §.§ Synthetic meteor data



We verified our software modules with the help of synthetic meteor data. A meteor trail is generated by providing a number of parameters which i) define the dynamic and photometric properties of a meteor, ii) define the geometric relation between the observers and the meteor, and iii) projects the luminous path of the meteor to the image plane of the given camera system. Table <ref> summarises the initial conditions used to generate synthetic meteor trails. 

Once the meteor path is generated in space, the corresponding meteor trail is projected on the image plane. The trail is created by convolving a 2D Gaussian curve imitating the motion of a point-like light source on a given camera system. The method is based on the photometric model in <cit.> implemented in reverse. The peak intensity value of each meteor is kept constant while the standard deviation of the Gaussian PSF is set equal to one pixel. The brightness of each synthetic meteor is normally distributed along the meteor trail. The peak brightness also varies between each meteor and resembles different shape curves <cit.>. The position of the peak along the meteor trail in our sample follows a normal distribution with its centre at n_t/2 and a standard deviation of n_t/10, where n_t is  the number of time intervals. The meteor trails in one of the stations were chopped periodically to simulate the effect of the rotating shutter. The starting point of a meteor at time t_0 is placed randomly within a shutter break and ranges between 0 s and 0.06 s. As an example, a meteor with t_0=0 will receive light directly, while a meteor with t_0=0.06 occurs exactly at the time when the shutter is located in front of the lens. For the parameters of the inner orientation of the camera, typical values for the SPOSH camera were used. The pointing of both cameras was chosen to slightly deviate from an optical axis parallel to the zenith. The distance between the two stations was set at 55.6 km. All synthetic meteors in our simulations occur at the same time, i.e. time information is not relevant. The position of the radiant is given in equatorial coordinates system as RA and Dec.



An image with random noise was generated using the noise properties of the SPOSH images that are the same size as the meteor image. Additionally, 20 2D-Gaussian PSF simulating the stellar sources in the image were distributed randomly in the FOV of the camera. The image with the synthetic stars was added to the noise image. The position and brightness levels of the stars were kept fixed for all synthetic images. The light of the synthetic stars in the noise image represent the remaining light due to scintillations, visible in the SPOSH difference image data. Finally, the noise image was added to the meteor image, creating the input for our algorithm. The software uses the images of each synthetic meteor as input and calculates its radiant position, speed, brightness, and heliocentric orbit.

We generated a dataset of synthetic meteor trails considering different geometric configurations. We present results for two types of synthetic data. For the first, we created 208 synthetic meteors with random positions and directions with respect to the location of the cameras, and then used our program to estimated their radiants, speeds, and magnitudes. Forty-seven of the synthetic meteors had convergence angles Q ≤ 10^∘ yielding large errors in the radiant determination. These meteors were therefore excluded from the procedure. For the remaining 161 synthetic meteor pairs with Q > 10^∘ we determined the radiant position (RA and Dec) and the speed (V) and computed their estimated uncertainties as the difference between the a priori value and that calculated from the code. We describe the statistical dispersion of the probability distributions by calculating the median absolute deviation (MAD), which is statistically a more robust measure for asymmetric distributions than the standard deviation. The distributions of the estimated uncertainties for RA and Dec were centred at zero with a median deviation of 0.11^∘ and 0.07^∘, respectively (Fig. <ref>). The statistical properties of the estimated uncertainties are shown graphically in Figure <ref>. The propagated uncertainties had a median value of 0.33^∘ for RA and 0.16^∘ for Dec (Fig. <ref>). The propagated and estimated errors are in good agreement, which implies a realistic stochastic model (Table <ref>). The distribution of the estimated errors for the speed appears to be offset from zero with a median of 0.24 km s^-1 and a median deviation of 0.51 km s^-1. For calculating the meteor magnitudes, we used a subset of the data consisting of 185 un-shuttered meteors of which the individual residual is ≤ 0.5 m. The estimated uncertainties had a median value of +0.01 m and a median deviation of 0.03 m.







A second set of 208 synthetic meteors was then created with the same radiant point for all meteors placed at the local zenith to simulate a meteor shower observed by the two cameras. One hundred seventy of the meteors had a convergence angle Q > 10^∘. Nine radiants with large estimated uncertainties were excluded from the procedure. As for the first set of synthetic data, the estimated errors for RA and Dec also have a zero median but slightly lower median deviation of 0.05^∘ and 0.03^∘ , respectively (Fig. <ref>). These reduce to 0.03^∘ and 0.03^∘ when considering only 126 meteors occurring > 40^∘ from the local horizon. The dispersion of the propagated uncertainties for RA and Dec are similar to those computed for the first synthetic dataset. The median propagated uncertainty was 0.37^∘ for RA and 0.27^∘ for Dec (Fig. <ref>). Figure <ref> shows the relation between the position of a meteor in the image and the estimated errors in radiant. Since the pointing of the camera and the radiant of the simulated shower was set to the local zenith, the angular separation between the radiant and position of a meteor in the image corresponds to the elevation angle of the meteor. The median for the estimated errors in the speed for these 161 meteors was 0.12 km s^-1. The median of the estimated errors for 192 un-shuttered synthetic meteors, for which the residuals in magnitude Δmag were < 0.5m, was +0.01 m with a median deviation of 0.02 m. 







 §.§ Real meteor data



We applied our software to three hours of double-station SPOSH image data acquired during an observing campaign held in Greece (i.e. at ∼37^∘N latitude ) on 12 August 2015 from 23 to 02 UT. The cameras were pointing to the zenith with their x-axis orientated to the north. The baseline between the two sites was 51.5 km. We reduced 177 meteor image pairs and determined their trajectories, velocities and heliocentric orbits. 



We focus on a 20×20^∘ area centred at RA=46^∘, Dec=58^∘, close to the nominal radiant position of the Perseids (Fig. <ref>). To distinguish between Perseid and non-Perseid meteors, we performed a classification based on radiant position and speed as follows: 132 meteors were found to radiate from within this area. We determine the radiant of the Perseid shower as the median value of these radiants: RA=45.96^∘ and Dec=57.77^∘. We assume that most of the meteors are Perseids and we find the 1σ uncertainty in RA and Dec to be 3.29^∘ and 2.27^∘, respectively. The geocentric speeds V_G have a median value at 58.97 km s^-1 and a median propagated error of 1.21 km s^-1 (Fig. <ref>). We classify all meteors with speeds closer to this median speed than four times the median propagated error, as Perseids. In this way, we identified 71 meteors belonging to the Perseids meteor shower (Fig. <ref>). Their median speed V_G was found to be 59.58 km s^-1 with a median deviation of the observed uncertainties of 0.48 km s^-1. Statistical properties are given in Table <ref>. As the aim of our example is to demonstrate the successful usage of our software using real data, we neglected the effect of radiant drift.







We calculated the magnitudes of the 71 meteors identified as Perseids from the un-shuttered images. We defined this magnitude to be the brightest value obtained for the light curve of each meteor. The magnitude distribution index r for the shower was found to be 2.10 ±0.10 (Fig. <ref>). The mean value is slightly above the upper limit of the range 1.87<r<2.01 given by <cit.> during the years 1988-1994 but within the range 1.86<r<2.12 found by <cit.> for the years 1991-1994. These observations cover both the pre- and post activity of the meteor shower that coincided with the perihelion passage of the parent comet in 1992. The index was computed from 61 meteors brighter than +0 m. The light curve of a bright, double-flaring Perseid was computed following the method described in Section <ref> (Fig. <ref>). The time step size dt was 0.006 and 0.01 seconds for the un-shuttered and shuttered camera respectively. The gaps between the points represent the time intervals with no information owing to the rotating shutter.





Heliocentric orbits were computed for all 177 meteors in our sample. The median orbital elements of the 71 Perseid meteors were compared with the orbits found in five studies <cit.> and the orbit of comet 109P/Swift-Tutle, parent comet of the Perseids <cit.> (Table <ref>). In general, we find a good agreement in the radiant, speed, and orbital elements. One exception is in the semi-major axis, which is known to be very sensitive to variations in the velocity of a meteoroid <cit.>.



§ CONCLUSIONS


We have presented SPOSH-Red, a software package for the data reduction of double-station meteor image data acquired by the SPOSH camera. The software extracts information about trajectories, heliocentric orbits and brightness levels of recorded meteors. 

The software was tested for simulated and real meteor data. We simulated different geometric configurations between a meteor shower and two observing sites. We suggest that such simulations can be used to assess the quality of the derived meteor trajectories for different camera network configurations and predicted meteor shower or outbursts. We expect that the results will greatly contribute to the planning of observing campaigns by finding the best location and orientation between the camera stations and predicted radiant position.

The software presented in this paper was developed to reduce data acquired by the SPOSH camera system. In the future we plan to provide a more generic version of the software package that can handle image datasets recorded by different camera systems.

The real meteor data used in this work is part of a large dataset comprising eight years of Perseids observations using SPOSH. The accumulated observing period spans more than 20 days within the activity period of the Perseids. During this period > 15,000 single meteors have been recorded from both stations. The reduction software opens the opportunity to analyse the available unique SPOSH meteor data. A full analysis of the data focussing on the Perseid meteor shower will be presented in a forthcoming paper.



§ ACKNOWLEDGMENT

The authors would like to thank the reviewer and editor for providing us with useful comments and detailed suggestions for the manuscript.

aa 


