Multi-path error weakening method for multi-frequency multi-mode signal

Abstract

The invention discloses a multi-frequency multi-mode signal multi-path error weakening method, and relates to the technical field of precision positioning, and the method comprises the steps: constructing a multi-frequency multi-mode non-combination precision single-point positioning observation equation, separating multi-path errors in pseudo-range and phase observation values through the equation, and forming residual data; performing quality control on the residual error, and eliminating an abnormal value; feature extraction and judgment are carried out based on residual data, wherein time repeatability features, spatial distribution features and trend complexity features are extracted; selecting a correction model according to a feature judgment result, wherein the correction model comprises a fixed star daily filtering model, a multi-path hemispherical graph model and a trend surface analysis multi-path hemispherical graph model; and generating a multi-path error correction amount through a selected model, correcting an original observation value, and outputting a positioning result. The method solves the problem that the positioning precision is reduced due to multi-path errors, can effectively weaken the multi-path errors, improves the positioning precision and stability, adapts to a complex environment, and reduces data processing redundancy.

Description

Technical Field

[0001] The present invention relates to the field of precision positioning technology, and in particular to a method for reducing multi-path errors of multi-frequency and multi-mode signals. Background Art

[0002] In GNSS precision positioning technology, multipath error is one of the main factors affecting positioning accuracy. This error arises from the superposition of interference between reflected waves and direct waves generated when satellite signals encounter buildings, the ground, or other reflectors during propagation, leading to deviations in receiver observations (such as pseudorange and phase). Especially in Precise Point Positioning (PPP) applications, such as surveying and mapping, geological disaster monitoring, and autonomous driving, positioning results have strict requirements for millimeter to centimeter-level accuracy. The cumulative effect of multipath error can significantly reduce positioning reliability. For example, in urban canyons, forests, or environments with significant multipath effects, pseudorange observations may produce decimeter-level deviations, while phase observations may introduce cycle slips or ambiguity resolution errors, which in turn affect the convergence speed and stability of the final position solution.

[0003] Existing methods for mitigating multipath errors primarily include sidereal filtering (SF) and multipath hemispheric mapping (MHM) models. The SF model exploits the periodicity of satellite orbits to correct errors through temporal repeatability, but its effectiveness is limited in dynamic environments (such as mobile receivers or changing reflectors) because the errors may exhibit non-repeatable characteristics. The MHM model, based on the statistical residual distribution of spatial grids, is suitable for static scenarios. However, in low-elevation areas or when signal obstruction is frequent, the sparse grid data leads to inaccurate corrections. Furthermore, existing methods are mostly targeted at a single frequency band or system (e.g., the GPS L1 band only), while modern GNSS has evolved to multi-frequency and multi-mode integration (e.g., GPS, GLONASS, Galileo, and BDS-3 in parallel). Signal frequency band differences (e.g., hardware delay bias) are not fully compensated, resulting in incomplete error separation. Existing technologies also lack adaptive mechanisms: Model selection is fixed and cannot be dynamically switched based on error characteristics (such as temporal repeatability, spatial distribution, or trend complexity). This reduces correction effectiveness in complex multipath scenarios, with residual variance reductions potentially limited to 10-15%. Positioning Dilution of Precision (PDOP) can easily exceed thresholds, impacting practicality. Therefore, a method is urgently needed that can adapt to multi-frequency and multi-mode signals, dynamically identify the dominant error characteristics, and optimize model selection to improve the robustness and accuracy of multipath error mitigation. Summary of the Invention

[0004] The present invention overcomes the problem in the prior art of multipath errors leading to decreased positioning accuracy, and provides a method for reducing multipath errors of multi-frequency and multi-mode signals, which can effectively reduce multipath errors, improve positioning accuracy and stability, adapt to complex environments, and reduce data processing redundancy.

[0005] In order to achieve the above object, the present invention adopts the following scheme: A multi-path error reduction method for multi-frequency and multi-mode signals comprises the following steps: Construct a multi-frequency, multi-mode, non-combined precise point positioning observation equation. The observation equation includes a combination of multi-frequency signals from multiple global navigation satellite systems. The observation equation is used to separate the multipath errors in pseudorange and phase observations. Quality control is performed on the separated pseudorange and phase residuals to eliminate abnormal residual data. Feature extraction and judgment of residual data after quality control, including: The temporal repeatability characteristics of the residual data were extracted by calculating the autocorrelation coefficient of the residual sequence. When the autocorrelation coefficient was greater than 0.7, it was determined to be dominated by temporal repeatability. The spatial distribution characteristics of the residual data were extracted by statistically analyzing the residual density in the elevation and azimuth grids. When the residual density in the grid with an elevation angle greater than 60° was less than 50% of the residual density in the grid with an elevation angle less than 30°, it was determined to be dominated by spatial distribution. The residual trend was fitted by a cubic polynomial and the goodness of fit R was calculated. 2 The value and F test value extract the trend complexity characteristics. When R 2 When the value is greater than or equal to 0.3 and the F test is significant at the 95% confidence level, it is determined to be dominated by trend complexity; The multipath error correction model is selected based on the feature judgment results. If it is determined that time repeatability is dominant, the sidereal day filter SF model is selected; if it is determined that spatial distribution is dominant, the multipath hemispherical map MHM model is selected; if it is determined that trend complexity is dominant, the trend surface analysis multipath hemispherical map TMHM model is selected; The multipath error correction of pseudorange observations and phase observations is generated through the selected correction model, the original observations are corrected and the positioning results are output.

[0006] Preferably, the specific steps of constructing the multi-frequency multi-mode non-combined precise point positioning observation equation include: The multi-frequency signal combination of multiple global navigation satellite systems includes GPS three-frequency, GLONASS two-frequency, Galileo five-frequency and BDS-3 five-frequency signals. The basic observation equations of pseudorange and phase are: in, and are the original pseudorange observation value and phase observation value of the receiver r for the i-th frequency point of satellite s, is the geometric distance, and are the receiver and satellite clock errors, is the mapping function, T r is the zenith tropospheric delay; , is the i-th frequency value of satellite s, is the first frequency value, is the slant ionospheric delay of the first frequency; 、 、 、 are the hardware delay deviations between the receiver and the satellite for pseudorange and phase, Fuzziness The wavelength, and are pseudorange and phase multipath errors, respectively; and All are observation noise; When separating multipath errors, the inter-band deviation compensation is performed on the pseudorange and phase observations at different frequencies based on the frequency band differences of the satellite and receiver hardware delay deviations; and the multi-frequency observations are weightedly fused according to the signal-to-noise ratio of each frequency signal and the satellite elevation angle.

[0007] Preferably, the specific steps of performing quality control include: Based on the normal distribution assumption, the mean μ is calculated for the pseudorange residual and phase residual in each elevation and azimuth grid. I and standard deviation σ I , the 3σ rule is used to identify outliers, and the judgment rules are: ,in, represents the i-th residual value in grid I; For the identified outliers, the F test is used to verify their significant impact on the grid residual variance. The specific steps are as follows: the residual data set in the grid is divided into the normal data set x I,1 and the abnormal data set x I,2 , and calculate the sample variance of the two data sets and , construct the F test statistic: If F is greater than the parameter α, take F of 0.05 α (n I,1 −1,n I,2 −1) value, it is determined that the outliers have a significant impact on the residual variance; For the residual data in the low elevation angle area with an elevation angle less than 15°, the abnormality judgment threshold is adjusted to: , where the satellite motion state dynamic factor , E s is the satellite altitude angle; For the grids that have eliminated outliers, the sliding window interpolation method of the residual mean of adjacent grids is used to fill the data gaps. The window size is 3×3 grids, and the interpolation formula is: , where N adj is the number of adjacent grids, μ I,jis the mean of residuals of adjacent grids; For the corrected residual data, the residual distribution is statistically analyzed by satellite system partition. If the residual variance of one system exceeds twice the global variance, the secondary quality control of the data of that system is triggered.

