A direct retrieval method of GNSS precipitable water vapor based on spatial sliding window grid

Abstract

The application discloses a GNSS precipitable water vapor direct inversion method based on a spatial sliding window grid, and comprises the following steps: obtaining historical reanalysis data of a research area, and calculating zenith wet delay sequence and atmospheric precipitable water vapor sequence; adopting a spatial sliding window strategy to perform grid processing on the research area; according to the zenith wet delay sequence, the atmospheric precipitable water vapor sequence and the grid processing result, a periodic regression model containing seasonal harmonic terms is used to construct a spatial sliding window direct grid model; obtaining real-time GNSS zenith wet delay of a target station, and obtaining atmospheric precipitable water vapor prediction results through grid interpolation. The application does not need real-time meteorological parameters, effectively captures the spatial inhomogeneity and time variability of atmospheric water vapor through the spatial sliding window, and is suitable for high-precision real-time inversion in a meteorological data-deficient area.

Description

Technical Field

[0001] This invention belongs to the field of GNSS meteorology, and in particular relates to a direct inversion method for GNSS precipitable water based on a spatial sliding window grid. Background Technology

[0002] Global Navigation Satellite System (GNSS) meteorology has become the mainstream technique for retrieving atmospheric precipitable water. Traditional GNSS water vapor retrieval methods typically employ a two-step approach: first, the total zenith delay is obtained through precise point positioning or other GNSS data processing techniques; then, a model (such as the Saastamonin model) is used to separate the zenith static delay, yielding the zenith wet delay; finally, the zenith wet delay is converted into atmospheric precipitable water using transformation parameters (usually related to the weighted average temperature). To obtain these transformation parameters, existing techniques largely rely on real-time observational data from ground-based meteorological sensors (such as barometers and thermometers) or estimate the weighted average temperature using global empirical models (such as global pressure-temperature series models). In addition, some studies have attempted to establish direct empirical regression models between zenith wet delay and atmospheric precipitable water, or to use artificial intelligence methods such as deep learning for end-to-end retrieval.

[0003] However, existing technologies still have the following prominent problems: First, traditional two-step methods and direct empirical models generally rely on real-time meteorological parameters or external numerical weather prediction data, making them difficult to apply in areas with scarce meteorological data (such as uninhabited areas, oceans, and polar regions). Furthermore, data latency and hardware costs limit the potential for real-time high-frequency inversion. Second, while existing global grid models can provide transformation parameters in the absence of on-site meteorological observations, their spatial resolution is limited, making it difficult to capture the spatial heterogeneity of atmospheric water vapor caused by complex terrain and local climate characteristics. Third, static linear models or global spherical harmonic function models fail to effectively fit the nonlinear time-varying characteristics of atmospheric precipitable water under the influence of seasonal variations and atmospheric circulation cycle fluctuations. While deep learning methods can learn complex nonlinear relationships, they require massive amounts of training data and expensive computational resources, have poor physical interpretability, and are not conducive to operational deployment. Summary of the Invention

[0004] The purpose of this invention is to provide a direct inversion method for GNSS precipitable water based on a spatial sliding window grid, so as to solve the technical problems in the prior art, such as the reliance on real-time meteorological parameters for PWV inversion, the lack of local adaptability of global grid models, and the difficulty of static models to capture spatiotemporal changes.

[0005] To achieve the above objectives, this invention provides a method for direct inversion of GNSS precipitable water based on a spatial sliding window grid, comprising: Historical reanalysis data of the study area are obtained, and the zenith wet delay sequence and atmospheric precipitable water sequence are calculated based on the historical reanalysis data. Based on the zenith wet delay sequence and the atmospheric precipitable water sequence, a spatial sliding window direct grid model is constructed using a spatial sliding window strategy and a periodic regression model containing seasonal harmonic terms. The real-time GNSS zenith wet delay of the target station is obtained. Based on the real-time GNSS zenith wet delay and the spatial sliding window direct grid model, the atmospheric precipitable water prediction result of the target station is obtained by grid interpolation.

[0006] Preferably, the process of obtaining the zenith wet delay sequence and the atmospheric precipitable water sequence includes: Based on the specific humidity data and air pressure data in the historical reanalysis data, the atmospheric precipitable water sequence is obtained by vertical integration calculation. Based on the geopotential, temperature, specific humidity, and relative humidity data in the historical reanalysis data, the total zenith retardation sequence is obtained by integrating the atmospheric refractive index layer by layer, and the wet zenith retardation sequence is calculated based on the total zenith retardation sequence and the zenith static retardation sequence.

[0007] Preferably, the process of calculating the atmospheric precipitable water sequence includes: The calculation was performed using specific humidity data and air pressure data through a vertical integration formula, where the lower limit of integration is the surface air pressure and the upper limit is the top atmospheric pressure. The trapezoidal rule was used to perform numerical calculations between each pressure layer. The process of calculating the zenith wet delay sequence includes: The total zenith delay is obtained by integrating the atmospheric refractive index layer by layer. The zenith static delay is calculated using the Sastamonin model based on surface air pressure, station latitude, and ellipsoidal height. The wet zenith delay is obtained by subtracting the zenith static delay from the total zenith delay.

