A Bayesian Approach to Track Multiple Extended Targets Using Particle Filter for Nonlinear System

To track multiple extended targets for the nonlinear system, this paper employs the idea of the particle filter to track kinematic states and shape formation of extended targets. First, the Bayesian framework is proposed for multiple extended targets to jointly estimate multiple extended target state and association hypothesis. Furthermore, a joint proposal distribution is defined for the multiple extended target state and association hypothesis. Then, the Bayesian framework of multiple extended target tracking is implemented by the particle filtering which could release the high computational burden caused by the increase in the number of extended targets and measurements. Simulation results show that the proposed multiple extended target particle filter has superior performance in shape estimation and improves the performance of the position estimation in the situation that there are spatially closed extended targets.


Introduction
In most multitarget tracking applications it is assumed that each target produces at most one measurement per time step.This is reasonable when their extension is assumed to be neglected in comparison with sensor resolution and error.However, with increasing sensor resolution capabilities this assumption is no longer valid, e.g., in short-range applications or for maritime surveillance where different scattering centers of the object under consideration may give rise to several distinct detection varying from scan to scan, both in its number and the relative origin location.
Extended target tracking has attracted much attention in the last decade.Gilholm et al. [1,2] present an approach for tracking the extended target under the assumption that the target extent is modeled by a spatial probability distribution and each target-related measurement is randomly and independently drawn from this spatial distribution.However, in real application the spatial probability is difficult to obtain.Baum and Hanebeck present the random hypersurface model [3], which describes the distribution of measurement sources, used to track elliptic targets [4] and more general shapes [5] as well.Another approach to elliptic target modeling is the random matrix framework in [6], where the measurement noise is assumed to be a zero-mean normally distributed random vector with variance depending on the current object extension and the Bayesian recursion is derived for the target extension and kinematical states of extended target.The target kinematical states are modeled using a Gaussian distribution, while the ellipsoidal target extension is modeled using an inverse Wishart distribution.In [7], Koch and Feldmann track group targets under kinematical constraints by making use of random matrices.Modifications and improvements to the Gaussian-inverse Wishart model of [6] are suggested in [8].Measurements of target downrange extent are used to aid track retention in [9].Zea et al. track elongated extended objects using splines in [10].
Those above-mentioned methods are to track single extended target.There are two kinds of approaches for tracking multiple extended targets.One is based on random finite sets.For example, in [11], Mahler presents an extension of the probability hypothesis density (PHD) filter [12] to handle extended targets.Orguner et al. propose a cardinalized probability hypothesis density (CPHD) filter for extended targets [13], and Lian et al. propose unified CPHD filters for extended targets and unresolved targets [14].Based on the Gaussian-mixture PHD filter [15], Granstrom presents the extended target GMPHD (ET-GMPHD) filter for extended target tracking in [16] and in [17] describes much more details and extensive investigations of the methodology.The ET-GMPHD filter still has the problem that its performance would obviously descend in the situation where two or more extended targets are spatially close [16,17].PHD implementations such as the sequential Monte Carlo PHD (SMCPHD) filter and the GMPHD filter do not include object identities.In some cases, we need to know the track continuity of objects in order for postprocessing, such as behavior of objects and activity recognition.There have been studies on the track continuity in PHD implementations [18][19][20][21].The Gaussian process measurement model is integrated in the recently developed labeled multi-Bernoulli filter for extended objects in [22].However, in theory this kind of methods needs to consider all possible partitions of the measurement set and thus is computation-intensive though some algorithms have been proposed to consider only a subset of all possible partitions [16,17].The other is those approaches which extend the data association algorithm applied in tracking point targets to track multiple extended targets.In [23,24], the random matrix model [6] and the random hypersurface model [3] are, respectively, integrated into a Probabilistic Multi-Hypothesis Tracking (PMHT) framework to track multiple extended targets.In [25] combined with the partitioning measurements algorithm, the JPDA (joint probabilistic data association) is applied to describe the association between measurement clusters and extended targets, and then the algorithm to track multiple extended targets is formed.However, those methods only consider the linear Gaussian dynamic system.
There are already single extended target tracking methods for nonlinear systems.In [26], a Rao-Blackwellised particle filter, which estimates the linear part of the state space with a Kalman filter and the nonlinear states using a particle filter, is proposed for single target tracking.Lan and Rong Li linearize the nonlinear measurement function to utilize the effectiveness and efficiency of the random matrix approach [27].To the best of our knowledge, almost all existing multiple extended target tracking methods are not specifically for nonlinear system, and to solve nonlinear system some of those methods apply approximations; for example, in [28] the ET-GMPHD filter applies UKF to handle nonlinear measurement model.But, the performance of those algorithms degrades rapidly as the nonlinearities become more severe.Sequential Monte Carlo, i.e., particle filtering, has become a practical numerical technique to approximate the Bayesian tracking recursions for nonlinear and/or non-Gaussian models, due to its efficiency, simplicity, flexibility, ease of implementation, and modeling success over a wide range of challenging applications [29].Vermaak et al. solve the data association problem within the context of particle filtering in tracking an extended target [30] where extended target is modeled as a set of fixed point features, and in multiple points target tracking [31].
To track multiple extended targets for the nonlinear system, this paper employs the idea of the particle filtering to jointly estimate the measurement-target data association and states of extended targets and then proposes a multiple extended target tracking algorithm.To do this, first, the Bayesian framework is proposed for multiple extended targets.Furthermore, the joint proposal distribution is defined for the multiple extended target combined state and association hypothesis, which could be obtained by the measurement model and the extended target dynamic evolution model.
The rest of the paper is organized as follows.We briefly describe the multiple extended target tracking problem in Section 2. Section 3 presents the proposed particle filter for multiple extended targets.The Bayesian framework of multiple extended targets and the joint proposal distribution are also discussed in this section.In Section 4, simulation results are given to demonstrate the performance of the proposed filter.Section 5 presents concluding remarks and outlines future research directions.

