© S. Jeanne et al. 2019

1. Introduction

Most of the fireball networks, such as the European Fireball Network (Oberst et al. 1998) and the Desert Fireball Network (Bland et al. 2012) use fish-eye lenses to survey the entire sky. However, the small focal length and the heavy distortion of these lenses induce reduction problems. In addition, external errors are difficult to determine, so that errors on orbits and strewn fields are poorly constrained. This problem was first addressed by Ceplecha (1987), taking advantage of the fact that trajectories are almost linear, and that they can be determined based on observations from two different stations. Borovička et al. (1995) subsequently improved on this methodology by using a least squares method to fit all the observations at the same time. This procedure is at the heart of all the methods used thereafter by the different fireball networks (Bland et al. 2012; Weryk et al. 2008; Trigo-Rodriguez et al. 2007), but the geometrical distortion and the error estimations differ among the projects. All networks have to deal with low elevation errors (h <  15°), which constitute a serious problem as most of the final trajectories are seen low on the horizon. We developed a new method, inspired by Hivon et al. (2017) and based on a nine degree odd polynomial distortion function. This method allows us to perform measurements at a 5° elevation with an accuracy of 1 arcmin or 0.1 pixel. As orbit accuracy strongly depends on measurement error, a reliable error estimate is needed. While the internal error is more easily determined than the external error, the latter is critical to the final result. Borovička et al. (1995) elected to estimate it based on the quality of the trajectory fit. In this paper, we propose a new external error estimation method based on observation errors between stations being independent from one another. This makes it possible to use several stations to perform a Gaussian analysis and therefore a global analysis of external errors.

The internal and external errors are then propagated on trajectories and velocities. The latter task is rendered difficult as observations are never reproducible (Egal 2017). However, such estimates are needed to compute the orbital parameters of meteoroids. In particular the measurement of the velocities is critical: a small error in velocity leads to significant errors in semi-major axis and eccentricity. In consequence, poorly constrained uncertainties in the trajectories and velocities of fireballs can lead to major errors in orbital element estimates. In the present paper, we detail the analytical method we developed to determine trajectories and their confidence intervals for fireballs seen by cameras that show all types of optical aberrations. Such aberrations include asymmetrical defects, for example staggered or slanted lenses, and cameras installed on structures that vibrate with the wind, or dilate with temperature.

FRIPON network data are used to the test our method. This fireball network was funded by Agence Nationale pour la Recherche (ANR) in 2014 (Colas et al. 2012). It now comprises 95 cameras (Colas et al., in prep.), and has started expanding to Europe (Gardiol et al. 2016; Koschny et al. 2017; Nedelcu et al. 2018). All the optical observations are made with the same hardware configuration (Colas et al. 2015). The heart of the network is a digital camera using a fish-eye lens. As the mesh of the network is 80 km, an event can be detected by many cameras, which can be used to estimate external errors on the measurements. The method described in the present paper is the cornerstone of the FRIPON reduction pipeline.

2. Astrometrical analysis

2.1. FRIPON network and data

All the observations of the FRIPON fireball network are made with the same hardware configuration (Colas et al. 2015): both for the fish-eye lens (f = 1.25 mm) and for the detector (Sony ICX445). Fish-eye lenses are known to produce heavy distortion, and we were unable to obtain the optical design from the lens manufacturer. Furthermore we found a strong heterogeneity among the distortion figures in the same set of lenses, and among the positions of the chip inside the camera. Hence, we developed our own distortion model to be fitted on each camera, as we needed to compute calibrations with background stars (Ceplecha 1987; Bland et al. 2012). Our objective is to get a good calibration for elevations superior to 10°. Indeed, with a characteristic distance between two cameras of about 80 km, a meteor at 20 km altitude will be seen with an elevation down to 10° by six or seven cameras; but on average only one or two cameras will see it at an elevation superior to 20°. So a calibration from zenith to 10° allows us to observe the whole trajectory of a large majority of bolides. Other networks, such as the Desert Fireball Network (Devillepoix et al. 2019) use better resolution cameras with larger spacing, but the need for low altitude astrometry is the same. As we used a video acquisition frame rate of 30 fps, the limiting magnitude is only zero, so there are almost no stars visible on a single frame. In order to detect enough calibration stars, longer exposure frames are needed. A pause of even five seconds is not sufficient to efficiently sample the field of view, as only a few dozen stars are measured. Hence, the measurement needs to be repeated and this is done every ten minutes so that the acquisition is disturbed as little as possible. With such a strategy, for a clear night, a few thousand stars are measured per camera and they fill the whole field of view as shown below in Sect. 5.2 (Fig. 4).

As the camera mount is stable, this calibration can be used over several nights in the case of a partly cloudy sky. Photometric corrections (flatfield, bias, and dark) are not used as such calibration images are not obtainable with remote fish-eye cameras. Instead, we computed the median of all the calibration images of the night, which automatically removes the stars as they move with the sky. This produces a map of the hot pixels and other static features related to the camera itself, which can in turn be subtracted from each frame. All the frames (calibration and detection) are then processed with SExtractor (Bertin & Arnouts 1996) to produce pixel positions (x, y). These calibrations are processed with Scamp (Bertin et al. 2006) to get astrometrical and photometrical plate solutions. As Scamp was not programmed for fish-eye observations, it only works at elevations above 30°, where the distortion is lower. This solution serves as the first order for our global model. The output of Scamp is a file associating objects found on the CCD (x, y) with topocentric positions (Alt, Az) and celestial positions for the observation epoch (Ra, Dec).

2.2. Reference frames and geometrical parameters

