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

Authors: Yun Wang; Guo-ping Hu; Hao Zhou
Journal: Journal of Sensors
 (2016)

---

## 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.


---

## Body

## 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 Section2; 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 setGk at time k, on condition that there is a certain measurement set Zk=Ztt=1k. This process is completed by transforming the prediction part Dk∣k-1· and the update part Dk∣k· of Gk.Firstly, at timek state set of group targets can be represented as (1)Gk=ξkii=1Nξ,k,ξki≜γki,xki,Xki,where Nξ,k is the unknown number of groups, ξk(i) is the state which needs to be estimated of the ith group, measurement rate γk(i)>0, kinematic state xk(i)∈Rnx, and extension state Xk(i)∈S++d. · denotes the potential of set, such that Gk=NG,k.The motion model of the target can be represented as(2)xki=Fk∣k-1xk-1i+vki,where Fk∣k-1=Fk∣k-1⊗Id, Fk∣k-1 is the state transition matrix, Id is a unit matrix of d×d dimension, and vk(i) is white Gaussian noise with zero mean and covariance Δk∣k-1(i)=Qk∣k-1⊗Xk(i). Qk∣k-1 is the covariance matrix of the state space with single dimension.Measurement set of group targets at timek can be represented as (3)Zk=zkjj=1NZ,k,where NZ,k=Zk denotes the number of measurements. So, measurement model of the target can be represented as (4)zkj=Hkxki+wkj,where Hk=Hk⊗Id; Hk is the measurement transition matrix of the state space with single dimension. For example, while d=2, Hk=1,0. wk(j) is white Gaussian noise with zero mean and covariance Xk(i)+R. While measurement noise covariance R is small enough compared to extension state Xk(i), the impact of R can be ignored [2] and covariance turns into Xk(i).Then, assuming the generation of measurement is a process of Poisson point and the measurement rate isγk at time k, hence the predicted probability density function of γk is submitted to Gamma distribution. Therefore, measurement rate is modeled by GAM· in order to obtain the posterior probability density function of γk [8]: (5)pγk∣Zk=pγk∣Zk-1pNZ,k∣γk=βk∣k-1αk∣k-1γkαk∣k-1-1e-βk∣k-1·γkΓαk∣k-1·γkNZ,ke-γkNZ,k!=GAMγk;αk∣k-1+NZ,k,βk∣k-1+1×Γαk∣k-1+NZ,kβk∣k-1αk∣k-1Γαk∣k-1βk∣k-1+1αk∣k-1+NZ,kNZ,k!=GAMγk;αk∣k,βk∣k×Lγαk∣k-1,βk∣k-1,NZ,k,where predicted likelihood function Lγ· is a negative binomial distribution and αk∣k=αk∣k-1+NZ,k,βk∣k=βk∣k-1+1. The expectation and variance of Gamma distribution are Eγk=αk∣k/βk∣k,Varγk=αk∣k/βk∣k2. Assume that the predicted probability density function is a Gamma distribution presented as 0<1/ηk<1 which is a form of index with a fading factor, where αk+1∣k=αk∣k/ηk,βk+1∣k=βk∣k/ηk. The length of efficient window in each prediction is we=1-1/ηk-1, and the expectation of predicted γk remains the same; namely, Eγk+1∣k=αk∣k/βk∣k=Eγk, and variance increases ηk times at each turn; namely, Varγk+1∣k=αk+1∣k/βk+1∣k2=ηkαk∣k/βk∣k2=ηkVarγk.Therefore,N· and IW· are proposed to describe xk and Xk separately [3], and the posterior probability density function is a GGIW distribution; namely, (6)pξk∣Zk=pγk∣Zkpxk∣Xk,ZkpXk∣Zk=GAMγk;αk∣k,βk∣kNxk;mk∣k,Pk∣k⊗Xk×IWXk;vk∣k,Vk∣k=GGIWξk;ζk∣k,where parameters of GGIW component are ζk∣k=αk∣k,βk∣k,mk∣k,Pk∣k,vk∣k,Vk∣k and covariance Pk∣k⊗Xk∈S+nx with Pk∣k∈S+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 βFA,k·S while the monitored d area is S.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] (7)Dk∣k-1ξk=Dk∣k-1sξk+Dkbξk=∫pSξk-1pk∣k-1ξk∣ξk-1×Dk-1∣k-1ξk-1dξk-1+Dkbξk,where pS· is survival probability of the targets, pk∣k-1· is state transition density, Dk∣k-1sξk is the predicted PHD of survival targets, and Dkb· 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] that(8)Dk∣k-1sξk=∫pSξk-1pk∣k-1ξk∣ξk-1×Dk-1∣k-1ξk-1dξk-1=pS∑j=1Jk-1∣k-1wk-1∣k-1j∫GAMγk-1;αk-1∣k-1j,βk-1∣k-1jpk∣k-1γγk∣γk-1dγk-1︸Measurementrate×∫Nxk-1;mk-1∣k-1j,Pk-1∣k-1j⊗Xk×pk∣k-1xxk∣xk-1,Xkdxk-1︸Kinematics×∫IWXk-1;vk-1∣k-1j,Vk-1∣k-1j×pk∣k-1XXk∣Xk-1dXk-1︸Extension,where αk-1∣k-1(j) and βk-1∣k-1(j) denote scalar shape parameter and scalar inverse scale parameter of estimation at time k-1 separately, mk-1∣k-1j and Pk-1∣k-1j denote mean and covariance of estimation group kinematic state at time k-1 separately, and vk-1∣k-1j and Vk-1∣k-1j separately denote degree of freedom and scalar inverse scale matrix that estimate group extension state at time k-1 separately.The predicted PHD can be obtained after integrating the formula above:(9)Dk∣k-1sξk=∑j=1Jk-1∣k-1wk∣k-1jGAMγk;αk∣k-1j,βk∣k-1j×Nxk;mk∣k-1j,Pk∣k-1j⊗Xk×IWXk;vk∣k-1j,Vk∣k-1j,where wk∣k-1(j)=pSwk-1∣k-1(j).


#### 2.2.2. The Birth PHD

The birth PHD can be obtained in accordance with assumption 7 in [9]:(10)Dkbξk=∑j=1Jb,kwb,kjGAMγk;αb,kj,βb,kj×Nxk;mb,kj,Pb,kj⊗Xk×IWXk;vb,kj,Vb,kj.Therefore, Dk∣k-1ξk equates to the sum of the predicted PHD of survival targets (9) and the birth PHD (10) and contains a total of Jk∣k-1=Jk-1∣k-1+Jb,k GGIW components.


#### 2.2.3. Predicted Cardinality Distribution of Group Targets

(11) P k ∣ k - 1 n = ∑ j ′ = 0 n p b n - j × ∑ l = j ∞ l i j ! l - j ! P k - 1 ∣ k - 1 l P S j 1 - P S l - i ,where pb(n-j) denotes the probability for (n-j) new targets to appear between scans k-1 and k.


### 2.3. Update Steps of GGIW-CPHD Filter

Update part of PHD can be obtained according to PHD theory [8]:(12)Dk∣kξk=LZkξkDk∣k-1ξk.Likelihood functionLZkξk is given as (13)LZkξk≜1-1-e-γkpDξk+e-γkpDξk×∑p∠Zkωp∑W∈pγkWdWLγαk∣k-1,βk∣k-1,W∏zk∈Wϕzkξkλkckzk,where estimation of measurement rate γk is the expectation of Gamma distribution; namely, γ^k=αk∣k/βk∣k. λk≜βFA,k·S denotes the mean number of cluttered measurements, and ckzk=1/S. p∠ZK denotes the cluster consequence of thepth nonempty measurement subset partitioned by measurement set ZK. W∈p denotes the Wth measurement subset in the cluster consequence of thepth partition; W denotes measurement number in subsetW. ωp and dW are nonnegative normalization coefficients referring to partition p and subsetW, specified as(14)ωp=∏W∈pdW∑p′∠Zk∏W′∈p′dW′,(15)dW=δW,1+Dk∣k-1pDαk∣kβk∣kWe-αk∣k/βk∣kLγ·∏zk∈Wϕzk·λkckzk,where Lγ· denotes the likelihood function of measurement rate and ϕzk· 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ϕzkξk can be obtained as (16)ϕzkξk=pzk∣ξk=pzk∣xk,Xk=Nzk;Hk⊗Idxk,Xk,where ξk∣k(j) denotes thejth component of GGIW estimated at time k; namely, ξk∣k(j)=αk∣kj,βk∣kj,mk∣kj,Pk∣kj,vk∣kj,Vk∣kj. In addition, we can have [3] (17)P^k∣kj=Pk∣kj⊗Vk∣kjvk∣kj+s-sd-2,X^k∣kj=Vk∣kjvk∣kj-2d-2,where P^k∣kj denotes the covariance matrix of estimation connected to the probability density function of xk,Pk∣kj denotes covariance matrix of estimation connected to marginalized probability density function of xk, and X^k∣kj is the mean estimation value of Xk obeying inverse Wishart distribution.In accordance with the fact whether the target can be detected, the update part of PHD can be partitioned as(18)Dk∣kξk=Dk∣kNDξk+∑p∠Zk∑W∈pωk∣kp,j,WDk∣kDξkj,W,where Dk∣kNDξk is the updated part of the target in leakage alarm and Dk∣kDξkj,W 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 as(19)Dk∣kNDξk=∑j=1Jk∣k-1wk∣kjGAMγk;αk∣kj,βk∣kj×Nxk;mk∣kj,Pk∣kj×IWXk;vk∣kj,Vk∣kj,where wk∣k(j)=1-1-e-αk∣k(j)/βk∣k(j)pDjwk∣k-1(j).


