*Article* **Reweighted Robust Particle Filtering Approach for Target Tracking in Automotive Radar Application**

**Qisong Wu 1,2,\*, Lingjie Chen 1, Yanping Li 1, Zijun Wang 1, Shuai Yao <sup>1</sup> and Hao Li <sup>3</sup>**


**Abstract:** In view of the decline of filtering accuracy caused by measured outliers in target tracking application, a novel reweighted robust particle filter is proposed to acquire accurate state estimates in an automotive radar system. To infer the importance of each entry in the multidimensional contaminated measurement vector, we employ a weight vector, which follows a Gamma distribution, to model the measured noise and carry out accurate state estimates. Additionally, the particle filter method is employed to perform approximate posterior inference of state estimates in the nonlinear model. The Cramer–Rao lower bound is provided for the performance evaluation of the proposed method. Both simulation and experimental results demonstrate the superiorities of the proposed algorithm over other robust solutions.

**Keywords:** outlier-robust; particle filter; tracking; automotive radar

#### **Citation:** Wu, Q.; Chen, L.; Li, Y.; Wang, Z.; Yao, S.; Li, H. Reweighted Robust Particle Filtering Approach for Target Tracking in Automotive Radar Application. *Remote Sens.* **2022**, *14*, 5477. https://doi.org/10.3390/ rs14215477

Academic Editors: Zhihuo Xu, Jianping Wang and Yongwei Zhang

Received: 7 September 2022 Accepted: 20 October 2022 Published: 31 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

In the field of autonomous driving, automotive radars have become the core sensors due to their all-time and all-weather environmental perception capabilities. Effective target tracking is essential for the autonomous driving systems to formulate correct and timely driving strategies. By clustering the point cloud data of Millimeter Waveform (MMW) radar, the target can be detected and located in the automotive application scenario. However, the occurrence of measurement abnormalities in real-world applications is a common phenomenon that may arise due to reasons such as unexpected disturbances, temporary sensor failures, and calibration errors [1]. A typical example is to track a walking pedestrian on the road through an MMW radar sensor in action. The occlusion of roadside obstacles, such as parked cars, walls, and so on, can corrupt the range measurement, and glint, speckle, and phase noise can contaminate radar data [2–4]. The input of most tracking algorithms is the result of target detection. If the target is suddenly lost or temporarily occluded, the target trajectory will be interrupted [5]. Moreover, the clutter cannot be completely filtered out of radar measurements can lead to false detections and causes measurement abnormalities [6]. In addition, the non-line-of-sight (NLOS) multipath propagation of the transmitted electromagnetic wave is also one of a source of measurement abnormalities [7]. The Kalman filter (KF) offers an optimal estimation method for the linear dynamic model with a Gaussian state and measurement noise [8]. However, owing to the complex electromagnetic propagation and changeable environment, the measurement noise is a typical non-Gaussian in the MMW radar sensor and would result in measurement abnormalities, also named outliers. The presence of these outliers may degrade the performance of standard KF, thus yielding poor estimates of the state of interest [9].

A number of efforts have been devoted to addressing the issue. Outliers rejection methods are used in robust state estimation and object tracking [10]. A robust filter approach in [11] has been developed by viewing measurements as clearly distinguished outliers from

normal values using Projection Statistics (PS). In [12–14], classical recursive approaches are converted into batch-mode regression form and iteratively reweighted least squares algorithms are used to estimate the states. A measurement-weighting-based whitening procedure is designed to distinguish outliers from normal measurements to improve the state estimation accuracy [15]. Another class of approaches is based on the heavy-tailed noise model, such as Student's *t* distributions and Laplace distributions [16–19]. A scalar weight, which follows a Gamma distribution, is introduced for each observed sample such that the variance matrices of measurement noise can be reassessed in the state estimation [20,21]. Filtering and smoothing algorithms, which are much more robust yet only slightly more involved than the standard inference algorithm, are derived by assigning Wishart distribution to the noise covariance matrix [22]. By choosing the inverse Wishart priors, the measurement noise covariance matrices can be inferred based on the variational Bayes (VB) approach [23].

In this paper, a reweighted robust particle filter is proposed to acquire more accurate state estimation with robustness to outliers in the nonlinear model. Motivated by the aforementioned methods in [20,21], we generalize models to the more comprehensive case of multidimensional outlier measurements. Specifically, a weight vector, which follows conditional independent Gamma distribution, is employed to reassess the contribution of each entry in the multidimensional outlier measurement. The proposed method can be regarded as a reweighted robust filter. The particle filter (PF) method is utilized to perform the approximate posterior inference of state estimates. The Cramer–Rao bound of the proposed method is also provided for performance analysis. Both simulations and experiments demonstrate the effectiveness of the proposed method.

Notations: We use lower-case (upper-case) bold characters to denote vectors (matrices). (·)*<sup>T</sup>* denotes the transpose of a matrix or vector. diag(**x**) represents a diagonal matrix that uses the elements of **x** as its diagonal elements. *p*(·) denotes the probability density function (pdf). *x* ∼ N (*a*, *b*) denotes that random variable *x* follows a Gaussian distribution with mean *<sup>a</sup>* and variance *<sup>b</sup>*. **<sup>x</sup>** ∈ R*<sup>N</sup>* denotes a real-value vector **<sup>x</sup>** with the dimension of *<sup>N</sup>* and **<sup>x</sup>** ∈ R*N*×*<sup>M</sup>* denotes a real-value matrix **<sup>x</sup>** with the dimension of *<sup>N</sup>* × *<sup>M</sup>*.

#### **2. Proposed Model**

*2.1. Problem Formulation*

A typical nonlinear discrete-time state-space model can be expressed as

$$\begin{aligned} \boldsymbol{\theta}\_{k} &= \boldsymbol{f}(\boldsymbol{\theta}\_{k-1}) + \mathbf{u}\_{k} \\ \mathbf{z}\_{k} &= \boldsymbol{h}(\boldsymbol{\theta}\_{k}) + \mathbf{n}\_{k} \end{aligned} \tag{1}$$

where *<sup>θ</sup><sup>k</sup>* ∈ R*<sup>N</sup>* and **<sup>z</sup>***<sup>k</sup>* ∈ R*<sup>M</sup>* represent the state vector and measurement vector at time step *k*, *f*(·) and *h*(·) denote nonlinear state and measurement functions, respectively, **<sup>u</sup>***<sup>k</sup>* ∈ R*<sup>N</sup>* and **<sup>n</sup>***<sup>k</sup>* ∈ R*<sup>M</sup>* are the zero-mean Gaussian white process noises satisfying **u***<sup>k</sup>* ∼ N (**0**, **Q**) and **n***<sup>k</sup>* ∼ N (**0**, **R**), and **Q** and **R** are the process noise covariance matrix and measurement noise covariance matrix, respectively. The model in Equation (1) is just regarded to be a typical nonlinear system model, and the system state can be estimated by a nonlinear Kalman filter [24].

