Enhanced Redundant Measurement-Based Kalman Filter for Measurement Noise Covariance Estimation in INS/GNSS Integration

Abstract

Adaptive Kalman filters (AKF) have been widely applied to the inertial navigation system (INS)/global navigation satellite system (GNSS) integrated navigation system. However, the traditional AKF methods suffer from the problems of filtering instability or covariance underestimation, especially when the GNSS measurement disturbances occur. In this paper, an enhanced redundant measurement-based AKF is developed to improve the filtering performance. The scheme is based on the mutual difference sequence derived from the redundant measurement of INS. By using the mutual difference sequence, the measurement noise covariance can be estimated without being affected by the inaccuracy estimates, hence avoiding the risk of filtering divergence. In addition, the kernel density estimation is used to estimate the GNSS measurement noise’s probability density to detect whether the Gaussian properties of the measurement noise are maintained. When the noise statistics are far from Gaussian distribution, the difference sequence will be modeled as an autoregressive process using the Burg’s method. The real variance of the difference sequence can then be updated relying on the autoregressive model in order to avoid the covariance underestimation. A field experiment was carried out to evaluate the performance of the proposed method. The test results demonstrate that the proposed method can effectively mitigate the GNSS measurement disturbances and improve the accuracy of the navigation solution.

1. Introduction

The integrated navigation system composed of the inertial navigation system (INS) and the global navigation satellite system (GNSS) is being widely used to provide the complete navigation states (i.e., position, velocity and orientation) due to their complementary characteristics [1]. In the integrated system, INS as a self-contained system can provide high-bandwidth, low short-term noise outputs, but its output accuracy degrades over time; GNSS solutions have long-term stability in an open sky environment, but its update rate is low and the solution accuracy suffers from degradation when the ideal navigation environment changes. Therefore, the fusion of the INS and GNSS information using a Kalman filter (KF) can combine their advantages so as to prevent the drifting of INS solution and smooth the GNSS position output [2]. Furthermore, considering the different architectures for INS/GNSS integration, nonlinear filters such as extended Kalman filter (EKF), unscented Kalman filter (UKF) and particle filter (PF) are commonly used to deal with the nonlinearity of the integrated system. The method applied by the EKF algorithm for approximating these nonlinearities is the Taylor series. In the EKF method, the nonlinear function is linearized through calculating the term of the first order of the Taylor series about the state point. However, both the EKF and KF method can only obtain the optimal estimates by the criterion of minimum mean square error when the statistical properties of the system parameters are known and the measurement noises are in Gaussian distribution [3,4]. The UKF method utilizes the unscented transformation to parameterize the first two moments of the random state which undergo a nonlinear transformation. This kind of transformation does not require that the noise distribution belongs to Gaussian. Moreover, it can obtain more accurate mean and covariance than that achieved by nonlinear function approximation because of the notion that it should be easier to approximate a distribution than to linearize a nonlinear function [5]. Although the UKF algorithm seems to have considerable merit compared with the EKF method, the unknown statistical properties of noise is still a stumbling block to its optimal estimation [6]. Different from the UKF method, which is based on deterministic sampling, the PF method applies a large number of random particles to approximately represent the states’ posterior density distribution [7]. The state estimates can be then obtained from the posterior density distribution by using the Bayesian filtering theory [8]. Nonetheless, due to the large number of samples, the PF method usually requires heavy computation burden. Meanwhile, the degeneracy problem and sample impoverishment always haunt PF [9]. Therefore, the computation cost and its own flaws limit the applications of PF in real-time navigation.

In practical integration navigation applications, especially in some challenging environments such as urban canyons, foliage covered areas and open pits, the GNSS signals are subject to interference. Therefore, the measurement noise statistic characteristics are uncertain and time varying. Besides, the uncertainty involved in the system state model also has a serious influence on the process noise covariance. Considering the noise covariance variation, the conventional KF and nonlinear filters are usually unable to be directly used to estimate the navigation states. The inadequate prior assumption of the noise covariance will significantly influence the localization performance of the INS/GNSS integrated system and potentially will lead to filter divergence [10,11,12]. To improve the performance of the integrated system, the concept of adaptive Kalman filter (AKF) was employed to deal with the inadequate noise covariance [13,14]. Generally, these methods fall into four categories: Bayesian methods, maximum likelihood estimations, covariance matching approaches, and correlation techniques [15,16]. The basic idea of Bayesian methods is to obtain a recursion equation for the posterior probability density function (PDF) of the state from the measurement and then estimate the parameters by integrating the PDF. In the maximum likelihood estimation, the likelihood function is maximized to obtain the noise statistic estimates. Although both of them are theoretically sound, they require heavy computations and have been falling out of application in favor of INS/GNSS integration. Covariance matching approaches balance the trade-off between the computation and the performance. Yet the core idea behind these methods is intuitive. Thus, the matching method’s convergence is not always accurate and gives biased covariance estimates [17]. Correlation techniques are fruitful in noise covariance estimation, where a set of linear equations relating the components of the noise statistics to the autocorrelation of the innovation or residual sequences are established. Through solving these linear equations, the measurement or process noise covariance can be obtained. One of the most popular correlation methods is innovation-based adaptive estimation (IAE), in which an estimate of the noise covariance is provided by a linear transformation of the innovation covariance [18,19].

