A Midcourse Guidance Method Combined with Trajectory Prediction for Antinear-Space-Gliding Vehicles

,


Introduction
With a large airspace, a high velocity, and strong maneuverability, near-space hypersonic vehicles have brought about new challenges [1][2][3] in the precise collision of antinearspace targets. For a high-velocity near-space-gliding vehicle (NSGV), the proportional navigation (PN) and its variants methods [4][5][6] of interceptors are prone to overload saturation caused by the divergence of the line-of-sight (LOS) angular rate at the end of the interception moment, which will result in a failure to intercept or a full deflection of the actuator. The midcourse guidance of interceptors can effectively reduce the overload at the terminal stage, if a certain intersection angle constraint is satisfied during shift handover at the predicted point of interception [7], such as reverse flight conditions. Moreover, the available overload from the interceptor is limited by the flight conditions in near space; reducing the full-range overload and ensuring the convergence of the terminal overload are necessary. Thus, a midcourse guidance method combined with trajectory prediction, which is subject to multiple constraints such as the seeker field of view, the overload limit, and the intersection-angle constraints, is one of the ways to effectively intercept hypersonic vehicles and provides the optimal flight platform for terminal accurate collision. Various approaches have been proposed to fulfill this capability.
Reference [8] classifies trajectory prediction into three types: the state estimation model, the kinetic model, and the machine-learning model, and then the definition and basic process of trajectory prediction are introduced for the four-dimensional trajectory prediction of aircraft. But, due to the complexity of the maneuver methods of near-space-gliding aircrafts and the fact that they cannot be represented as a single maneuver model, the interacting multiple model (IMM) containing multiple maneuver models is widely used for trajectory tracking of such targets [9]. A novel interacting multiple model (NIMM) filter [10] is presented that combines two different algorithms, the adaptive fading Kalman filter and the maximum correntropy Kalman filter, to track targets in hybrid situations. Target tracking was performed for the maneuverable model in two dimensions with bearing-only tracking (BOT) method and using interacting multiple model (IMM) filter in [11]. In addition, the stability of the trajectory-tracking method faces problems when the maneuver mode changes. Closed-form approximate solutions [12] are derived to predict 3D atmospheric gliding trajectory of hypersonic glide vehicle accurately under high-maneuver flight condition where the lateral maneuver range can be up to 4000 km. Aiming at the problem of gliding near-space hypersonic vehicle (NSHV) trajectory prediction, a trajectory prediction method [13] based on aerodynamic acceleration empirical mode decomposition (EMD) is proposed. A real-time three-dimensional constrained trajectory prediction method [14] for hypersonic glide vehicle based on 3D flight corridor is developed. To address the difficulty of maintaining high prediction accuracy and short prediction time simultaneously in maneuver trajectory prediction, a maneuver trajectory prediction method [15], that is based on a layered strategy, which combines long-term maneuver unit prediction and short-term maneuver trajectory prediction, was proposed. These methods improve the accuracy and convergence of trajectory prediction by identifying the target guidance law, by establishing a maneuver model, and by judging the flight intention.
Over the past few decades, numerous efforts have been devoted to improving the accuracy and optimality of the midcourse guidance for interceptors with close target velocities. References [16] proposed a new optimal guidance law and verified guidance accuracy and robustness for the time-varying velocity. Various new guidance methods have been proposed to improve the ability to intercept highspeed maneuver targets, including mainly improving target prediction accuracy [17,18], considering more constraints [19,20], and applying the numerical-programming algorithm [21] to reduce deviation. However, it is well known that nonlinear optimal control problems for intercepting maneuver targets is complex and difficult to solve in real time, which can be formulated through a variational calculus approach, accounting for nonlinear system dynamics and accounting for hard terminal constraints. To obtain the analytical closed-loop optimal guidance law, the concept of zero-effort errors that was first defined in Ref [22] must be introduced. To lessen the control effort, the gain is optimized by a novel particle swarm optimization (NPSO) algorithm. The novel reentry guidance law [23] with constrained impact angle is proposed for hypersonic weapons using a novel particle swarm optimization (NPSO) algorithm for less control effort. A hybrid particle swarm optimization-(PSO-) gauss pseudomethod (GPM) algorithm [24] with a fast convergence and high precision is proposed to deal with the reentry trajectory optimization problem [22] must be introduced. The characteristics of NSGV targets undermine the advantages in velocity and overload from interceptors, so the midcourse guidance methods not only require the requirements of real-time, multiconstraints, and high precision to be met but also needs to be combined with predictions of the trajectory to improve the interception feasible area.
In this paper, the midcourse guidance method takes NSGVs with various typical maneuver modes as the interception target and guides the interceptor to a feasible interception area at a certain intersection angle. The trajectory prediction method is used to provide terminal constraints for midcourse guidance for hypersonic targets in a space nearby. The midcourse guidance adopts an improved method based on the zero-error-miss/zero-error-velocity (ZEM/ZEV) guidance law to achieve multiple terminals with intersection angle restrictions. The main contributions of this study are as follows: (1) A series of target maneuver models is constructed with strong adaptability and clear physical meaning, based on prior information from the flight mission. Then, the IMM algorithm is improved by introducing the fading factor and by modifying the Markov probability transition matrix, so that improvements in these filtering methods enhance the trajectorytracking ability of NSGVs (2) The characteristic maneuver parameters of the trajectory-prediction method are periodically updated by the adaptive grid combined with the maneuver model to improve the accuracy of the trajectory prediction. The prediction results of target NSGVs at different times are obtained within the conditions of a numerical integration and provide the terminal constraints for the midcourse guidance (3) An improved ZEM/ZEV midcourse method combined with dynamic update of predicted trajectory is applied to meet the constraints of three-dimensional intersection angle and collision conditions and reduce the change rate of terminal overload. The guidance commands are solved in a closed loop according to an updated prediction of the trajectory, which provides the appropriate conditions for an interception to induce an accurate collision with the interceptor in the terminal guidance stage

