A combined filtering method for weakening GNSS multipath errors

Abstract

The present invention discloses a combined filtering method for weakening GNSS multipath errors, comprising the following steps: (1) reading the original coordinate sequence of GNSS observation data; (2) decomposing the original coordinate sequence into a series of intrinsic mode function terms and a remainder term by using ensemble empirical mode decomposition; (3) subdividing the IMF terms and the remainder term into noise terms, transition terms and useful terms by using classification indices k1 and k2; (4) discarding the noise terms and performing wavelet filtering on the transition terms; (5) reconstructing the filtered transition terms of the useful terms to obtain the multipath error sequence of the current day; (6) processing the original coordinate sequence according to the above steps, taking the multipath error sequence extracted on the first day as a reference signal, and calculating the amplitude ratio coefficient between the reference signal sequence and the multipath error sequences of the remaining days; (7) restoring the amplitude of the multipath error sequence and using it as an error correction model, subtracting it from the original coordinate sequence to obtain the corrected coordinate sequence; (8) outputting the coordinate sequence after multipath error correction.

Description

Technical Field

[0001] The present invention belongs to the field of surveying and mapping, and relates to a combined filtering method for weakening GNSS multipath errors. Background Art

[0002] The GNSS deformation monitoring technology has the advantages of high sampling frequency, real-time high precision, no need for line of sight, full automation, all-weather, etc., and can be applied to the real-time monitoring of civil engineering structures and natural disasters. GNSS deformation monitoring is usually carried out in the form of short-baseline static measurement. The carrier phase differential technology can basically eliminate the errors with strong correlations such as satellite and receiver clock errors, ephemeris errors, and atmospheric delays. However, the deformation monitoring environment is often relatively complex. The multipath effect caused by reflectors near the measuring station is not only unavoidable, but also difficult to weaken by differential, directly leading to a decrease in GNSS positioning accuracy. In addition, the noise generated by the measuring station receiver itself has a wide frequency domain distribution range and obvious Gaussian white noise characteristics, which is mixed with the multipath error, further increasing the difficulty of multipath error processing.

[0003] When the environment around the GNSS measuring station changes little or remains unchanged, the multipath effect shows strong weekly repeatability. Based on this characteristic, using the sidereal day filter to extract multipath errors and model repeatability from GNSS monitoring data for consecutive days can effectively weaken the multipath errors, with fast processing efficiency and high cost performance. Among them, the extraction of multipath errors is very critical in the whole processing process. The main existing extraction methods include the wavelet analysis (Wavelet) series, empirical mode decomposition (EMD) series, etc. Wavelet analysis (Wavelet) shows obvious advantages in multipath error processing due to its multi-resolution analysis and simultaneous localization in time and frequency domains. However, its effect is limited by the determination of parameters such as wavelet basis, decomposition level, threshold, and its estimation method, and there is no unified supporting theory for parameter selection; Ensemble Empirical Mode Decomposition (EEMD) is an improved method in the EMD series. It effectively improves the mode mixing problem by adding Gaussian white noise, performing multiple EMDs, and mean processing, and inherits the characteristic of being able to adaptively decompose the signal into several Intrinsic Mode Function (IMF) terms and one residue term, and can effectively and quickly extract the useful components in the signal, showing unique advantages in multipath error processing. However, its effect is limited by the amount of white noise added. Too much or too little white noise will lead to unsatisfactory final results.

[0004] As the main error source of GNSS deformation monitoring, the multipath effect seriously affects the accuracy of deformation monitoring. With the popularization of GNSS technology in deformation monitoring, there is an urgent need for a multipath error processing method that is convenient for computer operation, has a rigorous theory, and high cost performance to achieve real-time high-precision GNSS monitoring. Summary of the Invention

[0005] To solve the technical problems mentioned in the above background art, the present invention provides a combined filtering method for weakening GNSS multipath errors, which can more effectively weaken multipath errors in view of the deficiencies in the multipath error processing of Wavelet and EEMD.

[0006] To achieve the above technical objectives, the technical solution of the present invention is as follows:

[0007] A combined filtering method for weakening GNSS multipath errors, comprising the following steps:

[0008] Step 1: Read the original coordinate sequence of GNSS observation data, and start the calculation after reading.

[0009] Step 2: Use Ensemble Empirical Mode Decomposition (EEMD) to decompose the original coordinate sequence into a series of Intrinsic Mode Function (IMF) terms and a residue.

