A method for real-time, rapid and high-precision positioning of beidou navigation satellites based on FY satellites

Abstract

The application discloses a Beidou navigation satellite real-time fast high-precision positioning method based on a Fengyun satellite, first, real-time atmospheric precipitable water and land surface temperature data of the Fengyun satellite are collected, a water vapor conversion coefficient is obtained by using the land surface temperature data, and the atmospheric precipitable water is converted into zenith wet delay through the water vapor conversion coefficient; then, real-time electron density profile data of the Fengyun satellite are collected to obtain total electron content in the zenith direction of the occultation point, and a slant total electron content is obtained by using a Shepard interpolation method and a mapping relationship algorithm; further, a TimeMixer time sequence prediction model is constructed and trained, the zenith wet delay of the current moment is obtained by using the trained TimeMixer model, and the zenith wet delay calculated by the Fengyun satellite is combined with the Beidou receiver position, so that the Beidou receiver preferentially adopts the zenith wet delay calculated by the Fengyun satellite to perform PPP positioning calculation. The method significantly improves the performance of a low-cost Beidou receiver and realizes fast high-precision positioning.

Description

Technical Field

[0001] This invention relates to the field of global satellite navigation meteorology technology, specifically a real-time, rapid, and high-precision positioning method for BeiDou navigation satellites based on Fengyun satellites. Background Technology

[0002] Precise Point Positioning (PPP) technology, a key development direction in the field of satellite navigation and positioning, enables high-precision positioning at the centimeter to millimeter level using only a single receiver. With the full completion of the BeiDou-3 global satellite navigation system, the BeiDou system now possesses the capability to provide global high-precision positioning services. PPP technology, as one of the core technologies in satellite navigation, can achieve positioning accuracy at centimeter to decimeter levels at any location globally, without relying on a ground-based reference station network, and has broad application prospects. However, traditional PPP technology faces a critical bottleneck in practical applications: long convergence time during positioning. This severely limits its promotion in real-time dynamic application scenarios. In the PPP positioning model, ZTD (Zenith Tropospheric Delay) and ionospheric delay are key error sources affecting positioning accuracy and convergence speed. This is because satellite signals are delayed by the troposphere and ionosphere when passing through the atmosphere. In general tropospheric models, ZTD (Zenith Wet Delay) is divided into ZWD (Zenith Wet Delay) and ZHD (Zenith Hydrostatic Delay). ZHD is well correlated with parameters such as receiver elevation, atmospheric pressure, and receiver latitude and longitude, and can be accurately calculated using the Saastamoinen model. However, ZWD is greatly affected by meteorological factors, making calculation and estimation complex. In ionospheric delay models, STEC (Slant Total Electron Content) is the physical carrier of ionospheric delay and is linearly correlated with parameters such as receiver clock error and carrier phase integer ambiguity in the observation file, making it difficult to separate directly. Traditional methods typically estimate these parameters as unknowns in real time during positioning, which greatly increases the dimensionality of the parameter space and prolongs the convergence time. Therefore, improving the convergence speed of real-time PPP is an important research direction in global satellite navigation meteorology.

[0003] Currently, there are roughly three main technical approaches to accelerate PPP convergence. The first is multi-system fusion positioning technology. This technology combines data from multiple navigation systems such as GPS, GLONASS, and Galileo. By increasing the number of visible satellites and improving satellite geometry, it enhances the reliability and convergence speed of positioning solutions. However, differences in signal systems, coordinate frames, and time systems between different satellite systems introduce additional errors, and the system deviation calibration accuracy of low-cost receivers is not as good as that of geodetic receivers. The second is the use of external atmospheric constraints. This method introduces prior tropospheric and ionospheric delay products, treating parameters that would otherwise need to be estimated in PPP solutions, especially ZWD and ionospheric delay, as known quantities or strong constraints, reducing the number of unknown parameters and accelerating the convergence process. However, this method typically uses IGS (International GNSS Service) data, and the low spatial resolution of IGS products leads to a significant decrease in the accuracy of the calculated usable values, failing to meet real-time positioning requirements. The third category is precise modeling methods based on meteorological data. These methods utilize the physical relationship between meteorological observation data and tropospheric delay to establish accurate meteorological-tropospheric models, reducing the number of tropospheric parameters that need to be estimated in PPP calculations and accelerating the convergence process. The core of this approach is to accurately model the ZTD (Zero-Delay Time) using meteorological data provided by ground weather stations and radiosondes, particularly the highly variable ZWD (Zero-Wave Time). However, this method heavily relies on ground weather station data, resulting in limited spatial resolution and timeliness in areas with sparsely distributed ground weather stations.