The Equations of Midcourse Guidance Problem in NSGV Interception
The purpose of this paper is to perform a high-precision collision interception against NSGVs based on predictions of the trajectory combined with multiconstrained midcourse guidance. To this end, interceptor and NSGVrelated models, such as the dynamic model, the midcourse guidance model, and the trajectory-prediction model, need to be established.

2
International Journal of Aerospace Engineering

Dynamic Model of Interceptor and Ballistic
Missile. An interceptor with a dual-pulse engine was investigated in this paper, and the following assumptions are made to obtain the appropriate guidance equations: (1) The engine exhaust velocity and the normalized mass burn rate are constants (2) A nonrotating sphere model without J 2 perturbation is used to describe the relationship between the state parameters, the aerodynamic acceleration, and other variables Based on the above assumptions, the three-dimensional point-mass equations of motion in near space can be expressed as follows: where r i and v i are state vectors; x ib is the unit vector defining the longitudinal direction of the rocket body; y ib is the normal direction of the rocket body; z ib is the lateral direction of the rocket body; g ðr i Þ is a function of gravitational acceleration with r i as the independent variable; A i , N i , and Z i are axial aerodynamic force, normal aerodynamic force, and lateral aerodynamic force, respectively; T i is the thrust magnitude; and m i is the mass of the aerocraft. When the subscript i in Equation (1) is taken as M, it represents the interceptor, and when it is taken as T, it represents NSGV. The thrust magnitude T i is given by the following: where V ex is the engine exhaust velocity; _ m is the normalized mass burn rate of the interceptor; s m is the cross-sectional area of the engine nozzle; s T ig is a signum function whose form will be given later; P ðh M Þ is a function of atmospheric pressure with flight altitude of the interceptor h M as the independent variable. The aerocraft mass m i follows the following differential equation: The axial, normal, and lateral aerodynamic forces A i , N i , and Z i are given as follows: where S i is the aerodynamic reference area; q i is the dynamic pressure; C iA , C iN , and C iZ are the coefficient of axial force, the coefficient of normal force, and the coefficient of lateral force, respectively. x ib , y ib , and z ib can be expressed as the pitch angle φ i , the yaw angleψ i , and the roll angle γ i in inertial reference frame as follows: where φ i , ψ i , and γ i are the three Euler angles representing the rotation of the missile body relative to the inertial reference frame [25]. The interceptor adopts a two-stage solid motor and the interception flight process of interception includes a boost phase and a glide phase. Additionally, s T ig is a signum function representing the working state of the engine, which is expressed as follows: In the boost phase, thrust and aerodynamic force together produce a large acceleration to change the flight trajectory; in the glide phase, the available acceleration gener-ated by an aerodynamic force is limited by the flight conditions [2], which can be expressed as follows: where a i = ½ðT i − A i Þ × x ib + N i × y ib + Z i × z ib /m i ðtÞ; C iA:max ðMa, HÞ, C iN:max ðMa, HÞ, and C iZ:max ðMa, HÞ represent the maximum coefficient of axial aerodynamic force, the maximum coefficient of normal aerodynamic force, and 3 International Journal of Aerospace Engineering the maximum coefficient of lateral aerodynamic force with the corresponding Mach number and altitude, respectively.