Model Description and Problem Formulation
Considering a nonlinear system, the state evolution of the extended target  could be modeled by where x   is the state of the extended target  at the current moment , f  is the target state evolution function, and ^  is the process noise.The measurement model of the extended target  could be expressed by where h  is the measurement generation function and   is the measurement noise.One extended target could generate multiple measurements which contain the shape information of the extended target, and thus the extended target tracking is usually to estimate the kinematic state and shape parameters.In this paper, the extended target state x k consists of the kinematic state x   and the shape parameters x   at the current moment ; i.e., x  = [(x   ) T , (x   ) T ] T .The combined state is defined as the composition of the individual extended target state and is denoted by T at the current moment.Assume that there are   extended targets over the surveillance region, and individual extended targets evolve independently according to Markovian dynamic models (1) [17,31]; the dynamics of the combined state can be decomposed into a product of the individual extended targets dynamics: where (x   | x  −1 ) is the probability density function of the extended target  state transition.
At the current moment, the measurement set   consists of the measurements from extended targets and clutter measurements, and each extended target can generate more than one measurement at each time step.To deal with the data association problem, we introduce the association hypothesis denoted by where   is the number of clutter measurements,   is the number of measurements from extended targets, and a  is the association vector.The association vector a  is given by where  is the number of measurements at the current moment (note that  =   +   );    belongs to [0   ].    = 0 means that the measurement  generates from clutter, and    =  ̸ = 0 means that the measurement  stems from the extended target .
The association hypothesis   is unknown in most practical applications.Nevertheless, a likelihood model for the measurements, which is conditional on a known association hypothesis, could be obtained, and then the marginal likelihood model for measurements could be derived through marginalising over the association hypothesis uncertainty.Assume that the measurement is generated independently [2], the likelihood model conditional on a known association hypothesis could be given by where I 0 = { ∈ [1, . . ., ] :    = 0} is the subset of clutter measurements, and I = { ∈ [1, . . ., ] :    =  ̸ = 0} is the subset of measurements from extended targets.  is the clutter likelihood model, which is normally assumed to be uniform over the surveillance region.(z   | x   ) is the likelihood of the measurement  associated with the target .
For a given number of measurements originating from extended targets   , the number of ways of choosing a subset of   measurements from  measurements is !/   !( −   )!.The number of possible associations between a measurement and   extended targets is   since each measurement may originate from any extended target, and thus that of between   measurements and   extended targets is     .Therefore, for a given number of measurements from extended targets   the total number of association hypotheses is shown as follows: However, the number of measurements from extended targets   is unknown, and thus the total number of association hypotheses is given by If the total number of association hypotheses   is not large, the marginal likelihood could be obtained through marginalising over the association uncertainty, as shown below where (  |    ) is the prior distribution of the association hypothesis, which is assumed to be of the form For a known number of measurements from extended targets   , the prior probability of each association vector is assumed to be the same and thus could be given by The number of clutter measurements is assumed to follow a Poisson distribution, with the expected number   : Since each of the extended targets can generate more than one measurement at each time step,   measurements could stem from one extended target or more than one.The number of measurements generated by the extended target  at each time step is assumed to be a Poisson distributed random variable with rate   measurements per scan.The probability of the extended target  generating at least one measurement is 1 −  −  .Each target is detected with probability    , giving the effective probability of detection (1 −  −  )   .If the detected probability of each extended target is the same, denoted by   , and the mean number of measurements from each extended target is , the prior probability for   measurements from   extended targets is given by where     is the number of ways of choosing a subset of  elements from   and (  | ,    ) is the probability of   measurements generating from  extended targets.Since the number of measurements from each extended target is not the same and not fixed in the case of  > 1, the form of (  | ,    ) is complicated and consists of  −1   −1 components.The computational complexity for the conditional likelihood in (6) and the prior distribution of the association hypothesis in (10) is () and (min(  ,   )), respectively, and thus the computational complexity for the marginal likelihood in ( 9) is (  ( + min(  ,   ))).The number of association hypotheses increases exponentially with the increase in the number of measurements and extended targets.Thus, in many practical applications, the number of feasible association hypotheses   is prohibitively large, and the marginalisation cannot be performed explicitly.