The camera reference frame is a Cartesian coordinate system: O_z is the camera’s optical axis and O_x and O_y are defined by the orientation of the CCD chip. The goal of the calibration is to connect this reference frame with the Earth-fixed one International Terrestrial Reference Frame (ITRF), common to all cameras and necessary to compute the trajectory. Stars are used to fit a model to link these two reference frames as star positions are fixed into the celestial reference frame ICRF2 (Gaia Collaboration 2018), making it easy to convert their positions into the Earth-fixed reference frame ITRF. The geometrical parameters are the three angles describing the orientation of the camera (θ₁, θ₂, θ₃) in the Earth-fixed reference frame and the optical centre position (x₀, y₀) on the camera CCD.

At the beginning of the process, the geographical position of the camera is not used and no hypothesis is made on its orientation. The latter has to be determined to calibrate each new camera by exploring all possible directions and optical centres. The determination of the distortion function as described below allows us to convert a position in the image (x, y) into coordinates (Alt, Az) in the camera reference frame.

2.3. Distortion function

In addition to determining the geometrical parameters characteristic of each camera, its distortion function has to be modelled. In the CCD frame, the position of a star is described by the distance to the optical centre R and the argument A (Eq. (1)) where (x₀, y₀) is the position of the optical centre. In the first instance, we assumed that fish-eye lenses are perfectly symmetrical so that a radial distortion model is sufficient. In this case, the azimuth of a star (in the camera’s frame) is directly measured by the angle A.

$$\begin{array}{c}\hfill R=\sqrt{{(x-x_0)}^2+{(y-y_0)}^2}tan(A)=\frac{y-y_O}{x-x_0}·\end{array}$$(1)

The zenith distance z is given by a distortion function f(R) shown in Fig. 1. Kwon et al. (2014) show that it is possible to reach an accuracy of 0.4 pixel at 10° elevation with a standard equisolid fish-eye distortion model described by Eq. (2). However, this was not usable for lower elevation.

$$\begin{array}{c}\hfill R=k_{1}sin(k_2·z).\end{array}$$(2)

Fig. 1. Zenith distance of stars as a function of the distance R to the optical centre. Comparison of recorded stars (in black) with the equidistant model (in red).

Most lens manufacturers attempt to design lenses in which the distortion more or less follows the equisolid model (Bettonvil 2005). However, this is not always the case and, due to manufacturing processes, fish-eye lenses often deviate from this canonical projection. For a better accuracy Hivon et al. (2017) use an odd power law of sin(k₂ ⋅ z), the main term of the equisolid distortion model. With a very different setup (camera and lens), they reach an accuracy of 1 arcmin with a smaller pixel size of 3 arcmin, even at low elevation. We found this method unnecessarily complex, and followed Hughes et al. (2010) by using a standard odd polynomial function that can model any distortion feature (Eq. (3)). This model can also be understood as the modified development of a sin function. Our tests show that a nine-degree odd polynomial function is sufficient to model the distortion of FRIPON fish-eye camera as shown in Fig. 2.

$$\begin{array}{c}\hfill z=a_1·R+a_3·R^3+a_5·R^5+a_7·R^7+a_9·R^9.\end{array}$$(3)

Fig. 2. Residuals in elevation, using a nine degree odd polynomial function Eq. (3) to model the distortion. The standard deviation is 2.12 arcmin. Toulouse fish-eye camera.

2.4. Asymmetrical correction

Anisotropic systematic effects are present in some cameras. Such an effect is particularly visible for the Brussels FRIPON camera, see Fig. 3. In order to correct for it, we used an empirical correction formula according to Borovička et al. (1995):

$$\begin{array}{c}\hfill R=R_{\mathrm{old}}(1+K_{1}sin(A+\varphi )).\end{array}$$(4)

Fig. 3. Residuals in azimuth for the Brussels FRIPON camera fitted with a symmetrical distortion function (top) and with an asymmetrical correction (bottom). The standard deviation is 4.76′ without the asymmetrical correction. It goes down to 3.93′ with the asymmetrical correction.

This formula adds two new parameters to our model: an asymmetrical parameter K₁ and a phase ϕ. R_old is given by Eq. (1). The improvement provided by this correction is shown in Fig. 3. In this case, the asymmetrical factor K₁ is about 2.5 × 10⁻³.

2.5. Summary

Our astrometrical model uses 12 parameters that can be separated into four categories: the three angles describing the orientation of the camera (θ₁, θ₂, θ₃), the optical centre position (x₀,y₀), five parameters for the radial distortion function (a₁, a₃, a₅, a₇, a₉) and two parameters for the asymmetrical correction (K₁, ϕ) to account for chip inclination or lens misalignment. The median standard deviation of the positions of the stars calculated with this model is typically 2 arcmin or 0.2 pixel, as shown in Fig. 2, up to 4 arcmin in the worst cases. In the example described in Sect. 5, we show that, for a bright meteor, the internal error is of 0.7 arcmin. However, as the point spread function (PSF) with a FWHM of 1.5 pixel is under-sampled, it seems difficult to achieve an internal error better than 1/20 pixel, corresponding to 15 m at a 100 km distance, quite good enough for meteorite recovery. Another potential source of error is the motion of the fireball in a single frame. With an average velocity of 20 km s⁻¹ at a 100 km altitude, a meteor produces a 20 arcmin trail on a FRIPON CCD. With a 15 arcmin PSF, this produces only a slightly elongated feature, and induces no significant discrepancy in accuracy along and across the trajectory.

3. Data reduction

3.1. Trajectories

