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.

1. Introduction

1.1. 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 between them. Several approaches have emerged to improve the performance of conventional approaches in terms of tracking accuracy [9,10,11,12,13].

Recently, filters based on random finite sets (RFS) have been developed as alternative frameworks of the traditional algorithms [14,15,16,17] to estimate the multi-target number and states. In RFS formulations, both the multi-target state and observation sets are modelled as random finite sets. Given such models, the multi-target tracking problem is formulated under the Bayesian framework [14]. Since the optimal RFS-based Bayesian filter is computationally intractable, the probability hypothesis density (PHD) filter has been proposed as a sub-optimal Bayesian filter for multi-target tracking [15,16,17]. By propagating the first order moment (namely the intensity), the PHD filter can save much computational complexity in comparison to the optimal RFS-based Bayesian filter. Moreover, it avoids the combinatorial problem arising from data association, which is the bottleneck for conventional multi-target tracking approaches to estimate multi-target number and states.

However, the difficulty of the PHD filter is that it is intractable to derive close-form solutions to the PHD filter formulations. To solve such difficulties, Vo et al. [18] proposed the bootstrap SMC-PHD (BSMC-PHD) filter, which selects the transition density as the IS function. The modified versions of BSMC-PHD filter have emerged in [19,20,21]. As for those BSMC-PHD filters, the IS function selection simplifies the weight computation in the prediction stage of the SMC-PHD filter. In practice, such a selection leads to a few particles with large-valued weights in the update stage, when targets have nonlinear motion [22].

To improve the tracking performance of the conventional SMC-PHD filter, several efficient approaches have been presented. Morelande et al. proposed the Rao-Blackwellised SMC-PHD (RB-SMC-PHD) filter in [23]. The RB-SMC-PHD filter utilizes several auxiliary variables to define the IS function, thus making it suitable for the intensity approximation of the SMC-PHD filter with conditionally linear Gaussian models. Motivated by the auxiliary particle filter, Whiteley et al. proposed the auxiliary particle PHD (APHD) filter in [24], in which the auxiliary variable is pre-selected to minimize the variance of the IS weights. Although the APHD filter enhances the efficiency of the SMC-PHD filter in nonlinear scenario, it yields poor performance in case of severe nonlinearities or high process noise [22]. Inspired by the unscented particle filter, Yoon et al. [25] utilized the unscented information filter (UIF) to design the IS function of SMC-PHD (called USMC-PHD filter). Since such a design takes the current observations of targets into account, it is more stable than the BSMC-PHD filter. However, the drawback of this filter is that its performance has been influenced by the selection of the sigma-points of UIF. Ristic et al. [20] proposed a novel state estimation method (called IBSMC-PHD), rather than partitioning particles in ad-hoc manner. Their method groups the particles in the update stage, thus enjoying the computational efficiency.

1.2. 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.

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.

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.

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.

2. 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.

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 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 as , where and are the state and weight of the i-th particle at time , respectively. The posterior intensity at time can be modelled by where is the multi-target observations from time 1 to , and is the Dirac Delta function. Notice that . Given , 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 and , respectively. Then, according to [18], the predicted intensity can be formulated by these IS functions, where

In Equation (4), is calculated by , 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

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 . After resampling, we can use the clustering methods to extract number and states of targets.

3. 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 at time k is given as follows, where is the state transition function, and denotes the relationship between state and observation. and are the process and observation noises at time and k, respectively. Both and are assumed to be Gaussian noises with zero means, and their covariances are denoted as and . According to this model, the transition density and likelihood are subject to Gaussian distributions.

3.1. 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 where and are means of the survival and birth particles, respectively. and denotes the corresponding covariances of them.

Then, the problem of IS function design can be reduced to calculate and . Now, we discuss on how to calculate them. For simplicity, they are replaced by and , 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 is a function on the n-dimension , its Gaussian weight integral can be approximated by where and 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 and covariance according to cubature rules. Then, the predicted information state vector and matrix are estimated for the update step.

Let and be the previous state and covariance, respectively. According to Equations (12) and (13), the j-th cubature point can be estimated by

Then, we can calculate and using the following formulations: where is the transpose operator, and

Given Equations (15) and (16), the information forms of and are represented [28] by and where and are the information state and matrix, respectively.

