A Combined Method for MEMS Gyroscope Error Compensation Using a Long Short-Term Memory Network and Kalman Filter in Random Vibration Environments

Abstract

In applications such as carrier attitude control and mobile device navigation, a micro-electro-mechanical-system (MEMS) gyroscope will inevitably be affected by random vibration, which significantly affects the performance of the MEMS gyroscope. In order to solve the degradation of MEMS gyroscope performance in random vibration environments, in this paper, a combined method of a long short-term memory (LSTM) network and Kalman filter (KF) is proposed for error compensation, where Kalman filter parameters are iteratively optimized using the Kalman smoother and expectation-maximization (EM) algorithm. In order to verify the effectiveness of the proposed method, we performed a linear random vibration test to acquire MEMS gyroscope data. Subsequently, an analysis of the effects of input data step size and network topology on gyroscope error compensation performance is presented. Furthermore, the autoregressive moving average-Kalman filter (ARMA-KF) model, which is commonly used in gyroscope error compensation, was also combined with the LSTM network as a comparison method. The results show that, for the x-axis data, the proposed combined method reduces the standard deviation (STD) by 51.58% and 31.92% compared to the bidirectional LSTM (BiLSTM) network, and EM-KF method, respectively. For the z-axis data, the proposed combined method reduces the standard deviation by 29.19% and 12.75% compared to the BiLSTM network and EM-KF method, respectively. Furthermore, for x-axis data and z-axis data, the proposed combined method reduces the standard deviation by 46.54% and 22.30% compared to the BiLSTM-ARMA-KF method, respectively, and the output is smoother, proving the effectiveness of the proposed method.

1. Introduction

Fiber optic gyroscopes and laser gyroscopes have excellent performance, but they are too large and expensive for portable devices [1,2]. Micro-electro-mechanical-system (MEMS) gyroscopes have, in recent years, been used in low-cost inertial navigation systems (INS) due to their small size and low cost. However, the MEMS gyroscope has a significant error due to the manufacturing technology and structural composition [3,4]. The error of the MEMS gyroscope can be divided into deterministic error and random error. The deterministic error mainly refers to perturbation errors such as zero offsets and the scale factor, which can be corrected by a calibration test [5,6]. Random error refers to the random drift caused by uncertain factors, usually determined by the device’s accuracy level [7], with no precise repeatability. Therefore, it is difficult to accurately compensate for random error, which hinders the further improvement of MEMS gyroscope performance.

In MEMS gyroscope error compensation research, the MEMS gyroscope data are generally treated as time-series data. Scholars have proposed methods such as the autoregressive moving average (ARMA) model, the Allan variance (AV), the wavelet threshold (WT), the support vector machine (SVM), and the artificial neural network (ANN), and all of them have achieved excellent results [7,8,9,10,11,12,13,14]. Recently, various variants of the recurrent neural network (RNN), which has strong processing power for time-series data, have been shown to be superior to traditional methods in the research of error compensation in MEMS gyroscopes [15,16,17,18].

However, most of the research mentioned above has acquired data by placing the MEMS gyroscope in a static environment. In practical applications of the MEMS gyroscope, it is inevitably that it is affected by random vibration [19]. In random vibration environments, the MEMS gyroscope is interfered with by both internal device noise and external vibration noise [20], which dramatically affects the performance of the MEMS gyroscope. The degradation of performance in vibrating environments is a fatal problem for MEMS gyroscopes [21,22], so it is essential to research error compensation methods in random vibration environments.

Most of the current research on improving the performance of the MEMS gyroscope in random vibration environments is to fix the MEMS gyroscope on a vibration isolation platform [23,24,25]. However, this kind of method is not universal [26]. There is not much research based on time-series models—the windowed measurement error covariance (WMEC) method has been applied to compensate for the effects of the vibration environments [27], singular spectrum analysis (SSA) was proposed to remove the low-frequency vibration noise perturbations of MEMS accelerometers [28], and the third-order autoregressive (AR) model was used to estimate the Kalman filter to compensate for the MEMS gyroscope’s attitude angle error caused by random vibration [29].

Considering the dramatic perturbation of the MEMS gyroscope in random vibration environments, in this paper, a combined method of a long short-term memory (LSTM) network and Kalman filter is proposed for error compensation, with the Kalman smoother and expectation-maximization (EM) algorithm to dynamically adjust the predicted values of the LSTM network to improve the performance in error compensation. The main contributions of this paper are as follows:

(1)

The combination of LSTM network and Kalman filter is applied to MEMS gyroscope error compensation in random vibration environments;

(2)

The proper input data step and the network topology are explored, and the error compensation performance of the bidirectional LSTM (BiLSTM) network and other recurrent neural network (RNN) variants are compared;

(3)

In designing the Kalman filter, the EM algorithm is used to estimate the parameters. It is compared with the ARMA model, a parameter estimation method commonly used in research of the MEMS gyroscope error compensation problem.

