Multisensor Multi-Target Tracking Based on GM-PHD Using Out-Of-Sequence Measurements

In this paper, we study the issue of out-of-sequence measurement (OOSM) in a multi-target scenario to improve tracking performance. The OOSM is very common in tracking systems, and it would result in performance degradation if we used it inappropriately. Thus, OOSM should be fully utilized as far as possible. To improve the performance of the tracking system and use OOSM sufficiently, firstly, the problem of OOSM is formulated. Then the classical B1 algorithm for OOSM problem of single target tracking is given. Next, the random finite set (RFS)-based Gaussian mixture probability hypothesis density (GM-PHD) is introduced. Consequently, we derived the equation for re-updating of posterior intensity with OOSM. Implementation of GM-PHD using OOSM is also given. Finally, several simulations are given, and results show that tracking performance of GM-PHD using OOSM is better than GM-PHD using in-sequence measurement (ISM), which can strongly demonstrate the effectiveness of our proposed algorithm.


Introduction
Multi-target tracking has been studied widely in recent years. It has attracted a lot of attention and is a major research hotspot in the field of target tracking. Multi-target tracking is not a simple extension of single-target tracking; it involves many situations that will not occur in single-target tracking. For example, in a multi-target tracking scenario, the number of targets is time-varying because new targets may appear and original targets may die. Therefore, a multi-target tracking algorithm needs to estimate the number of targets and the state of each target, at the same time. Because of the existence of clutter, it cannot guarantee that every measurement comes from a real target.
Multi-sensor measurements at each sampling time correspond to many different targets (or clutters). The problem of correlating measurements with targets is called data association, which is the key point and difficulty of a multi-target tracking algorithm and has been the focus of multi-target tracking research. Today, joint probability data association (JPDA) and multiple hypothesis tracking (MHT) are regarded as representative algorithms to solve data association problems [1][2][3][4]. There are also many improved algorithms based on JPDA and MHT. For example, [5] studied the JPDA filter with unknown detection probability and clutter rate, and [6] studied the problem of passive sonar tracking using MHT algorithm. It is noted that key developments of MHT over the past 40 years have been already reviewed in [7].
Probability hypothesis density (PHD) filter is a multi-target tracking algorithm based on random finite set (RFS) that has emerged in recent years. It was first proposed by Mahler in the literature [8]. This algorithm avoids the problem of data association. Thus, it can effectively solve the problem of a large amount of computation in a classical data association algorithm. A comparison of the RFS approach and traditional multi-target tracking methods has been given in [9]. Under the framework of RFS, the set of states of targets is regarded as a set-valued state, and the set of observations is regarded as a set-valued observation. Thus, RFS is the natural representation of the states of targets and observations of the targets. After the state and observation are modeled as RFS, a multi-target tracking problem can be formulated in Bayes filtering. Because multi-target Bayesian filtering involves multiple integrals, it will result in computational intractability. The PHD filter can alleviate the computational intractability because it operates on the single-target state space. However, there is no closed-form solution for PHD filter, so it is a common way to use the sequential Mento Carlo method to approximate it [10]. Another common implementation of PHD is Gaussian mixture probability hypothesis density (GM-PHD) [11], in which multiple Gaussian components are used to approximate PHD. Under the assumption of the linear Gaussian multi-target model, the analytical solution of GM-PHD can be obtained. This filtering algorithm is also adopted in this paper. In addition, these three multi-target tracking approaches mentioned above have been widely used in sonar [12,13], autonomous vehicles [14,15], fusion [16][17][18], visual target tracking [19,20], and sensor networks [21,22]. There are also some multi-target tracking algorithms and related work in the framework of a particle filter; for example, [23] introduced an efficient particle filter for multi-target tracking which solved the problem of the curse of dimensionality, [24] proposed two sampling techniques based on group importance sampling, and the problem of model selection in the framework of particle filter is studied in [25,26].
In a multi-target tracking system, target information measured by sensors is transmitted to the fusion center by the wireless network. However, due to network delay and other reasons, it is inevitable that some of measurements may arrive at the fusion center later. In this situation, we call it out-of-sequence measurement (OOSM). OOSM is a widely studied problem, and the main difficulty is how to use past measurements to update the current state. The research of one-step-lag OOSM in single target tracking is mature. The representative works based on Kalman filter framework are A1, B1, and C1 algorithms [27][28][29], where A1 is the optimal algorithm, and B1 and C1 are sub-optimal algorithms. In recent years, a particle filter has been widely used in the tracking system. Thus, OOSM based on a particle filter framework has also been extensively studied. Y-bar Shalom gives the optimal algorithm based on a particle filter [30]; however, this method requires excessive computational resources. To overcome this drawback, an efficient particle filter using OOSM is proposed in [31].
Although the OOSM problem and the multi-target tracking problem may appear in the tracking system simultaneously, there is little related work on the combination of them. Moreover, to the best of our knowledge, there is no related work about multi-target tracking based on RFS using OOSM that has been before. Obviously, in the actual tracking system, using OOSM in multi-target tracking is very common. In view of this, a novel multi-target tracking algorithm based on RFS using OOSM is proposed. When OOSM arrives, it is used to re-update the current state, instead of being discarded directly.
The structure of this paper is as follows. Section 2 gives the formulations of the B1 algorithm which is used to solve the problem of single-target tracking when OOSM arrives, and the PHD filter is also introduced in this section. Section 3 presents the main contribution of this paper, namely the formulation of re-update of GM-PHD filter using OOSM. And the implementation of GM-PHD filter using OOSM is also given in this section. To demonstrate the effectiveness of our proposed algorithm, several simulation results are given in Section 4. Finally, the analysis and conclusion are given in Section 5.

