Gaussian mixture PHD filter for multi-sensor multi-target tracking with registration errors

Signal Processing 93 (2013) 86–99

Contents lists available at SciVerse ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Gaussian mixture PHD filter for multi-sensor multi-target tracking with registration errors Wenling Li a,n, Yingmin Jia a,b, Junping Du c, Fashan Yu d a

The Seventh Research Division and the Department of Systems and Control, Beihang University (BUAA), Beijing 100191, China Key Laboratory of Mathematics, Informatics and Behavioral Semantics (LMIB), Ministry of Education, SMSS, Beihang University (BUAA), Beijing 100191, China Beijing Key Laboratory of Intelligent Telecommunications Software and Multimedia, School of Computer Science and Technology, Beijing University of Posts and Telecommunications, Beijing 100876, China d School of Electrical Engineering and Automation, Henan Polytechnic University, Jiaozuo 454000, Henan, China b c

a r t i c l e in f o

abstract

Article history: Received 14 February 2012 Received in revised form 11 June 2012 Accepted 14 June 2012 Available online 11 July 2012

This paper studies the problem of multi-sensor multi-target tracking with registration errors in the formulation of random finite sets. The probability hypothesis density (PHD) recursion is applied by introducing the dynamics of the translational measurement bias into the associated intensity functions. Under the linear Gaussian assumptions on the bias dynamics, the Gaussian mixture implementation is used to give closed-form expressions. As the target state and the translational measurement bias are coupled through the likelihood in the update step, a two-stage Kalman filter is adopted to approximate the tractable form, which leads to a substantial reduction in computational complexity. Two numerical examples are provided to verify the effectiveness of the proposed filter. & 2012 Elsevier B.V. All rights reserved.

Keywords: Multi-target tracking Probability hypothesis density filter Registration errors

1. Introduction Registration errors compensation has been an important issue in multi-sensor data fusion systems regardless of the sensor measurements are processed in the centralized or distributed fashion. There are many kinds of sensor biases such as the translational bias in the state space, the rotational bias in the state space, the translational bias in the measurement space, the rotational bias in the measurement space and both translational and rotational biases [1]. The estimation of unknown translational measurement biases has received great attention [2], and it is this problem that we address in this paper. Note that it is vital to estimate these measurement biases

n

Corresponding author. Tel.: þ 86 10 8233 8683; fax: þ86 10 8231 6100. E-mail addresses: [email protected] (W. Li), [email protected] (Y. Jia), [email protected] (J. Du), [email protected] (F. Yu). 0165-1684/$ - see front matter & 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.sigpro.2012.06.030

as accurately as possible so that the multi-sensor measurements can be referenced to a common tracking coordinate frame [3]. To resolve this problem, many approaches have been proposed such as the least squares [4], the maximum likelihood [5] and the Kalman filtering [6], of which the Kalman filtering has been an attractive method for its efficiency. By stacking the sensor biases and the target states in a single vector, an augmented state Kalman filter (ASKF) can be applied to derive the bias estimates. However, the implementation of the ASKF might not be computationally feasible and numerical problems may arise especially for ill-conditioned systems [3]. To alleviate this disadvantage, for linear dynamics and measurement models, Friedland in [7] proposed a two-stage estimator by decoupling the bias estimates from target states, and it has been shown in [8] that this estimator is equivalent to the ASKF solution when a particular relationship between the initial parameters of two filters is satisfied. Similar ideas have been adopted in the literature [9–13] to develop various two-stage estimators. It should be mentioned that

W. Li et al. / Signal Processing 93 (2013) 86–99

the problem of measurement source uncertainty is not addressed in most of the existing approaches, which is often encountered in multi-target tracking. Although the data association techniques (e.g., the joint probabilistic data association (JPDA) and the multiple hypothesis tracking (MHT)) can be incorporated, the association results might be bad because of the effect of sensor biases. Recently, the finite set statistics (FISST) theory has been used to tackle the multi-target tracking problem which avoids the data association [14]. In the framework of FISST, the target states and measurements are modeled as two different random finite sets (RFSs), and as a consequence, the problem of tracking an unknown and time-varying number of targets in the clutter environment can be addressed in a natural manner. Moreover, the multi-target tracking can be formulated in a rigorous Bayesian framework by constructing the multi-target transition density and multi-target likelihood function. However, the optimal multi-target Bayes filter is generally intractable due to the existence of multiple set integrals and the combinatorial nature of the multi-target densities. To alleviate this intractability, the probability hypothesis density (PHD) filter has been proposed as a first order moment approximation to the multi-target posterior density [15]. It should be pointed out that the PHD recursion still requires solving multi-dimensional integrals. There are mainly two approaches to implement the PHD recursion, including the sequential Monte Carlo (SMC) [16,17] and the Gaussian mixture (GM) [18]. In the SMC-PHD filter, a large number of particles are used to approximate the multi-dimensional integrals, and therefore the main drawback is the high computational burden. In addition, some clustering techniques are required to extract the target state estimates, which might be often unreliable. To overcome these disadvantages, the GM-PHD filter was developed for linear Gaussian target dynamics and Gaussian birth model, in which the weights, means and covariance matrices are propagated analytically by the Kalman filter (KF). Specially, the nonlinear Kalman filter counterparts can be directly employed to deal with nonlinear target dynamics and measurement models. The convergence properties of two implementations were analyzed in [17,19]. As shown in [19], the GM-PHD filter can approximate the true PHD filter to any desired degree of accuracy under the linear Gaussian assumption of the dynamic model. Similar results have been extended to handle jump Markov models for tracking multiple maneuvering targets [20–23]. The particle and Gaussian mixture techniques have also been used to derive PHD smoothers [24–31]. In [32], the GM-PHD filter is extended to multi-sensor tracking system and the target state estimates are obtained sequentially at each sensor. However, the sensor registration errors are neglected. In [33], the problem of multi-sensor mis-registration with arbitrary biases has been addressed and a simplified version has been investigated in [34]. In [35], the translational measurement mis-registration is investigated in the PHD recursion and the multi-sensor SMC implementation is developed. Nevertheless, the proposed filter inherits the disadvantages of the SMC such as the high computational cost and unreliable clustering.

