A New Variational Bayesian-Based Kalman Filter with Unknown Time-Varying Measurement Loss Probability and Non-Stationary Heavy-Tailed Measurement Noise

Abstract

In this paper, a new variational Bayesian-based Kalman filter (KF) is presented to solve the filtering problem for a linear system with unknown time-varying measurement loss probability (UTVMLP) and non-stationary heavy-tailed measurement noise (NSHTMN). Firstly, the NSHTMN was modelled as a Gaussian-Student’s t-mixture distribution via employing a Bernoulli random variable (BM). Secondly, by utilizing another Bernoulli random variable (BL), the form of the likelihood function consisting of two mixture distributions was converted from a weight sum to an exponential product and a new hierarchical Gaussian state-space model was therefore established. Finally, the system state vector, BM, BL, the intermediate random variables, the mixing probability, and the UTVMLP were jointly inferred by employing the variational Bayesian technique. Simulation results revealed that in the scenario of NSHTMN, the proposed filter had a better performance than current algorithms and further improved the estimation accuracy of UTVMLP.

1. Introduction

Under the minimal mean square error criteria, the KF is the optimal estimator for the linear Gaussian state-space model [1,2]. KF has been widely employed in a variety of applications [3,4,5]. Unfortunately, in many practical applications, when the sensor produces intermittent faults, the actual measurement of the sensors may not be accurately represented by the KF measurement model [6,7]. If the random measurement loss occurs, the measurement of the sensors contains only pure noise. In this situation, the estimation accuracy of a typical KF will drop significantly or even diverge. Various filtering methods have been developed to address the measurement loss filtering issue, such as the intermittent KF (IKF) [8,9]. However, IKF has an important assumption: the measurement loss probability is known. In practical applications, the measurement loss probability is usually unknown and the IKF is no longer applicable in this case [7].

In order to address the filtering issues of the unknown measurement loss probability of the linear system, the first Bayesian Kalman filter and the second Bayesian Kalman filter were designed by estimating the process and posterior distribution of the measurement loss, respectively [10]. The above two filters, however, are no longer valid if the unknown measurement loss probability is time-varying. Recently, the variational Bayesian-based adaptive KF (VBAKF) was derived for a linear system with unknown time-varying measurement loss probability (UTVMLP) and both the system vector and UTVMLP are jointly estimated by introducing the variational Bayesian technique [7]. Additionally, VBAKF shows excellent performance in the context of white Gaussian measurement noise with known statistical characteristics. Unfortunately, in realistic engineering applications, measurement outliers may occur at various periods due to environmental changes and unreliable sensors, resulting in NSHTMN, i.e., when the system runs healthily, the measurement noise is the Gaussian-distributed, and when the time-varying measurement outliers erode the system, the measurement noise is heavy-tail-distributed [11,12]. In the scenario of NSHTMN, the estimation accuracy of VBAKF will drop sharply.

Recently, some mixture distribution-based algorithms have been presented to address NSHTMN, such as the Gaussian-Student’s t-mixture distribution-based KF (GSTKF) [13,14]. However, the filtering problem with UTVMLP and NSHTMN cannot directly solved by employing a mixture distribution, that is, under the scenario of UTVMLP and NSHTMN, the current likelihood function is a weighted sum of double-mixture distributions, which is an unclosed and unconjugated distribution that makes the Bayesian inference difficult to employ directly.

In this paper, a new variational Bayesian-based KF is presented to settle the filtering issue for linear discrete-time systems with UTVMLP and NSHTMN. Firstly, the Gaussian-Student’s t-mixture distribution with BM is employed to model the NSHTMN. Secondly, the form of the likelihood function is converted to an exponential product and constructs a new hierarchical Gaussian state-space model by utilizing BL. Thirdly, the variational Bayesian method is introduced to simultaneously estimate the system state vector, BM, BL, the intermediate random variables, the mixing probability, and the UTVMLP. Finally, a numerical simulation experiment reveals that the proposed filter has better estimation accuracy but is more time-consuming than existing filtering algorithms in the scenarios of NSHTMN and UTVMLP.

The contributions of this paper are as follows:

(a)

By employing a Bernoulli-distributed variable, the NSHTMN is modelled as a Gaussian-Student’s t-mixture distribution;

(b)

The measurement likelihood function is converted from the weight sum of two mixture distributions to an exponential product and a new hierarchical Gaussian state-space model is therefore derived;

(c)

The system state vector, UTVMLP, and the unknown variables are simultaneously estimated by utilizing the variational Bayesian technique;

(d)

Numerical simulation results indicate that the proposed filter has better performance than that of existing algorithms in the scenarios of NSHTMN and UTVMLP

2. Problem Formulation

Consider the linear stochastic system with the following state and measurement equations:

$$x_t=F_{t-1}{x}_{t-1}+e_{t-1}$$

(1)

$$y_t=b_t{H}_t{x}_t+g_t$$

(2)

where $x_t\in {ℝ}^m$ denotes the system state vector; $F_{t-1}\in {ℝ}^{m\times m}$ denotes the state transition matrix; $e_t\in {ℝ}^m$ represents the Gaussian-distributed white process noise vector with a zero mean value and covariance matrix $Q_t$; $y_t\in {ℝ}^n$ represents the measurement vector; $H_t\in {ℝ}^{n\times m}$ is the measurement matrix; $g_t\in {ℝ}^n$ is the white NSHTMN vector; and $t$ represents the index of discrete time. The phenomenon of measurement loss is described by introducing the identically distributed and mutually uncorrelated measurement loss defined as the Bernoulli random variable (BL) $b_t$, which is expressed by the following equations.

$$p\left(b_t=1\right)=\mathrm{E}\left[b_t\right]={\gamma }_t$$

(3)

$$p\left(b_t=0\right)=1-\mathrm{E}\left[b_t\right]=1-{\gamma }_t$$

(4)

where ${\gamma }_t\in \left[0,1\right]$ denotes the time-varying measurement loss probability. Note that the value of ${\gamma }_t$ is unknown in this paper. The initial Gaussian-distributed system state vector $x_0$ is the random vector with mean ${\widehat{x}}_{0|0}=0$ and covariance matrix $P_{0|0}$. Additionally, it is assumed that the initial system state vector $x_0$, the noise vectors $e_{t-1}$ and $g_t$, and the Bernoulli random variable $b_t$ are mutually independent.

It can be seen from Equations (1)–(4) that the ideal measurement was received by the sensor when $b_t=0$ and the measurement loss with UTVMLP occurred when $b_t=1$. Meanwhile, the measurement noise is NSHTMN due to measurement outliers, that is, when the system runs healthily, the measurement noise is Gaussian-distributed, and when measurement outliers erode the system, the measurement noise is heavy-tail-distributed. The NSHTMN and UTVMLP can result in estimation errors or even in filtering divergence. Therefore, a new variational Bayesian-based Kalman filter with NSHTMN and UTVMLP will be proposed.

3. Proposed Variational Bayesian-Based Kalman Filter

In this section, a new variational Bayesian-based Kalman filter is proposed to address the filtering issue for a linear system with NSHTMN and UTVMLP. Firstly, the Gaussian-Student’s t-mixture distribution is utilized to model the NSHTMN and the hierarchical form is derived. Secondly, by converting the measurement likelihood function into an exponential multiplication, a new hierarchical Gaussian state-space model is established. Thirdly, by using the variational Bayesian method, the system state and unknown variables are simultaneously estimated. Finally, the required mathematical expectations are given.

3.1. Gaussian-Student’s t-Mixture Distribution

