# A Novel Robust Student’st-Based Cubature Information Filter with Heavy-Tailed Noises

Authors: Yongtao Shui; Xiaogang Wang; Wutao Qin; Yu Wang; Baojun Pang; Naigang Cui
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

In this paper, a novel robust Student’st-based cubature information filter is proposed for a nonlinear multisensor system with heavy-tailed process and measurement noises. At first, the predictive probability density function (PDF) and the likelihood PDF are approximated as two different Student’s t distributions. To avoid the process uncertainty induced by the heavy-tailed process noise, the scale matrix of the predictive PDF is modeled as an inverse Wishart distribution and estimated dynamically. Then, the predictive PDF and the likelihood PDF are transformed into a hierarchical Gaussian form to obtain the approximate solution of posterior PDF. Based on the variational Bayesian approximation method, the posterior PDF is approximated iteratively by minimizing the Kullback-Leibler divergence function. Based on the posterior PDF of the auxiliary parameters, the predicted covariance and measurement noise covariance are modified. And then the information matrix and information state are updated by summing the local information contributions, which are computed based on the modified covariance. Finally, the state, scale matrix, and posterior densities are estimated after fixed point iterations. And the simulation results for a target tracking example demonstrate the superiority of the proposed filter.


---

## Body

## 1. Introduction

To obtain the reliable and precise information, multisensor systems have become more and more popular in a wide range of applications such as cooperative tracking, autonomous navigation, signal processing, and guidance [1, 2]. Therefore, how to best fuse the observations from different sensors has gained much attention in several decades. Kalman filter (KF), one of the most popular state estimation methods, has been used extensively in many multisensor estimation problems [3, 4]. However, the KF is a linear state estimation method, and it cannot get an optimal estimation result with non-Gaussian process and measurement noises. But in many engineering applications, such as tracking agile target in clutter with bearing sensors, the system is nonlinear and the heavy-tailed noises are presented [5, 6]. In this case, the performance of the classical KF degrades dramatically.To apply the KF in nonlinear system and improve the accuracy of nonlinear approximation, the first order linearization-based extended Kalman filter (EKF) [7], the unscented transformation-based unscented Kalman filtering (UKF) [8, 9], quadrature rule-based Gauss-Hermite filtering (GHF) [10], polynomial interpolation-based divided difference filtering (DDF) [11], and the spherical-radical rule-based cubature Kalman filter (CKF) [12, 13] have been proposed successively. Although the information filter (IF) is an algebraically equivalent form of Kalman filter, it is more suitable for multisensor estimation, since the IF can reduce the computational burden significantly and be decentralized easily [14]. Soon, the abovementioned filters were formulated in the information type and applied to the multisensor estimation, such as extended information filter (EIF) [15], unscented information filter (UIF) [16], and cubature information filter (CIF) [17]. Unfortunately, the performance of these information filters will decrease in the presence of heavy-tailed noises, since they all derived based on the Gaussian assumption.To tackle the heavy-tailed non-Gaussian noise, some Huber-based information filters have been proposed by minimizing a cost function which is a combinedl1 and l2 norm [18], such as Huber-based UIF [19] and Huber-based CIF [20]. Nevertheless, the Huber’s weight function cannot converge to zero even as the error approaches infinity. It indicates the measurements with large errors still will be employed, and the estimation performance will be deteriorated [5, 21]. Besides Huber’s technique, some maximum correntropy (MC) information filters such as MC-UIF [22] and MC-CIF [5] have been proposed by employing the maximum correntropy criterion as the optimality criterion. However, the estimation covariance cannot be developed due to the lack of theoretical basis, which will degrade the performance of MC Kalman filters.Since the Student’st distribution is a reasonable model of the heavy-tailed noise, some Student’s t-based Kalman filters (STKF) have been proposed by modeling the heavy-tailed noises as the Student’s t distribution [23–27]. To obtain the approximate solution of posterior estimations, the posterior probability density function (PDF) has also been approximated by the Student’s t distribution [21, 23], in which the growth of the dof parameter in STKF degrades the estimation accuracy obviously. To tackle with the difficulty in choosing scale matrix and dof parameter, variational Bayesian (VB) approach is employed to jointly estimate the state vector, scale matrix, and dof parameter [24, 25]. In [28], the inaccurate process noise covariance noise and measurement noise covariance matrix is estimated by choosing inverse Wishart priors. In [29], both skewed noise and heavy-tailed noise are considered, whereas Gaussian scale mixture distributions are used to model one-step prediction and likelihood PDFs. Despite the excellent work mentioned above, there have been relatively few results concerning the nonlinear multisensor estimation problem with heavy-tailed process, and measurement noises have not been considered, which serves as the main motivation of this paper.To cope with the nonlinear multisensor estimation problem with heavy-tailed process and measurement noises, a novel robust Student’st-based cubature information filter (NRSTCIF) is proposed in this paper. At first, the one-step predicted state vector and nominal covariance is propagated based on the cubature rule. Then, the predictive PDF and the likelihood PDF are transformed into a hierarchical Gaussian form to obtain the approximate solution of posterior PDF. Next, the posterior densities are approximated iteratively based on the variational Bayesian (VB) approximation method, where the information matrix and information state are updated by summing the local information contributions. Finally, the state, scale matrix, and posterior densities are estimated after fixed point iterations.The rest of the paper is organized as follows. In Sec. II, the problem is formulated. Sec. III gives the derivation of the proposed filter. The simulations are conducted in Sec. IV, and the conclusions are drawn in Sec. V.


