CEEMDAN-WPT combined filtering algorithm for weakening GNSS multi-path error

Abstract

The invention provides a CEEMDAN-WPT combined filtering algorithm for weakening GNSS multi-path errors, which comprises the steps of reading and resolving GNSS observation data, obtaining an original coordinate sequence, decomposing the original sequence into a plurality of intrinsic mode functions by using fully adaptive noise integrated empirical mode decomposition, and filtering the intrinsic mode functions by using a CEEMDAN-WPT algorithm. Dividing the signal into a useful signal part, a mixed signal part and a useless signal part according to a specific division criterion, and directly eliminating the useless part. For the mixed signal part, wavelet packet decomposition is adopted to carry out filtering processing, effective signals in the mixed signal part are extracted, reconstruction is carried out on the effective signals and the useful signal part, and therefore multi-path errors in the original coordinate sequence are extracted. And then, calculating the correlation between the coordinate sequences of the first day and the second day and the time offset thereof, and subtracting the extracted multi-path error from the coordinate sequence of the second day according to the offset to obtain a corrected coordinate sequence. And finally, the model is applied to coal mine structure deformation monitoring.

Description

Technical Field

[0001] The present invention relates to the field of deformation monitoring of coal mine structures, especially the technology for weakening multipath errors. Specifically, a CEEMDAN-WPT combined filtering algorithm for weakening GNSS multipath errors is proposed to improve the deformation monitoring accuracy of GNSS in complex environments. Background Technique

[0002] The GNSS deformation monitoring technology has the characteristics of high precision, all-weather, no need for line-of-sight, and high degree of automation. Compared with traditional monitoring means, such as ground survey and remote sensing technology, the GNSS technology can make up for the deficiencies in accuracy and timeliness in the field of deformation monitoring. Through the high-precision position information obtained by GNSS, not only the surface movement can be monitored in real time, but also the real-time deformation monitoring of targets such as buildings and bridges can be carried out. In GNSS deformation monitoring, it is usually in the short-baseline static mode. Through the carrier phase differential technology, most errors can be eliminated, such as atmospheric errors and satellite clock errors. However, since the multipath effect is caused by the reflection or scattering of the surrounding environment of the measuring station, it cannot be eliminated or weakened by the differential technology. Therefore, the multipath error becomes one of the main factors affecting the GNSS positioning accuracy, seriously affecting the accuracy and reliability of the monitoring data.

[0003] The multipath effect has a certain spatio-temporal repeatability. Especially in a relatively static measuring station environment, since the satellite orbit is repetitive and the surrounding environment of the measuring station remains basically unchanged, the multipath effect will show periodic and repetitive characteristics. Based on this characteristic, by using the daily repeatability of the multipath, the multipath error in the GNSS monitoring data can be extracted through sidereal day filtering and modeled for error correction.

[0004] In terms of multipath error processing, the selection of algorithms is crucial. Commonly used existing methods mainly utilize wavelet decomposition (Wavelet) and empirical mode decomposition (EMD), etc. Among them, complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) is an improved EMD algorithm. By gradually adding noise to the original signal and performing step-by-step averaging processing, it can effectively improve the problem of mode aliasing in EMD decomposition, effectively reduce the impact of white noise on the decomposition result, and adaptively decompose the signal into multiple intrinsic mode functions (IMFs) and one residue. It can quickly and effectively decompose and extract the effective part of the signal, and has unique advantages in the extraction of multipath error signals. However, it is easily affected by data mutations and outliers, resulting in unstable decomposition results; wavelet packet decomposition (Wavelet Packet) is a signal decomposition method based on wavelet transform. Compared with wavelet transform, wavelet packet decomposition can perform more detailed frequency band division and decomposition on the signal, which is suitable for the processing of multi-band signals. However, its effect is affected by the selection of parameters such as wavelet basis functions, decomposition levels, and threshold functions, and there is currently no unified theory to guide how to select the optimal parameters.

[0005] Since multipath error is the main error source in GNSS deformation monitoring, there is an urgent need for a method that can quickly and effectively process multipath error to achieve high-precision deformation monitoring of GNSS in complex environments. Summary of the Invention

[0006] The present invention aims to solve the multipath error problem mentioned in the above background technology, and proposes a multipath error weakening method based on the CEEMDAN-WPT combined filtering algorithm - the CEEMDAN-WPT combined filtering algorithm for weakening GNSS multipath error, aiming to improve the positioning accuracy of GNSS deformation monitoring. This algorithm combines the advantages of CEEMDAN in extracting the effective components of the signal and the fine processing ability of wavelet packet decomposition in multi-scale frequency band decomposition, and effectively weakens the influence of GNSS multipath effect on the monitoring accuracy.

