A Novel Adaptive Robust Cubature Kalman Filter for Maneuvering Target Tracking with Model Uncertainty and Abnormal Measurement Noises

Abstract

The features of measurement and process noise are directly related to the optimal performance of the cubature Kalman filter. The maneuvering target model’s high level of uncertainty and non-Gaussian mean noise are typical issues that the radar tracking system must deal with, making it impossible to obtain the appropriate estimation. How to strike a compromise between high robustness and estimation accuracy while designing filters has always been challenging. The H-infinity filter is a widely used robust algorithm. Based on the H-infinity cubature Kalman filter (HCKF), a novel adaptive robust cubature Kalman filter (ARCKF) is suggested in this paper. There are two adaptable components in the algorithm. First, an adaptive fading factor addresses the model uncertainty issue brought on by the target’s maneuvering turn. Second, an improved Sage–Husa estimation based on the Mahalanobis distance (MD) is suggested to estimate the measurement noise covariance matrix adaptively. The new approach significantly increases the robustness and estimation precision of the HCKF. According to the simulation results, the suggested algorithm is more effective than the conventional HCKF at handling system model errors and abnormal observations.

1. Introduction

Target tracking is estimating motion parameters such as the tracked target’s position and velocity in real time through measurement data [1]. It is a prerequisite for following tasks like target recognition and data fusion. Radar tracking is the primary method of target tracking. There are numerous methods, including infrared tracking, laser tracking, sonar tracking, and other types of tracking, due to the development of associated technologies, the diversity of the types of tracked objects, and the complexity of motion [2]. Target tracking plays an essential role in many fields, such as vehicle monitoring [3], space attitude perception [4], and missile guidance [5]. Many literature studies focus on efficiently tracking high-speed maneuvering targets with fast flight speeds and complex trajectories [6,7,8,9,10].

The Kalman filter, an ideal recursive estimator created for linear Gaussian systems, is the approach for target tracking that is most frequently used [11]. However, the system is frequently nonlinear in practical engineering applications, particularly the measurement equation. Scholars have proposed the extended Kalman filter (EKF), based on Taylor series expansion, and the unscented Kalman filter (UKF), based on sampling distance estimation, to suit the needs of practical nonlinear filtering better [12]. However, the EKF has disadvantages such as poor stability, imprecision, and sluggish target maneuver reaction. When faced with high-dimensional system states, the UKF is vulnerable to dimensional disasters, which cause a sharp reduction in filtering performance and cause it to fall behind the target [13]. Regarding numerical accuracy and operational stability, the cubature Kalman filter (CKF) [14] suggested by Arasaratnam is preferable to the algorithm mentioned above. The sigma point sampling approach and weight distribution of the UKF are optimized with the CKF using spherical integral and radial integral criteria, which resolves the dimension disaster issue. The CKF is a case of UKF when the free parameter equals zero. It offers a rigorous theoretical foundation for zero degrees of freedom in high-dimensional state estimation. However, maneuvering targets usually do not maintain only one motion mode, and the trajectories of maneuvering targets within the military domain exhibit a heightened level of complexity, compounded by an unpredictable operational environment, which often gives rise to indeterminate measurement noise. Therefore, the CKF, like other nonlinear filters, will suffer a severe decline in tracking accuracy when tracking targets under such conditions.

An H-infinity filter not based on the assumption of signal spectral characteristics is proposed in [15] to increase the robustness of navigation and tracking systems. The H-infinity filter is designed to minimize the impact of the worst disturbance on the estimation error by incorporating the H-infinity norm into the filter design. In [16,17], the H-infinity filter, which can be applied only to linear systems, is applied to nonlinear systems, maintaining the advantages of the CKF and H-infinity filter. However, various circumstances, including model uncertainty and external interference, will prevent information from fusing, lowering the H-infinity filter’s estimation accuracy.

Academics have presented some adaptive techniques to further enhance the filtering algorithm’s estimation performance. Zhou proposed a strong tracking filter (STF) that can maintain the residual sequence orthogonal to solve the system model uncertainty problem [18]. The STF method can be more robust when the system parameters or state change. The authors of [19,20] proposed an improved STF method based on a point estimation nonlinear Kalman filter and confirmed its effectiveness in navigation tracking systems. The Sage–Husa adaptive Kalman filter (SHAKF), based on maximum a posteriori estimation, is a straightforward and helpful solution for the uncertain noise in the measurement process [21]. The SHAKF algorithm can estimate the statistical characteristics of noise in real time and reduce the filtering divergence. However, ensuring that SHAKF will calculate the noise covariance matrix in a positive, definite manner can be challenging, which could result in filtering divergence. Performance will also decrease if measurement noise is mixed with non-Gaussian noise. The Mahalanobis distance (MD) is an index for anomaly statistical detection [22]. The robust estimate approach can be used for robust filtering when the observation value is abnormal, which can successfully lessen the impact of abnormal model deviation and abnormal measurement [23].

The research presented in this paper develops an adaptive fading factor and noise covariance estimation method with a robust strategy based on the traditional HCKF to solve the issues mentioned earlier. A low-cost STF calculation method is adopted to solve the problem of insufficient tracking accuracy caused by model errors. The MD-improved Sage–Husa estimation method is used to estimate abnormal or non-Gaussian measurement noise, further improving the robustness and tracking accuracy of the HCKF. The major contributions of this paper can be summarized as follows:

