A Robust SMC-PHD Filter for Multi-Target Tracking with Unknown Heavy-Tailed Measurement Noise

Abstract

In multi-target tracking, the sequential Monte Carlo probability hypothesis density (SMC-PHD) filter is a practical algorithm. Influenced by outliers under unknown heavy-tailed measurement noise, the SMC-PHD filter suffers severe performance degradation. In this paper, a robust SMC-PHD (RSMC-PHD) filter is proposed. In the proposed filter, Student-t distribution is introduced to describe the unknown heavy-tailed measurement noise where the degrees of freedom (DOF) and the scale matrix of the Student-t distribution are respectively modeled as a Gamma distribution and an inverse Wishart distribution. Furthermore, the variational Bayesian (VB) technique is employed to infer the unknown DOF and scale matrix parameters while the recursion estimation framework of the RSMC-PHD filter is derived. In addition, considering that the introduced Student- t distribution might lead to an overestimation of the target number, a strategy is applied to modify the updated weight of each particle. Simulation results demonstrate that the proposed filter is effective with unknown heavy-tailed measurement noise.

1. Introduction

Multi-target tracking (MTT) involves the estimate of target number and states from noisy measurements [1]. Influenced by clutter and detection uncertainty, the traditional MTT algorithms based on data association, such as joint probabilistic data association (JPDA) [2] and multiple hypotheses tracking (MHT) [3], are burdened with high computation complexity with the increase of target number and clutter number. Under the framework of finite set statistics (FISST), multi-target states and measurements are modeled as random finite sets (RFS) to avoid data association. Unfortunately, the MTT approaches based on RFS are still computationally intractable because of the multiple integrals. To obtain a computational tractable solution, the probability hypothesis density (PHD) filter [4] which propagates the first order moment approximation of multi-target posterior density is proposed and realizes the estimate of multi-target number and states. The two main implementations of the PHD filter are the Gaussian mixture PHD (GM-PHD) filter [5] and the sequential Monte Carle PHD (SMC-PHD) [6]. In recent years, the PHD filter has generated substantial interest in visual tracking [7], radar tracking [8], etc.

However, the above algorithms only achieve good performance with the known Gaussian measurement noise variance. When the measurement noise variance is unknown and inconsistent with the true variance, the performance of these algorithms decreases. In [9], Sarkka considers the unknown measurement noise variance problem under the linear Gaussian system and models noise variance as an inverse Gamma (IG) distribution. Then, the joint estimate of target state and noise variance is obtained by variational Bayesian (VB) technique. Since the inverse Gamma distribution is only suitable for the situation where the noise variance is diagonal, the inverse Wishart (IW) distribution is further used to model the full noise variance [10,11]. Following this, the VB technique was introduced into the PHD filter [12,13,14] for MTT.

In some applications, the measurement noise might be non-Gaussian. For example, due to the electromagnetic interference and the sensor’s unreliability, the measurement model is usually accompanied with outliers and the measurement noise can be treated as non-Gaussian noise with a heavy-tailed probability density function [15]. The outliers will result in the performance degradation of conventional tracking algorithms. In [15], the combination of a Laplacian distribution and a Gaussian distribution is used to describe the heavy-tailed measurement noise. In [16,17], the heavy-tailed measurement noise is characterized as a mixture of two Gaussian distributions. Due to the advantage of the heavy-tailed characteristic, Student-t distribution was widely used to deal with outliers in [18,19,20]. However, the relevant parameters of Student-t distribution need to be known in these applications. Similar to [9,10,11,12,13,14], the VB technique was employed in the situation where the Student-t distribution parameters are unknown in [21,22,23,24]. The joint estimate of target and unknown parameters is obtained in a linear system. However, the above approaches work worse in a nonlinear system. The nonlinear transformation algorithms, such as extended Kalman filter (EKF) and unscented Kalman filter (UKF), can be used to cope with the nonlinear system to some extent. Unfortunately, the EKF depends on the nature of nonlinearities and the UKF relies on the Gaussian characteristic of the system. They often give an unreliable estimate if the nonlinearities are severe or the system is non-Gaussian, which leads to target missing.

With heavy-tailed measurement noise, the robust PHD filter in a nonlinear system has been much less studied than that in a linear system. Compared with nonlinear transformations, the particle filter is more suitable for practical applications because it can deal with nonlinear and non-Gaussian dynamics. Motivated by [23], we extend the Student-t distribution into SMC-PHD filter and propose a robust SMC-PHD (RSMC-PHD) filter with unknown heavy-tailed measurement noise. In the proposed filter, the Student-t distribution where the degrees of freedom (DOF) and the scale matrix are both unknown is used to model the measurement noise. The Gamma distribution and IW distribution are respectively used to model the DOF and the scale matrix. Then, the augmented state of the target involves the target state, DOF, and the scale matrix. Further, the target state is estimated by particle approximations and the DOF and the scale matrix are updated by the VB technique within the SMC-PHD filter framework. In addition, considering the problem that the introduced Student-t distribution might lead to an overestimation of the target number, we apply a strategy to modify the updated weight of each particle. Simulation results show that the proposed RSMC-PHD filter can achieve higher estimate accuracy than the compared filter under unknown heavy-tailed measurement noise.

This paper is organized as follows. Section 2 gives a brief overview of the tracking model, measurement noise model, PHD filter, and SMC-PHD filter. Section 3 presents the detailed implementation of the proposed filter. Section 4 evaluates the performance of the proposed filter under two scenarios. Finally, Section 5 displays the conclusions.

2. Background

2.1. Tracking Model and Measurement Noise Model

Consider a state space system given by

$$x_k=f(x_{k-1})+w_k$$

(1)

$$z_k=h(x_k)+{\epsilon }_k$$

(2)

where $x_k$ is referred to as the target state and $z_k$ is referred to as the measurement at time $k$. $w_k$ is the process noise assumed to be Gaussian white noise with covariance matrix $Q_k$. $f(\cdot )$ and $h(\cdot )$ are the dynamic motion model function and measurement model function, respectively. ${\epsilon }_k$ is the measurement noise.

The target dynamic motion model function is given by

$$f\left(x_k\left|x_{k-1}\right.\right)=\mathcal{N}\left(x_k;f\left(x_{k-1}\right),Q_k\right)$$

(3)

where $\mathcal{N}\left(x;\widehat{x},P\right)$ denotes a Gaussian probability density function with mean $\widehat{x}$ and covariance matrix $P$. Assume the measurement noise has a heavy tail and can be described by a Student-t distribution as [18]

$$\begin{array}{cc}p({\epsilon }_k)& =\mathrm{St}({\epsilon }_k;{\mathit{0}}_d,R_k,{\nu }_k)\\ & =\frac{\mathrm{\Gamma }\left(\frac{{\nu }_k+d}{2}\right)}{\mathrm{\Gamma }\left(\frac{{\nu }_k}{2}\right)}\frac{1}{{({\nu }_k\pi )}^{\frac{d}{2}}}\frac{1}{{\left|R_k\right|}^{\frac{1}{2}}}{\left(1+\frac{1}{{\nu }_k}{\epsilon }_k^T{R}_k{}^{-1}{\epsilon }_k\right)}^{-\frac{{\nu }_k+d}{2}}\end{array}$$

(4)

where $\mathrm{St}(\cdot )$ denotes the Student-t distribution, ${\mathit{0}}_d$ denotes $d$ dimension zeros vector, $d$ is the dimension of $z_k$, ${\nu }_k$ denotes DOF, $R_k$ is a symmetric, positive defined $d\times d$ covariance matrix and $\mathrm{\Gamma }(\cdot )$ denotes the Gamma function. Different from the Gaussian distribution, the Student-t distribution has a thicker tail away from the mean value. When ${\nu }_k\to \infty$, the Student-t distribution is equivalent to a Gaussian distribution. The probability density function in Equation (4) can be also formulated as a mixture of Gamma distribution and Gaussian distribution denoted as [25]

