# GM-PHD Filter Combined with Track-Estimate Association and Numerical Interpolation

Authors: Jinguang Chen; Bugao Xu; Lili Ma; Rui Sun
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

For the standard Gaussian mixture probability hypothesis density (GM-PHD) filter, the number of targets can be overestimated if the clutter rate is too high or underestimated if the detection rate is too low. These problems seriously affect the accuracy of multitarget tracking for the number and the value of measurements and clutters cannot be distinguished and recognized. Therefore, we proposed an improved GM-PHD filter to tackle these problems. Firstly, a track-estimate association was implemented in the filtering process to detect and remove false-alarm targets. Secondly, a numerical interpolation technique was used to compensate the missing targets caused by low detection rate. At the end of this paper, simulation results were presented to demonstrate the proposed GM-PHD algorithm is more effective in estimating the number and state of targets than the previous ones.


---

## Body

## 1. Introduction

Multitarget tracking uses multisensor integration to provide a variety of observation data through comprehensive processing and to obtain multitarget state estimation, which has a wide range of civilian and military applications. Mahler introduced the theory of random finite set into the problem of multitarget tracking [1, 2]. Based on this theoretical framework, Mahler presented the probability hypothesis density (PHD) to replace the multitarget posterior probability density for multitarget tracking [3]. The PHD is also known as “first-order statistical moment,” which can recursively estimate the number and state of targets, avoid the problem of data association, and effectively reduce the time complexity. However, the recursion of the PHD filtering involves multiple integrals which may not have the closed solution. Under nonlinear and non-Gaussian conditions, the PHD filter does not have its analytical form [4].Subsequently, several researchers enhanced the PHD filter by designing the particle PHD filter (also known as the SMC-PHD) suitable for nonlinear and non-Gaussian conditions [5, 6] and the Gaussian mixture PHD (GM-PHD) filter suitable for linear Gaussian condition [7, 8]. However, The PHD recursion propagates cardinality information with the mean of the cardinality distribution. It effectively approximates the cardinality distribution by Poisson distribution. Since the mean and the variance of Poisson distribution are equal, when the number of targets present is large, the estimated error of the PHD filter is large. To solve this problem, PHD filter can be extended to the cardinality probability hypothesis density (CPHD) filter as seen in [9]. The posterior intensity and cardinality posterior probability distribution of multitarget state sets in the CPHD filter can be obtained through recursion at the same time. Subsequently, The Gaussian mixture cardinality probability hypothesis density (GM-CPHD) filter was given in [10, 11]; and the sequential Monte Carlo cardinality probability hypothesis density (SMC-CPHD) filter was presented in [12].Since the PHD filter can only obtain the target state estimation, it cannot determine which track the state estimation belongs to. To solve this problem, Lin et al. determined whether the estimated state at current time belongs to a previous track through track-estimate association at adjacent time [13]. In [14], the multiple hypotheses tracking approach [15] was used to form tracks among targets in SMC-PHD filter after clutters were removed from measurements. In [16], a new multitarget tracking filter was proposed based on GM-PHD filter, in which the Gaussian items were labeled to distinguish estimates.In addition, the clutter rate and the detection rate have important influences on effectiveness of tracking algorithms. Juang and Burlina showed the connection between the clutter rate and the false alarm rate [17]. When the clutter rate and the detection rate are unknown, Mahler et al. devised the process of PHD or CPHD to adaptively learn the clutter rate and the detection profile during the filtering and obtained the recursive closed-analytical-formula of PHD or CPHD by mixing the Beta distribution and the Gaussian distribution [18]. It is well known that the GM-PHD filter can track the target accurately only when the clutter rate is low and the detection rate is high. In many practical applications, however, these conditions are hardly satisfied. An excessive clutter rate may cause the number of targets to be overestimated, while a low detection rate may cause the number of targets to be underestimated. All these conditions can lead to estimation errors or poor filtering performance.In this paper, we proposed a new algorithm (GM-PHD-TANI) by incorporating track-estimate association and numerical interpolation into GM-PHD filter for solving the multitarget tracking problems under the conditions of a high clutter rate and a low detection rate. The algorithm mainly consists of two stages: track-estimate association and numerical interpolation. In the first stage, most clutters can be removed by the track-estimate association to lessen the overestimate problem. In the second stage, undetected targets can be made up to some extent through numerical interpolations so that the underestimate problem can be alleviated.This paper is organized as follows: Section2 is the problem formulation, including single target filtering problem, the theory of random finite set, and multitarget tracking model; Section 3 includes the PHD filter and GM-PHD filter; Section 4 is the proposed GM-PHD-TANI algorithm; and the experimental results and analysis are given in Section 5.


## 2. Problem Description

### 2.1. Single Target Bayesian Filter

Assume that, in the state spaceRnx, xk is the state vector and nx is the dimension of the state vector. The state transition equation is(1)xk=fk∣k-1xk-1+ωk-1,where ωk-1~N(0,Qk) is a process noise with a Gaussian distribution whose mean is zero and covariance is Qk.In the measurement spaceRnz, the measurement vector is zk and the dimension of the measurement vector is nz. The measurement equation is(2)zk=gkxk+υk,where υk~N(0,Rk) is the measurement noise with a Gaussian distribution whose mean is zero and covariance is Rk.The probability density functionpkxk∣z1:k of the state xk, decided by measurements z1:k=(z1,…,zk), is the posterior probability density function. According to the prior distribution, p0(x0), the posterior probability density function can be deduced via the Bayes theory; that is,(3)pk∣k-1xk∣z1:k-1=∫fk∣k-1xk∣xpk-1x∣z1:k-1dx,pkxkz1:k=gkzk∣xkpk∣k-1xk∣z1:k-1∫gkzk∣xpk∣k-1x∣z1:k-1dx.p
k
·
∣
z
1
:
k contains the relevant information about the state xk at time k. The estimated information can be obtained by minimizing the mean square error of state estimation or by maximizing the posterior probability pk(xk∣z1:k).


### 2.2. Multitarget Bayesian Filter Based on RFS

