DETECTION OF A FAINT FAST-MOVING NEAR-EARTH ASTEROID USING THE SYNTHETIC TRACKING TECHNIQUE

As introduced by Shao et al. (2014), synthetic tracking is a technique that integrates a set of short-exposure frames to simulate the tracking of the telescope at a specific velocity, called the synthetic tracking velocity. In synthetic tracking, we search for an object in four-dimensional space (x, y, v_x, v_y) in a three-dimensional data cube. Figure 1 schematically shows the integration of the displaced frames according to the two-dimensional (2D) velocity needed to track an NEA. Note that the locations of the NEA on the CCD are aligned as a consequence of shifting the frames. For the purposes of a blind search for faint NEAs, one has to choose all feasible tracking velocities. The objective is to search in the synthetic tracking velocity space for every asteroid whose image after integration does not streak (no trailing loss). This non-streaked image from synthetic tracking is very close to the image obtained with a telescope that is actually tracking the asteroid. Because the integration is done in post-processing, synthetic tracking can generate non-streaked images for both the moving asteroid and the background stars.

Zoom In Zoom Out Reset image size

Instrumentally, synthetic tracking is enabled by the new generation cameras that provide both fast frame rate and low read noise. Our observation on the Palomar 200 inch telescope used an Andor EMCCD camera ⁴ operating with an EM gain of 200, making the read noise benign compared with the sky background noise even for a 16.7 Hz frame rate. The key data processing for the synthetic tracking technique is a search for signals (or bright spots) in the synthetically integrated images within a grid of tracking velocities. As mentioned in (Shao et al. 2014), this process is computationally intensive. To overcome this difficulty, we use GPUs to perform synthetic tracking at different velocities in parallel, thus enabling a nearly real-time data processing.

Synthetic tracking significantly improves the detection S/N and the accuracy of astrometry by avoiding streaked images for detecting fast-moving NEAs. An additional benefit of synthetic tracking is that it allows one to estimate the velocity of an NEA using only 30 s of data, which provides timely follow-up observations with an approximate target trajectory.

Depending on the precision at which the images are shifted, synthetic tracking may be implemented in two different ways, which we call the Integer Pixel Shift and Add (IPSA) and the Continuous Shift and Add (CSA) methods, respectively. Because of the computational expediency, we use IPSA synthetic tracking in the GPU-enabled search for faint NEAs and apply CSA in post-analysis. These methods are discussed below.

2.2.1. Synthetic Tracking with IPSA Method

The IPSA synthetic tracking integrates frame images after shifting each by an integer amount of pixels (i.e., no shift of image is done at a fraction of a pixel) according to the tracking velocity and then adding the signals collected on each of these frames. The quantity of interest here is the integrated intensity of the synthetically tracked image, I_IPSA(x, y, v_x, v_y), which is computed as

$$\begin{eqnarray} I_{\rm IPSA}(x,y,v_x,v_y) &=& \sum _{n = 0}^{{{N_{\rm f}}{-}1}} I_n (x + {\tt round}(n v_x), \; y \nonumber\\ && + {\tt round}(n\,v_y)), \end{eqnarray} \tag{ 1 }$$

where n = 0, 1, ⋅⋅⋅, N_f − 1, labels the frames. The coordinates (x, y) represent the location of a particular pixel on a CCD in (row, column) order, with vector (v_x, v_y) being the 2D synthetic tracking velocity. I_n(x, y) is the intensity of the pixel at (x, y) in the n-th frame. The ${\tt round}$ function rounds its argument to the nearest integer. Because the frames have zero bias (i.e., the background has been subtracted), in the case where round values are out of the bound of the frame index, zeros are filled in. To properly estimate the noise level, it is necessary to record the actual number of frames that contribute to a signal when its location is near frame boundary. The frames that contribute "filled zeros" to the signal should be excluded because the signal is outside of their appropriate boundaries. The advantage of the IPSA synthetic tracking is that it requires minimal number of arithmetic operations. It is suitable for an extensive search for NEAs, where many such operations are performed. We implemented the IPSA synthetic tracking on GPUs to search for NEA signals over different tracking velocities in parallel.

2.2.2. Synthetic Tracking with CSA Method

The CSA synthetic tracking method integrates short-exposure frames by shifting the frames according to the displacement determined by the velocity of tracking, but not limited to integer number of pixels. The images are displaced using the spectral interpolation or the Fourier space interpolation method (Boyd 2001; Zhai et al. 2011). The spectral interpolation method for shifting the images is based on the fact that the Fourier transforms of the original and shifted images are related by a linear phase factor in spatial frequencies.

Mathematically, if we shift an image I(x, y) by (Δx, Δy) in row and column, the displaced image is related to the original image via

$$\begin{equation} I^{\rm FT} (x + \Delta x, y + \Delta y) = {\rm FT}^{-1} \lbrace {\rm FT} \lbrace I(x, y) \rbrace e^{-2 \pi i (k_x \Delta x + k_y \Delta y)}\rbrace, \end{equation} \tag{ 2 }$$

where FT represents the Fourier transform

$$\begin{equation} {\rm FT}\lbrace I(x,y)\rbrace \equiv \sum _{x, y} I(x,y) e^{2\pi i (k_x x + k_y y)}, \end{equation} \tag{ 3 }$$

and FT⁻¹ is the inverse Fourier transform. For a tracking velocity vector (v_x, v_y), the CSA synthetic tracking image is then computed as the following expression:

$$\begin{equation} I_{\rm CSA}(x,y,v_x,v_y) = \sum _{n=0}^{{{N_{\rm f}}{-}1}} I_n^{\rm FT} (x + nv_x, y + n v_y). \end{equation} \tag{ 4 }$$

