# Effective Computational Approach for Prediction and Estimation of Space Object Breakup Dispersion during Uncontrolled Reentry

Authors: Deok-Jin Lee; Eun Jung Choi; Sungki Cho; Jung-Hyun Jo; Tae Soo No
Journal: International Journal of Aerospace Engineering
 (2018)

---

## Abstract

This paper provides an effective approach for the prediction and estimation of space debris due to a vehicle breakup during uncontrolled reentry. For an advanced analysis of the time evolution of space debris dispersion, new efficient computational approaches are proposed. A time evolution of the dispersion of space pieces from a breakup event to the ground impact time is represented in terms of covariance ellipsoids, and in this paper, two covariance propagation methods are introduced. First, a derivative-free statistical linear regression method using the unscented transformation is utilized for performing a covariance propagation. Second, a novel Gaussian moment-matching method is proposed to compute the estimation of the covariance of a debris dispersion by using a Gauss-Hermite cubature-based numerical integration approach. Compared to a linearized covariance propagation method such as the Lyapunov covariance equation, the newly proposed Gauss-Hermite cubature-based covariance computation approach could provide high flexibilities in terms of effectively representing an initial debris dispersion and also precisely computing the time evolution of the covariance matrices by utilizing a larger set of sigma points representing debris components. In addition, we also carry out a parametric study in order to analyze the effects on the accuracy of the covariance propagation due to modeling uncertainties. The effectiveness of the newly proposed statistical linear regression method and the Gauss-Hermite computational approach is demonstrated by carrying out various simulations.


---

## Body

## 1. Introduction

In the past 50 years, over 16,000 metric tons of man-made space objects and unexpected asteroids have entered the Earth’s atmosphere. Although most of them burn while entering the atmosphere, some of the debris have survived to impact the Earth’s surface and do expose a risk to people and property [1, 2]. To prevent these casualties, researchers have given considerable attention to the research area of reentry estimation of space objects. Usually, the starting altitude of a reentry object will be at around 120 km and then, breakup event happens at 78 km. However, the problem of estimation of space debris dispersion during the atmospheric reentry is still quite difficult due to unknown effects and uncertainties. Especially, the effects of atmospheric uncertainties are very difficult to determine since the atmospheric density and drag coefficient undergoes large fluctuations while decaying at that altitude [3]. In addition, the characteristics of uncontrolled space objects including the ballistic coefficients are largely insufficient to model the reentry trajectory. This is because the exact physical parameters of the objects like mass, area of cross section, shape, and dimensions are not available. Thus, the prediction of a reentry trajectory is generally estimated based on the statistical and stochastic method [4]. As shown in Figure 1, an estimated position and velocity vector of a reentering object is propagated by using a numerical integration model up to the breakup point at an altitude of 78 km. After the breakup occurred, the debris dispersion could be modeled by using various stochastic estimation approaches [5]. Several approaches for the estimation of the reentry time and impact locations have been proposed over the last decades [4–8]. For instance, Tardy and Kluever have estimated the states of orbital debris by using an optimal estimation method, and then Monte-Carlo simulation was carried out to predict the impact location with a final covariance matrix [5]. Reyhanoglu and Alvarado have proposed a Lyapunov equation-based covariance propagation method to predict the time evolution of debris trajectories [7]. In that method, a concept of positional probability ellipsoids is employed for the visualization of the time evolution of the debris dispersion. However, the Lyapunov equation-based covariance propagation method is based on the linearization of the translational equations of motion of the reentering space debris and the equations of the atmospheric reentry is subject to unknown uncertainties; thus, a high degree of neglected and truncated nonlinearities could lead to a degraded estimation of the debris dispersion [8, 9].Figure 1
Methodology for designing a reentry trajectory of an uncontrolled space object.In order to compensate for the drawbacks of the Lyapunov-based estimation of the debris dispersion, alternative solutions could be employed by increasing the order of the Taylor-series expansion of the nonlinear system or by using an advanced numerical technique [9, 10]. These efficient alternatives include the unscented transformation [11], the divided-difference numerical integration [12–14], and the cubature quadrature integration [15]. As discussed in [11], the unscented transformation is able to capture the higher-order moments caused by the nonlinear transform better than the Taylor-series-based approximations used in the Lyapunov-based covariance propagation in [5]. The unscented transformation-based estimation technique is called the sigma point estimator [16, 17] in the sense that the estimator works on the principle that a set of deterministically sampled sigma points is used to parameterize the mean and covariance of the Gaussian random variables without the linearization step. The classification of the sigma point method is also interpreted as a special case of a weighted statistical linear regression (SLR) [18]. Even though the statistical linear regression approach and the unscented transformation approach could represent the estimation of the debris reentry trajectories, the estimation accuracy of the debris dispersion could be limited, because they capture the debris dispersion with the second-moment covariance using a fixed small number of sigma points extracted from a predefined covariance ellipsoid. Therefore, it is necessary to design an advanced approach which could include more precise distribution information of the debris dispersion and take into account the nonlinearities due to uncertain nonlinear motion of the reentering debris.The main contributions of this paper are twofold. First, a precise and effective Gaussian cubature transformation-based [19, 20] estimation approach is proposed to compute the time evolution of the debris dispersion by using a moment-matching technique. The moment-matching formulation enables usage of many precise numerical integration methods such as Gauss-Hermite quadratures and cubature rules for multidimensional states [21]. In the Gauss-Hermite integration approach, a much larger set of sigma points and weights could be chosen compared to those of the Lyapunov covariance and the statistical regression approach such that with polynomial integrand the numerical approximation becomes exact and leads to more precise representation of the statistical moments (mean and covariance) of the distribution of a debris dispersion. Second, various parametric studies are carried out to analyze the effects of modeling uncertainties due to unpredictable atmospheric density and drag coefficients as well as wind effects near to the ground. For the analysis of the effects of the density model affecting the accuracy of the reentry prediction, four difference density models were used in the prediction of reentry trajectory. Specifically, NRLMSISE-00 [22], CIRA-72 [23], U.S. standard 1976, and exponential models [24] are compared with several different values of drag coefficient describing the unknown shape of breakup debris. During the atmospheric reentry, a time-varying empirical wind model is utilized to take into account a realistic environment, and it is based on HWM-07 (Horizon Wind Model 2007) [24, 25].The remainder of this paper is organized as follows: Section2 provides a description of the formulation of the reentry translational equations of motion. Section 3 provides two new approaches for the estimation and prediction of the debris dispersion due to the breakup during the reentry which is introduced by propagating the covariances of positional probabilistic ellipsoids. Section 4 shows performance comparisons of the proposed estimation approaches, that is, the statistical linear regression with the unscented transformation and the Gauss-Hermite cubature-based numerical approximation are demonstrated. Lastly, Section 5 provides a discussion of the simulation results.


## 2. Reentry Equations of Motion

### 2.1. Coordinate System

Letx=x1,x2,x3T denote the topocentric horizon coordinate system at the breakup instant. The origin of this frame is at θ0,ϕ0 on the Earth’s surface, where θ0,ϕ0 denotes the initial longitude and latitude. The x1x2 plane is the local horizon, which is the plane tangent to the sphere at the origin. x3 is normal to this plane directed towards the zenith (i.e., the coordinate x3 is the geometric altitude). The x1 axis is directed eastward and the x2 axis points north. The x1x2x3 frame is also referred to as the ENU (east-north-up) coordinate [24] as shown in Figure 2.Figure 2
ENU (east-north-up) coordinate system for reentry equations of motion.


### 2.2. Nonlinear Translational Equations of Motion