Target-Maneuver Trajectory Prediction.
The trajectory prediction of near-space targets mainly solves the problem of predicting the state of a future target flight under the conditions of uncertainty in the measurement information and unknown target maneuvers. Let Z k = fZ 0 , Z 1 , LZ k g represents the target measurement sequence from time 0 (corresponding tot 0 ) to time k (corresponding to t k ) obtained by the ground-based radar; letX kÎ R 6 represents the estimation of target state X T k = ½r T T ðt k Þ, v T T ðt k Þ T at the time t k ; let _ XðtÞ = f T ðXðtÞ, UðtÞ, tÞ represents the target dynamic equation corresponding to (1), in which UðtÞ represents the target control function to be inferred. Thus far, the expected output of the target trajectory prediction problem is the target state predictionXðt p Þ from t k to a certain time t p in the future, determined by Z k , or a target state prediction trajectoryXðtÞðt k ≤ t ≤ t p Þ during a period.
Considering the uncertainty of measurement information Z k and target maneuver UðtÞ, the common practice in predicting the trajectory is to infer the density of the target control function pðUðtÞjZ k Þ through known measurement information and then to calculate the target state prediction trajectories through one or more high-probability control functions U M i The trajectory-prediction method constructs the target maneuver model set M k according to the maneuver mode of glider (NSGVs) takes the data sequence measured by radar as input Z k , estimates the target stateX k , and infer the model probability μ i k = pðM i k ðω i k ÞjZ k Þ using the control quantity U M i k ðω i k Þ ðtÞ determined by each model M i k ðω i k Þ and then modifies the parameters ω i k of the maneuver model according to the model probability. Finally, the predicted flight state of each model at time t p is calculated through (8), in which the one with the maximum model probability used as the predicted intercept point for midcourse guidance. The detailed design process of the trajectory prediction problem is as follows.
To infer the maneuver law of the target, the maneuver model set M k is designed to express different maneuver modes, whose element M i k ðω i k Þ represents the ith model at time k and ω i k is the maneuver parameter that needs to be inferred. Then, the expression of the maneuver mode set at time k is as follows: where N represents the total number of maneuver models used. Since the radar cannot directly obtain the state X k of NSGVs through measurement, it is necessary to estimate the target stateX k|M i According to N different parameters of maneuver model ω i k , N different target state estimation, and the model probability Therefore, the estimation of target state can be calculated aŝ To estimate and predict the target flight trajectory more accurately, the highest probability model parameter is used to gradually modify other low-probability model parameters ω i k so that the trajectory cluster calculated by the maneuver model set M k gradually becomes closer to the observed trajectory. The updated function of the maneuver parameters expressed by the model probability can be expressed as follows: The control function U M i k ðω i k Þ ðtÞ calculated with the maneuver model parameter ω i k is integrated according to the dynamic differential equation in (8) to obtain a prediction of the target trajectory. The target state prediction of t p during interception is expressed as follows: where t k represents the time when the prediction starts; that is, the time corresponding to the current time k. A series of target-predicted trajectory envelopes from t k to t p under each maneuver mode M i k ðω i k Þ can be expressed as follows: In summary, the NSGV trajectory-prediction process includes the design of a maneuver model set, the evaluation of the target state estimation and model probability, the adjustment of the maneuver model parameters, and the calculation of the target prediction-trajectory envelope.

