Method for small area fusion of gnss and era5 tropospheric zenith delay

Abstract

The present application relates to the technical field of GNSS satellite navigation positioning and GNSS meteorology, in particular to a method for fusing GNSS and ERA5 tropospheric zenith delay in a small area, which comprises correcting the geopotential height of four adjacent ERA5 grid positions of a GNSS station to the ellipsoidal height of the GNSS station, and interpolating the ERA5 tropospheric zenith delay of the four adjacent grid points to the GNSS station based on the inverse distance weighting method. The present application corrects the systematic deviation of ERA5 ZTD, makes up for the deficiency of traditional interpolation methods which excessively rely on the distribution density of stations, supports the flexible construction of ZTD grid with different spatial resolutions to balance the calculation load and spatial details, effectively eliminates the systematic error and restores the smoothed local micro-meteorological details through the fine statistical modeling of the residual grid, and realizes the regional tropospheric delay modeling with high precision and high spatial resolution.

Description

Technical Field

[0001] This invention relates to the fields of GNSS satellite navigation and positioning and GNSS meteorology, specifically a method for small-area fusion of GNSS and ERA5 tropospheric zenith delay. Background Technology

[0002] Global Navigation Satellite System (GNSS) has become the core of modern spatiotemporal information infrastructure, playing a crucial role in high-precision applications such as Precise Point Positioning (PPP), network-based Real-Time Kinematics (RTK), and crustal deformation monitoring. However, the refraction delay of GNSS signals as they pass through the troposphere, known as tropospheric zenith delay (ZTD), is a major source of error limiting positioning accuracy, convergence speed, and the reliability of elevation calculations. Furthermore, as a tracer of atmospheric water vapor, high-precision ZTD products are also key parameters in GNSS meteorology for retrieving precipitable water vapor (PWV) and for assimilation in numerical weather prediction (NWM). Therefore, constructing high-precision, high-spatiotemporal-resolution, and uniformly covered regional tropospheric zenith delay models has become a common and cutting-edge research topic in the fields of GNSS navigation and positioning and meteorology.

[0003] Currently, the technical approaches for regional tropospheric zenith delay modeling can be mainly summarized into three categories: The first category is empirical models based on meteorological parameters, including classic Hopfield, Saastamoinen, and GPT series models. The bottleneck of this type of method lies in its limited accuracy and poor timeliness. It primarily fits the climatic mean of historical meteorological data, failing to capture drastic changes in short-term extreme weather. The centimeter-level model accuracy is insufficient to meet the high-precision positioning requirements at the millimeter level. The second category is the inversion of tropospheric zenith delay based on meteorological parameters provided by NWM reanalysis data, represented by the ERA5 meteorological product released by the ECMWF. While its advantage lies in providing globally seamless grid data, its drawback is the existence of significant systematic bias. Limited by the terrain smoothing effect and physical parameterization scheme of numerical models, there is often a fixed deviation of several millimeters or even centimeters between the ERA5 tropospheric zenith delay grid and the actual GNSS tropospheric zenith delay, which is more severe in complex terrain areas. The third category is mathematical interpolation models based on GNSS observation networks. This involves spatial interpolation using high-precision GNSS tropospheric zenith delay products calculated from reference stations. A drawback of this method is its extreme dependence on station distribution. While accuracy is acceptable in densely populated urban areas, interpolation errors increase dramatically with distance in sparsely populated or unevenly distributed mountainous regions and ocean edges, making it difficult to construct high-precision continuous tropospheric zenith delay grids.

[0004] Given the limitations of a single technical approach, fusing GNSS tropospheric zenith delay (TZD) and ERA5 tropospheric zenith delay (ERA5) offers a new research perspective for resolving the aforementioned contradictions. However, although the fusion strategy theoretically offers complementary advantages, it still faces extremely severe challenges in practical applications, mainly in the following aspects: First, existing fusion methods handle the spatiotemporal nonstationarity of systematic biases in a coarse manner. The deviation between ERA5 and GNSS tropospheric zenith delays is not a fixed constant but exhibits complex nonlinear variations with geographical location, elevation, and season. Existing simple linear regression or fixed deviation subtraction methods assume that the deviation is uniformly distributed within the region, which is seriously inconsistent with reality, resulting in large systematic errors remaining in local areas of the corrected tropospheric zenith delay grid. Second, the non-uniformity of the spatial distribution of GNSS stations severely restricts the modeling accuracy of the residual grid. The core of the fusion technology lies in modeling the difference between GNSS and ERA5 tropospheric zenith delays, which is usually referred to as the residual. However, since GNSS stations are often concentrated in plains or towns, the spatial distribution of the calculated residual sample points is extremely uneven. When traditional interpolation algorithms construct a residual grid covering the entire area based on these sparse and uneven sample points, they often exhibit oversmoothing or edge distortion, failing to accurately reproduce the spatial details of the residuals. Summary of the Invention

[0005] The purpose of this invention is to provide a method for small-area fusion of GNSS and ERA5 tropospheric zenith delay to solve the problems mentioned in the background art.

