Two implementations of marginal distribution Bayes filter for nonlinear Gaussian models

Int. J. Electron. Commun. (AEÜ) 69 (2015) 1297–1304

Contents lists available at ScienceDirect

International Journal of Electronics and Communications (AEÜ) journal homepage: www.elsevier.com/locate/aeue

Two implementations of marginal distribution Bayes filter for nonlinear Gaussian models Zong-xiang Liu ∗ , Li-juan Li, Wei-xin Xie, Liang-qun Li ATR Key Laboratory, College of Information and Engineering, Shenzhen University, Shenzhen 518060, China

a r t i c l e

i n f o

Article history: Received 23 November 2014 Accepted 13 May 2015 Keywords: Multi-target tracking Bayes filter Marginal distribution Existence probability Nonlinear multi-target models

a b s t r a c t The marginal distribution Bayes (MDB) filter is an efficient approach for tracking an unknown and timevarying number of targets in the presence of clutter, noise, data association uncertainty, and detection uncertainty. This filter propagates the marginal distributions and existence probabilities of each target in the filter recursion, and it admits a closed-form solution for a linear Gaussian multi-target model. However, this closed-form solution is not general enough to accommodate nonlinear multi-target models. In this paper, we propose two implementations of the MDB filter to accommodate nonlinear multi-target models. The first is the first-order Taylor approximation MDB (FTA-MDB) filter which is based on the linearization technique of nonlinear function, and the second is the unscented transform MDB (UT-MDB) filter which is based on the unscented transform technique. Simulation results demonstrate that the proposed implementations are better on multiple targets tracking than the UK-PHD filter. © 2015 Elsevier GmbH. All rights reserved.

1. Introduction Multi-target tracking (MTT) has become increasingly important in many applications that range from traffic monitoring to automated surveillance. The primary objective of MTT is to detect individual targets in the surveillance region of interest and estimate their states according to a sequence of noisy and cluttered measurements collected by sensors [1–4]. The most efficient technique for MTT is the multi-target Bayes filter, which propagates joint posterior distribution of the multi-target state [1,2]. However, such propagation is computationally intensive because of the high dimensionality of the multi-target state space [2,3]. With the use of the Bayesian framework to propagate the posterior intensity of multiple targets recursively, the probability hypothesis density (PHD) filter provides a numerically tractable solution to this problem [2,3]. Two numerical solutions, namely, sequential Monte Carlo (SMC) [5] and Gaussian mixtures (GM) [6], have been developed for the PHD filter. Several extensions of the PHD filter have also been proposed to improve its performance [7–12]. The efficiency, simplicity of implementation, and success of the PHD filter in avoiding the combinatorial problem that arises from data association has made it more appealing than other filters.

∗ Corresponding author. Tel.: +86 755 26732055; fax: +86 755 26732049. E-mail address: [email protected] (Z.-x. Liu). http://dx.doi.org/10.1016/j.aeue.2015.05.007 1434-8411/© 2015 Elsevier GmbH. All rights reserved.

However, the PHD filter cannot distinguish multiple distinct targets when they are closely spaced. To efficiently distinguish closely spaced targets according to a sequence of measurements, Liu and Xie [13] proposed the marginal distribution multi-target Bayes (MDB) filter. Instead of maintaining the joint state distribution, this filter jointly propagates the marginal distributions and existence probabilities of each target in the filter recursion. A closed-form solution of the proposed filter for linear and Gaussian models is also presented to deal with an unknown and variable number of targets. In addition, a new approach on Gaussian Bayes filtering to compute the marginal function of a global function has been proposed in Ref. [14]. Although the closed-form recursion of the MDB filter is efficient for linear Gaussian models, it is not general enough to accommodate nonlinear Gaussian models. In this paper, we propose two implementations of the MDB filter to accommodate nonlinear multi-target models. The first is based on the linearization technique, while the second is based on the unscented transform technique. Simulation results are presented to demonstrate the capability of the proposed methods. The remainder of this paper is organized as follows: Section 2 briefly introduces the marginal distribution multi-target Bayes filter. Section 3 proposes two implementations of the MDB filter to accommodate nonlinear Gaussian models using the linearization technique and unscented transform technique. Section 4 evaluates the performance of the proposed methods. Section 5 states the drawn conclusions.

1298

Z.-x. Liu et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 1297–1304

2. Marginal distribution multi-target Bayes filter

The MDB filter combines the M distributions to a distribution as follows:

Marginal distribution multi-target Bayes filter, which was proposed by Liu and Xie [13], propagates the marginal distributions and existence probabilities of each target in the filter recursion. In this filter, the marginal distributions are used to model the uncertainties of individual target states caused by the dynamic uncertainty and measurement uncertainty of the target, while the existence probabilities of individual targets are used to characterize the randomness of target appearance and disappearance. Let Nt−1 denote the number of targets at time step t − 1; let xk,t−1 , k = 1, . . ., Nt−1 denote the states of individual targets at time step t − 1; let fk ( xk,t−1 | y 1:t−1 ), k = 1, . . ., Nt−1 denote probability density distributions of individual targets at time step t − 1; let pk,t−1 , k = 1, . . ., Nt−1 denote existence probabilities of individual targets at time step t − 1; and let y 1:t−1 = { y 1 , . . ., y t−1 } denote all observations up to time step t − 1. The prediction step of the MDB filter is given by



