Variational Bayesian Based Adaptive Shifted Rayleigh Filter for Bearings-Only Tracking in Clutters

Abstract

This paper considers bearings-only target tracking in clutters with uncertain clutter probability. The traditional shifted Rayleigh filter (SRF), which assumes known clutter probability, may have degraded performance in challenging scenarios. To improve the tracking performance, a variational Bayesian-based adaptive shifted Rayleigh filter (VB-SRF) is proposed in this paper. The target state and the clutter probability are jointly estimated to account for the uncertainty in clutter probability. Performance of the proposed filter is evaluated by comparing with SRF and the probability data association (PDA)-based filters in two scenarios. Simulation results show that the proposed VB-SRF algorithm outperforms the traditional SRF and PDA-based filters especially in complex adverse scenarios in terms of track continuity, track accuracy and robustness with a little higher computation complexity.

1. Introduction

Bearings-only target tracking is to estimate the current position and velocity of a target using only the noise-corrupted bearing measurements from one or multiple observer platforms. It is an important tracking problem that arises in both military and civilian applications, such as underwater sonar tracking, bistatic radar, air traffic control and computer vision [1].

Because of the intrinsic nonlinearities in the measurement models, it is difficult to acquire an optimal solution of this problem. Several suboptimal algorithms have been developed for bearings-only tracking in the literature. The extended Kalman filter (EKF) [2] in the Cartesian coordinate system is an early attempt. However, it is easy to diverge. To improve the stability of the EKF, the bearings-only tracking problem was formulated in modified polar coordinates, resulting in the modified polar coordinate EKF (MPEKF) [3]. However, it requires good initialization to guarantee convergence. The well-known pseudo-linear estimator (PLE) [4] was also developed to solve the bearings-only tracking problem. However, it gives a biased estimate at long ranges. In recent years, some bias compensation techniques were developed to improve the performance of PLE [5,6,7]. In addition, more sophisticated nonlinear Kalman filtering algorithms, such as unscented Kalman filter (UKF) [8,9], cubature Kalman filter (CKF) [1] and particle filter (PF) [10,11] were applied for bearings-only target tracking. PF can provide good performance but at the price of heavy computation load. Noteworthily, Clark et al. [12,13] proposed a novel shifted Rayleigh filter (SRF) for bearings-only tracking, which is still based on the approximation of conditional expectations but with novel feature which is performing a calculation to exploit the essential structure of the nonlinearities in a new way. It is shown to exhibit similar performance to the PF in certain challenging scenarios with much lower computational complexity [14].

However, these algorithms do not consider the impact of clutter which makes the bearings-only tracking problem more intractable. Apparently, the classical treatment of clutter in target tracking problem can be extended to the bearings-only tracking problem. For example, Reference [15] integrated the maximum entropy fuzzy probabilistic data association (MEFPDA) with the square-root cubature Kalman filter (SCKF) to deal with the clutter in bearings-only tracking. Clark et al. also included the effect of clutter measurements into the SRF algorithm in [12]. However, in the classical algorithms, the clutter probability is usually assumed known and constant, which maybe time-varying or hard to determine in advance especially in adverse scenarios. Use of incorrect value of the clutter probability may lead to track accuracy degradation even track loss. A straight-forward idea is to account for the unknown clutter probability in the process of estimation of the state. That is, we need to solve the problem of bearings-only target tracking with uncertain parameter of clutter probability.

As we know, the Bayesian approach is the most general approach of solving the problem with uncertain parameters. However, it is not trivial to get the analytical solution for most Bayesian approaches due to complex nonlinear probability density function or high dimension of integration. Recently, the variational Bayesian (VB)-based adaptive filters [16,17,18] have drawn extensive attention, which utilize a new simpler, analytically tractable distribution to approximate the true posterior distribution so that avoiding direct calculations of complex integrals. Its adaptive strategy has a strong ability of tracking time-varying parameters. Therefore, we adopt the VB method to jointly estimate the target state and the clutter probability within the framework of SRF in this paper for bearings-only tracking in clutters. By establishing a conjugate exponential model for clutter probability and data association indicator, the proposed filter is derived using the iterative filtering framework. The tracking performance of the proposed VB-SRF is evaluated by comparing with SRF, PDA-SCKF [15] and MEFPDA-SCKF [15] via two simulation examples. It shows that the proposed filter outperforms the traditional SRF and two PDA-SCKF-based filters in complex mismatched scenarios in terms of track continuity and track accuracy but at the cost of higher computation complexity.

The remainder of the paper is organized as follows. Section 2 gives the problem formulation. In Section 3, the variational Bayesian filtering is described. Section 4 derives the VB-based adaptive SRF. Section 5 provides simulation results and performance evaluation of the proposed approach, followed by the conclusions in Section 6.

2. The Shifted Rayleigh Filter Algorithm

2.1. The Bearing Model

Considering the shifted Rayleigh filter (SRF) for bearings-only tracking in $R^2$, the state equation and the measurement equation are described as [13]:

$$\begin{array}{cc}& {\mathbf{x}}_k={\mathbf{F}}_{k-1}{\mathbf{x}}_{k-1}+{\mathbf{u}}_{k-1}^s+{\mathbf{v}}_{k-1}\hfill \end{array}$$

(1)

$$\begin{array}{c}\begin{array}{cc}{\mathbf{y}}_k\hfill & ={\mathbf{H}}_k{\mathbf{x}}_k+{\mathbf{u}}_k^m+{\mathbf{w}}_k\hfill \\ {\mathbf{b}}_k\hfill & =\Pi \left({\mathbf{y}}_k\right)\hfill \end{array}\begin{array}{c}\end{array}\hfill \end{array}$$

(2)

where, ${\mathbf{x}}_k$ is the state vector which describes the position and velocity of the target; ${\mathbf{b}}_k$ is the noisy bearing measurement. $\Pi$ denotes the projection of the plane onto the unit circle. That is taking a 2-vector ${\mathbf{y}}_k$ into its normalized form ${{y}}_k/\left|\right|{\mathbf{y}}_k\left|\right|$. Then ${\mathbf{b}}_k={(sin{\theta }_k,cos{\theta }_k)}^T$, where ${\theta }_k$ is the bearing of the target position relative to the sensor platform. ${\mathbf{F}}_{k-1}$ and ${\mathbf{H}}_k$ are the state transition matrix and measurement matrix, respectively. ${\mathbf{u}}_{k-1}^s$ and ${\mathbf{u}}_k^m$ are the inputs to the system to increase the versatility of the model, for example, to reflect known perturbations to the dynamics and changes in sensor location; ${\mathbf{v}}_{k-1}$ is the Gaussian process noise with zero mean and covariance ${\mathbf{Q}}_v$, and ${\mathbf{w}}_k$ is the Gaussian measurement noise with zero mean and covariance ${\mathbf{Q}}_w$ and independent of ${\mathbf{v}}_k$.

