Wavelet noise reduction multi-path extraction method for improving threshold function

Abstract

The invention discloses a wavelet noise reduction multi-path extraction method for improving a threshold function, which can be used for effectively extracting and weakening multi-path errors of a global navigation satellite system (GNSS), solves the problem that signals of wavelet coefficients are discontinuous or distorted at a preset threshold, improves the traditional wavelet noise reduction effect and multi-path error modeling precision, and improves the accuracy of multi-path error modeling. Therefore, the influence of multi-path errors on GNSS positioning and orbit determination precision is effectively weakened, a feasible solution is provided for realizing high-precision positioning and orbit determination in a complex environment, and a solid foundation is laid for establishing 1mm space-time reference in the future.

Description

Technical Field

[0001] The present invention relates to the field of real-time GNSS precision positioning, particularly in urban canyon environments, where GNSS satellite signal quality is severely disrupted, making it difficult to fix GNSS ambiguities, which seriously affects the usability of positioning results. This invention provides a wavelet denoising method with an improved threshold function, addressing the problem of signal discontinuity or distortion at a preset threshold value. This method improves the wavelet denoising effect and the accuracy of multipath error modeling, effectively reducing the impact of multipath error and addressing the discontinuity and divergence of GNSS positioning in urban canyon environments. This method provides a feasible solution for achieving high-precision positioning in complex environments. Background Art

[0002] Global Navigation Satellite System (GNSS) positioning serves a wide range of applications. Growing demand is pushing GNSS positioning accuracy to the millimeter level and beyond. While significant progress has been made in high-precision GNSS positioning due to the reduction of errors in other systems, multipath error remains a major source of error because it cannot be eliminated using differential techniques.

[0003] At present, the methods for reducing multipath errors are mainly divided into three categories: (1) Selecting favorable observation conditions, relying on an open environment. In actual engineering applications, it is difficult to choose an environment that is completely unaffected by multipath. (2) Using GNSS hardware equipment, relying on GNSS antennas to reduce multipath errors. For example, measurement receiver antennas such as choke receiver antennas, single orthogonal dual linear polarization receiver antennas, and annular slot receiver antennas. These receiver antennas mainly rely on the development of antenna technology and will increase equipment costs. (3) Using GNSS data processing software based on algorithms or mathematical models to reduce multipath errors. The most commonly used are the sidereal day filter method with time-repetitive multipath effects and the semi-spherical model with spatial repetitive multipath effects.

[0004] The semi-spherical model utilizes a clock-synchronized, single-difference model with two antennas on a single machine. Under the condition that the reflection source and the satellite are relatively stationary, the multipath effect is dependent only on the satellite's elevation and azimuth, and is independent of the specific satellite. Therefore, the multipath effect is divided into grid points based on the satellite's elevation and azimuth. The multipath residuals within the grid points are averaged as the multipath correction value in that direction and modeled. This model can correct for multipath effects under the same conditions. However, to meet positioning accuracy requirements, the multipath correction model typically requires a finer grid, and high-order spherical harmonic models are computationally complex, making them difficult to implement.

[0005] Sidereal day filtering uses low-pass filtering to extract multipath errors from the residuals and construct a multipath correction model. Subsequently, by calculating the repetition period of the satellite's orbit, the multipath error corresponding to that satellite's observation period is corrected. Wavelet technology is commonly used for GNSS noise reduction and multipath error extraction. In wavelet transforms, the choice of threshold function and threshold directly impacts the quality of signal noise reduction and the accuracy of multipath error extraction. If the threshold is too small, the denoised signal will still contain noise; conversely, if the threshold is too large, important signal features will be filtered out, causing bias.

[0006] Wavelet threshold functions and wavelet thresholds are rules for modifying wavelet coefficients. Different functions reflect different strategies for processing wavelet coefficients. The most commonly used threshold functions are hard thresholding and soft thresholding. The hard thresholding method causes oscillations in the reconstructed signal due to discontinuity of wavelet coefficients at the pre-set threshold. While the soft thresholding method provides better continuity for wavelet coefficients, there is a fixed difference between the wavelet coefficients processed by the soft thresholding function and the original wavelet coefficients, which may cause some distortion in the reconstructed signal amplitude. Therefore, resolving the problem of discontinuity or distortion of wavelet coefficients at the pre-set threshold is crucial for effectively mitigating multipath errors. Summary of the Invention