[0006] To address the aforementioned technical problems, this invention provides the following technical solution: a method for fusing GNSS and ERA5 tropospheric zenith delay in a small area, the method comprising: S1. Obtain observation data from GNSS stations in the specified time period, construct observation equations and perform parameter calculations, and extract the estimated tropospheric zenith delay at the station. S2. Obtain ERA5 reanalysis meteorological data within a specified time period. The meteorological data includes air temperature, air pressure, and specific humidity. Calculate the regional grid tropospheric zenith delay based on the integral method. S3. Correct the geopotential height of the four adjacent ERA5 grid locations of the GNSS station to the ellipsoidal height of the GNSS station, and interpolate the ERA5 tropospheric zenith delay of the four adjacent grid points to the GNSS station based on the inverse distance weighting method. S4. Based on the GNSS tropospheric zenith delay estimate obtained in step S1 and the ERA5 tropospheric zenith delay obtained in step S3, calculate the tropospheric zenith delay residual sequence of the station location. S5. Based on the accuracy and computational efficiency requirements of the actual application scenario, a grid with a specific spatial resolution is adaptively constructed within the target area: if it is used to capture small-scale atmospheric changes, a high-resolution dense grid is constructed; if it is used to explore large-area atmospheric changes, a low-resolution sparse grid is constructed. S6. Spatial interpolation of the ERA5 tropospheric zenith delay is performed based on the inverse distance weighting method to generate the ERA5 tropospheric zenith delay at the corresponding resolution of each target grid point in step S5. S7. Using the tropospheric zenith delay residual sequence obtained in step S4, generate the tropospheric zenith delay residual of the corresponding resolution grid in step S5 based on Ordinary Kriging interpolation. S8. The tropospheric zenith delay residual grid obtained in step S7 is superimposed and compensated onto the corresponding ERA5 tropospheric zenith delay grid obtained in step S6 to generate a high-precision regional tropospheric zenith delay fusion product.

[0007] This invention relates to a method for establishing a regional fusion GNSS and ERA5 tropospheric zenith delay model, the specific process of which is as follows: First, acquire observation data from regional GNSS CORS stations within a specified time period, construct observation equations and perform parameter calculations, and extract the ZTD estimate at the station. Second, acquire ERA5 reanalysis meteorological data within a specified time period, specifically including air temperature, air pressure, and specific humidity, and calculate the regional grid ZTD based on the integral method. Third, correct the geopotential height of the four adjacent ERA5 grid locations of the GNSS station to the geodetic height of the GNSS station, and interpolate the ERA5 ZTD of the four adjacent grid points to the GNSS station using the IDW method. Fourth, based on the GNSS ZTD obtained in the first step and the ERA5 ZTD at the GNSS station obtained in the third step... The fifth step involves adaptively constructing a grid with a specific spatial resolution within the target area, based on the accuracy and computational efficiency requirements of the actual application scenario: if the actual application scenario is to capture small-scale atmospheric changes, a high-resolution dense grid is set to improve local accuracy; if the actual application scenario is to explore large-area atmospheric changes, a low-resolution sparse grid is constructed to reduce computational load and improve processing efficiency. The sixth step involves spatial interpolation of the ERA5 ZTD based on the IDW method to generate the ERA5 ZTD with a specific spatial resolution at each target grid point in the fifth step. The seventh step involves using the ZTD residual sequence of the GNSS station locations obtained in the fourth step to generate the ZTD residual of the grid with a specific spatial resolution in the fifth step based on Ordinary Kriging interpolation. The eighth step involves superimposing and compensating the ZTD residual grid at each target grid point obtained in the seventh step onto the corresponding ERA5 ZTD obtained in the sixth step to generate a high-precision regional ZTD grid product.

[0008] Furthermore, the specific method for extracting the estimated tropospheric zenith delay at the station in S1 includes: Construct ionospherically-free combined observation equations using dual-frequency or multi-frequency GNSS observations to eliminate the first-order ionospheric delay in pseudorange and carrier measurements. By projecting the slant path delay to the zenith direction using the tropospheric mapping function, the corresponding tropospheric zenith delay estimate is obtained. Specifically, the combined observations of the dual-frequency ionosphere-free system in this invention can be expressed as follows: In the formula, P is the corresponding pseudorange observation in the dual-frequency non-ionospheric combined observation, and L is the phase observation; Indicates the corresponding observation satellite number, Indicates the corresponding receiver number; This indicates the combined observation frequency points of the dual-frequency non-ionospheric combination; This represents the geometric distance the signal travels from the observation satellite to the receiver. The speed at which light travels in a vacuum; This represents the clock bias of the receiver numbered r; This represents the clock bias of the observation satellite numbered 's'. The tropospheric slant path delay is the corresponding distance between the receiver numbered r and the observation satellite numbered s. For the observation frequency point The corresponding carrier wavelength; The combined hardware delay at the receiver end for pseudorange-free ionospheric combination; The combined hardware delay at the satellite end of the pseudorange-free ionosphere combination; Phase delay at the receiver end for phase-free ionosphere combination; The phase delay at the satellite end of the phase-free ionosphere combination; The combined ambiguity corresponding to the phase delay between the receiver and satellite ends in the phase-free ionospheric combination; This refers to the observational noise in the combined pseudorange and phase-free ionospheric observations. This refers to the multipath effect in pseudorange and phase observations without ionosphere; Among them, tropospheric oblique path delay It can be represented as: In the formula, ZHD is the tropospheric zenith static delay, ZWD is the tropospheric zenith wet delay, and AP is atmospheric pressure. E is a function of geodetic latitude B and elevation h; These represent the elevation angle and azimuth angle of the observed satellite, respectively. The horizontal gradient is in the north-south direction. The horizontal gradient is in the east-west direction; , Let E represent the mapping functions for the dry and wet components of the observed satellite at an elevation angle of E, respectively. This is the tropospheric horizontal gradient mapping function for observing the satellite's elevation angle E.

