Tracking the Turn Maneuvering Target Using the Multi-Target Bayes Filter with an Adaptive Estimation of Turn Rate

Abstract

Tracking the target that maneuvers at a variable turn rate is a challenging problem. The traditional solution for this problem is the use of the switching multiple models technique, which includes several dynamic models with different turn rates for matching the motion mode of the target at each point in time. However, the actual motion mode of a target at any time may be different from all of the dynamic models, because these models are usually limited. To address this problem, we establish a formula for estimating the turn rate of a maneuvering target. By applying the estimation method of the turn rate to the multi-target Bayes (MB) filter, we develop a MB filter with an adaptive estimation of the turn rate, in order to track multiple maneuvering targets. Simulation results indicate that the MB filter with an adaptive estimation of the turn rate, is better than the existing filter at tracking the target that maneuvers at a variable turn rate.

1. Introduction

Target tracking has been discussed in many articles due to its military and civil applications, which range from threat warnings, to intelligent surveillance and situational awareness [1,2,3,4,5,6,7,8,9]. Maneuvering target tracking is the most essential ingredient of target tracking and has attracted the attention of many researchers. A number of efficient tracking algorithms for maneuvering targets have been developed and designed in the past few decades [10,11,12,13,14,15,16,17,18]. The interacting multiple model (IMM) algorithms were independently developed in [10,11], in order to track the maneuvering target in systems with Markov-switching coefficients, and in air traffic control, respectively. The mode-set adaptive IMM algorithm and multiple model method with variable structure were designed in [12,13], to improve the performance of IMM algorithms when tracking a maneuvering target. By combining the IMM method with the joint probabilistic data association (JPDA) and multiple hypothesis tracking (MHT) techniques, respectively, the IMM-JPDA filter and the IMM-MHT filter were developed in [14,15], in order to track multiple maneuvering targets. By applying the switching multiple models technique to the probability hypothesis density (PHD) filter and multi-target Bayes (MB) filter, respectively, Pasha developed a PHD filter to track maneuvering targets in the presence of clutter and noise [16], and Liu designed a MB filter for multiple maneuvering target tracking, in the case of low detection probability [17].

As mentioned above, the existing methods for maneuvering target tracking apply the IMM approach, or the switching multiple models technique, to the tracking filter. In these methods, a finite set of dynamic models are used each time. Because the motion mode space of a target is continuous, a sufficiently large set of dynamic models is usually required to cover the range of possible motion modes of the target. Such a large set is impractical because an increase in the number of dynamic models, also leads to an increase in the computational load. Additionally, it is worth noting that the actual motion mode of a target at any given time, may be different from all of the dynamic models, even if a sufficiently large set of dynamic models are used by the filter.

When tracking the target that maneuvers at a variable turn rate, a limited set of dynamic models with different turn rates are usually used by the filter [12,16]. Since the turn rate of a target at any given time is unknown and random, its actual turn rate at any given time may be different from the turn rates in dynamic models. This difference causes the filter to provide an inaccurate state estimation of the target at this time.

To track the target that maneuvers at a random turn rate, we establish a formula for computing the turn rate of the target. This formula solves the estimation issue of the turn rate by using the state vector of a target at a previous time, and its measurement at the current time. Applying the estimation method of the turn rate to the MB filter, we present the MB filter with an adaptive estimation of the turn rate. Its performance is demonstrated by the simulation results.

2. A Brief Description of Pasha’s PHD Filter

Pasha’s Gaussian mixture PHD filter is applied to the tracking of multiple maneuvering targets in systems with linear Gaussian jump Markov system models, and is used as the comparison object in this paper. We will first give a brief description of this filter. A simplified version of Pasha’s Gaussian mixture PHD filter is composed of the following four steps:

Step 1: Prediction

Let $v_{k-1}({\mathit{x}}_{k-1},r_{k-1})={\displaystyle \sum _{i=1}^{N_{k-1}}{w}_{i,k-1}}(r_{i,k-1})N({\mathit{x}}_{i,k-1};{\mathit{m}}_{i,k-1}(r_{i,k-1}),{\mathit{P}}_{i,k-1}(r_{i,k-1}))$ denote the posterior intensity at time $k-1$, where $N_{k-1}$ is the number of Gaussian terms at time $k-1$ and $N(\cdot ;\mathit{m},\mathit{P})$ is a Gaussian distribution with mean vector $\mathit{m}$, and covariance matrix $\mathit{P}$; $r_{i,k-1}$, ${\mathit{x}}_{i,k-1}$, $w_{i,k-1}(r_{i,k-1})$, ${\mathit{m}}_{i,k-1}(r_{i,k-1})$, and ${\mathit{P}}_{i,k-1}(r_{i,k-1})$ are the model label, state vector, weight, mean vector, and covariance matrix of Gaussian term $i$, respectively. The predicted posterior intensity is given by:

$$v_{k|k-1}({\mathit{x}}_k,r_k)={\displaystyle \sum _{i=1}^{N_{k-1}}{\displaystyle \sum _{r_{i,k}=1}^{M_r}{w}_{i,k|k-1}}}(r_{i,k})N({\mathit{x}}_{i,k};{\mathit{m}}_{i,k|k-1}(r_{i,k}),{\mathit{P}}_{i,k|k-1}(r_{i,k}))$$

(1)

where $M_r$ is the number of models used, and $w_{i,k|k-1}(r_{i,k})$, ${\mathit{m}}_{i,k|k-1}(r_{i,k})$, and ${\mathit{P}}_{i,k|k-1}(r_{i,k})$ are given by:

$$w_{i,k|k-1}(r_{i,k})=p_{S,k}{t}_{k|k-1}(r_{i,k}|r_{i,k-1})w_{i,k-1}(r_{i,k-1})$$

(2)

$${\mathit{m}}_{i,k|k-1}(r_{i,k})={\mathit{F}}_{k-1}(r_{i,k}){\mathit{m}}_{i,k-1}(r_{i,k-1})$$

(3)

$${\mathit{P}}_{i,k|k-1}(r_{i,k})={\mathit{Q}}_{k-1}(r_{i,k})+{\mathit{F}}_{k-1}(r_{i,k}){\mathit{P}}_{i,k-1}(r_{i,k-1}){\mathit{F}}_{k-1}^{\mathrm{T}}(r_{i,k})$$

(4)