[0008] Preferably, the specific steps of extracting the time repeatability feature include: The autocorrelation coefficients of pseudorange residual and phase residual sequences are calculated using sliding windows. The window length is set to an integer multiple of the satellite orbit repetition period, and the sliding step is 1 / 10 of the data sampling interval. A hypothesis test is performed on the autocorrelation coefficient of each window. If the autocorrelation coefficient is greater than 0.7 within the 95% confidence interval, the residual in the window is judged to have the dominant characteristic of time repeatability; According to the wavelength difference of signals in different frequency bands, the autocorrelation coefficients are subjected to frequency band weighted fusion, and the weight w i The calculation formula is: Among them, λ i is the wavelength of the signal in the i-th frequency band, n sys is the total number of frequency bands of the current satellite system; If the same satellite is judged to be time-repeatability dominated in multiple consecutive windows and the range of its orbital elevation angle variation is less than 10°, the residual sequence of the satellite is marked as time-repeatability dominated as a whole.

[0009] Preferably, the specific operation mode of the sidereal day filter model includes: Calculate the orbit repetition period T based on the satellite ephemeris data orbs , and the residual sequence is T orbs Segment alignment, controlling the alignment error to no more than 5% of the data sampling interval; The mean of the pseudorange and phase residuals within each orbital period is calculated, and the abnormal segments that deviate from the mean by more than 2 standard deviations are eliminated. The inter-band deviation of the mean sequence is corrected for the hardware delay bias of the multi-band signal. The mean sequence is smoothed using a sliding window with a length of 3 orbital periods. The weights within the window are dynamically allocated according to the satellite elevation angle. The weight allocation formula is: ,in, is the average altitude angle of the satellite in the kth orbital period; The smoothed mean sequence is used as the correction value of the sidereal day filter model, matched to the original observation value by timestamp, and the corrected pseudorange and phase observation values ​​are output.

[0010] Preferably, the specific steps of extracting spatial distribution features include: Dynamically divide the grid and dynamically adjust the grid's azimuth resolution according to the satellite elevation angle range. The specific adjustment rules are as follows: , where ΔAk is the adjusted azimuth resolution, ΔA0 is the reference azimuth resolution, and E k is the center value of the elevation angle interval of the kth layer, and the elevation angle range is divided into 5° intervals; the pseudorange residuals and phase residuals in each grid are statistically analyzed for density ρ I , the density calculation formula is: , where N I is the number of effective residuals in grid I, A I is the spherical area of ​​the grid; For high elevation angle areas with elevation angles greater than 60°, the residual density is corrected to , density compensation factor , E s is the satellite elevation angle; if the corrected average residual density of the high-elevation area grid is less than 50% of the average residual density of the low-elevation area grid with an elevation angle of less than 30°, it is determined to be spatially dominant; The residual density of all frequency bands of the same satellite system is jointly analyzed. If more than 70% of the frequency bands meet the spatial distribution dominance condition, the system is globally marked as spatial distribution dominance.

[0011] Preferably, the specific operation mode of the multipath hemispherical graph model includes: For each dynamically divided grid, the mean of pseudorange residuals and phase residuals is calculated respectively. For multi-band signals, the residual means are weighted and fused according to the wavelength of the frequency band. For grids without valid data, the residual means of the adjacent three layers of elevation grids are used for interpolation, and the interpolated grid data are spatially smoothed using a Gaussian kernel function with a standard deviation of 2°. When generating the correction value, the grid residual mean is dynamically weighted according to the real-time satellite elevation angle to enhance the correction contribution of the low elevation angle area; the weighted grid residual mean is mapped to the original observation value according to the real-time satellite position, and the corrected pseudorange and phase observation values ​​are output.

[0012] Preferably, the specific steps of extracting trend complexity features include: For the residual data in each elevation and azimuth grid, linear, quadratic and cubic polynomials are used to perform trend fitting in turn, and the optimal order AIC value is selected by the Akaike Information Criterion. , where k is the number of model parameters, L is the model likelihood value, and the polynomial order with the smallest AIC value is selected as the final fitting model; Calculates R for a polynomial fit of a selected order 2 The trend significance is verified by F test, and the trend complexity feature judgment rule is: R 2 ≥0.3, the F test statistic is significant at the 95% confidence level; For multi-band signals, the residual trend is weighted and fused according to the wavelength difference; if the change of the trend fitting coefficient of the same satellite in more than 5 consecutive grids is less than 10%, the residual trend of the satellite is globally smoothed; when R is satisfied, 2 When ≥0.3 and the F test is significant, it is determined to be dominated by trend complexity; if multiple order models meet the conditions at the same time, the highest order model is selected.

[0013] As a preferred embodiment, the specific operation mode of the trend surface analysis multi-path hemispherical graph model includes: For the residual data in each grid, the trend surface is fitted, and the model expressions of linear, quadratic and cubic trend surfaces are: in, is the multipath residual, A s is the azimuth, E s is the altitude angle, c0 to c9 are the polynomial coefficients of the trend surface; In view of the hardware delay deviation of multi-band signals, the trend surface coefficients are compensated for frequency band differences; the trend surface coefficients are dynamically weighted according to the real-time motion status of the satellite, and for grids without valid data, bilinear interpolation of the trend surface coefficients of adjacent grids is used to fill the gaps; the optimized trend surface coefficients are mapped to the observation values ​​according to the real-time satellite position, and the corrected pseudorange and phase data are output.

[0014] Preferably, when the residual data simultaneously satisfies multiple judgment conditions of temporal repeatability dominance, spatial distribution dominance, or trend complexity dominance, the multipath error correction model is selected according to the following priority rules: The model selection priority is defined as: temporal repeatability dominant > trend complexity dominant > spatial distribution dominant; If the residual data satisfies both the temporal repeatability-dominant and other feature-dominant conditions, the sidereal day filter model is preferred; if the residual data satisfies both the trend complexity-dominant and spatial distribution-dominant conditions, the trend surface analysis multipath hemispherical graph model is preferred. When multiple feature-dominant conditions coexist, calculate the residual variance contribution corresponding to each feature, and the contribution ,in, is the residual variance corresponding to feature f, is the total residual variance; The model corresponding to the feature with the highest contribution is selected. If the difference in contribution is less than 10%, the combined model correction mechanism is enabled: Perform joint correction on multiple feature models with similar contribution, and the comprehensive correction amount Δ comb The calculation formula is: , where Δ SF is the correction amount of the SF model, Δ TMHMis the correction amount of the TMHM model, Δ MHM is the correction amount of the MHM model, ω1, ω2, ω3 are the contribution C of each model f Proportional weight distribution, and satisfy ; During the correction process, the precision reduction factor (PDOP) and residual variance of the positioning results are monitored in real time. If the PDOP value exceeds 3 or the residual variance decreases by less than 20% compared with the pre-correction value, the system dynamically switches to the model with the second highest contribution. For conflict scenarios that trigger model switching multiple times, the timestamp, satellite system, feature contribution, and positioning error data are recorded.

[0015] The present invention has at least the following beneficial effects: (1) through dynamic feature extraction and adaptive model selection, multipath error is effectively weakened, the accuracy and stability of precise single-point positioning are improved, and the positioning residual variance can be significantly reduced in complex environments; (2) by utilizing multi-frequency signal combination, the error separation accuracy is improved through inter-band deviation compensation and weighted fusion based on signal-to-noise ratio and elevation angle, and the multi-system compatibility is enhanced; (3) through quality control steps including outlier identification, F test verification, dynamic threshold adjustment and data interpolation, the reliability of residual data is ensured, abnormal interference is reduced, and the accuracy of subsequent feature extraction is improved; (4) the feature extraction mechanism covers the diversity of multipath errors, and combined with the correction model, the robustness of the system in a changing environment is improved; (5) through priority rules and combined model mechanisms, feature conflicts are handled and the correction effect is optimized; PDOP and residual variance are monitored in real time and the model is switched dynamically to enhance fault tolerance and avoid positioning failure. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 This is a principle flow chart of the multi-frequency and multi-mode signal multi-path error reduction method provided by the present invention. DETAILED DESCRIPTION

[0017] The present invention will be described in further detail below in conjunction with the accompanying drawings so that those skilled in the art can implement the invention with reference to the description.

