A navigation and positioning method based on high-precision ZTD interpolation

Abstract

This application relates to the field of satellite navigation and positioning methods. Specifically, it relates to a navigation and positioning method based on high-precision ZTD interpolation, which can improve the accuracy of GNSS navigation and positioning to a certain extent. The navigation and positioning method includes: determining the ZHD elevation correction formula and the ZWD elevation correction formula; determining the first number of ERA5 grid points around the GNSS station, obtaining the ERA5ZHD and ERA5ZWD at each ERA5 grid point, obtaining the ERA5ZWD of the GNSS station height and the ERA5ZHD of the GNSS station height; based on the acquisition of the ERA5ZHD at the GNSS station location and the ERA5ZWD at the GNSS station location, obtaining the GNSS ZHD and GNSSZWD; obtaining the GNSSZHD and GNSSZWD at the first position on the mean elevation plane and performing interpolation at any position within the mean elevation plane through spherical harmonic functions to obtain the GNSSZHD and GNSSZWD at the second position on the mean elevation plane; based on the ZHD elevation correction formula and the ZWD elevation correction formula, respectively obtaining the GNSSZHD and GNSSZWD at the second height, and then obtaining the GNSSZTD at the second height for subsequent precise positioning and navigation.

Description

Technical Field

[0001] This application relates to the field of satellite navigation and positioning methods. Specifically, it relates to a navigation and positioning method based on high-precision ZTD interpolation. Background Technique

[0002] Tropospheric delay is one of the most important error sources in space geodetic surveying, and it is also an important factor restricting the navigation and positioning of high-precision Global Navigation Satellite System (GNSS). During the propagation of radio signals, the tropospheric delay of radio signals from the zenith direction to the horizon direction can reach 2m to 20m, seriously restricting the ambiguity convergence speed and positioning accuracy of GNSS precise point positioning and medium- and long-distance baseline differential positioning.

[0003] In the existing technologies of GNSS data processing, it is usually necessary to weaken the tropospheric delay. Generally, the zenith tropospheric delay (ZTD) is estimated through an empirical model, and then subsequent processing is carried out. Among them, ZTD is the average value obtained by projecting the slant path delay values in different directions to the vertical direction using a mapping function. According to the cause, ZTD can be further divided into zenith hydrostatic delay (ZHD) and zenith wet delay (ZWD). ZHD accounts for 90% or even more of the total atmospheric delay, and generally can be accurately estimated through surface pressure using an empirical model. While ZWD is less than 10% of the total atmospheric delay, but it seriously affects the accurate value of ZTD. Therefore, when obtaining ZTD, it is necessary to accurately calculate ZHD and ZWD.

[0004] However, since ZWD is greatly affected by the spatio-temporal distribution of atmospheric water vapor, ZTD cannot be directly obtained by simply adding the model values of ZHD and ZWD. Thus, it is difficult to accurately estimate ZTD using the model, which in turn leads to the accuracy of GNSS navigation and positioning. Summary of the Invention

[0005] In order to improve the accuracy of GNSS navigation and positioning, this application provides a navigation and positioning method based on high-precision ZTD interpolation.

[0006] The embodiments of this application are implemented as follows:

[0007] The embodiments of this application provide a navigation and positioning method based on high-precision ZTD interpolation, including:

[0008] Determine the ZHD elevation correction formula and the ZWD elevation correction formula, where the ZHD elevation correction formula is used to describe the relationship between ZHDs corresponding to different heights; the ZWD elevation correction formula is used to describe the relationship between ZWDs corresponding to different heights;

[0009] Determine the first quantity of ERA5 grid points near the GNSS site, obtain the ZHD and ZWD at each ERA5 grid point, and through the ZWD elevation correction formula, unify the ZHD at each ERA5 grid point to the ZHD at the GNSS site height, denoted as the GNSS site height ZHD; through the ZWD elevation correction formula, unify the ZWD at each ERA5 grid point to the ZWD at the GNSS site height, denoted as the GNSS site height ZWD;

[0010] Based on the GNSS site height ZWD and the GNSS site height ZHD, obtain the ZHD at the GNSS site location and the ZWD at the GNSS site location, and based on the ZHD at the GNSS site location, the ZWD at the GNSS site location, and the ZTD data at the GNSS site location, obtain the GNSS ZHD and the GNSS ZWD;

