Aerodynamic Statistics-Based Trajectory Estimation of Hypersonic Maneuvering Target

Trajectory tracking and estimation of hypersonic glide vehicles (HGVs) is a very challenging issue in the defense systems. Insufficient knowledge about the HGV and inaccurate dynamic models for the accelerating HGV are the main challenges in this regard. In the present study, an integrated nonlinear Markov acceleration model is established to formulate the nonlinear dynamic characteristics of HGVs. Since the aerodynamic accelerations of the HGV are dominant and the corresponding aerodynamic coefficients are unknown, a statistics-based aerodynamic model is proposed. The proposed aerodynamic model is capable of providing primary information of the aerodynamic characteristics even without knowing the configuration of the HGV. Then, considering the maneuver mode of the vehicle, the iterative extended Kalman filter (IEKF) is applied to track the trajectory of the HGV by using the proposed model. Obtained results from the numerical simulation for the equilibrium glide mode and skip maneuver mode indicate that the proposed model can effectively improve the velocity estimation accuracy by about 40%-50% and acceleration estimation accuracy by about 20%-50% in the given examples.


I. INTRODUCTION
Studies show that flexible maneuvers can be carried out by hypersonic glide vehicles (HGVs) to increase their survival rate in real operations [1]. This capability originates from the high maneuverability of HGVs. As a result, the trajectory tracking and estimation of HGVs are more challenging for the defense systems, when the comparison is made with that of ballistic vehicles. Generally, trajectory tracking and estimation are fulfilled through a filter, which combines measured data received from the radar and that of the reference model [2]. Accordingly, the reference model has a significant impact on the estimation accuracy of the trajectory [3], [4]. This is especially more pronounced for state variables, which cannot be measured directly. For ballistic targets, the ballistic coefficient or accelerations are almost constant or slowly change so that these coefficients can be modelled as a constant or a first-order Markov process [5].
The associate editor coordinating the review of this manuscript and approving it for publication was Yuan Zhuang .
Consequently, the target state of the ballistic target, including the velocity, position and acceleration can be accurately estimated. However, maneuvers of the HGV rely on the aerodynamic forces within the atmosphere so that the drag and lift accelerations are prominent and may rapidly change [6]. Since the acceleration data of HGVs can be hardly detected directly by the radar of the defense system, establishing an appropriate reference model of the filter is essential to accurately estimate the instantaneous state of the HGV.
Almost all maneuvering target trajectory estimation methods are model-based. An excellent tracking model should describe the true process of target motion and the moving characteristics. Reviewing the literature indicates that filter models for trajectory tracking can be roughly divided into two categories, including the kinematic models and the dynamic models. In this regard, different kinematic models such as Singer model [7], second-order time-dependent model [8], constant turning model [9], semi-Markov model [10], jerk model [11] and adaptive current statistical model [12], have been proposed so far. Studies show that although these models have reasonable real-time performance, they have some inherent drawbacks such as the lack of characterization of innate dynamics and inability to consider the nonlinearity of the dynamics. In order to develop new dynamic models, the concept of maneuvering coefficients was introduced to describe the motion of a maneuvering missile [13]. Then the state augmentation method was adopted to estimate the ballistic and maneuvering coefficients [14]. Fan et al. [15] introduced the angle of attack and the bank angle as two control variables in the Gaussian models and applied the cubature Kalman filter to estimate the target state with fixed control variables. In the foregoing studies, key parameters formulating the accelerations of the HGV are modelled as constants with zero acceleration derivatives or linear Gauss-Markov processes. However, the dynamic pressure, which has a great impact on the accelerations of the HGV, has not been fully considered yet. Moreover, many studies were performed to track tactical ballistic missiles adopting spiraling motion during the reentering phase. Dubois-Matra and Bishop [16] implemented time derivatives in a Kalman filter system and developed a flight model to track a spiraling maneuvering reentry vehicle accompanied by decoys. Tapiero and Bishop [17] established a model considering stochastic uncertainty parameters from ballistic geometry for tracking spiraling reentry vehicles. More recently, Hough [18] proposed nonlinear Markov acceleration models with variable process noise amplitudes. He established the dynamics of the accelerations and decomposed affecting factors into several components, which can describe the dynamics of HGVs. However, there are still some unsolved challenges to be faced. First, due to insufficient knowledge about the noncooperative HGV, configuration parameters and aerodynamic coefficients are unknown to the defense system. Consequently, the model mismatch, which can significantly deteriorate the tracking accuracy, may occur. Second, the accelerations of the HGV cannot be measured directly through the radar. Consequently, the nonlinear dynamics coupled with an inaccurate estimation of the accelerations may magnify tracking errors of the vehicle state.
Inspired by the aforementioned analysis, in the present study, it is intended to propose a statistics-based nonlinear dynamic model to solve the trajectory tracking problem of HGVs. To this end, the dynamic model of the acceleration components related to the position, velocity, atmospheric parameters and statistics of expected manoeuvres, is established. Since the configuration of the target HGV is unknown, a statistical analysis is initially performed for a family of HGVs to acquire the initial information of the primary aerodynamic characteristics of HGVs. Finally, the IEKF algorithm is conducted to accomplish the trajectory tracking.

