A Student's t Mixture Probability Hypothesis Density Filter for Multi-Target Tracking with Outliers.

In multi-target tracking, the outliers-corrupted process and measurement noises can reduce the performance of the probability hypothesis density (PHD) filter severely. To solve the problem, this paper proposed a novel PHD filter, called Student's t mixture PHD (STM-PHD) filter. The proposed filter models the heavy-tailed process noise and measurement noise as a Student's t distribution as well as approximates the multi-target intensity as a mixture of Student's t components to be propagated in time. Then, a closed PHD recursion is obtained based on Student's t approximation. Our approach can make full use of the heavy-tailed characteristic of a Student's t distribution to handle the situations with heavy-tailed process and the measurement noises. The simulation results verify that the proposed filter can overcome the negative effect generated by outliers and maintain a good tracking accuracy in the simultaneous presence of process and measurement outliers.

1. Introduction

Multi-target tracking (MTT) plays an important role in many sensing systems, such as infrared, radar, sonar, etc., which uses the sensor data to jointly estimate the target state and the number of targets. Nowadays it is widely used in civilian and military applications such as air traffic control, remote sensing, ballistic missile guidance, and computer vision [1,2]. In MTT, the time-varying number of targets causes a problem in that the associations between state and measurement sets of targets are hard to know, which makes the traditional data association-based multi-target tracking methods problematic. In recent years, the random finite set (RFS) theory-based multi-target tracking filters, such as probability hypothesis density (PHD) filter [3], cardinalized PHD (CPHD) filter [4], multi-target multi-Bernoulli (MeMber) filter [1] and cardinality-balanced MeMBer (CBMeMBer) filter [5], have attracted much more attention since they can avoid the combinatorial problem that arises from data association. Moreover, some labeled RFS-based multi-Bernoulli filters [6,7] which can accommodate target tracks were proposed.

The focus of this paper is the PHD filter, which has relatively simple recursion, making it suitable for the applications demanding real time results. The PHD filter provides a tractable sub-optimal strategy for jointly estimating the number and the state of a variable number of targets by propagating the first-order statistical moment of multi-target posterior probability density in time. It has two basic implementations: the sequential Monte Carlo (SMC) method [8] and the Gaussian mixture (GM) method [9] which can solve the problem of computationally intractable multiple integrals involved in the PHD recursions. Compared to SMC implementation, GM implementation of the PHD filter has the advantages of simple state extraction and low computational cost, which is suitable for the requirement of real-time scenes. Moreover, some nonlinear extensions [9,10,11] and improvements [12,13,14] extend the scope of applications for the GM-PHD filter.

In MTT, noise, as the important part of measurement uncertainty, is an inevitable problem that reduces estimation accuracy of the PHD filter. Vo [15] thought that the setting of a reasonably large noise variance can accommodate noise interference in most situations. However, this method is only suitable for Gaussian noise. In real applications, it is hard for the measurement noise from sensor data to follow the Gaussian distribution because of electromagnetic interference or sensors’ own unreliability. Such measurements with outliers, which often express heavy-tailed character, degrade the performance of the PHD filter strikingly. What’s worse, in some real applications, such as tracking some agile targets with unreliable sensors, outliers may appear in not only the measurement model, but also the process model. This situation with simultaneous heavy-tailed process and measurement noises reduces the performance of the PHD filter severely, and can even make it break down. Although the SMC-PHD filter can deal with the problem to a certain degree, it has to pay a high computational cost, especially in high dimensions. For the GM-PHD filter, its foundational Gaussian approximation limits the capability to handle heavy-tailed non-Gaussian noises. Although Huber’s M-estimation theory [16] or variational Bayesian (VB) method [17] can be utilized to improve the performance of the GM-PHD filter against outliers in the measurement model, they both cannot handle the outliers in the process model. More importantly, the two methods above do not change the foundation of the Gaussian approximation-based GM-PHD filter. This means that the noise model still cannot match the outliers-corrupted process and measurement noises well, leading to biased estimates of the target state and the number of targets. Obviously, the Gaussian noise model cannot handle the heavy-tailed non-Gaussian noise, so how to model the heavy-tailed noise becomes the key point. As [18] said, Student’s t distribution, which has a heavy tail characteristic, is a good choice to match the heavy-tailed non-Gaussian noise. Under the Bayesian filtering framework, the Student’s t approximation-based closed form recursions are obtained for the linear system [18]. Further, the Student’s t approximation-based approach also can be used in the nonlinear system [19,20,21,22]. However, up to the present, Student’s t approximation-based approaches to approximate a PHD filter with heavy-tailed process and measurement noises do not exist.