The remainder of this paper is organized as follows: (1) Section 2 introduces the methods, including BiLSTM network, Kalman filter, ARMA-KF model, and EM-KF model, and gives the illustration of this paper proposed method; (2) Section 3 presents the experiment, results, and comparisons; and (3) the remaining sections are the conclusion, appendix, and references.

2. Method

2.1. Multi-Layer BiLSTM Network and Kalman Filter

The long short-term memory network is a variant of the recurrent neural network used to solve the gradient vanishing or gradient explosion problem of RNNs [30,31]. A detailed description of LSTM units can be found in references [15,16,17,18].

The basic LSTM network only considers the historical and current inputs and ignores future inputs [32]. Therefore, the LSTM network can perform the reverse operation, superimpose the forward and reverse information flows, and fully utilize the front and back inputs at the current time to improve the error compensation performance. In addition, the previous hidden layer’s output is used as the input of the following layer to explore the more in-depth features of the time-series data, thus enhancing the model’s nonlinear fitting ability. The multi-layer BiLSTM network information flow is shown in Figure 1.

The cell state of the Layer n BiLSTM network at time t can be presented as:

$$\left[\begin{array}{c}{i}_t^{\left(n\right)}\\ f_t^{\left(n\right)}\\ o_t^{\left(n\right)}\\ {\tilde{c}}_t^{\left(n\right)}\end{array}\right]=\left[\begin{array}{c}\sigma \\ \sigma \\ \sigma \\ tanh\end{array}\right] \left[\begin{array}{cc}{W}_{i,x}^{\left(n\right)}& W_{i,h}^{\left(n\right)}\\ W_{f,x}^{\left(n\right)}& W_{f,h}^{\left(n\right)}\\ W_{o,x}^{\left(n\right)}& W_{o,h}^{\left(n\right)}\\ W_{\tilde{c},x}^{\left(n\right)}& W_{\tilde{c},h}^{\left(n\right)}\end{array}\right] \left[\begin{array}{c}{h}_t^{\left(n-1\right)}\\ h_t^{\left(n\right)}\end{array}\right]$$

(1)

where $h_t^{\left(n-1\right)}$ is the hidden state of the Layer $n-1$ at time $t$. Each hidden state is composed of forward and reverse superposition. The related equations are denoted as follows:

$$h_t^{\left(n\right)}={\overrightarrow{h}}_t^{\left(n\right)}\oplus {\overleftarrow{h}}_t^{\left(n\right)}$$

(2)

The Kalman filter is an optimal state estimation method that can be applied to dynamic systems with random disturbances. It estimates the system state based on discrete measurement that contain noise [33,34]. Suppose the state–space model is built as:

$${\widehat{x}}_k=\Phi {\widehat{x}}_{k-1}+\Gamma {\omega }_{k-1}$$

(3)

$$y_k=H{x}_k+v_k$$

(4)

where $\Phi$ is the system state transition matrix, $\Gamma$ is the system noise-driven matrix, $H$ is the measurement matrix, ${\widehat{x}}_k$ is the system state vector, $y_k$ is the measurement vector, ${\omega }_k$ is the system noise vector, and $v_k$ is the measurement noise vector.

The noise of the system models and measurement models are assumed to have normal distribution in the Kalman filter, such that [35]:

$${\omega }_k\sim \mathcal{N}\left(0,Q\right)$$

(5)

$$v_k\sim \mathcal{N}\left(0,R\right)$$

(6)

where $Q$ is the covariance matrix of the system models and $R$ is the covariance matrix of measurement models.

The Kalman filter is composed of two-stage optimization. In the predicted stage, the current system state vector is predicted based on the system state vector at the previous time, such that:

$${\widehat{x}}_{k/k-1}=\Phi {\widehat{x}}_{k-1}$$

(7)

$$P_{k/k-1}=\Phi P_{k-1}{\Phi }^T+\Gamma Q{\Gamma }^T$$

(8)

where ${\widehat{x}}_{k/k-1}$ is the predicted value of the system state vector and $P_{k/k-1}$ is the predicted covariance matrix of the system state vector.

In the updated stage of the Kalman filter, the current system state vector is updated by using the measurement vector, such that:

$$K_k=P_{k/k-1}{H}^T{\left(H{P}_{k/k-1}{H}^T+R\right)}^{-1}$$

(9)

$${\widehat{x}}_k={\widehat{x}}_{k/k-1}+K_k\left(y_k-H{\widehat{x}}_{k/k-1}\right)$$

(10)

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

(11)

where $K_k$ is the Kalman filter gain matrix and $P_k$ is the updated covariance matrix of the system state vector.

2.2. Kalman Filter Design with ARMA Model

The autoregressive moving average model is the most widespread model used in time-series analysis, and it is derived and developed on the basis of the linear regression model. The ARMA model can be described as [36]:

$$x_t={\displaystyle \sum }_{i=1}^p{\phi }_i{x}_{t-i}-{\displaystyle \sum }_{j=1}^q{\theta }_j{\epsilon }_{t-j}+{\epsilon }_t, {\epsilon }_t\sim \mathcal{N}\left(0,{\delta }_{\epsilon }^2\right)$$

(12)

where $p$ and $q$ are the order of the ARMA model; ${\phi }_i$ and ${\theta }_j$ are coefficients that satisfy stationary and invertible conditions, respectively [37]; $\mathrm{and} {\epsilon }_t$ is white noise, which is an uncorrelated random variable with mean zero and constant variance. The model expresses that the measured values of the stochastic process $\left\{x_t\right\}$ at time $t$ are correlated with the previous $p$ measurements and the previous $q$ white noise.

The steps to design a Kalman filter using the ARMA model as follows: (1) test the stationarity and normality of the measurement data, (2) determine the model type according to the autocorrelation function and partial autocorrelation function, (3) determine the order and parameters of the model according to the Akaike information criterion (AIC) [38], and (4) perform adaptive testing of the designed model.

2.3. Kalman Filter Design with EM Algorithm

The expectation-maximization algorithm is an iterative method proposed by Shumway and Stoffer to compute maximum likelihood estimates based on incomplete data [39]. It is convergent and can identify parameters and states in the model [40]. Andrieu and Doucet introduced the EM algorithm for parameter estimation for linear state–space models was introduced [41]. The EM algorithm is an iterative numerical algorithm for computing the maximum likelihood estimation (MLE). The linear Gaussian state–space model used for the EM algorithm can be expressed as follows [42]:

$$x_k=\Phi x_{k-1}+{\omega }_{k-1}$$

(13)

$$y_k=H{x}_k+{\upsilon }_k$$

(14)

The conditional probability densities of the state equation and the measurement equation are obtained from Equations (13) and (14), respectively:

$$p\left(x_k|x_{k-1}\right)=\exp \left\{-\frac{1}{2}{\left[x_k-\Phi x_{k-1}\right]}^T{Q}^{-1}\left[x_k-\Phi x_{k-1}\right]\right\}{\left(2\pi \right)}^{-\frac{n}{2}}{\left|Q\right|}^{-\frac{1}{2}}$$

(15)

$$p\left(y_k|x_k\right)=\exp \left\{-\frac{1}{2}{\left[y_k-H{x}_k\right]}^T{R}^{-1}\left[y_k-H{x}_k\right]\right\}{\left(2\pi \right)}^{-\frac{m}{2}}{\left|R\right|}^{-\frac{1}{2}}$$

(16)

It is assumed that the likelihood of the system state data and the evolution of the states is Gaussian, which are defined by the following equations [43]:

$$p\left(Y,X|\Theta \right)=p\left(x_1\right){\displaystyle \prod }_{k=1}^{N}p\left(y_k|x_k\right){\displaystyle \prod }_{k=2}^{N}p\left(x_k|x_{k-1}\right)$$

(17)

where $Y$ is the measurement data, $X$ is the unknown system state data, and $\Theta$ is the parameter set of linear Gaussian state–space model. $\Theta$ can be represented as follows:

$$\Theta =\left\{\Phi ,H,Q,R\right\}$$

(18)

By taking the log of the likelihood we arrive at the following formula:

$$\ln p\left(Y,X|\Theta \right)=-\frac{1}{2}{\displaystyle \sum }_{k=2}^N\left\{\ln \left|Q\right|+{\left[x_k-\Phi x_{k-1}\right]}^T{Q}^{-1}\left[x_k-\Phi x_{k-1}\right]\right\}-\frac{1}{2}{\displaystyle \sum }_{k=1}^N\left\{\ln \left|R\right|+{\left[y_k-H{x}_k\right]}^T{R}^{-1}\left[y_k-H{x}_k\right]\right\}$$

(19)

Depending on the maximum likelihood method, the linear Gaussian state–space model can be identified through an EM algorithm [42]. The algorithm alternates between two steps—the E-step (expectation) and the M-step (maximization) [44]. In general, the likelihood density function based on the measurement data, denoted by $p\left(\Theta |Y\right)$, is called the posterior distribution of the measurement. The EM algorithm aims to compute the maximum likelihood estimation of $p\left(\Theta |Y\right)$. ${\Theta }_i$ is denoted as the estimate of the likelihood function at the beginning of the $i$th iteration.

In the E-step, the expectation for the conditional distribution of $\ln p\left(Y,X|\Theta \right)$ concerning $X$, is calculated such that:

$$\mathsf{\Omega }\left(\Theta |{\Theta }_i,Y\right)\stackrel{\wedge}{=}{E}_X\left\{\ln p\left(\Theta |Y,X\right)|{\Theta }_i,Y\right\}={\displaystyle \int }\left[\ln p\left(\Theta |Y,X\right)\right]\left(X|{\Theta }_i,Y\right)dX$$