=

k = 1, . . ., Nt−1

Nt|t−1

(2)

k=1

=



Nt−1 

fk (xk,t |y 1:t−1 )

k=1

N

∪ fk (xk,t )

k=1

=



pk,t|t−1

Nt−1  k N ∪ p

k=1

tributions fk (xk,t ), k = 1, . . ., N . In the update step, the MDB filter uses the Bayes rule to obtain the individual update distributions as i fk,t|t xk,t |y 1:t−1 , y i,t

=



pD,t pk,t|t−1 gk (y i,t |xk,t )fk (xk,t |y 1:t−1 )

c,t + pD,t

Nt|t−1 e=1

pe,t|t−1



ge (y i,t |xe,t )fe (xe,t |y 1:t−1 )dxe,t

i = 1, . . ., M

, (5)

where M is the number of observations at time step t; y t = {y1,t , . . ., yM,t } denote the all measurements at time step t; c,t is the clutter density; pD,t is the detection probability of the target; and gk ( y i,t | xk,t ) is the measurement likelihood. The existence probability of the   i xk,t |y 1:t−1 , y i,t is defined as its integral updated distribution fk,t|t



pik,t|t =



e=1

pe,t|t−1



ge (y i,t |xe,t )fe (xe,t |y 1:t−1 )dxe,t

where c is the normalized coefficient and is given by c =

M 





i fk,t|t xk,t |y 1:t−1 , y i,t dxk,t

i=1

=

y i,t ∈y t

pD,t pk,t|t−1





gk (y i,t |xk,t )fk (xk,t |y 1:t−1 )dxk,t

Nt|t−1

c,t + pD,t

pe,t|t−1



ge (y i,t |xe,t )fe (xe,t |y 1:t−1 )dxe,t



i fk,t|t xk,t |y 1:t−1 , y i,t dxk,t ,

i = 1, . . ., M

The existence probability of the combined distribution fk,t|t ( xk,t | y 1:t ) is defined as pk,t|t = 1 −

M





1 − pik,t|t ,

k = 1, . . ., Nt|t−1

(9)

i=1

Considering the case that several targets may not be detected at time step t, the posterior distributions of individual targets at time step t are given by



fk (xk,t |y 1:t ) =

  

 

pk,t|t fk,t|t xk,t |y 1:t + 1 − pD,t pk,t|t−1 fk (xk,t |y 1:t−1 ) pk,t|t + 1 − pD,t pk,t|t−1

k = 1, . . ., Nt|t−1

,

(10)

The existence probabilities of posterior distributions fk ( xk,t | y 1:t ), k = 1, . . ., Nt|t−1 are given by pk,t = 1 − (1 − pk,t|t )(1 − (1 − pD,t )pk,t|t−1 ),

k = 1, . . ., Nt|t−1 (11)

(4)

k=1

where pk , k = 1, . . ., N are the existence probabilities of new dis-



c,t + pD,t

(7)

(3)

k=1

where Nt|t−1 is the number of the extended prediction distributions, N is the number of new distributions, and Nt|t−1 = Nt−1 + N ; fk (xk,t ), k = 1, . . ., N are the distributions of new targets. The existence probabilities of the extended prediction distributions fk ( xk,t | y 1:t−1 ), k = 1, . . ., Nt|t−1 are given by pk,t|t−1

y i,t ∈y t

pD,t pk,t|t−1 gk (y i,t |xk,t )fk (xk,t |y 1:t−1 )

Nt|t−1

(8)

where pS,t is the survival probability of the target. In addition to existing targets, new targets may appear at time step t. To enable the MDB filter to track the new targets appearing at time step t, The prediction distributions fk ( xk,t | y 1:t−1 ), k = 1, . . ., Nt−1 are extended to include the distributions of new targets as

Nt|t−1

c

(1)

pk,t|t−1 = pS,t pk,t−1 ,



1

f (xk,t |xk,t−1 )fk (xk,t−1 |y 1:t−1 )dxk,t−1 ,

where fk ( xk,t | y 1:t−1 ), k = 1, . . ., Nt−1 are the individual prediction distributions, and f( xk,t | xk,t−1 ) denotes the Markov transition probability from state xk,t−1 at time step t − 1 to state xk,t at time step t. The existence probability pk,t|t−1 of the prediction distribution fk ( xk,t | y 1:t−1 ) is given by

fk (xk,t |y 1:t−1 )

M

e=1

k = 1, . . ., Nt−1



 1 i  fk,t|t xk,t |y 1:t−1 , y i,t c i=1



