Method for mitigating multipath errors in multi-frequency multi-mode signals

Abstract

The application discloses a multi-frequency multi-mode signal multipath error weakening method, and relates to the technical field of precise positioning. The method comprises the following steps: constructing a multi-frequency multi-mode non-combination precise point positioning observation equation, separating multipath errors in pseudo-range and phase observation values through the equation, and forming residual data; performing quality control on the residual data, and eliminating abnormal values; performing feature extraction and judgment based on the residual data, including extracting time repeatability features, spatial distribution features and trend complexity features; selecting a correction model according to the feature judgment result, including a constant star day filtering model, a multipath hemisphere diagram model and a trend surface analysis multipath hemisphere diagram model; and generating a multipath error correction value through the selected model, correcting original observation values and outputting a positioning result. The method solves the problem of positioning precision reduction caused by multipath errors, can effectively weaken multipath errors, improve positioning precision and stability, is suitable for complex environments and reduces data processing redundancy.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of precise positioning, and in particular to a multi-frequency multi-mode signal multi-path error weakening method. BACKGROUND

[0002] In GNSS precise positioning technology, multi-path error is one of the main factors affecting positioning accuracy. This error is caused by the superimposed interference of reflected waves and direct waves when satellite signals encounter buildings, ground or other reflecting bodies during transmission, resulting in deviations in receiver observations (such as pseudorange and phase). Especially in precise point positioning (PPP) applications, such as surveying, geological disaster monitoring and autonomous driving fields, the positioning result has strict requirements for millimeter to centimeter level accuracy, and the cumulative effect of multi-path error can significantly reduce the positioning reliability. For example, in urban canyons, forests or environments with significant multi-path effects, pseudorange observations may produce decimeter-level deviations, phase observations may introduce cycle slips or ambiguity resolution errors, and thus affect the convergence speed and stability of the final position solution.

[0003] In the prior art, multi-path error weakening methods mainly include stellar day filtering (SF) and multi-path hemisphere map (MHM) models. The SF model uses the periodic characteristics of satellite orbits to correct errors through time repetition, but the effect is limited in dynamic environments (such as moving receivers or changes in reflecting bodies), because the errors may exhibit non-repetitive characteristics. The MHM model is based on spatial grid statistical residual distribution and is suitable for static scenarios, but in low elevation areas or when signal blocking is frequent, the sparse grid data results in inaccurate correction. In addition, existing methods are mostly for a single frequency band or system (such as only GPS L1 frequency band), while modern GNSS has developed to multi-frequency multi-mode fusion (such as GPS, GLONASS, Galileo, BDS-3 in parallel), and the differences in signal frequency bands (such as hardware delay bias) are not fully compensated, resulting in incomplete error separation. The existing technology also lacks an adaptive mechanism: the model selection is fixed and cannot dynamically switch according to error characteristics (such as time repetition, spatial distribution or trend complexity), and in complex multi-path scenarios, the correction effect decreases and the residual variance may only be reduced by 10-15%, the positioning dilution of precision (PDOP) is prone to exceed the threshold, affecting the practicality. Therefore, there is an urgent need for a method that can adapt to multi-frequency multi-mode signals, dynamically identify error dominant characteristics and optimize model selection, to improve the robustness and accuracy of multi-path error weakening. SUMMARY

[0004] The present application overcomes the problem of reduced positioning accuracy caused by multi-path error in the prior art, and provides a multi-frequency multi-mode signal multi-path error weakening method, which can effectively weaken multi-path error, improve positioning accuracy and stability, adapt to complex environments and reduce data processing redundancy.

[0005] To achieve the above object, the present application adopts the following scheme:

[0006] The method for weakening multipath error of multi-frequency multi-mode signal comprises the following steps:

[0007] A multi-frequency multi-mode non-combination precise point positioning observation equation is constructed, the observation equation comprises a multi-frequency signal combination from multiple global navigation satellite systems, the multipath error in the pseudo-range observation value and the phase observation value is separated through the observation equation, and quality control is performed on the separated pseudo-range residual and phase residual to eliminate abnormal residual data;

[0008] Feature extraction and determination are performed on the residual data after quality control, including:

[0009] The time repetition feature of the residual data is extracted by calculating the autocorrelation coefficient of the residual sequence, when the autocorrelation coefficient is greater than 0.7, it is determined that the time repetition is dominant; the spatial distribution feature of the residual data is extracted by counting the residual density in the elevation and azimuth grid, when the residual density in the grid with an elevation greater than 60° is lower than 50% of the residual density in the grid with an elevation less than 30°, it is determined that the spatial distribution is dominant; the trend complexity feature is extracted by fitting the residual trend with a cubic polynomial and calculating the goodness of fit R 2 value and the F test value, when the R 2 value is greater than or equal to 0.3 and the F test is significant at a 95% confidence level, it is determined that the trend complexity is dominant;

[0010] According to the feature determination result, a multipath error correction model is selected, if it is determined that the time repetition is dominant, a sidereal day filter SF model is selected; if it is determined that the spatial distribution is dominant, a multipath hemisphere map MHM model is selected; if it is determined that the trend complexity is dominant, a trend surface analysis multipath hemisphere map TMHM model is selected;

[0011] The multipath error correction amount of the pseudo-range observation value and the phase observation value is generated through the selected correction model, the original observation value is corrected, and a positioning result is output.

[0012] As a preferred, the specific steps of constructing the multi-frequency multi-mode non-combination precise point positioning observation equation include:

[0013] The multi-frequency signal combination of the multiple global navigation satellite systems includes GPS three-frequency, GLONASS two-frequency, Galileo five-frequency and BDS-3 five-frequency signals, and the pseudo-range and phase basic observation equations are:

[0014]

[0015] wherein, and are the original pseudo-range observation value and the phase observation value of the receiver r to the i-th frequency point of the satellite s, respectively, is the geometric distance, and are the receiver and satellite clock errors, respectively, is the mapping function, T r is the zenith tropospheric delay; , is the satellite s the i-th frequency value, is the first frequency value, is the slant ionospheric delay for the first frequency; , , , are the receiver and satellite end hardware delay biases for pseudo-range and phase, respectively, is the ambiguity is the wavelength, and are the pseudo-range and phase multipath errors, respectively; and are the observation noises;

[0016] In separating the multipath errors, the inter-frequency biases of the pseudo-range and phase observations at different frequencies are compensated based on the frequency band differences of the satellite and receiver hardware delay biases. The multi-frequency observations are fused by weighting according to the signal-to-noise ratios of the signals at different frequencies and the satellite elevation angles.

[0017] As preferred, the specific steps of performing the quality control include:

[0018] Based on the normal distribution assumption, the mean μ I and the standard deviation σ I of the pseudo-range residuals and the phase residuals in each grid of the elevation and azimuth angles are calculated, respectively, and the abnormal values are identified by the 3σ rule, and the determination rule is: wherein, represents the i-th residual value in the grid I;

[0019] For the identified abnormal values, the significant influence of the abnormal values on the grid residual variance is verified by F-test, and the specific steps are as follows: the residual data set in the grid is divided into a normal data set x I,1 and an abnormal data set x I,2 , and the sample variances of the two data sets are calculated and , and the F-test statistic is constructed: If the F is greater than the F α (n I,1 −1,n I,2 −1) value taken as 0.05 of the parameter α, it is determined that the abnormal value has a significant influence on the residual variance;

[0020] For the residual data in the low-elevation area with an elevation angle < 15°, the abnormal determination threshold is adjusted to: wherein the satellite motion state dynamic factor , E s is the satellite elevation angle;

[0021] For the grid that eliminates outliers, a sliding window interpolation method is used to fill in the data gaps, with a window size of 3x3 grids, and the interpolation formula is: wherein N adj is the number of adjacent grids, and μ I,j is the mean residual of adjacent grids;

[0022] For the corrected residual data, the residual distribution is counted according to the satellite system partition, and if the residual variance of one of the systems exceeds twice the global variance, the secondary quality control of the data of the system is triggered.