[0018] like Figure 1 As shown, the multi-path error reduction method for multi-frequency multi-mode signals provided by the present invention includes the following steps: Construct a multi-frequency, multi-mode, non-combined precise point positioning observation equation. The observation equation includes a combination of multi-frequency signals from multiple global navigation satellite systems. The observation equation is used to separate the multipath errors in pseudorange and phase observations. Quality control is performed on the separated pseudorange and phase residuals to eliminate abnormal residual data. Feature extraction and judgment of residual data after quality control, including: The temporal repeatability characteristics of the residual data were extracted by calculating the autocorrelation coefficient of the residual sequence. When the autocorrelation coefficient was greater than 0.7, it was determined to be dominated by temporal repeatability. The spatial distribution characteristics of the residual data were extracted by statistically analyzing the residual density in the elevation and azimuth grids. When the residual density in the grid with an elevation angle greater than 60° was less than 50% of the residual density in the grid with an elevation angle less than 30°, it was determined to be dominated by spatial distribution. The residual trend was fitted by a cubic polynomial and the goodness of fit R was calculated. 2 The value and F test value extract the trend complexity characteristics. When R 2 When the value is greater than or equal to 0.3 and the F test is significant at the 95% confidence level, it is determined to be dominated by trend complexity; The multipath error correction model is selected based on the feature judgment results. If it is determined that time repeatability is dominant, the sidereal day filter SF model is selected; if it is determined that spatial distribution is dominant, the multipath hemispherical map MHM model is selected; if it is determined that trend complexity is dominant, the trend surface analysis multipath hemispherical map TMHM model is selected; The multipath error correction of pseudorange observations and phase observations is generated through the selected correction model, the original observations are corrected and the positioning results are output.

[0019] Multi-frequency, multi-mode, non-combined precise point positioning observation equations are constructed to separate multipath errors. This involves constructing the observation equations and separating multipath errors from pseudorange and phase observations. By using multi-frequency signal combinations from the Global Navigation Satellite System (GNSS) (such as GPS triple-frequency, GLONASS dual-frequency, Galileo quintuple-frequency, and BDS-3 quintuple-frequency), mathematical equations are used to describe the signal propagation process. This allows multipath errors to be separated from other system errors (such as clock bias and ionospheric delay) to form pseudorange and phase residuals. These residuals, which incorporate multipath effects and random noise, serve as input data for subsequent processing.

[0020] The observation equation is based on the uncombined precise point positioning (PPP) model, avoiding the use of traditional differential techniques and thus preserving independent information for each frequency band. The equation accounts for parameters such as the geometric distance between the satellite and the receiver, clock error, zenithal tropospheric delay, slant ionospheric delay, and hardware delay bias. Actual observations (pseudoranges and phases) are substituted into the equation to calculate the residual sequence. First, multi-system multi-frequency signal data is received. Second, the equation is applied to compensate for inter-band biases in the observations at each frequency point (for example, adjusting signal weights based on hardware delay differences). Finally, the residual data is separated. When these technical features are combined, the multi-frequency signal combination provides a richer dimension for error separation. The uncombined PPP model ensures that the residual directly reflects multipath effects, laying the data foundation for subsequent quality control.

[0021] Quality control is performed on the separated pseudorange and phase residuals to ensure data reliability and improve the accuracy of subsequent feature extraction by eliminating anomalous residuals. Quality control involves identifying, verifying, and eliminating anomalous data. Based on the assumption that the residual data follows a normal distribution, statistical methods (such as the 3σ rule) are used to identify outliers, and the significance of the outliers is verified using an F-test. In low-elevation-angle regions (e.g., elevation angles less than 15°), satellite signals are susceptible to dynamic interference, necessitating dynamic threshold adjustment (e.g., by introducing a satellite motion state factor that varies based on the satellite elevation angle).

[0022] Feature extraction and judgment are key to multipath model selection. These features are divided into three sub-features: temporal repeatability, spatial distribution, and trend complexity. Each sub-feature corresponds to a decision logic used to identify the dominant characteristics of multipath errors.

[0023] Temporal repeatability is extracted by calculating the autocorrelation coefficient (ACF) of the residual sequence. The periodicity of satellite orbits causes multipath errors to repeat in time (e.g., the orbital repetition period on which the sidereal day filter depends). When the autocorrelation coefficient is greater than 0.7 (the reference threshold), temporal repeatability is considered dominant, indicating that the error is characterized by a predominantly periodic pattern.

[0024] Spatial distribution feature extraction is performed by statistically analyzing the residual density within an elevation-azimuth grid. Multipath errors are affected by the position of reflectors and exhibit a non-uniform spatial distribution (for example, near-ground reflections lead to concentrated errors at low elevation angles). When the residual density for grids with elevation angles greater than 60° is less than 50% (the reference threshold) of the density for grids with elevation angles less than 30°, the spatial distribution is considered dominant, indicating that the error is primarily azimuth-dependent.

[0025] The trend complexity feature extraction is to fit the residual trend by a cubic polynomial and calculate the goodness of fit R 2 and F-test value. Multipath error may show nonlinear trend (such as complex changes caused by the movement of reflectors). 2 When the value is ≥ 0.3 (reference threshold) and the F test is significant at the 95% confidence level, the trend complexity is determined to be dominant, indicating that the error is mainly due to high-order changes. This includes fitting linear, quadratic, and cubic polynomials in sequence; selecting the optimal order (e.g., based on the AIC criterion); and calculating R. 2 The F-value verifies significance. When the features are combined, polynomial fitting captures nonlinear patterns, statistical tests ensure reliability, and multi-frequency weighting optimizes trend extraction.

[0026] During the overall feature extraction phase, three sub-features are run in parallel: temporal repeatability focuses on periodicity, spatial distribution focuses on azimuth dependence, and trend complexity handles nonlinear variations. Together, they independently determine the diversity of multipath errors and provide a basis for model selection.

[0027] A multipath error correction model is selected based on feature judgment results, and the optimal correction model is dynamically selected based on feature judgment. Feature-dominant results are mapped to corresponding models: the sidereal filter (SF) model is selected when temporal repeatability is dominant, the multipath hemispheric map (MHM) model is selected when spatial distribution is dominant, and the trend-based multipath hemispheric map (TMHM) model is selected when trend complexity is dominant. Each model targets specific error characteristics (e.g., SF leverages temporal repeatability, MHM leverages spatial distribution), and model selection ensures that the correction accurately matches the error pattern. The specific process includes: receiving feature judgment results and resolving conflicts according to a priority rule (e.g., temporal dominance > trend dominance > spatial dominance). If multiple features coexist, the contribution to the residual variance is calculated and the one with the highest contribution is selected. When combining features, priority rules are used to resolve judgment overlaps, contribution calculations quantify feature influences, and a multi-model combination mechanism enhances robustness. Model selection dynamically optimizes the correction strategy to avoid the limitations of a single model.

[0028] Corrections are generated and corrected using a selected correction model. The selected model is then applied to generate multipath error corrections, and the original observations are corrected to output a positioning result. This involves generating corrections, correcting observations, and outputting a positioning result. Models (SF / MHM / TMHM) are used to invert multipath error from residual data and apply this error as a negative offset to the original pseudorange and phase observations, mitigating the error. Finally, a high-precision position is calculated using the PPP algorithm. The specific process involves running the selected model (e.g., the SF model calculates time-averaged residuals, the MHM model uses grid-based mean values, and the TMHM model fits trend surface coefficients); generating a correction table (to store model outputs), matching the model outputs to the real-time satellite positions (altitude and azimuth) to obtain corrections; correcting the original observations (pseudorange and phase), and solving the uncombined PPP equation to output a positioning result. When all features are combined, the model generates real-time corrections, position matching ensures dynamic adaptability, and the correction steps are directly integrated into the positioning process, achieving end-to-end error mitigation.

[0029] This method achieves multipath error mitigation through a systematic process. By dynamically selecting the optimal correction model, it effectively suppresses the impact of multipath error on pseudorange and phase observations, thereby improving the stability and reliability of PPP positioning. The feature extraction judgment mechanism (temporal, spatial, and trend) accounts for the diversity of multipath effects, making the method applicable to different observation scenarios (such as urban canyons and open areas) and enhancing system robustness. The integration of quality control and feature extraction reduces redundant data processing, and the priority and combination mechanism of model selection avoids overfitting, ensuring efficient operation of the method in real-time or post-processing applications. Multi-frequency and multi-mode signal processing integrates the advantages of GPS, GLONASS, Galileo, and BDS-3 systems, comprehensively improving global navigation performance through frequency band weighting and bias compensation. Overall, this method achieves high practicality and versatility in the field of multipath error mitigation. It automatically identifies and selects the most appropriate model for multipath error correction based on acquired positioning data. This efficient and accurate method avoids the impact of human judgment errors on positioning accuracy, providing reliable support for high-precision positioning.