#### 2.3.2. Update Part of the Detected Target

When the target is detected, component parameters can be updated by measurement, specified as(20)Dk∣kDξkj,W=βFA,k-WLkj,W,γLkj,W,x,XGAMγk;αk∣kj,W,βk∣kj,W×Nxk;mk∣kj,W,Pk∣kj,W⊗Xk×IWXk;vk∣kj,W,Vk∣kj,W.The likelihood function of measurement rateLkj,W,γ is expressed as (21)Lkj,W,γ=Γαk∣kj,Wβk∣k-1jαk∣k-1jΓαk∣k-1jβk∣kj,Wαk∣kj,WW!.The likelihood function of measurementLkj,W,x,X is expressed as(22)Lkj,W,x,X=1πWWSk∣k-1j,Wd/2Vk∣k-1jvk∣k-1j/2Vk∣kjvk∣kj,W/2Γdvk∣kj,W/2Γdvk∣k-1j/2,where Sk∣k-1j,W denotes the covariance of innovation, V denotes determinant of V, W denotes the measurement number of subset W, and Γd· denotes multiple Gamma function. Assuming that there are p kinds of measurement partitions at k time, pp denotes the number of thepth partition subsets W. Therefore, the number of GGIW components in the update step of CPHD is Jk∣k=Jk∣k-1+Jk∣k-1∑p=1Ppp.


#### 2.3.3. Posterior Cardinality Distribution of Group Targets [9]

(23) p k ∣ k n = ∑ p ∠ Z ∑ W ∈ p Ψ p , W G k / k - 1 n 0 G F A 0 η W 0,1 / p ρ 1 n - p / n - p ! δ n ≥ p + G F A W 0 ρ n - p + 1 / n - p + 1 ! δ n ≥ p - 1 ∑ p ∠ Z ∑ W ∈ p Ψ p , W l p , W . G k / k - 1 n · and GFA· 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)ηW≜pk/k-1PD·GZW0·∏Z′pZz′/·pFAz′,lp,W=GFA0Gk/k-1pρηWp+GFAW0Gk/k-1p-1ρ,ΨP,W≜∏W′∈P-WηW′,ρ≜pk/k-11-PD·+PD·GZ0·.


## 2.1. Math Model of GGIW Distribution

The purpose of GGIW-CPHD is to estimate the group targets’ state setGk at time k, on condition that there is a certain measurement set Zk=Ztt=1k. This process is completed by transforming the prediction part Dk∣k-1· and the update part Dk∣k· of Gk.Firstly, at timek state set of group targets can be represented as (1)Gk=ξkii=1Nξ,k,ξki≜γki,xki,Xki,where Nξ,k is the unknown number of groups, ξk(i) is the state which needs to be estimated of the ith group, measurement rate γk(i)>0, kinematic state xk(i)∈Rnx, and extension state Xk(i)∈S++d. · denotes the potential of set, such that Gk=NG,k.The motion model of the target can be represented as(2)xki=Fk∣k-1xk-1i+vki,where Fk∣k-1=Fk∣k-1⊗Id, Fk∣k-1 is the state transition matrix, Id is a unit matrix of d×d dimension, and vk(i) is white Gaussian noise with zero mean and covariance Δk∣k-1(i)=Qk∣k-1⊗Xk(i). Qk∣k-1 is the covariance matrix of the state space with single dimension.Measurement set of group targets at timek can be represented as (3)Zk=zkjj=1NZ,k,where NZ,k=Zk denotes the number of measurements. So, measurement model of the target can be represented as (4)zkj=Hkxki+wkj,where Hk=Hk⊗Id; Hk is the measurement transition matrix of the state space with single dimension. For example, while d=2, Hk=1,0. wk(j) is white Gaussian noise with zero mean and covariance Xk(i)+R. While measurement noise covariance R is small enough compared to extension state Xk(i), the impact of R can be ignored [2] and covariance turns into Xk(i).Then, assuming the generation of measurement is a process of Poisson point and the measurement rate isγk at time k, hence the predicted probability density function of γk is submitted to Gamma distribution. Therefore, measurement rate is modeled by GAM· in order to obtain the posterior probability density function of γk [8]: (5)pγk∣Zk=pγk∣Zk-1pNZ,k∣γk=βk∣k-1αk∣k-1γkαk∣k-1-1e-βk∣k-1·γkΓαk∣k-1·γkNZ,ke-γkNZ,k!=GAMγk;αk∣k-1+NZ,k,βk∣k-1+1×Γαk∣k-1+NZ,kβk∣k-1αk∣k-1Γαk∣k-1βk∣k-1+1αk∣k-1+NZ,kNZ,k!=GAMγk;αk∣k,βk∣k×Lγαk∣k-1,βk∣k-1,NZ,k,where predicted likelihood function Lγ· is a negative binomial distribution and αk∣k=αk∣k-1+NZ,k,βk∣k=βk∣k-1+1. The expectation and variance of Gamma distribution are Eγk=αk∣k/βk∣k,Varγk=αk∣k/βk∣k2. Assume that the predicted probability density function is a Gamma distribution presented as 0<1/ηk<1 which is a form of index with a fading factor, where αk+1∣k=αk∣k/ηk,βk+1∣k=βk∣k/ηk. The length of efficient window in each prediction is we=1-1/ηk-1, and the expectation of predicted γk remains the same; namely, Eγk+1∣k=αk∣k/βk∣k=Eγk, and variance increases ηk times at each turn; namely, Varγk+1∣k=αk+1∣k/βk+1∣k2=ηkαk∣k/βk∣k2=ηkVarγk.Therefore,N· and IW· are proposed to describe xk and Xk separately [3], and the posterior probability density function is a GGIW distribution; namely, (6)pξk∣Zk=pγk∣Zkpxk∣Xk,ZkpXk∣Zk=GAMγk;αk∣k,βk∣kNxk;mk∣k,Pk∣k⊗Xk×IWXk;vk∣k,Vk∣k=GGIWξk;ζk∣k,where parameters of GGIW component are ζk∣k=αk∣k,βk∣k,mk∣k,Pk∣k,vk∣k,Vk∣k and covariance Pk∣k⊗Xk∈S+nx with Pk∣k∈S+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 βFA,k·S while the monitored d area is S.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] (7)Dk∣k-1ξk=Dk∣k-1sξk+Dkbξk=∫pSξk-1pk∣k-1ξk∣ξk-1×Dk-1∣k-1ξk-1dξk-1+Dkbξk,where pS· is survival probability of the targets, pk∣k-1· is state transition density, Dk∣k-1sξk is the predicted PHD of survival targets, and Dkb· 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] that(8)Dk∣k-1sξk=∫pSξk-1pk∣k-1ξk∣ξk-1×Dk-1∣k-1ξk-1dξk-1=pS∑j=1Jk-1∣k-1wk-1∣k-1j∫GAMγk-1;αk-1∣k-1j,βk-1∣k-1jpk∣k-1γγk∣γk-1dγk-1︸Measurementrate×∫Nxk-1;mk-1∣k-1j,Pk-1∣k-1j⊗Xk×pk∣k-1xxk∣xk-1,Xkdxk-1︸Kinematics×∫IWXk-1;vk-1∣k-1j,Vk-1∣k-1j×pk∣k-1XXk∣Xk-1dXk-1︸Extension,where αk-1∣k-1(j) and βk-1∣k-1(j) denote scalar shape parameter and scalar inverse scale parameter of estimation at time k-1 separately, mk-1∣k-1j and Pk-1∣k-1j denote mean and covariance of estimation group kinematic state at time k-1 separately, and vk-1∣k-1j and Vk-1∣k-1j separately denote degree of freedom and scalar inverse scale matrix that estimate group extension state at time k-1 separately.The predicted PHD can be obtained after integrating the formula above:(9)Dk∣k-1sξk=∑j=1Jk-1∣k-1wk∣k-1jGAMγk;αk∣k-1j,βk∣k-1j×Nxk;mk∣k-1j,Pk∣k-1j⊗Xk×IWXk;vk∣k-1j,Vk∣k-1j,where wk∣k-1(j)=pSwk-1∣k-1(j).