[0007] The applicant's research found that when using wavelet technology to extract multipath errors, in the wavelet transform, the choice of threshold function and threshold is directly related to the quality of signal noise reduction and the accuracy of multipath error extraction. The wavelet threshold function and wavelet threshold are the rules for correcting wavelet coefficients. Different functions reflect different strategies for processing wavelet coefficients. The most commonly used threshold functions are hard threshold and soft threshold. The hard threshold method will cause oscillation of the reconstructed signal because the wavelet coefficients are discontinuous at the pre-set threshold. The soft threshold method has better continuity for the wavelet coefficients, but there is a fixed difference between the wavelet coefficients processed by the soft threshold function and the original wavelet coefficients, which may cause certain distortion of the reconstructed signal amplitude.

[0008] In order to solve the above problems, a wavelet denoising method with an improved threshold function is proposed for multipath reduction. By combining the advantages of hard threshold and soft threshold, the original signal is retained as much as possible, new oscillations in the signal after denoising are avoided, and high-precision positioning is achieved in complex environments.

[0009] In order to achieve the above object, the present invention adopts the following technical solutions:

[0010] A wavelet denoising method with an improved threshold function for multipath reduction includes the following steps:

[0011] Step Ll: Process the observation data of the rover station in real-time kinematic (RTK) mode, and perform quality control on the observation data of the reference station and the rover station to eliminate the observation data with unqualified quality. Then, the qualified observation data is used to calculate the corresponding post-difference residual.

[0012] Step L2: Convert the double-difference residual containing only multipath error and observation noise to obtain single-difference residual without receiver clock error.

[0013] Step L3: Decompose the single-difference residual by wavelet transform to generate wavelet coefficients; use the improved wavelet threshold function and heuristic threshold to correct the wavelet coefficients; and reconstruct the corrected wavelet coefficients to generate a signal without noise, thereby achieving the effect of extracting multipath.

[0014] Step L4: Estimate the repetition period of the orbit by the orbit repetition time method. Then, the multipath error of the previous day is deducted from the single-difference residual according to the orbit repetition period. Finally, the single-difference residual after multipath correction is converted into double-difference observation value, and the precise position of the rover station is estimated.

[0015] Preferably, the observation data preprocessing mainly includes pseudorange single-point positioning, pseudorange observation value checking, and correlation analysis method for gross error detection, MW (Melbourne-Wubbena) combination and geometry-free (GF) combination for cycle slip detection.

[0016] (a) Pseudorange observation value checking refers to calculating the satellite-to-station distance, i.e., satellite-geodetic distance, using the known reference station coordinates and the satellite coordinates calculated using ephemeris, and comparing the corresponding pseudorange observation value with the satellite-geodetic distance. If the pseudorange observation value differs greatly from the satellite-geodetic distance, it is considered that the observation value contains gross error and is eliminated. Correlation analysis method for gross error detection refers to detecting possible gross error in the observation value based on the correlation analysis theory of gross error after obtaining the residual of the pseudorange observation value through single-point positioning. If it is detected that the observation value contains gross error, it is eliminated.

[0017] (b) Linear combination of observation values of different frequencies or different types to obtain some combined observation values with excellent properties, which are used to detect cycle slip.

[0018] The calculation formulas of MW combination and GF combination are as follows:

[0019]

[0020] wherein, and are the observation values after MW and GF combination, respectively. and where P1 and P2 are carrier phase observations at L1 and L2 frequencies, respectively; P1 and P2 are pseudorange observations at L1 and L2 frequencies, respectively; and f1 and f2 are L1 and L2 frequencies, respectively. If the difference between new observations between adjacent epochs is greater than a threshold (0.25 cycles for GF and 0.5 cycles for MW combinations), an ambiguity flag is added, marking the end of one ambiguity arc and the beginning of another. The new ambiguity arc ends when satellite tracking ends or another cycle-skip epoch ends.

[0021] Preferably, the double-difference observations can be used to deduct some common error terms such as ionospheric error, tropospheric error, satellite orbit error, and satellite clock error, so that the integer ambiguity and the final baseline vector can be resolved, and then the double-difference residual after verification can be obtained. The original pseudorange and carrier phase observation equations are as follows:

[0022]

[0023] Wherein, subscripts r and p represent the station and pseudorange respectively; superscripts s and Represent a certain satellite and carrier phase respectively; and Original carrier phase and pseudorange observation value; dt r and dt s are the receiver and satellite clock errors respectively; c and λ are the propagation speed of light in vacuum and the carrier phase wavelength respectively; and are the tropospheric and ionospheric delays, respectively; and are the satellite-to-earth distance and phase integer ambiguity respectively; and are the hardware delays of pseudorange and phase, respectively; and denote the multipath errors of pseudorange and carrier phase respectively; and represent the noise of pseudorange and carrier phase respectively.

