PHD filter for multi-target tracking with glint noise

Signal Processing 94 (2014) 48–56

Contents lists available at SciVerse ScienceDirect

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

PHD filter for multi-target tracking with glint noise Wenling Li a,n, Yingmin Jia a,b, Junping Du c, Jun Zhang 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), SMSS, Beihang University (BUAA), Beijing 100191, China c 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 Electronic and Information Engineering, Beihang University (BUAA), Beijing 100191, China b

a r t i c l e in f o

abstract

Article history: Received 24 January 2013 Received in revised form 30 May 2013 Accepted 9 June 2013 Available online 20 June 2013

This paper studies the problem of multi-target tracking with glint noise in the formulation of random finite set theory. By using the Student's t-distribution to describe the glint noise statistics, the probability hypothesis density (PHD) filter is extended via augmenting the target state and the noise parameters. To derive a closed-form expression for the extended PHD filter, the prior Gamma distribution for the noise parameters is adopted so that the predicted and the updated intensities can be represented by mixtures of Gaussian– Gamma terms. As the target state and the noise parameters are coupled in the likelihood functions, the variational Bayesian approach is applied to derive approximated distributions so that the updated intensity is represented in the same form as the predicted intensity and the resulting algorithm is recursive. Simulation results are provided via a numerical example to illustrate the effectiveness of the proposed filter. & 2013 Elsevier B.V. All rights reserved.

Keywords: Multi-target tracking Probability hypothesis density filter Glint noise Variational Bayesian

1. Introduction Multi-target tracking has received great attention due to its wide applications in military and civil fields. As new targets may appear or disappear randomly in the surveillance region, tracking multi-target involves estimating an unknown and time-varying number of targets from a given set of measurements with uncertain origin. Many tracking algorithms have been proposed based on data association strategies such as the nearest neighbor (NN), the strongest neighbor (SN), the joint probabilistic data association (JPDA) and the multiple hypothesis tracking (MHT) [1]. Due to its combinatorial nature, a high computational load is often required to resolve the data association problem in multi-target tracking algorithms. Two well-known alternative formulations that avoid data 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] (J. Zhang). 0165-1684/$ - see front matter & 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.sigpro.2013.06.012

association are symmetric measurement equations [2] and the random finite sets (RFS) approach [3]. In the RFS formulation, the target states and measurements are represented by two different random finite sets (RFSs), and as a consequence, the multi-target tracking problem can be addressed in a rigorous Bayesian estimation framework based on the finite set statistics (FISST) theory. By constructing the multi-target transition density and the multi-target likelihood function, the optimal multi-target Bayes filter can be derived. However, it 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 [4]. It is worth mentioning that the resulting PHD filter still requires solving multi-dimensional integrals and the integrals might be also intractable in many cases of interest. Two schemes have been proposed to implement the PHD filter explicitly including the sequential Monte Carlo (SMC) [5] and the Gaussian mixture (GM) [6]. A main

W. Li et al. / Signal Processing 94 (2014) 48–56

drawback of the SMC-PHD filter is high computational cost since a large number of particles have to be sampled. To overcome this disadvantage, the GM-PHD filter was developed for linear target dynamic and measurement models with Gaussian distributions, in which the weights, means and covariance matrices are propagated analytically by the Kalman filter (KF). Moreover, the nonlinear KF counterparts can be directly employed to deal with nonlinear target dynamics and measurement models. The convergence properties of two implementations were analyzed in [5,7]. Many extensions have been recently developed to address different tracking scenarios [8–20]. For multi-target tracking problems, the Gaussian distribution has commonly been used for representing the measurement noise statistics due to its mathematical simplicity and effectiveness. However, it is known that not all the real-world data can be modeled well by Gaussian distribution. In particular for radar tracking systems, changes in the target aspect toward the radar may cause irregular electromagnetic wave reflections and this gives rise to outliers or glint noise. It was found that glint noise has a heavy-tailed probability density function and conventional filtering algorithms are known to show unsatisfactory performance in the presence of glint noise [21]. The statistical properties of the glint noise and its mathematical models have been studied extensively in [22]. Specially, the Student's t-distribution has been used to model the glint noise in [23] while the mixture of Gaussian distributions has been used in [24]. In [25], the glint noise was modeled by the mixture of a Gaussian distribution and a Laplacian distribution. As the Student's t-distribution has been shown to be much less sensitive than the Gaussian distribution to outliers [26–28], it has been used as an image prior [29–31]. A robust multisensor fusion algorithm for target tracking applications has been proposed based on the Student's t-distribution [32]. Note that the use of the Student's t-distribution raises significant difficulties of tractability and the variational Bayesian approach has been employed to derive approximated distributions. Nevertheless, the glint noise modeled by the Student's t-distribution has not been addressed for multi-target tracking in the RFS formulation. In this paper, we attempt to apply the PHD filter to address the problem of multi-target tracking with glint noise. By modeling the glint noise as a Student's tdistribution, we show how the variational Bayesian approach can be used in the PHD filter to derive closedform expressions. Based on the prior Gamma distributions for the parameters of the Student's t-distribution, we propose a novel implementation to the PHD filter by representing the predicted and the updated intensities as the mixtures of Gaussian–Gamma distributions. This is inspired by the idea of the GM-PHD filter for Gaussian distributions. To guarantee the same form of the predicted and the updated intensities, the main difficulty encountered is the computation of the predicted likelihood since the target state and the noise parameters are coupled in the likelihood functions. This is overcome by the variational Bayesian approximation method. In addition, a heuristic dynamics has been used which simply spreads the previous posterior distribution by a scalar factor so

49

