A Redundant Measurement-Based Maximum Correntropy Extended Kalman Filter for the Noise Covariance Estimation in INS/GNSS Integration

Abstract

The resolution accuracy of the inertial navigation system/global navigation satellite system (INS/GNSS) integrated system would be degraded in challenging areas. This paper proposed a novel algorithm, which combines the second-order mutual difference method with the maximum correntropy criteria extended Kalman filter to address the following problems (1) the GNSS measurement noise estimation cannot be isolated from the state estimation and suffers from the auto-correlated statistic sequences, and (2) the performance of EKF would be degraded under the non-Gaussian condition. In detail, the proposed algorithm determines the possible distribution of the measurement noise by a kernel density function detection, then depending on the detection result, either the difference sequences–based method or an autoregressive correction algorithm’s result is utilized for calculating the noise covariance. Then, the obtained measurement noise covariance is used in MCEKF instead of EKF to enhance filter adaptiveness. Meanwhile, to enhance the numerical stability of the MCEKF, we adopted the Cholesky decomposition to calculate the matrix inverse in the kernel function. The road experiment verified that our proposed method could achieve more accurate navigation resolutions than the compared ones.

1. Introduction

To achieve satisfactory navigation accuracy in real-world applications, multiple subsystems are often used instead of a single system. The inertial navigation system (INS) is usually an essential component in the navigation system due to its self-contained nature. The global navigation satellite system (GNSS) can be installed easily, and its result error does not accumulate over time. Thus, an effective combination of the two systems will produce more satisfactory outcomes [1]. Specifically, the advantages of the INS include high data output frequency, strong anti-interference, and low short-time noise [2]. However, the main drawback of the INS is the output accuracy could not maintain a high level over time. The GNSS could maintain long-term stability under an ideal environment, but the data bandwidth is much lower than that of the INS. Moreover, the navigation solution would be degraded if the GNSS is in a complex environment, such as in an urban street, tunnel, dense building area, and so on [3]. Therefore, the aim of INS/GNSS is to conquer their individual shortcomings, i.e., to suppress the prolonged drift of the INS and to filter the solution of GNSS to a more stable and accurate one [4,5]. Xia et al. proposed a method for the IMU and automotive onboard sensors to estimate the yaw misalignment autonomously. Based on this method, the estimation accuracy can be improved hugely when the INS is integrated with other systems [6]. Gao et al. proposed an improved low-cost vehicle localization method based on onboard sensors, including IMU, wheel speed sensors, and steering wheel angle sensor. This method could keep a sub-meter-level during about 20 s GNSS outages and can be seen as an excellent approach [7].

To obtain the navigation solution, the Kalman filter (KF) and the related extensions are frequently employed in the integrated navigation system [8,9,10,11]. For linear problems with the assumption that the noise both in the process and in measurement follow the Normal distribution, by utilizing the minimum mean square error (MMSE) criterion, KF could obtain an optimal estimation [12]. For nonlinear problems, several effective modified filters have been most commonly employed since the 1970s [13].

The extended Kalman filter (EKF) adopts a straightforward linearization method, which locally linearizes the nonlinear system state transition equation or observation equation through the first-order Taylor expansion (TE) [14]. Different from the EKF, to get a higher-order approximation, the unscented transformation (UT) is employed by the unscented Kalman filter (UKF) [15]. For high-dimensional problems, the performance of the UKF would degrade. Thus, a method based on the cubature integration rules with Gaussian weighted factors is introduced by the cubature Kalman filter (CKF) [16]. PF and the Gaussian sums filter (GSF) could give a better estimation than others but at the cost of a higher computational burden [12]. In INS/GNSS integrated system, EKF is usually adopted as the fusion algorithm to balance the estimation accuracy and the computation load [13].

However, if the measurement noise or the process noise contains non-Gaussian components, the performance of the above filters would be severely degraded [17]. For example, the Middleton Class-A model, Middleton Class-B model, alpha-stable distribution model, and the mixture Gaussian model are frequently discussed non-Gaussian noises [18]. The Middleton Class-A noise describes the type of electromagnetic interference (EMI) often encountered in telecommunication applications, and the Middleton Class-B one usually denotes the man-made or nonmessage-bearing noise [19]. Ndo et al. proposed a method to deal with the problem that the transmission over channels is impaired by impulsive noise [20]. Based on the same parameters used in the Middleton Class-A model, they added a correlation parameter that enables reproducing the bursty nature of impulses, and the built optimum detector could obtain significant performance gains compared to the Middleton Class-A model. The alpha-stable model and the Gaussian mixture model can be seen as the parametric method. Peng et al. proposed a mixture of Alpha-stable (MAS) distributions for modeling the statistical property of Synthetic Aperture Radar (SAR) images, and the effectiveness is verified by the experiments on TerraSAR-X images for modeling and classification of SAR images [21]. Zhang et al. proposed an individual tree segmentation (ITS) optimization method based on the Gaussian mixture model for airborne LiDAR data. The experiment was conducted with two datasets of airborne LiDAR data, and the results showed that the proposed method well improved the under-segmentation in the ITS of high-density, multistoried mixed forests, and achieved a better segmentation accuracy [22].

In practical INS/GNSS applications, impulsive-type noise in the GNSS signal is relatively typical, which makes the measurement noise tend to be heavy-tailed [23]. For example, in the following environments, the GNSS outputs are particularly susceptible to interference: foliage-covered districts, urban overpasses, and tunnels [24]. To get a better navigation solution accuracy, the noise covariance of the measurement noise should be estimated timely. Therefore, many studies have proposed solutions aiming to suppress the heavy-tailed noise in GNSS signals. Therefore, based on the fact that the MMSE is no longer suitable for this type of noise, the corresponding robust adaptive filtering methods have been proposed to address the issue by introducing different metrics, such as the Huber-based algorithms [25]. However, these filters can achieve desirable results only under certain conditions [17]. The maximum correntropy Kalman filter (MCKF) can be employed to address non-Gaussian noise problems, which was inspired by the information learning theorem (ILT) [26,27]. Additionally, to employ the above-mentioned non-Gaussian models, many parameters should be determined either by training offline or based on the knowledge of the statistics of the noise. Fortunately, the MCC based on the kernel technique conquers this parametric problem well.

The MCKF is a combination of MCC and KF. Due to the ability of MCC to capture higher-order statistics of the noise, it can achieve a higher estimation accuracy than other filters if non-Gaussian noise exists. For different applications, to get a more accurate estimation, scholars proposed modified versions of MCKF. Qu et al. proposed a generalized maximum asymmetric correntropy criterion. Different from the classical Gaussian kernel function, they adopted a generalized asymmetric Gaussian density function. The proposed method could achieve more accurate prediction than the classical MCC in a delay prediction simulation of a satellite network [28]. Liao et al. proposed a dynamic self-tuning MCKF (DSTMCKF), which was used to improve positioning accuracy in wireless sensor networks (WSN). This algorithm employed the observation as a priori information together with the innovation to calculate the noise covariance effectively; then, the new noise covariance estimation could be taken into a standard MCKF procedure [29]. Wang et al. introduced the MCC to the extended Kalman filter (MCEKF) under the condition that measurement outliers exist for nonlinear systems. They combined the MCC and the weighted least squares (WLS) algorithm in EKF. The MCEKF performed better than the standard EKF and UKF in an object trajectory estimation simulation [30].

Usually, in an INS/GNSS integrated system, the process noise covariance matrix $Q$ can be seen as the uncertainties in the dynamic model, while the measurement noise covariance matrix $R$ can be regarded as an indicator of the accuracy of the observations. $Q$ and $R$ are set to predetermined values in many realistic applications [31]. However, in practice, they could not be constant, if a more accurate state estimation is desired, one needs to estimate $Q$ and $R$ as accurately and timely as possible [21]. Although the kernel function technique-based MCEKF already takes process noise and measurement noise into consideration, it can get better estimation accuracy when $R$ could be accurately obtained. Many adaptive algorithms have been put forward to estimate $R$. Ref. [32] proposed the Sage-Husa statistics estimator for obtaining a more accurate state estimation. Similar methods also include the strong tracking Kalman filter [33,34] and the $H\infty$ filter [35], to name a few. To sum up, these algorithms employ the sequences that contains the state vector to evaluate $R$.

However, the sequences of innovation and residual are not independent, and the state estimation error will influence the estimation of $R$. In order to address this problem, a few researchers proposed algorithms that only employ the sensors’ data to estimate the covariance matrix in a system with redundant measurements. Ref. [36] proposed a random weighting method (RWM) to fuse the data in a multisensor system. This method could achieve variances of an individual sensor; therefore, a good estimation accuracy could be obtained through optimal weighting distribution. However, the RWM only uses the second-order moments of the raw sensors’ data; the first-order information is not fully used. Based on the RWM idea, a randomly weighted CKF algorithm is developed [37], and an adaptive unscented Kalman filter is proposed, which is based on the maximum posterior and random weighting method [38]. However, these algorithms adjust the weights basically on the state-contained sequences simultaneously, except for the RWM. Thus, the drawback that the error in state estimation would influence the evaluation of the measurement noise covariance inevitably exists in these algorithms [39]. Velázquez proposed a variance estimation algorithm that could be implemented online. By using a Minimum Norm Quadratic Unbiased Estimation [40], Velázquez argued that the proposed algorithm could estimate the individual variances at each time instant. However, the prerequisite for applying this algorithm is that the number of observers must be more than twice the amount of signals to be measured [39].

