**Motion State Estimation of Target Vehicle under Unknown Time-Varying Noises Based on Improved Square-Root Cubature Kalman Filter**

#### **Shiping Song and Jian Wu \***

State Key Laboratory of Automotive Simulation and Control, Jilin University, Changchun 130022, China; songsp17@mails.jlu.edu.cn

**\*** Correspondence: wujian@jlu.edu.cn

Received: 17 March 2020; Accepted: 2 May 2020; Published: 4 May 2020

**Abstract:** In the advanced driver assistance system (ADAS), millimeter-wave radar is an important sensor to estimate the motion state of the target-vehicle. In this paper, the estimation of target-vehicle motion state includes two parts: the tracking of the target-vehicle and the identification of the target-vehicle motion state. In the unknown time-varying noise, non-linear target-vehicle tracking faces the problem of low precision. Based on the square-root cubature Kalman filter (SRCKF), the Sage–Husa noise statistic estimator and the fading memory exponential weighting method are combined to derive a time-varying noise statistic estimator for non-linear systems. A method of classifying the motion state of the target vehicle based on the time window is proposed by analyzing the transfer mechanism of the motion state of the target vehicle. The results of the vehicle test show that: (1) Compared with the Sage–Husa extended Kalman filtering (SH-EKF) and SRCKF algorithms, the maximum increase in filtering accuracy of longitudinal distance using the improved square-root cubature Kalman filter (ISRCKF) algorithm is 45.53% and 59.15%, respectively, and the maximum increase in filtering the accuracy of longitudinal speed using the ISRCKF algorithm is 23.53% and 29.09%, respectively. (2) The classification and recognition results of the target-vehicle motion state are consistent with the target-vehicle motion state.

**Keywords:** millimeter-wave radar; square-root cubature Kalman filter; Sage-Husa algorithm; target tracking; stationary and moving object classification

#### **1. Introduction**

Millimeter-wave radar is an important sensor that constitutes an advanced driver assistance system (ADAS). The estimation of the moving state of the target vehicle based on the on-board millimeter-wave radar is essential for predicting the future trajectory of the target vehicle and determining the degree of danger of the target vehicle to the ego vehicle. In this paper, the motion state estimation of the target vehicle includes target-vehicle tracking and target-vehicle motion state classification.

The motion state information of the target vehicle (relative radial distance, azimuth, and relative radial rate) measured by millimeter-wave radar is obtained from the polar coordinate system. However, in the process of target tracking, the target motion model is usually established in the Cartesian coordinate system. As can be seen from the radar target-tracking process, the state equation is linear and the measurement equation is non-linear. Since the measurement equation is non-linear, the target-tracking system based on the millimeter-wave radar is a non-linear system.

Extended Kalman filter (EKF) [1,2], unscented Kalman filter (UKF) [3,4], particle filter (PF) [5] and cubature Kalman filter (CKF) [6] are common non-linear filtering state estimation algorithms.

The basic idea of EKF is: the non-linear system is linearized by Taylor series expansion, and then Kalman filtering is performed. Inaccurate modeling of system noise and changes in model parameters

due to environmental factors will cause the decrease of the EKF estimation accuracy, and even the filter divergence will be caused. To ensure the accuracy and stability of EKF under unknown and time-varying conditions, many scholars carry out research on the adaptive extended Kalman filter algorithm, adaptive extended Kalman filter (AEKF) algorithm, such as Sage–Husa's maximum a posteriori estimation [7], fictitious noise compensating [8], dynamic bias decoupling estimation [9], etc. To solve the filtering divergence problem caused by system modeling errors, some scholars proposed an EKF with a suboptimal fading factor [10]. Meanwhile, to improve the estimation accuracy of the EKF algorithm, some scholars proposed an AEKF algorithm based on a neural network [11]. The adaptive learning characteristics of neural networks are used to identify the non-linear system model online, to overcome the influence of unmodeled dynamic characteristics of the filter. However, due to the following shortcomings of the EKF algorithm, its development in engineering practice is limited:

(1) The higher-order terms are truncated (with truncation error) and only the first-order terms retained during the Taylor series are expanded, and the accuracy of EKF is the first-order Taylor series;

(2) In many engineering applications, the Jacobian matrix of the measurement equation is difficult to solve.

The basic idea of UKF is: "It is easier to approximate the probability density distribution of non-linear functions than the approximation of non-linear functions" [12]. UKF uses the unscented transformation to approximate the posterior distribution of the state of the non-linear system. The most important part of the UKF algorithm is the sampling strategy. Different sampling strategies differ in the number, location, and corresponding weights of the extracted Sigma points [13]. Compared with EKF, UKF has the following two advantages:

(1) The accuracy of UKF can reach at least two orders under the condition of EKF, and UKF algorithm takes the same order of magnitude;

(2) In UKF, it is not necessary to calculate the Jacobian matrix of the measurement equation.

The above two advantages of UKF expand the application range of the EKF algorithm. However, certain sigma point weights ω < 0 will cause the covariance matrix to non-positive definite condition when the dimension is too high (N ≥ 4). This situation will lead to the following two effects: firstly, the filter value is not stable or even divergent; secondly, the dimension disaster problem will occur [14]. Therefore, some scholars through theoretical analysis and experiments have proved that UKF has high accuracy for low-dimensional (N ≤ 2) non-linear systems [15].

The basic idea of PF is to approximate the posterior probability density function of the system state by random particles. PF uses the particle mean value instead of the integral operation to obtain the minimum variance of state. With the increase in the number of particles, the probability density function of particles gradually approximates the probability density function of the real state. However, PF has the following two shortcomings, which restrict the development of PF [16]:

(1) the particle degradation problem;

(2) It is difficult to realize online estimation due to the large amount of computation.

