Mathematical Problems in Engineering

Volume 2017 (2017), Article ID 7309637, 13 pages

https://doi.org/10.1155/2017/7309637

## Ballistic Coefficient Estimation for Reentry Prediction of Rocket Bodies in Eccentric Orbits Based on TLE Data

^{1}Astronautics Research Group, University of Southampton, Highfield Campus, Southampton SO17 1BJ, UK^{2}Surrey Space Centre, University of Surrey, Guildford GU2 7XH, UK^{3}Department of Integrated System Engineering, Kyushu Institute of Technology, Kitakyushu, Japan

Correspondence should be addressed to David J. Gondelach

Received 30 June 2017; Accepted 14 November 2017; Published 10 December 2017

Academic Editor: Alessandro Gasparetto

Copyright © 2017 David J. Gondelach et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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 m^{2}/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.