Multitarget Tracking Algorithm Based on Adaptive Network Graph Segmentation in the Presence of Measurement Origin Uncertainty

To deal with the problem of multitarget tracking with measurement origin uncertainty, the paper presents a multitarget tracking algorithm based on Adaptive Network Graph Segmentation (ANGS). The multitarget tracking is firstly formulated as an Integer Programming problem for finding the maximum a posterior probability in a cost flow network. Then, a network structure is partitioned using an Adaptive Spectral Clustering algorithm based on the Nyström Method. In order to obtain the global optimal solution, the parallel A* search algorithm is used to process each sub-network. Moreover, the trajectory set is extracted by the Track Mosaic technique and Rauch–Tung–Striebel (RTS) smoother. Finally, the simulation results achieved for different clutter intensity indicate that the proposed algorithm has better tracking accuracy and robustness compared with the A* search algorithm, the successive shortest-path (SSP) algorithm and the shortest path faster (SPFA) algorithm.


Introduction
The purpose of multitarget tracking is to jointly estimate the number of targets and their state of motion from sensor data [1]. During the past decade, it has developed in a variety of directions, such as Air-Traffic Control [2], Marine Monitoring [3], Computer Vision [4], Autonomous Vehicle and Robot [5], etc. At present, multitarget tracking has achieved substantial advances [6]. However, the measurement origin uncertainty, for instance the unknown and time-varying number of targets, clutters, jamming signals and so forth seriously deteriorates the performance of the multitarget tracking system. Resolving the uncertainty of the measurement origin is a computationally expensive task which relied on the prior information about the target motion. To find the mapping from each measurement to its origin is often called measurement-to-track association or just data association [7,8].
Markov Chain Monte Carlo Data Association algorithm (MCMCDA) based on Bayesian Inference [9] and the Probability Hypothesis Density filter (PHD) based on finite set statistics (FISST) [10] have been proposed to cope with this problem of tracking multiple targets with measurement uncertainty. The MCMCDA algorithm can be viewed as a deferred-logic method since a track is decided based on the current and past measurements. It uses the Markov Chain Monte Carlo sampling instead of enumerating over all possible associations. In a PHD filter, it can estimate the target states by recursively computing the first-order moment of the multi-target state posterior probability distribution, without using the complex data association techniques. However, it was not designed to estimate the trajectories of targets. For this problem, Ba-Ngu Vo [11] proposed a newly labeled Random Finite Set (RFS) approach, known as the generalized labeled multi-Bernoulli (GLMB) filter, it can output trajectories and has a better performance in harsh environments.
Bar-Shalom [12] proposed an multidimensional assignment algorithm for solving the data association problem. In essence, the data association problem is converted to a combinatorial optimization problem under certain linear constraints where the total distance/benefit of assigning targets to measurements is minimized/maximized. There is a wide range of algorithms, such as Greedy algorithm [13], Genetic algorithm [14] and Lagrange relaxation theory [15], are used to find the sub-optimum solutions of the multidimensional assignment problem [16]. These approaches, while effective, need to solve the Non-deterministic Polynomial Complete (NPC) with a large amount of data. Goldberg [17] constructed an efficient min-cost flow framework. It has applied a scaling push relabel method to find the optimal solution. Under this framework, Zhang [18] formulated the multitarget data association problem as a maximum a posteriori (MAP) problem. It is mapped into the cost flow network and finds the global optimum solution by depending on the min-cost flow algorithm. An approach combining Dynamic Programming (DP) and Successive Shortest-Path algorithm (SSP) is presented by using Hidden Markov Model (HMM) in [19]. The multitarget tracking problem is formulated as an Integer Linear Programming (ILP) problem, and a greedy, successive shortest-path algorithm is used to reduce the runtime costs. For k = 1, this algorithm can obtain the global optimal solution. For k > 1, it only obtains the approximate solutions, where k is the unknown number of targets. The Shortest Path Faster algorithm (SPFA) is used to solve the Integer Programming problem of the min-cost flow network and quickly obtains the global optimal solution in [20]. The algorithm improves the robustness and tracking accuracy. In [21], an A* based tracking association algorithm is presented. The integer assumption is relaxed to the standard Linear Programming (LP) problem so that the global optimal solution can be obtained by the A* search algorithm.
The above-mentioned approaches are mainly focused on the object tracking based on video image. The available information of image targets are more than that of point targets. In addition, less clutter leads to simple network structure so that the aforementioned algorithms have a good tracking performance. For the problem of multiple point targets tracking in the presence of measurement origin uncertainty, the network may become more complicated that result in an enormous computational burden.
In this paper, a multitarget tracking algorithm based on adaptive network graph segmentation is proposed to address the problem of tracking multitarget with measurement uncertainty. Parallel search strategy is employed to solve the NP-complete problem. The network flow model of multitarget tracking is divided into different sub-graphs. The optimal trajectory is extracted by using the A* search algorithm. Our main contributions are: (1) a parallel network search framework is presented to cope with the multitarget tracking in the presence of measurement origin uncertainty; and (2) we proposed a adaptive spectral clustering algorithm based on the Nyström Method to obtain the network segmentation results for an unknown cluster number data set.
The rest of paper is organized as follows: the problem of multiple targets tracking is formulated as a cost flow network, and transform it into an Integer Programming problem in Section 2. In Section 3, the A* search algorithm is briefly reviewed. The original contribution of the paper is presented in Section 4, where we describe the adaptive spectral clustering algorithm based on the Nyström Method. Simulation results are provided in Section 5. Conclusions and possible extensions appear in Section 6.