In this paper, a novel implementation of the PHD filter is proposed based on Student’s t mixture approximation, intending to improve the estimation accuracy in terms of the target states and the target number in the presence of heavy-tailed process and measurement noises. The proposed approach models the process noise and the measurement noise as a Student’s t distribution, meanwhile, the multi-target prior intensity is approximated as a mixture of the Student’s t distributions. Then, the Student’s t mixture-approximated predicted intensity and posterior intensity are obtained through utilizing Student’s t approximation, forming a closed form recursion of the PHD filter. The Student’s t mixture implementation is proposed in RFS-based MTT algorithms for the first time. Compared to the GM case, it is a Student’s t-based implementation, which propagates a mixture of Student’s t components. Because it utilizes the heavy-tailed characteristic of Student’s t distribution, the proposed filter has better accuracy and robustness in MTT scenes with heavy-tailed process and measurement noises. Moreover, the proposed filter also has relatively low computing cost just like the GM-PHD filter. The above advantages of the proposed approach are verified by simulations designed in linear scenario and nonlinear scenario, respectively.

The remainder of this paper is organized as follows: Section 2 presents an overview of the PHD filter and some properties of the Student’s t distribution. Section 3 presents the proposed filter for linear system and extends the proposed filter to the nonlinear system. Simulation results are given in Section 4, and conclusions are drawn in Section 5.

2. Background

2.1. The PHD Filter

The random finite set (RFS) approach [1] provides a mathematically elegant treatment for the difficult problem that is how to estimate the time-varying number of targets and their states jointly in MTT scenarios. According to RFS theory, the collections of target states and measurements at time k can be represented as finite sets () and () where and are the number of targets and the number of measurements, respectively, and are the collections of all finite subsets of target states and measurements, respectively, and are the dimensions of the state space and the measurement space , respectively. Let: where is the RFS of survival targets at time k from target states and is the RFS of spawned targets at time k from target states , and is the RFS of birth targets at time k.

The RFS of measurements is defined as: where and are the RFS of measurements coming from target and clutter at time k, respectively.

Based on the RFS approach, the MTT problem was recast into a Bayesian filtering framework, but optimal Bayesian recursions are computationally intractable due to their sets integral. Then the PHD filter, which recursively propagates the first-order moment of multi-target posterior probability density, was proposed as an approximate solution in [3]. Let denote the posterior intensity function at time k − 1, then the prediction step of the PHD filter is given by where is transition density function in a Markov process from time k − 1 to time k, is the survival probability of each target at time k, is the intensity function of the RFS of targets spawning from previous state u, is the intensity function of birth targets at time k.

Given the predicted intensity , the update intensity can be given by: where is the probability of detection at time k, denotes the measurement likelihood function at time k, is the intensity function of clutter. The integration domain of the integral functions in (3) and (4) is the state space .

2.2. Review of the GM-PHD Filter

The GM-PHD filter approximates the intensity function of multi-target as Gaussian mixture components under the assumption of the linear Gaussian dynamical model and measurement model. In addition, it still needs to satisfy the following assumptions: (1) the survival and detection probabilities are state-independent; (2) the intensities of the birth and spawn RFSs are Gaussian mixtures of the form. Then a closed-form solution of the PHD filter can be obtained as follows.

Prediction: The Gaussian mixture formulation of the PHD recursion at time k − 1 is given as: where denotes that random vector x follows Gaussian distribution, and are the mean and the covariance matrix of the jth item in all Gaussian mixture components, respectively, is the corresponding weight, is the number of Gaussian components.

Assume that the intensities of survival targets, spawned targets and birth targets are Gaussian mixtures of the form, the predicted intensity is given by: where:

Update: Given the PHD predictor as Equation (6), the PHD updater is formulated as: where: in the above, and are the state transition matrix and observation matrix, respectively, and are the process and measurement noise covariance, respectively.

2.3. Student’s t Distribution

Let a random vector admit the Student’s t distribution; its probability density function (PDF) can then be expressed as [23]: where denotes Gamma function and . The above PDF abbreviated by with mean m, scale matrix P and degrees of freedom parameter . The corresponding covariance matrix can be calculated as () [23].

