# Denoising Method of MEMS Gyroscope Based on Interval Empirical Mode Decomposition

Authors: Yang Liu; Guangwu Chen; Zongshou Wei; Juhua Yang; Dongfeng Xing
Journal: Mathematical Problems in Engineering
 (2020)

---

## 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.


---

## Body

## 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 modesj, j + 1 and the input signal is the boundary point. Thereby, the first demarcation point j is determined; it means IMFi∼IMFj 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 IMFj+1∼IMFk are the mixed IMFs, and the remaining IMFk+1∼IMFn are signal IMFs. Finally, the noise dominated IMF1∼IMFj are removed directly, the signal dominated IMFk+1∼IMFn 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

Figure1 shows the design block diagram of error denoising method.Figure 1
Global design module diagram of denoising method for MEMS gyroscope.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 remainderrt 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:(1)yt=∑i=1nIMFi+rt.In equation (1), IMFi is the component of intrinsic mode function, and rt is the remainder. The IMF1∼IMFi components are distributed from high frequency to low frequency in turn.The gyro output data can be expressed as follows after IEMD:(2)yt=∑i=1jIMFi+∑j+1kIMFi+∑k+1nIMFi+rt.In equation (2), IMF1−IMFj are the noise IMFs, IMFj+1−IMFk are the mixed IMFs, IMFk+1−IMFn are the signal IMFs, and rt 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.(3)yt′=∑j+1kIMFi′+∑k+1nIMFi+rt.In equation (3), IMFi′ 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 domainX, the Bhattacharyya distance of the probability distributions P and Q is defined as follows:(4)DBP,Q=−lnBCP,Q.Among equation (4) for discrete probability distribution,(5)BCP,Q=∑x∈XPxQx.For continuous probability distribution,(6)BCP,Q=∫PxQxdx.The smallerDB 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:(7)Li=distancepdfxt,pdfBLIMFit.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:(8)θ=maxLj+1−Lj,j=1,2,…,n−1.The boundary between the uncorrelated mode and the correlated mode isj; that is, IMF1∼IMFj 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:(9)Rxt1,t2=Ext1,xt2.In equation (9), xt 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(10)yt=∑i=1jIMFi+∑j+1kIMFi+∑k+1nIMFi+rt.In equation (10), IMFi−IMFj are the noise IMFs, IMFj+1−IMFk are the mixed IMFs, IMFk+1−IMFn are the signal IMFs, and rt 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 Table1 [17], and select the suitable model type for each sequence.Table 1
The characteristics of time series model.
Model nameARpMAqARMAp,qAutocorrelation coefficient functionTrailingq step censoringTrailingPartial autocorrelation coefficient functionP step censoringTrailingTrailing


#### 3.2.2. Model Ordering

After selecting the suitable model, the BIC criteria are used to determine the order of the model:(11)BICp=Nlnσa2+plnN.In equation (11), σa2 is variance of the residual and N is data length of the residual sequence.Take the model orderp with the minimum BICp 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 Cφ; the value of the Cφ′’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(12)xt=φ1xt−1+⋯+φpxt−p+at−θ1at−1−⋯−θqat−q.Then, the expression of the residual sequenceat is(13)at=xt−φ1xt−1−⋯−φpxt−p+θ1at−1+⋯+θqat−q.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 at is the residual sequence.The model is tested for applicability based on the autocorrelation coefficient criteria, after calculating the residual sequenceat.The test formula for the autocorrelation coefficient criterion [18] is(14)Nρ^a,k1+2∑i=1k−1ρ^a,k2≤1.96.Among them,(15)ρ^a,k=∑i=k+1Nxtxt−k∑i=k+1Nxt2,k=1,2,…,.If the autocorrelation coefficientρ^a,k satisfies this formula, at 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(16)xt=φ1xt−1+φ2xt−2+φ3xt−3+φ4xt−4+at.Then, the model equation is transformed into a state space model; the obtained state equation and measurement equation are as follows:(17)xk=φk/k−1xk−1+ωk−1,(18)yk=Hkxk+vk.In equations (17) and (18), xk is the system state, yk is the system measurement, φk/k−1 is the state transition matrix φk/k−1=φA, ϕ=φ1φ2,…,φp, A=1p−1×p−10p−1×1, Hk is the measurement matrix Hk=101×p−1, ωk−1 is the system noise, vk is the measurement noise, and ωk−1, vk 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:(19)yt′=∑j+1kIMFi′+∑k+1nIMFi+rt.In equation (19), IMFi′ are the mixed IMFs filtered and modelled.