II. DYNAMIC MODELS OF HGVs
In order to establish the dynamic models of the HGV, the earth is assumed as a homogeneous sphere. Considering the rotation of the earth, the three-dimensional point-mass dynamics of the HGV are constructed in a semi-speed coordinate frame (VTC) [19]. It should be indicated that the sideslip angle of the HGV is neglected in all calculations.
where r, θ, φ, v, γ , ψ are radial distance from the centre of the Earth to the vehicle, longitude, latitude, velocity of the HGV, flight path angle and the velocity heading angle, respectively. Moreover, σ and g are the bank angle and the gravity acceleration, respectively. C γ and C ψ are additional terms caused by the rotation of the earth [20]. The magnitudes of the drag and lift accelerations (a D and a L ) are related to aerodynamic coefficients, reference area, mass and dynamic pressure [19]. These parameters can be mathematically expressed as the following: where S and m denote the reference area and mass of the vehicle, respectively. Furthermore, C D and C L are drag and lift coefficients of HGV, respectively. q = ρv 2 /2 is the dynamic pressure and ρ is the atmospheric density, which can be expressed as the exponential function of height, in the form below [19]: where ρ 0 =1.225 kg/m 3 is the atmospheric density at the sea level and h denotes the altitude of the target distance from the sea level. Meanwhile, h = r − R e , where R e is the radius of the earth. Finally, hs = 6700 m is a constant.

III. GLIDING DYNAMICS FOR TRACKING
In order to track the trajectory of the HGV, the dynamic model is a base for the state estimation. In this section, the dynamic model of acceleration components for the HGV is derived.

A. DYNAMICS OF HGVs
Radar reference coordinate system (RRCS): In this coordinate system, the origin coincides with the radar station position, while unit vectors of X-and Y-axes are located in the tangent plane of the earth's reference sphere, along East and North directions, respectively. Meanwhile, the unit vector of the Z-axis is perpendicular to the local horizontal plane. Accordingly, the RRCS is also called the local East-North-Up (ENU) geographic coordinate system. VOLUME 8, 2020 During the gliding phase, the HGV motion relative to the radar reference coordinate can be mathematically expressed as: where r and v represent the center distance vector and the velocity vector of HGV, respectively. Moreover, a denotes the acceleration of HGV originating from aerodynamic forces; ω e is the angular velocity of the earth; −ω e ×(ω e ×r) denotes centrifugal acceleration; −2ω e ×ṙ denotes Coriolis acceleration; The gravity acceleration of gravity can be expressed as: where µ is the gravitational parameter of the earth.   [21]. Where u v , u t and u c are unit vectors of each axis in the VTC coordinate. Furthermore, a d , a t and a c are acceleration components of the drag, turn and climb, respectively. These components can be expressed as follows: or, where γ V is the bank angle. The aerodynamic acceleration vector can be decomposed in the VTC coordinate system [22]. This can be mathematically expressed as the following: where κ d , κ t and κ c denote the drag, lateral and longitudinal aerodynamic parameters, respectively. The total lift acceleration and aerodynamic parameters are defined as a L = a 2 t + a 2 c and κ L = κ 2 t + κ 2 c , respectively. Corresponding to the components of the lift acceleration, the angle of attack is also decomposed into the lateral and longitudinal components.
Lift coefficient C L (α,Ma) is a function of the angle of attack and Mach number, where the Mach number can be calculated from the following expression: It is worth noting that the sound velocity cs is a highly nonlinear function [19]. Drag coefficient can be described in the form of the drag pole curve [23]: where C 0 D is the zero-lift drag coefficient when no lift is generated by the HGV and K is a constant.
The critical lift coefficient C * L is the lift coefficient with a maximum lift-to-drag ratio.