Problem Formulation
The multitarget data association problem is regarded as a cost flow network. Let Z = {z i } be a set of measurements, z i = {φ i , ϕ i , t i }, where φ i , ϕ i is the position in x and y-axes, respectively. t i is the time step of the measurement z i . T k = {z k 1 , z k 2 , . . . z k n } represents a trajectory. The set of trajectories is T = {T 1 , T 2 , . . . , T L }, and the number of trajectories L is unknown. The key of the data association is to compute the maximum a posteriori (MAP) estimate of T given the measurement set Z: Assume that tracks are independent from each other. The cost flow network framework of multitarget tracking is as follows: where P(T) is modeled as a Markov Chain. P s (z k 0 ) is the initialization probability of a track starting at z k 0 , P t (z k n ) is the termination probability of a track ending at z k n , and P l (z k j |z k i ) is the transition density from measurement z k i to measurement z k j . P(z i |T ) denotes the likelihood function of measurement z i , which represents a measurement being a true target or a false alarm. In this paper, all measurements are regarded as targets; this means that P(z i |T ) = 1, and the posterior probability of trajectory set T is calculated as follows: To take advantage of the concept of network flow in Network Optimization [22], the indictor variable f i,j is defined as the directed flow variable that from measurement z i to measurement z j . f s,i and f i,t represent the starting flow variable and terminated flow variable, respectively. Depending on the flow conversation method [22], for all nodes z i , the sum of flows arriving at node z i is equal to the sum of outgoing flows from node z i . For any track T k , it is satisfies the following equation: Moreover, the cost flow network must guarantee that only one node is represented by a target in a moment. Let the upper bound of the sum of outgoing flows from node z i is 1. For any node, the constraint is Taking into account that the target may appear or disappear from any location in the cost flow network, the source and sink node are introduced in [18], which is connected to each node, respectively. We transform the Network Optimization problem of Equation (4) into the IP problem, and the logarithm of the objective function can be rewritten as where c s,i is the cost of the flow from the source node to measurement z i , c i,j is the cost of the flow from measurement z i to measurement z j , and c j,t is the cost of the flow from measurement z j to the sink node. Figure 1 shows an example of the cost flow network. The IP problem of multitarget tracking with measurement uncertainty can be described as The cost can be defined as follows: Enter/exit edges Association edges For the IP problem of multitarget tracking, the traditional algorithms have a relatively higher computational cost due to the large number of network nodes. Hence, a parallel processing technology based on the A* search algorithm is presented to find the optimal solution.