[0009] Furthermore, the specific method for calculating the regional grid tropospheric zenith delay based on the integral method in S2 includes: Using ERA5 to reanalyze multi-layer pressure, temperature and specific humidity data in meteorological data, the atmospheric refractive index is calculated layer by layer. By integrating the atmospheric refractive index calculation results of each layer in the vertical direction, the tropospheric zenith delay of each grid point is obtained. The number of pressure layers in the ERA5 reanalysis meteorological data corresponding to each grid point is the same, and the number of pressure layers in the ERA5 reanalysis meteorological data corresponding to each grid point is preset. Specifically, in this invention, each grid point of ERA5 contains 37 atmospheric pressure layers, and the atmospheric pressure, geopotential height, temperature, and specific humidity of each layer are needed to calculate the atmospheric refractive index of different layers: In the formula, N is the refractive index of the troposphere; , and These are preset physical constants related to the atmospheric refractive index; TP represents dry air pressure; TP represents temperature. It is the water vapor pressure.

[0010] The formula for calculating the tropospheric zenith delay using the integral method is as follows: In the formula, ZTD represents the tropospheric zenith delay of the corresponding grid point; m represents the total number of layers of reanalysis meteorological data contained above the GNSS station; and These are the heights of the vth and v+1th layers in the ERA5 reanalysis meteorological data corresponding to the grid points above the GNSS station; and These are the atmospheric refractive indices of the vth and v+1th layers in the ERA5 reanalysis meteorological data corresponding to the grid points above the GNSS station.

[0011] Furthermore, the specific implementation of S3 includes: The geopotential heights of the ERA5 grid locations adjacent to the GNSS station are successively converted into orthographic height and ellipsoidal height, thus unifying them with the GNSS station elevation system; Specifically, the ERA5 product provides geopotential data based on the distribution of pressure layers, while GNSS station data uses an ellipsoidal height elevation system. Therefore, the elevations of both need to be unified to the same standard. First, the ERA5 product needs to be converted to geopotential height. The relationship between geopotential and geopotential height is as follows: In the formula, G represents the potential, GH represents the potential height, and g is the acceleration constant. The obtained potential height is converted to positive height using the following formula: In the formula, Indicates the latitude of the corresponding grid point. Indicates latitude The positive height corresponding to the potential height GH of the corresponding grid point. Indicates latitude The radius of curvature of the Earth at that location. Indicates latitude The normal gravity value on the rotating ellipsoid at that location. This represents the normal gravity value on the rotating ellipsoid at latitude 45°.

[0012] The basic relationship for converting orthographic height to ellipsoidal height using geoid difference is as follows: In the formula, EH is the height of the ellipsoid after transformation, OH is the orthographic height before transformation, and NB is the geoid difference.

[0013] Based on the function model corresponding to the inverse distance weighting method, the ERA5 tropospheric zenith delay of the four grid points adjacent to the GNSS station is interpolated to the GNSS station location according to the weights of the grid points. The relevant function model expression is as follows: In the formula, n is the number of neighboring grid points, n=4; The distance from the i-th grid point adjacent to the GNSS station to the GNSS station; is the weight of the i-th grid point adjacent to the GNSS station; x is a preset grid point weight factor; and Let i represent the meteorological elements at the GNSS station and the meteorological elements at the i-th grid point, respectively, i∈[1,4].

[0014] Furthermore, in the process of adaptively constructing a grid with a specific spatial resolution within the target area in S5, the grid spacing of the high-resolution dense grid is no greater than 0.1°; and the grid spacing of the low-resolution sparse grid is no less than 0.5°.

[0015] Furthermore, the specific method for generating the tropospheric zenith delay residual of the corresponding resolution grid in step S5 based on Ordinary Kriging interpolation in step S7 includes: Based on the variation function of the residual sample points and interpolation points in the tropospheric zenith delay residual sequence obtained in step S4, an Ordinary Kriging function is constructed. By solving for the weight coefficients of the sample points involved in the interpolation, the tropospheric zenith delay residual value of the points to be interpolated in the corresponding resolution grid is obtained.