[0024] The double difference combination of the observation values ​​is used to eliminate the ionospheric and tropospheric errors. The calculation formula is as follows:

[0025]

[0026] Among them, the superscript i represents the selected reference satellite, the superscript j represents any other satellite; the subscript u represents the base station, and the subscript v represents the mobile station; A double difference operator representing the inter-station difference and inter-satellite difference.

[0027] The baseline used has been iteratively processed, and the geometric distance Determined by satellite orbit and station coordinates Equation (3a) can be further expressed as follows:

[0028]

[0029] (a) First, use Kalman filtering to estimate the floating point ambiguity and The corresponding standard deviation is σ L1 and σ L2 .

[0030] (b) Sort all ambiguities in order of standard deviation and resolve the ambiguity starting with the smallest one. If the floating ambiguity and its standard deviation satisfy equations (5a), 5(b), (6a) and (6b), round the floating ambiguity to a fixed ambiguity.

[0031]

[0032] σ L1 <σ t (6a)

[0033] σ L1 <σ t (6b)

[0034] in:

[0035]

[0036] relrank takes 25 as the threshold, σ t The threshold value is 0.25 weeks.

[0037] (c) After searching and fixing, the fixed fuzziness is used as the weight of 10 9 The Kalman filter is updated with virtual observations, which searches and fixes the floating-point ambiguity with the smallest standard deviation in the list until all ambiguities are fixed.

[0038] (d) Perform Kalman filtering with fixed ambiguities to estimate the unfixed ambiguities, and search and fix the floating ambiguities again. This process is iterated until all ambiguities are fixed or no more ambiguities can be fixed.

[0039] (e) Use Kalman filter to round the fixed ambiguity data to the nearest integer (IR) to obtain the fixed solution of the station. Finally, the double difference residual can be expressed as:

[0040]

[0041] Preferably, the single difference is first constrained to zero mean so that the coefficient matrix between the single difference and the double difference is full rank. Then, the double difference residual containing multipath and observation noise is transformed to obtain the single difference residual without receiver clock error. The calculation formula is as follows:

[0042]

[0043] ω i =sin 2 (θ) (10)

[0044]

[0045] Here, assuming i=1 is the reference star, is the single difference residual, ω i is the weight factor, and θ is the altitude angle.

[0046] Then, we perform an inverse operation on the coefficient matrix to obtain the single-difference residual:

[0047]

[0048] Preferably, the wavelet basis function is discretized, and then the inner product operation is performed using the discretized wavelet basis function and the single difference residual to obtain the coefficients after wavelet decomposition. The calculation formula is as follows:

[0049]

[0050] Where s(t) is the single difference residual signal to be transformed, ψ is the wavelet basis function, a is the expansion factor, and b is the translation parameter.

[0051] Using binary a=2 m and b = 2 m The wavelet basis function is discretized as follows:

[0052]

[0053] Wherein, m and n are integers.

[0054] The preferred, improved wavelet threshold function solves the problem of hard threshold being discontinuous at a pre-set threshold, and the problem of soft threshold having a fixed difference between the modified wavelet coefficients and the original wavelet coefficients. The heuristic threshold is a combination of the universal threshold and the unbiased risk estimation threshold.

[0055] The calculation formula of the improved threshold function is as follows:

[0056]

[0057] in, is the quantized wavelet coefficient, λ is the threshold parameter, and exp is the exponential operator.

[0058] The heuristic threshold calculation formula is as follows:

[0059]

[0060] Where β and γ are the calculated intermediate variable values, j represents the number of decomposition layers of the signal, and N j is the length of the single-difference residual signal. If β<γ, the universal threshold is selected as the wavelet threshold; otherwise, the smaller value between the universal threshold and the unbiased risk estimation threshold is selected as the wavelet threshold.

[0061] The general threshold is calculated as follows:

[0062]

[0063] Among them, σ j is the noise standard deviation, MAD j is the median absolute deviation of the wavelet coefficients.

[0064] The calculation formula for the unbiased risk estimate threshold is as follows:

[0065]

[0066] in, is the kth (k=1…n) risk vector; is the square of the wavelet coefficient, sorted from small to large; m corresponds to the minimum value of the risk vector.

[0067] Preferably, the satellite orbit period is first estimated using Kepler's third law based on the Earth's gravitational constant and the square root of the semi-major axis of the satellite's orbit in the broadcast ephemeris. The satellite's orbit repetition period is then estimated using the correction parameter for the satellite's mean motion in the broadcast. Finally, the orbit repetition period is subtracted from the international day to obtain the daily advance.

