Gyro signal denoising method based on kalman filtering and visushrink threshold processing

Abstract

The present application relates to a gyro signal denoising method based on Kalman filtering and Visushrink threshold processing, comprising the following steps: establishing a time series model of a gyro signal; using adaptive anti-outlier Kalman filtering to denoise the gyro signal; using wavelet analysis to respectively perform Visushrink threshold processing on low-frequency components and high-frequency components of the Kalman filtered gyro signal; and performing wavelet reconstruction on the high-frequency and low-frequency signals of the threshold processed gyro signal. The present application provides an adaptive anti-outlier denoising scheme for gyro signals, which combines Kalman filtering and wavelets, can more effectively improve the accuracy of sensors and reduce errors.

Description

Technical Field

[0001] This invention belongs to the field of sensor signal processing, specifically relating to a gyroscope signal denoising method using Kalman filtering and Visushrink thresholding. Background Technology

[0002] MEMS (Microelectro Mechanical Systems) sensors, manufactured using microelectronics technology, are a new type of miniature inertial sensor widely used in low-cost strapdown inertial navigation systems. Before using a MEMS gyroscope to measure data, noise reduction processing is necessary to reduce its random drift error and improve sensor accuracy.

[0003] For denoising MEMS gyroscope signals, Kalman filtering or wavelet analysis is generally used. Kalman filtering is simple to calculate and easy to process, but in practice, due to the influence of various factors, the filter may diverge. Often, the divergence of the filter is suppressed by increasing the weight of the current measurement value. Moreover, interference data such as outliers in the signal can also affect the accuracy of the filter and affect the denoising effect.

[0004] Wavelet analysis is commonly used in signal processing. It involves decomposing the signal, thresholding its high-frequency components, and then reconstructing the low-frequency and processed high-frequency signals using wavelets to obtain a denoised signal. When a large amount of noise is located in the high-frequency components, wavelet denoising has a good noise suppression effect. However, the angular random walk and zero-bias instability in MEMS gyroscope signals often reside in the low-frequency components. Summary of the Invention

[0005] The technical problem of this invention is that MEMS gyroscope measurement data often have random drift errors, which need to be denoised. However, the accuracy of existing Kalman filtering methods is affected by interference data such as outliers in the signal. Existing wavelet denoising only performs thresholding on the high-frequency components of the signal, ignoring the angular random walk and zero-bias instability that are often in the low-frequency components of MEMS gyroscope signals.

[0006] The purpose of this invention is to solve the above problems by providing an adaptive anti-outlier denoising method for MEMS gyroscopes that integrates Kalman filtering and wavelet filtering. After performing adaptive anti-outlier Kalman filtering on the gyroscope signal, wavelet thresholding is then applied to both low-frequency and high-frequency components simultaneously, given a determined wavelet decomposition level. This method combines the advantages of Kalman filtering and wavelet denoising to better suppress gyroscope signal noise and improve the measurement accuracy of the gyroscope.

[0007] The technical solution of this invention is a gyroscope signal denoising method using Kalman filtering and Visushrink thresholding, comprising the following steps:

[0008] Step 1: Establish a time-series ARMA model of the MEMS gyroscope signal;

[0009] First, the trend and periodic terms are extracted from the original signal of the gyroscope. Then, a time series ARMA model is established from its residuals. The final prediction error FPE criterion is used to determine the order of the model and obtain the random drift error of the gyroscope.

[0010] Step 2: Use adaptive anti-outlier Kalman filtering to denoise the gyroscope signal;

[0011] Step 2.1: Establish the state equation and measurement equation of the Kalman filter based on the gyroscope random drift error model in Step 1;

[0012] Step 2.2: Initialize the adaptive anti-outlier Kalman filter. At the start of filtering, use the gyroscope error estimated by the gyroscope random drift error model to initialize the filter's state vector and covariance matrix.

[0013] Step 2.3: Perform adaptive anti-outlier Kalman filtering;

[0014] Step 3: Using wavelet analysis, threshold processing is performed on the low-frequency and high-frequency components of the Kalman-filtered gyroscope signal, respectively.

[0015] Step 3.1: Select the wavelet function;

[0016] Step 3.2: Determine the wavelet decomposition level, perform wavelet decomposition on the filtered gyroscope signal using the selected wavelet function, and calculate the peak value ratio of the high-frequency components in each layer;

[0017] Step 3.3: Perform Visushrink threshold processing on low-frequency and high-frequency signals;

[0018] Step 4: Perform wavelet reconstruction on the high-frequency and low-frequency signals of the gyroscope signal after threshold processing.

[0019] Furthermore, in step 1, the calculation model for the gyroscope's random drift error is as follows:

[0020] m k+1 =am k +ξ k+1

