This study considers the tracking and impact point prediction (IPP) of a re-entry space launch vehicle with unknown ballistic coefficients and air densities. The tracking of a re-entry vehicle (RV) (Costa 1994, Bar-Shalom et al. 2001, Farina et al. 2002) is an extremely important problem for safety assurance through a correct ground IPP. In the case of a space launch vehicle, the major objective is to transport and place satellites into orbit. Thus it is designed to be launched in accordance with an elaborate predetermined trajectory. Thus far, however, quite many cases of launch failures have occurred. If an abnormal thrust halt or an abrupt departure from a planned path happens, it is necessary to track the vehicle or predict the impact point in real time for safety concerns. As for a scenario of a true target trajectory, we assumed that the vehicle has launched normally with beacon mode tracking, but an abnormal situation occurs 170 seconds after lift-off. The vehicle’s thrust is terminated and no communication is available with ground stations. At that time the beacon mode tracking function is halted and the mode should be adjusted as not to lose the target.

Tracking an RV and its IPP are the subjects to estimate the state variables of the dynamic model of a ballistic target. The performance of an estimation filter depends on the target dynamic model, the measurement model, and the filter structure. The uncertainties normally found during target modeling are ballistic coefficients, mass of target, initial values of state variables, and aerodynamic reference area, etc. From among these, the estimation of the ballistic coefficient is crucial for precise tracking of an RV. These subjects have been studied vigorously over the past decades. Mehra (1971) performed a comparison of an RV’s tracking performance using non-linear filters including the extended Kalman filter (EKF). Siouris et al. (1997) proposed the extended interval Kalman filter for the precise tracking of an RV. Cardillo et al. (1999) introduced a new tracking filter using EKF for re-entry objects with uncertain drag. Also, Minvielle (2005) reviewed the RV tracking chronology and introduced several non-linear filters including EKF. According to previous research, various tracking algorithms using EKF were proposed by linearizing the non-linear system to estimate the ballistic coefficients as well as other kinetic states, such as position, velocity, and acceleration. Typically the constant ballistic coefficient (CBC) filter was introduced to estimate the ballistic coefficient β (Dwivedi et al. 2005), the constant drag parameter (CDP) filter to estimate the drag parameter α, which is defined as the inverse of β (Mehra 1971), and the gamma augmented (GA) filter to estimate the gamma (γ) which is defined as the function of β and air density ρ (Bar-Shalom et al. 2001). Recently, Ghosh & Mukhopadhyay (2011) proposed the gamma bootstrapped (GB) filter which estimates basic kinetic states in the main filter with bootstrapping the separately estimated γ and its derivative from the sub filter. If the target modeling is inaccurate in a single model filter, especially for maneuvering target tracking, an estimation error may increase or diverge. To overcome this disadvantage of a single model filter, Blom & Bar-Shalom (1988) introduced the interacting multiple model (IMM) filter algorithm. The IMM filter is known for superior performance in tracking maneuvering targets and has high applicability in a variety of fields. For tracking a ballistic target, Farina et al. (2002) proposed a 4-model IMM filter which consists of a 6-dimensional Kalman filter (KF), a 9-dimensional KF, and 2 different CBC model KFs, and Farrell (2008) proposed a 3-model IMM filter which had a constant axial force model, a ballistic acceleration model, and standard auto-correlated acceleration Singer model. Particularly, Yuan et al. (2010) made use of an IMM filter to estimate the thrusting state for the purpose of IPP.

In this paper, we propose the ballistic coefficient β model based IMM-EKF (β-IMM-EKF) using alternately chosen β models according to the range of real β that the target may have. The β-IMM-EKF is compared and evaluated with the two other ballistic coefficient model based filters, these being the GA and GB filter, since both the GA or GB filter has generally better performance than a CBC or CDP filter. The parameter γ, which is the ratio of air density and ballistic coefficient for the target model, is used as an estimated state variable in the GA and GB filter using EKF. In β-IMM-EKF, γ can be calculated from the estimated β and an air density model to compare its tracking performance with the GA or GB filter. The performance of these filters to track and predict the impact point in the re-entry phase are compared and assessed.

