Accurate Ground-based Near-Earth-Asteroid Astrometry Using Synthetic Tracking

Figure 2 provides an overview of the data processing. We start with the raw images and apply calibration data to subtract the dark frame and factor out the flat field responses estimated using twilight images. We then perform preprocessing to remove the cosmic ray events and bad pixel signals. Cosmic ray events are identified as signal spikes above random noise level localized in both temporal and spatial dimension. The least-squares fitting is used to estimate centroid of reference stars. A pre-estimated field distortion correction, modeled as low-order polynomial functions, is applied to the reference stars in pixel coordinates. The telescope pointing and the size of FOV are used to look up star catalog, the Gaia Data Release 1 (DR1) catalog. A planar triangle matching algorithm identifies stars in the field by matching congruent triangles from the field and the catalog with the shape determined by the relative distances between the stars. A set of identified stars enables us to solve for an affine mapping between the pixel coordinate and the position in the plane of sky, and thus the right ascension (R.A.) and declination (decl.). With this mapping, we can covert the pixel coordinate of the asteroid (after the field distortion correction) to R.A. and decl. in sky.

Zoom In Zoom Out Reset image size

For bright objects, it is possible to estimate their centroid using single short-exposure images. We thus can perform NEA astrometry with these centroids from individual frames similar to stellar astrometry because both moving target and the reference stars have compact PSFs in these short-exposure frames. However, in most of the cases, we would need to use the whole data cube to estimate the centroid because when the signal in each short-exposure frame is low, centroiding may not be reliable. Instead, we perform a least-squares fitting of "moving PSF" to the whole data cube using the following cost function:

$$\begin{eqnarray}C({v}_{x},{v}_{y},{x}_{c},{y}_{c},\alpha ,{I}_{0}) & \equiv & \displaystyle \sum _{x,y,t}| I(x,y,t)-\alpha P(x-X(t),\\ & & y-Y(t)){-{I}_{0}| }^{2}w(x,y,t)\end{eqnarray} \tag{ 1 }$$

where we parameterize a moving PSF, $P(x-X(t),y-Y(t))$, with P(x, y) being the PSF function and $(X(t),Y(t))$ represent the location of the object in frame t. We use a Moffat function (Moffat 1969) to model the PSF and determine the parameter using a bright star in the field. The Moffat PSF function includes Gaussian PSF as its special case and fits better than a Gaussian to a seeing limited PSF; but the difference between using a Moffat and Gaussian PSF is much less than 10 mas. w(x, y, t) is a weighting function and N is the total number of frames. To minimize the variance of the estimation, the weight can be chosen to be the inverse of the variance of the measured I(x, y, t), including photon shot noise and sky background, dark current and read noise according to the Gauss–Markov theorem (Luenberger 1969). We typically choose w = 1 for simplicity because the noise in I(x, y, t) usually is not a limiting factor. For most of the objects, the motion is linear and can be modeled as

$$\begin{eqnarray}X(t) & = & {x}_{c}+{v}_{x}(t-(N+1)/2)+{\epsilon }_{x}(t),\\ Y(t) & = & {y}_{c}+{v}_{y}(t-(N+1)/2)+{\epsilon }_{y}(t)\end{eqnarray} \tag{ 2 }$$

where (x_c, y_c) is the location of the object at the center of the integration time interval and (v_x, v_y) is the velocity of the linear motion. (_x(t), _y(t)) is the tracking error with respect to sidereal, which can be estimated, for example, by an average of the centroids of a few bright reference stars in each frame. The location (x_c, y_c) and velocity (v_x, v_y) are solved simultaneously using a least-squares fitting together with parameters α and I₀, which give photometry and background intensity. We note that this approach to integrate a data cube avoids streaked images and keeps the jitter effect from atmosphere and telescope pointing common between target and reference objects, thus achieving accuracy comparable with stellar astrometry, which we shall present in Section 4.1.

After the star identification using a triangular matching scheme, we have a mapping between pixel coordinates of the N_s stars (X, Y)_i, i = 1, 2, ⋯, N_s and the catalog positions, the R.A. and decl., (R.A., decl.)_i, i = 1, 2, ⋯, N_s. We first convert the sky positions (R.A., decl.)_i into positions (ξ, η)_i in the tangent plane and then solve for the following fitting

