A Measurement Set Partitioning Algorithm Based on CFSFDP for Multiple Extended Target Tracking in PHD Filter

. The extended target probability hypothesis density (ET-PHD) filter is a promising approach for multiple extended target tracking. One crucial problem of the ET-PHD filter is partitioning the measurement set. This paper proposes a partitioning algorithm based on clustering by fast search and find density peaks (CFSFDP). Firstly, we adopt CFSFDP algorithm to partition the measurement set and the field theory is introduced to determine the cutoff distance of the CFSFDP algorithm. Then, the cluster center of the CFSFDP algorithm is determined according to solved cutoff distance and measurement rate. Finally, as the CFSFDP algorithm cannot handle the case in which targets are spatially close, an improved sub-partitioning method is implemented. Simulation results show that the proposed algorithm has less computational complexity and stronger robustness than the existing algorithm without losing tracking performance.


Introduction
In target tracking, the target echo signal may be occupied in different range resolution cells. The sensor will receive multiple measurements from a target per scan and such target is referred to as an extended target [1]. In recent years, multiple extended target tracking has been hot research in the target tracking field [2], [3].
In multiple extended target tracking, the conventional algorithms, such as joint probabilistic data association (JPDA) and multiple hypothesis tracker (MHT), burden exploded computational complexity. Based on the probability hypothesis density (PHD) filter in [4], Mahler presented an extension of the PHD filter to multiple extended target tracking, which is called the extended target PHD (ET-PHD) filter [5]. The ET-PHD filter avoids the data association and reduces computational complexity. Two main implementations for the ET-PHD filter are given in [6], [7], called the extended target Gaussian Mixture PHD (ET-GM-PHD) filter and the extended target particle PHD (ET-P-PHD) filter, respectively.
The key challenge of the ET-PHD filter is partitioning the measurement set rapidly and correctly. Therefore, various partitioning algorithms have been proposed. In [6], K-means++ algorithm is applied to partition the measurement set, but it needs the prior information of cluster number, and the partitioning results will become worse in a cluttered scene. In [8], distance partitioning is adopted to partition the measurement set. To get the correct partition, a large number of distance thresholds are needed. With the increase in the measurements generated by a target, distance partitioning makes ET-PHD filter computationally intractable. Based on neural network and adaptive resonance theory (ART), a fast partitioning algorithm with fuzzy ART is proposed in [9]. Although it has fewer partitions and less computational complexity than K-means++ partitioning and distance partitioning, it is sensitive to vigilance thresholds and becomes rather poor in a densely cluttered scene. An effective partitioning method [10] is proposed using spectral clustering technique, where the clutter is eliminated from the measurement set. But this method also needs the prior information of cluster number. The affinity propagation (AP) clustering [11], [12] is introduced into the measurement set partitioning. Although the AP clustering reduces the computational complexity, it is influenced by the initial value of the similarity matrix. In [13], an adaptive partitioning approach is proposed. However, it needs distance partitioning to partition the measurement set first and burdens high computation. In [14], an improved fuzzy c-means (FCM) algorithm is proposed to partition the measurement set. This algorithm sets a probability threshold to estimate the target number and then partitions the measurement set. The partitions number is much fewer, but it correspondingly increases as the target number increases.
When targets are spatially close, the distance between measurements generated by the different targets is small. These measurements may be partitioned into the same cell. It will lead to the underestimation of the target number. As a remedy, the sub-partitioning method is adopted in [6], [9]. The sub-partitioning method solves the underestimation of target number to some extent. However, the above subpartitioning method will result in a large number of additional partitions when there are multiple pairs of close targets whose cells are merged wrongly. More additional partitions will lead to high computational complexity.
In [15], clustering by fast search and find density peaks (CFSFDP) algorithm is proposed for clustering. Compared with the conventional clustering algorithms, the CFSFDP algorithm needs neither the prior information of cluster number nor iteration and it can realize clustering rapidly. Therefore, this paper adopts the CFSFDP algorithm to partition the measurement set. The contributions are as follows: 1) Since the cutoff distance of the CFSFDP algorithm needs to be determined by user experience, the idea of data field is introduced to determine the cutoff distance.
2) The choice of the cluster center in CFSFDP algorithm is relied on the decision graph and it is not a realtime process for measurement set partitioning. A method which can choose the cluster center automatically by the solved cutoff distance and the measurement rate of the target is proposed.
3) An improved sub-partitioning method is proposed. The improved sub-partitioning method decreases the computation burden due to no additional partitions.
The remainder of this paper is organized as follows: Section 2 gives a brief introduction of the ET-PHD filter and the partitioning of the measurement set; the implementation of the proposed algorithm is provided in Sec. 3; numerical simulations are shown in Sec. 4; conclusions are presented in Sec. 5.