[0016] Specifically, the Ordinary Kriging function expression is as follows: In the formula, Points to be interpolated; The value is the tropospheric zenith delay residual for the point to be interpolated. For the q-th sample point in the tropospheric zenith delay residual sequence that participates in the interpolation, Let be the tropospheric zenith delay residual value at the q-th sample point in the tropospheric zenith delay residual sequence, and SN be the number of sample points in the tropospheric zenith delay residual sequence that participated in the interpolation. For sample points participating in interpolation The weighting coefficients are obtained by solving the following system of equations: in, The semivariogram values ​​are the w-th and q-th sample points involved in the interpolation of the tropospheric zenith delay residual sequence. The semivariogram values ​​are the w-th sample point participating in the interpolation and the grid points to be interpolated in the tropospheric zenith delay residual sequence. For a Lagrange daily number. For a spatial distance of... For any two points, the corresponding semi-variogram values ​​are calculated using an exponential function model, the specific expression of which is as follows: in, This is a scaling factor used to adjust the magnitude of the semivariogram value between two very close points; it is generally assigned a relatively small value, such as... ; For varying ranges, for small areas of atmosphere, it can make It equals the distance between the two farthest GNSS stations in the area, meaning that all GNSS stations in the small area are included in the modeling scope.

[0017] Furthermore, in the process of generating a high-precision regional tropospheric zenith delay fusion product in S8, the tropospheric zenith delay residual grid obtained by interpolation in step S7 is superimposed point by point with the corresponding ERA5 tropospheric zenith delay grid obtained in step S6; the output is a high-precision tropospheric zenith delay grid dataset with a unified spatial reference system.

[0018] Furthermore, the method is also applicable to real-time or near-real-time tropospheric zenith delay modeling, and supports dynamic updates of residual sequences and grid products.

[0019] Compared with existing technologies, the beneficial effects achieved by this invention are as follows: First, this invention calculates the regional GNSS station observations, extracts ZTD estimates, and combines the ZTD retrieved from ERA5 meteorological data to calculate the ZTD residual sequence at the GNSS station locations. Second, it constructs a custom resolution target grid adapted to different application scenarios within the study area, and uses Ordinary Kriging to construct the ZTD residual grid. Finally, the residual grid is superimposed and compensated onto the ERA5 ZTD grid to generate the final regional ZTD fusion product. This method corrects the systematic bias in ERA5 ZTD and compensates for the shortcomings of traditional interpolation methods that over-rely on station distribution density. It supports the flexible construction of ZTD grids with different spatial resolutions to balance computational load and spatial details. Through refined statistical modeling of the residual grid, it effectively eliminates systematic errors and restores smoothed local micro-meteorological details, achieving regional tropospheric delay modeling with both high accuracy and high spatial resolution. Attached Figure Description

[0020] The accompanying drawings are provided to further illustrate the invention and form part of the specification. They are used in conjunction with embodiments of the invention to explain the invention and do not constitute a limitation thereof. In the drawings: Figure 1 This is a flowchart illustrating the method for small-area fusion of GNSS and ERA5 tropospheric zenith delay according to the present invention. Figure 2 This is a schematic diagram of the spatial resolution of the adaptively constructed ZTD residual grid and the ERA5 ZTD grid within the target area in this embodiment of the invention; Figure 3 This is a comparison chart of the ZTD time series of the ADR2 station throughout 2023 in this embodiment of the invention. Detailed Implementation

[0021] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0022] Please see Figures 1-3 The present invention provides a technical solution: such as Figure 1 As shown, this embodiment provides a method for small-area fusion of GNSS and ERA5 tropospheric zenith delay, the method comprising: S1. Obtain observation data from GNSS stations in the specified time period, construct observation equations and perform parameter calculations, and extract the estimated tropospheric zenith delay at the station. The specific methods for extracting the estimated tropospheric zenith delay at the station in S1 include: Construct ionospherically-free combined observation equations using dual-frequency or multi-frequency GNSS observations to eliminate the first-order ionospheric delay in pseudorange and carrier measurements. By projecting the slant path delay to the zenith direction using the tropospheric mapping function, the corresponding tropospheric zenith delay estimate is obtained. Specifically, the combined observations of the dual-frequency ionosphere-free system in this invention can be expressed as follows: In the formula, P is the corresponding pseudorange observation in the dual-frequency non-ionospheric combined observation, and L is the phase observation; Indicates the corresponding observation satellite number, Indicates the corresponding receiver number; This indicates the combined observation frequency points of the dual-frequency non-ionospheric combination; This represents the geometric distance the signal travels from the observation satellite to the receiver. The speed at which light travels in a vacuum; This represents the clock bias of the receiver numbered r; This represents the clock bias of the observation satellite numbered 's'. The tropospheric slant path delay is the corresponding distance between the receiver numbered r and the observation satellite numbered s. For the observation frequency point The corresponding carrier wavelength; The combined hardware delay at the receiver end for pseudorange-free ionospheric combination; The combined hardware delay at the satellite end of the pseudorange-free ionosphere combination; Phase delay at the receiver end for phase-free ionosphere combination; The phase delay at the satellite end of the phase-free ionosphere combination; The combined ambiguity corresponding to the phase delay between the receiver and satellite ends in the phase-free ionospheric combination; This refers to the observational noise in the combined pseudorange and phase-free ionospheric observations. This refers to the multipath effect in pseudorange and phase observations without ionosphere; Among them, tropospheric oblique path delay It can be represented as: In the formula, ZHD is the tropospheric zenith static delay, ZWD is the tropospheric zenith wet delay, and AP is atmospheric pressure. E is a function of geodetic latitude B and elevation h; These represent the elevation angle and azimuth angle of the observed satellite, respectively. The horizontal gradient is in the north-south direction. The horizontal gradient is in the east-west direction; , Let E represent the mapping functions for the dry and wet components of the observed satellite at an elevation angle of E, respectively. This is the tropospheric horizontal gradient mapping function for observing the satellite's elevation angle E.

