CN110542406B - Improved Gyroscope Signal Denoising Method Based on EMD-MPF - Google Patents

Abstract

The invention discloses an EMD-MPF (empirical mode decomposition-multi-point filter) based improved gyroscope signal denoising method, which comprises the following steps of: decomposing a noisy signal of a gyroscope into an intrinsic mode function and a residual signal; the intrinsic mode function is selected in order through the two determined identification parameters, and the intrinsic mode function is divided into a noise intrinsic mode function, a mixed intrinsic mode function and an information intrinsic mode function; abandoning the noise intrinsic mode function, reserving the information intrinsic mode function, and carrying out noise reduction processing on the mixed intrinsic mode function; and performing signal reconstruction on the mixed intrinsic mode function subjected to noise reduction and the reserved information intrinsic mode function to obtain a denoised gyroscope signal. The method has the advantages of realizing accurate noise signal removal and improving navigation accuracy.

Description

Improved Gyroscope Signal Denoising Method Based on EMD-MPF

technical field

The invention relates to the field of signal denoising, in particular to an improved gyroscope signal denoising method based on EMD-MPF.

Background technique

Due to the development of Micro-Electro-Mechanical System (MEMS), smaller and cheaper accelerometers and gyroscopes can be manufactured to form a lower-cost Inertial Measurement Unit (IMU). It has the advantages of small size, light weight, low power consumption, and mass production, which makes the inertial navigation or the positioning and navigation technology combined with the inertial measurement unit develop rapidly, and is widely used in many fields such as transportation, measurement, aerospace and so on. The size of the error of the MEMS gyroscope will determine the actual performance of the system to a large extent, so the suppression of various errors of the gyroscope has important practical significance. The possible errors of the gyroscope include installation error, lever arm error, scale factor error, etc.; while the noise mainly includes two aspects: system noise and measurement noise. Among them, the random drift of the MEMS gyroscope is an important source of noise. The random drift has a great influence on the accuracy of the gyroscope itself and the entire navigation system. Due to the long-term operation of the device, the slight drift is accumulated or After the system information is merged, it will lead to a large error. Therefore, how to reduce noise has become one of the most important issues in MEMS gyro signal processing.

The de-noising methods of MEMS gyroscope signals are mainly divided into two categories: modeling de-noising methods and non-modeling de-noising methods. First, mature random drift models are widely used or ARMA is used for modeling, and then Kalman is used for modeling. Filtering or various improved Kalman filtering techniques to filter the noise in the signal; second, directly use the wavelet decomposition or empirical mode decomposition (Empirical Mode Decomposition, EMD) method to analyze the signal and remove the noise in it to get This extracts a valid signal.

Traditional denoising methods generally use low-pass, high-pass or band-pass filters to remove noise. For example, using a digital low-pass filter to eliminate high-frequency noise in a strapdown Inertial Navigation System (INS); using an FIR low-pass filter to filter and suppress the measurement noise of MEMS sensors; applying a high-pass filter for gyroscope jitter A phenomenon with very high attenuation in the frequency band and a linear phase response in the passband, filtering the corresponding signal jitter noise. However, the traditional denoising method is based on the classical filtering theory, which is only suitable for the case where the signal and the noise do not overlap, so the filtered frequency component of the actual signal still contains white noise. For modern filtering theory, Gaussian filtering, Wiener filtering, Kalman filtering, etc. are its representative methods. These methods all require the use of prior knowledge of signal and noise statistics. For example, the adaptive robust Kalman filtering method is used, and Using the weighted covariance of the innovative sequence to adjust the measurement noise covariance matrix to reduce the random noise in the gyroscope; applying a random weighting method based on the anti-outlier recursive least squares (RLS) adaptive filtering algorithm To reduce the noise component in the output of the MEMS gyroscope; use the adaptive extended Kalman filter method to compensate the sideslip angle estimation, thereby reducing the potentially large drift estimation error caused by the inertial sensor offset. These methods overcome the shortcomings of classical filtering theory, but in the process of use, they need to give a certain maximum standard, and need to know the prior statistical information of noise, and these methods are usually difficult or impossible to obtain in practice, generally only Test an optimal estimate instead. At the same time, these methods need to accurately establish the model, and the deviation of the model will also affect the final denoising effect.

On this basis, a non-modeling denoising method has been developed. Wavelet transform (WT) is a new transform analysis method. It inherits and develops the idea of localization of short-time Fourier transform, and at the same time overcomes the shortcomings of window size not changing with frequency, so it has good multi-resolution. The analysis feature can subdivide the time at high frequency and frequency at low frequency, so as to achieve effective analysis of signals in different frequency bands and times. The error characteristics of the gyroscope are analyzed by the improved wavelet transform method, and a new soft threshold function is established to suppress the measurement noise of the gyroscope. The second generation wavelet transform is applied to de-noise the MEMS gyroscope. However, the wavelet transform analysis method needs to choose different decomposition layers and wavelet bases according to different situations, and also needs to choose an appropriate threshold for denoising, which is very inconvenient to use; at the same time, the WT method is more suitable for linear and steady-state signals. Therefore, there is a problem of inaccurate removal of noise information, resulting in low navigation accuracy.

SUMMARY OF THE INVENTION

