Wavelet denoising multi-path extraction method with improved threshold function

Abstract

The application discloses a wavelet denoising multi-path extraction method with improved threshold function, which can be used for effective extraction and weakening of global navigation satellite system (GNSS) multi-path error, solves the problem of signal discontinuity or distortion of wavelet coefficients at a preset threshold, improves traditional wavelet denoising effect and multi-path error modeling precision, effectively weakens the influence of multi-path error on GNSS positioning and orbit determination precision, provides a feasible solution for realizing high-precision positioning and orbit determination in a complex environment, and lays a solid foundation for establishing a 1mm space-time reference in the future.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of real-time GNSS precise positioning, especially aiming at the environment of urban canyon, GNSS satellite signal quality is seriously disturbed, which makes GNSS ambiguity difficult to fix, thereby seriously affecting the usability of positioning results. The present application provides a wavelet denoising method for improving threshold function, solves the problem of signal discontinuity or distortion of wavelet coefficients at the preset threshold, improves the wavelet denoising effect and the precision of multipath error modeling, thereby effectively weakening the influence of multipath error, making up for the discontinuity and divergence of GNSS positioning in the environment of urban canyon, and providing a feasible solution for high-precision positioning in complex environment. BACKGROUND

[0002] Global Navigation Satellite System (GNSS) positioning serves a variety of applications. The increasing demand pushes GNSS positioning accuracy to millimeter level or even higher. Although the reduction of other system errors has made great progress in GNSS high-precision positioning, multipath error is still one of the main error sources because it cannot be eliminated by differential technology.

[0003] Currently, the methods for weakening multipath error mainly fall into three categories, (1) selecting favorable observation conditions, relying on open environment. In actual engineering applications, it is difficult to choose an environment completely free of multipath effects. (2) Using GNSS hardware equipment, relying on GNSS antennas to reduce multipath error. For example, choke ring receiver antennas, single orthogonal bilinear polarization receiver antennas, ring slot receiver antennas, and other measurement type receiver antennas. These receiver antennas mainly rely on the development of antenna technology, while increasing the cost of equipment. (3) Using GNSS data processing software based on algorithms or mathematical models to reduce multipath error, the most commonly used are the constant star-day filtering method which has time repetition characteristics of multipath effect and the half-sphere model which has spatial repetition characteristics of multipath effect.

[0004] The half-sphere model is a clock-synchronized single-difference model with two antennas. In the case of a relatively static reflection source and satellite, the relative position between the multipath reflection source and the antenna is constant. The multipath effect is only related to the satellite elevation angle and azimuth angle, and is independent of the specific satellite. Therefore, the multipath effect can be divided into grid points according to the satellite elevation angle and azimuth angle, and the average of the multipath residual within the grid point is taken as the multipath correction value in that direction and is modeled. The multipath effect in the same environment can be corrected through the model. However, in order to meet the demand for positioning accuracy, the multipath correction model needs to be divided into more fine grids, and the calculation complexity of high-order spherical harmonic model is high, so it is difficult to implement.

[0005] The star day filter is to extract the multipath error in the residual error by using low-pass filter to construct a multipath correction model. Then, the repetition period of the satellite orbit is calculated to correct the multipath error of the observation period corresponding to the satellite. Wavelet technology is usually used for GNSS noise reduction and multipath error extraction. In wavelet transform, the selection of threshold function and threshold directly relates to the quality of signal denoising and the accuracy of multipath error extraction. If the threshold is too small, the denoised signal still has noise; on the contrary, if the threshold is too large, important signal characteristics are filtered out, causing deviation.

[0006] The wavelet threshold function and the wavelet threshold are the rules for correcting the wavelet coefficients, and different functions represent different strategies for processing the wavelet coefficients. The most commonly used threshold functions are hard threshold and soft threshold. The hard threshold method is discontinuous at the preset threshold of the wavelet coefficient, which may cause oscillation of the reconstructed signal. The soft threshold method has good continuity for the wavelet coefficient, but there is a fixed difference between the wavelet coefficient processed by the soft threshold function and the original wavelet coefficient, which may cause distortion of the amplitude of the reconstructed signal. Therefore, it is very important to solve the problem of signal discontinuity or distortion of the wavelet coefficient at the preset threshold for effectively weakening the multipath error. SUMMARY

