A New Triangulation Algorithm for Positioning Space Debris

: A single electro-optical (EO) sensor used in space debris observation provides angle-only information. However, space debris position can be derived using simultaneous optical measurements obtained from two EO sensors located at two separate observation sites, and this is commonly known as triangulation. In this paper, we propose a new triangulation algorithm to determine space debris position, and its analytical expression of Root-Mean-Square (RMS) position error is presented. The simulation of two-site observation is conducted to compare the RMS positioning error of the proposed triangulation algorithm with traditional triangulation algorithms. The results show that the maximum RMS position error of the proposed triangulation algorithm is not more than 200 m, the proposed triangulation algorithm has higher positioning accuracy than traditional triangulation algorithms, and the RMS position error obtained in the simulation is nearly consistent with the analytical expression of RMS position error. In addition, initial orbit determination (IOD) is carried out by using the triangulation positioning data, and the results show that the IOD accuracy of two-site observation is signiﬁcantly higher than that of the single-site observation.


Introduction
Space debris is defunct artificial objects in space, and it began to accumulate in Earth's orbit immediately with the first launch of an artificial satellite, Sputnik 1, into orbit in October 1957. As of August 2021, more than 900,000 pieces of debris larger than 1 cm are estimated to be in orbit around the Earth [1]. According to some computer models, the amount of space debris has reached a tipping point, with enough currently in orbit to continually collide and create even more debris, raising the risk of spacecraft failures [2,3]. Due to the potential threat that space debris poses, space debris observation and orbit determination are increasingly needed.
Optical observation is now widely used in space debris observation because it is unrestricted by the distance to observation objects and is cost-efficient compared to radar, and several countries and organizations have carried out optical space debris observation [4][5][6][7][8][9][10][11][12][13][14][15][16]. Electro-optical (EO) sensors are the primary instrument for optical observation systems. Optical observation systems use EO sensors to collect images of space debris and extract angle information from the collected images [17]. Optical observation does not intrinsically provide range information, and therefore we cannot directly obtain space debris position by a single EO sensor; however, space debris position can be derived using simultaneous optical measurements obtained from two EO sensors located at two separate observation sites, and this is commonly known as triangulation [18,19]. Yanagisawa et al. propose a ground-based optical observation system for monitoring and positioning low earth orbit (LEO) objects, which consists of two arrays of EO sensors located at two longitudinally separate observation sites [20]. In the absence of observation error, the triangulation problem is trivial; however, when observation error is present, the line-of-sight (LOS) of two sensors will not generally meet, in which case it is necessary to find the best point of intersection as Remote Sens. 2021, 13, 4878 2 of 12 the most probable space debris position. There are several algorithms to solve this problem. For example, a commonly suggested algorithm [21] takes the midpoint of the shortest line between two lines of sight as the most probable space debris position; in [22], the most probable space debris position is taken as a point on the shortest line between two lines of sight, which subtend equal residual angles at sites. Obtaining the positioning errors in the triangulation solution is a complicated problem, and numerical solutions are usually sought, but it is difficult to fully reflect the nature of the dependence between errors and experimental parameters involved in numerical studies. In [23], the analytical expression for the Root-Mean-Square (RMS) position error of triangulation is presented.
Triangulation has many applications, and one is to improve the accuracy of shortarc angles-only initial orbit determination. Initial orbit determination (IOD) refers to determining the orbit of the space objects using observation data, which is the key to space debris cataloging. Affected by the relative velocity between space debris and the EO sensors, the observed arcs are generally too short, especially for ground-based observation of LEO space debris or space-based observation [24,25]. Since the observed arcs are too short and lack range information, the short-arc angles-only initial orbit determination results generally have large errors [26]. However, based on triangulation using two EO sensors, the angles-only IOD problem is converted to the IOD using 3D position data, and the IOD accuracy can be significantly improved.
To obtain better positioning results, it is desirable to find a triangulation algorithm that improves positioning results by introducing observation accuracy parameters. This paper addresses this problem by assuming that the space debris state observed by EO sensors is subject to Gaussian distribution, and derives the space debris position by the statistics method. In addition, the analytical expression for the RMS position error of the proposed triangulation algorithm is presented. In a series of simulations, the algorithm is extensively tested against other triangulation algorithms. The rest of this paper is organized as follows: in Section 2, a new triangulation algorithm is proposed, and its analytical expression of RMS position errors is presented. In Section 3, the simulation of two-site observation is conducted to compare the RMS positioning error of the proposed triangulation algorithm with traditional triangulation algorithms, and the results and discussion are presented in Sections 4 and 5, respectively. Finally, the conclusions are given in Section 6.