[0030] In another technical solution, the specific steps of constructing the multi-frequency multi-mode non-combined precise point positioning observation equation include: The multi-frequency signal combination of multiple global navigation satellite systems includes GPS three-frequency, GLONASS two-frequency, Galileo five-frequency and BDS-3 five-frequency signals. The basic observation equations of pseudorange and phase are: in, and are the original pseudorange observation value and phase observation value of the receiver r for the i-th frequency point of satellite s, is the geometric distance, and are the receiver and satellite clock errors, is the mapping function, T r is the zenith tropospheric delay; , is the i-th frequency value of satellite s, is the first frequency value, is the slant ionospheric delay of the first frequency; 、 、 、 are the hardware delay deviations between the receiver and the satellite for pseudorange and phase, Fuzziness The wavelength, and are pseudorange and phase multipath errors, respectively; and All are observation noise; When separating multipath errors, the inter-band deviation compensation is performed on the pseudorange and phase observations at different frequencies based on the frequency band differences of the satellite and receiver hardware delay deviations; and the multi-frequency observations are weightedly fused according to the signal-to-noise ratio of each frequency signal and the satellite elevation angle.

[0031] The multi-frequency signal combination covers GPS tri-frequency, GLONASS bi-frequency, Galileo quintuple-frequency, and BDS-3 quintuple-frequency signals, enhancing observation redundancy by integrating multi-system, multi-frequency data. Signals in different frequency bands differ in their sensitivity to multipath errors. For example, high-frequency signals are more susceptible to reflectors, while low-frequency signals are more significantly affected by ionospheric delay. The multi-frequency combination provides a complementary basis for error separation. In the equation construction, the geometric range parameters include corrections such as antenna phase center offset and relativistic effects, which require precise calculation using precise ephemeris and station coordinates. Clock error parameters are divided into receiver clock error and satellite clock error. Based on the IGS protocol, pseudorange hardware delay bias is absorbed to ensure synchronous elimination of clock errors and hardware bias. Zenith tropospheric delay is converted to slant-path delay using a mapping function, while slant ionospheric delay is linked to each frequency point using coefficients to absorb differential code bias between frequency bands.

[0032] Hardware delay bias occurs on both the receiver and satellite sides. Different frequency bands (such as the B1I and B2B bands of the BDS-3) may experience nanosecond-level variations, which are compensated for using pre-calibrated parameters or global averages. Weighted fusion of observations depends on the signal-to-noise ratio and satellite elevation angle: when the signal-to-noise ratio falls below a certain threshold (e.g., 30-35 dBHz), the weight is proportionally reduced. When the elevation angle falls below 15°, the weight is further reduced to 30%-60% of the original value to mitigate weak signal interference. During operation, the original observations are input and the precise ephemeris is used to calculate the geometric distance. Clock errors and hardware bias parameters are then substituted. Finally, pseudorange and phase residuals are output through frequency band compensation and dynamic weighting. This solution improves error separation accuracy and enhances data reliability. Multi-frequency fusion and dynamic weighting significantly enhance error separation accuracy and support multipath modeling in complex environments.

[0033] In another technical solution, the specific steps of performing quality control include: Based on the normal distribution assumption, the mean μ is calculated for the pseudorange residual and phase residual in each elevation and azimuth grid. I and standard deviation σ I , the 3σ rule is used to identify outliers, and the judgment rules are: ,in, represents the i-th residual value in grid I; For the identified outliers, the F test is used to verify their significant impact on the grid residual variance. The specific steps are as follows: the residual data set in the grid is divided into the normal data set x I,1 and the abnormal data set x I,2, and calculate the sample variance of the two data sets and , construct the F test statistic: If F is greater than the parameter α, take F of 0.05 α (n I,1 −1,n I,2 −1) value, it is determined that the outliers have a significant impact on the residual variance; For the residual data in the low elevation angle area with an elevation angle less than 15°, the abnormality judgment threshold is adjusted to: , where the satellite motion state dynamic factor , E s is the satellite altitude angle; For the grids that have eliminated outliers, the sliding window interpolation method of the residual mean of adjacent grids is used to fill the data gaps. The window size is 3×3 grids, and the interpolation formula is: , where N adj is the number of adjacent grids, μ I,j is the mean of residuals of adjacent grids; For the corrected residual data, the residual distribution is statistically analyzed by satellite system partition. If the residual variance of one system exceeds twice the global variance, the secondary quality control of the data of that system is triggered.

[0034] Based on the assumption of a normal distribution, the mean and standard deviation of the residuals within each elevation-azimuth grid are calculated. The 3σ rule is used to initially screen out outliers: any residual that deviates from the mean by more than three standard deviations is flagged as a suspected anomaly. For low-elevation areas (elevation angles <15°), where strong satellite motion increases signal fluctuations, a dynamic adjustment factor is introduced to expand the threshold range. This factor is negatively correlated with the elevation angle. For example, at an elevation angle of 10°, the threshold is relaxed to 3.5 standard deviations, and at an elevation angle of 5°, it is further relaxed to 4 standard deviations to avoid excessive exclusion of valid data.

[0035] Outliers are tested for significance using an F-test. The grid residuals are divided into a normal dataset and an outlier dataset, and the ratio of the sample variances between the two groups is calculated as the F-statistic. If this value exceeds the critical value of the F-distribution corresponding to the selected significance level (e.g., α = 0.05), the outlier is considered to have a significant impact on the overall variance and is removed. This removal process is performed iteratively: the residual with the largest deviation from the mean is removed each time, and the statistic is recalculated until the F-test is insignificant. This process ensures the scientific nature of outlier removal and avoids subjective misjudgment.

[0036] For grids without valid data, interpolation is performed using the mean residuals of adjacent 3×3 grids (for example, taking the average of the means of the eight surrounding grids). After interpolation, residual distribution statistics are calculated by satellite system. If the residual variance of a single system (such as GLONASS) exceeds twice the global variance, secondary quality control is triggered for that system: the anomaly detection and rejection process is re-executed. This mechanism effectively intercepts localized systemic anomalies and ensures overall data consistency. This quality control mechanism enhances data reliability. Its dynamic threshold and secondary quality control design significantly reduce false rejection rates, providing high-integrity residual data for subsequent feature extraction.

[0037] In another technical solution, the specific steps of extracting the time repeatability feature include: The autocorrelation coefficients of pseudorange residual and phase residual sequences are calculated using sliding windows. The window length is set to an integer multiple of the satellite orbit repetition period, and the sliding step is 1 / 10 of the data sampling interval. A hypothesis test is performed on the autocorrelation coefficient of each window. If the autocorrelation coefficient is greater than 0.7 within the 95% confidence interval, the residual in the window is judged to have the dominant characteristic of time repeatability; According to the wavelength difference of signals in different frequency bands, the autocorrelation coefficients are subjected to frequency band weighted fusion, and the weight w i The calculation formula is: Among them, λ i is the wavelength of the signal in the i-th frequency band, n sys is the total number of frequency bands of the current satellite system; If the same satellite is judged to be time-repeatability dominated in multiple consecutive windows and the range of its orbital elevation angle variation is less than 10°, the residual sequence of the satellite is marked as time-repeatability dominated as a whole.

[0038] The sliding window length is set to an integer multiple of the satellite orbit repetition period, with possible values ​​including 1, 2, or 3 times the period, depending on the data sampling density and orbital characteristics. The window sliding step is set to 1 / 10 of the data sampling interval, for example, 3 seconds for 30 seconds of sampling data, ensuring that the window overlap covers more than 90% of the data. The periodic operation of the satellite orbit causes multipath errors to show temporal repeatability. The strength of the periodicity is quantified by the autocorrelation coefficient (ACF) of the residual sequence within the window. During implementation, the ACF value is calculated for the pseudorange residual and phase residual in each window. If the ACF value is continuously above a threshold (for example, 0.6-0.8), it indicates that significant temporal repeatability exists in that period.

