Issue 
A&A
Volume 626, June 2019



Article Number  A105  
Number of page(s)  10  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/201935580  
Published online  20 June 2019 
Astrometric calibration for allsky cameras revisited
^{1}
INAF – Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese, Italy
email: dario.barghini@inaf.it
^{2}
Università degli Studi di Torino – Dipartimento di Fisica, Via Pietro Giuria 1, 10125 Torino, Italy
^{3}
Osservatorio Astronomico della Regione Autonoma della Val d’Aosta, Lignan 39, 11020 Nus AO, Italy
Received:
29
March
2019
Accepted:
17
May
2019
Context. Several fireball networks deploy allsky cameras for the observation of bright meteors and bolides. Because the field is heavily distorted, a dedicated astrometric reduction is needed. A precise computation of the astrometric solution is essential to determine reliable orbital elements of the parent body and to recover possible fragments on ground.
Aims. The purpose of this article is to assess the astrometric performances of this type of instruments, which is characterized by a wide field of view and small apertures. The currently available parametric models show a high level of complexity and generally suffer from parameter crosstalk and local minimum confinement if the initial estimates are not precisely provided. We address these issues here and propose a solution by adopting a new explicit parametrisation.
Methods. The mismatch between the optical axis and the local zenith direction requires a geometric description that includes two centres of symmetry that lie very close to each other on the focal plane, causing an unreliable estimate of the related parameters. The introduction of new physical coordinates overcomes these issues, allowing a direct and independent estimation. We assessed the performances of different centroiding algorithms in the experimental conditions of an undersampled point spread function of reference stars and saturated bolides on video records. We implemented the algorithm for an automatic identification of bright sources on calibration frames and subsequent correlation with catalogue positions using astrometric projections of increasing complexity.
Results. The algorithm and the new parametrisation of the astrometric solution are tested against real data from the PRISMA Italian fireball network and ensure good convergence properties for all cameras we tested so far. By processing astrometric data with a few months’ statistics, we can achieve a random projection indeterminacy of the order of 10 arcsec, which is negligible with respect to single measurement errors on the bolide position.
Key words: meteorites, meteors, meteoroids / astrometry / methods: data analysis / techniques: image processing
© ESO 2019
1. Introduction
The recovery of fresh meteorite samples and the determination of the orbital elements of their parent bodies can provide useful and relatively lowcost information about the formation and evolution of the solar system (Lauretta & McSween 2006, and references therein). A systematic monitoring of the sky can improve the efficiency in the recovery of samples that survived the impact with our atmosphere and fell to the ground. Collecting new meteorites is also fundamental to improve their classification scheme (Grady et al. 2014). To these purposes, many networks are already operating worldwide, for example, the European Fireball Network (Spurný et al. 2017), the Southern Ontario AllSky Meteor Network (Weryk et al. 2008), the Desert Fireball Network (Devillepoix et al. 2019), and the recently launched FRIPON (Colas et al. 2015) and PRISMA (Gardiol et al. 2016) networks.
To achieve the goal of a successful recovery, an accurate and wellconfined ascertainment of the strewn field is desirable. The first major driver is computing a precise astrometric solution. For an allsky camera, a dedicated astrometric reduction is needed because the optical distortions are strong. This is obtained using a parametric description of the altazimuth coordinates as a function of physical coordinates onto the focal plane (Ceplecha 1987; Borovička 1992; Borovička et al. 1995).
This paper is organised as follows. In Sect. 2 we outline the current status for the astrometry of an allsky camera device. Section 3 gives a brief overview of the technical specification and data format of the PRISMA network. In Sect. 4 we describe the implemented algorithm to determine the associations between bright source positions on the CCD and reference star catalogue, and we discuss the main weaknesses in the current astrometric model. In Sect. 5 we provide a new explicit parametrisation that solves the highlighted issues. In Sect. 6 we discuss the astrometric performances of our algorithm, taking into account effects of point spread function (PSF) undersampling on calibration frames and saturation on bright bolides in video records. Finally, we draw our conclusions in Sect. 7.
2. Current status
The definition of an absolute astrometric model for an allsky camera was first addressed by Ceplecha (1987), who empirically deduced a parametric description that takes into account the heavy optical distortion in the zenith distance direction approaching the horizon:
$$\begin{array}{c}\hfill z=U+Vr+S{e}^{\mathit{Dr}},\end{array}$$(1)
whereas the projection in the azimuth direction is spherical:
$$\begin{array}{c}\hfill a={a}_{0}+\mathrm{atan}\left(\frac{y{y}_{C}}{x{x}_{C}}\right)\xb7\end{array}$$(2)
In Eq. (1), V is the linear plate scale, S and D describe the radial distortion, and r is the distance from (x_{C}, y_{C}) calculated as $r=\sqrt{{(x{x}_{C})}^{2}+{(\mathit{y}{\mathit{y}}_{C})}^{2}}$. After Ceplecha (1987), refinements of this model have been proposed mainly by Borovička (1992) and Borovička et al. (1995). When the optical axis is not perfectly aligned with the local zenith direction, Eqs. (1–2) refer to projections coordinates (b, u) instead of to the horizontal celestial coordinates (a, z). If (E, ϵ) are the azimuth and zenith distance, respectively, of the optical centre O = (x_{O}, y_{O}) with respect to the zenith direction Z = (x_{z}, y_{z}) on the CCD, then (b, u) and (a, z) are related through a translation in spherical coordinates as
$$\begin{array}{c}\hfill \{\begin{array}{c}sin(aE)={\displaystyle \frac{sinbsinu}{sinz}}\hfill \\ cosz=cosucos\u03f5sinucosbsin\u03f5,\hfill \end{array}\end{array}$$(3)
to be combined with the expression for the azimuth cosine:
$$\begin{array}{c}\hfill cos(aE)=\frac{cosucoszcos\u03f5}{sinzsin\u03f5}\xb7\end{array}$$(4)
The explicit form of the transformation for the projection coordinates (b, u) can be given as
$$\begin{array}{c}\hfill \{\begin{array}{c}b={a}_{0}E+\mathrm{atan}\left({\displaystyle \frac{y{y}_{O}}{x{x}_{O}}}\right)\hfill \\ u=Vr+S({e}^{\mathit{Dr}}1),\hfill \end{array}\end{array}$$(5)
where r is calculated from the optical centre as in Eq. (1). Morevoer, Borovička (1992) theoretically deduced that U = −S. Therefore, the explicit transformations to compute (a, z) from (b, u) are
$$\begin{array}{c}\hfill \{\begin{array}{c}a=E+\mathrm{atan}\left({\displaystyle \frac{sinbsinu}{cosbsinucos\u03f5+cosusin\u03f5}}\right)\hfill \\ z=\mathrm{acos}(cosucos\u03f5cosbsinusin\u03f5).\hfill \end{array}\end{array}$$(6)
For our purpose, it would be highly desirable for any astrometric model to be invertible, in order to allow the computation of the expected position of a known source onto the focal plane. The inversion of Eq. (1) with respect to r is not trivial but exists, and it can be derived by means of the Lambert Ω function:
$$\begin{array}{c}\hfill r=\frac{zU}{V}\frac{1}{D}\mathrm{\Omega}[\frac{\mathit{DS}}{V}exp\left(D\frac{zU}{V}\right)]\xb7\end{array}$$(7)
Equation (7) has different and simpler forms when one or more parameters of V, S, D are zero. Borovička et al. (1995) also considered other minor corrections, mainly due to (1) heavier radial distortion for z > 75° and (2) optical plate misalignment with respect to the local horizon plane. These further effects are taken into account and are specifically addressed in Sect. 6.4 in a different way. Therefore, considering Eqs. (5)–(6), the resulting astrometric model consists of M = 8 parameters (a_{0}, x_{O}, y_{O}, E, ϵ, V, S, D) to be estimated. Borovička (1992) reported some issues and advice to manage the convergence of the estimation algorithm. Other authors reported similar issues in the determination of the Borovička model parameters; Bannister et al. (2013) proposed a polynomial representation for u rather than the exponential of Eq. (1), for the NMSU SkySentinel^{1} cameras. Equations (5)–(6) are strongly nonlinear, and some parameters are not independent from each other.
3. PRISMA network
The PRISMA^{2} project was born with the aim of realising a network for detecting meteors and fireball events on the Italian territory. PRISMA stands for “Prima Rete Italiana per la Sorveglianza sistematica di Meteore e Atmosfera”, which means “First Italian Network for Meteors and Atmosphere systematic Surveillance”. PRISMA currently uses the same technology as the FRIPON^{3} network, with a grid spacing of about 60–100 km, which is the optimal distance between two stations for triangulation. As a wide field of view (FOV) is desired, a CCD camera (6 mm diagonal, 1296 × 966 px) is equipped with a short focal lens objective (1.25 mm). The camera is controlled by the FREETURE opensource software (Audureau et al. 2014), installed on a Linux operating miniPC system connected via LAN to the camera itself. In order to detect the meteor trail with a suitable sampling rate, the camera operates at 30 fps and triggers using a simple frame difference method. Detection frames are saved locally and postprocessed by temporal crosscorrelation with respect to other nodes of the network. However, these data are not useful for performing the astrometric calibration of the instrument because no stars are detected with an exposure of 1/30 s. To this purpose, each 10 min, FREETURE performs a 5 s exposure, resulting as a set of equally timespaced images (named captures, following FREETURE nomenclature), which are used to compute the astrometric solution. An example of a capture is presented in Fig. 1 for the Pino Torinese station. These data, acquired as FITS files, are organised in directories by station and month of acquisition; the dataset is sorted by Julian date (JD) and divided into Julian days, so that images from the same night are listed in the same subset.
Fig. 1. Example of capture for the Pino Torinese camera acquired on 6 January 2017, 00:16:12 UT (5 s exposure), with the results of the catalogue association algorithm: red circles enclose the sources that are associated with catalogue positions (N_{*} = 191). The image is oriented as described in Sect. 4.1. 
4. Algorithm implementation
To provide the astrometric calibration of PRISMA cameras, we developed a dedicated reduction software, written in IDL^{4} v8.7. In our procedures, we made an extensive use of the IDL Astronomy User’s Library^{5} (Landsman 1993), which contains very many routines for the analysis of astronomical data, and of the nonlinear fitting routines (MPFIT) included in the MarkwardtIDL Library^{6} (Markwardt 2009), which implements the Levenberg–Marquardt optimisation algorithm (Levenberg 1944; Marquardt 1963; Moré 1978) applied to chisquare minimisation.
4.1. Source identification and catalogue correlation
The first task to be performed is the identification of the bright sources in each frame and the correlation of these sources with a reference catalogue, that is, the definition of the list of associations (x, y)→(a_{cat}, z_{cat}). The reference astrometric catalogue is built through a query to the SIMBAD astronomical database (Wenger et al. 2000). The HIPPARCOS (Perryman et al. 1997) and Tycho catalogues (Høg et al. 1997) are the main sources of stellar positions in SIMBAD; in particular, the Tycho2 catalogue (Høg et al. 2000) provides equatorial coordinates of stars with V_{T} < +9 with standard errors of about 7 mas. At most a few hundred stars are detectable with optimal conditions of sky brightness (corresponding to a limiting magnitude of about +4.5), therefore the definition of the complete astrometric solution cannot be performed at this level because the sky vault is poorly sampled. For this same reason, we do not need an accurate astrometric solution to build the list of associations because the sky vault is not too crowded with stars. We therefore use the following simplified astrometric model:
$$\begin{array}{c}\hfill \{\begin{array}{c}a={a}_{0}+\mathrm{atan}\left({\displaystyle \frac{y{y}_{C}}{x{x}_{C}}}\right)\hfill \\ z=F\mathrm{asin}\left({\displaystyle \frac{r}{R}}\right)\xb7\hfill \end{array}\end{array}$$(8)
Equation (8) only depends upon five parameters, taking into account only one centre of symmetry C = (x_{C}, y_{C}), and provides residuals in a − a_{cat} sin z_{cat} and z − z_{cat} as large as 2°. The five parameters (a_{0}, x_{C}, y_{C}, F, R) can be easily estimated starting from the frame size and the known approximate radial plate scale of about 10′/pix. Thus, convergence can easily be achieved even with a small number of stars per frame. For the Pino Torinese station data (and similarly in all the other tested cameras), the residual minimisation algorithm returns values of F = 1.9 ± 0.1 and R = 650 ± 20 px. These values are consistent with an equisolid projection, for which nominal values are F = 2 and R = 2f/d_{px} = 667 px, being f = 1.25 mm the focal length of PRISMA cameras and d_{px} = 3.75 μm the pixel size of the CCD. Furthermore, it is easily shown that F = V ⋅ R. This implies that PRISMA cameras adhere to equisolid projection within the tolerance of 2°. Because the agreement is not exact, F and R are not fixed to their nominal values in order to provide the minimum residuals for stars associations. We prefer this description for z over the one proposed by Bannister et al. (2013) because with only two parameters, it allows star associations up to the local horizon while a secondorder polynomial barely reaches z ≃ 70°, requiring higher orders to achieve the same performance.
For each data subset from one Julian Day, the main steps performed by the algorithm are listed below.

Images are rotated to match the reference frame with N direction corresponding to the positive abscissa axis (x) and E direction to the positive ordinate axis (y).

For each capture, Sun and Moon ephemerides are computed from the known JD of the frame and the geographic latitude and longitude (λ_{z}, ϕ_{z}) of the station. Frames with z_{⊙} > − 12° are excluded (nautical dawn), and if the Moon is above the horizon and its phase is greater then a given limit (i.e., ϕ_{☾} > 0.5), a circular mask of 150 px radius (∼25°) is applied to minimize the straylight effect.

From all the captures of the same night, a median flat frame is computed in order to normalise each frame with respect to the mean brightness spatial distribution. This is useful to perform automatic identification of bright sources because the sky background is scaled to a value of unity and the detection limits can be set independently of the mean sky brightness of the single frame. This is also done to treat stationary sources of noise, while nonstationary sources are cancelled through the application of a mask. Even if the camera specifications provided by the manufacturer report a negligible spatial noise, this operation provides a partial correction for differential pixel gain and offset because no flatfielding calibration is available during operations.

When the list of bright sources is obtained, we consider the reference star catalogue and compute expected position of each star onto the frame. To this purpose, inverse relations of the astrometric projection of Eqs. (8) can be simply provided as
$$\begin{array}{c}\hfill \{\begin{array}{c}{x}_{\mathrm{cat}}={x}_{C}+{r}_{\mathrm{cat}}cos({a}_{\mathrm{cat}}{a}_{0})\hfill \\ {y}_{\mathrm{cat}}={y}_{C}+{r}_{\mathrm{cat}}sin({a}_{\mathrm{cat}}{a}_{0}),\hfill \end{array}\end{array}$$(9)
where r_{cat} = R sin(z_{cat}/F). In this way, we can spatially match the two lists of (x, y) registered on the frame and the expected position (x_{cat}, y_{cat}). This algorithm is implemented in an iterative way. Starting from very bright stars, the magnitude limit in the catalogue is progressively increased as the indeterminacy on the model parameters becomes smaller. The correlation radius is also reduced. Iterations stop when the number of identified stars and the projection parameters values become stable with respect to the imposed tolerance limits.
In Fig. 2 we present the resulting statistics about daily sky vault sampling for the Pino Torinese station in 2017–2018 (323 astrometric nights). In Fig. 2a the visible minima at 80° and 130° are caused by antennas and the flag located nearby the station. Other features are also evident, for instance, related to the detection of αUMi (peak at a ≃ 360° in Fig. 2a and at z ≃ 45° in Fig. 2b) and the lack of bright closely circumpolar stars. The azimuth coverage shows a slight decreasing trend (likely due to the lights of Turin) with a mean value of nine stars/(day ⋅ deg^{2}). The zenith distance distribution shows a more evident trend, with a decreasing efficiency in detection of stars closer to the horizon. This is mainly due to artificial light pollution, together with the radial counting efficiency loss of the optical system. As a consequence, no stars are detected in the last 10° above the horizon for the Pino Torinese station. This limit is generally valid for PRISMA cameras located nearby populated cities, for which the horizon is polluted by incoming artificial lights. For cameras located in darker places, stars are detected down to the natural horizon of the site.
Fig. 2. Daily mean number of stars per square degree as a function of (a) azimuth and (b) zenith distance for the Pino Torinese station, derived from 2017–2018 statistics. 
4.2. Determining the astrometric solution
The algorithm described in Sect. 4.1 provides the list of associations and the list of parameter values for the simplified model of Eq. (8). Figure 3 shows an example of the distribution of residuals between observed and computed horizontal coordinates. Both plots show significant systematic deviations; in Fig. 3a azimuth residuals are distributed on a family of fixed phase sinusoids, with amplitude depending on the zenith distance. This systematic error in the azimuth direction is reflected also in the exceeding spread of zenith distance residuals, shown in Fig. 3b, with respect to the expected divergence at high values due to the simplified distortion model.
Fig. 3. Fit residuals for data of the night on 5 January 2017 of the Pino Torinese station: (a) azimuth and (b) zenith distance residuals, obtained using the simplified model of Eq. (8). 
The next step is the derivation of the complete astrometric solution, that is, the value determination of the projection parameters (a_{0}, x_{O}, y_{O}, E, ϵ, V, S, D). Any optimisation algorithm requires a proper estimate of the starting point to ensure convergence of the output parameters, especially for complex and nonlinear mathematical relations such as those of Eqs. (5)–(6). With respect to Eq. (8), the parameter a_{0} is the north direction, and its meaning is the same in both models, so that a good estimate is already provided. An estimate of V is given by F/R. This value can be used to compute estimates for (S, D) assuming that z ≃ u and performing an exponential regression with respect to residuals:
$$\begin{array}{c}\hfill \mathrm{\Delta}z=zVr=S({e}^{\mathit{Dr}}1).\end{array}$$(10)
For the remaining parameters (x_{O}, y_{O}, E, ϵ), the main complication is that Eq. (8) assumes only one centre of symmetry C, but we know from Eqs. (5)–(6) that two centres exist, that is, the optical centre O that (b, u) refers to, and the zenith direction Z that (a, z) refers to. Because C does not have a direct physical meaning, we expect that it will lie somewhere in between O and Z. Figure 4 shows the distribution of the positions C computed for the Pino Torinese station, obtained through the available statistics, and outlines that the two distributions of x_{C} and y_{C} values are correlated. The twodimensional histogram shows an elongated distribution with a tilt angle of ∼15°, which corresponds to the azimuth where the residuals are close to zero (see Fig. 3a). In other words, this is the phase of the sinusoid family, with an indeterminacy of ±180°. The reason for this behaviour is clear in Fig. 5, which summarises the geometry of the problem. If C ≡ O, for a star lying at an azimuth a = E, the correction implemented through the translation of (E, ϵ) is null in the azimuth direction, that is, a = b. However, as it is possible to verify a posteriori, C ≠ O even if O, C, and Z lie approximately on the same great circle. Therefore the problem cannot be completely solved by resolving the spherical triangles of Fig. 5 because this can only provide information about the orientation E; Borovička (1992) suggested, for example, an initial estimate of ϵ = 0.01 rad and (x_{O}, y_{O}) = (x_{Z}, y_{Z}). Our approach aims to determine an initial estimate for these parameters specifically tailored to each camera, in order to avoid the correlation between (E, ϵ) and (x_{O}, y_{O}). The explicit form of Eqs. (5)–(6) is given by means of one point O = (x_{O}, y_{O}) and its distance (E, ϵ) with respect to another point Z. These two parameter pairs are clearly mutually dependent as moving the position of O also changes its distance with respect to Z. In the framework of an optimisation algorithm, all the parameters are required to be independent of one another. The focus of the problem is now to determine whether a new parametrisation can be given for the astrometric model in order to remove this dependence.
Fig. 4. Distribution of computed positions of the projection centre (x_{C}, y_{C}) of the simplified model of Eq. (8) derived from 2017–2018 statistics for the Pino Torinese station: contours lines (red) and tilt direction of the distribution (white), corresponding to the azimuth E. 
Fig. 5. Geometrical illustration of the optical centre (O), the zenith direction (Z) on the CCD, and corresponding projection coordinates (b, u) and horizontal coordinates (a, z). 
5. New parametrisation of the astrometric solution
The requirement on the parameters to be mutually independent suggests that (x_{O}, y_{O}, E, ϵ) might be replaced with the new subset (x_{O}, y_{O}, x_{Z}, y_{Z}). A direct benefit is the possibility of independently determining an estimate of the coordinates (x_{Z}, y_{Z}). This can be done, for example, by analysing meridian crossings of many stars. From the catalogue associations list, the path of each star during the night can be extracted, sampled every 10 min, that is, each 2.5° of hour angle. The meridian crossings can be evaluated when the local sideral time is equal to the right ascension α of the star. The result is a set of coordinates (x_{m}, y_{m}) that describes the local meridian projected onto the CCD plate as a function of declination δ. An example is given in Fig. 6a; the interpolation of (δ, x_{m}) and (δ, y_{m}) allows us to retrieve the zenith direction (x_{z}, y_{z}) evaluated for δ equal to the geographic latitude λ_{z} of the observational site, as shown in Figs. 6b,c. At this point, it is possible to also determine an initial estimate for (x_{O}, y_{O}) by varying the position of O around Z and requiring the minimisation of residuals Δa ⋅ sin (z) and Δz. An example of the results of this processing is shown in Fig. 7. The explicit formulation of the new parametrisation is obtained by adding the definition of (E, ϵ) as a function of the other parameters. (E, ϵ) defines the distance of O with respect to Z:
$$\begin{array}{c}\hfill \{\begin{array}{c}E={a}_{0}+\mathrm{atan}\left({\displaystyle \frac{{y}_{O}{y}_{Z}}{{x}_{O}{x}_{Z}}}\right)\hfill \\ \u03f5=V{r}_{\u03f5}+S({e}^{D{r}_{\u03f5}}1),\hfill \end{array}\end{array}$$(11)
Fig. 6. Example of the determination of the local meridian and zenith positions with data of the night on 5 January 2017 of the Pino Torinese station: (a) local meridian (solid red line) and zenith position (blue dot) on the CCD, identified by the interpolation of the stellar trajectories (black lines); (b) declination vs. x pixels of the observed crossings of the local meridian (black dots), spline interpolation (red line) and evaluation of x_{Z} value (blue); (c) same as (b), but for y_{Z} determination. 
Fig. 7. Map of the residuals sum values (χ_{LAV}, see Eq. (12)) for fixed Z = (x_{Z}, y_{Z}) (and the other distortion parameters) and varying O = (x_{O}, y_{O}) as described in Sect. 5 for the night on 5 January 2017 of the Pino Torinese station. An evident minimum is visible at about 12 px distance from Z, indicating the initial guess for O, where χ_{LAV} ≃ 1. 
where ${r}_{\u03f5}=\sqrt{{({x}_{O}{x}_{Z})}^{2}+{({\mathit{y}}_{O}{\mathit{y}}_{Z})}^{2}}$. The method we propose here has the further advantage that it is also applicable to lowresolution systems, where the sensitivity to starting point settings is magnified.
Figures 8a,b show an example for the residual histograms for azimuth and zenith distance obtained with this new parametrisation; both histograms highlight that residuals are distributed around zero, with a full width at halfmaximum (FWHM) Γ_{a} ≃ 4′ and Γ_{z} ≃ 4.5′ respectively. Both distributions are characterised by a positive excess kurtosis of k_{a} = 4.3 and k_{z} = 2.8, respectively, indicating a deviation from the normal distribution towards fatter tails. These figures also show a small fraction (typically ≲1%) of outliers, that is, results of incorrect associations of the correlation algorithm described in Sect. 4; this is an expected drawback of an automated procedure and is mainly due to false positives of the source search algorithm. The nonGaussianity of residuals together with the presence of outliers suggests the introduction of a weightreduction scheme. A possible approach is to modify the L^{2}norm of the ordinary leastsquares (OLS) by adopting the L^{1}norm of the least absolute value (LAV) regression (Narula & Wellington 1982; Dielman 2005), which is an efficient estimator as well, meaning that it provides parameter estimates close to those from OLS:
$$\begin{array}{c}\hfill {\chi}_{\mathrm{LAV}}={\displaystyle \sum _{i=0}^{{N}_{\ast}}}\frac{a({x}_{i},{y}_{i},\mathit{p}){a}_{\mathrm{cat},i}sin({z}_{\mathrm{cat},i})}{{\sigma}_{a}}+{\displaystyle \sum _{i=0}^{{N}_{\ast}}}\frac{z({x}_{i},{y}_{i},\mathit{p}){z}_{\mathrm{cat},i}}{{\sigma}_{z}}\xb7\end{array}$$(12)
Fig. 8. Summary of results for the astrometric reduction of the data for the night on 5 January 2017 for the Pino Torinese station (N_{*} = 8300). For azimuth (left column) and zenith distance (right column), we show (a–b) bi and monodimensional residuals histograms (2D histograms are logstretched); (c–d) distorsion maps from the linear projection; and (e–f) random projection error maps. 
Since LAV estimates have low variance in the case of a heavy tailed residual distribution and are provided with control bias in case of large samples, this method is preferable over OLS minimisation to reduce the effects of vertical outliers. However, L^{1}norm does not protect against the effect of leverage points. As a consequence, it can lead to possible multiple solutions because the breakpoint of LAV estimates is the same as OLS (Rousseeuw & Leroy 2005; Wilcox 2010). Nevertheless, this eventuality is avoided in our case as accurate and independent initial estimates for model parameters are provided, as described above.
Figures 8c,d show the distortion maps, that is, the absolute differences between the complete astrometric model and the model with S = D = 0 and O ≡ Z. The distortion introduced by the nonlinear component can be up to many degrees. Figures 8e,f show the expected uncertainties due to indeterminacy of the projection parameters of Eqs. (5)–(6)–(11) returned by the minimisation algorithm, which is of the order of 10 arcsec and does not exceed the value of 30 arcsec up to z = 80°.
An estimator of the global correlation among the model parameters can be given by the value of the determinant of the inverse correlation matrix d = detR^{−1}, which increases with the correlation degree. The determinant of the correlation matrix is linked to the Mdimensional space volume occupied by the indirect measurements of projection parameters. This volume is equal to unity for an uncorrelated system and goes to zero for a totally correlated system. The new parametrisation provides a substantial reduction of d, computed on a daily basis for the Pino Torinese station data, whose mean value drops from 300 ± 10 to 40 ± 1.
The new parametrisation implemented in the algorithm gives stable results for the whole PRISMA network currently in operation. Figure 9 shows the random projection errors for a and z for all the 33 cameras. Values are normalised to N_{*} = M + 1, as explained in Sect. 6.3. All cameras exhibit similar performances so that it is possible to give a mean PRISMA random projection error, shown as red curves.
Fig. 9. Random projection errors in azimuth (a) and zenith distance (b), renormalised to N_{*} = M + 1 for the sample of 33 cameras tested of the PRISMA network. The σ values are mildly dependent on the nightsky quality of the observational site. The red curves highlight the mean performances of the network. 
6. Astrometric error
The indeterminacy introduced by the astrometric calibration has two main contributions: (1) the random component, driven by the measurement error on calibration sources, and (2) the systematic component due to possible model inadequacy, for example, related to the radial distortion parametrisation. To be able to quantify the first term, the performance of the centroiding algorithm must be considered with respect to the spatial sampling resolution. The occurrence of systematic residuals is addressed through a numerical correction. These results have to be compared with the positioning precision achievable for bolide images on video records.
6.1. Sampling resolution
Considering an effective wavelength of 500 nm and the nominal effective aperture of 3.91 mm, we obtain a size of ∼1.1′ for the diffraction Airy figure, corresponding to an FWHM Γ ≃ 0.5′. The 5 s of integration time on captures introduce a variable smearing up to 2′ at the southern horizon. The Nyquist sampling theorem applied to astrometry requires the PSF to be sampled, in one dimension, with at least 3–4 px. With a pixel size of about 10′, we expect to be in a heavy undersampled regime. This is confirmed by looking at the star images on the CCD, whose signal is often confined to one pixel and almost always up to four adjacent pixels, depending on the position of the star with respect to the CCD lattice. Figure 10a shows some examples of star shapes; in some cases, the signal is spread over more than four pixels, likely because of unquantified optical distortions.
Fig. 10. Examples of point spread functions for (a) calibration stars on capture images and (b) evolution of bolide images on video records (2017/05/30 at 21:09:26 UT detection, Rovigo station). During the bolide transit, fragmentation (tails ejecting from the main PSF) and flares (sudden high brightness) are visible. 
6.2. Centroiding algorithm
One important consequence of the heavy undersampling is related to the accuracy of the centroiding algorithm. The Astronomy User’s Library provides an IDL implementation of the IRAF DAOPHOT procedures (Stetson 1987), and in particular, the FIND source identification and centroiding algorithm. The current version is based on marginal distribution fitting, against the original derivative search algorithm (referred to as DERIV in this work). To assess the performance of this software in the regime of undersampled PSF together with other simple algorithms, we performed tests on synthetic images. A sample of stars was generated with positions that were uniformly distributed with respect to the pixel grid, adding Poisson noise and constant sky background. The signaltonoise ratio was deduced from real data: the mean value is around 25 and the minimum value, driven by the detection cutoff, is about 5. Figure 11 shows the comparison of accuracy for the considered algorithms as a function of fractional position within the pixel. Figure 11a refers to stars with Γ = 2 px, near to the Nyquist limit. The FIND (green line) and DERIV (purple line) algorithms perform well over the whole range, except for a slight bias, positive and negative, respectively, at the pixel borders (but almost confined to the 1σ dispersion band of simulated results). The simple barycentre (SBC; red line) shows the known biases related to the sky background and the computing box relative positioning. The blue line refers to a filtered barycentre (FBC), computed by preserving only pixels with values 3σ over the sky median; this correction introduces subpixel discretisations but reduces the SBC bias, which is almost coincident with the brightest pixel computation (BP; not shown). Figure 11b shows the same comparison for stars with Γ = 0.2 px, of the order of our experimental conditions. As long as the PSF is completely confined to a single pixel, all the considered algorithms are consistent with the BP results. Instead, when the PSF signal spreads over a neighbouring pixel, different behaviours are evident. The FIND algorithm shows a major issue, that is, an increasing bias, that at the pixel border is even higher than 1 px. The DERIV method shows similar features, but with a reduced amplitude. The SBC computation maintains similar bias as the previous case, whereas the correction implemented by the FBC provides the smallest absolute deviation and gives the expected results of null bias when the signal is equally divided between two adjacent pixels. The residual histograms for x, y prove to be nearly uniform regardless of the algorithm. As a final consideration, the random error on the single measurement of a star position can be estimated by means of the standard deviation of the sample residuals and varies between 0.2 and 0.4 px, depending on the algorithm. For both undersampled and wellsampled regimes, FBC shows the best performances in terms of accuracy.
Fig. 11. Comparison of subpixel positioning accuracy for different centroiding algorithms: (a) spatial sampling near to the Nyquist criterium (Γ = 2 px) and (b) undersampling regime (Γ = 0.2 px). 
6.3. Random projection error
The random projection error, that is, the uncertainty on the computed astrometric positions due to the model parameter indeterminacy, can be estimated by means of error propagation. The covariance matrix is provided as an output of the optimisation algorithm. The result is mainly driven by (1) the single measurement error (see Sect. 6.2) and (2) the size of the calibration sample N_{*}. In Sect. 5 we briefly discussed these results for a daily statistics (N_{*} = 8300). To determine the algorithm performances as a function of N_{*}, we computed the mean astrometric error for a ⋅ sin(z) and z for different sample sizes (daily and monthly statistics over 2017–2018). Results are shown in Fig. 12. Values are mildly scattered above the expected square root dependence. For N_{*} = M + 1, being M = 8 the length of the parameters vector, the corresponding σ values are consistent with the standard deviation of the sample (histograms in Figs. 8a,b). On a daily time base, the random projection error is some tens of arcsec and drops below 10 arcsec for a whole month. For the PRISMA cameras tested so far, the projection parameter values can be considered constant over several months because no systematic deviations are evident at this timescale. When astrometric data from a few months are analysed, the expected projection error reaches roughly 1 arcsec.
Fig. 12. Azimuth (a) and zenith distance (b) random projection errors as a function of the calibration sample size (N_{*}). For zenith distance, the data are presented for z < 60° (blue) and z > 60° (red) separately to outline the plate scale degradation close to the horizon, whereas azimuth data (green) do not show any systematic dependence. The points represent the daily and monthly data, and the straight lines show the square root dependence. 
6.4. Systematic residuals
Figure 13 shows an example of systematic residuals for azimuth and zenith distance projections compared to the random error component. This result is already corrected for the optical plate misalignment, requiring two more parameters: amplitude and phase of the sinusoidal correction (Borovička et al. 1995) with a period of 360°. However, as shown in Fig. 13a, a further contribution in azimuth is still clearly visible in this picture with a period of 180°, but often aperiodic, and modulated in amplitude. This could be related to optical aberrations affecting the fisheye lens and/or mechanical stresses due to the holding system (not given to be symmetrical in azimuth). Figure 13b shows the zenith distance residuals, and systematic effects are also visible in this case. The model of Eq. (5), introduced by Borovička (1992) and also referred to as fisheye transform (FET) in the literature, describes the radial distortion within the pixel accuracy for PRISMA cameras, given the plate scale degradation from 10′/px at zenith to almost 20′/px at the horizon. Borovička et al. (1995) introduced an additional exponential correction of argument r^{2} to provide a correction for this effect and improved the zenith distance accuracy for z > 75°. Because our cameras are in general able to detect stars up to z ≲ 80°, with some exceptions in the case of very dark observational sites, we prefer a numerical rather than a parametric approach. This is implemented for both a and z by tabulating systematic residuals, as from Fig. 13, and adding them over the analytical projection results from Eqs. (5)–(6)–(11). Hughes et al. (2010) considered several analytical lens models, among which the FET, highlighting that mostly due to manufacturing tolerances, an ad hoc correction is always required regardless of the nominal distortion model that is applied. The application of these models to our data confirms this conclusion, showing a significant systematic deviation in all cases, as we showed for the equisolid projection in Sect. 4.1. The FET distortion model provide the smallest systematic residuals and preserves a small number of parameters. It can be verified that an error in the geographical position of the camera of ∼30 m causes systematic residuals of the order of 1 arcsec, so that this contribution can be considered negligible in the context of this work.
Fig. 13. Azimuth (a) and zenith distance (b) systematic residuals (black line and dots) deduced from the analysis of the January 2017 data for the Pino Torinese station. The random projection error is also reported for comparison (green: 1σ, orange: 2σ, and red: 3σ). 
6.5. Final astrometric error
The results discussed in the previous sections about the astrometric error are to be compared to uncertainties in the determination of the bolide position on the CCD. In particular, the main requirement on the astrometric model is that the introduced indeterminacy (Sect. 6.3) has to be negligible with respect to the single measurement error of the bolide position. Figure 10b shows an example of bolide evolution during the visible trajectory through the atmosphere. Compared to the discussion of Sects. 6.1–6.2 for the PSF of stars, substantial differences are obvious. The magnitude of the meteor during the atmospheric phase is highly variable. Almost all bolide detections are triggered during the initial stage of increasing light emission. While the atmospheric friction slows the meteoroid down, the light curve reaches one or more maxima, depending on the physical evolution, and fragmentation may occur causing flares. When the bolide velocity drops below 3−4 km s^{−1}, the light emission suddenly ceases; see Pecina & Koten (2009), CampbellBrown et al. (2013) for a comprehensive overview of atmospheric dynamic and ablation of the meteoroid. This is reflected in a highly dynamic intensity range, from very faint to very bright and even saturated images for the brightest bolides. The PSF size is also variable and can reach dimensions of several pixels, therefore issues related to undersampling become less relevant.
As for the case of Sect. 6.2, we performed numerical simulations to evaluate the performances of different centroiding algorithms with respect to these experimental conditions. We dropped the SBC computation for its poor performances with respect to FBC and replaced it with the PSF fitting, which was not suitable for the case of the undersampled reference stars. Figure 14 shows the results for both simulated and real images as a function of the normalised PSF height h, defined as the counts of the central pixel of the fitted PSF divided by the saturation capacity of the CCD. For simulated images (with Γ = 2 px), we report the sample standard deviation, obtained by varying the relative subpixel position of the bolide. On real images, the dispersion is computed with respect to a mean trajectory that is obtained as spline smoothing (Reinsch 1967). PSF fitting and FBC computation prove to be the bestperforming algorithms in both simulated and real data. The FIND results show significantly higher dispersion values, increasing for h ≳ 1. DERIV exhibits general poorer performances over the whole range. This comparison is consistent between Figs. 14a and b: the PSF and FBC simulated dispersion drops continuously for increasing signaltonoise ratio, whereas FIND and DERIV reach a minimum value for h ≲ 1 and show opposite behaviours for h > 1. On real data, however, the error does not go below 0.5′ (∼1/20 px) even for PSF and FBC methods. The limiting factor can be related to the presence of strong asymmetries in the bolide images (see Fig. 10b) due to either fragmentation, flares, and PSF elongation caused by high angular speed of its image on the CCD. Therefore, when PSF and FBC results are considered, the performance degradation on saturated images is not an limiting issue at least up to h ≤ 4.
Fig. 14. Positioning precision as a function of the normalized PSF height for (a) simulated PSF and (b) real bolide images, i.e., 20 events for a total of 273 video records of PRISMA network. The simulation of bolide images is provided with Γ = 2 px (mean value obtained on the respective bolide data). The dashed line highlights the saturation level of h = 1. 
The value of 0.5′ on the single image can be considered as a lower limit for the measurement accuracy on the bolide position, and it is the reference value to assess the quality of the astrometric solution (Sects. 6.3–6.4). From Fig. 12 and accordingly, for the mean random projection errors from Fig. 9, this requires a calibration sample of the order of 10^{4} star associations. In general, a monthly statistics is sufficient to provide quality astrometry.
7. Conclusions
The astrometric calibration of an allsky camera requires an ad hoc treatment that differs from the classical rotation and plate scaling because of heavy field distortions. The absolute astrometric solution is usually derived from a large sample of calibration images. In this work, we describe a new parametrisation of the most widely used model that solves the issues of local minimum confinement in the parameter value determination that is caused by the strong correlation between parameters. These effects are magnified in the case of instruments with low spatial resolution, such as digital cameras equipped with a small CCD. The introduction of the zenith coordinates in instrumental units as explicit parameters allows an independent and reliable estimation by the analysis of meridian crossings, providing also an improvement of one order of magnitude in the correlation degree and ensuring a stable and robust determination of the astrometric solution. This new method has been tested, implemented, and validated on the PRISMA network cameras. The astrometric solution is obtained for each camera by an automated process that enables discarding incorrect source associations due to, for example, moon brightness or artificial lights.
PRISMA cameras are usually able to detect stars up to z ≲ 80° due to effects of light pollution at the horizon, except for some stations that are located in dark sites and far away from populated cities. We found systematic deviations in both a and z, similar to those reported by Borovička et al. (1995). We preferred a numerical correction for these effects mostly because of our higher altitude limit and also for their subpixel magnitude. As a consequence of the lack of reference stars for z ≳ 80°, the astrometric reduction of bolides very close to the horizon is biased for some cameras. This issue is partially overcome by the high node density adopted in recently launched meteor networks, such as FRIPON and PRISMA. Cameras that record the bolide close to the horizon provide by themselves less accurate data for triangulation purposes because the distance between object and observer is greater.
The effects of PSF undersampling on reference frames and saturation on bolide images were evaluated and tested against different centroiding algorithms. In all cases the FBC showed the best performances. To achieve a negligible astrometric error with respect to the single frame position error of a meteor, it is necessary to process at least one month of calibration images.
Acknowledgments
PRISMA is the Italian Network for Systematic surveillance of Meteors and Atmosphere. It is a collaboration initiated and coordinated by the Italian National Institute for Astrophysics (INAF) that counts members among research institutes, associations and schools. The complete list of PRISMA members is available here: http://www.prisma.inaf.it. PRISMA was partially funded by a 2016 Research and Education grant from Fondazione Cassa di Risparmio di Torino and by a 2016 grant from Fondazione Agostino De Mari (Savona). The initial FRIPON hardware and software has been developed by the FRIPONFrance core team under a French ANR grant (20142018).
References
 Audureau, Y., Marmo, C., & Bouley, S. 2014, in Proceedings of the 33th International Meteor Conference, Giron, France, 18–21 September 2014, eds. J. L. Rault, & P. Roggemans, 39 [Google Scholar]
 Bannister, S. M., Boucheron, L. E., & Voelz, D. G. 2013, PASP, 125, 1108 [NASA ADS] [CrossRef] [Google Scholar]
 Borovička, J. 1992, Publ. Inst. Czech. Acad. Sci., 79, 19 [Google Scholar]
 Borovička, J., Spurny, P., & Keclikova, J. 1995, A&AS, 112, 173 [NASA ADS] [Google Scholar]
 CampbellBrown, M. D., Borovička, J., Brown, P. G., & Stokan, E. 2013, A&A, 557, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ceplecha, Z. 1987, Bull. Astr. Inst. Czechosl., 38, 222 [Google Scholar]
 Colas, F., Zanda, B., & Vaubaillon, J. 2015, in Proceedings of the 34th International Meteor Conference, Mistelbach, Austria, 27–30 August 2015, eds. J.L.Rault, & P. Roggemans, 37 [Google Scholar]
 Devillepoix, H. A. R., Bland, P. A., Sansom, E. K., et al. 2019, MNRAS, 483, 5166 [NASA ADS] [CrossRef] [Google Scholar]
 Dielman, T. E. 2005, J. Stat. Comput. Simul., 75, 263 [CrossRef] [Google Scholar]
 Gardiol, D., Cellino, A., & Di Martino, M. 2016, in Proceedings of the 35th International Meteor Conference, Egmond, The Netherlands, 2–5 June 2016, 76 [Google Scholar]
 Grady, M., Pratesi, G., & Moggi Cecchi, V. 2014, Atlas of Meteorites (Cambridge University Press) [Google Scholar]
 Høg, E., Bässgen, G., Bastian, U., et al. 1997, A&A, 323, L57 [Google Scholar]
 Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27 [Google Scholar]
 Hughes, C., Denny, P., Jones, E., & Glavin, M. 2010, Appl. Opt., 49, 3338 [NASA ADS] [CrossRef] [Google Scholar]
 Landsman, W. B. 1993, in Astronomical Data Analysis Software and Systems II, eds. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, PASPC, 52, 246 [NASA ADS] [Google Scholar]
 Lauretta, D. S., & McSween, H. Y. 2006, Meteorites and the Early Solar System II (University of Arizona Press) [Google Scholar]
 Levenberg, K. 1944, Q. Appl. Math., 2, 164 [Google Scholar]
 Markwardt, C. B. 2009, in Astronomical Data Analysis Software and Systems XVIII, eds. D. A. Bohlender, D. Durand, & P. Dowler, PASPC, 411, 251 [NASA ADS] [Google Scholar]
 Marquardt, D. W. 1963, SIAM J. Appl. Math., 11, 431 [Google Scholar]
 Moré, J. J. 1978, Numerical Analysis (Springer), 105 [Google Scholar]
 Narula, S. C., & Wellington, J. F. 1982, Int. Stat. Rev., 50, 317 [CrossRef] [Google Scholar]
 Pecina, P., & Koten, P. 2009, A&A, 499, 313 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perryman, M. A. C., Lindegren, L., Kovalevsky, J., et al. 1997, A&A, 323, L49 [NASA ADS] [Google Scholar]
 Reinsch, C. H. 1967, Numer. Math., 10, 177 [CrossRef] [Google Scholar]
 Rousseeuw, P. J., & Leroy, A. M. 2005, Robust Regression and Outlier Detection (John Wiley& Sons), 589 [Google Scholar]
 Spurný, P., Borovička, J., Mucke, H., & Svoren, J. 2017, A&A, 605, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stetson, P. B. 1987, PASP, 99, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weryk, R. J., Brown, P. G., Domokos, A., et al. 2008, Earth Moon Planets, 102, 241 [NASA ADS] [CrossRef] [Google Scholar]
 Wilcox, R. R. 2010, Fundamentals of Modern Statistical Methods: Substantially Improving Power and Accuracy (Springer) [CrossRef] [Google Scholar]
All Figures
Fig. 1. Example of capture for the Pino Torinese camera acquired on 6 January 2017, 00:16:12 UT (5 s exposure), with the results of the catalogue association algorithm: red circles enclose the sources that are associated with catalogue positions (N_{*} = 191). The image is oriented as described in Sect. 4.1. 

In the text 
Fig. 2. Daily mean number of stars per square degree as a function of (a) azimuth and (b) zenith distance for the Pino Torinese station, derived from 2017–2018 statistics. 

In the text 
Fig. 3. Fit residuals for data of the night on 5 January 2017 of the Pino Torinese station: (a) azimuth and (b) zenith distance residuals, obtained using the simplified model of Eq. (8). 

In the text 
Fig. 4. Distribution of computed positions of the projection centre (x_{C}, y_{C}) of the simplified model of Eq. (8) derived from 2017–2018 statistics for the Pino Torinese station: contours lines (red) and tilt direction of the distribution (white), corresponding to the azimuth E. 

In the text 
Fig. 5. Geometrical illustration of the optical centre (O), the zenith direction (Z) on the CCD, and corresponding projection coordinates (b, u) and horizontal coordinates (a, z). 

In the text 
Fig. 6. Example of the determination of the local meridian and zenith positions with data of the night on 5 January 2017 of the Pino Torinese station: (a) local meridian (solid red line) and zenith position (blue dot) on the CCD, identified by the interpolation of the stellar trajectories (black lines); (b) declination vs. x pixels of the observed crossings of the local meridian (black dots), spline interpolation (red line) and evaluation of x_{Z} value (blue); (c) same as (b), but for y_{Z} determination. 

In the text 
Fig. 7. Map of the residuals sum values (χ_{LAV}, see Eq. (12)) for fixed Z = (x_{Z}, y_{Z}) (and the other distortion parameters) and varying O = (x_{O}, y_{O}) as described in Sect. 5 for the night on 5 January 2017 of the Pino Torinese station. An evident minimum is visible at about 12 px distance from Z, indicating the initial guess for O, where χ_{LAV} ≃ 1. 

In the text 
Fig. 8. Summary of results for the astrometric reduction of the data for the night on 5 January 2017 for the Pino Torinese station (N_{*} = 8300). For azimuth (left column) and zenith distance (right column), we show (a–b) bi and monodimensional residuals histograms (2D histograms are logstretched); (c–d) distorsion maps from the linear projection; and (e–f) random projection error maps. 

In the text 
Fig. 9. Random projection errors in azimuth (a) and zenith distance (b), renormalised to N_{*} = M + 1 for the sample of 33 cameras tested of the PRISMA network. The σ values are mildly dependent on the nightsky quality of the observational site. The red curves highlight the mean performances of the network. 

In the text 
Fig. 10. Examples of point spread functions for (a) calibration stars on capture images and (b) evolution of bolide images on video records (2017/05/30 at 21:09:26 UT detection, Rovigo station). During the bolide transit, fragmentation (tails ejecting from the main PSF) and flares (sudden high brightness) are visible. 

In the text 
Fig. 11. Comparison of subpixel positioning accuracy for different centroiding algorithms: (a) spatial sampling near to the Nyquist criterium (Γ = 2 px) and (b) undersampling regime (Γ = 0.2 px). 

In the text 
Fig. 12. Azimuth (a) and zenith distance (b) random projection errors as a function of the calibration sample size (N_{*}). For zenith distance, the data are presented for z < 60° (blue) and z > 60° (red) separately to outline the plate scale degradation close to the horizon, whereas azimuth data (green) do not show any systematic dependence. The points represent the daily and monthly data, and the straight lines show the square root dependence. 

In the text 
Fig. 13. Azimuth (a) and zenith distance (b) systematic residuals (black line and dots) deduced from the analysis of the January 2017 data for the Pino Torinese station. The random projection error is also reported for comparison (green: 1σ, orange: 2σ, and red: 3σ). 

In the text 
Fig. 14. Positioning precision as a function of the normalized PSF height for (a) simulated PSF and (b) real bolide images, i.e., 20 events for a total of 273 video records of PRISMA network. The simulation of bolide images is provided with Γ = 2 px (mean value obtained on the respective bolide data). The dashed line highlights the saturation level of h = 1. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.