[0011] Based on the ZWD elevation correction formula, unify the GNSS ZWD to the mean elevation surface; based on the ZHD elevation correction formula, unify the GNSS ZHD to the mean elevation surface, and obtain the GNSS ZHD and the GNSS ZWD at the first position on the mean elevation surface;

[0012] Interpolate at any position within the mean elevation surface through spherical harmonic functions to obtain the GNSS ZHD and the GNSS ZWD at the second position on the mean elevation surface, based on the ZWD elevation correction formula, unify the GNSS ZWD at the second position to the GNSS ZWD at the second height; based on the ZHD elevation correction formula, unify the GNSS ZHD at the second position to the GNSS ZHD at the second height, and obtain the GNSS ZTD at the second height.

[0013] In some embodiments, determine the empirical model coefficients in the ZHD elevation correction formula and the ZWD elevation correction formula through ZHD and ZWD;

[0014] Wherein, the ZHD and the ZWD are calculated based on the parameter set provided by ERA5.

[0015] In some embodiments, in the step of obtaining the ZHD of the GNSS site location and the ZWD of the GNSS site location based on the GNSS site height ZWD and the GNSS site height ZHD, the ZHD of the GNSS site location and the ZWD of the GNSS site location are obtained by means of bilinear interpolation.

[0016] In some embodiments, in the step of obtaining the GNSS ZHD and the GNSS ZWD based on the ZHD of the GNSS site location, the ZWD of the GNSS site location, and the ZTD of the GNSS site location,

[0017] Calculate the ratio of ZHD to ZWD in the GNSS ZTD at the GNSS site location to determine the proportionality coefficient, where the ratio is equivalent to the ratio of the ERA5 ZHD of the GNSS site location to the ERA5 ZWD of the GNSS site location;

[0018] Determine the GNSS ZHD and the GNSS ZWD of the GNSS site according to the proportionality coefficient.

[0019] In some embodiments, in the step of obtaining the GNSS ZHD and the GNSS ZWD based on the ZHD of the GNSS site location, the ZWD of the GNSS site location, and the ZTD of the GNSS site location,

[0020] Calculate the monthly average value of the ZWD at the GNSS site location and the monthly average value of the ZHD at the GNSS site location, and calculate the proportionality coefficient in the GNSS ZTD month by month.

[0021] In some embodiments, in the step of unifying the GNSS ZWD to the mean elevation surface based on the ZWD elevation correction formula, calculate the average value of the heights at all grid points in the ERA5 data, or calculate the average value of all the heights of the GNSS sites to obtain the mean elevation surface.

[0022] In some embodiments, in the step of interpolating at any position within the mean elevation surface by means of spherical harmonic functions to obtain the GNSS ZHD and the GNSS ZWD at the second position on the mean elevation surface, further comprising:

[0023] Obtain the maximum order of the spherical harmonic function;

[0024] Input the maximum order, the longitude and latitude of the first position, and the GNSS ZHD and the GNSS ZWD of the first position to unify the GNSS data to any position on the mean elevation surface by means of spherical harmonic functions.

[0025] Input the longitude and latitude of the second location to output the GNSS ZHD and GNSS ZWD of the second location

[0026] In some embodiments, in the step of obtaining the maximum order of the spherical harmonic function, it further includes:

[0027] Adopt the leave-one-out cross-validation method, fit the ZTD and ZWD of each station with spherical harmonic functions of different orders, and calculate the RMS of the deviation between it and the actual ZHD and ZWD of the station. The order with the minimum RMS is the optimal maximum order of the spherical harmonic function.

[0028] In some embodiments, in the step of unifying the ZWD at the second location to the ZWD at the second height based on the ZWD elevation correction formula, and unifying the ZHD at the second location to the ZHD at the second height based on the ZHD elevation correction formula, where the height of the average elevation plane is the first height, it further includes:

[0029] Substitute the first height, the second height, and the ZHD at the first height into the ZHD elevation correction formula to obtain the ZHD at the second height;

