- Email: [email protected]
Improved Gaussian mixture probability hypothesis density smoother
Author’s Accepted Manuscript Improved Gaussian Mixture Probability Hypothesis Density Smoother Xiangyu He, Guixi Liu
www.elsevier.com
PII: DOI: Reference:
S0165-1684(15)00284-4 http://dx.doi.org/10.1016/j.sigpro.2015.08.011 SIGPRO5890
To appear in: Signal Processing Received date: 2 May 2015 Revised date: 26 July 2015 Accepted date: 22 August 2015 Cite this article as: Xiangyu He and Guixi Liu, Improved Gaussian Mixture Probability Hypothesis Density Smoother, Signal Processing, http://dx.doi.org/10.1016/j.sigpro.2015.08.011 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Improved Gaussian Mixture Probability Hypothesis Density Smoother Xiangyu He, Guixi Liu School of Mechano-electronic Engineering, Xidian University, Xi’an, 710071, China Corresponding author. Tel.: +86 13693823108 E-mail address: [email protected] (X. He) Abstract: The Gaussian mixture probability hypothesis density (GM-PHD) smoother proposed recently is a closed-form solution to the forward-backward PHD smoother for the linear Gaussian model, it can yield better state estimates than the GM-PHD filter. However, for the standard GM-PHD smoother, when one or more targets disappear during forward filtering, the smoothed PHD will be adjusted improperly in the backward smoothing, thus leading to a target number misestimation problem. In this paper, an improved GM-PHD smoother is proposed to solve such a problem, in which a modified backward corrector is used to adjust the smoothed PHD. Simulated results show that the improved GM-PHD smoother is superior to the standard GM-PHD smoother in both the aspects of target state estimate and target number estimate so that this improved GM-PHD smoother will have an applicable potential in related fields. Keywords: Gaussian mixture, Probability hypothesis density, Filtering, Smoothing, Target tracking. 1. Introduction Multi-target tracking (MTT) is an important research issue with wide applications in both civilian and military fields, such as passive radar tracking, terrain vehicle tracking, sonar image tracking, multi-target visual tracking, etc. The purpose of MTT is to jointly estimate the number of targets and their states from a sequence of observation sets in the presence of clutter, data association uncertainty, detection uncertainty and noise. The classical methods to the MTT problem requires data association that operates in conjunction with filtering [1-3], such as nearest neighbor (NN) [4], joint probabilistic data association (JPDA) [5], probabilistic data association (PDA) [6], and multiple hypothesis tracking (MHT) [7], etc. The data association problem in the classical MTT algorithms results in a huge computational load due to its combinatorial nature, although the combinatoric setting can be effectively used to reduce the computation complexity of the data association algorithm under some particular conditions [8], MTT still remains challenges in theory and application. In recent years, the random finite set (RFS) theory [9] has been used by more and more researchers to tackle with the MTT problem. The probability hypothesis density (PHD) multi-target filter [10] that propagates the posterior intensity function of the RFS of targets in time avoids the combinatorial problem that arises from data association and has attracted considerable interest. The PHD filter is a suboptimal but computationally tractable alternative to the multi-target Bayes filter in RFS framework, it still requires solving multiple integrals that have no closed-form solutions in general. Sequential Monte Carlo implementation of the PHD (SMC-PHD) filter has paved the way for its application to realistic nonlinear non-Gaussian filtering problems [11, 12]. Besides, a closed-form solution to the PHD recursion has also been derived for linear target dynamic and measurement models, called the Gaussian mixture PHD (GM-PHD), which estimates 1
the PHD distribution as a mixture of Gaussian densities [13, 14]. Afterward, many extensions have been developed to solve different MTT problems [15-23]. The filtering is to recursively estimate the current target states based on all measurements up to the current time, while smoothing can yield better estimates than filtering by using more measurements at a later time. To improve the capability of PHD-based filters, a forward-backward PHD smoother has been proposed recently. Similar to the PHD filter, there are two major implementation methods of the forward-backward PHD smoother known as the SMC-PHD smoother [24, 25], and the Gaussian mixture PHD (GM-PHD) smoother [26, 27]. The backward smoothing recursion is the key step of the forward-backward PHD smoother. In the backward smoothing step, the backward corrector, which involves the innovation from forward filtering using measurements beyond the current time, is used to adjust the filtered PHD distribution at current time step. For the standard forward-backward PHD smoother, however, if one or more targets disappear after some time steps of interest, the state transition density of the disappeared targets will become quite smaller, the smoothed PHD distribution will be adjusted incorrectly at the time step after which targets disappear, thus leading to a target number misestimation problem in the PHD smoother. In this paper, an improved GM-PHD smoother is proposed to deal with the target number misestimation problem in the forward-backward PHD smoother. By using the estimated target number of the forward GM-PHD filter as the judgement of targets disappearing, the proposed algorithm modifies backward corrector formula with the modified survival probability. The numerical simulation indicates that the proposed algorithm is a good solution for the target number misestimation problem caused by the disappeared targets. The rest of this paper is organized as follows. In Section 2 a brief background of RFS-based MTT is provided. In addition, the forward-backward PHD smoother is reviewed and the existing drawback of the PHD smoother is illustrated. The proposed algorithm is elaborated in Section 3. In Section 4, the simulated results to compare the performance of our algorithm with that of the standard GM-PHD smoother are presented and discussed. Finally, some meaningful conclusions are drawn in Section 5. 2. Background 2.1. Multi-target tracking review In MTT problems, assume that at time step k, there are nk target states x1k ,, xknk in a state space and mk measurements z1k ,, zkmk received in an observation space Z . Since there is no sequential order on the respective collections of target states and measurements, they can be naturally represented as finite sets, i.e.
X k {x1k ,, xknk }
(1)
Z k {z1k ,, zkmk } Z
(2)
2
where X k and Z k are the target state and measurement sets with nk targets and mk measurements, respectively. Some measurements in Z k may be due to clutter, the number of clutter points is assumed to be Poisson distributed. Denote the multi-target posterior density function as pk k ( X k Z1k: ) , then, the prediction and update equations for the optimal multi-target Bayes filter are described by pk 1 k ( X k 1 Z1:k ) f k 1 k ( X k 1 X ) pk k ( X Z1:k ) s (dX )
(3)
pk 1 k 1 ( X k 1 Z1:k 1 )
g k 1 ( Z k 1 X k 1 ) pk 1 k ( X k 1 Z1:k 1 )
(4)
g k 1 ( Z k 1 X ) pk 1 k ( X Z1:k 1 ) s (dX )
where the parameter s takes the place of the Lebesgue measure [12], pk 1 k denotes the predicted multi-target density, f k 1 k and g k 1 are the multi-target transition density function and the multi-target likelihood function, respectively.
2.2. PHD smoother PHD smoother involves a forward multi-target filtering recursion using the standard PHD filter and then a backward smoothing recursion [25]. The PHD filter can be derived from the optimal multi-target Bayes filter using finite set statistics. Let vk k denote the filtered PHD at time step k, the PHD filter involves a prediction step and an update step that propagate the intensity function vk k recursively in time, i.e. vk 1 k ( x) k 1 k ( x) pS ,k 1 k f k 1 k ( x ), vk k ()
(5)
vk 1 k 1 ( x) vk 1 k ( x) Lk 1 (Z k 1; x)
(6)
where Lk 1 ( Z k 1 ; x) 1 pD,k 1
pD,k 1 g k 1 ( z x)
zZ k 1
where f k 1 k is the single-target transition density;
(7)
k 1 ( z ) pD,k 1 g k 1 ( z ), vk 1 k ()
,
denotes inner product, i.e., the operation of multiple integrals;
k 1 k denotes the PHD of spontaneously target birth; pS ,k 1 k and pD,k 1 are the survival probability and detection probability, respectively; k 1 denotes the intensity function of the clutter RFS, g k 1 is the single-target measurement likelihood. In the backward smoothing step, the smoothed PHD is propagated backward via the backward smoothing recursion [27], i.e.
3
vk t ( x) vk k ( x) Bk t ( x)
(8)
where Bk t (x) is the backward corrector, it can be recursively computed as follows Bk t ( x) 1 pS ,k 1 k pS ,k 1 k Bk 1 t Lk 1 ( Z k 1;), f k 1 k ( x)
(9)
Note that the backward corrector starts with Bt t ( x) 1 . Then the expected target number of the PHD smoother can be written as Nˆ smoother ,k t 1, vk t ( x)
(10)
As seen from Eq. (10), the expected target number of the PHD smoother depends on the smoothed PHD. 2.3. Drawback of the PHD smoother The performance of PHD smoother degrades when one or more targets disappear during tracking. Suppose that one target survives at time step k but disappears at time step k 1 during forward filtering in PHD smoother, note the state corresponding to the disappeared target as xkl at time step k , then the state xkl 1 of the disappeared target at time step k 1 will not exist, so target state transition density f k 1 k ( xkl ) of the disappeared target becomes quite smaller.
As seen from Eqs. (8) and (9), the backward corrector Bk t (x) at time step k are related to the state transition density of each target from time step k to time step k 1 , if one target disappears between the time step k and time step k 1 , the smoothed PHD vk t (x) at time step k from time step t will be adjusted incorrectly by the backward corrector Bk t (x) , then the estimated target number of the PHD smoother will be misestimated according to Eq. (10).
As discussed above, if one or more targets disappear during the forward filtering in PHD smoother, the smoothed PHD will be adjusted incorrectly in the backward smoothing step. Because the estimated target number is obtained by the integral of the smoothed PHD, the target number misestimation problem in PHD smoother will occur, and the smoothing performance of PHD smoother will degrade. 3. Improved GM-PHD smoother Since the GM-PHD smoother is a closed-form smoothing solution to the forward-backward PHD smoother for the linear Gaussian model, the target number misestimation problem described in Section 2.3 is also present in the GM-PHD smoother if one or more targets disappear during the forward filtering, to solve this problem, based on the GM-PHD smoother derived in [27], an improved GM-PHD smoother is proposed in this section. 3.1. Forward filtering step Under linear Gaussian multi-target assumptions, each target follows a linear Gaussian dynamic model and linear Gaussian measurement model, i.e.
4
f k 1 k ( x ) N x; Fk , Qk
(11)
g k 1 ( z x) N z; H k 1 x, Rk 1
(12)
where N x; m, P denotes a Gaussian density with mean m and covariance P . The parameter Fk is the state transition matrix, Qk is the process noise covariance; H k 1 is the observation matrix and Rk 1 is the observation noise covariance. The PHD of spontaneously target birth at time k 1 is a Gaussian mixture of the form
k 1 ( x) w(i,)k 1 N x; m(i,)k 1 , P(,ik)1 J ,k 1
(13)
i 1
Under linear Gaussian assumptions on the target models, if the initial PHD is a Gaussian mixture, then the forward filtering is performed by using the standard GM-PHD filtering algorithm [14]. Assume that the PHD for targets at time step k is a Gaussian mixture given by Jk k
vk k ( x) wk(i k) N x; mk(i k) , Pk(ik) i 1
(14)
Then, the predicted PHD for targets to time step k 1 is also a Gaussian mixture and is given by vk 1 k ( x) vS ,k 1 k ( x) k 1 ( x)
(15)
where Jk k
vS ,k 1 k ( x) pS wk(i k) N x; mS(i,)k 1 k , PS(,ik)1 k i 1
(16)
Suppose that the predicted PHD for targets at time step k 1 is a Gaussian mixture of the form
J k 1 k
vk 1 k ( x) wk(i)1 k N x; mk(i)1 k , Pk(i1) k i 1
(17)
Then the posterior PHD for targets at time step k 1 is also a Gaussian mixture with the form vk 1 k 1 ( x) (1 pD )vk(i)1 k ( x)
J k 1 k
(i ) (i ) (i ) wk 1 ( z ) N x; mk 1 k 1 ( z ), Pk 1 k 1
zZ k 1 i 1
(18)
where wk(i)1 ( z )
pD wk(i)1 k qk(i)1 ( z )
(19)
k 1 ( z ) pD lJk11 k wk(l)1 k qk(l)1 ( z )
Thus, the filtered PHD for targets at forward time step k 1 can be approximated by the weighted sum of Gaussian components: J k 1 k 1
vk 1 k 1 ( x) wk(i)1 k 1 N x; mk(i)1 k 1 , Pk(i1) k 1 i 1
(20)
By using the filtered PHD intensity at forward time step k 1 , the filtered estimated number of targets can be obtained by summing up the filtered weights, i.e. 5
J k 1 k 1
Nˆ k 1 wk(i)1 k 1
(21)
i 1
Then, the filtered estimated target number Nˆ k 1 is used as the judgement of targets disappearing during forward filtering. To solve the problem in PHD smoother, we use the filtered estimated number of targets during forward filtering to build a step function as follows
1 Nˆ k 1 Nˆ k 0 u ( Nˆ k 1 Nˆ k ) ˆ ˆ 0 N k 1 N k 0
(22)
where Nˆ k and Nˆ k 1 are the filtered estimated target number at forward time step k and time step k 1 , respectively. By using the step function defined in Eq. (22), the survival probability used in the backward smoothing recursion can be modified as follows
pS ,k k 1 pS u( Nˆ k 1 Nˆ k )
(23)
where pS ,k k 1 denotes the modified survival probability used to calculate the backward corrector at time step k . As seen from Eqs. (22) and (23), if one or more targets disappear between time step k and time step k 1 , the value of pS ,k k 1 is equal to 0; Otherwise the value of pS ,k k 1 is equal to the value of pS . 3.2. Backward smoothing step In the GM-PHD backward smoothing step, the backward corrector is used to calculate the smoothed PHD. Thus, evaluating the backward corrector is a key step. According to Eqs. (9) and (11), the backward corrector at time step k from time step k 1 can be modified with pS ,k k 1 as follows (see Appendix A for the detailed derivations) Bk k 1 ( x) 1 pS ,k k 1 pD pS ,k k 1 pD
NH
T k 1Fk , Rk 1 H k 1Qk H k 1
( z; x )
(24)
k 1 ( z ) pD rk 1 ( z )
zZ k 1
where rk 1 ( z ) N H k 1 , Rk 1 ( z;), vk 1 k ()
(25)
J k 1 k
wk(i)1 k N H i 1
(i ) T k 1 , Rk 1 H k 1Pk 1 k H k 1
( z; mk(i)1 k )
Further, if the backward corrector at time step k 1 from time step t has the form J k 1 t
Bk 1 t ( x) k(i)1 N C ( i ) ,D( i ) k(i)1 ; x i 1
k 1
k 1
(26)
Then, the generic backward corrector recursion derived in [27] at time step k from time step t can be modified as follows
6
Bk t ( x) 1 pS ,k k 1 pS ,k k 1 (1 pD ) Bk 1 t ( x) pS ,k k 1 pD
Bk1 t ( x; z )
zZ k 1
(27)
k 1 ( z ) pD rk 1 ( z )
1 t ( x; z ) are given, respectively, by 1 t ( x) and Bk where rk1 ( z ) is given by Eq. (25), Bk J k 1 t
Bk 1 t ( x) k(i)1 N C ( i ) ,D ( i ) k(i )1 ; x k 1
i 1
k 1
(28)
C k(i)1 Ck(i)1 Fk D k(i)1 Dk(i)1 Ck(i)1Qk Ck(i)1T
(29)
J k 1 t (i ) Bk1 t ( x; z ) k(i)1 N C( i ) ,D( i ) k 1 ; x k 1 k 1 z i 1
(30)
(i ) (i ) Ck 1 F C k 1 k H k 1
(i ) D k 1
0 Ck(i)1 Ck(i)1 Qk Rk 1 H k 1 H k 1
D (i ) k 1 0
(31)
T
According to Eqs. (8), (20) and (26), the smoothed PHD at time step k 1 from time step t is a Gaussian mixture, and is given by vk 1 t ( x) vk 1 k 1 ( x) Bk 1 t ( x) J k 1 t J k 1 k 1
i 1
(i ) ( j) (i , j ) (i ) k 1wk 1 k 1 k 1 k 1 ( k 1 ) j 1
~ (i , j ) ~ (i , j ) ( (i ) ), P N x; m k 1 k 1 k 1 k 1 k 1
(32)
where
k(i,1j )k 1 ( k(i )1 ) NC (i ) ,D(i ) C (i ) P( j ) k 1
( i )T k 1 k 1 k 1Ck 1
k 1
(i ) k 1
; mk( j 1) k 1
~ (i , j ) ( (i ) ) m( j ) K (i , j ) (i ) C (i ) m( j ) m k 1 k 1 k 1 k 1 k 1 k 1 k 1 k 1 k 1 k 1 k 1
(33)
~ Pk(i1, jk)1 I K k(i,1j )k 1Ck(i)1 Pk(j1) k 1
K k(i,1j )k 1 Pk(j1) k 1Ck(i)1T Ck(i)1Pk(j1) k 1Ck(i)1T Dk(i)1
(34) (35)
1
(36)
The detailed derivations of the generic backward corrector and the smoothed PHD can be found in [27]. With the modification of backward corrector recursion formula by the modified survival probability pS , the proposed algorithm can solve target number misestimation problem in the standard GM-PHD smoother. The pseudo code of the improved GM-PHD smoother is provided in Algorithm 1. In the next section, we analyze the performance of the proposed algorithm by comparing with that of the standard GM-PHD smoother. Algorithm 1. Pseudo code of the improved GM-PHD smoother.
7
Step 1. Initialize the PHD {w0(i ) , m0(i ) , P0(i ) }iJ01 . Step 2. Forward filtering (Same as GM-PHD filtering algorithm [14]) J
Given {wk(i)1 k 1, mk(i)1 k 1, Pk(i )1 k 1}i k11 k 1 and measurement set Z k ; J
Step 2.1. Compute the predicted PHD {wk(i k) 1 , mk(i k) 1 , Pk(ik)1}i k1k 1 . J
Step 2.2. Compute the filtered PHD {wk(i k) , mk(i k) , Pk(ik) }i k1k . Step 3. Compute the filtered estimated number of targets Jk k
Nˆ k round ( wk(i k) ) ; i 1
Step 4. Modify the survival probability if Nˆ k Nˆ k 1 0 u( Nˆ k Nˆ k 1 ) 1 ;
else u( Nˆ k Nˆ k 1 ) 0 ;
end pS ,k 1 k ps u( Nˆ k Nˆ k 1 ) ;
Step 5. Backward smoothing J
J
Given {wt(it) , mt(it) , Pt (ti ) }it1t , {{wt(i )l t l 1, mt(i )l t l 1, Pt (il) t l 1}i t 1l t l 1 }lk1t , measurement sets {Zt l }lk1t and the modified survival probabilities { pS ,t l 1 t l }lk1t ; Step 5.1. Compute the backward corrector Bt k (x) . J
Step 5.2. Compute the smoothed PHD {wt(ik) , mt(ik) , Pt (ki ) }it1k . 4. Simulation results To demonstrate the efficiency of our proposed improved GM-PHD smoother, we consider a two-dimensional scenario with an unknown and time varying number of targets observed in clutter. Each target has survival probability pS 0.99 and follows a linear Gaussian dynamic model. The target state vector xk [ px,k , p x,k , p y ,k , p y ,k ]T is a vector of planar position and velocity at time k, and the measurement is a noisy version of the position. The simulation environment was as follows: AMD A8-6600K APU with Radeon HD(tm) Graphics 3.9 GHz, 4 GB DDR3 1600 Memory, Windows 7, and MATLAB R2012a. The linear Gaussian dynamic model is used as
8
xk 1
1 0 0 0 1 0 0 x q k 0 0 1 k 0 0 0 1
(37)
where 1 s is the sampling period, the process noise qk is the independent and identically distributed zero-mean white Gaussian noise with covariance matrix 0.0625 0.125 Qk 0 0
0.0025 0.005 0.005 0.01
0.125
0
0
0.25
0
0
0 0
(38)
The linear-Gaussian measurement model is described as 1 0 0 0 z k 1 xk 1 rk 1 0 0 1 0
(39)
where the measurement noise rk 1 is the independent and identically distributed zero-mean white Gaussian noise with covariance matrix 0 0.25 Rk 1 0 0 . 25
(40)
New targets can appear spontaneously according to a Poisson point process with intensity function
k w N x; m(i ) , P(i ) 2
(41)
i 1
where m(1) [0,0,0,0]T , m( 2) [0,2,0,2]T , w 0.2 , P(1) P( 2) diag([5,1,5,1]) , and diag () denotes the diagonal matrix. In our simulation, targets can appear or disappear in the scene at any time. Each target is detected with probability pD 0.98 , the maximum number of Gaussian components is J max 100 . The pruning and merging thresholds are Tp 105 and U 5 , respectively. The clutter is modelled as a Poisson RFS with the uniform density in the surveillance
region, and the clutter density is k ( z) 3.3 105 .
Fig. 1. True tracks and measurements. 9
Fig. 2. True target tracks and estimations.
Fig. 3. True target tracks and estimations in x coordinate versus time.
Fig. 4. True target tracks and estimations in y coordinate versus time. Fig. 1 shows a simulated scenario with true target tracks and measurements for a duration of 50 time steps in the presence of the clutter, where the solid lines and the plus signs denote the true target tracks and the measurements, respectively. Fig. 2 shows the target position estimations superimposed on the true target tracks over 50 time steps. Fig. 3 and Fig. 4 show the individual x and y coordinates of the true target tracks and the estimated targets at each time step, respectively. As seen form Fig. 3 and Fig. 4, the motion times of target 1 is [1, 50] s, target 2, target 3 and target 4 disappear after time steps 25, 35, 45, respectively. It is obvious that the standard GM-PHD smoother incurs some errors in estimating the target
10
positions at time steps 25, 35, 45, while the improved GM-PHD smoother has no missed position estimations at those time steps.
Fig. 5. Target number estimations for standard GM-PHD smoother and improved GM-PHD smoother. Fig. 5 shows the target number estimates. As seen, at time steps 25, 35, 45, since one or more targets disappear after those time steps, the number of targets is underestimated by the standard GM-PHD smoother. A comparison shows that the target number misestimation problem discussed above in standard GM-PHD smoother can be well solved by the improved GM-PHD smoother. To further evaluate the performance of the proposed improved GM-PHD smoother, an appropriate metric, known as the optimal subpattern assignment (OSPA) distance [28] is employed as follows m 1 d OSPA X , Y min d c xi , y i n k i1
p
1
p c n m p
(42)
where the parameters are set to p=1 and c=100 m (meters) in our simulation. The average OSPA errors are obtained from 1000 Monte Carlo trials.
Fig.6. Average OSPA errors for improved GM-PHD smoother and standard GM-PHD smoother. Fig. 6 plots the average OSPA errors for improved GM-PHD smoother, standard GM-PHD smoother and GM-PHD filter over 1000 Monte Carlo runs. In this figure, it can be seen that the standard GM-PHD smoother has three high error peaks at the time steps where the target number is misestimated. And what is more, the maximum value of the OSPA
11
average error 50.25m can be obtained at time step 45. As a comparison, it can be found that the proposed algorithm has the satisfactory tracking performance and is better than the standard one. It is apparent that the improved GM-PHD smoother can yield better state estimates than the GM-PHD filter.
Fig. 7. Average OSPA errors for improved GM-PHD smoother with lags of 1 and 2 time steps. Fig. 7 shows the average OSPA errors for the improved GM-PHD smoother with lags of 1 and 2 time steps. We can see that the performance of the proposed algorithm is improved over the GM-PHD filter. It can also be observed that the higher the lag we use in smoother the better the estimations we get, but improved performance is achieved at the cost of additional computational load (i.e., for one MATLAB implementation, the GM-PHD filter requires 0.257s, while the improved GM-PHD smoother with the lag of 1 requires 0.471s, and the improved GM-PHD smoother with the lag of 2 requires 1.087s). 5. Conclusions To solve the target number misestimation problem in the standard GM-PHD smoother, an improved GM-PHD smoother is proposed in this paper, which can eliminate the target number misestimation problem effectively by the modification of backward corrector recursion with the modified survival probability. Simulation results show that the improved GM-PHD smoother achieves better performance compared with the standard GM-PHD smoother when one or more targets disappear during tracking, and can correct target state and target number estimates with promising results, which indicates that this proposed smoother has good application prospects. Acknowledgments This work was supported by the National Ministries Foundation of China (Grant No. Y42013040181), and Fundamental Research Funds for the Central Universities (Grant No. NSIY191414). Appendix A. By using the time variable k 1 instead of t in Eq. (9), we have Bk k 1 ( x) 1 pS ,k k 1 pS ,k k 1 Bk 1 k 1 Lk 1 ( Z k 1 ; ), f k 1 k ( x) 1 pS ,k k 1 pS ,k k 1 Lk 1 ( Z k 1 ; ), f k 1 k ( x)
Substituting Eqs. (7), (11), (12) and (17) into (A.1), the inner product in (A.1) can be written as 12
(A.1)
Lk 1 ( Z k 1 ; ), f k 1 k ( x) 1 pD
zZ k 1
1 pD
pD N z; H k 1 , Rk 1
k 1 ( z ) pD g k 1 ( z ), vk 1 k () p D N z; H k 1 , Rk 1
zZ k 1
J k 1 k
k 1 ( z ) pD N z; H k 1 , Rk 1 , w i 1
1 pD pD
zZ k 1
(i ) k 1 k
N ;m
(i ) k 1 k
,P
, N ; Fk x, Qk
k 1 ( z ) pD wk(i)1 k N z; H k 1 , Rk 1 , N ; mk(i)1 k , Pk(i1) k
N z; H k 1 Fk x, Rk 1 H k 1Qk H kT1
1 pD pD
J k 1 k
k 1 ( z ) pD wk(i)1 k N z; H k 1mk(i)1 k , Rk 1 H k 1 Pk(i1) k H kT1 i 1
NH
1 pD pD
zZ k 1
(i ) k 1 k
N z; H k 1 , Rk 1 , N ; Fk x, Qk J k 1 k i 1
zZ k 1
, N ; Fk x, Qk
T k 1Fk , Rk 1 H k 1Qk H k 1
z; x
J k 1 k
k 1 ( z ) pD wk(i)1 k N H ,R H P( i ) k 1 k 1 k 1 i 1
(A.2)
k 1 k
H kT1
z; m (i ) k 1 k
Then, substituting (A.2) into (A.1) gives Eqs. (18) and (19), respectively. References [1] Y. Bar-Shalom, T.E. Fortmann, Tracking and Data Association, Academic Press, San Diego, 1988. [2] S. Blackman, Multiple Target Tracking with Radar Applications, Artech House, Norwood, 1986. [3] S. Blackman, R. Popoli, Design and Analysis of Modern Tracking Systems, Artech House, Norwood, 1999. [4] R. A. Singer, J. J. Stein, An optimal tracking filter for processing sensor data of imprecisely determined origin in surveillance system, in: Proceedings of the 10th IEEE Conference on Decision and Control, Miami Beach, Florida, USA, December, 1971, pp. 171-175. [5] T. E. Fortmann, Y. Bar-Shalom, M. Scheffe, Sonar tracking of multiple targets using joint probabilistic data association, IEEE Journal of Oceanic Engineering 8, (3) (1983) 173-184. [6] Y. Bar-Shalom, E. Tse, Tracking in a cluttered environment with probabilistic data association, Automatica 11, (5) (1975) 451–460. [7] S. Blackman, Multiple hypothesis tracking for multiple target tracking, IEEE Transactions on Aerospace and Electronic Systems 19, (1) (2004) 5-18. [8] F. Gianfelici, D. Farina, An effective classification framework for brain-computer interfacing based on a combinatoric setting, IEEE Transactions on Signal Processing 60, (3) (2012) 1446-1459. [9] I. Goodman, R. Mahler, H. Nguyen, Mathematics of Data Fusion, Kluwer Academic Publishers, Dordrecht, 1997. [10] R. Mahler, Multitarget Bayes filtering via first-order multitarget moments, IEEE Transactions on Aerospace and Electronic Systems 39, (4) (2003) 1152-1178. [11] 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 of Information Fusion, Cairns, Australia, July 2003, pp. 792-799.
13
[12] B.N. Vo, S. Singh, A. Doucet, Sequential Monte Carlo methods for multi-target filtering with random finite sets, IEEE Transactions on Aerospace and Electronic Systems 41, (4) (2005) 1224-1245. [13] B.N. Vo, W. K. Ma, A closed form solution for the probability hypothesis density filter, in: Proceedings of 8th International Conference on Information Fusion, Philadelphia, USA, July 2005, pp. 856-863. [14] B.N. Vo, W.K. Ma, The Gaussian mixture probability hypothesis density filter, IEEE Transactions on Signal Processing 54, (11) (2006) 4091-4104. [15] 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. [16] 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. [17] 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. [18] 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. [19] H. Zhang, Z. Jing, S. Hu, Gaussian mixture CPHD filter with gating technique, Signal Processing 89 (8) (2009) 1521-1530. [20] H. Zhang, Z. Jing, S. Hu, Location of multiple emitters based on the sequential PHD filter, Signal Processing 90 (1) (2010) 34-43. [21] 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. [22] 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. [23] 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. [24] N. Nadarajah, T. Kirubarajan, T. Lang, M. McDonald, K. Punithakumar, Multitarget tracking using probability hypothesis density smoothing, IEEE Transactions on Aerospace and Electronic Systems 47, (4) (2011) 2344-2360. [25] R. Mahler, B.N. Vo, B.T. Vo, Forward-backward probability hypothesis density smoothing, IEEE Transactions on Aerospace and Electronic Systems 48, (1) (2012) 707-728. [26] B.N. Vo, B.T. Vo, R. Mahler, A closed form solution to the probability hypothesis density smoother, in: Proceedings of 13th Conference on Information Fusion, Edinburgh, UK, July 2010, pp. 1-8. [27] 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.
14
[28] 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.
Highlights ► An improved Gaussian mixture probability hypothesis density smoother is proposed. ► The survival probability used in the backward smoothing recursion is replaced by the modified survival probability. ► A modified backward corrector based on the original formulation of the Gaussian mixture probability hypothesis density smoother is presented. ► Our algorithm overcomes the existing drawbacks in standard Gaussian mixture probability hypothesis density smoother. ► The simulation results show our algorithm achieves better performance compared with standard Gaussian mixture probability hypothesis density smoother.
15
www.elsevier.com
PII: DOI: Reference:
S0165-1684(15)00284-4 http://dx.doi.org/10.1016/j.sigpro.2015.08.011 SIGPRO5890
To appear in: Signal Processing Received date: 2 May 2015 Revised date: 26 July 2015 Accepted date: 22 August 2015 Cite this article as: Xiangyu He and Guixi Liu, Improved Gaussian Mixture Probability Hypothesis Density Smoother, Signal Processing, http://dx.doi.org/10.1016/j.sigpro.2015.08.011 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Improved Gaussian Mixture Probability Hypothesis Density Smoother Xiangyu He, Guixi Liu School of Mechano-electronic Engineering, Xidian University, Xi’an, 710071, China Corresponding author. Tel.: +86 13693823108 E-mail address: [email protected] (X. He) Abstract: The Gaussian mixture probability hypothesis density (GM-PHD) smoother proposed recently is a closed-form solution to the forward-backward PHD smoother for the linear Gaussian model, it can yield better state estimates than the GM-PHD filter. However, for the standard GM-PHD smoother, when one or more targets disappear during forward filtering, the smoothed PHD will be adjusted improperly in the backward smoothing, thus leading to a target number misestimation problem. In this paper, an improved GM-PHD smoother is proposed to solve such a problem, in which a modified backward corrector is used to adjust the smoothed PHD. Simulated results show that the improved GM-PHD smoother is superior to the standard GM-PHD smoother in both the aspects of target state estimate and target number estimate so that this improved GM-PHD smoother will have an applicable potential in related fields. Keywords: Gaussian mixture, Probability hypothesis density, Filtering, Smoothing, Target tracking. 1. Introduction Multi-target tracking (MTT) is an important research issue with wide applications in both civilian and military fields, such as passive radar tracking, terrain vehicle tracking, sonar image tracking, multi-target visual tracking, etc. The purpose of MTT is to jointly estimate the number of targets and their states from a sequence of observation sets in the presence of clutter, data association uncertainty, detection uncertainty and noise. The classical methods to the MTT problem requires data association that operates in conjunction with filtering [1-3], such as nearest neighbor (NN) [4], joint probabilistic data association (JPDA) [5], probabilistic data association (PDA) [6], and multiple hypothesis tracking (MHT) [7], etc. The data association problem in the classical MTT algorithms results in a huge computational load due to its combinatorial nature, although the combinatoric setting can be effectively used to reduce the computation complexity of the data association algorithm under some particular conditions [8], MTT still remains challenges in theory and application. In recent years, the random finite set (RFS) theory [9] has been used by more and more researchers to tackle with the MTT problem. The probability hypothesis density (PHD) multi-target filter [10] that propagates the posterior intensity function of the RFS of targets in time avoids the combinatorial problem that arises from data association and has attracted considerable interest. The PHD filter is a suboptimal but computationally tractable alternative to the multi-target Bayes filter in RFS framework, it still requires solving multiple integrals that have no closed-form solutions in general. Sequential Monte Carlo implementation of the PHD (SMC-PHD) filter has paved the way for its application to realistic nonlinear non-Gaussian filtering problems [11, 12]. Besides, a closed-form solution to the PHD recursion has also been derived for linear target dynamic and measurement models, called the Gaussian mixture PHD (GM-PHD), which estimates 1
the PHD distribution as a mixture of Gaussian densities [13, 14]. Afterward, many extensions have been developed to solve different MTT problems [15-23]. The filtering is to recursively estimate the current target states based on all measurements up to the current time, while smoothing can yield better estimates than filtering by using more measurements at a later time. To improve the capability of PHD-based filters, a forward-backward PHD smoother has been proposed recently. Similar to the PHD filter, there are two major implementation methods of the forward-backward PHD smoother known as the SMC-PHD smoother [24, 25], and the Gaussian mixture PHD (GM-PHD) smoother [26, 27]. The backward smoothing recursion is the key step of the forward-backward PHD smoother. In the backward smoothing step, the backward corrector, which involves the innovation from forward filtering using measurements beyond the current time, is used to adjust the filtered PHD distribution at current time step. For the standard forward-backward PHD smoother, however, if one or more targets disappear after some time steps of interest, the state transition density of the disappeared targets will become quite smaller, the smoothed PHD distribution will be adjusted incorrectly at the time step after which targets disappear, thus leading to a target number misestimation problem in the PHD smoother. In this paper, an improved GM-PHD smoother is proposed to deal with the target number misestimation problem in the forward-backward PHD smoother. By using the estimated target number of the forward GM-PHD filter as the judgement of targets disappearing, the proposed algorithm modifies backward corrector formula with the modified survival probability. The numerical simulation indicates that the proposed algorithm is a good solution for the target number misestimation problem caused by the disappeared targets. The rest of this paper is organized as follows. In Section 2 a brief background of RFS-based MTT is provided. In addition, the forward-backward PHD smoother is reviewed and the existing drawback of the PHD smoother is illustrated. The proposed algorithm is elaborated in Section 3. In Section 4, the simulated results to compare the performance of our algorithm with that of the standard GM-PHD smoother are presented and discussed. Finally, some meaningful conclusions are drawn in Section 5. 2. Background 2.1. Multi-target tracking review In MTT problems, assume that at time step k, there are nk target states x1k ,, xknk in a state space and mk measurements z1k ,, zkmk received in an observation space Z . Since there is no sequential order on the respective collections of target states and measurements, they can be naturally represented as finite sets, i.e.
X k {x1k ,, xknk }
(1)
Z k {z1k ,, zkmk } Z
(2)
2
where X k and Z k are the target state and measurement sets with nk targets and mk measurements, respectively. Some measurements in Z k may be due to clutter, the number of clutter points is assumed to be Poisson distributed. Denote the multi-target posterior density function as pk k ( X k Z1k: ) , then, the prediction and update equations for the optimal multi-target Bayes filter are described by pk 1 k ( X k 1 Z1:k ) f k 1 k ( X k 1 X ) pk k ( X Z1:k ) s (dX )
(3)
pk 1 k 1 ( X k 1 Z1:k 1 )
g k 1 ( Z k 1 X k 1 ) pk 1 k ( X k 1 Z1:k 1 )
(4)
g k 1 ( Z k 1 X ) pk 1 k ( X Z1:k 1 ) s (dX )
where the parameter s takes the place of the Lebesgue measure [12], pk 1 k denotes the predicted multi-target density, f k 1 k and g k 1 are the multi-target transition density function and the multi-target likelihood function, respectively.
2.2. PHD smoother PHD smoother involves a forward multi-target filtering recursion using the standard PHD filter and then a backward smoothing recursion [25]. The PHD filter can be derived from the optimal multi-target Bayes filter using finite set statistics. Let vk k denote the filtered PHD at time step k, the PHD filter involves a prediction step and an update step that propagate the intensity function vk k recursively in time, i.e. vk 1 k ( x) k 1 k ( x) pS ,k 1 k f k 1 k ( x ), vk k ()
(5)
vk 1 k 1 ( x) vk 1 k ( x) Lk 1 (Z k 1; x)
(6)
where Lk 1 ( Z k 1 ; x) 1 pD,k 1
pD,k 1 g k 1 ( z x)
zZ k 1
where f k 1 k is the single-target transition density;
(7)
k 1 ( z ) pD,k 1 g k 1 ( z ), vk 1 k ()
,
denotes inner product, i.e., the operation of multiple integrals;
k 1 k denotes the PHD of spontaneously target birth; pS ,k 1 k and pD,k 1 are the survival probability and detection probability, respectively; k 1 denotes the intensity function of the clutter RFS, g k 1 is the single-target measurement likelihood. In the backward smoothing step, the smoothed PHD is propagated backward via the backward smoothing recursion [27], i.e.
3
vk t ( x) vk k ( x) Bk t ( x)
(8)
where Bk t (x) is the backward corrector, it can be recursively computed as follows Bk t ( x) 1 pS ,k 1 k pS ,k 1 k Bk 1 t Lk 1 ( Z k 1;), f k 1 k ( x)
(9)
Note that the backward corrector starts with Bt t ( x) 1 . Then the expected target number of the PHD smoother can be written as Nˆ smoother ,k t 1, vk t ( x)
(10)
As seen from Eq. (10), the expected target number of the PHD smoother depends on the smoothed PHD. 2.3. Drawback of the PHD smoother The performance of PHD smoother degrades when one or more targets disappear during tracking. Suppose that one target survives at time step k but disappears at time step k 1 during forward filtering in PHD smoother, note the state corresponding to the disappeared target as xkl at time step k , then the state xkl 1 of the disappeared target at time step k 1 will not exist, so target state transition density f k 1 k ( xkl ) of the disappeared target becomes quite smaller.
As seen from Eqs. (8) and (9), the backward corrector Bk t (x) at time step k are related to the state transition density of each target from time step k to time step k 1 , if one target disappears between the time step k and time step k 1 , the smoothed PHD vk t (x) at time step k from time step t will be adjusted incorrectly by the backward corrector Bk t (x) , then the estimated target number of the PHD smoother will be misestimated according to Eq. (10).
As discussed above, if one or more targets disappear during the forward filtering in PHD smoother, the smoothed PHD will be adjusted incorrectly in the backward smoothing step. Because the estimated target number is obtained by the integral of the smoothed PHD, the target number misestimation problem in PHD smoother will occur, and the smoothing performance of PHD smoother will degrade. 3. Improved GM-PHD smoother Since the GM-PHD smoother is a closed-form smoothing solution to the forward-backward PHD smoother for the linear Gaussian model, the target number misestimation problem described in Section 2.3 is also present in the GM-PHD smoother if one or more targets disappear during the forward filtering, to solve this problem, based on the GM-PHD smoother derived in [27], an improved GM-PHD smoother is proposed in this section. 3.1. Forward filtering step Under linear Gaussian multi-target assumptions, each target follows a linear Gaussian dynamic model and linear Gaussian measurement model, i.e.
4
f k 1 k ( x ) N x; Fk , Qk
(11)
g k 1 ( z x) N z; H k 1 x, Rk 1
(12)
where N x; m, P denotes a Gaussian density with mean m and covariance P . The parameter Fk is the state transition matrix, Qk is the process noise covariance; H k 1 is the observation matrix and Rk 1 is the observation noise covariance. The PHD of spontaneously target birth at time k 1 is a Gaussian mixture of the form
k 1 ( x) w(i,)k 1 N x; m(i,)k 1 , P(,ik)1 J ,k 1
(13)
i 1
Under linear Gaussian assumptions on the target models, if the initial PHD is a Gaussian mixture, then the forward filtering is performed by using the standard GM-PHD filtering algorithm [14]. Assume that the PHD for targets at time step k is a Gaussian mixture given by Jk k
vk k ( x) wk(i k) N x; mk(i k) , Pk(ik) i 1
(14)
Then, the predicted PHD for targets to time step k 1 is also a Gaussian mixture and is given by vk 1 k ( x) vS ,k 1 k ( x) k 1 ( x)
(15)
where Jk k
vS ,k 1 k ( x) pS wk(i k) N x; mS(i,)k 1 k , PS(,ik)1 k i 1
(16)
Suppose that the predicted PHD for targets at time step k 1 is a Gaussian mixture of the form
J k 1 k
vk 1 k ( x) wk(i)1 k N x; mk(i)1 k , Pk(i1) k i 1
(17)
Then the posterior PHD for targets at time step k 1 is also a Gaussian mixture with the form vk 1 k 1 ( x) (1 pD )vk(i)1 k ( x)
J k 1 k
(i ) (i ) (i ) wk 1 ( z ) N x; mk 1 k 1 ( z ), Pk 1 k 1
zZ k 1 i 1
(18)
where wk(i)1 ( z )
pD wk(i)1 k qk(i)1 ( z )
(19)
k 1 ( z ) pD lJk11 k wk(l)1 k qk(l)1 ( z )
Thus, the filtered PHD for targets at forward time step k 1 can be approximated by the weighted sum of Gaussian components: J k 1 k 1
vk 1 k 1 ( x) wk(i)1 k 1 N x; mk(i)1 k 1 , Pk(i1) k 1 i 1
(20)
By using the filtered PHD intensity at forward time step k 1 , the filtered estimated number of targets can be obtained by summing up the filtered weights, i.e. 5
J k 1 k 1
Nˆ k 1 wk(i)1 k 1
(21)
i 1
Then, the filtered estimated target number Nˆ k 1 is used as the judgement of targets disappearing during forward filtering. To solve the problem in PHD smoother, we use the filtered estimated number of targets during forward filtering to build a step function as follows
1 Nˆ k 1 Nˆ k 0 u ( Nˆ k 1 Nˆ k ) ˆ ˆ 0 N k 1 N k 0
(22)
where Nˆ k and Nˆ k 1 are the filtered estimated target number at forward time step k and time step k 1 , respectively. By using the step function defined in Eq. (22), the survival probability used in the backward smoothing recursion can be modified as follows
pS ,k k 1 pS u( Nˆ k 1 Nˆ k )
(23)
where pS ,k k 1 denotes the modified survival probability used to calculate the backward corrector at time step k . As seen from Eqs. (22) and (23), if one or more targets disappear between time step k and time step k 1 , the value of pS ,k k 1 is equal to 0; Otherwise the value of pS ,k k 1 is equal to the value of pS . 3.2. Backward smoothing step In the GM-PHD backward smoothing step, the backward corrector is used to calculate the smoothed PHD. Thus, evaluating the backward corrector is a key step. According to Eqs. (9) and (11), the backward corrector at time step k from time step k 1 can be modified with pS ,k k 1 as follows (see Appendix A for the detailed derivations) Bk k 1 ( x) 1 pS ,k k 1 pD pS ,k k 1 pD
NH
T k 1Fk , Rk 1 H k 1Qk H k 1
( z; x )
(24)
k 1 ( z ) pD rk 1 ( z )
zZ k 1
where rk 1 ( z ) N H k 1 , Rk 1 ( z;), vk 1 k ()
(25)
J k 1 k
wk(i)1 k N H i 1
(i ) T k 1 , Rk 1 H k 1Pk 1 k H k 1
( z; mk(i)1 k )
Further, if the backward corrector at time step k 1 from time step t has the form J k 1 t
Bk 1 t ( x) k(i)1 N C ( i ) ,D( i ) k(i)1 ; x i 1
k 1
k 1
(26)
Then, the generic backward corrector recursion derived in [27] at time step k from time step t can be modified as follows
6
Bk t ( x) 1 pS ,k k 1 pS ,k k 1 (1 pD ) Bk 1 t ( x) pS ,k k 1 pD
Bk1 t ( x; z )
zZ k 1
(27)
k 1 ( z ) pD rk 1 ( z )
1 t ( x; z ) are given, respectively, by 1 t ( x) and Bk where rk1 ( z ) is given by Eq. (25), Bk J k 1 t
Bk 1 t ( x) k(i)1 N C ( i ) ,D ( i ) k(i )1 ; x k 1
i 1
k 1
(28)
C k(i)1 Ck(i)1 Fk D k(i)1 Dk(i)1 Ck(i)1Qk Ck(i)1T
(29)
J k 1 t (i ) Bk1 t ( x; z ) k(i)1 N C( i ) ,D( i ) k 1 ; x k 1 k 1 z i 1
(30)
(i ) (i ) Ck 1 F C k 1 k H k 1
(i ) D k 1
0 Ck(i)1 Ck(i)1 Qk Rk 1 H k 1 H k 1
D (i ) k 1 0
(31)
T
According to Eqs. (8), (20) and (26), the smoothed PHD at time step k 1 from time step t is a Gaussian mixture, and is given by vk 1 t ( x) vk 1 k 1 ( x) Bk 1 t ( x) J k 1 t J k 1 k 1
i 1
(i ) ( j) (i , j ) (i ) k 1wk 1 k 1 k 1 k 1 ( k 1 ) j 1
~ (i , j ) ~ (i , j ) ( (i ) ), P N x; m k 1 k 1 k 1 k 1 k 1
(32)
where
k(i,1j )k 1 ( k(i )1 ) NC (i ) ,D(i ) C (i ) P( j ) k 1
( i )T k 1 k 1 k 1Ck 1
k 1
(i ) k 1
; mk( j 1) k 1
~ (i , j ) ( (i ) ) m( j ) K (i , j ) (i ) C (i ) m( j ) m k 1 k 1 k 1 k 1 k 1 k 1 k 1 k 1 k 1 k 1 k 1
(33)
~ Pk(i1, jk)1 I K k(i,1j )k 1Ck(i)1 Pk(j1) k 1
K k(i,1j )k 1 Pk(j1) k 1Ck(i)1T Ck(i)1Pk(j1) k 1Ck(i)1T Dk(i)1
(34) (35)
1
(36)
The detailed derivations of the generic backward corrector and the smoothed PHD can be found in [27]. With the modification of backward corrector recursion formula by the modified survival probability pS , the proposed algorithm can solve target number misestimation problem in the standard GM-PHD smoother. The pseudo code of the improved GM-PHD smoother is provided in Algorithm 1. In the next section, we analyze the performance of the proposed algorithm by comparing with that of the standard GM-PHD smoother. Algorithm 1. Pseudo code of the improved GM-PHD smoother.
7
Step 1. Initialize the PHD {w0(i ) , m0(i ) , P0(i ) }iJ01 . Step 2. Forward filtering (Same as GM-PHD filtering algorithm [14]) J
Given {wk(i)1 k 1, mk(i)1 k 1, Pk(i )1 k 1}i k11 k 1 and measurement set Z k ; J
Step 2.1. Compute the predicted PHD {wk(i k) 1 , mk(i k) 1 , Pk(ik)1}i k1k 1 . J
Step 2.2. Compute the filtered PHD {wk(i k) , mk(i k) , Pk(ik) }i k1k . Step 3. Compute the filtered estimated number of targets Jk k
Nˆ k round ( wk(i k) ) ; i 1
Step 4. Modify the survival probability if Nˆ k Nˆ k 1 0 u( Nˆ k Nˆ k 1 ) 1 ;
else u( Nˆ k Nˆ k 1 ) 0 ;
end pS ,k 1 k ps u( Nˆ k Nˆ k 1 ) ;
Step 5. Backward smoothing J
J
Given {wt(it) , mt(it) , Pt (ti ) }it1t , {{wt(i )l t l 1, mt(i )l t l 1, Pt (il) t l 1}i t 1l t l 1 }lk1t , measurement sets {Zt l }lk1t and the modified survival probabilities { pS ,t l 1 t l }lk1t ; Step 5.1. Compute the backward corrector Bt k (x) . J
Step 5.2. Compute the smoothed PHD {wt(ik) , mt(ik) , Pt (ki ) }it1k . 4. Simulation results To demonstrate the efficiency of our proposed improved GM-PHD smoother, we consider a two-dimensional scenario with an unknown and time varying number of targets observed in clutter. Each target has survival probability pS 0.99 and follows a linear Gaussian dynamic model. The target state vector xk [ px,k , p x,k , p y ,k , p y ,k ]T is a vector of planar position and velocity at time k, and the measurement is a noisy version of the position. The simulation environment was as follows: AMD A8-6600K APU with Radeon HD(tm) Graphics 3.9 GHz, 4 GB DDR3 1600 Memory, Windows 7, and MATLAB R2012a. The linear Gaussian dynamic model is used as
8
xk 1
1 0 0 0 1 0 0 x q k 0 0 1 k 0 0 0 1
(37)
where 1 s is the sampling period, the process noise qk is the independent and identically distributed zero-mean white Gaussian noise with covariance matrix 0.0625 0.125 Qk 0 0
0.0025 0.005 0.005 0.01
0.125
0
0
0.25
0
0
0 0
(38)
The linear-Gaussian measurement model is described as 1 0 0 0 z k 1 xk 1 rk 1 0 0 1 0
(39)
where the measurement noise rk 1 is the independent and identically distributed zero-mean white Gaussian noise with covariance matrix 0 0.25 Rk 1 0 0 . 25
(40)
New targets can appear spontaneously according to a Poisson point process with intensity function
k w N x; m(i ) , P(i ) 2
(41)
i 1
where m(1) [0,0,0,0]T , m( 2) [0,2,0,2]T , w 0.2 , P(1) P( 2) diag([5,1,5,1]) , and diag () denotes the diagonal matrix. In our simulation, targets can appear or disappear in the scene at any time. Each target is detected with probability pD 0.98 , the maximum number of Gaussian components is J max 100 . The pruning and merging thresholds are Tp 105 and U 5 , respectively. The clutter is modelled as a Poisson RFS with the uniform density in the surveillance
region, and the clutter density is k ( z) 3.3 105 .
Fig. 1. True tracks and measurements. 9
Fig. 2. True target tracks and estimations.
Fig. 3. True target tracks and estimations in x coordinate versus time.
Fig. 4. True target tracks and estimations in y coordinate versus time. Fig. 1 shows a simulated scenario with true target tracks and measurements for a duration of 50 time steps in the presence of the clutter, where the solid lines and the plus signs denote the true target tracks and the measurements, respectively. Fig. 2 shows the target position estimations superimposed on the true target tracks over 50 time steps. Fig. 3 and Fig. 4 show the individual x and y coordinates of the true target tracks and the estimated targets at each time step, respectively. As seen form Fig. 3 and Fig. 4, the motion times of target 1 is [1, 50] s, target 2, target 3 and target 4 disappear after time steps 25, 35, 45, respectively. It is obvious that the standard GM-PHD smoother incurs some errors in estimating the target
10
positions at time steps 25, 35, 45, while the improved GM-PHD smoother has no missed position estimations at those time steps.
Fig. 5. Target number estimations for standard GM-PHD smoother and improved GM-PHD smoother. Fig. 5 shows the target number estimates. As seen, at time steps 25, 35, 45, since one or more targets disappear after those time steps, the number of targets is underestimated by the standard GM-PHD smoother. A comparison shows that the target number misestimation problem discussed above in standard GM-PHD smoother can be well solved by the improved GM-PHD smoother. To further evaluate the performance of the proposed improved GM-PHD smoother, an appropriate metric, known as the optimal subpattern assignment (OSPA) distance [28] is employed as follows m 1 d OSPA X , Y min d c xi , y i n k i1
p
1
p c n m p
(42)
where the parameters are set to p=1 and c=100 m (meters) in our simulation. The average OSPA errors are obtained from 1000 Monte Carlo trials.
Fig.6. Average OSPA errors for improved GM-PHD smoother and standard GM-PHD smoother. Fig. 6 plots the average OSPA errors for improved GM-PHD smoother, standard GM-PHD smoother and GM-PHD filter over 1000 Monte Carlo runs. In this figure, it can be seen that the standard GM-PHD smoother has three high error peaks at the time steps where the target number is misestimated. And what is more, the maximum value of the OSPA
11
average error 50.25m can be obtained at time step 45. As a comparison, it can be found that the proposed algorithm has the satisfactory tracking performance and is better than the standard one. It is apparent that the improved GM-PHD smoother can yield better state estimates than the GM-PHD filter.
Fig. 7. Average OSPA errors for improved GM-PHD smoother with lags of 1 and 2 time steps. Fig. 7 shows the average OSPA errors for the improved GM-PHD smoother with lags of 1 and 2 time steps. We can see that the performance of the proposed algorithm is improved over the GM-PHD filter. It can also be observed that the higher the lag we use in smoother the better the estimations we get, but improved performance is achieved at the cost of additional computational load (i.e., for one MATLAB implementation, the GM-PHD filter requires 0.257s, while the improved GM-PHD smoother with the lag of 1 requires 0.471s, and the improved GM-PHD smoother with the lag of 2 requires 1.087s). 5. Conclusions To solve the target number misestimation problem in the standard GM-PHD smoother, an improved GM-PHD smoother is proposed in this paper, which can eliminate the target number misestimation problem effectively by the modification of backward corrector recursion with the modified survival probability. Simulation results show that the improved GM-PHD smoother achieves better performance compared with the standard GM-PHD smoother when one or more targets disappear during tracking, and can correct target state and target number estimates with promising results, which indicates that this proposed smoother has good application prospects. Acknowledgments This work was supported by the National Ministries Foundation of China (Grant No. Y42013040181), and Fundamental Research Funds for the Central Universities (Grant No. NSIY191414). Appendix A. By using the time variable k 1 instead of t in Eq. (9), we have Bk k 1 ( x) 1 pS ,k k 1 pS ,k k 1 Bk 1 k 1 Lk 1 ( Z k 1 ; ), f k 1 k ( x) 1 pS ,k k 1 pS ,k k 1 Lk 1 ( Z k 1 ; ), f k 1 k ( x)
Substituting Eqs. (7), (11), (12) and (17) into (A.1), the inner product in (A.1) can be written as 12
(A.1)
Lk 1 ( Z k 1 ; ), f k 1 k ( x) 1 pD
zZ k 1
1 pD
pD N z; H k 1 , Rk 1
k 1 ( z ) pD g k 1 ( z ), vk 1 k () p D N z; H k 1 , Rk 1
zZ k 1
J k 1 k
k 1 ( z ) pD N z; H k 1 , Rk 1 , w i 1
1 pD pD
zZ k 1
(i ) k 1 k
N ;m
(i ) k 1 k
,P
, N ; Fk x, Qk
k 1 ( z ) pD wk(i)1 k N z; H k 1 , Rk 1 , N ; mk(i)1 k , Pk(i1) k
N z; H k 1 Fk x, Rk 1 H k 1Qk H kT1
1 pD pD
J k 1 k
k 1 ( z ) pD wk(i)1 k N z; H k 1mk(i)1 k , Rk 1 H k 1 Pk(i1) k H kT1 i 1
NH
1 pD pD
zZ k 1
(i ) k 1 k
N z; H k 1 , Rk 1 , N ; Fk x, Qk J k 1 k i 1
zZ k 1
, N ; Fk x, Qk
T k 1Fk , Rk 1 H k 1Qk H k 1
z; x
J k 1 k
k 1 ( z ) pD wk(i)1 k N H ,R H P( i ) k 1 k 1 k 1 i 1
(A.2)
k 1 k
H kT1
z; m (i ) k 1 k
Then, substituting (A.2) into (A.1) gives Eqs. (18) and (19), respectively. References [1] Y. Bar-Shalom, T.E. Fortmann, Tracking and Data Association, Academic Press, San Diego, 1988. [2] S. Blackman, Multiple Target Tracking with Radar Applications, Artech House, Norwood, 1986. [3] S. Blackman, R. Popoli, Design and Analysis of Modern Tracking Systems, Artech House, Norwood, 1999. [4] R. A. Singer, J. J. Stein, An optimal tracking filter for processing sensor data of imprecisely determined origin in surveillance system, in: Proceedings of the 10th IEEE Conference on Decision and Control, Miami Beach, Florida, USA, December, 1971, pp. 171-175. [5] T. E. Fortmann, Y. Bar-Shalom, M. Scheffe, Sonar tracking of multiple targets using joint probabilistic data association, IEEE Journal of Oceanic Engineering 8, (3) (1983) 173-184. [6] Y. Bar-Shalom, E. Tse, Tracking in a cluttered environment with probabilistic data association, Automatica 11, (5) (1975) 451–460. [7] S. Blackman, Multiple hypothesis tracking for multiple target tracking, IEEE Transactions on Aerospace and Electronic Systems 19, (1) (2004) 5-18. [8] F. Gianfelici, D. Farina, An effective classification framework for brain-computer interfacing based on a combinatoric setting, IEEE Transactions on Signal Processing 60, (3) (2012) 1446-1459. [9] I. Goodman, R. Mahler, H. Nguyen, Mathematics of Data Fusion, Kluwer Academic Publishers, Dordrecht, 1997. [10] R. Mahler, Multitarget Bayes filtering via first-order multitarget moments, IEEE Transactions on Aerospace and Electronic Systems 39, (4) (2003) 1152-1178. [11] 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 of Information Fusion, Cairns, Australia, July 2003, pp. 792-799.
13
[12] B.N. Vo, S. Singh, A. Doucet, Sequential Monte Carlo methods for multi-target filtering with random finite sets, IEEE Transactions on Aerospace and Electronic Systems 41, (4) (2005) 1224-1245. [13] B.N. Vo, W. K. Ma, A closed form solution for the probability hypothesis density filter, in: Proceedings of 8th International Conference on Information Fusion, Philadelphia, USA, July 2005, pp. 856-863. [14] B.N. Vo, W.K. Ma, The Gaussian mixture probability hypothesis density filter, IEEE Transactions on Signal Processing 54, (11) (2006) 4091-4104. [15] 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. [16] 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. [17] 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. [18] 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. [19] H. Zhang, Z. Jing, S. Hu, Gaussian mixture CPHD filter with gating technique, Signal Processing 89 (8) (2009) 1521-1530. [20] H. Zhang, Z. Jing, S. Hu, Location of multiple emitters based on the sequential PHD filter, Signal Processing 90 (1) (2010) 34-43. [21] 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. [22] 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. [23] 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. [24] N. Nadarajah, T. Kirubarajan, T. Lang, M. McDonald, K. Punithakumar, Multitarget tracking using probability hypothesis density smoothing, IEEE Transactions on Aerospace and Electronic Systems 47, (4) (2011) 2344-2360. [25] R. Mahler, B.N. Vo, B.T. Vo, Forward-backward probability hypothesis density smoothing, IEEE Transactions on Aerospace and Electronic Systems 48, (1) (2012) 707-728. [26] B.N. Vo, B.T. Vo, R. Mahler, A closed form solution to the probability hypothesis density smoother, in: Proceedings of 13th Conference on Information Fusion, Edinburgh, UK, July 2010, pp. 1-8. [27] 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.
14
[28] 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.
Highlights ► An improved Gaussian mixture probability hypothesis density smoother is proposed. ► The survival probability used in the backward smoothing recursion is replaced by the modified survival probability. ► A modified backward corrector based on the original formulation of the Gaussian mixture probability hypothesis density smoother is presented. ► Our algorithm overcomes the existing drawbacks in standard Gaussian mixture probability hypothesis density smoother. ► The simulation results show our algorithm achieves better performance compared with standard Gaussian mixture probability hypothesis density smoother.
15