The purpose of the present invention is to provide an improved method for denoising gyroscope signals based on EMD-MPF in view of the above problems, so as to achieve the advantages of accurately removing noise signals and improving navigation accuracy.

To achieve the above purpose, the technical solution adopted in the embodiment of the present invention is:

An improved gyroscope signal denoising method based on EMD-MPF, comprising:

Decompose the noisy signal of the gyroscope into eigenmode functions and residual signals;

The order of the eigenmode function is selected by the determined two identification parameters, and the eigenmode function is divided into noise eigenmode function, mixed eigenmode function and information eigenmode function;

discarding the noise eigenmode function, retaining the information eigenmode function, and performing noise reduction processing on the mixed eigenmode function;

Signal reconstruction is performed on the mixed eigenmode function after noise reduction processing and the retained information eigenmode function, thereby obtaining a denoised gyroscope signal.

As a specific implementation manner of the embodiment of the present invention, decomposing the noisy signal of the gyroscope into the eigenmode function and the residual signal is as follows:

The noisy signal of the gyroscope is decomposed into several eigenmode functions and a residual signal.

As a specific implementation manner of the embodiment of the present invention, the reconstruction formula in the signal reconstruction of the mixed eigenmode function after noise reduction processing and the retained information eigenmode function is:

表示噪声信号；M₁和M₂分别为标识参数，

表示混合本征模态函数；h⁽ⁱ⁾(t)表示信息本征模态函数；r_L(t)表示残差信号，MPF表示改进的PF方法。in,

represents the noise signal; M ₁ and M ₂ are the identification parameters, respectively,

represents the mixed eigenmode function; h ⁽ⁱ⁾ (t) represents the information eigenmode function; r _L (t) represents the residual signal, and MPF represents the improved PF method.

As a specific implementation of the embodiment of the present invention, the MPF is an improved PF method, including:

Denote the probability distribution of predicted particle samples as p(x _i ), 1≤i≤N, and denote the probability distribution of particle samples x _T with weight w _i >w ₀ as p ₁ (x _i ), 1≤i≤T , which is p, the probability distribution of the remaining particle sample x _NT is recorded as q(x _i ), NT<i≤N, which is q, w _i is the weight, w ₀ is the weight threshold, N and T are both constants ;

Calculate the K-L divergence between different particle samples as follows:

For the particle sample with small weight, generate NT random numbers u _j , whose range is u _j ∈(w ₀ ,1], update the weight of the particle sample, D _KL is the KL divergence:

将x_i+1距权重最大的粒子的距离记作S，进行重采样产生新的粒子如下：Sort the particle samples x _NT of the updated weights from small to large weights, that is,

Denote the distance of x _i+1 from the particle with the largest weight as S, and perform resampling to generate a new particle as follows:

表示重采样产生的新粒子。Integrate the two parts of particles, normalize the weights, and complete the continuous iteration of the PF method, x _i represents the ith particle, x _i+1 represents the i+1th particle, and x _N-T+1 represents the N-Tth particle +1 particle, x _N means the Nth particle,

Represents new particles resulting from resampling.

As a specific implementation manner of the embodiment of the present invention, among the two determined identification parameters, the two identification parameters are M ₁ and M ₂ respectively, and selecting and determining the identification parameter M ₁ includes: using an autocorrelation function The variance determines the identification parameter M ₁ .

As a specific implementation manner of the embodiment of the present invention, the use of the autocorrelation function variance to determine the identification parameter M ₁ includes:

According to the eigenmode function signals of different orders, the corresponding autocorrelation functions are calculated;

Calculate the mean square error of Euclidean distance between two autocorrelation functions according to the autocorrelation function of the continuous eigenmode function, and obtain the calculation result MSED;

The identification parameter M ₁ is determined according to the MSED.

As a specific implementation of the embodiment of the present invention, the identification parameter M ₁ is associated with the numerical value k and the numerical value L,

The 1≤k≤[2L/3], [2L/3] represents the largest integer not exceeding 2L/3, thereby ensuring that the value of M ₁ is too large to L, and both k and L are natural numbers.

As a specific implementation manner of the embodiment of the present invention, determining the identification parameter M ₂ includes:

Calculates the Bhattacharya distance between the probability density function of the original signal and the probability density function of each eigenmode function, where 0≤BC≤1 and 0≤BD≤∞, BD is the Bhattacharya distance , BC is the Battacharya coefficient;

Calculate, based on the power spectral density function, the separability of the Battacaria distance data under different characteristics between the probability density function of the original signal and the probability density function of each of the eigenmode functions;

The identification parameter M ₂ is determined based on the separability.

As a specific implementation manner of the embodiment of the present invention, the eigenmode function is selected by the order of the two determined identification parameters, and the eigenmode function is divided into noise eigenmode function and mixed eigenmode function. modal functions and information eigenmode functions, including:

When 1<i<M ₁ , it is the noise eigenmode function;

When M ₁ <i < M ₂ , it is a mixed eigenmode function;

When M ₂ <i<L, it is the information eigenmode function, and i represents the order of the eigenmode function

L represents the total order of the eigenstate function.

Embodiments of the present invention have the following beneficial effects:

Below, through the accompanying drawings and embodiments, the noisy signal is decomposed into an eigenmode function and a residual signal, and then the eigenmode function is divided into a noise eigenmode function and a mixed eigenmode through two identification parameter eigenmode functions. function and information eigenmode function; then the noise eigenmode function is discarded, and the mixed eigenmode function and information eigenmode function after noise reduction are reconstructed to accurately remove the noise signal and achieve The advantage of improving navigation accuracy.