[0023] As preferred, the specific steps of extracting the time repetition feature include:

[0024] The sliding window is used to calculate the autocorrelation coefficient for the pseudo-range residual and phase residual sequences respectively, and 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;

[0025] The autocorrelation coefficient of each window is subjected to hypothesis testing, and if the autocorrelation coefficient is greater than 0.7 within the 95% confidence interval, it is determined that the residual in the window has a time repetition dominant feature;

[0026] According to the wavelength difference of signals of different frequency bands, the autocorrelation coefficients are weighted and fused, and the weight w i is calculated as: wherein λ i is the wavelength of the i-th frequency band signal, and n sys is the total number of frequency bands of the current satellite system;

[0027] If the same satellite is determined to be time repetition dominant in multiple consecutive windows, and the range of its orbit elevation angle change is less than 10°, the residual sequence of the satellite is marked as time repetition dominant as a whole.

[0028] As preferred, the specific operation mode of the star day filter model includes:

[0029] According to the satellite ephemeris data, the orbit repetition period T orbs is calculated, and the residual sequence is segmented and aligned according to T orbs , and the alignment error is controlled to be not more than 5% of the data sampling interval;

[0030] The mean value of the pseudo-range and phase residual in each orbit period is calculated respectively, and the abnormal segment deviating from the mean value more than 2 times of the standard deviation is eliminated; for the hardware delay deviation of the multi-band signal, the mean value sequence is corrected among bands; the mean value sequence is smoothed by using a sliding window with a length of 3 orbit periods, and the weight in the window is dynamically allocated according to the satellite elevation angle, and the weight allocation formula is: wherein, is the average elevation angle of the satellite in the kth orbit period.

[0031] The smoothed mean value sequence is taken as the correction quantity of the sidereal day filtering model, is matched to the original observation value according to the time stamp, and the corrected pseudo-range and phase observation value are output.

[0032] As preferred, the specific steps of extracting the spatial distribution feature include:

[0033] The grid is dynamically divided, and the azimuth resolution of the grid is dynamically adjusted according to the satellite elevation angle range, and the specific adjustment rule is: wherein, ΔA k is the adjusted azimuth resolution, ΔA0 is the reference azimuth resolution, E k is the center value of the kth elevation angle interval, and the elevation angle range is divided into 5° intervals; the density ρ I of the pseudo-range residual and the phase residual in each grid is counted respectively, and the density calculation formula is: wherein, N I is the number of effective residuals in the grid I, and A I is the spherical area of the grid;

[0034] For the high-elevation-angle region with an elevation angle > 60°, the residual density is corrected as , the density compensation factor , E s is the satellite elevation angle; if the average residual density of the grid in the high-elevation-angle region after correction is lower than 50% of the average residual density of the grid in the low-elevation-angle region with an elevation angle < 30°, it is determined that the spatial distribution is dominant.

[0035] The residual densities of all bands of the same satellite system are jointly analyzed, and if more than 70% of the bands meet the spatial distribution dominant condition, the system is globally marked as spatial distribution dominant.

[0036] As preferred, the specific operation mode of the multipath hemisphere diagram model includes:

[0037] The mean values of the pseudo-range residuals and the phase residuals are calculated for each dynamically divided grid, respectively; for multi-band signals, the residual mean values are weighted and fused according to the wavelengths of the bands; for grids without effective data, the residual mean values of adjacent 3 layers of height angle grids are used for interpolation, and the interpolated grid data is spatially smoothed by using a Gaussian kernel function with a standard deviation of 2°;

[0038] In generating the correction amount, the grid residual mean values are dynamically weighted according to the real-time elevation angles of the satellites to enhance the correction contribution of the low-elevation-angle area; the weighted grid residual mean values are mapped to the original observation values according to the real-time positions of the satellites, and the corrected pseudo-range and phase observation values are output.

[0039] As preferred, the specific steps of extracting the trend complexity feature include:

[0040] For the residual data in each elevation angle and azimuth angle grid, linear, quadratic and cubic polynomials are used in turn for trend fitting, and the Akaike Information Criterion (AIC) value is selected to select the optimal order, wherein k is the number of model parameters, L is the model likelihood value, and the polynomial order with the minimum AIC value is selected as the final fitting model;

[0041] The R 2 value of the polynomial fitting of the selected order is calculated, and the trend significance is verified by F test, and the trend complexity feature determination rule is: R 2 ≥ 0.3, and the F test statistic is significant at the 95% confidence level;

[0042] For multi-band signals, the residual trend is weighted and fused according to the wavelength difference; if the trend fitting coefficient of the same satellite changes by less than 10% in more than 5 consecutive grids, the residual trend of the satellite is globally smoothed; when R 2 ≥ 0.3 and F test is significant, it is determined that the trend complexity is dominant; if multiple order models meet the condition at the same time, the highest order model is selected.

[0043] As preferred, the specific operation mode of the trend surface analysis multi-path hemisphere map model includes:

[0044] For the residual data in each grid, the trend surface is fitted, and the model expressions of linear, quadratic and cubic trend surfaces are:

[0045]

[0046] wherein, is the multi-path residual, A s is the azimuth angle, E s is the elevation angle, and c0 to c9 are the polynomial coefficients of the trend surface;

[0047] The trend surface coefficients are compensated for the hardware delay deviation of the multi-band signal, the trend surface coefficients are dynamically weighted according to the real-time motion state of the satellite, and the bilinear interpolation of the trend surface coefficients of adjacent grids is used to fill the gap for the grids without effective data.

[0048] As preferred, when the residual data simultaneously satisfy multiple determination conditions in the time repetition dominance, the spatial distribution dominance or the trend complexity dominance, the multipath error correction model is selected by the following priority rules:

[0049] The model selection priority is defined as: time repetition dominance > trend complexity dominance > spatial distribution dominance;

[0050] If the residual data simultaneously satisfy the time repetition dominance and other feature dominance conditions, the star day filtering model is preferentially selected; if the trend complexity dominance and the spatial distribution dominance are simultaneously satisfied, the trend surface analysis multipath hemispherical map model is preferentially selected;

[0051] When multiple feature dominance conditions coexist, the residual variance contribution degrees corresponding to the features are calculated, and the contribution degree of the feature with the highest contribution degree is selected. wherein, is the residual variance corresponding to the feature f, is the total residual variance;

[0052] If the difference between the contribution degrees is less than 10%, a combined model correction mechanism is enabled:

[0053] Multiple feature models with similar contribution degrees are jointly corrected, and the comprehensive correction amount Δ comb The calculation formula is: wherein Δ SF is the correction amount of the SF model, Δ TMHM is the correction amount of the TMHM model, and Δ MHM is the correction amount of the MHM model, ω1, ω2, and ω3 are weights allocated in proportion to the contribution degrees C f of the models, and satisfy ;

[0054] In the correction process, the accuracy decay factor PDOP and the residual variance of the positioning result are monitored in real time, if the PDOP value exceeds 3 or the residual variance is insufficiently reduced by 20% compared with before correction, the model with the second highest contribution degree is dynamically switched to; for the conflict scene of multiple triggered model switching, the time stamp, satellite system, feature contribution degree and positioning error data are recorded.

[0055] The present application at least includes the following beneficial effects: (1) by dynamic feature extraction and adaptive model selection, effectively weakening multipath error, improving the precision and stability of precise point positioning, and significantly reducing positioning residual variance in complex environments; (2) by using multi-frequency signal combination, through frequency band deviation compensation and weighted fusion based on signal-to-noise ratio and elevation angle, improving error separation precision and enhancing multi-system compatibility; (3) by quality control steps, including outlier identification, F-test verification, dynamic threshold adjustment and data interpolation, ensuring the reliability of residual data, reducing abnormal interference and improving the accuracy of subsequent feature extraction; (4) the feature extraction mechanism covers the diversity of multipath error, combined with the correction model, to improve the robustness of the system in variable environments; (5) by priority rules and combination model mechanism, processing feature conflicts, optimizing correction effect; real-time monitoring of PDOP and residual variance and dynamically switching models to enhance fault tolerance and avoid positioning failure. BRIEF DESCRIPTION OF DRAWINGS