In the RFS (Random Finite Set) theory, both the value of an element and the cardinality of a set are random. In a multitarget tracking problem, the birth, spawning, and disappearance of a target and the total number and the states of all targets change over time. Because of a false alarm, miss detection, clutters, and measurements received by sensors are changeable over time as well. Thus, the state set and the measurement set of the targets are a set of random variables.Assume that, at timek, Xk={xk,1,…,xk,Nk} is the multitarget state set, Nk is the number of targets, Zk={zk,1,…,zk,Mk} is the measurement set, and Mk is the number of measurements. Among Xk and Zk, xk,i denotes the ith target’s state, and zk,j denotes the jth target’s measurement:(4)Xk=xk,1,…,xk,Nk⊂χ,Zk=xk,1,…,xk,Mk⊂Z,where χ⊆Rnx and Z⊆Rnz denote the state space and the measurement space, respectively. The measurement set includes clutters. Assume the probability density of clutter is ck(·), its number obeys Poisson distribution whose mean is λc, the state of a target continues from time k-1 to k with probability Ps(xk-1) and disappears with probability 1-Ps(xk-1), and the evolution of each target is independent. According to the RFS theory, the evolution of a single target can be described as a RFS S(xk-1). The new birth target set at time k is expressed as Γk, and the spawned target derived from the survival target is denoted as B(x). The state model of multiple targets is(5)Xk=⋃x∈Xk-1Sx∪⋃x∈Xk-1Bx∪Γk.The received measurement set is the mixture of the target RFS,T(xk), with the clutter RFS, Kk; that is,(6)Zk=⋃x∈XkTx∪Kk.LetZ1:k={Z1,…,Zk}. The multitarget posterior density function can be obtained recursively from the optimal multitarget Bayesian filtering through the following formulas: (7)pk∣k-1Xk∣Z1:k-1=∫fk∣k-1Xk∣Xpk-1X∣Z1:k-1μsdX,pkXk∣Z1:k=gkZk∣Xkpk∣k-1Xk∣Z1:k-1∫gkZk∣Xpk∣k-1X∣Z1:k-1μsdX,where μs is a proper measure factor in space χ [6]. Equations (7) involve multiple integrations in space χ, which are often difficult, sometimes even infeasible, to calculate. Their computational complexity dramatically increases with the increase of the target number. Thus, reasonable approximations become inevitable in practice.


## 2.1. Single Target Bayesian Filter

Assume that, in the state spaceRnx, xk is the state vector and nx is the dimension of the state vector. The state transition equation is(1)xk=fk∣k-1xk-1+ωk-1,where ωk-1~N(0,Qk) is a process noise with a Gaussian distribution whose mean is zero and covariance is Qk.In the measurement spaceRnz, the measurement vector is zk and the dimension of the measurement vector is nz. The measurement equation is(2)zk=gkxk+υk,where υk~N(0,Rk) is the measurement noise with a Gaussian distribution whose mean is zero and covariance is Rk.The probability density functionpkxk∣z1:k of the state xk, decided by measurements z1:k=(z1,…,zk), is the posterior probability density function. According to the prior distribution, p0(x0), the posterior probability density function can be deduced via the Bayes theory; that is,(3)pk∣k-1xk∣z1:k-1=∫fk∣k-1xk∣xpk-1x∣z1:k-1dx,pkxkz1:k=gkzk∣xkpk∣k-1xk∣z1:k-1∫gkzk∣xpk∣k-1x∣z1:k-1dx.p
k
·
∣
z
1
:
k contains the relevant information about the state xk at time k. The estimated information can be obtained by minimizing the mean square error of state estimation or by maximizing the posterior probability pk(xk∣z1:k).


## 2.2. Multitarget Bayesian Filter Based on RFS

In the RFS (Random Finite Set) theory, both the value of an element and the cardinality of a set are random. In a multitarget tracking problem, the birth, spawning, and disappearance of a target and the total number and the states of all targets change over time. Because of a false alarm, miss detection, clutters, and measurements received by sensors are changeable over time as well. Thus, the state set and the measurement set of the targets are a set of random variables.Assume that, at timek, Xk={xk,1,…,xk,Nk} is the multitarget state set, Nk is the number of targets, Zk={zk,1,…,zk,Mk} is the measurement set, and Mk is the number of measurements. Among Xk and Zk, xk,i denotes the ith target’s state, and zk,j denotes the jth target’s measurement:(4)Xk=xk,1,…,xk,Nk⊂χ,Zk=xk,1,…,xk,Mk⊂Z,where χ⊆Rnx and Z⊆Rnz denote the state space and the measurement space, respectively. The measurement set includes clutters. Assume the probability density of clutter is ck(·), its number obeys Poisson distribution whose mean is λc, the state of a target continues from time k-1 to k with probability Ps(xk-1) and disappears with probability 1-Ps(xk-1), and the evolution of each target is independent. According to the RFS theory, the evolution of a single target can be described as a RFS S(xk-1). The new birth target set at time k is expressed as Γk, and the spawned target derived from the survival target is denoted as B(x). The state model of multiple targets is(5)Xk=⋃x∈Xk-1Sx∪⋃x∈Xk-1Bx∪Γk.The received measurement set is the mixture of the target RFS,T(xk), with the clutter RFS, Kk; that is,(6)Zk=⋃x∈XkTx∪Kk.LetZ1:k={Z1,…,Zk}. The multitarget posterior density function can be obtained recursively from the optimal multitarget Bayesian filtering through the following formulas: (7)pk∣k-1Xk∣Z1:k-1=∫fk∣k-1Xk∣Xpk-1X∣Z1:k-1μsdX,pkXk∣Z1:k=gkZk∣Xkpk∣k-1Xk∣Z1:k-1∫gkZk∣Xpk∣k-1X∣Z1:k-1μsdX,where μs is a proper measure factor in space χ [6]. Equations (7) involve multiple integrations in space χ, which are often difficult, sometimes even infeasible, to calculate. Their computational complexity dramatically increases with the increase of the target number. Thus, reasonable approximations become inevitable in practice.


## 3. Gaussian Mixture Probability Hypothesis Density Filter