## 3.1. Interval Empirical Mode Decomposition

The gyro output data can be expressed as follows after traditional EMD:(1)yt=∑i=1nIMFi+rt.In equation (1), IMFi is the component of intrinsic mode function, and rt is the remainder. The IMF1∼IMFi components are distributed from high frequency to low frequency in turn.The gyro output data can be expressed as follows after IEMD:(2)yt=∑i=1jIMFi+∑j+1kIMFi+∑k+1nIMFi+rt.In equation (2), IMF1−IMFj are the noise IMFs, IMFj+1−IMFk are the mixed IMFs, IMFk+1−IMFn are the signal IMFs, and rt 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.(3)yt′=∑j+1kIMFi′+∑k+1nIMFi+rt.In equation (3), IMFi′ 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 domainX, the Bhattacharyya distance of the probability distributions P and Q is defined as follows:(4)DBP,Q=−lnBCP,Q.Among equation (4) for discrete probability distribution,(5)BCP,Q=∑x∈XPxQx.For continuous probability distribution,(6)BCP,Q=∫PxQxdx.The smallerDB 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:(7)Li=distancepdfxt,pdfBLIMFit.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:(8)θ=maxLj+1−Lj,j=1,2,…,n−1.The boundary between the uncorrelated mode and the correlated mode isj; that is, IMF1∼IMFj 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:(9)Rxt1,t2=Ext1,xt2.In equation (9), xt 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(10)yt=∑i=1jIMFi+∑j+1kIMFi+∑k+1nIMFi+rt.In equation (10), IMFi−IMFj are the noise IMFs, IMFj+1−IMFk are the mixed IMFs, IMFk+1−IMFn are the signal IMFs, and rt 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.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 domainX, the Bhattacharyya distance of the probability distributions P and Q is defined as follows:(4)DBP,Q=−lnBCP,Q.Among equation (4) for discrete probability distribution,(5)BCP,Q=∑x∈XPxQx.For continuous probability distribution,(6)BCP,Q=∫PxQxdx.The smallerDB 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:(7)Li=distancepdfxt,pdfBLIMFit.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:(8)θ=maxLj+1−Lj,j=1,2,…,n−1.The boundary between the uncorrelated mode and the correlated mode isj; that is, IMF1∼IMFj 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:(9)Rxt1,t2=Ext1,xt2.In equation (9), xt 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(10)yt=∑i=1jIMFi+∑j+1kIMFi+∑k+1nIMFi+rt.In equation (10), IMFi−IMFj are the noise IMFs, IMFj+1−IMFk are the mixed IMFs, IMFk+1−IMFn are the signal IMFs, and rt 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 Table1 [17], and select the suitable model type for each sequence.Table 1
The characteristics of time series model.
Model nameARpMAqARMAp,qAutocorrelation coefficient functionTrailingq step censoringTrailingPartial autocorrelation coefficient functionP step censoringTrailingTrailing


### 3.2.2. Model Ordering

After selecting the suitable model, the BIC criteria are used to determine the order of the model:(11)BICp=Nlnσa2+plnN.In equation (11), σa2 is variance of the residual and N is data length of the residual sequence.Take the model orderp with the minimum BICp 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 Cφ; the value of the Cφ′’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(12)xt=φ1xt−1+⋯+φpxt−p+at−θ1at−1−⋯−θqat−q.Then, the expression of the residual sequenceat is(13)at=xt−φ1xt−1−⋯−φpxt−p+θ1at−1+⋯+θqat−q.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 at is the residual sequence.The model is tested for applicability based on the autocorrelation coefficient criteria, after calculating the residual sequenceat.The test formula for the autocorrelation coefficient criterion [18] is(14)Nρ^a,k1+2∑i=1k−1ρ^a,k2≤1.96.Among them,(15)ρ^a,k=∑i=k+1Nxtxt−k∑i=k+1Nxt2,k=1,2,…,.If the autocorrelation coefficientρ^a,k satisfies this formula, at is white noise, and the corresponding model is applicable.