The NSHTMN vector can be modeled as the Gaussian-Student’s t-mixture distribution by employing another mixing-defined Bernoulli random variable (BM), ${\zeta }_t$, and the probability density function (PDF), $p\left(g_t\right)$ is given as

$$p\left(g_t\right)={\displaystyle \sum }_{{\zeta }_t=0}^1{\displaystyle \int }{\left[\mathrm{N}\left(g_t;0,R_t\right)\right]}^{{\zeta }_t}{\left[\mathrm{ST}\left(g_t;0,R_t,\mu \right)\right]}^{\left(1-{\zeta }_t\right)}p\left({\zeta }_t|{\phi }_t\right)p\left({\phi }_t\right)d{\phi }_t {\zeta }_t\in \left\{0,1\right\}$$

(5)

where $\mathrm{N}\left(x;0,\mathsf{\Sigma }\right)$ represents the Gaussian PDF with a zero mean vector and covariance matrix $\mathsf{\Sigma }$, and $\mathrm{ST}\left(x;0,\mathsf{{\rm Y}},\mathsf{\omega }\right)$ represents the student’s t-PDF with a zero mean vector, covariance matrix $\mathsf{{\rm Y}}$, and degree of freedom (dof) parameter $\mathsf{\omega }$. $R_t$ represents the covariance matrix of the nominal measurement noise. The PDF of the mixing probability ${\phi }_t$ and the probability mass function (PMF) of ${\zeta }_t$ are defined as follows, respectively.

$$p\left({\phi }_t\right)=\mathrm{Be}\left({\phi }_t;h_0,1-h_0\right)$$

(6)

$$p\left({\zeta }_t|{\phi }_t\right)={\phi }_t{}^{{\zeta }_t}{\left(1-{\phi }_t\right)}^{\left(1-{\zeta }_t\right)}$$

(7)

where $\mathrm{Be}\left(x;\mathsf{\sigma },\kappa \right)$ represents the Beta PDF with shape parameters $\mathsf{\sigma }$ and $\kappa$.

Due to the hierarchical properties of the student’s t-distribution, Equation (5) can be rewritten as such:

$$p\left(g_t\right)={\displaystyle \sum }_{{\zeta }_t=0}^1{{\displaystyle \int }}_0^{+\infty }{\displaystyle \int }{\left[\mathrm{N}\left(g_t;0,R_t\right)\right]}^{{\zeta }_t}{\left[\mathrm{N}\left(g_t;0,R_t/{\beta }_t\right)\right]}^{\left(1-{\zeta }_t\right)}p\left({\beta }_t\right)p\left({\zeta }_t|{\phi }_t\right)p\left({\phi }_t\right)d{\phi }_{t}d{\beta }_t$$

(8)

$$p\left({\beta }_t\right)=\mathrm{G}\left({\beta }_t,\frac{\mu }{2},\frac{\mu }{2}\right)$$

(9)

where $\mathrm{G}\left(x,a,b\right)$ represents the Gamma PDF with shape parameter $a$ and rate parameter $b$, and ${\beta }_t$ represents the intermediate random variable.

3.2. New Hierarchical Gaussian State-Space Model (HGSSM)

According to Equations (2)–(4), the measurement likelihood PDF is derived as Based on Equation (2), the following equation can be obtained.

$$\begin{array}{ll}p\left(y_t|x_t,{\gamma }_t\right)& ={\displaystyle {\displaystyle \sum }_{b_t=0}^1}p\left(y_t,b_t|x_t,{\gamma }_t\right)\\ & =p\left(y_t|x_t,b_t=1\right)p\left(b_t=1\right)+p\left(y_t|x_t,b_t=0\right)p\left(b_t=0\right)\\ & =\left(1-{\gamma }_t\right)p\left(y_t|x_t,b_t=1\right)+{\gamma }_{t}p\left(y_t|x_t,b_t=0\right)\end{array}$$

(10)

$$p\left(y_t|x_t,b_t=1\right)=p_{g_t}\left(y_t-H_t{x}_t\right)$$

(11)

$$p\left(y_t|x_t,b_t=0\right)=p_{g_t}\left(y_t\right)$$

(12)

where $p_{g_t}(·)$ represents the measurement noise PD. Substituting Equations (11) and (12) in Equation (10) results in

$$p\left(y_t|x_t,{\gamma }_t\right)=\left(1-{\gamma }_t\right)p_{g_t}\left(y_t-H_t{x}_t\right)+{\gamma }_t{p}_{g_t}\left(y_t\right)$$

(13)

Remark 1.

The measurement likelihood PDF in Equation (13) is an unclosed and unconjugated weighted sum form, and it is impossible to infer the system state vector and unknown parameters directly by utilizing the variational Bayesian. The weighted sum will then be converted into an exponential multiplication form to address this problem.

The PMF of BL $b_t$ is given as

$$p\left(b_t|{\gamma }_t\right)={\left(1-{\gamma }_t\right)}^{b_t}{\gamma }_t{}^{\left(1-b_t\right)}$$

(14)

Exploiting Equations (13) and (14), the measurement likelihood PDF is reformulated as

$$\begin{array}{ll}p\left(y_t|x_t,{\gamma }_t\right)& ={\displaystyle {\displaystyle \sum }_{b_t=0}^1}p\left(y_t|x_t,b_t\right)p\left(b_t|{\gamma }_t\right)\\ & ={\displaystyle {\displaystyle \sum }_{b_t=0}^1}\left[{\left(1-{\gamma }_t\right)}^{b_t}{\left[p_{g_t}\left(y_t-H_t{x}_t\right)\right]}^{b_t}{\gamma }_t{}^{\left(1-b_t\right)}{\left[p_{g_t}\left(y_t\right)\right]}^{\left(1-b_t\right)}\right]\\ & ={\displaystyle {\displaystyle \sum }_{b_t=0}^1}\left[{\left[p_{g_t}\left(y_t-H_t{x}_t\right)\right]}^{b_t}{\left[p_{g_t}\left(y_t\right)\right]}^{\left(1-b_t\right)}\right]p\left(b_t|{\gamma }_t\right)\end{array}$$

(15)

According Equation (15), the exponential multiplication-formed likelihood PDF $p\left(y_t|x_t,b_t\right)$ is given as follows.

$$p\left(y_t|x_t,b_t\right)={\left[p_{g_t}\left(y_t-H_t{x}_t\right)\right]}^{b_t}{\left[p_{g_t}\left(y_t\right)\right]}^{\left(1-b_t\right)}$$

(16)

Remark 2.

The variational Bayesian method must select the suitable conjugate-prior distributions for unknown variables. Therefore, the appropriate prior PDFs to construct a new HGSSM are selected.

The one-step predicted PDF $p\left(x_t|y_{1:t-1}\right)$ of system state vector $x_t$ is assumed as being Gaussian distributed as follows.

$$p\left(x_t|y_{1:t-1}\right)=\mathrm{N}\left(x_t;{\widehat{x}}_{t|t-1},P_{t|t-1}\right)$$

(17)

where ${\widehat{x}}_{t|t-1}$ represents the mean vector and $P_{t|t-1}$ represents the covariance matrix. Both ${\widehat{x}}_{t|t-1}$ and $P_{t|t-1}$ can be updated by the typical Kalman filter, which is given as

$${\widehat{x}}_{t|t-1}=F_{t-1}{\widehat{x}}_{t-1|t-1}$$

(18)

$$P_{t|t-1}=H_{t-1}{P}_{t-1|t-1}{H}_{t-1}^T+Q_{t-1}$$

