Improved GM-PHD Filter with Birth Intensity and Spawned Intensity Estimation Based on Trajectory Situation Feedback Control

: The Gaussian Mixture Probability Hypothesis Density (GM-PHD) ﬁlter can effectively track multiple targets in a single scenario. However, for GM-PHD, unknown target behavior, e.g., target birth or target intersection, produces difﬁculties in terms of accurate estimation. First of all, GM-PHD assumes the model parameters about the birth target are prior information, which results in the inability to detect the birth target that occurs at random in complex scenarios. Then, since the measurements generated by the intersected targets overlap each other, GM-PHD cannot distinguish these targets, resulting in a biased estimation of the state and number of targets. To solve these problems, this paper proposes an improved GM-PHD ﬁlter with a birth intensity and spawned intensity updating method based on the trajectory situation feedback. In the ﬁltering process, the trajectory initiation feedback formed by the rule-based correlation of Gaussian components is introduced to GM-PHD to adjust the birth intensity in real time, which is used to improve the detection of birth targets. Simultaneously, the analysis of trajectory situation is designed to determine the relative motion trend between targets. On this basis, the ﬁlter improves the recognition of the intersected targets by enhancing the spawned intensity. Simulation results demonstrate that the proposed algorithm achieves better performance on the state and number of targets in complex scenarios, and shows superiority to other GM-PHD ﬁlters. Validation, C.Z., Y.Z. and Z.L. Investigation, Z.L. Resources, T.Q.; Data Curation, T.Q.; Writing—original Draft Preparation, C.Z.; Writing—review & Editing, C.Z. and Z.L. Visualization, Y.Z.; Supervision, Z.L. Project Administration, T.Q.