[0030] Substitute the first height, the second height, and the ZWD at the first height into the ZWD elevation correction formula to obtain the ZWD at the second height.

[0031] The GNSS station ZTD at the second height is the sum of the ZHD at the second height and the ZWD at the second height.

[0032] Advantages of the present application: By performing elevation vertical correction on the separated GNSS ZHD and GNSS ZWD and spherical harmonic fitting of the specified elevation plane, it helps to achieve ZTD interpolation at any spatial position within the region, can conveniently and accurately obtain the ZTD at any spatial position, so that GNSS navigation and positioning are more accurate; further, by considering the influence of the region and time on the elevation correction model coefficients, calculate the elevation correction model coefficients of ZHD and ZWD, and then based on the elevation correction model coefficients of ZHD and ZWD, the ZHD and ZWD at different heights can be accurately converted; further, by calculating the monthly average of the ZWD at the GNSS station location and the monthly average of the ZHD at the GNSS station location, it not only helps to reduce the influence of the climate conditions of each month on the test data, but also helps to reduce the subsequent calculation workload. Description of the Drawings

[0033] To more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are some embodiments of the present application. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0034] Figure 1 It is a flowchart of a navigation and positioning method based on high-precision ZTD interpolation according to one or more embodiments of the present application;

[0035] Figure 2 It is a specific method flowchart of step 100 in the navigation and positioning method based on high-precision ZTD interpolation according to one or more embodiments of the present application;

[0036] Figure 3 It is a specific method flowchart of step 300 in the navigation and positioning method based on high-precision ZTD interpolation according to one or more embodiments of the present application;

[0037] Figure 4 It is a specific method flowchart of step 500 in the navigation and positioning method based on high-precision ZTD interpolation according to one or more embodiments of the present application. Detailed implementation manners

[0038] To make the purpose, implementation manners and advantages of the present application clearer, the following will clearly and completely describe the exemplary implementation manners of the present application with reference to the drawings in the exemplary embodiments of the present application. Obviously, the described exemplary embodiments are only some embodiments of the present application, rather than all embodiments.

[0039] It should be noted that the brief description of the terms in the present application is only for the convenience of understanding the following described implementation manners, rather than intending to limit the implementation manners of the present application. Unless otherwise stated, these terms should be understood in their ordinary and common meanings.

[0040] The terms "first", "second", "third", etc. in the specification, claims and the above drawings of the present application are used to distinguish similar or like objects or entities, and do not necessarily mean to limit a specific order or sequence, unless otherwise noted. It should be understood that such terms can be interchanged under appropriate circumstances.

[0041] The terms "comprising" and "having" and any variations thereof are intended to cover but not be exclusive of inclusion. For example, a product or device comprising a series of components does not necessarily have to be limited to all the components clearly listed, but may include other components not clearly listed or inherent to these products or devices.

[0042] Currently, the methods for obtaining GNSS ZTD generally adopt the following ways:

[0043] Using the ray tracing principle, the vertical profile meteorological data above the known GNSS stations and the nearby areas is utilized, and GNSS ZTD is calculated through the integration method. However, this acquisition method has a high cost and a long cycle for obtaining results.

[0044] Or GNSS ZTD is obtained based on the empirical models of meteorological parameters and non-meteorological parameters. Compared with the integration method, the method of obtaining GNSS ZTD through the empirical model greatly simplifies the calculation process of GNSS ZTD, but the calculation accuracy of GNSS ZTD also decreases accordingly.

[0045] Or the observation results of space geodetic techniques such as Very Long Baseline Interferometry (VLBI) and GNSS are utilized. Parameters such as the dry delay (ZHD) and mapping function in the tropospheric delay are regarded as prior values, and the wet delay (ZWD) is regarded as an unknown parameter for estimation. However, this method requires permanent GNSS stations to be set up in the research area, and for areas with sparse station distributions, this method does not have high application value.

[0046] Or the global satellite navigation positioning system is used to obtain high-precision GNSS ZTD estimates. This method has advantages such as all-weather, high time resolution, and high accuracy. However, there are problems that the number of GNSS stations is scarce and unevenly distributed, resulting in the GNSS ZTD data obtained being discretely distributed in space, which limits related research in terms of spatial resolution.