Description of drawings

1 is a flowchart of an improved method for denoising gyroscope signals based on EMD-MPF according to an embodiment of the present invention;

2 is a block diagram of an EMD-MPF denoising method according to an embodiment of the present invention;

3 is the static data of the X-axis of the MEMS gyroscope according to the embodiment of the present invention;

Fig. 4 is the EMD exploded view of the X-axis gyro data described in the embodiment of the present invention;

Fig. 5 is the autocorrelation function of each order decomposition signal according to the embodiment of the present invention;

Fig. 6 is the MSED distribution situation described in the embodiment of the present invention;

7 is a probability density function of a decomposed signal according to an embodiment of the present invention;

Fig. 8 is the distribution situation of Barthel distance according to the embodiment of the present invention;

Fig. 9 is the filtering result of the hybrid IMF using MPF according to the embodiment of the present invention;

Fig. 10 is the X-axis static test result comparison diagram according to the embodiment of the present invention;

FIG. 11 is the result of Allen variance analysis according to the embodiment of the present invention.

Detailed ways

The preferred embodiments of the present invention will be described below with reference to the accompanying drawings. It should be understood that the preferred embodiments described herein are only used to illustrate and explain the present invention, but not to limit the present invention.

The EMD-MPF denoising method for MEMS gyro signals in this embodiment. The original signal is decomposed into several eigenmode functions through empirical mode decomposition, and then the IMF is divided into noise IMF, mixed IMF and information IMF through analysis and order selection. denoised signal.

The EMD method can be an adaptive signal decomposition algorithm for nonlinear and non-stationary signals. It not only breaks through the limitations of Fourier transform, but also does not have the problem of preselecting wavelet basis functions like wavelet transform. It has good time-frequency resolution and Adaptive. Based on empirical mode decomposition, a fully data-driven method is used to reconstruct the signal part to achieve the effect of signal filtering; the time-frequency domain FOG signal denoising method based on empirical mode decomposition and Allan variance is used to effectively improve the navigation accuracy .

As shown in Figure 1, an improved gyroscope signal denoising method based on EMD-MPF includes:

S101: Decompose the noisy signal of the gyroscope into an eigenmode function and a residual signal;

S102: Select the order of the eigenmode function (IMF) through the determined two identification parameters, and divide the eigenmode function into a noise eigenmode function, a mixed eigenmode function and an information eigenmode function ;

S103: Discard the noise eigenmode function, retain the information eigenmode function, and perform noise reduction processing on the mixed eigenmode function;

S104: Perform signal reconstruction on the mixed eigenmode function after noise reduction processing and the retained information eigenmode function, so as to obtain a denoised gyroscope signal.

As an optional implementation manner of the embodiment of the present invention, decomposing the noisy signal of the gyroscope into the eigenmode function and the residual signal is as follows:

The noisy signal of the gyroscope is decomposed into several eigenmode functions and a residual signal.

As a specific implementation manner of the embodiment of the present invention, the reconstruction formula in the signal reconstruction of the mixed eigenmode function after noise reduction processing and the retained information eigenmode function is:

表示噪声信号；M₁和M₂分别为标识参数，

表示混合本征模态函数；h⁽ⁱ⁾(t)表示信息本征模态函数；r_L(t)表示残差信号，MPF表示改进的PF方法。in,

represents the noise signal; M ₁ and M ₂ are the identification parameters, respectively,

represents the mixed eigenmode function; h ⁽ⁱ⁾ (t) represents the information eigenmode function; r _L (t) represents the residual signal, and MPF represents the improved PF method.

As an optional implementation manner of the embodiment of the present invention, the MPF is an improved PF method, including:

Denote the probability distribution of predicted particle samples as p(x _i ), 1≤i≤N, and denote the probability distribution of particle samples x _T with weight w _i >w ₀ as p ₁ (x _i ), 1≤i≤T , which is p, the probability distribution of the remaining particle sample x _NT is recorded as q(x _i ), NT<i≤N, which is q, w _i is the weight, w ₀ is the weight threshold, N and T are both constants ;

Calculate the K-L divergence between different particle samples as follows:

The introduction of K-L divergence can better ensure the effectiveness and diversity of particle sample sets;

For particle samples with small weights, generate NT random numbers u _j , whose range is u _j ∈(w ₀ ,1], u _j is generated in the range (w ₀ ,1] using the random function random in the program, D _KL is the KL divergence, and the weight of the particle sample is updated:

将x_i+1距权重最大的粒子的距离记作S，进行重采样产生新的粒子如下：Sort the particle samples x _NT of the updated weights from small to large weights, that is,

Denote the distance of x _i+1 from the particle with the largest weight as S, and perform resampling to generate a new particle as follows:

表示重采样产生的新粒子。Integrate the two parts of particles, normalize the weights, and complete the continuous iteration of the PF method, x _i represents the ith particle, x _i+1 represents the i+1th particle, and x _N-T+1 represents the N-Tth particle +1 particle, x _N means the Nth particle,

Represents new particles resulting from resampling.

