Abstract
The microelectromechanical system (MEMS) gyroscope has low measurement accuracy and large output noise; the useful signal is often submerged in the noise. A new denoising method 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. Then, 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.
1. 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. Therefore, 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. Therefore, the neural network modelling method theoretically has the ability to approximate nonlinear functions with arbitrary precision and has high-speed parallel computing capability, but it has a complex network computing structure and is prone to overfitting problems. The 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–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. This 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–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. Thereby, the first demarcation point j is determined; it means 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 are the mixed IMFs, and the remaining are signal IMFs. Finally, the noise dominated are removed directly, the signal dominated are retained directly, and the mixed IMFs are modelled and filtered; then, the processed mixed IMFs and signal IMFs and the remainder are reconstructed. The final gyro signal is obtained.
2. Global Design Module Diagram of Denoising Method
Figure 1 shows the design block diagram of error denoising method.
This method is mainly the following three steps:(i)Step 1. Data decomposition and screening. Firstly, the original signal is decomposed into multiple IMFs and a residual by EMD algorithm; then, Bhattacharyya distance analysis and autocorrelation function analysis of each component are performed to divide the signal into noise IMFs, mixed IMFs, and signal IMFs.(ii)Step 2. Time series modelling of the mixed IMFs. After the stationarity testing and normality testing, a series of modelling steps are performed on the filtered mixed IMFs to establish a reasonable time series model for the components.(iii)Step 3. Filtering and reconstruction. Kalman filtering is performed on the basis of the model to remove the random error. Finally, the filtered mixed IMFs and the signal IMFs and the remainder are reconstructed, and the final result is outputted.
3. Denoising Method Based on Interval Empirical Mode Decomposition
3.1. Interval Empirical Mode Decomposition
The gyro output data can be expressed as follows after traditional EMD:
In equation (1), is the component of intrinsic mode function, and is the remainder. The components are distributed from high frequency to low frequency in turn.
The gyro output data can be expressed as follows after IEMD:
In equation (2), are the noise IMFs, are the mixed IMFs, are the signal IMFs, and is the remainder.
The 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), are the mixed IMFs after filtering and modelling.
3.1.1. Screening Noise IMFs Based on Bhattacharyya Distance
The 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. The 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,
The smaller 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:
The 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:
The boundary between the uncorrelated mode and the correlated mode is j; that is, are noise IMFs.
3.1.2. Screening Signal IMFs Based on Autocorrelation Function
The 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), is a random signal.
The 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. The 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. Therefore, 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), are the noise IMFs, are the mixed IMFs, are the signal IMFs, and 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.
3.2. 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 pre-processing.
3.2.1. Model Identification
The 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.
3.2.2. Model Ordering
After selecting the suitable model, the BIC criteria are used to determine the order of the model:
In equation (11), is variance of the residual and N is data length of the residual sequence.
Take the model order p with the minimum value as the applicable model order.
3.2.3. Model Parameter Estimation
The selected model parameter estimation algorithm is the least squares estimation method. The 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 ; the value of the ’s main diagonal element is used as a parameter estimate.
3.2.4. Model Applicability Test
After the model is completed, the applicability of the model is tested, and the white noise test criterion is adopted. This 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.
The mathematical expression of the model is
Then, the expression of the residual sequence is
In equations (12) and (13), p and q are the order of the AR part and the MA part, respectively; and are the estimated parameters of the AR part and the MA part, respectively; and is the residual sequence.
The model is tested for applicability based on the autocorrelation coefficient criteria, after calculating the residual sequence .
The test formula for the autocorrelation coefficient criterion [18] is
Among them,
If the autocorrelation coefficient satisfies this formula, is white noise, and the corresponding model is applicable.
3.3. 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
Then, 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), is the system state, is the system measurement, is the state transition matrix , , , is the measurement matrix , is the system noise, is the measurement noise, and , 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]. The optimal estimation value of each component is obtained using the same method, and finally the signal is reconstructed. The reconstruction formula of the signal is as follows:
In equation (19), are the mixed IMFs filtered and modelled.
4. Experimental Verification and Analysis
4.1. Static Data Denoising Experiment
As shown in Figure 2, the device used in this experiment is an inertial navigation rotating platform. The 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. The 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.
Then, decompose and screen the above original signal. As shown in Figure 4, after decomposition, 15 IMF components and one remainder R are obtained.
4.1.1. Analysis of Bhattacharyya Distance
Use the kernel density estimation method to obtain the PDF of the original signal and each component, and then, calculate the Bhattacharyya distance between them according to formula (7). As clearly shown in Figure 5, the slope between BLIMF6 and BLINF7 is the largest. Therefore, IMF1-IMF6 are the noise IMFs; that is, the first boundary point is .
4.1.2. Analysis of Autocorrelation Function
In order to determine the second demarcation point , 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 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 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 , 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 are all less than 0.01; nevertheless, the variances of the autocorrelation function of are significantly larger than 0.01 and increase exponentially. Therefore, 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 as an example). Firstly, the component signal is tested for zero-mean and normal distribution, as shown in Figures 8 and 9. After testing, it is concluded that the component satisfies the modelling conditions. The same method tests the in turn, all of which meet the modelling conditions.
(a)
(b)
Next, model identification and model ordering of component signals are performed, and the autocorrelation function curve and the partial correlation function curve of 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 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.
(a)
(b)
In the same way, 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). The reconstruction result is the final denoising signal. Three 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
The 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]. The definitions of the two evaluation indicators are as follows:
In equations (20) and (21), is the original signal; is the denoising signal; and N is the length of the signal.
The 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-to-noise 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-to-noise 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.
4.2. Allan Variance Comparison Analysis
Allan variance analysis is an effective method for gyro random error identification and noise characteristics analysis [21]. The 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. The Allan variance of scheme 3 proposed in this paper is the smallest and has been reduced to 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. The 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 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.
4.3. Dynamic Data Denoising Experiment
The experiment of dynamic data is closer to the actual application scene. Therefore, a dynamic data denoising experiment is set up. The experiment of dynamic data is divided into two parts. The 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. The 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. The 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. The 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.
The 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. The experimental results take the Z-axis output signal as an example. The 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. The standard deviation is reduced from 0.4876 to 0.0799, and the signal-to-noise ratio is increased from 13.2743 dB to 25.4896 dB.
4.4. Vehicle Experiment
The 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]. Therefore, the following vehicle experiment is designed.
The vehicle experimental device is shown in Figure 16. The experimental device is GPS/IMU combined positioning system, and the inertial unit is IMU-200A. The performance parameters of vehicle sensors are shown in Table 8.
The 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. Therefore, 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
The 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. The 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%. The maximum error is reduced to 3.8464 m and 5.7696 m, respectively.
(a)
(b)
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.
5. 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. The 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
The 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
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was supported in part by the National Natural Science Foundation of China (61863024 and 71761023) and the Gansu Provincial Science Guidance Plan of China (2020-61).