The unusual point in the measurement model is that the noise ${\mathbf{w}}_k$ is modelled as additive noise present in an “augmented” measurement, ${\mathbf{y}}_k$, of the Cartesian coordinates of target relative to the sensor platform, which is projected onto the plane to generate the actual bearing measurement ${\mathbf{b}}_k$. It is different with the traditional “angle plus white noise” model, expressed as

$${\theta }_k=arctan(d_1/d_2)+{ϵ}_k$$

(3)

where, $d_1,d_2$ are the components of the displacement vector ${\mathbf{d}}_k={\mathbf{H}}_k{\mathbf{x}}_k+{\mathbf{u}}_k^m$. ${ϵ}_k$ is the sensor noise with Gaussian distribution $\mathcal{N}(0,{\sigma }^2)$ and independent of the displacement vector ${\mathbf{d}}_k$.

Actually, as explained in [13], the shifted Rayleigh bearing model (2) can be related with the traditional model (3) by making a variant on the shifted Rayleigh noise model

$$\begin{array}{cc}{\mathbf{y}}_k^{\prime }\hfill & ={\mathbf{d}}_k+\left|\right|{\mathbf{d}}_k\left|\right|{\mathbf{e}}^{\prime }\hfill \\ {\mathbf{b}}_k^{\prime }\hfill & =\Pi \left({\mathbf{y}}_k^{\prime }\right)\hfill \end{array}$$

(4)

where, ${\mathbf{e}}^{\prime }$ is an $\mathcal{N}(\mathbf{0},{\sigma }^2{\mathbf{I}}_{2\times 2})$ distributed random variable. The only difference with model (2) is the noise term ${\mathbf{w}}_k$ used to construct the augmented measurement ${\mathbf{y}}_k$ is replaced by $\left|\right|{\mathbf{d}}_k\left|\right|{\mathbf{e}}^{\prime }$. $\left|\right|{\mathbf{d}}_k\left|\right|{\mathbf{e}}^{\prime }$ differs from ${\mathbf{w}}_k$, but has identical first and second moments, and is uncorrelated with ${\mathbf{d}}_k$.

The angular ${\theta }_k^{\prime }$ of the modified vector bearing ${\mathbf{b}}_k^{\prime }$ can be represented as

$${\theta }_k^{\prime }=arctan(d_1/d_2)+{ϵ}_k^{\prime }$$

(5)

where ${ϵ}_k^{\prime }$ is a zero mean random variable, restricted to $[-\pi ,\pi ]$, independent of ${\mathbf{d}}_k$, with density ${\alpha }_{\sigma }(·)$:

$${\alpha }_{\sigma }\left(\theta \right)=\frac{e^{-1/2{\sigma }^2}}{2\pi }\left(1+\sqrt{2\pi }\frac{cos\theta }{\sigma }{F}_{normal}\left(\frac{cos\theta }{\sigma }\right)e^{1/2{(cos\theta /\sigma )}^2}\right)$$

(6)

where, $F_{normal}(·)$ is the cumulative distribution function of a standard $\mathcal{N}(0,1)$ variable.

Note that ${\theta }^{\prime }$ given by (5), is very close to the bearing $\theta$ in the standard model (3). The only difference is that the densities of the noise terms used in their construction are ${\alpha }_{\sigma }\left(\theta \right)$ and the normal $\mathcal{N}(0,{\sigma }^2)$ density, respectively. Reference [13] plots the two densities for ${\sigma }^2\le 1$. It shows the two density functions are virtually indistinguishable.

Given the bearing model, the SRF is to calculate the estimates of the conditional mean and covariance of the target state ${\mathbf{x}}_k$, given measurements up to time k, ${\mathbf{b}}_{1:k}$. That is,

$${\widehat{\mathbf{x}}}_{k|k}=E\left[{\mathbf{x}}_k\right|{\mathbf{b}}_{1:k}],{\mathbf{P}}_{k|k}=cov\left[{\mathbf{x}}_k\right|{\mathbf{b}}_{1:k}]$$

(7)

The formulas for the SRF algorithm can be seen in [13].

2.2. The Treatment of Clutter

Accounting for the effects of clutter on the bearing measurements, we represent a cluttered bearing measurement as

$$z_k=(1-r_k){\theta }_k+r_k{U}_k$$

(8)

where $U_k$ is the bearing measurement of the clutter, which is assumed to be an uniform random variable on $[-\pi ,\pi ]$; ${\theta }_k$ is the bearing measurement of the actual target from the sensor, $r_k$ is defined as data association indicator at time k.

