MEMS IMU stochastic error modelling

ABSTRACT Low cost inertial measurement units (IMU) are comprised of micro-electro-mechanical systems (MEMS) gyroscopes and accelerometers. These inertial sensors have high random noise and time varying bias. Hence, an accurate stochastic error model is necessary to predict performance and also to implement an attitude and heading reference system or a navigator. The parameters of this stochastic model are classically obtained via Allan Variance analysis. In this paper, two modern approaches, based on autoregressive moving average model and Kalman Filter, respectively, are investigated and their performances are unveiled via sensitivity analysis, realistic simulations with typical MEMS parameters and also with experimental data. Basic equations for the estimation problem are developed and the solution sensitivity to data and parameters is discussed. It is shown that the deficiency with the online methods can be traced back to the difficulty in estimating the parameters B (bias instability) and K (rate random walk for gyros, acceleration random walk for accelerometers) separately under high measurement noise N (angular random walk for gyros, velocity random walk for accelerometers), which is typical of MEMS IMU.


Introduction
Due to its low cost, micro-electro-mechanical systems inertial measurement units (MEMS IMU) has been extensively used in the implementation of inertial navigation systems (INS). See, for instance, Kingston and Beard (2004), Li, Dempster, Li, Wang, and Rizos (2006) and Grewal, Weil, and Andrews (2007). The accuracy of these systems depends on the adequate error modelling of the IMU sensors. This error model is actually divided into two parts: (a) a deterministic or systematic model and (b) a stochastic model. The deterministic model encapsulates the biases, scale factors and misalignments errors. The parameters associated to the deterministic error model can be obtained by two different approaches: (1) static procedure in a laboratory environment, such as in Wang, Hao, Wei, and Wang (2010) and (2) identification based procedure, with a cost function which incorporates the unknown parameters. See Skog and Händel (2006) for details. In both cases, if the gyroscope does not measure the Earth's angular velocity, as is usually the case with MEMS sensors, then a rotating table is required for calibration.
Once the deterministic IMU model is obtained, the remaining error has stochastic nature and must be modelled accordingly. For a detailed definition of all the CONTACT Elder M. Hemerly hemerly@ita.br parameters in this stochastic model, see IEEEStd 952-1997IEEEStd 952- (1998, which also recommends the Allan Variance method for experimentally obtaining these parameters. This model can then be used, for instance, to implement INS-based navigators, such as in Hasan, Samsudin, and Ramli (2011) and Wang, Wang, Liang, Zhang, and Zhou (2013). A comprehensive example of the classical approach for stochastic error modelling, based on Allan Variance analysis, is Petkov and Slavov (2010). Experimental works conducted there indicated that MEMS gyroscope stochastic error model is mainly composed of: (a) bias instability (B), which represents a time varying bias modelled by a first-order Gauss-Markov process; (b) angular random walk (N), which is an angular error due to white noise in the angular rate and (c) rate random walk (K), representing a rate error due to white noise in angular acceleration. These parameters are also the most relevant for accelerometers, except for the change in physical meanings: N is velocity random walk and K is the acceleration random walk. In Vaccaro and Zaki (2012), the statistical properties of the estimates obtained via Allan Variance analysis are carefully investigated, by considering tacticalgrade MEMS gyros. These sensors are easier to tackle than those considered by Petkov and Slavov (2010), since an accurate model can be obtained with only 2 parameters: N (angular random walk) and K (rate random walk).
Despite its broad application to IMU stochastic error modelling, the classical approach, based on Allan Variance analysis, is not an automatic procedure. It does require some interaction with the user, for determining adequate asymptotes to the Allan graphics: accordingly to the IEEE standard (IEEEStd 952-1997(IEEEStd 952- , 1998, each parameter in the error model has a distinctive signature in the Allan graphics. For instance, the parameters N and K are associated with asymptotes having slopes −1/2 and 1/2, respectively. Some procedures have then been proposed to estimate automatically the parameters associated to the IMU stochastic error model. In Seong, Lee, and Park (2000), it is shown that a combination of several different noises can be described by a single ARMA model, by properly calculating the parameters of the equivalent correlated noise. It is argued that for ring laser gyroscope (RLG), by considering only N (angular random walk) and Q (quantization noise), the performance is similar to the Allan approach and the implementation is simpler.
In Saini, Rana, and Kuber (2010), the equivalent ARMA model representation proposed by Seong et al. (2000) is employed as departing point, but the input/output model is rewritten in the state space form. MEMS sensors are considered, instead of RLG; hence there are four parameters associated to the stochastic error model: bias instability (B), angular random walk (N), rate random walk (K) and rate ramp noise (R). The problem then boils down to a state estimation problem, where the dynamic matrix A, the output matrix C and the measurement intensity noise are unknown. The solution is implemented via the Expectation Maximization (EM) algorithm, which solves the ensuing Maximum Likelihood Estimation problem. It is concluded that the proposed online method avoids the storage of data and produce results within the error limits of Allan Variance analysis.
It should be noted that the number of parameters to be retained in the IMU stochastic model depends on the fabrication technology: Petkov and Slavov (2010) considers three parameters, Vaccaro and Zaki (2012), Seong et al. (2000) only two and Saini et al. (2010) employs four parameters. In general, the better the technology, the lesser the number of parameters required in the model.
In this work, MEMS gyros similar to those investigated in Petkov and Slavov (2010) are considered, that is, the most relevant parameters in the stochastic model are B, N and K. The goal is to investigate if the modern approaches for IMU stochastic modelling, such as those proposed by Seong et al. (2000) and Saini et al. (2010), do provide reliable alternative to the classical approach base on Allan Variance analysis, such as in Petkov and Slavov (2010). As a first step, the approaches in Seong et al. (2000) and Saini et al. (2010) are tailored to the three parameters model employed by Petkov and Slavov (2010).
The main contributions are: (1) detailed sensitivity analysis of the estimates to the data and underlying parameters, which was not carried out in Seong et al. (2000) and Saini et al. (2010); (2) comparative performance evaluation of classical and modern approaches, by using realistic simulations with typical MEMS gyros parameters, and experimental data. This has not been done by Petkov and Slavov (2010), Vaccaro and Zaki (2012), Seong et al. (2000), Saini et al. (2010) and (3) light is shed on the online methods deficiency in the present case, by tracing it back to the difficulty in estimating the parameters B and K separately under the typical high measurement noise N exhibited by MEMS IMU.

Equivalent ARMA model approach
In this paper, this method is applied to the MEMS IMU modelling, with the three main parameters B, N and K, already discussed in Section 1. However, at first the simpler case of Seong et al. (2000) is considered, where only the parameters N (angle random walk) and Q (quantization noise) is revisited, for correcting a typo error and preliminary performance evaluation. In order to obtain the equivalent ARMA model in this simpler case, recall that in the discrete time domain the contributions of the error sources in the gyro measurement are, respectively, where the standard deviations σ 1 and σ 2 depend on the parameters N and B as indicated in Saini et al. (2010). Then, by defining the measurement as the application of the equivalent ARMA model representation to the right hand side of (3), which comprises two different noises, produces the equivalent model, with just one noise, which is an ARMA model with structure (0,1). In (4), the values of the parameters e 0 and e 1 are obtained by imposing equal correlation coefficients to the right hand sides in (3) and (4). This produces the relations σ 2 1 + 2σ 2 2 = e 2 0 + e 2 1 and − σ 2 2 = e 0 e 1 , with the constraint that the polynomial C(q −1 ) = e 0 + e 1 q −1 , where q −1 is the backward shift operator, is stable, that is, has root inside the unit circle. The solution of (5) presented in Equation (9) of Seong et al. (2000) is incorrect. The correct solution is The simulation from Seong et al. (2000) is now reproduced with the corrected relations (6). In this case, the standard deviations associated to the angle random walk and the quantization noise are, respectively, σ 1 = 35.6 deg/s and σ 2 = 114.525 deg/s, and the sampling time is T = 0.01 s. With these values, the true parameters in the equivalent ARMA model (4) are, accordingly to (6), e 0 = 133.7 and e 1 = −98.1. Table 1 displays the errors in per cent, as functions of the number of readings, obtained from a Monte Carlo simulation. In these simulations, the ARMA model identification was carried out by using the armax() command from Matlab, since the identification problem is not the main issue in this paper.
From Table 1 it can be concluded that the equivalent ARMA modelling is quite efficient for the two parameters (N and Q) model: good performance is obtained even with a small number of measurements. The problem now is to see if this good performance is also exhibited for the MEMS IMU stochastic modelling, where there are three parameters (B, N and K), instead of two (N and Q, which is the case shown in Table 1). The problem here does not lie only on the number of parameters, but mainly on their physical meaning. See IEEEStd 952-1997 (1998) for details.

Equivalent ARMA model for MEMS IMU
Instead of (3), now the MEMS IMU model must comprise, as in Petkov and Slavov (2010), the contributions of bias instability, angle random walk and rate random walk, that is, the gyro measurement is given by In (8) the parameter α is obtained by discretizing the continuous time first-order Gauss-Markov process, that is and the standard deviations σ 1 , σ 2 and σ 3 are related to the parameters B, N and K as in Saini et al. (2010).
Since the rate random walk (10) is a non-stationary process, instead of taking (7) as being the gyro measurement, it is advisable to consider as measurement the difference With this proviso, from (7) to (10) and (12) the measurement equation is rewritten as the stationary process Hence the equivalent ARMA model has order (1,2), that is and by calculating the correlations, the relations between the standard deviations and the MA parameters in (14) are The MEMS IMU stochastic error modelling then boils down to the following problem: given the gyro measurements (12), estimate the parameters (e 0 , e 1 , e 2 ) in the ARMA model (14) and then get the standard deviations (σ 1 , σ 2 , σ 3 ) from (15). The constraints here are: (a) the parameter α must be positive, from (11), and (b) the polynomial C(q −1 ) = e 0 + e 1 q −1 + e 2 q −2 must be stable, that is, must have roots inside the unit circle.
The estimation of the parameter α in (15) comes directly from the identified ARMA model, due to (14). Hence, (15) corresponds to an algebraic system with three equations and three unknowns (σ 1 , σ 2 , σ 3 ), from which the parameters B, N and K of the stochastic error model are finally obtained.

Simulation results
Just to recall, given gyro measurements (12), the ARMA model (14) is identified, thereby providing estimates for α, e 0 , e 1 and e 2 . With these values, from (15) the standard deviations σ 1 , σ 2 and σ 3 are obtained, and then parameters which define the stochastic error model, namely B, N and K.
However, for performance analysis, it is better to proceed the other way round: generate simulated data with given values for B, N, K and α, then obtain e 0 , e 1 and e 2 from (15). These are the true parameters associated to the equivalent ARMA model. Hence, good performance is obtained if the ARMA identification algorithm provides estimates close enough to the true values of e 0 , e 1 and e 2 .
A typical MEMS gyro is considered, with parameters given by Petkov and Slavov (2010) as B = 1.55 × 10 −2 deg/s, N = 5.75 × 10 −2 deg/s 1/2 and K = 9 × 10 −4 deg/s 3/2 . By assuming correlation time equal to 100 s and sampling time equal to 0.01 s, the true values for e 0 , e 1 and e 2 can be obtained from (15) in the following way: remove e 0 and e 1 from the last two equations in (15) and substitute them into the first equation in (15). An algebraic equation of order 8 for e 2 is obtained, with even powers only. Hence, this equation can be reduced to order 4. A total of eight roots are produced, but only the positive solutions are of interest, and each one must be tested for the constraint that C(q −1 ) = e 0 + e 1 q −1 + e 2 q −2 is Hurwitz. Let e 0 , e 1 and e 2 be the viable solution.
Now, suppose that, as an ideal case, the identification algorithm produces exactly the true viable solution e 0 , e 1 and e 2 . Then, σ 2 is easily calculated by using the last equation in (15), and from the first two equations, σ 1 and σ 3 are given by From (16), the solutions for σ 1 and σ 3 depend critically on the matrix M condition number. In this particular example, β = 1/100 = 0.01 s and T = 0.01 s. From (11), α = 0.9999, which implies cond(M) = 10 9 . Hence, the system of Equations (16) is very badly conditioned and the parameters σ 1 and σ 3 cannot be determined separately with accuracy. This conclusion is completely independent of the identification algorithm accuracy, hence the number of gyro measurements. More precisely, increasing the number of readings will not change the scenario. Figure  1 shows the matrix M condition number as function of the gyro correlation time and the sampling time, both in seconds.
From Figure 1 it is concluded that matrix M is only well conditioned if the correlation time is small and the sampling time is large. This mathematical difficulty in obtaining σ 1 and σ 3 separately with accuracy can be physically explained by two different factors: (a) the rate random walk variance increases as the sampling time decreases and (b) for small values of β and T, from (11) it follows that α is close to 1; hence the first order Gauss-Markov process (8) does not differ too much from the random walk process (10). More precisely, the problem boils down to estimating two quite close parameters under high noise condition. Indeed, recall that for MEMS IMU the parameter N (angle random walk) is large; hence σ 2 is large in (9). Consequently, the measurement (12) is very noisy.
With the previously defined values of B, N, K, β and T, 100,000 readings were generated by using (7)-(12), which would correspond to an acquisition time of 16.167 minutes. By using the Matlab armax() command for identification, the only parameter which was correctly estimated was σ 2 , from which the estimateN = 0.0574 was obtained. This represents an error of only 0.212% with respect to the true value of N. The estimate of α waŝ α = −0.9783, which is close to 1, as expected, but the signal is incorrect. The values for B and K were incorrect by orders of magnitude. Hence, a comparison with the Allan   Variance analysis is necessary. The corresponding Allan graphic, for the same realization, is shown in Figure 2.
From Figure 2, the estimates obtained from Allan Variance analysis areN = 5.6 × 10 −2 ,B = 1.2 × 10 −2 and K = 1.1 × 10 −3 . Hence, the error in per cent are, respectively, 2.6, 22.5 and 18.2%. Therefore, the results from the Allan Variance analysis are good for N and reasonable for B and K, whereas the results from the equivalent ARMA modelling are very good for N and useless for B and K. It is then concluded that Allan Variance is a better approach in the present case.
Remark 2.1: A potential cause for the error in the estimation of α in the equivalent ARMA approach could be the presence of correlated noise in (14). Hence, trails were carried out by employing the instrumental variable identification algorithm, by means of the ivar() command from Matlab. However, the results were equally incorrect, perhaps due to the high noise content of the measurement signal (12), which is, as already pointed out, the usual case when MEMS gyros are employed.

Remark 2.2:
An additional attempt was made to improve the estimate of α by using a constrained identification scheme, since it is known from the outset that this parameter is positive, according to (11). Again, this procedure was not reliable, due to the high measurement noise.

Remark 2.3:
The sampling time was then increased to T = 1 s, which implies an acquisition time of 27.77 hours. The performance did not improve for the equivalent ARMA model. Several other simulations, varying the number of data points, sampling time and so on did not alter the conclusion: the equivalent ARMA model is not applicable for estimating simultaneously the parameters N, B and K associated to the stochastic error model of a MEMS gyro.

Simplified model
Due to the inadequate performance reported in Section 2.2, a particular case is now considered, by dropping the parameter K from the MEMS gyro stochastic error model. This case is still relevant and simplifies the problem considerably, since the nonstationarity caused by (10) is removed. With this simplification, the equivalent MA part of the ARMA model is as in Section 2, that is, C(q −1 ) = e 0 + e 1 q −1 , but instead of (5), the algebraic system of equations is now whose solution is much simpler than the system (15). By generating 100,000 readings with the true parameters values B, N, β and T, the MATLAB armax() command produced inadequate result for the estimate of α. The problem now is not caused by almost singularity in the algebraic equation system (17), but solely by the high measurement noise.
The sampling time was then increased from 0.02 to 0.1 s, for reducing the random walk variance. The error in per cent associated to the parameter β estimation, as function of the sampling time, is shown in Figure 3.
For T = 0.1 s, several realizations were considered and a typical estimation result wasN = 5.74 × 10 −2 ,B = 1.52 × 10 −2 andα = 0.9992, with in per cent errors 0.21, 2.17 and 0.02%, respectively. Hence, by considering simplified MEMS gyros model with only parameters B and N, it is concluded that the equivalent ARMA modelling is

Kalman filter approach
Another modern approach for IMU stochastic modelling is proposed by Saini et al. (2010). The departing point is the equivalent ARMA approach proposed by Seong et al. (2000), but Saini et al. (2010) formulates the problem in the continuous time domain and employs state space representation. By considering the model with parameters B, N, K and β as in Section 2, the equivalent ARMA model is given bÿ which can be rewritten in state space form and discretized via Euler method, with small enough sampling time T, as in Saini et al. (2010), according to where, for small enough T, the state noise covariance matrix is given by (14) in the sense that the rate random walk pole at z = 1 has not been removed by data differencing. This explains why (19) is a second order model, whereas (14) is a first-order system.

Remark 3.2:
Technically, the problem of estimating B, N, K and β given (19)-(20) constitutes an adaptive state estimation problem, that is, joint estimation of state and parameters. The main difficulty here arises from the following fact: in (20), the intensity of the measurement noise v(k) depends on the unknown parameter N. Hence, the measurement noise covariance matrix is not known. Therefore, the usual approach for implementing adaptive state estimation, namely the one base on state augmentation and extended Kalman Filter, cannot be applied directly.

Remark 3.3:
In Saini et al. (2010), the adaptive state estimation problem posed by (19)-(20) is solved by using the EM, see Saini et al. (2010) for details. Several sensitivities must be calculated, that is, the ensuing algorithm is much more complicated than the Kalman Filter.

Remark 3.4:
From all the cases considered in Section 2, that is, for MEMS IMU, the estimation of the angle random walk N can be easily accomplished, since it is the parameter which influences the output the most. As a matter of fact, a simple sample variance can provide accurate results. For instance, by using the same parameters B, N, K, β and T in Section 2.2, with 500, 1000 and 6000 data points the error in per cent associated to the estimation of N was 3.47, 3.02 and 0.63%, respectively. Hence, in 1 minute (6000*0.01 = 60 s) it is possible to obtain a good enough estimate for N. Indeed, the data acquired to estimate N has to correspond to a short amount of time; otherwise B and K can display their effects, thereby spoiling the estimate of N.
Based on Remarks 3.2-3.4, the following simplified and much simpler version of Saini et al. (2010) is investigated for estimating the parameters B, N, K and β associated to the MEMS gyro stochastic error model: Step 1: By using initial data, estimate N by calculating the sample mean variance.
Step 2: Given the previous estimate of N, the adaptive estimation problem posed by (19)-(20) can be solved by the Extended Kalman Filter with augmented state where x(k) is defined in (19). For implementing the EKF, it is also supposed that the unknown parameters B, K and β are modelled as constants.
With the previous assumptions, the linearized system to which the Kalman Filter is applied has dynamic and output matrices defined by Simulations where carried out by generating data with the same parameters B, N, K, β and T as in Section 2, and the performance was highly dependent on the realization. The explanation was found in the very badly conditioning of the observability matrix of the linearized model. Indeed, its condition number was of order 10 22 . Similar difficulty has already been encountered in Section 2.2, where it was concluded that the separate estimation of the parameters B and K is not reliable. In Saini et al. (2010) this observability issue is not even addressed.

Experimental results
Experimental data from a MEMS gyro was obtained, by using the mAHRS, from Innalabs. Only raw roll gyro data is considered, with 130,000 samples, acquired at 100 Hz. The Allan graphic is shown in Figure 4.
From Figure 4, the estimates for N, B and K obtained from Allan Variance analysis, by using the asymptote lines, areN = 2.95 × 10 −2 deg/s 1/2 ,B = 1.74 × 10 −2 deg/s andK = 5.2 × 10 −3 deg/s 3/2 . The estimate for N looks adequate, since the sample variance, calculated with the first 1000 data points, produces the rough estimate 2.41 × 10 −2 . These parameters have approximately the same magnitudes as those exhibited by the sensor employed in Petkov and Slavov (2010), hence are of the same MEMS class.
By using now the equivalent ARMA model as in Section 2.2, the estimated parameters turned out to be completely unreliable, for instanceN = 1.87 × 10 −1 deg/s 1/2 . At least for the estimate of N, this result can be easily explained as follows: (1) the Matlab armax() command produces a final prediction error (FPE) equal to 0.0625. Hence, if all the noise effects were to come only from N, the estimated value would beN = 2.5 × 10 −2 deg/s 1/2 , which is coherent with the valueN = 2.95 × 10 −2 deg/s 1/2 estimated via Allan Variance analysis; (2) in the equivalent ARMA modelling, the N estimate depends on σ 2 , which from (15) depends on the estimates of α, e 0 and e 2 . But these three parameters are incorrectly estimated, despite several attempts of providing different initial estimates for the ARMA polynomials, via the Matlab command idpol().
For the Kalman Filter approach, the estimate for N was obtained as in Remark 3.4, that is,N = 2.41 × 10 −2 deg/s 1/2 . Several tests were performed by varying Figure 4. Allan graphic and asymptote lines with inclinations +1/2, 0 and −1/2, for estimating N, B and K, respectively, by using experimental data from mAHRS-Innalabs roll gyro. the initial estimates for B and K. The values obtained via Allan Variance analysis were used as references. The conclusions were the following: (a) the estimates for B do not vary much around their initial estimates and (b) the estimates for K are more sensitive to the data, but converge to values which are dependent on the initial estimates. This suggests the experimental data are not exciting enough.

Conclusions
The stochastic error modelling problem for MEMS IMU has been addressed in this work, by investigating two recently proposed techniques (equivalent ARMA model and state estimation, respectively), which aims at replacing the classical procedure based on Allan Variance analysis. Typical parameters from MEMS gyros were used to generate realistic data, which were then used for performance evaluation. Sensitivity of the parameter estimates to data was carefully investigated. The findings were corroborated by extensive simulations, which indicated that for the general case of stochastic model with three parameters (B, N and K), the performances of both modern approaches are deficient, due to badly conditioned data. Physically, this difficulty arises from the necessity of separately estimating the parameters B and K, under high measurement noise caused by N. Experimental data were used, and the conclusions were the same.
If the parameter K is ignored in the MEMS IMU stochastic model, then the equivalent ARMA approach can provide adequate estimates, as shown in Section 2.3, if the sampling time is large enough. This procedure is simpler than the state estimation approach in Section 3.
For general MEMS IMU stochastic modelling, the classical approach based on the Allan Variance analysis is still a reliable one for estimating the four parameters B, N, K and β (associated to the parameter B by means of the first order Gauss-Markov model (8), as indicated by (11)). Additionally, since in this approach the user knows beforehand what are the expected asymptote slopes for each parameter in the stochastic error model, it is possible do decide which parameters should be kept in the model, given the experimental data. This important structural analysis is not straightforward in the other two approaches: (a) the order selection in the ARMA equivalent modelling is bound to be spoiled by bad data conditioning and (b) in the Kalman Filter approach, each structure demands a different filter, besides the observability problem.

Disclosure statement
No potential conflict of interest was reported by the author.