{'fraction_non_alphanumeric': 0.032714325058032026, 'fraction_numerical': 0.011921155132391706, 'mean_word_length': 4.546039711979053, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Machine Learning-Based Approach for Improved Orbit Predictions of LEO Space Debris With Sparse Tracking Data From a Single Station

Accurate orbit prediction (OP) of space debris is vital in space situation awareness (SSA) related tasks, such as space collision warnings. However, owing to the sparse and low precision observations, unknown geometrical and physical features of debris, and effects of incomplete force models, OP based on the orbital mechanics theory or physics-based OP of space debris suffers from rapid error growth over a long duration, limiting the period of validity of debris OP for precise space applications. Considering that the tracking arcs of a debris object over a single station often share a similar temporal and spatial distribution in the inertial space, the resultant OP errors possibly have a coherent relationship with the temporal and spatial distribution of tracking arcs. This article proposes a machine learning (ML)-based approach to model the underlying pattern of debris OP errors from historical observations and apply it to modify the future physics-based OP results. The approach includes three steps: constructing a historical OP error set, training an ML model to fit the historical OP error set, and correcting the future physics-based OP with ML-predicted orbital errors. The ensemble learning algorithm of boosting tree is studied as the primary ML method for the error modeling and predicting process. Experiments with three low-Earth-orbit objects, tracked by a single radar station, demonstrate that the trained ML models can capture more than 80% of the underlying pattern of the historical OP errors. More importantly, the errors of physics-based OP over the future seven days reduce from thousands of meters to hundreds or even tens of meters through the error correction with the learned error pattern, achieving at least 50% accuracy improvement. Such dramatic OP improvements show the promising potential of ML for enhanced SSA capability.


I. INTRODUCTION
Space debris includes nonfunctional payloads, abandoned launch vehicle stages, mission-related debris, and fragmentation debris. There are more than 23 000 orbital objects larger than 10 cm in diameter tracked by the U.S. space surveillance network (SSN), of which space debris comprises more than 95% [1]. The population of space debris is rising rapidly, and the threats from space debris are reported frequently, highlighted by the well-known collision between operational Iridium 33 and defunct Cosmos 2251 on February 10, 2009 [2]. The requirements for safe and sustainable space operations drive the need to establish strong space situation awareness (SSA) and space traffic management capabilities in this increasingly crowded space environment. While detecting and tracking space debris requires very costly space-based and ground-based hardware facilities to gather data, how to accurately determine and predict the trajectories of space debris using the limited space surveillance facilities is the next biggest challenge. Accurate orbital knowledge of space debris can be widely applied in SSA-related tasks such as collision warnings [3], maneuvering target tracking [4], [5], debris detumbling [6], and the future active debris removal [7], [8].
The North American Aerospace Defence Command two-line element (TLE) catalog has been the sole orbital data source of space debris for public access via https:// www.space-track.org. However, the accuracy of TLE propagations through the Simplified General Perturbations-4 (SGP4) propagator decreases dramatically with time due to the simplifications made in the perturbation force models [9], [10]. For example, the seven-day maximum error of TLE propagations of Larets orbit (altitude 691 km) can grow up to 8 km, when comparing the TLE propagations to the high-precision ephemerides [11]. Hence, the TLE data are becoming gradually insufficient for the growing demands of space applications requiring higher accuracy, such as debris collision warnings, since missing such a warning could result in a catastrophic space event.
Alternatively, accurate orbital information of a space object can be generated from an orbit determination (OD) and prediction (OP) process with the orbit observations relating to the space object's positions and velocities [12]. In the OD process, the optimal estimation of the space object's truth orbit is obtained by adjusting the initial state parameters and force model parameters to fit the available observations best [13]. With the orbit state and force model parameters estimated from the OD process, the OP process begins by propagating the orbital solutions to the future time through the computation of the equations of motion. Theoretically, approaches to OP can be broadly categorized as numerical, analytical, and semianalytic methods [14]. The numerical methods provide accurate predictions of the positions and velocities of space objects by integrating the equations of motion with high fidelity force models, such as the Cowell integrator and the Runge-Kutta 4 integrator [13], but come with the disadvantage of slow computing speed. The analytical methods, such as the SGP4 propagator, can be used to solve the equations of motion efficiently, but the orbit solutions are not accurate as only simplified force models are employed. The semianalytic methods produce reasonably accurate but quickly delivered OPs, allowing the numerical integration of the equations of mean elements to proceed with long step sizes of half a day to one day, while the short-periodic terms are recovered analytically. Typical semianalytic methods include the averaging-based semianalytic integrator [15], [16] and the multiscaling-based semianalytic integrator [17], [18]. All of these OP approaches are essentially physics-based, which employ a series of physical models to quantify the effects of perturbation forces on the space objects' orbits. Consequently, the accuracy of any such OP process heavily depends on the OD accuracy and the fidelity of perturbation force models, such as Earth gravity, third-body gravitational attractions, solid Earth and ocean tides, atmospheric drag, solar radiation pressure, etc.
In theory, accurate OP can only be achievable when a good knowledge of the orbit state at the start of OP and the temporal and spatial variations of the space environment is available [13], [14]. For the satellites, their orbit states can be precisely determined with centimeter-level and globally distributed observations [19], [20] and then propagated in their high precision OP (HPOP) system. However, for the debris orbit problem, various sources could bring in large errors in the OD and OP processing. In an ideal case, the difference between the real observation and its theoretical values should be caused by the noise that occurs in making the observations. Yet, the sensors or instruments employed in the debris surveillance are not completely accurate, which inevitably introduces some discrepancies or errors in the debris observations. More importantly, our knowledge of the temporal and spatial variations of the space environment is limited, resulting in the dynamic model errors [21], namely, the differences between the nominal model parameters and the real model parameters, e.g., the atmospheric density model. The accuracy of widely used atmospheric density models, such as Jacchia71, MSIS86, and DTM94, is usually estimated at a 15% level [13]. As a result, the atmospheric drag, as the most significant nongravitational perturbation acting on low-Earth-orbit (LEO) objects, impacts the OP results seriously, especially for the long-duration prediction.
Although more complicated force models can be taken into account, the orbit estimation error is still an important factor affecting the accuracy of the subsequent OP, because of the erroneous debris observations. Radar and optical tracking are the two most popular tracking techniques in current space surveillance [22]- [25]. The accuracy of the radar tracking is typically at an order of tens of meters/arcseconds for the range/angle (azimuth and elevation) observations, according to ground-based radar station parameters in [14,. While the accuracy of the optical observations is usually at the level of 5-10 arcsec according to the Groundbased Electro-Optical Deep Space Surveillance [14], which is equivalent to 25-50 m at a distance of 1000 km between the station and the space objects. Additionally, due to the limited tracking facilities, the accurate OD and OP of space debris also confronts two main challenges [26], [27]: 1) the tracking data are mostly sparse and ill-distributed, leading to the OD computation difficult to converge, for instance, there are only one or two passes of tracking data observed per day; and 2) the ballistic coefficient or area-to-mass ratio needed for computing the nonconservative forces is unknown for most tracked debris.
To address the inaccurate debris OP challenges, multiple research efforts have been made through the modification of classical physical models. The major one is to improve current atmospheric density models through calibration methods [28]- [30] or develop new models [31]- [33]. As for the problem of unknown ballistic coefficients, research efforts have been made to use long-term historical TLE sets to estimation them [5], [34]. Besides, to overcome the geometrical weakness of sparsely distributed tracking arcs, the TLE-predicted positions are used as pseudoobservations to constrain the OD process so that the OP accuracy improves [26], [27]. However, these methods are still far from satisfactory for space safety operations to some extent, since developing high-precision physical models is a very challenging task that requires high-quality and long-term observations with sufficient temporal and spatial resolutions.
Considering a single-station observation scenario, a piece of space debris may be tracked once or twice each day. The tracking arcs from the single station appear stable in terms of the temporal and spatial distribution in the inertial space. Thus, it is likely that the OP errors using the single-station observations follow some underlying pattern if the space environment does not change much. The longterm space surveillance has accumulated huge amounts of debris tracking data, which makes it possible to study the underlying pattern of debris OP errors in the short term and use it to improve the future OP accuracy. Fortunately, with the development of machine learning (ML) technique [35], [36], through which a mathematical model is automatically learned from past experiences without being explicitly programmed to perform the task, it provides a novel solution for the pattern recognition of debris OP errors.
Unlike the physics-based OP method, the ML method is essentially a data-driven modeling and predicting approach. It learns a function or mapping from labeled samples, known as the training data, in order to make predictions or decisions for a new instance. An ML model, which is trained by the ML algorithm, is a mathematical expression Y =f (X ) that represents the relationship between the input X and the output Y in the context of a problem. Various ML algorithms and tools have been developed and used in a wide variety of space applications, such as automatic space objects characterization [37], planning and optimization of spacecraft trajectories [38], [39], and orbit determination [40], [41]. These applications clearly show the excellent capability of ML in the orbit dynamics field, which prompts researchers to a keen interest of applying the ML method for accurate OPs of resident space objects. Peng and Bai studied the support vector machine (SVM), Gaussian Processes, and artificial neural network to improve the OP accuracy [42]- [44], together with the performance comparisons of these three ML algorithms on the OP accuracy improvement [45].
In this article, we propose a hybrid method to achieve accurate debris OP results with single-station sparse tracking data through integrating the ML method into the physicsbased OP system. With the evolution pattern of historical OP errors learned by the ML method, it is used to produce orbital error predictions for the future physics-based OPs, improving their accuracy via the error compensation or correction. A complete ML framework for improved physics-based OP results is implemented in three steps, including constructing the historical OP error set in a sliding OD/OP window using sparse tracking data, training an error model with the supervised ML algorithms, and correcting the future physics-based OPs using ML-predicted orbital errors. In the developed ML framework, the learned ML model can be treated as an error corrector to the physicsbased OP system. Thus, accurate physics-based OPs of space debris can be achieved over a long duration, if the learned ML model captures the majority of the error pattern.
The rest of this article starts with the description of debris OD and OP process using sparse tracking data in Section II. Next, the framework of the ML approach in modeling the debris OP error pattern is provided. In Section IV, the experimental results of the ML approach in improving the physics-based OP accuracy are presented. Finally, Section V concludes this article.