Problem Formulation
This section is arranged as follows. At first, the model of the system is given. Next, the problem of OOSM is formulated. Finally, the RFS based GM-PHD filter is summarized.

System Model
Consider a dynamic system as follows: where F k is state transition matrix to time t k+1 from time t k , x k is the state vector, and w k is the process noise for this interval. The F k has different forms when the target adopts different motions, and specific F k related to our experiment is given in the Section 4. The state of target at time t k is given by represent the velocities in the x,y coordinate, and (x, y) is the target location. And the process noise w k is always assumed as zero-mean, white Gaussian noise: where Q k = E{w k w T k }. The measurement equation of system is given as follows: where H k is the measurement matrix of sensor; it also has different forms when sensors are different, and v k is the measurement noise which is always assumed as zero-mean, white Gaussian noise: where R k = E{v k v T k }.

OOSM B1 One-Step-Lag Algorithm
The OOSM's problem can be formulated as follows. At time t k , we use the measurement Z k to update the target state to get x k ; subsequently, an earlier measurement's Z d with time stamp t d arrives, which satisfies the following in Equality (5), where Z d is called an one-step-lag OOSM. The goal is to re-estimate the current state x k with OOSM Z d .
Applying Equation (1), we can get the backward transition equation as follows: where x d is the state of target at time d, F(k, d) is the state transition matrix to time t k from time t d , and w k,d is the process noise during this interval. And then, we can get the backward transition equation from (5): where F d,k is the backward transition matrix, and it is also the inverse matrix of F k,d , which means For a one-lag measurement problem, in which t k−1 < t d < t k , the schematic diagram is shown in Figure 1.

Measurement arrives at fusion
center: Time stamp of measurement: Respectively represent measurements of sensor 1,2 Out of sequence measurement Based on the LMMSE criterion, using the basic equation of linear estimation, we can get the backward prediction of the state from t k to t d as follows: whereŵ k,d|k is the estimation of process noise of backward prediction, and it is assumed to be zero in the sub-optimal B1 algorithm, which means under this situationŵ k,d|k = 0, and based on this point, we can get the following equation:x Then the covariance matrix of backward state prediction can be calculated as follows: The related covariance matrix in the Equation (10) can be calculated as follows: Here, we just give the equations of the B1 algorithm, and the derivation of the equations in detail can be found in [27]. The covariance matrix of the measurement of backward state prediction is The covariance matrix between state x k and measurement z d is Then we can calculate the filter gain as follows: therefore, the re-update of current statex k|k using the OOSM z d iŝ

