A method for calculating regional zenith tropospheric delay

Abstract

The present invention belongs to the technical field of space geodesy, and particularly relates to a method for calculating the zenith tropospheric delay in a region. The method first selects GNSS stations that meet the accuracy requirements for calculating the zenith tropospheric delay in the target region according to the corresponding relationship between the calculation accuracy of the zenith tropospheric delay and the shortest distance of the stations; wherein, the corresponding relationship is obtained by performing simulation calculations on the target region; then, according to the observables of the selected GNSS stations, the zenith tropospheric delay values at each GNSS station are calculated; the zenith tropospheric delays at each GNSS station are extended to other reference surfaces, and after the extension, grid processing is performed to obtain the regional delay values on other reference surfaces; finally, the regional delay values on other reference surfaces are inversely extended to the ground surface to obtain the regional zenith tropospheric delay values. The present invention can obtain accurate calculation results of the zenith tropospheric delay in the target region based on a small number of GNSS stations.

Description

Technical Field

[0001] The present invention belongs to the technical field of space geodesy, and particularly relates to a method for calculating the zenith tropospheric delay in a region. Background Art

[0002] The troposphere causes propagation delays to space observation technologies that mainly use electromagnetic waves as detection means, such as global navigation satellite systems, very long baseline interferometry, satellite laser ranging, synthetic aperture radar interferometry, and other technologies. This kind of delay is affected by meteorological conditions such as temperature, pressure, and water vapor content, has strong time variability, varies with geographical location and propagation path, and cannot be eliminated by using multi-frequency combined observations, which has become one of the main bottlenecks restricting the accuracy of space geodesy technologies.

[0003] For the convenience of research, the tropospheric delay is usually divided into two parts: the dry component and the wet component. Among them, the dry component has strong regularity and has a strong correlation with the pressure, elevation, latitude, etc. of the area to be studied, and can be accurately calculated using a simple empirical model; the wet delay is mainly affected by elements such as water vapor and temperature in the profile above the area to be studied, and the changes of these profile elements are unpredictable, resulting in much more difficulty in accurately determining the wet delay than the dry delay. Therefore, the restriction of the troposphere on space geodesy technologies is mainly reflected in the wet delay, and the acquisition of high-precision tropospheric delay information will also significantly improve the service capabilities of space technologies such as GNSS.

[0004] Currently, there are usually the following methods for calculating tropospheric delay: the stratified integration method, the empirical model method, the estimation method, and the machine learning method. Among them, the stratified integration method gives a rigorous calculation formula for ZWD from a physical perspective, with the highest theoretical accuracy, and can be used to evaluate the accuracy of other methods. With the development of technologies such as water vapor radiometers, radiosondes, drones, and artificial satellites, the means of obtaining meteorological data are more diverse, providing rich data sources for the integration method. However, the cost of obtaining profile meteorological data is still relatively high. Although there are several authoritative international institutions, such as ECMWF and NCEP, which can regularly release spatial stratified grid meteorological data, the product production cycle is long, and the updates often lag by several days or even months, which is not conducive to the real-time calculation of tropospheric delay. The empirical model method has the lowest calculation cost and the highest calculation efficiency. Compared with the integration method, the empirical model method greatly simplifies the solution process. Especially for non-meteorological parameter empirical models, such models are often constructed based on multi-year meteorological data or authoritative meteorological products and can be pre-implanted into data processing software such as GNSS and VLBI, with strong universality. However, the empirical model method sacrifices the calculation accuracy of ZWD and cannot reflect the characteristics of ZWD under sudden and abnormal weather conditions. The tropospheric delay calculated based on the estimation method has the highest time resolution and can obtain a high-precision delay sequence throughout the year without interruption. However, the prerequisite for the estimation method to obtain accurate delay is that there are uniform and sufficient spatial geodetic stations such as GNSS in the study area. Moreover, to separate the coupling effect of position error and tropospheric delay, it is also necessary to perform static continuous observations for a sufficient long time. This makes it difficult to reflect the application effect of this method in areas with blank stations or areas where it is inconvenient to build stations. The machine learning method has great potential. Many scholars have confirmed through a large number of experiments that machine learning has incomparable advantages in complex non-linear approximation fields such as meteorological and hydrological analysis and forecasting. However, the machine learning method usually requires a large amount of historical data, and it is also difficult to reflect the effect of this method in areas with scarce data and insufficient observations. Summary of the Invention

[0005] The purpose of the present invention is to provide a method for calculating the regional zenith tropospheric delay to solve the problem that the existing technical methods cannot achieve high-precision tropospheric delay calculation results in areas with limited station resources.