FRIPON cameras record 30 images per second. As a meteor has a velocity between 11 and 72 km s⁻¹ when it enters into the atmosphere, the distance travelled within one frame is between 0.4 and 2.4 km, equivalent to 1–6 pixel at a 100 km distance. The acquisition computer timing is obtained from an Network Time Protocol (NTP) internet server. Typical median timing accuracy is about 20 ms (Barry et al. 2015), but it can be worse if the internet connection is off, and the computer therefore depends on its internal clock for timing. As a result, the camera acquisitions are not synchronised as those of Desert Fireball Network are (Sansom et al. 2019). This lack of temporal accuracy combined with the high speed of the meteors prevents us from computing their positions by a direct triangulation as in Sansom et al. (2019), hence we have not directly computed the position of the meteor as a function of time as Gural (2012) does. Instead, we computed its trajectory assuming a straight-line after Ceplecha (1987). The relevance of this assumption is discussed in Sect. 4.1. This assumption allows us to separate the space and time components of our measurements and to overcome the problem of temporal precision. But it is possible to measure an internal timing accuracy referenced to an average time as displayed in Table 1. In this example, the timing accuracy is of the order of 10 ms, which is better than expected, perhaps due to the fact that every computer uses the same NTP address.

A straight-line in a three dimensional Euclidean space is defined by four parameters. It is necessary to first fit the best straight-line from all camera observations. Ceplecha (1987) computed a straight line per couple of cameras and, finally, a weighted average line from all these lines. For each all-sky camera, a plane through the camera is defined as the plane that minimises the square sum of scalar products between a vector normal to this plane and the measured vectors pointing from the camera towards the fireball (Ceplecha 1987). As the camera positions are known, we measured the vectors (x_j, y_j, z_j) normal to these planes:

$$\begin{array}{cc}\hfill {\mathit{x}}_j& =\left({\displaystyle \sum _i^{n_j}{\mathit{a}}_{\mathit{ji}}{\mathit{b}}_{\mathit{ji}}}\right)\left({\displaystyle \sum _i^{n_j}{\mathit{b}}_{\mathit{ji}}{\mathit{c}}_{\mathit{ji}}}\right)-\left({\displaystyle \sum _i^{n_j}{\mathit{b}}_{\mathit{ji}}^2}\right)\left({\displaystyle \sum _i^{n_j}{\mathit{a}}_{\mathit{ji}}{\mathit{c}}_{\mathit{ji}}}\right)\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill {\mathit{y}}_j& =\left({\displaystyle \sum _i^{n_j}{\mathit{a}}_{\mathit{ji}}{\mathit{b}}_{\mathit{ji}}}\right)\left({\displaystyle \sum _i^{n_j}{\mathit{a}}_{\mathit{ji}}{\mathit{c}}_{\mathit{ji}}}\right)-\left({\displaystyle \sum _i^{n_j}{\mathit{a}}_{\mathit{ji}}^2}\right)\left({\displaystyle \sum _i^{n_j}{\mathit{b}}_{\mathit{ji}}{\mathit{c}}_{\mathit{ji}}}\right)\hfill \end{array}$$(6)

$$\begin{array}{cc}\hfill {\mathit{z}}_j& =\left({\displaystyle \sum _i^{n_j}{\mathit{a}}_{\mathit{ji}}^2}\right)\left({\displaystyle \sum _i^{n_j}{\mathit{b}}_{\mathit{ji}}^2}\right)-{\left({\displaystyle \sum _i^{n_j}{\mathit{a}}_{\mathit{ji}}{\mathit{b}}_{\mathit{ji}}}\right)}^2\hfill \end{array}$$(7)

where the n vectors (a_j, b_j, c_j) are unit vectors pointing from the cameras towards the meteor and n_j is the number of images taken by the jth camera. When a meteor is observed by n_c cameras, we have n_c planes defining $\left(\genfrac{}{}{0pt}{}{2}{n_{\mathrm{c}}}\right)$ intersections. Of course, these $\left(\genfrac{}{}{0pt}{}{2}{n_{\mathrm{c}}}\right)$ trajectories are different, because of measurement errors and the assumption of the straight-line model.

Another method, proposed by Borovička (1990), is a global least squares regression minimising the sum of squares of angular distances between the lines of vision and the trajectory. We chose the latter method, to which we add a specific weight attached to each camera (see Sects. 3.2 and 3.3), in order to increase the robustness of the results. Meteor positions over time were then calculated by projecting the lines of vision of each image onto this line. The modified sum of residuals is presented in the section below.

3.2. Modified sum of residuals

We use a least squares regression to fit the four parameters of the straight-line trajectory. If the measurements are not biased, it is possible to get the internal uncertainties. By making this assumption, the sum of squares of residuals is:

$$\begin{array}{c}\hfill S_1(\mathit{\theta })={\displaystyle \sum _j^{n_{\mathrm{c}}}\sum _i^{n_j}\frac{{ϵ}_{\mathit{ji}}{(\mathit{\theta })}^2}{{\sigma }_j^2}}\end{array}$$(8)

where n_c is the number of cameras and n_j the number of images taken by camera j. ϵ_ji(θ) is the difference between the measurement and the model, assuming a vector of parameters θ. σ_j is the standard deviation expected for camera j. In this case, the uncertainty on the vector of parameters θ decreases with the number of recorded images.

However, as mentioned above, the images from any given station may have systematic errors related to the station itself, for example in the case of structures that evolve with time and temperature or located at the top of masts that vibrate with the wind. In addition, some cameras may have an inclined or misaligned lens. Lastly, the shape and thickness of the spherical porthole may vary, especially low on the horizon. The cameras could not be tested individually in the laboratory, and only the focus was adjusted. For these reasons, a few cameras may have flaws, which result in a permanent bias, and as a result the images from the individual stations are not independent from one another.

Taking this into account, if we assume that the astrometrical measurement noise is negligible compared to the systematic errors, the least squares method can still be used, except that our measurements have to be redefined. It can be assumed that each station is independent of the others. This assumption is more acceptable than the images of a given camera being independent from one another. Hence, we considered each set of images from one camera as one unique measurement, instead of each individual image being treated as a measurement. In the following text we use the term “video” for the whole set of calibration images of any given station. We obtain this new sum of residuals S:

$$\begin{array}{c}\hfill S_2(\mathit{\theta })={\displaystyle \sum _j^{n_{\mathrm{c}}}\frac{E_j{(\mathit{\theta })}^2}{{\mathrm{\Sigma }}_j^2}}\end{array}$$(9)

where E_j(θ) is the o − c (observed minus calculated, equal to the difference between the video j measurement and the model), assuming ${\mathrm{\Sigma }}_j^2$ is the standard deviation expected for the video measurement j of parameters θ:

$$\begin{array}{cc}& E_j{(\mathit{\theta })}^2={\displaystyle \sum _i^{n_j}{ϵ}_{\mathit{ji}}{(\mathit{\theta })}^2}\hfill \end{array}$$(10)

$$\begin{array}{cc}& {\mathrm{\Sigma }}_j^2=n_j{s}_j^2\hfill \end{array}$$(11)

where s_j is the systematic error expected for camera j. In this case, the uncertainty on the vector of parameters θ decreases only with the number of stations, and not with the number of recorded images. This results from our assumption that the error is dominated by the systematic errors for each station, making the images from a given station dependent on one another, so that the addition of new images does not contribute new information.

Hence, there are two ways to compute a sum of squares of residuals. The first method applies in the case where our distortion model is perfect for each station and errors only come from normal noise, each image being considered as one independent measurement. The sum of squares of residuals based on such an assumption is described by Eq. (8). In the second case, in which images from one video are dependent on one another and the video becomes one unique independent measurement, their square residuals are summed as shown in Eq. (9). A more general case can be constructed by combining these two cases in a third sum of residuals S:

$$\begin{array}{c}\hfill S_3(\mathit{\theta })={\displaystyle \sum _j^{n_{\mathrm{c}}}\sum _i^{n_j}\frac{{ϵ}_{\mathit{ji}}{(\mathit{\theta })}^2}{{\sigma }_j^2+n_j{s}_j^2}·}\end{array}$$(12)

This new modified sum of squares of residuals allows us to combine the measurements from different cameras with different characteristics. Moreover, it allows us to highlight the transition between errors dominated by the internal noise and errors dominated by the bias of each the camera.

3.3. Estimation of σ_j and s_j

In order to compute our modified sum of residuals (Eq. (12)), we need to characterise the standard deviation of measurements σ_j and the confidence parameters s_j for each camera. In other words, we need to separate the astrometrical noise from the systematic errors.

For σ_j, we use the standard deviation of the measurements of the jth camera referenced to the plane defined by the best fit of these measurements. This plane is characterised by the position of camera j and its normal vector (x_j, y_j, z_j) described in Eqs. (5)–(7). As this plane is only computed with the measurements of camera j, it will hide the systematic errors related to that station, so that the standard deviation will be close to the real astrometrical noise for camera j.

For s_j, we needed to characterise systematic errors in the measurements, which depend on the flaws and imperfections of the station itself. However, they also depend of the path of the meteor inside the images. Indeed, a difference between reality and our distortion model will not lead to a global degradation of the measurements, but to localised systematic errors. For example, a water droplet on the spherical porthole will induce a systematic error if the light crosses this droplet. So the quality of measurements on the path of the meteor needs to be characterised. For each image i of the meteor, we computed the average $\overline{{\mathit{ϵ}}_{\star ji}}$ of errors on stars in an area around the meteor. This area extends both in space (a circle around the meteor image) and in time (clearly, stars measured only at the same time as the meteor cannot be used). We chose to use stars seen during the same night as the fireball within a radius of 100 pixels around the fireball. The choice of such an area is discussed in Sect. 4.2. The standard deviation to zero of these error averages is taken as s_j:

$$\begin{array}{c}\hfill s_j^2=\frac{1}{2n_j}{\displaystyle \sum _i^{n_j}{\overline{{\mathit{ϵ}}_{\star ji}}}^2.}\end{array}$$(13)

The dividing factor of two is here to convert an uncertainty in a two-dimensional space (the celestial sphere) to an uncertainty in a one-dimensional space: the angular distance at the meridian, defined by the projection of the straight-line trajectory on the celestial sphere.

4. Discussion

4.1. Straight-line trajectories

Our straight-line model may be objected to as being unrealistic. Indeed a fireball is under the influence of the Earth’s gravity, therefore its trajectory must be curved. This effect is all the more important given that the object is slow and that its trajectory is at a low angle to the ground. The minimum radius of curvature of an object undergoing constant acceleration g is ${\rho }_{\min }=\frac{V^2}{g}$. For a slow fireball with a velocity about 11 km s⁻¹, ρ_min = 12 300 km, or 1.93 Earth radius. Thus, in the worst case of a slow meteor, with a length of ∼100 km and a trajectory parallel to the ground, the maximum difference between a straight-line and the actual curved trajectory is ∼50 m. This is comparable to the limit of the precision that can be reached on trajectories computed by FRIPON as discussed in Sect. 2.5. So, in a large majority of cases, we would not be able to detect the meteor curvature. This is confirmed by Borovička (1990). After having analysed five long fireballs, these authors conclude that it is not possible to detect the meteor curvature from all-sky images. There are some exceptions, such as the two cases described by Sansom et al. (2019). If such bolides are observed by the FRIPON network, they will need to be reduced by a particular method, and not by the regular automatic method described here.

4.2. Estimation of s_j

Equation Eq. (13), giving s_j, depends on the $\overline{{\mathit{ϵ}}_{\star ji}}$ averages. Average $\overline{{\mathit{ϵ}}_{\star ji}}$ is the average of errors on stars in an area around the meteor position in the ith image of the jth camera. The choice of this area is arbitrary, so the estimator of s_j depends on the area size we consider to compute the $\overline{{\mathit{ϵ}}_{\star ji}}$ averages. A big area hides short-frequency variations, but a small area introduces uncertainty because of the smaller number of stars. That uncertainty is inversely proportional to the square of the number of stars in the chosen area. A defect that is localised in space (e.g. a default related to the fabrication of the spherical porthole or the presence of a water droplet) has to be large enough to contain a large number of stars. If not, we are unable to detect it. For all these reasons, we can only estimate long-time and large-space variations. In the present model, we use stars from the night of the fall within a 100 pixel radius around the position of the meteor. Variations under these values will not be considered.

