An enhanced particle filtering method for GMTI radar tracking

This paper investigates the problem of ground vehicle tracking with ground-moving target indicator (GMTI) radar. In practice, the movement of ground vehicles may involve several different maneuvering types (acceleration, deceleration, standstill, etc.). Consequently, the GMTI radar may lose measurements when the radial velocity of the ground vehicle is below a threshold, i.e., falling into the Doppler blind region. In this paper, to incorporate the information gathered from normal measurements and knowledge on the Doppler blindness constraint, we develop an enhanced particle filtering method for which the importance distributions are inspired by a recent noise-related Doppler blind (NRDB) filtering algorithm for GMTI tracking. Specifically, when constructing the importance distributions, the proposed particle filter takes the advantages of the efficient NRDB algorithm by applying the extended Kalman filter and its generalization for interval-censored measurements. In addition, the linearization and Gaussian approximations in the NRDB algorithm are corrected by the weighting process of the developed filtering method to achieve a more accurate GMTI tracking performance. The simulation results show that the proposed method substantially outperforms the existing methods for the GMTI tracking problem.


I. INTRODUCTION
In this paper, we consider the problem of ground vehicle tracking using a ground-moving target indicator (GMTI) radar that discriminates a moving target against the static background based on the Doppler effect [1]. GMTI radar is well suited for detecting targets moving on ground due to its wide-area, all-weather, day/night, and real-time capabilities [2]. Therefore, GMTI-based tracking has received a wide range of applications in continuously tracking of vehicles to support surveillance in different environments (e.g., battlefield and urban).
Due to its practical importance, there have been a number of methods developed in recent years to address various research issues associated with GMTI tracking. Kirubarajan et al. [3] proposed a variable structure interacting multiple model (VS-IMM) algorithm for GMTI tracking, considering there may be more than one state model to describe the target object's movement and the number of state models may change as well. In order to overcome the nonlinearity in the measurement model of the GMTI radar and the non-Gaussian posterior distribution of the state vector, a particle filtering approach was proposed in [4] and the simulation results showed an improved performance with reduced root mean square error (RMSE) than [3]. In addition, a new approach was proposed in [5] to improve on the performance of [4] in which the traditional particle filter was replaced with a more advanced unscented particle filter developed in [6].
The major limitation of the methods in [3][4][5] is that the GMTI measurements were assumed to be recorded at every time step in these studies, which is not the case in the real-world problems. In practice, no GMTI measurements are received if the magnitude of a target object's radial velocity drops below the minimum detectable velocity.
In order to address this GMTI tracking problem, the researchers in [7] incorporated a separate stop model into the VS-IMM framework as a complement to the existing maneuvering models. The incorporation of the stop model in the VS-IMM framework improves the tracking performance when no measurements are recorded due to vehicle stoppage. On the other hand, there are other studies such as [2,8,9] that applied Gaussian mixture tracking algorithms by propagating the Gaussian mixture approximation to the conditional distribution of the target state. In addition, a suitable state-dependent detection probability was introduced, which helps to determine the conditional distribution of the target state when no measurements are recorded. To solve the problem that negative weights may possibly arise in the Gaussian mixture approximation, an extra approximation stage was introduced in [10] to replace the resulting negative Gaussian mixture with one having strictly positive mixture weights, and thus improving on algorithmic stability. The algorithm in [11] applied the particle filtering method and treated each nondetection case as evidence. The corresponding likelihood function of the nondetection evidence was formulated and incorporated into the particle filtering procedure to update the target state probability distribution when no measurements are recorded.
Recently, Clark et al. [12] have developed a new Gaussian mixture model (GMM)-based noise-related Doppler blind (NRDB) filtering algorithm for GMTI tracking in which a GMM is used to approximate the posterior probability distribution of the state vector. This algorithm was carefully designed so that the Doppler blindness constraint of the GMTI radar can be efficiently dealt with under the Gaussian distribution assumption. The simulation results in [12] showed that this new approach outperformed some existing methods (such as [8]), especially when dealing with the scenario of no measurement. However, to ensure that the filter in [12] retains an analytically tractable form and to process signals efficiently, the following approximations have been made at each time step in this algorithm. 1) Approximation A: a standard mixture reduction technique (see, e.g., [13]) is used to avoid the exponential growth of the number of the Gaussian mixture components.
2) Approximation B: the posterior probability distribution of the state vector for each mode is approximated by a Gaussian distribution.
3) Approximation C: the measurement equation associated with the GMTI radar is linearized based on the first-order Taylor expansion, and a Gaussian distribution assumption is made when dealing with the Doppler blindness constraint.
While approximation A is the most commonly used approach to ensure the number of components is manageable in practice, approximations B and C are adopted in [12] to form a generalized extended Kalman filter (EKF) for the situation where a certain measurement is an interval-censored value, i.e., it falls into a known interval (the Doppler blind region). More specifically, when the radial velocity of the target is within the Doppler blind region, Clark et al. [12] have made use of the results in [14] to compute conditional mathematical expectations and covariance matrices so that knowledge on the Doppler blindness constraint can be efficiently utilized to update the state estimate.
The aim of this paper is to address the limitations of the algorithm in [12] and remove approximations A-C as much as possible. Because GMTI tracking is a nonlinear, non-Gaussian filtering problem, we will use a particle filtering approach, rather than the GMM-based algorithm with EKF-like filtering in [12], to estimate the state vector.
First, we will use the resampling method developed in [15] to replace the mixture reduction technique adopted in approximation A to ensure that the approximation errors are kept at a minimum level. In addition, we note that by applying the EKF and its generalized version, the GMM-based NRDB algorithm in [12] works very well in many scenarios. A distinguished feature of the method proposed in this paper is to fully take the advantages of the NRDB algorithm by treating the EKF and its generalized version as the importance distributions to generate particles in the proposed particle filter. The errors caused by the linearization and Gaussian approximations associated with the importance distributions (i.e., approximations B and C outlined above) are then corrected in the later stage of the particle filter via a suitable weighting process. This paper is divided into the following sections. Section II considers problem formulation of GMTI tracking and briefly introduces the Bayesian inference framework developed in [15]. The details of the proposed novel particle filtering algorithm for GMTI tracking are provided in Section III. A simulation study is given in Section IV that evaluates the numerical performance of the proposed method and compares its performance with other state-of-the-art methods. Finally, we close this paper with some concluding remarks in Section V.