[0004] The Fengyun-3 (FY-3) satellite, my country's second-generation polar-orbiting meteorological satellite, can broadcast PWV (Precipitable Water Vapor) and VTEC (Vertical Total Electron Content) products in real time, with spatiotemporal resolutions of 5 minutes and 1000 meters, respectively. While the Fengyun satellite possesses real-time capabilities, its orbital characteristics and observation conditions limit its data coverage, resulting in approximately a 15% data gap rate, making it difficult to achieve global, all-weather coverage. ERA5 (European Centre for Medium-Range Weather Forecasts Reanalysis v5) is the fifth-generation reanalysis data from the European Centre for Medium-Range Weather Forecasts, providing data services since 1940. It offers global... Hourly meteorological parameters can be used to calculate high-precision ZWD (Zero-Wide Dynamics). However, ERA5 data has a time delay of approximately 4-6 hours, making true real-time applications impossible. With the development of artificial intelligence, high-precision time series forecasting technology has matured. TimeMixer is a time series forecasting model that employs a multi-scale fusion design. By separating historical data from multi-scale time series, it demonstrates excellent performance in both long-term and short-term forecasting tasks. Using ERA5 reanalysis data to train the TimeMixer time series forecasting model can provide high-precision short-term ZWD forecasts, compensating for the low real-time performance of ERA5 data. Therefore, fully utilizing my country's independent Fengyun meteorological satellite system, internationally available ERA5 reanalysis data, and the latest artificial intelligence algorithms to develop a rapid PPP (Public-Private Partnership) positioning method adapted to domestically produced, low-cost BeiDou hardware has significant strategic importance and application value. Summary of the Invention

[0005] The purpose of this invention is to overcome the shortcomings of the prior art and propose a real-time, fast, and high-precision positioning method for BeiDou navigation satellites based on Fengyun satellites. By introducing multi-source data complementarity and advanced prediction models, stable and reliable high-precision positioning can be achieved in a short time, significantly improving the performance of low-cost BeiDou receivers, realizing fast and high-precision positioning, and providing technical support for the large-scale application of the BeiDou system in the civilian field.

[0006] To achieve the above objectives, the technical solution specifically adopted by the present invention is as follows:

[0007] A real-time, fast, and high-precision positioning method for BeiDou navigation satellites based on Fengyun satellites includes the following steps:

[0008] Step 1: Collect real-time atmospheric precipitable water and land surface temperature data from Fengyun satellites, obtain the water vapor conversion coefficient using the land surface temperature data, convert atmospheric precipitable water into zenith wet delay using the water vapor conversion coefficient, and broadcast the zenith wet delay to the Beidou receiver as a known parameter for PPP positioning calculation.

[0009] Step 2: Collect real-time electron density profile data of Fengyun satellite to obtain the total electron content in the zenith direction of the occultation point. Use Shepard interpolation method and mapping relationship algorithm to obtain the oblique total electron content. Broadcast the oblique total electron content to the Beidou receiver as the initial convergence value for PPP positioning solution.

[0010] Step 3: Collect historical barosphere data and real-time barosphere data for the most recent seven days of ERA5 respectively. Calculate the historical zenith wet delay using the direct integration method, and construct a zenith wet delay sample library and a prediction sample set by combining latitude, longitude, and time information.

[0011] Step 4: Construct a TimeMixer time series prediction model and train the TimeMixer time series prediction model using the zenith wet delay sample library; use the prediction sample set to obtain the zenith wet delay of the current moment through the trained TimeMixer model in a rolling prediction manner, and combine it with the BeiDou receiver position, use bilinear interpolation coupled with elevation correction to calculate the predicted zenith wet delay of the receiver position.

[0012] Step 5: The BeiDou receiver prioritizes using the zenith wet delay calculated by the Fengyun satellite for PPP positioning calculation. When Fengyun satellite data is missing, it switches to the zenith wet delay predicted by ERA5 for calculation. At the same time, it combines the total oblique electron content to complete the PPP positioning calculation, records the convergence time of PPP positioning, saves and outputs the positioning results.

[0013] Preferably, in step 1, the method for obtaining the water vapor conversion coefficient is as follows: using land surface temperature data to calculate the atmospheric weighted average temperature through the Bevis model, and then combining the atmospheric weighted average temperature to calculate the water vapor conversion coefficient.

[0014] Preferably, the Bevis model is a weighted average atmospheric temperature equal to the product of model coefficient a plus model coefficient b and land surface temperature, where model coefficient a takes the value of 70.2 and model coefficient b takes the value of 0.72.

[0015] Preferably, in step 1, the calculation process of the water vapor conversion coefficient is as follows: the water vapor conversion coefficient is equal to 10 to the power of negative 6 multiplied by the liquid water density multiplied by the water vapor gas constant, and then multiplied by the sum of the atmospheric refractive constant K_1 and the atmospheric refractive constant K_2 divided by the weighted average temperature of the atmosphere; wherein the liquid water density is taken as 10^3 kg / m^2, the water vapor gas constant is taken as 461.495 J / (kg⋅K), the atmospheric refractive constant K_1 is taken as 16.48 K / hPa, and the atmospheric refractive constant K_2 is taken as 377600 K^2 / hPa.

