Anti-Clutter Gaussian Inverse Wishart PHD Filter for Extended Target Tracking

The extended target Gaussian inverse Wishart probability hypothesis density (ET-GIW-PHD) filter overestimates the number of targets under high clutter density. The reason for this is that the source of measurements cannot be determined correctly if only the number of measurements is used. To address this problem, we proposed an anti-clutter filter with hypothesis testing, we take into account the number of measurements in cells, the target state and spatial distribution of clutter to decide whether the measurements in cell are clutter. Specifically, the hypothesis testing method is adopted to determine the origination of the measurements. Then, the likelihood functions of targets and clutter are deduced based on the information mentioned above, resulting in the likelihood ratio test statistic. Next, the likelihood ratio test statistic is proved to be subject to a chi-square distribution and a threshold corresponding to the confidence coefficient is introduced and the measurements below this threshold are considered as clutter. Then the correction step of ET-GIW-PHD is revised based on hypothesis testing results. Extensive experiments have demonstrated the significant performance improvement of our proposed method.


Introduction
Extended target tracking (ETT) draws lots of attention in recent years because of its wide range of applications in traffic control [1], autonomous driving [2][3][4], person tracking [5,6] and etc. [7][8][9][10][11]. Since one extended target generates more than one measurement per time step, its shape information can be obtained. Using this information, the kinematic state and extent of the target can be estimated simultaneously. The extent of the target including the size, shape, and orientation can be further used for target identification.
The difference between point target tracking and extended target tracking lies in the measurement model and hypotheses. Point target generates at most one measurement per time step, while the extended target generates multiple measurements. Many algorithms were proposed to track point target based on point target hypothesis, such as probability hypothesis density (PHD) filter [12] and cardinalized PHD (CPHD) filter [13]. Since the extended target violates one measurement hypothesis, the number of targets will be overestimated if point target tracking algorithms are directly used for tracking extended targets. To address this problem, Mahler proposed an extended target tracking algorithm based on the inhomogeneous Poisson Point Process model (PPP model) [14] in random finite sets (RFSs) frame, namely extended target probability hypothesis density (ET-PHD) [15], Jiang et al. [16] proposed a novel time-matching ET-PHD filter, a Gaussian mixture implementation of the ET-PHD, called the extended target Gaussian inverse Wishart probability hypothesis density (ET-GIW-PHD) filter, has been presented in [17]. However, ET-PHD only estimates the kinematic state of the target (such as position, velocity) and does not estimate the extent of the target. Therefore, this method cannot extract the shape of the target. Nevertheless, the estimation of the target extent is important because it can be used to classify target and improve tracking accuracy [18][19][20].
Measurement partition is an important step in ETT. In ETT, measurements are partitioned into several non-empty subsets, each subset contains measurements that are all from the same source, either a single target or a clutter source, the subset is defined as cell. In ETT, the increase of measurements gives rise to the quick increase of the set partitions, thus the partition algorithm should be designed to achieve tractable computational complexity. Distance partition [17] is the most widely used method. Modified Bayesian adaptive resonance theory (MB-ART) [21] can also achieve good performance. For more details about other partition algorithms, please see [22][23][24].
One of the most important works in extended target tracking is how to model the target extent. To address this problem, the stick model is used for bicycle and pedestrian tracking [25,26]. The object extension is represented by a symmetric positive definite (SPD) random matrix [27], namely a random matrix (RM) model. Feldmann et al. [28] adapted the RM model for the case when the sensor error cannot be ignored. Lan et al. [29] took into account time variation and distortion of target extension in RM frame. In order to handle irregular shapes, a random hypersurface model (RHM) is introduced in [30][31][32]. Gaussian Processes (GP) was used to represent the target shape and achieved good performance [33][34][35][36]. Since shape estimation is similar to curving fitting, Kaulbersch et al. [37] applied a curve fitting method for shape estimation. Granström et al. [38] proposed an extension model for specific sensor. Granström et al. [39] proposed extended target Gaussian inverse Wishart PHD (ET-GIW-PHD) filter to incorporate widely used RM model into PHD filter and approximate the estimated PHD with an unnormalized mixture of Gaussian inverse Wishart (GIW) distributions. Later, Granström et al. [40] proposed extended target Gamma Gaussian inverse Wishart PHD (ET-GGIW-PHD) filter to estimate the measurement rate and target state simultaneously. The combination of several RM model was used to model nonelliptic targets in [41,42]. As mentioned in [39], more experiments that test ET-GIW-PHD filters are needed, e.g., for data that contains more clutter than typical laser data does, this provides the general motivation for this paper. We found that the number of targets will be overestimated which degrades the final performance when severe clutters are partitioned into one cell in ET-GIW-PHD. More analyses are presented in Section 3. In this paper, we proposed an anti-clutter ET-GIW-PHD filter for better cardinality estimation performance.
The main contributions of this paper are twofold. First, the reason why ET-GIW-PHD overestimates the number of targets is discussed detailedly, and the probability of the measurement generated by clutter against different scenario parameters is presented. Second, in order to deal with the cardinality overestimation in ET-GIW-PHD, we proposed an anti-clutter ET-GIW-PHD filter which revises the correction step of ET-GIW-PHD with hypothesis testing. Hypothesis testing is introduced to determine the source of measurements in the cell, hypothesis testing results are integrated into the correction step in ET-GIW-PHD. In order to deal with the source of measurements correctly, the essential differences between the measurements of targets and clutter should be recognized. Since the variation of target state over time follows certain rules (motion model and shape transition model), the state of targets could be predicted while clutter could not. Then, the likelihood functions of targets and clutter are deduced. The likelihood functions are built based on not only the number of measurements but also the target state and spatial distribution of clutter. Since the likelihood ratio test statistic is proved to be subject to chi-square distribution, a threshold corresponding to the confidence coefficient is introduced, this threshold is used to determine the source of measurements in the cell. It worth note that the perfect sensor resolution is advocated as a theoretical hypothesis in this paper. In reality, the results in the Section 5 will be affected by the limited sensor resolution. Future work will tackle the sensor's limited resolution.
The rest of the paper is outlined as follows. Section 2 reviews the ET-GIW-PHD filter. Section 3 discusses the reason why ET-GIW-PHD overestimates the number of targets. Our anti-clutter ET-GIW-PHD is presented detailedly in Section 4. We conduct experiments in different simulation scenarios to demonstrate the effectiveness of our proposed approach in Section 5; Conclusion is drawn in Section 6.