As an optional implementation manner of the embodiment of the present invention, among the two determined identification parameters, the two identification parameters are M ₁ and M ₂ respectively, and selecting and determining the identification parameter M ₁ includes: using autocorrelation The function variance determines the identity parameter M ₁ .

As an optional implementation manner of the embodiment of the present invention, the use of the autocorrelation function variance to determine the identification parameter M ₁ includes:

According to the eigenmode function signals of different orders, the corresponding autocorrelation functions are calculated;

Calculate the mean square error of Euclidean distance between two autocorrelation functions according to the autocorrelation function of the continuous eigenmode function, and obtain the calculation result MSED;

The identification parameter M ₁ is determined according to the MSED.

As an optional implementation manner of the embodiment of the present invention, the identification parameter M ₁ is associated with the numerical value k and the numerical value L,

The 1≤k≤[2L/3], [2L/3] represents the largest integer not exceeding 2L/3, thereby ensuring that the value of M ₁ is too large to L, and both k and L are natural numbers.

As an optional implementation manner of the embodiment of the present invention, determining the identification parameter M ₂ includes:

Calculates the Bhattacharya distance between the probability density function of the original signal and the probability density function of each eigenmode function, where 0≤BC≤1 and 0≤BD≤∞, BD is the Bhattacharya distance , BC is the Battacharya coefficient;

Calculate, based on the power spectral density function, the separability of the Battacaria distance data under different characteristics between the probability density function of the original signal and the probability density function of each of the eigenmode functions;

The identification parameter M ₂ is determined based on the separability.

As an optional implementation manner of the embodiment of the present invention, the order selection is performed on the eigenmode function by using the determined two identification parameters, and the eigenmode function is divided into a noise eigenmode function and a mixed eigenmode function. Modal functions and informative eigenmode functions, including:

When 1<i<M ₁ , it is the noise eigenmode function;

When M ₁ <i < M ₂ , it is a mixed eigenmode function;

When M ₂ <i<L, it is the information eigenmode function, and i represents the order of the eigenmode function

L represents the total order of the eigenstate function.

In a specific application scenario,

1. Improved EMD-MPF method:

EMD is an adaptive spatiotemporal analysis method suitable for dealing with non-stationary nonlinear sequences. This method decomposes the signal into a series of eigenmode functions (IMFs) that do not deviate from the time domain through an iterative process of screening, so it has good time-frequency resolution and adaptability, and can perfectly reconstruct the original signal, while at the same time Has the potential to highlight fine geological formations in the signal that might otherwise be overlooked. At present, the EMD method has been successfully applied in the field of reducing signal noise. Among them, the eigenmode function must satisfy the following two requirements: (1) the number of extrema and the number of zero crossings are at most equal or different by one; (2) at any point, the mean value of the envelope defined by the local maxima and the The envelope defined by the local minimum is zero. For the IMF signal h ⁽ⁱ⁾ (t) and a residual signal r _L (t) obtained by decomposing the noisy signal x(t), the number of L can be reconstructed as follows:

The idea of Particle Filter (PF) is based on Monte Carlo methods, which uses particle sets to represent probability and can be used in any form of state space model. This method is a process of approximately representing the probability density function by finding a set of random samples propagating in the state space, replacing the integral operation with the sample mean, and then obtaining the minimum variance estimation of the system state. Because this method adds another layer of importance sampling idea to the Monte Carlo method, a set of samples can be used to approximate the posterior probability distribution of the system, and then the state of the nonlinear system can be estimated using this approximate result. The general steps for the PF method are as follows:

Step 1: Initialize the state quantity: the initial state is simulated in space by a large number of particles according to the prior probability distribution.

Step 2, prediction stage: According to the state transition equation, each particle gets a prediction particle.

Step 3: Correction phase: Evaluate the predicted particles. The particles whose state is closer to the real state are given a larger weight.

Step 4, resampling stage: Screen the particles according to their weights. In the screening process, not only the particles with larger weights, but also a part of the particles with smaller weights should be retained.

Step 5, filtering stage: Substitute the resampled particles into the new predicted particles of the state transition equation, that is, return to step 2 for looping.

The traditional EMD method decomposes the original signal to obtain a series of IMFs with frequencies from high to low. It is generally believed that the noise basically exists in the high frequency band, so some IMFs in the high frequency band are directly deleted and then recombined to obtain a denoised signal. Although the PF method can get rid of the restriction that the random quantity must satisfy the Gaussian distribution when solving the nonlinear filtering problem, the problem is that it needs a large number of samples to well approximate the posterior probability density of the system, that is, it needs to use a large number of samples. There are a large number of effective particles with diversity, and the more complex the actual environment, the more particles are required, which leads to the higher complexity of the algorithm.

2. Improved IMF order selection method:

The EMD method can decompose the noisy signal into several IMF components of different scales and a residual residue. Usually, it is known that the noise signal is mainly concentrated in the high frequency part, and the composition of the IMF can be generally divided into three components, that is, the noise IMF, the mixed IMF and the information IMF. Now, for signal denoising rather than pure complete reconstruction, the main problem is how to distinguish correlated IMFs from uncorrelated IMFs, that is, how to find the demarcation point of different IMF components is the key problem of this method for denoising. Based on the above, the noise signal is re-described as follows:

表示噪声信号；M₁和M₂分别为下文描述的两个分界点；