[0016] Preferably, in step 2, the Shepard interpolation method calculates the total electron content in the zenith direction at the receiver position. This total electron content in the zenith direction at the receiver position is equal to the sum of the products of the normalized weight function of each occultation point within the influence radius and the total electron content in the zenith direction at the corresponding occultation point. The normalized weight function is the basic weight function of a single occultation point divided by the sum of the basic weight functions of all occultation points within the influence radius. The basic weight function is the difference between the local influence radius of the Earth and the distance between the receiver and the occultation point, divided by the product of the local influence radius of the Earth and the distance between the receiver and the occultation point, and then squared. The local influence radius of the Earth is taken as 400 kilometers.

[0017] Preferably, in step 2, the mapping relationship between the total electron content in the zenith direction and the total electron content in the oblique direction is as follows: the total electron content in the oblique direction is equal to the total electron content in the zenith direction multiplied by 1, then divided by 1 minus the sum of the Earth's radius and the height of the lowest ionospheric layer, then divided by the sum of the Earth's radius, the height of the lowest ionospheric layer, and the height of the ionospheric thin layer from the bottom layer, and then multiplied by the square root of the cosine of the satellite's elevation angle relative to the receiver; wherein the height of the lowest ionospheric layer is taken as 50 kilometers, and the height of the ionospheric thin layer from the bottom layer is taken as an empirical value of 300-450 kilometers.

[0018] Preferably, in step 3, historical barosphere data of ERA5 and real-time barosphere data of ERA5 for the most recent seven days are collected, both of which include air pressure, water vapor pressure, and specific humidity.

[0019] Preferably, in step 3, the direct integration method is used to calculate the zenith wet delay, which is equal to 10 to the power of negative 6 multiplied by the sum of the products of the atmospheric wet refractive index of each layer and the corresponding layer height; wherein the atmospheric wet refractive index is equal to the atmospheric refractive constant K_3 multiplied by the water vapor pressure divided by the land surface temperature, plus the atmospheric refractive constant K_2 multiplied by the water vapor pressure divided by the square of the land surface temperature, the water vapor pressure is equal to the specific humidity multiplied by the air pressure and then divided by 0.622, and the atmospheric refractive constant K_3 is taken as 64.79 K / hPa.

[0020] Preferably, in step 4, the TimeMixer time series forecasting model includes a PDM past decomposable mixing module and an FMM future multi-predictor mixing module. The PDM module decomposes multi-scale time series and mixes seasonal and trend components, while the FMM module integrates multiple predictors to achieve short-term forecasting with zenith wet delay.

[0021] Preferably, in step 4, the rolling prediction method is as follows: the input time series containing N consecutive time nodes is used as the model input to predict the zenith wet delay of the next H time nodes; when new time node data is input, both the input time series and the predicted time series step one time node in the direction of the new data, and the zenith wet delay of the next H time nodes is predicted based on the updated input time series.

[0022] Preferably, in step 4, the bilinear interpolation method coupled with elevation correction is as follows: first, the geodetic height of the receiver is converted into a geopotential height that matches the ERA5 data using the EGM008 model; then, the zenith wet delay of the receiver's horizontal position is calculated using bilinear interpolation; and finally, the elevation is corrected using an elevation attenuation correction factor to obtain the predicted zenith wet delay of the receiver's position. The elevation attenuation correction factor is the negative power of the difference between the geopotential height of the station and the grid center divided by the global average water vapor elevation, where the global average water vapor elevation is 2500 meters.

[0023] Preferably, in step 4, when performing bilinear interpolation, the initial coordinates of the receiver are preferably based on the meter-level coarse positioning results output by the receiver's internal SPS standard positioning service module. If the number of satellites is insufficient or the receiver is powered on for the first time, the historical positioning coordinates stored in the receiver, the coordinates of external auxiliary positioning equipment, or the preset coordinates are used as the initial values.

[0024] Preferably, in step 1, the nearest zenith wet delay value is broadcast to the BeiDou receiver, following the principle of proximity.

[0025] This invention has the following characteristics and beneficial effects:

[0026] The method described in this invention effectively solves the problem of slow convergence speed in real-time PPP positioning using low-cost BeiDou receivers. By extracting ZWD parameters from Fengyun satellite observation data and using this ZWD data as known parameters of tropospheric delay in real-time PPP positioning, the convergence speed can be improved without reducing accuracy. Simultaneously, by extracting ionospheric delay parameters from Fengyun satellite observation data and using them as initial convergence values ​​for the estimated ionospheric parameters in real-time PPP positioning, the convergence speed of real-time positioning can be further improved. Using the ERA5 dataset and training with the TimeMixer model, a short-time ZWD prediction model is constructed, which can compensate for the limited temporal and spatial coverage of a single Fengyun data source. By coupling elevation correction with bilinear interpolation, accurate ZWD at any location is obtained, improving the availability and accuracy of real-time PPP positioning. Using the above algorithm, the slow convergence speed of traditional real-time PPP algorithms in low-cost hardware environments due to the difficulty in estimating tropospheric and ionospheric delay parameters can be effectively solved, while ensuring the real-time performance, accuracy, and stability of real-time positioning. Attached Figure Description