[0007] The present invention is realized through the following technical solutions:

[0008] A CEEMDAN-WPT combined filtering algorithm for weakening GNSS multipath error, the method includes the following steps:

[0009] Step 1: Read the original observation data from the GNSS receiver and perform preprocessing. This includes the solution of carrier phase differential technology to fix the ambiguity and obtain the coordinate sequence of each epoch. Finally, obtain the original coordinate sequence in the ENU direction during the ambiguity-fixed period;

[0010] Step 2: Use the CEEMDAN algorithm to decompose the original coordinate sequence. Through multiple iterative decompositions, the signal is decomposed into a series of Intrinsic Mode Function (IMF) terms and a residue.

[0011] Step 3: Classify the IMFs of the signal after CEEMDAN decomposition based on the division criteria of Permutation Entropy (PE) and energy value (E-value), and divide them into three parts: useful signals, mixed signals, and useless signals.

[0012] Step 4: Directly eliminate the classified useless signal part. For the mixed signal part, perform wavelet packet filtering to further decompose the signal into high-frequency and low-frequency parts, and process it according to certain threshold criteria to obtain the filtered signal.

[0013] Step 5: Reconstruct the signal obtained after wavelet packet decomposition and combine it with the useful signal part identified after CEEMDAN decomposition to obtain a multipath error model for subsequent error correction.

[0014] Step 6: Calculate the correlation between the coordinate sequences of adjacent two days, determine the time offset of the coordinate sequences between the two observation days, and determine the time offset at the maximum correlation for aligning the multipath errors in the current and reference coordinate sequences.

[0015] Step 7: Adjust the current coordinate sequence according to the calculated offset time, perform multipath error correction, and obtain the coordinate sequence after error correction.

[0016] Further, in Step 2, the CEEMDAN decomposition process is as follows:

[0017] (2.1) Add white noise w(t) with a mean of 0 and a variance of 1 and an adaptive white noise coefficient β to the original coordinate sequence x(t).

[0018] (2.2) Perform N times of EMD decomposition on the sequence x(t) + βw(t) after adding white noise, and obtain the average value of N first-order IMF components to get the first-order IMF component.

[0019] (2.3) Use the sequence x(t) to subtract the first-order IMF component IMF1(t) to obtain the first residual signal component r1(t).

[0020] (2.4) Add white noise again to the first residual component and calculate the second-order IMF component.

[0021] (2.5) Repeat steps (2.3) and (2.4) until no further decomposition is possible. After completing the entire decomposition step, the original sequence can be expressed as:

[0022]

[0023] Wherein, IMF j (t) is the j-th IMF component obtained by decomposition, and R(t) is the residue; j = 1, 2, …, N, where N is the total number of IMF components.

[0024] Further, in the third step, the division criteria are as follows:

[0025] (3.1) Calculate the permutation entropy PE of each IMF component, and perform normalization processing on it. Determine that the IMF components with entropy values greater than 0.65 are the useless signal parts, and set the criteria for useless signals by setting the demarcation point k1.

[0026] (3.2) The more useful signal components in the IMF component, the higher the energy value. Use the following formula to calculate the energy value E of the IMF component, and determine that the IMF components with energy values greater than 1 are the useful signal parts, and set the demarcation point k2.

[0027]

[0028] Where pe j is the PE value of the j-th IMF component. When the noise content of the first few high-frequency noise IMF components is large and the calculated PE values are similar, when E j is greater than 1, it indicates that the PE value of the j-th IMF component is smaller than half of the average PE value of the previous j - 1 IMF components, proving that starting from the j-th IMF component, the effective components dominate the signal.

[0029] (3.3) Eliminate the IMF components with PE values greater than k1 as useless signals, and divide the IMF components between k1 and k2 into the mixed signal part, indicating that this signal part contains both noise and useful signals and needs further processing. Determine the IMF components greater than the demarcation point k2 as the useful signal part.

[0030] Further, in the fourth step, the steps of wavelet packet filtering are as follows:

[0031] (4.1) Perform wavelet packet decomposition on the mixed signal part obtained from the third step to decompose the signal into multiple scales and multiple frequency bands.

[0032] (4.2) Estimate the noise standard deviation σ in the signal through the median absolute deviation (σ = median(abs(ω j,k )) / 0.6745, where ω j,k is the coefficient of the j-th layer and the k-th node of the wavelet packet), and determine the filtering threshold N is the signal length. This formula is usually called the "VisuShrink" method and is applicable to Gaussian white noise). According to the set threshold thr, the improved threshold function is used to process the part greater than the threshold, and the part less than the threshold is not processed. The expression of this improved threshold function is as follows. This threshold function combines the advantages of hard and soft thresholds, can more effectively remove noise, and retain the effective part of the signal:

[0033]

[0034] In the formula, ω represents the wavelet packet coefficient, and λ represents the threshold parameter. When |ω|≥λ, threshold processing is performed. While removing the noise components less than the threshold, the structure of the larger coefficients is retained; when |ω|≤λ, η(ω,λ) = 0, that is, the coefficients less than the threshold are directly set to zero, achieving noise suppression.

[0035] (4.3) Reconstruct the signal after wavelet packet decomposition and threshold processing to obtain the filtered mixed signal. This reconstructed signal is used as part of the multipath error correction model to provide a basis for subsequent error correction of the GNSS coordinate sequence.

[0036] Furthermore, in the fifth step, the signal reconstruction steps are as follows:

[0037] (5.1) The wavelet decomposition of the input signal x of length N j,k (m) is:

[0038]

[0039] In the formula, x j,k (n) is the wavelet packet coefficient series of the k (k = 0, 1, 2, …, 2 j -1)th sub-band in the jth layer of decomposition; h and g are the low-pass filter and high-pass filter respectively, and satisfy g(n) = (-1) n h(-n).

[0040] (5.2) The wavelet packet reconstruction algorithm is:

[0041]

[0042] Furthermore, in the sixth step, the correlation calculation is as follows:

[0043]

[0044] Use the formula to calculate the correlation r of the GNSS coordinate sequences in different time periods. Among them, x i and y i are respectively the ith values in the two time series; and They are the averages of two time series respectively.

[0045] By calculating the correlation between the coordinate sequences observed on two adjacent days, through cross-correlation analysis, the time offset corresponding to the maximum correlation is identified. According to the formula, the correlation coefficients at different time offsets are calculated, and the offset value with the maximum correlation coefficient is obtained. This time offset is the optimal alignment point between the two-day coordinate sequences.

[0046] Furthermore, in step seven, the multi-path error correction steps are as follows:

[0047] Adjust the coordinate sequence monitored on the second day according to the calculated offset time. For example, if the calculated offset time is 5 seconds, the current coordinate sequence is shifted forward or backward by 5 seconds as a whole to correct the multi-path error in the two sets of data. Based on the adjusted coordinate sequence, the data corresponding to the offset time in the multi-path error model is subtracted from the adjusted coordinate sequence point by point to eliminate the multi-path error. Finally, the coordinate sequence after error correction is obtained.

[0048] The beneficial effects brought by adopting the above technical solutions are as follows:

[0049] The present invention provides a method for processing multi-path errors based on CEEMDAN-WPT combined filtering. This method fully considers the deficiencies of the wavelet packet and CEEMDAN methods and optimizes their combination, strengthening the reliability of the wavelet packet and CEEMDAN in processing multi-path errors, and can more quickly and effectively extract and correct multi-path errors, playing a good role in practical engineering applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0050] Figure 1 is a flowchart of the CEEMDAN-WPT combined filtering algorithm for weakening GNSS multi-path errors;

[0051] Figure 2 is a coordinate sequence diagram of the original ENU direction obtained by solving for two consecutive days;

[0052] Figure 3 is the CEEMDAD decomposition diagram of the east direction on DOY202;

[0053] Figure 4 is the multi-path error diagram extracted on DOY202;

[0054] Figure 5 is the coordinate sequence diagram after multi-path error correction on DOY203. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0055] Hereinafter, the technical solutions of the present invention will be described in detail with reference to the accompanying drawings.

[0056] The present invention proposes a CEEMDAN-WPT combined filtering algorithm for weakening GNSS multipath errors, and the implementation flowchart is as Figure 1 shown.

[0057] The data of this implementation case is from the short-baseline GNSS observation data of a coal mine shaft tower, the baseline length is about 400m, the observation epoch interval is 5s, the cut-off elevation angle is 10°, and the original coordinate sequence obtained by carrier phase difference resolution mainly contains multipath errors and random noise.

[0058] Step 1: Calculate the ENU original coordinate sequences for two consecutive days (DOY202 - DOY203), and the obtained original coordinate sequences are as Figure 2 shown. This sequence contains a large amount of multipath errors and noise signals.

[0059] Step 2: Use CEEMDAN to decompose the original coordinate sequence, and obtain 11 IMF terms and 1 residue. Figure 3 is the CEEMDAN decomposition result in the east direction of DOY202. Each IMF component corresponds to signal and noise components of different frequencies.