Midcourse Guidance Model.
To unify and simplify the derivation of the guidance equation, the specific force 4 International Journal of Aerospace Engineering acceleration a is introduced to replace Equation (1), which can be compactly expressed as follows: The intermediate guidance problem [22] of the interceptor is used to solve the optimal control with constraints, and the performance index is usually designed to minimize the total overload. This paper adopts a performance index with a weight varying over time, which can be expressed as follows: where t f is the shift handover time for midcourse guidance and terminal guidance, generally speaking, thus ensuring the space-time relationship for intercepting a collision; shift handover time t f is equal to the target trajectory-prediction time t p ; and n is a nonnegative integer. This performance index minimizes the energy used to control the entire process with time-varying weights, which can effectively reduce overload in the whole process and ensures convergence of the overload, and this effect is more obvious with an increase in n. Generally, when n = 0, only the performance index that consumes the least energy is represented. The states of the interceptor at the current guidance time are directly provided by the navigation equipment, and the terminal state is calculated by the trajectory-prediction method according to the target information measured by the radar. The midguidance process of the interceptor is used to guide the interceptor from the current state to the predicted point of interception provided by the target trajectory in Equation (12) at time t p . In order to make the terminal guidance more conducive to an accurate collision with the target aircraft, the terminal condition ψ½r M ðt f Þ, v M ðt f Þ = 0 of the midcourse guidance, shown in Figure 1, must meet the collision triangle relationship and intercept the target at certain intersection angles η θ and η ψ : which have significant impacts on the distance of a miss in the terminal guidance. Then, the velocity vector of the interceptor at the terminal time should meet the following: According to the geometric relationship of the collision triangle [26], the expression of the position vector at the terminal time under the condition of zero control interception is as follows: When the missile target distance reaches the maximum detection distance R max of the interceptor seeker, the interceptor is considered to have completed the shift handover for intermediate and terminal guidance and directly enters the terminal guidance phase, so the shift handover distance constraint is as follows: where from the trajectory-prediction method in Equation (12), and As a result, the midcourse guidance model of the interceptor can be expressed as follows:

Trajectory Prediction by IMM-ACKF
According to the establishment of the maneuver mode set in Section 2.3, the IMM method is used to evaluate the model probability of each maneuver mode (see Sections 3.1 and 3.2), and then, the target state is estimated on the basis of Figure 1: Geometric relationship to be satisfied in the collision between an interceptor and a target.
When designing the set of models, converting the implicit mapping relationship of the flight mission into acceleration is necessary. The mapping relationships involved in the conversion process include drag acceleration velocity mapping [27,28], altitude velocity mapping [29], and lateral flight boundary [30,31]. In order to facilitate establishing a target maneuver model, the acceleration commands are uniformly expressed as UðtÞ = ð−D, L Y , L Z Þ T /m T , and the acceleration generated by each component composed of drag, lift, and laterals force can be sorted as follows: where R L represents the remaining range; V represents the current speed; V f is the expected terminal speed; γ is the current intersection angle; γ′ðVÞ is the function of the correlation coefficient of the altitude velocity mapping relationship [29]; g is the gravitational acceleration; N L/D is the lift-drag ratio; and Sigðδ Z Þ represents the lateral boundary turnover function, which is defined as follows: where δ Z is the lateral flight line of a sight error, Sig Last represents the symbolic calculation result of the last time, and δ c is the boundary threshold related to the lateral boundary parameter V S . In the design of the abovementioned maneuver models, the model parameters that can be selected include the correlation coefficient a 2 of the altitude-velocity mapping relationship, the lift-drag ratio parameter N L/D , and the parameter V S expressing the lateral boundary. The design parameter vector of the ith maneuver mode M i 0 ðω i 0 Þ can be selected as follows: For the design of the aircraft maneuver mode set at the initial time (k = 0), the components of the parameter vector take the minimum, middle, and maximum values within the feasible range, which are combined to form a total of N = 33 = 27 feasible maneuver modes.
When intercepting NSGVs, a high-power active-passive composite radar is generally used to monitor and track the target remotely. At this time, the radar can measure the dis-tanceR r from the target, the target azimuth λ T , and the highlow angle λ D and can obtain the change rate _ R r of the target distance from the Doppler principle. As shown in Figure 2, the radar measurement equation of NSGVs is given.
where Δr = ½Δr x , Δr y , Δr z T is the position vector of the target relative to the radar and Δv = ½Δv x , Δv y , Δv z T is the velocity vector of the target relative to the radar.
Some errors and noises are found between the actual radar observation and the real state of the target. Equations (22) and (25) can be expressed as follows: where e z = ½n R n D n T n _ R T is the simulated radar observation noise and the matrix parameters are closely related to the radar performance. Radar measurement error can usually be considered Gaussian white noise with zero mean. where Nð0, R ν Þ represents a random variable e z that obeys zero mean and covariance R ν . Generally, the observation noise of each channel is believed to be independent of each other, and the radar observation noise covariance matrix can be obtained as follows: The diagonal elements represent the covariance in the measurement error of each channel. The establishment of a dynamic equation is made on the basis of maneuver model modelling, and the establishment of a radar measurement coordinate system is made on the basis of a radar pseudoobservation measurement and a filtering equation. The result of the measurement in Equation (25) from time 0 to time k constitutes the target measurement information sequence Z k = fZ 0 , Z 1 ,⋯Z k g.