The governing equation of motion for an uncontrolled space object entering into the atmosphere perturbed by the atmospheric drag uncertainty but ignoring the lifting force can be described by the following equations with positionr and velocity v with their corresponding initial conditions rt0 and vt0 [5].
(1)ṙt=vt,v̇t=aDt+aGt−2ω×vt−ω×ω×rt+Ree^3+ξt,where r=x1,x2,x3T and v=v1,v2,v3T denote the position and velocity vectors, respectively, with respect to the ENU coordinate frame. Re=6378×103m is the Earth’s radius, e3=0,0,1T is the third standard unit vector in ℛ3 and ξ is the random acceleration vector due to modeling uncertainties and disturbances, and it is assumed that the noise vector is a zero-mean Gaussian process, EξtξtT=δt−sQt, where Qt is the positive definite process noise covariance matrix. ω=ω1,ω2,ω3T denotes the angular velocity vector of the ENU coordinates with respect to the ECEI (Earth-centered-Earth-fixed inertial) frame and each term of the angular velocity is computed by θ0,ϕ0 if the origin of the ENU topocentric horizon coordinate system is given ω1,ω2,ω3T=0,ωecosϕ0,ωesinϕ0T, where ωe=72.9217×10−6rad/s is the Earth’s rotation rate. aG denotes the gravitational acceleration, aG=ge^3, and based on the inverse square gravity model, the gravitational acceleration is calculated by
(2)g=geR2Re+x32,where ge denotes the gravitational acceleration at zero altitude. aD is the instantaneous acceleration due to the atmospheric density and it is assumed to be opposed to the direction of motion and proportional to the atmospheric density ρ and the velocity squared [24].
(3)aDt=−12βρvreltνrelt,β=msCdA,where β is the ballistic drag coefficient, Cd is the aerodynamic drag coefficient, A is the cross-sectional area of the satellite in a plane normal to the relative velocity vector, and ms is the mass of the reentering space object. vrel=v−w denotes the relative air velocity vector, where w is the wind velocity vector, and vrel=vrel=vrelTvrel is the magnitude of the relative air velocity vector. For the atmospheric density, an analytical exponential atmosphere model is utilized as
(4)ρx3=ρeexp−x3H,where ρe=1.752kg/m3 and H=6.7×103m.Now, after defining an augmented state vectorxt≡rTtvTtT, then the nonlinear differential equation could be written by
(5)ẋt=fxt+Gtηt,where
(6)fxt≡vtaDt+aG−2ω×vt−ω×ω×rt+Ree^3,and ηt≡03×1TξTtT is the augmented process noise vector, and Gt=03×3I3×3T.


### 2.3. State-Space Nonlinear Equations of Motion

It is assumed that the initial nominal state of a reentering space vehicle right before the breakup is given by propagating an orbital decay reentry trajectory until 78 km in ECEF (Earth-centered-Earth-fixed). After the coordinate transformation from the ECEF to ENU frame, the initial position of a breakup event is given by geometric altitudez0, longitude θ0, and latitude ϕ0. The initial position θ0,ϕ0 which locates the origin of the ENU frame on the Earth’s surface is given by
(7)x1t0,x2t0,x3t0=0,0,z0.Then, the nonlinear translational equations of motion in (1) can be rewritten in the form of a state-space representation [5] for numerical integration as follows:
(8)ẋ1t=v1t,ẋ2t=v2t,ẋ3t=v3t,v̇1t=−ρ2βvreltv1t−w1−2ω2v3t−ω3v2t+ωe2x1t+ξ1t,v̇2t=−ρ2βvreltv2t−w2−2ω3v1t−ω2ω3Re+x3t+ω32x2t+ξ2t,v̇3t=−ρ2βvreltv3t−w3−gx3+2ω2v1t−ω2ω3x2t+ω22x3t+Re+ζ3t,where w=w1,w2,w3T is the wind velocity vector.


## 2.1. Coordinate System

Letx=x1,x2,x3T denote the topocentric horizon coordinate system at the breakup instant. The origin of this frame is at θ0,ϕ0 on the Earth’s surface, where θ0,ϕ0 denotes the initial longitude and latitude. The x1x2 plane is the local horizon, which is the plane tangent to the sphere at the origin. x3 is normal to this plane directed towards the zenith (i.e., the coordinate x3 is the geometric altitude). The x1 axis is directed eastward and the x2 axis points north. The x1x2x3 frame is also referred to as the ENU (east-north-up) coordinate [24] as shown in Figure 2.Figure 2
ENU (east-north-up) coordinate system for reentry equations of motion.


## 2.2. Nonlinear Translational Equations of Motion

The governing equation of motion for an uncontrolled space object entering into the atmosphere perturbed by the atmospheric drag uncertainty but ignoring the lifting force can be described by the following equations with positionr and velocity v with their corresponding initial conditions rt0 and vt0 [5].
(1)ṙt=vt,v̇t=aDt+aGt−2ω×vt−ω×ω×rt+Ree^3+ξt,where r=x1,x2,x3T and v=v1,v2,v3T denote the position and velocity vectors, respectively, with respect to the ENU coordinate frame. Re=6378×103m is the Earth’s radius, e3=0,0,1T is the third standard unit vector in ℛ3 and ξ is the random acceleration vector due to modeling uncertainties and disturbances, and it is assumed that the noise vector is a zero-mean Gaussian process, EξtξtT=δt−sQt, where Qt is the positive definite process noise covariance matrix. ω=ω1,ω2,ω3T denotes the angular velocity vector of the ENU coordinates with respect to the ECEI (Earth-centered-Earth-fixed inertial) frame and each term of the angular velocity is computed by θ0,ϕ0 if the origin of the ENU topocentric horizon coordinate system is given ω1,ω2,ω3T=0,ωecosϕ0,ωesinϕ0T, where ωe=72.9217×10−6rad/s is the Earth’s rotation rate. aG denotes the gravitational acceleration, aG=ge^3, and based on the inverse square gravity model, the gravitational acceleration is calculated by
(2)g=geR2Re+x32,where ge denotes the gravitational acceleration at zero altitude. aD is the instantaneous acceleration due to the atmospheric density and it is assumed to be opposed to the direction of motion and proportional to the atmospheric density ρ and the velocity squared [24].
(3)aDt=−12βρvreltνrelt,β=msCdA,where β is the ballistic drag coefficient, Cd is the aerodynamic drag coefficient, A is the cross-sectional area of the satellite in a plane normal to the relative velocity vector, and ms is the mass of the reentering space object. vrel=v−w denotes the relative air velocity vector, where w is the wind velocity vector, and vrel=vrel=vrelTvrel is the magnitude of the relative air velocity vector. For the atmospheric density, an analytical exponential atmosphere model is utilized as
(4)ρx3=ρeexp−x3H,where ρe=1.752kg/m3 and H=6.7×103m.Now, after defining an augmented state vectorxt≡rTtvTtT, then the nonlinear differential equation could be written by
(5)ẋt=fxt+Gtηt,where
(6)fxt≡vtaDt+aG−2ω×vt−ω×ω×rt+Ree^3,and ηt≡03×1TξTtT is the augmented process noise vector, and Gt=03×3I3×3T.


## 2.3. State-Space Nonlinear Equations of Motion

It is assumed that the initial nominal state of a reentering space vehicle right before the breakup is given by propagating an orbital decay reentry trajectory until 78 km in ECEF (Earth-centered-Earth-fixed). After the coordinate transformation from the ECEF to ENU frame, the initial position of a breakup event is given by geometric altitudez0, longitude θ0, and latitude ϕ0. The initial position θ0,ϕ0 which locates the origin of the ENU frame on the Earth’s surface is given by
(7)x1t0,x2t0,x3t0=0,0,z0.Then, the nonlinear translational equations of motion in (1) can be rewritten in the form of a state-space representation [5] for numerical integration as follows:
(8)ẋ1t=v1t,ẋ2t=v2t,ẋ3t=v3t,v̇1t=−ρ2βvreltv1t−w1−2ω2v3t−ω3v2t+ωe2x1t+ξ1t,v̇2t=−ρ2βvreltv2t−w2−2ω3v1t−ω2ω3Re+x3t+ω32x2t+ξ2t,v̇3t=−ρ2βvreltv3t−w3−gx3+2ω2v1t−ω2ω3x2t+ω22x3t+Re+ζ3t,where w=w1,w2,w3T is the wind velocity vector.


## 3. Algorithms for Estimation of Debris Dispersion

In this section, two new approaches for the estimation of debris dispersion due to a space vehicle breakup are proposed. The estimation of debris dispersion is represented in terms of a covariance propagation which describes probabilistic ellipsoids from a nominal reentry trajectory. First, a statistical linear regression using the unscented transformation is introduced to compute the time evolution of the covariance information. Second, a Gaussian moment-matching method which calculates the first and second moments of the probabilistic distribution of the debris dispersion by using a Gauss-Hermite cubature numerical approximation is proposed. Figure3 depicts the difference between the conventional Lyapunov-based covariance approach and the unscented transformation-based covariance computation. In this work, detailed derivation for the covariance computation using the Lyapunov equation is not included but it can be referred to [5].Figure 3
(a) Lyapunov equation-based covariance propagation. (b) Statistical linear regression-based covariance propagation method—unscented transformation approach.
(a)
(b)


### 3.1. Statistical Linear Regression Approach for Estimation of Debris Dispersion