87

In this paper, we attempt to apply the GM-PHD filter to address the problem of multi-sensor multi-target tracking with registration errors. For clarity, we consider the linear target dynamics and measurement models in this work. Under the linear Gaussian assumptions on the translational measurement bias dynamics, the PHD recursion is applied by introducing the bias dynamics into the associated intensity functions. Then the Gaussian mixture implementation is used to give closed-form expressions. As the target state and the translational measurement biases are coupled through the likelihood in the update step, a suboptimal two-stage Kalman filter is adopted to jointly estimate the target state and the biases, which leads to a substantial reduction in computational complexity. Simulation results show that the proposed filter performs better than the standard GM-PHD filter without incorporating registration errors. The rest of this paper is organized as follows. The PHD recursion for multi-sensor multi-target tracking with registration errors is given in Section 2. The Gaussian mixture implementation to the PHD recursion is presented in Section 3. In Section 4, two numerical examples are provided to illustrate the effectiveness of the GM-PHD filter. Conclusion is drawn in Section Appendix A.

2. PHD recursion for multi-sensor multi-target tracking For multi-sensor tracking systems, the measurements received from each sensor should be transformed into a common coordinate system for merging, and the registration errors are often encountered. This is also referred to the estimation of unknown translational measurement biases. In this paper, we consider the following linear target dynamics and measurement models: xk ¼ F k1 xk1 þwk1

l

zlk ¼ Hlk xk þ bk þ vlk ,

ð1Þ

l ¼ 1; 2, . . . ,L

ð2Þ

where xk 2 Rn and zlk 2 Rm denote the target state and the lth sensor measurement, respectively. L is the number of l sensors. F k1 and Hk are the state transition matrix and l measurement matrix. wk1 and vk are zero-mean white Gaussian process noise and measurement noise with l l covariance matrices Q xk1 and Rk, respectively. bk is the translational measurement bias of the lth sensor. In the multi-target tracking scenario, targets might appear and disappear randomly. Thus, it is natural to model the set of target states as an RFS. Similarly, due to the presence of the clutter and the time-varying number of targets, the number of measurements received from each sensor might be random, the RFS can be also used to characterize the property of measurements. To be specific, an RFS X is a finite set valued random variable, which can be described by a discrete probability distribution and a family of joint probability densities. Specially, the discrete distribution characterizes the cardinality of X whereas an appropriate density characterizes the joint distribution of the elements in X. Assume that there are nk targets with

88

W. Li et al. / Signal Processing 93 (2013) 86–99 l

l

states xk,1 , . . . ,xk,nk in the surveillance region and mk measurements zlk,1 , . . . ,zlk,ml can be received from the lth k sensor at time step k, then the multi-target state and multi-target measurement can be represented as [14]

where pd is the detection probability of the lth sensor, l h ð9Þ is the single-target measurement likelihood of the lth sensor, klk ðÞ denotes the intensity of the clutter RFS of the lth sensor, and

X k 9fxk,1 , . . . ,xk,nk g  X

ð3Þ

F lk ðZ lk 9xk ,bk Þ ¼ 1pld

Z lk 9fzlk,1 , . . . ,zlk,ml g  Z l

ð4Þ

k

n

nk9k1 ðxk ,bk Þ ¼

nk19k1 ðxk1 ,bk1 9Z 1:k1 Þ dxk1 dbk1 þ gk ðxk ,bk Þ

ð5Þ

where ps is the target surviving probability. f x ð9Þ is the single-target transition density. f b ð9Þ is the bias transition density. bk9k1 ð9Þ and gk ðÞ denote the intensity of the spawned target RFS and the intensity of the spontaneously birth target RFS, respectively. Suppose that the biases are independent between sensors and the dynamics of each bias is Markovian, then f b ðbk 9bk1 Þ ¼

l

l

l

X zk 2Z k

½ps f x ðxk 9xk1 Þf b ðbk 9bk1 Þ þ bk9k1 ðxk ,bk 9xk1 ,bk1 Þ

L Y

þ

pld h ðzk 9xk ,bk Þ

R

0 0 klk ðzk Þ þ pld hl ðzk 9x0k ,b0k Þnk9k1 ðx0k ,b0k 9Z 1:L 1:k1 Þdxk dbk

p

l

where X  R and Z  R denote the state and the observation space, respectively. Then, the multi-target tracking can be formulated as a filtering process as follows: given the set of measurements Z 1:k ¼ fZ 11:k , . . . ,Z L1:k g from all sensors up to time k, the problem is to find the expectation of the posterior density function pðX k 9Z 1:k Þ. By using the finite set statistics theory, an optimal Bayesian recursion can be obtained in terms of the multitarget posterior density functions. However, this recursion involves multiple integrals and the multi-target posterior density functions are combinatorial, which makes it computationally intractable. To alleviate this intractability, inspired by the single-target tracking approaches, the propagation of the statistical moments associated with the posterior densities is adopted. The PHD recursion, which propagates the first order moment or the intensity function of multi-target random finite sets, provides a computationally cheaper alternative to the optimal multi-target Bayesian recursion [18]. Since the measurement biases should be estimated as well as the target states, similar to the derivation of the standard PHD recursion, we can obtain an extended version by T 1 augmenting the state vector ½xTk ,bk T , where bk ¼ ½ðbk ÞT , L T T . . . ,ðbk Þ  . Define the posterior intensity nk19k1 ðxk1 , bk1 9Z 1:k1 Þ at time k1, the predicted intensity nk9k1 ðxk ,bk 9Z 1:k1 Þ, and the posterior intensity nk9k ðxk ,bk 9Z 1:k Þ. For simplicity, ns9t ðxs ,bs 9Z 1:t Þ is shortly denoted by ns9t . The predicted intensity can be derived as Z

l

l

f b ðbk 9bk1 Þ

ð6Þ

