Robust Student's t based Stochastic Cubature Filter for Nonlinear Systems with Heavy-tailed Process and Measurement Noises

TL;DR

This paper introduces a robust Student's t based stochastic cubature filter for nonlinear systems with heavy-tailed noises, improving estimation accuracy and computational efficiency over existing methods.

Contribution

A novel stochastic Student's t spherical radial cubature rule is developed, enabling a robust filter that outperforms existing approaches in accuracy and efficiency.

Findings

Achieves higher estimation accuracy than existing filters.

Demonstrates robustness in heavy-tailed noise environments.

Offers computational efficiency superior to particle filters.

Abstract

In this paper, a new robust Student's t based stochastic cubature filter (RSTSCF) is proposed for nonlinear state-space model with heavy-tailed process and measurement noises. The heart of the RSTSCF is a stochastic Student's t spherical radial cubature rule (SSTSRCR), which is derived based on the third-degree unbiased spherical rule and the proposed third-degree unbiased radial rule. The existing stochastic integration rule is a special case of the proposed SSTSRCR when the degrees of freedom parameter tends to infinity. The proposed filter is applied to a manoeuvring bearings-only tracking example, where an agile target is tracked and the bearing is observed in clutter. Simulation results show that the proposed RSTSCF can achieve higher estimation accuracy than the existing Gaussian approximate filter, Gaussian sum filter, Huber-based nonlinear Kalman filter, maximum correntropy criterion based Kalman filter, and robust Student's t based nonlinear filters, and is computationally much more efficient than the existing particle filter.

Tables

Table 1. TABLE I: The implementation pseudocode of the proposed SSTSRCR.

Inputs: $\mu$, $𝚺$, $\nu$, $𝐠(\cdot )$, $n$, $N$.

1. Initialization: $I_s^3\left[𝐠\right]=0$.

for $l=1:N$

2. Generate a $n\times n$ matrix ${𝐔}_l$ of standard norm variables.

3. Obtain the required random orthogonal matrix ${𝐐}_l$ using the

$QR$ factorization: ${𝐔}_l={𝐐}_l{𝐑}_l$.

4. Draw the random variable ${\tau }_l$ from a Beta distribution:

${\tau }_l\sim \mathrm{Beta}(\frac{n+2}{2},\frac{\nu -2}{2})$.

5. Calculate the random quadrature point $r_{2,l}$:

$r_{2,l}=\sqrt{\frac{{\tau }_l}{1-{\tau }_l}}$.

6. Update $I_s^3\left[𝐠\right]$ at current iteration:

$I_s^3\left[𝐠\right]=I_s^3\left[𝐠\right]+\frac{1}{N}\{[1-\frac{n}{\left(\nu -2\right)r_{2,l}^2}]𝐠\left(\mu \right)+\frac{1}{2\left(\nu -2\right)r_{2,l}^2}\times$

$\sum _{i=1}^n[𝐠(\mu -r_{2,l}\sqrt{\nu 𝚺}{𝐐}_l{𝐞}_i)+𝐠(\mu +r_{2,l}\sqrt{\nu 𝚺}{𝐐}_l{𝐞}_i)]\}$.

end for

Outputs: $I\left[𝐠\right]\approx I_s^3\left[𝐠\right]$.

Table 2. TABLE II: The implementation pseudocode for one time step of the proposed RSTSCF.

Inputs: ${\widehat{𝐱}}_{k-1|k-1}$, ${𝐏}_{k-1|k-1}$, ${𝐳}_k$, ${𝐐}_{k-1}$, ${𝐑}_k$, ${\nu }_1$, ${\nu }_2$, ${\nu }_3$, ${𝐟}_{k-1}(\cdot )$,

${𝐡}_k(\cdot )$, $n$, $N$.

Time update:

1. ${\widehat{𝐱}}_{k|k-1}=\mathrm{SSTSRCR}({\widehat{𝐱}}_{k-1|k-1},{𝐏}_{k-1|k-1},{\nu }_3,{𝐟}_{k-1}(\cdot ),n,N)$.

2. ${𝐏}_{k|k-1}=\frac{{\nu }_3-2}{{\nu }_3}\mathrm{SSTSRCR}({\widehat{𝐱}}_{k-1|k-1},{𝐏}_{k-1|k-1},{\nu }_3,$

${𝐟}_{k-1}(\cdot ){𝐟}_{k-1}^T(\cdot ),n,N)-\frac{{\nu }_3-2}{{\nu }_3}\widehat{𝐱}{}_{k|k-1}\widehat{𝐱}{}_{k|k-1}{}^T+\frac{{\nu }_1\left({\nu }_3-2\right)}{\left({\nu }_1-2\right){\nu }_3}𝐐{}_{k-1}$.

Measurement update:

3. ${\widehat{𝐳}}_{k|k-1}=\mathrm{SSTSRCR}({\widehat{𝐱}}_{k|k-1},{𝐏}_{k|k-1},{\nu }_3,{𝐡}_k(\cdot ),n,N)$.

4. ${𝐏}_{k|k-1}^{zz}=\frac{{\nu }_3-2}{{\nu }_3}\mathrm{SSTSRCR}({\widehat{𝐱}}_{k|k-1},{𝐏}_{k|k-1},{\nu }_3,{𝐡}_k(\cdot ){𝐡}_k^T(\cdot ),$

$n,N)-\frac{{\nu }_3-2}{{\nu }_3}\widehat{𝐳}{}_{k|k-1}\widehat{𝐳}{}_{k|k-1}{}^T+\frac{{\nu }_2\left({\nu }_3-2\right)}{\left({\nu }_2-2\right){\nu }_3}𝐑{}_k$.

5. ${𝐏}_{k|k-1}^{xz}=\frac{{\nu }_3-2}{{\nu }_3}\mathrm{SSTSRCR}({\widehat{𝐱}}_{k|k-1},{𝐏}_{k|k-1},{\nu }_3,{𝐱}_k{𝐡}_k(\cdot ),n,N)$

