An Adaptive Real-Time Detection Algorithm for Dim and Small Photoelectric GSO Debris

Geosynchronous orbit (GSO) is the ideal orbit for communication, navigation, meteorology and other satellites, but the space of GSO is limited, and there are still a large number of space debris threatening the safety of spacecraft. Therefore, real-time detection of GSO debris is necessary to avoid collision accidents. Because radar is limited by transmitting power and operating distance, it is difficult to detect GSO debris, so photoelectric detection becomes the mainstream way to detect GSO debris. This paper presents an adaptive real-time detection algorithm for GSO debris in the charge coupled device (CCD) images. The main work is as follows: An image adaptive fast registration algorithm and an enhanced dilation difference algorithm are proposed. Combining with mathematical morphology, threshold segmentation and global nearest neighbor (GNN) multi-target tracking algorithm, the functions of image background suppression, registration, suspected target extraction and multi-target tracking are realized. The processing results of a large number of measured data show that the algorithm can detect dim geostationary earth orbit (GEO) and non-GEO debris in GSO belt stably and efficiently, and the processing speed meets the real-time requirements, with strong adaptive ability, and has high practical application value.


Introduction
The GSO belt is 42,164 kilometers away from the center of the earth, 400 kilometers in radial thickness and delimited by ±15 • in latitude [1]. It is an ideal orbit for communications, navigation, meteorology and other satellites [2,3]. According to data released by the U.S. Strategic Command, as of July 2018, there were about 1200 objects in GSO orbit and about 400 satellites in service [1,4]. In addition, according to data from the European Space Agency, as of January 2019, there were about 1950 spacecraft in normal working state in all Earth orbits. At the same time, it was estimated that about 129 million debris were larger than 1mm, while only about 22.3 thousand objects were catalogued in the database and could be continuously tracked. Although the density of debris in GSO is much less than that in low Earth orbit (LEO), there is no debris clearance mechanism similar to atmospheric damping in GSO; debris will persist for millions of years and pose a greater threat to spacecraft [5,6]. Maintaining space security and accurately monitoring space debris is a part of space situational awareness [7]. Accurate monitoring of GSO objects can judge the operational status of satellites, evaluate the probability of collision events, and contribute to celestial mechanics, help people understand orbit evolution and motion characteristics [8,9]. Therefore, it is meaningful to study the detection algorithm of GSO debris [10,11]. The detection methods of space objects mainly include radar detection and optical detection. Microwave radar has the advantages of high accuracy and Top-Hat transform algorithm, image registration based on the adaptive fast registration algorithm, stellar suppression based on the enhanced dilation difference algorithm, GSO debris enhancement and target segmentation based on inter-frame correlation and threshold segmentation technology. The target tracker uses GNN multi-target tracker and can estimate the inter-frame interval adaptively to improve the performance of the tracker.
The rest of the paper is organized as follows. In the second section, each component of GSO debris extractor is introduced, and the adaptive fast registration algorithm and the enhanced dilation difference algorithm proposed in this paper are emphatically introduced. In the third section, the GSO debris tracker is introduced. The principle of GNN tracker and the adaptive inter-frame interval estimation technique proposed in this paper are described. The fourth section contains the introduction of test data and the processing results of the algorithm. The fifth section gives the conclusions and the prospects of future work.

GSO Debris Extractor
This section will introduce each component of GSO debris extractor, including background suppression, image registration, stellar suppression, GSO debris enhancement and target segmentation.

Background Suppression Based on Morphological Top-Hat Transform
The background of CCD star image will be disturbed by various noises, including photon noise, environmental noise, readout noise and so on. These interferences seriously affect the effect of subsequent image processing steps, so they must be suppressed. In view of the strong image adaptive ability of mathematical morphology, parallel processing and fast computing speed [34], this paper uses morphological Top-Hat transform algorithm for background suppression. Top-Hat transform background suppression algorithm has been widely used because of its advantages in background suppression and real-time operation. The formula of top-hat transform is as follows: where f is the original input image and b is the structural element. f • b stands for open operation of image f using structural element b, which is equivalent to a non-linear low-pass filter. It can remove bright details smaller than the size of structural element, while preserving the overall gray value of the image, and the larger bright areas are almost unaffected. The background of the image can be estimated by using appropriate structure element b to perform open operation. The g obtained by subtraction is the target image that suppresses the background. Whether the selection of structural element is appropriate or not will have a great impact on the processing efficiency and effect of Top-Hat transformation. According to the principle, the size of structural element should be larger than the targets size. By analyzing the measured data, the CCD star image mainly includes GEO debris, non-GEO debris and stars besides background noise, as shown in Figure 1.
In the image data used in this paper, the size of GSO debris is between 2 × 2 and 20 × 20 pixels. In order to reduce the computational complexity as much as possible while guaranteeing better background suppression effect, the size of the structural element used in this algorithm is 25 × 25, and it is further decomposed into flat structural elements of 1 × 25 and 25 × 1 [35]. The erosion and dilation operations of flat structural elements involved in the open operation are simplified to the minimum and maximum filters with the support domain of structural elements as the window respectively, which greatly reduces the algebraic computation. In order to further improve the processing speed and take into account the portability of the algorithm, flat structure elements of 1 × 25 and 25 × 1 are decomposed into small size structure elements of 1 × 3 and 3 × 1 [36], as shown Figure 2. If the image size is m × n, the computational complexity of the background suppression algorithm is O(mn).  Figure 3 shows the CCD star images before and after background suppression. It can be seen that the background is effectively suppressed. remove bright details smaller than the size of structural element, while preserving the overall gray value of the image, and the larger bright areas are almost unaffected. The background of the image can be estimated by using appropriate structure element b to perform open operation. The g obtained by subtraction is the target image that suppresses the background.
Whether the selection of structural element is appropriate or not will have a great impact on the processing efficiency and effect of Top-Hat transformation. According to the principle, the size of structural element should be larger than the targets size. By analyzing the measured data, the CCD star image mainly includes GEO debris, non-GEO debris and stars besides background noise, as shown in Figure 1.   In the image data used in this paper, the size of GSO debris is between 2 2 × and 20 20 × pixels. In order to reduce the computational complexity as much as possible while guaranteeing better background suppression effect, the size of the structural element used in this algorithm is 25 25 × , and it is further decomposed into flat structural elements of 1 25 × and 25 1 × [35]. The erosion and dilation operations of flat structural elements involved in the open operation are simplified to the minimum and maximum filters with the support domain of structural elements as the window respectively, which greatly reduces the algebraic computation. In order to further improve the processing speed and take into account the portability of the algorithm, flat structure elements of1 25 × and 25 1 × are decomposed into small size structure elements of1 3 × and 3 1 ×

