Denoising Method of MEMS Gyroscope Based on Interval Empirical Mode Decomposition

-emicroelectromechanical system (MEMS) gyroscope has lowmeasurement accuracy and large output noise; the useful signal is often submerged in the noise. A new denoisingmethod of interval empirical mode decomposition (IEMD) is proposed. Firstly, the traditional EMD algorithm is used to decompose the signal into a finite number of intrinsic mode functions (IMFs). Based on the Bhattacharyya distance analysis and the characteristics of the autocorrelation function, a screening mechanism is proposed to divide IMFs into three categories: noise IMFs, mixed IMFs, and signal IMFs. -en, the traditional modelling filtering method is used to filter the mixed IMFs. Finally, the mixed IMFs after modelling and filtering and signal IMFs are reconstructed to obtain the denoised signal. In the experimental analysis, the static denoising experiment of the turntable, the Allan variance analysis, dynamic denoising experiment, and vehicle experiment are set up in this paper, which fully proves that the method has obvious advantages in denoising and greatly improves the quality of signal and the accuracy of the inertial navigation system solution.


Introduction
Compared with other types of gyroscopes, MEMS gyroscopes based on micro electromechanical system technology have the advantages of miniaturization, low price, low power consumption, and easy installation [1]. MEMS gyroscope has been widely used in many fields because of the great development space and development value. However, due to the characteristics of the components and the external environment, the measurement accuracy of the MEMS gyroscope is relatively low, the output noise is large, and the useful signal is often submerged in the noise. erefore, in order to improve the measurement accuracy, it is very important to find an effective denoising method for MEMS gyroscope.
In recent years, domestic and foreign scholars have conducted a lot of in-depth research on suppressing the random drift of MEMS gyroscopes. Among them, wavelet transform, neural network modelling, time series modelling, and empirical mode decomposition are the main methods.
In reference [2], two wavelet basis functions with the best denoising effect are identified and selected; the multi-stage filtering method which combines wavelet denoising and median filtering is used to denoise the signal. Reference [3] proposed a strong tracking self-feedback model based on recursive least square multi-wavelet decomposition and reconstruction, established a new soft threshold function, and combined it with improved median filter to denoise. According to the characteristics of hard threshold and soft threshold, a custom threshold function is proposed for wavelet threshold denoising in reference [4]. It can be seen that the difficulty of the wavelet filtering method lies in the selection of the wavelet basis function and the determination of the decomposition layer number and the threshold function, which lacks algorithm adaptability. In reference [5], an error model of deep simple recurrent unit recurrent neural network (SRU-RNN) is proposed to reduce the parameters that need to be determined during the training process, but the fixed learning speed and training data lead to the reduction of model accuracy and are prone to overfitting problem. In reference [6], a new method of segmented compensation for the temperature error of fiber-optic gyroscope is proposed to model the data in different temperature intervals, which can avoid overfitting problems and improve model accuracy. erefore, the neural network modelling method theoretically has the ability to approximate nonlinear functions with arbitrary precision and has highspeed parallel computing capability, but it has a complex network computing structure and is prone to overfitting problems. e method of timing analysis and establishing reasonable ARMA model for random drift error is the most widely used method, the model has high accuracy, and this method has achieved good results in gyro denoising. In references [7][8][9][10], the most suitable and reasonable model of random drift error is established after the data pre-processing and a series of testing, and then various filtering algorithms are carried out to achieve compensation for error. However, the premise of this method is that the sequence to be processed should be a stationary sequence, and a complicated pre-processing process is required for the non-stationary sequence. In addition, the empirical mode decomposition (EMD) is a new adaptive method for processing nonlinear and non-stationary signals. is method does not require any prior knowledge of the signal. It decomposes the original data into a finite number of intrinsic mode functions (IMFs) and a margin, which is an effective method for data stabilization and denoising [11][12][13][14]. Traditional EMD thinks that noise mainly exists in the high frequency component and directly removes it, but there is no certain screening criterion to denoise, especially for the signal with low signal-to-noise ratio; noise and useful signal are often mixed together. Although the method of directly removing the high frequency component can achieve good denoising effect, some useful signals will be lost at the same time.
Based on the above analysis, this paper adopts a new method named IEMD to deal with the random error of MEMS gyroscope. Firstly, the method decomposes the original signal into several IMF components by traditional EMD. And then, calculate the Bhattacharyya distance between the signal and every component, and the maximum slope between the adjacent two modes j, j + 1 and the input signal is the boundary point. ereby, the first demarcation point j is determined; it means IMF i ∼ IMF j are noise IMFs; secondly, the autocorrelation function of each component is analyzed, and the second demarcation point k is found according to the properties of autocorrelation function. It is determined that IMF j+1 ∼ IMF k are the mixed IMFs, and the remaining IMF k+1 ∼ IMF n are signal IMFs. Finally, the noise dominated IMF 1 ∼ IMF j are removed directly, the signal dominated IMF k+1 ∼ IMF n are retained directly, and the mixed IMFs are modelled and filtered; then, the processed mixed IMFs and signal IMFs and the remainder are reconstructed. e final gyro signal is obtained. Figure 1 shows the design block diagram of error denoising method.