$-\frac{{\nu }_3-2}{{\nu }_3}{\widehat{𝐱}}_{k|k-1}{\widehat{𝐳}}_{k|k-1}^T$.

6. ${\Delta }_k=\sqrt{{\left({𝐳}_k-{\widehat{𝐳}}_{k|k-1}\right)}^T{\left({𝐏}_{k|k-1}^{zz}\right)}^{-1}\left({𝐳}_k-{\widehat{𝐳}}_{k|k-1}\right)}$.

7. ${𝐊}_k={𝐏}_{k|k-1}^{xz}{\left({𝐏}_{k|k-1}^{zz}\right)}^{-1}$.

8. ${\widehat{𝐱}}_{k|k}={\widehat{𝐱}}_{k|k-1}+{𝐊}_k\left({𝐳}_k-{\widehat{𝐳}}_{k|k-1}\right)$.

9. ${𝐏}_{k|k}=\frac{\left({\nu }_3-2\right)\left({\nu }_3+{\Delta }_k^2\right)}{{\nu }_3\left({\nu }_3+m-2\right)}\left({𝐏}_{k|k-1}-{𝐊}_k{𝐏}_{k|k-1}^{zz}{𝐊}_k^T\right)$.

Outputs: ${\widehat{𝐱}}_{k|k}$, ${𝐏}_{k|k}$.

Table 3. TABLE III: ARMSEs and implementation times of the proposed filter and existing filters.

Filters	${\mathrm{ARMSE}}_{\mathrm{pos}}$ (km)	${\mathrm{ARMSE}}_{\mathrm{vel}}$ (km/min)	Time (ms)
SIF	22.85	0.83	45.46
ARUKF	$2.57\times {10}^{25}$	$2.41\times {10}^{24}$	$0.41$
MCCKF	26.86	0.88	0.13
RSTEF	$3.51\times {10}^6$	$6.51\times {10}^4$	$0.09$
3rd RSTUF	88.86	2.43	$0.53$
3rd RSTCF	29.90	1.12	$0.50$
5th RSTUF	16.29	0.65	$1.0$
RSTMCF	53.65	1.65	193.7
GSCKF	19.23	0.58	54.6
PF	7.89	0.30	773.0
RSTSCF	12.08	0.56	48.2

Equations

$${\mathbf{x}}_k={\mathbf{f}}_{k-1}({\mathbf{x}}_{k-1})+{\mathbf{w}}_{k-1}\text{(process equation)}$$

$${\mathbf{z}}_k={\mathbf{h}}_k({\mathbf{x}}_k)+{\mathbf{v}}_k\text{(measurement equation)},$$

$$p({\mathbf{w}}_k)=\mathrm{St}({\mathbf{w}}_k;0,{\mathbf{Q}}_k,{\nu }_1)$$

$$p({\mathbf{v}}_k)=\mathrm{St}({\mathbf{v}}_k;0,{\mathbf{R}}_k,{\nu }_2),$$

$${\hat{\mathbf{x}}}_{k|k-1}={\int }_{{\mathbb{R}}^n}{\mathbf{f}}_{k-1}({\mathbf{x}}_{k-1})\mathrm{St}({\mathbf{x}}_{k-1};{\hat{\mathbf{x}}}_{k-1|k-1},{\mathbf{P}}_{k-1|k-1},{\nu }_3)$$

$$d{\mathbf{x}}_{k-1}$$

$${\mathbf{P}}_{k|k-1}=\frac{{\nu }_3-2}{{\nu }_3}{\int }_{{\mathbb{R}}^n}{\mathbf{f}}_{k-1}({\mathbf{x}}_{k-1}){\mathbf{f}}_{k-1}^T({\mathbf{x}}_{k-1})\mathrm{St}({\mathbf{x}}_{k-1};$$

$${\hat{\mathbf{x}}}_{k-1|k-1},{\mathbf{P}}_{k-1|k-1},{\nu }_3)d{\mathbf{x}}_{k-1}-\frac{{\nu }_3-2}{{\nu }_3}{\hat{\mathbf{x}}}_{k|k-1}{\hat{\mathbf{x}}}_{k|k-1}^T$$

$$+\frac{{\nu }_1({\nu }_3-2)}{({\nu }_1-2){\nu }_3}{\mathbf{Q}}_{k-1},$$

$${\Delta }_k=\sqrt{({\mathbf{z}}_k-{\hat{\mathbf{z}}}_{k|k-1}{)}^T({\mathbf{P}}_{k|k-1}^{zz}{)}^{-1}({\mathbf{z}}_k-{\hat{\mathbf{z}}}_{k|k-1})}$$

$${\mathbf{K}}_k={\mathbf{P}}_{k|k-1}^{xz}({\mathbf{P}}_{k|k-1}^{zz}{)}^{-1}$$

$${\hat{\mathbf{x}}}_{k|k}={\hat{\mathbf{x}}}_{k|k-1}+{\mathbf{K}}_k({\mathbf{z}}_k-{\hat{\mathbf{z}}}_{k|k-1})$$

$${\mathbf{P}}_{k|k}=\frac{({\nu }_3-2)({\nu }_3+{\Delta }_k^2)}{{\nu }_3({\nu }_3+m-2)}({\mathbf{P}}_{k|k-1}-{\mathbf{K}}_k{\mathbf{P}}_{k|k-1}^{zz}{\mathbf{K}}_k^T),$$

$${\hat{\mathbf{z}}}_{k|k-1}$$