To illustrate the synthetic tracking technique, we apply CSA to observational data of asteroid 2009BL₂ collected on the Palomar 200 inch telescope on 2013 April 3 with a Sloan g' filter in front of the CCD. The short-exposure frames were taken at a frame rate of 2 Hz. Figure 2 shows two synthetic tracking images for tracking the background stars (left) and the asteroid 2009BL₂ (right), respectively, while using CSA synthetic tracking with a total of 960 frames of data.

Zoom In Zoom Out Reset image size

The bright background star has a visual magnitude of 16 and the asteroid has an apparent magnitude of 18.5. In the left image, the stars are tracked. The asteroid appears as a streak with surface brightness significantly less than that of the star. In the right image of Figure 2, the asteroid is tracked and the stars are streaked. The surface brightness of the 16th magnitude star is now comparable with that of 2009BL₂, which is of 18.5th magnitude. As illustrated, the trailing loss is ≈2.3 visual magnitudes (a factor of 8.5). Due to the trailing loss, some of the faint background stars in the left image can barely be seen in the right image of Figure 2. By synthetically tracking the stars and the asteroids, they have essentially the same point-spread function (PSF). To quantify the PSF, denoted as P(x, y), we adopt Moffat's PSF template (Moffat 1969)

$$\begin{equation} P(x, y) = \left(1 + \frac{x^2 + y^2}{R^2} \right)^{-\beta } \end{equation} \tag{ 5 }$$

to model the measured seeing-limited PSF. Here, the quantity R is a size scale and β specifies the rate of how fast the PSF falls off while moving away from the center. The FWHM of Moffat's PSF is $W = 2R\sqrt{2^{1/\beta } - 1}.$ Figure 3 displays the radial intensities of the synthetic tracking PSFs computed for the 16th magnitude star (star marker) and the asteroid 2009BL₂ (circle marker), respectively. The corresponding fitting curves using Moffat's PSF are also displayed. The synthetic tracking PSFs of the star and asteroid are very similar, with W of the asteroid 2009BL₂ being only 4% larger than that of the star.

Zoom In Zoom Out Reset image size

Because the error due to neglecting sub-pixel (pix) displacement is much smaller than the size of the PSF (∼3 pixels), the IPSA synthetic tracking achieves synthetic tracking images very close to that of the CSA synthetic tracking. The third curve (diamond marker) in Figure 3 shows the radial dependency of the intensity profile as a function of the distance from the center of the image of asteroid 2009BL₂ using IPSA. The FWHM is only 3% larger for the IPSA synthetic tracking (W = 3.21 pixels) compared with the CSA synthetic tracking image (W = 3.11 pixels). The corresponding S/N is degraded by less than 3%. Therefore, IPSA synthetic tracking is used in the GPU search; hereafter, reference to "synthetic tracking" implies the IPSA method unless otherwise stated. CSA is only used in post-analysis for generating refined PSFs of tracked objects.

To detect faint objects, we need to accurately estimate the background bias to avoid excessive false positives (underestimated background) or missed detections (overestimated background). In the data frames, a majority of the data points measures the sky background; only a small portion detects the light from the stars and asteroids. Therefore, we can estimate the background by removing the signals as outliers. This is done by starting with including all the data points present in the background data set and then iteratively removing the signals that are above the average of the background by a threshold, e.g., 5σ, with σ being is the standard deviation.

The steps characterizing the background estimation procedure are given as follows.

• 1.  
  Let y_i, with i ∈ {1, 2, ⋅⋅⋅, N}, represent all the data points available including both background and signal measurements. Also, let B denote the set of the data points representing the background to be determined. We initialize B to be the entire set of data points and estimate B iteratively.
• 2.  
  Compute the sample mean 〈y〉 and standard deviation of the data in set B, σ_y, as

  $$\begin{eqnarray} &&\langle y \rangle = {1 \over N(B)} \sum _{i \in B} y_i , \sigma _y = \left[{1 \over N(B)-1} \sum _{i \in B}\left(y - \langle y \rangle \right)^2 \right]^\frac{1}{2}, \nonumber\\ \end{eqnarray} \tag{ 6 }$$

  with N(B) being the number of elements in set B.
• 3.  
  Update set B according to the rule B = {i|y_i ⩽ 〈y〉 + ξσ_y}, where ξ is the threshold chosen for the estimation.
• 4.  
  Iterate this procedure by going back to Step 2) above until the process converges.

If the population of the background is large compared with that of the signal, this process converges fast. This is because the background statistics, the mean and standard deviation, can be easily established and are relatively stable. In general, ξ should be chosen according to the population of the sample. The larger ξ is, the faster the convergence can be reached. However, a larger threshold means a higher chance of including weak signals into the background. At the beginning we used ξ = 5 to estimate an overall pixel-independent background using 50 frames of data and then used ξ = 4 to estimate a pixel-dependent background with ∼2000 frames.

Co-moving PSF fitting optimizes the astrometry and velocity of the detected object by a least-squares fitting of multiple short-exposure frames to a PSF that moves at the velocity of the object. This assumes a priori knowledge of the PSF function. It is performed after a synthetic velocity search where the location and tracking velocity have been determined to certain accuracy depending on the velocity grid of the search. The estimations of the location and velocity from synthetic velocity search are used as initial conditions for the optimization routines. The mathematical formulation for the fitting is based on explaining the observed multiple short-exposure signals as a moving PSF. This can be expressed as a minimization of the following least-squares cost function:

$$\begin{eqnarray} && C(v_x, v_y, x_c, y_c, \alpha, I_0) \equiv \sum _{n=0}^{{{N_{\rm f}}{-}1}}\sum _{x,y} \left |\right. I_n(x,y) -\alpha P \left(x-x_c\right. \nonumber\\ &&\quad \left.-v_x \left(n-{{{N_{\rm f}}{-}1} \over 2}\right), y-y_c-v_y \left(n-{{{N_{\rm f}}{-}1} \over 2} \right) \right) {-}I_0 \left|\right. ^2, \nonumber\\ \end{eqnarray} \tag{ 7 }$$

where P(x, y) is the PSF function, (x_c, y_c) is the location of the object at the mid-epoch of all the frames, and (v_x, v_y) is the velocity of the moving object relative to camera frame (i.e., the required tracking velocity). Quantities α and I₀ are two extra fitting parameters specifying linear and constant levels with respect to the PSF.

We use a bright star in the field to determine PSF model parameters. We estimate the PSF P(x, y) by fitting the Moffat's PSF template to a nearby bright star in the field. Because we critically sampled the PSFs (Zhai et al. 2011), in principle we could reconstruct the PSF from the star image itself, which is especially important for micro-arcsecond astrometry. For mas-level astrometry, it is sufficient to use the model for the PSF given by Equation (5), which helps reduce the amount of numerical computations during centroid fitting.

Minimizing the cost function C(v_x, v_y, x_c, y_c, α, I₀) gives an estimate of the velocity of the asteroid (v_x, v_y) and the location of the object (x_c, y_c) at the mid-epoch of all the frames,

$$\begin{eqnarray} &&(\hat{v}_x, \hat{v}_y, \hat{x}_c, \hat{y}_c, \hat{\alpha }) = {\arg \min }_{v_x, v_y, x_c, y_c, \alpha, I_0} C (v_x. v_y, x_c, y_c, \alpha, I_0), \nonumber\\ \end{eqnarray} \tag{ 8 }$$

where the estimate $\hat{\alpha }$ measures the signal level and is used to compute the S/N. We adopted the Matlab lsqnonlin routine to perform this optimization. We note that the co-moving PSF fitting can be applied when the objects are so faint that neither the star nor the asteroid are detectable in a single frame. The proper convergence relies crucially on a priori knowledge of the location and velocity of the detected faint object from the synthetic tracking search. Using the estimated velocity $(\hat{v}_x, \hat{v}_y)$, we use the CSA synthetic tracking to obtain a refined PSF of the tracked object.

Because the PSF is larger than a CCD pixel, the matched filter technique (Turin 1960) is used to improve the S/N. The matched filter convolves the data with a low pass filter, whose impulse response profile matches the PSF and yields the optimal S/N. With estimate $\hat{\alpha }$ in hand, we express the detection S/N after applying the matched filter as

$$\begin{equation} {\rm S/N } = {\hat{\alpha }\sqrt{{N_{\rm f}}}\sqrt{\sum _{x,y} P(x,y)^2} \over \sigma _n} , \end{equation} \tag{ 9 }$$

where the sum is taken over the pixels and σ_n is the standard deviation of the noise detected by each pixel (assumed to be uniform across pixels). Factor $\sqrt{{N_{\rm f}}}$ is from integrating N_f independent frames. A convenient approximation that expresses the S/N in terms of N_ph and the FWHM of the PSF is

$$\begin{equation} {\rm S/N } = 0.6 {N_{\rm ph}\sqrt{{N_{\rm f}}} \over \sigma _n W} \,. \end{equation} \tag{ 10 }$$

We studied the sensitivity of astrometric solutions to the noise present in the system using simulations. To do this, we adopt Moffat's PSF with R ≈ 3.3, β ≈ 4, which was obtained from fitting the 16th magnitude star images in the field taken while observing 2009BL₂. The co-moving PSF fitting procedure is applied to simulated signals to estimate the velocity and astrometric position of the moving object. The simulation was performed for many different signal levels and three different noise levels. Figure 4 displays the astrometric error RMS as a function of the total number of photons collected for three different levels of background noise, parameterized by σ_n, the standard deviation of the background noise per pixel. The squares, triangles, and circles represent the three different noise levels, respectively. By inverting the Hessian matrix (Press et al. 1986) of the least-squares cost function (7) for Moffat's PSF and assuming Poisson statistics for photon detection, we derived the following empirical formula to assess the uncertainty for astrometric position:

$$\begin{equation} \sigma _{x,y} = \sqrt{{0.457 \over N_{\rm ph}} + {1.37 \sigma _n^2 W^2 \over N_{\rm ph}^2}} {W \over \sqrt{{N_{\rm f}}}}, \end{equation} \tag{ 11 }$$

where W is the FWHM of the PSF in pixels and N_ph is the total number of the photons detected per frame. Figure 4 shows that the estimated RMS values using simulation agree well with the empirical formula. For faint objects, the background noise, given by the second term under the square root in Equation (11), dominates. In this case, it is convenient to approximate Equation (11) as

$$\begin{equation} \sigma _{x,y} \approx {0.65 W \over {\rm S/N }} , \end{equation} \tag{ 12 }$$

where S/N for detection is given by Equation (9).

Zoom In Zoom Out Reset image size

Precision of the synthetic tracking velocity is related to the centroid precision via

$$\begin{equation} \sigma _{v_x,v_y} = {\sqrt{12} \, \sigma _{x,y} \over T_{\rm tot}}, \end{equation} \tag{ 13 }$$

where T_tot is the total time duration covered by all the short-exposure frames.

