Real-time dynamic relative positioning methods, systems, and devices based on state domain correction information

Abstract

This invention discloses a real-time dynamic relative positioning method, system, and device based on state domain correction information. The method involves: a server receiving real-time GNSS observation data streams from the entire network of reference stations and real-time precise orbit clock correction streams to perform real-time precise single-point positioning calculations, and broadcasting network-wide floating-point ambiguity information; then performing real-time phase deviation estimation and broadcasting it; further, combining the network-wide reference station observation data streams and real-time precise orbit clock errors to perform real-time precise single-point positioning fixed-solution calculations, and broadcasting network-wide atmospheric delay information; performing atmospheric modeling and broadcasting the model information; and a user terminal receiving the precise orbit clock correction information, phase deviation information, and model information broadcast by the server in real-time, generating virtual observation values ​​at its approximate coordinates, and performing real-time dynamic relative positioning calculations. This invention can significantly reduce the algorithm complexity of current precise positioning servers and is compatible with existing relative positioning user terminals while ensuring one-way communication and privacy protection.

Description

Technical Field

[0001] This invention belongs to the field of GNSS positioning technology, and particularly relates to a real-time dynamic relative positioning method, system and device based on state domain correction information. Background Technology

[0002] With the increasing intelligence and refinement of social production activities, the demand for real-time high-precision positioning and navigation is becoming increasingly urgent. The basic positioning service of the Global Navigation Satellite System (GNSS) has an accuracy at the meter level, while the integration of various augmentation information can improve positioning accuracy to the centimeter or even millimeter level. Point Positioning-Precise Point Positioning (PPP-RTK) and Network Real-Time Relative Positioning (NRTK) are two representative GNSS real-time high-precision positioning technologies. PPP-RTK technology processes ground network observation data at the server end, generating and broadcasting state domain correction (SSR) information in real time, including precise satellite orbit clock errors, pseudorange bias, phase bias (UPD), and atmospheric delay models, to enhance absolute user positioning. NRTK technology generates observation domain correction (OSR) information by real-time calculation of spatial correlation errors of the reference station network and broadcasts it to relative positioning users, allowing users to achieve RTK calculations without setting up a reference station.

[0003] In theory, PPP-RTK and NRTK technologies are equivalent. However, in practical applications, NRTK, benefiting from its earlier theoretical development, boasts a mature user algorithm model, engineering, and commercial applications. However, its data processing is more complex, requiring bidirectional communication between the user and server. This model heavily relies on the communication environment, failing to protect user privacy and imposing high bandwidth requirements. Compared to NRTK, PPP-RTK technology started later, and its server-side and user-side algorithms are less mature. However, it overcomes NRTK's shortcomings. For example, PPP-RTK has simpler data processing, uses a broadcast-style unidirectional communication model, theoretically has no upper limit on the number of users it can serve, and is more effective at protecting user privacy than NRTK. Furthermore, PPP-RTK uses smaller SSR (Search Server Message) data, thus requiring less bandwidth than NRTK.

[0004] When providing real-time, high-precision positioning services to a massive number of users, service providers want to retain the advantages of PPP-RTK server-side algorithms—simple algorithms, one-way communication, massive concurrency, and privacy protection—while also striving for compatibility with existing RTK user-side algorithm models. However, the reality is that the augmentation information from NRTK and PPP-RTK technologies differs fundamentally from the user algorithm model; SSR augmentation information cannot be directly applied to RTK calculations. To fully integrate the advantages of PPP-RTK and NRTK technologies, this invention proposes a scheme for generating server-side SSR products and performing RTK positioning based on these SSR products, providing users with a novel real-time, high-precision positioning service. Summary of the Invention

[0005] To address the significant differences between the augmented information and user algorithm models of NRTK and PPP-RTK technologies, which prevent SSR information from being directly applied to RTK calculations, this invention provides a real-time dynamic relative positioning method, system, and device based on state domain correction information. After receiving SSR augmented information broadcast by the server, the user converts the SSR into observations from a virtual reference station at the terminal and differentiates them with the actual GNSS observations received by the user, thereby enabling RTK calculations.

[0006] To achieve the above technical objectives, the present invention adopts the following technical solution: Firstly, a real-time dynamic relative positioning method based on state domain correction information is provided, including the following steps: S1: The server receives the GNSS observation data stream from the entire network of reference stations and the real-time precise orbit clock error correction stream in real time, performs concurrent real-time precise single-point positioning calculation for the entire network, and broadcasts the floating-point ambiguity information for the entire network. S2: The server performs real-time phase deviation estimation based on the floating-point ambiguity information of the entire network and broadcasts the phase deviation information. S3: The server performs real-time concurrent precise single-point positioning fixed solution calculation for the entire network based on the GNSS observation data stream of the entire network reference stations, the real-time precise orbit clock error correction stream, and the phase deviation, and broadcasts the atmospheric delay information of the entire network. S4: The server performs atmospheric modeling in real time based on the atmospheric delay information of the entire network and broadcasts the atmospheric delay model information of the entire network. S5: The user terminal receives the precise orbital clock error correction flow, phase deviation information and whole-network atmospheric delay model information broadcast by the server in real time, generates virtual observation values ​​at its own approximate coordinates, and performs real-time dynamic relative positioning calculations.

[0007] Furthermore, the real-time precision orbit clock error correction stream received in step S1 includes the satellite precision orbit, precision clock error, and pseudorange deviation correction.