## 2. Problem Formulation

Consider a discrete-time stochastic nonlinear system described by(1)xk=fxk−1+wk−1(2)zk=hxk+vkwhere xk∈Rn is the n-dimensional state, zk∈Rm is the m-dimensional measurement, and wk∈Rn and vk∈Rm are zero mean heavily tailed process noise and measurement noise with nominal covariance Qk and Rk, respectively.Sincevk is a heavy-tailed measurement noise, its distribution is approximated as
(3)pvk=Stvk;0,Rk,υwhere Stvk;0,Rk,υ denotes the Student’s t probability distribution function (PDF) of vk, with mean 0 and scale matrix Rk, υ is the dof parameter.Based on the measurement equation and the distribution of measurement noise, the likelihood PDFpzkxk can be formulated as
(4)pzkxk=Stzk;hxk,Rk,υwhere hxk is the mean vector of zk.Due to the heavy-tailed process noise, the predictive PDFpxkZk−1 is also approximated as a Student’s t distribution.
(5)pxkZk−1=Stxk;ux^k∣k−1,Tk,μwhere x^k∣k−1 is the mean vector, μ is the dof parameter, and Tk is the scale matrix which will be approximated as the inverse Wishart distribution in the next section.


## 3. Derivation of the Proportional Filter

### 3.1. Brief Review of Extended Information Filter

Based on the state model given by (1), the predicted information matrix Yk∣k−1 and information state y^k∣k−1 can be formulated as [16, 17].
(6)Yk∣k−1=Pk∣k−1−1=ΦYk−1∣k−1−1ΦT+Qk−1−1y^k∣k−1=Yk∣k−1x^k∣k−1=Yk−1∣k−1fx^k−1∣k−1where Φ is the state-transition matrix which can be computed based on (1).Then, the updated information matrixYk∣k and information state y^k∣k can be formulated as
(7)y^k∣k=y^k∣k−1+ikYk∣k=Yk∣k−1+Ikwhere ik and Ik are information state contribution and information matrix contribution, respectively. They are given by
(8)ik=HkTRk−1zk−hx^k∣k−1+Hkx^k∣k−1(9)Ik=HkTRk−1Hkwhere zk is the measurement at time k and Hk is the measurement matrix.


### 3.2. The Transformation of the Predictive PDF and the Likelihood PDF

To get the approximate solution of posterior PDFpxkZk, the predictive PDF pxkZk−1 and the likelihood PDF pzkxk given by (5) and (4) are written as an infinite mixture of Gaussian distribution [24], i.e.,
(10)pxkZk−1=∫Nxk;x^k∣k−1,Tk/ςkGςk;μ2,μ2dςk(11)pzkxk=∫Nxk;hxk,Rk/ηkGηk;υ2,υ2dηkwhere ςk and ηk are auxiliary parameters, Gω;a,b is the Gamma PDF of auxiliary parameter ω. Gω;a,b is defined by
(12)Gω;a,b=baΓaωa−1e−bωwhere Γ⋅ denotes the Gamma function. Then, lnGω;a,b can be formulated as
(13)lnGω;a,b=lnωa−1+lne−bω+lnbaΓa=a−1lnω−bω+ca,bSince the process noise is a zero mean vector with nominal covarianceQk, the state vector and its associated nominal covariance can be propagated by
(14)x^k∣k−1=12n⋅∑i=12nfχi,k−1Pk∣k−1=12n⋅∑i=12nfχi,k−1−x^k∣k−1fχi,k−1−x∧k∣k−1T+Qk−1where χi,k−1i=1,2,⋯,2n is the i-th cubature point generated by
(15)χi,k−1=xk−1+Sk−1ξi,i=1,2,⋯,2nwhere Sk−1 is the Cholesky decomposition of Pk−1 satisfying Pk−1=Sk−1⋅Sk−1T, and ξi is the i-th element of the following set
(16)n10⋮0,01⋮0,⋯,00⋮1,−10⋮0,0−1⋮0,⋯,00⋮−1To accurately describe the predictive PDF, the VB approximation method is applied to dynamically estimate the scale matrixTk rather than set nominal covariance Pk∣k−1 as Tk. Then, the PDF of Tk can be appropriated as
(17)pTk=IWTk;tk,Τk=Τktk/2Tk−tk+n+1/2exp−0.5trΤkTk−12n⋅tk/2Γntk/2where IWTk;tk,Τk denotes the inverse Wishart distribution, whose functional form will not change between prior and posterior, tk is the dof parameter, and Τk is inverse scale matrix. Γn⋅ and tr⋅ denote the n-variate gamma function and the trace operation, respectively.Then,lnIWTk;tk,Τk can be formulated as
(18)lnIWTk;tk,Τk=lnTk−tk+n+1/2+lnexp−0.5trΤkTk−1+lnΤktk/2−ln2n⋅tk/2Γntk/2=−tk+n+12lnTk−trΤkTk−12+cTAccording to the property of inverse Wishart distribution, we have(19)ETk−1=tk−n−1Τk−1fortk>n+1To capture the statistics ofTk, the nominal covariance Pk∣k−1 is chosen as the mean value of pTk. Exploiting (19), the inverse scale matrix Τk is given by
(20)Τk=τPk∣k−1where the nonnegative tuning parameter τ=tk−n−1.Based on the (10), (11), and (17), the predictive PDF pxkZk−1 and the likelihood PDFpzkxk can be transformed as follows:
(21)pxkςk,Zk−1=Nxk;x^k∣k−1,Tk/ςkpTk=IWTk;tk,Τkpςk=Gςk;μ2,μ2(22)pzkxk,ηk=Nxk;hxk,Rk/ηkpηk=Gηk;υ2,υ2Note that we consider the case that the heavy-tailed process noise statistics, i.e., the scale matrixTk, is unknown and time varying, which is usually the case in practical applications especially for target tracking [29]. Therefore, the scale matrix Tk of process noise is modeled as inverse Wishart distribution and estimated dynamically. For measurement noise, we assume that the noise statistics Rk is known based on the prior information of sensors; therefore, it can be chosen wisely rather than estimated.