(19)

In employing Equations (8), (9) and (16), the conditional likelihood PDF $p\left(y_t|x_t,b_t,\beta ,{\zeta }_t\right)$ is derived as

$$\begin{array}{ll}p\left(y_t|x_t,b_t,{\beta }_t,{\zeta }_t\right)& ={\left[\mathrm{N}{\left(y_t;H_t{x}_t,R_t\right)}^{{\zeta }_t}\mathrm{N}{\left(y_t;H_t{x}_t,R_t/{\beta }_t\right)}^{\left(1-{\zeta }_t\right)}\right]}^{b_t}\\ & \times {\left[\mathrm{N}{\left(y_t;0,R_t\right)}^{{\zeta }_t}\mathrm{N}{\left(y_t;0,R_t/{\beta }_t\right)}^{\left(1-{\zeta }_t\right)}\right]}^{\left(1-b_t\right)}\\ & =\mathrm{N}{\left(y_t;H_t{x}_t,R_t\right)}^{{\zeta }_t{b}_t}\mathrm{N}{\left(y_t;H_t{x}_t,R_t/{\beta }_t\right)}^{\left(1-{\zeta }_t\right)b_t}\\ & \times N{\left(y_t;0,R_t\right)}^{{\zeta }_t\left(1-b_t\right)}N{\left(y_t;0,R_t/{\beta }_t\right)}^{\left(1-{\zeta }_t\right)\left(1-b_t\right)}\end{array}$$

(20)

It can be seen from Equations (6)–(9), (13) and (20) that the measurement vector $y_t$ depends on system state vector $x_t$, intermediate random variable ${\beta }_t$, BM ${\zeta }_t$, BL $b_t$, mixing probability ${\phi }_t$, and measurement loss probability ${\gamma }_t$. The following joint-prior PDF must be calculated, i.e.,

$$\begin{array}{c}p\left(\mathsf{\Xi }|y_{1:t-1}\right)=p\left(x_t|y_{1:t-1}\right)p\left({\beta }_t\right)p\left(b_t|{\gamma }_t\right)p\left({\zeta }_t|{\phi }_t\right)p\left({\phi }_t\right)p\left({\gamma }_t|y_{1:t-1}\right),\\ \mathsf{\Xi }\triangleq \left\{x_t,{\beta }_t,b_t,{\zeta }_t,{\gamma }_t,\phi \right\}\end{array}$$

(21)

where the definitions of $p\left({\gamma }_t\right)$, $p\left({\zeta }_t|\phi \right)$, $p\left(\beta \right)$, $p\left(b_t|{\gamma }_t\right)$, and $p\left(x_t|y_{1:t-1}\right)$ are given in Equations (6), (7), (9), (14) and (17), respectively. Additionally, $p\left({\gamma }_t|y_{1:t-1}\right)$ denotes the prior PDF of the time-varying measurement loss probability, which can be assumed as the following Beta distribution.

$$\begin{array}{ll}p\left({\gamma }_t|y_{1:t-1}\right)& =\mathrm{Be}\left({\gamma }_t;{\widehat{\eta }}_{t|t-1},{\widehat{\delta }}_{t|t-1}\right)\\ & ={\displaystyle \int }p\left({\gamma }_{t-1}|y_{1:t-1}\right)p\left({\gamma }_t|{\gamma }_{t-1}\right)d{\gamma }_{t-1}\left(\mathrm{Bayes}’ \mathrm{theorem}\right)\end{array}$$

(22)

where the shape parameters ${\widehat{\eta }}_{t|t-1}$ and ${\widehat{\delta }}_{t|t-1}$ can be calculated by introducing the forgetting factor $\rho \in (0,1]$ as follows.

$${\widehat{\eta }}_{t|t-1}=\rho {\widehat{\eta }}_{t-1|t-1}$$

(23)

$${\widehat{\delta }}_{t|t-1}=\rho {\widehat{\delta }}_{t-1|t-1}$$

(24)

where ${\widehat{\eta }}_{t-1|t-1}$ and ${\widehat{\delta }}_{t-1|t-1}$ represent posterior shape parameters.

The proposed new HGSSM is comprised of Equations (14) and (17)–(24). System state vector $x_t$, intermediate random variable ${\beta }_t$, BM ${\zeta }_t$, BL $b_t$, mixing probability ${\phi }_t$, and measurement loss probability ${\gamma }_t$ will be simultaneously estimated by utilizing the variational Bayesian method.

3.3. Variational Bayesian Approximation of the Joint Posterior PDFs

Aiming at the estimation of the unknown variables of the new HGSSM, the joint posterior PDF $p\left(\mathsf{\Xi }|y_t\right)$ with $\mathsf{\Xi }\triangleq \left\{x_t,{\beta }_t,b_t,{\zeta }_t,{\phi }_t,{\gamma }_t\right\}$ is required to be solved. However, the analytical solution of $p\left(\mathsf{\Xi }|y_t\right)$ is not accessible. The variational Bayesian approach is therefore employed to determine an approximate PDF for $p\left(\mathsf{\Xi }|y_t\right)$ and to compute an approximate solution [15,16,17], i.e.,

$$p\left(\mathsf{\Xi }|y_{1:t}\right)\approx q_a\left(x_t\right)q_a\left({\beta }_t\right)q_a\left(b_t\right)q_a\left({\zeta }_t\right)q_a\left({\phi }_t\right)q_a\left({\gamma }_t\right)$$

(25)

where $\theta$ represents an arbitrary element of $\mathsf{\Xi }$ and $q_a\left(\theta \right)$ denotes the approximate PDF or PMF. By minimizing the Kullback–Leibler divergence (KLD) between $p\left(\mathsf{\Xi }|y_{1:t}\right)$ and $q_a\left(x_t\right)q_a\left({\beta }_t\right)q_a\left(b_t\right)q_a\left({\zeta }_t\right)q_a\left({\phi }_t\right)q_a\left({\gamma }_t\right)$, $q_a\left(\theta \right)$ can be obtained as follows.

$$\begin{array}{c}\left\{q_a\left(x_t\right),q_a\left({\beta }_t\right),q_a\left(b_t\right),q_a\left({\zeta }_t\right),q_a\left({\phi }_t\right),q_a\left({\gamma }_t\right)\right\}=\mathrm{argminKLD}\\ \left(q_a\left(x_t\right)q_a\left({\beta }_t\right)q_a\left(b_t\right)q_a\left({\zeta }_t\right)q_a\left({\phi }_t\right)q_a\left({\gamma }_t\right)\left|\right|p\left(\mathsf{\Xi }|y_{1:t}\right)\right)\end{array}$$

(26)

$$\mathrm{KLD}\left(q_a\left(A\right)\left|\right|p\left(A\right)\right)\triangleq {\displaystyle \int }{q}_a\left(A\right)\log \frac{q_a\left(A\right)}{p\left(A\right)}dA$$

(27)

where $\mathrm{KLD}\left(q_a\left(A\right)\left|\right|p\left(A\right)\right)$ represents the KLD between $q_a\left(A\right)$ and $p\left(A\right)$, and the optimal solution of Equation (26) can be calculated via the following formula [15,17].

$$\log q_a\left(\theta \right)={\mathrm{E}}_{{\mathsf{\Xi }}^{-\theta }}\left[\log p\left(\mathsf{\Xi },y_{1:t}\right)\right]+c_{\theta }$$

(28)

