Group Targets Tracking Using Multiple Models GGIW-CPHD Based on Best-Fitting Gaussian Approximation and Strong Tracking Filter

Gamma Gaussian inverse Wishart cardinalized probability hypothesis density (GGIW-CPHD) algorithm was always used to track group targets in the presence of cluttered measurements and missing detections. A multiple models GGIW-CPHD algorithm based on best-fitting Gaussian approximation method (BFG) and strong tracking filter (STF) is proposed aiming at the defect that the tracking error of GGIW-CPHD algorithm will increase when the group targets are maneuvering. The best-fitting Gaussian approximation method is proposed to implement the fusion of multiple models using the strong tracking filter to correct the predicted covariance matrix of the GGIW component. The corresponding likelihood functions are deduced to update the probability of multiple tracking models. From the simulation results we can see that the proposed tracking algorithmMM-GGIWCPHD can effectively deal with the combination/spawning of groups and the tracking error of group targets in the maneuvering stage is decreased.


Introduction
In many tracking applications such as maritime surveillance, multiple targets formation, and ground moving convoys, the tracking objects are usually constituted by a series of targets in a coordinated fashion according to a certain constraint in close space, which was called group targets tracking (GTT) problems [1].Usually, we track the group as a whole instead of tracking all members of the targets in the group while the number of group targets is large or the resolution of sensor is limited which makes it unrealistic and unnecessary to track the single target in the group.It should indicate that the tracking objects in this paper are mainly pointed to the group targets: they have a large number of members and the group has a large distant to the sensor which enables the measurement set be similar to a cluster rather than a geometry structure.It indicates that group targets' math model is similar to its extension model while the group has a large number of members, so as to the tracking methods in [2].
Koch [3] proposed to describe the group's entire motion with the kinematic state and extension state of the group, the extension state was modeled by random matrix, and Bayesian recursion of group targets was deduced on the assumption that the sensor measurement error can be approximately ignored compared to extension state error of the group.Feldmann et al. [4] improved the estimate accuracy of Bayesian recursive estimate algorithm by interactive multiple models (IMM) in consideration of the measurement error based on the frame of [3].The kinematic state and shape estimate can be acquired directly in the algorithms, but it is limited to single group targets tracking and no clutter environment.Mahler [5] deduced recursive equation of extension targets in GM-PHD framework firstly, modeled the amount of measurement which is submitted to Poisson distribution with a certain measurement rate.Granström et al. [6] proposed a measurement clustered algorithm based on distance partition and achieved tracking kinematic state of extension targets using ET-GM-PHD algorithm.In [7] the GIW-PHD filtering algorithm was proposed using random