In real automotive radar applications, typical representative scenarios are often encountered in which the target of interest may be obscured by other interferences or clutters, and the acquired measurement vector, such as range, velocity, and azimuth, in the automotive radar system would fall outside a valid range and become outliers, since the acquired measurement samples at this time directly come from interferences. The weighted robust Kalman filter method has been developed by introducing a scalar *wk* into the measurement noise covariance matrix **R** for the improved estimates of the updated *θ<sup>k</sup>* [21]. However, the main limitation of this work is that even if it allows the treatment of multidimensional measurement, a scalar weight is considered for the entire observation vector, thus equally treating all the entries regardless of which entry is corrupted. On the one hand, the entries

in the measurement vector represent diverse physical quantities respectively in the real radar system and have diverse signal-to-noise ratios. On the other hand, the variance of each entry is often influenced by the automotive radar system, such as the wider bandwidth radar system, resulting in a smaller range variance and smaller number of receiver channels, leading to a poorer azimuth measurement. Therefore, introducing a scaler in the measurement vector may be not an optimal way to treat them equally.

#### *2.2. Generative Model*

To address this issue, we introduce a novel Bayesian algorithm that treats the weights associated with each entry in the measurement vector data probabilistically. In particular, we introduce a weight vector **w***<sup>k</sup>* to learn the dynamical noise covariance matrix, and a reweighted robust particle filter method is proposed to acquire improved state-space estimates in a nonlinear model in the presence of contaminated measurements. We modify the measurement noise equation as follows by introducing a weight vector **w***k*,

$$\mathbf{n}\_k \sim \mathcal{N}\left(\mathbf{0}, \mathbf{W}\_k^{-1}\mathbf{R}\right),\tag{2}$$

where **<sup>W</sup>***<sup>k</sup>* ∈ R*M*×*<sup>M</sup>* is a diagonal matrix with the weight vector **<sup>w</sup>***<sup>k</sup>* = [*wk*,1, ··· , *wk*,*M*] *T* ∈ <sup>R</sup>*<sup>M</sup>* on its diagonal at time step *<sup>k</sup>*, and **<sup>W</sup>**−<sup>1</sup> *<sup>k</sup>* is the inverse of the matrix **W***k*. The Gamma prior distribution is introduced for the elements of **w***<sup>k</sup>* to ensure that they remain positive and is written by

