Cubature Information SMC-PHD for Multi-Target Tracking

In multi-target tracking, the key problem lies in estimating the number and states of individual targets, in which the challenge is the time-varying multi-target numbers and states. Recently, several multi-target tracking approaches, based on the sequential Monte Carlo probability hypothesis density (SMC-PHD) filter, have been presented to solve such a problem. However, most of these approaches select the transition density as the importance sampling (IS) function, which is inefficient in a nonlinear scenario. To enhance the performance of the conventional SMC-PHD filter, we propose in this paper two approaches using the cubature information filter (CIF) for multi-target tracking. More specifically, we first apply the posterior intensity as the IS function. Then, we propose to utilize the CIF algorithm with a gating method to calculate the IS function, namely CISMC-PHD approach. Meanwhile, a fast implementation of the CISMC-PHD approach is proposed, which clusters the particles into several groups according to the Gaussian mixture components. With the constructed components, the IS function is approximated instead of particles. As a result, the computational complexity of the CISMC-PHD approach can be significantly reduced. The simulation results demonstrate the effectiveness of our approaches.


Background
Multi-target tracking refers to sequential approximation of the states (positions, velocities, etc.), and the number of individual targets. It has been widely used in ground-moving-target tracking [1], visual tracking [2], and distribution fusion [3]. In multi-target tracking, both the state and observation sets of targets are time-varying. In practice, the associations between state and observation sets are always unknown, thus posing a challenge for multi-target tracking. The conventional approaches, such as the nearest-neighbour Kalman filter (NNKF) [4], Extended Kalman Filter (EKF) [5,6], multiple hypothesis tracking (MHT) [7], and joint probabilistic data association (JPDA) [8], are used to formulate the explicit associations between states and observations of targets. However, caused by the targets appearing and disappearing, the state and observation sets of the multi-target are uncertain. Such uncertainty costs high complexity in the conventional approaches on constructing the association

Our Work and Contributions
In this paper, we propose a cubature information SMC-PHD (CISMC-PHD) and its fast implementation (F-CISMC-PHD) approaches, which can be used to estimate the time-varying number and states of multi-target. As aforementioned, the disadvantage of conventional SMC-PHD filters is the tracking inefficiency in nonlinear scenario. To avoid such inefficiency, our CISMC-PHD approach applies the posterior intensity as the IS function. With such a selection, the current observations can be incorporated into the IS function design. Then, we utilize the cubature information filter (CIF) [26] with a gating method to calculate the IS function. Benefitting from tracking accuracy of CIF in high dimensional nonlinear case, the CISMC-PHD approach is capable of estimating the time-varying number and states of targets. To avoid initializing birth intensity in whole state space, a birth intensity initialization method is proposed for our CISMC-PHD approach. At last, we present the F-CISMC-PHD approach to reduce the computational complexity by considering groups of particles as Gaussian mixture components. These components are applied to approximate the IS functions instead of particles of the CISMC-PHD approach. Since the number of Gaussian components is much less than that of particles, the computational complexity can be greatly reduced. The main contributions of our work are listed as follows.
(1) We propose a novel IS function approximating method, which utilizes the CIF with a gating method to enhance the estimation accuracy of the SMC-PHD filter. Specifically, first, the posterior intensity is applied as the IS function of our CISMC-PHD approach, in order to incorporate the current observation set into the IS function approximation. Then, the gating method is integrated into the update step of the CIF for approximating the IS function. Benefitting from the most recent success of CIF in nonlinear state estimation, the tracking performance of the proposed CISMC-PHD approach is significantly enhanced. (2) We develop a method to initialize the birth intensity for the next tracking recursion. Since the current estimated targets (i.e., current survival targets) are not possible to be the birth targets at the next recursion, the observations of estimated targets are removed from the current observation set. Then, using an unbiased model, the remaining observations are mapped to state space for the birth intensity initialization. As such, the birth intensity can be adaptively initialized, making the target tracking more accurate and stable. (3) We develop a fast version of the CISMC-PHD approach (namely F-CISMC-PHD). We first consider each group of particles as a Gaussian mixture component. Then these components are used to approximate the IS functions of the CISMC-PHD approach. With the approximated IS functions, the particles can be sampled from these components for the intensity prediction and update steps. As a result, the computational complexity of the proposed CISMC-PHD approach can be significantly reduced.
The rest of this paper is organized as follows. In Section 2, a brief overview of the SMC-PHD filter is provided. Section 3 proposes our CISMC-PHD approach. A fast implementation of CISMC-PHD approach is presented in Section 4. Simulation results are demonstrated in Section 5, and Section 6 concludes this paper.

