Re-entry prediction of objects with low-eccentricity orbits based on mean ballistic coefficients

Abstract During re-entry objects with low-eccentricity orbits traverse a large portion of the dense atmospheric region almost every orbital revolution. Their perigee decays slowly, but the apogee decays rapidly. Because ballistic coefficients change with altitude, re-entry predictions of objects in low-eccentricity orbits are more difficult than objects in nearly circular orbits. Problems in orbit determination, such as large residuals and non-convergence, arise for this class of objects, especially in the case of sparse observations. In addition, it might be difficult to select suitable initial ballistic coefficient for re-entry prediction. We present a new re-entry prediction method based on mean ballistic coefficients for objects with low-eccentricity orbits. The mean ballistic coefficient reflects the average effect of atmospheric drag during one orbital revolution, and the coefficient is estimated using a semi-numerical method with a step size of one period. The method is tested using Iridium-52 which uses sparse observations as the data source, and ten other objects with low-eccentricity orbits which use TLEs as the data source. We also discuss the performance of the mean ballistic coefficient when used in the evolution of drag characteristics and orbit propagation. The results show that the mean ballistic coefficient is ideal for re-entry prediction and orbit propagation of objects with low-eccentricity orbits.


Introduction
Massive space objects cannot completely burn up during a re-entry; 10% to 40% of the mass may survive (Ailor et al. 2005), and the surviving components may pose a threat to humans, buildings, and the environment (Choi et al. 2017). Usually we call these objects risky re-entry objects. To mitigate these risks we need predictions of re-entry of massive space objects. Atmospheric drag is the dominant non-gravitational perturbation acting on re-entering objects. Precise area-to-mass ratio of the re-entry object and proper modeling of the drag characteristic are the key issues in calculating the acceleration due to atmospheric drag and predicting the re-entry time. The drag coefficient depends on many parameters, such as the object's shape, the surface material, and the composition and the temperature of the atmosphere. An object's drag coefficient may vary under different altitudes and solar activity conditions (Moe and Moe 2005). Usually the space object's exact shape, mass, attitude, and surface material are unknown; therefore, it is difficult to estimate the drag coefficient, cross-sectional area, and mass. Hence, it is reasonable to combine these parameters into the ballistic coefficient B to scale how much the space object suffers from atmospheric drag (Bowman 2002). B is defined by Equation (1): where C D , A and m are the drag coefficient, cross-sectional area, and mass, respectively. Atmospheric drag decelerates the re-entry object, thus causing a gradual decrease in eccentricity and semi-major axis. In this paper we classify re-entry objects according to their initial apogee altitude, H a_ini , as follows (here initial means ten days before re-entry, which is also the limitation to the perigee altitude): 2. Low-eccentricity orbit: 500 km ≤ H a_ini < 5000 km; 3. High-eccentricity orbit: H a_ini ≥ 5000 km.
According to this classification, the statistics for large uncontrolled objects that have re-entered from 2012 to 2018 are shown in Figure 1. The information of the re-entered objects is obtained from the space-track website. We see that objects with low-eccentricity orbits are not in the minority, exceeding one-sixth in the past seven years, and even 40% in 2018. The main difficulties in re-entry prediction are determining the orbit and modeling the atmosphere drag . The optimal method for re-entry prediction is to determine the precise orbit based on abundant observations, estimate the ballistic coefficient as a fitting parameter, and then propagate the orbit until re-entry. Xiong et al. carry out the re-entry predictions by mainly employing a numerical integration method based on observations acquired from the optical surveillance network. They obtain the same re-entry time from two precise orbits by adjusting the ballistic coefficient to improve the prediction accuracy, and they propose that the change in the drag coefficient with altitude should be considered when the mean altitude is below 180 km (Xiong et al. 2009).
Only a few countries, however, have the ability to acquire observations of the re-entry objects. Most re-entry predictions are performed based on TLE (Two Line Element) sets provided by United States Strategic Command. Pardini et al. used Satellite Re-entry Analysis Program (SATRAP) software, which uses TLEs as the input, to study the prediction errors under different solar activities and different atmospheric model conditions Anselmo 2001, 2008;Anselmo and Pardini 2013). Different methods have been developed in the earlier studies to improve TLE-based re-entry prediction. Saunders et al. estimated the ballistic coefficient by comparing the change in the semi-major axis according to TLE data with the change in semi-major axis due to the drag computed by propagation using an initial state from TLEs (Saunders et al. 2012;Sang et al. 2013;Mutyalarao and Sharma 2011;Lu and Hu 2017). Agueda et al. derived pseudo-observations from TLEs and then estimated the ballistic coefficient, solar radiation pressure coefficient, state vector, or a combination of these (Agueda et al. 2013;Lidtke et al. 2016). Gondelach et al. studied the re-entry predictions of GTO rocket bodies, and proposed that TLEs should be preprocessed to filter out outliers (Gondelach et al. 2017;Lemmens and Krag 2014).
When compared with objects in nearly circular orbits, we find that re-entry predictions of objects in loweccentricity orbits are more difficult. Because objects with low-eccentricity orbits traverse a large portion of the dense atmospheric region almost every orbital revolution, the drag coefficients change significantly in one revolution. In a practical test it was found that it might be very difficult to obtain small residuals for each pass, particularly if we adopted the strategy that estimated a single ballistic coefficient from multiple passes of observations. If we adopted the strategy that estimated multiple ballistic coefficients from observations, it might be difficult to select a suitable estimated result as the initial value for re-entry prediction, even though residuals for each pass could be small. Almost all of the earlier studies assume the ballistic coefficient is constant, which might not be a reasonable choice for objects with eccentric orbits. This paper presents a new method for re-entry prediction of objects with loweccentricity orbits. This new method is based on mean ballistic coefficients, which reflect the average effect of atmospheric drag during one orbital revolution. The method is developed in Section 2. Then in Section 3 the method is tested using the Iridium-52 satellite which uses sparse observations as the data source, and ten other objects with low-eccentricity orbits which use TLEs as the data source. Evolution of drag characteristics and orbit propagation using the mean ballistic coefficient are also discussed in Section 3. Finally, some conclusions are drawn in Section 4.