Description of the A* Search Algorithm
The A* search algorithm is a heuristic graph search method which guides the search process by using the characteristic information of the problem. For the min-cost flow problem, the A* algorithm searches from the origin node, calculates and estimates the cost of each extended node, chooses the extended node that has a minimum cost and stops it when the algorithm reaches the destination node. Assuming that an evaluation function f * (x) is designed to estimate the minimum cost from origin node s o through node x to destination node s d . The estimated minimum cost is calculated as follows: where g * (x) is the cost from origin node s o to the node x, h * (x) is the lower bound on the minimum cost from node x to destination node s d , and h(x) is the actual value from the node x to destination node . In order to ensure the optimality of the algorithm, Admissibility and Consistency conditions must be satisfied [23]. Admissibility condition: f * (x) never overestimates the true cost of a solution along the current path through. Consistency condition: A heuristic function h(x) is consistent if, for every node x and every successor x of x, the estimated cost of reaching the goal from x is no greater than that of the step cost of getting to x form x plus the estimated cost of reaching the goal from In the implementation of the A* search algorithm, two lists need to be built, named Open list and Close list. Open list is the set of nodes that have been calculated, and that are candidates for the selection of the next node. Close list is the set of nodes that have been selected, and that are not in Open list. In the initial stage, Open list contains a source node and Close list is empty. During the iterative process, the A* search algorithm calculates the evaluation function of each node in Open list chooses the node with the minimum cost and judges whether it is the termination node. If so, the algorithm is done. Otherwise, it extends all adjacent nodes and calculates the cost function of each node. If the solution exists, the A* search algorithm can guarantee obtaining the optimal solution [24].
The A* search algorithm can be expressed as follows: here, the Open list and Close list are denoted as O and C. E is the set of edges: Step 1 Initialization: Step 2 Node Selection: Step 3 Stop Rule: Step 4 Update f * (x j ) and g * (x j ): Step 2.

Multitarget Tracking Algorithm Based on Adaptive Network Graph Segmentation
Under the multitarget tracking environment, the A* search algorithm have two obvious defects that long running time and large storage space. To solve the problem, an adaptive spectral clustering algorithm is presented to segment the cost flow network, and the A* search algorithm is used to find the optimal track of segmented sub-graph. To take advantage of track mosaic technology, the optimal track of each sub-graph is combined, and the combined track is smoothed by the Rauch-Tung-Striebel smoother. The flow chart of the proposed algorithm is depicted in Figure 2.