[0068] Using Kepler's third law, the formula for the orbital period is as follows:

[0069]

[0070] Where T is the orbital period, r is the semi-major axis of the satellite orbit, and GM = 3.986004418 × 10 14 m 3 / s 2 is the Earth's gravitational constant.

[0071] After obtaining the orbital period, perform a multiplication operation on it to obtain the orbital repetition period. The formula is as follows:

[0072]

[0073] Where T0 is the orbital repetition period, n0 and n c are the mean motion and the correction to the mean motion, respectively.

[0074] Daily lead time (T a ), the formula is as follows:

[0075] T a =86400-T0 (27)

[0076] Preferably, the multipath error of the previous day is extracted using the estimated orbit repetition period and wavelet transform to correct the single difference residual.

[0077]

[0078] in, is the single-difference residual corrected for multipath error on the second day, is the original single-difference residual on the second day, is the multipath error on the first day.

[0079] Preferably, the single-difference residual after multipath error correction is first substituted into equation (11) to obtain the double-difference residual after multipath error correction. Then, a Kalman filter is used to estimate the floating-point ambiguity and the floating-point solution of the position. Finally, the floating-point ambiguity is fixed to obtain the accurate position solution after the ambiguity is fixed.

[0080] The present invention relates to a wavelet denoising method with an improved threshold function, which is used to effectively extract multipath errors. This method solves the problem of signal discontinuity or distortion of wavelet coefficients at a preset threshold, improves the wavelet denoising effect and the accuracy of multipath error modeling, thereby effectively reducing the impact of multipath errors on the positioning accuracy of the Global Navigation Satellite System (GNSS), and provides a feasible solution for achieving high-precision positioning in complex environments. BRIEF DESCRIPTION OF THE DRAWINGS

[0081] Figure 1 Flowchart of a wavelet denoising multipath extraction method with an improved threshold function DETAILED DESCRIPTION

[0082] The present invention proposes a wavelet denoising method with an improved threshold function for multipath reduction, processes the observation data of the rover in real-time kinematic positioning (RTK) mode, performs quality control on the observation data preprocessing stage of the base station and the rover, and eliminates observation data with unqualified observation quality; then uses the qualified observation data to calculate their corresponding posterior double difference residuals.

[0083] The double-difference residual containing only multipath error and observation noise is transformed to obtain the single-difference residual without receiver clock error.

[0084] The single difference residual is decomposed by wavelet transform to generate wavelet coefficients; the wavelet coefficients are corrected using the improved wavelet threshold function and heuristic threshold; the corrected wavelet coefficients are reconstructed to generate a noise-free signal, thereby achieving the effect of extracting multipath.

[0085] The orbit repetition time method is used to estimate the orbit repetition period. The multipath error from the previous day is then subtracted from the single-difference residuals of the orbit repetition period. Finally, the multipath-corrected single-difference residuals are converted into double-difference observations to estimate the precise position of the rover.

[0086] The following flowcharts illustrate the specific embodiments of the present invention in detail. The accompanying drawings are in very simplified form and are only used to clearly illustrate the implementation of the present invention.

[0087] The present invention provides a wavelet denoising method with an improved threshold function for multipath reduction. The specific implementation method is as follows: Figure 1 As shown:

[0088] Step S1: The observation data preprocessing mainly includes pseudorange single point positioning, pseudorange observation value verification and correlation analysis method for gross error detection, MW (Melbourne-Wubbena) combination and geometry-free combination (GF) for cycle slip detection.

[0089] (a) Pseudorange observation verification involves calculating the satellite-to-station distance, i.e., the satellite-to-earth distance, using the known reference station coordinates and the satellite coordinates calculated using the ephemeris. The pseudorange observation corresponding to the satellite-to-earth distance is then compared. If the pseudorange observation differs significantly from the satellite-to-earth distance, the observation is considered to contain a gross error and is discarded. Gross error detection using correlation analysis involves detecting possible gross errors in the pseudorange observations obtained through single-point positioning using the gross error theory of correlation analysis. Any observed values ​​found to contain a gross error are discarded.

[0090] (b) Linearly combining observations of different frequencies or types to obtain some combined observations with excellent characteristics. These new observations are used to detect cycle slips.

[0091] The calculation formulas for MW combination and GF combination are as follows:

[0092]

[0093] in, and are the observed values ​​after combining MW and GF respectively; and Where: P1 and P2 are carrier phase observations at L1 and L2 frequencies, respectively; P1 and P2 are pseudorange observations at L1 and L2 frequencies, respectively; f1 and f2 are L1 and L2 frequencies, respectively. If the difference between new observations between adjacent epochs is greater than a threshold (0.25 cycles for GF and 0.5 cycles for MW combinations), an ambiguity flag is added, marking the end of one ambiguity arc and the beginning of another. The new ambiguity arc ends at the end of satellite tracking or at the end of another cycle-skip epoch. The preprocessed observations from S1 are provided to S2.