$$\begin{array}{c}\hfill r_k=\left\{\begin{array}{cc}1\hfill & \mathrm{if}\mathrm{the}\mathrm{measurement}\mathrm{is}\mathrm{from}\mathrm{clutter}\hfill \\ 0\hfill & \mathrm{if}\mathrm{the}\mathrm{measurement}\mathrm{is}\mathrm{from}\mathrm{target}\hfill \end{array}\right.\end{array}$$

(9)

The prior of $r_k$ is assumed independent of the previous data associations and can be described as

$$\begin{array}{c}\hfill p\left(r_k\right)=\left\{\begin{array}{cc}\xi \hfill & \mathrm{if}{r}_k=1\hfill \\ 1-\xi \hfill & \mathrm{if}{r}_k=0\hfill \end{array}\right.\end{array}$$

(10)

Then the likelihood of the measurement is

$$\begin{array}{c}\hfill p\left(z_k\right|{\mathbf{x}}_k,r_k)=\left\{\begin{array}{cc}1/2\pi \hfill & if{r}_k=1\hfill \\ f\left({\theta }_k\right|{\mathbf{x}}_k)\hfill & if{r}_k=0\hfill \end{array}\right.\end{array}$$

(11)

It can be represented in two forms:

$$\begin{array}{}(12)& \hfill p\left(z_k\right|{\mathbf{x}}_k,r_k)& =\frac{1}{2\pi }{r}_k+f\left({\theta }_k\right|{\mathbf{x}}_k)(1-r_k)\hfill (13)& & ={\left(\frac{1}{2\pi }\right)}^{r_k}f{\left({\theta }_k\right|{\mathbf{x}}_k)}^{1-r_k}\hfill \end{array}$$

Note that the two forms are used for different purpose in the later section. $f\left({\theta }_k\right|{\mathbf{x}}_k)$ is the likelihood of the target bearing measurement, whose expression is given in Appendix A.

Here, we need to explain the reason for using this representation (8) of the cluttered measurement. The information carried in (8) is identical with the projection measurement ${\mathbf{z}}_k^b=[sin{z}_k,cos{z}_k]$. So the mean and covariance of the target state ${\mathbf{x}}_k$ conditioned on cluttered measurement ${\mathbf{z}}_k^b$, which are the aims needs to be achieved in the SRF, are equivalent with these conditioned on $z_k$. Furthermore, this representation (8) is more convenient for calculation. Thus, $z_k$ instead of ${\mathbf{z}}_k^b$ is adopted to represent the clutter measurement.

Based on the cluttered measurement model, given clutter probability $\xi$, the state estimate and covariance can be calculated as follows.

Suppose the density of ${\mathbf{x}}_{k-1}$ given measurements up to $k-1$, $p_{k-1|k-1}\left({\mathbf{x}}_{k-1}\right)$, is normal with mean ${\overline{\mathbf{x}}}_{k-1|k-1}$ and covariance ${\overline{\mathbf{P}}}_{k-1|k-1}$. The posterior density of state ${\mathbf{x}}_k$ conditioned on cluttered bearing measurements ${\mathbf{z}}_{1:k}$ is given as

$$\begin{array}{}& \hfill p\left({\mathbf{x}}_k\right|{\mathbf{z}}_{1:k})& =p\left({\mathbf{x}}_k\right|z_k,{\mathbf{x}}_{k-1}\sim p_{k-1|k-1}\left({\mathbf{x}}_{k-1}\right))\hfill (14)& & =q_k\left(0\right)p_{k|k}\left({\mathbf{x}}_k\right)+q_k\left(1\right)p_{k|k-1}\left({\mathbf{x}}_k\right)\hfill \end{array}$$

where $q_k\left(i\right)=p(r_k=i|z_k,{\mathbf{z}}_{1:k-1})$, $p_{k|k}\left({\mathbf{x}}_k\right)$ is the non-normal density of ${\mathbf{x}}_k$ conditioned on $r_k=0$ and $z_k$, or, equivalently, on ${\theta }_k$, and $p_{k|k-1}\left({\mathbf{x}}_k\right)$ is the density of ${\mathbf{x}}_k$ when there is no target measurement at time k. It is normal with mean ${\widehat{\mathbf{x}}}_{k|k-1}$ and covariance ${\mathbf{P}}_{k|k-1}$. Thus, the state estimate and covariance at time k can be obtained as

$$\begin{array}{}& \hfill {\overline{\mathbf{x}}}_{k|k}& =E\left[{\mathbf{x}}_k\right|z_k,{\mathbf{x}}_{k-1}\sim p_{k-1|k-1}\left({\mathbf{x}}_{k-1}\right)]\hfill (15)& & =q_k\left(0\right){\widehat{\mathbf{x}}}_{k|k}+q_k\left(1\right){\widehat{\mathbf{x}}}_{k|k-1}\hfill \end{array}$$

$$\begin{array}{}& \hfill {\overline{\mathbf{P}}}_{k|k}& =cov\left[{\mathbf{x}}_k\right|z_k,{\mathbf{x}}_{k-1}\sim p_{k-1|k-1}\left({\mathbf{x}}_{k-1}\right)]\hfill & & =q_k\left(0\right)({\mathbf{P}}_{k|k}+({\widehat{\mathbf{x}}}_{k|k}-{\overline{\mathbf{x}}}_{k|k}){({\widehat{\mathbf{x}}}_{k|k}-{\overline{\mathbf{x}}}_{k|k})}^T)\hfill (16)& & +q_k\left(1\right)({\mathbf{P}}_{k|k-1}+({\widehat{\mathbf{x}}}_{k|k-1}-{\overline{\mathbf{x}}}_{k|k}){({\widehat{\mathbf{x}}}_{k|k-1}-{\overline{\mathbf{x}}}_{k|k})}^T)\hfill \end{array}$$

where ${\widehat{\mathbf{x}}}_{k|k}$ and ${\mathbf{P}}_{k|k}$ are the state estimate and its covariance based on the actual target measurement. ${\widehat{\mathbf{x}}}_{k|k-1}$ and ${\mathbf{P}}_{k|k-1}$ are the predicted target state estimate and covariance. All of these can be obtained using the basic formulas of SRF.

The conditional densities $q_k\left(0\right)$ and $q_k\left(1\right)$ are given by the equations

$$\begin{array}{}(17)& \hfill q_k\left(0\right)& =c_k(1-\xi )f_k\left({\theta }_k\right|{\mathbf{z}}_{1:k-1})\hfill (18)& \hfill q_k\left(1\right)& =\frac{c_k\xi }{2\pi }\hfill \end{array}$$

where $c_k$ is the normalizing constant and $f_k\left({\theta }_k\right|{\mathbf{z}}_{1:k-1})$ is the density of the actual target bearing ${\theta }_k$ conditioned on measurements ${\mathbf{z}}_{1:k-1}$. The expression is given in Appendix B.

However, in a complex environment, the probability of clutter is time-varying or hard to determine in advance. In this case, the probability $\xi$ is unknown, so the above formulas are not applicable. In this paper, we resort the VB method to find the joint posterior density of ${\mathbf{x}}_k$ and $\xi$ so as to account for the uncertainty in clutter probability.

3. Variational Bayesian Filtering

In this section, we first review the conjugate exponential (CE) model, and then gives the variational Bayesian solution of the CE model.

3.1. Conjugate Exponential Model

Given measurements ${\mathbf{z}}_{1:k-1}$, the posterior of the system state $p\left({\mathbf{x}}_{k-1}\right|{\mathbf{z}}_{1:k-1})$ and the posterior of the parameter $p\left({\mathbf{r}}_{k-1}\right|{\mathbf{z}}_{1:k-1})$, we assume the complete-data likelihood in the exponential family:

$$p({\mathbf{x}}_k,{\mathbf{z}}_k|{\mathbf{r}}_k,{\mathbf{z}}_{1:k-1})=g\left({\mathbf{r}}_k\right)f({\mathbf{x}}_k,{\mathbf{z}}_k)e^{\varphi {\left({\mathbf{r}}_k\right)}^{T}u({\mathbf{x}}_k,{\mathbf{z}}_k)}$$

(19)

where, $\varphi \left({\mathbf{r}}_k\right)$ is the vector of natural parameters ${\mathbf{r}}_k$, u and f are known functions, and g is a normalization constant:

$$g{\left({\mathbf{r}}_k\right)}^{-1}=\int f({\mathbf{x}}_k,{\mathbf{z}}_k)e^{\varphi {\left({\mathbf{r}}_k\right)}^{T}u({\mathbf{x}}_k,{\mathbf{z}}_k)}d{\mathbf{x}}_{k}d{\mathbf{z}}_k$$

The parameter prior is conjugate to the complete-data likelihood:

$$p\left({\mathbf{r}}_k\right|{\alpha }_k^{-},{\beta }_k^{-})=h({\alpha }_k^{-},{\beta }_k^{-})g{\left({\mathbf{r}}_k\right)}^{{\beta }_k^{-}}{e}^{\varphi {\left({\mathbf{r}}_k\right)}^T{\alpha }_k^{-}}$$

(20)

where ${\alpha }_k^{-}$ and ${\beta }_k^{-}$ are hyperparameters of the prior, and h is a normalization constant. Note the prior $p\left({\mathbf{r}}_k\right|{\alpha }_k^{-},{\beta }_k^{-})$ is said to be conjugate to the likelihood $p({\mathbf{x}}_k,{\mathbf{z}}_k|{\mathbf{r}}_k)$ if and only if the posterior

$$p\left({\mathbf{r}}_k\right|{\alpha }_k,{\beta }_k)\propto p\left({\mathbf{r}}_k\right|{\alpha }_k^{-},{\beta }_k^{-})p({\mathbf{x}}_k,{\mathbf{z}}_k|{\mathbf{r}}_k)$$

is of the same parametric form as the prior. Then we call models that satisfy Equations (19) and (20) conjugate-exponential.

3.2. VB Approximation Method

Applying Bayes’ rule, we have the joint posterior of ${\mathbf{x}}_k$ and $r_k$ as

$$p({\mathbf{x}}_k,{\mathbf{r}}_k|{\mathbf{z}}_{1:k})\propto p({\mathbf{x}}_k,{\mathbf{z}}_k|{\mathbf{r}}_k,{\mathbf{z}}_{1:k})p\left({\mathbf{r}}_k\right|{\alpha }_k^{-},{\beta }_k^{-})$$

(21)

The analytic solution to (21) would be difficult to calculate. Here, we use the VB method to approximate the true posterior distribution with a product of tractable marginal posteriors [17].

$$p({\mathbf{x}}_k,{\mathbf{r}}_k|{\mathbf{z}}_{1:k})\approx Q({\mathbf{x}}_k,{\mathbf{r}}_k)=Q_x\left({\mathbf{x}}_k\right)Q_r\left({\mathbf{r}}_k\right)$$

(22)

where $Q_x\left({\mathbf{x}}_k\right)$ and $Q_r\left({\mathbf{r}}_k\right)$ are unknown approximating marginal densities of ${\mathbf{x}}_k$ and ${\mathbf{r}}_k$.

The basic idea of VB approximation is to minimize the Kullback- Leibler (KL) divergence between the approximating posterior and the true posterior:

$$\begin{array}{c}KL\left[Q_x\left({\mathbf{x}}_k\right)Q_r\left({\mathbf{r}}_k\right)\left|\right|p({\mathbf{x}}_k,{\mathbf{r}}_k|{\mathbf{z}}_{1:k})\right]=\int Q_x\left({\mathbf{x}}_k\right)Q_r\left({\mathbf{r}}_k\right)\times log\left(\frac{Q_x\left({\mathbf{x}}_k\right)Q_r\left({\mathbf{r}}_k\right)}{p({\mathbf{x}}_k,{\mathbf{r}}_k|{\mathbf{z}}_{1:k})}\right)d{\mathbf{x}}_{k}d{\mathbf{r}}_k\hfill \end{array}$$

(23)

Given the measurements ${\mathbf{z}}_{1:k}$, we can minimize the KL divergence with respect to the probability densities $Q_x\left({\mathbf{x}}_k\right)$ and $Q_r\left({\mathbf{r}}_k\right)$ in turn, while keeping the other fixed. Then, the following equations can be given as:

$$Q_x\left({\mathbf{x}}_k\right)\propto exp({\langle \ln p({\mathbf{x}}_k,{\mathbf{r}}_k,{\mathbf{z}}_k|{\mathbf{z}}_{1:k-1})\rangle }_{{\mathbf{r}}_k})$$

(24)

$$Q_r({\mathbf{r}}_k)\propto exp({\langle \ln p({\mathbf{x}}_k,{\mathbf{r}}_k,{\mathbf{z}}_k|{\mathbf{z}}_{1:k-1})\rangle }_{{\mathbf{x}}_k})$$

(25)

where ${\langle ·\rangle }_{{\mathbf{x}}_k}$ and ${\langle ·\rangle }_{{\mathbf{r}}_k}$ denote the expectations with respect to $Q_x\left({\mathbf{x}}_k\right)$ and $Q_r\left({\mathbf{r}}_k\right)$, respectively. Obviously, it is not an explicit solution since the distribution of each parameter is dependent on the other and neither distributions is known. The mechanism of VB method is to firstly give the initial values of the parameters and then use expectation-maximum (EM) algorithm to iteratively calculate $Q_x\left({\mathbf{x}}_k\right)$ and $Q_r\left({\mathbf{r}}_k\right)$ until convergence. For the above CE models, $Q_x\left({\mathbf{x}}_k\right)$ and $Q_r\left({\mathbf{r}}_k\right)$ can be obtained from the following procedure [19]:

(1)

The VB expectation step yields:

$$Q_x\left({\mathbf{x}}_k\right)\propto f({\mathbf{x}}_k,{\mathbf{z}}_k)e^{{\langle \varphi \left({\mathbf{r}}_k\right)\rangle }_{{\mathbf{r}}_k}^{T}u({\mathbf{x}}_k,{\mathbf{z}}_k)}=p\left({\mathbf{x}}_k\right|{\mathbf{z}}_k,{\langle \varphi \left({\mathbf{r}}_k\right)\rangle }_{{\mathbf{r}}_k})$$

(26)

(2)

The VB maximization step yields that $Q_r\left({\mathbf{r}}_k\right)$ is conjugate and of the form

$$Q_r\left({\mathbf{r}}_k\right)=h({\alpha }_k,{\beta }_k^{-}g{\left({\mathbf{r}}_k\right)}^{{\beta }_k}{e}^{\varphi {\left({\mathbf{r}}_k\right)}^T{\alpha }_k}$$

(27)

where, ${\alpha }_k$ and ${\beta }_k$ are the hyper-parameters, and

$$\begin{array}{}(28)& \hfill {\alpha }_k& ={\alpha }_k^{-}+{\langle u({\mathbf{x}}_k,{\mathbf{z}}_k)\rangle }_{{\mathbf{x}}_k}\hfill (29)& \hfill {\beta }_k& ={\beta }_k^{-}+n\hfill \end{array}$$

where n is the dimension of the measurement.

4. VB Based Adaptive Shifted Rayleigh Filter with Unknown Clutter Probability

Considering the system model (1) and the measurement model (8) described in Section 2, we adopt the VB method within the SRF framework to get the joint estimation of the target state and the clutter probability.

The core is to determine the posterior approximation $Q({\mathbf{x}}_k,r_k,\xi )$. Assume factorization $Q({\mathbf{x}}_k,r_k,\xi )\approx Q_x\left({\mathbf{x}}_k\right)Q(r_k,\xi )$, we can obtain $Q_x\left({\mathbf{x}}_k\right)$, $Q_r\left(r_k\right)$ and $Q_{\xi }\left({\xi }_k\right)$ at each time k through the following procedure.

(1)

Optimization of $Q_x\left({\mathbf{x}}_k\right)$ for fixed $Q(r_k,\xi )$.

First, by using the first form of the measurement likelihood (12), the complete-data likelihood is presented as

$$\begin{array}{}& \hfill p({\mathbf{x}}_k,z_k|r_k,{\mathbf{z}}_{1:k-1})& =p\left(z_k\right|{\mathbf{x}}_k,r_k)p\left({\mathbf{x}}_k\right|r_k,{\mathbf{z}}_{1:k-1})\hfill (30)& & =[\frac{1}{2\pi }{r}_k+f\left({\theta }_k\right|{\mathbf{x}}_k)(1-r_k)]\mathcal{N}({\mathbf{x}}_k;{\widehat{\mathbf{x}}}_{k|k-1},{\mathbf{P}}_{k|k-1})\hfill \end{array}$$

Then according to (24), we can get

$$\begin{array}{}& \hfill Q_x\left({\mathbf{x}}_k\right)& \propto exp\{{\langle lnp(z_k,{\mathbf{x}}_k|r_k,{\mathbf{z}}_{1:k-1})\rangle }_{r_k,\xi }\}\hfill & & =[\frac{1}{2\pi }{\langle r_k\rangle }_{r_k}+f\left({\theta }_k\right|{\mathbf{x}}_k)(1-{\langle r_k\rangle }_{r_k})]\mathcal{N}({\mathbf{x}}_k;{\widehat{\mathbf{x}}}_{k|k-1},{\mathbf{P}}_{k|k-1})\hfill (31)& & =\frac{1}{2\pi }\mathcal{N}({\mathbf{x}}_k;{\widehat{\mathbf{x}}}_{k|k-1},{\mathbf{P}}_{k|k-1}){\langle r_k\rangle }_{r_k}+f\left({\theta }_k\right|{\mathbf{x}}_k)\mathcal{N}({\mathbf{x}}_k;{\widehat{\mathbf{x}}}_{k|k-1},{\mathbf{P}}_{k|k-1})(1-{\langle r_k\rangle }_{r_k})\hfill (32)& & \approx \frac{1}{2\pi }\mathcal{N}({\mathbf{x}}_k;{\widehat{\mathbf{x}}}_{k|k-1},{\mathbf{P}}_{k|k-1}){\langle r_k\rangle }_{r_k}+\mathcal{N}({\mathbf{x}}_k;{\widehat{\mathbf{x}}}_{k|k},{\mathbf{P}}_{k|k})f\left({\theta }_k\right|{\mathbf{z}}_{1:k-1})(1-{\langle r_k\rangle }_{r_k})\hfill \end{array}$$

The approximation sign in (32) is because the following:

$$\begin{array}{}& \hfill f\left({\theta }_k\right|{\mathbf{x}}_k)\mathcal{N}({\mathbf{x}}_k;{\widehat{\mathbf{x}}}_{k|k-1},{\mathbf{P}}_{k|k-1})& =p\left({\theta }_k\right|{\mathbf{x}}_k,{\mathbf{z}}_{1:k-1})p\left({\mathbf{x}}_k\right|{\mathbf{z}}_{1:k-1})\hfill & & =p\left({\mathbf{x}}_k\right|{\theta }_k,{\mathbf{z}}_{1:k-1})p\left({\theta }_k\right|{\mathbf{z}}_{1:k-1})\hfill (33)& & \approx \mathcal{N}({\mathbf{x}}_k;{\widehat{\mathbf{x}}}_{k|k},{\mathbf{P}}_{k|k})f\left({\theta }_k\right|{\mathbf{z}}_{1:k-1})\hfill \end{array}$$

where $f\left({\theta }_k\right|{\mathbf{z}}_{1:k-1})$ is derived in Appendix B.

Comparing (32) with the posterior density (14) of SRF, the difference lies in the weights. Except for a normalization constant, the clutter probability $\xi$ used before in (14) has been replaced by ${\langle r_k\rangle }_{r_k}$ in (32), which is updated online.

(2)

Optimization of $Q(r_k,\xi )$ for fixed $Q_x\left({\mathbf{x}}_k\right)$.

We use the VB method again by factorizing $Q(r_k,\xi )\approx Q_r\left(r_k\right)Q_{\xi }\left(\xi \right)$. Assume the conjugate prior of $r_k$ is binomial distributed with parameter $\xi$. That is,

$$\begin{array}{cc}\hfill p\left(r_k\right|\xi )& ={\xi }^{r_k}{(1-\xi )}^{1-r_k}\hfill \end{array}$$

(34)

and $p\left(\xi \right)$ follows beta distribution with parameters ${\alpha }_1$ and ${\alpha }_2$.

$$p(\xi ;{\alpha }_1,{\alpha }_2)=\frac{1}{B({\alpha }_1,{\alpha }_2)}{\xi }^{{\alpha }_1-1}{(1-\xi )}^{{\alpha }_2-1}$$

(35)

where $B({\alpha }_1,{\alpha }_2)=\Gamma \left({\alpha }_1\right)\Gamma \left({\alpha }_2\right)/\Gamma ({\alpha }_1+{\alpha }_2)$.

To have the form of (19), the complete-data likelihood $p({\mathbf{x}}_k,z_k|r_k,{\mathbf{z}}_{1:k-1})$ is re-derived using the second form (13) of $p\left(z_k\right|{\mathbf{x}}_k,r_k)$ as

$$\begin{array}{}& \hfill p({\mathbf{x}}_k,z_k|r_k,{\mathbf{z}}_{1:k-1})& ={(1/2\pi )}^{r_k}{\left[f\left({\theta }_k\right|{\mathbf{x}}_k)\right]}^{1-r_k}\mathcal{N}({\mathbf{x}}_k;{\widehat{\mathbf{x}}}_{k|k-1},{\mathbf{P}}_{k|k-1})\hfill (36)& & =f\left({\theta }_k\right|{\mathbf{x}}_k)\mathcal{N}({\mathbf{x}}_k;{\widehat{\mathbf{x}}}_{k|k-1},{\mathbf{P}}_{k|k-1})exp\{r_k[ln(1/2\pi )-ln\left[f\left({\theta }_k\right|{\mathbf{x}}_k)\right]\}\hfill \end{array}$$

Rewriting the conjugate prior of $r_k$ in the form of (20), we can get

$$\begin{array}{}& \hfill p\left(r_k\right|\xi )& ={\xi }^{r_k}{(1-\xi )}^{1-r_k}\hfill (37)& & =(1-\xi )exp\left[r_{k}ln\left(\frac{\xi }{1-\xi }\right)\right]\hfill \end{array}$$

Then applying (27), $Q_r\left(r_k\right)$ can be obtained as

$$Q_r\left(r_k\right)=(1-{\eta }_k)exp\left[r_{k}ln\left(\frac{{\eta }_k}{1-{\eta }_k}\right)\right]$$

(38)

where, ${\eta }_k$ is the hyper-parameter and updated as

$$ln(\frac{{\eta }_k}{1-{\eta }_k})={\langle ln(\frac{\xi }{1-\xi })\rangle }_{\xi }+\langle {[ln(1/2\pi )-ln\left[f\left({\theta }_k\right|{\mathbf{x}}_k)\right]\rangle }_{{\mathbf{x}}_k}$$

(39)

Likewise, rewriting the conjugate prior of $\xi$ in the form of (20) , we can get

$$\begin{array}{}& p\left(\xi \right)\hfill & =\frac{1}{B({\alpha }_1,{\alpha }_2)}{\xi }^{{\alpha }_1-1}{(1-\xi )}^{{\alpha }_2-1}\hfill & & =\frac{1}{B({\alpha }_1,{\alpha }_2)}{(1-\xi )}^{{\alpha }_1+{\alpha }_2-2}{(1-\xi )}^{-({\alpha }_1-1)}{\xi }^{{\alpha }_1-1}\hfill (40)& & =\frac{1}{B({\alpha }_1,{\alpha }_2)}{(1-\xi )}^{{\alpha }_1+{\alpha }_2-2}exp[({\alpha }_1-1)ln(\frac{\xi }{1-\xi })]\hfill \end{array}$$

Then applying (27), $Q_{\xi }\left(\xi \right)$ can be obtained as

$$Q_{\xi }\left(\xi \right)=\frac{1}{B({\alpha }_1^{\prime },{\alpha }_2^{\prime })}{(1-\xi )}^{{\alpha }_1^{\prime }+{\alpha }_2^{\prime }-2}exp[({\alpha }_1^{\prime }-1)ln(\frac{\xi }{1-\xi })]$$

(41)

with hyper-parameters

$$\begin{array}{}(42)& \hfill {\alpha }_1^{\prime }& ={\alpha }_1+{\langle r_k\rangle }_{r_k}\hfill (43)& \hfill {\alpha }_2^{\prime }& ={\alpha }_2+n-{\langle r_k\rangle }_{r_k}\hfill \end{array}$$

where n is the dimension of the measurement.

According to the approximated posteriors of $Q_r\left(r_k\right)$ and $Q_{\xi }\left(\xi \right)$, ${\langle r_k\rangle }_{r_k}$ and ${\langle ln(\frac{\xi }{1-\xi })\rangle }_{\xi }$ can be obtained as:

$$\begin{array}{}(44)& \hfill {\langle r_k\rangle }_{r_k}& ={\eta }_k\hfill (45)& \hfill {\langle ln(\frac{\xi }{1-\xi })\rangle }_{\xi }& =\psi \left({\alpha }_1^{\prime }\right)-\psi \left({\alpha }_2^{\prime }\right)\hfill \end{array}$$

where $\psi (·)$ is the digamma function.

Taking expectation and covariance on the posterior $Q_x\left({\mathbf{x}}_k\right)$, the conditional mean and covariance of the target state can then be obtained. We summarize the entire filtering procedure of the VB-based SRF (VB-SRF) in Algorithm 1.

Algorithm 1 : VB-SRF.

(1) Initialization: ${\overline{\mathbf{x}}}_{0|0}$, ${\overline{\mathbf{P}}}_{0|0}$, ${\mathbf{Q}}_v$, ${\mathbf{Q}}_w$, ${\eta }_0$, ${\alpha }_{1,0}$, ${\alpha }_{2,0}$ (2) Prediction: ${\widehat{\mathbf{x}}}_{k|k-1}={\mathbf{F}}_{k-1}{\overline{\mathbf{x}}}_{k-1|k-1}+{\mathbf{u}}_{k-1}^s$ ${\mathbf{P}}_{k|k-1}={\mathbf{F}}_{k-1}{\overline{\mathbf{P}}}_{k-1|k-1}{\mathbf{F}}_{k-1}^T+{\mathbf{Q}}_v$ ${\mathbf{S}}_k={\mathbf{H}}_k{\mathbf{P}}_{k|k-1}{\mathbf{H}}_k^T+{\mathbf{Q}}_k^m$ ${\eta }_{k|k-1}=\rho {\eta }_{k-1},{\alpha }_{1,k|k-1}=\rho {\alpha }_{1,k-1},{\alpha }_{2,k|k-1}=\rho {\alpha }_{2,k-1}$ where $\rho$ is the scale factor and $0<\rho \le 1$. (3) Update: the update of VB-SRF utilizes iterate filtering framework.  (3.a) First set: ${\overline{\mathbf{x}}}_{k|k}^{\left(0\right)}={\widehat{\mathbf{x}}}_{k|k-1}$, ${\overline{\mathbf{P}}}_{k|k}^{\left(0\right)}={\mathbf{P}}_{k|k-1}$, ${\eta }_k^{\left(0\right)}={\eta }_{k|k-1}$, ${\alpha }_{1,k}^{\left(0\right)}={\alpha }_{1,k|k-1}$, ${\alpha }_{2,k}^{\left(0\right)}={\alpha }_{2,k|k-1}$  (3.b) Calculate state estimation and its covariance using SRF when the measurement is from the target: ${\mathbf{K}}_k={\mathbf{P}}_{k|k-1}{\mathbf{H}}_k^T{\mathbf{S}}_k^{-1}$ ${\epsilon }_k={\left({\mathbf{b}}_k^T{\mathbf{S}}_k^{-1}{\mathbf{b}}_k\right)}^{-1/2}{\mathbf{b}}_k^T{\mathbf{S}}_k^{-1}({\mathbf{H}}_k{\widehat{\mathbf{X}}}_{k|k-1}+{\mathbf{u}}_k^m)$ ${\gamma }_k={\left({\mathbf{b}}_k^T{\mathbf{S}}_k^{-1}{\mathbf{b}}_k\right)}^{-1/2}{\rho }_n\left({\epsilon }_k\right)$ ${\delta }_k={\left({\mathbf{b}}_k^T{\mathbf{S}}_k^{-1}{\mathbf{b}}_k\right)}^{-1/2}[2+{\epsilon }_k{\rho }_2\left({\epsilon }_k\right)-{\rho }_2^2{\epsilon }_k]$ ${\rho }_2\left({\epsilon }_k\right)=\frac{{\epsilon }_k{e}^{-{\epsilon }_k^2/2}+\sqrt{2\pi }({\epsilon }_k^2+1)F_{normal}\left({\epsilon }_k\right)}{e^{-{\epsilon }_k^2/2}+\sqrt{2\pi }\left({\epsilon }_k\right)F_{normal}\left({\epsilon }_k\right)}$ ${\widehat{\mathbf{x}}}_{k|k}=(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k){\widehat{\mathbf{x}}}_{k|k-1}-{\mathbf{K}}_k{\mathbf{u}}_k^m+{\gamma }_k{\mathbf{K}}_k{\mathbf{b}}_k$ ${\mathbf{P}}_{k|k}=(\mathbf{I}-{\mathbf{K}}_k{\mathbf{H}}_k){\mathbf{P}}_{k|k-1}+{\delta }_k{\mathbf{K}}_k{\mathbf{b}}_k{\mathbf{b}}_k^T{\mathbf{K}}_k^T$  (3.c) For $j=1:N$, iterate the following N (N denotes iterated times) steps: $•$ Calculate the fused state estimation and its covariance: ${\overline{\mathbf{x}}}_{k|k}^{\left(j\right)}=\frac{1}{2\pi c}{\eta }_k^{(j-1)}{\widehat{\mathbf{x}}}_{k|k-1}+\frac{1}{c}(1-{\eta }_k^{(j-1)})f\left({\theta }_k\right|{\mathbf{z}}_{1:k-1}){\widehat{\mathbf{x}}}_{k|k}$ ${\overline{\mathbf{P}}}_{k|k}^{\left(j\right)}=\frac{1}{2\pi c}{\eta }_k^{(j-1)}({\mathbf{P}}_{k|k}+({\widehat{\mathbf{x}}}_{k|k}-{\overline{\mathbf{x}}}_{k|k}){({\widehat{\mathbf{x}}}_{k|k}-{\overline{\mathbf{x}}}_{k|k})}^T)$ $+\frac{1}{c}(1-{\eta }_k^{(j-1)})f\left({\theta }_k\right|{\mathbf{z}}_{1:k-1})({\mathbf{P}}_{k|k-1}+({\widehat{\mathbf{x}}}_{k|k-1}-{\overline{\mathbf{x}}}_{k|k}^{\left(j\right)}){({\widehat{\mathbf{x}}}_{k|k-1}-{\overline{\mathbf{x}}}_{k|k}^{\left(j\right)})}^T)$   where $c=\frac{1}{2\pi }{\eta }_k^{(j-1)}+f\left({\theta }_k\right)(1-{\eta }_k^{(j-1)})$ is a normalization term, and $f\left({\theta }_k\right|{\mathbf{z}}_{1:k-1})$ can be obtained using (A6). $•$ Update parameters: $ln(\frac{{\eta }_k^{\left(j\right)}}{1-{\eta }_k^{(j)}})=\psi ({\alpha }_{1,k}^{(j-1)})-\psi ({\alpha }_{2,k}^{(j-1)})+ln(1/2\pi )-lnf({\theta }_k|{\overline{\mathbf{x}}}_{k|k}^{\left(j\right)})$ ${\alpha }_{1,k}^{\left(j\right)}={\alpha }_{1,k}^{(j-1)}+{\eta }_k^{\left(j\right)}$ ${\alpha }_{2,k}^{\left(j\right)}={\alpha }_{2,k}^{(j-1)}-{\eta }_k^{\left(j\right)}+1$ $•$ End for and set ${\overline{\mathbf{x}}}_{k|k}={\overline{\mathbf{x}}}_{k|k}^{\left(N\right)}$, ${\overline{\mathbf{P}}}_{k|k}={\overline{\mathbf{P}}}_{k|k}^{\left(N\right)}$, ${\eta }_k={\eta }_k^{\left(N\right)}$, ${\alpha }_{1,k}={\alpha }_{1,k}^{\left(N\right)}$, ${\alpha }_{2,k}={\alpha }_{2,k}^{\left(N\right)}$.

5. Simulation Results

To evaluate the performance of the VB-SRF algorithm, two scenarios which are almost the same with these in [12,13] are utilized. The differences lie in the clutter probability in scenario 1 and the sensor tracks in scenario 2 which were not detailed in [12]. The two scenarios are very representative. In scenario 1, a maneuvering sensor is used in order to satisfy the condition of observability in bearings-only tracking. In scenario 2, multiple distributed sensors with large noise variance are utilized, which make the problem more challenging. The tracking performance of the VB-SRF algorithm was compared with the SRF algorithm, the MEFPDA-SCKF algorithm and PDA-SCKF algorithm in terms of track loss, track accuracy and computation complexity. The track loss is declared when the track error is large enough that making the filter diverge. Root mean square (RMS) error is used to show the track accuracy. In addition, the computation complexity is reflected by the computation time of each filter. The simulation codes can be downloaded through Github [20].

5.1. Scenario 1

In scenario 1, a target moves along a horizontal track, with zero vertical displacement, according to a white noise acceleration model. The state of the target is represented as ${\mathbf{x}}_k=[x_{1,k},{\dot{x}}_{1,k}]$, where $x_{1,k}$ and ${\dot{x}}_{1,k}$ are the horizontal distance and velocity at time k, respectively. The observer platform follows an approximately parallel track at a constant average speed. The horizontal and vertical displacements of the platform $x_k^p={[x_{1,k}^p,x_{2,k}^p]}^T$ are governed by the equations:

$$\begin{array}{}(46)& \hfill x_{1,k}^p& =4k+{\tilde{x}}_{1,k}^p\hfill (47)& \hfill x_{2,k}^p& =20+{\tilde{x}}_{2,k}^p\hfill \end{array}$$

in which ${\tilde{x}}_{1,k}^p$ and ${\tilde{x}}_{2,k}^p$ are zero mean Gaussian white noise processes, both with variance $q=1$.

The measurement is the angle (in radians) of the line-of-sight of this target from the platform. The sensor noise is Gaussian white noise with variance ${\sigma }^2={\left(0.05\right)}^2{\mathrm{rad}}^2={2.86}^2{\mathrm{deg}}^2$. The true clutter probability is 0.8. Other parameters are detailed in [13]. The configuration of the observer platform and target is illustrated in Figure 1. The bearing measurement of the target is presented in Figure 2. It can be seen that there is no abrupt change.

The clutter probability used in SRF is set as $p_c=0.8$, which is the same with the true clutter probability. The initial values of VB-SRF parameters are ${\eta }_0=0.8$, ${\alpha }_{1,0}=2$, ${\alpha }_{2,0}=10$. The clutter density $\lambda$ is calculated using $-kln(1-p_c)/2\pi$ in MEFPDA-SCKF and PDA-SCKF. Figure 3 presents the RMS target position errors of the four filters using 1000 Monte Carlo runs. It can be seen that SRF and VB-SRF have comparable performance under the correct clutter probability. MEFPDA-SCKF and PDA-SCKF have a little better tracking accuracy than SRF and VB-SRF.

In challenging scenarios, the clutter probability maybe unknown or time varying. The pre-set parameters are probably inaccurate. So here we set mistuned clutter probabilities for the four filters: (1) $p_c={\eta }_0=0.7$; (2) $p_c={\eta }_0=0.5$; (3) $p_c={\eta }_0=0.3$. The percentages of track losses in 1000 Monte Carlo runs are given in

Table 1. The percentages of track losses of four filters in two scenarios.

	Scenario 1			Scenario 2
	${\mathit{p}}_{\mathit{c}}\mathbf{=}\mathbf{0.7}$	${\mathit{p}}_{\mathit{c}}\mathbf{=}\mathbf{0.5}$	${\mathit{p}}_{\mathit{c}}\mathbf{=}\mathbf{0.3}$	${\mathit{p}}_{\mathit{c}}\mathbf{=}\mathbf{0.667}$	${\mathit{p}}_{\mathit{c}}\mathbf{=}\mathbf{0.5}$	${\mathit{p}}_{\mathit{c}}\mathbf{=}\mathbf{0.3}$
VB-SRF	0	0	0.1%	0	0	0
SRF	0.9%	1.6%	2.9%	0	0	2.7%
MEFPDA-SCKF	0	0	0	13.5%	13.3%	14.2%
PDA-SCKF	0	0	0	20.1%	20.8%	22.4%

. Clearly, the VB-SRF algorithm outperforms the SRF algorithm due to fewer track losses. Moreover the proportion of track losses increases as the mistuning aggravates. Especially, there are 2.9% tracks are lost for SRF while only 0.1% tracks are lost for VB-SRF in the worst case. In addition, we can see that MEFPDA-SCKF and PDA-SCKF have no track loss in all the three mistuned cases.

The RMS position errors of the four filters with different mistuned clutter probabilities are shown in Figure 4. We only consider the runs without track loss. From the figure, we can see that the tracking accuracy of VB-SRF is slightly better than SRF when ${\eta }_0\ge 0.5$. However, the performance difference is not obvious. When $p_c={\eta }_0=0.3$, VB-SRF exhibits distinct superiority over SRF. It implies that VB-SRF is more robust than SRF. Meanwhile, MEFPDA-SCKF and PDA-SCKF show better tracking accuracy than SRF and VB-SRF under all the three mistuned cases. They are almost not affected by the mistuning. It shows MEFPDA-SCKF and PDA-SCKF are more accurate and robust than SRF and VB-SRF under this simple scenario.

Table 2. Computation time of the four filters with 100 Monte Carlo runs for two scenarios.

	Scenario 1	Scenario 2
VB-SRF	0.7406 s	1.0236 s
SRF	0.3690 s	0.5779 s
MEFPDA-SCKF	0.2066 s	0.3314 s
MEFPDA-SCKF	0.2092 s	0.3128 s

shows the computation time of the four filters with 100 Monte Carlo runs. It is clear that VB-SRF has the maximum time of computation, which is about 2 times of that of SRF. The computation time of MEFPDA-SCKF and PDA-SCKF are comparable and both smaller than SRF and VB-SRF.

On the whole, the VB-SRF algorithm outperforms the SRF in terms of track continuity, track accuracy and robustness especially in severely mismatched scenarios but with higher computation complexity. MEFPDA-SCKF and PDA-SCKF perform better than SRF and VB-SRF in all aspects. It illustrates that the PDA-SCKF-based strategy has superiority over the SRF-based strategy in handing the clutters in simple scenarios.

5.2. Scenario 2

For scenario 2, the aim is to track a single target from several drifting sonobuoys. A monitoring aircraft estimates the positions of the drifting sonobuoys by observing the direction of arrival of sensor transmissions. The sonobuoys track the position of the target by means of noisy bearings measurements.

The state is 12-dimensional:

$${\mathbf{x}}_k={[x_{0,k},{\dot{x}}_{0,k},y_{0,k},{\dot{y}}_{0,k},x_{1,k},y_{1,k},x_{2,k},y_{2,k},x_{3,k},y_{3,k},u_{1,k},u_{1,k}]}^T$$

(48)

the first four components represent the $(x,y)$ coordinates of the position and velocity of the target, the next six, the coordinates of the positions of the three sonobuoys, and the last two, those of the drift current effecting all three sonobuoys.

Six simultaneous measurements are made at each time step. Three of these are measurements of the bearing angles of the sonobuoys from the monitoring platform and they are uncluttered. Three are the bearing angles of the target from the sonobuoys, which are subject to clutter. The standard deviation of monitoring sensor bearing noise and sonobuoy sensor bearing noise are 0.8° and 16°, respectively. The true probability of clutter is set as 0.667. In addition, the bearing of the clutter is uniformly distributed over $[-\pi ,\pi ]$. Other simulation parameters can be referred to [12]. 200 Monte Carlo runs are performed to evaluate the performance of the proposed filter. Figure 5 shows the behaviour of the estimates of target and sonobuoy positions provided by both the SRF and VB-SRF, for a typical simulation.

The RMS target position errors of the four filters with correct clutter probability assumption are given in Figure 6. Compared with scenario 1, the differences between the four filters are more dramatic in scenario 2. This is probably because scenario 2 is more complex in which multiple sensors are used to observe the target and give abruptly changing and severely noise-corrupted bearing measurements of the target, shown in Figure 7. Thus, even a minor change in filtering strategy could result in large variations in performance. Meanwhile, seen from Figure 7, an abrupt change (almost from +180° to −180°) occurs in the target bearing measurement from sonobuoy sensor 3 at k = 62 s, which leads to several track losses shown in Table 1 and much larger position errors of MEFPDA-SCKF and PDA-SCKF. Whereas, SRF and VB-SRF are less affected by the abrupt bearing variation since the value of the projected measurement $b_k={(sin{\theta }_k,cos{\theta }_k)}^T$ is invariant when there is a 360° change in bearing ${\theta }_k$. They have comparable tracking accuracy. It is hard to decide which one is better.

To compare the performance under adverse scenarios, we set mistuned clutter probabilities as: (1) $p_c={\eta }_0=0.5$; (2) $p_c={\eta }_0=0.3$. In case (1), there is no track loss in SRF and VB-SRF and the RMS position errors of the two filters are shown in Figure 8a. We can see that the RMS position errors of SRF are slightly increased compared with the case with no mistuning, while the RMS position errors of VB-SRF remain almost unchanged. In case (2), as shown in Table 1, 2.7% tracks are lost for SRF while no track is lost for VB-SRF. Meanwhile, as can be seen in Figure 8b, VB-SRF has much smaller RMS errors than the SRF. For MEFPDA-SCKF and PDA-SCKF, unlike with scenario 1, they have higher percentage of track losses and larger RMS position errors than VB-SRF in both mistuned cases. It shows that VB-SRF is superior to PDA-SCKF-based algorithms in more challenging scenarios. In addition, from Table 2, we can see that the computation time of VB-SRF is twice of the SRF and three times of MEFPDA-SCKF and PDA-SCKF. Overall, all these reveal that the proposed VB-SRF algorithm has significant performance superiority in severely mismatched and complex cases at the cost of a little higher computation complexity.

6. Conclusions

Bearings-only target tracking in the presence of clutter is a difficult problem because of the nonlinearity of the measurement model, the measurement origin uncertainty and the observability of the target. The Shifted Rayleigh filter (SRF) is shown to exhibit good performance for bearings-only target tracking in certain challenging scenarios through exploiting the essential structure of the nonlinearities in a new way. However, the clutter probability is assumed known and constant in SRF, which may not match with the truth especially in adverse scenarios. Therefore, to handle the bearings-only target tracking in clutters with uncertain clutter probability, a variational Bayesian-based adaptive shifted Rayleigh filter (VB-SRF) is proposed in this paper. By establishing a conjugate exponential model of the clutter probability and the data association indicator, the approximated posterior probability densities of the target state and the clutter parameters are iteratively calculated using the VB expectation and maximization steps. Finally, joint estimation of the target state and the clutter probability are achieved in the framework of SRF. The tracking performance of the proposed filter is compared with SRF, PDA-SCKF and MEFPDA-SCKF via two simulation examples. It shows that the proposed filter outperforms the other three filters in terms of track continuity and track accuracy with a little higher computation complexity in complex adverse scenarios. In addition, it also reveals that the proposed VB-SRF exhibits better robustness than the traditional SRF.