where $\mathrm{E}\left(\theta \right)$ denotes the mathematical expectation operation, ${\mathsf{\Xi }}^{-\theta }$ signifies a grouping of all the components in $\mathsf{\Xi }$ apart from $\theta$, and the constant with regard to $\theta$ is denoted by $c_{\theta }$. Additionally, the fixed-point iteration technique is utilized to derive the approximate formation of $q_a\left(\theta \right)$ due to the fact that estimated parameters are mutually coupled.

The joint PDF $p\left(\mathsf{\Xi },y_{1:t}\right)$ in Equation (26) can be derived as

$$\begin{array}{ll}p\left(\mathsf{\Xi },y_{1:t}\right)& =\mathrm{N}{\left(y_t;H_t{x}_t,R_t\right)}^{{\zeta }_t{b}_t}\mathrm{N}{\left(y_t;H_t{x}_t,R_t/{\beta }_t\right)}^{\left(1-{\zeta }_t\right)b_t}\\ & \times \mathrm{N}{\left(y_t;0,R_t\right)}^{{\zeta }_t\left(1-b_t\right)}\mathrm{N}{\left(y_t;0,R_t/{\beta }_t\right)}^{\left(1-{\zeta }_t\right)\left(1-b_t\right)}\mathrm{N}\left(x_t;{\widehat{x}}_{t|t-1},P_{t|t-1}\right)\\ & \times \mathrm{G}\left({\beta }_t,\frac{{\mu }_t}{2},\frac{{\mu }_t}{2}\right){\phi }_t{}^{{\zeta }_t}{\left(1-{\phi }_t\right)}^{\left(1-{\zeta }_t\right)}{\left(1-{\gamma }_t\right)}^{b_t}{\gamma }_t{}^{\left(1-b_t\right)}\\ & \times \mathrm{Be}\left({\phi }_t;h_0,1-h_0\right)\mathrm{Be}\left({\gamma }_t;{\widehat{\eta }}_{t|t-1},{\widehat{\delta }}_{t|t-1}\right)p\left(y_{1:t-1}\right)\end{array}$$

(29)

Proposition 1.

Let $\theta =x_t$ and by using Equation (29) in (28), $q_a^{\left(s+1\right)}\left(x_t\right)$ can be updated as Gaussian, i.e.,

$$q_a^{\left(s+1\right)}\left(x_t\right)=\mathrm{N}\left(x_t;{\widehat{x}}_{t|t}^{\left(s+1\right)},P_{t|t}^{\left(s+1\right)}\right)$$

(30)

where $q_a^{\left(s+1\right)}(·)$ represents the approximate PDF in the $\left(s+1\right)th$ iteration, while the mean vector ${\widehat{x}}_{t|t}^{\left(s+1\right)}$ and the covariance matrix $P_{t|t}^{\left(s+1\right)}$ are assumed to be updated in accordance with the traditional Kalman filter as follows.

$${\widehat{x}}_{t|t}^{\left(s+1\right)}={\widehat{x}}_{t|t-1}+K_t^{\left(s+1\right)}\left(y_t-H_t{\widehat{x}}_{t|t-1}\right)$$

(31)

$$P_{t|t}^{\left(s+1\right)}={\widehat{P}}_{t|t-1}^{\left(s+1\right)}-K_t^{\left(s+1\right)}{H}_t{\widehat{P}}_{t|t-1}^{\left(s+1\right)}$$

(32)

$$K_t^{\left(s+1\right)}={\widehat{P}}_{t|t-1}^{\left(s+1\right)}{H}_t^T{\left(H_t{\widehat{P}}_{t|t-1}^{\left(s+1\right)}{H}_t^T+{\tilde{R}}_t^{\left(s+1\right)}\right)}^{-1}$$

(33)

where $K_t^{\left(s+1\right)}$ represents the Kalman gain matrix. The modified measurement noise covariance matrix at $\left(s+1\right)th$ iteration ${\tilde{R}}_t^{\left(s+1\right)}$ is formulated as

$${\tilde{R}}_t^{\left(s+1\right)}=\frac{R_t}{{\mathrm{E}}^{\left(s\right)}\left[{\zeta }_t\right]{\mathrm{E}}^{\left(s\right)}\left[b_t\right]+{\mathrm{E}}^{\left(s\right)}\left[1-{\zeta }_t\right]{\mathrm{E}}^{\left(s\right)}\left[b_t\right]{\mathrm{E}}^{\left(s\right)}\left[{\beta }_t\right]}$$

(34)

where ${\mathrm{E}}^{\left(s\right)}[·]$ represents the mathematical expectation of variables in the $sth$ iteration.

Proof: see Appendix A.

Proposition 2.

Let $\theta ={\beta }_t$ and by using Equation (29) in Equation (28), $q_a^{\left(s+1\right)}\left({\beta }_t\right)$ can be updated as Gamma, i.e.,

$$q_a^{\left(s+1\right)}\left({\beta }_t\right)=\mathrm{G}\left({\beta }_t;{\pi }_t^{\left(s+1\right)},{\nu }_t^{\left(s+1\right)}\right)$$

(35)

where the shape parameter ${\pi }_t^{\left(s+1\right)}$ and rate parameter ${\nu }_t^{\left(s+1\right)}$ are formulated as

$${\pi }_t^{\left(s+1\right)}=0.5\left(n{\mathrm{E}}^{\left(s+1\right)}\left[1-{\zeta }_t\right]+{\mu }_t\right)$$

(36)

$${\nu }_t^{\left(s+1\right)}=0.5\left[\mathrm{tr}\left(G_t^{\left(s+1\right)}{\mathbf{R}}_t^{-1}\right)+{\mu }_t\right]$$

(37)

where $n$ represents the dimension of the measurement vector, $\mathrm{tr}(·)$ represents the trace operation, and $G_t^{\left(s+1\right)}$ is defined as

$$\begin{array}{ll}{G}_t^{\left(s+1\right)}& ={\mathrm{E}}^{\left(s\right)}\left[1-{\zeta }_t\right]\{{\mathrm{E}}^{\left(s\right)}\left[b_t\right]\left[H_t{P}_{t|t}^{\left(s+1\right)}{H}_t^T+\left(y_t-H_t{\widehat{x}}_{t|t}^{\left(s+1\right)}\right){\left(y_t-H_t{\widehat{x}}_{t|t}^{\left(s+1\right)}\right)}^T\right]\\ & +{\mathrm{E}}^{\left(s\right)}\left[1-{\zeta }_t\right]{\mathrm{E}}^{\left(s\right)}\left[1-b_t\right]y_t{\left(y_t\right)}^T\}\end{array}$$

(38)

Proof: see Appendix B.

Proposition 3.

Let $\theta =b_t$ and by using Equation (29) in Equation (28), $q_a^{\left(s+1\right)}\left(b_t\right)$ can be updated as the Bernoulli distribution. The posterior probabilities $p\left(b_t=1\right)$ and $p\left(b_t=0\right)$ of BL $b_t$ are given as

$$p^{\left(s+1\right)}\left(b_t=1\right){∆}^{\left(s+1\right)}\exp \left\{{\mathrm{E}}^{\left(s\right)}\left[\log \left(1-{\gamma }_t\right)\right]-0.5\mathrm{tr}\left(C_t^{\left(s+1\right)}{R}_t^{-1}\right)+0.5n{\mathrm{E}}^{\left(s\right)}\left[1-{\zeta }_t\right]{\mathrm{E}}^{\left(s\right)}\left[\log \left({\beta }_t\right)\right]\right\}$$

(39)