### 3.3. Derivation of Measurement Update Based on VB Approximation Method

To estimate the statexk, the joint posterior PDF pxk,ςk,Tk,ηkZk is required to be approximated. By making use of the VB approximation method, we have
(23)pxk,ςk,Tk,ηkZk≈qxkqςkqTkqηkwhere q⋅ is the approximate posterior PDF, which could be obtained by minimizing the following function.
(24)J=KLqxkqςkqTk−1qηkpxk,ςk,Tk−1,ηkZkwhere KL denotes the Kullback-Leibler divergence function.Minimizing (24) can reach the following VB-marginal [30].
(25)lnqφ=EΘ−φlnpΘ,Zk+cφwhere Θ=xk,ςk,Tk,ηk is a set consists of xk,ςk,Tk and ηk, φ is an arbitrary element of Θ and Θ−φ denotes the set Θ without φ, and cφ is the constant related to φ.Since the PDF given by (21) and (22) are conditional independent, pΘ,Zk could be factored as
(26)pΘ,Zk=pzkxk,ηkpxkςk,Zk−1pςkpTkpηkSubstituting (21) and (22) into (26), we have
(27)pΘ,Zk=Nxk;hxk,Rk/ηkNxk;x^k∣k−1,Tk/ςk⋅Gςk;μ2,μ2IWTk;tk,ΤkGηk;υ2,υ2Then,lnpΘ,Zk can be expressed as [24].
(28)lnpΘ,Zk=m+υ2−1lnηk−υ2ηk−ηk2zk−hxkTRk−1zk−hxk+n+μ2−1lnςk−μ2ςk−ςk2xk−x∧k∣k−1TTk−1xk−x^k∣k−1−trΤkTk−12−n+tk+22lnTkBy updating one element ofΘ and keeping the other fixed at their last estimated values, the VB-marginals given by (25) can be solved iteratively.Letϕ=ςk and substituting (28) in (25) results in
(29)lnqi+1ςk=n+μ2−1lnςk−μ+trPkiEiTk−12ςk+cςwhere the superscript “i” represents the i-th iteration. Pki is given by
(30)Pki=Eixk−x^k∣k−1xk−x∧k∣k−1T=Eixk−x^k∣kixk−x∧k∣kiT+x^k∣ki−x^k∣k−1x∧k∣ki−x∧k∣k−1T=Pk∣ki+x^k∣ki−x^k∣k−1x∧k∣ki−x∧k∣k−1TAnd the initial value of iteration is set asx^k∣k0=x^k∣k−1, Pk∣ki=Pk∣k−1, and E0Tk−1=tk−n−1Tk−1.Contrasting (13) with (29), the approximate posterior PDF qi+1ςk can be given by
(31)qi+1ςk=Gςk;αki+1,βki+1where αki+1=n+μ/2, βki+1=μ+trPkiEiTk−1/2.Letϕ=ηk and substituting (28) in (25) results in
(32)lnqi+1ηk=m+υ2−1lnηk−υ+trRkiRk−12ηk+cηwhere Rki is given by
(33)Rki=Eizk−hxkzk−hxkT=Eizk−hx^k∣ki+hx^k∣ki−hxk⋅⋅T=zk−hx^k∣kizk−hx∧k∣kiT+Pzzi=zk−hx^k∣kizk−hx∧k∣kiT+12n⋅∑j=12nhχj,k∣ki−z^k∣kihχj,k∣ki−z∧k∣kiTThei-th cubature pointχj,k∣kii=1,2,⋯,2n and the predicted measurement z^k∣ki are computed as follows:
(34)χj,k∣ki=x^k∣ki+Sk∣kiξj,j=1,2,⋯,2nz^k∣ki=12n⋅∑j=12nhχj,k∣kiwhere Sk∣ki is the Cholesky decomposition of Pk∣ki satisfying Pk∣ki=Sk∣ki⋅Sk∣kiT.Contrasting (13) with (32), the approximate posterior PDF qi+1ηk can be given by
(35)qi+1ηk=Gηk;ϖki+1,ζki+1where ϖki+1=m+υ/2, ζki+1=υ+trRkiRk−1/2Letϕ=Tk and substituting (28) in (25) results in
(36)lnqi+1Tk=−n+tk+22lnTk−trΤk+Ei+1ςkPkiTk−12+cTContrasting (36) with (18), the approximate posterior PDF qi+1Tk can be given by
(37)qi+1Tk=IWTk;tki+1,Tki+1where tki+1=tk+1, Tki+1=Τk+Ei+1ςkPki, and Ei+1ςk=αki+1/βki+1.Letφ=xk and substituting (28) into (25) results in
(38)lnqi+1xk=−Ei+1ηkzk−hxkTRk−1zk−hxk2−Ei+1ςkxk−x∧k∣k−1TEi+1Tk−1xk−x^k∣k−12+cxwhere Ei+1ηk=ϖki+1/ζki+1, Ei+1Tk−1=tki+1−n−1Tki+1−1.According to the Bayesian theory, the posterior PDF can be approximated by(39)qi+1xk=pi+1zkxkpi+1xkZk−1∫pi+1zkxkpi+1xkZk−1dxkExploiting (38), pi+1zkxk and pi+1xkZk−1 can be written as
(40)pi+1zkxk=Nzk;hxk,R¯ki+1pi+1xkZk−1=Nxk;x^k∣k−1,P¯k∣k−1i+1where R¯ki+1=Rk/Ei+1ηk, P¯k∣k−1i+1=Ei+1Tk−1−1/Ei+1ςk.Since Gaussian distribution can approximate the posterior PDF more accurate than the Student’st distribution [24], the approximate posterior PDF qi+1xk can be formulated by
(41)qi+1xk=Nxk;x^k∣ki+1,Pk∣ki+1According to the basic equations of information filter, we have(42)Y¯k∣k−1i+1=P¯k∣k−1i+1−1=Ei+1ςkEi+1Tk−1yk∣k−1i+1=Y¯k∣k−1i+1x^k∣k−1Based on the linear error propagation, the cross covarianceP¯xz,k∣k−1i+1 is approximated by
(43)P¯xz,k∣k−1i+1≈P¯k∣k−1i+1Hki+1MultiplyingP¯k∣k−1i+1−1P¯k∣k−1i+1 and P¯k∣k−1i+1−1P¯k∣k−1i+1T on the right hand side of (8) and (9) yields
(44)iki+1=P¯k∣k−1i+1−1P¯k∣k−1i+1Hki+1TRk−1⋅zk−z^k∣k−1i+1+Hki+1P¯k∣k−1i+1TP¯k∣k−1i+1−Tx^k∣k−1(45)Iki+1=P¯k∣k−1i+1−1P¯k∣k−1i+1HkTR¯ki+1−1HkP¯k∣k−1i+1TP¯k∣k−1i+1−TUsing (43) to simplify (44) and (45) results in
(46)iki+1=P¯k∣k−1i+1−1P¯xz,k∣k−1i+1Rk−1⋅zk−z^k∣k−1i+1+P¯xz,k∣k−1i+1TP¯k∣k−1i+1−Tx^k∣k−1(47)Iki+1=P¯k∣k−1i+1−1P¯xz,k∣k−1i+1R¯ki+1−1P¯xz,k∣k−1i+1TP¯k∣k−1i+1−TThe predicted measurementz^k∣k−1i+1 and cross covariance P¯xz,k∣k−1i+1 are computed as follows:
(48)z^k∣k−1i+1=12n⋅∑i=12nhχj,k∣k−1i+1P¯xz,k∣k−1i+1=12n⋅∑i=12nχj,k∣k−1i+1−x^k∣k−1hχj,k∣k−1i+1−z∧k∣k−1i+1Twhere χj,k∣k−1i+1i=1,2,⋯,2n is the i-th cubature point generated by
(49)χj,k∣k−1i+1=x^k∣k−1+S¯k∣k−1i+1ξj,j=1,2,⋯,2nand S¯k∣k−1i+1 is the Cholesky decomposition of P¯k∣k−1i+1.Based on the information contributions given by (46) and(47), the information state and matrix are updated as
(50)y^k∣ki+1=yk∣k−1i+1+iki+1Yk∣ki+1=Y¯k∣k−1i+1+Iki+1Then, the estimated state and covariance can be recovered as(51)Pk∣ki+1=Yk∣ki+1−1x^k∣ki+1=Pk∣ki+1y^k∣ki+1After the iterations, the approximate PDFs given in (23) are updated as
(52)qxk≈Nxk;x^k∣kd,Pk∣kdqςk≈Gςk;αk∣kd,βk∣kdqηk≈Gηk;ϖk∣kd,ζk∣kdqTk=IWTk;tkd,Tkdwhere d is the number of iterations (NOI).The main advantage of the IF is its ability to fuse multisource measurements simply by adding the information contributions to the information matrix and state. SupposeN sensors is available, the measurement equations can be given by
(53)zκ,k=hκxk+vκ,k,κ=1,2,⋯,Nwhere the subscript κ denotes the κ-th sensor, vκ,k~Stvκ,k;0,Rκ,k,υ is the measurement noise.Then, the information matrix and information state can be updated by(54)y^k∣ki+1=yk∣k−1i+1+∑κ=1Niκ,ki+1Yk∣ki+1=Y¯k∣k−1i+1+∑κ=1NIκ,ki+1where the information contributions of κ-th sensor are computed by
(55)iκ,ki+1=P¯k∣k−1i+1−1P¯κ,xz,k∣k−1i+1R¯k−1zκ,k−z^κ,k∣k−1i+1+P¯κ,xz,k∣k−1i+1TP¯k∣k−1i+1−Tx^k∣k−1,κ=1,2,⋯,NIκ,ki+1=P¯k∣k−1i+1−1P¯κ,xz,k∣k−1i+1R¯κ,ki+1−1⋅P¯κ,xz,k∣k−1i+1TP¯k∣k−1i+1−T,κ=1,2,⋯,NTo better illustrate the proposed filtering algorithm, the computational procedures are summarized in Figure1.Figure 1
Diagram of the proposed NRSTCIF.