fk (xk,t |y 1:t−1 ) =



fk,t|t xk,t |y 1:t =

(6)

The prediction distribution fk ( xk,t | y 1:t−1 ) generates M updated   i xk,t |y 1:t−1 , y i,t , i = 1, . . ., M in the update step. distributions fk,t|t

3. Two implementations of the MDB filter for nonlinear Gaussian models An implementation of the MDB filter for linear Gaussian models has been presented in Ref. [13]. In this section, we propose two implementations of the MDB filter to accommodate nonlinear Gaussian models. These two implementations require the following assumptions: A1. Each target follows a nonlinear Gaussian dynamic model and the sensor has a nonlinear Gaussian measurement model, i.e., xk,t = ϕk (xk,t−1 ) + ut−1

(12)

y k,t = hk (xk,t ) + w t

(13)

where ϕk (·) and hk (·) are nonlinear state and observation functions, respectively; xk,t−1 and xk,t denote the state vectors at time steps t − 1 and t, respectively; and ut−1 and w t are zero mean Gaussian noises with covariance matrices Q t−1 and R t , respectively. A2. The new distributions at time step t are Gaussian distributions and are given by



N

fk (xk,t )

k=1

=



N

N(xk,t ; mk , P k )

k=1

(14)

Z.-x. Liu et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 1297–1304

where N is the number of new distributions; mk and P k are the known mean vector and covariance matrix of new distribution k, respectively. At the same time, we also assume that a known existence probability p is assigned to each new distribution, i.e., pk

= p ,

k = 1, . . ., N

(15)

Based on assumptions A1 and A2, we present two nonlinear implementations of the MDB filter using linearization and unscented transform techniques. 3.1. First-order Taylor approximation MDB (FTA-MDB) filter

3.1.1. Prediction step If the Gaussian distributions and existence probabilities of individual targets at time step t − 1 are as follows: fk (xk,t−1 |y 1:t−1 ) = N(xk,t−1 ; mk,t−1 , P k,t−1 ), pk,t−1 ,

k = 1, . . ., Nt−1

k = 1, . . ., Nt−1

where xk,t−1 is the state vector, mk,t−1 is the mean of Gaussian distribution k, and P k,t−1 is the corresponding covariance. Then by using first-order approximations wherever nonlinearities are encountered, the Gaussian approximations of individual prediction distributions can be given by k = 1, . . ., Nt−1

mk,t|t−1 = ϕk (mk,t−1 )

(19)

P k,t|t−1 = F k,t−1 P k,t−1 F Tk,t−1 + Q t−1

(20)

F k,t−1 =



∂ϕk (x)  ∂x x=m

(21)

k = 1, . . ., Nt−1

(22)

Nt|t−1

N(xk,t ; mk,t|t−1 , P k,t|t−1 ) =



k=1

Nt−1 

N(xk,t ; mk,t|t−1 , P k,t|t−1 )

pk,t|t−1

Nt|t−1 k=1

=



k=1

 N

∪ N(xk,t ; mk , P k )

k=1

pk,t|t−1

Nt−1  k N k=1

∪ p

(27)

H k,t =



∂hk (x)  ∂x x=m

(28) k,t|t−1

Ak = P k,t|t−1 H Tk,t (H k,t P k,t|t−1 H Tk,t + R t )

fk (xk,t |y 1:t−1 ) = N(xk,t ; mk,t|t−1 , P k,t|t−1 ),

k = 1, . . ., Nt|t−1

(25)

and that the existence probabilities of prediction distributions are given by k = 1, . . ., Nt|t−1

−1

(29)

mik,t|t = mk,t|t−1 + Ak · (y i,t − hk (mk,t|t−1 ))

(30)

P ik,t|t = (I − Ak · H k,t )P k,t|t−1

(31)

c=

y i,t ∈y t

pD,t pk,t|t−1 N(y i,t ; hk (mk,t|t−1 ), H k,t P k,t|t−1 H Tk,t c,t + pD,t

Nt|t−1 e=1

+ Rt )

(32)

pe,t|t−1 N(y i,t ; he (me,t|t−1 ), H e,t P e,t|t−1 H Te,t + R t )

The existence probability of the combined distribution fk,t|t ( xk,t | y 1:t ) is given by M



pk,t|t = 1 −



1 − pik,t|t ,

k = 1, . . ., Nt|t−1

(33)

i=1

where pik,t|t =

pD,t pk,t|t−1 N(y i,t ; hk (mk,t|t−1 ), H k,t P k,t|t−1 H Tk,t + R t ) c,t + pD,t

Nt|t−1 e=1

,

pe,t|t−1 N(y i,t ; he (me,t|t−1 ), H e,t P e,t|t−1 H Te,t + R t )

(34)

i = 1, . . ., M

Considering missing detections of targets, the posterior distributions of individual targets at time step t are given by

fk (xk,t |y 1:t ) =

pk,t|t fk,t|t xk,t |y 1:t



