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

Authors: JiaQing Wan; Qian Zhang; Hao Zhang; JiaNing Liu
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

Near-space-gliding vehicles have variable maneuver modes and dramatic changes in their ballistic parameters, which lead to a need to accurately predict an intercept point based on predictions of their trajectories. First, a trajectory prediction method builds a set of time-varying maneuver models based on flight missions combined with an adaptive grid to infer maneuver modes. An interactive multiple-model method of variable structure is proposed to identify the characteristics of the maneuver mode by introducing a fading factor and the modified Markov probability transfer matrix and then predict the trajectory through numerical integration. In the midcourse guidance method, the prediction of the target trajectory is introduced, and the zero-control interception manifold with intersection angle constraints is designed as the midcourse guidance terminal constraint. For the calculation of the starting time of the boost phase, the optimization solution satisfying the remaining flight time constraint is realized based on the Newton-Raphson iterative method. The analytical expression of a guidance command based on zero-error-miss/zero-error-velocity is established on the basis of the optimal control theory to provide an optimal flight path guiding an NSGV fly toward a point of interception. The simulation results show that the trajectory-prediction method has high prediction accuracy and strong convergency for typical maneuver modes, and the proposed midcourse guidance algorithm meets the requirements of the zero-effort intercept manifold with the intersection angle constraints, which is of important theoretical significance and acts as a reference value for intercepting high-velocity maneuver targets.


---

## Body

## 1. Introduction