In general, statistical linear error propagation is more accurate than the error propagation by first-order Taylor-series approximation. For a statistically linearized error propagation, consider a nonlinear mapping with a set of weighted sigma pointsXi,Wi,i=1,…,r that are selected to line on the principle component axes of the state covariance with a weight [18].
(9)Zi=hXi.Then, the first- and second-order statistics of the transformed points, mean, and covariance, are computed from the sigma pointXi and transformed sigma point Zi:
(10)x^≜∑i=1rWiXi,PXX≜∑i=1rWiXi−x^Xi−x^T,PXZ≜∑i=1rWiXi−x^Zi−z^T.Before applying a statistical linear regression method, it is necessary to transform the continuous nonlinear differential equation in (6) into a discrete-time nonlinear difference equation using an approximate method, such as the Euler method, or numerical Runge-Kutta approach. In this section, it is assumed that a nonlinear discrete-time equation is given in the form [11]
(11)xk=fkxk,ζk,tk,where xk∈ℜn is the n×1 state vector. It is assumed that the noise vector ζk∈ℜq is the q×1 state noise vector and is a zero-mean Gaussian process satisfying
(12)EζkζjT=δkjQk,∀k,j,where Qk is the process noise covariance matrix. In this paper, we utilize the unscented transformation mapping technique to derive the first- and second-order moment statistics. First, in order to generate a set of weighted sigma points, the state vector xk∈ℜn is redefined as an augmented state vector xka∈ℜn+q along with noise variables and the corresponding augmented covariance matrix on the diagonal is reconstructed by
(13)xka=xkζk,Pka=Pk00Qk,where xka is the augmented state vector with the dimension na=n+q. Then, the set of scaled symmetric sigma points Xi,kai=02na=Xi,kxTXi,kςTTi=02na with the augmented state vector and covariance matrix is constructed by [11].
(14)X0,ka=x^ka,Xi,ka=x^ka+na+λPkai,i=1,…,na,Xi,ka=x^ka−na+λPkai,i=na+1,…,2na,where λ=α2na+κ−na is a scaling parameter with the constant parameter 0≤α≤1 and κ providing an extra degree of freedom for fine-tuning the higher-order moment (κ=3−na for a Gaussian distribution). The term na+λPkai, is the ith column or row vector of the weighted square root of the scaled covariance matrix na+λPka. The corresponding weights for the mean and covariance are defined by
(15)Wim=λna+λ,i=0,Wim=12na+λ,i=1,…,2na,Wic=λna+λ+1−α2+β,i=0,Wic=12na+λ,i=1,…,2na,where β is a third parameter for incorporating extra higher-order effects, and β=2 is the optimal value for Gaussian distribution.As for the prediction step for the next timek+1, the predicted state estimate vector x^k+1∣k and its predicted covariance matrix Pk+1∣k are computed from the propagated sigma point vectors as
(16)Xk+1x=fXkx,Xkw,k,(17)x^k+1∣k=∑i=02naWimfXkx,Xkw,k,(18)Pk+1∣k=∑i=02naWicXi,k+1x−x^k+1∣kXi,k+1x−x^k+1∣kT,where Xi,kx is a weighted sigma point vector of the first n elements of the ith augmented sigma point vector Xi,ka and Xi,kw is a weighted sigma point vector of the next q elements of Xi,ka, respectively.


### 3.2. Gaussian Moment-Matching Approach for Estimation of Debris Dispersion

Even though the statistical linear regression (SLR) using the unscented transformation approach could represent the propagation of the debris dispersion due to a breakup event, the SLR-based approach estimates the debris dispersion with the second-moment covariance using a finite small number of sigma points. The statistical regression approach is not a truly global approximation and it does not work well with nearly singular covariances, that is, nearly deterministic systems with small covariances. For compensating the drawbacks, it is necessary to consider a more realistic approach which can precisely capture the initial debris dispersion and also take into account the neglected nonlinearities from the statistical linear regression. In this section, the Gaussian cubature transformation approach [10] is proposed to estimate the debris dispersion with a moment-matching technique by using a set of transformed sigma lattice points which represent breakup components in the debris dispersion (shown in Figure 4).Figure 4
In Gauss-Hermite cubature integration rule, the lattice points in dimensionn=1,2,3 are required to perform p product rule integral approximation, that is, pn cubature points. The color marks the weight of the cubature points [19].Consider the transformation of the statex and the transformed random variable y=gx+q where x~Nm,P and q~N0,Q, then the moment-matching-based Gaussian approximation to the joint distribution of x and y is represented as [19]
(19)xy~Nmμ^M,PCMCMTSM,where
(20)μ^M=∫gxNx∣m,Pdx,SM=∫gx−μ^Mgx−μ^MTNx∣m,Pdx+Q,CM=∫x−mgx−μ^MTNx∣m,Pdx,and x is assumed to be Gaussian with the distribution, x∼Nm,P, m is the mean, and P is the covariance.The debris dispersion can be modeled with the firstμ^M and second moments SM of the original distribution of x and the transformed distribution y=gx+q, and this leads to the problem of forming Gaussian approximation to (x, y) by directly approximating the integrals. The advantage of the moment-matching formation lies in that it enables usage of well-known numerical integration methods such as Gauss-Hermite quadrature [20], cubature rules [19], and central difference methods [13, 14]. In this section, the integrals in (20) can be evaluated practically by the Gauss-Hermite cubature-based numerical integration method as shown in Figures 4 and 5. By generalizing the change of variable idea, we can form approximation to multidimensional integrals of the form in (19). First, it is needed to decompose the covariance, P=PPT, where P is the Cholesky factor of the covariance. Now, a new integration variable is defined by [19]
(21)x=m+Pσ,where σ=σi1,…,in is a n-dimensional vector with one-dimensional unit sigma point σik at element k. Then, the Gaussian integral in (20) can be rewritten as
(22)∫gxNx∣m,Pdx=∫gm+PσNσ∣0,Idσ.Figure 5
Illustration of a fifth-order Gauss-Hermite cubature-based approximation to a nonlinear transformation. The covariance of the true distribution is presented by the dashed line and the solid line is the approximation [10].
(a)
Original
(b)
TransformedThe integral in (22) is formulated in terms of the multidimensional unit Gaussian Ng∣0,I and can be written as an iterated integral over one-dimensional Gaussian distributions. Each of the one-dimensional integrals can be approximated with the Gauss-Hermite quadrature as
(23)∫gm+PσNσ∣0,Idσ=∫⋯∫gm+PσNσ1∣0,1dσ1×⋯×Nσn∣0,1dσn≈∑i1,…,inWi1×⋯×Wingm+Pσi1,…,in,where Wik, k=1,…,n are simply the corresponding one-dimensional Gauss-Hermite weights with pth-order approximation, which can be calculated by
(24)Wi=p!p2Hp−1σi2,and Hpx is a Hermite polynomial of order p given by
(25)Hpx=−1pexpx22dpdxpexp−x22.The unit sigma pointsσi, i=1,…,p is calculated by finding the roots of the Hermite polynomial Hpx. The extension of the Gauss-Hermite quadrature rule to an n-dimensional cubature rule by using the product rule lattice approach yields a rather good numerical integration method. Based on the Gauss-Hermite cubature rule, the multidimensional weights are calculated as the products of one-dimensional weights
(26)Wi1,…,in=Wi1×⋯×Win=p!p2Hp−1σi12×⋯×p!p2Hp−1σin2.Also, the multidimensional unit sigma points can be given as Cartesian product of the one-dimensional unit sigma points(27)σi1,…,in=σi1⋮σin.It is noted that the number of sigma points required forn-dimensional integral with pth-order rule is pn, which could become unfeasible when the number of dimensions grows which is indicated in Figure 4.Based on the previous Gauss-Hermite cubature integration approach, an additive form of the multidimensional Gauss-Hermite cubature integral technique is used to represent the covariance prediction of a debris dispersion. It is assumed that the posterior density functionpn is known. First, the unit sigma point σi, i=1,…,p is calculated by finding the roots of the Hermite polynomial Hpx. Based on the Gauss-Hermite cubature rule, the multidimensional weights are calculated as the products of one-dimensional weights:
(28)Wi1,…,in=Wi1×⋯×Win=p!p2Hp−1σi12×⋯×p!p2Hp−1σin2,where one-dimensional weight Wik is given by using (26). Now, the sigma points can be generated by
(29)Xk∣ki=mk∣k+Pk∣kσi,i=1,…,pn,where mk∣k is the mean vector of the current debris state and Pk∣k is the covariance matrix representing the debris dispersion at the event of the breakup time, k, and the unit sigma points σi1,…,in are defined in (27). Next, the sigma points can be propagated by using the nonlinear discrete equation in (11) as
(30)Xk+1∣ki=fXk∣ki,k,i=1,…,pn.As for the prediction step for the next timek + 1, the predicted state vector mk+1∣k of the debris trajectory and its predicted covariance matrix Pk+1∣k representing the prediction of the debris dispersion are computed by using the propagated sigma point vectors in (30) as
(31)mk+1∣k=∑i=1pnWiXk+1∣ki,(32)Pk+1∣k=∑i=1pnWiXk+1∣ki−mk+1∣kXk+1∣ki−mk+1∣kT+Qk+1,where the weight Wi is defined in (28) and Qk+1 is the process noise covariance representing the uncertainties of truncated modeling errors and unknown errors.The multidimensional Gaussian-Hermite cubature-based covariance prediction approach which represents the debris dispersion and the distribution in time can be interpreted as a special form of a Monte-Carlo integration approach.