$${\mathbf{P}}_{k|k-1}^{zz}=\frac{{\nu }_3-2}{{\nu }_3}{\int }_{{\mathbb{R}}^n}{\mathbf{h}}_k({\mathbf{x}}_k){\mathbf{h}}_k^T({\mathbf{x}}_k)\mathrm{St}({\mathbf{x}}_k;{\hat{\mathbf{x}}}_{k|k-1},{\mathbf{P}}_{k|k-1},$$

$${\nu }_3)d{\mathbf{x}}_k-\frac{{\nu }_3-2}{{\nu }_3}{\hat{\mathbf{z}}}_{k|k-1}{\hat{\mathbf{z}}}_{k|k-1}^T+\frac{{\nu }_2({\nu }_3-2)}{({\nu }_2-2){\nu }_3}{\mathbf{R}}_k$$

$${\mathbf{P}}_{k|k-1}^{xz}=\frac{{\nu }_3-2}{{\nu }_3}{\int }_{{\mathbb{R}}^n}{\mathbf{x}}_k{\mathbf{h}}_k^T({\mathbf{x}}_k)\mathrm{St}({\mathbf{x}}_k;{\hat{\mathbf{x}}}_{k|k-1},{\mathbf{P}}_{k|k-1},{\nu }_3)$$

$$d{\mathbf{x}}_k-\frac{{\nu }_3-2}{{\nu }_3}{\hat{\mathbf{x}}}_{k|k-1}{\hat{\mathbf{z}}}_{k|k-1}^T.$$

$$I[\mathbf{g}]={\int }_{{\mathbb{R}}^n}\mathbf{g}(\mathbf{x})\mathrm{St}(\mathbf{x};\mu ,\mathbf{\Sigma },\nu )d\mathbf{x},$$

$$\mathrm{St}(\mathbf{x};\mu ,\mathbf{\Sigma },\nu )=\frac{\Gamma (\frac{\nu +n}{2})}{\Gamma (\frac{\nu }{2})}\frac{1}{\sqrt{|\nu \pi \mathbf{\Sigma }|}}\times$$

$${\left[1+\frac{1}{\nu }(\mathbf{x}-\mu {)}^{\mathrm{T}}{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mu )\right]}^{-\frac{\nu +n}{2}},$$

$$\mathbf{x}=\mu +\sqrt{\nu \mathbf{\Sigma }}\mathbf{y},$$

$$I[\mathbf{g}]={\int }_{{\mathbb{R}}^n}\mathbf{l}(\mathbf{y})(1+{\mathbf{y}}^{\mathrm{T}}\mathbf{y}{)}^{-\frac{\nu +n}{2}}d\mathbf{y},$$

$$\mathbf{l}(\mathbf{y})=\frac{\Gamma (\frac{\nu +n}{2})}{\Gamma (\frac{\nu }{2}){\pi }^{\frac{n}{2}}}\mathbf{g}(\mu +\sqrt{\nu \mathbf{\Sigma }}\mathbf{y}).$$

$$I[\mathbf{g}]$$

$$I[\mathbf{g}]={\int }_0^{+\infty }\mathbf{S}(r)r^{n-1}(1+r^2{)}^{-\frac{\nu +n}{2}}dr,$$

$$\mathbf{S}(r)={\int }_{U_n}\mathbf{l}(r\mathbf{s})d\sigma (s).$$

$$\int \mathbf{g}(\mathbf{x})p(\mathbf{x})d\mathbf{x}=E\left[\sum _{l=1}^N{w}_l\mathbf{g}({\mathbf{x}}_l)\right],$$

$${\mathbf{S}}_u^3(r)\approx \frac{A_n}{2n}\sum _{i=1}^n\left[\mathbf{l}(-r\mathbf{Q}{\mathbf{e}}_i)+\mathbf{l}(r\mathbf{Q}{\mathbf{e}}_i)\right],$$

$${\int }_0^{+\infty }\mathbf{S}(r)r^{n-1}(1+r^2{)}^{-\frac{\nu +n}{2}}dr\approx {\omega }_{r,1}\mathbf{S}(r_1)+{\omega }_{r,2}\mathbf{S}(r_2),$$

$${\omega }_{r,1}{r}_1^0+{\omega }_{r,2}{r}_2^0={\int }_0^{+\infty }{r}^0{r}^{n-1}(1+r^2{)}^{-\frac{\nu +n}{2}}dr$$

$${\omega }_{r,1}{r}_1^2+{\omega }_{r,2}{r}_2^2={\int }_0^{+\infty }{r}^2{r}^{n-1}(1+r^2{)}^{-\frac{\nu +n}{2}}dr$$

$${\int }_0^{+\infty }\mathbf{S}(r)r^{n-1}(1+r^2{)}^{-\frac{\nu +n}{2}}dr=E\left[{\omega }_{r,1}\mathbf{S}(r_1)+{\omega }_{r,2}\mathbf{S}(r_2)\right].$$

Full text

I Introduction

Nonlinear filtering has been playing an important role in many applications, such as target tracking, detection, signal processing, communication and navigation. Under the Bayesian estimation framework, the nonlinear filtering problem is addressed by calculating the posterior probability density function (PDF) recursively based on the nonlinear state-space model. Unfortunately, there is not a closed form solution for posterior PDF for nonlinear state-space model since a closed PDF for nonlinear mapping doesn’t exist [1]. As a result, there is not an optimal solution for nonlinear filtering problem, and an approximate approach is necessary to obtain a suboptimal solution. In general, the posterior PDF is approximated as Gaussian by assuming the jointly predicted PDF of the state and measurement vectors is Gaussian, and the resultant Gaussian approximate (GA) filter can provide tradeoffs between estimation accuracy and computational complexity [2, 3]. Up to present, many variants of the GA filter have been proposed using different Gaussian weighted integral rules [3]–[9]. However, in some engineering applications, such as tracking an agile target that is observed in clutter, the heavy-tailed process noise may be induced by severe manoeuvering and the heavy-tailed measurement noise may be induced by measurement outliers from unreliable sensors [10, 11, 12]. The performance of the GA filters may degrade for such engineering applications with heavy-tailed noises since they all model the process and measurement noises as Gaussian distributions so that they are sensitive to heavy-tailed non-Gaussian noises [11].