With a large airspace, a high velocity, and strong maneuverability, near-space hypersonic vehicles have brought about new challenges [1–3] in the precise collision of antinear-space targets. For a high-velocity near-space-gliding vehicle (NSGV), the proportional navigation (PN) and its variants methods [4–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 high-speed 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 trajectory-tracking 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


## 2. 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 NSGV-related models, such as the dynamic model, the midcourse guidance model, and the trajectory-prediction model, need to be established.


### 2.1. 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 J2 perturbation is used to describe the relationship between the state parameters, the aerodynamic acceleration, and other variablesBased on the above assumptions, the three-dimensional point-mass equations of motion in near space can be expressed as follows:(1)ṙi=vi,v̇i=gri+Ti−Ai⋅xib+Ni⋅yib+Zi⋅zib/mi,i=M,T,where ri and vi are state vectors; xib is the unit vector defining the longitudinal direction of the rocket body; yib is the normal direction of the rocket body; zib is the lateral direction of the rocket body; gri is a function of gravitational acceleration with ri as the independent variable; Ai, Ni, and Zi are axial aerodynamic force, normal aerodynamic force, and lateral aerodynamic force, respectively; Ti is the thrust magnitude; and mi 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 Ti is given by the following:
(2)Ti=sTig⋅Vex⋅ṁ−PhM⋅sm,i=M,0,i=T,where Vex is the engine exhaust velocity; ṁ is the normalized mass burn rate of the interceptor; sm is the cross-sectional area of the engine nozzle; sTig is a signum function whose form will be given later; PhM is a function of atmospheric pressure with flight altitude of the interceptor hM as the independent variable. The aerocraft mass mi follows the following differential equation:
(3)ṁi=sTigṁi=M,0i=T.mit0=mi.0.The axial, normal, and lateral aerodynamic forcesAi, Ni, and Zi are given as follows:
(4)Ai=CiAqiSi,Ni=CiNqiSi,Zi=CiZqiSii=M,T.where Si is the aerodynamic reference area; qi is the dynamic pressure; CiA, CiN, and CiZ are the coefficient of axial force, the coefficient of normal force, and the coefficient of lateral force, respectively. xib, yib, and zib can be expressed as the pitch angle φi, the yaw angleψi, and the roll angle γi in inertial reference frame as follows:
(5)xib=cosφicosψi,sinφicosψi,−sinψiT,yib=cosφisinψisinγi−cosγisinφi,cosγicosφi+sinψisinγisinφi,cosψisinγiT,zib=sinφisinγi+cosγicosφisinψi,cosγisinψisinφi−sinγicosφi,cosψicosγiT,i=M,T,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, sTig is a signum function representing the working state of the engine, which is expressed as follows:
(6)sTig=1,s.t.tig<t<tig+Tig0,s.t.other.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 generated by an aerodynamic force is limited by the flight conditions [2], which can be expressed as follows:
(7)Ai≤Aimax=qiSiCiA.maxMa,H,Ni≤Nimax=qiSiCiN.maxMa,H,Zi≤Zimax=qiSiCiZ.maxMa,H,ai2≤qSAimax⋅xib+Nimaxyib+Zimaxzib2/mit,i=M,T,where ai=Ti−Ai×xib+Ni×yib+Zi×zib/mit; CiA.maxMa,H, CiN.maxMa,H, and CiZ.maxMa,H represent the maximum coefficient of axial aerodynamic force, the maximum coefficient of normal aerodynamic force, and the maximum coefficient of lateral aerodynamic force with the corresponding Mach number and altitude, respectively.


### 2.2. 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. LetZk=Z0,Z1,LZk represents the target measurement sequence from time 0 (corresponding tot0) to time k (corresponding to tk) obtained by the ground-based radar; let X^kI^R6 represents the estimation of target state XkT=rTTtk,vTTtkT at the time tk; let Ẋt=fTXt,Ut,t represents the target dynamic equation corresponding to (1), in which Ut represents the target control function to be inferred. Thus far, the expected output of the target trajectory prediction problem is the target state prediction X^tp from tk to a certain time tp in the future, determined by Zk, or a target state prediction trajectory X^ttk≤t≤tp during a period.Considering the uncertainty of measurement informationZk and target maneuver Ut, the common practice in predicting the trajectory is to infer the density of the target control function pUtZk through known measurement information and then to calculate the target state prediction trajectories through one or more high-probability control functions UMkiωkit.
(8)Sset=X^MkiωkitX^Mkiωkit=X^k+∫tktfTXt,UMkiωkit,t,t∈tk,tp.The trajectory-prediction method constructs the target maneuver model setMk according to the maneuver mode of glider (NSGVs) takes the data sequence measured by radar as input Zk, estimates the target state X^k, and infer the model probability μki=pMkiωkiZk using the control quantity UMkiωkit determined by each model Mkiωki and then modifies the parameters ωki of the maneuver model according to the model probability. Finally, the predicted flight state of each model at time tp 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 setMk is designed to express different maneuver modes, whose element Mkiωki represents the ith model at time k and ωki is the maneuver parameter that needs to be inferred. Then, the expression of the maneuver mode set at time k is as follows:
(9)Mk=Mk1ωk1,Mk2ωk2,⋯,MkNωkN,where N represents the total number of maneuver models used. Since the radar cannot directly obtain the state Xk of NSGVs through measurement, it is necessary to estimate the target state X^k∣Mkiωki,Zk according to the target control function UMkiωkit determined by the maneuver model Mkiωki. According to N different parameters of maneuver model ωki, N different target state estimation, and the model probability μkiωki=pMkiωkiZk can be obtained. Therefore, the estimation of target state can be calculated as
(10)X^k=∑i=1NμkiωkiX^kMkiωki,Zk.To estimate and predict the target flight trajectory more accurately, the highest probability model parameterωk∗=argmaxωkiμkiωki is used to gradually modify other low-probability model parameters ωki so that the trajectory cluster calculated by the maneuver model set Mk gradually becomes closer to the observed trajectory. The updated function of the maneuver parameters expressed by the model probability can be expressed as follows:
(11)ωki=gωk∗,ωk−1i,μk−11,μk−12,⋯,μk−1N,i=1,2,⋯N.The control functionUMkiωkit calculated with the maneuver model parameter ωki is integrated according to the dynamic differential equation in (8) to obtain a prediction of the target trajectory. The target state prediction of tp during interception is expressed as follows:
(12)X^Mkiωkitp=X^k+∫tktpfTXt,UMkiωkit,t,where tk 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 tk to tp under each maneuver mode Mkiωki can be expressed as follows:
(13)Sset=X^Mk1ωk1t,X^Mk2ωk2t,⋯,X^MkNωkNttk≤t≤tp.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.


### 2.3. Midcourse Guidance Model

To unify and simplify the derivation of the guidance equation, the specific force accelerationa is introduced to replace Equation (1), which can be compactly expressed as follows:
(14)r¨M=fMt,r,v,a,rMt0=rM.0,vMt0=vM.0.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:
(15)minJ=∫t0tfaTa2tf−tndt,where tf 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 tf is equal to the target trajectory-prediction time tp; 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 tp. In order to make the terminal guidance more conducive to an accurate collision with the target aircraft, the terminal condition ψrMtf,vMtf=0 of the midcourse guidance, shown in Figure 1, must meet the collision triangle relationship and intercept the target at certain intersection angles ηθ and ηψ:
(16)ηθ=π−vMtf⋅vTtfvMtf⋅vTtf,ηψ=atanvMytfvMztf,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:
(17)ψ1=vMtf⋅vTtf−cosηθ,ψ2=vTtf×vMtf×vTtf⋅vTtf×rTtfvTtf×vMtf×vTtf⋅vTtf×rTtf−cosηψ.Figure 1
Geometric relationship to be satisfied in the collision between an interceptor and a target.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:
(18)ψ3=arccosR⋅vMtfR⋅vMtf+arccosR⋅vTtfR⋅vTtf−arccosvMtf⋅vTtf,ψ4=vMtf×R−vTtf×RR.When the missile target distance reaches the maximum detection distanceRmax 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:
(19)ψ5=R−Rmax,where R=rMtf−rTtf and rTtf,vTtfT are obtained from the trajectory-prediction method in Equation (12), and X^tp=tf=rTtf,vTtfT. As a result, the midcourse guidance model of the interceptor can be expressed as follows:
(20)minJ=∫t0tfaTa2tf−tndt,s.t.Equations6,13andψirMtf,vMtf=0,i=1,2,3,4,5.


## 2.1. 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 J2 perturbation is used to describe the relationship between the state parameters, the aerodynamic acceleration, and other variablesBased on the above assumptions, the three-dimensional point-mass equations of motion in near space can be expressed as follows:(1)ṙi=vi,v̇i=gri+Ti−Ai⋅xib+Ni⋅yib+Zi⋅zib/mi,i=M,T,where ri and vi are state vectors; xib is the unit vector defining the longitudinal direction of the rocket body; yib is the normal direction of the rocket body; zib is the lateral direction of the rocket body; gri is a function of gravitational acceleration with ri as the independent variable; Ai, Ni, and Zi are axial aerodynamic force, normal aerodynamic force, and lateral aerodynamic force, respectively; Ti is the thrust magnitude; and mi 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 Ti is given by the following:
(2)Ti=sTig⋅Vex⋅ṁ−PhM⋅sm,i=M,0,i=T,where Vex is the engine exhaust velocity; ṁ is the normalized mass burn rate of the interceptor; sm is the cross-sectional area of the engine nozzle; sTig is a signum function whose form will be given later; PhM is a function of atmospheric pressure with flight altitude of the interceptor hM as the independent variable. The aerocraft mass mi follows the following differential equation:
(3)ṁi=sTigṁi=M,0i=T.mit0=mi.0.The axial, normal, and lateral aerodynamic forcesAi, Ni, and Zi are given as follows:
(4)Ai=CiAqiSi,Ni=CiNqiSi,Zi=CiZqiSii=M,T.where Si is the aerodynamic reference area; qi is the dynamic pressure; CiA, CiN, and CiZ are the coefficient of axial force, the coefficient of normal force, and the coefficient of lateral force, respectively. xib, yib, and zib can be expressed as the pitch angle φi, the yaw angleψi, and the roll angle γi in inertial reference frame as follows:
(5)xib=cosφicosψi,sinφicosψi,−sinψiT,yib=cosφisinψisinγi−cosγisinφi,cosγicosφi+sinψisinγisinφi,cosψisinγiT,zib=sinφisinγi+cosγicosφisinψi,cosγisinψisinφi−sinγicosφi,cosψicosγiT,i=M,T,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, sTig is a signum function representing the working state of the engine, which is expressed as follows:
(6)sTig=1,s.t.tig<t<tig+Tig0,s.t.other.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 generated by an aerodynamic force is limited by the flight conditions [2], which can be expressed as follows:
(7)Ai≤Aimax=qiSiCiA.maxMa,H,Ni≤Nimax=qiSiCiN.maxMa,H,Zi≤Zimax=qiSiCiZ.maxMa,H,ai2≤qSAimax⋅xib+Nimaxyib+Zimaxzib2/mit,i=M,T,where ai=Ti−Ai×xib+Ni×yib+Zi×zib/mit; CiA.maxMa,H, CiN.maxMa,H, and CiZ.maxMa,H represent the maximum coefficient of axial aerodynamic force, the maximum coefficient of normal aerodynamic force, and the maximum coefficient of lateral aerodynamic force with the corresponding Mach number and altitude, respectively.


## 2.2. 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. LetZk=Z0,Z1,LZk represents the target measurement sequence from time 0 (corresponding tot0) to time k (corresponding to tk) obtained by the ground-based radar; let X^kI^R6 represents the estimation of target state XkT=rTTtk,vTTtkT at the time tk; let Ẋt=fTXt,Ut,t represents the target dynamic equation corresponding to (1), in which Ut represents the target control function to be inferred. Thus far, the expected output of the target trajectory prediction problem is the target state prediction X^tp from tk to a certain time tp in the future, determined by Zk, or a target state prediction trajectory X^ttk≤t≤tp during a period.Considering the uncertainty of measurement informationZk and target maneuver Ut, the common practice in predicting the trajectory is to infer the density of the target control function pUtZk through known measurement information and then to calculate the target state prediction trajectories through one or more high-probability control functions UMkiωkit.
(8)Sset=X^MkiωkitX^Mkiωkit=X^k+∫tktfTXt,UMkiωkit,t,t∈tk,tp.The trajectory-prediction method constructs the target maneuver model setMk according to the maneuver mode of glider (NSGVs) takes the data sequence measured by radar as input Zk, estimates the target state X^k, and infer the model probability μki=pMkiωkiZk using the control quantity UMkiωkit determined by each model Mkiωki and then modifies the parameters ωki of the maneuver model according to the model probability. Finally, the predicted flight state of each model at time tp 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 setMk is designed to express different maneuver modes, whose element Mkiωki represents the ith model at time k and ωki is the maneuver parameter that needs to be inferred. Then, the expression of the maneuver mode set at time k is as follows:
(9)Mk=Mk1ωk1,Mk2ωk2,⋯,MkNωkN,where N represents the total number of maneuver models used. Since the radar cannot directly obtain the state Xk of NSGVs through measurement, it is necessary to estimate the target state X^k∣Mkiωki,Zk according to the target control function UMkiωkit determined by the maneuver model Mkiωki. According to N different parameters of maneuver model ωki, N different target state estimation, and the model probability μkiωki=pMkiωkiZk can be obtained. Therefore, the estimation of target state can be calculated as
(10)X^k=∑i=1NμkiωkiX^kMkiωki,Zk.To estimate and predict the target flight trajectory more accurately, the highest probability model parameterωk∗=argmaxωkiμkiωki is used to gradually modify other low-probability model parameters ωki so that the trajectory cluster calculated by the maneuver model set Mk gradually becomes closer to the observed trajectory. The updated function of the maneuver parameters expressed by the model probability can be expressed as follows:
(11)ωki=gωk∗,ωk−1i,μk−11,μk−12,⋯,μk−1N,i=1,2,⋯N.The control functionUMkiωkit calculated with the maneuver model parameter ωki is integrated according to the dynamic differential equation in (8) to obtain a prediction of the target trajectory. The target state prediction of tp during interception is expressed as follows:
(12)X^Mkiωkitp=X^k+∫tktpfTXt,UMkiωkit,t,where tk 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 tk to tp under each maneuver mode Mkiωki can be expressed as follows:
(13)Sset=X^Mk1ωk1t,X^Mk2ωk2t,⋯,X^MkNωkNttk≤t≤tp.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.


## 2.3. Midcourse Guidance Model

To unify and simplify the derivation of the guidance equation, the specific force accelerationa is introduced to replace Equation (1), which can be compactly expressed as follows:
(14)r¨M=fMt,r,v,a,rMt0=rM.0,vMt0=vM.0.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:
(15)minJ=∫t0tfaTa2tf−tndt,where tf 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 tf is equal to the target trajectory-prediction time tp; 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 tp. In order to make the terminal guidance more conducive to an accurate collision with the target aircraft, the terminal condition ψrMtf,vMtf=0 of the midcourse guidance, shown in Figure 1, must meet the collision triangle relationship and intercept the target at certain intersection angles ηθ and ηψ:
(16)ηθ=π−vMtf⋅vTtfvMtf⋅vTtf,ηψ=atanvMytfvMztf,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:
(17)ψ1=vMtf⋅vTtf−cosηθ,ψ2=vTtf×vMtf×vTtf⋅vTtf×rTtfvTtf×vMtf×vTtf⋅vTtf×rTtf−cosηψ.Figure 1
Geometric relationship to be satisfied in the collision between an interceptor and a target.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:
(18)ψ3=arccosR⋅vMtfR⋅vMtf+arccosR⋅vTtfR⋅vTtf−arccosvMtf⋅vTtf,ψ4=vMtf×R−vTtf×RR.When the missile target distance reaches the maximum detection distanceRmax 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:
(19)ψ5=R−Rmax,where R=rMtf−rTtf and rTtf,vTtfT are obtained from the trajectory-prediction method in Equation (12), and X^tp=tf=rTtf,vTtfT. As a result, the midcourse guidance model of the interceptor can be expressed as follows:
(20)minJ=∫t0tfaTa2tf−tndt,s.t.Equations6,13andψirMtf,vMtf=0,i=1,2,3,4,5.


## 3. Trajectory Prediction by IMM-ACKF

According to the establishment of the maneuver mode set in Section2.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 an assumption of each maneuver mode (see Section 3.3) to predict the target trajectory.


### 3.1. Target Maneuver Mode and Status Measured

The target maneuver model setMk provides a series of feasible target control functions UMkiωkit for tracking the target state X^k∣Mkiωki,Zk. 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 Ut=−D,LY,LZT/mT, and the acceleration generated by each component composed of drag, lift, and laterals force can be sorted as follows:
(21)D=2RLV2−Vf2,LY=−γ′V⋅VD+gsinγ−V2R−gcosγ,LZ=SigδZ⋅D⋅NL/D2−LY2,where RL represents the remaining range; V represents the current speed; Vf 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; NL/D is the lift-drag ratio; and SigδZ represents the lateral boundary turnover function, which is defined as follows:
(22)Sig=1,δZ>δc,−1,δZ<−δc,SigLast,δZ≤δc,where δZ is the lateral flight line of a sight error, SigLast represents the symbolic calculation result of the last time, and δc is the boundary threshold related to the lateral boundary parameter VS. In the design of the abovementioned maneuver models, the model parameters that can be selected include the correlation coefficient a2 of the altitude–velocity mapping relationship, the lift-drag ratio parameter NL/D, and the parameter VS expressing the lateral boundary. The design parameter vector of the ith maneuver mode M0iω0i can be selected as follows:
(23)ωki=a2iNL/DiVSiT.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.
(24)Mk=Mkiωkiωki=a2iNL/DiVSiT,a2i∈a2min,a2mid,a2maxNL/Di∈NL/Dmini,NL/Dmidi,NL/DmaxiVSi∈VSmini,VSmidi,VSmaxi.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 distanceRrfrom the target, the target azimuth λT, and the high-low 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.
(25)Zr=RrλDλTṘr=Δx2+Δy2+Δz2arctanΔy/Δx2+Δz2−arctanΔz/ΔxΔxΔẋ+ΔyΔẏ+ΔzΔż/Δx2+Δy2+Δz2,where Δr=Δrx,Δry,ΔrzT is the position vector of the target relative to the radar and Δv=Δvx,Δvy,ΔvzT is the velocity vector of the target relative to the radar.Figure 2
Schematic diagram of a radar measurement of the target trajectory.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:
(26)Z=Zr+ez=Zr+nRnDnTnṘT,where ez=nRnDnTnṘ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.
(27)ez∼N0,Rν,where N0,Rν represents a random variable ez 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:
(28)Rν=σR20000σD20000σT20000σṘ2.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 pseudo-observation 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 Zk=Z0,Z1,⋯Zk.


### 3.2. Maneuver Target Tracking Method

According to the multiple maneuver modelsMk in Equation (24) from Section 3.1 and the measurement information sequence Zk of the target, the IMM method was used to evaluate the model probability of each maneuver mode, and the probability μki=pMkiωkiZk 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.Figure 3
Schematic diagram of variable structure IMM trajectory tracking.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).
(29)λk+1=max1,trNk+1trMk+1,Nk+1=V0,k+1−Hk+1QkHk+1T−κRk+1,Μk+1=Hk+1Pk+1/kxHk+1T,Hk+1=Pk+1/kxy,lT,Pk+1/kxy,l=∑iωiχk+1/k,il−X^k+1/khχk+1/k,il−y∧k+1lT,V0,k+1=γ1γ1T,k=0,ρV0,k+γk+1γk+1T1+ρ,k>0,where Qk,Rk+1 are the variances in system noise and measurement noise, respectively; Pk+1/kx,Pk+1/kl are the one-step prediction error covariance matrix without and with system process noise in CKF; ωi. is the weight of the volume point; χk+1/k,il,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; y^k+1l is the values of predicted measurement obtained by weighting χk+1/k,il; 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, Hk+1,Μk+1,Nk+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:
(30)Pk+1/k=λk+1Pk+1/kx+Qk.After introducing the fading factor,X^k+1 conforms to the distribution NX^k+1/k,Pk+1/k. The predicted state X^k+1/k obtained by one-step prediction and the one-step predicted state covariance matrix Pk+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 pXkZk,Mkjωkj is calculated based on the Gaussian hypothesis. The likelihood function value of each model is as follows:
(31)Λkj=12πn/2skj1/2exp−12vkjTskj−1vkj,where vkj=Zk−hXkk−1l, skj=Pkk−1zz, and Zk are the observation at time k, hXkk−1l is the observation prediction value before introducing the fading factor, and Pkk−1zz 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 modeli to model j is calculated as follows:
(32)μk−1ij=πk−1ijμk−1ic¯j,where c¯j is calculated as:
(33)c¯j=∑i=1Nπk−1ijμk−1i.The mixed state and mixed covariance of the corresponding model are estimated as follows:(34)X^k−1k−10j=∑i=1NX^k−1k−1iμk−1ij,Pk−1k−10j=∑i=1Nμk−1ijPk−1k−1i+X^k−1k−10j−X^k−1k−1i⋅X∧k−1k−10j−X∧k−1k−1iT.