(20)

In the M-step, $\mathsf{\Omega }\left(\Theta |{\Theta }_i,Y\right)$ is maximized to find ${\Theta }_{i+1}$ such that:

$$\mathsf{\Omega }\left({\Theta }_{i+1}|{\Theta }_i,Y\right)=arg\underset{\Theta }{\max }\left[\mathsf{\Omega }\left(\Theta |{\Theta }_i,Y\right)\right]$$

(21)

The E-step and the M-step are iterated until,

$$\Vert L\left({\Theta }_{i+1}\right)-L\left({\Theta }_i\right)\Vert <\tau$$

(22)

where $\tau$ is the predefined threshold. Equation (22) means that it has satisfied the convergence criterion. The specific process of designing a Kalman filter using the EM algorithm is given as follows:

The value of $\mathsf{\Omega }\left(\Theta |{\Theta }_i,Y\right)$ is determined by the following [45]:

$$E_X\left(x_k|Y\right)={\widehat{x}}_{k|N}$$

(23)

$$E_X\left(x_k{x}_{k-1}^T|Y\right)=P_{k,k-1|N}+{\widehat{x}}_{k|N}{\widehat{x}}_{k-1|N}^T$$

(24)

$$E_X\left(x_k{x}_k^T|Y\right)=P_{k|N}+{\widehat{x}}_{k|N}{\widehat{x}}_{k|N}^T$$

(25)

where ${\widehat{x}}_{k|N}$ is the smoothed value of the system state vector and $P_{k|N}$ is the smoothed covariance matrix of the system state vector. $P_{k,k-1|N}$ is initialized by:

$$P_{k,k-1|k}=\left(I-K_{k}H\right)\mathsf{\Phi }{P}_{k-1}$$

(26)

$$P_{k,k-1|N}=P_{k,k-1|k}+\left[P_{k|N}-P_k\right]P_{k|k}^{-1}{P}_{k,k-1|k}$$

(27)

${\widehat{x}}_{k|N}$ and $P_{k|N}$ can be obtained by smoothing the outputs of Kalman filter using backward-pass methods such as the Rauch–Tung–Striebel (RTS) smoother [46]. This method is summarized in the following equations:

$$J_{k-1}=P_{k-1}{\Phi }^T{P}_{k|k-1}^{-1}$$

(28)

$${\widehat{x}}_{k-1|N}={\widehat{x}}_{k-1}+J_{k-1}\left({\widehat{x}}_{k|N}-{\widehat{x}}_{k/k-1}\right)$$

(29)

$$P_{k-1|N}=P_{k-1}-J_{k-1}\left(P_{k|N}-P_{k|k-1}\right)J_{k-1}^T$$

(30)

Then, the model parameters are re-estimated by maximizing the $\mathsf{\Omega }\left(\Theta |{\Theta }_i,Y\right)$ over $\Theta$ using partial derivatives of $\mathsf{\Omega }\left(\Theta |{\Theta }_i,Y\right)$ and setting them to zero. Solving these equations yields the updated parameters (in the $i$th iteration) as follows:

$$\frac{\partial L(\Theta )}{\partial \Phi }=-{\displaystyle \sum }_{k=2}^N{Q}^{-1}\left(P_{k,k-1|N}+{\widehat{x}}_{k|N}{\widehat{x}}_{k|N}^T\right)+{\displaystyle \sum }_{k=2}^N{Q}^{-1}\Phi \left(P_{k-1|N}+{\widehat{x}}_{k-1|N}{\widehat{x}}_{k-1|N}^T\right)=0$$

(31)

$${\Phi }_{i+1}=\left({\displaystyle \sum }_{k=2}^N{P}_{k,k-1|N}+{\widehat{x}}_{k|N}{\widehat{x}}_{k|N}^T\right){\left({\displaystyle \sum }_{k=2}^N{P}_{k-1|N}+{\widehat{x}}_{k-1|N}{\widehat{x}}_{k-1|N}^T\right)}^{-1}$$

(32)

$$\frac{\partial L(\Theta )}{\partial H}=-{\displaystyle \sum }_{k=1}^N{R}^{-1}{y}_k{\widehat{x}}_{k|N}^T+{\displaystyle \sum }_{k=1}^N{R}^{-1}H\left(P_{k|N}+{\widehat{x}}_{k|N}{\widehat{x}}_{k|N}^T\right)=0$$

(33)

$$H_{i+1}=\left({\displaystyle \sum }_{k=1}^N{y}_k{\widehat{x}}_{k|N}^T\right){\left[{\displaystyle \sum }_{k=1}^N\left(P_{k|N}+{\widehat{x}}_{k|N}{\widehat{x}}_{k|N}^T\right)\right]}^{-1}$$