$$\begin{eqnarray}&&a{\xi }_{i}+b{\eta }_{i}+{X}_{0}={X}_{i}-\displaystyle \sum _{n,m,2\leqslant n+m\leqslant {N}_{d}}{C}_{n,m}^{X}{X}_{i}^{n}{Y}_{i}^{m}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&c{\xi }_{i}+d{\eta }_{i}+{Y}_{0}={Y}_{i}-\displaystyle \sum _{n,m,2\leqslant n+m\leqslant {N}_{d}}{C}_{n,m}^{Y}{X}_{i}^{n}{Y}_{i}^{m}\end{eqnarray} \tag{ 4 }$$

where N_d is the order of polynomial, for which we found it sufficient to have N_d = 5 for mas level calibration, and a, b, c, d defines the linear transform between the sky plane position and measurements in pixel coordinates, and ${C}_{n,m}^{X},{C}_{n,m}^{Y}$ are distortion model coefficients. (X₀, Y₀) is an offset between the origins of tangent plane and the origin of pixel coordinate, which would be zero if the telescope had zero pointing error. To improve the fitting accuracy, we can weight the fitting according to the inverse of the variances of uncertainties of X_i and Y_i from the centroiding fitting, which is effective when faint reference stars are included.

ST avoids streaked images by having exposures short enough so that the moving object does not streak in individual images,⁷ allowing us to achieve astrometry similar to stellar astrometry for NEAs.

To illustrate this, we estimate the centroid positions of a bright (apparent magnitude of ∼13.2) asteroid 1983 TB, observed on 2017 December 20, at distance of ∼0.09 au from the Earth, with respect to reference stars in each of the 1 Hz frames. Figure 3 shows the frame-to-frame standard deviations of the differential centroids between stars (blue plus sign for R.A. and red dot for decl. respectively) as function of the angular distances. The frame-to-frame standard deviations serve as a measure of random errors for integration of 1 s, the exposure time. Similarly, we plot also the frame-to-frame standard deviations of the differential centroid of asteroid 1983 TB with respect to reference stars in the same field after removing a linear motion of (0.596, −0.599) pixel/second consistent with JPL Horizon ephemeris for comparison, marked with the blue squares and red diamonds representing R.A. and decl., respectively. We can see that the random errors are comparable for the differential centroids between asteroid and stars and stellar differential centroids. The random errors are mainly due to atmospheric turbulence because both the asteroid and reference stars are brighter than 16th magnitude, thus increasing with the angular separation.

Zoom In Zoom Out Reset image size

In general, we estimate the centroid using the "moving PSF fitting" to the data cube as described in the Section 3.2. The "moving PSF fitting" to the data cube achieves essentially the same accuracy as averaging the estimated centroids from each short-exposure frames. Figure 4 displays the Allan deviations for averaging centroid estimated from each individual frames based on totally 300 frames after removing the linear motion of the asteroid (blue lines with asterisk and triangle markers for R.A. and decl., respectively) as function of integration time. Similarly, we can divide the 300 frames into a bunch of "sub data cubes" corresponding to a specific integration time. For example, we have 30 "sub data cubes" for 10 s integration with each "sub data cube" containing 10 frames. We apply "moving PSF fitting" to each of the 30 "sub data cube" to obtain 30 centroid estimates. We can then estimate the centroiding error as the standard deviation of the 30 centroid estimates after removing the linear motion of the asteroid. The red squares and diamonds in Figure 4 represent R.A. and decl. centroiding errors using "moving PSF fitting" for the corresponding integration time. The integration shows the expected the inverse of the square root of integration time behavior (gray line) because the errors are uncorrelated noises. The "moving PSF fitting" gives a similar performance to estimates from averaging centroids estimated using individual frames. Because we're using individual short-exposure frames, the centroiding of asteroids and stars has similar performance as shown in Figure 3, and the integration down follows the inverse of square root of integration time behavior; thus, ST yields NEA astrometry with accuracy similar to that of stellar astrometry. The "moving PSF fitting" to the whole data cube is useful when the asteroid is dim because centroiding individual short-exposure frames may become too noisy to converge reliably.

Zoom In Zoom Out Reset image size