Global Design Module Diagram of Denoising Method
is method is mainly the following three steps:

Interval Empirical Mode
Decomposition. e gyro output data can be expressed as follows after traditional EMD: (1)

Mathematical Problems in Engineering
In equation (1), IMF i is the component of intrinsic mode function, and r(t) is the remainder. e IMF 1 ∼ IMF i components are distributed from high frequency to low frequency in turn. e gyro output data can be expressed as follows after IEMD: (2) In equation (2), IMF 1 − IMF j are the noise IMFs, IMF j+1 − IMF k are the mixed IMFs, IMF k+1 − IMF n are the signal IMFs, and r(t) is the remainder. e noise IMFs can be directly removed, and the signal IMFs are directly retained, mainly processing the mixed IMFs, modelling and filtering them.
In equation (3), IMF i ′ are the mixed IMFs after filtering and modelling.

Screening Noise IMFs Based on Bhattacharyya
Distance. e probability density function (PDF) can reflect the difference between the signals [15]. Because of this, the kernel density estimation method is used to obtain the probability density function of the input signal and each component, and the different components are distinguished by calculating the similarity between them. e Bhattacharyya distance can be used to measure the distance between two PDFs and it is an effective method to the judge the similarity [16].
In the same definition domain X, the Bhattacharyya distance of the probability distributions P and Q is defined as follows: Among equation (4) for discrete probability distribution, For continuous probability distribution, e smaller D B is, the closer the probability distribution is, and the more relevant the modal component is to the input signal.
In this paper, the similarity between the modal component and the input signal is defined as follows: e distinction between correlated and uncorrelated modes can be determined by evaluating the slope of the distance between two adjacent modes and the input signal; the max slope is defined as follows: e boundary between the uncorrelated mode and the correlated mode is j; that is, IMF 1 ∼ IMF j are noise IMFs.

Screening Signal IMFs Based on Autocorrelation
Function. e second demarcation point, that is, the boundary of the mixed IMFs and the signal IMFs, is determined by the characteristics of the autocorrelation function; the signal IMFs can be screened out according to autocorrelation function characteristics.
Autocorrelation function: In equation (9), x(t) is a random signal. e characteristics of the autocorrelation function: for a data sequence containing random noise, the function value of the autocorrelation function is the largest at zero, and the function value of the other points rapidly declines to zero, showing a weak correlation. e data sequence is dominated by useful signals; although the autocorrelation function value is also the largest at zero, the function value of other points slowly declines to zero; there is a certain regular change, showing a strong correlation. erefore, based on this characteristic, the demarcation point of mixed IMFs and signal IMFs can be determined.
After the IEMD, the original signal can be expressed as In equation (10), IMF i − IMF j are the noise IMFs, IMF j+1 − IMF k are the mixed IMFs, IMF k+1 − IMF n are the signal IMFs, and r(t) is the remainder.
In the next processing, the noise IMFs can be directly removed, and the signal IMFs are directly retained, mainly processing the mixed IMFs, modelling and filtering them.

Time Series Modelling.
Reasonable modelling of the mixed IMFs obtained by the above screening is required. Since the IMFs obtained by EMD decomposition are a stationary signal, it is only necessary to perform the stationarity and normality test on the signal during preprocessing.

Model Identification.
e time series model mainly includes three types, namely, autoregressive model, moving average model, and autoregressive moving average model. Calculate the autocorrelation function and partial correlation function of the sequence; then, identify the model type according to the different statistical characteristics of the model shown in Table 1 [17], and select the suitable model type for each sequence.