GGIW-CPHD Group Targets
Tracking Algorithm 2.1.Math Model of GGIW Distribution.The purpose of GGIW-CPHD is to estimate the group targets' state set G  at time , on condition that there is a certain measurement set Z  = {Z  }  =1 .This process is completed by transforming the prediction part  |−1 (⋅) and the update part  | (⋅) of G  .
Firstly, at time  state set of group targets can be represented as =1 ,  ()   ≜ ( ()  , x ()  , X ()  ) , where  , is the unknown number of groups,  ()   is the state which needs to be estimated of the th group, measurement rate  ()   > 0, kinematic state x ()  ∈ R  x , and extension state X ()   ∈ S  ++ .| ⋅ | denotes the potential of set, such that |G  | =  , .
The motion model of the target can be represented as where F |−1 =  |−1 ⊗ I  ,  |−1 is the state transition matrix, I  is a unit matrix of  ×  dimension, and k ()  is white Gaussian noise with zero mean and covariance Δ ()  |−1 = Q |−1 ⊗ X ()   .Q |−1 is the covariance matrix of the state space with single dimension.
Measurement set of group targets at time  can be represented as where  , = |Z  | denotes the number of measurements.So, measurement model of the target can be represented as where H  =   ⊗I  ;   is the measurement transition matrix of the state space with single dimension.For example, while is white Gaussian noise with zero mean and covariance X ()   + R. While measurement noise covariance R is small enough compared to extension state X ()   , the impact of R can be ignored [2] and covariance turns into X ()   .Then, assuming the generation of measurement is a process of Poisson point and the measurement rate is   at time , hence the predicted probability density function of   is submitted to Gamma distribution.Therefore, measurement rate is modeled by GAM(⋅) in order to obtain the posterior probability density function of   [8]: where predicted likelihood function L  (⋅) is a negative binomial distribution and  | =  |−1 + , ,  | =  |−1 + 1.The expectation and variance of Gamma distribution are (  ) =  | / | , Var(  ) =  | / 2 | .Assume that the predicted probability density function is a Gamma distribution presented as 0 < 1/  < 1 which is a form of index with a fading factor, where  +1| =  | /  ,  +1| =  | /  .The length of efficient window in each prediction is   = (1 − 1/  ) −1 , and the expectation of predicted   remains the same; namely, ( +1| ) =  | / | = (  ), and variance increases   times at each turn; namely, Var( +1| ) =  +1| / 2 +1| =   ( | / 2 | ) =   Var(  ).Therefore, N(⋅) and IW(⋅) are proposed to describe x  and X  separately [3], and the posterior probability density function is a GGIW distribution; namely, where parameters of GGIW component are  | = ( | ,  | ,  | ,  | , V | ,  | ) and covariance ( | ⊗ X  ) ∈ S   + with  | ∈ S  + .Moreover, assuming that the amount of clutter generated at each time is submitted to Poisson distribution, the appearing positions of clutter are evenly distributed in the monitoring region.So, the mean number of clutter measurement is  , ⋅  while the monitored  area is .
The combination of group target refers to multiple subgroups uniting to form a large group.The spawning of group target refers to one group splitting into two or more subgroups.The combination and spawning of group targets will result in changes of the groups which has a significant impact on the state estimation.We choose the combination and spawning models in [8] based on the Kullback-Leibler criterion.

Prediction
Steps of GGIW-CPHD Filter.GGIW-CPHD filter depends on the basis of assumptions in [9].In accordance with the PHD theory, the predicted part of PHD is [9] where   (⋅) is survival probability of the targets,  |−1 (⋅) is state transition density,   |−1 (  ) is the predicted PHD of survival targets, and    (⋅) is the PHD of newborn targets.

The Predicted PHD of Survival Targets.
It can be concluded from assumptions 4 and 6 in [9] that where  The predicted PHD can be obtained after integrating the formula above: where

The Birth PHD.
The birth PHD can be obtained in accordance with assumption 7 in [9]: Therefore,  |−1 (  ) equates to the sum of the predicted PHD of survival targets (9) and the birth PHD (10) and contains a total of Journal of Sensors

Predicted Cardinality Distribution of Group Targets
where   ( − ) denotes the probability for ( − ) new targets to appear between scans  − 1 and .

Update Steps of GGIW-CPHD Filter.
Update part of PHD can be obtained according to PHD theory [8]: Likelihood function  Z  (  ) is given as where estimation of measurement rate   is the expectation of Gamma distribution; namely, γ =  | / | .  ≜  , ⋅ denotes the mean number of cluttered measurements, and   (z  ) = 1/.∠Z  denotes the cluster consequence of the pth nonempty measurement subset partitioned by measurement set Z  . ∈  denotes the th measurement subset in the cluster consequence of the pth partition; || denotes measurement number in subset W.   and   are nonnegative normalization coefficients referring to partition  and subset W, specified as where L  (⋅) denotes the likelihood function of measurement rate and  z  (⋅) denotes the likelihood function of measurement.It should be illustrated that three methods about measurement partition are proposed in [7]: distance partition, prediction partition, and EM partition, in which the last two methods were estimated on the basis of a relatively precise consequence of first step prediction of group extension state, but the condition cannot be satisfied while the targets are maneuvering.Therefore, the distance partition method [7] is used in this paper for measurement set partition.
According to the provided model, the likelihood function  z  (  ) can be obtained as where where P() | denotes the covariance matrix of estimation connected to the probability density function of x  ,  () | denotes covariance matrix of estimation connected to marginalized probability density function of x  , and X() | is the mean estimation value of X  obeying inverse Wishart distribution.
In accordance with the fact whether the target can be detected, the update part of PHD can be partitioned as where   | (  ) is the updated part of the target in leakage alarm and   | ( , ) is the update part of the target being detected.