that the predicted intensity is still a mixture of Gaussian– Gamma terms. Simulation results are provided to illustrate the effectiveness of the proposed filter. The rest of this paper is organized as follows. The problem of multi-target tracking with glint noise is formulated in Section 2. The implementation to the PHD filter is presented by applying the variational Bayesian approach in Section 3. In Section 4, a numerical example is provided to illustrate the effectiveness of the proposed filter. Conclusion is drawn in Section 5. 2. Problem formulation 2.1. Tracking model Consider the following linear target dynamics and measurement models: xk ¼ F k−1 xk−1 þ wk−1 ð1Þ zk ¼ H k xk þ vk

ð2Þ

where xk ∈Rn and zk ∈Rm denote the target state and the measurement vectors, respectively. F k−1 is the state transition matrix and Hk is the measurement matrix. The process noise wk−1 is assumed to be zero-mean white Gaussian with covariance matrix Q k−1 and the measurement noise vk is assumed to be distributed according to heavy-tailed m-dimensional Student's t-distribution. Specifically, the probability density functions of the target state and the measurement can be represented by px ðxk jxk−1 Þ ¼ N ðxk ; F k−1 xk−1 ; Q k−1 Þ

ð3Þ

pz ðzk jxk ; Rk ; r k Þ ¼ Sðzk ; H k xk ; Rk ; r k Þ

ð4Þ

where N ðx; x; PÞ denotes a Gaussian density with mean x and covariance matrix P. Sðz; z; Σ; νÞ denotes the probability density function of a Student's t-distribution with mean z, precision Σ and degree of freedom ν. Unlike the conventional tracking algorithms with known measurement noise statistics, the parameters Rk ¼ diagfr 1k ; …; r m k g and rk of the Student's t-distribution in (4) are assumed to be unknown and they are estimated together with the target state xk. It is worth mentioning that the Student's tdistribution reduces to the Gaussian distribution as ν-∞ and it becomes uninformative when ν-0. Based on the target dynamics and measurement models (1)–(2), the single target tracking can be formulated as a filtering problem for (3)–(4) in the Bayesian estimation framework. However, it should be pointed out that the use of the Student's t-distribution for measurement noise renders the Bayesian marginalization analytically intractable. To overcome this problem, the variational Bayesian approach can be applied by representing the Student's t-distribution as an infinite mixture of Gaussian terms [26], i.e., the m-dimensional Student's t-distribution can be considered as a two-level generation process
m þ ν " #−ðmþνÞ=2 Γ ðz−zÞT Σðz−zÞ 1=2 2
 Sðz; z; Σ; νÞ ¼ jΣj 1þ ν ν ðνπÞm=2 Γ 2 Z ∞
ν ν ð5Þ ¼ N ðz; z; ðsΣÞ−1 ÞG s; ; ds 2 2 0 R∞ where ΓðaÞ ¼ 0 ua−1 e−u du is the gamma function. κ Gðs; κ; θÞ ¼ ðθ =ΓðκÞÞsκ−1 e−θs is the probability density

50

W. Li et al. / Signal Processing 94 (2014) 48–56

function of a Gamma distribution in terms of a shape parameter κ and a inverse scale parameter θ. The scaling parameter s in (5) is a latent variable since it is not apparent in the probability density function of the Student's t-distribution. This decomposition of the Student's t-distribution into separate Gaussian and Gamma components allows convenient application of the variational Bayesian method. 2.2. PHD filter for multi-target tracking

where ps is the target surviving probability and px ð  j  Þ is the single-target transition density (3). pr ð  j  Þ is the transition density of the parameters Rk and rk. Skjk−1 ð  j  Þ and Bk ð  Þ denote the intensity of the spawned target RFS and the intensity of the spontaneously birth target RFS, respectively. After receiving the measurements from the sensor at time step k, the updated intensity can be derived as Dkjk ¼ ð1−pd ÞDkjk−1 þ ∑

zk ∈Z k

In the RFS formulation, the collection of individual targets and measurements are treated as set-valued variables so that the multi-target tracking can be cast in the Bayesian filtering. An RFS is a finite set valued random variable, in which the cardinality of the RFS is characterized by a discrete distribution whereas the joint distribution of the elements is described by an appropriate density. To model the target appearance and disappearance randomly at each time step, the RFS can be used to represent the set of the target states. As the number of measurements received from the sensor might be also random due to the presence of the clutter and the timevarying number of targets, the RFS can be adopted to represent the measurement set at each time step. To be specific, assume that there are nk targets with states xk;1 ; …; xk;nk in the surveillance region and mk measurements zk;1 ; …; zk;mk can be received at time step k, then the multi-target state and measurement can represent by [3] X k ≜fxk;1 ; …; xk;nk g⊂X

ð6Þ

Z k ≜fzk;1 ; …; zk;mk g⊂Z

ð7Þ

where X ⊂Rn and Z⊂Rm denote the state and the observation space, respectively. Then, the multi-target tracking can be formulated a Bayesian filtering as follows: given the set of measurements Z 1:k ¼ fZ 1 ; …; Z k g up to time k, the problem is to find the expectation of the posterior density function pðX k jZ 1:k Þ. Although an optimal Bayesian recursion can be derived in terms of the multi-target posterior density functions, this recursion involves multiple integrals and the multi-target posterior density functions are combinatorial, which makes it computationally intractable. To alleviate this intractability, the propagation of the first order moment or the intensity function of multi-target posterior density functions can be used, i.e., the PHD filter provides a computationally cheaper alternative to the optimal multi-target Bayesian recursion [6]. As the measurements noise parameters Rk and rk should be estimated as well as the target states, we can obtain an extended version of the standard PHD recursion by augmenting the state vector ðxk ; Rk ; r k Þ. Define the posterior intensity Dk−1jk−1 ðxk−1 ; Rk−1 ; r k−1 jZ 1:k−1 Þ, the predicted intensity Dkjk−1 ðxk ; Rk ; r k jZ 1:k−1 Þ, and the posterior intensity Dkjk ðxk ; Rk ; r k jZ 1:k Þ. For simplicity, Dsjt ðxs ; Rs ; r s jZ 1:t Þ is shortly denoted by Dsjt . The predicted intensity can be derived as Z ½ps px ðxk jxk−1 Þpr ðRk ; r k jRk−1 ; r k−1 Þ Dkjk−1 ¼ þSkjk−1 ðxk ; Rk ; r k jxk−1 ; Rk−1 ; r k−1 Þ Dk−1jk−1 dxk−1 dRk−1 dr k−1 þ Bk ðxk ; Rk ; r k Þ

