A GNSS multipath error modeling and compensation method, system and storage medium

Abstract

This invention discloses a GNSS multipath error modeling and compensation method, system, and electronic device, relating to the field of GNSS data processing technology. The method involves acquiring the coordinate time series output from the GNSS baseline calculation and the corresponding satellite elevation and azimuth angles at each epoch; pre-stripping the coordinate time series with time-domain periodic signals to obtain time-domain multipath components and residual sequences; converting the residual sequences into zenith distance and azimuth angles based on satellite elevation and azimuth angles to obtain spatial domain modeling input data; performing spherical harmonic function fitting modeling with the residual values ​​as observations to obtain spatial domain multipath components; generating fused multipath compensation values ​​based on the time-domain and spatial domain multipath components; and using the fused multipath compensation values ​​to compensate for multipath errors in the coordinate time series. This application organically fuses time-domain and spatial-domain multipath errors through residual concatenation, capturing the spatiotemporal coupling characteristics of multipath errors and improving the accuracy of multipath error compensation.

Description

Technical Field

[0001] This invention relates to the field of GNSS data processing technology, specifically to the field of multipath error modeling and compensation, and particularly to a dual-mode fusion GNSS multipath error modeling and compensation method, system, and storage medium based on residual concatenation mechanism. Background Technology

[0002] Global Navigation Satellite Systems (GNSS) are widely used in fields such as geodesy, deformation monitoring, and precision agriculture. In carrier phase high-precision positioning, multipath effects are one of the main error sources restricting further improvement in positioning accuracy. Multipath effects refer to the interference effect caused by the superposition of non-direct signals and direct signals at the antenna after the satellite signal is reflected, scattered, or diffracted by reflectors such as buildings, ground, and water surfaces before reaching the receiver antenna. This leads to distortion of the pseudorange correlation function and periodic disturbances in the carrier phase observations. For code pseudorange measurement, multipath causes a shift in the peak value of the correlation function, resulting in meter-level or even more severe errors. For carrier phase measurement, the multipath error amplitude is constrained by integer ambiguity, and its maximum value is usually no more than one-quarter of the carrier wavelength (approximately 5 cm for the L1 band). However, multipath residuals can make the covariance estimation for fixed ambiguity overly optimistic, increase the search space, increase the fixation failure rate, and affect the accuracy of baseline calculation. In short baseline scenarios, multipath becomes a major error source.

[0003] Existing multipath processing methods can be broadly categorized into three types: At the hardware level, choke antennas and other methods are used to reduce reflected signals, but these are ineffective against scattering and diffraction signals and are costly; at the observation combination level, ionospheric desiccation combinations are used to mitigate the multipath effect, but these are compromises and cannot fundamentally eliminate it; at the post-processing level, multipath simulation and compensation based on environmental maps are cumbersome to model and extremely sensitive to environmental changes. Furthermore, parameter estimation methods only fit the multipath error in a single direction, failing to fully utilize the two-dimensional spatial distribution characteristics formed by elevation and azimuth angles; sidereal filtering and spherical harmonic function fitting methods follow time and spatial domain routes respectively, and existing schemes often treat them as independent optional routes, using parallel or weighted methods, resulting in the inability to effectively capture the time-space coupling component in the multipath error. Moreover, spherical harmonic function fitting suffers from rank deficiency when used with double-difference observations; double-difference combinations eliminate the general reference for the multipath field, making the normal equation matrix singular. Existing schemes either manually specify a reference star or add other prior conditions, lacking universality.

[0004] Therefore, how to organically integrate time-domain and spatial-domain multipath error modeling methods to capture the spatiotemporal coupling characteristics of multipath errors and eliminate the rank deficiency problem of spherical harmonic fitting of double-difference observations is an urgent problem to be solved.

[0005] In terms of time-domain processing, sidereal-day filtering methods (including template methods and multi-period joint fitting methods) can effectively capture the periodic repetitive features of multipath errors. However, their modeling capabilities are limited to one-dimensional periodic signals on the time axis, and they cannot represent the spatial distribution features determined by the geometric position of the reflectors. Specifically, when there are multiple reflectors around the station, different satellites pass through different reflection regions at the same sidereal time position, resulting in significant differences in the amplitude and phase of their multipath errors. This spatial heterogeneity is averaged in the one-dimensional time-domain model, leading to the loss of some multipath information.

[0006] In terms of spatial domain processing, existing spatial domain modeling methods (such as the MHGM grid method and the covariance function method) can describe the distribution characteristics of multipath errors in the elevation-azimuth two-dimensional space. However, their input is usually the original double-difference residuals of the baseline solution, which simultaneously contain strong periodic signals in the time domain and distributed signals in the spatial domain. The presence of strong periodic signals will have a "stretching" effect on the spatial domain fitting surface—in the data-dense periodic peak region, the spherical harmonic surface is forcibly raised or lowered to fit the large-amplitude periodic signal, resulting in a decrease in the fitting accuracy of non-periodic spatial components and the appearance of pseudo-oscillations on the surface.

[0007] The above analysis shows that time-domain models excel at capturing periodic principal components but miss spatial heterogeneity, while spatial-domain models excel at capturing spatial distribution characteristics but are susceptible to interference from strong periodic signals. The two models are complementary in information capture, but current technology lacks a data processing mechanism for effectively transferring and fusing their information. Simply weighting and superimposing the outputs of the two independent models (using them in parallel) fails to capture the time-space coupling components in multipath errors because the two models independently model the same signal, and the residual information is not mutually utilized. Summary of the Invention

[0008] To address the aforementioned technical problems, this application discloses a GNSS multipath error modeling and compensation method, system, and electronic device; the GNSS multipath error modeling and compensation method specifically includes: Obtain the coordinate time series output by the GNSS baseline calculation, as well as the satellite elevation angle and azimuth angle corresponding to each epoch of the coordinate time series; The coordinate time series is subjected to time domain periodic signal pre-stripping processing to obtain time domain multipath components and residual sequences. The residual sequences are the sequences remaining after stripping periodic signals and trend terms from the coordinate time series. The residuals of each epoch in the residual sequence are converted into zenith distance and azimuth according to the satellite elevation angle and azimuth angle of the corresponding epoch, to obtain spatial domain modeling input data composed of residual values, zenith distance and azimuth angle; wherein, the zenith distance is determined according to the satellite elevation angle and the azimuth angle is determined according to the satellite azimuth angle. Using the residual values ​​in the spatial domain modeling input data as the observed values ​​and the zenith distance and azimuth angle as independent variables, spherical harmonic function fitting modeling is performed to obtain the spatial domain multipath components; Based on the obtained spatial domain multipath components, a secondary residual is calculated on the residual sequence. The residual stellar diurnal component in the secondary residual is fed back to the time domain to update the time domain multipath components and the spherical harmonic function fitting model is re-executed. Cross-domain iteration is performed until convergence. The sum of the converged time domain multipath components and spatial domain multipath components is used as the fused multipath compensation value. The fused multipath compensation value is used to perform multipath error compensation on the coordinate time series.

