A shrinkage probability hypothesis density filter for multitarget tracking

The target state is ${{\bf{x}}_{k}}={\left[{{x_{k}},{{\dot{x}}_{k}},{y_{k}},{{\dot{y}}_{k}}}\right]^{T}}$,where $({x_{k}},{y_{k}})$ and $({\dot{x}_{k}},{\dot{y}_{k}})$ are the position and velocity. Since there is no ordering on the respective collections of all target states, they can be naturally represented as a finite set. ${X_{k}}$ is the multitarget state-set at time step $k$, i.e., the set of unknown target states (which are also of unknown number).

In a single-target scenario, the state ${{\bf{x}}_{k}}$ is modeled as random vectors. Then, the multitarget state, including target motion, birth, spawning, can be described by RFS. For target motion, given multitarget state set ${X_{k-1}}$, each ${{\bf{x}}_{k-1}}\in{X_{k-1}}$ either survives at time step $k$ with probability ${e_{k\left|{k-1}\right.}}({{\bf{x}}_{k-1}})$, and its transition probability density from ${x_{k-1}}$ to $x_{k}$ is ${f_{k\left|{k-1}\right.}}({{\bf{x}}_{k}}\left|{{{\bf{x}}_{k-1}}}\right.)$ . Therefore, the target motion is modeled as the RFS ${S_{k\left|{k-1}\right.}}({{\bf{x}}_{k-1}})$. In the same way, when the RFS of target birth at time $k$ is modeled by ${\Gamma_{k}}$, and the RFS of targets spawning from a target with ${{\bf{x}}_{k-1}}$ is modeled by ${B_{k\left|{k-1}\right.}}({{\bf{x}}_{k-1}})$ , the multitarget stat ${X_{k}}$ is given by

The measurements are measurements of reflected power as in Ref8 . ${n_{k}}$ is white complex Gaussian noise with variance $\sigma_{0}^{2}$. In this section, assume that the intensity of all targets is ${I_{k}}$, the SNR for the targets is defined by

The method to deal with multiple targets with different or unknown SNRs is discussed in section 5.

Because of the assumptions made in the modeling process and the essential difference between an RFS and a random vector, the two conditions that should be satisfied to model TBD measurements by a RFS are the following: there is no target that generates more than one measurement vector, and no measurement vector is generated by more than one targetRef22 . In summary, the measurement of targets should be modeled by a ’nail-like’ model on the range-Doppler-Bearing maps (the $ijk$ cell is defined by coordinat ($r_{i},d_{j},b_{l}$)) as follow:

R, D and B are the size of a range, the Doppler and the bearing cell. $r_{k}=\sqrt{{{(x_{k})}^{2}}+{{(y_{k})}^{2}}}$, $d_{k}=\frac{{x_{k}\dot{x}_{k}+y_{k}\dot{y}_{k}}}{{\sqrt{{{(x_{k})}^{2}}+{{(y_{k})}^{2}}}}}$ and $b_{k}=\arctan\left({\frac{{y_{k}}}{{x_{k}}}}\right)$. $R_{k}^{i}=\left[{{r_{i}}-\frac{R}{2},{r_{i}}+\frac{R}{2}}\right)$, $D_{k}^{j}=\left[{{d_{j}}-\frac{D}{2},{d_{j}}+\frac{D}{2}}\right)$ and $B_{k}^{l}=\left[{{b_{l}}-\frac{B}{2},{b_{l}}+\frac{B}{2}}\right)$.

Then the measurement model is

This means the measurement of a target is like a nail on the range-Doppler-bearing maps. This ’nail-like’ model is similar to the point target model in Ref5 and Ref6 . At step k, the measurement provided by the sensor consists of $N={N_{r}}\times{N_{d}}\times{N_{b}}$ measurements $z_{k}^{ijl}$, where ${N_{r}}$, ${N_{d}}$ and ${N_{b}}$ are the number of range, Doppler and bearing cells.

The number of $z_{k}^{ijl}$ is a constant $N$. However, the number of element of an RFS should be random and Poisson distributed for a PHD filter. Because weak signal information should be preserved by TBD algorithms, we should make sure all measurements produced by targets are included in the RFS. Therefore, a threshold could be chosen as follows:

where the likelihood function for the target is

and that for noise is

where ${{\rm I}_{0}}\left({}\right)$ is the zero order Bessel function and $z_{k}^{ijl}$ is assumed positive.The ${p_{D}}({{\bf{x}}_{k}})$ is the probability detection. Because this algorithm is for TBD, it should be ensured that ${p_{D}}({{\bf{x}}_{k}})\approx 1$.

$Z_{k}$ is the observation-set consisting of all measurements collected by all sensors at time-step $k$, no matter the measurement from the targets or from the noise. If $z_{k}^{ijl}\geq{\theta_{k}}$, let ${z_{k}}=z_{k}^{ijl}$. Then, the measurement RFS $Z_{k}$ is constructed by a subset of the $z_{k}^{ijl}$ using a thresholding mechanism. The threshold can insure that ${p_{D}}({{\bf{x}}_{k}})\approx 1$ by (5). The threshold is a function of the SNR of the targets. The element in the RFS $Z_{k}$ is as follows:

and the likelihood function according

At the same time, the sensor could be modeled by

where the measurements of targets are modeled by RFS ${\Theta_{k}}({x_{k}})$ and the noise is modeled as RFS ${K_{k}}$. The elements in ${\Theta_{k}}({x_{k}})$ are the measurements produced by the targets. The elements in ${K_{k}}$ are the measurements which are produced by the noise and bigger than ${\theta_{k}}$ as well. It is shown that when the measurements of TBD are modeled by RFSs, multitarget TBD can be regarded as a kind of classification problem of target and noise measurements for one scan. This classification problem will be further analyzed in Section 4.

As proposed in Ref13 , the PHD prediction and update equations are presented as follows:

where ${D_{k\left|k\right.}}({\bf{x}}\left|{{Z_{1:k}}}\right.)$ the PHD is the density whose integral $\int_{S}{{D_{k\left|k\right.}}({\bf{x}}\left|{{Z_{1:k}}}\right.)d{\bf{x}}}$ on any region S of state space is $\hat{n}_{k}(S)=\int{\left|{X\cap S}\right|}{p_{k}}({X_{k}}\left|{{Z_{1:k}}}\right.)\delta X$. ${b_{k\left|{k-1}\right.}}({{\bf{x}}_{k}}\left|{{{\bf{x}}_{k-1}}}\right.)$ and ${\gamma_{k}}({{\bf{x}}_{k}}){\rm{}}$ denotes the intensity of ${B_{k\left|{k-1}\right.}}({{\bf{x}}_{k-1}})$ and ${\Gamma_{k}}$ at time $k$, and ${\kappa_{k}}({{\bf{z}}_{k}})$ is the intensity of ${K_{k}}$. $Z_{1:k}$ is the time-sequence of observation-sets.

Because the PHD propagation equations involve multiple integrals that have no computationally tractable closed form expressions, Sequential Monte Carlo (SMC) methods are used to approximate the PHD in Ref12 . Let the update PHD at $k-1$ step ${D_{k-1\left|{k{\rm{-1}}}\right.}}({\bf{x}}\left|{{Z_{1:k-1}}}\right.)$ be represented by a set of particles $\left\{{\omega_{k-1}^{(p)},{\bf{x}}_{k-1}^{(p)}}\right\}_{p=1}^{{L_{k-1}}}$ , as

The predicted particles are generated by

where ${q_{k}}(\bullet\left|{{\bf{x}}_{k-1}^{(p)}}\right.)$ and ${v_{k}}(\bullet)$ are proposal density. The predicted density is

where

Set ${p_{D}}({\bf{x}})\equiv 1$, the update density will be

where

According to the standard treatment of particle filter, resample $\left\{{\omega_{k}^{*(p)}/\hat{n}_{k},{\bf{x}}_{\left.k\right|k-1}^{(p)}}\right\}_{p=1}^{{L_{k-1}}+{J_{k}}}$ to get $\left\{{\omega_{k}^{(p)}/\hat{n}_{k},{\bf{x}}_{k}^{(p)}}\right\}_{p=1}^{{L_{k}}}$