II. DEBRIS OD AND OP USING SPARSE TRACKING DATA
A. Single-Station Radar Tracking Scenario As an important part of the SSN, the ground-based radar system is best suited for space object tracking in the LEO region, and some very high powered radars are also capable of tracking objects beyond the LEO region [1], [14], [23]. Considering a single ground radar station and one debris, a set of radar tracking data includes the pointing angle and two-way range measurement of the debris, where the ground radar station serves as the signal transmitter and receiver. Essentially, the radar pointing angle is obtained in the topocentric coordinate system of the station by measuring the direction of the maximum signal amplitude of a space object, while the radar ranging measurement is the distance computed from the two-way flight time of a radar signal fired from the station to the space object and radiated back to the station. Thus, given the radar station location and debris position in the Earth Centered Earth Fixed coordinate system, the radar observations, including the azimuth/elevation O Az/El and the distance O Dis of a debris object with respect to (WRT) the radar station, can be derived according to the observational model, which can be referred to [13,Sec. 6.2].
As the purpose of this article is to develop an ML framework for improving the debris OP accuracy under the condition of sparse tracking data, only a single ground radar station is considered. To reflect the realistic tracking scenario, strict configurations are assumed for the radar station by referring to the radar system parameters given in [14,, as follows: 1) the debris must be visible to the ground radar station; 2) the accuracy of radar angular/ranging observation is 100 arcsec/100 m, respectively; 3) the maximum two-way ranging distance is 5000 km; 4) the minimum elevation mask is 20 • .
In the experiments, the ground radar station is located in Wuhan, China, with the geodetic latitude of 30.542 • , longitude of 114.346 • , and altitude of 39.3 m. The laser ranging satellites in the LEO region are selected as representatives in the radar tracking scenario. Their precise orbits in the consolidated prediction format (CPF) [46] are downloaded from the International Laser Ranging Service (ILRS) data center 1 and used to simulate the radar observations. A typical CPF file usually contains 3-5 days of accurate satellite positions and these positions are generated by propagating the precisely determined satellite state in its HPOP system. In this article, only the satellite positions from the first day in each CPF file are used, as they are more accurate than those of later days. Then, a precise CPF orbit can be produced by combining the first-day positions from a series of CPF files. Finally, a precise OD process with high fidelity force models is implemented to fit the combined positions, in order to avoid the discontinuities that occur in making the precise CPF orbit. Following these procedures, a precise and smooth CPF orbit over long arcs is obtained, which will be used as the truth orbit to compute the radar observations according to the preset configurations.