Introduction
With the rapid development of sensors, multi-target tracking (MTT) technology based on sensors (e.g., radar, infrared, and sonar) plays an important role in both military and civil fields [1][2][3][4]. MTT used to mainly rely on radar detection and orientation estimation. With the development of sensors, optical imaging technology is gradually applied in target tracking, which filters the optical target according to the motion model. In general, whether radar or optical imaging, MTT involves the analysis of the measurements from sensors. Generally, most traditional MTT formulations involve explicit associations between measurements and targets. Multiple Hypotheses Tracking (MHT) [5] propagates the association hypotheses between measurements and targets in time. The Joint Probabilistic Data Association (JPDA) [6] uses observations weighted by their association probabilities. Both of them can achieve perfect tracking performance in a low cluster scenario. However, these association-based approaches generally suffer from the combinatorial computations arising from the association of dense measurements with targets.
In contrast, the MTT approaches based on RFS formulate the tracking problem in the Bayesian framework without dealing with the complex association between measurements and targets, and have been popular for use in airplane tracking, vehicle tracking, and so on. The Probability Hypothesis Density (PHD) realizes MTT by transmitting the posterior intensity to iterate the target state [7][8][9], and it is the most popular and typical filter among those MTT approaches based on RFS. The Gaussian Mixture PHD (GM-PHD) filter [10] has the characteristics of high filtering efficiency and simple state extraction, which is widely used in the MTT system in the linear Gaussian dynamic model.
To guarantee the stability of cardinality estimation during the propagation, the Cardinality PHD (CPHD) [11] simultaneously propagates the posterior intensity and cardinality distribution. However, PHD/CPHD assume the birth intensity and spawned intensity as prior information, which leads to the decline in the tracking performance when prior information mismatches the real scenario.
The birth intensity containing the initial states of the birth targets determines the performance of the PHD filter. To adaptively approximate the birth intensity, the Track Initialization Based on GM-PHD (TIB-GM-PHD) [12] proposes the birth intensity estimation method based on trajectory initiation with the constraints in position and velocity. However, it only considers the association between measurements at adjacent moments, which leads to an increase in the false alarm rate. In [13], the intensity is propagated in the form of formal set and temporary, respectively. The state of birth targets is estimated through the Multi-Bernoulli RFS in [14]. All the above methods can be applied in GM-PHD. Furthermore, on the basis of CPHD, a discrete kernel estimator in conjunction with the exponential weighted moving average scheme is introduced to adjust the birth intensity [15]. In [16], the birth intensity is approximated by the birth particles drawn from each of the measurements with the same weight, which is applied into the Cardinality Balanced Multi-Target Multi-Bernoulli (CBMeMBer) filter. Although these methods effectively adjust the birth intensity according to the target birth situation, most of these lead to an increase in algorithm complexity, which makes it difficult to achieve a balance between adaptability and realtime performance.
Since the measurements generated by the intersected targets overlap each other, GM-PHD cannot distinguish these targets and inaccurately estimates the number of targets. To solve this problem, many improved GM-PHD filters are proposed. The Collaborative Penalized Gaussian Mixture PHD (CPGM-PHD) [17] filter uses different identities of each target to co-penalize the weights of the targets. However, the penalized weight schemes adopted in it only consider the distance between measurements and Gaussian components, and rely on the closest measurement too much. Based on the assumption that the measurements from multiple closing targets can be considered as multiple measurements from a single extended target, the closing target tracking method based on random matrices is proposed in [18,19]. Meanwhile, the Fuzzy C-means (FCM) is adopted into the PHD filter to classify intersected targets [20]. However, these filters are premised on the known physical extension of targets, which limits their practical application. Recently, the Radiation Intensity GM-PHD (RIGM-PHD) [21] is proposed to track closing targets, which adds radiation information into the GM-PHD framework. However, RIGM-PHD is greatly influenced by noise.
To accurately estimate the state and number of the birth targets and intersected targets in a dense cluttered environment, this paper proposes a novel GM-PHD filter with a birth intensity and spawned intensity updating method based on the trajectory situation feedback. The innovations of this work can be further summarized as follows: (1) The target trajectories are formed on the basis of the discrete state estimation of GM-PHD, which is conducive to further analysis of the target's behavior in the scenario. (2) In the process of the proposed GM-PHD filtering, the birth intensity is adjusted adaptively according to the trajectory initiation of Gaussian components, which facilitates the accurate tracking of the unknown and time-varying birth targets. (3) Through the analysis and feedback of the trajectory situation, the proposed GM-PHD filter improves the identification of the intersected targets by enhancing spawned intensity, and maintains the target's state during the intersection period.
In order to verify the effectiveness of the proposed method, two improved methods based on GM-PHD, e.g., TIB-GM-PHD and CPGM-PHD, are selected for simulation comparison. The simulation experimental results show that the performance of the proposed algorithm is better than TIB-GM-PHD, CPGM-PHD and GM-PHD filters in both the birth scenario and intersection scenario.

Traditional PHD Filter
According to the RFS theory, the measurement sets and the states of targets at time k can be represented as , respectively. The prediction step and the update step are the basic processes of the PHD filter, and they can be denoted as: where p D,k , p s,k , f k|k−1 (x|ζ) and g k (.) are the detection probability, survival probability, transition probability density of target states and single target likelihood, respectively; κ k (z), γ k (x), β k|k−1 (•|ζ) are the clutter intensity, birth intensity and spawned intensity, respectively.