[0027] Figure 1 This is a flowchart of a BeiDou real-time PPP fast and high-precision method based on Fengyun ZWD and ERA5 prediction of ZWD in an embodiment of the present invention.

[0028] Figure 2 This is a flowchart of the ZWD prediction method based on ERA5 data and the TimeMixer model in an embodiment of the present invention.

[0029] Figure 3 This is a flowchart of the BeiDou real-time PPP fast high-precision positioning method based on wind and cloud electron density profile data in an embodiment of the present invention. Detailed Implementation

[0030] The present invention will now be described in detail with reference to specific embodiments. These embodiments will help those skilled in the art to further understand the present invention, but do not limit the invention in any way. It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other.

[0031] A real-time, fast, and high-precision positioning method for BeiDou navigation satellites based on Fengyun satellites, such as... Figure 1 As shown, it includes the following steps:

[0032] Step 1: Collect real-time atmospheric precipitable water and land surface temperature data from Fengyun satellites. Calculate the weighted average atmospheric temperature using the Bevis model based on the land surface temperature data. Then, calculate the water vapor conversion coefficient using the weighted average atmospheric temperature. Convert the atmospheric precipitable water into zenith wet delay using the water vapor conversion coefficient. Broadcast this zenith wet delay to the BeiDou receiver as a known parameter for PPP positioning calculation.

[0033] Specifically, PWV and land surface temperature data products are collected in real time from the Fengyun satellite website. By combining the water vapor conversion coefficient, PWV can be converted into ZWD. The mathematical expression is as follows:

[0034]

[0035]

[0036] In the formula, The water vapor conversion coefficient; Let be the density of liquid water, and take the value of . kg / ; The water vapor gas constant is taken as 461.495 J / (kg). K); and is the atmospheric refractive constant, where The value is 16.48 K / hPa. The value is 377600 / hPa; This is the atmospheric weighted average temperature, expressed in Kelvin (K).

[0037] In the above expression, Land surface temperature data products collected by Fengyun satellites can be used to calculate the temperature using the Bevis model, an atmospheric weighted average temperature model widely used in atmospheric weighted average temperature calculations. Its mathematical expression is as follows:

[0038]

[0039] In the formula, and The coefficients of the Bevis atmospheric weighted average temperature empirical model are as follows: =70.2, =0.72; This refers to the land surface temperature, measured in Kelvin (K).

[0040] The above formula can be used to obtain a high-precision ZWD. Since the resolution of Fengyun data is 1000M, the accuracy is extremely high.

[0041] It should be noted that, considering that the change in ZWD at this scale is negligible, the principle of proximity is adopted to broadcast the ZWD value of the point closest to the receiver to the BeiDou receiver.

[0042] Step 2: Collect real-time electron density profile data of Fengyun satellites to obtain the total electron content in the zenith direction of the occultation point. Use Shepard interpolation to calculate the total electron content in the zenith direction at the Beidou receiver position. Then, calculate the total electron content in the oblique direction through the mapping relationship between the total electron content in the zenith direction and the total electron content in the oblique direction. Broadcast the total electron content in the oblique direction to the Beidou receiver as the initial convergence value for PPP positioning solution.

[0043] Specifically, such as Figure 3 As shown, the electron density profile product is collected in real time from the Fengyun satellite website. This product includes the electron density from near-Earth to satellite orbital altitude, retrieved from the radio occultation signal, i.e., the VTEC at the occultation point. Since the VTEC value is a disordered discrete point in the horizontal direction, the Shepard interpolation method is used to calculate the VTEC at the receiver location. The Shepard interpolation method is suitable for fitting large-scale scattered data and can provide an efficient, stable, and reliable solution strategy for calculating the VTEC at the receiver location. Generally, VTEC is interpolated and STEC is calculated using a single-layer ionospheric model. The principle of the single-layer ionospheric model is to assume that the free electrons along the satellite signal propagation path are concentrated on a thin spherical shell at a certain height, which facilitates calculation. The mathematical expression for calculating the VTEC at the receiver location using the Shepard interpolation method is as follows:

[0044]

[0045] In the formula, It is the VTEC of the interpolation point. These are longitude and latitude, respectively. Indicates the first Normalized weight function for each occultation point; For the first VTEC at each occultation point. The weighting function is... The mathematical expression is as follows:

[0046]

[0047]