#### 3.2.2. Mode-Matched Filtering Through ACKF

The mixed state and mixed covariance (34) are used as input to the ACKF matched to Mkjωkj, which uses Zk to yield X^k∣Mkiωki,Zk≜X^k∣kj,Pk∣kj. The model probability is updated as follows:
(35)μkjωkj=Λkjc¯jv,where c is calculated as
(36)c=∑j=1NΛkjc¯j.


#### 3.2.3. State Estimate and Covariance Combination

State estimate and covariance combination are done as follows:(37)X^k∣k=∑i=1NX^k∣kiμki,Pk∣k=∑i=1NμkiPk∣ki+X^k∣ki−X^k∣k⋅X∧k∣ki−X∧k∣kT.The state and covariance estimation describe the form of a posteriori probability density functionpXkZk for the target state under a Gaussian assumption. The above state estimation X^k≜X^k∣k from time 0 to time k constitute the estimation of the target state sequence X∧k=X^0,X^1,⋯X^k.


### 3.3. 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 asNL/D,Vs) of the maneuver target are in a continuous range ωmin,ωmax that the model set parameters at time k are ωkL,ωkc,ωkR, and that the posterior probabilities corresponding to the three model set parameters at time k are μkL,μkc,μkR; 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.(38)ωk+1c=ωkLμkL+ωkCμkC+ωkRμkR.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 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 thresholdtg. When the maximum posterior probability of the model exceeds the threshold, adjust the parameters of the central model ωk+1C(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(39)ωk+1C=αωkC+1‐αμkLωkL+μkCωkC+μkRωkR.(3)
IfμkC=maxμkL,μkC,μkR, the intermediate model has the highest probability
(40)ωk+1L=ωk+1C−λkL/2,μkL<t1,ωk+1C−λkL,others,ωk+1R=ωk+1C+λkR/2,μkL<t1,ωk+1C+λkR,others,where t1 is the posterior probability threshold. If the posterior probability of the models on the left and right sides is less than t1, the model parameters do not match the current maneuver mode and thus reduce the interval from the grid center.(4)
IfμkL=maxμkL,μkC,μkR, then the right model keeps the original interval unchanged and the left model adjusts the model interval according to the threshold value
(41)ωk+1L=ωk+1C−λkL2,μkL<t2,ωk+1C−λkL,others,ωk+1R=ωk+1C+λkR.(5)
IfμkR=maxμkL,μkC,μkR, then the left model keeps the original interval unchanged and the right model is adjusted according to the threshold value
(42)ωk+1L=ωk+1C−λkL,ωk+1R=ωk+1C+λkR2,μkL<t3,ωk+1C+λkR,others.The above rules for parameter updates in the maneuver model constitute the correction law of model probability ingωk−1i,pMk−1jωk−1jZkj=1,2,⋯,N.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 timek 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 modeli transferring to model j at time k be πkij and satisfy ∑j=1Mπkij=1. The elements of the Markov probability transition matrix are modified by the following methods:
(43)πk,0ij=λkijπk−1ij;λkij=ΛkjΛkiξ,0≤ξ≤1,where ξ=0.5 is used to control the speed of adjustment. The corresponding probability transition matrix is transformed into the following:
(44)πk,1ij=πk,0ij∑j=1Mπk,0ij=λkijπk−1ij∑j=1Mλkijπk−1ij.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 thresholdTh=0.5 for the main diagonal element if πk,1ij<Th,0<Th<1, then πk,1ij=Th, which thus has the following form:
(45)πkij=πk,1ij−Th−πk,1ijπk,1ij∑j=1Mπk,1ij−πk,1ij.