GM-PHD Filter
In the linear Gaussian dynamic system, the closed solution of PHD filter can be calculated by Gaussian mixture model. The prior intensity of multi-targets in the prediction step is expressed as: where ω k−1 are the weight, mean and covariance of the i-th Gaussian component at time k − 1, and the birth intensity γ k (x) and the spawned intensity β k|k−1 (x) are, respectively, denoted as: where J γ,k , ω γ,k describes the state of the birth target; similarly, J β,k , ω β,k−1 is a model determining the state of the spawning target with the previous state ζ. γ k (x) and β k|k−1 (x) determining the ability of the filter to detect the birth and spawned targets, respectively, and these model parameters are considered as fixed prior information in the filtering process.
Then the predicted intensity at time k based on Equation (1) can be denoted as: where v S,k|k−1 (x), v β,k|k−1 (x), γ k (x) is the predicted intensity of the survived target, the spawned target and the birth target, respectively. Afterward, the posterior intensity v k (x) based on Equation (2) in the update step is derived with the measurement sets Z k , and it can be calculated as: where The posterior intensity is estimated by means of assigning different weights to the Gaussian components in a recursive formulation at every scan. Afterward, the components that have weak weights will be pruned and the Gaussian components with close distance are merged into one. Note that, after the component pruning and merging procedure, the components with high weight can finally be extracted as target states.
It should be mentioned that the conventional filter sets the birth intensity γ k (x) and spawned intensity β k|k−1 (x) to fixed values, which results in the inability to perceive the target variation adaptively in the scene, especially in the birth scenario or intersection scenario [22,23]. There is a common phenomenon that the GM-PHD filter is unable to extract the state of the birth target during the initial time or extract the state of the intersected targets during the intersecting time accurately. Figure 1 shows the simulation scenario where five targets are considered.   Figure 2 shows the position and cardinality estimation for GM-PHD, and red circles represent the estimation of GM-PHD. According to Figure 2, target 2, 4 and 5, whose initial positions are not at the target birth positions, have not been detected successfully at the beginning. Furthermore, due to the fact that the measurements from the intersected targets overlap during intersecting time, the number of targets is estimated inaccurately when target 1 and 3 cross and target 2 and 3 cross. In the birth scenario, once the appearance of the target does not obey the description of birth intensity, the Gaussian components generated by measurements from the birth targets will obtain low weights and be filtered as clutter. Section 3.1 proposes the adaptive birth intensity updating method based on the trajectory initiation of the Gaussian components with low weights. Furthermore, in the case where overlapping measurements from intersected targets cannot be recognized, the filter can regard multiple intersected targets as the spawning of a single target according to the relative motion trend between targets, so as to maintain the number of targets. Section 3.2 proposes the spawned intensity updating method based on trajectory situation feedback.

GM-PHD Filter Based on Trajectory Situation Feedback Control
The flowchart of the proposed approach is shown in Figure 3. The feedback layer is appended on the basis of the traditional GM-PHD filter, and the discrete and disordered target states are formed into continuous trajectories [24] in the feedback layer. According to the trajectory situation analysis, i.e., birth trajectory situation or intersecting trajectory situation, the feedback layer adaptively adjusts birth intensity or spawned intensity at the next scan.  The Gaussian component generated by the birth target normally obtains a low weight when it appears at the first time [25][26][27]. Therefore, the component is regarded as clutter and is pruned in the state extraction procedure. To solve this problem, the Gaussian components after the component pruning and merging procedure [20] are recorded and the rule-based trajectory initiation technique is applied. Once the birth trajectory initializes successfully, the birth intensity will be automatically updated at the next scan.
When multiple targets are close to each other or even overlap, the Gaussian components generated by them are viewed as single components in the prediction and update procedure [20] and are stateless, which results in inaccurate estimation of the intersected target numbers and states. Hence, the proposed algorithm evaluates intersecting trajectory situation after state extraction procedure, thereby increasing the weights of components in the intersecting area by the means of enhancing the spawned intensity of targets in the intersecting area at the next scan.