The motion model of an RV is assumed to be a point mass which is applied to the rotating ellipsoidal earth. The difference between the actual and dynamic model is modeled as process noise. The reference frames used for an RV’s dynamic model are the Earth-centered-inertial coordinate system (ECI-CS), the Earth-centered-earth-fixed coordinate system (ECEF-CS) and the East-north-up coordinate system (ENU-CS) (Regan 1984, Siouris et al. 1997, Li & Jilkov 2001).

In Fig. 1, ECI-CS, Ox_Iy_Iz_I, is a fixed inertial frame as its origin is the center of the earth. It is a right-handed system with the origin O. Here the axis Ox_I points to the direction of the line where the longitude crossing the launch point of the vehicle and the equator meet. The axis Oz_I points to the direction of the North Pole. ECEF-CS, Ox_Fy_Fz_F, which has a rotational motion with the earth as its origin, is on the same position with ECI-CS. The radar ENU-CS, Ox_Ny_Nz_N, which has the origin O_N, is fixed on the surface of the earth where the tracking radar is placed. The up direction of ENU-CS, O_Nz_N, is normal to the earth’s surface whereas O_Nx_N axis points east and O_Ny_N axis points north.

Dynamic models (Li & Jilkov 2001, Ghosh & Mukhopadhyay 2011) are used to estimate the tracking of an RV. The issues of the ballistic and re-entry phase can be analyzed with respect to gravity and drag only. We derive the acceleration model of an RV.

2.2.1 Gravity

The acceleration of RV by gravitational force (Regan 1984, Bar-Shalom et al. 2001) in the ECI-CS can be expressed as:

In Eq. (1), μ is a universal gravitational constant as μ=3.986006 × 10¹⁴ m³/sec² and p is the position vector of an RV. Eq. (1) can be expressed in the ENU-CS as:

whereas

is the angular velocity vector of the earth’s rotation,

is the Coriolis force and

is the centrifugal force. From Eqs. (1) and (2), the gravitational acceleration in the ENU-CS is

whereas Ф is the coordinate transformation matrix and r' is the distance between the radar and the center of the earth.

2.2.2 Drag

The acceleration by drag force is as:

whereas the ballistic coefficient (Regan 1984, Dodin et el. 2005) is defined as:

Note that m denotes the mass, C_D denotes the aerodynamic drag coefficient and S denotes the cross-sectional reference area of the RV. β is basically an unknown parameter to be estimated during the tracking of an unknown RV, but if all information of an RV is known, it is quite possible to directly calculate β. In Eq. (4) the air density ρ can be modeled as a function of altitude h as follows:

whereas the parameters ρ₀ and K of the target are given as in Table 1, based on the 1976 US standard atmosphere model

(United States Committee on Extension to the Standard Atmosphere 1976) . Eq. (6) is used to generate the true target trajectory.

2.2.3 The Acceleration Model of an RV

From Eqs. (1) through (4), the acceleration model of RV in re-entry phase is as:

In this section, we introduce several ballistic coefficient model-based estimation algorithms which are evaluated as distinguished tracking performances in previous papers. These are the GA filter and the GB filter. In addition, we designed the β-IMM-EKF to compare the tracking performance with that of the GA and GB filters.

A measurement model is designed and used commonly in filter algorithms. Measurements which are modeled in the radar ENU-CS are the relative range, elevation, and azimuth between the target and the tracking radar. The measurement equation is as:

whereas:

In Eq. (8) the measurement noise vector ν_k and its covariance R_k are white Gaussian sequences.

The GA filter is a type of ballistic coefficient based estimation algorithm, but it regards the variation of air density additionally which significantly affects the nonlinearity of the system. Thus the GA filter estimates γ, which is composed of the ballistic coefficient β and the air density ρ, while a ballistic coefficient augmented filter (Mehra 1971, Dwivedi et al. 2005) does only β. The basic idea of the GA filter is that the new state γ is simply augmented in the state equation as described in Fig. 2. The other procedures are the same as EKF (Grewal & Andrews 2001).

3.2.1 Gamma Model

The ballistic coefficient β is a function of the aerodynamic drag coefficient as shown in Eq. (5). Also, the air density model ρ in Eq. (6) may vary according to the altitude of the target. Thus β and ρ can be estimated and calculated using intricate numerical formulae (Regan 1984). A parameter γ can thus be defined as:

Here, ρ is modeled by the relation as in Eq. (6) and β is approximated by a linear function model of the altitude h (Bar-Shalom et al. 2001) as follows:

Through combining Eqs. (6, 10, 11) and taking the first derivative of γ, we obtain:

3.2.2 State Equations of the GA Filter

The state equation of the GA filter is as:

whereas,

Furthermore,

where

was obtained from Eq. (12).

The GB filter by acceleration and jerk models was proposed in (Ghosh & Mukhopadhyay 2011). The GB filter handles ballistic coefficient β through gamma estimation similarly to the GA filter, but it uses a particular method of bootstrapping. To implement such bootstrapping, the GB filter is divided into two parts, being a main filter and sub-filter as described in Fig. 3.

The main filter is used to estimate the kinematic information including target position, velocity and acceleration while the sub-filter estimates γ and

The estimated results from the separated sub-filter are then bootstrapped into the main filter. Ghosh & Mukhopadhyay (2011) introduced two models of acceleration and jerk for the GB filter and presented that the GB filter had better tracking performance than the GA filter. However, the jerk model may adversely affect the estimates of other

parameters, such as position, velocity, acceleration and jerk, if the rate of change of jerk is not estimated accurately. Thus, we adopt only the acceleration model filter to compare its performance with the other filters.

3.3.1 GB Main Filter State Equations

The GB main filter estimates the target position, velocity, and acceleration by modeling and integrating jerk terms. The gravity, drag and Coriolis force are taken into account in the model states. The state equation of the main filter is as follows:

Hereas:

Furthermore:

Here,

which are estimated in the sub-filter are bootstrapped into the main filter and λ is the latitude of the tracking radar’s

position. We used the EKF for non-linear estimation. The linearized state transition matrix, Ф_k , is required for error covariance propagation. This is the Jacobian matrix of Eq. (23) with respect to X_GB as follows:

whereas the updated value of

at time k is used.

3.3.2 GB Sub Filter

3.3.2.1 Measurement Equations of Sub Filter

At the prediction step, the states γ ,

are predicted in the sub filter and are expressed as

At the update step, the states in the sub-filter are updated by the pseudo measurements γ_c , which are computed using the magnitude of the drag acceleration measurements D as calculated in the main filter. The updated estimate

is used at the prediction step in the main filter. The following equations describe the pseudo measurements γ_c and the measured drag acceleration D. From Eq. (4) we can derive the following:

whereas:

3.3.2.2 State Equations of Sub Filter

Assuming that

is a constant with additive white Gaussian noise, we can adopt the KF algorithm. Thus the GB sub-filter can be derived through a constant velocity model as a quadratic form:

The state variables are updated using the pseudo measurements as in Eq. (30).

whereas η_γ(k) is the measurement noise assumed to be white Gaussian. The measurement noise covariance R_GB(k) of the sub-filter is computed using the Jacobian matrix H_γ.

Whereas:

The IMM filter (Blom & Bar-Shalom 1988), which combines probabilistically the estimated variables from multiple dynamic models is well known as an algorithm that has similar computation load with the generalized Pseudo Bayesian approach of the first order (GPB1) (Bar-Shalom et al. 2001) and has estimation accuracy with the generalized Pseudo Bayesian approach of the second order (GPB2) (Bar-Shalom et al. 2001). In this section, we introduce the IMM-EKF algorithm and derive the state equations of β-IMM-EKF. We also designed various β according to hypotheses or modes to verify the tracking performance of this filter.

3.4.1 The IMM-EKF Algorithm

A distinct feature of the IMM filter algorithm is the combination process where the estimated mean and covariance are produced through the computation of mode probabilities. Mode probability means the probability where each filter dynamic model matches the target dynamic model. The process noises

which mean the uncertainties of a dynamic model, are white Gaussian noises with a zero mean and a covariance of

The IMM filter algorithm is composed of the steps of interaction, prediction, update, and combination which are derived under Bayes’ rule. The flow of probability density of each step is represented in the following Fig. 4.

3.4.1.1 Interaction

At the interaction step as in Fig. 4 the estimated mean

and error covariance are obtained by the conditional probability density function (CPDF),