Particle filter (PF) is a common method to address non-Gaussian noises, in which the posterior PDF is approximated as a set of random samples with associated weights based on sequential Monte Carlo sampling technique [13]. The PF can model the process and measurement noises as arbitrary distributions, such as the Student’s t distributions for heavy-tailed non-Gaussian noises [14, 15]. However, the PF suffers from substantial computational complexities in high-dimensional problems because the number of particles increases exponentially with the dimensionality of the state [16]. Gaussian sum filter (GSF) is an alternative method to handle heavy-tailed non-Gaussian noises, where the heavy-tailed process and measurement noises are modelled as a finite sum of Gaussian distributions, and the posterior distribution is then approximated as a weighted sum of Gaussian distributions by running a bank of GA filters [17]–[19]. However, for the GSF, it is very difficult to model the heavy-tailed process and measurement noises accurately using finite Gaussian distributions since the heavy-tailed non-Gaussian noises are induced by the unknown manoeuvering or outliers, which may degrade the estimation performance of the GSF.

To solve the filtering problem of nonlinear state-space model with heavy-tailed non-Gaussian noises, the Huber-based nonlinear Kalman filter (HNKF) has been proposed by minimising a Huber cost function that is a combined $l_{1}$ and $l_{2}$ norm [20]. A larger number of variants of the HNKF have been derived based on a linearized or statistical linearized method, such as the Huber-based extended Kalman filter [21], the Huber-based divided difference filter [22], the Huber-based unscented Kalman filter [23], the nonlinear regression Huber Kalman filter [24] and the adaptively robust unscented Kalman filter (ARUKF) [25]. However, the influence function of the HNKFs don’t redescend, which may deteriorate the estimation performance of the HNKFs [12]. The maximum correntropy criterion based Kalman filter (MCCKF) has been proposed by maximising the correntropy of the predicted error and residual [26]–[29]. However, there is a lack of theoretical basis to develop the estimation error covariance matrix of the MCCKF, which may degrade the estimation accuracy [12].

A reasonable approach to improve the estimation performance is utilizing a Student’s t distribution to model the heavy-tailed non-Gaussian noise. The Student’s t distribution is a generalized Gaussian distribution but has heavier tails than the Gaussian distribution, which makes it more suitable for modelling the heavy-tailed non-Gaussian noise. A general framework of the robust Student’s t based nonlinear filter (RSTNF) has been proposed, in which the jointly predicted PDF of the state and measurement vectors is assumed to be Student’s t, and the posterior PDF is then approximated as Student’s t [30]. The heart of the RSTNF is how to calculate the Student’s t weighted integral, and the estimation accuracy of the associated RSTNF is determined by the employed numerical integral technique. Many variants of the RSTNF have been derived based on different numerical integral methods, such as the robust Student’s t based extended filter (RSTEF) using the first-order linearization [10], the robust Student’s t based unscented filter (RSTUF) using the unscented transform (UT) [30, 31], and the robust Student’s t based cubature filter (RSTCF) using the third-degree Student’s t spherical radial cubature rule (STSRCR) [32]. However, the existing Student’s t integral rules can only capture the third-degree or fifth-degree information of the Taylor series expansion for nonlinear approximation, which may result in limited estimation accuracy. Although the Monte Carlo approach can be used to calculate the Student’s t weighted integral, it has low accuracy and slow convergence when the integrand is not approximately constant and the number of random samples is finite [33]. Therefore, there is a great demand to develop more accurate numerical integral approach for the Student’s t weighted integral to further improve the estimation accuracy of the existing RSTNFs.

In this paper, the Student’s t weighted integral is decomposed into the spherical integral and the radial integral based on the spherical-radial transformation. A new stochastic STSRCR (SSTSRCR) is derived based on the third-degree unbiased spherical rule (USR) and the proposed third-degree unbiased radial rule (URR), from which a new robust Student’s t based stochastic cubature filter (RSTSCF) is obtained. The existing stochastic integration rule (SIR) [8] is a special case of the proposed SSTSRCR when the degrees of freedom (dof) parameter tends to infinity. The proposed SSTSRCR can achieve better approximation to the Student’s t weighted integral as compared with existing Student’s t integral rules. As a result, the proposed RSTSCF has higher estimation accuracy than the existing RSTNFs. The proposed filter and existing filters are applied to a manoeuvring bearings-only tracking example, where an agile target is tracked and the bearing is observed in clutter. Simulation results show that the proposed RSTSCF can achieve higher estimation accuracy than the existing GA filter, GSF, HNKF, MCCKF and RSTNFs, and is computationally much more efficient than the existing PF.

The remainder of this paper is organized as follows. In Section II, a general frame of the RSTNF is reviewed. In Section III, a new SSTSRCR is derived based on the proposed third-degree URR, from which a new RSTSCF is obtained. Also, the relationship between the proposed SSTSRCR and the existing SIR is revealed in Section III. In Section IV, the proposed filter is applied to a manoeuvring bearings-only tracking example and simulation results are given. Concluding remarks are drawn in Section V.

II Problem Statement

Consider the following discrete-time nonlinear stochastic system as represented by the state-space model [30]

[TABLE]

[TABLE]

where ${k}$ is the discrete time index, $\mathbf{x}_{k}\in\mathbb{R}^{n}$ is the state vector, $\mathbf{z}_{k}\in\mathbb{R}^{m}$ is the measurement vector, and $\mathbf{f}_{k-1}(\cdot)$ and $\mathbf{h}_{k}(\cdot)$ are known process and measurement functions respectively. $\mathbf{w}_{k}\in\mathbb{R}^{n}$ and $\mathbf{v}_{k}\in\mathbb{R}^{m}$ are heavy-tailed process and measurement noise vectors respectively, which are induced by process and measurement outliers, and their distributions are modelled as Student’s t distributions, i.e.,