(34)

$$\frac{\partial L(\Theta )}{\partial Q^{-1}}=\frac{N}{2}Q-\frac{1}{2}{\displaystyle \sum }_{k=1}^N\left(P_{k|N}+{\widehat{x}}_{k|N}{\widehat{x}}_{k|N}^T\right)+\Phi \left[\frac{1}{N}{\displaystyle \sum }_{k=2}^N\left(P_{k,k-1|N}+{\widehat{x}}_{k|N}{\widehat{x}}_{k-1|N}^T\right)\right]=0$$

(35)

$$Q_{i+1}=\frac{1}{N}\left({\displaystyle \sum }_{k=1}^N\left(P_{k|N}+{\widehat{x}}_{k|N}{\widehat{x}}_{k|N}^T\right)-{\Phi }_{i+1}{\displaystyle \sum }_{k=2}^N\left(P_{k,k-1|N}+{\widehat{x}}_{k|N}{\widehat{x}}_{k-1|N}^T\right)\right)$$

(36)

$$\frac{\partial L(\Theta )}{\partial R^{-1}}=\frac{N+1}{2}R-{\displaystyle \sum }_{k=1}^N\left(\frac{1}{2}{y}_k{y}_k^T-H{\widehat{x}}_{k|N}{y}_k^T+\frac{1}{2}H\left(P_{k|N}+{\widehat{x}}_{k|N}{\widehat{x}}_{k|N}^T\right)H^T\right)=0$$

(37)

$$R_{i+1}=\frac{1}{N+1}\left({\displaystyle \sum }_{k=1}^N{y}_k{y}_k^T\right)-H_{i+1}{\left(\frac{1}{N+1}{\displaystyle \sum }_{k=1}^N{y}_k{\widehat{x}}_{k|N}\right)}^T$$

(38)

In this paper, based on the EM algorithm, the proposed LSTM and Kalman filter combination method is illustrated in Figure 2.

3. Experiments and Results

In this section, the designed experiments and the analysis of the results are presented to verify the effectiveness of the proposed method.

3.1. Data Acquisition

The MSI320H MEMS Inertial Measurement Unit (IMU) was employed for experiments. This consists of a three-axis MEMS gyroscope and a three-axis MEMS accelerometer. The real picture and the gyroscope specifications of MSI320H are shown in Figure 3a and

Table 1. Specifications of MSI320H gyroscope.

GYRO	Input range	±1800°/s
	Bias instability (Allan variance)	36°/h
	Angular random walk (Allan variance)	0.4°/√h
	Bandwidth (−3 dB)	≥220 Hz
GENERAL	Sample rate	100$\sim$1000 Hz
	Weight	≤25 g
	Supply voltage	5.0 ± 0.5 V
	RS422 transmission bit rate	921600 bps
	Mechanical shock, any direction	≥20,000 g

, respectively. The MSI320H was fixed on the vibration table. A picture of the vibration table is shown in Figure 3b. The data acquisition procedure of the MSI320H is shown in Figure 3c. Data from the MSI320H was sent to the xPC via the RS422 communication interface with a Baud of 921,600 bps. The xPC decoded the gyroscope data and sent it to the host computer via the network cable. The MSI320H was preheated at room temperature with power for 20 minutes. Then, linear vibration experiments were performed. The vibration direction of the vibration table is the y-axis of the gyroscope, and the power spectral density (PSD) of the linear random vibration loads is shown in Figure 3d.

As illustrated in Figure 3d, the acceleration of the applied vibration loads can be expressed as follows:

$$a_v\left(t\right)={\xi }_v\sin \left({\omega }_{v}t\right), {\omega }_v\in \left[20·2\pi ,2000·2\pi \right]$$

(39)

where ${\xi }_v=6$ g. The sample rate was set to 200 Hz, and approximately 60,000 data were acquired in the random vibration environment. The variation of gyroscope x-axis signal affected by random vibration is shown in Figure 4, and the error of the gyroscope increased significantly in random vibration environments.

3.2. Comparison of BiLSTM and Other RNN Variants

The outliers of the acquired data in random vibration environments were eliminated using the Puata criterion [47]. Because the linear random vibration test’s direction was the y-axis of the gyroscope, the x-axis and z-axis directions were in a random vibration environment. For the consideration of model generality, we took the last 80% of the processed x-axis data as the training set and the first 20% of the x-axis and z-axis data were used as the testing set.

If the hidden layer structure of the designed network is too simple, it is not easy to characterize the time-series model of gyroscope data. Conversely, it increases the complexity of the network, reduces the learning speed of the network, and tends to fall into local minima during the learning process. With the above considerations, in this paper, the proposed BiLSTM network is shown in Figure 5. The dense layer transformed the high-dimensional stacked sequence of the hidden layer into an output sequence with the same shape as the input sequence. Moreover, considering a large number of network parameters, dropout was set in the dense layer to prevent the overfitting phenomenon [48].