[0039] A confidence interval is constructed using a 95% confidence level. If the ACF value is continuously greater than 0.7 (reference threshold) within the confidence interval, the residual in the current window is determined to have a dominant characteristic of temporal repeatability. For multi-band signals (such as GPS L1 / L2 / L5), weights are assigned based on wavelength differences: signals with longer wavelengths (such as BDS-3 B3I, with a wavelength of approximately 25 cm) are given higher weights, while signals with shorter wavelengths (such as GPS L1, with a wavelength of approximately 19 cm) are given lower weights. In the weight calculation formula, wavelength λ i The sum of the system's total wavelengths is used as the numerator, ensuring that the long-wavelength signal dominates the decision. The runtime calculates the ACF value in the time window, performs hypothesis testing, and then fuses the multi-band ACF results according to the weights to output a time-dominant decision flag.

[0040] The overall satellite marking mechanism relies on continuous window determination and orbital stability. If the same satellite is determined to be time-dominated in multiple consecutive windows (for example, 3-5 windows) and its orbital elevation angle variation range is less than 10° (calculated through ephemeris data), all residual sequences of the satellite are uniformly marked as time-dominated. During implementation, the standard deviation of the satellite elevation angle needs to be monitored. If the standard deviation is less than 1.5° (for example, geostationary orbit satellites), the global marking is directly triggered. This design ensures full-time temporal modeling for highly stable satellites (such as geosynchronous satellites) to avoid discontinuities introduced by segmented processing. The temporal feature extraction mechanism significantly improves the recognition accuracy of periodic multipath errors, and the multi-band weighting strategy enhances the robustness of judgment in complex signal environments.

[0041] The specific operation mode of the sidereal day filter model includes: Calculate the orbit repetition period T based on the satellite ephemeris data orbs , and the residual sequence is T orbs Segment alignment, controlling the alignment error to no more than 5% of the data sampling interval; The mean of the pseudorange and phase residuals within each orbital period is calculated, and the abnormal segments that deviate from the mean by more than 2 standard deviations are eliminated. The inter-band deviation of the mean sequence is corrected for the hardware delay bias of the multi-band signal. The mean sequence is smoothed using a sliding window with a length of 3 orbital periods. The weights within the window are dynamically allocated according to the satellite elevation angle. The weight allocation formula is: ,in, is the average altitude angle of the satellite in the kth orbital period; The smoothed mean sequence is used as the correction value of the sidereal day filter model, matched to the original observation value by timestamp, and the corrected pseudorange and phase observation values ​​are output.

[0042] First, the orbital repetition period is accurately calculated based on satellite ephemeris data. For example, the GPS satellite period is approximately 86,164 seconds (23 hours, 56 minutes, and 4 seconds). The residual sequence is segmentally aligned according to this period, with the alignment error controlled to within 5% of the data sampling interval. For example, a 1.5-second time deviation is allowed for 30-second sampling. Linear interpolation is used to adjust the residual timestamps during implementation to ensure strict alignment of the start times of each period segment. This is based on the principle that the consistency of the orbital period determines the temporal repeatability of multipath errors, and accurate alignment maximizes the effectiveness of model correction.

[0043] The arithmetic mean of pseudorange and phase residuals is calculated. If a segment of the residual sequence deviates from the mean by more than two standard deviations (for example, a phase residual fluctuation exceeds 5 cm), it is marked as an abnormal segment and excluded. For multi-band signals (such as Galileo E1 / E5a / E5b), inter-band hardware delay offsets are compensated by correcting the mean sequence based on precalibrated parameters (for example, the E5a band offset is +0.05 m). This process includes segmenting the residuals by period, calculating the mean within the segment, removing the out-of-limit data segments, and correcting the mean according to the band compensation table.

[0044] The smoothing process uses a sliding window with a length of 3 orbital periods. The weight of each period in the window is dynamically allocated according to the average satellite elevation angle. In the weight calculation formula, the sine square value of the satellite elevation angle ( ) as a weighting factor, for example, a weight of approximately 0.25 at an altitude angle of 30° and approximately 0.75 at an altitude angle of 60°. During implementation, the mean sequences of the three periods within the window are weighted and averaged, outputting a smoothed correction sequence. Finally, the correction is matched to the original observations by timestamp, for example, applying the correction from period k to the corresponding observation at period k+1. Orbital period alignment and anomalous segment removal significantly improve model accuracy, while the dynamic weighted smoothing mechanism effectively suppresses random noise interference.

[0045] In another technical solution, the specific steps of extracting spatial distribution features include: Dynamically divide the grid and dynamically adjust the grid's azimuth resolution according to the satellite elevation angle range. The specific adjustment rules are as follows: , where ΔA k is the adjusted azimuth resolution, ΔA0 is the reference azimuth resolution, and E k is the center value of the elevation angle interval of the kth layer, and the elevation angle range is divided into 5° intervals; the pseudorange residuals and phase residuals in each grid are statistically analyzed for density ρ I , the density calculation formula is: , where N I is the number of effective residuals in grid I, A I is the spherical area of ​​the grid; For high elevation angle areas with elevation angles greater than 60°, the residual density is corrected to , density compensation factor , E s is the satellite elevation angle; if the corrected average residual density of the high-elevation area grid is less than 50% of the average residual density of the low-elevation area grid with an elevation angle of less than 30°, it is determined to be spatially dominant; The residual density of all frequency bands of the same satellite system is jointly analyzed. If more than 70% of the frequency bands meet the spatial distribution dominance condition, the system is globally marked as spatial distribution dominance.

[0046] The grid division dynamically adjusts the azimuth resolution based on the satellite elevation angle range. The specific rules are as follows: the base azimuth resolution ΔA0 is usually set to 5° (optional range 3°-10°), and the azimuth resolution ΔA of the kth layer elevation angle interval is set to 5°. t According to the formula ΔA t = ΔA0 / cosE t Calculate, where E t Take the center value of that layer (for example, the center value of the elevation angle range of 10°-15° is 12.5°). The elevation angle range is divided into fixed intervals, with a recommended interval of 5° (for example, 0°-5°, 5°-10°, and finally 85°-90°). As the satellite elevation angle increases, the signal coverage decreases. Dynamically expanding the azimuth resolution maintains approximately the same grid surface area and avoids data sparseness caused by too small grids in high elevation areas. Implementation requires traversing all elevation angle levels, calculating the azimuth resolution layer by layer, and generating the grid structure.