B. Orbit Estimation Strategy
The equations of motion of an Earth-orbiting object in the Earth-centered inertial (ECI) coordinate system is given by a system of three second-order differential equations of [14] f (x, c, t ) =r = −μ r where r,ṙ, andr are the space object's position, velocity, and accelerate vectors in the ECI frame, respectively; r = ||r||, with || · || denoting the Euclidean norm of a vector; μ = GM is the Earth gravitational constant; a per is the perturbation acceleration vector caused by nonspherical gravity, thirdbody gravity attraction, atmospheric drag, solar radiation pressure, and so on; Γ is the thrust acceleration vector; f (x, c, t ) denotes the equations of motion, in which x is the space object's state (six-dimensional), x = [r,ṙ]; c is the force model parameter, which refers to the drag coefficient C d used as a critical parameter to compute the atmospheric drag effect; t is the time. For a LEO object, given the initial condition of (1), i.e., the initial state x 0 = x(t 0 ) and the drag coefficient C d , x t at a future time epoch t can be uniquely predicted with a numerical, analytical, or semianalytical orbit propagator. In practice, the true value of x 0 is never exactly known, but can be accurately estimated using a batch or filtering algorithm. Assume that n radar angular observations O Az/El,i , i = 1, 2, . . . , n, and m radar distance observations O Dis, j , j = 1, 2, . . . , m of a tracked object are available, with their corresponding weights of P Az/El,i and P Dis, j characterizing their relative accuracy. Let v Az/El,i , v Dis, j denote the residuals between real observations O Az/El,i , O Dis, j and their theoretical values. Given the initial condition ( where C Az/El,i and C Dis, j correspond, respectively, to the theoretical values of O Az/El,i and O Dis, j , which are computed by the space object's state prediction and the radar station location according to the radar angle and range observational models. After obtaining all the residuals of available observations, a weighted sum of squares of the residuals (WSSR) is given by Obviously, WSSR is a function of the initial condition. In the sense of batch least squares (BLS), the optimal estimation of the initial condition can be obtained by minimizing WSSR WRT the initial state and force model parameter, asx where δx 0 and δC d are the estimated corrections to the initial condition of x 0 and C d , respectively; B is the error equation matrix; P is the weight matrix of the observations, which is typically a diagonal matrix with its elements being the inverse of the observation variance; l is the observation residual vector; H t,t is the observation matrix, H t,t = ∂H ∂x t ; Fig. 1. Illustration of the tracking arc (orbit observations of one pass of space debris WRT radar station) availability for each OD process. and t,t 0 is the state transition matrix, containing the variational partial of x t WRT the initial condition, In our experiments, an orbit determination and analysis software for processing the angular and ranging observations is developed by Li et al. [47], [48]. It applies the Cowell numerical orbit propagation method for the state solutions in (1) and the BLS algorithm for the optimal estimations in (5). In the BLS OD process, three days' radar tracking arcs of a debris object are used to compute the optimal estimation (x 0 ,Ĉ d ), and the OD fit span is set starting at the epoch time of the first tracking arc in radar tracking scenario. The initial orbit state in the BLS OD system is computed from the latest TLE before the OD fit span. As two types of observations, O Az/El,i (two-dimensional) and O Dis, j (one-dimensional), are weighted in the BLS OD process, they provide different geometrical strengths when forming the normal equations (6). Thus, their relative weights should be properly assigned for obtaining the optimal estimation, regarding to the observation accuracy and the geometrical strengths. In our BLS OD process, a weight strategy of P Az/El = 1.0, P Dis = 10 −10 is adopted by referring to the citation [27], so that the accuracy of O Dis (in meters) is at the same order as that of O Az/El (in radians) when fitting the angle and range observations together (three-dimensional observations).
The BLS OD computation is an iterative process and converges when changes between the estimated initial state and its previous estimate are less than the preset values of 1.0 m for the position vectors and 0.01 m/s for the velocity vectors. After the OD convergence, an outlier detection procedure is performed to ensure no observation with gross error is used in the OD process. After the completion of OD computation, the OD time window will slide to the next one with the starting time epoch set at the beginning time of the next tracking arc, as illustrated in Fig. 1. Following this OD procedure, all the pairs of the estimations (x,Ĉ d ) of the available tracking arcs are obtained.
To reflect the orbit error pattern relative to the high fidelity dynamical system of CPF orbits in simulating the radar tracking scenario, a simplified dynamical system (SDS) is adopted in (1) for the debris OD and OP processing so that any variations of space environment could introduce errors in the OD and OP processing. In the SDS, the main perturbations acting on LEO objects are considered, including Earth gravity, Sun and Moon attraction, solar radiation pressure, and atmospheric drag. The Joint Gravity Model-3, truncated to 20/20 degrees/orders, is used to compute the Earth gravity effect. The Sun and Moon positions are provided by the DE406 planetary ephemeris. When computing the solar radiation pressure and drag effects, the area-tomass ratio of a space object is determined as BC/2.2, where BC is the TLE-derived ballistic coefficient [34]. Initial values of the force model parameters, including the drag coefficient C d and the solar radiation pressure coefficient C r , are set with C d = 2.2 and C r = 1.1. The NRLMSISE-00 atmospheric density model [49] is applied to compute the atmospheric mass density ρ, which will be then added into a random error from a Gaussian distribution with zero mean and 0.1ρ standard deviation, in order to reflect the effect of atmospheric uncertainty.

C. Construction of Historical OP Error Set
After obtaining the estimation (x,Ĉ d ) of each tracking arc following the above OD strategy, an orbit integration of (1) with (x,Ĉ d ) as initial values is performed in the SDS through the Cowell integrator, yielding the predicted state as wherex i is the estimated state of the i-th tracking arc, with its initial time at t i ;x j;i is the predicted state at a future time t j based onx i , and t j − t i is the prediction time length; φ is the implicit state solution of (1), i.e., state propagation [13], [21], [48].
In (8), the error of the predicted statex j;i WRT the truth orbit can be interpreted as a result of the mutual effect of the orbit estimation error and the dynamic model error. As the truth orbit of space debris is unknown in the real data processing, the estimationx j is used instead as the reference to assess the accuracy ofx j;i , expressed as where e j;i is the error ofx j;i , i.e., the OP error;x j is the state estimation at t j ; t j > t i + 3 days, as the OP starts from the terminal time of the three-day OD window. In (9), e j;i is computed in the ECI coordinate system, which will be then transformed into the RSW coordinate system [14] for orbital error pattern analysis, where R axis points from the Earth's center along the radius vector towards the space object, S axis points in the direction of the velocity vector and is perpendicular to the radius vector, and W axis is normal to the orbital plane. In this way, an OP error set in the sliding OD window is constructed by combining e j;i of all the pairs of state estimations, as illustrated in Fig. 2.

III. ML IN MODELING THE OP ERROR PATTERN
Due to the observation error and incomplete perturbation force models, the estimated orbit states of space debris still contain errors. This causes the subsequently predicted orbit to further deviate from its truth orbit as a result of the effect of dynamic model error. According to the orbital error propagation characterization [21], [48], there would be some patterns existed in the OP errors. Inspired by the ML theory, such patterns can be modeled through the data analysis based on large amounts of OP error samples. Essentially, the ML method learns a mathematical model Y =f (X ) between the input X and the output Y of the labeled samples and then makes predictions when a new input is given. In this section, we will first discuss the data preparation for learning an error model Y =f (X ) using the ML method, followed by the implementation of ML algorithms in modeling the evolution pattern of historical OP error set.

A. Data Preparation
The ML method works on the sample data (training data and testing data), which are independent and identically distributed, even if their probability distribution is not known beforehand. Hence, if a pattern is contained in these samples, it can be learned by the ML algorithm. For a specific problem, the success of the pattern learning strongly relies on data preparation.
Generally, the training data are composed of the input X ∈ R N and the output Y ∈ R. The purpose of a learning process is to find the underlying mapping relationship between X and Y if it exists, i.e., Y =f (X ). Therefore, the training data have to be well-prepared such that this underlying mapping relationship can be captured. It is always easy to determine the target variable Y . On the contrary, the choice of the learning variables in constructing X is not straightforward, as there are many latent features associated with Y . Theoretically, the more learning variables are included in X , the more accurate but complex the trained model would be. However, it would result in the overfitting problem [50], in which the trained model works well for the training data, but produce bad predictions for the unseen data. In this sense, a desirable model should be trained to avoid overfitting, but at the same time take into account the most important factors affecting the target variable Y .
As expressed in (1), (8), and (9), the OP error e j;i is computed by comparing the predicted statex j;i to the reference state, whilex j;i is obtained by numerically integrating the equations of motion with the given initial condition (x i ,Ĉ d,i ). In this process, the factors that directly contribute to e j;i should be included to construct X . Thus, at a specific epoch t j , X includes the learning variables as follows. 1) OP time length t j;i = t j − t i , t j;i ≤ 7 days, because the accuracy of state propagation accuracy decreases as t j;i increases. The longer t j;i is, the larger the magnitude of e j;i will be. Obviously, the target variable Y corresponding to X refers to the OP error as follows.
1) e j;i at time t j , expressed in the RSW coordinate system as e = [e A , e C , e R , e vA , e vC , e vR ] T , where e A , e C , and e R denote the along-track, cross-track, and radial position errors, respectively, and e vA , e vC , and e vR denote the along-track, cross-track, and radial velocity errors, respectively. It means that six error models need to be trained for the position and velocity components.
Given the learning and target variables defined above, a sample of the training data is constructed in the form of ( and Y ∈ e j;i . It means that each target variable corresponds to 26 learning variables in the training data so that the underlying error pattern can be sufficiently modeled. Hence, with a series of samples (X , Y ) being well-prepared, an underlying mapping relationship Y =f (X ) can be trained by the ML approach in the learning system. More importantly, if the testing data beyond the learning period is correlated to the training data, the trained ML model can be applied to the testing data to predict new errors, which are then used to correct the state predictions during the testing period, as illustrated in Fig. 3.
In Fig. 3, the ML-predicted orbital error, denoted as e ML , is computed by e ML =Ȳ =f (X ) (10) whereȲ is the error prediction corresponding to X of the training data or testing data with the trained ML model. Once e ML is obtained, the error correction is performed to modify the physics-based state predictionsx j;i by removing e ML , asx ML j;i =x j;i − e ECI ML (11) wherex ML j;i is the improved physics-based state predictions; and e ECI ML is the e ML term expressed in the ECI frame. The residual error e corresponding tox j;i after error correction is equivalent to where e is the original error ofx j;i . When assessing the performance of learned ML models in modifyingx j;i via error correction, the improvement percentage is quantified by the performance metric P ML as where rms e = 1 n n i=0 e 2 i is the root mean square (rms) value of e; rms e = 1 n n i=0 e 2 i is the rms value of e; and n is the number of training or testing samples. In (13), all the samples are included in computing P ML , giving an overall performance assessment of the ML model on the training or testing data. The higher the P ML value is, the better the trained ML model is obtained. In an ideal case, P ML can reach 100%, indicating all the underlying patterns contained in e are perfectly extracted by the trained ML model, namely e ML = e. However, P ML could have a negative value if the OP error pattern cannot be found by the trained ML model.
The learned error pattern is strongly dependent on the configuration of the learning variables in the sense that whether the learning variables can characterize the features that affect OP errors, as well as the temporal and spatial distribution of the training and testing data. Thus, the trained ML model can be applicable to the time periods if they share a similar space environment. To guarantee the period of validity of the trained ML model for future error prediction, in the experiments, the time periods of the training data and the testing data are formed by a sliding time window, as shown in Fig. 4.
As illustrated in Fig. 4, the training data over the past 30 days is used to train an ML model. The trained model is then used to predict the OP errors over the seven days following the training period. In this case, only the tracking data in the three days before the last tracking arc, as well as the learning and target variables in the past 30 days, is needed. The radar tracking data are processed in the standard physicsbased OD/OP process, and the OP errors are predicted by the trained ML model with the learning variables derived from the OD/OP results. The ML-predicted errors are then applied to correct the physics-based OPs.

B. ML Algorithm for Training the OP Error Model
Modeling the relationship between the input X and the output Y of the prepared OP error set is a multivariable regression problem. Of all the types of ML methods, the supervised ML methods appear suitable to learn a mapping from X to Y . In the work of [42], the SVM [51], [52] was studied to improve the satellite OP accuracy, but the errors of improved position predictions were still at an order of tens of kilometers, although three radar stations were employed in the satellite observation. Such an OP accuracy cannot be applicable to the precise space applications. In this paper, we apply the ensemble learning method based on the boosting algorithm instead to mine the underlying OP error pattern of space debris under a single-station tracking scenario. In this section, the principle of the boosting tree (BT) algorithm is outlined, and more details are referred to the citations.
The BT algorithm is a popular ensemble learning method that can be applied for both classification and regression problems [53]. It is a boosting method with the decision tree as a base learner, in which a strong learner is learned by combining all the predictions of multiple weak learners (base learners) through the additive model. One significant advantage of the BT algorithm is that each base learner focuses its effort on the most difficult samples to fit, generating an ensemble model with lower bias at the end of training process.
A BT model through the AdaBoost algorithm [54] is expressed in the form of where x is the input variable, and x = X in our experiments; T (x; m ) is the decision tree (base learner) with its parameter of m ; and M is the number of T (x; m ). Here, T (x; m ) refers to the regression tree, since a regression task is performed in our study. Given a set of training data, the BT model is trained using the forward stage-wise algorithm as where N is the number of the training samples, y i is the output variable corresponding to the input variable x i of the training data, and y i = Y i in this article according to its definition. The quadratic loss function L(y, f (x)) = [y − f (x)] 2 is used to determineˆ m in (16) by minimizing the loss of (17) where r = y − f m−1 (x) denotes the residual at step m − 1.
In (14)- (17), the BT algorithm is interpreted that the first learner T (x; 1 ) is used to fit the whole space of input data, and then the second learner T (x; 2 ) is applied to fit the residuals of the first fitting process, thus overcoming the drawback of the first learner. This fitting process is repeated until meeting the stopping criterion. Finally, a BT model, as expressed in (14), is obtained to approximate the target function Y =f (X ) by the sum of the predictions of each base learner.
The implementation of the BT algorithm is summarized in Algorithm 1.
In this article, the BT algorithm for training the OP error model is programmed in the MATLAB (R2018a) environment. In the implementation of the BT algorithm, the least squares boosting is adopted for solving the regression problem [55]. A grid search method is adopted to adjust the parameter values in the BT model. In the experiments, the size of the OP error set is very limited, less than 900 samples, which will be given in Section IV-B. In this case, an amount of 30 base learners, with the leaf size being minimum 8 and learning rate being 0.1, are determined to fit them in the learning process. We note that these settings are not unique. It is a trial-and-error process, as the larger the leaf size and/or the number of employed learners is, the slower the training speed will be.

IV. EXPERIMENTAL RESULTS
Three LEO space objects, Larets (altitude 691 km), Starlette (altitude 815-1100 km), and HY-2A (altitude 971 km), are tested in the experiments. Their precise CPF-format orbits are available in the ILRS data center [56]. In modeling the OP errors of the satellites, the BT algorithm is used as main ML algorithm for the training tasks. In addition, the performance of ML based on the support vector regression (SVR) in the error pattern learning is also presented, in order to validate our developed ML framework for improved OPs of LEO space debris through modeling and transferring the historical orbit error patterns in the single-station radar tracking scenario.

A. Accuracy of OD and OP Using Sparse Tracking Arcs
In the experiments, a 90-day single-station radar tracking scenario is simulated for the test satellites, and up to two tracking arcs are picked each day to determine the orbits. In this radar tracking scenario, a total of 150 radar tracking arcs are collected by the radar station for each satellite. The average time length of each tracking arc is 305.7 s for Larets orbit, 483.6 s for Starlette orbit, and 502.0 s for HY-2A orbit, which are very short when comparing to their orbital periods.
According to the orbit estimation strategy in Section II, the BLS OD processing using tracking arcs in a three-day fit span is performed to obtain the optimal state estimationx for the three satellites. Asx will be used as the reference to compute the subsequent seven-day OP errors, an accuracy assessment ofx is necessary to verify thatx is sufficiently precise as the approximation of the truth orbit (precise CPF-format orbit). The true estimation error, denoted as , is computed as the difference betweenx and the truth orbit, and then transformed into the RSW coordinate system. It indicates the correctness of the developed BLS OD processing using sparse radar observations. Considering the magnitude of seven-day OP errors at such selected orbit altitudes is usually at hundreds or even thousands of meters, the true estimation errors hardly affect the OP error pattern learning. Hence, the state estimationsx can be regarded as accurate enough in representing the truth orbit to compute the OP errors. This is important for the OP error computation of space debris, as only the state estimationsx are available as references for the orbital accuracy assessments.
With the state estimationsx as references, a 30-day OP error set is constructed for the learning process. Fig. 6 presents an example of the OP error set over the first 30 days (1st-30th day) of the Larets orbit.
In Fig. 6, the along-track position prediction error e A grows dramatically as the OP time length increases, with the maximum value of the whole OP errors approximating 7 km over the 7-day OP time span, while the cross-track and radial position prediction errors, e C and e R , scatter around the zero value and are all within [−50 m, 50 m]. Obviously, e A contributes to the most of the position prediction error, indicating the significance of reducing e A to improve the OP accuracy. Therefore, in the following learning process, e A will be used as the target variable to be learned by the ML algorithms.

B. OP Error Pattern Learning Using Training Data
The purpose of ML is to learn a desirable OP error model Y =f (X ), which not only fits well with the training data, but also makes good predictions with new inputs. To avoid the overfitting problem and generalize the learned OP error pattern to future epochs, the cross-validation (CV) is necessary for the learning process. It is a statistical method used to obtain an accurate estimate of the ML model, which generally results in a low bias [57].
A common CV is the holdout CV that splits the original dataset into two parts (training and testing) and uses the testing score as a generalization measure [36]. It usually works well when the original dataset is large. However, due to the small size (less than 900 samples) of the OP error set, the k-fold CV [58] is more suitable for the error pattern learning in this study, as it generates a less biased estimate of the ML model than the holdout CV methods. In the k-fold CV process, the original OP error set is randomly partitioned into k disjoint sets with equal-sized subsamples. Of the k subsamples, a model is trained using the k − 1 subsamples (training data) while the remaining single subsample (validation or test data) is used to assess the performance of the trained model. Then, the CV is repeated for k times, with each of the k subsamples being used exactly once as the validation data. The performance metric of the trained model is finally computed as the average test error over all folds.
With e A being the target variable to be modeled, an MLpredicted error e ML,A can be computed by the learned error model Y =f (X ) and then used to correct e A . Figs. 7-9 present the ML results of the BT and SVR methods on the training data for the test orbits. In each figure, the lower scatter plot shows the original error e A (blue cross) and ML-predicted error e ML,A (red dot), and the upper scatter plot presents the residual error e A (green cross) after the error correction. Here, e ML,A refers to BT-predicted error e BT,A or SVR-predicted error e SVR,A .
In Figs. 7-9, it is clear that the ML-predicted error cluster e ML,A and the original error cluster e A are located closely, with their difference e A being evenly distributed around the zero value within a boundary. Taking Larets orbit, for example, the worst e A is −6722.3 m. However, after the error correction with the trained BT/SVR models, the worst residual error e A reduces to −622.7 m/−671.8 m through removing e BT,A /e SVR,A from e A , respectively. A decrease of 90.7%/90.0% is obtained for the worst error by the trained BT/SVR models. To demonstrate the ML performance on the whole in Figs. 7-9, the P ML values of the trained BT and SVR models on the training data are computed and used as the main metric to assess the performance of the trained BT and SVR models, as given in Table I.
In Table I, the performance metric P ML is computed by rms e A and rms e A according to its definition. It can be seen that rms e A is significantly smaller than rms e A and the computed P ML value exceed 81% for all the trained BT and SVR models. Moreover, according to the metrics, the trained BT models show stronger learning capability than the trained SVR models on the same training data. Although such superiority of the trained BT models is small, it demonstrates the proposed BT method has the ability to train an ML model with lower bias. All in all, the developed ML framework is capable of capturing the dominating components of the underlying error pattern when only 30 days' training data are available. If the learned pattern agrees well with the realistic patterns, it can predict accurate orbital errors beyond the training period, which will be evaluated in the following section.     From the application perspective, one is more concerned with the performance of the trained ML model with new data that are not used in the model learning. In this case, a desirable ML model should be learned with good generalization capability, that is, producing accurate predictions for the unseen data. In terms of the OP accuracy improvement, after the trained error model is obtained from the historical OP error set, it can be generalized to future epochs to predict new error terms corresponding to the latest OPs. Thus, the improved OPs are achieved via error correction. Now, the BT and SVR models, learned by the first 30 days' training data in Figs. 7-9, are applied to predict the future seven days' orbital errors. As illustrated in Fig. 4, the physics-based OPs during this period are used as the testing data to input the models.Figs. 10-12 plot the ML-predicted error e ML,A , the original error e A , and the difference between e ML,A and e A , over the future seven days (31st-37th day) for these three objects. In each figure, the generalization capability of the trained models can be seen from the performance metrics P ML .
We observed that in Figs. 10-12, the error predicted by the trained BT and SVR models with the testing data follows the same growing trend of the original error over the time span of 31st-37th day. After correcting the predicted state with the ML-predicted error, the maximum absolute error of Larets/Starlette/HY-2A over seven days reduces from 2186.3 m/548.9 m/396.7 m to 465.7 m/154.4 m/103.5 m for the trained BT model and 707.0 m/145.5 m/145.6 m for the trained SVR model, respectively. The performance metric P ML of the trained BT/SVR model on the testing data reaches 79.7%/75.1% for Larets, 80.9%/78.5% for Starlette, and 75.2%/60.8% for HY-2A. It means that the ML-predicted errors compensate the majority of the real errors of the physics-based OPs during the future seven days. More importantly, it proves the effectiveness of the learned error pattern. Besides, it rationalizes that performances of the learned error models on the testing data are slightly worse than that on the training data, since these two datasets are from different time spans. Hence, the error pattern learned from the past period may be not exactly the same as that in the future period, because of the variations of the orbital environment such as the atmospheric density. Fortunately, the orbital environment does not change much within short terms, indicating that the orbital error evolutions during these two periods share the common pattern in most cases. As a result, the learned error pattern can predict the majority of the real OP errors in the future seven days, as revealed in Figs. 10-12. To further illustrate the effectiveness of the learned error models in predicting the future OP errors, the ML training and testing processings in a sliding time window of 50 days are performed. In this sliding time window, 100 combinations of the training and testing data are constructed, as illustrated in Fig. 4. That is, the combination of tracking arcs in one ML training and testing is different from other ML processings. For each dataset combination, an ML model is trained with the training data over the previous 30 days' time span, and applied to the testing data in the following seven days. Fig. 13 plots the rms values of the future seven days' OP errors before and after the error correction with ML-predicted errors for the 100 ML runs for Larets orbit.
In Fig. 13, the generalization capability of the trained ML models on the future testing data are well demonstrated by the 100 ML runs, where the original OP error e A is reduced dramatically using the trained BT and SVR models, especially for cases with large rms e A . However, for cases that the OP results are reasonably accurate over the first one or two days of predictions, the trained models show limited accuracy improvement, and sometimes it undermines the accurate OP results.   Fig. 13, the original OP errors of these five combinations are all at the level of tens of meters. As most of the OP errors in the training data are large, they have a dominant effect on the error pattern learning. Therefore, it is expected that most of the large OP errors will be better modeled, and the few small errors will be modeled less accurately. On the other hand, a minor accuracy degradation in the originally accurate OP results is a worthy sacrifice for the overall OP accuracy improvement.

D. Effect of Training Data Size on the Performance Metric
More training data usually contributes to a better ML model. This is conditional on the assumption that the whole dataset follows the same distribution or shares the common underlying pattern. In practice, additional data could hardly improve the ML performance if it brings uncorrelated information or its related space weather condition in the training period and testing period changes heavily. To investigate the effect of training data size on the performance metric P ML , a series of BT and SVR models are trained with different sizes of training data and then evaluated on the same testing data. Fig. 14 gives an example of Larets orbit for the trained BT model using a training data size of 10, 20, 30, 40, 50, and 60 days. In Fig. 14, e denotes the original OP error e A of the testing data; e BT|10days denotes the predicted error e BT,A using the BT model trained by 10 days' training data, and so on for e BT|20days , . . ., e BT|60days .
As revealed in Fig. 14, the effect of the training data size on the performance of the trained ML model is demonstrated by using 60, 50, 40, 30, 20, and 10 days' training data, with their corresponding performance metric P ML being 90.5%, 87.9%, 86.0%, 84.3%, 70.1%, and 52.6%, respectively. To demonstrate more clear variations of P ML , the training data size is varied from 10 to 60 days with a step size of 5 days. A total of 30 dataset combinations of training data and testing data are processed. Fig. 15 presents the average P ML values of the trained BT and SVR models with different training data sizes for Larets orbit.
In Fig. 15, P ML increases quickly when the data size increases from 10 to 25 days for both the trained BT and SVR models, but tends to stabilize once the size of the training data reaches and exceeds 30 days. This indicates the performance of the trained BT and SVR models will not be further improved when a sufficient size of training data is used in the training process. It also rationalizes that our choice of using 30 days' training data in the previous experiments is sufficient in achieving the expected ML performance.

V. CONCLUSION
Providing accurate OP information of space debris is vital in enhancing the SSA capabilities, such as the reliable space conjunction warning. This article proposes a hybrid method of integrating the ML method into the physics-based OP system for accurate OPs of space debris when only sparse tracking data from a single radar station is available. This method displays a good generalization capability that the trained ML model is able to predict the future physics-based OP errors, which can be used to modify the physics-based OP results.
This article has resulted in a complete computational framework that embeds the ML method into the debris BLS OD/OP system. Two ML algorithms, including the BT method and the SVR method, are implemented to model the underlying pattern of past OP errors. Experimental results of three LEO satellites have demonstrated that both the trained BT and SVR models are capable of capturing most features of the underlying error pattern. The predicted errors by the trained BT and SVR models agree well with their truth values. By applying the ML-predicted OP errors to correct the physics-based OPs, the 7-day OP errors reduce from hundreds or even thousands of meters to only tens of meters through the error correction with the learned error models, achieving at least 50% accuracy improvement.
The current research has assumed the sparseness of tracking data from a single station. In the next phase of the research, datasets from a regional network of several stations in China will be processed, in order to assess the achievable OP accuracy of tracked space objects through the developed ML computing framework. As part of the orbital knowledge, the orbital uncertainty or error covariance knowledge plays an important role in the precise space applications such as space collision probability computation and data association. The application of ML method to the uncertainty propagation problem of space debris will be also studied in future work, thus providing the user groups with whole orbital information.