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.
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 variables
Based on the above assumptions, the three-dimensional point-mass equations of motion in near space can be expressed as follows: where and are state vectors; is the unit vector defining the longitudinal direction of the rocket body; is the normal direction of the rocket body; is the lateral direction of the rocket body; is a function of gravitational acceleration with as the independent variable; , , and are axial aerodynamic force, normal aerodynamic force, and lateral aerodynamic force, respectively; is the thrust magnitude; and is the mass of the aerocraft. When the subscript in Equation (1) is taken as , it represents the interceptor, and when it is taken as , it represents NSGV. The thrust magnitude is given by the following: where is the engine exhaust velocity; is the normalized mass burn rate of the interceptor; is the cross-sectional area of the engine nozzle; is a signum function whose form will be given later; is a function of atmospheric pressure with flight altitude of the interceptor as the independent variable. The aerocraft mass follows the following differential equation:
The axial, normal, and lateral aerodynamic forces , , and are given as follows: where is the aerodynamic reference area; is the dynamic pressure; , , and are the coefficient of axial force, the coefficient of normal force, and the coefficient of lateral force, respectively. , , and can be expressed as the pitch angle , the yaw angle, and the roll angle in inertial reference frame as follows: where , , and 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, is a signum function representing the working state of the engine, which is expressed as follows:
In the boost phase, thrust and aerodynamic force together produce a large acceleration to change the flight trajectory; in the glide phase, the available acceleration generated by an aerodynamic force is limited by the flight conditions [2], which can be expressed as follows: where ; , , and 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. Let represents the target measurement sequence from time 0 (corresponding to) to time (corresponding to ) obtained by the ground-based radar; let represents the estimation of target state at the time ; let represents the target dynamic equation corresponding to (1), in which represents the target control function to be inferred. Thus far, the expected output of the target trajectory prediction problem is the target state prediction from to a certain time in the future, determined by , or a target state prediction trajectory during a period.
Considering the uncertainty of measurement information and target maneuver , the common practice in predicting the trajectory is to infer the density of the target control function through known measurement information and then to calculate the target state prediction trajectories through one or more high-probability control functions .
The trajectory-prediction method constructs the target maneuver model set according to the maneuver mode of glider (NSGVs) takes the data sequence measured by radar as input , estimates the target state , and infer the model probability using the control quantity determined by each model and then modifies the parameters of the maneuver model according to the model probability. Finally, the predicted flight state of each model at time is calculated through (8), in which the one with the maximum model probability used as the predicted intercept point for midcourse guidance. The detailed design process of the trajectory prediction problem is as follows.
To infer the maneuver law of the target, the maneuver model set is designed to express different maneuver modes, whose element represents the th model at time and is the maneuver parameter that needs to be inferred. Then, the expression of the maneuver mode set at time is as follows: where represents the total number of maneuver models used. Since the radar cannot directly obtain the state of NSGVs through measurement, it is necessary to estimate the target state according to the target control function determined by the maneuver model . According to different parameters of maneuver model , N different target state estimation, and the model probability can be obtained. Therefore, the estimation of target state can be calculated as
To estimate and predict the target flight trajectory more accurately, the highest probability model parameter is used to gradually modify other low-probability model parameters so that the trajectory cluster calculated by the maneuver model set gradually becomes closer to the observed trajectory. The updated function of the maneuver parameters expressed by the model probability can be expressed as follows:
The control function calculated with the maneuver model parameter is integrated according to the dynamic differential equation in (8) to obtain a prediction of the target trajectory. The target state prediction of during interception is expressed as follows: where represents the time when the prediction starts; that is, the time corresponding to the current time . A series of target-predicted trajectory envelopes from to under each maneuver mode can be expressed as follows:
In summary, the NSGV trajectory-prediction process includes the design of a maneuver model set, the evaluation of the target state estimation and model probability, the adjustment of the maneuver model parameters, and the calculation of the target prediction-trajectory envelope.
2.3. Midcourse Guidance Model
To unify and simplify the derivation of the guidance equation, the specific force acceleration is introduced to replace Equation (1), which can be compactly expressed as follows:
The intermediate guidance problem [22] of the interceptor is used to solve the optimal control with constraints, and the performance index is usually designed to minimize the total overload. This paper adopts a performance index with a weight varying over time, which can be expressed as follows: where 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 is equal to the target trajectory-prediction time ; and 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 . Generally, when , 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 . In order to make the terminal guidance more conducive to an accurate collision with the target aircraft, the terminal condition of the midcourse guidance, shown in Figure 1, must meet the collision triangle relationship and intercept the target at certain intersection angles and : which have significant impacts on the distance of a miss in the terminal guidance. Then, the velocity vector of the interceptor at the terminal time should meet the following:
According to the geometric relationship of the collision triangle [26], the expression of the position vector at the terminal time under the condition of zero control interception is as follows:
When the missile target distance reaches the maximum detection distance of the interceptor seeker, the interceptor is considered to have completed the shift handover for intermediate and terminal guidance and directly enters the terminal guidance phase, so the shift handover distance constraint is as follows: where and are obtained from the trajectory-prediction method in Equation (12), and . As a result, the midcourse guidance model of the interceptor can be expressed as follows:
3. Trajectory Prediction by IMM-ACKF
According to the establishment of the maneuver mode set in Section 2.3, the IMM method is used to evaluate the model probability of each maneuver mode (see Sections 3.1 and 3.2), and then, the target state is estimated on the basis of 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 set provides a series of feasible target control functions for tracking the target state . 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 , and the acceleration generated by each component composed of drag, lift, and laterals force can be sorted as follows: where represents the remaining range; represents the current speed; is the expected terminal speed; is the current intersection angle; is the function of the correlation coefficient of the altitude velocity mapping relationship [29]; is the gravitational acceleration; is the lift-drag ratio; and represents the lateral boundary turnover function, which is defined as follows: where is the lateral flight line of a sight error, represents the symbolic calculation result of the last time, and is the boundary threshold related to the lateral boundary parameter . In the design of the abovementioned maneuver models, the model parameters that can be selected include the correlation coefficient of the altitude–velocity mapping relationship, the lift-drag ratio parameter , and the parameter expressing the lateral boundary. The design parameter vector of the th maneuver mode can be selected as follows:
For the design of the aircraft maneuver mode set at the initial time (), 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 feasible maneuver modes.
When intercepting NSGVs, a high-power active-passive composite radar is generally used to monitor and track the target remotely. At this time, the radar can measure the distancefrom the target, the target azimuth , and the high-low angle and can obtain the change rate of the target distance from the Doppler principle. As shown in Figure 2, the radar measurement equation of NSGVs is given. where is the position vector of the target relative to the radar and is the velocity vector of the target relative to the radar.
Some errors and noises are found between the actual radar observation and the real state of the target. Equations (22) and (25) can be expressed as follows: where is the simulated radar observation noise and the matrix parameters are closely related to the radar performance. Radar measurement error can usually be considered Gaussian white noise with zero mean. where represents a random variable that obeys zero mean and covariance . Generally, the observation noise of each channel is believed to be independent of each other, and the radar observation noise covariance matrix can be obtained as follows:
The diagonal elements represent the covariance in the measurement error of each channel. The establishment of a dynamic equation is made on the basis of maneuver model modelling, and the establishment of a radar measurement coordinate system is made on the basis of a radar pseudo-observation measurement and a filtering equation. The result of the measurement in Equation (25) from time 0 to time constitutes the target measurement information sequence .
3.2. Maneuver Target Tracking Method
According to the multiple maneuver models in Equation (24) from Section 3.1 and the measurement information sequence of the target, the IMM method was used to evaluate the model probability of each maneuver mode, and the probability of each maneuver model was obtained. Because the maneuver mode of the target cannot be well determined in advance, the single-motion model cannot accurately describe the flight state of the target. The interactive multimodel algorithm uses the combination of multiple motion models to estimate the system state through effective weighted fusion. The specific principle of the combination of IMM [9] and adaptive culture Kalman filter (ACKF) is shown in Figure 3.
The maneuver model set covers all possible typical maneuver modes; the deviation from the real motion model of the target is large. To prevent divergence in the filtering process, the fading factor is introduced into the state prediction covariance matrix, the gain matrix is adjusted online and in real time, and the output residual sequence is forced to be orthogonal. The calculation process of fading factor is shown in Equation (29). where are the variances in system noise and measurement noise, respectively; are the one-step prediction error covariance matrix without and with system process noise in CKF; . is the weight of the volume point; are the volume point, state estimation, and measurement residual generated by one-step state prediction without considering fading factor; is the values of predicted measurement obtained by weighting ; is the measurement functions, as shown in Formula (25); is the weakening factor, which is generally selected according to experience; and is the forgetting factor. In the formula, is the calculation process parameter matrix. The covariance matrix is modified by the fading factor to obtain the one-step prediction state covariance matrix:
After introducing the fading factor, conforms to the distribution . The predicted state obtained by one-step prediction and the one-step predicted state covariance matrix 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 is calculated based on the Gaussian hypothesis. The likelihood function value of each model is as follows: where , , and are the observation at time , is the observation prediction value before introducing the fading factor, and is the observation prediction covariance matrix for ACKF before calculating the fading factor.
The probability-estimation process in IMM-ACKF is described as below.
3.2.1. Interaction/Mixing
Mixed probability of model to model is calculated as follows: where is calculated as:
The mixed state and mixed covariance of the corresponding model are estimated as follows:
3.2.2. Mode-Matched Filtering Through ACKF
The mixed state and mixed covariance (34) are used as input to the ACKF matched to , which uses to yield . The model probability is updated as follows: where is calculated as
3.2.3. State Estimate and Covariance Combination
State estimate and covariance combination are done as follows:
The state and covariance estimation describe the form of a posteriori probability density function for the target state under a Gaussian assumption. The above state estimation from time 0 to time constitute the estimation of the target state sequence .
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 as ) of the maneuver target are in a continuous range that the model set parameters at time are , and that the posterior probabilities corresponding to the three model set parameters at time are ; the model set of parameters at the time is determined as follows.
An adjustment in the model set center parameters is obtained by taking the model posterior probability as the weight coefficient.
The inaccuracy of the system model, the influence of observation errors, and the competition between models cause the maximum a posteriori probability of the model to not always appear in the model closest to the real model. Therefore, continuous jumps, shaking, and even divergence of the central model around the real model space are easily produced. A grid center adjustment is thus needed; (1)To prevent irregular jitter of the model set, we set the grid center change threshold . When the maximum posterior probability of the model exceeds the threshold, adjust the parameters of the central model (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 (3)If , the intermediate model has the highest probability where is the posterior probability threshold. If the posterior probability of the models on the left and right sides is less than , the model parameters do not match the current maneuver mode and thus reduce the interval from the grid center.(4)If , then the right model keeps the original interval unchanged and the left model adjusts the model interval according to the threshold value (5)If , then the left model keeps the original interval unchanged and the right model is adjusted according to the threshold value
The above rules for parameter updates in the maneuver model constitute the correction law of model probability in .
The value of the likelihood function of the degree of the model reflects the match between the model and the real motion pattern of the target. The value of the likelihood function at time for model is calculated according to (31). A larger value relative to the likelihood function value of other models indicates that the model better matches the real motion state of the target and thus has a greater probability of other models transitioning to this matching model.
Let the probability of model i transferring to model at time be and satisfy . The elements of the Markov probability transition matrix are modified by the following methods: where is used to control the speed of adjustment. The corresponding probability transition matrix is transformed into the following:
These methods may cause some values of the main diagonal elements of the matrix to become smaller and smaller; the probability of transferring the mismatched model to itself becomes smaller and smaller. However, when the target maneuvers to turn the mismatched model into a matching model, the model switching is slow. Therefore, we set a threshold for the main diagonal element if , then , which thus has the following form:
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. where and 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:
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: where is the velocity increment caused by gravity and 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: where and . 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
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
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 . 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 set : given the minimum, middle, and maximum values of each component of model parameter vector Equation (23), a total of feasible maneuver modes are formed(2)Calculate target state estimation and model probability : input radar measurement data sequence and perform the evaluation of target state and model probability using IMM-ACKF, Equations (37) and (35)(3)Online adjustment of model parameters : through the model probability 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 envelope : calculate the target prediction trajectory envelopes using numerical integration Equations (8)–(12) for each maneuver model and evaluate the accuracy of the prediction results in combination with the model probability obtained in Step (2)(5)Calculation of the predicted intercept point: according to the trajectory envelope calculated in Step (4), select the maximum probability of the model and the flight time of the interceptor and calculate the predicted intercept point at time 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 and the terminal speed of the interceptor as the quantity and the maximum terminal speed as the performance index , the individual optimal position , and the global optimal position are calculated according to the particle swarm optimization algorithm, and the iteration is carried out in the following form where is the random number used when updating particles; is the proportion of global optimization and particle optimization; and , iterating continuously until the accuracy requirements are met, i.e., , and Equation (50) outputs the guidance instruction (7)Calculate the start time of the boost phase : estimate the remaining flight time which satisfies and calculate which satisfies the remaining flight time constraint, i.e., . Newton-Raphson iterative method is an effective approach used to solve the problem. The iterative variable is formulated as follows:
Iterating continuously until the accuracy requirements are met, i.e., .
5. Results and Discussion
5.1. Simulation Results
In Section 3, the NSGV maneuver model set was established based on target intent, the radar observation model was used to track the trajectory of NSGVs, and the interactive multimodel variable structure filtering method predicted feasible target trajectories. Then, in Section 4, the midcourse guidance method satisfying the terminal constraint of intersection angle guided the interceptor to collide with the target. In this section, to evaluate the performance of the proposed guidance and 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 is discretely chosen as [-0.025, 0.0, 0.025]; the lift-drag ratio parameters refer to the variation law of ballistic characteristics, which is taken as in proper order; and the relevant parameter of the roll angle inversion profile 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 is set to . The threshold for the model probability corresponding to the change in the grid center is set to , and the probability thresholds of , and are set to . The probability of model transferring to model () at time 0 is set to , and the probability of model remain unchanged at time 0 is set to .
The four nominal observation parameters are from Equation (25); the mean square deviation () of the observation noise for the relative distance is ; the mean square deviation () of the observation noise for the relative distance change rate is ; and the mean square deviation () of the angle of observation noise of X and Y is 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.
(a)
(b)
The simulation results in Figure 5 show obvious deviations in trajectory prediction in the longitudinal plane caused by the fixed maneuver model. For comparison, the adaptive grid in Equation (38) updates the parameters of the maneuver model set according to the measurement of the target trajectory by the radar to bring the predicted trajectories integrated by the maneuver model set in Equation (24) close to the true trajectory. The deviation in predictions of the trajectory at different times is plotted in Figure 6.
(a)
(b)
(c)
(d)
Figure 6 illustrates the trajectory prediction deviation obtained by the proposed IMM-ACKF method under the observation conditions of 0 s, 5 s, 15 s, and 25 s. The blue line in the figure is the predicted trajectory deviation corresponding to the maneuver model with high probability, and the corresponding red line is the trajectory-prediction deviation corresponding to the maneuver model with low probability. Obviously, with the increase in observation time and the decrease in prediction time for the response, the amount of deviation in the prediction increases first and then decreases, which shows a certain convergence because of the flight intention and constraints of the terminal mission for NSGVs.
Figure 6(a) shows that the target trajectory of the NSGVs is not observed and that the prediction parameters are not updated, which is consistent with the parameters of the fixed maneuver model in Figure 5. Additionally, the fixed maneuver model set fails to cover the true trajectory, so the prediction deviation is larger than the adjusted maneuver model updated by IMM-ACKF trajectory prediction. In Figure 6(b), the trajectory predicted by some high-priority maneuver models has a large deviation compared with the low-priority models primarily because the real target overturned the pitch angle and the model set was not updated in time. Therefore, the set of trajectories predicted by the high-probability models still needs to be updated and adjusted many times according to the actual radar measurement data to reduce the amount of deviation from the prediction of the trajectory.
Table 1 shows the simulation results of trajectory prediction accuracy deviation calculated by the proposed IMM-ACKF method in the same prediction duration of 15 seconds under different observation intervals and different initial prediction times. Table 1 shows that the trajectory prediction error of the proposed method decreases with the increase of the observation duration. And the predicted trajectory generated by the maximum probability model is close to the average value predicted by all other models in Equation (24). Therefore, the predicted trajectory generated by the maximum probability model can be used as the predicted value of the terminal constraint of the midcourse guidance method.
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 seeker was 20 km, and the constraints of the track angles and were and in the reverse flight conditions, respectively. At the same time, the methods proportional navigation with the trajectory prediction method (PN with TP) and proportional navigation without trajectory prediction (PN without TP) were introduced for a comparison with the methods proposed in this paper.
The simulation curve of the interception flight parameters based on the prediction information of the target trajectory is shown in Figure 7. In Figure 8, the change in the position of the interceptor when considering predictions of the trajectory of NSGVs is shown in (a), and the trajectories of interception during terminal guidance phase using different midcourse guidance methods is shown in (b). Figure 8(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 .
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
The simulation result in Table 2, compared with PN using the TP method and PN without using the TP method, the deviations in the third and fourth constraints are too large to meet the expected accuracy because the intersection angle is not considered in the design of these two methods. The proposed improved ZEM/ZEV method satisfies multiple constraints composed of the control interception manifold and intersection track angle, including the five equations of 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.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this article.
Acknowledgments
This research was funded by the Natural Science Basic Research Program of Shaanxi Province (Program No. 2022JQ-061) and the China Postdoctoral Science Foundation (Program No. 2022M712588).