A. State and Measurement Models
In a standard tracking problem, the system model includes a state equation of the target object and a measurement model for the sensor observations. The state equation is used to describe the dynamics of the movement of the target vehicle and the measurements provide the information for state filtering. The movement of the ground target vehicle may involve several different maneuvering types (acceleration, deceleration, standstill, etc.), and its dynamics are usually described using a Markov jump multiple model: where x t = [x t , y t ,ẋ t ,ẏ t ] is the state vector consisting of the positions (x t , y t ) and velocities (ẋ t ,ẏ t ) in x and y directions (here we assume that the vehicle moves on ground with z t ≡ 0 andż t ≡ 0). Note that more than one movement mode are normally involved in the tracking problem with m t ∈ M, where M is the set of mode indexes. The transitions among different modes are described by the transition probabilities. In this paper, we consider a Markov-jump system in [16] where the transition probability p(m t = i|m t-1 = j) from mode j to i is assumed to be constant. In addition, F(m t ) and w t-1 (m t ) represent the state transition matrix and process noise that could depend on a particular movement mode m t . The standard GMTI radar measures the range, azimuth angle, and range rate (denoted as y r , y θ , and yṙ at time step t, respectively) of the ground vehicle relative to the position of the GMTI radar. Following [12], we assume that these measurements are noise-corrupted from actual values: where arctan2 denotes the four quadrant inverse tangent function, (x o,t , y o,t , z o,t ) and (ẋ o,t ,ẏ o,t ,ż o,t ) represent the position and velocity of the observer (GMTI radar) at time t, respectively, h(x t ) = [r t , θ t ,ṙ t ] T is the ideal measurement vector without the measurement noises, and n t = [n r t , n θ t , nṙ t ] T represents the measurement noise vector that is assumed to be Gaussian with zero-mean vector and a diagonal covariance matrix diag{σ 2 r , σ 2 θ , σ 2 r }. To reduce the nonlinearity of (2), we adopt the unbiased measurement transformation method (see, e.g., [17,18] for details) to convert the nonlinear measurement components r t and θ t into the Cartesian measurements.
To deal with the Doppler blindness region, let κ denote the minimum detectable velocity of the GMTI radar and let v r t denote the radial velocity of the target that is defined to be the target's velocity projected along the range direction. According to the properties of GMTI radars, no measurements will be detected if the target radial velocity is within the Doppler blind region. In addition, even when the target radial velocity v r t is outside the Doppler blind region, the target is not always detected but has a detection probability P D (see, e.g., [11]). The actual measurement z t takes the values in the set Z = {R 3 0 /}, where 0 / denotes a missing measurement, i.e.,