## 3.1. Brief Review of Extended Information Filter

Based on the state model given by (1), the predicted information matrix Yk∣k−1 and information state y^k∣k−1 can be formulated as [16, 17].
(6)Yk∣k−1=Pk∣k−1−1=ΦYk−1∣k−1−1ΦT+Qk−1−1y^k∣k−1=Yk∣k−1x^k∣k−1=Yk−1∣k−1fx^k−1∣k−1where Φ is the state-transition matrix which can be computed based on (1).Then, the updated information matrixYk∣k and information state y^k∣k can be formulated as
(7)y^k∣k=y^k∣k−1+ikYk∣k=Yk∣k−1+Ikwhere ik and Ik are information state contribution and information matrix contribution, respectively. They are given by
(8)ik=HkTRk−1zk−hx^k∣k−1+Hkx^k∣k−1(9)Ik=HkTRk−1Hkwhere zk is the measurement at time k and Hk is the measurement matrix.


## 3.2. The Transformation of the Predictive PDF and the Likelihood PDF

To get the approximate solution of posterior PDFpxkZk, the predictive PDF pxkZk−1 and the likelihood PDF pzkxk given by (5) and (4) are written as an infinite mixture of Gaussian distribution [24], i.e.,
(10)pxkZk−1=∫Nxk;x^k∣k−1,Tk/ςkGςk;μ2,μ2dςk(11)pzkxk=∫Nxk;hxk,Rk/ηkGηk;υ2,υ2dηkwhere ςk and ηk are auxiliary parameters, Gω;a,b is the Gamma PDF of auxiliary parameter ω. Gω;a,b is defined by
(12)Gω;a,b=baΓaωa−1e−bωwhere Γ⋅ denotes the Gamma function. Then, lnGω;a,b can be formulated as
(13)lnGω;a,b=lnωa−1+lne−bω+lnbaΓa=a−1lnω−bω+ca,bSince the process noise is a zero mean vector with nominal covarianceQk, the state vector and its associated nominal covariance can be propagated by
(14)x^k∣k−1=12n⋅∑i=12nfχi,k−1Pk∣k−1=12n⋅∑i=12nfχi,k−1−x^k∣k−1fχi,k−1−x∧k∣k−1T+Qk−1where χi,k−1i=1,2,⋯,2n is the i-th cubature point generated by
(15)χi,k−1=xk−1+Sk−1ξi,i=1,2,⋯,2nwhere Sk−1 is the Cholesky decomposition of Pk−1 satisfying Pk−1=Sk−1⋅Sk−1T, and ξi is the i-th element of the following set
(16)n10⋮0,01⋮0,⋯,00⋮1,−10⋮0,0−1⋮0,⋯,00⋮−1To accurately describe the predictive PDF, the VB approximation method is applied to dynamically estimate the scale matrixTk rather than set nominal covariance Pk∣k−1 as Tk. Then, the PDF of Tk can be appropriated as
(17)pTk=IWTk;tk,Τk=Τktk/2Tk−tk+n+1/2exp−0.5trΤkTk−12n⋅tk/2Γntk/2where IWTk;tk,Τk denotes the inverse Wishart distribution, whose functional form will not change between prior and posterior, tk is the dof parameter, and Τk is inverse scale matrix. Γn⋅ and tr⋅ denote the n-variate gamma function and the trace operation, respectively.Then,lnIWTk;tk,Τk can be formulated as
(18)lnIWTk;tk,Τk=lnTk−tk+n+1/2+lnexp−0.5trΤkTk−1+lnΤktk/2−ln2n⋅tk/2Γntk/2=−tk+n+12lnTk−trΤkTk−12+cTAccording to the property of inverse Wishart distribution, we have(19)ETk−1=tk−n−1Τk−1fortk>n+1To capture the statistics ofTk, the nominal covariance Pk∣k−1 is chosen as the mean value of pTk. Exploiting (19), the inverse scale matrix Τk is given by
(20)Τk=τPk∣k−1where the nonnegative tuning parameter τ=tk−n−1.Based on the (10), (11), and (17), the predictive PDF pxkZk−1 and the likelihood PDFpzkxk can be transformed as follows:
(21)pxkςk,Zk−1=Nxk;x^k∣k−1,Tk/ςkpTk=IWTk;tk,Τkpςk=Gςk;μ2,μ2(22)pzkxk,ηk=Nxk;hxk,Rk/ηkpηk=Gηk;υ2,υ2Note that we consider the case that the heavy-tailed process noise statistics, i.e., the scale matrixTk, is unknown and time varying, which is usually the case in practical applications especially for target tracking [29]. Therefore, the scale matrix Tk of process noise is modeled as inverse Wishart distribution and estimated dynamically. For measurement noise, we assume that the noise statistics Rk is known based on the prior information of sensors; therefore, it can be chosen wisely rather than estimated.