[0048] In the formula, Indicates the first The basic weight function for each occultation point; This is the local radius of influence, taken as 400 km. It is the distance between the receiver and the occultation point; To influence the number of occultation points within the radius, its value is related to the formula. The N values ​​are equal. Using the above method, the VTEC value for the corresponding receiver position can be obtained. Then, the STEC is calculated using the mapping relationship between VTEC and STEC, as shown in the following mathematical expression:

[0049]

[0050] In the formula, It is the height of the lowest layer of the ionosphere, taken as 50 km; It is the distance between the thin ionospheric layer and the bottom layer, with an empirical value of 300-450km; It is the Earth's radius; This indicates the altitude angle of the satellite relative to the receiver.

[0051] Through the above calculation process, a high-precision STEC value can be obtained, which is then saved and broadcast to the BeiDou receiver as the initial convergence value for PPP positioning solution.

[0052] Step 3: Collect historical barosphere data and real-time barosphere data for the most recent seven days of ERA5 respectively. Calculate the historical zenith wet delay using the direct integration method, and construct a zenith wet delay sample library and a prediction sample set by combining latitude, longitude, and time information.

[0053] Specifically, a large amount of barosphere dataset was collected from the ERA5 website. This dataset contains essential parameters for calculating ZWD, such as air pressure, vapor pressure, and specific humidity, and has a horizontal resolution of [missing information]. The vertical resolution is 37 layers. The ZWD is calculated using the direct integration method, and its mathematical expression is as follows:

[0054]

[0055] In the formula: It refers to the height of the top floor; It is the receiver elevation; It is the wet refractive index; It is the first The atmospheric wet refractive index of the layer; For the first The height of the floor. The calculation formula is as follows:

[0056]

[0057]

[0058] In the formula: , The atmospheric refractive constant, =377600 / hPa, and The same parameter is used in the middle. =64.79K / hPa; Land surface temperature; It refers to air pressure; It is the water vapor pressure; The specific humidity is given. The accurate ZWD can be obtained using the above formula.

[0059] Using the ZWD obtained from the above formula, a ZWD sample library is constructed by combining latitude and longitude ranges, annual day counts, and UCT time, and stored on a local server.

[0060] Understandably, historical barosphere data from ERA5 can be used to construct a zenith wet delay sample library, while real-time barosphere data from the most recent seven days of ERA5 can be used to construct a prediction sample set.

[0061] It should be noted that the atmospheric refractive constant mentioned in this embodiment... , , The values ​​are all classic and universally accepted values ​​in the field of satellite navigation meteorology. They are industry-recognized parameters that have been verified through long-term theoretical research, atmospheric physics experiments and engineering practice. The core basis for their values ​​comes from the classic definition of the atmospheric refraction physical model and the statistical calibration results of global atmospheric observations. They are also adapted to the physical scenario of tropospheric delay calculation in this invention.

[0062] in, , These two are the core atmospheric refractive constants used in calculating the water vapor conversion coefficient. These values ​​are internationally accepted standard values ​​for calculating the conversion relationship between zenith wet delay (ZWD) and precipitable water volume (PWV), and are adapted to the physical laws governing tropospheric water vapor refraction.

[0063] The values ​​of this set of constants are derived based on the ideal gas law and the physical mechanism of atmospheric refraction, accurately characterizing the refraction effect of water vapor molecules in the atmosphere on satellite electromagnetic signals, as well as the quantitative influence of meteorological elements such as temperature and air pressure on the refraction effect.

[0064] After calibration using atmospheric observation data from radiosondes and satellite remote sensing in different regions and climate conditions around the world, it has become the industry default value for calculating the PWV to ZWD conversion relationship in the fields of satellite navigation and meteorological remote sensing. There is no need for scenario-specific adjustments, and it can be directly adapted to the solution scenarios of converting Fengyun satellite PWV data to ZWD.

[0065] This is a specific refractive constant used for calculating atmospheric wet refractive index. This value is a classic constant for tropospheric wet refractive index models and is suitable for physical scenarios where the direct integration method is used to solve ZWD.

[0066] Its value is derived from the empirical model of atmospheric wet refractive index, which is obtained by fitting a large amount of measured atmospheric refractive index data. It accurately represents the quantitative relationship between water vapor pressure, land surface temperature and atmospheric wet refractive index, and is the core basic parameter for calculating ZWD by direct integration method.

[0067] When used in conjunction with K2, the two parameters characterize the variation of wet refractive index from two dimensions: the first term of water vapor pressure-temperature and the second term of water vapor pressure-temperature. This combination of values ​​is a common choice for organizations such as the European Centre for Medium-Range Weather Forecasts (ECMWF) and the International GNSS Service (IGS) when calculating atmospheric delay, and can be directly adapted to the scenario of calculating ZWD from ERA5 meteorological data.

[0068] Overall, the values ​​of the three atmospheric refraction constants do not require additional scenario-specific calibration. They directly adopt the classic and universal values ​​in the field, which not only conforms to the basic laws of atmospheric physics, but also ensures the consistency of accuracy and reliability of results for the two ZWD solution paths of Fengyun satellite data and ERA5 data in this invention.

