A Labeled GM-PHD Filter for Explicitly Tracking Multiple Targets

In this study, an explicit track continuity algorithm is proposed for multitarget tracking (MTT) based on the Gaussian mixture (GM) implementation of the probability hypothesis density (PHD) filter. Trajectory maintenance and multitarget state extraction in the GM-PHD filter have not been effectively integrated to date. To address this problem, we propose an improved GM-PHD filter. In this approach, the Gaussian components are classified and labeled, and multitarget state extraction is converted into multiple single-state extractions. This provides the identity label of the individual target and can shield against the negative effects of clutter in the prior density region on the estimates, thus realizing the integration of trajectory maintenance with state extraction in the GM-PHD filter. As no additional associated procedures are required, the overall real-time performance of the proposed filter is similar to or slightly lower than that of the basic GM-PHD filter. The results of numerical experiments demonstrate that the proposed approach can achieve explicit track continuity.


Background and Motivation
Multitarget tracking (MTT) jointly estimates the target state and the number of targets simultaneously from a sequence of measurements [1]. In MTT, the time-varying number of targets causes data associations between state and measurement sets of targets, which are uncertain. The key challenges encountered involve detection uncertainty, clutter, and data association uncertainty. Therefore, accurately extracting the states of targets and their trajectories has become a critical issue and a hot topic in the field of MTT. The most popular technologies for MTT are the joint probabilistic data association (JPDA) [2,3], multiple hypothesis tracking (MHT) [4], and random finite set (RFS) [1,5].
Since RFS avoids the traditional data association, it has attracted significant attention. Based on RFS and under the framework of Bayesian filtering, many approximations of Bayes multitarget filters have been proposed, mainly including a probability hypothesis density (PHD) filter [1,6,7], cardinalized PHD (CPHD) filter [8], and multi-Bernoulli (MeM-Ber) filter [9]. These filters are not multitarget trackers, which only estimate target states at individual time instants as opposed to multitarget trajectories. While these filters were not designed with the aim of estimating the trajectories of targets [10], they have been used in many applications [11][12][13].
To provide the information of each track, the notion of labeled RFS is introduced in [14,15]. The generalized labeled multi-Bernoulli (GLMB) filter [14,15] leads to an analytic solution to the Bayes multitarget tracker, which significantly improves the accuracy of multitarget state extraction and explicitly accommodates the estimation of target trajectories. Recently, labeled RFS-based filters and smoothers [10,[16][17][18][19] have been further developed, for generating track estimates, which is also the focus of this paper.
A Gibbs-GLMB filter [10] is a computational efficient implementation of the GLMB filter that combines the prediction and update into a single step. A forward-backward labeled multi-Bernoulli smoother in [17] is proposed for MTT, which can improve both the cardinality estimation and the state estimation, and the major computational complexity is linear with the number of targets. The trajectory Poisson multi-Bernoulli filter and the trajectory Poisson multi-Bernoulli mixture filter [18] have better filtering accuracy and real-time performance than those of the Gibbs-GLMB filter. A onetime step lagged Bayes multitarget smoother using the δ−GLMB density [19], which also inherently outputs targets trajectories, outperforms the Gibbs-GLMB filter on both the estimates and target number and state at the cost of higher computational complexity. However, when targets are in a dense clutter, or misdetected, generating the correct multitarget trajectories is difficult in [16][17][18][19], although the computational complexity problem has been improved in [17,18].
Thus, it is necessary to design a relatively low computing cost filter that can provide explicit trajectories whilst considering the intense clutter environment.

Brief Survey of Related Work
From the viewpoint of practical applicability, the PHD filter is suitable for applications demanding real-time results [20][21][22]. Unlike the full Bayesian filter, the PHD filter iteratively propagates the first moment of the multitarget posterior density. It is a promising technique in terms of computational complexity to solve the MTT problem. However, it does not explicitly accommodate the trajectories of targets [10,21,23]. Furthermore, the PHD filter is known to produce highly uncertain estimates of the target number [21,23]. In recent years, several modifications have been proposed for the trajectories of targets [22,[24][25][26][27][28][29][30].
In [29], an explicit track continuity algorithm for the sequence Monte Carlo PHD (SMC-PHD) filter [31] is proposed, which not only shields against the negative effects of clutter but provides the identities of the individual targets. However, when several newborn targets appear simultaneously in one newborn region, it is unable to extract the state estimates of different observations from these newborn targets. Compared to the SMC-PHD filter, the Gaussian mixture PHD (GM-PHD) filter [32] has the advantages of simple state extraction and low computational cost, which is more suitable for the requirement of real-time scenes. Thus, the focus of this study is the GM-PHD filter.
For the estimation of the target trajectory of the GM-PHD filter, recently proposed solutions involve integrating the track identity into the GM-PHD filter. In each iteration process, each Gaussian component is assigned a label to indicate its correlation with one track. These methods can be mainly classified into two groups. One group [24,25] involves running an additional association method based on the estimates outputted by the filter. The other solution [26][27][28][29] extends a label to the state estimate. However, when targets are in close proximity to each other, in a dense clutter, or misdetected, maintaining the correct multitarget trajectories is difficult, although the accuracy of the estimates of the target number has been increased significantly by improving multitarget state extraction strategies [33][34][35][36][37]. Thus, one can conclude that in the GM-PHD filter, track maintenance and multitarget state extraction have not been effectively integrated. Recently, a Gaussian mixture trajectory PHD filter was proposed in [30], which is able to provide trajectory estimates for live targets, without adding labels or tags. It adapts the estimator for the GMPHD filter for the sets of trajectories; thus, its tracking accuracy is the same as that of the PHD filter. When clutters appear in the high prior density region of the target, it is difficult to shield against the interference of clutters.
An instructive idea is that the Gaussian components approximating the posterior information of one target have the same identity label, which will not be changed throughout the overall filtering process. In the multitarget state extraction, only one state estimate can be extracted from the updated survival Gaussian components with the same label. Then, the explicit trajectories can be obtained by simply connecting the state estimates with the same label, which requires no additional track management. To achieve this concept, in this paper, we propose an explicit track continuity algorithm for the GM-PHD filter.

Main Contributions
The work in this paper is based on a published conference paper [38]. The major contributions of this paper are as follows: • First, each Gaussian component is assigned a label, and the identity labels of Gaussian components are divided into three classes; once it is determined as a confirmed component, its label will be unchanged throughout the whole filtering process. This is convenient for track maintenance; that is, state estimates with the same label at different times belong to the same track. • Second, based on the measurements selected by eliminating most clutter, the updated components are obtained, and their parameters are stored in four different query tables. Based only on these query tables, the multitarget state extraction is converted to multiple single-target state extractions providing the identity label, and the negative effects of clutter in the prior density region on the estimates can be prevented. Thus, we realize the integration of trajectory maintenance with state extraction in the GM-PHD filter. • Third, since the component updating process is based on the selected measurements and no additional associated process is required, the overall real-time performance of the proposed filter is similar to or slightly lower than that of the basic GM-PHD filter. The performance of the proposed approach is verified by simulations designed in a linear scenario.
The remainder of this paper is organized as follows. The technical background is provided in Section 2, while the specific designing of the improved GM-PHD filter is presented in Section 3. The quantitative experiments are described in Section 4, and concluding remarks are provided in Section 5.

GM-PHD Filter
In MTT, the state and observation sets are time-varying and unknown. The collection of the states of targets at time k is represented as an RFS X k = x where |X k | denotes the number of targets, and x (i) k denotes the state of an individual target. Similarly, the collection of available measurements is an RFS Z k = z where |Z k | is the number of measurements and z (i) k is a measurement from a target or due to a clutter. The PHD [6] is the first-order moment of multitarget posterior density.
Assuming that each target follows a linear Gaussian dynamical model and the sensor has a linear Gaussian measurement model [32], where N (·; m, P) denotes the Gaussian density with mean m and covariance P, F k−1 is the state transition matrix, Q k−1 is the process noise covariance, H k is the observation matrix, and R k is the measurement noise covariance. The birth intensity function of the new targets at time k [32] is γ,k , where J γ,k denotes the number of newborn Gaussian components, and w γ,k , and P (i) γ,k are the weight, mean, and covariance, respectively, of the ith newborn Gaussian component. Suppose that the survival probability p S,k = p S,k (x k ), the detection probability p D,k = p D,k (x k ). Based on these assumptions [32], the PHD filter can be efficiently approximated by mixed Gaussian components. w is the weight of Gaussian density and the posterior intensity D k−1|k−1 (x) at time k − 1 [32] is a Gaussian mixture form with J k−1 survival components as follows: The predicted intensity D k|k−1 (x) at time k is also a Gaussian mixture with J k|k−1 components calculated as follows: When a new set Z k arrives, the posterior intensity D k|k at time k [32] is updated and can be described as follows: where For more details on the GM-PHD filter, we refer the reader to the original study [32].
To reduce the computational load, the increasing Gaussian components must be pruned and merged. Then, the posterior intensity D k|k is represented as follows: The multiple-target state estimates can be extracted based on the threshold w th (this is generally 0.5) from D k|k (x). Figure 1 shows the overall procedure of the basic GM-PHD filter.

Problem Formulation
The basic GM-PHD filter in Figure 1 does not explicitly accommodate target trajectory estimation. Trajectory continuity and state estimates essentially influence and interact with each other. In MTT, detection uncertainty, false state estimates, and closely spaced targets have detrimental effects on track continuity. Figure 2 graphically depicts these phenomena. As shown in Figure 2, the clutter rate was relatively high, and the detection probability was 0.9 in the current tracking scenario. In Figure 2a, at times k = 7 and k = 8, the target was misdetected, and its posterior intensity was weakened. This caused its inability to obtain timely state estimates even if the target was redetected at times k = 9 and k = 10. Managing the trajectories of the state estimates at times k = 6 and k = 11 was a challenge. In Figure 2b, at times k = 60, k = 61, k = 62, and k = 65, there were false state estimates caused by clutter points appearing in the high-prior density region of the target. In Figure 2c, two targets were closely spaced, and their state estimates were nearby. It is difficult to identify the state estimate with which they should be associated. The above phenomena make it difficult to provide a correct target trajectory based on the basic GM-PHD filter.
A trajectory is a time-sequence of state estimates with the same label, such as {(x k−1 , ), (x k , ), (x k+2 , )}, wherein provides the means to identify the trajectory or track of one target. Furthermore, even if there is no state estimate at time k + 1, the correct track continuity can still be obtained. To achieve explicit trajectories based on the GM-PHD filter, we proposed an algorithm that labels all Gaussian components, which we called the labeled GM-PHD (LGM-PHD) filter.

The Proposed Multitarget Tracking Algorithm
Assuming the target birth intensity is known a priori, a "one-to-one" principle for the multitarget state extraction of the GM-PHD filter was proposed herein.
The "one-to-one" principle for the GM-PHD filter: In the multitarget state extraction, only one state estimate can be extracted from the updated Gaussian components with the same label, originating from the survival ones.
To abide by the "one-to-one" principle, each Gaussian component is assigned a label, and these labels are divided into three categories in the LGM-PHD filter. In this paper, the identity labels of the Gaussian components are divided into three classes, the labels of newborn Gaussian components, the labels of unconfirmed Gaussian components (approximating the posterior information of targets having their own unconfirmed trajectory), and the labels of confirmed components (approximating the posterior information of targets having their own confirmed trajectory). The attribute of a label is determined based on the value range as follows: where V un is a value far greater than the total number of targets that could possibly appear in this scenario, V new = V un , and V new V un . In general, V un = 200 and V new = 400. r max and r maxt are used to count the confirmed and unconfirmed trajectories, respectively. At time k = 0, r max is initialized as 0 and r maxt is initialized as V un . When a new confirmed trajectory is initiated, r max is incremented by one; when a new unconfirmed trajectory is initiated, r maxt is also incremented by one.
The overall procedure of the LGM-PHD filter contains seven steps-initialization, prediction, observation selection, updating, state extraction, pruning, and merging and capping-as shown in Figure 3. When compared with the basic GM-PHD filter in Figure 1, the submodules marked with " * " are either modified or proposed in this paper.

Initialization and Prediction
At time k, newborn targets appear spontaneously according to a Poisson point process with the intensity function where each newborn Gaussian component is assigned a label According to (15), the initial intensity D 0|0 is described as follows: w With this assumption, at time k − 1, the posterior probability density k−1 is the identity label of the jth survival Gaussian component. Then, according to (4), the priori density at time k is approximated by w predicted Gaussian component is labeled as follows: To obtain the time-sequence of states with the same label, the labels of conformed components remain unchanged throughout the overall filtering process. In addition, the labels of all types of components can be inherited, and only components with the same label can be merged.

Observation Selection
To achieve robust track maintenance, clutter should be eliminated. Considering the overall real-time performance of the filter, it is necessary to eliminate most clutter before the updating procedure, to select effective observations. Given the zero-mean-measurement in which · is the Euclidean distance and a is the confidence coefficient. In this gating algorithm, a = 6. Initializing the effective observation set Z k,e f as ∅, based on all the predicted Gaussian components at is obtained by the operations in Algorithm 1.

Updating
First, four tables U w , U m , U P , and U , all with size J k|k−1 × (1 + Z k,e f ), are created to store parameters w, m, P, and of each updated Gaussian component. They are uniformly called query table U.
Second, based on Z k,e f , the posterior intensity D k|k at time k can be described as k|k . Based on the ith predicted Gaussian component and z (j) k,e f , we obtain the parameter set w k|k−1 . Based on this parameter set, the elements in the ith row and the (j + 1)th column of U w , U m , U P , and U are obtained as follows: . Based on this parameter set, the elements in the ith row and the 1th column of U w , U m , U P , and U are obtained as follows: Third, four tables, U w , U m , U P , and U , are created, which are collectively called When a new confirmed trajectory or a new unconfirmed trajectory is established, its information is documented and stored in U.

State Extraction
In this section, by means of U, certain Gaussian components with the same label are associated with one effective observation, and if the maximum weight in these components satisfies some conditions, one state estimate is extracted.
Complete target status information contains state estimatex, track label , and occurrence time k. Thus, at time k, the state estimate setX k , the track label set L k,st , and the occurrence time set T k,st are all initialized as ∅. The specific steps in state extraction process are as follows.
Step 1: First, obtain the maximum value w in U w (:, 2 : end), where U w (:, 2 : end) represents the elements from the second column to the last column of each row in U w . U w (i, j) denotes the element in the jth column of the ith row in U w . Here, w = U w (i * , j * ) is the maximum, (i * , j * ) = argmax(U w (i, j)), ∀i = 1, . . . , J k|k−1 , and ∀j = 2, . . . , 1 + Z k,e f . Very rarely, w may be in two or more different locations in U w .
Next, if w < w b , jump to Step 4. w b is the threshold of the basic weight. In general, w b is set to 0.02. If w ≥ w b , based on the row and column numbers of w, obtain the unique , and initialize CN = ∅. CN is used to store the column numbers of observations based on which some tracks will be established or maintained. In most cases, the label set only has one element, i.e., |L w | = 1.
Step 2: For i = 1, . . . , |L w |, determine the attribute of (i) w , and execute the following actions for the Gaussian components with

•
First, obtainx byx = H k U m (rn (i) , cn (i) ). The row numbers of elements with Next, the attribute of (i) w will be used to determine which of the following three cases will be executed.
Case A: (i) w ≥ V new , the label of newborn components. As several newborn targets may appear simultaneously in one newborn region, newborn Gaussian components with the same label may contain information on several newborn targets. Therefore, several components with large weights in U w RN (i) w , 2 : (1 + Z k,e f ) can be extracted as the state estimates of different observations, and new tracks can be established accordingly. The specific steps are as follows: 1 The sequence number set CN w l of the different observations with large weights in U w RN (i) w , 2 : (1 + Z k,e f ) is obtained as follows: where w l is the threshold of large weights and uni(A) is an operator that obtains the unique elements in set A. In general, w l is set to 0.4. Additionally, CN = CN, CN w l . 2 If CN w l ≥ 1, for j1 = 1, . . . , CN w l , CN w l state estimates will be extracted for the large weight observations and CN w l new confirmed tracks will be established. The j1th target status information is obtained as follows: In (21), the label of the new confirmed track is given by r max = r max + 1, rn w l is the row number of the maximum weight in U w RN (i) w , cn (j1) w l and the single-target state estimate is H k U m rn w l , cn (j1) w l , which is stored in state estimate setX k and labeled as r max .
The information on the j1th new confirmed track is established by modifying U w , U m , U P , and U as follows: U m (end + 1, :) = U m rn w l , : , U P (end + 1, :) = U P rn w l , : .
where U(end, j) denotes the elements in the jth column of the last row in U, and U(end + 1, :) denotes adding one row at the bottom of U. If j1 > 1, it means at least one state estimate has been extracted at time k. In other words, two or more confirmed trajectories will be simultaneously established from one prior newborn region. Thus, one additional row should be added in U to establish the posterior information of the individual track. To shield against interference from closely spaced components with another label, only the elements in U rn w l , cn (j1) w l can be reserved. 3 Additionally, some components whose weights are between w s and w l in U w RN (i) w , 2 : (1 + Z k,e f ) should also be associated with different observations. The sequence number set CN w bl of these observations with small weights is given as follows: where w s is the threshold of small weights. In general, w s is set to 0.1.
where U(end, j) and U(end + 1, :) are the same as those in (22). Note that no state estimate is extracted here, since state estimates are only extracted from newborn components with large weight or from confirmed components. Storing the parameter set in U makes it easy to determine the observation associated with one updated component with the maximum weight, such as CN w l and CN w bl .
w ∈ L k−2,un , new target status information is obtained as follows: where L k−1,un is the label set of unconfirmed components with small weights at time k − 1, w b < w m < w l . In general, w m is set to 0.2.
The information on this new confirmed track is established by modifying U w and U as follows: Additionally, CN = CN, cn (i) . 2 Otherwise, (i) w is stored in L k,un , i.e., L k,un = L k,un , w < V un , the label of confirmed components. The single-target state estimate is extracted. 1 Preprocessing state estimate. If the observation of one target drifts from its real location, or one target is misdetected at time k − 1, the distribution of this target's posterior information will diverge, and the weights of components associated with its observation at time k will be weakened. Thus, if w < w m , rn w m = arg max( is provided in Section 3.2. Figure A1, the 1st column elements in U represent the updated information for misdetection at time k. When there is some deviation, i.e., w < w m , H k U m (rn w m , 1) is the optimal state estimate since U m (rn w m , 1) is the update of the posterior component with 2 Survival target status information is obtained as follows:
In the effective region, the weights of Gaussian components remain unchanged; however, they are reduced to a 1 w outside the region. Here, a 1 = 0.3. Dividing the effective region of the posterior distribution based on the weights can further shield against interference, such as detection uncertainty and closely spaced measurements originating from clutter or the other target.

•
According to the "one-to-one" principle, the elements in U RN (i) w , : are employed in the state extraction; then, they will be discarded as follows: Step 3: Based on the components associated with effective observations, whose sequence numbers are stored in CN, some tracks are established or maintained. Thus, after the i = 1, . . . , |L w | loop is completed, the elements in U(:, CN) will be discarded as follows: Next, return to Step 1.
Step 4: Here, nr U denotes the number of rows in U, U w , U m , U P , and U , which all have size nr U × (1 + Z k,e f ). Converting these matrices into 1 × (nr U (1 + Z k,e f )) column by column, the parameter set of modified Gaussian components is obtained as follows: First, the maximum weight in U w (:, 2 : end) is w = 0.8256, which is U w (5, 2), as shown in Scheme 1a. As w > w b , obtain U (5, 2) = 1 and initialize CN = ∅. Here, |L w | = 1,x = H k U m (5,2) and RN (i) w = [5,6]. As U (5, 2) < V un , it is the label of the confirmed component. Then, based on (29), survival target status information is obtained. 2]. Via (32) and (33), U w ( [5,6], :) = 0 and U w (:, 2) = 0, Scheme 1b is obtained. It is obvious that the interference of U w (5, 8) caused by z (7) k,e f with state estimation is prevented, which represents the clutter close to z (1) k,e f . Second, the maximum weight in U w (:, 2 : end) is w = 0.5262, which is U w (1, 5), as shown in Scheme 1b. As w > w b , obtain U (1, 5) = 401 and initialize CN = ∅. Here, (1,5) and RN (i) w = 1. As U (1, 5) > V new , it is the label of the newborn component. Based on (20), CN w l = 5, and CN = [CN, 5]. Then, information on one new target status is established via (21), and r max is incremented to 3. The information on this new confirmed track is established by (22).
Sixth, the maximum weight in U w (:, 2 : end) is w = 0.0003 in Scheme 1f. As w < w b , step 1 to step 3 of the state extraction at time k is completed.X k , L k,st , and T k,st are all updated.

Remark 2.
In the proposed filter, since the components approximating one survival target can be identified, the weight threshold for state extraction is reduced compared with the threshold in the basic GM-PHD filter. This ensures that the survival target can be tracked when the proper measurement error occurs, such as U w (7, 3) = 0.2777 in Scheme 1. Additionally, only one state estimate can be extracted from the survival Gaussian components with the same label. This ensures that the interference of clutter on the number of state estimates, such as U w (5, 8) = 0.4703 in Scheme 1, can be prevented.

Remark 3.
The behavior of labeling state estimates in the proposed filter is determined by employing the weights and labels of the updated components, which requires no additional computational processing. Furthermore, the updating process and query tables in the state extraction algorithm are all based on the selected observations. Thus, the real-time performance of the LGM-PHD filter is guaranteed.

Pruning
To reduce the computational cost, reducing the number of Gaussian components by pruning is necessary.
Step 1: Prune confirmed and unconfirmed components that have no state estimates for successive n no time steps. In our filter, n no = 6. 1 The label set L k,no of confirmed and unconfirmed components with no state estimates at time k is: 2 For i = 1, . . . , L k,no , check the occurrence number of k,no appears in L k1,no at each time k1 = k − n no + 1, . . . , k − 1, w (j) This procedure improves accuracy.
Step 2: Based on the truncation threshold w et = 10 −5 , a certain number of the components with the strongest weights are retained in I, Step 3: k|k j = I(1)) , among others.