B. ACCELERATION MODEL BASED ON NONLINEAR MARKOV PROCESSES
Maneuver accelerations of HGVs are affected by a few trajectory parameters and guidance commands, including the height, velocity, atmospheric parameters, aerodynamic characteristics of the HGV, angle of attack and the bank angle. Dubois-Matra [16] and Tapiero [17] investigated the time derivative of acceleration components to have a deep insight into the dynamics of acceleration.

1) COMPONENTS OF THE LIFT ACCELERATION MODEL
In order to simplify the calculations without loss of generality, parameters a t and a c are replaced with a j (a j {j=t, c}). Based on Eq. (8), the derivative of a j can be expressed in the form below:ȧ Assuming that lateral and longitudinal components of the angle of attack are independent of each other, the first-order Markov processes can be established as: where τ j is the time constant and α j,ctrl is a term that controls the angle of attack.
Based on Eqs. (3), (8), (9) and (13), each term in the right hand side of Eq. (12) can be rewritten as the following: where a j,ctrl is the acceleration guidance command. The velocity derivative can be calculated in the form below: Meanwhile, the derivative of the flight height can be calculated from the following expression: The lift acceleration models yield the following result: Eq. (19) shows that the lateral and longitudinal lift accelerations a t and a c are established by nonlinear first-order Markov processes. The term S 1 is determined by the vehicle's state and leads to the nonlinearity in the dynamic model of the lift accelerations. If S 1 = 0, the dynamic model of the lift accelerations follows the linear Markov process. Moreover, the terms S 2 and S 3 denote the maneuvering frequency corresponding to aerodynamic parameters and the bias of the lift acceleration induced by the change of the state variables, respectively. It should be indicated that in the term S 3 the aerodynamic coefficients should be known. Furthermore, the term S 4 stands for the variation of the lift acceleration caused by guidance commands. Considering that the guidance commands of the acceleration a j,ctrl are generally not obtained by the defense industry, it can be assumed as a j,ctrl = 0 [18]. The terms S 3 and S 4 change the mean value of the distribution of the lift acceleration.
In the acceleration components caused by the lifting force, only the time constant τ j and the partial derivative of aerodynamic parameters to Mach number ∂κ j ∂Ma are not directly correlated to the state variables of the HGV. It should be indicated that the time constant can be set empirically.
The partial derivative of aerodynamic parameters to Mach number is described in section IV.

2) DRAG ACCELERATION MODEL
According to the description of the drag pole curve in Eqs. (8)-(10), the total differential of the drag acceleration is mathematically expressed as: Moreover, where n is the lift-drag ratio of the HGVs, which can be obtained by distribution statistics. Then: The drag acceleration model for the estimation is described as: where, the autopilot time constant is τ L = (τ t + τ c )/2, and the total lift acceleration guidance command is a Lc = a 2 t,ctrl + a 2 c,ctrl . Eq. (24) shows that the drag acceleration is determined by the nonlinear and inhomogeneous differential equation, which includes the linear and quadratic terms of the drag acceleration (Term Q 1 ), the total lift acceleration (Term Q 2 ), variable coefficient depending on the altitude, velocity and partial derivatives of aerodynamic parameters to Mach number (Term Q 3 ) and the guidance command of the total lift acceleration (Term Q 4 ). In the drag acceleration component, the partial derivative of aerodynamic parameters to Mach number ∂κ d ∂Ma , ∂κ L ∂Ma and the lift drag ratio n of the HGVcan be obtained by statistical methods, which is described in section IV.

IV. AERODYNAMIC CHARACTERISTICS STATISTICS
In Eqs. (19) and (24), the partial derivative of the aerodynamic parameters to the Mach number and the lift-to-drag ratio should be known. However, such parameters are generally unavailable for a noncooperative vehicle. Moreover, VOLUME 8, 2020 considering the typical characteristics of a hypersonic gliding vehicle, the prior information about the configuration can be inferred. In this section, a parameterized model of the hypersonic technology vehicle 2 (HTV-2) is established. Then, the typical configurations are investigated to obtain the aerodynamic characteristics. Finally, the statistical analysis of the aerodynamic characteristics is performed to approximate the unknown partial derivative of the aerodynamic parameters to Mach number and the lift-to-drag ratio in nonlinear Markov maneuver accelerations.