Adaptive Fast Image Registration
Because GSO debris are typical dim and small targets, in order to improve the probability of detection, an effective way is to accumulate energy through multi-frame images to improve the signalto-noise ratio (SNR) of the target. In this paper, local SNR is adopted, which is defined as follows: where, T f and B f are the gray mean of the target area and background area respectively, and σ is the gray mean variance of the background area. In order to prevent the interference between target

Adaptive Fast Image Registration
Because GSO debris are typical dim and small targets, in order to improve the probability of detection, an effective way is to accumulate energy through multi-frame images to improve the signal-to-noise ratio (SNR) of the target. In this paper, local SNR is adopted, which is defined as follows: where, f T and f B are the gray mean of the target area and background area respectively, and σ is the gray mean variance of the background area. In order to prevent the interference between target and background gray statistics, they are separated by isolation area. The definitions of different regions are shown in Figure 4.
(a) (b) Figure 3. CCD star image before and after background suppression: (a) before background suppression; (b) after background suppression.

Adaptive Fast Image Registration
Because GSO debris are typical dim and small targets, in order to improve the probability of detection, an effective way is to accumulate energy through multi-frame images to improve the signalto-noise ratio (SNR) of the target. In this paper, local SNR is adopted, which is defined as follows: where, T f and B f are the gray mean of the target area and background area respectively, and σ is the gray mean variance of the background area. In order to prevent the interference between target and background gray statistics, they are separated by isolation area. The definitions of different regions are shown in Figure 4. The premise of multi-frame image processing is image registration. GSO debris detection uses ground-based telescope, and the tracking mode of the telescope is the target tracking mode. Because the GSO debris moves slowly relative to the ground-based telescope, the target tracking mode here actually means that the telescope field of view is fixed in the target sky area, that is, the telescope is stationary when it is working normally. Therefore, the movement between the two frames is basically caused by the motion of the stars. The motion of stars is mainly caused by the rotation of the earth, the rotation speed of the earth is basically constant, and because the stars are very far away from the telescope. So, on the CCD image, there is only the translational motion of the stars, and there is no rotation. Overall, for the registration of CCD star images of GSO only involves translation The premise of multi-frame image processing is image registration. GSO debris detection uses ground-based telescope, and the tracking mode of the telescope is the target tracking mode. Because the GSO debris moves slowly relative to the ground-based telescope, the target tracking mode here actually means that the telescope field of view is fixed in the target sky area, that is, the telescope is stationary when it is working normally. Therefore, the movement between the two frames is basically caused by the motion of the stars. The motion of stars is mainly caused by the rotation of the earth, the rotation speed of the earth is basically constant, and because the stars are very far away from the telescope. So, on the CCD image, there is only the translational motion of the stars, and there is no rotation. Overall, for the registration of CCD star images of GSO only involves translation transformation, and often lacks information such as exposure time, telescope position, and optical axis pointing. In view of this, an adaptive fast registration algorithm is proposed in this paper. The correlation coefficient is a statistical index to reflect the degree of close correlation between variables. Obviously, when the correlation coefficient of two images is the largest, it can be judged that two images have been registered. Therefore, our goal is to find a translation vector u shi f t = x shi f t , y shi f t to maximize the correlation coefficient between the shifted image I shi f t and the reference image I f ixed . The formula for calculating the correlation coefficient is as follows: where I shi f t u shi f t represents the shifted image after translation according to translation vector u shi f t , ρ u shi f t represents the corresponding correlation coefficient, E(·) represents the average value of the image. The maximum correlation coefficient can be measured by minimizing the sum of squared differences (SSD) of two images [37]. The formula for calculating the SSD is as follows: where SSD u shi f t represents the SSD corresponding to translation vector u shi f t . When the size of the image is large, the amount of calculation to get the SSD of the whole image is very large, which seriously reduces the processing speed. In order to improve the efficiency, we uniformly select n identical sub-images on the reference image and the shifted image, in which the sub-images on the reference image is fixed and unchanged, and the sub-images on the shifted image moves within a given search area, each moving step is one pixel. At each step, the SSD of the corresponding sub-images is calculated. After all SSDs corresponding to the search area are calculated, the translation vector corresponding to the minimum SSD is required. Let D denote the search area, I shi f t is as follows: where function Find(·) obtains the translation vector corresponding to SSD. Figure 5 shows how to calculate the translation vector of the i − th sub-image. For clarity, the size of D is assumed to be 5 × 5 and the translation vector u shi f t = (−1, 2) in Figure 5. Obviously, using this method, the computational complexity is only related to the number and size of selected sub-images, and the size of search area, which is independent of the size of the whole image, and the algorithm has good parallel processing characteristics. When the number of sub-images is n s , the size of each sub-image is m s × m s , and the size of search area is r s × r s , then, the computational complexity is O n s m 2 s r 2 s .  An intuitive way to obtain the translation vector of the whole image is to directly calculate the mean of the translation vector determined in all sub-images. This method is feasible in most cases, but because the size of the selected sub-image is much smaller than that of the whole image, the translation vector determined in the sub-image may have a larger error with the actual translation vector of the whole image. The results of the measured data processing also show that there are such problems. In order to overcome this problem, we noticed that the translation vector determined in most sub-images are basically the same as the translation vector of the whole image. Therefore, we classify all translation vectors determined in sub-images according to the 2-norm (Euclidean distance) of the difference between translation vectors, and calculate the translation vector of the whole image by using the class with the largest number of translation vectors. The implementation process is as follows: The first step is to calculate the Euclidean distance ( )  An intuitive way to obtain the translation vector of the whole image is to directly calculate the mean of the translation vector determined in all sub-images. This method is feasible in most cases, but because the size of the selected sub-image is much smaller than that of the whole image, the translation vector determined in the sub-image may have a larger error with the actual translation vector of the whole image. The results of the measured data processing also show that there are such problems. In order to overcome this problem, we noticed that the translation vector determined in most sub-images are basically the same as the translation vector of the whole image. Therefore, we classify all translation vectors determined in sub-images according to the 2-norm (Euclidean distance) of the difference between translation vectors, and calculate the translation vector of the whole image by using the class with the largest number of translation vectors. The implementation process is as follows: The first step is to calculate the Euclidean distance d shi f t between the translation vector u where · 2 denotes 2-norm.
In the second step, all translation vectors whose distances from u    shi f t is calculated as follows: where N shi f t and |·| represents the potential (the number of elements) of the set. Because the search step is one pixel, the accuracy of the translation vectors is one pixel in both horizontal and vertical directions, so the distance threshold T shi f t is set to √ 2. The third step is to perform the same operation on all translation vectors determined in sub-images.
In the fourth step, the translation vector of the whole image is calculated according to the U By using the above method, reliable translation vectors of the whole image can be obtained. However, considering the sudden change of exposure time in order to adapt to the change of targets' intensity when using telescopes to observe GSO debris, it is difficult to determine the search area. It is often necessary to give a larger area to satisfy different exposure times as much as possible, which leads to a significant increase in processing time. In addition, it should be noted that the exposure time will not change frequently on a large scale. Due to the instability of device state and other factors, it usually changes in a very small range, so we can reduce the search area to reduce processing time.
In order to take into account the above two situations, an adaptive search strategy is proposed in this paper. Initially, the translation vector u o is determined by using large search area D large . According to the exposure time, we set D large as a 100 × 100 rectangular area centered on the origin. The processing steps of subsequent frames are as follows: In the first step, the translation vector u k−1 determined in the previous frame is directly used to extract suspected targets. It can be considered that the search area in this step is D 1 with size 1 × 1, and the extraction algorithm will be introduced in Sections 2.3 and 2.4. If the number of suspected targets is less than the threshold T SusObj , the translation vector u k−1 is judged to be correct and the registration is completed. Also, make u k = u k−1 . Otherwise, the translation vector is judged to be wrong.
In the second step, if u k−1 is wrong, considering that the accuracy of the translation vector is one pixel in both horizontal and vertical directions, and the exposure time of the telescope often varies in a small range, then the search area is enlarged to 5 × 5 area centered on u k−1 and marked as D 2 . The translation vector u 2 is obtained by searching area D 2 , and the suspected targets are extracted by u 2 . If the number of suspected targets met the quantitative constraint, less than the threshold T SusObj and the translation vector satisfied the translation constraint R u , then the translation vector u 2 is judged to be correct and the registration is completed. Also, make u k = u 2 . Otherwise, the translation vector is judged to be wrong. How to determine the threshold T SusObj and translation vector constraint R u will be introduced later.
In the third step, if u 2 is wrong, considering that the exposure time may vary in a wide range in order to adapt to the change of target intensity, the change range is generally half of the original exposure time to twice the original exposure time. The search area is enlarged to 2|u k−1 _x| × 2 u k−1 _y area centered on u k−1 and recorded as D 3 , where, u k−1 _x and u k−1 _y represent x and y coordinates of u k−1 respectively. The translation vector u 3 is obtained by searching area D 3 , and the suspected targets are extracted by u 3 . Similarly, if the number of suspected targets is less than the threshold T SusObj and the translation vector satisfied the translation constraint R u , the translation vector u 3 is judged to be correct and the registration is completed. Also, make u k = u 3 . Otherwise, the translation vector is judged to be wrong.
In the fourth step, if u 3 is wrong, consider that the abnormal change of exposure time may be caused by abnormal reasons. The translation vector u 4 is obtained by adjusting the search area to D 4 = D large , and then suspected targets are extracted by u 4 . The threshold T SusObj is used to judge whether u 4 is correct or not. Note that the step does not need translation constraint R u , because the exposure time is likely to be abnormal at this time, so it is likely that the correct u 4 does not meet the translation constraint. If u 4 is correct, registration is completed and u k = u 4 . If u 4 is wrong, consider that it may be due to the change of the telescope field of view, or a large amount of thick clouds occlusion that cannot be registered, then stop searching and discard the frame. Also, make u k = u k−1 . The search strategy schematic is shown in Figure 6. Through the analysis of a large number of measured data, we found that the number of GSO debris in CCD star image generally does not exceed 15, and the number of false targets extracted from a single frame image by the target extractor in this paper generally does not exceed 10. Therefore, the threshold SusObj T of quantitative constraints is set to 40. For translation constraint u R , it should be noted that the processing of each frame image involves two registration processes, including the registration of the previous frame image and the current frame image, whose registration translation vector is marked as , the registration of the first two frames with the current frame image, and the registration displacement vector is marked as . Considering that the stellar translation is Through the analysis of a large number of measured data, we found that the number of GSO debris in CCD star image generally does not exceed 15, and the number of false targets extracted from a single frame image by the target extractor in this paper generally does not exceed 10. Therefore, the threshold T SusObj of quantitative constraints is set to 40. For translation constraint R u , it should be noted that the processing of each frame image involves two registration processes, including the registration of the previous frame image and the current frame image, whose registration translation vector is marked as u shi f t_1 , the registration of the first two frames with the current frame image, and the registration displacement vector is marked as u shi f t_2 . Considering that the stellar translation is generally greater than 5 pixels and less than 50 pixels during the exposure period of GSO debris optical observation, the first constraint of R u is 5 ≤ u shi f t_1 2 ≤ 50 and 10 ≤ u shi f t_2 2 ≤ 100. In addition, in order to improve robustness, we add the ratio of two translation vector modulus as the second constraint condition. Normally, the interval of changing exposure time to adapt to the change of target intensity is usually 0.5 to 2 times of the original exposure time. Considering that both horizontal and vertical accuracy of translation vectors are 1 pixel, the maximum value M_R of the ratio of modulus of two translation vectors can be determined by the following formula.
and determine the minimum value m_R of the ratio of modulus of two translation vectors by the following formula: Using the software of MATLAB, we get M_R ≈ 14.811, m_R ≈ 2.3409. The above image registration algorithm, combined with the specific application of GSO debris optical detection, takes into account most of the situations in image registration, has the characteristics of strong adaptive ability, fast computing speed and parallel processing.