[TABLE]

[TABLE]

where $\mathrm{St}(\cdot;\mathbf{\mu},\mathbf{\Sigma},v)$ denotes the Student’s t PDF with mean vector $\mathbf{\mu}$, scale matrix $\mathbf{\Sigma}$, and dof parameter $v$, $\mathbf{Q}_{k}$ and $\nu_{1}$ are the scale matrix and dof parameter of process noise respectively, and $\mathbf{R}_{k}$ and $\nu_{2}$ are the scale matrix and dof parameter of measurement noise respectively. The initial state vector $\mathbf{x}_{0}$ is also assumed to have a Student’s t distribution with mean vector $\mathbf{\hat{x}}_{0|0}$, scale matrix $\mathbf{P}_{0|0}$, and dof parameter $\nu_{3}$, and $\mathbf{x}_{0}$, $\mathbf{w}_{k}$ and $\mathbf{v}_{k}$ are assumed to be mutually uncorrelated.

To achieve the filtering estimation, a general framework of RSTNF is derived for the nonlinear system formulated in equations (1)-(4), where the jointly predicted PDF of the state and measurement vectors is assumed as Student’s t, then the posterior PDF of the state vector can be approximated as Student’s t [30]. The time update and measurement update of the recursive RSTNF are given as follows:

Time update

[TABLE]

[TABLE]

where $(\cdot)^{T}$ denotes the transpose operation, $\mathbf{\hat{x}}_{k|k-1}$ and $\mathbf{P}_{k|k-1}$ are respectively the mean vector and scale matrix of the one-step predicted PDF $p(\mathbf{x}_{k}|\mathbf{Z}_{k-1})$, $\mathbf{Z}_{k-1}=\{\mathbf{z}_{j}\}_{j=1}^{k-1}$ is the set of $k-1$ measurement vectors, and $\nu_{3}$ denotes the dof parameter of the filtering PDF.

Measurement update

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $(\cdot)^{\mathrm{-1}}$ denotes the inverse operation, $\mathbf{\hat{x}}_{k|k}$ and $\mathbf{P}_{k|k}$ are respectively the mean vector and scale matrix of the filtering PDF $p(\mathbf{x}_{k}|\mathbf{Z}_{k})$, $\mathbf{\hat{z}}_{k|k-1}$ and $\mathbf{P}_{k|k-1}^{zz}$ are respectively the mean vector and scale matrix of the likelihood PDF $p(\mathbf{z}_{k}|$ $\mathbf{Z}_{k-1})$, and $\mathbf{P}_{k|k-1}^{xz}$ is the cross scale matrix of state and measurement vectors, which are given by

[TABLE]

[TABLE]

[TABLE]

The recursive RSTNF is composed of the analytical computations in equations (7)-(10) and the Student’s t weighted integrals in equations (5)-(6) and (11)-(13). The key problem in the design of the RSTNF is calculating the nonlinear Student’s t weighted integrals formulated in equations (5)-(6) and (11)-(13), whose integrands are all of the form nonlinear function$\times$Student’s t PDF. Therefore, the numerical integral technique is required to implement the RSTNF, which determines the estimation accuracy of associated RSTNF. Next, to further improve the estimation accuracy of existing RSTNFs, a new SSTSRCR will be proposed, based on which a new RSTSCF can be obtained.

III Main Results

III-A Spherical-radial transformation

The Student’s t weighted integrals involved in the RSTNF can be written as the general form as follows

[TABLE]

where the Student’s t PDF is given by

[TABLE]

where $\Gamma(\cdot)$ and $|\cdot|$ denote the Gamma function and determinant operation respectively. To derive the SSTSRCR, the Student’s t weighted integral in equation (14) requires to be transformed into a spherical-radial integral form.

A change of variable is utilized as follows

[TABLE]

where $\sqrt{\mathbf{\Sigma}}$ is the square-root of scale matrix $\mathbf{\Sigma}$ satisfying $\mathbf{\Sigma}=\sqrt{\mathbf{\Sigma}}\sqrt{\mathbf{\Sigma}}^{T}$.

Substituting equation (16) in equations (14)-(15) and using the identity $|\sqrt{\nu\mathbf{\Sigma}}|=\sqrt{|\nu\mathbf{\Sigma}|}$ yields

[TABLE]

where $\mathbf{l}(\mathbf{y})$ is given by

[TABLE]

Define $\mathbf{y}=r\mathbf{s}$ with $\mathbf{s}^{T}\mathbf{s}=1$, then equation (17) can be rewritten as [34]

[TABLE]

where $\mathbf{s}=[\mathbf{s}_{1},\mathbf{s}_{2},\cdots,\mathbf{s}_{n}]^{T}$, $U_{n}=\{\mathbf{s}\in\mathbb{R}^{n}:s_{1}^{2}+s_{2}^{2}+\cdots+s_{n}^{2}=1\}$, and $\sigma(\mathbf{s})$ is the spherical surface measure or an area element on $U_{n}$.

According to equation (19), the Student’s t weighted integral in equation (14) can be decomposed into the radial integral

[TABLE]

and the spherical integral

[TABLE]

Next, a new third-degree SSTSRCR will be derived, in which the spherical and the radial integrals are respectively calculated by the third-degree USR (Section III. B below) and the third-degree URR (Section III. C below). Before deriving the third-degree SSTSRCR, the unbiased integral rule is firstly defined as follows.

Definition 1**.**

The integral rule $\int\mathbf{g}(\mathbf{x})p(\mathbf{x})d\mathbf{x}\approx\sum\limits_{l=1}^{N}w_{l}\mathbf{g}(\mathbf{x}_{l})$ is unbiased if and only if [33]

[TABLE]

where $\mathbf{x}_{l}$ and $w_{l}$ are respectively cubature points and corresponding weights, and $E[\cdot]$ denotes the expectation operation.