A. PARAMETRIC MODEL OF HTV-2-LIKE HGVs
It should be indicated that although the configuration of the noncooperative HGV cannot be obtained by the defending system, its configuration is subjected to some constraints in order to achieve specific tasks. Therefore, according to the prior knowledge of the HGV, it is possible to parameterize the typical HGV approximately to provide deep insight into aerodynamic characteristics. Fig. 2 shows that HTV-2 liked configuration, which has an appropriate tradeoff in aspects of high lift drag ration and thermal protection is considered as a baseline in the present study.  The configurations of the HGV similar to HTV-2 are mainly determined by the following design parameters: The full length, the height and the maximum width of the aircraft, the radius of the bottom arc and the coordinates of the shape control points in the bottom cross-section [24], [25]. Table 1 presents the contour parameters corresponding to the parametric modeling of the HGV. Moreover, Fig. 3 and Fig. 4 show the bottom cross-sections of the HGV obtained by  choosing different design parameters. It is worth noting that X and Y coincide with the axis of the aircraft body coordinate frame.  Table 2 shows that it is necessary to analyze the general aerodynamic performance for a family of the HGV configurations, which are generated by the optimal Latin hypercube design method [26]. Analogously, the aerodynamic characteristics of the parametric model set of HGVs are analyzed. For each configuration, the aerodynamic coefficients of the HGV are fully calculated under various flight conditions. The angle of attack varies from −10 • to 10 • , Mach number from 3 to 20, and altitude from 30km to 50km. For the space limitation, only the cases of Ma = 10 and α = 10 • are listed as follows.    The configurations in Table 2 show that (1)     linear function of the aerodynamic coefficients concerning the Mach number are close. Therefore, the mean value of the slope and intercept can be used as the approximation of uncertain parameters. Fig. 9 and Fig. 10 show that by studying the lift-to-drag ratio of the typical configurations, the following conclusions are drawn: (1) The maximum lift-to-drag ratio appears at about 10 degrees of the angle of attack. (2) The maximum lift-to-drag ratio of the HTV-2-like HGV is around 4, which is in accordance with the results of the existing literature [27].
(3) In general, the lift-drag ratio of the HGV positively VOLUME 8, 2020 correlates with the angle of attack. (4) The lift-to-drag ratios of HGV are basically independent of the Mach number.

C. AERODYNAMIC PARAMETER STATISTICS
The Kernel density estimation (KDE) technique [28] is utilized to determine the probability density function (PDF) quantitatively for each unknown parameter in the acceleration component. Each PDF is expressed in ρ e (x). The mean value of the parameters, the variance and the root means square error are expressed as µ x , P x and S x , respectively. These statistics are determined by numerical integrations of ρ e (x) for each parameter component:  The five configurations and flight state sets that are composed of different angles of attack and Mach numbers, define ρ e (x) for each unknown parameter components. Assuming that the flight states of the HGV are uniformly distributed, the statistical results of the probability density function at different angles of attack and Mach number are given in Fig. 11 and Fig. 12. It should be indicated that the aerodynamic coefficient and lift-drag ratio are asymmetrically distributed. For the lift coefficient, ρ e (Cl) is distributed between −0.5 and 1, and the maximum value of the PDF appears at 0.09. For the drag coefficient, ρ e (Cd) is distributed between 0 and 1 (even though drag is always negative), and the maximum value appears at 0.12. For the lift-to-drag ratio, ρ e (n) is distributed between −3 and 6. Table 3 gives the statistical results of aerodynamic characteristic parameters. Since the distribution of the drag coefficient is more concentrated near the maximum value, the variance and root mean square error of the drag coefficient are lower than lift coefficient. Moreover, Fig. 13 shows that the statistics of derivatives of aerodynamic coefficients to Mach numbers are also required. It is observed that the derivatives of the lift coefficient and drag coefficient to the Mach number are both negative (near zero), while the distribution of derivatives of the drag coefficients to Mach numbers is more concentrated. In summary, Table 4 presents the uncertain parameters such as derivatives of aerodynamic coefficients to Mach numbers. Moreover, the lift-to-drag ratio in Eqs. (19) and (24) are obtained by statistical methods.
What needs to be emphasized here is that the purpose of this paper is to study a novel model, which can effectively realize the trajectory tracking for the unpowered long-range HGVs. In order to have the ability to complete some specific tasks, the aerodynamic configuration of this kind of vehicles is restricted by many constraints, such as large lift-to-drag ratio, appropriate volume ratio, etc. Therefore, the methods of parametric modelling and aerodynamic statistics proposed in this paper show the applicability in solving the tracking problem of the specific HGVs like HTV-2.

