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

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


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], 2 Mathematical Problems in Engineering 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: Section 2 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.

Problem Description
2.1.Single Target Bayesian Filter.Assume that, in the state space R   ,   is the state vector and   is the dimension of the state vector.The state transition equation is where  −1 ∼ (0,   ) is a process noise with a Gaussian distribution whose mean is zero and covariance is   .
In the measurement space R   , the measurement vector is   and the dimension of the measurement vector is   .The measurement equation is where   ∼ (0,   ) is the measurement noise with a Gaussian distribution whose mean is zero and covariance is   .The probability density function   (  |  1: ) of the state   , decided by measurements  1: = ( 1 , . . .,   ), is the posterior probability density function.According to the prior distribution,  0 ( 0 ), the posterior probability density function can be deduced via the Bayes theory; that is, (⋅ |  1: ) contains the relevant information about the state   at time .The estimated information can be obtained by minimizing the mean square error of state estimation or by maximizing the posterior probability   (  |  1: ).

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 time ,   = { ,1 , . . .,  ,  } is the multitarget state set,   is the number of targets,   = { ,1 , . . .,  ,  } is the measurement set, and   is the number of measurements.Among   and   ,  , denotes the th target's state, and  , denotes the th target's measurement: where  ⊆ R   and  ⊆ R   denote the state space and the measurement space, respectively.The measurement set includes clutters.Assume the probability density of clutter is   (⋅), its number obeys Poisson distribution whose mean is   , the state of a target continues from time  − 1 to  with probability   ( −1 ) and disappears with probability 1 −   ( −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 ( −1 ).The new birth target set at time  is expressed as Γ  , and the spawned target derived from the survival target is denoted as ().The state model of multiple targets is The received measurement set is the mixture of the target RFS, (  ), with the clutter RFS,   ; that is, Let  1: = { 1 , . . .,   }.The multitarget posterior density function can be obtained recursively from the optimal multitarget Bayesian filtering through the following formulas: where   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.

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  is , and the first-order moment is denoted as V, which is a negative function called the intensity function.For a region  in space  ( ⊆ ), denotes the posterior intensity of the multitargets at time  − 1, then the predicted intensity at time  is where   (  ) is the intensity of the birth target, is the state transition probability density of a single target, and  ,|−1 ( −1 ) is the survival probability of a target from time  − 1 to .
The target number at one prediction step is where  Γ |−1 = ∫   () is the number of birth targets; the number of survival targets; and The PHD is updated as follows when the set of measurements   is received: where   (  ) is the intensity of clutter,   (  |   ) is the likelihood function of a single measurement, and  , (  ) is the probability of detection.Correspondingly, the target number is updated as follows: The Gaussian mixture form of a PHD is formed by Gaussian components, that is, { ()   ,  ()  ,  ()  }   =1 .In the linear Gaussian system, the posterior intensity at time  can be constructed with Gaussian components in a linear combination with different weights: where   is the number of Gaussian components, and  ()  ,  ()   , and  ()  , respectively, denote weight, state, and covariance of the th Gaussian component.
Assume that the posterior intensity of multitargets at time  − 1 is expressed in the following Gaussian mixture form: The PHD of multitargets at time  can be also expressed in a Gaussian mixture form: Then, the posterior PHD of multitargets at time  can be achieved as follows: where  −1 is the state transition matrix,  −1 is the covariance matrix of process noise,   is the measurement matrix of the sensors, and   is the covariance matrix of measurement noise.The PHD of birth targets and spawned targets are written in the following Gaussian mixture forms: where  , ,  () , ,  () , , and  () , are the numbers of Gaussian components, weights, mean, and covariance of birth targets, and  , ,  Based on the posterior intensity at time  − 1 in ( 14), the posterior intensity at time  can also be predicted by where   () is given in (18) and V ,|−1 () is the intensity of the survival target which can be expressed as where The intensity of the spawned target, V ,|−1 (), is ( According to (20), the predicted PHD can be written in a Gaussian mixture form: Then, the posterior intensity function at time  also follows a Gaussian mixture form: where  ∈   are the measurements obtained at time  and 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  Th (usually chosen as 0.5).If the weight is larger than  Th , the corresponding Gaussian component is regarded as a target; that is, Figure 1: Track-estimate association.
Figure 1 is an illustration of the track-estimate association from the time  − 1 to .In Figure 1, m |−1 ] T .If targets move with a constant speed, the prediction of the th track can be calculated by Let the prediction set after alignment be  = {m (1)  |−1 , m(2) |−1 , . . ., m() |−1 } and let a set of state estimations obtained via GM-PHD at time  be  = {m (1)   , m(2)  , . . ., m(  )  }.Furthermore, assume that the permissible maximum interruption steps of a track is , the minimum estimated number is  num , the distance threshold is  th , and temporary track queue Ω  = ⌀.
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 rate   < 1, some targets are missed in detection.If Δ   > , the association is not needed, and this track is terminated and pushed into the terminated track queue Ω end .If Δ   ≤ , state estimation m()  at time  should be associated with Ω −1 .
Step 2. Compute the distances among elements between set  and set , and the results are used to construct a matrix .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 , its element   denotes the location distance between the th target's state m() in set  and the th predicted state m() |−1 in set .It belongs to the balance assignment problem when  =   .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  is Suppose  is an incidence matrix, where   indicates the association between the th target's state m()  and the th predicted state m() |−1 .When   = 1, the two targets are associated with each other; when   = 0, they are not associated.An associated matrix can be obtained by solving the following optimization problem: where  min 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  be minus the smallest element in the row and get the result matrix   in which each row has at least one element equal to 0 and all elements satisfy    ≥ 0. Second, the same process is used to deal with columns of matrix   .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 matrix  according to the final matrix of the costs   .For the size of  and   , there are three kinds of cases for the number of target's states and predicted states:  <   ,  =   , and  >   .When  =   , they can be directly matched.When  <   , they can be matched because the number of state estimations m()  matches the dummy variable.Here, these target's states can be regarded as starting points, and new tracks Ω new = {m ()   } can be created.When  >   , they can be matched according to the number of predicted states m() |−1 .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:

Mathematical Problems in Engineering
Step 3. Calculate the Mahalanobis distance between matched elements in set  and set  [21].( m() |−1 , m()  ) 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 satisfies 𝑑( m(𝑖)
|−1 , m()  ) ≤ 2 2 (  ), where   denotes the dimensions of the state vector.If the Mahalanobis distance is within the confidence interval, put m()  into the association Ω ()  = {Ω 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: Ω  = {Ω  , Ω end }.Test the number  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 Ω  .

Numerical Interpolation.
Interpolation is an important approximate method of a discrete function.For function () and  + 1 points  0 ,  1 , . . .,   ∈ [, ] the corresponding function values are ( 0 ), ( 1 ), . . ., (  ).The interpolation is to estimate the function value at unknown points  ∈ [, ].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.
According to (34), any corresponding value can be calculated for  at interval [, ].
For predicted state   and its corresponding time   of each effective track Ω ()  obtained from the track-estimate association, we can interpolate at the discontinuous time  by using the following cubic spline interpolation: where  is the starting time of a track and  is the ending time of a track;  ≤ ( −  + 1).
After   ,   ,   , and   are available, the th track can be estimated as follows: In terms of a complete track, if the value at   in the interval [, ] 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.where  is the sampling period ( = 1 s) and    is Gaussian white noise whose covariance is  2  = diag([0.5,0.1]).The measurement equation is

Simulation
where    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   = 0.99, the OSPA distance parameters  = 70 and  = 2, the pruning threshold  prun = 10 −5 , the merging threshold  merg = 4, the maximum tolerating number of Gaussian components  max = 100, the threshold of weight  th = 0.5, the maximum interrupted steps of each track  = 3, the minimum permitted number of track estimates  num = 3, and the distance threshold  th equals the 90% confidence interval of  2 distribution.

Performance Verification for GM-PHD-TANI.
In this simulation, the detection rate is set to be   = 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 4 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.
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 trackestimate 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.

Influence of the Clutter Rate and the Detection Rate on the GM-PHD-TANI.
Let the detection rate be   = 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.
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 Figure 7.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.

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

Figure 2 :
Figure 2: The real and the estimated target trajectories of the GM-PHD.

Figure 3 :Figure 4 :
Figure 3: The real and the estimated trajectories of the GM-PHD-TANI.
)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, RFS X at the current time is obtained according to the conational GM-PHD filter, and the track set Ω −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.
) denotes the end time of tracks at time  − 1.It should be pointed out that the end time is not always equal to  − 1 because some tracks may not be associated with the estimates at previous time.Δ denotes a time interval from the end time of the th track to time ; that is, Δ = −  −1 ,  = 1, 2, ..., .Δ   is the interruption time of the th track.m(1)|−1, m(2)|−1 , . . ., m() |−1 is the predicting vector corresponding to the end time state of each track at time −1, and m(