[0094] Step S2: Using the double-difference observations, some common error terms such as ionospheric error, tropospheric error, satellite orbit error, and satellite clock error can be deducted, thereby resolving the integer ambiguity and the final baseline vector, and then obtaining the posterior double-difference residual. The original pseudorange and carrier phase observation equations are as follows:

[0095]

[0096] Wherein, subscripts r and p represent the station and pseudorange respectively; superscripts s and Represent a certain satellite and carrier phase respectively; and Original carrier phase and pseudorange observation value; dt r and dt s are the receiver and satellite clock errors respectively; c and λ are the propagation speed of light in vacuum and the carrier phase wavelength respectively; and are the tropospheric and ionospheric delays, respectively; and are the satellite-to-earth distance and phase integer ambiguity respectively; and are the hardware delays of pseudorange and phase, respectively; and denote the multipath errors of pseudorange and carrier phase respectively; and represent the noise of pseudorange and carrier phase respectively.

[0097] The double difference combination of the observation values ​​is used to eliminate the ionospheric and tropospheric errors. The calculation formula is as follows:

[0098]

[0099] Among them, the superscript i represents the selected reference satellite, the superscript j represents any other satellite; the subscript u represents the base station, and the subscript v represents the mobile station; A double difference operator representing the inter-station difference and inter-satellite difference.

[0100] The baseline used has been iteratively processed, and the geometric distance Determined by satellite orbit and station coordinates Equation (3a) can be further expressed as follows:

[0101]

[0102] (a) First, use Kalman filtering to estimate the floating point ambiguity and The corresponding standard deviation is σ L1 and σ L2 .

[0103] (b) Sort all ambiguities in order of standard deviation and resolve the ambiguity starting with the smallest one. If the floating ambiguity and its standard deviation satisfy equations (5a), 5(b), (6a) and (6b), round the floating ambiguity to a fixed ambiguity.

[0104]

[0105] σ L1 <σ t (6a)

[0106] σ L1 <σ t (6b)

[0107] in:

[0108]

[0109] relrank takes 25 as the threshold, σ t The threshold value is 0.25 weeks.

[0110] (c) After searching and fixing, the fixed fuzziness is used as the weight of 10 9 The Kalman filter is updated with virtual observations, which searches and fixes the floating-point ambiguity with the smallest standard deviation in the list until all ambiguities are fixed.

[0111] (d) Perform Kalman filtering with fixed ambiguities to estimate the unfixed ambiguities, and search and fix the floating ambiguities again. This process is iterated until all ambiguities are fixed or no more ambiguities can be fixed.

[0112] (e) Use the Kalman filter integer rounding (IR) method to resolve the fixed ambiguity data and obtain the fixed solution of the station. Finally, the double difference residual can be expressed as:

[0113]

[0114] Step S3: Convert double-difference residuals into single-difference residuals. First, constrain the single-difference residuals to zero mean, and obtain the conversion relationship between single difference and double difference as follows:

[0115]

[0116] ω i =sin 2 (θ) (10)

[0117]

[0118] Here, assuming i=1 is the reference star, is the single difference residual, ω i is the weight factor, and θ is the altitude angle.

[0119] Then, the coefficient matrix is ​​inverted to obtain the single-difference residual, providing S4.

[0120]

[0121] Step S4: Discretize the wavelet basis function, then use the discretized wavelet basis function and the single difference residual to perform inner product operation to obtain the coefficients after wavelet decomposition, and provide them to S5. The calculation formula is as follows:

[0122]

[0123] Where s(t) is the single difference residual signal to be transformed, ψ is the wavelet basis function, a is the expansion factor, and b is the translation parameter.

[0124] Using binary a=2 m and b = 2 m The wavelet basis function is discretized as follows:

[0125]

[0126] Wherein, m and n are integers.

[0127] Then the improved wavelet threshold function is used to correct the decomposed wavelet coefficients. The calculation formula is as follows:

[0128]

[0129] in, is the quantized wavelet coefficient, λ is the threshold parameter, and exp is the exponential operator.

[0130] Step S5: The improved wavelet threshold function solves the problem of discontinuity of the hard threshold at the pre-set threshold, and the problem of a fixed difference between the modified wavelet coefficients and the original wavelet coefficients in the soft threshold. The heuristic threshold is a combination of the universal threshold and the unbiased risk estimation threshold.

