Adaptive Measurement Partitioning Algorithm for a Gaussian Inverse Wishart PHD Filter that Tracks Closely Spaced Extended Targets

Use of the Gaussian inverse Wishart probability hypothesis density (GIW-PHD) filter has demonstrated promise as an approach to track an unknown number of extended targets. However, when targets of various sizes are spaced closely together and performing maneuvers, estimation errors will occur because measurement partitioning algorithms fail to provide the correct partitions. Specifically, the sub-partitioning algorithm fails to handle cases in which targets are of different sizes, while other partitioning approaches are sensitive to target maneuvers. This paper presents an improved partitioning algorithm for a GIW-PHD filter in order to solve the above problems. The sub-partitioning algorithm is improved by considering target extension information and by employing Mahalanobis distances to distinguish among measurement cells of different sizes. Thus, the improved approach is not sensitive to either differences in target sizes or target maneuvering. Simulation results show that the use of the proposed partitioning approach can improve the tracking performance of a GIW-PHD filter when target are spaced closely together.


Introduction
Multi-target tracking (MTT) typically assumes that each target can produce at most one measurement per scan.Many studies have been done on MTT based on random finite sets (RFSs), such as [1], [2], and this kind of tracking has been used in several fields, such as those in [3], [4].However, along with the development of modern and highresolution sensors, a target may have a larger volume than the sensor resolution cell has and produce more than one measurement per scan.Such a target is called an extended target.Extended target tracking (ETT) is a vibrant area of research and has received increasing attention in recent years [5][6][7][8][9][10].This is especially true for multiple extended target tracking (METT) [11][12][13][14][15][16][17].
An RFS-based observation model for MTT based on finite set statistics (FISST) was first introduced by Mahler.One RFS implementation is the probability hypothesis density (PHD) filter [1], which can estimate target states by computing the function of the first order moment over the target state space.The Gaussian mixture PHD (GM-PHD) filter, presented in [2], is an effective approach to MTT in which the target state is approximated with a Gaussian mixture method.
In contrast to typical single-point targets, target extension is considered by ETT approaches as an extended state.Random hypersurface models (RHMs) for ETT are introduced in [5][6][7].These methods assume that target measurements are produced from a certain hypersurface by random selection.Thus, the targets can be tracked and their shapes can be estimated by computing the parameters of a given hypersurface function.Another typical ETT approach is the random matrix (RM) model [8][9][10].Because target extension is usually treated as elliptical extension, RM models assume that it follows a Wishart distribution.Therefore, the target extension state can be estimated with an inverse Wishart distribution.
To track an unknown number of extended targets, a METT PHD framework is presented in [11] by Mahler.Granström et al. employ both this framework and the Gaussian mixture method, presenting an extended target GM-PHD (ET-GM-PHD) filter for METT [12], [13].This filter uses the distance partitioning (DP) and sub-partitioning (SP) methods to provide partitions, and then the partitions are used to update the Gaussian components by computing the likelihood function.The ET-GM-PHD filter can track an unknown number of extended targets effectively, but it cannot estimate target extension states (i.e., the geometrical target extension must be a predetermined value).To solve this problem, a GIW-PHD filter is presented in [14].This filter assumes that "the last estimated PHD is approximated with an unnormalized mixture of Gaussian inverse Wishart distributions".Therefore, the RM approach can be employed in ETT to estimate the target extension state of a METT PHD filter.More discussion on determining the parameters of a GIW-PHD filter can be found in [15], [16], and the implementation of a GIW-PHD filter in X-band marine radar is introduced in [17].However, when targets are spatially close and performing maneuvers, the target states will be incorrectly estimated.The reason for this is that the partitioning approaches applied fail to provide the correct partitions, which leads to the failure of the GIW-PHD filter in successive scans.This paper presents an extension of the GIW-PHD filter, and its main contributions are as follows: 1) We present an adaptive sub-partitioning (ASP) approach to a GIW-PHD filter in order to solve the partitioning problems that occur when targets are spaced closely together.First, the candidate extension information is computed by using the GIW components that are around the measurements.Then the extension information is used to improve SP by employing Mahalanobis distances.As a result, the proposed ASP will not be sensitive to either the differences in target extensions or target maneuvers.
2) Since the GIW-PHD filter cannot provide the identities of the estimated targets, we present a labeling approach and introduce label evolvement.
The paper is organized as follows.Section 2 introduces the GIW-PHD filter.Section 3 presents the ASP algorithm.The labeling approach for the GIW-PHD filter is presented in Sec. 4. Simulation results are presented in Sec. 5. Section 6 contains our conclusions.

