Modeling of Inertial Rate Sensor Errors Using Autoregressive and Moving Average (ARMA) Models

Abstract

In this chapter, a low-cost micro electro mechanical systems (MEMS) gyroscope drift is modeled by time series model, namely, autoregressive-moving-average (ARMA). The optimality of ARMA (2, 1) model is identified by using minimum values of the Akaike information criteria (AIC). In addition, the ARMA model based Sage-Husa adaptive fading Kalman filter algorithm (SHAFKF) is proposed for minimizing the drift and random noise of MEMS gyroscope signal. The suggested algorithm is explained in two stages: (i) an adaptive transitive factor (a1) is introduced into a predicted state error covariance for adaption. (ii) The measurement noise covariance matrix is updated by another transitive factor (a2). The proposed algorithm is applied to MEMS gyroscope signals for reducing the drift and random noise in a static condition at room temperature. The Allan variance (AV) analysis is used to identify and quantify the random noise sources of MEMS gyro signal. The performance of the suggested algorithm is analyzed using AV for static signal. The experimental results demonstrate that the proposed algorithm performs better than CKF and a single transitive factor based adaptive SHFKF algorithm for reducing the drift and random noise in the static condition.

1. Introduction

During the last two decades, low cost, small size, accurate and reliable navigation system development is a hot research in the modern navigation technology. The word navigation is a process of monitoring and controlling any moving object from one place to other. Inertial navigation system (INS) is a dead reckoning positioning method based on measurements and mathematical processing of the vehicle absolute acceleration and angular speed in order to estimate its attitude, speed and position related to difference[1, 2, 3, 4, 5]. INS technology is categorized into (i) gimbal INS and (ii) strap-down INS. In the early 1940s, a gimbal INS system was developed based on the mechanical inertial sensor (i.e., accelerometers and gyroscopes) for providing the navigation information [5]. Its accuracy was limited by mechanical inertial sensor errors. The main drawbacks of the gimbal INS system are its designed complexity and it requires synchronous servo motors, slip rings, control electronics, etc., for acquiring the navigation information. Because of these factors, the gimbals INS systems are used in low grade navigation applications [6]. In the early 1950, strap-down INS (SINS) was developed based on solid state inertial sensor [5]. The SINS is a self-contained navigation system that has been developed for providing the accurate navigation information (i.e., position, velocity and rotation information. It has three gyroscopes ad three accelerometers. In general, the operation principle of SINS follows the physical laws of motion equations. It is an emerging technology as compared to gimbal INS systems and it has significant features such as easy to design, lower cost of ownership, moderate manufacturing cost and also high reliability. SINS consist of an inertial measurement that includes 3-axis accelerometers and 3-axis gyroscopes, and a processing computer. IMU is a key device to the INS and has been widely used for measuring the rotation rate and acceleration of an object. In practice, SINS accuracy degrades due to internal and external errors of the inertial sensors. These errors are mainly caused due to fluctuation in temperature, pressure and internal electronics components of the sensor. Due to these factors, stochastic errors and drift errors are generated at the IMU output [7, 8].

With the recent development of modern navigation technology, inertial sensor based SINS technology have been characterized into three categories, (i) low accuracy (tactical applications), (ii) medium accuracy (navigation applications) and (iii) high accuracy (strategic navigation applications) sensor technology. The performance improvements of inertial sensors are decided by the inertial sensor errors [9]. Currently, the strap-down INS use (i) low-cost MEMS and (ii) precision fiber optic gyroscope. MEMS sensor has more attractive to manufacturers of navigation systems because of their small size, low cost, light weight, low power consumption and ruggedness [10]. However, MEMS sensors give poor performance in the highly dynamic environment. Hence, the reliability of MEMS-based INS navigation accuracy is limited. Because of these features MEMS have only been used for low-end navigation applications (i.e., commercial domain) [11].

In the recent years, MEMS devices have been developed and tested successfully for low-end accuracy applications [12, 13]. MEMS sensor operates for a long time under poor condition and it generates the noise due to internal circuits and electronics interferences of the MEMS sensor [14, 15, 16]. As a result, noise and drift are generated at the MEMS output. In general, drift error is affected by ambient temperatures and magnetic field effect [17, 18, 19]. Many studies have been reported for temperature error model of MEMS sensor to capture the temperature variation affects [20]. According to the IEEE standard specification, MEMS errors can be characterized into two categories, such as (i) deterministic errors and (ii) stochastic errors. Deterministic errors are due to scale factor errors, bias and misalignment errors [18, 19]. Several calibration methods have been developed for eliminating the bias errors, scale factor errors in the lab environments. Stochastic errors are due to quantization effect, temperature effect (random bias), random drift, and additive noise of MEM sensor. In the case of stochastic errors analysis, calibration techniques cannot be suitable because of randomness [21, 22, 23, 24]. This chapter concentrates on random errors modeling and random noise elimination techniques. The developments of random noise suppressing methods are helpful for improving the MEMS accuracy as well as SINS accuracy. In general stochastic error includes quantitation noise (QN), bias instability (BS), angle random walk (ARW), rate random walk (RRW) and rate ramp (RR) drift. With the extension of research, random noise and bias drift are the non-negligible errors in the MEMS sensor output. In this chapter, different signal processing techniques are developed to minimize the bias drift and random noise [25].

In time domain, Allan Variance (AV) is a popular technique has been widely used to identify and quantify different random noises present in the MEMS sensor [16, 26, 27]. In literature, several noise compensation techniques such as discrete wavelet transform (DWT), empirical mode decomposition (EMD) method and Forward linear prediction (FLP) methods have been developed and applied to MEMS sensors for filtering the high-frequency noise [28, 29]. These methods are not suitable when the sensor includes the correlated noise. Kalman filter (KF) is a most popular state estimation technique that has been used for minimizing the correlated noise of the MEMS sensor [30, 31, 32, 33, 34]. The priori knowledge of an initial values of the process and measurement noise covariance matrix are known exactly, when the KF become an optimal. However, in practice, these noise parameters may vary with time so that the performance of the KF can be degraded and then the filter become diverge.

To solve the divergence problems, Adaptive Kalman filter technique (AKF) is a better solution. The adaptation can be based on either (i) innovation based adaptive estimation AKF (IAE-AKF) or (ii) residual based estimation AKF (RAE-AKF) and also multiple model based AKF [34, 35]. Among the other methods, adaptive KF is developed using IAE. In general, an innovation sequence is defined as the difference between true and estimated values. In the IAE-AKF method, the measurement and process noise matrices are estimated based on innovation sequence and followed by sliding average window method. In real time, the selection of window size is a critical issue. Sage-Husa Adaptive KF is another version of adaptive KF that has been developed to improve the AKF performance by introducing a time varying estimator. In the SHAKF, using a time-varying noise estimator can be helpful in estimating the statistical characteristics of the uncertainty in the measurements in real time and mitigating the filter divergence. A further study on the SHAKF is developed based on adaptive factors for improving the filter practicability and optimality [23].

An adaptive fading Kalman filter (AFKF) was proposed for compensating the effect of the uncertainty in the measurements by transitive factor to the state error covariance (P). In AFKF, the state error covariance (P) is scaled with a single transitive factor for improving the filter variance and gain correction. When it is used for complex systems, the performance of AFKF degrades because of it may not be sufficient to use a single transitive factor for estimating the covariance matrix of the filter [24]. To overcome the difficulties of single transitive factor, multiple fading factors are used in AFKF. Because of that reason, authors are developed double transitive factor based SHAFKF that adapts both predicted state error covariance (P) and measurement noise covariance matrix (R) based on the innovation sequence. Although it has been successively applied to different domains, its performance for MEMS gyroscope sensor signal is not explored. The stochastic errors of MEMS gyroscope cannot be eliminated using calibration technique. It needs to be modeled before filtering the signal. Therefore, adaptive filtering techniques have been developed for minimizing the random noise from MEMS gyroscope system. In general, auto-regressive (AR), Moving Average (MA), and Auto-Regressive and Moving Average (ARMA) and Gauss-Markov model (GM) have been used for modeling stochastic signal [17]. Among these models, ARMA is a better choice for modeling MEMS gyroscope drift errors. In general, the ARMA modeling involved three steps as (i) randomness and stationary test (ii) selection of suitable time series model and (iii) estimation of model parameters. The unit root test and inverse sequence techniques have been used for checking the stationary of the signal. The model order is obtained by using auto correlation function (ACF) and partial auto correlation function (PACF). Moreover, Akaike Information Criterion also used to check the model order. The modified Yule-Walker method is used estimate the model parameters. Once an optimal ARMA model is defined, a suitable adaptive Kalman filter can be applied to minimize the drift of inertial sensors [14, 30].

In this chapter, we developed double transitive factors based on Sage-Husa adaptive fading Kalman filter (SHAFKF), namely SHAFKF-P Adaption and SHAFKF-R adaption. In addition, ARMA model is used to model the random drift errors of MEMS sensor. ARMA model based SHAFKF algorithm is developed and applied for minimizing the bias drift and random noise in the presence of MEMS gyroscope signal. The suggested algorithm is analyzed in two stages. In the first stage, the predicted state error covariance is adapted by a transitive factor, whereas, in the second stage, another transitive factor is scaled to the measurement noise covariance matrix (R). The efficiency of the algorithm is analyzed using Allan Variance technique.

The remainder of the paper is organized as follows. Section 2, explains the theory of ARMA models for MEMS gyroscope random noise analysis. The Allan Variance method is explained in Section 3. In Section 4, Conventional and adaptive Kalman filters are discussed based on innovation sequence. Section 5 explains the proposed algorithm based on double transitive factors. Designing state space model for ARMA (2, 1) model is presented in Section 6. Experimental results and static test analysis are explained in Section 7 and also followed by conclusions in Section 8.

2. Auto regressive and moving average (ARMA) model

In literature, several time series models have been widely used in many fields such as industry, science and engineering. Among the other model, auto regressive (AR) and moving average (MA) models have been most popular and since then widely used for forecasting [14, 15, 16]. The combination of AR and MA models has been used for inertial sensors error modeling. In this chapter, stationary ARMA model is proposed for characterizing the stochastic errors of the MEMS gyroscope signals. In general, the ARMA model is a combination of weighted sum of AR and MA model. The expression for the ARMA model with an order (p, q) is defined as

$$Y_n=\sum _{i=1}^p{\varnothing }_i{Y}_{n-i}\sum _{j=1}^q{\theta }_j{\epsilon }_{n-j}+{\epsilon }_n$$E1

where pand qare the AR and MA model orders, receptively. ${\epsilon }_n$is a sequence of independent and identical distributed random variable. $Y_n$is the measured time series data of MEMS gyroscope signal. ${\varnothing }_1$, ${\varnothing }_2$, ${\varnothing }_3$, …, ${\varnothing }_p$and ${\theta }_1$, ${\theta }_2$, ${\theta }_3$, …, ${\theta }_q$are the auto regressive (AR) and moving average (MA) coefficients, respectively. The MEME gyroscope sensor raw data is used to test the normality and zero mean of the time series data of MEMS gyroscope. In general, the skewness and Kurtosis should be 0 and 1 that tells that checking the zero mean and normal distributed data of the time series data of the sensors.

2.1 Time series model selection

In the time series analysis, several methods have been developed for selecting the order of the AR, MA and ARMA order. In general, auto-correlation function (ACF) and partial ACF (PACF) are the basic methods to select the model based on the characteristics of the ACF and PCF graphs as shown in Table 1. From Table 1, we observed that both ACF and PACF are tail off. In this chapter, ARMA (p, q) is suitable for modeling the MEMS Gyroscope data.

Model order	ACF	PACF
AR(p)	Tail off	Cut of at order P
MA(q)	Cut of at order q	Tail off
ARMA(p, q)	Tail off	Tail off

Table 1.

Determining the model and order of the MEMS gyro signal.

The samples autocorrelation function (ACF) is defined as

$$\mathit{ACF}=g_k=\frac{\frac{1}{N}{\sum }_{n=1}^{N-k}\left(Y_n-{\mu }_y\right)\left(\begin{array}{cc}{Y}_{n+k}-& {\mu }_y\end{array}\right)}{\frac{1}{N}{\sum }_{n=1}^N{\left(\begin{array}{cc}{Y}_n-& {\mu }_y\end{array}\right)}^2}$$E2

and the partial autocorrelation is expressed as

$$\mathit{\text{PACF}}=g_{\mathit{kk}}=\left\{\begin{array}{c}{g}_1\mathit{\text{if}}k=1\\ \frac{g_k{\sum }_{j=1}^{k-j}\left(g_{k-1}\right)\left(g_{k-j}\right)}{1-{\sum }_{j=1}^{k-j}\left(g_{k-1}\right)\left(g_k\right)}\mathit{\text{if}}k=2,3,\dots ,n,\end{array}\right.$$E3

where $k$is the lag and $g_k$is the sample autocorrelation. The ${\mu }_y$and $g_{\mathit{kk}}$are the samples mean and partial correlation at lag $k$.

This can also cross checked using Akaike Information Criterion (AIC) method. In this work, AIC values of the time series data are evaluated using Table 1. The model order is selected based on the minimum value of AIC.

The general expression of Akaike information criterion (AIC) is

$$\mathit{AIC}=log\left(\mathrm{\Theta }\right.\left[\begin{array}{cc}1+& \frac{2d_{\mathit{aic}}}{N_{\mathit{aic}}}\end{array}\right])$$E4

where $\mathrm{\Theta }$denotes the estimated residual of the model. $d_{\mathit{aic}}$and $N_{\mathit{aic}}$are the model order and the number of time series observation respectively.

2.2 Model parameter estimation

Suitable model parameters are estimated by using Yule-Walker, Burg, Unconstrained Least-Squares method and Levinson-Durbin methods. In general, for large data-set analysis, Yule-Walker and Unconstrained Least-Squares method are the better estimators.

3. Allan variance analysis

Allan variance (AV) is a popular time domain method has been widely used for identifying and quantifying random errors in the presence of inertial sensor [14]. Cluster based analysis is used to develop the AV technique. In the AV analysis, the IMU raw data can be divided into clusters with specified length, “$m$.” Let us take “$n$” measurements of gyroscope $\left(\omega \right)$, denote it by ${\omega }^{x\left[1\right]},{\omega }^{x\left[2\right]},{\omega }^{x\left[3\right]},\dots \dots {\omega }^{x\left[n\right]}$. The collected MEMS sensor data is sampled at rate of $\mathit{fs}$(samples per seconds). The set of samples called as cluster and denoted as “k_c”. The total number of clusters can be represented by “K,” (i.e., $K=\frac{n}{m}$). The measured date of the gyroscope can be written as

$${\omega }^{x\left[1\right]},{\omega }^{x\left[2\right]},{\omega }^{x\left[3\right]},\dots \dots {\omega }^{x\left[m\right]},{\omega }^{x\left[m+1\right]},\dots \dots {\omega }^{x\left[2m\right]},\dots \dots {\omega }^{x\left[n-m\right]}\dots {\omega }^{x\left[n\right]}$$

To calculate each clusters average is

$${\overline{\omega }}^{x⌈k_c⌉}=\sum _{i=1}^m{\overline{w}}^{x\left[k_c-1\right]m+i}$$E5

Here, $k_c$=1, 2, 3, …, K is the number of clusters.

The Allan variance is computed from two successive cluster averages for the specified correlation time which is defined as:

$${\sigma }^2\left({\tau }_m\right)=\frac{1}{2\left(K-1\right)}\sum _{k_c=1}^{K-1}(\left({\overline{\omega }}^{x\left[k_c+1\right]}\left(m\right)-{\overline{\omega }}^{x\left[k_c\right]}\left(m\right)){}^2\right)$$E6

where $k_c$=1, 2, 3, …, K, and ${\tau }_m=\raisebox{1ex}{$m$}\!\left/ \!\raisebox{-1ex}{$f_s$}\right.$is averaged period (or specified correlation time). The AV can be computed as a function of correlation times versus the Allan deviation plot are shown in Figure 1. The different contribution error sources are carried out simply by examining the slope of the AV plot. To extract the information on a specific source of error from the AV plot.

Figure 1.

Allan variance log-log plot.

There is a unique relationship between the Allan Variance (time domain) and the PSD (frequency domain) of the random process as:

$${\sigma }^2\left({\tau }_m\right)=4\underset{0}{\overset{\infty }{\int }}{S}_{\mathrm{\Omega }}\left(f\right)\frac{{\mathit{sin}}^4\left(\mathit{\pi fT}\right)}{{\left(\mathit{\pi fT}\right)}^2}$$E7

where $S_{\mathrm{\Omega }}\left(f_s\right)$is the power spectral density (PSD) of the random process and $\frac{{\mathit{sin}}^4\left(\mathit{\pi fT}\right)}{{\left(\mathit{\pi fT}\right)}^2}$is the transfer function of PSD.

The different random noise processes are characterized at various frequencies that are fitted by the AV method. The root Allan variance with each correlation time and slope are computed and presented in Table 2.

Noisy type	Units	Slope	Root Allan variance
Quantization noise (QN)	$°/\mathit{sec}$	−1	${\sigma }_{\mathit{QN}}\left(\tau \right)=\frac{\sqrt{3}\mathit{QN}}{\tau }$
Angle random walk (ARW)	$°/\sqrt{\mathit{hr}}$	−1/2	${\sigma }_{\mathit{ARW}}\left(\tau \right)=\frac{\mathit{ARW}}{\sqrt{\tau }}$
Bias instability (BS)	$°/\mathit{hr}$	0	${\sigma }_{\mathit{BS}}\left(\tau \right)=0.668\mathit{Bs}$
Rate random walk (RRW)	$°/\sqrt[3]{\mathit{hr}}$	1/2	${\sigma }_{\mathit{RRW}}\left(\tau \right)=\mathit{RRW}\sqrt{\frac{\tau }{3}}$
Rate ramp (RR)	$°/{\mathit{hr}}^2$	1	${\sigma }_{\mathit{RR}}\left(\tau \right)=\mathit{RRW}\frac{\tau }{\sqrt{2}}$

Table 2.

Allan variance analysis results.

4. Adaptive Kalman filtering

4.1 Conventional Kalman filter

The application of conventional Kalman filter (CKF) for the MEMS gyroscope requires a prior knowledge of dynamic process and measurement models. In addition, the process and measurement noise of the MEMS gyroscope. Considering a linear dynamic system, the state and measurement equations can be written as

$$x_k=A{x}_{k-1}+B{u}_k+w_k$$E8

$$z_k=H{x}_k+v_k$$E9

where $x_k$is the state vector at epoch k; $A$is the state transition matrix; $w_k$is the system (process) noise; $z_k$is the observation (measurement) at epoch k; $H$represents the observation matrix; and $v_k$is the measurement noise. Let us assume that the process $w_k$and measurement noises $\left(v_k\right)$are the Gaussian white noise with zero mean and finite variance that means that these are statistically independent from each other, following properties can be satisfied:

$$\mathrm{E}\left\{w_k\right\}=0,\mathrm{E}\left\{v_k\right\}=0$$E10

$$\mathrm{E}\left\{w_k{w}_k^T\right\}=Q_k$$E11

$$\mathrm{E}\left\{v_k{v}_k^T\right\}=R_k$$E12

Basically, the Kalman Filtering estimation algorithm comprises two steps, namely prediction and updating equations. The main Kalman Filtering equations are given below.

Prediction equations can be expressed as

$${\widehat{x}}_k^{-}=A{\widehat{x}}_{k-1}$$E13

$$P_k^{-}=A{P}_{k-1}{A}^T+Q_k$$E14

In the above equations, $A$is the state transition matrix and $A^T$denotes the transpose of $A$. $P_k^{-}$and $Q_k$represents prediction state error covariance and process noise covariance matrix at epoch k.

In the linear Kalman filter, the measurement updated equations are

$$K_k=P_k^{-}{H}^T{\left(H{P}_k^{-}{H}^T+R\right)}^{-1}$$E15

$${\widehat{x}}_k={\widehat{x}}_k^{-}+K_k\left(z_k-H{\widehat{x}}_k^{-}\right)$$E16

$$P_k=\left(I-K_{k}H\right)P_k^{-}$$E17

where ${\widehat{x}}_k$is the estimated state, $K_k$is the gain matrix and $P_k$is the estimated of state vector. $R$and $I$are the measurement noise covariance matrix and identify matrix respectively.

4.2 Innovation based adaptive estimation adaptive Kalman filter (IAE-AKF)

CKF requires a prior knowledge of the measurement and dynamic process models of MEMS IMU. In practice, statistical noise models of the process and measurement models are varying with time because of that the CKF would deprive optimality. To address this divergence, an adaptive KF (AKF) is a better solution. In the AKF, the adaptation can be carried out using three ways: (a) varying Qby trial and error until a stable state is estimated with fixed R[20]; (b) varying Rby keeping Qfixed; (c) varying Qand Rsimultaneously [21]. In the IAE-AKF algorithm, we selected the second adaption method is that varying the measurement noise covariance matrix (R) by keeping Q fixed based on innovation sequence $V_k$.

The innovation sequence is defined as the difference between true measurements and predicated measurements that can assume to be additional information to the filter. The innovation sequence is a zero-mean white Gaussian noise sequence, defined as

$$V_k=z_k-H{\widehat{x}}_k^{-}$$E18

The weighted innovation $K_k\left(z_k-H{\widehat{x}}_k^{-}\right)$acts as a correction to the predicted estimation ${\widehat{x}}_k^{-}$to form ${\widehat{x}}_k$. By substituting the measurement model (5) in (14), we get $V_k=H(x_k-{\widehat{x}}_k^{-}$). By taking variance on both sides of this, the theoretical covariance matrix of $V_k$is

$$C_{V_k}=H{P}_k^{-}{H}^T+R_k$$E19

The optimal estimation of covariance matrix of innovation sequence using average window method can be expressed as

$${\widehat{C}}_{V_k}=\frac{1}{D_s}\sum _{j=j0}^{D_s}{V}_j{V}_j^T$$E20

where $V_j$is the innovation sequence, $D_s$is the window size, $j0=k-D_s+1$is the first epoch. If the window size is too small, the measurement estimation covariance can be noisy; on the other hand, the estimation of measurement covariance will be smoother. Usually, window size is chosen empirically for statistical smoothing.

The estimated measurement noise covariance based on innovation sequence is

$${\widehat{R}}_k={\widehat{C}}_{V_k}-H{P}_k^{-}{H}^T$$E21

where ${\widehat{R}}_k$is the estimated measurement noise covariance matrix, $H$is the observation matrix, $P_k^{-}$is the prediction state error covariance and ${\widehat{C}}_{V_k}$is the estimated covariance matrix of innovation sequence.

4.3 Sage-Husa adaptive Kalman filter (SHAKF)

Sage-Husa AKF (SHAKF) is another class of adaptive filtering that uses a time-varying noise statistical estimator to proceed recursively. It is also used to reduce the sensor noise in the presence of MEMS IMU signals [16]. The linear dynamical process and measurement model equations can be written in the Eqs. (4) and (5).

The expectation and the covariance matrices of $w_k$and $v_k$are written as.

$$\mathrm{E}\left\{w_k\right\}={\widehat{q}}_k$$E22

$$\mathrm{E}\left\{v_k\right\}={\widehat{r}}_k$$E23

$$\mathrm{E}\left\{w_k{w}_k^T\right\}={\widehat{Q}}_k$$E24

$$\mathrm{E}\left\{v_k{v}_k^T\right\}={\widehat{R}}_k$$E25

where ${\widehat{Q}}_k$and ${\widehat{R}}_k$are the initial estimated process and measurement noise covariance matrices, respectively.

The time-varying noise statistic recursive estimator is given by:

$${\widehat{r}}_{k+1}=\left(1-d_k\right){\widehat{r}}_k+d_k\left(z_k-H{\widehat{x}}_k^{-}\right)$$E26

$${\widehat{R}}_{k+1}=\left(1-d_k\right){\widehat{R}}_k+d_k\left(V_k{V}_k^T-H{P}_k^{-}{H}^T\right)$$E27

$${\widehat{q}}_{k+1}=\left(1-d_k\right){\widehat{q}}_k+d_k\left(x_k-A{\widehat{x}}_k\right)$$E28

$${\widehat{Q}}_{k+1}=\left(1-d_k\right){\widehat{Q}}_k+d_k\left(K_k{V}_k{V}_k^T{K}_k^T+P_k-A{P}_{k-1}{A}^T\right)$$E29

where $d_k=\left(1-b_k\right)/\left(1-b_k^{k+1}\right)$is the amnestic factor, value range between 0 and 1. The Kalman filtering output signal and Sage-Husa self-adaptive Kalman filtering output signal are expressed in terms of following equations.

Prediction equations as.

$${\widehat{x}}_{k-1}^{-}=A{\widehat{x}}_{k-1}+{\widehat{q}}_k$$E30

$$P_{k-1}^{-}=A{P}_{k-1}{A}^T+{\widehat{Q}}_{k-1}$$E31

Measurement updated equations are equations:

$$K_k=P_k^{-}{H}^T{\left(H{P}_k^{-}{H}^T+R_k\right)}^{-1}$$E32

$${\widehat{x}}_k={\widehat{x}}_k^{-}+K_k\left(z_k-H{\widehat{x}}_k^{-}\right)$$E33

$$P_k=\left(I-K_{k}H\right)P_{k-1}^{-}$$E34

Here, an innovation sequence can be written as

$$V_k=z_k-H{\widehat{x}}_k-{\widehat{r}}_k$$E35

The $K_k$is the Kalman updated gain. $R_k$and I are the measurement noise covariance matrix and identity matrix respectively.

5. Proposed: Sage-Husa adaptive fading Kalman filter (SHAFKF) based on double transitive factors

Adaptive estimation methods have been developed for improving the CKF performance [22, 23]. In the AKF, covariance matching techniques is used to estimate the covariance matrix of the innovation or residual by fixing the values of Q. By using a scale factor in the AKF hence the performance filter was improved for estimating the state error covariance and also it improves the variance of the predicted state. Further, adaptive fading Kalman filters have been developed for improving the filter performance by introducing multiple adaptive scaling factors [24]. In the proposed algorithms, adaptive transitive factors based linear Adaptive Kalman filter algorithm is proposed also used for improving the MEMS gyroscope performance [25, 30]. However, a limited work has been reported the use of transitive factors in ARMA model based Sage-Husa KF. The proposed algorithm is explained in two cascaded stages. The predicted state error covariance Pis adapted in the stage one, whereas in the second stage, the measurement noise covariance Ris adapted by another transitive factor. The proposed scheme is shown in Sections 5.1 and 5.2, respectively.

5.1 Stage one: adaptation of predicted state error covariance (P)

In this stage, the predicted state error covariance is modified using an adaptive transitive factor. This stage is also termed as SHAFKF-P adaptation. The transitive factor is used to reduce the process noise of kinematic model based on the residual sequence.

The transitive factor $a_1\left(\mathrm{k}\right)$is evaluated as

$$a_1\left(\mathrm{k}\right)=\left\{\begin{array}{ll}1,& \mathit{tr}\left(C_{V_k}\right)>\mathit{tr}\left({\widehat{P}}_{\overline{v}k}\right)\\ \frac{\mathit{tr}\left({\widehat{C}}_{V_k}-R_k\right)}{\mathit{tr}\left(C_{V_k}-R_k\right)},& \mathit{\text{Otherwise}}\end{array}\right.$$E36

where tris the trace function and ${\widehat{P}}_{\overline{v}k}$is the estimated covariance matrix of the residual sequence expressed as

$${\widehat{P}}_{\overline{v}k}={\overline{V}}_k{\overline{V}}_k^T$$E37

The predicted state covariance ${\widehat{P}}_k^{-}$is updated as

$${\widehat{P}}_k^{-}=\frac{1}{a_1\left(k\right)}{\widehat{P}}_{k-1}^{-}$$E38

The SHAFKF-P adaptation algorithm, the predicted and estimated state error covariance are updated based on the SHAKF algorithm.

$${\widehat{P}}_k^{-}=\frac{1}{a_1\left(k\right)}{\widehat{P}}_{k-1}^{-}$$E39

$$C_{V_k}=H{P}_k^{-}{H}^T+R_k$$E40

The suboptimal state and update the measurement equations as

$$K_k=P_k^{-}{H}^T{\left(H{P}_k^{-}{H}^T+R_k\right)}^{-1}$$E41

$${\widehat{x}}_k={\widehat{x}}_k^{-}+K_k\left(z_k-H{\widehat{x}}_k^{-}\right)$$E42

$$P_k=\left(I-K_{k}H\right)P_{k-1}^{-}$$E43

where ${\widehat{x}}_k$indicates the suboptimal estimated state vector and $C_{V_k}$denotes the suboptimal covariance matrix of innovation the state vector. For optimal filter purpose, ${\widehat{x}}_k$and $P_k$are further passed to the next stage. The flowchart of the stage one of the algorithm is shown in Figure 2.

Figure 2.

Flow chart of the SHAFKF-P adaptation algorithm.

5.2 Stage two: adaptation of measurement noise covariance matrix (R)

The stage one algorithm requires prior knowledge of the state error vector and kinematic of model errors. To overcome this drawback and to eliminate the influence of the measurement noise disturbances, another transitive factor is introduced for updating the measurement noise covariance matrix (R). This stage is also termed as AUFKF-R adaptation.

In this stage, modified residual sequence is evaluated as the difference between measurement vector zkand the suboptimal estimated state (${\widehat{x}}_k$) evaluated using Eq. (45). Thus the modified residual sequence can be defined as

$${\overline{V}}_k=z_k-H{\widehat{x}}_k-{\widehat{r}}_k$$E44

Furthermore, using the suboptimal state error covariance ${\widehat{P}}_{\overline{v}k}$similar to Eq. (41), the estimated covariance matrix of the residual sequence can be written as

$$C_{\overline{v}k}=H{P}_k^{-}{H}^T+R_k$$E45

The suboptimal estimation of covariance matrix of residual sequence using the average window method is

$${\widehat{C}}_{\overline{v}k}=\frac{1}{D_s}\sum _{j=j0}^k{\overline{V}}_k\left(j\right){\overline{V}}_k^T\left(j\right)$$E46

The transitive factor $a_2\left(\mathrm{k}\right)$for the stage two is evaluated as

$$a_2\left(\mathrm{k}\right)=\left\{\begin{array}{ll}1,& \mathit{tr}\left(C_{\overline{v}k}\right)>\mathit{tr}\left({\widehat{C}}_{\overline{v}k}\right)\\ \frac{\mathit{tr}\left(C_{\overline{v}k}\right)}{\mathit{tr}\left({\widehat{C}}_{\overline{v}k}\right)},& \mathit{\text{Otherwise}}\end{array}\right.$$E47

In this algorithm, the measurement noise covariance matrix is scaled by a factor $a_2\left(\mathrm{k}\right)$. Thus Eq. (48) can be rewritten as

$$C_{\overline{v}k}=H{P}_k^{-}{H}^T+a_2\left(\mathrm{k}\right)R_k$$E48

The Kalman gain and state equations are updated as Eqs. (41)–(46). In this algorithm, measurement noise covariance matrix is multiplied by the adaptive transitive factor, $a_2\left(\mathrm{k}\right)$. If $a_2\left(\mathrm{k}\right)$large, $R_k$becomes larger, this helps to reduce the influence of uncertain measurement noise [23, 24]. The flow chart of the stage two, i.e., SHAFKF-R adaptation, is shown in Figure 3.

Figure 3.

Flow chart of the SHAFKF-R adaptation algorithm.

6. Designing state space model for ARMA (2, 1) model

The ARMA (p, q) model order is obtained using AIC method as in Table 3. The minimum values of AIC can be decided the optimal order of the ARMA (2, 1) is chosen. The ARMA (2, 1) model parameters such as Φ1 = −0.5422, Φ2 = −0.1204and θ1 = 0.1382are estimated based on the minimum AIC value, i.e., −5.7612. The parameters are tabulated in Table 3.

Model	${\mathit{\phi }}_{\mathbf{1}}$	${\mathit{\phi }}_{\mathbf{2}}$	${\mathit{\theta }}_{\mathbf{1}}$	AIC
AR(1)	−0.5422			−5.1769
AR(2)	−0.5422	−0.1204		−5.1832
MA(1)			−0.1382	−5.1860
ARMA(1, 1)	−0.5422		−0.1382	−5.5726
ARMA(2, 1)	−0.5422		−0.1382	−5.7612

Table 3.

Parameters estimation results with AIC values.

The ARMA (2, 1) model is used to approximate the MEMS Gyro sensor as:

$$Y_n={\phi }_1{Y}_{n-1}+{\phi }_2{Y}_{n-2}+{\theta }_1{\epsilon }_{n-1}+{\epsilon }_n$$E49

where Φis the AR coefficients and θis the MA model parameter, ${\epsilon }_n$is the system input Gaussian white noise with zero mean and variance ${\sigma }_n^2$. State-space representation of the optimal ARMA (2, 1) model is described as

$$X_k=\left[\begin{array}{cc}{\phi }_1& {\phi }_2\\ 1& 0\end{array}\right]\left[\begin{array}{c}{Y}_{n-1}\\ Y_{n-2}\end{array}\right]+\left[\begin{array}{cc}1& {\theta }_1\\ 0& 0\end{array}\right]W_k$$E50

$$Z_k=\left[\begin{array}{cc}1& 0\end{array}\right]\left[\begin{array}{c}{Y}_n\\ Y_{n-1}\end{array}\right]+V_k$$E51

where $W_k$= ${\left[\begin{array}{cc}{\epsilon }_k& {\epsilon }_{k-1}\end{array}\right]}^T$is the process noise. The initialize the state estimate ${\widehat{x}}_0$= ${\left[\begin{array}{cc}0& 0\end{array}\right]}^T$and state error covariance, ${\widehat{P}}_0$= Iare selected. In practice, the process noise covariance matrix and the measurement noise covariance matrix are assumed as $\left[\begin{array}{cc}{\sigma }_{W_k}^2& 0\\ 0& {\sigma }_{V_k}^2\end{array}\right]$. In the CKF, the process and measurement noise covariance matrices are constant whereas in the adaptive proposed algorithms, these parameters are changed iteratively.

7. Test results and discussion

The experimental setup consists of a single axis a prototype Xsens MTi 10 series MEMS sensor, turn table control unit, data acquisition board, and data processing computer. The MEMS gyroscope specification and test conditions of three single axis gyro sensor detailed results are reported in [36]. The experimental raw data is collected for 1 hour duration with sampling frequency at 100 Hz at room temperature. In the static condition, MEMS gyro is in zero rotation under the room temperature, for a more detailed specification of the Xsens MTi 100 series MEMS please refer to [36, 37].

7.1 Static performance test analysis

Three single-axis MEMS gyro sensor raw data are collected for 1 hour duration with sampling frequency at 100 Hz. The pre-processing methods are required to test the zero mean values for the sensor raw data before analyzing the Allan variance (AV) results [16]. Three single-axes of the MEMS Gyro sensor signals and corresponding AV results are plotted in Figure 4a and b respectively. From these figures, we see that the −1/2 slope indicates the angle random walk (ARW), which is a white noise characteristics. Bias instability (Bs) is due to internal and external electronic components of the sensor and is indicated at zero slope in log-log AV plot [16]. The three axes of MEMS IMU sensors are identified and quantified using AV analysis, which are presented in Table 2. From this table, we can observe that ARW and BI are the two that dominate noises in the presence of the MEMS sensor.

Figure 4.

(a) Three single axes of MEMS gyroscope raw signals and (b) corresponding Allan variance plot.

Conventional Kalman filter (CKF) algorithm is applied for minimizing the all three axis MEMS gyro static signal. In this experiment, the initial values of measurement and process noise covariance matrix are chosen as 0.098 and 0.0001 respectively. In practical application, these noise covariance matrices vary with time. In real-time, by adjusting the noise parameters are critical. The adaptive KF algorithm, an innovation sequence is used to adjust the noise parameters of process and measurement noise matrices ad it is followed by covariance matching principle. In the IAE-AKF algorithm, the window width selection is critical and can decide the filter optimality. In general, the window width is varied between 5 and 30. In this analysis, we observed that 15 samples of the window width is the optimal choice for statistical smoothing.

In the SHAKF algorithm, the innovation sequence is used to estimate the measurement noise covariance matrix and followed by sliding window average method. In addition, statistical noise estimator is used in the AKF frame work for updating the noise coefficients in each iteration recursively. The window width is 15 samples for statistical smoothing. The SHAKF algorithm results are plotted in Figures 5a–7a, respectively.

Figure 5.

(a) X-axis MEMS gyro signal and de-noised results using the SHAFKF algorithm and (b) corresponding Allan variance plot.

Figure 6.

(a) Y-axis MEMS gyro signal and de-noised results using the SHAFKF algorithm and (b) corresponding Allan variance plot.

Figure 7.

(a) Z-axis MEMS gyro signal and de-noised results using the SHAFKF algorithm and (b) corresponding Allan variance plot.

In the proposed approach, the predicted state error covariance is updated by one transitive factor whereas the measurement noise covariance matrix is updated using another transitive factors based on the residual sequence. The covariance matrix of residual sequence is estimated using sliding average window method. In this method, window width is chosen empirically as 15. In the first stage of the proposed algorithm (SHAFKF-P adaption), the transitive factor ($a_1$) is calculated in stage one. The measurement noise covariance matrix is scaled by an adaptive transitive factor ($a_2$) is in the second stage. The transitive factors are used to scale $R_k$and reciprocal to ${\widehat{P}}_k^{-}$for reducing the variance of uncertainty in the process model and measurements, respectively. The developed algorithm is also applied to X, Y and Z-axis MEMS gyroscope static signal. The test results of the proposed algorithm for X, Y and Z-axis data are shown in Figures 5a–7a, respectively. From these figures, it is observed that the angle random walk (ARW) and bias instability (Bs) noise are the dominated noise sources. The quantified noise coefficients are tabulated in the Tables 4–6, respectively. All the random noise and drift are quantified before and after applying the de-noising algorithm. The drift is also calculated before and after de-noising MEMS signal and tabulated in Tables 4–6, respectively. From these tables, it is observed that the ARW is reduced by 1000 and also Bs random noise is minimalized by order of 100 compared to the original value.

Methods	ARW ($\mathbf{°}/\sqrt{\mathit{hr}}$)	BS ($\mathbf{°}/\mathit{hr}$)	Drift ($\mathbf{°}/\mathit{hr}$)
MEMS raw data	165.115	8.775	1.758
CKF	103.235	7.459	1.362
IAE-AKF	24.858	3.496	0.859
SHAKF	4.228	2.296	0.0014
SHAFKF-P Adaption	1.279	0.690	0.00038
SHAFKF-R Adaption	0.331	0.421	0.00012

Table 4.

Allan variance and drift results of X-axis MEMS gyro using proposed scheme in static condition.

Methods	ARW ($\mathbf{°}/\sqrt{\mathit{hr}}$)	BS ($\mathbf{°}/\mathit{hr}$)	Drift ($\mathbf{°}/\mathit{hr}$)
MEMS raw data	33.0437	4.7297	1.7587
CKF	30.6510	1.5762	1.3624
SHKF	4.0098	1.2862	0.0859
IAE-SHAKF	3.4075	0.4371	0.00386
SHAFKF-P adaption	0.6570	0.2653	0.000562
SHAFKF-R adaption	0.442	0.150	0.000312

Table 5.

Allan variance and drift results of Y-axis MEMS gyro using proposed scheme in static condition.

Methods	ARW ($\mathbf{°}/\sqrt{\mathit{hr}}$)	BS ($\mathbf{°}/\mathit{hr}$)	Drift ($\mathbf{°}/\mathit{hr}$)
MEMS raw data	38.9222	9.8105	1.758
CKF	24.3805	8.3228	1.216
SHKF	7.0068	3.502	0.597
IAE-SHAKF	3.3229	3.071	0.0013
SHAFKF-P adaption	1.2516	0.832	0.00028
SHAFKF-R adaption	0.914	0.542	0.00014

Table 6.

Allan variance and drift results of Z-axis MEMS gyro using proposed scheme in static condition.

From these tables, it is evident that SHAFKF of R adaptation using transitive factor improves the performance of the algorithm. In this proposed algorithm, measurement noise covariance is scaled by the transitive factor. It ensures the variance is inversely proportional to the uncertainty of measurement. Due to this, SHAFK-R adaptation algorithm outperforms other algorithms.

In addition, we observed the Drift error for the MEMS gyroscope signals. Drift error is considered as one of the performance indicator of all the proposed algorithms. From Tables 4–6, it is observed that the proposed SHAFKF-R adaptation filter performs better than CKF, IAE-AKF SHAKF, and SHAFKF-P adaptation filters because of that the measurement noise covariance tunes by the adaptive transitive factor $a_2\left(\mathrm{k}\right)$to reduce the influence of uncertainty in measurement noise of the sensor.

8. Conclusions

In this chapter, the MEMS gyroscope drift is modeled by using ARMA (2, 1) for characterizing the MEMS gyro noise behavior. Moreover, ARMA-based linear Sage-Husa adaptive fading Kalman filter with double transitive factors is proposed. In the proposed algorithm, double adaptive transitive factors are used to update in the predicted state vector and measurement noise covariance matrix. The suggested algorithm is used to reduce the drift and random noise in the presence of MEMS gyroscope. From the AV analysis, the noise terms of ARW and Bs are reduced by order of 100. The proposed SHAFKF outperforms the CKF, IAE-AKF, and SHAKF algorithms in static case. It concludes that the SHAFKF algorithm is suitable for MEMS gyroscope signal drift minimization.