As an extension of the Bayes filtering under the framework of RFS, the PHD filter uses the PHD form of multitarget posterior probability density, that is, its first-order moment of estimation, to replace the posterior probability density for estimating target’s number and state recursively. In spaceχ, the probability distribution of RFS X is P, and the first-order moment is denoted as v, which is a negative function called the intensity function. For a region S in space χ(S⊆χ),(8)∫X∩SPdX=∫Svxdx.Ifvk-1∣k-1(xk-1∣Z1:k-1) denotes the posterior intensity of the multitargets at time k-1, then the predicted intensity at time k is(9)
v
k
∣
k
-
1
x
k
∣
Z
1
:
k
-
1
=
γ
k
x
k
+
∫
p
S
,
k
∣
k
-
1
x
k
-
1
f
k
∣
k
-
1
x
k
∣
x
k
-
1
+
β
k
∣
k
-
1
x
k
∣
x
k
-
1
v
k
-
1
∣
k
-
1
x
k
-
1
∣
Z
1
:
k
-
1
d
x
k
-
1
,where γk(xk) is the intensity of the birth target, βk∣k-1(xk∣xk-1) is the intensity of the spawned target, fk∣k-1(xk∣xk-1) is the state transition probability density of a single target, and pS,k∣k-1(xk-1) is the survival probability of a target from time k-1 to k.The target number at one prediction step is(10)Nk∣k-1=∫vk∣k-1xk∣Z1:k-1dx=Nk∣k-1Γ+Nk∣k-1S+Nk∣k-1B,where Nk∣k-1Γ=∫γk(x)dx is the number of birth targets; Nk∣k-1S=∫pS,k∣k-1(ξ)fk∣k-1(x∣ξ)vk-1∣k-1(ξ∣Z1:k-1)dξ is the number of survival targets; and Nk∣k-1B=∫βk∣k-1(x∣ξ)vk-1∣k-1(ξ∣Z1:k-1)dξ is the number of spawned targets.The PHD is updated as follows when the set of measurementsZk is received:(11)vk∣kxk∣Z1:k=1-pD,kxkvk∣k-1xk∣Z1:k-1+∑zk∈ZkpD,kxkgkzk∣xkvk∣k-1xk∣Z1:k-1κkzk+∫pD,kxkgkzk∣xkvk∣k-1xk∣Z1:k-1dxk,where κk(zk) is the intensity of clutter, gk(zk∣xk) is the likelihood function of a single measurement, and pD,k(xk) is the probability of detection. Correspondingly, the target number is updated as follows:(12)Nk=∫vk∣kxk∣Z1:kdx=Nk∣k-1-∫pD,kxvk∣k-1xk∣Z1:k-1dx+∑z∈Zk∫pD,kxgkz∣xvk∣k-1xk∣Z1:k-1dxκkz+∫pD,kξgkz∣ξvk∣k-1ξ∣Z1:k-1dξ.The Gaussian mixture form of a PHD is formed by Gaussian components, that is,{wk(i),mk(i),Pk(i)}i=1Jk. In the linear Gaussian system, the posterior intensity at time k can be constructed with Gaussian components in a linear combination with different weights:(13)vkx=∑i=1JkwkiNx;mki,Pki,where Jk is the number of Gaussian components, and wk(i), mk(i), and Pk(i), respectively, denote weight, state, and covariance of the ith Gaussian component.Assume that the posterior intensity of multitargets at timek-1 is expressed in the following Gaussian mixture form:(14)vk-1x=∑i=1Jk-1wk-1iNx;mk-1i,Pk-1i.The PHD of multitargets at time k can be also expressed in a Gaussian mixture form:(15)vk∣k-1x=vS,k∣k-1x+vB,k∣k-1x+γkx=PS,k∑j=1Jk-1wk-1jNx;mS,k-1j,PS,k-1j+∑j=1Jk-1∑l=1JB,kwk-1jwB,klNx;mB,k-1j,l,PB,k-1j,l+∑i=1Jγ,kwγ,kiNx;mγ,ki,Pγ,ki=∑l=1Jk∣k-1wk∣k-1iNx;mk∣k-1i,Pk∣k-1i.Then, the posterior PHD of multitargets at time k can be achieved as follows:(16)vkx=1-pD,kvk∣k-1x+∑z∈ZkvD,kx;z.


## 4. GM-PHD Filter with Track-Estimate Association and Numerical Interpolation

### 4.1. Implementation of GM-PHD

Assume that the state model of targets and the measurement model of sensors are linear Gaussian; that is,(17)pxk∣k-1xk-1=Nxk∣k-1;Fk-1xk-1,Qk-1,pzk∣xk=Nzk;Hkxk,Rk,where Fk-1 is the state transition matrix, Qk-1 is the covariance matrix of process noise, Hk is the measurement matrix of the sensors, and Rk is the covariance matrix of measurement noise. The PHD of birth targets and spawned targets are written in the following Gaussian mixture forms: (18)γkx=∑i=1Jγ,kwγ,kiNx;mγ,ki,Pγ,ki,(19)βk∣k-1x∣ξ=∑j=1Jβ,kwβ,kjNx;Fβ,k-1jξ+dβ,k-1j,Qβ,k-1j,where Jγ,k, wγ,k(i), mγ,k(i), and Pγ,k(i) are the numbers of Gaussian components, weights, mean, and covariance of birth targets, and Jβ,k, wβ,k(j), Fβ,k-1(j), dβ,k-1(j), Fβ,k-1(j)ξ+dβ,k-1(j), and Qβ,k-1(j) are the numbers of Gaussian components, weights, state transition matrix, noise disturbing items, mean, and covariance of spawned targets.Based on the posterior intensity at timek-1 in (14), the posterior intensity at time k can also be predicted by(20)vk∣k-1x=vS,k∣k-1x+vβ,k∣k-1x+γkx,where γk(x) is given in (18) and vS,k∣k-1(x) is the intensity of the survival target which can be expressed as(21)vS,k∣k-1x=pS,k∑j=1Jk-1wk-1jNx;mS,k∣k-1j,PS,k∣k-1j,where(22)mS,k∣k-1j=Fk-1mk-1j,PS,k∣k-1j=Qk-1+Fk-1Pk-1jFk-1T.The intensity of the spawned target, vβ,k∣k-1(x), is(23)vβ,k∣k-1x=∑j=1Jk-1∑i=1Jβ,kwk-1jwβ,kiNx;mβ,k∣k-1j,i,Pβ,k∣k-1j,i,mβ,k∣k-1j,i=Fβ,k-1imk-1j+dβ,k-1i,Pβ,k∣k-1j,i=Qβ,k-1i+Fβ,k-1iPk-1jFβ,k-1iT.According to (20), the predicted PHD can be written in a Gaussian mixture form:(24)vk∣k-1x=∑i=1Jk∣k-1wk∣k-1iNx;mk∣k-1i,Pk∣k-1i.Then, the posterior intensity function at time k also follows a Gaussian mixture form:(25)vkx=1-pD,kvk∣k-1x+∑z∈ZkvD,kx;z,where z∈Zk are the measurements obtained at time k and(26)vD,kx;z=∑j=1Jk∣k-1pD,kwk∣k-1jqkjzκkz+pD,k∑i=1Jk∣k-1wk∣k-1iqkiz·Nx;mk∣kjz,Pk∣kj,qkjz=Nz;Hkmk∣k-1j,HkPk∣k-1jHkT+Rk,mk∣kjz=mk∣k-1j+Kkjz-Hkmk∣k-1j,Pk∣kj=I-KkjHkPk∣k-1j,Kkj=Pk∣k-1jHkTHkPk∣k-1jHkT+Rk-1.A large number of the Gaussian components will bring a heavy burden into the computation. Usually, a merging and pruning stage is needed to deal with this problem [8]. In order to find targets in the Gaussian components, we set a weight threshold wTh (usually chosen as 0.5). If the weight is larger than wTh, the corresponding Gaussian component is regarded as a target; that is,(27)X^k=m~ki:w~ki>wTh,i=1,…,Jk.


