Abstract
Spent rocket bodies in geostationary transfer orbit (GTO) pose impact risks to the Earth’s surface when they reenter the Earth’s atmosphere. To mitigate these risks, reentry prediction of GTO rocket bodies is required. In this paper, the reentry prediction of rocket bodies in eccentric orbits based on only Two-Line Element (TLE) data and using only ballistic coefficient (BC) estimation is assessed. The TLEs are preprocessed to filter out outliers and the BC is estimated using only semimajor axis data. The BC estimation and reentry prediction accuracy are analyzed by performing predictions for 101 rocket bodies initially in GTO and comparing with the actual reentry epoch at different times before reentry. Predictions using a single and multiple BC estimates and using state estimation by orbit determination are quantitatively compared with each other for the 101 upper stages.
1. Introduction
Rocket bodies in geostationary transfer orbits (GTOs) have their apogee near geosynchronous altitude and their perigee within the Earth’s atmosphere. The atmospheric drag reduces the orbital energy of the rocket bodies and lowers the orbit until reentry occurs. Lunisolar perturbations speed up or slow down this process by changing the eccentricity of the orbit and raising or lowering the perigee altitude, which in extreme cases results in direct reentry without drag-induced decay. The reentry of spent rocket bodies is desirable because the deorbiting of these uncontrolled bodies prevents collisions with functional spacecraft and potential generation of new space debris. However, the reentry poses a risk to the Earth’s population because rocket bodies consist of components likely to survive the reentry and impact the Earth’s surface (such as propellant tanks) [1]. Therefore, to be able to mitigate any risks due to deorbiting, the reentry of rocket bodies needs to be predicted.
The major source of error in orbit prediction is the computation of the atmospheric drag [2]. The perturbing acceleration due to drag, , depends on the spacecraft’s drag coefficient, , area-to-mass ratio, , velocity with respect to the atmosphere, , and the atmospheric density, :The drag coefficient and the effective area-to-mass ratio depend on the object’s attitude, which is generally uncertain. The local atmospheric density, on the other hand, depends on the solar and geomagnetic activity, for which future values are unknown [3, 4]. In addition, the drag calculation is subject to inaccuracies in the atmospheric density model and possible mismodeling of the drag coefficient [1]. Finally, the velocity with respect to the atmosphere is uncertain, because the local wind speed is unknown.
For state-of-the-art reentry prediction, the accuracy of atmospheric density calculations can be improved by calibrating the density models using near real-time satellite tracking data [5–7]. In addition, the effective area can be computed by performing six degrees-of-freedom (6DoF) propagation to calculate the attitude of the rocket body [8]. Moreover, using the attitude and a physical model of the rocket body the drag coefficient can be computed [8, 9]. Furthermore, a wind model can be used to compute the horizontal wind speeds in the atmosphere [10].
When density correction models and 6DoF propagation techniques are not available (e.g., because the object details are unknown or the measurements necessary for density corrections are unavailable), the drag coefficient and area-to-mass ratio can be combined into one parameter called the ballistic coefficient () that can be estimated from orbital data. Such an estimated BC depends on the actual and area-to-mass ratio but also soaks up atmospheric density model errors and possibly other errors, for example, orbital data inaccuracies. More accurate orbital data and dynamical models therefore result in estimated BCs that are closer to the true BC [6].
The application of highly accurate models and orbital data is required for accurately predicting the impact point of reentering objects. Sufficiently accurate orbital data is, however, often not available and Two-Line Element sets (TLEs) provided by the United States Strategic Command are the only available data to perform reentry prediction. The accuracy of TLE data is, however, limited due to the application of simplified perturbation models (SGP4 and SDP4) [11, 12], especially for objects in GTOs [13, 14] and in orbits with high energy dissipation rates [15].
In this paper, the reentry prediction of rocket bodies in eccentric orbits based on only TLE data is assessed. Because attitude and density correction data are not directly available from TLEs, the predictions are carried out using 3DOF propagation and a standard empirical atmospheric density model. Different methods have been developed in the past to improve TLE-based reentry prediction by preprocessing TLE data and by estimating the BC, solar radiation pressure coefficient (SRPC), object state vector, or a combination of these. In this paper, reentry predictions using only an estimate for the BC are investigated. This approach is straightforward and can be used to obtain a first-order guess of the reentry date several weeks or months before reentry when accurate prediction of the impact point is not feasible due to uncertainties in future space weather predictions. In addition, reentry predictions using only BC estimates can easily be automated to perform daily predictions for many objects. Within this assumption (only BC estimation), the goal of this paper is to provide guidelines on how to estimate the BC to obtain the most accurate reentry predictions.
Ballistic Coefficient Estimation. For the estimation of the BC based on TLEs, several methods have been developed [16–20]. Saunders et al. [17] and Sang et al. [18] estimate the BC by comparing the change in semimajor axis according to TLE data with the change in semimajor axis due to drag computed by propagation using an initial state from TLEs. This method is straightforward and uses semimajor axis data from TLEs which are generally accurate. The methods by Saunders and Sang are almost equivalent, the main difference is that Sang computes a single BC estimate directly, where Saunders finds improved estimates by iteration. Gupta and Anilkumar [20], on the other hand, estimate the BC by minimizing the difference between apogee and perigee altitudes according to TLEs and propagation. This method is said to perform well for reentry prediction during the last phase of orbital decay. It is, however, more complex and requires the use of the eccentricity from TLEs which is generally less accurate than semimajor axis data. A method for estimating both the BC and initial eccentricity was developed by Sharma et al. [16] to improve reentry prediction of upper stages in GTO [21–23]. Here the eccentricity and BC are estimated by fitting the apogee altitude according to propagation to TLE apogee data using the response surface methodology. Finally, Dolado-Perez et al. [19] developed a method for estimating the BC and SRPC simultaneously. This is carried out by comparing the rate of change of the semimajor axis and eccentricity according to TLE data and propagation. The method assumes that the change in semimajor axis is due to both drag and SRP, which should improve the BC estimate. However, again less accurate eccentricity data from TLEs are used for the estimation. In addition, because the eccentricity is strongly affected by lunisolar perturbations, the changes in eccentricity due to drag and SRP are hard to observe. Finally, the methods by Sharma et al. [16] and Gupta and Anilkumar [20] estimate a single BC that is used for the purpose of reentry prediction. Saunders, Sang, and Dolado-Perez, on the other hand, estimate multiple BCs and subsequently take a statistical measure of the set as final estimate.
It should be noted that all these methods estimate a single, and thus fixed, ballistic coefficient. In reality, the BC, however, varies over time due to, for example, rotation of the object or changes in due to altering atmospheric conditions. Efforts can be made to predict the future variation of the BC [24] or assume a relation between the drag coefficient and the orbital regime [25], but this is beyond the scope of this paper.
State Estimation. To obtain an accurate state of the object for reentry prediction, state estimation can be carried out by orbit determination using pseudo-observations derived from TLE data. This approach is widely used and is described by, for example, Levit and Marshall [26], Vallado et al. [14], and Dolado-Perez et al. [19]. In this paper, state estimation will only be utilized for comparison.
TLE Preprocessing. TLE data is used for estimating the BC and state of an object; however, the quality of TLEs associated with an object is not homogeneous: sometimes low quality or even wrong TLEs are distributed. For this reason, preprocessing of TLEs is needed to identify outliers and TLEs of poor quality [27].
TLE Based Reentry Prediction Approach. The goal of this paper is to obtain accurate reentry predictions of decaying GTO rocket bodies using only an estimate for the BC and irrespective of TLE quality and availability. This is achieved by TLE preprocessing (see Lidtke et al. [27]) and enhancing the BC estimation for the purpose of reentry prediction. The main contributions of this work are as follows.(i)The estimation of the BC is tailored for reentry predictions by comparing the decay of the mean semimajor axis according to TLE data and according to a high-fidelity propagator considering all perturbations.(ii)The impact of the initial state used for BC estimation on the reentry prediction is shown.(iii)The performance of the method is assessed and improved based on predicting the reentry dates of 101 upper stages in highly eccentric orbits (all initially in GTO) and the sources of inaccurate predictions are analyzed.(iv)The good performance of using a single BC estimate versus the use of a median BC estimate and versus BC and state estimation is shown.
Because the considered rocket bodies are in highly eccentric orbits, all relevant perturbations (geopotential, lunisolar, drag, and SRP) are always considered during orbit propagation.
The methods used in this approach are discussed in the following section. After that the BC estimation and reentry prediction results using a single and multiple BC estimates are discussed.
2. Methods
The orbital propagator and BC and state estimation and TLE preprocessing methods used for TLE-based reentry prediction are discussed in the following.
2.1. Propagation Method
The orbital propagator used in this study is the Accurate Integrator for Debris Analysis (AIDA), a high-precision numerical propagator tailored for the analysis of space debris dynamics using up-to-date perturbation models. AIDA includes the following force models [28]: geopotential acceleration computed using the EGM2008 model (10 × 10), atmospheric drag modeled using the NRLMSISE-00 air density model, solar radiation pressure with dual-cone shadow model, and third body perturbations from Sun and Moon.
NASA’s SPICE toolbox (https://naif.jpl.nasa.gov/naif/index.html) is used both for Moon and Sun ephemerides (DE405 kernels) and for reference frame and time transformations (ITRF93 and J2000 reference frames and leap-seconds kernel). Solar and geomagnetic activity data (F10.7 and Ap indexes) are obtained from CelesTrak (http://www.celestrak.com/SpaceData/sw19571001.txt) and Earth orientation parameters from IERS (ftp://ftp.iers.org/products/eop/rapid/standard/finals.data). A wind model is not used, because the effect of wind generally cancels out over one orbital revolution [29] and the impact of neglecting wind is small compared to the effect of inaccuracies in atmospheric density modeling.
2.2. Ballistic Coefficient Estimation Method
The approach used for the estimation of the BC is based on the method for deriving accurate satellite BCs from TLEs proposed by Saunders et al. [17]. Several modifications were made to improve the method for the reentry prediction purpose. The BC estimation algorithm uses the data of two TLEs. The BC is estimated by comparing the change in semimajor axis according to two TLEs to the change in semimajor axis due to drag computed by accurate orbit propagation using an initial state derived from the first TLE (if not stated otherwise, states are obtained from TLEs using SGP4 to convert the TLE to an osculating state at the desired epoch and subsequently converting the state from the TEME to J2000 reference frame). Since short-periodic changes are removed from TLE data, the change in semimajor axis according to TLEs can be assumed to be purely the secular change caused by atmospheric drag (long-periodic variation of semimajor axis due to gravitational terms and SRP may be included in TLE data but are generally small compared to changes due to drag [30]). Therefore, any difference between the change in semimajor axis according to TLE data and due to drag computed by orbit propagation can be assumed to be caused by a wrong guess for the BC. The BC that gives the correct change in semimajor axis is obtained as follows:(1)Compute the change in semimajor axis between the two TLEs, , using the “mean” mean motion, , available in a TLE:(2)Take guess for value of the BC.(3)Propagate the orbit with the full dynamical model between the two TLE epochs and simultaneously compute where is the semilatus rectum, the true anomaly, and and the acceleration due to drag in radial and transverse direction, respectively.(4)Integrate over time to obtain the change in semimajor axis due to drag only, ,(5)Update the BC estimate value using the Secant method: where is the th BC estimate and .(6)Repeat the procedure from step 3 until convergence is reached.
The first guess, , for this method is taken from of the first TLE. The parameter in TLEs is an SGP4 drag-like coefficient and a BC value can be recovered from it: [31]. The second guess, , needed for the Secant method is computed by performing one propagation using the first guess and assuming a linear relation between the BC and : The convergence criterion is met when is less than km.
Several changes were made to the original method by Saunders. First, during the BC estimation process, it may happen that the object unexpectedly reenters during propagation. Such a reentry is generally the result of a too-high estimate for the BC. Therefore, the propagation is then repeated assuming a smaller value for BC, namely, 90% of the initial value. This prevents failure of BC estimation due to reentry but may require several iterations to sufficiently reduce the BC value.
By default forward propagation is applied for BC estimation, that is, taking the state at the earliest TLE and propagating it until the epoch of the latest TLE. In addition, also backward propagation was implemented starting from the latest TLE and propagating backward until the prior one. By propagating backward one prevents reentry occurring during propagation. This is especially useful when estimating the BC close to reentry where an inaccurate BC guess can easily cause unexpected reentry.
Furthermore, the change in semimajor axis due to drag (see (3)) is computed considering all perturbations during propagation. This is important because the effect of coupling between different perturbations cannot be neglected.
Finally, the average semimajor axis is computed from osculating data from AIDA to compare the change in semimajor axis with TLE data. This improves the estimation because the osculating data includes short-periodic variations whereas the mean TLE data does not [30].
Besides estimating the BC also the SRPC can be estimated. Dolado-Perez et al. [19] developed a method where the BC and SRPC are estimated simultaneously by comparing semimajor axis and eccentricity data from TLEs with the changes in semimajor axis and eccentricity due to drag, SRP, and conservative forces. This method was implemented and tested but was found to give aberrant results because in all test cases the effect of SRP was at least an order of magnitude smaller than the effect of drag. This resulted in an ill-conditioned system of equations and consequently aberrant SRPC estimates. Therefore, SRPC estimation was omitted and known area-to-mass ratio data was used to compute the SRPC for SRP perturbation computation, assuming the typical reflectivity coefficient value of .
2.3. State Estimation
The state estimation performed in this work is carried out by fitting accurate orbit propagation states to pseudo-observations derived from TLEs using nonlinear least-squares. This is a consolidated method widely used for offline (ground-based) orbit determination (OD) [32]. A five-day observation window with 21 pseudo-observations is used to estimate the state together with the BC. The initial state is located at the end of the observation period and is expressed in modified equinoctial elements [33]. The residuals minimized during least-squares optimization are expressed in Cartesian coordinates aligned with satellite coordinate system in radial, transverse, normal directions. More details on the algorithm and settings can be found in Gondelach et al. [34].
2.4. TLE Preprocessing
The TLEs have to be filtered because incorrect, outlying TLEs and entire sequences thereof could be present in the data from Space-Track, and using such aberrant TLEs in subsequent analyses would deteriorate the accuracy of the results. Filtering out aberrant, or incorrect, TLEs consists of a number of stages [27]; namely,(1)filter out TLEs that were published but subsequently corrected;(2)find large time gaps between TLEs because they hinder proper checking of TLE consistency;(3)identify single TLEs with inconsistent mean motion, as well as entire sequences thereof, using a sliding window approach;(4)filter out TLEs outlying in perigee radius;(5)filter out TLEs outlying in inclination;(6)filter out TLEs with negative as they cause incorrect SGP4 propagation.
TLEs with negative are filtered out, because they produce SGP4 propagations where the semimajor axis increases, which is not realistic for decaying orbits. More details on the applied filtering methods and results are discussed by Lidtke et al. [27].
3. Test Cases
To determine the quality of the BC estimates, the estimates were compared with BC values derived from in TLEs and with real object data. In addition, to measure accuracy of the reentry predictions, the error between the predicted and actual reentry date is computed. This error with respect to the time to reentry is calculated as follows: where is the predicted reentry date, the actual reentry date, and the epoch of the last TLE used for the prediction.
To test the reentry prediction performance, a set of 101 rocket bodies that reentered in the past 50 years was selected. This makes it possible to compare the predicted reentry date with the real one. The reentry dates were taken from satellite decay messages from the Space-Track.org website (https://www.space-track.org) that provides the decay date of space objects. It is worth mentioning that the exact reentry time is not known, because all decay times are at midnight (this can produce a bias in the calculated reentry prediction error when predictions are made close to the actual reentry). All upper stages were initially in GTOs, but their reentry dates, lifetimes, inclinations, and area-to-mass ratios differ significantly. To give an indication, the perigee altitude 180 days before reentry lies between 131 and 259 km and the eccentricity between 0.1 and 0.73. The number of TLEs available in the last 180 days before reentry varies from 45 to 543 and the area-to-mass ratio according to object data lies between 0.002 and 0.03 m2/kg.
In addition, all objects have been used to predict the reentry 10, 20, 30, 60, 90, and 180 days before the actual reentry date. Some of the 101 objects were not suitable for several reentry prediction tests, because they had no TLEs within a specific number of days before the reentry (e.g., last TLE is 90 days before reentry).
In real reentry prediction cases, the actual reentry date of the object is, of course, not known. Analyzing the results has therefore not only the goal to examine the quality of the reentry predictions but also the goal to define guidelines for real reentry prediction scenarios.
4. Results
4.1. Ballistic Coefficient Estimation
Figure 1 shows BC estimates and BCs from for object 28452 together with the perigee radius according to TLE data in the 180 days before reentry. For the left plots TLEs filtered on mean motion were used, whereas for the right plots the TLEs were filtered on mean motion and perigee radius. First of all, the trend of the BC estimates is similar to the trend of the BC from , but with an offset (note that in general it is however not true that BC estimates and BC from follow the same trend). This proves that a BC estimate is required to perform reentry prediction with a dynamical model different from SGP4/SDP4.
(a)
(b)
(c)
(d)
Besides, there is a clear relation between outliers in TLE perigee radius and estimated BC; an outlier in perigee radius results in an outlier in the BC estimates. More precisely, of the two TLEs that are used for BC estimation, the outlying TLE that is used to obtain the initial state for propagation results in an outlier in BC estimate. The other TLE is only used to compute the change in semimajor axis according to the TLEs and does not have such a strong effect. Therefore, it can be concluded that the BC estimate strongly depends on the initial state used in the estimation. Because the atmospheric drag depends largely on altitude, an incorrect value of the initial state that translates in an aberrant perigee height results in a poor BC estimate. The BC estimate compensates for the incorrect initial state such that the state and BC together give the correct decay in the estimation period. is strongly correlated to the perigee height and thus both BC estimate and depend on the initial state. This may explain why the BC estimate and in Figure 1 follow the same trend.
Figures 1(b) and 1(d) show the BC estimates and perigee radius after filtering the TLEs on outliers in perigee radius. The BC estimates improve, because outliers in BC estimate disappear when TLE outliers in perigee radius are removed. Nevertheless, there are still outliers in the BC estimates, which may be removed when also smaller outliers in perigee radius are filtered out.
To have a closer look at the dependency of the BC estimate on the perigee radius, the BC estimates are plotted against perigee radius according to TLE data for object 27808 in Figure 2, where the color indicates the epoch of the BC estimate. In Figure 2(c) one can observe a correlation between the BC estimates and perigee radii for estimates at similar epochs. For a set of BC estimates with similar epochs, the BC varies almost linearly with changing perigee radius. Figures 2(a) and 2(b) show that this relation is mainly due to noise in the perigee radius that is compensated by the BC estimates. If the TLE data were more accurate then the BC estimates would not vary as much and would be closer to the real BC.
(a)
(b)
(c)
This proves that to obtain a good single BC estimate the TLEs should be filtered on perigee radius, or on both semimajor axis and eccentricity. Another option to reduce the impact of outliers on the estimate is to compute multiple BC estimates and take the median of the estimates as the final BC estimate. The reentry prediction results using a single and a median BC estimate are discussed in the next two sections.
Besides, different epoch separations between the two TLEs used for BC estimation have been tested, namely, 2, 5, 10, and 20 days. A TLE separation of 10 days was found to be least sensitive to outliers and short-period effects, because the difference between mean and median of the estimates was the smallest and the dispersion in terms of standard deviation and median absolute deviation was small as well. Therefore, 10-day separation is used for BC estimation, which is in agreement with Saunders et al. [17].
Finally, BCs were estimated for the 101 test objects in the 180 days before reentry. It was found that 80% of the medians of the BC estimates were within the range of possible area-to-mass ratio (assuming ) according to physical object data taken from European Space Agency’s DISCOS database (https://discosweb.esoc.esa.int); see Figure 3. This gives confidence that the estimation method provides good results.
4.1.1. Reentry Prediction Using Single BC Estimate
The objective of this section is to show that, for reentry prediction using only a BC estimate, it is of fundamental importance to run the reentry predictions using the same state that is used for BC estimation.
As described in Section 2.2, two TLEs are needed for estimating the BC; thus to run the subsequent reentry prediction one can use the state of either one of the two TLEs. Now, consider the test case of predicting the reentry for 91 rocket bodies 30 days before reentry; that is, all reentry predictions start from the state of the TLE at 30 days (). In one case, and an older TLE () are used for BC estimation. BC is estimated by propagating from the state of backward to and the state of is also used for the reentry prediction. This case is labeled “older TLE, same state.” In the second case, the BC estimation is performed using and a newer TLE () by propagating backward from to . Here the state (of ) that is used for BC estimation is not equal to the state (of ) that is used for the reentry prediction. This case is called “newer TLE, different state.” Figure 4 shows the cumulative distributions of the reentry prediction errors and their 90%-confidence regions (the 90%-confidence region is the interval where the true cumulative distribution is located with 90% probability. The width of the interval depends on the number of samples and is computed using the Dvoretzky-Kiefer-Wolfowitz inequality [35]) for both cases. One can see that although newer information is used in the second case, the first case, which uses an older TLE but the same state, results in more accurate reentry predictions. The difference between the prediction results of the two cases is significant, because the corresponding 90%-confidence intervals only overlap for small prediction errors. The use of the newer TLE only gives more accurate reentry predictions if the same state is used for BC estimation and reentry prediction; see case “newer TLE, same state” in Figure 4. For completeness, Figure 4 also shows the case “older TLE, different state” that results in less accurate predictions compared to using the “same state.”
(a) Cumulative distributions and 90%-confidence regions of reentry prediction errors using only an estimate for BC for 91 objects 30 days before reentry
(b) Schematic diagram of BC estimation
Using the same state for BC estimation and reentry prediction gives better results, because the BC estimate is computed such that together with the state it gives the correct decay rate of the semimajor axis in the estimation period. Using that BC estimate with another state will generally not result in the correct decay rate and the reentry prediction is thus more likely to be less accurate. Therefore, the same initial state for BC estimation and reentry prediction should be applied.
The reentry predictions using a single BC estimate that are presented in the following sections are computed using the “older TLE, same state” approach such that the latest available TLE is used for the initial state.
4.1.2. Reentry Prediction Using Multiple BC Estimates
Instead of using a single estimate, one can compute multiple estimates and take the mean or median of the set that may better represent the average BC behavior. This approach was tested by estimating the BC for every TLE between 90 and 30 days and from 180 to 60 before reentry and use the median of the estimates for reentry prediction at 30 and 60 days before reentry, respectively. The prediction errors are shown in Figure 5. Compared with the predictions based on a single BC the results are significantly worse; the majority of the median-BC samples is outside the 90%-confidence interval of the single-BC error distribution. On average, the reentry predictions are 8% and 6% less accurate at 30 and 60 days before reentry, respectively.
(a) 30 days before reentry; median taken from BC estimates between 90 and 30 days before reentry
(b) 60 days before reentry; median taken from BC estimates between 120 and 60 days before reentry
It was found that especially for orbits with a high eccentricity and low inclination the predictions with median BC are less accurate. Figure 6 shows the prediction error against eccentricity with different markers for different inclinations at 60 days before reentry (similar results were found for 30 days). The results with median BC show a correlation between increasing eccentricity and increasing error, whereas with a single BC estimate this correlation is less strong. In addition, the majority of the inaccurate predictions with median BC at lower eccentricity corresponds to low inclination orbits ( deg). A possible cause for this is the TLE accuracy, because the accuracy of TLEs for objects in HEO, GTO, and orbits with low inclination is less than for other objects [36]. This is also shown in Figure 7 that shows the dispersion of the mean perigee data (the median absolute deviation of detrended perigee data (the mean perigee radius data was detrended by subtracting the moving median from the data; see Lidtke et al. [27])) against eccentricity. The dispersion of the perigee data, that is, the noise, increases with increasing eccentricity. A single BC estimate can compensate for such inaccuracies by soaking up the error. However, when using a median BC the individual TLE errors are averaged out and not compensated for, except for possible biases.
These results suggest that estimation of the perigee altitude or eccentricity is required in order to improve the perigee data and thus the BC estimation and reentry prediction. Indeed, Sharma et al. [16] developed a method for estimating both the BC and eccentricity with good reentry prediction results for upper stages in GTO.
4.1.3. Only BC versus Full State Estimation
The reentry predictions using only BC estimates are compared with those after full state estimation using OD. Figure 8(a) shows the reentry prediction results for 30 days before reentry after only BC estimation (orange) and after full state estimation (blue). Surprisingly, the results obtained after OD are not better than the predictions using only an estimate for the BC. The BC-only predictions are on average 0.6% better; however, this difference is not significant for the number of samples (notice that the cumulative distributions are well within each others 90%-confidence intervals). This outcome is opposite to what one would expect, because a state estimated using OD is supposed to be a better starting point for accurate orbit propagation than a state taken directly from TLE data using SGP4. To check if state estimation improves reentry predictions at all, a test was performed where after the state estimation the BC is reestimated using the new state estimate. The results are shown in Figure 8(b) and they are on average 0.4% better than using only an estimate for the BC; however, again this difference is not significant for the number of samples used. This indicates that state estimation has less impact on the reentry prediction accuracy than BC estimation.
(a) Prediction errors using only BC estimate and after OD
(b) Prediction errors using only BC estimate and after OD with subsequent BC reestimation
To assess whether an accurate state and BC estimate result in an accurate reentry prediction, the six objects with the lowest position residuals after state and BC estimation using OD at 30 days before reentry were analyzed. Table 1 shows their mean position residuals and reentry prediction errors before OD (i.e., only BC estimation) and after OD. The residuals after OD are all two orders of magnitude smaller than before OD. The state estimation thus improved the accuracy of the orbit in the 5-day observation period significantly with respect to only estimating the BC. However, just half of the corresponding reentry predictions improved and the highest prediction error is still 16.6%. This shows that a state and BC that give an accurate orbit in the past do not necessarily give an accurate reentry prediction.
This outcome may be the consequence of taking a fixed BC for prediction. Figures 1 and 2 show that the BC changes over time (possibly due to object attitude variation, changing drag coefficient [25] and atmospheric modeling errors [6]). These variations in the BC are not accounted for during reentry prediction and, therefore, even if the initial state is very accurate, the prediction may not be accurate.
4.1.4. 10 to 180 Days before Reentry
Finally, the reentry prediction results for 10, 20, 30, 60, 90, and 180 days before reentry using single BC estimates are shown in Figure 9 together with the cumulative distribution and 90%-confidence interval of all predictions. The predictions at 60 days before reentry are on average most accurate. The predictions at 10 and 20 days before reentry, on the other hand, are significantly less accurate than the overall result. It should, however, be noticed here that the given reentry epochs are only accurate within one day (as they are given at midnight) which can result in a 10% reentry prediction error 10 days before reentry even if the prediction is perfect. The fact that the short-term predictions are less accurate is possibly due to the fast-changing dynamics close to reentry. The local atmosphere changes largely and the BC can vary quickly at lower altitudes; see, for example, Figure 1. Assuming a constant value for the BC may therefore not be a good approximation and accurate computation of the atmospheric drag becomes difficult.
(a) All prediction errors and at 10, 20, and 30 days before reentry
(b) All prediction errors and at 60, 90, and 180 days before reentry
Overall, with 90% confidence, 62 to 72% of the predictions is within 10% error and 85 to 95% within 20% error. Using a single BC estimate one can thus obtain a first-order estimate of the reentry date irrespective of TLE quality and availability. More sophisticated methods, such as 6DoF propagation and density corrections, should subsequently be applied to accurately estimate the impact point of the reentering object.
5. Conclusion
The estimation of the BC is tailored for reentry predictions by comparing the decay of the mean semimajor axis according to TLE data with the decay of the average semimajor axis due to drag according to a high-fidelity propagator considering all perturbations. The BC estimation results show that the estimated BC depends strongly on the initial state because TLE outliers and noise in the perigee radius result in outliers and noise in BC estimates. Therefore, filtering TLEs on eccentricity or perigee radius is important. Because of the dependency on the initial state, it is important to use the same initial state for BC estimation and reentry prediction as inaccuracy in the state is absorbed by a single BC estimate such that they provide the correct decay of the semimajor axis. Taking the median of multiple BC estimates for predicting the reentry does not give good results, because the median BC is not related to the initial state. The accuracy of reentry predictions after state and BC estimation using OD are not significantly different from using only a single BC estimate. Moreover an accurate initial state and BC do not necessarily give accurate reentry predictions. Overall, using a single BC estimate 62 to 72% of the reentry predictions is within 10% error (with 90% confidence). These conclusions are based on reentry predictions using TLE data and are thus subject to their accuracy and availability that vary largely for different objects.
Besides using more accurate orbital data, the fixed-BC approach can be improved by using more accurate atmospheric density models and by applying a wind model to increase the accuracy of density and velocity calculations during both BC estimation and reentry prediction. Furthermore, if the accuracy of the orbital data is very low, estimation of the eccentricity or perigee radius could improve the predictions as they strongly affect the BC estimate and reentry prediction. However, if the drag coefficient or frontal area of the object changes over time, then the achievable accuracy using a fixed BC is limited. Knowledge of the object’s attitude and 6DoF propagation or a forecasting model for the BC could significantly reduce the reentry prediction error.
Appendix
Test Objects
Rocket bodies with the following NORAD catalog numbers were used for reentry prediction:
625, 2609, 7252, 7794, 8479, 9017, 9787, 9859, 10983, 11072, 11718, 11719, 12562, 12810, 13025, 13087, 13098, 13136, 13294, 13447, 13599, 13684, 13940, 14130, 14168, 14287, 14332, 14369, 14423, 14787, 14989, 15157, 15165, 15679, 16600, 18352, 18923, 19218, 19332, 19877, 20042, 20123, 20254, 20778, 20920, 21057, 21141, 21654, 21766, 21895, 21990, 22118, 22254, 22906, 22928, 22932, 22997, 23315, 23416, 23572, 23797, 23916, 24314, 24666, 24770, 24799, 24847, 25051, 25129, 25154, 25240, 25313, 25372, 25496, 25776, 26560, 26576, 26579, 26641, 27514, 27719, 27808, 28185, 28239, 28253, 28418, 28452, 28623, 28703, 29497, 32764, 36829, 37211, 37239, 37257, 37482, 37764, 37805, 37949, 39499, 40142
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was partly carried out within the European Space Agency project ITT AO/1-8155/15/D/SR, titled “Technology for Improving Re-Entry Predictions of European Upper Stages through Dedicated Observations.” The authors acknowledge Dr. Hugh G. Lewis of the University of Southampton (UoS), Dr. Camilla Colombo of Politecnico di Milano, and Dr. Tim Flohrer and Quirin Funke of the European Space Agency for their valuable contributions. In addition, the use of the IRIDIS High Performance Computing Facility and associated support services at UoS in the completion of this work are acknowledged. David J. Gondelach was funded by an EPSRC Doctoral Training Grant awarded by the Faculty of Engineering and the Environment of UoS. Aleksander A. Lidtke would like to acknowledge the funding he received from the Ministry of Education, Culture, Sports, Science and Technology of Japan. Roberto Armellin acknowledges the support received by the Marie Skłodowska-Curie Grant 627111 (HOPT, Merging Lie perturbation theory and Taylor Differential algebra to address space debris challenges).