$$p^{\left(s+1\right)}\left(b_t=0\right)={∆}^{\left(s+1\right)}\exp \left\{{\mathrm{E}}^{\left(s\right)}\left[\log \left({\gamma }_t\right)\right]-0.5\mathrm{tr}\left(D_t^{\left(s+1\right)}{R}_t^{-1}\right)+0.5n{\mathrm{E}}^{\left(s\right)}\left[1-{\zeta }_t\right]{\mathrm{E}}^{\left(s\right)}\left[\log \left({\beta }_t\right)\right]\right\}$$

(40)

where ${∆}^{\left(s+1\right)}$ represents the normalizing constant and the parameters $C_t^{\left(s+1\right)}$ and $D_t^{\left(s+1\right)}$ are, respectively, defined as

$$C_t^{\left(s+1\right)}=\left\{{\mathrm{E}}^{\left(s\right)}\left[{\zeta }_t\right]+{\mathrm{E}}^{\left(s\right)}\left[{\beta }_t\right]{\mathrm{E}}^{\left(s\right)}\left[1-{\zeta }_t\right]\right\}\left[H_t{P}_{t|t}^{\left(s+1\right)}{H}_t^T+\left(y_t-H_t{\widehat{x}}_{t|t}^{\left(s+1\right)}\right){\left(y_t-H_t{\widehat{x}}_{t|t}^{\left(s+1\right)}\right)}^T\right]$$

(41)

$$D_t^{\left(s+1\right)}=\left\{{\mathrm{E}}^{\left(s\right)}\left[{\zeta }_t\right]+{\mathrm{E}}^{\left(s\right)}\left[{\beta }_t\right]{\mathrm{E}}^{\left(s\right)}\left[1-{\zeta }_t\right]\right\}{y}_t{y}_t^T$$

(42)

Proof: see Appendix C.

Proposition 4.

Let $\theta ={\zeta }_t$ and by using Equation (29) in (28), $q_a^{\left(s+1\right)}\left({\zeta }_t\right)$ can be also updated as the Bernoulli distribution. The posterior probabilities $p\left({\zeta }_t=1\right)$ and $p\left({\zeta }_t=0\right)$ of BM ${\zeta }_t$ are given as

$$p^{\left(s+1\right)}\left({\zeta }_t=1\right)={\nabla }^{\left(s+1\right)}\exp \left\{{\mathrm{E}}^{\left(s\right)}\left[\log \left({\phi }_t\right)\right]-0.5\mathrm{tr}\left(V_t^{\left(s+1\right)}{R}_t^{-1}\right)\right\}$$

(43)

$$p^{\left(s+1\right)}\left({\zeta }_t=0\right)={\nabla }^{\left(s+1\right)}\exp \left\{{\mathrm{E}}^{\left(s\right)}\left[\log \left(1-{\phi }_t\right)\right]-0.5\mathrm{tr}\left(W_t^{\left(s+1\right)}{R}_t^{-1}\right)+0.5n{\mathrm{E}}^{\left(s\right)}\left[\log \left({\beta }_t\right)\right]\right\}$$

(44)

where ${\nabla }^{\left(s+1\right)}$ also represents the normalizing constant and the definitions of parameters $V_t^{\left(s+1\right)}$ and $W_t^{\left(s+1\right)}$ are, respectively, given as

$$\begin{array}{ll}{V}_t^{\left(s+1\right)}& ={\mathrm{E}}^{\left(s\right)}\left[b_t\right]\left[H_t{P}_{t|t}^{\left(s+1\right)}{H}_t^T+\left(y_t-H_t{\widehat{x}}_{t|t}^{\left(s+1\right)}\right){\left(y_t-H_t{\widehat{x}}_{t|t}^{\left(s+1\right)}\right)}^T\right]\\ & +{\mathrm{E}}^{\left(s\right)}\left[1-b_t\right]y_t{y}_t^T\end{array}$$

(45)

$$\begin{array}{ll}{W}_t^{\left(s+1\right)}& ={\mathrm{E}}^{\left(s\right)}\left[b_t\right]{\mathrm{E}}^{\left(s\right)}\left[{\beta }_t\right]\left[H_t{P}_{t|t}^{\left(s+1\right)}{H}_t^T+\left(y_t-H_t{\widehat{x}}_{t|t}^{\left(s+1\right)}\right){\left(y_t-H_t{\widehat{x}}_{t|t}^{\left(s+1\right)}\right)}^T\right]\\ & +{\mathrm{E}}^{\left(s\right)}\left[1-b_t\right]{\mathrm{E}}^{\left(s\right)}\left[{\beta }_t\right]y_t{y}_t^T\end{array}$$

(46)

Proof: see Appendix D.

Proposition 5.

Let $\theta ={\phi }_t$ and by using Equation (29) in (28), $q_a^{\left(s+1\right)}\left({\phi }_t\right)$ can be updated as the Beta distribution, i.e.,

$$q_a^{\left(s+1\right)}\left({\phi }_t\right)=\mathrm{Be}\left({\phi }_t;h_t^{\left(s+1\right)},d_t^{\left(s+1\right)}\right)$$

(47)

where the shape parameters $h_t^{\left(s+1\right)}$ and $d_t^{\left(s+1\right)}$ are formulated as

$$h_t^{\left(s+1\right)}={\mathrm{E}}^{\left(s+1\right)}\left[{\zeta }_t\right]+h_0$$

(48)

$$d_t^{\left(s+1\right)}={\mathrm{E}}^{\left(s+1\right)}\left[1-{\zeta }_t\right]+1-h_0$$

(49)

Proof: see Appendix E.

Proposition 6.

Let $\theta ={\gamma }_t$ and by using Equation (29) in Equation (28), $q_a^{\left(s+1\right)}\left({\gamma }_t\right)$ can be also updated as the Beta distribution, i.e.,

$$q_a^{\left(s+1\right)}\left({\gamma }_t\right)=\mathrm{Be}\left({\gamma }_t;{\widehat{\eta }}_t^{\left(s+1\right)},{\widehat{\delta }}_t^{\left(s+1\right)}\right)$$

(50)

where the definitions of shape parameters ${\widehat{\eta }}_t^{\left(s+1\right)}$ and ${\widehat{\delta }}_t^{\left(s+1\right)}$ are given as

$${\widehat{\eta }}_t^{\left(s+1\right)}={\mathrm{E}}^{\left(s+1\right)}\left[1-b_t\right]+{\widehat{\eta }}_{t|t-1}$$

(51)

$${\widehat{\delta }}_t^{\left(s+1\right)}={\mathrm{E}}^{\left(s+1\right)}\left[b_t\right]+{\widehat{\delta }}_{t|t-1}$$

(52)

Proof: see Appendix F.

3.4. Calculation of the Required Mathematical Expectations

The required mathematical expectations ${\mathrm{E}}^{\left(s+1\right)}\left[{\beta }_t\right]$, ${\mathrm{E}}^{\left(s+1\right)}\left[\log \left({\beta }_t\right)\right]$,${\mathrm{E}}^{\left(s+1\right)}\left[b_t\right]$, ${\mathrm{E}}^{\left(s+1\right)}\left[1-b_t\right]$, ${\mathrm{E}}^{\left(s+1\right)}\left[{\zeta }_t\right]$, ${\mathrm{E}}^{\left(s+1\right)}\left[1-{\zeta }_t\right]$, ${\mathrm{E}}^{\left(s+1\right)}\left[\log \left({\phi }_t\right)\right]$, ${\mathrm{E}}^{\left(s+1\right)}\left[\log \left(1-{\phi }_t\right)\right]$, ${\mathrm{E}}^{\left(s+1\right)}\left[\log \left({\gamma }_t\right)\right],$ and ${\mathrm{E}}^{\left(s+1\right)}\left[\log \left(1-{\gamma }_t\right)\right]$ in Section 3.3 are defined, respectively, as follows:

$${\mathrm{E}}^{\left(s+1\right)}\left[{\beta }_t\right]=\frac{{\pi }_t^{\left(s+1\right)}}{{\nu }_t^{\left(s+1\right)}}$$

(53)

$${\mathrm{E}}^{\left(s+1\right)}\left[\log \left({\beta }_t\right)\right]=-\Psi \left({\pi }_t^{\left(s+1\right)}\right)-\log \left({\nu }_t^{\left(s+1\right)}\right)$$

(54)

$${\mathrm{E}}^{\left(s+1\right)}\left[b_t\right]=\frac{p^{\left(s+1\right)}\left(b_t=1\right)}{p^{\left(s+1\right)}\left(b_t=1\right)p^{\left(s+1\right)}\left(b_t=0\right)}$$

(55)

$${\mathrm{E}}^{\left(s+1\right)}\left[1-b_t\right]=1-{\mathrm{E}}^{\left(s+1\right)}\left[b_t\right]$$

(56)

$${\mathrm{E}}^{\left(s+1\right)}\left[{\zeta }_t\right]=\frac{p^{\left(s+1\right)}\left({\zeta }_t=1\right)}{p^{\left(s+1\right)}\left({\zeta }_t=1\right)p^{\left(s+1\right)}\left({\zeta }_t=0\right)}$$

(57)

$${\mathrm{E}}^{\left(s+1\right)}\left[1-{\zeta }_t\right]=1-{\mathrm{E}}^{\left(s+1\right)}\left[{\zeta }_t\right]$$

(58)

$${\mathrm{E}}^{\left(s+1\right)}\left[\log \left({\phi }_t\right)\right]=\Psi \left(h_t^{\left(s+1\right)}\right)-\Psi \left(h_t^{\left(s+1\right)}+d_t^{\left(s+1\right)}\right)$$

(59)

$${\mathrm{E}}^{\left(s+1\right)}\left[\log \left(1-{\phi }_t\right)\right]=\Psi \left(d_t^{\left(s+1\right)}\right)-\Psi \left(h_t^{\left(s+1\right)}+d_t^{\left(s+1\right)}\right)$$

(60)

$${\mathrm{E}}^{\left(s+1\right)}\left[\log \left({\gamma }_t\right)\right]=\Psi \left({\eta }_t^{\left(s+1\right)}\right)-\Psi \left({\eta }_t^{\left(s+1\right)}+{\delta }_t^{\left(s+1\right)}\right)$$

(61)

$${\mathrm{E}}^{\left(s+1\right)}\left[\log \left(1-{\gamma }_t\right)\right]=\Psi \left({\delta }_t^{\left(s+1\right)}\right)-\Psi \left({\eta }_t^{\left(s+1\right)}+{\delta }_t^{\left(s+1\right)}\right)$$

(62)

where $\Psi (·)$ represents the digamma function [18].

The presented variational Bayesian-based Kalman filter with UTVMLP and NSHTMN consists of Equations (18), (19) and (30)–(62).

Table 1. The proposed variational Bayesian-based Kalman filter with UTVMLP and NSHTMN (one-time step).

Inputs: ${\widehat{x}}_{t-1|t-1}$, $P_{t-1|t-1}$, $Q_{t-1|t-1}$, $R_{t-1|t-1}$, $y_t$, $F_{t-1}$, $H_t$, $n$, $m$, ${\mu }_t$, $h_0$, ${\widehat{\eta }}_{t-1}$, ${\widehat{\delta }}_{t-1}$, $N_I$, $\varsigma$

Time update:

1. Obtain ${\widehat{x}}_{t|t-1}$ and $P_{t|t-1}$ utilizing Equations (18) and (19) (time update of typical Kalman filter).

Variational measurement update:

2. Initialization: ${\widehat{x}}_{t|t}^{\left(0\right)}=x_{t|t-1}$, $P_{t|t}^{\left(0\right)}=P_{t|t-1}$, ${\mathrm{E}}^{\left(0\right)}\left[{\beta }_t\right]=1$, ${\mathrm{E}}^{\left(0\right)}\left[\log \left({\beta }_t\right)\right]=0$, ${\mathrm{E}}^{\left(0\right)}\left[b_t\right]={\widehat{\eta }}_{t-1}/{\widehat{\delta }}_{t-1}$, ${\mathrm{E}}^{\left(0\right)}\left[1-b_t\right]=1-{\mathrm{E}}^{\left(0\right)}\left[b_t\right]$, ${\widehat{\eta }}_t^{\left(0\right)}={\widehat{\eta }}_{t-1}$, ${\widehat{\delta }}_t^{\left(0\right)}={\widehat{\delta }}_{t-1}$, ${\mathrm{E}}^{\left(0\right)}\left[{\zeta }_t\right]=1$, ${\mathrm{E}}^{\left(0\right)}\left[{\zeta }_t\right]=1-{\mathrm{E}}^{\left(0\right)}\left[{\zeta }_t\right]$, ${\mathrm{E}}^{\left(0\right)}\left[\log \left({\phi }_t\right)\right]=\Psi \left(h_0\right)-\Psi \left(1\right)$, ${\mathrm{E}}^{\left(s+1\right)}\left[\log \left(1-{\phi }_t\right)\right]=\Psi \left(1-h_0\right)-\Psi \left(1\right)$, ${\mathrm{E}}^{\left(s+1\right)}\left[\log \left({\gamma }_t\right)\right]=\Psi \left({\widehat{\eta }}_t^{\left(0\right)}\right)-\Psi \left({\widehat{\eta }}_t^{\left(0\right)}+{\widehat{\delta }}_t^{\left(0\right)}\right)$, ${\mathrm{E}}^{\left(s+1\right)}\left[\log \left(1-{\gamma }_t\right)\right]=\Psi \left({\widehat{\delta }}_t^{\left(0\right)}\right)-\Psi \left({\widehat{\eta }}_t^{\left(0\right)}+{\widehat{\delta }}_t^{\left(0\right)}\right)$

for$s=0:N_I-1$.

3. Update $q_a^{\left(s+1\right)}\left(x_t\right)$ by Equation (30).

4. Obtain ${\widehat{x}}_{t|t}^{\left(s+1\right)}$, $P_{t|t}^{\left(s+1\right)}$, and ${\tilde{R}}_t^{\left(s+1\right)}$ by utilizing Equations (31)–(34) (typical Kalman filter).

5. Update the Gamma-distributed $q_a^{\left(s+1\right)}\left({\beta }_t\right)$ by Equation (35).

6. Obtain ${\pi }_t^{\left(s+1\right)}$, ${\nu }_t^{\left(s+1\right)},$ and $G_t^{\left(s+1\right)}$ by utilizing Equations (36)–(38).

7. Update the Bernoulli-distributed $q_a^{\left(s+1\right)}\left(b_t\right)$.

8. Obtain $p^{\left(s+1\right)}\left(b_t=1\right)$, $p^{\left(s+1\right)}\left(b_t=0\right)$, $C_t^{\left(s+1\right)}$, and $D_t^{\left(s+1\right)}$ by utilizing Equations (39)–(42).

9. Obtain ${\mathrm{E}}^{\left(s+1\right)}\left[b_t\right]$ and ${\mathrm{E}}^{\left(s+1\right)}\left[1-b_t\right]$ by utilizing Equations (55) and (56).