GM-PHD Filter
As mentioned in the previous section, the multi-target tracking algorithm based on random sets effectively avoids the problem of data association. Moreover, it does not need to know the number of targets in advance. In a word, it brings a new idea to solve the problem of multi-target tracking.
In multi-target tracking scenarios, the birth and death of the target should also be considered except from the target motion. Thus, it is necessary not only to estimate the state of multiple targets, but also to estimate the number of targets. This means that the number of elements in the multi-target state set and measurement set is variable. So, in a multi-target situation, the states of targets and measurements can be represented by the following: where n(k) and m(k) are the number of targets and the number of measurements at time t k , respectively, and F (E s ) and F (E o ) are the multi-target state space and measurement space.
In the multi-target scenario, the birth and death of the target should also be considered separately from the target motion. At time t k−1 , each target either survives at time t k with the probability p s,k , or dies with probability 1 − p s,k . At time t k , a new target may arise either by newborn or by spawning from a target at time t k−1 . Therefore, the multi-target state of next time t k can be modeled by the following RFS: where S k (X k−1 ) is the RFS of multi-target state at time t k , which evolves from time t k−1 , B k (X k−1 ) is the RFS of the newborn target, and Γ k (X k−1 ) is RFS of targets spawned from the previous state. Due to the existence of clutters in measurement, the RFS of measurement can be modeled as: where Θ k (X k−1 ) is the RFS of measurement of targets, and K k (X k−1 ) is the RFS of measurement of clutters. The PHD filter is a typical RFS-based multi-target tracking algorithm. However, there is no closed form solution to PHD recursion, in general. But the PHD recursion has a closed-form solution in a certain class of situation, in which the multi-target model satisfies the linear Gaussian assumption. In this situation, the Gaussian mixture is used to approximate the intensity.
For clarity, the Gaussian mixture PHD (GM-PHD) filter is summarized in this paper, and readers can find more details in [11]. Under linear Gaussian assumption, the intensity at time t k−1 can be represented by the following equation: Then, the predicted intensity is given as follows: where v S,k|k−1 (x), v β,k|k−1 (x), v γ,k|k−1 (x) represent intensity of survival target, intensity of new target spawned by the previous target, and intensity of newborn target, respectively. These three intensities are also Gaussian mixture, and their specific expressions can be found in [11]. Next, we can get the posterior intensity, which is given as follows: where p D,k is the probability of detect, v D,k (x; z) is the update item of predicted intensity, and the expression of ω k|k can be found in [11]. Finally, the pruning and merging is performed, and readers can also find the details in [11]

Multi-Target OOSM Tracking Algorithm
In the above two sections, the problems of OOSM and multi-target tracking based on RFS were discussed. In fact, these two issues have been widely studied. And the issue of multi-target tracking using OOSM has also attracted a lot of attention recently. However, to the best of our knowledge, there is no relative work about RFS-based multi-target tracking using OOSM. In the following, a novel multi-target tracking algorithm based on RFS is proposed, in which OOSM is used to re-update the posterior intensity. As we can see from Equations (9) and (12), predicted covariance P k|k−1 is needed to calculate the retrodiction covariance P d|k . However, in the GM-PHD algorithm, after pruning and merging, there is only posterior information left in the Gaussian component. In this situation, OOSM cannot be used to re-update the posterior intensity. Therefore, when measurement arrives, the judgment of whether the measurement is a delay should be executed first. If the timestamp of measurement shows that the measurement is in sequence, the update process of the GM-PHD filter is performed. On the contrary, if the measurement is a delay, our proposed algorithm is performed to re-update the multi-target state at time t k . The flowchart of the proposed algorithm is given in Figure 2. Recalling the update equation of GM-PHD filter, substituting Equations (23) and (25) into update Equation (24), then the update equation can be rewritten as follows: where the |Z k | is the number of elements of set Z k . In Equation (26), the first item in the right side is the prior intensity, which is associated with the estimated statex k|k−1 . It also represents targets which are not detected, and we refer to it herein a 'prior target'. The second item represents detected targets, which we refer to herein as 'posterior target'.
The 'prior target' is the target which is not detected at time t k , but when OOSM arrives, it may contain information about the 'prior target'. Therefore, when OOSM arrives, the estimated statex k|k−1 associated with 'prior target' should be updated. Analogously, the 'posterior target' is also needed to re-update using OOSM, and the weights of Gaussian components also need to be modified. The diagram of the classification processing algorithm using OOSM is shown in Figure 3.  In the following, the derivation of the proposed algorithm is given.