To better solve the problem of poor filtering performance and even divergence in high-dimensional non-linear filtering estimation, Arasaratnam and Haykin proposed a third-degree spherical-radial rule CKF [17]. After CKF was proposed, it was widely used in target tracking [18] and navigation [19]. Compared with UKF, CKF has the following advantages:

(1) The UKF algorithm selects 2n+1 sampling points with different weights, while the CKF algorithm selects 2n sampling points with the same weight. The CKF algorithm has fewer sampling points than UKF, so the CKF algorithm takes less time than UKF.

(2) Since the weight coefficients of the sampling points of the CKF algorithm are all positive, the robustness of the CKF algorithm is high when the dimension of the observed variable is too high (N ≥ 4).

In conclusion, compared with EKF and UKF, CKF has higher estimation accuracy. During the operation of the standard CKF algorithm, the following two conditions should be met: (1) symmetry, (2) positive qualitative. However, in practical engineering, these two conditions are sometimes difficult

to meet. Therefore, scholars proposed the SRCKF algorithm based on the CKF algorithm [20]. SRCKF has the following two advantages. On the one hand, SRCKF avoids computing the square root of a matrix by directly calculating the square root of the predicted value of error covariance and the estimated value of error covariance. On the other hand, in the SRCKF algorithm, the symmetry and positive qualitative value of the error covariance matrix can always be guaranteed.

During the derivation of the CKF algorithm, it is generally assumed that the statistical characteristics of system noise and measurement noise are known [21]. However, in practice, the statistical characteristics of noise are often unknown and time-varying. The Sage–Husa estimator is often used to estimate the statistical characteristics of noise because of its simplicity and good real-time performance [22]. However, the conventional Sage–Husa estimator is suitable for estimating the statistical properties of constant coefficient noise in linear systems [23]. Based on the conventional Sage–Husa algorithm, an adaptive noise statistical estimator for non-linear systems is derived by using the cubature rule.

For time-varying noise statistics, the real-time updated data play a leading role, while the old data play a small role compared with the new data. Therefore, we should gradually reduce the weight of old data and increase the weight of new data. The exponential weighted attenuation method for fading memory is introduced to estimate time-varying noise. The exponential weighted attenuation method has the characteristic of remembering the past historical data, but the weighted coefficient of the old data is small [24].

The current international standard "ISO/DIS15622 Intelligent transportation systems-adaptive cruise control systems-performance requirements and test procedures" clearly states that adaptive cruise control (ACC) may ignore stationary targets or do not respond to stationary targets. At the same time, for full-speed ACC and autonomous emergency braking (AEB) systems, it is necessary to accurately identify the target-vehicle as a stationary target-vehicle or a moving target-vehicle.

The recognition of the target motion state has the following two functions for the ADAS. On the one hand, it can predict the future trajectory of the target-vehicle; On the other hand, it can determine the degree of danger of the target-vehicle to the ego-vehicle. Therefore, it is essential to study the classification of the target-vehicle motion state. In the literature [25], targets detected by radar are divided into high-altitude targets, stationary targets, moving targets, and road targets, but the basis for target classification is not discussed in detail. Therefore, a method of classifying the motion state of the target-vehicle based on a time window is proposed by analyzing the transfer mechanism of the motion state of the target-vehicle.

The motivation of writing the paper is as follows: (1) For the on-board millimeter-wave radar in the unknown and time-varying noise environment, the accuracy of a high-dimensional non-linear target tracking process is low. The ISRCKF algorithm based on SRCKF is proposed to accurately estimate the unknown and time-varying noise statistics. (2) To accurately predict the future trajectory of the target vehicle and determine the danger degree of the target vehicle to the ego vehicle. We present a classification method for moving objects and stationary objects based on the mechanism of moving state transfer in a time window. The vehicle test results show: (1) The filter accuracy of the ISRCKF algorithm is higher than that of SRCKF and SH-EKF. (2) The classification and recognition results of the target-vehicle's motion state are consistent with the target-vehicle's motion state.

The rest of the paper is organized as follows. Section 2, based on the Cartesian coordinate system of millimeter-wave radar, the target-vehicle motion state model is established; In Section 3, based on the SRCKF, an adaptive square-root cubature Kalman filter (ASRCKF) is derived. In Section 4, based on the analysis of the motion state and transfer principle of the target, a classification method of moving target and stationary target based on the motion state transfer mechanism in a time window is proposed. In Section 5, the algorithm is validated and its results are analyzed by establishing a vehicle test platform. Section 6 presents the conclusions.

#### **2. Motion Model of Target-Vehicle**

#### *2.1. Coordinate System*

It can be learned from the dynamics that the description of the same target motion state varies in different reference coordinate systems. Therefore, it is meaningful to clarify the coordinate system that describes the target motion.

The target measurement information (relative radial distance, azimuth, and relative radial speed) measured by the millimeter-wave radar is obtained from the millimeter-wave radar polar coordinate system. This paper studies the target tracking algorithm based on the millimeter-wave radar Cartesian coordinate system, to verify the performance and accuracy of the proposed algorithm in the target tracking process.

The three coordinate systems herein are respectively, the geodetic coordinate system *xoyozo*, on the horizontal ground, the vehicle motion coordinate system *xv yvzv* with its origin coinciding with the center of gravity vehicle, the millimeter-wave radar coordinate system *xroyr*, as shown in the Figure 1. The millimeter-wave radar is fixedly mounted on the front of the vehicle and the radar beam is aligned with the longitudinal axis of the vehicle. Therefore, the millimeter-wave radar coordinate system *xroyr* is parallel to the vehicle motion coordinate system *xv yvzv*. The radar *xr* axis direction is identical with the *xv* direction in the vehicle motion coordinate system. The radar *yr* axis direction is identical with the *yv* direction in the vehicle motion coordinate system.

**Figure 1.** Coordinate system.

#### *2.2. Equation of Motion*