10. Update the Bernoulli-distributed $q_a^{\left(s+1\right)}\left({\zeta }_t\right)$.

11. Obtain $p^{\left(s+1\right)}\left({\zeta }_t=1\right)$, $p^{\left(s+1\right)}\left({\zeta }_t=0\right)$, $V_t^{\left(s+1\right)},$ and $W_t^{\left(s+1\right)}$ by utilizing Equations (43)–(46).

12. Obtain ${\mathrm{E}}^{\left(s+1\right)}\left[{\zeta }_t\right]$ and ${\mathrm{E}}^{\left(s+1\right)}\left[1-{\zeta }_t\right]$ by utilizing Equations (57) and (58).

13. Update the Beta-distributed $q_a^{\left(s+1\right)}\left({\phi }_t\right)$ by Equation (47).

14. Obtain $h_t^{\left(s+1\right)}$ and $d_t^{\left(s+1\right)}$ by utilizing Equations (48) and (49).

15. Obtain ${\mathrm{E}}^{\left(s+1\right)}\left[\log \left({\phi }_t\right)\right]$ and ${\mathrm{E}}^{\left(s+1\right)}\left[\log \left(1-{\phi }_t\right)\right]$ by utilizing Equations (59) and (60).

16. Update the Beta-distributed $q_a^{\left(s+1\right)}\left({\gamma }_t\right)$ by Equation (50).

17. Obtain ${\widehat{\eta }}_t^{\left(s+1\right)}$ and ${\widehat{\delta }}_t^{\left(s+1\right)}$ by utilizing Equations (51) and (52).

18. Obtain ${\mathrm{E}}^{\left(s+1\right)}\left[\log \left({\gamma }_t\right)\right]$ and ${\mathrm{E}}^{\left(s+1\right)}\left[\log \left(1-{\gamma }_t\right)\right]$ by utilizing Equations (61) and (62).

19. If $\left(\parallel {\widehat{x}}_{t|t}^{\left(s+1\right)}-{\widehat{x}}_{t|t}^{\left(s\right)}\parallel /\parallel {\widehat{x}}_{t|t}^{\left(s\right)}\parallel \right)\le \varsigma$, the iteration stopped.

End for:

20. ${\widehat{x}}_{t|t}={\widehat{x}}_{t|t}^{\left(s\right)}$, $P_{t|t}=P_{t|t}^{\left(s\right)}$, $h_t=h_t^{\left(s\right)}$, $d_t=d_t^{\left(s\right)}$, ${\eta }_t={\eta }_t^{\left(s\right)}$, ${\delta }_t={\delta }_t^{\left(s\right)}$

Outputs: ${\widehat{x}}_{t|t}$, $P_{t|t}$, $h_t$, $d_t$, ${\eta }_t$, ${\delta }_t$, $h_t/\left(h_t+d_t\right)$, ${\eta }_t/\left({\eta }_t+{\delta }_t\right)$

describes the implementation of the proposed new KF.

4. Simulations

Aimed at demonstrating the superiority of the presented filter in the scenario with UTVMLP and NSHTMN, a numerical example is simulated. The process and measurement equations of the stochastic system are, respectively, given as [7]

$$x_t=\left[\begin{array}{cc}0.6& 0.4\\ 0.1& 0.9\end{array}\right]x_{t-1}+e_{t-1}$$

(63)

$$y_t=b_t\left[\begin{array}{cc}1& -2\end{array}\right]x_t+g_t$$

(64)

where the Gaussian process noise $e_{t-1}$ and the NSHTMN $g_t$ are given as [12]

$$e_{t-1}~\mathrm{N}\left(0,Q_t\right)$$

(65)

