# Multiple Maneuvering Target Tracking by Improved Particle Filter Based on Multiscan JPDA

Authors: Jing Liu; ChongZhao Han; Feng Han; Yu Hu
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

The multiple maneuvering target tracking algorithm based on a particle filter is addressed. The equivalent-noise approach is adopted, which uses a simple dynamic model consisting of target state and equivalent noise which accounts for the combined effects of the process noise and maneuvers. The equivalent-noise approach converts the problem of maneuvering target tracking to that of state estimation in the presence of nonstationary process noise with unknown statistics. A novel method for identifying the nonstationary process noise is proposed in the particle filter framework. Furthermore, a particle filter based multiscan Joint Probability Data Association (JPDA) filter is proposed to deal with the data association problem in a multiple maneuvering target tracking. In the proposed multiscan JPDA algorithm, the distributions of interest are the marginal filtering distributions for each of the targets, and these distributions are approximated with particles. The multiscan JPDA algorithm examines the joint association events in a multiscan sliding window and calculates the marginal posterior probability based on the multiscan joint association events. The proposed algorithm is illustrated via an example involving the tracking of two highly maneuvering, at times closely spaced and crossed, targets, based on resolved measurements.


---

## Body

## 1. Introduction

The problem of multiple maneuvering target tracking has received a considerable attention in recent years. Multiple model approaches are proposed to model highly maneuvering targets and the interacting multiple model (IMM) algorithm [1–3] is the most popular one among them. Moreover, in the presence of clutter, some of the sensor measurements may not relate to the target of interest. The problem of data association has to be solved. Effective approaches in a Bayesian framework include the probabilistic data association (PDA) approach [4, 5] and the integrated probabilistic data association (IPDA) approach [6–8] for tracking a single target in clutter, the JPDA approach [4, 9, 10], and the joint integrated probabilistic data association (JIPDA) approach [11] for tracking multiple targets in clutter. The PDA and JPDA approaches address the track maintenance issue, while the IPDA and JIPDA approaches tackle the issues related to track initiation and track termination as well as track maintenance.Different combinations of the IMM method and the PDA/JPDA/IPDA/JIPDA methods have been used to tackle multiple maneuvering target tracking problems. In [12], the IMM algorithm is combined with the PDA filter in a multiple sensor scenario proposing a combined IMMPDA algorithm. The IMMPDA algorithm has good tracking performance when the targets are sufficiently separated. In [13], an IMMJPDA-coupled filter is developed for situations where the measurements of two targets are unresolved during periods of close encounter. The filters corresponding to the individual targets are coupled through cross-target-covariance terms. In [14], JPDA is combined with a crude approximation of IMM, under the name IMMJPDA. In [15], an IMMJIPDA algorithm is used in situations with a small number of crossing targets and low clutter measurement density. In order to deal with more complex scenarios with a large number of targets, an IMM-LMIPDA algorithm has been proposed by integrating linear multitarget (LM) method with IMM-IPDA. The computational requirements of the proposed IMM-LMIPDA algorithm increase linearly in the number of targets and the number of measurements.Recently, several methods which significantly outperform IMMJPDA have been proposed to deal with the multiple maneuvering target tracking problem. In [16], an IMMJPDA uncoupled fixed-lag smoothing algorithm is proposed, which performs far better than IMMJPDA filtering. In [17, 18], the descriptor system approach [19] is used to develop a track-coalescence-avoiding version of IMMJPDA. In [20], multiscan information is incorporated in JPDA to tackle the data association problem associated with the multiple maneuvering target tracking. In [21], a new track-coalescence-avoiding IMMCPDA method is proposed to improve the situation of tracking closely spaced targets by characterizing the exact Bayesian filter, and by developing novel target tracking combinations of IMM and PDA.In the above methods, the data models are assumed Gaussian and weakly nonlinear, and the Kalman filter/extended Kalman filter (KF/EKF) is used to perform target state estimation. More recently, nonlinear filtering techniques have been gaining more attention. The most common one among them is the particle filter [22], a state estimate pdf sampling based algorithm.There are not many reported works concerning the particle filter based multiple maneuvering target tracking methods. In [23], the problem of maintaining tracks of multiple maneuvering targets from unassociated measurements is formulated as a problem of estimating the hybrid state of a Markov jump linear system from measurements made by a descriptor system with independent, identically distributed (i.i.d.) stochastic coefficients. The possible multiple motion models and transition probability matrices are assumed to be known in the proposed method. In [24], the marginal filtering distributions for each of the targets are approximated with particles. The posterior marginal probability is calculated based on the joint association hypotheses, which are examined in a single scan. However, single scan information may not be enough for tracking maneuvering targets, especially in the condition of tracking two closely spaced targets.In this work, a novel method is proposed for a multiple maneuvering target tracking, which is a combination of the particle filter based process noise identification algorithm and the particle filter based multiscan JPDA algorithm. The process noise identification process is effective in estimating both the maneuvering movement and the random acceleration of the target, using one general model for each of the maneuvering targets. The multiscan JPDA is effective in maintaining the tracks of multiple targets using multiple scan information. Compared with the single scan JPDA method, the multiscan JPDA method uses richer information, which results in better estimated probabilities.The rest of the paper is organized as follows. Section2 provides a brief introduction to particle filter, and the multiple maneuvering target tracking problem is defined in Section 3. The particle filter based process noise identification method for tracking highly maneuvering target is introduced in Section 4. The data association method based on multiscan JPDA is discussed in Section 5. The proposed multiple maneuvering target tracking algorithm, which is a combination of the above two algorithms, is introduced in Section 6. The simulation results are shown in Section 7 and the paper is summarized in Section 8.


## 2. Basic Theory of Particle Filter

To define the problem of tracking, consider a dynamic system represented by the state equation:(2.1)xk=f(xk-1,vk-1*),
where x is the state, f is a possibly nonlinear function, and v* is the known process noise with a zero mean Gaussian distribution. The objective of tracking is to recursively estimate xk from a sequence of measurements up to time step k, z1:k={z1,z2,…,zk}. The observation model is described as follows:
(2.2)zk=h(xk,nk),
where h is a possibly nonlinear function. n is the observation noise with a zero mean Gaussian distribution. From the Bayesian perspective, the tracking problem is to recursively calculate the posterior distribution p(xk∣z1:k).In this paper, a particle filter is considered to solve the state estimation problem due to its ability to tackle the nonlinear and non-Gaussian systems. The posterior distributionp(xk∣z1:k) is approximated by a set of particles with associated weights. The detailed introduction about particle filter algorithm can be found in [22]. The procedures associated with the standard particle filter are listed in the following.Algorithm 2.1 (particle filter algorithm).
(1)Initialization. Sample initial particles {x0i,i=1,…,H} from the initial posterior distribution p(x0) and set the weights w0i to 1/H, i=1,…,H. H is the number of particles.(2)Prediction. Particles at time step k-1, {xk-1i,i=1,…,H}, are passed through the system model (2.3) to obtain the predicted particles at time step k, {x^ki,i=1,…,H}:
(2.3)x^ki=f(xk-1i,vk-1*,i),
where vk-1*,i is a sample drawn from the known zero mean Gaussian distribution.(3)Update. Once the observation data, zk, is measured, evaluate the importance weight of each predicted particle and obtain the normalized weight for each particle:
(2.4)w^ki=p(zk∣x^ki),wki=w^ki∑i=1Hw^ki.
Thus at time step k we can obtain the estimate of the state, x~~k=∑i=1H‍wkix^ki.(4)Resample. Resample the discrete distribution {wki:i=1,…,H}, H times to generate particles {xkj:j=1,…,H}, so that for any j, Pr{xkj=x^ki}=wki. Set the weights wki to 1/H,  i=1,…,H, and move to Stage 2.