Because ${\kappa_{k}}({\bf{z}}_{k})$ can be presented by the number of noise samples ${\lambda_{k}}$ multiplying the clutter probability density, the update operator (19) turns into

Equation (21) is the likelihood function for the measurement of noise in the measurements set, which denotes the clutter probability density in theoretical filtering (Ref12 and Ref13 ), and the optimal choice of $\sigma$ should be the variance of noise , as shown in (20). However, for TBD applications, the SNR is extremely low. Hence, this setting of leads to the situation in which measurements generated by the target and noise are not separated correctly and heavy degradation of tracking performance is seen. In fact, the derivation of (20) depends on the hypothesis that the number of noise sample is Poisson distributed. This means the clutter follows a Poisson distribution and is independent of target-originated measurements. In this way, the PHD can be the ”best fit” approximations of the multitarget posterior probability Ref13 . The average number of noise samples included in RFS is

When $N$ is large and $\lambda_{k}$ is relatively small, the Poisson distribution could be approximated by the distribution of a number of noise samples. However, this approximation will be destroyed when become sufficiently large, which is the case in low SNR scenarios, as in Table 1.

When ${\sigma_{0}}$ in (20) cannot make PHD a sufficient approximation of the multitarget posterior probability, choosing the ”optimal” becomes the key problem for extending the PHD filter to be applied in TBD. As discussed in the following section, this selection can be considered as a shrinkage operation, and the optimal parameter for the shrinkage operation can be obtained via certain optimization procedures.

It is known that the Bayesian classifier is the optimal classifier when the posterior probability can be calculated. However, as indicated in Section 3.2, the PHD cannot be a sufficient approximation of multitarget posterior probability when the SNR of the targets is low. Therefore, another classification criterion, the Fisher class separability criterion, is introduced into our methods.

The separation between the two classes of objects defined by Fisher is the ratio of the distance between the centers of classes to the scattering of the classes Ref24 :

where ${\mu_{i}}$ represents the centers of classes and $d_{i}^{2}$ denotes the scattering of the classes. In our analysis, they are set to the mean and variance of the likelihood functions of the target and the noise, respectively.

For targets, the mean and variance are defined as follows:

For the noise, the mean and variance are given by

Therefore, the Fisher class separability is

Simulations indicated that the correlation between $S$ and the SNR of the targets is a positive, approximately linear relationship (as shown in Fig.1). Therefore, when the SNR is low, the Fisher class separability is so small that the targets and false alarms could not be distinguished. In other words, we should enlarge the Fisher class separability of the target class and the noise class when the SNR of the targets is low.

To enlarge the Fisher class separability between targets and noise, some kind of ”shrinkage” operation is applied to the likelihood function of the noise. That is, the parameter $\sigma_{0}^{2}$ in (26) and (27) is adjusted. If ${\sigma_{s}}\leq{\sigma_{0}}$ is chosen instead of $\sigma_{0}^{2}$ itself, then Fisher class separability becomes

Obviously, ${S_{s}}$ is larger than $S$ whenever ${\sigma_{s}}\leq{\sigma_{0}}$ . The following question naturally arises: how does one choose the ”optimal” ${\sigma_{s}}$ , which is the key parameter for the shrinkage process?

To demonstrate how influences the classification problem of target and noise measurements, we can define the Mahalanobis distance Ref24 of measurement $z_{k}$ using the center of the noise class and the target class, respectively:

and

Because the selection of the threshold is a difficult task for the threshold shrinkage algorithm Ref20 in the field of image denoising, the selection of the optimal parameter for the shrinkage-PHD filter involves similar problems. For the threshold shrinkage algorithm, an oversized threshold may result in the loss of signal. In contrast, too much noise may remain if an undersized threshold is used. In the same way, for the shrinkage-PHD filter, preserving the weak target information and reducing the noise by tracking are the major challenges to be addressed for TBD problems. Therefore, when considering the two sides, the selection criterion should be designed as follows:

The first constraint represents the ”shrinkage”. To preserve the weak target information and classify as many low SNR targets to the target class as possible, the Mahalanobis distance of measurement with the center of noise class should be enhanced; hence, the variance of the noise class should be reduced.