[0006] To solve the above technical problems, the present invention provides a method for calculating the regional zenith tropospheric delay, including the following steps:

[0007] 1) Select GNSS stations that meet the requirements for the calculation accuracy of the regional zenith tropospheric delay according to the corresponding relationship between the calculation accuracy of the zenith tropospheric delay and the shortest distance of the stations; wherein, the corresponding relationship is obtained through simulation calculation of the target area.

[0008] 2) Calculate the zenith tropospheric delay values at each GNSS station based on the observed values of the selected GNSS stations.

[0009] 3) Extend the zenith tropospheric delays at each GNSS station to other reference surfaces, and after extension, perform gridding processing to obtain the regional delay values on other reference surfaces.

[0010] 4) Reverse-extend the regional delay values on other reference surfaces to the ground surface to obtain the regional zenith tropospheric delay values.

[0011] The beneficial effects are as follows: The method of the present invention first needs to establish the correspondence between the calculation accuracy of the zenith tropospheric delay and the shortest distance of the stations. Then, this correspondence can be used to select appropriate GNSS stations in the target area. On this basis, combined with the tropospheric delay extension, reverse extension models, and gridding methods, the real-time and high-precision calculation of the regional tropospheric delay is finally realized, successfully solving the current situation of difficult or inaccurate acquisition of regional tropospheric delay information due to limited station resources, and can provide effective technical support for the research and application in the field of space geodesy such as GNSS precise point positioning, very long baseline interferometry, and spaceborne radar interferometry.

[0012] Further, the process of the simulation calculation includes: calculating the zenith tropospheric delay values of the grids in the area based on the vertical profile meteorological data provided by the hierarchical grid meteorological model as the reference values of the zenith tropospheric delay; using the data of the GNSS stations in the target area and the barometric data at the corresponding positions to calculate the regional zenith tropospheric delay values, and comparing them with the reference values of the zenith tropospheric delay to obtain the errors of each grid point; performing regression analysis on the errors of each grid point and the shortest distance between each grid point and the station, so as to obtain the correspondence between the calculation accuracy of the zenith tropospheric delay and the shortest distance of the station.

[0013] The beneficial effects are as follows: By using the above method, the correspondence between the calculation accuracy of the zenith tropospheric delay and the shortest distance of the station can be accurately and quickly obtained, and then this relationship can be used to accurately evaluate whether the GNSS stations in the area meet the accuracy requirements.

[0014] Further, it is necessary to combine the layered integration method to calculate the zenith tropospheric delay values of the grids in the area based on the vertical profile meteorological data provided by the hierarchical grid meteorological model. The calculation formula of the layered integration method is:

[0015]

[0016] In the formula, ZTD is the calculated zenith tropospheric delay value; surf and trop are the surface elevation and tropopause elevation respectively; k1, k'2, k3 are constant coefficients; p, p w, \(P\), \(e\), and \(T\) are the air pressure, water vapor pressure, and air temperature at each height level respectively.

[0017] Its beneficial effect is that the stratified integration method gives a rigorous calculation formula from a physical perspective and has the highest theoretical accuracy. Therefore, it is used to evaluate the accuracy of other methods.

[0018] Furthermore, when the existing GNSS stations in the target area cannot meet the accuracy requirements for calculating the zenith tropospheric delay in the target area, temporary GNSS stations need to be set up in the target area to meet the requirements.

[0019] Its beneficial effect is that the accuracy requirements for calculating the zenith tropospheric delay in the target area can be met by setting up temporary GNSS stations.

[0020] Furthermore, the other reference surface is the geoid.

[0021] Its beneficial effect is that using the geoid as the other reference surface can improve the accuracy of calculating the tropospheric delay in the region.

[0022] Furthermore, in step 3), a quadratic surface function is used for gridding.

[0023] Furthermore, the form of the quadratic surface function is:

[0024] f(B,L) = a0 + a1B + a2L + a3B 2 + a4L 2 + a5BL

[0025] In the formula, f(B,L) is the quadratic surface function with respect to the geodetic latitude B and the geodetic longitude L; a i is the function coefficient, i = 1, 2, …, 5.

[0026] Furthermore, the zenith tropospheric delay values at each GNSS station in step 2) include the zenith tropospheric dry delay value and the zenith tropospheric wet delay value at each GNSS station, which are respectively:

[0027]

[0028] ZWD surf = ZTD surf - ZHD surf