[0056] Figure 1 The principle flowchart of the multi-frequency multi-mode signal multipath error weakening method provided by the present application. DETAILED DESCRIPTION

[0057] The present application will be further described in detail below with reference to the accompanying drawings, so that those skilled in the art can implement it according to the description.

[0058] As shown in the drawings, the multi-frequency multi-mode signal multipath error weakening method provided by the present application includes the following steps: Figure 1

[0059] Constructing a multi-frequency multi-mode non-combination precise point positioning observation equation, the observation equation including multi-frequency signal combination from multiple global navigation satellite systems, separating multipath error in pseudo-range observation value and phase observation value through the observation equation, and performing quality control on the separated pseudo-range residual and phase residual to eliminate abnormal residual data;

[0060] Extracting and determining the features of the residual data after quality control, including:

[0061] Extracting the time repeatability feature of the residual data by calculating the autocorrelation coefficient of the residual sequence, and determining that time repeatability is dominant when the autocorrelation coefficient is greater than 0.7; extracting the spatial distribution feature of the residual data by counting the residual density in the elevation and azimuth grid, and determining that spatial distribution is dominant when the residual density in the grid with an elevation greater than 60° is less than 50% of the residual density in the grid with an elevation less than 30°; extracting the trend complexity feature by fitting the residual trend with a cubic polynomial and calculating the goodness of fit R 2 value and F-test value, and determining that trend complexity is dominant when the R 2 value is greater than or equal to 0.3 and the F-test is significant at a 95% confidence level.​

[0062] According to the feature judgment result, a multipath error correction model is selected, if it is judged that time repetition is dominant, a sidereal filter SF model is selected, if it is judged that spatial distribution is dominant, a multipath hemisphere map MHM model is selected, and if it is judged that trend complexity is dominant, a trend surface analysis multipath hemisphere map TMHM model is selected.

[0063] The multipath error correction amount of the pseudo-range observation value and the phase observation value is generated through the selected correction model, the original observation value is corrected, and the positioning result is output.

[0064] The multi-frequency multi-mode non-combination precise point positioning observation equation is constructed for separating the multipath error, including constructing the observation equation, and separating the multipath error in the pseudo-range and phase observation value. The multi-frequency signal combination (such as GPS three-frequency, GLONASS two-frequency, Galileo five-frequency, and BDS-3 five-frequency) of the global navigation satellite system (GNSS) is used, the signal propagation process is described by mathematical equations, the multipath error is separated from other system errors (such as clock error, ionospheric delay), and pseudo-range residuals and phase residuals are formed. These residuals contain multipath effects and random noise, which are input data for subsequent processing.

[0065] The observation equation is based on the non-combination precise point positioning (PPP) model, avoids using traditional differential technology, and thus retains the independent information of each frequency band. The parameters such as the geometric distance between the satellite and the receiver, the clock error, the zenith tropospheric delay, the slant ionospheric delay, and the hardware delay bias are considered in the equation. The residual sequence is calculated by substituting the actual observation value (pseudo-range and phase). First, the multi-system multi-frequency signal data is received; second, the equation is applied to compensate for the inter-frequency bias of the observation value of each frequency band (for example, the signal weight is adjusted according to the hardware delay difference); and finally, the residual data is separated. When the technical features are matched, the multi-frequency signal combination provides more abundant error separation dimensions, and the non-combination PPP model ensures that the residuals directly reflect the multipath effect, laying a data foundation for subsequent quality control.

[0066] Quality control is performed on the separated pseudo-range residuals and phase residuals to ensure data reliability, and the accuracy of subsequent feature extraction is improved by eliminating abnormal residuals. Quality control includes identifying, verifying, and eliminating abnormal data. Based on the assumption that the residual data is subject to normal distribution, statistical methods (such as the 3σ rule) are used to identify outliers, and the significance of abnormal values is verified by F test. For low-elevation areas (such as an elevation angle less than 15°), satellite signals are easily affected by dynamic interference, so the threshold needs to be dynamically adjusted (for example, a satellite motion state factor that changes according to the satellite elevation angle is introduced).

[0067] Feature extraction is the key of multipath model selection, which is divided into three sub-features: temporal repetitiveness, spatial distribution and trend complexity. Each sub-feature corresponds to a decision logic, which is used to identify the dominant characteristics of multipath error.

[0068] Temporal repetitiveness feature extraction is to extract temporal repetitiveness by computing the autocorrelation coefficient (ACF) of residual series. The periodic operation of satellite orbit leads to the repetition of multipath error in time (e.g. the orbit repetition period dependent on the star day filter). When the autocorrelation coefficient is greater than 0.7 (reference threshold), it is determined that temporal repetitiveness is dominant, indicating that the error is mainly in a periodic pattern.

[0069] Spatial distribution feature extraction is to extract spatial distribution characteristics by counting the residual density in the elevation-azimuth grid. Multipath error is affected by the position of the reflector, and presents uneven distribution in space (e.g. near-ground reflection leads to dense error density at low elevation). When the residual density of the grid with elevation greater than 60° is less than 50% of the grid density with elevation less than 30° (reference threshold), it is determined that spatial distribution is dominant, indicating that the error is mainly in azimuth dependence.

[0070] Trend complexity feature extraction is to extract trend complexity by fitting the residual trend with a cubic polynomial and computing the goodness of fit R 2 and F-test value. Multipath error may present a nonlinear trend (e.g. complex changes caused by reflector movement). When R 2 ≥ 0.3 (reference threshold) and F-test is significant at 95% confidence level, it is determined that trend complexity is dominant, indicating that the error is mainly in high-order change. Including fitting linear, quadratic and cubic polynomials in turn; select the optimal order (e.g. based on AIC criterion); compute R 2 and F value to verify significance. When the three features are combined, the polynomial fitting captures the nonlinear pattern, the statistical test ensures the reliability, and the multi-frequency weighting optimizes the trend extraction.

[0071] In the overall feature extraction section, the three sub-features run in parallel: temporal repetitiveness focuses on periodicity, spatial distribution pays attention to azimuth dependence, and trend complexity handles nonlinear change. They cooperate to cover the diversity of multipath error through independent decision, providing the basis for model selection.

[0072] According to the feature judgment result, the multipath error correction model is selected, and the optimal correction model is dynamically selected based on the feature judgment. The feature dominant result is mapped to the corresponding model: the time repeatability dominant selects the sidereal filter (SF) model, the spatial distribution dominant selects the multipath hemisphere map (MHM) model, and the trend complexity dominant selects the trend surface analysis multipath hemisphere map (TMHM) model. Each model is aimed at specific error characteristics (such as SF using time repeatability, MHM using spatial distribution), and the model selection ensures that the correction accurately matches the error mode. The specific process includes: receiving the feature judgment result, processing conflicts according to priority rules (for example, time dominant > trend dominant > spatial dominant); if multiple features coexist, calculate the residual variance contribution degree, and select the highest contributor. When each feature is matched, the priority rule solves the judgment overlap, the contribution degree quantifies the feature influence, and the multi-model combination mechanism enhances the robustness. The model selection dynamically optimizes the correction strategy, avoiding the limitations of a single model.

[0073] The correction amount is generated and corrected by the selected correction model, and the multipath error correction amount is generated by the selected model, and the original observation value is corrected to output the positioning result. It involves generating correction, correcting observation and outputting positioning. The model (SF / MHM / TMHM) is used to invert the multipath error value from the residual data, which is applied as a negative compensation to the original pseudorange and phase observation value, thereby weakening the error influence, and finally calculating the high-precision position through the PPP algorithm. The specific process is: running the selected model (such as SF model calculating time average residual, MHM model based on grid mean, TMHM model fitting trend surface coefficient); generate correction table (store model output value), match satellite real-time position (elevation angle, azimuth angle) to obtain correction; correct the original observation value (pseudorange and phase), solve the non-combined PPP equation to output the positioning result. When each feature is matched, the model generates real-time correction, and the position matching ensures dynamic adaptability, and the correction step is directly integrated into the positioning process to realize end-to-end error weakening.