表示混合IMF部分；h⁽ⁱ⁾(t)表示信息IMF部分；r_L(t)表示剩余残差信号。in,

represents the noise signal; M ₁ and M ₂ are the two demarcation points described below;

represents the mixed IMF part; h ⁽ⁱ⁾ (t) represents the information IMF part; r _L (t) represents the residual residual signal.

2.1. About the determination of the dividing point M1:

The current problem is how to determine the two demarcation points M ₁ and M ₂ . For the first parameter, use the variance of the autocorrelation function to determine. In practical applications, the actual signal components in the noise signal must not be known. In other words, the noise in a noisy signal can never be fully known and removed in theory or in practice. Therefore, the parameters can neither be determined by minimizing the mean square error of the target signal and the noise signal nor by obtaining the respective signal-to-noise ratio of each component. The reconstructed signal of the IMF is as follows:

The autocorrelation function can determine the correlation degree of the signal in different time domains. According to the IMF signals of different orders, the corresponding autocorrelation functions are calculated as follows:

Then, the mean square error of Euclidean distance between the two is calculated according to the autocorrelation function of the continuous IMF, where Pea is the correlation coefficient, and the calculation result of the last formula (5) is abbreviated as MSED, and the calculation is as follows:

According to the statistical characteristics of random noise, the autocorrelation function of random noise is the largest only at the zero point, and its autocorrelation function will quickly decay to the vicinity of the zero value at other times. The use of Euclidean distance is a common classification basis, and the introduction of the variance of the autocorrelation function not only considers the dimensionless situation under the autocorrelation function, but also considers the time domain between different components by relying on the autocorrelation function. To a certain extent, the correlation on the Euclidean distance reduces the interference of different components in the Euclidean distance under a single feature.

Finally, the first indicator parameter M ₁ is obtained, which is determined as follows:

Among them, the value of k is set between 1 and [2L/3], and [2L/3] represents the largest integer not exceeding 2L/3. This setting is to ensure that the value of M ₁ is too large to L.

Through the above analysis, according to the autocorrelation characteristics, if it is the IMF component of the noise main body, the random noise has a maximum value at the zero point, and will rapidly decrease to zero approximately in other time domains; if it is the IMF component of mixed information, then the random noise will Larger fluctuations also occur in other time domains. At the same time, according to the difference between IMF components in the continuous time domain calculated and described by MSED, the smaller the MSED value, the smaller the corresponding h ⁽ⁱ⁾ (t) component; There is a lot of useful information for MEMS gyroscopes, so the corresponding h(t) in them has high energy. In summary, the calculated value of MSED can be used as the first representation parameter.

2.2. About the determination of the dividing point M2:

The demarcation point M ₂ is mainly used for identifying and distinguishing information IMFs and mixed IMFs. To determine the second identification parameter, the Bhattacharyya distance (BD) between the probability density function (PDF) of the original signal and the PDF of each IMF is introduced. The Bhattacharyya distance and the Bhattacharyya coefficient are named after A. Bhattacharya, a statistician who worked at the Indian Institute of Statistics in the 1930s. Some scholars have proposed to use Euclidean distance as a separation marker, but Euclidean distance will be significantly affected by the sample dimension; the Mahalanobis distance method considers the relationship between various features and eliminates the difference between variables. The correlation interference between the two is an indirect distinction method, and the Mahalanobis distance may exaggerate the effect of the variable with small changes, thus reducing the accuracy of the dividing point. The Bavarian distance is often used in the probability distribution between two discrete samples to judge the separability between measurement classes, which not only eliminates the interference of the correlation between variables, but also is a more direct method of distinction. The Battacharia distance between discrete probability distributions p and q is defined as follows:

BD(p,q)=-ln[BC(p,q)] (9),

Among them, BC is called the Bhattacharyya coefficient; and 0≤BC≤1 and 0≤BD≤∞ are satisfied.

The power spectral density function defines how the power of a signal or time series is distributed with frequency and is a good description of the characteristics of the data. Therefore, the separability of the Battacharia distance data between the original signal and the PDF of each IMF under different characteristics, is calculated as:

Analyzing the above formula, the more noise mixed in the IMF, the lower the separability compared with the original signal; if it is a "pure" information IMF, the separability will be higher compared with the original signal . Therefore, the selection method of _M2 is as follows:

There is a point of view, if the Batacharya distance between the two PDFs of the IMF component and the original signal is smaller, the separability will be smaller, that is, the more noise is contained in the IMF component; on the contrary, the difference between the two PDFs will be smaller. The larger the Bhattacharya distance, the smaller the similarity and the more separable, that is, the less noise contained in the IMF component. Therefore, the next point of the maximum value of the Bhattacharya distance is regarded as the dividing point, and the mixed IMF component and the information IMF component can be effectively distinguished. Usually, the IMF order determined by the demarcation point M ₂ is higher than the order determined by the demarcation point M ₁ , but there may be a situation that most of the IMF orders decomposed by the noisy signal contain noise , resulting in M ₁ behind M ₂ , so formula (10) is improved as follows to ensure that the IMF order determined by the identification parameter M ₂ is higher than the order determined by M ₁ :

To sum up, the two identification parameters M ₁ and M ₂ decompose the IMF component obtained by the EMD method into three parts, wherein, when 1<i<M ₁ , it is the noise IMF component and is discarded; M ₁ <i< When M ₂ , it is considered as a mixed IMF component, and the following improved particle filter is used for further noise reduction processing; when M ₂ <i<L, it is considered as a time-information IMF component and retained.