[0047] To address the above problems, this application takes into account the influence of region and time on the coefficients of the elevation correction model, and uses the ZHD and ZWD of ERA5 grid points to calculate the elevation correction model coefficients that are more suitable for the ZHD and ZWD in this area. Secondly, by applying the ratio of ERA5 ZHD at the GNSS station location to ERA5 ZWD at the GNSS station location to GNSS, the separation of GNSS wet and dry delays is realized. Finally, the elevation vertical correction and spherical harmonic fitting of the specified elevation surface are performed on the separated GNSS ZHD and GNSS ZWD to achieve ZTD interpolation at any spatial position within the region. Through the above method, the GNSS ZTD at any spatial position within the region can be accurately obtained, which helps to improve the convenience and universality of accurately obtaining GNSS ZTD, is applicable to the acquisition of high-precision ZTD in large ranges and large height differences, and is more conducive to enabling GNSS to achieve accurate navigation and positioning based on a large amount of high-precision ZTD data. The following specifically elaborates on the navigation and positioning method based on high-precision ZTD interpolation:

[0048] Figure 1 This is a flowchart of a navigation and positioning method based on high-precision ZTD interpolation in an embodiment. As Figure 1 shown, the navigation and positioning method based on high-precision ZTD interpolation specifically includes:

[0049] In step 100, determine the ZHD elevation correction formula and the ZWD elevation correction formula, where both the ZHD elevation correction formula and the ZWD elevation correction formula include ERA5 data.

[0050] Among them, since the magnitude of the ZTD value is closely related to the height, and as the height increases, the ZTD value continuously decreases, it is necessary to use the ZHD elevation correction formula and the ZWD elevation correction formula to complete the spatial matching of ERA5 ZTD and GNSS ZTD.

[0051] It should be noted that since the European Centre for Medium-Range Weather Forecasts reanalysis data ERA5 can provide a continuously updated numerical description of the climate, and the ERA5 data is easily available and has a large amount of data, the ERA5 data values include the estimation of atmospheric parameters such as air temperature, pressure, and wind force in different altitude regions, as well as surface parameters such as rainfall, soil moisture content, and sea wave height. Therefore, in some embodiments, the empirical model coefficients in the ZHD elevation correction formula and the ZWD elevation correction formula are determined through ERA5 ZHD and ERA5 ZWD; and ZHD and ZWD are calculated based on the parameter group provided by ERA5, where the parameter group provided by ERA5 includes atmospheric pressure, atmospheric weighted average temperature, precipitable water vapor (PWV), station latitude, station height, and liquid water density, etc. As Figure 2 shown, the specific determination process of the ZWD elevation correction formula and the ZHD elevation correction formula is as follows:

[0052] In step 110, calculate ERA5 ZHD using the atmospheric pressure data, and the formula is as follows (Equation 1):

[0053]

[0054] In the formula, P represents the surface atmospheric pressure (unit: hPa), φ represents the station latitude (unit: rad), and H represents the station height (unit: km).

[0055] In step 120, calculate ERA5 ZWD using the atmospheric pressure data, and the formula is as follows (Equation 2):

[0056]

[0057] In the formula, T m represents the weighted average temperature (unit: K), and ρ represents the liquid water density (unit: kg / m 3), R v represents the specific gas constant for water vapor (461.495 J / kg / K), k'2 and k3 are physical constants, and the value of k'2 is generally 16.48 · K / hPa, and the value of k3 is generally (3.776 ± 0.014) × 10 5 K 2 / hPa.

[0058] Step 130: Determine the model coefficients in the ZHD elevation correction formula and the ZWD elevation correction formula based on the known ERA5 ZWD and ERA5 ZHD.

[0059] Among them, substitute ERA5 ZHD at different heights into Equation (3a) to calculate the empirical model coefficient α, and substitute ERA5 ZWD at different heights into Equation (3b) to calculate the empirical model coefficient β. It should be noted that the values of α and β are related to the area to be studied.

[0060]

[0061]

[0062] Among them, ZHD h1 , ZWD h1 and ZHD h2 , ZWD h2 respectively represent the dry delay and wet delay corresponding to heights h1 and h2. The units of h1 and h2 are both m, and i represents the month number.