Target Motion Model and Measurement Model
In an MTT system, we assume the unknown number of targets is N k , and the target states set is i is referred to as the ith target state. The measurement set obtained at time k is Z k = {z k 1 ,…, z k M k }. z k j is referred to as the jth measurement and M k is the number of measurements. The target dynamic motion model is defined as where f() is the dynamic motion model function and w k i is Gaussian white noise with covariance matrix Q k . Each target state evolves according to the same dynamic motion model independent of the other targets.
The number of measurements generated by each target at each time is a Poisson distributed random variable with rate (x k i ) [16]. () is referred to as the measurement rate of the target and it is a known non-negative function defined over the target state space. The target measurement model is defined as where h() is the measure model function and e k j is Gaussian white noise with covariance matrix R k . Each target is assumed to give rise to measurements independently of the other targets.

ET-PHD Filter
The aim of multiple extended target tracking is to obtain an estimate of the target states X k = {x k 1 ,…,x k N k }. given the measurement set Z k = {z k 1 ,…, z k M k } at time k. It can be achieved under the ET-PHD framework by propagating the predicted PHD and updated PHD. The predicted PHD and updated PHD can be expressed as [5], [6] , where x is the target state. D kk -1 (x) and D kk (x) are predicted PHD and updated PHD, respectively.  kk -1 (xζ) and  k (x) are the spawn target PHD and birth target PHD, respectively. p S,k (ζ) is the probability of target survival and L Zk (x) is the pseudo-likelihood function, which can be expressed as ( ) , ( ) , (5) p D,k (x) is the detection probability. λ k is the average clutter number and c k (z k ) is the spatial distribution of the clutter. p Z k denotes the pth non-empty partition subset of Z k and W  p denotes the Wth cell of the partition p. W denotes the number of measurements in cell W.  Zk (x) = g k (z k x) is the likelihood function for a single target.  p and d W denote the coefficients for each partition p and cell W shown as where  i,j is the Kronecker delta and

The Description of Partitioning the Measurement Set
As shown in (5), all the possible partitions of measurement set Z k are needed in updated PHD. Assuming at time k, the measurement set Z k contains three individual measurements, i.e., Then all the possible partitions can be described as where p i denotes the ith partition and W i j denotes the jth cell of partition i. It can be observed that as the number of measurements increases, the number of possible partitions grows sharply. When all the partitions are used for ET-PHD update, it will make ET-PHD filter computationally intractable. To reduce computational complexity, various measurement set partitioning algorithms have been proposed in [6][7][8][9][10][11][12][13][14]. They consider the most likely subset of all possible partitions to update the ET-PHD filter. The essence of the measurement set partitioning is the clustering problem. Distance partitioning, Kmeans++ partitioning, FCM partitioning and AP partitioning can be regarded as the clustering algorithms. In contrast, fuzzy ART can be regarded as a classification algorithm. Various partitioning algorithms have their own advantage and limitation. Although the above algorithms reduce the computational complexity to some extent, they still have the disadvantage of high computational complexity with the number of clutter and measurements increasing.

The Measurement Set Partitioning
Algorithm Based on CFSFDP

CFSFDP Algorithm
As a clustering algorithm, CFSFDP algorithm calculates the local density  i of each data point and its distance  i from points of higher density. The  i and  i can be expressed as [15], [17] c ( ) where d ij is the distance between data point i and data point j, d c is the cutoff distance. X(x) = 1 if x  0 and X(x) = 0 otherwise. In CFSFDP algorithm, the cluster centers usually have a higher local density, and they are at a relatively large distance from any points with higher local density. The outliers have a lower local density, and they are at a relatively large distance from any points with higher local density. After the cluster centers have been found and outliers have been removed, each remaining point is assigned to the same cluster as its nearest neighbor of higher density.