## 3.3. Derivation of Measurement Update Based on VB Approximation Method

To estimate the statexk, the joint posterior PDF pxk,ςk,Tk,ηkZk is required to be approximated. By making use of the VB approximation method, we have
(23)pxk,ςk,Tk,ηkZk≈qxkqςkqTkqηkwhere q⋅ is the approximate posterior PDF, which could be obtained by minimizing the following function.
(24)J=KLqxkqςkqTk−1qηkpxk,ςk,Tk−1,ηkZkwhere KL denotes the Kullback-Leibler divergence function.Minimizing (24) can reach the following VB-marginal [30].
(25)lnqφ=EΘ−φlnpΘ,Zk+cφwhere Θ=xk,ςk,Tk,ηk is a set consists of xk,ςk,Tk and ηk, φ is an arbitrary element of Θ and Θ−φ denotes the set Θ without φ, and cφ is the constant related to φ.Since the PDF given by (21) and (22) are conditional independent, pΘ,Zk could be factored as
(26)pΘ,Zk=pzkxk,ηkpxkςk,Zk−1pςkpTkpηkSubstituting (21) and (22) into (26), we have
(27)pΘ,Zk=Nxk;hxk,Rk/ηkNxk;x^k∣k−1,Tk/ςk⋅Gςk;μ2,μ2IWTk;tk,ΤkGηk;υ2,υ2Then,lnpΘ,Zk can be expressed as [24].
(28)lnpΘ,Zk=m+υ2−1lnηk−υ2ηk−ηk2zk−hxkTRk−1zk−hxk+n+μ2−1lnςk−μ2ςk−ςk2xk−x∧k∣k−1TTk−1xk−x^k∣k−1−trΤkTk−12−n+tk+22lnTkBy updating one element ofΘ and keeping the other fixed at their last estimated values, the VB-marginals given by (25) can be solved iteratively.Letϕ=ςk and substituting (28) in (25) results in
(29)lnqi+1ςk=n+μ2−1lnςk−μ+trPkiEiTk−12ςk+cςwhere the superscript “i” represents the i-th iteration. Pki is given by
(30)Pki=Eixk−x^k∣k−1xk−x∧k∣k−1T=Eixk−x^k∣kixk−x∧k∣kiT+x^k∣ki−x^k∣k−1x∧k∣ki−x∧k∣k−1T=Pk∣ki+x^k∣ki−x^k∣k−1x∧k∣ki−x∧k∣k−1TAnd the initial value of iteration is set asx^k∣k0=x^k∣k−1, Pk∣ki=Pk∣k−1, and E0Tk−1=tk−n−1Tk−1.Contrasting (13) with (29), the approximate posterior PDF qi+1ςk can be given by
(31)qi+1ςk=Gςk;αki+1,βki+1where αki+1=n+μ/2, βki+1=μ+trPkiEiTk−1/2.Letϕ=ηk and substituting (28) in (25) results in
(32)lnqi+1ηk=m+υ2−1lnηk−υ+trRkiRk−12ηk+cηwhere Rki is given by
(33)Rki=Eizk−hxkzk−hxkT=Eizk−hx^k∣ki+hx^k∣ki−hxk⋅⋅T=zk−hx^k∣kizk−hx∧k∣kiT+Pzzi=zk−hx^k∣kizk−hx∧k∣kiT+12n⋅∑j=12nhχj,k∣ki−z^k∣kihχj,k∣ki−z∧k∣kiTThei-th cubature pointχj,k∣kii=1,2,⋯,2n and the predicted measurement z^k∣ki are computed as follows:
(34)χj,k∣ki=x^k∣ki+Sk∣kiξj,j=1,2,⋯,2nz^k∣ki=12n⋅∑j=12nhχj,k∣kiwhere Sk∣ki is the Cholesky decomposition of Pk∣ki satisfying Pk∣ki=Sk∣ki⋅Sk∣kiT.Contrasting (13) with (32), the approximate posterior PDF qi+1ηk can be given by
(35)qi+1ηk=Gηk;ϖki+1,ζki+1where ϖki+1=m+υ/2, ζki+1=υ+trRkiRk−1/2Letϕ=Tk and substituting (28) in (25) results in
(36)lnqi+1Tk=−n+tk+22lnTk−trΤk+Ei+1ςkPkiTk−12+cTContrasting (36) with (18), the approximate posterior PDF qi+1Tk can be given by
(37)qi+1Tk=IWTk;tki+1,Tki+1where tki+1=tk+1, Tki+1=Τk+Ei+1ςkPki, and Ei+1ςk=αki+1/βki+1.Letφ=xk and substituting (28) into (25) results in
(38)lnqi+1xk=−Ei+1ηkzk−hxkTRk−1zk−hxk2−Ei+1ςkxk−x∧k∣k−1TEi+1Tk−1xk−x^k∣k−12+cxwhere Ei+1ηk=ϖki+1/ζki+1, Ei+1Tk−1=tki+1−n−1Tki+1−1.According to the Bayesian theory, the posterior PDF can be approximated by(39)qi+1xk=pi+1zkxkpi+1xkZk−1∫pi+1zkxkpi+1xkZk−1dxkExploiting (38), pi+1zkxk and pi+1xkZk−1 can be written as
(40)pi+1zkxk=Nzk;hxk,R¯ki+1pi+1xkZk−1=Nxk;x^k∣k−1,P¯k∣k−1i+1where R¯ki+1=Rk/Ei+1ηk, P¯k∣k−1i+1=Ei+1Tk−1−1/Ei+1ςk.Since Gaussian distribution can approximate the posterior PDF more accurate than the Student’st distribution [24], the approximate posterior PDF qi+1xk can be formulated by
(41)qi+1xk=Nxk;x^k∣ki+1,Pk∣ki+1According to the basic equations of information filter, we have(42)Y¯k∣k−1i+1=P¯k∣k−1i+1−1=Ei+1ςkEi+1Tk−1yk∣k−1i+1=Y¯k∣k−1i+1x^k∣k−1Based on the linear error propagation, the cross covarianceP¯xz,k∣k−1i+1 is approximated by
(43)P¯xz,k∣k−1i+1≈P¯k∣k−1i+1Hki+1MultiplyingP¯k∣k−1i+1−1P¯k∣k−1i+1 and P¯k∣k−1i+1−1P¯k∣k−1i+1T on the right hand side of (8) and (9) yields
(44)iki+1=P¯k∣k−1i+1−1P¯k∣k−1i+1Hki+1TRk−1⋅zk−z^k∣k−1i+1+Hki+1P¯k∣k−1i+1TP¯k∣k−1i+1−Tx^k∣k−1(45)Iki+1=P¯k∣k−1i+1−1P¯k∣k−1i+1HkTR¯ki+1−1HkP¯k∣k−1i+1TP¯k∣k−1i+1−TUsing (43) to simplify (44) and (45) results in
(46)iki+1=P¯k∣k−1i+1−1P¯xz,k∣k−1i+1Rk−1⋅zk−z^k∣k−1i+1+P¯xz,k∣k−1i+1TP¯k∣k−1i+1−Tx^k∣k−1(47)Iki+1=P¯k∣k−1i+1−1P¯xz,k∣k−1i+1R¯ki+1−1P¯xz,k∣k−1i+1TP¯k∣k−1i+1−TThe predicted measurementz^k∣k−1i+1 and cross covariance P¯xz,k∣k−1i+1 are computed as follows:
(48)z^k∣k−1i+1=12n⋅∑i=12nhχj,k∣k−1i+1P¯xz,k∣k−1i+1=12n⋅∑i=12nχj,k∣k−1i+1−x^k∣k−1hχj,k∣k−1i+1−z∧k∣k−1i+1Twhere χj,k∣k−1i+1i=1,2,⋯,2n is the i-th cubature point generated by
(49)χj,k∣k−1i+1=x^k∣k−1+S¯k∣k−1i+1ξj,j=1,2,⋯,2nand S¯k∣k−1i+1 is the Cholesky decomposition of P¯k∣k−1i+1.Based on the information contributions given by (46) and(47), the information state and matrix are updated as
(50)y^k∣ki+1=yk∣k−1i+1+iki+1Yk∣ki+1=Y¯k∣k−1i+1+Iki+1Then, the estimated state and covariance can be recovered as(51)Pk∣ki+1=Yk∣ki+1−1x^k∣ki+1=Pk∣ki+1y^k∣ki+1After the iterations, the approximate PDFs given in (23) are updated as
(52)qxk≈Nxk;x^k∣kd,Pk∣kdqςk≈Gςk;αk∣kd,βk∣kdqηk≈Gηk;ϖk∣kd,ζk∣kdqTk=IWTk;tkd,Tkdwhere d is the number of iterations (NOI).The main advantage of the IF is its ability to fuse multisource measurements simply by adding the information contributions to the information matrix and state. SupposeN sensors is available, the measurement equations can be given by
(53)zκ,k=hκxk+vκ,k,κ=1,2,⋯,Nwhere the subscript κ denotes the κ-th sensor, vκ,k~Stvκ,k;0,Rκ,k,υ is the measurement noise.Then, the information matrix and information state can be updated by(54)y^k∣ki+1=yk∣k−1i+1+∑κ=1Niκ,ki+1Yk∣ki+1=Y¯k∣k−1i+1+∑κ=1NIκ,ki+1where the information contributions of κ-th sensor are computed by
(55)iκ,ki+1=P¯k∣k−1i+1−1P¯κ,xz,k∣k−1i+1R¯k−1zκ,k−z^κ,k∣k−1i+1+P¯κ,xz,k∣k−1i+1TP¯k∣k−1i+1−Tx^k∣k−1,κ=1,2,⋯,NIκ,ki+1=P¯k∣k−1i+1−1P¯κ,xz,k∣k−1i+1R¯κ,ki+1−1⋅P¯κ,xz,k∣k−1i+1TP¯k∣k−1i+1−T,κ=1,2,⋯,NTo better illustrate the proposed filtering algorithm, the computational procedures are summarized in Figure1.Figure 1
Diagram of the proposed NRSTCIF.