III-B Unbiased spherical rule

The Stewart’s method is employed to construct the third-degree USR. If $\mathbf{Q}$ is a random orthogonal matrix drawn with a Haar distribution from the set of all matrices in the orthogonal group, the third-degree USR can be constructed as [33, 35]

[TABLE]

where $A_{n}=\frac{2\pi^{\frac{n}{2}}}{\Gamma(\frac{n}{2})}$ is the surface area of the unit sphere, and $\mathbf{e}_{i}$ denotes the $i$-th column of an $n\times{n}$ unit matrix. To produce a random orthogonal matrix $\mathbf{Q}$, a $n\times n$ matrix $\mathbf{U}$ of standard norm variables is first generated, then the required random orthogonal matrix $\mathbf{Q}$ is obtained based on the $QR$ factorization, i.e., $\mathbf{U}=\mathbf{Q}\mathbf{R}$ [35].

Next, a new third-degree URR will be proposed for the radial integral in equation (20).

III-C Unbiased radial rule

Generally, the monomials $S(r)=1$, $S(r)=r$, $S(r)=r^{2}$, and $S(r)=r^{3}$ need to be matched to derive the third-degree URR. However, only monomials $S(r)=1$ and $S(r)=r^{2}$ need to be matched for the third-degree URR since the USR and the resultant STSRCR are fully symmetry. Thus, two points $\{r_{1},\omega_{r,1}\}$ and $\{r_{2},\omega_{r,2}\}$ are sufficient to design the third-degree URR, where one point is used to match monomials $S(r)=1$ and $S(r)=r^{2}$ and the other point is employed to retain unbiasedness. That is to say, the third-degree URR can be written as

[TABLE]

where $\{r_{1},\omega_{r,1}\}$ and $\{r_{2},\omega_{r,2}\}$ satisfy the following equations

[TABLE]

[TABLE]

[TABLE]

Since there are three equations and four variables in equations (25)-(27), there is one free variable. In order to get the STSRCR with the minimum number of points, $r_{1}$ is chosen as the free variable and set to zero.

Theorem 1**.**

If $r_{1}=0$ and the PDF of random variable $r_{2}$ is

[TABLE]

where $\mathrm{B}(\cdot,\cdot)$ denotes the beta function, then the third-degree URR is given by

[TABLE]

Proof.

Firstly, a general integral $\int_{0}^{+\infty}r^{l}r^{n-1}(1+r^{2})^{-\frac{\nu+n}{2}}dr$ is calculated to obtain the right-hand parts in equations (25)-(26). Making a change of variable via $t=r^{2}$ results in

[TABLE]

where $\mathrm{B}(\cdot,\cdot)$ denotes the beta function.

Substituting equation (30) in equations (25)-(26), we have

[TABLE]

[TABLE]

Utilizing the identities $\Gamma(a+1)=a\Gamma(a)$ and $\mathrm{B}(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}$ in equation (32) yields

[TABLE]

Employing $r_{1}=0$ in equation (33) yields

[TABLE]

Substituting equation (34) in equation (31) results in

[TABLE]

Utilizing $r_{1}=0$ and equations (34)-(35), the expectation of the third-degree radial rule with respect to $p(r_{2})$ is written as

[TABLE]

Using equations (28) and (30), we have

[TABLE]

[TABLE]

Substituting equations (37)-(38) in equation (36) yields

[TABLE]

With $r_{1}=0$, equations (34)-(35) and (39), the third-degree URR can be formulated as equation (29). ∎

It is very difficult to directly generate random samples from $p(r_{2})$ since $p(r_{2})$ is not a special PDF. To solve this problem, Theorem 2 is presented as follows.

Theorem 2**.**

If random variable $\tau=\frac{r_{2}^{2}}{1+r_{2}^{2}}$, then random variable $\tau$ obeys the Beta distribution, i.e.,

[TABLE]

where $\mathrm{Beta}(\cdot;\alpha,\beta)$ denotes the beta PDF with parameters $\alpha$ and $\beta$.

Proof.

Since $\tau=\frac{r_{2}^{2}}{1+r_{2}^{2}}$ and $r_{2}\in[0,+\infty)$, random variable $\tau\in[0,1)$. According to $\tau=\frac{r_{2}^{2}}{1+r_{2}^{2}}$, $r_{2}$ is formulated as

[TABLE]

Employing the transformation theorem and equation (41), the PDF of random variable $\tau$ is given by

[TABLE]

where $p_{r_{2}}(\cdot)$ denotes the PDF of $r_{2}$ and $c^{\prime}(\tau)$ denotes the derivative of $c(\tau)$ with respect to $\tau$ given by

[TABLE]

Substituting equations (28), (41) and (43) in equation (42) obtains

[TABLE]

which proves the theorem. ∎

III-D Stochastic STSRCR

A theorem is first presented to derive the unbiased STSRCR.

Theorem 3**.**

If the spherical and radial rules are unbiased, then the resultant STSRCR is also unbiased.

Proof.

If the spherical and radial rules are given by

[TABLE]

then the STSRCR can be formulated as

[TABLE]

where $\mathbf{s}_{i}$ and $w_{s,i}$ are respectively cubature points and weights of the spherical rule, and $r_{j}$ and $w_{r,j}$ are respectively quadrature points and weights of the radial rule.

Since the spherical and radial rules are unbiased, we obtain

[TABLE]

Using equation (47) yields

[TABLE]

Since the set $\{\mathbf{s}_{i},\,w_{s,i}\}_{i=1}^{N_{s}}$ is independent of the set $\{\mathbf{r}_{j},\,w_{r,j}\}_{j=1}^{N_{r}}$, we have

[TABLE]

which proves the theorem. ∎

Using Theorems 1-3 obtains

[TABLE]

By employing the Monte Carlo approach, the right-hand parts of equation (50) can be approximated as