[0009] Preferably, obtaining the coordinate time series output from the GNSS baseline solution includes data preprocessing and quality assessment: Perform consistency checks on the input data, examining the temporal continuity of the data, sampling intervals, and observation quality indicators; Filter low elevation angle satellite observations based on the set elevation angle cutoff angle; Remove gross errors and outlier observations, and check the rationality of the station coordinates.

[0010] Preferably, the step of performing time-domain periodic signal pre-stripping processing on the coordinate time series to obtain time-domain multipath components and residual sequences includes: Convert the timestamps of each epoch in the coordinate time series to the sidereal time reference time of the corresponding GNSS system, so that the multi-day data are aligned on the sidereal day time axis. For each component of the coordinate time series, the following joint observation equation is established for multi-period joint fitting: The formula is as follows:

[0011] Where Y(t) is the coordinate time series, and t is the continuous cumulative time. Let k be the sidereal day period of the GNSS system, and k be the harmonic order. For constant terms, The slope of the linear trend. , Here, N represents the k-th harmonic coefficient, and N is the highest truncation order in the time domain. It is random noise; The parameters of the joint observation equation are obtained by least squares solution, and the time-domain multipath component is extracted from the fitting results. The formula is:

[0012] and the residual sequence The formula is: .

[0013] Preferably, the highest truncation order N in the time domain is determined by a successive test method: The fitting residuals for each coordinate component under different harmonic orders are calculated. Starting from the first harmonic, the order is gradually increased. When the residual decreases gradually after increasing a certain order, the process stops. The order at this point is taken as the highest truncation order N in the time domain.

[0014] Preferably, the step of using the residual values ​​from the spatial domain modeling input data as observed values ​​and the zenith distance and azimuth angle as independent variables to perform spherical harmonic function fitting modeling to obtain spatial domain multipath components specifically involves: The zenith distance is obtained by converting the satellite elevation angle using z = 90° - E, and the azimuth angle is obtained directly from the satellite azimuth angle, where E is the satellite elevation angle. The spherical harmonic function fitting model adopts the following spherical harmonic function expansion model. The formula is:

[0015] Where z is the zenith distance and α is the azimuth angle. To fully normalize the Legendre accompanying polynomial, , Let N be the spherical harmonic coefficients to be estimated, N be the highest cutoff order of the spherical harmonics, and m be the order within the order, with m ranging from 0 to n. The summation starts from n=1 to N, and zero-order terms n=0 are removed.

[0016] Preferably, the spherical harmonic function fitting modeling further includes introducing zenith virtual constraints and Tikhonov regularization: Construct the zenith virtual constraint equations. At the zenith direction (i.e., at zenith distance z=0), based on the Legendre polynomial property, retain only the harmonic terms where m=0 to construct the zenith virtual constraint equations. The formula is:

[0017] The zenith virtual constraint equation is given extremely high weight. To establish the zenith direction as the absolute zero reference for multiple paths; Tikhonov regularization is introduced to construct a regularization matrix R, which is a diagonal matrix with n (n+1) diagonal elements, and the regularization parameter λ is set. The spherical harmonic coefficient vector X is solved using the joint equation method, as follows:

[0018] Wherein, B1 is the design matrix for fitting the spherical harmonic function, P1 is the observation weight matrix, L1 is the observation vector, and B2 is the design matrix for the zenith virtual constraint. The joint equation is obtained by minimizing the objective function, and the formula is as follows: , where V1=B1·X-L1.

[0019] Preferably, the process of solving the spherical harmonic coefficient vector X also includes gross error detection: After each round of solving, the fitting residual is calculated, and the gross error observations are weighted according to the post-hoc unit weight mean square error, updating the observation weight matrix. Then, the combined equations are solved again, and the above process is repeated until the solution converges.

[0020] Preferably, the method further includes outputting model files and quality control information: Output a time-domain compensation file, a spatial-domain compensation file, and a fusion compensation configuration file, wherein the time-domain compensation file records the mapping between the sidereal time position and the time-domain multipath component, the spatial-domain compensation file records the spherical harmonic coefficients and truncation order, and the fusion compensation configuration file records the fusion method and quality control parameters; The output fitted residual sequence and standard deviation are used as quality control indicators. When the fused residual exceeds the preset threshold, an abnormal alarm is output.

[0021] A GNSS multipath error modeling and compensation system includes an acquisition module, a time-domain processing module, a spatial-domain processing module, a fusion module, and a compensation module, specifically: The acquisition module acquires the coordinate time series output by the GNSS baseline calculation, as well as the satellite elevation angle and azimuth angle corresponding to each epoch of the coordinate time series; The time domain processing module performs time domain periodic signal pre-stripping processing on the coordinate time series to obtain time domain multipath components and residual sequences. The residual sequences are the sequences remaining after stripping periodic signals and trend terms from the coordinate time series. The spatial domain processing module converts the residuals of each epoch in the residual sequence into zenith distance and azimuth according to the satellite elevation angle and azimuth of the corresponding epoch, thereby obtaining spatial domain modeling input data composed of residual values, zenith distance and azimuth. The module then uses the residual values ​​in the spatial domain modeling input data as observation values ​​and the zenith distance and azimuth as independent variables to perform spherical harmonic function fitting modeling, thereby obtaining spatial domain multipath components. The fusion module generates a fused multipath compensation value based on the time-domain multipath component and the spatial-domain multipath component. The spatial-domain multipath component is obtained by fitting the residual sequence with a spherical harmonic function, so that the time-domain processing and the spatial-domain processing form a unidirectional data flow fusion relationship through residual concatenation. Before generating the fused multipath compensation value, a secondary residual is calculated on the residual sequence based on the initially obtained spatial-domain multipath component. The residual sidereal periodic component in the secondary residual is fed back to the time domain to update the time-domain multipath component and the spherical harmonic function fitting model is re-executed. Cross-domain iterative refinement is performed until convergence. The sum of the converged time-domain multipath component and the spatial-domain multipath component is used as the fused multipath compensation value. The compensation module uses the fused multipath compensation value to perform multipath error compensation on the coordinate time series.

[0022] An electronic device includes: a memory and a processor, the memory and the processor being communicatively connected to each other, the memory storing computer instructions, and the processor executing the computer instructions to perform any of the GNSS multipath error modeling and compensation methods described above.

[0023] Compared with the prior art, the technical solution of this application has the following technical effects: This invention uses sidereal time alignment combined with multi-period joint fitting to process coordinate time series, which can completely separate periodic multipath components and residual sequences within the observation data, extract periodic multipath interference in a unified manner, simplify the data basis for subsequent spatial domain modeling, reduce signal interference in the spatial fitting process, and provide regular input data for multipath error processing in a domain-specific collaborative manner.