The adaptive moment estimation (Adam) optimization algorithm was used to update the network parameters [49]. The Adam algorithm uses default parameters. The activation function of the Dense layer is rectified linear unit (ReLU). The specifications used for network training are illustrated in

Table 2. Specifications used for network training.

The output dimension of dense layer	1
Activation function of dense layer	ReLU
Dropout rate	0.5
Batch size	256
Training epoch	50
Learning rate	0.001

. Moreover, the root mean square error (RMSE) was used as the loss function. The network training was performed on Tensorflow 2.0.0 and Keras 2.3.1 over the Ubuntu 16.04-LTS-x86 64 operating system. The heterogeneous computing platform was equipped with Intel Xeon E5-1620 and GeForce RTX-2080Ti GPUs.

In order to verify the performance of BiLSTM for gyroscope error compensation in random vibration environments, proper values for the input data step size, the number of hidden units, and the number of hidden layers was first explored using the x-axis testing set. Subsequently, training was performed using the identified values. The BiLSTM network results were compared with the LSTM network, gated recurrent unit (GRU) network, and bidirectional GRU (BiGRU) network using the x-axis and z-axis testing sets, respectively.

As shown in

Table 3. Performance with varying values of the input data step.

Number of Hidden Layers	Number of Hidden Units	Input Data Step	STD (°/s)	Time/Epoch
10	64	5	0.1551	23 s
10	64	10	0.1481	38 s
10	64	15	0.1483	60 s
10	64	20	0.1346	82 s
10	64	25	0.1368	98 s
10	64	30	0.1501	115 s

,

Table 4. Performance with varying values of the number of hidden units.

Number of Hidden Layers	Number of Hidden Units	Input Data Step	STD (°/s)	Time/Epoch
10	8	20	0.1504	82 s
10	16	20	0.1459	82 s
10	32	20	0.1559	81 s
10	64	20	0.1346	82 s
10	128	20	0.1326	81 s
10	256	20	0.1468	88 s

and

Table 5. Performance with varying values of the number of hidden layers.

Number of Hidden Layers	Number of Hidden Units	Input Data Step	STD (°/s)	Time/Epoch
1	128	20	0.1493	10 s
2	128	20	0.1513	19 s
3	128	20	0.1597	27 s
4	128	20	0.1557	34 s
5	128	20	0.1658	42 s
6	128	20	0.1505	50 s
7	128	20	0.1470	58 s
8	128	20	0.1459	66 s
9	128	20	0.1542	75 s
10	128	20	0.1326	81 s
11	128	20	0.1405	90 s
12	128	20	0.1393	98 s

, when the input data step size and the number of hidden layers are more extensive, the training time per epoch will be longer. So we needed to make a trade-off between the results and the computational performance. According to the results, the best results were obtained when taking the input data step of 20, the number of hidden units of 128, and the number of hidden layers of 10. Although it does not indicate that this is the optimal parameter for the network, it will be the proper value to be obtained considering the computational resources.

The results are shown in Figure 6 and Figure 7 and

Table 6. Comparison of raw data, BiLSTM, LSTM, BiGRU, and GRU standard deviation values for x-axis data.

x-axis	STD (°/s)	Percentage
Raw data	0.2493	$-$
BiLSTM	0.1326	53.19%
LSTM	0.1543	61.89%
BiGRU	0.1501	60.21%
GRU	0.1604	64.34%

and

Table 7. Comparison of raw data, BiLSTM, LSTM, BiGRU, and GRU standard deviation values for z-axis data.

z-axis	STD (°/s)	Percentage
Raw data	0.2400	$-$
BiLSTM	0.1353	56.38%
LSTM	0.1550	64.58%
BiGRU	0.1504	62.67%
GRU	0.1574	65.58%

. Figure 6 shows the training loss within 50 epochs, and convergence was achieved for all networks. Table 6 and Table 7 show that the standard deviations of the BiLSTM network results for the x-axis and z-axis were reduced by 46.81% and 43.63%, respectively, compared to the raw data, proving that the BiLSTM network is feasible for the application in the research of MEMS gyroscope error compensation. Furthermore, compared with the results of the LSTM network, the BiGRU network, and the GRU network, the standard deviation values of BiLSTM results in the x-axis were reduced by 14.06%, 11.66%, and 17.33%, respectively, and the standard deviation values of BiLSTM results in the z-axis were reduced by 12.71%, 10.04%, and 14.04%, respectively. This indicates that the error compensation performance of the BiLSTM network is better than these three networks.

3.3. Comparison of LSTM-EM-KF and LSTM-ARMA-KF

In this section, the raw data and BiLSTM network results of the x-axis and z-axis are used as the measurement, respectively. The Kalman filter parameters are estimated by the ARMA model and EM algorithm, and the filter results are compared.