Model Ordering.
After selecting the suitable model, the BIC criteria are used to determine the order of the model:

Mathematical Problems in Engineering
In equation (11), σ 2 a is variance of the residual and N is data length of the residual sequence.
Take the model order p with the minimum BIC(p) value as the applicable model order.

Model Parameter Estimation.
e selected model parameter estimation algorithm is the least squares estimation method. e least squares estimation method is actually an unbiased estimation; the basic idea is to obtain random vector ϕ and calculate the auto-covariance function of ϕ, named C(φ); the value of the C(φ) ′ 's main diagonal element is used as a parameter estimate.

Model Applicability Test.
After the model is completed, the applicability of the model is tested, and the white noise test criterion is adopted. is criterion is to use different statistics to check whether the model residual is white noise, and the residual sequence is closer to white noise; the model accuracy is better. e mathematical expression of the model is en, the expression of the residual sequence a t is In equations (12) and (13), p and q are the order of the AR part and the MA part, respectively; φ 1 , . . . , φ p and θ 1 , . . . , θ q are the estimated parameters of the AR part and the MA part, respectively; and a t is the residual sequence. e model is tested for applicability based on the autocorrelation coefficient criteria, after calculating the residual sequence a t . e test formula for the autocorrelation coefficient criterion [18] is Among them, If the autocorrelation coefficient ρ a,k satisfies this formula, a t is white noise, and the corresponding model is applicable.

Kalman Filtering Based on the AR Model.
Kalman filtering is performed on the above established model. Taking the AR (4) model as an example, the model expression is en, the model equation is transformed into a state space model; the obtained state equation and measurement equation are as follows: In equations (17) and (18), x k is the system state, y k is the system measurement, φ k/k− 1 is the state transition matrix , ω k− 1 is the system noise, v k is the measurement noise, and ω k− 1 , v k are usually Gaussian white noise sequence with zero mean.
According to the five recursive formulas and the parameters of Kalman filtering, the optimal estimation of the signal can be obtained [19]. e optimal estimation value of each component is obtained using the same method, and finally the signal is reconstructed. e reconstruction formula of the signal is as follows: In equation (19), IMF i ′ are the mixed IMFs filtered and modelled. Figure 2, the device used in this experiment is an inertial navigation rotating platform. e performance parameters of the inertial measurement unit 3DM-10A are shown in Table 2. Take the X-axis data outputted from the MEMS gyroscope as the experimental object. Firstly, the device is turned on and preheated for 5 minutes to stabilize the working state, and then the signal output from the gyro is sampled and collected. e sampling frequency is 100 Hz, the sampling time is about 1 h, and the static X-axis data from the gyroscope is shown in Figure 3.

Static Data Denoising Experiment. As shown in
en, decompose and screen the above original signal. As shown in Figure 4, after decomposition, 15 IMF components and one remainder R are obtained.  (7). As clearly shown in Figure 5, the slope between BLIMF6 and BLINF7 is the largest. erefore, IMF1-IMF6 are the noise IMFs; that is, the first boundary point is j � 6.