Stellar Suppression Based on the Enhanced Dilation Difference Algorithm
The stellar suppression algorithm based on enhanced dilation difference is an improved algorithm proposed in this paper to overcome the shortcomings of dilation difference algorithm. GSO debris detection generally adopts target tracking mode, so the translation of all stars in CCD images is basically the same, and the motion mode of GSO debris is different from that of stars whether GEO or non-GEO objects, so we can suppress stars by using inter-frame difference technology after image registration. There are two main problems to be solved in the suppression of stars by inter-frame difference technology. Firstly, using adaptive fast image registration algorithm, the accuracy of translation vector in horizontal and vertical directions is 1 pixel. Therefore, if absolute difference is used directly, there will be residual marks on the edge, even if the translation vector obtained is accurate, because there are many stars, there are small differences in their translation, and the difference results will also leave edge remnants. In particular, when GSO observations are performed, the camera usually does not have a mechanical shutter, which results in a serious smear of high-brightness stars, and the edge remnants left by the smear will seriously affect the subsequent target segmentation effect. Secondly, the intensity of stars varies in time because of the interference of various internal and external factors, that is, the intensity of the same star in the first and second frames is generally inconsistent, especially the intensity of high-brightness stars is often very different, which leads to the difference results leaving the remnants of stars. The remnants of high-brightness stars can be severe and can even lead to the inability to segment targets. For the first problem, we can use the dilation difference algorithm to solve [38]. Dilation difference first dilates the subtracted image using structural element b. Considering the accuracy of the registration algorithm, this paper chooses b as a flat structure element of 3 × 3. Let I 0 (x, y) be the intensity of the image at (x, y) after background suppression of the current frame, and I 1 (x, y) be the intensity of the image at (x, y) after background suppression of the previous frame and already registered with the current frame. The dilation difference formula is as follows: where I d_01 (x, y) represents the intensity of the image at (x, y) of the dilation difference result, " ⊕ " is the dilation operator. Dilation difference algorithm solves the problem of stellar edge remnants caused by scale mismatch. For the second problem, we need to solve the problem of stellar remnants caused by the difference of stellar intensity. The specific reason is that the intensity of the star in the current frame is stronger than that in the previous frame. In order to solve this problem, we propose multiplying an intensity enhancement factor η after the dilation operation of the previous image.
It should be pointed out that the same GSO debris will not be in the same position in the previous frame and the current frame image after registration because of the difference between the motion mode of GSO debris and the stars, while the enhancement factor only changes the gray value of the image, so the enhanced dilation difference operation will not suppress the GSO debris. The enhanced dilation difference formula is as follows: where I ed_01 (x, y) represents the intensity of the image at (x, y) after the enhanced dilation difference operation, and the enhancement factor η generally ranges from 1.5 to 2. In order to determine the range of η, the average number of false alarms and missed detections per frame of 496 images in seven image sequences is counted when η = 1 : 0.1 : 3, as shown in Figure 7. Thus, when 1.5 ≤ η ≤ 2, the number of false alarms and missed detection is relatively small. In addition, according to [39][40][41], the magnitude of a variable star will not change more than 0.75 in 5 to 15 seconds of exposure time, that is, the brightness of the same variable star is less than twice that of the previous frame. In order to show the effect of the algorithm, the result image and the corresponding gray distribution of the same high-brightness star after the operation of absolute difference, dilation difference and enhanced dilation difference are given, as shown in Figure 8.