The second constraint refers to the fixed probability of target loss. The reason for this term is that for TBD problems, false alarms can be reduced by integration over time. Meanwhile, if target loss occurred, the integration of the information of targets is broken off. Hence, the cost of target loss is much larger than the cost of false alarms, which is the main difference between TBD and common radar detection. Moreover, for the PHD filter, a missed detection can result in loss of the track Ref18 . Therefore, target loss should be fixed as $\beta$. When the $\beta$ is set, a critical value $z_{s}$ is determined.

The third constraint means that the measurements $z_{k}$ determined by the second constraint should be classified as targets. When the SNR is low, ${S_{s0}}({z_{s}},{\sigma_{s}})$ always intersects ${S_{s1}}({z_{s}})$ , so the optimal solution $\sigma_{s}^{M}$ to (32) is deduced, as indicated in Fig.2(a). With the increasing SNR, ${S_{s0}}({z_{s}},{\sigma_{s}})$ becomes much larger than ${S_{s1}}({z_{s}})$ and thus we choose $\sigma_{s}^{M}={\sigma_{0}}$, as shown in Fig.2(b).

Note that although the shrinkage operation could highlight the target measurement and reduce the possibility of losing the target during tracking, the noise might be incorrectly increased unavoidably. To minimize the influence of shrinkage on the classification of the noise measurement, the maximum $\sigma_{s}$ , which is consistent with the above constraints, is obtained. The optimal solutions to (32) under different SNRs are listed in Table 2:

The update operator of the shrinkage-PHD filter can be summarized as follows:

where $\sigma_{s}^{M}$ should be set according to Table 2.

The optimal sub-pattern assignment (OSPA) distance Ref15 between the estimated and true multitarget state is adopted here to estimate error.

The OSPA distance is defined as follows Ref15 . Let ${d^{(c)}}({\bf{q}},{\bf{y}}):=\min(c,\left\|{{\bf{q}}-{\bf{y}}}\right\|)$ for ${\bf{q}},{\bf{y}}\in{\Re^{d}}$, and ${\Pi_{p}}$ denote the set of permutations on $\left\{{1,2,\ldots,p}\right\}$ for any positive integer $p$. $c$ is the cut-off value. Then, for $Q=\left\{{{{\bf{q}}_{1}},\cdots,{{\bf{q}}_{m}}}\right\}$ and $Y=\left\{{{{\bf{y}}_{1}},\cdots,{{\bf{y}}_{n}}}\right\}$,

if $m\leq n$; and${\bar{d}^{(c)}}\left({Q,Y}\right):={\bar{d}^{(c)}}\left({Y,Q}\right)$ if $m>n$; and ${\bar{d}^{(c)}}\left({Q,Y}\right)={\bar{d}^{(c)}}\left({Y,Q}\right)=0$ if $m=n=0$.

For example, when $m=1$ and $n=2$, $Q_{1}=\left\{q_{1}\right\}$ and $Y_{1}=\left\{y_{1},y_{2}\right\}$, the ${\bar{d}^{(c)}}\left({Q_{1},Y_{1}}\right)=0.5*(min(min(|q_{1}-y_{1}|,c),min(|q_{1}-y_{2}|,c))+c)$.

Consider range cells in the interval $[80000,90000]$ with meters as the unit and Doppler cells in the interval $[-400,-150]$ with meters/second as the unit. ${N_{r}}=200$, ${N_{d}}=10$ and ${N_{b}}=1$. The size of a range and Doppler cell is as follows: R=50 m and D=25 m/s. ${L_{k}}\equiv 2000$ and ${J_{k}}\equiv 800$ . The time interval T is equal to 1 s.

At time step 1, the first target appears with the initial target state ${\left[{\begin{array}[]{*{20}{c}}{89000}&{-200}&0&0\\ \end{array}}\right]^{T}}$. At step 10, a secondary target spawns from the first target, moving at -300 m/s in the x-direction. The standard deviation of the noise is 0.25.

In the first example, the SNRs of both targets are 9 dB. The measurements are shown in Fig. 3. Tracking with the shrinkage-PHD filter, we plot the PHD of the position and velocity versus time in Fig. 4 and 5. It is shown that the PHD of targets can be integrated over time by simultaneous detection and tracking. Observe that both trajectories are automatically initiated and tracked. In step 3, the first target can be tracked stably. At step 11, just 1 s after target spawning, the second target is already initiated and tracked.