[0060] Step 3: According to the IMF results obtained by CEEMDAN decomposition, calculate the permutation entropy PE and the energy value E, and determine the standard demarcation points k1 and k2 for signal division. Table 1 shows the calculation results of the PE values and E values of the IMF components in the east direction of DOY202.

[0061] Table 1 PE values and E values corresponding to each IMF component in the east direction of DOY202

[0062] IMF component 1 2 3 4 5 6 Permutation entropy PE 0.94 0.78 0.57 0.42 0.31 0.24 Energy value E - 0.20 0.49 0.81 1.17 1.51 IMF component 7 8 9 10 11 12 Permutation entropy PE 0.19 0.17 0.16 0.15 0.14 0.12 Energy value E 1.76 1.86 1.82 1.75 1.75 1.57

[0063] According to the permutation entropy in the table, it can be seen that the permutation entropy value of the 3rd-order IMF component is less than 0.65 for the first time, so k1 = 3; the E value of the 5th-order IMF component is greater than 1, so k2 = 5 is determined. Therefore, IMF1 - 2 are classified as the useless signal part, IMF3 - 5 are classified as the mixed signal part, and IMF6 - 11 are the useful signal part.

[0064] Step 4: After removing the useless signal part, perform wavelet packet filtering on the mixed signal part (IMF3 - 5) to remove the noise part in the mixed signal. The mixed signal after wavelet packet filtering is reconstructed with the useful signal part (IMF6 - 11) to obtain the multipath error. The extracted multipath error results are as Figure 4 shown.

[0065] Step 5. Calculate the correlation between the two-day coordinate sequences and the maximum correlation epoch interval. The results are shown in Table 2. Since the receiver sampling interval is 5 seconds, the calculated maximum correlation epoch intervals are all multiples of 5 seconds, with a few seconds deviation from the actual interval.

[0066] Table 2 Correlation between the coordinate sequences of DOY202 - DOY203 and the maximum correlation epoch interval table

[0067]

[0068] Step 6. Perform time offset on the data of the second day according to the maximum correlation epoch time, and subtract the multipath error extracted on the first day to finally correct the multipath error. Figure 5 Table 3 shows the comparison results before and after the multipath error correction. The RMS and improvement rate in the table are used to represent the improvement degree of the multipath error. The higher the improvement rate value, the more obvious the correction effect. The improvement rate calculation formula is:

[0069]

[0070] Table 3 RMS values and improvement rate table before and after the multipath error correction on DOY203

[0071]

Claims

1. A CEEMDAN-WPT combined filtering algorithm for reducing GNSS multipath error, the method comprising the following steps: Step 1: Read the original observation data from the GNSS receiver and perform preprocessing, including solving the carrier phase difference technology, to fix the ambiguity, obtain the coordinate sequence of each epoch, and finally obtain the original coordinate sequence of the ENU direction in the ambiguity fixed period; Step 2: Use the CEEMDAN algorithm to decompose the original coordinate sequence into a series of intrinsic mode function terms IMF and a remainder term through multiple iterative decompositions; Step 3: Classify the signal IMF after CEEMDAN decomposition based on the division criteria of permutation entropy and energy value, and divide it into three parts: useful signal, mixed signal and useless signal; Step 4: directly remove the useless signal part after classification, and perform wavelet packet filtering on the mixed signal part, further decompose the signal into high-frequency and low-frequency parts, and process it according to a certain threshold criterion to obtain the filtered signal; Step 5: Reconstruct the signal obtained after wavelet packet decomposition and combine it with the useful signal part identified after CEEMDAN decomposition to obtain the multipath error model for subsequent error correction. Step 6: Calculate the correlation between the coordinate sequences of two adjacent days, determine the time offset of the coordinate sequence between the two observation days, and determine the time offset at the maximum correlation to align the multipath errors in the current and reference coordinate sequences. Step 7: Adjust the current coordinate sequence according to the calculated offset time, perform multipath error correction, and obtain a coordinate sequence after error correction.

2. The CEEMDAN-WPT combined filtering algorithm for reducing GNSS multipath error according to claim 1, characterized in that: In step 2, the CEEMDAN decomposition process steps are as follows: (2.1) Add white noise w(t) with mean 0 and variance 1 to the original coordinate sequence x(t), and the adaptive white noise coefficient β, (2.2) Perform N EMD decompositions on the sequence x(t)+βw(t) after adding white noise, obtain the average value of N first-order IMF components, and obtain the first-order IMF component. (2.3) Use the sequence x(t) to subtract the first-order IMF component IMF1(t) to obtain the first residual signal component r1(t), (2.4) Add white noise to the first residual component again and calculate the second-order IMF component, (2.5) Repeat steps (2.3) and (2.4) until no further decomposition is possible. After completing the entire decomposition step, the original sequence can be expressed as: In the formula, IMF j (t) is the jth IMF component obtained by decomposition, R(t) is the remainder; j = 1, 2, ..., N, N is the total number of IMF components.

