Space Object Detection Algorithm Leveraging Absolute Photo-Detector Calibration

This paper introduces a new space object detection algorithm that is designed to process image data taken from astronomical telescopes for the purposes of finding sources of optical radiation in space. Specifically, the algorithm is designed to find unresolvable space objects or objects that possess an angular size that is too small to appear as anything, but a point source as viewed through the telescope conducting the search. The proposed approach involves calibrating the image data into units of photoelectrons and then executing an estimation algorithm to compute the strength of the hypothetical sources in the image. A Likelihood Ratio Test (LRT) is then implemented to make a determination if the hypothetical sources are classified as space objects or not. The proposed algorithm is demonstrated to achieve a higher probability of detecting unresolvable objects than the matched filter, which is still the state-of-the-art approach for finding optical sources in astronomical images. The new approach involves a pre-processing step where the amplitude of the optical source in a given test location is estimated under the hypothesis that at optical source exists at that location. The median filter is used to estimate the background level in the vicinity of the test location. Once these parameters are estimated, A likelihood ratio test is used to determine whether an object is present at the test location. The new algorithm is tested against the matched filter detector using two sets of measured short exposure data of the star Polaris and two stars in its vicinity taken with a small telescope. Receiver Operating Characteristic (ROC) curves are produced for the two detection schemes showing that the new algorithm out-performs the old one with a difference of 10 percent in the probability of detection, which is demonstrated to be statistically significant in these experiments with confidence as high as 90%.


Introduction
Detection of objects in space is an ever increasingly important function of astronomical telescopes and the astronomical community. Useful space situational awareness is important for planetary protection from asteroids as well as mission planning for earth orbiting satellite programs and safety of manned space missions. The future need for space object detection will only grow as access to space becomes more affordable. The economics of the future space economy are still in doubt; however, so making space surveillance as affordable and cost effective remains a priority [1,2].
Design and construction of space surveillance equipment remains expensive, whether those assets are on the ground or in space itself. The research explored in this paper is designed to economize space surveillance by improving the ability of existing sensors to accomplish the space object detection mission. For this reason, this paper will not discuss new advances in device or optics technology for accomplishing space surveillance but will instead only discuss research that can improve an existing system's ability to detect objects in space through image processing. Additionally, with the advent of radar systems capable of detecting objects in low earth orbit, the research discussed in this paper will focus on improving the detection of objects in geosynchronous and geostationary orbit, which represent one of the most crowded regions where satellites operate.
Image processing algorithms for accomplishing space object detection utilize different phenomenology to accomplish the detection. Some types of algorithms are designed to process long exposure images of objects in space while the telescope is rate tracking, thus making all stars appear as point sources and objects not moving at the sidereal rate appear as streaks of light. This type of algorithm attempts to accomplish streak detection, which suffers efficiency problems as the light from the desired object is spread over many pixels, thus decreasing its visibility [3]. A more efficient method of telescope operation is to fix the telescope and allow the stars to streak across the detector during a long exposure. In this mode of operation, objects in geosynchronous or geostationary orbit will not streak and will have their photons integrate in fewer pixels, thus achieving a higher signal to noise ratio in those detectors [4]. For this reason, this research will endeavor to detect objects that appear as point sources to the observing telescope. A related area of research involves how to process collections of observations to make detection decisions [5,6], but this research focuses on how to process a single observation as multiple observations are not always available in every case.
Detection of objects that appear as point sources has been a subject of research for decades. MIT Lincoln Labs developed the some of the first algorithms for detecting such objects with what is called by some as the point detector [7]. This type of algorithm endeavors to find a space object by looking at a single pixel and normalizing its response by the local standard deviation of the measured intensities around the test location after subtracting the local background. The resulting signal is compared to a threshold set to achieve the desired probability of false alarm. A similar approach has been implemented with great success using a matched filter approach with the same kind of background subtraction and normalization scheme [4,8]. Other attempts have been made at improving space object detection over what is possible with a matched filter, by employing collections of matched filters or exploiting the statistics of space object occurrences [9.10] These techniques however do not improve space object detection in cases where the image field is not under sampled and rely on the basic matched filter approach, just with more filters. The matched filter approach is a state-of-the-art approach for accomplishing space object detection and will be used as the method to compare to in this research.
The proposed space object detection algorithm explored in this paper uses a Log-Likelihood Ratio Test (LRT) based on Poisson statistics similar to that proposed by Polig [3]. The main difference in what is proposed here is that the brightness of the hypothetical object to be detected is estimated from the data and the natural logarithm is not approximated by a Taylor series expansion. The proposed technique also requires the use of a calibration algorithm [11] to convert the digital counts collected by the photodetector into units of photoelectrons that exhibit Poisson statistics. The digital counts are discrete and non-negative, but do not display Poisson statistics due to the detector gain and bias present in the system. This research will demonstrate how absolute calibration can be used in conjunction with an LRT to improve the ability of a telescope to detect objects in space over what can be achieve through a matched filter approach. Section 2 of this paper introduces the hybrid estimation and detection algorithm. Section 3 showcases the data sets used in this study. Section 4 displays the results obtained from applying the new algorithm and the matched filter to the data sets. Conclusions of the research are discussed in Section 5 as well as future research possibilities.