[0024] This invention constructs a spherical harmonic function spatial model based on zenith distance and azimuth angle, and performs the solution with zenith virtual constraints and regularization processing. It can establish a multipath distribution mapping covering the entire satellite observation space domain, make up for the spatial error characteristics missing in the time domain processing, and at the same time ensure the stability of the solution process and eliminate the equation singularity problem in spatial modeling.

[0025] This invention employs a residual concatenation combined with cross-domain iteration processing logic to achieve interconnection and linkage between the time domain and spatial domain modeling processes. It feeds back and corrects the error components that are not fully extracted by the two types of models, improves the ability to capture the spatiotemporal coupling characteristics of multi-path errors, and allows the error components fitted by the two types of models to complement and improve each other.

[0026] This invention outputs standardized compensation files and supporting quality control information, corrects coordinate timing based on the multipath compensation values ​​obtained through fusion, fully realizes integrated multipath error correction, and simultaneously retains model parameters and error assessment basis, facilitating subsequent data verification and dynamic model maintenance, and is suitable for long-term continuous GNSS monitoring operations.

[0027] The above description is only an overview of the technical solution of this application. In order to better understand the technical means of this application and implement it in accordance with the contents of the specification, and to make the above and other objects, features and advantages of this application more obvious and understandable, the preferred embodiments of this application are described in detail below with reference to the accompanying drawings.

[0028] The above and other objects, advantages and features of this application will become more apparent to those skilled in the art from the following detailed description of specific embodiments in conjunction with the accompanying drawings. Attached Figure Description

[0029] To more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings used in the following description of the embodiments will be briefly introduced. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0030] Figure 1 This is a system module framework diagram provided in an embodiment of the present invention; Figure 2 This is a flowchart illustrating a GNSS multipath error modeling and compensation method provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of a time-domain component pre-stripping process provided in an embodiment of the present invention; Figure 4 This is a schematic diagram of a sky map fitted with a spherical harmonic function according to an embodiment of the present invention. Figure 5 This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present invention. Detailed Implementation

[0031] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. In the following description, specific details such as specific configurations and components are provided merely to help fully understand the embodiments of this application. Therefore, those skilled in the art should understand that various changes and modifications can be made to the embodiments described herein without departing from the scope and spirit of this application. In addition, for clarity and brevity, descriptions of known functions and structures are omitted in the embodiments.

[0032] It should be understood that the phrase "an embodiment" or "this embodiment" throughout the specification means that a specific feature, structure, or characteristic related to the embodiment is included in at least one embodiment of this application. Therefore, "an embodiment" or "this embodiment" appearing throughout the specification does not necessarily refer to the same embodiment. Furthermore, these specific features, structures, or characteristics can be combined in any suitable manner in one or more embodiments.

[0033] Furthermore, reference numerals and / or letters may be repeated in different examples within this application. Such repetition is for the purpose of simplification and clarity and does not in itself indicate a relationship between the various embodiments and / or settings discussed.

[0034] In this article, the term "and / or" is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can mean: A exists alone, B exists alone, and A and B exist simultaneously. The term " / and" in this article describes another type of relationship between related objects, indicating that two relationships can exist. For example, A / and B can mean: A exists alone, and A and B exist alone. In addition, the character " / " in this article generally indicates that the related objects before and after it are in an "or" relationship.

[0035] In this article, the term "at least one" is merely a description of the relationship between related objects, indicating that there can be three relationships. For example, "at least one of A and B" can mean: A exists alone, A and B exist simultaneously, or B exists alone.

[0036] It should also be noted that, in this document, relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion. like Figure 1 As shown, Figure 1 This is a system module framework diagram provided in the embodiments of this application, such as... Figure 1 As shown, the framework adopts a modular hierarchical architecture, with input interface, data management, preprocessing and quality control, dual-branch parallel modeling, normal equation solving and result output functional units arranged sequentially from top to bottom along the data flow direction, forming a series-parallel hybrid pipeline structure.

[0037] On the data input side, the system receives three types of data: coordinate time series, double-difference residual files, and multi-system observation data conforming to the RINEX 2 / 3 standard.

[0038] Input data enters the data management unit, which includes two functions: memory pool allocation and file reading. These functions are used to schedule and load underlying data resources, providing a unified data foundation for subsequent processing.

[0039] After data management, the observation data is transferred to the preprocessing and quality control unit, which sequentially performs three processes: gross error detection, elevation angle filtering, and consistency detection, in order to eliminate abnormal observations, filter out low-quality data, and complete data consistency verification.

[0040] The clean data that has passed quality control further branches into two parallel modeling branches, forming the characteristic dual-branch structure of this framework.

[0041] The left branch is the sidereal-sunday filtering modeling unit, which includes three processes: multi-period joint fitting, construction of the design matrix, and least squares solution, used to achieve multipath error modeling in the time domain. The right branch is the spherical harmonic function fitting modeling unit, which includes three processes: zenith constraint, regularization, and iterative weight reduction, used to achieve multipath error modeling in the spatial domain.

[0042] The modeling results of the two branches are submitted to the normal equation solving unit for joint processing. This unit includes three functions: constructing the joint normal equation, matrix inversion, and solving the variance and covariance, which are used to perform unified calculation of the parameters of the two branches and evaluate the accuracy.

[0043] The output unit generates two types of correction products: time-domain compensation files and spatial-domain compensation files, and outputs quality control statistics, thus providing usable products in both the time and space domains for precise correction of multipath errors. This framework integrates the two types of models into a unified solution system through a fusion approach of parallel time-domain and spatial-domain modeling and joint solution of normal equations.

[0044] The GNSS multipath error modeling and compensation method provided in this application is designed for static observation scenarios and can be applied to GNSS continuous operation monitoring projects such as reservoir dam deformation monitoring, slope and landslide surface displacement monitoring, and engineering control networks. The input data consists of the NEU coordinate time series output from the GNSS baseline calculation and the elevation and azimuth angles of visible satellites at each epoch. After dual-mode fusion modeling involving pre-stripping of the time-domain periodic signal and fitting of the spatial-domain spherical harmonic function, the fused multipath compensation value is output and used to compensate the coordinate time series. The method can be executed by electronic devices, and the software system GDAmm supports local private deployment and offline operation.

[0045] It should be noted that multipath error has the following characteristics: high environmental dependence - the multipath effect is most severe in densely built-up urban canyons, relatively smaller but not non-existent in open environments, and exhibits obvious periodicity over water due to water surface reflection; significant spatiotemporal correlation - adjacent stations may experience similar but not identical multipath environments, and the same station may experience similar multipath conditions at the same time of day; complex frequency characteristics - including low-frequency components (caused by changes in satellite geometry, with periods up to several hours), mid-frequency components (caused by the movement of reflectors around the antenna, such as water surface ripples), and high-frequency components (caused by dynamic scatterers such as vehicles and pedestrians); nonlinear superposition characteristics - the total multipath error of multiple reflected signals is a nonlinear superposition of the contributions of each reflected signal, which may result in constructive or destructive interference.