ð9Þ For the multi-sensor PHD filter, both the iteratedcorrector PHD filter and the PHD filter defined by (8) have been investigated in [36]. Specially, the quality of the iterated PHD is affected by the sensor ordering at moderate and low probability of detection [37]. However, the improvement in performance is minor at a high probability of detection for the product multi-sensor PHD filter [36]. In addition, the approximation in (8) is valid only for relatively large numbers of targets [38]. Note that the PHD recursion given by (5)–(8) operates on the single-target state space and avoids the explicit problem of data association. However, it does not admit tractable solutions in general due to the multi-dimensional integrals. In the following section, the Gaussian mixture implementation is used to give closed-form expressions. 3. GM-PHD filter with registration errors Before we present the GM-PHD filter for multi-sensor multi-target tracking with registration errors, the following lemmas are adopted for further development [18]. Lemma 1. Given F, d, Q, m, and P of compatible dimensions and that Q and P are positive definite, then Z N ðx; F x þd,Q ÞN ðx; m,PÞ dx ¼ N ðx; Fm þ d,Q þFPF T Þ ð10Þ Lemma 2. Given H, R, m, and P of compatible dimensions and that R and P are positive definite, then N ðz; Hx,RÞN ðx; m,PÞ ¼ qðzÞN ðx; m,PÞ

ð11Þ

where qðzÞ ¼ N ðz; Hm,R þ HPHT Þ

ð12Þ

m ¼ m þKðzHmÞ

ð13Þ

P ¼ ðIKHÞP

ð14Þ

K ¼ PHT ðHPHT þ RÞ1

ð15Þ

l¼1 l

l

l

where f b ðbk 9bk1 Þ is the bias transition density of the lth l sensor. As in [35], the bias bk can be modeled as a first order Gauss–Markov process, i.e., l

l

l

l

f b ðbk 9bk1 Þ ¼ N ðbk ; bk1 ,Blk1 Þ

ð7Þ

After receiving the measurements from all sensors at time k, the updated intensity can be derived as in [35]

nk9k ¼ F Lk ðZ Lk 9xk ,bLk Þ    F 1k ðZ 1k 9xk ,b1k Þnk9k1

To derive Gaussian mixture implementations of the PHD recursion, the intensities of the birth and spawning random finite sets are assumed to be of the following forms:

gk ðxk ,bk Þ ¼

J g,k L Y X

l

j,l

wjg,k N ð½xk ; bk ; ½mjg,k ; bg,k ,½P jg,k ; Bj,l g,k Þ

j¼1l¼1

ð8Þ

ð16Þ

W. Li et al. / Signal Processing 93 (2013) 86–99

bk9k1 ðxk ,bk 9xk1 ,bk1 Þ ¼

Jb,k L Y X

89

Update step: Given that the predicted intensity can be represented as the form of l

i

l

i i,l wib,k N ð½xk ; bk ; ½F ix,k xk1 þ dx,k ; F i,l b,k bk1 ,½Q x,k ; Q b,k Þ

i¼1l¼1

Jk9k1

nk9k1 ¼

L X Y

l

j,l

wjk9k1 N ð½xk ; bk ; ½mjk9k1 ; bk9k1 ,½Pjk9k1 ; Bj,l Þ k9k1

j¼1l¼1

ð17Þ

ð30Þ

j

where J g,k , wjg,k , mjg,k , bg,k , Bjg,k and P jg,k are given parameters that determine the shape of the birth intensity. i J b,k , wib,k , F ix,k , dx,k F ib,k , Q ix,k and Q ib,k are given parameters that determine the shape of the spawning intensity. It is worth mentioning here that the intensity of the translational measurement biases are introduced into the intensities of the birth and spawning random finite sets, whereas this is not done in the standard PHD filter [18] and the extended PHD filter in [35]. The extended PHD recursion (5)–(8) can be carried out as follows (see Appendix for the detailed derivations). Prediction step: Given that the posterior intensity is a Gaussian mixture

nk19k1 ¼

J k1 Y L X

then the posterior intensity is updated sequentially for each sensor as X n1k9k ¼ ð1p1d Þnk9k1 þ n1d,k ðxk ; zk Þ ð31Þ zk 2Z 1k

Jk9k1

n1d,k ðxk ; zk Þ ¼

then the predicted intensity is also a Gaussian mixture with the form

nk9k1 ¼ ns,k9k1 þ nb,k9k1 þ gk ðxk ,bk Þ

l

j,l

ðzk Þ ¼ wj,1 k

p1d wjk9k1 qj,1 ðzk Þ k P J k9k1 k1k ðzÞ þ p1d t ¼ 1 wtk9k1 qt,1 ðzk Þ k j,1

ð33Þ

ð34Þ

qj,1 ðzk Þ ¼ N ðzk ; z^ k9k1 ,H1k Pjk9k1 ðH1k ÞT þ R1k þ Bj,1 Þ k k9k1

ð35Þ

j,1 mj,1 ðz Þ ¼ mjk9k1 þ K j,1 ðzk z^ k9k1 Þ k k9k k

ð36Þ

j,1 j,1 z^ k9k1 ¼ H1k mjk9k1 þ bk9k1

ð37Þ

P j,1 ¼ ðIK j,1 H1k ÞP jk9k1 k k9k

ð38Þ

¼ P jk9k1 ðH1k ÞT ½H1k Pjk9k1 ðH1k ÞT þ R1k 1 K j,1 k

ð39Þ

ð19Þ

where gk ðxk ,bk Þ is given by (16), and J k1 Y L X

j,1

wj,1 ðzk ÞN ðxk ; mj,1 ðz Þ,Pj,1 ÞN ðbk ; bk9k ,Bj,1 Þ k k9k k k9k k9k

N ðbk ; bk9k1 ,Bj,l Þ k9k1

j¼1l¼1

ns,k9k1 ¼ ps

L X Y

j¼1l¼2

j,l

ð18Þ

ð32Þ

where

wjk1 N ð½xk1 ; bk1 ; ½mjk19k1 ; bk19k1 ,

Þ ½P jk19k1 ; Bj,l k19k1

nld,k ðxk ; zk Þ, l ¼ 2, . . . ,L

zk 2Z lk

l

l

X

nlk9k ¼ ð1pld Þnl1 k9k þ

j,l