[0074] This method realizes multi-path error mitigation through systematic steps, effectively suppresses the influence of multi-path error on pseudorange and phase observations by dynamically selecting the optimal correction model, thereby improving the stability and reliability of PPP positioning. The feature extraction decision mechanism (time, space, trend) covers the diversity of multi-path effects, making the method applicable to different observation scenarios (such as urban canyons, open areas), enhancing system robustness. The combination 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 multi-mode signal processing integrates the advantages of GPS, GLONASS, Galileo, BDS-3, etc. systems, through frequency band weighting and bias compensation, to comprehensively improve global navigation performance. Overall, this method achieves high practicality and universality in the field of multi-path error mitigation. For the obtained positioning data, the most suitable model can be automatically identified and selected for multi-path error correction, achieving high efficiency and accuracy, avoiding the influence of human judgment errors on positioning accuracy, and providing reliable support for high-precision positioning.

[0075] In another technical solution, the specific steps of constructing a multi-frequency multi-mode non-combined precise point positioning observation equation include:

[0076] 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, and the basic observation equations for pseudorange and phase are:

[0077]

[0078] wherein, and are the original pseudorange and phase observation values of the receiver r to the satellite s at the i-th frequency point, is the geometric distance, and are the receiver and satellite clock errors, respectively, is the mapping function, T r is the zenith troposphere delay; , is the i-th frequency value of the satellite s, is the first frequency value, is the first frequency slant ionospheric delay; , , , are the receiver and satellite end hardware delay biases of pseudorange and phase, respectively, is the wavelength of ambiguity , and are the pseudorange and phase multi-path errors, respectively; and All of these are observation noise;

[0079] When separating multipath errors, frequency band deviation compensation is performed on pseudorange and phase observations at different frequency points based on the frequency band difference of the hardware delay deviation between the satellite and the receiver; and multi-frequency observations are weighted and fused according to the signal-to-noise ratio of each frequency signal and the satellite elevation angle.

[0080] The multi-frequency signal combination covers GPS three-frequency, GLONASS two-frequency, Galileo five-frequency, and BDS-3 five-frequency signals, enhancing observation redundancy by fusing multi-system, multi-band data. Different frequency bands exhibit varying sensitivities to multipath errors; for example, high-frequency signals are more susceptible to reflectors, while low-frequency signals are more significantly affected by ionospheric delay. Multi-frequency combinations provide complementary error separation mechanisms. In equation construction, geometric distance parameters include corrections for antenna phase center offset and relativistic effects, requiring precise calculation using accurate ephemeris and station coordinates. Clock bias parameters are divided into receiver clock bias and satellite clock bias. Based on the IGS protocol, pseudorange hardware delay bias is absorbed, ensuring that clock errors and hardware biases are eliminated synchronously. Zenith tropospheric delay is converted to slant path delay using a mapping function, while slant ionospheric delay is correlated with each frequency point through coefficients and absorbs inter-band differential code bias.

[0081] Hardware delay bias is divided into receiver-side and satellite-side biases. For different frequency bands (such as the B1I and B2B bands of BDS-3), nanosecond-level differences may exist, which are compensated for using pre-calibrated parameters or global averages. Weighted fusion of observations relies on signal-to-noise ratio (SNR) and satellite elevation angle: when the SNR is below a certain threshold (e.g., 30-35 dBHz), the weights are reduced proportionally; when the elevation angle is below 15°, the weights are further reduced to 30%-60% of their original value to suppress weak signal interference. During runtime, the original observations are first input, and precise ephemeris is used to calculate the geometric distance. Then, clock errors and hardware bias parameters are substituted, and finally, pseudorange and phase residuals are output through frequency band compensation and dynamic weighting. This scheme improves error separation accuracy and enhances data reliability. Multi-frequency fusion and dynamic weighting significantly improve error separation accuracy and support multipath modeling in complex environments.

[0082] In another technical solution, the specific steps for performing quality control include:

[0083] Based on the normal distribution assumption, the mean μ is calculated for the pseudorange residual and phase residual within each elevation and azimuth grid. I and standard deviation σ I The 3σ rule is used to identify outliers, and the judgment rule is as follows: ,in, This represents the i-th residual value within grid I;

[0084] For the identified outliers, the significance of their influence on the grid residual variance is verified by F-test. The specific steps are as follows: the residual data set in the grid is divided into a normal data set x I,1 and an abnormal data set x I,2 , and the sample variances of the two data sets are calculated and , and the F-test statistic is constructed: If F is greater than the F α (n I,1 −1,n I,2 −1) value with an α of 0.05, it is determined that the outliers have a significant influence on the residual variance.

[0085] For the residual data in the low-elevation-angle area with an elevation angle < 15°, the abnormality judgment threshold is adjusted to: wherein the satellite motion state dynamic factor , and E s is the satellite elevation angle.

[0086] For the grid with the outliers removed, the adjacent grid residual mean sliding window interpolation method is used to fill in the data gaps, and the window size is 3x3 grids. The interpolation formula is: wherein N adj is the number of adjacent grids, and μ I,j is the adjacent grid residual mean.

[0087] For the corrected residual data, the residual distribution is statistically analyzed according to the satellite system partition. If the residual variance of one system exceeds twice the global variance, the secondary quality control of the data of the system is triggered.

[0088] Based on the normal distribution assumption, the mean and standard deviation of the residual in each elevation-angle-azimuth-angle grid need to be calculated. The 3σ rule is used for initial screening of outliers: if the deviation of the residual from the mean exceeds 3 times the standard deviation, it is marked as a suspected abnormality. For the low-elevation-angle area (elevation angle < 15°), due to the large signal fluctuation caused by the intense satellite motion, a dynamic adjustment factor is introduced to expand the threshold range. The factor is negatively correlated with the elevation angle. For example, when the elevation angle is 10°, the threshold is relaxed to 3.5 times the standard deviation, and when the elevation angle is 5°, it is further relaxed to 4 times the standard deviation, avoiding excessive removal of valid data.

[0089] The significance of the outliers needs to be verified by F-test. The grid residual is divided into a normal data set and an abnormal data set, and the ratio of the sample variances of the two groups is calculated as the F statistic. If the value exceeds the critical value of the F distribution corresponding to the selected significance level (such as α=0.05), it is determined that the outliers have a significant influence on the overall variance, and are removed. The removal operation is performed in an iterative manner: each time the residual with the largest deviation from the mean is removed, the statistic is recalculated until the F-test is not significant. This process ensures the scientificity of the outlier removal and avoids subjective misjudgment.

[0090] For grids without valid data, the residual mean of adjacent 3x3 grids is interpolated (for example, the average of the mean of the 8 surrounding grids). After interpolation, the residual distribution is counted by satellite system partition, and if the residual variance of a single system (such as GLONASS) exceeds 2 times the global variance, a secondary quality control for the system is triggered: the abnormal detection and rejection process is re-executed. This mechanism effectively intercepts local systematic abnormalities and ensures the overall consistency of the data. The quality control mechanism enhances the reliability of the data, and the dynamic threshold and secondary quality inspection design significantly reduce the false rejection rate, providing high-integrity residual data for subsequent feature extraction.

[0091] In another technical solution, the specific steps of extracting the time repetition feature include:

[0092] The autocorrelation coefficients of the pseudo-range residual and phase residual sequences are calculated using a sliding window, 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;

[0093] Hypothesis testing is performed on the autocorrelation coefficients of each window, and if the autocorrelation coefficient is greater than 0.7 within a 95% confidence interval, it is determined that the residual in the window has a time repetition dominant feature;

[0094] According to the wavelength difference of signals of different frequency bands, the autocorrelation coefficients are weighted and fused, and the weight w i The calculation formula is: Where λ i is the wavelength of the i-th frequency band signal, and n sys is the total number of frequency bands of the current satellite system;