[0007] The applicant found that in the wavelet transform, the selection of threshold function and threshold directly relates to the quality of signal denoising and the accuracy of multipath error extraction. The wavelet threshold function and the wavelet threshold are the rules for correcting the wavelet coefficients, and different functions represent different strategies for processing the wavelet coefficients. The most commonly used threshold functions are hard threshold and soft threshold. The hard threshold method is discontinuous at the preset threshold of the wavelet coefficient, which may cause oscillation of the reconstructed signal. The soft threshold method has good continuity for the wavelet coefficient, but there is a fixed difference between the wavelet coefficient processed by the soft threshold function and the original wavelet coefficient, which may cause distortion of the amplitude of the reconstructed signal.

[0008] In order to solve the above problems, a wavelet denoising method with improved threshold function is proposed for multipath weakening, which combines the advantages of hard threshold and soft threshold to preserve the original signal as much as possible, avoid new oscillation of the denoised signal and achieve high-precision positioning in complex environment.

[0009] In order to achieve the above purpose, the application adopts the following technical scheme:

[0010] A wavelet denoising method with improved threshold function for multipath weakening comprises 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 P1and P2are the pseudorange observations on L1 and L2 bands, respectively; and f1and f2are the L1 and L2 frequencies, respectively. If the difference between the new observations is greater than a threshold (0.25 cycles for GF, 0.5 cycles for MW combination) between adjacent epochs, a flag is added, which means that one ambiguity arc ends and another ambiguity arc begins. The new ambiguity arc will end at the end of satellite tracking or another cycle slip epoch.

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

[0022]

[0023] wherein the subscripts r and p represent the station and the pseudorange, respectively; the superscripts s and represent a satellite and a carrier phase, respectively; and the original carrier phase and pseudorange observations; 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 geodetic distance and the phase integer ambiguity, respectively; and are the hardware delays of the pseudorange and the phase, respectively; and represent the multipath errors of the pseudorange and the carrier phase, respectively; and represent the noises of the pseudorange and the carrier phase, respectively.

[0024] The double difference combination of the observations is used to eliminate the ionospheric and tropospheric errors, and the calculation formula is as follows:

[0025]

[0026] wherein the superscripts i and j represent the selected reference satellite and any other satellite, respectively; the subscripts u and v represent the reference station and the rover station, respectively; represent the double difference operators of the inter-station difference and the inter-satellite difference.

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

[0028]

[0029] (a) First, estimate the float ambiguities with Kalman filter and The corresponding standard deviation is σ L1 and σ L2 .

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

[0031]

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

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

[0034] where:

[0035]

[0036] relrank takes 25 as the threshold for σ t and 0.25 week as the threshold.

[0037] (c) After searching and fixing, use the fixed ambiguities as virtual observations with weight 10 9 to update the Kalman filter, and search and fix the float ambiguity with the smallest standard deviation in the search and fix list until all ambiguities are fixed.

[0038] (d) Use the fixed ambiguities to estimate the unfixed ambiguities with Kalman filter, and search and fix the float ambiguities again. This process will be iterated until all ambiguities are fixed or no ambiguity can be fixed.

[0039] (e) Use Kalman filter with Integer Rounding (IR) to solve the data with fixed ambiguities to obtain the fixed solution of the station. Finally, the double difference residuals can be expressed as:

[0040]

[0041] Preferably, firstly, zero-mean constraint is performed on single difference, so that the coefficient matrix between single difference and double difference is full rank, then double difference residual containing multipath and observation noise is converted, so as to obtain single difference residual without receiver clock error, and the calculation formula is as follows:

[0042]

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

[0044]

[0045] Wherein, it is assumed that i=1 is a reference star, is single difference residual, ω i is a weight factor, and θ is an elevation angle.

[0046] Then, inverse operation is performed on the coefficient matrix, so as to obtain single difference residual:

[0047]