[0008] Furthermore, in step S1, precise point positioning (PPP) is performed using satellite observations from multiple frequencies and systems; the PPP process is as follows: Precise orbit and clock bias products are obtained by combining broadcast ephemeris data with satellite precise orbit, precise clock bias, and pseudorange bias corrections. Then, parameterized estimates are performed on receiver 3D coordinates, clock bias, ionospheric delay, tropospheric wet delay (ZWD), and ambiguity. ; The parameters in the formula can be expressed as: ; In the formula, and They represent GNSS systems and satellites, respectively. The value can be G, E, or C, where G, E, and C represent GPS, Galileo, and BDS-3, respectively. and These represent the "observed value minus the calculated value" for the pseudorange and phase at the i-th frequency, respectively. For satellite The unit vector of the line of sight direction. For the receiver's three-dimensional coordinates, For receiver clock bias estimation, For the physical clock bias at the receiver end, and for The ionosphere-free combination coefficient of the system, where , , for The frequency of the i-th carrier wave in the system, for System-to-system deviation between the system and GPS and They are respectively System satellite Original and estimated values ​​of first-frequency tilted ionospheric delay. for The ionospheric delay mapping factor at the i-th frequency of the system, , and They are respectively The pseudorange hardware delay at the i-th frequency receiver and satellite ends of the system. and These represent the carrier phase hardware delays at the receiver and satellite ends for the i-th frequency, respectively. For tropospheric wet delay, For satellite Wet delay projection factor, for Estimation of the inter-frequency deviation parameter at the k-th frequency of the system. , For satellite Phase ambiguity estimation at frequency i. For integer ambiguity, Let be the carrier wavelength of the i-th frequency. During the calculation, Kalman filtering is used to estimate the above state parameters to obtain the floating-point ambiguity information of each reference station, and then broadcast it through the NRTIP protocol.

[0009] Furthermore, the method for real-time phase deviation estimation in step S2 is as follows: The method for real-time phase deviation estimation in step S2 is as follows: S21: Perform a linear combination of the original floating-point ambiguities of each reference station; for example, for satellite s and receiver r, the combined floating-point ambiguities are expressed as: ; In the formula, the subscript For different types of linear combinations, For linear combinations Ambiguity, The floating-point ambiguity at the j-th frequency obtained from multi-frequency precise single-point positioning estimation; For linear combinations The integer part of the ambiguity, The decimal part of the corresponding ambiguity; S22: Obtain the fractional part of the linear combination ambiguity by rounding; since the carrier wavelength of the linear combination is relatively long and less affected by other errors, it can be obtained by rounding. and Suppose that n reference stations observe a total of m satellites, then the following equation can be obtained: ; In the formula, This represents the total number of satellites at the nth reference station. The phase deviation observations for quasi-linear combinations correspond to the fractional part of the ambiguity of the linear combination. and Let m and n represent column vectors with dimensions m and n, respectively, where all elements are 1. and Let m and n represent the identity matrices of dimensions m and n, respectively. Indicates the Kronecker product. and Let n and m represent the receiver and satellite respectively. Phase deviation of quasi-linear combinations; S23: Introduce a reference satellite; to eliminate the rank deficiency in the above equation, the satellite observed by the most reference stations will be used. The phase deviation is fixed to zero, i.e., the following constraint is introduced: ; S24: Calculate the original frequency phase deviation; after determining the reference satellite, estimate the phase deviation of the linear combination of each satellite using the least squares method; then obtain the original phase deviation of each frequency of satellite s using the following formula and broadcast it via the NRTIP protocol: ; In the formula, , , These represent the estimated original phase deviations at frequencies 1, 2, and 3, respectively.

[0010] Furthermore, the method for real-time precise single-point positioning fixed solution calculation in step S3 is as follows: using the phase deviation information broadcast by the server, the ambiguity is fixed during precise single-point positioning calculation by means of inter-satellite difference and phase deviation correction, thereby realizing precise single-point positioning fixed solution calculation. At the same time, fixed solution ionospheric delay and tropospheric wet delay information are generated and broadcast through the NRTIP protocol.

[0011] Furthermore, the method for atmospheric modeling in step S4 is as follows: S41: A single-star oblique ionospheric delay model is established based on a second-order polynomial model. This model can be expressed as: ; In the formula, These are values ​​from the ionospheric delay model. and Satellites The difference between the latitude and longitude of the puncture point and the latitude and longitude of the reference center point. These are the coefficients of the ionospheric delay model; when enough ground stations observe satellite s, the least squares method can be used to fit the model coefficients. ; S42: Delayed estimation in the ionosphere Subtracting ionospheric delay model values Computational modeling residuals ,Right now: ; based on A grid model of ionospheric delay residuals is constructed. After considering the spatial correlation of atmospheric delay, an inverse distance-weighted method is used to interpolate the ionospheric residuals to grid points near the puncture point. This method can be expressed as: ; In the formula, This represents the model residual value of the q-th grid point of satellite s. This represents the distance from the p-th ionospheric puncture point to the q-th grid point, where P represents the number of puncture points surrounding the grid point. S43: The server sequentially establishes polynomial and residual grid models for the tilted ionospheric delay of all satellites within the entire reference station network, and broadcasts the model coefficients and grid values ​​to user terminals within the network. After receiving this information, the user terminal first obtains the model value by combining its own approximate coordinates and model coefficients, and then calculates the model residual based on the surrounding four residual grid points, thus obtaining a high-precision ionospheric delay. This process can be represented as: ; In the formula, , and These represent the model value, grid value, and final reference value of the ionospheric delay of satellite S, respectively. S44: When modeling tropospheric delay, a modified polynomial model is used to construct a tropospheric wet delay background model: ; In the formula, For tropospheric wet delay, Here, h represents the tropospheric model coefficients, and h is the ellipsoidal height of the station. After fitting the model coefficients using the least squares method, the tropospheric delay residual grid model is established using the same method as in step S42. After receiving the tropospheric model coefficients and the corresponding residual grid, the user terminal can calculate the tropospheric wet delay reference value. .

