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

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.

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. Section 2 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: 𝑥_𝑘=𝑓𝑥_𝑘−1,𝑣^∗_𝑘−1, (2.1) where 𝑥 is the state, 𝑓 is a possibly nonlinear function, and 𝑣^∗ is the known process noise with a zero mean Gaussian distribution. The objective of tracking is to recursively estimate 𝑥_𝑘 from a sequence of measurements up to time step 𝑘, 𝑧_1∶𝑘={𝑧₁,𝑧₂,…,𝑧_𝑘}. The observation model is described as follows: 𝑧_𝑘=ℎ𝑥_𝑘,𝑛_𝑘, (2.2) where ℎ is a possibly nonlinear function. 𝑛 is the observation noise with a zero mean Gaussian distribution. From the Bayesian perspective, the tracking problem is to recursively calculate the posterior distribution 𝑝(𝑥_𝑘∣𝑧_1∶𝑘).

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 distribution 𝑝(𝑥_𝑘∣𝑧_1∶𝑘) 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 {𝑥^𝑖₀,𝑖=1,…,𝐻} from the initial posterior distribution 𝑝(𝑥₀) and set the weights 𝑤^𝑖₀ to 1/𝐻, 𝑖=1,…,𝐻. 𝐻 is the number of particles.

(2)  Prediction. Particles at time step 𝑘−1, {𝑥^𝑖_𝑘−1,𝑖=1,…,𝐻}, are passed through the system model (2.3) to obtain the predicted particles at time step 𝑘, {𝑥^𝑖_𝑘,𝑖=1,…,𝐻}: 𝑥^𝑖_𝑘=𝑓𝑥^𝑖_𝑘−1,𝑣^∗,𝑖_𝑘−1, (2.3) where 𝑣^∗,𝑖_𝑘−1 is a sample drawn from the known zero mean Gaussian distribution.

(3)  Update. Once the observation data, 𝑧_𝑘, is measured, evaluate the importance weight of each predicted particle and obtain the normalized weight for each particle: 𝑤^𝑖_𝑘=𝑝𝑧_𝑘∣𝑥^𝑖_𝑘,𝑤^𝑖_𝑘=𝑤^𝑖_𝑘∑^𝐻_𝑖=1𝑤^𝑖_𝑘. (2.4) Thus at time step 𝑘 we can obtain the estimate of the state, 𝑥_𝑘=∑^𝐻_𝑖=1𝑤^𝑖_𝑘𝑥^𝑖_𝑘.

(4)  Resample. Resample the discrete distribution {𝑤^𝑖_𝑘∶𝑖=1,…,𝐻}, 𝐻 times to generate particles {𝑥^𝑗_𝑘∶𝑗=1,…,𝐻}, so that for any 𝑗, Pr{𝑥^𝑗_𝑘=𝑥^𝑖_𝑘}=𝑤^𝑖_𝑘. Set the weights 𝑤^𝑖_𝑘 to 1/𝐻,  𝑖=1,…,𝐻, and move to Stage 2.

3. Multiple Maneuvering Target Tracking Problem Formulation

The number of maneuvering targets (𝑀) 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 𝑥_𝑚,𝑘, where 𝑚 denotes the 𝑚th target and 𝑘 denotes time step 𝑘. The combined state, 𝑥_𝑘=(𝑥_1,𝑘,…,𝑥_𝑀,𝑘), is the concatenation of the individual target states. The individual targets are assumed to evolve independently according to a general motion model, 𝑥_𝑘=𝑓𝑥_𝑘−1,𝑢_𝑘−1,𝑣^∗_𝑘−1, (3.1) where 𝑢 is the maneuver acceleration and 𝑣^∗ is the process noise.