## 3.1. Statistical Linear Regression Approach for Estimation of Debris Dispersion

In general, statistical linear error propagation is more accurate than the error propagation by first-order Taylor-series approximation. For a statistically linearized error propagation, consider a nonlinear mapping with a set of weighted sigma pointsXi,Wi,i=1,…,r that are selected to line on the principle component axes of the state covariance with a weight [18].
(9)Zi=hXi.Then, the first- and second-order statistics of the transformed points, mean, and covariance, are computed from the sigma pointXi and transformed sigma point Zi:
(10)x^≜∑i=1rWiXi,PXX≜∑i=1rWiXi−x^Xi−x^T,PXZ≜∑i=1rWiXi−x^Zi−z^T.Before applying a statistical linear regression method, it is necessary to transform the continuous nonlinear differential equation in (6) into a discrete-time nonlinear difference equation using an approximate method, such as the Euler method, or numerical Runge-Kutta approach. In this section, it is assumed that a nonlinear discrete-time equation is given in the form [11]
(11)xk=fkxk,ζk,tk,where xk∈ℜn is the n×1 state vector. It is assumed that the noise vector ζk∈ℜq is the q×1 state noise vector and is a zero-mean Gaussian process satisfying
(12)EζkζjT=δkjQk,∀k,j,where Qk is the process noise covariance matrix. In this paper, we utilize the unscented transformation mapping technique to derive the first- and second-order moment statistics. First, in order to generate a set of weighted sigma points, the state vector xk∈ℜn is redefined as an augmented state vector xka∈ℜn+q along with noise variables and the corresponding augmented covariance matrix on the diagonal is reconstructed by
(13)xka=xkζk,Pka=Pk00Qk,where xka is the augmented state vector with the dimension na=n+q. Then, the set of scaled symmetric sigma points Xi,kai=02na=Xi,kxTXi,kςTTi=02na with the augmented state vector and covariance matrix is constructed by [11].
(14)X0,ka=x^ka,Xi,ka=x^ka+na+λPkai,i=1,…,na,Xi,ka=x^ka−na+λPkai,i=na+1,…,2na,where λ=α2na+κ−na is a scaling parameter with the constant parameter 0≤α≤1 and κ providing an extra degree of freedom for fine-tuning the higher-order moment (κ=3−na for a Gaussian distribution). The term na+λPkai, is the ith column or row vector of the weighted square root of the scaled covariance matrix na+λPka. The corresponding weights for the mean and covariance are defined by
(15)Wim=λna+λ,i=0,Wim=12na+λ,i=1,…,2na,Wic=λna+λ+1−α2+β,i=0,Wic=12na+λ,i=1,…,2na,where β is a third parameter for incorporating extra higher-order effects, and β=2 is the optimal value for Gaussian distribution.As for the prediction step for the next timek+1, the predicted state estimate vector x^k+1∣k and its predicted covariance matrix Pk+1∣k are computed from the propagated sigma point vectors as
(16)Xk+1x=fXkx,Xkw,k,(17)x^k+1∣k=∑i=02naWimfXkx,Xkw,k,(18)Pk+1∣k=∑i=02naWicXi,k+1x−x^k+1∣kXi,k+1x−x^k+1∣kT,where Xi,kx is a weighted sigma point vector of the first n elements of the ith augmented sigma point vector Xi,ka and Xi,kw is a weighted sigma point vector of the next q elements of Xi,ka, respectively.


## 3.2. Gaussian Moment-Matching Approach for Estimation of Debris Dispersion

Even though the statistical linear regression (SLR) using the unscented transformation approach could represent the propagation of the debris dispersion due to a breakup event, the SLR-based approach estimates the debris dispersion with the second-moment covariance using a finite small number of sigma points. The statistical regression approach is not a truly global approximation and it does not work well with nearly singular covariances, that is, nearly deterministic systems with small covariances. For compensating the drawbacks, it is necessary to consider a more realistic approach which can precisely capture the initial debris dispersion and also take into account the neglected nonlinearities from the statistical linear regression. In this section, the Gaussian cubature transformation approach [10] is proposed to estimate the debris dispersion with a moment-matching technique by using a set of transformed sigma lattice points which represent breakup components in the debris dispersion (shown in Figure 4).Figure 4
In Gauss-Hermite cubature integration rule, the lattice points in dimensionn=1,2,3 are required to perform p product rule integral approximation, that is, pn cubature points. The color marks the weight of the cubature points [19].Consider the transformation of the statex and the transformed random variable y=gx+q where x~Nm,P and q~N0,Q, then the moment-matching-based Gaussian approximation to the joint distribution of x and y is represented as [19]
(19)xy~Nmμ^M,PCMCMTSM,where
(20)μ^M=∫gxNx∣m,Pdx,SM=∫gx−μ^Mgx−μ^MTNx∣m,Pdx+Q,CM=∫x−mgx−μ^MTNx∣m,Pdx,and x is assumed to be Gaussian with the distribution, x∼Nm,P, m is the mean, and P is the covariance.The debris dispersion can be modeled with the firstμ^M and second moments SM of the original distribution of x and the transformed distribution y=gx+q, and this leads to the problem of forming Gaussian approximation to (x, y) by directly approximating the integrals. The advantage of the moment-matching formation lies in that it enables usage of well-known numerical integration methods such as Gauss-Hermite quadrature [20], cubature rules [19], and central difference methods [13, 14]. In this section, the integrals in (20) can be evaluated practically by the Gauss-Hermite cubature-based numerical integration method as shown in Figures 4 and 5. By generalizing the change of variable idea, we can form approximation to multidimensional integrals of the form in (19). First, it is needed to decompose the covariance, P=PPT, where P is the Cholesky factor of the covariance. Now, a new integration variable is defined by [19]
(21)x=m+Pσ,where σ=σi1,…,in is a n-dimensional vector with one-dimensional unit sigma point σik at element k. Then, the Gaussian integral in (20) can be rewritten as
(22)∫gxNx∣m,Pdx=∫gm+PσNσ∣0,Idσ.Figure 5
Illustration of a fifth-order Gauss-Hermite cubature-based approximation to a nonlinear transformation. The covariance of the true distribution is presented by the dashed line and the solid line is the approximation [10].
(a)
Original
(b)
TransformedThe integral in (22) is formulated in terms of the multidimensional unit Gaussian Ng∣0,I and can be written as an iterated integral over one-dimensional Gaussian distributions. Each of the one-dimensional integrals can be approximated with the Gauss-Hermite quadrature as
(23)∫gm+PσNσ∣0,Idσ=∫⋯∫gm+PσNσ1∣0,1dσ1×⋯×Nσn∣0,1dσn≈∑i1,…,inWi1×⋯×Wingm+Pσi1,…,in,where Wik, k=1,…,n are simply the corresponding one-dimensional Gauss-Hermite weights with pth-order approximation, which can be calculated by
(24)Wi=p!p2Hp−1σi2,and Hpx is a Hermite polynomial of order p given by
(25)Hpx=−1pexpx22dpdxpexp−x22.The unit sigma pointsσi, i=1,…,p is calculated by finding the roots of the Hermite polynomial Hpx. The extension of the Gauss-Hermite quadrature rule to an n-dimensional cubature rule by using the product rule lattice approach yields a rather good numerical integration method. Based on the Gauss-Hermite cubature rule, the multidimensional weights are calculated as the products of one-dimensional weights
(26)Wi1,…,in=Wi1×⋯×Win=p!p2Hp−1σi12×⋯×p!p2Hp−1σin2.Also, the multidimensional unit sigma points can be given as Cartesian product of the one-dimensional unit sigma points(27)σi1,…,in=σi1⋮σin.It is noted that the number of sigma points required forn-dimensional integral with pth-order rule is pn, which could become unfeasible when the number of dimensions grows which is indicated in Figure 4.Based on the previous Gauss-Hermite cubature integration approach, an additive form of the multidimensional Gauss-Hermite cubature integral technique is used to represent the covariance prediction of a debris dispersion. It is assumed that the posterior density functionpn is known. First, the unit sigma point σi, i=1,…,p is calculated by finding the roots of the Hermite polynomial Hpx. Based on the Gauss-Hermite cubature rule, the multidimensional weights are calculated as the products of one-dimensional weights:
(28)Wi1,…,in=Wi1×⋯×Win=p!p2Hp−1σi12×⋯×p!p2Hp−1σin2,where one-dimensional weight Wik is given by using (26). Now, the sigma points can be generated by
(29)Xk∣ki=mk∣k+Pk∣kσi,i=1,…,pn,where mk∣k is the mean vector of the current debris state and Pk∣k is the covariance matrix representing the debris dispersion at the event of the breakup time, k, and the unit sigma points σi1,…,in are defined in (27). Next, the sigma points can be propagated by using the nonlinear discrete equation in (11) as
(30)Xk+1∣ki=fXk∣ki,k,i=1,…,pn.As for the prediction step for the next timek + 1, the predicted state vector mk+1∣k of the debris trajectory and its predicted covariance matrix Pk+1∣k representing the prediction of the debris dispersion are computed by using the propagated sigma point vectors in (30) as
(31)mk+1∣k=∑i=1pnWiXk+1∣ki,(32)Pk+1∣k=∑i=1pnWiXk+1∣ki−mk+1∣kXk+1∣ki−mk+1∣kT+Qk+1,where the weight Wi is defined in (28) and Qk+1 is the process noise covariance representing the uncertainties of truncated modeling errors and unknown errors.The multidimensional Gaussian-Hermite cubature-based covariance prediction approach which represents the debris dispersion and the distribution in time can be interpreted as a special form of a Monte-Carlo integration approach.