Although the IAE method can improve the accuracy of the INS/GNSS integration navigation solutions by using the estimated noise covariance, the stability of the IAE method is compromised. As well-known, the IAE method is based on the innovation, which is easily influenced by both the contaminated measurements as well as inaccurate estimates. Therefore, the estimated covariance using IAE has the risk to decline in accuracy and even results in filtering divergence, especially in the case where the statistical properties of both measurement noise and process noise are needed to be estimated [20]. Aiming for solving this problem, Zhang et al. [21] proposed a redundant measurement-based adaptive estimation (RMAE) for the INS/GNSS integrated navigation system, in which the INS solutions are used as redundant measurements to estimate the measurement noise covariances. As the estimates of the measurement noise covariance are only based on the measurements, the RMAE method is immune to the state errors. In view of this advantage, Ge et al. [22] combined the RMAE method with the algorithms for process noise covariance estimation so as to avoid the coupling problem of estimating the measurement and process noise covariance simultaneously. In the INS/GNSS navigation, however, the GNSS measurement noise may no longer be in a Gaussian distribution and the measurements are correlated in time or even involve many random outliers, especially in multipath afflicted environments. Thus, the measurement noise covariance obtained using RMAE may be underestimated due to ignoring the autocorrelation when the mutual difference sequence is non-Gaussian. To improve the position accuracy and reliability, Dhital et al. [23] assumed that the GNSS measurement noise follows a heavy-tailed distribution in the challenging environments such as deep canyons and adjusted the measurement noise covariance using a scalar based on the acceleration consistency between the GNSS and INS. Ghaleb et al. [24] modeled the correlated measurement residual sequence as a first-order autoregressive process AR(1), and then used the Yule-Walker method to estimate the measurement noise covariance. The measurement residual sequence collected in estimation, however, neglects the steady-state error of GNSS and INS. Therefore, the measurement residual sequence is not capable of representing the measurement noises properly. Moreover, the randomness testing based on the improper measurement residuals may lead to the misjudgment of the optimality of KF.

In this work, we focus on the measurement noise covariance estimation. Based on the RMAE method, an enhanced redundant measurement-based adaptive Kalman filter is developed, which extends RMAE to the case with correlated measurements. Furthermore, the abnormal observations are also considered. The random weighting method is introduced into the measurement covariance estimation to control outliers’ precision levels, which overcomes the limitation of the RMAE method; that is, all samples from historical epochs make an equivalent contribution to the measurement covariance estimation. The improved RMAE method can accurately achieve the varying measurement noise covariance without relying on states estimation. Thus, the coupling problem of covariance estimation can be solved by combining the enhanced RMAE method with other schemes for the estimation of process noise covariance.

The remainder of this paper is organized as follows: Section 2 describes the INS/GNSS tightly coupled scheme briefly; Section 3 introduces the implementation of the enhanced redundant measurement-based adaptive Kalman filter; the field test results are discussed in Section 4, and Section 5 summarizes this paper with a conclusion.

2. INS/GNSS Tightly Coupled Integration

According how GNSS measurements are used in the integration algorithm, the way of corrections applied to the inertial navigation solution, and whether the GNSS user equipment is aided, the INS/GNSS integration architecture can be mainly divided into three groups: loosely coupled, tightly coupled and ultra-tightly coupled integration systems [25].

The loosely coupled system has a simple structure, where the outputs of the GNSS navigation filter are directly input as measurements to the integration filter. However, the errors of the filtering outputs are time correlated, so the usage of a cascaded architecture may disrupt the state estimation of KF. Different from the loosely coupled integration, the GNSS navigation filter and the integration filter are combined into one in the tightly coupled integration to get rid of the statistical problem of cascading. Furthermore, this architecture has a suppression effect on GNSS multipath and still works even when the number of visible satellites drops to one [26]. The main feature of an ultra-tightly coupled system is that the INS is used to aid GNSS receiver acquisition and tracking. In the ultra-tight coupling, the corrected INS solution is used to generate the commands of the GNSS receiver’s numerically controlled oscillator. Thus, compared with tight coupling, the ultra-tightly coupled integration can operate with a lower tracking bandwidth. Nevertheless, this architecture is usually very complex and computationally intensive, and a new interface between GNSS receiver and INS need to be designed. Hence, the tightly coupled INS/GNSS integration is widely preferred for its good robustness as well as easy implementation.

2.1. Dynamic Model

In the tightly coupled integration, the system state model can be established by combining an INS and a GNSS state error propagation. The states of INS consist of the attitude error φⁿ, velocity error δVⁿ and position error (δL, δλ, δh), together with the gyroscope drift ε^b and accelerometer error ∇^b, while the states additionally considered in GNSS are the GNSS user equipment’s clock offset δt_u and its drift δt_ru.

The attitude error propagation in the local-level frame, which can also be stated as the navigation (n-)frame, is formulated as

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

(1)

where ${\mathit{\omega }}_{ie}^n$ is the earth rate expressed in n-frame, ${\mathit{\omega }}_{en}^n$ denotes the rate of n-frame with respect to the earth-centered, earth-fixed (e-)frame, $C_b^n$ is the direction cosine matrix from the body (b-)frame to n-frame.

The transition function of the velocity error vector in the n-frame can be described as

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

(2)

where ${\mathit{f}}^n$ and $\delta {\mathit{g}}^n$ represent the specific force and gravitational error, respectively.

The position error propagation is shown as follows