[0047] The residual density ρ of each grid I t Defined as the number of effective residuals N t The area of ​​the grid tennis surface is A t The ratio (ρ t = N t / A t For high elevation angle areas with elevation angles greater than 60°, a density compensation factor α = 1 + 0.5 × (1 - Eˢ / 90°) is introduced for correction (e.g. E s =70°). The condition for determining spatial dominance is that the corrected average density in high-elevation angle areas (>60°) is less than 50% of the average density in low-elevation angle areas (<30°). The operation process includes calculating the density of all grids, correcting the high-elevation angle density, calculating the average density for each region, and comparing the calculated average density with the threshold.

[0048] The residual density of all frequency bands within a satellite system (such as GPS or BDS-3) is calculated. If more than 70% (with an optional threshold of 60%-80%) of the frequency bands meet the spatial distribution dominance criteria, all data from that system is marked as spatially dominant. Implementation requires calculating the percentage of density compliance for each frequency band by system partition. Once a global flag is triggered, it is automatically applied to all frequency point data for that system. This design avoids misjudgments of single frequency bands and improves system-level modeling consistency. Dynamic gridding and density correction significantly enhance spatial feature characterization capabilities, and system-level joint analysis enhances the reliability of judgments in multi-band scenarios.

[0049] The specific operation of the multipath hemispheric graph model includes: For each dynamically divided grid, the mean of pseudorange residuals and phase residuals is calculated respectively. For multi-band signals, the residual means are weighted and fused according to the wavelength of the frequency band. For grids without valid data, the residual means of the adjacent three layers of elevation grids are used for interpolation, and the interpolated grid data are spatially smoothed using a Gaussian kernel function with a standard deviation of 2°. When generating the correction value, the grid residual mean is dynamically weighted according to the real-time satellite elevation angle to enhance the correction contribution of the low elevation angle area; the weighted grid residual mean is mapped to the original observation value according to the real-time satellite position, and the corrected pseudorange and phase observation values ​​are output.

[0050] For each dynamically divided grid, the arithmetic mean of the pseudorange residual and phase residual is calculated. For multi-band signals (such as Galileo five-band), the residual mean is weighted and fused according to the wavelength of the band: the band with longer wavelength (such as BDS-3B3I, wavelength 25.2cm) is given a higher weight, and the band with shorter wavelength (such as GPS L1, wavelength 19.0cm) is given a lower weight. In the weight distribution formula, the wavelength λ i The total wavelength of the system is used as the numerator and the total wavelength of the system is used as the denominator. During implementation, the residual mean is stored by frequency band and the weighted average is calculated as the final grid output value.

[0051] A layered interpolation strategy is used for grids with no valid data (e.g., due to lack of satellite coverage). For grids without valid data (e.g., due to lack of satellite coverage), the mean of the residuals from three adjacent elevation grid layers (e.g., the current layer ± 1 layer) is used for interpolation. For example, for a grid with an empty elevation of 20°-25°, the mean of the 15°-20° and 25°-30° layers is used for interpolation. After interpolation, spatial smoothing is performed using a Gaussian kernel function. The recommended kernel standard deviation is 2° (optional range: 1°-3°). The smoothing range covers the surrounding 3×3 grid area. The implementation process includes identifying the empty grid location, retrieving data from adjacent layers, calculating the interpolated value, and applying Gaussian kernel convolution for smoothing.

[0052] Dynamically adjust the grid residual contribution weight according to the real-time satellite elevation angle: elevation angle E s When the angle is less than 30°, the weight increases to 1.2-1.5 times, Es When the angle is >60°, the weight is reduced to 0.8-1.0. During implementation, the weighted grid residual mean is mapped to the observation value based on the satellite's real-time azimuth and elevation angles (for example, an azimuth of 120.3° matches a 120°-125° grid). When outputting correction values, the phase observations must be multiplied by the wavelength of the corresponding frequency band to convert them into distance units. Multi-frequency weighted fusion optimizes spatial model accuracy, layered interpolation and smoothing ensure data integrity, and dynamic weighting enhances correction effectiveness in low-elevation areas.

[0053] In another technical solution, the specific steps of extracting trend complexity features include: For the residual data in each elevation and azimuth grid, linear, quadratic and cubic polynomials are used to perform trend fitting in turn, and the optimal order AIC value is selected by the Akaike Information Criterion. , where k is the number of model parameters, L is the model likelihood value, and the polynomial order with the smallest AIC value is selected as the final fitting model; Calculates R for a polynomial fit of a selected order 2 The trend significance is verified by F test, and the trend complexity feature judgment rule is: R 2 ≥0.3, the F test statistic is significant at the 95% confidence level; For multi-band signals, the residual trend is weighted and fused according to the wavelength difference; if the change of the trend fitting coefficient of the same satellite in more than 5 consecutive grids is less than 10%, the residual trend of the satellite is globally smoothed; when R is satisfied, 2 When ≥0.3 and the F test is significant, it is determined to be dominated by trend complexity; if multiple order models meet the conditions at the same time, the highest order model is selected.

[0054] For the residual data within each elevation-azimuth grid, a linear (first-order), quadratic, and cubic polynomial trend fit is performed, in sequence. The order is chosen based on the Akaike Information Criterion (AIC), whose calculation depends on the number of model parameters k and the likelihood value L. Lower AIC values ​​indicate better model fit. During implementation, the AIC values ​​for the third-order model should be calculated separately. For example, the AIC value for a cubic polynomial with parameter k=10 may be lower than that for a linear model with parameter k=3. Multipath errors require higher-order polynomials, and the AIC criterion balances fitting accuracy with model complexity. Optionally, when the number of residuals is insufficient (e.g., <20 data points per grid), a lower-order model is forced to avoid overfitting.

[0055] Goodness of fit and significance test are the core of determining trend dominance. Calculate the coefficient of determination R of the selected order polynomial fit 2 (value range 0-1), when R 2The fit is considered valid when the value is ≥0.3 (reference threshold 0.25-0.35). Simultaneously, an F test is performed to verify the significance of the trend: at a 95% confidence level, if the F statistic exceeds the critical value, it is considered significant. The implementation process includes: 1) performing polynomial regression on the grid residuals; 2) calculating R 2 value; 3) Check the F distribution table to verify the significance; 4) When R 2 ≥0.3 and the F test is significant, it is marked as trend complexity dominated. For example, a grid three-dimensional fitting R 2 =0.32 and the F value is significant, the trend dominant flag is triggered.

[0056] For signals in different frequency bands (such as GPS L1 / L2 / L5), the residual trend is weighted and fused according to the wavelength ratio: the frequency band with longer wavelength has a higher weight (such as the BDS-3 B3I weight is set to 0.4, and the GPS L1 weight is 0.2). If the trend fitting coefficient of the same satellite in more than 5 consecutive grids changes by less than 10% (for example, the difference in the azimuth coefficient c1 between adjacent grids is ≤0.05), a global smoothing process is performed on the satellite residual, and the trend coefficients of adjacent grids are fused using the moving average method. During implementation, the coefficient change rate is monitored and R is recalculated after smoothing is triggered. 2 Polynomial fitting and the AIC criterion significantly improve the ability to model nonlinear errors, the multi-band weighting strategy effectively integrates system differences, and the global smoothing mechanism enhances large-scale spatial consistency.

[0057] The specific operation mode of the trend surface analysis multi-path hemispherical graph model includes: For the residual data in each grid, the trend surface is fitted, and the model expressions of linear, quadratic and cubic trend surfaces are: in, is the multipath residual, A s is the azimuth, E s is the altitude angle, c0 to c9 are the polynomial coefficients of the trend surface; In view of the hardware delay deviation of multi-band signals, the trend surface coefficients are compensated for frequency band differences; the trend surface coefficients are dynamically weighted according to the real-time motion status of the satellite, and for grids without valid data, bilinear interpolation of the trend surface coefficients of adjacent grids is used to fill the gaps; the optimized trend surface coefficients are mapped to the observation values ​​according to the real-time satellite position, and the corrected pseudorange and phase data are output.

[0058] The model expression contains three levels: 1) Linear trend surface only contains azimuth A s and altitude angle E s 2) The quadratic trend surface increases 、 and Cross terms; 3) Further introduction of the third trend surface 、 Coefficients c0 to c9 characterize the trend surface. For example, if c1 > 0, multipath error increases with increasing azimuth. The corresponding expression is invoked based on the selected order. For example, a quadratic model requires the calculation of six coefficients, c0-c5. Polynomial surfaces are used to fit spatially varying error trends, and higher-order terms capture the characteristics of complex reflective environments.

[0059] Frequency band difference compensation and dynamic weighting optimize coefficient accuracy. To address multi-band hardware delay variations, frequency band compensation is applied to trend surface coefficients: for example, a +0.1m offset is added to the c0 term for the BDS-3 B2b band. Real-time satellite motion (e.g., velocity > 0.1° / s) triggers dynamic weighting: high-order terms are weighted down (e.g., the cubic term is weighted ×0.7) at high speeds, and restored at low speeds. The implementation process includes: 1) pre-storing compensation tables by frequency band; 2) acquiring satellite angular velocity in real time; and 3) adjusting coefficient weights based on motion. For example, a simplified linear model is used for high-speed GPS satellites, while a full cubic model is used for stationary receivers.

[0060] Data gap processing and real-time mapping achieve efficient correction. For grids without valid data, bilinear interpolation of the trend surface coefficients of the adjacent four grids (adjacent in the same layer azimuth + co-located in the upper and lower layers) is used to fill the gaps. For example, the coefficient of the missing grid = 0.25 × (upper left + upper right + lower left + lower right coefficients). In the correction stage, the optimized coefficients are substituted into the trend surface equation, and the real-time azimuth and altitude angles of the satellite (such as A s =45.3°, E s =28.7°), output multipath error value Δ TMHM During implementation, a coefficient-position mapping table is established, and corrections are updated on a second-by-second basis and superimposed on the original observations. A multi-order trend surface structure significantly enhances the modeling capabilities of complex reflection scenarios, motion state weighting effectively suppresses dynamic errors, and bilinear interpolation ensures correction continuity across the entire airspace.