Triangulation Algorithm
In Figure 1, the two EO sensors located at two arbitrary observation sites are denoted by A and B, respectively. The LOS vectors from the EO sensors to space debris are denoted by L A and L B , respectively, and the line segment MN is the common perpendicular of them. Figure 2 illustrates the right ascension and declination (RADEC) of a LOS vector in the topocentric coordinate systems. With RADEC measurements, it is possible to calculate the LOS vector: The object in triangulation is to derive the space debris position X T from the space debris state X θ observed by two EO sensors. We address this problem by assuming that the space debris state observed by EO sensors is subject to Gaussian distribution, with mean value X θ = [α A , δ A , α B , δ B ], which are RADEC measurements of EO sensors, and variance σ θ , which is the uncertainty of RADEC measurements, and it can be assumed that RMS errors σ α A , σ δ A , σ α B and σ δ B are uncorrelated. According to the traditional triangulation algorithms, it is generally believed that the actual space debris position is near the line segment MN, and therefore it can be assumed that the position of point M is the mean value of the space debris position determined by the EO sensor A with the assistance of EO sensor B. Similarly, the position of point N is regarded as the mean  The object in triangulation is to derive the space debris position T X from the space debris state X θ observed by two EO sensors. We address this problem by assuming that the space debris state observed by EO sensors is subject to Gaussian distribution, with mean value       To determine the space debris position, the positions of points M and N and their variances need to be calculated. From Figure 1, according to the triangle law of vector addition, in geocentric coordinate systems, the expressions for X M and X N with respect to LOS vectors L A and L B are given by: where L A and L B can be calculated by Equation (1); geometry. Since line segment MN is the common perpendicular of L A and L B , we can obtain equation: Substituting Equations (2) and (3) into Equation (4) yields: Solving Equation (5) and after further simplification, we obtain: where → AB = X B − X A . Then, substituting Equations (6) and (7) into Equations (2) and (3), respectively, we obtain the expressions for X M and X N , respectively: Using Equation (1), Equations (8) and (9) can be further simplified to give the expressions for X M and X N with respect to variables α A , δ A , α B , and δ B , respectively.
According to the law of propagation error, the covariance matrix of point M can be calculated: Similarly, the covariance matrix of point N can be calculated: We now have two Gaussian distributions-the position of point M is subject to Gaussian distribution: and the position of point N is subject to Gaussian distribution: Both of them are the actual position estimates of space debris, and our goal is to obtain a new distribution that is the configuration for which both estimates are most likely, and Remote Sens. 2021, 13, 4878 5 of 12 the mean of the new distribution is the most probable space debris position given all the information we have. According to Bayes' theorem, the new distribution is given by: The new distribution is also a Gaussian distribution, and its mean value and variance are given by: where X T is the most probable space debris position obtained by using the proposed triangulation algorithm. The summary of the above triangulation algorithm (Algorithm 1) is listed below.

Algorithm 1. Summary of the proposed triangulation algorithm
Input: The two EO sensor positions X A and X A ; The two EO sensors' observation RMS errors σ α A ,σ δ A ,σ α B and σ δ B ; The simultaneous optical measurements X θ . Output: Space debris position X T . 1: Compute X M and X N from Equations (8) and (9). 2: Compute P M and P N from Equations (10) and (11). 1: Compute X T from Equation (15).

Positioning Error
According to previous studies, the positioning error of triangulation is not only related to the observation error of EO sensors but also depends on the angles θ A and θ B between the LOS vectors and the baseline of two EO sensors. In order to analyze the RMS position error of the proposed triangulation algorithm, we rewrite Equation (15) as: where K = P M −1 (P M −1 + P N −1 ) −1 is a weighting matrix and I is an identity matrix, and then the space debris position error is given by: where R M and R N are the position errors of points M and N, respectively, and the vectors → r M and → r N are their direction vectors. In order to obtain the expressions of RMS position error, we have to derive the expressions of R M and R N with respect to θ A and θ B .
Form Figure 3, the two EO sensors located at two arbitrary observation sites are denoted by A and B, respectively, and d is the distance between them. The actual position of space debris is denoted by T. The dθ A is the LOS vector error in plane θ A from site A, and dϕ A is the LOS vector error out of plane θ A from site A. Similarly, the dθ B is the LOS vector error in plane θ B from site B, and dϕ B is the LOS vector error out of plane θ B from site B. Form Figure 3, the position error of point M in triangulation is given by:  The RMS position error of point M in triangulation is then given by: It can be assumed that the LOS vector errors A dθ and A dϕ are uncorrelated.
Therefore, Equation (20) becomes: Similarly, the RMS position error of point N in triangulation is given by: From Figure 3, according to the law of sines, we obtain: and it is obvious that the error A dθ manifests itself in a change in the range B ρ , while the angle B θ is unaffected; therefore: Since A dϕ is small, The RMS position error of point M in triangulation is then given by: It can be assumed that the LOS vector errors dθ A and dϕ A are uncorrelated. Therefore, Equation (20) becomes: Similarly, the RMS position error of point N in triangulation is given by: From Figure 3, according to the law of sines, we obtain: and it is obvious that the error dθ A manifests itself in a change in the range ρ B , while the angle θ B is unaffected; therefore: Since dϕ A is small, Therefore, Substituting Equations (25) and (27) into Equation (20) yields: Similarly, Therefore, the expressions of RMS position errors can be derived by Equations (28), (29) and (18) Therefore, Equation (30) can be rewritten as: where → r MN = L A × L B because line segment MN is the common perpendicular of L A and L B .