[TABLE]

where $N$ denotes the number of random samples, and $\mathbf{Q}_{l}$ is a random orthogonal matrix, and $r_{2,l}$ is drawn randomly from $p(r_{2})$. The form $I_{s}^{3}[\mathbf{g}]$ denotes the proposed third-degree SSTSRCR, and the implementation pseudocode of the proposed SSTSRCR is shown in Table I.

According to the Monte Carlo approach, $I_{s}^{3}[\mathbf{g}]$ converges to $I[\mathbf{g}]$ when $N$ tends to infinity, i.e.,

[TABLE]

Thus, the proposed SSTSRCR provides asymptotically exact integral evaluations when $N$ tends to infinity. A new RSTSCF can be obtained by employing the proposed SSTSRCR to calculate the Student’s t weighted integrals involved in the RSTNF, and the implementation pseudocode for one time step of the proposed RSTSCF is shown in Table II, where $\mathrm{SSTSRCR}(\cdot)$ denotes the proposed SSTSRCR algorithm. The proposed SSTSRCR can achieve better approximation to the Student’s t weighted integral as compared with existing Student’s t integral rules. As a result, the proposed RSTSCF has higher estimation accuracy than the existing RSTNFs.

Remark 1**.**

The Monte Carlo approach can be also used to calculate the Student’s t weighted integral, and it provides asymptotically exact integral evaluations when the number of random samples tends to infinity. However, it has low accuracy and slow convergence when the integrand is not approximately constant and the number of random samples is finite [33]. Fortunately, the proposed SSTSRCR is at least exact up to third-degree polynomials for any number of random samples, and it can capture more and more higher-degree moment information as the number of random samples increases.

III-E Relationship between the proposed SSTSRCR and the existing SIR [8]

Theorem 4**.**

The proposed SSTSRCR will degrade to the existing SIR when the dof parameter $\nu\rightarrow+\infty$, i.e.

[TABLE]

where the right-hand side of the equation (53) is the SIR for the Gaussian weighted integral, and $\rho_{l}$ is drawn randomly from $p(\rho_{l})$ that is given by

[TABLE]

Proof.

Make a change of variable as follows

[TABLE]

Substituting equation (55) in equation (51) results in

[TABLE]

Taking the limit of equation (56) when the dof parameter $\nu\rightarrow+\infty$, we can obtain equation (53).

Using the transformation theorem and equation (55), the PDF of random variable $\rho_{l}$ is given by

[TABLE]

where $p_{r_{2}}(\cdot)$ denotes the PDF of $r_{2}$ , and $c^{\prime}(\rho_{l})$ denotes the derivative of $c(\rho_{l})$ with respect to $\rho_{l}$ given by

[TABLE]

Substituting equations (28), (55) and (58) in equation (57), we obtain

[TABLE]

Utilizing the identity $\mathrm{B}(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}$, the first limit in equation (59) can be formulated as

[TABLE]

Using the property of Gamma function $\lim\limits_{t\rightarrow+\infty}\frac{\Gamma(t+\alpha)}{\Gamma(t)t^{\alpha}}=1$ in equation (60) gives

[TABLE]

The second limit in equation (59) can be reformulated as

[TABLE]

where the functions $s(\nu)$ and $d(\nu)$ are given by

[TABLE]

[TABLE]

Using the identity $\lim\limits_{t\rightarrow+\infty}\left(1+\frac{1}{t}\right)^{t}=e$ and equations (63)-(64), equation (62) can be rewritten as

[TABLE]

Substituting equations (61) and (65) in (59), we can obtain (54), which proves the theorem. ∎

Considering that the Student’s t PDF turns into the Gaussian PDF as the dof parameter $\nu\rightarrow+\infty$, we obtain

[TABLE]

According to the Theorem 4 and equation (66), we can conclude that the proposed SSTSRCR with $\nu\rightarrow+\infty$ can be utilized to calculate the Gaussian weighted integral. Thus, the proposed SSTSRCR is a generalized SIR, which can calculate not only the Gaussian weighted integral but also the Student’s t weighted integral.

IV Simulation study

In this simulation, the superior performance of the proposed RSTSCF as compared with existing filters is shown in the problem of manoeuvring bearing-only tracking observed in clutter. The target moves according to the continuous white noise acceleration motion model [8]

[TABLE]

where $\mathbf{x}_{k}=[x_{k}\;y_{k}\;\dot{x}_{k}\;\dot{y}_{k}]$, and $x_{k}$, $y_{k}$, $\dot{x}_{k}$ and $\dot{y}_{k}$ denote the cartesian coordinates and corresponding velocities respectively; $\mathbf{F}$ and $\mathbf{G}$ denote respectively the state transition matrix and noise matrix given by

[TABLE]

where $\Delta{t}=1\mathrm{min}$ is the sampling interval, and $\mathbf{I}_{2}$ is the two dimensional identity matrix, and $\mathbf{0}_{2\times 1}$ is the two dimensional zero vector, and $\mathbf{\Gamma}=[0.5\Delta{t}^{2}\;\;\Delta{t}]^{T}$.

The target is observed by an angle sensor installed in a manoeuvring platform. If the platform is located at $(x_{k}^{p},y_{k}^{p})$ at time $k$, then the measurement model is given by

[TABLE]

where $z_{k}$ is the angle between the target and the platform at time $k$. Outlier corrupted process and measurement noises are generated according to [10, 12, 30]

[TABLE]

[TABLE]

where $\mathrm{w.p.}$ denotes “with probability”, the nominal process noise covariance matrix is $\mathbf{\Sigma}_{w}=10^{-6}\mathbf{I}_{2}\mathrm{km}^{2}/\mathrm{min}^{2}$, and the nominal measurement noise variance is $\Sigma_{v}=(0.02\mathrm{rad})^{2}$. Process and measurement noises, which are generated according to equations (70)-(71), have heavier tails.