Method
The mean ballistic coefficient estimation method and the reentry prediction method are laid out in this section. In this paper we use classical orbital elements σ (a, e, i, Ω, ω, M) to describe the orbital motion, where (a, e, i, Ω, ω, M) are semi-major axis, eccentricity, inclination, right ascension of the ascending node, argument of perigee, and mean anomaly, respectively.

State Estimation
State estimation can be carried out based on observations or TLEs. When multiple passes of observations are used, problems such as large residuals and non-convergence may arise for the re-entering objects with low-eccentricity orbits, especially when there are sparse observations. Therefore, we perform state estimation based on a single pass of observations using a numerical method. Most widely used methods for determining a precise orbit are suitable for this purpose. A weighted least-squares algorithm is used in this paper when determine the precise orbit based on observations. The estimation results are state vectors, (⃗ r,⃗ r), which can easily be transformed to osculating orbital elements. Ballistic coefficients are estimated in some cases, for example, when the observations are obtained in a altitude region below 500 km. We also use corresponding mean orbital elements. The mean orbits used in this paper are obtained by removing the short periodic variations from the osculating orbits (Liu 2000) as given by Equation (2): where σ is the mean orbital element, including secular terms and long periodic terms, σ is the osculating orbital element, and ∆σs are short periodic terms. An analytical scheme for calculating mean orbital elements was proposed by Kozai (Kozai 1959) and Brouwer (Brouwer 1959). More details for short periodic terms have been discussed in (Liu 2000;Kozai 1959;Brouwer 1959;Zhong and Gurfil 2013). When TLE is used, the state can be obtained by converting TLE to an osculating state at the desired epoch using SGP4 (Vallado et al. 2006), or determining the orbit based on pseudo-observations derived from TLE (Vallado et al. 2013;Levit and Marshall 2011).