B. Bayesian Inference for Recursive Filtering
Let X t = {x 1 , . . ., x t } and Z t = {z 1 , . . ., z t } denote the sequences of the state vectors and measurement vectors up to time step t. With the state space model (1)-(3), Bayesian inference can be drawn to derive the posterior distribution p(X t , m t |Z t ) at each time step t, which in turn can be used to determine both the mode probability p(m t |Z t ) and the mode-conditioned probability p(X t |m t , Z t ). In this subsection, we briefly summarize the general recursive Bayesian filtering for the formulated estimation problem; see [15] for a detailed description.
Given the measurement sequence Z t = {z 1 , . . ., z t }, the recursive Bayesian filtering obtains the probability distribution p(X t , m t |Z t ) at each time step t using the previous distribution p(X t-1 , m t-1 |Z t-1 ) at t -1. First, we note that given p(X t-1 , m t-1 |Z t-1 ), p(X t-1 , m t |Z t-1 ) can be obtained by law of total probability (termed mode interaction in [15]): On the basis of the mode interaction in (4) and the state transition probability distribution p(x t |X t-1 , m t , Z t-1 ) defined in (1), the predicted distribution p(X t , m t |Z t-1 ) can be obtained: Finally, following the Bayesian theorem, the posterior distribution of the state vector can be derived from the predicted distribution by taking into account the current measurement z t : where p(z t |x t ) is the likelihood function that is the probability of z t conditional on the state vector x t . Its definitions are given by the corresponding measurement model, (2) and (3). Hence, following the general Bayesian statistical inference procedure, the distribution p(X t-1 , m t-1 |Z t-1 ) at time step t -1 can be updated to form the posterior distribution p(X t , m t |Z t ) at time step t. Furthermore, the marginal distribution p(x t , m t |Z t ) can be obtained from p(X t , m t |Z t ), upon which the state and mode estimations at time instance t can be worked out.