Adaptive Spectral Clustering
The segmentation problem of graph structure represented by the dissimilarity degree between nodes is the combinatorial optimization problem, which is NP-hard. The general solution is to consider the continuous relaxation form for this problem. Spectral clustering is an unsupervised learning method based on graph theory [25] for arbitrary image shapes. It uses the eigenvalue decomposition of graph matrix to build a spectral mapping space of the original data set, and the new space is partitioned by the K-means algorithm.
Let G 1 (V, E, W) be an undirected graph that transforms from the directed graph G 0 (V 0 , E 0 ) with vertex set V = (v 1 , v 2 , . . . , v N ) and edge set E = (e ij ) i,j=1,2,...,m . We assume that the edge-weighted adjacency matrix of undirected graph W = (w ij ) i,j=1,2,...,m , which is also called affinity matrix, is nonnegative, w ij = w ji ≥ 0. If w ij = 0 that means the node v i and v j are not connected. Here, the weight of two nodes w ij is the cost value c ij in the directed graph model G 0 (V 0 , E 0 ). The sets Ψ 1 ....Ψ k are the subsets of the graph. For sets Ψ 1 ....Ψ k , Ψ 1 Ψ 2 .. Ψ k = V and Ψ i Ψ j = Ø, i = j. Let cut(Ψ 1 , Ψ 2 , . . . , Ψ k ) be the sum of the cuts between sets Ψ 1 , Ψ 2 , . . . , Ψ k : whereΨ i denotes the complement V\Ψ i . The purpose of the spectral clustering is to find the sets Ψ 1 , Ψ 2 , . . . , Ψ k such that MNcut is minimized. The objective function is given as follows: In [26], eigenvectors are clustered in the subspace that is generated by the first k eigenvectors of normalized Laplacian matrix. The degree matrix of the graph D and the normalized Laplacian matrix L sym are defined as Unfortunately, W grows as the square of the number of elements in the grouping problem, and it quickly becomes infeasible to fit W in memory. Hence, an adaptive spectral clustering based on the Nyström Method is proposed to reduce the complexity of the time and space. The Nyström Method is a technique for finding numerical approximations to eigenfunction problems.
Assuming that randomly sample n points from vertex set V and the number of the remaining points is N − n. Now, partition the affinity matrix W as where A ∈ R n×n represents the sub-block of weights among the random samples, B ∈ R (N−n)×n contains the weights between the sample points and the rest of points, and C ∈ R (N−n)×(N−n) denotes the weights matrix between all of the remaining points. LetŪ represent the approximate eigenvectors of W:Ū The approximation of W, which we denoteŴ, can be written aŝ From Equations (17) and (19), C is approximated by B T A −1 B. Therefore, calculating the affinity matrix between remaining points is avoided. It is noteworthy that the columns ofŪ are not orthogonal. We need to orthogonalizeŪ. If A is a positive definite matrix, A −1/2 represents the symmetric positive definite square root of A. Let Q = A + A −1/2 BB T A −1/2 , and diagonalize Q as then the affinity matrix W is diagonalized by U v and Λ s . Without calculating B T A −1 B, a simple approach is proposed to calculated the row sums ofŴ: where a r , b r ∈ R N−n denote the rows sum of A and B. b c ∈ R n denotes column sum of B.
The normalized A and B can be obtained byd. The elements of the normalized A and B are given by The number of clusters is generally determined by human experience and background knowledge. Next, the relationship between the spectrum of the weight matrix and the number of clusters can be obtained by analyzing the affinity matrix of graph.
If v i and v j belong to the same class, then w ij = 1, otherwise w ij = 0. According to perturbation analysis of spectral clustering [27], a permutation matrix P always exists to make elements of any node set V in the sequence belongs to a class after the PV transformation. The affinity matrix W is a block diagonal matrix that consists of n k all-1 matrices. Thus, the elements of Laplacian matrix L sym is divided into k matrix blocks. For each matrix blockL n i , λ k = n i , k = 1, 2, . . . n i − 1, where λ k is the eigenvalues of matrix blockL n i , and λ n i = 0. In light of Matrix Theory [28], the union of the eigenvalues of real symmetric matrix equals that of the block diagonal matrix that consists of these real symmetric matrices. Therefore, the eigenvalues of Laplacian matrix is the union of the eigenvalues of k matrices. The eigenvalues of Laplacian matrix consists of (n − k)(n i − 1) nonzero eigenvalues and k zero eigenvalues. The number of clusters is the number of zero eigenvalues of Laplacian matrixL sym .

k-Short Paths Algorithm
In the cost flow network of multitarget tracking, there may be multiple paths that are directed and paths are edge-and node-disjoint, which means that any two paths cannot share the same edge and node, and a path visits a node in the sub-graph at most once. Here, a path represents a possible track. We reformulated the multitarget tracking problem as an edge-and node-disjoint k-shortest paths problem on a directed acyclic graph. In order to obtain the k-shortest paths, the segmented graph is transformed into a undirected graph firstly, and then the A* search algorithm is adopted to find the single shortest path. If the maximum iteration count is not reached, remove nodes except source node and sink node on the single shortest path, and search for the next path until reaching the maximum iteration count.

Track Mosaic
To segment an undirected graph in the multitarget tracking environment using the adaptive spectral clustering algorithm, it may arise over segmentation. For example, a trajectory may be divided into several segments. In order to obtain the integral track, the mosaic technique is employed to deal with these segmented trajectories. Let T i and T j are two trajectories of different sub-graphs. x io and x jo are the initial position, and x id and x jd are the terminal position, respectively. The Euclidean distance d D is used to decide whether to mosaic two trajectories. If d D < τ, T i and T j are mosaicked, where τ is the mosaic threshold.

Rauch-Tung-Striebel Smoother
A Rauch-Tung-Striebel (RTS) smoother [29] is used to smooth the extracted tracks. The RTS smoother consists of two parts. The first part is the Kalman filter, which calculates the state of target at each time and estimates the corresponding covariance matrix. The second part is backward recursion [30]. In this process, target state and the covariance matrix are taken as inputs to obtain the smoother output.