R

pd pz ðzk jxk ; Rk ; r k ÞDkjk−1 pd pz ðzk jx′k ; R′k ; r′k ÞDkjk−1 dx′k dR′k dr′k ð9Þ

where pd is the detection probability, pz ð  j  Þ is the singletarget measurement likelihood (4), κk ð  Þ denotes the intensity of the clutter RFS. The aim of this paper is to develop a closed-form solution to the PHD recursion (8)–(9). As the measurement noise is modeled by the Student's t-distribution, the PHD recursion does not admit tractable solutions in general due to the multi-dimensional integrals. In the following section, the variational Bayesian approach is adopted to give approximate expressions. Before we end this section, the variational Bayesian approach is briefly reviewed for deriving a recursive filter of (1)–(2).

2.3. Variational Bayesian approach For linear dynamical systems (1)–(2) with probability density functions (3)–(4), the aim of the Bayesian filtering is to derive the posterior distribution pðxk ; Rk ; r k jZ k Þ where Z k ¼ fz1 ; …; zk g is the cumulative set of measurements up to time step k. The recursive solution consists of the predictive and the update steps, i.e., Z pðxk ; Rk ; r k jZ k−1 Þ ¼ pðxk ; Rk ; r k jxk−1 ; Rk−1 ; r k−1 Þ pðxk−1 ; Rk−1 ; r k−1 jZ k−1 Þ dxk−1 dRk−1 dr k−1 ð10Þ

pðxk ; Rk ; r k jZ k Þ ¼ R

pðzk jxk ; Rk ; r k Þpðxk ; Rk ; r k jZ k−1 Þ pðzk jxk ; Rk ; r k Þpðxk ; Rk ; r k jZ k−1 Þ dxk dRk dr k ð11Þ

It can be observed that the recursion (10)–(11) reverts to the conventional results when the measurement noise is modeled as Gaussian density with known covariance matrix and a closed-form solution for the first two moments can be obtained in the form of the Kalman filter. To derive a closedform expression for (10)–(11) with unknown Student's tdistribution noise parameters, the conjugate prior distributions for Rk and rk are used. To be specific, suppose that the posterior density can be approximated by pðxk−1 ; Rk−1 ; r k−1 jZ k−1 Þ ¼ N ðxk ; x^ k−1jk−1 ; P k−1jk−1 Þ m

ð8Þ

κk ðzk Þ þ

∏ Gðr lk ; αk−1jk−1 ; βk−1jk−1 ÞGðr k ; γ k−1jk−1 ; ηk−1jk−1 Þ

l¼1

ð12Þ

W. Li et al. / Signal Processing 94 (2014) 48–56

where the diagonal elements of the matrices Rk and rk are assumed to be distributed according to Gamma distributions. Then the predictive density can be derived by substituting (12) into (10) Z pðxk ; Rk ; r k jZ k−1 Þ ¼ px ðxk jxk−1 Þpr ðRk ; r k jRk−1 ; r k−1 Þ N ðxk ; x^ k−1jk−1 ; P k−1jk−1 Þ m

 ∏ Gðr lk ; αk−1jk−1 ; βk−1jk−1 Þ Gðr k ; γ k−1jk−1 ; ηk−1jk−1 Þ dxk−1 dRk−1 dr k−1 m

¼ N ðxk ; x^ kjk−1 ; P kjk−1 Þ ∏ Gðr lk ; αkjk−1 ; βkjk−1 Þ l¼1

ð13Þ

where a heuristic approach is adopted as in [32] to generate the predicted parameters of the Gamma distributions x^ kjk−1 ¼ F k−1 x^ k−1jk−1 P kjk−1 ¼

F k−1 P k−1jk−1 F Tk−1

Then the approximated densities can be derived by Q x ðxk Þ ¼ N ðxk ; x^ kjk ; P kjk Þ

ð25Þ

m

Q R ðRk Þ ¼ ∏ Gðr lk ; αlkjk ; βlkjk Þ

ð26Þ

Q r ðr k Þ ¼ Gðr k ; γ kjk ; ηkjk Þ

ð27Þ

l¼1

where

l¼1

Gðr k ; γ kjk−1 ; ηk−1jk−1 Þ

51

xkjk ¼ xkjk−1 þ K k ðzk −H k xkjk−1 Þ

ð28Þ

P kjk ¼ ðI−K k H k ÞP kjk−1

ð29Þ

b k Þ−1 −1 K k ¼ P kjk−1 H Tk ½H k P kjk−1 H Tk þ ðb skR b k ¼ diag R

( 1 αkjk

; …;

β1kjk

αm kjk

ð30Þ

) ð31Þ

βm kjk

ð14Þ þ Q k−1

ð15Þ

αkjk−1 ¼ ρα αk−1jk−1

ð16Þ

βkjk−1 ¼ ρβ βk−1jk−1

ð17Þ

γ kjk−1 ¼ ργ γ k−1jk−1

ð18Þ

ηkjk−1 ¼ ρη ηk−1jk−1

ð19Þ