3. The CEEMDAN-WPT combined filtering algorithm for reducing GNSS multipath error according to claim 1, characterized in that: In step 3, the classification criteria are as follows: (3.1) Calculate the permutation entropy PE of each IMF component and normalize it. Determine the IMF component with an entropy value greater than 0.65 as the useless signal part. By setting the demarcation point k1, the standard of the useless signal is defined. (3.2) The more useful signal components there are in the IMF component, the higher the energy value. The energy value E of the IMF component is calculated using the following formula. The IMF component with an energy value greater than 1 is determined to be the useful signal part, and the demarcation point k2 is set. Where pe j is the PE value of the jth IMF component. When the noise content of the current several high-frequency noise IMF components is large and the calculated PE values ​​are similar, when E j When it is greater than 1, it indicates that the PE value of the jth IMF component is smaller than half of the average PE value of the previous j-1 IMF components, proving that starting from the jth IMF component, the effective component dominates the signal. (3.3) The IMF components with PE values ​​greater than k1 are eliminated as useless signals, and the IMF components between k1 and k2 are classified as mixed signal parts, indicating that the signal part contains both noise and useful signals and needs further processing. The IMF components greater than the dividing point k2 are determined as useful signal parts.

4. The CEEMDAN-WPT combined filtering algorithm for reducing GNSS multipath error according to claim 1, characterized in that: In step 4, the steps of wavelet packet filtering are as follows: (4.1) Perform wavelet packet decomposition on the mixed signal obtained in step 3, and decompose the signal into multi-scale and multi-band. (4.2) The noise standard deviation σ in the signal is estimated by the median absolute deviation to determine the filtering threshold N is the signal length. According to the set threshold thr, the improved threshold function is used to process the part greater than the threshold, and the part less than the threshold is not processed. The expression of the improved threshold function is as follows. This threshold function combines the advantages of soft and hard thresholds and can remove noise more effectively while retaining the valid part of the signal: Where ω represents the wavelet packet coefficient, and λ represents the threshold parameter. When |ω|≥λ, threshold processing is performed to remove noise components less than the threshold while retaining the structure of larger coefficients; when |ω|≤λ, η(ω,λ)=0, that is, the coefficients less than the threshold are directly set to zero, achieving noise suppression. (4.3) The signal after wavelet packet decomposition and threshold processing is reconstructed to obtain a filtered mixed signal. This reconstructed signal is used as part of the multipath error correction model to provide a basis for the error correction of the subsequent GNSS coordinate sequence.

5. The CEEMDAN-WPT combined filtering algorithm for reducing GNSS multipath error according to claim 1, characterized in that: In step 5, the steps of wavelet packet decomposition and reconstruction are as follows: (5.1) Input signal x of length N j,k The wavelet decomposition of (m) is: Where x j,k (n) is the kth (k=0,1,2,…,2 j -1) sub-band wavelet packet coefficient series; h and g are low-pass filter and high-pass filter respectively, and satisfy g(n)=(-1) n h(-n), (5.2) The wavelet packet reconstruction algorithm is: x j,k (m)=2[∑ n h(m-2n)x j+1,2k-1 (n)+∑ n g(m-2n)x j+1,2k (n)]。 6. The CEEMDAN-WPT combined filtering algorithm for reducing GNSS multipath error according to claim 1, characterized in that: In step 6, the correlation is calculated as follows: The correlation r of GNSS coordinate sequences in different time periods is calculated using the formula, where x i and i are the i-th values ​​in the two time series respectively; and are the average values ​​of the two time series, By calculating the correlation between the coordinate sequences observed on two adjacent days, the time offset corresponding to the maximum correlation is identified through cross-correlation analysis. The correlation coefficient under different time offsets is calculated according to the formula, and the offset value with the largest correlation coefficient is obtained. This time offset is the optimal alignment point between the coordinate sequences of the two days.

7. The CEEMDAN-WPT combined filtering algorithm for reducing GNSS multipath error according to claim 1, characterized in that: In step seven, the multipath error correction steps are as follows: Adjust the currently monitored coordinate sequence according to the calculated offset time. Based on the adjusted coordinate sequence, the data corresponding to the offset time in the multipath error model is subtracted point by point from the adjusted coordinate sequence to eliminate the multipath error. Finally, the error-corrected coordinate sequence is obtained.

Images