Mean Ballistic Coeflcient Estimation
How to estimate the ballistic coefficient of the re-entry object has been the focus of many studies. In most of these studies the ballistic coefficient is considered as a constant. This might not be sufficient for objects with low-eccentricity orbits, because the drag coefficients change with the altitude. When the space object's attitude is stable, the change of the ballistic coefficient is mainly caused by the drag coefficient. The drag coefficient depends on many parameters, such as the object's shape, the surface material, and the composition and temperature of the atmosphere. The drag coefficient of a plate whose normal aligns to the flight direction can be calculated from Equation (3) (Pilinski et al. 2011): where s is the speed ratio, defined as the ratio of the relative velocity of the re-entry object with respect to the atmospheric molecule and the most probable molecular speed; T k,r is the temperature of the reflected molecule; T i is the kinetic temperature carried to the surface by the incoming molecule; and erf() denotes the error function defined as Equation (4): To deal with the drag characteristic we propose to carry out re-entry prediction by employing the mean ballistic coefficient defined in this paper. The mean ballistic coefficient B reflects the average effect of atmospheric drag in one orbital revolution, which is defined by Equation (5): where ρ is the atmospheric density. Atmospheric drag has a long-term effect on orbital variations. The change in the semi-major axis due to the drag effect during one revolution ∆a can be defined by Equation (6): where E is the eccentric anomaly. Obviously, the mean ballistic coefficient can be easily estimated using the seminumerical method, which integrates the mean orbital elements with a step size that is a multiple of the orbital period. The following perturbations are considered: the central body perturbation, J 2 and J 3 zonal harmonics of the Earth potential, and perturbation due to atmospheric drag. Modeling higher order perturbations is much more difficult and complex for a semi-numerical method (Wu 2011). In addition the long-term impact of other perturbations is small compared with atmospheric drag (Picone et al. 2002). Thus, more accurate perturbations are deemed to be unnecessary. For convenience the mean ballistic coefficient is written as B in the following sections. Suppose that N sets of mean orbital elements (σ 1 , σ 2 , · · · , σ N ) have been obtained in the state estimation stage, and the corresponding epochs are (t 1 , t 2 , · · · , t N ). Then the procedure to estimate the mean ballistic coefficient is as follows: 1. Set the last orbital element σ N and the initial ballistic coefficient B 0 as the initial values of the integration; B 0 could be one of the estimated results in observation case, or set as 12.741621B * in TLE case, where B * is the drag item in TLE. 2. Propagate backward by integrating the following equation, Equation (7), with a step size of one period (Liu and Alford 1980). By propagating backward one prevents re-entry occurring during propagation.
where ρ is the atmospheric density, computed by the MSIS-90 atmospheric model, v is the relative velocity of the re-entry object with respect to the atmosphere, f is the true anomaly, ω E is the Earth rotation rate, µ is the Earth gravitational constant, r E is the Earth equatorial radius, n is the mean motion, p = a(1 − e 2 ), and J 2 and J 3 are the second and third gravity harmonics of the Earth's field. 3. When one integration step is finished, compare the updated time with the epoch of the mean orbital elements. If the difference is within half a step size, compute the difference between the computed semimajor axis and the "real" semi-major axis, ∆a k , and then take the partial derivative about the ballistic coefficient, ∂a k ∂B . Some researches estimate the ballistic coefficient by employing the semi-major axis and eccentricity simultaneously. But the eccentricity is generally less accurate than the semi-major axis (Hejduk et al. 2013). Therefore, we consider only the change in semi-major axis. 4. After the first orbital element σ 1 is compared, compute the correction value for B using Equation (8): 5. If ∆B is less than the convergence threshold, the estimation process is finished. If not, correct the coefficient B as follows: B ′ = B + ∆B, repeat step 1 and begin a new iteration.

Re-entry time prediction
To predict the re-entry time we set the last osculating orbital element σ N and the estimated mean ballistic coefficient B as the initial values, and then propagate the orbit until reentry. We use a 10th-order fixed step size KSG integrator for numerical integration. Dynamical and geometrical models employed in the re-entry prediction are shown in Table 1. Because atmospheric densities from different atmospheric models may be quite different from each other (Bruinsma et al. 2004), the atmospheric models used to estimate the ballistic coefficient and predict the re-entry time should be consistent. Here we use the MSIS-90 atmospheric model.

Case studies
To measure the prediction accuracy we use the error between the predicted and real re-entry time, which is calculated using Equation (9): where δ is the relative error, t real is the real re-entry time, t pred is the predicted re-entry time, and t 0 is the epoch of the initial state. For the objects in nearly circular orbits, the relative error generally is 10% to 20% and sometimes larger than 30%. The error mainly depends on the attitude of the re-entering objects and the space environment during reentry (Choi et al. 2017;Pardini and Anselmo 2008;Geul et al. 2018). Iridium-52 and ten other objects with low-eccentricity orbits which re-entered in 2017 are used as test objects.
Re-entry predictions of Iridium-52 are carried out based on sparse observations with time intervals larger than ten hours; and re-entry predictions of other objects are carried out based on TLEs accessed from the space-track website.