[0048] Preferably, the wavelet base function is discretized, then inner product operation is performed on the discretized wavelet base function and single difference residual, so as to obtain wavelet decomposition coefficient, and the calculation formula is as follows:

[0049]

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

[0051] The wavelet base function is discretized by using binary a=2 m and b=2 m , and the formula is as follows:

[0052]

[0053] Wherein, m and n are integers respectively.

[0054] Preferably, the improved wavelet threshold function solves the problem that the hard threshold value is discontinuous at the preset threshold value, and the problem that the soft threshold value has a fixed difference between the corrected wavelet coefficient and the original wavelet coefficient. The heuristic threshold value is a combination of a general threshold value and an unbiased risk estimation threshold value.

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

[0056]

[0057] Wherein, For quantizing the wavelet coefficients, λ 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 decomposition level of the signal, 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 calculation formula of the universal threshold is as follows:

[0062]

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

[0064] The calculation formula of the unbiased risk estimation threshold is as follows:

[0065]

[0066] where, is the kth (k = 1…n) risk vector; is the square of the wavelet coefficients sorted in descending order; m corresponds to the minimum value of the risk vector.

[0067] Preferably, first, the period of the satellite orbit is estimated according to the Earth's gravitational constant and the square root of the semi-major axis of the satellite orbit in the broadcast ephemeris using Kepler's third law. Then, the orbit repeat period of the satellite is estimated according to the correction of the mean motion of the satellite in the broadcast. Finally, the orbit repeat period is subtracted from the international day to obtain the time in advance per day.

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

[0069]

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

[0071] After obtaining the orbit period, a multiple operation is performed on it to obtain the orbit repeat period, and 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] The daily advance time (T a ) is given by the following equation:

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

[0076] Preferably, the previous day's multipath error is extracted using the estimated orbital repetition period and wavelet transform to correct the single-difference residuals.

[0077]

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

[0079] Preferably, the single-difference residual corrected for multipath error is first substituted into equation (11) to obtain the double-difference residual corrected for multipath error. Then, the Kalman filter is used to estimate the float ambiguity and the float solution of the position, and finally the float ambiguity is fixed to obtain the precise position solution after fixing the ambiguity.

[0080] The present application relates to a wavelet denoising method for improving threshold function, which is used for effectively extracting multipath error, solving the problem of signal discontinuity or distortion at the preset threshold of wavelet coefficient, improving the wavelet denoising effect and the precision of multipath error modeling, thereby effectively weakening the influence of multipath error on the positioning precision of global navigation satellite system (GNSS), and providing a feasible solution for realizing high-precision positioning in a complex environment. BRIEF DESCRIPTION OF DRAWINGS

[0081] Figure 1 A flowchart of a wavelet denoising multipath extraction method for improving threshold function DETAILED DESCRIPTION

[0082] The present application proposes a wavelet denoising method for improving threshold function for multipath weakening, which is used for processing the observation data of a rover station in real-time kinematic (RTK) mode, performing quality control on the observation data of the reference station and the rover station in the observation data preprocessing stage, and eliminating the observation data with unqualified quality; and then calculating the corresponding post-difference residuals of the qualified observation data.

[0083] The double-difference residuals containing only multipath error and observation noise are converted to obtain single-difference residuals without receiver clock error.

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

[0085] The orbit repetition time method is used to estimate the repetition period of the orbit. Then, the single difference residual is deducted from the multipath error of the previous day according to the orbit repetition period. Finally, the single difference residual after the multipath correction is converted into double difference observations, and the precise position of the rover station is estimated.

[0086] The following flow chart details the specific embodiments of the present application. The accompanying drawings are very simplified and only serve to clearly illustrate the purpose of the present application.

[0087] The present application provides a wavelet denoising method for improving the threshold function for multipath weakening, and the specific implementation method is as shown in Figure 1 .

[0088] S1 step: pre-processing of observation data 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) for cycle slip detection.

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

[0090] (b) Linear combination of observation values of different frequencies or different types to obtain some combined observation values with excellent properties. These new observation values are used to detect cycle slips.

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

[0092]