with ρα ; ρβ ; ργ and ρη being spread factors taken values in ð0; 1. As shown in [28–32], the latent variable sk in Student's t-distribution in (5) should be incorporated to facilitate the variational Bayesian approach so that the posterior density can be approximated by pðxk ; Rk ; r k jZ k Þ≃Q x ðxk ÞQ R ðRk ÞQ r ðr k Þ

αlkjk ¼ 12 þ αlkjk−1

ð32Þ

s k ½zk −H k xkjk−1 ½zk −H k xkjk−1 T þ H k P kjk H Tk g βlkjk ¼ βlkjk−1 þ12 trfb ð33Þ γ kjk ¼ 12 þ γ kjk−1 ηkjk ¼ ηkjk−1 −

ð34Þ

  1 Γ′ðak Þ 1þ −ln bk −b sk 2 Γðak Þ

ð35Þ

ak ¼

γ kjk 1 þ 2ηkjk 2

bk ¼

b k ½zk −Hk xkjk−1 ½zk −H k xkjk−1 T þ H k P kjk H T g γ kjk trfR k þ 2ηkjk 2

ð36Þ

ð37Þ

ð20Þ

where the approximate densities can be formed by minimizing the Kullback–Leibler (KL) divergence between the approximation and the true posterior  KL½Q x ðxk ÞQ R ðRk ÞQ r ðr k Þ∥pðxk ; Rk ; r k Z k Þ Z ¼ Q x ðxk ÞQ R ðRk ÞQ r ðr k Þ   Q ðx ÞQ ðR ÞQ ðr Þ ð21Þ ln x k R k r k dxk dRk dr k pðxk ; Rk ; r k jZ k Þ Using the calculus of variations for minimizing the KLdivergence with respect to the probability densities while keeping the other fixed yields Z  Q R ðRk ÞQ r ðr k Þln pðxk ; Rk ; r k ; zk jZ k−1 Þ dRk dr k Q x ðxk Þ∝exp ð22Þ Z  Q x ðxk ÞQ r ðr k Þln pðxk ; Rk ; r k ; zk jZ k−1 Þ dxk dr k Q R ðRk Þ∝exp ð23Þ Z  Q x ðxk ÞQ R ðRk Þln pðxk ; Rk ; r k ; zk jZ k−1 Þ dxk dRk Q r ðr k Þ∝exp ð24Þ

b sk ¼

ak bk

ð38Þ

3. Main results To derive a closed-form solution to the PHD recursion (8)–(9), the intensities of the birth and spawning RFSs are assumed to be of the following forms: J b;k

m

Bk ðxk ; Rk ; r k Þ ¼ ∑ ∏ wjb;k N ðxk ; mjb;k ; P jb;k Þ j¼1l¼1

Gðr lk ; αj;l ; βj;l ÞGðr k ; γ jk ; ηjk Þ k k

ð39Þ

Skjk−1 ðxk ; Rk ; r k jxk−1 ; Rk−1 ; r k−1 Þ J p;k

m

i

¼ ∑ ∏ wip;k N ðxk ; F ix;k xk−1 þ dx;k ; Q ix;k Þ i¼1l¼1

pðr lk jr lk−1 Þpðr k jr k−1 Þ wjb;k ,

mjb;k ,

P jb;k ,

ð40Þ αj;l , k

βj;l , k

γ jk

ηjk

where J b;k , and are given parameters that determine the shape of the birth intensity. i J p;k , wip;k , F ix;k , dx;k and Q ix;k are given parameters that determine the shape of the spawning intensity.

52

W. Li et al. / Signal Processing 94 (2014) 48–56

Similar to the Gaussian mixture implementations to the standard PHD filter [6], the PHD filter (8)–(9) can be carried out as follows. Prediction step: Given that the posterior intensity at time step k−1 is J k−1

where J kjk−1

m

Dd;k ðxk ; zk Þ ¼ ∑ ∏ wjk ðzk ÞN ðxk ; mjkjk ðzk Þ; P jkjk Þ j¼1l¼1

; βj;l ÞGðr k ; γ jkjk ; ηjkjk Þ Gðr lk ; αj;l kjk kjk

ð59Þ

m

Dk−1jk−1 ¼ ∑ ∏ wjk−1 N ðxk−1 ; mjk−1jk−1 ; P jk−1jk−1 Þ

wjk ðzk Þ ¼

j¼1l¼1

; βj;l ÞGðr k ; γ jk−1jk−1 ; ηjk−1jk−1 Þ Gðr lk ; αj;l k−1jk−1 k−1jk−1

ð41Þ

then the predicted intensity is Dkjk−1 ¼ Ds;kjk−1 þ Dp;kjk−1 þ Bk ðxk ; Rk ; r k Þ

ð42Þ

where Bk ðxk ; Rk ; r k Þ is given by (39), and

pd wjkjk−1 qjk ðzk Þ J

t t κk ðzÞ þ pd ∑tkjk−1 ¼ 1 wkjk−1 qk ðzk Þ

ð60Þ

 1 1 qjk ðzk Þ ¼ exp − lnjP jkjk−1 j− trfðP jkjk−1 Þ−1 P jkjk g 2 2 1 j j j −1 − trfðP kjk−1 Þ ðmkjk −mkjk−1 Þðmjkjk −mjkjk−1 ÞT g 2 m

J k−1

m

Ds;kjk−1 ¼ ps ∑ ∏

j¼1l¼1

; βj;l ÞGðr k ; γ js;kjk−1 ; ηjs;kjk−1 Þ Gðr lk ; αj;l s;kjk−1 s;kjk−1 J k−1

J p;k

ln βj;l −ln Γðαj;l Þ þ ðαj;l −1Þln αj;l  þ ∑ ½αj;l kjk−1 kjk−1 kjk−1 kjk−1 kjk−1