The observation vector 𝑧_𝑘 is composed of multiple sensor measurements {𝑧_𝑗,𝑘,𝑗=1,…,𝑁_𝑘}, where 𝑁_𝑘 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 Sections 4 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 𝑥_𝑘=𝑓𝑥_𝑘−1,𝑣_𝑘−1, (4.1) with an adequate accuracy, where 𝑣 is the equivalent noise that quantifies the error of this model in describing the target motions, in particular, maneuvers. The statistics of this noise 𝑣, 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 vector 𝑣_𝑘−1 is chosen as the state of the noise system. The observation vector is 𝑧_𝑘, which is the same as in (2.2). The observation equation is defined in 𝑧_𝑘=ℎ𝑥_𝑘,𝑛_𝑘=ℎ𝑓𝑥_𝑘−1,𝑣_𝑘−1,𝑛_𝑘. (4.2) The aim of the proposed method is to estimate the posterior distribution of the process noise, 𝑝(𝑣_𝑘−1∣𝑧_1∶𝑘). According to the Bayesian inference theory, 𝑝𝑣_𝑘−1∣𝑧_1∶𝑘=𝑝𝑧_𝑘∣𝑣_𝑘−1𝑝𝑣_𝑘−1∣𝑧_1∶𝑘−1𝑝𝑧_𝑘∣𝑧_1∶𝑘−1, (4.3) where 𝑝(𝑧_𝑘∣𝑧_1∶𝑘−1) is a normalizing constant and 1/𝑝(𝑧_𝑘∣𝑧_1∶𝑘−1) is defined as Υ. So we can obtain 𝑝𝑣_𝑘−1∣𝑧_1∶𝑘=Υ⋅𝑝𝑧_𝑘∣𝑣_𝑘−1𝑝𝑣_𝑘−1∣𝑧_1∶𝑘−1. (4.4) The noise 𝑣_𝑘−1 may be from random accelerations, sudden maneuvers, or both, and there is no information about the distribution of 𝑣_𝑘−1. 𝑣_𝑘−1 is not dependent on the previous measurements 𝑧_1∶𝑘−1, which results in 𝑝𝑣_𝑘−1∣𝑧_1∶𝑘=Υ⋅𝑝𝑧_𝑘∣𝑣_𝑘−1𝑝𝑣_𝑘−1. (4.5) 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 about 𝑣_𝑘−1, it is reasonable to assume that 𝑣_𝑘−1 is uniformly distributed in the range of 𝑈(−𝑑,𝑑), where 𝑈 denotes a uniform distribution and 𝑑 is the known process noise bound accounting for the maximin uncertain dynamics. According to the Monte Carlo theory, 𝑝(𝑣_𝑘−1) can be represented by 𝐻 samples, {𝑣^𝑗_𝑘−1,𝑗=1,…,𝐻}, from 𝑈(−𝑑,𝑑). Consider𝑝𝑣_𝑘−1=1𝐻^𝐻_𝑗=1𝛿𝑣_𝑘−1−𝑣^𝑗_𝑘−1. (4.6) The number of process noise samples (𝐻) is proportional to the magnitude of the maximum uncertain dynamics (‖𝑑‖). The posterior distribution of 𝑣_𝑘−1 can be represented as 𝑝𝑣_𝑘−1∣𝑧_1∶𝑘=Υ𝐻^𝐻_𝑗=1𝑝𝑧_𝑘∣𝑣^𝑗_𝑘−1𝛿𝑣_𝑘−1−𝑣^𝑗_𝑘−1=Υ𝐻^𝐻_𝑗=1𝜉^𝑗_𝑘⋅𝛿𝑣_𝑘−1−𝑣^𝑗_𝑘−1, (4.7) where 𝜉^𝑗_𝑘=𝑝(𝑧_𝑘∣𝑣^𝑗_𝑘−1) is defined as the weight assigned to the 𝑗th process noise sample 𝑣^𝑗_𝑘−1. The process noise samples {𝑣^𝑗_𝑘−1,𝑗=1,…,𝐻} are then resampled according to {𝜉^𝑗_𝑘,𝑗=1,…,𝐻} based on the principle thatPr{𝑣^𝑖_𝑘−1=𝑣^𝑗_𝑘−1}=𝜉^𝑗_𝑘, where {𝑣^𝑖_𝑘−1,𝑖=1,…,𝐻} are the process noise samples obtained from resampling. The obtained resampled process noise samples are approximately distributed according to the posterior distribution 𝑝(𝑣_𝑘−1∣𝑧_1∶𝑘).

In order to calculate 𝜉^𝑗_𝑘, the likelihood function 𝑝(𝑧_𝑘∣𝑣^𝑗_𝑘−1) is expanded based on the resampled state particles at time step 𝑘−1, {𝑥^𝑖_𝑘−1,𝑖=1,⋯,𝐻}. Consider𝑝𝑧_𝑘∣𝑣^𝑗_𝑘−1=^𝐻_𝑖=1𝑝𝑧_𝑘∣𝑥^𝑖_𝑘−1,𝑣^𝑗_𝑘−1𝑝𝑥^𝑖_𝑘−1∣𝑣^𝑗_𝑘−1. (4.8) Since 𝑥^𝑖_𝑘−1 and 𝑣^𝑗_𝑘−1 are independent, 𝑝(𝑥^𝑖_𝑘−1∣𝑣^𝑗_𝑘−1)=𝑝(𝑥^𝑖_𝑘−1).