[0010] Step 3: Use classification indicators k 1 and k 2 to subdivide the IMF terms and the residue into noise terms, transition terms, and useful terms.

[0011] Step 4: Discard the noise terms, and perform Wavelet filtering on the transition terms to obtain the filtered transition terms.

[0012] Step 5: Reconstruct the useful terms and the filtered transition terms to obtain the multipath error sequence.

[0013] Step 6: Take the multipath error sequence of the first day as the reference signal sequence, and calculate the amplitude ratio coefficient between the reference signal sequence and the multipath error sequences of the remaining days.

[0014] Step 7: Restore the amplitudes of the multipath error sequences of the remaining days, and use them as the error correction model to subtract from the corresponding original coordinate sequences to obtain the corrected coordinate sequences.

[0015] Step 8: End the calculation and output the coordinate sequences after multipath error correction.

[0016] Preferably, in the above Step 2, the EEMD processing procedure is as follows:

[0017] 21) Add Gaussian white noise n(t) to the original coordinate sequence y(t) to obtain the coordinate sequence to be processed y m (t). The magnitude of the added Gaussian white noise n(t) is determined by the standard deviation ratio σ between it and the original coordinate sequence y(t).

[0018] 22) Perform EMD processing on y m (t) to obtain a set of IMF terms and a residue. The EMD processing steps are as follows:

[0019] (221) Let y m (t) = R i-1 (t), i = 1;

[0020] (222) Find all the extreme points on R i-1 (t). Use cubic spline interpolation to fit the sequences of minimum and maximum points respectively, and obtain the upper and lower envelopes of R i-1 (t);

[0021] (223) Calculate the mean sequence m i (t) of the upper and lower envelopes, and calculate the difference h i-1 (t) between it and R i (t) = R i-1 (t) - m i (t);

[0022] (224) Judge whether h i (t) meets the following conditions: ① The number of extreme points and zeros of h i (t) is equal or differs by only 1; ② The upper and lower envelopes of h i (t) are locally symmetric about the time axis. If it meets the conditions, h i (t) = IMF i (t), otherwise, let R i-1 (t) = h i (t), and jump to step (222) until the conditions are met;

[0023] (225) Calculate R i (t) = R i-1 (t) - IMF i (t), and judge whether R i (t) is monotonic or the number of its extreme points is less than 2. If R i (t) is monotonic or the number of its extreme points is less than 2, the decomposition ends, let n = i, otherwise, let i = i + 1, jump to step (222), until the conditions are met, and obtain a decomposition result consisting of n IMF terms and 1 residue term. The expression is as follows:

[0024]

[0025] In the formula, is the i-th IMF term obtained by EMD, is the residue; i = 1, 2,..., n, and n is the total number of IMF terms;

[0026] 23) Repeat steps 21) and 22) N times to obtain N groups of IMF terms and residues;

[0027] 24) Calculate the mean values of N groups of IMF terms and the residue term respectively. The n mean values of the IMF terms and 1 mean value of the residue term obtained are used as the final results, that is, the n mean values of the IMF terms and 1 mean value of the residue term obtained by EEMD are as shown in the following formula:

[0028]

[0029] In the formula, j = 1, 2, …, N, where N is the total number of EMD repetitions; IMF i (t) is the mean value of the i-th IMF term, and R n (t) is the mean value of the residue term.

[0030] Preferably, in the third step, the determination step of the classification index k 1 is as follows:

[0031] 31) Through EEMD processing, the coordinate sequence reconstructed starting from the mean value of the k-th IMF term is expressed as follows:

[0032]

[0033] 32) Calculate the squared Euclidean distance between two continuously reconstructed coordinate sequences using the continuous root mean square error criterion, and the calculation is as follows:

[0034]

[0035] In the formula, k = 1, 2, 3,...., n - 1; L is the length of the coordinate sequence, that is, the number of epochs;

[0036] 33) Calculate the CMSE value according to the value range of k. The k value at the first occurrence of the minimum value of the CMSE value is used as the starting point of the transition term, that is, k = k 1 :

[0037]

[0038] In the formula, argfirstlocalmin(·) is the k value corresponding to obtaining the first minimum value; is the largest integer not exceeding , and n is the total number of mean values of the IMF terms.

[0039] Preferably, in the third step, the determination step of the classification index k 2 is as follows:

[0040] 41) According to the EEMD decomposition result, calculate the product ET k of the energy density of the k-th IMF mean term and its average period, and the formula is as follows:

[0041]

[0042] Wherein, E k is the energy density of the k-th IMF mean term, is the average period of the k-th IMF mean term, and the calculation formulas are as follows:

[0043]

[0044]

[0045] Wherein, k = 1, 2, 3,...., n, O k is the number of extreme points of the k-th IMF term; L is the length of the k-th IMF term, that is, the number of epochs;

[0046] 42) Calculate the effective coefficient R k of the k-th IMF term, and the expression is as follows:

[0047]

[0048] Wherein, k = 2, 3,...., n;

[0049] 43) Calculate R k successively according to the value range of k, and take the k value when R k > 3 appears for the first time as the starting point of the useful term, that is, k = k 2 .

[0050] Preferably, in the fourth step, the processing steps of wavelet filtering are as follows:

[0051] 51) Perform wavelet decomposition on the transition term sequence to obtain a decomposition result composed of multiple high-frequency and low-frequency terms;

[0052] 52) Perform soft thresholding on the high-frequency term sequence, and do not process the low-frequency terms. The soft threshold function soft(·) is expressed as follows:

[0053]

[0054] Wherein, x t is the amplitude of the high-frequency term sequence, t = 1, 2, 3,...., L, L is the sequence length, that is, the number of epochs; sgn(·) is the step function, when x t > 0, the result is taken as 1, when x t < 0, the result is taken as -1, and when x t = 0, it is taken as 0; λ is the threshold, and the calculation is as follows:

[0055] 53) Reconstruct the low-frequency term and the high-frequency term after soft thresholding to obtain the filtered transition term.

[0056] In the fifth step, the reconstructed GNSS multipath error sequence is expressed as follows:

[0057]

[0058] In the formula, wavelet is wavelet filtering process; IMF i (t) is the i-th IMF term obtained by EEMD, and R n (t) is the remainder term; k 1 and k 2 are the indicators for distinguishing the noise term from the transitional term and the transitional term from the useful term respectively. Preferably, in the sixth step, the amplitude ratio coefficient a r is calculated as follows:

[0059]

[0060] In the formula, average(·) is the mean processing; r(t) is the reference signal sequence, is the reconstructed multipath error sequence.

[0061] Preferably, in the seventh step, the calculation of the multipath error amplitude recovery is as follows:

[0062]

[0063] In the formula, s(t) is the multipath error after amplitude recovery; a r is the amplitude ratio coefficient; is the reconstructed multipath error sequence.

[0064] Beneficial effects brought by adopting the above technical solutions:

[0065] The present invention fully considers the physical characteristics of the multipath effect in GNSS deformation monitoring and the deficiencies of the Wavelet and EEMD methods in multipath error processing, and introduces a combined classification criterion on the basis of the EEMD-Wavelet combined filtering. This method strengthens the reliability of EEMD and Wavelet in multipath error processing, can not only extract and correct the multipath error more effectively, but also has a simple process and a unified algorithm, and can adapt to actual engineering. Description of the Drawings

[0066] Figure 1 is the implementation flowchart of a combined filtering method for weakening GNSS multipath error;

[0067] Figure 2 is the coordinate sequence diagram in the elevation direction for 3 days;

[0068] Figure 3 is the EEMD decomposition result diagram;

[0069] Figure 4 is the waveform diagram of the filtered transitional term;

[0070] Figure 5 is the multipath error sequence diagram of the first day;

[0071] Figure 6 is the coordinate sequence diagram before and after multipath error correction. Specific implementation manner

[0072] The technical solution of the present invention will be described in detail below in conjunction with the accompanying drawings.

[0073] The present invention proposes a combined filtering method for weakening GNSS multipath error, and the implementation process is as Figure 1 shown.

[0074] The data of this implementation case is from the ultra-short baseline GNSS observation data of Curtin University. The baseline length is 6.15 m. The receiver used is Trimble-NETR9. The observation epoch interval is 30 s. The cut-off elevation angle is 10°. The original coordinate sequence obtained by resolving the carrier phase difference mainly contains multipath error and random noise.

[0075] Step 1: The computer reads the original coordinate sequence of 2500 epochs in the elevation direction at the same time period for three consecutive days, as Figure 2 shown, and the original coordinate sequence of the first day is taken to elaborate on the following steps.