ET-GIW-PHD Review
In ET-GIW-PHD, both predicted PHD and corrected PHD can be approximated as an unnormalized mixture of Gaussian inverse Wishart distributions. Let ζ k = {x k , X k } be the sufficient statistics of the GIW components at time which contains kinematical state x k and extension state X k which is mathematically described by a symmetric and positively definite (SPD) random matrix. The iterative formulae for ξ k are obtained in [39]. More implementation details, such as pruning and merging, can also be found in [39]. Prediction: where p s (·) is the probability of survival, p k+1|k (ξ k+1 |ξ k ) is the state transition density, D b k+1 (·) is the birth PHD, new target spawning is omitted [39]. Correction: The corrected PHD D k|k (ξ k ) can be summarized as: where p∠Z k means that the measurement sets Z k are partitioned into non-empty cells, W ∈ p means that the cell W is in the partition p. D ND k|k (ξ k ) handles the undetected target case, because D k+1|k (ξ k+1 ) is approximated as an unnormalized mixture of Gaussian inverse Wishart distributions, it is given by where J k|k−1 is the number of components of predicted PHD, w (j) k|k is the weight of GIW component. N (x; m, P) means that a vector x is subject to Gaussian distribution with mean m and covariance P, IW (X; v, V) means that a matrix is subject to inverse Wishart distribution with degree of freedom v and inverse scale matrix V. ⊗ is the Kronecker product. D D k|k (ξ k , W) handles the detected target case, which is given by can be obtained by where presents the likelihood of the jth GIW component given the measurements of the Wth cell, ω p is the weight of pth partition, p (j) D is the detection probability of jth GIW component, γ (j) is the expected number of measurements generated by jth GIW component, λ c is the mean number of clutter measurements, c k is the spatial distribution of the clutter over the surveillance volume, δ i,j is the Kronecker delta, |W| is the the number of measurements in the Wth cell, S (j,W) k|k−1 is innovation factor, Γ d (·) is the multivariate Gamma function.