[0021] In the formula m k ,m k+1 They represent t respectively k ,t k+1 The gyroscope error at time t, where α is the autoregressive coefficient and ξ is the gyroscope error at time t. k+1 For t k+1 White noise at any given moment.

[0022] Preferably, step 2.3 specifically includes:

[0023] After completing t k After a filter update at time k+1, calculate the information sequence ε(ε1, ε2, ..., ε) for the previous k+1 times. k+1 variance Where ε k+1 For t k+1 Information at time, calculate ||ε k+1 |-|E[ε]|| and together with In comparison, where E[·] represents the mean of the matrix, β is a constant, and β∈(0,1);

[0024] if This data is considered interference data, and t is corrected. k+1 Kalman filter gain K at time step k+1 , αK k+1 Assigned to K k+1 α is a constant, α∈(0,1), recalculate t k+1 System status updates in real time Covariance Update P k+1|k+1 and filter output Enter the next filtering loop;

[0025] if Then modify t k One-step prediction of the system state vector at time step Will Assign to And recalculate t k+1 Information at time ε k+1 Then recalculate t k+1 System status updates in real time Covariance Update P k+1|k+1 and filter output Enter the next filtering loop;

[0026] if Then the filtering process will begin in the next time step.

[0027] Furthermore, the formula for calculating the peak value ratio of the high-frequency components in step 3.2 is as follows:

[0028]

[0029] In the formula, j is the number of decomposition layers, and N j Let W be the number of high-frequency components in the j-th layer, and max(|W j |) represents the maximum absolute value of the high-frequency component located in the j-th layer, |W j，i | represents the absolute value of the i-th high-frequency component located in the j-th layer, J j is the peak ratio of the j-th layer.

[0030] When J j ≤ρ and J j+1 When the value is greater than ρ, the number of decomposition layers L = j is determined, where ρ is a threshold for distinguishing whether a high-frequency component in a certain layer contains noise information, and ρ ∈ (0, 1).

[0031] In step 3.3, thresholding is performed on the low-frequency signal of layer L and the high-frequency signal of each layer using the Visushrink threshold. in |c1| represents the absolute value of the wavelet coefficients in the first layer of wavelet decomposition, M is the signal length, and a soft thresholding function is used.

[0032]

[0033] In the formula, r represents each wavelet coefficient.

[0034] Compared with the prior art, the beneficial effects of the present invention include:

[0035] 1) This invention uses an adaptive anti-outlier Kalman filter, which, compared to the classic Kalman filter, adds an anti-outlier step to avoid the influence of outliers in the signal on the filtering and improves the filtering accuracy;

[0036] 2) This invention performs thresholding on low-frequency signals during wavelet denoising while ensuring no signal distortion. Compared to traditional wavelet denoising, it better suppresses low-frequency noise.

[0037] 3) After adaptive anti-outlier Kalman filtering, the present invention performs wavelet processing on the signal. Compared with the previous Kalman filtering and wavelet fusion method, the algorithm is simpler and more effective.

[0038] 4) This invention provides an adaptive anti-outlier denoising scheme that integrates Kalman filtering and wavelet filtering for MEMS gyroscope signals, which can more effectively improve the accuracy of the sensor and reduce errors. Attached Figure Description

[0039] The present invention will be further described below with reference to the accompanying drawings and embodiments.

[0040] Figure 1 This is a flowchart illustrating the method of the present invention.

[0041] Figure 2 The following is an analysis diagram of the gyroscope x-axis noise reduction effect in the example.

[0042] Figure 3 The diagram shows the noise reduction effect of the gyroscope's y-axis in the example.

[0043] Figure 4The diagram shows the noise reduction effect of the gyroscope z-axis in the example. Detailed Implementation

[0044] like Figure 1 As shown, the gyroscope signal denoising method using Kalman filtering and Visushrink thresholding specifically includes the following steps:

[0045] Step 1: Establish a time-series ARMA model of the MEMS gyroscope signal:

[0046] The trend and periodic terms are extracted from the original signal of the MEMS gyroscope, and a time series ARMA model is established from its residuals. The final prediction error FPE criterion is used to determine the order of the model, and the random drift error model of the MEMS gyroscope is obtained as AR(1) model, which is expressed as:

[0047] m k+1 =am k +ξ k+1

[0048] In the formula m k m k+1 They represent t respectively k Time and t k+1 The gyroscope error at time t, where α is the autoregressive coefficient and ξ is the gyroscope error at time t. k+1 For t k+1 Time-based white noise;

[0049] Step 2: Use Kalman filtering to denoise the gyroscope signal;

[0050] Step 2.1: Based on the AR(1) model obtained in Step 1, establish the state equation and measurement equation of the Kalman filter, specifically as follows:

[0051]