$$\begin{aligned} \mathbf{z}\_k | \boldsymbol{\theta}\_{k'} \mathbf{W}\_k &\sim \mathcal{N} \left( h(\boldsymbol{\theta}\_k), \mathbf{W}\_k^{-1} \mathbf{R} \right) \\ \boldsymbol{\theta}\_k | \boldsymbol{\theta}\_{k-1} &\sim \mathcal{N} (f(\boldsymbol{\theta}\_{k-1}), \mathbf{Q}) \\ w\_{k, \text{m}} &\sim \text{Gamma} (w\_{k, \text{m}} | v\_{\text{m}} / 2, v\_{\text{m}} / 2), \end{aligned} \tag{3}$$

where *m* = 1, ··· , *M*, *vm* is the *m*th hyperparameter in the vector **v** for the Gamma distribution. According to Equation (3), the modified model has the capability of dynamically learning the measurement noise covariance and adjusting the weights of each entry in the measurement vector for improved estimates of the state vector with effectively alleviating outliers.

#### *2.3. Reweighted Robust Particle Filter*

In this subsection, we develop a reweighted robust particle filter using statistical inference over the aforementioned model. The particle filter, utilizing the Monte Carlo method, is selected to provide approximate posterior estimation in the nonlinear model [25]. The state estimation generally focuses on the expectation calculation of posterior distribution,

$$\mathbb{E}[\mathcal{g}(\theta)|\mathbf{z}\_{1:k}] = \int \mathbb{g}(\theta) p(\theta|\mathbf{z}\_{1:k}) d\theta \tag{4}$$

where *g*(*θ*) is the function of the variable *θ*, and *p*(*θ*|**z**1:*k*) represents the posterior probability density function for *θ* given the measurement vector data **z**1:*k*. However, the posterior distribution *p*(*θ*|**z**1:*k*) is often difficult to analytically acquire due to the manipulation of the multiple integration. The importance sampling (IS) is often used to perform the posterior distribution by introducing a nominal distribution *π*(*θ*|**z**1:*k*) that often utilizes a Gaussian distribution,

$$\begin{split} \mathbb{E}[\mathcal{X}(\theta)|\mathbf{z}\_{1:k}] &= \int \mathcal{g}(\theta) \frac{p(\theta|\mathbf{z}\_{1:k})}{\pi(\theta|\mathbf{z}\_{1:k})} \pi(\theta|\mathbf{z}\_{1:k}) d\theta \\ &= \int \mathcal{g}(\theta) \frac{p(\mathbf{z}\_{1:k}|\theta) p(\theta)}{p(\mathbf{z}\_{1:k}) \pi(\theta|\mathbf{z}\_{1:k})} \pi(\theta|\mathbf{z}\_{1:k}) d\theta \\ &= \frac{\int \mathcal{g}(\theta) \frac{p(\mathbf{z}\_{1:k}|\theta) p(\theta)}{\pi(\theta|\mathbf{z}\_{1:k})} \pi(\theta|\mathbf{z}\_{1:k}) d\theta}{\int \frac{p(\mathbf{z}\_{1:k}|\theta) p(\theta)}{\pi(\theta|\mathbf{z}\_{1:k})} \pi(\theta|\mathbf{z}\_{1:k}) d\theta}. \end{split} \tag{5}$$

The particle weight *ω*(*θ*) is defined as

$$
\omega(\theta) = \frac{p(\mathbf{z}\_{1:k}|\theta)p(\theta)}{\pi(\theta|\mathbf{z}\_{1:k})}.\tag{6}
$$

The Equation (5) can be rewritten by introducing the *ω*(*θ*):

$$\begin{split} \mathbb{E}[\mathcal{g}(\boldsymbol{\theta})|\mathbf{z}\_{1:k}] &= \frac{\int \mathcal{g}(\boldsymbol{\theta})\omega(\boldsymbol{\theta})\pi(\boldsymbol{\theta}|\mathbf{z}\_{1:k})d\boldsymbol{\theta}}{\int \omega(\boldsymbol{\theta})\pi(\boldsymbol{\theta}|\mathbf{z}\_{1:k})d\boldsymbol{\theta}}\\ &\approx \frac{\sum\_{s=1}^{S} \mathcal{g}\left(\boldsymbol{\theta}^{(s)}\right)\omega(\boldsymbol{\theta}^{(s)})}{\sum\_{s=1}^{S} \omega(\boldsymbol{\theta}^{(s)})} \end{split} \tag{7}$$

where *s* = 1, ... , *S*, and *S* is the total number of sampling particles. According to the Markov properties of the model, the particle weights *ω*(*s*) *<sup>k</sup>* can be deduced in the recursion form

$$
\omega\_k^{(s)} \propto \frac{p\left(\mathbf{z}\_k | \boldsymbol{\theta}\_k^{(s)}\right) p\left(\boldsymbol{\theta}\_k^{(s)} | \boldsymbol{\theta}\_{k-1}^{(s)}\right)}{\pi \left(\boldsymbol{\theta}\_k^{(s)} | \boldsymbol{\theta}\_{0:k-1}^{(s)}, \mathbf{z}\_{1:k}\right)} \omega\_{k-1}^{(s)}.\tag{8}
$$

Generally, *<sup>π</sup>*(*θk*|*θ*0:*k*−1, **<sup>z</sup>**1:*k*) = *<sup>p</sup>*(*θk*|*θk*−1) is often used in the importance sampling, and thus the recursion in Equation (8) can be reduced into

$$
\omega\_k^{(s)} \propto p\left(\mathbf{z}\_k | \boldsymbol{\theta}\_k^{(s)}\right) \boldsymbol{\omega}\_{k-1}^{(s)}.\tag{9}
$$

According to Equation (9), it is obvious that the problem eventually converts to computation of *p*(**z***k*|*θk*), which can be achieved from the marginal likelihood utilizing Bayesian inference. Considering the fact that **W***<sup>k</sup>* is independent of *θk*, we have

$$\begin{split} &p(\mathbf{z}\_k|\boldsymbol{\theta}\_k) \\ &= \int p(\mathbf{z}\_k|\boldsymbol{\theta}\_{k'}\mathbf{W}\_k) p(\mathbf{W}\_k|\boldsymbol{\theta}\_k) d\mathbf{W}\_k \\ &= \prod\_{m=1}^M \int p(z\_{k,m}|\boldsymbol{\theta}\_{k'}\boldsymbol{w}\_{k,m}) p(\boldsymbol{w}\_{k,m}) d\boldsymbol{w}\_{k,m} \\ &= \prod\_{m=1}^M \int\_0^\infty \mathcal{N}\left(z\_{k,m} \mid \boldsymbol{\mu}\_{k,m}, \frac{r\_{m,m}}{r\nu\_{k,m}}\right) \text{Gam}\left(w\_{k,m} \mid \frac{v\_m}{2}, \frac{v\_m}{2}\right) d\boldsymbol{w}\_{k,m} \\ &= \prod\_{m=1}^M \frac{\Gamma\left(\frac{v\_m}{2}+\frac{1}{2}\right)}{\Gamma\left(\frac{v\_m}{2}\right)} \left(\frac{\lambda\_m}{\pi \upsilon\_m}\right)^{\frac{1}{2}} \left[1 + \frac{\lambda\_m(z\_{k,m}-\mu\_{k,m})^2}{\upsilon\_m}\right]^{-\frac{\mu\_m}{2}-\frac{1}{2}} \end{split} \tag{10}$$

where *μk*,*<sup>m</sup>* is the *m*th element of the vector *h*(*θk*); *λ<sup>m</sup>* = *r*−<sup>1</sup> *<sup>m</sup>*,*<sup>m</sup>* is the *m*th element on the diagonal of matrix **R**<sup>−</sup>1. The entire reweighted robust particle filter (RR-PF) method is summarized as Algorithm 1.

**Algorithm 1** Reweighted Robust Particle Filter **Input:** {*θ* (*s*) *<sup>k</sup>*−<sup>1</sup> : *<sup>s</sup>* <sup>=</sup> 1, . . . , *<sup>S</sup>*} **Output:** *θk*, {*θ* (*s*) *<sup>k</sup>* : *s* = 1, . . . , *S*} 1: Draw new samples: *θ*˜(*s*) *<sup>k</sup>* ∼ *p θk*|*θ* (*s*) *k*−1 ,*s* = 1, . . . , *S* 2: Calculate particles weights: *ω*˜ (*s*) *<sup>k</sup>* in Equations (9) and (10) 3: Perform resampling with *θ*˜(*s*) *<sup>k</sup>* and *ω*˜ (*s*) *<sup>k</sup>* to acquire *θ* (*s*) *<sup>k</sup>* in Equation (7) 4: Estimate system state *θk* <sup>=</sup> <sup>∑</sup>*<sup>S</sup> <sup>s</sup>*=<sup>1</sup> *θ* (*s*) *<sup>k</sup>* /*S* 5: **return** *θk*, {*θ* (*s*) *<sup>k</sup>* : *s* = 1, . . . , *S*}

#### *2.4. Algorithm Analysis*

In this subsection, detailed analysis is provided to illustrate the reason why the proposed method is more robust to the outlier compared with the conventional heavytail multivariate Student's *t*-distributed (CMSt) [26,27]. For comparison, we provide the probability density function (pdf) of the CMSt as

$$p(\mathbf{z}\_k|\boldsymbol{\theta}\_k) \propto \frac{\Gamma\left(\frac{\upsilon}{2} + \frac{M}{2}\right)}{\Gamma\left(\frac{\upsilon}{2}\right)} \left[1 + \frac{(\mathbf{z}\_k - \boldsymbol{\mu})^T \mathbf{R}^{-1} (\mathbf{z}\_k - \boldsymbol{\mu})}{\upsilon}\right]^{-\frac{M}{2} - \frac{\upsilon}{2}}.\tag{11}$$

We take 2D pdfs with *M* = 2 as an example to illustrate the roles of the robustness to outliers. In this example, **z***<sup>k</sup>* = [*xk*, *yk*] *<sup>T</sup>*, *μ* = [0, 0] *<sup>T</sup>*, **R**−<sup>1</sup> = diag([1, 1]), and **v** = [*v*, *v*] *<sup>T</sup>* = [0.5, 0.5] *<sup>T</sup>* are shown in Figure 1. It is observed that both pdfs exhibit heavy-tailed properties, which indicates their robustness to outliers. The differences lie in that there exist quite prominent ridges along the *x*-axis and *y*-axis in the proposed pdf, which makes the model more adapted to the changes in the independent variable in the direction perpendicular to the ridges than CMSt.

**Figure 1.** Pdfs of the proposed model and CMSt in the 2D example. (**a**) Pdf plot of the proposed model. (**b**) Pdf plot of the CMSt.

In order to intuitively show the role of the proposed pdf in the particle filter, we illustrate the evolution of particle weights. In the particle filter, at (*k* − 1)th time step, the particles are acquired with an identical weight to estimate *θk*−1 based on the resampling technique, shown in Figure 2a. It is observed that these particle balls have the same sizes, and thus have identical weights, although the particle unreliability is different, as shown in the middle surface in Figure 2a. It is also found that the blue particle ball with a small value in the color bar denotes that this particle is highly reliable, however, the yellow ball with a large value suggests an unreliable particle. At the *k*th time step, when the new

measurement vector data *z<sup>k</sup>* are provided, particle weights in the *k*th time step will evolve according to Equation (9). The evolved result of particle weights in Equation (10) is shown on the top surface, and the result in Equation (11) is provided in the bottom surface.

**Figure 2.** The evolution of particle weights. (**a**) 3D plots of evolution of particle weights. (**b**) slice plots of evolution of particle weights. (Particle weights are denoted by the sizes of their corresponding particle balls, and the particle unreliability is depicted and illustrated by the colors.)

To clearly illustrate the roles of the pdfs in Equation (10) compared with that in Equation (11), we provide slice plots of the pdfs the along *y*-axis. It is observed that the plot of the sharper pdf in the proposed method would encourage the contributions from the reliable particle, and it would weaken the particle weight concerning these unreliable particles, shown in the top plot in Figure 2b. It is found that the blue particle ball has a higher weight due to its high reliability in the current time step, whereas the yellow one has a lower weight in the proposed model compared with those in the bottom pdf of the CMSt. In conclusion, the proposed method encourages the contributions from the particle weights of reliable particles and weakens the contributions from the unreliable particles by discriminatively treating the observed vector data. This could be the reason why the proposed method is more robust to outliers and acquires improved state estimates compared with the CMSt method.

#### *2.5. Cramer–Rao Bound for the Model*

To further explore the performance bound of the proposed method, CRLB analysis is provided in this section. We first define a latent parameter set concerning the random parameter set p(**Θ**); the mean-square error matrix is generally served as a measure of error and is given by

$$\Sigma\_{\mathfrak{c}} = \text{MSE}(\boldsymbol{\Theta}) \triangleq \mathbb{E}\left[ \left( \boldsymbol{\Theta} - \boldsymbol{\Theta} \right) \left( \boldsymbol{\Theta} - \boldsymbol{\Theta} \right)^{T} \right],\tag{12}$$

where **Θ** are the latent variable parameter collection, **Θ**ˆ denote estimates of the corresponding variables, and **E**(*x*) denotes the expectation operator for the variable *x*.

The Bayesian Information Matrix (BIM), which is commonly employed in timeinvariant statistical models, is often used for CRLB analysis and can be written by

$$\mathbf{J} = -\mathbf{E}\left\{ \bigtriangledown \big[ \big[ \big[ \big[ \big[ \big[ \big[ \Theta \big[ \big] \big] \big] \big] \big] \big] \big] \bigtriangleq \mathbf{E}\left\{ \big[ \big[ \big[ \big[ \Theta \big] \big] \big] \big( \Theta \big) \big] \big[ \big[ \big[ \big[ \Theta \big] \big] \big] \big} \big( \Theta \big) \big] \big) \big} \right\} \right) \right\} \tag{13}$$

where *p*(**Θ**,**Z**) is the joint probability distribution, and **<sup>Θ</sup>** denotes the first-order partial derivative operator concerning the variable vector **Θ**, which can be expressed by

$$
\nabla \clubsuit = \left[ \frac{\partial}{\partial \Theta\_1}, \frac{\partial}{\partial \Theta\_2}, \dots \right]^T \tag{14}
$$

According to the Bayesian CRLB analysis, we have

$$
\Sigma\_{\mathbf{c}} \gtrless \mathbf{J}^{-1}.\tag{15}
$$

The inequality suggests that the matrix on the left is non-negative.

According to posterior Cramer–Rao bounds analysis, we have the recursive sequence **J***<sup>k</sup>* of posterior information submatrices for estimating the state vector *θ<sup>k</sup>* at the *k*th time step [28]

$$\mathbf{J}\_k = \mathbf{D}\_k^{22} - \mathbf{D}\_k^{21} \left(\mathbf{J}\_{k-1} + \mathbf{D}\_k^{11}\right)^{-1} \mathbf{D}\_k^{12} \tag{16}$$

where

$$\mathbf{D}\_{k}^{11} = \mathbf{E} \left[ -\nabla \boldsymbol{\theta}\_{k-1} \, \bigtriangledown \boldsymbol{\theta}\_{k-1}^{T} \, \ln p(\boldsymbol{\theta}\_{k} | \boldsymbol{\theta}\_{k-1}) \right],\tag{17}$$

$$\mathbf{D}\_k^{12} = \mathbf{E} \left[ -\nabla \mathfrak{a}\_{k-1} \, \bigtriangledown \mathfrak{f}\_{\theta\_k}^T \ln p(\theta\_k | \theta\_{k-1}) \right],\tag{18}$$

$$\mathbf{D}\_{k}^{21} = \mathbb{E}\left[ -\left. \nabla \boldsymbol{\theta}\_{k} \right| \boldsymbol{\nabla}\_{\boldsymbol{\theta}\_{k-1}}^{T} \ln p(\boldsymbol{\theta}\_{k}|\boldsymbol{\theta}\_{k-1}) \right] = (\mathbf{D}\_{k}^{12})^{T}.\tag{19}$$

**D**<sup>22</sup> *<sup>k</sup>* can be decomposed into

$$\mathbf{D}\_k^{22} = \mathbf{D}\_k^{22,a} + \mathbf{D}\_k^{22,b},\tag{20}$$

where

$$\mathbf{D}\_{k}^{22,\mathfrak{a}} = \mathbf{E} \left[ -\nabla \mathfrak{g}\_{k} \, \bigcirc\_{}^{T}\_{\mathfrak{G}\_{k}} \ln p(\mathfrak{G}\_{k}|\mathfrak{G}\_{k-1}) \right],\tag{21}$$

$$\mathbf{D}\_{k}^{22,b} = \mathbf{E}\left[ -\left. \nabla \boldsymbol{\theta}\_{k} \, \big|\, \nabla\_{\boldsymbol{\theta}\_{k}}^{T} \ln p(\mathbf{z}\_{k}|\boldsymbol{\theta}\_{k}) \right] \right] \tag{22}$$

The initial information submatrix can be obtained from the prior probability distribution *p*(*θ*0)

$$\mathbf{J}\_0 = \mathbf{E}\left[-\left.\nabla\theta\_0 \,\nabla\right|\_{\theta\_0} \ln p(\theta\_0)\right],\tag{2.3}$$

where *p*(*θ*0) denotes the initial prior distribution of the variable vector *θ*0.

Using the nonlinear discrete-time state-space model in Equations (1) and (10), we have

$$-\ln p(\boldsymbol{\theta}\_{k}|\boldsymbol{\theta}\_{k-1}) = \frac{1}{2}\left\{ \left[\boldsymbol{\theta}\_{k} - f(\boldsymbol{\theta}\_{k-1})\right]^{T} \mathbf{Q}^{-1} [\boldsymbol{\theta}\_{k} - f(\boldsymbol{\theta}\_{k-1})] \right\} + \text{const},\tag{24}$$

$$-\ln p(\mathbf{z}\_k|\boldsymbol{\theta}\_k) = \sum\_{m=1}^{M} \left(\frac{\upsilon\_m}{2} + \frac{1}{2}\right) \ln \left\{ 1 + \frac{\lambda\_m [\mathbf{z}\_{k,m} - h(\boldsymbol{\theta}\_k)\_m]^2}{\upsilon\_m} \right\} + \text{const.} \tag{25}$$

Then,

$$\mathbf{D}\_{k}^{11} = \mathbf{E}\_{\theta\_{k-1}} \left\{ \left[ \nabla \boldsymbol{\theta}\_{k-1} f^{T}(\boldsymbol{\theta}\_{k-1}) \right] \mathbf{Q}^{-1} \left[ \nabla \boldsymbol{\theta}\_{k-1} f^{T}(\boldsymbol{\theta}\_{k-1}) \right]^{T} \right\} \tag{26}$$

$$\mathbf{D}\_{k}^{12} = -\mathbf{E}\_{\theta\_{k-1}} \left\{ \left[ \nabla \boldsymbol{\theta}\_{k-1} \boldsymbol{f}^{T} (\boldsymbol{\theta}\_{k-1}) \right] \right\} \mathbf{Q}^{-1} \tag{27}$$

$$\mathbf{D}\_k^{21} = \left[\mathbf{D}\_k^{12}\right]^T,\tag{28}$$

$$\mathbf{D}\_{k}^{22} = \mathbf{Q}^{-1} + \mathbf{E}\_{\theta\_{k}} \Big[ \bigtriangledown \theta\_{k} \bigtriangledown^{T}\_{\theta\_{k}} \ln p(\mathbf{z}\_{k}|\theta\_{k}) \Big]. \tag{29}$$

It is observed that it is difficult to carry out these integrals, the sequential Monte Carlo method is used to perform the approximations, and we define the variables **Λ**11, **Λ**12, and **Λ**<sup>22</sup> as

$$\mathbf{A}^{11}(\boldsymbol{\theta}\_{k-1}) = \left[ \bigtriangledown \boldsymbol{\theta}\_{k-1} \boldsymbol{f}^T(\boldsymbol{\theta}\_{k-1}) \right] \mathbf{Q}^{-1} \left[ \bigtriangledown \boldsymbol{\theta}\_{k-1} \boldsymbol{f}^T(\boldsymbol{\theta}\_{k-1}) \right]^T,\tag{30}$$

$$\Lambda^{12}(\theta\_{k-1}) = -\bigtriangledown \theta\_{k-1} f^T(\theta\_{k-1})\_\prime \tag{31}$$

$$\mathbf{A}^{ZZ}(\theta\_k) = \bigtriangledown \bigotimes\_{\theta\_k} \bigtriangledown p(\mathbf{z}\_k | \theta\_k). \tag{32}$$

Therefore, we can rewrite Equations (26)–(29) as follows:

$$\begin{split} \mathbf{D}\_{k}^{11} &= \int \mathbf{A}^{11}(\boldsymbol{\theta}\_{k-1}) p(\boldsymbol{\theta}\_{k-1}) d\boldsymbol{\theta}\_{k-1} \\ &= \int p(\mathbf{z}\_{k-1}) \left\{ \int \mathbf{A}^{11}(\boldsymbol{\theta}\_{k-1}) p(\boldsymbol{\theta}\_{k-1}|\mathbf{z}\_{k-1}) d\boldsymbol{\theta}\_{k-1} \right\} d\mathbf{z}\_{k-1} \end{split} \tag{33}$$

$$\begin{split} \mathbf{D}\_{k}^{12} &= \int \mathbf{A}^{12}(\theta\_{k-1}) p(\theta\_{k-1}) d\theta\_{k-1} \\ &= \int p(\mathbf{z}\_{k-1}) \left\{ \int \mathbf{A}^{12}(\theta\_{k-1}) p(\theta\_{k-1}|\mathbf{z}\_{k-1}) d\theta\_{k-1} \right\} d\mathbf{z}\_{k-1}, \end{split} \tag{34}$$

$$\begin{split} \mathbf{D}\_{k}^{22} &= \mathbf{Q}^{-1} + \int \mathbf{A}^{22}(\theta\_{k}) p(\theta\_{k}) d\theta\_{k} \\ &= \mathbf{Q}^{-1} + \int p(\mathbf{z}\_{k}) \left\{ \int \mathbf{A}^{22}(\theta\_{k}) p(\theta\_{k}|\mathbf{z}\_{k-1}) d\theta\_{k} \right\}. \end{split} \tag{35}$$

$$=\mathbf{Q}^{-1} + \int p(\mathbf{z}\_k) \left\{ \int \boldsymbol{\Lambda}^{22}(\boldsymbol{\theta}\_k) p(\boldsymbol{\theta}\_k | \mathbf{z}\_k) d\boldsymbol{\theta}\_k \right\} d\mathbf{z}\_k. \tag{2.1}$$
  $\text{L'can be amnovximated by the mean of samples after having a series of } \mathbf{Q}$ 

These integrals can be approximated by the mean of samples after having a series of samples corresponding to the density. We perform *T* Monte Carlo samples with the time steps of *K* for the calculation of the expectation instead of the integral operator over **z***<sup>k</sup>* and obtain the measurement sample matrix **Z**1:*T*,1:*K*. **z***t*,*k*−<sup>1</sup> denotes the *t*th measurement sample at the (*k* − 1)th time step. Consequently, Equations (33)–(35) can be calculated by

$$\mathbf{D}\_{k}^{11} \approx \frac{1}{T} \sum\_{t=1}^{T} \left\{ \int \Lambda^{11}(\boldsymbol{\theta}\_{k-1}) p(\boldsymbol{\theta}\_{k-1}|\mathbf{z}\_{t,k-1}) d\boldsymbol{\theta}\_{k-1} \right\} \tag{36}$$

$$\mathbf{D}\_{k}^{12} \approx \frac{1}{T} \sum\_{t=1}^{T} \left\{ \int \Lambda^{12}(\boldsymbol{\theta}\_{k-1}) p(\boldsymbol{\theta}\_{k-1}|\mathbf{z}\_{t,k-1}) d\boldsymbol{\theta}\_{k-1} \right\} \tag{37}$$

$$\mathbf{D}\_{k}^{22} \approx \mathbf{Q}^{-1} + \frac{1}{T} \sum\_{t=1}^{T} \left\{ \int \Lambda^{22} (\theta\_k, \mathbf{z}\_k) p(\theta\_k | \mathbf{z}\_{t,k}) d\theta\_k \right\} \tag{38}$$

We further perform Monte Carlo samples over the variable *θ<sup>k</sup>* instead of an integral operator over *θ<sup>k</sup>* at the *t*th sample, and formulas inside the braces in Equations (36)–(38) can be expressed as

$$\mathbf{D}\_{t,k}^{11} \approx \frac{1}{L} \sum\_{l=1}^{L} \boldsymbol{\Lambda}^{11} (\boldsymbol{\theta}\_{t,k-1}^{(l)}) \tag{39}$$

$$\mathbf{D}\_{t,k}^{12} = \left(\mathbf{D}\_{t,k}^{21}\right)^T \cong \frac{1}{L} \sum\_{l=1}^{L} \Lambda^{12}(\boldsymbol{\theta}\_{t,k-1}^{(l)}) \tag{40}$$

$$\mathbf{D}\_{t,k}^{22} \approx \frac{1}{L} \sum\_{l=1}^{L} \Lambda^{22} (\boldsymbol{\theta}\_{t,k'}^{(l)} \mathbf{z}\_{t,k}) \tag{41}$$

where *θ* (*l*) *<sup>t</sup>*,*k*−<sup>1</sup> and *<sup>θ</sup>* (*l*) *<sup>t</sup>*,*<sup>k</sup>* are, respectively, samples drawn from the posterior distributions of *<sup>p</sup>*(*θt*,*k*−1|*zt*,*k*−1) and *<sup>p</sup>*(*θt*,*k*|*zt*,*k*), and *<sup>L</sup>* is the number of samples. According to the approximate inference method, the Markov Chain Monte Carlo (MCMC) method is utilized to estimate *<sup>θ</sup>t*,*<sup>k</sup>* with the state transition probability *<sup>p</sup>*(*θk*|*θk*−1) in Equation (24) and the likelihood function *p*(**z***k*|*θk*) in Equation (25). The following Equations (33)–(35) and (39)– (41), can be rewritten by

$$\mathbf{D}\_{k}^{11} \approx \frac{1}{T} \sum\_{t=1}^{T} \left\{ \frac{1}{L} \sum\_{l=1}^{L} \Lambda^{11} (\boldsymbol{\theta}\_{t,k-1}^{(l)}) \right\} \tag{42}$$

$$\mathbf{D}\_k^{12} \approx \left(\mathbf{D}\_k^{21}\right)^T \cong \frac{1}{T} \sum\_{t=1}^T \left\{ \frac{1}{L} \sum\_{l=1}^L \Lambda^{12} (\boldsymbol{\theta}\_{t,k-1}^{(l)}) \right\} \tag{43}$$

$$\mathbf{D}\_k^{22} \approx \mathbf{Q}^{-1} + \frac{1}{T} \sum\_{t=1}^T \left\{ \frac{1}{L} \sum\_{l=1}^L \boldsymbol{\Lambda}^{22} (\boldsymbol{\theta}\_{t,k'}^{(l)} \mathbf{z}\_{t,k}) \right\} \tag{44}$$

When the number of samples *L* is large enough, the estimated results would converge to true values. Finally, we arrive at the resulting Monte-Carlo-based posterior Cramer–Rao bound (MC-PCRB).

#### **3. Simulations and Experiments**

In this simulation, we take target tracking in the automotive radar system as an example. The state process model with the constant velocity, which is a typically linear equation, can be expressed by

$$\boldsymbol{\theta}\_{k} = f(\boldsymbol{\theta}\_{k-1}) + \mathbf{u}\_{k} = \mathbf{A}\boldsymbol{\theta}\_{k-1} + \mathbf{u}\_{k} \tag{45}$$

where *θ<sup>k</sup>* = [*xk*, *yk*, *x*˙*k*, *y*˙*k*], containing the two-dimensional location (*xk*, *yk*) and the corresponding velocities (*x*˙*k*, *y*˙*k*) of the vehicle. The state transition matrix **A** can be written by

$$\mathbf{A} = \begin{bmatrix} 1 & 0 & \Delta T & 0 \\ 0 & 1 & 0 & \Delta T \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{46}$$

where Δ*T* denotes the sampling interval of the radar. The nonlinear measurement model in the automotive radar system can be denoted by,

$$
\begin{bmatrix} R\_k \\ \phi\_k \\ V\_k \end{bmatrix} = h \begin{pmatrix} \begin{bmatrix} \mathbf{x}\_k \\ \mathbf{y}\_k \\ \dot{\mathbf{x}}\_k \\ \dot{\mathbf{y}}\_k \end{bmatrix} \\
\end{bmatrix} = \begin{bmatrix} \sqrt{\mathbf{x}\_k^2 + \mathbf{y}\_k^2} \\ \arctan\left(\frac{\mathbf{x}\_k}{\mathbf{y}\_k}\right) \\ \frac{\mathbf{x}\_k \mathbf{x}\_k + \mathbf{y}\_k \mathbf{y}\_k}{\sqrt{\mathbf{x}\_k^2 + \mathbf{y}\_k^2}} \end{bmatrix}. \tag{47}
$$

where *Rk* denotes the measured range of the target at the time step *k*, *φ<sup>k</sup>* represents the measured spatial information, and *Vk* is its measured radial velocity.

#### *3.1. Performance Testing on Simulated Data*

In this simulation, a spiral trajectory is considered, and the spiral trajectory can be generated by a radius reduction circular motion. Assume that the radar continuously acquires the observed measurement vector, and the measurements containing outliers are generated according to the mixture model

$$\mathbf{z}\_{k}|\boldsymbol{\theta}\_{k}, \boldsymbol{c}\_{k,m} \sim \begin{cases} \mathbf{U}(\mathbf{z}\_{k}|a\_{m}, b\_{m}), & \text{if } c\_{k,m} = 0\\ \mathcal{N}(\mathbf{z}\_{k}|h(\boldsymbol{\theta}\_{k}), \mathbf{R}), & \text{if } c\_{k,m} = 1 \end{cases} \tag{48}$$

where *ck*,*<sup>m</sup>* = 0 indicates a normal measurement from the *m*th dimension at the *k*th time step clutter and *U*(·) is a uniform distribution, while *ck*,*<sup>m</sup>* = 1 indicates the measurement vector is a normal one, and noises follow Gaussian distribution with the measurement noise covariance matrix **R** = diag([4, 0.01, 1]). The processed noise covariance matrix is assumed to be **Q** = diag([0.1, 0.1, 0.01, 0.01]). Multidimensional outliers caused by interferences or clutters are assumed to be drawn from the uniform distribution with the following parameters: [*a*1, *b*1]=[10, 90] in the range entry in the measurement vector, [*a*2, *b*2]=[−*π*/2, *π*/2] in the azimuth entry, and [*a*3, *b*3]=[−10, 10] in the velocity entry. The probability of outliers is assumed to be constant over time, that is, *ck*,*<sup>m</sup>* is a Bernoulli random variable with the probability of *p*(*ck*,*<sup>m</sup>* = 0) = 0.2. In the particle filter method, the number of particles *S* is 500, and hyperparameters **v** are set to be [10, 10, 10] *<sup>T</sup>*. The total step number *K* is 1000 with the sampling interval Δ*T* of 0.01 s. All the results are implemented on the PC with CPU (Intel Core i7-10700K) using Matlab R2020a.