## 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 Table1 [17], and select the suitable model type for each sequence.Table 1
The characteristics of time series model.
Model nameARpMAqARMAp,qAutocorrelation coefficient functionTrailingq step censoringTrailingPartial autocorrelation coefficient functionP step censoringTrailingTrailing


## 3.2.2. Model Ordering

After selecting the suitable model, the BIC criteria are used to determine the order of the model:(11)BICp=Nlnσa2+plnN.In equation (11), σa2 is variance of the residual and N is data length of the residual sequence.Take the model orderp with the minimum BICp 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 Cφ; the value of the Cφ′’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(12)xt=φ1xt−1+⋯+φpxt−p+at−θ1at−1−⋯−θqat−q.Then, the expression of the residual sequenceat is(13)at=xt−φ1xt−1−⋯−φpxt−p+θ1at−1+⋯+θqat−q.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 at is the residual sequence.The model is tested for applicability based on the autocorrelation coefficient criteria, after calculating the residual sequenceat.The test formula for the autocorrelation coefficient criterion [18] is(14)Nρ^a,k1+2∑i=1k−1ρ^a,k2≤1.96.Among them,(15)ρ^a,k=∑i=k+1Nxtxt−k∑i=k+1Nxt2,k=1,2,…,.If the autocorrelation coefficientρ^a,k satisfies this formula, at 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(16)xt=φ1xt−1+φ2xt−2+φ3xt−3+φ4xt−4+at.Then, the model equation is transformed into a state space model; the obtained state equation and measurement equation are as follows:(17)xk=φk/k−1xk−1+ωk−1,(18)yk=Hkxk+vk.In equations (17) and (18), xk is the system state, yk is the system measurement, φk/k−1 is the state transition matrix φk/k−1=φA, ϕ=φ1φ2,…,φp, A=1p−1×p−10p−1×1, Hk is the measurement matrix Hk=101×p−1, ωk−1 is the system noise, vk is the measurement noise, and ωk−1, vk 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:(19)yt′=∑j+1kIMFi′+∑k+1nIMFi+rt.In equation (19), IMFi′ are the mixed IMFs filtered and modelled.


## 4. Experimental Verification and Analysis

### 4.1. Static Data Denoising Experiment

As shown in Figure2, 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.Figure 2
Inertial platform test system.Table 2
Performance parameters of 3DM-10A.
ParametersValueUnitGyro bias stability0.03°/sGyro random walk0.03°/(s1/2)Accelerometer bias stability0.03°/sAccelerometer random walk0.03°/(s1/2)Resolution±0.01°/sFigure 3
Static data from gyroX-axis.Then, decompose and screen the above original signal. As shown in Figure4, after decomposition, 15 IMF components and one remainder R are obtained.Figure 4
The IMF components and remainder of the original signal.


#### 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 j=6.Figure 5
Bhattacharyya distance between input signals and each component.


#### 4.1.2. Analysis of Autocorrelation Function