### 4.2. Track-Estimate Association

By associating a single target track at the previous time with the state of the target at the current time, the overestimation of the target number in the GM-PHD can be restrained. First, RFSX~k at the current time is obtained according to the conational GM-PHD filter, and the track set Ωk-1 at the previous time is already generated, and then track-estimate association can be done. Due to the randomness of clutter, the error state generated by the clutter also has randomness, so we can eliminate the error of state through the existing track.Figure1 is an illustration of the track-estimate association from the time k-1 to k. In Figure 1, m~k-1(1),m~k-1(2),…,m~k-1(Lk-1) are the end nodes of tracks in Ωk-1 formed at time k-1. The end node is the same as a state vector m~k-1(i)=[xk-1(i),x˙k-1(i),yk-1(i),y˙k-1(i)]T where (xk-1(i),yk-1(i)) denote the location of the ith target and (x˙k-1(i),y˙k-1(i)) denote the velocity of the ith target. (tk-11,tk-12,…,tk-1Lk-1) denotes the end time of tracks at time k-1. It should be pointed out that the end time is not always equal to k-1 because some tracks may not be associated with the estimates at previous time. Δtki denotes a time interval from the end time of the ith track to time k; that is, Δtki=k-tk-1i, i=1,2,…,m. Δtki is the interruption time of the ith track. m~k∣k-1(1),m~k∣k-1(2),…,m~k∣k-1(l) is the predicting vector corresponding to the end time state of each track at time k-1, and m~k∣k-1(i)=[xk∣k-1(i),x˙k∣k-1(i),yk∣k-1(i),y˙k∣k-1(i)]T. If targets move with a constant speed, the prediction of the ith track can be calculated by(28)xk∣k-1i=xk-1i+Δtki·x˙k-1i,yk∣k-1i=yk-1i+Δtki·y˙k-1i,x˙k∣k-1i=x˙k-1i,y˙k∣k-1i=y˙k-1i.Figure 1
Track-estimate association.Let the prediction set after alignment beB={m~k∣k-1(1),m~k∣k-1(2),…,m~k∣k-1(l)} and let a set of state estimations obtained via GM-PHD at time k be C={m~k(1),m~k(2),…,m~k(Nk)}. Furthermore, assume that the permissible maximum interruption steps of a track is S, the minimum estimated number is anum, the distance threshold is dth, and temporary track queue Ωk=⌀.Step 1.
Judge whether a track is associated with a target or not. For a target, track should be continuous strictly; however, because of the detection ratepD<1, some targets are missed in detection. If Δtki>S, the association is not needed, and this track is terminated and pushed into the terminated track queue Ωend. If Δtki≤S, state estimation m~k(i) at time k should be associated with Ωk-1.Step 2.
Compute the distances among elements between setB and set C, and the results are used to construct a matrix D. At this time, the data association is done via Hungary algorithm so that the whole distance cost will be smallest. In the matrix of the costs D, its element Dij denotes the location distance between the ith target’s state m~k(i) in set C and the jth predicted state m~k∣k-1(j) in set B. It belongs to the balance assignment problem when l=Nk. Otherwise, it belongs to the unbalanced assignment problem. For an unbalanced assignment problem, it can be converted to a balance assignment problem via the importing dummy variable (such as zero or infinite) and turning the matrix into a square matrix. So we only need to consider the balance assignment. In this case, the matrix of the costs D is(29)D=DijNk×Nk=D11D12⋯D1NkD21D22⋯D2Nk⋮⋮⋮DNk1DNk2⋯DNkNk.Suppose H is an incidence matrix, where Hij indicates the association between the ith target’s state m~k(i) and the jth predicted state m~k∣k-1(j). When Hij=1, the two targets are associated with each other; when Hij=0, they are not associated. An associated matrix can be obtained by solving the following optimization problem:(30)zmin=minHij⁡∑i=1n∑j=1nDijHij,where zmin denotes the minimum overall benefit.The solution of (30) may be achieved through the Hungary algorithm [19, 20], which is described briefly as below. First, let each row of the matrix of the costs D be minus the smallest element in the row and get the result matrix D′ in which each row has at least one element equal to 0 and all elements satisfy Dij′≥0. Second, the same process is used to deal with columns of matrix D′. We get a new matrix of the costs from the original, in which there is at least one 0 element in each column and each row. Repeat the above two steps until there is only one nonzero element in each column and each row, which is the optimal solution of the assignment problem.We can obtain the matching matrixu according to the final matrix of the costs D′′. For the size of l and Nk, there are three kinds of cases for the number of target’s states and predicted states: l<Nk, l=Nk, and l>Nk. When l=Nk, they can be directly matched. When l<Nk, they can be matched because the number of state estimations m~k(i) matches the dummy variable. Here, these target’s states can be regarded as starting points, and new tracks Ωnew={m~k(i)} can be created. When l>Nk, they can be matched according to the number of predicted states m~k∣k-1(j). In this case, the track is interrupted if the predicted state is matched with a dummy variable according to the cost matrix. This track is regarded as an effective track: Ωk={Ωk,Ωk-1j}.Step 3.
Calculate the Mahalanobis distance between matched elements in setB and set C [21]. d(m~k∣k-1(i),m~k(j)) denotes the square of the distance judging whether its value is located in the corresponding confidence region of χ2 distribution or not.The Mahalanobis distance satisfiesd(m~k∣k-1(i),m~k(j))≤2χ2(nx), where nx denotes the dimensions of the state vector. If the Mahalanobis distance is within the confidence interval, put m~k(j) into the association Ωk(i)={Ωk-1(i),m~k(j)}, i≤m,j≤n. Otherwise, take target’s state estimation m~k(j) as a starting point and create a new track Ωnew={Ωnew,m~k(j)}. The tracks in this step are merged as one set: Ωk={Ωk,Ωnew}.Step 4 (pruning the tracks).
When a target exists in a certain time, the track should remain a continuous time. If the time is short, then this track is regarded as an invalid one.Merge all the tracks together:Ωk={Ωk,Ωend}. Test the number a of estimated targets in each track. If it is smaller than one predefined threshold, keep the track. Otherwise, the track is eliminated. At last, we get the final track set Ωk.