The forward millimeter-wave radar is installed in the middle of the front bumper and fixed to the vehicle body. As the vehicle travels, the millimeter-wave radar detects the target in front of the vehicle, mainly referring to the target-vehicles. These target-vehicles have no vertical motion or small moving speed in the vertical direction, so only the movement of the target-vehicles in the XY plane needs to be considered. Since the target-vehicles motion state has the characteristics of small mobility and low speed, a constant velocity model is established based on the millimeter-wave radar Cartesian coordinate system to describe the motion state of the front target-vehicles.

The model of the constant velocity of the target-vehicle can be obtained from the millimeter-wave radar Cartesian coordinate system:

$$\begin{cases} \mathbf{x}(k+1) = \mathbf{x}(k) + \dot{\mathbf{x}}(k) \ast T + \frac{1}{2} \mathbf{w}\_{\mathbf{x}}(k) \ast T^2 \\ \dot{\mathbf{x}}(k+1) = \dot{\mathbf{x}}(k) + \mathbf{w}\_{\mathbf{x}}(k) \ast T \\ \mathbf{y}(k+1) = \mathbf{y}\_{(k)} + \dot{\mathbf{y}}(k) \ast T + \frac{1}{2} \mathbf{w}\_{\mathbf{y}}(k) \ast T^2 \\ \dot{\mathbf{y}}(k+1) = \dot{\mathbf{y}}(k) + \mathbf{w}\_{\mathbf{y}}(k) \ast T \end{cases} \tag{1}$$

where (*x*(*k* + 1), *y*(*k* + 1)) represents the longitudinal distance and the lateral distance from the forward target-vehicles in the millimeter-wave radar Cartesian coordinate system at *k* + 1, respectively. . *x*(*k* + 1), . *y*(*k* + 1) is the longitudinal velocity and the lateral velocity of the target-vehicles relative to the millimeter-wave radar Cartesian coordinate system movement at *k* + 1. *T* is the sampling time.

The model of the target motion state is as follows:

$$X(k+1) = A \ast X(k) + B \ast w(k)\tag{2}$$

The target motion state vector *X* = [*x*, . *x*, *y*, . *y*] *<sup>T</sup>* here is designed to describe the motion state of the target-vehicle in the millimeter-wave radar Cartesian coordinate system; where, *A* is the state transition matrix, *B* noise-driven matrix, and *w*(*k*) the measurement noise at k;


The target motion state information of millimeter-wave radar in a polar coordinate system is transformed to the front target vehicle in the millimeter-wave radar Cartesian coordinate system by coordinate transformation, the conversion formula is:

$$\begin{cases} \mathbf{x} = \mathbf{R} \ast \cos \theta \\ \mathbf{y} = \mathbf{R} \ast \sin \theta \end{cases} \tag{3}$$

⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

where (*x*, *y*) is the position of the target-vehicle in the millimeter-wave radar Cartesian coordinate system θ is the azimuth.

The measurement equation of the target-vehicle in a millimeter-wave radar Cartesian coordinate system:

$$Z(k+1) = \begin{cases} \mathcal{R}(k+1) = \sqrt{\mathbf{x}^2(k+1) + \mathbf{y}^2(k+1)} + \nu\_1(k+1) \\\mathcal{V}(k+1) = \sqrt{\dot{\mathbf{x}}^2(k+1) + \dot{\mathbf{y}}^2(k+1)} + \nu\_2(k+1) \end{cases} \tag{4}$$

where *vi*(*k*), *i* = 1, 2 are the observation noises at k.

Equations (2) and (4) are the state equation and the measurement equation, respectively. As can be seen from the system model, the state equation is linear and the measurement equation is non-linear. Since the measurement equation is non-linear, the target tracking based on the millimeter-wave radar Cartesian coordinate system in a non-linear system.

#### *2.3. Parameters of the Millimeter-Wave Radar*

There are two main frequency bands for automobile millimeter-wave radar: 77 GHz and 24 GHz. The 24 GHz millimeter-wave radar is usually installed on the side of the vehicle, the detection range is small and mainly used for blind spot detection (BSD), lane change assistance (LCA)**.** The 77 GHz millimeter-wave radar has a large detection range and is usually installed in front of the vehicle. It is mainly used for forward collision warning (FCW) and AEB. In order to verify the accuracy of the algorithm's estimation of the motion state of the target vehicle ahead, a 77 GHz millimeter-wave radar is used.

Table 1 shows the specific technical specifications of the on-board millimeter-wave radar provided by the supplier. Its main technical parameters such as detection range, measurement accuracy, and resolution are usually provided in its specification. A 77 GHz millimeter-wave radar is used herein provided by a supplier. The radar data rate is 40 ms. There are long-distance radar and medium-distance radar working states, which are subject to the driving speed of the ego-vehicle. When the travel speed of ego-vehicle is greater than 80 km/h and the millimeter-wave radar is in the long-distance radar working mode, it achieves the farthest detection distance of 120 m, and the angle detection range of around 30◦. When the travel speed of the ego-vehicle is less than 80 km/h and the millimeter-wave radar is in the medium-distance radar working mode, it reaches the farthest detection distance of 100 m, and the angle detection range of around 50◦.


**Table 1.** Main technical parameters of the millimeter-wave radar.

Remark: -1 The horizontal angle range is negative when the target is on the left side of the radar and positive when on the right side. -2 The relative speed is positive when the target is far from the radar, and negative when close to the radar.

#### **3. Adaptive Square-Root Cubature Kalman Filter**

This section is based on SRCKF. The Sage–Husa noise statistic estimator is extended from the linear constant noise statistic estimator to the non-linear time-varying noise statistic estimator.

#### *3.1. Cubature Rule*

CKF is derived from Bayesian filter theory in the Gaussian field. Under the framework of Bayesian filter theory in Gaussian domain, the non-linear filter problem can be summarized as a weighted integral in the form of "*non-linear* × *Gaussian density*". that is:

$$I(f) = \int\_{\Omega} f(\mathbf{x}) \omega(\mathbf{x}) d\mathbf{x} \tag{5}$$

where *<sup>f</sup>*(*x*) is a non- linear function, <sup>Ω</sup> <sup>⊆</sup> *<sup>R</sup><sup>n</sup>* is the region of integration, <sup>ω</sup>(*x*) <sup>=</sup> exp <sup>−</sup>*xT<sup>x</sup>* is a Gaussian density and satisfies the non-negativity condition in the entire region.

CKF uses the third-degree spherical-radial rule to calculate the non-linear filter problem [17]. Based on the third-degree spherical-radial rule, a set of 2n equal weight sampling points is used to approximate the integral value. that is:

$$I(f) \approx \sum\_{i=1}^{2n} a\_i \mathfrak{f}(\xi\_i) \tag{6}$$

where <sup>ξ</sup>*<sup>i</sup>* <sup>=</sup> <sup>√</sup> *<sup>n</sup>*[**1**]*<sup>i</sup>* is the cubature point, <sup>ω</sup>*<sup>i</sup>* = <sup>1</sup> <sup>2</sup>*<sup>n</sup>* is the corresponding weight, *i* = 1, ... , 2*n*, [**1**]*<sup>i</sup>* is the no.i element of the cubature points set:

$$[\mathbf{1}]\_i = \left\{ \left[ \begin{array}{c} 1 \\ 0 \\ \vdots \\ 0 \end{array} \right] \left[ \begin{array}{c} 0 \\ 1 \\ \vdots \\ 0 \end{array} \right] \dots \left[ \begin{array}{c} 0 \\ 0 \\ \vdots \\ 1 \end{array} \right] \left[ \begin{array}{c} -1 \\ 0 \\ \vdots \\ 0 \end{array} \right] \left[ \begin{array}{c} 0 \\ -1 \\ \vdots \\ 0 \end{array} \right] \dots \left[ \begin{array}{c} 0 \\ 0 \\ \vdots \\ -1 \end{array} \right] \tag{7}$$

#### *3.2. Square-Root Cubarure Kalman Filter*

In CKF, the error covariance matrix needs to satisfy two conditions: (1) symmetry (2) positive qualitative. It is important to ensure these two points in the algorithm iteration process to improve the robustness of the filter. The robustness of algorithm is defined as free of failure that the parameters in the algorithm violate the constraints, which makes the algorithm unable to continue running.

*Sensors* **2020**, *20*, 2620

In order to ensure the symmetry and positive qualitative of the error covariance matrix, Arasaratnam proposed the SRCKF algorithm based on the CKF algorithm [17]. The SRCKF algorithm process is as follows:

(1) Initialization

*<sup>x</sup>*ˆ*<sup>k</sup>* is the state variable, *<sup>P</sup><sup>k</sup>* is the error covariance matrix, *<sup>Q</sup><sup>k</sup>* is the process noise, *<sup>R</sup><sup>k</sup>* is the measurement noise.

(2) Evaluate the cubature points:

$$X\_{j,k} = S\_k \mathbb{1}\_j + \mathfrak{k}\_k; \ j = 1, 2, \dots, 2n \tag{8}$$

where ξ*<sup>j</sup>* is the cubature points set, as shown below:

$$\mathfrak{E}\_{\flat} = \begin{cases} \sqrt{n}[\mathbf{1}]\_i & i = 1, 2, \dots, n \\ -\sqrt{n}[\mathbf{1}]\_i & i = n+1, n+2, \dots, 2n \end{cases}$$

where [**1**] is the unit matrix.

(3) Spread cubature points and calculate the state prediction:

$$\mathbf{X}\_{j,k+1}^\* = f(\mathbf{X}\_{j,k}, \boldsymbol{\mu}\_k) \tag{9}$$

$$\overline{\mathbf{x}}\_{k+1} = \frac{1}{2n} \ast \sum\_{j=1}^{2n} \mathbf{X}\_{j,k+1}^{\ast} \tag{10}$$

(4) Estimate the square-root factor of the predicted error covariance

$$\overline{\mathbf{S}}\_{k+1} = \operatorname{Trid} \Big( \left[ \mathbf{x}\_{k+1}^\* \mathbf{S}\_{Q,k+1} \right] \Big) \tag{11}$$

where *Tria*() represents a triangular operation. *B* = *Tria*(*A*) means: *B* is a matrix of *A<sup>T</sup>* upper triangular matrix obtained *QR* decomposition.

χ∗ *<sup>k</sup>*+<sup>1</sup> and *SQ*,*k*+<sup>1</sup> are as follows:

$$\mathbf{x}\_{k+1}^{\*} = 1/\left(\sqrt{2n}\right) \begin{bmatrix} \mathbf{X}\_{1,k+1}^{\*} - \overline{\mathbf{x}}\_{k+1} \ \mathbf{X}\_{2,k+1}^{\*} - \overline{\mathbf{x}}\_{k+1} \cdot \cdots \ \mathbf{X}\_{2n,k+1}^{\*} - \overline{\mathbf{x}}\_{k+1} \end{bmatrix}$$

$$\mathbf{S}\_{Q,k+1} = \text{chol}\left(\mathbf{Q}\_{k+1}\right)$$

where chol() is Cholesky decomposition.

(5) Recalculate the cubature points

$$X\_{j,k+1} = \overline{\mathbf{S}}\_{k+1} \underline{\mathbf{x}}\_{j} + \overline{\mathbf{x}}\_{k+1} \tag{12}$$

(6) Spread cubature points

$$Z\_{j,k+1} = h(\mathbf{X}\_{j,k+1}, \mathbf{u}\_{k+1}) \tag{13}$$

(7) Observation Prediction:

$$\overline{z}\_{k+1} = \frac{1}{2n} \sum\_{j=1}^{2n} Z\_{j,k+1} \tag{14}$$

(8) Estimate the square-root of the innovation covariance matrix:

$$\mathbf{S}\_{\rm ZZ,k+1} = \operatorname{Trid} \Big( \left[ \mathbf{y}\_{k+1} \; \mathbf{S}\_{\mathbf{R},k+1} \right] \Big) \tag{15}$$

where γ*k*+<sup>1</sup> and *S<sup>R</sup>* are as follows:

$$\mathbf{y}\_{k+1} = 1/\left(\sqrt{2n}\right) \begin{bmatrix} \mathbf{Z}\_{1,k+1} - \mathbf{\overline{z}}\_{k+1} \ \mathbf{Z}\_{2,k+1} - \mathbf{\overline{z}}\_{k+1} \cdot \cdots \ \mathbf{\overline{z}}\_{2u,k+1} - \mathbf{\overline{z}}\_{k+1} \end{bmatrix}$$

$$\mathbf{S}\_{\mathbf{R},k+1} = \operatorname{chol}(\mathbf{R}\_{k+1})$$

(9) Estimate the cross-covariance matrix

$$P\_{\mathbf{X}\mathbf{Z},k+1} = \mathbf{x}\_{k+1}\mathbf{y}\_{k+1}^T \tag{16}$$

where χ*k*+<sup>1</sup> is as follows:

$$\chi\_{k+1} = 1/\left(\sqrt{2n}\right) \left[ \mathbf{X}\_{1,k+1} - \overline{\mathbf{x}}\_{k+1}, \mathbf{X}\_{2,k+1} - \overline{\mathbf{x}}\_{k+1}, \dots, \mathbf{X}\_{2n,k+1} - \overline{\mathbf{x}}\_{k+1} \right]$$

(10) Estimate the Kalman gain

$$\mathbf{K}\_{k+1} = (\mathbf{P}\_{\mathbf{X}\mathbf{Z},k+1}/\mathbf{S}\_{\mathbf{z}\mathbf{z},k+1}^T) / \mathbf{S}\_{\mathbf{z}\mathbf{z},k+1} \tag{17}$$

(11) Estimate the updated state

$$
\mathfrak{X}\_{k+1} = \overline{\mathfrak{x}}\_{k+1} + \mathbb{K}\_{k+1} (z\_{k+1} - \overline{z}\_{k+1}) \tag{18}
$$

(12) Estimate the square-root factor of the corresponding error covariance

$$\mathbf{S}\_{k+1} = \operatorname{Train}([\mathbf{x}\_{k+1} - \mathbf{K}\_{k+1} \mathbf{y}\_{k+1} \mathbf{K}\_{k+1} \mathbf{S}\_{\mathbf{R}}]) \tag{19}$$

#### *3.3. Improved Square-Root Cubature Kalman Filter*

The on-board millimeter-wave radar is driving in the roads, and the statistical parameters of measurement noise are often unknown and time-varying. If the values of measurement noise and process noise are not consistent with the actual noise statistics, the filter will diverge. Therefore, constructing the ASRCKF for online estimation of unknown noise statistics is of great significance for improving the accuracy of filters. Sage-Husa is often used to estimate the statistical characteristics of noise online because of its simplicity and good real-time performance. However, the conventional Sage–Husa noise statistical estimator is suitable for the constant coefficient noise statistical characteristic estimation of linear systems. In the literature [26], the Sage–Husa noise statistic estimator in linear system is extended to the non-linear Sage–Husa noise statistic estimator. Therefore, based on the literature [26], we combine the Sage–Husa noise statistical estimator and cubature rules to derive a time-varying noise statistical estimator suitable for nonlinear systems.

The noise statistical estimator based on Sage–Husa in a non-linear system is:

Recursive formula for the process noise means:

$$\hat{\eta}\_{k+1} = \frac{1}{k+1} [k \ast \eta\_k + \hat{\mathbf{x}}\_{k+1} - \frac{1}{2n} \sum\_{i=1}^{2n} f\_k(\mathbf{X}\_{i,k|k})] \tag{20}$$

Recursive formula for the measurement noise means:

$$\hat{\mathbf{r}}\_{k+1} = \frac{1}{k+1} [k \ast \hat{\mathbf{r}}\_k + \mathbf{z}\_{k+1} - \frac{1}{2n} \sum\_{i=1}^{2n} h(\mathbf{X}\_{i,k+1|k})] \tag{21}$$

The process noise variance is written in the form of a recursive estimation formula as:

$$\begin{split} \hat{\mathbf{Q}}\_{k+1} &= \frac{1}{k+1} [\mathbf{k} \ast \hat{\mathbf{Q}}\_{k} + \mathbf{K}\_{k+1} \mathbf{v}\_{k+1} \mathbf{v}\_{k+1}^{T} \mathbf{K}\_{k+1}^{T} + \mathbf{P}\_{k+1} \\ &- (\frac{1}{2\overline{n}} \sum\_{i=1}^{2n} \mathbf{X}\_{i,k+1|k}^{\star} \mathbf{X}\_{i,k+1|k}^{\star T} - \hat{\mathbf{x}}\_{k+1|k} \hat{\mathbf{x}}\_{k+1|k}^{T}) \big] \end{split} \tag{22}$$

The measurement noise variance is written in the form of a recursive estimation formula as:

$$\mathcal{R}\_{k+1} = \frac{1}{k+1} [\mathbf{k} \ast \mathbf{\hat{R}}\_k + \mathbf{v}\_{k+1} \mathbf{v}\_{k+1}^T - (\frac{1}{2n} \sum\_{i=1}^{2n} \mathbf{Z}\_{i,k+1|k} \mathbf{Z}\_{i,k+1|k}^T - \mathbf{\hat{z}}\_{k+1|k} \mathbf{\hat{z}}\_{k+1|k}^T)] \tag{23}$$

Including: ˆ*q*, ˆ*r*, *Q*ˆ , *R*ˆ are the maximum a posterior of *q*, *r*, *Q*, *R* respectively.

As can be seen from (20), (21), (22), (23), the weight coefficient of each factor in *k* + 1 moments, every factor of the weighted coefficient of *q*ˆ, *r*ˆ, *Q*ˆ , *R*ˆ is 1/(k + 1). For time-varying noise, the role of new data should be increased, while the role of old data should be gradually forgotten. Therefore, a different weighting factor should be multiplied for each factor in the system noise. That is: the weighting coefficient of new data should be greater than that of the old data. On the basis of the constant noise statistic estimator, the time-varying noise statistic estimator suitable for SRCKF is deduced by using the method of fading memory index weighting.

Weighted coefficient: <sup>λ</sup>*<sup>i</sup>* <sup>=</sup> <sup>λ</sup>*i*−<sup>1</sup> <sup>∗</sup> *<sup>b</sup>*, 0 <sup>&</sup>lt; *<sup>b</sup>* <sup>&</sup>lt; 1, <sup>5</sup>*k*+<sup>1</sup> *<sup>i</sup>*=<sup>1</sup> λ*<sup>i</sup>* = 1.

Therefore, the weighted index of fading memory is:

$$\begin{cases} \lambda\_i = d\_k b^{i-1} \\ \quad d\_k = \frac{1-b}{1-b^k} \end{cases} \tag{24}$$

where b is the forgetting factor, and its value range is usually between 0.95 and 0.99. A time-varying noise statistical estimator is obtained by replacing λ*k*+1−*<sup>i</sup>* with the weight factor of 1/(*k* + 1) in the constant noise statistical estimator:

Recursive formula for the process noise means:

$$\mathfrak{q}\_{k+1} = (1 - d\_{k+1})\mathfrak{q}\_k + d\_{k+1} [\mathfrak{x}\_{k+1|k+1} - \frac{1}{2n} \sum\_{i=1}^{2n} f\_k(\mathbf{X}\_{i,k|k})] \tag{25}$$

Recursive formula for the measurement noise means:

$$\mathfrak{H}\_{k+1} = (1 - d\_{k+1})\mathfrak{H}\_k + d\_{k+1} [z\_{k+1} - \frac{1}{2n} \sum\_{i=1}^{2n} h\_{k+1} \{ \mathbf{X}\_{i,k+1|k} \} ] \tag{26}$$

The process noise variance is written in the form of a recursive estimation formula as:

$$\begin{split} \hat{\mathbf{Q}}\_{k+1} &= (1 - d\_{k+1})\hat{\mathbf{Q}}\_{k} + d\_{k+1} [\mathbf{K}\_{k+1} \boldsymbol{\varepsilon}\_{k+1} \boldsymbol{\varepsilon}\_{k+1}^T \mathbf{K}\_{k+1}^T + \mathbf{P}\_{k+1} \\ &- (\frac{1}{2u} \sum\_{i=1}^{2u} \mathbf{X}\_{i,k+1|k}^\* \mathbf{X}\_{i,k+1|k}^{\*T} - \mathbf{\hat{x}}\_{k+1|k} \boldsymbol{\hat{\mathbf{x}}}\_{k+1|k}^T)] \end{split} \tag{27}$$

where

$$\mathbf{z}\_k = \mathbf{z}\_k - \mathbf{\hat{z}}\_{k|k-1}$$

The measurement noise variance is written in the form of a recursive estimation formula as:

$$\hat{\mathcal{R}}\_{k+1} = (1 - d\_{k+1})\hat{\mathcal{R}}\_k + d\_{k+1} [\varepsilon\_{k+1} \varepsilon\_{k+1}^T - (\frac{1}{2n} \sum\_{i=1}^{2n} Z\_{i,k+1|k} \mathbf{Z}\_{i,k+1|k}^T - \mathfrak{z}\_{k+1|k} \hat{\mathfrak{z}}\_{k+1|k}^T)] \tag{28}$$

#### **4. Classification of Target-Vehicle Motion State**

Due to the influence of the ego-vehicle speed sensor and millimeter-wave radar measurement error, the direct use of the current moment of the ego vehicle and target vehicle motion relationship to identify the target-vehicle motion state, will lead to the vibration and even inaccurate motion state classification results. A method of classifying the motion state of the target vehicle based on time window is proposed by analyzing the transfer mechanism of the motion state of the target-vehicle. According to the absolute velocity of the target vehicle, the motion state of the target vehicle is divided into stationary target vehicle, moving target vehicle, oncoming target vehicle, start-stop target vehicle, and unclassified target vehicle.

(1) Unclassified target vehicle: the motion state of the target vehicle obtained at the initial moment of radar is the unclassified target vehicle;

(2) Stationary target vehicle: the target vehicle whose absolute speed stays near zero for a long time;

(3) Moving the target vehicle in the same direction: the movement direction of the target vehicle is the same as that of the ego vehicle;

(4) Oncoming target vehicle: the movement direction of the target vehicle is opposite to that of the ego vehicle;

(5) Start-stop target vehicle: the speed of the moving target vehicle (or the oncoming target vehicle) is reduced to near zero.

Since the velocity measured by the on-board millimeter-wave radar is the relative motion velocity of the target vehicle relative to the ego vehicle. Therefore, the absolute velocity of the target vehicle relative to the geodetic coordinate system can be deduced:

$$V\_1 = V\_{\mathbf{r}} + V\_{\mathbf{v}} \tag{29}$$

where:

*V*1: The absolute velocity of the target vehicle;

*Vr*: The relative velocity of the target vehicle;

*Vv*: The speed of the ego vehicle.

Figure 2 shows the flow chart of movement state transfer of the target vehicle.

**Figure 2.** The flow chart of movement state transfer of the target vehicle.

According to the absolute speed of the target vehicle, the determination of the target motion state is mainly influenced by the following two factors: first, the measurement error of the ego vehicle speed sensor; Second, millimeter-wave radar speed error. Due to the above two measurement errors, the stationary target may also return a non-zero velocity value. Therefore, it is essential to determine the appropriate threshold value to judge the target motion state. The influence of velocity sensor measurement error and millimeter-wave radar measurement error is considered. In this paper, the reference ranges of the velocity threshold are [ 0.5 ∼ 1 m/s] and [−1 ∼ −0.5 m/s]. Because the fluctuation range of the ego-vehicle velocity is 0.3 m/s. Therefore, the threshold value of the velocity of the stationary target-vehicle is set to ±0.5*m*/*s*. The moving state transition rules of the target vehicle are as follows:

(1) The motion state of the target vehicle obtained during the initial operation of the radar is unclassified;

(2) The absolute speed of the target vehicle is between [−0.5 ∼ 0.5*m*/*s*] for *N* consecutive cycles. The target motion state transition can be divided into the following four conditions:





(3) When the absolute speed of the target-vehicle is greater than 0.5m/s for *N* consecutive cycles, the target motion state transition has the following four conditions:





(4) When the absolute speed of the target-vehicle is less than –0.5 m/s for *N* consecutive cycles, the target motion state transition has the following four conditions:





(5) By recording the time *T* when the target vehicle is recognized as the start-stop motion state, the transition relationship between the start-stop motion state and the stationary state is identified.



As the length of time window *N* is longer, the delay of target-vehicle motion state recognition is more serious. The value of *M* has a significant influence on decision-making and control of the vehicle. In this paper, *N* = 3, *M* = 2*s*.

Because the target vehicle has inertia, there is no sudden change in the speed of the target-vehicle. In the process of state transfer between moving target vehicles in the same direction and moving target vehicles in the opposite direction, it is necessary to go through the start-stop motion state.

#### **5. Experiment and Discussion**

#### *5.1. Construction of the Experimental Platform*

As shown in the Figure 3, to check the accuracy and reliability of the proposed algorithm, the Sagitar vehicle is applied herein in the test. Figure 3 (a) is the Sagitar test vehicle platform; (b) is the experimental platform communication.

**Figure 3.** Experimental vehicle. (**a**) Test platform equipment, (**b**) experimental platform communication.