To address the problems for INS/GNSS in unfriendly areas: (1) the method for estimating $R$ usually involves the estimation of the state; hence the estimation result would be contaminated if the state estimation is not accurate; and (2) the performance of the MCEKF would be degraded if non-Gaussian noise exists, we proposed a Redundant measurement-based Maximum Correntropy Extended Kalman Filter (RMCEKF). The key contributions are listed as follows:

• To make full use of the first- and second-order measurement statistics, the proposed method adopts the enhanced SOMD algorithm to estimate $R$ in real time. The RMCEKF is the first algorithm that combines the estimation of $R$ with observations only. Moreover, to get better numerical stability for the MCEKF, the matrix inverse operation is calculated by the Cholesky decomposition.

• To verify the effectiveness and performance of the RMCEKF, real road experiment data were used. The result showed that RMCEKF could achieve better estimation accuracy than the compared ones.

The ensuing sections of this paper are arranged in the subsequent manner. The INS/GNSS tightly coupled integration and the standard EKF solution procedure are reviewed in Section 2. In Section 3, the MCEKF and the SOMD are described. Meanwhile, the drawback of the MCEKF is analyzed. Then, in Section 4, the novel RMCEKF will be given. A real road test and the precision comparisons of different algorithms are given in Section 5. The conclusions are given in Section 6.

2. INS/GNSS Tightly Coupled Integration and EKF

There are generally three groups for INS/GNSS integration systems [41]. The difference is the way how the GNSS signals are employed in the integration procedure. In the loosely coupled architecture, the GNSS data are used as the measurements of KF. Although redundancy and simplicity could be achieved, the disadvantage is that the requirement of visibility of at least four satellites makes the loosely coupled method only suitable to some ideal environments. Xu et al. constructed a loosely coupled system by employing the BeiDou global navigation satellite system (BDS-3). Through two experiments with different INS in a complex urban environment for 15 min, they verified their method could improve the positioning accuracy by over 60% in every direction [42]. Zhai et al. proposed a vision-aided loosely coupled navigation system based on semantic information. They verified the algorithm using the KITTI dataset and argued that the loosely coupled navigation system could achieve satisfying estimation accuracy in GNSS-challenged environments [43]. Chiang et al. proposed a method based on the Non-Holonomic Constraint (NHC) to improve the INS/GPS vehicle navigation solution accuracy. They applied the algorithm with both loosely coupled and tightly coupled schemes to test the performance. One field test result showed that in poisoning, compared to the loosely coupled strategy, the tightly coupled and the proposed method with abnormal GPS measurement rejection could achieve an improvement of 33% and 61%, respectively [44].

In a tightly coupled integration system, the pseudo-range and pseudo-range rate of GNSS is used as elements in KF’s state vector. This algorithm is widely adopted for its anti-multipath ability. For the ultra-tightly coupled mechanism, the INS is utilized to aim the GNSS receiver for the loop track. Thus, INS can be seen as a part of the GNSS receiver. This architecture can enhance the anti-jamming ability. However, this architecture is complex because it requires intensive computation resources and hardware redesign [45]. Thus, in this study, we adopt the tightly coupled mechanism as the navigation solution method. As the loosely coupled scheme is widely used in real applications, such as in a low-cost scenario, our proposed algorithm could be deployed by modifying the measurement function similar to [44].

2.1. Dynamic Model

We adopt the kinetic equation mechanism for the navigation system discussed in [46].

The expression for the propagation of attitude error in the navigation frame is as follows:

(1)${\dot{\mathsf{\phi }}}^n=-\left({\omega }_{ie}^n+{\omega }_{en}^n\right)\times {\mathsf{\phi }}^n+\left(\delta {\omega }_{ie}^n+\delta {\omega }_{en}^n\right)-C_b^n{\epsilon }^b$

where the superscript $n$ and $b$ denote the navigation frame and the body frame, respectively, the subscript $e$ represents the earth-centered, earth-fixed frame, ${\mathsf{\phi }}^n$ is the attitude error, ${\omega }_{ie}^n$ denotes the earth rate, ${\omega }_{en}^n$ is the rate of $n$-frame with respect to $e$-frame, $C_b^n$ represents the direction cosine matrix transmitted from the $b$-frame to the $n$-frame, and ${\epsilon }^b$ is the gyroscope drift.

The position error can be given as:

(2)$\left[\begin{array}{c}\delta \dot{L}\\ \delta \dot{\lambda }\\ \delta \dot{h}\end{array}\right]=\left[\begin{array}{c}\frac{\delta V_y^n}{R_M+h}-\frac{V_y^n}{{\left(R_M+h\right)}^2}\delta h\\ \frac{\delta V_x^n}{\left(R_N+h\right)\cos L}-\frac{V_x^n}{{\left(R_N+h\right)}^2\cos L}\delta h+\frac{V_x^n\cdot \tan L}{\left(R_N+h\right)\cos L}\delta L\\ V_z^n\end{array}\right]$

where $[\delta L\delta \lambda \delta h]$ stands for the position error vector, which includes the error in latitude, longitude, and height are listed. $V_x^n$ denotes the velocity in the east direction, $V_y^n$ in the north direction, and $V_z^n$ in the up direction. The meridian and transverse radius of curvature are denoted by $R_M$ and $R_N$, respectively.

The velocity error transition function is derived as:

(3)$\delta {\dot{\mathit{V}}}^n=-f^n\times f^n+\delta {\mathit{V}}^n\times \left(2{\omega }_{ie}^n+{\mathbf{\omega }}_{en}^n\right)+V^n\times \left(2\delta {\mathbf{\omega }}_{ie}^n+\delta {\mathbf{\omega }}_{en}^n\right)+\delta g^n+C_b^n{\mathsf{\nabla }}^b$

in which the velocity error is denoted by $\delta {\mathit{V}}^n$, $f^n$ denotes the specific force error, $\delta g^n$ is the gravitational error, and ${\mathsf{\nabla }}^b$ represents the accelerometer error.

The gyroscope error model is expressed by:

(4)$\left\{\begin{array}{l}{G}_{\epsilon }=G_{\epsilon }{}_b+G_{\epsilon }{}_r\\ {\dot{G}}_{\epsilon }=0\\ {\dot{G}}_{\epsilon r}=-\frac{1}{C_{tg}}{G}_{\epsilon r}+{\omega }_g\end{array}\right.$

in which $G_{\epsilon }{}_b$ is the gyroscope bias, $G_{\epsilon }{}_r$ stands for the first-order Gauss-Markov (GM) process, $C_{tg}$ denotes coefficient of the correlation time, and ${\omega }_g$ represents the noise.

It can be assumed to be the first-order Markov process for the model of the accelerometer too. Additionally, we assume that the three axes follow the same model [46]:

(5)${\dot{\mathsf{\nabla }}}^b=-\frac{1}{T_a}{\mathsf{\nabla }}^b+{\omega }_a$

where $T_a$ represents correlation time coefficient and ${\omega }_a$ is assumed as the Gaussian white noise.

The stochastic model of $\delta t_u$ and $\delta t_{tu}$ in this paper are also considered as the first-order Markov process:

(6)$\begin{array}{l}\delta {\dot{t}}_u=\delta t_{ru}+{\omega }_u\\ \delta {\dot{t}}_{ru}=-\beta \delta t_{ru}+{\omega }_{ru}\end{array}$

where $\delta t_u$ denotes the GNSS receiver’s clock offset, $\delta t_{ru}$ is the drift of $\delta t_u$, ${\omega }_u$ and ${\omega }_{tu}$ represent the Gaussian white noise item, and $-\beta$ stands for the correlation time factor.

2.2. Measurement Model

The observation function considered in the system can be given by the difference between the pseudo-range and pseudo-range rate calculated by INS and the pseudo-range and pseudo-range rate output by the GNSS receiver. First, the pseudo-range of INS-driven and GNSS together with the difference between them are given as [3]:

(7)$\left\{\begin{array}{l}{\rho }_I^j=\sqrt{{(x_I-x_{sa}^j)}^2+{(y_I-y_{sa}^j)}^2+{(z_I-z_{sa}^j)}^2}(j=1,2,3,4)\hfill \\ {\rho }_G^j=\sqrt{{\left(x-x_{sa}^j\right)}^2+{\left(y-y_{sa}^j\right)}^2+{\left(z-z_{sa}^j\right)}^2}-\delta t_u-v_{\rho }^j(j=1,2,3,4)\hfill \\ \delta {\rho }_j={\rho }_I^j-{\rho }_G^j=e_{j1}\delta x+e_{j2}\delta y+e_{j3}\delta z+\delta t_u+v_{\rho }^j\hfill \end{array}\right.$

in which $[x_I,y_I,z_I]$ represents the coordinate of INS in the $e$-frame. For the GNSS receiver, $[x,y,z]$ is true position in the $e$-frame, $[x_{sa}^j,y_{sa}^j,z_{sa}^j]$ denotes the position of the $j\mathrm{th}$ satellite, $v_{\rho }^j$ is the measurement noise of the $j\mathrm{th}$ satellite.

Subsequently, the difference in pseudo-range rate between the two systems is expressed as follows:

(8)$\left\{\begin{array}{l}\begin{array}{l}{r}_j=\sqrt{{\left(x-x_{sa}^j\right)}^2+{\left(y-y_{sa}^j\right)}^2+{\left(z-z_{sa}^j\right)}^2}\\ {\dot{\rho }}_I^j=e_{j1}({\dot{x}}_I-{\dot{x}}_{sa}^j)+e_{j2}({\dot{y}}_I-{\dot{y}}_{sa}^j)+e_{j3}({\dot{z}}_I-{\dot{z}}_{sa}^j)+e_{j1}\delta \dot{x}+e_{j2}\delta \dot{y}+e_{j3}\delta \dot{z}\end{array}\hfill \\ {\dot{\rho }}_G^j=e_{j1}(\dot{x}-{\dot{x}}_{sa}^j)+e_{j2}(\dot{y}-{\dot{y}}_{sa}^j)+e_{j3}(\dot{z}-{\dot{z}}_{sa}^j)+\delta t_{ru}+v_{\dot{\rho }}^j\hfill \\ \delta {\dot{\rho }}_j={\dot{\rho }}_I^j-{\dot{\rho }}_G^j=e_{j1}\delta \dot{x}+e_{j2}\delta \dot{y}+e_{j3}\delta \dot{z}-\delta t_{ru}-v_{\dot{\rho }}^j\hfill \end{array}\right.$

where $v_{\dot{\rho }}^j$ is the measurement noise, and $e_j$ can be given by:

(9)$\left\{\begin{array}{c}\frac{\partial {\rho }^j}{\partial x}=\frac{x-x_s^j}{r_j}=e_{j1}\\ \frac{\partial {\rho }^j}{\partial y}=\frac{y-y_s^j}{r_j}=e_{j2}\\ \frac{\partial {\rho }^j}{\partial z}=\frac{z-z_s^j}{r_j}=e_{j3}\end{array}\right.$

In summary, we can draw the measurement function as follows:

(10)$\left\{\begin{array}{l}\delta {\rho }_j={\rho }_I^j-{\rho }_G^j\hfill \\ \delta {\dot{\rho }}_j={\dot{\rho }}_I^j-{\dot{\rho }}_G^j\hfill \end{array}\right.$

2.3. Navigation Filtering with the Standard EKF

An EKF algorithm is applied for the integrated system to obtain the navigation solution. The model can be expressed as:

(11)$\left\{\begin{array}{l}{x}_k=A_{k,k-1}{x}_{k-1}+{{\rm M}}_k{v}_k^P\\ z_k={{\rm B}}_k{x}_k+v_k^R\end{array}\right.$

in which $x={\left[{\mathsf{\phi }}^n,\delta {\mathit{V}}^n,\delta L,\delta \lambda ,\delta h,{\epsilon }^b,{\mathsf{\nabla }}^b,\delta t_u,\delta t_{ru}\right]}^T\in {ℝ}^{17}$, $A_{k,k-1}$ is the state transition matrix, and $v_k^P$ denotes the process noise. Both $A_{k,k-1}$ and $v_k^P$ can be calculated by stacking (1)–(6). ${{\rm M}}_k$ represents the noise-driving matrix. $z=\left[\delta {\rho }_j,\delta {\dot{\rho }}_j\right]$, ${{\rm B}}_k$ denotes the measurement matrix, also can be linearized by the first-order TE, i.e., ${{\rm B}}_k=\partial b/\partial {\widehat{x}}_{k/k-1}$, $b$ is the nonlinear measuring function obtained from (7)–(10), $v_k^R$ is the measurement noise. The detail form of $A_{k,k-1}$ and ${{\rm B}}_k$ can be found in [45]. The statistics of $v_k^P$ and $v_k^R$ are assumed:

(12)$\left\{\begin{array}{l}E[v_t^P{v}_{\tau }^P]=Q_k{\delta }_{t,\tau }\\ E[v_t^R{v}_{\tau }^R]=R_k{\delta }_{t,\tau }\end{array}\right.$

in which ${\delta }_{t,\tau }$ stands for the Kronecker delta function.

Typically, a two-step process is used in an EKF:

Step Prediction:

(13)$\left\{\begin{array}{l}{\widehat{x}}_{k/k-1}=A_{k,k-1}{\widehat{x}}_{k-1}\hfill \\ P_{k/k-1}=A_{k,k-1}{P}_{k-1}{A}_{k,k-1}^T+{{\rm M}}_k{Q}_{k-1}{{\rm M}}_k^T\hfill \end{array}\right.$

in which ${\widehat{x}}_{k-1}$ is the estimated state and $P_{k-1}$ denotes the estimation error covariance at the time instant $k-1$. ${\widehat{x}}_{k/k-1}$ and $P_{k/k-1}$ represent the prediction of ${\widehat{x}}_{k-1}$ and $P_{k-1}$, respectively. $Q_{k-1}$ stands for the process noise covariance.

Step Measurement correction (updating):

(14)$\left\{\begin{array}{l}{K}_k=P_{k/k-1}{{\rm B}}_k^T{\left({{\rm B}}_k{P}_{k/k-1}{{\rm B}}_k^T+R_k\right)}^{-1}\\ {\widehat{x}}_k={\widehat{x}}_{k/k-1}+K_k\left(z_k-{{\rm B}}_k{\widehat{x}}_{k/k-1}\right)\\ P_k=\left(I-K_k{{\rm B}}_k\right)P_{k/k-1}\end{array}\right.$

in which the Kalman gain is denoted by $K_k$, $I$ represents the identity matrix. As can be seen in (14) that if $R_k$ is not accurate, the estimation would not be reliable. To level up the precision of the navigation solution, $R_k$ should be estimated timely.

3. Maximum Correntropy Extended Kalman Filter and Second Order Mutual Difference Method

3.1. Maximum Correntropy EKF

To further enhance the capabilities of the conventional ILT, the correntropy was originally proposed to use both the statistical and time domain information from the observer(s) [47]. Based on ILT, MCKF was proposed aiming to address the non-Gaussian problems and achieved promising results [48]. MCKF employs the MCC, which is more robust compared to the MMSE. Later, MCEKF was proposed in state estimation for the nonlinear system under outliers; the authors confirmed that the result obtained by MCEKF had much better than the traditional EKF [30].

Considering two stochastic variables $\mathrm{A},\mathrm{B}\in ℝ$, to obtain the correntropy is to calculate the following equation:

(15)$V(\mathrm{A},\mathrm{B})=E[{\kappa }_{ernel}(\mathrm{A},\mathrm{B})]={\displaystyle \int {\kappa }_{ernel}(\alpha ,\beta )d{F}_{\mathrm{A}\mathrm{B}}(\alpha ,\beta )}$

in which $E[\cdot ]$ represents the mathematic expectation, ${\kappa }_{ernel}(\cdot ,\cdot )$ is the kernel function satisfying the Mercer theory, and $F_{\mathrm{A}\mathrm{B}}(\alpha ,\beta )$ is the joint probability distribution of the two variables.

As it is effective and versatile, the Gaussian kernel is usually the widely adopted one, the kernel function is expressed as:

(16)${\kappa }_{ernel}(\alpha ,\beta )=G_{\sigma }(\alpha -\beta )=\exp \left(-\frac{{(\alpha -\beta )}^2}{2{\sigma }^2}\right)$

in which $G_{\sigma }(\cdot )$ represents the Gaussian kernel function, $\sigma$ stands for the kernel size. The details of the selection procedure of our study are given in the experiment section.

In real applications, only limited numbers of data could be used. Moreover, it is not easy to calculate the joint probability density function of random variables. So, the correntropy is usually calculated by the arithmetic operation:

(17)$\widehat{V}(X,Z)=\frac{1}{L}\sum _{i=1}^L{\kappa }_{ernel}\left(e_i\right)$

in which $e_i=x_i-z_i,i=1\dots ,L$, $L$ represents the length of the sampling window.

MCEKF can be seen as a combination of MCC and EKF. When the MCEKF is applied to the INS/GNSS integrated system, the core steps can be listed as follows.

Different from [30], we first calculate the $P_{k/k-1}$ and $R_k$ by the Cholesky decomposition as follows:

(18)$\left\{\begin{array}{l}{P}_{k/k-1}={\Theta }_p{\Theta }_p^T\\ R_k={\Theta }_r{\Theta }_r^T\end{array}\right.$

Then, $e_x$ and $e_z$ are used in the Gaussian kernel function and calculated as:

(19)$\begin{array}{l}{e}_x={\Theta }_p^{-1}\left({\widehat{x}}_k-A_{k,k-1}{\widehat{x}}_{k-1/k-1}\right)\\ e_z={\Theta }_r^{-1}\left(z_k-{{\rm B}}_k{\widehat{x}}_{k/k}\right)\end{array}$

Recalling (16) and (17), the Gaussian kernel correntropy and the cost function are written as [20]:

(20)$J_{MCC}=G_{\sigma }\left(‖e_x‖\right)+G_{\sigma }\left(‖e_z‖\right)$

By solving the following equation, ${\widehat{x}}_{k/k}$ can be obtained:

(21)$\frac{\partial J_{\mathrm{MCC}}}{\partial {\widehat{x}}_{k|k}}=\frac{1}{{\sigma }^2}\left[G_{\sigma }\left(‖e_x‖\right)P_{k/k-1}^{-1}\left({\widehat{x}}_k-A_{k,k-1}{\widehat{x}}_{k-1/k-1}\right)-G_{\sigma }\left(‖e_z‖\right){{\rm B}}_k^T{R}_k^{-1}\left(z_k-{{\rm B}}_k{\widehat{x}}_{k/k}\right)\right]=0$

Note that ${\widehat{x}}_{k/k-1}=A_{k,k-1}{\widehat{x}}_{k-1/k-1}$, if the sampling frequency is high enough, then ${\widehat{x}}_{k/k}={\widehat{x}}_{k/k-1}$ [30], Equation (21) should be rearranged as:

(22)$P_{k/k-1}^{-1}{\widehat{x}}_{k|k}-P_{k/k-1}^{-1}{\widehat{x}}_{k-1/k-1}=L_k{{\rm B}}_k^T{R}_k^{-1}\left(z_k-{{\rm B}}_k{\widehat{x}}_{k/k}\right)$

where:

(23)$L_k=\frac{G_{\sigma }\left(‖e_z‖\right)}{G_{\sigma }\left(‖e_x‖\right)}$

in which ${\widehat{x}}_{k|k}\approx {\widehat{x}}_{k-1/k-1}$ for $L_k$ [30]. So, in MCEKF, $K_k$ in (14) can be replaced as follows:

(24)${\tilde{K}}_k={\left(P_{k/k-1}^{-1}+L_k{{\rm B}}_k{P}_{k/k-1}{{\rm B}}_k^T\right)}^{-1}{L}_k{{\rm B}}_k^T{R}_k^{-1}$

Then, ${\widehat{x}}_{k|k}$ and $P_{k/k}$ can be calculated as:

(25)$\begin{array}{l}{\widehat{x}}_{k|k}={\widehat{x}}_{k|k-1}+{\tilde{K}}_k\left(z_k-{{\rm B}}_k{\widehat{x}}_{k/k}\right)\\ P_{k/k}=\left(I-{\tilde{K}}_k{{\rm B}}_k\right)P_{k/k-1}{\left(I-{\tilde{K}}_k{{\rm B}}_k\right)}^T+{\tilde{K}}_k{R}_k{\tilde{K}}_k^T\end{array}$

It can be found in (25) that ${\widehat{x}}_{k|k}$ are shown at both sides, in MCC based Kalman filters, a fixed-point iteration way can be employed for solving the equation [26]. The fixed-point iterative algorithm in our prosed RMCEKF is calculated by:

(26)${\widehat{x}}_{k|k}^{l+1}={\widehat{x}}_{k|k-1}^l+{\tilde{K}}_k^l\left(z_k-{{\rm B}}_k{\widehat{x}}_{k/k-1}\right)$

where ${\widehat{x}}^l{}_{k|k-1}$ is the result at iteration $l$, ${\tilde{K}}_k^l$ and its factor $L_k^l$ are calculated as:

(27)${\tilde{K}}_k^l={\left(P_{k/k-1}^{-1}+L_k^l{{\rm B}}_k{P}_{k/k-1}{{\rm B}}_k^T\right)}^{-1}{L}_k^l{{\rm B}}_k^T{R}_k^{-1}$

(28)$L_k^l=\frac{G_{\sigma }\left(‖e_z‖\right)}{G_{\sigma }\left(‖e_x^l‖\right)},e_x^l={\Theta }_p^{-1}\left({\widehat{x}}_{k|k}^l-A_{k,k-1}{\widehat{x}}_{k-1}\right)$

To terminate the iteration, a termination condition is defined as:

(29)$\frac{‖{\widehat{x}}_{k|k}^{l+1}-{\widehat{x}}_{k|k}^l‖}{‖{\widehat{x}}_{k|k}^l‖}\le \epsilon$

where $\epsilon \in {ℝ}^{+}$ is a constant small value. When the termination condition is reached, then ${\widehat{x}}_{k|k}={\widehat{x}}_{k|k}^{l+1}$, $K_k={\tilde{K}}_k^l$, and calculate Equation $P_{k/k}$ by (25). Otherwise, set $l=l+1$, continue the next iteration. Other details of MCEKF can be found in [30].

As an important issue, the global convergence analysis of the fixed-point iteration can be found in Ref. [26], so it is omitted here. Similar derivation also can be found in Ref. [17], in which a student’s t kernel is adopted. The key idea is that if the kernel bandwidth is greater than a threshold value, then by Banach fixed-point Theorem, the fixed-point iteration algorithm will surely converge to a unique fixed point. Additionally, Ref. [49] verified that given an appropriate kernel width, the MCEKF could obtain good estimation when the standard EKF even failed to converge under different measurement noises.

It can be found in (23) that the MCEKF takes both the innovation and the residual into the Gaussian kernel. Although the kernel technique can provide higher robustness when the non-Gaussian noises exist. Unfortunately, if the noise encountered by the filter is time-varying and contains many impulsive elements or heavy-tailed properties, the performance of the MCEKF would not be as good as desired [17]. The limitations of MCEKF are given in the next subsection.

3.2. Limitations of MCEKF

In real INS/GNSS applications, If the state estimation error is greater than expected led by some unknown reason, even though the measurement is accurate, the iteration would not take this situation into consideration. Vice versa, if the measured signal is disturbed by impulsive noise or the noise covariance of the measurement violates the preset assumption, the iterative process still could not adjust the strategy to concentrate on the contaminated observation. Based on the above analysis, a method is desired that if it could isolate the noise variance estimate from the state estimate as considerably as possible. Fortunately, the SOMD algorithm mentioned above can meet this requirement.

In the following subsection, the SOMD and its enhanced version which is adopted by this paper will be described for INS/GNSS integrated navigation application.

3.3. Second Order Mutual Difference Algorithm for Estimation of Noise Variance

In the INS/GNSS integrated system, there are two independent measurements for the same signals, i.e., the pseudo-range (s) provided by the GNSS receiver and that calculated by the INS in an indirect way. Let $z_{GN}\left(t\right)$ and $z_{IN}\left(t\right)$ be the measurements obtained from GNSS and INS, respectively, the following redundant measurement model can be made

(30)$\left\{\begin{array}{l}{z}_{GN}\left(t\right)=z_{true}\left(t\right)+{\epsilon }_{GN}\left(t\right)+v_{GN}\left(t\right)\\ z_{IN}\left(t\right)=z_{true}\left(t\right)+{\epsilon }_{IN}\left(t\right)+v_{IN}\left(t\right)\end{array}\right.$

in which $z_{true}\left(t\right)$ is the true signal value, ${\epsilon }_{GN}\left(t\right)$ and ${\epsilon }_{IN}\left(t\right)$ denote the unknown bias item which can be seen as the Gaussian random walk process [39], $t$ is the discrete time epoch, $v_{GN}\left(t\right)$ and $v_{IN}\left(t\right)$ are uncorrelated discrete noise, the statistics of them can be given as

(31)$\left\{\begin{array}{l}E[v_{GN}\left(t\right)v_{GN}^T\left(l\right)]=R_{GN}{\delta }_{t,l}\\ E[v_{IN}\left(t\right)v_{IN}^T\left(l\right)]=R_{IN}{\delta }_{t,l}\end{array}\right.$

in which ${\delta }_{t,l}$ is the Kronecker delta function.

The first order-self-difference (FOSD) of individual observation can be calculated by

(32)$\left\{\begin{array}{l}\mathsf{\Delta }{z}_{GN}(t)=z_{GN}(t)-z_{GN}(t-1)\\ =\left[z_{true}\left(t\right)-z_{true}\left(t-1\right)\right]+\left[{\epsilon }_{GN}\left(t\right)-{\epsilon }_{GN}\left(t-1\right)\right]+\left[v_{GN}\left(t\right)-v_{GN}\left(t-1\right)\right]\\ \mathsf{\Delta }{z}_{IN}(t)=z_{IN}(t)-z_{IN}(t-1)\\ =\left[z_{true}\left(t\right)-z_{true}\left(t-1\right)\right]+\left[{\epsilon }_{IN}\left(t\right)-{\epsilon }_{IN}\left(t-1\right)\right]+\left[v_{IN}\left(t\right)-v_{IN}\left(t-1\right)\right]\end{array}\right.$

The second order-self-difference sequences of $z_{GN}\left(t\right)$ and $z_{IN}\left(t\right)$ can be calculated as

(33)$\begin{array}{l}\mathsf{\Delta }{z}_{GI}(t)=\mathsf{\Delta }{z}_{GN}(t)-\mathsf{\Delta }{z}_{IN}(t)\\ \approx \left[v_{GN}\left(t\right)-v_{GN}\left(t-1\right)\right]-\left[v_{IN}\left(t\right)-v_{IN}\left(t-1\right)\right]\end{array}$

where the items about ${\epsilon }_{GN}$ and ${\epsilon }_{IN}$ are neglected because they are usually much smaller in numerical than the noise [50]. The statistics of noise of $\mathsf{\Delta }{z}_{GN}(t)$ and $\mathsf{\Delta }{z}_{IN}(t)$ could be obtained on the assumption that they are independent and zero-mean:

(34)$\left\{\begin{array}{l}E[v_{GN}\left(t\right)v_{IN}^T\left(t\right)]=0\\ E[v_{GN}\left(t\right)v_{GN}^T\left(t-1\right)]=0\\ E[v_{IN}\left(t\right)v_{IN}^T\left(t-1\right)]=0\end{array}\right.$

and obviously, $E[\mathsf{\Delta }{z}_{GI}(t)]=0$.

The mathematical expectation of $E\left[\mathsf{\Delta }{z}_{GI}(t)\mathsf{\Delta }{z}_{GI}^T(t)\right]$ can be calculated as

(35)$\begin{array}{l}E\left[\mathsf{\Delta }{z}_{GI}(t)\mathsf{\Delta }{z}_{GI}^T(t)\right]=\\ E[v_{GN}\left(t\right)v_{GN}^T\left(t\right)]+E[v_{GN}\left(t-1\right)v_{GN}^T\left(t-1\right)]\\ +E[v_{IN}\left(t\right)v_{IN}^T\left(t\right)]+E[v_{IN}\left(t-1\right)v_{IN}^T\left(t-1\right)]\end{array}$

In real applications, the error accumulation of INS is significantly less than that of GNSS in a short time interval. So, the last two items in (35) can be neglected, $R_{GN}$ can be got by

(36)$R_{GN}=\left(E[v_{GN}\left(t\right)v_{GN}^T\left(t\right)]+E[v_{GN}\left(t-1\right)v_{GN}^T\left(t-1\right)]\right)/2=E\left[\mathsf{\Delta }{z}_{GI}(t)\mathsf{\Delta }{z}_{GI}^T(t)\right]/2$

where $E\left[\mathsf{\Delta }{z}_{GI}(t)\mathsf{\Delta }{z}_{GI}^T(t)\right]$ can be determined by:

(37)$E\left[\mathsf{\Delta }{z}_{GI}(t)\mathsf{\Delta }{z}_{GI}^T(t)\right]=\frac{1}{L}{\displaystyle \sum _{l=0}^{L-1}\mathsf{\Delta }{z}_{GI}(t-l)}\mathsf{\Delta }{z}_{GI}^T(t-l)$

in which $L$ stands for the sample window size.

The SOMD could obtain satisfying performance if the measurement sequences are Gaussian and independent. However, when the GNSS receiver is in an undesirable environment, the two assumptions that (1) the noise distribution conforms to a Gaussian distribution and (2) the measurement difference sequences are uncorrelated would hold anymore.

Although the MCEKF could resist non-Gaussian noise, especially in the outlier’s form, it still suffers from performance degradation when heavy-tailed noise exists [17]. Unfortunately, in INS/GNSS applications, heavy-tailed noise usually exists in challenging environments [23]. Additionally, if the noise arises in a correlated form, the MCEKF would not obtain an accurate estimation as expected.

To address the two problems, we adopt a two-step methodology that has already been proposed by our team. For simplicity, the core ideas are presented as follows, and the details can be found in [50].

First, we use the method established by kernel density estimation (KDE) to detect the potential distribution of the measurement noise. KDE is a powerful tool in that it can evaluate the probability density with a moderate hypothesis about the concerned noise statistics.

Second, if the noise distribution is not Gaussian, then Burg’s algorithm is employed to evaluate $R_{GN}$ by the mutual difference observations instead of SOMD.

4. RMCEKF for the INS/GNSS Integration System

To illustrate the above ideas, the flowchart of RMCEKF is given in Figure 1. The dark green block is a modified MCEKF version proposed by this study, in which the Cholesky decomposition technique is utilized for calculating $P_{k/k-1}$ and $R_k^{}$ and enhancing the numerical stability. The light green block stands for the proposed method; only GNSS signal is available so that the enhanced redundant measurement could be employed.

When the GNSS signal is available, the SOMD sequences are calculated, then the KDE based method is used to make a decision whether the sequences calculated by the SOMD follow the Gaussian distribution or not. Based on the decision, either the standard SOMD algorithm or the Buger’s method is adopted to estimate $R_{GNSS}$ of the GNSS. Finally, $R_{GNSS}$ is employed by the MCEKF in an iterative way to estimate the state ${\hat{\mathrm{x}}}_{k|k}$.

The detail of the RMCEKF is summarized in a pseudocode form as follows.

Initialization: the Gaussian kernel width parameter, $\sigma$, and the fixed-point termination parameter, $\epsilon$, are chosen. ${\widehat{x}}_{0|0}$ and $P_{0|0}^{}$ are set. $Q$ and $R$ are obtained by the INS ‘s and GNSS receiver’s datasheet, respectively (Algorithm 1).

Algorithm 1 The pseudocode of the RMCEKF.

Input: ${\widehat{x}}_{0|0}$, $P_{0|0}^{}$, $Q_{k-1}$, $R_k=diag\left(R_{IN},R_{GN}\right)$, $\sigma$, $\epsilon$.

Prediction:

1. ${\widehat{x}}_{k/k-1}=A_{k,k-1}{\widehat{x}}_{k-1}$.

2. $P_{k/k-1}=A_{k,k-1}{P}_{k-1}{A}_{k,k-1}^T+{{\rm M}}_k{Q}_{k-1}{{\rm M}}_k^T$.

Measurement update (when GNSS signal is available):

1.  Calculate the $\mathsf{\Delta }{z}_{GI}(k)$ by (32)–(33).

2.  Gaussian distribution of $\mathsf{\Delta }{z}_{GI}(k)$ test by KDE method [36]:

if $\mathsf{\Delta }{z}_{GI}(k)\sim Gaussian$

Calculate ${\tilde{R}}_{GN}$ by (36).

Else

${\tilde{R}}_{GN}=\mathrm{Burg}’\mathrm{s}\mathrm{mehtod}(\mathsf{\Delta }{z}_{GI}(k))$.

End

$R_k=diag\left(R_{IN},{\tilde{R}}_{GN}\right)$.

3. ${\Theta }_p=chol\left(P_{k/k-1}\right)$, ${\Theta }_r=chol\left(R_k\right)$.

4. $e_x={\Theta }_p^{-1}\left({\widehat{x}}_k-A_{k|k-1}{\widehat{x}}_{k-1/k-1}\right)$, $e_z={\Theta }_r^{-1}\left(z_k-{{\rm B}}_k{\widehat{x}}_{k/k}\right)$.

5. $L_k^l=\frac{G_{\sigma }\left(‖e_z‖\right)}{G_{\sigma }\left(‖e_x^l‖\right)},e_x^l={\Theta }_p^{-1}\left({\widehat{x}}_{k|k}^l-A_{k,k-1}{\widehat{x}}_{k-1}\right)$.

6. ${\tilde{K}}_k^l={\left(P_{k/k-1}^{-1}+L_k^l{{\rm B}}_k{P}_{k/k-1}{{\rm B}}_k^T\right)}^{-1}{L}_k^l{{\rm B}}_k^T{R}_k^{-1}$.

7. ${\widehat{x}}_{k|k}^{l+1}={\widehat{x}}_{k|k-1}^l+{\tilde{K}}_k^l\left(z_k-{{\rm B}}_k{\widehat{x}}_{k/k-1}\right)$

8. if $\frac{‖{\widehat{x}}_{k|k}^{l+1}-{\widehat{x}}_{k|k}^l‖}{‖{\widehat{x}}_{k|k}^l‖}\le \epsilon$

go to 9.

Else

$l=l+1$, go to 3.

End

9. ${\widehat{x}}_{k|k}={\widehat{x}}_{k|k}^{l+1}$.

$P_{k/k}=\left(I-{\tilde{K}}_k^l{{\rm B}}_k\right)P_{k/k-1}{\left(I-{\tilde{K}}_k^l{{\rm B}}_k\right)}^T+{\tilde{K}}_k^l{R}_k{\tilde{K}}_k^l$.

Outputs: ${\widehat{x}}_{k|k}$, $P_{k/k}$.

5. Experiment and Results

To verify the RMCEKF, we use the data from a road experiment conducted by our team in Tianjin City, China. An INS and a differential GNSS receiver were used in the trial. In Table 1, the primary parameters of the INS are presented. The ground truth is obtained in a post-processing way by the NovAtel Inertial Explorer. The horizontal accuracy of post-processed positioning can reach as high as 1 cm, while the vertical accuracy can reach up to 2 cm.

In Figure 2, the trajectory of this road trial is given. The ground truth is marked by the yellow line, the GNSS outputs is denoted by the red one. It can be seen that there are mainly two intervals, in which the GNSS signals fluctuated due to the poor environment. The challenges include tree foliage in blue blocks and dense buildings in green blocks.

The standard EKF, the MCEKF [30], the ERMAKF [50], the EIAKF [51], and the proposed RMCEKF are used as comparisons.

Scholars adopted several kernel size values in their simulations or experiments to demonstrate the performance of the MCC-based algorithms [26]. Additionally, we tried several kernel sizes among $\sigma =5$, $\sigma =15$, $\sigma =20$, $\sigma =30$, and $\sigma =50$ to find the most appropriate one for the INS/GNSS application. The experiment results showed that when $\sigma =30$, the estimation error is the smallest. For simplicity, we only showed the experiment result when $\sigma =30$.

The parameters used in MCEKF and RMCEKF are the same, i.e.,

(38)$\left\{\begin{array}{l}{\sigma }_{MCEKF}={\sigma }_{RMCEKF}=30\\ {\epsilon }_{MCEKF}={\epsilon }_{RMCEKF}={10}^{-4}\end{array}\right..$

The sliding windows size used in SOMD is set to 50 [36].

The root mean square error (RMSE) is utilized for evaluating the precision of the methods in the following manner.

(39)$RMSE=\sqrt{\frac{1}{L}{\displaystyle \sum _{i=1}^L{\left({\widehat{x}}_i-x_i\right)}^2}}$

in which ${\widehat{x}}_i$ and $x_i$ stand for the navigation solution and the ground truth at the $i\mathrm{th}$ epoch, respectively. $L$ denotes the length of the interval.

Figure 3 shows the acceleration information.

5.1. Estimation Accuracy Comparision in Foliage Environment

It can be found from Figure 4a that the positioning precision of RMCEKF is leveled up by about 17% and 24% compared to ERMAKF and MCEKF. The velocity accuracy of RMCEKF is improved by 39% and 35% compared to standard EKF and MCEKF and 10% worse than ERMAKF. As for the attitude resolution accuracy, it is improved by about 52% and 33% compared to the ERMAKF and MCEKF, respectively.

To furtherly analyze the estimation accuracy of the compared methods, Figure 4b gives the estimation error in longitude, latitude, and height, respectively. For longitude estimation, ERMAKF and RMCEKF performed better than others because of their ability to estimate the GNSS noise variance in line, while the standard EKF and EIAKF reached a maximum error of about 20 m in the interval [90 s, 100 s], while the MCEKF obtained a better accuracy than the EKF and EIAKF because the MCC was used instead of the MMSE. For latitude estimation, RMCEKF performed neck and neck as our previously proposed ERMAKF, the reason why RMCEKF performed not and in longitude and height is one of our future works. For height resolution, RMCEKF performed much better than the compared ones. Figure 4c,d give the details of the velocity and attitude estimation result; ERMAKF and RMCEKF performed better than others; the reason is similar to that in the poisoning estimation.

The RMSE of the navigation errors for different methods is shown in Table 2. The bold type denotes the highest accuracy among the algorithms. The accuracy obtained from the proposed RMCEKF has improved by 23% and 21% in longitude and height estimation compared to the second-best one. MCEKF performed the best in latitude estimation; this is deserved to study and could be a direction for furtherly improving the positioning accuracy. RMCEKF performed in velocity estimation not and in positions, mainly because the measurement noise covariance calculated from the enhanced SOMD is overestimated. For attitude estimation, the proposed RMCEKF performed the best of all. For the sake of simplicity, the comparison of the measurement noise covariance estimation results from GNSS is omitted; furthermore, there is no mechanism for adapting the measurement noise covariance for MCEKF, so readers, please refer to [50] for the detailed comparison result.

5.2. Estimation Accuracy Comparision in Dense Building Area

Compared to Section 5.1, RMCEKF performed even better compared to other ones. Specifically, as can be seen from Figure 5a, for positioning accuracy, it is improved by 35% and 40% compared to ERMAKF and MCEKF. For velocity estimation, the accuracy is refined by about 43% and 53% compared to ERMAKF and MCEKF. It is improved by 35% and 23% in attitude accuracy compared to ERMACK and MCEKF.

In Figure 5b–d, the estimation errors comparison in the dense building environment at intervals [260 s, 320 s] is given. The reason why ERMAKF and RMCEKF performed better than the others is similar to that in Section 5.1. However, it is worth to emphasis that RMCEKF performed much better in this environment than in the foliage environment because the measurement noises in this environment exhibited a heavy-tailed distribution characteristic, which is consistent with the insightful analysis in [52]. In that study, the authors concluded that in the dense urban area (Shanghai, China), the multipath-contaminated line-of-sight (LOS) signal reception follows the Gamma distribution, and non-line-of-sight (NLOS) signal reception follows the exponential distribution, respectively. Clearly, both are heavy-tailed styles.

The RMSE of the navigation errors for different methods in the dense building area is shown in Table 3. The bold type denotes the highest accuracy among the algorithms. It can be found that all the solutions provided by the proposed RMCEKF maintained the lowest error except for the north velocity item, which was only about 1% worse than the best one (MCEKF). Specifically, in positioning error compared to ERMAKF, our method achieved an improvement of about 45% (1.4543 m vs. 2.6509 m), 1% (1.4129 m vs. 1.4191 m), 37% (2.1926 m vs. 3.4597 m) in longitude, latitude, and height estimations, respectively, which resulted in about 43% improvement in 3-D positioning error.

5.3. Discussion

In the first challenging environment of the road experiment, our proposed method could achieve the best estimation accuracy among the compared ones in 3-D positioning, but in the latitude item, there is still space for improvement. Additionally, the estimation of velocity is not as good as ERMAKF and MCEKF. Thus, one of our future works should be focused on this issue. In the second dense building area, RMCEKF performed much better than in the first environment. This is consistent with the argument that the MCC performs better than the MMSE if the noises follow the Non-Gaussian distribution, especially in a heavy-tailed form.

This paper focuses on the tightly coupled mechanism. However, if the GNSS receiver could just provide the position and the velocity information instead of the pseudo-range and pseudo-range rate data, then only a loosely coupled integrated system could be employed. As the measurement function is different from that of the tightly coupled system, the calculations of SOMD should be revised accordingly.

Currently, the development of autonomous driving technology and machine learning is getting faster. The cameras and LiDARs are essential sensors in addition to GNSS in these fields [53]. These sensors also may generate non-Gaussian measurement information [54]. Therefore, the method proposed in this paper can be combined with existing work by other scholars to improve their performance. Additionally, this is still a very meaningful direction for our future work.

6. Conclusions

In this study, to address the challenges of the INS/GNSS integrate system caused by complex environments, a novel Redundant measurement-based Maximum Correntropy Extended Kalman Filter (RMCEKF) is presented. This novel method used the independent measurement from INS and GNSS to evaluate the noise covariance of the GNSS signal, then MCEKF was employed instead of the standard EKF. A real road experiment verified the methodology, and the positioning accuracy in the plant area and dense building area were improved by 17% and 35% compared to ERMAKF, respectively.

Author Contributions

Conceptualization, H.Z. and D.W.; methodology, D.W. and H.H.; software, D.W; validation, H.Z., D.W., H.H. and B.G.; formal analysis, H.Z. and D.W.; investigation, H.Z. and D.W.; resources, H.Z. and D.W.; data curation, H.Z.; writing—original draft preparation, D.W.; writing—review and editing, H.Z., D.W., H.H. and B.G.; visualization, D.W.; supervision, H.Z.; project administration, H.Z.; funding acquisition, H.Z. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

We thank all the reviewers and the editors for the insightful comments and suggestions that hugely help to improve the quality of the paper.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Diagram of the RMCEKF.

Figure 2. The ground trajectory of the road experiment.

Figure 3. The raw 200 Hz acceleration data.

Figure 4. Estimation errors comparison in foliage environment at interval [40 s, 120 s]. (a) Comparison of 3-D poisoning, velocity, and attitude for different methods. (b) Position errors are estimated by different algorithms. (c) Velocity errors are estimated by different algorithms. (d) Attitude errors are estimated by different algorithms.

Figure 5. Estimation errors comparison in building environment at interval [260 s, 320 s]. (a) Comparison of 3-D poisoning, velocity, and attitude for different methods. (b) Position errors for different methods. (c) Velocity errors. (d) Attitude errors for different methods.

Table 1

The key parameters of IMU-ISA-100C and GNSS.

	Gyroscope Parameters	Accelerator Parameters
Measurement range	±495 deg/s	±10 g
Bias stability(In-run)	≤0.05 deg/hr	≤100 μg
Non-linearity (Scale factor)	≤100 ppm	≤100 ppm
Random walk	0.012 deg/√hr	≤100 μg/√Hz
Raw Data Rate	200 Hz	200 Hz
	Horizontal Position Accuracy of GNSS (RMS)
Single point L1/L2	1.2 m
SBAS	60 cm
DGPS	40 cm

Table 2

Navigation Errors for Compared Algorithms in Foliage Environment.

Algorithm	Longitude(m)	Latitude(m)	Height(m)	East Velocity(m/s)	North Velocity(m/s)	Up Velocity(m/s)	Pitch(°)	Roll(°)	Yaw(°)
Standard EKF	5.4944	1.2624	5.6871	0.6415	0.2246	0.2221	0.2785	0.2165	0.0485
EIAKF	6.4422	1.8996	4.5233	0.6779	0.3255	0.1585	0.2784	0.2645	0.0454
ERMAKF	5.0846	1.9428	4.1101	0.2884	0.2544	0.1006	0.2783	0.3518	0.0312
MCEKF	5.1309	1.1844	5.3577	0.6000	0.2009	0.2183	0.2616	0.2032	0.0483
RMCEKF	3.9347	2.6159	3.2446	0.3112	0.2588	0.1792	0.1710	0.1294	0.0497

Table 3

Navigation Errors for Compared Algorithms in Dense Building Area.

Algorithm	Longitude(m)	Latitude(m)	Height(m)	East Velocity(m/s)	North Velocity(m/s)	Up Velocity(m/s)	Pitch(°)	Roll(°)	Yaw(°)
Standard EKF	3.1620	1.7765	4.1326	0.3017	0.1075	0.1870	0.1653	0.0830	0.1347
EIAKF	2.9980	4.7459	3.2645	0.1961	0.1597	0.1011	0.1921	0.0888	0.1837
ERMAKF	2.6209	1.4191	3.4597	0.2152	0.1456	0.1037	0.1620	0.1902	0.0868
MCEKF	2.8755	1.7279	3.6706	0.2758	0.1020	0.1733	0.1573	0.0750	0.1361
RMCEKF	1.4543	1.4129	2.1926	0.0853	0.1031	0.0783	0.0884	0.0538	0.1346

References

1. Chen, H.; Wu, W.; Zhang, S.; Wu, C.; Zhong, R. A GNSS/LiDAR/IMU Pose Estimation System Based on Collaborative Fusion of Factor Map and Filtering. Remote Sens.; 2023; 15, 790. [DOI: https://dx.doi.org/10.3390/rs15030790]

2. Xu, Y.; Wang, K.; Jiang, C.; Li, Z.; Yang, C.; Liu, D.; Zhang, H. Motion-Constrained GNSS/INS Integrated Navigation Method Based on BP Neural Network. Remote Sens.; 2023; 15, 154. [DOI: https://dx.doi.org/10.3390/rs15010154]

3. Li, Z.; Zhang, H.; Zhou, Q.; Che, H. An Adaptive Low-Cost INS/GNSS Tightly-Coupled Integration Architecture Based on Redundant Measurement Noise Covariance Estimation. Sensors; 2017; 17, 2032. [DOI: https://dx.doi.org/10.3390/s17092032] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28872629]

4. Fang, W.; Jiang, J.; Lu, S.; Gong, Y.; Tao, Y.; Tang, Y.; Yan, P.; Luo, H.; Liu, J. A LSTM Algorithm Estimating Pseudo Measurements for Aiding INS during GNSS Signal Outages. Remote Sens.; 2020; 12, 256. [DOI: https://dx.doi.org/10.3390/rs12020256]

5. Li, X.; Chen, W.; Chan, C.; Li, B.; Song, X. Multi-sensor fusion methodology for enhanced land vehicle positioning. Inf. Fusion; 2019; 46, pp. 51-62. [DOI: https://dx.doi.org/10.1016/j.inffus.2018.04.006]

6. Xia, X.; Xiong, L.; Huang, Y.; Lu, Y.; Gao, L.; Xu, N.; Yu, Z. Estimation on IMU yaw misalignment by fusing information of automotive onboard sensors. Mech. Syst. Signal Process.; 2022; 162, 107993. [DOI: https://dx.doi.org/10.1016/j.ymssp.2021.107993]

7. Gao, L.; Xiong, L.; Xia, X.; Lu, Y.; Yu, Z.; Khajepour, A. Improved Vehicle Localization Using On-Board Sensors and Vehicle Lateral Velocity. IEEE Sens. J.; 2022; 22, pp. 6818-6831. [DOI: https://dx.doi.org/10.1109/JSEN.2022.3150073]

8. Groves, P.D. Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems; Artech House: London, UK, 2008; pp. 363-387.

9. Xia, X.; Hashemi, E.; Xiong, L.; Khajepour, A. Autonomous Vehicle Kinematics and Dynamics Synthesis for Sideslip AngleEstimation Based on Consensus Kalman Filter. IEEE Trans. Control. Syst. Technol.; 2022; 31, pp. 179-192. [DOI: https://dx.doi.org/10.1109/TCST.2022.3174511]

10. Xiong, L.; Xia, X.; Lu, Y.; Liu, W.; Gao, L.; Song, S.; Yu, Z. IMU-based automated vehicle body sideslip angle and attitude estimation aided by GNSS using parallel adaptive Kalman filters. IEEE Trans. Veh. Technol.; 2020; 69, pp. 10668-10680. [DOI: https://dx.doi.org/10.1109/TVT.2020.2983738]

11. Liu, W.; Xia, X.; Xiong, L.; Lu, Y.; Gao, L.; Yu, Z. Automated Vehicle Sideslip Angle Estimation Considering Signal Measurement Characteristic. IEEE Sens. J.; 2021; 21, pp. 21675-21687. [DOI: https://dx.doi.org/10.1109/JSEN.2021.3059050]

12. Duník, J.; Biswas, S.K.; Dempster, A.G.; Pany, T.; Closas, P. State estimation methods in navigation: Overview and application. IEEE Aerosp. Electron. Syst. Mag.; 2020; 35, pp. 16-31. [DOI: https://dx.doi.org/10.1109/MAES.2020.3002001]

13. Khodaparast, J. A Review of Dynamic Phasor Estimation by Non-Linear Kalman Filters. IEEE Access; 2022; 10, pp. 11090-11109. [DOI: https://dx.doi.org/10.1109/ACCESS.2022.3146732]

14. Bar-Shalom, Y.; Li, X.R.; Kirubarajan, T. Estimation with Applications to Tracking and Navigation: Theory Algorithms and Software; John Wiley & Sons: Hoboken, NJ, USA, 2004; pp. 381-394.

15. Julier, S.J.; Uhlmann, J.K. Unscented filtering and nonlinear estimation. Proc. IEEE.; 2004; 92, pp. 401-422. [DOI: https://dx.doi.org/10.1109/JPROC.2003.823141]

16. Arasaratnam, I.; Haykin, S. Cubature kalman filters. IEEE Trans. Autom. Control.; 2009; 54, pp. 1254-1269. [DOI: https://dx.doi.org/10.1109/TAC.2009.2019800]

17. Huang, H.; Zhang, H. Student’s t-Kernel-Based Maximum Correntropy Kalman Filter. Sensors; 2022; 22, 1683. [DOI: https://dx.doi.org/10.3390/s22041683]

18. Shongwe, T.; Vinck, A.J.H.; Ferreira, H.C. A Study on Impulse Noise and Its Models. SAIEE Afr. Res. J.; 2015; 106, pp. 119-131. [DOI: https://dx.doi.org/10.23919/SAIEE.2015.8531938]

19. Middleton, D. Non-Gaussian noise models in signal processing for telecommunications: New methods and results for class A and class B noise models. IEEE Trans. Inf. Theory; 1999; 45, pp. 1129-1149. [DOI: https://dx.doi.org/10.1109/18.761256]

20. Ndo, G.; Labeau, F.; Kassouf, M. A markov-middleton model for bursty impulsive noise: Modeling and receiver design. IEEE Trans. Power Deliv.; 2013; 28, pp. 2317-2325. [DOI: https://dx.doi.org/10.1109/TPWRD.2013.2273942]

21. Peng, Y.; Chen, J.; Xu, X.; Pu, F. SAR Images Statistical Modeling and Classification Based on the Mixture of Alpha-Stable Distributions. Remote Sens.; 2013; 5, pp. 2145-2163. [DOI: https://dx.doi.org/10.3390/rs5052145]

22. Zhang, Z.; Wang, J.; Li, Z.; Zhao, Y.; Wang, R.; Habib, A. Optimization Method of Airborne LiDAR Individual Tree Segmentation Based on Gaussian Mixture Model. Remote Sens.; 2022; 14, 6167. [DOI: https://dx.doi.org/10.3390/rs14236167]

23. Guangcai, W.; Xu, X.; Zhang, T. MM estimation-based robust cubature Kalman filter for INS/GPS integrated navigation system. IEEE Trans. Instrum. Meas.; 2020; 70, pp. 1-11. [DOI: https://dx.doi.org/10.1109/TIM.2020.3021224]

24. Li, Y.; Mi, J.; Xu, Y.; Li, B.; Jiang, D.; Liu, W. A Robust Adaptive Filtering Algorithm for GNSS Single-Frequency RTK of Smartphone. Remote Sens.; 2022; 14, 6388. [DOI: https://dx.doi.org/10.3390/rs14246388]

25. Chang, L.; Hu, B.; Chang, G.; Li, A. Huber-based novel robust unscented Kalman filter. IET Sci. Meas. Technol.; 2012; 6, pp. 502-509. [DOI: https://dx.doi.org/10.1049/iet-smt.2011.0169]

26. Chen, B.; Liu, X.; Zhao, H.; Principe, J.C. Maximum correntropy Kalman filter. Automatica; 2017; 76, pp. 70-77. [DOI: https://dx.doi.org/10.1016/j.automatica.2016.10.004]

27. Principe, J.C. Information Theoretic Learning: Renyi’s Entropy and Kernel Perspectives; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2010; pp. 123-126.

28. Qu, H.; Wang, M.; Zhao, J.; Zhao, S.; Li, T.; Yue, P. Generalized Asymmetric Correntropy for Robust Adaptive Filtering: A Theoretical and Simulation Study. Remote Sens.; 2022; 14, 3677. [DOI: https://dx.doi.org/10.3390/rs14153677]

29. Liao, T.; Hirota, K.; Wu, X.; Shao, S.; Dai, Y. A Dynamic Self-Tuning Maximum Correntropy Kalman Filter forWireless Sensors Networks Positioning Systems. Remote Sens.; 2022; 14, 4345. [DOI: https://dx.doi.org/10.3390/rs14174345]

30. Wang, S.Y.; Yang, J.Z. State Estimation under Outliers by the Maximum Correntropy Extended Kalman Filter. Proceedings of the 2021 60th Annual Conference of the Society of Instrument and Control Engineers of Japan (SICE); Tokyo, Japan, 8–10 September 2021; pp. 1426-1431.

31. Wang, D.; Zhang, H.; Ge, B. Adaptive Unscented Kalman Filter for Target Tacking with Time-Varying Noise Covariance Based on Multi-Sensor Information Fusion. Sensors; 2021; 21, 5808. [DOI: https://dx.doi.org/10.3390/s21175808]

32. Zhang, S. An Adaptive Unscented Kalman Filter for Dead Reckoning Systems. Proceedings of the 2009 International Conference on Information Engineering and Computer Science; Wuhan, China, 19–20 December 2009.

33. Lu, D.H.; Wang, J.Q.; Xiong, K.; Hou, B.W.; He, Z.M. Strong tracking Kalman filter for non-Gaussian observation. Kongzhi Lilun Yu Yingyong/Control Theory Appl.; 2019; 36, pp. 1997-2004.

34. Xie, S.B.; Zhang, Y.; Wen, K.R.; Zhang, S.; Liu, Z.M.; Qi, N.M. Motion estimation for non-cooperative target based on strong tracking cubature Kalman filter. Jilin Daxue Xuebao J. Jilin Univ.; 2021; 51, pp. 1482-1489.

35. Wang, M.; Dong, X.; Qin, C.; Liu, J. Adaptive H-infinity Kalman filter based random drift modeling and compensation method for ring laser gyroscope. Measurement; 2021; 167, 108170. [DOI: https://dx.doi.org/10.1016/j.measurement.2020.108170]

36. Gao, S.; Zhong, Y.; Li, W. Random weighting method for multisensory data fusion. IEEE Sens. J.; 2011; 11, pp. 1955-1961. [DOI: https://dx.doi.org/10.1109/JSEN.2011.2107896]

37. Gao, S.; Hu, G.; Zhong, Y. Windowing and random weighting-based adaptive unscented Kalman filter. Int. J. Adapt. Control Signal Process.; 2015; 29, pp. 201-223. [DOI: https://dx.doi.org/10.1002/acs.2467]

38. Gao, Z.; Mu, D.; Wei, W.; Zhong, Y.; Gu, C. Adaptive unscented Kalman filter based on maximum posterior estimation and random weighting. Aerosp. Sci. Technol.; 2017; 71, pp. 12-24. [DOI: https://dx.doi.org/10.1016/j.ast.2017.08.020]

39. Huang, H.; Zhang, H.; Jiang, L. An Optimal Fusion Method of Multiple Inertial Measurement Units Based on Measurement Noise Variance Estimation. IEEE Sens.; 2023; 23, pp. 2693-2706. [DOI: https://dx.doi.org/10.1109/JSEN.2022.3229475]

40. Velazquez, J.R. Analysis and Development of Algorithms for Data Fusion in Sensor Arrays. Ph.D. Thesis; Universit Montpellier: Montpellier, France, 2020.

41. Srinivas, P.; Kumar, A. Overview of architecture for GPS-INS integration. Proceedings of the 2017 Recent Developments in Control, Automation & Power Engineering; Noida, India, 26–27 October 2017.

42. Xu, X.; Nie, Z.; Wang, Z.; Wang, B.; Du, Q. Performance Assessment of BDS-3 PPP-B2b/INS Loosely Coupled Integration. Remote Sens.; 2022; 14, 2957. [DOI: https://dx.doi.org/10.3390/rs14132957]

43. Zhai, R.; Yuan, Y. A Method of Vision Aided GNSS Positioning Using Semantic Information in Complex Urban Environment. Remote Sens.; 2022; 14, 869. [DOI: https://dx.doi.org/10.3390/rs14040869]

44. Chiang, K.W.; Duong, T.T.; Liao, J.K. The performance analysis of a real-time integrated INS/GPS vehicle navigation system with abnormal GPS measurement elimination. Sensors; 2013; 13, pp. 10599-10622. [DOI: https://dx.doi.org/10.3390/s130810599]

45. Zhao, Y. Performance evaluation of Cubature Kalman filter in a GPS/IMU tightly-coupled navigation system. Signal Proc.; 2016; 119, pp. 67-79. [DOI: https://dx.doi.org/10.1016/j.sigpro.2015.07.014]

46. Kong, X.; Nebot, E.M.; Durrant-Whyte, H. Development of a nonlinear psi-angle model for large misalignment errors and its application in INS alignment and calibration. Proceedings of the 1999 IEEE International Conference on Robotics and Automation; Detroit, MI, USA, 10–15 May 1999; pp. 1430-1435.

47. Santamaría, I.; Pokharel, P.P.; Principe, J.C. Generalized correlation function: Definition, properties, and application to blind equalization. IEEE Trans. Signal Process.; 2006; 54, pp. 2187-2197. [DOI: https://dx.doi.org/10.1109/TSP.2006.872524]

48. Xiong, W.; Schindelhauer, C.; So, H.C.; Wang, Z. Maximum Correntropy Criterion for Robust TOA-Based Localization in NLOS Environments. Circuits Syst. Signal Process.; 2021; 40, pp. 6325-6339.

49. Mohiuddin, S.M.; Qi, J. Maximum Correntropy Extended Kalman Filtering for Power System Dynamic State Estimation. Proceedings of the IEEE Power & Energy Society General Meeting (PESGM); Atlanta, GA, USA, 4–8 August 2019; pp. 1-5. [DOI: https://dx.doi.org/10.1109/PESGM40551.2019.8973525]

50. Ge, B.S.; Zhang, H.; Fu, W.X.; Yang, J.B. Enhanced Redundant Measurement-Based Kalman Filter for Measurement Noise Covariance Estimation in INS/GNSS Integration. Remote Sens.; 2020; 12, 3500. [DOI: https://dx.doi.org/10.3390/rs12213500]

51. Ghaleb, F.; Zainal, A.; Rassam, M.; Abraham, A. Improved vehicle positioning algorithm using enhanced innovation-based adaptive Kalman filter. Pervasive Mob. Comput.; 2017; 40, pp. 139-155. [DOI: https://dx.doi.org/10.1016/j.pmcj.2017.06.008]

52. Chen, X.; Yu, J.M.; Yu, W.; Truong, T.K. GPS L1CA/BDS B1I Multipath Channel Measurements and Modeling for Dynamic Land Vehicle in Shanghai Dense Urban Area. IEEE Trans. Veh. Technol.; 2020; 69, pp. 14247-14263. [DOI: https://dx.doi.org/10.1109/TVT.2020.3038646]

53. Xia, X.; Meng, Z.; Han, X.; Li, H.; Tsukiji, T.; Xu, R.; Zhang, Z.; Ma, J. Automated Driving Systems Data Acquisition and Processing Platform. arXiv; 2022; arXiv: 2211.13425

54. Liu, W.; Quijano, K.; Crawford, M.M. YOLOv5-Tassel: Detecting Tassels in RGB UAV Imagery With Improved YOLOv5 Based on Transfer Learning. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2022; 15, pp. 8085-8094. [DOI: https://dx.doi.org/10.1109/JSTARS.2022.3206399]