The resmapled particles at time step 𝑘−1, {𝑥^𝑖_𝑘−1,𝑖=1,…,𝐻}, should be assigned with the same and equal weights, 1/𝐻. We can then obtain 𝑝𝑥^𝑖_𝑘−1=1𝐻. (4.9)

To calculate 𝑝(𝑧_𝑘∣𝑥^𝑖_𝑘−1,𝑣^𝑗_𝑘−1), define 𝜇_𝑘^𝑖,𝑗 as the intermediate particle, 𝜇_𝑘^𝑖,𝑗=𝑓𝑥^𝑖_𝑘−1,𝑣^𝑗_𝑘−1,𝑖=1,…𝐻,𝑗=1,…,𝐻 (4.10) and expand 𝑝(𝑧_𝑘∣𝑥^𝑖_𝑘−1,𝑣^𝑗_𝑘−1) based on 𝜇_𝑘^𝑖,𝑗, 𝑝𝑧_𝑘∣𝑥^𝑖_𝑘−1,𝑣^𝑗_𝑘−1=^𝐻𝐻_𝑝=1𝑞=1𝑝𝑧_𝑘∣𝜇_𝑘^𝑝,𝑞,𝑥^𝑖_𝑘−1,𝑣^𝑗_𝑘−1𝑝𝜇_𝑘^𝑝,𝑞∣𝑥^𝑖_𝑘−1,𝑣^𝑗_𝑘−1. (4.11) Since 𝑥^𝑖_𝑘−1 and 𝑣^𝑗_𝑘−1 are known and 𝜇_𝑘^𝑝,𝑞 is obtained from a purely deterministic relationship in (4.10), we obtain 𝑝𝜇_𝑘^𝑝,𝑞∣𝑥^𝑖_𝑘−1,𝑣^𝑗_𝑘−1=𝑝𝑧1,for𝑝=𝑖and𝑞=𝑗,0,for𝑝≠𝑖or𝑞≠𝑗,_𝑘∣𝑥^𝑖_𝑘−1,𝑣^𝑗_𝑘−1=𝑝𝑧_𝑘∣𝜇_𝑘^𝑖,𝑗. (4.12)(4.13) Combining (4.9) and (4.13) with (4.8) results in 𝑝𝑧_𝑘∣𝑣^𝑗_𝑘−1=^𝐻_𝑖=1𝑝𝑧_𝑘∣𝜇_𝑘^𝑖,𝑗⋅1𝐻. (4.14) At each time step, the process noise samples are drawn from a uniform distribution 𝑈(−𝑑,𝑑). Each process noise sample 𝑣^𝑗_𝑘−1 is evaluated and assigned its corresponding weight 𝜉^𝑗_𝑘. 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 {𝑥^𝑖_𝑘∶𝑖=1,…,𝐻} are then obtained based on the resampled process noise samples {𝑣^𝑖_𝑘−1∶𝑖=1,…,𝐻} 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 𝑘−1, draw process noise samples {𝑣^𝑗_𝑘−1∶𝑗=1,…,𝐻} from a uniform distribution 𝑈(−𝑑,𝑑).
(2)  Calculate the intermediate particles {𝜇_𝑘^𝑖,𝑗∶𝑖=1,…,𝐻;𝑗=1,…,𝐻} according to (4.10).
(3)  Calculate the process noise sample weights {𝜉^𝑗_𝑘∶𝑗=1,…,𝐻} via (4.14).
(4)  Resample the discrete distribution {𝜉^𝑗_𝑘∶𝑗=1,…,𝐻}, 𝐻 times to generate the resampled process noise samples {𝑣^𝑖_𝑘−1∶𝑖=1,…,𝐻}, so that for any 𝑖, Pr{𝑣^𝑖_𝑘−1=𝑣^𝑗_𝑘−1}=𝜉^𝑗_𝑘. Set the weights 𝜉^𝑗_𝑘 to 1/𝐻,𝑖=1,…,𝐻.
(5)  Calculate the predicted particles at time step 𝑘, {𝑥^𝑖_𝑘,𝑖=1,…,𝐻}, 𝑥^𝑖_𝑘=𝑓𝑥^𝑖_𝑘−1,𝑣^𝑖_𝑘−1. (4.15)(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, 𝐻∗𝐻 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 𝐻 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 {𝑥^𝑖_𝑘−1,𝑖=1,…,𝐻} are replaced by 𝑥_𝑘−1, the estimate of the state at time step 𝑘−1, which results in a simplified version of (4.10),𝜇^𝑗_𝑘=𝑓𝑥_𝑘−1,𝑣^𝑗_𝑘−1, (4.16) and 𝑝(𝑧_𝑘∣𝑣^𝑗_𝑘−1) is expanded directly on 𝜇^𝑗_𝑘, 𝑝𝑧_𝑘∣𝑣^𝑗_𝑘−1=^𝐻_𝜏=1𝑝𝑧_𝑘∣𝜇^𝜏_𝑘𝑝𝜇^𝜏_𝑘∣𝑣^𝑗_𝑘−1. (4.17) Using the similar derivation process with (4.12), we can obtain 𝑝𝑧_𝑘∣𝑣^𝑗_𝑘−1=𝑝𝑧_𝑘∣𝜇^𝑗_𝑘. (4.18)
In the simplified version of the proposed algorithm, the number of intermediate particles is reduced to 𝐻, 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 bound 𝑑 preliminarily, which could be estimated based on the maximin uncertain dynamics. For most of the manned air vehicles, the process noise bound 𝑑 could be determined according to the flight envelop of the aircraft. For the unmanned air vehicles with unexpected maneuvers, 𝑑 is set large enough to cover the unexpected maneuvers. If our knowledge about 𝑑 is imprecise, tracking process may fail when coping with the maneuvers that are not covered by the process noise bound 𝑑.

5. Multiple Target Tracking Using Particle Filter Based Multiscan JPDAF

The number of targets (𝑀) 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 𝑥_𝑚,𝑘, where 𝑚 denotes the 𝑚th target and 𝑘 denotes time step 𝑘. The combined state, 𝑥_𝑘=(𝑥_1,𝑘,…,𝑥_𝑀,𝑘), is the concatenation of the individual target states. The individual targets are assumed to evolve independently according to the Markovian dynamic models 𝑝_𝑚(𝑥_𝑚,𝑘∣𝑥_{𝑚,𝑘−1}). The observation vector 𝑧_𝑘 is composed of multiple sensor measurements {𝑧_𝑗,𝑘,𝑗=1,…,𝑁_𝑘}, where 𝑁_𝑘 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 targets 𝑝_𝑚(𝑥_𝑚,𝑘∣𝑧_1∶𝑘),𝑚=1⋯𝑀, and these distributions are approximated with particles, {𝑥^𝑖_𝑚,𝑘,𝑖=1⋯𝐻,𝑚=1⋯𝑀}, and their associated weights {𝑤^𝑖_𝑚,𝑘,𝑖=1⋯𝐻,𝑚=1⋯𝑀}, as in 𝑝_𝑚𝑥_𝑚,𝑘∣𝑧_1∶𝑘=^𝐻_𝑖=1𝑤^𝑖_𝑚,𝑘𝛿𝑥_𝑚,𝑘−𝑥^𝑖_𝑚,𝑘. (5.1) 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 𝜆_{𝑘−𝐿+1∶𝑘}, where 𝐿 denotes the length of the multiple scan sliding window. The multiple scan joint association events are mutually exclusive, and they form a complete set Λ_{𝑘−𝐿+1∶𝑘} [11]. 𝜆_{𝑘−𝐿+1∶𝑘} is composed by the association vectors at each scan in the sliding window, 𝜆_{𝑘−𝐿+1∶𝑘}=(𝜃_{𝑘−𝐿+1},𝜃_{𝑘−𝐿+2},…,𝜃_𝑘). The elements of the association vector at time step 𝑘, 𝜃_𝑘=(𝜁_1,𝑘,…,𝜁_𝑗,𝑘,…,𝜁_{𝑁𝑘,𝑘}) are given by 𝜁_𝑗,𝑘=0,if𝑧_𝑗,𝑘{}isduetoclutter,𝑚∈1⋯𝑀,if𝑧_𝑗,𝑘isfromtarget𝑚. (5.2) The heart of the new algorithm is to find the posterior probability for the joint association event of multiple scans. That is to calculate 𝑝(𝜆_{𝑘−𝐿+1∶𝑘}∣𝑧_1∶𝑘) and it can be written as 𝑝𝜆_{𝑘−𝐿+1∶𝑘}∣𝑧_1∶𝑘∝𝑝𝑧_𝑘⋯𝑧_{𝑘−𝐿+1}∣𝜆_{𝑘−𝐿+1∶𝑘},𝑧_{1∶𝑘−𝐿}𝑝𝜆_{𝑘−𝐿+1∶𝑘}∣𝑧_{1∶𝑘−𝐿}∝𝑝𝑧_𝑘⋯𝑧_{𝑘−𝐿+1}∣𝜆_{𝑘−𝐿+1∶𝑘},𝑧_{1∶𝑘−𝐿}𝑝𝜆_{𝑘−𝐿+1∶𝑘}, (5.3) where the conditioning of 𝜆_{𝑘−𝐿+1∶𝑘} 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 𝑝𝑧_𝑘⋯𝑧_{𝑘−𝐿+1}∣𝜆_{𝑘−𝐿+1∶𝑘},𝑧_{1∶𝑘−𝐿}=^𝐿_𝑠=1^{𝑁𝑘−𝐿+𝑠}_𝑗=1𝑝𝑧_{𝑗,𝑘−𝐿+𝑠}∣𝜆_{𝑘−𝐿+1∶𝑘},𝑧_{1∶𝑘−𝐿}. (5.4) To reduce the notation, the index of the scan 𝑠 in the sliding window is denoted by 𝑘_𝑠=𝑘−𝐿+𝑠. Then we can obtain 𝑝𝑧_𝑘⋯𝑧_{𝑘−𝐿+1}∣𝜆_{𝑘−𝐿+1∶𝑘},𝑧_{1∶𝑘−𝐿}=^𝐿_𝑠=1^𝑁𝑘𝑠_𝑗=1𝑝𝑧_{𝑗,𝑘𝑠}∣𝜆_{𝑘−𝐿+1∶𝑘},𝑧_{1∶𝑘−𝐿}=^𝐿_{𝑠=1𝑗∈𝐼𝑠0,𝑘}𝑝_𝑐𝑧_{𝑗,𝑘𝑠}⋅_{𝑗∈𝐼𝑘𝑠}𝑝𝑧_{𝑗,𝑘𝑠}∣𝑥_{𝜁𝑠𝑗,𝑘,𝑘𝑠}=^𝐿_𝑠=1(𝑉)^{−𝐶𝑘𝑠}⋅_{𝑗∈𝐼𝑘𝑠}𝑝𝑧_{𝑗,𝑘𝑠}∣𝑥_{𝜁𝑠𝑗,𝑘,𝑘𝑠}, (5.5) where 𝐼_0,𝑘𝑠={𝑗∈{1,…,𝑁_𝑘𝑠}∶𝜁_{𝑗,𝑘𝑠}=0} and 𝐼_𝑘𝑠={𝑗∈{1,…,𝑁_𝑘𝑠}∶𝜁_{𝑗,𝑘𝑠}≠0} are, respectively, the subsets of measurement indices corresponding to clutter measurements and measurements from the targets being tracked, on scan 𝑘_𝑠. 𝑝_𝑐 denotes the clutter likelihood model, which is assumed to be uniform over the volume of the surveillance area 𝑉. The volume of the surveillance area could be calculated as per 𝑉=2𝜋𝑅_max, where 𝑅_max is the maximum range of the sensor. 𝐶_𝑘𝑠 is defined as the number of clutter measurements.

The joint association prior 𝑝(𝜆_{𝑘−𝐿+1∶𝑘}), can be calculated as in (5.6) according to [4, 16], 𝑝𝜆_{𝑘−𝐿+1∶𝑘}=^𝐿_𝑠=1𝐶_𝑘𝑠!𝜀𝑁_𝑘𝑠!^𝑀_𝑚=1𝑃_𝐷^{𝛿𝑚(𝜃𝑘𝑠)}1−𝑃_𝐷^{1−𝛿𝑚(𝜃𝑘𝑠)}, (5.6) where 𝜀 is a “diffuse” prior [16] and 𝑃_𝐷 is the detection probability. 𝛿_𝑚(𝜃_𝑘𝑠) is a binary variable and set to one if the 𝑚th target is assigned with a measurement in the event 𝜃_𝑘𝑠.

The posterior probability for the joint association event of multiple scans is obtained as 𝑝𝜆_{𝑘−𝐿+1∶𝑘}∣𝑧_1∶𝑘∝𝑝𝜆_{𝑘−𝐿+1∶𝑘}^𝐿_𝑠=1(𝑉)^{−𝐶𝑘𝑠}⋅_{𝑗∈𝐼𝑘𝑠}𝑝𝑧_{𝑗,𝑘𝑠}∣𝑥_{𝜁𝑠𝑗,𝑘,𝑘𝑠}. (5.7) The posterior probability that the 𝑗th measurement is associated with the 𝑚th target on scan 𝑘, 𝛽_𝑗𝑚, is calculated by summing over the probabilities of the corresponding joint association events via 𝛽_𝑗𝑚=𝑝𝜁_𝑗,𝑘=𝑚∣𝑧_1∶𝑘=_{𝜆𝑘−𝐿+1∶𝑘∈Λ𝑘−𝐿+1∶𝑘∶𝜁𝑗,𝑘=𝑚}𝑝𝜆_{𝑘−𝐿+1∶𝑘}∣𝑧_1∶𝑘. (5.8)

These approximations can, in turn, be used in (5.9) to approximate the target likelihood according to [24], 𝑝_𝑚𝑧_𝑘∣𝑥_𝑚,𝑘=𝛽_0𝑚+^𝑁𝑘_𝑗=1𝛽_𝑗𝑚𝑝𝑧_𝑗,𝑘∣𝑥_𝑚,𝑘, (5.9) where 𝛽_0𝑚 is the posterior probability that the 𝑚th target is undetected. Finally, setting the new importance weights to 𝑤^𝑖_𝑚,𝑘∝𝑤^𝑖_{𝑚,𝑘−1}𝑝_𝑚𝑧_𝑘∣𝑥^𝑖_𝑚,𝑘𝑝_𝑚𝑥^𝑖_𝑚,𝑘∣𝑥^𝑖_{𝑚,𝑘−1}𝑞_𝑚𝑥^𝑖_𝑚,𝑘∣𝑥^𝑖_{𝑚,𝑘−1𝐻𝑖=1},𝑤^𝑖_𝑚,𝑘=1, (5.10) where 𝑞_𝑚(𝑥^𝑖_𝑚,𝑘∣𝑥^𝑖_{𝑚,𝑘−1}) is the proposal distribution, which is used to generate the predicted particles. Equation (5.10) leads to the sample set {𝑤^𝑖_𝑚,𝑘,𝑥^𝑖_𝑚,𝑘}^𝐻_𝑖=1 being approximately distributed according to the marginal filtering distribution 𝑝_𝑚(𝑥_𝑚,𝑘∣𝑧_1∶𝑘).

A summary of the particle filter based multiple scan JPDA filter algorithm is presented in what follows. Assuming that the sample sets {𝑤^𝑖_{𝑚,𝑘−1},𝑥^𝑖_{𝑚,𝑘−1}}^𝐻_𝑖=1,𝑚=1⋯𝑀, are approximately distributed according to the corresponding marginal filtering distributions at the previous time step 𝑝_𝑚(𝑥_{𝑚,𝑘−1}∣𝑧_1∶𝑘−1),𝑚=1⋯𝑀, the algorithm proceeds as follows at the current time step.

Algorithm 5.1 (particle filter based multiscan JPDA filter). (1)  For 𝑚=1⋯𝑀,𝑖=1⋯𝐻, generate predicted particles for the target states 𝑥^𝑖_𝑚,𝑘∼𝑞_𝑚(𝑥^𝑖_𝑚,𝑘∣𝑥^𝑖_{𝑚,𝑘−1}).
(2)  For 𝑚=1⋯𝑀, calculate 𝑥_𝑚,𝑘, the preapproximation of 𝑥_𝑚,𝑘, which is to be substituted into (5.7) to calculate the posterior probability of the joint association event in multiple scan, 𝑝(𝜆_{𝑘−𝐿+1∶𝑘}∣𝑧_1∶𝑘). Consider𝑥_𝑚,𝑘=^𝐻_𝑖=1𝑤^𝑖_{𝑚,𝑘−1}𝑥^𝑖_𝑚,𝑘. (5.11)
(3)  Enumerate all the valid joint measurement to target association events in the sliding window 𝑘−𝐿+1∶𝑘 to form the set Λ_{𝑘−𝐿+1∶𝑘}.
(4)  For each 𝜆_{𝑘−𝐿+1∶𝑘}∈Λ_{𝑘−𝐿+1∶𝑘}, compute the posterior probability of the joint association event in multiple scan, 𝑝(𝜆_{𝑘−𝐿+1∶𝑘}∣𝑧_1∶𝑘), via (5.7).
(5)  For 𝑚=1⋯𝑀,𝑗=1⋯𝑁_𝑘, compute the marginal association posterior probability, 𝛽_𝑗𝑚, via (5.8).
(6)  For 𝑚=1⋯𝑀,𝑖=1⋯𝐻, compute the target likelihood, 𝑝_𝑚(𝑧_𝑘∣𝑥^𝑖_𝑚,𝑘), via (5.9).
(7)  For 𝑚=1⋯𝑀,𝑖=1⋯𝐻, compute and normalize the particle weights, via (5.10).
(8)  Resample the discrete distribution {𝑤^𝑖_𝑚,𝑘∶𝑖=1,…,𝐻}, 𝐻 times to generate particles {𝑥^𝑗_𝑚,𝑘∶𝑗=1,…,𝐻}, so that for any 𝑗, Pr{𝑥^𝑗_𝑚,𝑘=𝑥^𝑖_𝑚,𝑘}=𝑤^𝑖_𝑚,𝑘. Set the weights 𝑤^𝑖_𝑚,𝑘 to 1/𝐻,𝑖=1,…,𝐻, 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 Section 4) 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 𝑘−1, for 𝑚=1⋯𝑀,(a)draw process noise samples {𝑣^𝑗_{𝑚,𝑘−1}∶𝑗=1,…,𝐻} from a uniform distribution 𝑈(−𝑑,𝑑),(b)calculate the intermediate particles {𝜇^𝑗_𝑚,𝑘∶𝑗=1,…,𝐻} according to 𝜇^𝑗_𝑚,𝑘=𝑓𝑥_{𝑚,𝑘−1},𝑣^𝑗_{𝑚,𝑘−1}, (6.1)(c)calculate the process noise sample weights {𝜉^𝑗_𝑚,𝑘∶𝑗=1,…,𝐻} via (6.2) and normalize each weight: 𝜉^𝑗_𝑚,𝑘=𝑝𝑧^𝑚_𝑘∣𝜇^𝑗_𝑚,𝑘, (6.2) where 𝑧^𝑚_𝑘 are the measurements that are close to 𝜇^𝑗_𝑚,𝑘 obtained from the nearest neighbor method,(d)resample the discrete distribution {𝜉^𝑗_𝑚,𝑘∶𝑗=1,…,𝐻}, 𝐻 times to generate the new process noise samples {𝑣^𝑖_{𝑚,𝑘−1}∶𝑖=1,…,𝐻}, so that for any 𝑖, Pr{𝑣^𝑖_{𝑚,𝑘−1}=𝑣^𝑗_{𝑚,𝑘−1}}=𝜉^𝑗_𝑚,𝑘. Set the weights 𝜉^𝑗_𝑚,𝑘 to 1/𝐻,𝑖=1,…,𝐻,(e)obtain the predicted particles {𝑥^𝑖_𝑚,𝑘∶𝑖=1,…,𝐻} at time step 𝑘 from the new process noise samples via 𝑥^𝑖_𝑚,𝑘=𝑓𝑥^𝑖_{𝑚,𝑘−1},𝑣^𝑖_{𝑚,𝑘−1}. (6.3)
(2)  At time step 𝑘, for 𝑚=1⋯𝑀, calculate 𝑥_𝑚,𝑘, the postapproximation of 𝑥_𝑚,𝑘, 𝑥_𝑚,𝑘=^𝐻_𝑖=1𝜉^𝑖_𝑚,𝑘𝑥^𝑖_𝑚,𝑘, (6.4) where {𝜉^𝑖_𝑚,𝑘∶𝑖=1,…,𝐻} are the process noise samples' weights and are calculated based on the measurements at the current time step (𝑘). However, in (5.11) 𝑥_𝑚,𝑘 is calculated based on {𝑤^𝑖_{𝑚,𝑘−1}∶𝑖=1,…,𝐻}, which rely on the measurements at the previous time step (𝑘−1). As a result, in maneuvering target tracking, 𝑥_𝑚,𝑘 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 𝑥_𝑚,𝑘=∑^𝐻_𝑖=1𝑤^𝑖_𝑚,𝑘𝑥^𝑖_𝑚,𝑘.  𝑥_𝑚,𝑘 is the estimate of the true state 𝑥_𝑚,𝑘.

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 in 𝑥-𝑦 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 𝑋_𝑘=Φ𝑋_𝑘−1+Γ𝑣_𝑘−1 (7.1) is adopted during the whole tracking process, where in (7.1), ΦΦ=^𝑏00Φ_𝑏,Φ_𝑏=1Δ𝑇Δ𝑇²201Δ𝑇,001Γ=𝐼_6×6. (7.2) Matrix Φ is the transition matrix and Δ𝑇 is the sample interval. 𝑋_𝑘=[𝑝_𝑥,𝛾_𝑥,𝑎_𝑥,𝑝_𝑦,𝛾_𝑦,𝑎_𝑦]^𝑇_𝑘 is the state vector; 𝑝_𝑥,𝛾_𝑥, and 𝑎_𝑥 denote, respectively, the position, velocity, and acceleration of the moving object along the 𝑥-axis of Cartesian frame and, 𝑝_𝑦,𝛾_𝑦, and 𝑎_𝑦 along the 𝑦-axis. The equivalent process noise, 𝑣_𝑘−1=[𝑣_𝑝𝑥,𝑣_𝛾𝑥,𝑣_𝑎𝑥,𝑣_𝑝𝑦,𝑣_𝛾𝑦,𝑣_𝑎𝑦]^𝑇_𝑘−1, with unknown statistics is required to be identified. The bound of the process noise (𝑑), which accounts for the uncertain dynamics, is chosen as {20m,20m/s,10m/s²,20m,20m/s,10m/s²}. 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/s²,10m,10m/s,5m/s²}.