For the performance evaluation of estimated target location, the root of mean square error (RMSE) and normalized mean square error (NMSE) are defined as performance indexes

$$\text{RMSE}\_{k} = \sqrt{\left(\langle \mathbf{x}\_{k} \rangle - \mathbf{x}\_{k}^{\text{TRUE}}\right)^{2} + \left(\langle y\_{k} \rangle - y\_{k}^{\text{TRUE}}\right)^{2}} \tag{49}$$

$$\text{NMSE} = \frac{1}{K} \sum\_{k=1}^{K} \frac{\text{RMSE}\_k}{\sqrt{\left(x\_k^{\text{TRUE}}\right)^2 + \left(y\_k^{\text{TRUE}}\right)^2}} \tag{50}$$

where *θk* = [*xk*,*yk*,*x*˙*k*,*y*˙*k*] *<sup>T</sup>* denotes the estimated state vector at the step time *k*. The performance comparisons of tracking results with the conventional particle filter (PF) and other robust solutions, including the weighted robust PF (WR-PF) and the outlier robust extended KF (OR-EKF), are provided in Figure 3. It is observed that the estimated trajectory is the closest to the truth one compared with the WR-PF and OR-EKF methods, and thus the proposed method has the smallest NMSE. Note that the results in the PF method are not provided due to quite a sensitivity to outliers, yielding significantly poor performances.