Algorithm Design
The algorithms utilized in this paper for processing the image data comprise both an estimation scheme for the model parameters as well as a detection algorithm. The first is the estimation scheme, which is used to recover the background level as well as the hypothetical amplitude of the source as the test location. The second is the Likelihood Ratio Test (LRT) used to make the detection decision. The estimation step begins with the use of the median filter to determine the level of the background in the local area around the test location. The next step in the estimation procedure involves the estimation of the amplitude of a star that might exist at the test location through the use of an Expectation-Maximization approach. This new algorithm is compared to that of a baseline approach previously reported in the literature [8].

Estimation Algorithm
The number of photoelectrons measured at each detector in the array is assumed to be a Poisson random variable with a mean equal to the sum of the mean number of photons from the objects in the field of view and the background.
The data from these sources is modeled as separate Poisson random variables such that the data from the objects, d k (x,y) has the mean show in Equation (1).
In this Equation, (x,y) are coordinates in the detector plane in an M-by-M pixel area centered at the test location (the test location has coordinates (M/2+1,M/2+1), γ k is the number of photoelectrons from optical sources in the k'th frame of image data out of a total of J frames and h k is the PSF in the k'th frame. The second set of data. d B (x, y,k) has a mean equal to the background level as shown in Equation (2).
In this equation, B k is the background level in the k'th frame of data, which is assumed to be uniform in the field of view on average and is computed using a median filter. The measured data, , is also Poisson and is equal to the sum of all both data sets as shown in Equation (3).
An Expectation-Maximization approach is adopted for computing the amplitude of an unknown object from a single frame of camera data [12][13][14]. This procedure begins by identifying the joint log-likelihood of the two data sets, d k (x,y) and d B (x, y,k). Since these data are assumed to be Poisson and both statistically independent from one another as well as statistically independent from pixel to pixel, the joint log-likelihood, L(g k ), can be expressed using the following equation, which assumes that the PSF's sum to one and the data in each pixel in every frame are statistically independent from one another: The derivative of L(g k ), with respect to γk is computed as shown in Equation (5).
A second derivative yields a strictly negative result, implying that a solution for Equation (5) when the derivative is set equal to zero will produce a maximum of the log-likelihood function. The solution for γ k that maximizes L(g k ) is shown in Equation (6). The conditional mean of the estimate given the measured data is computed to show the form of the estimate obtained from the Expectation-Maximization process.
The expectation step shown in Equation (6) is necessary since the data described in Equation (1) is not directly observable but is contained within the observed data, , as shown in Equation (3). The conditional expectation has been solved by Shepp and Vardi and is substituted into Equation (6) to yield the following result [15].
In this equation a variable with the old superscript denotes a quantity computed with estimates obtained from the previous iteration of the EM algorithm. γ k new denotes the new estimate obtained from the EM algorithm for the number of photons present in the k'th data frame. In this way the algorithm is implemented by the user providing a starting estimate which is updated from iteration to iteration via Equation (7). One obvious way to initialize the estimate for γ k is to subtract the median value, B k , from the image, , , and sum the result over all pixels.

Likelihood Ratio Test
The test for determining whether an object is present at the test location is an LRT. To derive the form of the test, the probability of the data in the M-by-M pixel area around the test location must be surmised for both the hypothesis that a source is present at the test location as well as the hypothesis that a source is not present at the test location. The test is only conducted if the estimated value of γ k new exceeds a threshold of γ thresh photons. The natural log of the probability of the data given there is an optical source at the test location is shown in Equation (8).
The natural log of the probability of the data given there is no optical source present is shown in Equation (9).
Subtracting the natural logarithms in Equation (9) from that found in Equation (8) yields the natural logarithm of the LRT, Q. If Q is greater than a threshold, t, we say that H 1 is true and if it is not, we say that H 0 is true.
Q is the test statistic for the new algorithm and is notably different than the test statistic used when employing a matched filter detector, which is denoted as C and is shown in Equation (11) for reference. 1 1

Experimental Data
The data utilized in this research is taken with an 8-inch reflecting telescope with a 2-meter focal length that has a Cassegrain design. The aperture is stopped down with a mask that admits light in a 2-inch diameter circle within the entrance pupil of the optic at a position away from the secondary mirror. This creates an optical system that has a 2-inch diameter with a 2-meter focal length, achieving a F# of 40. The average wavelength of the light entering the telescope is approximately 500 nm. The CCD camera attached to the telescope is from a Roper Cascade 512B scientific camera which has a 16 micro-meter pixel pitch.
Two sets of image data were taken containing a collection of optical sources with integration times of 10 and 100 ms. The telescope was pointed at Polaris and regions of the image were identified as containing optical sources and other regions were shown to possess no discernable optical sources. Figure  1 shows an average of 1000 frames of data taken at 100 ms integration time. Polaris is the brightest object in the field of view.  Two Regions were chosen to provide data where weak optical sources reside. A region was also chosen where there appear to be no optical sources present in the average data. Figure 2 shows an area that corresponds to the location of Polaris B. It is not visible in Figure 1, but its image can be viewed when isolated by itself as shown in Figure 2. Figure 2 represents an average of 1000 frames of image data. Figure 3 shows a typical individual frame. Figure 4 shows the 1000 frame average of an area near Polaris A that appears to be free of optical sources. Figures 5 and 6 show another area used for testing where a weak optical source appears to be present in the 10ms data set. All 5 of these data sets were calibrated using an algorithm known as the SANUC calibration method. This method utilized 1000 frames of background data taken of the sky area shown in Figure 4. Using 1000 frames of the 100ms data and 1000 frames of the 10ms data, this calibration technique was able to compute that the median pixel gain was 20.7 digital counts per electron, the average sky background photo-electron count was 3.0 photo electrons and the digital bias level was 146 counts.  All the data frames were calibrated into units of photoelectrons by subtracting off the local mean, setting a negative value to be zero and then dividing the remaining digital counts by the gain of 20.7. A value of 3.0 photoelectrons were then added to all the pixel values to account for the background light.
This calibrated data was then processed by both the matched filter detector and the proposed detection algorithm with the results shown in the next section. In order to accomplish this a shape for the Point Spread Function (PSF) otherwise known as the system impulse response function must be given. Figure 7 shows the shape of the PSF used in this study that was obtained from using Polaris A as a guide star. This star is chosen because it is close to the optical source data and is comparatively bright.

Results
1000 frames of image data collected on an object in the vicinity of Polaris A as shown in Figure 6 was processed by the algorithm proposed in Section 2 of this paper. This data contains an object centered at the middle of the image which is the test location for this study. Since we know there is an object present (as seen in the 1000 frame average at that location) the output of the LRT defined in Section 2 serves to characterize the behavior of the detection algorithm when an object is present at the test location   Figure 8 shows the histogram of outputs from the LRT when processing the data of the dim object shown in Figure 5. This data set of 1000 LRT outputs from equation (10) was fit to a Gaussian distribution using the MATLAB © NORMFIT function. The mean and standard deviation are computed as 2082.5 and 486.2 respectively and a lower bound on the mean was found to be 2074.5 and a lower bound on the standard deviation was computed to be 480.7.
These bounds were computed with 70% confidence, meaning that there is a 70 percent chance that the mean is above the reported bound and the standard deviation is below the bound. The same procedure is followed with data collected of the area shown in Figure 4. Figure 9 shows a histogram of the output of the algorithm when processing an area that appears to not have any optical sources present in the long-term average. This data set was also fit to a Gaussian distribution using the MATLAB © NORMFIT function. The mean and standard deviation are computed as 1419.7 and 356.5 respectively and an upper bound on the mean was found to be 1425. 6 and an upper bound on the standard deviation was computed to be 360.8.
These statistics were used to compute the probability of detection and probability of false alarm of the proposed algorithm. Using a normal distribution, the probability of detection is computed by choosing a threshold and then integrating the probability of the normal distribution above the threshold. This is mathematically equivalent to taking one minus the Cumulative Distribution Function (CDF) at the threshold value.   Figure 6 of the object in Figure 5.
A lower bound on the probability of detection can also be computed by taking the 70% confidence bound values reported earlier and using the CDF method just described to compute the probability of detection using the lower bound of the mean and standard deviation of the detector outputs. By choosing the lower limit of the mean and lower limit on the standard deviation, the probability of detection is bounded using these 'worst case' values. If the mean was higher and the standard deviation was higher, then the probability of detection for a given threshold (over the mean value) would be lower than the that computed using the nominal values for the mean and standard deviation.
The same process is utilized for computing the probability of false alarm using the outputs from the 1000 frames of data of the empty sky fed into the proposed algorithm. In this case because the sky is empty, we are computing the probability of false alarm. To obtain the 'worst case' for this statistic the upper bound on the mean and upper bound on the standard deviation are used to compute the upper bound on the probability of false alarm.   Figure 6 of the object shown in Figure 5.
With a lower bound on the probability of detection and an upper bound on the false alarm rate computed for a given threshold, the Receiver Operating Characteristic (ROC) curve can be graphed by plotting the probability of detection vs. the probability of false alarm for each threshold value. The range of threshold values was chosen to be between zero and 18,000, which spans the entire range of possible outputs of the proposed algorithm in this case with much greater than 99.9999% probability.
This same procedure is applied to the output of the matched filter algorithm. The difference is that instead of computing the 'worst case' performance of the algorithm with 70% confidence, now the 'best case' is computed with 70% confidence again using the MATLAB © NORMFIT function. Figure 10 shows the histogram obtained from the matched filter output when the images used to generate Figure 6 are input. The mean and standard deviation are computed as .0491 and .04 respectively and an upper bound on the mean was found to be .0498 and a lower bound on the standard deviation was computed to be .0397.  The same procedure is followed with data collected of the area shown in Figure 4. Figure 11 shows a histogram of the output of the algorithm when processing an area that appears to not have any optical sources present in the long-term average. This data set was also fit to a Gaussian distribution using the MATLAB © NORMFIT function. The mean and standard deviation are computed as -.0021 and .03 respectively and a lower bound on the mean was found to be -.0026 and a lower bound on the standard deviation was computed to be .0295. Figure 12 shows the ROC curves produced by the proposed LRT detection algorithm for the nominal case and the 'worst case' estimates of the mean and standard deviations of the outputs. This figure also shows plots of the ROC curve for the nominal values of mean and standard deviations of the outputs obtained from the matched filter detector as well as the ROC computed using the 'best case' with 70 percent confidence. The 70 percent confidence interval was chosen as this is the point at which the 'worst case' of the proposed algorithm and the 'best case' of the matched filter are approximately equal. This signifies that the proposed detection algorithm achieves approximately a 5 percent improvement in the probability of detection with at least 70 percent confidence.
The exact same procedure is used to compute the performance of the algorithms on the data shown in Figure 3 of Polaris B. The same set of empty space is used for computing the probability of false alarm. Figure 13 shows the ROC curves for nominal values of the mean and standard deviation of the detector outputs as well as the 'worst case' ROC for the proposed detector and the 'best case' ROC for the matched filter detector. The Pfa region between 0 and 0.1 is shown here on the x-axis as this is the region in which the new detector has the greatest performance advantage over the matched filter and lower false alarm probabilities are generally desirable for practical operation. The red curve showing the nominal performance is just over the line representing the 90% lower bound of performance.
The same set of data without employing the photo-calibration step was processed in the exact same way as before and produced the ROC curves seen in Figure 14, showing how important the absolute calibration step is.

Conclusions
The results presented in the previous section show that the combination of a calibration algorithm together with an LRT based on Poisson statistics can improve the probability of detecting space objects over what is achievable with a matched filter. The improvement is as much as 10 percent and is shown to be statistically significant with 90% confidence. The application of the new technique does not require expensive of time-consuming calibration work, since the sky background is used to compute the calibration parameters. All that is required is that two sets of data taken at different integration times be available to complete the calibration step. Failure to employ the calibration step produces results that are inferior to the matched filter as shown in Figure 14. This is due to the reality that although the photodetector is measuring a Poisson random variable, the system gain and bias skew the statistics in such a way as to reduce the effectiveness of an LRT that is based on Poisson statistics. Calibrating the data back into units of photoelectrons alleviates this problem.
The method proposed in this paper does require knowledge of the shape of the PSF as does the matched filter algorithm. Future work could include a sensitivity analysis to determine which algorithm requires more exact knowledge of the PSF. Another topic for future exploration could include an algorithm that adaptively determines the PSF while detecting objects in space. Such a technique would potentially be more robust as it would deal with the fluctuations in atmospheric seeing that inevitably occur during long observation periods. It would also deal with telescope aberrations that are field dependent and change across the sensor field of view, which is a serious problem for synoptic survey telescopes.