where $p_{S,k}$ is the survival probability; ${\mathit{F}}_{k-1}(r_{i,k})$ and ${\mathit{Q}}_{k-1}(r_{i,k})$ are the state transition and process noise covariance matrices, respectively, of model $r_{i,k}$; and $t_{k|k-1}(r_{i,k}|r_{i,k-1})$ is the Markov transition probability from model $r_{i,k-1}$ to model $r_{i,k}$.

Step 2: Update

If the predicted posterior intensity is given by Equation (1), then the updated posterior intensity is given by:

$$\begin{array}{cc}{v}_{k|k}({\mathit{x}}_k,r_k)=& {\displaystyle \sum _{i=1}^{N_{k-1}}{\displaystyle \sum _{r_{i,k}=1}^{M_r}{w}_{i,k|k}}}(r_{i,k})N({\mathit{x}}_{i,k};{\mathit{m}}_{i,k|k-1}(r_{i,k}),{\mathit{P}}_{i,k|k-1}(r_{i,k}))\\ & +{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \sum _{i=1}^{N_{k-1}}{\displaystyle \sum _{r_{i,k}}^{M_r}{w}_{i,k|k}^j}}}(r_{i,k})N({\mathit{x}}_{i,k};{\mathit{m}}_{i,k|k}^j(r_{i,k}),{\mathit{P}}_{i,k|k}^j(r_{i,k}))\end{array}$$

(5)

where:

$$w_{i,k|k}(r_{i,k})=(1-p_{D,k})w_{i,k|k-1}(r_{i,k})$$

(6)

$${\mathit{m}}_{i,k|k}^j(r_{i,k})={\mathit{m}}_{i,k|k-1}(r_{i,k})+{\mathit{A}}_{i,k}(r_{i,k})\cdot ({\mathit{y}}_{j,k}-\mathit{H}(r_{i,k}){\mathit{m}}_{i,k|k-1}(r_{i,k}))$$

(7)

$${\mathit{P}}_{i,k|k}^j(r_{i,k})=(\mathit{I}-{\mathit{A}}_{i,k}(r_{i,k})\cdot \mathit{H}(r_{i,k})){\mathit{P}}_{i,k|k-1}(r_{i,k})$$

(8)

$${\mathit{A}}_{i,k}(r_{i,k})={\mathit{P}}_{i,k|k-1}(r_{i,k}){\mathit{H}}^{\mathrm{T}}(r_{i,k}){[\mathit{H}(r_{i,k}){\mathit{P}}_{i,k|k-1}(r_{i,k}){\mathit{H}}^{\mathrm{T}}(r_{i,k})+\mathit{R}(r_{i,k})]}^{-1}$$

(9)

$$w_{i,k|k}^j(r_{i,k})=\frac{p_{D,k}{w}_{i,k|k-1}(r_{i,k})N({\mathit{y}}_{j,k};\mathit{H}(r_{i,k}){\mathit{m}}_{i,k|k-1}(r_{i,k}),\mathit{H}(r_{i,k}){\mathit{P}}_{i,k|k-1}(r_{i,k}){\mathit{H}}^{\mathrm{T}}(r_{i,k})+\mathit{R}(r_{i,k}))}{{\lambda }_{c,k}+p_{D,k}{\displaystyle \sum _{e=1}^{N_{k-1}}{\displaystyle \sum _{r_{e,k}=1}^{M_r}{w}_{e,k|k-1}(r_{e,k})N({\mathit{y}}_{j,k};\mathit{H}(r_{e,k}){\mathit{m}}_{e,k|k-1}(r_{e,k}),\mathit{H}(r_{e,k}){\mathit{P}}_{e,k|k-1}(r_{i,k}){\mathit{H}}^{\mathrm{T}}(r_{e,k})+\mathit{R}(r_{e,k}))}}}$$

(10)

where $M_r$ is the number of measurements at time $k$, ${\mathit{y}}_{j,k}$ denotes a measurement at time $k$; $\mathit{H}(r_{i,k})$ and $\mathit{R}(r_{i,k})$ are the observation matrix and covariance matrix of observation noise, respectively; and $\mathit{I}$, ${\lambda }_{c,k}$, and $p_{D,k}$ denote the identity matrix, clutter rate, and detection probability, respectively.

Step 3: Generation of the Birth Intensity

The birth intensity is generated from the measurements at time $k$ and is given by:

$$\gamma (x_k,r_k)={\displaystyle \sum _{j=1}^{M_k}{w}_{\gamma ,k}^j}N({\mathit{x}}_{j,k};{\mathit{m}}_{\gamma ,k}^j(r_{j,k}),{\mathit{P}}_{\gamma ,k}^j(r_{j,k}))$$

(11)

where ${\mathit{m}}_{\gamma ,k}^j$ is taken from measurement ${\mathit{y}}_{j,k}={\left[\begin{array}{cc}{r}_{x,k}^j& r_{y,k}^j\end{array}\right]}^{\mathrm{T}}$; and ${\mathit{m}}_{\gamma ,k}^j={\left[\begin{array}{cccc}{r}_{x,k}^j& 0& r_{y,k}^j& 0\end{array}\right]}^{\mathrm{T}}$, $r_{j,k}=1$, $w_{\gamma ,k}^j={\rho }_r$, and ${\mathit{P}}_{\gamma ,k}^j={\mathit{P}}_{\gamma }$ where ${\rho }_r$ and ${\mathit{P}}_{\gamma }$ are the known parameter and covariance matrix, respectively.

Step 4: Combination of the Updated Posterior Intensity and Birth Intensity

The posterior intensity at time $k$ is obtained by the combination of the updated posterior intensity in Equation (5), and the birth intensity in Equation (11), which is given by:

$$\begin{array}{cc}\hfill v_k({\mathit{x}}_k,r_k)& =v_{k|k}({\mathit{x}}_k,r_k)+\gamma ({\mathit{x}}_k,r_k)\hfill \\ & ={\displaystyle \sum _{i=1}^{N_{k-1}}{\displaystyle \sum _{r_{i,k}=1}^{M_r}{w}_{i,k|k}}}(r_{i,k})N({\mathit{x}}_{i,k};{\mathit{m}}_{i,k|k-1}(r_{i,k}),{\mathit{P}}_{i,k|k-1}(r_{i,k}))\hfill \\ & +{\displaystyle \sum _{j=1}^{M_k}{\displaystyle \sum _{i=1}^{N_{k-1}}{\displaystyle \sum _{r_{i,k}}^{M_r}{w}_{i,k|k}^j}}}(r_{i,k})N({\mathit{x}}_{i,k};{\mathit{m}}_{i,k|k}^j(r_{i,k}),{\mathit{P}}_{i,k|k}^j(r_{i,k}))+{\displaystyle \sum _{j=1}^{M_k}{w}_{\gamma ,k}^j}N({\mathit{x}}_{j,k};{\mathit{m}}_{\gamma ,k}^j,{\mathit{P}}_{\gamma ,k}^j)\hfill \end{array}$$