wjk1 N ð½xk ; bk ; ½mjs,k9k1 ; bs,k9k1 ,

j¼1l¼1

Þ ½P js,k9k1 ; Bj,l s,k9k1

nb,k9k1 ¼

J b,k J k1 X L Y X

ð20Þ

wjk1 wib,k N ðxk ; mj,i ,Pj,i Þ b,k9k1 b,k9k1

j,1

l j,i,l N ðbk ; bb,k9k1 ,Bj,i,l Þ b,k9k1

j,1

j,1

bk9k ¼ bk9k1 þ Bj,1 ðz z^ k9k1 Þ k9k1 k

j¼1i¼1l¼1

ð21Þ

ð40Þ

Bj,1 ¼ Bj,1 Bj,1 ½H1k P jk9k1 ðH1k ÞT þBj,1 þR1k 1 Bj,1 k9k k9k1 k9k1 k9k1 k9k1

mjs,k9k1 ¼ F k1 mjk19k1

ð22Þ

ð41Þ

P js,k9k1 ¼ F k1 P jk19k1 F Tk1 þQ xk1

ð23Þ

Remark 1. It should be pointed out that the proposed two-stage Kalman filtering process is not optimal in the updated step. This idea is used to maintain the same form of the updated posterior intensity at each time step, i.e., j,l the state estimates mjk9k and the bias estimates bk9k can be expressed in a decoupled form. Another feature of the proposed filter is that the nonlinear target dynamics and measurement models can be addressed by using nonlinear filtering techniques such as the extended Kalman filter (EKF), the unscented Kalman filter (UKF) and the cubature Kalman filter (CKF).

i

¼ F ix,k mjk19k1 þ dx,k mj,i b,k9k1

ð24Þ

¼ F ix,k P jk19k1 ðF ix,k ÞT þ Q ix,k P j,i b,k9k1

ð25Þ

j,l

j,l

bs,k9k1 ¼ bk19k1

ð26Þ

¼ Bj,l þBlk1 Bj,l s,k9k1 k19k1

ð27Þ

j,i,l

j,l

bb,k9k1 ¼ F i,l b,k bk19k1

ð28Þ

j,l i,l i,l T Bj,i,l ¼ F i,l b,k Bk19k1 ðF b,k Þ þ Q b,k b,k9k1

ð29Þ

Remark 2. It is worth noting that the pruning scheme is required after the updated step since the number of Gaussian components increases without bound as time progresses. A simple pruning procedure has been provided by truncating components that have weak weights

90

W. Li et al. / Signal Processing 93 (2013) 86–99

to mitigate this problem. Interested readers are referred to [18] for the details.

The spontaneous birth RFS is Poisson with intensity

gk ðxk ,bk Þ ¼ 0:1 Remark 3. It is expected that the filtering with sensor registration technique can be extended to the Gaussian mixture cardinalized PHD (GM-CPHD) filter to improve the accuracy of multi-target state estimates. Another possible application is to implement the cardinality balanced multi-target multi-Bernoulli (CBMeMBer) recursion as shown in [39].

2 Y 2 X

l

2 Y

N ðxk ; xk1 ,Q x,k Þ

l¼1

4. Simulation results l

l

N ðbk ; bk1 ,Q lb,k Þ In this section, we present two numerical examples, including the non-maneuvering and maneuvering multitarget tracking, to illustrate the effectiveness of the proposed filter and compare its performance with that of the standard GM-PHD filter without introducing sensor biases. Example 1. Consider a two-dimensional scenario with an unknown and time-varying number of targets. The state is denoted by xk ¼ ðpx,k , p_ x,k ,py,k , p_ y,k ÞT , where ðpx,k ,py,k Þ represents the Cartesian coordinates in the horizontal plane and ðp_ x,k , p_ y,k Þ represents its velocities. The target dynamics is described by the following nearly constant velocity model: 2 3 1 T 0 0 60 1 0 07 6 7 þwk1 ð42Þ xk ¼ 6 7x 4 0 0 1 T 5 k1 0

0

ð45Þ

where m1g,k ¼ ð10 000; 0,20 000; 0ÞT , m2g,k ¼ ð0; 0,30 000; 0ÞT , j,1 j,2 P jg,k ¼ diagf106 ,104 ,106 ,104 g, bg,k ¼ ð500, p=360ÞT , bg,k ¼ j,l 2 T 2 ð400, p=120Þ and Bg,k ¼ diagf100 ,ð0:02p=180Þ g (j,l¼1,2). The intensity of the Poisson RFS of spawn births is given by

bk9k1 ðxk ,bk 9xk1 ,bk1 Þ ¼ 0:05

0

j,l

N ðxk ; mjg,k ,P jg,k ÞN ðbk ; bg,k ,Bj,l g,k Þ

j¼1l¼1

1

where T is the sampling period, and wk1 is zero-mean white Gaussian noise with covariance 2 4 3 T T3 0 0 4 2 6 3 7 6 T T2 0 0 7 6 7 26 2 7 Qk ¼ s 6 ð43Þ T4 T3 7 0 0 6 7 4 2 4 5 3 0 0 T2 T 2 with s ¼ 5. Two sensors are used to generate the range and the bearing measurements 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 ðpx,k slx Þ2 þ ðpy,k sly Þ2 7 6 l zlk ¼ 4 5 þbk þ vlk , l ¼ 1; 2 ð44Þ arctan½ðpx,k slx Þ=ðpy,k sly Þ l

where ðslx ,sly Þ is the location of the lth sensor. bk is the translational measurement bias and the measurement l noise vk is assumed to be zero-mean white Gaussian with l Rk. In the simulations, the first sensor is located at (2000, 15 000) m and the second sensor is located at (18 000, 15 000) m. The biases for two sensors are set to ð500, p=360Þ and ð400, p=120Þ, respectively. The covariance matrix of the measurement noise is R1k ¼ R2k ¼ diagf1002 ,ðp=180Þ2 g. Specially, the CKF is used to handle nonlinear measurements [40]. The number of targets is time-varying due to target appearance and disappearance in the scene at any time.

ð46Þ 4

4