• A low-cost strong tracking fading factor is proposed. This method avoids the calculation of the Jacobian matrix and additional sampling operations, significantly reducing the algorithm’s time complexity. At the same time, it also solves the problem of STF failure caused by the large difference in the order of magnitude of each dimension of the measurement system.
• An MD-based method is proposed to correct the noise covariance of the Sage–Husa estimate. Sage–Husa estimation decreases measuring system error by estimating measurement noise in real time, but this method can easily lead to non-positive definite matrices and filter divergence. The MD approach is used in this study to construct the expansion factor, which is used to correct the measurement noise updated by the Sage–Husa estimate, limit the interference of outliers on the filter results, and improve the algorithm’s robustness.
• For the composite problems of model mismatch and measurement noise anomaly that may occur in maneuvering target tracking, a new robust adaptive cubature Kalman filter is proposed. The new method, which is based on the HCKF algorithm, introduces the two new, enhanced methods mentioned above, and suppresses the influence of system mutation and non-Gaussian noise. The simulation results show that the proposed method achieves a similar impact as the interacting-multiple model CKF (IMMCKF) algorithm in the face of a sudden change in maneuvering target motion. It also exhibits improved robustness and tracking accuracy in the presence of anomalous measurement noise compared with the CKF, HCKF, and Sage–Husa CKF.

The rest of this article is organized as follows. Section 2 describes the system model and the problems to be solved. Section 3 gives the principle and process steps of the traditional HCKF algorithm. Section 4 introduces an improved adaptive robust cubature Kalman filter for high-speed maneuvering target tracking. In Section 5, the proposed algorithm’s effectiveness is verified with various simulation experiments. Finally, the conclusion is followed in Section 6.

2. Problem Formulation

Consider the following nonlinear discrete system:

$$X_k=f\left(X_{k-1}\right)+v_{k-1}$$

(1)

$$Z_k=h\left(X_k\right)+w_k$$

(2)

where $X_k\in {\mathbf{R}}^n$ is the system state quantity at time k; $Z_k\in {\mathbf{R}}^m$ is the system measurement at time k; $f\left(·\right)$ and $h\left(·\right)$ are the nonlinear functions of the system state equation and the measurement equation, respectively; $v_{k-1}$ and $w_k$ represent the process noise and measurement noise of the system; the mean is zero; and the variances are $Q_{k-1}$ and $R_k$, respectively.

According to the classic CKF, the system must be accurately modeled, and the process and measurement noise must be Gaussian white noise with well-known statistical properties. The process state equation will unavoidably vary in the radar monitoring system due to the maneuvering target’s easily modifiable trajectory. The measurement noise is typically non-standard Gaussian noise due to external environment disturbance, and its statistical features cannot be recognized in real time. The traditional filter algorithm will not be able to estimate the state accurately in this situation, and divergence issues may even arise.

The problems solved in this paper are described as follows:

• The high-speed maneuvering target has uncertain reentry motion, which is the main reason for the significant decrease in tracking accuracy. An improved algorithm is needed to solve the model error caused by maneuvering motion.
• The statistical characteristics of non-Gaussian noise cannot be accurately estimated in real time, and it is necessary to improve the system’s resistance to unknown noise as much as possible.
• Due to parameter settings, the proposed improvement approaches would interact with one another during implementation, necessitating specific designs to avoid system interference.

3. The H-Infinity Cubature Kalman Filter

The H-infinity filter is based on game theory, which aims to minimize the estimation error for all disturbances with bounded energy. It is a particular form of Kalman filter [15]. For the nonlinear discrete system proposed in the previous section, the H-infinity filter design idea is that when $P_k$, $Q_{k-1}$, and $R_k$ reach the upper limit, there is a specific cost function to minimize [24]:

$$J_{\mathrm{\infty }}=\frac{\sum _{k=1}^N|X_k-{\widehat{x}}_k{|}_{P_k^{-1}}^2}{|X_0-{\widehat{x}}_0{|}_{P_0^{-1}}^2+\sum _{k=1}^N\left(|v_k{|}_{Q_k^{-1}}^2+|w_k{|}_{R_k^{-1}}^2\right)}$$

(3)

where $P_k$ represents the posterior covariance matrix and $P_0$ and $X_0$ are the initial covariance matrix and the system’s initial state, respectively. The standard symbol operations used in the formula are as follows: $|x{|}_A^2=x^{T}Ax$.

The purpose of the H-infinity filter is to estimate ${\widehat{x}}_k$ minimized $J_{\mathrm{\infty }}$. In general, the analytical solution of the optimal H-infinity filter problem is challenging to obtain, so the suboptimal solution is generally sought [25]. In the worst case, the boundary of $J_{\infty }$ satisfies:

$$\sup \left(J_{\mathrm{\infty }}\right)<{\gamma }^2$$

(4)

where sup represents the upper bound and $\gamma$ is a positive scalar parameter that limits the estimation error of the system due to uncertainty. Based on Equation (4), the designer must find ${\widehat{x}}_k$ such that $\sup \left(J_{\mathrm{\infty }}\right)<{\gamma }^2$ holds for any perturbation in $v_k$, $w_k$, and $X_0$.