[0095] If the same satellite is determined to be time repetition dominant in multiple consecutive windows, and its orbit height angle change range is less than 10°, the residual sequence of the satellite is marked as time repetition dominant.

[0096] The sliding window length is set to an integer multiple of the satellite orbit repetition period, and the possible value selection includes 1 times, 2 times or 3 times the period, which depends on the data sampling density and orbit 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, to ensure that the window overlap rate covers more than 90% of the data. The periodic operation of the satellite orbit causes the multipath error to be repetitive in time, and the autocorrelation coefficient (ACF) of the residual sequence in the window is used to quantify the periodicity. In practice, the ACF value is calculated for the pseudo-range residual and phase residual of each window, and if the ACF is continuously higher than the threshold (for example, 0.6-0.8), it indicates that there is significant time repetition in this period.

[0097] A confidence interval is constructed with a 95% confidence level. If the ACF value is greater than 0.7 (reference threshold) within the confidence interval, it is determined that the current window residual has a time repetition dominant feature. For multi-band signals (such as GPS L1 / L2 / L5), weights are assigned according to wavelength differences: signals with longer wavelengths (such as BDS-3 B3I, wavelength about 25 cm) are given higher weights, and signals with shorter wavelengths (such as GPS L1, wavelength about 19 cm) are given lower weights. In the weight calculation formula, the wavelength λ i As a molecule, the total wavelength sum is the denominator, ensuring that long-wavelength signals dominate the determination results. In specific operation, the ACF value is calculated in a window, hypothesis testing is performed, multi-band ACF results are fused according to weights, and finally a time-dominant determination flag is output.

[0098] The satellite overall marking mechanism relies on continuous window determination and orbit stability. If the same satellite is determined to be time repetition dominant in multiple consecutive windows (for example, 3-5 windows), and the range of its orbit elevation angle change is less than 10° (calculated through ephemeris data), the overall residual sequence of the satellite is marked as time dominant. In implementation, the satellite elevation angle standard deviation needs to be monitored. If the standard deviation is less than 1.5° (for example, geostationary satellite), the global marking is triggered directly. This design ensures that time modeling is implemented for highly stable satellites (such as geostationary satellites) in the entire period, avoiding discontinuity introduced by segmented processing. The time feature extraction mechanism significantly improves the identification accuracy of periodic multipath errors, and the multi-band weighting strategy enhances the determination robustness in complex signal environments.

[0099] The specific operation mode of the sidereal filter model includes:

[0100] According to the satellite ephemeris data, the orbit repetition period T orbs of the satellite is calculated orbs The residual sequence is segmented and aligned according to T

[0101] The mean value of the pseudorange and phase residual in each orbit period is calculated, and the abnormal segment deviating from the mean value by more than 2 times the standard deviation is removed. For the hardware delay bias of multi-band signals, the mean value sequence is corrected for inter-band bias. A sliding window with a length of 3 orbit periods is used to smooth the mean value sequence, and the weight is dynamically allocated according to the satellite elevation angle in the window. The weight allocation formula is: wherein, is the average elevation angle of the satellite in the kth orbit period.

[0102] The smoothed mean value sequence is used as the correction quantity of the sidereal filter model, and is matched to the original observation value according to the timestamp to output the corrected pseudorange and phase observation value.

[0103] Firstly, the orbit repeat period is calculated accurately according to the satellite ephemeris data, for example, the GPS satellite period is about 86164 seconds (23 hours 56 minutes 4 seconds). The residual sequence is segmented and aligned according to the period, and the alignment error is controlled within 5% of the data sampling interval, for example, 1.5 seconds of time deviation is allowed for 30 seconds of sampling data. In implementation, the linear interpolation method is used to adjust the residual timestamp, and the starting time of each period segment is strictly aligned. The principle is that the consistency of the orbit period determines the time repeatability of the multipath error, and accurate alignment can maximize the model correction effect.

[0104] The arithmetic mean of the pseudo-range and phase residual is calculated, and if a certain segment of residual sequence deviates from the mean value by more than 2 times the standard deviation (for example, the 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), the inter-band hardware delay bias is compensated: based on the pre-calibration parameters (such as E5a band bias +0.05m), the mean value sequence is corrected. The operation process includes: segmenting the residual according to the period, calculating the mean value within the segment, excluding the out-of-limit data segment, and correcting the mean value according to the frequency band compensation table.

[0105] The smoothing process uses a sliding window with a length of 3 orbit periods, and the weight of each period in the window is dynamically allocated according to the average elevation angle of the satellite. In the weight calculation formula, the sine square value of the satellite elevation angle ( ) is used as the weight factor, for example, when the elevation angle is 30°, the weight is about 0.25, and when the elevation angle is 60°, the weight is about 0.75. In implementation, the mean value sequence of the three periods in the window is weighted and averaged according to the weight, and the smoothed correction quantity sequence is output. Finally, it is matched to the original observation value according to the timestamp, for example, the correction quantity of the kth period is applied to the observation value at the corresponding time of the k+1 period. The alignment of the orbit period and the exclusion of the abnormal segment significantly improve the model accuracy, and the dynamic weighted smoothing mechanism effectively suppresses the random noise interference.

[0106] In another technical solution, the specific steps of extracting the spatial distribution feature include:

[0107] The grid is dynamically divided, and the azimuth resolution of the grid is dynamically adjusted according to the satellite elevation angle range. The specific adjustment rule is: , wherein ΔA k is the adjusted azimuth resolution, ΔA0 is the reference azimuth resolution, E k is the center value of the kth layer elevation angle interval, and the elevation angle range is divided into 5° intervals; the pseudo-range residual and phase residual in each grid are counted respectively, and the density ρ I is calculated. The density calculation formula is: , wherein N I is the number of effective residuals in the grid I, and A I is the spherical area of the grid;

[0108] For high elevation area with elevation > 60°, the residual density is corrected as , the density compensation factor , E s is the satellite elevation angle; if the average residual density of the high elevation area after correction is lower than 50% of the average residual density of the low elevation area with elevation < 30°, it is determined as spatial distribution dominant;

[0109] For all frequency bands of the same satellite system, if more than 70% of the frequency bands meet the spatial distribution dominant condition, the whole system is marked as spatial distribution dominant.

[0110] The grid is divided according to the dynamic adjustment of the azimuth resolution based on the satellite elevation angle range, and the specific rules are as follows: the reference azimuth resolution ΔA0 is usually set to 5° (optional range 3°-10°), and the azimuth resolution ΔA t of the kth elevation angle interval is calculated as ΔA t = ΔA0 / cosE t , where E t is the center value of the layer (for example, the center value of the elevation angle interval 10°-15° is 12.5°). The elevation angle range is divided at fixed intervals, and the recommended interval is 5° (for example, 0°-5°, 5°-10°, and so on to 85°-90°). As the satellite elevation angle increases, the signal coverage range decreases, and dynamically expanding the azimuth resolution can keep the grid spherical area approximately equal, avoiding the data sparseness caused by the small grid in the high elevation area. When implementing, all elevation angle levels need to be traversed, and the azimuth resolution is calculated layer by layer to generate the grid structure.

[0111] The residual density ρ t of each grid I is defined as the ratio of the effective residual number N t to the grid spherical area A t (ρ t = N t / A t ). For the high elevation area with elevation > 60°, a density compensation factor α = 1 + 0.5×(1 - Eˢ / 90°) is introduced for correction (for example, when E s =70°, α≈1.11). The condition for determining spatial distribution dominance is that the average density of the high elevation area (> 60°) after correction is lower than 50% of the average density of the low elevation area with elevation < 30°. The running process includes calculating the density of all grids, correcting the high elevation density, calculating the average density in different areas, and comparing the threshold value.

[0112] The residual density of all frequency bands of the same satellite system (such as GPS or BDS-3) is counted, and if more than 70% (optional threshold 60%-80%) of the frequency bands meet the spatial distribution dominant condition, the system data is marked as spatial dominant. When implementing, the system is divided into zones, and the proportion of each frequency band density that meets the standard is calculated. After triggering the global marking, it is automatically applied to all frequency point data of the system. This design avoids single frequency band misjudgment and improves the consistency of system-level modeling. Dynamic grid and density correction significantly improve the spatial feature representation capability, and system-level joint analysis enhances the reliability of multi-band scene judgment.