[0063] Step 140: Substitute the α and β determined in Step 130 into Equation (3a) and Equation (3b) to obtain the ZHD elevation correction formula and the ZWD elevation correction formula.

[0064] It should be noted that for most areas to be studied, the value of α can be 5.225, and the value of β can be 2000. The above ZHD elevation correction formula and ZWD elevation correction formula fully consider influencing factors such as spatial position and season, ensuring the accuracy and robustness of the method. The ZHD elevation correction formula describes the relationship between ZHD corresponding to different heights, and the ZWD elevation correction formula describes the relationship between ZWD corresponding to different heights.

[0065] In step 200, determine the ERA5 grid points around the GNSS station at the first height. The number of the ERA5 grid points is the first quantity. Obtain the ERA5 ZHD and ERA5 ZWD at each ERA5 grid point. Based on the height at each of the ERA5 grid points, through the ZHD elevation correction formula and the ZWD elevation correction formula, unify the ERA5 ZHD and ERA5 ZWD at each ERA5 grid point into the ERA5 ZHD and ERA5 ZWD at the first height respectively.

[0066] Among them, the ERA5 grid is a planar grid with a certain height. ERA5 provides global planar data with a spatial resolution of 0.25°. Therefore, each ERA5 grid generally has four grid points, that is, the first quantity is usually four. It can be understood that it is necessary to find four ERA5 grid points near the GNSS station, record the parameter groups at each ERA5 grid point, and obtain the ERA5 ZHD and ERA5 ZWD at each grid point through the above formulas (1) and (2).

[0067] After that, since the ERA5 grid also has a certain height and the height of the GNSS station is generally different from the height of the ERA5 grid points, through the above ZHD elevation correction formula (3a) and ZWD elevation correction formula (3b), unify the ERA5 ZHD and ERA5 ZWD at the ERA5 grid into the ERA5 ZHD at the first height and the ERA5 ZWD at the first height respectively, and obtain the ERA5 ZHD and ERA5 ZWD at the height of the GNSS station.

[0068] In step 300, based on the first quantity of the ERA5 ZHD at the first height and the first quantity of the ERA5 ZWD at the first height, obtain the ERA5 ZHD at the GNSS station location and the ERA5 ZWD at the GNSS station location by means of bilinear interpolation; calculate the proportionality coefficient, and obtain the GNSS ZHD and GNSS ZWD at the GNSS station location according to the relationship between the proportional relationship and the GNSS ZTD at the GNSS station location, where the GNSS ZTD at the GNSS station location is directly obtained by the GNSS station.

[0069] Among them, bilinear interpolation, also known as bilinear interpolation, is an interpolation algorithm in numerical analysis. Bilinear interpolation is an extension of linear interpolation for interpolation functions with two variables. Its core idea is to perform linear interpolation in two directions respectively. In this embodiment, by performing linear interpolation on the ERA5 ZHD and ERA5 ZWD at four grid points in the longitude of the GNSS station and then performing linear interpolation in the latitude of the GNSS station, at this time, the purpose of unifying and matching the height, accuracy, and latitude at the GNSS station with those at the ERA5 grid is achieved. It can be understood that the ERA5 ZHD at the GNSS station position and the ERA5 ZWD at the GNSS station position obtained at this time are the ZHD and ZWD using the GNSS station position information.

[0070] In some embodiments, as Figure 3 shown, in the process of determining the ERA5 ZHD at the GNSS station position and the ERA5 ZWD at the GNSS station position, since the time resolution of the data volume of ERA5 ZHD and ERA5 ZWD is in hours. For example, when data is collected every hour, the cumulative data volume for one month is 24×30. Therefore, in order to reduce the computational load in the subsequent process, the process of determining the GNSS ZHD and GNSS ZWD of the GNSS station specifically includes:

[0071] Step 310, by adding the corresponding data values of each month and then dividing by the data volume of that month, calculate the monthly average value of ERA5 ZHD at the GNSS station position and the monthly average value of ERA5 ZWD at the GNSS station position.