In our simulation scenario, the initial positions of the target and the platform are respectively $(3\mathrm{km},3\mathrm{km})$ and $(0\mathrm{km},0\mathrm{km})$. The target moves at a constant speed of 180 knots (1 knot is $1.852\mathrm{km}/\mathrm{h}$) with a course of $-135.4^{\circ}$. The platform moves at a constant speed of 50 knots with a initial course of $-80^{\circ}$, and the course reaches $146^{\circ}$ at time $k=15\mathrm{min}$ by executing a manoeuvre [8]. The initial estimation error covariance matrix is $\bm{P}_{0|0}=\mathrm{diag}[16\mathrm{km}^{2}\,16\mathrm{km}^{2},4\mathrm{km}^{2}/\mathrm{min}^{2},$ $4\mathrm{km}^{2}/\mathrm{min}^{2}]$, and the initial state estimate $\bm{\hat{x}}_{0|0}$ is chosen randomly from $\mathrm{N}(\bm{x}_{0},\bm{P}_{0|0})$, where $\bm{x}_{0}$ denotes the initial true state.

In this simulation, the stochastic integration filter (SIF) [8], the ARUKF with free parameter $\kappa=0$ [25], the MCCKF with kernel size $\sigma=5$ [29], the RSTEF [10], the 3rd-degree RSTUF with free parameter $\kappa=3-n$ [30], the 3rd-degree RSTCF [32], the fifth-degree RSTUF [31], the robust Student’s t based Monte Carlo filter (RSTMCF), the Gaussian sum-cubature Kalman filter (GSCKF) [18], the PF [13], and the proposed RSTSCF are tested. In the RSTMCF, the Student’s t weighted integral is calculated using the conventional Monte Carlo approach with $10000$ random samples. In the GSCKF, the process and measurement noises are modelled as $p(\mathbf{w}_{k})=\sum\limits_{i=1}^{5}\alpha_{i}\mathrm{N}\left(\mathbf{w}_{k};\mathbf{0},\lambda_{i}\mathbf{\Sigma}_{w}\right)$ and $p(v_{k})=\sum\limits_{i=1}^{5}\alpha_{i}\mathrm{N}\left(v_{k};0,\lambda_{i}\Sigma_{v}\right)$, where the weights $\alpha_{1}=0.8$ and $\alpha_{2}=\alpha_{3}=\alpha_{4}=\alpha_{5}=0.05$, and the scale parameters $\lambda_{1}=1$, $\lambda_{2}=50$, $\lambda_{3}=100$, $\lambda_{4}=500$ and $\lambda_{5}=1000$. Moreover, to prevent the computational complexity of the GSCKF increasing exponentially as the time, the posterior distribution is approximated as a weighted sum of five Gaussian terms with the highest weights. In the PF, the process and measurement noises are modelled as Student’s t distributions, and the number of particle is chosen as $10000$. In the existing RSTEF, 3rd-degree RSTUF, 3rd-degree RSTCF, fifth-degree RSTUF, RSTMCF, PF and the proposed RSTSCF, the dof parameters are all chosen as $\nu_{1}=\nu_{2}=\nu_{3}=5$ and the scale matrices are all set as $\mathbf{Q}_{k}=\mathbf{\Sigma}_{w}$ and ${R}_{k}={\Sigma}_{v}$. In the SIF and the proposed RSTSCF, the number of random samples is selected as $N=100$. The proposed filter and existing filters are coded with MATLAB and the simulations are run on a computer with Intel Core i7-3770 CPU at 3.40 GHz.

To compare the performances of the proposed filter and existing filters, the RMSEs and the averaged RMSEs (ARMSEs) of the position and velocity are chosen as performance metric. The RMSE and ARMSE in position are respectively defined as

[TABLE]

[TABLE]

where $M=1000$ denotes the number of Monte Carlo runs, and $T=100\mathrm{min}$ denotes the simulation time, and $(x_{k}^{s},y_{k}^{s})$ and $(\hat{x}_{k}^{s},\hat{y}_{k}^{s})$ respectively denote the true and estimated positions at the $s$-th Monte Carlo run. Similar to the RMSE and ARMSE in position, we can also formulate the RMSE and ARMSE in velocity.

The RMSEs and ARMSEs of position and velocity from the proposed filter and existing filters are respectively shown in Fig. 1–Fig. 2 and Table III. The implementation times of the proposed filter and existing filters in single step run are given in Table III. Note that the existing ARUKF and RSTEF diverge in the simulation, as shown in Fig. 1–Fig. 2 and Table III.

It is seen from Fig. 1–Fig. 2 and Table III that the RMSEs and ARMSEs of the proposed RSTSCF are smaller than the existing SIF, ARUKF, MCCKF, RSTEF, 3rd-degree RSTUF, 3rd-degree RSTCF, 5th-degree RSTUF, RSTMCF, GSCKF but larger than the existing PF. Furthermore, it can be also seen from Table III that the implementation time of the proposed RSTSCF are greater than the existing SIF, ARUKF, MCCKF, RSTEF, 3rd-degree RSTUF, 3rd-degree RSTCF, 5th-degree RSTUF but significantly smaller than the existing RSTMCF, GSCKF and PF. Therefore, the proposed RSTSCF has better estimation accuracy than the existing SIF, ARUKF, MCCKF, RSTNFs and GSCKF, and is computationally much more efficient than the existing PF.

V Conclusion

In this paper, a new SSTSRCR was derived based on the third-degree USR and the proposed third-degree URR, from which a new RSTSCF was obtained. The existing SIR is a special case of the proposed SSTSRCR when the dof parameter tends to infinity. Simulation results for a manoeuvring bearings-only tracking example illustrated that the proposed RSTSCF can achieve higher estimation accuracy than the existing GA filter, GSF, HNKF, MCCKF and RSTNFs, and is computationally much more efficient than the existing PF.