Extended GM-PHD Filter For Multitarget Tracking in ... - IEEE Xplore
Extended GM-PHD Filter For Multitarget Tracking in Nonlinear/Non-Gaussian System Ming Lei, Zhongliang Jing and Peng Dong School of Aeronautics and Astronautics, Shanghai Jiaotong University, Shanghai Email: {mlei, zljing, dongpengkty}@sjtu.edu.cn Abstract—The Gaussian mixture probability hypothesis density (GM-PHD) filter involves the joint estimation of number of targets as well as their individual states in linear/nonlinear Gaussian system, however theoretically not suit for dynamics with nonGaussian noise. In this paper we show that for the state transition and likelihood both with non-Gaussian distribution, the prior and posterior PHD still can be formulated by the weighted Gaussian sum and propagating the Gaussian mixtures separately over time, in this way, state estimation in a nonlinear/non-Gaussian system can be approximately recast as state estimation in a set of parallel nonlinear/Gaussian systems, moreover, the reduced rank scaled unscented/ensemble transform variational (RSEV) filtering [11] is applied to each individual nonlinear/Gaussian system for an improved accuracy of estimation of Gaussian pdf. In addition, an implementation of the proposed algorithm is proposed by combining the closed-form recursions with a strategy for pruning/merging to the number of Gaussian components to increase efficiency.
I.
I NTRODUCTION
Multitarget tracking involves the joint estimation of a timevarying number of targets as well as their individual states from a sequence of sets of noisy measurements/clutters with uncertain data association [1]. Analogous to the single-target Bayesian filtering for being an optimal framework for singletarget state estimation, by using the random finite set (RFS) to represent the collections of targets and measurements, the finite set statistics (FISST) was adopted to construct an extended optimal Bayesian framework for multitarget tracking [2], [3]. The optimal multitarget Bayes filter under the framework of FISST, requires to evaluate combinatorial sums of multiple integrals of high dimensionality with a prohibitively large number of combinations, thus the approach is computationally intractable [3], [4]. The probability hypothesis density (PHD) filter is proposed as a suboptimal but computationally tractable alternative to the RFS-based Bayes multitarget estimation [3], [5]. Aim at the implementation of PHD filtering, a Gaussian mixture PHD recursion, called GM-PHD filter, has been suggested for jointly estimating the time-varying number of targets and corresponding state RFS from the observation RFS in linear/nonlinear Gaussian multitarget system [6], [7]. Based on the Gaussian assumption of state noise and observation noise, the prior PHD is approximated by a sum of weighted Gaussian components, and with the weights, means, and covariances being updated separately by the Kalman filter/Extended Kalman filter/Unscented Kalman filter, the posterior PHD is proved also in form of a weighted Gaussian sum and therefore has an approximate close-form solution to the original problem.
The GM-PHD filter theoretically can not processes the dynamics with non-Gaussian noise, e.g., the nonlinear/nonGaussian system in Eqs.(2), through some other approximations such as the SMC implementations [7], [8] being suggested, however shown that they are still suffer from intensive computation due to the combinatorial nature of the densities, especially when the number of targets is large. In this paper, focusing on the multitarget tracking in nonGaussian system, as described in Eqs.(2), we approximately solve the problem through the weighted Gaussian sum, or Gaussian mixture model [9, ch.8], here use the extended Gaussian mixture PHD (EGM-PHD) for being different from the original GM-PHD. It is shown that for the dynamics and observation both with arbitrary distribution, the prior and posterior PHD still can be formulated by the weighted Gaussian sum and propagating the Gaussian mixtures separately over time. In this way, state estimation in a nonlinear/non-Gaussian system can be approximately recast as state estimation in a set of parallel nonlinear/Gaussian systems. In order to achieve an improved accuracy of nonlinear estimation, the reduced rank scaled unscented/ensemble transform variational (RSEV) filtering in [10], [11] can be applied to each individual nonlinear/Gaussian system to evaluate the mean and covariance of the corresponding Gaussian pdf. Concretely, consider that the pdfs 𝑝 (u𝑘 ) and 𝑝 (v𝑘 ) are of the Gaussian mixture model of the non-Gaussian dynamical and observation noise as shown in Eq.(2), the pdfs of state transition and likelihood at time 𝑘, i.e., 𝑝(x𝑘 ∣x𝑘−1 ) and 𝑝(y𝑘 ∣x𝑘 ), can be approximated by 𝑛𝑢𝑘 and 𝑛𝑣𝑘 Gaussian distributions, see Eqs.(3a)(3b) for details. In addition, like the weighted Gaussian sum presented in [12], [13], when the specified number of Gaussian terms of dynamics and observation noise more than 1 (however a common choice is more than 1 of the terms for the pdf approximation), the computation problems of increasing Gaussian terms as time progresses will be more severe than ever. To reduce the computational cost, [6], [7] suggested a pruning procedure that “discarding those with weights below some preset threshold, or by keeping only a certain number of components with the strongest weights”, and merging that “some of the Gaussian components are so close together that they could be accurately approximated by a single Gaussian.” We also adopt this scheme for computation reduction. For performance illustration, we demonstrate the ability of proposed EGM-PHD to estimate the correct number of targets and their trajectories in high-density clutter. The performance
of the new EGM-PHD filter is benchmarked against the reality specification and its efficacy is shown in terms of its ability to pick up the correct number tracks and the improved accuracy of the trajectory estimates. The rest of this paper is organized as follows. Section II provides a summary of multitarget models with arbitrary noise distributions, and PHD filter in the recursive Bayes equations. In section III, we then proceed to consider the multitarget state estimation problem in nonlinear/non-Gaussian systems. Based on the framework of weighted Gaussian sum and the reduced rank scaled unscented/Ensemble transform variational filtering, we derive the reduced rank scaled unscented/ensemble transform variational-based extended Gaussian mixture PHD (RSEV-EGM-PHD) filter as an extension version to the original GM-PHD filter. For convenience, we also outline the major procedures in the RSEV-EGM-PHD equipped with pruning/merging reduction. An example is then given in section IV to illustrate the details in implementing the RSEV-EGM-PHD. Finally we conclude this paper in section V. II.
P ROBLEM F ORMULATION
Let’s consider the state estimation problem in multitarget scenario. The orderless collections of target states and measurements at time 𝑘, can be represented as finite sets, 𝑛𝑥 ×𝑁𝑘 𝑘 x𝑘 = {x𝑘,𝑖 }𝑁 , 𝑖=1 ∈ ℝ
(1a)
y𝑘 =
(1b)
𝑘 {y𝑘,𝑖 }𝑀 𝑖=1
𝑛𝑦 ×𝑀𝑘
∈ℝ
.
where 𝑁𝑘 and 𝑀𝑘 are unknown and variable with time. Similar to that in single-target filtering, the uncertainty is represented by random vector, the originals uncertainty in a multitarget system is characterized by modelling the multitarget state x𝑘 and multiple target measurement y𝑘 as random finite sets (RFS). An RFS x𝑘 is a finite-set-valued random variable. Then the dynamics can be modeled by x𝑘 = ℳ𝑘,𝑘−1 (x𝑘−1 ) + u𝑘−1 , y𝑘 = ℋ𝑘 (x𝑘 ) + v𝑘 ,
(2a) (2b)
where the transition operator ℳ𝑘,𝑘−1 and the observation operator ℋ𝑘 are both possibly nonlinear. The dynamical and observation noise, in terms of u𝑘−1 and v𝑘 respectively, are non-Gaussian and their pdfs, 𝑝(u ) and 𝑝(v𝑘 ), are known and ∑𝑛𝑢𝑘 𝑘𝑢 approximately 𝑝(u𝑘 ) ≈ 𝛼 𝑖=1 𝑘,𝑖 𝑁 (u𝑘 ; 0, Q𝑘,𝑖 ), 𝑝(v𝑘 ) ≈ ∑𝑛𝑣𝑘 𝑣 𝑖=1 𝛼𝑘,𝑖 𝑁 (v𝑘 ; 0, R𝑘,𝑖 ). Then pdf of transition and likelihood, i.e., 𝑝(x𝑘 ∣x𝑘−1 ) and 𝑝(y𝑘 ∣x𝑘 ), can also be approximated by Gaussian’s with same number of terms, such that 𝑛𝑢 𝑘
𝑝(x𝑘 ∣x𝑘−1 ) ≈
∑
𝑢 𝛼𝑘,𝑖 𝑁 (x𝑘 ; ℳ𝑘,𝑘−1 (x𝑘−1 ) , Q𝑘−1,𝑖 ) ,
𝑖=1 𝑛𝑣 𝑘
𝑝(y𝑘 ∣x𝑘 ) ≈
∑
(3a) 𝑣 𝛼𝑘,𝑖 𝑁 (y𝑘 ; ℋ𝑘 (x𝑘 ) , R𝑘,𝑖 ) ,
(3b)
𝑖=1 𝑢 denotes the weight associated where 𝛼𝑘,𝑖 ∑𝑛𝑢𝑘 𝑢 𝑢 satisfies 𝛼𝑘,𝑖 ∈ [0, 1] and 𝑖=1 𝛼𝑘,𝑖 = 1. 𝑣 definite of 𝛼𝑘,𝑖 .
with Eq.(3a) and Similarly for the
A. Multitarget filtering represented by random finite sets The multitarget state x𝑘 , including target motion, birth, spawning, can be described by RFS. For target motion, given multitarget state RFS x𝑘−1 , each target with x𝑘−1,𝑖 ∈ x𝑘−1 survives at time 𝑘 with surviving probability 𝑝𝑆,𝑘 (x𝑘−1 ), its transition pdf is specified by 𝑝(x𝑘 ∣x𝑘−1 ) in Eq.(3) and therefore the target motion is modeled as the RFS 𝑆(x𝑘 ∣x𝑘−1 ). Similarly, the RFS of target birth at time 𝑘 is modeled by Γ𝑘 , the RFS of targets spawning from a target with x𝑘−1 is modeled by 𝐵(x𝑘 ∣x𝑘−1 ). Summarily the RFS of multitarget state x𝑘 is consist of ⎡ ⎤ ⎡ ⎤ ∪ ∪ ∪ ∪ 𝑋𝑘 = ⎣ 𝑆(x𝑘 ∣x𝑘−1 )⎦ ⎣ 𝐵(x𝑘 ∣x𝑘−1 )⎦ Γ𝑘 . x𝑘−1
x𝑘−1
(4) For the RFS measurement model, suppose that a given target x𝑘 ∈ ℝ𝑛𝑥 ×𝑁𝑘 is either detected with probability 𝑝𝐷,𝑘 (x𝑘 ) or missed with probability 1−𝑝𝐷,𝑘 (x𝑘 ), conditional on detection, the pdf of observed measurement y𝑘 from x𝑘 is modeled by 𝑝(y𝑘 ∣x𝑘 ) in Eq.(3b), then count on either {y𝑘 , ⋅ ⋅ ⋅ } when the target is detected or {0} when the target is not detected, an RFS of the obtained measurements that generated by each state x𝑘 is Θ𝑘 (x𝑘 ). Meanwhile, a set 𝐾𝑘 of false measurements or clutter, may be companied with the target originated measurements. If some of the observations y𝑘,𝑖 ’s are due to clutter, it follows a clutter probability density 𝑐𝑘 , and the number of clutter points are assumed to be Poisson distributed with a mean of 𝜆𝑐 . Thus, given a multitarget state x𝑘 , multitarget measurement y𝑘 can be formed by an union of target generated measurements and clutter, ] [ ∪ ∪ Θ𝑘 (x𝑘 ) . (5) 𝑌 𝑘 = 𝐾𝑘 x𝑘
Then RFSs of multitarget evolution and observation described in Eq.(4) and (5), can be formulated in multitarget transition density 𝑝(x𝑘 ∣x𝑘−1 ) and likelihood 𝑝(y𝑘 ∣x𝑘 ) [5]. B. The probability hypothesis density (PHD) filter PHD filtering propagates a first-order statistics of the posterior multitarget state by appropriate assumptions [1], [5]. Particularly, let 𝑣(x𝑘 ∣Y1:𝑘 ) and 𝑣(x𝑘 ∣Y1:𝑘−1 ) be the respective intensities associated with the multitarget posterior pdf 𝑝(x𝑘 ∣Y1:𝑘 ) and the multitarget predicted pdf 𝑝(x𝑘 ∣Y1:𝑘−1 ), by using the finite set statistics (FISST), PHD recursion can be propagated in time via PHD prediction and PHD update [5]: 𝑣(x𝑘 ∣Y1:𝑘−1 ) ∫ = 𝑝𝑆,𝑘 (x𝑘 )𝑝(x𝑘 ∣x𝑘−1 )𝑣(x𝑘−1 ∣Y1:𝑘−1 )𝑑x𝑘−1 + ∫ 𝑣 𝛾 (x𝑘 ) + 𝑣 𝛽 (x𝑘 ∣x𝑘−1 )𝑣(x𝑘−1 ∣Y1:𝑘−1 )𝑑x𝑘−1 ,
(6a)
𝑣(x𝑘 ∣Y1:𝑘 ) = [1 − 𝑝𝐷,𝑘 (x𝑘 )]𝑣(x𝑘 ∣Y1:𝑘−1 )+ ∑ 𝑝𝐷,𝑘 (x𝑘 )𝑝(y𝑘 ∣x𝑘 )𝑣(x𝑘 ∣Y1:𝑘−1 ) ∫ , 𝜅𝑘 (y𝑘 ) + 𝑝𝐷,𝑘 (x𝑘 )𝑝(y𝑘 ∣x𝑘 )𝑣(x𝑘 ∣Y1:𝑘−1 )𝑑x𝑘 y𝑘 ∈Y𝑘
(6b)
where Y1:𝑘 = {y1 , ⋅ ⋅ ⋅ , y𝑘 } is historical observations up to and including time 𝑘. the ∫posterior intensity 𝑣(x𝑘 ∣Y1:𝑘 ) is the density whose integral 𝑆 𝑣(x𝑘 ∣Y1:𝑘 ))𝑑x𝑘 on any region 𝑆 of RFS gives the expected number of elements of 𝑋 that ∫ ˆ𝑘 = ∣𝑋𝑘 ∩ 𝑆∣𝜇𝑠 (𝑑𝑋𝑘 ). 𝑣 𝛽 (x𝑘 ∣x𝑘−1 ) and are in 𝑆, i.e., 𝑁 𝑣 𝛾 (x𝑘 ) denotes the intensities of 𝐵(x𝑘 ∣x𝑘−1 ) and Γ𝑘 at time 𝑘 in Eq.(4), 𝜅𝑘 (y𝑘 ) is the intensity of the clutter RFS 𝐾𝑘 in Eq.(5) and equals 𝜆𝑐 𝑐𝑘 . 𝑝𝑆,𝑘 (⋅) and 𝑝𝐷,𝑘 (⋅) are the survival and detection probabilities.
∑𝑛𝑥𝑎 𝑘−1 𝑥𝑎 𝑥𝑎 ∈ [0, 1] and 𝑖=1 𝛼𝑘−1,𝑖 = 1. Then the prior where 𝛼𝑘−1,𝑖 intensity 𝑣(x𝑘 ∣Y1:𝑘−1 ) to time 𝑘 is also approximately a Gaussian mixture, and accordingly which consists of three terms 𝑣 𝑆 (x𝑘 ∣Y1:𝑘−1 ), 𝑣 𝛽 (x𝑘 ∣Y1:𝑘−1 ) and 𝑣 𝛾 (x𝑘 ) respectively, to the existing targets, the spawned targets and the spontaneous births, such that 𝑣(x𝑘 ∣Y1:𝑘−1 ) = 𝑣 𝑆 (x𝑘 ∣Y1:𝑘−1 ) + 𝑣 𝛽 (x𝑘 ∣Y1:𝑘−1 ) + 𝑣 𝛾 (x𝑘 ) , (10) where
III.
RSEV-EGM-PHD FILTERING FOR NONLINEAR / NON -G AUSSIAN MULTITARGET TRACKING
𝑢 𝑛𝑥𝑎 𝑘−1 𝑛𝑘
A. Nonlinear/non-Gaussian multitarget model Nonlinear/non-Gaussian multitarget model includes certain assumptions on the birth, death, and detection of targets, are summarized: (1) the survival and detection probabilities are both state independent, i.e., 𝑝𝑆,𝑘 (x𝑘 ) = 𝑝𝑆,𝑘 , 𝑝𝐷,𝑘 (x𝑘 ) = 𝑝𝐷,𝑘 .
𝑣 𝑆 (x𝑘 ∣Y1:𝑘−1 ) = 𝑝𝑆,𝑘
(7a) (7b)
(2) the intensities of the birth and spawn RFSs are Gaussian mixtures of the form
= 𝑝𝑆,𝑘
∑∑
∫ 𝑝(x𝑘 ∣x𝑘−1 )𝑣(x𝑘−1 ∣Y1:𝑘−1 )𝑑x𝑘−1 ,
∫ 𝑢 𝑥𝑎 𝛼𝑘,𝑗 𝛼𝑘−1,𝑖 𝑁 (x𝑘 ; ℳ𝑘,𝑘−1 (x𝑘−1 ) , Q𝑘−1,𝑗 ) ×
𝑖=1 𝑗=1
) ( 𝑁 x𝑘−1 ; x𝑎𝑘−1,𝑖 , P𝑎𝑥𝑥,𝑘−1,𝑖 𝑑x𝑘−1 , 𝑢 𝑛𝑥𝑎 𝑘−1 𝑛𝑘
≈ 𝑝𝑆,𝑘
∑∑ 𝑖=1 𝑗=1 𝑛𝑆 𝑘
= 𝑝𝑆,𝑘
∑
) ( 𝑢 𝑥𝑎 𝛼𝑘,𝑗 𝛼𝑘−1,𝑖 𝑁 x𝑘 ; x𝑆𝑘,(𝑖,𝑗) , P𝑆𝑥𝑥,𝑘,(𝑖,𝑗) ,
) ( 𝑆 𝛼𝑘,𝑠 𝑁 x𝑘 ; x𝑆𝑘,𝑠 , P𝑆𝑥𝑥,𝑘,𝑠 ,
𝑠=1
(11)
𝛾
𝑛𝑘 ∑
𝛾 𝛼𝑘,𝑖 𝑁 (x𝑘 ; x𝛾𝑘,𝑖 , P𝛾𝑥𝑥,𝑘,𝑖 ) ,
x𝑆𝑘,(𝑖,𝑗)
P𝑆𝑥𝑥,𝑘,(𝑖,𝑗)
where and the mean and covariance of the background evaluated by SUT [14], that is, to propagate 𝑖=1 forward the 𝑖-th indexed analysis statistics (i.e., x𝑎𝑘−1,𝑖 and 𝛽 𝑛𝑘 𝑎 𝑥𝑎 ∑ , here 𝑖 ∈ [1, 𝑛𝑘−1 ] ) at time 𝑘 − 1, through the 𝑗P 𝛽 𝛼𝑘,𝑖 𝑁 (x𝑘 ; ℳ𝑘,𝑘−1 (x𝑘−1 ) + 𝑑𝛽𝑘−1,𝑖 , Q𝛽𝑘−1,𝑖 ) , th𝑥𝑥,𝑘−1,𝑖 𝑣 𝛽 (x𝑘 ∣x𝑘−1 ) = indexed nonlinear/Gaussian system x𝑘 = ℳ𝑘,𝑘−1 (x𝑘−1 ) + 𝑖=1 𝑁(u𝑘−1,𝑗 ; 0, Q𝑘−1,𝑗 ) , 𝑗 ∈ [1, 𝑛𝑢𝑘 ], then to estimate the first(8b) two moment of transferred statistics. 𝛾 , x𝛾𝑘,𝑖 , P𝛾𝑥𝑥,𝑘,𝑖 , 𝑖 = 1, ⋅ ⋅ ⋅ , 𝑛𝛾𝑘 are given model where 𝑛𝛾𝑘 , 𝛼𝑘,𝑖 𝑆 𝑢 𝑥𝑎 Moreover, 𝑛𝑆𝑘 = 𝑛𝑢𝑘 𝑛𝑥𝑎 𝑘−1 , 𝛼𝑘,𝑠 = 𝛼𝑘,𝑗 𝛼𝑘−1,𝑖 with the parameters that determine the shape of the birth intensity; integer index 𝑠 being a representation of the two-dimensional 𝛽 𝑥𝑎 similarly 𝑛𝛽𝑘 , 𝛼𝑘,𝑖 , ℳ𝑘,𝑘−1 , 𝑑𝛽𝑘−1,𝑖 , Q𝛽𝑘−1,𝑖 , 𝑖 = 1, ⋅ ⋅ ⋅ , 𝑛𝛽𝑘 deindex (𝑖, 𝑗), e.g., 𝑠 = 𝑖 + 𝑛𝑥𝑎 𝑘−1 (𝑗 − 1), 1 ≤ 𝑖 ≤ 𝑛𝑘−1 and termine the shape of the spawning intersity of a target with 1 ≤ 𝑗 ≤ 𝑛𝑢𝑘 . Similarly the previous state x𝑘−1 . ∫ 𝛽 𝛾 𝛾 𝑣 (x ∣Y ) = 𝑣 𝛽 (x𝑘 ∣x𝑘−1 )𝑣(x𝑘−1 ∣Y1:𝑘−1 )𝑑x𝑘−1 , 𝑘 1:𝑘−1 In the second assumption, x𝑘,𝑖 , 𝑖 = 1, ⋅ ⋅ ⋅ , 𝑛𝑘 are the peaks of the spontaneous birth intensity in (8a). These points 𝛽 𝑛𝑥𝑎 ∫ 𝑘−1 𝑛𝑘 ∑ ∑ 𝛽 have the highest local concentrations of expected number 𝛽 𝛽 𝑥𝑎 𝛾 = 𝛼 𝛼 𝑘,𝑗 𝑘−1,𝑖 𝑁 (x𝑘 ; ℳ𝑘,𝑘−1 (x𝑘−1 ) + 𝑑𝑘−1,𝑗 , Q𝑘−1,𝑗 ) of spontaneous births and represent. The weight 𝛼𝑘,𝑖 given 𝑖=1 𝑗=1 the expected number of new targets originating from x𝛾𝑘,𝑖 . × 𝑁 (x𝑘−1 ; x𝑎𝑘−1,𝑖 , P𝑎𝑥𝑥,𝑘−1,𝑖 )𝑑x𝑘−1 , Analogously for the spawning intensity in Eq.(8b). 𝑣 𝛾 (x𝑘 ) =
(8a)
𝛽 𝑛𝑥𝑎 𝑘−1 𝑛𝑘
B. RSEV-based extended GM-PHD (RSEV-EGM-PHD) Filter The multitarget model described with Gaussian mixture model for nonlinear/non-Gaussian system given in III-A, the procedure of RSEV-based extended GM-PHD filtering of Eqs.(6a)(6b) yields a approximated solution called RSEV-EGMPHD recursion, we split the recursion into 3 steps as below: 1) Propagation step: Suppose that at time 𝑘 − 1, the posterior intensity (analogous to pdf) 𝑣(x𝑘−1 ∣Y1:𝑘−1 ) can be approximated in terms of 𝑛𝑥𝑎 𝑘−1 Gaussian distributions of the form 𝑛𝑥𝑎 𝑘−1
𝑣(x𝑘−1 ∣Y1:𝑘−1 ) ≈
∑ 𝑖=1
) ( 𝑥𝑎 𝛼𝑘−1,𝑖 𝑁 x𝑘−1 ; x𝑎𝑘−1,𝑖 , P𝑎𝑥𝑥,𝑘−1,𝑖 , (9)
≈
∑∑ 𝑖=1 𝑗=1
𝑛𝐵 𝑘
=
∑ 𝑠=1
) ( 𝛽 𝑥𝑎 𝐵 𝛼𝑘,𝑗 𝛼𝑘−1,𝑖 𝑁 x𝑘 ; x𝐵 𝑘,(𝑖,𝑗) , P𝑥𝑥,𝑘,(𝑖,𝑗) ,
) ( 𝐵 𝐵 𝛼𝑘,𝑠 𝑁 x𝑘 ; x𝐵 𝑘,𝑠 , P𝑥𝑥,𝑘,𝑠 , (12)
𝐵 Analogous to Eq.(11), x𝐵 𝑘,(𝑖,𝑗) and P𝑥𝑥,𝑘,(𝑖,𝑗) are the mean and covariance of spawn evaluated by propagating forward the 𝑖-th indexed analysis statistics at time 𝑘 − 1, through the 𝑗-th indexed nonlinear/Gaussian system of the spawn x𝑘 = ℳ𝑘,𝑘−1 (x𝑘−1 ) + 𝑑𝛽𝑘−1,𝑗 + 𝑁(u𝑘−1,𝑗 ; 0, Q𝛽𝑘−1,𝑗 ) , 𝑗 ∈ [1, 𝑛𝛽𝑘 ]. 𝛽 𝑥𝑎 𝛽 𝐵 𝑥𝑎 There 𝑛𝐵 𝑘 = 𝑛𝑘 𝑛𝑘−1 , 𝛼𝑘,𝑠 = 𝛼𝑘,𝑗 𝛼𝑘−1,𝑖 with the index 𝑠 to represent the index (𝑖, 𝑗), e.g., 𝑠 = 𝑖 + 𝑛𝑥𝑎 𝑘−1 (𝑗 − 1), 1 ≤ 𝑖 ≤ 𝛽 and 1 ≤ 𝑗 ≤ 𝑛 . 𝑛𝑥𝑎 𝑘−1 𝑘
2) Filtering step: Suppose that the prior (or background) intensity for time 𝑘 is a Gaussian mixture with the form 𝑥𝑏
𝑣(x𝑘 ∣Y1:𝑘−1 ) =
𝑛𝑘 ∑
𝑥𝑏 𝛼𝑘,𝑖 𝑁 (x𝑘 ; x𝑏𝑘,𝑖 , P𝑏𝑥𝑥,𝑘,𝑖 ) ,
(13)
𝑖=1
since the posterior intensity approximately consists of a misdetection term and ∣Y𝑘 ∣ detection terms for each measurement y 𝑘 ∈ Y𝑘 , ∑ 𝑣(x𝑘 ∣Y1:𝑘 ) = (1−𝑝𝐷,𝑘 )𝑣(x𝑘 ∣Y1:𝑘−1 )+ 𝑣𝐷 (x𝑘 ∣Y1:𝑘 ) , y𝑘 ∈Y𝑘
(14) to determine the term 𝑣𝐷,𝑘 (x𝑘 ∣Y1:𝑘 ), according to Bayes’ rule and the new observation y𝑘 ∈ Y𝑘 , we update prior intensity 𝑣(x𝑘 ∣Y1:𝑘−1 ) to posterior 𝑣𝐷,𝑘 (x𝑘 ∣Y1:𝑘 ) by using Eqs.(13)(3b) and yields 𝑣𝐷 (x𝑘 ∣Y1:𝑘 ) ∝ 𝑝𝐷,𝑘 𝑝(y𝑘 ∣x𝑘 )𝑣(x𝑘 ∣Y1:𝑘−1 ) 𝑛𝑥𝑏 𝑛𝑣 𝑘 𝑘
≈ 𝑝𝐷,𝑘
∑∑ 𝑖=1 𝑗=1 𝑥𝑏
= 𝑝𝐷,𝑘 ×
= 𝑝𝐷,𝑘 ×
𝑝𝐷,𝑘 𝒰𝑖,𝑗 ∑𝑛𝑥𝑏 ∑𝑛𝑣𝑘 𝑘
𝐷 𝛼𝑘,𝑠 =
, 𝜅𝑘 (y𝑘 ) + 𝑝𝐷,𝑘 𝑖=1 𝑗=1 𝒰𝑖,𝑗 ) ( 𝑥𝑏 𝑣 = 𝛼𝑘,𝑖 𝛼𝑘,𝑗 𝑁 y𝑘 ; ℋ𝑘 (x𝑏𝑘,𝑖 ), P𝑏𝑦𝑦,𝑘,𝑖 + R𝑘,𝑗 , 𝑛𝐷 𝑘
𝑥𝑏 𝑣 𝛼𝑘,𝑖 𝛼𝑘,𝑗 𝑁
(
y𝑘 ; ℋ𝑘 (x𝑏𝑘,𝑖 ), P𝑏𝑦𝑦,𝑘,𝑖
𝑖=1 𝑗=1 𝑁 (x𝑘 ; x𝑏𝑘,𝑖 , P𝑏𝑥𝑥,𝑘,𝑖 )𝑁 (y𝑘 ; ℋ𝑘 (x𝑘 ), R𝑘,𝑗 ) 𝑁 (y𝑘 ; ℋ𝑘 (x𝑏𝑘,𝑖 ), P𝑏𝑦𝑦,𝑘,𝑖 + R𝑘,𝑗 ) 𝑥𝑏
𝑥𝑏 𝑣 𝑥𝑏 Analogously, if we let 𝑛𝐷 𝑘 = 𝑛𝑘 𝑛𝑘 , 𝑠 = 𝑖 + 𝑛𝑘 (𝑗 − 1), or 𝑣 𝑥𝑏 𝑣 𝑠 = 𝑗 +𝑛𝑘 (𝑖−1), and 1 ≤ ∫ 𝑖 ≤ 𝑛𝑘 , 1 ≤ 𝑗 ≤ 𝑛𝑘 . Then consider the discarded constant 𝑝𝐷,𝑘 𝑝(y𝑘 ∣x𝑘 )𝑣(x𝑘 ∣Y1:𝑘−1 )𝑑x𝑘 by using the symbol “∝” in the first line of Eq.(15), we can compute the normalized weight with index 𝑠 and the corresponding intensity through
) ( 𝑥𝑏 𝑣 𝛼𝑘,𝑖 𝛼𝑘,𝑗 𝑁 x𝑘 ; x𝑏𝑘,𝑖 , P𝑏𝑥𝑥,𝑘,𝑖 𝑁 (y𝑘 ; ℋ𝑘 (x𝑘 ), R𝑘,𝑗 ) ,𝒰𝑖,𝑗
𝑣
𝑛𝑘 𝑛𝑘 ∑ ∑
where a singular value decomposition is applied to the co𝑥 𝑏 variance matrix (H𝑥𝑘 S𝑏𝑥,𝑘,𝑖 )𝑇 R−1 𝑘,𝑗 (H𝑘 S𝑥,𝑘,𝑖 ) (Eq.(17c)), and 2 𝑛 the eigenvalues {𝜎𝑘,𝑖 }𝑖=1 and the eigenvectors {e𝑘,𝑖 }𝑛𝑖=1 are ˆ 𝑘 = diag(𝜎 2 , ⋅ ⋅ ⋅ , 𝜎 2 ) is a sorted in descending order, D 𝑘,1 𝑘,ℓ𝑘 ℓ𝑘 × ℓ𝑘 diagonal matrix formed with the first ℓ𝑘 -th bigger ˆ 𝑘 = [e𝑘,1 , ⋅ ⋅ ⋅ , e𝑘,ℓ ] is a 𝑛 × ℓ𝑘 eigenvectors eigenvalues, V 𝑘 matrix. The definitions of H𝑥𝑘 , S𝑏𝑥,𝑘,𝑖 , Iℓ𝑘 ×ℓ𝑘 , as well as the explaining of approximation by truncation can be seen in [11].
+ R𝑘,𝑗
)
𝑣𝐷 (x𝑘 ∣Y1:𝑘 ) ≈
∑
(18a) (18b)
) ( 𝐷 𝐷 𝛼𝑘,𝑠 𝑁 x𝑘 ; x𝐷 𝑘,𝑠 , P𝑥𝑥,𝑘,𝑠 .
(18c)
𝑠=1
,
3) Statistics abstraction: Then let’s return to Eq.(14). By substituting Eqs.(18c)(13) into Eq.(14) we have
𝑣
𝑛𝑘 𝑛𝑘 ∑ ∑
𝑥𝑎
𝑥𝑏 𝑣 𝛼𝑘,𝑖 𝛼𝑘,𝑗 𝑁 (y𝑘 ; ℋ𝑘 (x𝑏𝑘,𝑖 ), P𝑏𝑦𝑦,𝑘,𝑖 + R𝑘,𝑗 )
𝑣(x𝑘 ∣Y1:𝑘 ) ≈
𝑖=1 𝑗=1 𝐷 𝑁 (x𝑘 : x𝐷 𝑘,(𝑖,𝑗) , P𝑥𝑥,𝑘,(𝑖,𝑗) ) ,
) ( 𝑥𝑎 𝑥𝑎 𝛼𝑘,𝑠 𝑁 x𝑘 ; x𝑥𝑎 𝑘,𝑠 , P𝑥𝑥,𝑘,𝑠 ,
(19)
𝑠=1
(15) P𝑏𝑦𝑦,𝑘,𝑖
where is the projection covariance of the Gaussian random variable with the background mean x𝑏𝑘,𝑖 and covariance P𝑏𝑥𝑥,𝑘,𝑖 , and can be computed by SEVF [11]. The x𝐷 𝑘,(𝑖,𝑗) and are the posterior first-two moments of x with the P𝐷 𝑘 𝑥𝑥,𝑘,(𝑖,𝑗) 𝑖-th indexed prior 𝑁(x𝑘 ; x𝑏𝑘,𝑖 , P𝑏𝑥𝑥,𝑘,𝑖 ) of x𝑘 , and the 𝑗-th indexed observer y𝑘 = ℋ𝑘 (x𝑘 )+𝑁 (v𝑘,𝑗 ; 0, R𝑘,𝑗 ) , 𝑗 ∈ [1, 𝑛𝑣𝑘 ]. Concretely, consider the nonlinear/Gaussian components in Eq.(3), and follow the filtering step of the RSEV filtering pre𝐷 sented in [11], we compute the analysis x𝐷 𝑘,(𝑖,𝑗) and P𝑥𝑥,𝑘,(𝑖,𝑗) by ′ x𝐷 𝑘,(𝑖,𝑗) = arg min 𝐽𝑘 (x) , subject to 𝐽𝑘 (x) ≥ 0 , (16a) x
𝑏 ℓ𝑘 ×ℓ𝑘 −1 ˆ 𝑇 ˆ ˆ P𝐷 ) (V𝑘 ) (S𝑏𝑥,𝑘,𝑖 )𝑇 , 𝑥𝑥,𝑘,(𝑖,𝑗) = S𝑥,𝑘,𝑖 V𝑘 (D𝑘 + I (16b) where 1 𝐽𝑘 (x) = (x − x𝑏𝑘,𝑖 )𝑇 (P𝑏𝑥𝑥,𝑘,𝑖 )−1 (x − x𝑏𝑘,𝑖 ) 2 )𝑇 1( y𝑘 − ℋ𝑘 (x𝑏𝑘,𝑖 ) − H𝑥𝑘 (x − x𝑏𝑘,𝑖 ) + 2 ( ) 𝑏 𝑥 𝑏 × R−1 𝑘,𝑗 y𝑘 − ℋ𝑘 (x𝑘,𝑖 ) − H𝑘 (x − x𝑘,𝑖 ) , (17a)
𝐽𝑘′ (x) = (P𝑏𝑥𝑥,𝑘,𝑖 )−1 (x − x𝑏𝑘,𝑖 ) − (H𝑥𝑘 )𝑇 R−1 𝑘,𝑗 ( ) 𝑏 𝑥 𝑏 × y𝑘 − ℋ𝑘 (x𝑘,𝑖 ) − H𝑘 (x − x𝑘,𝑖 ) , ˆ 𝑘V ˆ 𝑇 ≈ (H𝑥 S𝑏 )𝑇 R−1 (H𝑥 S𝑏 ) , ˆ 𝑘D V 𝑘 𝑘 𝑥,𝑘,𝑖 𝑘 𝑥,𝑘,𝑖 𝑘,𝑗
𝑛𝑘 ∑
where 𝑥𝑏 𝐷 𝑛𝑥𝑎 𝑘 = 𝑛𝑘 + ∣Y𝑘 ∣𝑛𝑘 , } { 𝑥𝑏 𝐷 𝑛𝑥𝑎 𝑥𝑏 𝑛𝑘 𝐷 ∣Y𝑘 ∣𝑛𝑘 𝑘 , {x𝑥𝑎 𝑘,𝑠 }𝑠=1 = {x𝑘,𝑖 }𝑖=1 , {x𝑘,𝑖 }𝑖=1 } { 𝑥𝑎 𝑥𝑏 𝑛𝑘 𝑛𝑘 ∣Y𝑘 ∣𝑛𝐷 𝑥𝑏 𝐷 𝑘 {P𝑥𝑎 , } = {P } , {P } 𝑥𝑥,𝑘,𝑠 𝑠=1 𝑥𝑥,𝑘,𝑖 𝑖=1 𝑥𝑥,𝑘,𝑖 𝑖=1 } { 𝑥𝑎 𝑥𝑏 𝐷 𝑥𝑎 𝑛𝑘 𝑥𝑏 𝑛𝑘 𝐷 ∣Y𝑘 ∣𝑛𝑘 {𝛼𝑘,𝑠 . }𝑠=1 = {(1 − 𝑝𝐷,𝑘 )𝛼𝑘,𝑖 }𝑖=1 , {𝛼𝑘,𝑖 }𝑖=1
(20a) (20b) (20c) (20d)
Given the Gaussian mixture intensities 𝑣𝐷 (x𝑘 ∣Y1:𝑘−1 ) and 𝑣(x𝑘 ∣Y1:𝑘 ), the corresponding the mean of the updated numˆ (x𝑘 ∣Y1:𝑘 ) of tarˆ (x𝑘 ∣Y1:𝑘−1 ) and predicted number 𝑁 ber 𝑁 gets can be obtained by summing up the appropriate weights, 𝑥𝑎
ˆ (x𝑘 ∣Y1:𝑘 ) = (1 − 𝑝𝐷,𝑘 )𝑁 ˆ (x𝑘 ∣Y1:𝑘−1 ) + 𝑁
𝑛𝑘 ∑ ∑
𝑥𝑎 𝛼𝑘,𝑠 ,
y𝑘 ∈Y𝑘 𝑠=1
(21a) ˆ (x𝑘 ∣Y1:𝑘−1 ) 𝑁
⎛
ˆ (x𝑘−1 ∣Y1:𝑘−1 ) ⎝𝑝𝑆,𝑘 =𝑁
𝛽
𝑢
𝑛𝑘 ∑ 𝑗=1
𝑢 𝛼𝑘,𝑗
+
𝑛𝑘 ∑ 𝑗=1
⎞ 𝛽 ⎠ 𝛼𝑘,𝑗
𝐵
+
𝑛𝑘 ∑
𝐵 𝛼𝑘,𝑠 .
𝑠=1
(21b) (17b) (17c)
𝑥𝑎 ˆ (x𝑘−1 ∣Y1:𝑘−1 ) = ∑𝑛𝑘−1 𝛼𝑥𝑎 where 𝑁 𝑖=1 𝑘−1,𝑖 is target number estimated by the last recursion.
C. Outline of RSEV-EGM-PHD with pruning/merging
particularly, the collection of projection covariance, 𝑛𝑥𝑏 𝑛𝑆 𝑛𝐵 𝑘 𝑘 𝑘 , consists of {P𝑆𝑦𝑦,𝑘,𝑠 }𝑠=1 , {P𝐵 {P𝑏𝑦𝑦,𝑘,𝑖 }𝑖=1 𝑦𝑦,𝑘,𝑠 }𝑠=1
This section we focus on the procedures after the RSEVEGM-PHD being initialized. The initialization of the RSEVEGM-PHD will be discussed in IV. We assume that at time 𝑘 − 1, the posterior intensity 𝑣(x𝑘−1 ∣Y1:𝑘−1 ) is approximated by 𝑞 weighted Gaussian sum and show 𝑞 ∑ 𝑣(x𝑘−1 ∣Y1:𝑘−1 ) = 𝑏𝑘−1,𝑠 𝑁 (x𝑘−1 ; 𝒵𝑘−1,𝑠 , Φ𝑘−1 ) . (22) 𝑠=1
𝑛𝛾
𝑘 and {P𝛾𝑦𝑦,𝑘,𝑠 }𝑠=1 .
2) Filtering step: ∙
∙
With the available observation y𝑘 ∈ Y𝑘 , update the prior moment x𝑏𝑘,𝑗 and P𝑏𝑥𝑥,𝑘,𝑗 (𝑗 = 1, ⋅ ⋅ ⋅ , 𝑛𝑥𝑏 𝑘 ) to 𝐷 𝐷 the posterior ones x𝑘,𝑠 and P𝑥𝑥,𝑘,𝑠 (𝑠 = 1, ⋅ ⋅ ⋅ , 𝑛𝐷 𝑘 , 𝑣 𝑛𝐷 = 𝑛𝑥𝑏 𝑘 𝑘 × 𝑛𝑘 ) for each Gaussian distribution 𝑁 (x𝑘 ; x𝑏𝑘,𝑗 , P𝑏𝑥𝑥,𝑘,𝑗 ) in each 𝐻-sub-system according formulaes in Eq.(16).
∙
𝑥𝑏 Update the prior weights 𝛼𝑘,𝑗 ’s (𝑗 = 1, ⋅ ⋅ ⋅ , 𝑛𝑥𝑏 𝑘 ) to 𝐷 the posterior ones 𝛼𝑘,𝑠 ’s (𝑠 = 1, ⋅ ⋅ ⋅ , 𝑛𝐷 𝑘 ) according to Eq.(18).
∙
For the term related with no observations in Eq.(14), the mean and covariance are kept the same with that in 𝑛𝑥𝑏 𝑥𝑏 𝑥𝑏 𝑘 {𝛼𝑘,𝑠 , x𝑏𝑘,𝑠 , P𝑏𝑥𝑥,𝑘,𝑠 }𝑠=1 except the weights 𝛼𝑘,𝑠 being 𝑥𝑏 replaced by (1 − 𝑝𝐷,𝑘 )𝛼𝑘,𝑠 .
∙
Form the posterior intensity in Eq.(14) with
Then for each Gaussian distribution 𝑁 (x𝑘−1 ; 𝒵𝑘−1,𝑠 , Φ𝑘−1 ) (𝑠 = 1, ⋅ ⋅ ⋅ , 𝑞), there is a set of 2ℓ𝑘−1 + 1 sigma points 2ℓ𝑘−1 2ℓ𝑘−1 𝑠 𝑠 }𝑖=0 , with the associated weights {𝑊𝑘−1,𝑖 }𝑖=0 . {𝒳𝑘−1,𝑖 The birth targets at time 𝑘 with intensity approximated by 𝛾
𝛾
𝑣 (x𝑘 ) =
𝑛𝑘 ∑ 𝑟=1
𝛾 𝛼𝑘,𝑟 𝑁 (x𝑘 ; x𝛾𝑘,𝑟 , P𝛾𝑥𝑥,𝑘,𝑟 ) ,
(23)
for each Gaussian distribution 𝑁 (x𝑘 ; x𝛾𝑘,𝑟 , P𝛾𝑥𝑥,𝑘,𝑟 ) (𝑟 = 1, ⋅ ⋅ ⋅ , 𝑛𝛾𝑘 ), there is a set of 2ℓ′𝑘 + 1 pair of weights and sigma ′ 𝑟 𝑟 2ℓ𝑘 points {𝑊𝑘,𝑖 , 𝒳𝑘,𝑖 }𝑖=0 , which are only used for computing the projection covariance P𝛾𝑦𝑦,𝑘,𝑟 . The procedure in the next recursion is outlined as follows: 1) Propagation step: ∙
∙
∙
∑𝑛𝑢𝑘
𝑢 𝑁 (u𝑘 ; 0, Q𝑘,𝑖 ) of Given the pdf 𝑝(u𝑘 ) ≈ 𝑖=1 𝛼𝑘,𝑖 the dynamical noise, divide the original dynamical system into 𝑛𝑢𝑘 𝑆-sub-systems described by Eq.(3a), and 𝑛𝛽𝑘 𝛽-sub-systems described by Eq.(8b) respectively. 2ℓ
𝑘−1 𝑠 Propagate each set of sigma points {𝒳𝑘−1,𝑖 }𝑖=0 forward through the 𝑆-sub-systems, and evaluate the background mean x𝑆𝑘,𝑠 and covariance P𝑆𝑥𝑥,𝑘,𝑠 as well as the projection covariance P𝑆𝑦𝑦,𝑘,𝑠 according to the formulae in [11], with index 𝑠 = 1, ⋅ ⋅ ⋅ , 𝑛𝑆𝑘 , 𝑛𝑆𝑘 = 𝑞 × 𝑛𝑢𝑘 . Then, update the weights 𝑏𝑘−1,𝑖 ’s 𝑆 𝑢 to 𝛼𝑘,𝑠 = 𝛼𝑘,𝑗 𝑏𝑘−1,𝑖 , with 𝑠 = 𝑖 + 𝑞 × (𝑗 − 1), 𝑖 = 1, ⋅ ⋅ ⋅ , 𝑞, 𝑗 = 1, ⋅ ⋅ ⋅ , 𝑛𝑢𝑘 .
Similarly propagate each set of sigma points 2ℓ𝑘−1 𝑠 {𝒳𝑘−1,𝑖 }𝑖=0 forward through the 𝛽-sub-systems, 𝑛𝐵
𝐵 𝐵 𝐵 𝐵 𝑘 and evaluate {𝛼𝑘,𝑠 , x𝐵 𝑘,𝑠 , P𝑥𝑥,𝑘,𝑠 , P𝑦𝑦,𝑘,𝑠 }𝑠=1 , 𝑛𝑘 = 𝛽 𝛽 𝐵 = 𝛼𝑘,𝑗 𝑏𝑘−1,𝑖 , with index 𝑞 ×𝑛𝑘 , and the weights 𝛼𝑘,𝑠 𝑠 = 𝑖 + 𝑞 × (𝑗 − 1), 𝑖 = 1, ⋅ ⋅ ⋅ , 𝑞, 𝑗 = 1, ⋅ ⋅ ⋅ , 𝑛𝛽𝑘 .
∙
For term 𝑣 𝛾 (x𝑘 ), we regard each set of sigma ′ 𝑟 2ℓ𝑘 points {𝒳𝑘,𝑖 }𝑖=0 as the background samples and therefore keep the same mean and covariance for 𝑟-th component. We only use the each set to compute the square root matrix S𝛾𝑦,𝑘,𝑟 , thus obtain P𝛾𝑦𝑦,𝑘,𝑟 = S𝛾𝑦,𝑘,𝑟 (S𝛾𝑦,𝑘,𝑟 )𝑇 . Finally we have 𝑛𝛾
𝛾 𝑘 , x𝛾𝑘,𝑠 , P𝛾𝑥𝑥,𝑘,𝑠 , P𝛾𝑦𝑦,𝑘,𝑠 }𝑠=1 . {𝛼𝑘,𝑠
∙
𝑥𝑏
𝑣(x𝑘 ∣Y1:𝑘−1 ) ≈
𝑠=1
𝑥𝑎
𝑣(x𝑘 ∣Y1:𝑘 ) ≈
𝑛𝑘 ∑
𝑥𝑎 𝛼𝑘,𝑠 𝑁 (x𝑘 : x𝑎𝑘,𝑠 , P𝑎𝑥𝑥,𝑘,𝑠 ) . (25)
𝑠=1
3) Reduction with pruning/merging: To prune/merge the 𝑣(x𝑘 ∣Y1:𝑘 ) by a set of 𝑝 (here 𝑝 < 𝑛𝑥𝑎 𝑘 ) Gaussian distributions 𝑣˜(x𝑘 ∣Y1:𝑘 ) =
𝑝 ∑
𝑏𝑘,ℓ 𝑁 (x𝑘 ; 𝒵𝑘,ℓ , Φ𝑘,ℓ ) ,
(26)
ℓ=1
based on the pruning/merging procedure [7], first should specify parameters such as a truncation threshold 𝑇 𝑡𝑟𝑢𝑛𝑐 , and a merging threshold 𝑈 𝑚𝑒𝑟𝑔𝑒 , then perform reduction recursively. IV.
R ESULTS OF NUMERICAL SIMULATION
This section assesses the performance of the RSEV-EGMPHD filter using numerical simulations. Consider a twodimensional scenario with an unknown and time varying number of targets observed in clutter over the surveillance region [9 × 104 , 2.6 × 105 ] × [−3 × 104 , 1.2 × 105 ] (in 𝑚). The target state takes the augmented form x𝑘 = [z𝑇𝑘 , 𝜔𝑘 ]𝑇 , where z𝑘 = [𝑝𝑥,𝑘 , 𝑝˙𝑥,𝑘 , 𝑝¨𝑥,𝑘 , 𝑝𝑦,𝑘 , 𝑝˙𝑦,𝑘 , 𝑝¨𝑦,𝑘 ]𝑇 of target state consists of position (𝑝𝑥,𝑘 , 𝑝𝑦,𝑘 ), velocity (𝑝˙𝑥,𝑘 , 𝑝˙𝑦,𝑘 ) and acceleration (¨ 𝑝𝑥,𝑘 , 𝑝¨𝑦,𝑘 ). 𝜔𝑘 is the turn-rate. The state dynamics are given by z𝑘 = 𝐹 (𝜔𝑘−1 )z𝑘−1 + u𝑧,𝑘−1 , 𝜔𝑘 = 𝜔𝑘−1 + u𝜔,𝑘−1 ,
(27a) (27b)
where the non-Gaussian process noise u𝑧,𝑘 will be specified later. u𝜔,𝑘 ∼ 𝑁 (⋅; 0, 𝜎𝜔2 ), the standard deviation 𝜎𝜔 = 𝜋/360 rad/s.
Form the prior intensity 𝑛𝑘 ∑
∑𝑛𝑣𝑘 𝑣 Given the pdf 𝑝(v𝑘 ) ≈ 𝑖=1 𝛼𝑘,𝑖 𝑁 (v𝑘 ; 0, R𝑘,𝑖 ) of the observation noise, divide the original observation system into 𝑛𝑣𝑘 𝐻-sub-systems described by Eq.(3b).
𝑥𝑏 𝛼𝑘,𝑠 𝑁 (x𝑘 ; x𝑏𝑘,𝑠 , P𝑏𝑥𝑥,𝑘,𝑠 ) ,
(24)
The IMM filter in [15] is adapted for target tracking and the model set consists of 3 sub-models, the nearly-constant velocity (CV) model, the coordinated-turn (CT) model and
the nearly-constant acceleration (CA) model, the definitions of their transition matrix and covariance of process noise are omitted here, to see details in [16]. The model probability is initialized by 𝑝𝐶𝑉,0 = 𝑝𝐶𝐴,0 = 𝑝𝐶𝑇,0 = 1/3. The transition probability matrix of] sub-models is set as a [ 0.8 0.1 0.1 const by 𝜋𝑘∣𝑘−1 = 0.1 0.8 0.1 . The standard deviation 0.1 0.1 0.8 of noise covariances of CV, CA and CT models is initialized by same value, [𝜎𝐶𝑉,𝑥 , 𝜎𝐶𝑉,𝑦 ]𝑇 = [𝜎𝐶𝐴,𝑥 , 𝜎𝐶𝐴,𝑦 ]𝑇 = [𝜎𝐶𝑇,𝑥 , 𝜎𝐶𝑇,𝑦 ]𝑇 = [0.01, 0.01]𝑇 , the turn rate 𝜔0 = 0.001. The parameters involved by the RSEV filtering in [11] are set with: 𝛼 = 1, 𝜆 = −2 and 𝛽 = 2. Assume that simulation time is 400s and sample time is set as 𝑇 = 5s, so the total simulation steps is 80. There are 4 targets moving in the surveillance region and each one can be described by CV, CT and CA model movement or their combination. The designed trajectories consist of 4 different tracks: Target 1 performs a clock-wise turn flying with the (1) initial state x0 = [78000, 800, 1, 97000, 140, −4, 0]𝑇 , the turn moving follows a CT model movement with a specified turn-rate 𝜔𝑘 = 0.008 during steps [1, 76], then appends with a CV model movement during left steps; Target 2 performs (2) a counter-clock-wise turn flying with initial state x0 = 𝑇 [230000, −800, −1, −23000, 120, 4, 0] , the turn flying also follows a CT movement however with a changed turn-rate 𝜔𝑘 = −0.008 during steps [1, 76], then appends with a CV movement and keeps its moving mode to the end; Target 3 is spawned from target 4 at the step 40 and fulfils a CV model movement with an initial velocity and accelerate [𝑝𝑥,𝑘 , 𝑝˙𝑥,𝑘 , 𝑝𝑦,𝑘 , 𝑝˙ 𝑦,𝑘 ] = [350, 0, 350, 0.1], it disappears at the step 64; Target 4 performs a CA model movement and (4) initialized by x0 = [100000, 400, 0.1, 12000, 150, 0.1, 0]𝑇 , its existing covers steps [16, 80]. The 4 targets have same survival probability 𝑝𝑆,𝑘 = 0.98, detection probability 𝑝𝐷,𝑘 = 0.98. The non-Gaussian process noise u𝑘 in Eq.(2a) or then in [u𝑇𝑧,𝑘 , u𝜔,𝑘 ]𝑇 in Eq.(27), is assumed as a weighted sum of 2 student t-distribution, i.e., [u𝑇𝑧,𝑘 , u𝜔,𝑘 ]𝑇 = 0.2𝜎𝑄1 ,𝑘 𝑇 (100) + 0.8𝜎𝑄2 ,𝑘 𝑇 (120), the standard deviation of noise variances of the 2 components are 𝜎𝑄1 ,𝑘 = diag{1, 0.1, 0.01, 1, 0.1, 0.01, 0.0005} and 𝜎𝑄2 ,𝑘 = diag{1.5, 0.15, 0.015, 1.5, 0.15, 0.015, 0.0005} respectively. The observations modeled by the ℋ𝑘 in Eq.(2b) in polar coordinates consists of range, bearing and doppler measurements, is formulated by ⎤ ⎡ (𝑝2𝑥,𝑘 + 𝑝2𝑦,𝑘 )1/2 ⎦ + v𝑘 atan(𝑝𝑦,𝑘 /𝑝𝑥,𝑘 ) y𝑘 = ⎣ 2 2 1/2 (𝑝𝑥,𝑘 𝑝˙ 𝑥,𝑘 + 𝑝𝑦,𝑘 𝑝˙ 𝑦,𝑘 )/(𝑝𝑥,𝑘 + 𝑝𝑦,𝑘 ) Each component of measurement noise vector v𝑘 is subjected student t-distribution with a free degree 𝜈 = 100 and the standard deviation of noise covariance 𝜎𝑅,𝑘 = diag{𝜎𝑟,𝑘 , 𝜎𝜃,𝑘 , 𝜎𝑟,𝑘 ˙ } = diag{100, 𝜋/9, 1}, that is v𝑘 = 𝜎𝑅,𝑘 𝑇 (100). The Gaussian intensity with a mixture form (2) 𝑣𝑘𝛾 (x) =0.3𝒩 (x : x(1) 𝛾 , P𝛾 ) + 0.3𝒩 (x : x𝛾 , P𝛾 )+
0.4𝒩 (x : x(3) 𝛾 , P𝛾 ) ,
x 10
4
10
Y coordinate (m)
8
target 1: born at k=1 dies at k=80 target 3: born at k=40 dies at k=64
6
4
2
0
target 2: born at k=1 dies at k=80
target 4: born at k=16 dies at k=80
−2 1
1.2
1.4
1.6
1.8
2
2.2
X coordinate (m) Fig. 1.
x 10
5
True trajectories of 4 targets in the surveillance region. (1)
(2)
where x𝛾 = −x𝛾 = [103 , 102 , 1, 103 , 102 , 1, 10−2 ]𝑇 , (3) x𝛾 = [10, 1, 0.1, 10, 1, 0.1, 0.01]𝑇 and P𝛾 = diag([104 , 103 , 102 , 104 , 103 , 102 , 0.1]𝑇 ) is used to approximate the spontaneous births that defined in Eq.(8a) in the (1) (2) (3) vicinity of x𝛾 , x𝛾 and x𝛾 . Similarly, the RFS 𝛽𝑘∣𝑘−1 (𝜁) of target spawned from a existing target with previous state 𝜁 is modeled by a Gaussian and with intensity (1)
(2)
𝑣𝑘𝛽 (x∣𝜁) =0.44𝒩 (x : 𝜁, Q𝛽 ) + 0.22𝒩 (x : 𝜁, Q𝛽 )+ (3)
0.34𝒩 (x : 𝜁, Q𝛽 ) , (1)
(2)
where Q𝛽 = diag{[103 , 102 , 1, 103 , 102 , 1, 10−2 ]𝑇 }, Q𝛽 = (1) (3) (1) 10Q𝛽 and Q𝛽 = 5Q𝛽 . Also, the intensity of existing target is approximated by (1)
(2)
𝑣𝑘𝑆 (x) =0.3𝒩 (x : x𝑆 , Q𝑆 ) + 0.3𝒩 (x : x𝑆 , Q𝑆 )+ (3)
(1)
0.4𝒩 (x : x𝑆 , Q𝑆 ) , (2)
(3)
x𝑆 = − x𝑆 = 2x𝑆 = [102 , 10, 0.1, 102 , 10, 0.1, 10−3 ]𝑇 , Q𝑆 =diag{[10002 , 5002 , 102 , 10002 , 5002 , 102 , 10−2 ]𝑇 } . The detected measurements are immersed in clutter that can be modelled as a Poisson RFS 𝐾𝑘 with intensity 𝜅𝑘 (y) = 𝜆𝑐 𝑉 𝑢(y) , where 𝑢(⋅) is the uniform density over the surveillance region, 𝑉 = 2.55 × 106 m2 is the “volume” of the surveillance region, and 𝜆𝑐 = 3.922 × 106 m−2 is the average number of clutter returns per unit volume (i.e., 10 clutter returns over the surveillance region 𝑉 ). The true target trajectories are plotted in Fig.1, there are totally 4 targets. At time 𝑘 = 1, targets 1 and 2 appear and begin to travel in turn lines with specified turn rate, then switch to a straight line mode at time 𝑘 = 76 and exists to the end. Target 4 appears at time 𝑘 = 16 and travel with straight acceleration mode until time 𝑘 = 80. Target 3 spawns from target 4 at time 𝑘 = 40 and fulfils a straight line mode until time 𝑘 = 64. The real change of target number with time can also be seen from Fig.5.
5
2.4
x 10
2400 Reality tracks Estimate tracks
2200
2.2
2000
D pc , OSPA measure (m)
X coordinate (m)
2
1.8
1.6
1.4
1800 1600 1400 1200 1000
1.2 800 1
600
0.8 0
50
100
150
200
400 0
250
10
20
30
simulation time (s)
Fig. 2.
40
50
60
70
80
time step (k)
X-axis position comparison between estimates and realities.
Fig. 4.
OSPA measure with parameter 𝑐 = 3𝑘𝑚 and 𝑝 = 2.
4
12
x 10
4
Reality tracks Estimate tracks
Reality num Estimate num
3.8
10
3.6
8
target num
Y coordinate (m)
3.4 6
4
2
3.2 3 2.8 2.6
0
2.4 −2 2.2 −4 0
50
100
150
200
2
250
simulation time (s)
Fig. 3.
Y-axis position comparison between estimates and realities.
The pruning/merging parameters are specified by 𝑇 𝑡𝑟𝑢𝑛𝑐 = 10 and 𝑈 𝑚𝑒𝑟𝑔𝑒 = 5. It was shown in simulations that the performance of pruning/merging reduction may not be too sensible when specifying threshold values with a small change within a bound. −4
0
10
20
30
40
50
60
70
80
time step (k)
Fig. 5.
Comparison of the target number estimates and true ones.
𝑑(𝑐) (x, y) = min(𝑐, 𝑑(x, y)) is the distance 𝑑 between x and y cut off at 𝑐, Π𝑌 is the set of all possible permutations of 𝑌 , 𝑑 is the Euclidean distance in this work. Fig.4 shows the OSPA measure (cutoff 𝑐 = 3𝑘𝑚 and 𝑝 = 2) for the RSEV-EGMPHD filtering, we observe that in most steps of simulation, the OSPA keeps lower level and fluctuates around an amplitude 1.4𝑘𝑚, but in certain steps there are several jumps with higher amplitude, accordingly there are total 4 mis-estimated target number appeared during the all 80 steps of simulation, as shown in Fig.5, it indicates the right ratio of estimated target number is as high as 95%.
The estimation results, shown in Figs.2 and 3, indicate that both in x-anis and y-axis directions, the RSEV-EGM-PHD filter exhibits good tracking performance. In order to performs an optimal assignment of target estimates to the real states, we adopt the optimal sub-pattern assignment (OSPA) measure [17] for a further evaluation. Let 𝑌 be the set of true target states and 𝑋 be the set of target estimates, with cardinalities 𝑛 and Theoretically for nonlinear Bayes estimation in practice, 𝑚 respectively, the OSPA measure is defined by the RSEV filtering has been shown better performance than the [ ( )] 1/𝑝 𝑚 ∑ 1 KF-like filters such as EKF, UKF, EnKF (see details in [11]), 𝑐 (𝑐) 𝑝 𝑝 𝐷𝑝 (𝑋, 𝑌 ) = 𝑑 (x𝑖 , y𝜋(𝑖) ) + 𝑐 (𝑛 − 𝑚) , the same is expected in PHD filtering with nonlinear dynamics min 𝑛 𝜋∈Π𝑌 𝑖=1 and observations, the performance improvement between the KF-like based PHD filter and the RSEV based PHD filter if 𝑚 ≤ 𝑛 and 𝐷𝑝𝑐 (𝑋, 𝑌 ) = 𝐷𝑝𝑐 (𝑌, 𝑋) otherwise. Here
surely is explicit and therefore not be focused in this work. Through many efforts being made for estimation improvement, such as the thresholds choice for pruning/merging reduction, the parameters choice for the RSEV filtering, the initial state specified for targets maneuvering, the sub-models initialization for IMM, etc, however some extra errors are still introduced with too many parameters prior in propagation, how to optimally specify the related parameters to match the “gap” between the dynamics with arbitrary distribution of noise and the filtering algorithm, is still an important aspect of PHD-like filters application in practice. V.
C ONCLUSION
The Gaussian mixture PHD (GM-PHD) filter was suggested as a closed-form solutions to the PHD recursion in linear/nonlinear Gaussian multitarget tracking [7], theoretically it can not process the dynamics with non-Gaussian noise, e.g., the nonlinear/non-Gaussian system in Eqs.(2). In this paper, focusing on the multitarget tracking in nonlinear/non-Gaussian system in Eqs.(2), we extend the original GM-PHD through using the weighted Gaussian sum [9, ch.8] to process the dynamics with arbitrary distribution of noise (or non-Gaussian noise). Under the weighted Gaussian sum formed state transition and likelihood in Eqs.(3), and linear Gaussian assumptions, we have derived closed-form recursions for the weights, means, and covariances of the Gaussian components of the posterior intensity. Particularly, multitarget state estimation in a nonlinear/non-Gaussian system can be approximately recast as state estimation in a set of parallel nonlinear/Gaussian systems. Moreover for an improved accuracy of nonlinearity estimation, the reduced rank scaled unscented/ensemble transform variational (RSEV) filtering [10] is applied to each individual nonlinear/Gaussian system to evaluate the mean and covariance of the corresponding Gaussian pdf. Thus the suggested approach, the RSEV-based extended Gaussian mixture PHD (RSEV-EGM-PHD), is termed for being different from the original GM-PHD. In addition, an implementation of the RSEV-EGM-PHD has been proposed by combining the closed-form recursions with the pruning/merging to reduce the growing number of components. Simulations have demonstrated the capability of the RSEV-EGM-PHD to track an unknown and time-varying number of targets under the student distribution noise, detection uncertainty and clutters. One of possible future research on the RSEV-EGM-PHD is the efficiency and simplicity in implementation, which also the requirement of possible application such as to tracking in sensor networks. ACKNOWLEDGMENTS This research was supported by the Chinese Natural Science Foundation (granted No.61271317 and No.61175028). R EFERENCES [1] Y. Bar-Shalom and T. E. Fortmann, Tracking and Data Association. San Diego, CA: Academic, 1988. [2] L. D. Stone, C. A. Barlow, and T. L. Corwin, Bayesian Multiple Target Tracking. Boston: Artech House, 1999.
[3] Mahler, R., Random set theory for target tracking and identification. In D. Hall and J. Llinas (Eds)., Data Fusion Handbook,Boca Raton, FL: CRC Press, 2001. [4] Vo, B-N., S. Singh and A. Doucet, Sequential Monte Carlo implementation of the PHD filter for multitarget tracking, Proc. 6th Int. Conf. Information Fusion, Vol.2,July 2003, pp.792-799. [5] R. Mahler, Multi-target Bayes filtering via first-order multi-target moments, IEEE Trans. Aerosp. Electron. Syst., vol.39, no.4, pp.1152C1178, 2003. [6] Vo, B-N., and Ma, W. K. A closed-form solution for the probability hypothesis density filter. In Proceedings of the 7th International Conference on Information Fusion, vol.2, July 2005, 25 28. [7] Vo, B-N., Ma W K. The Gaussian mixture probability hypothesis density filter. IEEE Transactions on Signal Processing, 54(11):40914104,2006. [8] Vo, B-N., S. Singh, and A. Doucet,, Sequential Monte Carlo methods for multi-target filtering with random finite sets, IEEE Trans. Aerosp. Electron. Syst., vol. 41, no. 4, pp. 1224 C1245, 2005. [9] B. D. O. Anderson, J. B. Moore, Optimal Filtering, Prentice-Hall, 1979. [10] Ming Lei, Zhongliang Jing and Shiqiang Hu, Scaled Unscented Transform-based Variational Optimality Filter, 15th International Conference on Information Fusion. Singapore, 9-12th, July 2012. [11] Ming Lei, Christophe Baehr, Unscented/Ensemble Transform-based Variational Filter, Physica D: Nonlinear Phenomena, Vol.246, No.1, March 2013. [12] H. W. Sorenson and D. L. Alspach, Recursive Bayesian estimation using Gaussian sum, Automatica, vol.7, pp.465-479, 1971. [13] D. L. Alspach and H. W. Sorenson, Nonlinear Bayesian estimation using Gaussian sum approximations, IEEE Trans. Autom. Control, vol. AC17, no. 4, pp. 439C448, 1972. [14] S. J. Julier, The scaled unscented transformation, The Proceedings of the American Control Conference, Anchorage, AK, 2004. [15] Li, X.R.; Jilkov, V.P., Survey of maneuvering target tracking. Part V. Multiple-model methods,IEEE Trans. on Aerospace alectronic Systems, vol.41, no.4, pp.1255-1321, 2005. [16] X.Rong Li, Survey of Maneuvering Target Tracking. Part1:Dynamic Models, IEEE Trans. on Aerospace and Electronic Systems,Vol.39, no.4,pp: 1337-1364, 2003. [17] Schuhmacher, Dominic, B-T. Vo, and B-N. Vo., A consistent metric for performance evaluation of multi-object filters, IEEE Trans. on Signal Processing, vol.56, no.8, pp.3447-3457,2008.
I.
I NTRODUCTION
Multitarget tracking involves the joint estimation of a timevarying number of targets as well as their individual states from a sequence of sets of noisy measurements/clutters with uncertain data association [1]. Analogous to the single-target Bayesian filtering for being an optimal framework for singletarget state estimation, by using the random finite set (RFS) to represent the collections of targets and measurements, the finite set statistics (FISST) was adopted to construct an extended optimal Bayesian framework for multitarget tracking [2], [3]. The optimal multitarget Bayes filter under the framework of FISST, requires to evaluate combinatorial sums of multiple integrals of high dimensionality with a prohibitively large number of combinations, thus the approach is computationally intractable [3], [4]. The probability hypothesis density (PHD) filter is proposed as a suboptimal but computationally tractable alternative to the RFS-based Bayes multitarget estimation [3], [5]. Aim at the implementation of PHD filtering, a Gaussian mixture PHD recursion, called GM-PHD filter, has been suggested for jointly estimating the time-varying number of targets and corresponding state RFS from the observation RFS in linear/nonlinear Gaussian multitarget system [6], [7]. Based on the Gaussian assumption of state noise and observation noise, the prior PHD is approximated by a sum of weighted Gaussian components, and with the weights, means, and covariances being updated separately by the Kalman filter/Extended Kalman filter/Unscented Kalman filter, the posterior PHD is proved also in form of a weighted Gaussian sum and therefore has an approximate close-form solution to the original problem.
The GM-PHD filter theoretically can not processes the dynamics with non-Gaussian noise, e.g., the nonlinear/nonGaussian system in Eqs.(2), through some other approximations such as the SMC implementations [7], [8] being suggested, however shown that they are still suffer from intensive computation due to the combinatorial nature of the densities, especially when the number of targets is large. In this paper, focusing on the multitarget tracking in nonGaussian system, as described in Eqs.(2), we approximately solve the problem through the weighted Gaussian sum, or Gaussian mixture model [9, ch.8], here use the extended Gaussian mixture PHD (EGM-PHD) for being different from the original GM-PHD. It is shown that for the dynamics and observation both with arbitrary distribution, the prior and posterior PHD still can be formulated by the weighted Gaussian sum and propagating the Gaussian mixtures separately over time. In this way, state estimation in a nonlinear/non-Gaussian system can be approximately recast as state estimation in a set of parallel nonlinear/Gaussian systems. In order to achieve an improved accuracy of nonlinear estimation, the reduced rank scaled unscented/ensemble transform variational (RSEV) filtering in [10], [11] can be applied to each individual nonlinear/Gaussian system to evaluate the mean and covariance of the corresponding Gaussian pdf. Concretely, consider that the pdfs 𝑝 (u𝑘 ) and 𝑝 (v𝑘 ) are of the Gaussian mixture model of the non-Gaussian dynamical and observation noise as shown in Eq.(2), the pdfs of state transition and likelihood at time 𝑘, i.e., 𝑝(x𝑘 ∣x𝑘−1 ) and 𝑝(y𝑘 ∣x𝑘 ), can be approximated by 𝑛𝑢𝑘 and 𝑛𝑣𝑘 Gaussian distributions, see Eqs.(3a)(3b) for details. In addition, like the weighted Gaussian sum presented in [12], [13], when the specified number of Gaussian terms of dynamics and observation noise more than 1 (however a common choice is more than 1 of the terms for the pdf approximation), the computation problems of increasing Gaussian terms as time progresses will be more severe than ever. To reduce the computational cost, [6], [7] suggested a pruning procedure that “discarding those with weights below some preset threshold, or by keeping only a certain number of components with the strongest weights”, and merging that “some of the Gaussian components are so close together that they could be accurately approximated by a single Gaussian.” We also adopt this scheme for computation reduction. For performance illustration, we demonstrate the ability of proposed EGM-PHD to estimate the correct number of targets and their trajectories in high-density clutter. The performance
of the new EGM-PHD filter is benchmarked against the reality specification and its efficacy is shown in terms of its ability to pick up the correct number tracks and the improved accuracy of the trajectory estimates. The rest of this paper is organized as follows. Section II provides a summary of multitarget models with arbitrary noise distributions, and PHD filter in the recursive Bayes equations. In section III, we then proceed to consider the multitarget state estimation problem in nonlinear/non-Gaussian systems. Based on the framework of weighted Gaussian sum and the reduced rank scaled unscented/Ensemble transform variational filtering, we derive the reduced rank scaled unscented/ensemble transform variational-based extended Gaussian mixture PHD (RSEV-EGM-PHD) filter as an extension version to the original GM-PHD filter. For convenience, we also outline the major procedures in the RSEV-EGM-PHD equipped with pruning/merging reduction. An example is then given in section IV to illustrate the details in implementing the RSEV-EGM-PHD. Finally we conclude this paper in section V. II.
P ROBLEM F ORMULATION
Let’s consider the state estimation problem in multitarget scenario. The orderless collections of target states and measurements at time 𝑘, can be represented as finite sets, 𝑛𝑥 ×𝑁𝑘 𝑘 x𝑘 = {x𝑘,𝑖 }𝑁 , 𝑖=1 ∈ ℝ
(1a)
y𝑘 =
(1b)
𝑘 {y𝑘,𝑖 }𝑀 𝑖=1
𝑛𝑦 ×𝑀𝑘
∈ℝ
.
where 𝑁𝑘 and 𝑀𝑘 are unknown and variable with time. Similar to that in single-target filtering, the uncertainty is represented by random vector, the originals uncertainty in a multitarget system is characterized by modelling the multitarget state x𝑘 and multiple target measurement y𝑘 as random finite sets (RFS). An RFS x𝑘 is a finite-set-valued random variable. Then the dynamics can be modeled by x𝑘 = ℳ𝑘,𝑘−1 (x𝑘−1 ) + u𝑘−1 , y𝑘 = ℋ𝑘 (x𝑘 ) + v𝑘 ,
(2a) (2b)
where the transition operator ℳ𝑘,𝑘−1 and the observation operator ℋ𝑘 are both possibly nonlinear. The dynamical and observation noise, in terms of u𝑘−1 and v𝑘 respectively, are non-Gaussian and their pdfs, 𝑝(u ) and 𝑝(v𝑘 ), are known and ∑𝑛𝑢𝑘 𝑘𝑢 approximately 𝑝(u𝑘 ) ≈ 𝛼 𝑖=1 𝑘,𝑖 𝑁 (u𝑘 ; 0, Q𝑘,𝑖 ), 𝑝(v𝑘 ) ≈ ∑𝑛𝑣𝑘 𝑣 𝑖=1 𝛼𝑘,𝑖 𝑁 (v𝑘 ; 0, R𝑘,𝑖 ). Then pdf of transition and likelihood, i.e., 𝑝(x𝑘 ∣x𝑘−1 ) and 𝑝(y𝑘 ∣x𝑘 ), can also be approximated by Gaussian’s with same number of terms, such that 𝑛𝑢 𝑘
𝑝(x𝑘 ∣x𝑘−1 ) ≈
∑
𝑢 𝛼𝑘,𝑖 𝑁 (x𝑘 ; ℳ𝑘,𝑘−1 (x𝑘−1 ) , Q𝑘−1,𝑖 ) ,
𝑖=1 𝑛𝑣 𝑘
𝑝(y𝑘 ∣x𝑘 ) ≈
∑
(3a) 𝑣 𝛼𝑘,𝑖 𝑁 (y𝑘 ; ℋ𝑘 (x𝑘 ) , R𝑘,𝑖 ) ,
(3b)
𝑖=1 𝑢 denotes the weight associated where 𝛼𝑘,𝑖 ∑𝑛𝑢𝑘 𝑢 𝑢 satisfies 𝛼𝑘,𝑖 ∈ [0, 1] and 𝑖=1 𝛼𝑘,𝑖 = 1. 𝑣 definite of 𝛼𝑘,𝑖 .
with Eq.(3a) and Similarly for the
A. Multitarget filtering represented by random finite sets The multitarget state x𝑘 , including target motion, birth, spawning, can be described by RFS. For target motion, given multitarget state RFS x𝑘−1 , each target with x𝑘−1,𝑖 ∈ x𝑘−1 survives at time 𝑘 with surviving probability 𝑝𝑆,𝑘 (x𝑘−1 ), its transition pdf is specified by 𝑝(x𝑘 ∣x𝑘−1 ) in Eq.(3) and therefore the target motion is modeled as the RFS 𝑆(x𝑘 ∣x𝑘−1 ). Similarly, the RFS of target birth at time 𝑘 is modeled by Γ𝑘 , the RFS of targets spawning from a target with x𝑘−1 is modeled by 𝐵(x𝑘 ∣x𝑘−1 ). Summarily the RFS of multitarget state x𝑘 is consist of ⎡ ⎤ ⎡ ⎤ ∪ ∪ ∪ ∪ 𝑋𝑘 = ⎣ 𝑆(x𝑘 ∣x𝑘−1 )⎦ ⎣ 𝐵(x𝑘 ∣x𝑘−1 )⎦ Γ𝑘 . x𝑘−1
x𝑘−1
(4) For the RFS measurement model, suppose that a given target x𝑘 ∈ ℝ𝑛𝑥 ×𝑁𝑘 is either detected with probability 𝑝𝐷,𝑘 (x𝑘 ) or missed with probability 1−𝑝𝐷,𝑘 (x𝑘 ), conditional on detection, the pdf of observed measurement y𝑘 from x𝑘 is modeled by 𝑝(y𝑘 ∣x𝑘 ) in Eq.(3b), then count on either {y𝑘 , ⋅ ⋅ ⋅ } when the target is detected or {0} when the target is not detected, an RFS of the obtained measurements that generated by each state x𝑘 is Θ𝑘 (x𝑘 ). Meanwhile, a set 𝐾𝑘 of false measurements or clutter, may be companied with the target originated measurements. If some of the observations y𝑘,𝑖 ’s are due to clutter, it follows a clutter probability density 𝑐𝑘 , and the number of clutter points are assumed to be Poisson distributed with a mean of 𝜆𝑐 . Thus, given a multitarget state x𝑘 , multitarget measurement y𝑘 can be formed by an union of target generated measurements and clutter, ] [ ∪ ∪ Θ𝑘 (x𝑘 ) . (5) 𝑌 𝑘 = 𝐾𝑘 x𝑘
Then RFSs of multitarget evolution and observation described in Eq.(4) and (5), can be formulated in multitarget transition density 𝑝(x𝑘 ∣x𝑘−1 ) and likelihood 𝑝(y𝑘 ∣x𝑘 ) [5]. B. The probability hypothesis density (PHD) filter PHD filtering propagates a first-order statistics of the posterior multitarget state by appropriate assumptions [1], [5]. Particularly, let 𝑣(x𝑘 ∣Y1:𝑘 ) and 𝑣(x𝑘 ∣Y1:𝑘−1 ) be the respective intensities associated with the multitarget posterior pdf 𝑝(x𝑘 ∣Y1:𝑘 ) and the multitarget predicted pdf 𝑝(x𝑘 ∣Y1:𝑘−1 ), by using the finite set statistics (FISST), PHD recursion can be propagated in time via PHD prediction and PHD update [5]: 𝑣(x𝑘 ∣Y1:𝑘−1 ) ∫ = 𝑝𝑆,𝑘 (x𝑘 )𝑝(x𝑘 ∣x𝑘−1 )𝑣(x𝑘−1 ∣Y1:𝑘−1 )𝑑x𝑘−1 + ∫ 𝑣 𝛾 (x𝑘 ) + 𝑣 𝛽 (x𝑘 ∣x𝑘−1 )𝑣(x𝑘−1 ∣Y1:𝑘−1 )𝑑x𝑘−1 ,
(6a)
𝑣(x𝑘 ∣Y1:𝑘 ) = [1 − 𝑝𝐷,𝑘 (x𝑘 )]𝑣(x𝑘 ∣Y1:𝑘−1 )+ ∑ 𝑝𝐷,𝑘 (x𝑘 )𝑝(y𝑘 ∣x𝑘 )𝑣(x𝑘 ∣Y1:𝑘−1 ) ∫ , 𝜅𝑘 (y𝑘 ) + 𝑝𝐷,𝑘 (x𝑘 )𝑝(y𝑘 ∣x𝑘 )𝑣(x𝑘 ∣Y1:𝑘−1 )𝑑x𝑘 y𝑘 ∈Y𝑘
(6b)
where Y1:𝑘 = {y1 , ⋅ ⋅ ⋅ , y𝑘 } is historical observations up to and including time 𝑘. the ∫posterior intensity 𝑣(x𝑘 ∣Y1:𝑘 ) is the density whose integral 𝑆 𝑣(x𝑘 ∣Y1:𝑘 ))𝑑x𝑘 on any region 𝑆 of RFS gives the expected number of elements of 𝑋 that ∫ ˆ𝑘 = ∣𝑋𝑘 ∩ 𝑆∣𝜇𝑠 (𝑑𝑋𝑘 ). 𝑣 𝛽 (x𝑘 ∣x𝑘−1 ) and are in 𝑆, i.e., 𝑁 𝑣 𝛾 (x𝑘 ) denotes the intensities of 𝐵(x𝑘 ∣x𝑘−1 ) and Γ𝑘 at time 𝑘 in Eq.(4), 𝜅𝑘 (y𝑘 ) is the intensity of the clutter RFS 𝐾𝑘 in Eq.(5) and equals 𝜆𝑐 𝑐𝑘 . 𝑝𝑆,𝑘 (⋅) and 𝑝𝐷,𝑘 (⋅) are the survival and detection probabilities.
∑𝑛𝑥𝑎 𝑘−1 𝑥𝑎 𝑥𝑎 ∈ [0, 1] and 𝑖=1 𝛼𝑘−1,𝑖 = 1. Then the prior where 𝛼𝑘−1,𝑖 intensity 𝑣(x𝑘 ∣Y1:𝑘−1 ) to time 𝑘 is also approximately a Gaussian mixture, and accordingly which consists of three terms 𝑣 𝑆 (x𝑘 ∣Y1:𝑘−1 ), 𝑣 𝛽 (x𝑘 ∣Y1:𝑘−1 ) and 𝑣 𝛾 (x𝑘 ) respectively, to the existing targets, the spawned targets and the spontaneous births, such that 𝑣(x𝑘 ∣Y1:𝑘−1 ) = 𝑣 𝑆 (x𝑘 ∣Y1:𝑘−1 ) + 𝑣 𝛽 (x𝑘 ∣Y1:𝑘−1 ) + 𝑣 𝛾 (x𝑘 ) , (10) where
III.
RSEV-EGM-PHD FILTERING FOR NONLINEAR / NON -G AUSSIAN MULTITARGET TRACKING
𝑢 𝑛𝑥𝑎 𝑘−1 𝑛𝑘
A. Nonlinear/non-Gaussian multitarget model Nonlinear/non-Gaussian multitarget model includes certain assumptions on the birth, death, and detection of targets, are summarized: (1) the survival and detection probabilities are both state independent, i.e., 𝑝𝑆,𝑘 (x𝑘 ) = 𝑝𝑆,𝑘 , 𝑝𝐷,𝑘 (x𝑘 ) = 𝑝𝐷,𝑘 .
𝑣 𝑆 (x𝑘 ∣Y1:𝑘−1 ) = 𝑝𝑆,𝑘
(7a) (7b)
(2) the intensities of the birth and spawn RFSs are Gaussian mixtures of the form
= 𝑝𝑆,𝑘
∑∑
∫ 𝑝(x𝑘 ∣x𝑘−1 )𝑣(x𝑘−1 ∣Y1:𝑘−1 )𝑑x𝑘−1 ,
∫ 𝑢 𝑥𝑎 𝛼𝑘,𝑗 𝛼𝑘−1,𝑖 𝑁 (x𝑘 ; ℳ𝑘,𝑘−1 (x𝑘−1 ) , Q𝑘−1,𝑗 ) ×
𝑖=1 𝑗=1
) ( 𝑁 x𝑘−1 ; x𝑎𝑘−1,𝑖 , P𝑎𝑥𝑥,𝑘−1,𝑖 𝑑x𝑘−1 , 𝑢 𝑛𝑥𝑎 𝑘−1 𝑛𝑘
≈ 𝑝𝑆,𝑘
∑∑ 𝑖=1 𝑗=1 𝑛𝑆 𝑘
= 𝑝𝑆,𝑘
∑
) ( 𝑢 𝑥𝑎 𝛼𝑘,𝑗 𝛼𝑘−1,𝑖 𝑁 x𝑘 ; x𝑆𝑘,(𝑖,𝑗) , P𝑆𝑥𝑥,𝑘,(𝑖,𝑗) ,
) ( 𝑆 𝛼𝑘,𝑠 𝑁 x𝑘 ; x𝑆𝑘,𝑠 , P𝑆𝑥𝑥,𝑘,𝑠 ,
𝑠=1
(11)
𝛾
𝑛𝑘 ∑
𝛾 𝛼𝑘,𝑖 𝑁 (x𝑘 ; x𝛾𝑘,𝑖 , P𝛾𝑥𝑥,𝑘,𝑖 ) ,
x𝑆𝑘,(𝑖,𝑗)
P𝑆𝑥𝑥,𝑘,(𝑖,𝑗)
where and the mean and covariance of the background evaluated by SUT [14], that is, to propagate 𝑖=1 forward the 𝑖-th indexed analysis statistics (i.e., x𝑎𝑘−1,𝑖 and 𝛽 𝑛𝑘 𝑎 𝑥𝑎 ∑ , here 𝑖 ∈ [1, 𝑛𝑘−1 ] ) at time 𝑘 − 1, through the 𝑗P 𝛽 𝛼𝑘,𝑖 𝑁 (x𝑘 ; ℳ𝑘,𝑘−1 (x𝑘−1 ) + 𝑑𝛽𝑘−1,𝑖 , Q𝛽𝑘−1,𝑖 ) , th𝑥𝑥,𝑘−1,𝑖 𝑣 𝛽 (x𝑘 ∣x𝑘−1 ) = indexed nonlinear/Gaussian system x𝑘 = ℳ𝑘,𝑘−1 (x𝑘−1 ) + 𝑖=1 𝑁(u𝑘−1,𝑗 ; 0, Q𝑘−1,𝑗 ) , 𝑗 ∈ [1, 𝑛𝑢𝑘 ], then to estimate the first(8b) two moment of transferred statistics. 𝛾 , x𝛾𝑘,𝑖 , P𝛾𝑥𝑥,𝑘,𝑖 , 𝑖 = 1, ⋅ ⋅ ⋅ , 𝑛𝛾𝑘 are given model where 𝑛𝛾𝑘 , 𝛼𝑘,𝑖 𝑆 𝑢 𝑥𝑎 Moreover, 𝑛𝑆𝑘 = 𝑛𝑢𝑘 𝑛𝑥𝑎 𝑘−1 , 𝛼𝑘,𝑠 = 𝛼𝑘,𝑗 𝛼𝑘−1,𝑖 with the parameters that determine the shape of the birth intensity; integer index 𝑠 being a representation of the two-dimensional 𝛽 𝑥𝑎 similarly 𝑛𝛽𝑘 , 𝛼𝑘,𝑖 , ℳ𝑘,𝑘−1 , 𝑑𝛽𝑘−1,𝑖 , Q𝛽𝑘−1,𝑖 , 𝑖 = 1, ⋅ ⋅ ⋅ , 𝑛𝛽𝑘 deindex (𝑖, 𝑗), e.g., 𝑠 = 𝑖 + 𝑛𝑥𝑎 𝑘−1 (𝑗 − 1), 1 ≤ 𝑖 ≤ 𝑛𝑘−1 and termine the shape of the spawning intersity of a target with 1 ≤ 𝑗 ≤ 𝑛𝑢𝑘 . Similarly the previous state x𝑘−1 . ∫ 𝛽 𝛾 𝛾 𝑣 (x ∣Y ) = 𝑣 𝛽 (x𝑘 ∣x𝑘−1 )𝑣(x𝑘−1 ∣Y1:𝑘−1 )𝑑x𝑘−1 , 𝑘 1:𝑘−1 In the second assumption, x𝑘,𝑖 , 𝑖 = 1, ⋅ ⋅ ⋅ , 𝑛𝑘 are the peaks of the spontaneous birth intensity in (8a). These points 𝛽 𝑛𝑥𝑎 ∫ 𝑘−1 𝑛𝑘 ∑ ∑ 𝛽 have the highest local concentrations of expected number 𝛽 𝛽 𝑥𝑎 𝛾 = 𝛼 𝛼 𝑘,𝑗 𝑘−1,𝑖 𝑁 (x𝑘 ; ℳ𝑘,𝑘−1 (x𝑘−1 ) + 𝑑𝑘−1,𝑗 , Q𝑘−1,𝑗 ) of spontaneous births and represent. The weight 𝛼𝑘,𝑖 given 𝑖=1 𝑗=1 the expected number of new targets originating from x𝛾𝑘,𝑖 . × 𝑁 (x𝑘−1 ; x𝑎𝑘−1,𝑖 , P𝑎𝑥𝑥,𝑘−1,𝑖 )𝑑x𝑘−1 , Analogously for the spawning intensity in Eq.(8b). 𝑣 𝛾 (x𝑘 ) =
(8a)
𝛽 𝑛𝑥𝑎 𝑘−1 𝑛𝑘
B. RSEV-based extended GM-PHD (RSEV-EGM-PHD) Filter The multitarget model described with Gaussian mixture model for nonlinear/non-Gaussian system given in III-A, the procedure of RSEV-based extended GM-PHD filtering of Eqs.(6a)(6b) yields a approximated solution called RSEV-EGMPHD recursion, we split the recursion into 3 steps as below: 1) Propagation step: Suppose that at time 𝑘 − 1, the posterior intensity (analogous to pdf) 𝑣(x𝑘−1 ∣Y1:𝑘−1 ) can be approximated in terms of 𝑛𝑥𝑎 𝑘−1 Gaussian distributions of the form 𝑛𝑥𝑎 𝑘−1
𝑣(x𝑘−1 ∣Y1:𝑘−1 ) ≈
∑ 𝑖=1
) ( 𝑥𝑎 𝛼𝑘−1,𝑖 𝑁 x𝑘−1 ; x𝑎𝑘−1,𝑖 , P𝑎𝑥𝑥,𝑘−1,𝑖 , (9)
≈
∑∑ 𝑖=1 𝑗=1
𝑛𝐵 𝑘
=
∑ 𝑠=1
) ( 𝛽 𝑥𝑎 𝐵 𝛼𝑘,𝑗 𝛼𝑘−1,𝑖 𝑁 x𝑘 ; x𝐵 𝑘,(𝑖,𝑗) , P𝑥𝑥,𝑘,(𝑖,𝑗) ,
) ( 𝐵 𝐵 𝛼𝑘,𝑠 𝑁 x𝑘 ; x𝐵 𝑘,𝑠 , P𝑥𝑥,𝑘,𝑠 , (12)
𝐵 Analogous to Eq.(11), x𝐵 𝑘,(𝑖,𝑗) and P𝑥𝑥,𝑘,(𝑖,𝑗) are the mean and covariance of spawn evaluated by propagating forward the 𝑖-th indexed analysis statistics at time 𝑘 − 1, through the 𝑗-th indexed nonlinear/Gaussian system of the spawn x𝑘 = ℳ𝑘,𝑘−1 (x𝑘−1 ) + 𝑑𝛽𝑘−1,𝑗 + 𝑁(u𝑘−1,𝑗 ; 0, Q𝛽𝑘−1,𝑗 ) , 𝑗 ∈ [1, 𝑛𝛽𝑘 ]. 𝛽 𝑥𝑎 𝛽 𝐵 𝑥𝑎 There 𝑛𝐵 𝑘 = 𝑛𝑘 𝑛𝑘−1 , 𝛼𝑘,𝑠 = 𝛼𝑘,𝑗 𝛼𝑘−1,𝑖 with the index 𝑠 to represent the index (𝑖, 𝑗), e.g., 𝑠 = 𝑖 + 𝑛𝑥𝑎 𝑘−1 (𝑗 − 1), 1 ≤ 𝑖 ≤ 𝛽 and 1 ≤ 𝑗 ≤ 𝑛 . 𝑛𝑥𝑎 𝑘−1 𝑘
2) Filtering step: Suppose that the prior (or background) intensity for time 𝑘 is a Gaussian mixture with the form 𝑥𝑏
𝑣(x𝑘 ∣Y1:𝑘−1 ) =
𝑛𝑘 ∑
𝑥𝑏 𝛼𝑘,𝑖 𝑁 (x𝑘 ; x𝑏𝑘,𝑖 , P𝑏𝑥𝑥,𝑘,𝑖 ) ,
(13)
𝑖=1
since the posterior intensity approximately consists of a misdetection term and ∣Y𝑘 ∣ detection terms for each measurement y 𝑘 ∈ Y𝑘 , ∑ 𝑣(x𝑘 ∣Y1:𝑘 ) = (1−𝑝𝐷,𝑘 )𝑣(x𝑘 ∣Y1:𝑘−1 )+ 𝑣𝐷 (x𝑘 ∣Y1:𝑘 ) , y𝑘 ∈Y𝑘
(14) to determine the term 𝑣𝐷,𝑘 (x𝑘 ∣Y1:𝑘 ), according to Bayes’ rule and the new observation y𝑘 ∈ Y𝑘 , we update prior intensity 𝑣(x𝑘 ∣Y1:𝑘−1 ) to posterior 𝑣𝐷,𝑘 (x𝑘 ∣Y1:𝑘 ) by using Eqs.(13)(3b) and yields 𝑣𝐷 (x𝑘 ∣Y1:𝑘 ) ∝ 𝑝𝐷,𝑘 𝑝(y𝑘 ∣x𝑘 )𝑣(x𝑘 ∣Y1:𝑘−1 ) 𝑛𝑥𝑏 𝑛𝑣 𝑘 𝑘
≈ 𝑝𝐷,𝑘
∑∑ 𝑖=1 𝑗=1 𝑥𝑏
= 𝑝𝐷,𝑘 ×
= 𝑝𝐷,𝑘 ×
𝑝𝐷,𝑘 𝒰𝑖,𝑗 ∑𝑛𝑥𝑏 ∑𝑛𝑣𝑘 𝑘
𝐷 𝛼𝑘,𝑠 =
, 𝜅𝑘 (y𝑘 ) + 𝑝𝐷,𝑘 𝑖=1 𝑗=1 𝒰𝑖,𝑗 ) ( 𝑥𝑏 𝑣 = 𝛼𝑘,𝑖 𝛼𝑘,𝑗 𝑁 y𝑘 ; ℋ𝑘 (x𝑏𝑘,𝑖 ), P𝑏𝑦𝑦,𝑘,𝑖 + R𝑘,𝑗 , 𝑛𝐷 𝑘
𝑥𝑏 𝑣 𝛼𝑘,𝑖 𝛼𝑘,𝑗 𝑁
(
y𝑘 ; ℋ𝑘 (x𝑏𝑘,𝑖 ), P𝑏𝑦𝑦,𝑘,𝑖
𝑖=1 𝑗=1 𝑁 (x𝑘 ; x𝑏𝑘,𝑖 , P𝑏𝑥𝑥,𝑘,𝑖 )𝑁 (y𝑘 ; ℋ𝑘 (x𝑘 ), R𝑘,𝑗 ) 𝑁 (y𝑘 ; ℋ𝑘 (x𝑏𝑘,𝑖 ), P𝑏𝑦𝑦,𝑘,𝑖 + R𝑘,𝑗 ) 𝑥𝑏
𝑥𝑏 𝑣 𝑥𝑏 Analogously, if we let 𝑛𝐷 𝑘 = 𝑛𝑘 𝑛𝑘 , 𝑠 = 𝑖 + 𝑛𝑘 (𝑗 − 1), or 𝑣 𝑥𝑏 𝑣 𝑠 = 𝑗 +𝑛𝑘 (𝑖−1), and 1 ≤ ∫ 𝑖 ≤ 𝑛𝑘 , 1 ≤ 𝑗 ≤ 𝑛𝑘 . Then consider the discarded constant 𝑝𝐷,𝑘 𝑝(y𝑘 ∣x𝑘 )𝑣(x𝑘 ∣Y1:𝑘−1 )𝑑x𝑘 by using the symbol “∝” in the first line of Eq.(15), we can compute the normalized weight with index 𝑠 and the corresponding intensity through
) ( 𝑥𝑏 𝑣 𝛼𝑘,𝑖 𝛼𝑘,𝑗 𝑁 x𝑘 ; x𝑏𝑘,𝑖 , P𝑏𝑥𝑥,𝑘,𝑖 𝑁 (y𝑘 ; ℋ𝑘 (x𝑘 ), R𝑘,𝑗 ) ,𝒰𝑖,𝑗
𝑣
𝑛𝑘 𝑛𝑘 ∑ ∑
where a singular value decomposition is applied to the co𝑥 𝑏 variance matrix (H𝑥𝑘 S𝑏𝑥,𝑘,𝑖 )𝑇 R−1 𝑘,𝑗 (H𝑘 S𝑥,𝑘,𝑖 ) (Eq.(17c)), and 2 𝑛 the eigenvalues {𝜎𝑘,𝑖 }𝑖=1 and the eigenvectors {e𝑘,𝑖 }𝑛𝑖=1 are ˆ 𝑘 = diag(𝜎 2 , ⋅ ⋅ ⋅ , 𝜎 2 ) is a sorted in descending order, D 𝑘,1 𝑘,ℓ𝑘 ℓ𝑘 × ℓ𝑘 diagonal matrix formed with the first ℓ𝑘 -th bigger ˆ 𝑘 = [e𝑘,1 , ⋅ ⋅ ⋅ , e𝑘,ℓ ] is a 𝑛 × ℓ𝑘 eigenvectors eigenvalues, V 𝑘 matrix. The definitions of H𝑥𝑘 , S𝑏𝑥,𝑘,𝑖 , Iℓ𝑘 ×ℓ𝑘 , as well as the explaining of approximation by truncation can be seen in [11].
+ R𝑘,𝑗
)
𝑣𝐷 (x𝑘 ∣Y1:𝑘 ) ≈
∑
(18a) (18b)
) ( 𝐷 𝐷 𝛼𝑘,𝑠 𝑁 x𝑘 ; x𝐷 𝑘,𝑠 , P𝑥𝑥,𝑘,𝑠 .
(18c)
𝑠=1
,
3) Statistics abstraction: Then let’s return to Eq.(14). By substituting Eqs.(18c)(13) into Eq.(14) we have
𝑣
𝑛𝑘 𝑛𝑘 ∑ ∑
𝑥𝑎
𝑥𝑏 𝑣 𝛼𝑘,𝑖 𝛼𝑘,𝑗 𝑁 (y𝑘 ; ℋ𝑘 (x𝑏𝑘,𝑖 ), P𝑏𝑦𝑦,𝑘,𝑖 + R𝑘,𝑗 )
𝑣(x𝑘 ∣Y1:𝑘 ) ≈
𝑖=1 𝑗=1 𝐷 𝑁 (x𝑘 : x𝐷 𝑘,(𝑖,𝑗) , P𝑥𝑥,𝑘,(𝑖,𝑗) ) ,
) ( 𝑥𝑎 𝑥𝑎 𝛼𝑘,𝑠 𝑁 x𝑘 ; x𝑥𝑎 𝑘,𝑠 , P𝑥𝑥,𝑘,𝑠 ,
(19)
𝑠=1
(15) P𝑏𝑦𝑦,𝑘,𝑖
where is the projection covariance of the Gaussian random variable with the background mean x𝑏𝑘,𝑖 and covariance P𝑏𝑥𝑥,𝑘,𝑖 , and can be computed by SEVF [11]. The x𝐷 𝑘,(𝑖,𝑗) and are the posterior first-two moments of x with the P𝐷 𝑘 𝑥𝑥,𝑘,(𝑖,𝑗) 𝑖-th indexed prior 𝑁(x𝑘 ; x𝑏𝑘,𝑖 , P𝑏𝑥𝑥,𝑘,𝑖 ) of x𝑘 , and the 𝑗-th indexed observer y𝑘 = ℋ𝑘 (x𝑘 )+𝑁 (v𝑘,𝑗 ; 0, R𝑘,𝑗 ) , 𝑗 ∈ [1, 𝑛𝑣𝑘 ]. Concretely, consider the nonlinear/Gaussian components in Eq.(3), and follow the filtering step of the RSEV filtering pre𝐷 sented in [11], we compute the analysis x𝐷 𝑘,(𝑖,𝑗) and P𝑥𝑥,𝑘,(𝑖,𝑗) by ′ x𝐷 𝑘,(𝑖,𝑗) = arg min 𝐽𝑘 (x) , subject to 𝐽𝑘 (x) ≥ 0 , (16a) x
𝑏 ℓ𝑘 ×ℓ𝑘 −1 ˆ 𝑇 ˆ ˆ P𝐷 ) (V𝑘 ) (S𝑏𝑥,𝑘,𝑖 )𝑇 , 𝑥𝑥,𝑘,(𝑖,𝑗) = S𝑥,𝑘,𝑖 V𝑘 (D𝑘 + I (16b) where 1 𝐽𝑘 (x) = (x − x𝑏𝑘,𝑖 )𝑇 (P𝑏𝑥𝑥,𝑘,𝑖 )−1 (x − x𝑏𝑘,𝑖 ) 2 )𝑇 1( y𝑘 − ℋ𝑘 (x𝑏𝑘,𝑖 ) − H𝑥𝑘 (x − x𝑏𝑘,𝑖 ) + 2 ( ) 𝑏 𝑥 𝑏 × R−1 𝑘,𝑗 y𝑘 − ℋ𝑘 (x𝑘,𝑖 ) − H𝑘 (x − x𝑘,𝑖 ) , (17a)
𝐽𝑘′ (x) = (P𝑏𝑥𝑥,𝑘,𝑖 )−1 (x − x𝑏𝑘,𝑖 ) − (H𝑥𝑘 )𝑇 R−1 𝑘,𝑗 ( ) 𝑏 𝑥 𝑏 × y𝑘 − ℋ𝑘 (x𝑘,𝑖 ) − H𝑘 (x − x𝑘,𝑖 ) , ˆ 𝑘V ˆ 𝑇 ≈ (H𝑥 S𝑏 )𝑇 R−1 (H𝑥 S𝑏 ) , ˆ 𝑘D V 𝑘 𝑘 𝑥,𝑘,𝑖 𝑘 𝑥,𝑘,𝑖 𝑘,𝑗
𝑛𝑘 ∑
where 𝑥𝑏 𝐷 𝑛𝑥𝑎 𝑘 = 𝑛𝑘 + ∣Y𝑘 ∣𝑛𝑘 , } { 𝑥𝑏 𝐷 𝑛𝑥𝑎 𝑥𝑏 𝑛𝑘 𝐷 ∣Y𝑘 ∣𝑛𝑘 𝑘 , {x𝑥𝑎 𝑘,𝑠 }𝑠=1 = {x𝑘,𝑖 }𝑖=1 , {x𝑘,𝑖 }𝑖=1 } { 𝑥𝑎 𝑥𝑏 𝑛𝑘 𝑛𝑘 ∣Y𝑘 ∣𝑛𝐷 𝑥𝑏 𝐷 𝑘 {P𝑥𝑎 , } = {P } , {P } 𝑥𝑥,𝑘,𝑠 𝑠=1 𝑥𝑥,𝑘,𝑖 𝑖=1 𝑥𝑥,𝑘,𝑖 𝑖=1 } { 𝑥𝑎 𝑥𝑏 𝐷 𝑥𝑎 𝑛𝑘 𝑥𝑏 𝑛𝑘 𝐷 ∣Y𝑘 ∣𝑛𝑘 {𝛼𝑘,𝑠 . }𝑠=1 = {(1 − 𝑝𝐷,𝑘 )𝛼𝑘,𝑖 }𝑖=1 , {𝛼𝑘,𝑖 }𝑖=1
(20a) (20b) (20c) (20d)
Given the Gaussian mixture intensities 𝑣𝐷 (x𝑘 ∣Y1:𝑘−1 ) and 𝑣(x𝑘 ∣Y1:𝑘 ), the corresponding the mean of the updated numˆ (x𝑘 ∣Y1:𝑘 ) of tarˆ (x𝑘 ∣Y1:𝑘−1 ) and predicted number 𝑁 ber 𝑁 gets can be obtained by summing up the appropriate weights, 𝑥𝑎
ˆ (x𝑘 ∣Y1:𝑘 ) = (1 − 𝑝𝐷,𝑘 )𝑁 ˆ (x𝑘 ∣Y1:𝑘−1 ) + 𝑁
𝑛𝑘 ∑ ∑
𝑥𝑎 𝛼𝑘,𝑠 ,
y𝑘 ∈Y𝑘 𝑠=1
(21a) ˆ (x𝑘 ∣Y1:𝑘−1 ) 𝑁
⎛
ˆ (x𝑘−1 ∣Y1:𝑘−1 ) ⎝𝑝𝑆,𝑘 =𝑁
𝛽
𝑢
𝑛𝑘 ∑ 𝑗=1
𝑢 𝛼𝑘,𝑗
+
𝑛𝑘 ∑ 𝑗=1
⎞ 𝛽 ⎠ 𝛼𝑘,𝑗
𝐵
+
𝑛𝑘 ∑
𝐵 𝛼𝑘,𝑠 .
𝑠=1
(21b) (17b) (17c)
𝑥𝑎 ˆ (x𝑘−1 ∣Y1:𝑘−1 ) = ∑𝑛𝑘−1 𝛼𝑥𝑎 where 𝑁 𝑖=1 𝑘−1,𝑖 is target number estimated by the last recursion.
C. Outline of RSEV-EGM-PHD with pruning/merging
particularly, the collection of projection covariance, 𝑛𝑥𝑏 𝑛𝑆 𝑛𝐵 𝑘 𝑘 𝑘 , consists of {P𝑆𝑦𝑦,𝑘,𝑠 }𝑠=1 , {P𝐵 {P𝑏𝑦𝑦,𝑘,𝑖 }𝑖=1 𝑦𝑦,𝑘,𝑠 }𝑠=1
This section we focus on the procedures after the RSEVEGM-PHD being initialized. The initialization of the RSEVEGM-PHD will be discussed in IV. We assume that at time 𝑘 − 1, the posterior intensity 𝑣(x𝑘−1 ∣Y1:𝑘−1 ) is approximated by 𝑞 weighted Gaussian sum and show 𝑞 ∑ 𝑣(x𝑘−1 ∣Y1:𝑘−1 ) = 𝑏𝑘−1,𝑠 𝑁 (x𝑘−1 ; 𝒵𝑘−1,𝑠 , Φ𝑘−1 ) . (22) 𝑠=1
𝑛𝛾
𝑘 and {P𝛾𝑦𝑦,𝑘,𝑠 }𝑠=1 .
2) Filtering step: ∙
∙
With the available observation y𝑘 ∈ Y𝑘 , update the prior moment x𝑏𝑘,𝑗 and P𝑏𝑥𝑥,𝑘,𝑗 (𝑗 = 1, ⋅ ⋅ ⋅ , 𝑛𝑥𝑏 𝑘 ) to 𝐷 𝐷 the posterior ones x𝑘,𝑠 and P𝑥𝑥,𝑘,𝑠 (𝑠 = 1, ⋅ ⋅ ⋅ , 𝑛𝐷 𝑘 , 𝑣 𝑛𝐷 = 𝑛𝑥𝑏 𝑘 𝑘 × 𝑛𝑘 ) for each Gaussian distribution 𝑁 (x𝑘 ; x𝑏𝑘,𝑗 , P𝑏𝑥𝑥,𝑘,𝑗 ) in each 𝐻-sub-system according formulaes in Eq.(16).
∙
𝑥𝑏 Update the prior weights 𝛼𝑘,𝑗 ’s (𝑗 = 1, ⋅ ⋅ ⋅ , 𝑛𝑥𝑏 𝑘 ) to 𝐷 the posterior ones 𝛼𝑘,𝑠 ’s (𝑠 = 1, ⋅ ⋅ ⋅ , 𝑛𝐷 𝑘 ) according to Eq.(18).
∙
For the term related with no observations in Eq.(14), the mean and covariance are kept the same with that in 𝑛𝑥𝑏 𝑥𝑏 𝑥𝑏 𝑘 {𝛼𝑘,𝑠 , x𝑏𝑘,𝑠 , P𝑏𝑥𝑥,𝑘,𝑠 }𝑠=1 except the weights 𝛼𝑘,𝑠 being 𝑥𝑏 replaced by (1 − 𝑝𝐷,𝑘 )𝛼𝑘,𝑠 .
∙
Form the posterior intensity in Eq.(14) with
Then for each Gaussian distribution 𝑁 (x𝑘−1 ; 𝒵𝑘−1,𝑠 , Φ𝑘−1 ) (𝑠 = 1, ⋅ ⋅ ⋅ , 𝑞), there is a set of 2ℓ𝑘−1 + 1 sigma points 2ℓ𝑘−1 2ℓ𝑘−1 𝑠 𝑠 }𝑖=0 , with the associated weights {𝑊𝑘−1,𝑖 }𝑖=0 . {𝒳𝑘−1,𝑖 The birth targets at time 𝑘 with intensity approximated by 𝛾
𝛾
𝑣 (x𝑘 ) =
𝑛𝑘 ∑ 𝑟=1
𝛾 𝛼𝑘,𝑟 𝑁 (x𝑘 ; x𝛾𝑘,𝑟 , P𝛾𝑥𝑥,𝑘,𝑟 ) ,
(23)
for each Gaussian distribution 𝑁 (x𝑘 ; x𝛾𝑘,𝑟 , P𝛾𝑥𝑥,𝑘,𝑟 ) (𝑟 = 1, ⋅ ⋅ ⋅ , 𝑛𝛾𝑘 ), there is a set of 2ℓ′𝑘 + 1 pair of weights and sigma ′ 𝑟 𝑟 2ℓ𝑘 points {𝑊𝑘,𝑖 , 𝒳𝑘,𝑖 }𝑖=0 , which are only used for computing the projection covariance P𝛾𝑦𝑦,𝑘,𝑟 . The procedure in the next recursion is outlined as follows: 1) Propagation step: ∙
∙
∙
∑𝑛𝑢𝑘
𝑢 𝑁 (u𝑘 ; 0, Q𝑘,𝑖 ) of Given the pdf 𝑝(u𝑘 ) ≈ 𝑖=1 𝛼𝑘,𝑖 the dynamical noise, divide the original dynamical system into 𝑛𝑢𝑘 𝑆-sub-systems described by Eq.(3a), and 𝑛𝛽𝑘 𝛽-sub-systems described by Eq.(8b) respectively. 2ℓ
𝑘−1 𝑠 Propagate each set of sigma points {𝒳𝑘−1,𝑖 }𝑖=0 forward through the 𝑆-sub-systems, and evaluate the background mean x𝑆𝑘,𝑠 and covariance P𝑆𝑥𝑥,𝑘,𝑠 as well as the projection covariance P𝑆𝑦𝑦,𝑘,𝑠 according to the formulae in [11], with index 𝑠 = 1, ⋅ ⋅ ⋅ , 𝑛𝑆𝑘 , 𝑛𝑆𝑘 = 𝑞 × 𝑛𝑢𝑘 . Then, update the weights 𝑏𝑘−1,𝑖 ’s 𝑆 𝑢 to 𝛼𝑘,𝑠 = 𝛼𝑘,𝑗 𝑏𝑘−1,𝑖 , with 𝑠 = 𝑖 + 𝑞 × (𝑗 − 1), 𝑖 = 1, ⋅ ⋅ ⋅ , 𝑞, 𝑗 = 1, ⋅ ⋅ ⋅ , 𝑛𝑢𝑘 .
Similarly propagate each set of sigma points 2ℓ𝑘−1 𝑠 {𝒳𝑘−1,𝑖 }𝑖=0 forward through the 𝛽-sub-systems, 𝑛𝐵
𝐵 𝐵 𝐵 𝐵 𝑘 and evaluate {𝛼𝑘,𝑠 , x𝐵 𝑘,𝑠 , P𝑥𝑥,𝑘,𝑠 , P𝑦𝑦,𝑘,𝑠 }𝑠=1 , 𝑛𝑘 = 𝛽 𝛽 𝐵 = 𝛼𝑘,𝑗 𝑏𝑘−1,𝑖 , with index 𝑞 ×𝑛𝑘 , and the weights 𝛼𝑘,𝑠 𝑠 = 𝑖 + 𝑞 × (𝑗 − 1), 𝑖 = 1, ⋅ ⋅ ⋅ , 𝑞, 𝑗 = 1, ⋅ ⋅ ⋅ , 𝑛𝛽𝑘 .
∙
For term 𝑣 𝛾 (x𝑘 ), we regard each set of sigma ′ 𝑟 2ℓ𝑘 points {𝒳𝑘,𝑖 }𝑖=0 as the background samples and therefore keep the same mean and covariance for 𝑟-th component. We only use the each set to compute the square root matrix S𝛾𝑦,𝑘,𝑟 , thus obtain P𝛾𝑦𝑦,𝑘,𝑟 = S𝛾𝑦,𝑘,𝑟 (S𝛾𝑦,𝑘,𝑟 )𝑇 . Finally we have 𝑛𝛾
𝛾 𝑘 , x𝛾𝑘,𝑠 , P𝛾𝑥𝑥,𝑘,𝑠 , P𝛾𝑦𝑦,𝑘,𝑠 }𝑠=1 . {𝛼𝑘,𝑠
∙
𝑥𝑏
𝑣(x𝑘 ∣Y1:𝑘−1 ) ≈
𝑠=1
𝑥𝑎
𝑣(x𝑘 ∣Y1:𝑘 ) ≈
𝑛𝑘 ∑
𝑥𝑎 𝛼𝑘,𝑠 𝑁 (x𝑘 : x𝑎𝑘,𝑠 , P𝑎𝑥𝑥,𝑘,𝑠 ) . (25)
𝑠=1
3) Reduction with pruning/merging: To prune/merge the 𝑣(x𝑘 ∣Y1:𝑘 ) by a set of 𝑝 (here 𝑝 < 𝑛𝑥𝑎 𝑘 ) Gaussian distributions 𝑣˜(x𝑘 ∣Y1:𝑘 ) =
𝑝 ∑
𝑏𝑘,ℓ 𝑁 (x𝑘 ; 𝒵𝑘,ℓ , Φ𝑘,ℓ ) ,
(26)
ℓ=1
based on the pruning/merging procedure [7], first should specify parameters such as a truncation threshold 𝑇 𝑡𝑟𝑢𝑛𝑐 , and a merging threshold 𝑈 𝑚𝑒𝑟𝑔𝑒 , then perform reduction recursively. IV.
R ESULTS OF NUMERICAL SIMULATION
This section assesses the performance of the RSEV-EGMPHD filter using numerical simulations. Consider a twodimensional scenario with an unknown and time varying number of targets observed in clutter over the surveillance region [9 × 104 , 2.6 × 105 ] × [−3 × 104 , 1.2 × 105 ] (in 𝑚). The target state takes the augmented form x𝑘 = [z𝑇𝑘 , 𝜔𝑘 ]𝑇 , where z𝑘 = [𝑝𝑥,𝑘 , 𝑝˙𝑥,𝑘 , 𝑝¨𝑥,𝑘 , 𝑝𝑦,𝑘 , 𝑝˙𝑦,𝑘 , 𝑝¨𝑦,𝑘 ]𝑇 of target state consists of position (𝑝𝑥,𝑘 , 𝑝𝑦,𝑘 ), velocity (𝑝˙𝑥,𝑘 , 𝑝˙𝑦,𝑘 ) and acceleration (¨ 𝑝𝑥,𝑘 , 𝑝¨𝑦,𝑘 ). 𝜔𝑘 is the turn-rate. The state dynamics are given by z𝑘 = 𝐹 (𝜔𝑘−1 )z𝑘−1 + u𝑧,𝑘−1 , 𝜔𝑘 = 𝜔𝑘−1 + u𝜔,𝑘−1 ,
(27a) (27b)
where the non-Gaussian process noise u𝑧,𝑘 will be specified later. u𝜔,𝑘 ∼ 𝑁 (⋅; 0, 𝜎𝜔2 ), the standard deviation 𝜎𝜔 = 𝜋/360 rad/s.
Form the prior intensity 𝑛𝑘 ∑
∑𝑛𝑣𝑘 𝑣 Given the pdf 𝑝(v𝑘 ) ≈ 𝑖=1 𝛼𝑘,𝑖 𝑁 (v𝑘 ; 0, R𝑘,𝑖 ) of the observation noise, divide the original observation system into 𝑛𝑣𝑘 𝐻-sub-systems described by Eq.(3b).
𝑥𝑏 𝛼𝑘,𝑠 𝑁 (x𝑘 ; x𝑏𝑘,𝑠 , P𝑏𝑥𝑥,𝑘,𝑠 ) ,
(24)
The IMM filter in [15] is adapted for target tracking and the model set consists of 3 sub-models, the nearly-constant velocity (CV) model, the coordinated-turn (CT) model and
the nearly-constant acceleration (CA) model, the definitions of their transition matrix and covariance of process noise are omitted here, to see details in [16]. The model probability is initialized by 𝑝𝐶𝑉,0 = 𝑝𝐶𝐴,0 = 𝑝𝐶𝑇,0 = 1/3. The transition probability matrix of] sub-models is set as a [ 0.8 0.1 0.1 const by 𝜋𝑘∣𝑘−1 = 0.1 0.8 0.1 . The standard deviation 0.1 0.1 0.8 of noise covariances of CV, CA and CT models is initialized by same value, [𝜎𝐶𝑉,𝑥 , 𝜎𝐶𝑉,𝑦 ]𝑇 = [𝜎𝐶𝐴,𝑥 , 𝜎𝐶𝐴,𝑦 ]𝑇 = [𝜎𝐶𝑇,𝑥 , 𝜎𝐶𝑇,𝑦 ]𝑇 = [0.01, 0.01]𝑇 , the turn rate 𝜔0 = 0.001. The parameters involved by the RSEV filtering in [11] are set with: 𝛼 = 1, 𝜆 = −2 and 𝛽 = 2. Assume that simulation time is 400s and sample time is set as 𝑇 = 5s, so the total simulation steps is 80. There are 4 targets moving in the surveillance region and each one can be described by CV, CT and CA model movement or their combination. The designed trajectories consist of 4 different tracks: Target 1 performs a clock-wise turn flying with the (1) initial state x0 = [78000, 800, 1, 97000, 140, −4, 0]𝑇 , the turn moving follows a CT model movement with a specified turn-rate 𝜔𝑘 = 0.008 during steps [1, 76], then appends with a CV model movement during left steps; Target 2 performs (2) a counter-clock-wise turn flying with initial state x0 = 𝑇 [230000, −800, −1, −23000, 120, 4, 0] , the turn flying also follows a CT movement however with a changed turn-rate 𝜔𝑘 = −0.008 during steps [1, 76], then appends with a CV movement and keeps its moving mode to the end; Target 3 is spawned from target 4 at the step 40 and fulfils a CV model movement with an initial velocity and accelerate [𝑝𝑥,𝑘 , 𝑝˙𝑥,𝑘 , 𝑝𝑦,𝑘 , 𝑝˙ 𝑦,𝑘 ] = [350, 0, 350, 0.1], it disappears at the step 64; Target 4 performs a CA model movement and (4) initialized by x0 = [100000, 400, 0.1, 12000, 150, 0.1, 0]𝑇 , its existing covers steps [16, 80]. The 4 targets have same survival probability 𝑝𝑆,𝑘 = 0.98, detection probability 𝑝𝐷,𝑘 = 0.98. The non-Gaussian process noise u𝑘 in Eq.(2a) or then in [u𝑇𝑧,𝑘 , u𝜔,𝑘 ]𝑇 in Eq.(27), is assumed as a weighted sum of 2 student t-distribution, i.e., [u𝑇𝑧,𝑘 , u𝜔,𝑘 ]𝑇 = 0.2𝜎𝑄1 ,𝑘 𝑇 (100) + 0.8𝜎𝑄2 ,𝑘 𝑇 (120), the standard deviation of noise variances of the 2 components are 𝜎𝑄1 ,𝑘 = diag{1, 0.1, 0.01, 1, 0.1, 0.01, 0.0005} and 𝜎𝑄2 ,𝑘 = diag{1.5, 0.15, 0.015, 1.5, 0.15, 0.015, 0.0005} respectively. The observations modeled by the ℋ𝑘 in Eq.(2b) in polar coordinates consists of range, bearing and doppler measurements, is formulated by ⎤ ⎡ (𝑝2𝑥,𝑘 + 𝑝2𝑦,𝑘 )1/2 ⎦ + v𝑘 atan(𝑝𝑦,𝑘 /𝑝𝑥,𝑘 ) y𝑘 = ⎣ 2 2 1/2 (𝑝𝑥,𝑘 𝑝˙ 𝑥,𝑘 + 𝑝𝑦,𝑘 𝑝˙ 𝑦,𝑘 )/(𝑝𝑥,𝑘 + 𝑝𝑦,𝑘 ) Each component of measurement noise vector v𝑘 is subjected student t-distribution with a free degree 𝜈 = 100 and the standard deviation of noise covariance 𝜎𝑅,𝑘 = diag{𝜎𝑟,𝑘 , 𝜎𝜃,𝑘 , 𝜎𝑟,𝑘 ˙ } = diag{100, 𝜋/9, 1}, that is v𝑘 = 𝜎𝑅,𝑘 𝑇 (100). The Gaussian intensity with a mixture form (2) 𝑣𝑘𝛾 (x) =0.3𝒩 (x : x(1) 𝛾 , P𝛾 ) + 0.3𝒩 (x : x𝛾 , P𝛾 )+
0.4𝒩 (x : x(3) 𝛾 , P𝛾 ) ,
x 10
4
10
Y coordinate (m)
8
target 1: born at k=1 dies at k=80 target 3: born at k=40 dies at k=64
6
4
2
0
target 2: born at k=1 dies at k=80
target 4: born at k=16 dies at k=80
−2 1
1.2
1.4
1.6
1.8
2
2.2
X coordinate (m) Fig. 1.
x 10
5
True trajectories of 4 targets in the surveillance region. (1)
(2)
where x𝛾 = −x𝛾 = [103 , 102 , 1, 103 , 102 , 1, 10−2 ]𝑇 , (3) x𝛾 = [10, 1, 0.1, 10, 1, 0.1, 0.01]𝑇 and P𝛾 = diag([104 , 103 , 102 , 104 , 103 , 102 , 0.1]𝑇 ) is used to approximate the spontaneous births that defined in Eq.(8a) in the (1) (2) (3) vicinity of x𝛾 , x𝛾 and x𝛾 . Similarly, the RFS 𝛽𝑘∣𝑘−1 (𝜁) of target spawned from a existing target with previous state 𝜁 is modeled by a Gaussian and with intensity (1)
(2)
𝑣𝑘𝛽 (x∣𝜁) =0.44𝒩 (x : 𝜁, Q𝛽 ) + 0.22𝒩 (x : 𝜁, Q𝛽 )+ (3)
0.34𝒩 (x : 𝜁, Q𝛽 ) , (1)
(2)
where Q𝛽 = diag{[103 , 102 , 1, 103 , 102 , 1, 10−2 ]𝑇 }, Q𝛽 = (1) (3) (1) 10Q𝛽 and Q𝛽 = 5Q𝛽 . Also, the intensity of existing target is approximated by (1)
(2)
𝑣𝑘𝑆 (x) =0.3𝒩 (x : x𝑆 , Q𝑆 ) + 0.3𝒩 (x : x𝑆 , Q𝑆 )+ (3)
(1)
0.4𝒩 (x : x𝑆 , Q𝑆 ) , (2)
(3)
x𝑆 = − x𝑆 = 2x𝑆 = [102 , 10, 0.1, 102 , 10, 0.1, 10−3 ]𝑇 , Q𝑆 =diag{[10002 , 5002 , 102 , 10002 , 5002 , 102 , 10−2 ]𝑇 } . The detected measurements are immersed in clutter that can be modelled as a Poisson RFS 𝐾𝑘 with intensity 𝜅𝑘 (y) = 𝜆𝑐 𝑉 𝑢(y) , where 𝑢(⋅) is the uniform density over the surveillance region, 𝑉 = 2.55 × 106 m2 is the “volume” of the surveillance region, and 𝜆𝑐 = 3.922 × 106 m−2 is the average number of clutter returns per unit volume (i.e., 10 clutter returns over the surveillance region 𝑉 ). The true target trajectories are plotted in Fig.1, there are totally 4 targets. At time 𝑘 = 1, targets 1 and 2 appear and begin to travel in turn lines with specified turn rate, then switch to a straight line mode at time 𝑘 = 76 and exists to the end. Target 4 appears at time 𝑘 = 16 and travel with straight acceleration mode until time 𝑘 = 80. Target 3 spawns from target 4 at time 𝑘 = 40 and fulfils a straight line mode until time 𝑘 = 64. The real change of target number with time can also be seen from Fig.5.
5
2.4
x 10
2400 Reality tracks Estimate tracks
2200
2.2
2000
D pc , OSPA measure (m)
X coordinate (m)
2
1.8
1.6
1.4
1800 1600 1400 1200 1000
1.2 800 1
600
0.8 0
50
100
150
200
400 0
250
10
20
30
simulation time (s)
Fig. 2.
40
50
60
70
80
time step (k)
X-axis position comparison between estimates and realities.
Fig. 4.
OSPA measure with parameter 𝑐 = 3𝑘𝑚 and 𝑝 = 2.
4
12
x 10
4
Reality tracks Estimate tracks
Reality num Estimate num
3.8
10
3.6
8
target num
Y coordinate (m)
3.4 6
4
2
3.2 3 2.8 2.6
0
2.4 −2 2.2 −4 0
50
100
150
200
2
250
simulation time (s)
Fig. 3.
Y-axis position comparison between estimates and realities.
The pruning/merging parameters are specified by 𝑇 𝑡𝑟𝑢𝑛𝑐 = 10 and 𝑈 𝑚𝑒𝑟𝑔𝑒 = 5. It was shown in simulations that the performance of pruning/merging reduction may not be too sensible when specifying threshold values with a small change within a bound. −4
0
10
20
30
40
50
60
70
80
time step (k)
Fig. 5.
Comparison of the target number estimates and true ones.
𝑑(𝑐) (x, y) = min(𝑐, 𝑑(x, y)) is the distance 𝑑 between x and y cut off at 𝑐, Π𝑌 is the set of all possible permutations of 𝑌 , 𝑑 is the Euclidean distance in this work. Fig.4 shows the OSPA measure (cutoff 𝑐 = 3𝑘𝑚 and 𝑝 = 2) for the RSEV-EGMPHD filtering, we observe that in most steps of simulation, the OSPA keeps lower level and fluctuates around an amplitude 1.4𝑘𝑚, but in certain steps there are several jumps with higher amplitude, accordingly there are total 4 mis-estimated target number appeared during the all 80 steps of simulation, as shown in Fig.5, it indicates the right ratio of estimated target number is as high as 95%.
The estimation results, shown in Figs.2 and 3, indicate that both in x-anis and y-axis directions, the RSEV-EGM-PHD filter exhibits good tracking performance. In order to performs an optimal assignment of target estimates to the real states, we adopt the optimal sub-pattern assignment (OSPA) measure [17] for a further evaluation. Let 𝑌 be the set of true target states and 𝑋 be the set of target estimates, with cardinalities 𝑛 and Theoretically for nonlinear Bayes estimation in practice, 𝑚 respectively, the OSPA measure is defined by the RSEV filtering has been shown better performance than the [ ( )] 1/𝑝 𝑚 ∑ 1 KF-like filters such as EKF, UKF, EnKF (see details in [11]), 𝑐 (𝑐) 𝑝 𝑝 𝐷𝑝 (𝑋, 𝑌 ) = 𝑑 (x𝑖 , y𝜋(𝑖) ) + 𝑐 (𝑛 − 𝑚) , the same is expected in PHD filtering with nonlinear dynamics min 𝑛 𝜋∈Π𝑌 𝑖=1 and observations, the performance improvement between the KF-like based PHD filter and the RSEV based PHD filter if 𝑚 ≤ 𝑛 and 𝐷𝑝𝑐 (𝑋, 𝑌 ) = 𝐷𝑝𝑐 (𝑌, 𝑋) otherwise. Here
surely is explicit and therefore not be focused in this work. Through many efforts being made for estimation improvement, such as the thresholds choice for pruning/merging reduction, the parameters choice for the RSEV filtering, the initial state specified for targets maneuvering, the sub-models initialization for IMM, etc, however some extra errors are still introduced with too many parameters prior in propagation, how to optimally specify the related parameters to match the “gap” between the dynamics with arbitrary distribution of noise and the filtering algorithm, is still an important aspect of PHD-like filters application in practice. V.
C ONCLUSION
The Gaussian mixture PHD (GM-PHD) filter was suggested as a closed-form solutions to the PHD recursion in linear/nonlinear Gaussian multitarget tracking [7], theoretically it can not process the dynamics with non-Gaussian noise, e.g., the nonlinear/non-Gaussian system in Eqs.(2). In this paper, focusing on the multitarget tracking in nonlinear/non-Gaussian system in Eqs.(2), we extend the original GM-PHD through using the weighted Gaussian sum [9, ch.8] to process the dynamics with arbitrary distribution of noise (or non-Gaussian noise). Under the weighted Gaussian sum formed state transition and likelihood in Eqs.(3), and linear Gaussian assumptions, we have derived closed-form recursions for the weights, means, and covariances of the Gaussian components of the posterior intensity. Particularly, multitarget state estimation in a nonlinear/non-Gaussian system can be approximately recast as state estimation in a set of parallel nonlinear/Gaussian systems. Moreover for an improved accuracy of nonlinearity estimation, the reduced rank scaled unscented/ensemble transform variational (RSEV) filtering [10] is applied to each individual nonlinear/Gaussian system to evaluate the mean and covariance of the corresponding Gaussian pdf. Thus the suggested approach, the RSEV-based extended Gaussian mixture PHD (RSEV-EGM-PHD), is termed for being different from the original GM-PHD. In addition, an implementation of the RSEV-EGM-PHD has been proposed by combining the closed-form recursions with the pruning/merging to reduce the growing number of components. Simulations have demonstrated the capability of the RSEV-EGM-PHD to track an unknown and time-varying number of targets under the student distribution noise, detection uncertainty and clutters. One of possible future research on the RSEV-EGM-PHD is the efficiency and simplicity in implementation, which also the requirement of possible application such as to tracking in sensor networks. ACKNOWLEDGMENTS This research was supported by the Chinese Natural Science Foundation (granted No.61271317 and No.61175028). R EFERENCES [1] Y. Bar-Shalom and T. E. Fortmann, Tracking and Data Association. San Diego, CA: Academic, 1988. [2] L. D. Stone, C. A. Barlow, and T. L. Corwin, Bayesian Multiple Target Tracking. Boston: Artech House, 1999.
[3] Mahler, R., Random set theory for target tracking and identification. In D. Hall and J. Llinas (Eds)., Data Fusion Handbook,Boca Raton, FL: CRC Press, 2001. [4] Vo, B-N., S. Singh and A. Doucet, Sequential Monte Carlo implementation of the PHD filter for multitarget tracking, Proc. 6th Int. Conf. Information Fusion, Vol.2,July 2003, pp.792-799. [5] R. Mahler, Multi-target Bayes filtering via first-order multi-target moments, IEEE Trans. Aerosp. Electron. Syst., vol.39, no.4, pp.1152C1178, 2003. [6] Vo, B-N., and Ma, W. K. A closed-form solution for the probability hypothesis density filter. In Proceedings of the 7th International Conference on Information Fusion, vol.2, July 2005, 25 28. [7] Vo, B-N., Ma W K. The Gaussian mixture probability hypothesis density filter. IEEE Transactions on Signal Processing, 54(11):40914104,2006. [8] Vo, B-N., S. Singh, and A. Doucet,, Sequential Monte Carlo methods for multi-target filtering with random finite sets, IEEE Trans. Aerosp. Electron. Syst., vol. 41, no. 4, pp. 1224 C1245, 2005. [9] B. D. O. Anderson, J. B. Moore, Optimal Filtering, Prentice-Hall, 1979. [10] Ming Lei, Zhongliang Jing and Shiqiang Hu, Scaled Unscented Transform-based Variational Optimality Filter, 15th International Conference on Information Fusion. Singapore, 9-12th, July 2012. [11] Ming Lei, Christophe Baehr, Unscented/Ensemble Transform-based Variational Filter, Physica D: Nonlinear Phenomena, Vol.246, No.1, March 2013. [12] H. W. Sorenson and D. L. Alspach, Recursive Bayesian estimation using Gaussian sum, Automatica, vol.7, pp.465-479, 1971. [13] D. L. Alspach and H. W. Sorenson, Nonlinear Bayesian estimation using Gaussian sum approximations, IEEE Trans. Autom. Control, vol. AC17, no. 4, pp. 439C448, 1972. [14] S. J. Julier, The scaled unscented transformation, The Proceedings of the American Control Conference, Anchorage, AK, 2004. [15] Li, X.R.; Jilkov, V.P., Survey of maneuvering target tracking. Part V. Multiple-model methods,IEEE Trans. on Aerospace alectronic Systems, vol.41, no.4, pp.1255-1321, 2005. [16] X.Rong Li, Survey of Maneuvering Target Tracking. Part1:Dynamic Models, IEEE Trans. on Aerospace and Electronic Systems,Vol.39, no.4,pp: 1337-1364, 2003. [17] Schuhmacher, Dominic, B-T. Vo, and B-N. Vo., A consistent metric for performance evaluation of multi-object filters, IEEE Trans. on Signal Processing, vol.56, no.8, pp.3447-3457,2008.