**Figure 3.** Performance comparisons on simulated data. (**a**) Estimated target locations. (**b**) Estimated RMSE.

The performances versus outlier ratios with *p*(*ck*,*<sup>m</sup>* = 0) are also examined. The outlier ratio varies from 0% to 50%, and all other parameters are kept in this simulation. It is reasonable that the NMSEs generally decreases with decreases in the outlier ratios. It is observed that the proposed method has the lowest NMSE across all ratios compared with other methods, as shown in Figure 4.

**Figure 4.** Estimated NMSE versus outlier ratios.

In addition, we also introduce average Mahalanobis distance as the performance index to evaluate the performance of the entire estimated state-space vector and defined it as

$$\mathbf{D}\_{\mathbf{M}} = \frac{1}{K} \sum\_{k=1}^{K} \sqrt{\left( \langle \theta\_k \rangle - \theta\_k^{\mathrm{TRUE}} \right) \mathbf{Q}^{-1} \left( \langle \theta\_k \rangle - \theta\_k^{\mathrm{TRUE}} \right)} \tag{51}$$

where *θ<sup>k</sup>* = [*xk*, *yk*, *x*˙*k*, *y*˙*k*] is the true state vector containing the two-dimensional location (*xk*, *yk*) and the corresponding velocities (*x*˙*k*, *y*˙*k*) of the vehicle. *θk* denotes the estimated vector of *θk*, and **Q** is the process noise covariance matrix provided in the simulations.