[0012] Furthermore, step S5 generates virtual observations at its own approximate coordinates as follows: S51: The user determines the approximate coordinates of the current epoch and the receiver clock error through single-point positioning, and sets the approximate coordinates as the position of the virtual reference station; S52: Perform tidal correction on the virtual reference station coordinates, then subtract the receiver clock error from the observation time to obtain the signal reception time, and then calculate the initial satellite coordinates assuming the satellite-to-ground distance is 0. Iterate the satellite-to-ground distance based on the initial satellite coordinates until convergence, thus obtaining the signal transmission time and the precise satellite coordinates and clock error. S53: Calculate the satellite elevation and azimuth angles, then calculate the tropospheric dry delay, tropospheric projection coefficient, and antenna phase winding correction terms related to the elevation and azimuth angles based on the corresponding model. Additionally, obtain the ionospheric delay, tropospheric wet delay, phase deviation, and pseudorange deviation correction terms based on the received server broadcast SSR information. Finally, obtain the pseudorange and carrier phase observations of the virtual reference station; these observations are expressed as: ; In the formula, the distance between the satellite and the ground is... Represented as: ; In the formula, Indicates satellite, , Let these represent the virtual pseudorange and carrier phase observations at the j-th frequency, respectively. For the precise clock difference of the satellite at the time of signal transmission, This is a relativistic correction term. For the receiver's approximate clock bias, This is a reference value for ionospheric delay. Let j be the ionospheric delay mapping factor at frequency j. For satellite Wet delay projection factor, This is a reference value for tropospheric wet delay. For the pseudorange hardware delay at the j-th frequency satellite end, This is a phase winding correction term. Let be the unit vector in the line-of-sight direction from the receiver to the satellite. The position of the satellite's center of mass at the moment of signal transmission. The coordinates of the receiver antenna reference point at the moment of signal reception. These are the PCO and PCV correction values ​​for the satellite antenna. These are the receiver antenna PCO and PCV correction values. These are tidal load and Earth's rotation error, respectively. Estimate the original phase deviation at frequency j. Among these parameters, and Using the results of SPP, , , , Various PCO and PCV correction values ​​are calculated using a priori error models. , , , , The pseudorange deviation is obtained based on the information broadcast by the server.

[0013] Furthermore, the method for real-time dynamic relative positioning calculation in step S5 is as follows: the user performs real-time dynamic relative positioning calculation based on the observations of the recovered virtual reference station and the actual received GNSS observations. This process includes: S54: In satellite By performing inter-station differencing between actual and virtual observations, we can obtain single-difference inter-station observations: ; in, Indicates on the real site With virtual sites The results of the single-difference analysis between stations were obtained. It is an inter-station single difference operator. This represents the single-difference pseudorange observation value between stations at frequency i. Let be the wavelength of the carrier wave at frequency i. This represents the inter-station single-difference carrier phase observation value at the i-th frequency. For single-difference satellite-to-ground distance, For single-difference receiver clock bias, and These are the single-difference ionospheric delay and the tropospheric wet delay, respectively. Let i be the ionospheric delay mapping factor at frequency i. and These represent the pseudorange and phase hardware delay at the single-difference receiver at the i-th frequency, respectively. For single-difference ambiguity, and These are single-difference pseudorange and phase observation noise, respectively; based on this, the satellite... , Inter-station single-difference observations are quantified into inter-satellite differences to obtain double-difference observations: ; in, It is a double difference operator. and These are the double-difference pseudorange and double-difference carrier phase observations, respectively. It is a double-difference satellite distance. For double-difference ambiguity, and To address the double-difference pseudorange and phase observation noise, Kalman filtering is used to solve the above equation, and the double-difference ambiguity is fixed, thus enabling real-time dynamic relative positioning.

[0014] Secondly, a real-time dynamic relative positioning system based on state domain correction information is provided, including a server and a user terminal; The server includes a real-time PPP floating-point solution module, a real-time UPD estimation module, a real-time fixed solution module, and a real-time atmospheric modeling module; The real-time PPP floating-point solution module is used to receive the GNSS observation data stream of the entire network reference station and the real-time precise orbit clock error correction stream in real time, perform concurrent real-time precise single-point positioning solution for the entire network, and broadcast the floating-point ambiguity information of the entire network. The real-time UPD estimation module is used to receive the floating-point ambiguity information of the entire network in real time, perform real-time phase deviation estimation, and broadcast the phase deviation information. The real-time fixed solution module is used to receive GNSS observation data streams from the entire network of reference stations, real-time precise orbit clock error correction streams, and phase deviations in real time, perform network-wide concurrent real-time precise single-point positioning fixed solution calculations, and broadcast network-wide atmospheric delay information. The real-time atmospheric modeling module is used to receive atmospheric delay information from the entire network in real time, perform atmospheric modeling, and broadcast the atmospheric delay model information from the entire network. The user end includes a real-time virtual observation recovery and RTK module, which is used to receive the precise orbital clock error correction flow, phase deviation information and whole-network atmospheric delay model information broadcast by the server in real time, generate virtual observations at its own approximate coordinates, and perform real-time dynamic relative positioning calculations.

[0015] Thirdly, a real-time dynamic relative positioning device based on state domain correction information is provided, including a memory and a processor. The memory stores a computer program, and when the computer program is executed by the processor, the processor enables the processor to implement the real-time dynamic relative positioning method based on state domain correction information as described above.