[0052] Where A = a, B = 1, H = 1, V k+1 Let R and X be the estimation error of the AR(1) model. k+1 X k t k+1 t k The state quantity at time Z k+1 For t k+1 Measurement of time, η k The system process noise has a variance of Q, and its value is ξ as described in step 1. k+1 The variance of white noise, Q and R are uncorrelated, and the gyroscope error m is used as the filter input;

[0053] Step 2.2: Adaptive anti-outlier Kalman filter initialization. At the start time t0 of filtering, the filter's state vector is initialized using the gyroscope signal error m estimated by the ARMA model. The sum and variance matrix P0 are as follows:

[0054]

[0055]

[0056] Where E[·] represents the mean, (·) T This is the transpose of the matrix;

[0057] Step 2.3: Perform adaptive anti-outlier Kalman filtering;

[0058] Step 3: Determine the wavelet basis db2;

[0059] Determine the wavelet decomposition level: Perform wavelet decomposition on the filtered gyroscope signal using the selected db2 wavelet, and calculate the peak-to-peak ratio of the high-frequency components at each level, i.e.

[0060]

[0061] Where j is the number of decomposition levels, N j Let W be the number of high-frequency components in the j-th layer, and max(|W j |) represents the maximum absolute value of the high-frequency component located in the j-th layer, |W j,i | represents the absolute value of the i-th high-frequency component located in the j-th layer, J j The peak ratio of the j-th layer; when J j ≤ρ and J j+1 When the value is greater than ρ, the number of decomposition layers L = j is determined, where ρ is a threshold for distinguishing whether a high-frequency component in a certain layer contains noise information, and ρ ∈ (0, 1).

[0062] Thresholding is performed on the low-frequency signal of layer L and the high-frequency signal of each layer using the Visushrink threshold.

[0063]

[0064]

[0065] Where |c1| represents the absolute value of the wavelet coefficients in the first layer of wavelet decomposition, and M is the signal length. A soft thresholding function is used in this embodiment.

[0066]

[0067] Where r represents each wavelet coefficient.

[0068] Step 4: Perform wavelet reconstruction on the processed high-frequency and low-frequency signals using the waverec function in MATLAB software.

[0069] In step 2, the filtering process of the adaptive anti-outlier Kalman filter is as follows:

[0070] The first step is to predict the state in one step:

[0071]

[0072] The second step involves predicting the covariance matrix in one step:

[0073] P k+1|k =AP k|k A T +BQB T

[0074] The third step is to calculate the filter gain matrix:

[0075] K k+1 =P k+1|k H T HP k+1|k H T +R] -1

[0076] Step 4, calculate the information:

[0077]

[0078] Step 5, Status Update:

[0079]

[0080] Step 6, Covariance Update:

[0081] P k+1|k+1 =[I n -K k+1|k H]P k+1|k

[0082] Step 7, filter output at time k+1:

[0083]

[0084] In the formula For t k The one-step prediction of the system state vector at time t. t k , t k+1 The time update of the system state vector at time P. k+1|k For t k The one-step prediction of the system covariance matrix at time e k+1 For tk+1 Information about time, P k|k P k+1|k+1 t k , t k+1 System covariance update at time K k+1 For t k+1 Kalman filter gain at time I n Let be an n-order identity matrix, where n is the dimension of the filter state vector, [·] -1 It is the inverse of the matrix. This is the filter output at time k+1.

[0085] The adaptive anti-wildness value step: after completing t k After a filter update at time k+1, calculate the information sequence ε(ε1, ε2, ... ε) for the previous k+1 times. k+1 The variance S k+1 :

[0086]

[0087] And compare:

[0088]

[0089] In the formula, |·| represents the absolute value. If equation (1) is satisfied, then this data is considered to be interference data, and aK is... k+1 Assigned to K k+1 In the example where a∈(0,1), a = 0.5 is taken, and the state update is performed again in step 5.

[0090] If equation (1) is not satisfied, continue to compare whether it is satisfied:

[0091]

[0092] In the example, β = 0.7 is taken. If equation (2) is satisfied, then... Assign to g∈(0,1), in this example g=0.5, return to the fourth step to recalculate the information; if neither equation (1) nor equation (2) is satisfied, then enter the next time step of the filtering process loop.

[0093] This was verified through experiments, the details of which are as follows:

[0094] 1. The data in this experiment comes from the actual measurement data of a certain type of MEMS gyroscope static base. The static time is 1 hour and the sampling frequency is 200HZ.

[0095] 2. The autoregressive coefficients α of the gyroscope's x, y, and z axes are 0.0092, 0.0299, and 0.0080, respectively, and their corresponding white noise variances are 3.8595 × 10⁻⁶.-7 (rad 2 / s 2 ), 4.6831×10 -7 (rad 2 / s 2 ), 3.6487×10 -7 (rad 2 / s 2 ).

[0096] 3. The initial conditions for the Kalman filters of the gyroscope's x, y, and z axes are set as follows:

[0097]

[0098]