Analysis of Autocorrelation Function.
In order to determine the second demarcation point k, that is, the boundary between the mixed IMFs and the signal IMFs, the autocorrelation function of each order IMF component is calculated and normalized, as shown in Figure 6. It can be concluded from the characteristics of the autocorrelation function that the autocorrelation function value of IMF 1 − IMF 10 is the largest at the zero point, and the other points are rapidly attenuated to zero, which can be initially determined as it contains noise, while the function value of IMF 11 − IMF 15 has a certain regularity change, which can be determined as it is dominated by the signal. In order to further accurately determine the second demarcation point k, the variance threshold method is used to verify it. As shown in Figure 7, the variance of the IMF autocorrelation function of each order is calculated. According to the variance threshold method, the variances of the autocorrelation function of IMF 1 − IMF 10 are all less than 0.01; nevertheless, the variances of the autocorrelation function of IMF 11 − IMF 15 are significantly larger than 0.01 and increase exponentially. erefore, IMF 6 − IMF 10 are mixed IMFs that can be accurately determined; that is, the next step modelling and filtering are required. A reasonable model is established for the filtered component signals (because of limited space, this article only takes IMF 6 as an example). Firstly, the IMF 6 component signal is tested for zero-mean and normal distribution, as shown in Figures 8  and 9. After testing, it is concluded that the IMF 6 component satisfies the modelling conditions. e same method tests the IMF 8 − IMF 10 in turn, all of which meet the modelling conditions.
Next, model identification and model ordering of component signals are performed, and the autocorrelation function curve and the partial correlation function curve of IMF 6 are drawn, respectively, as shown in Figure 10. It can be seen from the figure that the autocorrelation function is tailed and the partial autocorrelation function presents a p-step censored. According to Table 1, the model is an AR model. At the same time, according to the BIC criterion, as shown in Figure 11, the number of BIC functions with the minimum value is selected as the optimal order, and the optimal order of the model of IMF 6 is 4, so the model is judged as the AR (4) model. After the model is determined, the parameters of the model are estimated by least squares estimation method, and finally the model applicability test is performed.
In the same way, IMF 6 − IMF 10 are modelled, respectively. After a series of modelling steps and applicability tests, the final modelling results are shown in Table 3.
Finally, Kalman filtering is performed on the above model, and the parameters are set and updated. After each component filter is updated, the signals are reconstructed according to equation (3). e reconstruction result is the final denoising signal. ree different schemes are adopted to denoise the same set of data, respectively: (i) Scheme 1: direct modelling method (ii) Scheme 2: traditional EMD method (iii) Scheme 3: the method in this article e denoising results of each scheme are shown in Figure 12.
Additionally, the root mean square error (RMSE) and signal-to-noise ratio (SNR) are introduced as evaluation criteria [20]. e definitions of the two evaluation indicators are as follows: In equations (20) and (21), x i is the original signal; x i ′ is the denoising signal; and N is the length of the signal. e calculation results are shown in Table 4. It can be seen from Table 5 that the RMSE of the original signal is 0.3497°/s. After three different denoising schemes, the RMSE of the original signal is reduced to different degrees, and the signal-to-noise ratio is improved to  different degrees. Scheme 3 proposed in this paper is the most effective, the RMSE is minimum, and the signal-tonoise ratio is the largest. Compared with the original signal, scheme 1, and scheme 2, the RMSE of scheme 3 decreased by 71%, 48%, and 36%, respectively, and the signal-to-noise ratio increased by 48%, 28%, and 6%. It is proven that scheme 3 proposed in this paper has effective denoising effect and has obvious advantages compared with scheme 1 and scheme 2. However, the improvement of the signal-tonoise ratio of scheme 3 is not significantly improved compared with scheme 2, but as can be seen from Figure 9, the integrity and smoothness of the signal processed by scheme 3 are better, the signal is smoother, and many peaks are reduced. In order to further verify the effect of scheme 3 proposed in this paper, Allan's variance comparison analysis experiment is calculated.

Allan Variance Comparison Analysis.
Allan variance analysis is an effective method for gyro random error identification and noise characteristics analysis [21]. e method can identify multiple different types of random errors in different time domains and can classify the error term into five error terms, quantization noise, angular random walk, zero offset instability, angular rate random walk, and speed ramp, and the error coefficients can be analyzed quantitatively [22,23]. In order to further verify the validity and applicability of the proposed method, the Allan variance double logarithm graph of the signals obtained by the three schemes is plotted. As shown in Figure 13, it can be seen that the Allan variance of original signal is the largest, and the Allan variances of scheme 1, scheme 2, and scheme 3 are sequentially reduced.
e Allan variance of scheme 3 proposed in this paper is the smallest and has been reduced to 10 1 magnitude.
In addition, the error term coefficients of the original signal and the signals processed by the three schemes are, respectively, obtained by the fitting method and recorded in Table 5. e parameters Q, N, B, K, and R in the table are, respectively, quantified noise coefficient, angular random walk coefficient, zero offset instability coefficient, angular rate random walk coefficient, and speed ramp coefficient (for    Mathematical Problems in Engineering convenience of operation, the unit°/s is converted to°/h). It can be obtained from the table that the coefficients of each error term are reduced to different degrees through three kinds of denoising schemes, and the coefficients of each error term of the signal are the smallest after scheme 3 processing; the error coefficients Q, N, B, K, and R of the signal processed by scheme 3 are reduced by 70%, 83%, 43%, 35%, and 34%, respectively, which fully demonstrates that the method proposed in this paper has the best denoising effect.