The tail behavior of the Student’s t distribution is very much influenced by the degrees of freedom parameter . The smaller is, the heavier the tail is, and vice versa. In addition, the Student’s t distribution reduces to the Gaussian as tends to infinity, thus includes it as a special case. A number of convenient properties are shared by both and can be easily derived. As for the PDF of affine mappings of Student’s t variables [23] we have that , with appropriate A and b, admits:

The degrees of freedom parameter remains unaltered. Turning to random vectors and that are jointly Student’s t distributed with:

The marginal PDF of can be computed by applying a linear transformation with to (18):

From (18) and (19), the conditional PDF is also a Student’s t as follows: with:

It can be seen that the conditional mean (21) corresponds to the Gaussian conditional mean. The matrix parameter (22) is a scaled version of the conditional covariance in the Gaussian case, which is recovered as tends to infinity. In contrast to the Gaussian, depends on . Also, the degrees of freedom parameter in (23) increases. Above properties of Student’s t distribution like (17) and (20) are the foundations of Lemmas 1 and 2 (in Section 3.3) that are the key of our proposed approach.

3. Student’s t Mixture PHD Recursion

Differing from the GM-PHD recursion, the proposed approach approximates the intensity of multi-target RFS as the Student’s t mixture of the form, which is propagated in the closed-form recursions. The following gives the novel algorithm for linear system at first and then extends it to nonlinear system.

3.1. Basic Assumptions for Linear Model

As for the linear system, some foundational assumptions are given as follows.

Given that the process and measurement noises admit the zero-mean t distribution with scale matrix where

The survival and detection probabilities are state independent, i.e.:

The intensities of the birth and spawn RFSs are Student’s t mixtures of the form:

Assumption 1 is given according the affine transition nature of Student’s t variables [19]. Assumption 2 is commonly used in GM-PHD filters. Assumption 3 is given in the presence of process and measurement noises with heavy tails, so it is reasonable in actual applications, i.e., when using unreliable sensors or sensors suffer from electromagnetic interference while tracking some agile targets, where both process noise and measurement noise are prone to show a heavy-tailed character. For illustrating this, a one-dimensional heavy-tailed noise distribution and Gaussian noise distribution are given as Figure 1.

Figure 1

Illustration of heavy-tailed noise distribution and Gaussian noise distribution.

3.2. Student’s t Mixture PHD Recursion

Based on the Assumptions 1–3, the closed solution to the PHD recursion (3) and (4) is presented as two Propositions. The proposed Propositions show how the Student’s t components of the posterior intensity are analytically propagated to the next time, analogous to the GM-PHD recursion.

Suppose that Assumptions 1–3 hold and that the posterior intensity at time k − 1 is a Student’s t mixture of the form:

Then, the predicted intensity at time k is also a Student’s t mixture and is given by: where:

Suppose that Assumptions 1–3 hold and that the predicted intensity for time k is a Student’s t mixture of the form:

Then, the posterior intensity at time k is also a Student’s t mixture and is given by: where:

Aforementioned Propositions 1 and 2 can be established by applying the following approximated results for Student’s t functions.

Given that the jointly PDF of the current state and one-step ahead state vectors is Student’s t and F, Q, m and P of appropriate dimensions and that Q and P are positive definite:

Given that the jointly PDF of the state and measurement vectors is Student’s t and H, R, m, P of appropriate dimensions and that R and P are positive definite: where:

Lemmas 1 and 2 are derived based on the properties of Student’s t distribution listed in Section 2.3. The detailed derivations of Lemmas 1 and 2 can be seen [19].

Proposition 1 is established by substituting (25), (27) and (29)–(31) into the PHD prediction (3), and replacing integrals of the form (51) by appropriate Student’s t as given by Lemma 1. Similarly, Proposition 2 is established by substituting (26), (28) and (39) into the PHD update (4), and then replacing integrals of the form (51) and product of Student’s t of the form (52) by appropriate Student’s t as given by Lemmas 1 and 2 respectively. The concrete proofs of Propositions 1 and 2 are given in Appendix B.

From (45), we know that the degree of freedom will increase infinitely with recursion performing, which results in that the Student’s t mixture degrades the Gaussian mixture according to [18]. It means that the robustness against outliers for the Student’s t mixture PHD filter will be lost with time going by. To solve the problem, we adapts the moment matching approach [18] to obtain the correction of posterior intensity at time k given as (59)–(61). A pseudocode of main process of the proposed algorithm is given by Table A1 in Appendix A:

Table A1

Pseudocode for the STM-PHD filter.