Update Part of the Target in Leakage Alarm
. While the target is in leakage alarm, namely, the emerged target has not been detected, component parameters cannot be updated by measurement at the moment, so it can be represented by prediction of the target state, specified as where

Update
The likelihood function of measurement rate L (,),  is expressed as The likelihood function of measurement L (,),x,X  is expressed as where  (,) |−1 denotes the covariance of innovation, || denotes determinant of , || denotes the measurement number of subset , and Γ  (⋅) denotes multiple Gamma function.Assuming that there are  kinds of measurement partitions at  time, |  | denotes the number of the pth partition subsets .Therefore, the number of GGIW components in the update step of CPHD is

Posterior Cardinality Distribution of Group Targets [9]
| () /−1 (⋅) and   (⋅) denote the predicted probability generating function of the state and the probability generating function of the false alarms; the coefficients are defined as follows: (24)

Best-Fitting Gaussian (BFG) Approximation Method.
BFG approximation method [15] is a function estimation method of IMM on the condition of jump Markov linear systems (JMLS).The basic thought is to approximate state transition equation and covariance matrix of process noise with a BFG distribution and keep the target's mean value and variance of prediction state to be the same in the approximation model, so as to achieve rational transform from multiple models to single model.Therefore, it can be obtained according to BFG approximating method.
We take the following JMLS model into account: where   denotes the target motion model in time interval [ − 1, ),  −1 (  ) and  −1 (  ) denote the state transition matrix in model   and k −1 (  ) denotes white Gauss noise whose mean value is zero and covariance matrix is BFG approximation is to replace the above formulas with a BFG distribution: where k −1 ∼ (0, Σ −1 ); Σ −1 is the process noise variance after approximation.Formulas (25) and (26) in JMLS are replaced by factor A and factor B, respectively.The key of the BFG approximation is to find out rational Φ −1 and Σ −1 to make the following formulas established: According to the conclusion in [14], set  −1 ≜ { −1 | }, Θ −1 = cov{ −1 | }, and we can get where  is the number of tracking models,  , is the occurrences probability of model  at  time, and   represents the model transition matrix.

Prediction Steps of Multiple Models GGIW-CPHD.
According to the theory described above, prediction part of GGIW-PHD  |−1 (  ) at time  includes the GGIW-PHD birth intensity    (  ) at time  and prediction part of posterior GGIW-PHD intensity   |−1 (  ) at time  − 1.

𝐷 𝑏 𝑘 (𝜉 𝑘
). One-step prediction of parameter in    (  ) can be obtained by the prior knowledge.

𝐷 𝑠
|−1 (  ).One-step prediction of parameter in   |−1 (  ) can be specified as where  Therefore, it can be obtained according to BFG approximating method.
As prediction probability density function of measurement rate is a Gamma distribution, an index fading factor 1/ −1 is introduced specified as where  −1, > 1 and the length of effective prediction window of the th model is length The prediction part of kinematic state: where The th model of prediction probability The prediction part of extension state: (34)

Fixed One-Step Prediction Covariance Using STF.
It should be indicated that it is available to adjust the filter gain online by STF on the basis of measurement set partition in order to prevent the problem of tracking model mismatched which eventually makes the discrepancy sequence orthogonal.The adaptation of STF is mainly reflected in the identification of prediction variance, namely, fixing prediction covariance with a fading factor [16].In this paper, where z(,) |−1 is innovation,  is a weakening factor, and  ≥ 1.The single dimension state space process in [8] is transformed to all dimensional state space process due to the using of BFG approximation method.Therefore, Y |−1 can be obtained referring to the all dimension process method in [3] and X () |−1 is the jth component of one-step prediction of extension state.