The performance results versus outlier ratios are also examined in Figure 5, the proposed method has the lowest Mahalanobis distance across all ratios, and thus has better performance than those in the WR-PF and OR-EKF methods. Note that the results in the PF method are not provided due to quite a sensitivity to outliers, yielding significantly poor performances.

**Figure 5.** Estimated Mahalanobis distance versus outlier ratios.

In the following simulations, diverse types of outliers in the multidimensional measurements are considered to examine the performances of the proposed method. First, we consider only range outliers in the multidimensional measurements, that is to say, the azimuth and velocity entries are abnormal values with *p*(*ck*,2) = *p*(*ck*,3) = 0. It is found that both PF-based methods have better performances than that in the KF-based method, shown in Figure 6. The main reason may be the nonlinear approximation in the measurement function in the EKF method. These two PF-based outlier-robust methods can effectively address the range outlier issue and acquire comparable NMSEs, even though the NMSEs in the WR-PF method are a little bit smaller than those in the proposed method in the target location estimation, shown in Figure 6a. The possible reason could be that the entry of range is much larger than the other two entries in the multidimensional measurements, and thus the scalar weight in the WR-PF method is learned and estimated to accurately

model the measurement noise variances, particularly in the noise variance of the range entry in the multidimensional measurements, such that the lower NMSEs in the target locations can be obtained. However, to evaluate the performance of the entire estimated state-space vector in the filtering methods, the proposed method has the lowest average Mahalanobis distance across outlier ratios, shown in Figure 6b. Therefore, it suggests that the proposed method has better estimation performances in the entire state-space vector estimation by treating the weights associated with each entry in the measurement vector data probabilistically.