[0113] The specific operation mode of the multipath hemisphere diagram model includes:

[0114] For each dynamically divided grid, the mean of the pseudorange residual and the phase residual is calculated respectively. For multi-band signals, the residual mean is weighted and fused according to the frequency band wavelength. For grids with no valid data, the residual mean of adjacent 3 layers of elevation angle grids is used for interpolation, and the interpolated grid data is spatially smoothed using a Gaussian kernel function with a standard deviation of 2°.

[0115] When generating the correction, the grid residual mean is dynamically weighted according to the real-time elevation angle of the satellite 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 position of the satellite, and the corrected pseudorange and phase observation values are output.

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

[0117] The empty data grid uses a hierarchical interpolation strategy. For grids with no valid data (e.g. not covered by satellites), the residual mean of adjacent three layers of elevation angle grids (e.g. current layer ±1 layer) is interpolated. For example, the empty grid of elevation angle 20°-25°, the mean of 15°-20°, 25°-30° layer is interpolated. After interpolation, spatial smoothing is performed using a Gaussian kernel function, and the standard deviation of the kernel function is recommended to be set to 2° (optional range 1°-3°), and the smoothing range covers the surrounding 3x3 grid area. The implementation process includes identifying the position of the empty grid, retrieving the data of the adjacent layers, calculating the interpolated value, and applying Gaussian kernel convolution smoothing.

[0118] The grid residual contribution weight is dynamically adjusted according to the real-time elevation angle of the satellite: the elevation angle Es Weight is raised to 1.2-1.5 times at 30°, E s Weight is reduced to 0.8-1.0 times at 60°. In implementation, the weighted grid residual mean is mapped to the observation value according to the real-time azimuth and elevation angle of the satellite (for example, the azimuth 120.3° matches the 120°-125° grid). When outputting the correction value, the phase observation value needs to be multiplied by the corresponding frequency band wavelength to convert to distance units. Multi-frequency weighting fusion optimizes the spatial model precision, hierarchical interpolation and smoothing processing ensures data integrity, and dynamic weight enhances the correction effect in the low-elevation-angle area.

[0119] In another technical solution, the specific steps of extracting the trend complexity feature include:

[0120] For the residual data in each elevation-azimuth grid, linear, quadratic, and cubic polynomials are used in turn for trend fitting, and the Akaike Information Criterion (AIC) is used to select the optimal order, where k is the number of model parameters, and L is the model likelihood value. The polynomial order with the smallest AIC value is selected as the final fitting model;

[0121] The R 2 value of the selected order polynomial fitting is calculated, and the trend significance is verified by F test. The trend complexity feature determination rule is: R 2 ≥ 0.3, and the F test statistic is significant at the 95% confidence level;

[0122] For multi-band signals, the residual trend is weighted and fused according to the wavelength difference. If the trend fitting coefficient of the same satellite changes by less than 10% in more than 5 consecutive grids, the residual trend of the satellite is globally smoothed. When R 2 ≥ 0.3 and F test is significant, it is determined that the trend complexity is dominant. If multiple order models meet the conditions at the same time, the highest order model is selected.

[0123] For the residual data in each elevation-azimuth grid, linear (first order), quadratic, and cubic polynomials are used in turn for trend fitting. The order selection is based on the Akaike Information Criterion (AIC), which depends on the number of model parameters k and the likelihood value L. The smaller the AIC value, the better the model fitting effect. In implementation, the AIC values of three order models need to be calculated, for example, the AIC value of a cubic polynomial with k=10 parameters may be lower than that of a linear model with k=3 parameters. Multipath errors need to be represented by high-order polynomials, and the AIC criterion balances the fitting precision and model complexity. Alternatively, when the number of residuals is insufficient (e.g., the number of data points in the grid is less than 20), a low-order model is forced to avoid overfitting.

[0124] Fitting goodness and significance test are the core of determining the trend dominance. The determination coefficient R 2(range 0-1), the fit is considered valid when R 2 ≥ 0.3 (reference threshold 0.25-0.35). F-test is performed simultaneously to verify the significance of the trend: at 95% confidence level, the trend is considered significant if the F-statistic exceeds the critical value. The implementation process includes: 1) polynomial regression on grid residuals; 2) R 2 value calculation; 3) F-distribution table lookup to verify significance; 4) mark as trend complexity dominant when R 2 ≥ 0.3 and F-test is significant. For example, a certain grid has a cubic fit R 2 = 0.32 and significant F-value, the trend dominant flag is triggered.

[0125] For different frequency band signals (such as GPS L1 / L2 / L5), the residual trend is weighted and fused according to the wavelength ratio: the longer the wavelength, the higher the weight (for example, 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 between the adjacent grid azimuth angle coefficients c1 is less than or equal to 0.05), then the global smoothing processing is performed on the satellite residuals, and the moving average method is used to fuse the adjacent grid trend coefficients. When implementing, monitor the coefficient change rate, and re-calculate R 2 and F-value after triggering smoothing. Polynomial fitting and AIC criterion significantly improve the modeling capability of non-linear errors, the multi-band weighting strategy effectively fuses system differences, and the global smoothing mechanism enhances the spatial consistency in a large range.

[0126] The specific operation mode of the trend surface analysis multipath semisphere model includes:

[0127] For the residual data in each grid, the trend surface is fitted, and the model expression of linear, quadratic and cubic trend surfaces is as follows:

[0128]

[0129] wherein, is the multipath residual, A s is the azimuth angle, E s is the elevation angle, and c0 to c9 are the polynomial coefficients of the trend surface.

[0130] For the hardware delay bias of multi-band signals, the trend surface coefficients are compensated for frequency band differences; according to the real-time motion state of the satellite, the trend surface coefficients are dynamically weighted, and for the grids without effective data, the bilinear interpolation of the adjacent grid trend surface coefficients is used to fill the gap; the optimized trend surface coefficients are mapped to the observation values according to the real-time position of the satellite, and the corrected pseudo-range and phase data are output.

[0131] The model expression contains three levels: 1) the linear trend surface only contains the azimuth angle A s and the elevation angle Es a linear term; 2) quadratic trend surface increases , and a cross term; 3) cubic trend surface further introduces , higher order terms. Coefficients c0-c9 represent the trend surface shape, for example, c1>0 indicates the multipath error increases as the azimuth angle increases. According to the selected order, the corresponding expression is invoked, for example, six coefficients c0-c5 are calculated for the quadratic model. The polynomial surface is used to fit the spatially varying error trend, and the high order terms capture the complex reflection environment characteristics.

[0132] Frequency band difference compensation and dynamic weighting optimize coefficient accuracy. For multi-band hardware delay bias, frequency band compensation is applied to the trend surface coefficients: for example, the BDS-3 B2b frequency band adds a +0.1m offset to the c0 term. Real-time satellite motion state (such as speed >0.1° / s) triggers dynamic weighting: when moving at high speed, the weight of high order terms (such as the weight of cubic terms x 0.7) is reduced, and when moving at low speed, the weight is restored. The implementation process includes: 1) prestore the compensation table according to the frequency band; 2) obtain the satellite angular velocity in real time; 3) adjust the coefficient weight according to the motion state. For example, a GPS satellite moving at high speed uses a simplified linear model, and a stationary receiver enables a complete cubic model.

[0133] Data gap processing and real-time mapping achieve efficient correction. For grids without valid data, bilinear interpolation is used to fill in the trend surface coefficients of the adjacent 4 grids (adjacent azimuth angle on the same layer + same position on the upper and lower layers). For example, the missing grid coefficient = 0.25 x (left-up + right-up + left-down + right-down coefficients). In the correction stage, the optimized coefficients are substituted into the trend surface equation, and the real-time azimuth and elevation angles of the satellite (such as A s =45.3°, E s =28.7°) are input to output the multipath error value Δ TMHM . In implementation, a coefficient-position mapping table is established, and the correction amount is updated by seconds and superimposed on the original observation value. The multi-order trend surface structure significantly enhances the modeling capability of complex reflection scenarios, the motion state weighting effectively suppresses dynamic errors, and the bilinear interpolation guarantees the correction continuity of full airspace coverage.