which is computed using the transition probabilities of each hypothesis as in the following:

Here, the mode probability

and the state transition probability π_ij from mode j to mode i are defined as follows:

3.4.1.2 Prediction

At the prediction step as in Fig. 4, the predicted state and error covariance of

is calculated as the manner of EKF:

whereas:

Here

is the Jacobian matrix that uses partial derivatives of the non-linear function in Eq. (57) with respect to X.

3.4.1.3 Update

At the update step in Fig. 4, the estimate and error covariance of

are calculated as EKF:

whereas:

Mode probability

which means the probability of the hypothesis

is true, can be calculated as follows:

Here the likelihood function,

is calculated using measuring residual

and innovation covariance

as follows through the condition of linear Gaussian.

Here, n is chosen as 3 for the dimension of the system.

3.4.1.4 Combination

At the combination step, one combined estimate and error covariance are computed with the estimates

the error covariance

and the each mode probabilities

which are obtained from the multiple modes according to the Gaussian mixture model as follows:

3.4.2 The Ballistic Coefficient Model based IMM-EKF

In this paper we propose the β-IMM-EKF, which tracks a non-linear system such as an RV, with modeling of the Markov chain for N-hypotheses. This algorithm is to be compared with GA and GB filter algorithms for tracking performance. The dynamic model for β-IMM-EKF is designed as follows:

whereas:

Hereas:

Here, i means the i-th hypothesis of the i-th mode target dynamic model. Each model has different ballistic coefficients in different modes. Thus the β term in Eqs. (59) through (61) for β-IMM-EKF should be seen instead of γ as in Eq. (10), which is used in either the GA or GB filter.

3.4.3 Modeling of β

The tracking performance of β-IMM-EKF depends on how well the multiple β models are selected, and this might be a demerit of the IMM-EKF algorithm. In other words, without knowledge of the area of the model parameters to be estimated, it may produce poor performance when applying the IMM filter algorithm for tracking. In our application, however, it can be assumed that the vicinity of true β values is known in advance as the true RV model was selected from a known space launch vehicle (SLV). Generally, the ballistic coefficients of the SLV are the factors to be obtained from the design and development phase of the vehicle. This means that the major ballistic coefficients are in accordance with the different flight stages for an SLV and are known prior to launch. Even though the exterior shape of the vehicle changes owing to the abnormal conditions during flight, β values can be expected to have a neighboring range in their predetermined values. Thus the β-IMM-EKF algorithm can be used in this application.

Fig. 5 shows the true β and γ of the nominal trajectory. In Fig. 5a the range of true β value is about 750-1,320 kg/m², which were acquired by computer simulations during

the flight range from 300 seconds to 482 seconds. Fig. 5b shows the gamma variation through the computation of (10) with the corresponding β. In this paper, the performance comparison of β-IMM-EKF with the GA or GB filter is executed with respect to the parameter γ. The gamma value can be easily calculated with the estimated β using Eqs. (6) and (10). We can observe that γ has a significant rise from approximately 400 seconds; likewise, β also shows a rapid change at that time. This occurs due to the influence of a rapid increase in air density according to the close approach of an RV to the surface of the earth.

In Table 2, we selected three hypotheses or modes since the variation of the β value in Fig. 5a is not very wide from the minimum to maximum values, and since the computation time can be shortened by small numbers of β models. Thus we designed three ballistic coefficient models with the values of 800, 1,200, and 1,600 kg/m² respectively for the three hypotheses based on the range of true β as in Fig. 5a.

Fig. 6 shows the change of mode probabilities at each mode after tracking simulation for the true RV model. Since observability is not built, three modes have the same averaged values until approximately 350 seconds of flight time. Considering the entire flight during the re-entry phase, the mode probabilities begin to actively change after about 350 seconds when observability is formed. Thus, after about 350 seconds, each mode has its own variations of probability according to the flight conditions which are related with Eq. (5).

From approximately 380 to 420 seconds, the probability of Mode 2 rises to about 0.6 while that of Mode 1 and 3 have

0.25 and 0.15 each. From 420 to 440 seconds, the probability of Mode 1 rises to about 0.7. After about 450 seconds, the probability of Mode 3 rises to about 0.5 until the RV impacts the surface of the earth.