Determination of d c and the Cluster Center
In this paper, we adopt CFSFDP algorithm to partition the measurement set. It can be seen from Sec. 3.1 that d c has an impact on the calculation of the local density. In the conventional CFSFDP algorithm, d c needs to be determined according to user experience. If d c is too large, each measurement has a higher density. The measurements generated by different targets may be clustered together, leading to the underestimation of target number. In contrast, if d c is too small, each measurement has a lower density. The measurements generated by the same target may be divided into multiple clusters, leading to the overestimation of target number. To avoid the blindness of empirical experience d c , data field method is introduced to determine d c adaptively.
In the data field, suppose that there is a dataset X ={x 1 , x 2 ,…, x n }in a data space . The potential function of an arbitrary data point x i is defined as [18] 1 where m j is the mass of x j , and  is an impact factor. K(x) is a unit potential function, and expresses the law how a data point diffuses its data contribution. x j -x i  is the distance between point i and point j. Since  has a great impact on the distribution of data field, it is adaptively solved by minimizing the Shannon entropy of the potential as follows is the Shannon entropy, and In [19], the above method is introduced to extract the optimal cutoff distance of the CFSFDP algorithm when the local density is calculated by a Gaussian kernel. Similar to [19], we use (12) to determine the cutoff distance. Moreover, it can be seen from (2) that the measurements generated by a target are following Gaussian distribution. For Gaussian distribution, most data stochastically is distributed inside the interval between the expectation plus threefold variances and the expectation minus threefold variances [20]. So d c in our paper is determined by The objective function in (13) is a one-dimensional nonlinear function and it is solved by gradient descent method in this paper.
After d̂c is determined,  i and  i of each measurement can be obtained by (9) and (10). Then, the cluster center is selected manually according to the decision graph of  i and  i [15]. However, manual selection of the cluster center with the decision graph is not a real-time process for measurement set partitioning per time scan. Considering the cluster centers have a higher local density and higher distance, z k i can be treated as a cluster center if  i >  c and  i >  c .  c and  c are thresholds of density and distance, respectively.
We can see that  c reflects the distance between a cluster center and the adjacent measurement. Since d̂c can perfectly reflect the distribution of each measurement in the measurement set, we set  c = d̂c in this paper.  c denotes the number of measurements near the cluster center, that is, the number of measurements generated by a target.
Assume (x) is a constant, i.e., (x) = . Let event A  n denotes there are n measurements generated by a target with , and its probability is expressed as The probability of A  n is shown in Fig. 1 under different . In Fig. 1, it is clear that under a constant , P(A  n ) approaches 0 when n is smaller. This indicates that the extended target is less likely to generate fewer measurements per scan. When  becomes larger, n also becomes larger. Thus, we set a probability threshold p th1 in this paper and the corresponding measurements are regarded as originated-target measurements when the cumulative probability of P(A  i ) is greater than p th1 . Therefore,  c can be defined as th1 0 arg min ( ) According to (15), we can obtain  c adaptively under different .

Partitioning the Measurement Set
After the cluster center is determined, the measurement set Z k can be partitioned by the CFSFDP algorithm. The detailed steps are shown in Tab. 1.

Input:
, p th1 , and  Output: Step 1: Determine the cutoff distance d c by (13) and calculate i  and i  by (9) and (10).
Step 2: Determine  c and c Step 3: Step 4: Remove all the clutter and assign the remaining measurements to the cluster.
Step 5: Get the partition Remark1: When the measurements generated by a target are relatively dispersed, these measurements may be partitioned into multiple cells with the above partitioning algorithm. The wrong cells will lead to the overestimation of the target number after ET-PHD filter update. To solve the overestimation of the target number, a merge operation is added after Step 5. Letting d m̃ñ denote the distance between cluster center z k m̃ and cluster center z k ñ , if d m̃ñ < d̂c is satisfied, z k m̃ and z k ñ are more likely to be generated by the same target. We merge the corresponding cell W̃m̃ and cell