Backward State Prediction
For simplicity, we assumed that there is no targets newborning and spawning from exist targets during time t d to t k , which means the number of targets at time t k and time t d are the same. The multi-target state at time t d can be represented by the following: where S d|k (x) is the RFS of multi-target state at time t d which evolve from time t k , X k|k−1 is the RFS of undetected targets' state, and X D,k is the RFS of detected targets' state. The intensity at time t d can be divided into two items. One is retrodicted from the undetected intensity (1 − p D,k )v k|k−1 (x). The other is retrodicted from the detected intensity ∑ z∈Z k v D,k (x; z). Applying retrodiction Equations (9) and (10) of the B1 algorithm, the intensity retrodicted from the undetected intensity is given as follows: where where P ww(i) k,d|k−1 is the ith Gaussian component's covariance of process noise from time t d to time t k , and P xw(i) k,d|k is the ith Gaussian component's cross-covariance between state and process noise. They can be represented as follows: Substituting Equations (31) and (32) into (30), we can rewrite Equation (30) as follows: Analogously, we can easily get the intensity retrodicted from the detected intensity: where ω (i) where P (i) k|k , P ww,(i) k,d|k , and P xw,(i) k,d|k can be derived from Equations (10)- (12). And the pseudo-code for backward state prediction is given in Algorithm 1.

Algorithm 1 Backward state prediction
, and the out of sequence measurement Z d Step 1. Backward prediction

Re-Update with OOSM
It has been mentioned in Equation (26) that intensity v k (x) can be divided into two items. Therefore, when OOSM arrives, it is used for re-updating these two items, respectively, by using OOSM, and we can get the following re-update equation: where the first item represents the undetected intensity at time t d , and the second item represents the detected intensity at time t d . And intensities v D,k|k−1,d (x; z), v D,k|k,d (x; z) are given by: where K d|k can be derived from the Equations (13) and (15). The pseudo-code is given in Algorithm 2. And weights ω k|k,d are given in the next section.

Weights Correction
In the above section, we get the intensities that are re-updated with OOSM. These two intensities are re-updated using the same measurement Z d , which is different from Equation (39). Weights of the intensities need to be corrected, and corrected weights are given as follows: Prior target: Posterior target: where For completeness, the pseudo-code for GM-PHD using OOSM is given in Algorithm 3.

Simulations
In this section, the distance which is used to measure the performance of multi-target tracking is introduced first. Then another two algorithms are given to compare with the proposed algorithm.
Finally, several simulation scenarios are given to prove the effectiveness of our proposed algorithm. And simulation results show that the proposed GM-PHD using OOSM has almost the same performance as the GM-PHD using in-sequence measurement (ISM).

Distance
It is also an important issue to measure the tracking performance properly, which directly leads to the choice of the tracking algorithm. Improper criterion will result in tracking accuracy that cannot meet the requirement. Root mean square error (RMSE) is used to measure the distance error between estimate position and real position in a single-target tracking situation. Analogously, in a multi-target tracking scenario, RMSE is also an indicator of tracking performance. However, there is more criterion than that. For example, the accurate rate of the track associate is also an important indicator, but it does not work because the GM-PHD filter based on RFS avoids using the data association technique. Thus, Mahler [32] and Vo [33] proposed three general criteria to measure the effectiveness of multi-target tracking. The essence of these criteria for multi-target tracking performance is to measure the distance between the real state set and estimate state set.
Optimal subpattern assignment (OSPA) distance is an improvement of Wasserstein distance. The definition of OSPA distance is given as follows: where d (c) = min(c, ||x − y||), || · || represents the Euler norm, c is the truncation distance, |X| is the cardinality of set X, and Π n represents all permutations of {0, 1, · · · , n}. The OSPA distance can be divided into location error and cardinality error which, respectively, correspond to first item and second item ofd p (X, Y). Thus, it can evaluate the performance of multi-target filtering more comprehensively. It also serves the shortcomings of Wasserstein distance. And it has definition even when one of the two sets is empty or all the two sets are empty. Parameter c and p affect the weights of cardinality error and location error. How to choose the values of c and p can be refer to [33]. In general, the OSPA distance is currently recognized as an indicator of multi-target filtering.

Compared Algorithm
To prove the effectiveness of the proposed algorithm, two other algorithms are given here for comparison.
The first is GM-PHD using ISM. Under this circumstance, the order of arrival of two sensors at the fusion center is the same as that of sampling. Then, the sequential fusion of asynchronous measurement based on the GM-PHD filter is performed in fusion center. In another situation, the OOSM is discarded directly. Estimation of targets' states is performed using only in sequential measurements. This situation is called Frame Drop here. The representations of these two situations are shown in Figures 4 and 5.
Respectively represent measurements of sensor 1,2 Measurement arrives at fusion center: Time stamp of measurement: In sequence measurement