given {wk−1(i),mk−1(i),Pk−1(i)}i=1J_k−1, the measurement set Z_k and the degree of freedom parameters υ_1,υ_2,υ_3 .step1.(prediction for birth targets)        i=0.        for j=1,…,J_γ,k        i:=i+1.        wk|k−1(i)=wγ,k(j),   mk|k−1(i)=mγ,k(j),   Pk|k−1(i)=Pγ,k(j),        end        for j=1,…,J_β,k               for l=1,…,J_k−1                    i:=i+1.                    wk|k−1(i)=wk−1(l)wβ,k(j),              mk|k−1(i)=Fβ,k−1(j)mk−1(l)+dβ,k−1(j),                    Pk|k−1(i)=Fβ,k−1(j)Pk−1(l)(Fβ,k−1(j))^T+υ_1(υ_3−2)υ_3(υ_1−2)Qβ,k−1,                end         endstep2. (prediction for existing targets)        for j=1,…,J_k−1               i:=i+1.               wk|k−1(i)=p_S,kwk−1(j),               mk|k−1(i)=Fk−1(j)mk−1(j),   Pk|k−1(i)=Fk−1Pk−1(j)Fk−1T+υ_1(υ_3−2)υ_3(υ_1−2)Qk−1,        end        J_k|k−1=i.step3.(construction of PHD update components)        for j=1,…,J_k|k−1               ηk|k−1(j)=H_kmk|k−1(j),             Pzz,k|k−1(j)=HkPk|k−1(j)HkT+υ_2(υ_3−2)υ_3(υ_2−2)Rk,               Kk(j)=Pk|k−1(j)HkT[Pzz,k|k−1(j)]^−1,     Pk|k(j)=Pk|k−1(j)−Kk(j)HkPk|k−1(j),        endstep4.(update)        for j=1,…,J_k|k−1               wk(j)=(1−p_D,k)wk|k−1(j),               mk(j)=mk|k−1(j),     Pk(j)=Pk|k−1(j),        end        l=0.        for each z∈Z_k               l:=l+1.               for j=1,…,J_k|k−1                    Δ^(lJ_k|k−1+j)=(z−ηk|k−1(j))^T(Pzz,k|k−1(j))^−1(z−ηk|k−1(j))                    wk(lJ_k|k−1+j)=p_D,kwk|k−1(j)St(z;ηk|k−1(j),Sk(j),υ_3),                    m^′k(lJ_k|k−1+j)=mk|k−1(j)+Kk(j)(z−ηk|k−1(j)),                    P^′k(lJ_k|k−1+j)=υ_3+(Δ^(lJ_k|k−1+j))^2υ_3+d_zPk|k(j),               end               wk(lJ_k|k−1+j):=wk(lJ_k|k−1+j)κ_k(z)+∑i=1J_k|k−1wk(lJ_k|k−1+i),   for j=1,…,J_k|k−1.        end        υ^′_3=υ_3+d_z,      J_k=lJ_k|k−1+J_k|k−1.step5.(moment matching)        for j=1,…,J_k               mk(j)=m^′k(j),       Pk(j)=υ^′_3(υ_3−2)υ_3(υ^′_3−2)P^′k(j),        endoutput{wk(i),mk(i),Pk(i)}i=1J_k

The degree of freedom parameters for the process noise model, measurement noise model and the initial multi-target state intensity are different in general. As in the recursion performed from (31) to (39), there is a problem that how to select the degree of freedom parameter between and as the degree of freedom parameter of predicted multi-target intensity to be propagated in the next recursive step. A valid method [18] is that choose the minimized value between and to be propagated. The same problem, existing in the recursion from (39) to (40), also can be handled by selecting the minimized value between and . For simplicity, this paper assumes that the degree of freedom parameters are equal.

Analogous to the derivation of GM implementation of the CPHD filter in [24], the proposed Student’s t mixture implementation shown as in Propositions 1 and 2 can also be used in the CPHD filter, forming the corresponding Student’s t mixture CPHD recursion.

3.3. Implementation Issues

Like the GM-PHD filter, the Student’s t mixture PHD (STM-PHD) filter also suffers from the computation problem that the number of Student’s t components increases endlessly with recursive time, so a pruning procedure and a merging procedure are necessary for the STM-PHD filter. The concrete procedures are similar to the procedures in the GM-PHD filter (readers can refer to [9]). The different point is that the Student’s t components give the scale matrix not the covariance matrix, so the covariance matrix, calculated by , should be used in the merging procedure of the STM-PHD recursion. The estimated number of targets is obtained by summing up the weights of all the Student’s t components. This procedure for the STM-PHD filter is no different from the GM-PHD filter. Again, the state extraction is also analogous to the GM-PHD filter as to select the means of the Student’s t components that have weights greater than some threshold (generally set as 0.5 [9]).