The cubature Kalman filter is a nonlinear filter based on the third-order spherical-radial cubature rule proposed by Arasaratnam [14]. The CKF regards the estimation of nonlinear equations as the estimation of a probability distribution based on Bayesian estimation, which dramatically simplifies the nonlinear problem and can be well applied to high-dimensional nonlinear system filtering [26]. In order to apply H-infinity filtering to nonlinear systems, a cubature Kalman filter (HCKF) based on H-infinity filtering is proposed [16,27]. The algorithm steps are as follows:

Prediction Update:

(1)

The system’s initial predictive value and covariance matrix are set to ${\widehat{x}}_0$ and $P_0$. Calculate 2n cubature points $X_{i,cub}$ from ${\widehat{x}}_{k-1}$ and $P_{k-1}$:

$$X_{i,cub}=S_{k-1}{\xi }_i+{\widehat{x}}_{k-1}$$

(5)

where $S_{k-1}$ is obtained with Cholesky decomposition of $P_{k-1}$. The cubature points ${\mathsf{\xi }}_i=\sqrt{n}{\left[I_n,-I_n\right]}_i$ and $I_n$ represent the n-dimensional unit matrix.

The obtained cubature points $X_{i,cub}$ are propagated through the nonlinear state function:

$$X_{i,k|k-1}^{*}=f\left(X_{i,cub}\right)$$

(6)

(2)

The prior prediction state value ${\widehat{x}}_{k|k-1}$ and the prediction covariance matrix $P_{k|k-1}$ are calculated.

$${\widehat{x}}_{k|k-1}=\omega \sum _{i=1}^{2n}{X}_{i,k|k-1}^{*}$$

(7)

$$P_{k|k-1}=\omega \sum _{i=1}^{2n}{X}_{i,k|k-1}^{*}{\left(X_{i,k|k-1}^{*}\right)}^T-{\widehat{x}}_{k|k-1}{\widehat{x}}_{k|k-1}^T+Q_{k-1}$$

(8)

where $\omega =1/2n$ is the cubature points’ weight.

Measurement Update:

(3)

The cubature points are calculated again according to the prior prediction ${\widehat{x}}_{k|k-1}$ and the prediction covariance matrix $P_{k|k-1}$ obtained in the previous step, and ${\widehat{z}}_{k|k-1}$ is obtained by nonlinear measurement function propagation.

$$P_{k|k-1}=S_{k|k-1}{S}_{k|k-1}^T$$

(9)

$$X_{i,cub}^{*}=S_{k|k-1}{\xi }_i+{\widehat{x}}_{k|k-1}$$

(10)

$$Z_{i,k|k-1}^{*}=h\left(X_{i,cub}^{*}\right)$$

(11)

$${\widehat{z}}_{k|k-1}=\omega \sum _{i=1}^{2n}{Z}_{i,k|k-1}^{*}$$

(12)

(4)

The autocorrelation covariance matrix $P_{zz}$ and the cross-correlation covariance matrix $P_{xz}$ are as follows:

$$P_{zz}=\omega \sum _{i=1}^{2n}{Z}_{i,k|k-1}^{*}{\left(Z_{i,k|k-1}^{*}\right)}^T-{\widehat{z}}_{k|k-1}{\widehat{z}}_{k|k-1}^T+R_k=P_{ee}+R_k$$

(13)

$$P_{xz}=\omega \sum _{i=1}^{2n}{X}_{i,k|k-1}^{*}{\left(Z_{i,k|k-1}^{*}\right)}^T-{\widehat{x}}_{k|k-1}{\widehat{z}}_{k|k-1}^T$$

(14)

(5)

Calculate the Kalman gain $K_k$ according to $P_{zz}$ and $P_{xz}$:

$$K_k=P_{xz}{P}_{zz}^{-1}$$

(15)

(6)

Finally, the update state ${\widehat{x}}_k$ and update covariance matrix $P_k$ of the system are calculated:

$${\widehat{x}}_k={\widehat{x}}_{k|k-1}+K_k\left(Z_k-{\widehat{z}}_{k|k-1}\right)$$

(16)

$$P_k=P_{k|k-1}-\left[P_{xz},P_{k|k-1}\right]R_{e,k}^{-1}{\left[P_{xz},P_{k|k-1}\right]}^T$$

(17)

where

$$R_{e,k}=\left[\begin{array}{cc}{P}_{zz}+R_k& P_{xz}^T\\ P_{xz}& P_{k|k-1}-{\gamma }^{2}I\end{array}\right]$$

The threshold $\gamma$ determines the robustness of the filter. When it approaches infinity, the HCKF degenerates into the ordinary CKF. In other words, the Kalman filter is a particular case of an H-infinity filter when its performance boundary is infinite. The adversary in Kalman filtering is considered irrelevant from the game theory standpoint. The probability distribution function of the noise is known when designing a Kalman filter, and it may be used to create a statistically optimal state estimation. At the same time, the adversary does not alter this probability distribution function to degrade the performance of the estimated state. Therefore, while the Kalman filter minimizes the variance of the estimation error, it does not guarantee the estimation error in the worst case; that is, it does not ensure the boundary of the cost function in Equation (3). The H-infinity filter equation has a more unambiguous interpretation than the Kalman filter equation. It is a filter that evaluates the worst-case scenario. The opponent can maximize the estimating error by just using the infinite $X_0$, $w_k$, and $v_k$, which makes the game unfair. Thus, when defining $J_{\infty }$, $X_0-{\widehat{x}}_0$, $w_k$, and $v_k$ are in the denominator. Even if using infinite $X_0$, $w_k$, and $v_k$ increases the estimation error, the $J_{\infty }$ may not increase because the denominator also increases. Because the form of $J_{\infty }$ prohibits the adversary from utilizing cruel means to maximize the estimation error, the architecture of the H-infinity filter is robust.

