A Cost-Effective Vehicle Localization Solution Using an Interacting Multiple Model−Unscented Kalman Filters (IMM-UKF) Algorithm and Grey Neural Network

Abstract

In this paper, we propose a cost-effective localization solution for land vehicles, which can simultaneously adapt to the uncertain noise of inertial sensors and bridge Global Positioning System (GPS) outages. First, three Unscented Kalman filters (UKFs) with different noise covariances are introduced into the framework of Interacting Multiple Model (IMM) algorithm to form the proposed IMM-based UKF, termed as IMM-UKF. The IMM algorithm can provide a soft switching among the three UKFs and therefore adapt to different noise characteristics. Further, two IMM-UKFs are executed in parallel when GPS is available. One fuses the information of low-cost GPS, in-vehicle sensors, and micro electromechanical system (MEMS)-based reduced inertial sensor systems (RISS), while the other fuses only in-vehicle sensors and MEMS-RISS. The differences between the state vectors of the two IMM-UKFs are considered as training data of a Grey Neural Network (GNN) module, which is known for its high prediction accuracy with a limited amount of samples. The GNN module can predict and compensate position errors when GPS signals are blocked. To verify the feasibility and effectiveness of the proposed solution, road-test experiments with various driving scenarios were performed. The experimental results indicate that the proposed solution outperforms all the compared methods.

1. Introduction

Accurate and reliable vehicle ego-position is important and necessary information in more and more Intelligent Transportation System (ITS) applications [1,2,3]. The most popular technique is Global Positioning System (GPS), which can provide satisfactory localization performance in open areas [4,5], but in modern urban environments, more and more tall buildings or overpasses may affect the GPS signals and cause the failure of GPS. To improve the GPS localization performance, it is usually integrated with Inertial Navigation System (INS), which is a self-contained system and is not affected by external disturbances [4]. For land vehicles, the intent is to have a low-cost localization system [6,7], so in general only low-cost inertial sensors based on microelectromechanical systems (MEMS) are affordable enough [8,9]. In order to further lower the cost of vehicle localization systems, research efforts have recently been made to investigate the applicability of reduced inertial sensor systems (RISS) [10,11]. Usually, RISS involves a single-axis gyroscope and two-axis accelerometers. In RISS mechanization, pitch and roll are calculated from accelerometers instead of gyroscopes, and vehicle velocity is calculated from the forward speed derived from wheel speed sensors instead of accelerometers.

Although MEMS-based inertial sensors are portable and low-cost, their measurements often suffer from large and uncertain noises [8], which can seriously affect the localization performance. Unscented Kalman filter (UKF) has been widely used to fuse GPS and inertial sensor data in vehicle localization problems [12,13]. It can essentially provide derivative-free high order approximations of nonlinear models [14,15]. However, one of the limitations of UKF is the necessity to have priori statistical information of the process and measurement noises [16]. They are specified in the form of process noise covariances and measurement noise covariances in the filtering algorithm, normally expressed as the Q matrix and R matrix, respectively. Usually, the specification of the R matrix can be directly derived from the accuracy characteristics of the measurement devices, while the specification of the Q matrix is often determined by a trial-and-error approach and considered as a constant [17]. In a number of practical situations, due to the high uncertainty of the MEMS inertial sensor noise, the covariances are variable and difficult to determine. An improper specification will degrade the performance of the filter or even cause divergence. The adaptive filter algorithm has been considered as one of the strategies to adjust the covariance matrices through scale factors [18,19]. However, the approaches for determining the scale factors heavily rely on personal experience or computer simulation [20]. Besides, even small changes of the scale factor can greatly affect the final performance of the filter.

Alternatively, a structural adaptation approach called Interacting Multiple Model (IMM) makes it possible for a set of models with different characteristics to run in parallel [21,22]. IMM algorithms have already been applied to vehicle localization, and are usually used to represent the possible vehicle driving patterns with a set of models, which are generally established according to different maneuvering or driving conditions [21,22,23]. The IMM algorithm has shown better results than conventional switching schemes, because a smooth transition from one model to another is achieved [24]. Different from the common applications, we envisioned that the IMM algorithm can provide a soft switching among the filters designed for different noise levels and contribute to adapt to the uncertain noise of MEMS inertial sensors. To the authors’ knowledge, this application has seldom been evaluated or discussed in the existing literature.

Besides, in order to further compensate the position errors during GPS outages, artificial intelligence (AI) methods have attracted researchers’ interest due to their abilities of modeling and predicting for nonlinear system [25]. Usually, AI methods are used to model position errors when GPS signal is available and predict position errors during GPS outages. These AI-based approaches include Back-Propagation Neural Network (BPNN) [26], Support Vector Machine [27], etc. These approaches generally require a large number of training samples to achieve good generalization performance. However, to meet the real-time requirement of vehicle localization, a sliding window with a certain size is usually adopted to select training samples [28,29]. When there is a limited amount of samples, these approaches have a high probability of being affected [30].

Grey system theory can make full use of the historical data sequence information and is characterized by modeling with insufficient data [31]. However, grey system models have some drawbacks due to the lack of feedback, self-learning, and self-adaption. BPNN is one of the most popular learning algorithms and it can approximate an arbitrary nonlinear function with satisfactory precision [32]. Through mapping a grey system model to a BPNN, we can get a grey neural network (GNN) [33], which sufficiently exploits the advantages of both grey system theory and BPNN. We also develop a GNN module to model and predict position errors in this paper. To the authors’ knowledge, GNN is applied to vehicle localization for the first time.

This paper presents a cost-effective localization solution to adapt to uncertain inertial sensor noise and compensate position errors during GPS outages. The proposed vehicle localization solution integrates low-cost GPS, MEMS-RISS, and in-vehicle sensors. Here, the in-vehicle sensors specifically refer to wheel speed sensors and steering angle sensors. A novel IMM-UKF algorithm is proposed by introducing several UKFs designed for different noise levels into the framework of an IMM algorithm. Two IMM-UKFs are utilized to work in parallel when GPS is available. One fuses all the sensors including GPS, in-vehicle sensors, and RISS, while the other only fuses those which are not affected by the GPS-denied environments. The differences between the state vectors of the two IMM-UKFs are considered as the output of training samples. Meanwhile, the measurements of inertial sensors are considered as the corresponding input. A GNN module is adopted to establish the model and thus the position errors can be predicted and compensated during GPS outages. The proposed localization solution has been extensively evaluated in road-test experiments.