**Figure 6.** Estimated errors versus only range outliers in multidimensional measurements. (**a**) NMSE comparisons. (**b**) Mahalanobis distance comparisons.

Then, we consider only azimuth outliers in the multidimensional measurements, and the other two entries are normal values. The proposed method has the lowest estimated NMSE in both average Mahalanobis distance and NMSE compared with those in the OR-EKF and WR-PF methods, shown in Figure 7. Finally, Figure 8 shows the performance comparisons between the proposed method and the other two robust solutions in the only velocity outliers case, and it is observed that the proposed method has clear performance improvements and has the lowest NMSE and average Mahalanobis distance across normal measurement ratios.

**Figure 7.** Estimated errors versus only azimuth outliers in multidimensional measurements. (**a**) NMSE comparisons. (**b**) Mahalanobis distance comparisons.

**Figure 8.** Estimated errors versus only radial velocity outliers in multidimensional measurements. (**a**) NMSE comparisons. (**b**) Mahalanobis distance comparisons.

Table 1 shows the average executive comparisons of these robust solutions. It is observed that the executive time of the proposed method is comparable to these PF-based methods, whereas the proposed method has the best performance among them. It is reasonable that the OR-EKF method has the smallest executive time among all these methods owing to its analytical formulation, whereas it generally has slightly poorer performance than the PF-based robust solutions such as the WR-PF and RR-PF methods.