Sub-partitioning
When targets are spatially close, the distance between measurements generated by the different targets will be small. These measurements may be partitioned into the same cell which will lead to the underestimation of the target number. As a remedy, the sub-partitioning method is taken in [6], [9] to solve the problem. Assuming the measurements generated from different targets are independent, the likelihood function for the number of targets N W̃j corresponding to the cell W̃j in partition p can be expressed as [6]   If N̂Wj > 1, the cell W̃j is regarded as a wrong cell and the measurements in W̃j are generated by n W̃j = ceil(N̂Wj) targets. ceil(x) denotes the smallest integer more than x. An additional partition p l = p is added into the list of partitions and the cell W̃j in partition p l is split into n W̃j smaller cells   1 j W n s s W    by K-means method. The original cell W̃j in p l is deleted. It can be seen that the above method solves the underestimation of the target number by identifying the wrong cell by N̂Wj > 1 and adding additional partitions. Therefore, multiple additional partitions are added when there are multiple wrong cells in one partition or there are wrong cells in multiple partitions. The added additional partitions will lead to high computational complexity.
To reduce the computational complexity, an improved sub-partitioning method is proposed. We can see from where n  min denotes the minimum number of measurements when the cumulative probability of the event A  i is larger than p th2 . For the partition p = {W̃j} Ñ j=1 , we calculate min min If N  min > 1, the W̃j is treated as a wrong cell from N  = ceil(N  min ) targets. As the classification result of Kmeans method depends on the choice of initial centers, the CFSFDP algorithm is used to split the wrong cell W̃j. Subpartitioning can be denoted as Since the number of cluster centers is adaptively determined in CFSFDP algorithm, the estimate number N̂ may not always equal N  when the measurements generated by a target are relatively dispersed or clustered. Thus, a modification operation is implemented in this paper. If N̂ equals N  , N̂ centers are reserved as cluster centers of subpartitioning directly. If N̂ is larger than N  , the center which has the largest  i is selected as the first cluster center of sub-partitioning from N̂ centers. Then the distances between the first center and the N̂ -1 remaining centers are calculated, and N  -1 centers with larger distances are selected as the remaining cluster centers of sub-partitioning. If N̂ is smaller than N  , other N  -N̂ centers should be added. The distance factor is considered for the new centers. We calculate the sum of distances between each measurement in W̃j and N̂ centers, and N  -N̂ measurements with larger distance sum are selected as added cluster centers of sub-partitioning. After the above steps, all N  cluster centers are selected.
When all the wrong cells in p are split, the wrong cells are removed and the split smaller cells are added into p. It can be seen that the improved sub-partitioning method makes use of the fact that the number of measurements generated by a target is Poisson distribution, which can be used to identify the wrong cells. The number of wrong cells will decrease, and the improved sub-partitioning method does not need to add the new partition which can reduce computational complexity.
Remark 2: It is worth noting that the d̂c should be recalculated when CFSFDP algorithm is performed to split cell W̃j into N  smaller cells.

Computational Complexity Analysis
To demonstrate the performance of the proposed algorithm, we adopt K-means++ partitioning [6], distance partitioning [8], fuzzy ART partitioning [9], and FCM partitioning [14] as the compared algorithms. All the algorithms are implemented under the framework of ET-GM-PHD filter. In sub-partitioning, the compared algorithms adopt the sub-partitioning method in [6] and the proposed algorithm adopts the improved sub-partitioning method in this paper. For the computational complexity of update in ET-GM-PHD filter, the filter is updated by all cells in all partitions. The computational complexity of the five algorithms is approximated as . n x is the dimension of the target state vector and J kk -1 is the number of predicted Gaussian components. While P of fuzzy ART partitioning and FCM partitioning is much fewer than that of distance partitioning and K-means++ partitioning. For FCM partitioning, the clutter is eliminated and p i  in each partition is much fewer than that of fuzzy ART partitioning. In contrast, the proposed algorithm has the only one partition and the clutter is eliminated. Its computational complexity is much less than that of FCM partitioning.