The Multiple Extended Target Particle Filter
In the condition that the combined state evolution (   |   −1 ) and the marginal likelihood (  |    ) have been known, and a suitable proposal distribution is given, the idea of the particle filter could be applied in tracking multiple extended targets.However, the computational complexity is (  ( + min(  ,   ))) at each time step ( is the number of particles per step) and quickly becomes infeasible for the numbers of extended targets and measurements.
The most of computational burden is generated by marginalising over the association uncertainty in computing the marginal likelihood.To avoid the expensive marginalisation computation, we apply the idea of the particle filtering to jointly estimate the multiple extended target combined state    and the data association  and define the joint proposal distribution for the combined state and association hypothesis (   ,   |   −1 ,  −1 ,   ), from which the combined state and the association hypotheses are sampled.

The Framework of the Multiple Extended Target Particle
Filter.In the Bayesian framework, each update of the joint probability density (   ,   |   ) is preceded by a prediction step and a filtering step.The predicted joint density (   ,   |  −1 ) is calculated by integration where (  |    ) is the prior distribution of the association hypothesis and (   |   −1 ) is the transition density of the combined state, which describes the state evolution dynamic model.According to Bayes' rule, the joint probability density function (   ,   |   ) is given by where (  |    ,   ) is the conditional likelihood as in (6).The Monte Carlo implementation of the formulas ( 14) and (15) could be described as follows: assuming that a set of samples with weights { ,() −1 ,  () −1 ,  () −1 }  =1 is available, approximately distributed according to the marginal filtering distribution at the previous time step (  −1 ,  −1 |  −1 ), at the current time step, the set of new samples for the multiple extended target combined state and association hypothesis are generated from a suitably joint proposal distribution The particle weights are given by The resulting sample set { ,()  ,  ()  ,  ()  }  =1 is approximately distributed by the joint posterior probability distribution at the current time step (   ,   |   ), and resampling could be applied to avoid degeneracy of the particle sets.

The Joint Proposal Distribution.
According to the multiplication theorem on probability, the joint proposal distribution for the multiple extended target combined state and association hypothesis can be written as We assume that the combined state and the association hypothesis at the current time step are independent of the association hypothesis at the previous time step.This assumption is reasonable due to only a weak temporal dependence between the association hypotheses.Thus the joint proposal distribution can be written as a product of the proposal for the association hypothesis and the proposal for the combined state