[0061] In another technical solution, when the residual data simultaneously meets multiple judgment conditions among time repeatability dominance, spatial distribution dominance, or trend complexity dominance, the multipath error correction model is selected according to the following priority rules: The model selection priority is defined as: temporal repeatability dominant > trend complexity dominant > spatial distribution dominant; If the residual data satisfies both the temporal repeatability-dominant and other feature-dominant conditions, the sidereal day filter model is preferred; if the residual data satisfies both the trend complexity-dominant and spatial distribution-dominant conditions, the trend surface analysis multipath hemispherical graph model is preferred. When multiple feature-dominant conditions coexist, calculate the residual variance contribution corresponding to each feature, and the contribution ,in, is the residual variance corresponding to feature f, is the total residual variance; The model corresponding to the feature with the highest contribution is selected. If the difference in contribution is less than 10%, the combined model correction mechanism is enabled: Perform joint correction on multiple feature models with similar contribution, and the comprehensive correction amount Δ comb The calculation formula is: , where Δ SF is the correction amount of the SF model, Δ TMHM is the correction amount of the TMHM model, Δ MHM is the correction amount of the MHM model, ω1, ω2, ω3 are the contribution C of each model f The weight of the proportional distribution, and satisfy ; During the correction process, the precision reduction factor (PDOP) and residual variance of the positioning results are monitored in real time. If the PDOP value exceeds 3 or the residual variance decreases by less than 20% compared with the pre-correction value, the system dynamically switches to the model with the second highest contribution. For conflict scenarios that trigger model switching multiple times, the timestamp, satellite system, feature contribution, and positioning error data are recorded.

[0062] When residual data simultaneously meet the conditions of temporal repeatability, spatial distribution, or trend complexity, a three-level priority is established: temporal repeatability > trend complexity > spatial distribution. For example, if both temporal repeatability (autocorrelation coefficient > 0.8) and spatial distribution (high elevation angle density < 50% of low elevation angle density) are detected, the sidereal day filter (SF) model is preferred. The principle is that temporal repeatability, derived from the periodicity of satellite orbits, is more stable than variations in the spatial reflection environment, and prioritizing corrections can cover more significant error sources. A decision flag matrix is ​​established, triggering the first model that meets the conditions in order of priority.

[0063] The residual variance contribution calculation is used to resolve the priority conflicts at the same level. The contribution is defined as the ratio of the residual variance corresponding to a single feature to the total variance. For example, the time feature variance is 0.15, total variance When is 0.30, we get If the difference between the contribution of two features is less than 10% (e.g., time contribution is 42% and trend contribution is 38%), the combined model is used: the correction value Δ of SF, TMHM, and MHM models SF , Δ TMHM , Δ MHMWeighted fusion is performed based on the contribution ratio (for example, the SF weight is 0.42, the TMHM weight is 0.38, and the MHM weight is 0.20). The operation process includes decomposing the residual components of each feature, calculating the variance ratio, and weighting them to generate a comprehensive weight.

[0064] Dynamic switching and conflict monitoring ensure real-time reliability. During the correction process, the positioning precision dilution factor (PDOP) and residual variance are continuously monitored. If PDOP is greater than 3 (threshold range 2.5-4.0) or the residual variance decreases by less than 20% compared to before correction (for example, only by 15%), the system automatically switches to the model with the second highest contribution. For scenarios where switching is frequently triggered (e.g., ≥3 switches within 10 minutes), the timestamp, satellite system number, feature contribution, and positioning error data are recorded to form a conflict log for offline analysis. A real-time monitoring thread is established to interrupt the current model when the trigger condition is met and load the backup model coefficient library. Priority rules and contribution weighting significantly improve the rationality of decision-making in complex scenarios. The dynamic switching mechanism enhances the system's fault tolerance, and the conflict log supports long-term model optimization.

[0065] It should be noted that although the steps are described above in a specific order, this does not necessarily mean that the steps must be performed in this specific order. In fact, some of these steps can be performed concurrently or even in a different order, as long as the required functions can be achieved. The number of devices and processing scales described here are intended to simplify the description of the present invention. Applications, modifications, and variations of the present invention will be apparent to those skilled in the art.

[0066] Although the embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the description and implementation methods. They can be fully applied to various fields suitable for the present invention. For those familiar with the art, additional modifications can be easily implemented. Therefore, without departing from the general concept defined by the claims and the scope of equivalents, the present invention is not limited to the specific details and illustrations shown and described herein.

Claims

1. A method for reducing multipath errors of multi-frequency and multi-mode signals, characterized in that: The following steps are involved: Construct a multi-frequency, multi-mode, non-combined precise point positioning observation equation. The observation equation includes a combination of multi-frequency signals from multiple global navigation satellite systems. The observation equation is used to separate the multipath errors in pseudorange and phase observations. Quality control is performed on the separated pseudorange and phase residuals to eliminate abnormal residual data. Feature extraction and judgment of residual data after quality control, including: The temporal repeatability characteristics of the residual data were extracted by calculating the autocorrelation coefficient of the residual sequence. When the autocorrelation coefficient was greater than 0.7, it was determined to be dominated by temporal repeatability. The spatial distribution characteristics of the residual data were extracted by statistically analyzing the residual density in the elevation and azimuth grids. When the residual density in the grid with an elevation angle greater than 60° was less than 50% of the residual density in the grid with an elevation angle less than 30°, it was determined to be dominated by spatial distribution. The residual trend was fitted by a cubic polynomial and the goodness of fit R was calculated. 2 The value and F test value extract the trend complexity characteristics. When R 2 When the value is greater than or equal to 0.3 and the F test is significant at the 95% confidence level, it is determined to be dominated by trend complexity; The multipath error correction model is selected based on the feature judgment results. If it is determined that time repeatability is dominant, the sidereal day filter SF model is selected; if it is determined that spatial distribution is dominant, the multipath hemispherical map MHM model is selected; if it is determined that trend complexity is dominant, the trend surface analysis multipath hemispherical map TMHM model is selected; The multipath error correction of pseudorange observations and phase observations is generated through the selected correction model, the original observations are corrected and the positioning results are output.

2. The multi-path error reduction method for multi-frequency and multi-mode signals according to claim 1, characterized in that: The specific steps for constructing the multi-frequency, multi-mode, non-combined precise point positioning observation equations include: The multi-frequency signal combination of multiple global navigation satellite systems includes GPS three-frequency, GLONASS two-frequency, Galileo five-frequency and BDS-3 five-frequency signals. The basic observation equations of pseudorange and phase are: in, and are the original pseudorange observation value and phase observation value of the receiver r for the i-th frequency point of satellite s, is the geometric distance, and are the receiver and satellite clock errors, is the mapping function, T r is the zenith tropospheric delay; , is the i-th frequency value of satellite s, is the first frequency value, is the slant ionospheric delay of the first frequency; 、 、 、 are the hardware delay deviations between the receiver and the satellite for pseudorange and phase, Fuzziness The wavelength, and are pseudorange and phase multipath errors, respectively; and All are observation noise; When separating multipath errors, the inter-band deviation compensation is performed on the pseudorange and phase observations at different frequencies based on the frequency band differences of the satellite and receiver hardware delay deviations; and the multi-frequency observations are weightedly fused according to the signal-to-noise ratio of each frequency signal and the satellite elevation angle.