Update: We use the observation set to update the predicted and in the current step. In order to construct the associations between the observation set and predicted observation , a gating method is applied to extract the associated observations. With the extracted observations, we can finally obtain and covariance .

We denote as the predicted observation, computed by where and

Utilizing the predicted observation 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 is the innovation of the j-th observation (), expressed by

In our scenario, is a two-dimension vector, and 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 by where is the covariance matrix of the predicted observation , and is the matrix inversion. is the threshold of the gate. According to Equation (27), the innovation follows the Chi-square distribution. Thus, can be determined by the dimension of and association probability. Commonly, the square root of is known as the number of Sigma. Literature [30] proved that the number of Sigma gates ranging from 3 to 5 (corresponding to ) can guarantee the true observation lying inside the gate with “enough” probability (), when the dimension of is less than three. In this paper, we select the number of Sigma gates being to 4 (corresponding to ). When the dimension of is less than three, such a selection guarantee that the association probability .

Then, we concentrate on computing of Equation (27). Let be the cross covariance matrix between observation and state space. According to the linear error propagating of [31], can be rewritten as where is the linearized matrix.

Obviously, can be approximated by

With the achieved , according to [32], can be calculated by

Substituting the achieve into Equation (27), can be extracted from the current observation set .

With the extracted observation set , the information state vector and matrix are represented as:

Given information state and matrix , posterior state can be reconstructed based on Equation (18) :

Moreover, the posterior covariance is recovered based on Equation (19):

If there is no observation that lies inside the gate (), we approximate and in the following,

Substituting the above obtained and into Equations (10) and (11), we can approximate the IS functions of survival and birth targets for our CISMC-PHD approach.

3.2. 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 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 be the state of the i-th target in the estimate state set . consists of position and velocity, while () consists of the bearing angle and range. We define the distance between and () as where 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 is the error of the range-dimension (known as a prior). l is the confidence level, commonly selected from . Here, we use , 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. Let be the observations of birth targets, the removing procedure is summarized in Table 1.

Table 1

Removing observations generated by the estimated targets.

– Input: Observation set Z_k, and estimated state set Xke – Output: Observation set of birth targets Z˜_k - For: j=1,2,⋯,|Xke|, 1 Compute the distance d_i,j between xk,ie (xk,ie∈Xke) and z_k,j for each z_k,i∈Z_k by Equation (37). 2 Extract z_k,i satisfying i={i|d_i,j<=3·σ_r}. 3 Remove z_k,i from Z_k to obtain Z˜_k. - End For - Return Z˜_k.

Step2. Estimate the birth target components.

Once is obtained, we turn to estimate the birth target components (the mean of the i-th target state vector and its corresponding covariance ) by the unbiased model of [37].

Let , we map into state space denoted by . and can be computed by and . is a biased comparison factor, where , as a prior, is the error of bearing angle . According to [37], , the mean of the i-th birth target state, can be estimated as

The covariance can be approximated by where , as a prior, is the standard deviation of velocity. In Equation (40), the following exists,

Finally, we can construct the new-born targets as , where 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, , where N is the number particles for each target, and is defined in Section 2. On this basis, these birth particles become survival particles at time . That is to say, the IS functions of these particles at time 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.

3.3. 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 of the i-th particle consists of two parts,

In Equation (42), denotes that there is no observation, and can be computed by Equation (7). For state estimation, we aggregate of particle weights corresponding to observation ,

According to Equations (42) and (43), if is generated by the clutter, then the likelihood may be small, leading to low value of . However, if is generated by the target, then may be large due to the large value of . Thus, setting certain threshold for , we can assign particles that satisfy to the j-th target. In this paper, we set , the same as [21]. Then, and 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.

Table 2

The CISMC-PHD filter at time k.