[0134] In another technical solution, when the residual data simultaneously satisfies multiple determination conditions in time repeatability dominance, spatial distribution dominance, or trend complexity dominance, the multipath error correction model is selected through the following priority rules:

[0135] The model selection priority is defined as: time repeatability dominance > trend complexity dominance > spatial distribution dominance;

[0136] If the residual data meets the time repetition dominant and other characteristic dominant conditions, the sidereal filter model is preferred; if the trend complexity dominant and spatial distribution dominant are met, the trend surface analysis multi-path hemispherical model is preferred;

[0137] When multiple characteristic dominant conditions coexist, the residual variance contribution degree of each characteristic is calculated, and the contribution degree , wherein is the residual variance corresponding to the characteristic f, is the total residual variance;

[0138] The model corresponding to the characteristic with the highest contribution degree is selected, and if the difference in contribution degree is less than 10%, the combined model correction mechanism is enabled:

[0139] The multiple characteristic models with similar contribution degrees are jointly corrected, and the comprehensive correction amount Δ comb The calculation formula is: , wherein Δ SF is the correction amount of the SF model, Δ TMHM is the correction amount of the TMHM model, and Δ MHM is the correction amount of the MHM model, ω1, ω2, ω3 are weights allocated in proportion to the contribution degree C f of each model, and satisfy ;

[0140] In the correction process, the accuracy decay factor PDOP and the residual variance of the positioning result are monitored in real time, and if the PDOP value exceeds 3 or the residual variance decreases by less than 20% compared with before correction, the model with the second highest contribution degree is dynamically switched to. For conflict scenarios that trigger model switching multiple times, record the timestamp, satellite system, characteristic contribution degree, and positioning error data.

[0141] When the residual data meets multiple conditions of time repetition dominant, spatial distribution dominant, or trend complexity dominant, set a three-level priority: time repetition dominant > trend complexity dominant > spatial distribution dominant. For example, if time repetition dominant (autocorrelation coefficient > 0.8) and spatial distribution dominant (high elevation density < low elevation density 50%) are detected at the same time, the sidereal filter (SF) model is preferred. The principle is that time repetition is derived from satellite orbit periodicity, which is more stable than spatial reflection environment changes, and the error source covered is more significant. A determination flag matrix is established, and the first model that meets the conditions is triggered in order of priority.

[0142] Residual variance contribution degree calculation is used to solve the same priority conflict. The contribution degree is defined as the ratio of the residual variance corresponding to a single characteristic to the total variance . For example, when the time characteristic variance is 0.15 and the total variance is 0.30, the contribution degree is If the difference of two feature contribution is less than 10% (e.g. time contribution 42%, trend contribution 38%), the combined model is enabled: SF, TMHM, MHM model correction amount Δ SF , Δ TMHM , Δ MHM are weighted by the contribution ratio, for example, SF weight is 0.42, TMHM weight is 0.38, and MHM weight is 0.20. The running process includes decomposing each feature residual component, calculating the variance ratio, and generating a comprehensive weight by weighting.

[0143] Dynamic switching and conflict monitoring guarantee real-time reliability. During the correction process, the positioning dilution of precision (PDOP) and residual variance are continuously monitored: if PDOP > 3 (threshold range 2.5-4.0) or the residual variance is reduced by less than 20% (e.g. only reduced by 15%) compared to before correction, automatically switch to the model with the next highest contribution. For scenarios that frequently trigger switching (e.g. switch > 3 times within 10 minutes), record the timestamp, satellite system number, feature contribution, and positioning error data to form a conflict log for offline analysis. Establish a real-time monitoring thread, interrupt the current model when the trigger condition is met, and load the backup model coefficient library. The priority rules and contribution weighting significantly improve the rationality of decision-making in complex scenarios, the dynamic switching mechanism enhances the fault tolerance of the system, and the conflict log supports long-term model optimization.

[0144] It should be noted that although the above describes the steps in a specific order, it does not mean that the steps must be performed in the above specific order, in fact, some of the steps can be executed concurrently, or even in reverse order, as long as the desired function can be achieved. The number of devices and the processing scale described here are used to simplify the description of the invention, and the application, modification and change of the invention are obvious to those skilled in the art.

[0145] Although the embodiments of the present application have been disclosed as above, it is not limited to the application listed in the specification and embodiments, and can be fully applied to various fields suitable for the present application, and additional modifications can be easily realized by those skilled in the art, therefore, the present application is not limited to specific details and the figures shown and described herein, without departing from the general concept defined by the claims and equivalent scope.

Claims

1. A method of mitigating multipath errors in a multi-frequency multi-mode signal, characterized by, The method comprises the following steps: Constructing a multi-frequency multi-mode non-combination precise point positioning observation equation, the observation equation comprising a combination of multi-frequency signals from multiple global navigation satellite systems, separating the multipath error in the pseudorange observation value and the phase observation value through the observation equation, performing quality control on the separated pseudorange residual and phase residual, and eliminating abnormal residual data; Extracting and determining the characteristics of the residual data after quality control, comprising: The time repetition characteristic of the residual data is extracted by calculating the autocorrelation coefficient of the residual sequence, and when the autocorrelation coefficient is greater than 0.7, it is determined that the time repetition is dominant; the spatial distribution characteristic of the residual data is extracted by counting the residual density in the elevation and azimuth grid, and when the residual density in the grid with an elevation greater than 60° is less than 50% of the residual density in the grid with an elevation less than 30°, it is determined that the spatial distribution is dominant; the trend complexity characteristic is extracted by fitting the residual trend with a cubic polynomial and calculating the goodness of fit R 2 value and F test value, and when the R 2 value is greater than or equal to 0.3 and the F test is significant at a 95% confidence level, it is determined that the trend complexity is dominant; Selecting a multipath error correction model according to the characteristic determination result, if the time repetition is dominant, a sidereal filter SF model is selected, if the spatial distribution is dominant, a multipath hemisphere map MHM model is selected, and if the trend complexity is dominant, a trend surface analysis multipath hemisphere map TMHM model is selected; Generating the multipath error correction of the pseudorange observation value and the phase observation value through the selected correction model, correcting the original observation value, and outputting the positioning result.

2. The method of claim 1, wherein the method is implemented in a mobile station. The specific steps for constructing the multi-frequency multi-mode non-combination precise point positioning observation equation comprise: The combination of multi-frequency signals of multiple global navigation satellite systems comprises GPS three-frequency, GLONASS two-frequency, Galileo five-frequency and BDS-3 five-frequency signals, and the pseudorange and phase basic observation equation is: where, and are the raw pseudorange and phase observations of the receiver r to the i-th frequency of the satellite s, respectively, is the geometric range, and are the receiver and satellite clock biases, respectively, is the mapping function, T r is the zenith tropospheric delay; , is the i-th frequency value of the satellite s, is the first frequency value, is the slant ionospheric delay at the first frequency; , , , are the receiver and satellite hardware delay biases at the pseudorange and phase, respectively, is the ambiguity of the wavelength, and are the pseudorange and phase multipath errors, respectively; and are the observation noises; When separating the multipath error, the inter-band bias of the pseudorange and phase observation values of different frequency points is compensated based on the frequency band difference of the satellite and the receiver hardware delay bias; and the multi-frequency observation values are weighted and fused according to the signal-to-noise ratio of each frequency point signal and the satellite elevation angle.