[0029] In the formula, ZHD surf , ZWD surf are respectively the zenith tropospheric dry and wet delay values at the GNSS station; P surf is the surface air pressure value; B and H are respectively the geodetic latitude and the geodetic height at the GNSS station; ZTD surf are respectively the zenith tropospheric delay values at the GNSS station.

[0030] Its beneficial effects are as follows: By combining the surface air pressure value and the surface GNSS station information, the regional zenith tropospheric delay value can be accurately calculated.

[0031] Further, the formula for the extension in step 3) is:

[0032]

[0033]

[0034] In the formula, ZHD g , ZWD g are the dry and wet zenith tropospheric delay values after being extended to the geoid respectively; P surf is the surface air pressure value; B is the geodetic latitude at the GNSS station; H norm is the orthometric height to be extended; N is the geoid undulation at the position to be extended; ZWD surf is the wet zenith tropospheric delay value at the surface.

[0035] Further, the formula for the reverse extension in step 4) is:

[0036] ZHD surf = ZHD g ×(1 - 0.0000226×H norm ) 5.225

[0037]

[0038] In the formula, ZHD surf , ZWD surf are the dry and wet zenith tropospheric delay values at the surface respectively. Description of the Drawings

[0039] Figure 1 is the flowchart of the method for calculating the regional zenith tropospheric delay of the present invention;

[0040] Figure 2 is the distribution map of GNSS stations in the region in the embodiment of the method of the present invention;

[0041] Figure 3 is the geoid map of the region in the embodiment of the method of the present invention;

[0042] Figure 4 is the pre - calculation accuracy map of the regional zenith tropospheric delay in the embodiment of the method of the present invention;

[0043] Figure 5 is the distribution map of the shortest distance between grid points and stations in the embodiment of the method of the present invention. Detailed Embodiment

[0044] In order to make the objectives, technical solutions and advantages of the present invention clearer, the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0045] Method Embodiment:

[0046] As Figure 1 shown, the overall implementation process of the method of the present invention is as follows:

[0047] Step 1: Collect data in the target area, including the three-dimensional position information of GNSS stations, geoid and surface elevation information in the target area.

[0048] Collect available GNSS station data in the target area. In this embodiment, the central European region (45°N - 55°N, 5°E - 25°E) is taken as an example, and the three-dimensional coordinates of 31 GNSS stations, the geoid with a resolution of 1°×1° in the region, and the surface elevation information with a resolution of 10′×10′ are collected. For relevant situations, see Figure 2 、 Figure 3 . Among them, Figure 2 the triangles in

[0049] are the positions of GNSS stations, the letters are the names of the stations, and the terrain undulation is represented by color blocks with different depths of color.

[0050] 1) According to the vertical profile meteorological data provided by the stratified grid meteorological model, and combining with the stratified integration method, calculate the zenith tropospheric delay values of the grids in the region, and use this as reference data (i.e., the zenith tropospheric delay reference value). The calculation formula of the stratified integration method is:

[0051]

[0052] In the formula, ZTD is the calculated zenith tropospheric delay value; surf and trop are the surface elevation and tropopause elevation respectively; k1, k'2, k3 are constant coefficients, which are equal to 77.6890K / hPa, 22.13K / hPa, 3.739×105K 2 / hPa; p, p w , T are the air pressure, water vapor pressure, and temperature of each height layer respectively. In this embodiment, the stratified grid meteorological data at UTC 00:00 every day in 2021 provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) is used to calculate the zenith tropospheric delay reference value with a resolution of 15′×15′.

[0053] 2) According to the methods in Steps 3 to 5 below, use the GNSS stations in the area and the barometric pressure values at the corresponding positions to complete the calculation of the tropospheric delay in the area, and subtract the tropospheric delay reference value to obtain the simulation preview accuracy. Moreover, all existing stations need to be included. After obtaining the preview accuracy, stations can be appropriately deleted according to requirements. Each time a station is deleted, a new simulation needs to be performed and the preview accuracy obtained. If the accuracy drops significantly, the deleted stations need to be restored. The preview accuracy in this embodiment is based on the annual mean absolute error (MAE) of each grid point in the area, and the specific accuracy distribution is shown in Figure 4 . Since there are some areas in this embodiment where the MAE is significantly larger (such as the southeast corner), if the calculation accuracy of the tropospheric delay is required to be better than 20 mm, then no more stations can be deleted in this embodiment.

[0054] 3) Conduct a regression analysis on the mean absolute error of each grid point and the shortest distance between each grid point and the station, and establish a linear function model. The distribution of the shortest distance between the grid points and the stations in this embodiment is shown in Figure 5 , and the established linear function model is:

[0055] A = 7.002×10 -2 D + 9.041

[0056] where A is the MAE, in millimeters; D is the shortest distance between each grid point and the stations in the area, in kilometers. The correlation coefficient between A and D is 0.8765. In this embodiment, if MAE < 20 mm is required, then according to the regression model, the distance D ≈ 165.978 km, that is, it is required that the distance of each position in the study area from the discrete stations does not exceed 165.978 km. In this way, appropriate stations can be selected, or temporary stations can be added in the areas that do not meet the requirements.

[0057] Step 3: According to the observables of the selected GNSS stations and the surface barometric pressure observables, calculate the zenith tropospheric delay values at each GNSS station (i.e., at the surface).

[0058] The calculation formulas for the zenith tropospheric dry and wet delays at the GNSS station are as follows:

[0059]

[0060] ZWD surf = ZTD surf - ZHD surf

[0061] where ZHD surf , ZWD surf are the zenith tropospheric dry and wet delay values at the surface respectively; ZTD surf is the total zenith tropospheric delay value, estimated from the GNSS pseudorange observation data; Psurf is the surface air pressure value, obtained from the barometric sensors of GNSS stations or nearby meteorological stations; B and H are the geodetic latitude and geodetic height of the GNSS station respectively. Taking January 4, 2021, UTC 00:00 as an example, Table 1 lists the coordinates, air pressure, and calculated zenith tropospheric dry and wet delay values of 29 GNSS stations in this embodiment (two stations in the area had no observation data at this moment).

[0062] Table 1 Coordinates, air pressure observations of GNSS stations in the survey area, and calculated zenith tropospheric dry and wet delay values

[0063]

[0064]

[0065] Step 4: Extend the zenith tropospheric delay at each GNSS station to the geoid (i.e., the surface with zero orthometric height), and after extension, perform gridding processing to obtain the regional delay value on the geoid.

[0066] 1) Extend the zenith tropospheric dry and wet delay values obtained in Step 3 to obtain the extended delay values at the GNSS stations. The extension functions are respectively:

[0067]

[0068]

[0069] In the formula, ZHD g and ZWD g are the zenith tropospheric dry and wet delay values after extension to the geoid respectively; ZWD surf is the surface wet delay value; H norm is the orthometric height to be extended; N is the geoid undulation at the position to be extended, which can be obtained according to the geoid information collected in this area.

[0070] 2) After extension, grid the discrete dry and wet delay values respectively. The gridding adopts the form of a quadratic surface function:

[0071] f(B,L) = a0 + a1B + a2L + a3B 2 + a4L 2 + a5BL

[0072] In the formula, B and L are the geodetic latitude and geodetic longitude of the discrete points respectively; a i(i = 1, 2, …, 5) are the coefficients of the function to be solved. By substituting the geodetic latitude and geodetic longitude of the GNSS station into the above formula, the coefficients can be obtained according to the least squares method. Accordingly, the discrete tropospheric delay values on the geoid can be transformed into a grid form. Taking January 4, 2021, UTC 00:00 as an example, the coefficients to be solved in this embodiment are shown in Table 2.

[0073] Table 2 Coefficients of the tropospheric dry and wet delay grid functions

[0074] component <![CDATA[a0]]> <![CDATA[a1]]> <![CDATA[a2]]> <![CDATA[a3]]> <![CDATA[a4]]> <![CDATA[a5]]> dry delay 15249.887 27.373 0.365 -553.902 5.878 -0.747 wet delay 1025.519 26.275 0.050 -45.440 0.519 -0.524

[0075] Step five, reverse extend the regional delay value on the geoid to the surface of the earth to obtain the regional zenith tropospheric delay value.

[0076] 1) Reverse extend the tropospheric dry and wet delays on the geoid. The formula is:

[0077] ZHD surf = ZHD g × (1 - 0.0000226 × H norm ) 5.225

[0078]

[0079] 2) Add the reversed dry and wet delay values to obtain the final regional tropospheric delay value.

[0080] It should be noted that in step three of this embodiment, when performing the extension process, it is extended to the geoid because through multiple experiments, it is found that the overall calculation accuracy is relatively high when extended to the geoid. Of course, if it is found that the accuracy is higher when extended to other reference surfaces, it can be extended to other reference surfaces.