[0072] Step 320, calculate the ratio of ERA5 ZHD at the GNSS station position to ERA5 ZWD at the GNSS station position month by month with reference to formula (4), and the proportionality coefficient k m can be determined. This proportionality coefficient can be used to characterize the ratio of GNSS ZHD to GNSS ZWD in GNSS ZTD. Since ERA5 ZHD and ERA5 ZWD are related to factors such as water vapor content, rainfall, and weather humidity in a month, calculating the monthly average value of ERA5 ZHD at the GNSS station position and the monthly average value of ERA5 ZWD at the GNSS station position can reduce the influence of climate factors on ERA5 ZHD and ERA5 ZWD within a month.

[0073]

[0074] In the formula, ZHD ERA5 and ZWD ERA5 respectively represent the dry delay and wet delay of the GNSS station position retrieved by ERA5, and m represents the month.

[0075] Step 330: Separate the GNSS ZTD according to formulas (5) and (6), that is, determine the GNSS ZHD and GNSS ZWD of the GNSS station according to the proportionality coefficient km. It can be understood that the acquisition of GNSS ZHD and GNSS ZWD is obtained by separating the GNSS ZTD based on ERA5 data.

[0076]

[0077]

[0078] In the formula, and respectively represent the ZTD, ZHD, and ZWD of the GNSS station, k m represents the proportionality coefficient, m represents the month, and m = 1, 2, 3,..., 12.

[0079] It should be noted that in steps 100 - 300, since the number of GNSS ZTD is limited, only a limited number of GNSS ZHD and GNSS ZWD can be obtained. Therefore, in order to accurately obtain the GNSS ZHD and GNSS ZWD at any location, the following steps are also required.

[0080] In step 400, based on the ZWD elevation correction formula, unify the GNSS ZWD at the GNSS station location to the average elevation plane to obtain the GNSS ZHD at the first location on the average elevation plane; based on the ZHD elevation correction formula, unify the GNSS ZHD at the GNSS station location to the average elevation plane to obtain the GNSS ZWD at the first location on the average elevation plane.

[0081] It can be understood that when the spherical harmonic function fits the site data to the planar data, the influence of height is not considered, and the ERA5 grid is a curved surface with a certain height. Therefore, in order to better unify the ERA5 grid curved surface to the plane, it is necessary to calculate the average elevation plane. In some embodiments, the average elevation plane is obtained by averaging the heights at all grid points in ERA5. In other embodiments, the average elevation plane is obtained by averaging all the heights of the GNSS stations. It should be noted that since the amount of data of the GNSS stations is small, the effect of obtaining the average elevation plane by averaging the heights at all grid points in ERA5 will be better.

[0082] After obtaining the mean elevation surface, since the height of the mean elevation surface is known, the GNSS ZWD obtained in step 300 above can be unified to the mean elevation surface through the ZWD elevation correction formula to obtain the GNSS ZHD at the first position on the mean elevation surface; through the ZHD elevation correction formula, the GNSS ZHD obtained in step 300 above is unified to the mean elevation surface to obtain the GNSS ZHD at the first position on the mean elevation surface.

[0083] In step 500, interpolation is performed at any position within the mean elevation surface using spherical harmonic functions to obtain the GNSS ZHD and GNSS ZWD at the second position on the mean elevation surface. Based on the ZWD elevation correction formula and the ZHD elevation correction formula, the GNSS ZWD at the second position and the GNSS ZHD at the second position are respectively unified to the GNSS ZWD at the second height and the GNSS ZHD at the second height; wherein, the second height is different from the height of the mean elevation surface, and the second height is the GNSS station height; the GNSS ZHD at the second height and the GNSS ZWD at the second height are added together to obtain the GNSS ZTD at the second height.

[0084] Wherein, when the height of the mean elevation surface is the third height, the heights of both the first position and the second position are the third height, and the second height is the GNSS station height required in the actual application process.

[0085] It should be noted that during the process of interpolation at any position within the mean elevation surface using spherical harmonic functions, as Figure 4 shown, it specifically includes the following steps:

[0086] Step 510, obtain the maximum order of the spherical harmonic function;

[0087] Wherein, the calculation formula of the spherical harmonic function is as shown in formula (7):

[0088]

[0089] Wherein, M and N respectively represent the maximum degree and order of the spherical harmonic function. And the calculation formula of a ij is formula (8), and the calculation formula of b ij is formula (9):