Moreover, if the systematic error for a camera is close to zero, the expected value of $\overline{{\mathit{ϵ}}_{\star ji}}$ will be $\frac{{\sigma }_{\star }}{\sqrt{n_{\star }}}$, where σ_⋆ is the expected standard deviation on the measurements of stars and n_⋆ is the number of stars used to compute $\overline{{\mathit{ϵ}}_{\star ji}}$. This value is a limit, under which we cannot estimate s_j. In other words, we cannot determine if a camera is better than this confidence limit s_limit. In the case for which a camera is actually better than its confidence limit s_limit, uncertainties on the meteor trajectory are overestimated.

These problems might be reduced by increasing the density of sources used for camera calibration. In order to increase the number of stars recorded, we might increase the number of long exposure images per night (at present, the long exposure images are taken every ten minutes to avoid long acquisition disruptions), or we might increase the exposure time of these images (but a five-second exposure is already at the limit of sky saturation in some locations, as many camera are located in the centre of cities). In some cases (e.g. partly cloudy nights), the data for several nights are used, based on the hypothesis that cameras do not move.

4.3. Advantage and disadvantage of the method

The modified least squares method, described in Sect. 3.2 is not perfect, and some objections may be raised. Indeed, this method assumes that each camera is independent of the others. Moreover, it also assumes that the station characteristics, expressed by the term s_j, follow a normal distribution. The first assumption (independence of stations) is solid: each station is several dozens of kilometers away from the others. But the assumption of a normal distribution is much less so. Each camera has its own particularities, its own installation, with its own sensitivity to temperature, wind, etc. So, for a meteor seen by a particular set of only a few cameras, it is almost impossible to describe the distribution law of the errors precisely.

Moreover, the computation of a confidence interval depends on parameters s_j and σ_j for each camera. These parameters are not straightforward, especially s_j, and errors can be made, cf. Sect. 4.2. For all these reasons, there will always remain some level of arbitrary choice, and our confidence intervals have to be viewed with critical eyes. Hence, our uncertainties should be taken as orders of magnitude, not as exact values.

Nevertheless, this method allows us to compute meteor trajectories from a heterogeneous group of cameras and to get an idea of the uncertainties on these trajectories. These uncertainties can be propagated in the calculation of velocity, orbit, strewn field, etc. This method is useful for a large meteor survey system such as FRIPON or in cases where videos from motley sources have to be used (survey cameras, amateur videos, etc.). Moreover, this method allows us to compute trajectories and confidence intervals consistent with one another, hence to compute statistics on the populations of meteoroids that cross the path of the Earth.

5. Application to a fireball

5.1. Characterisation of cameras

As an example of results obtained with the method presented above, we have chosen a fireball event which occurred between Italy, Switzerland, and France at 00:06 UT on August 4, 2017. This fireball was recorded by three of the FRIPON and PRISMA (Gardiol et al. 2016) cameras: St Lupicin, Sutrieu, and Lignan. The stars used to calibrate each camera are displayed in Fig. 4. This figure shows that all cameras are not equivalent. Residuals between computed star positions and our predictions are larger for the Lignan camera than for the others. Moreover, for this station, residuals do not seem independent of their location in the image. This shows that our distortion model, presented in Sect. 2, is not able to fully describe the Lignan camera. This is our reason for choosing this particular event: it allows us to demonstrate how our method behaves in the face of measurements from motley stations.

Fig. 4. Superposition of all calibration stars observed during one night used to compute distortion parameters for the St Lupicin, Sutrieu, and Lignan cameras (from top to bottom). X and Y are the positions of the stars in the images (pixel). Errors between the distortion model and measurements are represented by a colour scale (arcmin). The positions of the meteor of 2017-08-04 UT 00:06 are shown in dark blue.

For these three cameras, we computed estimates of confidence parameters s_j and standard deviations σ_j, with the methods described in Sect. 3.3. These estimates are displayed in Table 2.

It is interesting to compare the astrometrical noise term ${\sigma }_j^2$ and the systematic error term $n_j{s}_j^2$. Data from the Lignan camera are clearly in the case $n_j{s}_j^2\gg {\sigma }_j^2$. All the images are strongly dependent on one another. In this case the full video should be considered as a unique measurement (cf. Eq. (9)).

For the St Lupicin and Sutrieu cameras, $n_j{s}_j^2$ is only about ten times larger than ${\sigma }_j^2$. But, for these two stations, s_j is near to the confidence limit s_limit discussed in Sect. 4.2. There are not enough stars in our calibration sample to estimate the divergence between the distortion model and reality under this limit. This means that these two stations probably have a smaller systematic error than is indicated by their confidence parameters s_j. But we do not have enough stars in our calibration sample to prove this. Conversely, for the Lignan station, s_j is clearly higher than its limit. This describes the divergence between our distortion model and that camera.

5.2. Computing fireball trajectory

We used the least squares method, with the modified sum of residuals described by Eq. (12) with the parameters shown in Table 2. The calculated straight-line trajectory starts at 46.01898° N, 6.77494° E, and at an altitude of 72.96 km. Its ends at 46.28917° N, 7.04590° E, and at an altitude of 43.49 km. The residuals of this trajectory are shown in Fig. 5. The standard deviations of these residuals are 0.71, 0.51 and 1.78 arcmin for, respectively, St Lupicin, Sutrieu, and Lignan. These values are consistent with those displayed in Table 2. The residuals from the Lignan camera clearly do not follow a Gaussian law, but they are compatible with the systematic error s_j expected for this camera. The confidence intervals are shown in Fig. 6. At the top of the trajectory, the 2σ confidence interval is an ellipse with a 220 m semi-major axis and a 50 m semi-minor axis. For the 1σ confidence interval, these values fall to 100 m and 30 m.