wjk−1 N ðxk ; mjs;kjk−1 ; P js;kjk−1 Þ

l¼1

m 1 ln βj;l −ln Γðαj;l Þ þ ðαj;l −1Þln αj;l  þ lnjP jkjk j− ∑ ½αj;l kjk kjk kjk kjk kjk 2 l¼1

ð43Þ

þγ jkjk−1 ln ηjkjk−1 −ln Γðγ jkjk−1 Þ þ ðγ jkjk−1 −1Þln γ j;l kjk−1  nþ2 j j j j j;l −γ kjk ln ηkjk þ ln Γðγ kjk Þ−ðγ kjk −1Þln γ kjk þ 2

m

Dp;kjk−1 ¼ ∑ ∑ ∏ wjk−1 wip;k N ðxk ; mj;i ; P j;i Þ p;kjk−1 p;kjk−1 j¼1i¼1l¼1

; βj;i;l ÞGðr k ; γ j;i ; ηj;i Þ Gðr lk ; αj;i;l p;kjk−1 p;kjk−1 p;kjk−1 p;kjk−1 mjs;kjk−1 P js;kjk−1

¼ F k−1 mjk−1jk−1

P j;i p;kjk−1

ð45Þ

¼ F k−1 P jk−1jk−1 F Tk−1

mj;i p;kjk−1

¼ F ix;k mjk−1jk−1

ð44Þ

þ

þ Q k−1

i dx;k

¼ F ix;k P jk−1jk−1 ðF ix;k ÞT

þ

ð46Þ ð47Þ

Q ix;k

ð61Þ

j mjkjk ðzk Þ ¼ mjkjk−1 þ K jk ðzk −z^ kjk−1 Þ

ð62Þ

j z^ kjk−1 ¼ H k mjkjk−1

ð63Þ

P jkjk ¼ ðI−K jk H k ÞP jkjk−1

ð64Þ j

b Þ−1 −1 skR K jk ¼ P jkjk−1 H Tk ½H k P jkjk−1 H Tk þ ðb k j

ð48Þ b j ¼ diag R k

8 <αj;1 kjk :βj;1

; …;

9 = αj;m kjk

ð65Þ

ð66Þ

; βj;m kjk

αj;l ¼ ρα αj;l s;kjk−1 k−1jk−1

ð49Þ

¼ ρβ βj;l βj;l s;kjk−1 k−1jk−1

ð50Þ

αj;l ¼ 12 þ αj;l kjk kjk−1

ð67Þ

γ js;kjk−1 ¼ ργ γ jk−1jk−1

ð51Þ

j j j βj;l ¼ βj;l þ 12 trfb s k ½zk −z^ kjk−1 ½zk −z^ kjk−1 T þ H k P jkjk H Tk g kjk kjk−1

ð68Þ

ηjs;kjk−1 ¼ ρη ηjk−1jk−1

ð52Þ

γ jkjk ¼ 12 þ γ jkjk−1

ð69Þ

¼ ρiα αj;l αj;i;l p;kjk−1 k−1jk−1

ð53Þ

ηjkjk

βj;i;l ¼ ρiβ βj;l p;kjk−1 k−1jk−1

ð54Þ

γ j;i p;kjk−1

ρiγ γ jk−1jk−1

ð55Þ

ηj;i ¼ ρiη ηjk−1jk−1 p;kjk−1

ð56Þ

¼

kjk

1 ¼ ηjkjk−1 − 2

ajk ¼

j

bk ¼

Update step: Given that the predicted intensity at time step k−1 is J kjk−1

m

Dkjk−1 ¼ ∑ ∏

j¼1l¼1

j b sk ¼

wjkjk−1 N ðxk ; mjkjk−1 ; P jkjk−1 Þ

; βj;l ÞGðr k ; γ jkjk−1 ; ηjkjk−1 Þ Gðr lk ; αj;l kjk−1 kjk−1

zk ∈Z k

2ηjkjk γ jkjk 2ηjkjk ajk j

bk

þ

1þ

Γ′ðajk Þ Γðajk Þ

# −ln

j j bk −b sk

1 2

ð70Þ

ð71Þ

j

þ

b ½zk −zj trfR ½zk −zjkjk−1 T þ H k P jkjk H Tk g k kjk−1 2

ð72Þ

ð73Þ

ð57Þ

then the posterior intensity is updated as Dkjk ¼ ð1−pd ÞDkjk−1 þ ∑ Dd;k ðxk ; zk Þ

γ jkjk

"

ð58Þ

Remark 1. Note that the computation of qjk ðzk Þ in (61) is different from that of the GM-PHD filter due to the application of the Student's t-distribution. In fact, qjk ðzk Þ is the measurement likelihood function N ðzk ; zjkjk−1 ; H k P jkjk−1 H TK þ Rk Þ in [6]

W. Li et al. / Signal Processing 94 (2014) 48–56 −1000

while the variational Bayesian approach is used to calculate qjk ðzk Þ in (61) (see Appendix for detailed derivations).

−2000 −3000

Y coordinate (m)

Remark 2. As the number of Gaussian components increases without bound as time progresses, the pruning scheme is required after the updated step and a simple pruning procedure has been provided by truncating components that have weak weights [6]. 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) and the unscented Kalman filter (UKF).

−8000 2000

0

1