Maneuver Target Tracking
Method. According to the multiple maneuver models M k in Equation (24) from Section 3.1 and the measurement information sequence Z k of the target, the IMM method was used to evaluate the model probability of each maneuver mode, and the probability μ i k = pðM i k ðω i k ÞjZ k Þ of each maneuver model was obtained. Because the maneuver mode of the target cannot be well determined in advance, the single-motion model cannot accurately describe the flight state of the target. The interactive multimodel algorithm uses the combination of multiple motion models to estimate the system state through effective weighted fusion. The specific principle of the combination of IMM [9] and adaptive culture Kalman filter (ACKF) is shown in Figure 3.
The maneuver model set covers all possible typical maneuver modes; the deviation from the real motion model of the target is large. To prevent divergence in the filtering process, the fading factor is introduced into the state prediction covariance matrix, the gain matrix is adjusted online and in real time, and the output residual sequence is forced to be orthogonal. The calculation process of fading factor is shown in Equation (29).
where Q k , R k+1 are the variances in system noise and measurement noise, respectively; P x k+1/k , P l k+1/k are the one-step prediction error covariance matrix without and with system process noise in CKF; ω i . is the weight of the volume point; χ l k+1/k,i ,X k+1/k , γ k+1 are the volume point, state estimation, and measurement residual generated by one-step state prediction without considering fading factor;ŷ l k+1 is the values of predicted measurement obtained by weighting χ l k+1/k,i ; h is the measurement functions, as shown in Formula (25); κ ≥ 1 is the weakening factor, which is generally selected according to experience; and 0:95 ≤ ρ ≤ 0:995 is the forgetting factor. In the formula, H k+1 , Μ k+1 , N k+1 is the calculation process parameter matrix. The covariance matrix is modified by the fading factor to obtain the one-step prediction state covariance matrix:

International Journal of Aerospace Engineering
After introducing the fading factor,X k+1 conforms to the distribution NðX k+1/k , P k+1/k Þ. The predicted stateX k+1/k obtained by one-step prediction and the one-step predicted state covariance matrix P k+1/k with a fading factor are filtered as the input of CKF to posteriorly estimate the state quantity, to update the covariance matrix, and to form an adaptive culture Kalman filter with the fading factor. Thus far, the posterior probability density function pðX k jZ k , M j k ðω j k ÞÞ is calculated based on the Gaussian hypothesis. The likelihood function value of each model is as follows: where v j k = Z k − hðX l kjk−1 Þ, s j k = P zz kjk−1 , and Z k are the observation at time k, hðX l kjk−1 Þ is the observation prediction value before introducing the fading factor, and P zz kjk−1 is the observation prediction covariance matrix for ACKF before calculating the fading factor.
The probability-estimation process in IMM-ACKF is described as below.
3.2.1. Interaction/Mixing. Mixed probability of model i to model j is calculated as follows: where c j is calculated as: The mixed state and mixed covariance of the corresponding model are estimated as follows: where c is calculated as The state and covariance estimation describe the form of a posteriori probability density function pðX k jZ k Þ for the target state under a Gaussian assumption. The above state estimationX k ≜X k|k from time 0 to time k constitute the estimation of the target state sequenceX k = fX 0 ,X 1 ,⋯X k g.