[0131] The calculation formula of the improved threshold function is as follows:

[0132]

[0133] in, is the quantized wavelet coefficient, λ is the threshold parameter, and exp is the exponential operator.

[0134] The heuristic threshold calculation formula is as follows:

[0135]

[0136] Where β and γ are the calculated intermediate variable values, j represents the number of decomposition layers of the signal, and N j is the length of the single-difference residual signal. If β<γ, the universal threshold is selected as the wavelet threshold; otherwise, the smaller value between the universal threshold and the unbiased risk estimation threshold is selected as the wavelet threshold.

[0137] The general threshold is calculated as follows:

[0138]

[0139] Among them, σ j is the noise standard deviation, MAD j is the median absolute deviation of the wavelet coefficients.

[0140] The calculation formula for the unbiased risk estimate threshold is as follows:

[0141]

[0142] in, is the kth (k=1…n) risk vector; is the square of the wavelet coefficient, sorted from small to large; m corresponds to the minimum value of the risk vector.

[0143] Finally, the corrected wavelet coefficients are recombined to generate a noise-free signal, thereby achieving the effect of extracting multipath errors and providing it to S7.

[0144] Step S6: Using the Earth's gravitational constant and the square root of the satellite's semi-major axis in the broadcast ephemeris, Kepler's third law is used to estimate the satellite's orbital period. The satellite's orbital repetition period is then estimated using the correction for the satellite's mean motion in the broadcast. Finally, the orbital repetition period is subtracted from the international day to obtain the daily advance.

[0145] Using Kepler's third law, the formula for the orbital period is as follows:

[0146]

[0147] Where T is the orbital period, r is the semi-major axis of the satellite orbit, and GM = 3.986004418 × 10 14 m 3 / s 2 is the Earth's gravitational constant.

[0148] After obtaining the orbital period, perform a multiplication operation on it to obtain the orbital repetition period. The formula is as follows:

[0149]

[0150] Where T0 is the orbital repetition period, n0 and n c are the mean motion and the correction to the mean motion, respectively.

[0151] Daily lead time (T a ) is provided to S7, and the formula is as follows:

[0152] T a =86400-T0 (27)

[0153] Step S7: Use the estimated orbit repetition period and wavelet transform to extract the multipath error of the previous day and correct the single difference residual.

[0154]

[0155] in, is the single-difference residual corrected for multipath error on the second day, is the original single-difference residual on the second day, is the multipath error on the first day. The single-difference residual corrected for the multipath error is provided to S8.

[0156] Step S8: Substitute the single-difference residual after multipath error correction into equation (11) to obtain the double-difference residual after multipath error correction.

[0157] Step S9: Use Kalman filtering to estimate the floating-point ambiguity and the floating-point solution of the position, and finally fix the floating-point ambiguity to obtain the precise position solution after the ambiguity is fixed.

[0158] The above description is merely a preferred embodiment of the present invention and does not limit the present invention in any way. Any person skilled in the art who, without departing from the scope of the present invention, makes any equivalent substitution, modification, or other changes to the technical solution and technical content disclosed in the present invention shall be deemed to be within the scope of the present invention and still fall within the scope of protection of the present invention.

Claims

1. A wavelet denoising multipath extraction method with an improved threshold function, characterized in that: The main steps are as follows: (1) Step L1: When processing the observation data of the rover in the GNSS real-time kinematic positioning (RTK) mode, the multipath effect must generally be considered. To this end, the observation data of the base station and the rover are first preprocessed to eliminate data with unqualified observation quality. Then, the qualified observation data are used for kinematic positioning processing, and their corresponding posterior double-difference residuals are calculated. The double-difference residuals are considered to contain only multipath errors and observation noise. (2) Step L2: Convert the double-difference residual containing only multipath error and observation noise to obtain the single-difference residual without receiver clock error. (3) Step L3: Decompose the single-difference residual using wavelet transform to generate wavelet coefficients; The improved wavelet threshold function and heuristic threshold proposed in this patent are used to correct the wavelet coefficients; the corrected wavelet coefficients are reconstructed to generate a noise-free signal, thereby obtaining a high-precision multipath error and achieving the purpose and effect of extracting multipath. (4) Step L4: Calculate the orbit repetition period using the orbit repetition time method; then, based on the orbit repetition period, use the single-difference residual to deduct the multipath error of the previous day; Finally, the single-difference residuals after multipath correction are converted into double-difference observations to estimate the precise position of the mobile station.