$$\{\begin{array}{c}{g}_t~\mathrm{N}\left(0,R_t\right) t\in \left[1,100\right]\left(\mathrm{Gaussian}\right)\\ g_t~\{\begin{array}{c}\mathrm{N}\left(0,R_t\right)\mathrm{w}.\mathrm{p}.=0.98\\ \mathrm{N}\left(0,500R_t\right)\mathrm{w}.\mathrm{p}.=0.02\end{array}t\in \left[101,200\right]\left(\mathrm{slightly} \mathrm{heavy}-\mathrm{tailed}\right)\\ g_t~\{\begin{array}{c}\mathrm{N}\left(0,R_t\right)\mathrm{w}.\mathrm{p}.=0.95\\ \mathrm{N}\left(0,500R_t\right)\mathrm{w}.\mathrm{p}.=0.05\end{array}t\in \left[201,300\right]\left(\mathrm{moderately} \mathrm{heavy}-\mathrm{tailed}\right)\\ g_t~ \mathrm{N}\left(0,R_t\right) t\in \left[301,400\right]\left(\mathrm{Gaussian}\right)\end{array}$$

(66)

where $\mathrm{w}.\mathrm{p}.$ represents “with probability”. The true process noise covariance matrix $Q_t$ with parameter $M=1$ and the nominal measurement noise covariance matrix $R_t$ with parameter $N=150 {\mathrm{m}}^2$ are set as

$$Q_t=M\left[\begin{array}{cc}1& 1\\ 0& 1\end{array}\right]$$

(67)

$$R_t=N$$

(68)

The real UTVMLP is set as

$$p\left({\gamma }_t\right)=\{\begin{array}{c}0.1 t\in \left[1,100\right]\\ 0.15 t\in \left[101,200\right]\\ 0.3 t\in \left[201,300\right]\\ 0.1 t\in \left[301,400\right]\end{array}$$

(69)

From Equations (66)–(69), it can be seen that the measurement noise and UTVMLP are divided into four stages, as shown in

Table 2. The stages of measurement and the corresponding measurement noise and UTVMLP.

Measurement Stage	Measurement Noise	UTVMLP
Stage 1, time 1 s~100s	$g_t~\mathrm{N}\left(0,R_t\right)$ (Gaussian)	0.1 (slight loss)
Stage 2, time 101 s~200 s	$g_t~\{\begin{array}{c}\mathrm{N}\left(0,R_t\right)\mathrm{w}.\mathrm{p}.=0.98\\ \mathrm{N}\left(0,500R_t\right)\mathrm{w}.\mathrm{p}.=0.02\end{array}$ (slightly heavy-tailed)	0.15 (slight loss)
Stage 3, time 201 s~300 s	$g_t~\{\begin{array}{c}\mathrm{N}\left(0,R_t\right)\mathrm{w}.\mathrm{p}.=0.95\\ \mathrm{N}\left(0,500R_t\right)\mathrm{w}.\mathrm{p}.=0.05\end{array}$ (moderately heavy-tailed)	0.3 (moderate loss)
Stage 4, time 301 s~400 s	$g_t~ \mathrm{N}\left(0,R_t\right)$ (Gaussian)	0.1 (slight loss)

. The remaining system parameters are as follows: the sampling interval $∆k=0.01 \mathrm{s}$ and the total simulation time $T=400 \mathrm{s}$. The proposed filter is compared with the typical Kalman filter (KF) [2]; the existing variational Bayesian-based adaptive KF with UTVMLP (VBAKF) [7]; the existing Gaussian-Student’s t-mixture distribution-based KF (GSTKF) with Gaussian process noise [14]; and the existing IKF with known real measurement loss probability [8]. The parameters of VBAKF are selected as $p_0=0.5$, ${\alpha }_0=5$, ${\beta }_0=5$, $\rho =1-\exp \left(-5\right)$, and $N_m=10$. The parameters of GSTKF are selected as $v_t=5$ and $e_0=0.85$. The parameters of the proposed filter are given as $\rho =0.99$, $\mu =5$, $h_0=0.85$, $N_I=10$, ${\eta }_0=5$, and ${\delta }_t=5$, $\varsigma ={10}^{-16}$. All filters are programmed with MATLAB R2018a and run on a computer with Intel Core i5-6300HQ CPU at 2.30 GHz and 8 GB of RAM.

Aimed at evaluating the performances in the estimation of the system state vector of all the algorithms, the root-mean square error (RMSE) and the averaged root-mean square error (AGRMSE) are utilized as performance indicators. The definitions of RMSE and AGRMSE of the system state are given as

$${\mathrm{RMSE}}_x=\sqrt{\frac{1}{M_c}{\displaystyle \sum }_{r=1}^{M_r}\left({\left(x_t^r-{\widehat{x}}_t^r\right)}^2+{\left(y_t^r-{\widehat{y}}_t^r\right)}^2\right)}$$

(70)

$${\mathrm{AGRMSE}}_{\mathrm{x}}=\sqrt{\frac{1}{M_{c}T}{\displaystyle \sum }_{t=1}^T{\displaystyle \sum }_{r=1}^{M_r}\left({\left(x_t^r-{\widehat{x}}_t^r\right)}^2+{\left(y_t^r-{\widehat{y}}_t^r\right)}^2\right)}$$

(71)

where $\left(x_t^r,y_t^r\right)$ and $\left({\widehat{x}}_t^r,{\widehat{y}}_t^r\right)$ denote the actual and estimated system state at the $j\mathrm{th}$ Monte Carlo run and discrete-time $t$, respectively. $M_c=500$ represents the total Monte Carlo run time.

Different from the proposed algorithm and VBAKF, the KF, IKF and GSTKF do not estimate UTVMLP. Although IKF can also address the filtering problem with measurement loss, IKF is based on the assumption that the measurement loss probability is known. Therefore, only VBAKF and the proposed algorithm participate in the comparison of the UTVMLP estimation performance.

Figure 1 shows the ${\mathrm{RMSE}}_x$s of the proposed filter and the existing filters over 500 times of the Monte Carlo run. Additionally, the ${\mathrm{AGRMSE}}_{\mathrm{x}}$s and SSRTs of different filters are listed in Table 2. It can be seen from Figure 1 and

Table 3. ${\mathrm{AGRMSE}}_{\mathrm{x}}$s and single-step running times (SSRT) of different filters.

Filters	KF	IKF	VBAKF	GSTKF	The Proposed Filter
${\mathrm{AGRMSE}}_{\mathrm{x}}$ in Stage 1	10.1336	4.4326	4.7705	8.7683	4.7704
${\mathrm{AGRMSE}}_{\mathrm{x}}$ in Stage 2	25.4009	11.1049	18.7603	14.4343	5.0660
${\mathrm{AGRMSE}}_{\mathrm{x}}$ in Stage 3	58.0561	17.5274	35.8194	28.0354	6.3044
${\mathrm{AGRMSE}}_{\mathrm{x}}$ in Stage 4	24.9129	4.4448	27.0252	9.5892	4.5081
${\mathrm{AGRMSE}}_{\mathrm{x}}$ in all stages	29.6259	9.3774	21.5944	15.2068	5.1623
SSRT (ms)	0.0276	0.0925	0.1166	0.1491	0.2989

that in the contexts of UTVMLP and NSHTMN, when the measurement is the Gaussian measurement noise and there is slight loss probability, as shown in stages 1 and 4, the estimation accuracy of the proposed filter is close to the IKF with true loss probability and the performance of the proposed algorithm is better than the other algorithms. We can also find that the proposed algorithm still has better performance than the existing algorithms when the measurement has heavy-tailed measurement noise and larger measurement loss probability, as shown in stages 2 and 4. In addition, the proposed algorithm has longer SSRT and higher computational complexity than the existing filters, which can be observed from Table 3.

Figure 2 shows the curves of the true and estimated UTVMLPs of VBAKF and the proposed filter over 500 times of the Monte Carlo run. Obviously, the NSHTMN has a great influence on the filtering performance of VBAKF and the proposed filter has better UTVMLP estimation accuracy than VBAKF in the scenario of NSHTMN.

Figure 3 and Figure 4 show the ${\mathrm{RMSE}}_x$s and the estimated UTVMLPs of the proposed filter with shape parameter $\mu =3,4,5,6$ over the 500 Monte Carlo run, respectively. The corresponding SSRTs of the proposed filter with $\mu =3,4,5,6$ are 0.2991, 0.2983, 0.2993, and 0.2989. It can be seen that the proposed filter with different shape parameters has better performance than current algorithms in the system state and UTVMLP estimations. Moreover, the degree of freedom parameter $\mu$ has little influence on the estimation accuracy and time complexity of the proposed algorithm, and the recommended value of $\mu$ is therefore set as 5.

Figure 5 and Figure 6 show the ${\mathrm{RMSE}}_x$s and the estimated UTVMLPs of the proposed filter with forgetting factor $\rho =0.93,0.95,0.97,0.99$ over the 500 Monte Carlo run, respectively. The corresponding SSRT of the proposed filter with $\rho =0.93,0.95,0.97,0.99$ is approximately equal to 0.2990. We can find that the proposed filter with $\rho =0.99$ has the best performance in the system state and UTVMLP estimations, and the value of $\rho$ has little effect on calculation complexity. Therefore, the recommended value of $\rho$ is 0.99.

Figure 7 shows the ${\mathrm{AGRMSE}}_{\mathrm{x}}$s of the proposed filter and the current algorithms with the iteration number $N_I=1,2,\cdots 10$. We can see from Figure 7 that the proposed filter has a smaller ${\mathrm{AGRMSE}}_{\mathrm{x}}$ than the existing filters when $N_I\ge 3$ and the proposed filter converges faster than the existing filters. However,

Table 4. SSRT of the proposed filter with the iteration number $N_I=1,2,\cdots 10$.

Iteration Number${\mathit{N}}_{\mathit{I}}$	1	2	3	4	5
SSRT (ms)	0.0327	0.0550	0.0894	0.1307	0.1570
Iteration Number${\mathit{N}}_{\mathit{I}}$	6	7	8	9	10
SSRT (ms)	0.1883	0.2243	0.2451	0.2835	0.2989

shows that the setting of $N_I$ has a huge impact on the time consumption of the proposed filter and the SSRT increases with the increase of $N_I$. Therefore, considering time consumption and estimation accuracy, the recommended value of $N_I$ ranges from 4 to 10.

5. Conclusions

In this paper, a new VB-based KF is presented to address the filtering issue with UTVMLP and NSHTMN. The system state vector, BM, BL, the intermediate random variables, the mixing probability, and the UTVMLP are simultaneously inferred by introducing the variational Bayesian technique. Simulation results illustrated that the proposed filter has a better performance than existing filters in the estimations of the system state vector and UTVMLP.

6. Future Work

The environmental factors in practical applications may be more complicated than this paper illustrated. Apart from non-stationary heavy-tailed measurement noise and unknown loss probability, the system process noise may also present non-stationary heavy-tailed distribution. In terms of measurement, random delays in measurement will also appear. Therefore, in our future research, we will further consider the factors of non-stationary heavy-tailed process noise and random measurement delay based on the theoretical content of this paper. Additionally, we will design a non-linear filtering method with lower computational complexity to verify the effectiveness in the real world.