## 4. Simulation Results

In this section, first, various parametric studies are carried out to analyze the effects of uncertainties in the atmospheric density and drag coefficient along with wind effects. Then, the performance of the proposed estimation techniques are investigated by comparing two covariance propagation methods; the unscented transformation approach and the Gauss-Hermite cubature integration method.


### 4.1. Parametric Study for Analysis of Effects of Modeling Uncertainties

Atmospheric drag most strongly influences the motion of reentry objects near the Earth. When it comes to determining accurately atmospheric drag, the values of atmospheric density and drag coefficient are the main factors to the trajectory in the sense that the information on the specific shape and quantity of objects are generally unknown. In light of this reason, the different values of the drag coefficient are selected to see the how the drag coefficient affects the reentry point and the probabilistic ellipsoid of debris dispersion caused by the breakup event in this study. In this study, four different atmospheric density models including CIRA-72, Exponential model, U.S. standard model, and NRLMSISE-00 are used to analyze the effects of the atmospheric density model onto the reentry trajectory. Figure6 shows density profiles corresponding to each density model as a function of altitude. Since wind also influences the reentry trajectory near to the ground, it is necessary to design a mathematical wind profile model. In this study, HWM-07 (Horizontal Wind Model 2007) [24] is employed to consider a more accurate reentry model. Figure 7 shows the velocity of eastward and westward winds with respect to the history of altitude.Figure 6
Density profiles for density models.Figure 7
Wind profiles for different altitudes.
(a)
(b)In order to analyze the effects of the density model and the drag coefficient, two simulation examples are investigated. The first example considers the lifetime prediction of uncontrolled space objects with different atmospheric density models and drag coefficients. The satellite under consideration has the following orbit parameters:a=6678km, e=1.0×10−5, i=28.5°, ω=0°, and Ω=0°, and epoch time is defined as August 1, 2015. Figure 8 shows that the orbital decay of the space object until 78 km is different with respect to atmospheric models. Using the U.S. standard model, the decay of altitude is faster than any other density models. On the other hand, the NRLSMSISE-00 model predicts the decay of altitude of the reentry object with the slowest drop among the three models.Figure 8
Decay of altitude with different atmospheric models.Table1 presents the reentry latitude as well as longitude of different atmospheric models at the 78 km altitude and reentry time. Compared to the STK LTP (lifetime prediction) tool [26], it is shown that the results of reentry time are similar. Figure 9 indicates the orbital decay of space objects with a different drag coefficient. As the drag coefficient increases, it is shown that the space object falls faster.Table 1
Terminal data of lifetime prediction with different density models.
Density model
Reentry latitude
Reentry longitude
Reentry time
STK LTP tool
CIRA-72
−38.936°
−345.39°
Aug. 26 (01:15)
Aug. 25 (18:21)
U.S. standard
−27.306°
−224.34°
Aug. 28 (15:44)
Aug. 27 (09:02)
NRLMSISE-00
−8.302°
−319.91°
Aug. 20 (05:31)
Aug. 19 (03:22)Figure 9
Decay of altitude with different drag coefficients.On the other hand, a second example is illustrated to show the effects of different density models and drag coefficient onto the debris nominal trajectory as well as the time evolution of the debris dispersion from 78 km at the breakup event to impact ground instance. For the simulation study, the number of 200 samples is drawn to model the initial debris breakup dispersion with the assumption of a Gaussian distribution. It is assumed that the initial breakup state is given byxt0=0078,000m709.9m/s0m/s−124.0m/s with the breakup altitude at 78 km. The total simulation time was about 32 min. An initial breakup dispersion is generated with an incremental speed Δv uniformly added to the nominal velocity vt0 in every direction by creating a spherical grid of 32 × 32 = 1024 points, that is, vi,j=vt0+Δvi,j [7]. So, with each 32 points for longitude and latitude, respectively, initial velocities in the initial debris dispersion are created by
(33)vi,j=vt0+Δvcosϕicosθjcosϕisinθjsinϕi,i,j=1,2,…,32.The values of each latitude and longitude are taken withϕi∈−π/2,π/2 and θj∈0,2π, respectively, where i,j=1,2,…,32 and Δv=100m/s is an incremental speed.Figures10 and 11 show the positional ellipsoid trajectory of breakup dispersion with different atmospheric density models and drag coefficients. As can be seen, even though each case of breakup event occurred at the same position, the final location of the ground impact point is very different depending on the type of the density model and the drag coefficient. Figures 12 and 13 show the final location and size of the positional ellipsoid on the ground in detail. For the case of different atmospheric models, the CIRA-72 model generated the largest impact area of dispersion debris, while the NRLSMSISE-00 model showed a smaller one, which indicates that the CIRA-72 model is subject to a lot of uncertainty errors compared to that of the NRLMSISE-00 models. Table 2 describes the footprint statistics including impact areas. For the case of a different drag coefficient, the size of the final positional ellipsoid of debris dispersion becomes smaller as the values of drag coefficients increase. It is indicated that a higher drag coefficient value leads to a faster fall to the ground with a short evolution. As can be shown in Table 3, the positional ellipsoidal area of CD=1.5 is almost twice larger than that of CD=2.5. It is seen that a properly estimated drag coefficient is one of the main factors in the propagation of the probabilistic positional ellipsoid.Figure 10
Time evolution of probability ellipsoid with different density models.Figure 11
Time evolution of probability ellipsoid with different drag coefficients.Figure 12
Impact footprint (density models).Figure 13
Impact footprint (drag coefficients).Table 2
Footprint statistics for various atmospheric density models.
Density model
Impact latitude
Impact longitude
Impact area
CIRA-72
−3.007°
−240.2°
1522 km2
U.S. standard
−27.306°
−224.34°
619.2 km2
NRLMSISE-00
−8.302°
−319.91°
913.8 km2Table 3
Footprint statistics for various drag coefficients.
Drag coefficient
Impact latitude
Impact longitude
Impact area
CD = 1.5
12.23°
−137.3°
1093 km2
CD = 2.0
−7.452°
−347.3°
916.5 km2
CD = 2.5
−11.96°
−118.3°
547.3 km2


### 4.2. Statistical Computation of Debris Dispersion with Covariance Estimation