+ (1 − pD,t ) pk,t|t−1 N(xk,t ; mk,t|t−1 , P k,t|t−1 )

pk,t|t + (1 − pD,t ) pk,t|t−1

,

(35)

k = 1, . . ., Nt|t−1

The existence probabilities of posterior distributions fk ( xk,t | y 1:t ), k = 1, . . ., Nt|t−1 are given by











1 − 1 − pD,t pk,t|t−1 ,

k = 1, . . ., Nt|t−1 (36)

(24)

k=1

pe,t|t−1 N(y i,t ; he (me,t|t−1 ), H e,t P e,t|t−1 H e,t + R t )

e=1

(23)

3.1.2. Update step Assuming that the prediction distributions at time step t are given by

pk,t|t−1 ,

c,t + pD,t

y i,t ∈y t

pk,t = 1 − 1 − pk,t|t

The existence probabilities of extended prediction distributions are given by



pD,t pk,t|t−1 N(yi,t ; hk (mk,t|t−1 ), H k,t P k,t|t−1 H Tk,t + Rt )N(xk,t ; mik,t|t , P ik,t|t ) Nt|t−1 T



Extending the prediction distributions in (18) to include the new distributions in (14), the extended prediction distributions are as follows:



1 c

k,t−1

The existence probabilities of prediction distributions N( xk,t ; mk,t|t−1 , P k,t|t−1 ), k = 1, . . ., Nt−1 are given by pk,t|t−1 = pS,t pt−1 ,

=



(18)

where mk,t|t−1 and P k,t|t−1 are the mean vector and covariance matrix of prediction distribution k, respectively, and are given by

where



fk,t|t xk,t |y 1:t

(16) (17)

fk (xk,t |y 1:t−1 ) = N(xk,t ; mk,t|t−1 , P k,t|t−1 ),

Approximating the nonlinear measurement models by using first-order approximations wherever nonlinearities are encountered, the combined distribution fk,t|t ( xk,t | y 1:t ) in (7) is as follows:

where

Similar to the extended Kalman filter (EKF), a nonlinear approximation to the MDB recursion is proposed based on applying first-order Taylor approximations of nonlinear functions ϕk (·) and hk (·). The FTA-MDB filter consists of the following steps.

1299

(26)

Similar to the recursion of MDB filter for linear Gaussian models in Ref. [13], a recursion of FTA-MDB filter for nonlinear Gaussian models comprises, in addition to the above prediction and update steps, filtering distribution approximation step and multi-target state extraction step in Ref. [13]. 3.2. Unscented transform MDB (UT-MDB) filter This section considers an implementation of the MDB filter to accommodate nonlinear Gaussian models by applying the unscented transform (UT) technique. Similar to the unscented Kalman filter (UKF), we use UT technique to obtain a nonlinear approximation to the MDB recursion in this section. The prediction step and update step of the UT-MDB filter are as follows:

1300

Z.-x. Liu et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 1297–1304

3.2.1. Prediction step If the Gaussian distributions and existence probabilities of individual targets at time step t − 1 are as follows:

Applying the UT technique wherever nonlinearities are encountered, the combined distribution fk,t|t ( xk,t | y 1:t ) in (7) can be given by

fk (xk,t−1 |y 1:t−1 ) = N(xk,t−1 ; mk,t−1 , P k,t−1 ),

fk,t|t xk,t |y 1:t

k = 1, . . ., Nt−1

k = 1, . . ., Nt−1

pk,t−1 ,

(37)



N(xk,t ; mk,t|t−1 , P k,t|t−1 ),

k = 1, . . ., Nt−1

y k,t|t−1 =

j=0



ϕk (xk,j ) − mk,t|t−1

T

(40)



xk,j = mk,t−1 +

xk,j+n = mk,t−1 −

(n + )P k,t−1



 j

(n + )P k,t−1

j = 1, . . ., n

 j

wj+n

(42)

j = 1, . . ., n

k = 1, . . ., Nt−1

(43)

Extending the prediction distributions in (39) to include the new distributions in (14), the extended prediction distributions are as follows:

Nt|t−1

N(xk,t ; mk,t|t−1 , P k,t|t−1 )

k=1

Nt−1 

N(xk,t ; mk,t|t−1 , P k,t|t−1 )

k=1

N

∪ N(xk,t ; mk , P k )

k=1

(44)

The existence probabilities of the extended prediction distributions are given by



pk,t|t−1

Nt|t−1 k=1

=



pk,t|t−1

Nt−1  k N k=1

∪ p

(45)

k=1

3.2.2. Update step Assuming that the prediction distributions at time step t are given by fk (xk,t |y 1:t−1 ) = N(xk,t ; mk,t|t−1 , P k,t|t−1 ),

k = 1, . . ., Nt|t−1

(46)

and that the existence probabilities of the prediction distributions are given by pk,t|t−1 ,

k = 1, . . ., Nt|t−1

+ Rt