where Q x,k ¼ diagf10 ,400,10 ,400g and Q lb,k ¼ diagf4002 ,ð0:02p=180Þ2 gðl ¼ 1; 2Þ. The intensity of the clutter RFS is assumed to be

klk ðzk Þ ¼ llc Uðzk Þ, l ¼ 1; 2

ð47Þ

l

where lk ¼ 10 (l ¼1,2) and UðÞ is the uniform density over the surveillance region. l The bias bk is modeled by a first order Gauss–Markov process with transition density l

l

l

l

f b ðbk 9bk1 Þ ¼ N ðbk ; bk1 ,Blk1 Þ, Blk1

2

l ¼ 1; 2 2

ð48Þ

where ¼ diagf2 ,ð0:01np=180Þ gðl ¼ 1; 2Þ. The true target trajectories are shown in Fig. 1. To be specific, target 1 starts at time k¼1 with initial position at (10 000, 20 000) m and ends at time k¼100; target 2 is spawned from target 1 at time k¼30 and ends at time k¼70; target 3 starts at time k¼5 with initial position at (0, 30 000) m and ends at time k¼100; target 4 is spawned from target 3 at time k¼40 and ends at time k¼80. In the simulations, the survival probability and the detection probability are set to ps ¼0.99 and pd ¼0.98, respectively. The pruning threshold is taken as T Th ¼ 0:001, the merging threshold U Th ¼ 5, the weight threshold wTh ¼ 0:5 and the maximum number of Gaussian terms J max ¼ 100 (see [18] for the meanings of these parameters). The criterion known as optimal subpattern assignment (OSPA) metric is used for performance evaluation since the OSPA metric captures the differences in cardinality and individual elements between two finite sets [41]. To verify the performance of the proposed filter, the simulation results are obtained from 100 Monte Carlo runs. The position estimates of the proposed filter (shortly denoted as ‘GM-PHD-RE’) and the standard GM-PHD filter for one trial are shown in Fig. 2, the simulation results suggest that the GM-PHD-RE filter provides more accurate tracking performance for almost all the time. This is expected since the transitional measurement biases have been estimated and incorporated in the proposed filter. In Fig. 3, the Monte Carlo averages of the OSPA distance with p ¼2 and c¼1000 are shown versus time. It can be seen that the GM-PHD and the GM-PHD-RE filters produce average errors of approximately 580 m and 200 m, respectively. These results also suggest that the GM-PHDRE filter outperforms the standard GM-PHD filter. In addition, the true and the estimated target numbers are

W. Li et al. / Signal Processing 93 (2013) 86–99

91

x 104

4

True tracks Target birth Target death

3.8 3.6

Y coordinate (m)

3.4 3.2 3 2.8 2.6 2.4 2.2 2 0

0.2

0.4

0.6

0.8 1 1.2 X coordinate (m)

1.4

1.6

1.8

2 x 104

Fig. 1. True target trajectories.

4

x 104 True tracks GM−PHD GM−PHD−RE

3.8 3.6

Y coordinate (m)

3.4 3.2 3 2.8 2.6 2.4 2.2 2 0

0.2

0.4

0.6

0.8 1 1.2 X coordinate (m)

1.4

1.6

1.8

2 x 104

Fig. 2. Position estimates with low clutter rate.

shown in Fig. 4, which indicates that the cardinality statistics can be estimated accurately with some penalty of time delay.

In order to investigate the behavior of the GM-PHD-RE filter to clutter rates, we study a tracking scenario l with high clutter rate lk ¼ 30 (l¼ 1,2). The performance

92

W. Li et al. / Signal Processing 93 (2013) 86–99

700 GM−PHD GM−PHD−RE

OSPA (c=1000, p=2)

600

500

400

300

200

100 0

10

20

30

40

50

60

70

80

90

100

Time (s) Fig. 3. Performance comparison with low clutter rate.

4 Truth GM−PHD−RE

Estimated target number

3.5

3

2.5

2

1.5

1 0

10

20

30

40

50

60

70

80

90

100

Time (s) Fig. 4. True and estimated target numbers against time with low clutter rate.

comparison results are shown in Fig. 5, which indicates that the GM-PHD-RE filter achieves higher accuracy than the GM-PHD filter. Another tracking scenario with two

1

2

sensors having different clutter rates (lk ¼ 10 and lk ¼ 30) is carried out as shown in Fig. 6 and similar conclusions can be drawn.

W. Li et al. / Signal Processing 93 (2013) 86–99

93

600 GM−PHD GM−PHD−RE

550

OSPA (c=1000, p=2)

500 450 400 350 300 250 200 150

0

10

20

30

40

50 Time (s)

60

70

80

90

100

Fig. 5. Performance comparison with high clutter rate.

650 GM−PHD GM−PHD−RE

600 550

OSPA (c=1000, p=2)

500 450 400 350 300 250 200 150 0

10

20

30

40

50

60

70

80

90

100

Time (s) Fig. 6. Performance comparison with different clutter rates.

Example 2. In this example, we consider the multiple maneuvering targets tracking by using jump Markov models. For simplicity, the same surveillance region and

measurement equations as in the non-maneuvering target tracking are used. However, two sensors are located at (12 000, 14 000) m and ð2000; 14 000Þ m, respectively.

94

W. Li et al. / Signal Processing 93 (2013) 86–99

To indicate the target maneuvers, the target dynamics is described by the following coordinated turn model: 2 sinðoTÞ oTÞ 3 1 0  1cosð o o 6 7 6 0 cosðoTÞ 0 sinðoTÞ 7 7xk1 þ wk1 ðoÞ xk ¼ 6 ð49Þ sinðoTÞ 6 0 1cosðoTÞ 1 7 4 5 o o 0 sinðoTÞ 0 cosðoTÞ where o denotes the turn rate, and wk1 ðoÞ is zero-mean white Gaussian noise with covariance 2 4 3 T T3 0 0 6 43 2 7 6 T T2 0 0 7 62 7 2 7 Q k ðoÞ ¼ s ðoÞ6 ð50Þ 6 0 0 T4 T3 7 6 4 2 7 4 5 3 0 0 T2 T 2