2. The wavelet denoising multipath extraction method with an improved threshold function as claimed in claim 1 is used for multipath reduction, characterized in that: The quality control of the pre-processing stage in step L1 includes: The observation data preprocessing mainly includes pseudorange single point positioning, pseudorange observation value verification and correlation analysis method for pseudorange data gross error detection, and MW (Melbourne-Wubbena) combination and geometry-free (GF) combination for phase data cycle slip detection. (a) Pseudorange observation verification involves calculating the satellite-to-station distance, i.e., the satellite-to-earth distance, using the known reference station coordinates and the satellite coordinates calculated using the ephemeris. The pseudorange observation corresponding to the satellite-to-earth distance is then compared. If the pseudorange observation differs significantly from the satellite-to-earth distance, the observation is considered to contain a gross error and is discarded. Gross error detection using correlation analysis involves detecting possible gross errors in the pseudorange observations obtained through single-point positioning using the gross error theory of correlation analysis. Any observed values ​​found to contain a gross error are discarded. (b) Linearly combining observations of different frequencies or types to obtain some combined observations with certain excellent properties. These combined observations can be effectively used to detect cycle slips. The calculation formulas for MW combination and GF combination are as follows: in, and are the observation values ​​after combining MW and GF respectively; and Where P1 and P2 are carrier phase observations at L1 and L2 frequencies, respectively; P1 and P2 are pseudorange observations at L1 and L2 frequency bands, respectively; and f1 and f2 are L1 and L2 frequencies, respectively. If the difference between combined observations between adjacent epochs exceeds a threshold (0.25 cycles for GF and 0.5 cycles for MW), an ambiguity flag is added, marking the end of one ambiguity arc and the beginning of a new one. The new ambiguity arc ends when satellite tracking ends or another cycle-skip epoch is encountered, and this cycle repeats until the ambiguities of all data arcs are flagged.

3. The wavelet denoising method with an improved threshold function as claimed in claim 1 is used for multipath reduction, characterized in that: The double difference residual in step L1 is calculated as follows: Using double-difference observations, some common error terms such as ionospheric error, tropospheric error, satellite orbit error, and satellite clock error can be deducted, so that the whole-cycle ambiguity and baseline vector can be solved more accurately, and then the double-difference residual can be obtained. Assume that the original pseudorange and carrier phase observation equations are as follows: Among them, the subscripts r, p and Represent the measuring station, pseudorange and carrier phase respectively; the superscript s represents a satellite; and are the original carrier phase and pseudorange observation values; dt r and dt s are the receiver and satellite clock errors respectively; c and λ are the propagation speed of light in vacuum and the carrier phase wavelength respectively; and are the tropospheric and ionospheric delays, respectively; and are the satellite-to-earth distance and phase integer ambiguity respectively; and are the hardware delays of pseudorange and phase, respectively; and denote the multipath errors of pseudorange and carrier phase respectively; and represent the noise of pseudorange and carrier phase respectively. The double difference combination of the observation values ​​is used to eliminate the ionospheric and tropospheric errors. The calculation formula is as follows: Among them, the superscript i represents the selected reference satellite, the superscript j represents any other satellite; the subscript u represents the base station, and the subscript v represents the mobile station; A double difference operator representing the inter-station difference and inter-satellite difference. The baseline used has been iteratively processed, and the geometric distance Determined by satellite orbit and station coordinates Equation (3a) can be further expressed as follows: (a) First, use Kalman filtering to estimate the floating point ambiguity and The corresponding standard deviation is σ L1 and σ L2 . (b) Sort all ambiguities in order of standard deviation and resolve the ambiguity starting with the smallest one. If the floating ambiguity and its standard deviation satisfy equations (5a), 5(b), (6a) and (6b), round the floating ambiguity to a fixed ambiguity. s L1 <s t (6a) s L1 <s L (6b) in: relrank takes 25 as the threshold, σ t The threshold value is 0.25 weeks. (c) After searching and fixing, the fixed fuzziness is used as the weight of 10 9 The Kalman filter is updated with virtual observations, which searches and fixes the floating-point ambiguity with the smallest standard deviation in the list until all ambiguities are fixed. (d) Perform Kalman filtering with fixed ambiguities to estimate the unfixed ambiguities, and search and fix the floating ambiguities again. This process is iterated until all ambiguities are fixed or no more ambiguities can be fixed. (e) Use the Kalman filter integer rounding (IR) method to resolve the fixed ambiguity data and obtain the fixed solution of the station. Finally, the double difference residual can be expressed as:

4. The wavelet denoising method with an improved threshold function as claimed in claim 1 is used for multipath reduction, characterized in that: The calculation of converting the double difference residual in step L2 into a single difference residual is as follows: First, the single difference is constrained to zero mean so that the coefficient matrix between the single difference and the double difference is full rank. Then, the double difference residual containing multipath and observation noise is transformed to obtain the single difference residual without receiver clock error. The calculation formula is as follows: oh i =sin 2 (i) (10) Here, assuming i=1 is the reference star, is the single difference residual, ω i is the weight factor, and θ is the altitude angle. Then, we perform an inverse operation on the coefficient matrix to obtain the single-difference residual:

5. The wavelet denoising method with an improved threshold function as claimed in claim 1 is used for multipath reduction, characterized in that: The wavelet transform in step L3 decomposes the single difference residual to generate the wavelet coefficients as follows: First, the wavelet basis function is discretized, and then the inner product operation is performed using the discretized wavelet basis function and the single difference residual to obtain the coefficients after wavelet decomposition. The calculation formula is as follows: Where s(t) is the single difference residual signal to be transformed, ψ is the wavelet basis function, a is the expansion factor, and b is the translation parameter. Using binary a=2 m and b = 2 m The wavelet basis function is discretized as follows: Wherein, m and n are integers.

6. The wavelet denoising method with an improved threshold function as claimed in claim 1 is used for multipath reduction, characterized in that: In step L3, the improved wavelet threshold function and the heuristic threshold are used to correct the decomposed wavelet coefficients, and the calculation is as follows: The improved wavelet threshold function solves the problem of hard threshold being discontinuous at the pre-set threshold, and the problem of soft threshold having a fixed difference between the modified wavelet coefficients and the original wavelet coefficients. The heuristic threshold is a combination of the universal threshold and the unbiased risk estimation threshold. The calculation formula of the improved threshold function is as follows: in, is the quantized wavelet coefficient, λ is the threshold parameter, and exp is the exponential operator. The heuristic threshold calculation formula is as follows: Where β and γ are the calculated intermediate variable values, j represents the number of decomposition layers of the signal, and N j is the length of the single-difference residual signal. If β<γ, the universal threshold is selected as the wavelet threshold; otherwise, the smaller value between the universal threshold and the unbiased risk estimation threshold is selected as the wavelet threshold. The general threshold is calculated as follows: Among them, σ j is the noise standard deviation, MAD j is the median absolute deviation of the wavelet coefficients. The calculation formula for the unbiased risk estimate threshold is as follows: in, is the kth (k=1…n) risk vector; is the square of the wavelet coefficient, sorted from small to large; m corresponds to the minimum value of the risk vector.

7. The wavelet denoising method with an improved threshold function as claimed in claim 1 is used for multipath reduction, characterized in that: The calculation of the orbit repetition period estimation in step L4 is as follows: First, Kepler's third law is used to estimate the satellite's orbital period based on the Earth's gravitational constant and the square root of the semi-major axis of the satellite's orbit as reported in the broadcast ephemeris. The satellite's orbital repetition period is then estimated using the correction for the satellite's mean motion as reported in the broadcast. Finally, the daily advance is calculated by subtracting the orbital repetition period from the international day. Using Kepler's third law, the formula for the orbital period is as follows: Where T is the orbital period, r is the semi-major axis of the satellite orbit, and GM = 3.986004418 × 10 14 m 3 / s 2 is the Earth's gravitational constant. After obtaining the orbital period, perform a multiplication operation on it to obtain the orbital repetition period. The formula is as follows: Where T0 is the orbital repetition period, n0 and n c are the mean motion and the correction to the mean motion, respectively. Daily lead time (T a ), the formula is as follows: T a =86400-T0 (27)。 8. The wavelet denoising method with an improved threshold function as claimed in claim 1 is used for multipath reduction, characterized in that: In the L4 step, the multipath error extracted the previous day is deducted from the single difference residual. The calculation formula is as follows: The multipath error of the previous day is extracted using the estimated orbit repetition period and wavelet transform to correct the single difference residual. in, is the single-difference residual corrected for multipath error on the second day, is the original single-difference residual on the second day, is the multipath error on the first day.

9. The wavelet denoising method with an improved threshold function as claimed in claim 1 is used for multipath reduction, characterized in that: The step L4 converts the corrected single-difference residuals into double-difference observations and estimates the precise position of the mobile station, including: First, the single-difference residual after multipath error correction is substituted into Equation (11) to obtain the double-difference residual after multipath error correction. Then, a Kalman filter is used to estimate the floating-point ambiguity and the floating-point solution of the position. Finally, the floating-point ambiguity is fixed to obtain the accurate position solution after the ambiguity is fixed.

Images