Iridium-52
Iridium-52 (NORAD No. 25169) was deployed into a 760 km circular orbit on 18 February 1998. According to the TLEs, Iridium-52 was deorbit maneuvered on 16 October 2018, and when the deorbit maneuver finished on 24 October 2018, the perigee altitude was 162 km and the apogee altitude was 668 km. Most of the Iridium satellites were deorbited by lowering their perigee altitudes. That is why there are a large number of objects with low-eccentricity orbits which have re-entered since 2018. The initial perigee altitude of Iridium-52 was about 153 km and initial apogee altitude was about 612 km. According to the message from the spacetrack website, the re-entry time of Iridium-52 was 2018-11-05 21:20(UTC), which is consistent with our last prediction result. Hence, we adopt 2018-11-05 21:20(UTC) as the real re-entry time. Re-entry predictions of Iridium-52 are carried out based on three successive passes of ground-based radar observations, with time intervals larger than ten hours since 2018-10-27. Most passes last 5-7 minutes with measurement frequency at 1 Hz. Take the computation on 2018-10-30 as example, the observations information is shown in Table  2. State vectors and ballistic coefficients are estimated in orbit determination computation, the ballistic coefficients relative to the observations are also shown in Table 2. We   can see the orbit determination problem stated previously, that we estimate multiple ballistic coefficients relative to observations to obtain small residuals, it is difficult to select a suitable estimated result as the initial value for re-entry prediction.
The mean ballistic coefficient and the re-entry time results are shown in Table 3 and Figure 2. To understand clearly the performance of the mean ballistic coefficient, the real environmental indices are used in the atmospheric density calculation, thus reducing the errors caused by inaccuracies in the atmosphere. We can see that the mean ballistic coefficient of Iridium-52 varies between 0.0251 and 0.0256, and the mean value is 0.0255. The variability rate of the mean ballistic coefficient is less than 1%, which is much smaller than the variations of the drag coefficient in one revolution. Eight of the ten prediction results have relative errors less than 1%, and the maximum relative error is 3.25%. This shows the efficiency and accuracy of the mean ballistic coefficient used in re-entry prediction.

The ten test objects
We also use ten other objects with low-eccentricity orbits which re-entered in 2017 as test objects. The difference is that TLEs are used as the data source for these ten objects. Five TLEs over five days are used to estimate the mean ballistic coefficients for each object, and the TLEs are not preprocessed. The states are obtained by converting TLE to an osculating state at the TLE epoch using SGP4, and then transformed to osculating orbital elements and mean orbital elements. The information and real re-entry time of these ten objects are shown in Table 4. The real re-entry times are obtained from the TIP (Tracking and Impact Prediction) messages from the space-track website. If no TIP messages are released, the last predicted result is adopted as the real re-entry time. Re-entry predictions are carried out once a day over ten days before re-entry for each object, and the maximum relative errors are shown in Figure 3. All the re-entry predictions are carried out with real data, which means that the predicted environmental indices are used and the inaccuracy in the atmosphere is involved. Clearly, the maximum relative errors for eight objects are less than 15%, and the other two objects have maximum relative errors of 35.6% (NORAD No. 41726) and 42.5% (NORAD  No. 09980), respectively. It is also seen that the maximum relative error has no apparent correlation with initial eccentricity (or initial apogee altitude). Poor accuracy of TLEs might be the primary reason for the large relative error for these two objects.

Mean Ballistic Coeflcient Used in the Evolution of Drag Characteristics
Even within the same revolution the ballistic coefficient of the object with a low-eccentricity orbit is not constant, and that makes it difficult to select initial values for re-entry prediction. The ballistic coefficient derived from observations tracked in a certain altitude might be difficult to propagate to the other altitudes. Therefore, the traditional method for re-entry prediction, which estimates the ballistic coefficients from orbit determination computations, is limited to objects with low-eccentricity orbits. Suppose that the attitude of the object is stable, it is reasonable to agree that the change of the ballistic coefficient is caused by the drag coefficient. Use CZ-3B 3rd stage (NORAD No. 38253), which is one of the ten re-entered objects in 2017, as the test case. The initial osculating orbital element of CZ-3B 3rd stage is shown in Table 5.
To simplify the modeling we suppose that the object in the orbit of CZ-3B 3rd stage is a flat plate, whose normal aligns to the flight direction. Calculate C D of the flat plate on August 8 for one revolution using Equation (3); the result is shown in Figure 4a. The mean drag coefficient C D is 2.167 as calculated by Equation (10). Simulate the drag characteristic (defined as the product of C D and atmospheric density) of the flat plate for one day, three days, and seven days later, respectively, as shown in Figure 4b, Figure 4c, and Figure 4d. The value of the polar axis is lg(C D × ρ).
(10) Figure 4a shows how the drag coefficient changes with the altitude in one revolution. The drag coefficient is about 2.15 at perigee, and about 3.24 at apogee. Use the mean drag coefficient (2.167) to calculate the subsequent drag characteristics. The results agree with the results calculated using the real coefficient. The obvious difference occurs at the apogee. The real drag coefficient at the apogee is obviously different from the mean drag coefficient; however, because the density at apogee is three orders smaller than the density at perigee, the difference can be neglected. The error due to the mean ballistic coefficient is smaller than 0.2% after three days and smaller than 0.3% after seven days. But the error due to the mean ballistic coefficient of the flat plate in the nearly circular orbit is greater by about 3% after seven days, which is one order greater than the flat plate in a low-eccentricity orbit.   Therefore, we conclude that the mean ballistic coefficients are quite suitable for re-entry predictions of objects with low-eccentricity orbits. The change with altitude in the drag coefficient should be considered, if more accurate predictions are required for the objects with nearly circular orbits.

