Space Debris Detection and Positioning Technology Based on Multiple Star Trackers

Abstract

This paper focuses on the opportunity to use multiple star trackers to help space situational awareness and space surveillance. Catalogs of space debris around Earth are usually based on ground-based measurements, which rely on data provided by ground-based radar observations and ground-based optical observations. However, space-based observations offer new opportunities because they are independent of the weather and the circadian rhythms to which the ground system is subjected. Consequently, space-based optical observations improve the possibility of space debris detection and cataloging. This work deals with a feasibility study of an innovative strategy, which consists of the use of a star sensor with a dedicated algorithm to run directly on board. This approach minimizes the impact on the original mission of the satellite, and on this basis, it has also the function of space debris monitoring. Therefore, theoretically, every satellite with a star tracker can be used as a space surveillance observer. In this paper, we propose a multi-star space debris detecting and positioning method with constant geocentric observation. Using the multi-star tracker joint positioning method, the angle measurement data of the star tracker is converted into the spatial coordinates of the target. In addition, the Gaussian MMSE difference correction algorithm is used to realize the target positioning of multiple optical observations, and the spatial target position information of the multi-frame images is fused, thus completing the solution of the orbit equation. The simulation results show that the proposed method is sufficient to detect and position space debris. It also demonstrates the necessity and feasibility of cooperative network observation by multiple star trackers.

1. Introduction

With the increasingly active space launch activities of various countries, the number of spacecraft worldwide is on a surge. The ensuing mass of space debris threatens the operational safety of orbiting spacecraft [1]. A complete catalog of space debris is useful for current and future missions, as unexpected impacts may compromise the mission of the impacted operative platforms [2]. Given the necessity of many interactions among international stakeholders [3], cooperation, responsibility, liability, and the obligations to register space objects under the existing space treaties are fundamental concerns of space situational awareness (SSA) [4].

At present, space debris is mainly detected by ground-based observation and then cataloged, recorded, and entered into a database by major aerospace nations. There are two types of ground-based observations: ground-based radar observations and ground-based optical observations. More than 50 radars are currently used worldwide for space target and space debris monitoring. More and more efforts have been made to extend the network of available radar systems devoted to the control of space [5]. Ground-based radar detection methods are not limited by factors such as size and weight. A radar with a large aperture antenna and large transmitting power is usually chosen to obtain high detection accuracy and increase the detection distance. However, since the target signal loss is proportional to the fourth power of distance, the detection range of the radar is usually limited to low orbit. In addition, due to the effects of atmospheric transmission jitter, ionospheric scintillation, astronomical refraction error, and atmospheric attenuation limitations, the radar can only operate in lower frequency bands, which will limit the detection accuracy of ground-based radar. At present, there are still great difficulties for ground-based radar to identify small-sized space debris [6].

Compared with ground-based radar detection, ground-based optical detection can obtain a longer observation distance because the signal attenuation is proportional to the square of the distance, thus enabling the observation of large-size space debris in higher orbit [7]. With regard to the ground-based optical observations, ESA has installed a 1 m Schmidt telescope on the Canary Islands with a limiting detection size of about 15 cm in order to characterize the GEO environment [8]. Simulated results are reported, where it is claimed that the accuracy of the satellite’s estimated position is better than 5 arc seconds, and the satellite’s tracking accuracy is around 1–3 arc seconds. Other ongoing missions focused on ground-based optical observations are described in [9]. However, due to the limitation of lighting conditions and other factors, the observation efficiency of ground-based optical detection methods is low, and clouds, fog, atmospheric pollution, urban glow or glow at full moon, and the thermal effect of the atmosphere will increase the background noise, thus reducing the detection signal-to-noise ratio. The atmospheric turbulence also limits the observation accuracy and limiting resolution, thus limiting the detection capability of ground-based optical observation for small-sized space debris.

In recent years, on the basis of research on ground-based monitoring means, some countries have further developed space-based detection means to improve space debris monitoring and warning capabilities, especially the observation and warning capabilities for small-sized space debris with wide distribution, large numbers, and high collision threats [10]. For example, the United States, Canada, and Europe have successively launched plans to establish a space-based space target surveillance system [11].

A pathfinder satellite with the Space Based Visible (SBV) on-board was launched in 2010 in a sun-synchronous orbit to successfully detect objects as small as 1 m³ in the geosynchronous belt. The pathfinder satellite and the geosynchronous space situational awareness program are part of the U.S. Air Force space-based surveillance system program [12]. In Europe, ESA supported and is still supporting several research programs for space-based SSA, such as the space-based space surveillance service [13], the optical In-situ Monitor project [14], A summary of ESA activities related to space-based optical observations of space debris and NEOs is reported in [15].

All of the above missions are designed with dedicated payloads, but most satellites already have one or more-star trackers for attitude determination, which, if properly utilized, can constitute a zero-cost space surveillance system. Using the star as the reference frame for attitude measurement, the star trackers can output the vector direction of the star in the coordinates of the star trackers, thus providing highly accurate measurement data for spacecraft attitude control and astronomical navigation. The star map taken by the star trackers contains not only stars but also space targets, such as space debris [16]. In the starry background, space debris can also reflect part of the energy from the Sun and has similar luminous characteristics to the stars. When the signal-to-noise ratio is high enough, space debris can be captured by star trackers, and all its information is contained in the point target of the star map. Combined with the principle of extracting star points and star map recognition of star trackers [17], it can distinguish space targets from stars. In simple terms, the identification of all the stars in the star map is accomplished with the help of star trackers, and the targets are naturally filtered out. In fact, in multi-star tracker missions, one or more sensors are usually not used to determine the attitude. Therefore, the images acquired in these sensors can be used for surveillance operations. In this way, unlike most space-based missions designed specifically for space situational awareness, no specific attitude maneuver is required [18].