In order to determine the second demarcation pointk, 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 IMF1−IMF10 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 IMF11−IMF15 has a certain regularity change, which can be determined as it is dominated by the signal.Figure 6
Autocorrelation functions of each IMF.In order to further accurately determine the second demarcation pointk, 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 IMF1−IMF10 are all less than 0.01; nevertheless, the variances of the autocorrelation function of IMF11−IMF15 are significantly larger than 0.01 and increase exponentially. Therefore, IMF6−IMF10 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 IMF6 as an example). Firstly, the IMF6 component signal is tested for zero-mean and normal distribution, as shown in Figures 8 and 9. After testing, it is concluded that the IMF6 component satisfies the modelling conditions. The same method tests the IMF8−IMF10 in turn, all of which meet the modelling conditions.Figure 7
Variances of ACF of each IMF.Figure 8
Normality test of IMF6.Figure 9
Mean and mean square error changing chart. (a) Mean change. (b) Variance of mean square error.
(a)(b)Next, model identification and model ordering of component signals are performed, and the autocorrelation function curve and the partial correlation function curve ofIMF6 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 IMF6 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.Figure 10
Model identification of IMF6. (a) Autocorrelation coefficient graph. (b) Partial autocorrelation coefficient graph.
(a)(b)Figure 11
BIC criterion for ordering determination.In the same way,IMF6−IMF10 are modelled, respectively. After a series of modelling steps and applicability tests, the final modelling results are shown in Table 3.Table 3
Modelling results ofIMF6−IMF10.
Component signalModelφ1φ2φ3φ4φ5φ6IMF6AR (4)−3.47994.4906−2.53360.5232IMF7AR (6)−2.42961.5645−0.0048−0.0099−0.0673−0.0528IMF8AR (4)−2.43961.44870.4232−0.4322IMF9AR (4)−2.09451.09460.09500.0950IMF10AR (4)−2.00821.00820.0088−0.0088Finally, 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 articleThe denoising results of each scheme are shown in Figure12.Figure 12
Denoising results of each scheme.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:(20)RMSE=1n∑i=1Nxi−xi′2,(21)SNR=10log∑i=1Nxi2∑i=1Nxi−xi′2.In equations (20) and (21), xi is the original signal; xi′ is the denoising signal; and N is the length of the signal.The calculation results are shown in Table4.Table 4
Denoising results comparison of each scheme.
Data/schemeRMSE (°/s)SNR (dB)Original signal0.349711.9505Scheme 10.193016.4550Scheme 20.155321.6059Scheme 30.100122.9756It can be seen from Table5 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.Table 5
Comparison of signal error term coefficients of each scheme.
Data/SchemeQ (°)N (°/h1/2)B (°/h)K (°/h3/2)R (°/h2)Original signal0.00270.010918.895026.01616.0252Scheme 10.00230.008714.761421.35085.0599Scheme 20.00200.005614.131721.97755.2616Scheme 30.00080.001810.714116.79333.9505


### 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 Figure13, 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 101 magnitude.Figure 13
Allan variances comparison of each scheme.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 Table5. 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 theX-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.Figure 14
Rotating data denoising results comparison of each scheme.Table 6
Denoising results comparison of rotating data.
Data/schemeRMSE (°/s)SNR (dB)Original signal0.019810.9428Scheme 10.019316.2687Scheme 20.019120.9575Scheme 30.008332.4631It can be seen from Table7 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.Table 7
Denoising results comparison of swinging data.
Data/schemeRMSE (°/s)SNR (dB)Original signal0.487613.2743Scheme 10.199821.1113Scheme 20.089723.1142Scheme 30.079925.4896The 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 theZ-axis output signal as an example. The results are shown in Figure 15 and Table 7.Figure 15
Swinging data denoising results comparison of each scheme.It can be seen from Table8 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.Table 8
Performance parameters of vehicle sensors.
ParametersValueUnitUpdate rate300HzGyro bias stability0.007°/sGyro random walk2.4°/s/HzGyro resolution0.007mgAccelerometer bias stability0.2°/sAccelerometer random walk0.2°/s/HzAccelerometer resolution0.33mg


### 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 Figure16. 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.Figure 16
Vehicle experimental device.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 paperThe longitude and latitude error results of the three schemes are shown in Figure17, 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.Figure 17
Comparison of position errors for each scheme. (a) Longitude error. (b) Latitude error.
(a)(b)Table 9
Position errors’ parameter of each scheme.
MeanStandard deviationMax errorScheme 1Longitude error (m)9.24549.722834.0624Latitude error (m)5.36006.438923.038Scheme 2Longitude error (m)2.68003.219511.5190Latitude error (m)2.75862.24007.6928Scheme 3Longitude error (m)1.37931.11201.6800Latitude error (m)2.06891.68005.7696Based 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.


## 4.1. Static Data Denoising Experiment

As shown in Figure2, 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.Figure 2
Inertial platform test system.Table 2
Performance parameters of 3DM-10A.
ParametersValueUnitGyro bias stability0.03°/sGyro random walk0.03°/(s1/2)Accelerometer bias stability0.03°/sAccelerometer random walk0.03°/(s1/2)Resolution±0.01°/sFigure 3
Static data from gyroX-axis.Then, decompose and screen the above original signal. As shown in Figure4, after decomposition, 15 IMF components and one remainder R are obtained.Figure 4
The IMF components and remainder of the original signal.