– Input: Birth particle set {xk-1,b(i),Pk-1,b(i),wk-1,b(i)}i=1J_k-1, survival particle set {xk-1,s(i),Pk-1,s(i),wk-1,s(i)}i=1L_k-1, and current observation set Z_k – Output: Target number M_k, and estimated state set Xke 1 Calculate the IS function q(xk(i)|xk-1,s(i),Z_1:k) of Equation (10) by Equations (33) and (34), and draw particles {xk|k(i),wk|k-1(i)}i=1L_k-1 for survival targets by Equations (3) and (4), where mk-1,s(i)=xk-1,s(i). 2 Approximate the IS function q(xk,b(i)|xk-1,b(i),Z_1:k) of Equation (10) using Equations (33) and (34), and draw particle set {xk|k(i),wk|k-1(i)}i=L_k-1+1L_k-1+J_k-1 for birth particles by Equations (3) and (4), where mk-1,b(i)=xk-1,b(i). 3 Calculate wk|k(i) by Equation (6) for resampling and wk|k(i,j) by Equation (42) for estimating, using the particle set {xk|k(i),wk|k-1(i)}i=1L_k-1+J_k-1. 4 Compute W^j by Equation (43) with wk|k(i,j) 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 Xke by Equation (44) with wk|k(i,j) and xk|k,j(i), where wk|k(i,j) and xk|k,j(i) belong to the j-th group of step 4. In addition, target number can be approximated by M_k=∣Xke∣. 6 Resample particles {xk|k(i),Pk|k(i),wk|k(i)}i=1L_k-1+J_k-1 to obtain {xk,s(i),Pk,s(i),wk,s(i)}i=1L_k for the next recursion, where L_k=[∑i=1L_k-1+J_kwk|k(i)], 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 {xk,b(i),Pk,b(i),wk,b(i)}i=1J_k when B_k is given. 9 Return M_k and Xke.

4. 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.

4.1. 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 is the number of birth components. and denote the mean and covariance of the i-th birth component, respectively. The posterior intensity of Equation (1) can be approximated by where and are the mean and covariance of the i-th survival component, respectively. is its corresponding weight, is the number of survival components, and is the weight of the i-th birth component. Here, .

Commonly, birth components at time become survival components at time k. We combine birth and survival components into one set. That is , where 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 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 where . In addition, and 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 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 of Equation (43) can be used to assign the particles into the j-th group, where the mean and covariance can be calculated by Equations (44) and (45) in Section 3.3. Given , and , we model the j-th group as a Gaussian mixture component , where we set . With such a selection, it can guarantee that the groups with large have enough number of sampling particles. Note that the group with small 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 for , subject to . In this paper, we set . The construction procedure of the target components is summarized in Table 3.

Table 3

The construction of target components.