GSO Debris Enhancement and Target Segmentation Based on Inter-Frame Correlation and Threshold Segmentation Technology
This section contains two processing steps, GSO debris enhancement and target segmentation, but because the two operations are closely related, we did not discuss them separately for ease of understanding. In view of the weak SNR of GSO debris on CCD image, there are two main methods to improve the SNR of the target. The first method is to prolong the exposure time. This strategy is indeed used in the GSO debris observation task. However, it should be noted that exposure time is limited. When the exposure time is too long, it will cause the saturation overflow of CCD sensor pixels. In addition, the long exposure time only improves the SNR of the GEO debris, but has little effect on the non-GEO debris. The second method is to accumulate energy by using multiple frames. This paper uses the inter-frame correlation processing algorithm to achieve this goal. Firstly, the current frame and the previous frame, and the current frame and the first two frames are used to obtain two differential images Then, multiplying the two differential images can enhance the GSO debris and further suppress the stars:

GSO Debris Enhancement and Target Segmentation Based on Inter-Frame Correlation and Threshold Segmentation Technology
This section contains two processing steps, GSO debris enhancement and target segmentation, but because the two operations are closely related, we did not discuss them separately for ease of understanding. In view of the weak SNR of GSO debris on CCD image, there are two main methods to improve the SNR of the target. The first method is to prolong the exposure time. This strategy is indeed used in the GSO debris observation task. However, it should be noted that exposure time is limited. When the exposure time is too long, it will cause the saturation overflow of CCD sensor pixels. In addition, the long exposure time only improves the SNR of the GEO debris, but has little effect on the non-GEO debris. The second method is to accumulate energy by using multiple frames. This paper uses the inter-frame correlation processing algorithm to achieve this goal. Firstly, the current frame and the previous frame, and the current frame and the first two frames are used to obtain two differential images I ed_01 and I ed_02 respectively by using the enhanced dilation difference algorithm.
Then, multiplying the two differential images can enhance the GSO debris and further suppress the stars: I multip_0 (x, y) = I ed_01 (x, y) * I ed_02 (x, y), where I multip_0 represents the result of the current frame processed by inter-frame correlation. At this point, we get the background and stars are effectively suppressed, the GSO debris are effectively enhanced image, then we can extract the suspected GSO debris on I multip_0 .
In this paper, threshold segmentation technology is used to extract GSO debris. Threshold segmentation is mainly divided into global threshold segmentation and local threshold segmentation. Generally speaking, the performance of segmentation algorithm based on local threshold is better than that of global threshold, but the complexity of the algorithm is much higher than that of global threshold. In view of the fact that the background and stars in the CCD star image have been effectively suppressed based on the preprocessing algorithm mentioned above, in order to improve the processing speed, we choose to segment and extract GSO debris using global thresholds. Firstly, the segmentation threshold is calculated as follows: where, E I multip_0 represents the mean of the image, SD I multip_0 represents the standard deviation of the image, and κ is the adjustment coefficient, which is determined by the minimum SNR of the detectable target. For dim and small GSO debris segmentation, in order to improve detection probability, κ is usually set as small as possible, but at the same time, it also increases the false alarm probability when extracting suspected targets, which challenges how to efficiently complete GSO debris track association. After the segmentation threshold is determined, the image is segmented by the threshold, and the binary image I binary_0 is obtained as follows: For the binary image I binary_0 , the connected area is searched according to the 8-connected objects, and each connected component is labeled. Finally, the centroid of each suspected target area is calculated. Generally, there are three kinds of suspected targets: GEO debris, non-GEO debris and false targets. The image of an GEO debris is a circular dim region. The image of a non-GEO debris appears as a dim strip or ellipse. False targets are usually incompletely suppressed stars, and their images are long strips, but their lengths or directions are generally different from those of non-GEO debris. There are two methods for calculating the centroid of an object image: fitting method and modified moment method. Fitting method has high positioning accuracy, but its calculation is huge, and its accuracy depends heavily on the accuracy of the fitting template. The shape of the extracted suspected object image is different, so it is difficult to determine the fitting template. Modified moment method has good accuracy for centroid location of dim target and low computational complexity. Considering the complexity and accuracy of the algorithm, the modified moment method is used to locate the centroid of the suspected target. The calculation formula is as follows: where (centroid_x, centroid_y) represents the centroid coordinate of the target, D is the target area, that is, the connected area of the suspected target, and I(x, y) represents the image used to calculate the centroid of the target. In order to avoid background interference, this paper calculates the centroid of the target using the image after background suppression. When the centroid of all suspected targets is extracted, all steps of GSO debris extractor are completed. It should be noted that we only extract the centroid information of the suspected target, but not the size, gray level, shape and other information, which is often an important basis for judging whether the target is a false alarm. This is because through the measured data testing, we found that the GSO debris extractor used in this paper can suppress background and stars very well. Among all the suspected targets extracted from a single frame image, the false targets are generally about 2 to 3, and in extreme cases, there will be no more than 10, so there is no need to eliminate false targets. In addition, because GSO debris are dim and small targets, if false alarm removal is performed, the GSO debris is likely to be deleted by mistake. Although we do not extract other information of suspected targets, the algorithm has good expansibility, and it is easy to extract other information in parallel while calculating the centroid.