V. ITERATIVE EXTENDED KALMAN FILTER
Since the dynamics of acceleration components are established, the accelerations can be estimated by a filter. The acceleration components are augmented in the state variables. Then, the IEKF method is utilized to estimate the states of the trajectory.

A. STATE VECTOR DIFFERENTIAL EQUATION
The IEKF with nine state variables is used to estimate the 9 × 1 state vector and 9 × 9 covariance.
The elements in the state vector X include the HGV's position components x, y, z and the velocity components v x , v y , v z described in the RRCS. Moreover, the acceleration components a d , a t and a c are resolved in the VTC coordinate frame.
The aerodynamic acceleration vector is conversed to the RRCS coordinate system, which is mathematically expressed as the following: where T RRCS VTC is the coordinate transformation matrix from the VTC coordinate frame to RRCS coordinate frame of the radar station. T RRCS VTC is described in terms of v [21]: where v = ẋ 2 +ẏ 2 +ż 2 and v g = ẋ 2 +ẏ 2 are the 'total velocity' and 'in-plane velocity', respectively. The continuous state equation of the estimation model is described as:Ẋ = f (X) + w (29) where w is a process noise, which is assumed to be the zero-mean Gauss white noise. The nonlinear function f(X) based on equations (19) and (24) can be expressed as: By discretizing the nonlinear function f(X), the following results are obtained: where, t is the sampling interval, which is correlated to the tracking data rate of the radar. Moreover, F(X k ) is the Jacobian matrix of f(X k , t k ), which is relative to X k .

B. RADAR MEASUREMENT MODEL
The states of the HGVs are estimated in a mixed coordinate system. In other words, the position and velocity of the HGVs are described in the ENU rectangular coordinate system, and the measurement equation is established in the radar spherical coordinate system. Considering the gliding trajectory of the HGV, a radarbased measurement model is proposed in the present study. The radar is assumed to measure the relative range and two gimbal angles.
where, R, A and E denote the relative range, the azimuth angle and the elevation angle, respectively. VOLUME 8, 2020 Consequently, the nonlinear measurement equation is as follows: where x n represents the state in Eq (26). Moreover, n R , n A and n E are radar measurement noises, assuming that they are the zero-mean white noise. Partial derivatives of h(x n ) with respect to the state variables specify the 3 × 9 measurement Jacobi matrix: Measurement errors are approximated by the Gaussian random variable v n with 3 × 3 variance matrix, which is described as the following:

C. ITERATIVE EXTENDED KALMAN FILTER UPDATES AND PREDICTIONS
For the conventional EKF method, the Taylor expansion highorder items of h(x n ) are omitted in the process of linearizing the observation equation [18]. During each iteration step, error will be magnified once, which leads to the worse effect of the maneuvering target state estimation or the divergence of filtering.
Moreover, the IEKF requires the algorithm to linearize the observation equation at each observation time and repeat the filtering stage to obtain the optimal estimation. The specific steps are as follows: 1) In the k-time, according to the observed data Z k , the EKF method is used to obtain the state predictionX k/k−1 and the prediction covariance matrix P k/k−1 .
2) Select the relative error between the state estimation and the observation as the variable. Eqs. (37) shows that this error is also known as the residual. Then, set a reasonable threshold according to the tolerable relative positioning error. In the iteration cycle of each step, the iteration operation is stopped when the parameter is less than the threshold value. It should be indicated that the number of iterations is expressed as N . Each initial iteration is determined by a prior estimation of the state and covariance of the previous prediction, i.e.X 1 k/k = X k/k−1 and P 1 k = P k/k−1 . The filtering gain K N k , state updatê X N k/k and filtering error covariance matrix P N k are obtained by using the observation data at k time and the EKF updating method.
3) Let X k/k =X N k/k , P k = P N k ,and then use the EKF method to predict the updates by introducing the observation data at k + 1 time. Then, repeat step 2 until the end of the filtering.
From the above steps, it is observed that the IEKF method improves the correctness of the filter covariance with nonlinear measurements. The residuals are evaluated using a nonlinear measurement function: The IEKF algorithm is an approximate method for the nonlinear system. Its iteration process is equivalent to applying the Gauss-Newton method to the iterative optimization of the least-squares residual problem.