### 4.3. Numerical Interpolation

Interpolation is an important approximate method of a discrete function. For functionf(x) and n+1 points x0,x1,…,xn∈[a,b] the corresponding function values are f(x0),f(x1),…,f(xn). The interpolation is to estimate the function value at unknown points x∈[a,b]. There are many general interpolation methods, such as neighboring interpolation, linear interpolation, Lagrange interpolation, Newton interpolation, cubic spline interpolation, and Hermite interpolation. Among them, the spline interpolation uses a special piecewise, low-order polynomial spline to realize a small interpolation error [22]. Therefore, the cubic spline interpolation method is used in this paper.Assume that, in the interval[a,b], there are nodes:(31)x:a=x0<x1<⋯<xn=b,y:y0,y1,…,yn.Spline curve S(x) denotes a spline function interpolating function f(x):(32)Sx=S0x,x∈x0,x1S1x,x∈x1,x2⋮Sn-1x,x∈xn-1,xn.Equation (32) must meet some conditions: (a) on each small interval (xi-1,xi),i=1,2,…,n, S(x) is a respective cubic polynomial function; (b) its derivative S′(xi) and 2-order derivative S′′(xi) are both continuous at node xi(i=1,2,…,n-1); (c) let S(xi)=yi, i=1,2,…,n-1; then n cubic piecewise polynomials can be written as(33)Six=ai+bix-xi+cix-xi2+dix-xi3,where ai,bi,ci, and di denote the unknown coefficients. By solving the matrix or (33), the values of ai,bi,ci, and di can be obtained as follows:(34)ai=yi,bi=yi+1-yihi-hi2mi-hi6mi+1-mi,ci=mi2,di=mi+1-mi6hi,where hi=xi+1-xi denotes the length of step and mi=S′′(xi) denotes the quadratic differential value, i=1,2,…,n-1. According to (34), any corresponding value can be calculated for x at interval [a,b].For predicted statemi and its corresponding time ti of each effective track Ωk(j) obtained from the track-estimate association, we can interpolate at the discontinuous time t by using the following cubic spline interpolation:(35)x:a=t0<t1<⋯<tn=b,y:m0,m1,…,mn,where a is the starting time of a track and b is the ending time of a track; n≤(b-a+1).Afterai, bi, ci, and di are available, the jth track can be estimated as follows:(36)Sijx=aij+bijx-xi+cijx-xi2+dijx-xi3,i=1,2,…,n-1.In terms of a complete track, if the value at ti in the interval [a,b] is missing, it can be computed according to the specific expression (36). After the numerical interpolation is carried on for all tracks, each track becomes more complete; that is, the undetected targets are compensated. Thus, this is a way to solve the underestimate problem.


## 4.1. Implementation of GM-PHD

Assume that the state model of targets and the measurement model of sensors are linear Gaussian; that is,(17)pxk∣k-1xk-1=Nxk∣k-1;Fk-1xk-1,Qk-1,pzk∣xk=Nzk;Hkxk,Rk,where Fk-1 is the state transition matrix, Qk-1 is the covariance matrix of process noise, Hk is the measurement matrix of the sensors, and Rk is the covariance matrix of measurement noise. The PHD of birth targets and spawned targets are written in the following Gaussian mixture forms: (18)γkx=∑i=1Jγ,kwγ,kiNx;mγ,ki,Pγ,ki,(19)βk∣k-1x∣ξ=∑j=1Jβ,kwβ,kjNx;Fβ,k-1jξ+dβ,k-1j,Qβ,k-1j,where Jγ,k, wγ,k(i), mγ,k(i), and Pγ,k(i) are the numbers of Gaussian components, weights, mean, and covariance of birth targets, and Jβ,k, wβ,k(j), Fβ,k-1(j), dβ,k-1(j), Fβ,k-1(j)ξ+dβ,k-1(j), and Qβ,k-1(j) are the numbers of Gaussian components, weights, state transition matrix, noise disturbing items, mean, and covariance of spawned targets.Based on the posterior intensity at timek-1 in (14), the posterior intensity at time k can also be predicted by(20)vk∣k-1x=vS,k∣k-1x+vβ,k∣k-1x+γkx,where γk(x) is given in (18) and vS,k∣k-1(x) is the intensity of the survival target which can be expressed as(21)vS,k∣k-1x=pS,k∑j=1Jk-1wk-1jNx;mS,k∣k-1j,PS,k∣k-1j,where(22)mS,k∣k-1j=Fk-1mk-1j,PS,k∣k-1j=Qk-1+Fk-1Pk-1jFk-1T.The intensity of the spawned target, vβ,k∣k-1(x), is(23)vβ,k∣k-1x=∑j=1Jk-1∑i=1Jβ,kwk-1jwβ,kiNx;mβ,k∣k-1j,i,Pβ,k∣k-1j,i,mβ,k∣k-1j,i=Fβ,k-1imk-1j+dβ,k-1i,Pβ,k∣k-1j,i=Qβ,k-1i+Fβ,k-1iPk-1jFβ,k-1iT.According to (20), the predicted PHD can be written in a Gaussian mixture form:(24)vk∣k-1x=∑i=1Jk∣k-1wk∣k-1iNx;mk∣k-1i,Pk∣k-1i.Then, the posterior intensity function at time k also follows a Gaussian mixture form:(25)vkx=1-pD,kvk∣k-1x+∑z∈ZkvD,kx;z,where z∈Zk are the measurements obtained at time k and(26)vD,kx;z=∑j=1Jk∣k-1pD,kwk∣k-1jqkjzκkz+pD,k∑i=1Jk∣k-1wk∣k-1iqkiz·Nx;mk∣kjz,Pk∣kj,qkjz=Nz;Hkmk∣k-1j,HkPk∣k-1jHkT+Rk,mk∣kjz=mk∣k-1j+Kkjz-Hkmk∣k-1j,Pk∣kj=I-KkjHkPk∣k-1j,Kkj=Pk∣k-1jHkTHkPk∣k-1jHkT+Rk-1.A large number of the Gaussian components will bring a heavy burden into the computation. Usually, a merging and pruning stage is needed to deal with this problem [8]. In order to find targets in the Gaussian components, we set a weight threshold wTh (usually chosen as 0.5). If the weight is larger than wTh, the corresponding Gaussian component is regarded as a target; that is,(27)X^k=m~ki:w~ki>wTh,i=1,…,Jk.


