A Robust Kalman Filter Based on the Pearson Type VII-Inverse Wishart Distribution: Symmetrical Treatment of Time-Varying Measurement Bias and Heavy-Tailed Noise

Abstract

This paper introduces a novel robust Kalman filter designed to leverage symmetrical properties within the Pearson Type VII-Inverse Wishart (PVIW) distribution, enhancing state estimation accuracy in the presence of time-varying biases and non-stationary heavy-tailed (NSHT) noise. The filter includes a shape parameter from the normal distribution and an extra variable from the Gamma distribution, which are used to symmetrically adjust the average and variation measures of the data to fit better under difficult noise conditions. To deal with unknown noise that changes over time, the filter uses the Inverse Wishart distribution to model and estimate the scale matrix deviations, making it easier to adapt to changes. The filter also uses a technique called Variational Bayesian to estimate both the state and the parameters at the same time. The results from simulations show that this new filter greatly improves the accuracy and strength of the estimation compared to the usual Kalman filters that assume a normal distribution, especially when there is non-stationary heavy-tailed noise. The main objective is to improve estimation in signal processing and control systems where heavy-tailed noise is prevalent.

1. Introduction

The Kalman filter is a well-established method for recursive state estimation and is widely utilized in information processing domains [1,2,3]. However, its performance heavily depends on Gaussian noise assumptions and precise system models. In practical applications, heavy-tailed noise and unknown biases often affect measurement data, significantly degrading the traditional Kalman filters’ performance [4,5,6]. Particularly in the case of heavy-tailed noise, filtering methods utilizing Gaussian assumptions produce substantial estimation errors [7,8,9,10,11]. Several improved Kalman filters (KFs) have been proposed to help with measurement outliers [12,13,14]. However, the performance of these filters could be better when faced with measurement noise that has pronounced heavy-tailed characteristics due to their lack of modeling for such noise. As a result, several Kalman filters have been introduced that incorporate modeling based on the Student’s t distribution lately.

To address this challenge, researchers have developed a family of filtering methods founded on the Student’s t distribution [15,16,17,18]. Compared to the distribution known as Gaussian, the Student’s t distribution features heavier tails, making it more effective at modeling heavy-tailed noise in natural systems [19]. These methods approximate the posterior probability density function using the Student’s t distribution and recursively estimate joint state vectors and distribution parameters. However, in dynamic environments, the statistical properties of measurement noise often exhibit time-varying characteristics, posing significant challenges to Student’s t-distribution filters with fixed parameters. To address this issue, ref. [20] proposed a robust interactive multiple-model filter, innovatively introducing Gamma and inverse Wishart (IW) distributions to the control degree of freedom (DOF) and the scale matrices, respectively, achieving adaptive estimation of time-varying noise characteristics. Building upon this, ref. [21] proposed the RSTKF method, effectively addressing heavy-tailed measurement noise in linear systems by combining measurement differences, state augmentation strategies, and variational Bayesian methods.

Although Student’s t-distribution-based approaches have made substantial advancements for controlling heavy-tailed noise, their set DOF parameters still limit the filter’s adaptability to varying heavy-tailed characteristics. This limitation has motivated academics to examine the Pearson Type VII distribution, which offers an improved and adjustable parametric structure. By removing the need for “equal DOF parameters”, the Pearson VII distribution provides more flexibility in adapting to random heavy-tailed noise than the Student’s t distribution. Making use of this benefit, ref. [22] suggested a Kalman filter based on the Gaussian–Pearson Type VII (GPTV) mixed distribution, achieving adaptive estimation of unknown time-varying measurement noise through Beta-Bernoulli distribution judgment factors and dual Gamma distribution parameter modeling. Ref. [23] used variational Bayesian inference to jointly estimate states, noise covariance, and distribution parameters, further developing the sliding window-based PTVSWAKF technique.

Nevertheless, these approaches mainly depend on the assumption of zero-mean measurement noise, which is often unrealistic in practical applications. These approaches’ estimating accuracy drastically declines when systems display unpredictable and time-varying measurement biases. Researchers have proposed various solutions to this issue: Ref. [24] developed a filtering algorithm based on incremental measurement equations, while ref. [25] designed an online state estimation scheme with robust mechanisms to detect and mitigate measurement biases through recursive Bayesian inference. For handling unknown measurement biases, ref. [26] developed a strong filtering algorithm employing a Gaussian distribution, while ref. [27] developed the STIWRKF technique centered on the Student’s t Inverse Wishart distribution, jointly estimating system states, bias parameters, and noise covariance through VB methods. Ref. [28] proposed the PRKF method based on Student’s t-distributed process noise, achieving adaptive real-time estimation of time-varying process biases. Despite their excellent performance in specific scenarios, Student’s t-distribution-based filtering methods still struggle to fully adapt to varying heavy-tailed characteristics due to fixed DOF parameters. At the same time, their gradual mean bias iteration strategy fails to respond quickly to dramatic noise changes in complex scenarios.

Most algorithms that rely on fixed degree of freedom (DOF) parameters face challenges when dealing with time-varing heavy-tailed measurement noise. This is because imprecise a priori parameter selection for the Student’s t distribution may lead to incorrect noise modeling, affecting the algorithm’s precision. In contrast, the Pearson Type VII distribution provides more freedom in parameter selection without the restriction that ’two DOF parameters must be equal’, improving its adaptability to random heavy-tailed noise. The robustness of the Pearson Type VII distribution in modeling heavy-tailed noise makes it an essential tool for constructing robust Kalman filtering methods [29,30]. Unfortunately, these existing algorithms suffer from varying degrees of degradation in estimating heavy-tailed measurement noise (HTMN) with not clear and dynamic biases. Most of these algorithms rely on the assumption that measurement noise is zero-centered. Although a few bias estimation methods have been presented in the literature [31,32], their algorithms’ accuracy still needs to be raised.

This paper proposes a novel Pearson Type VII-Inverse Wishart (PVIW) distribution to accurately model HTMN with unknown and time-varying biases to overcome these limitations. This new distribution is able to accommodate not only symmetric noises, such as Gaussian distributions, but also a variety of time-varying heavy-tailed noise environments, demonstrating its adaptability and robustness under different noise conditions. We present an innovative hierarchical representation of the PVIW distribution that combines Gaussian, inverse Wishart, and normal Gamma distributions. Additionally, we introduce a robust Kalman filter developed through the joint learning of state vectors and auxiliary random variables, utilizing variational Bayesian (VB) methods. This approach enables dynamic adjustment of the Pearson Type VII distribution’s mean and scale matrix, effectively adapting to non-stationary heavy-tailed (NSHT) noise characteristics. Simulation experiments validate the superiority and effectiveness of PVIW.

The primary contributions of this study are as follows:

• To address the issues caused by NSHT noise, we suggest a novel distribution dubbed the PVIW distribution. This distribution enhances the Pearson Type VII distribution by incorporating two auxiliary parameters that change its form and scale. Our findings reveal that the PVIW distribution outperforms typical Kalman filters and the resilient Student’s t distribution Kalman filter (RSTKF) in effectively modeling NSHT noise;
• We optimise the Pearson VII distribution by a Gamma distribution to enhance the state estimation accuracy and robustness of the filter under time-varying bias and non-stationary heavy-tailed noise, and propose an innovative time-varying bias estimation method;
• We develop an adaptive Kalman filter (PVIW) optimized for NSHT measurement noise. The filter performs successfully in simulations and experiments, confirming its usefulness and superiority in actual applications.

2. Preliminaries

2.1. State-Space Model

The general filtering problems in the environments with outlier interference are usually based on the following state-space model:

$$\left\{\begin{array}{c}{X}_k=f\left(X_{k-1}\right)+{\omega }_k\hfill \\ Y_k=h\left(X_k\right)+v_k\hfill \end{array}\right.,$$

(1)

where the vector $X_k\in {\mathbb{R}}^n$ represents the state of the system at time k, while the vector $Y_k\in {\mathbb{R}}^m$ denotes the measurement data at that time. The vector ${\omega }_k\in {\mathbb{R}}^n$ corresponds to the Gaussian process noise, and the vector $v_k\in {\mathbb{R}}^m$ stands for the NSHT measurement noise. The function f represents the dynamic model, while the function h denotes the measurement model. With these equations, we can predict the state of the system and update our estimates based on new measurements, which is the core of how the Kalman filter works.

2.2. PVIW Distribution

The impact of the mysterious HTMN poses a significant challenge to the performance of traditional filtering methods, as shown in Figure 1. The graph vividly demonstrates the complex interaction between measurement noise and outliers, with their magnitudes varying from very small to large. Figure 2 presents the contrast of the probability density functions (PDFs) for HTMN, the Gaussian distribution, the Pearson Type VII distribution, and the Student’s t distribution. The Pearson Type VII distribution is distinguished by its extreme tail heaviness in the probability distribution, which provides it with an outstanding ability to decrease the effect of heavy-tailed noise. However, the heavy-tailed characteristics of HTMN often need to be clarified in practical situations, potentially resulting in inadequate representation by the conventional Student’s t distribution. Therefore, there is an urgent requirement for a more adaptable and flexible approach to address this uncertainty, enhancing the filter’s efficiency and robustness in the face of complex noise environments.

In real-world scenarios, measurement noise often shows irregular and heavy-tailed traits, mainly due to time-varying outliers. To tackle this issue, this research proposes a novel PVIW distribution model to deal with the dynamic nature of anomalous measurement noise. The PVIW distribution model integrates a Pearson Type VII distribution with an IW distribution, resulting in a composite model that accurately captures the statistical features of noise. The mathematical expression of the PVIW distribution is described in the following form:

$$PVIW(x;\epsilon ,\Sigma ,\delta ,\varphi ,\sigma ,u,U)=\int PV(x;\epsilon +\delta ,\Sigma +R,\varphi ,\sigma )IW(R;u,U)dR,$$

(2)

in which x is a random variable, $\epsilon$ is the mean, $\Sigma$ is the scale matrix, $\varphi$ and $\sigma$ are the first and second DOF parameters, and $\delta$ is the term that adjusts the mean. R is the term that adjusts the scale matrix $\Sigma$. PV follows a Pearson Type VII distribution and u is the DOF parameter of R. Besides $\varphi$ following a Gamma distribution, the mean adjustment parameters $\delta$ are assumed to follow a Gaussian distribution. They can be represented in the following form:

$$\begin{array}{cc}\hfill p\left(\varphi \right)& =G(\varphi ;a,b),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill p\left(\sigma \right)& =G(\sigma ;c,d),\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill p\left(\delta \right)& =N(\delta ;g,h).\hfill \end{array}$$

(5)

Here, $G(x;a,b)$ follows a Gamma distribution, where a is the shape parameter and b is the scale parameter, and x is a random variable. $N(\delta ;g,h)$ follows a normal distribution, where g is the mean and h is the variance. The Pearson Type VII distribution is characterized as the convolution of a Gaussian distribution and the Gamma distribution [30]. It is described by the following form:

$$PV(x;\epsilon ,\Sigma ,\varphi ,\sigma )=\int N(x;\epsilon ,\Sigma /\kappa )G(\kappa ;\varphi ,\sigma )d\kappa .$$

(6)

We must select its prior distribution to adaptively estimate and update the measurement model’s unknown scale matrix $\Sigma$. According to [33], the IW distribution is suitable for $\Sigma$ as a conjugate prior:

$$\begin{array}{cc}\hfill p(\Sigma )& =IW(\Sigma ;u,U)\hfill \\ \hfill & ={\displaystyle \frac{{\left|U\right|}^{u/2}{|\Sigma |}^{-(u+m+1)/2}exp\left\{\mathrm{tr}\left(U{\Sigma }^{-1}\right)/2\right\}}{2^{mu/2}{\Gamma }_m(u/2)}},\hfill \end{array}$$

(7)

where u is the DOF parameter, and U is the scale matrix. When $u>m+1$, $E\left[R^{-1}\right]=(u-m-1)U^{-1}$, where m is the dimension of U.

According to the additivity of the normal distribution, the following can be calculated:

$$N(x;\delta +\mu ,\Sigma +R)=\int N(x;\eta ,R)N(\eta ;\delta +\mu ,\Sigma )d\eta .$$

(8)

Using Equations (2), (6) and (8), it can be derived that:

$$\begin{array}{cc}\hfill PVIW(x;\epsilon ,\Sigma ,\delta ,\varphi ,\sigma ,u,T)& =\int \int \int N(x;\eta ,R/\kappa )N(\eta ;\epsilon +\delta ,\Sigma /\kappa )\hfill \\ \hfill & \times G(\kappa ;\varphi ,\sigma )IW(R;u,U)d\eta d\kappa dR.\hfill \end{array}$$

(9)

Utilizing Equation (9), the proposed PVIW distribution can be expressed as the Gaussian hierarchical structural model, as shown in the equation:

$$\left\{\begin{array}{c}p\left(x\right|\eta ,R,\kappa )=N(x;\eta ,R/\kappa )\hfill \\ p(\eta ,\kappa |\delta ,\varphi ,\sigma )=N(\eta ;\epsilon +\delta ,\Sigma /\kappa \left)G\right(\kappa ;\varphi ,\sigma )\hfill \\ p\left(\varphi \right)=G(\varphi ;a,b)\hfill \\ p\left(\sigma \right)=G(\sigma ;c,d)\hfill \\ p\left(\delta \right)=N(g,h)\hfill \\ p\left(R\right)=IW(R;u,U)\hfill \end{array}.\right.$$

(10)

3. State Estimation for Uncertain Heavy-Tailed Noise

3.1. Conjugate Prior Distributions and Time Update

The one-step estimation of the state equation could be represented as a Gaussian distribution of the following form:

$$p\left(X_k\right|Y_{1:k-1})=N(X_k;{\widehat{X}}_{k|k-1},P_{k|k-1}),$$

(11)

where ${\widehat{X}}_{k|k-1}$ is the one-step predicted value at time $k-1$, and $P_{k|k-1}$ is the appropriate covariance matrix. It can be determined using the following variant of the Kalman filter equation:

$$X_{k|k-1}=f\left({\widehat{X}}_{k-1|k-1}\right),$$

(12)

$$P_{k|k-1}=F_{k-1}{P}_{k-1|k-1}{F}_{k-1}^T+Q_k,$$

(13)

where ${\widehat{X}}_{k-1|k-1}$ and $P_{k-1|k-1}$ are the state estimation vector and the corresponding covariance matrix at time $k-1$, respectively. Matrix $Q_k$ denotes the covariance of the state noise, while $F_{k-1}$ represents the Jacobian matrix derived from function f. The measurement noise is characterized using the PVIW distribution:

$$p\left({\upsilon }_k\right)=PVIW({\upsilon }_k;0,{\Sigma }_k,{\delta }_k,{\varphi }_k,{\sigma }_k,u_k,U_k),$$

(14)

where ${\Sigma }_k$ and $U_k$ are the scale matrices, and ${\varphi }_k$ and ${\sigma }_k$ are the two-degree-of-freedom parameters, and ${\delta }_k$ is the third-degree-of-freedom parameter. Using Equations (1) and (14), the likelihood probability density function $p\left(Y_k\right|X_k)$ can be written as follows:

$$p\left(Y_k\right|X_k)=PVIW(Y_k;h\left(X_k\right),{\Sigma }_k,{\varphi }_k,{\sigma }_k,{\delta }_k,u_k,U_k).$$

(15)

Utilizing Equations (10) and (15), $p(Y_k\mid X_k)$ can be expressed as the Gaussian hierarchical structural model, as shown in the equation:

$$\left\{\begin{array}{c}p\left(Y_k\right|X_k,R_k,{\eta }_k,{\kappa }_k)=N(Y_k;f\left(x_k\right)+{\eta }_k,R_k/{\kappa }_k)\hfill \\ p({\eta }_k,{\kappa }_k|{\delta }_k,{\varphi }_k,{\sigma }_k)=N({\eta }_k;{\delta }_k,{\Sigma }_k/{\kappa }_k)G({\kappa }_k;{\varphi }_k,{\sigma }_k)\hfill \\ p\left({\varphi }_k\right)=G({\varphi }_k;a_k,b_k)\hfill \\ p\left({\sigma }_k\right)=G({\sigma }_k;c_k,d_k)\hfill \\ p\left({\delta }_k\right)=N({\delta }_k;{\delta }_0,y{I}_m)\hfill \\ p\left(R_k\right)=IW(R_k;u_k,U_k)\hfill \end{array},\right.$$

(16)

where ${\delta }_0$ is the initial vector and $y{I}_m$ is the covariance matrix, and y is the restricted random variable $y{I}_m$. ${\delta }_k$ serves as the adjustment factor for the variable $h\left(x_k\right)$.

Figure 3 displays the suggested algorithm’s Gaussian hierarchical structural model. Here, variational Bayesian inference is used to estimate the parameters and states, where ${\lambda }_k$, $\kappa$, ${\varphi }_k$, ${\sigma }_k$ and $R_k$ are the parameters to be estimated.

3.2. Approximate Posterior Based on Variational Bayesian

In nonlinear systems, the variational Bayesian approach has become a popular approximation strategy. This is because the joint posterior probability density function and the parameter to be estimated often make it challenging to obtain an exact analytical solution. The core idea of this approach is to identify a parametric distribution q that closely approximates the true posterior distribution. Subsequently, this approximation is improved by minimizing the Kullback–Leibler divergence (KL divergence), or relative entropy, between the two distributions. The relative entropy, denoted as $D_{KL}(q\parallel p)$, measures the loss of information when the actual distribution is p and the assumed distribution is q. Its definition formula is as follows:

$$D_{KL}(q\parallel p)=\sum _{x}q\left(x\right)log{\displaystyle \frac{q\left(x\right)}{p\left(x\right)}}.$$

(17)

Since it is more challenging to compute $D_{KL}(q\parallel p)$ directly, variational Bayesian methods usually minimize $D_{KL}(q\parallel p)$ indirectly by maximizing the Evidence Lower Bound (ELBO). The ELBO can be expressed as:

$$\mathrm{ELBO}\left(q\right)={\mathbb{E}}_q[logp(X,Z)-D_{KL}(q\parallel p\left(Z\right))],$$

(18)

where X represents the observed data, Z represents the latent variables, $p(X,Z)$ is the joint probability distribution, and $p\left(Z\right)$ is the prior distribution. The first term of the ELBO, ${\mathbb{E}}_q[logp(X,Z)]$, represents the expectation of the log-likelihood under the approximate distribution q. The second term, $D_{KL}(q\parallel p\left(Z\right))$, represents the relative entropy between the approximate distribution and the prior distribution. Maximizing the ELBO makes it possible to keep a good fit for the observed data while making the approximate distribution q close to the actual posterior distribution p.

Under uncertain heavy-tailed noise conditions, the joint posterior distribution of the system’s state, the two DOF, and the auxiliary parameters can be expressed as:

$$p(X_k,{\eta }_k,{\kappa }_k,{\varphi }_k,{\sigma }_k,{\delta }_k,R_k|Y_{1:k})=p\left({\Theta }_k\right|Y_{1:k}),$$

(19)

where ${\Theta }_k\triangleq \{X_k,{\eta }_k,{\kappa }_k,{\varphi }_k,{\sigma }_k,{\delta }_k,R_k\}$.

The approximate distribution of $p\left({\Theta }_k\right|Y_{1:k})$ can be written as:

$$p\left({\Theta }_k\right|Y_{1:k})\approx q\left(X_k\right)q({\eta }_k,{\kappa }_k)q\left({\varphi }_k\right)q\left({\sigma }_k\right)q\left({\delta }_k\right)q\left(R_k\right),$$

(20)

in which $q\left(X_k\right)$, $q({\eta }_k,{\kappa }_k)$, $q\left({\varphi }_k\right)$, $q\left({\sigma }_k\right)$, $q\left({\delta }_k\right)$, and $q\left(R_k\right)$ are the approximate posterior distributions of the variables $X_k$, ${\eta }_k$, ${\kappa }_k$, ${\varphi }_k$, ${\sigma }_k$, ${\delta }_k$, $R_k$. The problem is equivalent to

$$logq\left(\theta \right)={\mathbb{E}}_{{\Theta }_k(-\theta )}[logp({\Theta }_k,Y_{1:k})]+C_{\theta },$$

(21)

where $\theta$ is any element in ${\Theta }_k\triangleq \{X_k,{\eta }_k,{\kappa }_k,{\varphi }_k,{\sigma }_k,{\delta }_k,R_k\}$, $C_{\theta }$ is a constant related to $\theta$, and ${\Theta }_k(-\theta )$ represents all variables in ${\Theta }_k$ except $\theta$.

The following is an expression for the joint posterior probability density function:

$$\begin{array}{cc}\hfill p({\Theta }_k,Y_{1:k})& =p\left(Y_k\right|X_k,{\eta }_k,R_k,{\kappa }_k)p\left(X_k\right|Y_{1:k-1})\hfill \\ \hfill & \times q\left({\varphi }_k\right)p\left({\delta }_k\right)p\left(R_k\right)p\left(Y_{1:k-1}\right)\hfill \\ \hfill & \times p({\eta }_k,{\kappa }_k|{\delta }_k,{\varphi }_k,{\sigma }_k).\hfill \end{array}$$

(22)

Using Equations (11), (16) and (22), we can obtain the following:

$$\begin{array}{cc}\hfill p({\Theta }_k,Y_{1:k})& =N(Y_k;f\left(x_k\right)+{\eta }_k,R_k/{\kappa }_k)N(X_k;X_{k|k-1},P_{k|k-1})\hfill \\ \hfill & \times N({\eta }_k;{\delta }_k,{\Sigma }_k/{\kappa }_k)G({\kappa }_k;{\varphi }_k,{\sigma }_k)G({\varphi }_k;a_k,b_k)\hfill \\ \hfill & \times G({\sigma }_k;c_k,d_k)N({\delta }_k;{\delta }_0,y{I}_m)IW(R_k;u_k,U_k)p\left(Y_{1:k-1}\right).\hfill \end{array}$$

(23)

Since Equation (23) has no closed-form analytical solution, it cannot be solved directly. The fixed-point iteration approach can be used to tackle this problem. During the iteration process, when updating any variational parameters, the values of PDFs obtained from the previous iteration must be utilized to ensure the stability and convergence of the iteration. This method continuously adjusts the parameters using stepwise approximation until a solution satisfies the condition.

Proposition 1.

Let $\theta =X_k$. Using Equations (21) and (23), the state estimation $q^{(i+1)}\left(X_k\right)$ follows a normal distribution:

$$q^{(i+1)}\left(X_k\right)=N(X_k;{\widehat{X}}_{k|k}^{(i+1)},P_{k|k}^{(i+1)}),$$

(24)

where ${\widehat{X}}_{k|k}^{(i+1)}$ is the updated estimate and $P_{k|k}^{(i+1)}$ is the updated covariance matrix given by:

$${\widehat{X}}_{k|k}^{(i+1)}={\widehat{X}}_{k|k-1}+K_k^{(i+1)}(Y_k-{\tilde{Y}}_k^i),$$

(25)

$$P_{k|k}^{(i+1)}=(I_m-K_k^{(i+1)}{H}_k)P_{k|k-1},$$

(26)

$$K_k^{(i+1)}=P_{k|k-1}{H}_k^T{(H_k{P}_{k|k-1}{H}_k^T+{\tilde{R}}_k^{\left(i\right)})}^{-1}.$$

(27)

Here, $K_k^{(i+1)}$ is the Kalman gain matrix, $H_k$ is the Jacobian matrix of h, and ${\tilde{Y}}_k$ and ${\tilde{R}}_k$ are intermediate variables defined by:

$${\tilde{Y}}_k^{\left(i\right)}=h\left(X_k\right)+E^{\left(i\right)}\left[{\eta }_k\right],$$

(28)

$${\tilde{R}}_k^{\left(i\right)}={\displaystyle \frac{{E^{\left(i\right)}\left[R_k^{-1}\right]}^{-1}}{E^{\left(i\right)}\left[{\kappa }_k\right]}}.$$

(29)

Proof. 

See Appendix A. □

Proposition 2.

Let $\theta =\left\{{\eta }_k,{\kappa }_k\right\}$, using Equations (21) and (23), yields the result that $q^{(i+1)}({\eta }_k,{\kappa }_k)$ satisfies the normal Gamma distribution, which can be expressed as:

$$\begin{array}{cc}\hfill q^{(i+1)}({\eta }_k,{\kappa }_k)& =N({\eta }_k;{\widehat{\mu }}_{k|k}^{(i+1)},{\displaystyle \frac{{\Sigma }_{k|k}^{(i+1)}}{E^{\left(i\right)}\left[{\kappa }_k\right]}})\hfill \\ \hfill & \times G({\kappa }_k,A^{(i+1)},B^{(i+1)}).\hfill \end{array}$$

(30)

The parameters are updated as follows:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\hfill W_k^{(i+1)}& ={\Sigma }_{k|k-1}{({\Sigma }_{k|k-1}+{\widehat{R}}_k^{(i+1)})}^{-1}\hfill \\ \hfill {\widehat{\mu }}_{k|k}^{(i+1)}& =E^i\left[{\delta }_k\right]+W_k^{(i+1)}(y_k-h\left({\widehat{x}}_{k|k}^{(i+1)}\right)-E^i\left[{\delta }_k\right])\hfill \\ \hfill {\Sigma }_{k|k}^{(i+1)}& =(I_n-W_k^{(i+1)}){\Sigma }_{k|k-1}\hfill \\ \hfill A_k^{(i+1)}& =m/2+E^{\left(i\right)}\left[{\varphi }_k\right]\hfill \\ \hfill B_k^{(i+1)}& =E\left[{\sigma }_k\right]+0.5\{\mathrm{tr}\left(H_k{P}_{k|k}^{(i+1)}{H}_k^T{\left({\widehat{R}}_k^{(i+1)}\right)}^{-1}\right)\hfill \\ \hfill & +tr(({\tilde{Y}}_k^{i+1}-E^{(i+1)}\left[{\delta }_k\right]){({\tilde{Y}}_k^{i+1}-E^{(i+1)}\left[{\delta }_k\right])}^T{({\Sigma }_{k|k-1}+{\widehat{R}}_k^{(i+1)})}^{-1}\}\hfill \end{array},\right.\end{array}$$

(31)

Here,

$${\tilde{Y}}_k^{(i+1)}=Y_k-h\left({\widehat{x}}_{k|k}^{(i+1)}\right),{\widehat{R}}_k^{(i+1)}={\left(E^{(i+1)}\left[R_k^{-1}\right]\right)}^{-1}.$$

(32)

Proof. 

See Appendix B. □

Proposition 3.

Let $\theta ={\varphi }_k$, using Equations (21) and (23), yields that $q^{(i+1)}\left({\varphi }_k\right)$ satisfies the Gamma distribution, which can be expressed as:

$$q^{(i+1)}\left({\varphi }_k\right)=G({\varphi }_k,a_k^{(i+1)},b_k^{(i+1)}),$$

(33)

where the parameters $a_k^{(i+1)}$ and $b_k^{(i+1)}$ are updated as:

$$\begin{array}{cc}\hfill a_k^{(i+1)}& =a+0.5-E^{\left(i\right)}\left[{\varphi }_k\right],\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill b_k^{(i+1)}& =b-1-E^{(i+1)}[log{\kappa }_k]-E^{\left(i\right)}[log{\sigma }_k].\hfill \end{array}$$

(35)

Let $\theta ={\sigma }_k$, the distribution $q^{(i+1)}\left({\sigma }_k\right)$ also follows a Gamma distribution:

$$q^{(i+1)}\left({\sigma }_k\right)=G({\sigma }_k,c_k^{(i+1)},d_k^{(i+1)}),$$

(36)

with parameters $c_k^{(i+1)}$ and $d_k^{(i+1)}$ updated as:

$$\begin{array}{cc}\hfill c_k^{(i+1)}& =c+E^{\left(i\right)}\left[{\varphi }_k\right],\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill d_k^{(i+1)}& =d+E^{(i+1)}\left[{\kappa }_k\right].\hfill \end{array}$$

(38)

Proof. 

See Appendix C. □

Proposition 4.

Let $\theta ={\delta }_k$, using Equations (21) and (23), yields that $q^{(i+1)}\left({\delta }_k\right)$ satisfies the Gaussian distribution, which can be expressed as:

$$q^{(i+1)}\left({\delta }_k\right)=N({\delta }_k;{\delta }_k^{(i+1)},P_{\delta }^{(i+1)}),$$

(39)

where ${\delta }_k^{(i+1)}$ and $P_{\delta }^{(i+1)}$ can be represented as:

$${\delta }_k^{(i+1)}={\delta }_0+K_{\delta }^{(i+1)}(E^{\left(i\right)}\left[{\eta }_k\right]-{\delta }_0),$$

(40)

$$P_{\delta }^{(i+1)}=(I_m-K_{\delta }^{(i+1)})y.$$

(41)

where $K_{\delta }^{(i+1)}$ can be expressed as:

$$\begin{array}{cc}\hfill K_{\delta }^{(i+1)}& =y{[y{I}_m+{\tilde{R}}_k^{\left(i\right)}]}^{-1},\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill {\tilde{R}}_k^{\left(i\right)}& ={\displaystyle \frac{\Sigma }{E^{\left(i\right)}\left[{\kappa }_k\right]}}.\hfill \end{array}$$

(43)

Proof. 

See Appendix D. □

Proposition 5.

Let $\theta =R_k$, using Equations (21) and (23), the distribution of $q^{(i+1)}\left(R_k\right)$ follows an Inverse Wishart distribution:

$$q^{(i+1)}\left(R_k\right)=IW(R_k;u_k^{(i+1)},U_k^{(i+1)}).$$

(44)

The parameters $u_k^{(i+1)}$ and $U_k^{(i+1)}$ are updated as follows:

$$\begin{array}{cc}\hfill {\widehat{u}}_{k|k}^{(i+1)}& ={\widehat{u}}_{k|k-1}+1,\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill {\widehat{U}}_{k|k}^{(i+1)}& ={\widehat{U}}_{k|k-1}+D_k^{(i+1)}.\hfill \end{array}$$

(46)

where $D_k^{(i+1)}$ is given by:

$$\begin{array}{cc}\hfill D_k^{(i+1)}& =E^{(i+1)}\left[\kappa \left((Y_k-h\left(X_k\right)-{\eta }_k){(Y_k-h\left(X_k\right)-{\eta }_k)}^T\right)\right].\hfill \end{array}$$

(47)

Proof. 

See Appendix E. □

In order to realize the fixed-point iteration, the required expectation values need to be calculated with the following formulas.

Using Equation (24), the following equation can be obtained:

$$E^{(i+1)}[Y_k-h\left(X_k\right)]=Y_k-h\left(X_k^{(i+1)}\right),$$

(48)

$$\begin{array}{cc}\hfill E^{(i+1)}\left[(Y_k-h\left(X_k\right)){(Y_k-h\left(X_k\right))}^T\right]& =(Y_k-h({\widehat{X}}_{k|k}^{(i+1)}))(Y_k-h({\widehat{X}}_{k|k}^{(i+1)}){)}^T\hfill \\ \hfill & +H_k{P}_{k|k}^{(i+1)}{H}_k^T,\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill E^{(i+1)}\left[(Y_k-h\left(X_k\right)-{\eta }_k){(Y_k-h\left(X_k\right)-{\eta }_k)}^T\right]& =(Y_k-h({\widehat{X}}_{k|k}^{(i+1)})-E^{(i+1)}({\eta }_k))\hfill \\ \hfill & \times (Y_k-h({\widehat{X}}_{k|k}^{(i+1)})-E^{(i+1)}({\eta }_k){)}^T\hfill \\ \hfill & +H_k{P}_{k|k}^{(i+1)}{H}_k^T.\hfill \end{array}$$

(50)

Using Equation (30), we obtain the expectation associated with ${\eta }_k$, ${\kappa }_k$:

$$\begin{array}{cc}\hfill E^{\left(i\right)}\left[{\eta }_k\right]& ={\widehat{\mu }}_{k|k}^{(i+1)},\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill E^{(i+1)}\left[{\kappa }_k\right]& ={\displaystyle \frac{A^{(i+1)}}{B^{(i+1)}}},\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill E^{(i+1)}[log{\kappa }_k]& =-log\left(B^{(i+1)}\right)+\Psi \left(A^{(i+1)}\right).\hfill \end{array}$$

(53)

in which $\Psi$ is a polygamma function.

Using Equations (33) and (36), we obtain the expectations associated with ${\varphi }_k$ and ${\sigma }_k$:

$$\begin{array}{cc}\hfill E^{\left(i\right)}\left[{\varphi }_k\right]& ={\displaystyle \frac{a_k^{(i+1)}}{b_k^{(i+1)}}},\hfill \end{array}$$

(54)

$$\begin{array}{cc}\hfill E^{\left(i\right)}\left[{\sigma }_k\right]& ={\displaystyle \frac{c_k^{(i+1)}}{d_k^{(i+1)}}},\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill E^{\left(i\right)}[log{\sigma }_k]& =\Psi \left(c_k\right)-log\left(d_k\right).\hfill \end{array}$$

(56)

Using Equation (39), we obtain the expectations associated with ${\delta }_k$:

$$\begin{array}{cc}\hfill E^{(i+1)}\left[{\delta }_k\right]& ={\delta }_k^{(i+1)},\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill E^{(i+1)}\left[{\delta }_k{\delta }_k^T\right]& ={\delta }_k^{(i+1)}{\left({\delta }_k^{(i+1)}\right)}^T+P_{\delta }^{(i+1)}.\hfill \end{array}$$

(58)

Using Equation (44), we obtain the expectations associated with $R_k$:

$$E^{(i+1)}\left[R_k^{-1}\right]={\widehat{u}}_{k|k}^{(i+1)}{\left[{\widehat{U}}_{k|k}^{(i+1)}\right]}^{-1}.$$

(59)

Thus, the robust Kalman filtering method based on the Gamma Pearson Type VII mixture distribution is completed, as shown in

Table 1. One time step of the proposed KF.

Inputs: ${\widehat{x}}_{k-1|k-1},P_{k-1|k-1},F_k,H_k,y_k,Q_k,R_k,{\Sigma }_{k|k},a_k,b_k,{\delta }_k,c_k,d_k,y$

Time Update:

Calculate ${\widehat{X}}_{k|k-1},P_{k|k-1}$ by (12), (13)

1. Initial $E^0\left[{\varphi }_k\right]={\displaystyle \frac{a}{b}},E^0\left[{\sigma }_k\right]={\displaystyle \frac{c}{d}},E^0\left[{\kappa }_k\right]={\displaystyle \frac{E^0\left[{\varphi }_k\right]}{E^0\left[{\sigma }_k\right]}},E^0\left[{\delta }_k\right]={\delta }_0,E^0\left[{\delta }_k{\delta }_k^T\right]={\delta }_0{\delta }_0^T+y{I}_m$

for $i=0:N_I-1$ do

Update ${\widehat{X}}_{k|k}^{(i+1)},P_{k|k}^{(i+1)},K_k^{(i+1)}$ by (24)–(28)

Calculate $E^{(i+1)}[Y_k-h\left(X_k\right)],E^{(i+1)}\left[(Y_k-h\left(X_k\right)){(Y_k-h\left(X_k\right))}^T\right]$ by (48), (49)

Calculate $E^{(i+1)}\left[{\kappa }_k\right]$ by (30), (52)

Calculate $E^{\left(i\right)}\left[{\varphi }_k\right]$ by (33), (54)

Calculate $E^{\left(i\right)}\left[{\sigma }_k\right]$ by (36), (55)

Calculate ${\delta }_k^{(i+1)},P_{\delta }^{(i+1)},K_{\delta }^{(i+1)}$ by (40)–(44)

Calculate $E^{(i+1)}\left[{\delta }_k\right]$ by (39), (57)

Calculate $E^{(i+1)}\left[{\delta }_k{\delta }_k^T\right]$ by (58)

Calculate $E^{(i+1)}\left[(Y_k-h\left(X_k\right)-{\eta }_k){(Y_k-h\left(X_k\right)-{\eta }_k)}^T\right]$ by (50)

end for

${\widehat{X}}_{k|k}={\widehat{X}}_{k|k}^i,P_{k|k}=P_{k|k}^i,{\widehat{\eta }}_{k|k}={\widehat{\eta }}_{k|k}^{(i+1)},{\Sigma }_{k|k}={\Sigma }_{k|k}^{(i+1)},{\widehat{u}}_{k|k}={\widehat{u}}_{k|k}^{(i+1)},{\widehat{U}}_{k|k}={\widehat{U}}_{k|k}^{(i+1)}$

Outputs: ${\widehat{X}}_{k|k},P_{k|k},{\widehat{\eta }}_{k|k},{\Sigma }_{k|k},{\widehat{u}}_{k|k},{\widehat{U}}_{k|k}$

.

3.3. Discussion on the Proposed Algorithm

Remark 1.

The algorithm proposed in this study is suitable for handling nonlinear measurement systems and can effectively cope with the challenges of linear systems. The flexibility and robustness of the algorithm make it ideal for dealing with systems characterized by heavy-tailed noise with offsets. It can be combined with different tracking methods, such as solid tracking, input estimation or multi-modeling algorithms [31,32,33], to estimate rapidly changing measurements and potentially nonlinear characteristics. This paper aims to investigate and improve the algorithm’s ability to handle unknown heavy-tailed noise without addressing the above aspects.

Remark 2.

When ${\varphi }_k$ and ${\sigma }_k$ are equal, and ${\delta }_k$ is a given fixed value, the proposed algorithm becomes [27]. The algorithm presented in this paper has additional DOF, allowing for flexible adjustment of ${\Sigma }_k$. Consequently, the algorithm exhibits greater robustness in handling heavy-tailed noise.

Remark 3.

The filter we have developed is mainly used in the field of signal processing and control systems, where heavy-tailed noise can have a significant impact on system performance. The filter has been validated by simulation and has been applied in maneuvering target tracking. In addition, it can be widely used in areas such as autopilot control and navigation.

Remark 4.

The traditional Kalman filter, designed for Gaussian noise, struggles with heavy-tailed noise. The Pearson Type VII distribution, with its adjustable degrees of freedom, provides greater flexibility in adapting to heavy-tailed noise compared to the Student’s t distribution. The inverse Wishart distribution further enhances adaptability by modeling the scale matrix of the Pearson Type VII distribution. By integrating these distributions, we aim to develop a robust Kalman filter capable of handling non-stationary heavy-tailed noise and time-varying biases. This approach is particularly suited for natural systems and real-world applications where heavy-tailed noise is prevalent, though it may still face challenges with rapidly changing noise characteristics and unknown biases.

4. Simulation

4.1. Linear System Simulation

We created and carried out a number of simulation tests to guarantee the algorithm’s wide applicability and completely validate the suggested algorithm’s efficacy. In the experiments, we compared the newly proposed Kalman filtering algorithm with a variety of existing algorithms, including RSTKF [21], NSMRKF [34], GPTVMAKF [22], RORFS [35], and AETVKF [27], as well as Kalman filtering with known actual error covariance matrix (KFTNCM). Through computer simulations, we comparatively analyze the performance of these filters in handling outliers.

The equations of state and measurement are as follows:

$$\left\{\begin{array}{cc}\hfill X_k& =F_k{X}_{k-1}+W_k\hfill \\ \hfill Z_k& =H_k{X}_k+V_k\hfill \end{array}.\right.$$

(60)

Here, the state vector is $X_k=x_k$. $x_k$ are the position coordinates, the coefficients of the state turn coefficients and measurement equations are $F_k=0.5$, $H_k=1$, and the sampling time is $\Delta T=0.01$ s. The measurement noise exhibits heavy-tailed characteristics, whereas the state noise follows a Gaussian distribution. The outcomes of the filtered solution are simulated and compared using the root mean square error (RMSE) and the average root mean square error (ARMSE) metrics. The metrics are specified as follows:

$$\begin{array}{cc}\hfill RMSEx& =\sqrt{{\displaystyle \frac{1}{M_c}}\sum _{r=1}^{M_c}{(x_k^r-{\widehat{x}}_k^r)}^2},\hfill \end{array}$$

(61)

$$\begin{array}{cc}\hfill ARMSEx& =\sqrt{{\displaystyle \frac{1}{M_{c}N}}\sum _{k=1}^N\sum _{r=1}^{M_c}{(x_k^r-{\widehat{x}}_k^r)}^2}.\hfill \end{array}$$

(62)

where $x_k^r$ and ${\widehat{x}}_k^r$ represent the true and estimated state results at time k for the r-th Monte Carlo simulation, respectively. $M_c=500$ denotes the total number of Monte Carlo simulations, and $N=40$ is the total sampling time. The initial state vector is $X_0=100$, and the covariance matrix is $P_0=1000$. From a Gaussian distribution, such as ${\widehat{x}}_{0|0}\sim N(x_0,P_0)$, the initial state estimate is chosen at random. The status and measurement noises are generated as follows:

$$W_k\sim N(0,Q_0),$$

(63)

$$V_k\sim \left\{\begin{array}{c}N({\mu }_k,R_0)\mathrm{w}.\mathrm{p}.1-p_1\hfill \\ N({\mu }_k,100R_0)\mathrm{w}.\mathrm{p}.p_1\hfill \end{array},\right.$$

(64)

the measurement noise’s real mean vector, which is provided by

$${\mu }_k\sim \left\{\begin{array}{c}100\le k\le N/4\hfill \\ 20N/4<k\le N/2\hfill \\ 30N/2<k\le 3N/4\hfill \\ 103N/4<k\le N\hfill \end{array},\right.$$

(65)

where $w.p.$ stands for “with probability”, $p_1$ represents the probability of outliers occurring within the sampling period. The nominal and the measured state noise covariance matrices are set as $Q_0=R_0=100$. The following are the pertinent parameters of the suggested algorithm: $a_0=25$, $b_0=5$, $c_0=25$, $d_0=5$, ${\delta }_0=0$, $y=20$, ${\Sigma }_{0|0}=100$. The parameters for the comparative algorithms are set up as follows: (1) RSTKF: The DOF parameter is set to 5; (2) NSMRKF: The tuning parameter is set to 5, with the first prior DOF parameter $v_1$ at 0.5 and the second prior DOF parameter $v_2$ at 5, and the mixing probability is 0.85; (3) GPTVMAKF: The Gaussian distribution shape parameter is 5, and the DOF parameter for the Student’s t distribution is also 5, with a mixing probability of 0.85; (4) RORFS: The DOF parameter for the Student’s t distribution is 5, and the scale parameter is set to $1-exp(-4)$; (5) AETVKF: The parameter settings are identical to those of RORFS. By setting the maximum iteration count to $N_m=20$ for both the suggested and current filtering algorithms, the fixed-point iterations are guaranteed to converge. All of the algorithms are built in MATLAB 2019a and run on a PC with a 2.50 GHz Intel Core i7-6500U CPU sourced from Intel Corporation, Santa Clara, CA, USA.

First, we contrasted the suggested algorithm’s estimated accuracy with that of current robust filtering techniques when the outlier probability $p_1$ is 0.1. The measurement noise indicated different heavy-tailed noise characteristics in each phase. Figure 4 displays the RMSE of the state variables for the proposed and current robust filtering technologies.

Table 2. ARMSEs and run time of KFS when $p_1=0.1$.

	ARMSE					Time (ms)
Filters	(1:100)	(101:200)	(201:300)	(301:400)	(1:400)
RORFS	9.82	11.30	11.92	9.92	10.74	0.32
AETVKF	8.95	9.06	9.06	9.21	9.07	0.61
The Proposed Filter	8.12	8.17	8.18	8.24	8.18	1.21
KFTNCM	10.13	15.23	21.21	10.16	14.19	0.06
RSTKF	10.51	16.42	23.77	10.55	15.31	0.91
GPTVMAKF	10.01	13.08	15.61	10.06	12.19	0.92
NSMRKF	10.23	15.62	22.18	10.26	14.57	1.73

lists the corresponding average root mean square error (ARMSE) and single-step execution time for each phase of the state vector. The suggested filter has a lower RMSE and ARMSE than a number of comparable algorithms, as shown in Figure 4 and Table 2. In each phase where the mean changes, the proposed filtering algorithm maintains a stable overall variation in filtering accuracy without significant fluctuations. In terms of the execution time of a single step of the filtering, the proposed filtering algorithm takes longer because it involves the estimation of the scale matrix, the mean, etc. More parameters are involved in the estimation, so the filtering takes longer, and a comparison of the time of each filtering algorithm is given in Table 2.

Secondly, we contrast how quickly the filtering methods converge. Among the filtering algorithms, only the AETVKF algorithm and the proposed filtering algorithm can estimate the noise mean. Here, the two are selected for the comparative study. The findings can be found in Figure 5. As observed in the figure, in the first stage of the simulation (0–10 s), the proposed filtering algorithm converges at 0.5 s. The comparative algorithm AETVKF algorithm converges at about 1.5 s. In the second simulation stage (10–20 s), when the noise mean value jumps from 10 to 20, the proposed algorithm can adapt very quickly and converge at 10.5 s, while the comparative algorithm AETVKF converges at about 12 s. In the third stage of simulation (20–30 s), when the noise mean value jumps from 20 to 30, the proposed filtering algorithm can adapt very quickly and converges with 20.7 s, while the comparison algorithm AETVKF converges at about 23.6 s; in the fourth stage of simulation, the noise mean value jumps from 30 to 10; at this time, the proposed filtering algorithm can adapt very quickly. It converges at the moment of 30.7 s, while the comparison algorithm AETVKF converges at about 23.7 s. The proposed filtering algorithm converges faster and ensures minimum filtering error at all simulation stages.

As seen in Table 2, for $p_1=0.1$, the computational complexity of the proposed filter is similar to that of the algorithms such as RSTKF, GPTVMAKF, etc., but lower than the existing NSMRKF. In the following, the results of state estimation using different iteration numbers are investigated by setting the iteration number as $N_m=1,2,\cdots ,20$, respectively, for simulation, with the outcomes depicted in Figure 6.

Figure 6 clearly demonstrates that the filtering algorithm we proposed can reach a stable state in a relatively short period of time. Not only is its convergence speed rapid, but it also performs better than other methods in terms of state estimate accuracy. This finding confirms that our algorithm can provide higher precision in state estimation when dealing with unknown heavy-tailed noise, thus offering a significant advantage in practical applications.

We further test the robustness of the suggested technique by varying the covariance matrix of measurement noise. The nominal measurement error covariance matrix is set to $0.01R_0$, $0.1R_0$, $R_0$, $2R_0$, $10R_0$, $20R_0$, respectively, and the RMSE of the proposed filtering algorithm is shown in Figure 7. The ARMSE statistics are shown in

Table 3. ARMSEs of the proposed KF with different $R_0$ matrices.

R	0.01 R0	0.1 R0	R0	2 R0	5 R0	10 R0	20 R0
ARMSEs (1:400)	14.03	9.94	8.18	8.23	8.94	9.64	10.29

, from which it can be seen that the ARMSE result of the proposed filtering algorithm is minimized when the nominal error covariance matrix is set to $R_0$. From Table 3, it can be observed that the ARMSE of the suggested filter is minimized when the noise is heavy-tailed. The aforementioned findings suggest that the proposed filter can be well fitted in the case of location of measurement errors.

For simulation, ${\delta }_0$ = 0, 1, 2, 5, 10, 20, 30 are chosen, and Figure 8 displays the RMSEs findings.

According to

Table 4. ARMSEs of the proposed KF with different ${\delta }_0$.

${\mathit{\delta }}_0$	0	1	2	5	10	20	30
ARMSEs (1:400)	8.23	8.19	8.18	8.19	8.20	8.22	8.27

, the simulation results of the suggested filtering method are very similar for a range of ${\delta }_0$ values, indicating that the process may self-adapt to varying initial ${\delta }_0$ values.

Figure 9 displays the RMSE findings for the ${\Sigma }_0=10,20,50,100,200,300,500$ that were chosen for simulation.

From Figure 9, it is evident that the RMSE results remain close as ${\Sigma }_0$ varies, which suggests that the proposed filtering algorithm possesses a certain level of adaptability. This allows it to adjust to varying values of ${\Sigma }_0$ to optimize the outcomes of state estimation.

The simulation results for $y=1,2,5,10,15,20,25$ are displayed in Figure 10.

It can be observed from the figure that the RMSE results are close when y takes different values, indicating that the proposed filtering algorithm has a certain degree of adaptability, enabling it to adapt to different values of y to optimize the state estimation results.

Lastly, we evaluate the accuracy of the suggested filter estimate against the current robust filter in the case of the probability of outlier $p_1=0.4$. In this example, the state and measurement noise have robust heavy-tailed distributions for all periods. Figure 11 shows the stability of the proposed filtering algorithm, which is well adapted to heavy-tailed and general noise. The ARMSEs computations’ outcomes are displayed in

Table 5. ARMSEs OF KFS WHEN $p_1=0.4$.

Filters	ARMSEs (1:400)
RORFS	11.2
AETVKF	11.0
The Proposed Filter	9.20
KFTNCM	13.4
RSTKF	14.3
GPTVMAKF	11.9
NSMRKF	13.7

.

The simulation outcomes demonstrate that the proposed filtering algorithm adapts effectively to heavy-tailed noise across various probability scenarios, exhibiting more stability and robustness than prior algorithms.

4.2. Nonlinear System Simulation

In this section, we initially employ a distinct classical maneuvering target tracking model to assess the performance of the proposed algorithm. To demonstrate the superiority of our algorithm in terms of precision evaluation, we have selected several commonly used filtering algorithms for comparison, such as RSTKF [21], NSMRKF [34], the Extended Kalman Filter (EKF) and so on, through simulation. The state transition equation for the system is given by:

$$x_k=\left[\begin{array}{cc}{I}_2& \Delta T{I}_2\\ 0_2& I_2\end{array}\right]x_{k-1}+w_k.$$

(66)

The state vector $x_k={[p_{x,k},p_{y,k},v_{x,k},v_{y,k}]}^T$, where $(p_{x,k},p_{y,k})$ represents the position of the target, and $(v_{x,k},v_{y,k})$ denotes the velocity. The sampling time is $\Delta T=0.01$ s. The process noise $w_k$ follows a normal distribution $N(0,Q_k)$, with the covariance matrix $Q_k$ defined as:

$$Q_k=\left[\begin{array}{cc}{\displaystyle \frac{\Delta T^3}{3}}{I}_2& {\displaystyle \frac{\Delta T^2}{2}}{I}_2\\ {\displaystyle \frac{\Delta T^2}{2}}{I}_2& \Delta T{I}_2\end{array}\right].$$

(67)

The measurement data consist of the range and bearing collected by sensors in a cluttered environment, and the measurement equation is expressed as:

$$y_k={\left[\begin{array}{c}\sqrt{p_{x,k}^2+p_{y,k}^2}\\ \mathrm{atan}2(p_{y,k},p_{x,k})\end{array}\right]}^T+{\nu }_k.$$

(68)

Here, the measurement noise ${\nu }_k$ is non-stationary, which can be represented as:

$${\nu }_k\sim \left\{\begin{array}{cc}N(b_1,d_1),\hfill & \mathrm{w}.\mathrm{p}.1-p_1\hfill \\ N(b_1,100d_1),\hfill & \mathrm{w}.\mathrm{p}.p_1\hfill \end{array}\right.,$$

(69)

with

$$b_1=\left[\begin{array}{c}5\\ 0.002\end{array}\right],$$

$$d_1=\left[\begin{array}{cc}100& 0.015\\ 0.015& {10}^{-5}\end{array}\right],$$

$$p_1=0.1.$$

The initial state is assumed to follow a normal distribution $N(x_0,P_0)$, with the initial mean $x_0={[100,1,300,-4]}^T$ and covariance matrix $P_0=\mathrm{diag}\left([25,4,25,4]\right)$, where the units are meters (m) for position and meters per second (m/s) for velocity, and the covariance matrix $P_0$ is a diagonal matrix with the square of the standard deviations for the positions and velocity components.

The results of the filtering solutions are also modeled and compared using the root mean square error (RMSE) and the average root mean square error (ARMSE) metrics. The specific metrics are as follows:

$${\mathrm{RMSE}}_{\mathrm{pos}}=\sqrt{{\displaystyle \frac{1}{M_c}}\sum _{r=1}^{M_c}\left[{(x_k^r-{\widehat{x}}_k^r)}^2+{(y_k^r-{\widehat{y}}_k^r)}^2\right]},$$

(70)

$${\mathrm{ARMSE}}_{\mathrm{pos}}=\sqrt{{\displaystyle \frac{1}{M_{c}N}}\sum _{k=1}^N\sum _{r=1}^{M_c}\left[{(x_k^r-{\widehat{x}}_k^r)}^2+{(y_k^r-{\widehat{y}}_k^r)}^2\right]}.$$

(71)

In these equations, $x_k^r$ and $y_k^r$ are the real coordinates of the time k in the r-th Monte Carlo simulation, while ${\widehat{x}}_k^r$ and ${\widehat{y}}_k^r$ are the corresponding estimated coordinates. $M_c=500$ denotes the total number of Monte Carlo simulations, and $N=40$ is the total sampling time. The RMSE for velocity, denoted as ${\mathrm{RMSE}}_{\mathrm{vel}}$, is similar to the RMSE for position ${\mathrm{RMSE}}_{\mathrm{pos}}$, and thus is not listed here.

We compared the estimation accuracy of the proposed algorithm with that of current robust filtering techniques when the outlier probability $p_1$ is 0.1. The measurement noise exhibited different heavy-tailed noise characteristics in each phase. Figure 12 shows the RMSE of the position state for the proposed and current robust filtering techniques, while Figure 13 shows the RMSE of the velocity state.

Table 6. ARMSEs of filters.

Filters	ARMSEpos	ARMSEvel
RSTKF	13.26	3.90
EKF	13.88	4.78
NSMRKF	19.65	7.82
Proposed Filter	8.31	3.23

lists the average root mean square error (ARMSE) for position and velocity. The simulation results indicate that the proposed filtering algorithm outperforms other algorithms in terms of estimation accuracy for both position and velocity.

5. Conclusions

In this study, a unique Pearson Type VII-Inverse Wishart (PVIW) distribution is successfully introduced with the aim of effectively modeling non-stationary heavy-tailed noise. The PVIW distribution is constructed by incorporating adjustable variables in a hierarchical Gaussian format. On this basis, we propose a novel Kalman filter. The filter employs a variational Bayesian method to jointly infer the state vector, mean vector, and other auxiliary random variables, which realizes the accurate estimation of the system state. Simulation results show that our proposed filter provides higher estimation accuracy in the presence of NSHT measurement noise compared to existing Kalman filters and multi-model filters. In addition, the proposed filtering algorithm exhibits good fast convergence properties. Looking ahead, we believe that the PVIW distribution and the corresponding Kalman filter have great potential for applications in dealing with a wider range of non-stationary noises and complex systems. Future work will explore the application of the filter to multi-sensor fusion and nonlinear systems, and investigate how the algorithm can be further optimized to improve its applicability and robustness to real-world engineering problems. In addition, we will also consider potential improvements to the filter to enhance its performance under extreme noise conditions and explore its applications in other fields.