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.
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: 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.
2. Problem Description
2.1. Single Target Bayesian Filter
Assume that, in the state space , is the state vector and is the dimension of the state vector. The state transition equation iswhere is a process noise with a Gaussian distribution whose mean is zero and covariance is .
In the measurement space , the measurement vector is and the dimension of the measurement vector is . The measurement equation iswhere is the measurement noise with a Gaussian distribution whose mean is zero and covariance is .
The probability density function of the state , decided by measurements , is the posterior probability density function. According to the prior distribution, , the posterior probability density function can be deduced via the Bayes theory; that is,
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 .
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 time , is the multitarget state set, is the number of targets, 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 and 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 to with probability and disappears with probability , 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 . 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 . 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.
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 is , and the first-order moment is denoted as , which is a negative function called the intensity function. For a region in space ,
If denotes the posterior intensity of the multitargets at time , then the predicted intensity at time iswhere is the intensity of the birth target, is the intensity of the spawned target, is the state transition probability density of a single target, and is the survival probability of a target from time to .
The target number at one prediction step iswhere is the number of birth targets; is the number of survival targets; and is the number of spawned targets.
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, . 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 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:
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,where is the state transition matrix, 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 , , , , , and 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 time in (14), the posterior intensity at time can also be predicted bywhere is given in (18) and is the intensity of the survival target which can be expressed aswhereThe intensity of the spawned target, , 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 (usually chosen as 0.5). If the weight is larger than , the corresponding Gaussian component is regarded as a target; that is,
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 at the current time is obtained according to the conational GM-PHD filter, and the track set 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.
Figure 1 is an illustration of the track-estimate association from the time to . In Figure 1, are the end nodes of tracks in formed at time . The end node is the same as a state vector where denote the location of the th target and denote the velocity of the th target. denotes the end time of tracks at time . It should be pointed out that the end time is not always equal to 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, , . is the interruption time of the th track. is the predicting vector corresponding to the end time state of each track at time , and . If targets move with a constant speed, the prediction of the th track can be calculated by
Let the prediction set after alignment be and let a set of state estimations obtained via GM-PHD at time be . Furthermore, assume that the permissible maximum interruption steps of a track is , the minimum estimated number is , the distance threshold is , 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 , some targets are missed in detection. If , the association is not needed, and this track is terminated and pushed into the terminated track queue . If , state estimation at time should be associated with .
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 in set and the th predicted state 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 isSuppose is an incidence matrix, where indicates the association between the th target’s state and the th predicted state . When , the two targets are associated with each other; when , they are not associated. An associated matrix can be obtained by solving the following optimization problem:where 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 . 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 matches the dummy variable. Here, these target’s states can be regarded as starting points, and new tracks can be created. When , they can be matched according to the number of predicted states . 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: .
Step 3. Calculate the Mahalanobis distance between matched elements in set and set [21]. denotes the square of the distance judging whether its value is located in the corresponding confidence region of distribution or not.
The Mahalanobis distance satisfies , where denotes the dimensions of the state vector. If the Mahalanobis distance is within the confidence interval, put into the association , . Otherwise, take target’s state estimation as a starting point and create a new track . The tracks in this step are merged as one set: .
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: . 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 .
4.3. Numerical Interpolation
Interpolation is an important approximate method of a discrete function. For function and points the corresponding function values are . 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.
Assume that, in the interval , there are nodes:Spline curve denotes a spline function interpolating function :Equation (32) must meet some conditions: (a) on each small interval , is a respective cubic polynomial function; (b) its derivative and 2-order derivative are both continuous at node ; (c) let , ; then cubic piecewise polynomials can be written aswhere , and denote the unknown coefficients. By solving the matrix or (33), the values of , and can be obtained as follows:where denotes the length of step and denotes the quadratic differential value, . 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; .
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.
5. Simulation
Suppose there are four targets appearing at different times in a monitoring region . The initial target states are , , , and . Assume that the state vector is , where and are the position and velocity of the th target. The target’s motion equation is given as follows:where is the sampling period ( s) and is Gaussian white noise whose covariance is . The measurement equation iswhere is Gaussian white noise whose mean is zero and covariance is . Other experimental parameters are set as follows: the survival probability of each target , the OSPA distance parameters and , the pruning threshold , the merging threshold , the maximum tolerating number of Gaussian components , the threshold of weight , the maximum interrupted steps of each track , the minimum permitted number of track estimates , and the distance threshold equals the 90% confidence interval of distribution.
5.1. Performance Verification for GM-PHD-TANI
In this simulation, the detection rate is set to be ; 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 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 be . 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.
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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by National Natural Science Foundation of China (61201118), the China Postdoctoral Science Foundation (2013M532020), the Scientific Research Program funded by Shaanxi Provincial Education Department (14JK1304), and the Xi’an Polytechnic University Young Scholar Backbone Supporting Plan.