where T¼1 is the sampling period, and wk−1 is zero-mean white Gaussian noise with covariance matrix 2 4 3 T T3 0 0 6 43 2 7 6 T T2 0 0 7 62 7 7 ð75Þ Q k−1 ¼ 6 6 0 0 T4 T3 7 6 4 2 7 5 4 3 0 0 T2 T 2 The range and the bearing measurements are generated by a radar 2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 ðpx;k −sx Þ2 þ ðpy;k −sy Þ2 5 þ vk ð76Þ zk ¼ 4 arctan½ðpx;k −sx Þ=ðpy;k −sy Þ where ðsx ; sy Þ is the location of the radar. The measurement noise vk is assumed to be Student's t-distribution with Rk ¼ diagf0:05; 28:66g and rk ¼10. In the simulations, the radar is located at ð4000; −8000Þ. Specially, the EKF is used to handle nonlinear measurements. The number of targets is time-varying due to target appearance and disappearance in the scene at any time. For simplicity, the target spawn is not considered in this example. The intensity of the spontaneous birth RFS is 3

2

Bk ðxk ; Rk ; r k Þ ¼ 0:1 ∑ ∏

j¼1l¼1

N ðxk ; mjγ;k ; P jγ;k ÞGðr lk ; αj;l ; βj;l ÞGðr k ; γ jk ; ηjk Þ k k ð77Þ

where

m1γ;k

T

¼ ð4000; 0; −5000; 0Þ ,

m3γ;k ¼ ð3000; 0; −6000; 0ÞT , αj;l k

βj;l k

γ jk

m2γ;k

¼ ð2000; 0; −4000; 0Þ ,

P jγ;k ¼ diagf106 ; ηjk

T

104 ; 106 ; 104 g,

¼ 10, ¼ 4000, ¼ 100 and ¼ 3 ðj ¼ 1; 2; 3; l ¼ 1; 2Þ. The target trajectories are shown in Fig. 1. Specifically, target 1 starts at time k ¼1 with initial position at ð4000; −5000Þ and ends at time k ¼100; target 2 starts at time

2500

3000

3500

4000

4500

5000

5500

600

X coordinate (m) Fig. 1. True and estimated target trajectories.

800 700 600

OSPA (c=1000, p=2)

0

−5000

−7000

500

55 50

400 45

300

40 35

200

50

60

70

80

90

100

100 0

0

10

20

30

40

50

60

70

80

90

100

Time (s) Fig. 2. Performance comparison with respect to OSPA when r k ¼ 10.

800 700 600

OSPA (c=1000, p=2)

0

−4000

−6000

4. Simulation results Consider a two-dimensional scenario with an unknown and time-varying number of targets. The target state is denoted by xk ¼ ðpx;k ; p_ x;k ; py;k ; p_ y;k ÞT , where ðpx;k ; py;k Þ and ðp_ x;k ; p_ y;k Þ denote the target position and velocity, respectively. The target dynamics is described by the nearly constant velocity model 2 3 1 T 0 0 60 1 0 07 6 7 xk ¼ 6 þ wk−1 ð74Þ 7x 4 0 0 1 T 5 k−1

53

500

50 45

400

40

300

35 30

200

50

60

50

60

70

80

90

100

100 0

0

10

20

30

40

70

80

90

10

Time (s) Fig. 3. Performance comparison with respect to OSPA when r k ¼ 1000.

k ¼5 with initial position at ð2000; −4000Þ and ends at time k¼ 85; target 3 starts at time k¼25 with initial position at ð3000; −6000Þ and ends at time k ¼90; target

54

W. Li et al. / Signal Processing 94 (2014) 48–56

4 starts at time k ¼10 with initial position at ð2000; −4000Þ and ends at time k¼100. In the simulations, the survival probability and the detection probability are set to ps ¼0.99 and pd ¼0.98, respectively. The spreading factor ρα ¼ ρβ ¼ ργ ¼ ρη ¼ ρ

in the predicted intensity is taken to be 0.75 and two fixed point iteration steps is adopted to generate the updated intensity. As in [6], 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. To compare the tracking performance, 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 [33]. For simplicity of notation, the proposed filter is shortly denoted by VB-PHD-EKF in the figures. The position estimates of the VB-PHD-EKF for one trial are shown in Fig. 1, the simulation results suggest that the proposed VB-PHD-EKF filter can provide accurate tracking performance for almost all the time. In Fig. 2, the averages of the OSPA distance over 100 Monte Carlo runs with p ¼2 and c ¼1000 are presented versus time. It can be shown that the GM-PHD-EKF and the VB-PHD-EKF filters produce average errors of approximately 64 and 52, respectively. These results also suggest that the VB-PHD-EKF filter outperforms the GM-PHD-EKF filter. This is expected since the measurement noise is modeled as heavy-tailed Student's t-distribution. However, as shown in Fig. 3, the performance of the proposed VB-PHD-EKF is almost

5 4.5

3 2.5 2 1.5 1 0.5 0

10

20

30

40

50

60

70

80

90

100

Time (s) Fig. 4. True and estimated target numbers against time.

95 90

Average of OSPA w.r.t. Time

Estimated target number

4 3.5

85 80 75 70 65 60 55 50 45 0.6

0.65

0.7

0.75

0.8

0.85

0.9

0.95

The value of ρ

Fig. 5. Averaged OSPA with respect to ρ.

1

W. Li et al. / Signal Processing 94 (2014) 48–56