## 3.1. Target Maneuver Mode and Status Measured

The target maneuver model setMk provides a series of feasible target control functions UMkiωkit for tracking the target state X^k∣Mkiωki,Zk. 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 Ut=−D,LY,LZT/mT, and the acceleration generated by each component composed of drag, lift, and laterals force can be sorted as follows:
(21)D=2RLV2−Vf2,LY=−γ′V⋅VD+gsinγ−V2R−gcosγ,LZ=SigδZ⋅D⋅NL/D2−LY2,where RL represents the remaining range; V represents the current speed; Vf 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; NL/D is the lift-drag ratio; and SigδZ represents the lateral boundary turnover function, which is defined as follows:
(22)Sig=1,δZ>δc,−1,δZ<−δc,SigLast,δZ≤δc,where δZ is the lateral flight line of a sight error, SigLast represents the symbolic calculation result of the last time, and δc is the boundary threshold related to the lateral boundary parameter VS. In the design of the abovementioned maneuver models, the model parameters that can be selected include the correlation coefficient a2 of the altitude–velocity mapping relationship, the lift-drag ratio parameter NL/D, and the parameter VS expressing the lateral boundary. The design parameter vector of the ith maneuver mode M0iω0i can be selected as follows:
(23)ωki=a2iNL/DiVSiT.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.
(24)Mk=Mkiωkiωki=a2iNL/DiVSiT,a2i∈a2min,a2mid,a2maxNL/Di∈NL/Dmini,NL/Dmidi,NL/DmaxiVSi∈VSmini,VSmidi,VSmaxi.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 distanceRrfrom the target, the target azimuth λT, and the high-low 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.
(25)Zr=RrλDλTṘr=Δx2+Δy2+Δz2arctanΔy/Δx2+Δz2−arctanΔz/ΔxΔxΔẋ+ΔyΔẏ+ΔzΔż/Δx2+Δy2+Δz2,where Δr=Δrx,Δry,ΔrzT is the position vector of the target relative to the radar and Δv=Δvx,Δvy,ΔvzT is the velocity vector of the target relative to the radar.Figure 2
Schematic diagram of a radar measurement of the target trajectory.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:
(26)Z=Zr+ez=Zr+nRnDnTnṘT,where ez=nRnDnTnṘ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.
(27)ez∼N0,Rν,where N0,Rν represents a random variable ez 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:
(28)Rν=σR20000σD20000σT20000σṘ2.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 pseudo-observation 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 Zk=Z0,Z1,⋯Zk.