(12)

After this combination, the Gaussian terms whose weight is less than threshold $\tau$ are pruned, and the posterior intensity, which is composed of the remaining Gaussian terms, is propagated to the next time step. Those Gaussian terms whose weight is greater than 0.5 are picked as the output of the filter at time $k$.

3. Estimation of Turn Rate

In this Section, we will estimate the turn rate of a maneuvering target by using its state vector at time $k-1$ and its position measurement at time $k$. Figure 1 shows a maneuvering target with turn rate ${\omega }_k$ which moves from point $O^e$ at time $k-1$, to point $E$ at time $k$. Let ${\mathit{X}}_{k-1}={\left[\begin{array}{cccc}{x}_{k-1}& {\dot{x}}_{k-1}& y_{k-1}& {\dot{y}}_{k-1}\end{array}\right]}^{\mathrm{T}}$ and ${\mathit{X}}_k={\left[\begin{array}{cccc}{x}_k& {\dot{x}}_k& y_k& {\dot{y}}_k\end{array}\right]}^{\mathrm{T}}$ denote the state vectors of the target at times $k-1$ and $k$, respectively, where $(x_{k-1},y_{k-1})$ and $(x_k,y_k)$ denote the position coordinates of the target at times $k-1$ and $k$, respectively, and $({\dot{x}}_{k-1},{\dot{y}}_{k-1})$ and $({\dot{x}}_k,{\dot{y}}_k)$ denote its velocities at times $k-1$ and $k$, respectively. Obviously, in the x-y Cartesian coordinate system, the Cartesian coordinates of points $O^e$ and $E$ are $(x_{k-1},y_{k-1})$ and $(x_k,y_k)$, respectively.

We introduce an ${\mathrm{x}}^{\mathrm{e}}$ − ${\mathrm{y}}^{\mathrm{e}}$ Cartesian coordinate system, whose origin is located at point $O^e$. The introduced ${\mathrm{x}}^{\mathrm{e}}$ − ${\mathrm{y}}^{\mathrm{e}}$ Cartesian coordinate system is shown in Figure 2. The coordinate transformation from position coordinates $\mathit{X}$ in the x-y coordinate system, to position coordinates ${\mathit{X}}^c$ in the ${\mathrm{x}}^{\mathrm{e}}$ − ${\mathrm{y}}^{\mathrm{e}}$ coordinate system, is given by:

$${\mathit{X}}^c=\left[\begin{array}{cc}\cos {\alpha }_{k-1}& \sin {\alpha }_{k-1}\\ -\sin {\alpha }_{k-1}& \cos {\alpha }_{k-1}\end{array}\right]\times \left(\mathit{X}-\left[\begin{array}{c}{x}_{k-1}\\ y_{k-1}\end{array}\right]\right)$$

(13)

where ${\alpha }_{k-1}$ is given by:

$${\alpha }_{k-1}=arc\cos \frac{{\dot{x}}_{k-1}}{\sqrt{{\dot{x}}_{k-1}^2+{\dot{y}}_{k-1}^2}}$$

(14)

Similarly, the vector transformation from velocity vector $\mathit{V}$ in the x − y coordinate system, to velocity vector ${\mathit{V}}^c$ in the ${\mathrm{x}}^{\mathrm{e}}$ − ${\mathrm{y}}^{\mathrm{e}}$ coordinate system, is given by:

$${\mathit{V}}^c=\left[\begin{array}{cc}\cos {\alpha }_{k-1}& \sin {\alpha }_{k-1}\\ -\sin {\alpha }_{k-1}& \cos {\alpha }_{k-1}\end{array}\right]\times \mathit{V}$$

(15)

Using Equations (13) and (15) to address position coordinates $\mathit{X}={\left[\begin{array}{cc}{x}_{k-1}& y_{k-1}\end{array}\right]}^{\mathrm{T}}$ and velocity vector $\mathit{V}={\left[\begin{array}{cc}{\dot{x}}_{k-1}& {\dot{y}}_{k-1}\end{array}\right]}^{\mathrm{T}}$ in state vector ${\mathit{X}}_{k-1}={\left[\begin{array}{cccc}{x}_{k-1}& {\dot{x}}_{k-1}& y_{k-1}& {\dot{y}}_{k-1}\end{array}\right]}^{\mathrm{T}}$, respectively, we obtain the state vector of the target in the ${\mathrm{x}}^{\mathrm{e}}$ − ${\mathrm{y}}^{\mathrm{e}}$ Cartesian coordinate system, as:

$${\mathit{X}}_{k-1}^e=\left[\begin{array}{c}{x}_{k-1}^e\\ {\dot{x}}_{k-1}^e\\ y_{k-1}^e\\ {\dot{y}}_{k-1}^e\end{array}\right]=\left[\begin{array}{c}0\\ \sqrt{{\dot{x}}_{k-1}^2+{\dot{y}}_{k-1}^2}\\ 0\\ 0\end{array}\right]$$

(16)

Since the target moves at turn rate ${\omega }_k$, from time $k-1$ to time $k$, the state transition matrix of the target motion is given by:

$$\mathit{F}({\omega }_k)=\left[\begin{array}{cccc}1& \frac{\sin {\omega }_k\Delta t_k}{{\omega }_k}& 0& \frac{\cos {\omega }_k\Delta t_k-1}{{\omega }_k}\\ 0& \cos {\omega }_k\Delta t_k& 0& -\sin {\omega }_k\Delta t_k\\ 0& \frac{1-\cos {\omega }_k\Delta t_k}{{\omega }_k}& 1& \frac{\sin {\omega }_k\Delta t_k}{{\omega }_k}\\ 0& \sin {\omega }_k\Delta t_k& 0& \cos {\omega }_k\Delta t_k\end{array}\right]$$