[0093] wherein, and are the observation values after MW and GF combination, respectively; and P1and P2are the pseudorange observations on L1 and L2 bands, respectively; f1and f2are the L1 and L2 frequencies, respectively. If the difference between the new observations is greater than a threshold (0.25 cycles for GF, 0.5 cycles for MW combination) between adjacent epochs, a flag is added, which means that one ambiguity arc ends and another ambiguity arc begins. The new ambiguity arc will end at the end of satellite tracking or another cycle slip epoch. The S1 pre-processed observations are provided to S2.

[0094] S2 step: 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 integer ambiguity and the final baseline vector can be solved, and then the posteriori double difference residuals can be obtained. The original pseudorange and carrier phase observation equations are as follows:

[0095]

[0096] wherein the subscripts r and p represent the station and the pseudorange, respectively; the superscripts s and represent a satellite and a carrier phase, respectively; and original carrier phase and pseudorange observations; 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 geodetic distance and the phase integer ambiguity, respectively; and are the hardware delays of the pseudorange and the phase, respectively; and represent the multipath errors of the pseudorange and the carrier phase, respectively; and represent the noises of the pseudorange and the carrier phase, respectively.

[0097] Using the observations, the ionospheric and tropospheric errors are eliminated by double difference combination, and the calculation formula is as follows:

[0098]

[0099] wherein the superscripts i represent the selected reference satellite, and the superscripts j represent any other satellite; the subscripts u represent the reference station, and the subscripts v represent the rover station; represent the double difference operators of the inter-station difference and the inter-satellite difference.

[0100] The baseline used has been processed by iteration, and the geometric distance Determination from satellite orbits and station coordinates Equation (3a) is further expressed as follows:

[0101]

[0102] (a) First, estimate the float ambiguities using Kalman filter and The corresponding standard deviation is σ L1 and σ L2 .

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

[0104]

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

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

[0107] where:

[0108]

[0109] relrank takes 25 as the threshold for σ t and 0.25 cycle as the threshold.

[0110] (c) After searching and fixing, use the fixed ambiguities as virtual observations with weight 10 9 in the Kalman filter to update the float ambiguities. The float ambiguity with the smallest standard deviation in the searching and fixing list is searched and fixed until all ambiguities are fixed.

[0111] (d) Estimate the unfixed ambiguities using the Kalman filter with the fixed ambiguities, and search and fix the float ambiguities again. This process is iterated until all ambiguities are fixed or no ambiguity can be fixed.

[0112] (e) Use the Kalman filter Integer Rounding (IR) method to solve the data with fixed ambiguities to obtain the fixed solution of the station. Finally, the double-difference residuals can be expressed as:

[0113]

[0114] S3 step: double difference residual is converted into single difference residual, first zero mean constraint is carried out to single difference residual, and conversion relationship of single difference and double difference is obtained as follows:

[0115]

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

[0117]

[0118] Wherein, it is assumed that i=1 is a reference star, is single difference residual, ω i is weight factor, and θ is elevation angle.

[0119] Then, inverse operation is carried out to coefficient matrix, and single difference residual is obtained, and S4 is provided.

[0120]

[0121] S4 step: wavelet base function is discretized, then inner product operation is carried out to discretized wavelet base function and single difference residual, and coefficient after wavelet decomposition is obtained, and S5 is provided, and calculation formula is as follows:

[0122]

[0123] Wherein, s(t) is single difference residual signal needing transformation, ψ is wavelet base function, a is expansion factor, and b is translation parameter.

[0124] Binary a=2 m and b=2 m are adopted to discretize wavelet base function, and formula is as follows:

[0125]

[0126] Wherein, m and n are integers respectively.

[0127] Then, improved wavelet threshold function is adopted to modify wavelet coefficient after decomposition, and calculation formula is as follows:

[0128]

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

[0130] S5 step: the improved wavelet threshold function solves the problem of hard threshold that is not continuous at the preset threshold, and the problem of soft threshold that has a fixed difference between the corrected wavelet coefficient and the original wavelet coefficient. 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] Wherein, is the quantized wavelet coefficient, λ is the threshold parameter, and exp is the exponential operator.

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