## 3.2. Maneuver Target Tracking Method

According to the multiple maneuver modelsMk in Equation (24) from Section 3.1 and the measurement information sequence Zk of the target, the IMM method was used to evaluate the model probability of each maneuver mode, and the probability μki=pMkiωkiZk 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.Figure 3
Schematic diagram of variable structure IMM trajectory tracking.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).
(29)λk+1=max1,trNk+1trMk+1,Nk+1=V0,k+1−Hk+1QkHk+1T−κRk+1,Μk+1=Hk+1Pk+1/kxHk+1T,Hk+1=Pk+1/kxy,lT,Pk+1/kxy,l=∑iωiχk+1/k,il−X^k+1/khχk+1/k,il−y∧k+1lT,V0,k+1=γ1γ1T,k=0,ρV0,k+γk+1γk+1T1+ρ,k>0,where Qk,Rk+1 are the variances in system noise and measurement noise, respectively; Pk+1/kx,Pk+1/kl are the one-step prediction error covariance matrix without and with system process noise in CKF; ωi. is the weight of the volume point; χk+1/k,il,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; y^k+1l is the values of predicted measurement obtained by weighting χk+1/k,il; 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, Hk+1,Μk+1,Nk+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:
(30)Pk+1/k=λk+1Pk+1/kx+Qk.After introducing the fading factor,X^k+1 conforms to the distribution NX^k+1/k,Pk+1/k. The predicted state X^k+1/k obtained by one-step prediction and the one-step predicted state covariance matrix Pk+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 pXkZk,Mkjωkj is calculated based on the Gaussian hypothesis. The likelihood function value of each model is as follows:
(31)Λkj=12πn/2skj1/2exp−12vkjTskj−1vkj,where vkj=Zk−hXkk−1l, skj=Pkk−1zz, and Zk are the observation at time k, hXkk−1l is the observation prediction value before introducing the fading factor, and Pkk−1zz 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 modeli to model j is calculated as follows:
(32)μk−1ij=πk−1ijμk−1ic¯j,where c¯j is calculated as:
(33)c¯j=∑i=1Nπk−1ijμk−1i.The mixed state and mixed covariance of the corresponding model are estimated as follows:(34)X^k−1k−10j=∑i=1NX^k−1k−1iμk−1ij,Pk−1k−10j=∑i=1Nμk−1ijPk−1k−1i+X^k−1k−10j−X^k−1k−1i⋅X∧k−1k−10j−X∧k−1k−1iT.


### 3.2.2. Mode-Matched Filtering Through ACKF

The mixed state and mixed covariance (34) are used as input to the ACKF matched to Mkjωkj, which uses Zk to yield X^k∣Mkiωki,Zk≜X^k∣kj,Pk∣kj. The model probability is updated as follows:
(35)μkjωkj=Λkjc¯jv,where c is calculated as
(36)c=∑j=1NΛkjc¯j.


### 3.2.3. State Estimate and Covariance Combination

State estimate and covariance combination are done as follows:(37)X^k∣k=∑i=1NX^k∣kiμki,Pk∣k=∑i=1NμkiPk∣ki+X^k∣ki−X^k∣k⋅X∧k∣ki−X∧k∣kT.The state and covariance estimation describe the form of a posteriori probability density functionpXkZk for the target state under a Gaussian assumption. The above state estimation X^k≜X^k∣k from time 0 to time k constitute the estimation of the target state sequence X∧k=X^0,X^1,⋯X^k.


## 3.2.1. Interaction/Mixing

Mixed probability of modeli to model j is calculated as follows:
(32)μk−1ij=πk−1ijμk−1ic¯j,where c¯j is calculated as:
(33)c¯j=∑i=1Nπk−1ijμk−1i.The mixed state and mixed covariance of the corresponding model are estimated as follows:(34)X^k−1k−10j=∑i=1NX^k−1k−1iμk−1ij,Pk−1k−10j=∑i=1Nμk−1ijPk−1k−1i+X^k−1k−10j−X^k−1k−1i⋅X∧k−1k−10j−X∧k−1k−1iT.


## 3.2.2. Mode-Matched Filtering Through ACKF

The mixed state and mixed covariance (34) are used as input to the ACKF matched to Mkjωkj, which uses Zk to yield X^k∣Mkiωki,Zk≜X^k∣kj,Pk∣kj. The model probability is updated as follows:
(35)μkjωkj=Λkjc¯jv,where c is calculated as
(36)c=∑j=1NΛkjc¯j.


## 3.2.3. State Estimate and Covariance Combination

State estimate and covariance combination are done as follows:(37)X^k∣k=∑i=1NX^k∣kiμki,Pk∣k=∑i=1NμkiPk∣ki+X^k∣ki−X^k∣k⋅X∧k∣ki−X∧k∣kT.The state and covariance estimation describe the form of a posteriori probability density functionpXkZk for the target state under a Gaussian assumption. The above state estimation X^k≜X^k∣k from time 0 to time k constitute the estimation of the target state sequence X∧k=X^0,X^1,⋯X^k.