### 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 j=6.Figure 5
Bhattacharyya distance between input signals and each component.


### 4.1.2. Analysis of Autocorrelation Function

In order to determine the second demarcation pointk, 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 IMF1−IMF10 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 IMF11−IMF15 has a certain regularity change, which can be determined as it is dominated by the signal.Figure 6
Autocorrelation functions of each IMF.In order to further accurately determine the second demarcation pointk, 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 IMF1−IMF10 are all less than 0.01; nevertheless, the variances of the autocorrelation function of IMF11−IMF15 are significantly larger than 0.01 and increase exponentially. Therefore, IMF6−IMF10 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 IMF6 as an example). Firstly, the IMF6 component signal is tested for zero-mean and normal distribution, as shown in Figures 8 and 9. After testing, it is concluded that the IMF6 component satisfies the modelling conditions. The same method tests the IMF8−IMF10 in turn, all of which meet the modelling conditions.Figure 7
Variances of ACF of each IMF.Figure 8
Normality test of IMF6.Figure 9
Mean and mean square error changing chart. (a) Mean change. (b) Variance of mean square error.
(a)(b)Next, model identification and model ordering of component signals are performed, and the autocorrelation function curve and the partial correlation function curve ofIMF6 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 IMF6 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.Figure 10
Model identification of IMF6. (a) Autocorrelation coefficient graph. (b) Partial autocorrelation coefficient graph.
(a)(b)Figure 11
BIC criterion for ordering determination.In the same way,IMF6−IMF10 are modelled, respectively. After a series of modelling steps and applicability tests, the final modelling results are shown in Table 3.Table 3
Modelling results ofIMF6−IMF10.
Component signalModelφ1φ2φ3φ4φ5φ6IMF6AR (4)−3.47994.4906−2.53360.5232IMF7AR (6)−2.42961.5645−0.0048−0.0099−0.0673−0.0528IMF8AR (4)−2.43961.44870.4232−0.4322IMF9AR (4)−2.09451.09460.09500.0950IMF10AR (4)−2.00821.00820.0088−0.0088Finally, 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 articleThe denoising results of each scheme are shown in Figure12.Figure 12
Denoising results of each scheme.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:(20)RMSE=1n∑i=1Nxi−xi′2,(21)SNR=10log∑i=1Nxi2∑i=1Nxi−xi′2.In equations (20) and (21), xi is the original signal; xi′ is the denoising signal; and N is the length of the signal.The calculation results are shown in Table4.Table 4
Denoising results comparison of each scheme.
Data/schemeRMSE (°/s)SNR (dB)Original signal0.349711.9505Scheme 10.193016.4550Scheme 20.155321.6059Scheme 30.100122.9756It can be seen from Table5 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.Table 5
Comparison of signal error term coefficients of each scheme.
Data/SchemeQ (°)N (°/h1/2)B (°/h)K (°/h3/2)R (°/h2)Original signal0.00270.010918.895026.01616.0252Scheme 10.00230.008714.761421.35085.0599Scheme 20.00200.005614.131721.97755.2616Scheme 30.00080.001810.714116.79333.9505


## 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 j=6.Figure 5
Bhattacharyya distance between input signals and each component.


## 4.1.2. Analysis of Autocorrelation Function