The Proposal Distribution for the Association Vector.
If the association vector a has been generated from the proposal distribution, the number of clutter measurements   and measurements from extended target   could be calculated, and thus we define the proposal distribution for the association vector instead of the association hypothesis.
Based on the notion of a soft-gating of measurements, the proposal for the association vector is assumed to take a form that factorises sequentially over the individual target associations, which is similar in spirit to the single extended target tracking developed in [30], as shown below where  is the number of measurements and (   |  1:−1  ,    , z   ) is the proposal distribution for the th component of the association vector, the probability of which should be high if the measurement is close to the associated extended target and becomes low as the distance between the measurement and the associated extended target increases.To meet this requirement, we define the proposal distribution for the th component using Bayes' rule where (   =  |  1:−1  , x   ) is defined by the prior probability distribution of the th component as where  0 is the prior probability of the measurement  from the clutter and set to some fixed value,    is the detection probability of the extended target , and   is the mean number of measurements from the extended target .
Since the extended target could generate more than one measurement, for any given measurement, it may be generated from any extended target or the clutter.Thus, we could assume that the prior probability distribution of the th component ( 22) is independent of the association hypotheses for those 1 :  − 1th measurements; i.e., ( . Therefore, the proposal distribution for the th component of the association vector (   |  1:−1  ,    , z   ) is also independent of  1:−1  and thus can be represented as (   |    , z   ).

The Proposal Distribution for the Multiple Extended
Target Combined State.As assumed in Section 2, the individual extended targets evolve independently according to Markovian dynamic models in (3), and thus the proposal distribution function for the multiple extended target combined state factorises over the proposal for individual extended target The simplest state proposal distribution is taken to be the state transition function However, the new particle set generated from the simplest state proposal distribution is seriously dependent on the target evolution model and does not work well in some situations due to taking no account of the current measurements.An alternative method is to select the optimal importance function as the proposal distribution, which minimises the variance of the importance weights where (  | x   ) is the likelihood function of the extended target .For linear and Gaussian models, the optimal importance function is also the Gaussian function, and its analytic form could be obtained.For models with nonlinearity and/or non-Gaussian noise, it is generally not possible to obtain a closed-form expression for the optimal proposal distribution.

Implementation of the Multiple Extended Target Particle
Filter.If the particle filter discussed above is directly implemented, most particles will have the following phenomenon: the resulting joint particle has good estimates for some target' states, and bad estimates for other.However, the weight of this particle is small for penalty of the poor target's states.Therefore, during the resampling stage, this particle would be ignored and the component with good estimates is also abandoned by suppressing the poor estimates.This leads to a great decline of the filter performance.In this case, to maintain the same estimation accuracy, the number of particles required will increase exponentially with the increase of the combined state dimension, which increases as the number of targets increases.
To solve the above problem, the proposed method is as follows.Firstly, the generation and evolution of the multiple extended target combined state particles are decomposed into that of the individual target state particles, and the particle set of each target is resampled according to the weights associated with it.This allows each target to retain particles whose state is better, avoiding the problem generated from the combined state particles evolution directly.The feasibility and rationality are discussed below.
Since the proposal distribution for the th component of the association vector (   |    , z   ), as discussed in Section 3.2.1, is independent of  1:−1  , according to (3), ( 19), (20), and (23), the particle weight represented in (17) can be decomposed into weights related to the individual extended target where ( ()  ) is the prior probability that produces  = }.If a particle has bad estimates for the extended target , the weight associated with this extended target  ,()  is small for the penalty of the poor state and causes the small weight of this particle  ()   .In this case, even with good state estimates for some targets, this particle would be ignored, and the component with good estimates is also abandoned during the resampling stage.
The joint proposal distribution for the multiple extended target combined state and association hypothesis can be written as a product of that for the association vector and that for the combined state, as shown in (18), and the proposal distribution function for the extended target combined state can factorise over the proposal for individual extended target, as shown in (23).Therefore, new particle set of the extended target combined state can be independently sampled from those proposal distributions of individual target at  step.Furthermore, as shown in (26), the particle weight can be decomposed into weights related to the individual extended target, and thus the generation and evolution of extended target combined state particles can be decomposed into that of individual target state particles.Therefore, the particle set of individual extended targets could be resampled individually according to those weights  ,()  , and then the particle set of the multiple extended target combined state is resampled according to weights  ()   .Assuming that the sample set { ,() −1 ,  () −1 }  =1 , which is approximately distributed by the marginal posterior probability (  −1 |  −1 ), is obtained by marginalizing the sample set { ,() −1 ,  () −1 ,  () −1 }  =1 of the joint posterior probability distribution (  −1 ,  −1 |  −1 ), the implementation steps of the proposed multiple extended target particle filter are as follows.
Not only the kinematic state but also the shape information can be obtained from the multiple measurements generated from the extended target.According to the target state evolution model, the measurement source model, and the sensor measurement model, the proposed multiple extended target particle filter can estimate the kinematic state and shape parameter of each extended target.For example, the random hypersurface model can be used to describe the extended target measurement source model [5].