Analysis of ET-GIW-PHD
In ET-GIW-PHD, the calculation of w (j,W) k|k is important. If the measurements in Wth cell are generated by clutter, w (j,W) k|k is expected to be smaller than the pruning threshold, then the corresponding component will be eliminated and the clutter will be eliminated.
In Equation (5), w (j,W) k|k contains two parts, one is the weight of the pth partition, denoted by ω p , the other is the weight of Wth cell in partition. Without loss of generality, only one partition is considered for clarity, therefore ω p = 1. Substituting Equation (6) into Equation (5), we arrive at From Equation (9), the numerator is a part of denominator, the measurements of Wth cell is used to correct each GIW component, then Λ can be given based on some prior parameters, such as p D , γ, λ c and c k (for brevity, the subscript and superscript are omitted here).
If the measurements in the cell are generated by clutter, the likelihood Λ (j,W) k of each GIW component will be very small since clutter does not obey the kinematic and extent model of target. If the number of clutter measurements in the cell is equal to one, then |W| = 1, δ |W|,1 = 1, k|k−1 will be much smaller than 1 because the likelihood Λ (j,W) k achieves a small value mentioned above and other parameters can be considered as constants, the value of w (j,W) k|k will be close to 0 and is smaller than the pruning threshold, then the corresponding component will be eliminated and the clutter is eliminated. However, if the number of clutter measurements in the cell is more than one, then |W| = 1, δ |W|,1 = 0, Equation (9) is the normalization process. Although can still take a large value. In this case, ghost targets will emerge and the number of targets will be overestimated. Further details on numerical implementation can be found in Section 5.
According to the analysis above, if the measurement in the cell is clutter, (9) should be added by 1. Otherwise, it should be added by 0 and the clutter can be eliminated. However, from Equation (9), if the cell contains only one measurement, ψ l,W is added by 1, it means that the cell contains only one measurement is considered as clutter in ET-GIW-PHD. Otherwise, it is considered as a target if the cell contains more than one measurement. In fact, this assumption can be violated under strong clutter. The criterion whether measurements in the cell are generated by clutter based on only the number of measurements can be erroneous. A simple numerical calculation is shown below to illustrate this point.
In ET-GIW-PHD, the probability of the measurements of the cell generated by clutter is obtained based on the Bayesian theorem, see Equation (10). Note that, only the number of measurement is considered in this calculation.
where Z W presents the measurements in cell, n W is the number of measurements in cell, C and T mean clutter and target respectively, P(Z W ⊂ C) and P(Z W ⊂ T) are the prior information.
The number of measurements generated by the target is subject to Poisson distribution with Poisson rate γ, the detection probability is p d , then where C n m = n! m!(m−n)! denotes the combinatorial number of the events that m out of n. Remark: p d is not equal to p D in Equation (6). p d is the probability that one measurement generated by target or clutter is detected, while p D is the probability that an extended target will generate a measurement set [15]. p D can be derived if p d is already known.
The clutter measurements are assumed to be uniformly distributed over the surveillance area, and the number of clutter is subject to Poisson distribution with Poisson rate λ c . So we have When the number of measurements is 1, the probability of the measurement in the cell generated by clutter is shown in Table 1 with different γ and λ c . In this simulation, the prior information is set to 0.5, then P(Z W ⊂ C) = P(Z W ⊂ T) = 0.5. Table 1. The probability of the measurement generated by clutter.
From Table 1 we can see that when γ = 10, λ c = 35 and p d = 0.9, P(Z W ⊂ C|n W = 1) = 0.013. Although the cell contains only one measurement, the probability of the measurement in the cell generated by clutter is close to 0. Consequently, the criterion of ET-GIW-PHD does not work well in this case. when γ = 10 and λ c = 5, P(Z W ⊂ C|n W = 1) = 1. In this case, the clutter is distinguished correctly based on the criterion of ET-GIW-PHD. In summary, the determination whether measurements are generated by clutter based on only the number of measurements can be erroneous.