(17)

where $\Delta t_k=t_k-t_{k-1}$ is the interval between times $k$ and $k-1$. Using the state transition matrix in Equation (17), we obtain the state vector of the target in the ${\mathrm{x}}^{\mathrm{e}}$ − ${\mathrm{y}}^{\mathrm{e}}$ Cartesian coordinate system at time $k$, as:

$${\mathit{X}}_k^e=\left[\begin{array}{c}{x}_k^e\\ {\dot{x}}_k^e\\ y_k^e\\ {\dot{y}}_k^e\end{array}\right]=\mathit{F}({\omega }_k){\mathit{X}}_{k-1}^e=\left[\begin{array}{c}\frac{\sqrt{{\dot{x}}_{k-1}^2+{\dot{y}}_{k-1}^2}\sin {\omega }_k\Delta t_k}{{\omega }_k}\\ \sqrt{{\dot{x}}_{k-1}^2+{\dot{y}}_{k-1}^2}\cos {\omega }_k\Delta t_k\\ \frac{\sqrt{{\dot{x}}_{k-1}^2+{\dot{y}}_{k-1}^2}(1-\cos {\omega }_k\Delta t_k)}{{\omega }_k}\\ \sqrt{{\dot{x}}_{k-1}^2+{\dot{y}}_{k-1}^2}\sin {\omega }_k\Delta t_k\end{array}\right]$$

(18)

Based on Equation (18), the position coordinates of point $E$ in the ${\mathrm{x}}^{\mathrm{e}}$ − ${\mathrm{y}}^{\mathrm{e}}$ Cartesian coordinate system are given by:

$$\left[\begin{array}{c}{x}_k^e\\ y_k^e\end{array}\right]=\left[\begin{array}{c}\frac{\sqrt{{\dot{x}}_{k-1}^2+{\dot{y}}_{k-1}^2}\sin {\omega }_k\Delta t_k}{{\omega }_k}\\ \frac{\sqrt{{\dot{x}}_{k-1}^2+{\dot{y}}_{k-1}^2}(1-\cos {\omega }_k\Delta t_k)}{{\omega }_k}\end{array}\right]$$

(19)

and the angle ${\beta }_k$ in Figure 2 is given by:

$${\beta }_k=arc\tan \frac{y_k^e}{x_k^e}=arc\tan \frac{1-\cos {\omega }_k\Delta t_k}{\sin {\omega }_k\Delta t_k}$$

(20)

We assume that a sensor observes the position of the target and we use ($r_x$, $r_y$) to denote the position measurement of the target in the x − y Cartesian coordinate system at time $k$. Obviously, the position measurement of the target at time $k$ is its position coordinates at this time, if no measurement error appears in the measurement of the sensor. Therefore, we have:

$$\left[\begin{array}{c}{r}_x\\ r_y\end{array}\right]=\left[\begin{array}{c}{x}_k\\ y_k\end{array}\right]\mathrm{or}\left[\begin{array}{c}{r}_x^e\\ r_x^e\end{array}\right]=\left[\begin{array}{c}{x}_k^e\\ y_k^e\end{array}\right]$$

(21)

where ($r_x^e$,$r_y^e$) is given by:

$$\left[\begin{array}{c}{r}_x^e\\ r_y^e\end{array}\right]=\left[\begin{array}{cc}\cos {\alpha }_{k-1}& \sin {\alpha }_{k-1}\\ -\sin {\alpha }_{k-1}& \cos {\alpha }_{k-1}\end{array}\right]\times \left(\left[\begin{array}{c}{r}_x\\ r_y\end{array}\right]-\left[\begin{array}{c}{x}_{k-1}\\ y_{k-1}\end{array}\right]\right)$$

(22)

Replacing ($x_k^e$, $y_k^e$) in Equation (20) with ($r_x^e$, $r_y^e$) in Equation (22), we obtain:

$$arc\tan \frac{r_y^e}{r_x^e}=arc\tan \frac{1-\cos {\omega }_k\Delta t_k}{\sin {\omega }_k\Delta t_k}$$

(23)

Let $c=\frac{r_y^e}{r_x^e}$ and $\varphi ={\omega }_k\Delta t_k$, we have:

$$(c^2+1){\cos }^2\varphi -2\cos \varphi -c^2+1=0$$

(24)

Solving Equation (24), we obtain:

$${\omega }_k=\{\begin{array}{l}\frac{1}{\Delta t_k}arc\cos \frac{{(r_x^e)}^2-{(r_y^e)}^2}{{(r_x^e)}^2+{(r_y^e)}^2},if{r}_y^e\ge 0\\ -\frac{1}{\Delta t_k}arc\cos \frac{{(r_x^e)}^2-{(r_y^e)}^2}{{(r_x^e)}^2+{(r_y^e)}^2},if{r}_y^e<0\end{array}$$

(25)

where:

$$\{\begin{array}{l}{r}_x^e=(r_x-x_{k-1})\cos {\alpha }_{k-1}+(r_y-y_{k-1})\sin {\alpha }_{k-1}\\ r_y^e=-(r_x-x_{k-1})\sin {\alpha }_{k-1}+(r_y-y_{k-1})\cos {\alpha }_{k-1}\end{array}$$

(26)

Obviously, if no measurement error appears in the position measurement of the target, ${\omega }_k$ in Equation (25) is its turn rate from time $k-1$ to time $k$. Otherwise, we use ${\omega }_k$ in Equation (25) as the estimation of its turn rate from time $k-1$ to time $k$. Thus, by using the state vector of a target at a previous time, and its position measurement at the current time, we may estimate its turn rate at the current time.

4. MB Filter with an Adaptive Estimation of Turn Rate

In [19,20], Liu presented the MB filter to track multiple targets in the presence of clutter and noise. In this section, we apply the estimation method of the turn rate to the MB filter in [20], in order to develop the MB filter with an adaptive estimation of the turn rate. This filter consists of the following four steps:

Step 1: Prediction

In this step, we predict the marginal distribution and existence probability of each target at the current time, according to its marginal distribution, existence probability, and the estimation of the turn rate at a previous time.

Let $N({\mathit{x}}_{i,k-1};{\mathit{m}}_{i,k-1},{\mathit{P}}_{i,k-1})$, ${\rho }_{i,k-1}$, and ${\omega }_{i,k-1}$ denote the marginal distribution of target $i$, its existence probability, and the estimation of the turn rate at time $k-1$, respectively, where $i=1,2,\cdots ,N_{k-1}$ $N_{k-1}$ is the target number and $N(\cdot ;\mathit{m},\mathit{P})$ is a Gaussian distribution with mean vector $\mathit{m}$ and covariance matrix $\mathit{P}$. The predicted marginal distribution, existence probability, and turn rate of each target at time $k$ are given by:

$$N({\mathit{x}}_{i,k};{\mathit{m}}_{i,k|k-1},{\mathit{P}}_{i,k|k-1}),i=1,2,\cdots ,N_{k-1}$$

(27)

$${\rho }_{i,k|k-1}=p_{S,k}{\rho }_{i,k-1},i=1,2,\cdots ,N_{k-1}$$

(28)

$${\omega }_{i,k|k-1}={\omega }_{i,k-1},i=1,2,\cdots ,N_{k-1}$$

(29)

where $p_{S,k}$ is the survival probability, and ${\mathit{m}}_{i,k|k-1}$ and ${\mathit{P}}_{i,k|k-1}$ are given by:

$${\mathit{m}}_{i,k|k-1}=\mathit{F}({\omega }_{i,k-1}){\mathit{m}}_{i,k-1}$$

(30)

$${\mathit{P}}_{i,k|k-1}={\mathit{Q}}_{i,k-1}+\mathit{F}({\omega }_{i,k-1}){\mathit{P}}_{i,k-1}{\mathit{F}}^{\mathrm{T}}({\omega }_{i,k-1})$$

(31)

where T denotes the transpose, ${\mathit{Q}}_{i,k-1}$ is the covariance of process noise, and $\mathit{F}({\omega }_{i,k-1})$ is given by:

$$\mathit{F}({\omega }_{i,k-1})=\left[\begin{array}{cccc}1& \frac{\sin ({\omega }_{i,k-1}\Delta t_k)}{{\omega }_{i,k-1}}& 0& -\frac{1-\cos ({\omega }_{i,k-1}\Delta t_k)}{{\omega }_{i,k-1}}\\ 0& \cos ({\omega }_{i,k-1}\Delta t_k)& 0& -\sin ({\omega }_{i,k-1}\Delta t_k)\\ 0& \frac{1-\cos ({\omega }_{i,k-1}\Delta t_k)}{{\omega }_{i,k-1}}& 1& \frac{\sin ({\omega }_{i,k-1}\Delta t_k)}{{\omega }_{i,k-1}}\\ 0& \sin ({\omega }_{i,k-1}\Delta t_k)& 0& \cos ({\omega }_{i,k-1}\Delta t_k)\end{array}\right]$$

(32)

Step 2: Estimation of Turn Rate

In this step, we use the measurements at time $k$, and marginal distribution of target $i$ at time $k-1$, to estimate its turn rate from time $k-1$ to time $k$.

Let ${\mathit{y}}_{j,k}={\left[\begin{array}{cc}{r}_{x,k}^j& r_{y,k}^j\end{array}\right]}^{\mathrm{T}}$ denote the measurement at time $k$, where $j=1,2,\cdots ,M_k$, $M_k$ is the number of measurements, and $r_{x,k}^j$ and $r_{y,k}^j$ denote the $x$ and $y$ components of measurement ${\mathit{y}}_{j,k}$; and let ${\mathit{m}}_{i,k-1}={\left[\begin{array}{cccc}{\eta }_{x,k-1}^i& {\dot{\eta }}_{x,k-1}^i& {\eta }_{y,k-1}^i& {\dot{\eta }}_{y,k-1}^i\end{array}\right]}^{\mathrm{T}}$ denote the state vector of target $i$ at time $k-1$, where $({\eta }_{x,k-1}^i,{\eta }_{y,k-1}^i)$ denotes the position coordinates of target $i$ and $({\dot{\eta }}_{x,k-1}^i,{\dot{\eta }}_{y,k-1}^i)$ denotes its velocities. According to Equation (25), the turn rate of target $i$ that corresponds with measurement ${\mathit{y}}_{j,k}$, is given by:

$${\omega }_{k,e}^{i,j}=\frac{\mathrm{sgn}(y_{i,j}^e)}{t_k-t_{k-1}}\arccos \frac{{(x_{i,j}^e)}^2-{(y_{i,j}^e)}^2}{{(x_{i,j}^e)}^2+{(y_{i,j}^e)}^2},i=1,2,\cdots ,N_{k-1},j=1,2,\cdots ,M_k$$

(33)

where

$$\left[\begin{array}{c}{x}_{i,j}^e\\ y_{i,j}^e\end{array}\right]=\left[\begin{array}{cc}\cos {\alpha }_{i,k-1}& \sin {\alpha }_{i,k-1}\\ -\sin {\alpha }_{i,k-1}& \cos {\alpha }_{i,k-1}\end{array}\right]\times \left(\left[\begin{array}{c}{r}_{x,k}^j\\ r_{y,k}^j\end{array}\right]-\left[\begin{array}{c}{\eta }_{x,k-1}^i\\ {\eta }_{y,k-1}^i\end{array}\right]\right)$$

(34)

$${\alpha }_{i,k-1}=arc\cos \frac{{\dot{\eta }}_{x,k-1}^i}{\sqrt{{({\dot{\eta }}_{x,k-1}^i)}^2+{({\dot{\eta }}_{y,k-1}^i)}^2}}$$

(35)

$$\mathrm{sgn}(y_{i,j}^e)=\{\begin{array}{ll}1,& if{y}_{i,j}^e\ge 0\\ -1,& if{y}_{i,j}^e<0\end{array}$$

(36)

Considering that the turn rate of a target is generally within a known range, the turn rate of target $i$ that corresponds with measurement ${\mathit{y}}_{j,k}$, may be given by:

$${\omega }_k^{i,j}=\{\begin{array}{ll}{\omega }_{\max }& if{\omega }_{k,e}^{i,j}\ge {\omega }_{\max }\\ {\omega }_{k,e}^{i,j}& if-{\omega }_{\max }<{\omega }_{k,e}^{i,j}<{\omega }_{\max }\\ -{\omega }_{\max }& if{\omega }_{k,e}^{i,j}\le -{\omega }_{\max }\end{array}$$

(37)

where ${\omega }_{\max }$ is the maximal turn rate.

Step 3: Update

In this step, we use the marginal distribution $N({\mathit{x}}_{i,k-1};{\mathit{m}}_{i,k-1},{\mathit{P}}_{i,k-1})$, predicted existence probability ${\rho }_{i,k|k-1}$, turn rate estimation ${\omega }_k^{i,j}$, and measurement ${\mathit{y}}_{j,k}$, to obtain the updated marginal distribution, existence probability, and turn rate. The updated distribution and existence probability of target $i$ that corresponds with measurement ${\mathit{y}}_{j,k}$, are as follows:

$$N({\mathit{x}}_{i,k};{\mathit{m}}_k^{i,j},{\mathit{P}}_k^{i,j}),i=1,2,\cdots ,N_{k-1},j=1,2,\cdots ,M_k$$

(38)

$${\rho }_k^{i,j}=\frac{p_{D,k}{\rho }_{i,k|k-1}N({\mathit{y}}_{j,k};{\mathit{H}}_k{\mathit{m}}_{k|k-1}^{i,j},{\mathit{H}}_k{\mathit{P}}_{k|k-1}^{i,j}{\mathit{H}}_k^{\mathrm{T}}+{\mathit{R}}_k)}{{\lambda }_{c,k}+p_{D,k}{\displaystyle \sum _{e=1}^{N_{k-1}}{\rho }_{e,k|k-1}N({\mathit{y}}_{j,k};{\mathit{H}}_k{\mathit{m}}_{k|k-1}^{e,j},{\mathit{H}}_k{\mathit{P}}_{k|k-1}^{e,j}{\mathit{H}}_k^{\mathrm{T}}+{\mathit{R}}_k)}},i=1,2,\cdots ,N_{k-1},j=1,2,\cdots ,M_k$$

(39)

where:

$${\mathit{m}}_{k|k-1}^{i,j}=\mathit{F}({\omega }_k^{i,j}){\mathit{m}}_{i,k-1}$$

(40)

$${\mathit{P}}_{k|k-1}^{i,j}={\mathit{Q}}_{i,k-1}+\mathit{F}({\omega }_k^{i,j}){\mathit{P}}_{i,k-1}{\mathit{F}}^{\mathrm{T}}({\omega }_k^{i,j})$$

(41)

$${\mathit{m}}_k^{i,j}={\mathit{m}}_{k|k-1}^{i,j}+{\mathit{A}}_{i,j}\cdot ({\mathit{y}}_{j,k}-{\mathit{H}}_k{\mathit{m}}_{k|k-1}^{i,j})$$

(42)

$${\mathit{P}}_k^{i,j}=(\mathit{I}-{\mathit{A}}_{i,j}\cdot {\mathit{H}}_k){\mathit{P}}_{k|k-1}^{i,j}$$

(43)

$${\mathit{A}}_{i,j}={\mathit{P}}_{k|k-1}^{i,j}{\mathit{H}}_k^{\mathrm{T}}{[{\mathit{H}}_k{\mathit{P}}_{k|k-1}^{i,j}{\mathit{H}}_k^{\mathrm{T}}+{\mathit{R}}_k]}^{-1}$$

(44)

where ${\mathit{H}}_k$ and ${\mathit{R}}_k$ are the observation matrix and covariance matrix of the observation noise, respectively, $\mathit{I}$ denotes the identity matrix, ${\lambda }_{c,k}$ denotes the clutter rate, and $p_{D,k}$ denotes the detection probability.

Equations (27) and (38) indicate that a marginal distribution $N({\mathit{x}}_{i,k-1};{\mathit{m}}_{i,k-1},{\mathit{P}}_{i,k-1})$ at time $k-1$ generates a prediction distribution $N({\mathit{x}}_{i,k};{\mathit{m}}_{i,k|k-1},{\mathit{P}}_{i,k|k-1})$ and $M_k$ update distributions $N({\mathit{x}}_{i,k};{\mathit{m}}_k^{i,j},{\mathit{P}}_k^{i,j})$, $j=1,2,\cdots ,M_k$ at time $k$. To decrease the computational load, we merge these $M_k+1$ distributions to a singular distribution. The merging procedures are as follows:

Let $N({\mathit{x}}_{i,k};{\mathit{m}}_k^{i,j},{\mathit{P}}_k^{i,j})$, $j=1,2,\cdots ,M_k+1$ denote the $M_k+1$ distributions where $N({\mathit{x}}_{i,k};{\mathit{m}}_k^{i,M_k+1},{\mathit{P}}_k^{i,M_k+1})=N({\mathit{x}}_{i,k};{\mathit{m}}_{i,k|k-1},{\mathit{P}}_{i,k|k-1})$. Similarly, let ${\rho }_k^{i,j}$, $j=1,2,\cdots ,M_k+1$ and ${\omega }_k^{i,j}$, $j=1,2,\cdots ,M_k+1$ denote the corresponding existence probabilities and turn rates, respectively, where ${\rho }_k^{i,M_k+1}={\rho }_{i,k|k-1}$ and ${\omega }_k^{i,M_k+1}={\omega }_{i,k|k-1}$. We first find the index of the maximal existence probability from ${\rho }_k^{i,j}$, $j=1,2,\cdots ,M_k+1$, namely:

$$q=\arg \underset{j\in \{1,\cdots ,M_k+1\}}{\max }\{{\rho }_k^{i,j}\}$$

(45)

We then use the distribution, existence probability, and turn rate with index $q$, as the marginal distribution of target $i$, its existence probability, and the turn rate at time $k$, respectively, namely:

$$N({\mathit{x}}_{i,k};{\mathit{m}}_{i,k},{\mathit{P}}_{i,k})=N({\mathit{x}}_{i,k};{\mathit{m}}_k^{i,q},{\mathit{P}}_k^{i,q}),i=1,2,\cdots ,N_{k-1}$$

(46)

$${\rho }_{i,k}={\rho }_k^{i,q},i=1,2,\cdots ,N_{k-1}$$

(47)

$${\omega }_{i,k}={\omega }_k^{i,q},i=1,2,\cdots ,N_{k-1}$$

(48)

Step 4: Generation of New Target Distribution and the Output of the Filter

In this step, we use the measurement at time $k$ to generate the marginal distribution of the new target as:

$$N({\mathit{x}}_{j,k};{\mathit{m}}_{\gamma ,k}^j,{\mathit{P}}_{\gamma ,k}^j),j=1,2,\cdots ,M_k$$

(49)

where ${\mathit{m}}_{\gamma ,k}^j$ is from measurement ${\mathit{y}}_{j,k}={\left[\begin{array}{cc}{r}_{x,k}^j& r_{y,k}^j\end{array}\right]}^{\mathrm{T}}$, and ${\mathit{m}}_{\gamma ,k}^j={\left[\begin{array}{cccc}{r}_{x,k}^j& 0& r_{y,k}^j& 0\end{array}\right]}^{\mathrm{T}}$ and ${\mathit{P}}_{\gamma ,k}^j={\mathit{P}}_{\gamma }$ where ${\mathit{P}}_{\gamma }$ is a known covariance matrix. Meanwhile, we designate parameter ${\rho }_{\gamma }$ as the existence probability of the new target and assign 0 as its turn rate, namely:

$${\rho }_{\gamma ,k}^j={\rho }_{\gamma },j=1,2,\cdots ,M_k$$

(50)

$${\omega }_{\gamma ,k}^j=0,j=1,2,\cdots ,M_k$$

(51)

We then combine the marginal distributions of the existing targets in Equation (46), with those of new targets in Equation (49), to form the marginal distributions of individual targets at time $k$, as:

$${\left\{N({\mathit{x}}_{i,k};{\mathit{m}}_{i,k},{\mathit{P}}_{i,k})\right\}}_{i=1}^{N_k}={\left\{N({\mathit{x}}_{i,k};{\mathit{m}}_{i,k},{\mathit{P}}_{i,k})\right\}}_{i=1}^{N_{k-1}}\cup {\left\{N({\mathit{x}}_{j,k};{\mathit{m}}_{\gamma ,k}^j,{\mathit{P}}_{\gamma ,k}^j)\right\}}_{j=1}^{M_k}$$

(52)

where $N_k=N_{k-1}+M_k$. The corresponding existence probabilities and turn rates of individual targets at time $k$, are as follows:

$${\left\{{\rho }_{i,k}\right\}}_{i=1}^{N_k}={\left\{{\rho }_{i,k}\right\}}_{i=1}^{N_{k-1}}\cup {\left\{{\rho }_{\gamma ,k}^j\right\}}_{j=1}^{M_k}$$

(53)

$${\left\{{\omega }_{i,k}\right\}}_{i=1}^{N_k}={\left\{{\omega }_{i,k}\right\}}_{i=1}^{N_{k-1}}\cup {\left\{{\omega }_{\gamma ,k}^j\right\}}_{j=1}^{M_k}$$

(54)

After this combination, we prune the targets whose existence probability ${\rho }_{i,k}$ is less than threshold $\tau$, and propagate the marginal distributions, existence probabilities, and turn rates of the remaining targets to the next time step. Those targets whose existence probability ${\rho }_{i,k}$ is greater than 0.5 are picked as the output of the filter at time $k$.

5. Simulation Results

In this section, we use an example to reveal the tracking performance of the MB filter with an adaptive estimation of the turn rate for multiple maneuvering targets. In this example, Pasha’s PHD filter [16] is used as the comparison object, and the OSPA distance [21], with parameters $c=50\mathrm{m}$ and $p=2$, is used as the measure. The covariance matrix ${\mathit{Q}}_{i,k-1}$, observation matrix ${\mathit{H}}_k$, and covariance matrices ${\mathit{R}}_k$ and ${\mathit{P}}_{\gamma }$ used in the experiment, are as follows:

$${\mathit{Q}}_{i,k-1}=\mathit{Q}({\sigma }_v)=\left[\begin{array}{cccc}\frac{\Delta t_k^4}{4}& \frac{\Delta t_k^3}{2}& 0& 0\\ \frac{\Delta t_k^3}{2}& \Delta t_k^2& 0& 0\\ 0& 0& \frac{\Delta t_k^4}{4}& \frac{\Delta t_k^3}{2}\\ 0& 0_k& \frac{\Delta t_k^3}{2}& \Delta t_k^2\end{array}\right]{\sigma }_v^2$$

(55)

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

(56)

$${\mathit{R}}_k=\left[\begin{array}{cc}{\sigma }_w^2& 0\\ 0& {\sigma }_w^2\end{array}\right]$$

(57)

$${\mathit{P}}_{\gamma }=\left[\begin{array}{cccc}2500& 0& 0& 0\\ 0& 625& 0& 0\\ 0& 0& 2500& 0\\ 0& 0& 0& 625\end{array}\right]$$

(58)

where ${\sigma }_v$ and ${\sigma }_w$ denote the standard deviations of noises.

Three coordinated turn models with different turn rates are used in Pasha’s PHD filter. The state transition and covariance matrices for models $r_{i,k}=1$, $r_{i,k}=2$, and $r_{i,k}=3$ are given by ${\mathit{F}}_{k-1}(r_{i,k}=1)=\mathit{F}({\omega }_k=0^{\circ }{\mathrm{s}}^{-1})$, ${\mathit{Q}}_{k-1}(r_{i,k}=1)=\mathit{Q}({\sigma }_v=1{\mathrm{ms}}^{-2})$, ${\mathit{F}}_{k-1}(r_{i,k}=2)=\mathit{F}({\omega }_k=5^{\circ }{\mathrm{s}}^{-1})$, ${\mathit{Q}}_{k-1}(r_{i,k}=2)=\mathit{Q}({\sigma }_v=3{\mathrm{ms}}^{-2})$, ${\mathit{F}}_{k-1}(r_{i,k}=3)=\mathit{F}({\omega }_k=-5^{\circ }{\mathrm{s}}^{-1})$ and ${\mathit{Q}}_{k-1}(r_{i,k}=3)=\mathit{Q}({\sigma }_v=3{\mathrm{ms}}^{-2})$, respectively. The Markov transition probabilities among different motion models are given by:

$$[t_{k|k-1}(r_{i,k}|r_{i,k-1})]=\left[\begin{array}{ccc}0.9& 0.05& 0.05\\ 0.05& 0.9& 0.05\\ 0.05& 0.05& 0.9\end{array}\right]$$

(59)

$\mathit{H}(r_{i,k})$ and $\mathit{R}(r_{i,k})$ of Pasha’s PHD filter are given by $\mathit{H}(r_{i,k})={\mathit{H}}_k$ and $\mathit{R}(r_{i,k})={\mathit{R}}_k$.

Example 1.

Five targets are considered in this example. Targets 1, 2, 3, and 4 appear at $t=1\mathrm{s}$, $t=1\mathrm{s}$, $t=3\mathrm{s}$, and $t=3\mathrm{s}$, respectively, and disappear at $t=70\mathrm{s}$. Target 5 appears at $t=5\mathrm{s}$ and disappears at $t=60\mathrm{s}$. Each target changes its turn rate at $t=15\mathrm{s}$, $t=30\mathrm{s}$, $t=40\mathrm{s}$, and $t=55\mathrm{s}$, respectively. The initial positions and moving trajectories of these five targets are shown in Figure 3.

We use parameters $\Delta t_k=1\mathrm{s}$, ${\sigma }_v=0{\mathrm{m}/\mathrm{s}}^2$, ${\sigma }_w=1\mathrm{m}$, $p_{S,k}=1.0$, ${\lambda }_{c,k}=1.25\times {10}^{-5}{\mathrm{m}}^{-2}$, and $p_{D,k}=0.95$ to generate the simulation measurement. Figure 4 shows the simulation measurement for a trial. Setting the parameters of the proposed filter to $\Delta t_k=1\mathrm{s}$, ${\sigma }_v=1\mathrm{m}/\mathrm{s}$, $p_{S,k}=0.6$, ${\lambda }_{c,k}=1.25\times {10}^{-5}{\mathrm{m}}^{-2}$, $p_{D,k}=0.95$, $\tau =0.001$, ${\sigma }_w=2\mathrm{m}$, ${\rho }_{\gamma }=0.1$, and ${\omega }_{\max }=6{}^{\mathrm{o}}/\mathrm{s}$, and the parameters of Pasha’s filter to $\Delta t_k=1\mathrm{s}$, ${\sigma }_v=1\mathrm{m}/\mathrm{s}$, $p_{S,k}=1.0$, ${\lambda }_{c,k}=1.25\times {10}^{-5}{\mathrm{m}}^{-2}$, $p_{D,k}=0.95$, $\tau =0.001$, ${\sigma }_w=2\mathrm{m}$, and ${\rho }_{\gamma }=0.1$, respectively, we use the proposed filter and Pasha’s filter to address the simulation measurements for 150 trials. The experimental results are shown in Figure 5. Based on these experimental results, the proposed filter performs better than Pasha’s filter, most of the time. Two factors are responsible for this result. The first factor is the difference between the actual motion mode of a target, and the dynamic model used by the filter. This difference causes Pasha’s filter to provide an inaccurate state estimation of the target. The proposed filter reduces this difference by estimating the turn rate of the target at each given time. The second factor is the filter’s memory. Due to the poor memory of Pasha’s filter, it is prone to discarding the information of a target from the posterior intensity, and cannot provide its state estimation if the target is not detected by a sensor at each point in time. In contrast to Pasha’s filter, the proposed filter provides the state estimation of a missed target, due to its sufficient memory of the target. The effect of the filter’s memory on the OSPA distance has been discussed in detail in [20]. As shown in Figure 5, a peak appears at $t=60\mathrm{s}$, because the proposed filter furnishes the state estimation of target 5 at its disappearing time. According to the definition of OSPA distance in [21], the OSPA distance is used to measure the similarity between two different sets. The excessive or deficient state estimation of the target will be punished with the cutoff distance.

To reveal the effect of the clutter rate and detection probability on the tracking performance of the proposed filter, we use different clutter rates and detection probabilities to generate simulation measurements, and use the proposed filter and Pasha’s filter to address the simulation measurements, respectively, for 150 trials.

Table 1. OSPA distance (m) at five clutter rates.

${\mathit{\lambda }}_{\mathit{c}\mathbf{,}\mathit{k}}$ ($\mathbf{\times }{\mathbf{10}}^{\mathbf{-}\mathbf{6}}{\mathbf{m}}^{\mathbf{-}\mathbf{2}}$)	3.125	6.25	9.375	12.5	15.63
Proposed filter	5.18	5.61	6.05	6.32	6.48
Pasha’s filter	13.03	14.37	15.22	15.97	16.61

and

Table 2. OSPA distance (m) at five detection probabilities.

${\mathit{p}}_{\mathit{D}\mathbf{,}\mathit{k}}$	1.0	0.95	0.9	0.85	0.8
Proposed filter	5.04	6.32	7.79	11.67	12.95
Pasha’s filter	8.78	15.97	21.65	26.14	30.40

show the results obtained from different clutter rates and detection probabilities. Table 1 suggests that an increase in the clutter rate leads to a larger OSPA distance for both the proposed filter and Pasha’s filter, but the proposed filter performs better at each clutter rate than Pasha’s filter. A similar conclusion is also reached from the result in Table 2. A decrease in the detection probability enlarges the OSPA distance of the proposed filter and Pasha’s filter, but the proposed filter obtains a smaller OSPA distance than Pasha’s filter at each detection probability.

The performance time is also an important measure for the performance of a filter.

Table 3. Performance time (s) at five clutter rates.

${\mathit{\lambda }}_{\mathit{c}\mathbf{,}\mathit{k}}$ ($\mathbf{\times }{\mathbf{10}}^{\mathbf{-}\mathbf{6}}{\mathbf{m}}^{\mathbf{-}\mathbf{2}}$)	3.125	6.25	9.375	12.5	15.63
Proposed filter	3.14	4.96	7.70	10.70	14.37
Pasha’s filter	1.08	1.98	3.67	6.33	10.49

displays the required time of a trial for the proposed filter and Pasha’s filter, at different clutter rates. Based on Table 3, the proposed filter requires more time than Pasha’s filter for each trial, because the proposed filter is used to estimate the turn rate of the target at each time step, and this estimation requires a number of calculations.

6. Conclusions

In this study, the formula for calculating the turn rate of a maneuvering target is derived. Based on this formula, we may estimate the turn rate of a target by using its state vector at a previous time, and its measurement at the current time. Applying the estimation method of the turn rate to the MB filter, we present a MB filter with an adaptive estimation of the turn rate, to track multiple targets maneuvering at a random turn rate. Based on simulation experimental data, we test the performance of the proposed filter, by comparing it with Pasha’s filter. The experimental results suggest that the proposed filter is better than Pasha’s filter at tracking the targets that maneuver at a variable turn rate.