[0135]

[0136] Wherein, β and γ are the calculated intermediate variable values, j represents the decomposition level of the signal, 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 calculation formula of the universal threshold is as follows:

[0138]

[0139] Wherein, σ j is the noise standard deviation, and MAD j is the median absolute deviation of the wavelet coefficient.

[0140] The calculation formula of the unbiased risk estimation threshold is as follows:

[0141]

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

[0143] Finally, the corrected wavelet coefficient is recombined to generate a noise-free signal, thereby achieving the effect of extracting the multi-path error and providing S7.

[0144] S6 step: according to the earth's gravitational constant and the square root parameter of the satellite orbit in the broadcast ephemeris, the period of the satellite orbit is estimated by using Kepler's third law. Then, according to the correction amount parameter of the average motion of the satellite in the broadcast, the orbit repetition period of the satellite is estimated. Finally, the orbit repetition period and the international day difference are subtracted to obtain the daily advance time.

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

[0146]

[0147] where T is the orbit period, r is the semi-major axis of the satellite orbit, and GM = 3.986004418 x 10 14 m 3 / s 2 is the gravitational constant of the earth.

[0148] After obtaining the orbit period, the multiple operation is performed on it to obtain the orbit repeat period, and the formula is as follows:

[0149]

[0150] where T0 is the orbit repeat period, and n0 and n c are the mean motion and the correction of the mean motion, respectively.

[0151] The daily advance time (T a ) is provided to S7, and the formula is as follows:

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

[0153] S7 step: The multipath error of the previous day is extracted by using the estimated orbit repeat period and wavelet transform, and the single-difference residual to be corrected is obtained.

[0154]

[0155] where, is the single-difference residual to be corrected by the multipath error of the second day, is the original single-difference residual of the second day, is the multipath error of the first day. The single-difference residual to be corrected by the multipath error is provided to S8.

[0156] S8 step: The single-difference residual to be corrected by the multipath error is substituted into equation (11) to obtain the double-difference residual to be corrected by the multipath error.

[0157] S9 step: The Kalman filter is used to estimate the float ambiguity and the float solution of the position, and finally the float ambiguity is fixed to obtain the precise position solution after the ambiguity is fixed.

[0158] The above merely describes the preferred embodiments of the present application and does not limit the present application in any way. Any person skilled in the art can make any form of equivalent replacement or modification to the technical solutions and technical contents disclosed by the present application without departing from the scope of the technical solutions of the present application, and such changes still belong to the protection scope of the present application.

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 GNSS real-time dynamic positioning (RTK) mode, the multipath effect must generally be considered. Therefore, firstly, the observation data of the base station and the rover are preprocessed to remove data with unqualified observation quality. Then, the qualified observation data is used for dynamic positioning processing to calculate their corresponding post-hoc double-difference residuals. The double-difference residuals are considered to contain only multipath error 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; use the improved wavelet threshold function and heuristic threshold to correct the wavelet coefficients; reconstruct the corrected wavelet coefficients to generate a noise-free signal, thereby obtaining high-precision multipath error and achieving the purpose and effect of multipath extraction. (4) Step L4: Calculate the repetition period of the track using the track repetition time method; then, based on the track repetition period, use the single-difference residual to deduct the multipath error of the previous day; finally, convert the single-difference residual after multipath correction into double-difference observations to estimate the precise location of the rover station.