Simulation Results
In this section, we compare the performance of the proposed multiple extended target particle filter with the ET-RHM-GMPHD filter and the ET-RHM-SMCPHD filter.The ET-RHM-GMPHD filter uses the framework of the Gaussian-Mixture probability hypothesis density (GMPHD) filter to avoid the data association [17] and applies the random hypersurface model describing the measurement source distribution for convex-star extended target [5] and employs the particle filter to deal with the nonlinear measurement model in update step.Details for the ET-RHM-GMPHD filter are referred to in [32].The ET-RHM-SMCPHD filter is a sequential Monte Carlo implementation of the extended target probability hypothesis density (ET-PHD) filter based on RHM.

Simulation Setup.
For illustration purposes, we consider a two-dimensional scenario over the surveillance region [−1000, 1000] × [−1000, 1000] (in m).The extended target state is , where [  ,   ] T is the position vector of the centroid, [ ẋ  , ẏ  ] T is the velocity vector, and x   is the shape parameter vector.The shape parameter vector x   , as defined in [5], is [a 0  ,   ] T , which consists of Fourier coefficients of the Fourier series expansion of the radial function for the extended target boundary.Here   is set to 5.
The nonlinear system which has the linear Gaussian dynamics and the nonlinear measurement model is considered.Each extended target follows the linear Gaussian dynamics where k   is process noise with zero-mean and covariance matrix   .
] , 0.03I 11 ), and  = 1 is the sampling period, and  = 1 is the standard deviation of the process noise, and I  denotes the  ×  identity matrix, and 0  is the  ×  zero matrix.The measurement sources are uniformly drawn from the surface of the extended target (the squared scaling factor   is uniformly distributed on the interval [0, 1] [5]).
As illustrated in Figure 1, the measurement model consists of the measurement source model which is described by the random hypersurface model (RHM) and sensor model [5].The measurement z  is a noisy version of the position of measurement source y  , but there is serious nonlinearity between shape parameters of the irregular extended target and measurements.Each extended target has multiple measurement sources, so multiple measurements can be generated.
Based on (33), the following new measurement equation ℎ * (⋅, ⋅, ⋅, ⋅) is obtained: which maps the state   , the measurement noise  , , and the scaling factor  , to the pseudomeasurement 0. The obtained final measurement equation ( 34) is applied to update the kinematic state and shape parameter estimates [5,33] and is used as the measurement model for the star-convex shape target.The measurement noise is set to zero-mean Gaussian with covariance matrix Σ V = diag(0.2,0.2).The scaling factor  is set to a Gaussian function with mean 0.5 and variance 0.02.In each simulation, clutters are generated with a Poisson rate of 10 clutter measurements and each extended target generates measurements with a Poisson rate 10 measurements per scan.The detection probability of the extended target is set to 0.9, and the number of particles per steps is set to 1000.

