Adaptive Markov IMM Based Multiple Fading Factors Strong Tracking CKF for Maneuvering Hypersonic-Target Tracking

Abstract

Hypersonic targets have complex motion states and high maneuverability. The traditional interactive multi-model (IMM) has low tracking accuracy and a slow convergence speed. Therefore, this paper proposes a strong tracking cubature Kalman filter (CKF) adaptive interactive multi-model (AIMM) based on multiple fading factors. Firstly, this paper analyzes the structure of the CKF algorithm, introduces the fading factor of the strong tracking algorithm into the covariance matrix of the time update and measurement update, and adjusts the filter gain online and in real time, which can reduce the decline infilter accuracy caused by model mismatch. Secondly, Singer model, “current” statistical (CS) model, and Jerk model are selected in the model set of IMM and introduced singular value decomposition (SVD) decomposition to solve the problem that Cholesky decomposition cannot be performed in the CKF due to the model dimension expansion. Last, an adaptive algorithm for the Markov matrix in the IMM is proposed. The transition probability was adaptively modified by the value of the model likelihood function to enhance the proportion of matching models. The simulation results show that the proposed algorithm enhanced the proportion of matching models in the IMM and improved the tracking accuracy by 16.51% and the convergence speed by 37.5%.

1. Introduction

A hypersonic glide target (HGT) [1] has a flight Mach number greater than 5. It has the characteristics of a fast flight speed, strong penetration capability, strong mobility, and high accuracy. It adopts an aerodynamic shape with a large lift-drag ratio, and its trajectory is complex and changeable, so it is difficult to achieve effective defense. Therefore, tracking and locating the target have always been difficult and hot issues in the field of target tracking.

With the continuous deepening of the research on the characteristics of maneuvering targets, experts and scholars have proposed a variety of single-model motion models to fit maneuvering targets: A. W. B. et al. proposed the constant velocity (CV) model and constant acceleration (CA) model in 1973 [2]. The Singer model, also known as the time-dependent model, was first developed by R Singer and proposed in 1970 [3]. Professor Hongren Zhou of China put forward the CS model in the 1980s [4]. In 1997, Indian scholar Kishore proposed the Jerk model [5], with reference to the Singer model, expanded the dimension of traditional acceleration, and proposed the concept of acceleration. Some experts and scholars have improved these models alone, but for highly maneuvering targets such as hypersonic targets, it is difficult for a single model to match their motion states. When the difference between motion states is large, it will lead to filtering divergence and tracking failure [6].

To solve the above problems, in 1965, Magill proposed multiple models (MM) [7] to improve the limitation of single models. On the basis of the multi-model algorithm, Bar Shalom and Blom proposed the IMM algorithm by using the pseudo-Bayesian algorithm, and used the Markov transformation matrix to overcome the filtering divergence caused by the fixed model [8]. It tracks the uniform motion process, significantly reduces the noise, and has high tracking accuracy and a fast convergence speed [9]. Since the IMM algorithm was proposed, scholars have applied it to the target tracking field and conducted a significant amount of research. For the tracking of hypersonic gliding targets, Ref. [10] proposed a self-adaptive algorithm of maneuvering frequency in the Singer model, which uses filtering information to calculate and adjust the maneuvering frequency in the Singer model. Compared with the IMM Singer algorithm, it has high tracking accuracy and requires a small amount of calculation. The nonlinear filtering algorithm plays an important role in the field of target tracking, such as multi-base radar positioning [11] and nonlinear filtering vehicle adaptive navigation [12]. Since the measurement model of hypersonic target radar is also nonlinear, Ref. [13] proposed the use of the unscented Kalman filter (UKF) [14] combined with the IMM algorithm to realize the tracking filtering of hypersonic targets. For the problem of model interaction, Ref. [15] proposed an improved AIMM algorithm with Markov parameters. The tracking accuracy is improved, but it is only applicable to slow maneuvering targets with prior information. On this basis, Ref. [16] proposed an adaptive variable structure interactive multi-model (AVSIMM), which improves the tracking accuracy of hypersonic targets.

To sum up, there are three main problems in the IMM tracking algorithm for highly maneuvering targets such as hypersonic targets: The first is, the selection of models. The conventional CV, CA, and other models are weak in mobility, and cannot match the motion state of hypersonic targets. The second is the filtering algorithm. Most of the measurements of hypersonic targets are nonlinear measurement equations, which are not suitable for the linear Kalman filter (KF) in the IMM model; at the same time, the accuracy of the filtering algorithm needs to be improved. The third is the model interaction problem. The Markov matrix of the conventional IMM algorithm is fixed, the model probability interaction adaptability is weak, and the matching accuracy is low.

On the basis of previous algorithms, this paper analyzed the structure of the CKF [17] and proposes a strong tracking CKF algorithm with multiple fading factors to improve the filtering accuracy. To solve the problem that Cholesky decomposition fails due to the different model dimensions in the CKF algorithm, an improved strong tracking CKF algorithm based on SVD [18] is proposed. The aim of the algorithm is to improve the standard Markov matrix to improve the adaptability of the algorithm. The simulation results show that, compared with the single-model algorithm and the standard IMM algorithm, the proposed algorithm can significantly improve the tracking accuracy and convergence speed of the hypersonic target.

2. Modeling of Hypersonic Target Motion Characteristics

2.1. Motion Equation of Hypersonic Target

In the process of establishing the motion equation of hypersonic targets, the following coordinate systems need to be considered (Figure 1):

(1)