– Input: The group weight W^j, particle {xk,j(i),wki,j}, ,j=0,1,…,|Z^_k|, and i=1,2,…,L_k-1+J_k – Output: The target components {Gk(a),m_k,a,P_k,p}a=1N_s,k. 1 Initiate the weight of the target component Gk(1)=0,N_s,k=0. 2 Normalize W^j by W˜^j=W^jmax_jW^j 3 for: j=0,1,…,|Z_k| if W˜^j>T_g – Add the number of target components N_s,k=N_s,k+1. – Substitute {wk(i,j)} and {xk(i)} into Equation (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 Gk(N_s,k)=W˜^j. 4 Return the target components {Gk(a),m_k,a,P_k,p}a=1N_s,k.

The target components of Table 3 refer to the survival target components. The birth components and target state estimation can be achieved by Section 3.2 and Section 3.3, respectively.

4.2. 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 .

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 , where is the dimension of the particle state. The computational complexity of the CIF-based IS functions is . Besides, the computational complexity of the PHD update step is , where M is the number of observations. In addition, the resampling step of our CISMC-PHD approach consumes computational complexity. Thus, the computational complexity of our CISMC-PHD approach in total is .

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 , where is the number of target components, . Assuming that in the PHD prediction step, target components generate particles. The computational complexity of the PHD update step is . Besides, the Gaussian target component forming consumes . Combining the computational complexity of these steps together, the computational complexity of the F-CISMC-PHD approach is . In practice, , 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 , where and are the computational complexity of the update and resampling steps. Obviously, . 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.

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.

5.1. Simulation Scenarios

In our simulations, we use the nonlinear scenario, the same as [36]. Let be the target states, represented by . In this paper, is the position, is the velocity, and α is the turn rate. With the above definitions, we model the nonlinear dynamic equation as where second (s) is the sampling interval. In addition, is the noise, defined by . We denote as covariance matrices of . In this paper, we set meter/second (m/s), rad. The initial states, appearing times, and disappearing times of targets are listed in Table 4.

Table 4

Initial states of the targets.

Target	State	Appearing(s)	Disappearing(s)
1	[320,5,320,5,0]	1	40
2	[400,-5,400,5,0]	8	50
3	[375,5,375,-5,0]	25	70
4	[400,5,325,-5,0]	59	70
5	[325,-5,375,5,0]	59	70

Besides, the observation model is given by where is the observation noise defined by , and is the covariance. Here, we set rad and meter (m). We also assume that clutters are uniformly distributed in the detection region, where the angle range is rad, and distance range is 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.

Figure 2

Ground-truth trajectories of five targets with the clutter number setting to be 10. The target trajectories are depicted by circle-solid lines with different colors, while the asterisks denote clutters.

For parameters, we set the gating threshold according to [32], the probability of detection and survival are , and , 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.

5.2. 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 and be two RFSs, where m and n are numbers of elements in and , respectively. Supposing that represents the set of permutations of , the first order OSPA metric can be rewritten by where , is a cut-off factor, and is the distance between and . In this paper, we set in accordance with [43], and use the Euclidean distance to compute .

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 Figure 3b,c, while most of points are not covered with the ground-truth trajectories.

Figure 3

Estimated trajectories of three approaches with clutter number being 10. The estimated trajectories are represented with the red (light) points, while the true trajectories are with the black (dark) solid line. (a) Trajectories of IBSMC-PHD; (b) Trajectories of CISMC-PHD; (c) Trajectories of F-CISMC-PHD.

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.

Figure 4

OSPA distances of the IBSMC-PHD, CISMC-PHD and F-CISMC-PHD approaches with clutter number being 10.

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.

Figure 5

Estimated numbers and RMSEs of IBSMC-PHD, CISMC-PHD and F-CISMC-PHD approaches with clutter number being 10. (a) Estimated numbers of the three approaches; (b) RMSEs of the three approaches.

Table 5

Averaged Estimation Errors and Computational Times per 100 particles.

Approaches	OSPA (m)	RMSE	Time (s)
IBSMC-PHD	46.20	0.44	0.02
CISMC-PHD	25.48	0.24	0.2
F-CISMC-PHD	31.93	0.30	0.06

5.3. 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.

Figure 6

Estimated error of the IBSMC-PHD, CISMC-PHD and F-CISMC-PHD approaches along with the clutter number changing from 1 to 30. (a) Averaged OSPA distances; (b) Averaged RMSEs of estimated number.

5.4. 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.

	p_d=0.92			p_d=0.94
	IBSMC-PHD	CISMC-PHD	F-CISMC-PHD	IBSMC-PHD	CISMC-PHD	F-CISMC-PHD
OSPA(m)	52.13	34.40	40.21	43.55	27.39	29.82
RMSE	0.54	0.37	0.47	0.47	0.32	0.37
	p_d=0.96			p_d=0.98
	IBSMC-PHD	CISMC-PHD	F-CISMC-PHD	IBSMC-PHD	CISMC-PHD	F-CISMC-PHD
OSPA(m)	43.55	27.39	29.82	39.35	23.20	23.02
RMSE	0.42	0.27	0.31	0.35	0.22	0.24

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 to . Thus, the F-CISMC-PHD approach is suitable for the high probabilities of detection.

5.5. 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 to . 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 to . 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.

Table 7

Tracking Performance over different .

	σ_θ=1.5π180			σ_θ=2π180
	IBSMC-PHD	CISMC-PHD	F-CISMC-PHD	IBSMC-PHD	CISMC-PHD	F-CISMC-PHD
OSPA(m)	49.62	27.32	29.66	50.40	27.64	29.99
RMSE	0.43	0.33	0.37	0.47	0.34	0.38
	σ_θ=2.5π180			σ_θ=3π180
	IBSMC-PHD	CISMC-PHD	F-CISMC-PHD	IBSMC-PHD	CISMC-PHD	F-CISMC-PHD
OSPA(m)	51.04	27.43	30.74	52.56	28.02	31.3
RMSE	0.50	0.35	0.39	0.45	0.35	0.39

6. 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.