The situation is, however, very different when we estimate centroid using streaked images as illustrated in Figure 5. The left plot shows detection of a well-tracked PSF (unstreaked) of width W, in the sense of full-width-half-maximum, and a streaked image with streak length L. The light signal detected by the streaked image is spread over an area of a factor 1 + aL/W larger than the unstreaked image, where a is of order 1 depending on the detailed shape of the PSF. Supposing that we integrate all of the pixel intensities contained in the image, the total photon signal and the corresponding shot noise from the object that is imaged would be the same for both streaked and unstreaked images. But the variance of the total noise is larger for the streaked image because of the extra sky background and detector noises (by a factor ∼1 + aL/W in variance) from the larger detecting area, leading to a degraded S/N, which can be as large as $\sim \sqrt{1+{aL}/W}$ when limited by background noise. When estimating centroid using streaked images, the degradation of astrometry across the streak is caused by this degradation of S/N. Along the streak, the degradation is far greater because the centroid position of a streak along the streak is mainly determined by the signals at both ends and is insensitive to the central portion of the streak. Therefore, in the whole streak, only the two ends, or approximately an area of an unstreaked PSF, contains a signal for determining centroid along the streak, leading to a degradation in the total signal by a factor of ∼1 + aL/W. Along the streak, another major error of centroiding is caused by the fact that jitter effect from the atmosphere and telescope pointing becomes non-common between streaked and unstreaked images, which is illustrated by the right plot in Figure 5. Because the centroid of the streaked image is only sensitive to the signals at both ends, only the jitter effects at the beginning and end of the integration affect the centroiding along the streak, while for a well-tracked reference star, the jitter effects during the whole integration affect its centroiding. This non-common effect contributes an extra error with variance ∝L/W to the centroiding along the streak.

Zoom In Zoom Out Reset image size

We now show the discussed degradation of astrometric precision from using the streaked images, simulated by combining our short-exposure images, compared with the precision of ST. The differential centroid residuals of asteroid 1983 TB with respect to a nearby reference star are calculated by removing its linear motion. The performance of astrometry from centroiding-streaked images clearly shows degradation with respect to ST, especially along the streak as displayed in Figure 6, where we plot residuals from using ST (red dots), fitting with a two-dimensional elliptical Gaussian PSF (2d-G, black cross), and fitting with a streaked-Gaussian PSF (s-G, blue squares) together with the corresponding 1983 TB images of different streak lengths. At timescale of 3 s, the centroiding error mainly comes from atmospheric jittering effect. As we increase the integration time from 3 to 6 s, the spread of the residuals shrinks as we integrate down jitter errors. Further increasing integration time, the performance along the streak degrades due to the fact that the jitter effects from atmosphere and telescope pointing are no longer common between the asteroid and the nearby reference star, as we just discussed.

Zoom In Zoom Out Reset image size

Figure 7 shows centroiding performance along and across the streak using a two-dimensional (2D) elliptical Gaussian PSF, streaked-Gaussian PSF Vereš et al. (2012), compared with ST method for asteroids 1983 TB (left) and 2003 EB50 (right) respectively. Asteroid 1983 TB moves at ∼038 per second, so it takes about 5 s for the streak length L to be comparable with size of the PSF, W ∼ 2'', at which the precision of centroiding using streaked images starts to degrade significantly. Asteroid 2003 EB50 moves at ∼ 01 per second, so it takes more than 10 s to see the degradation of accuracy. At 1 Hz, errors are dominated by the air turbulences because the asteroids and reference stars are all bright. The reference star for 2003 EB50 is ∼25'' away, and the reference star for 1983 TB is ∼68'' from the asteroid. A closer reference star is the reason why the error at integration time of 1 s is smaller for 2003 EB50. As expected, for long steaks, streaked-Gaussian PSF fitting gives a better performance than a 2D elliptical Gaussian PSF fitting (Vereš et al. 2012) for centroiding along the streak, because the streaked-Gaussian PSF is a higher fidelity model of the intensity distribution of a streak. Across the streak, the atmospheric effect is common between the asteroid and reference star; therefore, the performances are comparable until the degraded S/N from the increased sky background starts to degrade the performance. For asteroid 1983 TB, the sky background is ∼5 photon per second per pixel. As we integrated the air turbulence down to less than 15 mas at about 30 s, the sky background is 150 photon, becoming comparable with the average pixel intensity of the asteroid signal, which is a few hundreds of photons. Therefore, sky background noise starts to degrade the S/N, thus increasing the the centroiding error. Similarly, this is also true for asteroid 2003 EB50. The degradation of S/N due to extra sky background noise happens at about 50 s integration because it moves slower.

Zoom In Zoom Out Reset image size