In this section, a detailed analysis of the estimation of debris dispersion due to the breakup event during the reentry is made by using one of the proposed covariance propagation methods, the unscented transformation technique. Note that the positional variation at the breakup instant is small enough to be neglected but the variation of the debris velocity becomes relatively big. Therefore, it is assumed that an initial breakup dispersion could be generated by adding an incremental speedδvi,j to the initial nominal velocity vt0 in every direction, vi,j=vt0+δvi,j like in (33). Now, based on the initial breakup model with an incremental velocity dispersion as above, an initial debris dispersion at the event of the breakup time t0 is remodeled with a Gaussian distribution before covariance propagation for the time evolution of the debris dispersion in (34).
(34)Pt0=106diag10−1010−1010−100.0025m/s20.0025m/s20.0053m/s2.Note that the initial breakup dispersion generated with an incremental speedΔv uniformly added to the nominal velocity vt0 in every direction could enter the inside of the initial covariance and thus the initial debris dispersion at the breakup event is simulated as shown in Figure 14.Figure 14
Initial estimate of debris dispersion.Since the probabilistic distribution of the debris dispersion is Gaussian, the center of the ellipsoid becomes the breakup point and the dispersion is represented by an ellipsoid. Uncertainty due to unknown and neglected acceleration terms is characterized by using the process noise covariance, that is,EηtηςT≡Qtδτ−ς in (12).
(35)Qt=106diag00010−1610−1610−16.It is assumed that the breakup altitude is 78 km and the initial breakup position is located atxt0=0078,000m709.9m/s0m/s−124.0m/s as seen in Figure 15. The total simulation time was about 32 minutes. Based on the initial covariance in (34), the process noise covariance matrix in (35), and other initial conditions, the initial covariance estimate of the debris dispersion could be propagated in time by using the one of the covariance propagation techniques such as the Gauss-Hermite cubature or the unscented transformation approach. Before propagating the covariance of the debris dispersion, the unscented transformation technique needs to generate a set of sigma points which also represents a set of initial debris fragments at the breakup event. Now, a set of scaled symmetric sigma points Xi,kai=02na=Xi,kxTXi,kςTTi=02na with the augmented state vector and covariance matrix is constructed by
(36)X0,0a=x^0a,Xi,0a=x^0a+na+λP0ai,i=1,…,na,Xi,0a=x^0a−na+λP0ai,i=na+1,…,2na,where the term, na+λP0ai, is the ith column or row vector of the weighted square root of the scaled covariance matrix na+λP0a, and the augmented initial covariance is given by
(37)x0a=x0ζ0,P0a=Pt000QtΔt.Figure 15
Breakup point and impact footprint.After the generation of the initial sigma points, the covariance is propagated by using (16), (17), and (18).Figure16 depicts the 3-dimensional positional probability ellipsoids corresponding to a confidence interval of 99.99% with a sequence of time instants. The breakup event during the reentry happened at the altitude of 78 km, and the predicted trajectory and the positional covariance ellipsoid of the debris fragments are well depicted with the usage of the proposed breakup dispersion estimation method. As can be shown, the positional probability ellipsoid representing the debris dispersion increases quickly within 2.4 minutes, and after 5 minutes the ellipsoid increases gradually until it impacts the ground. This phenomenon is understood by the fact that after 3 minutes from the breakup event the motion becomes a sharp fall and becomes nearly constant in motion.Figure 16
Dispersion probability ellipsoids for a sequence of time and impact footprint on the ground.This phenomenon is verified in Figure17 which shows the plots of the standard deviation of the positional covariance matrices in each direction over time. As can be expected from the time evolution of the debris dispersion in Figure 16, the variance sharply increases within 3 minutes but gradually changes after that time, that is, the falling motion 3 minutes after the breakup event.Figure 17
Magnitude of positional covariances in each direction (x, y, and z).For a better analysis of the first few minutes, a detailed view is illustrated in Figure18 which depicts the change of the initial debris dispersion in terms of the ellipsoid from the breakup event to 100 seconds. The covariance ellipsoid contains the debris fragments which are broken into all directions with the magnitude of the initial covariance matrix. For the first 100 seconds, the dispersion grows fast and becomes 10 times larger than the initial breakup ellipsoid. This indicates that right after the breakup event there are lots of uncertainties acting on the motion of the instantaneous debris dispersion and also the motion is highly nonlinear. A commonly used metric for covariance matrices is the square root of the trace of the covariance ellipsoid as seen in Figure 17.Figure 18
Initial dispersion of ellipsoids after the breakup event within 100 seconds.The simulation is made with different ballistic coefficients using the empirical density and wind models. Figure19 shows the plot of flight time versus ballistic coefficient. In Figure 19, for the variance of the ballistic coefficient in terms of altitude, an empirical density profile obtained from the MSISE-00 density model [22] is used. In addition, in order to include the effects of wind near the ground, a wind model is developed by using a piecewise cubic spline method through a curve fitting of real wind data. Figure 20 represents the wind profiles used in the simulation in the eastern and northern direction as a function of different altitudes. The wind model used in this work is based on HWM-07 (Horizon Wind Model 2007) [7, 24]. Figure 21 illustrates the effects of the wind strength onto the time variation of the dispersion trajectory. As can be seen, if the strength of the wind is stronger, then it takes a longer time until the debris arrives the ground and impact it.Figure 19
Plot of the empirical ballistic coefficient as a function of flight time [7].Figure 20
Simulated eastern and northern wind profiles for different altitudes.
(a)
(b)Figure 21
Effects of the wind to the time evolution of the decay of the debris altitude.


### 4.3. Performance Comparison for Estimation of Debris Dispersion

In this section, the performance of the proposed covariance methods for estimating the debris dispersion due to the reentry breakup is investigated in terms of the time evolution of the positional probabilistic ellipsoid. The covariance propagation methods include the unscented transformation-based statistical linear regression and the Gauss-Hermite cubature-based numerical integration methods described in Section3.For the simulation test, it is also assumed that the breakup altitude is 78 km and the initial breakup position is located atxt0=0078,000m709.9m/s0m/s−124.0m/s. The initial breakup dispersion is modeled as a Gaussian distribution with the initial covariance ellipsoid as given in (34); the process noise covariance matrix representing modeling uncertainties is also given by Qt=106diag00010−1610−1610−16 in (35). Now, based on the initial debris dispersion covariance matrix Pt0 and the process noise matrix Qt, a set of sigma points representing a set of breakup components can be generated by using (14) for the unscented transformation (UT) method and using (29) for the Gauss-Hermite cubature (GHC) method. After the generation of the initial sigma points for the initial debris dispersion, the covariance is propagated by using (16), (17), and (18) for the UT method and (30), (31), and (32) for the GHC approach. The total simulation time used was 32 minutes.Figure22 depicts the estimated nominal trajectory from the breakup event to the ground impact instant, and also the footprint of the ground impact location. After the breakup event of the space debris during reentry, it is assumed that the simulation starts at t=0 at the event, and then the nominal trajectory in terms of the mean value of all the debris fragments and the 3-dimensional positional ellipsoid are propagated until the debris fragments reach the ground. Based on the result in Figure 22, it is seen that the time evolution of the nominal estimates of the debris dispersion have different final locations but the overall trajectories are similar to each other. The nominal trajectory generated by the Gauss-Hermite cubature method has a little bit faster fall into the ground. The locational difference between the final footprint impact locations of two different methods was about 2 km in the x direction, and 10 m in the y direction, that is, δxtfinal=1.8980km0.00860km0km. It is seen from the nominal trajectory plot that the uncertainties of the motion of the debris dispersion are captured more precisely by utilizing the Gauss-Hermite cubature technique compared to the unscented transformation method, and this is because the GHC method adopts a much larger set of sigma lattice points at the initial debris dispersion which leads to a more precise estimation of the debris dispersion. In the simulation, it is also seen that the reentry motion is a nonlinear motion right after the breakup event, but from the altitude near the ground the motion has relatively small nonlinearities due to a vertical falling motion.Figure 22
Performance comparison between unscented transformation and Gauss-Hermite cubature methods in nominal trajectory prediction and ground impact footprint.Figures23 and 24 depict the profiles of the positional probability ellipsoids corresponding to a confidence interval of 99.99% for a sequence of time instants generated from the proposed covariance propagation methods, statistical linear regression, and the Gauss-Hermite moment-matching methods, respectively. In general, since the Gauss-Hermite cubature approach could include more debris breakup distribution points with p=36=729 sigma points, it could generate a more accurate trajectory of the core debris as well as the covariance estimate of the breakup dispersion, while the unscented transformation method generates the initial breakup components with p=2×6+1=13 sigma points. In addition, the performance of the proposed approach is investigated by checking the time evolution of the dispersion ellipsoid shown in Figures 23 and 24 where the time evolution of the positional probability ellipsoids can be seen to almost match between the statistical linear regression in Figure 23 and the Gauss-Hermite cubature-based numerical integration shown in Figure 24. In fact, the close agreement of these techniques can be justified by comparing the positional covariance matrices from both of the methods in each direction. It is seen that the size of the overall debris dispersion is similar to each other, 0.6 km in the x direction and 0.01 km in the y direction. Figure 25 shows the plots of the sigma values for the positional covariance matrices for both of the methods overtime. The difference of the sigma in the x direction is large but the other directions have similar values. The Gauss-Hermite cubature method could provide a little bit more accurate estimates of the nominal trajectory and the positional ellipsoid which could affect the calculation of the risk causality, but it also requires more computational power compared to the unscented transformation method. Both methods could satisfy the computational procedures required for real-time computation in a sequential estimation, while a Monte-Carlo-based approach uses offline computational procedures.Figure 23
Time evolution of debris dispersion ellipsoids and impact footprint on the ground using the statistical linear regression with the unscented transformation method.Figure 24
Time evolution of debris dispersion ellipsoids and impact footprint on the ground using the Gauss-Hermite cubature integration method.Figure 25
Comparison of the magnitude of positional covariances in each direction (x, y, and z).