Simulation Procedure
This section introduces the simulation procedure of space debris observation by triangulation. The purposes of this simulation are, on one hand, to compare the RMS positioning error of the proposed triangulation algorithm with traditional triangulation algorithms, and on the other hand, to evaluate the IOD errors using triangulation positioning data. Based on optical measurements from Jilin Space Tracking Base of Changchun Observatory, China, a ground-based simulation is conducted from Jan. 13, 2020, to Feb. 8, 2020. A total of 702 pieces of TLE-cataloged LEO space debris in the near-circular path are selected as the reference population, and the mean altitude distribution of these space debris is shown in Figure 4. The triangulation for positioning LEO space debris has the minimum error when the baseline between two EO sensors is about 1600 km [27]. Therefore, we choose the Changchun Observatory, China, and the Sheshan Observatory in Shanghai, China as the simulated observation sites where the EO sensors are set (the baseline between them is 1456.36 km). The simulation procedure is organized as follows: 1.
Based on optical measurements from Jilin Space Tracking Base of Changchun Observatory, precise orbit determination (POD) is performed on the selected LEO space debris to obtain their precise ephemeris; 2.
Based on the ephemeris obtained in step 1, the simultaneous observation arc consists of a series of consecutive simultaneous RADEC pseudo-measurements generated for the two simulation observation sites, and each RADEC pseudo-measurement is corrupted with Gaussian noise of 9". The length and sampling interval of the simultaneous observation arc are set at 30 s and 1 s, respectively. In addition, the detection constraints are also considered in the simulation, including the visibility of space debris to the observation sites, but not the EO sensors' performance and weather conditions; 3.
Based on the simultaneous RADEC pseudo-measurements obtained in step 2, the space debris position can be calculated by the proposed triangulation algorithm and two traditional triangulation algorithms, respectively. The positioning results are compared with the ephemeris, and the RMS position errors of the above three triangulation algorithms are evaluated; 4.
The initial orbit of selected LEO space debris is calculated using the obtained positioning results in step 3, the initial orbit results are compared with the reference orbit obtained in step 1, and IOD errors are evaluated. Due to the length of the detection arc, the Herrick-Gibbs method and differential correction are used in IOD processes [28].
debris to obtain their precise ephemeris; 2. Based on the ephemeris obtained in Step 1, the simultaneous observation arc consists of a series of consecutive simultaneous RADEC pseudo-measurements generated for the two simulation observation sites, and each RADEC pseudo-measurement is corrupted with Gaussian noise of 9″. The length and sampling interval of the simultaneous observation arc are set at 30s and 1s, respectively. In addition, the detection constraints are also considered in the simulation, including the visibility of space debris to the observation sites, but not the EO sensors' performance and weather conditions; 3. Based on the simultaneous RADEC pseudo-measurements obtained in step 2, the space debris position can be calculated by the proposed triangulation algorithm and two traditional triangulation algorithms, respectively. The positioning results are compared with the ephemeris, and the RMS position errors of the above three triangulation algorithms are evaluated; 4. The initial orbit of selected LEO space debris is calculated using the obtained positioning results in step 3, the initial orbit results are compared with the reference orbit obtained in step 1, and IOD errors are evaluated. Due to the length of the detection arc, the Herrick-Gibbs method and differential correction are used in IOD processes [28].