In addition, for the scene with high clutter density, gating strategy is always used to reduce the computing cost for the GM-PHD filter. It is easy to know that the gating strategy is also suitable for the STM-PHD filter. The core principle to select the measurement can be expressed as: where denotes the reduced set of the measurement at time k, is the innovation covariance matrix corresponding to the jth predicted measurement, and is the gate threshold. In contrast to the Gaussian mixture case, the judging variable for the Student’s t mixture case follows an F-distribution not a chi-squared distribution. It means that the gate threshold for the STM-PHD filter, which can be chosen from F-distribution table, is different from the GM-PHD filter at the same probability regions.

3.4. Extension to Nonlinear Model

This section considers the situation that process and measurement models are nonlinear. The models in Assumption 1 change to become: where fk and h are known nonlinear functions, w−1 and v are additional noises, which follow zero-mean Student’s t distribution with scale matrix Q−1 and R, respectively.

Different from the GM-PHD filter to cope with nonlinear problems via giving the numerical solution to Gaussian integrals, the core problem of the STM-PHD filter for nonlinear systems is how to compute the Student’s t integrals. Some researchers have given numerical solutions to Student’s t integrals based on Taylor linearization, unscented transform or cubature rule [19,21,22]. Compared with Taylor linearization and unscented transform, the cubature rule-based filter is a derivative-free and vigorous method [25], so this paper utilizes the cubature rule to extend the STM-PHD filter to nonlinear system according to [22]. The proposed algorithm for handle nonlinearity is given as Table A2 in Appendix A.

Table A2

Pseudocode for the cubature rule based STM-PHD filter.

given {wk−1(i),mk−1(i),Pk−1(i)}i=1J_k−1, the measurement set Z_k and the degree of freedom parameters υ_1,υ_2,υ_3 .step1.(prediction for birth targets)        follow Step 1. of table.        for j=1,…,i

- use the cubature rule suitable for Student′s t integral (see [22]) with mk−1(j) and Pk−1(j)

to obtain a set of cubature points and weights, denoted by {xk−1(l),μ(l)}l=1L                    - compute                    xk|k−1(l):=f_β,k(xk−1(l)),        l=1,…,L,                    mk|k−1(j)=∑l=1Lμ(l)xk|k−1(l),                    Pk|k−1(j)=υ_3−2υ_3∑l=1Lμk−1(l)(xk|k−1(l)−mk|k−1(j))(xk|k−1(l)−mk|k−1(j))^T+υ_1(υ_3−2)υ_3(υ_1−2)Qβ,k−1,        endstep2. (prediction for existing targets)        for j=1,…,J_k−1               i:=i+1.               wk|k−1(i)=p_S,kwk−1(j),               - use the cubature rule suitable for Student′s t integral with mk−1(j) and Pk−1(j)                to obtain a set of cubature points and weights, denoted by {xk−1(l),μ(l)}l=1L               - compute               xk|k−1(l):=f_k(xk−1(l)),        l=1,…,L,               mk|k−1(j)=∑l=1Lμk−1(l)xk|k−1(l),               Pk|k−1(j)=υ_3−2υ_3∑l=1Lμ^(l)(xk|k−1(l)−mk|k−1(j))(xk|k−1(l)−mk|k−1(j))^T+υ_1(υ_3−2)υ_3(υ_1−2)Qk−1,        end        J_k|k−1=i.step3.(construction of PHD update components)        for j=1,…,J_k|k−1               - compute               zk|k−1(l):=h_k(xk|k−1(l)),        l=1,…,L,               ηk|k−1(j)=∑l=1Lμ(l)zk|k−1(l),               Pzz,k|k−1(j)=υ_3−2υ_3∑l=1Lμ^(l)(zk|k−1(l)−ηk|k−1(j))(zk|k−1(l)−ηk|k−1(j))^T+υ_2(υ_3−2)υ_3(υ_2−2)Rk,               Pxz,k|k−1(j)=υ_3−2υ_3∑l=1Lμ^(l)(xk|k−1(l)−mk|k−1(j))(zk|k−1(l)−ηk|k−1(j))^T               Kk(j)=Pxz,k|k−1(j)[Pzz,k|k−1(j)]^−1,               Pk|k(j)=Pk|k−1(j)−Pxz,k|k−1(j)[Pzz,k|k−1(j)]^−1[Pxz,k|k−1(j)]^T,        endstep4.(update) andstep5.(moment matching)        follow Step 4. and Step 5. of table to obtain {wk(i),mk(i),Pk(i)}i=1J_k.output {wk(i),mk(i),Pk(i)}i=1J_k.