Mean Ballistic Coeflcient Used in Orbit Propagation
The position error in the orbit propagation for over 24 hours might be larger than several kilometers, even tens of kilometers for re-entering objects in elliptical orbits. Large position errors cause tracking problems for observer stations. Therefore, we are also interested in whether the mean ballistic coefficient could improve the accuracy of orbit propagation.
We use the observations of Iridium-52 to test this problem. To show the effect of the ballistic coefficients obtained from different combinations of observations, we select observations tracked from a low altitude region and a high altitude region. The observations are shown in Table 6. Two orbits are derived from these five passes of observations with the combination of pass 1, 2 and the combination of pass 3, 4, and 5 respectively. The ballistic coefficients are estimated in the computation of the determination of the orbit, and then the orbits are propagated using the results of the orbit determination, the estimated ballistic coefficient, and the mean ballistic coefficient, marked as "EBC-L", "MBC-L," "EBC-H," and "MBC-H," respectively. "EBC-L" means that the orbit is propagated based on the state result and the ballistic coefficient result derived from observations obtained from the low altitude region, "MBC-L" means that the orbit is propagated based on the state result derived from observations obtained from the low altitude region and the mean ballistic coefficient result. It is similar to "EBC-H" and "MBC-H". Finally, the propagated orbits are compared with the observations obtained over the next three days, then we get the position errors. The position errors are shown in Figure 5.
The ballistic coefficient derived from observations obtained from the high altitude region, namely passes 1 and 2, is 0.0281. The mean ballistic coefficient is 0.0255, obtained from Table 2. The total position errors are shown in Figure  5 using red marks. We can see that the position error using estimated ballistic coefficient is 107 km one day later and 939 km three days later. It can be concluded that the large error in the ballistic coefficient causes the large position error. The ballistic coefficient derived from observations obtained from the low altitude region, namely passes 3, 4, and 5, is 0.0249, which is close to the mean ballistic coefficient. The position errors for this case are marked in blue in Figure 5. We can see that "EBC-L" is much more accurate than "EBC-H," but less accurate than "MBC-L." "MBC-H" is a little less accurate than "MBC-L." The difference might be caused by the difference in the initial state. It can be con-cluded that the mean ballistic coefficient is optimal for orbit propagation. The ballistic coefficient derived from observations obtained from the low altitude region is next, and the third is the ballistic coefficient derived from observations obtained from the high altitude region.

Conclusions
Because objects with low-eccentricity orbits traverse a large portion of the dense atmospheric region almost every revolution, the drag coefficients change significantly in one orbital revolution, thus making it more difficult to obtain predictions for re-entry. We propose a new method for reentry prediction based on the mean ballistic coefficient, which reflects the average effect of atmospheric drag during one revolution. Then the method is tested using Iridium-52 and ten other objects with low-eccentricity orbits. The variability rate of the mean ballistic coefficient of Iridium-52 is less than 1%, which is much smaller than the variation of the drag coefficient in one revolution. When the mean ballistic coefficients are used in re-entry predictions, the prediction results can be improved.
We also discuss the evolution of the drag characteristics and orbit propagation using the mean ballistic coefficient. We find that the mean ballistic coefficients are quite suitable for re-entry predictions and orbit propagation of those objects with low-eccentricity orbits. When the mean ballistic coefficient is unavailable, an alternate choice can be to estimate the ballistic coefficient from the observations tracked from the low altitude region.