[0008] Preferably, in the periodic regression model used to construct the spatial sliding window direct grid model, atmospheric precipitable water is expressed as the sum of the following terms: The terms include the intercept term, the product of the zenith wet delay and the linear transformation coefficient, the annual cycle cosine term, the annual cycle sine term, the semi-annual cycle cosine term, and the semi-annual cycle sine term, where the annual cycle term and the semi-annual cycle term both use the accumulated days of the year as independent variables.

[0009] Preferably, the method further includes: performing gridding processing on the study area using the spatial sliding window strategy; The process of gridding the study area using the aforementioned spatial sliding window strategy includes: The study area is divided into grids based on the resolution of the historical reanalysis data; Determine the size and step size of the sliding window; Starting from the boundary of the study area, multiple sliding windows are defined sequentially, with each sliding window covering multiple grid points.

[0010] Preferably, the process of constructing the spatial sliding window direct grid model includes: For each sliding window, aggregate the zenith wet delay sequence and atmospheric precipitable water sequence data of all grid points within the window; The parameters of the periodic regression model are fitted using the least squares method to obtain the model parameters for each sliding window; Assign the model parameters of each sliding window to the center grid point of the corresponding sliding window; The model parameters of all sliding window center grid points are reorganized to generate multiple coefficient grid diagrams.

[0011] Preferably, the plurality of coefficient grid plots include an intercept grid, a linear transformation coefficient grid, an annual periodic cosine coefficient grid, an annual periodic sine coefficient grid, a semi-annual periodic cosine coefficient grid, and a semi-annual periodic sine coefficient grid, and the spatial resolution of each grid plot is consistent with the resolution of the historical reanalysis data.

[0012] Preferably, the process of obtaining the real-time GNSS zenith wet delay of the target station includes: Obtain the real-time total zenith delay of the target station; Obtain the predicted surface air pressure value of the target station, and calculate the zenith static delay based on the predicted surface air pressure value; The real-time GNSS wet zenith delay is calculated based on the real-time total zenith delay and the zenith static delay.

[0013] Preferably, the process of obtaining the atmospheric precipitable water prediction result by grid interpolation based on the real-time GNSS zenith wet delay and the spatial sliding window direct grid model includes: Based on the latitude and longitude coordinates of the target station, the model parameters for the corresponding location are obtained by interpolation from multiple coefficient grid diagrams; Obtain the cumulative day of the year corresponding to the current observation time; The real-time GNSS zenith wet delay, the annual cumulative day, and the model parameters are substituted into the periodic regression model to calculate the predicted atmospheric precipitable water.

[0014] Preferably, the method further includes evaluating the accuracy of the spatial sliding window direct grid model; The process of evaluating the accuracy of the spatial sliding window direct grid model: Using historical reanalysis data from the test period, the zenith wet delay of the test period is input into the spatial sliding window direct grid model to obtain the predicted atmospheric precipitable water. The predicted atmospheric precipitable water is then compared with the atmospheric precipitable water of the test period, and the root mean square error and bias are calculated to verify internal consistency. Using GNSS observation zenith wet delay as input, the atmospheric precipitable water inverted by the spatial sliding window direct grid model is compared with the GNSS atmospheric precipitable water inverted by traditional methods and the radiosonde atmospheric precipitable water, to verify external accuracy. The zenith wet delay, calculated in real time using precise single-point positioning technology, is used as input to evaluate the inversion accuracy of the spatial sliding window direct grid model in real-time scenarios.

[0015] Compared with the prior art, the present invention has the following advantages and technical effects: First, by pre-establishing a coefficient grid based on historical reanalysis data, this invention enables real-time atmospheric precipitable water inversion without the need for any ground meteorological sensors or external numerical weather forecast data during the real-time application phase, thus overcoming the bottleneck of applicability in areas with scarce meteorological data.

[0016] Second, the present invention uses a spatial sliding window strategy for grid modeling. Each window aggregates data from neighboring grid points and obtains location-specific transformation coefficients through least squares fitting. This effectively captures the spatial heterogeneity of atmospheric water vapor caused by changes in elevation, latitude, and topography, and significantly improves the local adaptability of complex regions.

[0017] Third, this invention introduces annual and semi-annual harmonic terms into the regression model, which accurately characterizes the seasonal time-varying features of the conversion efficiency and makes up for the shortcomings of static linear models in fitting the nonlinear time variation of atmospheric precipitable water. At the same time, the model has clear physical interpretability, high computational efficiency, and is easy to promote in business applications.

[0018] Fourth, verified by various accuracy evaluation indicators, the present invention exhibits excellent inversion performance in terms of internal consistency, external accuracy, and real-time precise single-point positioning scenarios. It is highly consistent with radiosonde data and traditional methods, and can meet the high-timeliness application requirements such as extreme weather monitoring and flash flood forecasting. Attached Figure Description