[0069] Step 4: Construct a TimeMixer time series prediction model and train the TimeMixer time series prediction model using a zenith wet delay sample library; use the prediction sample set to obtain the zenith wet delay of the current moment through the trained TimeMixer model in a rolling prediction manner; combine the BeiDou receiver position and use bilinear interpolation coupled with elevation correction to calculate the predicted zenith wet delay of the receiver position.

[0070] Specifically, such as Figure 2 As shown, TimeMixer consists of two main modules:

[0071] The PDM (Past-Decomposable-Mixing) module is responsible for decomposing multi-scale sequences and extracting seasonal and trend components from two dimensions, thereby continuously aggregating micro-level seasonal information and macro-level trend information. First, historical data is pooled using average pooling. Perform downsampling By considering multiple scales, a set of multi-scale time series is ultimately obtained:

[0072]

[0073] in , which indicates the first A sequence of scales, where Indicates the first The length of the time step in the scale; Indicates the length of the historical observation sequence; Indicates the number of variables; the lowest level sequence. It is the input sequence, containing the finest temporal variations, while the highest-level sequence... This represents macroscopic changes. Then, in this embodiment, these multi-scale sequences are projected into deep features through an embedding layer. This process can be formalized as follows:

[0074]

[0075] Through the above design, a multi-scale representation of the input sequence is obtained in this embodiment. Next, stacked PDM modules are used to mix historical data at different scales. For the... Layer, input is The PDM process can be formalized as follows:

[0076]

[0077] in, It is the total number of floors, and Indicates having The mixed past representation of the channel.

[0078] The FMM (Future-Multipredictor-Mixing) module further integrates multiple predictors to leverage complementary predictive capabilities from multi-scale observations, thereby improving the accuracy of future sequence predictions. The process is formulated as follows:

[0079]

[0080] in, This indicates the final prediction. This indicates the predicted time series length. Through the above design, TimeMixer can effectively extract key historical information from decomposed multi-scale observations and predict future ZWD data trends accordingly.

[0081] Furthermore, real-time barosphere data from the ERA5 website for the past week is collected as preparation for constructing the ZWD prediction sample set, while the latest barosphere data is continuously collected to supplement it. Then, through the process in step three, ZWD is obtained, and its value is used to construct the ZWD prediction sample set. Unlike the ZWD sample library used for training, the data in the prediction sample set is updated in real time to achieve rolling prediction. Taking time series as an example, the rolling prediction mechanism is as follows:

[0082]

[0083]

[0084] When the latest data is input

[0085]

[0086]

[0087] in, The input data is a time series; The output is the time series of the predicted data; The length of the input sequence; The length of the output sequence; This represents the time series of the latest input data. That is, when the latest data is input, and The entire forecast is incremented by one time unit to the left, ensuring that each prediction is based on all available information up to the current moment, thereby improving the accuracy of the prediction.

[0088] Finally, the predicted ZWD is interpolated to obtain the accurate ZWD for the corresponding receiver position. Bilinear interpolation, as a natural adaptation tool for grid data, significantly outperforms commonly used algorithms such as Kriging interpolation in terms of computational time complexity. By introducing an elevation correction factor, the problem of elevation inconsistencies can be systematically corrected. Furthermore, its accuracy fully meets the requirements of precise positioning, and it is simple, stable, and applicable. Since the ERA5 elevation system is geopotential height while the receiver's elevation system is geodetic height, the receiver elevation needs to be converted to geopotential height before interpolation. This elevation conversion can be achieved using the EGM008 model, which is common knowledge in the field and will not be elaborated further.

[0089] Then, the ZWD (Zero-Wide Position) for the corresponding receiver horizontal position is calculated using bilinear interpolation. This method is a classic algorithm and will not be elaborated further. For the initial receiver coordinates required for interpolation, the meter-level coarse positioning result output by the receiver's internal SPS (Standard Positioning Service) module is preferred. In scenarios with insufficient satellite count or initial power-on, historical positioning coordinates stored in the receiver, external auxiliary positioning equipment, or preset receiver coordinates can be used as initial values. The coarse positioning error has minimal impact on ZWD interpolation, meeting the initial value accuracy requirements for PPP convergence.

[0090] Finally, the horizontal interpolation results are corrected for elevation, and the calculation method is as follows:

[0091]

[0092]

[0093] In the formula, This represents the ZWD value of the receiver position after elevation correction. The ZWD value before elevation correction; Elevation attenuation correction factor; The difference in geopotential between the station and the center of the grid. The global average water vapor elevation is taken as 2500m.

[0094] Using the above method, a high-precision ZWD of the corresponding receiver location can be obtained, saved, and broadcast to the BeiDou receiver. When the Fengyun data is missing, ERA5 data is switched to enter the receiver positioning calculation process.