We conducted a 12 hr blind search using the Caltech HIgh-speed Multi-color camERA (CHIMERA) ⁵ at the Palomar 200 inch telescope, with the Wynne correctors, over two nights on 2013 September 11–12. CHIMERA uses two Andor iXon3 888 EMCCDs with (1024 × 1024) pixel detectors, allowing readout at a 10 MHz data rate in two colors with an effective read noise of ≪1e⁻ with EM gain applied. We used this instrument to take images at 16.7 Hz with EM gain of 200 to avoid excessive read noise. This allowed us to scan over the sky at 5 ''/s rate continuously instead of slewing and stopping the telescope repeatedly. Each scan was along the R.A. direction and lasted for approximately 1 hr to cover a field of size (53 × 0046). The scan rate was chosen so that each object would remain in the FOV = 2.7 arcmin (') (or 0046) for about 30 s. We binned (2×2) CCD pixels to have (512×512) frame pixels. No additional optical filters were applied. During each night, we divided the allocated 6 hr observation into three two hour pairs, with each pair scanning over the same field twice, to cover three different fields. To facilitate data processing, we took calibration data sets including the shortest exposure (10 μs) frames allowed for estimating frame bias and twilight frames for estimating the flat field responses. At a future time, when we will be able to process the data in real time, the second one hour observation would only be conducted for follow-up observations of a detected object.

The entire data process contains four main steps: preprocessing, detecting and removing objects from frames, synthetic tracking velocity search, and post analysis to refine the results. Figure 5 displays a flow chart showing the relationships between the steps.

Zoom In Zoom Out Reset image size

In pre-processing, we adjust data frames using calibration data of the A2D bias and flat-field responses. A pixel-dependent sky background is estimated and subtracted to make the data frames zero-biased. We remove the bright pixel data caused by cosmic ray events by setting the values to zero. In the second step, we detect stationary objects (stars and galaxies) as well as bright asteroids and remove their signals from each frame by setting the relevant pixel values to zero. The data is then passed to the synthetic tracking velocity search, where the bright signals above the detection threshold in the synthetically tracked image at each of the grid velocities are detected. The results are then reported for the post-analysis to refine the detection S/N and astrometry using optimization.

In pre-processing, we first deal with known features in the measurements. To enable asteroid detection, the A2D bias map is subtracted from each of the data images. All the frames are then divided by the flat-field responses, estimated using twilight frames, to compensate for pixel-dependent throughput. Next, we remove the cosmic ray events, which generate high counts that only stay in a single frame. Cosmic ray events are detected by first zeroing out the bright stars that have a signal level comparable with the cosmic ray events in a single frame. (We differentiate the bright stars from the cosmic ray events by integrating over hundreds of frames so that they are much brighter than the extra counts due to cosmic ray events, which only appear in one frame.) After zeroing the signals from the bright stars, the differences between consecutive frames are used to identify highly varying signals on a frame-to-frame basis and we attribute them to cosmic ray events.

We then estimate a pixel-dependent background intensity and an average background noise level using about 2000 frames of data. To do this, we first estimate a uniform background (pixel-independent) and noise level by computing the sample mean and standard deviation of background data in a small number (e.g., 50) of frames. It is convenient to divide all the data points into two parts, background data and signals, assuming signals are above the background for a specified threshold. The background data set is estimated by iteratively removing signals above the sample mean by 5σ. (See Section 3.1 for the details of the algorithm.) Using this initial estimation of the background and noise levels, we are able to eliminate most of the strong signals from each frame.

For the rest of the data, we apply the same algorithm to the data measured by each pixel (∼2000 frames) and estimate a pixel-dependent background using a threshold of 4σ. After removing signals from each of the frame, we further detect and remove background stars that have S/N above 4 after 30 s of integration (the average observation time for each object in the field for the chosen scan rate) by performing a synthetic tracking at the sidereal rate. (Note that S/N = 4 is only used for estimating background and has nothing to do with the detection S/N threshold 7). After removing bright objects and faint stars, the field-dependent background is estimated by taking the average over the background data (typically more than 1000 data points after removing signals) for each pixel. We also compute the frame-to-frame variation of the background for each pixel and then take the average over all the pixels as an estimate for the noise level of background. The estimated sky background is subtracted from each frame to have zero-biased frames.

Before searching for faint asteroids, we detect and remove the bright objects from each frame. With the estimated background noise level, it is straightforward to detect all the bright objects in each of the zero-biased data frames with a single frame S/N above our threshold of 7. We then synthetically track at the sidereal rate to detect faint static objects with S/N above 7 after 30 s of integration. The detected objects are passed to post-analysis for further identification. For example, a very bright asteroid may be detected if its single frame signal is above the specified detection threshold. Also, if an asteroid does not move much during the observation period, it may be detected as a static object.⁶ We remove the signals of the detected objects by setting the values at the relevant pixels to zero. Before passing to the next step, we apply the matched filter by convolving each data frame with the PSF profile to minimize computational operations during the extensive synthetic tracking velocity search.

The synthetic tracking velocity search examines the data over a two-dimensional tracking velocity grid. To avoid trailing loss, the grid spacing should be no more than the size of the PSF divided by the integration time, i.e., the speed at which the motion is less than the size of the PSF over the integration time (30 s in our case). An IPSA synthetic tracking is performed for each grid velocity. Because the frames have zero background, we fill in zeros for the missing data at the boundaries of frames that are displaced. If the synthetic tracking image shows a signal above the noise level by the detection threshold of 7, we report this signal level together with four numbers (x, y, v_x, v_y), where (x, y) specifies the location of the detected signal and (v_x, v_y) is the grid velocity at which the synthetic tracking yields the signal. The post-analysis uses this information in optimization schemes to refine the detection S/N and compute the astrometry. Because data processing for the synthetic tracking search is independent between different tracking velocities, it can be easily speeded up by implementing parallel computing.

We have implemented the search using NVIDIA's Tesla K20c GPU⁷ to accelerate the faint object search. The K20c is based on the Kepler architecture equipped with 2496 CUDA processing cores and 5 GB of GDDR5 RAM. Peak single-precision processing performance is 3.52 TFLOPS (10¹² floating operations per second). Our performance is currently limited by the memory bandwidth, which is 209 GB/s. The search software is implemented in C/C++. To process a 90 s data cube of 1,500 frames of (512×512) images searching over a (100 × 100) synthetic tracking velocity grid covering a velocity range of ±12 °/day in both R.A. and decl., the average GPU processing time is under 90 s. Our current grid spacing is 1 pixel/integration time, which is finer than needed with a PSF size of ∼3 pix. Using a coarser grid, this performance allows real-time processing. Our detection threshold is set to S/N = 7 to have less than 1% false positive probability per 30 s of data. See the next subsection for details.

After the synthetic-tracking-velocity search detects a signal (S/N ⩾ 7), we refine the detection using a least-squares fitting of a co-moving PSF to all the relevant data frames as described in the algorithm in Section 3. The least-squares fitting yields an optimal trajectory (average position and velocity) of the moving object relative to the camera frames. In the same fashion, we obtain solutions of the trajectories of stellar objects to compute the relative astrometry of the asteroid with respect to them. The fitting results also yield the signal strength above the background, which is described by αP(x, y), with α determined by the fitting procedure. The detection S/N can be then computed using Equation (9). This approach yields the same S/N as the matched filter scheme (Gural et al. 2005; Shucker & Stuart 2013) if a predetermined template filter velocity matches that of the asteroid. A powerful feature of our approach is that we search in the tracking velocity space for the faint object using parallel computing and then optimize the tracking velocity using a co-moving PSF fitting.

We estimate the false positive probability for initial detection as

$$\begin{equation} P_{\rm false-alarm}^{\rm detection} = N_x N_y N_{v_x} N_{v_y} \times {1 \over 2} {\tt erfc} \left({\rm S/N } \over \sqrt{2} \right), \end{equation} \tag{ 14 }$$

where ${\tt erfc}$ is the Gaussian complimentary error function, N_x, N_y are the dimensions of the CCD, and $N_{v_x}, N_{v_y}$ are dimensions of the synthetic-tracking-velocity search grid along R.A. and decl., respectively. This is derived by assuming the background fluctuates independently with Gaussian statistics so that the total false detection probability is the number of trials, which is the product of the total number of pixels and the total number of grid velocities $N_x N_y N_{v_x} N_{v_y}$, multiplying the false detection probability of one trial by $\rm erfc({\rm S/N}/\sqrt{2})/2$.

We choose S/N = 7 as the threshold for N_x = N_y = 512, $N_{v_x}=N_{v_y}=100$, to have a false alarm rate of less than 1% (∼0.34%) per 30 s data. For our 12 hr data, we expect about 0.34%× 12 × 3600/30 ∼ 5 false detections. In practice, artifacts from charge leakages or diffractions of bright objects yield many more false positives, which so far are rejected by human-involved inspections.

To further reduce the false-positive probability in post-analysis, we actually broke the data frames into a few segments (typically four as shown in Figure 7) to test the signal level of each segment separately. This ensures that the signal is not from a transient event that has not been properly handled in the data preparation steps for the synthetic-tracking-velocity search, and yields a much lower false-positive rate. Therefore, we had zero false positive for the two night search after the post-analysis.

Upon a detection, for confirmation, the false-positive rate becomes

$$\begin{equation} P_{\rm false-alarm}^{\rm confirm} = \Delta N_x \Delta N_y \Delta N_{v_x} \Delta N_{v_y} \times {1 \over 2} {\tt erfc} \left( \frac{\rm S/N }{\sqrt{2}} \right), \end{equation} \tag{ 15 }$$

where (ΔN_x, ΔN_y), and $(\Delta N_{v_x }, \Delta N_{v_y})$, respectively, represent uncertainties in the predicted location and velocity of the asteroid from the first observation. Now the velocity search space only needs to cover the uncertainty of the estimated synthetic velocity $(\Delta N_{v_x }, \Delta N_{v_y})$, which typically does not exceed one grid spacing of the tracking velocity search grid, which is much smaller than the detection search space of $N_{v_x} N_{v_y}$ grid points. Therefore, the false-positive probability of confirmation is significantly smaller than that of detection. The uncertainties in determining the location, on the other hand, grow linearly in time. Therefore, it is important to follow up with the second observation not too long from the initial detection so that the positional uncertainty will be smaller than the FOV. For an object detected at S/N = 7 with 1'' seeing, the estimated velocity uncertainty is approximately 025 day⁻¹. So, if we want the object to be within the FOV with 3σ confidence, the follow up should be scheduled no more than 1 hr for a 4'' FOV. This time window grows linearly with the detection S/N and the size of FOV.

In post-processing, we also go through the list of objects that are removed in the pre-processing (as described in Section 4.3) and estimate their velocities relative to the sidereal background to determine whether the motion is statistically significant using Equation (13). This will detect bright asteroids that have S/N ⩾ 7 with a single short-exposure frame as well as those asteroids whose streaked image still yields S/N above 7 (after the trailing loss).

We detected a faint asteroid using approximately 30 s data taken with the Palomar 200 inch telescope on 2013 September 11. The object emerged with a signal at S/N ∼ 15 in the co-added image after we performed an IPSA synthetic tracking at a velocity of approximately 632 day⁻¹ along the R.A. It was confirmed in a second data set acquired 77 minutes later at S/N ∼ 15 and with a similar velocity, which provides another 30 s of observation. We now provide the details of the detection, confirmation, astrometry, and photometry.

5.1.1. Initial Detection

The faint object was first observed at 2013-09-11T15:50:48 UTC for about 30 s. The synthetic-tracking-velocity search detected its signal at S/N ∼ 15 after performing an IPSA synthetic tracking on 524 frames, taken with an exposure time of 0.06 s, at a tracking velocity ∼638 day⁻¹ along the R.A.

Figure 6 shows the IPSA synthetic images tracking the asteroid and background stars, respectively, using 524 frames. The left image, tracking the asteroid, shows the asteroid at the middle near the bottom as pointed by the arrow with the surface brightness equivalent to a star of ≈23 mag. The detection S/N is 14.7, which, according to Equation (14), gives a statistical false positive probability of less than one part in 10³⁰. This asteroid moved ∼7'' during the 30 s observation. If a traditional 30 s exposure were used, the surface brightness of the streaked image would be approximately S/N ≈ 4 (apparent magnitude 25), falling below the detection threshold of S/N = 7. The right image in Figure 6 shows what a 30 s exposure(tracking sky background) would have yielded. The asteroid cannot be identified from the background fluctuations due to trailing loss.

Zoom In Zoom Out Reset image size

To further ensure that the signal did not come from a transient event, we broke the data into four segments, with each lasting about 7.8 s (the actual data of observations lasted 31.2 s) and found that the signal appeared in all four segments of data. The left four images in Figure 7 display the synthetic tracking images of the asteroid at the four different epochs corresponding to the four segments of data. Here we have arranged the images so that the background stars are at approximately the same locations to show the motion of the asteroid. The corresponding S/Ns for all the data segments are displayed.

Zoom In Zoom Out Reset image size

Using the co-moving PSF fitting, we refined the astrometric solutions by computing the trajectories of the asteroid and background stars relative to the camera frames. The relative position of objects can be estimated as the difference between their trajectories. We computed the average of the relative positions between the asteroid and background stars as the relative astrometry. We are able to identify background stars S1, S2, and S3 in the 2MASS catalog (Cutri et al. 2003) to be 2MASS J00500913+000035, 2MASS J00500710-0000242, and 2MASS J00500583-0000584, respectively. Using these, we were able to calibrate our plate scale to be 0324 pixel⁻¹ and that our camera frame is rotated clockwise relative to the Equatorial Coordinate System (ECS) by an angle ∼076.

The calibration enabled us to compute both relative and absolute astrometry. The right plot in Figure 7 shows the locations of the asteroid relative to the background star S3 (2MASS J00500583-0000584). The estimated uncertainty of astrometry due to noise is ≈60 mas for the 30 s observations using Equation (12) (the seeing is ≈13). The estimated velocity of the object is [637, −077 day⁻¹] with precision 016 day⁻¹ according to Equation (13) in the camera frame, giving [638, −068 day⁻¹] in (RA, decl.) after the 076 counterclockwise rotation. Because the asteroid is near the edge of the field, we estimated that the field distortion gives a correction to the velocity mainly along the decl. by 046⁻¹. With the correction, we estimate the velocity to be [638, −022 day⁻¹]. At the mid-epoch of the first observation, UTC 20130911T15:51:03, we found that the asteroid was at R.A. 0:50:4.39 decl. −0:0:55.5, (− 2166, 294) relative to the star S3 (2MASS J00500583-0000584). Field distortion is our main source of error, which we did not systematically calibrate. We estimate that the absolute astrometry is accurate to ∼05 even though the precision of asteroid centroid is ∼60 mas.

We perform photometry using the data from the four epochs shown in Figure 7. The mean and standard deviation of the photon fluxes were estimated to be (207 ± 18) photon/s. This gives an apparent magnitude of 23.1 ± 0.1 for the asteroid using a system throughput calibration from observing a known star of 9.1 mag, which yielded a flux of ∼8.1 × 10⁷ photon/s at the detector.

5.1.2. Confirmation

It is useful to estimate the expected number of detections for our two night search using the asteroid population distribution from Harris (Harris 2011). To do this, we divide the H-magnitude range from 28 through 31 into six bins with an equal bin size of 0.5 magnitude. The populations of the asteroids in all the bins from Harris are converted into densities of the asteroids (per AU³) using a method from P. Chodas (2013, private communication). For each H-magnitude bin, we determine the maximal distance for the asteroids to be observable according to the limiting magnitude of 24 for the Palomar 200 inch telescope using 30 s integration time. The total observation volume is estimated as the volume of a cone ∼ (max-distance)³ × (total field)/3. The total field covered over the two nights is 6° × 53  ×  0046 ≈ 1.5 deg² ≈ 4.5 × 10⁻⁴ rad². With this, we are ready to estimate the expected number of detections for each H-magnitude bin. For example, the density for the asteroid in bin at H = 29 is ∼1.8 × 10⁶/AU³ and the maximum distance allowed to observe an asteroid of H = 29 is ∼0.09 AU; we thus estimated 1.8 × 10⁶/AU³ × (0.09 AU)³ × 4.5 × 10⁻⁴/3 ∼ 0.2 NEA to be detected in this bin. Rounding to the nearest tenth, we estimated 0.2 asteroids per bin for all six bins, ⁸ giving a total of 1.2 asteroids for two nights. Even though this is a very crude estimation, it tells us that Harris's distribution is consistent with what we found.

Using the location and velocity of the NEA estimated based on the data of the first detection, we predicted that the asteroid will be at coordinates (R.A. = 0:51:36, decl. = −0:01:37) with an uncertainty of ±77 minutes × 016 day⁻¹ × 3 ∼ ±00255 at a confidence level of 3σ. We indeed detected a signal in the predicted field 77 minutes later. The size of the predicted field, where the asteroid was expected to be, was 0.17% of the total field that we searched over two nights. According to Harris (2011), we expected to detect only ∼1.2 faint asteroids over the two observing nights. If we scale the probability linearly with the area of the observed field, the chance for detecting a different asteroid in the expected field was 0.0017 × 1.2 ∼ 0.002. Therefore, the second observation is a confirmation of the first detection with 0.002 false positive probability, which is further reduced by the consistency of estimated velocities.

We identified the reference stars S5 and S8 in the 2MASS catalog as 2MASS J00512128-0000319 (S5) and 2MASS J00511743+0000477 (S8) for the second observation of the asteroid. Based on these, we again estimated the plate scale to be 0324 per pixel and determined a clockwise rotation between our camera frame and the ECS of 065.

Using the 30 s data from the second observation, we estimated the velocity of the asteroid to be [651, 0003 day⁻¹], which is within 2σ of the velocity of the initial detection, [638, −023 day⁻¹] (the uncertainty is σ_v = 016 day⁻¹). Assuming a uniform distribution of asteroid velocity within a circle of radius v₀, which is roughly the velocity of this asteroid, i.e., v₀  ∼  64 day⁻¹, the chance for a random asteroid to be within 2σ of the velocity of the initial detection is $\sim \pi (2\sigma _v)^2/(\pi v_0^2) \sim (2\sigma _v /v_i)^2 \sim 0.0025$. With the consistency of velocity, we conclude that the probability for the two detected signals being from two different asteroids is 0.002 × 0.0025 = 5 × 10⁻⁶.

The photon fluxes for the two detections are 207 ± 18 photons s⁻¹ (first) and 213 ± 21 photons s⁻¹ (second, 77 minutes later), respectively. This consistency in the photon fluxes shows that the asteroid did not change its brightness significantly after rotating over 77 minutes. This further reduces the false positive probability that the two detections are for two different two asteroids. Based on both observations, the apparent magnitude of the asteroid is 23.1 ± 0.1.

Using star S5 (2MASS J00512128-0000319) as a reference, the astrometry of the asteroid is RA 0:51:25.85 decl. −0:0:59.35 at UTC 20130911:T17:08:10, the mid-epoch of the second observation. The proper motion of stars during one night of observing time is negligible, thus the angular distances between stars can be assumed to be constant. Figure 8 shows the sidereal background where we conducted the two observations. Reference stars (S1, S2, S3, S4), and (S5, S6, S7, S8) can be easily identified in the field for the first and second observations of the asteroid. As discussed previously, these reference stars offer a calibration of our camera frame relative to the ECS. With two observations, we found that the asteroid moved 1219'' along R.A. and −4'' along decl. during 4626 s, providing an estimate of proper motion of [632, −002 day⁻¹] with an accuracy of 0005 day⁻¹.

Zoom In Zoom Out Reset image size

Assuming an asteroid's typical speed of ∼10 km s⁻¹ relative to the Earth, we estimate that the distance of the object is approximately at a distance d ∼ (10 km s⁻¹)/(632 day⁻¹) ∼8× 10⁶ km ∼ 20 lunar distances, giving an H magnitude of 29. Assuming an albedo of 0.15 (albedo of the moon), we estimate the size of the asteroid to be 8 m.

A discovery of an asteroid requires three observations to uniquely determine its orbit. Our two night survey was designed mainly for detecting faint objects, not discovery, because the expected outcome is only a couple of asteroids. We were not able to claim a discovery of the NEA because it was only observed at two different epochs on the same night. Even though detecting a single asteroid does not provide much statistics to improve the models for NEA population distribution, the consistency between the actual and expected number of detections in the H-magnitude range of 28–31 gives us confidence on the available distribution of the NEA population (Harris 2011). A 10 times longer survey would have the potential of reducing the uncertainties, which is currently a factor of two to three, in Harris's population model of the 10 m class NEAs.

On 2013 April 3 we observed two known asteroids, 2009BL₂ and 2013FQ₁₀. We compare the astrometry obtained by using synthetic tracking and traditional techniques relying on long-exposure images, which is simulated by co-adding short-exposure images without a prior shifting. We demonstrate the improvement of astrometric precision resulting from the improved S/N and from cancellation of the effects due to atmosphere and imprecise telescope pointing that nominally impact relative astrometric measurements. Because an asteroid moves at approximately a constant velocity during the 15 minutes of observation, we compute the relative astrometry between the asteroid and background stars and compare it with a motion with constant velocity to determine errors in the astrometry.

The top and middle plots in Figure 9 show the de-trended temporal variation of the location of a 16th magnitude background star (empty dots/squares) and the location of asteroid 2009BL₂ (solid dots/squares) (of apparent magnitude 18.5) relative to the camera frames (telescope pointing), for declination (top) and right ascension (middle), respectively. In the field, the star locations vary due to atmosphere effect and imprecise telescope pointing with an RMS of 100 mas. We de-trended (by removing an overall constant and a linear trend from the data) both star locations (there is a small drift in the tracking) and asteroid locations to show the common temporal variation of the locations. The bottom plot in Figure 9 displays the de-trended relative astrometry of the asteroid with respect to the 16th magnitude star, which is mostly linear (not shown, de-trended to see the residual) with less than 20 mas residual RMS, which is much smaller than the 100 mas level variation. This shows the cancellation of the common temporal variations shown in the top and middle plots between the asteroid and the star. However, this would not happen if we co-added the images to synthesize a long exposure.

Zoom In Zoom Out Reset image size

We fit the relative astrometry of the asteroid with respect to the background stars to a linear motion model and use the RMS of the residuals as a measure of the astrometry error. Figure 10 shows the astrometry error RMS as a function of the integration time. The three pairs of curves show the RMSs of the astrometry errors for using synthetic tracking with short-exposures (dots) and a single exposure with the duration being the integration time and processed using a conventional centroid fitting to a 2D Gaussian PSF (diamond) or a Gaussian trail function (squares) as described in (Vereš et al. 2012). The solid and empty markers denote declination and right ascension, respectively. Compared with the estimated photon-noise-limited RMS error (dot-dashed line) using Equation (7), as the integration time increases, the synthetic tracking astrometry error is close to be photon-noise-limited,⁹ and decreases as the inverse of the square root of the integration time. However, the astrometric errors for the traditional long exposure approach does not decrease as the integration time goes beyond 30 s because the image of the asteroid is streaked by more than 1'' (size of the seeing-limited PSF) and the error is no longer dominated by the photon noise, but by the effects due to atmosphere and imprecise telescope pointing, which is no longer common between the streaked asteroid images and the non-streaked background star images.

Zoom In Zoom Out Reset image size

Figure 11 displays the residuals of relative astrometry, after fitting to a linear motion, for 80 s integration time. This calculation is done for the cases of using synthetic tracking (circles), long exposure with a 2D Gaussian fitting (diamonds), and long exposure with a Gaussian trail fitting (squares), respectively. The synthetic tracking using short-exposures enables us to achieve mas-level astrometric precision, close to being photon and background noise limited as in relative stellar astrometry (Boss et al. 2009). On the other hand, by using long exposures, the errors due to atmosphere and imprecise telescope pointing are uncommon between the asteroid and stars and thus lead to tens of mas errors in relative astrometry. These are the result of using a 200 inch (5 m) telescope; the atmospheric effect for using 1 m class telescopes is expected to be even larger.

Zoom In Zoom Out Reset image size

We obtained similar results for the asteroid 2013FQ₁₀, which we observed for 300 s. Figure 12, which is the same plot as Figure 10 for 2013FQ₁₀, shows the improvement of astrometry as we integrate more 2 Hz frames. The residual errors are mainly due to photon and background noise as estimated using the empirical formula Equation (12). The error becomes single-digit mas after integrating over 100 s. Again, if we were to use long-exposure (i.e., 30 s) images, even with a state-of-the-art fitting of a Gaussian trail (square markers), the error would still be as large as tens of mas and grow with the integration time.

Zoom In Zoom Out Reset image size

Synthetic tracking is especially valuable for detecting very small and fast-moving asteroids. The vast majorities of these objects are so small and move so rapidly that follow-up observations have to be planned from the start. Finding the object one hour later when it has moved 1000'' is possible only because our single observation has a coarse velocity measurement. For a really faint object with S/N ∼ 7, a coarse velocity measurement with uncertainty ∼03 day⁻¹ would not let another observatory find it easily one day later. In the recent observing run, we scanned each part of the sky twice. However, a large portion of that time was wasted. Ideally, we would have the GPU software running to detect the NEAs in near real time so that the second confirming observation could be made within the time frame (a few hours) later before the uncertainty of the predicted location of the asteroid becomes larger than the FOV. Once confirmed, we should spend significantly more than 30 s to get astrometry approaching 50 mas so that three observations spaced a few days apart would let us derive an orbit where the object could be observed again at the next apparition (J. Giorgini 2013, private communication). Objects with H ∼ 28–30 moving at 6°day⁻¹ are detectable only on medium-to-large telescopes with synthetic tracking cameras, and even then, only for a few weeks or so before they are too faint to be detected. Therefore, it is crucial to have a strategy as discussed above to discover the object within its time window for observation.

The astrometric accuracy of the synthetic tracking technique enables us to determine the area-to-mass ratio of an NEA by measuring the effect of the solar radiation pressure on its orbit (Micheli et al. 2014). As an example, we simulated an NEA with a semi-major axis of 1.2 AU that comes within ∼1 lunar distance of the Earth. This NEA would be brighter than 24 mag for ∼60 days. Astrometry during these 60 days would show the difference between a purely gravitational orbit and one that is influenced by solar radiation pressure on a 10 m sized NEA with a density of 1 gm cm⁻³. The left panel in Figure 13 shows the apparent magnitude of the object from 100 days before to 100 days after the closest approach. The top right panel in Figure 13 shows how the NEA's orbit would change if solar radiation pressure were turned on. However, the modified orbit with solar pressure could be partially mimicked by a different orbit without solar pressure. The bottom right panel in Figure 13 shows the difference between the best-fit zero-solar-pressure orbit to the solar pressure orbit, which is the astrometric signal due to solar pressure that cannot be explained by a pure gravitation-driven motion. ¹⁰ If the astrometric position of the NEA during the "close" encounter and the subsequent ∼50 days are measured to <10 mas at ∼20 epochs concentrating on the ∼14 days when the NEA is closest to Earth, the area–mass ratio of this NEA can be measured at an S/N ∼ 7.

Zoom In Zoom Out Reset image size

The measurement of area–mass ratio by optical astrometry is made easier when the NEA passes very close to the Earth. Micheli et al. (2014) measured the area–mass ratio of 2011MD with astrometry that was only accurate to ∼02–04 because this NEA came within a couple of Earth radii to the Earth. With more precise astrometry in our example above, measuring the area–mass ratio becomes possible for NEAs with much larger impact parameters. With the additional radar ranging and size data, milliarcsecond optical astrometry may be able to measure mass and density for a significant number of NEAs per year.

Because astrometric measurements of NEA positions are absolute position measurements, to achieve a 1 mas accuracy we need a star catalog with accuracy better than 1 mas. Currently the most useful catalog for asteroid astrometry is the USNO's UCAC catalog of ∼50 million stars, accurate to ∼50 mas, anchored by the stars in the it Hipparcos and it Tycho catalogs. ESA's it Gaia mission will begin science operations in 2014 June. The mission's ultimate goal is to achieve ∼10 μas accuracy. The release of an initial catalog, expected in the fall of 2015, will provide an accuracy better than 1mas for all ground-based NEA astrometry.