[0016] The beneficial effects of this invention are: the method proposed in this invention fully integrates the technical features of PPP-RTK and NRTK, and can meet the needs of real-time precise positioning services such as one-way communication, protection of user privacy, massive concurrency, and support for user relative positioning calculation without changing the current PPP-RTK server-side algorithm model and RTK user-side algorithm model and significantly reducing the cost of technology updates. Attached Figure Description

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

[0018] Figure 1 This is a flowchart of the real-time dynamic relative positioning method based on state domain correction information provided in an embodiment of the present invention; Figure 2 This is a schematic diagram of a real-time dynamic relative positioning system based on state domain correction information provided in an embodiment of the present invention; Figure 3 This is a schematic diagram of the double-difference residual sequence of virtual and real observations of five-frequency pseudorange and phase at WRRN station BDS-3 on April 5, 2025, provided by an embodiment of the present invention; wherein, (a)-(e) represent the double-difference residual distribution of five-frequency pseudorange of BDS-3 B1I / B3I / B1C / B2a / B2b, respectively, and (f)-(j) represent the double-difference residual distribution of five-frequency carrier phase, respectively; Figure 4 The present invention provides a 24-hour continuous GPS / Galileo / BDS-3 three-system single-frequency RTK positioning error sequence for the RGLN station based on virtual observations, wherein (a), (b), and (c) represent the positioning errors in the east, north, and elevation directions, respectively. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be described in detail below. Obviously, the described embodiments are merely some embodiments of this invention, and not all embodiments. Based on the embodiments of this invention, all other implementation methods obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0020] This invention provides a real-time dynamic relative positioning method based on state domain correction information, referencing... Figure 1 As shown, it includes the following steps: S1, the server receives real-time GNSS observation data streams from the entire network of reference stations and real-time precise orbit clock correction streams, performs concurrent real-time precise point positioning (PPP) calculations for the entire network, and broadcasts the entire network's floating-point ambiguity information. This step uses satellite observations from multiple frequencies and systems for precise point positioning (PPP). The PPP process is as follows: Precise orbit and clock bias products are obtained by combining broadcast ephemeris and precise orbit, clock bias, and pseudorange bias corrections. Then, parameterized estimates are performed on the receiver's three-dimensional coordinates, clock bias, ionospheric delay, tropospheric wet delay (ZWD), and ambiguity. ; The parameters in the formula can be expressed as: ; In the formula, and They represent GNSS systems and satellites, respectively. The value can be G, E or C, where G, E and C represent GPS (GPS system), Galileo (Galileo system) and BDS-3 (BeiDou system) respectively. and These represent the "observed value minus the calculated value" for the pseudorange and phase at the i-th frequency, respectively. For satellite The unit vector of the line of sight direction. For the receiver's three-dimensional coordinates, For receiver clock bias estimation, For the physical clock bias at the receiver end, and for The ionosphere-free combination coefficient of the system, where , , for The frequency of the i-th carrier wave in the system, for System-to-system deviation between the system and GPS and They are respectively System satellite Original and estimated values ​​of first-frequency tilted ionospheric delay. for The ionospheric delay mapping factor at the i-th frequency of the system, , and They are respectively The pseudorange hardware delay at the i-th frequency receiver and satellite ends of the system. and These represent the carrier phase hardware delays at the receiver and satellite ends for the i-th frequency, respectively. For tropospheric wet delay, For satellite Wet delay projection factor, for Estimation of the inter-frequency deviation parameter at the k-th frequency of the system. , For satellite Phase ambiguity estimation at frequency i. For integer ambiguity, Let be the carrier wavelength of the i-th frequency. During the calculation, Kalman filtering is used to estimate the above state parameters to obtain the floating-point ambiguity information of each reference station, and then broadcast it through the NRTIP protocol.

[0021] S2, the server receives the floating-point ambiguity information of the entire network in real time, performs real-time phase deviation estimation, and broadcasts the phase deviation information. The method for real-time phase deviation estimation is as follows: S21: Perform a linear combination of the original floating-point ambiguities of each reference station; for example, for satellite s and receiver r, the combined floating-point ambiguity can be expressed as: ; In the formula, the subscript For different types of linear combinations, For linear combinations Ambiguity, The floating-point ambiguity at the j-th frequency obtained from multi-frequency precise single-point positioning estimation; For linear combinations The integer part of the ambiguity, This represents the decimal part of the corresponding ambiguity.

[0022] S22: Obtain the fractional part of the linear combination ambiguity by rounding. Since the carrier wavelength of the linear combination is relatively long and less affected by other errors, it can be obtained by rounding. and Suppose that n reference stations observe a total of m satellites, then the following equation can be obtained: ; In the formula, This represents the total number of satellites at the nth reference station. The phase deviation observations for quasi-linear combinations correspond to the fractional part of the ambiguity of the linear combination. and Let m and n represent column vectors with dimensions m and n, respectively, where all elements are 1. and Let m and n represent the identity matrices of dimensions m and n, respectively. Indicates the Kronecker product. and Let n and m represent the receiver and satellite respectively. Phase deviation of a linear combination.

[0023] S23: Introduce a reference satellite; to eliminate the rank deficiency in the above equation, the satellite observed by the most reference stations will be used. The phase deviation is fixed to zero, i.e., the following constraint is introduced: ; S24: Calculate the original frequency phase deviation. After determining the reference satellite, the phase deviation of the linear combination of each satellite can be estimated using the least squares method. Then, the original phase deviation of each frequency of satellite 's' can be obtained using the following formula and broadcast via the NRTIP protocol: ; In the formula, , , These represent the estimated original phase deviations at frequencies 1, 2, and 3, respectively.