4. Adaptive Robust H-Infinity Cubature Kalman Filter

The HCKF requires that all parameters and noise statistics of the system model must be known, which is almost impossible in actual maneuvering target tracking. The trajectory of a maneuvering target usually cannot be represented by a single kinematic model, and the measurement system’s noise is usually heavy tailed. The above problems will lead to a decrease in the estimation accuracy of the filter.

In order to compensate for the modeling error and improve the ability of the filter to deal with abnormal measurements, we introduce the STF and the Sage–Husa noise estimator method based on Section 3 and propose a new adaptive robust cubature Kalman filter (ARCKF). The STF corrects the prediction error covariance matrix by introducing an adaptive fading factor to compensate for the modeling error and maintain the excellent performance of the filter [28]. The Sage–Husa method is based on maximum a posteriori estimation, which has a simple structure and good real-time performance. It is often used to estimate the statistical characteristics of noise [29].

4.1. Fading Factor Based on Strong Tracking Filter

The strong tracking filter [18] is an algorithm based on EKF. This algorithm is based on the orthogonality principle of output residual sequence. It introduces the adaptive fading factor into the state prediction covariance matrix, which solves the problem that the EKF cannot converge when the model error is significant. The strong tracking filter mainly modifies the prediction covariance matrix $P_{k|k-1}$ by satisfying the following two equality conditions:

$$E\left[\left(X_k-{\widehat{x}}_k\right){\left(X_k-{\widehat{x}}_k\right)}^T\right]=min$$

(18)

$$E\left({\epsilon }_{k+j}{\epsilon }_k^T\right)=0,\left(k=\mathrm{0,1},\dots ;j=\mathrm{1,2},\dots \right)$$

(19)

where ${\widehat{x}}_k$ is the state estimation at time k, ${\epsilon }_k=Z_k-{\widehat{z}}_{k|k-1}$ is the output residual at time k, and ${\widehat{z}}_{k|k-1}$ is the observation prediction at time k. Equations (18) and (19) are sufficient conditions for STF. These two conditions are called the orthogonality principle. The filtering performance index under the minimum variance estimation criterion is given by Equation (18), which states that the output of the filter must meet the minimal mean square error; Equation (19) shows that the residual sequence is orthogonal. That is, the residual sequences at different times are not related. When the filter has a model mismatch problem, the estimated value of the system state deviates from the actual state of the system, which is reflected in the residual sequence. By adjusting the Kalman gain in real time to meet Equation (19), it can ensure that the residual sequence always has Gaussian white noise characteristics and improve the tracking ability of the system; when the system model is accurate, Equation (19) is satisfied, and the filter is equivalent to the original filter.

Equation (19) is the core of the orthogonality principle. Equation (18) represents the design standard the original filter must meet. When the original filter is combined with the criteria of Equation (19), it exhibits the features of the STF. The STF modifies the prediction covariance matrix to update the Kalman gain using the adaptive fading factor${\lambda }_k$. The following is the computation method [28]:

$${\lambda }_k=\max \left({\lambda }_0,1\right)$$

(20)

$${\lambda }_0=\frac{tr\left[N_k\right]}{tr\left[M_k\right]}$$

(21)

$$N_k=V_k-C_k{Q}_{k-1}{C}_k^T-R_k$$

(22)

$$M_k=C_k{F}_k{P}_{k-1}{F}_k^T{C}_k^T$$

(23)