For the true trajectory of an RV, the specifications of a small low earth orbit satellite launch vehicle, the Korea SLV-I (KSLV-I), was used in the computer simulation. The true target trajectory assumes that the vehicle is safely launched at first. The predetermined thrust duration is 230 seconds but the thrust is assumed to be cut off due to an abnormal condition at Lift-off +170 seconds. Then the RV enters the coast phase. For simple analysis, the conditions of body tumbling or splitting of the RV are not considered. The launch point is set at longitude 126 degrees east and latitude 34 degrees north. The launch direction is 170 degrees from the North Pole. The results are 482 seconds of total flight time, 427 km of down range, and 113 km at the highest altitude for the nominal trajectory. The radar station for tracking is located at longitude 127 degrees east, latitude 19 degrees north. The skin mode tracking begins at Lift-off 190 seconds when the relative distance between target and radar is about 1,594 km. The tracking filter receives measurements from the sensor at a 16 Hz cycle in the form of range, elevation, and azimuth (r,θ,ψ) in polar coordinates. The measurement noise is white Gaussian with standard deviation (25 m, 0.05 deg, 0.05 deg). Figs. 7a-c represents the true target trajectory with respect to the East, North, and Up axes, respectively. The trajectory was generated in the

radar ENU-CS after proper coordinate transformations were performed from the ENU-CS referenced at the launch point. The entire trajectory is represented by the solid black line, and the tracking interval is expressed by the scattered red circles on the trajectory.

Fig. 8 shows the estimation root mean square error (RMSE) in terms of position, velocity and γ for each algorithm. In Figs. 8a and b, the position and velocity of RMSE, both the β-IMM-EKF and GA filter generally has better estimation performance than the GB filter during the tracking time from 190 seconds to 482 seconds. The GB filter cannot reduce errors less than about 190m of position and 30 m/s of velocity, whereas the GA filter maintained an error trend similar to that of the β-IMM-EKF until about 420 seconds. In the enlarged section of Figs. 8a and b from 400 to 500 seconds, the GA filter shows a rapid wave of position error up to 240 m at 420 seconds and twice oscillations of its velocity error occurred from about 0 to 18 m/s. This result is because the GA filter cannot follow the fast change of γ related to the ballistic coefficient and air density according to Fig. 5 and Eq. (6). In the case of the GB filter, the estimation of jerk terms as in Eq. (23) has a vital influence on overall system performance, but an accurate estimation of jerk terms, as derivatives of acceleration terms, seems not to be completely achieved. This means that the GB filter shows limitations in comparison of that of the β-IMM-EKF or GA filter. Thus, the β-IMM-EKF achieved a stable and more accurate result than both the GA and GB filter. In Fig. 8c of gamma RMSE, three models of GA, GB, and β-IMM-EKF show good performance until 400 seconds. This result comes from very low air density by Eq. (6). However, in the enlarged section of Fig. 8c, the GB filter shows the largest error in the results after about 440 seconds amongst the three filters since it does not obtain full estimation performance due to the jerk terms. Between the GA filter and β-IMM-EKF, the GA filter shows a smaller number of errors at some parts than the β-IMM-EKF from about 430 to 475 seconds. But the error trend of the GA filter is unstable with oscillation, while that of the β-IMM-EKF is stable. Overall, the β-IMM-EKF has more improved and stable performance than that of the other two filters.

Table 3 shows the processing time per scan for each algorithm based on our simulation environments. As for the run time speed per one scan, the GA filter proved to be the fastest. This result can be expected regarding the dimensions of each filter. As described, the state equations in Sections 3.2-3.4, the dimension of the GA filter is 7 whereas that of the GB filter is 9 and that of the β-IMM-EKF is 18. Thus the complexity of the β-IMM-EKF is the highest. The GB filter actually requires considerable computation of the Jacobian matrix when its sub-filter also works to calculate pseudo measurements. The complexity of the β-IMM-EKF depends on the number of models. Since we propose three models of the β-IMM-EKF in this paper, it has three times the complexity of each 6-dimension EKF.