[0024] Step 3: The server receives real-time GNSS observation data streams, real-time precise orbit clock errors, and phase deviations from the entire network of reference stations. It then performs concurrent real-time precise single-point positioning fixed-solution calculations for the entire network and broadcasts the network-wide atmospheric delay information. The method for performing real-time PPP fixed-solution calculations is as follows: using the phase deviation information broadcast by the server, ambiguities are fixed during PPP calculations through inter-satellite differences and phase deviation corrections, thus achieving PPP fixed-solution calculations. Simultaneously, fixed-solution ionospheric delay and tropospheric wet delay information are generated and broadcast via the NRTIP protocol.

[0025] Step 4: The server receives the fixed-solution atmospheric delay information from the entire network in real time, performs atmospheric modeling, and broadcasts the atmospheric delay model information for the entire network. The atmospheric modeling method is as follows: S41: A single-star oblique ionospheric delay model is established based on a second-order polynomial model. This model can be expressed as: ; In the formula, These are values ​​from the ionospheric delay model. and Satellites The difference between the latitude and longitude of the puncture point and the latitude and longitude of the reference center point. These are the coefficients of the polynomial model. When enough ground stations observe satellite s, the model coefficients can be fitted using the least squares method. .

[0026] S42: Delayed estimation in the ionosphere Subtracting ionospheric delay model values Computational modeling residuals ,Right now: ; based on A grid model of ionospheric delay residuals is further constructed. Taking into account the spatial correlation of atmospheric delay, an inverse distance-weighted method is used to interpolate the ionospheric residuals to grid points near the puncture point. This method can be expressed as: ; In the formula, This represents the model residual value of the q-th grid point of satellite s. This represents the distance from the p-th ionospheric puncture point to the q-th grid point, where P represents the number of puncture points around the grid point.

[0027] S43: The server will sequentially establish polynomial models and residual grids for the tilted ionospheric delays of all satellites within the entire reference station network, and broadcast the model coefficients and grid values ​​to users within the network. After receiving this information, users first obtain the model values ​​by combining their own approximate coordinates and model coefficients, and then calculate the model residuals based on the surrounding four residual grid points, thus obtaining a high-precision ionospheric delay. This process can be represented as:

[0028] ;

[0029] In the formula, , and These represent the model value, grid value, and final reference value of the ionospheric delay of satellite S, respectively.

[0030] S44: When modeling tropospheric delay, a modified polynomial model is used to construct a tropospheric wet delay background model: ; In the formula, For tropospheric wet delay, Here, h represents the tropospheric model coefficients, and h is the ellipsoidal height of the station. After fitting the model coefficients using the least squares method, the tropospheric delay residual grid model can be established using the same method as in step S42. After receiving the tropospheric model coefficients and residual grid, the user can obtain the tropospheric wet delay (ZWD) reference value using a method similar to step S43. .

[0031] Step 5: The user terminal receives the precise orbital clock correction flow, phase deviation information, and overall atmospheric delay model information broadcast by the server in real time. It then generates virtual observations at its own approximate coordinates and performs real-time dynamic relative positioning calculations. The method for generating virtual observations at its own approximate coordinates in this step is as follows:

[0032] S51: The user determines the approximate coordinates of the current epoch and the receiver clock error through single-point positioning, and sets the approximate coordinates as the position of the virtual reference station.

[0033] S52: Perform tidal correction on the virtual reference station coordinates, then subtract the receiver clock error from the observation time to obtain the signal reception time, and then calculate the initial satellite coordinates assuming the satellite-to-ground distance is 0. Iterate the satellite-to-ground distance based on the initial coordinates until convergence, thus obtaining the signal transmission time and the precise satellite coordinates and clock error.

[0034] S53: Calculate the satellite elevation angle and azimuth angle, then calculate corrections related to elevation and azimuth angles such as tropospheric interference delay, tropospheric projection coefficient, and antenna phase winding based on the corresponding model. Additionally, obtain corrections for ionospheric delay, ZWD, phase deviation, and pseudorange deviation based on the received server broadcast SSR information. Finally, obtain the pseudorange and carrier phase observations of the virtual reference station. This observation can be expressed as: ; In the formula, the distance between the satellite and the ground is... This can be further expressed as: ; In the formula, Indicates satellite, , Let these represent the virtual pseudorange and carrier phase observations at the j-th frequency, respectively. For the precise clock difference of the satellite at the time of signal transmission, This is a relativistic correction term. For the receiver's approximate clock bias, This is a reference value for ionospheric delay. Let j be the ionospheric delay mapping factor at frequency j. For satellite Wet delay projection factor, This is a reference value for tropospheric wet delay. For the pseudorange hardware delay at the j-th frequency satellite end, This is a phase winding correction term. Let be the unit vector in the line-of-sight direction from the receiver to the satellite. The position of the satellite's center of mass at the moment of signal transmission. The coordinates of the receiver antenna reference point at the moment of signal reception. These are the PCO and PCV correction values ​​for the satellite antenna. These are the receiver antenna PCO and PCV correction values. These are tidal load and Earth's rotation error, respectively. Estimate the original phase deviation at frequency j. Among these parameters, and Using the results of SPP, , , , Various PCO and PCV correction values ​​are calculated using a priori error models. , , , , The pseudorange deviation is obtained based on the information broadcast by the server.