Fig. 5. Residuals of the best trajectory in arcmin (top) and metres (bottom) from cameras of Lignan (red), St Lupicin (green), and Sutrieu (blue). Meteor of 2017-08-04 UT 00:06.

Fig. 6. Top: lines of vision from cameras to the meteor of 2017-08-04 UT 00:06, projected in a plane normal to the best trajectory found. This trajectory crosses the plane at coordinates (0, 0). Bottom: zoom on the centre showing random trajectories drawn within 1σ confidence interval (blue) and in 2σ confidence interval (green).

This example was chosen to demonstrate our statistical method. We emphasise that if the lines of view for two stations are close, as for St Lupicin and Sutrieu in the present example, the results are not significantly affected. This example also demonstrates the effects of systematic biases as in the Lignan case. Last, it indicates that an accuracy as high as ∼20 m might be reached in the case of bright events seen by many cameras.

5.3. Photometry and velocity

Based on the astrometrical data reduction presented in Sect. 3, a detected calibration star of the calibration frames (X, Y, F) can be connected with a star from a catalogue (Ra, Dec, Mag). In order to convert the arbitrary fluxes computed by SExtractor (Bertin & Arnouts 1996) with the magnitudes of a photometric catalogue, we used the Hipparcos catalogue (Bessell 2000), as its filter passband is similar to the CCD one. As these calibration frame are made with a 5 s exposure time, we have to apply a correction factor of 1/150 for the flux of the meteor on detection images made with a 30 frame per second acquisition rate. Finally, to compare the light curves originating from different cameras, we used the absolute meteor photometric system, which corresponds to the magnitude of the same meteor as seen from 100 km without atmospheric absorption. The atmospheric absorption is modelled with the simplified law (Eq. (14)):

$$\begin{array}{c}\hfill F=F_z-\kappa ·m\end{array}$$(14)

where F is the flux from a source, F_z is the flux from this source at the zenith, m is the air mass ($m=\frac{1}{cos(z)}$), and κ is the extinction coefficient.

To measure the extinction coefficient k, it is necessary to compute for each star, the flux of a magnitude zero (F₀) star located at the same elevation, hence with the same absorption (Eq. (15)), and then to compute the linear regression of F₀ measured as a function of the air mass m (Eq. (16)),

$$\begin{array}{cc}\hfill F_0& =F·{10}^{\frac{\mathrm{Mag}}{2.5}}\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill F_0& =F_{0z}-\kappa ·m\hfill \end{array}$$(16)

where F_0z is the flux from a zero magnitude star at the zenith. Finally, the magnitude of a fireball in the Hipparcos photometric system is determined as:

$$\begin{array}{c}\hfill {\mathrm{Mag}}_{\mathrm{fireball}}=2.5·{log}_{10}(\frac{F_{0z}-\kappa ·m}{F_{\mathrm{fireball}}}·\frac{t_{\mathrm{s}}}{t_{\mathrm{l}}})\end{array}$$(17)

where F_fireball is the flux measured from the fireball, t_s is the exposure time of short exposure images, and t_l is the exposure time of long exposure images. The absolute magnitude of a fireball is defined as its magnitude at a distance of 100 km (Eq. (18)).

$$\begin{array}{c}\hfill A{\mathrm{Mag}}_{\mathrm{fireball}}={\mathrm{Mag}}_{\mathrm{fireball}}-5·{log}_{10}\left(\frac{d}{100\mathrm{km}}\right)·\end{array}$$(18)

We also determine the position of the meteor as a function of time by projecting the lines of vision on its previously computed trajectory. Altitude, velocity, and magnitude are shown in Figs. 7–9. The fireball magnitude brutally increases at a 45 km altitude, shortly before it disappears. Moreover, the velocity at the last visible point is 10 km s⁻¹. According to Brown et al. (2013), these values indicate that the body was totally disintegrated. This result will be confirmed below in Sect. 5.4 and Table 3. Meteoroids that survive their journey through the atmosphere, can be seen at lower altitudes (under 30 km). In addition, their magnitude does not increase suddenly, but gradually with their deceleration (Campbell-Brown et al. 2013).

Fig. 7. Absolute magnitude as function of altitude for the fireball of 2017-08-04 UT 00:06 seen by the Lignan (red), St Lupicin (green), and Sutrieu (blue) cameras. The absolute magnitude of a fireball is its magnitude at 100 km.

Fig. 8. Altitude of the fireball of 2017-08-04 UT 00:06 as a function of time, seen by the Lignan (red), St Lupicin (green), and Sutrieu (blue) cameras. Time starts at 2017-08-04 UT 00:06:18.159.

Fig. 9. Velocity as function of altitude for the fireball of 2017-08-04 UT 00:06 seen by the Lignan (red), St Lupicin (green), and Sutrieu (blue) cameras. Dots are point to point measurements. Lines are linear approximation on groups of five points.

5.4. Deceleration model

We estimated (18.59 ± 0.17) km s⁻¹ as the velocity at the top of the visible trajectory by using a linear approximation on the first five points of each video. In order to estimate the velocity at the top of the atmosphere, we used a dynamic model from Bronshten (1983), Eqs. (19) and (20), which describes the deceleration and the ablation of a meteoroid in the atmosphere. This model takes four parameters: the pre-entry velocity V_e (our unknown variable), a deceleration parameter A, an ablation parameter B, and a shape change coefficient μ.