[0023] S2. Obtain ERA5 reanalysis meteorological data within a specified time period. The meteorological data includes air temperature, air pressure, and specific humidity. Calculate the regional grid tropospheric zenith delay based on the integral method. The specific methods for calculating the tropospheric zenith delay of the regional grid based on the integral method in S2 include: Using ERA5 to reanalyze multi-layer pressure, temperature and specific humidity data in meteorological data, the atmospheric refractive index is calculated layer by layer. By integrating the atmospheric refractive index calculation results of each layer in the vertical direction, the tropospheric zenith delay of each grid point is obtained. The number of pressure layers in the ERA5 reanalysis meteorological data corresponding to each grid point is the same, and the number of pressure layers in the ERA5 reanalysis meteorological data corresponding to each grid point is preset. Specifically, in this invention, each grid point of ERA5 contains 37 atmospheric pressure layers, and the atmospheric pressure, geopotential height, temperature, and specific humidity of each layer are needed to calculate the atmospheric refractive index of different layers: In the formula, N is the refractive index of the troposphere; , and These are preset physical constants related to the atmospheric refractive index; TP represents dry air pressure; TP represents temperature. It is the water vapor pressure.

[0024] The formula for calculating the tropospheric zenith delay using the integral method is as follows: In the formula, ZTD represents the tropospheric zenith delay of the corresponding grid point; m represents the total number of layers of reanalysis meteorological data contained above the GNSS station; and These are the heights of the vth and v+1th layers in the ERA5 reanalysis meteorological data corresponding to the grid points above the GNSS station; and These are the atmospheric refractive indices of the vth and v+1th layers in the ERA5 reanalysis meteorological data corresponding to the grid points above the GNSS station.

[0025] S3. Correct the geopotential height of the four adjacent ERA5 grid locations of the GNSS station to the ellipsoidal height of the GNSS station, and interpolate the ERA5 tropospheric zenith delay of the four adjacent grid points to the GNSS station based on the inverse distance weighting method. The specific implementation methods of S3 include: The geopotential heights of the ERA5 grid locations adjacent to the GNSS station are successively converted into orthographic height and ellipsoidal height, thus unifying them with the GNSS station elevation system; Specifically, the ERA5 product provides geopotential data based on the distribution of pressure layers, while GNSS station data uses an ellipsoidal height elevation system. Therefore, the elevations of both need to be unified to the same standard. First, the ERA5 product needs to be converted to geopotential height. The relationship between geopotential and geopotential height is as follows: In the formula, G represents the potential, GH represents the potential height, and g is the acceleration constant. The obtained potential height is converted to positive height using the following formula: In the formula, Indicates the latitude of the corresponding grid point. Indicates latitude The positive height corresponding to the potential height GH of the corresponding grid point. Indicates latitude The radius of curvature of the Earth at that location. Indicates latitude The normal gravity value on the rotating ellipsoid at that location. This represents the normal gravity value on the rotating ellipsoid at latitude 45°.

[0026] The basic relationship for converting orthographic height to ellipsoidal height using geoid difference is as follows: In the formula, EH is the height of the ellipsoid after transformation, OH is the orthographic height before transformation, and NB is the geoid difference.

[0027] Based on the function model corresponding to the inverse distance weighting method, the ERA5 tropospheric zenith delay of the four grid points adjacent to the GNSS station is interpolated to the GNSS station location according to the weights of the grid points. The relevant function model expression is as follows: In the formula, n is the number of neighboring grid points, n=4; The distance from the i-th grid point adjacent to the GNSS station to the GNSS station; is the weight of the i-th grid point adjacent to the GNSS station; x is a preset grid point weight factor; and Let i represent the meteorological elements at the GNSS station and the meteorological elements at the i-th grid point, respectively, i∈[1,4].