Target Tracking Setup
Consider the multiple extended target tracking scenarios over the surveillance region [-1000,1000][-1000,1000](m 2 ). The target state vector x k = [x k ,ẋ k ,y k ,ẏ k ] T consists of position [x k ,y k ] T and velocity [ẋ k ,ẏ k ] T . In order to verify the performance of the proposed algorithm, two scenarios are used. 50 Monte Carlo simulations are performed for each scenario. All the simulations are implemented in MATLAB 2015a with a PC (CPU Core i5 4GHz and 5GB RAM). The results are presented in terms of the target estimate number and optimal sub-pattern assignment (OSPA) [22] metric with the parameters p = 2 and c = 100.
The parameters of two scenarios are shown in Tab. 2. Clutter is uniformly distributed over the surveillance region. I 2 denotes 2  2 identity matrix.  v = 2 m/s 2 is the standard deviation of the process noise, and  ε = 20 m is the standard deviation of measurement noise. In gradient descent method, iteration termination threshold is ε = 0.001, step length is ξ = 0.001, and the maximum number of iteration is  max = 100.
The birth PHD is described as where P  = diag([100,900,100,900]), and m  i is initialized according to different tracking scenarios. The spawn target is not considered in the simulation.    Figure 3 gives the corresponding partitioning results with the proposed algorithm at times k = 11 and k = 25 in Scenario 1. At times k = 11 and k = 25, there are two targets and three targets, respectively. It can be seen that the proposed algorithm can divide the measurement set into two and three cells correctly.  Figure 4 shows the estimate of the target number and OSPA in Scenario 1. It is clear that the estimate accuracy of the proposed algorithm approximately equals that of distance partitioning, fuzzy ART partitioning, and FCM partitioning, but it outperforms that of the K-means++ partitioning. This is because K-means++ partitioning often fails to obtain the correct partitions, and it partitions the measurements generated by a target into multiple cells [23], leading to the overestimation of the target number and larger OSPA.  Figure 5 shows the estimate of target number and OSPA in Scenario 2. It can be seen that the K-means++ partitioning still performs worse than the other algorithms when targets are spatially close. Influenced by clutter and sub-partitioning, the target number is still underestimated for distance partitioning and fuzzy ART partitioning. The target number estimate of FCM partitioning has similar accuracy to that of the proposed algorithm. The proposed algorithm splits the wrong cells by CFSFDP algorithm, and it can get more accurate sub-cells, and has smaller OSPA than FCM partitioning. Figure 6 shows the average run time of each step for five algorithms in two scenarios. In Scenario 1, the average run time of five algorithms is 5.07 s, 3.08 s, 0.37 s, 0.12 s, and 0.02 s, respectively. It can be seen that the average run time of five algorithms increases with the increase in target number. However, the average run time of the proposed algorithm is much less than that of the compared algorithms. In Scenario 2, the average run time of five algorithms is 17.19 s, 7.07 s, 0.49 s, 0.20 s, and 0.02 s, respectively. When targets are spatially close, the number of partitions for distance partitioning will increase due to more distance thresholds, which lead to the average run time increasing sharply. Meanwhile, for the compared algorithms, the number of partitions with wrong cells will increase. More additional partitions are added in sub-partitioning, which leads to higher computation. In contrast, the proposed algorithm need not add the new partition in subpartitioning and the average run time is approximately unchanged.

Simulation Result Analysis
To further validate the performance of the proposed algorithm, in Scenario 1, the impact of different  k and different measure noise is analyzed. Figure 7 shows the estimate of the target number and OSPA for distance partitioning, fuzzy ART partitioning, FCM partitioning, and the proposed algorithm under different  k . It can be seen that at a higher clutter rate, the estimate error of fuzzy ART partitioning algorithm increases. The results are based on the fact that the detection threshold in fuzzy ART partitioning might be unreasonable in a densely cluttered scene. Thus, the target number will be overestimated, and OSPA becomes larger. Thus, the estimate error of distance partitioning also increases. Since the clutter is eliminated by elliptical gating in FCM partitioning, the identification of clutter is related to the target predicted location. When the target predicted location is not accurate, the clutter may be considered as the measurement generated by targets and is used for updating. This will increase the estimate error. The proposed algorithm removes the clutter by local density and distance, and it is not affected by the target predicted location. With the increase of  k , the proposed algorithm has approximately the same estimate performance and the estimate accuracy is the highest than the other algorithms.   Figure 8 shows the estimate of target number and OSPA of the proposed algorithm under different measurement noise. From Fig. 8, it can be seen that the estimate of target number keeps nearly the same. This is because the estimate of the target number relies on the accuracy of the partitioned measurement set, and the proposed algorithm could achieve the correct partition. With the measurement noise increasing, the estimate error of target state becomes larger which leads to the increase of OSPA. In general, the proposed algorithm has good performance under different  k and different measure noise, which reflects the strong robustness of the proposed algorithm.

Conclusions
In this paper, a measurement set partitioning algorithm based on CFSFDP is proposed for multiple extended target tracking in PHD filter. Based on the advantage of the CFSFDP algorithm, the measurement set can be partitioned rapidly and correctly. To avoid the blindness of empirical experience cutoff distance and cluster center, the Shannon entropy of potential for measurement set is used to solve the cutoff distance adaptively. Then, the cluster center is determined by the solved cutoff distance and measurement rate. In addition, an improved sub-partitioning method is proposed to deal with the case that the targets are spatially close. As the simulation results show, compared with the other partitioning algorithms, the proposed algorithm has lower computational complexity. Moreover, the proposed algorithm has stronger robustness under some challenging scenarios, such as high clutter rate and different measurement noise.
In the future works, we plan to apply the proposed algorithm to extended target Gaussian inverse wishart PHD (ET-GIW-PHD) filter in [13], [24]. Compared with ET-GM-PHD filter, the ET-GIW-PHD filter can estimate the target extension states. The proposed algorithm can also reduce the computation complexity of the ET-GIW-PHD filter.