## 4.2. Track-Estimate Association

By associating a single target track at the previous time with the state of the target at the current time, the overestimation of the target number in the GM-PHD can be restrained. First, RFSX~k at the current time is obtained according to the conational GM-PHD filter, and the track set Ωk-1 at the previous time is already generated, and then track-estimate association can be done. Due to the randomness of clutter, the error state generated by the clutter also has randomness, so we can eliminate the error of state through the existing track.Figure1 is an illustration of the track-estimate association from the time k-1 to k. In Figure 1, m~k-1(1),m~k-1(2),…,m~k-1(Lk-1) are the end nodes of tracks in Ωk-1 formed at time k-1. The end node is the same as a state vector m~k-1(i)=[xk-1(i),x˙k-1(i),yk-1(i),y˙k-1(i)]T where (xk-1(i),yk-1(i)) denote the location of the ith target and (x˙k-1(i),y˙k-1(i)) denote the velocity of the ith target. (tk-11,tk-12,…,tk-1Lk-1) denotes the end time of tracks at time k-1. It should be pointed out that the end time is not always equal to k-1 because some tracks may not be associated with the estimates at previous time. Δtki denotes a time interval from the end time of the ith track to time k; that is, Δtki=k-tk-1i, i=1,2,…,m. Δtki is the interruption time of the ith track. m~k∣k-1(1),m~k∣k-1(2),…,m~k∣k-1(l) is the predicting vector corresponding to the end time state of each track at time k-1, and m~k∣k-1(i)=[xk∣k-1(i),x˙k∣k-1(i),yk∣k-1(i),y˙k∣k-1(i)]T. If targets move with a constant speed, the prediction of the ith track can be calculated by(28)xk∣k-1i=xk-1i+Δtki·x˙k-1i,yk∣k-1i=yk-1i+Δtki·y˙k-1i,x˙k∣k-1i=x˙k-1i,y˙k∣k-1i=y˙k-1i.Figure 1
Track-estimate association.Let the prediction set after alignment beB={m~k∣k-1(1),m~k∣k-1(2),…,m~k∣k-1(l)} and let a set of state estimations obtained via GM-PHD at time k be C={m~k(1),m~k(2),…,m~k(Nk)}. Furthermore, assume that the permissible maximum interruption steps of a track is S, the minimum estimated number is anum, the distance threshold is dth, and temporary track queue Ωk=⌀.Step 1.
Judge whether a track is associated with a target or not. For a target, track should be continuous strictly; however, because of the detection ratepD<1, some targets are missed in detection. If Δtki>S, the association is not needed, and this track is terminated and pushed into the terminated track queue Ωend. If Δtki≤S, state estimation m~k(i) at time k should be associated with Ωk-1.Step 2.
Compute the distances among elements between setB and set C, and the results are used to construct a matrix D. At this time, the data association is done via Hungary algorithm so that the whole distance cost will be smallest. In the matrix of the costs D, its element Dij denotes the location distance between the ith target’s state m~k(i) in set C and the jth predicted state m~k∣k-1(j) in set B. It belongs to the balance assignment problem when l=Nk. Otherwise, it belongs to the unbalanced assignment problem. For an unbalanced assignment problem, it can be converted to a balance assignment problem via the importing dummy variable (such as zero or infinite) and turning the matrix into a square matrix. So we only need to consider the balance assignment. In this case, the matrix of the costs D is(29)D=DijNk×Nk=D11D12⋯D1NkD21D22⋯D2Nk⋮⋮⋮DNk1DNk2⋯DNkNk.Suppose H is an incidence matrix, where Hij indicates the association between the ith target’s state m~k(i) and the jth predicted state m~k∣k-1(j). When Hij=1, the two targets are associated with each other; when Hij=0, they are not associated. An associated matrix can be obtained by solving the following optimization problem:(30)zmin=minHij⁡∑i=1n∑j=1nDijHij,where zmin denotes the minimum overall benefit.The solution of (30) may be achieved through the Hungary algorithm [19, 20], which is described briefly as below. First, let each row of the matrix of the costs D be minus the smallest element in the row and get the result matrix D′ in which each row has at least one element equal to 0 and all elements satisfy Dij′≥0. Second, the same process is used to deal with columns of matrix D′. We get a new matrix of the costs from the original, in which there is at least one 0 element in each column and each row. Repeat the above two steps until there is only one nonzero element in each column and each row, which is the optimal solution of the assignment problem.We can obtain the matching matrixu according to the final matrix of the costs D′′. For the size of l and Nk, there are three kinds of cases for the number of target’s states and predicted states: l<Nk, l=Nk, and l>Nk. When l=Nk, they can be directly matched. When l<Nk, they can be matched because the number of state estimations m~k(i) matches the dummy variable. Here, these target’s states can be regarded as starting points, and new tracks Ωnew={m~k(i)} can be created. When l>Nk, they can be matched according to the number of predicted states m~k∣k-1(j). In this case, the track is interrupted if the predicted state is matched with a dummy variable according to the cost matrix. This track is regarded as an effective track: Ωk={Ωk,Ωk-1j}.Step 3.
Calculate the Mahalanobis distance between matched elements in setB and set C [21]. d(m~k∣k-1(i),m~k(j)) denotes the square of the distance judging whether its value is located in the corresponding confidence region of χ2 distribution or not.The Mahalanobis distance satisfiesd(m~k∣k-1(i),m~k(j))≤2χ2(nx), where nx denotes the dimensions of the state vector. If the Mahalanobis distance is within the confidence interval, put m~k(j) into the association Ωk(i)={Ωk-1(i),m~k(j)}, i≤m,j≤n. Otherwise, take target’s state estimation m~k(j) as a starting point and create a new track Ωnew={Ωnew,m~k(j)}. The tracks in this step are merged as one set: Ωk={Ωk,Ωnew}.Step 4 (pruning the tracks).
When a target exists in a certain time, the track should remain a continuous time. If the time is short, then this track is regarded as an invalid one.Merge all the tracks together:Ωk={Ωk,Ωend}. Test the number a of estimated targets in each track. If it is smaller than one predefined threshold, keep the track. Otherwise, the track is eliminated. At last, we get the final track set Ωk.