In the simulations, three motion models corresponding to different turn rates are used. Model 1 is the coordinated turn model with turn rate o ¼ 0J =s and sð0Þ ¼ 5. Model 2 is the coordinated turn model with clockwise turn rate of o ¼ 4J =s and sð4Þ ¼ 20. Model 3 is the coordinated turn model with counterclockwise turn rate of o ¼ 4J =s and sð4Þ ¼ 20. The switching between three models is governed by a first order Markov chain with known transition probability matrix 2 3 0:8 0:1 0:1 6 P ¼ 4 0:1 0:8 0:1 7 ð51Þ 5 0:1 0:1 0:8

The true target trajectories are shown in Fig. 7. To be specific, target 1 starts at time k¼1 with initial position at (10 000, 20 000) m and ends at time k¼100; target 2 is spawned from target 1 at time k¼30 and ends at time k¼ 70; target 3 starts at time k¼5 with initial position at (0, 30 000) m and ends at time k¼100; target 4 is spawned from target 3 at time k¼ 40 and ends at time k¼80. To deal with the multiple model estimation in the framework of the RFS, the best-fitting Gaussian (BFG) approximation approach is applied to derive the GM-PHD-RE filter. The purpose of the BFG approximation is to express the dynamics of the jump Markov linear system (JMLS) (49) by a linear Gaussian model xk þ 1 ¼ Fk xk þwk

ð52Þ

where wk is a zero-mean white Gaussian random vector with covariance matrix Sk , i.e., wk  N ð0, Sk Þ. In other words, we want to replace the JMLS (given by (49) and referred to as ‘‘A’’) with a single BFG distribution (given by (52) and referred to as ‘‘B’’). Fk and Sk are determined such that the distribution of xk has the same mean and covariance under each model, i.e., Efxk 9Ag ¼ Efxk 9Bg

ð53Þ

Covfxk 9Ag ¼ Covfxk 9Bg

ð54Þ

As stated in [23], the matrices Fk and Sk can be determined by pk þ 1,r ¼

M X

pir pk,i

ð55Þ

i¼1

x 104 3.2 True tracks Target birth Target death

3

Y coordinate (m)

2.8

2.6

2.4

2.2

2

1.8

1.6 −4000

−2000

0

2000

4000

6000

8000

X coordinate (m) Fig. 7. True target trajectories.

10000

12000

14000

16000

W. Li et al. / Signal Processing 93 (2013) 86–99

Fk ¼

M X

pk þ 1,r F rk

ð56Þ

95

Sk ¼ Yk þ 1 Fk Yk FTk

ð58Þ

ek þ 1 ¼ Fk ek

ð59Þ

r¼1

pk þ 1,r ½F rk ðYk þ k Tk Þ½F rk T

ee

þQ rk Fk k Tk FTk

ee

ð57Þ

r¼1

3.2

r Fk

is the system transition matrix, pir is the where transition probability of the Markov chain, pk þ 1,r is the

x 104 True tracks GM−PHD GM−PHD−RE

3 2.8 Y coordinate (m)

2.6 2.4 2.2 2 1.8 1.6 −4000 −2000

0

2000

4000

6000

8000 10000 12000 14000 16000

X coordinate (m) Fig. 8. Position estimates with low clutter rate.

700 GM−PHD GM−PHD−RE

600

OSPA (c=1000, p=2)

Yk þ 1 ¼

M X

500

400

300

200

100 0

10

20

30

40

50

60

70

Time (s) Fig. 9. Performance comparison with low clutter rate.

80

90

100

96

W. Li et al. / Signal Processing 93 (2013) 86–99

probability of the event that model r is in effect during the sampling period ½k,k þ 1Þ. ek and Yk þ 1 are auxiliary variables. It has been shown in [23] that the BFG-based GM-PHD filter outperforms the multiple model GM-PHD filter without interacting. To implement the filtering algorithm, the settings for other parameters are identical to those in

Example 1. The position estimates of the proposed filter and the BFG-based GM-PHD filter in [23] for one trial are shown in Fig. 8. In Fig. 9, the Monte Carlo averages of the OSPA metric are shown versus time. The performance l comparisons with high clutter rate lk ¼ 30 (l ¼1,2) and 1 2 different clutter rates (lk ¼ 10 and lk ¼ 30) are presented in Figs. 10 and 11, respectively. The simulation results

700 GM−PHD GM−PHD−RE

OSPA (c=1000, p=2)

600

500

400

300

200

100 0

10

20

30

40

50

60

70

80

90

100

Time (s) Fig. 10. Performance comparison with high clutter rate.

700 GM−PHD GM−PHD−RE

OSPA (c=1000, p=2)

600

500

400

300

200

100 0

10

20

30

40

50 60 Time (s)

70

80

Fig. 11. Performance comparison with different clutter rates.

90

100

W. Li et al. / Signal Processing 93 (2013) 86–99

97

4 Truth GM−PHD−RE

Estimated target number

3.5

3

2.5

2

1.5

1 0

10

20

30

40

50

60

70

80

90

100

Time (s) Fig. 12. True and estimated target numbers against time with low clutter rate.

demonstrate that the GM-PHD-RE filter performs better than the GM-PHD filter without incorporating sensor biases. Similarly, the true and estimated target number are shown in Fig. 12, which is consistent with the conclusion for non-maneuvering target tracking.

main difficulty is the computation of n1d,k ðxk ; zk Þ. It follows from (8) and (9) that 1

n1d,k ðxk ; zk Þ ¼

R

p1d h ðzk 9xk ,bk Þnk9k1

0 0 k1k ðzk Þ þ p1d h1 ðzk 9x0k ,b0k Þnk9k1 ðx0k ,b0k 9Z 1:L 1:k1 Þdxk dbk

ðA:1Þ 5. Conclusion In this paper, the GM-PHD filter is applied for tracking an unknown and time-varying number of targets with multiple sensors, in which the translational measurement biases are considered and estimated. To derive a decoupled form of the state estimates and the bias estimates in the updated posterior intensity, a suboptimal two-stage Kalman filtering process is used. The proposed filter can be extended to track maneuvering targets that follow Markovian switching. Simulation results show that the proposed filter outperforms the standard GM-PHD filter without incorporating measurement biases.