[0076] Step 2: Use EEMD to decompose the original coordinate sequence with default parameters to obtain 10 IMF terms and 1 residue term, Figure 3 as shown in the EEMD decomposition result;

[0077] Step 3: According to the EEMD decomposition result, use the CMSE criterion and the effective coefficient R defined based on energy density and average period to calculate the classification indexes k 1 and k 2 , and the calculation results are shown in Table 1.

[0078] Table 1. CMSE values and energy effective coefficient R

[0079] Number of Imf terms k 1 2 3 4 5 6 7 8 9 10 CMSE value 11.38 3.75 2.99 3.39 1.85 2.29 1.37 0.45 1.28 0.03 R value - 0.32 0.57 1.73 0.28 0.90 3.46 0.95 12.52 5.39

[0080] In the table, the CMSE value of the 3rd-order IMF term first appears as a minimum value, k 1 = 3, the R value of the 6th-order IMF term is greater than 3, k 2 = 6, IMF1-2 terms are noise terms, IMF3-6 terms are transitional terms, and IMF7-10 and the residue term are useful terms.

[0081] Step 4: Eliminate the noise terms, perform wavelet filtering on the transition terms with default parameters to obtain the filtered transition terms; as Figure 4 shown.

[0082] Step 5: Reconstruct the filtered transition terms and the useful terms to obtain the multipath error;

[0083] Step 6: Use the multipath error on the first day as the reference signal, and then calculate the amplitude ratio coefficients a r , a r respectively with the multipath errors extracted in the following two days. The obtained values are as Figure 5 shown in Table 2;

[0084] Table 2. a r Value