hk (xk,j ) − y k,t|t−1

T

(49)

P ik,t|t = P k,t|t−1 − Ak P k,yy ATk

y i,t ∈y t

pD,t pk,t|t−1 N(y i,t ; y k,t|t−1 , P k,yy )



(50)

Nt|t−1

c,t + pD,t

j = 1, . . ., n

j = 1, . . ., n

pk,t|t−1 = pS,t pt−1 ,

T

pe,t|t−1 N(y i,t ; y e,t|t−1 , P e,yy )

e=1

The existence probabilities for the prediction distributions N( xk,t ; mk,t|t−1 , P k,t|t−1 ), k = 1, . . ., Nt−1 are given by





wj xk,j − mk,t|t−1

The used sigma points in (49) include 2 points of xk,0 , n points of xk,j , and n points of xk,j+n with xk,0 = mk,t|t−1



xk,j+n = mk,t|t−1 −

1 = , 2(n + )

=

2n

xk,j = mk,t|t−1 +

1 , 2(n + )

hk (xk,j ) − y k,t|t−1

mik,t|t = mk,t|t−1 + Ak (y i,t − y k,t|t−1 )

(41)

 w0 = n+



P k,xy =

c=

where n denotes the dimension of the state vector,  is the scaling factor, and (·)j denotes the column vector of the matrix. The weights of sigma points are given by

wj =



j=0

The used sigma points in (40) include 2 points of xk,0 , n points of xk,j , and n points of xk,j+n with





wj hk (xk,j ) − y k,t|t−1

Ak = P k,xy P −1 k,yy

+ Q t−1

j=0

xk,0 = mk,t−1

wj hk (xk,j )

j=0

wj ϕk (xk,j )

wj ϕk (xk,j ) − mk,t|t−1

2n

j=0 2n

2n

P k,t|t−1 =

(48)

pe,t|t−1 N(y i,t ; y e,t|t−1 , P e,yy )

e=1

where

P k,yy =

2n 

c,t + pD,t

y i,t ∈y t

(39)

where

pD,t pk,t|t−1 N(yi,t ; yk,t|t−1 , P k,yy )N(xk,t ; mik,t|t , P ik,t|t ) Nt|t−1

1 c

=

(38)

Then the Gaussian approximations of the individual prediction distributions can be given by

mk,t|t−1 =



(47)

(n + )P k,t|t−1



 j

(n + )P k,t|t−1

j = 1, . . ., n

 j

(51)

j = 1, . . ., n

where n denotes the dimension of the state vector, (·)j denotes the column vector of a matrix, and the weights of sigma points are given by (42). The existence probability of the combined distribution fk,t|t ( xk,t | y 1:t ) is given by pk,t|t = 1 −

M





1 − pik,t|t ,

k = 1, . . ., Nt|t−1

(52)

i=1

where pik,t|t =

pD,t pk,t|t−1 N(y i,t ; y k,t|t−1 , P k,yy ) c,t + pD,t

Nt|t−1 e=1

pe,t|t−1 N(y i,t ; y e,t|t−1 , P e,yy )

,

i = 1, . . ., M

(53)

Considering missing detections of targets, the posterior distributions of individual targets at time step t are given by fk (xk,t |y 1:t ) =





pk,t|t fk,t|t xk,t |y 1:t + (1 − pD,t ) pk,t|t−1 N(xk,t ; mk,t|t−1 , P k,t|t−1 ) pk,t|t + (1 − pD,t ) pk,t|t−1

k = 1, . . ., Nt|t−1

,

(54)

The existence probabilities of the posterior distributions fk ( xk,t | y 1:t ), k = 1, . . ., Nt|t−1 are given by



pk,t = 1 − 1 − pk,t|t









1 − 1 − pD,t pk,t|t−1 ,

k = 1, . . ., Nt|t−1 (55)

Z.-x. Liu et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 1297–1304

400

300

200

100

y (m)

Similar to the recursion of FTA-MDB filter for nonlinear Gaussian models, a recursion of UT-MDB filter for nonlinear Gaussian models comprises, in addition to the above prediction and update steps, filtering distribution approximation step and multi-target state extraction step in Ref. [13]. Note that the FTA-MDB filter requires the computation of Jacobians, and is only applicable to the differentiable state and measurement models, whereas the UT-MDB filter completely avoids the differentiation requirement and is even applicable to models with discontinuities.

0

4. Simulation results

-100

In order to demonstrate the performance of the proposed FTAMDB and UT-MDB filters, we select the UK-PHD filter as a contesting comparison object. The selected filter that was proposed by Vo in Ref. [6] is a state-of-the-art filter. In the simulation, we consider a 2D nonlinear tracking scenario where targets move at a constant turn model with varying turn rate in the surveillance region [−400 m, 400 m] × [−400 m, 400 m]. The state of each target consists of its position, velocity, and turn rate, and is represented T by xk,t = [xk,t , x˙ k,t , yk,t , y˙ k,t , ωk,t ] . The used nonlinear dynamic model is given by (12) with