[0019] The accompanying drawings, which form part of this application, are used to provide a further understanding of this application. The illustrative embodiments and descriptions of this application are used to explain this application and do not constitute an undue limitation of this application. In the drawings: Figure 1 This is a schematic diagram of the method flow according to an embodiment of the present invention; Figure 2 This is a schematic diagram of the spatial sliding window algorithm according to an embodiment of the present invention; Figure 3 This is a spatial distribution diagram of RMSE and Bias of the SSW-DGM model in the fitting and prediction periods according to an embodiment of the present invention. Figure 4 Site-specific RMSE, Bias, and maximum absolute residual maps for 40 GNSS stations in this embodiment of the invention; Figure 5 The spiral bar chart for deriving PWV from SSW-DGM of GNSS-ZWD observations in this embodiment of the invention; Figure 6 This is a scatter plot showing the correlation between the PWV inverted from the representative GNSS station SSW-DGM and the ERA5-PWV in this embodiment of the invention. Figure 7 This is a correlation diagram between SSW-DGM inverted PWV and co-located radiosonde PWV in an embodiment of the present invention. Figure 8 This is a monthly average RMSE, Bias, and maximum residual map of PWV inversion for all GNSS stations using SSW-DGM in this embodiment of the invention. Figure 9 This is a comparison chart of the statistical performance indicators of the SSW-DGM model in four seasons according to an embodiment of the present invention. Figure 10 This is a spatial relationship diagram between the BJFS GNSS station and the surrounding SSW-DGM grid nodes in an embodiment of the present invention. Figure 11 This paper presents a correlation analysis and residual distribution diagram of the SSW-DGM inversion of BJFS sites based on PPP and traditional methods, as shown in this embodiment of the invention. Detailed Implementation

[0020] It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0021] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0022] like Figure 1 As shown, this embodiment provides a method for direct inversion of GNSS precipitable water based on a spatial sliding window grid, including: Historical reanalysis data of the study area were obtained, and the zenith wet delay sequence and atmospheric precipitable water sequence were calculated based on the historical reanalysis data. Based on the zenith wet delay sequence and atmospheric precipitable water sequence, a spatial sliding window direct grid model was constructed using a spatial sliding window strategy and a periodic regression model containing seasonal harmonic terms. The real-time GNSS zenith wet delay of the target station is obtained. Based on the real-time GNSS zenith wet delay and the spatial sliding window direct grid model, the atmospheric precipitable water prediction result of the target station is obtained by grid interpolation.

[0023] Furthermore, this embodiment acquires ERA5 reanalysis data of the study area, including geopotential, temperature, specific humidity, and relative humidity of the barosphere. Zenith wet delay (ZWD) and precipitable water vapor (PWV) are calculated through integration. The conversion relationship between ZWD and PWV is established, and a periodic regression model incorporating seasonal harmonic terms is constructed. A spatial sliding window strategy is used to grid the study area. The least squares method is used to fit the model parameters within each sliding window, constructing a spatial sliding window direct grid model (SSW-DGM) and generating a high-resolution coefficient grid map. Using the constructed SSW-DGM model, real-time GNSS-ZWD and spatiotemporal information from the target station are input, and the predicted precipitable water vapor is directly output through grid interpolation. The accuracy of the inversion model is evaluated using accuracy evaluation indicators to verify the model's internal consistency, external accuracy, and real-time applicability.

[0024] This embodiment calculates ZWD and PWV by acquiring ERA5 reanalysis data; establishes a ZWD-PWV direct regression model including seasonal harmonic terms; uses a spatial sliding window strategy to grid the study area, generating location-specific transformation coefficients for each grid node; constructs a spatial sliding window direct grid model (SSW-DGM) using the least squares method to generate a high-resolution coefficient grid map; inputs real-time GNSS-ZWD and directly inverts PWV through grid interpolation; and evaluates accuracy using RMSE, bias, and correlation coefficients. This embodiment does not require real-time meteorological parameters, effectively capturing spatial heterogeneity and temporal variability. In the Loess Plateau validation, the internal RMSE is 0.23 mm, highly consistent with traditional methods and radiosonde data (R>0.98). 91.4% of the results meet the standards in the real-time PPP scenario, making it suitable for high-precision real-time PWV inversion in areas with scarce meteorological data.

[0025] Furthermore, the process of obtaining the zenith wet delay sequence and the atmospheric precipitable water sequence includes: Based on the specific humidity and air pressure data in the historical reanalysis data, the atmospheric precipitable water sequence is obtained by vertical integration calculation. Based on the geopotential, temperature, specific humidity, and relative humidity data in the historical reanalysis data, the total zenith retardation sequence is obtained by integrating the atmospheric refractive index layer by layer. Based on the total zenith retardation sequence and the zenith static retardation sequence, the wet zenith retardation sequence is calculated.

[0026] Furthermore, the process of calculating the atmospheric precipitable water sequence includes: The calculation was performed using specific humidity data and air pressure data through a vertical integration formula, where the lower limit of integration is the surface air pressure and the upper limit is the top atmospheric pressure. The trapezoidal rule was used to perform numerical calculations between each pressure layer. The process of calculating the zenith wet delay sequence includes: The total zenith delay is obtained by integrating the atmospheric refractive index layer by layer. The zenith static delay is calculated using the Sastamonin model based on surface air pressure, station latitude, and ellipsoidal height. The wet zenith delay is obtained by subtracting the zenith static delay from the total zenith delay.