Background: The Gaussian Inverse Wishart PHD Filter
The GIW-PHD filter is an important implementation of the extended target PHD framework.The target states at time k are defined as where N x,k is the unknown number of targets, and x k (i and X k (i) are the kinematic and the extension states, respectively.

Update:
The detection intensity is where the notation pZ k .indicates that the summation is taken of all partitions p in the current set of measurements is the missed detection PHD, and the detection D can be approximated by a mixture of Gaussian inverse Wishart distributions as follows: X are the mean and covariance of the jth Gaussian distribution, respectively.Here,  denotes the Kronecker product of the matrices.
are the degrees of freedom and the scale matrix of the jth GIW distribution, respectively.The updated component weight is given by where d W denotes the probability intensity of cell W, and w p is the weight of a partition p. γ (j) is the mean number of the measurements produced by a target, and β FA,k is a rate parameter to determine the clutter measurements per surveillance volume per scan.δ W,1 is the Kronecker delta and p D is the detection probability.L k (j,W) denotes the likelihood between the jth components and the cell W.
Since we have not made original contributions in this section, only the known basic methods and functions are introduced.The details of the GIW-PHD filter can be found in [14], and the measurement partitioning methods are discussed in Sec. 3.

Adaptive Sub-partitioning Algorithm
Measurement partitioning is an integral part of an extended target PHD filter, because incorrect partitions will lead directly to estimation error.However, the difficulty of partitioning closely spaced targets is still an unresolved issue that leads to serious errors in GIW PHD filtering (see [14], section VI).This section proposes an ASP algorithm to improve the estimation performance of a GIW PHD filter when targets are closely spaced.

Problems and Key Methods
Several partitioning approaches to extended target PHD filtering have been proposed, such as distance partitioning (DP) [12], SP [13], mixture partitioning (MP) [18], expectation maximum partitioning (EMP), and prediction partitioning (PP) [14].However, MP, EMP, and PP are all sensitive to target maneuvers because of their uses of predicted target positions.DP is not sensitive to target maneuvers, but it fails to handle closely spaced targets correctly.SP is an additional partitioning that follows DP to divide the measurement cells produced by multiple targets.Thus, SP is not sensitive to target maneuvers and can handle cases in which targets are closely spaced.However, when targets are closely spaced and their extensions are different, SP will provide incorrect partitions because a clustering error occurs in K-means++ (an algorithm employed in SP) [19].Therefore, the typical partitioning approaches encounter the following problems: 1) no partitioning approach can handle cases in which targets are both closely spaced and performing maneuvers; 2) the use of additional partitioning approaches can indeed improve the performance of the GIW-PHD filter (DP-SP, EMP and PP are used in [14]), but doing so will require more computations.
The SP algorithm is comprised of two steps: detecting and dividing.Suppose    denotes a measurement cell of a partition obtained by DP.First, SP will detect whether W should be divided into sub-cells with where N ̂ denotes the possible number of targets.Pois(•)and • denote the Poisson distribution and the number of elements in a set, respectively.
If N ̂  2, the cell W should be divided into N ̂ sub-cells by the K-means++ algorithm as follows: Step 1: Choose the initial centers using the method in [19].
Step 2: For each j  {1,…, N ̂}, assign each z  to the cluster with the shortest distance D , j between z  and c j .
Thus, N ̂ sub-cells W ͠ j are obtained.
Step 3: Update the centers with Step 4: Repeat Steps 2 and 3 until no W ͠ j values change any longer.
The distance D , j is generally considered a Euclidean distance, meaning Using a Euclidean distance is similar to using a straight line to divide the measurements.However, if target extensions differ sufficiently, such a line may converge to the incorrect position, as shown in Fig. 1.
In contrast to Euclidean distances, Mahalanobis distances are scale-invariant distances that can employ covariance matrices to eliminate the influence differences in target sizes.A key innovation of the proposed approach is employing Mahalanobis distances to improve the K-means++ algorithm.When this is done, (6) becomes where S j denotes the corresponding covariance matrix of W ͠ j .Note that the parameter S j is not given, but it can be calculated using the inverse scale matrix provided by the GIW-PHD filter.The primary method of determining S j can be summarized as follows: Find the predicted components around the cell W. The extension matrices {X̃i} of these components can be employed to calculate the candidates for {S j }.For each subcell W ͠ j , a corresponding S j can be initialized using these candidates.