3.3.1. Estimating Kalman Filter Parameters Using the ARMA Model

When modeling time-series data using the ARMA model, the time-series data must meet stationarity and normality requirements. Therefore, a polynomial fitting method was used to eliminate the trend term before modeling. In this paper, the stationarity was tested using the run test, and the normality was tested by calculating the skewness, $\xi ,$ and the kurtosis, $\upsilon$.

According to the test, after eliminating the trend term, the raw data and BiLSTM network results of the x-axis and z-axis met the stationarity and normality requirements. The test process is shown in Appendix A. Moreover, as illustrated in Figure A1, the autocorrelation function and partial autocorrelation function diagrams exhibit trailing properties, and all models can be identified as an ARMA $\left(p,q\right)$ model.

The next step is to determine the order of the model. If the order is increased, the identified model will be more realistic, but the computational difficulty will also increase as the order increases. Therefore, the maximum order was set to 3, which means the maximum value of $p$ and $q$ was set to 3. Furthermore, the Akaike information criterion (AIC) was used for determining model order. Determining the model order process and the Durbin–Watson test results are shown in Appendix B.

For x-axis raw data, the model identified is identified as ARMA(3,3):

$$x_n=-0.3126x_{n-1}+0.8168x_{n-2}+0.1520x_{n-3}+{\epsilon }_n+0.6885{\epsilon }_{n-1}-0.3241{\epsilon }_{n-2}-0.0366{\epsilon }_{n-3}$$

(40)

For x-axis BiLSTM network results, the model is identified as ARMA(3,3):

$$x_n=0.8505x_{n-1}+0.8891x_{n-2}-0.7745x_{n-3}+{\epsilon }_n+0.1354{\epsilon }_{n-1}-0.8476{\epsilon }_{n-2}-0.0815{\epsilon }_{n-3}$$

(41)

For z-axis raw data, the model identified is identified as ARMA(1,3):

$$x_n=0.8593x_{n-1}+{\epsilon }_n-0.51915{\epsilon }_{n-1}+0.0328{\epsilon }_{n-2}-0.0236{\epsilon }_{n-3}$$

(42)

For z-axis BiLSTM network results, the model identified is identified as ARMA(3,3):

$$x_n=2.0000x_{n-1}-1.2268x_{n-2}+0.2031x_{n-3}+{\epsilon }_n-1.0446{\epsilon }_{n-1}-0.0809{\epsilon }_{n-2}-0.1086{\epsilon }_{n-3}$$

(43)

where $x_n$ is the output of the ARMA model, ${\epsilon }_n$ is the driving white noise (with mean, 0, and variance, ${\widehat{\delta }}_{\epsilon }^2$). The Kalman filter parameters are presented in

Table 8. Kalman filter parameters of all measurements.

	$\mathbf{\Phi }$	${\Gamma }$	$\mathit{H}$	$\mathit{Q}$	$\mathit{R}$