Maneuver Mode Set Adjustment
Strategy. When many model parameters need to be designed, the scale of the model set is limited to ensure operation efficiency. Under the conditions of the same model set scale, to improve the coverage of the model set to various maneuver modes, an adaptive grid is introduced to adjust the maneuver model parameters in real time. In the standard IMM algorithm, the transition probability is preset and cannot be adjusted in real time according to posterior information. However, due to the uncertainty of target mobility and prior information, the preset Markov matrix cannot well reflect the switching between target real-motion modes, resulting in an increase in estimation errors. Therefore, bridging the gaps between the weights of each model and their real values according to the likelihood function is necessary [30].
Assuming that the current maneuver model parameters (such as N L/D , V s ) of the maneuver target are in a continuous range ½ω min , ω max that the model set parameters at time k are ½ω L k , ω c k , ω R k , and that the posterior probabilities corresponding to the three model set parameters at time k are ½μ L k , μ c k , μ R k ; the model set of parameters at the time k + 1 is determined as follows.
An adjustment in the model set center parameters is obtained by taking the model posterior probability as the weight coefficient.
The inaccuracy of the system model, the influence of observation errors, and the competition between models cause the maximum a posteriori probability of the model to not always appear in the model closest to the real model. Therefore, continuous jumps, shaking, and even divergence 8 International Journal of Aerospace Engineering of the central model around the real model space are easily produced. A grid center adjustment is thus needed; (1) To prevent irregular jitter of the model set, we set the grid center change threshold t g . When the maximum posterior probability of the model exceeds the threshold, adjust the parameters of the central model ω C

k+1
(2) To prevent excessive drastic changes in the center of the grid and to reduce the total distance of the grid center when moving, a smoothing factor is introduced

<
: where t 1 is the posterior probability threshold. If the posterior probability of the models on the left and right sides is less than t 1 , the model parameters do not match the current maneuver mode and thus reduce the interval from the grid center.
(4) If μ L k = max fμ L k , μ C k , μ R k g, then the right model keeps the original interval unchanged and the left model adjusts the model interval according to the threshold value The above rules for parameter updates in the maneuver model constitute the correction law of model probability in The value of the likelihood function of the degree of the model reflects the match between the model and the real motion pattern of the target. The value of the likelihood function at time k for model j is calculated according to (31). A larger value relative to the likelihood function value of other models j indicates that the model better matches the real motion state of the target and thus has a greater probability of other models transitioning to this matching model.
Let the probability of model i transferring to model j at time k be π ij k and satisfy ∑ M j=1 π ij k = 1. The elements of the Markov probability transition matrix are modified by the following methods: where ξ = 0:5 is used to control the speed of adjustment. The corresponding probability transition matrix is transformed into the following: These methods may cause some values of the main diagonal elements of the matrix to become smaller and smaller; the probability of transferring the mismatched model to itself becomes smaller and smaller. However, when the target maneuvers to turn the mismatched model into a matching model, the model switching is slow. Therefore, we set a threshold T h = 0:5 for the main diagonal element if π ij k,1 < T h , 0 < T h < 1, then π ij k,1 = T h , which thus has the following form:

Improved ZEM/ZEV Midcourse Guidance Method
Guiding the interceptor to the zero-control interceptor manifold can effectively reduce terminal guidance, which is one of the most effective midcourse guidance methods. The midcourse guidance algorithm based on ZEM-ZEV mainly studies the relationship between guidance command and zero-effort-miss distance and zero-effort-velocity. The specific derivation is as follows.
4.1. Principle of ZEM/ZEV. The guidance law can be obtained by solving the two-point boundary value problem of Equation (15), and the Lagrange multiplier is introduced to combine the terminal constraint with the performance index to construct the corresponding Hamiltonian function. International Journal of Aerospace Engineering (6) Calculation of interceptor closed-loop guidance command: since the midcourse guidance method needs to calculate ZEM and ZEV in Equation (50) but the terminal constraint lacks speed, it needs to be calculated using an additional iterative method. Taking the flight time t p and the terminal speed v M ðt f Þ of the interceptor as the quantity and the maximum terminal speed as the performance index J = −v M ðt f Þ, the individual optimal position p i,j , and the global optimal position g j are calculated according to the particle swarm optimization algorithm, and the iteration is carried out in the following form where r 1,2 is the random number used when updating particles; c is the proportion of global optimization and particle optimization; and ε = ½t p , v M ðt f Þ T , iterating continuously until the accuracy requirements are met, i.e., kg j − g j−1 k ≤ kΔεk, and Equation (50) outputs the guidance instruction a (7) Calculate the start time of the boost phase t ig : estimate the remaining flight time t go ðt ig Þ which satisfies ðv T ðt + t go ðt ig ÞÞ − v M ðt + t go ðt ig ÞÞÞ · ðr T ðt + t go ðt ig ÞÞ + R − r M ðt + t go ðt ig ÞÞÞ ≥ 0 and calculate t ig which satisfies the remaining flight time constraint, i.e., kΔt go ðt ig Þk = kt f − t − t go ðt ig Þk = 0. Newton-Raphson iterative method is an effective approach used to solve the problem. The iterative variable is formulated as follows:

12
International Journal of Aerospace Engineering Iterating continuously until the accuracy requirements are met, i.e., kt ig,j − t ig,j−1 k ≤ kΔεk.

Simulation Results.
In Section 3, the NSGV maneuver model set was established based on target intent, the radar observation model was used to track the trajectory of NSGVs, and the interactive multimodel variable structure filtering method predicted feasible target trajectories. Then, in Section 4, the midcourse guidance method satisfying the terminal constraint of intersection angle guided the interceptor to collide with the target. In this section, to evaluate the performance of the proposed guidance and trajectoryprediction method, several cases are applied.

Performance of Trajectory Prediction Method.
To verify the proposed trajectory prediction method, the true trajectory maneuver is generated using the method mentioned in [32], which starts from 35 s to 78 s. The initial height of target gliding segment is about 22 km, the terminal height is 14 km, and the gliding distance is about 100 km. In target trajectory prediction, the constant coefficients used in control law response to drag profile deviations are unknown, as well as the coefficients in the roll angle command equation and the angle-of-attack command equation. The effect of the above coefficients on the target trajectory is estimated approximately by the proposed model adjustment method mentioned in Step 3 of section 4.2. According to the set of maneuver models for the designed target in Equation (24)   14 International Journal of Aerospace Engineering and the mean square deviation (1σ) of the angle of observation noise n D , n T of X and Y is 0:1 ∘ in the radar observation noise covariance matrix of Equation (28). After estimating the state of the target and model probability correction, the target trajectory obtained by trajectory numerical integration and extrapolation are shown in Figure 5.
The simulation results in Figure 5 show obvious deviations in trajectory prediction in the longitudinal plane caused by the fixed maneuver model. For comparison, the adaptive grid in Equation (38) updates the parameters of the maneuver model set according to the measurement of the target trajectory by the radar to bring the predicted trajectories integrated by the maneuver model set in Equation (24) close to the true trajectory. The deviation in predictions of the trajectory at different times is plotted in Figure 6. Figure 6 illustrates the trajectory prediction deviation obtained by the proposed IMM-ACKF method under the observation conditions of 0 s, 5 s, 15 s, and 25 s. The blue line in the figure is the predicted trajectory deviation corresponding to the maneuver model with high probability, and the corresponding red line is the trajectory-prediction deviation corresponding to the maneuver model with low probability. Obviously, with the increase in observation time and the decrease in prediction time for the response, the amount of deviation in the prediction increases first and then decreases, which shows a certain convergence because of the flight intention and constraints of the terminal mission for NSGVs. Figure 6(a) shows that the target trajectory of the NSGVs is not observed and that the prediction parameters are not updated, which is consistent with the parameters of the fixed maneuver model in Figure 5. Additionally, the fixed maneuver model set fails to cover the true trajectory, so the prediction deviation is larger than the adjusted maneuver model updated by IMM-ACKF trajectory prediction. In Figure 6(b), the trajectory predicted by some high-priority maneuver models has a large deviation compared with the low-priority models primarily because the real target overturned the pitch angle and the model set was not updated in time. Therefore, the set of trajectories predicted by the high-probability models still needs to be updated and adjusted many times according to the actual radar measurement data to reduce the amount of deviation from the prediction of the trajectory. Table 1 shows the simulation results of trajectory prediction accuracy deviation calculated by the proposed IMM-ACKF method in the same prediction duration of 15 seconds under different observation intervals and different initial prediction times. Table 1 shows that the trajectory prediction error of the proposed method decreases with the increase of the observation duration. And the predicted trajectory generated by the maximum probability model is close to the average value predicted by all other models in Equation (24). Therefore, the predicted trajectory generated by the maximum probability model can be used as the predicted value of the terminal constraint of the midcourse guidance method.