[0095] Step 5: The BeiDou receiver prioritizes using the zenith wet delay calculated by the Fengyun satellite for PPP positioning calculation. When Fengyun satellite data is missing, it switches to the zenith wet delay predicted by ERA5 for calculation. At the same time, it combines the total oblique electron content to complete the PPP positioning calculation.

[0096] Specifically, the data switching mechanism is as follows:

[0097] The BeiDou receiver performs three core checks on the real-time received Fengyun ZWD data: completeness, spatiotemporal matching, and numerical reasonableness. Only if all three checks pass is the Fengyun ZWD data deemed valid; otherwise, it is directly judged as missing / invalid data, triggering an automatic data source switching mechanism. The specific judgment criteria are as follows:

[0098] Completeness: The Fengyun ZWD data contains latitude and longitude, time series, ZWD values, etc., without any missing or garbled characters; Spatiotemporal matching: The time difference between the data time series and the current observation epoch of the receiver is less than 30 seconds, and the distance between the latitude and longitude of the data and the actual location of the receiver is less than 1000m; Numerical reasonableness: The Fengyun ZWD values ​​are within the normal reasonable range of 0.01~0.5m of global tropospheric wet delay, without any abnormally large or small values.

[0099] When the Fengyun ZWD data is determined to be valid, the receiver substitutes it as a known fixed parameter into the BeiDou PPP positioning model. At the same time, the Fengyun STEC is used as a priori initial convergence value for ionospheric delay to constrain the entire process, reducing the dimension of unknown parameters and constraining the iteration range of ionospheric parameters. When the Fengyun ZWD data is determined to be invalid, only the ZWD of the tropospheric delay term in the positioning solution is replaced with the ZWD predicted by ERA5, while other input parameters do not change with the data source of ZWD.

[0100] Furthermore, the convergence time of PPP positioning is recorded, and the positioning results are saved and output.

[0101] In this invention, the sole criterion for PPP positioning convergence is: the change in the epoch of each component of the three-dimensional coordinates of the BeiDou receiver must be less than 0.01m for three consecutive epochs. The convergence time is recorded with the moment this criterion is achieved as the endpoint. Taking a positioning process as an example, when the BeiDou receiver starts the PPP calculation, an epoch counter is simultaneously activated. Each time the receiver completes an epoch update, it automatically calculates the change in the epoch of each component of the three-dimensional coordinates, monitoring in real time whether the convergence criterion is met. Once met, the counter stops immediately. The receiver's data update frequency is once per second. Multiplying the update frequency by the number of epochs in the counter gives the convergence time for this positioning. The recorded convergence time is used as a core field and bound to key data such as the positioning calculation time series, receiver latitude and longitude, ZWD value, STEC value, and positioning accuracy. It is stored in the receiver's local calculation log in a unified format for easy subsequent analysis.

[0102] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely preferred examples and are not intended to limit the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A real-time, rapid, and high-precision positioning method for BeiDou navigation satellites based on Fengyun satellites, characterized in that, Includes the following steps: Step 1: Collect real-time atmospheric precipitable water and land surface temperature data from Fengyun satellites, obtain the water vapor conversion coefficient using the land surface temperature data, convert atmospheric precipitable water into zenith wet delay using the water vapor conversion coefficient, and broadcast the zenith wet delay to the Beidou receiver as a known parameter for PPP positioning calculation. Step 2: Collect real-time electron density profile data of Fengyun satellite to obtain the total electron content in the zenith direction of the occultation point. Use Shepard interpolation method and mapping relationship algorithm to obtain the oblique total electron content. Broadcast the oblique total electron content to the Beidou receiver as the initial convergence value for PPP positioning solution. The mapping relationship between the total electron content in the zenith direction and the total electron content in the oblique direction is as follows: the total electron content in the oblique direction is equal to the total electron content in the zenith direction multiplied by 1, then divided by 1 minus the sum of the Earth's radius and the height of the lowest ionospheric layer, then divided by the sum of the Earth's radius, the height of the lowest ionospheric layer, and the height of the ionospheric thin layer from the bottom layer, and then multiplied by the square root of the cosine of the satellite's elevation angle relative to the receiver; where the height of the lowest ionospheric layer is taken as 50 kilometers, and the height of the ionospheric thin layer from the bottom layer is taken as an empirical value of 300-450 kilometers; Step 3: Collect historical barosphere data and real-time barosphere data for the most recent seven days of ERA5 respectively. Calculate the historical zenith wet delay using the direct integration method, and construct a zenith wet delay sample library and a prediction sample set by combining latitude, longitude, and time information. Step 4: Construct a TimeMixer time series prediction model and train the TimeMixer time series prediction model using the zenith wet delay sample library; use the prediction sample set to obtain the zenith wet delay of the current moment through the trained TimeMixer model in a rolling prediction manner, and combine it with the BeiDou receiver position, use bilinear interpolation coupled with elevation correction to calculate the predicted zenith wet delay of the receiver position. The TimeMixer time series forecasting model includes a PDM past decomposable mixing module and an FMM future multi-predictor mixing module. The PDM module decomposes multi-scale time series and mixes seasonal and trend components, while the FMM module integrates multiple predictors to achieve short-term forecasting with zenith wet delay. Step 5: The BeiDou receiver prioritizes using the zenith wet delay calculated by the Fengyun satellite for PPP positioning calculation. When Fengyun satellite data is missing, it switches to the zenith wet delay predicted by ERA5 for calculation. At the same time, it combines the total oblique electron content to complete the PPP positioning calculation, records the convergence time of PPP positioning, saves and outputs the positioning results.