$$\{\begin{array}{l}\delta \dot{L}=\frac{\delta V_y^n}{R_M+h}-\frac{V_y^n}{{(R_M+h)}^2}\delta h\\ \delta \dot{\mathit{\lambda }}=\frac{\delta V_x^n}{(R_N+h)\cos L}-\frac{V_x^n}{{(R_N+h)}^2\cos L}\delta h+\frac{V_x^n·\tan L}{(R_N+h)\cos L}\delta L\\ \delta \dot{h}=V_z^n\end{array}$$

(3)

where δL, δλ, δh are the position error components in the latitude, longitude and height. $V_x^n$, $V_y^n$, $V_z^n$ are the velocity components in the east, north, and up direction. R_M and R_N denote the meridian and transverse radius of curvature, respectively.

The error model of gyroscope is assumed to be the combination between zero-offset bias and first-order Gauss-Markov (GM) process, thus

$$\{\begin{array}{l}\mathit{\epsilon }={\mathit{\epsilon }}^b+{\mathit{\epsilon }}^r\\ {\dot{\mathit{\epsilon }}}^b=0\\ {\dot{\mathit{\epsilon }}}^r=-\frac{1}{T_g}{\mathit{\epsilon }}^r+{\mathit{\omega }}_g\end{array}$$

(4)

where ${\mathit{\epsilon }}^b$ and ${\mathit{\epsilon }}^r$ denote gyroscope bias and first-order GM process, respectively. $T_g$ is the correlation time, ${\mathit{\omega }}_g$ is the driving noise.

The stochastic model for accelerometer error is modeled as

$${\dot{\mathit{\nabla }}}^b=-\frac{1}{T_a}{\mathit{\nabla }}^b+{\mathit{w}}_a$$

(5)

where $T_a$ is the correlation time and ${\mathit{w}}_a$ is Gaussian white noise.

For the state errors of GNSS, the propagation functions are

$$\{\begin{array}{l}\delta {\dot{t}}_u=\delta t_{ru}+w_u\\ \delta {\dot{t}}_{ru}=-\frac{1}{T_{ru}}\delta t_{ru}+w_{ru}\end{array}$$

(6)

where T_ru is the correlation times, w_u and w_ru are Gaussian white noises.

The error propagation model of the measurements can be represented as follows

$$\{\begin{array}{l}\delta {\rho }_i={\rho }_{Gi}-{\rho }_{Ii}\\ \delta {\dot{\rho }}_i={\dot{\rho }}_{Gi}-{\dot{\rho }}_{Ii}\end{array}$$

(7)

where ${\rho }_{Gi}$ and ${\dot{\rho }}_{Gi}$ are the pseudo-range and pseudo-range rate, respectively. ${\rho }_{Ii}$ and ${\dot{\rho }}_{Ii}$ represent the INS-driven pseudo-range and pseudo-range rate.

The pseudo-range ${\rho }_{Gi}$ and pseudo-range rate ${\dot{\rho }}_{Gi}$ can be calculated using the broadcasted ephemeris from the satellites, whose models are as follows

$$\{\begin{array}{l}{\rho }_{Gi}=r_i+\left(\delta t_u-\delta T_i\right)c+{\mathsf{\Delta }}_{ioni}+{\mathsf{\Delta }}_{tropi}+v_{\rho i}\\ {\dot{\rho }}_{Gi}={\dot{r}}_i+\left(\delta {\dot{t}}_u-\delta {\dot{T}}_i\right)c+v_{\dot{\rho }i}\\ r_i=\sqrt{{\left(x-x_i\right)}^2+{\left(y-y_i\right)}^2+{\left(z-z_i\right)}^2}\\ {\dot{r}}_i=\frac{x-x_i}{r_i}\frac{\partial x}{\partial t}+\frac{y-y_i}{r_i}\frac{\partial y}{\partial t}+\frac{z-z_i}{r_i}\frac{\partial z}{\partial t}\end{array}$$

(8)

where $r_i$ denotes the true geometric distance between the receiver and satellite, and ${\dot{r}}_i$ is the true range rate. $\delta T_i$ and $\delta {\dot{T}}_i$ are the satellite clock offset and drift which can be directly obtained from the broadcasted ephemeris. c is the speed of light. ${\mathsf{\Delta }}_{ion}$ and ${\mathsf{\Delta }}_{trop}$ are propagation delays of ionosphere and troposphere, respectively. These propagation errors can be mitigated by the Klobuchar model and troposphere delay models, such as Hopfield and Egnos models [27,28]. $v_{\rho i}$ and $v_{\dot{\rho }i}$ are the measurement noise. (x,y,z) denotes the true position of the receiver in e-frame and (x_i,y_i,z_i) is the position of the satellite, which can also be obtained from the ephemeris.

The INS-driven pseudo-range ${\rho }_{Ii}$ and pseudo-range rate are represented as

$$\{\begin{array}{l}{\rho }_{Ii}=\sqrt{{\left(x_I-x_i\right)}^2+{\left(y_I-y_i\right)}^2+{\left(z_I-z_i\right)}^2}\\ {\dot{\rho }}_{Ii}=\frac{x_I-x_i}{{\rho }_{Ii}}\frac{\partial x_I}{\partial t}+\frac{y_I-y_i}{{\rho }_{Ii}}\frac{\partial y_I}{\partial t}+\frac{z_I-z_i}{{\rho }_{Ii}}\frac{\partial z_I}{\partial t}\end{array}$$

(9)

where (x_I,y_I,z_I) denotes the inertial measurement unit’s position in the e-frame, which can be converted from the INS position solution.

2.2. Navigation Fusion with EKF

The non-linearity of the measurement model limits the application of KF. Thus, an EKF method is used to fuse the outputs of the INS and GNSS system. In the EKF method, the non-linear models are linearized applying Taylor series expansion, which makes the measurement models suitable for the linear filtering conditions. The implementation of EKF can be described as follows

Time Updating

$$\{\begin{array}{l}{\hat{X}}_{k/k-1}=F_{k,k-1}{\hat{X}}_{k-1}\\ P_{k/k-1}=F_{k,k-1}{P}_{k-1}{F}_{k,k-1}^T+{\mathbf{\Gamma }}_{k-1}{Q}_{k-1}{\mathbf{\Gamma }}_{k-1}^T\end{array}$$

(10)

where ${\hat{X}}_{k-1}$ is the state estimate at time k−1 and X = [φⁿ, δV_n, δL, δλ, δh, ε^b, ε^r, ∇^b, δt_u, δt_ru], P_k−1 is the estimation covariance. ${\hat{X}}_{k/k-1}$ and P_k/k−1 denote the state prediction and prediction covariance. F_k,k−1 is the state transition matrix. ${\Gamma }_{k-1}$ and Q_k−1 are the noise-driven matrix and process noise covariance matrix, respectively.

Measurement Updating

$$\{\begin{array}{l}{\widehat{X}}_k={\widehat{X}}_{k/k-1}+K_k\left(Z_k-H_k{\widehat{X}}_{k/k-1}\right)\\ K_k=P_{k/k-1}{H}_k^T{(H_k{P}_{k/k-1}{H}_k^T+R_k)}^{-1}\\ P_k=\left(I-K_k{H}_k\right)P_{k/k-1}\end{array}$$

(11)

where Z_k denotes the difference between the GNSS pseudo-range (rate) and the INS-driven pseudo-range (rate). H_k is the measurement matrix and ${\mathbf{H}}_{\mathrm{k}}=\partial \mathbf{h}/\partial {\widehat{\mathbf{X}}}_{k/k-1}$, in which h represents the nonlinear measurement function. K_k denotes the Kalman gain and R_k is the measurement noise covariance matrix.

For the INS/GNSS tightly coupled system, the EKF method is the problem. Although it can solve the nonlinear filtering problem and obtain better navigation results, the time-varying measurement noise covariance and abnormal observations caused by the disturbances of challenging environments may result in inaccurate solutions or even divergence problems. Therefore, the measurement noise covariance used in the filter should be adapted to that of the real-time environments.

3. Enhanced Redundant Measurement-Based AKF

The RMAE method is immune to the inaccurate estimates and performs well in the case where the measurement stochastic noises are Gaussian. Nonetheless, the assumption of the white noise for GNSS measurement is not reasonable, especially in the harsh environment where the multipath errors are correlated over time. In addition, the outliers involved in the measurements will also limit the accuracy of the covariance estimation. Thus, an enhanced redundant measurement-based AKF scheme is developed in this section. The enhancement consists of three parts: a measurement-based robust adaptive noise covariance estimation for Gaussian noise, a switching strategy between the RMAE method and the Burg’s algorithm according to the statistics analysis of the second-order mutual difference sequences, and the measurement noise covariance estimation for non-Gaussian noise.

3.1. Redundant Measurement-Based Robust Adaptive Measurement Noise Covariance EstimationSubsection

The observation equations expressed in (8) and (9) show that there are two separate observations available for the same measurement values. GNSS provides the pseudo-range and pseudo-range rate directly, while INS provides them in an implicit manner. Based on the RMAE method, the measurements Z_GNSS and Z_INS from GNSS and INS can respectively be modeled as [29]

$$\{\begin{array}{c}{Z}_{GNSS}\left(k\right)=Z\left(k\right)+f_{GNSS}\left(k\right)+V_{GNSS}\left(k\right)\\ Z_{INS}\left(k\right)=Z\left(k\right)+f_{INS}\left(k\right)+V_{INS}\left(k\right)\end{array}$$

(12)

where Z(k) denotes the true measurement at time epoch k, f_GNSS(k) and f_INS(k) are the steady-state errors of the measurements, V_GNSS(k) and V_INS(k) are mutually independent Gaussian noises.

Subtracting every two adjacent measurements, the difference sequences of the GNSS and INS measurement yield as

$$\{\begin{array}{cc}\hfill \Delta Z_{GNSS}\left(k\right)& =Z_{GNSS}\left(k\right)-Z_{GNSS}\left(k-1\right)\hfill \\ & =Z\left(k\right)-Z\left(k-1\right)+f_{GNSS}\left(k\right)-f_{GNSS}\left(k-1\right)+V_{GNSS}\left(k\right)-V_{GNSS}\left(k-1\right)\\ \hfill \Delta Z_{INS}\left(k\right)& =Z_{INS}\left(k\right)-Z_{INS}\left(k-1\right)\hfill \\ & =Z\left(k\right)-Z\left(k-1\right)+f_{INS}\left(k\right)-f_{INS}\left(k-1\right)+V_{INS}\left(k\right)-V_{INS}\left(k-1\right)\end{array}$$

(13)

Then, the second-order mutual difference sequence of the measurements can be calculated as

$$\nabla Z\left(k\right)=\mathsf{\Delta }{Z}_{GNSS}\left(k\right)-\mathsf{\Delta }{Z}_{INS}\left(k\right)\approx V_{GNSS}\left(k\right)-V_{GNSS}\left(k-1\right)-\left[V_{INS}\left(k\right)-V_{INS}\left(k-1\right)\right]$$

(14)

where the steady item f of the measurement errors is much smaller than the noise, so it can be neglected after subtraction in (14).

The two separate measurement noises are uncorrelated and zero-mean,

$$\{\begin{array}{l}E\left[V_{GNSS}\left(k\right)V_{INS}^T\left(k\right)\right]=0\\ E\left[V_{GNSS}\left(k\right)V_{GNSS}^T\left(k-1\right)\right]=0\\ E\left[V_{INS}\left(k\right)V_{INS}^T\left(k-1\right)\right]=0\\ E\left[\nabla Z\left(k\right)\right]=0\end{array}$$

(15)

where E[ ] is the expectation operator.

Thus, the autocorrelation of the second-order mutual difference sequences can be written as

$$\begin{array}{ll}E\left[\nabla Z\left(k\right)\nabla Z^T\left(k\right)\right]& =E\left[V_{GNSS}\left(k\right)V_{GNSS}^T\left(k\right)\right]+E\left[V_{GNSS}\left(k-1\right)V_{GNSS}^T\left(k-1\right)\right]+\\ & E\left[V_{INS}\left(k\right)V_{INS}^T\left(k\right)\right]+E\left[V_{INS}\left(k-1\right)V_{INS}^T\left(k-1\right)\right]\end{array}$$

(16)

On the other side, the INS error accumulation is much smaller than the GNSS error in a short term. Therefore, the measurement noise covariance of GNSS can be approximated as

$$R_{GNSS}=E\left[V_{GNSS}{V}_{GNSS}^T\right]\approx E\left[\nabla Z\left(k\right)\nabla Z^T\left(k\right)\right]/2$$

(17)

The calculation of the autocorrelation E[∇Z(k)∇Z^T(k)] can utilize a limited number of samples of the second-order mutual difference sequence,

$$E\left[\nabla Z\left(k\right)\nabla Z^T\left(k\right)\right]=\frac{1}{M}{\displaystyle \sum _{m=0}^{M-1}\nabla Z\left(k-m\right)\nabla Z^T\left(k-m\right)}$$

(18)

where M denotes the size of the window, which is chosen by experience as it implies a trade-off between the smoothness and sharpness of the covariance estimation.

The RMAE method is decoupled from state estimation errors and suitable for the case without any random outliers involved in the GNSS measurement. However, if the abnormal observations occur, the samples from different sequences are likely to be contaminated and cannot accurately describe the statistical properties of the measurement noise, thus leading to the bias for measurement noise covariance estimation. To improve the robustness of the RMAE method, the random weights concept [30] is applied to (18), the autocorrelation can be then obtained using random weighting estimation as

$$E\left[\nabla Z\left(k\right)\nabla Z^T\left(k\right)\right]=\frac{1}{M}{\displaystyle \sum _{m=0}^{M-1}{w}_m\nabla Z\left(k-m\right)\nabla Z^T\left(k-m\right)}$$

(19)

where w_m denotes the random weighting factor from the Dirichlet distribution, which meets the condition ${\displaystyle \sum _{m=0}^{M-1}{w}_m}=1$. The details on their determination can be found in [31].

3.2. Detection of GNSS Measurement Noise Statistics

Since the GNSS measurements are likely to be disturbed by the surroundings around the GNSS receiver, the statistical properties of the measurement noise may be far from the assumption and change to be non-Gaussian as well [32]. Before the stage of covariance estimation for the correlated noise, a detection to determine whether the measurement noise is Gaussian or not is introduced first in this subsection.

Considering the disturbance of the challenging environment would contaminate the observation’s statistical properties significantly, it is not easy to find out which family of the distribution the observation is from when the disturbance occurs. Unlike the parametric estimation assuming that the observed data belong to a given parametric family of distributions, the kernel density estimation (KDE) as a nonparametric estimation just assumes that the probability density of distribution exists [33,34]. So, the KDE method can be used to estimate the probability density under a less rigid assumption about the measurement noise statistics.

For the sake of testing accuracy, instead of innovation or residual which is also influenced by the inaccurate estimates, the second-order mutual difference sequence that relies only on the measurement is used to test the statistical properties of the measurement noise. Suppose that the independent samples ∇Z₁, ∇Z₂, …, ∇Z_N are extracted from the second-order mutual difference sequence, the PDF of the data sequence can be estimated by KDE method as

$${\widehat{f}}_n\left(z\right)=\frac{1}{Nh}{\displaystyle \sum _{i=1}^{N}K\left(\frac{z-\nabla Z_i}{h}\right)}$$

(20)

where N is the size of sampling window, h is called bandwidth, K(⋅) denotes the kernel function.

Generally, the global estimation accuracy of ${\widehat{f}}_n\left(z\right)$ is measured with the mean integrated square error (MISE), which is defined by

$$MISE\left[{\widehat{f}}_n\left(z\right)\right]=E{{\displaystyle \int \left\{{\widehat{f}}_n\left(z\right)-f_n\left(z\right)\right\}}}^{2}dz$$

(21)

Theories and experiments have verified that there is very little influence on the biases of MISE when the different kernel functions are selected in the context of sufficient sample size [35]. Therefore, the kernel function is commonly chosen as the Gaussian kernel function in practice, which is defined as

$$K\left(u\right)=\frac{1}{\sqrt{2\pi }}{e}^{-\frac{u^2}{2}}$$

(22)

Along with the Gaussian kernel function, the optimal bandwidth that minimizes the MISE can be obtained by the rule of thumb [33]

$$h_{opt}=1.06\sigma N^{-1/5}$$

(23)

where σ is the variance of the observed sequence, which can be obtained by the adaptive estimation described in Section 3.1.

According to the definition of PDF, the probability of the data falls around z can be written as

$$f\left(z\right)=\frac{dF\left(z\right)}{dz}=\underset{h\to 0}{\lim }\frac{F\left(z+h/2\right)-F\left(z-h/2\right)}{h}$$

(24)

where F(z) denotes the cumulative distribution function.

Equation (24) indicates that f(z) can represent the probability of each data point that falls around z in the second-order mutual difference sequence. Furthermore, the data’s credibility can be measured through the probability as well. As illustrated in Figure 1, the dashed lines represent the individual bumps placed at the samples from the difference sequence interval ${\left\{z\right\}}_{t-1}^{t-n}$ and the solid line is the density estimate which is constructed though adding each bump up. When the current datum sample z_t shown by the red circle appears at the edge of the kernel density estimate, the probability of this point is very small. So, the red datum point is not likely to be from the same PDF of the previous data ${\left\{z\right\}}_{t-1}^{t-n}$ and it implies that the GNSS measurement may be contaminated by the surroundings. If z_t represented by the blue circle falls in the interior of the kernel density estimate, it is determined that this point belongs to the same PDF as the previous data ${\left\{z\right\}}_{t-1}^{t-n}$. The credibility of the current datum can directly be calculated by (20) in real time. If the measurement noise is normal, the credibility of each datum will satisfy the condition that ${\widehat{f}}_n\left(z_t\right)>K\left(3\sigma \right)$; Otherwise the measurement disturbance is inferred to occur and the Gaussian properties of the measurement noise may change when the condition does not hold continuously.

3.3. Measurement Noise Covariance Estimation for Non-Gaussian Noise

The RMAE method is sensitive to the correlated errors induced by the environmental noises that affect the GNSS receiver. Ignoring this autocorrelation of the measurement, especially in the multipath environment, can lead to underestimation of the measurement noise covariances. Affected by these underestimates, the filter will trust the measurement’s correctness more and cause serious implications on the filtering performance. Therefore, when the changes of the Gaussian properties of the measurement noise are detected, the second-order mutual difference sequence should be considered as a corrected sequence rather than a Gaussian noise process. Ghaleb et al. [24] used an AR(1) process to describe the corrected noise. However, the autoregressive model order also has to be estimated according to the practical noise process. Furthermore, the Yule-Walker method is unstable and may lead to poor parameter estimation in some special cases [36].

The Burg’s method is a widely employed algorithm for autoregressive modelling as it can provide reliable parameter estimates and a stable model on the contrary to the Yule-Walker and covariance method. Unlike the Yule-Walker method which can be expressed in a single set of equations, the Burg’s method operates on the difference sequence in iterations. Once the autoregressive model order p of the difference sequence is estimated using Akaike’s criterion [37], the Burg’s recursion can be implemented as follows [38,39].

Step1: Initialize the AR coefficients $a_0^{(0)}=1$ as well as the forward prediction error vector ${\mathbf{f}}_0^{(0)}$ and the backward prediction error vector ${\mathbf{b}}_0^{(0)}$ using the difference sequence, namely,

$$f_0^{(0)}=b_0^{(0)}=\left[\nabla Z_1,\nabla Z_2,\cdots ,\nabla Z_N\right]$$

(25)

where N denotes the number of difference sequence samples.

Step2: Calculate the reflection coefficient

$$k_l=-2\frac{{\displaystyle \sum _{n=l}^N{f}_n^{(l-1)}}{b}_{n-1}^{(l-1)}}{{\displaystyle \sum _{n=l}^N\left[{\left(f_n^{(l-1)}\right)}^2+{\left(b_{n-1}^{(l-1)}\right)}^2\right]}}$$

(26)

where l is the stage.

Step3: Update the prediction errors and the AR coefficients a_m via the Levinson-Durbin recursion.

$$\{\begin{array}{l}{f}_n^{(l)}=f_n^{(l-1)}+k_l{b}_{n-1}^{(l-1)}\\ b_n^{(l)}=b_{n-1}^{(l-1)}+k_l{f}_{n-1}^{(l-1)}\begin{array}{cc}& n=l,l+1,\cdots ,N\end{array}\\ a_m^{(l)}=a_m^{(l-1)}+k_l{a}_{l-m}^{(l-1)}\begin{array}{cc}& m=1,2,\cdots ,l\end{array}\\ a_l^{(l)}=k_l\end{array}$$

(27)

where m and n are the symbols for counting.

Step4: Repeat steps from 2 to 3 until l = p.

After the iterations, the AR model that describes the correlation of the measurement noise can be constructed as

$$\nabla Z_N={\displaystyle \sum _{i=1}^p{a}_i\nabla Z_{N-i}}+v$$

(28)

where v is a zero-mean white noise process.

According to Ghaleb in [24], the actual variance of the difference sequence can be then obtained as follows

$${\sigma }_{\nabla Z}^2=\frac{{\displaystyle \sum _{m=0}^{M-1}\nabla Z\left(k-m\right)\nabla Z^T\left(k-m\right)}}{M\left(1-{\displaystyle \sum _{i=1}^p{a}_i{\beta }_i}\right)}$$

(29)

where β_i is the autocorrelation coefficient of the difference sequence.

3.4. Implementation of the Enhanced AKF

To illustrate the above descriptions, the flowchart of the proposed enhanced redundant measurement-based AKF method is shown in Figure 2. The main steps are as follows: After the states and covariance initialization, the time-propagation is operated to predict the states using the knowledge of the a priori processing noise covariance of the INS/GNSS integrated system. Once the measurements from INS and GNSS are available, the second-order mutual difference observations will be formed. Then, the KDE-based detection approach will be applied to identify the GNSS measurement noise statistics. If the statistical properties of the measurement noise are stable and belong to the family of Gaussian distribution, the robust RMAE method will be selected. Otherwise, the measurement noise covariance will be expanded and estimated by the Burg’s method using the stream of the mutual difference observations. Finally, the updated measurement noise covariance is adopted to optimally update the state estimates during the measurement-update phase of EKF.

4. Test Results

4.1. Implementation of the Enhanced AKF

In order to validate the proposed algorithm, a field experiment was carried out in Tianjin City, China. The NovAtel IMU-ISA-100C device, whose bias offsets of the gyroscope and the accelerometer are 0.5 °/h and 1250 μg and the angle random walk of the gyroscope and the velocity random walk of the accelerometer are respectively 0.012°/$\sqrt{\mathsf{h}}$ and 100 μg/$\sqrt{\mathrm{Hz}}$, was used to collect the raw inertial data at the rate of 100 Hz. The NovAtel GNSS receiver applied in the experiment can output the pseudo-range and pseudo-range rate data at the rate of 1 Hz. In our test, about 30 min data was collected and its post-processing positioning results processing by the NovAtel Inertial Explorer using the differential GNSS measurements in the Tightly Coupled mode were used as the ground truth. Figure 3 shows the trajectory for the test. It is clear that the GNSS outputs show a fluctuation and even deviates from the road due to the environmental disturbances when the vehicle went through the foliage or buildings.

4.2. Performance Evaluation in Foliage Environment

The navigation performance of the enhanced redundant measurement-based AKF (ERMAKF) was evaluated in comparison with the standard EKF method, an enhanced innovation-based adaptive Kalman filter (EIAKF) [24], and the simple redundant measurement-based AKF (SRMAKF) which ignored the autocorrelation of the measurement disturbance. The vehicle went through the foliage during 40–120 s periods and the navigation errors with respect to the ground truth are shown in Figure 4, Figure 5 and Figure 6. Figure 7 shows the noise standard deviations of the pseudo-range $\left(\sqrt{r_1}\sim \sqrt{r_4}\right)$ and pseudo-range rate $\left(\sqrt{r_5}\sim \sqrt{r_8}\right)$ for each filter.

As shown in Figure 4, the navigation errors of the INS/GNSS tightly integrated system using the standard EKF algorithm increased significantly in the period of 50–120 s, when the vehicle was driving under the dense vegetation. The maximum error in longitude even reaches about 18 m. This can be explained by the multipath effect, which will contaminate the GNSS measurements and lead to large navigation errors. Figure 7 shows the fixed measurement noise standard deviations applied in the standard EKF. The priori values of noise covariance will lower the filtering accuracy due to the overweighed abnormal measurements. Although the EIAKF method can turn the measurement noise covariance adaptively, its noise covariance estimates may not match that of the actual measurement. This is because the estimated measurement noise covariance in the EIAKF method is based on the innovation sequence, which is also influenced by the inaccurate state estimates. In addition, the measurement residual sequence used in the EIAKF method is directly obtained by differencing the GNSS measurements and the “dead reckoning” results, which ignores the steady-state errors. It can be seen from the estimated measurement noise standard deviations of the EIAKF method in Figure 7 that the standard deviations do not update in a timely manner due to the improper difference sequence applied to the randomness testing which may lead to the misjudgment for the filtering optimality. Thus, the filtering errors of the EIAKF algorithm, as shown in Figure 4, Figure 5 and Figure 6, fluctuate quite a lot. Different from the EIAKF algorithm, the SRMAKF method estimates the measurement noise covariance relying on the redundant mutual difference sequence, which is not only immune to the state estimation error, but also considers the steady-state error of the measurement systems. Therefore, the navigation errors of the SRMAKF method are decreased compared with the EIAKF algorithm, although it causes an increase in the height error. The reason for the decrease in the height accuracy is that the altitude solutions from the GNSS are usually unsteady and even correlated in time due to the imperfect tropospheric delay model and obstructed environments. As can be seen from Figure 7, such correlation will result in the underestimation of the noise standard deviations. Compared with the EIAKF algorithm, ERMAKF takes the measurement autocorrelation into full consideration. When the measurement autocorrelation is detected using the KDE method, the Burg’s method is used to estimate the parameter of the AR model and the variance of the difference sequence is then obtained according to the AR coefficients. Figure 8, Figure 9 and Figure 10 give insight about the navigation accuracy of each algorithm in terms of the root mean square error (RMSE), which is defined as

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}{\displaystyle \sum _{k=1}^N{\left({\widehat{X}}_k-X_k\right)}^2}}$$

(30)

where ${\widehat{X}}_k$ and $X_k$ are the filtering and reference states at time k, respectively.

As shown in Figure 8, the longitude solution accuracy is improved by 26.3% and 46.9% when the SRMAKF scheme is used, in comparison with the Standard EKF and EIAKF solution, respectively. It is clear that the ERMAKF scheme achieves the best performance in the longitude component, whose accuracy is further improved by 36.0% compared to the SRMAKF scheme due to the reduction in the negative influence of the measurement correlation on the noise covariance estimation. As explained before, the positioning RMSE of the EIAKF method has an increase because of the unreliable difference sequence applied in the measurement noise covariance estimation. Although the SRMAKF method can improve the longitude accuracy, the height solution deteriorates when the GNSS altitude is strongly correlated. In comparison, the height solution applying the ERMAKF scheme outperforms that of the Standard EKF and EIAKF method by 56.0% and 50.0%, respectively. The improvement in the latitude accuracy is not that significant, but it still maintains an acceptable accuracy.

The RMSE of the velocity accuracy is presented in Figure 9. It can be seen that compared with the Standard EKF solutions, the SRMAKF method improves the velocity solutions by 62.8%, 26.7% and 43.5% in the east, north and up direction, respectively. Similarly, SRMAKF obtains improvements of 68.0%, 52.2% and 18.7% compared to the EIAKF method. Moreover, the velocity accuracy can be further improved when the ERMAKF scheme is used. It should be noted that although the altitude solved by SRMAKF includes more errors, the velocity accuracy it obtains in the up direction is not the same and shows better accuracy. This is because the noise covariance of the pseudo-range and pseudo-range rate are estimated separately and these estimates will determine the measurement’s weights in filtering. Even though the accuracy of the pseudo-range and pseudo-range rates is interlinked with each other at the output from the receiver’s ranging processor, the following navigation processing will break this direct correspondence of accuracy between the two observations due to the changes in the estimated covariance while filtering. It is not surprising that the ERMAKF scheme performs better than the SRMAKF method. This is due to the fact that the SRMAKF method is not able to deal with the measurement correlations and it may underestimate the noise covariance.

As can be seen from Figure 10, although the high-precision attitude is achievable through the standard EKF method, the ERMAKF method can reduce the RMSE further, especially for the pitch and roll component. The disturbances of the abnormal GNSS measurements are rejected through de-weighting the observations, and thus the ERMAKF solutions are improved by 47.4% for the pitch component and 21.4% for the roll component compared to the standard EKF method. The discussions of the EIAKF and SRMAKF algorithms are similar as aforementioned. In general, the attitude performance using the ERMAKF scheme is satisfactory.

4.3. Performance Evaluation in Obstructed Areas

Dense buildings constitute a harsh environment, where the GNSS signals are susceptible to the disturbances. As shown in Figure 3, the positioning quality of GNSS was degraded when the vehicle passed that buildings during 260–320 s periods. In this dense building environment, the ERMAKF method can also be applied to mitigate the GNSS measurement errors. To quantify the filtering errors concisely, the RMSE in 3-D space is listed in

Table 1. 3-Dimensional Navigation Errors for Different Methods.

Algorithm	Position Error (m)	Velocity Error (m/s)	Attitude Error (°)
Standard EKF	5.52	0.38	0.22
EIAKF	6.24	0.29	0.25
SRMAKF	6.12	0.38	0.24
ERMAKF	4.60	0.17	0.18

.

Table 1 illustrates the filtering performance for different methods. It shows that although the adaptive strategies are used in the INS/GNSS integrated system, the navigation accuracy may not be improved in some special scenarios without considering the model’s influence and noise’s correlation. The EIAKF method can estimate the time-varying measurement noise covariance applying the innovation, but the mismatched error models also have a negative influence on the estimation, hence leading to the decrease in filtering accuracy. The SRMAKF method is immune to the state errors. However, the correlation between the measurements will cause the underestimation of the measurement covariance. Thus, the filtering accuracy of the SRMAKF method can still be improved. The ERMAKF method is an enhanced version of SRMAKF, in which the measurement’s autocorrelation is considered and the Burg’s method is employed to model this autocorrelation. Through classifying the measurement noise and updating their statistics, the positioning error decreases from 6.12 m to 4.60 m while applying the ERMAKF method. The velocity errors are respectively reduced by 41.4% and 55.3% compared to EIAKF and SRMAKF when the ERMAKF method is used. Meanwhile, the attitude errors obtained by EIAKF also appeared to decline, falling by 25.0% compared to the SRMAKF method. The results demonstrate that the proposed enhanced method can provide more consistent navigation solutions than the other state-of-the-art algorithms.

5. Conclusions

This paper presents a scheme of the mitigation of GNSS measurement disturbances, and the methods of RMAE and Burg’s are used to estimate the time-varying measurement noise covariance based on the statistic property of the mutual difference sequence. The scheme mainly includes the following steps: first, to obtain the mutual difference observations by using two separate measurements from INS and GNSS; second, to detect via the KDE method, whether the measurement noise statistics belong to a Gaussian distribution, and then determine the autoregressive model applying the Burg’s recursion when the Gaussian properties of the difference sequence noise are not maintained; finally, to obtain the navigation results through fusing the sensors data with the timely updating measurement noise covariance.

In the field experiment, the fusing results demonstrate that the proposed enhanced redundant measurement-based AKF method outperforms the compared methods in general. The positioning accuracy is improved by 66.0%, 20.7% and 50.0%, respectively, for longitude, latitude, and height components compared to the conventional innovation-based adaptive Kalman filter in foliage environment. However, there is a slight decrease in the yaw accuracy due to the windowing estimation errors. In addition, the ERMAKF method also obtains an improvement for the 3-dimensional navigation accuracy in the harsh environment such as dense buildings.

The current research mainly focuses on the feasibility of the proposed enhanced method and only the time-varying measurement noise covariance is inspected. In the future, the processing noise covariance estimation will be considered for improving the navigation precision, especially for the low-cost integrated system.