Acknowledgments This work was supported by the National 973 Program (2012CB821200) and the NSFC (61134005, 60921001, 90916024, 91116016). Appendix A. Derivation of the predicted and updated intensities By applying Lemma 1, the predicted intensity nk9k1 can be derived by substituting (7) and (16)–(18) into the PHD prediction (5). For the updated intensity n1k9k , the

where the likelihood function is 1

1

h ðzk 9xk ,bk Þ ¼ N ðzk ; H1k xk þbk ,R1k Þ

ðA:2Þ

Substituting (30) into the numerator of (A.1) yields J k9k1

L X Y

1

p1d wjk9k1 N ðzk ; H1k xk þ bk ,R1k ÞN ðxk ; mjk9k1 ,P jk9k1 Þ

j¼1l¼1 l

j,l

N ðbk ; bk9k1 ,Bj,l Þ k9k1

ðA:3Þ 1

It can be seen that the target state xk and the bias bk are coupled in the likelihood function, it is difficult to 1 derive the decoupled form of xk and bk in (A.3). To overcome this difficulty, the idea is to develop a twostage estimator as follows. First, fix the bias terms to obtain the state estimates. Second, consider the bias term and obtain the bias estimates. Specifically, by applying Lemma 2, we have 1

N ðzk ; H1k xk þ bk ,R1k ÞN ðxk ; mjk9k1 ,P jk9k1 Þ 1

,Pj,1 ÞN ðzk ; H1k mjk9k1 þbk ,H1k Pjk9k1 ðH1k ÞT þ R1k Þ ¼ N ðxk ; mj,1 k9k k9k ðA:4Þ where 1

¼ mjk9k1 þ K j,1 ðzk H1k mjk9k1 bk Þ mj,1 k k9k

ðA:5Þ

98

W. Li et al. / Signal Processing 93 (2013) 86–99

P j,1 ¼ ðIK j,1 H1k ÞPjk9k1 k k9k

ðA:6Þ

¼ P jk9k1 ðH1k ÞT ½H1k P jk9k1 ðH1k ÞT þ R1k 1 K j,1 k

ðA:7Þ

Note that there are still bias terms in the state estimates (A.5), to derive an explicit expression, we adopt 1 j,1 the suboptimal solution, i.e., bk ¼ bk9k1 . Then, this reduces to (36). Similarly, we can derive the bias estimates 1

j,1

N ðzk ; H1k mjk9k1 þ bk ,H1k P jk9k1 ðH1k ÞT þ R1k ÞN ðbk ; bk9k1 ,Bj,1 Þ k9k1 j,1

¼ N ðzk ; H1k mjk9k1 þbk9k1 ,H1k P jk9k1 ðH1k ÞT þ R1k þBj,1 Þ k9k1 1

j,1

Þ N ðbk ; bk9k ,Bj,1 k9k j,1 bk9k

ðA:8Þ Bj,1 k9k

and are given by (40) and (41), where respectively. Following the same line, the denominator of n1d,k ðxk ; zk Þ in (A.1) can be obtained, and the updated intensity nlk9k (l ¼ 2, . . . ,L) can be derived sequentially for each sensor. References [1] J.J. Sudano, A least square algorithm with covariance weighting for computing the translational and rotational errors between two radar sites, in: Proceedings of the IEEE Aerospace and Electronics Conference, Dayton, OH, 1993, pp. 383–387. [2] N. Nabaa, R.H. Bishop, Solution to a multisensor tracking problem with sensor registration errors, IEEE Transactions on Aerospace and Electronic Systems 35 (1) (1999) 354–363. [3] X. Lin, Y. Bar-Shalom, Multisensor target tracking performance with bias compensation, IEEE Transactions on Aerospace and Electronic Systems 42 (3) (2006) 1139–1149. [4] Y.F. Zhou, H. Leung, B. Martin, Sensor alignment with earthcentered earth-fixed coordinate system, IEEE Transactions on Aerospace and Electronic Systems 35 (2) (1999) 410–416. [5] N. Okello, B. Ristic, Maximum likelihood registration for multiple dissimilar sensors, IEEE Transactions on Aerospace and Electronic Systems 39 (3) (2003) 1074–1083. [6] J.G. Herrero, J. Portas, J. Corredera, On-line multi-sensor registration for data fusion on airport surface, IEEE Transactions on Aerospace and Electronic Systems 43 (1) (2007) 356–370. [7] B. Friedland, Treatment of bias in recursive filtering, IEEE Transactions on Automatic Control 14 (4) (1969) 359–367. [8] M.B. Ignagni, An alternate derivation and extension of Friedland’s two-stage Kalman estimator, IEEE Transactions on Automatic Control 26 (3) (1981) 746–750. [9] A.T. Alouani, T.R. Rice, W.D. Blair, A two-stage filter for state estimation in the presence of dynamical stochastic bias, in: Proceedings of the 1992 American Control Conference, Chicago, IL, USA, June 1992, pp. 1784–1788. [10] B.A. Van Doorn, H.A.P. Blom, Systematic error estimation in multisensor fusion systems, in: Proceedings of SPIE Conference on Signal and Data Processing of Small Targets, vol. 1954, Orlando, FL, April 1993, pp. 1–6. [11] J.Y. Keller, M. Darouach, Optimal two-stage Kalman filter in the presence of random bias, Automatica 33 (9) (1997) 1745–1748. [12] D. Haessig, B. Friedland, Separate-bias estimation with reducedorder Kalman filters, IEEE Transactions on Automatic Control 43 (7) (1998) 983–987. [13] C.S. Hsieh, F.C. Chen, Optimal solution of the two-stage Kalman estimator, IEEE Transactions on Automatic Control 44 (1) (1999) 194–199. [14] R. Mahler, Statistical Multisource-Multitarget Information Fusion, Artech House, Norwood, MA, 2007. [15] R. Mahler, Multitarget Bayes filtering via first-order multitarget moments, IEEE Transactions on Aerospace and Electronic Systems 39 (4) (2003) 1152–1178. [16] B.N. Vo, S. Singh, A. Doucet, Sequential Monte Carlo implementation of the PHD filter for multi-target tracking, in: Proceedings of the 6th International Conference on Information Fusion, Cairns, Australia, 2003, pp. 792–799.