Earth-centered, Earth-fixed (ECEF) coordinate system: The origin $O$ is located at the Earth’s center, $OX$ points to the zero degrees longitude in the equatorial plane, the $OZ$ axis points to the north pole, and the $OY$ axis is perpendicular to the $OX$ and $OZ$ axes in the equatorial plane. The unit vectors on the three $O-XYZ$ axes are expressed as $i$, $j$, and $k$.

(2)

Longitude, latitude, altitude (LLA) coordinate system: This is the geodetic coordinate system. The coordinate origin is located at the geo-center $O$, and the longitude $\mu$, latitude $\lambda$, and altitude $h$ of the target location.

(3)

Aircraft position coordinate system: The origin $O$ is located at the center of the Earth. The $Ox$ axis is the connecting line from the center of the Earth to the hypersonic target, and it is positive upward. The $Oy$ axis is perpendicular to the axis and is in the equatorial plane. The due east direction is positive, and the $Oz$ axis is determined by the other two axes and the right-hand law. The $O-xyz$ coordinate system can be obtained by combining the $O-XYZ$ coordinate system with the latitude and longitude of the hypersonic target position.

(4)

East-north-up (ENU): This is the radar coordinate system. The radar position $R$ is the origin of the coordinate system, and the $RE$, $RN$, and $RU$ axes point to the east, north, and sky directions perpendicular to each other. It can be used to represent the pitch angle $a$ and azimuth angle $e$ of the hypersonic target.

Assuming the Earth is a homogeneous sphere with rotation [19], the aerodynamic force and gravity of the hypersonic target are converted into the same coordinate system through the flight heading angle, aerodynamic angle, attitude angle, and other important angles to obtain the six-degree-of-freedom motion equation of the hypersonic target. Due to the complex structure of the equation, it is simplified when analyzing the scene, and the target trajectory model is obtained as follows [20]:

$$\{\begin{array}{l}\dot{r}=V\sin \gamma \\ \dot{\mu }=\frac{V\cos \gamma \cos \xi }{r\cos \lambda }\\ \dot{\lambda }=\frac{V\cos \gamma \sin \xi }{r}\\ \dot{V}=-\frac{D}{m}-g\sin \gamma +w_e^{2}r\cos \lambda \left(\sin \gamma \cos \lambda -\cos \gamma \sin \lambda \sin \xi \right)\\ \dot{\gamma }=\frac{L\cos \sigma }{mV}+\left(\frac{V}{r}-\frac{g}{V}\right)\cos \gamma +2w_e\cos \lambda \cos \xi \\ +\frac{w_e^{2}r}{V}\cos \lambda \left(\cos \gamma \cos \lambda +\sin \gamma \sin \lambda \sin \xi \right)\\ \dot{\xi }=\frac{L\sin \sigma }{mV\cos \gamma }-\frac{V}{r}\cos \gamma \cos \xi \tan \lambda +\\ 2w_e\left(\tan \gamma \cos \lambda \sin \xi -\sin \lambda \right)-\frac{w_e^{2}r}{V\cos \gamma }\sin \lambda \cos \lambda \cos \xi \end{array}$$

(1)

where $r$ is the distance between the center of mass of the target aircraft and the geo-center, $\mu$ is the geocentric longitude, and $\lambda$ is the geocentric latitude. These three parameters describe the location of the target. $V$ is the target velocity, $\gamma$ is the trajectory inclination angle, and $\xi$ is the trajectory deflection angle. These three parameters describe the target velocity. $D$ and, $L$ are the aerodynamic resistance and aerodynamic lift, $\sigma$ is the inclination angle, $w_e$ is the value of the Earth’s rotation angle velocity, and $g$ is the local gravity acceleration.

In the main movement process of the hypersonic target, the complex aerodynamic force is the main factor to be considered. The common flight aerodynamic lift $L$ and drag $D$ are as follows:

$$\{\begin{array}{l}L=C_{L}qS\\ D=C_{D}qS\\ q=0.5\rho V^2\end{array}$$

(2)

where $C_L$ is the lift coefficient, and $C_D$ is the drag coefficient. These two coefficients are usually related to the flight Mach number and angle of attack of the target. $S$ is the aerodynamic reference area of the hypersonic target, and $\rho$ is the air density at the flight altitude. The atmospheric index model is used to describe this as follows:

$$\rho ={\rho }_0{e}^{-\frac{h}{{\rho }_{\beta }}}$$

(3)

where ${\rho }_0=1.225 \mathrm{k}\mathrm{g}/{\mathrm{m}}^3$ is the standard sea-level atmospheric density, and ${\rho }_{\beta }$ is 6700 m [21].

This paper took Falcon Hypersonic Technology Vehicle 2 (HTV-2) as the research object. Since the aerodynamic parameters of HTV-2 are similar to those of Currency Aerodynamic Vehicle H (CAV-H) [22], the aerodynamic parameters of CAV-H can be cited [23], and the least square method can be used to identify the model parameters:

$$\{\begin{array}{l}{C}_D=C_{D0}+C_{D1}{\alpha }^2+C_{D2}\exp (C_{D3}Ma)\\ C_L=C_{L0}+C_{L1}\alpha +C_{L2}\exp (C_{L3}Ma)\end{array}$$

(4)

where $Ma$ refers to the flight Mach number, which is the ratio of the flight speed to the local sound speed. Other coefficients are as follows:

$$\begin{array}{l}{C}_{D0}=0.024, C_{D1}=7.24\times {10}^{-4}, C_{D2}=0.406, C_{D3}=-0.323\\ C_{L0}=-0.2317, C_{L1}=0.0513, C_{D2}=0.2945, C_{D3}=-0.1028\end{array}$$

(5)

2.2. Measurement Model

The measurement of the single-station radar includes the pitch angle $a$, azimuth angle $e$, and relative distance $r$. According to the geometric relationship shown in Figure 2, the following systematic observation equations can be established:

$${\mathit{Z}}_k=\left[\begin{array}{l}a\\ e\\ r\end{array}\right]=\left[\begin{array}{l}\arctan (\frac{z}{\sqrt{x^2+y^2}})\\ \arctan (\frac{y}{x})\\ \sqrt{x^2+y^2+z^2}\end{array}\right]+{\mathit{v}}_k$$

(6)

where ${\mathit{v}}_k$ is the measurement noise, ${\mathit{v}}_k~N(0,{\mathit{R}}_k)$, and ${\mathit{R}}_k$ is the noise covariance matrix.

The radar measurement equation is established in the ENU coordinate system, while the HTV-2 ballistic reference coordinate system in the previous text is the LLA coordinate system. Therefore, coordinate transformation is required for the hypersonic target state information, as shown in Figure 3.

Assuming that the coordinates of the hypersonic target are $(\mu ,\lambda ,h)$ in the LLA coordinate system and $(X,Y,Z)$ in the ECEF coordinate system, the conversion relationship is as follows:

$$\{\begin{array}{l}X=(N+h)\cos (\mu )\cos (\lambda )\\ Y=(N+h)\cos (\mu )\sin (\lambda )\\ Z=(N(1-e^2)+h)\sin (\mu )\end{array}$$

(7)

$$\{\begin{array}{l}{e}^2=\frac{a^2-b^2}{a^2}\\ N=\frac{a}{\sqrt{1-e^2{\sin }^2(\mu )}}\end{array}$$

(8)

where $N$ is the radius of the reference ellipsoid, $e$ is the eccentricity of the ellipsoid, $a$ is the long half-axis of the Earth, and $b$ is the short half-axis of the Earth.

Similarly, the LLA coordinates $({\mu }_0,{\lambda }_0,h_0)$ of the radar can be converted into ECEF coordinates $(X_0,Y_0,Z_0)$ by Equations (6) and (7). Therefore, the position information of the hypersonic target in the ENU coordinate system can be obtained:

$$\left[\begin{array}{l}x\\ y\\ z\end{array}\right]=\left[\begin{array}{l}e\\ n\\ u\end{array}\right]=T_{\mathrm{ECEF}}^{\mathrm{ENU}}\left[\begin{array}{l}X-X_0\\ Y-Y_0\\ Z-Z_0\end{array}\right]$$

(9)

$$T_{\mathrm{ECEF}}^{\mathrm{ENU}}=\left[\begin{array}{ccc}-\sin ({\lambda }_0)& \cos ({\lambda }_0)& 0\\ -\sin ({\mu }_0)\cos ({\lambda }_0)& -\sin ({\mu }_0)\sin ({\lambda }_0)& \sin ({\mu }_0)\\ \cos ({\mu }_0)\cos ({\lambda }_0)& \cos ({\mu }_0)\sin ({\lambda }_0)& \sin ({\mu }_0)\end{array}\right]$$

(10)

3. Strong Tracking CKF Based on Multiple Fading Factors

3.1. CKF Algorithm

Taking the hypersonic target motion model as the state equation and the single-station radar observation model as the measurement equation, the following nonlinear system model is formed:

$$\{\begin{array}{l}{x}_k=f(x_{k-1})+w_{k-1}\\ z_k=h(x_k)+v_k\end{array}$$

(11)

where:

$$x_k={\left[\begin{array}{cccccccccccc}x& \dot{x}& \ddot{x}& \stackrel{⃛}{x}& y& \dot{y}& \ddot{y}& \stackrel{⃛}{y}& z& \dot{z}& \ddot{z}& \stackrel{⃛}{z}\end{array}\right]}^T\in R^n$$

(12)

It refers to the flight state of the hypersonic target at time $k$, including the position, speed, acceleration, and other $n$ dimensional state information of the target. $z_k={[a e r]}^T$ is the observation vector of the radar system at time $k$, representing the measurement element of the radar. Additionally, $f(\cdot )$ is the process function and $h(\cdot )$ is the measurement function of the nonlinear tracking system. $w_k$ and $v_k$ are uncorrelated zero-mean Gaussian white noise, and the corresponding covariance matrices are $Q_k$ and $R_k$.

CKF is a nonlinear Gaussian filtering algorithm, whose core is to select the sampling points through the spherical radial cubature criterion. The traditional calculation process is as follows [17].

3.2. Multiple Fading Factor Strong Tracking CKF

The strong tracking filter (STF) can adjust the state prediction covariance matrix $P_{k+1|k}$ in real time by introducing the changing fading factor ${\lambda }_k$, it can keep the filtering residual matrix orthogonal, in order to deal with the problem of the sudden change of the target motion state and improve the robustness. The core idea is to adjust the gain matrix $K_k$ online to meet the following two conditions:

(1)

$E[({\mathit{x}}_{k+1}-{\hat{\mathit{x}}}_{k+1}){({\mathit{x}}_{k+1}-{\hat{\mathit{x}}}_{k+1})}^T]=\min$;

(2)

$\begin{array}{cc}E[{\mathit{e}}_{k+1+j}{\mathit{e}}_{k+1}^T]=0& \begin{array}{cc}k=0,1,2,\cdots & j=1,2,\cdots \end{array}\end{array}$, where ${\mathit{e}}_{k+1}={\mathit{z}}_{k+1}-{\hat{\mathit{z}}}_{k+1|k}$ is the residual at time $k+1$. That is, it is required that the residual sequences should remain orthogonal to each other at all times.

The traditional STF algorithm is based on expend Kalman filter (EKF). In this paper, the STF algorithm and CKF algorithm were combined to establish the STF algorithm for nonlinear systems. The core structure of the algorithm is the calculation of the residual:

$${\widehat{x}}_k={\widehat{x}}_{k|k-1}+K_k{e}_k$$

(13)

$$e_k=z_k-h\left({\widehat{x}}_{k|k-1}\right)$$

(14)

The fading factor ${\lambda }_{k+1}$ is calculated as follows:

$${\lambda }_k=\max \left(1,\frac{\mathrm{tr}\left[N_k\right]}{\mathrm{tr}\left[M_k\right]}\right)$$

(15)

$$N_k=V_{0,k}-H_k{Q}_{k-1}{H}_k^T-\beta R_k$$

(16)

$$M_k=H_k{P}_{k|k-1}^x{H}_k^T$$

(17)

$$H_k=\left(P_{k|k-1}^{xz,l}\right){\left(P_{k|k-1}^l\right)}^{-1}$$

(18)

$$P_{k|k-1}^{xz,l}=\frac{1}{m}{\displaystyle \sum _{i=1}^m{x}_{k|k-1}^{i,l}{h}^T\left(x_{k|k-1}^{i,l}\right)-{\widehat{x}}_{k|k-1}{h}^T\left({\widehat{x}}_{k|k-1}\right)}$$

(19)

$$V_{0,k}=\{\begin{array}{l}{e}_1{e}_1^T,k=0\hfill \\ \frac{\epsilon V_{0,k}+e_k{e}_k^T}{1+\epsilon },k>0\hfill \end{array}$$

(20)

where $P_{k|k-1}^x$ is the variance matrix of the state prediction error when the fading factor is not introduced and the system noise variance matrix is not considered; $P_{k|k-1}^l$ is the state prediction error variance matrix considering the system noise variance matrix; $\beta \ge 1$ is the weakening factor; and $\epsilon$ is the forgetting factor, usually taken as $0.95\le \epsilon \le 0.995$.

The nonlinear system variables of hypersonic targets are multi-dimensional and complex, so it is difficult to determine the optimal fading factor, and the accuracy of a single fading factor is limited. Therefore, a strong tracking cubature Kalman filter (STCKF) algorithm based on multiple fading factors is introduced:

(1)

In the time update process, the symmetric fading factor is introduced into the state prediction error covariance matrix, and $P_k^{-}$ is changed to the following formula:

$$P_k^{-}={\lambda }_{k+1}[\frac{1}{m}{\displaystyle \sum _{i=1}^m\left({\chi }_{i,k}^{-}-{\widehat{x}}_k^{-}\right){\left({\chi }_{i,k}^{-}-{\widehat{x}}_k^{-}\right)}^T]}{\lambda }_{k+1}{}^T+Q_{k-1}$$

(21)

(2)

In the measurement update process, a fading factor is introduced into the innovation covariance matrix and the measurement autocovariance matrix:

$$P_z={\lambda }_{k+1}[\frac{1}{m}{\displaystyle \sum _{i=1}^m{{\rm Z}}_{i,k}^{-}{\left({{\rm Z}}_{i,k}^{-}\right)}^T-{\widehat{z}}_k^{-}{\left({\widehat{z}}_k^{-}\right)}^T]}+R_k$$

(22)

$$P_{xz}={\lambda }_{k+1}[\frac{1}{m}{\displaystyle \sum _{i=1}^m{X}_{i,k}^{-}{\left({{\rm Z}}_{i,k}^{-}\right)}^T}-{\widehat{x}}_k^{-}{\left({\widehat{z}}_k^{-}\right)}^T]$$

(23)

Compared with [10], the algorithm in this paper adds a fading factor in the time updating process, which can adjust the prediction covariance matrix in real time and the gain matrix more obviously. With a slight increase in the amount of calculation, the algorithm accuracy can be improved.

In nonlinear systems, the system noise covariance matrix $Q_k$ is often time-varying and uncertain [24]. A suboptimal unbiased MAP constant noise estimator can be designed according to the maximum a posteriori estimation principle and the exponential weighting method. The fading factor is added to the filter to adjust the system noise, and the following recursive formula is given:

$${\widehat{Q}}_{k+1}=(1-d_{k+1}){\widehat{Q}}_{k-1}+d_{k+1}[K_{k+1}{e}_{k+1}{e}_{k+1}{}^T{K}_{k+1}{}^T+P_{k+1}-{\displaystyle \sum _{i=1}^m\frac{1}{m}\left({\chi }_{i,k}^{-}-{\widehat{x}}_k^{-}\right){\left({\chi }_{i,k}^{-}-{\widehat{x}}_k^{-}\right)}^T]}$$

(24)

where $d_k=(1-b)/(1-b^k)$, b is the noise fading factor, and the value range is $0<b<1$. Through multiple simulations, it can be seen that a better tracking filtering effect can be obtained when b is 0.89.

Substitute the selected b value into the formula to acquire:

$$Q_k=(1-d_k)Q_{k-1}+d_k[K_k{e}_k{e}_k{}^T{K}_k{}^T\hspace{1em}]k\ge 2$$

(25)

Add $Q_k$ to the multi fading factor strong tracking algorithm, the filtering value of the target is closer to the real value and the filtering accuracy is improved.

4. AIMM Algorithm

4.1. IMM Algorithm

The IMM algorithm is mainly composed of an interactor, multiple filters, a state estimation fusion device, and a model probability estimator. The main calculation process is as follows (Figure 4):

The traditional IMM algorithm has four calculation steps, it can refer to the literature [8].

4.2. AIMM Algorithm

In the IMM algorithm, model switching is a Markov process, so it is very dependent on the Markov matrix of the model transformation. In the standard IMM algorithm, the Markov matrix is given in advance. Due to the uncertainty of prior information such as the mobility of the target, it is generally determined based on experience. Therefore, for highly maneuvering targets such as hypersonic targets, the fixed Markov matrix cannot sufficiently reflect the real motion state of the target and the switching process between the models. Therefore, based on traditional IMM algorithm, we proposed to use the posterior information to adaptively modify the Markov matrix.

The model likelihood function value ${\mathsf{\Lambda }}_j\left(k\right)$ contains the residual and residual covariance between the predicted state and the real state at time $k$, so it can sufficiently reflect the matching degree between model $j$ and the current motion state. The degree of model matching is proportional to the value of the likelihood function. Therefore, based on the model likelihood function, this paper proposes a transfer probability matrix adaptive algorithm based on the likelihood function. The calculation process is as follows:

$${\phi }_{ij}\left(k\right)={(\frac{{\mathsf{\Lambda }}_j\left(k\right)}{{\mathsf{\Lambda }}_i\left(k\right)})}^{\gamma }$$

(26)

$$P_{ij}\left(k\right)=P_{ij}^{\prime }\left(k\right)/{\displaystyle \sum _{j=1}^u{P}_{ij}^{\prime }\left(k\right)}$$

(27)

$$P_{ij}^{\prime }\left(k\right)={\phi }_{ij}\left(k\right)P_{ij}\left(k-1\right)$$

(28)

where ${\phi }_{ij}\left(k\right)$ is the ratio of the model likelihood function, which is used as the adjustment factor to adjust the model conversion rate, and $\gamma$ is the correction rate, $0\le \gamma \le 1$, which is used to control the speed of correction. It can be seen from Equation (24) that at time $k$, when the likelihood function index ${\mathsf{\Lambda }}_j\left(k\right)$ is greater than ${\mathsf{\Lambda }}_i\left(k\right)$, model $j$ has a higher matching degree with the real motion than model $i$, and the probability from model $j$ to model $i$ increases after the correction. Therefore, the model with a larger likelihood function value accounts for a larger proportion in the next moment. On the contrary, the proportion of the model with a small likelihood function decreases in the next moment. Therefore, it can be concluded that this method enhances the proportion of matching models by using the measurement information to adaptively modify the transfer probability matrix.

However, it can be found through substitution calculation that Formula (25) has defects. During iterative calculation, some elements of the main diagonal may become larger and larger, and the transition probability between models may be reduced. Therefore, when the model state needs to be switched, the efficiency will be reduced. Therefore, the thresholds $M$ and $N$ are set for the main diagonal element and the non-main diagonal element, respectively, and the following improvements are made.

If the main diagonal element is too large, i.e., $P_{ii}\left(k\right)>M$, let $P_{ii}^{\prime }\left(k\right)=M$, and redistribute the non-main diagonal elements of its row according to the original ratio:

$$P_{ij}^{″}\left(k\right)=P_{ij}\left(k\right)\left[1-\frac{M-P_{ii}\left(k\right)}{{\displaystyle \sum _{j=1}^m{P}_{ii}\left(k\right)-P_{ii}^{″}\left(k\right)}}\right]$$

(29)

Similarly, if the non-main diagonal element is too small, i.e., $P_{ij}\left(k\right)<N(i\ne j)$, let $P_{ij}^{\prime }\left(k\right)=N(i\ne j)$, and redistribute the other elements of its row according to the original proportion (refer to Formula (27)).

The correction of the elements by Equation (27) still satisfies ${\displaystyle \sum _{j=1}^m{P}_{ij}\left(k\right)=1}$ and $0\le P_{ij}\left(k\right)\le 1$, and can enhance the stability and adaptability of the model transformation.

4.3. SVD Algorithm

The selection of the IMM model set needs to meet the highest matching degree between the model and the motion state at any time. The existing models include the CV model, CA model, constant turning (CT) model, Singer model, CS model, and Jerk model, which are applicable to different scenarios. Among them, the CV and CA models are applicable to the non-maneuvering or weak maneuvering movement state of the target, CT models are applicable to the turning movement of the target, the Singer and CS models can describe the normal movement, and the Jerk model, as the model with the highest order, can describe the acceleration of the target, which is mainly used to track strong maneuvering targets. Literature [25] proposed the single model to target tracking, but the hypersonic target has high maneuverability, which is difficult to tracking. Therefore, this paper selected three models with high maneuverability, namely, the Singer model, CS model, and Jerk model, to build the IMM model set.

The state equation of the target motion model can be expressed as

$$x_k=F_1{x}_{k-1}+w_{k-1}$$

(30)

where $x_{k-1}$ is the motion state of the target at the previous time, $x_k$ is the motion state of the current target, and $F_1$ is the state transition matrix.

Singer model: Singer model is time related, that is, the target acceleration can be set as a first-order time correlation function with the mean value of 0, and its corresponding function can be expressed as:

$$R_a\left(\tau \right)=E\left[a\left(t\right)a\left(t+\tau \right)\right]={\sigma }_a{}^2{e}^{-\alpha \left|\tau \right|}\begin{array}{cc}& \left(\alpha \ge 0\right)\end{array}$$

(31)

where $a(t)$ is the target acceleration, ${\sigma }_a{}^2$ is the variance, $\alpha$ is the maneuvering frequency. The probability density function can be used to solve ${\sigma }_a{}^2$:

$${\sigma }_a{}^2={\displaystyle {\int }_{-a_{\max }}^{a_{\max }}P\left(a\right)a^{2}da=\frac{a_{\max }{}^2}{3}\left(1+4P_{\max }-P_0\right)}$$

(32)

The state transition matrix and state noise covariance matrix in one dimension are as follows:

$$F_1=\left[\begin{array}{ccc}1& T& (\alpha T-1+e^{-\alpha T})/{\alpha }^2\\ 0& 1& (1-e^{-\alpha T})/\alpha \\ 0& 0& e^{-\alpha T}\end{array}\right]$$

(33)

In the formula, $\alpha$ refers to the maneuvering frequency, which is usually selected according to experience: generally, 1, 1/20, and 1/60 are taken according to the maneuverability of the target. The state noise covariance matrix can be found in [26].

CS model: The state transition matrix and state noise covariance matrix are the same as in the Singer model. Compared with the Singer model, the CS model has a control input matrix, and its discrete state equation can be expressed as:

$$x_k=F_1{x}_{k-1}+U_k\overline{a}+w_{k-1}$$

(34)

where $\overline{a}$ is the average value of acceleration, and the one-step prediction of the acceleration state can be taken as the average value of the current acceleration, i.e., ${\overline{a}}_x={\ddot{x}}_k^{-}$; according to the average value of the “current” acceleration, the variance of the current acceleration can be obtained:

$${\sigma }_a^2=\{\begin{array}{l}(4-\pi )/\pi ·{\left(a_{\max }-\overline{a}\right)}^2, \overline{a}>0\\ (4-\pi )/\pi ·{\left(\overline{a}-a_{\min }\right)}^2, \overline{a}<0\end{array}$$

(35)

where $a_{\max }$ and $a_{\min }$ are the maximum and minimum accelerations set for each axis, which are generally set according to prior conditions.

Jerk model: Based on the Singer model, zero-mean colored noise is introduced to express the acceleration. The exponential autocorrelation function of target acceleration $j\left(t\right)$ is:

$$r_j\left(\tau \right)=E\left[j\left(t\right)j\left(t+\tau \right)\right]={\sigma }_j{}^2{e}^{-\alpha \left|\tau \right|}\begin{array}{cc}& \left(\alpha \ge 0\right)\end{array}$$

(36)

where, ${\sigma }_j{}^2$ is the variance of acceleration of maneuvering target, and ${\alpha }_{\mathrm{m}}$ is the maneuvering frequency. By Laplace transformation of the above equation, we can acquire:

$$R\left(s\right)=Laplace\left\{R_j\left(\tau \right)\right\}=\frac{-2{\alpha }_{\mathrm{m}}{\sigma }_j{}^2}{\left(s-{\alpha }_{\mathrm{m}}\right)\left(s+{\alpha }_{\mathrm{m}}\right)}=H\left(s\right)H\left(-s\right)V\left(s\right)$$

(37)

where,

$$H\left(s\right)=1/\left(s+a\right)$$

(38)

$$V\left(s\right)=2{\alpha }_{\mathrm{m}}{\sigma }_j{}^2$$

(39)

Then change the above equation in time domain to obtain the differential equation as follows:

$$j\left(t\right)=-\alpha j\left(t\right)+v\left(t\right)$$

(40)

where the correlation function of white noise a is $v\left(t\right)$:

$$r_v\left(\tau \right)=2\alpha {\sigma }_j{}^2\delta \left(\tau \right)$$

(41)

The discrete state transition matrix and state noise covariance matrix are as follows:

$$F_2=\left[\begin{array}{cccc}1& T& \frac{T^2}{2}& p_1\\ 0& 1& T& q_1\\ 0& 0& 1& r_1\\ 0& 0& 0& s_1\end{array}\right]$$

(42)

The elements of the state noise covariance matrix can be found in [27].

Since the Jerk model is of the fourth order and the Singer and CS models are of the third order, it is necessary to expand the dimension of their state transition matrix and state noise variance matrix. All zero rows will appear in the traditional dimension expansion method, which leads to the fact that Cholesky in the CKF algorithm cannot be decomposed due to the loss of the positive definiteness of the matrix. Therefore, in order to enhance the stability of the calculation, SVD decomposition was considered to replace Cholesky decomposition in this paper. The specific process is as follows:

Cholesky decomposition of the covariance matrix in the time update process was rewritten into SVD decomposition:

$$P_{k-1}=U_{k-1}\left[\begin{array}{cc}{S}_{k-1}& 0\\ 0& 0\end{array}\right]U_{k-1}^T$$

(43)

where $S_{k-1}$ is a diagonal matrix, $S_{k-1}=\mathrm{diag}\left[s_1,s_2,\cdots s_n\right]$, and $s_i$ is the singular value of the matrix. The cubature point can be expressed as:

$$X_{i,k-1}={\widehat{x}}_{k-1}+U_{k-1}\sqrt{S_{k-1}}{\xi }_i, i=1,2,\cdots m$$

(44)

where ${\xi }_i=\sqrt{m/2}{\left[1\right]}_i$.

The measurement update is similar to the time update.

$$P_{k-1}^{-}=U_{k-1}^{-}\left[\begin{array}{cc}{S}_{k-1}^{-}& 0\\ 0& 0\end{array}\right]U_{k-1}^{-}{}^T$$

(45)

$$X_{i,k-1}={\widehat{x}}_{k-1}+U_{k-1}^{-}\sqrt{S_{k-1}^{-}}{\xi }_i, i=1,2,\cdots m$$

(46)

This algorithm takes into account the representation of the Jerk model in the IMM model set and the reasonable decomposition of the covariance matrix in the CKF algorithm and enhances the stability of the algorithm.

5. Simulation and Analysis

5.1. Simulation Environment Settings

(1)

Construction of real trajectory of HTV-2 trajectory

In this paper, the trajectory model of hypersonic target is simulated, and the real water drift trajectory model is constructed. Using the characteristic data of HTV-2 given by the United States, the mass of the aircraft was 907 kg, the aerodynamic reference area $S=0.4837{ \mathrm{m}}^2$, and the gravity acceleration $g$ was taken as $9.8 \mathrm{m}/{\mathrm{s}}^2$. In order to obtain the motion state of the hypersonic target, the initial test state of the target was set as follows: the geocentric distance was 6426 km, the longitude and latitude were 0 degrees, the initial velocity of the target was 17 Mach, the target trajectory inclination angle was $0.25 \mathrm{rad}$, and the trajectory deflection angle was $-0.42 \mathrm{rad}$. The angle of attack setting followed that presented in [28]:

$$\alpha =\{\begin{array}{ll}{\alpha }_{\max },\hfill & V_1<V\le V_e\hfill \\ {\alpha }_{\mathrm{mid}}+{\alpha }_{\mathrm{bal}}\sin \left[\frac{\left(V-V_{\mathrm{mid}}\right)\pi }{V_1-V_2}\right],\hfill & V_2\le V\le V_1\hfill \\ {\alpha }_{\max (K)}\hfill & V_f\le V<V_2\hfill \end{array}$$

(47)

$$\{\begin{array}{l}{\alpha }_{\mathrm{mid}}=\left({\alpha }_{\max }+{\alpha }_{\max (K)}\right)/2\hfill \\ {\alpha }_{\mathrm{bal}}=\left({\alpha }_{\max }-{\alpha }_{\max (K)}\right)/2\hfill \\ V_{\mathrm{mid}}=\left(V_1+V_2\right)/2\hfill \end{array}$$

(48)

where the maximum angle of attack ${\alpha }_{\max }$ was set to 20 degrees, and the angle of attack ${\alpha }_{\max (K)}$ was set to 11 degrees at the maximum lift–drag ratio. The target motion trajectory can be obtained through the fourth-order Runge–Kutta method simulation, as shown in the Figure 5.

(2)

Radar simulation condition setting

The coordinates in the radar longitude, latitude, altitude coordinate system were $({37.45}^{°},-{15}^{°},0 \mathrm{m})$, the detection range of the single-station radar was 400 km, the standard deviations of the angle measurement and range were 500 m and 0.5 degrees, respectively, the radar sampling period was 1 s, and the tracking period was 710 s–1060 s. The tracking track is shown in the red track part in Figure 5.

5.2. Algorithm Simulation and Analysis

(1)

Algorithm simulation verification

The strong tracking algorithm is verified against the state mutation of maneuvering target (Figure 6):

As can be seen from Figure 7, the target has a state mutation at time 40 and 135. The RMSE of the traditional IMM-CKF is obvious when the state mutation occurs, and the tracking accuracy is reduced. The strong tracking algorithm based on multiple fading factors proposed in this paper can improve the tracking accuracy by adjusting the gain matrix in real time, so that it can reduce the impact of state mutation residual. Therefore, the following simulation results can be obtained by applying the algorithm in this paper to the strong maneuver phase of the hypersonic target trajectory.

As shown in the Figure 8, three algorithms are used to simulate and verify the strong maneuvering segment of the hypersonic target. It can be seen from the trajectory tracking graph of X, Y, Z axes that the three algorithms have achieved effective tracking. The algorithm proposed in this paper has the highest tracking accuracy in the three axes. It is obvious from the z-axis tracking process that the hypersonic target has a strong maneuver in about 150 s. The algorithm proposed in this paper can reduce the impact of state mutation and improve the tracking accuracy. According to the strong maneuvering phase of the hypersonic target, the algorithm in this paper is compared with the recent literature:

Compared with the reference [13], this paper proposed adaptive IMM. The comparative simulation results are as follows (Figure 9):

From the experimental simulation results, it can be seen that the tracking accuracy of the algorithm in reference [13] is lower than AIMM in this paper, and the mean square error variance is large. Take the Jerk model as an example, it can be seen from the change trend of the trace value of the position gain matrix in the figure that at time t = 158–164 s, the trace value of the position gain matrix has a sudden change, and the model state has changed; The gain matrix $K_k$ trace value of the algorithm in this paper changes less than that of the algorithm in reference [13], so the adaptability of the model matching degree is enhanced.

(2)

Hypersonic target tracking simulation

Take the weak maneuver phase as an example, the hypersonic target tracking algorithm is simulated. The initial value of the hypersonic target was set as $x_0$ = [−719,000, 4270, 1, 0.2, 233,000, −1530, −0.6, 0.3, 199,000, 1030, −5.9, 0.5], and the target maneuvering frequency was set as $\alpha =1/20$. In the CS model, the maximum and minimum accelerations of the $x$ axis, $y$, axis and $z$ axis were, respectively, $\pm 700 \mathrm{m}/{\mathrm{s}}^2$, $\pm 500 \mathrm{m}/{\mathrm{s}}^2$, and $\pm 200 \mathrm{m}/{\mathrm{s}}^2$, and the observation noise covariance matrix was set as $R=\mathrm{diag}[\begin{array}{ccc}{0.5}^2& {0.5}^2& {500}^2\end{array}]$. In the IMM algorithm, the initial probability matrix of the model was $u_1=\left[\begin{array}{ccc}1/3& 1/3& 1/3\end{array}\right]$, and the initial Markov transition probability matrix was

$$P_1=\left[\begin{array}{ccc}0.9& 0.05& 0.05\\ 0.05& 0.9& 0.05\\ 0.05& 0.05& 0.9\end{array}\right]$$

(49)

After conducting 100 Monte Carlo simulations, the performance of the algorithm was measured by the root mean square error (RMSE) and the mean value of the RMSE.

$$RMS{E}_i=\sqrt{\frac{1}{N}{\displaystyle \sum _{k=1}^N{\left(X^i-{\widehat{X}}^i\right)}^2+{\left(Y^i-{\widehat{Y}}^i\right)}^2}+{\left(Z^i-{\widehat{Z}}^i\right)}^2}$$

(50)

where $N$ refers to the number of Monte Carlo simulations, $X^i,Y^i,Z^i$ refer to the real state of each axis in the $i$ simulation, and ${\widehat{X}}^i,{\widehat{Y}}^i,{\widehat{Z}}^i$ refer to the filter estimation value of each axis. The simulation results are as follows.

The tracking trajectory of the target is shown in Figure 6. It can be seen that the three algorithms could effectively track the hypersonic target. In order to compare the tracking accuracy of the three algorithms, the position and speed root mean square errors of each axis were simulated, as shown in Figure 10 and Figure 11. Additionally, the comparison table of the root mean square error and mean value of the three algorithms was obtained through calculation, as shown in

Table 1. Table of the mean RMSEs of the three algorithms.

Algorithm		IMM-CKF	IMM-STCKF	AIMM-STCKF
Position/m	$x$	449.2	447.9	381.9
	$y$	213.7	207	171.3
	$z$	164	158.2	97.1
	$r$	537.3	530.2	448.6
Velocity/(m/s)	$v_x$	52.7	37.9	19.4
	$v_y$	56.2	56.3	27.7
	$v_z$	51.9	49.5	41.4
	$v$	104.1	92.7	66.9

. The model transition probability of IMM-CKF is shown in Figure 12a, and the model transition probability of the algorithm in this paper is shown in Figure 12b.

As defenders, the velocity tracking error and position tracking error are the main indicators for HTV-2 tracking performance [29,30,31]. According to Figure 10 and Figure 11 and Table 1, the RMSE of the x, y, z positions and speed can be found:

(1)

The IMM-CKF and IMM-STCKF algorithms converged faster because the setting error of the initial value of the axis position and speed was small. Due to the large error of the initial value, the convergence time of the two algorithms was about 40 s. In contrast, the AIMM-STCKF algorithm has model adaptability, so the convergence time was significantly faster. Overall, it was about 25 s, and the convergence speed was increased by about 37.5%.

(2)

The AIMM-STCKF algorithm proposed in this paper significantly improved the accuracy of the algorithm after adding the Markov adaptive matrix. Among the three algorithms, the root mean square error of the speed and position was the smallest, the tracking accuracy was the highest and the variance was small. Overall, the position tracking accuracy was about 16.51% higher than that of the IMM-STCKF algorithm, and the speed tracking accuracy was about 35.74%.

(3)

Since there was no obvious mutation in the state of the hypersonic target, the accuracy of the strong tracking algorithm was low. In contrast, the strong tracking algorithm based on adaptive model interaction could accelerate the model transformation and enhance the adaptive matching ability, so the convergence speed and tracking accuracy were better.

By comparing the probability transfer diagram of the model, the following can be found:

(1)

As the maneuverability of the hypersonic target decreased, the probability of the Singer model with weak maneuverability gradually increased. As the maneuverability of the hypersonic target increased, the probability of the Jerk model with strong maneuverability gradually increased.

(2)

It can be found through comparison that the algorithm in this paper enhanced the proportion of matching models in the IMM algorithm, reduced the probability variance in the model, and made the model transformation more robust.

6. Concluding Remarks

A strong tracking adaptive IMM algorithm based on CKF was proposed to solve the problem of the low tracking accuracy or even failure of the traditional IMM algorithm under the high maneuverability of hypersonic targets.

Firstly, the six-degree-of-freedom equation was established for the dynamics of HTV-2. The jump glide trajectory of HTV-2 was obtained by using the fourth-order Runge–Kutta method and the angle of attack model, and the radar measurement model was established according to the coordinate axis transformation relationship. After analyzing the structure of CKF, the gain matrix was adjusted in real time by adding multiple pairs of fading factors, and the strong tracking filter was combined with CKF to improve the filtering accuracy when the model is out of order. Aiming at the problem that the model dimension in the IMM is expanded, which makes Cholesky decomposition impossible in the CKF algorithm, SVD decomposition was proposed to replace Cholesky decomposition. Aiming at the shortcomings of the IMM in estimating the Markov matrix, an IMM with an adaptive transition probability was proposed. The transition probability was adaptively modified by posterior information to enhance the proportion of matching models.

The multi-fading factor strong tracking algorithm proposed in this paper adjusts the filtering gain in real time and improves the tracking accuracy. The adaptive interaction multi-model proposed in this paper improves the matching degree of the motion state of the hypersonic target and speeds up the convergence of the algorithm. The SVD decomposition of CKF in this paper improves the stability of the algorithm. The simulation results show that the proposed algorithm enhanced the proportion of matching models in the IMM and improved the tracking accuracy by 16.51% and the convergence speed by 37.5%. It has engineering application value.