$$\begin{array}{cc}& \frac{\mathrm{d}V}{\mathrm{d}t}=-\frac{1}{2}A{\rho }_{\mathrm{atm}}{V}^2\frac{s}{m}\hfill \end{array}$$(19)

$$\begin{array}{cc}& \frac{\mathrm{d}m}{\mathrm{d}t}=-\frac{1}{2}B{\rho }_{\mathrm{atm}}{V}^{3}s\hfill \end{array}$$(20)

$$\begin{array}{cc}& s=m^{\mu }.\hfill \end{array}$$(21)

With c_d= drag coefficient, c_h= heat-transfer coefficient, H= destruction, enthalpy by kg, ρ_atm= gas density, M= meteoroid mass, M_e= pre-entry mass, S= middle section area, S_e= pre-entry middle section area, μ= shape change coefficient, we have:

$$\begin{array}{c}\hfill A=\frac{c_{\mathrm{d}}{S}_{\mathrm{e}}}{M_{\mathrm{e}}},B=\frac{c_{\mathrm{h}}{S}_{\mathrm{e}}}{H{M}_{\mathrm{e}}},m=\frac{M}{M_{\mathrm{e}}},\mathrm{and}s=\frac{S}{S_{\mathrm{e}}}·\end{array}$$

These four parameters can be reduced to three if only the deceleration of the meteoroid is considered, and not its mass loss (Sansom et al. 2017). After Turchak & Gritsevich (2014) Eqs. (19) and (20) can be rewritten as follows:

$$\begin{array}{cc}\hfill \frac{\mathrm{d}V}{\mathrm{d}t}& =-\frac{1}{2}A{\rho }_{\mathrm{atm}}{V}^2\exp \left(\frac{B(1-\mu )}{2A}(V^2-V_{\mathrm{e}}^2)\right)\hfill \end{array}$$(22)

$$\begin{array}{cc}\hfill m& =\exp \left(\frac{B}{2A}(V^2-V_{\mathrm{e}}^2)\right).\hfill \end{array}$$(23)

The deceleration of a meteoroid is then described by only three parameters: V_e, A, and B(1 − μ). As a result, only B(1 − μ) can be determined by observing the evolution of the velocity of a meteoroid, and not B and μ individually. This new equation system gives us an explicit expression of the normalised mass m as a function of velocity. So the final mass (Eq. (24)) can be computed.

$$\begin{array}{c}\hfill m_{\mathrm{final}}=\exp (-\frac{B{V}_{\mathrm{e}}^2}{2A})·\end{array}$$(24)

These three parameters can be converted into the Turchak & Gritsevich (2014) system:

$$\begin{array}{c}\hfill V_{\mathrm{e}}\alpha =\frac{{\rho }_0{h}_{0}A}{2sin(\gamma )}\beta =(1-\mu )\frac{B{V}_{\mathrm{e}}^2}{2A}\end{array}$$

with ρ₀ the gas density at sea level, h₀ the scale height and γ the angle between horizon and trajectory.

We computed confidence intervals in the three-dimension parameter space using the least squares method. The atmospheric gas density ρ_atm was obtained from the empirical NRLMSISE-00 model. The pertinence of this atmospheric model for meteor science is discussed by Lyytinen & Gritsevich (2016). Results are presented in Table 3. According to this deceleration model, there is about 300 m s⁻¹ difference between the velocity at the top of the atmosphere and the first observed velocity. This illustrates the fact that a direct measurement of the meteor velocity cannot be used to estimate its pre-entry velocity (Egal 2017).

It is not possible to measure the pre-entry mass of the meteoroid, we can only estimate it by making some assumptions. If we assume that the meteoroid is spherical, that its drag coefficient c_d ≃ 1, and that its density is ∼3.5 g cm⁻³, we find a pre-entry mass of 17 g.

5.5. Orbit

In order to propagate the geometrical and velocity uncertainties of our calculation of orbital elements, we used a Monte-Carlo method. We drew random meteoroid positions and directions from the 2σ confidence interval presented in Sect. 5.2. Their velocities were uniformly drawn between 18.842 and 18.968 km s⁻¹ (cf. Sect. 5.4 and Table 3). The results we obtained are shown in Fig. 10 and Table 4.

Fig. 10. Meteoroid orbits (black) from the fireball of 2017-08-04 UT 00:06. Initial position and direction, and initial velocity are uniformly drawn from the 2σ confidence intervals of, respectively, trajectory geometry (see Sect. 5.2) and dynamical parameters (see Sect. 5.4). The orbits of planets Mercury (brown), Venus (green), Earth (blue), Mars (red), and Jupiter (orange) are also shown.

This meteoroid may have originated from the Alinda asteroid family, which is an important source of near-Earth objects (Milani et al. 1989). The objects in this family are in 3:1 resonance with Jupiter, and they have a semi-major axis close to 2.5 AU. Their eccentricities are modified by gravitational interaction with Jupiter, until a close encounter with an inner planet ejects them from their resonance. Known objects from this group have an eccentricity between 0.3 and 0.75, which is comparable to the results we obtained for this meteoroid.

6. Conclusion

We have designed a new method to compute the confidence intervals of meteor trajectories obtained from fish-eye cameras with video rate acquisition speed. This method is an evolved version of the classical least squares method. It consists in applying a least squares regression on a modified sum of residuals. Each camera involved in an event is assigned a specific weight, which depends on the astrometrical measurement noise, the number of images, and the accuracy of its calibration. Our modified sum of residuals is a cross between two cases. In the first case, we assumed that the cameras are perfectly calibrated, so that each image can be considered as one independent measurement, yielding the sum of residuals described by Eq. (8). In the second case, we assumed that residuals are dominated by artefacts due to the cameras themselves. In this case, the images from one camera are perfectly dependent on one another, so that adding or removing images does not add or remove information. In this latter case, it is not the images but the video itself that constitutes one independent measurement; the sum of residuals is described by Eq. (9).