Results
. In this section, four star-convex extended targets are tracked, and their trajectories are shown in Figure 2. The targets 1 and 2 are born at 1 s.The target 3 is spawned from target 1 at 6 s and disappears at 15 s.The target 4 is born at 20 s.The size of the targets 3 and 4 is about half of the targets 1 and 2, as shown in Figures 3 and 4. The parameters of the shape are a priori set to Gaussian with mean [2, 0 1×10 ]  and covariance matrix 0.03I 11 , i.e., an uncertain circle with radius 2.
To evaluate the of the proposed filter, the root mean squared errors (RMSE) of the position estimates are compared with the corresponding results produced by the ET-RHM-GMPHD filter and the ET-RHM-SMCPHD filter.Because Fourier coefficients with lower indices encode information about the coarse features of the shape and Fourier coefficients with higher indices give information about finer details, their RMSE could not describe the shape estimation error well.In this paper, the Jaccard distance between the estimated shape and the true shape is applied to evaluate the estimated shape of the extended target, as done in [28,34].
The Jaccard distance between two shapes is derived from the Jaccard distance between the samples, defined as [34] where (A) is the region formed by geometric shape A; () is the area of region .In [28], we have proved that the Jaccard distance accords with three conditions of the metric.The Jaccard distance reflects the difference in the shape contour, size, position, and direction of the two shapes.In this paper, the center position of the estimated shape is moved to the center of the actual shape, and then the Jaccard distance is used to measure the shape contour and size and the direction of the estimated shape performance.
The shape estimation of the proposed multiple extended target filter and the two compared filters at some moment are shown as Figures 5 and 6.Figures 7 and 8 of the position estimates and the Jaccard distance between the estimated shape and the true shape, respectively, for the ET-RHM-GMPHD filter and the ET-RHM-SMCPHD and the proposed multiple extended target filter.Figure 7 shows that the proposed filter has better position estimates than the ET-RHM-GMPHD filter and the ET-RHM-SMCPHD filter at the initial several moments of the new target tracking and in the situation that the two extended targets are spatially close to each other, and slightly better in other situation.The reason is as follows: it is difficult to distinguish measurements from those spatially close extended targets using the measurement set partitioning algorithm, on which the performance of the ET-RHM-GMPHD filter and the ET-RHM-SMCPHD depends; however association hypothesis of the proposed particle filter is extracted from its proposal distribution function, and the estimation error caused by the uncertainty of the association hypothesis is small.As shown in Figures 5, 6, and 8, the proposed filter significantly outperforms the ET-RHM-GMPHD filter and the ET-RHM-SMCPHD in the shape estimation performance.The reason for this is that the function, which describes the relationship between measurements and the shape parameters, is severely nonlinear, and the proposed multiple extended target filter, which jointly estimates extended targets and association hypothesis using particle filtering techniques, is specifically presented for the nonlinear non-Gaussian system, and furthermore the estimation error caused by the uncertainty of the association hypothesis is small.

Conclusion
For the nonlinear system, this paper proposed a multiple extended target particle filter which applies the idea of particle filter to jointly estimate the kinematic state and shape parameters of extended targets.Firstly, the Bayesian framework for multiple extended targets is proposed to jointly estimate multiple extended target state and data association.Furthermore, the joint proposal distribution function of multiple extended target combining state and association hypothesis is defined according to the extended target measurement source model, the measurement model, and the target state evolution model.Then, the particle implementation is applied under the Bayesian framework for multiple extended targets, and to improve the performance the generation and evolution of the joint target state particles are decomposed into the generation and evolution of the particles of individual target state.Simulation results show that the proposed multiple extended target particle filter achieves the less error of the shape estimates and enhances the position tracking accuracy at the initial several moments of the new target tracking and when there are spatially close extended targets.

Figure 7 :Figure 8 :
Figure 7: The mean RMSE of the position estimates.
(31)e measurement noise  , , and the scaling factor   , to the measurement z  , , where z  , is the th measurement from the extended target ,    = [   ,    ] T is the position vector of the centroid, and (  , ) = [cos(  , ), sin(  , )] T , and R  (  , ) = [1,cos(  , ), sin(  , ), ..., cos(    , ), sin(    , )].The term   , denotes the angle between the vector from the center to the measurement source y  , and the -axis, but it is unknown in(31).However, we can substitute the unknown value of   , by a proper point estimate.In case of isotropic measurement noise, a proper point estimate φ , is given by the angle between the vector from the current center estimate μ ,−1 to the measurement z  , and the -axis.To reduce the influence of φ, on   , algebraic manipulations are performed on (31) as follows: ,  , ,  , ,  , )