[0046] Multipath errors exhibit various characteristics over time: diurnal periodicity is the most significant temporal feature, stemming from the sidereal repeatability of GNSS satellite orbits, which forms the physical basis of sidereal filtering methods; semi-diurnal periodicity and higher-order harmonics arise from the non-sinusoidal changes in zenith and azimuth angles caused by the ellipticity of satellite orbits, the beat frequency effect resulting from the superposition of multipath contributions from different satellites, and the inherent periodicity of ground reflection characteristics; slow trend changes are caused by seasonal changes in the environment surrounding the station (such as seasonal changes in vegetation and construction activities); and random rapid changes are caused by dynamic scatterers (vehicle traffic, pedestrian movement, and water surface waves), which do not exhibit sidereal repeatability and cannot be corrected using periodic-based methods.

[0047] Multipath error exhibits spatial characteristics related to satellite geometry and reflector location: azimuth dependence reflects the influence of reflector location distribution around the station, with multipath errors typically more severe in densely built-up areas; elevation dependence – satellite signals at elevation angles travel longer low-altitude paths with more reflectors, making multipath effects more pronounced than those at higher elevation angles; spatial continuity – the multipath environment of the same station remains relatively stable over a short period, and multipath errors show a certain spatial correlation when adjacent stations are close together and have similar environments.

[0048] Based on the above description, this invention proposes a GNSS multipath error modeling and compensation method, such as... Figure 2 As shown, it includes the following steps: S201, Obtain the coordinate time series output by the GNSS baseline calculation, as well as the satellite elevation angle and azimuth angle corresponding to each epoch of the coordinate time series.

[0049] Specifically, the GNSS baseline solution outputs the coordinate time series Y(t) of the three NEU components (North, East, and Zenith), and also outputs the visible satellite elevation angle E and azimuth angle A corresponding to each epoch. The coordinate time series contains continuous observation data over multiple days, and the sampling interval can be customized.

[0050] In one possible implementation, S201 also includes data preprocessing and quality assessment: performing consistency checks on the input data, checking the temporal continuity of the data, sampling intervals, and observation quality indicators; filtering low-elevation-angle satellite observations based on a set elevation angle cutoff angle (usually 10° to 15°), because low-elevation-angle signals are more susceptible to multipath interference and have a lower signal-to-noise ratio; removing gross errors and anomalous observations, and checking the rationality of station coordinates to ensure that the time series of the NEU component fluctuates within a reasonable range.

[0051] S202, perform time-domain periodic signal pre-stripping processing on the coordinate time series to obtain time-domain multipath components and residual sequences. The residual sequences are the sequences remaining after stripping periodic signals and trend terms from the coordinate time series.

[0052] For details, see Figure 3 , Figure 3 This is a schematic diagram of a time-domain component pre-stripping process provided by an embodiment of the present invention; the periodicity of multipath error originates from the periodicity of satellite orbital motion, exhibiting a clear periodicity on the sidereal day time scale. The timestamps of each epoch in the coordinate time series are converted to the corresponding sidereal time reference time of the GNSS system, aligning the multi-day data on the sidereal day time axis, where Tsid is the sidereal day period length of the corresponding GNSS system (86164 seconds for GPS, 8 sidereal days for GLONASS regression, 10 sidereal days for Galileo, and 7 sidereal days for BDS MEO). For each component of the coordinate time series, the following joint observation equation is established for multi-period joint fitting:

[0053] Where Y(t) is the coordinate time series, and t is the continuous cumulative time. Here, k represents the sidereal day period length corresponding to the GNSS system, and k is the harmonic order. For constant terms, The slope of the linear trend. , denoted as the kth harmonic coefficient, and N is the highest cutoff order (typically 1 to 3). This is random noise. To ensure the purity of the periodic function, the trend term (constant term) is... With linear drift term The multipath periodic function is included in the same observation equation for joint estimation, avoiding the propagation of truncation error caused by "detrending before fitting". The multipath DC component has been incorporated into the constant term c0.

[0054] It should be noted that there are two implementation paths for the sidereal daily filtering method: the template method directly uses the calculated values ​​from the same time of the previous day or several days as the estimation template for multipath error. This method is simple and intuitive but sensitive to noise and outliers, and requires the multipath environment to remain strictly unchanged between two days; changes in ambient temperature will alter the reflection characteristics and cause the template to fail. The method based on periodic function fitting, according to the Fourier series expansion theorem, any periodic function satisfying the Dirichlet conditions can be expanded into a superposition of a series of harmonic components. By jointly fitting the periodic function parameters using observation data from multiple periods, a more robust and accurate multipath estimate is obtained. This invention adopts the latter approach.

[0055] Due to the constant term The DC component A0 of the multipath model and the DC component A0 are mathematically completely collinear (their superposition remains a constant). To make the normal equations invertible, they must be merged—discarding A0 in the multipath model and absorbing its effect into the constant parameter of the trend term. middle.

[0056] The sidereal-day filtering method is applicable under the following conditions: the multipath reflection environment around the station remains stable within consecutive sidereal days, the satellite orbital period repeatability is sufficiently accurate, and the influence of dynamic scatterers can be ignored or smoothly suppressed. When significant environmental changes occur (such as the addition of new buildings, vegetation growth, or construction activities), the time repeatability assumption fails, and the filtering effect decreases.

[0057] It should be noted that the orbital periods and return periods of different GNSS systems are different: GPS system has an orbital period of approximately 11h 58min and a return period of 1 sidereal day; GLONASS has an orbital period of approximately 11h 15min 44s and a return period of 8 sidereal days; Galileo has an orbital period of approximately 14h 4min 45s and a return period of 10 sidereal days; BDS MEO has an orbital period of approximately 12h 55min 23s and a return period of 7 sidereal days; BDS GEO / IGSO has an orbital period of 24h and is stationary relative to the ground with no return period.

[0058] The parameters are obtained by solving using least squares. And extract the time-domain multipath component from the fitting results. =

[0059] The formula for the residual sequence is: The physical meaning of the residual sequence R(t) is: the multipath error component remaining after removing the strong periodic signal and linear trend in the time domain, which is mainly composed of multipath and random noise distributed in the spatial domain; N is the optimal truncation order in the time domain.

[0060] In one possible implementation, the highest cutoff order N is determined by a step-by-step verification method: the fitting residuals for each coordinate component are calculated at different harmonic orders, starting with the first harmonic and gradually increasing the order. The process stops when the residual decreases at a certain order and the increase in order becomes gradual; this order is then taken as the highest cutoff order N. For most static monitoring stations, harmonics of orders 1-3 are sufficient to cover the main multipath periodic components.