3. Improved PF method:

The traditional resampling process is that in each iteration of the algorithm, the particles are resampled according to the particle weights, and the particles with smaller weights are discarded and replaced by particles with larger weights. This resampling process can better solve the original degradation phenomenon in particle filtering, but it also brings about the problem of sample exhaustion. So combined with KL divergence to improve, an innovative resampling process is proposed. The weight threshold of the resampling process is set to w ₀ , the number of particles corresponding to the weight higher than the threshold is T, and the total number of particles is N, then the traditional resampling process is to retain the particles with a weight higher than w ₀ , and then resample according to the weight Generate NT particles. This destroys the validity and diversity of particle sample sets. Now, the specific steps of the improved resampling process are as follows:

1) Denote the probability distribution of predicted particle samples as p(x _i ), 1≤i≤N; denote the probability distribution of particle samples x _T with weight w _i >w ₀ as p ₁ (x _i ), 1≤i ≤T; denote the probability distribution of the remaining particle samples x _NT as q(x _i ), NT<i≤N.

2) Calculate the K-L divergence between different particle samples as follows:

3) For the particle sample with small weight, generate NT random numbers u _j , whose range is u _j ∈(w ₀ ,1], and update the weight of the particle sample:

将x_i+1距权重最大的粒子的距离记作S，进行重采样产生新的粒子如下：4) Sort the updated weight particle samples x _NT from small to large, that is

Denote the distance of x _i+1 from the particle with the largest weight as S, and perform resampling to generate a new particle as follows:

5) Integrate the two parts of particles, normalize the weights, and complete the iterative operation of the PF method.

The improved resampling process not only re-evaluates the weight of some particles, but also uses the K-L divergence to increase the order of the re-evaluation process, so that the random distribution difference during the resampling process is minimized, and the particle size is maintained as much as possible. Effectiveness and diversity, thereby increasing the efficiency of the resampling process.

4. The noise reduction system process of the improved method:

Based on the above improvement process, the specific flow of the EMD-MPF denoising method proposed in this embodiment is described in Figure 2, and the execution process is divided into the following steps:

Step 1: Decompose the noisy signal x(t) of the MEMS gyroscope according to the EMD method to obtain several IMFs and a residual signal.

Step 2: Perform order selection on the IMF through two identification parameters M ₁ and M ₂ , and divide the IMF into three parts: noise IMF, mixed IMF and information IMF.

Step 3: Discard the noise IMF; keep the information IMF; use the improved PF method to denoise the mixed IMF.

Step 4: Perform signal reconstruction on the mixed IMF and information IMF after noise reduction processing to obtain the final complete denoising result. The signal reconstruction formula is as follows:

In the improved gyro signal denoising method based on EMD-MPF, the noise in the noise signal cannot be completely known and removed in theory or practice. Therefore, the parameters can neither be determined by minimizing the mean square error of the target signal and the noise signal nor by obtaining the respective signal-to-noise ratio of each component. Therefore, according to the statistical characteristics of random noise, the autocorrelation function of random noise is the largest only at the zero point, and its autocorrelation function will rapidly decay to the vicinity of the zero value at other times. The use of Euclidean distance is a common classification basis, and the introduction of the variance of the autocorrelation function not only considers the dimensionless situation under the autocorrelation function, but also considers the time domain between different components by relying on the autocorrelation function. To some extent, it reduces the interference of Euclidean distance by different components under a single feature; according to the autocorrelation characteristics, if it is the IMF component of the main noise, the random noise has a maximum value at zero, and at other The time domain will rapidly drop to zero approximately; if it is the IMF component of mixed information, the random noise will fluctuate greatly in other time domains. After modification on the basis of the above-mentioned related contents, it is possible to quickly distinguish the IMF components mainly existing in noise and the information IMFs without knowing any actual signal components, so as to determine the index parameters between the two.

To determine the second identification parameter, the Bhattacharyya distance (BD) between the probability density function (PDF) of the original signal and the PDF of each IMF is introduced. Barthesian distance is often used in the probability distribution between two discrete samples to judge the separability between measurement classes, which not only excludes the interference of the correlation between variables, but also is a more direct method of distinction. If the Battacharya distance between the two PDFs of the IMF component and the original signal is smaller, the separability is smaller, that is, the more noise is contained in the IMF component; on the contrary, the Battacharya distance between the two PDFs is smaller. The larger the elegant distance, the smaller the similarity and the more separable, that is, the less noise contained in the IMF component. Therefore, the next point of the maximum value of the Bhattacharya distance is regarded as the dividing point, and the mixed IMF component and the information IMF component can be effectively distinguished.

Combined with the improvement of KL divergence, an innovative resampling process is proposed. The weight threshold of the resampling process is set to w ₀ , the number of particles corresponding to the weight higher than the threshold is T, and the total number of particles is N, then the traditional resampling process is to retain the particles with a weight higher than w ₀ , and then resample according to the weight Generate NT particles. This destroys the validity and diversity of particle sample sets.

This implementation is verified by experiments as follows:

The experimental environment is an inertial navigation rotation platform. The inertial measurement unit of the platform is a full attitude measurement sensing device IMU-300A. The device consists of a three-axis MEMS gyroscope, a three-axis MEMS accelerometer and a three-axis magnetoresistive type. The magnetometer consists of three parts. The device is fixed on a two-axis electric turntable, and the system can measure the absolute angular rate and acceleration values on the three axes X, Y, and Z respectively.

The effective position of the data in the experiment is determined by the calibration standard of the experimental platform. The experimental system adopts a dual-axis turntable. The system includes: a heading reference system and an angle sensor, which can accurately control the angular rotation of the turntable. The system can be calibrated automatically, and provides relatively accurate reference attitude angle information through the output attitude angle of the heading reference system.

In order to more clearly and comprehensively verify the effectiveness of the method proposed in this embodiment, two parts of the experimental content are set up later, the first part uses the static data of the MEMS gyroscope; the second part uses the dynamic data. In order to reflect the signal denoising effect of the improved method, the following four methods were used to analyze and compare the results. Method 1: using traditional particle filtering; method 2: using traditional EMD method; method 3: using improved PF method; method Fourth: use the improved gyroscope denoising method based on MED-MPF proposed in this embodiment.

3.1. Static data test of turntable

In the test preparation stage, after power-on, the experimental system equipped with the gyroscope is first reset to zero, and then it is preheated for about 5 minutes. The data sampling frequency is 100Hz. After the system is stabilized, the raw data output by the gyroscope will be collected; as shown in Figure 3 As shown, the static data of the gyro (take the signal data of the X-axis as an example) is collected, the data length is 49428, and the duration is about 8 minutes.

First, the noise signal is decomposed using the EMD method, and the 15th-order IMF signal and a signal residual as shown in Figure 4 are obtained.

According to the 15th-order IMF signal and the first-order residual signal obtained by decomposing, the autocorrelation function of each order IMF signal is calculated respectively, and the result is shown in Fig. 5. Simply analyze the autocorrelation function. The autocorrelation function describes the degree of correlation between the values of the random signal at any different time, and is the cross-correlation of the signal itself; because the random noise has the largest autocorrelation coefficient at the zero point, while the gyro signal There is a certain correlation between the time before and after the random noise, so that the association of random noise will not be interrupted immediately. The autocorrelation of the time before and after also reflects the correlation of random noise in the time domain to a certain extent. Therefore, the autocorrelation of the continuous IMF signal is used. The mean square error of the distance between the functions can judge the correlation between random noise and useful signal to a certain extent, so as to achieve the purpose of distinguishing the IMF order mainly containing random noise.

According to the results obtained in Figure 5 and the method for determining the M1 parameter, the results are shown in Figure 6.

The probability density function gives the probability of the signal taking different amplitudes, and describes the distribution of signal amplitudes well. The location of the binary demarcation point to distinguish the mixed IMF from the information IMF. The probability density function relationship between each order IMF and the original signal PDF is shown in Figure 7, and the determination result of the M2 identification point is shown in Figure 8.

According to the selection results of the boundary points in Figures 7 and 9, it can be seen that before the 5th-order IMF signal, it is basically a noise IMF, which can be directly discarded; after the 10th-order IMF signal is basically a large number of noise-free signals; , between the 5th order and the 10th order, that is, the selected mixed IMF part, the particle filtering method improved in this embodiment is used for this part of the signal, and the comparison result shown in Figure 10 is obtained, and finally the filtered The result is reconstructed with the information-dominated IMF signal to obtain the final processing result.

The filtering results and the information-dominated IMF signal are reconstructed, and compared and analyzed according to the four methods mentioned above. For the other three methods, the noise reduction amplitude is larger. In order to compare the experimental results more intuitively, Allen variance analysis was performed subsequently, and the signal-to-noise ratio and mean square error of the respective methods were calculated respectively.

Allen variance is a commonly used method for analyzing gyro noise. This method can easily characterize and identify the statistical characteristics of various error sources in the time domain. Therefore, the quantitative analysis of Allen variance can further verify this implementation. The effectiveness and applicability of the method for gyro denoising. By comparing the results in FIG. 11 , it can be found that using the method of this embodiment, the errors obtained by the Allen variance analysis are reduced to varying degrees.

RMSE can well reflect the fluctuation of sample data. The smaller the value, the smaller the fluctuation and the less noise it contains. Compared with the original data, the RMSE value of the method proposed in this example is about 91% lower than that of the original data, about 89% lower than that of method 1, about 81% lower than that of method 2, and about 63% lower than that of method 3, and the effect is very obvious. The signal-to-noise ratio more directly reflects the noise content in the signal. Generally, the larger the signal-to-noise ratio, the smaller the noise mixed in the signal. Compared with the original signal, the signal-to-noise of the method proposed in this embodiment is increased by 20.96dB, compared with the method 1, increased by 19.01dB, compared with the method 2, increased by 14.62dB, and improved by 8.67dB compared with the method 3, indicating that the method is suitable for the MEMS gyro signal. The denoising effect is more obvious. The analysis results are shown in Table 1.

Table 1 Error analysis of experimental results.

Aiming at the problem of random noise in the output of the MEMS gyroscope, this embodiment proposes a new denoising method combining EMD and improved PF. Analyze multiple IMF components obtained by EMD decomposition, obtain noise IMF, mixed IMF and information IMF through the calculation of two decomposition points, and then use the improved PF method to denoise the mixed IMF. The information IMF is reconstructed to obtain the final denoised signal. In the experiment part, the static data and dynamic data were used respectively, and the results were compared between the original data and the four denoising methods. noise, which greatly improves the signal quality.

