© 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 low-cost 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 All-Sky 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 well-confined ascertainment of the strewn field is desirable. The first major driver is computing a precise astrometric solution. For an all-sky camera, a dedicated astrometric reduction is needed because the optical distortions are strong. This is obtained using a parametric description of the alt-azimuth 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 all-sky 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 all-sky 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:

(1)

whereas the projection in the azimuth direction is spherical:

(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 . 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

(3)

to be combined with the expression for the azimuth cosine:

(4)

The explicit form of the transformation for the projection coordinates (b, u) can be given as

(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

(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:

(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₀, 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¹ cameras. Equations (5)–(6) are strongly non-linear, and some parameters are not independent from each other.

3. PRISMA network

The PRISMA² 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³ 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 open-source software (Audureau et al. 2014), installed on a Linux operating mini-PC 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 post-processed by temporal cross-correlation 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 time-spaced 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⁴ v8.7. In our procedures, we made an extensive use of the IDL Astronomy User’s Library⁵ (Landsman 1993), which contains very many routines for the analysis of astronomical data, and of the non-linear fitting routines (MPFIT) included in the Markwardt-IDL Library⁶ (Markwardt 2009), which implements the Levenberg–Marquardt optimisation algorithm (Levenberg 1944; Marquardt 1963; Moré 1978) applied to chi-square 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 Tycho-2 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:

(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₀, 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 second-order 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.

1. 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).

2. 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 stray-light effect.

3. 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 non-stationary 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 flat-fielding calibration is available during operations.

4. 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

   (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²). 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₀, 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 non-linear mathematical relations such as those of Eqs. (5)–(6). With respect to Eq. (8), the parameter a₀ 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:

(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 two-dimensional 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 (xC, yC) 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:

(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 xZ value (blue); (c) same as (b), but for yZ determination.

Fig. 7. Map of the residuals sum values (χLAV, see Eq. (12)) for fixed Z = (xZ, yZ) (and the other distortion parameters) and varying O = (xO, yO) 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 . The method we propose here has the further advantage that it is also applicable to low-resolution 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 half-maximum (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 non-Gaussianity of residuals together with the presence of outliers suggests the introduction of a weight-reduction scheme. A possible approach is to modify the L²-norm of the ordinary least-squares (OLS) by adopting the L¹-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:

(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 mono-dimensional residuals histograms (2D histograms are log-stretched); (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¹-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 non-linear 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⁻¹, which increases with the correlation degree. The determinant of the correlation matrix is linked to the M-dimensional 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 night-sky 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 signal-to-noise ratio was deduced from real data: the mean value is around 25 and the minimum value, driven by the detection cut-off, 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 sub-pixel 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 well-sampled regimes, FBC shows the best performances in terms of accuracy.

Fig. 11. Comparison of sub-pixel 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 fish-eye 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 fish-eye 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² 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⁻¹, the light emission suddenly ceases; see Pecina & Koten (2009), Campbell-Brown 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 sub-pixel 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 best-performing 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 signal-to-noise 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⁴ star associations. In general, a monthly statistics is sufficient to provide quality astrometry.

7. Conclusions

The astrometric calibration of an all-sky 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 sub-pixel 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 FRIPON-France core team under a French ANR grant (2014-2018).

References

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 (xC, yC) 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 xZ value (blue); (c) same as (b), but for yZ determination.

In the text

	Fig. 7. Map of the residuals sum values (χLAV, see Eq. (12)) for fixed Z = (xZ, yZ) (and the other distortion parameters) and varying O = (xO, yO) 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 mono-dimensional residuals histograms (2D histograms are log-stretched); (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 night-sky 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 sub-pixel 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