4. Simulations and Results

To illustrate the performance of the proposed filter, simulation examples are designed to compare with standard GM-PHD filter in linear and nonlinear scenarios, respectively. To compare the performance of two filters, we choose the Optimal Sub-pattern Assignment (OSPA) distance as the metric, which can comprehensively measure the cardinality and localization errors [26]. The OSPA distance is defined as follows. Let for , and denotes the set of permutations on for any . For , and arbitrary finite subsets and belong to , where :

If , and if ; and if . is the order parameter that determines the sensitivity to outliers and is the cut-off parameter that determines the relative weighting of the penalties assigned to cardinality and localization errors. The details to choose the parameters and can be seen in [26]. In our simulation examples, we set and .

4.1. Linear Scenario

Consider a two-dimensional scenario where there are twelve targets over region × during the interval of 100 s. Assuming no target spawning and each target moves as a constant velocity model similar to [9] with: where and . The state of each target consist of position and velocity at time k. Their corresponding initial state and life time of each target are given as Table 1.

Table 1

A list of initial target states.

Target Index	Life Time (s)	Initial States (m, m/s, m, m/s)
#1	(1, 70)	[0, 0, 0, −10]
#2	(1, 100)	[400, −10, −600, 5]
#3	(1, 70)	[−800, 20, −200, −5]
#4	(20, 100)	[400, −7, −600, −4]
#5	(20, 100)	[400, −2.5, −600, 10]
#6	(20, 100)	[0, 7.5, 0, −5]
#7	(40, 100)	[−800, 12, −200, 7]
#8	(40, 100)	[−200, −3, 800, −10]
#9	(60, 100)	[−800, 3, −200, 15]
#10	(60, 100)	[−200, −3, 800, −15]
#11	(80, 100)	[0, −20, 0, −15]
#12	(80, 100)	[−200, 15, 800, −5]

The noisy measurement model is the same as [9] with:

The process and measurement noises with heavy tails are given as (66) and (67):

For the process and measurement noises in (66) and (67), about one percent of process and measurement noise values are drawn from Gaussian with severely high covariance. This percentage is also called contaminated rate which can be denoted by [27].

Assuming no spawned target and birth targets appear spontaneously according to a Poisson point process with intensity function: where , , and , and . The true trajectories of each target are shown in Figure 2, while Figure 3 plots these trajectories with Gaussian measurements and heavy-tailed measurements over time (not plot clutter in figure). From Figure 3, it is can be seen that the individual heavy-tailed measurements obviously bias the true position compared with the corresponding Gaussian measurements, which may degrade the estimation accuracy.

Figure 2

True trajectories of each target.

Figure 3

Measurements and true target positions versus time: (a) in x coordinate; (b) in y coordinate.

The detection probability and target survival probability are and , respectively. Truncated threshold, merged threshold and the maximum Student’s t components related to pruning and merging process are , and , respectively. For simplification, set .

To evaluate the performance of the STM-PHD filter, we compare it with GM-PHD filter over 100 Monte Carlo (MC) trails with fixed clutter density. Under the uniform distribution assumption, the clutter density can be given by clutter rate λc with the relationship = λc/V. In this simulation, we set λc = 20 (giving an average of 20 clutter returns per scan). Figure 4 and Figure 5 respectively show the estimated cardinality and the OSPA distance for two filters. The result in Figure 4 shows that the STM-PHD filter provides a noticeable improvement in terms of cardinality estimation accuracy compared with the GM-PHD filter, although some biased cardinality estimates appear for the STM-PHD filter.

Figure 4

Comparison of cardinality estimation of two filters with fixed clutter rate (λc = 20).

Figure 5

Comparison of OSPA distance of two filters with fixed clutter rate (λc = 20).

In Figure 5 the OSPA distance of the STM-PHD filter is lower than that of the GM-PHD filter. Especially after 40 s more targets appear, difference of the OSPA distance between two filters is more noticeable. The main reason is that the STM-PHD filter has more accurate cardinality estimation.