x-axis raw data	$\left[\begin{array}{cccc}-0.3126& 0.8168& 0.1520& 0\\ 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right]$	$\left[\begin{array}{cccc}1& 0.6885& -0.3241& -0.0366\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$	$\left[\begin{array}{cccc}1& 0& 0& 0\end{array}\right]$	$\left[\begin{array}{cccc}0.0370& 0& 0& 0\\ 0& 0.0370& 0& 0\\ 0& 0& 0.0370& 0\\ 0& 0& 0& 0.0370\end{array}\right]$	$0.0604$
x-axis BiLSTM	$\left[\begin{array}{cccc}0.8505& 0.8891& -0.7745& 0\\ 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right]$	$\left[\begin{array}{cccc}1& 0.1354& -0.8476& -0.0815\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$	$\left[\begin{array}{cccc}1& 0& 0& 0\end{array}\right]$	$\left[\begin{array}{cccc}0.0036& 0& 0& 0\\ 0& 0.0036& 0& 0\\ 0& 0& 0.0036& 0\\ 0& 0& 0& 0.0036\end{array}\right]$	$0.0175$
z-axis raw data	$\left[\begin{array}{cccc}0.8593& 0& 0& 0\\ 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right]$	$\left[\begin{array}{cccc}1& -0.5191& 0.0328& 0.0236\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$	$\left[\begin{array}{cccc}1& 0& 0& 0\end{array}\right]$	$\left[\begin{array}{cccc}0.0405& 0& 0& 0\\ 0& 0.0405& 0& 0\\ 0& 0& 0.0405& 0\\ 0& 0& 0& 0.0405\end{array}\right]$	$0.0596$
z-axis BiLSTM	$\left[\begin{array}{cccc}2.0000& -1.2268& 0.2031& 0\\ 1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right]$	$\left[\begin{array}{cccc}1& -1.0446& 0.0809& 0.1086\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$	$\left[\begin{array}{cccc}1& 0& 0& 0\end{array}\right]$	$\left[\begin{array}{cccc}0.0043& 0& 0& 0\\ 0& 0.0043& 0& 0\\ 0& 0& 0.0043& 0\\ 0& 0& 0& 0.0043\end{array}\right]$	$0.0183$

. The value of R is the variance of the measurement. The initial value of the Kalman filter is set as follows: $x_1=\left[0;0;0;0\right]$, and $P_1$ is the fourth-order identity matrix.

3.3.2. Estimating Kalman Filter Parameters Using the EM Algorithm

When using the EM algorithm to estimate the Kalman filter parameters, only the iteration convergence conditions and initial parameters need to be set. The M-step convergence constant $\tau$ in Equation (22) was set to 0.1. The Kalman filter’s initial values were set to $x_1=0 \mathrm{and} P_1=1$, and the Kalman filter’s initial parameters were set to ${\Phi }_1=1$, $H_1=1$, $Q_1=1$, and $R_1=1$. The change of the log-likelihood function during the iteration of the EM algorithm is shown in Figure 8. The parameter estimation results are presented in

Table 9. The parameter estimation results.

	$\mathbf{\Phi }$	$\mathit{H}$	$\mathit{Q}$	$\mathit{R}$
x-axis raw data	0.9723	0.1350	0.0016	0.1181
x-axis BiLSTM	0.9738	0.3422	0.0016	0.2852
z-axis raw data	0.9471	0.1327	0.0065	0.1155
z-axis BiLSTM	0.9530	0.2983	0.0044	0.2430

.

3.3.3. Kalman Filtering Results

The results are illustrated in

Table 10. Comparison of raw data, BiLSTM, ARMA-KF, EM-KF, BiLSTM-ARMA-KF, and BiLSTM-EM-KF standard deviation values for x-axis data.

x-axis	STD (°/s)	Percentage
Raw data	0.2493	$-$
BiLSTM	0.1326	53.19%
ARMA-KF	0.1618	64.90%
EM-KF	0.0943	37.83%
BiLSTM-ARMA-KF	0.1201	48.17%
BiLSTM-EM-KF	0.0642	25.75%

and

Table 11. Comparison of raw data, BiLSTM, ARMA-KF, EM-KF, BiLSTM-ARMA-KF, and BiLSTM-EM-KF standard deviation values for z-axis data.

z-axis	STD (°/s)	Percentage
Raw data	0.2400	$-$
BiLSTM	0.1353	56.38%
ARMA-KF	0.1853	77.21%
EM-KF	0.1098	45.75%
BiLSTM-ARMA-KF	0.1233	51.38%
BiLSTM-EM-KF	0.0958	39.92%

. For the x-axis data, the standard deviation of the BiLSTM-EM-KF results was reduced by 51.58% and 31.92% compared to the BiLSTM network and EM-KF, respectively. For the z-axis data, the standard deviation of the BiLSTM-EM-KF results was reduced by 29.19% and 12.75% compared to the BiLSTM network and EM-KF, respectively. Therefore, the combined method proposed in this paper can be demonstrated to improve the gyroscope error compensation performance of the BiLSTM network and EM algorithm. Moreover, compared with BiLSTM-ARMA-KF results, the standard deviation of the BiLSTM-EM-KF was reduced by 46.54% and 22.30% in x-axis and z-axis, respectively. It indicates the proposed method’s superior performance to that of BiLSTM-ARMA-KF. Furthermore, according to Figure 9, the curves of BiLSTM-EM-KF results are smoother, which proves that the proposed combined method is effective.

4. Conclusions

In this paper, a combined method of an LSTM network and Kalman filter is proposed for MEMS gyroscope error compensation in random vibration environments. Through the results, the following conclusions were obtained:

(1)

After exploring proper input data step size and network topology, the network was trained, and the test results showed that the BiLSTM network outperformed the LSTM network, the GRU network, and the BiGRU network in gyroscope error compensation;

(2)

Combining the BiLSTM network with the EM-KF method can improve their gyroscopic error compensation performance;

(3)

In the classical gyroscope error compensation method, the ARMA-KF method, tedious data testing and model checking are required. In contrast, the EM-KF method only needs to set the initial parameters and the convergence value, which is much easier to apply. Moreover, the ARMA-KF method parameters cannot be updated through the filtering process, which means that satisfactory results cannot be obtained if the parameters are not defined correctly before the filtering process. From the filtering results, compared with BiLSTM-ARMA-KF, the standard deviation of the BiLSTM-EM-KF results were 46.54% and 22.30% lower, in x-axis and z-axis, and the output curve was smoother, which proves the effectiveness of the proposed method in this paper.

Future work should include conducting dynamic field experiments to obtain MEMS gyroscope outputs, as well as combining neural networks with more state-of-the-art Kalman filter methods for MEMS gyroscope error compensation and using fiber optic gyroscopes or laser gyroscopes as benchmarks for comparison.