Time Complexity
The proposed algorithm is parallel processing in that each sub-graph uses the A* search algorithm to obtain the shortest path. The time complexity of the proposed algorithm is mainly composed of two parts, which are the time complexity of the adaptive spectral clustering algorithm and the worst time complexity of searching sub-graph. The time complexity of the adaptive spectral clustering algorithm is O(n × (N − n)) + O(n 3 ) + O(KN I). O(n × (N − n)), O(n 3 ) and O(KN I) are the time complexity of calculating degreed, orthogonal eigenvectors U v and K-means clustering algorithm, respectively. N is the number of nodes in the undirected graph, n is the sample points, K is the number of trajectories and I is the iteration number of K-means algorithm. The worst time complexity of searching sub-graph is O(n 2 max ), and that n max is less than N. The proposed algorithm is adopted to calculate the k-shortest path in O(n × (N − n)) + O(n 3 ) + O(KN I). The A* search algorithm, SPFA and SSP are recognized as three effective methods to solve the k-shortest path problem. The time complexity of these algorithms are O(N 2 ), O(NE) and O(KN 2 ), respectively. E is the number of edge in the undirected graph. The complexity of these algorithms are primarily related to the value of K and N. For multitarget tracking systems, that is, N is large, and the time complexity of the proposed algorithm is far less than that of the algorithms mentioned above.

Experimental Results
In this section, the proposed algorithm was tested in different tracking scenarios. The optimal subpattern assignment (OSPA) [31] metric is used for performance assessment. The experiments have been performed on a computer with an Intel G840 2.8 GHz CPU and 4 GB of memory.

Clustering Quality Evaluation
To evaluate the clustering quality of the adaptive spectral clustering algorithm, the external quality measure (F-score) [32] and (MCR) [33] are used. F-score is used in spam detection for documentation as an overall assessment performance that combines the precision and recall ideas from information retrieval. F-score is defined as follows: where P f (i, j) = N ij /N j represents the precision of the cluster j for the given class i. R f (i, j) = N ij /N i represents the recall of the cluster j for the given class i. N i is the number of the members of class i. N ij is the number of numbers of class i in cluster j. The MCR is given by where C m is the number of misclassified targets. C t is the total number of targets. Figure 3 displays the comparison of classification results with different sampling rates. Since the sampling rate is 1%, clustering result is seriously affected. When the sampling rate exceeds 1%, it is clear that the classification results outperform that shown in Figure 3b. The clustering performance versus sampling rate are shown in Figure 4. It can be noted that F-score increases with increasing sample rate. Once the sample rate exceeds 30%, the data set can be accurately segmented. This is because the similarity matrix between the sampling points can be approximated by that of all points. In Figure 4b, the sampling rate is equal to 1%, MCR is 0.18. The MCR of other sampling rate is equal to 0. Based on the two evaluation criteria above, we choose the sampling rate as 10% in this paper.