GSO Debris Tracker
Track Association of GSO debris belongs to the field of multi-target tracking. There are mainly two kinds of multi-target tracking algorithms, one is based on RFS and the other is based on data association. Multi-target tracker based on RFS includes unlabeled RFS type and labeled RFS type. Unlabeled RFS type has the advantage of fast computation, but it does not have the ability to distinguish target tracks. Strictly speaking, it can be considered that the unlabeled label type is a target validator, and subsequent processing steps are needed to obtain target tracks. Labeled RFS type has the ability to distinguish target tracks, but it has a large amount of computation. In addition, all RFS trackers require that the target's birth model be known prior, while in GSO debris tracking applications, the birth model is unknown. Multi-target tracking algorithm based on data association includes GNN, JPDA, MHT and their variants. In complex application scenarios, the tracking effect of MHT and JPDA is better than that of GNN, but the real-time performance of GNN is better than that of MHT and JPDA. Considering the simplicity of GSO debris motion, approximately stationary or uniformly low-speed linear motion, and the GSO debris extractor used in this paper can suppress background and stars very well, and the extraction results contain few false alarms. In this case, GNN multi-target tracker can complete GSO debris tracking task very well. Therefore, this paper uses GNN multi-target tracker to track GSO debris.

Adaptive Inter-Frame Interval Estimation
Before introducing the principle of GNN multi-target tracker, there is still a problem to be solved. Any tracker needs a relatively accurate time input, that is, the inter-frame interval is known. It is not necessary to consider this problem when the image sequence is of equal interval, but sometimes the exposure time is changed during GSO debris observation, which results in inconsistent inter-frame intervals. When the tracking error caused by time error is larger than the maximum error allowed by the tracker, it will cause the failure of target-track association at that time. The whole image sequence shows that the track is discontinuous, that is, the same target has multiple tracks. Therefore, it is necessary to extract the information related to the inter-frame interval from the image. Through analysis, the inter-frame interval can be represented by the translation vector obtained during image registration. Specifically, the inter-frame interval is the sum of charge transfer time and exposure time. The charge transfer time is much smaller than the exposure time, so the inter-frame interval can be approximated to the exposure time. When the telescope is in target tracking mode, the stars will form stripes on the CCD image. The stripes are caused by the rotation of the earth during exposure time, and the rotation speed of the earth is basically constant. Therefore, the exposure time can be estimated by the length of the stripes. The length of the stripes is positively correlated with the modulus of the translation vector, so the modulus of the translation vector can represent the exposure time. In summary, the modulus of the translation vector can indirectly represent the inter-frame interval.
In this paper, the translation vector modulus of 5 seconds exposure is taken as the reference modulus. The ratio of translation vector modulus obtained by adaptive fast registration algorithm to reference modulus is used as an adaptive estimation of inter-frame interval.