Dynamic Data Denoising Experiment.
e experiment of dynamic data is closer to the actual application scene. erefore, a dynamic data denoising experiment is set up. e experiment of dynamic data is divided into two parts. e first part is the experiment of MEMS gyroscope rotating at a constant speed. Both the main axis and the pitch axis of the turntable rotate at a constant speed of 10°/s, where the main axis rotates in the clockwise direction and the pitch axis rotates in the counter clockwise direction. e experimental results take the X-axis output signal as an example. In the experiment, three different schemes proposed in Section 3.1 are still used to compare the denoising results. e results are shown in Figure 14 and Table 6. It can be seen from Table 7 that the denoising result of scheme 3 proposed in this paper is significantly better than scheme 1 and scheme 2. e standard deviation is reduced from 0.0198 to 0.0083, and the signal-to-noise ratio is increased from 10.9428 dB to 32.4631 dB. e second part of the experiment is the MEMS gyroscope swinging motion experiment. In this part, the main axis of the turntable and the pitch axis are both oscillated at a rate of 10°/s. e experimental results take the Z-axis output signal as an example. e results are shown in Figure 15 and Table 7.
It can be seen from Table 8 that the denoising result of scheme 3 proposed in this paper is significantly better than scheme 1 and scheme 2. e standard deviation is reduced     Mathematical Problems in Engineering from 0.4876 to 0.0799, and the signal-to-noise ratio is increased from 13.2743 dB to 25.4896 dB.

Vehicle Experiment.
e ultimate purpose of gyro denoising is to improve the navigation accuracy of inertial navigation system, so the position information of the inertial navigation solution can directly reflect the advantages and disadvantages of the algorithm [24]. erefore, the following vehicle experiment is designed. e vehicle experimental device is shown in Figure 16. e experimental device is GPS/IMU combined positioning system, and the inertial unit is IMU-200A. e performance parameters of vehicle sensors are shown in Table 8. e vehicle experimental environment is relatively complicated and there are many interference factors, especially in the satellite occlusion area; the inertial navigation solution data will have serious error divergence. erefore, in this part of the experiment, select the data within 20 s of the satellite occlusion area under the vehicle linear motion state, and set different solutions to solve the problem: (i) Scheme 1: calculate the original data of the gyro outputted by the inertial navigation system (ii) Scheme 2: calculate the gyro data processed by traditional EMD (iii) Scheme 3: calculate the gyro data processed by the method proposed in this paper e longitude and latitude error results of the three schemes are shown in Figure 17, and the corresponding mean, standard deviation, and maximum error are compared as shown in Table 9. Scheme 3 proposed in this paper optimizes the position information to the maximum extent. e mean of the longitude error and the latitude error is reduced to 1.3793 m and 2.0689 m, respectively; compared with scheme 2, it decreased by 48.5% and 25%; the standard deviation is reduced to 1.1120 m and 1.6800 m; compared with scheme 2, it decreased by 65.5% and 25%. e maximum error is reduced to 3.8464 m and 5.7696 m, respectively.
Based on the above analysis, in the complex vehicle test environment, the position information of the inertial navigation calculation is significantly optimized. Compared with the traditional EMD method, the error divergence of inertial navigation system is further suppressed, and the inertial navigation calculation accuracy is improved.

Conclusions
Based on the traditional EMD algorithm, this paper proposes a screening mechanism, which combines Bhattacharyya distance analysis and the characteristics of autocorrelation function to classify the IMF into three categories, namely, noise IMFs, mixed IMFs, and signal IMFs. And then the mixed IMFs are modelled and filtered; finally, the signal is reconstructed. Static experiment and dynamic experiment on the turntable are set to verify the  performance of the algorithm.
e experimental results show that the denoising effect of the proposed method is better than that of the traditional EMD denoising method. In addition, the vehicle experiment in complex environment is designed; the results show that the accuracy of inertial navigation position calculation is still significantly optimized. It fully proves that the method has obvious denoising effect, greatly improves the signal quality, improves the accuracy of the inertial navigation calculation, and has certain guiding significance for engineering applications.

Data Availability
e data used to support the findings of this study have not been made available because the data also form part of an ongoing study at this time.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.