To evaluate the performance of the proposed filter sufficiently, a simulation is executed over 100 MC trials with different contamination rates from to . Then the time averaged OSPA distance of the STM-PHD filter and the GM-PHD filter, respectively, are shown in Figure 6.

Figure 6

Comparison of OSPA distance of two filters with different contaminated rate.

From Figure 6, it can be seen that the time averaged OSPA distance of the STM-PHD filter is lower than that of the GM-PHD filter overall. The time averaged OSPA distances of two filters increases with the increasing contamination rate. Remarkably, the gap of OSPA distance between the STM-PHD filter and the GM-PHD filter changes wider from to . It means that the STM-PHD filter has strong robustness against the negative effect of outliers, especially for high contamination rates. This is due to the fact the Student’s t noise model in the proposed approach can match the heavy-tailed non-Gaussian noise well. On the contrary, the Gaussian-based GM approach matches such a non-Gaussian noise worse and worse with the increasing of contaminated rate. Additionally, at , the OSPA distance for the STM-PHD filter is the same as the GM-PHD filter. It indicates that the STM-PHD filter and the GM-PHD filter have the same tracking performance when outliers do not exist.

To further evaluate the performance of the proposed filter, a simulation is performed over 100 MC trails with different clutter rates from λc = 0 to λc = 50. The resulting time averaged OSPA distances of the proposed filter and the GM-PHD filter are shown in Figure 7. It can be seen that the time averaged OSPA distances of two filters increase with the increasing clutter rate and the time averaged OSPA distance of the STM-PHD filter is always lower than that of the GM-PHD filter under different clutter rates. This means that the STM-PHD filter generally outperforms the GM-PHD filter when outliers exist, no matter what the clutter rate is.

Figure 7

Comparison of OSPA distance for two filters with different clutter rate.

In addition, the computational cost for the STM-PHD filter lies at the same level as that of the GM-PHD filter for the linear system. Running on a computer with an Intel(R) Core(TM) i5-4570 CPU at 3.2 GHz, the average computing times per execution of the GM-PHD filter and the STM-PHD filter with different clutter rate are given in Table 2.

Table 2

Average computing time with different clutter rates.

	Clutter Rate
	0	10	20	30	40	50
GM-PHD	0.9217 s	0.9917 s	1.0782 s	1.1253 s	1.1938 s	1.2068 s
STM-PHD	0.9146 s	0.9961 s	1.0780 s	1.1417 s	1.2576 s	1.2486 s

4.2. Nonlinear Scenario

In this example, we assume a maximum of ten targets appears on the observation region × and a nearly constant turn state model and nonlinear bearings and range measurement model are considered according to [9]. The state consists of position and velocity as well as the turn rate . The state model is given by: where:

The noisy measurement model with range and bearing measurement is given by:

Like the linear scenario, the outliers contaminated process and measurement noises can be given by: with , .

In the simulation, we assume no spawned target and that the birth target is Poisson with intensity: where , , and , and . (The unit of distance, angle and time in this paper are meter, radian and second, respectively.)

The initial target states are given by Table 3 and the true trajectories of each target are shown as Figure 8. In addition, Figure 9 plots corresponding measurements with Gaussian noise and heavy-tailed noise respectively over time (not plot clutter in figure). In Figure 9, it also shows the results analogous to the linear case as shown in Figure 3.

Table 3

A list of initial target states.

Target Index	Life Time (s)	Initial States (m, m/s, m, m/s, rad/s)
#1	(1, 100)	[1000, −10, 1500, −10, 2π/(180 × 8)]
#2	(10, 100)	[−250, 20, 1000, 3, −2π/(180 × 3)]
#3	(10, 100)	[−1500, 11, 250, 10, −2π/(180 × 2)]
#4	(10, 66)	[−1500, 43, 250, 0, 0]
#5	(20, 80)	[250, 11, 750, 5, 2π/(180 × 4)]
#6	(40, 100)	[−250, −12, 1000, −12, 2π/(180 × 2)]
#7	(40, 100)	[1000, 0, 1500, −10, 2π/(180 × 4)]
#8	(40, 80)	[250, −50, 750, 0, −2π/(180 × 4)]
#9	(60, 100)	[1000, −50, 1500, 00, −2π/180 × 4]
#10	(60, 100)	[250, −40, 750, 25, 2π/(180 × 4)]

Figure 8

True trajectories of each target.

Figure 9

Measurements and true target positions versus time: (a) in x coordinate; (b) in y coordinate.

