Real-time PPP positioning method considering troposphere anisotropy
By introducing the solar elevation angle parameter into real-time PPP, an improved wet mapping function and stochastic model are constructed, solving the problem that the influence of tropospheric anisotropy is not considered in traditional PPP, and achieving high-precision and fast positioning results.
Patent Information
- Application Number
- CN202511971889.2
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2025-12-25
- Publication Date
- 2026-03-20
AI Technical Summary
Existing real-time precise point positioning (PPP) technology fails to effectively account for atmospheric anisotropy caused by solar radiation when dealing with tropospheric delay, resulting in decreased positioning accuracy and prolonged convergence time, which is particularly pronounced in high-latitude regions and under low solar altitude angle conditions.
By introducing the solar elevation angle parameter, an improved wet mapping function and stochastic model that take into account the anisotropy of the troposphere are constructed. Parameter estimation is performed using Kalman filtering to correct the solar radiation influence that is ignored in the traditional model.
It improves the positioning accuracy and convergence speed of real-time PPP, especially in areas and conditions where traditional models are weak, meeting the requirements of high-precision real-time positioning, and has the characteristics of environmental robustness and controllable computational burden.
Smart Images
Figure CN121703849A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of satellite navigation and positioning technology, and in particular to a real-time PPP positioning method that takes into account tropospheric anisotropy. Background Technology
[0002] Real-time precise point positioning technology, as a core development direction for high-precision positioning in Global Navigation Satellite Systems (GNSS), boasts advantages such as independent operation of a single receiver, wide coverage, and positioning accuracy ranging from centimeters to decimeters. It has been widely applied in surveying engineering, autonomous driving, crustal deformation monitoring, meteorological services, and disaster emergency response. This technology relies on real-time GNSS observation data and precise orbital clock error products, achieving high-precision positioning through rigorous error models and parameter estimation.
[0003] In handling tropospheric delay errors, current real-time PPP commonly employs empirical mapping functions (such as GMF, VMF1, etc.) to project the zenith delay onto the satellite direction and uses a stochastic model constructed based on the satellite elevation angle to weight observations (e.g., Chinese authorized patent, publication number CN115051749B, publication date 2022-11-04, disclosing an automatic network topology design method and related equipment for satellite networking). However, these traditional methods only consider the influence of the satellite elevation angle, neglecting the anisotropic characteristics of the troposphere itself, especially the spatial differences in atmospheric thermal and moisture distribution caused by solar radiation. Recent studies have shown that solar radiation causes directional changes in the atmospheric refractive index in the troposphere, thus introducing systematic biases in GNSS signal propagation.
[0004] In summary, the shortcomings of the existing technologies mentioned above are mainly reflected in the following three aspects: 1) The construction of the mapping function only takes the satellite elevation angle as the core parameter, without taking into account the modulation effect of the sun's azimuth and solar radiation intensity on the vertical structure and horizontal distribution of the atmosphere. It cannot accurately describe the tropospheric delayed projection law under anisotropic conditions, resulting in a significant increase in delayed projection error in areas with strong solar radiation. 2) The weight allocation of the stochastic model depends only on the satellite elevation angle and observation geometry, which fails to reflect the impact of atmospheric anisotropy caused by solar radiation on the accuracy of the observations. This results in a mismatch between the statistical characteristics of the weight matrix and the actual observation error, reducing the optimality of the parameter estimation. 3) Under special conditions such as high latitude regions and low solar altitude angle, the spatiotemporal distribution of solar radiation varies more significantly, the tropospheric anisotropy effect is stronger, and the systematic error of traditional models is further amplified, resulting in a significant decrease in positioning accuracy and a significant increase in convergence time, making it difficult to meet the application requirements of high-precision real-time positioning. Summary of the Invention
[0005] To address the aforementioned technical problems, the present invention provides a real-time PPP positioning method that takes into account tropospheric anisotropy, the method comprising the following steps: S01. Receive real-time GNSS observation data, broadcast ephemeris data and precision correction data, and perform preprocessing on the above data, including clock slip and cycle slip detection and repair, data quality control and antenna phase correlation error, to obtain the preprocessed dataset. S02. Calculate the solar altitude angle after atmospheric refraction correction based on the preprocessed datasets and the geographical location and observation time of the corresponding stations for each dataset. S03. Based on the solar elevation angle and the satellite elevation angle, the traditional empirical wet mapping function is anisotropically corrected to obtain an improved wet mapping function that takes into account the anisotropy of the troposphere. S04. Based on the solar elevation angle, anisotropic correction is performed on the traditional stochastic model based on satellite elevation angle to obtain an improved stochastic model for weighting observations. S05. Using the improved wet mapping function and the improved stochastic model, a real-time PPP function model that takes into account the anisotropy of the troposphere is constructed, and Kalman filtering is used for real-time parameter estimation. Finally, the high-precision positioning coordinates of the station are output.
[0006] Preferably, step S02, which involves calculating the solar altitude angle corrected for atmospheric refraction based on the preprocessed datasets and the corresponding station geographical location and observation time, includes the following steps: S21. Calculate solar declination based on the season. Calculated according to the season: ; in: ; S22. Based on the station longitude, Coordinated Universal Time (UTC), and solar angle parameters determined by the season, calculate true solar time, which includes: S221, Longitude correction as follows: ; in, 15 represents the longitude of the station, and 15 represents the longitude span of the Earth's rotation per hour (360° / 24h).
[0007] S222, Time Difference Correction is as follows: ; in, For solar angle parameters; S22, the complete formula for calculating true solar time is: ; in, This represents the true solar time of the station, expressed in hours; it is the fundamental parameter for subsequent calculations of the solar hour angle and solar altitude angle. The station uses Coordinated Universal Time, expressed in hours. This is a longitude correction, expressed in hours; used to correct for deviations between the longitude of the station and Greenwich Mean Time. This is a time difference correction, in hours; the difference between mean solar time and true solar time is caused by the Earth's orbital eccentricity and the obliquity of the ecliptic. S23. Calculate the solar hour angle based on the true solar time, using the following formula: ,when <0 or When the value is greater than 24, modulo 24 processing is required; S24. Based on the station latitude, the solar declination, and the solar hour angle, calculate the original solar altitude angle using the following formula: ,in, The solar altitude angle, The latitude of the measuring station. The solar declination, Solar hour angle; S25. Perform atmospheric refraction correction on the original solar altitude angle, using the following formula: ; in, It is the solar altitude angle after atmospheric refraction correction. This is the original solar altitude angle before correction. It is the angle of atmospheric refraction: The formula for calculating the atmospheric refraction angle is: This yields the final solar altitude angle used for modeling.
[0008] Preferably, the improved wet mapping function in step S03 is obtained in the following way: S31. Obtain the basic form of the traditional GMF wet mapping function: ; in, Here, represents the satellite elevation angle, which is the elevation angle of the GNSS satellite relative to the ground plane of the station. This value is used to characterize the distribution of tropospheric wet delay near the zenith, and varies with station latitude and observation season. The wet delay projection pattern used to adapt to mid-altitude angle satellites is related to the vertical distribution of water vapor at the station. Used to correct wet delay projection bias for low-elevation-angle satellites; S32. Construct an anisotropic correction function The formula is as follows: ; in, The solar altitude angle, The optimal solar altitude angle is usually 45°. The solar altitude angle affects the width parameter, which is usually taken as 15°. The width parameter is affected by the satellite elevation angle, and is usually taken as 30°; S33. The traditional GMF wet mapping function is combined with the anisotropy correction function to form the improved wet mapping function, as shown in the following formula: ; in, This is the anisotropy intensity coefficient, with a value ranging from 0.01 to 0.05, which can be adjusted according to the region and season.
[0009] Preferably, the anisotropy correction function It is the product of a Gaussian function based on the solar elevation angle and a decreasing function based on the satellite elevation angle. Its expression includes adjustable parameters, which include at least the optimal solar elevation angle, the solar elevation angle influence width, and the satellite elevation angle influence width.
[0010] Preferably, an anisotropic intensity coefficient is introduced into the improved wet mapping function expression. An anisotropy intensity coefficient is introduced to control the overall magnitude of the anisotropy correction. Adaptive adjustments are made based on at least one of the following factors: station latitude, season, or local time, wherein: The formula for adjusting the latitude of the S100 station is as follows: ; in, As the baseline value, This is the latitude adjustment factor, with a value ranging from 0.3 to 0.5; S101. Seasonal adjustments, the formula is as follows: ; in This is a seasonal adjustment factor, with a value ranging from 0.1 to 0.2. For reference, the summer solstice is usually chosen. S102. Time adjustment, the formula is as follows: ; in, For local time, This is the daily variation coefficient, with a value range of 0.05-0.1.
[0011] Preferably, the improved stochastic model in step S04 includes the following steps: S41. Obtain the traditional satellite elevation angle-based observation weighting function. ,in, For the prior variance, This is an empirical constant, and its range is 1-3 mm. This is an empirical constant, and its range is 3-5 mm; S42. Construct an anisotropic weight adjustment function. It uses the solar altitude angle obtained in step S02 as the independent variable to adjust the weights of the observed values according to the anisotropic effects caused by solar radiation. The anisotropy weight adjustment ranges from 0.1 to 0.3. The optimal solar altitude angle is 45°. To account for the influence of the solar altitude angle on the width, we take 20°. S43. The traditional observation weight function is modified using the anisotropic weight adjustment function to obtain the improved stochastic model, as shown in the following formula: ,in, The improved observation weights are used for constructing the weight matrix for real-time PPP positioning. Based on weights, For anisotropic adjustment functions, and for functions incorporating solar altitude angles. Satellite elevation angle Quantify the correction of observation weights by tropospheric anisotropy; The variance of the corresponding observed values is: ; in, Let V be the variance of the observed values.
[0012] Preferably, the anisotropic weight adjustment function in step S43 is a Gaussian function with the solar altitude angle as the variable, and its expression includes adjustable parameters, which include at least: anisotropic weight adjustment range, optimal solar altitude angle, and solar altitude angle influence width.
[0013] Preferably, step S04 further includes a multi-system weighting strategy, which includes GPS system, BDS system, Galileo system and GLONASS system, and the relative weights between each system are as follows: ,in, Based on weights, The system characteristic factors are as follows: GPS system BDS system Galileo system and GLONASS system ; By introducing a system characteristic factor into the observation weights of the improved stochastic model, relative weights are assigned between observations from different navigation satellite systems.
[0014] Preferably, in step S05, when using Kalman filtering for real-time parameter estimation, the station's three-dimensional coordinates, receiver clock error, zenith wet delay, and integer ambiguity are used as state parameters, as shown in the following formula: ; in, Let k be the state vector at time k. Here is the state transition matrix. For process noise, For the observation vector, To design the matrix, To observe noise.
[0015] The present invention has at least the following beneficial effects: 1. This method effectively improves the positioning accuracy of real-time PPP, especially in areas and directions where traditional models are weak. By introducing the solar elevation angle parameter and constructing an anisotropic correction function, this method finely characterizes the directional differences in atmospheric refractive index caused by solar radiation, thereby suppressing systematic errors introduced by neglecting anisotropy.
[0016] 2. The improved stochastic model, by incorporating solar elevation angle information, makes the weight allocation of observations more closely reflect the error characteristics caused by actual atmospheric anisotropy, thereby optimizing the weight matrix of the observation equation for parameter estimation. This improves the efficiency and stability of filtering estimation, helping to meet the real-time requirements of autonomous driving, disaster emergency response, and other applications that require rapid acquisition of reliable location information.
[0017] 3. By designing adaptive adjustment strategies for latitude-dependent, seasonal, and time-dependent parameters, the model can automatically adapt to the dynamic changes in solar radiation intensity in different geographical regions, seasons, and within a day, demonstrating strong environmental robustness.
[0018] 4. The existing PPP processing framework only adds a solar altitude angle calculation module, which increases the algorithm complexity to a limited extent and keeps the computational burden under control. Furthermore, the designed multi-system fusion weighting strategy that supports GPS, BDS, Galileo and GLONASS further leverages the observation advantages of each satellite system and enhances the reliability and availability of the positioning service. Attached Figure Description
[0019] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying 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.
[0020] Figure 1 This is a flowchart of a real-time PPP positioning method that takes into account tropospheric anisotropy, provided in Embodiment 1 of the present invention. Detailed Implementation
[0021] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0022] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or server that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or devices.
[0023] Example 1 This embodiment provides a real-time PPP positioning method that takes into account tropospheric anisotropy. The method includes the following steps: Figure 1 As shown: S01. Receive real-time GNSS observation data, broadcast ephemeris data, and precise correction data, and perform preprocessing on the above data, including clock slip and cycle slip detection and repair, data quality control, and antenna phase correlation error, to obtain a preprocessed dataset.
[0024] Specifically, the system synchronously receives real-time observation data (including pseudorange and carrier phase) from GNSS receivers, satellite broadcast ephemeris, and real-time state-space representation (SSR) precise orbit and clock error correction data acquired via data link. Subsequently, the raw observation data undergoes systematic preprocessing to ensure its quality and consistency. This includes: jointly detecting and repairing clock slips at the receiver or satellite using a Melbourne-Wübbena (MW) wide-lane combination and a Geometry-Free (GF) geometry-free combination; reliably detecting and repairing cycle slips based on ionospherically free combination observations combined with phase smoothing pseudorange technology; performing data quality control, such as removing low-quality observations with a signal-to-noise ratio below 36 dB-Hz; and finally, correcting receiver antenna phase center deviation and variations using the igs14.atx standard file published by the International GNSS Service (IGS), and eliminating antenna phase winding errors caused by signal polarization direction changes due to satellite motion through a physical model.
[0025] S02. Based on the preprocessed datasets and the corresponding station locations and observation times, calculate the solar altitude angle corrected for atmospheric refraction, including the following steps: S21. Calculate solar declination based on the season. Calculated according to the season: ; in: ; S22. Based on the station longitude, Coordinated Universal Time (UTC), and solar angle parameters determined by the season, calculate true solar time, which includes: S221, Longitude correction as follows: ; in, 15 represents the longitude of the station, and 15 represents the longitude span of the Earth's rotation per hour (360° / 24h).
[0026] S222, Time Difference Correction is as follows: ; in, For solar angle parameters; S22, the complete formula for calculating true solar time is: ; in, This represents the true solar time of the station, expressed in hours; it is the fundamental parameter for subsequent calculations of the solar hour angle and solar altitude angle. The station uses Coordinated Universal Time, expressed in hours. This is a longitude correction, expressed in hours; used to correct for deviations between the longitude of the station and Greenwich Mean Time. This is a time difference correction, in hours; the difference between mean solar time and true solar time is caused by the Earth's orbital eccentricity and the obliquity of the ecliptic. S23. Calculate the solar hour angle based on true solar time. The formula is as follows: ,when <0 or When the value is greater than 24, modulo 24 processing is required; S24. Based on the station latitude, solar declination, and solar hour angle, calculate the original solar altitude angle using the following formula: ,in, The solar altitude angle, The latitude of the measuring station. The solar declination, Solar hour angle; S25. Atmospheric refraction correction is applied to the original solar altitude angle, using the following formula: ; in, It is the solar altitude angle after atmospheric refraction correction. This is the original solar altitude angle before correction. It is the angle of atmospheric refraction: The formula for calculating the atmospheric refraction angle is: This yields the final solar altitude angle used for modeling.
[0027] Specifically, the solar declination is calculated based on the annual day corresponding to the observation date to characterize the sun's north-south position in the equatorial coordinate system. Next, combining the station's longitude, Coordinated Universal Time (UTC), and the solar angle parameters derived from the annual day, the local true solar time is calculated. This step converts standard UTC into local time reflecting the sun's true apparent motion by introducing longitude corrections and time difference corrections caused by Earth's orbital eccentricity and obliquity of the ecliptic. Then, the solar hour angle is calculated using this true solar time to determine the sun's east-west position relative to the local meridian. Based on this, the station's latitude, the calculated solar declination, and the solar hour angle are substituted into spherical trigonometry formulas to calculate the original solar altitude angle without atmospheric correction. Finally, to obtain an accurate value suitable for precise modeling, atmospheric refraction correction is applied to this original altitude angle to eliminate the light bending effect caused by the vertical gradient of atmospheric density. This ensures that the final solar altitude angle accurately characterizes the solar radiation intensity at the station's location and a specific time, providing reliable key input parameters for subsequent construction of anisotropic mapping functions and stochastic models.
[0028] Taking the high-latitude station MAC1 as an example: Station location: 54.5°S, 158.94°E (MAC1 station), Date: March 8, 2023 (DOY=67, Southern Hemisphere autumn), Time: 08:00 UTC (approximately 18:00 local time), System used: GPS+Galileo+BDS (GEC combination), Meteorological conditions: Clear, stable atmosphere. Then, calculate the solar altitude angle: (I) Calculation of Solar Declination Calculate the annual day parameters using DOY : ; Substituting into the solar declination formula:
[0029] The calculation yields: (In the Southern Hemisphere during autumn, the sun's declination is negative.) (ii) Calculation of True Solar Time Longitude correction:
[0030] Time difference correction (based on) calculate)
[0031] The calculation yields:
[0032] True solar time:
[0033] (III) Calculation of solar hour angle
[0034] (iv) Calculation of solar altitude angle Station latitude (South latitude is negative): The calculation yields: (At low solar altitude angles, anisotropy effects are significant) (v) Atmospheric refraction correction: Atmospheric refraction angle:
[0035] Corrected solar altitude angle: .
[0036] S03. Based on the solar elevation angle and the satellite elevation angle, the traditional empirical wet mapping function is anisotropically corrected to obtain an improved wet mapping function that takes into account the tropospheric anisotropy. The improved wet mapping function described above is obtained in the following way: S31. Obtain the basic form of the traditional GMF wet mapping function: ; in, Here, represents the satellite elevation angle, which is the elevation angle of the GNSS satellite relative to the ground plane of the station. This value is used to characterize the distribution of tropospheric wet delay near the zenith, and varies with station latitude and observation season. The wet delay projection pattern used to adapt to mid-altitude angle satellites is related to the vertical distribution of water vapor at the station. Used to correct wet delay projection bias for low-elevation-angle satellites; S32. Construct an anisotropic correction function The formula is as follows: ; in, The solar altitude angle, The optimal solar altitude angle is usually 45°. The solar altitude angle affects the width parameter, which is usually taken as 15°. The width parameter is affected by the satellite elevation angle, and is usually taken as 30°; S33. Combining the traditional GMF wet mapping function with the anisotropic correction function, an improved wet mapping function is formed, as shown in the following formula: ; in, This is the anisotropy intensity coefficient, with a value ranging from 0.01 to 0.05, which can be adjusted according to the region and season.
[0037] Specifically, based on the satellite elevation angle, the station's geographical location, and the annual day, the basic wet delay projection value (i.e., the basic form of the traditional GMF wet mapping function) is calculated using the traditional GMF wet mapping function model. Secondly, based on the solar elevation angle calculated in step S02 and the current satellite elevation angle, an explicit anisotropy correction function is constructed. As shown in the disclosure document, this function is a product of a Gaussian function and a decreasing function. Its core lies in using both the solar elevation angle and the satellite elevation angle as independent variables: the Gaussian term quantifies the modulation effect of solar radiation intensity (characterized by the solar elevation angle) on the degree of atmospheric anisotropy, with the greatest correction effect when the solar elevation angle is close to the preset optimal value; the decreasing term reflects the influence of the satellite elevation angle, i.e., satellite signals at low elevation angles traverse longer and more complex atmospheric paths, and are more significantly affected by anisotropy, thus requiring a larger correction weight. Finally, the anisotropy correction function is multiplied by an adjustable anisotropy intensity coefficient and linearly combined with the result of the traditional GMF wet mapping function to "form an improved wet mapping function." By using mathematical modeling, the anisotropic physical effects caused by solar radiation are transformed into deterministic corrections to the existing projection function, enabling the final model to respond simultaneously to changes in satellite geometric position and solar radiation state.
[0038] Taking a satellite with an elevation angle of E=20° as an example (satellites with low elevation angles are more affected by anisotropy): ; Based on the GMF model parameters (latitude-dependent), take , , :
[0039] Anisotropic correction function Adopting the modified function form from the patent: ; set up , , : (Low solar altitude angle, moderate correction) Improved wet mapping function: ; Adaptive adjustment based on high latitude: The calculation yields: It can be seen that the wet delay projection value increases by approximately 1.3%, reflecting the additional delay caused by anisotropy under low solar altitude angle conditions.
[0040] S04. Based on the solar elevation angle, the traditional stochastic model based on the satellite elevation angle is anisotropically corrected to obtain an improved stochastic model for weighting observations. Furthermore, the improved stochastic model described above includes the following steps: S41. Obtain the traditional satellite elevation angle-based observation weighting function. ,in, For the prior variance, This is an empirical constant, and its range is 1-3 mm. This is an empirical constant, and its range is 3-5 mm; S42. Construct an anisotropic weight adjustment function. It uses the solar altitude angle obtained in step S02 as the independent variable to adjust the weights of the observed values according to the anisotropic effects caused by solar radiation. The anisotropy weight adjustment ranges from 0.1 to 0.3. The optimal solar altitude angle is 45°. To account for the influence of the solar altitude angle on the width, we take 20°. S43. The traditional observation weight function is modified using the anisotropic weight adjustment function to obtain an improved stochastic model, as shown in the following formula: ,in, The improved observation weights are used for constructing the weight matrix for real-time PPP positioning. Based on weights, For anisotropic adjustment functions, and for functions incorporating solar altitude angles. Satellite elevation angle Quantify the correction of observation weights by tropospheric anisotropy; The variance of the corresponding observed values is: ; in, Let V be the variance of the observed values.
[0041] In step S43 above, the anisotropic weight adjustment function is a Gaussian function with the solar altitude angle as the variable. Its expression includes adjustable parameters, which include at least: the anisotropic weight adjustment range, the optimal solar altitude angle, and the influence width of the solar altitude angle.
[0042] Specifically, a traditional observation weighting function based on satellite elevation angle is used as the base model. This function typically represents a decreasing function of the satellite elevation angle, reflecting the increased observation noise caused by the lengthening of the signal propagation path. Subsequently, an anisotropic weighting adjustment function with the solar elevation angle as the core independent variable is constructed. This function directly uses the solar elevation angle, calculated in step S02 and corrected for atmospheric refraction, as input to quantify the impact of tropospheric anisotropy caused by solar radiation on the accuracy of the observations. Finally, this adjustment function is combined with the traditional weighting function to generate the final improved stochastic model, thus considering both satellite geometry and atmospheric anisotropy factors during the weighting process.
[0043] Furthermore, the anisotropy weight adjustment function described above specifically adopts a Gaussian function form with the solar altitude angle as the variable. This function includes three key adjustable parameters: the anisotropy weight adjustment magnitude, used to control the overall strength of the correction; the optimal solar altitude angle, representing the solar position where the anisotropy effect is most significant; and the influence width of the solar altitude angle, used to describe the sensitivity of the effect to changes in the solar altitude angle. By reasonably setting these parameters, this function can accurately characterize the following physical law: when the solar altitude angle is close to the optimal value, the atmospheric anisotropy effect is strongest, and the observed values are most affected; therefore, its weight should be appropriately reduced. Conversely, when the solar altitude angle deviates significantly, the anisotropy effect weakens, and its influence on the weight also decreases. A detailed example is provided below: Traditional elevation angle weighting, ,set up , , : ; Anisotropic weight adjustment, ,set up , , : ; Improved stochastic models ; Corresponding variance: ; In conclusion, the observation weights were increased by approximately 6.2%, reflecting a decrease in observation accuracy under low solar altitude angle conditions, and therefore a weight reduction process is necessary.
[0044] S05. By using the improved wet mapping function and the improved stochastic model, a real-time PPP function model that takes into account the anisotropy of the troposphere is constructed, and Kalman filtering is used for real-time parameter estimation. Finally, the high-precision positioning coordinates of the station are output.
[0045] When using Kalman filtering for real-time parameter estimation, the station's three-dimensional coordinates, receiver clock error, zenith wet delay, and integer ambiguity are used as state parameters, as shown in the following formula: ; in, Let k be the state vector at time k. Here is the state transition matrix. For process noise, For the observation vector, To design the matrix, To observe noise.
[0046] Specifically, Taking MAC1 station, DOY64-70 as an example, Traditional model RMS: N=7.10cm, E=5.01cm, U=15.80cm Improved model RMS: N=6.86cm, E=4.66cm, U=13.95cm Improvement rate: ; ; Convergence time Convergence time of the traditional model: 18.94 minutes in the U direction; Improved model convergence time: 16.76 minutes in the U direction; Shortened convergence time: ; Verification calculation, , , The improvement was statistically significant (P<0.05).
[0047] Example 2 Based on the above embodiment one, the anisotropic correction function It is the product of a Gaussian function based on the solar elevation angle and a decreasing function based on the satellite elevation angle. Its expression contains adjustable parameters, including at least the optimal solar elevation angle, the solar elevation angle influence width, and the satellite elevation angle influence width.
[0048] Furthermore, an anisotropic intensity coefficient is introduced into the improved wet mapping function expression. An anisotropy intensity coefficient is introduced to control the overall magnitude of the anisotropy correction. Adaptive adjustments are made based on at least one of the following factors: station latitude, season, or local time, wherein: The formula for adjusting the latitude of the S100 station is as follows: ; in, As the baseline value, This is the latitude adjustment factor, with a value ranging from 0.3 to 0.5; S101. Seasonal adjustments, the formula is as follows: ; in This is a seasonal adjustment factor, with a value ranging from 0.1 to 0.2. For reference, the summer solstice is usually chosen. S102. Time adjustment, the formula is as follows: ; in, For local time, This is the daily variation coefficient, with a value range of 0.05-0.1.
[0049] Specifically, an anisotropy correction function is constructed by multiplying the solar elevation angle Gaussian function by the satellite elevation angle decreasing function. The core adjustable parameters of this function provide the basic framework for describing the anisotropic effects under different astronomical and geometric conditions. Subsequently, an anisotropy intensity coefficient is introduced before this correction function and dynamically adjusted through a set of adaptive parameter strategies: for example, adjusting the coefficient baseline value according to the station latitude, seasonally modulating it according to the day of year (DOY), and incorporating local time to reflect diurnal variations. In this way, the final improved wet mapping function not only finely characterizes the spatial distribution of the anisotropic effect through function, but also enables the model to intelligently respond to environmental changes at different geographical locations, seasons, and times of day through adaptive coefficient adjustment. This achieves high-precision, adaptive modeling of tropospheric anisotropy effects and compensates for them in the final delayed projection calculation.
[0050] Example 3 Based on the above embodiment one, step S04 further includes a multi-system weighting strategy, which includes the GPS system, BDS system, Galileo system, and GLONASS system, and the relative weights between each system are as follows: ,in, Based on weights, The system characteristic factors are as follows: GPS system BDS system Galileo system and GLONASS system ; By introducing a system characteristic factor into the observation weights of the improved stochastic model, relative weights can be assigned between observations from different navigation satellite systems.
[0051] Specifically, using the improved wet mapping function and improved stochastic model constructed in the previous steps, an ionospheric-free combined observation equation that takes into account tropospheric anisotropy is constructed. The improved wet mapping function is directly used to calculate the tropospheric delay along the signal propagation path. When entering the parameter estimation stage, a Kalman filter algorithm is used to incorporate parameters such as the station's three-dimensional coordinates, receiver clock error, zenith wet delay (ZWD), and integer ambiguity of each satellite into the state vector for real-time recursive estimation. The most crucial step is that in step S05, when constructing the observation update stage of the Kalman filter, the improved stochastic model established in step S04 must be concretized into an observation weight matrix (or covariance matrix). This means that the weights of the observation values (pseudorange and phase) of each satellite at each epoch are no longer solely determined by the satellite elevation angle, but are dynamically determined by the satellite elevation angle and the real-time calculated solar elevation angle through the improved stochastic model and the anisotropic weight adjustment function. This reflects the spatial differences in observation accuracy caused by atmospheric anisotropy under current solar radiation conditions, enabling the Kalman filter parameter estimation process to more reasonably balance the contributions of observation information from different directions, and ultimately output station coordinates with higher accuracy and faster convergence.
[0052] For example, for GPS satellites (E=20°, h=15.835°): ; For Galileo satellites (same elevation angle): ; The Galileo observations are weighted by 10%, reflecting their higher observation accuracy.
[0053] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.
Claims
1. A real-time PPP positioning method that takes into account tropospheric anisotropy, characterized in that, The method includes the following steps: S01. Receive real-time GNSS observation data, broadcast ephemeris data and precision correction data, and perform preprocessing on the above data, including clock slip and cycle slip detection and repair, data quality control and antenna phase correlation error, to obtain the preprocessed dataset. S02. Calculate the solar altitude angle after atmospheric refraction correction based on the preprocessed datasets and the geographical location and observation time of the corresponding stations for each dataset. S03. Based on the solar elevation angle and the satellite elevation angle, the traditional empirical wet mapping function is anisotropically corrected to obtain an improved wet mapping function that takes into account the anisotropy of the troposphere. S04. Based on the solar elevation angle, anisotropic correction is performed on the traditional stochastic model based on satellite elevation angle to obtain an improved stochastic model for weighting observations. S05. Using the improved wet mapping function and the improved stochastic model, a real-time PPP function model that takes into account the anisotropy of the troposphere is constructed, and Kalman filtering is used for real-time parameter estimation. Finally, the high-precision positioning coordinates of the station are output.
2. The real-time PPP positioning method considering tropospheric anisotropy according to claim 1, characterized in that, Step S02, which involves calculating the solar altitude angle after atmospheric refraction correction based on the preprocessed datasets and the corresponding station geographical location and observation time, includes the following steps: S21. Calculate solar declination based on the season. Calculated according to the season: ； in: ； S22. Based on the station longitude, Coordinated Universal Time (UTC), and solar angle parameters determined by the season, calculate true solar time, which includes: S221, Longitude correction as follows: ； in, 15 represents the longitude of the station, and 15 represents the longitude span of the Earth's rotation per hour (360° / 24h). S222, Time Difference Correction is as follows: ； in, For solar angle parameters; S22, the complete formula for calculating true solar time is: ； in, This represents the true solar time of the station, expressed in hours; it is the fundamental parameter for subsequent calculations of the solar hour angle and solar altitude angle. The station uses Coordinated Universal Time, expressed in hours. This is a longitude correction, expressed in hours; used to correct for deviations between the longitude of the station and Greenwich Mean Time. This is a time difference correction, in hours; the difference between mean solar time and true solar time is caused by the Earth's orbital eccentricity and the obliquity of the ecliptic. S23. Calculate the solar hour angle based on the true solar time, using the following formula: ,when <0 or When the value is greater than 24, modulo 24 processing is required; S24. Based on the station latitude, the solar declination, and the solar hour angle, calculate the original solar altitude angle using the following formula: ,in, The solar altitude angle, The latitude of the measuring station. The solar declination, Solar hour angle; S25. Perform atmospheric refraction correction on the original solar altitude angle, using the following formula: ； in, It is the solar altitude angle after atmospheric refraction correction. This is the original solar altitude angle before correction. It is the angle of atmospheric refraction: The formula for calculating the atmospheric refraction angle is: This yields the final solar altitude angle used for modeling.
3. The real-time PPP positioning method considering tropospheric anisotropy according to claim 1, characterized in that, The improved wet mapping function in step S03 is obtained in the following way: S31. Obtain the basic form of the traditional GMF wet mapping function: ； in, Here, represents the satellite elevation angle, which is the elevation angle of the GNSS satellite relative to the ground plane of the station. This value is used to characterize the distribution of tropospheric wet delay near the zenith, and varies with station latitude and observation season. The wet delay projection pattern used to adapt to mid-altitude angle satellites is related to the vertical distribution of water vapor at the station. Used to correct wet delay projection bias for low-elevation-angle satellites; S32. Construct an anisotropic correction function The formula is as follows: ； in, The solar altitude angle, The optimal solar altitude angle is usually 45°. The solar altitude angle affects the width parameter, which is usually taken as 15°. The width parameter is affected by the satellite elevation angle, and is usually taken as 30°; S33. The traditional GMF wet mapping function is combined with the anisotropy correction function to form the improved wet mapping function, as shown in the following formula: ； in, This is the anisotropy intensity coefficient, with a value ranging from 0.01 to 0.05, which can be adjusted according to the region and season.
4. The real-time PPP positioning method considering tropospheric anisotropy according to claim 3, characterized in that, The anisotropic correction function It is the product of a Gaussian function based on the solar elevation angle and a decreasing function based on the satellite elevation angle. Its expression includes adjustable parameters, which include at least the optimal solar elevation angle, the solar elevation angle influence width, and the satellite elevation angle influence width.
5. A real-time PPP positioning method considering tropospheric anisotropy according to any one of claims 3 or 4, characterized in that, An anisotropy intensity coefficient is introduced into the improved wet mapping function expression. An anisotropy intensity coefficient is introduced to control the overall magnitude of the anisotropy correction. Adaptive adjustments are made based on at least one of the following factors: station latitude, season, or local time, wherein: The formula for adjusting the latitude of the S100 station is as follows: ； in, As the baseline value, This is the latitude adjustment factor, with a value ranging from 0.3 to 0.5; S101. Seasonal adjustments, the formula is as follows: ； in This is a seasonal adjustment factor, with a value ranging from 0.1 to 0.
2. For reference, the summer solstice is usually chosen. S102. Time adjustment, the formula is as follows: ； in, For local time, This is the daily variation coefficient, with a value range of 0.05-0.
1.
6. The real-time PPP positioning method considering tropospheric anisotropy according to claim 1, characterized in that, The improved stochastic model in step S04 includes the following steps: S41. Obtain the traditional satellite elevation angle-based observation weighting function. ,in, For the prior variance, This is an empirical constant, and its range is 1-3 mm. This is an empirical constant, and its range is 3-5 mm; S42. Construct an anisotropic weight adjustment function. It uses the solar altitude angle obtained in step S02 as the independent variable to adjust the weights of the observed values according to the anisotropic effects caused by solar radiation. The anisotropy weight adjustment ranges from 0.1 to 0.
3. The optimal solar altitude angle is 45°. To account for the influence of the solar altitude angle on the width, we take 20°. S43. The traditional observation weight function is modified using the anisotropic weight adjustment function to obtain the improved stochastic model, as shown in the following formula: ,in, The improved observation weights are used for constructing the weight matrix for real-time PPP positioning. Based on the weights, For anisotropic adjustment functions, and for functions incorporating solar altitude angles. Satellite elevation angle Quantify the correction of observation weights by tropospheric anisotropy; The variance of the corresponding observed values is: ； in, Let be the variance of the observed values.
7. A real-time PPP positioning method considering tropospheric anisotropy according to claim 6, characterized in that, The anisotropic weight adjustment function in step S43 is a Gaussian function with the solar altitude angle as the variable. Its expression includes adjustable parameters, which include at least: anisotropic weight adjustment range, optimal solar altitude angle, and solar altitude angle influence width.
8. A real-time PPP positioning method considering tropospheric anisotropy according to claim 1, characterized in that, Step S04 further includes a multi-system weighting strategy, which includes GPS system, BDS system, Galileo system and GLONASS system, and the relative weights between each system are as follows: ,in, Based on the weights, The system characteristic factors are as follows: GPS system BDS system Galileo system and GLONASS system ; By introducing a system characteristic factor into the observation weights of the improved stochastic model, relative weights are assigned between observations from different navigation satellite systems.
9. A real-time PPP positioning method considering tropospheric anisotropy according to claim 1, characterized in that, In step S05, when using Kalman filtering for real-time parameter estimation, the station's three-dimensional coordinates, receiver clock error, zenith wet delay, and integer ambiguity are used as state parameters, as shown in the following formula: ； in, Let k be the state vector at time k. Here is the state transition matrix. For process noise, For the observation vector, To design the matrix, To observe noise.