A track-while-scan (TWS) radar is positioned at the origin of the plane. The measurement equation is as follows: 𝑍_𝑘=ℎ𝑋_𝑘+𝑛_𝑘, (7.3) where 𝑍_𝑘=[𝑧₁,𝑧₂]_𝑘 is the observation vector. 𝑧₁ is the distance between the radar and the target, and 𝑧₂ is the bearing angle. The measurement noise 𝑛_𝑘=[𝑛_𝑧1,𝑛_𝑧2]_𝑘 is a zero mean Gaussian white noise process with standard deviations of 20m(𝜎_𝑧1) and 0.01rad(𝜎_𝑧2). Resolution of the sensor is selected after from [29] (twice of the standard deviations of the measurement noise). The sampling interval is Δ𝑇=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 from 1000 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.

Figure 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 a 2.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, 𝑣_𝑘−1=[𝑣_𝑝𝑥,𝑣_𝛾𝑥,𝑣_𝑎𝑥,𝑣_𝑝𝑦,𝑣_𝛾𝑦,𝑣_𝑎𝑦]^𝑇_𝑘−1, is a zero mean Gaussian white noise process with standard deviations of 1m(𝜎_𝑣𝑝𝑥), 1m/s(𝜎_𝑣𝛾𝑥), 20m/s²(𝜎_𝑣𝑎𝑥), 1m(𝜎_𝑣𝑝𝑦), 1m/s(𝜎_𝑣𝛾𝑦), and 20m/s²(𝜎_𝑣𝑎𝑦).