**Table 1.** Executive time comparisons.


#### *3.2. Cramer–Rao Bound Analysis of the Proposed Algorithm*

To further illustrate the estimation performance of the algorithm, we provide the Cramer–Rao Bound for the nonlinear system to serve as the reference for the estimation error. We present the universal MC-PCRB method above for discrete-time state-space models to solve the iterative Cramer–Rao bound at each step in Equation (16). Figure 9 shows the comparison results among these methods. The outlier ratios are assigned to be *p*(*ck*,1) = *p*(*ck*,2) = *p*(*ck*,3) = 0.2. It is observed that the RMSEs in the outlier robust methods generally decrease with the increases in the time step, except in the conventional PF method. The proposed method is generally close to the PCRLB, and thus has the lowest RMSE. The proposed method has higher RMSEs than the WR-PF method in the early stage because of more parameters to learn. In theory, the proposed method would approach the PCRLB, as long as the number of particles approaches infinity while tracking under a stable state.

**Figure 9.** Bounds comparisons of estimated errors.

#### *3.3. Experimental Results Based on Real Data*

In this subsection, real data collected by the automotive millimeter wave (MMW) radar with linear modulated continuous waveform (LFMCW) are used to verify the algorithm's effectiveness. In the automotive radar application, pedestrian tracking based on the MMW radar is often an issue of great importance. Therefore, we consider pedestrian tracking in these experiments. A 77 GHz millimeter wave radar, whose parameters are shown in Table 2, is used to carry out experiments. An LFMCW waveform with a bandwidth of 4 GHz is used for the high range resolution of 0.0375 m, and the sampling rate of 6.25 MHz is used to acquire raw data. Two-dimensional (2D) Fourier Transform (FT) is firstly implemented on the raw data, and thus the constant false alarm ratio (CFAR) technique is used to acquire the estimates of the range and velocity of the walking pedestrian. Finally, the digital beamforming technique is utilized to obtain the azimuth of the pedestrian. To this end, the measurement vector is acquired for further target tracking. During the data collection of the pedestrian, several artificial blocks are place along the tracking path, and then some outlier measurements occur in the acquired measurement vector.



In the first experiment trial, the pedestrian walked from the location (1.5, 2.5) to the location (−1.5, 2.1) along a straight line. There is an artificial blocker from the location (−0.3, 1.5) to the location (0.3, 1.5). The sketch map is shown in Figure 10. The measurement vector data have some outliers from the 30th to the 37th time step because of crossing the blocker. The acquired measurement vector data are given in Figure 11. It is observed that the measurement vectors involve several obvious outliers at several time steps due to the artificial occlusion or clutters reflected from the complex environment. The tracking results and the corresponding RMSEs are, respectively, shown in Figure 12a,b. It demonstrates that the estimated trajectory in the proposed method is much closer to the true one, and thus it has the lowest RMSE compared with other methods. It is found that there exist a few fluctuations in the interval from the 20th time step and the 35th time step; a possible reason may result from the measurement precision issue of pedestrian walking trajectory in the real experiment, and the true trajectory of the walking pedestrian may be around the preassigned straight line, with some inevitable fluctuations when the pedestrian walked along the straight line in the natural walking state. Although a few fluctuations appear, the RMSEs of the proposed method reach a pretty small level, and they are almost the lowest across time steps compared with other robust solutions.

**Figure 10.** Experiment environment sketch map.

**Figure 11.** Raw data plots in three-dimensional entries in the measurement vector.

**Figure 12.** Performance comparisons of tracking results in the real data. (**a**) Tracking results based on real data. (**b**) RMSE comparisons.

In the second experiment trial, the slightly complex trajectory, which involves curvilinear motion, turning around, and sharp turning, is considered to verify the effectiveness of the proposed method. The acquired measurement vector data are provided in Figure 13, and it is observed that the measurement vectors have obvious outliers at a few time steps due to the clutters reflected from the complex electromagnetic environment. The tracking results and the corresponding RMSEs are, respectively, shown in Figure 14a,b. They demonstrate that the estimated trajectory in the proposed method is much closer to the true one and almost has the lowest RMSE compared with other robust methods.

**Figure 13.** Raw data plots in three-dimensional entries in the measurement vector.

**Figure 14.** Performance comparisons of tracking results in the real data. (**a**) Tracking results based on real data. (**b**) RMSE comparisons.

#### **4. Conclusions**

In this paper, we consider the nonlinear filtering problem in the presence of contaminated measurements. A novel reweighted robust particle filter method was proposed for accurate target tracking in the automotive radar system. We generalized nonlinear filtering to the more comprehensive case of multidimensional outlier measurements and enhanced state vector estimates by probabilistically learning a weight vector for each entry in the multidimension measurement vector. The particle filter technique was adopted to make an approximate posterior in this nonlinear system. We also illustrated the reason that the proposed method is more robust to outliers. Finally, the Cramer–Rao lower bound of the proposed method was also provided. Both simulations and experiments verified the superiority of the proposed method over other robust solutions.

**Author Contributions:** Conceptualization, Q.W.; methodology, Q.W., L.C., Y.L. and S.Y.; softwarep, L.C.; validation,; formal analysis, Q.W.; investigation, Q.W., Z.W. and H.L.; writing—original draft preparation, Q.W., L.C. and Y.L.; writing—Q.W. and L.C.; supervision, S.Y. and H.L.; funding acquisition, Q.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work was supported in part in part by Key Laboratory of Underwater Acoustic Countermeasure Technology under Grant No. 2022JCJQLB03305, National Natural Science Foundation under Grants No. 91938203, in part by Fundamental Research Funds for the Central University under Grant No. 2242022k30016, in part by National Defense Basis Scientific Research program of China under Grant No. JCKY2019110C143, and in part by Science and Technology on Near-Surface Detection Laboratory Pre-research Fund under Grant No. 6142414200505.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