To evaluate the performance of the CKF based STM-PHD filter to cope with nonlinear problem, we compare it with CKF based GM-PHD filter [11] over 100 Monte Carlo (MC) trails with fixed clutter rate λc = 20. Figure 10 and Figure 11 respectively show the estimated cardinality and the OSPA distance of two filters. From Figure 10, it can be seen that the STM-PHD filter is superior to the GM-PHD filter in terms of cardinality estimation accuracy, although the STM-PHD filter has cardinality bias when the number of targets increases. The main reason for generating cardinality bias is that some necessary approximations for coping with nonlinear problems induce errors.

Figure 10

Comparison of cardinality estimation of two filters with fixed clutter rate (λc = 20).

Figure 11

Comparison of OSPA distance of two filters with fixed clutter rate (λc = 20).

In Figure 11 the OSPA distances of the STM-PHD filter and the GM-PHD filter are at the same level before 65 s and later on the OSPA distance of the STM-PHD filter is obviously lower than that of the GM-PHD filter. This result matches the result in Figure 10 well. It indicates that the difference of OSPA distance between two filters is due to the difference of cardinality estimation accuracy.

To evaluate the performance of the STM-PHD filter sufficiently, simulation is performed over 100 MC trials with different contaminated rate from to . Then the time averaged OSPA distances of the STM-PHD filter and the GM-PHD filter are shown in Figure 12.

Figure 12

Comparison of OSPA distance of two filters with different contamination rates.

Analogous to the linear scenario, it can be seen that the time averaged OSPA distance of the STM-PHD filter is lower than that of the GM-PHD filter at almost all contaminated rates. The gap of the time averaged OSPA distance between the STM-PHD filter and the GM-PHD filter also changes widely as the contamination rates increase. Nevertheless, the gap is not noticeable like in the linear scenario. This is due to the relatively big approximation error induced in the nonlinear scenario. The result indicates that the STM-PHD filter still outperforms the GM-PHD filter to cope with outliers for nonlinear systems.

To further evaluate the performance of the proposed filter, 100 MC trails are performed from λc = 0 to λc = 50. The time averaged OSPA distances versus varying clutter rate for the STM-PHD filter and the GM-PHD filter are shown in Figure 13. It can be seen that the time averaged OSPA distance of the STM-PHD filter is lower than that of the GM-PHD filter at different clutter rates, and the trend for two filters goes up with the increase of clutter rate. Generally speaking, the STM-PHD filter is superior to the GM-PHD filter but the superiority is not noticeable like in the linear scenario. This is due to the fact the approximation error induced in a nonlinear scenario is bigger than that in a linear scenario. Nevertheless, the results still indicate that the STM-PHD filter is valid to handle the outliers.

Figure 13

Comparison of OSPA distance of two filters with different clutter rates.

Again, the computing time for each filter under the different clutter rate is given in Table 4. It shows that the STM-PHD filter has a higher computing cost compared to the GM-PHD filter and the higher the clutter rate is, much more time is consumed for the STM-PHD filter. The reason is that computing nonlinear Student’s t integrals is more complex than computing nonlinear Gaussian integrals.

Table 4

Average computing time with different clutter rate.

	Clutter Rate
	0	10	20	30	40	50
GM-PHD	1.0142 s	1.3780 s	1.8481 s	2.4988 s	3.3321 s	3.9051 s
STM-PHD	1.8534 s	3.3179 s	5.4081 s	8.5213 s	12.9163 s	16.6113 s

5. Conclusions

To solve the problem that the process and measurement outliers degrade the performance of PHD filter, this paper proposes a Student’s t mixture-based PHD filter, where the prior multi-target intensities are approximated as a mixture of the Student’s t components to be propagated in time and the Student’s t mixture approximated posterior multi-target intensity are obtained through a Student’s t approximation-based recursion. The proposed approach can efficiently suppress the negative impact caused by the process and measurement outliers, while still maintaining good estimation accuracy. The advantages of the proposed filter are verified sufficiently by simulation. The simulation results also imply that the proposed approach outperforms the GM-based approach in scenarios where outliers appear due to electromagnetic interference or sensors’ own unreliability. In addition, the proposed approach is also suitable for the CPHD filter. However, the proposed approach has relatively high computational cost for nonlinear systems. In further study, we will try to improve this. Moreover, how to combine the proposed Student’s t mixture implementation with the multi-Bernoulli filters (including CBMeMBer filter and labeled RFS based multi-Bernoulli filter) is also our focus in the next work.