For the second example, we conducted simulations to compare the performance of the shrinkage-PHD filter (shrinkage-PHDF) the PHD filter (PHDF) Ref16 , and the classical multitarget particle filtering method (MPF) Ref8 in tracking accurately and consistently, and 50 Monte Carlo runs for the same scenario as used in the first example were performed. The states of targets are estimated by expectation-maximization (EM) algorithm. In addition, the maximum possible number of targets is limited and known for MPF, but unlimited and unknown for the PHD filter and the shrinkage-PHDF.

To measure the estimation quantitatively, the estimation errors in terms of Monte Carlo averaged OSPA distance are shown in Figs. 6 and 7. For position estimation, the cut-off value c is five times the size of the range cell. For velocity estimation, the cut-off value c is two times the size of the velocity cell. It is shown that the shrinkage-PHDF estimates target position and velocity on all tracks with higher estimated precision than the MPF, because MPF requires a modeling setup to accommodate varying numbers of targets, which is very difficult to accomplish in reality.It also shows that the performance of shrinkage-PHDF is better than PHDF, which is discussed further below.

To compare the shrinkage-PHDF and the PHDF clearly, Figs. 8 and 9 show the comparison of the mean OSPA for varying SNRs. The SNR of the targets is set to known values between 6 and 13 dB. When the value is between 6 and 11 dB, the shrinkage-PHDF performs significantly better than the PHDF. When the value is large than 11 dB, they perform similarly. This finding illustrates the limitation of ordinary the PHDF for TBD: the PHD cannot sufficiently approximate the multitarget posterior probability when the SNR is low. The superiority of the shrinkage-PHDF takes advantage of the shrinkage operation is indicated, especially for low SNR scenarios. It is clear that for multitarget TBD, the shrinkage-PHDF performs better than the PHDF and MPF, especially when the SNR is low.

In this example, the targets with different SNRs (the first target is 8dB, the second target is 9dB and the last one is 10dB). At time step1, the first target appears as in the first example. At step 7,the second target spawns form the first one, and the velocity in the x-direction is -300m/s. The third target new-birth far from the previous two targets at the time step 13 with the intital target state ${\left[{\begin{array}[]{*{20}{c}}{89000}&{-250}&0&0\\ \end{array}}\right]^{T}}$. A comparison of the estimation errors of the position and the velocity of the three algorithms is shown in Fig. 10 and 11. In this scene, the targets vary not only the value of SNRs but also the average number of targets. Because the SNR of the first target is lower than which in section 6.2, the convergence velocity of the estimation first target is slower. After time step 7, comparing Figs. 6 and 10, Fig.7 and 11, the varying SNRs and the number of the targets did not significantly influence the performance of the PHDF and Shrinkage-PHDF. In this example, the performance of MPF is better than that in the first example. The reason is the intensity of the targets is the most important characteristic to differentiate different targets for MPF method. In this way, the SNRs of all targets should be known for MPF method.

The PHD filter and the shrinkage-PHD filter are compared in Fig. 12 and 13. The sense is similar with that in the section 6.2, but the SNRs of the targets are unknown and the SNR of the second target changes from 9dB to 10dB. It is obvious that the performance of both algorithms is similar to what was presented in Figs. 6 and 7. After the 11th time step, the performance is even better than which in the section 6.2, because the value of the SNR of the second target is increased. Therefore, whether the SNRs of multiple targets are known has not significant influence on the performance of PHDF and Shrinkage-PHDF. In addition, MPF from Ref8 does not work when the value of the SNRs of the targets are unknown.

Table 3 shows the computation resources needed for the PHDF, the shrinkage-PHDF and MPF. All simulations were run on a PC with an Intel (R) Core TM2 Duo CPU E4500 @ 2.20 GHz processor. The time listed is Monte Carlo-averaged, for a run in a scenario similar to that given in Section 6.2. The amount of calculation time required for the shrinkage-PHDF and the PHDF is much less than for MPF. The calculated amount time required for the shrinkage-PHDF is slightly greater than that required for the PHD filter.