2. The method as described in claim 1, characterized in that, Quality control in the pretreatment stage of step L1 includes: Preprocessing of observation data mainly includes pseudorange single-point positioning, pseudorange observation value verification, and correlation analysis for gross error detection of pseudorange data, as well as cycle slip detection of phase data using MW combination and GF combination. (a) Pseudorange observation check refers to calculating the distance from the satellite to the station, i.e., the satellite-to-ground distance, using the known coordinates of the reference station and the satellite coordinates obtained from the ephemeris. The pseudorange observation corresponding to the satellite-to-ground distance is then compared. If the pseudorange observation differs significantly from the satellite-to-ground distance, the observation is considered to contain gross errors and is discarded. Correlation analysis gross error detection refers to detecting potential gross errors in the observations based on the gross error theory of correlation analysis after obtaining the residuals of the pseudorange observations through single-point positioning. If a gross error is detected, the observation is discarded. (b) By linearly combining observations of different frequencies or types, some combined observations with certain excellent characteristics can be obtained. These combined observations can be used efficiently to detect cycle slips. The calculation formulas for MW and GF combinations are as follows: (1a) (1b) in, and These are the observed values ​​after combining MW and GF, respectively; and These are carrier phase observations at frequencies L1 and L2, respectively. and These are pseudorange observations in the L1 and L2 frequency bands, respectively; and These are the L1 and L2 frequencies, respectively. If the difference between combined observations in adjacent epochs is greater than a threshold, where the threshold for the GF combination is 0.25 cycles and the threshold for the MW combination is 0.5 cycles, an ambiguity flag is added. This means that one ambiguity arc ends and another new ambiguity arc begins. The new ambiguity arc will end when satellite tracking ends or when another cycle jump epoch is encountered, and this process repeats until the ambiguity of all data arcs is flagged.

3. The method as described in claim 1, characterized in that, The double-difference residual calculation in step L1 is as follows: Double-difference observations can be used to subtract common error terms such as ionospheric error, tropospheric error, satellite orbit error, and satellite clock error, thereby enabling more accurate calculation of integer ambiguity and baseline vector, and ultimately obtaining the posterior double-difference residuals. Assume the original pseudorange and carrier phase observation equations are as follows: （2a） (2b) Among them, subscript , and These represent the station, pseudorange, and carrier phase, respectively; superscript Represents a specific satellite; and These are the original carrier phase and pseudorange observations; and These are the receiver and the satellite clock bias, respectively. and These represent the speed of light in a vacuum and the carrier phase wavelength, respectively. and These are the tropospheric and ionospheric delays, respectively. and These are the satellite-to-ground distance and phase integer ambiguity, respectively. and These are the hardware delays for pseudorange and phase, respectively; and These represent the pseudorange and carrier phase multipath errors, respectively. and These represent the noise of the pseudorange and the carrier phase, respectively. The ionospheric and tropospheric errors are eliminated by combining observations using a double-difference method, and the calculation formula is as follows: (3a) (3b) Among them, superscript Indicates the selected reference satellite, indicated by the superscript. Represents any other satellite; subscript Indicates the base station, subscript Indicates a mobile station; The double difference operator representing the difference between stations and between satellites; The baseline used has been iteratively processed, geometric distance The theoretical geometric distance can be calculated from satellite orbit parameters and station coordinates. Instead, equation (3a) is further expressed as follows: (4a) (4b) (a) First, the floating-point ambiguity is estimated using Kalman filtering. and The corresponding standard deviation is and ; (b) Sort all ambiguities according to their standard deviations and start unambiguity resolution from the smallest; if the floating-point ambiguity and its standard deviation satisfy equations (5a), 5(b), (6a) and (6b), round the floating-point ambiguity to the fixed ambiguity. (5a) (5b) (6a) (6b) in: (7a) (7b) Take 25 as the threshold. Take 0.25 weeks as the threshold; (c) After searching and fixing, the fixed fuzziness is used as the weight. The virtual observations are used to update the Kalman filter, and the floating-point ambiguities with the smallest standard deviation in the list are searched and fixed until all ambiguities are fixed. (d) Use Kalman filtering with fixed ambiguity to estimate unfixed ambiguity, and search and fix floating-point ambiguity again; this process will be iterated until all ambiguities are fixed, or no more ambiguities can be fixed. (e) The data with fixed ambiguity is solved using the Kalman filter with Integer Rounding (IR) method to obtain the fixed solution for the station; finally, the double-difference residual can be expressed as: (8a) (8b) The pseudorange and carrier phase double difference residuals are obtained.