Merging and Capping
In the basic GM-PHD filter, some of the Gaussian components that are very close together will be approximated by a single Gaussian [32]. This means that, in the merging procedure, the posterior information of targets in close proximity will interfere. This will result in the loss of posterior information, approximated by components with small weights, of some targets. To avoid such loss, in the proposed filter, only the components with the same label can be merged.
First, let d Th be the merging threshold. In general, d Th = 4. Set J k = 0, and I1 = {i1} J k|k i1=1 . The specific steps in merging process are as follows.
Step 1: k|k is obtained.
Step 2: The sequence number set I2 is given as follows: Step 3: Then, a merged Gaussian component is given as follows: Step 4: Execute I1 = I1 − I2. If I1 is not empty, skip to Step 1.
Step 5: The output of merging procedure is obtained as follows: w by those of the J max1 Gaussian components with the largest weights, and the posterior intensity for propagating is obtained as follows: w . In this paper, J max1 = 150.

Steps in the Proposed Algorithm
For completeness, we summarized the key steps of the LGM-PHD filter as follows.
Step 1: Prediction. The predicted components w are obtained via (4) and their labels are obtained via (16). Then, obtain the parameter set of the predicted components.
Step 2: Observation selection. The effective observation set Z k,e f is obtained by eliminating the clutter in the low prior region via Algorithm 1.
Step 3: Updating. Based on Z k,e f , the predicted components are updated via (5) to (12), and the parameters w, m, P, and associated with each component are stored in U and U, respectively.
Step 4: State extraction. As per the steps in 3.4, based on U and U, the multitarget state extraction is converted into multiple single-target state extractions that can provide the identity label, and the occurrence time for each estimate, i.e.,X k , L k,st , and T k,st are obtained. Then, converting the four U values into 1 × (nr U (1 + Z k,e f )) column by column, the parameter set w Step 5: Pruning. As per the steps in Section 3.5, w After Step 4, based onX k , L k,st , and T k,st , explicit track maintenance is achieved by connecting the state estimates with the same label in sequence.