At time step 1, target one starts at location [−310,310] in 𝑥-𝑦 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 Δ𝑇=1s and it is assumed that the probability of detection 𝑃_𝐷=0.9 for the radar. For generating measurements in simulations, the clutter is assumed uniformly distributed with density 1×10⁻⁶/m².

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/s²,10m,10m/s,5m/s²}.

In the simulations carried out, the length of the multiple scan sliding window (𝐿) varies from 1 to 3. The corresponding particle filter based multiscan JPDA methods with 𝐿=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 (𝐿). 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.

Figure 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 (𝐿=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 Figure 7 in 𝑥-𝑦 plane. The velocities of the targets are shown in Figure 8, where 𝑉𝑋 and 𝑉𝑌 represent the velocity, respectively, in 𝑥 and 𝑦 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 Section 7.1). The sampling interval (Δ𝑇), the probability of detection (𝑃_𝐷), 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 (𝐿) is chosen equal to 3. The bound of the process noise (𝑑) 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/s²,10m,10m/s,5m/s²}.

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 from 1000 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.

Table 4 

Influence of particle number in the performance of the proposed algorithm for tracking multiple maneuvering target.

Figure 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 of 46s~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 Table 4, 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.

Acknowledgments

This work was supported by the Foundation for Innovative Research Groups of the National Natural Science Foundation of China (60921003), the Natural Science Foundations of China (nos. 61104051, 61074176, and 61004087), and the Program for New Century Excellent Talents in University.