Birth Intensity Estimation Method
Since the target state is continuous at multiple scans, the rule-based trajectory initiation can be applied in updating birth intensity. Based on this technique, as soon as the Gaussian components after the component pruning procedure at n consecutive scans are correlated successfully, the birth trajectory is considered to be initialized and the birth intensity γ k (x) will be updated according to the current merged Gaussian component. v p,k represent n Gaussian components after the component pruning and merging procedure at time k − n + 1, and k − n + 2, . . ., k, respectively, L represents a certain component in the posterior intensity. In order to indicate whether v p,k−n+1 , the deviation distance about them is expressed as: where r min , r max and t are the minimum rate of status change, maximum rate of status change and interval time, and d k−n+1,k−n+2 denotes the deviation distance from the normal rate region. Once the averaged rate between v p,k−n+2 falls in the normal region [r min , r max ], they are considered related and d k−n+1,k−n+2 equals the closest distance to the normal region.
Generally, the covariance of the zero-mean Gaussian white noise is R, and the normalized square is denoted as: k−n+1,k−n+2 is random variables with χ 2 distribution and the degree of freedom is p. In consequence, the deviation threshold ε can be determined in the χ 2 distribution table.
If the merged Gaussian components v < ε, the birth intensity composed of the above n components will be confirmed. Based on the least square technique, the intensity of birth target is estimated as: where, ω γ * ,k is the weight of birth target, and where m and P are the mean and covariance of the Gaussian components.
In summary, the birth trajectory can be initialized by repeating the above processes, and the birth intensity containing the intensity of all birth trajectories at the current scan is rewritten as: where J γ * ,k is the number of birth targets. The implementation details of the birth intensity estimation method are summarized by the pseudo-code in Algorithm 1.
Let N and M denote the respective number of targets and measurements. N γ and N β are the number of birth targets and spawned targets, respectively, and N g is the number of Gaussian components after the component pruning and merging procedure. As is known to all, the computational complexity of GM-PHD is O(M(N + N γ + NN β )). Compared with GM-PHD, TIB-GM-PHD requires additionally calculating the association between the measurements except from surviving targets to initiate birth trajectory, thus its computational complexity is approximately described as O(M(N + N γ + NN β + M)). In contrast, instead of using measurements, the proposed birth intensity estimation method initiates birth trajectory according to the components after the pruning and merging procedure. Since the trajectory initiation as described in Equations (10) and (11) requires computing associations between components at n consecutive scans, the computational complexity of the proposed method is O(M(N + NN β ) + nN g 2 ). Generally, M is much greater than N g . Hence, the proposed method does not have much negative impact on the real-time property of tracking.
Algorithm 1 Pseudo-code for the birth intensity estimation.
, minimum rate of status change r min , maximum rate of status change r max , interval time t and deviation threshold ε. Set: I = ∅, γ k+1 = ∅ and m = 0.