The ego-vehicle is equipped with a 77 GHz millimeter-wave radar as shown in (a) and (b) of Figure 3. The 77 GHz millimeter-wave radar is installed directly above the license plate on the front of the ego-vehicle. The 77 GHz millimeter-wave radar is provided with two-way controller area network (CAN). One CAN is connected to the vehicle gateway. The other is connected to MicroAutoBoxII through CAN. The AR023Z-1080p camera is installed in the bracket above the ego-vehicle to synchronously record the scene of the vehicle.

We use the light detection and ranging (lidar) produced by Ibeo to detect the motion state information of the target vehicle and take it as the ground truth, which is used to verify the performance of the target tracking algorithm of the millimeter-wave radar. In order to make the sensors smarter, Ibeo will provide point cloud processing algorithm software for the lidar. At present, the algorithm provided by Ibeo supports target recognition and tracking. The motion state information of the target vehicle includes longitudinal distance and longitudinal speed. As shown in Figure 4a, the ego-vehicle is equipped with two lidars, which are lidar-1 and lidar-2. The Ibeo TrackingBox produced by Ibeo is responsible for data fusion of the two lidars. Table 2 shows the specific technical specifications of the four scan levels lidar provided by Ibeo.


**Table 2.** Main technical parameters of the lidar.

The ego-vehicle is equipped with MicroAutoBoxII and connected to the vehicle gateway port via the CAN to obtain the longitudinal velocity and steering wheel angle information of the vehicle measured by the electronic stability controller (ESC). MicroAutoBoxII is connected to the Ibeo TrackingBox through CAN and obtains the target vehicle movement status after data fusion.

The motion state estimation algorithm of the target vehicle is written in the environment of MatlabR2018a/Simulink in the host computer. The automatic code generation software provided by the dSPACE company is used to download the code to the MicroAutoBoxII1401/ 1505/1507 rapid prototyping controller through user datagram protocol (UDP) to run.

In order to exclude the randomness of the experiments, we conducted multiple sets of experiments. As shown in Figure 4, four test environments are selected, and three groups of experiments are carried out for each test environment.

(**a**) (**b**)

**Figure 4.** Experimental environment. (**a**) underground parking, (**b**) tunnel, (**c**) campus, (**d**) expressways.

#### *5.2. Analysis of Experiment Results*

Since the test results for the 12 groups of experiments are similar, one set of test data is selected for discussion and analysis. Figure 5 to Figure 6 show the experimental results of the effect of verifying the ISRCKF target tracking algorithm.

**Figure 5.** Target-vehicle longitudinal distance time history curve.

**Figure 6.** Target-vehicle longitudinal velocity time history curve.

The measurement data of the target vehicle are indicated by the red circle, without the ISRCKF algorithm optimization and containing noise (R-M). The solid blue line represents the ground truth of the motion state of the target-vehicle measured by lidar. The measurement data of the target vehicle are indicated by the magenta dashed line, with the ISRCKF algorithm optimization. The measurement data of the target vehicle are indicated by the green solid line, with the SH-EKF algorithm optimization. The measurement data of the target vehicle are indicated by the cyan dot-dash line, with the SRCKF algorithm optimization.