[0028] S4. Based on the GNSS tropospheric zenith delay estimate obtained in step S1 and the ERA5 tropospheric zenith delay obtained in step S3, calculate the tropospheric zenith delay residual sequence of the station location. S5. Based on the accuracy and computational efficiency requirements of the actual application scenario, adaptively construct a grid with a specific spatial resolution within the target area: if used to capture small-scale atmospheric changes, construct a high-resolution dense grid; if used to explore large-area atmospheric changes, construct a low-resolution sparse grid; such as... Figure 2 The diagram shows the spatial resolution of the adaptively constructed ZTD residual grid and the ERA5 ZTD grid within the target area. In step S5, during the adaptive construction of a grid with a specific spatial resolution within the target area, the grid spacing of the high-resolution dense grid is no greater than 0.1°; and the grid spacing of the low-resolution sparse grid is no less than 0.5°.

[0029] S6. Spatial interpolation of the ERA5 tropospheric zenith delay is performed based on the inverse distance weighting method to generate the ERA5 tropospheric zenith delay at the corresponding resolution of each target grid point in step S5. S7. Using the tropospheric zenith delay residual sequence obtained in step S4, generate the tropospheric zenith delay residual of the corresponding resolution grid in step S5 based on Ordinary Kriging interpolation. The specific methods for generating the tropospheric zenith delay residual based on Ordinary Kriging interpolation in step S5 in S7 include: Based on the variation function of the residual sample points and interpolation points in the tropospheric zenith delay residual sequence obtained in step S4, an Ordinary Kriging function is constructed. By solving for the weight coefficients of the sample points involved in the interpolation, the tropospheric zenith delay residual value of the points to be interpolated in the corresponding resolution grid is obtained.

[0030] Specifically, the Ordinary Kriging function expression is as follows: In the formula, Points to be interpolated; The value is the tropospheric zenith delay residual for the point to be interpolated. This refers to the q-th sample point in the tropospheric zenith delay residual sequence that participates in the interpolation. SN is the tropospheric zenith delay residual value of the qth sample point in the tropospheric zenith delay residual sequence, and SN is the number of sample points in the tropospheric zenith delay residual sequence that participate in interpolation. For sample points participating in interpolation The weighting coefficients are obtained by solving the following system of equations: in, The semivariogram values ​​are the w-th and q-th sample points involved in the interpolation of the tropospheric zenith delay residual sequence. The semivariogram values ​​are the w-th sample point participating in the interpolation and the grid points to be interpolated in the tropospheric zenith delay residual sequence. For a Lagrange daily number. For a spatial distance of... For any two points, the corresponding semi-variogram values ​​are calculated using an exponential function model, the specific expression of which is as follows: in, This is a scaling factor used to adjust the magnitude of the semivariogram value between two very close points; it is generally assigned a relatively small value, such as... ; For varying ranges, for small areas of atmosphere, it can make It equals the distance between the two farthest GNSS stations in the area, meaning that all GNSS stations in the small area are included in the modeling scope.

[0031] S8. The tropospheric zenith delay residual grid obtained in step S7 is superimposed and compensated onto the corresponding ERA5 tropospheric zenith delay grid obtained in step S6 to generate a high-precision regional tropospheric zenith delay fusion product.

[0032] In the process of generating a high-precision regional tropospheric zenith delay fusion product in step S8, the tropospheric zenith delay residual grid obtained by interpolation in step S7 is superimposed point by point with the corresponding ERA5 tropospheric zenith delay grid obtained in step S6; the output is a high-precision tropospheric zenith delay grid dataset with a unified spatial reference system.

[0033] The method described in this embodiment is also applicable to real-time or near-real-time tropospheric zenith delay modeling, and supports dynamic updates of residual sequences and grid products.