The proposed strategy has also been discussed in previous literature. In [19], a convolutional filter-based centroid estimation image preprocessing was first used. Then, stars were removed using the pyramidal star identification technique, and finally, the spatial target was tracked by combining the center-of-mass estimation and the Kalman filter. However, this work did not address the orbiting of space targets. The initial possibility of using a star tracker for space target detection was discussed in [20]. The idea is to use convolutional neural networks, which is interesting but not feasible with the current technology of spacecraft applications. In [21], a faint streak detection algorithm based on a streak-like spatial filter was proposed, but it took 2–3 h for 4000 × 4000-pixel images. Considering the limited computing power of spacecraft on-board computers, none of these techniques is designed for efficient implementation. Therefore, research on novel observation methods to improve the cost-effectiveness of observation is the main research direction at present.

In this paper, the configuration of a space debris monitoring network based on multiple star trackers is proposed and studied for detection feasibility. Unlike the previous algorithms, the proposed algorithm uses satellite equipped star trackers to detect and extract space targets from star maps taken during star tracker attitude fixation and mines the available observation information of space targets from a large amount of data to achieve space target surveillance, thus saving the cost of launching dedicated surveillance satellites. The space target is observed by a single satellite star tracker and the target orientation data are given. The maximum projection time information method is used to achieve target detection, and then the goniometric data of a single optical platform are converted into the spatial coordinates of the target by a Gaussian minimum mean square error difference correction algorithm through a space-based multi-optical platform to achieve target localization. Subsequently, the multi-frame image spatial target position information is fused to complete the solution of orbital equations, followed by space debris orbit determination for accurate safety prediction of spacecraft in-orbit operation.

By applying the performance advantages of star trackers to space target detection, the advantages of star trackers in terms of high accuracy and high update rate are thus effectively exploited in the field of space detection. Different from the tracking and observation mode of specialized monitoring satellites, the space debris monitoring network based on multiple star trackers proposed in this paper relies on the natural intersection between the field of view of star trackers and the target to capture the target. Although each continuous observation is short, the number of star trackers is large, which can realize the intermittent short arc observation of space targets under different observation geometry. This kind of network does not require enhancing the detection ability of the system but increases the detection opportunities of the system so that more targets can be detected by the system at the same time. In fact, it improves the visible turn state between the target and the platform and provides ideas and engineering practice for the expansion and application of star trackers.

This paper is structured as follows: Section 2 presents the comprehensive analysis of space debris. Section 3 presents the star tracker model. Section 4 focuses on the target identification method between star trackers. Section 5 presents the localization and initial orbiting of space targets, presenting both typical and advanced operations required for the specific task of this work. Section 6 reports the results obtained from the feasibility and characterization experiments, the results of the simulation validation, and a discussion of the related computational effort. Finally, concluding remarks are given in Section 7.

2. Comprehensive Analysis of Space Debris

2.1. Observability Analysis of Space Debris

The observation of space debris by the star tracker should consider the influence of the Earth and the Sun: when the space debris is blocked by the Earth, the satellite cannot detect the space debris; when the space debris is within the Earth’s shadow region, the satellite cannot observe it either; when the satellite observes the debris against the Sun’s light, the observation cannot be achieved. The analysis of the visible situation can be divided into the following three cases:

(1)

Consider the effect of ground shadow

Space debris is captured by the CCD camera of the tracking satellite only when it is illuminated by the Sun. As space debris moves around the Earth, the Earth is irradiated by daylight and creates a ground shadow (as shown in Figure 1). Whether the space debris is within the ground shadow can be expressed by the irradiation factor F [22], which is defined as