Figures 5 and 6 show the time history curves of the longitudinal distance and longitudinal velocity of the target-vehicle relative to the ego-vehicle. It can be concluded from the figures that the data fluctuation of the ISRCKF algorithm during tracking the target is small, and the target tracking performance of the ISRCKF algorithm is significantly better than the SH-EKF and SRCKF algorithms. Through references [27–29], it is found that the root mean square error (RMSE) can be used to quantitatively analyze the filtering accuracy. In order to quantitatively analyze the filtering accuracy of the ISRCKF, SH-EKF, and SRCKF algorithms, we counted the RMSE of the 12 groups of experiments. and the results are shown in Figures 7 and 8. The test environment of the statistical value for RMSE of groups 1, 2, and 3 is the underground parking. The test environment of the statistical value for RMSE of groups 4, 5, and 6 is the tunnel. The test environment of the statistical value for RMSE of groups 7, 8, and 9 is the campus. The test environment for the statistical value for RMSE of groups 10, 11, and 12 is the expressway.

**Figure 7.** MSRE curves of longitudinal velocity for various filters.

It can be concluded from Figures 7 and 8 that in the driving environment where the system noise is unknown and time-variant, the RMSE of the longitudinal speed and longitudinal distance of target tracking using SRCKF and SH-EKF algorithm is larger than that of ISRCKF algorithm. Compared with SH-EKF and SRCKF algorithms, the maximum increase in filtering accuracy of longitudinal distance using the ISRCKF algorithm is 45.53% and 59.15%, respectively, and the maximum increase in filtering accuracy of longitudinal speed using the ISRCKF algorithm is 23.53% and 29.09%, respectively. Through the above analysis, it can be concluded that using ISRCKF algorithm to track can effectively suppress the divergence of target tracking, thereby reducing the tracking error and improving the tracking accuracy.

As shown in Figure 9, the time consumption of different algorithms is measured by the mean time of the algorithm running once in the MicroAutoBoxII.

**Figure 8.** MSRE curves of longitudinal distance for various filters.

**Figure 9.** Time mean of the algorithm running once.

From Figure 9, it can be concluded that the target tracking algorithm proposed herein takes 0.0053 s on average, slightly up that of SH-EKF and SRCKF. SH-EKF and SRCKF take 0.0051 s and 0.0047 s on average, respectively.

The vehicle test results show that the ISRCKF algorithm has the highest accuracy when the time consumption is not increased much compared with SH-EKF and SRCKF.

In order to ensure the safety and rigorous of the experiment, we conducted 10 groups of experiments on campus. During the test, the driving process of the target vehicle was divided into the following stages: parking, accelerated reverse, braking to stop, starting acceleration, and decelerating to stop. Since the test results of the 10 groups are similar, we choose one of the data for analysis.

Figure 10 is the time history curve of the target-vehicle's absolute velocity. The meanings of -1 , -<sup>2</sup> , -<sup>3</sup> , and -<sup>4</sup> in Figure 10 represent the stationary, moving, start-stop, oncoming movement states in Figure 2, respectively. We have used the same colors to indicate the same motion state in Figures 2 and 10.

Figure 11 is the transition process of the target vehicle's motion state. The initial motion state of the target vehicle is the unclassified motion state. As the reversing speed of the target vehicle increases at 2.3 s, the motion state of the target-vehicle changes from the stationary motion state to the oncoming motion state. As the reversing speed of the target vehicle decreases, the motion state of the target vehicle changes from the oncoming motion state to the start-stop motion state (at 14.5 s). As the time T of the start-stop motion state increases, the status of the target vehicle changes from the start-stop motion state to the stationary state at the time of 16.5 s. As the forward speed of the target vehicle increases, the motion state of the target vehicle changes from the stationary state to the same motion state (at 19.4 s). As the forward speed of the target vehicle decreases, the motion state of the target vehicle changes from the same motion state to the start-stop motion state at 28.4 s.

**Figure 10.** Target-vehicle absolute velocity time history curve.

**Figure 11.** The movement state transfer of the target-vehicle.

In the two sports periods of 14.5 s ∼ 16.5 s and 28.4 s ∼ 30 s, the classification method proposed classifies target vehicles as start-stop targets. The classification method shows that it has a certain memory effect on the motion state of the target vehicle. The test results show that the classification and recognition results are consistent with the motion state of the target vehicle.

#### **6. Conclusions**

In the process of on-board millimeter-wave radar target tracking, the characteristics of unknown time-varying noise cannot be accurately counted, which will cause the filter accuracy to decline or even diverge. In response to this question, based on the SRCKF, the Sage–Husa noise statistic estimator and the fading memory exponential weighting method are combined to derive a time-varying noise statistic estimator for non-linear systems. ISRCKF can effectively overcome the problem of low accuracy and even divergence of SRCKF filtering-varying noise. When the motion relation between the ego vehicle and the target vehicle in the current is used to identify the motion state of the target vehicle, there exists the problem that the classification result of motion state is vibration or even inaccurate. A method of classifying the motion state of the target vehicle based on the time window is proposed by analyzing the transfer mechanism of the motion state of the target vehicle. According to the absolute velocity of the target vehicle, the motion state of the target vehicle is divided into stationary target vehicle, moving target vehicle, oncoming target vehicle, start-stop target vehicle, and unclassified target vehicle. Because the target vehicle has inertia, there is no sudden change in the speed of the target vehicle. In this classification method, the target vehicle needs to go through the start-stop motion state during the state transfer between the same motion and the reverse motion. Since the starting and stopping motion state of the target vehicle is different from the stationary motion state of the target vehicle, this method reflects that the target vehicle has a certain memory effect on its motion state. The results of the vehicle test show that: (1) the accuracy of the ISRCKF algorithm is significantly improved; (2) the classification and recognition results of the target vehicle motion state are consistent with the target vehicle motion state.

**Author Contributions:** Funding acquisition, J.W.; Methodology, S.S.; Writing—original draft, S.S.; Writing—review and editing, S.S. and J.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Research on Construction and Simulation Technology of Hardware in Loop Testing Scenario for Self-driving Electric Vehicle of China, grant number 2018YFC0105103 and the National Natural Science Foundation of China, grant number U1564211.

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

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).