$$\mathrm{St}({\epsilon }_k;\mathit{0},R_k,{\nu }_k)={\displaystyle \underset{0}{\overset{\infty }{\int }}\mathcal{N}\left({\epsilon }_k;\mathit{0},R_k/u_k\right)}\mathrm{G}\left(u_k;\frac{{\nu }_k}{2},\frac{{\nu }_k}{2}\right)d{u}_k$$

(5)

where $\mathrm{G}\left(\cdot \right)$ denotes the gamma distribution, $u_k$ is an auxiliary variable. The likelihood of the measurement can be written as

$$g_k(z_k\left|x_k,R_k,{\nu }_k\right.)=\mathrm{St}(z_k;h(x_k),R_k,{\nu }_k)$$

(6)

According to Equation (5), we can obtain

$$g_k(z_k\left|x_k,R_k,{\mu }_k\right.)=\mathcal{N}\left(z_k;h(x_k),R_k/u_k\right)$$

(7)

where $u_k$ is distributed as

$$p(u_k\left|{\nu }_k\right.)=\mathrm{G}\left(u_k;{\nu }_k/2,{\nu }_k/2\right)$$

(8)

In this paper, the DOF ${\nu }_k$ and scale matrix $R_k$ are both unknown. Considering the conjugate prior to Gaussian distribution with unknown variance is the IW distribution, we model $R_k$ as the IW distribution. Similarly, the Gamma distribution is assumed to be the conjugate prior to the DOF ${\nu }_k$. Then the joint probability density of $R_k$, $u_k$, and ${\nu }_k$ can be expressed as

$$\begin{array}{cc}p\left(R_k,u_k,{\nu }_k\right)& =\mathrm{IW}\left(R_k;{\upsilon }_k,V_k\right)\mathrm{G}\left(u_k;{\nu }_k/2,{\nu }_k/2\right)\\ & \times \mathrm{G}\left({\nu }_k;{\alpha }_k,{\beta }_k\right)\end{array}$$

(9)

where $\mathrm{IW}(\cdot )$ denotes the IW distribution, ${\upsilon }_k$, $V_k$ are the sufficient statistics of $R_k$, ${\alpha }_k$ and ${\beta }_k$ are the sufficient statistics of ${\nu }_k$.

2.2. PHD Filter

Assume that $X_k=\{x_k^1,\cdots ,x_k^{N_k}\}$ is the multi-target states set and $Z_k=\{z_k^1,\cdots ,z_k^{M_k}\}$ is the measurements set at time $k$. $x_k^n$ is referred to as the nth target state and $z_k^m$ is referred to as the mth measurement. $N_k$ and $M_k$ are the number of targets and measurements, respectively. Then the Bayesian filter recursion can be expressed as

$$p_{k\left|k-1\right.}(X_k\left|Z_{1:k-1}\right.)={\displaystyle \int f(X_k\left|X_{k-1}\right.)}{p}_{k-1}(X_{k-1}\left|Z_{1:k-1}\right.){\mu }_{s}d{X}_{k-1}$$

(10)

$$p_k(X_k\left|Z_{1:k}\right.)=\frac{g_k(Z_k\left|X_k\right.)p_{k\left|k-1\right.}(X_k\left|Z_{1:k-1}\right.)}{{\displaystyle \int g_k(Z_k\left|X_k\right.)p_{k\left|k-1\right.}(X_k\left|Z_{1:k-1}\right.){\mu }_{s}d{X}_k}}$$

(11)

where $g_k(\cdot )$ is the likelihood function. $p_{k\left|k-1\right.}(\cdot )$ and $p_k(\cdot )$ denote multi-target predicted density and multi-target posterior density, respectively. ${\mu }_s$ is the reference measure.

The recursion Equations (10) and (11) involve multiple integrals, which are computationally intractable. To alleviate the computational intractability, the PHD filter provides a sub-optimal solution by propagating the first-order statistical moment of the multiple posterior densities. Without the spawning targets, the predicted and updated equations of the PHD filter are given by [5]

$$D_{k\left|k-1\right.}(x_k)={\displaystyle \int p_{S,k}f(x_k\left|x_{k-1}\right.)D_{k-1}(x_{k-1})d{x}_{k-1}}+{\gamma }_k(x_k)$$

(12)

$$D_k(x_k)=(1-p_{D,k})D_{k\left|k-1\right.}(x_k)+{\displaystyle \sum _{z\in Z_k}\frac{p_{D,k}{g}_k(z\left|x_k\right.)D_{k\left|k-1\right.}(x_k)}{{\kappa }_k(z)+{\displaystyle \int p_{D,k}{g}_k(z\left|x_k\right.)D_{k\left|k-1\right.}(x_k)d{x}_k}}}$$

(13)

where $D_{k\left|k-1\right.}(x_k)$ and $D_{k\left|k\right.}(x_k)$ are the predicted intensity and posterior intensity, respectively. ${\gamma }_k(x_k)$ is the birth targets intensity, $p_{S,k}$ is the target survival probability, $p_{D,k}$ is the detection probability. ${\kappa }_k(z)=\lambda c(z)$ is the clutter intensity with average clutter number $\lambda$ and spatial distribution $c(z)$.

2.3. SMC-PHD Filter

Assume that the multi-target posterior intensity can be denoted by weighted particles set at time $k-1$, i.e., $D_{k-1}(x_{k-1}\left|Z_{1:k-1}\right.)={\displaystyle {\sum }_{i=1}^{L_{k-1}}{\omega }_{k-1}^{(i)}\delta (x_{k-1}-x_{k-1}^{(i)})}$. $L_{k-1}$ denotes the particle number at time $k-1$, $x_{k-1}^{(i)}$ is the ith particle state and ${\omega }_{k-1}^{(i)}$ is the corresponding weight. Then, the SMC-PHD filter can be summarized by the following steps.

Step 1 Prediction