III. THE PROPOSED ALGORITHM
A. The General Structure Overall, the enhanced particle filtering algorithm developed in this paper is based on (4)- (6), and the filtered state vector at each time step t is derived by the corresponding marginal posterior distribution p(x t , m t |Z t ). First, we focus on (6). When there is a measurement recorded at t, p(z t |x t ) in (6) is given by (2). However, when there is no measurement (i.e., z t = 0 /), the likelihood function p(z t = 0 /|x t ) is given by, according to (3): In general, there is no exact analytical solution for the posterior distribution p(X t , m t |Z t ) for state and mode estimation. Recently, Clark et al. [12] have proposed an efficient GMM-based NRDB algorithm to approximate p(X t , m t |Z t ) on the basis of approximations A-C as outlined in the previous section. In order to address the limitations of the algorithm in [12] and to obtain more accurate state estimation, a novel particle filtering approach is developed as follows.
It is well-known that the quality and efficiency of any particle filters depend on the quality of the chosen importance distribution. Ideally, the importance distribution should be 1) close to the posterior distribution of the state vector and the mode, i.e., p(X t , m t | Z t ), and 2) easy to sample from. The standard choice of importance distribution is the state transition distribution p(x t |X t-1 , m t , Z t-1 ) given in (1). The main drawback of such a choice is that it does not take into account the current measurement z t (see, e.g., [16]) that reveals the state information.
Essentially the new filter to be developed in this paper includes two key elements. First, the EKF and its generalization for interval-censored measurements in [12] are used to construct the importance distributions. In comparison with the standard choice of importance distribution, i.e., the state transition distribution at each state mode, the sampling efficiency is improved for each mode of the multiple models where the information about the current measurement z t has been incorporated. Second, all three approximations in the NRDB algorithm are addressed by utilizing the advantages of particle filtering. This includes 1) the standard mixture reduction technique is avoided by the resampling method and 2) the linearization and Gaussian approximations are corrected through the weighting process in the proposed particle filter.
Specifically, at time t -1, suppose that for each mode s ∈ M, N weighted particles are allocated to the corresponding mode-matched filter where δ (·) is a Dirac delta function. From (4), the mode interaction is given by: Assuming the number of the state modes to be M, from (9), we can see that N × M particles are required to represent p(X t-1 , m t = r|Z t-1 ). The increasing number of particles from N to N × M for probability representation at each time step will make the number of particles grow in an exponential way as time t becomes large.
Blom and Bloem [15] developed a resampling method to address this problem, where the particles are resampled from the above distribution, conditional on the model m t : Because N is usually taken as a large number, the error caused by the above resampling process is small. In contrast, the number of modes in the NRDB algorithm is usually small or medium, and hence the standard mixture reduction technique used in [12] may not be able to approximate well a multimodal posterior distribution.
Equations (5), (6), and (10) can be applied to derive the posterior distribution. See [15] for details. As mentioned earlier, this posterior is analytically intractable, and hence Monte Carlo methods are usually used and new particles need to be drawn. In the particle filtering, this step is undertaken using an importance distribution.
Specifically, for each mode r, importance sampling is used to estimate the desired posterior distribution. Now suppose that a new sample x r,k t of the current state vector x t is drawn from the importance distribution The obtained new state vector samples for time step t and the corresponding X together with the corresponding weights w r,k t , are used to represent the posterior distribution at time step t: upon which the marginal distribution p(x t , m t = r|Z t ) for the current state vector can be obtained. According to whether GMTI measurements are recorded or not, different importance distributions q r (x t |X r,k t−1 , Z t ) will be used in the proposed method. This will be investigated in detail below.