### 2.2.2. The Birth PHD

The birth PHD can be obtained in accordance with assumption 7 in [9]:(10)Dkbξk=∑j=1Jb,kwb,kjGAMγk;αb,kj,βb,kj×Nxk;mb,kj,Pb,kj⊗Xk×IWXk;vb,kj,Vb,kj.Therefore, Dk∣k-1ξk equates to the sum of the predicted PHD of survival targets (9) and the birth PHD (10) and contains a total of Jk∣k-1=Jk-1∣k-1+Jb,k GGIW components.


### 2.2.3. Predicted Cardinality Distribution of Group Targets

(11) P k ∣ k - 1 n = ∑ j ′ = 0 n p b n - j × ∑ l = j ∞ l i j ! l - j ! P k - 1 ∣ k - 1 l P S j 1 - P S l - i ,where pb(n-j) denotes the probability for (n-j) new targets to appear between scans k-1 and k.


## 2.2.1. The Predicted PHD of Survival Targets

It can be concluded from assumptions 4 and 6 in [9] that(8)Dk∣k-1sξk=∫pSξk-1pk∣k-1ξk∣ξk-1×Dk-1∣k-1ξk-1dξk-1=pS∑j=1Jk-1∣k-1wk-1∣k-1j∫GAMγk-1;αk-1∣k-1j,βk-1∣k-1jpk∣k-1γγk∣γk-1dγk-1︸Measurementrate×∫Nxk-1;mk-1∣k-1j,Pk-1∣k-1j⊗Xk×pk∣k-1xxk∣xk-1,Xkdxk-1︸Kinematics×∫IWXk-1;vk-1∣k-1j,Vk-1∣k-1j×pk∣k-1XXk∣Xk-1dXk-1︸Extension,where αk-1∣k-1(j) and βk-1∣k-1(j) denote scalar shape parameter and scalar inverse scale parameter of estimation at time k-1 separately, mk-1∣k-1j and Pk-1∣k-1j denote mean and covariance of estimation group kinematic state at time k-1 separately, and vk-1∣k-1j and Vk-1∣k-1j separately denote degree of freedom and scalar inverse scale matrix that estimate group extension state at time k-1 separately.The predicted PHD can be obtained after integrating the formula above:(9)Dk∣k-1sξk=∑j=1Jk-1∣k-1wk∣k-1jGAMγk;αk∣k-1j,βk∣k-1j×Nxk;mk∣k-1j,Pk∣k-1j⊗Xk×IWXk;vk∣k-1j,Vk∣k-1j,where wk∣k-1(j)=pSwk-1∣k-1(j).


## 2.2.2. The Birth PHD

The birth PHD can be obtained in accordance with assumption 7 in [9]:(10)Dkbξk=∑j=1Jb,kwb,kjGAMγk;αb,kj,βb,kj×Nxk;mb,kj,Pb,kj⊗Xk×IWXk;vb,kj,Vb,kj.Therefore, Dk∣k-1ξk equates to the sum of the predicted PHD of survival targets (9) and the birth PHD (10) and contains a total of Jk∣k-1=Jk-1∣k-1+Jb,k GGIW components.


## 2.2.3. Predicted Cardinality Distribution of Group Targets

(11) P k ∣ k - 1 n = ∑ j ′ = 0 n p b n - j × ∑ l = j ∞ l i j ! l - j ! P k - 1 ∣ k - 1 l P S j 1 - P S l - i ,where pb(n-j) denotes the probability for (n-j) new targets to appear between scans k-1 and k.


## 2.3. Update Steps of GGIW-CPHD Filter

Update part of PHD can be obtained according to PHD theory [8]:(12)Dk∣kξk=LZkξkDk∣k-1ξk.Likelihood functionLZkξk is given as (13)LZkξk≜1-1-e-γkpDξk+e-γkpDξk×∑p∠Zkωp∑W∈pγkWdWLγαk∣k-1,βk∣k-1,W∏zk∈Wϕzkξkλkckzk,where estimation of measurement rate γk is the expectation of Gamma distribution; namely, γ^k=αk∣k/βk∣k. λk≜βFA,k·S denotes the mean number of cluttered measurements, and ckzk=1/S. p∠ZK denotes the cluster consequence of thepth nonempty measurement subset partitioned by measurement set ZK. W∈p denotes the Wth measurement subset in the cluster consequence of thepth partition; W denotes measurement number in subsetW. ωp and dW are nonnegative normalization coefficients referring to partition p and subsetW, specified as(14)ωp=∏W∈pdW∑p′∠Zk∏W′∈p′dW′,(15)dW=δW,1+Dk∣k-1pDαk∣kβk∣kWe-αk∣k/βk∣kLγ·∏zk∈Wϕzk·λkckzk,where Lγ· denotes the likelihood function of measurement rate and ϕzk· 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ϕzkξk can be obtained as (16)ϕzkξk=pzk∣ξk=pzk∣xk,Xk=Nzk;Hk⊗Idxk,Xk,where ξk∣k(j) denotes thejth component of GGIW estimated at time k; namely, ξk∣k(j)=αk∣kj,βk∣kj,mk∣kj,Pk∣kj,vk∣kj,Vk∣kj. In addition, we can have [3] (17)P^k∣kj=Pk∣kj⊗Vk∣kjvk∣kj+s-sd-2,X^k∣kj=Vk∣kjvk∣kj-2d-2,where P^k∣kj denotes the covariance matrix of estimation connected to the probability density function of xk,Pk∣kj denotes covariance matrix of estimation connected to marginalized probability density function of xk, and X^k∣kj is the mean estimation value of Xk obeying inverse Wishart distribution.In accordance with the fact whether the target can be detected, the update part of PHD can be partitioned as(18)Dk∣kξk=Dk∣kNDξk+∑p∠Zk∑W∈pωk∣kp,j,WDk∣kDξkj,W,where Dk∣kNDξk is the updated part of the target in leakage alarm and Dk∣kDξkj,W 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 as(19)Dk∣kNDξk=∑j=1Jk∣k-1wk∣kjGAMγk;αk∣kj,βk∣kj×Nxk;mk∣kj,Pk∣kj×IWXk;vk∣kj,Vk∣kj,where wk∣k(j)=1-1-e-αk∣k(j)/βk∣k(j)pDjwk∣k-1(j).


### 2.3.2. Update Part of the Detected Target

When the target is detected, component parameters can be updated by measurement, specified as(20)Dk∣kDξkj,W=βFA,k-WLkj,W,γLkj,W,x,XGAMγk;αk∣kj,W,βk∣kj,W×Nxk;mk∣kj,W,Pk∣kj,W⊗Xk×IWXk;vk∣kj,W,Vk∣kj,W.The likelihood function of measurement rateLkj,W,γ is expressed as (21)Lkj,W,γ=Γαk∣kj,Wβk∣k-1jαk∣k-1jΓαk∣k-1jβk∣kj,Wαk∣kj,WW!.The likelihood function of measurementLkj,W,x,X is expressed as(22)Lkj,W,x,X=1πWWSk∣k-1j,Wd/2Vk∣k-1jvk∣k-1j/2Vk∣kjvk∣kj,W/2Γdvk∣kj,W/2Γdvk∣k-1j/2,where Sk∣k-1j,W denotes the covariance of innovation, V denotes determinant of V, W denotes the measurement number of subset W, and Γd· denotes multiple Gamma function. Assuming that there are p kinds of measurement partitions at k time, pp denotes the number of thepth partition subsets W. Therefore, the number of GGIW components in the update step of CPHD is Jk∣k=Jk∣k-1+Jk∣k-1∑p=1Ppp.