4. The method as described in claim 1, characterized in that, The calculation for converting the double-difference residual into a single-difference residual in step L2 is as follows: First, a zero-mean constraint is applied to the single-difference residual to ensure that the coefficient matrix between the single-difference and double-difference residuals is full rank. Then, the double-difference residual containing multipath and observation noise is transformed to obtain the single-difference residual without receiver clock bias. The calculation formula is as follows: （9） （10） （11） Where, assuming As the reference star, For single-difference residuals, As a weighting factor, The elevation angle; Then, perform the inverse operation on the coefficient matrix to obtain the single-difference residuals: （12） (13a) (13b) The pseudorange and carrier phase single-difference residuals are obtained.

5. The method as described in claim 1, characterized in that, The wavelet transform in step L3 decomposes the single-difference residual, and the wavelet coefficients are calculated as follows: First, the wavelet basis functions are discretized. Then, the discretized wavelet basis functions and the single-difference residuals are used to perform an inner product operation to obtain the coefficients after wavelet decomposition. The calculation formula is as follows: （14） in The single-difference residual signal that needs to be transformed For wavelet basis functions, It is the expansion factor. These are translation parameters; Use binary and The wavelet basis functions are discretized as follows: （15） in, and They are all integers.

6. The method as described in claim 1, characterized in that, In step L3, the decomposed wavelet coefficients are corrected using an improved wavelet threshold function and a heuristic threshold, calculated as follows: The improved wavelet thresholding function solves the problem of discontinuity of hard threshold at a pre-set threshold and the problem of fixed difference between the modified wavelet coefficient and the original wavelet coefficient in soft threshold; the heuristic threshold is a combination of general threshold and unbiased risk estimation threshold. The improved threshold function is calculated using the following formula: （16） in, To quantize wavelet coefficients, Threshold parameter, For exponentiation operator; The heuristic threshold calculation formula is as follows: （17） （18） （19） in and For the calculated intermediate variable values, Indicates the number of decomposition layers of the signal. The length of the single-difference residual signal; if If a general threshold is selected as the threshold for the wavelet, then the smaller of the general threshold and the unbiased risk estimation threshold is selected as the threshold for the wavelet. The formula for calculating the general threshold is as follows: （20） （21） in, The standard deviation of noise. The median absolute deviation of the wavelet coefficients; The formula for calculating the unbiased risk estimation threshold is as follows: （22） （23） （24） in, For the first A risk vector, , It is a positive integer; The squares of the wavelet coefficients are sorted in ascending order. This corresponds to the minimum value of the risk vector.

7. The method as described in claim 1, characterized in that, The calculation of the orbital repetition period estimation in step L4 is as follows: First, 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, Kepler's third law is used to estimate the period of the satellite's orbit. Then, based on the correction parameter of the satellite's average motion in the broadcast, the orbital repetition period of the satellite is estimated. Finally, the difference between the orbital repetition period and the International Day is used to obtain the daily advance time. Using Kepler's third law, the formula for the orbital period is as follows: （25） in, For orbital period, And the semi-major axis of the satellite orbit, The gravitational constant of Earth; After obtaining the orbital period, perform a multiplication operation on it to obtain the orbital repetition period, as shown in the following formula: （26） in, It is the orbital repetition period. and These are the mean motion and the correction to the mean motion, respectively. Daily advance time The formula is as follows: （27） in, It's the time to arrive earlier each day.

8. The method as described in claim 1, characterized in that, In step L4, the multipath error extracted from the previous day is subtracted from the single-difference residual. The calculation formula is as follows: The estimated orbit repetition period and wavelet transform are used to extract the multipath error of the previous day, which is then used to correct the single-difference residual. (28) in, This is the single-difference residual after multipath error correction on the second day. This is the original single-difference residual for the second day. This represents the multipath error for the first day.

9. The method as described in claim 4, characterized in that, Step L4, which converts the corrected single-difference residuals into double-difference observations and estimates the precise location of the rover, includes: First, substitute the single-difference residual after multipath error correction into equation (11) to obtain the double-difference residual after multipath error correction; then, use Kalman filtering to estimate the floating-point ambiguity and the floating-point solution of the position; finally, fix the floating-point ambiguity to obtain the precise position solution after ambiguity fixation.

Images