Results
In this section, the results for the RMS position error and IOD error of space debris described above in Section 3 are presented. In the 27-day simulation period, we obtained a total of 18,954 simultaneous observation arcs of the 702 pieces of TLE-cataloged LEO space debris, and each space object has more than 20 simultaneous observation arcs. Figure  5 shows the RMS position error using the proposed triangulation algorithm and two traditional triangulation algorithms, respectively. The RMS position error refers to the Root-Mean-Square position error of all sampling points on the simultaneous observation arcs. In addition, the analytical RMS position error calculated by Equation (32) is also plotted in Figure 5.
IOD is carried out by using the triangulation positioning data obtained in step 3. We calculate the total of 18954 IOD cases for 702 pieces of TLE-cataloged LEO space debris, and compare the results with the reference orbit obtained in step 1. Figures 6-8 show the IOD error, and we use orbit elements' semi-major axis (SMA), eccentricity, and inclination to characterize IOD errors. space debris, and each space object has more than 20 simultaneous observation arcs. Figure 5 shows the RMS position error using the proposed triangulation algorithm and two traditional triangulation algorithms, respectively. The RMS position error refers to the Root-Mean-Square position error of all sampling points on the simultaneous observation arcs. In addition, the analytical RMS position error calculated by Equation (32) is also plotted in Figure 5. IOD is carried out by using the triangulation positioning data obtained in step 3. We calculate the total of 18954 IOD cases for 702 pieces of TLE-cataloged LEO space debris, and compare the results with the reference orbit obtained in step 1. Figures 6-8 show the IOD error, and we use orbit elements' semi-major axis (SMA), eccentricity, and inclination to characterize IOD errors.   IOD is carried out by using the triangulation positioning data obtained in step 3. We calculate the total of 18954 IOD cases for 702 pieces of TLE-cataloged LEO space debris, and compare the results with the reference orbit obtained in step 1. Figures 6-8 show the IOD error, and we use orbit elements' semi-major axis (SMA), eccentricity, and inclination to characterize IOD errors.

Discussion
For LEO space debris in the simulation, Figure 5 shows that the maximum RMS position error of two traditional triangulation algorithms exceeds 250 m; however, the maximum RMS position error of the proposed triangulation algorithm is not more than 200 m. It is obvious that the RMS position errors of the proposed triangulation algorithm are better than those of the other two traditional triangulation algorithms. Therefore, the introduction of observation accuracy parameters into the triangulation algorithm can effectively improve the positioning accuracy. For the proposed triangulation algorithm, the RMS position error of the numerical solution is slightly higher than that of the analytical solution obtained by using Equation (32). The difference between them is mainly due to the approximate calculation in the analytic expression derivation of RMS position error and the uncertainty of RADEC measurements in numerical simulation. Ignoring these effects, the RMS position errors of the analytical solution and the numerical solution are in

Discussion
For LEO space debris in the simulation, Figure 5 shows that the maximum RMS position error of two traditional triangulation algorithms exceeds 250 m; however, the maximum RMS position error of the proposed triangulation algorithm is not more than 200 m. It is obvious that the RMS position errors of the proposed triangulation algorithm are better than those of the other two traditional triangulation algorithms. Therefore, the introduction of observation accuracy parameters into the triangulation algorithm can effectively improve the positioning accuracy. For the proposed triangulation algorithm, the RMS position error of the numerical solution is slightly higher than that of the analytical solution obtained by using Equation (32). The difference between them is mainly due to the approximate calculation in the analytic expression derivation of RMS position error and the uncertainty of RADEC measurements in numerical simulation. Ignoring these effects, the RMS position errors of the analytical solution and the numerical solution are in good agreement, and both of them can conclude that the RMS position error increases as the space debris altitude increases.
From Equation (30), it is not difficult to find that the triangulation positioning error not only depends on the observation accuracy but also depends on angles θ A and θ B between LOS vectors and the baseline of two EO sensors. It is noted that the RMS position error is infinite for θ A + θ B = 180 • , which occurs when the two LOS vectors are parallel. Therefore, ground-based triangulation is only suitable for positioning LEO space debris, but this limitation is avoidable in space-based triangulation.
From Figures 6-8, the RMS errors of SMA, eccentricity, and inclination are about 4.36 km, 0.0077, and 0.053 • , respectively. The IOD success rate is 99.4%. For single-site angles-only IOD, with the double-r method, the RMS errors of SMA, eccentricity, and inclination are about 550 km, 0.08, and 0.45 • , respectively, and the IOD success rate is only 35.1% [26]. Due to the constraint of range information, the IOD accuracy of triangulation is improved significantly. Compared with the single-site IOD, the two-site IOD is more stable and reliable.

Conclusions
In this paper, a new triangulation algorithm for positioning space debris is presented, and its analytical expression of RMS position errors is derived. We evaluate the new algorithm in comparison to traditional triangulation algorithms, and the simulation results not only validate our theoretical derivation but also show that the positioning results of the new triangulation algorithm are better than those of the traditional triangulation algorithms. IOD is carried out by using the triangulation positioning data, and the results show that the two-site IOD using triangulation positioning data is more stable and reliable than the single-site IOD using angles-only data. Triangulation is useful for positioning space debris and the IOD. However, to obtain precise positioning results, the configuration between the two observation sites and the observed space debris should be considered. Since IOD is the key to cataloging space debris, optical observation systems using triangulation would significantly contribute to cataloging space debris in the future.