### 2.3.3. Posterior Cardinality Distribution of Group Targets [9]

(23) p k ∣ k n = ∑ p ∠ Z ∑ W ∈ p Ψ p , W G k / k - 1 n 0 G F A 0 η W 0,1 / p ρ 1 n - p / n - p ! δ n ≥ p + G F A W 0 ρ n - p + 1 / n - p + 1 ! δ n ≥ p - 1 ∑ p ∠ Z ∑ W ∈ p Ψ p , W l p , W . G k / k - 1 n · and GFA· 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)ηW≜pk/k-1PD·GZW0·∏Z′pZz′/·pFAz′,lp,W=GFA0Gk/k-1pρηWp+GFAW0Gk/k-1p-1ρ,ΨP,W≜∏W′∈P-WηW′,ρ≜pk/k-11-PD·+PD·GZ0·.


## 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 as(19)Dk∣kNDξk=∑j=1Jk∣k-1wk∣kjGAMγk;αk∣kj,βk∣kj×Nxk;mk∣kj,Pk∣kj×IWXk;vk∣kj,Vk∣kj,where wk∣k(j)=1-1-e-αk∣k(j)/βk∣k(j)pDjwk∣k-1(j).


## 2.3.2. Update Part of the Detected Target

When the target is detected, component parameters can be updated by measurement, specified as(20)Dk∣kDξkj,W=βFA,k-WLkj,W,γLkj,W,x,XGAMγk;αk∣kj,W,βk∣kj,W×Nxk;mk∣kj,W,Pk∣kj,W⊗Xk×IWXk;vk∣kj,W,Vk∣kj,W.The likelihood function of measurement rateLkj,W,γ is expressed as (21)Lkj,W,γ=Γαk∣kj,Wβk∣k-1jαk∣k-1jΓαk∣k-1jβk∣kj,Wαk∣kj,WW!.The likelihood function of measurementLkj,W,x,X is expressed as(22)Lkj,W,x,X=1πWWSk∣k-1j,Wd/2Vk∣k-1jvk∣k-1j/2Vk∣kjvk∣kj,W/2Γdvk∣kj,W/2Γdvk∣k-1j/2,where Sk∣k-1j,W denotes the covariance of innovation, V denotes determinant of V, W denotes the measurement number of subset W, and Γd· denotes multiple Gamma function. Assuming that there are p kinds of measurement partitions at k time, pp denotes the number of thepth partition subsets W. Therefore, the number of GGIW components in the update step of CPHD is Jk∣k=Jk∣k-1+Jk∣k-1∑p=1Ppp.


## 2.3.3. Posterior Cardinality Distribution of Group Targets [9]

(23) p k ∣ k n = ∑ p ∠ Z ∑ W ∈ p Ψ p , W G k / k - 1 n 0 G F A 0 η W 0,1 / p ρ 1 n - p / n - p ! δ n ≥ p + G F A W 0 ρ n - p + 1 / n - p + 1 ! δ n ≥ p - 1 ∑ p ∠ Z ∑ W ∈ p Ψ p , W l p , W . G k / k - 1 n · and GFA· 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)ηW≜pk/k-1PD·GZW0·∏Z′pZz′/·pFAz′,lp,W=GFA0Gk/k-1pρηWp+GFAW0Gk/k-1p-1ρ,ΨP,W≜∏W′∈P-WηW′,ρ≜pk/k-11-PD·+PD·GZ0·.


## 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:(25)xk=Fk-1rkxk-1+Gk-1rkvk-1rk,where rk denotes the target motion model in time interval [k-1,k), Fk-1rk and Gk-1rk denote the state transition matrix in model rk and vk-1rk denotes white Gauss noise whose mean value is zero and covariance matrix is Qk-1rk. Model transition matrix is πij=Prrk=j∣rk-1=i.BFG approximation is to replace the above formulas with a BFG distribution:(26)xk=Φk-1xk-1+vk-1,where vk-1~N0,Σk-1; Σk-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 Φk-1 and Σk-1 to make the following formulas established: (27)Exk-1∣A=Exk-1∣B,covxk-1∣A=covxk-1∣B.According to the conclusion in [14], set μk-1≜Exk-1∣A,Θk-1=covxk-1∣A, and we can get(28)Φk-1=∑r=1Mpk,rFk-1r,Σk-1=Θk-Φk-1Θk-1Φk-1T,Θk=∑r=1Mpk,rFk-1rΘk-1+μk-1μk-1TFk-1rT+Gk-1rQk-1rGk-1rT-Φk-1μk-1μk-1TΦk-1T,μk=Φk-1μk-1,pk,r=∑i=1Mπirpk-1,r,where M is the number of tracking models, pk,r is the occurrences probability of model r at k time, and πir represents the model transition matrix.


### 3.2. Prediction Steps of Multiple Models GGIW-CPHD

According to the theory described above, prediction part of GGIW-PHDDk∣k-1ξk at time k includes the GGIW-PHD birth intensity Dkbξk at time k and prediction part of posterior GGIW-PHD intensity Dk∣k-1sξk at time k-1.


#### 3.2.1.Dkbξk

One-step prediction of parameter inDkbξk can be obtained by the prior knowledge.


#### 3.2.2.Dk∣k-1sξk

One-step prediction of parameter inDk∣k-1sξk can be specified as(29)Dk∣k-1sξk=∑j=1Jk-1∣k-1wk∣k-1jGAMγk;αk∣k-1j,βk∣k-1jNxk;mk∣k-1j,Pk∣k-1j⊗Xk×IWXk;vk∣k-1j,Vk∣k-1j,where αk∣k-1(j) and βk∣k-1(j), respectively, represent one-step prediction of scalar shape parameter and scalar inverse scale parameter of thejth component of the measurement rate, mk∣k-1j and Pk∣k-1j⊗Xk, respectively, represent Gauss mean value and one-step prediction of covariance of thejth component, and vk∣kj and Vk∣kj, respectively, represent inverse Wishart degree of freedom parameter and one-step prediction of inverse scale matrix of thejth 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 factor1/ηk-1 is introduced specified as (30)αk∣k-1j=∑r=1Mpk∣k-1,rjηk-1,rαk-1∣k-1j,βk∣k-1j=∑r=1Mpk∣k-1,rjηk-1,rβk-1∣k-1j,where ηk-1,r>1 and the length of effective prediction window of the rth model is lengthk-1,r=1/1-1/ηk-1,r=ηk-1,r/ηk-1,r-1.The prediction part of kinematic state:(31)mk∣k-1j=Φk∣k-1jmk-1∣k-1j,Pk∣k-1j=Φk∣k-1jPk-1∣k-1jΦk∣k-1jT+Σk∣k-1j,Φk∣k-1j=∑r=1Mpk∣k-1,rjFk-1r,where(32)Σk∣k-1j=∑r=1Mpk∣k-1,rjFk-1rPk-1∣k-1j+mk-1∣k-1jmk-1∣k-1jTFk-1rT+Gk-1rQk-1rGk-1rT-Φk∣k-1jmk-1∣k-1jmk-1∣k-1jTΦk∣k-1jT-Φk∣k-1jPk-1∣k-1jΦk∣k-1jT.Therth model of prediction probability pk∣k-1,r(j) is (33)pk∣k-1,rj=∑i=1Mπirpk-1,ij.The prediction part of extension state:(34)vk∣k-1j=∑r=1Mpk∣k-1,rje-Ts/τk-1,rvk-1∣k-1j,Vk∣k-1j=vk∣k-1j-d-1vk-1∣k-1j-d-1Vk-1∣k-1j.


#### 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, Pk∣k-1(j,W) is fixed by a fading factor ϖk(j,W) of STF where Pk∣k-1(j,W) represents one-step prediction covariance which means thejth component of theWth measurement subset. The GGIW of survival targets is mainly modified and the specific calculation formula is (35)Pk∣k-1j,W=ϖkj,WΦk∣k-1j,WPk-1∣k-1j,WΦk∣k-1j,WT+Σk∣k-1j,W,ϖkj,W=ηkj,W,ηkj,W>11,ηkj,W≤1,ηkj,W=trGkj,WtrLkj,W,Gkj,W=N0j,Wk-ξYk∣k-1jW-HΣk∣k-1j,WHT,Lkj,W=HΦk∣k-1j,WPk-1∣k-1j,WΦk∣k-1j,WTHT,N0j,Wk=E1W∑zki∈Wzki-Hmk∣k-1j,W1W∑zki∈Wzki-Hmk∣k-1j,WT,where z~k(j,W)=1/W∑zk(i)∈Wzk(i)-Hmk∣k-1(j,W) 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, Yk∣k-1(j)=Xk∣k-1(j) can be obtained referring to the all dimension process method in [3] and Xk∣k-1(j) is thejth 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].(36)Dk∣kξk=Dk∣kNDξk+∑p∠Zk∑W∈pDk∣kDξk,W.