Anti-Clutter ET-GIW-PHD
The difference between ET-GIW-PHD and anti-clutter ET-GIW-PHD is how to determine the source of measurements in the cell, specifically, the difference is how to obtain d W in Equation (6). Using only the number of measurements in ET-GIW-PHD to determine whether the measurements in the cell is the target or not may be erroneous. In contrast, our anti-clutter ET-GIW-PHD uses hypothesis testing to deal with this problem. The number of measurements, the kinematic state and extent state of target and clutter spatial distribution are taken into account to obtain the likelihood ratio test statistic.
There are two hypotheses: where Z W = {z 1 , z 2 , ..., z n W } is the measurements of Wth cell, n W is the number of the measurements, C and T represent clutter and target respectively. The likelihood ratio test statistic for hypothese is given by where L(Z W |T) and L(Z W |C) are the likelihood to measure the set Z W given Z W ⊂ T and Z W ⊂ C respectively and L(Z W |T) and L(Z W |C) will be presented later.
Because log(·) is monotony increase, log(η) is also the test statistic for hypothesis H 0 versus H 1 . If these measurements are generated by the target, L(Z W |T) will achieve a large value and L(Z W |C) will be small. Consequently, the statistics log(η) will grow to a large value. Using a threshold, we can distinguish between targets and clutter. Specifically, if log(η) is greater than the threshold, the measurements in the cell is considered to be generated by targets. Otherwise, these measurements are considered to be clutter. The expression of L(Z W |T) and L(Z W |C) are given below, the setting of the threshold is discussed.
If the measurements are generated by a target, different extent models lead to a different expression of L(Z W |T). In this paper, L(Z W |T) is deduced based on the model in [27], where I d is an unit matrix with d dimension, X k is the extension of target at time k, H k is the 1D observation matrix.
The clutters are assumed to be uniformly distributed over the surveillance area [27], then where β FA = λ c c k , λ c is the mean number of clutter measurements, c k is the spatial distribution of the clutter over the surveillance volume. Substitute Equations (17) and (18) into Equation (16), we have where Because z j is subject to Gaussian distribution with mean (H k ⊗ I d )x k and covariance X k , where G is subject to chi-square distribution with degree of freedom n W . In Equation (21), γ, λ c and β FA are priori known, the volume of the target extension is proportional to |X k |, the size of the target could be assumed to be unchanged, then D could be considered as a constant. The confidence coefficient is set to α and a threshold is introduced (denoted by g), suppose hypothesi H 1 is true, then where From Equation (23), the probability of log(η) < g is 1 − α. Generally, α is set to be a value close to 1 and log(η) < g is a small probability event. If log(η) < g is satisfied, hypothesi H 1 should be rejected. Finally, we have If log(η) < g, the measurements are generated by clutter, then If log η ≥ g,the measurements are generated by targets, then The pseudo-code for anti-clutter ET-GIW-PHD is illustrated in Table 2. The difference between ET-GIW-PHD and anti-clutter ET-GIW-PHD lies in correction step. Pseudo-code for anti-clutter ET-GIW-PHD filter correction is shown in Table 3, pseudo-code for other steps (prediction, prune and merge etc) can be found in [39]. Table 3. Pseudo-code for anti-clutter ET-GIW-PHD filter correction.