[0027] Furthermore, this embodiment calculates ERA5_PWV using meteorological data, and the formula is as follows: In the formula, q represents the specific humidity (g / kg), and g is the acceleration due to gravity (m / s²). 2 ), where p is the air pressure (hpa), The density of liquid water (kg / m³) P surf Surface air pressure, P top The integral is performed numerically on the pressure layer provided by ERA5, representing the top atmospheric pressure. The zenith total time delay (ZTD) was calculated using ERA5 meteorological data, obtained by integrating the atmospheric refractive index N layer by layer: In the formula, h top The height of the upper level h surf The height of the lower layer, where the refractive index is... N The following is given by the Smith-Weintraub equation: In the formula, k 1 = 77.604 (K / hpa), k 2 = 64.79 (K / hpa), k 3 = 377600.0 (K2 / hpa), P dThe partial pressure of dry air (hPa) T Thermodynamic temperature (K), e The water vapor pressure (hPa) is obtained by passing through the saturated water vapor pressure. P sa relative humidity RH and wet q The calculation yields the following formula: Zenith static delay was calculated using the Saastamoinen model. ZHD : In the formula, P s φ represents the surface air pressure (hPa) and the latitude of the station (rad). h The height of the ellipsoid is in km. Calculate the zenith wet delay (ZWD): Through the above process, ZWD and PWV time series data for the training period (2016-2019) were generated for all ERA5 grid points in the study area.

[0028] Furthermore, in the periodic regression model used to construct the spatial sliding window direct grid model, atmospheric precipitable water is expressed as the sum of the following terms: The terms include the intercept term, the product of the zenith wet delay and the linear transformation coefficient, the annual cycle cosine term, the annual cycle sine term, the semi-annual cycle cosine term, and the semi-annual cycle sine term, where the annual cycle term and the semi-annual cycle term both use the accumulated days of the year as independent variables.

[0029] Furthermore, this embodiment establishes a periodic regression model for the direct conversion of ZWD-PWV, as shown in the following formula: In the formula, PWV and ZWD The units are all in mm. a 0 is the intercept. a 1 represents the primary linear transformation coefficient. doy For years, A 1. B 1 is the coefficient for the annual periodic term. A 2. B 2 represents the coefficient for the semi-annual periodic term; This model captures harmonic terms, where aThe coefficient 1 acts as a dimensionless constant Π that changes dynamically, eliminating the need for explicit calculation of Tm and thus bypassing the dependence on meteorological parameters in the traditional two-step method.

[0030] Furthermore, the method also includes: using a spatial sliding window strategy to grid the study area; The process of gridding the study area using a spatial sliding window strategy includes: The study area is divided into grids based on the resolution of historical reanalysis data; Determine the size and step size of the sliding window; Starting from the boundary of the study area, multiple sliding windows are defined sequentially, with each sliding window covering multiple grid points.

[0031] Furthermore, in this embodiment, the study area is divided into grids according to the ERA5 data resolution (0.25°×0.25°), covering all grid points in the study area; The sliding window size was determined to be 5×5 grid points, corresponding to a spatial range of approximately 1.25°×1.25°, to ensure that the window contains enough samples for parameter solving, while maintaining sensitivity to local atmospheric characteristics; Starting from the top left corner of the study area, define the first window G1, which covers 5×5 grid points; Set the sliding step size to 5 grid points to ensure that the center points of adjacent windows do not overlap, thus achieving complete spatial coverage.

[0032] Furthermore, the process of constructing a spatial sliding window direct grid model includes: For each sliding window, aggregate the zenith wet delay sequence and atmospheric precipitable water sequence data of all grid points within the window; The parameters of the periodic regression model are fitted using the least squares method to obtain the model parameters for each sliding window. Assign the model parameters of each sliding window to the center grid point of the corresponding sliding window; The model parameters of all sliding window center grid points are reorganized to generate multiple coefficient grid diagrams.

[0033] Furthermore, multiple coefficient grid plots include intercept grids, linear transformation coefficient grids, annual cycle cosine coefficient grids, annual cycle sine coefficient grids, semi-annual cycle cosine coefficient grids, and semi-annual cycle sine coefficient grids, with the spatial resolution of each grid plot consistent with the resolution of the historical reanalysis data.

[0034] Furthermore, in this embodiment, for each sliding window, the ZWD and PWV time series data of all 25 grid points in the aggregation window are guaranteed; The model parameters of this window are solved using least squares estimation. a0, a 1, A 1, B 1, A 2, B 2) Determine the optimal model parameters based on minimizing the sum of squared predicted residuals for all grid points and all time points within the window; Assign the solved parameters to the center grid points of the window; Move the window 5 grid points eastward and repeat the parameter estimation process until the current latitude zone is completed; Move the window 5 grid points south to the next latitude zone and repeat the above process; Iterate until the entire study area is covered, then reorganize the parameters at the center points of all windows to form six high-resolution coefficient grid diagrams, which are respectively... a 0 grid, a 1 grid A 1 grid B 1 grid A 2 grid and B Two grids, with each grid having the same spatial resolution as ERA5.

[0035] Furthermore, the process of obtaining the real-time GNSS zenith wet delay of the target station includes: Obtain the real-time total zenith delay of the target station; Obtain the predicted surface air pressure value of the target station, and calculate the zenith static delay based on the predicted surface air pressure value; The real-time GNSS wet zenith delay is calculated based on the real-time total zenith delay and the zenith static delay.

[0036] Furthermore, the process of obtaining atmospheric precipitable water prediction results through grid interpolation based on the real-time GNSS zenith wet delay and spatial sliding window direct grid model includes: Based on the latitude and longitude coordinates of the target station, the model parameters for the corresponding location are obtained by interpolation from multiple coefficient grid maps; Obtain the cumulative day of the year corresponding to the current observation time; By substituting real-time GNSS zenith wet delay, annual cumulative day, and model parameters into a periodic regression model, the predicted atmospheric precipitable water content is obtained.