$$V_k=\left\{\begin{array}{cc}{\epsilon }_1{\epsilon }_1^T& \left(k=1\right)\\ \frac{\rho V_k+{\epsilon }_k{\epsilon }_k^T}{1+\rho }& \left(k>1\right)\end{array}\right.$$

(24)

where $F_k$ and $C_k$ are the Jacobian matrices after the first-order Taylor expansion of the state function and the measurement function, respectively; $\rho \left(0<\rho \le 1\right)$ represents the forgetting factor, usually 0.95. The computational complexity of Jacobian matrix calculation is high, and this method, using first-order Taylor expansion approximation, only applies to weakly nonlinear systems, which significantly limits the application scenarios of the STF. Moreover, the fading factor calculates the ratio of the eigenvalues of each dimension element of the measurement residual matrix to the theoretical output value of the filter. This calculation method brings problems to the measurement system with significant differences in the order of magnitude of each dimension. For example, the angle residual in the radar tracking system is much smaller than the distance residual. When the maneuvering target is mainly maneuvering at the angle, the small change in the angle residual may lead to the failure of the fading factor. Therefore, it is necessary to improve the calculation of the fading factor to better cope with the above problems.

Applying statistical linearization to the nonlinear observation equation [30], we can obtain

$$z_k=H_k\left(x_k-{\widehat{x}}_{k|k-1}\right)+{\widehat{z}}_{k|k-1}+v_k$$

(25)

where $H_k=P_{xz}^T{P}_{k|k-1}^{-1}$ is the statistical regression equation.

According to the theorem in [18], the sufficient condition for the orthogonality of the output residual sequence is:

$$P_{k|k-1}{H}_k^T-K_k{V}_k=0$$

(26)

Using the definition of the cross-correlation covariance matrix $P_{xz}$, we can obtain:

$$\begin{array}{cc}\hfill P_{xz}& =E\left[\left(x_k-{\widehat{x}}_{k|k-1}\right){\left(z_k-{\widehat{z}}_{k|k-1}\right)}^T\right]\hfill \\ & =E\left[\left(x_k-{\widehat{x}}_{k|k-1}\right){\left(H_k\left(x_k-{\widehat{x}}_{k|k-1}\right)+v_k\right)}^T\right]\hfill \\ & =E\left[\left(x_k-{\widehat{x}}_{k|k-1}\right){\left(x_k-{\widehat{x}}_{k|k-1}\right)}^T{H}_k^T\right]=P_{k|k-1}{H}_k^T\hfill \end{array}$$

(27)

Equation (26) is rewritten as:

$$P_{xz}-K_k{V}_k=P_{xz}-P_{xz}{P}_{zz}^{-1}{V}_k=0$$

(28)

$$P_{zz}-V_k=0$$

(29)

It is known that the autocorrelation covariance matrix has the form

$$\begin{array}{cc}\hfill P_{zz}& =E\left[\left(z_k-{\widehat{z}}_{k|k-1}\right){\left(z_k-{\widehat{z}}_{k|k-1}\right)}^T\right]\hfill \\ & =E\left[\left(H_k\left(x_k-{\widehat{x}}_{k|k-1}\right)+v_k\right){\left(H_k\left(x_k-{\widehat{x}}_{k|k-1}\right)+v_k\right)}^T\right]\hfill \\ & =H_k{P}_{k|k-1}{H}_k^T+R_k\hfill \end{array}$$

(30)

The fading factor ${\lambda }_k$ is introduced to correct $P_{k|k-1}$ to ensure that the performance index of Equation (29) is satisfied. Therefore, Equation (29) is further rewritten as:

$${\lambda }_k\left(P_{zz}-R_k\right)+R_k=V_k$$

(31)

Redefining $M_k$ and $N_k$ in accordance with Equation (31),

$$N_k=V_k-R_k$$

(32)

$$M_k=P_{zz}-R_k=P_{ee}$$

(33)

where $P_{ee}$ was defined in Equation (13).

Therefore, rather than using Equations (22) and (23), we can determine the fading factor ${\lambda }_k$ using Equations (32) and (33). The original method of calculating the fading factor is to compare the trace of the measurement residual to the theoretical output value (see Equation (21)). Algorithm failure may happen when measurement systems have significant variances in the order of magnitude of each dimension. In order to ensure that the fading factor has the same sensitivity to different measurement dimensions, the maximum ratio of different measurement dimensions is considered the fading factor. Then Equation (21) is changed to

$${\lambda }_0=\max \left\{\frac{{\left[diag\left(N_k\right)\right]}_i}{{\left[diag\left(M_k\right)\right]}_i}\right\},i=\mathrm{1,2},\dots ,m$$

(34)

4.2. Robust Adaptive Strategy for Measurement Noise Estimator

The adaptive method for estimating the $R_k$ matrix is a scaling method based on residual covariance. It is calculated based on the difference between the actual measured value obtained at time k and its estimated value. The recursive formula of the measurement noise is expressed as [31]:

$$R_k=\left(1-d_{k-1}\right)R_{k-1}+d_{k-1}\left({\epsilon }_k{\epsilon }_k^T-P_{ee}\right)$$

(35)

where

$$d_k=\frac{1-b}{1-b^{k+1}},0<b<1$$

(36)

$d_k$ is a weighted coefficient to measure the influence of noise on the observed value, and b represents the forgetting factor. The larger the value is, the stronger the influence of the last measurement is. However, if the value of b is too small, the predicted noise will oscillate. The forgetting factor b in this method is often set between 0.95 and 0.99 based on numerical simulation attempts or experience (the frequency of changes in noise statistics). The noise estimation in Equation (36) still needs to satisfy the Gaussian distribution, and its resistance to abnormal noise needs to be increased. In order to further improve the robust performance of the HCKF, the square of the Mahalanobis distance from $Z_k$ to ${\widehat{z}}_{k|k-1}$ is considered the criterion [32]. The formula is as follows:

$$M{D}_k^2={\left(\sqrt{{\left(Z_k-{\widehat{z}}_{k|k-1}\right)}^T{\left(P_{zz}\right)}^{-1}\left(Z_k-{\widehat{z}}_{k|k-1}\right)}\right)}^2={\epsilon }_k^T{\left(P_{ee}+R_k\right)}^{-1}{\epsilon }_k$$

(37)

If the actual criterion index satisfies $M{D}_k^2>{\chi }_{n,\alpha }^2$, then the observed value $Z_k$ is marked as an outlier, and the expansion factor $\mu$ also amplifies the covariance of the observed noise; that is,

$${\tilde{R}}_k=\mu R_k$$

(38)

at this time should satisfy

$${MD}_k^2={{\epsilon }_k^T\left(P_{ee}+{\tilde{R}}_k\right)}^{-1}{\epsilon }_k={\chi }_{n,\alpha }^2$$

(39)

where ${\chi }_{n,\alpha }^2$ is the statistical detection threshold conforming to the chi-square distribution and ${\chi }_{n,\alpha }^2={\chi }^2\left(n\right)$. The problem of solving ${\mu }_k$ can be transformed into the following nonlinear equation:

$$f\left({\mu }_k\right)={{\epsilon }_k^T\left(P_{ee}+{\tilde{R}}_k\right)}^{-1}{\epsilon }_k-{\chi }_{n,\alpha }^2=0$$

(40)

Considering the use of a Newton iterative method to solve ${\mu }_k$ in Equation (40), the recursive relationship can be expressed as follows:

$${\mu }_k\left(i+1\right)={\mu }_k\left(i\right)+\frac{M{D}_k^2\left(i\right)-{\chi }_{n,\alpha }^2}{{{\epsilon }_k^T\left({\tilde{P}}_{zz}\left(i\right)\right)}^{-1}{R}_k{\left({\tilde{P}}_{zz}\left(i\right)\right)}^{-1}{\epsilon }_k}$$

(41)

where i represents the number of iterations, ${\tilde{P}}_{zz}\left(i\right)=P_{ee}+{\mu }_k\left(i\right)R_k$, we set the initial value of the iteration to ${\mu }_k\left(0\right)=1$, and the iteration ends only when the decision condition ${\mu }_k\left(i\right)\le {\chi }_{n,\alpha }^2$ is satisfied. It can be considered that $\alpha$ is set to 0.99, which means that the filter’s efficiency is 99%. According to the look-up table method, the chi-square distribution ${\chi }_{\mathrm{2,0.99}}^2$ of 2 degrees of freedom can be determined to be 9.21.

It should be noted that there are many reasons for the abnormal state prediction value. The system’s modeling error and abnormal measurement noise will affect the final state prediction value. The wrong use of the correction matrix ${\tilde{R}}_k$ is likely to lead to misjudgment, so that the STF algorithm and the robust estimation method of measurement noise work simultaneously. However, it is easy to lead to the divergence of the filter. Therefore, the calculation of the adaptive fading factor in Equation (34) uses the original $R_k$instead of the modified ${\tilde{R}}_k$ to prevent the ${\tilde{R}}_k$ based on the Mahalanobis distance correction from affecting the calculation in the presence of kinematic model errors.

Remark 1.

Adding a fading factor and adjusting measurement noise are performed using a scaling coefficient and positive definite matrix operation. Therefore, both the new approach and the conventional HCKF may guarantee the positive definiteness of the error covariance matrix. Positive definiteness may be lost in practice due to errors brought about by arithmetic operations carried out on a digital computer with a finite word length. The number-sensitive operations that are most likely to destroy the positive definiteness of the covariance include matrix square-rooting [see (5) and (9)] and matrix inversion [see (15) and (17)]. When working with high-precision measurement systems, numerically ill-conditioned is another possibility that could lead to a non-positive definite covariance matrix. The most widely utilized methods to lessen the negative consequences that could cause unstable or even divergent behavior include swapping out the CKF algorithm for SCKF and substituting Cholesky decomposition with SVD decomposition. This paper does not examine this section of the content because it is outside the purview of the work.

In summary, we give the adaptive robust correction strategy of the HCKF, and the specific process steps of the algorithm are shown in Figure 1.

5. Numerical Examples and Analysis

In order to verify the robustness and accuracy of the ARCKF algorithm proposed in this paper, three possible situations in maneuvering target tracking are simulated and compared with the estimation results of the CKF [14], HCKF [27], and Sage–Husa cubature Kalman filter (SHCKF) [29]. It should be noted that in order to maintain robustness, the conventional SHCKF method cannot concurrently perceive process noise and measurement noise. Therefore, we consider improving the SHCKF with the method in Reference [8]. Although this approach trades some noise perception accuracy, it can guarantee that the algorithm will maintain robustness and perform comparative experiments more effectively. In the experiment, 100 Monte Carlo simulations were run to reflect the tracking capability of these algorithms as accurately as possible.

In the simulation experiment, the radar position is set as the origin of the rectangular coordinate system, and the measurement information of the radar is distance and azimuth; then, the corresponding measurement equation is

$$Z_k=\left[\begin{array}{c}{r}_k\\ {\theta }_k\end{array}\right]=\left[\begin{array}{c}\sqrt{x^2+y^2}\\ \arctan \frac{y}{x}\end{array}\right]+w_k$$

(42)

The state vector of the system is $X={\left[x,v_x,y,v_y,\omega \right]}^T$, $x,y$ is the position of the maneuvering target in x direction and y direction, $v_x,v_y$ is the velocity of the maneuvering target in x direction and y direction, and $\omega$ is the turning rate. In most cases, the maneuvering target always maintains uniform motion during the cruise phase, so the system model considering the filtering algorithm is set to the uniform motion model:

$$F=\left[\begin{array}{ccccc}1& T& 0& 0& 0\\ 0& 1& 0& 0& 0\\ 0& 0& 1& T& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right]$$

(43)

The relevant parameters of the simulation experiment are given in

Table 1. Parameters for simulation.

Parameter	Corresponding Value
Number of Monte Carlo simulations	100
Discrete sampling period	T = 1 s
Process noise intensities	$q_1=0.1{\mathrm{m}}^2{\mathrm{s}}^{-3},q_2=1.75\times {10}^{-4}{\mathrm{s}}^{-3}$
Initial process noise covariance matrix	$Q_{k-1}=diag\left(\left[q_1{M}_1,q_1{M}_1,q_{2}T\right]\right),M_1=\left[T^3/3,T^2/2;T^2/2,T\right]$
Measurement noise intensities	${\sigma }_r=10\mathrm{m},{\sigma }_{\theta }=3.1\mathrm{m}\mathrm{r}\mathrm{a}\mathrm{d}$
Initial measurement noise covariance matrix	$R_k=diag\left(\left[{\sigma }_r^2,{\sigma }_{\theta }^2\right]\right)$
Initial state of the maneuvering target	$X_0={\left[1000\mathrm{m},200\mathrm{m}/\mathrm{s},8000\mathrm{m},10\mathrm{m}/\mathrm{s},5°/\mathrm{s}\right]}^T$
Initial state covariance matrix	$P_0=diag\left(\left[100{\mathrm{m}}^2,10{\mathrm{m}}^2/{\mathrm{s}}^2,100{\mathrm{m}}^2,10{\mathrm{m}}^2/{\mathrm{s}}^2,0.001\mathrm{r}\mathrm{a}{\mathrm{d}}^2/{\mathrm{s}}^2\right]\right)$

. The simulation results select the root-mean-square error (RMSE) as the performance index. Position RMSE and velocity RMSE are defined as:

$$\begin{array}{c}RMS{E}_{pos}^k=\sqrt{\frac{1}{N}{\displaystyle \sum _{j=1}^N}{\left(x_k^j-{\widehat{x}}_k^j\right)}^2+{\left(y_k^j-{\widehat{y}}_k^j\right)}^2}\\ RMS{E}_{vel}^k=\sqrt{\frac{1}{N}{\displaystyle \sum _{j=1}^N}{\left({v_x}_k^j-{{\widehat{v}}_x}_k^j\right)}^2+{\left({v_y}_k^j-{{\widehat{v}}_y}_k^j\right)}^2}\end{array}$$

(44)

where N is the number of Monte Carlo simulations, j is the jth Monte Carlo simulation, and k is the simulation time. The values $\left(x_k^j,y_k^j\right)$ and $\left({v_x}_k^j,{v_y}_k^j\right)$ represent the true position and velocity of the maneuvering target, respectively; $\left({\widehat{x}}_k^j,{\widehat{y}}_k^j\right)$ and $\left({{\widehat{v}}_x}_k^j,{{\widehat{v}}_y}_k^j\right)$ represent the estimated position and velocity of the filter, respectively.

For filters, consistency is equally as crucial as accuracy because it can help us gauge the algorithm’s robustness. The average normalized estimation error square (ANEES), the consistency analysis index, can be calculated as follows [7]:

$$ANEE{S}_k=\frac{1}{nN}\sum _{i=1}^N{\tilde{x}}_k^T{P}_k^{-1}{\tilde{x}}_k$$

(45)

where ${\tilde{x}}_k$ and $P_k$ are the state estimation error and covariance matrix at time k, respectively. If $ANEE{S}_k\in \left[l_b,u_b\right]$, the filter is considered consistent, where $l_b$ and $u_b$ are the lower and upper bounds of the acceptance interval, respectively. If $ANEE{S}_k<l_b$, then the filter is regarded as “pessimistic” (lacking confidence) because the posterior error covariance is very high compared with the true value; if $ANEE{S}_k>u_b$, then the filter is considered “optimistic” (overconfident), and the covariance $P_{k|k}$ is too small. Different probability intervals can change the upper and lower bounds, but the overall acceptability range is close to 1.

5.1. The System Model Does Not Match

A single kinematic model usually cannot express the motion characteristics of maneuvering targets, which may lead to the loss of targets in the target tracking system when maneuvering orbit changes occur. The interacting-multiple model is a widely used and effective method of dealing with model mismatch issues. As a result, we incorporated IMMCKF as one of the comparison algorithms in this experiment. In this case, the maneuvering target maintains an alternating motion state in 0~250 s, and the periods of uniform motion are 0~70 s, 121~145 s, 161~200 s, and 226~250 s; the periods of maneuver turning are 71~120 s, 146~160 s, and 201~225 s. Among these, the uniform motion section satisfies Equation (43), the maneuvering turning section’s motion equation satisfies Equation (46), and the maneuvering turning rate $\omega =5°/\mathrm{s}$. Figure 2 displays the maneuvering target’s trajectory.

$$F_{ct}=\left[\begin{array}{ccccc}1& \frac{\sin \omega T}{\omega }& 0& \frac{\cos \omega T-1}{\omega }& 0\\ 0& \cos \omega T& 0& -\sin \omega T& 0\\ 0& \frac{1-\cos \omega T}{\omega }& 1& \frac{\sin \omega T}{\omega }& 0\\ 0& \sin \omega T& 0& \cos \omega T& 0\\ 0& 0& 0& 0& 1\end{array}\right]$$

(46)

The RMSE of the position estimate and velocity estimation for various filtering algorithms in Case 5.1 is shown in Figure 3. The filtering algorithm’s model and the initial stage motion are uniform. As a result, these algorithms can produce accurate tracking results. The filtering algorithm’s model and the target’s actual motion diverge significantly when the target maneuvers. The CKF is entirely divergent at this moment and cannot trace the target. The HCKF, based on the CKF, reduces the impact of outliers by a particular function. Although the filter’s robustness is much better than the CKF’s, it still cannot perform effective tracking during the model mismatch stage. The improved SHCKF outperforms the algorithm above at tracking performance. The estimation accuracy of the system noise is nevertheless decreased to prevent algorithm divergence. As a result, when the target changes its movement suddenly, it cannot be rectified in time, and accurate tracking requires some time. By introducing an adaptive fading factor, the ARCKF fully utilizes the useful information in the residual sequence and increases the algorithm’s resistance to uncertainty modeling. The RMSE of the speed estimation only slightly increases when the state changes abruptly, and the RMSE of the position estimation almost wholly ignores the model mismatch. The ARCKF can achieve tracking accuracy similar to the IMMCKF. However, the fundamental model and the Markov transition matrix significantly impact IMMCKF performance. The IMMCKF’s filtering accuracy will be significantly decreased if the trajectory of the maneuvering target does not follow the fundamental model.

5.2. Non-Gaussian Measurement Noise

In general, the measurement system is susceptible to abnormal noise pollution. The noise distribution cannot meet the Gaussian distribution and usually presents a heavy-tailed distribution. Contaminated Gaussian distribution noise can be expressed as [22]:

$$w_k~\left(1-\delta \right)N\left(0,R_n\right)+\delta \left(0,R_p\right)$$

(47)

where $\delta \in \left[\mathrm{0,1}\right]$ denotes the proportion of contaminated noise, $R_n$ denotes the standard measurement noise error covariance matrix, and $R_p$ denotes the contaminated noise covariance matrix, which can be any symmetric distribution. If $R_p$ is also a Gaussian distribution with a large standard deviation, the actual noise is also called Gaussian mixed noise, in which the likelihood of abnormal noise grows with an increase in the weight $\delta$ of the polluted noise $R_p$, and the variance of $R_p$ determines the magnitude of the observation deviation produced by the abnormal noise. In this experiment, we assume that $R_u=50R_n$, and the pollution ratio is 0.2.

Figure 4 shows the RMSE of the position and velocity estimation of different filtering algorithms in Case 5.2. In this case, the HCKF has better filtering accuracy than CKF, converging more quickly during the early filtering stage than CKF; the filtering accuracy of the SHCKF and ARCKF is significantly higher than that of the above two algorithms, but the convergence speed of the SHCKF in the initial stage of filtering is the slowest among the algorithms, which is most obvious in the RMSE of position estimation. The ARCKF retains the advantage that HCKF has the fastest convergence speed in the initial filtering stage. The ARCKF can update the measurement noise covariance in real time, which improves the filtering accuracy of HCKF, and the tracking effect is better than the other three filtering algorithms. Figure 5 shows the estimation consistency of the various filtering methods in Case 5.2. Only the consistency of the CKF is severely impacted by non-Gaussian noise. The other algorithms offer good consistency, but only the ARCKF can do so at the initial filtering stage.

5.3. Variable Noise Covariance

When the measurement system tracks the target, sometimes the signal is unstable, which will cause the noise covariance error to change. In this case, we mainly discuss when the noise covariance matrix suddenly increases significantly. The measurement noise is assumed to satisfy the Gaussian distribution, but it will sometimes change:

$$R_{k+1}=\left\{\begin{array}{cc}{R}_n& t<50s\\ {10R}_n& t\ge 50s\end{array}\right.$$

(48)

where $R_n$ represents the standard measurement noise error covariance matrix, and t represents the simulation time.

Figure 6 shows the RMSE of the position and velocity estimation of the different filtering algorithms in Case 5.3. When the noise covariance matrix becomes larger, the robustness of the four algorithms decreases to a certain extent. However, the SHCKF and ARCKF use the process of adaptively updating the measurement noise covariance. Therefore, these two algorithms have the best resistance to abnormal noise interference after the noise covariance matrix changes. However, when each state changes, the SHCKF always takes a long time to make the algorithm converge. In the face of frequent and complex state changes, the filtering algorithm’s divergence probability will be greatly improved. While using Sage–Husa estimation to update $R_k$, the ARCKF introduces the MD method to correct the $R_k$ matrix that is easy to diverge, improving the algorithm’s robustness. Figure 7 shows the estimation consistency of the various filtering methods in Case 5.3. The first stage’s actual noise covariance is identical to that of $R_0$. These algorithms are relatively consistent, while the CKF and ARCKF can keep consistency more quickly. In the second stage, the measurement noise increases, seriously impairing the consistency of the CKF. Over time, the SHCKF can return to good consistency, but only the HCKF and ARCKF can consistently maintain good consistency.

6. Conclusions

An adaptive robust cubature Kalman filter is proposed in this study to address the problems of system model uncertainty and abnormal measurement noise in maneuvering target tracking. Based on the traditional HCKF, the algorithm introduces a simplified adaptive fading factor to solve the problem of system model uncertainty. Sage–Husa estimation is used to update $R_k$ in real time to significantly increase the estimation accuracy of the HCKF for the time-varying measurement noise covariance matrix. When an observation value is heavily tailed or anomalous, the MD method includes a scaling factor to lower the observation weight, further enhancing the algorithm’s robustness. The simulation results demonstrate that the proposed ARCKF method has good filtering and estimation performance in tracking moving targets. It can produce results for the model mismatch problem comparable to those of the IMM method, and it also performs better than the CKF, HCKF, and SHCKF algorithms in terms of robustness to non-Gaussian measurement noise.