GNN Multi-Target Tracking Algorithm
This section introduces the principle of GNN multi-target tracker and gives a comparison with other trackers in GSO debris tracking applications. In the following description, X o = x o_1 , · · · , x o_m is used to represent the track table of the previous moment, where m is the number of targets. X = x 1 , · · · , x m represents the predicted track table at the current time, P 1 , · · · , P m represents the corresponding prediction error covariance matrix. Linear multiple-target motion model X = FX o + V = X + V is adopted, where X = x 1 , · · · , x q is the real track table at the current time, F is the state transition matrix, V is the process noise with zero mean and covariance matrix Q; then the transfer function can be expressed as f x |x o = N Q x − Fx o . Z = {z 1 , · · · , z n } is used to represent the suspected target state table extracted at the current time, that is, the observation at the current time. The observation model is Z = H X + W, where H is the observation matrix, W is the observation noise with zero mean and covariance matrix R; then the likelihood function can be expressed as f z x = N R z − H x .
Considering the Euclidean distance between observation z and prediction x , the uncertainty R of z and the uncertainty P of x , referring to [42] (p. 95), the associated distance d z, x between observation and track is defined by Bayesian method as follows: GNN multi-target tracker is based on global associated distance. All possible associated distance between observation and prediction tracks are considered as follows: where π : {1, · · · , m} → {1, · · · , n} is any possible one-to-one mapping. Therefore, the optimal association problem can be solved by solving π, which minimizes the following formula: The mapping π obtained here is the optimal association between observation and prediction tracks determined by GNN multi-target tracking algorithm.
GNN multi-target tracking algorithm has low complexity and good tracking effect when the motion model is relatively simple. Let Dnum, Tnum and Fnum represent the number of observations, tracks and false alarms respectively, and assume that there is no target missed detection. Then the complexity of GNN tracking algorithm is O C Tnum Dnum Tnum! = O(Dnum!/Fnum!). Since the motion form of GSO debris is relatively simple, GNN algorithm is adopted in this paper. In order to show the rationality of using GNN tracker in GSO debris tracking application, we record the processing time of GSO debris tracking using GNN tracker and other trackers in seconds, such as Table 1. It should be noted that according to the processing results of the measured data, all the tracking algorithms used in the table can achieve better tracking results. There is no missing track in all tracking algorithms. There are three and one false tracks in the tracking results of PHD and CPHD algorithms, respectively. There are no false tracks in the tracking results of other algorithms. When the target track is established, the tracking algorithm based on RFS and MHT has a delay of 3-6 frames, and the GNN and JPDA algorithms have a delay of 3 frames. Seq_I-Seq_VII image sequences were tested. The numbers in parentheses represent the number of images in the corresponding image sequences. JointGLMB and GLMB are tracking algorithms based on labeled RFS. The contents in parentheses indicate the implementation method adopted. SMC represents sequential Monte Carlo method and GMS represents Gauss mixture method. CBMeMBer, CPHD and PHD are tracking algorithms based on unlabeled RFS, and their processing results need follow-up processing. In order to complete the test of the RFS tracker, we first analyze the data to obtain the birth model, and then use it as the input of the RFS tracker. This is obviously unreasonable in practical application. However, it should be noted that the unlabeled RFS tracker has a better processing speed advantage, and it is necessary to further study it, but this is beyond the content of this paper. MHT, JPDA and GNN are all tracking algorithms based on data association. When tracking applications are complex, MHT has the best performance, followed by JPDA and GNN. However, because the motion mode of target in GSO debris tracking application is relatively simple, and the target extractor used in this paper can suppress background and stars very well, there are a few false alarms in the extracted results. Therefore, there is no significant difference in the effectiveness of the three trackers. According to the above information, considering the tracking effect and real-time performance of the algorithm, it is reasonable to select GNN multi-target tracker for GSO debris tracking.

Introduction of Measured Data
In order to test the algorithm comprehensively, we collected seven groups of image sequences by two telescopes, named Seq_I-Seq_VII. Seq_I-Seq_III was collected by telescope I, and Seq_IV-Seq_VII was collected by telescope II. The observation sites of telescope I and II are different, and their parameters are identical. The number of frames in each sequence is shown in Table 1. The exposure time of a single frame is usually between 5 and 10 seconds, so the frame rate is 0.1 to 0.2 frames per second. The images are all 16-bit gray-scale images of 2K × 2K, and the single frame data is 8MB. We annotated ground-truth labels based on the public version of USSTRATCOM catalogue [44] and combined with manual screening. The characteristics of each image sequence are described as follows: Seq_I, which contains two frames of irrelevant images, is located in the first frame and in the middle of the sequence, to test the adaptability of the registration algorithm in case of abnormality. It also contains two targets with very close distance. As shown in Figure 9a, the distance between two targets is less than 10 pixels, which can test the algorithm's ability to distinguish close targets.
Seq_VII is a typical data of telescope II, which can be used to test the performance of the algorithm under normal conditions.