## 4.3. Numerical Interpolation

Interpolation is an important approximate method of a discrete function. For functionf(x) and n+1 points x0,x1,…,xn∈[a,b] the corresponding function values are f(x0),f(x1),…,f(xn). The interpolation is to estimate the function value at unknown points x∈[a,b]. There are many general interpolation methods, such as neighboring interpolation, linear interpolation, Lagrange interpolation, Newton interpolation, cubic spline interpolation, and Hermite interpolation. Among them, the spline interpolation uses a special piecewise, low-order polynomial spline to realize a small interpolation error [22]. Therefore, the cubic spline interpolation method is used in this paper.Assume that, in the interval[a,b], there are nodes:(31)x:a=x0<x1<⋯<xn=b,y:y0,y1,…,yn.Spline curve S(x) denotes a spline function interpolating function f(x):(32)Sx=S0x,x∈x0,x1S1x,x∈x1,x2⋮Sn-1x,x∈xn-1,xn.Equation (32) must meet some conditions: (a) on each small interval (xi-1,xi),i=1,2,…,n, S(x) is a respective cubic polynomial function; (b) its derivative S′(xi) and 2-order derivative S′′(xi) are both continuous at node xi(i=1,2,…,n-1); (c) let S(xi)=yi, i=1,2,…,n-1; then n cubic piecewise polynomials can be written as(33)Six=ai+bix-xi+cix-xi2+dix-xi3,where ai,bi,ci, and di denote the unknown coefficients. By solving the matrix or (33), the values of ai,bi,ci, and di can be obtained as follows:(34)ai=yi,bi=yi+1-yihi-hi2mi-hi6mi+1-mi,ci=mi2,di=mi+1-mi6hi,where hi=xi+1-xi denotes the length of step and mi=S′′(xi) denotes the quadratic differential value, i=1,2,…,n-1. According to (34), any corresponding value can be calculated for x at interval [a,b].For predicted statemi and its corresponding time ti of each effective track Ωk(j) obtained from the track-estimate association, we can interpolate at the discontinuous time t by using the following cubic spline interpolation:(35)x:a=t0<t1<⋯<tn=b,y:m0,m1,…,mn,where a is the starting time of a track and b is the ending time of a track; n≤(b-a+1).Afterai, bi, ci, and di are available, the jth track can be estimated as follows:(36)Sijx=aij+bijx-xi+cijx-xi2+dijx-xi3,i=1,2,…,n-1.In terms of a complete track, if the value at ti in the interval [a,b] is missing, it can be computed according to the specific expression (36). After the numerical interpolation is carried on for all tracks, each track becomes more complete; that is, the undetected targets are compensated. Thus, this is a way to solve the underestimate problem.


## 5. Simulation

Suppose there are four targets appearing at different times in a monitoring region[-100,100]m×[-100,100]m. The initial target states are x11=[0,2.6,0,-1.2]T, x12=[0,0,0,0]T, x13=[-10,0,0,0]T, and x14=[0,0,-5,0]T. Assume that the state vector is xki=[xki,x˙ki,yki,y˙ki]T, where (xki,yki) and (x˙ki,y˙ki) are the position and velocity of the ith target. The target’s motion equation is given as follows:(37)
x
k
i
=
1
T
0
0
0
1
0
0
0
0
1
T
0
0
0
1
x
k
-
1
i
+
T
2
2
0
T
0
0
T
2
2
0
T
ω
k
i
,where T is the sampling period (T=1 s) and ωki is Gaussian white noise whose covariance is σω2=diag⁡([0.5,0.1]). The measurement equation is(38)zki=10000010xki+υki,where υki is Gaussian white noise whose mean is zero and covariance is συ2=diag⁡([0.5,0.5]). Other experimental parameters are set as follows: the survival probability of each target Ps=0.99, the OSPA distance parameters c=70 and p=2, the pruning threshold Tprun=10-5, the merging threshold Umerg=4, the maximum tolerating number of Gaussian components Jmax⁡=100, the threshold of weight wth=0.5, the maximum interrupted steps of each track S=3, the minimum permitted number of track estimates anum=3, and the distance threshold dth equals the 90% confidence interval of χ2 distribution.


### 5.1. Performance Verification for GM-PHD-TANI