[0035] The method for real-time dynamic relative positioning calculation in this step is as follows: the user performs real-time dynamic relative positioning calculation based on the observation values ​​of the recovered virtual reference station and the actual received GNSS observation values. This process can be represented as: S54: In satellite By performing inter-station differencing between actual and virtual observations, we can obtain single-difference inter-station observations: ; in, Indicates on the real site With virtual sites The results of the single-difference analysis between stations were obtained. It is an inter-station single difference operator. This represents the single-difference pseudorange observation value between stations at frequency i. Let be the wavelength of the carrier wave at frequency i. This represents the inter-station single-difference carrier phase observation value at the i-th frequency. For single-difference satellite-to-ground distance, For single-difference receiver clock bias, and These are the single-difference ionospheric delay and the tropospheric wet delay, respectively. Let i be the ionospheric delay mapping factor at frequency i. and These represent the pseudorange and phase hardware delay at the single-difference receiver at the i-th frequency, respectively. For single-difference ambiguity, and These are single-difference pseudorange and phase observation noise, respectively; based on this, the satellite... , Inter-station single-difference observations are quantified into inter-satellite differences to obtain double-difference observations: ; in, It is a double difference operator. and These are the double-difference pseudorange and double-difference carrier phase observations, respectively. It is a double-difference satellite distance. For double-difference ambiguity, and This addresses the noise from double-difference pseudorange and phase observations. By performing inter-satellite single-difference on top of inter-station single-difference, receiver clock errors and hardware delays are eliminated, and ionospheric and tropospheric delay errors are further reduced, resulting in highly accurate double-difference observations. Real-time dynamic relative positioning can be achieved by using Kalman filtering to solve the above equation and fixing the double-difference ambiguity.

[0036] This invention also provides a real-time dynamic relative positioning system based on state domain correction information, including a server and a user terminal; The server includes a real-time PPP floating-point solution module, a real-time UPD estimation module, a real-time fixed solution module, and a real-time atmospheric modeling module; The real-time PPP floating-point solution module is used to receive the GNSS observation data stream of the entire network reference station and the real-time precise orbit clock error correction stream in real time, perform concurrent real-time precise single-point positioning solution for the entire network, and broadcast the floating-point ambiguity information of the entire network. The real-time UPD estimation module is used to receive the floating-point ambiguity information of the entire network in real time, perform real-time phase deviation estimation, and broadcast the phase deviation information. The real-time fixed solution module is used to receive GNSS observation data streams from the entire network of reference stations, real-time precise orbit clock error correction streams, and phase deviations in real time, perform network-wide concurrent real-time precise single-point positioning fixed solution calculations, and broadcast network-wide atmospheric delay information. The real-time atmospheric modeling module is used to receive atmospheric delay information from the entire network in real time, perform atmospheric modeling, and broadcast the atmospheric delay model information from the entire network. The user end includes a real-time virtual observation recovery and RTK module, which is used to receive the precise orbital clock error correction flow, phase deviation information and whole-network atmospheric delay model information broadcast by the server in real time, generate virtual observations at its own approximate coordinates, and perform real-time dynamic relative positioning calculations.

[0037] Furthermore, this embodiment of the invention also provides a real-time dynamic relative positioning device based on state domain correction information, including a memory and a processor. The memory stores a computer program, and when the computer program is executed by the processor, the processor implements the real-time dynamic relative positioning method based on state domain correction information as described above.

[0038] Example 1

[0039] This example illustrates the accuracy of virtual observations.

[0040] Figure 3 This figure shows the double-difference residual sequence of virtual and real pseudorange and phase observations from the WRRN station BDS-3 five-frequency system on April 5, 2025. Theoretically, both results should be white noise. As can be seen from the figure, the double-difference results of the pseudorange observations vary within a range of about 2m, exhibiting white noise characteristics overall, while the carrier phase observations remain within ±0.2 cycles for most of the time, which meets the requirement of fixed RTK double-difference ambiguity.

[0041] Example 2

[0042] This example illustrates the positioning accuracy of RTK using virtual observations.

[0043] Figure 4 The figure shows the GPS / Galileo / BDS-3 three-system single-frequency RTK positioning error sequence of the RGLN station over 24 hours based on virtual observations, with initialization performed every 4 hours during positioning. As can be seen from the figure, the results for each time period exhibit instantaneous convergence, and after convergence, the results achieve positioning performance with horizontal accuracy better than 5 cm and vertical accuracy better than 10 cm.

[0044] Examples 1 and 2 together demonstrate that the recovered virtual observations can support user RTK calculations and that the positioning results can achieve instantaneous convergence.

[0045] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0046] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0047] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0048] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0049] Although embodiments of the present invention have been shown and described above, it is understood that the above embodiments are exemplary and should not be construed as limiting the present invention. Those skilled in the art can make changes, modifications, substitutions and variations to the above embodiments within the scope of the present invention.

Claims

1. A real-time dynamic relative positioning method based on state domain correction information, characterized in that, Includes the following steps: S1: The server receives the GNSS observation data stream from the entire network of reference stations and the real-time precise orbit clock error correction stream in real time, performs concurrent real-time precise single-point positioning calculation for the entire network, and broadcasts the floating-point ambiguity information for the entire network. S2: The server performs real-time phase deviation estimation based on the floating-point ambiguity information of the entire network and broadcasts the phase deviation information. S3: The server performs real-time concurrent precise single-point positioning fixed solution calculation for the entire network based on the GNSS observation data stream of the entire network reference stations, the real-time precise orbit clock error correction stream, and the phase deviation, and broadcasts the atmospheric delay information of the entire network. S4: The server performs atmospheric modeling in real time based on the atmospheric delay information of the entire network and broadcasts the atmospheric delay model information of the entire network. S5: The user terminal receives the precise orbital clock error correction flow, phase deviation information and whole-network atmospheric delay model information broadcast by the server in real time, generates virtual observation values ​​at its own approximate coordinates, and performs real-time dynamic relative positioning calculations.

2. The real-time dynamic relative positioning method based on state domain correction information according to claim 1, characterized in that, The real-time precision orbit clock error correction stream received in step S1 includes the satellite precision orbit, precision clock error, and pseudorange deviation correction.