Processing Results of Measured Data
All the tests in this paper are completed on Lenovo M410-PC with Intel i7-7700CPU@3.6GHz (8 CPUs), 8GB DDR4 RAM by using the software of MATLAB R2019(a). In order to show the detection ability for dim and small targets, we give the time-varying SNR of all targets detected in Seq_III, as shown in Figure 10, and the SNR of targets after image preprocessing, as shown in Figure 11. Seq_II, which contains more targets, is used to test the performance of the algorithm in the case of a large number of targets. Also, a new non-GEO debris enters the field of view, which can test the acquisition ability of the algorithm for the new-emerging target.
Seq_III is the typical data of Telescope I. It contains a very dim GEO debris, which can be used to test the limit detection ability of the algorithm. As shown in Figure 9b, the SNR of the target is very small. It also contains a non-GEO debris whose direction of motion is basically the same as that of stars, but whose velocity is smaller than that of stars. As shown in Figure 9c, the stripe direction of the target is basically the same as that of stars, but the stripe length is shorter than that of stars, which can be used to test the ability of the algorithm to detect targets with similar stellar motions.
Seq_IV image quality is poor. There are multiple frame dropouts and interference from fast moving bright object. As shown in Figure 9d, the interference target has a long bright smear, which can be used to test the performance of the algorithm in complex situations.
Seq_V, there is a change of field of view, testing the performance of the algorithm when the field of view changes. In addition, the exposure time changed once. As shown in Figure 9e, due to the change of exposure time, the length of the stellar stripe in frame 16 is shorter than that in frame 17, which can be used to test the adaptive ability of the algorithm when the exposure time changes.
Seq_VI, in which there is strong moonlight interference. As shown in Figure 9f, moonlight causes strong background noise, which can be used to test the effectiveness of background suppression algorithm.
Seq_VII is a typical data of telescope II, which can be used to test the performance of the algorithm under normal conditions.

Processing Results of Measured Data
All the tests in this paper are completed on Lenovo M410-PC with Intel i7-7700CPU@3.6GHz (8 CPUs), 8GB DDR4 RAM by using the software of MATLAB R2019(a). In order to show the detection ability for dim and small targets, we give the time-varying SNR of all targets detected in Seq_III, as shown in Figure 10, and the SNR of targets after image preprocessing, as shown in Figure 11.  Figure 10. SNR of targets before image preprocessing. Figure 10. SNR of targets before image preprocessing.  Figure 11. SNR of targets after image preprocessing. Figures 10 and 11 show that the proposed image processing algorithm can effectively improve the SNR of the target, and then increase the detection probability of dim and small targets. Figure 12 shows the receiver operating characteristic curves (ROC) of the proposed algorithm and the twice frame differential algorithm [1]. The method of drawing the ROC curve is referred to [45]. Figure 11. SNR of targets after image preprocessing. Figures 10 and 11 show that the proposed image processing algorithm can effectively improve the SNR of the target, and then increase the detection probability of dim and small targets. Figure 12 shows the receiver operating characteristic curves (ROC) of the proposed algorithm and the twice frame differential algorithm [1]. The method of drawing the ROC curve is referred to [45].  Figure 12 shows that the proposed algorithm has better detection performance than the twice frame differential algorithm. Table 2 shows the number of GEO debris and non-GEO debris detected in seven groups of sequences, and the mean of SNR enhancement multiples of all targets in the sequence after preprocessing. Detection accuracy Figure 12. The receiver operating characteristic curves (ROC) of the proposed algorithm and the twice frame differential algorithm. Figure 12 shows that the proposed algorithm has better detection performance than the twice frame differential algorithm. Table 2 shows the number of GEO debris and non-GEO debris detected in seven groups of sequences, and the mean of SNR enhancement multiples of all targets in the sequence after preprocessing.  Table 2 shows that the algorithm can detect not only GEO debris, but also non-GEO debris, and at least can improve the SNR of detected targets by about 200 multiples. Table 3 gives the minimum SNR and missed detection accuracy of all detected GSO debris in seven sequences, in which SNR_min represents the minimum SNR of detected targets and Omi_R represents the missed detection accuracy. The definition of Omi_R in this paper is as follows: where M_num denotes the number of frames that the target has not been detected, E_num denotes the total number of frames in which the target exists.  Table 3 shows that the minimum SNR of the target detected by this algorithm is not more than 5.91 dB. At the same time, it should be noted that the missed detection accuracy of some targets is very high, so we relax the condition of judging the target track loss in the tracking algorithm. Only when the target is missed eight times in a row can it be judged that the track is lost. In addition, when the minimum SNR of different detectable targets is almost the same, the target miss detection accuracy is significantly different. This is due to the time-varying brightness of GSO objects. Because of the rotation and irregular shape of the object, the brightness of the object changes with time, and the change rate will be very large [21,22]. The intense change of brightness results in discontinuous detection results of targets, which leads to the above problems. For example, let the minimum SNR of the target detected by the algorithm be 5.9 dB. At k-time, the SNR of target A and target B are both 7 dB and detected. At k + 1, the target A changes bright, the SNR increases, the target B changes dark, and the SNR is less than 5.9 dB, which results in that A can be detected and B cannot be detected. That is to say, the minimum detectable SNR of A and B is 7 dB, but the missed detection accuracy is different. The number of false targets detected in each frame of seven image sequences is shown in Figure 13.  Table 3 shows that the minimum SNR of the target detected by this algorithm is not more than 5.91 dB. At the same time, it should be noted that the missed detection accuracy of some targets is very high, so we relax the condition of judging the target track loss in the tracking algorithm. Only when the target is missed eight times in a row can it be judged that the track is lost. In addition, when the minimum SNR of different detectable targets is almost the same, the target miss detection accuracy is significantly different. This is due to the time-varying brightness of GSO objects. Because of the rotation and irregular shape of the object, the brightness of the object changes with time, and the change rate will be very large [21,22]. The intense change of brightness results in discontinuous detection results of targets, which leads to the above problems. For example, let the minimum SNR of the target detected by the algorithm be 5.9 dB. At k-time, the SNR of target A and target B are both 7 dB and detected. At k + 1, the target A changes bright, the SNR increases, the target B changes dark, and the SNR is less than 5.9 dB, which results in that A can be detected and B cannot be detected. That is to say, the minimum detectable SNR of A and B is 7 dB, but the missed detection accuracy is different. The number of false targets detected in each frame of seven image sequences is shown in Figure 13.  Figure 13 shows that the GSO debris extractor can extract no more than 10 false targets per frame, generally less than 4. This is because the extractor can suppress background and stars very well, and effectively enhance GSO debris. The processing time of each sequence is shown in Table 4.  Table 4 shows that the single frame processing time of the proposed algorithm is not more than 2.112 s, while the exposure time of GSO debris observed with telescopes is generally more than 5 seconds, so the algorithm designed in this paper fully meets the real-time requirements of GSO debris detection tasks. In addition, there are differences in the processing time of single frame in Table 4, Figure 13. The number of false targets detected in each frame of seven image sequences. Figure 13 shows that the GSO debris extractor can extract no more than 10 false targets per frame, generally less than 4. This is because the extractor can suppress background and stars very well, and effectively enhance GSO debris. The processing time of each sequence is shown in Table 4.  Table 4 shows that the single frame processing time of the proposed algorithm is not more than 2.112 s, while the exposure time of GSO debris observed with telescopes is generally more than 5 seconds, so the algorithm designed in this paper fully meets the real-time requirements of GSO debris detection tasks. In addition, there are differences in the processing time of single frame in Table 4, which is caused by the different stability of exposure time in image sequence. If the exposure time is stable, there is no need to use large-scale search frequently for registration, and the processing speed will be fast. However, if the exposure time changes frequently, it will lead to frequent use of large-scale search for registration, resulting in increased processing time. Finally, a representative Seq_III detection result image is given, as shown in Figure 14. Among them, color lines represent tracks, and different tracks are represented by different colors. The target with obvious track is non-GEO debris, otherwise it is GEO debris. The green rectangle represents the position of the target whose track has been established, and the red dot represents the position of the suspected target. Obviously, in the detection result of Seq_III, target 1, 2, 3, 4, 6, 9 is GEO debris and 5, 7, 8 is non-GEO debris. scale search for registration, resulting in increased processing time. Finally, a representative Seq_III detection result image is given, as shown in Figure 14. Among them, color lines represent tracks, and different tracks are represented by different colors. The target with obvious track is non-GEO debris, otherwise it is GEO debris. The green rectangle represents the position of the target whose track has been established, and the red dot represents the position of the suspected target. Obviously, in the detection result of Seq_III, target 1, 2, 3, 4, 6, 9 is GEO debris and 5, 7, 8 is non-GEO debris.