In this simulation, the detection rate is set to bePD=0.85; the number of clutters in the observation area is uniformly distributed and obeys the Poisson distribution with a mean of 25. The real and estimated trajectories via GM-PHD are illustrated in Figure 2. There are some obvious error estimations, such as points (−23.19, 5.959), (−19.05, 0.4054), and (−14.7, −3.663). This kind of errors is the overestimation problem. It actually comes from the clutters because these locations deviate from the real trajectories in a large degree. On the other hand, some estimations are leakage according to the real trajectories, and some potential intermittent points are present in the intervals from point (8.713, 31.31) to (17.65, 32.05) and from point (0, 0) to (4.903, −2.088) and so forth. The numbers of estimated targets in these time intervals are leakage, which is the underestimation error due to the low detection rate. The real and estimated trajectories via the GM-PHD-TANI filter are shown in Figure 3. It can be seen that the overestimation errors caused by clutters in Figure 2 have been overcome after the track-estimate association, and some of the leakage estimations caused by the low detection rate have been compensated effectively as well.Figure 2
The real and the estimated target trajectories of the GM-PHD.Figure 3
The real and the estimated trajectories of the GM-PHD-TANI.Figure4 shows the true target number and the target numbers estimated by the GM-PHD and GM-PHD-TANI filters. The number of estimated targets via the GM-PHD after the track-estimate association becomes lower because of the removal of some false alarms. Moreover, the number of estimated targets via the GM-PHD after the numerical interpolation becomes larger because of the compensation for undetected targets. Overall, the new processing method (GM-PHD-TANI) can improve the tracking accuracy. The OSPA distances of the GM-PHD and GM-PHD-TANI filters are given in Figure 5. It is easy to know that the GM-PHD-TANI has lowered errors (OSPA distances) as opposed to the GM-PHD.Figure 4
The numbers of the real and the estimated targets.Figure 5
OSPA distance.The same scenario as given above is repeated independently for 50 times and the average computational time is calculated. The computing times for GM-PHD and GM-PHD-TANI are 1.3053 s and 1.5284 s, respectively. GM-PHD-TANI adds only 0.2231 s to the computations of track-estimate association and numerical interpolation. It should be noted that there may be delays in the stage of numerical interpolation, and the length of the delayed time depends on the number of points interpolated.


### 5.2. Influence of the Clutter Rate and the Detection Rate on the GM-PHD-TANI

Let the detection rate bePd=0.9. The simulation was performed 100 times by using clutter numbers varying from 5 to 50. The average OSPA distances of the GM-PHD and the GM-PHD-TANI filters were calculated and displayed in Figure 6. As the average clutter number increases, the OSPA distance using GM-PHD increases much more sharply than that using GM-PHD-TANI. Hence, the GM-PHD-TANI filter is more robust than the GM-PHD when dealing with high clutter rates.Figure 6
OSPA distance using different clutter rates.Given that the average clutter number is 10, another simulation is repeated 100 times by varying the detection rate from 0.75 to 1. The corresponding OSPA distance curves are plotted in Figure7. The OSPA distances of the two algorithms have declining trends as the detection rate rises. However, the OSPA distances of the GM-PHD-TANI filter are much lower than those of the GM-PHD at all the detection rates.Figure 7
OSPA distance using different detection rates.


## 5.1. Performance Verification for GM-PHD-TANI

In this simulation, the detection rate is set to bePD=0.85; the number of clutters in the observation area is uniformly distributed and obeys the Poisson distribution with a mean of 25. The real and estimated trajectories via GM-PHD are illustrated in Figure 2. There are some obvious error estimations, such as points (−23.19, 5.959), (−19.05, 0.4054), and (−14.7, −3.663). This kind of errors is the overestimation problem. It actually comes from the clutters because these locations deviate from the real trajectories in a large degree. On the other hand, some estimations are leakage according to the real trajectories, and some potential intermittent points are present in the intervals from point (8.713, 31.31) to (17.65, 32.05) and from point (0, 0) to (4.903, −2.088) and so forth. The numbers of estimated targets in these time intervals are leakage, which is the underestimation error due to the low detection rate. The real and estimated trajectories via the GM-PHD-TANI filter are shown in Figure 3. It can be seen that the overestimation errors caused by clutters in Figure 2 have been overcome after the track-estimate association, and some of the leakage estimations caused by the low detection rate have been compensated effectively as well.Figure 2
The real and the estimated target trajectories of the GM-PHD.Figure 3
The real and the estimated trajectories of the GM-PHD-TANI.Figure4 shows the true target number and the target numbers estimated by the GM-PHD and GM-PHD-TANI filters. The number of estimated targets via the GM-PHD after the track-estimate association becomes lower because of the removal of some false alarms. Moreover, the number of estimated targets via the GM-PHD after the numerical interpolation becomes larger because of the compensation for undetected targets. Overall, the new processing method (GM-PHD-TANI) can improve the tracking accuracy. The OSPA distances of the GM-PHD and GM-PHD-TANI filters are given in Figure 5. It is easy to know that the GM-PHD-TANI has lowered errors (OSPA distances) as opposed to the GM-PHD.Figure 4
The numbers of the real and the estimated targets.Figure 5
OSPA distance.The same scenario as given above is repeated independently for 50 times and the average computational time is calculated. The computing times for GM-PHD and GM-PHD-TANI are 1.3053 s and 1.5284 s, respectively. GM-PHD-TANI adds only 0.2231 s to the computations of track-estimate association and numerical interpolation. It should be noted that there may be delays in the stage of numerical interpolation, and the length of the delayed time depends on the number of points interpolated.


## 5.2. Influence of the Clutter Rate and the Detection Rate on the GM-PHD-TANI

Let the detection rate bePd=0.9. The simulation was performed 100 times by using clutter numbers varying from 5 to 50. The average OSPA distances of the GM-PHD and the GM-PHD-TANI filters were calculated and displayed in Figure 6. As the average clutter number increases, the OSPA distance using GM-PHD increases much more sharply than that using GM-PHD-TANI. Hence, the GM-PHD-TANI filter is more robust than the GM-PHD when dealing with high clutter rates.Figure 6
OSPA distance using different clutter rates.Given that the average clutter number is 10, another simulation is repeated 100 times by varying the detection rate from 0.75 to 1. The corresponding OSPA distance curves are plotted in Figure7. The OSPA distances of the two algorithms have declining trends as the detection rate rises. However, the OSPA distances of the GM-PHD-TANI filter are much lower than those of the GM-PHD at all the detection rates.Figure 7
OSPA distance using different detection rates.


## 6. Conclusions

In multitarget tracking problems, the clutter rate and the detection rate can affect the tracking results greatly. The GM-PHD filter can exhibit a good tracking performance only when the clutter rate is low and the detection rate is high. This paper improves the GM-PHD filter via track-estimate association and numerical interpolation. The track-estimate association is used to remove most clutters to lessen the overestimate problem, while the numerical interpolation is used to compensate undetected targets so that the underestimate problem can be alleviated. The experimental results show that the improved algorithm has a higher tracking accuracy for multitargets than the GM-PHD filter. However, the real-time performance of this improved algorithm is lowered slightly because of the computation costs for the track-estimate association and the numerical interpolation.



---

*Source: 629023-2015-07-12.xml*