Do
Step 2 and Step 3 of the K-means++ algorithm.
Calculate the covariance matrix Ψ j for each sub-cell W ͠ j and assign each X ̂i to its most similar Ψ j .For example, suppose X ̂1 and X ̂2 are the covariance matrices of the large and small targets, respectively, in Fig. 1; Ψ 1 and Ψ 2 denote the covariance matrices of cluster 1 and cluster 2, respectively, in Fig. 2 (c).{X ̂1, X ̂2} will be assigned to {Ψ 1 , Ψ 2 } (i.e., S 2 = X ̂1 and S 1 = X ̂2) because of the similarity of their extensions (the details are shown in Sec.3.2).Thus, the partitioning result will converge to the extension similarity maximum.Figure 2 shows an example the partitioning process:  (a) shows the result of the first loop.The triangles are the initial centers, which are obviously inaccurate.
The original K-means++ result may converge to the result shown in Fig. 1, because the use of ( 6) is like using a straight line to divide measurements, and the situation in Fig. 1 demonstrates the best convergence result. However, as shown in (b), Equation (7) ensures that the proposed approach converges to an elliptical extension.Thus, the result in (b) is more similar to two elliptical extensions than that of (a).
 (c) shows the results of the third loop.As we would intuitively expect, the measurements are divided almost correctly with the exceptions of several points.The extension of cluster 2 is similar to that found in X ̂1, and the extension of cluster 1 is similar to that found in X ̂2.Hence, the covariance matrices in (7) become S 1 = X ̂2 and S 2 = X ̂1.
 Finally, the proposed algorithm converges to the result shown in (d).The proposed method can ensure the coverage of results whose clusters are similar to those of X ̂1 and X ̂2.This is why the proposed method can achieve a better performance than the algorithm originally used in the GIW-PHD filter can.
Remark: the proposed method applies only to the situation that the number of elements in {X ̂j} is equal to N ̂.If this condition are not established, the number of targets may be estimated incorrectly in the last time step (i.e., the extension information of the predicted components may be incorrect).Thus, the original K-means++ algorithm should be used in this step.

Initialization
denotes the corresponding d-dimensional set of predicted positions.The states of the predicted components around W can be defined as where  is the maximum threshold of DP.M includes the predicted positions around the cell W, and  includes the corresponding extension matrices.X ̂ikk -1 is the extension matrix and can be calculated by (9) where d denotes the number of dimensions of the physical space [9].V i kk -1 and v i kk -1 are the predicted inverse scale matrix and the inverse Wishart degrees of freedom, respectively.
ASP is only used when M= N ̂.The calculation of the initial centers can give results equal to those of the original K-means++.In addition, to improve the accuracy of the initialization, the centers can be obtained by translating the corresponding , where | 1 k k m  and z are the mean positions of M and W, respectively.c j indicates the relative position and removes the influence of target maneuvers.Thus, the initialization of the covariance matrix is

Modification of K-means++
In the modified K-means ++ algorithm, the dividing process remains the same as in Step 2 and Step 3 in Sec.3.1, but the distance D , j should be calculated by (7) instead.Then the covariance matrix Ψ j of each sub-cell W ͠ j can be calculated by S j is updated based on the similarity of Ψ j and X ̂ikk -1 .

Suppose
denotes a permutation of the elements of  (e.g.,     where   tr  denotes the trace of the matrix.  denotes the

Track Maintenance
The GIW-PHD filter cannot provide the trajectories of individual extended targets.This section proposes a track maintenance approach to a GIW-PHD filter based on a component labeling technique [20].The proposed approach consists of the following steps:

1) Prediction:
Suppose that at time k -1 (k  2), the labels of the GIW components are denoted by where J k -1 is the number of GIW components at time k -1 and that l (j) k-1 denotes the label of the jth GIW component.The labels of predicted GIW components can be written as where l (i) β denotes the labels of the birth GIW components and J β is the number of birth components.

2) Update:
Suppose that Ñ k denotes the total number of all cells in all partitions (i.e., , where Ñ p , k is the number of partitions at time k, and Ñ (i) W,k is the number of the cells of the ith partition).According to (2), there would now be Ñ k + 1 times the number of predicted GIW components.Thus, each predicted component gives rise to Ñ k + 1 corresponding updated components.The corresponding labels can be expressed as where

3) Pruning & merging:
When the GIW components with low weights are pruned, the corresponding labels and their attributes should be also pruned.Note that if several components