[0037] Furthermore, in this embodiment, the real-time total zenith delay (ZTD) of the target GNSS station can be obtained in real time through precise point positioning (PPP) technology or network RTK technology. The GPT3 model was used to obtain the predicted surface air pressure at the station location, and the zenith static delay ZHD was calculated. Then, the zenith wet delay ZWD was obtained. Based on the latitude and longitude coordinates of the station, the model parameters for that location are obtained from a pre-established grid map of six coefficients using bilinear interpolation. a 0, a 1, A 1, B 1, A 2, B 2); Obtain the year-day corresponding to the current observation time. doy ZWD, doy The coefficients obtained by interpolation are substituted into the periodic regression model formula to directly calculate the output PWV value.

[0038] Furthermore, the method also includes accuracy evaluation of the spatial sliding window direct grid model; The process of evaluating the accuracy of a spatial sliding window direct grid model: Using historical reanalysis data from the test period, the zenith wet delay of the test period is input into the spatial sliding window direct grid model to obtain the predicted atmospheric precipitable water. The predicted atmospheric precipitable water is then compared with the atmospheric precipitable water of the test period, and the root mean square error and bias are calculated to verify internal consistency. Using GNSS observation zenith wet delay as input, the atmospheric precipitable water inverted by the spatial sliding window direct grid model is compared with the GNSS atmospheric precipitable water inverted by traditional methods and the radiosonde atmospheric precipitable water, to verify external accuracy. The zenith wet delay calculated in real time using precise single-point positioning technology is used as input to evaluate the inversion accuracy of the spatial sliding window direct grid model in real-time scenarios.

[0039] Furthermore, the accuracy evaluation in this embodiment specifically includes the following steps: Internal consistency verification was performed using retained test period (2020) ERA5 data. The ERA5_ZWD was input into the SSW-DGM model, the predicted PWV was calculated, and compared with the ERA5_PWV. The RMSE and Bias were calculated to evaluate the mathematical fidelity of the model. External accuracy verification: Using GNSS observation ZWD as input, the PWV inverted by SSW-DGM is compared with the GNSS-PWV inverted by traditional methods and the PWV of the same radiosonde, to evaluate the reliability of the model on actual observation data. Real-time applicability verification: Using the ZWD calculated in real time by PPP technology as input, the inversion accuracy and stability of SSW-DGM in real-time scenarios are evaluated. The formula for calculating the accuracy evaluation index is as follows: In the formula, xi For reference only. yi For predicted values, and These represent the average values ​​of their respective samples. N This represents the number of valid samples.

[0040] As a preferred implementation method, this embodiment provides a direct inversion method for GNSS precipitable water based on a spatial sliding window grid, including the following steps: Step 1: Obtain ERA5 reanalysis data for the study area, including geopotential, temperature, specific humidity, and relative humidity of the barosphere. Calculate the Zenith Wet Delay (ZWD) and Precipitable Water Vapor (PWV) through integration. Step 2: Establish the conversion relationship between ZWD and PWV, and construct a periodic regression model that includes seasonal harmonic terms; Step 3: Grid the study area using a spatial sliding window strategy; Step 4: Fit the model parameters within each sliding window using the least squares method to construct the Spatial Sliding Window Direct Grid Model (SSW-DGM) and generate a high-resolution coefficient grid diagram. Step 5: Using the constructed SSW-DGM model, input the real-time GNSS-ZWD and spatiotemporal information of the target station, and directly output the atmospheric precipitable water prediction results through grid interpolation; Step 6: Use accuracy evaluation metrics to evaluate the accuracy of the inversion model and verify its internal consistency, external accuracy, and real-time applicability.

[0041] Further, step 1 includes: Step 1-1: Extract geopotential, temperature, specific humidity, and relative humidity data for 37 pressure layers (from 1000 hPa to 1 hPa) from the ERA5 data. Calculate PWV using specific humidity q via vertical integration. In the formula, Take 1000 kg / m³, g Take 9.80665 m / s², q Indicates specific moisture content (g / kg). g Acceleration due to gravity (m / s²) 2 ), P The pressure is in hPa. P surfSurface air pressure, P top For the top atmospheric pressure, the integral is calculated numerically between each pressure layer using the trapezoidal rule. Steps 1-2: For each pressure layer, first calculate the water vapor pressure e using relative humidity RH and saturated water vapor pressure. P sa Calculate, and calculate again. q : Then, the atmospheric refractive index of each layer was calculated using the Smith-Wynter equation. N : In the formula, k 1 = 77.604 (K / hpa), k 2 = 64.79 (K / hpa), k 3 = 377600.0 (K2 / hpa), P d The partial pressure of dry air (hPa) T Let N be the thermodynamic temperature (K). The total zenith delay (ZTD) is calculated by integrating the refractive index layer by layer. ᵢ Δh represents the average refractive index between adjacent layers. ᵢ This represents the height difference between two adjacent floors.