The rest of the paper is structured as follows. In the next section, we briefly describe the overview of the proposed localization solution. Then, in Section 3, we explain the detailed implementation of IMM-UKF algorithm. The design of GNN module is presented in Section 4. Section 5 presents the setup, scenarios, and results of experimental validation. Finally, a conclusion is presented in Section 6.

2. Overview of the Proposed Solution

The whole mechanism and functionality of the proposed localization solution is illustrated in Figure 1, which can be divided into two parts, i.e., the sensor part and the fusion part. In the sensor part, GPS, in-vehicle sensors, and MEMS-RISS are all cost-effective ones. For the MEMS-RISS, one single-axis gyroscope is mounted with its sensitive axis aligned with the vertical axis of the vehicle and measures the rotation rate of yaw angle in the body frame, denoted as ${\omega }_k^z$. Two accelerometers are mounted along the longitudinal and lateral axes of the vehicle, and the longitudinal and lateral accelerations of the vehicle are measured respectively, denoted as $a_k^x$ and $a_k^y$. For the in-vehicle sensors, since more and more vehicles are equipped with Antilock Brake System (ABS) and Electronic Stability Program (ESP), the information about vehicle forward speed and steering angle can be directly obtained from the controller area network (CAN) bus [22]. The wheel speed sensor can provide longitudinal velocity $v_k^{wh}$, while steering angle sensor can derive lateral velocity ${\widehat{v}}_k^y$.

In the fusion part, improving the robustness against the uncertain noise and bridging the GPS outages are the main priorities. Specifically, the proposed IMM-UKF algorithm contains three UKFs designed for different noise characteristics and the output is a weighted sum of the three individual UKFs. The proposed algorithm can adapt to a wide variation of inertial sensor noise. When GPS is available, IMM-UKF1 fuses the information from GPS, in-vehicle sensors, and MEMS-RISS. Because GPS can provide direct position and velocity observations, more accurate vehicle positions can be achieved by IMM-UKF1. Meanwhile, IMM-UKF2 is introduced to work in parallel with IMM-UKF1. The IMM-UKF2 only fuse the information of in-vehicle sensors and MEMS-RISS. In order to bridge GPS outages, a GNN module is employed to establish the model of position errors. The difference between the state vectors of the two IMM-UKFs at each epoch is transferred to the GNN module as the desired output. While the RISS output is fed to the GNN as the corresponding input at the same epoch. Considering the balance between model accuracy and computation efficiency, a sliding window with a certain window size is considered for sample selection [27]. The GNN parameters are continuously updated till the occurrence of GPS outages to ensure the predicting precision.

When satellite signals are blocked, the absence of GPS observations causes the invalidation of IMM-UKF1 and it is removed from the system. However, the IMM-UKF2 can still efficiently execute the fusion of in-vehicle sensors and RISS. As shown in Figure 2, the GNN module can predict the position errors with the input of current RISS measurements. Due to the adaptation of IMM-UKF, even if the uncertain noise causes changes on the statistical properties, the proposed solution will not be affected and can still maintain the performance. Thus, accurate vehicle positions can be obtained even when the localization system suffers from GPS outages and uncertain noise of MEMS inertial sensors simultaneously.

3. Proposed IMM-UKF Algorithm

Since the inertial sensor noise is highly uncertain, a fixed value of noise statistics can lead to poor filter performance and even result in filter divergence. Thus, it is advisable to use IMM, which can represent the noise behavior with different characteristics and provide a soft switching among these noise characteristics. We study the IMM-UKF algorithm to adapt to the uncertain noise of inertial sensors. The details about IMM-UKF are shown in the following paragraphs.

3.1. Motion Model

The nonlinear motion model for the RISS involving attitude, velocity, and position states is presented in this section. When the vehicle is moving, the forward accelerometer measures the forward vehicle acceleration as well as the component due to gravity. Similarly, the transversal accelerometer measures the normal component of the vehicle acceleration as well as the component due to gravity. Thus, the pitch angle can be calculated by removing the vehicle acceleration derived from the wheel speed sensor measurements from the forward accelerometer measurements, while the roll angle can be calculated by compensating the transversal accelerometer measurements for the normal component of acceleration. The equations can be expressed as [34]:

$$\begin{array}{c}{\mathrm{p}}_{\mathrm{k}}={\sin }^{-1}\frac{a_k^x-{\dot{v}}_k^{wh}}{g}\\ {\mathrm{r}}_{\mathrm{k}}=-{\sin }^{-1}\frac{a_k^y+{\omega }_k^z{v}_k^{wh}}{\mathrm{g}\cos p_k}\end{array}$$

(1)

where the subscript $k$ represents the time step, $p_k$ and $r_k$ are pitch and roll angle, respectively.${\dot{v}}_k^{wh}$ is the differentiation of $v_k^{wh}$, $g$ denotes the acceleration due to gravity.

Note that $v_k^{wh}$ is derived from the wheel speed sensor rather than the longitudinal accelerometer. This is because any uncompensated accelerometer bias will introduce an error to the speed during the integration. However, the speed derived from the wheel speed sensor avoids the integration. Besides, ${\dot{v}}_k^{wh}$ at each time step can be calculated as:

$${\dot{v}}_k^{wh}=\frac{v_k^{wh}-v_{k-1}^{wh}}{dt}$$

(2)

where $dt$ is the time interval between $v_{k-1}^{wh}$ and $v_k^{wh}$. Since the output frequency of vehicle speed is 100 Hz, $dt$ is 0.01 s here.

When calculating the azimuth angle, both the Earth’s rotation and the change of orientation of the local-level frame are taken into consideration [35]. Thus, the calculation of azimuth angle $A_k$ is:

$$A_k={\tan }^{-1}\frac{U^E}{U^N}+\left({\omega }^e\sin {\phi }_{k-1}\right)\Delta t+\frac{v_{k-1}^{wh}\sin A_{k-1}\cos p_{k-1}\tan {\phi }_{k-1}}{R_N+h_{k-1}}\Delta t$$

(3)

where:

• $U^E=\sin A_{k-1}\cos p_{k-1}\cos {\gamma }_k^z-\left(\cos A_{k-1}\cos r_{k-1}+\sin A_{k-1}\sin p_{k-1}\sin r_{k-1}\right)\sin {\gamma }_k^z$
• $U^N=\cos A_{k-1}\cos p_{k-1}\cos {\gamma }_k^z-\left(-\sin A_{k-1}\cos r_{k-1}+\cos A_{k-1}\sin p_{k-1}\sin r_{k-1}\right)\sin {\gamma }_k^z$
• ${\gamma }_k^z={\omega }_k^z\Delta t$, ${\omega }^e$ is the Earth’s rotation rate,
• ${\phi }_k$ is the latitude of the vehicle position,
• $h_k$ is the altitude,
• $R_N$ is the normal radius of curvature of the Earth.

Since the vehicle does not jump off the ground during common driving maneuvers [36], the vertical velocity can be presumed to be zero. Thus, the relationship between the vehicle’s velocity in the body frame and in the local-level frame (navigation frame) can be calculated as:

$$\left[\begin{array}{c}{v}_k^E\\ v_k^N\\ v_k^U\end{array}\right]=R_{b,l}\left[\begin{array}{c}{v}_k^x\\ v_k^y\\ 0\end{array}\right]$$

(4)

where $v_k^E$ is the velocity component along the east direction, $v_k^N$ is the velocity component along the north direction, $v_k^U$ is the velocity component along the up direction. $v_k^x$ is the velocity component along the forward longitudinal direction, and it can be calculated from the longitudinally aligned accelerometer as $v_k^x=v_{k-1}^x+a_k^x\Delta t$. $v_k^y$ is the velocity component along the transversal direction, and it can be calculated from the laterally aligned accelerometer as $v_k^y=v_{k-1}^y+a_k^y\Delta t$. $R_{b,l}$ is the rotation matrix that transforms from the vehicle body frame to the local-level frame, and is given as:

$$R_{b,l}=\left[\begin{array}{ccc}\sin A_k\cos p_k& \cos A_k\cos r_k+\sin A_k\sin p_k\sin r_k& \cos A_k\sin r_k-\sin A_k\sin p_k\cos r_k\\ \cos A_k\cos p_k& -\sin A_k\cos r_k+\cos A_k\sin p_k\sin r_k& -\sin A_k\sin r_k-\cos A_k\sin p_k\cos r_k\\ \sin p_k & -\cos p_k\sin r_k& \cos p_k\cos r_k\end{array}\right]$$

Expanding Equation (4), we can get:

$$\begin{gathered} v_k^E=v_k^x\sin A_k\cos p_k+v_k^y\left(\cos A_k\cos r_k+\sin A_k\sin p_k\sin r_k\right) \\ {v}_k^N=v_k^x\cos A_k\cos p_k+v_k^y\left(-\sin A_k\cos r_k+\cos A_k\sin p_k\sin r_k\right) \\ {v}_k^U=v_k^x\sin p_k-v_k^y\cos p_k\sin r_k \end{gathered}$$

(5)

Then, the position calculation can be expressed as:

$$\begin{array}{c}{\phi }_k={\phi }_{k-1}+\frac{v_k^N}{R_M+h_k}\Delta t\\ {\lambda }_k={\lambda }_{k-1}+\frac{v_k^E}{\left(R_N+h_k\right)\cos {\phi }_k}\Delta t\\ h_k=h_{k-1}+v_k^U\Delta t\end{array}$$

(6)

where ${\lambda }_k$ is the longitude of the vehicle position, $R_M$ is the meridian radius of curvature of the Earth. Based on Equations (1), (3), (5) and (6), the discrete-time system state equation can be presented as:

$$X_k=f\left(X_{k-1},u_k\right)+W_k$$

(7)

where $X_k$ and $u_k$ represent the state vector and the input vector respectively, $W_k$ is the corresponding system noise vector, $f(·)$ is the nonlinear system function. $X_k$, $u_k$, and $f(·)$ can be described as:

$$X_k=\left[\begin{array}{ccc}\begin{array}{ccc}{\phi }_k& {\lambda }_k& h_k\end{array}& \begin{array}{ccc}{v}_k^E& v_k^N& v_k^U\end{array}& \begin{array}{ccc}{p}_k& r_k& A_k\end{array}\end{array}\right]\prime$$

(8)

$$u_k=\left[\begin{array}{cc}\begin{array}{cc}{v}_k^{od}& a_k^x\end{array}& \begin{array}{cc}{a}_k^y& {\omega }_k^z\end{array}\end{array}\right]\prime$$

(9)