[0081] In summary, according to the above steps, the present invention can obtain accurate zenith tropospheric delay results in the target area based on a small number of GNSS stations and surface meteorological observation data. Through simulation, it judges whether the distribution of stations in the area meets the accuracy requirements, and when the station positions do not meet the requirements, appropriate temporary stations are set according to the regression function relationship between the accuracy and the shortest distance index of the stations. Thus, the problem of difficult or inaccurate acquisition of tropospheric delay information caused by the shortage of stations can be successfully solved, which can provide effective technical support for the research and application in the field of space geodesy such as GNSS precise point positioning, very long baseline interferometry, and spaceborne radar interferometry.

Claims

1. A method for calculating the regional zenith tropospheric delay, characterized in that, It includes the following steps: 1) Select GNSS stations that meet the requirements for the calculation accuracy of zenith tropospheric delay in the target area according to the corresponding relationship between the calculation accuracy of zenith tropospheric delay and the shortest distance from the station; the following simulation calculations are performed on the target area to obtain the corresponding relationship: Calculate the zenith tropospheric delay value of the grid in the area based on the vertical profile meteorological data provided by the stratified grid meteorological model as the reference value of zenith tropospheric delay; Use the data of GNSS stations in the target area and the barometric data at the corresponding positions to calculate the zenith tropospheric delay value of the area, and compare it with the reference value of zenith tropospheric delay to obtain the error of each grid point; Perform regression analysis on the error of each grid point and the shortest distance between each grid point and the station to obtain the corresponding relationship; 2) Calculate the zenith tropospheric delay value including the zenith tropospheric dry delay value and the zenith tropospheric wet delay value at each GNSS station according to the observed values of the selected GNSS stations; 3) Extend the zenith tropospheric delay at each GNSS station to other reference surfaces, and perform gridding processing after extension to obtain the regional delay value on other reference surfaces; 4) Reverse extend the regional delay value on other reference surfaces to the ground surface to obtain the regional zenith tropospheric delay value.

2. The method for calculating the regional zenith tropospheric delay according to claim 1, characterized in that, It is necessary to combine the stratified integration method to calculate the zenith tropospheric delay value of the grid in the area according to the vertical profile meteorological data provided by the stratified grid meteorological model. The calculation formula of the stratified integration method is: wherein, ZTD is the calculated zenith tropospheric delay value; surf and trop are the surface elevation and tropopause elevation respectively; k1, k'2, k3 are constant coefficients; p, p w , and T are the air pressure, water vapor pressure, and air temperature at each altitude level respectively.

3. The regional zenith tropospheric delay calculation method according to claim 1, characterized in that, When the existing GNSS stations in the target area cannot meet the requirements for the calculation accuracy of zenith tropospheric delay in the target area, it is necessary to set up temporary GNSS stations in the target area to meet the requirements.

4. The method for calculating the regional zenith tropospheric delay according to claim 1, characterized in that, The other reference surface is the geoid.

5. The regional zenith tropospheric delay calculation method according to claim 1, characterized in that, In step 3), a quadratic surface function is used for gridding processing.

6. The method for calculating the regional zenith tropospheric delay according to claim 5, wherein The form of the quadratic surface function is: f(B,L) = a0 + a1B + a2L + a3B 2 + a4L 2 + a5BL where f(B,L) is a quadratic surface function with respect to geodetic latitude B and geodetic longitude L; a i is the function coefficient, and i = 1, 2, …, 5.

7. The method for calculating the regional zenith tropospheric delay according to claim 1, characterized in that The calculation formulas for the zenith tropospheric dry delay value and the zenith tropospheric wet delay value at each GNSS station are: ZWD surf = ZTD surf - ZHD surf Wherein, ZHD surf and ZWD surf are respectively the dry and wet zenith tropospheric delay values at the GNSS station; P surf is the surface air pressure value; B and H are respectively the geodetic latitude and geodetic height at the GNSS station; ZTD surf are respectively the zenith tropospheric delay values at the GNSS station.

8. The method for calculating the regional zenith tropospheric delay according to claim 4, wherein The extension formula in step 3) is: wherein, ZHD g , ZWD g are respectively the zenith tropospheric dry and wet delay values after being extended to the geoid; P surf is the surface air pressure value; B is the geodetic latitude at the GNSS station; H norm is the orthometric height to be extended; N is the geoid undulation at the position to be extended; ZWD surf is the zenith tropospheric wet delay value at the surface.

9. The method for calculating the regional zenith tropospheric delay according to claim 8, characterized in that, The reverse extension formula in step 4) is: ZHD surf = ZHD g × (1 - 0.0000226 × H norm ) 5.225 In the formula, ZHD surf and ZWD surf are the dry and wet zenith tropospheric delay values of the ground surface, respectively.

Images