Spawned Intensity Estimation Method
When targets are close to each other, the filter will consider the overlapping components generated by multi-targets as a single component. For this reason, the intersecting trajectory situation should be fed back to the prediction and update procedure at the next scan. Once the Gaussian component is too close to the intersection area at the next scan, it is considered as the component formed by multiple targets, and the weight of this component is increased by each intersecting trajectory according to enhancing the spawned intensity. Ultimately, the weight of the component is approximately equal to the number of all intersecting trajectories.
Therefore, the algorithm concerned with the spawned intensity v β,k|k−1 (x) is used to decrease the impact of the intersection according to trajectory situation feedback. Let's e,k−1 ) as the components generated by targets after the state extraction procedure at time k − 1, and these targets in successive scans can form the corresponding trajectory η Each trajectory state can be represented by the current component in Gaussian form, and KL divergence [28] is used to express the relation between trajectories. The KL diver- k−1 is calculated as: where m In fact, the spawned intensity can be determined not only from its own trajectory but also from other trajectories, especially the neighboring trajectories. Therefore, the spawned intensity β k|k−1 (x|ζ) includes two parts, i.e., the original spawned intensity in Equation (5) and the enhanced spawned intensity calculated according to the closeness of other trajectories, and it is rewritten as: where J k−1 is the number of targets at time k − 1, and h Correspondingly, the predicted spawned intensity v β,k|k−1 (x) in Equation (6) will be updated as: where, m As discussed above, even though multiple components are mistaken for a single component due to coincidence, the weight of the single component will be enhanced by other trajectories and should be approximately equal to the number of intersection trajectories. The expected number of intersecting targets will be estimated accurately according to the weight of the component.
It is worth noting that the weight of the components that are close to each other but can still be distinguished will be much greater than 1 in the state extraction procedure, which leads to one target being estimated as multiple targets. The distinguishable component not only obtains a large weight given by the original trajectory but is also overly enhanced by other trajectories. Here, the component that has been pruned and merged needs to determine which trajectory it belongs to before the state extraction procedure, and the corresponding trajectory should be deleted as long as the target state is extracted once from the component.
Suppose v where m and P are the mean and covariance. Let I k−1 be the summary of all trajectories at time k − 1, then the trajectory η (l) k−1 to which v (i) p,k belongs meets the condition: where U is the divergence threshold between component and trajectory. For trajectory η (l) k−1 , the summary of trajectories adjacent to it can be denoted as: where T is trajectory closeness threshold. Generally, the number of targets should be less than or be equal to the number of trajectories in the intersection area. Therefore, as long as the state of the target is extracted once from the component v The implementation details of the spawned intensity estimation method are summarized by the pseudo-code in Algorithm 2.
Step 1. prediction for birth and existing targets (details omitted, as seen in Algorithm 1 and Ref. [10]).
Step 5. state extraction. 18: for i = 1, · · · , J k do 19: if ω if p = ∅ then 24: 25: for j = 1, · · · , round(ω Let N and M denote the respective number of targets and measurements. N γ and N β are the number of birth targets and spawned targets, respectively. The computational complexity of GM-PHD is known as O(M(N + N γ + NN β )). The CPGM-PHD utilizes the one-to-one correspondence between each measurement and each target to collaboratively penalize and renormalize the weight of targets, which almost doubles the computational complexity of GM-PHD, hence the complexity of it is described as O (2M(N + N γ + NN β )). For the proposed spawned intensity estimation method, the additional computation is required to calculate the closeness of trajectories and the enhancement of the spawned intensity, respectively, thus its computational complexity is O(M(N + N γ + NN β ) + N 2 (1 + M)).

Simulation Results and Discussion
This section presents simulation results from two scenarios to illustrate the performance of the proposed PHD comparing to GM-PHD (described in Ref. [10]), TIB-GM-PHD (described in Ref. [12]), CPGM-PHD (described in Ref. [17]). In two simulation experiments, the surveillance region in the two-dimensional scenario is [−500, 500] (m) × [−500, 500] (m), and the target state is [p x,k , p y,k ,ṗ x,k ,ṗ y,k ] T , where (p x,k , p y,k ) is the position and (ṗ x,k ,ṗ y,k ) is the velocity. For each target, its detected probability p D,k = 0.98, the survival probability p S,k = 0.99.
The transition probability density of the target follows linear Gaussian dynamics where I 2 and ⊗ denote the n × n identity and Kronecker product, respectively. ∆t = 1 s is the time step, and σ v = 5 (m/s 2 ) is the standard deviation of process noise. The clutter is modeled as a Poisson RFS with the intensity κ k (z) = λ c Vu(z). In the simulation, the clutter rate is λ c = 5.0 × 10 −5 m −2 and the area of the surveillance region is V = 1.0 × 10 6 m 2 . u(z) is the uniform density over the surveillance region. Further, if the length required to initiate the birth trajectory is too long it will increase the computational complexity, and, on the contrary, too short will lead to a high false alarm, thus the length is set to n = 3. Meanwhile, there may be bias in the acquisition of measurements from the sensor. To ensure that components and trajectories can be accurately correlated in the state estimation procedure, the values of the measurement-trajectory threshold U and the trajectory closeness threshold T are slightly large, and they are set to U = 5 and T = 0.8, respectively.
To effectively evaluate the performance of different PHD filters, 200 Monte Carlo runs are performed, and the optimal sub-pattern assignment (OSPA) criteria and cardinality estimate criteria [29] are adopted. The OSPA parameters are fixed at p = 1, c = 150.

Target Birth Scenario
In this scenario, there are ten birth targets with a random initial state in the twodimensional plane. The trajectories of the ten targets are shown in Figure 4, and the detailed parameters of targets are shown in Table 1.  As seen in Table 1, only target 1, 3, 4 and 9 are born at time 1 s, and at different locations. The rest of the targets are born one by one during the tracking process. In order to verify the tracking performance of the proposed PHD for unknown birth targets, the initial birth intensity does not include all initial birth target states as shown in Table 1. Specifically, the initial birth intensity γ k (x) and the spawned intensity β k|k−1 (x|ζ) are modeled as: where m  Figure 5. Figure 5a is the position estimation of GM-PHD. Figure 5b is the position estimation of TIB-GM-PHD. Figure 5c is the position estimation of CPGM-PHD. Figure 5d is the position estimation of the proposed PHD. According to Figure 5a,c, due to the lack of the adaptive-update of the birth intensity, the birth intensity of GM-PHD and CPGM-PHD is not in accordance with the state of the birth targets at time 1-18 s, which results in the loss of the birth target position estimation. Unlike them, it can be seen in Figure 5b that TIB-GM-PHD has good tracking performance for target 1, 3, 4 and 9. However, simple trajectory initiation conditions of TIB-GM-PHD lead to target 9 being underestimated at the beginning. Compared with TIB-GM-PHD, the proposed PHD considers prior information at multiple scans, which brings the better tracking performance of birth targets, as shown in Figure 5d.   Figures 6 and 7 show the OSPA distance and cardinality estimate for GM-PHD, TIB-GM-PHD, CPGM-PHD and the proposed PHD filter, respectively. During the period 1-18 s, the OSPA and cardinality estimate demonstrates that the proposed PHD has the best performance in the estimation of the birth target's state and number compared with the other three PHDs. In Figure 6, the peaks at around 15 s, 25 s, 35 s and 65 s indicate that the birth targets appear. For these targets born in the tracking process, the proposed PHD also has the lowest error rate in the state estimation. Meanwhile, Figure 7 shows that the cardinality estimate curve of the proposed PHD is closer to the truth curve with the help of the trajectory initiation feedback.  Considering that the clutter rate λ c and the trajectory initiation length n have a great influence on the performance of the proposed PHD, the effects of λ c on the different PHDs and the effects of n on the proposed PHD are verified in the scenario 1, respectively. Figure 8 illustrates the averaged OSPA distance and the sum of the cardinality estimate error of four PHDs with different clutter rates, over 200 Monte Carlo runs. The clutter rate varies from 1.0 × 10 −6 m −2 to 10.0 × 10 −6 m −2 . As shown in Figure 8, the estimation errors of four PHDs for the target state and number increase with the clutter rate. However, the performances of GM-PHD and CPGM-PHD are at the same level and worse than TIB-GM-PHD and the proposed PHD. Moreover, since the initial conditions of the proposed PHD involve more prior information, it estimates the target state and number more accurately.  Figure 9 illustrates the averaged OSPA distance and the sum of the cardinality estimate error of the proposed PHD with different trajectory initiation length, over 200 Monte Carlo runs. In Figure 9, there are three performance curves of the proposed PHD with different detection probabilities p D . It was observed that the tracking performance of the proposed PHD improves with the increase in detection probability. When the detection probabilities are 0.90 and 0.95, the best performances are located at n = 3 and n = 4, respectively. This is because the initiation conditions for the birth targets become more stringent as the trajectory initiation length grows. Once the measurements from birth targets at a certain scan are lost, the birth target fails to be detected.

Target Intersection Scenario
In this scenario, there are five targets in the two-dimensional plane. The trajectories of five targets are shown as Figure 10, and the detailed parameters of the targets are shown in Table 2.  As seen in Figure 10 and Table 2, to verify the effective tracking of the intersected targets, target 1 and 3, target 2 and 4 are born at the same time but at two different locations, and they will meet at some point. During the movement of targets, target 2 and 4 meet at 28 s, target 1 and 3 meet at 65 s, and target 2 and 3 meet at 89 s. Specifically, the initial birth intensity γ k (x) and the spawned intensity β k|k−1 (x|ζ) are modelled as: where m  Figure 11. Figure 11a is the position estimation of GM-PHD. Figure 11b is the position estimation of TIB-GM-PHD. Figure 11c is the position estimation of CPGM-PHD. Figure 11d is the position estimation of the proposed PHD. From Figure 11a-d, it is obvious that GM-PHD has the worst performance in the complex scenario, regardless of the birth scenario or intersection scenario. Figure 11b shows that TIB-GM-PHD successfully tracks birth target 5 in a timely manner but fails to maintain the position estimation for intersected targets.
By contrast, CPGM-PHD can distinguish the intersected targets but has poor performance in the tracking of the birth targets, as shown in Figure 11c. From Figure 11d, compared with these three filters, the proposed PHD stably estimates the position of the birth and intersected targets due to the model parameters updating according to the trajectory situation feedback.   Figures 12 and 13 show the OSPA distance and cardinality estimate for GM-PHD, TIB-GM-PHD, CPGM-PHD and the proposed PHD filter, respectively. In Figure 12, the periods at around 23 s-38 s, 60 s-70 s and 85 s-94 s represent the intersection of targets and the peak at around 25 s indicates the birth of target 5. During the intersection period, GM-PHD and TIB-GM-PHD cannot distinguish measurements from the intersected targets, which leads to a rapid increase in OSPA distance. Moreover, the performance of the proposed PHD is better than CPGM-PHD on the basis of the trajectory analysis. Furthermore, at 25 s the birth target appears, the OSPA distance of all the filters rises rapidly except TIB-GM-PHD and the proposed filter. According to Figure 13, the cardinality estimation of all filters are lower than the truth during intersection periods. However, the cardinality estimate curve of the proposed PHD is closer to the truth curve.   Figure 14 illustrates the averaged OSPA distance and the sum of the cardinality estimate error of four PHDs with different clutter rates, over 200 Monte Carlo runs. The clutter rate varies from 1.0 × 10 −6 m −2 to 10.0 × 10 −6 m −2 . As seen in Figure 14, the tracking performance of four PHDs for the intersected targets decreases as the clutter rate increases. At the period of high clutter rate, the TIB-GM-PHD performs the worst because of false initiation of dense clutters as birth targets. In contrast, the proposed PHD has the best performance in the estimation of the target state and number.  Figure 15 illustrates the averaged OSPA distance and the sum of cardinality estimate error of the proposed PHD with different trajectory initiation length, over 200 Monte Carlo runs. In Figure 15, there are three performance curves of the proposed PHD with different detection probabilities p D . It suggests that the trajectory initiation length is not significant in the estimation of the intersected target. Moreover, the increase in detection probability can improve the tracking performance of the intersected target.  Further, in order to evaluate the proposed approach more intuitively, the averaged OSPA distance, the sum of the cardinality estimates error and the averaged running time per 100 samples (with MATLAB 2016b, Intel i5 1.8 GHz processor with 8G RAM) between filters are compared and are given in Table 3. Compared with GM-PHD, TIB-GM-PHD and CPGM-PHD, the performance of the proposed filter in OSPA distance improves 42.36%, 38.00%, and 18.74%, respectively; and improves 58.70%, 50.54% and 24.17%, respectively, in cardinality error. However, due to the additional calculation about the track initiation and the trajectories situation analysis, the averaged running time of proposed filter is slightly higher than GM-PHD and TIB-GM-PHD but is still lower than CPGM-PHD significantly.

Conclusions
This paper proposes a novel GM-PHD filter with a birth intensity and spawned intensity updating method based on the trajectory situation feedback. Firstly, considering the appearance of the birth target in the scene, the proposed method adjusts the birth intensity in time according to the trajectory initiation formed by the rule-based correlation of the Gaussian components, which facilitates the accurate tracking of the birth targets. In addition, by adding the feedback of the trajectory situation analysis, the proposed method perceives the intersection of targets in time, and improves the identification of the intersected targets through enhancing spawned intensity. Experimental results show that the proposed method has superior tracking performance in the scenario where there are both birth targets and intersected targets. However, the cost of trajectory initiation and additional analysis of the trajectories leads to the increased computational complexity in this algorithm. Moreover, for the situation where the target dies at the intersection period, since the estimation of the state and number of targets relies on the feedback of the previous state, the proposed approach cannot perceive the death of the target and maintains this error during the subsequent tracking.