3. The multi-path error reduction method for multi-frequency and multi-mode signals according to claim 1, wherein: Specific steps in performing quality control include: Based on the normal distribution assumption, the mean μ is calculated for the pseudorange residual and phase residual in each elevation and azimuth grid. I and standard deviation σ I , the 3σ rule is used to identify outliers, and the judgment rules are: ,in, represents the i-th residual value in grid I; For the identified outliers, the F test is used to verify their significant impact on the grid residual variance. The specific steps are as follows: the residual data set in the grid is divided into the normal data set x I,1 and the abnormal data set x I,2 , and calculate the sample variance of the two data sets and , construct the F test statistic: If F is greater than the parameter α, take F of 0.05 α (n I,1 −1,n I,2 −1) value, it is determined that the outliers have a significant impact on the residual variance; For the residual data in the low elevation angle area with an elevation angle less than 15°, the abnormality judgment threshold is adjusted to: , where the satellite motion state dynamic factor , E s is the satellite altitude angle; For the grids that have eliminated outliers, the sliding window interpolation method of the residual mean of adjacent grids is used to fill the data gaps. The window size is 3×3 grids, and the interpolation formula is: , where N adj is the number of adjacent grids, μ I,j is the mean of residuals of adjacent grids; For the corrected residual data, the residual distribution is statistically analyzed by satellite system partition. If the residual variance of one system exceeds twice the global variance, the secondary quality control of the data of that system is triggered.

4. The multi-path error reduction method for multi-frequency and multi-mode signals according to claim 1, wherein: The specific steps of extracting time-repeated features include: The autocorrelation coefficients of pseudorange residual and phase residual sequences are calculated using sliding windows. The window length is set to an integer multiple of the satellite orbit repetition period, and the sliding step is 1 / 10 of the data sampling interval. A hypothesis test is performed on the autocorrelation coefficient of each window. If the autocorrelation coefficient is greater than 0.7 within the 95% confidence interval, the residual in the window is judged to have the dominant characteristic of time repeatability; According to the wavelength difference of signals in different frequency bands, the autocorrelation coefficients are subjected to frequency band weighted fusion, and the weight w i The calculation formula is: Among them, λ i is the wavelength of the signal in the i-th frequency band, n sys is the total number of frequency bands of the current satellite system; If the same satellite is judged to be time-repeatability dominated in multiple consecutive windows and the range of its orbital elevation angle variation is less than 10°, the residual sequence of the satellite is marked as time-repeatability dominated as a whole.

5. The multi-path error reduction method for multi-frequency and multi-mode signals according to claim 4, characterized in that: The specific operation mode of the sidereal day filter model includes: Calculate the orbit repetition period T based on the satellite ephemeris data orbs , and the residual sequence is T orbs Segment alignment, controlling the alignment error to no more than 5% of the data sampling interval; The mean of the pseudorange and phase residuals within each orbital period is calculated, and the abnormal segments that deviate from the mean by more than 2 standard deviations are eliminated. The inter-band deviation of the mean sequence is corrected for the hardware delay bias of the multi-band signal. The mean sequence is smoothed using a sliding window with a length of 3 orbital periods. The weights within the window are dynamically allocated according to the satellite elevation angle. The weight allocation formula is: ,in, is the average altitude angle of the satellite in the kth orbital period; The smoothed mean sequence is used as the correction value of the sidereal day filter model, matched to the original observation value by timestamp, and the corrected pseudorange and phase observation values ​​are output.

6. The multi-path error reduction method for multi-frequency and multi-mode signals according to claim 1, wherein: The specific steps of extracting spatial distribution features include: Dynamically divide the grid and dynamically adjust the grid's azimuth resolution according to the satellite elevation angle range. The specific adjustment rules are as follows: , where ΔA k is the adjusted azimuth resolution, ΔA0 is the reference azimuth resolution, and E k is the center value of the elevation angle interval of the kth layer, and the elevation angle range is divided into 5° intervals; the pseudorange residuals and phase residuals in each grid are statistically analyzed for density ρ I , the density calculation formula is: , where N I is the number of effective residuals in grid I, A I is the spherical area of ​​the grid; For high elevation angle areas with elevation angles > 60°, the residual density is corrected to , density compensation factor , E s is the satellite elevation angle; if the corrected average residual density of the high-elevation area grid is less than 50% of the average residual density of the low-elevation area grid with an elevation angle of less than 30°, it is determined to be spatially dominant; The residual density of all frequency bands of the same satellite system is jointly analyzed. If more than 70% of the frequency bands meet the spatial distribution dominance condition, the system is globally marked as spatial distribution dominance.

7. The multi-path error reduction method for multi-frequency and multi-mode signals according to claim 6, characterized in that: The specific operation of the multipath hemispheric graph model includes: For each dynamically divided grid, the mean of pseudorange residuals and phase residuals is calculated respectively. For multi-band signals, the residual means are weighted and fused according to the wavelength of the frequency band. For grids without valid data, the residual means of the adjacent three layers of elevation grids are used for interpolation, and the interpolated grid data are spatially smoothed using a Gaussian kernel function with a standard deviation of 2°. When generating the correction value, the grid residual mean is dynamically weighted according to the real-time satellite elevation angle to enhance the correction contribution of the low elevation angle area; the weighted grid residual mean is mapped to the original observation value according to the real-time satellite position, and the corrected pseudorange and phase observation values ​​are output.

8. The multi-path error reduction method for multi-frequency and multi-mode signals according to claim 1, characterized in that: The specific steps of extracting trend complexity features include: For the residual data in each elevation and azimuth grid, linear, quadratic and cubic polynomials are used to perform trend fitting in turn, and the optimal order AIC value is selected by the Akaike Information Criterion. , where k is the number of model parameters, L is the model likelihood value, and the polynomial order with the smallest AIC value is selected as the final fitting model; Calculates R for a polynomial fit of a selected order 2 The trend significance is verified by F test, and the trend complexity feature judgment rule is: R 2 ≥0.3, the F test statistic is significant at the 95% confidence level; For multi-band signals, the residual trend is weighted and fused according to the wavelength difference; if the change of the trend fitting coefficient of the same satellite in more than 5 consecutive grids is less than 10%, the residual trend of the satellite is globally smoothed; when R is satisfied, 2 When ≥0.3 and the F test is significant, it is determined to be dominated by trend complexity; if multiple order models meet the conditions at the same time, the highest order model is selected.

9. The method for reducing multipath errors of multi-frequency and multi-mode signals according to claim 8, wherein: The specific operation mode of the trend surface analysis multi-path hemispherical graph model includes: For the residual data in each grid, the trend surface is fitted, and the model expressions of linear, quadratic and cubic trend surfaces are: in, is the multipath residual, A s is the azimuth, E s is the altitude angle, c0 to c9 are the polynomial coefficients of the trend surface; In view of the hardware delay deviation of multi-band signals, the trend surface coefficients are compensated for frequency band differences; the trend surface coefficients are dynamically weighted according to the real-time motion status of the satellite, and for grids without valid data, bilinear interpolation of the trend surface coefficients of adjacent grids is used to fill the gaps; the optimized trend surface coefficients are mapped to the observation values ​​according to the real-time satellite position, and the corrected pseudorange and phase data are output.

10. The multi-path error reduction method for multi-frequency and multi-mode signals according to claim 1, characterized in that: When the residual data simultaneously meets multiple criteria of temporal repeatability dominance, spatial distribution dominance, or trend complexity dominance, the multipath error correction model is selected according to the following priority rules: The model selection priority is defined as: temporal repeatability dominant > trend complexity dominant > spatial distribution dominant; If the residual data satisfies both the temporal repeatability-dominant and other feature-dominant conditions, the sidereal day filter model is preferred; if the residual data satisfies both the trend complexity-dominant and spatial distribution-dominant conditions, the trend surface analysis multipath hemispherical graph model is preferred. When multiple feature-dominant conditions coexist, calculate the residual variance contribution corresponding to each feature, and the contribution ,in, is the residual variance corresponding to feature f, is the total residual variance; The model corresponding to the feature with the highest contribution is selected. If the difference in contribution is less than 10%, the combined model correction mechanism is enabled: Perform joint correction on multiple feature models with similar contribution, and the comprehensive correction amount Δ comb The calculation formula is: , where Δ SF is the correction amount of the SF model, Δ TMHM is the correction amount of the TMHM model, Δ MHM is the correction amount of the MHM model, ω1, ω2, ω3 are the contribution C of each model f Proportional weight distribution, and satisfy ; During the correction process, the precision reduction factor (PDOP) and residual variance of the positioning results are monitored in real time. If the PDOP value exceeds 3 or the residual variance decreases by less than 20% compared with the pre-correction value, the system dynamically switches to the model with the second highest contribution. For conflict scenarios that trigger model switching multiple times, the timestamp, satellite system, feature contribution, and positioning error data are recorded.

Images