We now present asteroid astrometry from observations taken since 2017 June. Gaia DR1 (Gaia Collaboration et al. 2016) is used for our data reduction, but it does not have proper motions except for the Tycho stars. We estimate the proper motion for the stars that are in both Gaia DR 1 and UCAC4 (http://ad.usno.navy.mil/ucac/readme_u4v5) catalogs, assuming a linear motion between epoch J2015 (Gaia) and J2000 (UCAC4). Because the accuracy of UCAC4 is about ∼50 mas, so our proper motion is only accurate to a few mas per year. Propagating from Gaia's 2015 astrometry to about 2.5 years into 2017, we have a propagation error generally less than 10 mas.

Figure 8 shows astrometry of asteroid 2005 UP156 observed on 2017 June 14 with the left plot displays respectively the R.A. and decl. from our instrument (red squares) on top of the ephemeris from the JPL Horizon System (https://ssd.jpl.nasa.gov/?horizons) (blue curve). The right plot displays the residuals, the difference between our astrometry and JPL Horizon system solution #274. Each of the data point is derived from integrating 100 × 1 Hz frames (100 s). The spread of our measurements have standard deviations of ∼10 mas. There is clear bias around 30 mas, which we believe is due to the uncertainty of the JPL Horizon System solution ∼100 mas at the 3σ level. 2005 UP156 is a kilometer size asteroid and was at ∼0.16 au from Earth during the observation with an apparent magnitude of ∼14.7, moving at a slow speed of ∼2'' per minute.

Zoom In Zoom Out Reset image size

Figure 9 shows similar plots for a faster moving asteroid, 2010 NY65 observed on 2017 June 29. This asteroid has a size ∼200 m and was at ∼0.04 au away from the Earth during the observation with apparent magnitude ∼16.8. It moves at speed of 02 per second.

Zoom In Zoom Out Reset image size

Again, the standard deviations (spread of our measurements) are ∼14 mas, and there are biases consistent with the JPL Horizon System solution uncertainty of solution #83. Each of the measurements was obtained from integrating 100 × 1 Hz frames. It is possible to further average down random errors to have better precision; however, we are limited by systematic errors. Figure 10 displays astrometry residuals of 2005 UP156 from observations over five weeks, where each data point represents the average value of all the observations within one night with the standard deviations of data within a night shown as the error bar. The spread of daily mean astrometry however, is still ∼10 mas, not smaller from averaging data over one night; this suggests systematic errors. The consistency of the overall constant bias over five weeks is unlikely due to our systematic errors; it is likely to be the prediction uncertainty of the JPL Horizon system, whose 3σ is ∼ 100 mas.

Zoom In Zoom Out Reset image size

It is usually sufficient to use only third-order polynomials to correct the field distortion to 10 mas over the 19' × 16' field. Our major error comes from chromatic distortion effect. A typical field distortion calibration has 40–50 mas root-mean-square (rms) of errors with 50–200 reference stars in the field. Near the center of the field over a region about 6' × 6', we can correct to 10–20 mas. In case of a wide range of stellar spectra, our field distortion over the whole FOV can be as large as ∼100 mas. Currently, we mitigate this effect by restricting reference stars according to colors, e.g., using the difference of magnitudes in B and V bands. Applying an R or V bandpass filter (∼ 120 nm) helps us reduce field distortion calibration error to less than 20 mas over the whole 19' × 16' field, and better than 10 mas over the center 6' × 6' region.

Since 2017 June, we have observed a dozen asteroids listed in Table 1 and reported more than 400 data points to the Minor Planet Center (https://minorplanetcenter.net/). Each of our data points represent an integration over 300 s of data. The rms of our astrometry residuals is about 50 mas as shown in Table 2, where we have obtained the rms of NEA astrometry residuals for Pan-STARRS, which has the smallest rms among the major NEA surveying facilities, from the reference Vereš et al. (2017). Because we believe a major contribution to the rms is the uncertainty of the predicted asteroid ephemeris from the JPL Horizon System, just like in the case of Figure 10, we also compute the standard deviations of our data points over each night, exactly like we did for the spread in Figure 9, which we believe is closer to our accuracy for most of the data points.

The left plot in Figure 11 displays the histogram of the standard deviations of our astrometry residuals of reported data within a night, and the right plot shows the corresponding cumulative distribution. We do not know why the decl. has slightly better performance. More than 85% of the data have accuracy better than 50 mas. The largest residuals can come from centroiding errors due to confusion field or random errors for faint targets. For bright targets (mag < 17), our accuracy are typically better than 20 mas with the best ∼10 mas.

Zoom In Zoom Out Reset image size

Most of the Pan-STARRS astrometry statistics comes from slowly moving asteroids, probably main belt asteroids, with an integration time of 45 s. According to Figure 2 in Vereš et al. (2017), as the speed of asteroid motion increases, the degradation becomes clear. In contrast, our accuracy is independent of how fast the target moves in sky with the ST technique.