Comparison and Analysis of Midcourse Guidance
Methods. We designed a multiconstraint interception task based on predictions of the trajectory to verify the effectiveness of the algorithm. According to the target state information obtained by the trajectory-tracking method based on the guidance law, the experiment comparing interceptions under nominal conditions was designed: the interceptor was launched at 0 s, and the whole flight time lasted 75 s, which is consistent with the time predicted by the trajectory. The maximum detection range of the seeker R max was 20 km, and the constraints of the track angles η θ and η ψ were 7°and −3°in the reverse flight conditions, respectively. At the same time, the methods proportional navigation with the trajectory prediction method (PN with TP) and proportional navigation without trajectory prediction (PN without TP) were introduced for a comparison with the methods proposed in this paper. The simulation curve of the interception flight parameters based on the prediction information of the target trajectory is shown in Figure 7. In Figure 8, the change in the position of the interceptor when considering predictions of the trajectory of NSGVs is shown in (a), and the trajectories of interception during terminal guidance phase using different midcourse guidance methods is shown in (b). Figure 8 Table 2, it indicates that the proposed IMM-ACKF trajectory prediction method can bring the results of the traditional PN method closer to the target, thus leading to a perfect interception after passing the handover point shown in green cross in (a) and (b). In Figure 7(a) and (b), the overload command of PN with the TP method continues to increase at the end of the flight time, which makes meeting the guidelines for the available overload from the interceptor difficult. However, the overload convergence of the proposed ZEM/ZEV method is gradually enhanced with an increase in weight coefficient n, ensuring that the end overload

16
International Journal of Aerospace Engineering command approaches zero. Figure 7(c) and (d) show that the intersection angles η θ and η ψ of the two directions gradually reach the expected values according to the proposed improved ZEM/ZEV method and satisfy the multiple terminal constraints of the interception in Equation (20). In addition, Figure 7(a), (b), and (f) show that the proposed midcourse guidance method makes the terminal overload quickly converge to zero with an increase in n but reduces the terminal velocity. Therefore, the weight coefficient n must be designed to take into account specific problems and can generally be considered as n = 1. The simulation result in Table 2, compared with PN using the TP method and PN without using the TP method, the deviations in the third and fourth constraints are too large to meet the expected accuracy because the intersection angle is not considered in the design of these two methods. The proposed improved ZEM/ZEV method satisfies multiple constraints composed of the control interception manifold and intersection track angle, including the five equations of ψ½r M ðt f Þ, v M ðt f Þ = 0 established by the terminal velocity direction and the terminal position vector, to guide the interceptor to accurately collide with the target NSGVs.

Conclusions
According to prior information about the target, this paper constructed a maneuver model based on the target intent. Using the improved variable structure IMM-ACKF method, the adaptability of the trajectory prediction method to the target maneuver mode was improved and prediction of target trajectory set was obtained.
To intercept maneuvering glider, we introduced target trajectory prediction parameters into the midcourse guidance method and designed a zero-control intercept manifold with intersection angle constraints as the midcourse guidance terminal constraint, providing a new way to intercept maneuver targets in near space. The simulation results verify that the trajectory prediction method based on the typical target maneuver model proposed in this paper has high accuracy and strong robustness. However, in an actual attack-defense confrontation environment, the hypersonic glider will increase the number of penetration maneuvers, flight avoidances, and detours to reduce the probability of being intercepted. Therefore, in a subsequent analysis of the prior information statistics and to establish the set of models, further consideration and improvement in the set of targets under penetration conditions will be needed.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this article.