$$x_k^{(i)}\sim \left\{\begin{array}{l}{q}_k(\cdot \left|x_{k-1}^{(i)},Z_k\right.),i=1,\dots ,L_{k-1}\hfill \\ p_k(\cdot \left|Z_k\right.),i=L_{k-1}+1,\dots ,L_{k-1}+J_k\hfill \end{array}\right.$$

(14)

$${\omega }_{k\left|k-1\right.}^{(i)}\sim \left\{\begin{array}{l}{\omega }_{k-1}^{(i)}\frac{p_{S,k}{f}_{k\left|k-1\right.}(x_k^{(i)}\left|x_{k-1}^{(i)}\right.)}{q_k(x_k^{(i)}\left|x_{k-1}^{(i)},Z_k\right.)},i=1,\dots ,L_{k-1}\hfill \\ \frac{1}{J_k}\frac{{\gamma }_k(x_k^{(i)})}{p_k(x_k^{(i)}\left|Z_k\right.)},i=L_{k-1}+1,\dots ,L_{k-1}+J_k\hfill \end{array}\right.$$

(15)

$$D_{k\left|k-1\right.}(x_k\left|Z_{1:k-1}\right.)={\displaystyle \sum _{i=1}^{L_{k-1}+J_k}{\omega }_{k\left|k-1\right.}^{(i)}\delta (x_k-x_k^{(i)})}$$

(16)

where $q_k\left(\cdot \left|x_{k-1}^{(i)},Z_k\right.\right)$ is the proposal density of survival targets and $p_k\left(\cdot \left|Z_k\right.\right)$ is the proposal density of birth targets. $J_k$ is the particles number of birth targets.

Step 2 Update

$${\omega }_k^{(i)}=(1-p_{D,k}){\omega }_{k\left|k-1\right.}^{(i)}+{\displaystyle \sum _{z\in Z_k}\frac{p_{D,k}{g}_k(z\left|x_k^{(i)}\right.)}{{\kappa }_k(z)+C_k(z)}}{\omega }_{k\left|k-1\right.}^{(i)}$$

(17)

$$C_k(z)={\displaystyle \sum _{j=1}^{L_{k-1}+J_k}{p}_{D,k}{g}_k(z\left|x_k^{(i)}\right.){\omega }_{k\left|k-1\right.}^{(i)}}$$

(18)

$$D_k(x_k\left|Z_{1:k}\right.)={\displaystyle \sum _{i=1}^{L_{k-1}+J_k}{\omega }_k^{(i)}\delta (x_k-x_k^{(i)})}$$

(19)

Step 3 Resampling

Target number is estimated by ${\widehat{N}}_k=\mathrm{ceil}\left({\displaystyle {\sum }_{i=1}^{L_{k-1}+J_k}{\omega }_k^{(i)}}\right)$, and resample ${\left\{x_k^{(i)},{\omega }_k^{(i)}\right\}}_{i=1}^{L_{k-1}+J_k}$. $\mathrm{ceil}(x)$ denotes the smallest integer more than $x$, $L_k=N_{particle}\cdot {\widehat{N}}_k$, $N_{particle}$ is the number of particles per expected target.

Step 4 State extraction

Extract the target states $X_k=\{{\widehat{x}}_k^1,\cdots ,{\widehat{x}}_k^{{\widehat{N}}_k}\}$ from the resampling particles set, where ${\widehat{x}}_k^1,\cdots ,{\widehat{x}}_k^{{\widehat{N}}_k}$ denote the estimated target states.

3. The Robust SMC-PHD Filter

3.1. VB Approximation

The aforementioned SMC-PHD filter performs well under Gaussian measurement noise. However, the outliers in measurements will have a negative effect on the SMC-PHD filter. To solve the problem, the Student-t distribution is introduced to model the heavy-tailed measurement noise. Since the parameters of the Student-t distribution are unknown, the joint posterior probability density of target state $x_k$, scale matrix $R_k$, auxiliary variable $u_k$, and DOF ${\nu }_k$ need to be estimated. Therefore, the augmented target state can be expressed as

$${\ddot{x}}_k=\left(x_k,R_k,u_k,{\nu }_k\right)$$

(20)

Since there is no explicit analytical solution existing for the joint posterior density $p({\ddot{x}}_k\left|Z_{1:k}\right.)$, we utilize the VB inference [9] to obtain the approximate solution. As for an approximate method, VB inference is widely used in searching for a factored approximate posterior density as follows

$$\begin{array}{cc}p({\ddot{x}}_k\left|Z_{1:k}\right.)& =p(x_k,R_k,u_k,{\nu }_k\left|Z_{1:k}\right.)\\ & \approx q({\ddot{x}}_k\left|Z_{1:k}\right.)\triangleq q_x(x_k)q_R(R_k)q_u(u_k)q_v({\nu }_k)\end{array}$$

(21)

where $q({\ddot{x}}_k\left|Z_{1:k}\right.)$ is the approximate joint posterior density, $q_x(\cdot )$, $q_R(\cdot )$, $q_u(\cdot )$,and $q_{\nu }(\cdot )$ are the approximate posterior densities. Equation (21) has the implicit assumption that the target state and the Student-t distribution parameters are independent. By minimizing the Kullback–Leibler (KL) divergence between $p({\ddot{x}}_k\left|Z_{1:k}\right.)$ and $q({\ddot{x}}_k\left|Z_{1:k}\right.)$, the estimate ${\widehat{q}}_x(x_k)$, ${\widehat{q}}_R(R_k)$, ${\widehat{q}}_u(u_k)$, and ${\widehat{q}}_{\nu }({\nu }_k)$ can be easily obtained by

$$\left\{{\widehat{q}}_x,{\widehat{q}}_R,{\widehat{q}}_u,{\widehat{q}}_{\nu }\right\}=\underset{q_x,q_R,q_u,q_{\nu }}{\arg \min }KL\left[q(\left\{{\ddot{x}}_k\right\}\left|Z_{1:k}\right.)‖p(\left\{{\ddot{x}}_k\right\}\left|Z_{1:k}\right.)\right]$$

(22)

where $KL(q(x)\left|p(x)\right.)\triangleq {\displaystyle \int q(x)\log \left(\frac{q(x)}{p(x)}\right)dx}$. The optimal solution of Equation (22) can be obtained by [9,21,26]

$$\ln \widehat{q}(\varphi )=E_{{\ddot{x}}_k^{(-\varphi )}}\left[\ln p\left({\ddot{x}}_k,Z_k\left|Z_{1:k-1}\right.\right)\right]+c_{\varphi }$$

(23)

where $\varphi$ denotes an arbitrary element in ${\ddot{x}}_k$, $E[\cdot ]$ denotes the expectation, ${\ddot{x}}_k^{(-\varphi )}$ denotes the elements set in ${\ddot{x}}_k$ except for $\varphi$, and $c_{\varphi }$ is an irrelevant constant for element $\varphi$.

3.2. The Update of Unknown Student-t Distribution Parameters

Based on the assumption that the target state is independent of the Student-t distribution parameters, the joint posterior density of the augmented target state can be factorized as

$$p({\ddot{x}}_k\left|Z_{1:k}\right.)=p(R_k,u_k,{\nu }_k\left|Z_{1:k}\right.)p(x_k\left|Z_{1:k}\right.)$$

(24)

where the unknown Student-t distribution parameters are integrated out. Based on (24), the target state can be estimated by particle samples and the Student-t distribution parameters can be updated by the VB approximation method.

In this subsection, we derive the update of Student-t distribution parameters for each particle. At time $k$, let ${\ddot{x}}_k=\{{\ddot{x}}_k^{(i)}{\}}_{i=1}^N=\{x_k^{(i)},R_k^{(i)},u_k^{(i)},{\nu }_k^{(i)}{\}}_{i=1}^N$ is the augmented target state where $R_k^{(i)}$, $u_k^{(i)}$, ${\nu }_k^{(i)}$ are the Student-t distribution parameters with respect to the ith particle, and $N$ is the number of particles. Assume the predicted density of ${\ddot{x}}_k^{(i)}$ at time $k$ can be expressed as

$$\begin{array}{cc}p({\ddot{x}}_k^{(i)}\left|Z_{1:k-1}\right.)& ={\omega }_{k\left|k-1\right.}^{(i)}\delta (x_k-x_k^{(i)})\mathrm{IW}(R_k^{(i)};{\upsilon }_{k\left|k-1\right.}^{(i)},V_{k\left|k-1\right.}^{(i)})\\ & \times \mathrm{G}\left(u_k^{(i)};\frac{{\nu }_k^{(i)}}{2},\frac{{\nu }_k^{(i)}}{2}\right)\mathrm{G}({\nu }_k^{(i)};{\alpha }_{k\left|k-1\right.}^{(i)},{\beta }_{k\left|k-1\right.}^{(i)})\end{array}$$

(25)

where ${\upsilon }_{k\left|k-1\right.}^{(i)}$, $V_{k\left|k-1\right.}^{(i)}$, ${\alpha }_{k\left|k-1\right.}^{(i)}$ and ${\beta }_{k\left|k-1\right.}^{(i)}$ are the predicted sufficient statistics with respect to $R_k^{(i)}$ and ${\nu }_k^{(i)}$. Then the joint density of $p({\ddot{x}}_k^{(i)},z_k\left|Z_{1:k-1}\right.)$ can be expressed as

$$\begin{array}{cc}p({\ddot{x}}_k^{(i)},z_k\left|Z_{1:k-1}\right.)=& {\omega }_{k\left|k-1\right.}^{(i)}\delta (x_k-x_k^{(i)})\mathcal{N}\left(z_k;h(x_k^{(i)}),\frac{R_k^{(i)}}{u_k^{(i)}}\right)\\ & \times \mathrm{IW}(R_k^{(i)};{\upsilon }_{k\left|k-1\right.}^{(i)},V_{k\left|k-1\right.}^{(i)})\mathrm{G}\left(u_k^{(i)};\frac{{\nu }_k^{(i)}}{2},\frac{{\nu }_k^{(i)}}{2}\right)\\ & \times \mathrm{G}({\nu }_k^{(i)};{\alpha }_{k\left|k-1\right.}^{(i)},{\beta }_{k\left|k-1\right.}^{(i)})\end{array}$$

(26)

where $z_k$ is the measurement with respect to the target represented by particle $x_k^{(i)}$. Based on Equation (26), $\ln p({\ddot{x}}_k^{(i)},z_k\left|Z_{1:k-1}\right.)$ can be expressed as

$$\begin{array}{cc}\ln p({\ddot{x}}_k^{(i)},z_k\left|Z_{1:k-1}\right.)& =-\frac{1}{2}\ln \left|R_k^{(i)}\right|+\frac{d}{2}\ln u_k^{(i)}-\frac{u_k^{(i)}}{2}tr\left({(R_k^{(i)})}^{-1}{\Xi }_k^{(i)}\right)-\frac{{\upsilon }_{k\left|k-1\right.}^{(i)}}{2}\ln \left|R_k^{(i)}\right|\\ & -\frac{1}{2}tr\left({(R_k^{(i)})}^{-1}{V}_{k\left|k-1\right.}^{(i)}\right)+\frac{{\nu }_k^{(i)}}{2}\ln \frac{{\nu }_k^{(i)}}{2}-\ln \Gamma (\frac{{\nu }_k^{(i)}}{2})+(\frac{{\nu }_k^{(i)}}{2}-1)\ln u_k^{(i)}\\ & -\frac{{\nu }_k^{(i)}}{2}{u}_k^{(i)}+({\alpha }_{k\left|k-1\right.}^{(i)}-1)\ln {\nu }_k^{(i)}-{\beta }_{k\left|k-1\right.}^{(i)}{\nu }_k^{(i)}+c_{{\ddot{x}}_k^{(i)}}\end{array}$$

(27)

where ${\Xi }_k^{(i)}=(z_k-h(x_k^{(i)})){(z_k-h(x_k^{(i)}))}^T$ denotes the innovation covariance matrix. $tr(\cdot )$ denotes the trace operation. $c_{{\ddot{x}}_k^{(i)}}$ denotes an irrelevant constant for variables $x_k^{(i)}$, $R_k^{(i)}$, $u_k^{(i)}$, and ${\nu }_k^{(i)}$.

Letting $\varphi =R_k^{(i)}$ and using Equation (27) in Equation (23), we can obtain

$$\ln \widehat{q}(R_k^{(i)})=-\frac{{\upsilon }_{k\left|k-1\right.}^{(i)}+1}{2}\ln \left|R_k^{(i)}\right|-\frac{1}{2}\left\{tr\left({(R_k^{(i)})}^{-1}\left[V_{k\left|k-1\right.}^{(i)}+u_k^{(i)}{\Xi }_k^{(i)}\right]\right)\right\}+c_{R_k^{(i)}}$$

(28)

The posterior distribution $\widehat{q}(R_k^{(i)})$ of $R_k^{(i)}$ is an IW distribution, i.e.,

$$\widehat{q}(R_k^{(i)})=\mathrm{IW}(R_k^{(i)};{\upsilon }_k^{(i)},V_k^{(i)})$$

(29)

with parameters

$${\upsilon }_k^{(i)}={\upsilon }_{k\left|k-1\right.}^{(i)}+1$$

(30)

$$V_k^{(i)}=V_{k\left|k-1\right.}^{(i)}+E[u_k^{(i)}]{\Xi }_k^{(i)}$$

(31)

Letting $\varphi ={\nu }_k^{(i)}$ and using Equation (27) in Equation (23), we can obtain

$$\begin{array}{cc}\ln \widehat{q}({\nu }_k^{(i)})& =\frac{{\nu }_k^{(i)}}{2}\ln \frac{{\nu }_k^{(i)}}{2}-\ln \mathrm{\Gamma }(\frac{{\nu }_k^{(i)}}{2})+\frac{{\nu }_k^{(i)}}{2}\ln u_k^{(i)}\\ & -\frac{{\nu }_k^{(i)}}{2}{u}_k^{(i)}+({\alpha }_{k\left|k-1\right.}^{(i)}-1)\ln {\nu }_k^{(i)}-{\beta }_{k\left|k-1\right.}^{(i)}{\nu }_k^{(i)}+c_{{\nu }_k^{(i)}}\end{array}$$

(32)

Using Stirling’s approximation in [27]

$$\ln \mathrm{\Gamma }(\frac{{\nu }_k^{(i)}}{2})\approx \frac{{\nu }_k^{(i)}-1}{2}\ln \frac{{\nu }_k^{(i)}}{2}-\frac{{\nu }_k^{(i)}}{2}$$

(33)

And $\ln \widehat{q}({\nu }_k^{(i)})$ can be expressed as

$$\begin{array}{cc}\ln \widehat{q}({\nu }_k^{(i)})& =({\alpha }_{k\left|k-1\right.}^{(i)}+\frac{1}{2}-1)\ln {\nu }_k^{(i)}\\ & -({\beta }_{k\left|k-1\right.}^{(i)}-\frac{1}{2}(1+E[\ln u_k^{(i)}]-E[u_k^{(i)}])){\nu }_k^{(i)}\end{array}$$

(34)

The posterior distribution $\widehat{q}({\nu }_k^{(i)})$ of ${\nu }_k^{(i)}$ is a Gamma distribution, i.e.,

$$\widehat{q}({\nu }_k^{(i)})=\mathrm{G}({\nu }_k^{(i)};{\alpha }_k^{(i)},{\beta }_k^{(i)})$$

(35)

with parameters

$${\alpha }_k^{(i)}={\alpha }_{k\left|k-1\right.}^{(i)}+\frac{1}{2}$$

(36)

$${\beta }_k^{(i)}={\beta }_{k\left|k-1\right.}^{(i)}-\frac{1}{2}(1+E[\ln u_k^{(i)}]-E[u_k^{(i)}])$$

(37)

Letting $\varphi =u_k^{(i)}$ and using Equation (27) in Equation (23), we can obtain

$$\ln \widehat{q}(u_k^{(i)})=(\frac{{\nu }_k^{(i)}+d}{2}-1)\ln u_k^{(i)}-\frac{1}{2}[{\nu }_k^{(i)}+tr\left({(R_k^{(i)})}^{-1}{\Xi }_k^{(i)}\right)]u_k^{(i)}+c_{u_k^{(i)}}$$

(38)

The posterior distribution $\widehat{q}(u_k^{(i)})$ of $u_k^{(i)}$ is a Gamma distribution, i.e.,

$$\widehat{q}(u_k^{(i)})=\mathrm{G}(u_k^{(i)};\frac{{\nu }_k^{(i)}+d}{2},\frac{1}{2}[{\nu }_k^{(i)}+tr\left({(R_k^{(i)})}^{-1}{\Xi }_k^{(i)}\right)])$$

(39)

Then we can obtain

$$E[u_k^{(i)}]=a_k^{(i)}/b_k^{(i)}$$

(40)

$$E[\ln u_k^{(i)}]=\phi (a_k^{(i)})-\ln b_k^{(i)}$$

(41)

where $a_k^{(i)}=\frac{E[{\nu }_k^{(i)}]+d}{2}$, $b_k^{(i)}=\frac{1}{2}\left\{E[{\nu }_k^{(i)}]+tr\left({(R_k^{(i)})}^{-1}{\Xi }_k^{(i)}\right)\right\}$. Since the parameters in Equations (30), (31), (36), (37), (40) and (41) are coupled, it can be solved iteratively by the fixed point method in [21].

3.3. The Implementation of RSMC-PHD Filter

In this subsection, we introduce the Student-t distribution into the SMC-PHD filter and present an implementation of the robust SMC-PHD filter with unknown heavy-tailed measurement noise. Assume that the multi-target posterior intensity is represented by a weighted particles set ${\{{\omega }_{k-1}^{(i)},x_{k-1}^{(i)},R_{k-1}^{(i)},u_{k-1}^{(i)},{\nu }_{k-1}^{(i)}\}}_{i=1}^{L_{k-1}}$ at time $k-1$, i.e.,

$$\begin{array}{cc}{D}_{k-1}({\ddot{x}}_{k-1}\left|Z_{1:k-1}\right.)& ={\displaystyle \sum _{i=1}^{L_{k-1}}{\omega }_{k-1}^{(i)}\delta (x_{k-1}-x_{k-1}^{(i)})}\mathrm{IW}(R_{k-1}^{(i)};{\upsilon }_{k-1}^{(i)},V_{k-1}^{(i)})\\ & \times \mathrm{G}\left(u_k^{(i)};\frac{{\nu }_{k-1}^{(i)}}{2},\frac{{\nu }_{k-1}^{(i)}}{2}\right)\mathrm{G}({\nu }_{k-1}^{(i)};{\alpha }_{k-1}^{(i)},{\beta }_{k-1}^{(i)})\end{array}$$

(42)

Assume $R_k$, $u_k$, and ${\nu }_k$ are independent and the predicted intensity at time $k$ can be expressed as

$$\begin{array}{cc}{D}_{k\left|k-1\right.}({\ddot{x}}_k\left|Z_{1:k-1}\right.)& ={\displaystyle \sum _{i=1}^{L_{k-1}+J_k}{\omega }_{k\left|k-1\right.}^{(i)}\delta (x_k-x_k^{(i)})}\mathrm{G}\left(u_k^{(i)};\frac{{\nu }_k^{(i)}}{2},\frac{{\nu }_k^{(i)}}{2}\right)\\ & \times \mathrm{IW}(R_k^{(i)};{\upsilon }_{k\left|k-1\right.}^{(i)},V_{k\left|k-1\right.}^{(i)})\mathrm{G}({\nu }_k^{(i)};{\alpha }_{k\left|k-1\right.}^{(i)},{\beta }_{k\left|k-1\right.}^{(i)})\end{array}$$

(43)

where ${\omega }_{k\left|k-1\right.}^{(i)}$ and $x_k^{(i)}$ are given in Equations (14) and (15). There are $L_{k-1}$ particles for survival targets and $J_k$ particles for birth targets. Generally, the dynamic of measurement noise is unknown and we use the similar heuristic dynamic in [9] for the predicted sufficient statistics. For the particles of survival targets $i=1,2,\dots ,L_{k-1}$, the predicted sufficient statistics are given by

$${\upsilon }_{k\left|k-1\right.}^{(i)}=\rho ({\upsilon }_{k-1}^{(i)}-d-1)+d+1,V_{k\left|k-1\right.}^{(i)}=\rho V_{k-1}^{(i)}$$

(44)

$${\alpha }_{k\left|k-1\right.}^{(i)}=\rho {\alpha }_{k-1}^{(i)},{\beta }_{k\left|k-1\right.}^{(i)}=\rho {\beta }_{k-1}^{(i)}$$

(45)

where $\rho \in (0,1]$ is the decreasing factor. For the particles of birth targets $i=L_{k-1}+1,L_{k-1}+2,\dots ,L_{k-1}+J_k$, the predicted sufficient statistics are given by

$${\upsilon }_{k\left|k-1\right.}^{(i)}={\upsilon }_{\gamma ,k}^{(i)},V_{k\left|k-1\right.}^{(i)}=V_{\gamma ,k}^{(i)}$$

(46)

$${\alpha }_{k\left|k-1\right.}^{(i)}=\rho {\alpha }_{k-1}^{(i)},{\beta }_{k\left|k-1\right.}^{(i)}=\rho {\beta }_{k-1}^{(i)}$$

(47)

where ${\upsilon }_{\gamma ,k}^{(i)}$, $V_{\gamma ,k}^{(i)}$, ${\alpha }_{\gamma ,k}^{(i)}$ and ${\beta }_{\gamma ,k}^{(i)}$ are sufficient statistics with respect to the birth targets.

The posterior intensity can be expressed as

$$\begin{array}{cc}{D}_k({\ddot{x}}_k\left|Z_{1:k}\right.)& ={\displaystyle \sum _{i=1}^{L_{k-1}+J_k}{\omega }_k^{(i)}\delta (x_k-x_k^{(i)})}\mathrm{G}\left(u_k^{(i)};\frac{{\nu }_k^{(i)}}{2},\frac{{\nu }_k^{(i)}}{2}\right)\\ & \times \mathrm{IW}(R_k^{(i)};{\upsilon }_k^{(i)},V_k^{(i)})\mathrm{G}({\nu }_k^{(i)};{\alpha }_k^{(i)},{\beta }_k^{(i)})\end{array}$$

(48)

where ${\omega }_k^{(i)}$ is given in (17). It is worth mentioning that the likelihood function is $g_k(z\left|{\ddot{x}}_k\right.)=\mathrm{St}(z;h(x_k^{(i)}),R_k^{(i)},{\nu }_k^{(i)})$. Then, the sufficient statistics ${\upsilon }_k^{(i)}$, $V_k^{(i)}$, ${\alpha }_k^{(i)}$, and ${\beta }_k^{(i)}$ are updated by

$$R_k^{(i)(j)}=V_k^{(i)(j-1)}/({\upsilon }_k^{(i)}-d-1)$$

(49)

$${\tilde{R}}_k^{(i)(j)}=R_k^{(i)(j)}/u_k^{(i)(j-1)}$$

(50)

$${\gamma }_k^{(i)(j)}=tr({\Xi }_k^{(i)}/{\tilde{R}}_k^{(i)(j)})$$

(51)

$$a_k^{(i)(j)}=({\nu }_k^{(i)(j-1)}+d)/2$$

(52)

$$b_k^{(i)(j)}=({\nu }_k^{(i)(j-1)}+{\gamma }_k^{(i)(j)})/2$$

(53)

$$u_k^{(i)(j)}=a_k^{(i)(j)}/b_k^{(i)(j)}$$

(54)

$${\beta }_k^{(i)(j)}={\beta }_{k\left|k-1\right.}^{(i)}-\frac{1}{2}-\frac{\phi (a_k^{(i)(j)})-\ln (b_k^{(i)(j)})}{2}+\frac{u_k^{(i)(j)}}{2}$$

(55)

$$V_k^{(i)(j)}=V_{k\left|k-1\right.}^{(i)}+{\Xi }_k^{(i)}{u}_k^{(i)(j)}$$

(56)

$${\nu }_k^{(i)(j)}={\alpha }_k^{(i)}/{\beta }_k^{(i)(j)}$$

(57)

where ${(\cdot )}^{(i)(j)}$ denotes the update of parameter ${(\cdot )}^{(i)}$ at the jth iteration. The initial values are given by $V_k^{(i)(\mathit{0})}=V_{k\left|k-1\right.}^{(i)}$, ${\upsilon }_k^{(i)}={\upsilon }_{k\left|k-1\right.}^{(i)}+1$, ${\alpha }_k^{(i)}={\alpha }_{k\left|k-1\right.}^{(i)}+1/2$, ${\beta }_k^{(i)(\mathit{0})}={\beta }_{k\left|k-1\right.}^{(i)}$, ${\nu }_k^{(i)(\mathit{0})}=$${\alpha }_k^{(i)}/{\beta }_k^{(i)(\mathit{0})}$.

It should be noted that the innovation covariance matrix ${\Xi }_k^{(i)}=(z_k-h(x_k^{(i)})){(z_k-h(x_k^{(i)}))}^T$ for particle $x_k^{(i)}$ is needed in each iteration. Since there is no explicit association between measurements and targets, the measurement origin is ambiguous. Considering the fact that the measurement with larger likelihood value is more likely generated by a potential target, we choose the measurement with the largest likelihood value to calculate the innovation covariance matrix as follows

$${\Xi }_k^{(i)}=(z_k^l-h(x_k^{(i)})){(z_k^l-h(x_k^{(i)}))}^T$$

(58)

where $l=\underset{p\in I}{\arg \max }{g}_k^p({\ddot{x}}_k^{(i)})$, $I=\left\{1,\dots ,\left|Z_k\right|\right\}$. $g_k^p({\ddot{x}}_k^{(i)})$ is the normalized likelihood function value given by

$$g_k^p({\ddot{x}}_k^{(i)})=\frac{p_{D,k}{g}_k(z_k^p\left|{\ddot{x}}_k^{(i)}\right.)}{{\kappa }_k(z_k^p)+C_k(z_k^p)}$$

(59)

$$C_k(z_k^p)={\displaystyle \sum _{j=1}^{L_{k-1}+J_k}{p}_{D,k}{g}_k(z_k^p\left|{\ddot{x}}_k^{(i)}\right.){\omega }_{k\left|k-1\right.}^{(i)}}$$

(60)

3.4. The Modification of Particles Weight

The estimate of target number relies on the updated weights of particles in the SMC-PHD filter. In our paper, the weight is updated by Equation (17) for each particle. The likelihood function is $g_k(z\left|{\ddot{x}}_k\right.)=\mathrm{St}(z;h(x_k^{(i)}),R_k^{(i)},{\nu }_k^{(i)})$, and the updated weight of the ith particle is given by

$${\omega }_k^{(i)}=(1-p_{D,k}){\omega }_{k\left|k-1\right.}^{(i)}+{\displaystyle \sum _{p=1}^{\left|Z_k\right|}{g}_k^p({\ddot{x}}_k^{(i)})}{\omega }_{k\left|k-1\right.}^{(i)}\propto \mathrm{St}(z_k^p;h(x_k^{(i)}),R_k^{(i)},{\nu }_k^{(i)})$$

(61)

Figure 1 gives the illustration of the measurement likelihood function. We can see that the Student-t distribution has a heavy-tailed characteristic and it can deal with outliers well by increasing the weights of outliers. However, the weight of clutter will increase at the same time. Hence, some clutter may be regarded as target-generated measurements, resulting in the overestimate of the target number at the current time.

To solve this problem, a strategy is adopted to modify the updated weight in this subsection. In the PHD filter, one target can produce at most one measurement and only one measurement can be attributed to one target per time step [4], named as the ‘one-to-one’ assumption. For each particle, the normalized measurement likelihood value $g_k^p({\ddot{x}}_k^{(i)})$ can be obtained by Equations (59) and (60). Let $J=\{p=1,\dots ,\left|Z_k\right|\left|g_k^p({\ddot{x}}_k^{(i)})\ge \sigma \right.\}$, $\sigma$ is a preset threshold. If $\left|J\right|\ge 2$, it indicates there might be multiple measurements associated with particle $x_k^{(i)}$. In this case, the $g_k^p({\ddot{x}}_k^{(i)})$ is modified as

$${\tilde{g}}_k^p({\ddot{x}}_k^{(i)})=\left\{\begin{array}{c}{g}_k^p({\ddot{x}}_k^{(i)}),p=l or p\notin J\\ g_k^q({\ddot{x}}_k^{(i)}),p\ne l \& p\in J\end{array}\right.$$

(62)

where $q=\underset{p\in I}{\arg \min }{g}_k^p({\ddot{x}}_k^{(i)})$, $\&$ is a logic and operation and $or$ is a logic or operation. In our paper, the measurement with the largest normalized likelihood value is regarded as the corresponding measurement generated from the target corresponding to particle $x_k^{(i)}$, and other measurements are regarded as clutter. The proposed strategy retains the maximum likelihood value and penalizes the other larger likelihood value to hold the ‘one-to-one’ assumption.

Then the modified weight can be expressed as

$${\omega }_k^{(i)}=(1-p_{D,k}){\omega }_{k\left|k-1\right.}^{(i)}+{\displaystyle \sum _{p=1}^{\left|Z_k\right|}{\tilde{g}}_k^p({\ddot{x}}_k^{(i)})}{\omega }_{k\left|k-1\right.}^{(i)}$$

(63)

The robust SMC-PHD filter is summarized in Algorithm 1.

Remark 1:

Although the heavy-tailed characteristic of the Student-t distribution can deal with outliers, the target might still miss when the weights of outliers are lower than the extracted threshold and the target number would be underestimated. In our paper, the weight modification of the particles is mainly aimed to deal with the target number overestimate problem caused by clutter at the current time.

Algorithm 1. The robust SMC-PHD filter.

Step 1 Prediction for $i=1,2,\dots ,L_{k-1}$     $x_k^{(i)}\sim q_k(\cdot \left|x_{k-1}^{(i)},Z_k\right.),{\omega }_{k\left|k-1\right.}^{(i)}\sim {\omega }_{k-1}^{(i)}\frac{p_{S,k}{f}_{k\left|k-1\right.}(x_k^{(i)}\left|x_{k-1}^{(i)}\right.)}{q_k(x_k^{(i)}\left|x_{k-1}^{(i)},Z_k\right.)}$     $\begin{array}{l}{\upsilon }_{k\left|k-1\right.}^{(i)}=\rho ({\upsilon }_{k-1}^{(i)}-d-1)+d+1,V_{k\left|k-1\right.}^{(i)}=\rho V_{k-1}^{(i)}\\ {\alpha }_{k\left|k-1\right.}^{(i)}=\rho {\alpha }_{k-1}^{(i)},{\beta }_{k\left|k-1\right.}^{(i)}=\rho {\beta }_{k-1}^{(i)}\end{array}$ end for for $i=L_{k-1}+1,L_{k-1}+2,\dots ,L_{k-1}+J_k$     $x_k^{(i)}\sim p_k(\cdot \left|Z_k\right.),{\omega }_{k\left|k-1\right.}^{(i)}\sim \frac{1}{J_k}\frac{{\gamma }_k(x_k^{(i)})}{p_k(x_k^{(i)}\left|Z_k\right.)}$     ${\upsilon }_{k\left|k-1\right.}^{(i)}={\upsilon }_{\gamma ,k}^{(i)},V_{k\left|k-1\right.}^{(i)}=V_{\gamma ,k}^{(i)}$     ${\alpha }_{k\left|k-1\right.}^{(i)}={\alpha }_{\gamma ,k}^{(i)},{\beta }_{k\left|k-1\right.}^{(i)}={\beta }_{\gamma ,k}^{(i)}$ end for

Step 2 Update For $i=1,2,\dots ,L_{k-1}+J_k$         $\begin{array}{l}{R}_k^{(i)}=V_{k\left|k-1\right.}^{(i)}/({\upsilon }_{k\left|k-1\right.}^{(i)}-d-1),{\nu }_k^{(i)}={\alpha }_{k\left|k-1\right.}^{(i)}/{\beta }_{k\left|k-1\right.}^{(i)}\\ g_k(z\left|{\ddot{x}}_k^{(i)}\right.)=\mathrm{St}(z;h(x_k^{(i)}),R_k^{(i)},{\nu }_k^{(i)})\end{array}$         $\begin{array}{l}{g}_k^p({\ddot{x}}_k^{(i)})=\frac{p_{D,k}{g}_k(z_k^p\left|{\ddot{x}}_k^{(i)}\right.)}{{\kappa }_k(z_k^p)+C_k(z_k^p)}\\ C_k(z_k^p)={\displaystyle \sum _{j=1}^{L_{k-1}+J_k}{p}_{D,k}{g}_k(z_k^p\left|{\ddot{x}}_k^{(i)}\right.){\omega }_{k\left|k-1\right.}^{(i)}}\end{array}$         $I=\left\{1,\dots ,\left|Z_k\right|\right\},J=\left\{p=1,\dots ,\left|Z_k\right|\left|g_k^p({\ddot{x}}_k^{(i)})\ge \sigma \right.\right\}$         $l=\underset{p\in I}{\arg \max }{g}_k^p({\ddot{x}}_k^{(i)}),q=\underset{p\in I}{\arg \min }{g}_k^p({\ddot{x}}_k^{(i)})$         ${\tilde{g}}_k^p({\ddot{x}}_k^{(i)})=\left\{\begin{array}{c}{g}_k^p({\ddot{x}}_k^{(i)}),p=l or p\notin J\\ g_k^q({\ddot{x}}_k^{(i)}),p\ne l \& p\in J\end{array}\right.$         ${\omega }_k^{(i)}=(1-p_{D,k}){\omega }_{k\left|k-1\right.}^{(i)}+{\displaystyle \sum _{p=1}^{\left|Z_k\right|}{\tilde{g}}_k^p({\ddot{x}}_k^{(i)})}{\omega }_{k\left|k-1\right.}^{(i)}$ Perform iteration initialization, set the iteration number $\tau$, and update the parameters of measurement noise by Equations (49)–(60). end for

Step 3 Resampling Step 4 State extraction     Step 3 and step 4 are the same as that in the standard SMC-PHD filter.

4. Simulation Results

To reveal the performance of the proposed robust SMC-PHD (RSMC-PHD) filter, we compare the Gaussian (Normal) Gamma inverse Wishart Gamma PHD (NGIWG-PHD) filter [23] and the SMC-PHD filter [6] with the proposed filter. For each simulation, 100 Monte Carlo trials are performed. The target number and optimal sub-pattern assignment (OSPA) [28] distance are chosen as the metric, where the parameters of OSPA are $p=2$ and $c=100$.

4.1. The Choice of Threshold $\sigma$

In Section 3.4, threshold $\sigma$ is utilized to distinguish the target-generated measurements from clutter. In the SMC-PHD filter, the underlying target relies on the updated weight of each particle, and the likelihood value in Equation (59) can directly reflect the contribution of a given measurement or clutter to the underlying target. In most cases, the target-generated measurements usually have the largest likelihood value for each particle. If $\sigma$ is larger, some clutter might be still regarded as target-generated measurements, and the overestimate of the target number cannot be solved. By contrast, if $\sigma$ is smaller, the likelihood values of some clutter which contribute a small proportion to the updated weight will be modified. However, this modification will have no obvious effect on the updated weight of the particle except for unnecessary complexity. As we all know, $N_{particle}$ is the number of particles for each expected target in the SMC-PHD filter and it is used to extract the weighted particles to approximate the posterior intensity. To some extent, $1/N_{particle}$ can approximately reflect the average contribution of each particle for an underlying target. Therefore, we choose $\sigma =1/N_{particle}$ as the threshold.

4.2. Linear Scenario

The target state of each target is $x_k={[x_k,{\dot{x}}_k,y_k,{\dot{y}}_k]}_{}^T$, consisting of the position ${[x_k,y_k]}^T$ and velocity ${[{\dot{x}}_k,{\dot{y}}_k]}_{}^T$. The target motion model is given by

$$x_k=F{x}_{k-1}+G{w}_k$$

(64)

where

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

(65)

$$G={\left[\begin{array}{cccc}\Delta & 1& 0& 0\\ 0& 0& \Delta & 1\end{array}\right]}^T$$

(66)

$w_k\sim \mathcal{N}(\cdot ;0_2,Q_k)$, $Q_k=\mathrm{diag}({[{\sigma }_w^2,{\sigma }_w^2]}^T)$, ${\sigma }_w=5 { \mathrm{m}/\mathrm{s}}^2$ is the standard deviation of the process noise, $\Delta =1 \mathrm{s}$ is the sampling period.

The observation model is formulated by

$$z_k=H{x}_k+{\epsilon }_k$$

(67)

where

$$H=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\end{array}\right]$$

(68)

The heavy-tailed measurement noise corrupted by outliers is represented as [21,23]

$${\epsilon }_k\sim \left\{\begin{array}{l}\mathcal{N}({\epsilon }_k;{\mathit{0}}_2,R_{\mathit{0}}) 1-\eta \\ \mathcal{N}({\epsilon }_k;{\mathit{0}}_2,R_n) \eta \end{array}\right.$$

(69)

where $\eta$ denotes the probability of the measurement noise outlier. $R_0=\mathrm{diag}({[{10}^2,{10}^2]}^T)$ is the nominal noise variance, and $R_n=20R_0$.

The model of the birth target is a Poisson RFS with intensity

$${\gamma }_k(x)={\displaystyle {\sum }_{j=1}^{3}0.1\mathcal{N}(x;m_{\gamma }^{(j)},P_{\gamma })}$$

(70)

where $m_{\gamma }^{(1)}={[1200,0,-400,0]}^T$, $m_{\gamma }^{(2)}={[600,0,700,0]}^T$, $m_{\gamma }^{(3)}={[600,0,-300,0]}^T$,$P_{\gamma }=\mathrm{diag}({[100,100,100,100]}^T)$. The survival probability is $p_{S,k}=0.99$ and the detection probability is $p_{D,k}=0.98$. The surveillance region is $[0,2000] \mathrm{m}\times [-1000,1000] \mathrm{m}$. The clutter is uniformly distributed over the surveillance region, and the average clutter number is $\lambda =10$. For the particles of birth targets, $V_{\gamma ,k}^{(i)}=\mathrm{diag}({[{10}^2,{10}^2]}^T)$, ${\upsilon }_{\gamma ,k}^{(i)}=4$, ${\alpha }_{\gamma ,k}^{(i)}=3$, ${\beta }_{\gamma ,k}^{(i)}=1$. The number of iteration is $\tau =5$, the decreasing factor is $\rho =0.98$, and the particles number is $N_{particle}=800$.

The trajectories of the targets are shown in Figure 2, where circles are the beginning positions and triangles are the ending positions.

To illustrate the performance of the proposed filter under the Gaussian measurement noise, $\eta$ is set to 0. Figure 3 and Figure 4 show the estimate of target number and OSPA distance of the three filters, respectively. From Figure 3, we can see that all the filters achieve unbiased estimate of the target number and have a comparable estimate accuracy under Gaussian measurement noise. This is because the SMC-PHD filter utilizes the true measurement noise covariance. For the NGIWG-PHD filter and the RSMC-PHD filter, the DOF and scale matrix of the Student-t distribution can be adjusted automatically and the target tracking under Gaussian measurement noise realized. Figure 4 shows the OSPA distance of the three filters. Obviously, all the filters have an approximate OSPA distance.

To verify the multi-tracking performance under heavy-tailed measurement noise, $\eta$ is set to 0.2. The remaining parameters are the same as above.

Figure 5 and Figure 6 show the estimate of target number and OSPA distance of the three filters, respectively. Under heavy-tailed measurement noise, the SMC-PHD filter will diverge because the assumed measurement noise is not matched with the true measurement noise. For the NGIWG-PHD filter and the RSMC-PHD filter, the two filters will produce a weight which is non-negligible for filter update when there are outliers in measurements. Thus, the estimate accuracy of the target number of the NGIWG-PHD filter and the RSMC-PHD filter is higher than that of the SMC-PHD filter. From Figure 6, it can be observed that the RSMC-PHD filter has the best performance on OSPA distance.

To verify the validity of the weight modification method, we compare the proposed RSMC-PHD filter with the proposed robust SMC-PHD filter without weight modification (RSMC-PHD(M1)). Figure 7 and Figure 8 show the mean absolute error (MAE) of the target number and the average OSPA (AOSPA) distance for the RSMC-PHD(M1) filter and the RSMC-PHD filter under different clutter number. The MAE of the target number and AOSPA are given by

$${\mathrm{MAE}}_{\mathrm{number}}=\frac{1}{M_{C}T}{\displaystyle \sum _{j=1}^{M_C}{\displaystyle \sum _{k=1}^T\left|N_k^j-{\widehat{N}}_k^j\right|}}$$

(71)

$$\mathrm{AOSPA}=\frac{1}{M_{C}T}{\displaystyle \sum _{j=1}^{M_C}{\displaystyle \sum _{k=1}^T{\mathrm{OSPA}}_k^j}}$$

(72)

where $M_C$ is the number of Monte Carle trials, $T$ is the simulation tracking time, $N_k^j$ and ${\widehat{N}}_k^j$ are the true and the estimated target numbers at time $k$ in the jth Monte Carle trial, ${\mathrm{OSPA}}_k^j$ is the corresponding OSPA distance.

Figure 7 illustrates that the MAE estimate of target number for the RSMC-PHD filter is lower than that of the RSMC-PHD(M1) filter. The more the clutter number, the higher is the probability that the clutter is regarded as the target-generated measurement, resulting in an overestimate of the target number at the current time. From Figure 8, it can be seen that the RSMC-PHD filter also outperforms the RSMC-PHD(M1) filter on AOSPA distance.

4.3. Nonlinear Scenario

To further evaluate the performance of the proposed filter, we consider a nonlinear scenario. For the NGIWG-PHD filter, EKF is used to handle nonlinearity.

The target state transition function is given by

$$x_k=F_{ct}{x}_{k-1}+G_{ct}{w}_k$$

(73)

where

$$F_{ct}=\left[\begin{array}{ccccc}1& \sin \omega \Delta /\omega & 0& -(1-\cos \omega \Delta )/\omega & 0\\ 0& \cos \omega \Delta & 0& -(1-\cos \omega \Delta )/\omega & 0\\ 0& (1-\cos \omega \Delta )/\omega & 1& \sin \omega \Delta /\omega & 0\\ 0& \sin \omega \Delta /\omega & 0& \cos \omega \Delta & 0\\ 0& 0& 0& 0& 1\end{array}\right]$$

(74)

$$G_{ct}={\left[\begin{array}{ccccc}{\Delta }^2/2& \Delta & 0& 0& 0\\ 0& 0& {\Delta }^2/2& \Delta & 0\\ 0& 0& 0& 0& 1\end{array}\right]}^T$$

(75)

The target state is $x_k={[x_k,{\dot{x}}_k,y_k,{\dot{y}}_k,\omega ]}^T$, $\omega$ is the turn rate, $w_k\sim \mathcal{N}(\cdot ;{\mathit{0}}_2,Q_k)$, $Q_k=\mathrm{diag}({[{\sigma }_x^2,{\sigma }_y^2,{\sigma }_{\omega }^2]}^T)$, ${\sigma }_x^{}={\sigma }_y^{}=5 { \mathrm{m}/\mathrm{s}}^2$, ${\sigma }_{\omega }^{}=5(\pi /180) \mathrm{rad}/\mathrm{s}$, $\Delta =1 \mathrm{s}$.

The observation model is formulated by

$$z_k=\left[\begin{array}{c}\sqrt{x_k^2+y_k^2}\\ \arctan (y_k/x_k)\end{array}\right]+{\epsilon }_k$$

(76)

Assume the nominal measurement noise variance is $R_{\mathit{0}}=\mathrm{diag}({[{10}^2,{\left(\pi /180 \right)}^2]}^T)$, and the true measurement noise is represented as (69) with $R_n=20R_{\mathit{0}}$ and $\eta =0.2$.

The model of birth target is a Poisson RFS with intensity

$${\gamma }_k(x)={\displaystyle {\sum }_{j=1}^{3}0.02\mathcal{N}(x;m_{\gamma }^{(j)},P_{\gamma })}$$

(77)

where $P_{\gamma }=\mathrm{diag}({[900,900,900,900,{(6\times \pi /180)}^2]}^T)$, $m_{\gamma }^{(1)}={[-400,0,600,0,0]}^T$, $m_{\gamma }^{(2)}={[-250,0,250,0,0]}^T$, and $m_{\gamma }^{(3)}={[250,0,450,0,0]}^T$. The surveillance region is $\left[0,\pi \right]\mathrm{rad}\times \left[0,1500\right]\mathrm{m}$. The clutter is uniformly distributed over the surveillance region and the average clutter number is $\lambda =10$. For the particles of birth targets,$V_{\gamma ,k}^{(i)}=\mathrm{diag}({[{10}^2,{\left(\pi /180 \right)}^2]}^T)$, ${\upsilon }_{\gamma ,k}^{(i)}=4$, ${\alpha }_{\gamma ,k}^{(i)}=3$, ${\beta }_{\gamma ,k}^{(i)}=1$. The remaining parameters are the same as that in Section 4.2.

The trajectories of the targets are shown in Figure 9, where circles are the beginning positions and triangles are the ending positions.

Figure 10 and Figure 11 show the estimate of target number and OSPA distance for the two filters, respectively. It indicates that the NGIWG-PHD filter performs worse in the nonlinear system. This is because the NGIWG-PHD filter uses EKF to linearize the nonlinear model, and the estimate error of the target state will become large. Under the influence of heavy-tailed measurement noise, the target will miss and the target number will be underestimated. On the contrary, the RSMC-PHD filter is not subject to nonlinear constraint, and has a better performance in the nonlinear system.

5. Conclusions

In this paper, a robust SMC-PHD filter with unknown heavy-tailed measurement noise is proposed. The proposed filter uses the Student-t distribution to model the unknown measurement noise. The IW distribution is used to model the unknown scale matrix parameter. Meanwhile, the Gamma distribution is used to model the unknown DOF parameter. Then, the unknown noise parameters are updated by the VB technique within the SMC-PHD filter framework. In addition, a strategy is applied to solve the overestimate of the target number induced by the Student-t distribution. Simulation results under the linear system and nonlinear system show that the proposed filter has a better performance than the compared filters with the unknown heavy-tailed measurement noise, especially in a nonlinear system.