[0090]

[0091]

[0092] In the formula, φ and λ respectively represent the latitude and longitude of the station or grid point, and Pn(x) represents the Legendre function, and its function expression is as shown in the following formula (10):

[0093]

[0094] In the formula, represents the integer part of, then formula (7) can be converted to formula (11):

[0095]

[0096] In formula (11), A represents the coefficient to be estimated of the spherical harmonic function model, which can be determined by the finally determined spherical harmonic function model itself.

[0097] It can be seen from the above analysis that in the process of using the spherical harmonic function, when it is necessary to determine the optimal maximum degree and order, the leave-one-out cross-validation method can be adopted, that is, each time one GNSS site is left out, and the GNSS ZHD and GNSS ZWD of the remaining GNSS sites are used to fit the GNSS ZHD and GNSS ZWD of this site through spherical harmonic functions of different orders, and the root mean square (RMS) of the deviation between the above fitting values and the actual GNSS ZHD and GNSS ZWD of this GNSS site is calculated. The order with the smallest root mean square is the optimal maximum order of the spherical harmonic function.

[0098] Step 520, after obtaining the maximum order of the spherical harmonic function, input the maximum order, the longitude, latitude of the first position, and the GNSS ZHD and GNSS ZWD of the first position into the spherical harmonic function model, so as to unify the GNSS data to any position on the mean elevation surface through the spherical harmonic function.

[0099] Step 530, input the longitude and latitude of the second position into the spherical harmonic function model, and output the GNSS ZHD and GNSS ZWD of the second position through the operation of the spherical harmonic function model.

[0100] In step 600, based on performing multiple interpolations on the mean elevation surface, obtain the GNSS ZTD at different heights of the GNSS site for GNSS navigation and positioning.

[0101] The beneficial effects of the embodiments in this part are as follows. By considering the influence of region and time on the coefficients of the elevation correction model, the coefficients of the elevation correction model for ZHD and ZWD are calculated. Then, based on the coefficients of the elevation correction model for ZHD and ZWD, the ZHD and ZWD at different heights can be accurately converted. Further, by applying the ratio of ERA5 ZHD and ERA5 ZWD to GNSS, the separation of ZHD and ZWD in GNSS can be achieved efficiently and accurately, which helps to accurately obtain GNSS ZHD and GNSS ZWD. Further, by performing elevation vertical correction and spherical harmonic fitting on the specified elevation surface for the separated GNSS ZHD and GNSS ZWD, the interpolation of ZTD at any spatial position within the region can be realized, and the GNSS ZTD at any spatial position can be obtained conveniently and accurately. Furthermore, it helps the GNSS station to perform precise navigation and positioning based on the GNSS ZTD at multiple positions.

[0102] For the sake of convenience in explanation, the above description has been made in conjunction with specific embodiments. However, the above discussion in some embodiments is not intended to be exhaustive or to limit the embodiments to the specific forms disclosed above. According to the above teachings, various modifications and variations can be obtained. The selection and description of the above embodiments are for better explaining the principles and practical applications, so that those skilled in the art can better use the embodiments and various different modified embodiments suitable for specific use considerations.

Claims