-200

⎡

sin(ωk,t T )

1

0

ωk,t ⎢ ⎢ ⎢0 cos(ωk,t T ) 0 ⎢ ⎢ 1 − cos(ωk,t T ) ϕk (xk,t−1 ) = ⎢ 1 ⎢0 ωk,t ⎢ ⎢ sin(ωk,t T ) 0 ⎣0 0

0

−

1 − cos(ωk,t T ) ωk,t

0

−sin(ωk,t T )

0⎥

sin(ωk,t T ) ωk,t cos(ωk,t T )

0

0

⎡ T 4  2 /4 T 3  2 /2 0 0 v v ⎢ T 3  2 /2 T 2  2 0 0 v v ⎢ ⎢ Q t−1 = ⎢ 0 T 4 v2 /4 T 3 v2 /2 ⎢ 0 ⎢ ⎣ 0 0 T 3 v2 /2 T 2 v2 0

0

0

⎤

0

⎥ ⎥ ⎥ ⎥ ⎥x 0 ⎥ k,t−1 ⎥ ⎥ 0⎦ 1

⎤

0



⎥ ⎥ ⎥ 0 ⎥ ⎥ ⎥ 0 ⎦ 0

⎢ ⎢

⎛

hk (xk,t ) = ⎢ ⎢

⎣ arccos ⎝ 

xk,t − xs (xk,t − xs )2 + (yk,t − ys )2



⎞⎥ 2 ⎥ ⎥ , Rt = r  2 ⎠⎥ ⎦

dp,c (X, Y ) =

1 n

min 

∈

n

m i=1

1/p dc (xi , y (i) )p + c p (n − m)

T4

T6 T5

-400 -400

-300

-200

-100

0

100

200

300

400

x (m) Fig. 1. True target trajectories.



where dc (x, y) = min(c, d(x, y)) is the cutoff distance, n represents the set {1, . . ., n}, and 1≤ p < ∞ is the OSPA metric parameter. If m > n, the definition is dp,c ( X, Y ) = dp,c ( Y , X). In our simulations, the parameters p and c are set to p = 2 and c = 10 m. There are six targets in the simulations. Targets 1 and 2 appear at t = 1 s, and continue to exist for the remaining times. Targets 3 and 4 appear at t = 3 s, and continue to exist for the remaining times. Targets 5 and 6 appear at t = 5 s, and disappear at t = 65 s and t = 60 s, respectively. All targets change their turn rates from ωk,t−1 to ωk,t = − ωk,t−1 at t = 26 s. The initial states of the individual targets are as follows:

10 (m s−1 ), −0.045 (rad s−1 )]



where  r = 2 m,  = 0.003 rad, and [xs , ys is the position vector of the sensor. To assess the performances of the FTA-MDB, UT-MDB and UKPHD filters, we use the optimal subpattern assignment (OSPA) distance, which was recently developed by Schuhmacher et al. in [15], as the measure. For two arbitrary finite subsets X = {x1 , · · · , xm } and Y = {y1 , · · · , yn }, if m ≤ n, the OSPA distance is defined as

 

T3

T

x2,1 = [−340 (m), 0 (m s−1 ), −200 (m),

2 T 2 ω

]T

Initial position Target trajectory

-300

10 (m s−1 ), −0.045 (rad s−1 )]

⎤ (xk,t − xs )2 + (yk,t − ys )2

T1 T2