[17] B.N. Vo, S. Singh, A. Doucet, Sequential Monte Carlo methods for Bayesian multi-target filtering with random finite sets, IEEE Transactions on Aerospace and Electronic Systems 41 (4) (2005) 1224–1245. [18] B.N. Vo, W.K. Ma, The Gaussian mixture probability hypothesis density filter, IEEE Transactions on Signal Processing 54 (11) (2006) 4091–4104. [19] D.E. Clark, B.N. Vo, Convergence analysis of the Gaussian mixture PHD filter, IEEE Transactions on Signal Processing 55 (4) (2007) 1204–1211. [20] B.N. Vo, A. Pasha, H.D. Tuan, A Gaussian mixture PHD filter for nonlinear jump Markov models, in: Proceedings of the 45th IEEE Conference on Decision and Control, San Diego, CA, USA, December 13–15, 2006, pp. 3162–3167. [21] K. Punithakumar, T. Kirubarajan, A. Sinha, Multiple model probability hypothesis density filter for tracking maneuvering targets, IEEE Transactions on Aerospace and Electronic Systems 44 (1) (2008) 87–98. [22] S.A. Pasha, B.N. Vo, H.D. Tuan, W.K. Ma, A Gaussian mixture PHD filter for jump Markov system models, IEEE Transactions on Aerospace and Electronic Systems 45 (3) (2009) 919–936. [23] W. Li, Y. Jia, Gaussian mixture PHD filter for jump Markov models based on best-fitting Gaussian approximation, Signal Processing 91 (4) (2011) 1036–1042. [24] N. Nandakumaran, K. Punithakumar, T. Kirubarajan, Improved multi-target tracking using probability hypothesis density smoothing, in: Proceedings of the SPIE Conference on Signal and Processing of Small Targets, vol. 6699, August 2007, pp. 1–6. [25] N. Nandakumaran, R. Tharmarasa, T. Lang, T. Kirubarajan, Gaussian mixture probability hypothesis density smoothing with multistatic sonar, in: Proceedings of the SPIE Conference on Signal Processing, Sensor Fusion and Target Recognition, vol. 6968, March 2008, pp. 1–6. [26] N. Nandakumaran, T. Kirubarajan, Maneuvering target tracking using probability hypothesis density smoothing, in: Proceedings of the SPIE Conference on Signal Processing, Sensor Fusion and Target Recognition, vol. 7336, March 2009, pp. 1–6. [27] W. Li, Y. Jia, J. Du, F. Yu, Gaussian mixture PHD smoother for jump Markov models in multiple maneuvering targets tracking, in: Proceedings of the 2011 American Control Conference, San Francisco, CA, USA, June 29–July 1, 2011, pp. 3024–3029. [28] R. Mahler, B.T. Vo, B.N. Vo, The forward–backward probability hypothesis density smoother, in: Proceedings of the 13th International Conference on Information Fusion, Edinburgh, UK, 2010, pp. 1–8. [29] R. Mahler, B.T. Vo, B.N. Vo, Multi-target forward–backward smoothing with the probability hypothesis density, IEEE Transactions on Aerospace and Electronic Systems 48 (1) (2012) 707–728. [30] B.N. Vo, B.T. Vo, R. Mahler, A closed form solution to the probability hypothesis density smoother, in: Proceedings of the 13th International Conference on Information Fusion, Edinburgh, UK, 2010, pp. 1–8. [31] B.N. Vo, B.T. Vo, R. Mahler, Closed form solutions to forward– backward smoothing, IEEE Transactions on Signal Processing 60 (1) (2012) 2–17. [32] B.N. Vo, S. Singh, W.K. Ma, Tracking multiple speakers with random sets, in: Proceedings of the International Conference on Acoustics, Speech and Signal Processing, Montreal, Canada, 2004, pp. 357–360. [33] R. Mahler, A. El-Fallah, Bayesian unified registration and tracking, in: Proceedings of the SPIE Conference on Signal Processing, Sensor Fusion and Target Recognition, vol. 8050, May 2011, pp. 1–11. [34] B. Ristic, D. Clark, Particle filter for joint estimation of multi-object dynamic state and multi-sensor bias, in: Proceedings of the 37th International Conference on Acoustics, Speech and Signal Processing, Kyoto, Japan, March 2012. [35] F. Lian, C. Han, W. Liu, H. Chen, Joint spatial registration an multi-target tracking using an extended probability hypothesis density filter, IET Radar, Sonar and Navigation 5 (4) (2011) 441–448. [36] R. Mahler, Approximate multisensor CPHD and PHD filters, in: Proceedings of the 13th International Conference on Information Fusion, Edinburgh, UK, 2010, pp. 1–8. [37] S. Nagappa, D. Clark, On the ordering of the sensors in the iteratedcorrector probability hypothesis density (PHD) filter, in: Proceedings of the SPIE Conference on Signal Processing, Sensor Fusion and Target Recognition, vol. 8050, May 2011, pp. 1–6.

W. Li et al. / Signal Processing 93 (2013) 86–99

[38] N.T. Pham, W. Huang, S.H. Ong, Multiple sensor multiple object tracking with GMPHD filter, in: Proceedings of the 10th International Conference on Information Fusion, Quebec, Que, 2007, pp. 1–7. [39] B.T. Vo, B.N. Vo, A. Cantoni, The cardinality balanced multi-target multi-Bernoulli filter and its implementations, IEEE Transactions on Signal Processing 57 (2) (2009) 409–423.

99

[40] I. Arasaratnam, S. Haykin, Cubature Kalman filters, IEEE Transactions on Automatic Control 54 (6) (2009) 1254–1269. [41] D. Schuhmacher, B.T. Vo, B.N. Vo, A consistent metric for performance evaluation of multi-object filters, IEEE Transactions on Signal Processing 56 (8) (2008) 3447–3457.