VI. NUMERICAL SIMULATION EXPERIMENTS
In order to verify the proposed method, a simulation is conducted to track the HGV. It is worth noting that the HGV weighs 970 kg and has a reference area of 1.02 m 2 . The trajectories of HGVs are generally divided into two forms, including equilibrium glide mode and skip manoeuvre mode [29]. Considering the realistic environment of the pursuit-evasion confrontation, both flight modes are tested.
It is assumed that the ground-based radar is located in the east longitude 165 • and north latitude 12 • . Moreover, the range of azimuth, the range of elevation and the sampling rate are A ∈ [−180 • , 180 • ), E ∈ [0 • , 20 • ) and 1 Hz, respectively. Furthermore, the standard deviation of the radial distance measurement and the standard deviation of azimuth and elevation measurements are 100 m and 2 × 10 −4 rad, respectively.
For the comparison with the proposed acceleration model, three other acceleration models, including the linear Gauss-Markov processes [30], adaptive Current Statistical (CS) model [31] and the constant acceleration model [29], are utilized in the test case. The statistical average results are obtained through 100 Monte Carlo simulations. It should be indicated that root-mean-square error (RMSE) and average root-mean-square error (ARMSE) are utilized as evaluation metrics: where, M and S denote the total number of Monte-Carlo simulations and the total steps of the simulation, respectively. Moreover, X i k andX i k are the real value and the estimation value in time k t at i th simulation, respectively.    plane approaches zero. The trajectory of the equilibrium gliding HGV is designed in the VTC coordinate frame from the angle of attack and bank as two control variables, which is discussed in section II. The initial altitude, velocity and path angle are 90 km, 6500 m/s and 0 • , respectively. Fig. 14 shows the height and velocity distributions of the HGV during gliding. In this mode, the path angle of the HGV is small and its rate approaches zero. Fig. 15 depicts the geometric relationship between the radar and the target trajectory. The range of the equilibrium glide HGV detected by radar is more than 1200 km.   By applying the proposed model, the comparisons of EKF, UKF and IEKF on the position and velocity estimation for equilibrium glide trajectory are shown in Fig. 16 and Fig. 17. The simulation results are extracted at 15 seconds interval. Therefore, the details of the lines and points in the diagram can be clearly presented.
For the equilibrium glide scenario, Table 5 shows the ARMSEs and the runtime completed 100 Monte Carlo simulations of the three algorithms. It depicts that the average  tracking precision of the IEKF algorithm is similar to that of the UKF algorithm, which is better than the conventional EKF algorithm. ARMSEs on position and velocity of the IEKF algorithm decrease by 10% and 33%, respectively, compared with the EKF algorithm. Although the tracking accuracy of the IEKF algorithm is slightly lower than that of the UKF algorithm, the runtime is reduced by 52%. In summary, the IEKF algorithm exhibits satisfactory accuracy and real-time performance.
For the equilibrium glide trajectory, the RMSEs on position and velocity estimation in the X, Y and Z-axis applying the proposed model are shown in Fig. 18 and Fig. 19, respectively. These results evident that the proposed model and IEKF algorithm exhibit stable tracking performance for equilibrium glide targets. Fig. 20 shows RMSEs of the estimations on the position, velocity, lift acceleration and drag acceleration for the four filter models, respectively. Fig. 20(a) shows that there is no significant difference in the accuracy of the position estimation for these four estimation models when HGV is in the equilibrium glide mode (during 700 s to 950 s). This phenomenon is reasonable since the position of the HGV can be measured directly. The RMSE of the position estimation by the proposed model can converge to about 50 m. Moreover, Fig. 20 (b) illustrates that the velocity estimation error of the proposed acceleration model is significantly lower than that of the other three models. Since the velocity cannot be measured by the radar, the estimation of the velocity relies on the estimation accuracy of the accelerations. Furthermore, Fig. 20 (c) shows that the total lift acceleration by the proposed model obtains the highest accuracy. In Fig. 20 (d), the estimation of the drag acceleration of the proposed model does not show much higher accuracy. However, it steadily converges to a moderate value. Table 6 presents that the ARMSE of the proposed model is significantly better than that of the other three models. It is found that the position estimation accuracy, velocity estimation accuracy and acceleration estimation accuracy are improved by about 7%, 50% and 20%, respectively.