#### 3.3.1. Update Part of Targets in Leakage Alarm

D k ∣ k N D · is the update part of PHD while the targets are in leakage alarm, specified as (37)Dk∣kND·=∑j=1Jk∣k-1wk∣kjGAMγk;αk∣kj,βk∣kj×Nxk;mk∣kj,Pk∣kj×IWXk;vk∣kj,Vk∣kj,(38)wk∣kj=1-1-e-αk∣kj/βk∣kjpDjwk∣k-1j,(39)αk∣kj=αk∣k-1j,βk∣kj=βk∣k-1j,mk∣kj=mk∣k-1j,Pk∣kj=Pk∣k-1j,vk∣kj=vk∣k-1j,Vk∣kj=Vk∣k-1j.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 measurementz-kW and ZkW scattering matrix of theWth measurement subset are (40)z-kW=1W∑zki∈Wzki,ZkW=∑zki∈Wzki-z-kWzki-z-kWT.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(41)αk∣kj,W=αk∣k-1j+W,βk∣kj,W=βk∣k-1j+1,mk∣kj,W=mk∣k-1j+Kk∣k-1j,Wεk∣k-1j,W,Pk∣kj,W=Pk∣k-1j-Kk∣k-1j,WSk∣k-1j,WKk∣k-1j,WT,vk∣kj,W=vk∣k-1j+W,Vk∣kj,W=Vk∣k-1j+Nk∣k-1j,W+ZkW,where Kk∣k-1j,W represents gain matrix, εk∣k-1j,W represents the innovation vector, Nk∣k-1j,W represents the innovation matrix, and Sk∣k-1j,W represents the innovation covariance. Specified calculation formulas are (42)Kk∣k-1j,W=Pk∣k-1jHkTSk∣k-1j,W-1,εk∣k-1j,W=z-kW-Hkmk∣k-1j,Nk∣k-1j,W=Sk∣k-1j,W-1εk∣k-1j,Wεk∣k-1j,WT,Sk∣k-1j,W=HkPk∣k-1jHkT+Xk∣k-1jW.In addition, the weight of thejth GGIW component of the subset W is (43)wk∣kj,W=ωpdWe-αk∣kj,W/βk∣kj,Wαk∣kj,Wβk∣kj,WβFA,kWpDjLkj,Wwk∣k-1j,where(44)Lkj,W=Lkj,W,γ×Lkj,W,x,X,dW=δW,1+∑j=1Jk∣k-1e-αk∣kj,W/βk∣kj,Wαk∣kj,Wβk∣kj,WβFA,kWpDjLkj,Wwk∣k-1j,where ωp 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:(45)pk∣kn=pk∣k,rjnpk,rj∑r=1Mpk∣k,rjnpk,rj.


### 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 modelr can be assumed as pk,r(j): (46)pk,rj=1M,where M 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(47)pk,rj=pk∣k-1,rjLk,rj,W∑r=1Mpk∣k-1,rjLk,rj,W,where Lk,rj,W is the likelihood function of therth model and can be specified as (48)Lk,rj,W=Γαk∣k,rj,Wβk∣k-1,rjαk∣k-1,rjΓαk∣k-1,rjβk∣k,rj,Wαk∣k,rj,WW!×1πWWSk∣k-1,rj,Wd/2Vk∣k-1,rjvk∣k-1,rj/2Vk∣k,rj,Wvk∣k,rj,W/2Γdvk∣k,rj,W/2Γdvk∣k-1,rj/2,where W represents the measurement number of subsetW and S and V denote the determinant of matrices S andV.One-step prediction parameters of therth model are (49)αk∣k-1,rj=αk-1∣k-1jηk-1,r,βk∣k-1,rj=βk-1∣k-1jηk-1,r,vk∣k-1,rj=e-Ts/τk-1,rvk-1∣k-1j,Vk∣k-1,rj=vk∣k-1,rj-d-1vk-1∣k-1j-d-1Vk-1∣k-1j.Update formulas parameters of therth model are (50)αk∣k,rj,W=αk∣k-1,rj+W,βk∣k,rj,W=βk∣k-1,rj+1,vk∣k,rj,W=vk∣k-1,rj+W,Vk∣k,rj,W=Vk∣k-1,rj+Nk∣k-1,rj,W+ZkW,Nk∣k-1,rj,W=Sk∣k-1,rj,W-1εk∣k-1,rj,Wεk∣k-1,rj,WT,Sk∣k-1,rj,W=HkFk-1rPk-1∣k-1jFk-1rT+Gk-1rQk-1rGk-1rTHkT+Xk∣k-1jW,εk∣k-1,rj,W=z-kW-HkFk-1rmk∣k-1j.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.


## 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:(25)xk=Fk-1rkxk-1+Gk-1rkvk-1rk,where rk denotes the target motion model in time interval [k-1,k), Fk-1rk and Gk-1rk denote the state transition matrix in model rk and vk-1rk denotes white Gauss noise whose mean value is zero and covariance matrix is Qk-1rk. Model transition matrix is πij=Prrk=j∣rk-1=i.BFG approximation is to replace the above formulas with a BFG distribution:(26)xk=Φk-1xk-1+vk-1,where vk-1~N0,Σk-1; Σk-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 Φk-1 and Σk-1 to make the following formulas established: (27)Exk-1∣A=Exk-1∣B,covxk-1∣A=covxk-1∣B.According to the conclusion in [14], set μk-1≜Exk-1∣A,Θk-1=covxk-1∣A, and we can get(28)Φk-1=∑r=1Mpk,rFk-1r,Σk-1=Θk-Φk-1Θk-1Φk-1T,Θk=∑r=1Mpk,rFk-1rΘk-1+μk-1μk-1TFk-1rT+Gk-1rQk-1rGk-1rT-Φk-1μk-1μk-1TΦk-1T,μk=Φk-1μk-1,pk,r=∑i=1Mπirpk-1,r,where M is the number of tracking models, pk,r is the occurrences probability of model r at k time, and πir represents the model transition matrix.


## 3.2. Prediction Steps of Multiple Models GGIW-CPHD

According to the theory described above, prediction part of GGIW-PHDDk∣k-1ξk at time k includes the GGIW-PHD birth intensity Dkbξk at time k and prediction part of posterior GGIW-PHD intensity Dk∣k-1sξk at time k-1.


### 3.2.1.Dkbξk

One-step prediction of parameter inDkbξk can be obtained by the prior knowledge.


### 3.2.2.Dk∣k-1sξk

One-step prediction of parameter inDk∣k-1sξk can be specified as(29)Dk∣k-1sξk=∑j=1Jk-1∣k-1wk∣k-1jGAMγk;αk∣k-1j,βk∣k-1jNxk;mk∣k-1j,Pk∣k-1j⊗Xk×IWXk;vk∣k-1j,Vk∣k-1j,where αk∣k-1(j) and βk∣k-1(j), respectively, represent one-step prediction of scalar shape parameter and scalar inverse scale parameter of thejth component of the measurement rate, mk∣k-1j and Pk∣k-1j⊗Xk, respectively, represent Gauss mean value and one-step prediction of covariance of thejth component, and vk∣kj and Vk∣kj, respectively, represent inverse Wishart degree of freedom parameter and one-step prediction of inverse scale matrix of thejth 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 factor1/ηk-1 is introduced specified as (30)αk∣k-1j=∑r=1Mpk∣k-1,rjηk-1,rαk-1∣k-1j,βk∣k-1j=∑r=1Mpk∣k-1,rjηk-1,rβk-1∣k-1j,where ηk-1,r>1 and the length of effective prediction window of the rth model is lengthk-1,r=1/1-1/ηk-1,r=ηk-1,r/ηk-1,r-1.The prediction part of kinematic state:(31)mk∣k-1j=Φk∣k-1jmk-1∣k-1j,Pk∣k-1j=Φk∣k-1jPk-1∣k-1jΦk∣k-1jT+Σk∣k-1j,Φk∣k-1j=∑r=1Mpk∣k-1,rjFk-1r,where(32)Σk∣k-1j=∑r=1Mpk∣k-1,rjFk-1rPk-1∣k-1j+mk-1∣k-1jmk-1∣k-1jTFk-1rT+Gk-1rQk-1rGk-1rT-Φk∣k-1jmk-1∣k-1jmk-1∣k-1jTΦk∣k-1jT-Φk∣k-1jPk-1∣k-1jΦk∣k-1jT.Therth model of prediction probability pk∣k-1,r(j) is (33)pk∣k-1,rj=∑i=1Mπirpk-1,ij.The prediction part of extension state:(34)vk∣k-1j=∑r=1Mpk∣k-1,rje-Ts/τk-1,rvk-1∣k-1j,Vk∣k-1j=vk∣k-1j-d-1vk-1∣k-1j-d-1Vk-1∣k-1j.