Remark 4.
The proposed filter does not specifically define track disappearance, and each track is left as it is, whether it has disappeared or been misdetected, except that no state estimate appears in seven successive time steps. When a target has been redetected, as long as confirmed Gaussian components for this target are still available for propagation, this target can still have clues to be tracked, and the state estimate will be extracted from these components only if its posterior information satisfies the necessary conditions.
The linear target-originated measurement is given as follows: where υ x,k and υ y,k represent mutually independent zero-mean Gaussian white noise with standard deviations of σ x = 10 m and σ y = 10 m , respectively. The clutter is uniformly distributed over the region −1000 1000 m × −1000 1000 m with an average rate of λ points per scan, i.e., κ = λ/2000 2 , for the basic GM-PHD filter. To ensure a fair comparison, the basic GM-PHD filter employs the gating algorithm that is used in Gibbs-GLMB filter, and G th = 9. The maximum number of Gaussian components is J max = 100 in the proposed filter and the basic GM-PHD filter. The optimal subpattern assignment (OSPA) metric [40] is used to evaluate the tracking performance of different filters, and the order and cut-off are set as p = 1 and c = 100, respectively. The presented result is the average of 250 independent repeated experiments but with independently generated observations for each trial. All experiments were run in MATLAB R2013b on a computer with a 3.3 GHz CPU and 8 GB RAM.
Example 1: Detailed performance comparison under λ = 100 and p D,k = 0.9. As shown in Figures 4-6, the number of clutter points is λ = 100, and the detection probability is p D,k = 0.9. Figure 4 shows the track estimates obtained from the proposed filter for a single trial [38]. The average OSPA errors and computing time for 250 sample runs of each filter are shown in Figures 5 and 6.    Figure 4a shows that the proposed filter can achieve correct trajectory maintenance in MTT in the presence of dense clutter. As the proposed filter does not supplement the state estimate for the misdetected target, there are discontinuities in the tracks. However, this does not affect the trajectory continuity, and the redetected targets can be tracked in time.
Additionally, as shown in Figure 4b, one false track is around H k m (2) γ,k , one false track is around H k m (3) γ,k and two false tracks are around H k m (4) γ,k . This is because the proposed filter has no special procedure to account for the false estimates around the newborn locations. Similarly, they do not affect the continuity of other established tracks, and will disappear over time.
In Figure 5, the average OSPA distances of the proposed filter, the Gibbs-GLMB filter and the basic GM-PHD filter were 22.522 m, 20.497 m, and 31.990 m, respectively. Compared with the basic GM-PHD filter, the average tracking performance of the proposed filter and the Gibbs-GLMB filter increased by 29.60% and 35.92%, respectively. Clearly, the tracking performance of the LGM-PHD filter was slightly lower than that of the Gibbs-GLMB filter. This is mainly caused by the phenomena shown in Figure 4b. However, this did not affect the survival target tracking.
As shown in Figure 6, the average computing times in seconds of a single run (100 time steps) of our MATLAB implementations were 4.60 s (the proposed filter), 30.21 s (the Gibbs-GLMB filter), and 4.96 s (the basic GM-PHD filter), respectively. Compared with the basic GM-PHD filter, the real-time performance of the proposed filter and the Gibbs-GLMB filter increased by 7.26% and −509%, respectively. As expected, the real-time performance of the proposed filter was better than that of the Gibbs-GLMB filter.
Example 2: Performance comparison under various clutter rates for p D,k = 0.95. In this example, the performance of the proposed filter was evaluated with various clutter rates. The clutter rate λ changed from 1 to 100, and the probability of detection was set to p D,k = 0.95. Figure 7 shows the comparison results in terms of the average OSPA errors for 250 sample runs of the different algorithms under different clutter rates, and their computing times are shown in Figure 8.  From Figures 7 and 8, we can see that the proposed filter achieved a better OSPA distance than the basic GM-PHD filter and kept a similar real-time performance as the clutter rate increased. Compared with the basic GM-PHD filter, with λ = 20, the average tracking performances of the proposed filter and the Gibbs-GLMB filter increased by 15.75% and 33.14%, respectively; with λ = 100, they increased by 31.96% and 37.14%, respectively. The advantage of the proposed filter can be fully demonstrated in the dense clutter environment.
Example 3: Performance comparison under various probabilities of detection for λ = 100. In this example, various probabilities of detection were employed to evaluate the effectiveness of the proposed filter. The probability of detection p D,k changed from 0.8 to 1, and the clutter rate was kept unchanged at λ = 100. The comparison results of the average OSPA errors for 250 sample runs of the different algorithms under different probabilities of detection are shown in Figure 9. Figure 10 shows the corresponding computing times of each filter.  Clearly, the proposed filter achieved a better OSPA distance than the basic GM-PHD filter and maintained a similar real-time performance at each probability of detection. As the probability of detection increased, the OSPA distance of the proposed filter decreased. Compared with the basic GM-PHD filter, with p D,k = 0.84, the average tracking performances of the proposed filter and the Gibbs-GLMB filter increased by 23.03% and 29.13%, respectively; with p D,k = 0.96, they increased by 33.21% and 35.79%, respectively.

Conclusions
The paper proposed a computationally efficient LGM-PHD filter that can explicitly accommodate the estimation of target trajectories in MTT. Based on the proposed query matrices, state estimates with identity labels and occurrence times were extracted, and no additional associated procedures were required. Three numerical and simulations analyses demonstrated that the LGM-PHD filter could effectively improve the tracking performance as compared with the basic GM-PHD filter, which was slightly lower than that of the Gibbs-GLMB filter. As the proposed filter could maintain a similar to or slightly lower real-time performance than that of the basic GM-PHD filter, it could be an alternative filter explicitly accommodating target trajectories in real-time processing applications. We should point out that the LGM-PHD filter is not suitable for low detection scenarios. In our future work, we will investigate multisensor technology [12,13,39] for improving tracking performance.
Author Contributions: Y.G. performed the simulations and wrote the paper; D.J. and C.Z. revised this paper and offered some useful suggestions in methodologies and simulations; S.G. offered revision during the whole process. All authors have read and agreed to the published version of the manuscript. .