Performance Analysis
In this subsection, we consider two scenarios for multitarget tracking. There are two types of dense clutter areas. Inside Type I dense clutter area, clutter points are uniformly distributed in the surveillance region. While Type II dense clutter area is elliptical and the position of its clutter points follows a 2D Gaussian distribution, whose mean is target position at k time and the standard deviations are determined by the major axes of the ellipse. The measurements are obtained from radar which located at [0, 0] m. The measurement model is described as where T is the state variable. The measurement noise v i t is zero-mean Gaussian random vector with covariance matrix R = diag([δ 2 , δ 2 ]), and δ = 10 m in the scenarios 1 and 2.
The target motion can be modeled by combination of constant turn rate (CT) motion and constant velocity (CV) motion [34]. The dynamistic model of the target is described as follows: where F i t is the state transition matrix of the target i at the time t. Under the assumption of CV motion, it is defined as follows: where T denotes the sampling time. In CT motion, it is defined as where ω i t denotes the turn rate. The process noise w i t is zero-mean Gaussian random vector with covariance matrix  [27,30,30,28] s, respectively. The mosaic threshold τ = 10.
In this case, the Type I clutter intensity is κ 1 (z c ) = m Ic V I = 5 6.3×10 7 = 0.79 × 10 −7 , where m Ic is the expected number of clutter measurements in this Type I dense clutter area, and V I is the surveillance region. For the Type II dense clutter, the expected number of clutter measurements in this Type II dense clutter area is 10, and V I I = 4.54 × 10 5 m 2 is the 'volume' of the Type II dense clutter surveillance region. Therefore, the Type II clutter intensity is κ 2 (z c ) = m I Ic V I I = 10 π×600×400 = 1.33 × 10 −5 . We perform a total of 100 Monte-Carlo runs to obtain the average optimal subpattern assignment (OSPA) distance [35] and average of the estimated number of targets. Figure 5 shows the average OSPA distances of the proposed algorithm, the A* search algorithm, SSP and SPFA. It is observed that the average OSPA distance of the proposed algorithm is approximately equal to that of the A* search algorithm, which are better than that of SPFA and SSP. Figure 6 shows an average of the estimated number of targets of the proposed algorithm, the A* search algorithm, SSP and SPFA. It can be seen that the proposed algorithm presents an considerable performance in the estimation of the target numbers. To make a comparison, set the Type II clutter intensity κ 2 (z c ) = 1.33 × 10 −6 , 6.33 × 10 −6 , 19.89 × 10 −6 , 26.53 × 10 −6 , 39.79 × 10 −6 , respectively. The average OSPA distances of the proposed algorithm, the A* search algorithm, SSP and SPFA versus the Type II clutter intensity are shown in Figure 7.  The average OSPA distances of the proposed algorithm have a relatively small difference than that of other algorithms at κ 2 (z c ) = 1.33 × 10 −6 . The average OSPA distance of other algorithms increase rapidly with increasing the clutter intensity. The average OSPA distance of the proposed algorithm is basically no significant change, still at a lower value. When κ 2 (z c ) = 39.79 × 10 −6 , it is obvious that the average OSPA distances of SSP and SPFA is about 40 and 50 times of that of the proposed algorithm, respectively.  Figures 8 and 9. The average OSPA distances of the proposed algorithm, the A* search algorithm, SSP and SPFA are given in Figure 10. The Type I clutter intensity is κ 1 (z c ) = 5 9×10 7 = 0.56 × 10 −7 , and the Type II clutter intensity are κ 2 (z c ) = 1.33 × 10 −6 , 6.33 × 10 −6 , 19.89 × 10 −6 , gaihao26.53 × 10 −6 , 39.79 × 10 −6 , respectively.  As seen from Figure 8, the SSP and SPFA have a larger average OSPA distances. In Figure 9, it is apparent that the average target number estimation of the proposed algorithm is exactly the same as the number of the true targets. In Figure 10, as expected, there is an overall increase of OSPA distances with with increasing the clutter intensity. The average OSPA distance of the proposed algorithm is at a lower value. When κ 2 (z c ) = 39.79 × 10 −6 , the average OSPA distances of SSP and SPFA is about 45 and 70 times of that of the proposed algorithm, respectively.

Run Time
The average running time of different Type II clutter intensity is shown in Figure 11. It clear that the running time of the A* search algorithm, SSP, and SPFA increase exponentially. The running time of the proposed algorithm is growing slowly. When κ 2 (z c ) = 2.222 × 10 −7 , the running time of the A* search algorithm is about 14 times that of the proposed algorithm.

Conclusions
In this paper, we have presented a novel data association framework for multitarget tracking with measurement uncertainty that estimates unknown number and states of targets using the continuous multi-frame data. The multitarget tracking problem was formulated as network flow optimization problem for finding k-shortest paths, and an adaptive spectral clustering algorithm was used to segment the network structure. The optimal solution of each sub-graph can be obtained by the A* search algorithm. Experiment results indicate that the proposed algorithm is helpful in improving the accuracy of track extraction and can reduce the computational complexity. Future work will focus on tracking multitargets with low detection probability.
Author Contributions: T.M. and S.G. conceived and designed the algorithm; T.M. performed the experiments, analyzed the data, and drafted the paper; C.C. and X.S. revised the manuscript.