Conclusions and Prospects
In this paper, an adaptive real-time detection algorithm for dim and small photoelectric GSO debris was proposed. It can effectively improve the detection ability of dim and small GSO debris. It has strong adaptability and does not require exposure time, telescope position and optical axis pointing information. Moreover, the real-time performance of the algorithm is strong, and it can fully meet the requirements of real-time processing ability in the application of GSO debris photoelectric detection. In addition, the algorithm has good parallelism. The algorithm consists of GSO debris extractor and GSO debris tracker. In the process of GSO debris extraction, an adaptive fast registration algorithm is proposed. The adaptability and real-time performance of the algorithm are effectively improved. In addition, in view of the shortcoming of the dilation difference algorithm, an enhanced dilation difference algorithm is proposed to improve the suppression effect for stars. For GSO debris tracker, an adaptive inter-frame interval technique is proposed to improve the performance. Using seven groups of measured data in different situations, the algorithm is tested comprehensively. The test results show that the proposed algorithm has the ability of self-adaptive detection of dim and small photoelectric GSO debris, and its processing speed fully meets the realtime requirement.
Some prospects for future work can be given. Firstly, because the algorithm has good parallel processing characteristics, it can be transplanted to embedded systems to further improve the processing speed. Secondly, in view of the theoretical advantages of RFS tracker, it is meaningful to find the method to solve the problem of unknown birth model in GSO debris tracking application with RFS tracker. Thirdly, considering that ISAL can image GSO debris, it can combine GSO optical detection system with ISAL imaging system to improve the recognition of GSO debris.

Conclusions and Prospects
In this paper, an adaptive real-time detection algorithm for dim and small photoelectric GSO debris was proposed. It can effectively improve the detection ability of dim and small GSO debris. It has strong adaptability and does not require exposure time, telescope position and optical axis pointing information. Moreover, the real-time performance of the algorithm is strong, and it can fully meet the requirements of real-time processing ability in the application of GSO debris photoelectric detection. In addition, the algorithm has good parallelism. The algorithm consists of GSO debris extractor and GSO debris tracker. In the process of GSO debris extraction, an adaptive fast registration algorithm is proposed. The adaptability and real-time performance of the algorithm are effectively improved. In addition, in view of the shortcoming of the dilation difference algorithm, an enhanced dilation difference algorithm is proposed to improve the suppression effect for stars. For GSO debris tracker, an adaptive inter-frame interval technique is proposed to improve the performance. Using seven groups of measured data in different situations, the algorithm is tested comprehensively. The test results show that the proposed algorithm has the ability of self-adaptive detection of dim and small photoelectric GSO debris, and its processing speed fully meets the real-time requirement.
Some prospects for future work can be given. Firstly, because the algorithm has good parallel processing characteristics, it can be transplanted to embedded systems to further improve the processing speed. Secondly, in view of the theoretical advantages of RFS tracker, it is meaningful to find the method to solve the problem of unknown birth model in GSO debris tracking application with RFS tracker. Thirdly, considering that ISAL can image GSO debris, it can combine GSO optical detection system with ISAL imaging system to improve the recognition of GSO debris.