Simulation
In this section, the effectiveness of ET-GIW-PHD and anti-clutter ET-GIW-PHD were tested. Two scenarios with multiple targets were established. The surveillance area was set as [−1000 m, 1000 m] × [−1000 m, 1000 m], then c k is 2.5 × 10 −7 under the assumption that clutter is uniformly distributed over the surveillance area. We set totally 100 time steps and the sampling time is 1 s.
In the first scenario, four targets moving along different lines were generated: In the second scenario, two targets were born at (−1000 m, 300 m) and (−1000 m, −300 m), respectively at k = 0 (k is time step). Next, they moved close gradually and then moved in parallel before separating. The birth intensity in the second scenario is where m 1 = [−1000, 300, 25, −25], m 2 = [−1000, −300, 25,25], 100, 100]). The true trajectories of two scenario are shown in Figure 1.  The dynamic and measurement model are shown below. The target kinematic state is denoted as x = [r x , r y ,ṙ x ,ṙ y ], where r x andṙ x is the position and velocity in the x direction, likewise of y direction. The time evolution of kinematic state given by where x (j) k is the target state of jth target at time k, w k is the process noise of jth target and is the Gaussian white noise with zero mean and covariance Q (j) k is the kinematic state transition matrix of jth target, given by where t is the sampling time and Ω represents the acceleration error, t = 1 s and Ω = 5 m/s 2 in this simulation.
In this simulation, the major and minor axes are 20 m and 15 m respectively for all targets. The major axis was aligned with the direction of motion of the target and the extent of these targets remained unchanged.
The measurement model can be expressed as where z (j) k is the measurements generated by the jth target at time k, q k is the measurement noise and is the Gaussian white noise with zero mean and covariance R k , H k ⊗ I d is the observation matrix, given by where In our experiment, the confidence coefficient α of anti-clutter ET-GIW-PHD is set to 0.99. A distance partition algorithm [17] is used for both filters, a measurement partition that contains several cells can be obtained for a given distance threshold. Clutter Poisson rate λ c is set to 35, then clutter density λ c c k is 8.75 × 10 −6 (the clutter density in this paper is higher than that of related references, such as [18,21]). The expected number of measurements generated by targets γ is set to 15. The probability of survival p s and the detection probability p D are assumed to be state independent and set to 0.99 and 0.98, respectively. The probability p d is set to 0.99.
The OSPA distance is defined by where m < n, κ k = {x k , ...,x (n) k } is the estimated RFS, ∏ n is the assignment results which assign κ toκ, p means p − norm, c is the penalty cost for cardinality mismatch. In this simulation, c = 60 and p = 2.
ET-GIW-PHD and anti-clutter ET-GIW-PHD are applied to two scenarios mentioned above for performance evaluation. The trajectories generated by these two methods are presented in Figures 2 and 3. To further verify the analysis in Section 3, the calculation of w (j,W) k|k in Equation (9) at k = 40 (k is time step) in scenario 1 is shown below. The partition result at k = 40 is given firstly in Figure 4. From Figure 4 we can see that the measurements of four targets are correctly clustered, and two clutter (marked with arrows in Figure 4) are incorrectly partitioned into one cell.
k|k−1 is denoted as ψ j,W for jth GIW component in the Wth cell, then From simulation results, the number of components of predicted PHD is 14 at k = 40, then J k|k−1 = 14, ψ j,W of the clutter cell (marked with arrows in Figure 4) is obtained and shown in Table 4.
The likelihood Λ (j,W) k of each GIW component in this cell is very small since clutter does not obey the kinematic and extent model of target, therefore ψ j,W achieve small value as shown in Table 4. Because the number of measurement in this cell is two, Equation (37) is represent as Equation (38) is a normalization process, w (j,W) k|k of the clutter cell is shown in Table 5.  Although ψ j,W is small, w (j,W) k|k may achieve a large value (w (10,W) k|k = 0.99) and results in a ghost target. At k = 40, the estimated number of targets was 5 while true number is 4. That is, the number of targets was overestimated.
To test the influence of the clutter density on tracking performance, ET-GIW-PHD and anti-clutter ET-GIW-PHD were tested under different numbers of clutter modeled as Poisson distribution with Poisson rate λ c . The clutter measurements are assumed to be uniformly distributed over the surveillance area. The OSPA distance of these two filters under different Poisson rate λ c is shown in Figures 5 and 6.
As we can see from Figures 5 and 6, when λ c is small, ET-GIW-PHD achieves good performance. However, as λ c increases, the performance of ET-GIW-PHD degrades significantly. In contrast, our anti-clutter ET-GIW-PHD achieves superior performance with varying λ c , which demonstrates that anti-clutter ET-GIW-PHD is more robust to clutter than ET-GIW-PHD. The results of cardinality estimation are shown in Figures 7 and 8.   Figures 7 and 8 we can see that the cardinality estimation error of ET-GIW-PHD increases as the λ c grows. That is, ET-GIW-PHD cannot avoid the overestimation of cardinality under high clutter density. When clutter density is small, the clutter spreads apart. Thus, it is unlikely to partition more than one clutter into one cell. In the presence of severe clutter, the probability that multiple clutter being partitioned into one cell increases, and thus ET-GIW-PHD could overestimate the cardinality. However, our anti-clutter ET-GIW-PHD uses not only the number of measurement, but also target state and spatial distribution of clutter for better cardinality estimation performance. Using hypothesis testing, the measurements can be distinguished more correctly. Therefore, a better tracking performance can be achieved. Extensive experiments have demonstrated the effectiveness of anti-clutter ET-GIW-PHD.

Conclusions
In this paper, we propose an anti-clutter ET-GIW-PHD filter which revises the correction step of ET-GIW-PHD with hypothesis testing for better tracking performance under severe clutter. Our anti-clutter ET-GIW-PHD adopts a hypothesis testing method to distinguish between measurements from targets and clutter, hypothesis testing results are incorporated into the correction step. Specifically, likelihood functions are built to incorporate the number of measurements, the target state, and clutter spatial distribution in anti-clutter ET-GIW-PHD, the source of measurements in the cell is determined more correctly. Compared with ET-GIW-PHD, our method improves the cardinality estimation accuracy and achieves better tracking performance. The effectiveness of our method has been demonstrated by extensive experiments.
Author Contributions: Y.H. conceived the idea of the study and conducted the research. Y.H. and L.W. analyzed most of the data and wrote the initial draft of the paper. X.W. and W.A. contributed to finalizing this paper.
Funding: This work was supported by the National Natural Science Foundation of China (Nos.61605242).