identical to that of the GM-PHD-EKF when the degree of freedom of Student's t-distribution is taken to be 1000. This is reasonable since the Student's t-distribution reduces to the Gaussian distribution as the degree of freedom r k -∞. In addition, it can be observed that high peaks of the OSPA are avoided for the proposed VB-PHDEKF when new targets appear. This might due to the fact that the updated weight wjkjk is calculated in different ways and the correct target state with low weight has been deleted or merged for the GM-PHD-EKF in the pruning and merging step. As shown in Fig. 4, the proposed filter can derive more accurate target numbers, especially when new target emerges at time steps 5 and 10. To illustrate the choice of the spreading factor ρ in the predicted intensity, the averaged OSPA with different values of ρ is shown in Fig. 5. It is suggested that it is better to choose ρ in the interval ½0:7; 0:9. To evaluate the computational requirement of the proposed algorithm, the averaged CPU time is computed in MATLAB 7.1 on a 2.80-GHz 4 CPU Pentium-based computer operating under Windows XP (Professional). The proposed VB-PHD-EKF consumed approximately 3.17 s per sample run over 100 time steps and the GM-PHD-EKF consumed approximately 2.84 s. This is due to the fact that two fixed point iteration steps are used and a complex computation of qjk ðzk Þ is required to generate the updated intensity. 5. Conclusion In this paper, the PHD filter is applied for tracking an unknown and time-varying number of targets in glint noise environment. Based on a Student's t-distribution model for the glint noise statistics, the PHD filter is extended by augmenting the target state and the noise parameters and a novel implementation to the extended PHD filter is developed by using the variational Bayesian approach. The proposed filter is carried out by representing the predicted and the updated intensities as mixtures of Gaussian–Gamma terms. It is expected that the proposed approach can be used to the cardinalized PHD (CPHD) filter to improve the accuracy of multi-target state estimates. Simulation results show that the proposed filter outperforms the GM-PHD filter for glint noise.

55

m

 ∏ Gðr lk ; αj;l ; βj;l ÞGðr k ; γ jkjk−1 ; ηjkjk−1 Þ dxk dRk dr k kjk−1 kjk−1 l¼1

ðA:1Þ It should be pointed out that the above integral cannot be computed explicitly and the variational Bayesian approach provides a way to approximate it. To be specific, by comparing ((11) and A.1), qjk ðzk Þ can be treated as the predictive likelihood function similar to the denominator in (9). By variational Bayesian approach, the predictive log-likelihood function can be written as [26] ln qjk ðzk Þ ¼ KL½Q x ðxk ÞQ R ðRk ÞQ r ðr k Þ∥pðxk ; Rk ; r k jZ k Þ þ Lj

ðA:2Þ

where   Z pðxk ; Rk ; r k ; zk jZ k−1 Þ Lj ¼ dxk dRk dr k Q x ðxk ÞQ R ðRk ÞQ r ðr k Þln Q x ðxk ÞQ R ðRk ÞQ r ðr k Þ ¼ Efln pðxk ; Rk ; r k ; zk jZ k−1 Þg−Efln Q x ðxk ÞQ R ðRk ÞQ r ðr k Þg ðA:3Þ As the first term on the right hand side of (A.2) can be minimized by the variational iteration, we can get qjk ðzk Þ≃expðLj Þ

ðA:4Þ j

To derive L explicitly, we have j