### 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, Pk∣k-1(j,W) is fixed by a fading factor ϖk(j,W) of STF where Pk∣k-1(j,W) represents one-step prediction covariance which means thejth component of theWth measurement subset. The GGIW of survival targets is mainly modified and the specific calculation formula is (35)Pk∣k-1j,W=ϖkj,WΦk∣k-1j,WPk-1∣k-1j,WΦk∣k-1j,WT+Σk∣k-1j,W,ϖkj,W=ηkj,W,ηkj,W>11,ηkj,W≤1,ηkj,W=trGkj,WtrLkj,W,Gkj,W=N0j,Wk-ξYk∣k-1jW-HΣk∣k-1j,WHT,Lkj,W=HΦk∣k-1j,WPk-1∣k-1j,WΦk∣k-1j,WTHT,N0j,Wk=E1W∑zki∈Wzki-Hmk∣k-1j,W1W∑zki∈Wzki-Hmk∣k-1j,WT,where z~k(j,W)=1/W∑zk(i)∈Wzk(i)-Hmk∣k-1(j,W) 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, Yk∣k-1(j)=Xk∣k-1(j) can be obtained referring to the all dimension process method in [3] and Xk∣k-1(j) is thejth component of one-step prediction of extension state.


## 3.2.1.Dkbξk

One-step prediction of parameter inDkbξk can be obtained by the prior knowledge.


## 3.2.2.Dk∣k-1sξk

One-step prediction of parameter inDk∣k-1sξk can be specified as(29)Dk∣k-1sξk=∑j=1Jk-1∣k-1wk∣k-1jGAMγk;αk∣k-1j,βk∣k-1jNxk;mk∣k-1j,Pk∣k-1j⊗Xk×IWXk;vk∣k-1j,Vk∣k-1j,where αk∣k-1(j) and βk∣k-1(j), respectively, represent one-step prediction of scalar shape parameter and scalar inverse scale parameter of thejth component of the measurement rate, mk∣k-1j and Pk∣k-1j⊗Xk, respectively, represent Gauss mean value and one-step prediction of covariance of thejth component, and vk∣kj and Vk∣kj, respectively, represent inverse Wishart degree of freedom parameter and one-step prediction of inverse scale matrix of thejth 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 factor1/ηk-1 is introduced specified as (30)αk∣k-1j=∑r=1Mpk∣k-1,rjηk-1,rαk-1∣k-1j,βk∣k-1j=∑r=1Mpk∣k-1,rjηk-1,rβk-1∣k-1j,where ηk-1,r>1 and the length of effective prediction window of the rth model is lengthk-1,r=1/1-1/ηk-1,r=ηk-1,r/ηk-1,r-1.The prediction part of kinematic state:(31)mk∣k-1j=Φk∣k-1jmk-1∣k-1j,Pk∣k-1j=Φk∣k-1jPk-1∣k-1jΦk∣k-1jT+Σk∣k-1j,Φk∣k-1j=∑r=1Mpk∣k-1,rjFk-1r,where(32)Σk∣k-1j=∑r=1Mpk∣k-1,rjFk-1rPk-1∣k-1j+mk-1∣k-1jmk-1∣k-1jTFk-1rT+Gk-1rQk-1rGk-1rT-Φk∣k-1jmk-1∣k-1jmk-1∣k-1jTΦk∣k-1jT-Φk∣k-1jPk-1∣k-1jΦk∣k-1jT.Therth model of prediction probability pk∣k-1,r(j) is (33)pk∣k-1,rj=∑i=1Mπirpk-1,ij.The prediction part of extension state:(34)vk∣k-1j=∑r=1Mpk∣k-1,rje-Ts/τk-1,rvk-1∣k-1j,Vk∣k-1j=vk∣k-1j-d-1vk-1∣k-1j-d-1Vk-1∣k-1j.


## 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, Pk∣k-1(j,W) is fixed by a fading factor ϖk(j,W) of STF where Pk∣k-1(j,W) represents one-step prediction covariance which means thejth component of theWth measurement subset. The GGIW of survival targets is mainly modified and the specific calculation formula is (35)Pk∣k-1j,W=ϖkj,WΦk∣k-1j,WPk-1∣k-1j,WΦk∣k-1j,WT+Σk∣k-1j,W,ϖkj,W=ηkj,W,ηkj,W>11,ηkj,W≤1,ηkj,W=trGkj,WtrLkj,W,Gkj,W=N0j,Wk-ξYk∣k-1jW-HΣk∣k-1j,WHT,Lkj,W=HΦk∣k-1j,WPk-1∣k-1j,WΦk∣k-1j,WTHT,N0j,Wk=E1W∑zki∈Wzki-Hmk∣k-1j,W1W∑zki∈Wzki-Hmk∣k-1j,WT,where z~k(j,W)=1/W∑zk(i)∈Wzk(i)-Hmk∣k-1(j,W) 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, Yk∣k-1(j)=Xk∣k-1(j) can be obtained referring to the all dimension process method in [3] and Xk∣k-1(j) is thejth 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].(36)Dk∣kξk=Dk∣kNDξk+∑p∠Zk∑W∈pDk∣kDξk,W.


### 3.3.1. Update Part of Targets in Leakage Alarm

D k ∣ k N D · is the update part of PHD while the targets are in leakage alarm, specified as (37)Dk∣kND·=∑j=1Jk∣k-1wk∣kjGAMγk;αk∣kj,βk∣kj×Nxk;mk∣kj,Pk∣kj×IWXk;vk∣kj,Vk∣kj,(38)wk∣kj=1-1-e-αk∣kj/βk∣kjpDjwk∣k-1j,(39)αk∣kj=αk∣k-1j,βk∣kj=βk∣k-1j,mk∣kj=mk∣k-1j,Pk∣kj=Pk∣k-1j,vk∣kj=vk∣k-1j,Vk∣kj=Vk∣k-1j.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 measurementz-kW and ZkW scattering matrix of theWth measurement subset are (40)z-kW=1W∑zki∈Wzki,ZkW=∑zki∈Wzki-z-kWzki-z-kWT.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(41)αk∣kj,W=αk∣k-1j+W,βk∣kj,W=βk∣k-1j+1,mk∣kj,W=mk∣k-1j+Kk∣k-1j,Wεk∣k-1j,W,Pk∣kj,W=Pk∣k-1j-Kk∣k-1j,WSk∣k-1j,WKk∣k-1j,WT,vk∣kj,W=vk∣k-1j+W,Vk∣kj,W=Vk∣k-1j+Nk∣k-1j,W+ZkW,where Kk∣k-1j,W represents gain matrix, εk∣k-1j,W represents the innovation vector, Nk∣k-1j,W represents the innovation matrix, and Sk∣k-1j,W represents the innovation covariance. Specified calculation formulas are (42)Kk∣k-1j,W=Pk∣k-1jHkTSk∣k-1j,W-1,εk∣k-1j,W=z-kW-Hkmk∣k-1j,Nk∣k-1j,W=Sk∣k-1j,W-1εk∣k-1j,Wεk∣k-1j,WT,Sk∣k-1j,W=HkPk∣k-1jHkT+Xk∣k-1jW.In addition, the weight of thejth GGIW component of the subset W is (43)wk∣kj,W=ωpdWe-αk∣kj,W/βk∣kj,Wαk∣kj,Wβk∣kj,WβFA,kWpDjLkj,Wwk∣k-1j,where(44)Lkj,W=Lkj,W,γ×Lkj,W,x,X,dW=δW,1+∑j=1Jk∣k-1e-αk∣kj,W/βk∣kj,Wαk∣kj,Wβk∣kj,WβFA,kWpDjLkj,Wwk∣k-1j,where ωp 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:(45)pk∣kn=pk∣k,rjnpk,rj∑r=1Mpk∣k,rjnpk,rj.


