- Email: [email protected]
Particle filters for maneuvering target tracking problem
ARTICLE IN PRESS
Signal Processing 86 (2006) 195–203 www.elsevier.com/locate/sigpro
Particle filters for maneuvering target tracking problem Yihua Yua,, Qiansheng Chengb a
School of Information Engineering, Beijing University of Posts and Telecommunications, Beijing 100088, People’s Republic of China Department of Information Science, School of Mathematical Sciences, Peking University, Beijing 100871, People’s Republic of China
b
Received 8 August 2003; received in revised form 15 March 2005; accepted 15 March 2005 Available online 8 August 2005
Abstract In this paper, we address the target tracking problem for the case of maneuvering target, including single target and multiple target tracking. We propose a suitable model to characterize the maneuvering acceleration and develop a state space model to describe the maneuvering target tracking problem. Algorithms based on particle filters are developed. Simulations are given to demonstrate the performance of the procedures. r 2005 Elsevier B.V. All rights reserved. Keywords: Target tracking; Particle filter; Data association
1. Introduction Mobility tracking is an important issue for cellular radio network [1,2]. It enables efficient network control and provides the ability to offer additional services. Similar issues appear in many other civilian and military applications [3,4]. All these problems are related in that they can be described by state space models. Particle filters [5,6], which can be applied to estimate non-linear and non-Gaussian dynamic process, have received much attention in recent years. Many solutions for target tracking based on particle filters have been proposed [7,4]. Corresponding author.
E-mail addresses: [email protected] (Y. Yu), [email protected] (Q. Cheng).
Multitarget tracking deals with the state estimation of a number of moving targets. The extension of the particle filters to multitarget tracking has progressively received attention only in the last 10 years [3,8,9]. But these methods are mainly devoted to the random-acceleration target tracking and little work is related to the case of maneuvering target. Since the acceleration cannot be observed directly, maneuvering acceleration is usually difficult to track. In this paper, we focus on the maneuvering target tracking problem, including single target and multiple target tracking. We propose a suitable model to characterize the maneuvering acceleration and develop a state space model to describe the maneuvering target tracking problem. Particle filtering approach is applied and tracking algorithms are derived.
0165-1684/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2005.03.021
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
196
The rest of this paper is organized as follows. In Section 2, we describe our maneuvering acceleration model and motion models of target. In Sections 3 and 4, we present our algorithms for single target tracking and multiple target tracking, respectively. Simulation results are provided in Section 5 and a brief conclusion is given in Section 6.
2. Models Target tracking on a two-dimensional plane can be modelled as follows: let st ¼ ðst1 ; st2 ÞT be position vector, vt ¼ ðvt1 ; vt2 ÞT be velocity vector and at ¼ ðat1 ; at2 ÞT be acceleration vector, the target is supposed to evolve in the following way: 1 ! ! ! 0 Dt2 st st1 I 2 Dt I 2 I 2A at , ¼ þ@ 2 vt 0 I2 vt1 Dt I 2 (1) where Dt is the sample period, and I 2 is 2 2 identity matrix. If no maneuver exists, at is random acceleration and usually considered as the state noise. When some unknown maneuver exists, at is maneuvering acceleration. In this case, we characterize at by the following way: at ¼ at1 þ Dat ,
(2)
where Dat denotes the change of acceleration between time t 1 and t. We model Dat as Dat ¼ ut þ et ,
(3)
where t is a zero mean continuous random variable and ut is a discrete random variable. The choice of et and ut is dependent on different applications. For example, suppose the change of acceleration between two successive time points is not too large, we can choose et Nð0; 0:1I 2 Þ and ut 2 X fð0; 0ÞT ; ð1; 0ÞT ; ð1; 1ÞT ; ð0; 1ÞT ; ð1; 1ÞT , ð1; 0ÞT ; ð1; 1ÞT ; ð0; 1ÞT ; ð1; 1ÞT g. The distribution of ut on X is determined by some prior information.
From (1)–(3), we obtain the motion equation of maneuvering target as 0 1 2 0 1 0 1 Dt st B I 2 Dt I 2 2 I 2 C st1 CB B C B C CB v C B vt C ¼ B @ A B0 I2 Dt I 2 C@ t1 A @ A at at1 0 0 I2 1 1 0 2 0 2 Dt Dt B 2 I2 C B 2 I2 C C C B B C C B þ B Dt I Cut þ B ð4Þ B Dt I 2 Cet . 2 A A @ @ I2 I2 Basically, the observations that are related to the target can be described as yt ¼ hðst Þ þ Zt ,
(5)
where the observation noise Zt are characterized by its density and hðÞ is determined by different applications, such as the distance [4] relative to a reference point p ¼ ðp1 ; p2 ÞT , hðst Þ ¼ kst pk2 , and the direction of the target [7,4] relative to the reference point st2 p2 hðst Þ ¼ tan1 . st1 p1
3. Single target tracking 3.1. Description of the algorithm The state equation (4) and the observation equation (5) form a dynamic system for the maneuvering target tracking problem. Our goal is to estimate st ; vt ; at based on all available observations y1:t ¼ ðy1 ; . . . ; yt Þ up to time t. In this section, we develop an algorithm based on auxiliary particle filter [10] to this problem. From the Bayes rule, pðyt jst ; vt ; at Þpðst ; at ; vt jy1:t1 Þ pðyt jy1:t1 Þ pðyt jst Þpðst ; at ; vt jy1:t1 Þ . ¼ pðyt jy1:t1 Þ
pðst ; vt ; at jy1:t Þ ¼
ð6Þ
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
where the value set of k is f1; . . . ; MNg. Obviously, ^ t ; vt ; at jy1:t Þ is the marginal density of pðs pðst ; vt ; at ; kjy1:t Þ. In order to generate random samples from pðst ; vt ; at ; kjy1:t Þ, we choose the following importance sample density:
Moreover, pðst ; vt ; at jy1:t1 Þ Z ¼ pðst ; vt ; at ; st1 ; vt1 ; at1 , ut jy1:t1 Þ dst1 dvt1 dat1 dut Z ¼ pðst ; vt ; at jst1 ; vt1 ; at1 ; ut ; y1:t1 Þ
qðst ; vt ; at ; kjy1:t Þ / pðyt jmðkÞ t Þ ðkÞ ðkÞ ðkÞ ðkÞ pðst ; vt ; at jsðkÞ t1 ; vt1 ; at1 ; ut Þwt1 .
pðut jst1 ; vt1 ; at1 ; y1:t1 Þ pðst1 ; vt1 ; at1 jy1:t1 Þ dst1 dvt1 dat1 dut Z ¼ pðst ; vt ; at jst1 ; vt1 ; at1 ; ut Þpðut Þ
ð7Þ
where we utilize the Markovian structure of the state space model, and ut is independent on st1 ; vt1 ; at1 ; y1:t1 . Suppose a set of weighted random samples ðiÞ ðiÞ ðiÞ N fðsðiÞ from pðst1 ; vt1 ; t1 ; vt1 ; at1 Þ; wt1 gi¼1 ðjÞ M at1 jy1:t1 Þ and fut gj¼1 from pðut Þ are available, then (7) can be expressed as pðst ; vt ; at jy1:t1 Þ N X M 1 X ðiÞ ðiÞ ðjÞ ðiÞ pðst ; vt ; at jsðiÞ t1 ; vt1 ; at1 ; ut Þwt1 M i¼1 j¼1
9
MN 1 X ðkÞ ðkÞ ðkÞ ðkÞ pðst ; vt ; at jsðkÞ t1 ; vt1 ; at1 ; ut Þwt1 , M k¼1
^ t ; vt ; at jy1:t Þ / pðyt jst Þ pðs MN X
ðkÞ ðkÞ ðkÞ ðkÞ pðst ; vt ; at jsðkÞ t1 ; vt1 ; at1 ; ut Þwt1 .
From (11), Z qðkjy1:t Þ / pðyt jmðkÞ t Þ ðkÞ ðkÞ ðkÞ ðkÞ pðst ; vt ; at jsðkÞ t1 ; vt1 ; at1 ; ut Þwt1 dst dvt dat ðkÞ ¼ pðyt jmðkÞ t Þwt1 ,
then ðkÞ ðkÞ ðkÞ qðst ; vt ; at ; kjy1:t Þ / pðst ; vt ; at jsðkÞ t1 ; vt1 ; at1 ; ut Þ
qðkjy1:t Þ. The procedure to generate random samples from qðst ; vt ; at ; kjy1:t Þ can be implemented by the following way: (1) Generate N random samples fkðjÞ gN j¼1 from qðkjy1:t Þ. (2) For j ¼ 1; . . . ; N, generate a random sample ðjÞ ðjÞ ðsðjÞ according to the density t ; v t ; at Þ ðjÞ
ðjÞ
ðjÞ
ðjÞ
ðk Þ ðk Þ ðk Þ ðk Þ pðst ; vt ; at jst1 ; vt1 ; at1 ; ut Þ.
ð9Þ
k¼1
Next we develop an algorithm based on auxiliary particle filter to generate random samples of pðst ; vt ; at jy1:t Þ from (9). We add an auxiliary random variable k in (9) and get a new density pðst ; vt ; at ; kjy1:t Þ / pðyt jst Þ ðkÞ ðkÞ ðkÞ ðkÞ pðst ; vt ; at jsðkÞ t1 ; vt1 ; at1 ; ut Þwt1 ,
Dt2 ðkÞ Dt2 ðkÞ at1 þ u . 2 2 t
ð8Þ
where we reset the indices of the summation. From (6) and (8), an approximation of pðst ; vt ; at jy1:t Þ is given by
ð11Þ
ðkÞ ðkÞ ðkÞ ðkÞ Here, we let mðkÞ t ¼ Efst jst1 ; vt1 ; at1 ; ut g. From (4), we have ðkÞ ðkÞ mðkÞ t ¼ st1 þ Dt vt1 þ
pðst1 ; vt1 ; at1 jy1:t1 Þ dst1 dvt1 dat1 dut ,
197
ð10Þ
ðjÞ ðjÞ ðjÞ If each random sample ðsðjÞ t ; vt ; at ; k Þ is assigned a weight
wðjÞ t /
ðjÞ pðstðjÞ ; vtðjÞ ; aðjÞ t ; k jy1:t Þ ðjÞ qðstðjÞ ; vtðjÞ ; aðjÞ t ; k jy1:t Þ
¼
pðyt jsðjÞ t Þ ðjÞ
Þ pðyt jmðk Þ t
,
then we get the weighted random samples ðjÞ ðjÞ ðjÞ ðjÞ N fðsðjÞ t ; vt ; at ; k Þ; wt gj¼1 which are distributed ðjÞ ðjÞ ðjÞ N from pðst ; vt ; at ; kjy1:t Þ and fðsðjÞ t ; vt ; at Þ; wt gj¼1 which are distributed from pðst ; vt ; at jy1:t Þ.
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
198
Algorithm 1 (Single target tracking algorithm). ðiÞ ðiÞ ðiÞ N ðiÞ ðiÞ ðiÞ ½fðsðiÞ t ; vt ; at Þ; wt gi¼1 ¼ STT ½fðst1 ; vt1 ; at1 Þ; N wðiÞ t1 gi¼1 ðiÞ ðiÞ For each sample ðsðiÞ t1 ; vt1 ; at1 Þ, i ¼ 1; . . . ; N, generate M random samples futðjÞ gM j¼1 from pðut Þ, and denote ðkÞ ðkÞ ðkÞ ðkÞ NM fðsðkÞ t1 ; vt1 ; at1 ; ut Þ; wt1 gk¼1
mðkÞ t and fkðjÞ gN j¼1
For k ¼ 1; . . . ; NM, compute qðkjy1:t Þ. Generate N random samples from the set f1; . . . ; NMg with the probability proportional to qðkjy1:t Þ. For j ¼ 1; . . . ; N, generate one random sample ðjÞ
ðk Þ ðstðjÞ ; vtðjÞ ; aðjÞ t Þ from the density pðst ; vt ; at jst1 ; ðjÞ
ðjÞ
ðjÞ
wðjÞ t ¼
ðjÞ
ðk Þ Þ. pðyt jsðjÞ t Þ=pðyt jmt
3.2. Delayed-weight estimation The delayed-weight estimation method [11] is a simple extension of the particle filters. Since the state space model is highly correlated, the future observations often contain information about the current state. Hence, a delayed estimate is usually more accurate than the concurrent estimate. The basic idea of the delayed-weight estimation is as follows. Suppose at time t þ D, D40, the ðjÞ N N samples fxðjÞ 1:tþD gj¼1 and weights fwtþD gj¼1 , where x1:tþD ¼ ðx1 ; . . . ; xtþD Þ, are properly generated from pðx1:tþD jy1:tþD Þ, then we have the samples ðjÞ N fxtðjÞ gN j¼1 and weights fwtþD gj¼1 which are properly distributed from pðxt jy1:tþD Þ. With these samples and weights, we obtain a delayed estimate of the quantity of interest at time t as , N N X X ðjÞ ðjÞ ðjÞ Efhðxt Þjy1:tþD g hðxt ÞwtþD wtþD . j¼1
4.1. Multiple target tracking model Consider a system of M targets. Let xm t ¼ m m ðsm t ; vt ; at Þ, m ¼ 1; . . . ; M, be the state of the M targets at time t. The state of each target evolves in the following way: m m sm t ¼ st1 þ Dt vt1 þ
ðiÞ ðiÞ ðjÞ ðiÞ N M 9fðsðiÞ t1 ; vt1 ; at1 ; ut Þ; wt1 gi¼1 j¼1
Þ ðk Þ ðk Þ Þ and assign weight vðk t1 ; at1 ; ut
4. Multiple target tracking
j¼1
The delayed-weight estimation method incurs no additional computational cost, but it requires some extra memory for storing fðxðjÞ tþ1 ; . . . ; ðjÞ xtþD ÞgN . j¼1
Dt2 m at , 2
(12)
m m vm t ¼ vt1 þ Dt at ,
(13)
m m am t ¼ at1 þ Dat .
(14)
We model Dam t in the same way as in Section 2, i.e., m m Dam Let xt ¼ ðx1t ; . . . ; xM and t ¼ ut þ e t . t Þ 1 M ut ¼ ðut ; . . . ; ut Þ. The observation vector collected at time t is t denoted by yt ¼ ðy1t ; . . . ; yN t Þ, where N t is the number of observations collected at time t. Because of the effect of observation noise and false alarm, the origin of each observation is not known. It may come from one target or from false alarm. The false alarms are assumed to be uniformly distributed in the observation area. Their number is assumed to arise from Poisson density of parameter lV , where V is the volume of the observation area, and l is the number of false alarms per unit volume. The most common traditional approaches for multitarget tracking consist of two sequential steps: (1) data association, and (2) target state estimation. 4.2. Data association The most commonly used method for data association is probably the joint probability data association filter (JPDAF) [12,13]. In this paper, we take advantage of JPDAF for data association. As we do not know the origin of each observation, one has to introduce the vector At to describe the associations between the observations and the targets. Each component Ant , n ¼ 1; . . . ; N t , is a random variable that takes its values among f0; . . . ; Mg. Ant ¼ 0 indicates that ynt is associated with the false alarm and Ant ¼ mðma0Þ
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
indicates that ynt is associated with the mth target. In the case of Ant ¼ mðma0Þ, ynt is a realization of the following observation equation: ynt
¼
hðxm t Þ
þ
Zm t .
n bnm t ¼ pðAt ¼ mÞ;
and an approximation of pðxm t jy1:t Þ is given by jy Þ pbðxm 1:t t /
The JPDAF begins with the gating of the observations. Only the observations that are inside an ellipsoid around the predicted state of one target are kept. Then the probabilities of each data association n ¼ 1; . . . ; N t ; m ¼ 0; . . . ; M,
are estimated. For more details of JPDAF, we refer to [13,14].
199
9
N X R 1 X m m;ðiÞ m;ðjÞ pðyt jxm Þwm;ðiÞ t Þpðxt jxt1 ; ut t1 R i¼1 j¼1 NR X 1 m;ðkÞ m;ðkÞ pðyt jxm Þ pðxm Þwm;ðkÞ t t jxt1 ; ut t1 , R k¼1
ð16Þ
where we reset the indices of the summation. The following procedure can be used to generate the weighted random samples fxm;ðiÞ ; wm;ðiÞ gN t t i¼1 from (16) which are distributed from pðxm t jy1:t Þ, m ¼ 1; . . . ; M. Algorithm 2 (Multitarget tracking algorithm). m;ðiÞ m;ðiÞ N M M ; wm;ðiÞ gN ½fxm;ðiÞ t t i¼1 m¼1 ¼ MTT½fxt1 ; wt1 gi¼1 m¼1 :
4.3. State estimation With the output of data association, we estimate the state of each target. Our goal is to track M posterior density pðxm t jy1:t Þ, m ¼ 1; . . . ; M. From the Bayes rule, m pðyt jxm t Þpðxt jy1:t1 Þ pðxm , t jy1:t Þ ¼ pðyt jy1:t1 Þ
and pðxm t jy1:t1 Þ Z ¼ pðxm t ; xt1 ; ut jy1:t1 Þ dxt1 dut Z ¼ pðxm t jxt1 ; ut Þpðxt1 ; ut jy1:t1 Þ dxt1 dut Z m m m m m m ¼ pðxm t jxt1 ; ut Þpðxt1 ; ut jy1:t1 Þ dxt1 dut Z m m m m m m ¼ pðxm t jxt1 ; ut Þpðxt1 jy1:t1 Þpðut Þ dxt1 dut . Here, we utilize the fact the xm t is only related to m xm t1 and ut and independent on other elements in xt1 and ut from Eqs. (12)–(14). Suppose a set of random samples and weights m;ðiÞ N m fxm;ðiÞ t1 ; wt1 gi¼1 from pðxt1 jy1:t1 Þ is available at m;ðjÞ R time t 1, and fut gj¼1 are generated from pðum t Þ, then N X R 1 X m;ðiÞ m;ðjÞ pðxm jy Þ pðxm Þwm;ðiÞ 1:t1 t t jxt1 ; ut t1 , R i¼1 j¼1
(15)
For m ¼ 1; . . . ; M, i ¼ 1; . . . ; N, generate R m random samples fum;ðjÞ gR t j¼1 from pðut Þ and denote m;ðkÞ NR Þ; wm;ðkÞ fðxm;ðkÞ t1 ; ut t1 gk¼1 m;ðjÞ N R 9fðxm;ðiÞ Þ; wm;ðiÞ t1 ; ut t1 gi¼1 j¼1 .
For m ¼ 1; . . . ; M; k ¼ 1; . . . ; NR, generate one random sample xm;ðkÞ from the density t m;ðkÞ m;ðkÞ pðxm jx ; u Þ. t t t1 For each xm;ðkÞ , assign a weight t m;ðkÞ ¼ wm;ðkÞ Þ, wm;ðkÞ t t1 pðyt jxt m ¼ 1; . . . ; M; k ¼ 1; . . . ; NR.
ð17Þ
For m ¼ 1; . . . ; M, generate N random samples from the set fxm;ðkÞ gNR t k¼1 according to the density m;ðkÞ m;ðkÞ pðxt Þ / wt , and denote the results as m;ðiÞ fxm;ðiÞ ; wm;ðiÞ gN is reset t t i¼1 , where each weight wt as 1=N. The computation of the weight (17) in Algorithm 2 needs the output of the data association. From the total probability theorem with the event ðAnt ¼ mÞ, n ¼ 1; . . . ; N t , pðyt jxm t Þ¼
Nt X
n n pðyt jxm t ; At ¼ mÞpðAt ¼ mÞ
n¼1
¼
Nt X n¼1
nm n pðynt jxm t ; At ¼ mÞbt .
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
200
)1=2
5. Experiments 2
st2 ðgÞÞ þðst2 ðgÞ b
5.1. Single target tracking With the sample period Dt ¼ 1s in (1), we obtain the motion model of the target. The observation includes two elements: the relative distance and the relative direction from the target to a reference point p ¼ ðp1 ; p2 ÞT , 0
1 kst pk2 B st2 p2 C yt ¼ @ A þ Zt . tan1 st1 p1
(18)
In this example, we choose p ¼ ð1; 1ÞT and 2
0 Zt Nðð00Þ; ðð2:0Þ 0 ð0:05Þ2 ÞÞ. Before the simulation, we artificially give the maneuvering accelerations, the actual initial position s0 ¼ ð0; 0ÞT and the actual initial velocity v0 ¼ ð0; 0ÞT of the target. With these settings, we get the simulated position st and velocity vt from (1) and the simulated observation data yt from (18), where t ¼ 1; . . . ; L and L ¼ 100 is taken. With the simulated data y1:L ¼ fy1 ; . . . ; yL g, we apply Algorithm 1 to estimate the simulated position, velocity and acceleration of the target, with N ¼ 500. The distribution of pðut Þ is set to 1 pðut Þ ¼ 45 if ut ¼ ð0; 0ÞT , and pðut Þ ¼ 58 otherwise. In Algorithm 1, M ¼ 5 samples are generated from pðut Þ every time. The results of the estimation are shown in Fig. 1 and the results of the delayed-weight estimation with D ¼ 10 are shown in Fig. 2. Since the position of target can be directly computed from the noisy observations. Fig. 1(f) presents the results of direct computation. But the velocity and acceleration cannot be directly computed. We perform G ¼ 100 runs of the simulations. All the computations are performed on Pentium 4 1.7 GHz PC with Matlab 6.0. The performance measure is root mean squared error (RMS) computed as follows [15]: ( L G 1 X 1 X 2 ðRMSÞ ¼ ½ðst1 ðgÞ b st1 ðgÞÞ2 L t¼1 G g¼1
,
where ðst1 ðgÞ; st2 ðgÞÞT is the simulated target position and ðb st1 ðgÞ; b st2 ðgÞÞT is the estimated target position at time t of the gth simulation. The results of performance are RMS ¼ 2:3977 for Algorithm 1, RMS ¼ 1:8737 for the delayed-weight estimation method, and RMS ¼ 2:7536 for the direct computation, respectively. 5.2. Multitarget tracking In this simulation study, the proposed MTT algorithm is implemented with five targets and Dt ¼ 1 s. An observation produced by the mth target is generated according to 0 1 ksm t pk2 m B s p2 C yt ¼ @ A þ Zt , tan1 t2 sm p 1 t1 2 0 where p ¼ ð1; 1ÞT and Zt N 00 ; ð1:0Þ . 0 ð0:04Þ2 The number of false alarm follows the Poisson distribution with mean lV ¼ 2. These false alarms are assumed to be independent and uniformly distributed within the observation volume V. With these settings, the states of the five targets and the observations are simulated with t ¼ 1; . . . ; L and L ¼ 100. The trajectories of the targets are represented in Fig. 3, where each ‘‘ ’’ indicates the starting position and each ‘‘ ’’ indicates the ending position of one target. From the simulated observations, the proposed MTT algorithm is implemented to estimate the states of the targets. The results are presented in Fig. 3(a). The delayed-weight estimation method is also applied here and the results are presented in Fig. 3(b). If two targets are in the approximately same position in the two-dimensional plane at the same time, it is usually difficult to distinguish them. Since our algorithm estimates the velocity of target, we can take advantage of the velocity information to distinguish them. Moreover, if the velocities of these two targets are also approximate, we can take advantage of the acceleration
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
201
10 200
5 0
150
–5 100 –10 50
–15
0
–20 –100
–50
0
50
100
150
(a)
0
10
20
30
40
50
60
70
80
90 100
0
10
20
30
40
50
60
70
80
90 100
(b) 15
2 1.5
10 1 0.5
5
0 0
–0.5 –1
–5
–1.5 –10
–2 0
10
20
30
40
50
60
70
80
90 100
(c)
(d) 300
2.5 2
250
1.5 200
1 0.5
150
0 100
–0.5 –1
50
–1.5 0
–2 0
(e)
10
20
30
40
50
60
70
80
90 100
–200 –150 –100
–50
0
50
100
150
(f)
Fig. 1. The simulated data (solid line) and estimated data by Algorithm 1 (dashed line): (a) target position; (b) x-coordinate of target velocity; (c) y-coordinate of target velocity; (d) x-coordinate of target acceleration; (e) y-coordinate of target acceleration; (f) the result of direct computation.
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
202
10 200
5
150
0
100
–5
50
–10
0 –100
–50
0
50
100
–15
150
0
10
20
30
40
50
60
70
80
90
100
0
10
20
30
40
50
60
70
80
90 100
(b)
(a) 12
2
10
1.5
8 1
6 4
0.5
2
0
0
–0.5
–2
–1
–4 –1.5
–6
–2
–8 0
10
20
30
40
50
60
70
80
90
100 (d)
(c) 2 1.5 1 0.5 0 –0.5 –1 –1.5 0
10
20
30
40
50
60
70
80
90 100
(e) Fig. 2. The simulated data (solid line) and estimated data by delayed-weight estimation method (dashed line): (a) target position; (b) xcoordinate of target velocity; (c) y-coordinate of target velocity; (d) x-coordinate of target acceleration; (e) y-coordinate of target acceleration.
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203 800
800
700
700 4
600
500
400
400 5
300 200
3
200
0
3 2
0
1 –300 –200–100
5
300
100
2
0 –100
4
600
500
100
203
100 200 300 400 500 600 700
(a)
–100
1 –300 –200–100
0
100 200 300 400 500 600 700
(b)
Fig. 3. The simulated positions (solid line) and the estimated positions (dashed line) of the multitarget tracking with five targets by (a) Algorithm 2 and (b) delayed-weight estimation method.
information. But if the positions, velocities and accelerations of these two targets are all approximate, we cannot distinguish them any more. 6. Conclusions In this paper, algorithms are presented for the target tracking problem for the case of the maneuvering target. A suitable model is proposed to characterize the maneuvering acceleration. Algorithms based on particle filters are developed. Finally, two simulations are given to demonstrate the performance of the procedures. References [1] M. Hellebrandt, R. Mathar, Location tracking of mobiles in cellular radio networks, IEEE Trans. Veh. Technol. 48 (5) (September 1999) 1158–1162. [2] T. Liu, P. Bahl, I. Chlamtac, Mobility modeling, location tracking, and trajectory prediction in wireless ATM networks, IEEE J. Select. Areas Commun. 16 (6) (August 1998) 922–936. [3] D. Avitzour, Stochastic simulation Bayesian approach to multitarget tracking, IEE Proc.-F 142 (2) (April 1995) 41–44. [4] F. Gustafsson, F. Gunnarsson, N. Bergman, U. Forssell, J. Jansson, R. Karlsson, P. Nordlund, Particle filters for positioning, navigation and tracking, IEEE Trans. Signal Process. 50 (2) (February 2002) 425–437.
[5] A. Doucet, S.J. Godsill, C. Andrieu, On sequential Monte Carlo sampling methods for Bayesian filtering, Statist. Comput. 10 (3) (July 2000) 197–208. [6] J.S. Liu, R. Chen, Sequential Monte Carlo methods for dynamic systems, J. Amer. Statist. Assoc. 93 (443) (September 1998) 1032–1044. [7] N.J. Gordon, D.J. Salmond, A.F.M. Smith, Novel approach to nonlinear/non-Gaussian Bayesian state estimation, IEE Proc.-F 140 (2) (April 1993) 107–113. [8] N. Gordon, A hybrid bootstrap filter for target tracking in clutter, IEEE Trans. Aerosp. Electron. Syst. 33 (1) (January 1997) 353–358. [9] C. Hue, J.L. Cadre, P. Perez, Sequential Monte Carlo methods for multiple target tracking and data fusion, IEEE Trans. Signal Process. 50 (2) (February 2002) 309–325. [10] M.K. Pitt, N. Shephard, Filtering via simulation: auxiliary particle filters, J. Amer. Statist. Assoc. 94 (446) (June 1999) 590–599. [11] R. Chen, X.D. Wang, J.S. Liu, Adaptive joint detection and decoding in flat-fading channels via mixture Kalman filtering, IEEE Trans. Inform. Theory 46 (6) (September 2000) 2079–2094. [12] Y. Bar-Shalom, T.E. Fortmann, Tracking and Data Association, Academic Press, New York, 1988. [13] T.E. Fortmann, Y. Bar-Shalom, M. Scheffe, Sonar tracking of multiple targets using joint probabilistic data association, IEEE J. Oceanic Eng. OE-8 (July 1983) 173–184. [14] B. Zhou, N.K. Bose, Multitarget tracking in clutter: fast algorithm for data association, IEEE Trans. Aerosp. Electron. Syst. 29 (2) (April 1993) 352–363. [15] A. Doucet, N.J. Gordon, V. Krishnamurthy, Particle filters for state estimation of jump Markov linear systems, IEEE Trans. Signal Process. 49 (3) (March 2001) 613–624.
Signal Processing 86 (2006) 195–203 www.elsevier.com/locate/sigpro
Particle filters for maneuvering target tracking problem Yihua Yua,, Qiansheng Chengb a
School of Information Engineering, Beijing University of Posts and Telecommunications, Beijing 100088, People’s Republic of China Department of Information Science, School of Mathematical Sciences, Peking University, Beijing 100871, People’s Republic of China
b
Received 8 August 2003; received in revised form 15 March 2005; accepted 15 March 2005 Available online 8 August 2005
Abstract In this paper, we address the target tracking problem for the case of maneuvering target, including single target and multiple target tracking. We propose a suitable model to characterize the maneuvering acceleration and develop a state space model to describe the maneuvering target tracking problem. Algorithms based on particle filters are developed. Simulations are given to demonstrate the performance of the procedures. r 2005 Elsevier B.V. All rights reserved. Keywords: Target tracking; Particle filter; Data association
1. Introduction Mobility tracking is an important issue for cellular radio network [1,2]. It enables efficient network control and provides the ability to offer additional services. Similar issues appear in many other civilian and military applications [3,4]. All these problems are related in that they can be described by state space models. Particle filters [5,6], which can be applied to estimate non-linear and non-Gaussian dynamic process, have received much attention in recent years. Many solutions for target tracking based on particle filters have been proposed [7,4]. Corresponding author.
E-mail addresses: [email protected] (Y. Yu), [email protected] (Q. Cheng).
Multitarget tracking deals with the state estimation of a number of moving targets. The extension of the particle filters to multitarget tracking has progressively received attention only in the last 10 years [3,8,9]. But these methods are mainly devoted to the random-acceleration target tracking and little work is related to the case of maneuvering target. Since the acceleration cannot be observed directly, maneuvering acceleration is usually difficult to track. In this paper, we focus on the maneuvering target tracking problem, including single target and multiple target tracking. We propose a suitable model to characterize the maneuvering acceleration and develop a state space model to describe the maneuvering target tracking problem. Particle filtering approach is applied and tracking algorithms are derived.
0165-1684/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2005.03.021
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
196
The rest of this paper is organized as follows. In Section 2, we describe our maneuvering acceleration model and motion models of target. In Sections 3 and 4, we present our algorithms for single target tracking and multiple target tracking, respectively. Simulation results are provided in Section 5 and a brief conclusion is given in Section 6.
2. Models Target tracking on a two-dimensional plane can be modelled as follows: let st ¼ ðst1 ; st2 ÞT be position vector, vt ¼ ðvt1 ; vt2 ÞT be velocity vector and at ¼ ðat1 ; at2 ÞT be acceleration vector, the target is supposed to evolve in the following way: 1 ! ! ! 0 Dt2 st st1 I 2 Dt I 2 I 2A at , ¼ þ@ 2 vt 0 I2 vt1 Dt I 2 (1) where Dt is the sample period, and I 2 is 2 2 identity matrix. If no maneuver exists, at is random acceleration and usually considered as the state noise. When some unknown maneuver exists, at is maneuvering acceleration. In this case, we characterize at by the following way: at ¼ at1 þ Dat ,
(2)
where Dat denotes the change of acceleration between time t 1 and t. We model Dat as Dat ¼ ut þ et ,
(3)
where t is a zero mean continuous random variable and ut is a discrete random variable. The choice of et and ut is dependent on different applications. For example, suppose the change of acceleration between two successive time points is not too large, we can choose et Nð0; 0:1I 2 Þ and ut 2 X fð0; 0ÞT ; ð1; 0ÞT ; ð1; 1ÞT ; ð0; 1ÞT ; ð1; 1ÞT , ð1; 0ÞT ; ð1; 1ÞT ; ð0; 1ÞT ; ð1; 1ÞT g. The distribution of ut on X is determined by some prior information.
From (1)–(3), we obtain the motion equation of maneuvering target as 0 1 2 0 1 0 1 Dt st B I 2 Dt I 2 2 I 2 C st1 CB B C B C CB v C B vt C ¼ B @ A B0 I2 Dt I 2 C@ t1 A @ A at at1 0 0 I2 1 1 0 2 0 2 Dt Dt B 2 I2 C B 2 I2 C C C B B C C B þ B Dt I Cut þ B ð4Þ B Dt I 2 Cet . 2 A A @ @ I2 I2 Basically, the observations that are related to the target can be described as yt ¼ hðst Þ þ Zt ,
(5)
where the observation noise Zt are characterized by its density and hðÞ is determined by different applications, such as the distance [4] relative to a reference point p ¼ ðp1 ; p2 ÞT , hðst Þ ¼ kst pk2 , and the direction of the target [7,4] relative to the reference point st2 p2 hðst Þ ¼ tan1 . st1 p1
3. Single target tracking 3.1. Description of the algorithm The state equation (4) and the observation equation (5) form a dynamic system for the maneuvering target tracking problem. Our goal is to estimate st ; vt ; at based on all available observations y1:t ¼ ðy1 ; . . . ; yt Þ up to time t. In this section, we develop an algorithm based on auxiliary particle filter [10] to this problem. From the Bayes rule, pðyt jst ; vt ; at Þpðst ; at ; vt jy1:t1 Þ pðyt jy1:t1 Þ pðyt jst Þpðst ; at ; vt jy1:t1 Þ . ¼ pðyt jy1:t1 Þ
pðst ; vt ; at jy1:t Þ ¼
ð6Þ
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
where the value set of k is f1; . . . ; MNg. Obviously, ^ t ; vt ; at jy1:t Þ is the marginal density of pðs pðst ; vt ; at ; kjy1:t Þ. In order to generate random samples from pðst ; vt ; at ; kjy1:t Þ, we choose the following importance sample density:
Moreover, pðst ; vt ; at jy1:t1 Þ Z ¼ pðst ; vt ; at ; st1 ; vt1 ; at1 , ut jy1:t1 Þ dst1 dvt1 dat1 dut Z ¼ pðst ; vt ; at jst1 ; vt1 ; at1 ; ut ; y1:t1 Þ
qðst ; vt ; at ; kjy1:t Þ / pðyt jmðkÞ t Þ ðkÞ ðkÞ ðkÞ ðkÞ pðst ; vt ; at jsðkÞ t1 ; vt1 ; at1 ; ut Þwt1 .
pðut jst1 ; vt1 ; at1 ; y1:t1 Þ pðst1 ; vt1 ; at1 jy1:t1 Þ dst1 dvt1 dat1 dut Z ¼ pðst ; vt ; at jst1 ; vt1 ; at1 ; ut Þpðut Þ
ð7Þ
where we utilize the Markovian structure of the state space model, and ut is independent on st1 ; vt1 ; at1 ; y1:t1 . Suppose a set of weighted random samples ðiÞ ðiÞ ðiÞ N fðsðiÞ from pðst1 ; vt1 ; t1 ; vt1 ; at1 Þ; wt1 gi¼1 ðjÞ M at1 jy1:t1 Þ and fut gj¼1 from pðut Þ are available, then (7) can be expressed as pðst ; vt ; at jy1:t1 Þ N X M 1 X ðiÞ ðiÞ ðjÞ ðiÞ pðst ; vt ; at jsðiÞ t1 ; vt1 ; at1 ; ut Þwt1 M i¼1 j¼1
9
MN 1 X ðkÞ ðkÞ ðkÞ ðkÞ pðst ; vt ; at jsðkÞ t1 ; vt1 ; at1 ; ut Þwt1 , M k¼1
^ t ; vt ; at jy1:t Þ / pðyt jst Þ pðs MN X
ðkÞ ðkÞ ðkÞ ðkÞ pðst ; vt ; at jsðkÞ t1 ; vt1 ; at1 ; ut Þwt1 .
From (11), Z qðkjy1:t Þ / pðyt jmðkÞ t Þ ðkÞ ðkÞ ðkÞ ðkÞ pðst ; vt ; at jsðkÞ t1 ; vt1 ; at1 ; ut Þwt1 dst dvt dat ðkÞ ¼ pðyt jmðkÞ t Þwt1 ,
then ðkÞ ðkÞ ðkÞ qðst ; vt ; at ; kjy1:t Þ / pðst ; vt ; at jsðkÞ t1 ; vt1 ; at1 ; ut Þ
qðkjy1:t Þ. The procedure to generate random samples from qðst ; vt ; at ; kjy1:t Þ can be implemented by the following way: (1) Generate N random samples fkðjÞ gN j¼1 from qðkjy1:t Þ. (2) For j ¼ 1; . . . ; N, generate a random sample ðjÞ ðjÞ ðsðjÞ according to the density t ; v t ; at Þ ðjÞ
ðjÞ
ðjÞ
ðjÞ
ðk Þ ðk Þ ðk Þ ðk Þ pðst ; vt ; at jst1 ; vt1 ; at1 ; ut Þ.
ð9Þ
k¼1
Next we develop an algorithm based on auxiliary particle filter to generate random samples of pðst ; vt ; at jy1:t Þ from (9). We add an auxiliary random variable k in (9) and get a new density pðst ; vt ; at ; kjy1:t Þ / pðyt jst Þ ðkÞ ðkÞ ðkÞ ðkÞ pðst ; vt ; at jsðkÞ t1 ; vt1 ; at1 ; ut Þwt1 ,
Dt2 ðkÞ Dt2 ðkÞ at1 þ u . 2 2 t
ð8Þ
where we reset the indices of the summation. From (6) and (8), an approximation of pðst ; vt ; at jy1:t Þ is given by
ð11Þ
ðkÞ ðkÞ ðkÞ ðkÞ Here, we let mðkÞ t ¼ Efst jst1 ; vt1 ; at1 ; ut g. From (4), we have ðkÞ ðkÞ mðkÞ t ¼ st1 þ Dt vt1 þ
pðst1 ; vt1 ; at1 jy1:t1 Þ dst1 dvt1 dat1 dut ,
197
ð10Þ
ðjÞ ðjÞ ðjÞ If each random sample ðsðjÞ t ; vt ; at ; k Þ is assigned a weight
wðjÞ t /
ðjÞ pðstðjÞ ; vtðjÞ ; aðjÞ t ; k jy1:t Þ ðjÞ qðstðjÞ ; vtðjÞ ; aðjÞ t ; k jy1:t Þ
¼
pðyt jsðjÞ t Þ ðjÞ
Þ pðyt jmðk Þ t
,
then we get the weighted random samples ðjÞ ðjÞ ðjÞ ðjÞ N fðsðjÞ t ; vt ; at ; k Þ; wt gj¼1 which are distributed ðjÞ ðjÞ ðjÞ N from pðst ; vt ; at ; kjy1:t Þ and fðsðjÞ t ; vt ; at Þ; wt gj¼1 which are distributed from pðst ; vt ; at jy1:t Þ.
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
198
Algorithm 1 (Single target tracking algorithm). ðiÞ ðiÞ ðiÞ N ðiÞ ðiÞ ðiÞ ½fðsðiÞ t ; vt ; at Þ; wt gi¼1 ¼ STT ½fðst1 ; vt1 ; at1 Þ; N wðiÞ t1 gi¼1 ðiÞ ðiÞ For each sample ðsðiÞ t1 ; vt1 ; at1 Þ, i ¼ 1; . . . ; N, generate M random samples futðjÞ gM j¼1 from pðut Þ, and denote ðkÞ ðkÞ ðkÞ ðkÞ NM fðsðkÞ t1 ; vt1 ; at1 ; ut Þ; wt1 gk¼1
mðkÞ t and fkðjÞ gN j¼1
For k ¼ 1; . . . ; NM, compute qðkjy1:t Þ. Generate N random samples from the set f1; . . . ; NMg with the probability proportional to qðkjy1:t Þ. For j ¼ 1; . . . ; N, generate one random sample ðjÞ
ðk Þ ðstðjÞ ; vtðjÞ ; aðjÞ t Þ from the density pðst ; vt ; at jst1 ; ðjÞ
ðjÞ
ðjÞ
wðjÞ t ¼
ðjÞ
ðk Þ Þ. pðyt jsðjÞ t Þ=pðyt jmt
3.2. Delayed-weight estimation The delayed-weight estimation method [11] is a simple extension of the particle filters. Since the state space model is highly correlated, the future observations often contain information about the current state. Hence, a delayed estimate is usually more accurate than the concurrent estimate. The basic idea of the delayed-weight estimation is as follows. Suppose at time t þ D, D40, the ðjÞ N N samples fxðjÞ 1:tþD gj¼1 and weights fwtþD gj¼1 , where x1:tþD ¼ ðx1 ; . . . ; xtþD Þ, are properly generated from pðx1:tþD jy1:tþD Þ, then we have the samples ðjÞ N fxtðjÞ gN j¼1 and weights fwtþD gj¼1 which are properly distributed from pðxt jy1:tþD Þ. With these samples and weights, we obtain a delayed estimate of the quantity of interest at time t as , N N X X ðjÞ ðjÞ ðjÞ Efhðxt Þjy1:tþD g hðxt ÞwtþD wtþD . j¼1
4.1. Multiple target tracking model Consider a system of M targets. Let xm t ¼ m m ðsm t ; vt ; at Þ, m ¼ 1; . . . ; M, be the state of the M targets at time t. The state of each target evolves in the following way: m m sm t ¼ st1 þ Dt vt1 þ
ðiÞ ðiÞ ðjÞ ðiÞ N M 9fðsðiÞ t1 ; vt1 ; at1 ; ut Þ; wt1 gi¼1 j¼1
Þ ðk Þ ðk Þ Þ and assign weight vðk t1 ; at1 ; ut
4. Multiple target tracking
j¼1
The delayed-weight estimation method incurs no additional computational cost, but it requires some extra memory for storing fðxðjÞ tþ1 ; . . . ; ðjÞ xtþD ÞgN . j¼1
Dt2 m at , 2
(12)
m m vm t ¼ vt1 þ Dt at ,
(13)
m m am t ¼ at1 þ Dat .
(14)
We model Dam t in the same way as in Section 2, i.e., m m Dam Let xt ¼ ðx1t ; . . . ; xM and t ¼ ut þ e t . t Þ 1 M ut ¼ ðut ; . . . ; ut Þ. The observation vector collected at time t is t denoted by yt ¼ ðy1t ; . . . ; yN t Þ, where N t is the number of observations collected at time t. Because of the effect of observation noise and false alarm, the origin of each observation is not known. It may come from one target or from false alarm. The false alarms are assumed to be uniformly distributed in the observation area. Their number is assumed to arise from Poisson density of parameter lV , where V is the volume of the observation area, and l is the number of false alarms per unit volume. The most common traditional approaches for multitarget tracking consist of two sequential steps: (1) data association, and (2) target state estimation. 4.2. Data association The most commonly used method for data association is probably the joint probability data association filter (JPDAF) [12,13]. In this paper, we take advantage of JPDAF for data association. As we do not know the origin of each observation, one has to introduce the vector At to describe the associations between the observations and the targets. Each component Ant , n ¼ 1; . . . ; N t , is a random variable that takes its values among f0; . . . ; Mg. Ant ¼ 0 indicates that ynt is associated with the false alarm and Ant ¼ mðma0Þ
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
indicates that ynt is associated with the mth target. In the case of Ant ¼ mðma0Þ, ynt is a realization of the following observation equation: ynt
¼
hðxm t Þ
þ
Zm t .
n bnm t ¼ pðAt ¼ mÞ;
and an approximation of pðxm t jy1:t Þ is given by jy Þ pbðxm 1:t t /
The JPDAF begins with the gating of the observations. Only the observations that are inside an ellipsoid around the predicted state of one target are kept. Then the probabilities of each data association n ¼ 1; . . . ; N t ; m ¼ 0; . . . ; M,
are estimated. For more details of JPDAF, we refer to [13,14].
199
9
N X R 1 X m m;ðiÞ m;ðjÞ pðyt jxm Þwm;ðiÞ t Þpðxt jxt1 ; ut t1 R i¼1 j¼1 NR X 1 m;ðkÞ m;ðkÞ pðyt jxm Þ pðxm Þwm;ðkÞ t t jxt1 ; ut t1 , R k¼1
ð16Þ
where we reset the indices of the summation. The following procedure can be used to generate the weighted random samples fxm;ðiÞ ; wm;ðiÞ gN t t i¼1 from (16) which are distributed from pðxm t jy1:t Þ, m ¼ 1; . . . ; M. Algorithm 2 (Multitarget tracking algorithm). m;ðiÞ m;ðiÞ N M M ; wm;ðiÞ gN ½fxm;ðiÞ t t i¼1 m¼1 ¼ MTT½fxt1 ; wt1 gi¼1 m¼1 :
4.3. State estimation With the output of data association, we estimate the state of each target. Our goal is to track M posterior density pðxm t jy1:t Þ, m ¼ 1; . . . ; M. From the Bayes rule, m pðyt jxm t Þpðxt jy1:t1 Þ pðxm , t jy1:t Þ ¼ pðyt jy1:t1 Þ
and pðxm t jy1:t1 Þ Z ¼ pðxm t ; xt1 ; ut jy1:t1 Þ dxt1 dut Z ¼ pðxm t jxt1 ; ut Þpðxt1 ; ut jy1:t1 Þ dxt1 dut Z m m m m m m ¼ pðxm t jxt1 ; ut Þpðxt1 ; ut jy1:t1 Þ dxt1 dut Z m m m m m m ¼ pðxm t jxt1 ; ut Þpðxt1 jy1:t1 Þpðut Þ dxt1 dut . Here, we utilize the fact the xm t is only related to m xm t1 and ut and independent on other elements in xt1 and ut from Eqs. (12)–(14). Suppose a set of random samples and weights m;ðiÞ N m fxm;ðiÞ t1 ; wt1 gi¼1 from pðxt1 jy1:t1 Þ is available at m;ðjÞ R time t 1, and fut gj¼1 are generated from pðum t Þ, then N X R 1 X m;ðiÞ m;ðjÞ pðxm jy Þ pðxm Þwm;ðiÞ 1:t1 t t jxt1 ; ut t1 , R i¼1 j¼1
(15)
For m ¼ 1; . . . ; M, i ¼ 1; . . . ; N, generate R m random samples fum;ðjÞ gR t j¼1 from pðut Þ and denote m;ðkÞ NR Þ; wm;ðkÞ fðxm;ðkÞ t1 ; ut t1 gk¼1 m;ðjÞ N R 9fðxm;ðiÞ Þ; wm;ðiÞ t1 ; ut t1 gi¼1 j¼1 .
For m ¼ 1; . . . ; M; k ¼ 1; . . . ; NR, generate one random sample xm;ðkÞ from the density t m;ðkÞ m;ðkÞ pðxm jx ; u Þ. t t t1 For each xm;ðkÞ , assign a weight t m;ðkÞ ¼ wm;ðkÞ Þ, wm;ðkÞ t t1 pðyt jxt m ¼ 1; . . . ; M; k ¼ 1; . . . ; NR.
ð17Þ
For m ¼ 1; . . . ; M, generate N random samples from the set fxm;ðkÞ gNR t k¼1 according to the density m;ðkÞ m;ðkÞ pðxt Þ / wt , and denote the results as m;ðiÞ fxm;ðiÞ ; wm;ðiÞ gN is reset t t i¼1 , where each weight wt as 1=N. The computation of the weight (17) in Algorithm 2 needs the output of the data association. From the total probability theorem with the event ðAnt ¼ mÞ, n ¼ 1; . . . ; N t , pðyt jxm t Þ¼
Nt X
n n pðyt jxm t ; At ¼ mÞpðAt ¼ mÞ
n¼1
¼
Nt X n¼1
nm n pðynt jxm t ; At ¼ mÞbt .
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
200
)1=2
5. Experiments 2
st2 ðgÞÞ þðst2 ðgÞ b
5.1. Single target tracking With the sample period Dt ¼ 1s in (1), we obtain the motion model of the target. The observation includes two elements: the relative distance and the relative direction from the target to a reference point p ¼ ðp1 ; p2 ÞT , 0
1 kst pk2 B st2 p2 C yt ¼ @ A þ Zt . tan1 st1 p1
(18)
In this example, we choose p ¼ ð1; 1ÞT and 2
0 Zt Nðð00Þ; ðð2:0Þ 0 ð0:05Þ2 ÞÞ. Before the simulation, we artificially give the maneuvering accelerations, the actual initial position s0 ¼ ð0; 0ÞT and the actual initial velocity v0 ¼ ð0; 0ÞT of the target. With these settings, we get the simulated position st and velocity vt from (1) and the simulated observation data yt from (18), where t ¼ 1; . . . ; L and L ¼ 100 is taken. With the simulated data y1:L ¼ fy1 ; . . . ; yL g, we apply Algorithm 1 to estimate the simulated position, velocity and acceleration of the target, with N ¼ 500. The distribution of pðut Þ is set to 1 pðut Þ ¼ 45 if ut ¼ ð0; 0ÞT , and pðut Þ ¼ 58 otherwise. In Algorithm 1, M ¼ 5 samples are generated from pðut Þ every time. The results of the estimation are shown in Fig. 1 and the results of the delayed-weight estimation with D ¼ 10 are shown in Fig. 2. Since the position of target can be directly computed from the noisy observations. Fig. 1(f) presents the results of direct computation. But the velocity and acceleration cannot be directly computed. We perform G ¼ 100 runs of the simulations. All the computations are performed on Pentium 4 1.7 GHz PC with Matlab 6.0. The performance measure is root mean squared error (RMS) computed as follows [15]: ( L G 1 X 1 X 2 ðRMSÞ ¼ ½ðst1 ðgÞ b st1 ðgÞÞ2 L t¼1 G g¼1
,
where ðst1 ðgÞ; st2 ðgÞÞT is the simulated target position and ðb st1 ðgÞ; b st2 ðgÞÞT is the estimated target position at time t of the gth simulation. The results of performance are RMS ¼ 2:3977 for Algorithm 1, RMS ¼ 1:8737 for the delayed-weight estimation method, and RMS ¼ 2:7536 for the direct computation, respectively. 5.2. Multitarget tracking In this simulation study, the proposed MTT algorithm is implemented with five targets and Dt ¼ 1 s. An observation produced by the mth target is generated according to 0 1 ksm t pk2 m B s p2 C yt ¼ @ A þ Zt , tan1 t2 sm p 1 t1 2 0 where p ¼ ð1; 1ÞT and Zt N 00 ; ð1:0Þ . 0 ð0:04Þ2 The number of false alarm follows the Poisson distribution with mean lV ¼ 2. These false alarms are assumed to be independent and uniformly distributed within the observation volume V. With these settings, the states of the five targets and the observations are simulated with t ¼ 1; . . . ; L and L ¼ 100. The trajectories of the targets are represented in Fig. 3, where each ‘‘ ’’ indicates the starting position and each ‘‘ ’’ indicates the ending position of one target. From the simulated observations, the proposed MTT algorithm is implemented to estimate the states of the targets. The results are presented in Fig. 3(a). The delayed-weight estimation method is also applied here and the results are presented in Fig. 3(b). If two targets are in the approximately same position in the two-dimensional plane at the same time, it is usually difficult to distinguish them. Since our algorithm estimates the velocity of target, we can take advantage of the velocity information to distinguish them. Moreover, if the velocities of these two targets are also approximate, we can take advantage of the acceleration
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
201
10 200
5 0
150
–5 100 –10 50
–15
0
–20 –100
–50
0
50
100
150
(a)
0
10
20
30
40
50
60
70
80
90 100
0
10
20
30
40
50
60
70
80
90 100
(b) 15
2 1.5
10 1 0.5
5
0 0
–0.5 –1
–5
–1.5 –10
–2 0
10
20
30
40
50
60
70
80
90 100
(c)
(d) 300
2.5 2
250
1.5 200
1 0.5
150
0 100
–0.5 –1
50
–1.5 0
–2 0
(e)
10
20
30
40
50
60
70
80
90 100
–200 –150 –100
–50
0
50
100
150
(f)
Fig. 1. The simulated data (solid line) and estimated data by Algorithm 1 (dashed line): (a) target position; (b) x-coordinate of target velocity; (c) y-coordinate of target velocity; (d) x-coordinate of target acceleration; (e) y-coordinate of target acceleration; (f) the result of direct computation.
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203
202
10 200
5
150
0
100
–5
50
–10
0 –100
–50
0
50
100
–15
150
0
10
20
30
40
50
60
70
80
90
100
0
10
20
30
40
50
60
70
80
90 100
(b)
(a) 12
2
10
1.5
8 1
6 4
0.5
2
0
0
–0.5
–2
–1
–4 –1.5
–6
–2
–8 0
10
20
30
40
50
60
70
80
90
100 (d)
(c) 2 1.5 1 0.5 0 –0.5 –1 –1.5 0
10
20
30
40
50
60
70
80
90 100
(e) Fig. 2. The simulated data (solid line) and estimated data by delayed-weight estimation method (dashed line): (a) target position; (b) xcoordinate of target velocity; (c) y-coordinate of target velocity; (d) x-coordinate of target acceleration; (e) y-coordinate of target acceleration.
ARTICLE IN PRESS Y. Yu, Q. Cheng / Signal Processing 86 (2006) 195–203 800
800
700
700 4
600
500
400
400 5
300 200
3
200
0
3 2
0
1 –300 –200–100
5
300
100
2
0 –100
4
600
500
100
203
100 200 300 400 500 600 700
(a)
–100
1 –300 –200–100
0
100 200 300 400 500 600 700
(b)
Fig. 3. The simulated positions (solid line) and the estimated positions (dashed line) of the multitarget tracking with five targets by (a) Algorithm 2 and (b) delayed-weight estimation method.
information. But if the positions, velocities and accelerations of these two targets are all approximate, we cannot distinguish them any more. 6. Conclusions In this paper, algorithms are presented for the target tracking problem for the case of the maneuvering target. A suitable model is proposed to characterize the maneuvering acceleration. Algorithms based on particle filters are developed. Finally, two simulations are given to demonstrate the performance of the procedures. References [1] M. Hellebrandt, R. Mathar, Location tracking of mobiles in cellular radio networks, IEEE Trans. Veh. Technol. 48 (5) (September 1999) 1158–1162. [2] T. Liu, P. Bahl, I. Chlamtac, Mobility modeling, location tracking, and trajectory prediction in wireless ATM networks, IEEE J. Select. Areas Commun. 16 (6) (August 1998) 922–936. [3] D. Avitzour, Stochastic simulation Bayesian approach to multitarget tracking, IEE Proc.-F 142 (2) (April 1995) 41–44. [4] F. Gustafsson, F. Gunnarsson, N. Bergman, U. Forssell, J. Jansson, R. Karlsson, P. Nordlund, Particle filters for positioning, navigation and tracking, IEEE Trans. Signal Process. 50 (2) (February 2002) 425–437.
[5] A. Doucet, S.J. Godsill, C. Andrieu, On sequential Monte Carlo sampling methods for Bayesian filtering, Statist. Comput. 10 (3) (July 2000) 197–208. [6] J.S. Liu, R. Chen, Sequential Monte Carlo methods for dynamic systems, J. Amer. Statist. Assoc. 93 (443) (September 1998) 1032–1044. [7] N.J. Gordon, D.J. Salmond, A.F.M. Smith, Novel approach to nonlinear/non-Gaussian Bayesian state estimation, IEE Proc.-F 140 (2) (April 1993) 107–113. [8] N. Gordon, A hybrid bootstrap filter for target tracking in clutter, IEEE Trans. Aerosp. Electron. Syst. 33 (1) (January 1997) 353–358. [9] C. Hue, J.L. Cadre, P. Perez, Sequential Monte Carlo methods for multiple target tracking and data fusion, IEEE Trans. Signal Process. 50 (2) (February 2002) 309–325. [10] M.K. Pitt, N. Shephard, Filtering via simulation: auxiliary particle filters, J. Amer. Statist. Assoc. 94 (446) (June 1999) 590–599. [11] R. Chen, X.D. Wang, J.S. Liu, Adaptive joint detection and decoding in flat-fading channels via mixture Kalman filtering, IEEE Trans. Inform. Theory 46 (6) (September 2000) 2079–2094. [12] Y. Bar-Shalom, T.E. Fortmann, Tracking and Data Association, Academic Press, New York, 1988. [13] T.E. Fortmann, Y. Bar-Shalom, M. Scheffe, Sonar tracking of multiple targets using joint probabilistic data association, IEEE J. Oceanic Eng. OE-8 (July 1983) 173–184. [14] B. Zhou, N.K. Bose, Multitarget tracking in clutter: fast algorithm for data association, IEEE Trans. Aerosp. Electron. Syst. 29 (2) (April 1993) 352–363. [15] A. Doucet, N.J. Gordon, V. Krishnamurthy, Particle filters for state estimation of jump Markov linear systems, IEEE Trans. Signal Process. 49 (3) (March 2001) 613–624.