L ¼ Efln pðxk ; Rk ; r k ; zk jZ k−1 Þg−Efln Q x ðxk ÞQ R ðRk ÞQ r ðr k Þg ¼ Efln pðxk ; Rk ; r k jZ k−1 Þpðzk jxk ; Rk ; r k ;k−1 Þg −Efln Q x ðxk ÞQ R ðRk ÞQ r ðr k Þg ( m

; βj;l Þ ¼ E ln N ðxk ; mjkjk−1 ; P jkjk−1 Þ ∏ Gðr lk ; αj;l kjk−1 kjk−1 l¼1

Gðr k ; γ jkjk−1 ; ηjkjk−1 ÞSðzk ; H k xk ; Rk ; r k Þ (

o )

m

; βj;l ÞGðr k ; γ jkjk ; ηjkjk Þ −E ln N ðxk ; mjkjk ; P jkjk Þ ∏ Gðr lk ; αj;l kjk kjk l¼1

(

m

; βj;l Þ ¼ E ln N ðxk ; mjkjk−1 ; P jkjk−1 Þ þ ∑ ln Gðr lk ; αj;l kjk−1 kjk−1 l¼1

ðA:5Þ þln Gðr k ; γ jkjk−1 ; ηjkjk−1 Þ þln Sðzk ; H k xk ; Rk ; r k Þ−ln N ðxk ; mjkjk ; P jkjk Þ ) m

; βj;l Þ−ln Gðr k ; γ jkjk ; ηjkjk Þ − ∑ ln Gðr lk ; αj;l kjk kjk Acknowledgments

l¼1

n 1 1 ln 2π− lnjP jkjk−1 j− trfðP jkjk−1 Þ−1 P jkjk g 2 2 2 1 j j −1 − trfðP kjk−1 Þ ðmkjk −mjkjk−1 Þðmjkjk −mjkjk−1 ÞT g 2

¼− This work was supported by the National Basic Research Program of China (973 Program, 2012CB821200, 2012CB821201), the NSFC (61203044, 61134005, 60921001, 90916024, 91116016) and the Beijing Natural Science Foundation (4132040). Appendix A Derivation of qjk ðzk Þ: By substituting (57) into (9), we can get qjk ðzk Þ ¼

Z

m

ln βj;l −ln Γðαj;l Þ þ ∑ ½αj;l kjk−1 kjk−1 kjk−1 l¼1

−1Þln αj;l −αj;l þðαj;l kjk−1 kjk−1 kjk−1 þγ jkjk−1 ln ηjkjk−1 −ln Γðγ jkjk−1 Þ þ ðγ jkjk−1 −1Þln γ j;l −γ jkjk−1 kjk−1 þ

m n 1 n ln 2π þ lnjP jkjk j þ − ∑ ½αj;l ln βj;l kjk 2 2 2 l ¼ 1 kjk

−ln Γðαj;l Þ þ ðαj;l −1Þln αj;l  kjk kjk kjk pz ðzk jxk ; Rk ; r k ÞN ðxk ; mjkjk−1 ; P jkjk−1 Þ

þ γ jkjk ln ηjkjk −ln Γðγ jkjk Þ þ ðγ jkjk −1Þln γ j;l þ γ jkjk þαj;l kjk kjk

o

56

W. Li et al. / Signal Processing 94 (2014) 48–56

1 1 lnjP jkjk−1 j− trfðP jkjk−1 Þ−1 P jkjk g 2 2 1 j −1 − trfðP kjk−1 Þ ðmjkjk −mjkjk−1 Þðmjkjk −mjkjk−1 ÞT g 2

¼−

m

ln βj;l −ln Γðαj;l Þ þ ðαj;l −1Þ ln αj;l  þ ∑ ½αj;l kjk−1 kjk−1 kjk−1 kjk−1 kjk−1 l¼1

m 1 þ lnjP jkjk j− ∑ ½αj;l ln βj;l −ln Γðαj;l Þ þ ðαj;l −1Þln αj;l  kjk kjk kjk kjk kjk 2 l¼1

þγ jkjk−1 ln ηjkjk−1 −ln Γðγ jkjk−1 Þ þ ðγ jkjk−1 −1Þln γ j;l kjk−1 −γ jkjk ln ηjkjk þ ln Γðγ jkjk Þ−ðγ jkjk −1Þln γ j;l þ kjk

nþ2 2

ðA:6Þ

Then, taking exponentials on both sides of (A.5) yields qjk ðzk Þ in (61). References [1] Y. Bar-Shalom, X.R. Li, Multitarget-Multisensor Tracking: Principles and Techniques, YBS Publishing, Storrs, CT, 1995. [2] E.W. Kamen, Multiple target tracking based on symmetric measurement equations, IEEE Transactions on Automatic Control 37 (3) (1992) 371–374. [3] R. Mahler, Statistical Multisource-Multitarget Information Fusion, Artech House, Norwood, MA, 2007. [4] R. Mahler, Multitarget Bayes filtering via first-order multitarget moments, IEEE Transactions on Aerospace and Electronic Systems 39 (4) (2003) 1152–1178. [5] 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. [6] B.N. Vo, W.K. Ma, The Gaussian mixture probability hypothesis density filter, IEEE Transactions on Signal Processing 54 (11) (2006) 4091–4104. [7] D.E. Clark, B.N. Vo, Convergence analysis of the Gaussian mixture PHD filter, IEEE Transactions on Signal Processing 55 (4) (2007) 1204–1211. [8] B.T. Vo, B.N. Vo, A. Cantoni, Analytic implementations of the cardinalized probability hypothesis density filter, IEEE Transactions on Signal Processing 55 (7) (2007) 3553–3567. [9] 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. [10] 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. [11] 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. [12] 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.

[13] 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. [14] 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. [15] 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. [16] H. Zhang, Z. Jing, S. Hu, Gaussian mixture CPHD filter with gating technique, Signal Processing 89 (8) (2009) 1521–1530. [17] H. Zhang, Z. Jing, S. Hu, Location of multiple emitters based on the sequential PHD filter, Signal Processing 90 (1) (2010) 34–43. [18] W. Li, Y. Jia, J. Du, F. Yu, Gaussian mixture PHD filter for multi-sensor multi-target tracking with registration errors, Signal Processing 93 (1) (2013) 86–99. [19] K. Granstrom, U. Orguner, A PHD filter for tracking multiple extended targets using random matrices, IEEE Transactions on Signal Processing 60 (11) (2012) 5657–5671. [20] K. Granstrom, C. Lundquist, U. Orguner, Extended target tracking using a Gaussian mixture PHD filter, IEEE Transactions on Aerospace and Electronic Systems 48 (4) (2012) 3268–3286. [21] G.A. Hewer, R.D. Martin, J. Zeh, Robust preprocessing for Kalman filtering of glint noise, IEEE Transactions on Aerospace and Electronic Systems 23 (1) (1987) 120–128. [22] J. Kim, M. Tandale, P.K. Menon, Particle filter for ballistic target tracking with glint noise, Journal of Guidance, Control,and Dynamics 33 (6) (2010) 1918–1921. [23] B.H. Borden, M.L. Mumford, A statistical glint radar cross section target model, IEEE Transactions on Aerospace and Electronic Systems 19 (5) (1983) 781–785. [24] I. Bilik, J. Tabrikian, Maneuvering target tracking in the presence of glint using the nonlinear Gaussian mixture Kalman filter, IEEE Transactions on Aerospace and Electronic Systems 46 (1) (2010) 246–262. [25] W. Wu, Target tracking with glint noise, IEEE Transactions on Aerospace and Electronic Systems 29 (1) (1993) 174–185. [26] C. Bishop, Pattern Recognition and Machine Learning, Springer Verlag, New York, 2006. [27] T. Arne, A. Hanssenb, R.E. Hansenc, F. Godtliebsen, EM-estimation and modeling of heavy-tailed processes with the multivariate normal inverse Gaussian distribution, Signal Processing 85 (8) (2005) 1655–1673. [28] M.E. Tipping, N.D. Lawrence, Variational inference for Student-t models: robust Bayesian interpolation and generalised component analysis, Neurocomputing 69 (1–3) (2005) 123–141. [29] G. Chantas, N. Galatsanos, A. Likas, M. Saunders, Variational Bayesian image restoration based on a product of t-distributions image prior, IEEE Transactions on Image Processing 17 (10) (2008) 1795–1805. [30] D.G. Tzikas, A.C. Likas, N.P. Galatsanos, Variational Bayesian sparse kernel-based blind image deconvolution with Student's-t priors, IEEE Transactions on Image Processing 18 (4) (2009) 753–764. [31] J. Gai, R.L. Stevenson, Studentized dynamical system for robust object tracking, IEEE Transactions on Image Processing 20 (1) (2011) 186–199. [32] H. Zhu, H. Leung, Z. He, A variational Bayesian approach to robust sensor fusion based on Student-t-distribution, Information Sciences 221 (1) (2013) 201–214. [33] 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.