Abstract
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 algorithm MM-GGIW-CPHD can effectively deal with the combination/spawning of groups and the tracking error of group targets in the maneuvering stage is decreased.
1. 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 matrix in consideration of the extension state estimate which was dismissed in the anterior literatures. The likelihood function of multitargets was deduced and the recursive process of the extension state was able to achieve. GGIW-PHD filtering algorithm was proposed to model measurement rate using Gamma distribution and established the merging/deriving judgment standard and math model of group targets based on Kullback-Leibler Divergence (KL-DIV) [8]. In [9] the GGIW-CPHD filtering algorithm was proposed to achieve better tracking accuracy than GGIW-PHD filtering algorithm, but the tracking error will increase apparently once the group targets were maneuvering as group motion was formed by singer model.
Multiple models could track maneuvering targets effectively [10]. In [11, 12] multiple models PHD is applied to multiple maneuvering targets tracking where model interacting had not been brought into filtering process. To solve this problem, BFG approximation method [13] was brought in [14] where the prediction results of multiple models were approximated to single model in the prediction stage of PHD. It was pointed out in [15] that models probability had not been updated in [14] which resulted in the prior probability of models solidified.
Therefore, the multiple models based on best-fitting Gaussian approximation and strong tracking filter will be introduced to GGIW-CPHD algorithm to track maneuvering group targets with unknown number in this paper. In order to improve the tracking accuracy when group targets maneuver, the best-fitting Gaussian approximation method is used to implement the fusion of multiple models in the CPHD prediction stage and a fading factor of strong tracking filter [16] is used to correct the prediction covariance matrix of the GGIW component.
The rest of the paper is organized as follows: GGIW-CPHD algorithm is introduced in Section 2; the MM-GGIW-CPHD algorithm proposed in this paper is introduced in Section 3; Simulations are given in Section 4; Conclusions are drawn in Section 5.
2. 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 at time , on condition that there is a certain measurement set . This process is completed by transforming the prediction part and the update part of .
Firstly, at time state set of group targets can be represented as where is the unknown number of groups, is the state which needs to be estimated of the th group, measurement rate , kinematic state , and extension state . denotes the potential of set, such that .
The motion model of the target can be represented as where , is the state transition matrix, is a unit matrix of dimension, and is white Gaussian noise with zero mean and covariance . is the covariance matrix of the state space with single dimension.
Measurement set of group targets at time can be represented as where denotes the number of measurements. So, measurement model of the target can be represented as where ; is the measurement transition matrix of the state space with single dimension. For example, while , . is white Gaussian noise with zero mean and covariance . While measurement noise covariance is small enough compared to extension state , the impact of can be ignored [2] and covariance turns into .
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 in order to obtain the posterior probability density function of [8]: where predicted likelihood function is a negative binomial distribution and . The expectation and variance of Gamma distribution are . Assume that the predicted probability density function is a Gamma distribution presented as which is a form of index with a fading factor, where . The length of efficient window in each prediction is , and the expectation of predicted remains the same; namely, , and variance increases times at each turn; namely, .
Therefore, and are proposed to describe and separately [3], and the posterior probability density function is a GGIW distribution; namely, where parameters of component are and covariance with . 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.
2.2. 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, is state transition density, is the predicted PHD of survival targets, and is the PHD of newborn targets.
2.2.1. The Predicted PHD of Survival Targets
It can be concluded from assumptions 4 and 6 in [9] thatwhere and denote scalar shape parameter and scalar inverse scale parameter of estimation at time separately, and denote mean and covariance of estimation group kinematic state at time separately, and and separately denote degree of freedom and scalar inverse scale matrix that estimate group extension state at time separately.
The predicted PHD can be obtained after integrating the formula above:where .
2.2.2. The Birth PHD
The birth PHD can be obtained in accordance with assumption 7 in [9]:Therefore, equates to the sum of the predicted PHD of survival targets (9) and the birth PHD (10) and contains a total of GGIW components.
2.2.3. Predicted Cardinality Distribution of Group Targets
where denotes the probability for new targets to appear between scans and .
2.3. Update Steps of GGIW-CPHD Filter
Update part of PHD can be obtained according to PHD theory [8]:
Likelihood function is given as where estimation of measurement rate is the expectation of Gamma distribution; namely, . denotes the mean number of cluttered measurements, and . denotes the cluster consequence of the pth nonempty measurement subset partitioned by measurement set . 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 aswhere denotes the likelihood function of measurement rate and 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 can be obtained as where denotes the jth component of GGIW estimated at time ; namely, . In addition, we can have [3] where denotes the covariance matrix of estimation connected to the probability density function of denotes covariance matrix of estimation connected to marginalized probability density function of , and is the mean estimation value of 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.
2.3.1. 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 aswhere .
2.3.2. Update Part of the Detected Target
When the target is detected, component parameters can be updated by measurement, specified as
The likelihood function of measurement rate is expressed as
The likelihood function of measurement is expressed aswhere 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 .
2.3.3. Posterior Cardinality Distribution of Group Targets [9]
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:
3. Multiple Models GGIW-CPHD Algorithm of Group Tracking Based on BFG and STF
3.1. 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 , and denote the state transition matrix in model and denotes white Gauss noise whose mean value is zero and covariance matrix is . Model transition matrix is .
BFG approximation is to replace the above formulas with a BFG distribution:where ; 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 and to make the following formulas established:
According to the conclusion in [14], set , and we can getwhere is the number of tracking models, is the occurrences probability of model at time, and represents the model transition matrix.
3.2. Prediction Steps of Multiple Models GGIW-CPHD
According to the theory described above, prediction part of GGIW-PHD at time includes the GGIW-PHD birth intensity at time and prediction part of posterior GGIW-PHD intensity at time .
3.2.1.
One-step prediction of parameter in can be obtained by the prior knowledge.
3.2.2.
One-step prediction of parameter in can be specified aswhere and , respectively, represent one-step prediction of scalar shape parameter and scalar inverse scale parameter of the jth component of the measurement rate, and , respectively, represent Gauss mean value and one-step prediction of covariance of the jth component, and and , respectively, represent inverse Wishart degree of freedom parameter and one-step prediction of inverse scale matrix of the jth component.
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 is introduced specified as where and the length of effective prediction window of the th model is .
The prediction part of kinematic state:where
The th model of prediction probability is
The prediction part of extension state:
3.2.3. 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, is fixed by a fading factor of STF where 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 calculation formula is where is innovation, is a weakening factor, and . 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, can be obtained referring to the all dimension process method in [3] and is the jth component of one-step prediction of extension state.
3.3. 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].
3.3.1. Update Part of Targets in Leakage Alarm
is the update part of PHD while the targets are in leakage alarm, specified as 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.
3.3.2. Update Part of the Detected Targets
Firstly, in update steps, kinematic measurement 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 where represents gain matrix, represents the innovation vector, represents the innovation matrix, and represents the innovation covariance. Specified calculation formulas are
In addition, the weight of the jth GGIW component of the subset is wherewhere can be calculated by formula (14). Therefore, it can be deduced that the update of PHD can be calculated by formula (20) where the weight can be calculated by formula (43), and Gamma parameters of measurement rate, Gauss parameters of kinematic state, and inverse Wishart parameters of extension state can be calculated by formulas (41)~(42) [17].
The updated cardinality distribution is specified:
3.4. Update of Models Probability
Here, the conducting of updating is divided into two parts.
3.4.1. Newborn Target Component
Lacking the posterior information, as model probability of newborn target component is mainly influenced by posterior information, the probability of model can be assumed as : where is the number of tracking models in model sets.
3.4.2. Existing Target Component
Model probability is mainly updated by likelihood function, specified as where is the likelihood function of the rth model and can be specified as where represents the measurement number of subset W and and denote the determinant of matrices and V.
One-step prediction parameters of the rth model are
Update formulas parameters of the rth model 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.
4. 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 . 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 beenIn CT models, there have beenwhere denotes turning angular velocity of CT models and sets . Observation model of group targets refers to formula (4) where measurement equation . In addition, transformation matrix of model probability is
Observation location of sensor in the simulation was assumed as , the generation of clutter is a process of Poisson point which obeys a Poisson distribution with mean value , and is the clutter distribution and uniformly distributed in monitoring area. The detection probability of group targets is , and survival probability of the target is . The real extension state of the ith group targets is , 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 ; set the real originate position of the target to be the random set of each newborn target’s components; variances are . The extension state of newborn group obeys inverse Wishart distribution and . The measurement rate of the newborn targets state set is ; parameters of Gamma distribution are . The fading factors of Gamma distribution prediction of the 3 models set in this paper are . In addition, in the clipping and merging of GGIW components, set the GGIW components number maximally as , pruning threshold is , merging threshold is , and the picked state threshold is .
Figure 1 describes motion status of 3 group targets in clutter, group 1 is initialed with kinematic state , extension state , and measurement rate . Group 1 is generated at time , is extinct at time , and is motioned in CT model where during 1~9 s, during 10~17 s, and during 18~30 s. Group 2 is initialed with kinematic state , extension state , and measurement rate . Group 2 is generated at time , is extinct at time , is motioned in CV model during 3~9 s, is motioned in CT model where during 10~17 s, and is motioned in CV model during 18~32 s. Group 3 is generated at time , is extinct at time , is motioned in CV model during 13~16 s, is motioned in CT model where during 17~21 s, and is motioned in CT model where during 22~25 s. Specified measurements are presented in Figure 1, where the black dots represent the group targets’ real measurements, and black “” represents the clutter.
4.2. 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: and estimation state finite sets , where and separately represent the real and estimation number of group targets. While , the order OSPA distance is defined [20]:where measurement rate ; kinematic state ; extension state represents the absolute value, 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.
5. 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.
Competing Interests
The authors declare that they have no competing interests.