## 3.3.1. Update Part of Targets in Leakage Alarm

D k ∣ k N D · is the update part of PHD while the targets are in leakage alarm, specified as (37)Dk∣kND·=∑j=1Jk∣k-1wk∣kjGAMγk;αk∣kj,βk∣kj×Nxk;mk∣kj,Pk∣kj×IWXk;vk∣kj,Vk∣kj,(38)wk∣kj=1-1-e-αk∣kj/βk∣kjpDjwk∣k-1j,(39)αk∣kj=αk∣k-1j,βk∣kj=βk∣k-1j,mk∣kj=mk∣k-1j,Pk∣kj=Pk∣k-1j,vk∣kj=vk∣k-1j,Vk∣kj=Vk∣k-1j.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 measurementz-kW and ZkW scattering matrix of theWth measurement subset are (40)z-kW=1W∑zki∈Wzki,ZkW=∑zki∈Wzki-z-kWzki-z-kWT.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(41)αk∣kj,W=αk∣k-1j+W,βk∣kj,W=βk∣k-1j+1,mk∣kj,W=mk∣k-1j+Kk∣k-1j,Wεk∣k-1j,W,Pk∣kj,W=Pk∣k-1j-Kk∣k-1j,WSk∣k-1j,WKk∣k-1j,WT,vk∣kj,W=vk∣k-1j+W,Vk∣kj,W=Vk∣k-1j+Nk∣k-1j,W+ZkW,where Kk∣k-1j,W represents gain matrix, εk∣k-1j,W represents the innovation vector, Nk∣k-1j,W represents the innovation matrix, and Sk∣k-1j,W represents the innovation covariance. Specified calculation formulas are (42)Kk∣k-1j,W=Pk∣k-1jHkTSk∣k-1j,W-1,εk∣k-1j,W=z-kW-Hkmk∣k-1j,Nk∣k-1j,W=Sk∣k-1j,W-1εk∣k-1j,Wεk∣k-1j,WT,Sk∣k-1j,W=HkPk∣k-1jHkT+Xk∣k-1jW.In addition, the weight of thejth GGIW component of the subset W is (43)wk∣kj,W=ωpdWe-αk∣kj,W/βk∣kj,Wαk∣kj,Wβk∣kj,WβFA,kWpDjLkj,Wwk∣k-1j,where(44)Lkj,W=Lkj,W,γ×Lkj,W,x,X,dW=δW,1+∑j=1Jk∣k-1e-αk∣kj,W/βk∣kj,Wαk∣kj,Wβk∣kj,WβFA,kWpDjLkj,Wwk∣k-1j,where ωp 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:(45)pk∣kn=pk∣k,rjnpk,rj∑r=1Mpk∣k,rjnpk,rj.


## 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 modelr can be assumed as pk,r(j): (46)pk,rj=1M,where M 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(47)pk,rj=pk∣k-1,rjLk,rj,W∑r=1Mpk∣k-1,rjLk,rj,W,where Lk,rj,W is the likelihood function of therth model and can be specified as (48)Lk,rj,W=Γαk∣k,rj,Wβk∣k-1,rjαk∣k-1,rjΓαk∣k-1,rjβk∣k,rj,Wαk∣k,rj,WW!×1πWWSk∣k-1,rj,Wd/2Vk∣k-1,rjvk∣k-1,rj/2Vk∣k,rj,Wvk∣k,rj,W/2Γdvk∣k,rj,W/2Γdvk∣k-1,rj/2,where W represents the measurement number of subsetW and S and V denote the determinant of matrices S andV.One-step prediction parameters of therth model are (49)αk∣k-1,rj=αk-1∣k-1jηk-1,r,βk∣k-1,rj=βk-1∣k-1jηk-1,r,vk∣k-1,rj=e-Ts/τk-1,rvk-1∣k-1j,Vk∣k-1,rj=vk∣k-1,rj-d-1vk-1∣k-1j-d-1Vk-1∣k-1j.Update formulas parameters of therth model are (50)αk∣k,rj,W=αk∣k-1,rj+W,βk∣k,rj,W=βk∣k-1,rj+1,vk∣k,rj,W=vk∣k-1,rj+W,Vk∣k,rj,W=Vk∣k-1,rj+Nk∣k-1,rj,W+ZkW,Nk∣k-1,rj,W=Sk∣k-1,rj,W-1εk∣k-1,rj,Wεk∣k-1,rj,WT,Sk∣k-1,rj,W=HkFk-1rPk-1∣k-1jFk-1rT+Gk-1rQk-1rGk-1rTHkT+Xk∣k-1jW,εk∣k-1,rj,W=z-kW-HkFk-1rmk∣k-1j.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.


## 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 modelr can be assumed as pk,r(j): (46)pk,rj=1M,where M 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(47)pk,rj=pk∣k-1,rjLk,rj,W∑r=1Mpk∣k-1,rjLk,rj,W,where Lk,rj,W is the likelihood function of therth model and can be specified as (48)Lk,rj,W=Γαk∣k,rj,Wβk∣k-1,rjαk∣k-1,rjΓαk∣k-1,rjβk∣k,rj,Wαk∣k,rj,WW!×1πWWSk∣k-1,rj,Wd/2Vk∣k-1,rjvk∣k-1,rj/2Vk∣k,rj,Wvk∣k,rj,W/2Γdvk∣k,rj,W/2Γdvk∣k-1,rj/2,where W represents the measurement number of subsetW and S and V denote the determinant of matrices S andV.One-step prediction parameters of therth model are (49)αk∣k-1,rj=αk-1∣k-1jηk-1,r,βk∣k-1,rj=βk-1∣k-1jηk-1,r,vk∣k-1,rj=e-Ts/τk-1,rvk-1∣k-1j,Vk∣k-1,rj=vk∣k-1,rj-d-1vk-1∣k-1j-d-1Vk-1∣k-1j.Update formulas parameters of therth model are (50)αk∣k,rj,W=αk∣k-1,rj+W,βk∣k,rj,W=βk∣k-1,rj+1,vk∣k,rj,W=vk∣k-1,rj+W,Vk∣k,rj,W=Vk∣k-1,rj+Nk∣k-1,rj,W+ZkW,Nk∣k-1,rj,W=Sk∣k-1,rj,W-1εk∣k-1,rj,Wεk∣k-1,rj,WT,Sk∣k-1,rj,W=HkFk-1rPk-1∣k-1jFk-1rT+Gk-1rQk-1rGk-1rTHkT+Xk∣k-1jW,εk∣k-1,rj,W=z-kW-HkFk-1rmk∣k-1j.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 isTs=1s. Kinematic state of group targets contains position and velocity; namely, xk=xk,yk,x˙k,y˙k. 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:(51)xk=Fxk-1+Gvk-1.In CV model, there has been(52)FCV=10Ts0010Ts00100001,G=Ts2200Ts22Ts00Ts.In CT models, there have been(53)FCT=10sinΩTsΩ-1-cosΩTsΩ011-cosΩTsΩsinΩTsΩ00cosΩTs-sin⁡ΩTs00sinΩTscosΩTs,G=Ts2200Ts22Ts00Ts,where Ω denotes turning angular velocity of CT models and sets Ω=±0.5rad/s. Observation model of group targets refers to formula (4) where measurement equation H=10. In addition, transformation matrix of model probability is (54)Π=0.90.050.050.050.90.050.050.050.9.Observation location of sensor in the simulation was assumed asx0,y0=0,0m, the generation of clutter is a process of Poisson point which obeys a Poisson distribution with mean value λk=10, and ckzk is the clutter distribution and uniformly distributed in monitoring area. The detection probability of group targets is pD=0.99, and survival probability of the target is pS=0.95. The real extension state of theith group targets is Xki=Rkidiagai2bi2RkiT, where Rki represents a rotation matrix which is related to the angle between long axis of the ellipse and coordinate axis, ai represents long axis of the ellipse, and bi represents short axis of the ellipse. The number of mixture components of the newborn targets’ intension function is Jb,k=3,wb,k(j)=0.1; set the real originate position of the target to be the random set of each newborn target’s components; variances are Pb,k(j)=diag1002452. The extension state of newborn group obeys inverse Wishart distribution and vb,k(j)=7,Vb,k(j)=diag11. The measurement rate of the newborn targets state set is γb,k=10,15,20; parameters of Gamma distribution are αb,k=10,20,30,βb,k=1,1,1. The fading factors of Gamma distribution prediction of the 3 models set in this paper are ηk=5,6,6. In addition, in the clipping and merging of GGIW components, set the GGIW components number maximally as Jmax=100, pruning threshold is TTh=10-5, merging threshold is UTh=102, and the picked state threshold is ETh=0.5.Figure1 describes motion status of 3 group targets in clutter, group 1 is initialed with kinematic state [-254m,509m,0m/s,-50m/s]T, extension state a1,b1=10,5, and measurement rate γk1=20. Group 1 is generated at time k=1, is extinct at time k=30, and is motioned in CT model where Ω=π/16rad/s during 1~9 s, Ω=-π/8rad/s during 10~17 s, and Ω=π/8rad/s during 18~30 s. Group 2 is initialed with kinematic state [-248m,232m,41.4m/s,0m/s]T, extension state a2,b2=15,6, and measurement rate γk2=10. Group 2 is generated at time k=3, is extinct at time k=23, is motioned in CV model during 3~9 s, is motioned in CT model where Ω=-π/8rad/s during 10~17 s, and is motioned in CV model during 18~32 s. Group 3 is generated at time k=13, is extinct at time k=28, is motioned in CV model during 13~16 s, is motioned in CT model where Ω=-π/10rad/s during 17~21 s, and is motioned in CT model where Ω=π/8rad/s 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.Figure 1
Trajectories of group targets.