Simulations Results
The simulation scenario designed here is in the 2-D plane. In the detection area of [−800, 800] m ×[−800, 800] m , there are four targets appearing, all of which adopt the constant velocity (CV) model. We recall the dynamic equation as follows: x k+1 = F k x k + w k , when the motion of target is CV, the transition matrix and the covariance of process noise are given as follows, respectively: where q v = 1(m/s 2 ), T is sampling time, and we set it to 5s. Each target has survival probability p S,k = 0.99, and each target is detected with probability p D,k = 0.98. The observation model has been given in Equation (3). The measurement matrix and covariance of measurement noise are given as follows: where σ x = σ y = 5 m . Measurements are immersed in clutter, and the clutter can be modeled as a Possion RFS K k with intensity: where u(z) is the uniform density over the surveillance area, V is the volume of the surveillance region, and λ c is the average number of clutter returns per unit volume.
The algorithm uses the Gaussian component pruning method in the GM-PHD filter. Pruning threshold is T w = 10 −5 , merging threshold is U = 4 m, and the maximum number of Gaussian component is J max = 200. Gaussian component with weights over 0.5 is considered as the estimation of state of target.
Simulation results are averaged by 100 Monte Carlo experiments. The real trajectory of each simulation is the same while measurements are regenerated each time. The real trajectory of targets and measurements of 100 simulations are shown in Figure 6. The tracking results of three algorithms are shown in Figure 7. It is obvious that the tracking performance of GM-PHD using OOSM and GM-PHD using ISM is nearly the same, while the performance of Frame Drop GM-PHD has decreased significantly.   A detailed analysis of the performance of tracking algorithms is given below. The real number of targets at different times and the number of targets estimated by different algorithms in the presence of clutter are given in Figure 8a. And the estimated number and real number in the absence of clutter is given in Figure 9a.     From Figure 8a we can see that in the presence of clutter, the result of target number estimation of GM-PHD using OOSM is even more accurate than GM-PHD using ISM. Meanwhile, the estimation of the number of GM-PHD discarding OOSM has a great deviation. And as time goes by, the algorithm even begin to diverge. The more detailed data can be seen in Figure 8b. From Figure 9, we can see that these three algorithms have almost the same performance in the estimation of number.
The simulation results of OSPA distance used to measure the results of multi-target tracking are also given as follows.   From Figure 10b, we can see that under OSPA distance criteria, our proposed algorithm which can process OOSM has almost the same performance as that of GM-PHD using ISM in the absence of clutter. However, as we all know, the clutter exists in most tracking scenario. From Figure 10a, we can see that our proposed algorithm using OOSM has better performance than GM-PHD using ISM in the presence of clutter. At the same time, in the presence of clutter situation, the performance of GM-PHD has decreased evidently when OOSMs are discarded directly. It is easy to understand that discarding OOSM directly leads to a reduction of information and consequently results in performance degradation. Compared with GM-PHD using ISM, GM-PHD using OOSM can be seen as doubling the sampling interval while the sampling information of each sample double. Due to targets are all adopting CV model, the maneuverability is not strong. Therefore, compared with the increase of information obtained by each sampling, the increase of sampling interval has little effect on the estimation accuracy. Thus, our proposed algorithm which re-updates intensity using OOSM outperforms GM-PHD using ISM in performance. Where clutter = 0 represents there is no clutter in the tracking region, and clutter = 4 represents that the average number of clutter in tracking region is 4 at each sampling time. The detailed data is shown in Table 1 which can adequately demonstrate the effectiveness of our proposed algorithm.

Conclusions
The issue of OOSM appears frequently in tracking systems. To better use the OOSM in multi-target tracking scenarios, a novel algorithm was designed in this paper. Firstly, the GM-PHD filter is introduced, and then the B1 algorithm which can solve the OOSM problem in a single-target tracking scenario is also given. Based on this, we used OOSM to re-update the posterior intensity of GM-PHD. Under the framework of the Kalman filter, we derived closed-form recursion for the intensity which is re-updated by the OOSM. Implementation of GM-PHD using OOSM was also given. Finally, several simulations were given to show that the tracking performance of our proposed algorithm is almost the same as that of GM-PHD using ISM in the absence of clutter, and the tracking performance of our proposed algorithm is even better than GM-PHD using ISM in the presence of clutter. This strongly proves the effectiveness of our algorithm.