Update Steps of Multiple Models GGIW-CPHD.
In accordance with the fact whether targets can be detected, the updating part of PHD can be partitioned as two parts [17,18].(39) The formulas above indicate that the parameters cannot be updated by measurement while the targets are in leakage alarm, but its' update parameters can be replaced by one-step prediction.

Update Part of the Detected Targets.
Firstly, in update steps, kinematic measurement z   and    scattering matrix of the Wth measurement subset are Update steps are also needed to be conducted in all dimension state space which has been used in prediction steps of BFG approximating method.Specified update steps are It can be analyzed that the model probability is updated by using the likelihood function of different models from each GGIW component in the process of updating, namely, introducing the newest measurement information to update the model probability.

Simulations
4.1.Scene Settings.The cluttered group targets motion scene is two dimensions.Monitoring area is [−1000, 1000] m × [−1000, 1000] m and the number of group targets is unknown and time-varying, the simulation time is 30 s, and the sampling interval of measurement data is   = 1 s.Kinematic state of group targets contains position and velocity; namely, The time evolution model of group targets' kinematic state refers to formula (2); three motion models are applied to describe the targets motion process in the simulation scene of this paper including 1 CV model and 2 CT models [19].The two kinds of models are both satisfied with the following state transform equation: In CV model, there has been In CT models, there have been where Ω denotes turning angular velocity of CT models and sets Ω = ±0.5 rad/s.Observation model of group targets refers to formula (4) where measurement equation Observation location of sensor in the simulation was assumed as ( 0 ,  0 ) = (0, 0) m, the generation of clutter is a process of Poisson point which obeys a Poisson distribution with mean value   = 10, and   (  ) is the clutter distribution and uniformly distributed in monitoring area.The detection probability of group targets is   = 0.99, and survival probability of the target is   = 0.95.The real extension state of the ith group targets is    =    diag([ 2   2  ])( ()  )  , where    represents a rotation matrix which is related to the angle between long axis of the ellipse and coordinate axis,   represents long axis of the ellipse, and   represents short axis of the ellipse.The number of mixture components of the newborn targets' intension function is  , = 3,  ).The measurement rate of the newborn targets state set is  , = 10, 15, 20; parameters of Gamma distribution are  , = 10, 20, 30,  , = 1, 1, 1.The fading factors of Gamma distribution prediction of the 3 models set in this paper are   = 5, 6, 6.In addition, in the clipping and merging of GGIW components, set the GGIW components number maximally as  max = 100, pruning threshold is  ℎ = 10 −5 , merging threshold is  ℎ = 10 2 , and the picked state threshold is  ℎ = 0.5.