## 4. Simulation

To demonstrate the feasibility of the proposed algorithm, a target tracking problem is investigated in this section. The target moves in a plane with an unknown speed and unknown constant turn rate, it is observed by multiple radars in clutter. Due to the rapid motion and clutter, the heavy-tailed noises are introduced.The state equation is given by(56)xk=1sinωk−1Tωk−10cosωk−1T−1ωk−100cosωk−1T0−sinωk−1T001−cosωk−1Tωk−11sinωk−1Tωk−100sinωk−1T0cosωk−1T000001xk+wk−1where the state vector xkis defined as xk=xk,ẋk,yk,ẏk,ωkT, ωk is the turn rate, T=1s is the sampling interval, and wk−1 is the process noise generated by
(57)wk−1=N0,Qk−1withprobability1−pN0,50Qk−1withpoobabilitypwhere p=0.05 is probability of outlier, Qk−1 is given by
(58)Qk−1=T3/3T2/2000T2/2T00000T3/3T2/2000T2/20000001.75×10−3TThe measurement equation ofj-th radar is given by
(59)zj,k=xk−xj,s2+yk−yj,s2atan2yk−yj,s,xk−xj,s+vkwhere xj,s,yj,s is the position of j-th radar, vk is the measurement noise generated by
(60)vk=N0,Rkwithprobability1−pN0,100Rkwithprobabilitypand Rk is given by
(61)Rk=100m200100mrad2The locations of the four sensors are given as(62)x1,s=−30000my1,s=−5000mx2,s=0my2,s=35000mx3,s=30000my3,s=−5000mx4,s=0my4,s=−40000mThe initial state estimationx^0∣0 is generated randomly based on the true initial value x0 and initial covariance P0 given by
(63)x0=1000m,300m/s,1000m,0,−3∘/sTP0=diag100m2,10m2/s2,100m2,10m2/s2,100mrad2/s2TIn the simulation, the dof parameter, tuning parameter, and NOI of the proposed filter are set asμ=υ=5, τ=5, and d=10, respectively. The root mean square error (RMSE) is chosen as the performance metric. Based on 1000 independent Monte Carlo simulation runs, Figure 2 compares the performance of cubature information filter (CIF) [17], Huber-based cubature information filter (HCIF) [19], Student’s t-based cubature information filter (STCIF) [23], and the proposed filter. It can be seen that the CIF has the biggest RMSEs among four filters. This is not surprise, since the CIF is essentially a minimizer of a weighted least-squares criterion, which is sensitive to heavy-tailed noises. In addition, one can also observe that the RMSEs of the proposed filter are smaller than the existing HCIF and STCIF, which indicates the superiority of the proposed filter.Figure 2
Tracking performance of four filters; (a) Position estimation error, (b) Velocity estimation error, (c) Turn rate estimation error.
(a)(b)(c)The average computational costs of the various filters for each Monte-Carlo run are given in Table1. It can be seen that the computational time of the NRSTCIF is the longest, and that of the CIF is the shortest. It takes the NRSTCIF about 5 times as long as the CIF. Therefore, the algorithm proposed in this paper does not increase the computation time dramatically.Table 1
Computational time comparisons.
FilterCIFHCIFSTCIFNRSTCIFTime(s)0.110.190.240.5To clarify the relation between estimation accuracy and NOI, another 1000 independent Monte Carlo simulation runs with different NOIs are conducted based on the proposed filter, and the simulation results are summarized in Figure3. It can be seen that the averaged RMSEs decrease as the NOI increase when NOI=1,2,and3, but the averaged RMSEs are almost same when NOI ≥4. Therefore, it is unnecessary to choose the NOI as big as possible, and the result is an important reference for NOI choosing.Figure 3
Averaged RMSEs with different NOIs.Figure4 shows the RMSEs of the CIF, HCIF, and the proposed filter with different probabilities of outlier p=0.01,0.02,⋯,0.1. It can be seen the proposed filter is superior to the CIF and HCIF for all values of p. And the averaged RMSEs of CIF and HCIF increase more rapidly than the proposed NRSTCIF as the increase of p, which indicate the superiority of the NRSTCIF is more and more significant. In other words, the proposed NRSTCIF is more promising when the probability of outlier is larger.Figure 4
Performance comparison of three filters with different probabilities of outlier. (a) Position estimation error. (b) Velocity estimation error. (c) Turn rate estimation error.
(a)(b)(c)


## 5. Conclusions

In this paper, a novel robust Student’st-based cubature information filter is proposed to handle the heavy-tailed process and measurement noises in nonlinear multisensor estimation. Based on the cubature rule, the one-step predicted state and the nominal covariance are propagated. To get the approximate solution of the posterior PDF, the predictive PDF and the likelihood PDF described by Student’s t distribution are transformed into a hierarchical Gaussian form. To approximate the joint posterior PDF, the VB approximation method is utilized. After fixed point iterations, the state, scale matrix, and posterior densities are estimated. And the simulation results demonstrate the superiority of the proposed filter.



---

*Source: 7075037-2020-07-17.xml*