## 3. Multiple Maneuvering Target Tracking Problem Formulation

The number of maneuvering targets (M) to be tracked is assumed as fixed and known, where each target track has been initiated, and our objective is to maintain the tracks. Each target is parameterized by a state xm,k, where m denotes the mth target and k denotes time step k. The combined state, xk=(x1,k,…,xM,k), is the concatenation of the individual target states. The individual targets are assumed to evolve independently according to a general motion model,
(3.1)xk=f(xk-1,uk-1,vk-1*),
where u is the maneuver acceleration and v* is the process noise.The observation vectorzk is composed of multiple sensor measurements {zj,k,j=1,…,Nk}, where Nk is the total number of measurements. It is assumed that there are no unresolved measurements (i.e., measurements associated with two or more targets simultaneously); any measurement is either associated with a single target or caused by clutter. Clutter is modeled as independently and identically distributed (IID) with an uniform spatial distribution over the surveillance area.In this paper, the multiple maneuvering target tracking problem is divided into two parts: single maneuvering target tracking and multiple target tracking, which are introduced in Sections4 and 5, respectively. The particle filter based process noise identification method for tracking a highly maneuvering target is introduced in Section 4. The data association method based on multiscan JPDA is discussed in Section 5. Finally, in Section 6, the combination of the two methods is combined to deal with the multiple maneuvering target tracking problem.


## 4. Maneuvering Target Tracking Based on Process Noise Identification