x1,1 = [−380 (m), 0 (ms−1 ), −200 (m),

where v and  ω are standard deviations of process noises; and T = 1 s is the sampling period. A sensor located at [0 m, −100 m]T observes the targets in the surveillance region. The nonlinear observation model is given by (13) with

⎡

1301

T

x3,3 = [−350 (m), 8 (m s−1 ), −300 (m), 10 (m s−1 ), −0.01 (rad s−1 )]

T

x4,3 = [−300 (m), 8 (m s−1 ), −300 (m), 10 (m s−1 ), 0.01 (rad s−1 )]

T

x5,5 = [−100 (m), 10 (m s−1 ), −340 (m), 10 (m s−1 ), 0.035 (rad s−1 )]

T

x6,5 = [0 (m), 10 (m s−1 ), −300 (m), 10 (m s−1 ), 0.035 (rad s−1 )]

T

The true trajectories of the six targets are shown in Fig. 1. The measurements are generated with the survival probability pS,t = 1.0, the detection probability pD,t = 0.98, the clutter density c,t = 0.0071 m−1 rad−1 (20 clutter returns over the surveillance region per scan period), the standard deviations of the process

1302

Z.-x. Liu et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 1297–1304

7

400

300

FTA-MDB filter UT-MDB filter UK-PHD filter

6

200

OSPA Distance (m)

5

y (m)

100

0

-100

4

3

-200

2

-300

1

-400 -400

-300

-200

-100

0

100

200

300

400

x (m)

0 0

N(xk,t ; mk , P k ),

k = 1, . . ., 6

where m1 = [−380 (m), 0 (m s−1 ), −200 (m), 0 (m s−1 ), 0 (rad s−1 )]

T

20

30

40

50

60

70

t (s)

Fig. 2. Measurement data.

noises v = 1 m s−2 and  ω = 0.001 rad s−2 , and the standard deviations of the observation noises  r = 2 m and  = 0.003 rad. Fig. 2 shows the measurement data for an experiment. A 1D view of the sensor measurements along with the true target trajectories is shown in Fig. 3. To deal with the simulation data, we set the corresponding parameters of the FTA-MDB, UT-MDB and UK-PHD filters to pS,t = 1.0, pD,t = 0.98, c,t = 0.0071 m−1 rad−1 , v = 1 m s−2 ,  ω = 0.1 rad s−2 ,  r = 2 m,  = 0.003 rad,  = 0,  = 10−3 and U = 4. The new state distributions at time step t are given by

10

Fig. 4. OSPA distances of the three filters.

m2 = [−340 (m), 0 (m s−1 ), −200 (m), 0 (m s−1 ), 0 (rad s−1 )]

T

m3 = [−350 (m), 0 (m s−1 ), −300 (m), 0 (m s−1 ), 0 (rad s−1 )]

T

m4 = [−300 (m), 0 (m s−1 ),

− 300 (m),

0 (m s−1 ), 0 (rad s−1 )]T

m5 = [−100 (m), 0 (m s−1 ), −340 (m), 0 (m s−1 ), 0 (rad s−1 )]T

m6 = [0 (m), 0 (m s−1 ), −300 (m), 0 (m s−1 ), 0 (rad s−1 )]

T 2

and P k = (diag([ 5 20 5 20 0.1])) , k = 1, . . ., 6. The existence probabilities of the six new distributions are set to pk = 0.1,

Fig. 3. 1D view of true tracks and measurement data.

k = 1, . . ., 6

We use FTA-MDB, UT-MDB, and UK-PHD filters to deal with the measurement data, respectively. Fig. 4 shows the OSPA distances given by these three filters for a Monte Carlo run, and Fig. 5 shows the average OSPA distances of these three filters for 100 Monte Carlo runs. Figs. 4–5 show that these three filters can detect multiple targets and estimate their states in the presence of clutter, noise, data association uncertainty, and detection uncertainty. As shown in Fig. 4, the peaks in the plot of the OSPA distance indicate that the filter may provide either excessive state estimates or insufficient state estimates and is penalized with the cutoff distance. According to the average OSPA distance shown in Fig. 5, the proposed FTA-MDB

Z.-x. Liu et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 1297–1304

1303

Table 1 Average OSPA distance (m) for 100 Monte Carlo runs. Filter

v = 0.1 m s−2  ω = 0.001 rad s−2

v = 1 m s−2  ω = 0.001 rad s−2

v = 1 m s−2  ω = 0.1 rad s−2

v = 3 m s−2  ω = 0.001 rad s−2

v = 3 m s−2  ω = 0.1 rad s−2

FTA-MDB UT-MDB UK-PHD

2.5414 2.5615 2.6748

2.5378 2.5550 2.6778

2.5469 2.5537 2.6767

2.5098 2.5304 2.6590

2.5462 2.5606 2.6909

Table 2 Average running time (s) for a Monte Carlo run. Filter

v = 0.1 m s−2  ω = 0.001 rad s−2

v = 1 m s−2  ω = 0.001 rad s−2

v = 1 m s−2  ω = 0.1 rad s−2

v = 3 m s−2  ω = 0.001 rad s−2

v = 3 m s−2  ω = 0.1 rad s−2

FTA-MDB UT-MDB UK-PHD

2.9939 4.4933 4.9423

2.8117 4.2242 4.5778

3.0387 4.5445 5.0630

2.8095 4.2850 4.6069

2.7467 4.0853 4.4208

and UT-MDB filters provide more accurate state estimates than the UK-PHD filter. The main reason for this may be that the UK-PHD filter propagates the first-order moment of the joint multi-target posterior density, whereas the FTA-MDB and UT-MDB filters propagate marginal density distributions of individual targets. Thus, the target state of the UK-PHD filter is more easily affected by clutter and other target states than the target states of the FTA-MDB and UT-MDB filters. Fig. 5 also shows that the FTA-MDB filter is slightly better than the UT-MDB filter on multiple targets tracking because its average OSPA distance is slightly smaller than that of the UTMDB filter at most times. Due to the fact that these three filters need a few time steps to confirm the appearance of new targets, large OSPA distances emerge between t = 1 s and t = 8 s in the plots of these three filters. 4.1. Effect of process noise To show the effect of the process noise on the tracking performance of these three filters, we use different standard deviations of process noises to generate measurement data, and run 100 Monte Carlo trials for each filter to obtain the average OSPA distance. The average OSPA distance shown in Table 1 reveals the following results: (1) these three filters are insensitive to the standard deviation of the process noise; (2) the average OSPA distances of

FTA-MDB filter UT-MDB filter UK-PHD filter

Average OSPA Distance (m)

5 4.5 4 3.5

Similar to the PHD filter, the computational complexity of the MDB filter is O(KM), where K is the number of targets in the scene, and M is the number of observations in the current measurement set [13]. Table 2 shows the average performing times of these three filters for different standard deviations of process noises. As shown in Table 2, the average performing times of the FTAMDB and UT-MDB filters are smaller than the average performing time of the UK-PHD filter because the FTA-MDB and UT-MDB filters propagate less Gaussian terms than the UK-PHD filter in the filter recursion. Table 2 also shows that the average performing time of the FTA-MDB filter is less than that of the UT-MDB filter. This is due to the fact that the UT-MDB filter uses several sigma points to handle a Gaussian term in the filter recursion, whereas the FTAMDB filter uses one point to deal with it. The UT approximation of the MDB filter requires more computation than its first-order Taylor approximation. 5. Conclusions

Acknowledgement This work was supported by National Natural Science Foundation of China (Nos. 61271107 and 61301074), Shenzhen Basic Research Project (No. JCYJ20140418095735618) and Defense Advance Research Fund Project (91400C800501140C80340).

3 2.5

2 0

4.2. Computational complexity

In this paper, we propose two implementations of the MDB filter to accommodate nonlinear multi-target models by applying linearization and unscented transform techniques. These two implementations propagate the marginal distributions and existence probabilities of each target in the filter recursion, and are applicable to tracking multiple targets in the presence of clutter, noise, data association uncertainty, and detection uncertainty. Based on the simulation data, we compare these two implementations with the UK-PHD filter according to the OSPA distance. The simulation results show that both the FTA-MDB and UT-MDB filters track multiple targets better as well as require less computation than the UK-PHD filter in the presence of clutter, noise and a variable number of targets.

6 5.5

the FTA-MDB and UT-MDB filters are consistently smaller than the average OSPA distance of the UK-PHD filter.

10

20

30

40

50

60

70

References

t (s) Fig. 5. Average OSPA distances of the three filters.

[1] Vermaak J, Godsill SJ, Perez P. Monte Carlo filtering for multi-target tracking and data association. IEEE Trans Aerosp Electron Syst 2005;41(1):309–31.

1304

Z.-x. Liu et al. / Int. J. Electron. Commun. (AEÜ) 69 (2015) 1297–1304

[2] Mahler R. Statistical multisource-multitarget information fusion. Norwood, MA: Artech House; 2007. [3] Mahler R. Multitarget Bayes filtering via first-order multitarget moments. IEEE Trans Aerosp Electron Syst 2003;39(4):1152–78. [4] Li L-Q, Xie W-X. Bearings-only maneuvering target tracking based on fuzzy clustering in a cluttered environment. AEU – Int J Electron Commun 2014;68(2):130–7. [5] Vo BN, Doucet SS. Sequential Monte Carlo methods for multi-target filtering with random finite sets. IEEE Trans Aerosp Electron Syst 2005;41(4):1224–45. [6] Vo BN, Ma WK. The Gaussian mixture probability hypothesis density filter. IEEE Trans Signal Process 2006;54(11):4091–104. [7] Ouyang C, Ji HB, Guo ZQ. Extensions of the SMC-PHD for jump Markov systems. Signal Process 2012;92:1422–30. [8] Pasha SA, Vo BN, Tuan HD, Ma WK. A Gaussian mixture PHD filter for jump Markov system models. IEEE Trans Aerosp Electron Syst 2009;45(3):919–36.

[9] Liu ZX, Xie WX, Wang P, Yu Y. A sequential GM-based PHD filter for a linear Gaussian system. Sci China Inf Sci 2013;56(10), 102302(1–10). [10] Wang Y, Jing ZL, Hu SQ, Wu JJ. Detection-guided multi-target Bayesian filter. Signal Process 2012;92:564–74. [11] Liu ZX, Xie WX, Wang P, Yu Y. A Gaussian mixture PHD filter with the capability of information hold. Acta Electron Sin 2013;41(8):1603–8 [in Chinese]. [12] Vo BT, Vo BN, Cantoni A. Analytic implementations of the cardinalized probability hypothesis density filter. IEEE Trans Signal Process 2007;55(7):3553–67. [13] Liu ZX, Xie WX. Multi-target Bayesian filter for propagating marginal distribution. Signal Process 2014;105:328–37. [14] Li B, Wu N, Wang H, Tseng PH, Kuang JM. Gaussian message passing-based cooperative localization on factor graph in wireless networks. Signal Process 2015;111:1–12. [15] Schuhmacher D, Vo BT, Vo BN. A consistent metric for performance evaluation of multi-object filters. IEEE Trans Signal Process 2008;56(8):3447–57.