## 3.3. 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 asNL/D,Vs) of the maneuver target are in a continuous range ωmin,ωmax that the model set parameters at time k are ωkL,ωkc,ωkR, and that the posterior probabilities corresponding to the three model set parameters at time k are μkL,μkc,μkR; 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.(38)ωk+1c=ωkLμkL+ωkCμkC+ωkRμkR.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 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 thresholdtg. When the maximum posterior probability of the model exceeds the threshold, adjust the parameters of the central model ωk+1C(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(39)ωk+1C=αωkC+1‐αμkLωkL+μkCωkC+μkRωkR.(3)
IfμkC=maxμkL,μkC,μkR, the intermediate model has the highest probability
(40)ωk+1L=ωk+1C−λkL/2,μkL<t1,ωk+1C−λkL,others,ωk+1R=ωk+1C+λkR/2,μkL<t1,ωk+1C+λkR,others,where t1 is the posterior probability threshold. If the posterior probability of the models on the left and right sides is less than t1, the model parameters do not match the current maneuver mode and thus reduce the interval from the grid center.(4)
IfμkL=maxμkL,μkC,μkR, then the right model keeps the original interval unchanged and the left model adjusts the model interval according to the threshold value
(41)ωk+1L=ωk+1C−λkL2,μkL<t2,ωk+1C−λkL,others,ωk+1R=ωk+1C+λkR.(5)
IfμkR=maxμkL,μkC,μkR, then the left model keeps the original interval unchanged and the right model is adjusted according to the threshold value
(42)ωk+1L=ωk+1C−λkL,ωk+1R=ωk+1C+λkR2,μkL<t3,ωk+1C+λkR,others.The above rules for parameter updates in the maneuver model constitute the correction law of model probability ingωk−1i,pMk−1jωk−1jZkj=1,2,⋯,N.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 timek 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 modeli transferring to model j at time k be πkij and satisfy ∑j=1Mπkij=1. The elements of the Markov probability transition matrix are modified by the following methods:
(43)πk,0ij=λkijπk−1ij;λkij=ΛkjΛkiξ,0≤ξ≤1,where ξ=0.5 is used to control the speed of adjustment. The corresponding probability transition matrix is transformed into the following:
(44)πk,1ij=πk,0ij∑j=1Mπk,0ij=λkijπk−1ij∑j=1Mλkijπk−1ij.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 thresholdTh=0.5 for the main diagonal element if πk,1ij<Th,0<Th<1, then πk,1ij=Th, which thus has the following form:
(45)πkij=πk,1ij−Th−πk,1ijπk,1ij∑j=1Mπk,1ij−πk,1ij.


## 4. 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.
(46)H=12aTatf−tn+λrT⋅vM+λvT⋅grM+a,where λr and λV are the costate vectors corresponding to the relative position vector and the relative velocity vector, respectively. By solving the optimal control problem using the variational method, the control equation and adjoint equation can be obtained as follows:
(47)∂H∂a=atf−tn+λv=0∂H∂rM=−λ̇r=0∂H∂vM=−λ̇v=λr⟹λr=λrf,λv=λvf+tf−t⋅λrf,a=−tf−tn⋅λv.According to the guidance Equation (14) and the optimal control Equation (47), the analytical form of the costate vector and the current flight state based on the control variables and the terminal position and velocity of the interceptor can be obtained:
(48)rf=Δrg−tf−tn+3n+3λrf−tf−tn+2n+2λvf+tf−t⋅vMt+rMt,vf=Δvg−tf−tn+2n+2λrf−tf−tn+1n+1λvf+vMt,where Δvg=∫ttfgMrMτdτ is the velocity increment caused by gravity and Δrg=∫ttf∫tsgrMτdτds is the position increment caused by gravity. By solving the state constraint Equation (48), the Lagrange multipliers satisfying the midcourse guidance model of the interceptor Equation (20) can be obtained:
(49)λrf=−n+22n+3tf−tn+3⋅rZEM+n+1n+2n+3tf−tn+2⋅vZEM,λvf=n+1n+2n+3tf−tn+2⋅rZEM−n+1n+22tf−tn+1⋅vZEM,where rZEM=rf−rMt+tf−t⋅vMt+Δrg and vZEM=vf−tf−t⋅vMt+Δvg. The physical meaning of ZEM is the position deviations caused by the interceptor flight to the terminal time under zero overload conditions; similarly, the physical meaning of ZEV is the position deviations under the same conditions. After substituting Equations (48) and (49) into Equation (41), the ZEM/ZEV feedback guidance law under the performance index Equation (15) can be obtained as
(50)a=n+2n+3tf−tn+2⋅rZEM−n+1n+2tf−t⋅vZEV.This method takes the minimum full-time overload with a time-varying weight as the performance indicator. While effectively reducing full-time overload, it ensures that the algorithm has an overload convergence, which is perfect from a theoretical point of view. Its core is to modify the predicted deviation in relative position and relative velocity based on the idea of correcting predictions to adapt to the terminal constraint equations(51)ψirMtf,vMtf=0;i=1,2,3,4,5.To effectively reduce overload in the whole process and to consider the convergence of overload, we add time-varying weights to the performance indicators.


### 4.2. Midcourse Guidance Scheme Combined with Trajectory Prediction