In order to determine the second demarcation pointk, 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 IMF1−IMF10 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 IMF11−IMF15 has a certain regularity change, which can be determined as it is dominated by the signal.Figure 6
Autocorrelation functions of each IMF.In order to further accurately determine the second demarcation pointk, 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 IMF1−IMF10 are all less than 0.01; nevertheless, the variances of the autocorrelation function of IMF11−IMF15 are significantly larger than 0.01 and increase exponentially. Therefore, IMF6−IMF10 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 IMF6 as an example). Firstly, the IMF6 component signal is tested for zero-mean and normal distribution, as shown in Figures 8 and 9. After testing, it is concluded that the IMF6 component satisfies the modelling conditions. The same method tests the IMF8−IMF10 in turn, all of which meet the modelling conditions.Figure 7
Variances of ACF of each IMF.Figure 8
Normality test of IMF6.Figure 9
Mean and mean square error changing chart. (a) Mean change. (b) Variance of mean square error.
(a)(b)Next, model identification and model ordering of component signals are performed, and the autocorrelation function curve and the partial correlation function curve ofIMF6 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 IMF6 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.Figure 10
Model identification of IMF6. (a) Autocorrelation coefficient graph. (b) Partial autocorrelation coefficient graph.
(a)(b)Figure 11
BIC criterion for ordering determination.In the same way,IMF6−IMF10 are modelled, respectively. After a series of modelling steps and applicability tests, the final modelling results are shown in Table 3.Table 3
Modelling results ofIMF6−IMF10.
Component signalModelφ1φ2φ3φ4φ5φ6IMF6AR (4)−3.47994.4906−2.53360.5232IMF7AR (6)−2.42961.5645−0.0048−0.0099−0.0673−0.0528IMF8AR (4)−2.43961.44870.4232−0.4322IMF9AR (4)−2.09451.09460.09500.0950IMF10AR (4)−2.00821.00820.0088−0.0088Finally, 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 articleThe denoising results of each scheme are shown in Figure12.Figure 12
Denoising results of each scheme.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:(20)RMSE=1n∑i=1Nxi−xi′2,(21)SNR=10log∑i=1Nxi2∑i=1Nxi−xi′2.In equations (20) and (21), xi is the original signal; xi′ is the denoising signal; and N is the length of the signal.The calculation results are shown in Table4.Table 4
Denoising results comparison of each scheme.
Data/schemeRMSE (°/s)SNR (dB)Original signal0.349711.9505Scheme 10.193016.4550Scheme 20.155321.6059Scheme 30.100122.9756It can be seen from Table5 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.Table 5
Comparison of signal error term coefficients of each scheme.
Data/SchemeQ (°)N (°/h1/2)B (°/h)K (°/h3/2)R (°/h2)Original signal0.00270.010918.895026.01616.0252Scheme 10.00230.008714.761421.35085.0599Scheme 20.00200.005614.131721.97755.2616Scheme 30.00080.001810.714116.79333.9505


## 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 Figure13, 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 101 magnitude.Figure 13
Allan variances comparison of each scheme.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 Table5. 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 theX-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.Figure 14
Rotating data denoising results comparison of each scheme.Table 6
Denoising results comparison of rotating data.
Data/schemeRMSE (°/s)SNR (dB)Original signal0.019810.9428Scheme 10.019316.2687Scheme 20.019120.9575Scheme 30.008332.4631It can be seen from Table7 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.Table 7
Denoising results comparison of swinging data.
Data/schemeRMSE (°/s)SNR (dB)Original signal0.487613.2743Scheme 10.199821.1113Scheme 20.089723.1142Scheme 30.079925.4896The 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 theZ-axis output signal as an example. The results are shown in Figure 15 and Table 7.Figure 15
Swinging data denoising results comparison of each scheme.It can be seen from Table8 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.Table 8
Performance parameters of vehicle sensors.
ParametersValueUnitUpdate rate300HzGyro bias stability0.007°/sGyro random walk2.4°/s/HzGyro resolution0.007mgAccelerometer bias stability0.2°/sAccelerometer random walk0.2°/s/HzAccelerometer resolution0.33mg


## 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 Figure16. 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.Figure 16
Vehicle experimental device.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 paperThe longitude and latitude error results of the three schemes are shown in Figure17, 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.Figure 17
Comparison of position errors for each scheme. (a) Longitude error. (b) Latitude error.
(a)(b)Table 9
Position errors’ parameter of each scheme.
MeanStandard deviationMax errorScheme 1Longitude error (m)9.24549.722834.0624Latitude error (m)5.36006.438923.038Scheme 2Longitude error (m)2.68003.219511.5190Latitude error (m)2.75862.24007.6928Scheme 3Longitude error (m)1.37931.11201.6800Latitude error (m)2.06891.68005.7696Based 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.



---

*Source: 3019152-2020-12-08.xml*