For the FRIPON observation programme, which manages more than 150 cameras and covers a large area, it is impossible to have a perfect calibration of each camera of the network. As a result, assuming each image to be independent from the others as in Eq. (8) may lead to an underestimation of the uncertainties in our trajectories. Conversely, assuming that all the images from a station are perfectly dependent on one another, as in Eq. (9), may result in under-using the available information, and hence overestimating the uncertainties. Our new method is a compromise between these two cases, which are combined to construct our modified sum of residuals: Eq. (12). This allows us to compute realistic confidence intervals for meteor trajectories, even when videos from cameras with different precision and properties are combined. Our results hence have to be looked at with a critical eye, as our method is not derived from a rigorous mathematical demonstration in a general case (which will probably be impossible, due to the complexity of our system), but from a mathematical interpolation between two simpler cases.

As an example of our new method, we present the trajectory of a meteor fallen on August 4, 2017 at 00h06 UT. This meteor was observed by three cameras in Lignan, Sutrieu, and St Lupicin. For the Lignan camera, the non-efficiency of our distortion model is clearly visible, unlike for the other two. The best trajectory we computed reflects this property: the residuals are consistent with Gaussian noise for the Sutrieu and St Lupicin cameras, but not for the Lignan one. In addition, the confidence intervals we computed also reflect our knowledge of the cameras. Our data are less constrained by the observations from the Lignan camera than by data from the other two. This shows that our new method is now ready to be included in the FRIPON pipeline to automatically provide an efficient analysis of meteors observed each night.

Last, our analysis demonstrates that the measurement accuracy for the majority of FRIPON network stations is better than 1 arcmin or 0.1 pixel after implementation of the proper calibration procedures as described here, and that better precision can be expected for bright fireballs. Our method would yield even better results with a more evolved hardware as the FRIPON cameras are limited by the under-sampling of the PSF.

Acknowledgments

We want to thank all FRIPON (www.fripon.org) and Vigie-Ciel (www.vigie-ciel.org) teams for fostering and maintaining the fireball network, especially, G. Quitté, L. Jahan (Toulouse), E. Ravier (Sutrieu), and M. Schieffer (Saint-Lupicin). FRIPON was funded by ANR grant N.13-BS05-0009-03. Vigie-Ciel is coordinated by A. Steinhausser. It contributed to FRIPON network funding through the 65 Millions d’Observateurs project, run by the Muséum National d’Histoire Naturelle and funded by the French Investissements d’Avenir FRIPON data are hosted by Service Informatique Pythéas (SIP) – OSU Institut Pythéas – UMS 3470 (CNRS and Aix-Marseille Université) and at Integrated Data and Operation. Center (IDOC) (Support for IDOC is provided by CNRS and CNES). Database and interface are designed by CeSAM (Centre de donnéeS Astrophysiques de Marseille), LAM (Laboratoire d’Astrophysique de Marseille) UMR 7326.

References

All Tables

All Figures

	Fig. 1. Zenith distance of stars as a function of the distance R to the optical centre. Comparison of recorded stars (in black) with the equidistant model (in red).
In the text

	Fig. 2. Residuals in elevation, using a nine degree odd polynomial function Eq. (3) to model the distortion. The standard deviation is 2.12 arcmin. Toulouse fish-eye camera.
In the text

	Fig. 3. Residuals in azimuth for the Brussels FRIPON camera fitted with a symmetrical distortion function (top) and with an asymmetrical correction (bottom). The standard deviation is 4.76′ without the asymmetrical correction. It goes down to 3.93′ with the asymmetrical correction.
In the text

Fig. 4. Superposition of all calibration stars observed during one night used to compute distortion parameters for the St Lupicin, Sutrieu, and Lignan cameras (from top to bottom). X and Y are the positions of the stars in the images (pixel). Errors between the distortion model and measurements are represented by a colour scale (arcmin). The positions of the meteor of 2017-08-04 UT 00:06 are shown in dark blue.

In the text

	Fig. 5. Residuals of the best trajectory in arcmin (top) and metres (bottom) from cameras of Lignan (red), St Lupicin (green), and Sutrieu (blue). Meteor of 2017-08-04 UT 00:06.
In the text

	Fig. 6. Top: lines of vision from cameras to the meteor of 2017-08-04 UT 00:06, projected in a plane normal to the best trajectory found. This trajectory crosses the plane at coordinates (0, 0). Bottom: zoom on the centre showing random trajectories drawn within 1σ confidence interval (blue) and in 2σ confidence interval (green).
In the text

	Fig. 7. Absolute magnitude as function of altitude for the fireball of 2017-08-04 UT 00:06 seen by the Lignan (red), St Lupicin (green), and Sutrieu (blue) cameras. The absolute magnitude of a fireball is its magnitude at 100 km.
In the text

	Fig. 8. Altitude of the fireball of 2017-08-04 UT 00:06 as a function of time, seen by the Lignan (red), St Lupicin (green), and Sutrieu (blue) cameras. Time starts at 2017-08-04 UT 00:06:18.159.
In the text

	Fig. 9. Velocity as function of altitude for the fireball of 2017-08-04 UT 00:06 seen by the Lignan (red), St Lupicin (green), and Sutrieu (blue) cameras. Dots are point to point measurements. Lines are linear approximation on groups of five points.
In the text

	Fig. 10. Meteoroid orbits (black) from the fireball of 2017-08-04 UT 00:06. Initial position and direction, and initial velocity are uniformly drawn from the 2σ confidence intervals of, respectively, trajectory geometry (see Sect. 5.2) and dynamical parameters (see Sect. 5.4). The orbits of planets Mercury (brown), Venus (green), Earth (blue), Mars (red), and Jupiter (orange) are also shown.
In the text