[0042] Steps 1-3: Calculate the zenith static delay ZHD using the Saastamoinen model: In the formula, P s φ represents the surface air pressure (hPa) and the latitude of the station (rad). h The height of the ellipsoid is in km. Steps 1-4: Calculate the zenith wet delay (ZWD): Through the above process, ZWD and PWV time series data for the training period (2016-2019) were generated for all ERA5 grid points in the study area.

[0043] Further, step 2 includes: Step 2-1: Establish a periodic regression model for the direct conversion of ZWD to PWV, as shown in the following formula: In the formula, the units of PWV and ZWD are both mm, a0 is the intercept, and a1 is the principal linear transformation coefficient. doy This is the accumulated days of the year, with a value range of 1-365. A 1. B 1 is the coefficient for the annual periodic term. A 2. B 2 represents the coefficient for the semi-annual periodic term; The physical significance of this model lies in: through a 1. The dimensionless constant Π is implicitly expressed as a function of time, without the need for explicit calculation of Tm; Further, step 3 includes: Step 3-1: Divide the study area into grids at ERA5 resolution, totaling 2,225 grid points; Step 3-2: Determine the sliding window size as a 5×5 grid. The choice of window size needs to balance two factors: (1) the number of samples within the window is large enough to ensure the stability of parameter estimation; (2) the window should not be too large to maintain sensitivity to local atmospheric characteristics. The 5×5 window contains 25 grid points, each grid point has data at 35,064 time points, for a total of 875,000 sample points, which is sufficient to robustly solve for the 6 parameters; Step 3-3: Starting from the top left corner (northeast and west) of the study area, define the first window G1, covering the grid points in rows 1-5 and columns 1-5; Further, step 4 includes: Step 4-1: Aggregate all ZWD and PWV time series data of 25 grid points within window G1 to construct the design matrix X and observation vector Y; Step 4-2: Solve for the model parameters using the least squares method: Where, β=[ a 0, a 1, A 1, B 1, A 2, B 2]ᵀ; Step 4-3: Assign the solved parameters to the center grid point (3rd row, 3rd column) of window G1; Step 4-4: Move the window 5 grid points eastward, define window G2, covering the grid points of rows 1-5 and columns 6-10. Repeat steps 4-1 to 4-3; Steps 4-5: When the window reaches the easternmost point of the current latitude zone, the window moves 5 grid points south and returns to the westernmost point to begin processing the next latitude zone; Steps 4-6: Iterate until the entire study area is covered. Since the sliding step size is equal to the window size, the center points of adjacent windows do not overlap, ensuring that each grid point is uniquely assigned a set of parameters.

[0044] Ultimately, six coefficient grid diagrams are generated, each containing the same number of nodes as the original ERA5 grid, with each node storing a coefficient value.

[0045] Further, step 5 includes: Step 5-1: For the target GNSS station, process the raw GNSS observation data using PPP technology to obtain the ZTD in real time. PPP processing uses precise ephemeris and clock error products provided by IGS, with a calculation interval of 5 minutes. Step 5-2: Using the GPT3 model, predict the surface air pressure based on the latitude, longitude, elevation, and time information of the station. P s Then, the ZHD is calculated using the Saastamoinen model, and the ZWD is obtained from the ZTD. Step 5-3: Based on the station's latitude and longitude (λ, φ), obtain the model parameters for that location using bilinear interpolation in the six-coefficient grid map. For any coefficient C (e.g., ... a 0、 a (e.g., 1), assuming the station is located within a rectangle formed by four grid points (λ1, φ1), (λ2, φ1), (λ1, φ2), and (λ2, φ2), where λ1 < λ < λ2 and φ1 < φ < φ2, then the interpolation coefficients are: in, C 11 , C 21 , C 12 , C 22 These are the coefficient values ​​for the four grid points; Step 5-4: Calculate the cumulative day based on the date of the observation. doy ZWD, doy The six coefficients obtained by interpolation are substituted into the model formula to directly calculate the PWV value. Further, step 6 includes: Step 6-1: Validation using 2020 ERA5 data. For 2,225 grid points, input ERA5_ZWD into the SSW-DGM model, calculate the predicted PWV, and compare it with ERA5_PWV; like Figure 3As shown, the model exhibits excellent internal consistency in both the fitting period (2016-2019) and the prediction period (2020). The average RMSE during the fitting period was 0.24 mm, and the average RMSE during the prediction period was 0.21 mm. The average bias was 0.00 mm and -0.05 mm in the fitting and prediction periods, respectively, close to zero. Spatial distribution analysis shows that the RMSE in the southeast region (approximately 0.4-0.5 mm) is slightly higher than that in the northwest region. This is because the southeast is adjacent to the humid monsoon region, with a higher average PWV value, and the absolute error index naturally increases with the magnitude of the estimated variable. For the locations of 40 GNSS stations, such as... Figure 4 As shown, using ERA5_ZWD as input, the PWV inverted by SSW-DGM is compared with ERA5_PWV. The average RMSE is 0.230 mm, the average bias is -0.001 mm, and the maximum absolute residual is usually controlled within 1 mm. Step 6-2: Use ZTD products from 40 GNSS stations provided by CMONOC in 2020. ZHD is precisely calculated using ERA5 surface pressure, resulting in ZWD = ZTD_GNSS - ZHD_ERA5: like Figure 5 As shown, a spiral bar chart is used to display the error distribution of 40 GNSS stations. The PWV (PWV_DGM_GNSS) retrieved by SSW-DGM is compared with two reference values: (1) Traditional GNSS-PWV (PWV_GNSS_Tm): Using the standard conversion factor method, Tm is obtained from the numerical integration of the ERA5 profile, such as... Figure 6 As shown in the scatter plot of representative sites, the correlation coefficient between SSW-DGM inverted PWV and ERA5-PWV always exceeds 0.98, indicating a strong linear relationship. (2) Radiosonde PWV (PWV_RO): Obtained by integrating the profiles of 8 co-located radiosonde stations, such as... Figure 7 As shown, compared with the corresponding radiosonde PWV, the correlation coefficient also exceeds 0.98, proving that the model-driven conversion is highly synchronized with the actual atmospheric water vapor changes detected by physical sensors.