## 4.1. Parametric Study for Analysis of Effects of Modeling Uncertainties

Atmospheric drag most strongly influences the motion of reentry objects near the Earth. When it comes to determining accurately atmospheric drag, the values of atmospheric density and drag coefficient are the main factors to the trajectory in the sense that the information on the specific shape and quantity of objects are generally unknown. In light of this reason, the different values of the drag coefficient are selected to see the how the drag coefficient affects the reentry point and the probabilistic ellipsoid of debris dispersion caused by the breakup event in this study. In this study, four different atmospheric density models including CIRA-72, Exponential model, U.S. standard model, and NRLMSISE-00 are used to analyze the effects of the atmospheric density model onto the reentry trajectory. Figure6 shows density profiles corresponding to each density model as a function of altitude. Since wind also influences the reentry trajectory near to the ground, it is necessary to design a mathematical wind profile model. In this study, HWM-07 (Horizontal Wind Model 2007) [24] is employed to consider a more accurate reentry model. Figure 7 shows the velocity of eastward and westward winds with respect to the history of altitude.Figure 6
Density profiles for density models.Figure 7
Wind profiles for different altitudes.
(a)
(b)In order to analyze the effects of the density model and the drag coefficient, two simulation examples are investigated. The first example considers the lifetime prediction of uncontrolled space objects with different atmospheric density models and drag coefficients. The satellite under consideration has the following orbit parameters:a=6678km, e=1.0×10−5, i=28.5°, ω=0°, and Ω=0°, and epoch time is defined as August 1, 2015. Figure 8 shows that the orbital decay of the space object until 78 km is different with respect to atmospheric models. Using the U.S. standard model, the decay of altitude is faster than any other density models. On the other hand, the NRLSMSISE-00 model predicts the decay of altitude of the reentry object with the slowest drop among the three models.Figure 8
Decay of altitude with different atmospheric models.Table1 presents the reentry latitude as well as longitude of different atmospheric models at the 78 km altitude and reentry time. Compared to the STK LTP (lifetime prediction) tool [26], it is shown that the results of reentry time are similar. Figure 9 indicates the orbital decay of space objects with a different drag coefficient. As the drag coefficient increases, it is shown that the space object falls faster.Table 1
Terminal data of lifetime prediction with different density models.
Density model
Reentry latitude
Reentry longitude
Reentry time
STK LTP tool
CIRA-72
−38.936°
−345.39°
Aug. 26 (01:15)
Aug. 25 (18:21)
U.S. standard
−27.306°
−224.34°
Aug. 28 (15:44)
Aug. 27 (09:02)
NRLMSISE-00
−8.302°
−319.91°
Aug. 20 (05:31)
Aug. 19 (03:22)Figure 9
Decay of altitude with different drag coefficients.On the other hand, a second example is illustrated to show the effects of different density models and drag coefficient onto the debris nominal trajectory as well as the time evolution of the debris dispersion from 78 km at the breakup event to impact ground instance. For the simulation study, the number of 200 samples is drawn to model the initial debris breakup dispersion with the assumption of a Gaussian distribution. It is assumed that the initial breakup state is given byxt0=0078,000m709.9m/s0m/s−124.0m/s with the breakup altitude at 78 km. The total simulation time was about 32 min. An initial breakup dispersion is generated with an incremental speed Δv uniformly added to the nominal velocity vt0 in every direction by creating a spherical grid of 32 × 32 = 1024 points, that is, vi,j=vt0+Δvi,j [7]. So, with each 32 points for longitude and latitude, respectively, initial velocities in the initial debris dispersion are created by
(33)vi,j=vt0+Δvcosϕicosθjcosϕisinθjsinϕi,i,j=1,2,…,32.The values of each latitude and longitude are taken withϕi∈−π/2,π/2 and θj∈0,2π, respectively, where i,j=1,2,…,32 and Δv=100m/s is an incremental speed.Figures10 and 11 show the positional ellipsoid trajectory of breakup dispersion with different atmospheric density models and drag coefficients. As can be seen, even though each case of breakup event occurred at the same position, the final location of the ground impact point is very different depending on the type of the density model and the drag coefficient. Figures 12 and 13 show the final location and size of the positional ellipsoid on the ground in detail. For the case of different atmospheric models, the CIRA-72 model generated the largest impact area of dispersion debris, while the NRLSMSISE-00 model showed a smaller one, which indicates that the CIRA-72 model is subject to a lot of uncertainty errors compared to that of the NRLMSISE-00 models. Table 2 describes the footprint statistics including impact areas. For the case of a different drag coefficient, the size of the final positional ellipsoid of debris dispersion becomes smaller as the values of drag coefficients increase. It is indicated that a higher drag coefficient value leads to a faster fall to the ground with a short evolution. As can be shown in Table 3, the positional ellipsoidal area of CD=1.5 is almost twice larger than that of CD=2.5. It is seen that a properly estimated drag coefficient is one of the main factors in the propagation of the probabilistic positional ellipsoid.Figure 10
Time evolution of probability ellipsoid with different density models.Figure 11
Time evolution of probability ellipsoid with different drag coefficients.Figure 12
Impact footprint (density models).Figure 13
Impact footprint (drag coefficients).Table 2
Footprint statistics for various atmospheric density models.
Density model
Impact latitude
Impact longitude
Impact area
CIRA-72
−3.007°
−240.2°
1522 km2
U.S. standard
−27.306°
−224.34°
619.2 km2
NRLMSISE-00
−8.302°
−319.91°
913.8 km2Table 3
Footprint statistics for various drag coefficients.
Drag coefficient
Impact latitude
Impact longitude
Impact area
CD = 1.5
12.23°
−137.3°
1093 km2
CD = 2.0
−7.452°
−347.3°
916.5 km2
CD = 2.5
−11.96°
−118.3°
547.3 km2


## 4.2. Statistical Computation of Debris Dispersion with Covariance Estimation

In this section, a detailed analysis of the estimation of debris dispersion due to the breakup event during the reentry is made by using one of the proposed covariance propagation methods, the unscented transformation technique. Note that the positional variation at the breakup instant is small enough to be neglected but the variation of the debris velocity becomes relatively big. Therefore, it is assumed that an initial breakup dispersion could be generated by adding an incremental speedδvi,j to the initial nominal velocity vt0 in every direction, vi,j=vt0+δvi,j like in (33). Now, based on the initial breakup model with an incremental velocity dispersion as above, an initial debris dispersion at the event of the breakup time t0 is remodeled with a Gaussian distribution before covariance propagation for the time evolution of the debris dispersion in (34).
(34)Pt0=106diag10−1010−1010−100.0025m/s20.0025m/s20.0053m/s2.Note that the initial breakup dispersion generated with an incremental speedΔv uniformly added to the nominal velocity vt0 in every direction could enter the inside of the initial covariance and thus the initial debris dispersion at the breakup event is simulated as shown in Figure 14.Figure 14
Initial estimate of debris dispersion.Since the probabilistic distribution of the debris dispersion is Gaussian, the center of the ellipsoid becomes the breakup point and the dispersion is represented by an ellipsoid. Uncertainty due to unknown and neglected acceleration terms is characterized by using the process noise covariance, that is,EηtηςT≡Qtδτ−ς in (12).
(35)Qt=106diag00010−1610−1610−16.It is assumed that the breakup altitude is 78 km and the initial breakup position is located atxt0=0078,000m709.9m/s0m/s−124.0m/s as seen in Figure 15. The total simulation time was about 32 minutes. Based on the initial covariance in (34), the process noise covariance matrix in (35), and other initial conditions, the initial covariance estimate of the debris dispersion could be propagated in time by using the one of the covariance propagation techniques such as the Gauss-Hermite cubature or the unscented transformation approach. Before propagating the covariance of the debris dispersion, the unscented transformation technique needs to generate a set of sigma points which also represents a set of initial debris fragments at the breakup event. Now, a set of scaled symmetric sigma points Xi,kai=02na=Xi,kxTXi,kςTTi=02na with the augmented state vector and covariance matrix is constructed by
(36)X0,0a=x^0a,Xi,0a=x^0a+na+λP0ai,i=1,…,na,Xi,0a=x^0a−na+λP0ai,i=na+1,…,2na,where the term, na+λP0ai, is the ith column or row vector of the weighted square root of the scaled covariance matrix na+λP0a, and the augmented initial covariance is given by
(37)x0a=x0ζ0,P0a=Pt000QtΔt.Figure 15
Breakup point and impact footprint.After the generation of the initial sigma points, the covariance is propagated by using (16), (17), and (18).Figure16 depicts the 3-dimensional positional probability ellipsoids corresponding to a confidence interval of 99.99% with a sequence of time instants. The breakup event during the reentry happened at the altitude of 78 km, and the predicted trajectory and the positional covariance ellipsoid of the debris fragments are well depicted with the usage of the proposed breakup dispersion estimation method. As can be shown, the positional probability ellipsoid representing the debris dispersion increases quickly within 2.4 minutes, and after 5 minutes the ellipsoid increases gradually until it impacts the ground. This phenomenon is understood by the fact that after 3 minutes from the breakup event the motion becomes a sharp fall and becomes nearly constant in motion.Figure 16
Dispersion probability ellipsoids for a sequence of time and impact footprint on the ground.This phenomenon is verified in Figure17 which shows the plots of the standard deviation of the positional covariance matrices in each direction over time. As can be expected from the time evolution of the debris dispersion in Figure 16, the variance sharply increases within 3 minutes but gradually changes after that time, that is, the falling motion 3 minutes after the breakup event.Figure 17
Magnitude of positional covariances in each direction (x, y, and z).For a better analysis of the first few minutes, a detailed view is illustrated in Figure18 which depicts the change of the initial debris dispersion in terms of the ellipsoid from the breakup event to 100 seconds. The covariance ellipsoid contains the debris fragments which are broken into all directions with the magnitude of the initial covariance matrix. For the first 100 seconds, the dispersion grows fast and becomes 10 times larger than the initial breakup ellipsoid. This indicates that right after the breakup event there are lots of uncertainties acting on the motion of the instantaneous debris dispersion and also the motion is highly nonlinear. A commonly used metric for covariance matrices is the square root of the trace of the covariance ellipsoid as seen in Figure 17.Figure 18
Initial dispersion of ellipsoids after the breakup event within 100 seconds.The simulation is made with different ballistic coefficients using the empirical density and wind models. Figure19 shows the plot of flight time versus ballistic coefficient. In Figure 19, for the variance of the ballistic coefficient in terms of altitude, an empirical density profile obtained from the MSISE-00 density model [22] is used. In addition, in order to include the effects of wind near the ground, a wind model is developed by using a piecewise cubic spline method through a curve fitting of real wind data. Figure 20 represents the wind profiles used in the simulation in the eastern and northern direction as a function of different altitudes. The wind model used in this work is based on HWM-07 (Horizon Wind Model 2007) [7, 24]. Figure 21 illustrates the effects of the wind strength onto the time variation of the dispersion trajectory. As can be seen, if the strength of the wind is stronger, then it takes a longer time until the debris arrives the ground and impact it.Figure 19
Plot of the empirical ballistic coefficient as a function of flight time [7].Figure 20
Simulated eastern and northern wind profiles for different altitudes.
(a)
(b)Figure 21
Effects of the wind to the time evolution of the decay of the debris altitude.