[0085] Number DAY1 - DAY2 DAY1 - DAY3 <![CDATA[a r > 0.96 0.97

[0086] Step 7: Restore the amplitude of the multipath error, use this as the error correction model, and subtract it from the original coordinate sequence to achieve multipath error correction. Figure 6 Table 3 shows the corrected coordinate sequence and its RMSE respectively:

[0087] Table 3. RMSE of the coordinate sequence after multipath error correction

[0088]

[0089] The embodiments are only used to illustrate the technical idea of the present invention, and the protection scope of the present invention cannot be limited thereby. Any modification made on the basis of the technical solution according to the technical idea proposed by the present invention falls within the protection scope of the present invention.

Claims

1. A combined filtering method for weakening GNSS multipath errors, characterized in that, it includes the following steps: Step 1: Read the original coordinate sequence of GNSS observation data, and start the calculation after reading; Step 2: Use Ensemble Empirical Mode Decomposition (EEMD) to decompose the original coordinate sequence into a series of Intrinsic Mode Function (IMF) terms and a residue; Step 3. Use classification indicators k 1 and k 2 to subdivide the IMF terms and the residual term into noise terms, transitional terms, and useful terms; Step 4: Discard the noise terms, and perform wavelet filtering (Wavelet) on the transitional terms to obtain the filtered transitional terms; Step 5: Reconstruct the useful terms and the filtered transitional terms to obtain the multipath error sequence; Step 6: Take the multipath error sequence of the first day as the reference signal sequence, and calculate the amplitude ratio coefficient between the reference signal sequence and the multipath error sequences of the other days; Step 7: Restore the amplitudes of the multipath error sequences of the other days, and use this as the error correction model to subtract from the corresponding original coordinate sequences to obtain the corrected coordinate sequences; Step 8: End the calculation and output the coordinate sequences after multipath error correction; In the said Step 2, the EEMD processing procedure is as follows: 21) Add Gaussian white noise n(t) to the original coordinate sequence y(t) to obtain the coordinate sequence y m (t) to be processed. The magnitude of the added Gaussian white noise n(t) is determined by the standard deviation ratio σ between it and the original coordinate sequence y(t); 22) For y m (t) is processed by EMD to obtain a set of IMF terms and a residue. The EMD processing steps are as follows: (221) Let y m (t) = R i-1 (t), i = 1; (222) Search for R i-1 All the extreme points on (t), and use cubic spline interpolation to fit the sequences of minimum points and maximum points respectively to obtain R i-1 The upper and lower envelopes of (t); (223) Calculate the mean sequence m of the upper and lower envelope lines i (t), and take the difference h i-1 (t) with R i (t) = R i-1 (t) - m i (t); (224) Determine h i (t) satisfies the following: ① The number of extreme points and zeros of h i (t) is equal or differs by only 1; ② The upper and lower envelopes of h i (t) are locally symmetric about the time axis; if satisfied, h i (t) = IMF i (t), otherwise, let R i-1 (t) = h i (t), and jump to step (222) until the condition is satisfied; (225) Calculate R i R(t) = R i-1 (t) - IMF i (t), and determine whether R i (t) is monotonic or the number of its extreme points is less than 2. If R i (t) is monotonic or the number of its extreme points is less than 2, the decomposition ends, let n = i. Otherwise, let i = i + 1, jump to step (222) until the condition is satisfied, and obtain a decomposition result consisting of n IMF terms and 1 residual term. The expression is as follows: In the formula, is the i-th IMF item obtained from EMD, is the remainder item; i = 1, 2, …, n, where n is the total number of IMF items; 23) Repeat Steps 21) and 22) N times to obtain N groups of IMF terms and residues; 24) Calculate the mean values of the N groups of IMF terms and residues respectively, and take the mean values of the n IMF terms and the mean value of 1 residue obtained as the final results, that is, the mean values of the n IMF terms and the mean value of 1 residue obtained by EEMD, as shown in the following formula: where \(j = 1, 2, \ldots, N\), \(N\) is the total number of EMD repetitions; \(IMF_i(t)\) i is the mean value of the \(i\)-th IMF term, and \(R(t)\) n is the mean value of the remainder term; In step 3, the classification index k 1 is determined as follows: 31) The coordinate sequence reconstructed starting from the mean value of the k-th IMF term through EEMD processing is expressed as follows: 32) Calculate the squared Euclidean distance between two consecutive reconstructed coordinate sequences using the continuous root mean square error criterion, as follows: ​ In the formula, k = 1, 2, 3,...., n - 1; L is the length of the coordinate sequence, that is, the number of epochs; 33) Calculate the CMSE value according to the value range of k. The k value at the first occurrence of the minimum value of the CMSE value is used as the starting point of the transition term, that is, k = k 1 : where argfirstlocalmin(·) is the value of k corresponding to obtaining the first minimum; is the largest integer not exceeding , and n is the total number of the means of the IMF terms; In the third step, the classification index k 2 is determined as follows: 41) Calculate the product ET of the energy density of the k-th IMF mean term and its average period according to the EEMD decomposition result k , and the formula is as follows: where, E k is the energy density of the k-th IMF mean term, is the average period of the k-th IMF mean term, and their calculation formulas are as follows: where k = 1, 2, 3,...., n, O k is the number of extreme points of the k-th IMF term; L is the length of the k-th IMF term, that is, the number of epochs; 42) Calculate the effective coefficient R of the k-th IMF term k , and the expression is as follows: In the formula, k = 2, 3,...., n; 43) Calculate R successively according to the value range of k k , and take the k value when R k > 3 appears for the first time as the starting point of the useful term, that is, k = k 2 ; In the said Step 4, the processing steps of wavelet filtering are as follows: 51) Perform wavelet decomposition on the transitional term sequence to obtain a decomposition result composed of a group of multiple high-frequency and low-frequency terms; 52) Perform soft thresholding on the high-frequency term sequence, and do not process the low-frequency terms. The soft threshold function soft(·) is expressed as follows: where x t is the amplitude of the high-frequency term sequence, t = 1, 2, 3,...., L, where L is the sequence length, i.e., the number of epochs; sgn(·) is the step function, and when x t > 0, the result is taken as 1, when x t < 0, the result is taken as -1, and when x t = 0, the result is taken as 0; λ is the threshold, which is calculated as follows: 53) Reconstruct the low-frequency terms and the high-frequency terms after soft thresholding to obtain the filtered transitional terms; In the fifth step, the reconstructed GNSS multipath error sequence is expressed as follows: where wavelet is wavelet filtering process; IMF i (t) is the i-th IMF term obtained by EEMD, R n (t) is the remainder term; k 1 and k 2 are the indicators for distinguishing the noise term from the transitional term and the transitional term from the useful term, respectively. In the sixth step, the amplitude proportionality coefficient a r is calculated as follows: where average(·) is for mean processing; r(t) is the reference signal sequence, is the multi-path error sequence obtained by reconstruction; In the said Step 7, the calculation of restoring the multipath error amplitude is as follows: where \(s(t)\) is the multipath error after amplitude recovery; \(a\) r is the amplitude proportionality coefficient; is the multipath error sequence obtained by reconstruction.

Images