[0061] S203, convert the residuals of each epoch in the residual sequence into zenith distance and azimuth according to the satellite elevation angle and azimuth angle of the corresponding epoch, to obtain spatial domain modeling input data composed of residual values, zenith distance and azimuth angle.

[0062] Specifically, for static stations, multipath error depends only on the signal incident direction (elevation angle E, azimuth angle A). The residual value R(t) of each epoch is matched with the elevation angle Ei and azimuth angle Ai of all visible satellites at that epoch, transforming them to obtain the zenith distance zi = 90° - Ei and azimuth angle αi = Ai, forming the input data for space domain modeling. At each epoch t, ​​R(t) reflects the combined contribution of the multipath effect of all visible satellites at that epoch.

[0063] In one possible implementation, the zenith distance is obtained by converting the satellite elevation angle using z = 90° - E, and the azimuth angle is obtained directly from the satellite azimuth angle, where E is the satellite elevation angle.

[0064] S204, using the residual values ​​in the spatial domain modeling input data as the observed values ​​and the zenith distance and azimuth angle as independent variables, perform spherical harmonic function fitting modeling to obtain the spatial domain multipath components.

[0065] Specifically, a fully normalized spherical harmonic function expansion model is adopted:

[0066] Spherical harmonic functions are orthogonal and complete basis functions defined on a sphere. Considering the satellite sky view as a hemisphere centered at the station, the multipath error is a scalar field on this hemisphere, which can be fitted with spherical harmonic functions for high-precision 3D spatial unfolding and modeling. Spherical harmonic function fitting modeling is based on spatial position modeling, unaffected by daily satellite time offsets, and completely avoids the problem of multi-system periodic differences. Introducing normalization is crucial; otherwise, the coefficients of higher-order terms will be extremely small, leading to severe ill-conditioned normal equations. Here, z is the zenith distance, and α is the azimuth angle. To fully normalize the Legendre adjoint polynomial, Cnm and Snm are the spherical harmonic coefficients to be estimated, N is the highest truncation order of the spherical harmonics (generally taken as 2 to 4), and m is the degree within the order, with m ranging from 0 to n. The summation starts from n=1 to N, and the zero-order term n=0 is removed (the zero-order term represents the deviation of the full spatial constant, which has been absorbed by the ambiguity parameter in the double difference). The highest truncation order N of the spherical harmonics and the highest truncation order N of the time-domain harmonics are two independent parameters, and they can take different values ​​(the time-domain N is generally 1 to 3, and the spherical harmonic N is generally 2 to 4).

[0067] The design matrix B1 for fitting the spherical harmonic function is constructed as follows: For each epoch t, ​​the zenith distance zi(t) and azimuth angle αi(t) of all visible satellites in that epoch are substituted into the spherical harmonic basis function. For each triple consisting of the residual value, zenith distance, and azimuth angle, a row vector of the design matrix B1 is constructed. The observation vector L1 consists of the residual values ​​R(t) of each epoch. The observation weight matrix P1 is based on elevation angle weighting and can use a sine function model, assigning higher weights to satellites with high elevation angles and lower weights to satellites with low elevation angles. The dimension of the coefficient vector X is... (Because the n=0 term was removed), the regularization matrix R has a dimension of P×P, and the value of its diagonal element Rn,m depends only on the order n and is independent of the series m, i.e. Rn,m=n(n+1).

[0068] The spherical harmonic function fitting method is applicable when the multipath error is continuous and separable in the satellite orientation space and can be accurately described by a finite-order spherical harmonic function. When the spatial distribution of reflectors is relatively uniform and regular, the spatial variation of multipath error is relatively smooth, and the spherical harmonic function fitting effect is good. However, when there are strong scatterers or the spatial distribution is discontinuous, the multipath error may show abrupt changes or irregular distribution, affecting the fitting effect.

[0069] In one possible implementation, see Figure 4 , Figure 4 This is a schematic diagram of a sky map fitted with a spherical harmonic function according to an embodiment of the present invention; the spherical harmonic function fitting modeling also includes introducing a zenith virtual constraint and Tikhonov regularization: constructing a zenith virtual constraint, at the zenith direction, i.e., at the zenith distance z=0, based on the properties of Legendre polynomials, only the harmonic terms with m=0 are retained, and constructing the zenith virtual constraint equation:

[0070] Assign extremely high weights to P2 (10 8 ~10 10The zenith direction is established as the absolute zero reference for multipath propagation. Since the zenith direction signal exhibits almost no multipath effect, this constraint mathematically establishes an absolute zero reference for the sphere. Tikhonov regularization is introduced—to suppress high-frequency oscillations caused by holes in the sky map data, a penalty term for the surface roughness of the spherical harmonic coefficients is added to the objective function. A diagonal regularization matrix R (with diagonal elements of n(n+1)) is constructed, and a regularization parameter λ is set (ranging from 0.01 to 10, with smaller values ​​for good data quality and larger values ​​for poor data quality). The spherical harmonic coefficient vector X is solved using the following joint equation:

[0071] Where B1 is the design matrix for fitting the spherical harmonic function, P1 is the observation weight matrix, L1 is the observation vector, and B2 is the design matrix for the zenith virtual constraint. The joint normal equation is derived from minimizing the objective function. We obtain V1 = B1·X - L1.

[0072] In one possible implementation, the process of solving the spherical harmonic coefficient vector X also includes gross error detection: after each round of solving, the fitting residual is calculated, the gross error observation is identified based on the post-hoc unit weight error and its weight is reduced, the observation weight matrix P1 is updated and the joint normal equation is solved again, and the above process is repeated until the solution converges, and finally the spherical harmonic coefficient vector X is output.

[0073] It should be noted that when using double-difference observation residuals as the spatial domain modeling object, the construction method of design matrix B1 is adjusted accordingly: for each double-difference observation, the Legendre function values ​​at the non-reference star and reference star positions are filled into the corresponding two columns of the design matrix, and the two columns are subtracted to obtain the row vector of the design matrix corresponding to that double-difference observation. At this time, the double-difference residuals... This involves multipath coupling between a non-reference star and a reference star. Since the double-difference combination eliminates the general reference of the multipath field, the normal equations exhibit rank deficiency. In this scenario, the aforementioned zenith virtual constraint not only serves as a physical prior but also mathematically eliminates the rank deficiency and establishes an absolute zero reference for the sphere. Although the solved spherical harmonic coefficients are input with double differences, their physical meaning is already a single-star multipath spatial field with the zenith as the absolute zero reference.

[0074] S205: After completing the sequential process of time domain stripping, residual spatial domain modeling, and fusion compensation value construction for the first time, this step uses a reverse feedback mechanism to refine the parameters between the two domains, upgrading the unidirectional data flow to bidirectional coupling.

[0075] The quadratic residual calculation, based on the initial spatial domain modeling results, calculates the quadratic residual sequence:

[0076] in R(t) represents the multipath estimate obtained from the initial spatial domain modeling, and R2(t) represents the first-order residual generated in step S204. The physical meaning of R2(t) is the remaining components that were not captured after fitting the spatial domain model.

[0077] Residual periodic backtesting is performed on R2(t) to check for any remaining sidereal periodic components: Calculate R2(t) at the sidereal daily period frequency Power spectral density at and with a residual Power spectral density at the same frequency Comparison. Define the cycle residual rate:

[0078] when When the value is ≥ 0.1, it is considered that a significant periodic component still remains in R2(t), indicating that the initial time-domain fitting is insufficient (some periodic signals were not captured). Captured and absorbed into the spherical harmonic coefficients by the spatial domain model, resulting in and There are non-orthogonal couplings between them.

[0079] Time-domain parameter feedback refinement, when the cycle residual rate When ≥ 0.1, perform the following feedback refinement: Extract the residual periodic component from R2(t): Perform a bandpass filter on R2(t) at the sidereal daily periodic frequency to extract the residual periodic signal. ; Will Superimposed on the initial time-domain estimate, updating the time-domain parameters:

[0080] Equivalently, for the parameters of the joint observation equation { , Corrections are needed:

[0081]

[0082] Where Δa k Δb k Depend on The Fourier coefficients are determined.

[0083] Regenerate residuals based on the updated time-domain parameters:

[0084] With R （2） Using (t) as input, re-perform the spatial domain spherical harmonic function fitting to obtain the updated result. (z, α).

[0085] Convergence criterion: Define cross-domain iteration convergence conditions. Convergence is determined when both of the following conditions are met: Rate of change of time domain parameters:

[0086] Where θ is the time-domain parameter vector {a k , b k}, where Δθ is the parameter change between two consecutive iterations. The convergence threshold (empirically taken as 0.01). Secondary residual variance reduction rate:

[0087] The variance reduction threshold (empirically taken as 0.05) indicates that the spatial domain model can no longer extract more effective information from the updated residuals.

[0088] The maximum number of iterations is limited to 3 to prevent overfitting. If convergence is not achieved after exceeding the maximum number of iterations, the result of the last iteration is used, and a convergence status indicator is output.

[0089] The final fusion after iteration termination: After iteration convergence, the fusion compensation value is constructed using the final time domain parameters and spatial domain parameters.

[0090] Where j is the iteration number at convergence. At this point, re-performing the orthogonality verification is expected to significantly reduce the cross-correlation coefficient |ρ| (because the iterative refinement has incorporated the periodic component into the time domain model, while the spatial domain model retains only the pure spatial correlation component).

[0091] The technical significance of cross-domain iterative refinement: This mechanism allows the parameters of the time-domain model and the spatial-domain model to no longer be determined independently, but rather refined through residual feedback. The initial concatenation may result in insufficient time-domain stripping, leading to the leakage of periodic components into the spatial-domain model. Iterative refinement recovers these leaked periodic components into the time-domain model, enabling each model to achieve optimal fitting within its preferred feature subspace. This bidirectional coupling mechanism is impossible to achieve with existing techniques that only involve single-domain processing.

[0092] S206, Generate a fused multipath compensation value based on the time-domain multipath component and the spatial-domain multipath component, wherein the fused multipath compensation value is the sum of the time-domain multipath component and the spatial-domain multipath component.

[0093] Specifically, for any epoch t to be processed and visible satellite i, the time domain compensation value is... Based on the current sidereal time position, the spatial domain compensation value is calculated from the time domain model. (z,α) is calculated from the spherical harmonic model based on the current satellite elevation and azimuth angles; fused compensation value. (t,z,α)= + (z,α) are superimposed according to the residual series relationship, and used as the final multipath correction number.

[0094] It should be noted that the aforementioned residual concatenation mechanism means that the time-domain processing does not directly output the final compensation value, but instead outputs the residual R(t) and the time-domain components. The spatial domain spherical harmonic function is fitted to model the spatial distribution on R(t), and the spatial domain components are output. This fusion method differs from the parallel approach of weighted superposition of two independent estimates. Instead, it achieves an organic fusion through residual concatenation—the time domain removes strong periodic signals, while the spatial domain fits spatially correlated components onto the residuals, capturing the spatiotemporal coupling characteristics of multipath errors. Early removal of strong periodic principal components in the time domain prevents the spherical harmonic surface from being stretched by large periodic signals, while the spatial domain performs fine modeling on low-amplitude residuals.

[0095] S207, use the fused multipath compensation value to perform multipath error compensation on the coordinate time series.

[0096] Specifically, the fusion compensation value As the final multipath correction, it is subtracted from the original NEU coordinate time series to complete the dual-mode fusion compensation.

[0097] In one possible implementation, the method further includes outputting model files and quality control information: outputting a time-domain compensation file, a spatial-domain compensation file, and a fusion compensation configuration file. The time-domain compensation file records the sidereal position and the time-domain multipath components. The mapping is as follows: the spatial domain compensation file records the spherical harmonic coefficients and truncation order; the fusion compensation configuration file records the fusion method and quality control parameters. The output fitted residual sequence and standard deviation are used as quality control indicators to evaluate the quality of the fusion model. When the fused residual exceeds a preset threshold (usually set to 3 times the unit weight mean square error), an anomaly alarm is output, prompting the user to check the data quality or adjust the model parameters.

[0098] It should be noted that time-domain and spatial-domain processing complement each other at the application level: the sidereal-diurnal filtering method is suitable for handling static multipath effects with obvious sidereal-diurnal periodic characteristics, while the spherical harmonic function fitting method is suitable for handling the spatial distribution characteristics of multipath effects and is applicable to single-day data analysis or scenarios with rapid environmental changes. Users can flexibly choose to use the time-domain method independently, the spatial-domain method independently, or the dual-mode fusion method according to data characteristics and application requirements to maximize the advantages of both methods. The system defaults to the fusion mode to obtain better compensation results.

[0099] Based on the description of the above embodiment of the GNSS multipath error modeling and compensation method, this invention also discloses a GNSS multipath error modeling and compensation system. The system includes an acquisition module, a time-domain processing module, a spatial-domain processing module, a fusion module, and a compensation module.

[0100] Specifically, the acquisition module is used to acquire the coordinate time series output by the GNSS baseline calculation, as well as the satellite elevation angle and azimuth angle corresponding to each epoch of the coordinate time series; the time domain processing module is used to perform time domain periodic signal pre-stripping processing on the coordinate time series to obtain time domain multipath components and residual sequences, wherein the residual sequence is the sequence remaining after stripping the periodic signal and trend term from the coordinate time series; the spatial domain processing module is used to convert the residuals of each epoch in the residual sequence into zenith distance and azimuth angle according to the satellite elevation angle and azimuth angle of the corresponding epoch, to obtain spatial domain modeling input data composed of residual values, zenith distance, and azimuth angle, and to perform spherical harmonic function fitting modeling with the residual values ​​in the spatial domain modeling input data as observation values ​​and the zenith distance and azimuth angle as independent variables to obtain spatial domain multipath components; the fusion module is used to generate fused multipath compensation values ​​based on the time domain multipath components and the spatial domain multipath components; the compensation module is used to perform multipath error compensation on the coordinate time series using the fused multipath compensation values.

[0101] In one possible implementation, the time-domain processing module converts the timestamps of each epoch in the coordinate time series into the sidereal time reference time of the corresponding GNSS system, establishes a joint observation equation for multi-period joint fitting, and uses least squares to solve and extract the time-domain multipath components and residual sequences.

[0102] In one possible implementation, the time-domain processing module determines the highest truncation order N by a step-by-step verification method, and gradually increases the order for each coordinate component starting from the first harmonic, stopping when the residual decreases gradually.

[0103] In one possible implementation, the spatial domain processing module employs a fully normalized spherical harmonic function expansion model: The highest truncation order N of the spherical harmonics and the highest truncation order N of the time domain are independent parameters; the zenith distance is obtained by converting the satellite elevation angle according to z=90°-E, and the azimuth angle is obtained directly from the satellite azimuth angle.

[0104] In one possible implementation, the spatial domain processing module also constructs a zenith virtual constraint equation and assigns it a very high weight to establish the zenith direction as the absolute zero reference for multiple paths; at the same time, it introduces Tikhonov regularization, constructs a diagonal regularization matrix and sets regularization parameters, and solves the spherical harmonic coefficient vector according to the joint method equation.

[0105] In one possible implementation, during the process of solving the spherical harmonic coefficient vector, gross error detection is performed based on the apologetic unit weight mean square error. After reducing the weight of the gross error observations and updating the observation weight matrix, the joint normal equation is solved again, and the iteration continues until convergence.

[0106] In one possible implementation, the acquisition module further includes data preprocessing and quality assessment functions, which perform consistency checks on the input data, filter low-elevation-angle satellite observations according to the set elevation angle cutoff angle, eliminate gross errors and abnormal observations, and check the rationality of the station coordinates.

[0107] In one possible implementation, the system further includes a result output module for outputting a time-domain compensation file, a spatial-domain compensation file, a fusion compensation configuration file, and outputting the fitted residual sequence and standard deviation as quality control indicators. When the fused residual exceeds a preset threshold, an abnormal alarm is output.

[0108] It should be noted that the system adopts a modular architecture to implement the GNSS multipath error compensation software system GDAmm, including: a data management module (reading coordinate time series and residual files), a preprocessing and quality control module (elevation angle filtering, gross error detection, and consistency detection), a time domain component pre-stripping module (calling multi-period joint fitting to complete time domain pre-stripping and residual generation), a spherical harmonic function fitting and modeling module (establishing spherical harmonic function observation equations, zenith virtual constraints, and Tikhonov regularization matrix construction), a normal equation solving module (normal equation analysis and solving, matrix inversion and variance-covariance propagation), a dual-mode fusion modeling module (time domain to spatial domain data relay control and fusion compensation value generation), and a result output module (writing compensation model files and quality control statistics). The software supports multiple systems including GPS, GLONASS, Galileo, and BDS, supports satellite system selection and single-satellite rejection, and allows customization of elevation angle cutoff angle, sampling interval, and fitting order. Data formats support RINEX 2.x and 3.x versions and are compatible with the file formats output by mainstream GNSS data processing software. The software supports both command-line mode and graphical interface mode, and can be deployed locally and run offline without transmitting data externally.

[0109] This invention also provides a computer storage medium, which is a memory device in an electronic device used to store programs and data. It is understood that the computer storage medium here can include both the built-in storage medium in the electronic device and extended storage media supported by the electronic device. The computer storage medium provides storage space that stores the operating system of the electronic device. Furthermore, the storage space also stores one or more instructions suitable for loading and execution by the processor 301. These instructions can be one or more computer programs (including program code). It should be noted that the computer storage medium here can be a high-speed random access memory (RAM) or a non-volatile memory, such as at least one disk storage device; optionally, it can also be at least one computer storage medium located remotely from the aforementioned processor.

[0110] Figure 5 This is a schematic diagram of the structure of an electronic device provided in an embodiment of the present invention. The device includes: at least one processor 301, such as a central processing unit, at least one memory 302, and at least one bus 303.

[0111] The memory 302 can store program instructions, and the processor 301 can be used to call the program instructions to execute the GNSS multipath error modeling and compensation method.

[0112] Those skilled in the art will understand that all or part of the steps in the various methods of the above embodiments can be implemented by a program instructing related hardware. The program can be stored in a computer-readable storage medium, including read-only memory, random access memory, programmable read-only memory, erasable programmable read-only memory, programmable read-only memory, electronically erasable rewritable read-only memory, read-only optical disk, solid-state drive or other optical disk storage, disk storage, magnetic tape storage, or any other computer-readable medium capable of carrying or storing data.

[0113] It should be noted that all steps in the embodiments of the present invention are performed under legal and compliant conditions, that is, all steps in the embodiments of the present invention are performed under authorized conditions.

[0114] The above description discloses only preferred embodiments of the present invention and should not be construed as limiting the scope of the present invention. Therefore, equivalent variations made in accordance with the claims of the present invention are still within the scope of the present invention.

Claims

1. A GNSS multipath error modeling and compensation method, characterized in that, include: Obtain the coordinate time series output by the GNSS baseline calculation, as well as the satellite elevation angle and azimuth angle corresponding to each epoch of the coordinate time series; The coordinate time series is subjected to time domain periodic signal pre-stripping processing to obtain time domain multipath components and residual sequences. The residual sequences are the sequences remaining after stripping periodic signals and trend terms from the coordinate time series. The residuals of each epoch in the residual sequence are converted into zenith distance and azimuth according to the satellite elevation angle and azimuth angle of the corresponding epoch, to obtain spatial domain modeling input data composed of residual values, zenith distance and azimuth angle; wherein, the zenith distance is determined according to the satellite elevation angle and the azimuth angle is determined according to the satellite azimuth angle. Using the residual values ​​in the spatial domain modeling input data as the observed values ​​and the zenith distance and azimuth angle as independent variables, spherical harmonic function fitting modeling is performed to obtain the spatial domain multipath components; Based on the obtained spatial domain multipath components, a secondary residual is calculated on the residual sequence. The residual stellar diurnal component in the secondary residual is fed back to the time domain to update the time domain multipath components and the spherical harmonic function fitting model is re-executed. Cross-domain iteration is performed until convergence. The sum of the converged time domain multipath components and spatial domain multipath components is used as the fused multipath compensation value. The fused multipath compensation value is used to perform multipath error compensation on the coordinate time series.

2. The method according to claim 1, characterized in that, The acquisition of the coordinate time series output from the GNSS baseline solution includes data preprocessing and quality assessment: Perform consistency checks on the input data, examining the temporal continuity of the data, sampling intervals, and observation quality indicators; Filter low elevation angle satellite observations based on the set elevation angle cutoff angle; Remove gross errors and outlier observations, and check the rationality of the station coordinates.

3. The method according to claim 1, characterized in that, The step of performing time-domain periodic signal pre-stripping on the coordinate time series to obtain time-domain multipath components and residual sequences includes: Convert the timestamps of each epoch in the coordinate time series to the sidereal time reference time of the corresponding GNSS system, so that the multi-day data are aligned on the sidereal day time axis. For each component of the coordinate time series, the following joint observation equation is established for multi-period joint fitting: The formula is as follows: Where Y(t) is the coordinate time series, and t is the continuous cumulative time. Let k be the sidereal day period of the GNSS system, and k be the harmonic order. For constant terms, The slope of the linear trend. , Here, N represents the k-th harmonic coefficient, and N is the highest truncation order in the time domain. It is random noise; The parameters of the joint observation equation are obtained by least squares solution, and the time-domain multipath component is extracted from the fitting results. The formula is: and the residual sequence The formula is: 。 4. The method according to claim 2, characterized in that, The highest truncation order N in the time domain is determined by a successive test method: The fitting residuals for each coordinate component under different harmonic orders are calculated. Starting from the first harmonic, the order is gradually increased. When the residual decreases gradually after increasing a certain order, the process stops. The order at this point is taken as the highest truncation order N in the time domain.

5. The method according to claim 1, characterized in that, The method involves using the residual values ​​from the spatial domain modeling input data as observed values, and the zenith distance and azimuth angle as independent variables to perform spherical harmonic function fitting modeling, thereby obtaining the spatial domain multipath components. Specifically: The zenith distance is obtained by converting the satellite elevation angle using z = 90° - E, and the azimuth angle is obtained directly from the satellite azimuth angle, where E is the satellite elevation angle. The spherical harmonic function fitting model adopts the following spherical harmonic function expansion model. The formula is: Where z is the zenith distance and α is the azimuth angle. To fully normalize the Legendre accompanying polynomial, , Let N be the spherical harmonic coefficients to be estimated, N be the highest cutoff order of the spherical harmonics, and m be the order within the order, with m ranging from 0 to n. The summation starts from n=1 to N, and zero-order terms n=0 are removed.

6. The method according to claim 4, characterized in that, The spherical harmonic function fitting modeling also includes the introduction of zenith virtual constraints and Tikhonov regularization: Construct the zenith virtual constraint equations. At the zenith direction (i.e., at zenith distance z=0), based on the Legendre polynomial property, retain only the harmonic terms where m=0 to construct the zenith virtual constraint equations. The formula is: The zenith virtual constraint equation is given extremely high weight. To establish the zenith direction as the absolute zero reference for multiple paths; Tikhonov regularization is introduced to construct a regularization matrix R, which is a diagonal matrix with n (n+1) diagonal elements, and the regularization parameter λ is set. The spherical harmonic coefficient vector X is solved using the joint equation method, as follows: Wherein, B1 is the design matrix for fitting the spherical harmonic function, P1 is the observation weight matrix, L1 is the observation vector, and B2 is the design matrix for the zenith virtual constraint. The joint equation is obtained by minimizing the objective function, and the formula is as follows: , where V1=B1·X-L1.

7. The method according to claim 5, characterized in that, The process of solving the spherical harmonic coefficient vector X also includes gross error detection: After each round of solving, the fitting residual is calculated, and the gross error observations are weighted according to the post-hoc unit weight mean square error, updating the observation weight matrix. Then, the combined equations are solved again, and the above process is repeated until the solution converges.

8. The method according to claim 1, characterized in that, The method also includes outputting model files and quality control information: Output a time-domain compensation file, a spatial-domain compensation file, and a fusion compensation configuration file, wherein the time-domain compensation file records the mapping between the sidereal time position and the time-domain multipath component, the spatial-domain compensation file records the spherical harmonic coefficients and truncation order, and the fusion compensation configuration file records the fusion method and quality control parameters; The output fitted residual sequence and standard deviation are used as quality control indicators. When the fused residual exceeds the preset threshold, an abnormal alarm is output.

9. A GNSS multipath error modeling and compensation system, applicable to the GNSS multipath error modeling and compensation method according to any one of claims 1-8, characterized in that, It includes an acquisition module, a time domain processing module, a spatial domain processing module, a fusion module, and a compensation module, specifically: The acquisition module acquires the coordinate time series output by the GNSS baseline calculation, as well as the satellite elevation angle and azimuth angle corresponding to each epoch of the coordinate time series; The time domain processing module performs time domain periodic signal pre-stripping processing on the coordinate time series to obtain time domain multipath components and residual sequences. The residual sequences are the sequences remaining after stripping periodic signals and trend terms from the coordinate time series. The spatial domain processing module converts the residuals of each epoch in the residual sequence into zenith distance and azimuth according to the satellite elevation angle and azimuth of the corresponding epoch, thereby obtaining spatial domain modeling input data composed of residual values, zenith distance and azimuth. The module then uses the residual values ​​in the spatial domain modeling input data as observation values ​​and the zenith distance and azimuth as independent variables to perform spherical harmonic function fitting modeling, thereby obtaining spatial domain multipath components. The fusion module generates a fused multipath compensation value based on the time-domain multipath component and the spatial-domain multipath component. The spatial-domain multipath component is obtained by fitting the residual sequence with a spherical harmonic function, so that the time-domain processing and the spatial-domain processing form a unidirectional data flow fusion relationship through residual concatenation. Before generating the fused multipath compensation value, a secondary residual is calculated on the residual sequence based on the initially obtained spatial-domain multipath component. The residual sidereal periodic component in the secondary residual is fed back to the time domain to update the time-domain multipath component and the spherical harmonic function fitting model is re-executed. Cross-domain iterative refinement is performed until convergence. The sum of the converged time-domain multipath component and the spatial-domain multipath component is used as the fused multipath compensation value. The compensation module uses the fused multipath compensation value to perform multipath error compensation on the coordinate time series.

10. An electronic device, characterized in that, include: The system includes a memory and a processor, which are interconnected. The memory stores computer instructions, and the processor executes the computer instructions to perform a GNSS multipath error modeling and compensation method according to any one of claims 1 to 8.

Images