Simulation Results
The extended targets are modeled by , where X k is a symmetric positive definite (SPD) matrix, A is the rotation matrix determined by the motion model, and R is the Gaussian measurement noise.The number of measurements of each target follows the Poisson distribution.
The parameters of simulated scenarios are given as 1 s, =4000 4000 m , 6.25 10 where T s is the sensor-scanning interval. denotes the surveillance volume with a rate parameter (i.e., the Poisson mean of clutter measurements is ). R k and Q k are the covariances of process noise and measurement noise, respectively.
The parameters of the GIW-PHD filter are where w 0 is the weight of the birth GIW component.γ (j) is the mean number of measurements produced by each target each time and v 0 , V 0 and P 0 are inverse Wishart degrees of freedom, the inverse scale matrix, and the Gaussian covariance of generated GIW components, respectively.Measurements are partitioned by DP with d = {1,3,5,10,15,…, 30} as the distance thresholds.Then, either SP or ASP is used to divide cells containing more than one target into sub-cells.Thus, the partitioning approaches are called DP-SP and DP-ASP.

Extended Target OSPA (ET-OSPA)
According to [21], the Optimal SubPattern Assignment (OSPA) distance is employed to evaluate the performance of a PHD filter.For METT, the ET-OSPA is defined as where denotes the extension error and {X̃i} and {Ỹ π(i) } are the estimated and true extension matrices, respectively.The first term on the right hand side of (20b), namely the localization error, is equivalent to original function of the OSPA algorithm.The second term is the added extension error.Therefore, the error function d (c)  involves both localization and extension errors.

Results
The target tracks simulated in this study are shown in Fig. 3.The targets moved closer together from 0-20 s and then moved linearly together from 21-40 s.They began circular motion at 41 s and moved uniformly again from 61-80 s.Finally, they separated from each other at 81 s.The edges of the targets separated from each other by 2 m.The shape matrices were X (1) = diag( [20,4]) and X (2) = diag( [10,2]).
Figure 4 shows the average results of 100 Monte Carlo (MC) runs.The ET-OSPA distance of each approach is shown in (a).The ET-OSPA values from 20-80 s are higher when DP-SP is used than they are when other approaches are used.This means that DP-SP cannot provide the correct partitions when targets with different extensions are closely spaced.The results of using DP-SP with EMP and PP are better than those of using just DP-SP because EMP and PP consider the target extension information and can provide the correct partitions.However, the ET-OSPA values increase significantly during the times 20-25 s and 40-60 s.The reason for this is that PP and EMP are sensitive to target maneuvers.The ET-OSPA values of the proposed DP-ASP were clearly lower than those of DP-SP with EMP and PP when the targets were both closely spaced and performing maneuvers.Equation (10) insured that DP-ASP was not sensitive to target maneuvers, while the use of target extension information made DP-ASP insensitive to differences in target extensions.Therefore, DP-ASP could provide more correct partitions than other partitioning approaches could when the targets were closely spaced.This conclusion is also evident in Fig. 4(b), which shows the sums of the weights of each approach.
Figure 5 shows the average time costs of making 100 MC runs for each partitioning approach.The use of more partitioning approaches required a greater number of computations.Thus, the values of DP-SP with EMP and PP were significantly larger than those of the other two approaches.The time costs of DP-ASP and DP-SP were roughly the same because the two approaches provided the same number of partitions.Note that providing a greater number of correct partitions can reduce the number of GIW components.Thus, the time costs of DP-ASP are sometimes lower than those of DP-SP.
Figure 6 shows the association results of a single trial, and Figure 6(a) shows the results of the GIW-PHD filter using DP-SP with EMP and PP.When the targets were maneuvering, many incorrect and excrescent tracks occurred.The reason for this is that the partitioning approaches failed to provide accurate partitions.Then, the inaccurate partitions led to estimation errors in the filter, which ultimately resulted in association problems.The association results when using DP-ASP, shown in (b), were relatively better than those in (a) (i.e., most of the tracks were correct, and the number of excrescent tracks was lower than that of (a)).This means that even though an error may occur during an individual scan, DP-ASP can generally provide more accurate partitions than other approaches when targets are closely spaced.

Conclusions
This paper proposes a measurement partitioning algorithm for a GIW PHD filter, called ASP, which can provide better partitioning results than other approaches can when targets are closely spaced.In addition, a track maintenance approach is included in the GIW PHD filter.In the future works, we plan to apply ASP to other METT particle filters as was done in [22], [23].ASP is a promising approach for dividing particle swarms.Thus, the performance of METT particle filters can be improved when targets are closely spaced.


has the largest weight, the label of the merged component is equal to the label of   | n k k  .
Results of DP-SP with EMP and PP.Results of DP-ASP.