### 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.Figure 2
Single simulation of GGIW-CPHD algorithm.Figure 3
Single simulation of MM-GGIW-CPHD algorithm.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:Gk=ξk(j)j=1Nξ,k,ξk(j)≜γk(j),xk(j),Xk(j) and estimation state finite sets G^k=ξ^k(i)i=1N^ξ,k,ξ^k(i)≜γ^k(i),x^k(i),X^k(i), where Nξ,k and N^ξ,k separately represent the real and estimation number of group targets. While Nξ,k≤N^ξ,k, the p order OSPA distance is defined [20]:(55)d-pc=1N^ξ,kminπ∈∏N^ξ,k∑j=1Nξ,kdcξkj,ξ^kjp+cpN^ξ,k-Nξ,k1/p,where measurement rate dcγk(j),γ^k(i)=mincX,γk(j)-γ^k(i),cγ=35,p=2; kinematic state dcxk(j),x^k(i)=mincx,xk(j)-x^k(i)2,cx=60,p=2; extension state dcXk(j),X^k(i)=mincX,Xk(j)-X^k(i)F,cX=120,p=2;· represents the absolute value, ·2 represents 2-norm of vector, and ·F represents F-norm of vector.Figures4 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.Figure 4
OSPA distance of GGIW-CPHD algorithm.Figure 5
OSPA distance of MM-GGIW-CPHD algorithm.Figures6 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.Figure 6
OSPA distance of GGIW-CPHD algorithm.Figure 7
OSPA distance of MM-GGIW-CPHD algorithm.Figures8 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.Figure 8
OSPA distance of GGIW-CPHD algorithm.Figure 9
OSPA distance of MM-GGIW-CPHD algorithm.Figures10 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.Figure 10
Group number of GGIW-PHD algorithm.Figure 11
Group number of MM-GGIW-PHD 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.


## 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 isTs=1s. Kinematic state of group targets contains position and velocity; namely, xk=xk,yk,x˙k,y˙k. 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:(51)xk=Fxk-1+Gvk-1.In CV model, there has been(52)FCV=10Ts0010Ts00100001,G=Ts2200Ts22Ts00Ts.In CT models, there have been(53)FCT=10sinΩTsΩ-1-cosΩTsΩ011-cosΩTsΩsinΩTsΩ00cosΩTs-sin⁡ΩTs00sinΩTscosΩTs,G=Ts2200Ts22Ts00Ts,where Ω denotes turning angular velocity of CT models and sets Ω=±0.5rad/s. Observation model of group targets refers to formula (4) where measurement equation H=10. In addition, transformation matrix of model probability is (54)Π=0.90.050.050.050.90.050.050.050.9.Observation location of sensor in the simulation was assumed asx0,y0=0,0m, the generation of clutter is a process of Poisson point which obeys a Poisson distribution with mean value λk=10, and ckzk is the clutter distribution and uniformly distributed in monitoring area. The detection probability of group targets is pD=0.99, and survival probability of the target is pS=0.95. The real extension state of theith group targets is Xki=Rkidiagai2bi2RkiT, where Rki represents a rotation matrix which is related to the angle between long axis of the ellipse and coordinate axis, ai represents long axis of the ellipse, and bi represents short axis of the ellipse. The number of mixture components of the newborn targets’ intension function is Jb,k=3,wb,k(j)=0.1; set the real originate position of the target to be the random set of each newborn target’s components; variances are Pb,k(j)=diag1002452. The extension state of newborn group obeys inverse Wishart distribution and vb,k(j)=7,Vb,k(j)=diag11. The measurement rate of the newborn targets state set is γb,k=10,15,20; parameters of Gamma distribution are αb,k=10,20,30,βb,k=1,1,1. The fading factors of Gamma distribution prediction of the 3 models set in this paper are ηk=5,6,6. In addition, in the clipping and merging of GGIW components, set the GGIW components number maximally as Jmax=100, pruning threshold is TTh=10-5, merging threshold is UTh=102, and the picked state threshold is ETh=0.5.Figure1 describes motion status of 3 group targets in clutter, group 1 is initialed with kinematic state [-254m,509m,0m/s,-50m/s]T, extension state a1,b1=10,5, and measurement rate γk1=20. Group 1 is generated at time k=1, is extinct at time k=30, and is motioned in CT model where Ω=π/16rad/s during 1~9 s, Ω=-π/8rad/s during 10~17 s, and Ω=π/8rad/s during 18~30 s. Group 2 is initialed with kinematic state [-248m,232m,41.4m/s,0m/s]T, extension state a2,b2=15,6, and measurement rate γk2=10. Group 2 is generated at time k=3, is extinct at time k=23, is motioned in CV model during 3~9 s, is motioned in CT model where Ω=-π/8rad/s during 10~17 s, and is motioned in CV model during 18~32 s. Group 3 is generated at time k=13, is extinct at time k=28, is motioned in CV model during 13~16 s, is motioned in CT model where Ω=-π/10rad/s during 17~21 s, and is motioned in CT model where Ω=π/8rad/s 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.Figure 1
Trajectories of group targets.


## 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.Figure 2
Single simulation of GGIW-CPHD algorithm.Figure 3
Single simulation of MM-GGIW-CPHD algorithm.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:Gk=ξk(j)j=1Nξ,k,ξk(j)≜γk(j),xk(j),Xk(j) and estimation state finite sets G^k=ξ^k(i)i=1N^ξ,k,ξ^k(i)≜γ^k(i),x^k(i),X^k(i), where Nξ,k and N^ξ,k separately represent the real and estimation number of group targets. While Nξ,k≤N^ξ,k, the p order OSPA distance is defined [20]:(55)d-pc=1N^ξ,kminπ∈∏N^ξ,k∑j=1Nξ,kdcξkj,ξ^kjp+cpN^ξ,k-Nξ,k1/p,where measurement rate dcγk(j),γ^k(i)=mincX,γk(j)-γ^k(i),cγ=35,p=2; kinematic state dcxk(j),x^k(i)=mincx,xk(j)-x^k(i)2,cx=60,p=2; extension state dcXk(j),X^k(i)=mincX,Xk(j)-X^k(i)F,cX=120,p=2;· represents the absolute value, ·2 represents 2-norm of vector, and ·F represents F-norm of vector.Figures4 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.Figure 4
OSPA distance of GGIW-CPHD algorithm.Figure 5
OSPA distance of MM-GGIW-CPHD algorithm.Figures6 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.Figure 6
OSPA distance of GGIW-CPHD algorithm.Figure 7
OSPA distance of MM-GGIW-CPHD algorithm.Figures8 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.Figure 8
OSPA distance of GGIW-CPHD algorithm.Figure 9
OSPA distance of MM-GGIW-CPHD algorithm.Figures10 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.Figure 10
Group number of GGIW-PHD algorithm.Figure 11
Group number of MM-GGIW-PHD 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.



---

*Source: 7294907-2016-12-25.xml*