3. The method of claim 1, wherein the method is implemented in a mobile station. The specific steps for performing quality control comprise: Based on the normal distribution assumption, the mean μ and the standard deviation σ of the pseudo-range residual and the phase residual in each grid of the elevation and azimuth are calculated respectively I and the standard deviation σ I The 3σ rule is used to identify abnormal values, and the judgment rule is: wherein, represents the i-th residual value in the grid I; For the identified outliers, the significance of their influence on the grid residual variance is verified by F-test. The specific steps are as follows: the residual data set in the grid is divided into a normal data set x I,1 and an abnormal data set x I,2 , the sample variances of the two data sets are calculated and , and the F-test statistic is constructed: If F is greater than the F α (n I,1 −1,n I,2 −1) value with α=0.05, it is determined that the outliers have a significant influence on the residual variance. For the residual data of low elevation angle area with elevation angle < 15°, the abnormality determination threshold is adjusted as: wherein the satellite motion state dynamic factor , E s is the satellite elevation angle; For the grid with abnormal values, the sliding window interpolation method with the mean of adjacent grid residuals is used to fill the data gaps, and the window size is 3x3 grid. The interpolation formula is: where N is the number of adjacent grids, μ is the mean of adjacent grid residuals. adj I,j ​​ For the corrected residual data, the residual distribution is statistically analyzed according to the satellite system, if the residual variance of one of the systems exceeds 2 times the global variance, secondary quality control of the data of the system is triggered.

4. The method of claim 1, wherein the method is implemented in a mobile station. The specific steps for extracting the time repetition characteristic comprise: The autocorrelation coefficients of the pseudorange residual and phase residual sequences are calculated by using a sliding window, 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; The autocorrelation coefficients of each window are subjected to hypothesis testing, if the autocorrelation coefficient is greater than 0.7 within the 95% confidence interval, it is determined that the residual in the window has a time repetition dominant characteristic; According to the wavelength difference of different frequency band signals, the autocorrelation coefficients are weighted and fused by frequency band, and the weight w i The calculation formula is: Wherein, λ i is the wavelength of the i-th frequency band signal, n sys is the total number of frequency bands of the current satellite system; If the same satellite is determined to have a time repetition dominant characteristic in a plurality of consecutive windows, and the range of the orbit elevation angle change is less than 10°, the residual sequence of the satellite is marked as a time repetition dominant characteristic.

5. The method of claim 4, wherein, The specific operation mode of the sidereal filter model comprises: According to the satellite ephemeris data, the orbit repetition period T is calculated orbs The residual sequence is segmented and aligned according to T orbs The alignment error is controlled to be less than 5% of the data sampling interval. The mean value of the pseudo-range and phase residual in each orbit period is calculated respectively, and the abnormal section deviating from the mean value more than 2 times of the standard deviation is eliminated; for the hardware delay deviation of the multi-band signal, the mean value sequence is corrected among the bands; the mean value sequence is smoothed by using a sliding window with a length of 3 orbit periods, and the weight in the window is dynamically distributed according to the satellite elevation angle, and the weight distribution formula is: wherein, is the average elevation angle of the satellite in the kth orbit period. The smoothed mean sequence is taken as the correction amount of the sidereal filter model, is matched to the original observation value according to the time stamp, and the corrected pseudorange and phase observation values are output.

6. The method of claim 1, wherein, The specific steps for extracting the spatial distribution characteristic comprise: Dynamic grid partitioning, dynamically adjust the azimuth resolution of the grid according to the satellite elevation angle range, the specific adjustment rule is: Where, ΔA k is the adjusted azimuth resolution, ΔA0 is the reference azimuth resolution, E k is the center value of the kth layer elevation angle interval, and the elevation angle range is divided into 5° intervals; The density ρ I of the pseudo-range residual and the phase residual in each grid is counted respectively, and the density calculation formula is: Where, N I is the number of effective residuals in the grid I, and A I is the spherical area of the grid. For high elevation angle area with elevation angle > 60°, the residual density is corrected as , the density compensation factor , E s is the satellite elevation angle; if the corrected average residual density of the high elevation angle area grid is lower than 50% of the average residual density of the low elevation angle area grid with elevation angle < 30°, it is determined that the spatial distribution is dominant; All frequency band residual densities of the same satellite system are jointly analyzed, if more than 70% of the frequency bands meet the spatial distribution dominant condition, the system is globally marked as a spatial distribution dominant.

7. The method of claim 6, wherein the method is implemented in a mobile station. The specific operation mode of the multipath hemisphere map model comprises: For each dynamic grid, the mean of the pseudo-range residual and the phase residual is calculated respectively; for multi-band signals, the residual mean is weighted and fused according to the wavelength of the frequency band; for the grid without effective data, the residual mean of the adjacent 3 layers of the height angle grid is used for interpolation, and the interpolated grid data is spatially smoothed by using the Gaussian kernel function with a standard deviation of 2°; In the generation of the correction amount, the grid residual mean is dynamically weighted according to the real-time height angle of the satellite 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 position of the satellite, and the corrected pseudo-range and phase observation values are output.

8. The method of claim 1, wherein, The specific steps of extracting the trend complexity feature include: For each residual data in the grid of elevation and azimuth, linear, quadratic and cubic polynomials are used to fit the trend, and the Akaike information criterion (AIC) is used to select the optimal order, where k is the number of model parameters, and L is the model likelihood value. The polynomial order with the minimum AIC value is selected as the final fitting model. R value of the selected order polynomial fitting is calculated, and the trend significance is verified by F test, and the trend complexity feature determination rule is: R 2 ≥ 0.3, and the F test statistic is significant at the 95% confidence level; 2 ≥ 0.3, and the F test statistic is significant at the 95% confidence level; For multi-band signals, the residual trends are weighted and fused according to the wavelength difference; if the trend fitting coefficient of the same satellite in more than 5 consecutive grids changes less than 10%, the residual trend of the satellite is globally smoothed; when R 2 ≥0.3 and the F test is significant, it is determined that the trend complexity is dominant; if multiple order models simultaneously meet the conditions, the highest order model is selected.

9. The method of claim 8, wherein the method is implemented in a mobile station. The specific operation mode of the trend surface analysis multipath hemispherical chart model includes: For the residual data in each grid, the trend surface is fitted, and the model expression of the linear, quadratic and cubic trend surface is used as: wherein is a multipath residual, A s is an azimuth angle, E s is an elevation angle, c0to c9are polynomial coefficients of a trend surface; For the hardware delay bias of multi-band signals, the trend surface coefficient is compensated for the difference between frequency bands; according to the real-time motion state of the satellite, the trend surface coefficient is dynamically weighted, and for the grid without effective data, the bilinear interpolation of the trend surface coefficient of the adjacent grid is used to fill the gap; the optimized trend surface coefficient is mapped to the observation value according to the real-time position of the satellite, and the corrected pseudo-range and phase data are output.

10. The method of claim 1, wherein, When the residual data meets multiple judgment conditions in time repeatability dominance, spatial distribution dominance or trend complexity dominance, the multipath error correction model is selected by the following priority rules: The model selection priority is defined as: time repeatability dominance > trend complexity dominance > spatial distribution dominance; If the residual data meets the time repeatability dominance and other feature dominance conditions at the same time, the constant day filtering model is preferred; if the trend complexity dominance and the spatial distribution dominance are met at the same time, the trend surface analysis multipath hemispherical chart model is preferred; When multiple feature dominant conditions coexist, the residual variance contribution degree of each feature is calculated, and the contribution degrees wherein, is the residual variance corresponding to the feature f, is the total residual variance; The model corresponding to the feature with the highest contribution degree is selected, and if the difference in contribution degree is less than 10%, the combined model correction mechanism is enabled: The multiple feature models with similar contribution degrees are jointly corrected, and a comprehensive correction amount Δ comb The calculation formula is: Wherein Δ 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 weights proportionally distributed according to the contribution degrees C f of the models, and satisfy ; During the correction process, the accuracy decay factor PDOP and the residual variance of the positioning result are monitored in real time, if the PDOP value exceeds 3 or the residual variance decreases by less than 20% compared with before correction, the model with the second highest contribution degree is dynamically switched to; for the conflict scene of multiple triggered model switching, the time stamp, satellite system, feature contribution degree and positioning error data are recorded.

Images