## 4.3. Performance Comparison for Estimation of Debris Dispersion

In this section, the performance of the proposed covariance methods for estimating the debris dispersion due to the reentry breakup is investigated in terms of the time evolution of the positional probabilistic ellipsoid. The covariance propagation methods include the unscented transformation-based statistical linear regression and the Gauss-Hermite cubature-based numerical integration methods described in Section3.For the simulation test, it is also assumed that the breakup altitude is 78 km and the initial breakup position is located atxt0=0078,000m709.9m/s0m/s−124.0m/s. The initial breakup dispersion is modeled as a Gaussian distribution with the initial covariance ellipsoid as given in (34); the process noise covariance matrix representing modeling uncertainties is also given by Qt=106diag00010−1610−1610−16 in (35). Now, based on the initial debris dispersion covariance matrix Pt0 and the process noise matrix Qt, a set of sigma points representing a set of breakup components can be generated by using (14) for the unscented transformation (UT) method and using (29) for the Gauss-Hermite cubature (GHC) method. After the generation of the initial sigma points for the initial debris dispersion, the covariance is propagated by using (16), (17), and (18) for the UT method and (30), (31), and (32) for the GHC approach. The total simulation time used was 32 minutes.Figure22 depicts the estimated nominal trajectory from the breakup event to the ground impact instant, and also the footprint of the ground impact location. After the breakup event of the space debris during reentry, it is assumed that the simulation starts at t=0 at the event, and then the nominal trajectory in terms of the mean value of all the debris fragments and the 3-dimensional positional ellipsoid are propagated until the debris fragments reach the ground. Based on the result in Figure 22, it is seen that the time evolution of the nominal estimates of the debris dispersion have different final locations but the overall trajectories are similar to each other. The nominal trajectory generated by the Gauss-Hermite cubature method has a little bit faster fall into the ground. The locational difference between the final footprint impact locations of two different methods was about 2 km in the x direction, and 10 m in the y direction, that is, δxtfinal=1.8980km0.00860km0km. It is seen from the nominal trajectory plot that the uncertainties of the motion of the debris dispersion are captured more precisely by utilizing the Gauss-Hermite cubature technique compared to the unscented transformation method, and this is because the GHC method adopts a much larger set of sigma lattice points at the initial debris dispersion which leads to a more precise estimation of the debris dispersion. In the simulation, it is also seen that the reentry motion is a nonlinear motion right after the breakup event, but from the altitude near the ground the motion has relatively small nonlinearities due to a vertical falling motion.Figure 22
Performance comparison between unscented transformation and Gauss-Hermite cubature methods in nominal trajectory prediction and ground impact footprint.Figures23 and 24 depict the profiles of the positional probability ellipsoids corresponding to a confidence interval of 99.99% for a sequence of time instants generated from the proposed covariance propagation methods, statistical linear regression, and the Gauss-Hermite moment-matching methods, respectively. In general, since the Gauss-Hermite cubature approach could include more debris breakup distribution points with p=36=729 sigma points, it could generate a more accurate trajectory of the core debris as well as the covariance estimate of the breakup dispersion, while the unscented transformation method generates the initial breakup components with p=2×6+1=13 sigma points. In addition, the performance of the proposed approach is investigated by checking the time evolution of the dispersion ellipsoid shown in Figures 23 and 24 where the time evolution of the positional probability ellipsoids can be seen to almost match between the statistical linear regression in Figure 23 and the Gauss-Hermite cubature-based numerical integration shown in Figure 24. In fact, the close agreement of these techniques can be justified by comparing the positional covariance matrices from both of the methods in each direction. It is seen that the size of the overall debris dispersion is similar to each other, 0.6 km in the x direction and 0.01 km in the y direction. Figure 25 shows the plots of the sigma values for the positional covariance matrices for both of the methods overtime. The difference of the sigma in the x direction is large but the other directions have similar values. The Gauss-Hermite cubature method could provide a little bit more accurate estimates of the nominal trajectory and the positional ellipsoid which could affect the calculation of the risk causality, but it also requires more computational power compared to the unscented transformation method. Both methods could satisfy the computational procedures required for real-time computation in a sequential estimation, while a Monte-Carlo-based approach uses offline computational procedures.Figure 23
Time evolution of debris dispersion ellipsoids and impact footprint on the ground using the statistical linear regression with the unscented transformation method.Figure 24
Time evolution of debris dispersion ellipsoids and impact footprint on the ground using the Gauss-Hermite cubature integration method.Figure 25
Comparison of the magnitude of positional covariances in each direction (x, y, and z).


## 5. Summary and Conclusion

In this paper, for advanced analysis of the time evolution of space debris dispersion due to the breakup during reentry, two effective computational approaches were proposed to estimate the statistical distribution of the debris dispersion. The estimated debris dispersion is represented by using the prediction of positional probability ellipsoids for the visualization of the results. First, the time evolution of the covariance propagation was computed by using a statistical linear regression using the unscented transformation. Second, a novel covariance estimation technique was proposed by utilizing the Gaussian moment-matching method with a Gauss-Hermite cubature-based numerical integration approach. Compared to other covariance propagation methods, the newly proposed Gauss-Hermite cubature-based covariance computation could not only represent more exact initial dispersion, but also precisely calculate the time evolution of the debris dispersion with a large set of sigma points representing debris components. The Gauss-Hermite cubature approach does not require any linearization of the nonlinear translation equations of motion of a reentering debris, which leads to precisely taking into account nonlinearities exposed on the motion of a reentering debris. Furthermore, we also carried out a parametric study in order to analyze the effects on the accuracy of the covariance propagation due to modeling uncertainties. For the detailed parametric analysis, four types of density models and drag coefficients were used in the computation of the lifetime prediction and the covariance computation. In the simulation studies, it was shown that the newly proposed statistical linear regression method and the Gauss-Hermite computational approach are in close agreement to each other, but the Gauss-Hermite covariance propagation method provides a more precise estimate of the debris dispersion of the breakup fragments.



---

*Source: 6824978-2018-07-02.xml*