B. Particle Filtering with Measurements Recorded
When the measurement at time step t is recorded, we have z t = y t given by (2). The EKF is applied to construct the importance distribution for every X r,k t−1 . As in [6], we assume that in the local region near the state vectorx r,k t−1 , the mode-conditioned distribution p(x t-1 |m t = r, Z t-1 ) is approximated as a Gaussian distribution N (x t−1 |x r,k t−1|t−1 , P r,k t−1|t−1 ) with meanx r,k t−1|t−1 and covariance matrix P r,k t−1|t−1 . The EKF estimates the Gaussian approximation of p(x t |m t = r, Z t ) in a local region, which is taken as the importance distribution q r (x t |X r,k t−1 , Z t ) for the new particle generation. The EKF scheme follows the standard Kalman filtering procedure by approximating the nonlinear measurement model via its first-order Taylor expansion. It is divided into the following prediction and correction steps: CORRECTIONỹ wherex r,k t−1|t−1 and P r,k t−1|t−1 represent the initial mean and covariance matrix associated with the particleX r,k t−1 . Q r t−1 and R t are the covariance matrices of the process noise w t-1 (m t ) with m t = r in (1) and the measurement noise n t in (2), respectively. F r t is F (m t ) evaluated at m t = r. H r,k t is the gradients of h(x t ) evaluated at the point of x r,k t|t−1 : The obtained Gaussian distribution N (x t |x r,k t|t , P r,k t|t ) is taken as the importance distribution q r (x t |X r,k t−1 , Z t ) for generating new state particles. Clearly, the measurement information is incorporated for constructing N (x t |x r,k t|t , P r,k t|t ) so that the generated particles are more likely in the high measurement likelihood region.
It should be noted that (16) involves the inverse of a matrix for every particleX r,k t−1 , and thus potentially it requires a substantial amount of computational time. In order to reduce the computation cost, we set every initial covariance P r,k t−1|t−1 (k = 1, . . ., N) to be equal to a preset matrix τ diag{1, 1,τ 0 , τ 0 } with τ and τ 0 being two tuning parameters. In addition, the gradient of h(x t ) is evaluated at a common value, 1 N N k=1x r,k t|t−1 ,which are in turn used to approximate H r,k t for each k. In doing so, only one matrix inverse is needed for estimating the importance distribution for each particle.

C. Particle Filtering Without Measurements
When there is no measurement recorded at a time step t, i.e., z t = 0 /, either the target radial velocity is within the Doppler blindness constraint region (|v r t | ≤ κ) or the target is not detected (with probability 1 -P D ). For each mode m t = r, the mode-conditioned probability of the state vector x t follows: . (20) Similar to the scenario that the measurements are recorded, we use a mixture of two Gaussian distributions to approximate the local distribution of p(x t |m t = r, Z t ). This Gaussian mixture distribution is taken as the importance distribution for generating a new state particle. This will be discussed in detail below.
For the first term in (20), the probability p(x t |m t = r, Z t-1 ) represents the predicted distribution based on the previous measurements Z t-1 . From the approximated initial Gaussian distribution N (x t−1 |x r,k t−1|t−1 , P r,k t−1|t−1 ) near the particlex r,k t−1 (as in the previous subsection), the predicted distribution can be approximated as a Gaussian distribution with the meanx r,k t|t−1 and covarianceP r,k t|t−1 . The recursive formulae are identical to the ones given previously.
We now turn to consider the second term of (20). Note that now the only information we know is that the target radial velocity v r t falls into a given interval, i.e., |v r t | ≤ κ. This is termed an interval-censored problem in the literature. To deal with this interval-censored problem, we follow [12] and use the results in [14] to work out the mathematical expectation and variance conditional on the interval-censored measurement v r t (see Appendix A for the detailed results) as follows.
First, we approximate the target radial velocity by the first-order Taylor expansion around the predicted statê x r,k t|t−1 as: From (21), the conditional meanx r,k t|t and covariance P r,k t|t of p(x t |m t = r, Z t−1 , |v r t | ≤ κ) can be calculated using the formulae in Appendix A: This procedure is a generalization of the EKF to the problem of state estimation based on an interval-censored measurement. On the basis of this generalized EKF, the distribution p(x t |m t = r, Z t−1 , |v r t | ≤ κ) can be approximated as N(x t |x r,k t|t , P r,k t|t ). In addition, based on (21), the probability that the measured range rate is within the Doppler blind zone, i.e., p(|v r t |≤ κ|m t = r, Z t−1 ), is approximated by Finally, on the basis of the above analysis, we choose the importance distribution for each particleX r,k t−1 as:

REMARKS
1) The above importance distribution is a Gaussian mixture distribution, and thus it is straightforward to draw particles from this distribution.
2) It can be seen from the construction of the importance distribution (29) that, when no measurement is recorded, the relevant information is incorporated by taking into account the Doppler blindness constraint |v r t | ≤ κ. The constructed importance distribution is a local approximation of the desired distribution p(x t |m t = r, Z t ) so that the particles of p(x t |m t = r, Z t ) can be generated from the importance distribution.
3) Similar to the scenario that the measurements are recorded, we set the initial covariance P r,k 1, ..., N) to be equal and each H r,k t is approximated by the first-order gradient of the radial velocity evaluated at the averaged predicted means. In doing so, only oneσ r,k and one K r,k t are needed to calculate for all the particles, and hence the computational cost is reduced.

D. Correcting the Approximation Errors
It can be seen from the previous subsections that several approximations were made when deriving the importance distributions. The corresponding approximation errors can be corrected in the proposed particle filtering algorithm through a weighting operation.
Suppose that we have drawn a new state vector particle x r,k t from the importance distribution q r (x t |X r,k t−1 , Z t ) and hence formed a new particle X r,k t = [X r,k t−1 , x r,k t ]. The corresponding weight for each new particle X r,k t is given by the following ratio: . (30) Each particle X r,k t is assigned a suitable weight w r,k t . This results in an estimate of the posterior distribution, (11). Two things are worth noting from (30). 1) In some practical problems, the state transition probability distribution given by the state equation (1) is rank-deficient and p(x r,k t |X r,k t−1 , m t = r, Z t ) could not be calculated directly. We apply the QR decomposition to 2) We have used a novel importance distribution that is estimated from either the EKF or the generalized EKF to incorporate the GMTI measurement information, and therefore the weight in (30) is different from that used in [15]. In [15], the importance distribution is chosen as the system transition distribution p(x t |x t-1 , m t = r) and hence w r,k t ∝ p(m t = r|Z t−1 )p(z t |x r,k t ).

IV. SIMULATION STUDY
In this section, we use a simulation study to evaluate the numerical performance of the developed filter. A simulation study was carried out for a single target-tracking scenario in which a GMTI sensor mounted on an airborne platform was employed to track the motion of a moving ground vehicle. The target moving eastbound started at a constant speed of 10 m/s maintained for 180 s before it accelerated at a rate of 1 m/s 2 up to a speed of 25 m/s. After travelling at this constant speed for 180 s, the target started to decelerate for 25 s until it came to a standstill. The target remained stationary for 60 s, before accelerating again to a speed of 15 m/s. Target speed as a function of time is plotted in Fig. 1. The same simulation scenario was also investigated in [12].
The sensor platform travelled northbound at a constant speed of 120 m/s and at an altitude of 10 km. The moving sensor platform took noisy measurements (including the 3-dimensional range, azimuth angle, and range rate) of the ground-moving target every 5 s. The movements of the target and sensor platform are displayed in Fig. 2.

A. The State and Measurement Models
Because the vehicle moved in several different maneuvering styles, we applied a multiple models scheme to describe the movement of the target. Following [12,19], we define a low-intensity nearly constant velocity (LINCV) model and a high-intensity nearly constant velocity (HINCV) model. In addition, a stop model based on the nearly constant position model in [20] is used to model the scenario of vehicle stopping.
To better compare the proposed algorithm with the GMM-based NRDB algorithm in [12], we used the same method to discretize the continuous-time model of vehicle dynamics and assumed the following general form of the state equation for every model (see, e.g., [21], for other forms of the discretized state equation in the literature): where m t = {1, 2, 3} represents the model index: m t = 1 for the LINCV model, m t = 2 for the HINCV model, and m t = 3 for the stop mode. A Gaussian distribution is assumed for the 2 × 1 noise vector w t-1 (m t ) with zero means and mode-related covariance matrix diag{σ 2 x (m t ), σ 2 y (m t )} with 1) σ 2 x (m t ) = σ 2 y (m t ) = 0.05 2 (m/s 2 ) 2 for m t = 1; 2) σ 2 x (m t ) = σ 2 y (m t ) = 0.5 2 (m/s 2 ) 2 for m t = 2; and 3) σ 2 x (m t ) = σ 2 y (m t ) = 0.005 2 (m/s) 2 for m t = 3.
The state transition matrix F(m t ) and constant matrix G(m t ) for m t = {1, 2} are defined as: where the evolution of the position is modeled by adding nonzero process noise as in [20]. In addition, the velocity is forced to be zero as in [22,23]. The relevant model transition probabilities were calculated according to the method described in [7]. The resulting state mode transition matrix is given by: The measurement model, as discussed in Section II, includes two scenarios: measurements recorded and not recorded. As illustrated in Fig. 3, when a measurement is recorded and the target's radial velocity is outside the Doppler blindness constraint region (-κ, + κ) (here we followed [12] and set κ to be 3 m/s), the corresponding likelihood function is Gaussian as given by (2). In the simulation study, we followed [12] and set the parameters of the measurement equation σ r and σ θ as 20 m and 0.001 rad, respectively. On the other hand, when no measurement is recorded, the measurement model is given by (7). Finally, throughout the simulation study, the two tuning parameters τ and τ 0 were set as 1 and 0.1, respectively.