In this section, the particle filter based process noise identification method is proposed for tracking highly maneuvering target. The general motion model of a maneuvering target is described by (3.1). In the equivalent-noise approach [25–27], it is assumed that the general motion model (3.1) that describes target motions can be simplified to
(4.1)xk=f(xk-1,vk-1),
with an adequate accuracy, where v is the equivalent noise that quantifies the error of this model in describing the target motions, in particular, maneuvers. The statistics of this noise v, nonstationary in general, are not known.In this section, a novel method is proposed for process noise identification. The process noise is modeled as a dynamic system. The noise vectorvk-1 is chosen as the state of the noise system. The observation vector is zk, which is the same as in (2.2). The observation equation is defined in
(4.2)zk=h(xk,nk)=h[f(xk-1,vk-1),nk].
The aim of the proposed method is to estimate the posterior distribution of the process noise, p(vk-1∣z1:k). According to the Bayesian inference theory,
(4.3)p(vk-1∣z1:k)=p(zk∣vk-1)p(vk-1∣z1:k-1)p(zk∣z1:k-1),
where p(zk∣z1:k-1) is a normalizing constant and 1/p(zk∣z1:k-1) is defined as Υ. So we can obtain
(4.4)p(vk-1∣z1:k)=Υ·p(zk∣vk-1)p(vk-1∣z1:k-1).
The noise vk-1 may be from random accelerations, sudden maneuvers, or both, and there is no information about the distribution of vk-1. vk-1 is not dependent on the previous measurements z1:k-1, which results in
(4.5)p(vk-1∣z1:k)=Υ·p(zk∣vk-1)p(vk-1).
According to the Bayesian theory, the prior distribution of parameters, on which there is no information, could be considered as a uniform distribution. In this paper, the prior distribution of the process noise is considered as a uniform distribution, and the posterior distribution of the process noise is then obtained via a Monte Carlo deviation process. The general procedure of the particle filter based process noise identification algorithm is as follows.Initially a number of process noise samples are generated from a uniform distribution since there is no information on the process noise distribution due to the uncertain dynamics. Each of the process noise samples is assigned with a weight based on the current measurement and the particles from the previous time step. According to Monte Carlo theory, the process noise samples and their associated weights approximate the true posterior distribution of the process noise. The process noise samples are then resampled according to their associate weights. The samples with large weights are duplicated, while the samples with small weights are removed. The resampled process noise samples distribute approximately according to the true posterior distribution of the process noise.Since there is no information aboutvk-1, it is reasonable to assume that vk-1 is uniformly distributed in the range of U(-d,d), where U denotes a uniform distribution and d is the known process noise bound accounting for the maximin uncertain dynamics. According to the Monte Carlo theory, p(vk-1) can be represented by H samples, {v^k-1j,j=1,…,H}, from U(-d,d). Consider(4.6)p(vk-1)=1H∑j=1Hδ(vk-1-v^k-1j).
The number of process noise samples (H) is proportional to the magnitude of the maximum uncertain dynamics (∥d∥). The posterior distribution of vk-1 can be represented as
(4.7)p(vk-1∣z1:k)=ΥH∑j=1Hp(zk∣v^k-1j)δ(vk-1-v^k-1j)=ΥH∑j=1Hξkj·δ(vk-1-v^k-1j),
where ξkj=p(zk∣v^k-1j) is defined as the weight assigned to the jth process noise sample v^k-1j. The process noise samples {v^k-1j,j=1,…,H} are then resampled according to {ξkj,j=1,…,H} based on the principle thatPr{vk-1i=v^k-1j}=ξkj, where {vk-1i,i=1,…,H} are the process noise samples obtained from resampling. The obtained resampled process noise samples are approximately distributed according to the posterior distribution p(vk-1∣z1:k).In order to calculateξkj, the likelihood function p(zk∣v^k-1j) is expanded based on the resampled state particles at time step k-1, {xk-1i,i=1,⋯,H}. Consider(4.8)p(zk∣v^k-1j)=∑i=1Hp(zk∣xk-1i,v^k-1j)p(xk-1i∣v^k-1j).
Since xk-1i and v^k-1j are independent, p(xk-1i∣v^k-1j)=p(xk-1i).The resmapled particles at time stepk-1, {xk-1i,i=1,…,H}, should be assigned with the same and equal weights, 1/H. We can then obtain
(4.9)p(xk-1i)=1H.To calculatep(zk∣xk-1i,v^k-1j), define μki,j as the intermediate particle,
(4.10)μki,j=f(xk-1i,v^k-1j),i=1,…H,j=1,…,H
and expand p(zk∣xk-1i,v^k-1j) based on μki,j,
(4.11)p(zk∣xk-1i,v^k-1j)=∑p=1H∑q=1H[p(zk∣μkp,q,xk-1i,v^k-1j)p(μkp,q∣xk-1i,v^k-1j)].
Since xk-1i and v^k-1j are known and μkp,q is obtained from a purely deterministic relationship in (4.10), we obtain
(4.12)p(μkp,q∣xk-1i,v^k-1j)={1,forp=iandq=j,0,forp≠iorq≠j,(4.13)p(zk∣xk-1i,v^k-1j)=p(zk∣μki,j).
Combining (4.9) and (4.13) with (4.8) results in
(4.14)p(zk∣v^k-1j)=∑i=1H‍p(zk∣μki,j)·1H.
At each time step, the process noise samples are drawn from a uniform distribution U(-d,d). Each process noise sample v^k-1j is evaluated and assigned its corresponding weight ξkj. A resampling procedure is then used to redistribute the process noise samples, from which the process noise samples with large weights are replicated while the samples with small weights are eliminated.The standard particle filter procedure for state estimation follows next. The predicted particles{x^ki:i=1,…,H} are then obtained based on the resampled process noise samples {vk-1i:i=1,…,H} through the dynamic model (4.1). The predicted particles are updated and resampled as in the conventional particle filter algorithm.The complete algorithm including the process noise estimation and state estimation parts is summarized below.Algorithm 4.1 (particle filter based process noise identification).
(
1
)  At time step k-1, draw process noise samples {v^k-1j:j=1,…,H} from a uniform distribution U(-d,d).
(
2
)  Calculate the intermediate particles {μki,j:i=1,…,H;j=1,…,H} according to (4.10).
(
3
)  Calculate the process noise sample weights {ξkj:j=1,…,H} via (4.14).
(
4
)  Resample the discrete distribution {ξkj:j=1,…,H}, H times to generate the resampled process noise samples {vk-1i:i=1,…,H}, so that for any i, Pr{vk-1i=v^k-1j}=ξkj. Set the weights ξkj to 1/H,i=1,…,H.
(
5
)  Calculate the predicted particles at time step k, {x^ki,i=1,…,H},
(4.15)x^ki=f(xk-1i,vk-1i).(6) and (7) are the same with the stages (3) and (4) in Algorithm 2.1.Algorithm 4.2 (simplification of the proposed algorithm).
In the proposed algorithm, at each iteration,H*H intermediate particles are calculated through the permutation of particles and process noise samples in (4.10) and evaluated via (4.13). This increases the computation burden and the algorithm runs slowly compared to the conventional particle filter, which are based on H particles (Algorithm 2.1). It is observed that at each time step, after resampling the particles focus in some smaller area, a large portion of particles are with the same value (Algorithm 2.1, Step (4)). To simplify the algorithm, it is assumed that the particles (Algorithm 2.1, Step (4)) are less variable compared with the process noise samples (4.10). In (4.10), particles {xk-1i,i=1,…,H} are replaced by x~~k-1, the estimate of the state at time step k-1, which results in a simplified version of (4.10),(4.16)μkj=f(x~~k-1,v^k-1j),
and p(zk∣v^k-1j) is expanded directly on μkj,
(4.17)p(zk∣v^k-1j)=∑τ=1Hp(zk∣μkτ)p(μkτ∣v^k-1j).
Using the similar derivation process with (4.12), we can obtain
(4.18)p(zk∣v^k-1j)=p(zk∣μkj).
In the simplified version of the proposed algorithm, the number of intermediate particles is reduced toH, which reduces the computation burden and increases the computing speed. More importantly, the performance of the algorithm with the simplification procedure is verified through simulation study. In the following sections, the particle filter based process noise identification algorithm refers to the simplified version.Discussion.
It is required to set the process noise boundd preliminarily, which could be estimated based on the maximin uncertain dynamics. For most of the manned air vehicles, the process noise bound d could be determined according to the flight envelop of the aircraft. For the unmanned air vehicles with unexpected maneuvers, d is set large enough to cover the unexpected maneuvers. If our knowledge about d is imprecise, tracking process may fail when coping with the maneuvers that are not covered by the process noise bound d.


## 5. Multiple Target Tracking Using Particle Filter Based Multiscan JPDAF

The number of targets (M) to be tracked is assumed as fixed and known, where each target track has been initiated, and our objective is to maintain the tracks. Each target is parameterized by a state xm,k, where m denotes the mth target and k denotes time step k. The combined state, xk=(x1,k,…,xM,k), is the concatenation of the individual target states. The individual targets are assumed to evolve independently according to the Markovian dynamic models pm(xm,k∣xm,k-1). The observation vector zk is composed of multiple sensor measurements {zj,k,j=1,…,Nk}, where Nk is the total number of measurements. It is assumed that there are no unresolved measurements (i.e., measurements associated with two or more targets simultaneously); any measurement is either associated with a single target or caused by clutter. Clutter is modeled as independently and identically distributed (IID) with uniform spatial distribution over the surveillance area.In the particle filter based JPDA algorithm, the distributions of interest are the marginal filtering distributions for each of the targetspm(xm,k∣z1:k),m=1⋯M, and these distributions are approximated with particles, {x^m,ki,i=1⋯H,m=1⋯M}, and their associated weights {wm,ki,i=1⋯H,m=1⋯M}, as in
(5.1)pm(xm,k∣z1:k)=∑i=1Hwm,kiδ(xm,k-x^m,ki).
At each time step, when the new observation vector arrives, the marginal filtering distributions for each of the targets are updated through the Bayesian sequential estimation recursions [24].In the standard single scan JPDA framework, a track is updated with a weighted sum of the measurements which could have reasonably originated from the target in track. The only information the standard JPDA algorithm uses is the measurements on the present scan and the state vectors. If more scans of measurements are used, additional information is available resulting in better computed probabilities [4]. Since the tracking systems are unable to store all of the measurements from all the scans, a Bayesian tracking system can at best rely on a sliding window of scans. In [28], the single scan JPDA filter is extended to the multiple scan JPDA filter for tracking multiple targets. And in this work the multiple scan JPDA filter is developed in a particle filter framework.The multiple scan JPDA calculation examines multiple scan joint association events [28]. The measurement to target the association event of a multiple scan is defined as λk-L+1:k, where L denotes the length of the multiple scan sliding window. The multiple scan joint association events are mutually exclusive, and they form a complete set Λk-L+1:k [11]. λk-L+1:k is composed by the association vectors at each scan in the sliding window, λk-L+1:k=(θk-L+1,θk-L+2,…,θk). The elements of the association vector at time step k, θk=(ζ1,k,…,ζj,k,…,ζNk,k) are given by
(5.2)ζj,k={0,ifzj,kisduetoclutter,m∈{1⋯M},ifzj,kisfromtargetm.
The heart of the new algorithm is to find the posterior probability for the joint association event of multiple scans. That is to calculate p(λk-L+1:k∣z1:k) and it can be written as
(5.3)p(λk-L+1:k∣z1:k)∝p(zk⋯zk-L+1∣λk-L+1:k,z1:k-L)p(λk-L+1:k∣z1:k-L)∝p(zk⋯zk-L+1∣λk-L+1:k,z1:k-L)p(λk-L+1:k),
where the conditioning of λk-L+1:k on the history of measurements before the sliding window has been eliminated.The distribution of the measurements in the sliding window based on a specific association event is given by(5.4)p(zk⋯zk-L+1∣λk-L+1:k,z1:k-L)=∏s=1L[∏j=1Nk-L+sp(zj,k-L+s∣λk-L+1:k,z1:k-L)].
To reduce the notation, the index of the scan s in the sliding window is denoted by ks=k-L+s. Then we can obtain
(5.5)p(zk⋯zk-L+1∣λk-L+1:k,z1:k-L)=∏s=1L[∏j=1Nksp(zj,ks∣λk-L+1:k,z1:k-L)]=∏s=1L[∏j∈I0,kspc(zj,ks)·∏j∈Iksp(zj,ks∣xζj,ks,ks)]=∏s=1L[(V)-Cks·∏j∈Iksp(zj,ks∣xζj,ks,ks)],
where I0,ks={j∈{1,…,Nks}:ζj,ks=0} and Iks={j∈{1,…,Nks}:ζj,ks≠0} are, respectively, the subsets of measurement indices corresponding to clutter measurements and measurements from the targets being tracked, on scan ks. pc denotes the clutter likelihood model, which is assumed to be uniform over the volume of the surveillance area V. The volume of the surveillance area could be calculated as per V=2πRmax, where Rmax is the maximum range of the sensor. Cks is defined as the number of clutter measurements.The joint association priorp(λk-L+1:k), can be calculated as in (5.6) according to [4, 16],
(5.6)p(λk-L+1:k)=∏s=1L[Cks!εNks!∏m=1M(PD)δm(θks)(1-PD)1-δm(θks)],
where ε is a “diffuse” prior [16] and PD is the detection probability. δm(θks) is a binary variable and set to one if the mth target is assigned with a measurement in the event θks.The posterior probability for the joint association event of multiple scans is obtained as(5.7)p(λk-L+1:k∣z1:k)∝p(λk-L+1:k)∏s=1L[(V)-Cks·∏j∈Iksp(zj,ks∣xζj,ks,ks)].
The posterior probability that the jth measurement is associated with the mth target on scan k, βjm, is calculated by summing over the probabilities of the corresponding joint association events via
(5.8)βjm=p(ζj,k=m∣z1:k)=∑{λk-L+1:k∈Λk-L+1:k:ζj,k=m}p(λk-L+1:k∣z1:k).These approximations can, in turn, be used in (5.9) to approximate the target likelihood according to [24],
(5.9)pm(zk∣xm,k)=β0m+∑j=1Nkβjmp(zj,k∣xm,k),
where β0m is the posterior probability that the mth target is undetected. Finally, setting the new importance weights to
(5.10)wm,ki∝wm,k-1ipm(zk∣x^m,ki)pm(x^m,ki∣xm,k-1i)qm(x^m,ki∣xm,k-1i),∑i=1H‍wm,ki=1,
where qm(x^m,ki∣xm,k-1i) is the proposal distribution, which is used to generate the predicted particles. Equation (5.10) leads to the sample set {wm,ki,x^m,ki}i=1H being approximately distributed according to the marginal filtering distribution pm(xm,k∣z1:k).A summary of the particle filter based multiple scan JPDA filter algorithm is presented in what follows. Assuming that the sample sets{wm,k-1i,xm,k-1i}i=1H,m=1⋯M, are approximately distributed according to the corresponding marginal filtering distributions at the previous time step pm(xm,k-1∣z1:k-1),m=1⋯M, the algorithm proceeds as follows at the current time step.Algorithm 5.1 (particle filter based multiscan JPDA filter).
(
1
)  For m=1⋯M,i=1⋯H, generate predicted particles for the target states x^m,ki~qm(x^m,ki∣xm,k-1i).
(
2
)  For m=1⋯M, calculate x~m,k, the preapproximation of xm,k, which is to be substituted into (5.7) to calculate the posterior probability of the joint association event in multiple scan, p(λk-L+1:k∣z1:k). Consider(5.11)x~m,k=∑i=1Hwm,k-1ix^m,ki.
(
3
)  Enumerate all the valid joint measurement to target association events in the sliding window k-L+1:k to form the set Λk-L+1:k.
(
4
)  For each λk-L+1:k∈Λk-L+1:k, compute the posterior probability of the joint association event in multiple scan, p(λk-L+1:k∣z1:k), via (5.7).
(
5
)  For m=1⋯M,j=1⋯Nk, compute the marginal association posterior probability, βjm, via (5.8).
(
6
)  For m=1⋯M,i=1⋯H, compute the target likelihood, pm(zk∣x^m,ki), via (5.9).
(
7
)  For m=1⋯M,i=1⋯H, compute and normalize the particle weights, via (5.10).
(
8
)  Resample the discrete distribution {wm,ki:i=1,…,H}, H times to generate particles {xm,kj:j=1,…,H}, so that for any j, Pr{xm,kj=x^m,ki}=wm,ki. Set the weights wm,ki to 1/H,i=1,…,H, and move to Stage 1.


## 6. An Algorithm for Tracking Multiple Maneuvering Targets

The proposed multiple maneuvering target tracking algorithm is a combination of the particle filter based process noise identification algorithm (proposed in Section4) and the particle filter based multiscan JPDA algorithm (proposed in Section 5). In the proposed algorithm, firstly the particle filter based process noise identification algorithm is used to estimate the maneuvering movement for each target. Then the particles of each target model are propagated to the next time step based on the new distributed process noise samples to obtain the predicted particles. In the update process, each predicated particle of one target model is assigned with a weight based on the multiscan JPDA process. The steps of the proposed multiple maneuvering target tracking algorithm are listed in the following.Algorithm 6.1 (multiple maneuvering target tracking).
(
1
)  At time step k-1, for m=1⋯M,(a)
draw process noise samples{v^m,k-1j:j=1,…,H} from a uniform distribution U(-d,d),(b)
calculate the intermediate particles{μm,kj:j=1,…,H} according to
(6.1)μm,kj=f(x~~m,k-1,v^m,k-1j),(c)
calculate the process noise sample weights{ξm,kj:j=1,…,H} via (6.2) and normalize each weight:
(6.2)ξm,kj=p(zkm∣μm,kj),
where zkm are the measurements that are close to μm,kj obtained from the nearest neighbor method,(d)
resample the discrete distribution{ξm,kj:j=1,…,H}, H times to generate the new process noise samples {vm,k-1i:i=1,…,H}, so that for any i, Pr{vm,k-1i=v^m,k-1j}=ξm,kj. Set the weights ξm,kj to 1/H,i=1,…,H,(e)
obtain the predicted particles{x^m,ki:i=1,…,H} at time step k from the new process noise samples via
(6.3)x^m,ki=f(xm,k-1i,vm,k-1i).
(
2
)  At time step k, for m=1⋯M, calculate x~m,k, the postapproximation of xm,k,
(6.4)x~m,k=∑i=1H‍ξm,kix^m,ki,
where {ξm,ki:i=1,…,H} are the process noise samples' weights and are calculated based on the measurements at the current time step (k). However, in (5.11) x~m,k is calculated based on {wm,k-1i:i=1,…,H}, which rely on the measurements at the previous time step (k-1). As a result, in maneuvering target tracking, x~m,k calculated via (6.4) is closer to the true state than that via (5.11) since the information from the current time step is considered in (6.4).
(
3
)–(8) are the same with the stages (3)–(8) in Algorithm 5.1.
(
3
)
Thus we can obtain x~~m,k=∑i=1Hwm,kix^m,ki.  x~~m,k is the estimate of the true state xm,k.


## 7. Simulation Results and Analysis

This section consists of three parts. In the first and second part, the particle filter based process noise identification algorithm and the particle filter based multiscan JPDA algorithm are, respectively, used to track single maneuvering targets and two slow-maneuvering targets in clutter. In the third part, the combination of the two algorithms are used to track two highly maneuvering, at times closely spaced and crossed, targets.


### 7.1.  Single Maneuvering Target Tracking

The simulation study using nearly coordinate turn model is performed. The maneuvering target tracking is done by setting up a 2D flight path inx-y plane, which is similar to the path considered in [16]. At time step 1, the target starts at location [-310,310] in Cartesian coordinates in meters with initial velocity (in m/s) [10,-400]. The following trajectory is considered: a straight line with constant velocity between 1 and 17 s, a coordinated turn (0.09 rad/s) between 17 and 34 s, a straight line with constant velocity between 34 and 51 s, a coordinated turn (0.09 rad/s) between 51 and 68 s, and a straight line with constant velocity between 68 and 100 s.In the particle filter with process noise identification, a general model(7.1)Xk=ΦXk-1+Γvk-1
is adopted during the whole tracking process, where in (7.1),
(7.2)Φ=[Φb00Φb],Φb=[1ΔTΔT2201ΔT001],Γ=I6×6.
Matrix Φ is the transition matrix and ΔT is the sample interval. Xk=[px,γx,ax,py,γy,ay]kT is the state vector; px,γx, and ax denote, respectively, the position, velocity, and acceleration of the moving object along the x-axis of Cartesian frame and, py,γy, and ay along the y-axis. The equivalent process noise, vk-1=[vpx,vγx,vax,vpy,vγy,vay]k-1T, with unknown statistics is required to be identified. The bound of the process noise (d), which accounts for the uncertain dynamics, is chosen as {20m,20m/s,10m/s2,20m,20m/s,10m/s2}. The number of the process noise samples is equal to the number of particles, which is set as 500. The algorithm is initialized with Gaussian around the initial state of the true target, and the standard deviation of the Gaussian distribution is chosen as {10m,10m/s,5m/s2,10m,10m/s,5m/s2}.A track-while-scan (TWS) radar is positioned at the origin of the plane. The measurement equation is as follows:(7.3)Zk=h(Xk)+nk,
where Zk=[z1,z2]k is the observation vector. z1 is the distance between the radar and the target, and z2 is the bearing angle. The measurement noise nk=[nz1,nz2]k is a zero mean Gaussian white noise process with standard deviations of 20m(σz1) and 0.01rad(σz2). Resolution of the sensor is selected after from [29] (twice of the standard deviations of the measurement noise). The sampling interval is ΔT=1s.The particle filter based process noise identification algorithm is compared to the IMM filter and the regularized particle filter. The IMM filter consists of three extended Kalman filter (EKF) with different motion models. The details regarding these models may be found in [16]. The initial model probabilities and the mode switching probability matrix are set the same values as in [16]. For the regularised particle filter, Epanechnikov kernel is chosen as the rescaled kernel density, which is the same as in [30]. Moreover, the proposed algorithm is compared to its complete version in the same simulation setup.The simulation results are obtained from1000 Monte Carlo runs except the results from the complete version of the particle filter based multiscan JPDA algorithm, which is run for 100 times. Figure 1 shows the true trajectory of the maneuvering target. The root mean-square errors (RMSEs) in position at each time step, respectively, using the four methods are shown in Figure 2, where RPF and PFPNI represent the regularised particle filter algorithm and the particle filter based process noise identification algorithm. The performance of the methods is also compared via the global RMSE (in position), the tracking loss rate (TLR), and the executing time (ET), which are listed in Table 1. The tracking loss rate (TLR) is defined as the ratio of the number of simulations, in which the target is lost in track, to the total number of simulations carried out. The target is defined as lost in track when its global RMSE in position is larger than the ten times of the magnitude of the standard deviation in position. The executing time (ET) is the CPU time needed to execute one time step in MATLAB 7.1 on a 3 GHz (Mobile) Pentium IV operating under Windows 2000.Table 1
Performance comparison.
RMSE (m)
ET (s)
TLR
IMM
28.34
0.0239
3.7%
RPF
16.66
0.8708
0
PFPNI (simplified)
13.67
0.4503
0
PFPNI (complete)
13.29
190.47
0Figure 1
True trajectory of the single maneuvering target.Figure 2
RMSE in position using IMM, RPF, PFPNI (simplified), and PFPNI (complete) algorithms.From the simulation results, it can be seen that the simplified version of the proposed particle filter based process noise identification algorithm outperforms the IMM filter and the regularised particle filter algorithm, with computing time within the limits of practically realizable systems. Moreover, the proposed algorithm needs neither the possible multiple motion models nor the transition probability matrices, which makes it a more general algorithm for maneuvering target tracking. From the simulation results, it can also be seen that the complete version of the particle filter based process noise identification algorithm is not suitable for practical application due to the long computing time, though it gains a2.8% increase in accuracy (RMSE) compared with its simplified version.


### 7.2. Multiple Target Tracking

The simulation is carried out for tracking two slow-maneuvering targets in clutter. The Wiener process acceleration model (7.1) is chosen as the motion model for the two targets. The process noise, vk-1=[vpx,vγx,vax,vpy,vγy,vay]k-1T, is a zero mean Gaussian white noise process with standard deviations of 1m(σvpx), 1m/s(σvγx), 20m/s2(σvax), 1m(σvpy), 1m/s(σvγy), and 20m/s2(σvay).At time step 1, target one starts at location[-310,310] in x-y Cartesian coordinates in meters with the initial velocity (in m/s) [10,-400]. Target two starts at location (in m) [-310,-19000] with the initial velocity (in m/s) [10,400]. The length of each simulation run is 50 seconds. A TWS radar is positioned at the origin of the plane, whose details are provided in Section 7.1.The sampling interval isΔT=1s and it is assumed that the probability of detection PD=0.9 for the radar. For generating measurements in simulations, the clutter is assumed uniformly distributed with density 1×10-6/m2.In the particle filter based multiscan JPDA methods, each target model is assigned with 500 particles. The algorithm is initialized with Gaussians around the initial states of the true targets, and the standard deviations of the two Gaussian distributions are chosen equally as{10m,10m/s,5m/s2,10m,10m/s,5m/s2}.In the simulations carried out, the length of the multiple scan sliding window (L) varies from 1 to 3. The corresponding particle filter based multiscan JPDA methods with L=1,2,3 are utilized to track two slow-maneuvering targets individually. A comparison to the standard JPDA filter is also studied on the same simulation setup. The simulation results are obtained from 1000 Monte Carlo runs. Figure 3 shows the true trajectories of the two targets and Figure 4 shows the distance between the two targets along time, through which we can see that the two targets reach the smallest distance at time step 25. The RMSEs in position for the two targets are, respectively, shown in Figures 5 and 6, where PFJPDA represents the particle filter based JPDA filter with its following number in the bracket (e.g., (2.1)) denoting the scan number (L). The performance of the four methods is also compared in Table 2. The swap rate (SR) is defined as the ratio of the number of simulations, in which the two targets swap, to the total number of simulations.Table 2
Performance comparison.
RMSE (m)
ET (s)
TLR
SR
JPDA
T1: 20.24, T2: 22.87
0.07
24%
0
PFJPDA(1)
T1: 13.21, T2: 11.36
1.28
0
0
PFJPDA(2)
T1: 12.92, T2: 11.01
1.54
0
0
PFJPDA(3)
T1: 12.75, T2: 10.67
1.8
0
0Figure 3
True trajectories of the two slow-maneuvering targets.Figure 4
Distance between the two targets (- distance;-* zoom in distance).Figure 5
Target 1: RMSE in position using standard JPDA filter, particle filter based multiscan JPDA algorithms with a scan number equals to 1, 2, and 3, respectively.Figure 6
Target 2: RMSE in position using standard JPDA filter, particle filter based multiscan JPDA algorithms with a scan number equals to 1, 2, and 3, respectively.Compared with the standard JPDA method (based on extended Kalman filter), the particle filter based JPDA methods (single scan and multiple scan) are much more accurate and robust, at the cost of longer computing time. This verifies that when dealing with nonlinear problem (nonlinear observation equation) and large random acceleration (large process noise), the performance of particle filter is better than extended Kalman filter using local linearization.In the comparison between the three particle filter based JPDA methods (L=1,2,3, resp.), from Figures 5 and 6 it can be seen that there is no significant deference between the three methods except around time step 25, when the two targets are very close to each other (Figure 4). As the scan number increases, the RMSE (in position) of the corresponding algorithm decreases significantly around time step 25 (Figures 5 and 6). This verified that the particle filter based multiscan JPDA method provides a better performance especially when the targets are very close. The additional information of more scans improve the association probabilities in such critical situation, resulting in lower estimation errors (RMSE) and larger robustness (tracking loss rate).


### 7.3. Multiple Maneuvering Target Tracking

We now consider tracking two highly maneuvering targets in clutter. The true trajectories are shown in Figure7 in x-y plane. The velocities of the targets are shown in Figure 8, where VX and VY represent the velocity, respectively, in x and y coordinate. The distance between the two targets is shown in Figure 9. Target one is with the same setup as in the example proposed in Section 7.1. Target two starts at location [-310,-20310] in Cartesian coordinates in meters with the initial velocity (in m/s) [10,400]. Its trajectory is a straight line with constant velocity between 1 and 17 s, a coordinated turn (0.09rad/s) between 17 and 34 s, and a straight line constant velocity between 34 and 100 s.Figure 7
True trajectories of the two maneuvering targets.Figure 8
True velocities of the two maneuvering targets.Figure 9
Distance between the two targets (- distance;-* zoom in distance).A TWS radar is positioned at the origin of the plane (refer to Section7.1). The sampling interval (ΔT), the probability of detection (PD), and the clutter density are the same as in Section 7.2.In the proposed method, each target model is assigned with 1000 particles. The length of the multiple scan sliding window (L) is chosen equal to 3. The bound of the process noise (d) is the same as in Section 7.1. The algorithm is initialized with Gaussians around the initial states of the true targets, and the standard deviations of the two Gaussian distributions are chosen equally as {10m,10m/s,5m/s2,10m,10m/s,5m/s2}.In this work, the proposed method is compared to an IMMJPDA method and a multiscan IMMJPDA method, whose details could be found in [16, 20], respectively.The simulation results are obtained from1000 Monte Carlo runs. The RMSEs (resp. in position and velocity) at each time step for the two targets are shown in Figures 10, 11, 12, and 13, where MS-IMMJPDA and PFPNI-PFMSJPDA represent, respectively, the multiscan IMMJPDA algorithm and the proposed algorithm. The performance comparison is listed in Table 3. Moreover, for the proposed algorithm, the influence of the particle number in its performance is studied and simulations are carried out based on different sample sizes. The results are listed in Table 4.Table 3
Performance comparison.
RMSE in position (m)
RMSE in velocity (m/s)
ET (s)
TLR
SR
IMMJPDA
T1: 26.04, T2: 31.28
T1: 39.48, T2: 31.98
0.75
8.7%
13.9%
MS-IMMJPDA
T1: 25.30, T2: 30.26
T1: 39.22, T2: 31.76
3.78
0
4.5%
PFPNI-PFMSJPDA
T1: 14.42, T2: 18.14
T1: 34.44, T2: 30.58
4.0
0
3.7%Table 4
Influence of particle number in the performance of the proposed algorithm for tracking multiple maneuvering target.
Particle number
RMSE in position (m)
ET (s)
TLR
200
NA
NA
100%
500
T1: 17.37, T2: 21.22
2.9
35.9%
1000
T1: 14.42, T2: 18.14
4.0
0
2000
T1: 13.75, T2: 17.93
20.2
0Figure 10
Target 1: RMSE in position using IMMJPDA filter, multiple scan IMMJPDA filter, and the proposed method.Figure 11
Target 2: RMSE in position using IMMJPDA filter, multiple scan IMMJPDA filter, and the proposed method.Figure 12
Target 1: RMSE in velocity using IMMJPDA filter, multiple scan IMMJPDA filter, and the proposed method.Figure 13
Target 2: RMSE in velocity using IMMJPDA filter, multiple scan IMMJPDA filter, and the proposed method.The simulation results show that the proposed algorithm is more accurate than the multiple scan IMMJPDA filter and IMMJPDA filter, though it takes longer computing time. However, the three algorithms are implemented using Matlab. The computing time is considerably reduced when coded in C++. For the proposed method, it takes about 0.8 second for one time step in C++.In the scene of simulation, at the period of46s~52s, the two targets are very near to each other, which easily results in track swap. In the IMMJPDA method, the judgement of which measurements belong to which target depends on the information from the current scan. If the measurements from two targets are very close, it is hard to distinguish the different targets based on the measurements from a single scan, which leads to a track swap. In the proposed method and the multiple scan JPDA method, the information from several previous scans is combined with the information from current scan to calculate the association probabilities. The decision of which measurements belong to which target based on the information from multiple scans will be more accurate than single scan method, which reduces the swap rate effectively.From Table4, it can be seen that when the number of particles is increased, the performance of the proposed algorithm also increases. But when the number of particles exceeds some threshold (1000), the increase in the performance slows down.


## 7.1.  Single Maneuvering Target Tracking

The simulation study using nearly coordinate turn model is performed. The maneuvering target tracking is done by setting up a 2D flight path inx-y plane, which is similar to the path considered in [16]. At time step 1, the target starts at location [-310,310] in Cartesian coordinates in meters with initial velocity (in m/s) [10,-400]. The following trajectory is considered: a straight line with constant velocity between 1 and 17 s, a coordinated turn (0.09 rad/s) between 17 and 34 s, a straight line with constant velocity between 34 and 51 s, a coordinated turn (0.09 rad/s) between 51 and 68 s, and a straight line with constant velocity between 68 and 100 s.In the particle filter with process noise identification, a general model(7.1)Xk=ΦXk-1+Γvk-1
is adopted during the whole tracking process, where in (7.1),
(7.2)Φ=[Φb00Φb],Φb=[1ΔTΔT2201ΔT001],Γ=I6×6.
Matrix Φ is the transition matrix and ΔT is the sample interval. Xk=[px,γx,ax,py,γy,ay]kT is the state vector; px,γx, and ax denote, respectively, the position, velocity, and acceleration of the moving object along the x-axis of Cartesian frame and, py,γy, and ay along the y-axis. The equivalent process noise, vk-1=[vpx,vγx,vax,vpy,vγy,vay]k-1T, with unknown statistics is required to be identified. The bound of the process noise (d), which accounts for the uncertain dynamics, is chosen as {20m,20m/s,10m/s2,20m,20m/s,10m/s2}. The number of the process noise samples is equal to the number of particles, which is set as 500. The algorithm is initialized with Gaussian around the initial state of the true target, and the standard deviation of the Gaussian distribution is chosen as {10m,10m/s,5m/s2,10m,10m/s,5m/s2}.A track-while-scan (TWS) radar is positioned at the origin of the plane. The measurement equation is as follows:(7.3)Zk=h(Xk)+nk,
where Zk=[z1,z2]k is the observation vector. z1 is the distance between the radar and the target, and z2 is the bearing angle. The measurement noise nk=[nz1,nz2]k is a zero mean Gaussian white noise process with standard deviations of 20m(σz1) and 0.01rad(σz2). Resolution of the sensor is selected after from [29] (twice of the standard deviations of the measurement noise). The sampling interval is ΔT=1s.The particle filter based process noise identification algorithm is compared to the IMM filter and the regularized particle filter. The IMM filter consists of three extended Kalman filter (EKF) with different motion models. The details regarding these models may be found in [16]. The initial model probabilities and the mode switching probability matrix are set the same values as in [16]. For the regularised particle filter, Epanechnikov kernel is chosen as the rescaled kernel density, which is the same as in [30]. Moreover, the proposed algorithm is compared to its complete version in the same simulation setup.The simulation results are obtained from1000 Monte Carlo runs except the results from the complete version of the particle filter based multiscan JPDA algorithm, which is run for 100 times. Figure 1 shows the true trajectory of the maneuvering target. The root mean-square errors (RMSEs) in position at each time step, respectively, using the four methods are shown in Figure 2, where RPF and PFPNI represent the regularised particle filter algorithm and the particle filter based process noise identification algorithm. The performance of the methods is also compared via the global RMSE (in position), the tracking loss rate (TLR), and the executing time (ET), which are listed in Table 1. The tracking loss rate (TLR) is defined as the ratio of the number of simulations, in which the target is lost in track, to the total number of simulations carried out. The target is defined as lost in track when its global RMSE in position is larger than the ten times of the magnitude of the standard deviation in position. The executing time (ET) is the CPU time needed to execute one time step in MATLAB 7.1 on a 3 GHz (Mobile) Pentium IV operating under Windows 2000.Table 1
Performance comparison.
RMSE (m)
ET (s)
TLR
IMM
28.34
0.0239
3.7%
RPF
16.66
0.8708
0
PFPNI (simplified)
13.67
0.4503
0
PFPNI (complete)
13.29
190.47
0Figure 1
True trajectory of the single maneuvering target.Figure 2
RMSE in position using IMM, RPF, PFPNI (simplified), and PFPNI (complete) algorithms.From the simulation results, it can be seen that the simplified version of the proposed particle filter based process noise identification algorithm outperforms the IMM filter and the regularised particle filter algorithm, with computing time within the limits of practically realizable systems. Moreover, the proposed algorithm needs neither the possible multiple motion models nor the transition probability matrices, which makes it a more general algorithm for maneuvering target tracking. From the simulation results, it can also be seen that the complete version of the particle filter based process noise identification algorithm is not suitable for practical application due to the long computing time, though it gains a2.8% increase in accuracy (RMSE) compared with its simplified version.


## 7.2. Multiple Target Tracking

The simulation is carried out for tracking two slow-maneuvering targets in clutter. The Wiener process acceleration model (7.1) is chosen as the motion model for the two targets. The process noise, vk-1=[vpx,vγx,vax,vpy,vγy,vay]k-1T, is a zero mean Gaussian white noise process with standard deviations of 1m(σvpx), 1m/s(σvγx), 20m/s2(σvax), 1m(σvpy), 1m/s(σvγy), and 20m/s2(σvay).At time step 1, target one starts at location[-310,310] in x-y Cartesian coordinates in meters with the initial velocity (in m/s) [10,-400]. Target two starts at location (in m) [-310,-19000] with the initial velocity (in m/s) [10,400]. The length of each simulation run is 50 seconds. A TWS radar is positioned at the origin of the plane, whose details are provided in Section 7.1.The sampling interval isΔT=1s and it is assumed that the probability of detection PD=0.9 for the radar. For generating measurements in simulations, the clutter is assumed uniformly distributed with density 1×10-6/m2.In the particle filter based multiscan JPDA methods, each target model is assigned with 500 particles. The algorithm is initialized with Gaussians around the initial states of the true targets, and the standard deviations of the two Gaussian distributions are chosen equally as{10m,10m/s,5m/s2,10m,10m/s,5m/s2}.In the simulations carried out, the length of the multiple scan sliding window (L) varies from 1 to 3. The corresponding particle filter based multiscan JPDA methods with L=1,2,3 are utilized to track two slow-maneuvering targets individually. A comparison to the standard JPDA filter is also studied on the same simulation setup. The simulation results are obtained from 1000 Monte Carlo runs. Figure 3 shows the true trajectories of the two targets and Figure 4 shows the distance between the two targets along time, through which we can see that the two targets reach the smallest distance at time step 25. The RMSEs in position for the two targets are, respectively, shown in Figures 5 and 6, where PFJPDA represents the particle filter based JPDA filter with its following number in the bracket (e.g., (2.1)) denoting the scan number (L). The performance of the four methods is also compared in Table 2. The swap rate (SR) is defined as the ratio of the number of simulations, in which the two targets swap, to the total number of simulations.Table 2
Performance comparison.
RMSE (m)
ET (s)
TLR
SR
JPDA
T1: 20.24, T2: 22.87
0.07
24%
0
PFJPDA(1)
T1: 13.21, T2: 11.36
1.28
0
0
PFJPDA(2)
T1: 12.92, T2: 11.01
1.54
0
0
PFJPDA(3)
T1: 12.75, T2: 10.67
1.8
0
0Figure 3
True trajectories of the two slow-maneuvering targets.Figure 4
Distance between the two targets (- distance;-* zoom in distance).Figure 5
Target 1: RMSE in position using standard JPDA filter, particle filter based multiscan JPDA algorithms with a scan number equals to 1, 2, and 3, respectively.Figure 6
Target 2: RMSE in position using standard JPDA filter, particle filter based multiscan JPDA algorithms with a scan number equals to 1, 2, and 3, respectively.Compared with the standard JPDA method (based on extended Kalman filter), the particle filter based JPDA methods (single scan and multiple scan) are much more accurate and robust, at the cost of longer computing time. This verifies that when dealing with nonlinear problem (nonlinear observation equation) and large random acceleration (large process noise), the performance of particle filter is better than extended Kalman filter using local linearization.In the comparison between the three particle filter based JPDA methods (L=1,2,3, resp.), from Figures 5 and 6 it can be seen that there is no significant deference between the three methods except around time step 25, when the two targets are very close to each other (Figure 4). As the scan number increases, the RMSE (in position) of the corresponding algorithm decreases significantly around time step 25 (Figures 5 and 6). This verified that the particle filter based multiscan JPDA method provides a better performance especially when the targets are very close. The additional information of more scans improve the association probabilities in such critical situation, resulting in lower estimation errors (RMSE) and larger robustness (tracking loss rate).


## 7.3. Multiple Maneuvering Target Tracking

We now consider tracking two highly maneuvering targets in clutter. The true trajectories are shown in Figure7 in x-y plane. The velocities of the targets are shown in Figure 8, where VX and VY represent the velocity, respectively, in x and y coordinate. The distance between the two targets is shown in Figure 9. Target one is with the same setup as in the example proposed in Section 7.1. Target two starts at location [-310,-20310] in Cartesian coordinates in meters with the initial velocity (in m/s) [10,400]. Its trajectory is a straight line with constant velocity between 1 and 17 s, a coordinated turn (0.09rad/s) between 17 and 34 s, and a straight line constant velocity between 34 and 100 s.Figure 7
True trajectories of the two maneuvering targets.Figure 8
True velocities of the two maneuvering targets.Figure 9
Distance between the two targets (- distance;-* zoom in distance).A TWS radar is positioned at the origin of the plane (refer to Section7.1). The sampling interval (ΔT), the probability of detection (PD), and the clutter density are the same as in Section 7.2.In the proposed method, each target model is assigned with 1000 particles. The length of the multiple scan sliding window (L) is chosen equal to 3. The bound of the process noise (d) is the same as in Section 7.1. The algorithm is initialized with Gaussians around the initial states of the true targets, and the standard deviations of the two Gaussian distributions are chosen equally as {10m,10m/s,5m/s2,10m,10m/s,5m/s2}.In this work, the proposed method is compared to an IMMJPDA method and a multiscan IMMJPDA method, whose details could be found in [16, 20], respectively.The simulation results are obtained from1000 Monte Carlo runs. The RMSEs (resp. in position and velocity) at each time step for the two targets are shown in Figures 10, 11, 12, and 13, where MS-IMMJPDA and PFPNI-PFMSJPDA represent, respectively, the multiscan IMMJPDA algorithm and the proposed algorithm. The performance comparison is listed in Table 3. Moreover, for the proposed algorithm, the influence of the particle number in its performance is studied and simulations are carried out based on different sample sizes. The results are listed in Table 4.Table 3
Performance comparison.
RMSE in position (m)
RMSE in velocity (m/s)
ET (s)
TLR
SR
IMMJPDA
T1: 26.04, T2: 31.28
T1: 39.48, T2: 31.98
0.75
8.7%
13.9%
MS-IMMJPDA
T1: 25.30, T2: 30.26
T1: 39.22, T2: 31.76
3.78
0
4.5%
PFPNI-PFMSJPDA
T1: 14.42, T2: 18.14
T1: 34.44, T2: 30.58
4.0
0
3.7%Table 4
Influence of particle number in the performance of the proposed algorithm for tracking multiple maneuvering target.
Particle number
RMSE in position (m)
ET (s)
TLR
200
NA
NA
100%
500
T1: 17.37, T2: 21.22
2.9
35.9%
1000
T1: 14.42, T2: 18.14
4.0
0
2000
T1: 13.75, T2: 17.93
20.2
0Figure 10
Target 1: RMSE in position using IMMJPDA filter, multiple scan IMMJPDA filter, and the proposed method.Figure 11
Target 2: RMSE in position using IMMJPDA filter, multiple scan IMMJPDA filter, and the proposed method.Figure 12
Target 1: RMSE in velocity using IMMJPDA filter, multiple scan IMMJPDA filter, and the proposed method.Figure 13
Target 2: RMSE in velocity using IMMJPDA filter, multiple scan IMMJPDA filter, and the proposed method.The simulation results show that the proposed algorithm is more accurate than the multiple scan IMMJPDA filter and IMMJPDA filter, though it takes longer computing time. However, the three algorithms are implemented using Matlab. The computing time is considerably reduced when coded in C++. For the proposed method, it takes about 0.8 second for one time step in C++.In the scene of simulation, at the period of46s~52s, the two targets are very near to each other, which easily results in track swap. In the IMMJPDA method, the judgement of which measurements belong to which target depends on the information from the current scan. If the measurements from two targets are very close, it is hard to distinguish the different targets based on the measurements from a single scan, which leads to a track swap. In the proposed method and the multiple scan JPDA method, the information from several previous scans is combined with the information from current scan to calculate the association probabilities. The decision of which measurements belong to which target based on the information from multiple scans will be more accurate than single scan method, which reduces the swap rate effectively.From Table4, it can be seen that when the number of particles is increased, the performance of the proposed algorithm also increases. But when the number of particles exceeds some threshold (1000), the increase in the performance slows down.


## 8. Conclusions

A new algorithm is proposed for multiple maneuvering target tracking in the particle filter framework. In order to track the highly maneuvering target, the particle filter based process noise identification method is proposed to estimate the equivalent process noise induced by both maneuvering and random acceleration. Compared with the multiple model based methods for maneuvering target tracking, only one general model is adopted in the proposed method. In order to tackle the data association problem in multiple maneuvering target tracking, the particle filter based multiscan JPDA filter is adopted. Compared with the single scan JPDA method, the multiscan JPDA method uses richer information, which results in better estimated probabilities.



---

*Source: 372161-2012-12-24.xml*