[0034] This embodiment selects GNSS ZTD products from 36 GNSS stations (10 modeling stations and 26 verification stations) in the Netherlands (3.5°E-8.0°E, 50.5°N-53.5°N) in 2023 and ERA5 grid ZTD products from the same period for modeling and verification. Figure 3As shown. Its characteristic is that it includes the following steps: S1. Acquire observation data from GNSS stations in the specified time period, construct observation equations and perform parameter calculations, and extract the estimated tropospheric zenith delay (ZTD) value at the station. For example, for the AMST station (52.39°N, 4.84°E), at 12:00 UTC on May 1, 2023, the estimated GNSS ZTD value extracted by PPP calculation is 2417.6 mm. S2. Obtain ERA5 reanalysis meteorological data within a specified time period. The meteorological data includes air temperature, air pressure, and specific humidity. Calculate the regional grid tropospheric zenith delay based on the integral method. S3. Correct the geopotential heights of the four adjacent ERA5 grid locations of the GNSS station to the ellipsoidal height of the GNSS station, and interpolate the ERA5 tropospheric zenith delay of the four adjacent grid points to the GNSS station based on the inverse distance weighting method. For example, for the AMST station, at 12:00 UTC on 2023-05-01, the ERA5 ZTD value is 2411.2 mm; S4. Based on the GNSS tropospheric zenith delay estimate obtained in step S1 and the ERA5 tropospheric zenith delay obtained in step S3, calculate the tropospheric zenith delay residual sequence for the station location. For the AMST station mentioned above, the ZTD residual calculated at 12:00 UTC on 2023-05-01 is 6.4 mm, i.e., 2417.6 mm - 2411.2 mm; S5. Based on the accuracy and computational efficiency requirements of the actual application scenario, a grid with a specific spatial resolution is adaptively constructed within the target area: If used to capture small-scale atmospheric changes, a high-resolution dense grid is constructed to improve local accuracy. In this embodiment, the target grid resolution is set to 0.1°×0.1°, which is higher than the 0.25°×0.25° of the ERA5 product, providing denser residual compensation information and improving the interpolation accuracy of the user end; if used to explore large-area atmospheric changes, a low-resolution sparse grid of 0.5°×0.5° is constructed to reduce computational load and improve processing efficiency. S6. Spatial interpolation of the ERA5 tropospheric zenith delay is performed based on the inverse distance weighting method to generate the ERA5 tropospheric zenith delay at the corresponding resolution of each target grid point in step S5. S7. Using the tropospheric zenith delay residual sequence obtained in step S4, generate the tropospheric zenith delay residual of the corresponding resolution grid in step S5 based on Ordinary Kriging interpolation. S8. The tropospheric zenith delay residual grid obtained in step S7 is superimposed and compensated onto the corresponding ERA5 tropospheric zenith delay grid obtained in step S6 to generate a high-precision regional tropospheric zenith delay fusion product. Taking the verification station ADR2 in the target area as an example, at 12:00 UTC on May 1, 2023, the ERA5 ZTD value at this location was 2401.1 mm, while the measured GNSS ZTD true value was 2406.1 mm, resulting in a systematic bias of 5.0 mm. The residual correction obtained by the Ordinary Kriging interpolation method was 3.6 mm. After superimposing this residual onto the ERA5 ZTD, the fused ZTD value obtained was 2404.7 mm. The results show that, compared with the ERA5 ZTD, the deviation between the fused ZTD value after residual correction and the true value was reduced from 5.0 mm to 1.4 mm, effectively correcting the inherent systematic bias of ERA5 under normal weather conditions, and verifying the high accuracy and reliability of the regional ZTD model.

[0035] In the embodiments provided by this invention, the method for establishing a regional fusion GNSS and ERA5 tropospheric zenith delay model comprises the following steps: First, acquiring observation data from regional GNSS CORS stations within a specified time period, constructing observation equations and performing parameter calculations, and extracting ZTD estimates at the stations; Second, acquiring ERA5 reanalysis meteorological data within a specified time period, specifically including air temperature, air pressure, and specific humidity, and calculating the regional grid ZTD based on the integral method; Third, correcting the geopotential height of the four adjacent ERA5 grid locations of the GNSS station to the geodetic height of the GNSS station, and interpolating the ERA5 ZTD of the four adjacent grid points of the GNSS station to the GNSS station based on the IDW method; Fourth, based on the GNSS ZTD obtained in the first step and the ERA5 values ​​at the GNSS station obtained in the third step... The fifth step involves adaptively constructing a grid with a specific spatial resolution within the target area, based on the accuracy and computational efficiency requirements of the actual application scenario: if the actual application scenario is to capture small-scale atmospheric changes, a high-resolution dense grid is set to improve local accuracy; if the actual application scenario is to explore large-area atmospheric changes, a low-resolution sparse grid is constructed to reduce computational load and improve processing efficiency. The sixth step involves spatial interpolation of the ERA5 ZTD based on the IDW method to generate the ERA5 ZTD with a specific spatial resolution at each target grid point in the fifth step. The seventh step involves using the ZTD residual sequence of the GNSS station location obtained in the fourth step to generate the ZTD residual of the grid with a specific spatial resolution in the fifth step based on OrdinaryKriging interpolation. The eighth step involves superimposing and compensating the ZTD residual grid at each target grid point obtained in the seventh step onto the corresponding ERA5 ZTD obtained in the sixth step to generate a high-precision regional ZTD grid product. This method corrects the systematic bias in ERA5 ZTD and makes up for the shortcomings of traditional interpolation methods that rely too much on station distribution density. It supports the flexible construction of ZTD grids with different spatial resolutions to balance computational load and spatial details. Through refined statistical modeling of the residual grid, it effectively eliminates systematic errors and restores smoothed local micro-meteorological details, achieving regional tropospheric delay modeling with both high accuracy and high spatial resolution.

[0036] It should be noted that, in this document, relational terms such as "first" and "second" are used only to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such process, method, article, or apparatus.