A Brief Overview of The SMC-PHD Filter
In this section, we review the basic idea of the SMC-PHD filter in detail. The main notations used in this section are defined as follows.
x k The state of a dynamic target at time k z k The observation of a dynamic target at time k The intensity of birth target at time k L k The number of the survival particles at time k J k The number of the birth particles at time k p s (·) The survival probability of target The detected probability of target π(·, ·) The IS function of birth intensity q(·|·) The IS function of survival intensity The SMC-PHD filter, motivated by the particle filter, is a sequential implementation of the PHD filter. In the SMC-PHD filter, the posterior intensity can be represented by a set of random samples of state vector x k with associated weights, which are usually called particles. By substituting these particles into the recursion of the PHD filter, the multi-dimensional integrals can be replaced by summations of the particles, which is computationally tractable.
More specifically, we define the particle set at time k − 1 as {x k−1|k−1 are the state and weight of the i-th particle at time k − 1, respectively. The posterior intensity at time k − 1 can be modelled by where Z 1:k−1 is the multi-target observations from time 1 to k − 1, and δ(·) is the Dirac Delta function. Notice that Z k = {z k,1 , z k,2 , . . . , z k,M }. Given D k−1|k−1 (x k−1 |Z 1:k−1 ), the implementation of the SMC-PHD filter consists of the prediction and update stages.

Prediction:
We first denote IS functions for the survival and birth targets as q(x k , Z k ), respectively. Then, according to [18], the predicted intensity can be formulated by these IS functions, In Equation (4), J k is calculated by J k = ν γ(x)dx, where ν is the particle number of each birth target. Equations (2)-(4) can be then used to predict the states and weights of the particles.
Update: In this stage, we obtain the posterior intensity by updating Equation (2). Then, we have the following update strategy, where w (i) κ(·) denotes the clutter intensity, and Equations (2)-(6) include the main procedure of SMC-PHD at one recursion. Commonly, to avoid the degeneracy of particles, the resampling strategy is utilized to resample particle set {x . After resampling, we can use the clustering methods to extract number and states of targets.

The CISMC-PHD Approach for Multi-Target Tracking
In this section, we present our CISMC-PHD approach. To be more specific, we propose a novel IS function approximation algorithm incorporating the CIF and a gating method in Section 3.1. Then, Section 3.2 develops a method to initialize the birth intensity. Finally, Section 3.3 introduces the state extraction method for state estimation.
The nonlinear dynamic model of the target with state x k at time k is given as follows, Process model: Observation model: where φ(·) is the state transition function, and ϕ(·) denotes the relationship between state and observation. v k−1 and w k are the process and observation noises at time k − 1 and k, respectively. Both v k−1 and w k are assumed to be Gaussian noises with zero means, and their covariances are denoted as Q k−1 and R k . According to this model, the transition density f (x k |x k−1 ) and likelihood g k (z k |x k ) are subject to Gaussian distributions.

The IS Function Approximation Algorithm
As mentioned in Section 2, most of the conventional SMC-PHD filters utilize the transition density function as the IS functions, resulting in great tracking error for targets with nonlinear dynamics. A novel IS function approximation algorithm, incorporating the CIF with a gating method, is presented to improve the tracking accuracy.
In our approach, we select the IS functions of Equations (3) and (4) as Birth IS: where m k . Now, we discuss on how to calculate them. For simplicity, they are replaced by m k and P k , respectively. Here we use the CIF and gating methods to estimate them.
Before introducing the CIF method, we review the cubature rules [27]. The cubature rules are used to approximate the Gaussian weight integral. Assuming c(x) is a function on the n-dimension R n , its Gaussian weight integral can be approximated by where and [1] j is the j-th vector of the set According to Equation (12), the cubature rules can be used to compute the multi-dimension integrals in the prediction and update steps of the CIF method.
Prediction: In this step, we first predict the state m k|k−1 and covariance P k|k−1 according to cubature rules. Then, the predicted information state vector y k|k−1 and matrix Y k|k−1 are estimated for the update step.
Let m k−1 and P k−1 be the previous state and covariance, respectively. According to Equations (12) and (13), the j-th cubature point χ k−1,j can be estimated by Then, we can calculate m k|k−1 and P k|k−1 using the following formulations: where (·) T is the transpose operator, and Given Equations (15) and (16), the information forms of m k|k−1 and P k|k−1 are represented [28] by and where y k|k−1 and Y k|k−1 are the information state and matrix, respectively. Update: We use the observation set Z k to update the predicted y k|k−1 and Y k|k−1 in the current step. In order to construct the associations between the observation set Z k and predicted observation z k|k−1 , a gating method is applied to extract the associated observations. With the extracted observations, we can finally obtain m k and covariance P k .
We denote z k|k−1 as the predicted observation, computed by where and Utilizing the predicted observation z k|k−1 of Equation (20), the error cross covariance matrix of state and observation can be evaluated by With the above obtained parameters, we can calculate the state contribution and its corresponding information matrix as and where µ j is the innovation of the j-th observation z j (z j ∈ Z k ), expressed by In our scenario, z j is a two-dimension vector, and µ j follows a two-dimension Gaussian distribution.
In practice, the observation set may contain large clutters. The existence of these clutters cannot only degenerate the estimation accuracy, but also increase the computational complexity. Recently, several gating technologies have been proposed to remove the clutters from the observation set [29,30]. Inspired by [30], we utilize the gating technology to reduce the influence of clutters.
Intuitively, observations far away from the predicted observation are subject to be generated by clutters. These observations must be removed from the observation set. With the gating technology, the left observations can be represented bŷ where P zz k|k−1 is the covariance matrix of the predicted observation z k|k−1 , and (·) −1 is the matrix inversion. T h is the threshold of the gate. According to Equation (27), the innovation µ j follows the Chi-square distribution. Thus, T h can be determined by the dimension of µ j and association probability. Commonly, the square root of T h is known as the number of Sigma. Literature [30] proved that the number of Sigma gates ranging from 3 to 5 (corresponding to T h = 9 − 25) can guarantee the true observation lying inside the gate with "enough" probability (≥ 0.971), when the dimension of µ j is less than three. In this paper, we select the number of Sigma gates being to 4 (corresponding to T h = 16). When the dimension of µ j is less than three, such a selection guarantee that the association probability ≥ 0.998.
Then, we concentrate on computing P zz k|k−1 of Equation (27). Let P mz k|k−1 be the cross covariance matrix between observation and state space. According to the linear error propagating of [31], P mz k|k−1 can be rewritten as P mz where H k is the linearized matrix.
Obviously, H k can be approximated by With the achieved H k , according to [32], P zz k|k−1 can be calculated by Substituting the achieve P zz k|k−1 into Equation (27),Ẑ k can be extracted from the current observation set Z k .
With the extracted observation setẐ k , the information state vector y k and matrix Y k are represented as: Given information state y k and matrix Y k , posterior state m k can be reconstructed based on Equation (18) : Moreover, the posterior covariance P k is recovered based on Equation (19): If there is no observation that lies inside the gate (Ẑ k = ∅), we approximate m k and P k in the following, Substituting the above obtained m k and P k into Equations (10) and (11), we can approximate the IS functions of survival and birth targets for our CISMC-PHD approach.

The Birth Intensity Initialization Method
According to Equation (2), the birth intensity has large effect on the posterior intensity estimation. Targets may "born at anywhere" of the state space. In other words, birth intensity γ(x) may cover the whole state space, which is rather exhaustive [25]. To avoid such a disadvantage, observation-driven birth intensity initiation methods were proposed [20,33,34]. Inspired by these methods, an adaptive birth intensity initialization method is proposed for the CISMC-PHD approach. Instead of initializing birth intensity across the whole state space, the proposed method of CISMC-PHD approach utilizes the current observations and estimated targets to initialize the birth intensity at the next recursion. Compared with the conventional SMC-PHD filters, our method can initialize the birth intensity without knowing it as a prior.
The implementation of our method consists of two steps. First, in order to initialize the birth intensity, we remove observations generated by the estimated targets That is because the current survival targets cannot be new-born targets at the next recursion. Second, we use the remaining observations to estimate the birth target components, which can be used to calculate the birth intensity. With these two steps, the birth intensity can be initialized for the next recursion.
Step1. Remove observations generated by the estimated targets.
In the basic PHD filter, it is assumed that each target can yield at most one observation [35]. According to this assumption, each target has one and only one corresponding observation. Influenced by the noises and clutters, the observation generated by the target may appear around the target. In other words, observations around the target has the large probability to be generated by the same target. Therefore, the birth target state set can be estimated by removing states of estimated targets from the multi-target state.
Here, we adopt the bearing and range tracking model [36] to illustrate the birth intensity initialization method of our CISMC-PHD approach. Let x e k,i be the state of the i-th target in the estimate state set X e k . x e k,i consists of position and velocity, while z k,j (z k,j ∈ Z k ) consists of the bearing angle and range. We define the distance between x e k,i and z k,j (z k,j ∈ Z k ) as where (·) r denotes the range-dimension element, and | · | is the absolution value.
We follow the way of [30] to select the certain threshold for Equation (37), where σ r is the error of the range-dimension (known as a prior). l is the confidence level, commonly selected from l = 3 ∼ 5. Here, we use l = 3, which can guarantee that the associated probability equals to 0.997. With Equations (37) and (38), we can remove the observations associated with the estimated targets. LetZ k be the observations of birth targets, the removing procedure is summarized in Table 1.
Compute the distance d i,j between x e k,i (x e k,i ∈ X e k ) and z k,j for each z k,i ∈ Z k by Equation (37).
Step2. Estimate the birth target components.
OnceZ k is obtained, we turn to estimate the birth target components (the mean of the i-th target state vector m (i) k,b and its corresponding covariance P (i) k,b ) by the unbiased model of [37]. Letz k,i ∈Z k , we mapz k,i into state space denoted byz c k,i = [p x k,i , p y k,i ] T . p x k,i and p y k,i can be computed by p x k,i = β −1 θ r k,i cos θ k,i and p y k,i = β −1 θ r k,i sin θ k,i . β θ = σ θ is a biased comparison factor, where σ θ , as a prior, is the error of bearing angle θ k,i . According to [37], m (i) k , the mean of the i-th birth target state, can be estimated as The covariance can be approximated by where σ v , as a prior, is the standard deviation of velocity. In Equation (40), the following exists, Finally, we can construct the new-born targets as {m k,i , P k,i } , where N k,b = |Z c k | is considered as the number of birth targets. We use Equation (11) to sample states of birth particles. The weights of these birth particles are initialized with the same values, w , where N is the number particles for each target, and γ(x) is defined in Section 2. On this basis, these birth particles become survival particles at time k + 1. That is to say, the IS functions of these particles at time k + 1 can be computed by the method of Section 3.1. Notice that the new-born target in this section may contain clutters, and these clutters can be removed in the resampling step of the CISMC-PHD approach.
The proposed initialization method may cause overestimation of targets. To overcome the issue of overestimation, some advanced methods, such as [33,34], etc., may be incorporated for initialization of our approach. It is an interesting future work.

State Estimation
In multi-target tracking, it is rather important to estimate the target number and states. As for the state estimation, clustering methods, are commonly used in SMC-PHD filters [16,18]. However, they are subject to biased estimation [21]. Ristic et al. [21] proposed an method that clusters the particles into several groups at the update stage. In this paper, we intend to adopt the method of [38] for state estimation, which is an improved method of [21]. There are also several alternative methods, such as Zhao's method [39] and MEAP method [40], which have better estimation performance.
According to Equation (6), the updated weight w i k of the i-th particle consists of two parts, In Equation (42), j = 0 denotes that there is no observation, and C z k,j can be computed by Equation (7). For state estimation, we aggregate W j of particle weights corresponding to observation z k,j , According to Equations (42) and (43), if z k,j is generated by the clutter, then the likelihood g k (z k,j |x (i) k ) may be small, leading to low value of W j . However, if z k,j is generated by the target, then W j may be large due to the large value of g k (z k,j |x (i) k ). Thus, setting certain threshold W th for W j , we can assign particles {x k } that satisfy W j > W th to the j-th target. In this paper, we set W th = 0.5, the same as [21]. Then, x k,j and P k,j can be calculated in the following Given Equations (44) and (45), the states of multi-target can be finally output. We summarize our approach at time k in Table 2.
, and current observation set Z k -Output: Target number M k , and estimated state set X e k 1 Calculate the IS function q(x (10) by Equations (33) and (34), and draw particles {x for survival targets by Equations (3) and (4), where m (10) using Equations (33) and (34), and draw particle set {x for birth particles by Equations (3) and (4), where m Calculate w (i) k|k by Equation (6) for resampling and w (i,j) k|k by Equation (42) for estimating, using the particle set {x Compute W j by Equation (43) with w (i,j) k|k calculated by step 3 and assign particles into the corresponding group by W th to estimate the target states and number.

5
Estimate the state set X e k by Equation (44) k|k and x (i) k|k,j belong to the j−th group of step 4. In addition, target number can be approximated by M k =| X e k |. 6 Resample particles {x k|k ], and [·] denotes the nearest integer. 7 Remove the observations of survival targets using Table 1, and estimate the birth components by Equations (39) and (40) to obtain the birth component set B k .
8 Draw particle from Equation (11) to obtain the birth particles {x Return M k and X e k .

A Fast Approach For The CISMC-PHD Filter
In this section, we focus on reducing the computational complexity of our CISMC-PHD approach. Section 4.1 presents a fast implementation of the CISMC-PHD approach, namely the F-CISMC-PHD approach. Then, we analyze the computational complexity of the CISMC-PHD and F-CISMC-PHD approaches in Section 4.2. The framework of the improved approach is illustrated in Figure 1. Figure 1. Framework of the F-CISMC-PHD approach. In the IS function Approximating step, we utilize the survival and birth components to approximate the IS functions. Then the predicted particles are generated and updated in the PHD prediction and update step to achieve the posterior particles. By clustering the particles into several groups, the Gaussian Mixture components (namely survival components) can be constructed in the survival components construction step. Meanwhile, we also apply the survival components to estimate the birth components in the birth components construction step. These components are used to approximate the IS functions in the next iteration.

The Fast CISMC-PHD Filter
In this section, we introduce the F-CISMC-PHD approach, which can save the computational time of the CISMC-PHD approach. Inspired by [41], we consider the particle groups of targets as the Gaussian mixture components. Recall that N k−1,b is the number of birth components. m where m Commonly, birth components at time k − 1 become survival components at time k. We combine birth and survival components into one set. That is {G k−1 denote the weight, mean and covariance of the i-th combined component.
With the combined components, we use the CIF method of Section 3.1 to approximate the IS function of each component. Here, we utilize q(x|m (i) k−1 , Z k ) to represent the IS function of the i-th component. On this basis, the j-th predicted particle, which is generated by the i-th component, can be represented by k|k−1 and P (i) k|k−1 are the predicted mean and covariance of i-th component, respectively, which can be computed by Equations (15) and (16).
· denotes the nearest floor integer, and G (i) k−1 · N is the number of particles generated by the i-th component. According to Equations (47) and (48), the number of predicted particles is Then, the predicted weight of Equation (48) is substituted into Equation (42) for particle grouping and state estimation. Recall that, weight W j of Equation (43) can be used to assign the particles into the j-th group, where the mean x k,j and covariance P k,j can be calculated by Equations (44) and (45) in Section 3.3. Given W j , x k,j and P k,j , we model the j-th group as a Gaussian mixture component With such a selection, it can guarantee that the groups with large W j have enough number of sampling particles. Note that the group with small W j has the large probability to be generated by the clutter, and such a group should be neglected to avoid the waste of computational time. We set a certain threshold T g for {G In this paper, we set T g = 0.1. The construction procedure of the target components is summarized in Table 3.
The target components of Table 3 refer to the survival target components. The birth components and target state estimation can be achieved by Sections 3.2 and 3.3, respectively.  (44) to approximate the mean x k,N s,k and covariance P k,N s,k . − Save the weight of the current component as G Return the target components {G

Computational Complexity
In this section, we analyze the computational complexity of the CISMC-PHD, F-CISMC-PHD and conventional SMC-PHD approaches. For justice, we adopt the same state estimation and birth target initialization methods for the three approaches. The computational complexity of the three approaches on state estimating and birth target initializing is the same, when they have same particle numbers and observations. Thus, it can be neglected for the computational complexity comparing. In addition, we select the multinomial resampling method as the resampling methods of the CISMC-PHD and conventional SMC-PHD approaches. The particle numbers of these approaches are equal to N p .
We begin with the computational complexity analysis of the CISMC-PHD approach. Incorporating the CIF and gating methods into the SMC-PHD approach, the CISMC-PHD approach can achieve a good estimation accuracy in nonlinear target tracking. As mentioned in Section 3.1, each particle is applied to approximate the IS functions. The computational complexity of the CIF and gating method per particle is nearly O(n 3 d ), where n d is the dimension of the particle state. The computational complexity of the CIF-based IS functions is O(N p · n 3 d ). Besides, the computational complexity of the PHD update step is O(N p · M), where M is the number of observations. In addition, the resampling step of our CISMC-PHD approach consumes O(N p ) computational complexity. Thus, the computational complexity of our CISMC-PHD approach in total is O(N p · n 3 d + N p · M + N p ). Then, we turn to the analysis of the computational complexity of the F-CISMC-PHD approach. Since the F-CISMC-PHD approach adopts the target components to compute the CIF-based IS functions, the computational complexity of the CIF-based IS functions is O(N G · n 3 d ), where N G is the number of target components, N G N p . Assuming that in the PHD prediction step, N G target components generate N p particles. The computational complexity of the PHD update step is O(N p · M). Besides, the Gaussian target component forming consumes O(M + 1). Combining the computational complexity of these steps together, the computational complexity of the F-CISMC-PHD approach is O(N G · n 3 d + N p · M + M + 1). In practice, M N p , the computational complexity of F-CISMC-PHD filter is much less than the CISMC-PHD filter.
The conventional SMC-PHD approach uses the transitional density as the IS function. Compared with the CISMC-PHD and F-CISMC-PHD approaches, it does not need to the IS function computing. Hence, the computational complexity of the conventional SMC-PHD approach can be approximated as O(N p · M + N p ), where O(N p · M) and O(N p ) are the computational complexity of the update and resampling steps. Obviously, O(N p · M + N p ) < O(N G · n 3 d + N p · M + M + 1) < O(N p · n 3 d + N p · M + N p ). Thus, the computational complexity of the conventional SMC-PHD approach is smaller than the CISMC-PHD and F-CISMC-PHD approaches.
With the above discussion, we can conclude that the conventional SMC-PHD approach has the lowest computational complexity, and the CISMC-PHD approach has the highest computational complexity. However, the estimation accuracy of the conventional SMC-PHD approach is the lowest, and the estimation accuracy of the CISMC-PHD approach is the highest. The F-CISMC-PHD approach can make a trade-off between the computational complexity and estimation accuracy. Such a conclusion can be observed in Section 5.

Simulation Results
In this section, we validate the tracking performance of the proposed CISMC-PHD and F-CISMC-PHD approaches. In Section 5.1, a nonlinear simulation scenario composed of five targets is constructed. Then, Section 5.2 compares the estimation results of the IBSMC-PHD [20], proposed CISMC-PHD and F-CISMC-PHD approaches, in terms of the optimal subpattern assignment (OSPA) metric [42] and Root Mean Square Error (RMSE). At last, we compare the simulation results of all three approaches with different numbers of clutters and detection probabilities, to validate the effectiveness of our approaches.

Simulation Scenarios
In our simulations, we use the nonlinear scenario, the same as [36]. Let x be the target states, represented by x = [p x , v x , p y , v y , α] T . In this paper, (p x , p y ) is the position, (v x , v y ) is the velocity, and α is the turn rate. With the above definitions, we model the nonlinear dynamic equation as where T = 1 second (s) is the sampling interval. In addition, ε k−1 is the noise, defined by ε k−1 ∼ N(ε k−1 ; 0, Q). We denote Q = diag(σ 2 x,ε , σ 2 y,ε , σ 2 α ) as covariance matrices of ε k−1 . In this paper, we set σ x,ε = σ y,ε = 1 meter/second 2 (m/s 2 ), σ α = π/180 rad. The initial states, appearing times, and disappearing times of targets are listed in Table 4. Besides, the observation model is given by where η k is the observation noise defined by η k ∼ N(η k ; 0, R), and R = diag(σ 2 θ , σ 2 r ) is the covariance. Here, we set σ θ = π/180 rad and σ r = 1 meter (m). We also assume that clutters are uniformly distributed in the detection region, where the angle range is (0, π/2) rad, and distance range is (0, 1000) m. Trajectories of the five targets in both scenarios are shown in Figure 2, where the clutter number is set to be 10 for each scenario. For parameters, we set the gating threshold T h = 16 according to [32], the probability of detection and survival are p d (x k ) = 0.95, and p s (x k ) = 0.99, in accordance with [36]. The particle numbers for each birth target and survival targets are 5 and 100, respectively. All of the simulations are run in a computer with MATLAB 2015a, and i5 3.2 GHz processor with 4GB RAM.

Comparison of Estimation Accuracy on Certain Number of Clutters
To compare the estimation accuracy, we adopt the first order OSPA and RMSE as the metric. Here, we discuss the first order OSPA distance in the following. Let X = {x 1 , . . . , x n } and Y = {y 1 , . . . , y n } be two RFSs, where m and n are numbers of elements in X and Y, respectively. Supposing that Ω n represents the set of permutations of {1, 2, . . . , n}, the first order OSPA metric can be rewritten byd where d c (x, y) = min(c, d(x, y)), c > 0 is a cut-off factor, and d(x, y) is the distance between x and y.
In this paper, we set c = 150 in accordance with [43], and use the Euclidean distance to compute d(x, y). Then, we conducted 500 Monte Carlo runs for multi-target tracking with the IBSMC-PHD, our CISMC-PHD, and F-CISMC-PHD approaches. The estimated trajectories of all three approaches are demonstrated in Figure 3. We can observe that most of estimated points are covered with the true trajectories in Figures 3b,c, while most of points are not covered with the ground-truth trajectories. Furthermore, Figure 4 depicts the OSPA distances of all three approaches. In these figures, the OSPA distances of the proposed CISMC-PHD and F-CISMC-PHD approaches are smaller than the BSMC-PHD approach. Note that large OSPA distance denotes large tracking error. Thus, the proposed CISMC-PHD and F-CISMC-PHD approaches have the smaller tracking error than the IBSMC-PHD approach. We can also observe that the OSPA distances of the proposed CISMC-PHD and F-CSMC-PHD approaches have four peaks at time 2, 9, 26, and 60. That is to say, targets may appear at these times.  We also plot the estimated numbers and corresponding RMSEs of all three approaches in Figure 5. As seen from Figure 5a, the estimated number of the proposed CISMC-PHD approach is closest to the ground truth among three approaches, thus enjoying the lowest RMSE in Figure 5b. In addition, the numerical results of the averaged OSPA distances and RMSEs are listed in Table 5, demonstrating that our CISMC-PHD approach achieves the best estimation on numbers and states. According to Table 5, the F-CIMS-PHD approach can make a compromise between the computational time and estimation accuracy.

Comparison of Estimation Accuracy on Various Numbers of Clutters
To validate the influence of clutters on multi-target tracking, the IBSMC-PHD, CISMC-PHD and F-CISMC-PHD approaches were implemented with 500 Monte Carlo simulations alongside the clutter number changing from 1 to 30. The results are illustrated in Figure 6a,b. From this figure, we can see that the averaged OSPA distances of all three approaches are enhanced, when the clutter number increases from 1 to 30. Among these approaches, the CISMC-PHD approach has the smallest averaged OSPA distance and RMSE.

Comparison of Estimation Accuracy over Different Probabilities of Detection
In this section, we compare the estimation accuracy at different detection probabilities (varying from 0.92 to 0.98). Here, the clutter number is chosen to be 10, and 500 Monte Carlo simulations are run for the comparison. Table 6 reports the OSPA distances and RMSEs of the IBSMC-PHD, CISMC-PHD and F-CISMC-PHD approaches. Table 6. Tracking Performance over different detection probabilities.  From Table 6, we can observe that the estimated accuracy of the F-CISMC-PHD approach get close to the CISMC-PHD approach, when the probabilities of detection increase from 0.92 to 0.98. Thus, the F-CISMC-PHD approach is suitable for the high probabilities of detection.

Comparison of Estimation Accuracy at Challenging Nonlinear Scenarios
In this section, we compare the estimation accuracy of the IBSMC-PHD, CISMC-PHD and F-CISMC-PHD approaches at challenging nonlinear scenarios. For the challenging nonlinear scenarios, the standard deviation σ θ varies from 1.5π 180 to 3π 180 . We implement each approach with 500 Monte Carlo simulations.
The averaged accuracy, evaluated by OSPA and RMSE, is listed in Table 7. Table 7 indicates that the estimation accuracy of all three approaches decreases, when σ θ increases from 1.5π 180 to 3π 180 . Compared with the IBSMC-PHD approach, the OSPA distances and RMSEs of the CISMC-PHD and F-CISMC-PHD approaches are smaller at all four values of σ θ . It means that the estimation accuracy of the CISMC-PHD and F-CISMC-PHD approaches is more stable than the IBSMC-PHD approach in the challenging nonlinear scenarios.

Conclusions
In this paper, we have proposed the CISMC-PHD and F-CISMC-PHD approaches, which can estimate the time-varying number and states of multi-target nonlinear tracking. In our CISMC-PHD approach, a novel IS function approximation method is presented, which incorporates a gating method into the CIF method. To initiate the birth intensity of the next recursion, we use the current observations and estimated states to estimate the birth target components. In addition, we also present a fast implementation of the CISMC-PHD approach, namely F-CISMC-PHD, to reduce the time complexity of the CISMC-PHD approach. By clustering the particles into several groups, the target components can be obtained by representing the groups as Gaussian mixture components. Utilizing these components to approximating the IS functions, the computational time can be reduced magnificently. The simulation results demonstrate that the proposed CISMC-PHD and F-CISMC-PHD approaches outperform the conventional BSMC-PHD approach.
This paper concentrates on improving efficiency and accuracy of the conventional SMC-PHD filter. It simply utilizes the multinomial resampling method as the resampling method. Other resampling methods may be integrated in the CISMC-PHD approach for the future work. Besides, the F-CISMC-PHD approach is only suitable for the high probability of detection and small number of clutters. Study on improving the tracking performance of the F-CISMC-PHD approach at low probability of detection and large number of clutters may be seen as another direction of the future work.