[0099]

[0100] in P represents the initial state vectors of the filters along the x, y, and z axes, respectively. 0x P 0y P 0z These are the initial covariances of the filters on the x, y, and z axes, respectively.

[0101] 4. Based on step 3, for this gyroscope signal, when the peak value ratio of the gyroscope error signal, J j ≤ρ,J j+1 When ρ > 0.05, determine the number of decomposition levels for the x, y, and z axis data as L. x =L y =L z =13.

[0102] 5. The x, y, and z axis data of the gyroscope were subjected to classical Kalman filtering, Visushrink denoising, and the proposed denoising method, respectively. The gyroscope noise figure was estimated by allan variance and the methods were compared.

[0103] The comparison of noise reduction effects of the gyroscope's x, y, and z axes through computer simulation is shown in the following figure. Figure 2 , Figure 3 , Figure 4 As shown, the noise figures of the gyroscope's x, y, and z axes are compared in Tables 1, 2, and 3, respectively.

[0104] Table 1 Comparison and Analysis of Gyroscope X-Axis Noise Figure

[0105]

[0106] Table 2 Comparison and Analysis of Gyroscope Y-axis Noise Figure

[0107]

[0108] Table 3 Comparison and Analysis of Gyroscope Z-Axis Noise Figure

[0109]

[0110] As can be seen from Table 1-3, the quantization noise figure, angle random walk, and zero-bias instability of the proposed method are all lower than those of Kalman filtering and wavelet denoising. In particular, the zero-bias instability of the gyroscope's x, y, and z axes is improved by 31.0%, 29.3%, and 30.5% respectively compared with Kalman filtering, and by 2.4%, 12.1%, and 12.4% respectively compared with traditional wavelet denoising.

[0111] Table 4 shows the variance comparison of the gyroscope's x, y, and z axes after denoising. As can be seen from Table 4, the variance after denoising using this scheme is smaller than that of Kalman filtering and wavelet denoising.

[0112] Table 4. Comparison and Analysis of Variance of Gyroscope x, y, and z Axes

[0113]

Claims

1. A method for gyro signal denoising by Kalman filtering and Visushrink thresholding, characterized in that, The method comprises the following steps: Step 1: establishing a time series ARMA model of the gyro signal; First, the trend term and the periodic term in the original gyro signal are extracted, and then a time series ARMA model is established for the residual error, the order of the model is determined by using the final prediction error (FPE) criterion, and the random drift error of the gyro is obtained; Step 2: denoising the gyro signal by using an adaptive anti-outlier Kalman filter; Step 2.1: establishing the state equation and the measurement equation of the Kalman filter according to the gyro random drift error model in step 1; Step 2.2: initializing the adaptive anti-outlier Kalman filter, and initializing the state vector and the covariance matrix of the filter by using the gyro error estimated by the gyro random drift error model at the start time of the filtering; Step 2.3: performing the adaptive anti-outlier Kalman filtering; After a one-step filter update at time t k , compute the variance of the information sequence , , at time , where is the information at time and compare it to , , where denotes the mean of the matrix, β is a constant, ;​​ If , the data is considered to be interference data, and the Kalman filter gain at time is modified to , is a constant, , the system state update and the covariance update at time are recalculated, and the filter output is entered into the next time filtering cycle; if Then correct One-step prediction of the system state vector at time step ,Will Assign to , And recalculate Information of time Then recalculate System status updates in real time Covariance Update and filter output Then proceed to the next filtering cycle; If then the filtering process cycle is entered for the next time instant; Step 3: performing threshold processing on the low-frequency component and the high-frequency component of the gyro signal filtered by the Kalman filter by using wavelet analysis; Step 3.1: selecting a wavelet function; Step 3.2: determining the wavelet decomposition level, and performing wavelet decomposition on the gyro signal output after filtering by using the selected wavelet function, and calculating the peak ratio of each layer of high-frequency component; The calculation formula of the peak ratio of the high-frequency component is: ； In the formula The number of decomposition layers, For the first The number of high-frequency components in the layer. For the position located at the The maximum absolute value of the high-frequency components of the layer. For the position located at the Layer The absolute value of each high-frequency component, For the first Peak ratio of the layer; When and , the number of decomposition layers is determined , wherein is a threshold value for distinguishing whether the high-frequency component of a certain layer contains noise information, ; Step 3.3: performing Visushrink threshold processing on the low-frequency signal and the high-frequency signal; right The low-frequency signals of each layer and the high-frequency signals of each layer are thresholded using the Visushrink threshold. ,in The absolute values ​​of the wavelet coefficients in the first layer of wavelet decomposition are given. M The length of the signal is determined by a soft threshold function. ； In the formula r is each wavelet coefficient; Step 4: performing wavelet reconstruction on the high-frequency and low-frequency signals of the gyro signal after threshold processing.

Images