Based on an estimate of the target state and adjustment of the parameters of the set of models in the first two sections, the future target prediction state is obtained using numerical integration of Equation (12) to obtain the NSGV target prediction-trajectory envelope of the NSGVs. The terminal constraint of the midcourse guidance method is obtaining a matched predicted point of interception in the target predicted trajectory envelope Sset. The calculation flow for the predicted intercept point based on predictions of the trajectory is shown in Figure 4. The detailed calculation steps are as follows:
(1)
Design maneuver model setMk: given the minimum, middle, and maximum values of each component of model parameter vector Equation (23), a total of N feasible maneuver modes are formed(2)
Calculate target state estimationX^k and model probability μkiωki: input radar measurement data sequence Zk and perform the evaluation of target state X^k and model probability μkiωki using IMM-ACKF, Equations (37) and (35)(3)
Online adjustment of model parametersωki: through the model probability μkiωki obtained in Step (2), combined with the adaptive grid adjustment rules Equations (40)–(42), the maneuver model parameters can be adjusted online(4)
Calculation of target prediction trajectory envelopeSset: calculate the target prediction trajectory envelopes using numerical integration Equations (8)–(12) for each maneuver model Mkiωki and evaluate the accuracy of the prediction results in combination with the model probability μkiωki obtained in Step (2)(5)
CalculationX^tp of the predicted intercept point: according to the trajectory envelope Sset calculated in Step (4), select the maximum probability of the model μkiωki and the flight time tp of the interceptor and calculate the predicted intercept point X^Mkiωkitp at time tp using numerical integration of Equation (12)(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 tp and the terminal speed vMtf of the interceptor as the quantity and the maximum terminal speed as the performance index J=−vMtf, the individual optimal position pi,j, and the global optimal position gj are calculated according to the particle swarm optimization algorithm, and the iteration is carried out in the following form
(52)vi,j+1=w⋅vi,j+c1⋅r1×pi,j−εi,j+c2⋅r2×gj−εi,j,εi,j+1=εi,j+vi,j,where r1,2 is the random number used when updating particles; c is the proportion of global optimization and particle optimization; and ε=tp,vMtfT, iterating continuously until the accuracy requirements are met, i.e., gj−gj−1≤Δε, and Equation (50) outputs the guidance instruction a(7)
Calculate the start time of the boost phasetig: estimate the remaining flight time tgotig which satisfies vTt+tgotig−vMt+tgotig·rTt+tgotig+R−rMt+tgotig≥0 and calculate tig which satisfies the remaining flight time constraint, i.e., Δtgotig=tf−t−tgotig=0. Newton-Raphson iterative method is an effective approach used to solve the problem. The iterative variable is formulated as follows:(53)tig,j=tig,j−1−Δtgo,jtigΔṫgo,j−1tig,j=1,2⋯.Figure 4
Logic diagram of midcourse guidance considering predictions of the trajectory for maneuver NSGVs.Iterating continuously until the accuracy requirements are met, i.e.,tig,j−tig,j−1≤Δε.


## 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.
(46)H=12aTatf−tn+λrT⋅vM+λvT⋅grM+a,where λr and λV are the costate vectors corresponding to the relative position vector and the relative velocity vector, respectively. By solving the optimal control problem using the variational method, the control equation and adjoint equation can be obtained as follows:
(47)∂H∂a=atf−tn+λv=0∂H∂rM=−λ̇r=0∂H∂vM=−λ̇v=λr⟹λr=λrf,λv=λvf+tf−t⋅λrf,a=−tf−tn⋅λv.According to the guidance Equation (14) and the optimal control Equation (47), the analytical form of the costate vector and the current flight state based on the control variables and the terminal position and velocity of the interceptor can be obtained:
(48)rf=Δrg−tf−tn+3n+3λrf−tf−tn+2n+2λvf+tf−t⋅vMt+rMt,vf=Δvg−tf−tn+2n+2λrf−tf−tn+1n+1λvf+vMt,where Δvg=∫ttfgMrMτdτ is the velocity increment caused by gravity and Δrg=∫ttf∫tsgrMτdτds is the position increment caused by gravity. By solving the state constraint Equation (48), the Lagrange multipliers satisfying the midcourse guidance model of the interceptor Equation (20) can be obtained:
(49)λrf=−n+22n+3tf−tn+3⋅rZEM+n+1n+2n+3tf−tn+2⋅vZEM,λvf=n+1n+2n+3tf−tn+2⋅rZEM−n+1n+22tf−tn+1⋅vZEM,where rZEM=rf−rMt+tf−t⋅vMt+Δrg and vZEM=vf−tf−t⋅vMt+Δvg. The physical meaning of ZEM is the position deviations caused by the interceptor flight to the terminal time under zero overload conditions; similarly, the physical meaning of ZEV is the position deviations under the same conditions. After substituting Equations (48) and (49) into Equation (41), the ZEM/ZEV feedback guidance law under the performance index Equation (15) can be obtained as
(50)a=n+2n+3tf−tn+2⋅rZEM−n+1n+2tf−t⋅vZEV.This method takes the minimum full-time overload with a time-varying weight as the performance indicator. While effectively reducing full-time overload, it ensures that the algorithm has an overload convergence, which is perfect from a theoretical point of view. Its core is to modify the predicted deviation in relative position and relative velocity based on the idea of correcting predictions to adapt to the terminal constraint equations(51)ψirMtf,vMtf=0;i=1,2,3,4,5.To effectively reduce overload in the whole process and to consider the convergence of overload, we add time-varying weights to the performance indicators.


## 4.2. Midcourse Guidance Scheme Combined with Trajectory Prediction

Based on an estimate of the target state and adjustment of the parameters of the set of models in the first two sections, the future target prediction state is obtained using numerical integration of Equation (12) to obtain the NSGV target prediction-trajectory envelope of the NSGVs. The terminal constraint of the midcourse guidance method is obtaining a matched predicted point of interception in the target predicted trajectory envelope Sset. The calculation flow for the predicted intercept point based on predictions of the trajectory is shown in Figure 4. The detailed calculation steps are as follows:
(1)
Design maneuver model setMk: given the minimum, middle, and maximum values of each component of model parameter vector Equation (23), a total of N feasible maneuver modes are formed(2)
Calculate target state estimationX^k and model probability μkiωki: input radar measurement data sequence Zk and perform the evaluation of target state X^k and model probability μkiωki using IMM-ACKF, Equations (37) and (35)(3)
Online adjustment of model parametersωki: through the model probability μkiωki obtained in Step (2), combined with the adaptive grid adjustment rules Equations (40)–(42), the maneuver model parameters can be adjusted online(4)
Calculation of target prediction trajectory envelopeSset: calculate the target prediction trajectory envelopes using numerical integration Equations (8)–(12) for each maneuver model Mkiωki and evaluate the accuracy of the prediction results in combination with the model probability μkiωki obtained in Step (2)(5)
CalculationX^tp of the predicted intercept point: according to the trajectory envelope Sset calculated in Step (4), select the maximum probability of the model μkiωki and the flight time tp of the interceptor and calculate the predicted intercept point X^Mkiωkitp at time tp using numerical integration of Equation (12)(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 tp and the terminal speed vMtf of the interceptor as the quantity and the maximum terminal speed as the performance index J=−vMtf, the individual optimal position pi,j, and the global optimal position gj are calculated according to the particle swarm optimization algorithm, and the iteration is carried out in the following form
(52)vi,j+1=w⋅vi,j+c1⋅r1×pi,j−εi,j+c2⋅r2×gj−εi,j,εi,j+1=εi,j+vi,j,where r1,2 is the random number used when updating particles; c is the proportion of global optimization and particle optimization; and ε=tp,vMtfT, iterating continuously until the accuracy requirements are met, i.e., gj−gj−1≤Δε, and Equation (50) outputs the guidance instruction a(7)
Calculate the start time of the boost phasetig: estimate the remaining flight time tgotig which satisfies vTt+tgotig−vMt+tgotig·rTt+tgotig+R−rMt+tgotig≥0 and calculate tig which satisfies the remaining flight time constraint, i.e., Δtgotig=tf−t−tgotig=0. Newton-Raphson iterative method is an effective approach used to solve the problem. The iterative variable is formulated as follows:(53)tig,j=tig,j−1−Δtgo,jtigΔṫgo,j−1tig,j=1,2⋯.Figure 4
Logic diagram of midcourse guidance considering predictions of the trajectory for maneuver NSGVs.Iterating continuously until the accuracy requirements are met, i.e.,tig,j−tig,j−1≤Δε.


## 5. Results and Discussion

### 5.1. Simulation Results

In Section3, 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 trajectory-prediction method, several cases are applied.


### 5.2. 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), the coefficient of altitude velocity mapping a2 is discretely chosen as [-0.025, 0.0, 0.025]; the lift-drag ratio parameters L/D refer to the variation law of ballistic characteristics, which is taken as 1.82.83.8 in proper order; and the relevant parameter of the roll angle inversion profile VS in Equation (24) is chosen as [800 m/s, 1300 m/s, and 1800 m/s] in turn. These three parameters can be cross-combined to form the corresponding maneuver model set, which is divided according to the maneuver model parameters. The model set is composed of 27 maneuver models based on trajectory tracking. At the same time, the initial probability of each maneuver model μ0j is set to 1/27. The threshold for the model probability tg corresponding to the change in the grid center is set to 1/25, and the probability thresholds of t1, t2 and t3 are set to 1/25. The probability of model i transferring to model j (j≠i) at time 0 π0ij is set to 0.3%, and the probability of model i remain unchanged at time 0 π0ii is set to 92.2%.The four nominal observation parameters areRrλDλTṘrT from Equation (25); the mean square deviation (1σ) of the observation noise for the relative distance nR is 100m/s; the mean square deviation (1σ) of the observation noise for the relative distance change rate nṘ is 10m/s; and the mean square deviation (1σ) of the angle of observation noise nD,nT 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.Figure 5
Comparison of the trajectory-prediction results between the proposed method and fixed maneuver model. (a) Two-dimensional curve comparison. (b) Comparison of spatial curves.
(a)(b)The simulation results in Figure5 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
The amount of deviation in the prediction of the trajectory corresponding to different time lengths. (a) No observation and tracking. (b) Track and identify for 5 s. (c) Track and identify for 15 s. (d) Track and identify for 25 s.
(a)(b)(c)(d)Figure6 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.Figure6(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.Table1 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.Table 1
Accuracy analysis of trajectory prediction method (the terminal time is 75 s).
Start prediction time (s)35405060Duration of observation (s)051525Duration of prediction (s)15151515Mean deviation of all model (km)1.92251.20591.19761.1209Maximum probability model’s deviation (km)——1.29241.22431.0124


### 5.3. 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 seekerRmax 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 Figure7. 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(a) and (b) depict normal overload and lateral overload during terminal guidance phase. Combine with the results in 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 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.Figure 7
Changes in the flight parameters for interception using different midcourse guidance methods under multiple constraints based on predictions of the trajectory. (a) Normal overload change curve. (b) Lateral overload change curve. (c) Curve of the track angle curve. (d) Curve of the track declination angle. (e) Relative distance curve. (f) Velocity curve.
(a)(b)(c)(d)(e)(f)Figure 8
Comparison of results considering predictions of trajectories of NSGVs under PN and the proposed guidance methods. (a) Trajectories of interception using different midcourse guidance methods. (b) Trajectories of interception during terminal guidance phase using different midcourse guidance methods. (c) Normal overload change curve during terminal guidance phase. (d) Lateral overload change curve during terminal guidance phase.
(a)(b)(c)(d)Table 2
Simulated comparison of the results for an interception of the maneuver target under multiple guidance methods.
Guidance methodsψ1(-)ψ2(-)ψ3(°)ψ4(m/s)PN with TP-0.05050.132620.032550.0214PN without TP-0.0121-0.0327-12.671119.0337ZEM/ZEV (n=0)-0.0003-9.286e-50.0110.2371ZEM/ZEV (n=1)0.00028.982e-50.0050.3249ZEM/ZEV (n=2)0.00071.7e−4-0.0020.2495The simulation result in Table2, 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 ψrMtf,vMtf=0 established by the terminal velocity direction and the terminal position vector, to guide the interceptor to accurately collide with the target NSGVs.


## 5.1. Simulation Results

In Section3, 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 trajectory-prediction method, several cases are applied.


## 5.2. 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), the coefficient of altitude velocity mapping a2 is discretely chosen as [-0.025, 0.0, 0.025]; the lift-drag ratio parameters L/D refer to the variation law of ballistic characteristics, which is taken as 1.82.83.8 in proper order; and the relevant parameter of the roll angle inversion profile VS in Equation (24) is chosen as [800 m/s, 1300 m/s, and 1800 m/s] in turn. These three parameters can be cross-combined to form the corresponding maneuver model set, which is divided according to the maneuver model parameters. The model set is composed of 27 maneuver models based on trajectory tracking. At the same time, the initial probability of each maneuver model μ0j is set to 1/27. The threshold for the model probability tg corresponding to the change in the grid center is set to 1/25, and the probability thresholds of t1, t2 and t3 are set to 1/25. The probability of model i transferring to model j (j≠i) at time 0 π0ij is set to 0.3%, and the probability of model i remain unchanged at time 0 π0ii is set to 92.2%.The four nominal observation parameters areRrλDλTṘrT from Equation (25); the mean square deviation (1σ) of the observation noise for the relative distance nR is 100m/s; the mean square deviation (1σ) of the observation noise for the relative distance change rate nṘ is 10m/s; and the mean square deviation (1σ) of the angle of observation noise nD,nT 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.Figure 5
Comparison of the trajectory-prediction results between the proposed method and fixed maneuver model. (a) Two-dimensional curve comparison. (b) Comparison of spatial curves.
(a)(b)The simulation results in Figure5 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
The amount of deviation in the prediction of the trajectory corresponding to different time lengths. (a) No observation and tracking. (b) Track and identify for 5 s. (c) Track and identify for 15 s. (d) Track and identify for 25 s.
(a)(b)(c)(d)Figure6 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.Figure6(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.Table1 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.Table 1
Accuracy analysis of trajectory prediction method (the terminal time is 75 s).
Start prediction time (s)35405060Duration of observation (s)051525Duration of prediction (s)15151515Mean deviation of all model (km)1.92251.20591.19761.1209Maximum probability model’s deviation (km)——1.29241.22431.0124


## 5.3. 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 seekerRmax 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 Figure7. 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(a) and (b) depict normal overload and lateral overload during terminal guidance phase. Combine with the results in 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 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.Figure 7
Changes in the flight parameters for interception using different midcourse guidance methods under multiple constraints based on predictions of the trajectory. (a) Normal overload change curve. (b) Lateral overload change curve. (c) Curve of the track angle curve. (d) Curve of the track declination angle. (e) Relative distance curve. (f) Velocity curve.
(a)(b)(c)(d)(e)(f)Figure 8
Comparison of results considering predictions of trajectories of NSGVs under PN and the proposed guidance methods. (a) Trajectories of interception using different midcourse guidance methods. (b) Trajectories of interception during terminal guidance phase using different midcourse guidance methods. (c) Normal overload change curve during terminal guidance phase. (d) Lateral overload change curve during terminal guidance phase.
(a)(b)(c)(d)Table 2
Simulated comparison of the results for an interception of the maneuver target under multiple guidance methods.
Guidance methodsψ1(-)ψ2(-)ψ3(°)ψ4(m/s)PN with TP-0.05050.132620.032550.0214PN without TP-0.0121-0.0327-12.671119.0337ZEM/ZEV (n=0)-0.0003-9.286e-50.0110.2371ZEM/ZEV (n=1)0.00028.982e-50.0050.3249ZEM/ZEV (n=2)0.00071.7e−4-0.0020.2495The simulation result in Table2, 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 ψrMtf,vMtf=0 established by the terminal velocity direction and the terminal position vector, to guide the interceptor to accurately collide with the target NSGVs.


## 6. 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.



---

*Source: 4528803-2022-12-28.xml*