3. The real-time dynamic relative positioning method based on state domain correction information according to claim 1, characterized in that, Step S1 uses satellite observations from multiple frequencies and systems to perform precise point positioning. The precise point positioning process includes: By combining broadcast ephemeris data with satellite precise orbit, precise clock bias, and pseudorange bias corrections, precise orbit and clock bias products are obtained, and parameterized estimates are performed on receiver 3D coordinates, clock bias, ionospheric delay, tropospheric wet delay, and ambiguity. ； The parameters in the formula can be expressed as: ； In the formula, and They represent GNSS systems and satellites, respectively. The value can be G, E, or C, where G, E, and C represent GPS, Galileo, and BDS-3, respectively. and These represent the "observed value minus the calculated value" for the pseudorange and phase at the i-th frequency, respectively. For satellite unit vector of line of sight direction For the receiver's three-dimensional coordinates, For receiver clock bias estimation, For the physical clock bias at the receiver end, and for The system's non-ionospheric combination coefficient, where , , for The frequency of the i-th carrier wave in the system, for System-to-system deviation between the system and GPS and They are respectively System satellite Original and estimated values ​​of first-frequency tilted ionospheric delay. for The ionospheric delay mapping factor at the i-th frequency of the system, , and They are respectively The pseudorange hardware delay at the i-th frequency receiver and satellite ends of the system. and These represent the carrier phase hardware delays at the receiver and satellite ends for the i-th frequency, respectively. For tropospheric wet delay, For satellite Wet delay projection factor, for Estimation of the inter-frequency deviation parameter at the k-th frequency of the system. , For satellite Phase ambiguity estimation at frequency i. For integer ambiguity, Let be the carrier wavelength of the i-th frequency. During the calculation, Kalman filtering is used to estimate the above state parameters to obtain the floating-point ambiguity information of each reference station, and then broadcast it through the NRTIP protocol.

4. The real-time dynamic relative positioning method based on state domain correction information according to claim 1, characterized in that, The method for real-time phase deviation estimation in step S2 is as follows: S21: Perform linear combination of the original floating-point ambiguities of each reference station; for satellite s and receiver r, the combined floating-point ambiguity is expressed as: ； In the formula, the subscript For different types of linear combinations, For linear combinations Ambiguity, The floating-point ambiguity at the j-th frequency obtained from multi-frequency precise single-point positioning estimation; For linear combinations The integer part of the ambiguity, The decimal part of the corresponding ambiguity; S22: Obtain the fractional part of the linear combination ambiguity by rounding; assuming n base stations observe a total of m satellites, the following equation can be obtained: ； In the formula, This represents the total number of satellites at the nth reference station. The phase deviation observations for quasi-linear combinations correspond to the fractional part of the ambiguity of the linear combination. and Let m and n represent column vectors with dimensions m and n, respectively, where all elements are 1. and Let m and n represent the identity matrices of dimensions m and n, respectively. Indicates the Kronecker product. and Let n and m represent the receiver and satellite respectively. Phase deviation of quasi-linear combinations; S23: Introduce a reference satellite; the satellite observed by the most reference stations will be included. The phase deviation is fixed to zero, i.e., the following constraint is introduced: ； S24: Calculate the original frequency phase deviation; after determining the reference satellite, estimate the phase deviation of the linear combination of each satellite using the least squares method; then obtain the original phase deviation of each frequency of satellite s using the following formula and broadcast it via the NRTIP protocol: ； In the formula, , , These represent the estimated original phase deviations at frequencies 1, 2, and 3, respectively.

5. The real-time dynamic relative positioning method based on state domain correction information according to claim 1, characterized in that, The method for real-time precise single-point positioning fixed solution calculation in step S3 is as follows: using the phase deviation information broadcast by the server, the ambiguity is fixed by inter-satellite difference and phase deviation correction during precise single-point positioning calculation, thereby realizing precise single-point positioning fixed solution calculation. At the same time, fixed solution ionospheric delay and tropospheric wet delay information are generated and broadcast through the NRTIP protocol.

6. The real-time dynamic relative positioning method based on state domain correction information according to claim 1, characterized in that, The method for atmospheric modeling in step S4 is as follows: S41: A single-star oblique ionospheric delay model is established based on a second-order polynomial model. This model can be expressed as: ； In the formula, These are values ​​from the ionospheric delay model. and Satellites The difference between the latitude and longitude of the puncture point and the latitude and longitude of the reference center point. These are the coefficients of the ionospheric delay model; S42: Delayed estimation in the ionosphere Subtracting ionospheric delay model values Computational modeling residuals ,Right now: ； based on A grid model of ionospheric delay residuals is constructed. Taking into account the spatial correlation of atmospheric delay, an inverse distance-weighted method is used to interpolate the ionospheric residuals to grid points near the puncture point. This method is expressed as follows: ； In the formula, This represents the model residual value of the q-th grid point of satellite s. This represents the distance from the p-th ionospheric puncture point to the q-th grid point, where P represents the number of puncture points surrounding the grid point. S43: The server sequentially establishes polynomial and residual grid models for the tilted ionospheric delay of all satellites within the entire reference station network, and broadcasts the model coefficients and grid values ​​to user terminals within the network. After receiving this information, the user terminal first obtains the model value by combining its own approximate coordinates and model coefficients, and then calculates the model residual based on the surrounding four residual grid points, thus obtaining a high-precision ionospheric delay. This process is represented as: ； In the formula, , and These represent the model value, grid value, and final reference value of the ionospheric delay of satellite S, respectively. S44: When modeling tropospheric delay, a modified polynomial model is used to construct a tropospheric wet delay background model: ； In the formula, For tropospheric wet delay, Here, h represents the tropospheric model coefficients, and h is the ellipsoidal height of the station. After fitting the model coefficients using the least squares method, the tropospheric delay residual grid model is established using the same method as in step S42. After receiving the tropospheric model coefficients and the corresponding residual grid, the user terminal can calculate the tropospheric wet delay reference value. .