$$F=\{\begin{array}{c}1\\ 0\end{array}\begin{array}{c}\\ \end{array}\begin{array}{c}Spacedebrisoutsidethegroundshadow\\ Spacedebriswithinthegroundshadow\end{array}$$

Assuming that the sunlight is parallel (${\overrightarrow{R}}_{SOL}$ is the sunlight unit vector), i.e., the ground shadow is a cylindrical ground shadow model, ${\overrightarrow{r}}_D$ is the vector of space debris relative to the center of the earth, we can obtain D and H

$$D={\overrightarrow{r}}_D\cdot {\overrightarrow{R}}_{SOL}$$

(1)

$$H^2={\left|{\overrightarrow{r}}_D\right|}^2-D^2$$

(2)

When $D\le 0$, the space debris is illuminated by sunlight and $F=1$.

When $D>0$, if $H^2\ge R^2$, the space debris is in the sunlight, F = 1; if $H^2<R^2$, the space debris is in the earth’s shadow, $F=0$.

(2)

Considering the effect along and against the sunlight

When the satellite is tracked against the Sun to observe the debris, the background light is so strong that the observation results will be blurred and invalid. Considering such factors, we define the tracking sun angle θ as the angle between the tracking satellite and the line connecting the space debris and the sun (as shown in Figure 2), and according to the performance index of the CCD camera, θ should be greater than a certain value (e.g., 80°) for effective observation.

(3)

Considering the effect of earth–atmosphere radiation

Since weak energy is available for the star point and space debris acquired by the star tracker, the star tracker cannot work properly when there are too many stray lights in the field of view after the local gas light enters the optical system, so the space debris detection based on the star tracker can only function from low orbit to high orbit or in the same orbit direction, as shown in Figure 3.

2.2. Analysis of Space Debris Target Characteristics

Targets in star maps usually exhibit point-like characteristics without structure and texture features. Star trackers detect the sunlight reflected from the space target, and the brightness of the resulting space target image varies to some extent due to the different reflectivity of different parts of the space target and the different observation phase angles at different observation moments. In general, the factors affecting the brightness of the detected target include the size of the target, the distance between the detector and the target, the target material properties, and the observation azimuth angle.

According to theoretical analysis and experimental verification, the calculation formula of spatial target brightness characteristics is [23,24]:

$$M_V=-26.7-2.5\log (A\rho F(\phi ))+5.0\log (R)$$

(3)

where A represents the lateral area of the target, $\rho$ denotes the reflection coefficient of the target, R denotes the distance between the target and the observer, $\phi$ represents the sun–target–observer phase angle, and F is the phase function, which is related to the current situation and direction of the target.

Assuming the target reflectance $\rho =0.1$ and the observed dislocation angle $\phi =0$, the calculation can be performed according to the typical spherical target brightness characteristics analysis shown in

Table 1. Analysis of brightness characteristics of typical spherical targets.

Target Diameter D	Diffuse Target Distance R = 200 km	Diffuse Target Distance R = 100 km	Diffuse Target Distance R = 50 km
1 m	+3.8 Mv	+2.3 Mv	+0.8 Mv
0.5 m	+5.3 Mv	+3.81 Mv	+2.31 Mv
0.3 m	+6.63 Mv	+4.92 Mv	+3.42 Mv
0.2 m	+7.3 Mv	+5.8 Mv	+4.3 Mv
0.1 m	+8.8 Mv	+7.3 Mv	+5.8 Mv
0.05 m	+9.5 Mv	+8.0 Mv	+6.5 Mv

.

The observation limit magnitude of a traditional star tracker is 5~6 magnitudes. Combining with the above table, we find that the star tracker can complete the observation of a target with a diameter over 0.5 m at 200 km, a target with a diameter over 0.2 m at 100 km, and a target with a diameter over 0.1 m at 50 km. In summary, it is clear that star trackers have the ability to observe space debris near its own orbit.

3. Star Tracker Observation Model

In the process of observing space debris, assuming that the vector of debris observed by the star tracker in the geocentric coordinate system is $\overrightarrow{\rho }$, the basic relationship between the instantaneous geocentric position vector of debris ${\overrightarrow{r}}_{deb}$ and the position vector of the satellite in which the star tracker is located is as follows [25]:

$$\overrightarrow{\rho }={\overrightarrow{r}}_{deb}-{\overrightarrow{r}}_{sat}=\lfloor \begin{array}{c}{x}_{deb}-x_{sat}\\ y_{deb}-y_{sat}\\ z_{deb}-z_{sat}\end{array}\rfloor$$

(4)

where ${\overrightarrow{r}}_{deb}={\left[\begin{array}{ccc}{x}_{deb}& y_{deb}& z_{deb}\end{array}\right]}^T$ and ${\overrightarrow{r}}_{sat}={\left[\begin{array}{ccc}{x}_{sat}& y_{sat}& z_{sat}\end{array}\right]}^T$ are the coordinates in the inertial system of the debris and star tracker, respectively, and the spatial polar coordinates of the observed debris $\overrightarrow{\rho }$ are expressed as follows:

$$\overrightarrow{\rho }=\rho \left[\begin{array}{c}\cos {\delta }_{star}\cos {\alpha }_{star}\\ \cos {\delta }_{star}\sin {\alpha }_{star}\\ \sin {\delta }_{star}\end{array}\right]$$

(5)

where $\rho$ is the distance from the debris to the observation star tracker, $\alpha$ and $\delta$ the rectascension and declination of the debris. The angle subscript “star” indicates the observation vector of the rectascension and declination by the star tracker.

$${\overrightarrow{r}}_{deb}=r\left[\begin{array}{c}\cos {\delta }_{earth}\cos {\alpha }_{earth}\\ \cos {\delta }_{earth}\sin {\alpha }_{earth}\\ \sin {\delta }_{earth}\end{array}\right]$$

(6)

The angle subscript “earth” indicates the rectascension and declination of the debris vector observed under the center of the celestial sphere, that is, the center of the earth.

After the coordinate transformation,

$$\left[\begin{array}{c}\cos {\delta }_{star}\cos {\alpha }_{star}\\ \cos {\delta }_{star}\sin {\alpha }_{star}\\ \sin {\delta }_{star}\end{array}\right]=k{A}_1\left[\begin{array}{c}\cos {\delta }_{earth}\cos {\alpha }_{earth}\\ \cos {\delta }_{earth}\sin {\alpha }_{earth}\\ \sin {\delta }_{earth}\end{array}\right]+r_{sat}$$

(7)

The relationship between the direction vector of the debris under the celestial sphere, that is, under the center of the earth, and the direction vector of the debris by satellite observation is the coordinate transformation from the coordinate system of the star tracker to the coordinate system of the geocentric coordinate system, where k is the scaling factor, $A_1$ is the coordinate system rotation matrix, $r_{sat}$ is the vector of the observation satellite (quasi-satellite optical observation platform) in the geocentric coordinate system, that is, the translation vector from the star tracker coordinate system to the geocentric coordinate system [5]. Since translation and scaling do not change the directionality of the vector, it can be derived

$$\left[\begin{array}{c}\cos {\delta }_{star}\cos {\alpha }_{star}\\ \cos {\delta }_{star}\sin {\alpha }_{star}\\ \sin {\delta }_{star}\end{array}\right]=A_1\left[\begin{array}{c}\cos {\delta }_{earth}\cos {\alpha }_{earth}\\ \cos {\delta }_{earth}\sin {\alpha }_{earth}\\ \sin {\delta }_{earth}\end{array}\right]$$

(8)

Then, after the observation vector under the geocentric observation point is rotated, it can be transformed into the observation vector under the satellite observation point. The rotation matrix is the rotation transformation matrix from the geocentric coordinate system to the star tracker coordinate system.

Obtained by the reversibility of the rotation matrix

$$\left[\begin{array}{c}\cos {\delta }_{earth}\cos {\alpha }_{earth}\\ \cos {\delta }_{earth}\sin {\alpha }_{earth}\\ \sin {\delta }_{earth}\end{array}\right]=A_1{}^{-1}\left[\begin{array}{c}\cos {\delta }_{star}\cos {\alpha }_{star}\\ \cos {\delta }_{star}\sin {\alpha }_{star}\\ \sin {\delta }_{star}\end{array}\right]$$

(9)

It can be concluded from the above equation that any satellite observation vector of the debris can be rotated to obtain the observation vector of the debris based on the center of the earth, and this property is defined as the invariability of the geocentric observation vector.

With the observation model of the star tracker, the rectascension and declination of the observed debris can be obtained, and it can be obtained according to Formulas (4) and (5)

$${\alpha }_{star}=\arctan \frac{y_{deb}-y_{sat}}{x_{deb}-x_{sat}}$$

(10)

$${\delta }_{star}=\arctan \frac{z_{deb}-z_{sat}}{\sqrt{{\left(x_{deb}-x_{sat}\right)}^2+{\left(y_{deb}-y_{sat}\right)}^2}}$$

(11)

Then the observation equation is derived:

$${\alpha }_d=\arctan \frac{y_{deb}-y_{sat}}{x_{deb}-x_{sat}}+{\nu }_{\alpha }$$

(12)

$${\delta }_d=\arctan \frac{z_{deb}-z_{sat}}{\sqrt{{\left(x_{deb}-x_{sat}\right)}^2+{\left(y_{deb}-y_{sat}\right)}^2}}+{\nu }_{\delta }$$

(13)

The instantaneous direction of the moving target in the space (rectascension, declination), brightness, velocity, and motion trajectory can be obtained by a single star tracker. From the joint calculation of the instantaneous pointing and the motion trajectory, the orbital inclination of the target can be calculated. Because of the high pointing accuracy of the star tracker, the more precise orbit inclination calculated with its data will be used as the basis for the main target recognition and classification index. Because the brightness information is closely related to the observation angle and the angle of the sun, even if certain space debris can simultaneously be observed by two star trackers, its brightness information will also vary greatly due to the influence of the observation angle and the sun angle [26]. Therefore, it is only a reference value and can be used to evaluate the target size information after the completion of orbit determination. The speed is the relative speed, which is the projection of the space debris velocity vector on the satellite optical observation platform detector. Because of the different angles of the two satellite observation positions, the relative velocity will vary greatly and is not significant for the classification and recognition by the two satellite observations.

4. Target Identification Method between Star Trackers

The proposed target identification method in this paper is combined with the star tracker technology. First, the proposed maximum projection method is used for star detection, where the median frame image obtained during star detection is removed from the maximum projection image, and the stars in the image can be removed. Next, the target detection method based on the projection time information is used to detect the moving targets using the trajectory of the moving targets as well as the continuity of the time information, while eliminating the false targets. After the detection of the targets in frame groups is completed, the target information array of each frame group can be obtained, and then the inter-frame target trajectory association method is used for the trajectory association of the target information to form a complete target trajectory. Finally, the identification method of targets observed by the same orbital arc segment as well as different orbital arc segments by using the star tracker is analyzed.

4.1. Maximum Projection Method

For k frames of $M\times N$ astronomical images, denote the data samples as $f(i,j,k)$, where $i=1,\dots ,M$, $j=1,\dots ,N$, $k=1,\dots ,K$. Assume that the set of images satisfies the following conditions:

The data samples taken from the image sequence are mutually independent Gaussian random variables and the noise can be considered Gaussian white noise.

If the moving target does not exist, the data sample can be regarded as a temporally stable but spatially unstable data model. This means that all images are spatially registered and there is no motion clutter.

The maximum projection of K samples is expressed as follows:

$$z(i,j)=\max [f(i,j,1)+f(i,j,1),\dots ,f(i,j,K)]$$

(14)

The above equation implies that the 3D data are projected onto the 2D plane. To describe this process, we assume that there are three consecutive frames of images, as shown in Figure 4a–c below, in which the information of the image is the same except for the moving target and some perturbations. After the maximum projection, the motion of the target is connected, as shown in Figure 4d below. Most importantly, in addition to the maximum projection information, we can obtain the time scale corresponding to the maximum value generated at each pixel position, as shown in Figure 4e below. It can be used to estimate the velocity of the target’s motion. Figure 5 shows the image after the projection of the maximum value.

The median frame image obtained during stellar detection is removed from the maximum projection image to remove the stars from the image.

The maximum projection image after background removal is binarized by adaptive threshold segmentation.

In addition to the maximum projection frame, it is also possible to obtain the time scale corresponding to the maximum value on each pixel during projection, i.e., the time information of the maximum projection frame $\max posi(i,j)$. Binarize the maximum projected image. After being separated from the image noise, and then multiplied with the time marker information, the time information frame of the maximum projected image can be obtained (as shown in Figure 6). The formula is as follows:

$$\begin{array}{l}BWz(i,j)=\{\begin{array}{l}0, {{\displaystyle }}_{}{{\displaystyle }}_{}x<threshold\\ 1, {{\displaystyle }}_{}{{\displaystyle }}_{}z(i,j)\le threshold\end{array}\\ timetag(i,j)=\max posi(i,j)\times BWz(i,j)\end{array}$$

(15)

From Figure 7, it can be seen that in the temporal information frame of the maximum projection, the moving targets are connected domains where the information with the same time markers are clustered together, while the stars are connected domains where the different temporal information are clustered together. After stellar removal, the time information and position information of the moving target can be used for target detection.

4.2. Target Detection Based on Projection Time Information

On the maximum projected frame after star removal, only candidate moving targets remain, which include stars that are not removed, and false target points composed of isolated noise, as shown below. In this paper, we propose a target detection method based on temporal information, which uses the trajectory of moving targets as well as the continuity of temporal information to detect moving targets while removing false targets. The principle of the algorithm is shown in Figure 8 and Figure 9. Figure 9 shows the process of setting a tracking window based on speed to detect possible objects, and the specific steps are as follows.

(1) The search starts with the candidate target $P3$ with time number T = 3, and the search radius is set to R (R is the maximum motion speed threshold of the moving target). When the candidate target points of T = 2 and T = 4 are found in the search range, the distance $D_{23}$ between T = 3 and T = 2 is calculated, and the distance $D_{34}$ between T = 3 and T = 4 is calculated. If $\left|D_{23}-D_{34}\right|\le D_{th}$, (where $D_{th}$ is the distance threshold, because the target is moving at a uniform speed, so $D_{th}$ is set to a very small value), then P2, P3, and P4 are judged to be a candidate trajectory; if the distance condition is not satisfied, the trajectory is judged as false trajectory.

(2) Determine the angle of the candidate trajectories $P2$, $P3$, and $P4$ that meet the distance judgment conditions, and calculate the angle between the straight line composed of $P2$ and $P3$ and the straight line composed of $P3$ and $P4$ and denote it as …. If $Angl{e}_{234}\le Angl{e}_{th}$ (where $Angl{e}_{th}$ is the angle threshold. Because the target is moving in a straight line, it is set to a small value), then confirm that $P2$, $P3$, and $P4$ are a true track, if not, it is judged as a false trajectory.

(3) The estimated velocity of the moving target is obtained by calculating the mean value of $D_{23}$, $D_{34}$: $V=(D_{23}+D_{34})/2$. The forward and backward searches are performed on the trajectories of $P2$, $P3$, and $P4$ to obtain the estimated positions $P{1}_e$ and $P{5}_e$ of $P1$ and $P5$. A tracking window of a certain size is established with the estimated position as the center, and the detection of targets within the window is performed by binarization and center-of-mass localization to obtain all candidate targets within the window. If a candidate target, $P$, satisfies $\left|P-P_e\right|\le P_{th}$ (where $P_e$ is the estimated position $P{1}_e$, and $P_{th}$ is the position error threshold), then $P$ is determined to be the true position of $P1$ or $P5$. In this way, the position of $P1$, $P5$ can be detected.

(4) Use the time information to judge the time markers of $P1$–$P5$ for the trajectory again and declare it as the real trajectory if it is the relation of contact increment.

4.3. Identification of Targets Observed in the Same Orbital Arc

The effectiveness analysis of information obtained from a single star tracker shows that the magnitude brightness and velocity components of the debris are different with the observation position when the debris is observed with respect to different star trackers, making it impossible to be used as the debris identification characteristic between the star trackers. The debris identifies the feature value. However, the geocentric observation vector of the debris can be transformed by the observation vector of the star tracker and the rotation matrix between the geocentric coordinate system and the coordinate system of the star tracker based on the invariance of the geocentric observation of the debris, and the calculation of the rotation matrix can be calculated by star image recognition and attitude of the star tracker. Thus, the identification of the debris image trajectory observed by the star tracker can be transformed into the similarity of the geocentric observation vector. The debris observed by different star trackers is identified by comparing the number of similar points of the motion trajectories of the debris observed by the star tracker.

The vector C is defined as follows:

$$b=\frac{1}{\sqrt{{\alpha }^2+{\delta }^2}}\left[\begin{array}{c}\alpha \\ \delta \end{array}\right]$$

(16)

$\alpha$, $\delta$ are the rectascension and declination of the debris, respectively.

The similarity criterion uses the Bhattacharyya coefficient, and Bhattacharyya is defined as follows:

$$\widehat{\rho }(b{}^{\prime })\equiv \rho \left[b,b{}^{\prime }\right]={\displaystyle \sum _{m=1}^n\sqrt{b(m)b^{\prime }(m)}}$$

(17)

It can be seen from the definition of the Bhattacharyya coefficient that the closer the two modes are, the larger the value $\widehat{\rho }(b{}^{\prime })$ is, and the above formula can be transformed into the following:

$$\xi (b{}^{\prime })=1-\widehat{\rho }(b{}^{\prime })$$

(18)

The more similar the measurement mode structure to the standard mode structure is, the smaller the value $\xi (b{}^{\prime })$ is. When $\xi (b^{\prime })<\epsilon$, it can be judged as a point in a similar trajectory.

Figure 10 shows the simulating orbital trajectory of space debris in 100 different motion velocities. The 100 sets of debris are identified by similarity criteria, and the error of 20 s of arc is added to 100 sets of data to form another set of observation data of the star tracker. After 100 measurements, the recognition rate is above 85%.

4.4. Identification of Targets Observed in the Different Orbital Arc

The identification of targets observed in the different orbital arc segments is done based on the existing multiple star trackers in the respective orbital arc segments, and the star trackers can complete the solution of the positioning information to identify the same debris targets in different orbital arc segments space debris. The observation correlation of different orbital arcs is low, and the similarity cannot be used to extract the same debris targets. The initial orbital parameters of space debris can only be solved by the positioning information of each arc segment, and then the similarity of orbital parameters can be used to identify the target in different orbital arcs. The orbit determination algorithm will be described in the next section.

5. Positioning and Orbit Determination of Space Targets

In this section, a joint position method using multiple star trackers is proposed to achieve the target position of multiple optical observations using a Gaussian minimum mean square error difference correction algorithm, fusing the spatial target position information of multiple frames of images and completing the solution of the orbital equations. The schematic diagram of multiple star trackers observing debris simultaneously is shown in Figure 11:

From Equations (10) and (11), it can be seen that the right ascension and declination information of the observed debris through the star tracker has a high order nonlinear relationship with the debris orbital position information, which cannot be obtained with high accuracy by using the traditional geometric solution. From the observation equation, a high-precision static equation solving algorithm is introduced to solve the problem to filter out the noise to the maximum extent and achieve the optimal estimation result.

5.1. Gaussian Minimum Mean Square Error Correction Positioning Algorithm

The optimal solution of the nonlinear static equation can be determined with the Gaussian minimum mean square error differential correction algorithm [26]. The loss equation determined by the Wahba minimum mean square error is as follows:

$$J=\frac{1}{2}{{\displaystyle \sum _{i=1}^N\left({\tilde{b}}_i-h_i\left(\widehat{x}\right)\right)}}^T{W}_i\left({\tilde{b}}_i-h_i\left(\widehat{x}\right)\right)$$

(19)

The optimization of the loss equation can be attributed to the minimization of the redundant mean square error with the weight. Where, $i$ is a single measurement label, and $N$ is the total number of available measurements. $W=W^T<0$ is a weight coefficient matrix to assess the relative importance of each measurement. Estimated by the minimum variance estimation criterion, the optimal weight coefficient matrix $W$ is the inverse matrix of the measurement covariance matrix $R^{-1}$.

The Gaussian minimum mean square error differential correction algorithm uses an iterative approximation to find accurate minimum mean square estimates.

Assuming the estimate of the current state is ${\widehat{x}}_k$, then the next target state estimates ${\widehat{x}}_{k+1}$ can be expressed as follows:

$${\widehat{x}}_{k+1}={\widehat{x}}_k+\Delta x$$

(20)

If $\Delta x$ is small enough, $h(\widehat{x})$ can be linearized with Taylor’s first-order expansion linearization.

$$h_i\left(\widehat{x}\right)\approx h_i\left({\widehat{x}}_k\right)+H_{i,k}\Delta x$$

(21)

where $H_{i,k}$ is the Jacobian matrix of the measurement model.

$$H_{i,k}={\left[\frac{\partial h_i}{\partial x}\right]}_{{\widehat{x}}_k}$$

(22)

Then the redundancy error of the state estimation can be expressed as follows:

$$\begin{array}{l}\Delta b_{i,k+1}={\tilde{b}}_i-h_i\left({\widehat{x}}_{k+1}\right)\\ {{\displaystyle }}^{}\approx {\tilde{b}}_i-h_i\left({\widehat{x}}_k\right)-H_{i,k}\Delta x\\ {{\displaystyle }}^{}=\Delta b_{i,k}-H_{i,k}\Delta x\end{array}$$

(23)

The measurement vector and the prediction vector can be expressed as follows:

$$\tilde{Y}=\left[\begin{array}{c}{\tilde{b}}_1\\ {\tilde{b}}_2\\ ⋮\\ {\tilde{b}}_N\end{array}\right],{\widehat{Y}}_k=\left[\begin{array}{c}{h}_1\left({\widehat{x}}_k\right)\\ h_2\left({\widehat{x}}_k\right)\\ ⋮\\ h_N\left({\widehat{x}}_k\right)\end{array}\right]$$

(24)

Then the measurement error and sensitivity matrix can be expressed as follows:

$$\Delta Y_k=\tilde{Y}-{\widehat{Y}}_k=\left[\begin{array}{c}{\tilde{b}}_1-h_1\left({\widehat{x}}_k\right)\\ {\tilde{b}}_2-h_2\left({\widehat{x}}_k\right)\\ ⋮\\ {\tilde{b}}_N-h_N\left({\widehat{x}}_k\right)\end{array}\right],A_k=\left[\begin{array}{c}{H}_1\\ H_2\\ ⋮\\ H_N\end{array}\right]$$

(25)

To minimize the loss Equation of (19), that is, the optimal correction amount is searched to minimize the linearized prediction redundancy mean square error:

$$J_{k+1}\approx \frac{1}{2}{\left(\Delta Y_k-A_k\Delta x\right)}^{T}W\left(\Delta Y_k-A_k\Delta x\right)$$

(26)

Solve the above formula, the optimal approximate solution is

$$\Delta x=P_k{A}_k^{T}W\Delta Y_k$$

(27)

where $P_k$ is the covariance matrix,

$$P_k={\left(A_k^{T}W{A}_k\right)}^{-1}$$

(28)

Therefore, ${\widehat{x}}_k$ of Formula (20) will be updated to be ${\widehat{x}}_{k+1}$.

An initial estimation state ${\widehat{x}}_0$ is required to start the entire Gaussian minimum mean square error differential correction algorithm, and an iterative condition is required to stop the entire algorithm. The iteration condition can be expressed by:

$$\delta J\equiv \frac{\left|J_k-J_{k-1}\right|}{J_k}<\frac{\epsilon }{\Vert W\Vert }$$

(29)

Orbital position solution, the Equations (12) and (13) can be transformed as follows:

$$\tan \alpha =\frac{y_{deb}-y_{sat}}{x_{deb}-x_{sat}}+{\nu }_{\alpha }$$

(30)

$$\tan \delta =\frac{z_{deb}-z_{sat}}{\sqrt{{\left(x_{deb}-x_{sat}\right)}^2+{\left(y_{deb}-y_{sat}\right)}^2}}+{\nu }_{\delta }$$

(31)

${\nu }_{\alpha }$, ${\nu }_{\delta }$ is the tangent error of the rectascension and declination caused by the measurement error of the satellite optical observation platform, respectively.

Let vector $Y={[\tan \alpha {{\displaystyle }}^{}\tan \delta ]}^T$, then we have:

$$Y={[\tan \alpha {{\displaystyle }}^{}\tan \delta ]}^T=\left[\begin{array}{c}\frac{y_{deb}-y_{sat}}{x_{deb}-x_{sat}}\\ \frac{z_{deb}-z_{sat}}{\sqrt{{\left(x_{deb}-x_{sat}\right)}^2+{\left(y_{deb}-y_{sat}\right)}^2}}\end{array}\right]$$

(32)

The measured sensitivity matrix is:

$$\begin{array}{l}H=\frac{\partial h_{}}{\partial r_{deb}}\\ {{\displaystyle }}^{}=\left[\begin{array}{ccc}\frac{\partial h_1}{\partial x_{deb}}& \frac{\partial h_1}{\partial y_{deb}}& \frac{\partial h_1}{\partial z_{deb}}\\ \frac{\partial h_2}{\partial x_{deb}}& \frac{\partial h_2}{\partial y_{deb}}& \frac{\partial h_2}{\partial z_{deb}}\end{array}\right]\\ {{\displaystyle }}^{}=\left[\begin{array}{ccc}-\frac{y_{deb}-y_{sat}}{{\left(x_{deb}-x_{sat}\right)}^2}& \frac{1}{x_{deb}-x_{sat}}& 0\\ -\left(z_{deb}-z_{sat}\right)\left(x_{deb}-x_{sat}\right)A^3& -\left(z_{deb}-z_{sat}\right)\left(y_{deb}-y_{sat}\right)A^3& A\end{array}\right]\end{array}$$

(33)

where $A=\frac{1}{\sqrt{{\left(x_{deb}-x_{sat}\right)}^2+{\left(y_{deb}-y_{sat}\right)}^2}}$, vector $Y$ is a nonlinear equation about the debris vector $r_{deb}={[x_{deb }{{\displaystyle }}^{}{y}_{deb }{{\displaystyle }}^{}{z}_{deb }]}^T$. In order to solve the debris vector, three or more measurement vectors are required, i.e., more than three satellites are needed to observe the debris.

From the measurement sensitivity matrix, the error transfer function for the right ascension and declination is obtained as follows:

$$\begin{array}{l}\left[\begin{array}{c}\Delta \alpha \\ \Delta \delta \end{array}\right]=H\left[\begin{array}{c}\Delta x_{rdeb}\\ \Delta y_{rdeb}\\ \Delta z_{rdeb}\end{array}\right]\\ {{\displaystyle }}^{}{{\displaystyle }}^{}=\left[\begin{array}{ccc}-\frac{y_{deb}-y_{sat}}{{\left(x_{deb}-x_{sat}\right)}^2}& \frac{1}{x_{deb}-x_{sat}}& 0\\ -\left(z_{deb}-z_{sat}\right)\left(x_{deb}-x_{sat}\right)A^3& -\left(z_{deb}-z_{sat}\right)\left(y_{deb}-y_{sat}\right)A^3& A\end{array}\right]\left[\begin{array}{c}\Delta x_{rdeb}\\ \Delta y_{rdeb}\\ \Delta z_{rdeb}\end{array}\right]\end{array}$$

(34)

The smaller the $x_{deb}-x_{sat}$, the larger the $y_{deb}-y_{sat}$, the larger the error in calculating the declination, and the smaller the $A$, the larger the $z_{deb}-z_{sat}$, the larger the error in calculating the right ascension.

Since the observation quantities of multiple star trackers are needed to perform the solution, the observation time synchronization for the star tracker observation quantities is particularly important. The instantaneous time error of the observation quantities includes the satellite-borne clock error, the clock transmission–reception errors, and the tracker information acquisition time errors. In terms of the tracker information acquisition time error, the error is on the order of 10⁻⁶ s, the satellite-borne clock error is about 10⁻¹¹ s, and the clock transmission reception error increases with the increase in distance, these time errors bring great difficulties to the time synchronization of the observation quantity. However, for space debris, no matter how many satellites observe, the spatial position of debris at the same time is unique, i.e., geocentric observation invariance. By comparing the motion trajectories of debris observed by different star trackers, the right ascension and declination of the motion trajectories are compared, and the two sets of measurements with the closest declination and declination are taken as a time synchronization point, and two time synchronization points in a set of motion trajectories are selected as the starting and ending points of the trajectory observations. Since the debris observed in the star tracker image is a very short orbital arc segment, and its motion trajectory is a straight line motion, the starting and ending points of the trajectory observation are re-fitted differentially according to the customized time to obtain a new measurement quantity, which can reduce the error caused by the measurement, and also solve the difficulty of not being able to perform high precision orbit determination due to the insufficient measurement quantity. The arc segment with an observation duration of less than 180 s is defined as a short arc. In order to highlight the effect of short arc correlation, trajectory data with an observation duration of 50~120 s are selected for correlation. The revisit time of general space targets is 20~50 min. The longer the trajectory interval that can be associated, the higher the utilization rate of data. The maximum interval in the simulation is 10 h, which can meet most of the requirements of interrupted trajectory association.

5.2. A Model for Determining Orbital Parameters of Space Debris Based on Space-Based Observations

Accurate evaluation of orbital precision is very difficult because, first, there is no way to know the true criterion for precision assessment, the true value of the orbit, and secondly, it is often difficult to obtain higher precision orbits that can be used as a criterion. Therefore, the evaluation of orbital precision can only be an approximately statistical result under certain assumptions. The orbital precision can be evaluated by comparing the instantaneous orbital roots, the average orbital roots, or the spatial debris position and velocity components of the time series [27].

In this paper, referring to the Gooding initial orbit determination method of space-based optical measurement, three sets of angular measurement data and the position information of the measurement platform are used to estimate the position and velocity of the space target. In addition, the space coordinate position is converted into six orbital elements using the star tracker coordinates.

The orbital elements of elliptical motion are commonly expressed in terms of six Kepler orbital elements (as shown in Figure 12), namely: orbital half-length diameter $a$, orbital eccentricity $e$, argument of perigee $\omega$, orbital inclination $i$, longitude of ascending node $\Omega$ and mean anomaly $M$.

$M$ is an orbital element describing the position of the satellite in orbit:

$$M=n(t-\tau )$$

(35)

where $n=\sqrt{\frac{\mu }{a^3}}$, $n$ is often called the translation and $\tau$ is the moment of perigee of the debris.

The angle between the directional diameter $\stackrel{\rightharpoonup }{r}$ and the perigee is called the true anomaly $f$. The orbital equation and energy integral are then converted into:

$$\begin{array}{l}r=\frac{a\left(1-e^2\right)}{1+e\cos f}\\ v^2=\mu \left(\frac{2}{r}-\frac{1}{a}\right)\end{array}$$

(36)

The spatial coordinate position is transformed into the six orbital elements. The coordinates $\stackrel{\rightharpoonup }{r}$ and $\stackrel{\rightharpoonup }{v}$ at a given moment are known, and the six orbital elements are calculated as follows:

①

Calculate $a$

$$a=\frac{\mu r}{2\mu -r{v}^2}$$

(37)

②

Calculate $e\sin E,e\cos E,e,E$

$$\begin{array}{l}e\sin E=\frac{\stackrel{\rightharpoonup }{r}\cdot \stackrel{\rightharpoonup }{v}}{n{a}^2}\\ e\cos E=1-\frac{r}{a}\end{array}$$

(38)

③

Calculate $i,\Omega$

$$\begin{array}{l}\stackrel{\rightharpoonup }{h}=\stackrel{\rightharpoonup }{r}\times \stackrel{\rightharpoonup }{v}=\sqrt{\mu p}\left[\begin{array}{c}\sin \Omega \sin i\\ -\cos \Omega \sin i\\ \cos i\end{array}\right]\\ \cos i=\frac{h_3}{\sqrt{\mu p}}\\ \tan \Omega =-\frac{h_1}{h_2}\end{array}$$

(39)

④

Calculate $\omega$

The true anomaly angle $f$ can be calculated according to the following formula.

$$f=E+2{\tan }^{-1}\left(\frac{e\sin E}{1+\sqrt{1-e^2}-e\cos E}\right)$$

(40)

The eccentric anomaly $u$ can be calculated according to the following formula:

$$\left[\begin{array}{c}r\cos u\\ r\sin u\\ 0\end{array}\right]=R_1(i)R_3(\Omega )\stackrel{\rightharpoonup }{r}$$

(41)

$$\omega =u-f$$

(42)

⑤

Calculate mean anomaly $M$

$$M=E-e\sin E$$

(43)

By calculating the positioning algorithm in the previous section, a higher accuracy position determination of the target can be accomplished. If the spatial positions of the target at two moments can be obtained, the initial orbital parameters of the space debris can be calculated by Equations (37)–(43). If there are multiple measurements of space debris in one orbital arc, the initial orbital elements can be further obtained with higher accuracy by calculating the average orbital elements.

The location information of each orbital arc segment of space debris can be solved through the positioning algorithm in the previous section. On each small orbital arc, the motion trajectory can be considered as a uniform linear motion. Since the positioning information has a certain error, if the position between two points is directly divided by the elapsed time, the introduced speed error will be larger. About this point, we calculate the average speed of the arc through the speed fitting of each small arc, and the error will be minimized.

6. Positioning Algorithm Analysis

The satellite tool kit (STK) platform was used to build the observation model. The observation platform was located in the solar synchronous orbit, and the orbit parameters refer to the observation data of the union of American Yousi scientists. The space debris position vector $r_{deb}$ was generated by simulation to traverse the low orbit and high orbit. The instantaneous position height was from 6800 km to 40,000 km, and the orbit inclination was 0°~2°. The orbital inclination range of the observation platform was 96.8°~98.5°, and the eccentricity was 0. The positions of three star trackers were randomly generated according to the space debris position vector, and the observation distance from the space debris was guaranteed. A total of 200 cooperative positioning of space debris were carried out, and the positioning errors of star trackers at different heights were calculated at the observation distances of 200 km, 100 km, and 50 km. In order to simulate the attitude determination error of the star tracker and the position error of the platform itself in the task, the position errors of 10 m and 100 m were introduced into the position of the star tracker in the geocentric coordinate system according to the performance of the existing star tracker, and the pointing accuracy of the star tracker was set to 2″, which is more in line with the practical application. The code had been verified on computers with an Intel Core i7-6500u (dual cores) CPU up to 3.1 GHz, 8 GB Ram, and operating system Windows 10. The simulation results are shown in Figure 13 and Figure 14.

It can be seen from the above curve data that the positioning accuracy of the debris also changed with the change in the positioning accuracy of the satellite itself. When the positioning accuracy of the satellite itself was 10 m, the positioning accuracy of the debris was also about 10 m, as shown in Figure 13. It can be found from Figure 14 that when the positioning accuracy of the satellite itself was 100 m, the positioning accuracy of the debris was about 100 m. For the positioning of space debris with a plurality of optical sensors, the positioning accuracy varied with the positioning accuracy of the observation satellite itself, and the positioning accuracy of the space debris was limited by the positioning accuracy of the satellite itself.

By calculating the positioning algorithm, a higher accuracy position determination of the target can be accomplished. If the spatial positions of the target at two moments can be obtained, the initial orbital parameters of the space debris can be calculated by Equations (37)–(43). If there are multiple measurements of space debris in one orbital arc, the initial orbital elements can be further obtained with higher accuracy by calculating the average orbital elements. The space debris was in different orbits, and the positioning accuracy of space debris changed due to the influence of the positioning accuracy of the observation star tracker itself, and the positioning accuracy of the high-orbit satellite itself was poor compared with that of the low-orbit satellite. The following simulation calculations were completed by using 20 sets of position measurements in one orbital arc segment to solve the orbital elements, the simulation results are shown in Figure 15 and Figure 16.

As shown in Figure 15 and Figure 16, the orbital parameters were estimated by using the arc segments distributed at different positions and the 20-position information of each arc segment with a position error of 10 m, where there were 20 sets of measurements for each arc segment. The real orbital parameters were $a=20,000,e=0.003,i=40,\Omega =40,w=70,M=0$. The initial orbital parameters estimated using one arc segment had high accuracy. Figure 15 indicates that the maximum error of orbital half-length diameter was 0.5 km; the maximum error of orbit eccentricity was 1.5 × 10⁻⁴, and the error of orbit inclination was below 6 × 10⁻⁶⁰; the error of longitude of ascending node was also below 1 × 10⁻⁴⁰. With the help of orbital parameter estimation, the initial orbital parameters can be used as feature vectors to accomplish target identification for different orbital arc segments. From Figure 16, it can be seen that the longitude of orbit parameters still maintained high calculation accuracy as the orbit height decreased, The real orbital parameters were $a=7000,e=0.002,i=40,\Omega =40,w=70,M=0$, the maximum error of orbital half-length diameter was below 0.1 km, the maximum error of orbit eccentricity was 6 × 10⁻⁴, and the error of orbit inclination and the error of longitude of the ascending node were both in the order of 1 × 10⁻⁴⁰.

It can be seen from the positioning accuracy of space debris, the positioning accuracy of the star tracker itself largely determined the observation accuracy of the debris. Only by improving the positioning accuracy of the satellite itself, can the positioning accuracy of the observation target (that is, the space debris) be achieved, and finally the high-precision orbit determination of the target can be completed.

7. Conclusions

This paper presented an innovative strategy for space debris detection using multiple star trackers for space-based observations. It was demonstrated that the typical optical and electronic characteristics of current star trackers are sufficient for space debris detection. This paper demonstrated in detail the positioning method and the observability of space target debris by a multi-star tracker network. The method takes advantage of the high accuracy of star trackers to achieve observable positioning of space debris smaller than 10 cm. Through the study, the conclusion of space debris geocentric observation vector invariance was derived, and this conclusion was used as a breakthrough to solve the problem of efficient identification of multi-star tracker targets with the same orbit identification rate greater than 80%. Then, a positioning algorithm and an orbit determination algorithm were further proposed. The results of the algorithms in this paper are encouraging and the proposed method is sufficient to detect and position space debris. The recognition of targets in different orbit arcs is completed through the initial orbit parameters. The final positioning accuracy was better than 15 m, and the orbit determination accuracy was better than 1000 m at high orbit and better than 100 m at low orbit. This method has the potential to successfully determine the orbit of space debris, and this research result can provide a reference for the orbit identification of space debris as well as early warning and avoidance. The current engineering star tracker is only used as an attitude sensor application, and the combination of the measured data and simulation software can provide a reference for subsequent engineering applications. With the development of new star trackers with a small field of view and high limiting magnitude, further space debris positioning studies will be conducted with a real star map tracker.