B. Tracking Performance Analysis
Based on the state and measurement models, the proposed algorithm was applied in the tracking of the moving target in the set scenario and the corresponding tracking performance was evaluated. As in [12], we assumed the initial target state distribution followed a single Gaussian prior density obtained by the single-point track initialization algorithm [24,25] for the proposed approach and for the other algorithms used for comparison below.
The maneuvering types of the moving target at different time points are shown in Fig. 4a, and the estimated model probabilities using the proposed algorithm are given in Fig. 4b. It can be seen that in the majority of time instances, the model that corresponded to the actual movement type had the largest probability. Hence, the maneuvering characteristics of the moving target could be correctly reflected by the proposed algorithm in the simulation study. Next, we analyze the tracking accuracy of the proposed method and draw a comparison with some other state-of-the-art approaches, including the multiple model version of the noise-related Doppler blind mixture filter (NRDB-MM) [12] and the multiple model particle filtering (MMPF) approach. Here the multiple-model-based approaches were chosen for a fair comparison because our proposed approach considers multiple state models. In particular, the NRDB-MM algorithm in [12] was used as the benchmark algorithm for comparison purposes, whereas the MMPF approach is a widely used algorithm in the literature (see, e.g., [4,11,26]). The parameter for the NRDB-MM was set exactly the same as in [12]. For the two particle filtering-based methods, i.e., the proposed method and the MMPF, the number of particles corresponding to each state model (denoted as N) was chosen at different levels so that we can assess its impact on the accuracy of the estimation.
For each approach, 100 Monte-Carlo simulation experiments were carried out and the RMSEs averaged over the 100 simulations at different time instances were calculated. We first focus on the average RMSE obtained by the different methods, as displayed in Fig. 5. It can be seen from Fig. 5 that: 1) The proposed method achieved the best performance with the smallest errors during most of the time instances than its counterparts, especially for the sample indexes between 80 and 92 when the target stopped.
2) Compared with the MMPF approach, the proposed method was less sensitive to the number of particles. Fig. 5 shows that the errors for the MMPF method were substantially increased when the number of particles reduced from N = 2500 to N = 1000. In contrast, similar performances were obtained for the proposed method with different particle sizes.
There are several reasons for these differences. First, in comparison with the NRDB-MM algorithm, the proposed approach adopts the particle filtering-based method, and therefore it avoids the approximations made in the GMM-based NRDB approach in [12], leading to a performance improvement over the NRDB-MM algorithm. Second, compared the MMPF approach, a novel sampling method is proposed to incorporate either the measurement information (when measurements are recorded) or knowledge on the Doppler blindness constraint (when no measurements are recorded) for the construction of the importance distributions. Effective particles can thus be sampled from the constructed importance distributions, which leads to a more accurate state estimation even with a relatively small particle size.
Finally, we focus on vehicle tracking when the vehicle was standstill and hence was within the Doppler blindness constraint region. This is a challenging scenario because no measurements were recorded during the vehicle stopping interval. For this end, we considered the same parameter settings as did in [12], including different detection probabilities P D and the standard deviations of range rate measurement σṙ . A comprehensive evaluation was performed under these parameter settings to compare the different approaches.
Under each parameter setting, the RMSEs during the time interval when the vehicle stopped for different approaches were calculated. Tables I and II display the average RMSEs over the 100 simulation experiments. It can be seen from the two tables that the proposed approach outperformed the other methods with the smallest tracking errors. This shows that for the scenario where the vehicle stopped and no measurements were recorded, the proposed approach could most efficiently utilize knowledge on the Doppler blindness constraint for vehicle tracking.