B. MODE 2: SKIP GLIDE MODE
In this mode, the skip gliding trajectory of the HGV, including the initial altitude of 90 km and initial velocity of 6500 m/s is designed. Moreover, the angle of attack and bank are    considered as two control variables. Fig. 21 shows the height and velocity of the HGV during gliding. Due to the influence of the curvature of the earth, the ground-based radar can only   observe the HGV above the local horizon. Fig. 22 depicts the geometric relationship between the radar and the skip glide target trajectory. It is worth noting that the range of the skip glide HGV detected by radar is more than 1700 km and the observable time is about 300 s.
By applying the proposed model, the comparisons of EKF, UKF and IEKF on the position and velocity estimation for skip glide trajectory are shown in Fig. 23 and Fig. 24. The simulation results are extracted at 15 seconds interval. For the skip glide scenario, Table 7 shows the ARMSEs and the runtime of the three algorithms. The ARMSEs on position and velocity of the IEKF algorithm decrease by 11% and 16%, respectively, compared with the EKF algorithm. Indeed, using UKF algorithm can improve the tracking accuracy of skip glide trajectory to some extent, however, the computational burden of the IEKF algorithm is much less than the UKF algorithm. For skip glide trajectory, the RMSEs on position and velocity estimation in the X, Y and Z-axis applying the proposed model are shown in Fig. 25 and Fig. 26, respectively. Since the HGV performs skip maneuver between 800s and 900s, the RMSE of position and velocity estimation increases to about 100 m and 50m/s respectively. And then, the filter converges quickly, the proposed model achieves satisfied estimation precision. Fig. 27 shows RMSEs on the position estimation, velocity estimation, lift acceleration estimation and drag acceleration estimation of the four models. Fig. 27 (a) and Fig. 27 (b) illustrate that the position and velocity estimation accuracy in the abovementioned four filter models is close to each other, and the proposed model shows a reasonable superiority in the estimation accuracy. Moreover, the RMSE of the position estimation by the proposed model can converge to about 30 m and the RMSE of the velocity estimation converges to about 20 m/s. Furthermore, Fig. 27 (c) and Fig. 27 (d) show RMSEs of the acceleration estimation of the four models. It is observed that the tracking precision on the accelerations of the proposed model is more reasonable than the other three models.  Table 8 presents the tracking results quantitatively. It is found that the proposed model improves the position estimation accuracy, velocity estimation accuracy and acceleration estimation accuracy by about 10%, 40% and 50%, respectively.
Consequently, the Gauss model describes the nonmaneuvering motion of the target, which is significantly different from the real HGV motion. While the linear Markov model is suitable for the maneuver motion between constant velocity and constant acceleration, its estimation precision decreases greatly if the target performs a type of strong maneuver. For the CS model, although it can describe the current acceleration accurately, the limitations of the linear Markov model still exist.

VII. CONCLUSION
Motion modes of HGV are complex and uncertain. Therefore, the conventional estimation models cannot describe the HGV motion characteristics accurately. In the present study, the nonlinear characteristic of the gliding dynamics is considered. Then, differentials of the maneuvering acceleration components are modelled by analyzing the statistics of typical maneuvers and configurations. Moreover, an integrated nonlinear Markov acceleration model is formulated to match the nonlinear dynamics of the HGV. In order to obtain the primary information of the aerodynamic characteristics of the HGV, a class of HGV configurations is analyzed and a statistical model for the aerodynamic characteristics is established. Finally, IEKF is used to estimate the trajectory with the proposed model. The simulations for two maneuver modes of HGV are performed. It is concluded that the proposed model is able to improve the trajectory estimation accuracy effectively.
SHUO TANG received the B.S., M.S., and Ph.D. degrees from Northwestern Polytechnical University, Xi'an, China, in 1982China, in , 1985China, in , and 1988, respectively. He is currently a Professor with the School of Astronautics, Northwestern Polytechnical University, Xi'an. His current research interests include aerospace vehicle design, flight dynamics, flight simulation, and virtual prototype technology.