Analysis.
The MM-GGIW-CPHD algorithm in this paper is mainly compared with GGIW-CPHD algorithm in [9].Single simulation results of the two algorithms are presented in Figures 2 and 3. We can see from the simulation results that the two algorithms can track multiple group targets well on the whole and can deal with the merging and spawning factors.It can be deduced that the tracking error of GGIW-CPHD algorithm increases apparently while the group targets are maneuvering.However, MM-GGIW-CPHD can make full use of the posterior measurement information in tracking process which can decrease the tracking error while the group targets are maneuvering.
200 orders of Monte-Carlo experiments are conducted.OSPA missing distance is used to estimate the performance of the two algorithms; specified calculation formulas are real state finite sets: G  = { estimation number of group targets.While  , ≤ N, , the  order OSPA distance is defined [20]: where measurement rate  () ( ,  X = 120,  = 2; ‖ ⋅ ‖ represents the absolute value, ‖ ⋅ ‖ 2 represents 2-norm of vector, and ‖ ⋅ ‖  represents -norm of vector. Figures 4 and 5 represent the OSPA distance and the corresponding standard deviation (dotted line) of kinematical states.We can know from the simulation results that while the newborn group targets are emerging or the group targets are maneuvering, the OSPA distances of kinematic state of GGIW-CPHD algorithm and MM-GGIW-CPHD algorithm will both increase.Through the comparison of the results we can find that while the algorithm is stable, the OSPA distance of kinematic state of MM-GGIW-CPHD algorithm is less than that of GGIW-CPHD algorithm on the whole.The fluctuation range of estimation results is also smaller in the MM-GGIW-CPHD algorithm due to the interactive fusion of different models and the introducing of STF filter.
Figures 6 and 7 represent the OSPA distance and the corresponding standard deviation (dotted line) of extension states.We can know from the comparison of simulation results that the OSPA distance of extension state of MM-GGIW-CPHD algorithm is close to that of GGIW-CPHD algorithm on the whole, but the fluctuation range of estimation results of the MM-GGIW-CPHD algorithm is smaller mainly due to the interactive fusion of different models to estimate the extension state that make the stability of the MM-GGIW-CPHD algorithm better.
Figures 8 and 9 represent the measurement rate estimation (solid line) and the corresponding standard deviation (dotted line) of the two algorithms.We can know from the comparison of simulation results that the estimation accuracy of MM-GGIW-CPHD algorithm is better than that of GGIW-CPHD, and the former algorithm has a smaller standard deviation of estimation results and a better robustness.Figures 10 and 11 represent the number estimation of the group targets (solid line) and the corresponding standard deviation (dotted line) of the two algorithms.According to PHD theory, the sum weight of GGIW components is the estimation number of the group targets.We can know from the comparison of simulation results that the estimation results of MM-GGIW-CPHD algorithm are closer to the real number of the group targets than GGIW-CPHD algorithm.What is more, the standard deviation of MM-GGIW-CPHD algorithm is smaller than that of GGIW-CPHD algorithm.
In order to contrast the calculation amount of MM-GGIW-CPHD algorithm and GGIW-CPHD algorithm, we take the comparison with the CPU time used in each step into consideration.In 100 time Monte Carlo simulations, the average time each step needed in single simulation of GGIW-CPHD filter is 10.35 s, as to MM-GGIW-CPHD filter proposed in this paper 13.69 s is needed.As the introduction of the multiple models filter, the calculating time of the MM-GGIW-CPHD algorithm will increase accordingly.With the applying of BFG approximation method in this paper, the calculating amount decreases in a certain degree due to the approximation of multiple models to single model.

Conclusions
The MM-GGIW-CPHD algorithm based on best-fitting Gaussian approximation and strong tracking filter was proposed to solve the large error problem in multiple maneuvering group targets tracking.The best-fitting Gaussian approximation method is proposed to implement the fusion of multiple models using the strong tracking filter to correct the predicted covariance matrix of the GGIW component.The corresponding likelihood functions are deduced to update the probability of multiple tracking models using the newest measurement information.We can see from the simulations that the tracking performance of the MM-GGIW-CPHD algorithm is better than the GGIW-CPHD algorithm in the estimation of kinematic state, extension state, and the number of group targets.
−1 denote scalar shape parameter and scalar inverse scale parameter of estimation at time  − 1 separately,  () −1|−1 and  () −1|−1 denote mean and covariance of estimation group kinematic state at time  − 1 separately, and V −1 separately denote degree of freedom and scalar inverse scale matrix that estimate group extension state at time  − 1 separately.
denotes the jth component of GGIW estimated at time ; namely,  In addition, we can have[3]

1
, respectively, represent one-step prediction of scalar shape parameter and scalar inverse scale parameter of the jth component of the measurement rate,  () |−1 and  () |−1 ⊗ X  , respectively, represent Gauss mean value and one-step prediction of covariance of the jth component, and V () | and  () | , respectively, represent inverse Wishart degree of freedom parameter and one-step prediction of inverse scale matrix of the jth component.

(
,) |−1 is fixed by a fading factor  (,)  of STF where  (,) |−1 represents one-step prediction covariance which means the jth component of the Wth measurement subset.The GGIW of survival targets is mainly modified and the specific calcu-

1 ;
set the real originate position of the target to be the random set of each newborn target's components; variances are  () , = diag([100 2 45 2 ]).The extension state of newborn group obeys inverse Wishart distribution and V = diag([1 1]
Part of the Detected Target.When the target is detected, component parameters can be updated by measurement, specified as