$$f\left(·\right)=\left[\begin{array}{c}{\phi }_{k-1}+\frac{v_k^N}{R_M+h_k}\Delta t\\ {\lambda }_{k-1}+\frac{v_k^E}{\left(R_N+h_k\right)\cos {\phi }_k}\Delta t\\ h_{k-1}+v_k^U\Delta t\\ v_k^x\sin A_k\cos p_k+v_k^y(\cos A_k\cos r_k+\sin A_k\sin p_k\sin r_k\\ v_k^x\cos A_k\cos p_k+v_k^y\left(-\sin A_k\cos r_k+\cos A_k\sin p_k\sin r_k\right)\\ v_k^x\sin p_k-v_k^y\cos p_k\sin r_k\\ {\sin }^{-1}\frac{a_k^x-{\dot{v}}_k^x}{g}\\ -{\sin }^{-1}\frac{a_k^y+{\omega }_k^z{v}_k^x}{g\cos p_k}\\ {\tan }^{-1}\left(\frac{U^E}{U^N}\right)+\left({\omega }^e\sin {\phi }_{k-1}\right)\Delta t+\frac{v_{k-1}^E\tan {\phi }_{k-1}}{R_N+h_{k-1}}\Delta t\end{array}\right]$$

(10)

3.2. Observation Model

As shown in Figure 1, the observation information comes from two sources, i.e., the in-vehicle sensors and the GPS. The observation equation of the in-vehicle sensors can be established as:

$$Z_k^1=h_1\left(X_k\right)+n_1=\left[\begin{array}{c}{v}_k^E\cos A_k+v_k^N\sin A_k+n_{vx}\\ -v_k^E\sin A_k+v_k^N\cos A_k+n_{vy}\end{array}\right]$$

(11)

where $Z_k^1=\left[\begin{array}{cc}{v}_k^{wh}& {\widehat{v}}_k^y\end{array}\right]\prime$, ${\widehat{v}}_k^y$ is the vehicle lateral velocity derived from front wheel steering angle data. $n_1=\left[\begin{array}{cc}{n}_{vx}& n_{vy}\end{array}\right]\prime$ is the corresponding observation noise vector.

In order to estimate the lateral velocity from front wheel steering angle data in real time, we adopt the simple but effective bicycle model [37,38]. Moreover, to achieve more accurate estimation, the influences of roll and pitch angles are also considered. Assume that inner tires and outer tires have the same tire cornering stiffnesses and tire slip angles, the equations for the lateral motion of the vehicle can be established according to Newton’s law of motion [34], described as:

$$m\left({\dot{v}}_k^y+{\omega }_k^z{v}_k^{wh}\right)=2F_k^{sf}+2F_k^{sr}$$

(12)

where $m$ is the mass of the vehicle. $F_k^{sf}$ and $F_k^{sr}$ are the front-tire lateral force and the rear-tire lateral force, respectively. The tire slip is usually small and the tire lateral forces $F_{sf}$ and $F_{sr}$ can usually be approximated by a linear function [22], expressed as:

$$F_k^{sf}=C_{af}{\alpha }_k^f,F_k^{sr}=C_{ar}{\alpha }_k^r$$

(13)

where $C_{af}$ and $C_{ar}$ are the front tire cornering stiffness and rear tire cornering stiffness, respectively. ${\alpha }_k^f$ and ${\alpha }_k^r$ are the front-tire slip angle and rear-tire slip angle respectively, and they can be described as:

$${\alpha }_k^f={\delta }_k^f-\frac{{\widehat{v}}_k^y-a{\omega }_k^z}{v_k^{wh}},{\alpha }_k^r=\frac{b{\omega }_k^z-{\widehat{v}}_k^y}{v_k^{wh}}$$

(14)

where ${\delta }_k^f$ denotes the front-wheel steering angle, $a$ and $b$ are the distance between the center of gravity (CoG) and the front axle and the distance between the CoG and the rear axle, respectively.

Substituting Equations (13) and (14) into Equation (12), we can obtain:

$${\dot{v}}_k^y=\frac{{\beta }_1}{v_k^{wh}}{\widehat{v}}_k^y+\left(\frac{{\beta }_2}{v_k^{wh}}-v_k^{wh}\right){\omega }_k^z+{\beta }_3{\delta }_k^f$$

(15)

where ${\beta }_1=\frac{-2\left(C_{af}+C_{ar}\right)}{m}$, ${\beta }_2=\frac{2\left(C_{ar}b-C_{ar}a\right)}{m}$, ${\beta }_3=\frac{2C_{af}}{m}$.

Thus, the ${\widehat{v}}_k^y$ can be calculated as:

$${\widehat{v}}_k^y={\widehat{v}}_{k-1}^y+{\dot{v}}_k^y\Delta t$$

(16)

Furthermore, the observation equation of the GPS measurements is:

$$Z_k^2=h_2\left(X_k\right)+n_2=\left[\begin{array}{c}{\phi }_k+n_{\phi }\\ {\lambda }_k+n_{\lambda }\\ h_k+n_h\\ v_k^E+n_{vE}\\ v_k^N+n_{vN}\\ v_k^U+n_{vU}\end{array}\right]$$

(17)

where $Z_k^2=[\begin{array}{cccccc}{\phi }_k^G& {\lambda }_k^G& h_k^G& v_k^{GE}& v_k^{GN}& v_k^{GU}\end{array}]\prime$. ${\phi }_k^G$, ${\lambda }_k^G$, and $h_k^G$ are the latitude, longitude, and altitude output by GPS, respectively. $v_k^{GE}$, $v_k^{GN}$, and $v_k^{GU}$ are the east, north and up velocity measured by GPS, respectively. $n_2=[\begin{array}{cccccc}{n}_{\phi }& n_{\lambda }& n_h& n_{vE}& n_{vN}& n_{vU}\end{array}]\prime$ is the corresponding observation noise vector.

3.3. Implementation of the Proposed Algorithm

In our study, the IMM-UKF approach contains three UKFs with different Q matrices, as shown in Figure 3. Using the system state equation and measurement equation described above, we can execute the recursive procedure of the proposed IMM-UKF algorithm, which can be described in four parts [22,23]:

(1) Interaction

The individual filter estimation $X_{k-1}^i$ of the ith UKF (i = 1,2,3) is mixed with the predicted model probability ${\mu }_{k-1}^i$ and the Markov transition probability ${\pi }_{ji}$, i.e., the probability of the transition from state j to state i, to give:

$${\mu }_{k,k-1}^i={{\displaystyle \sum }}_{j=1}^3{\pi }_{ji}{\mu }_{k-1}^j(i=1,2,3)$$

(18)

The mixing weight is given by:

$${\mu }_{k-1}^{j|i}=\frac{{\pi }_{ji}{\mu }_{k-1}^j}{{\mu }_{k,k-1}^i}(i,j=1,2,3)$$

(19)

The mixing of the state estimates ${\overline{X}}_{k-1}^i$ can be computed as:

$${\overline{X}}_{k-1}^i={{\displaystyle \sum }}_{j=1}^3{\mu }_{k-1}^{j|i}{X}_{k-1}^j(i=1,2,3)$$

(20)

The mixing of the covariance ${\overline{P}}_{k-1}^i$ is given as:

$${\overline{P}}_{k-1}^i={{\displaystyle \sum }}_{j=1}^3{\mu }_{k-1}^{j|i}\left\{P_{k-1}^j+\left[{\overline{X}}_{k-1}^i-X_{k-1}^j\right]\left[{\overline{X}}_{k-1}^i-X_{k-1}^j\right]\prime \right\}(i=1,2,3)$$

(21)

(2) Specific Filtering

Using the mixing state and covariance obtained in the interaction step, each UKF predicts and updates the model state and covariance individually. Since the specification of Q matrix depends on the noise characteristics of inertial sensors [39], UKF1 is designed for high-level noise with Q¹, UKF2 is designed for medium-level noise with Q², and UKF3 is designed for low-level noise with Q³. The execution of the ith UKF (i = 1,2,3) can be described as follows:

Step 1: Calculate the Sigma Points

The Cholesky factorization is utilized in obtaining the sigma points, which is numerical efficient and stable, given by:

$$\{\begin{array}{c}{\xi }_{k-1}^{iq}={\overline{X}}_{k-1}^{i}q=0\\ {\xi }_{k-1}^{iq}={\overline{X}}_{k-1}^i+\sqrt{\left(n+\eta \right)}{\left\{chol\left({\overline{P}}_{k-1}^i\right)\right\}}_q^{\prime }q=1,2,\dots ,n\\ {\xi }_{k-1}^{iq}={\overline{X}}_{k-1}^i-\sqrt{\left(n+\eta \right)}{\left\{chol\left({\overline{P}}_{k-1}^i\right)\right\}}_{q-n}^{\prime }q=n+1,n+2,\dots ,2n\end{array}$$

(22)

where $n$ is the dimension of state vector $X$, $\eta ={\alpha }_1^2\left(n+{\alpha }_2\right)-n$ is a scaling parameter ${\alpha }_1$ determines the spread of the sigma points around ${\overline{X}}_{k-1}^i$ and is usually set to a small positive value. ${\alpha }_2$ is a secondly scaling parameter. $\left\{chol\left({\overline{P}}_{k-1}^i\right)\right\}\prime$ is the lower-triangular matrix of the Cholesky factorization of ${\overline{P}}_{k-1}^i$, the subscript q means the qth column.

Step 2: Time Propagation

$${\xi }_{k,k-1}^{iq}=f\left({\xi }_{k-1}^{iq},u_{k-1}\right)q=0,1,\dots ,2n$$

(23)

$${\widehat{X}}_{k,k-1}^i={\displaystyle \sum _{q=0}^{2n}}{\omega }_q^{\left(m\right)}{\xi }_{k,k-1}^{iq}$$

(24)

$$P_{k.k-1}^i={\displaystyle \sum }_{q=0}^{2n}{\omega }_q^{\left(c\right)}\left[{\xi }_{k,k-1}^{iq}-{\widehat{X}}_{k,k-1}^i\right]·{\left[{\xi }_{k,k-1}^{iq}-{\widehat{X}}_{k,k-1}^i\right]}^{\prime }+Q^i$$

(18)

$${\zeta }_{k,k-1}^{iq}=h\left({\xi }_{k,k-1}^{iq}\right)q=0,1,\dots ,2n$$

(26)

$${\widehat{Z}}_{k,k-1}^i={\displaystyle \sum }_{q=0}^{2n}{\omega }_q^{\left(m\right)}{\zeta }_{k,k-1}^{iq}$$

(27)

where the weighting factors are calculated as:

$$\{\begin{array}{c}{\omega }_0^{\left(m\right)}=\frac{\eta }{n+\eta }\\ {\omega }_0^{\left(c\right)}=\frac{\eta }{n+\eta }+\left(1-{{\alpha }_1}^2+{\alpha }_3\right)\\ {\omega }_q^{\left(m\right)}={\omega }_q^{\left(m\right)}=\frac{1}{2\left(n+\eta \right)}q=1,2,\dots ,2n\end{array}$$

${\alpha }_3$ is used to incorporate prior knowledge of the distribution of ${\overline{X}}_{k-1}^i$ and is optimally set to 2 for Gaussian distributions. Note that $h(·)$ is the combination of $h_1(·)$ and $h_2(·)$ in IMM-UKF1, while it is equal to $h_1(·)$ in IMM-UKF2.

Step 3: Measurement Update

$$P_{ZZ}^i={\displaystyle \sum }_{q=0}^{2n}{\omega }_q^{\left(c\right)}\left[{\zeta }_{k,k-1}^{iq}-{\widehat{Z}}_{k,k-1}^i\right]·{\left[{\zeta }_{k,k-1}^{iq}-{\widehat{Z}}_{k,k-1}^i\right]}^{\prime }+R$$

(28)

$$P_{XZ}^i={\displaystyle \sum }_{i=0}^{2n}{\omega }_q^{\left(c\right)}\left[{\xi }_{k,k-1}^{iq}-{\widehat{X}}_{k,k-1}^i\right]·{\left[{\zeta }_{k,k-1}^{iq}-{\widehat{Z}}_{k,k-1}^i\right]}^{\prime }$$

(29)

$$K_k^i=P_{XZ}^i{\left(P_{ZZ}^i\right)}^{-1}$$

(30)

$$X_k^i={\widehat{X}}_{k,k-1}^i+K_k^i\left(Z_k-{\widehat{Z}}_{k,k-1}^i\right)$$

(31)

$$P_k^i=P_{k.k-1}^i-K_k^i{P}_{ZZ}^i{\left(K_k^i\right)}^{\prime }$$

(32)

Note that $Z_k$ is the combination of $Z_k^1$ and $Z_k^2$ in IMM-UKF1, while it is equal to $Z_k^1$ in IMM-UKF2.

(3) Model probability Update

Under Gaussian statistics assumption, the likelihood for the observation can be calculated from the innovation vector $v_k^i$ and its covariance $s_k^i$ as follows:

$${\mathsf{\Lambda }}_k^i=\frac{\exp \left\{-\left(\frac{1}{2}\right){\left(v_k^i\right)}^{\prime }{\left(s_k^i\right)}^{-1}{v}_k^i\right\}}{\sqrt{\left|2\pi s_k^i\right|}}(i=1,2,3)$$

(33)

where $v_k^i=Z_k-{\widehat{Z}}_{k,k-1}^i$, $s_k^i=P_{ZZ}^i$

Then, the model probability update is calculated as:

$${\mu }_k^i=\frac{{\mu }_{k,k-1}^i{\mathsf{\Lambda }}_k^i}{{{\displaystyle \sum }}_{j=1}^3{\mu }_{k,k-1}^j{\mathsf{\Lambda }}_k^j}(i=1,2,3)$$

(34)

(4) Estimation Fusion

Finally, the combined state $X_k$ can be calculated as:

$$X_k={{\displaystyle \sum }}_{i=1}^3{\mu }_k^i{X}_k^i$$

(35)

Since the proposed IMM-UKF can adapt to a wide variation of inertial sensor noise, the vehicle localization system is robust enough to achieve an accurate position output when facing uncertain inertial sensor noise.

4. Design the GNN Module

Considering the uncertain noise of MEMS inertial sensors and varied driving situations, it is very difficult to establish appropriate functions or equations to describe the dynamic behaviors of RISS position errors. The grey system theory requires only a limited amount of data to estimate the behavior of unknown systems. Through combing grey system theory with neural network, the predicting precision can be raised undoubtedly when the training samples are not sufficient. Therefore, the GNN module is developed here to predict the future position errors using the current available inertial sensor data.

For land vehicle applications, the horizontal localization performance is generally the main concern [7,40]. Thus, the latitude and longitude components of the difference between two state vectors associated with IMM-UKF1 and IMM-UKF2 are selected as the outputs. Since the vehicle maneuverability can affect the position errors [41], the longitudinal acceleration $a_k^x$ and the yaw rate ${\omega }_k^z$ measured by RISS are considered as the corresponding inputs. In actual implementation, two separate GNNs are designed in parallel for the position errors along latitude and longitude, respectively. It is worthwhile to mention here that the two GNNs have similar designing process. For simplicity, we choose the latitude component to show how to establish the GNN while the other one for the longitude component can be processed similarly. The GNN for the latitude component is developed as follows [42]:

(1) Construct the original data series:

$$\begin{array}{c}{z}_t^{\left(0\right)}=X_t^{IMM-UKF1}\left[1\right]-X_t^{IMM-UKF2}\left[1\right]\\ y_t^{1\left(0\right)}=a_t^x\\ y_t^{2\left(0\right)}={\omega }_t^z\end{array}t=1,2,\dots ,N$$

(36)

where $X_t^{IMM-UKF1}\left[1\right]$ and $X_t^{IMM-UKF2}\left[1\right]$ are the latitude component of the state vector of IMM-UKF1 and IMM-UKF2 at time step $t$, respectively. $N$ is the length of the sliding window and can be adjusted according to the length of the assumed GPS-outage time.

(2) Take accumulated generating operation (AGO) on $z_t^{\left(0\right)}$, $y_t^{1\left(0\right)}$, and $y_t^{2\left(0\right)}$, respectively. Then the AGO sequence can be obtained as:

$$\begin{array}{c}{Z}^{\left(1\right)}=\left(z_1^{\left(1\right)},z_2^{\left(1\right)},\dots ,z_N^{\left(1\right)}\right)\\ y^{1\left(1\right)}=\left(y_1^{1\left(1\right)},y_2^{1\left(1\right)},\dots \dots ,y_N^{1\left(1\right)}\right)\\ y^{2\left(1\right)}=\left(y_1^{2\left(1\right)},y_2^{2\left(1\right)},\dots ,y_N^{2\left(1\right)}\right)\end{array}$$

(37)

where $z_t^{\left(1\right)}={\displaystyle \sum }_{i=1}^t{z}_i^{\left(0\right)}$, $y_t^{1\left(1\right)}={\displaystyle \sum }_{i=1}^t{y}_i^{1\left(0\right)}$, $y_t^{2\left(1\right)}={\displaystyle \sum }_{i=1}^t{y}_i^{2\left(0\right)}$, $t=1,2,\dots ,N$.

(3) Form the whitening differential equation according to grey system theory:

$$\frac{d{z}_t^{\left(1\right)}}{dt}+b_1{z}_t^{\left(1\right)}=b_2{y}_t^{1\left(1\right)}+b_3{y}_t^{2\left(1\right)}$$

(38)

The solution of Equation (38) can be obtained as:

$${\widehat{z}}_t^{\left(1\right)}=\left(z_1^{\left(1\right)}-\frac{b_2}{b_1}{y}_t^{1\left(1\right)}-\frac{b_3}{b_1}{y}_t^{2\left(1\right)}\right)e^{-b_{1}t}+\frac{b_2}{b_1}{y}_t^{1\left(1\right)}+\frac{b_3}{b_1}{y}_t^{2\left(1\right)}$$

(39)

Let $d=\frac{b_2}{b_1}{y}_t^{1\left(1\right)}+\frac{b_3}{b_1}{y}_t^{2\left(1\right)}$, and Equation (39) can be transformed to:

$${\widehat{z}}_t^{\left(1\right)}=\left(\left(z_1^{\left(1\right)}-d\right)·\frac{e^{-b_{1}t}}{1+e^{-b_{1}t}}+d·\frac{1}{1+e^{-b_{1}t}}\right)·\left(1+e^{-b_{1}t}\right)=\left(\left(z_1^{\left(1\right)}-d\right)-z_1^{\left(1\right)}·\frac{1}{1+e^{-b_{1}t}}+2d·\frac{1}{1+e^{-b_{1}t}}\right)·\left(1+e^{-b_{1}t}\right)$$

(40)

(4) Map Equation (40) to an expanded BPNN. Thus we can obtain the GNN with two input variables and one output variable, as shown in Figure 4.

Here, t is the time step and can also be treated as one hidden input. ω₁₁, ω₂₁, ω₂₂, ω₂₃, ω₃₁, ω₃₂, ω₃₃ are network weighting values. LA, LB, LC, LD are the four layers of GNN respectively. The learning process of GNN is as follows:

Step 1: Initialize the network parameters and weighting values.

Let $\frac{2b_2}{b_1}=u_1$ and $\frac{2b_3}{b_1}=u_2$, then the network initial weighting value can be represented as: ${\omega }_{11}=b_1$, ${\omega }_{21}=-z_1^{\left(1\right)}$, ${\omega }_{21}=u_1$, ${\omega }_{23}=u_2$, ${\omega }_{31}={\omega }_{32}={\omega }_{33}=1+e^{-b_{1}t}$

Step 2: Calculate the output of each layer at each time step.

• LA layer: $o_a={\omega }_{11}t$
• LB layer: $o_b=\frac{1}{1+e^{-{\omega }_{11}t}}$
• LC layer: $o_{c1}=o_b{\omega }_{21}$, $o_{c2}=y_t^{1\left(1\right)}{o}_b{\omega }_{22}$, $o_{c3}=y_t^{2\left(1\right)}{o}_b{\omega }_{23}$
• LD layer: $o_d=o_{c1}{\omega }_{31}+o_{c2}{\omega }_{32}+o_{c1}{\omega }_{33}-\theta$

where $\theta$ is the threshold value and can be calculated as $\theta =(1-e^{-b_{1}t})(d-z_1^{\left(1\right)})$.

Step 3: Calculate the errors between the forecast and expectation, and then adjust the weighting values and the threshold value.

The error of each level can be calculated as:

LD layer error: ${\delta }_d=o_d-z_t^{\left(1\right)}$

LC layer error: ${\delta }_{c1}={\delta }_{c2}={\delta }_{c3}={\delta }_d\left(1+e^{-{\omega }_{11}t}\right)$

LB layer error: ${\delta }_b=\frac{1}{1+e^{-{\omega }_{11}t}}\left(1-\frac{1}{1+e^{-{\omega }_{11}t}}\right)\left({\omega }_{21}{\delta }_{c1}+{\omega }_{22}{\delta }_{c2}+{\omega }_{23}{\delta }_{c3}\right)$

The weighting values can be adjusted as:

${\omega }_{21}=-z_1^{\left(1\right)}$, ${\omega }_{22}={\omega }_{22}-{\mu }_1{\delta }_{c2}{o}_b$, ${\omega }_{23}={\omega }_{23}-{\mu }_2{\delta }_{c3}{o}_b$, ${\omega }_{11}={\omega }_{11}+o_{a}t{\delta }_b$

where ${\mu }_1$ and ${\mu }_2$ are learning rates, which are defined previously.

The threshold value can be adjusted as:

$$\mathsf{\theta }=\left(1+e^{-{\omega }_{11}t}\right)\left(\frac{{\omega }_{22}}{2}{y}_t^{1\left(1\right)}+\frac{{\omega }_{23}}{2}{y}_t^{2\left(1\right)}-z_1^{\left(1\right)}\right)$$

Step 4: Back to Step 2, re-adjust the weighting values and the threshold value until GNN is convergent.

GNN has a rapid convergence rate. Usually, the optimal weighting values and threshold value can be achieved after adjusted twice. After the GNN is convergent, it can be utilized to efficiently predict the corresponding position error, which is calculated by taking the inverse accumulated generation operation (IAGO) operation on ${\widehat{z}}_t^{\left(1\right)}$.

$${\widehat{z}}_1^{\left(0\right)}={\widehat{z}}_1^{\left(1\right)}$$

$${\widehat{z}}_t^{\left(0\right)}={\widehat{z}}_t^{\left(1\right)}-{\widehat{z}}_{t-1}^{\left(1\right)},t=2,\dots ,N$$

5. Experiments and Results

5.1. Equipment and Road Trajectories

To evaluate the localization performance of the proposed solution, several experiments were conducted on a Chery TIGGO5 SUV (Chery Automobile Co., Ltd., Wuhu, China). Since the vehicle was equipped with ABS and ESP, the information about steering angle and forward speed could be directly obtained from the in-vehicle CAN bus. Besides, a low-cost NovAtel Superstar II GPS receiver (NovAtel, Calgary, AL, Canada) with 1 Hz rate and MEMSIC MEMS-based IMU VG440CA-200 sampled at 100 Hz were installed. The RISS data used in this research is from the one vertical gyroscope and two horizontal accelerometers of the full six-degree-of-freedom (6-DoF) IMU VG440CA-200. For the MEMS-based inertial sensors, the gyroscope has a bias stability of 10°/h and angle random walk of 4.5°/√h, while each accelerometer has bias stability of 1 mg and velocity random walk of 1 m/s/√h. The accuracies of other sensors (1σ) are 0.05 m/s and 3 m for the GPS velocity and position, 0.05 m/s for the wheel speed sensor, and 4° for the steering angle, respectively. Moreover, an accurate and reliable NovAtel SPAN-CPT system was used as a reference for quantitative comparison. The horizontal position accuracy of SPAN-CPT system was 0.01 m in absence of GPS outages and 0.02 m during 10 s outage.

Several road-test experiments were carried out along different trajectories using the setup described above. One of the trajectories was on the Fifth Ring Road in Beijing, which was a typical urban scenario with real GPS-denied environments in some parts. Besides, a series of typical driving maneuvers, such as lane-changes, accelerations and decelerations etc., were conducted according to actual driving conditions. It is worthwhile to mention here that, in this paper, that the position errors denote the horizontal Euclidean distance error between the estimated position and the corresponding reference, which is the main concern for land vehicle localization.

5.2. Test 1: Performance Evaluation of the Proposed Localization Solution in Trajectories 1

The trajectory was shown in Figure 5. Straight portions and curves were considered when selecting outages in this trajectory. Since some periods of real GPS outages were shorter than 45 s, the selected outages were all extended to 45 s for convenient comparison.

In this test, the overall performance of the proposed localization solution was evaluated. As shown in Figure 1 and Figure 2, our proposed solution fuses the information from GPS, MEMS-RISS, and in-vehicle sensors utilizing IMM-UKF and compensates position errors utilizing GNN during GPS outages. Thus, the proposed solution is termed as IMM-UKF-GNN. In order to highlight the advantages of our proposed methodology, three other methods are also conducted for comparison: (1) General UKF without any compensation during GPS outages, termed as UKF; (2) IMM-UKF without any compensation during GPS outages, termed as IMM-UKF; (3) IMM-UKF with Radial Basis Function (RBF) compensation during GPS outages, termed as IMM-UKF-RBF. Note that both general UKF and IMM-UKF have the same motion model and measurement model described in Section 3. The difference is that the UKF method only has a constant Q matrix. In other words, the general UKF can be treated as one of the UKFs in the IMM-UKF. Both the UKF method and the IMM-UKF method can only execute the measurement update associated with in-vehicle sensors during GPS outages and without any further compensation. Since RBF has been widely regarded as the most remarkable ANN during the past decades [28], it is elected to compare with GNN in this paper. The RBF module was designed with the same inputs and outputs as GNN. Besides, the same 45 s sliding window was also utilized to train the RBF module. The learning procedures of both GNN and RBF continue as long as the GPS signal is available. In case of GPS outages, the trained RBF and GNN module are utilized to predict and compensate the position errors. In the absence of GPS outages, all the four methods can provide an accurate position output. Therefore, we focus on the comparisons of the performances among the four methods during GPS outages.

Table 1. Maximum position errors during GPS outages.

Outage Number	Maximum Error (m)
	UKF	IMM-UKF	IMM-UKF-RBF	IMM-UKF-GNN
1	26.32	24.62	5.24	3.91
2	55.20	46.94	16.21	15.82
3	64.62	56.76	17.92	10.29
4	53.64	51.74	17.20	9.28
5	59.98	51.76	12.28	7.83
6	26.13	23.83	11.66	8.21

and

Table 2. RMS position errors during GPS outages after inserting biases.

Outage Number	RMS Error (m)
	UKF	IMM-UKF	IMM-UKF-RBF	IMM-UKF-GNN
1	6.31	5.42	1.22	1.18
2	13.84	12.96	7.69	7.41
3	21.41	16.43	4.72	3.02
4	16.54	14.34	5.86	2.84
5	18.41	14.23	2.81	1.79
6	7.95	7.81	3.71	2.31

give a quantitative comparison of the maximum and RMS position errors among the four methods described above during the six GPS outages, respectively. The highlighted columns correspond to the least errors achieved by the proposed localization solution. From the tables, it can be determined that the IMM-UKF method outperforms the UKF method. When GPS outages occur, the position errors will accumulate rapidly due to the uncertain noise of MEMS inertial sensors, the IMM-UKF can adapt to the uncertain noise and thus mitigating the error accumulation. On average, the IMM-UKF method achieves 9.9% and 13.6% improvements on maximum error and RMS error over the UKF method, respectively. However, both methods cannot ensure the localization accuracy and reliability during GPS outages.

The localization results also show that the methods with compensation can achieve much smaller errors than those without. Since the RBF and GNN can mimic the latest position errors, these errors can be removed from corresponding position components and thereby improving the localization accuracy during GPS outages. Furthermore, due to the advantages of GNN with respect to insufficient modeling information, the maximum error of the proposed IMM-UKF-GNN solution is 28.5% lower than that of the IMM-UKF-RBF method on average. When it comes to the RMS error, the proposed solution achieves a 28.1% lower value than IMM-UKF-RBF method.

In order to directly show the localization results of different methods, three representative outages, i.e., outages 1, 3, and 4, were chosen to show the trajectories. Outages 1 and 4 correspond to the portion of the trajectory for the vehicle moving along curves, illustrated in Figure 6 and Figure 7. During outage 1, the percentage improvement of the proposed solution in maximum error is found to improve by 25.4%, 84.1%, and 85.1% against IMM-UKF-RBF method, IMM-UKF method and UKF method respectively, while the percentage improvement in RMS error is found to improve by 3.3%, 78.2%, and 81.5% respectively. For outage 4, the proposed solution effectively reduces the maximum error by 46.0%, 82.1%, and 82.7% against IMM-UKF-RBF method, IMM-UKF method, and UKF method, and reduces the RMS error by 51.5%, 80.2%, and 82.8%, respectively.

Figure 8 gives the localization results for outage 3, which belongs to a typical straight road. The percentage improvement of the proposed solution in maximum error is found to improve by 42.6%, 81.9%, and 84.1% against IMM-UKF-RBF method, IMM-UKF method and UKF method respectively and the percentage improvement in RMS error is 36.0%, 81.6%, and 85.9% respectively.

5.3. Test 2: Further Evaluation of the Proposed Localization Solution

In order to further test the adaption to uncertain noise of the proposed solution, we inserted biases into the inertial sensor data during the periods of GPS outages. The biases are modeled by the first-order Gauss-Markov process. The correlation time is defined as 100 s and the standard deviation of the white noise associated with the process is 10 mg for the accelerometers and 100°/h for the gyroscope. After inserting biases, the statistical properties of inertial sensor errors were dramatically changed during simulated outages, and the Q matrix should be updated correspondingly, which was not capable for the general UKF. Thus, the inaccurate Q matrix would cause performance degradation in the UKF method. However, the proposed IMM-UKF was envisioned to be adaptive to the inserted biases.

Table 3. Maximum position errors during GPS outages after inserting biases.

Outage Number	Maximum Error (m)
	UKF	IMM-UKF	IMM-UKF-RBF	IMM-UKF-GNN
1	31.79	25.77	6.84	4.35
2	62.37	47.23	17.38	16.13
3	71.31	57.20	21.52	11.77
4	61.17	53.63	21.62	12.85
5	67.44	54.14	14.16	8.15
6	32.62	25.39	13.98	9.71

and

Table 4. RMS position errors during GPS outages after inserting biases.

Outage Number	RMS Error (m)
	UKF	IMM-UKF	IMM-UKF-RBF	IMM-UKF-GNN
1	7.36	6.03	1.84	1.45
2	15.58	13.21	8.26	7.73
3	22.47	17.19	5.45	3.53
4	17.28	14.57	6.34	3.22
5	20.02	15.09	2.66	1.66
6	9.33	8.10	4.38	2.89

show the results of maximum and RMS position errors among the four methods during the six GPS outages after inserting biases. The proposed IMM-UKF-GNN solution can still achieve the best maximum and RMS position errors.

Comparing Table 1 and Table 3, it can be determined that, after the statistical properties of inertial sensor errors are intentionally changed, the maximum error of the methods with IMM-UKF only increases 1.8 m on average while the increase of the UKF method is 6.8 m. The increase of the maximum position errors among the four methods is also shown in Figure 9. Besides, comparing Table 2 and Table 4 it can be found that the increase of RMS error is 0.47 m on average for the methods with IMM-UKF, while the increase is 1.26 m for the UKF method. The increase of the maximum position errors among the four methods is also depicted in Figure 10. Thus, it can be concluded that, when facing the same situation of uncertain inertial sensor noise, the IMM-UKF can achieved better performance than the general UKF.

6. Conclusions

This paper has presented a cost-effective vehicle localization solution, which can simultaneously address uncertain noises of MEMS inertial sensors and GPS outages. The proposed IMM-UKF fuses information from low-cost GPS, MEMS-RISS, and in-vehicle sensors. Three UKFs with different covariances are developed to cover a wide variation of inertial sensor noise. Then, an accurate estimation of vehicle positions can be obtained when GPS is available. Meanwhile, another IMM-UKF is developed to execute the measurement update associated with in-vehicle sensors. The difference between the state vector of the two IMM-UKFs are modeled by a GNN module. When GPS outages occur, the latest updated GNN module can predict and compensate position errors. Thus, the proposed solution can achieve accurate localization even without GPS observations.

The proposed localization solution has been successfully implemented and tested with real road-test trajectories. Through comparison with other three representative localization methods, it can be concluded that the research fulfills the basic aim of proposing a cost-effective vehicle localization solution which can maintain relatively good performance when facing uncertain inertial sensor noises and GPS outages simultaneously.