1. A navigation and positioning method based on high-precision ZTD interpolation, characterized in that Including: Determine the ZHD elevation correction formula and the ZWD elevation correction formula, where both the ZHD elevation correction formula and the ZWD elevation correction formula include ERA5 data; Determine the ERA5 grid points around the GNSS site at the first height, the number of the ERA5 grid points being the first quantity, obtain the ERA5 ZHD and ERA5 ZWD at each ERA5 grid point, and based on the height at each ERA5 grid point, unify the ERA5 ZHD and ERA5 ZWD at each ERA5 grid point into the ERA5 ZHD and ERA5 ZWD at the first height respectively through the ZHD elevation correction formula and the ZWD elevation correction formula; Based on the first quantity of ERA5 ZHD at the first height and the first quantity of ERA5 ZWD at the first height, obtain the ERA5 ZHD and ERA5 ZWD at the GNSS site location through bilinear interpolation; calculate the proportionality coefficient, and obtain the GNSS ZHD and GNSS ZWD at the GNSS site location according to the relationship between the proportionality relationship and the GNSS ZTD at the GNSS site location, where the GNSS ZTD at the GNSS site location is directly obtained by the GNSS site; Based on the ZWD elevation correction formula and the ZHD elevation correction formula, unify both the GNSS ZWD at the GNSS site location and the GNSS ZHD at the GNSS site location to the mean elevation surface, and obtain the GNSS ZHD at the first position on the mean elevation surface and the GNSS ZWD at the first position; Interpolate at any position within the mean elevation surface through spherical harmonic functions to obtain the GNSS ZHD and GNSS ZWD at the second position on the mean elevation surface, and based on the ZWD elevation correction formula and the ZHD elevation correction formula, unify the GNSS ZWD at the second position and the GNSS ZHD at the second position to the GNSS ZWD at the second height and the GNSS ZHD at the second height respectively, where the second height is different from the height of the mean elevation surface, and the second height is the height of the GNSS site; add the GNSS ZHD at the second height and the GNSS ZWD at the second height to obtain the GNSS ZTD at the second height; Based on obtaining the GNSS ZTD at different heights of the GNSS site after multiple interpolations on the mean elevation surface, perform GNSS navigation positioning.

2. The navigation and positioning method based on high-precision ZTD interpolation according to claim 1, wherein In the step of determining the ZHD elevation correction formula and the ZWD elevation correction formula, it further includes: Fit the ERA5 ZHD and ERA5 ZHD at different heights to determine the empirical model coefficients in the ZHD elevation correction formula and the ZWD elevation correction formula; Wherein, both the ERA5 ZHD and the ERA5 ZWD are calculated from the parameter group provided by ERA5.

3. The navigation and positioning method based on high-precision ZTD interpolation according to claim 2, wherein The ZHD elevation correction formula is as follows: The ZHD elevation correction formula is: Among them, ZHD h1 , ZWD h1 and ZHD h2 , ZWD h2 respectively represent the corresponding dry delay and wet delay at heights h1 and h2, i represents the number of months; both α and β are the empirical model coefficients.

4. The navigation and positioning method based on high-precision ZTD interpolation according to claim 1, wherein: The proportionality coefficient is the ratio of ERA5 ZHD at the GNSS station location to ERA5 ZWD at the GNSS station location, and this ratio is equivalent to the ratio of ZHD to ZWD in GNSS ZTD at the GNSS station location.

5. The navigation and positioning method based on high-precision ZTD interpolation according to claim 4, characterized in that: Calculate the monthly mean of ERA5 ZWD at the GNSS station location and the monthly mean of ERA5 ZHD at the GNSS station location, and calculate the proportionality coefficient month by month.

6. The navigation and positioning method based on high-precision ZTD interpolation according to claim 1, characterized in that In the step of obtaining GNSS ZHD and GNSS ZWD at the GNSS station location according to the relationship between the proportional relationship and GNSS ZTD at the GNSS station location, it further includes: calculating GNSS ZHD and GNSS ZWD at the GNSS station location according to the following formula Among them, and represent the ZTD, ZHD, and ZWD of the GNSS station respectively, k m represents the proportionality coefficient, m represents the month, and m = 1, 2, 3,..., 12.

7. The navigation and positioning method based on high-precision ZTD interpolation according to claim 1, characterized in that Before the step of unifying GNSS ZWD to the mean elevation surface based on the ZWD elevation correction formula, it further includes: Calculate the average value of the heights at all grid points in the ERA5 data, or calculate the average value of all the heights of the GNSS stations, to obtain the mean elevation surface.

8. The navigation and positioning method based on high-precision ZTD interpolation according to claim 1, wherein, In the step of interpolating at any position within the mean elevation surface through spherical harmonic functions to obtain GNSS ZHD and GNSS ZWD at the second position on the mean elevation surface, it further includes: Obtain the maximum degree of the spherical harmonic function; Input the maximum degree, the longitude, latitude of the first position, and GNSS ZHD and GNSS ZWD at the first position, so as to unify GNSS data to any position on the mean elevation surface through spherical harmonic functions; Input the longitude and latitude of the second position to output GNSS ZHD and GNSS ZWD at the second position.

Images