This embodiment has the following advantages:

1. Compared with the traditional EMD method, a more detailed way of selecting the dividing point is proposed, and the hierarchy is divided into three parts, which is more targeted for random noise;

2. The second advantage is that the PF algorithm is improved in combination with K-L divergence, and an innovative resampling process is proposed, which keeps the effectiveness and diversity of the particle sample set;

3. The third advantage: Empirical mode decomposition decomposes the original signal into several eigenmode functions, and then divides the IMF into noise IMF, mixed IMF and information IMF through analysis and order selection. Finally, by eliminating noise and improving particle filtering method to obtain the final denoised signal.

Finally, it should be noted that the above descriptions are only preferred embodiments of the present invention, and are not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, for those skilled in the art, the The technical solutions described in the foregoing embodiments may be modified, or some technical features thereof may be equivalently replaced. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims (3)

1. a gyroscope signal denoising method based on EMD-MPF improvement, is characterized in that, comprises: Decompose the noisy signal of the gyroscope into eigenmode functions and residual signals; The order of the eigenmode function is selected by the determined two identification parameters, and the eigenmode function is divided into noise eigenmode function, mixed eigenmode function and information eigenmode function; discarding the noise eigenmode function, retaining the information eigenmode function, and performing noise reduction processing on the mixed eigenmode function; performing signal reconstruction on the mixed eigenmode function after noise reduction processing and the retained information eigenmode function, thereby obtaining a denoised gyroscope signal; The reconstruction formula in the signal reconstruction of the mixed eigenmode function after noise reduction processing and the retained information eigenmode function is:

in,

represents the noise signal; M ₁ and M ₂ are the identification parameters, respectively,

represents the mixed eigenmode function; h ⁽ⁱ⁾ (t) represents the information eigenmode function; r _L (t) represents the residual signal, and MPF represents the improved PF method; The MPF is an improved PF method, including: Denote the probability distribution of predicted particle samples as p(x _i ), 1≤i≤N, and denote the probability distribution of particle samples x _T with weight w _i >w ₀ as p ₁ (x _i ), 1≤i≤T , which is p, the probability distribution of the remaining particle sample x _NT is recorded as q(x _i ), NT<i≤N, which is q, w _i is the weight, w ₀ is the weight threshold, N and T are both constants ; Calculate the K-L divergence between different particle samples as follows:

For the particle sample with small weight, generate NT random numbers u _j , whose range is u _j ∈(w ₀ ,1], update the weight of the particle sample, D _KL is the KL divergence:

Sort the particle samples x _NT of the updated weights from small to large weights, that is,

Denote the distance of x _i+1 from the particle with the largest weight as S, and perform resampling to generate a new particle as follows:

Integrate the two parts of particles, normalize the weights, and complete the continuous iteration of the PF method, x _i represents the ith particle, x _i+1 represents the i+1th particle, and x _N-T+1 represents the N-Tth particle +1 particle, x _N means the Nth particle,

Represents new particles generated by resampling; Among the two determined identification parameters, the two identification parameters are M ₁ and M ₂ respectively, and selecting and determining the identification parameter M ₁ includes: determining the identification parameter M ₁ by using the variance of the autocorrelation function; The use of the variance of the autocorrelation function to determine the identification parameter M ₁ includes: According to the eigenmode function signals of different orders, the corresponding autocorrelation functions are calculated; Calculate the mean square error of Euclidean distance between two autocorrelation functions according to the autocorrelation function of the continuous eigenmode function, and obtain the calculation result MSED; Determine the identification parameter M ₁ according to the MSED; The identification parameter M ₁ is associated with the value k and the value L, The 1≤k≤[2L/3], [2L/3] represents the largest integer not exceeding 2L/3, thereby ensuring that the value of M ₁ is too large to L, and both k and L are natural numbers; Determine the identification parameter M ₂ , including: Calculates the Bhattacharya distance between the probability density function of the original signal and the probability density function of each eigenmode function, where 0≤BC≤1 and 0≤BD≤∞, BD is the Bhattacharya distance , BC is the Battacharya coefficient; Calculate, based on the power spectral density function, the separability of the Battacaria distance data under different characteristics between the probability density function of the original signal and the probability density function of each of the eigenmode functions; The identification parameter M ₂ is determined based on the separability. 2. the improved gyroscope signal denoising method based on EMD-MPF according to claim 1, is characterized in that, described by the noisy signal of gyroscope is decomposed into eigenmode function and residual signal as: The noisy signal of the gyroscope is decomposed into several eigenmode functions and a residual signal. 3. the improved gyroscope signal denoising method based on EMD-MPF according to claim 1, is characterized in that, described by the two identification parameters determined to carry out order selection to the eigenmode function, the eigenmode The functions are divided into noise eigenmode functions, mixed eigenmode functions and information eigenmode functions, including: When 1<i<M ₁ , it is the noise eigenmode function; When M ₁ <i < M ₂ , it is a mixed eigenmode function; When M ₂ <i<L, it is the information eigenmode function, and i represents the order of the eigenmode function L represents the total order of the eigenstate function.

Images