C. Computational Costs
In this subsection, we briefly evaluate the computational cost of the different methods. The execution was performed on a PC with 3.40 GHz processing speed and 8.00 GB memory using Matlab 2013(a). For the scenario under investigation, the vehicle   moved 700 s and the measurement was taken every 5 s, so it was required that the computation time for processing one sample be less than 5 s. In total there were 700/5 = 140 sample time indexes (time steps). For each method considered in the simulation study, the execution time and the average processing time for each sample were recorded, as displayed in Table III. We can see that the computational costs of all the algorithms for processing one sample are far less than 5 s and fast enough to make the tracking algorithm run in practice. Further reduction of the computation time can be achieved by implementing the algorithm in a C/C++ programming environment.

V. CONCLUSIONS
In this paper, a novel particle filtering approach for GMTI tracking has been developed. The proposed particle filter makes use of the EKF and its generalized version in [12] to construct importance distributions for efficient generation of particles. In comparison with the GMM-based NRDB algorithm recently developed in [12], the three approximations outlined in the Introduction section can be satisfactorily addressed by the proposed algorithm, so that more accurate state filtering can be obtained. On the other hand, comparing to general-purpose particle filters, where the importance distribution is usually chosen as the state transition distribution, the information on both the measurements and the Doppler blindness constraint can be directly incorporated in the proposed method when forming the importance distributions, as a result, substantially enhancing the quality of the particle filter. The numerical results in the Monte Carlo simulation have verified that the proposed method outperformed the existing methods with an acceptable computational cost.

APPENDIX A. MEAN AND VARIANCE CONDITIONAL ON AN INTERVAL-CENSORED MEASUREMENT
We summarize the results on the mathematical expectation and variance conditional on some interval-censored measurements below; see [14] for details. Consider a measurement m given by m = q T x + w with w ∼ N (w|0, σ 2 ). Assume that the prior distribution of the state vector x is Gaussian, N (x|x 0 , P 0 ). Let A denote an interval A = [a, b].

APPENDIX B. RANK-DEFICIENT PROCESS NOISE
In many practical problems, the state transition probability distribution given by the state equation is rank-deficient (see, e.g., [21]). Consider the following state transition equation: where x t is an n-dimensional state vector. G t is an n × q matrix with n > q and w t-1 (m t ) ∈ R q . In this case, x t is constrained in a vector space with the deficient rank q < n. As a consequence, if particles x i t are sampled from an importance distribution in R n , they may fall outside the required vector space and the corresponding state transition probability then could not be calculated. Clearly, an accept-reject method where samples x i t outside the given space are discarded is very inefficient.
In order to solve this rank-deficient problem, we perform the QR decomposition for the matrix G t as: where W t = [ t , 0] T . t is a matrix with rank q and U t is an orthogonal matrix. We define: . (44) From (44), we can see that by transforming x t to s t , the state transition function is split into a stochastic part s 1,t and a deterministic part s 2,t . The deterministic part s 2,t can be directly worked out from the previous state value. Therefore, we apply the developed method only to generate particles that are related to the stochastic part s 1,t . The obtained particles for s 1,t and s 2,t are then transformed back to the original scale.