The run time results represent that it seems to be approximately proportional to the dimensions. That is, the GA filter is somewhat faster than that of the GB filter, but rather had similar results since the difference is only 2 dimensions between the two. The run time results of the β-IMM-EKF showed almost 3 times more than that of the GA filter as its dimension is 3 times larger. Judging from the results, β-IMM-EKF is the slowest but the differences are insignificant, being 0.021 to 0.327 msec among the 3 filters. Thus all 3 filters are applicable for real time tracking.

The comparisons of IPP are conducted at 5 points. These are the points at 350 and 400 seconds where drag has a small effect due to light air density; and the points at 425 and 450 seconds where drag has a significant effect at less than approximately 35 km altitude. Only the prediction step for each filter is performed without measurements. Lastly, the impact point is shown through full tracking from 190 to 482 seconds.

Fig. 9 shows the results of the predicted impact points obtained from each starting time by performing 100 Monte Carlo runs for each algorithm. In Figs. 9a and b, both the

GA and GB filter show an error radius up to about 15 km around the true impact point while the β-IMM-EKF shows less than a few kilometers. This means that observability for the GA or GB filter is not attained sufficiently at about 100 km altitude of the RV prior to prediction beginning, while β-IMM-EKF has better results than the others owing to the selection of proper modeling. In Figs. 9c, and d, the GA and GB filter was closer to the true IPP. This is due to the rising of observability to estimate how γ is formed. Here the tracking performance starts to influence the accuracy of the IPP. In Fig. 9e the results of the final IPPs are shown through the full tracking of each filter. Though all the filters showed good results of IPP around the true at less than 1km, β-IMM-EKF still showed better results than either the GA or GB filter. As a result in Figs. 9a-e, β-IMM-EKF shows noticeably better prediction performance amongst the 3 filters.

Table 4 shows the values of circular error probable (CEP) based on the results of the IPPs for the three algorithms. The CEPs are computed by mean square error between the true impact point (TIP) and the predicted impact point (PIP) of the target. Each value represents the maximum radius from TIP as the center of the circle where the 50% of PIPs are located. In all cases, we can find that the β-IMM-EKF showed the best results against the others.

A comparative assessment of tracking performance and IPP of an RV has been performed for the GA filter, GB filter, and the proposed β-IMM-EKF, which are ballistic coefficient model based algorithms. To verify practical performance of tracking, the RV model has been applied with actual specifications of a space launch vehicle. With this RV model as a true target, the tracking and IPPs through Monte Carlo simulations were performed for those three algorithms.

Concerning the RMSE of position, velocity and gamma, the β-IMM-EKF produced the best tracking performance amongst the 3 filters. Comparing the results between the GA filter and β-IMM-EKF, both showed a very similar

equivalent performance during certain sections. However, some unstable oscillations appeared in the GA filter while no oscillations in occurred in the β-IMM-EKF when the RV approached the surface of the earth.

Concerning the IPP, β-IMM-EKF still showed the best results among the 3 filters. The IPP was performed without measurement updates for five cases. At an altitude of about 100 km, the GA and GB filter both have larger prediction errors than the β-IMM-EKF since observability was not attained. When the influence of air density increases and observability is acquired from approximately 35 km altitude, the predicted impact points of the GA or GB filter started to point closer to the true impact point. But in any case, the β-IMM-EKF had better performance than the other two. Thus IPP after β-IMM-EKF tracking is more confident and less affected by the altitude of an RV.

Concerning the algorithm execution time per scan, the GA filter was the fastest whereas the β-IMM-EKF was the slowest. This result was expected due to the smaller dimension of the GA filter than the others. However, this does not mean that either the GB filter or β-IMM-EKF is inapplicable to a real environment. As identified in Table 3, the maximum differences amongst the three filters are less than 1 msec. Thus any of the 3 filters in real time are able to be used in real tracking. Moreover, the run time per scan also depends on the capability of the computer.

In conclusion, the proposed β-IMM-EKF showed better performance than the other two filters in comparison of RMSE and IPP. In regards to the run time per scan, such is minor since the difference between the filters is quite small. As discussed in Subsection 3.4.3, β-IMM-EKF can offer stable and good tracking performance through proper modeling of β. Thus, if the pre-determined range of the ballistic coefficients is known, it can be expected that it is able to design a distinguished good estimator with the β-IMM-EKF for applications of precise tracking of an RV.