[0046] like Figure 8 and Figure 9 As shown, the time stability of the model was analyzed using monthly and seasonal error statistics. The RMSE was higher in summer (June-August), averaging about 1.7 mm, and decreased to about 0.7 mm in winter (December-February). This is because the water vapor content is high and highly variable in summer. However, the bias fluctuated around zero in all seasons, indicating that the seasonal harmonic term correctly modeled the annual periodicity and prevented systematic seasonal drift.

[0047] Step 6-3: Taking the BJFS site as an example, the model is trained using ERA5 data from 2020 to 2023, and the ZWD is calculated in real time using the original GNSS data from 2025 through PPP technology.

[0048] like Figure 10 As shown, the spatial relationship between the BJFS site and the surrounding SSW-DGM grid nodes is illustrated.

[0049] like Figure 11 As shown, the correlation coefficient between the PPP-based SSW-DGM inversion (PWV_DGM_PPP) and the high-precision radiosonde PWV is 0.978, comparable to 0.985 for the traditional method (PWV_GNSS_Tm). The residual histogram shows that 91.4% of the DGM inversion PWV values ​​fall within the acceptable threshold (92.0% for the traditional method). The residual distribution center is close to zero, verifying that even driven by real-time PPP estimation, the model does not introduce significant systematic bias.

[0050] Step 6-4: The formula for calculating the accuracy evaluation index is as follows: In the formula, x i Value, y i Measured value and These represent the average values ​​of their respective samples, and N represents the number of valid samples.

[0051] This embodiment achieves true meteorological parameter-free PWV inversion by pre-establishing a coefficient grid based on historical ERA5 data, eliminating the need for any real-time meteorological sensor data or external NWP data transmission during the application phase. It is particularly suitable for areas with scarce meteorological data, eliminating data latency bottlenecks and hardware costs.

[0052] Traditional grid models typically treat each grid point as an isolated estimation unit, making them susceptible to local outliers. Sliding windows, however, effectively smooth out local anomalies by aggregating neighboring data while preserving regional trends. More importantly, the window contains grid points at different elevations and latitudes, allowing least squares estimation to "learn" the local decline rate of transformation efficiency, implicitly considering the influence of topography and latitude on water vapor distribution.

[0053] like Figure 2As shown, this embodiment introduces a spatial sliding window algorithm, which takes into account the influence of elevation and latitude changes on the ZWD-PWV conversion efficiency of the model coefficients. This effectively captures the atmospheric spatial heterogeneity in complex terrain, avoids the lack of local adaptability of global grid models, and maintains stable accuracy in rugged terrain areas.

[0054] This embodiment effectively captures the seasonal time-varying characteristics of water vapor conversion efficiency by embedding annual and semi-annual harmonic terms into the regression model, thus making up for the shortcomings of traditional static linear models in effectively fitting the nonlinear characteristics of PWV affected by seasonal variations and atmospheric circulation cycle fluctuations.

[0055] Unlike "black box" machine learning models, the model in this embodiment has clear physical interpretability. Key coefficients a 1 acts as a dynamic proxy for the reciprocal of the weighted average temperature Tm, with the harmonic term capturing the seasonal variations in thermodynamics. The high correlation (>0.98) between the model and standard physical methods demonstrates that it has mathematically internalized the physical relationship between ZWD and PWV. In real-time scenario validation based on PPP, the model maintains comparable accuracy to traditional methods (correlation coefficient 0.978 vs 0.985), and 91.4% of the inversion results show consistent agreement with the reference values, proving the feasibility of the model in time-critical applications such as flash flood nowcasting and disaster early warning systems.

[0056] It should be noted that the calculation of atmospheric precipitable water volume (PWV) in this invention is not limited to the methods described above, and corresponding data products can also be selected. Similarly, the acquisition of GNSS ZTD in this invention is not limited to PPP technology; other GNSS data processing technologies such as network RTK and relative positioning can also be used.

[0057] In real-time PPP scenarios, although the ZWD derived from PPP includes orbital / clock error and ambiguity convergence noise, SSW-DGM effectively converts it into PWV without amplifying the error. Compared to traditional methods that require external meteorological data, SSW-DGM eliminates the data latency bottleneck (NWP data is typically delayed by several hours) and hardware costs (isotropic barometers), which is highly advantageous for time-critical applications such as flash flood nowcasting and disaster warning.