[0037] Finally, it should be noted that the above descriptions are merely preferred embodiments of the present invention and are not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for small-area fusion of GNSS and ERA5 tropospheric zenith delay, characterized in that, The method includes: S1. Obtain observation data from GNSS stations in the specified time period, construct observation equations and perform parameter calculations, and extract the estimated tropospheric zenith delay at the station. S2. Obtain ERA5 reanalysis meteorological data within a specified time period. The meteorological data includes air temperature, air pressure, and specific humidity. Calculate the regional grid tropospheric zenith delay based on the integral method. S3. Correct the geopotential height of the four adjacent ERA5 grid locations of the GNSS station to the ellipsoidal height of the GNSS station, and interpolate the ERA5 tropospheric zenith delay of the four adjacent grid points to the GNSS station based on the inverse distance weighting method. S4. Based on the GNSS tropospheric zenith delay estimate obtained in step S1 and the ERA5 tropospheric zenith delay obtained in step S3, calculate the tropospheric zenith delay residual sequence of the station location. S5. Based on the accuracy and computational efficiency requirements of the actual application scenario, a grid with a specific spatial resolution is adaptively constructed within the target area: if it is used to capture small-scale atmospheric changes, a high-resolution dense grid is constructed; if it is used to explore large-area atmospheric changes, a low-resolution sparse grid is constructed. S6. Spatial interpolation of the ERA5 tropospheric zenith delay is performed based on the inverse distance weighting method to generate the ERA5 tropospheric zenith delay at the corresponding resolution of each target grid point in step S5. S7. Using the tropospheric zenith delay residual sequence obtained in step S4, generate the tropospheric zenith delay residual of the corresponding resolution grid in step S5 based on Ordinary Kriging interpolation. S8. The tropospheric zenith delay residual grid obtained in step S7 is superimposed and compensated onto the corresponding ERA5 tropospheric zenith delay grid obtained in step S6 to generate a high-precision regional tropospheric zenith delay fusion product.

2. The method for small-area fusion of GNSS and ERA5 tropospheric zenith delay according to claim 1, characterized in that: The specific methods for extracting the estimated tropospheric zenith delay at the station in S1 include: Construct ionospherically-free combined observation equations using dual-frequency or multi-frequency GNSS observations to eliminate the first-order ionospheric delay in pseudorange and carrier measurements. By projecting the slant path delay to the zenith direction using the tropospheric mapping function, the corresponding tropospheric zenith delay estimate is obtained.

3. The method for small-area fusion of GNSS and ERA5 tropospheric zenith delay according to claim 1, characterized in that: The specific methods for calculating the tropospheric zenith delay of the regional grid based on the integral method in S2 include: Using ERA5 to reanalyze multi-layer pressure, temperature and specific humidity data in meteorological data, the atmospheric refractive index is calculated layer by layer. By integrating the atmospheric refractive index calculation results of each layer in the vertical direction, the tropospheric zenith delay of each grid point is obtained. The number of pressure layers in the ERA5 reanalysis meteorological data corresponding to each grid point is the same, and the number of pressure layers in the ERA5 reanalysis meteorological data corresponding to each grid point is preset.

4. The method for small-area fusion of GNSS and ERA5 tropospheric zenith delay according to claim 1, characterized in that: The specific implementation methods of S3 include: The geopotential heights of the ERA5 grid locations adjacent to the GNSS station are successively converted into orthographic height and ellipsoidal height, thus unifying them with the GNSS station elevation system; Based on the function model corresponding to the inverse distance weighting method, the ERA5 tropospheric zenith delay of the four grid points adjacent to the GNSS station is interpolated to the GNSS station location according to the weights of the grid points. The relevant function model expression is as follows: In the formula, n is the number of neighboring grid points, n=4; The distance from the i-th grid point adjacent to the GNSS station to the GNSS station; is the weight of the i-th grid point adjacent to the GNSS station; x is a preset grid point weight factor; and Let i represent the meteorological elements at the GNSS station and the meteorological elements at the i-th grid point, respectively, i∈[1,4].

5. The method for small-area fusion of GNSS and ERA5 tropospheric zenith delay according to claim 1, characterized in that: In step S5, during the adaptive construction of a grid with a specific spatial resolution within the target area, the grid spacing of the high-resolution dense grid is no greater than 0.1°; and the grid spacing of the low-resolution sparse grid is no less than 0.5°.

6. The method for small-area fusion of GNSS and ERA5 tropospheric zenith delay according to claim 1, characterized in that: The specific methods for generating the tropospheric zenith delay residual based on Ordinary Kriging interpolation in step S5 in S7 include: Based on the variation function of the residual sample points and interpolation points in the tropospheric zenith delay residual sequence obtained in step S4, an Ordinary Kriging function is constructed. By solving for the weight coefficients of the sample points involved in the interpolation, the tropospheric zenith delay residual value of the points to be interpolated in the corresponding resolution grid is obtained.

7. The method for small-area fusion of GNSS and ERA5 tropospheric zenith delay according to claim 1, characterized in that: In the process of generating a high-precision regional tropospheric zenith delay fusion product in step S8, the tropospheric zenith delay residual grid obtained by interpolation in step S7 is superimposed point by point with the corresponding ERA5 tropospheric zenith delay grid obtained in step S6; the output is a high-precision tropospheric zenith delay grid dataset with a unified spatial reference system.

8. The method for small-area fusion of GNSS and ERA5 tropospheric zenith delay according to any one of claims 1-7, characterized in that: The method is also applicable to real-time or near-real-time tropospheric zenith delay modeling, and supports dynamic updates of residual sequences and grid products.

Images