2. The method according to claim 1, characterized in that, In step 1, the method for obtaining the water vapor conversion coefficient is as follows: using land surface temperature data to calculate the atmospheric weighted average temperature through the Bevis model, and then combining the atmospheric weighted average temperature to calculate the water vapor conversion coefficient.

3. The method according to claim 2, characterized in that, The Bevis model states that the atmospheric weighted average temperature equals the product of the model coefficient a plus the model coefficient b and the land surface temperature, where model coefficient a takes a value of 70.2 and model coefficient b takes a value of 0.

72.

4. The method according to claim 3, characterized in that, In step 1, the calculation process of the water vapor conversion coefficient is as follows: the water vapor conversion coefficient is equal to 10 to the power of negative 6 multiplied by the liquid water density multiplied by the water vapor gas constant, and then multiplied by the sum of the atmospheric refractive constant and the atmospheric refractive constant divided by the atmospheric weighted average temperature; where the liquid water density is taken as kg / , the water vapor gas constant is taken as 461.495 J / (kgK), the atmospheric refractive constant is taken as 16.48 K / hPa, and the atmospheric refractive constant is taken as 377600 / hPa.

5. The method according to claim 1, characterized in that, In step 2, the Shepard interpolation method calculates the total electron content in the zenith direction at the receiver location. The total electron content in the zenith direction at the receiver location is equal to the sum of the products of the normalized weight function of each occultation point within the influence radius and the total electron content in the zenith direction of the corresponding occultation point. The normalized weight function is the basic weight function of a single occultation point divided by the sum of the basic weight functions of all occultation points within the influence radius. The basic weight function is the difference between the local influence radius of the Earth and the distance between the receiver and the occultation point, divided by the product of the local influence radius of the Earth and the distance between the receiver and the occultation point, and then squared. The local influence radius of the Earth is taken as 400 kilometers.

6. The method according to claim 1, characterized in that, In step 3, historical barosphere data of ERA5 and real-time barosphere data of ERA5 for the past seven days are collected, both of which include air pressure, water vapor pressure, and specific humidity.

7. The method according to claim 6, characterized in that, In step 3, the direct integration method is used to calculate the zenith wet delay, which is equal to 10 to the power of negative 6 multiplied by the sum of the products of the atmospheric wet refractive index of each layer and the corresponding layer height. The atmospheric wet refractive index is equal to the atmospheric refractive constant multiplied by the water vapor pressure divided by the land surface temperature, plus the atmospheric refractive constant multiplied by the water vapor pressure divided by the square of the land surface temperature. The water vapor pressure is equal to the specific humidity multiplied by the air pressure and then divided by 0.

622. The atmospheric refractive constant is taken as 64.79 K / hPa.

8. The method according to claim 1, characterized in that, In step 4, the rolling prediction method is as follows: the input time series containing N consecutive time nodes is used as the model input to predict the zenith wet delay of the next H time nodes; when new time node data is input, both the input time series and the predicted time series step one time node in the direction of the new data, and the zenith wet delay of the next H time nodes is predicted based on the updated input time series.

9. The method according to claim 1, characterized in that, In step 4, the bilinear interpolation method coupled with elevation correction is as follows: first, the earth height of the receiver is converted into a potential height that matches the ERA5 data through the EGM008 model; then, the zenith wet delay of the receiver's horizontal position is calculated using the bilinear interpolation method; and finally, the elevation is corrected by the elevation attenuation correction factor to obtain the predicted zenith wet delay of the receiver position. The elevation attenuation correction factor is the negative of the natural constant e, which is the difference between the geopotential height between the station and the grid center divided by the global average water vapor elevation, which is taken as 2500 meters.

10. The method according to claim 1, characterized in that, In step 4, when performing bilinear interpolation, the initial coordinates of the receiver are preferentially based on the meter-level coarse positioning results output by the receiver's internal SPS standard positioning service module. When the number of satellites is insufficient or when the receiver is powered on for the first time, the historical positioning coordinates stored in the receiver, the coordinates of external auxiliary positioning equipment, or the preset coordinates are used as the initial values.

11. The method according to claim 1, characterized in that, In step 1, the nearest zenith wet delay value is broadcast to the BeiDou receiver, based on the principle of proximity.

Images