[0058] The above are merely preferred embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for direct inversion of GNSS precipitable water based on a spatial sliding window grid, characterized in that, include: Historical reanalysis data of the study area are obtained, and the zenith wet delay sequence and atmospheric precipitable water sequence are calculated based on the historical reanalysis data. Based on the zenith wet delay sequence and the atmospheric precipitable water sequence, a spatial sliding window direct grid model is constructed using a spatial sliding window strategy and a periodic regression model containing seasonal harmonic terms. The real-time GNSS zenith wet delay of the target station is obtained. Based on the real-time GNSS zenith wet delay and the spatial sliding window direct grid model, the atmospheric precipitable water prediction result of the target station is obtained by grid interpolation.

2. The method according to claim 1, characterized in that, The process of obtaining the zenith wet delay sequence and the atmospheric precipitable water sequence includes: Based on the specific humidity data and air pressure data in the historical reanalysis data, the atmospheric precipitable water sequence is obtained by vertical integration calculation. Based on the geopotential, temperature, specific humidity, and relative humidity data in the historical reanalysis data, the total zenith retardation sequence is obtained by integrating the atmospheric refractive index layer by layer, and the wet zenith retardation sequence is calculated based on the total zenith retardation sequence and the zenith static retardation sequence.

3. The method according to claim 1, characterized in that, The process of calculating the atmospheric precipitable water sequence includes: The calculation was performed using specific humidity data and air pressure data through a vertical integration formula, where the lower limit of integration is the surface air pressure and the upper limit is the top atmospheric pressure. The trapezoidal rule was used to perform numerical calculations between each pressure layer. The process of calculating the zenith wet delay sequence includes: The total zenith delay is obtained by integrating the atmospheric refractive index layer by layer. The zenith static delay is calculated using the Sastamonin model based on surface air pressure, station latitude, and ellipsoidal height. The wet zenith delay is obtained by subtracting the zenith static delay from the total zenith delay.

4. The method according to claim 1, characterized in that, In the periodic regression model used to construct the spatial sliding window direct grid model, atmospheric precipitable water is expressed as the sum of the following terms: The terms include the intercept term, the product of the zenith wet delay and the linear transformation coefficient, the annual cycle cosine term, the annual cycle sine term, the semi-annual cycle cosine term, and the semi-annual cycle sine term, where the annual cycle term and the semi-annual cycle term both use the accumulated days of the year as independent variables.

5. The method according to claim 1, characterized in that, The method further includes: performing gridding processing on the study area using the spatial sliding window strategy; The process of gridding the study area using the aforementioned spatial sliding window strategy includes: The study area is divided into grids based on the resolution of the historical reanalysis data; Determine the size and step size of the sliding window; Starting from the boundary of the study area, multiple sliding windows are defined sequentially, with each sliding window covering multiple grid points.

6. The method according to claim 1, characterized in that, The process of constructing the spatial sliding window direct grid model includes: For each sliding window, aggregate the zenith wet delay sequence and atmospheric precipitable water sequence data of all grid points within the window; The parameters of the periodic regression model are fitted using the least squares method to obtain the model parameters for each sliding window; Assign the model parameters of each sliding window to the center grid point of the corresponding sliding window; The model parameters of all sliding window center grid points are reorganized to generate multiple coefficient grid diagrams.

7. The method according to claim 6, characterized in that, The multiple coefficient grid plots include intercept grids, linear transformation coefficient grids, annual period cosine coefficient grids, annual period sine coefficient grids, semi-annual period cosine coefficient grids, and semi-annual period sine coefficient grids, with the spatial resolution of each grid plot being consistent with the resolution of the historical reanalysis data.

8. The method according to claim 1, characterized in that, The process of obtaining the real-time GNSS zenith wet delay of the target station includes: Obtain the real-time total zenith delay of the target station; Obtain the predicted surface air pressure value of the target station, and calculate the zenith static delay based on the predicted surface air pressure value; The real-time GNSS wet zenith delay is calculated based on the real-time total zenith delay and the zenith static delay.

9. The method according to claim 1, characterized in that, The process of obtaining the atmospheric precipitable water prediction result by grid interpolation based on the real-time GNSS zenith wet delay and the spatial sliding window direct grid model includes: Based on the latitude and longitude coordinates of the target station, the model parameters for the corresponding location are obtained by interpolation from multiple coefficient grid diagrams; Obtain the cumulative day of the year corresponding to the current observation time; The real-time GNSS zenith wet delay, the annual cumulative day, and the model parameters are substituted into the periodic regression model to calculate the predicted atmospheric precipitable water.

10. The method according to claim 1, characterized in that, The method also includes evaluating the accuracy of the spatial sliding window direct grid model; The process of evaluating the accuracy of the spatial sliding window direct grid model: Using historical reanalysis data from the test period, the zenith wet delay of the test period is input into the spatial sliding window direct grid model to obtain the predicted atmospheric precipitable water. The predicted atmospheric precipitable water is then compared with the atmospheric precipitable water of the test period, and the root mean square error and bias are calculated to verify internal consistency. Using GNSS observation zenith wet delay as input, the atmospheric precipitable water inverted by the spatial sliding window direct grid model is compared with the GNSS atmospheric precipitable water inverted by traditional methods and the radiosonde atmospheric precipitable water, to verify external accuracy. The zenith wet delay, calculated in real time using precise single-point positioning technology, is used as input to evaluate the inversion accuracy of the spatial sliding window direct grid model in real-time scenarios.

Images