7. The real-time dynamic relative positioning method based on state domain correction information according to claim 1, characterized in that, Step S5 involves generating virtual observations at its approximate coordinates as follows: S51: The user determines the approximate coordinates of the current epoch and the receiver clock error through single-point positioning, and sets the approximate coordinates as the position of the virtual reference station; S52: Perform tidal correction on the virtual reference station coordinates, then subtract the receiver clock error from the observation time to obtain the signal reception time, and then calculate the initial satellite coordinates assuming the satellite-to-ground distance is 0. Iterate the satellite-to-ground distance based on the initial satellite coordinates until convergence, thus obtaining the signal transmission time and the precise satellite coordinates and clock error. S53: Calculate the satellite elevation and azimuth angles, then calculate the tropospheric dry delay, tropospheric projection coefficient, and antenna phase winding correction terms related to the elevation and azimuth angles based on the corresponding model. Additionally, obtain the ionospheric delay, tropospheric wet delay, phase deviation, and pseudorange deviation correction terms based on the received server broadcast state domain correction information. Finally, obtain the pseudorange and carrier phase observations of the virtual reference station; these observations are expressed as: ； In the formula, the distance between the satellite and the ground is... Represented as: ； In the formula, Indicates satellite, , Let these represent the virtual pseudorange and carrier phase observations at the j-th frequency, respectively. For the precise clock difference of the satellite at the time of signal transmission, This is a relativistic correction term. For the receiver's approximate clock bias, This is a reference value for ionospheric delay. Let j be the ionospheric delay mapping factor at frequency j. For satellite Wet delay projection factor, This is a reference value for tropospheric wet delay. For the pseudorange hardware delay at the j-th frequency satellite end, This is a phase winding correction term. Let be the unit vector in the line-of-sight direction from the receiver to the satellite. The position of the satellite's center of mass at the moment of signal transmission. The coordinates of the receiver antenna reference point at the moment of signal reception. These are the PCO and PCV correction values ​​for the satellite antenna. These are the receiver antenna PCO and PCV correction values. These are tidal load and Earth's rotation error, respectively. This is the estimated value of the original phase deviation at frequency j.

8. The real-time dynamic relative positioning method based on state domain correction information according to claim 1, characterized in that, Step S5 involves performing real-time dynamic relative positioning calculations: the user performs real-time dynamic relative positioning calculations based on the recovered virtual reference station observations and the actual received GNSS observations. This process includes: S54: In satellite By performing inter-station differencing between actual and virtual observations, we can obtain single-difference inter-station observations: ； in, Indicates on the real site With virtual sites The results of the single-difference analysis between stations were obtained. It is an inter-station single difference operator. This represents the single-difference pseudorange observation value between stations at frequency i. Let be the wavelength of the carrier wave at frequency i. This represents the inter-station single-difference carrier phase observation value at the i-th frequency. For single-difference satellite-to-ground distance, For single-difference receiver clock bias, and These are single-difference ionospheric delay and tropospheric wet delay, respectively. Let i be the ionospheric delay mapping factor at frequency i. and These represent the pseudorange and phase hardware delay at the single-difference receiver at the i-th frequency, respectively. For single-difference ambiguity, and These are single-difference pseudorange and phase observation noise, respectively; based on this, the satellite... , Inter-station single-difference observations are quantified into inter-satellite differences to obtain double-difference observations: ； in, It is a double difference operator. and These are the double-difference pseudorange and double-difference carrier phase observations, respectively. The satellite distance is double-difference. For double-difference ambiguity, and To address the double-difference pseudorange and phase observation noise, Kalman filtering is used to solve the above equation, and the double-difference ambiguity is fixed, thus enabling real-time dynamic relative positioning.

9. A real-time dynamic relative positioning system based on state domain correction information, characterized in that, Including server-side and user-side; The server includes a real-time PPP floating-point solution module, a real-time UPD estimation module, a real-time fixed solution module, and a real-time atmospheric modeling module; The real-time PPP floating-point solution module is used to receive the GNSS observation data stream of the entire network reference station and the real-time precise orbit clock error correction stream in real time, perform concurrent real-time precise single-point positioning solution for the entire network, and broadcast the floating-point ambiguity information of the entire network. The real-time UPD estimation module is used to receive the floating-point ambiguity information of the entire network in real time, perform real-time phase deviation estimation, and broadcast the phase deviation information. The real-time fixed solution module is used to receive GNSS observation data streams from the entire network of reference stations, real-time precise orbit clock error correction streams, and phase deviations in real time, perform network-wide concurrent real-time precise single-point positioning fixed solution calculations, and broadcast network-wide atmospheric delay information. The real-time atmospheric modeling module is used to receive atmospheric delay information from the entire network in real time, perform atmospheric modeling, and broadcast the atmospheric delay model information from the entire network. The user end includes a real-time virtual observation recovery and RTK module, which is used to receive the precise orbital clock error correction flow, phase deviation information and whole-network atmospheric delay model information broadcast by the server in real time, generate virtual observations at its own approximate coordinates, and perform real-time dynamic relative positioning calculations.

10. A real-time dynamic relative positioning device based on state domain correction information, comprising a memory and a processor, wherein the memory stores a computer program, characterized in that, When the computer program is executed by the processor, the processor causes the processor to implement the method as described in any one of claims 1 to 8.

Images