A method for tropospheric delay correction taking into account non-isotropy

Abstract

The present invention belongs to the field of GNSS meteorology, and discloses a tropospheric delay correction method that takes non-isotropy into account. In view of the problem that most traditional methods are bound to isotropy or anisotropy when obtaining SPD, which limits the prediction accuracy of SPD, the present invention introduces the theory of non-isotropy and improves the calculation method of SPD. The method starts with the observation data between the measuring station and the satellite that can be obtained, and separates the unmodelable error and the modelable error from the residual through the PPP positioning algorithm; then, the non-isotropic tropospheric delay correction value is calculated by introducing a method based on least squares to solve a linear equation group, and a solution that minimizes the error is found to deal with the situation where the equation is redundant or cannot be accurately solved. Finally, the obtained non-isotropic tropospheric delay correction value is added to the SPD calculated by the mapping function to obtain the corrected SPD. The present invention is conducive to improving the prediction accuracy of SPD.

Description

Technical Field

[0001] The invention belongs to the field of GNSS meteorology and relates to a tropospheric delay correction method taking into account non-isotropy. Background Art

[0002] Tropospheric delay refers to the change in the path of electromagnetic wave signals emitted by Global Navigation Satellite System (GNSS) satellites when passing through the neutral atmosphere, which causes the signals to be refracted and delayed. It is an important error source that affects the accuracy of GNSS positioning. In order to improve the accuracy of GNSS and accelerate the convergence speed of Precise Point Positioning (PPP), it is necessary to obtain high-precision tropospheric slant path delay (SPD).

[0003] The current methods for calculating SPD mainly include ray tracing, mapping function and parameter estimation. Among them, the ray tracing method has the highest accuracy, but because its data source depends on local meteorological data, and is limited by the fact that meteorological data cannot be obtained and updated in real time, this method cannot obtain SPD in real time. The mapping function method and the SPD with horizontal gradient correction respectively reflect the characteristics of isotropy and anisotropy, but they are different from the SPD calculated by the most accurate ray tracing method. Therefore, real-time or predicted acquisition of high-precision SPD has become the key to high-precision real-time positioning.

[0004] At present, there are two main ways to predict SPD. One is to use neural networks or deep learning to train existing data to achieve prediction. The other is to perform parameter estimation modeling based on the SPD law and combined with a large amount of data.

[0005] For example, Patent Document 1 discloses a global tropospheric delay modeling method based on deep learning. This method uses a deep learning tropospheric atmospheric delay meteorological parameter model, and then regenerates the meteorological parameters required for tropospheric delay calculation, and finally realizes the prediction of SPD with anisotropic characteristics. Patent Document 2 discloses a tropospheric delay prediction method for a global navigation satellite system. This method obtains the parameters of the zenith dry delay model and the zenith wet delay model corresponding to each grid point by least squares fitting, and improves the elevation and prediction time according to the exponential change of the tropospheric delay in the zenith direction with the altitude, and obtains the zenith total delay model based on altitude and time. Finally, the zenith total delay model is multiplied by the mapping function to obtain the final SPD with isotropic characteristics.

[0006] It can be seen that the current SPD prediction can meet some real-time positioning needs. However, the characteristics of the SPD obtained by the current prediction are mostly constrained by isotropy or anisotropy, which limits the accuracy of SPD prediction. At the same time, when using neural network training, due to the influence of wild values, multiple trainings are required to achieve the expected optimal results.

[0007] Related Literature

[0008] Patent Document 1 China invention patent application publication number: CN111382507A, publication date: 2020.07.07;

[0009] Patent document 2 Chinese invention patent application publication number: CN115980809A, publication date: 2023.04.18. Summary of the invention

[0010] In view of the above technical problems, the present invention proposes a tropospheric delay correction method taking into account non-isotropy. The method introduces the theory of non-isotropy and helps to improve the prediction accuracy of SPD by improving the calculation method of SPD.

[0011] In order to achieve the above object, the present invention adopts the following technical scheme:

[0012] A method for correcting tropospheric delay taking into account non-isotropy comprises the following steps:

[0013] Step 1. Calculate the phase residual and separate the non-isotropic tropospheric delay correction from the residual;

[0014] A ground station is randomly selected, and the GNSS continuous observation data of the reference station is used to determine the precise location of the ground station through the SINEX file to obtain the PPP ionosphere-free combined observation value.

[0015] The obtained PPP ionosphere-free combined observation value is subtracted from the pseudorange and the known value to obtain the phase residual, and the residual of the PPP ionosphere-free combined observation value at the ground station is denoised using the EDM model;

[0016] A sample satellite is randomly selected, a sampling time is preset, and a satellite whose altitude angle differs from that of the sample satellite by no more than 1° is selected during the sampling time as a satellite at the same altitude angle as the sample satellite;

[0017] The phase residuals of the two selected satellites at the same altitude angle are subtracted to construct the residual difference; since the two satellites are at the same altitude angle, after the two satellites are subtracted, only the difference in the tropospheric delay mapping function model error caused by the different azimuth angles remains, that is, the difference in the non-isotropic tropospheric delay correction values ​​on the paths of the two satellites;

[0018] Step 2. Calculate the non-isotropic tropospheric delay correction using the least squares method;

[0019] The formula of the difference of the residuals obtained in step 1 is formed into a linear system of equations;

[0020] Determine whether the coefficient matrix of the linear equation system is rank-deficient. If it is rank-deficient, change the observation station. If it is full-rank, use the least square method to calculate the linear equation and obtain the optimal solution of the corresponding azimuth non-isotropic tropospheric delay correction number.

[0021] Step 3. Calculate the corrected tropospheric slant delay;

[0022] The optimal solution of the non-isotropic tropospheric delay correction calculated in step 2 is compared with the tropospheric slant delay calculated by the mapping function method. Add them together to get the corrected tropospheric slant delay.

[0023] The present invention has the following advantages:

[0024] As described above, the present invention relates to a method for correcting tropospheric delays that takes into account non-isotropic tropospheric delays. The method first uses the GNSS continuous observation data of the reference station to determine the precise position of the ground station through the SINEX file to obtain the PPP ionosphere-free combined observation value, and then subtracts the observation value from the pseudorange and the known value (satellite clock error, satellite orbit error, etc.) to obtain the residual (here it is considered to be the non-isotropic tropospheric delay correction value), and then uses the EDM model to denoise the residual of the PPP ionosphere-free combined observation value at the station, and then randomly selects a sample satellite, and presets Sampling time, select a satellite whose altitude angle differs from the sample satellite by no more than 1° within the sampling time, make a difference between the phase residuals of the two selected satellites, construct the difference of the residuals, and obtain a linear equation system. Then introduce a linear equation system solving algorithm based on least squares to calculate the non-isotropic tropospheric delay correction value, find a solution that minimizes the error, deal with the situation where the equation is redundant or cannot be solved accurately, so that the process reduces the influence of wild values ​​compared to the neural network, and finally add the obtained non-isotropic tropospheric delay correction value to the SPD calculated by the mapping function to obtain the corrected SPD. This method breaks through the constraints of isotropy and anisotropy of tropospheric delay, introduces the theory of non-isotropy, and improves the prediction accuracy of SPD by improving the calculation method of SPD. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 is a flow chart of a tropospheric delay correction method taking into account non-isotropy in an embodiment of the present invention;

[0026] Figure 2 This is a schematic diagram of the distribution of IGS stations in Europe;

[0027] Figure 3 Schematic diagram of the non-isotropic tropospheric delay correction value of the BRUX station on June 1, 2023. Figure 3 (a), (b), (c), and (d) show the solution results of 0-6 hours, 7-12 hours, 13-18 hours, and 18-24 hours respectively;

[0028] Figure 4 RMS histogram of the European reference station on 1 June 2023 to take into account the non-isotropic tropospheric delay correction;

[0029] Figure 5 STD histogram of the European reference stations on 1 June 2023, corrected for the non-isotropic tropospheric delay. DETAILED DESCRIPTION

[0030] Glossary:

[0031] The non-isotropy of tropospheric delay means that at a certain altitude angle, the tropospheric delays at certain azimuths can be considered to be approximately equal, but at other azimuths the tropospheric delays have large differences.

[0032] The non-isotropic tropospheric delay correction is the difference between the true value of SPD (i.e., the tropospheric slant delay calculated by the ray tracing method) and the tropospheric slant delay calculated by the mapping function method.

[0033] The present invention is further described in detail below with reference to the accompanying drawings and specific embodiments:

[0034] This embodiment proposes a tropospheric delay correction method that takes non-isotropy into account. This method addresses the problem that most traditional methods are constrained to isotropy or anisotropy when obtaining SPD, thereby limiting the prediction accuracy of SPD. This method innovatively introduces the theory of non-isotropy and improves the prediction accuracy of SPD by improving the calculation method of SPD.

[0035] The method of the present invention starts with the observation data between the station and the satellite that can be obtained, separates the unmodeled error and the modeled error from the residual through the PPP positioning algorithm, and then solves the SPD correction value to obtain the approximate true value of the SPD. The unmodeled error here is, for example, white noise, and the modeled error refers to the non-isotropic tropospheric delay correction number.

[0036] like Figure 1 As shown, the tropospheric delay correction method taking into account the non-isotropic nature of the present embodiment comprises the following steps:

[0037] Step 1. Calculate the phase residual and separate the non-isotropic tropospheric delay correction from the residual.

[0038] A ground station is randomly selected, and the GNSS continuous observation data of the European reference station is used to determine the precise position of the station through the SINEX file to obtain the PPP ionosphere-free combined observation value.

[0039] The observed value is subtracted from the pseudorange and known values ​​(satellite clock error, satellite orbit error, etc.) to obtain the phase residual.

[0040] like Figure 1 As shown in the blue box on the left, the calculation process of the phase residual is as follows:

[0041] First, the GNSS carrier ionosphere-free combined observation values ​​are listed through precise ephemeris and combined with SINEX station location information. The equation is shown in formula (1):

[0042] (1)

[0043] In the formula, is the geometric distance between the station and the satellite, s is the satellite, c is the speed of light in vacuum, and are the receiver and satellite clock errors, represents the tropospheric delay, is the wavelength, N is the carrier integer ambiguity, is the remaining error term, including the tropospheric delay mapping function model error and the remaining error terms that cannot be modeled , as shown in formula (2):

[0044] (2)

[0045] In the non-differential positioning data processing, precise ephemeris is used to determine the satellite orbital position, and precise clock error products are used to eliminate the error effect of the satellite clock. The present invention introduces the continuous observation data (including post-orbital error) and precise clock error provided by the international GNSS service to eliminate the satellite clock error. The error equation of the combined carrier phase without ionospheric error and orbit error can be rewritten as:

[0046] (3)

[0047] In the formula, Represents the residual, that is, the error term in formula (1) .

[0048] The EDM model is further used to denoise the residuals of the PPP ionosphere-free combined observations at the station to weaken the correlation between the various errors in the residuals. Here, the EDM model is relatively conventional and will not be described in detail.

[0049] A sample satellite is randomly selected, and the sampling time is τ seconds. During the sampling time, a satellite whose altitude angle differs from that of the sample satellite by no more than 1° is selected as the satellite at the same altitude angle as the sample satellite.

[0050] For example, the value of τ is 30 seconds. Of course, in this embodiment, τ is not limited to the above 30 seconds.

[0051] The same altitude angle here is just an approximation that when the altitude angle difference does not exceed 1°, the two are at the same altitude angle.

[0052] The phase residuals of the two satellites at the same selected altitude angle are subtracted to construct the residual difference.

[0053] Since the satellites have the same elevation angle, the If the two satellites are the same, then the residual can be considered to contain only the error term caused by ignoring the non-isotropic characteristics of the tropospheric delay. Since the elevation angles of the two satellites are approximately the same, the only remaining error after the two satellites are subtracted is the difference in the tropospheric delay mapping function model error caused by the different azimuths, that is, the difference in the non-isotropic tropospheric delay correction values ​​on the paths of the two satellites.

[0054] Take τ seconds as the sampling time, select a satellite whose altitude angle does not differ from the sample satellite by more than 1° within the sampling time, make a difference between the phase residuals of the two selected satellites, and construct the residual difference. The calculation formula is as follows: As shown:

[0055] (4)

[0056] In the formula, represents the difference of residuals, , Respectively represent the residuals of the i-th and j-th satellites, , They represent the tropospheric delay mapping function model errors of the i-th and j-th satellites, , They represent the remaining error terms of the i-th and j-th satellites that cannot be modeled, respectively, i, j = 1, 2, ..., n, i and j are not equal, and n represents the number of satellites with the same altitude angle.

[0057] Due to the error that cannot be modeled is very small, so it is ignored. Since the SPD calculated by the mapping function for the satellites is the same when the altitude angle is the same, it is considered that the residual approximately only contains the error term caused by ignoring the non-isotropic characteristics of the tropospheric delay, that is, the difference between the non-isotropic tropospheric delay correction values ​​on the paths of the two satellites, as shown in the formula and As shown;

[0058] (5)

[0059] Right now:

[0060] (6)

[0061] In the formula, , are the non-isotropic tropospheric delay correction values ​​of the i-th and j-th satellites, respectively. It represents the difference between the non-isotropic tropospheric delay correction values ​​of the i-th and j-th satellites.

[0062] The proposal of non-isotropy of tropospheric delay is based on the assumption that the tropospheric delay of the mapping function is isotropic. The non-isotropy of tropospheric delay refers to the change of SPD caused by the change of azimuth at the same altitude.

[0063] Since the SPD calculated by the mapping function is isotropic, that is, the variation of SPD with azimuth is ignored, the delay correction value of the tropospheric delay taking into account the non-isotropic angle is defined as the true value of SPD at the same altitude angle and The difference between .

[0064] The calculation process of the non-isotropic tropospheric delay correction is as follows:

[0065] Firstly, the tropospheric slant delay is calculated by using the ray tracing method and the mapping function method, and the calculated results are recorded as SPD and ; Among them, the SPD calculated by the ray tracing method is used as the reference value.

[0066] Then the non-isotropic tropospheric delay correction value is calculated, which is the reference value SPD of the tropospheric oblique delay at the same altitude angle and The difference between them is shown in formula (7):

[0067] (7)

[0068] in, is the non-isotropic tropospheric delay correction value; SPD is the reference value, which is calculated by the ray tracing method; is the tropospheric slant delay calculated by the mapping function method.

[0069] Step 2. Calculate the non-isotropic tropospheric delay correction using the least squares method.

[0070] The difference in the non-isotropy values ​​between two different satellites in formula (6) is regarded as an equation in which two unknowns are subtracted to obtain a known number, that is, the non-isotropy value of each satellite at different time points is regarded as an unknown number.

[0071] At the same time, ignoring the short-term changes of tropospheric delay, assuming that the non-isotropic value of a certain satellite remains unchanged within six hours, within this time interval, using the satellites observed by the selected station, the satellites that meet the condition that the elevation angle of all satellites selected within the 30s sampling time differs from the elevation angle of the sample satellite by no more than 1°, and the phase residual of its observation data is subtracted from the phase residual of the random sample satellite observation data to construct a linear equation system, which is expressed as:

[0072] (8)

[0073] In the formula, represents the non-isotropic tropospheric delay correction value of the first satellite, represents the non-isotropic tropospheric delay correction value of the i-th satellite, It represents the difference between the non-isotropic tropospheric delay correction values ​​of the i-th and j-th satellites.

[0074] In formula (8), the azimuth angles of all satellite positions are Rounding to the nearest integer, and is an integer, and equation (8) is written as The form is ; where A is the coefficient matrix.

[0075] For example, if the azimuth angle of the i-th satellite is 181.3 degrees, then its azimuth angle is approximately 181 degrees. .

[0076] If there are enough satellites with approximately the same elevation angle, the coefficient matrix A of the listed equations will have a full rank property, and the linear equations can be solved to obtain the optimal solution of the non-isotropic tropospheric delay correction value corresponding to the azimuth angle. , and the solution is It has a unique solution; if the rank is deficient, change the ground station and repeat steps 1-2 above.

[0077] Step 3. Calculate the corrected tropospheric slant delay.

[0078] The optimal solution of the non-isotropic tropospheric delay correction calculated in step 2 is compared with the tropospheric slant delay calculated by the mapping function method. Add them together to get the corrected tropospheric slant delay.

[0079] Taking the i-th satellite as an example, the non-isotropic tropospheric delay correction value of the i-th satellite will be obtained Add , that is, the corrected tropospheric slant delay of the i-th satellite is obtained.

[0080] Next, the SPD at the satellite orbit observed by the station is calculated by using the satellite positions in the GNSS observation data and the NWM data through a ray tracing algorithm, and the accuracy of the method of the present invention is evaluated using the standard deviation and mean error.

[0081] (9)

[0082] (10)

[0083] Among them, SPD R To use the SPD data of the Radiate reference, consider it as the true value, SPD M is the SPD data calculated using the tropospheric delay correction method taking into account the non-isotropy proposed in the present invention, and N is the sample size.

[0084] At the same time, in the process of solving the non-isotropic value of tropospheric delay, the resolution of the altitude angle is set to 5°, and the same solution process is repeated in each interval. Since the non-isotropy of altitude angles above 50° has little effect on the SPD value, and the quality of satellite observation data below 10° is low, a cutoff altitude angle of 10° is set in the PPP solution.

[0085] Therefore, this experiment only considers the tropospheric delay with an altitude angle range of 10°-50° and an azimuth angle range of 0°-360°.

[0086] Taking Europe as the research area, the non-isotropic correction method based on linear equations is used to model this area. Figure 2 The IGS stations shown in the figure use equations (9) and (10) as error evaluation indicators.

[0087] Taking the BRUX station as an example, the station-based non-isotropic tropospheric delay correction value was constructed.

[0088] The present invention uses the GNSS observation data of June 1, 2023 as a sample, with a solution interval of 6 hours, and the altitude and azimuth resolutions are 1° and 5° respectively. The non-isotropic tropospheric delay correction value of that day is as follows: Figure 3 shown.

[0089] Figure 3 In the figure, the polar axis has an altitude range of 10°-50°, an azimuth range of 0°-360°, and the color of the scattered points is the optimal solution for the non-isotropic value at that location, with a maximum value of 0.4m and a minimum value of -0.4m. Figure 3(a), (b), (c), and (d) show the solution results for 0-6 o'clock, 7-12 o'clock, 13-18 o'clock, and 18-24 o'clock, respectively.

[0090] The method of the present invention only corrects the non-isotropic value at the satellite trajectory. The blank area in the figure indicates that no satellite passed on that day. The results show that the non-isotropic value of tropospheric delay is highly correlated with the altitude angle, and there are also fluctuations in the azimuth angle domain.

[0091] 17 stations were selected as reference stations in Europe to verify the SPD correction accuracy of the proposed method on June 1, 2023. The STD and RMS histograms at each station are as follows Figure 4 and Figure 5 shown.

[0092] Figure 4 and Figure 5 In the figure, the horizontal axis is the measuring station, and the vertical axis represents the RMS and STD values ​​(in meters). The colors of the bar graph represent the results of different time periods: purple represents UTC 0-6, blue represents UTC 6-12, orange represents UTC 12-18, and green represents UTC 18-24.

[0093] from Figure 4 and Figure 5 It can be seen that compared with the SPD calculated by Radiate, the internal and external accuracy of the present invention are both 6cm to 14cm, and the external accuracy is better overall. In order to further compare the accuracy differences between different stations, the average RMS and average STD values ​​of each IGS station on the day were counted, as shown in Table 1 below.

[0094] Table 1 Statistics of average RMS and average STD of European reference stations on June 1, 2023

[0095]

[0096] As can be seen from Table 1, the average RMS and average STD values ​​of each station are distributed between 0.07 m and 0.12 m. Among them, the RMS and STD values ​​of stations bor1 and ganp are relatively low, which are 0.0769 m, 0.0755 m and 0.0761 m, 0.0733 m, respectively, showing a higher solution accuracy. Relatively speaking, the error of graz station is the highest, with RMS and STD values ​​of 0.1170 m and 0.1060 m, respectively, indicating that the solution results of this station fluctuate greatly. The error values ​​of the remaining stations are between the above extreme values.

[0097] Overall, the average STD value of each station is slightly lower than the RMS value, reflecting the good stability of the solution results and the high internal accuracy. However, there are still certain accuracy differences between different stations, which may be related to the geographical location of the station, the observation environment and the quality of the observation data of the station. But in general, the present invention has greatly improved the accuracy of SPD.

[0098] Compared with the traditional method of predicting and obtaining SPD, the method of the present invention first uses the GNSS continuous observation data of the reference station, determines the precise position of the station through the SINEX file, obtains the PPP ionosphere-free combined observation value, and subtracts the observation value from the pseudorange and the known value (satellite clock error, satellite orbit error, etc.) to obtain the phase residual (here it is considered to be the non-isotropic tropospheric delay correction value), then uses the EDM model to denoise the residual of the PPP ionosphere-free combined observation value at the station, and then randomly selects a sample satellite, takes τ seconds as the sampling time, selects a satellite whose altitude angle does not differ from that of the sample satellite by no more than 1° within the sampling time, subtracts the phase residuals of the two selected satellites, constructs the residual difference, and obtains a linear equation group, and then introduces a linear equation group solving algorithm based on the least squares method to model the difference of the non-isotropic value of SPD, reduces the influence of wild values ​​compared with the neural network, and adds the obtained non-isotropic tropospheric delay correction value to the SPD calculated by the mapping function to obtain the corrected SPD. The method of the present invention breaks through the constraints of isotropy and anisotropy of tropospheric delay, proposes a theory of non-isotropy of tropospheric delay, and improves the calculation accuracy of SPD.

[0099] Of course, the above description is only a preferred embodiment of the present invention, and the present invention is not limited to the above embodiments. It should be noted that all equivalent substitutions and obvious deformation forms made by any technician familiar with the field under the guidance of this specification fall within the essential scope of this specification and should be protected by the present invention.

Claims

1. A method for correcting tropospheric delay taking into account non-isotropy, characterized in that: The steps include: Step 1. Randomly select a ground station, use the GNSS continuous observation data of the reference station, determine the precise location of the ground station through the SINEX file, and obtain the PPP ionosphere-free combined observation value; The obtained PPP ionosphere-free combined observation value is subtracted from the pseudorange and the known value to obtain the phase residual, and the residual of the PPP ionosphere-free combined observation value at the ground station is denoised using the EDM model; A sample satellite is randomly selected and a sampling time is preset. During the sampling time, a satellite whose altitude angle differs from that of the sample satellite by no more than 1° is selected as a satellite at the same altitude angle as the sample satellite. The phase residuals of the two selected satellites at the same altitude angle are subtracted to construct the residual difference; since the two satellites are at the same altitude angle, the only difference left after the two satellites are subtracted is the difference in the tropospheric delay mapping function model error caused by the different azimuth angles, that is, the difference in the non-isotropic tropospheric delay correction values ​​on the paths of the two satellites; Step 2. The formula of the difference of the residuals obtained in step 1 is formed into a linear system of equations; Determine whether the coefficient matrix of the linear equation system is rank-deficient. If it is rank-deficient, change the observation station. If it is full-rank, use the least square method to calculate the linear equation and obtain the optimal solution of the corresponding azimuth non-isotropic tropospheric delay correction number. Step 3. Add the optimal solution of the non-isotropic tropospheric delay correction calculated in step 2 to the tropospheric slant delay calculated by the mapping function method to obtain the corrected tropospheric slant delay.

2. The method for correcting tropospheric delay taking into account non-isotropy according to claim 1, characterized in that: In step 1, the calculation process of the phase residual is as follows: First, the GNSS carrier ionosphere-free combined observation values ​​are listed through precise ephemeris and combined with SINEX station location information. The equation is shown in formula (1): （1） In the formula, is the geometric distance between the station and the satellite, s is the satellite, c is the speed of light in vacuum, and are the receiver and satellite clock errors, represents the tropospheric delay, is the wavelength, N is the carrier integer ambiguity, is the remaining error term, including the tropospheric delay mapping function model error and the remaining error terms that cannot be modeled , as shown in formula (2): （2） Eliminate satellite clock errors by introducing the post-orbit errors and precise clock errors provided by international GNSS services and orbit error, the error equation of the ionosphere-free combined carrier phase can be rewritten as: （3） In the formula, Represents the phase residual, that is, the error term in formula (1) .

3. The method for correcting tropospheric delay taking into account non-isotropy according to claim 1, characterized in that: In step 1, the calculation formula of the residual difference is as follows: As shown: （4） In the formula, represents the difference of residuals, , Respectively represent the residuals of the i-th and j-th satellites, , They represent the tropospheric delay mapping function model errors of the i-th and j-th satellites, , They represent the remaining error terms that cannot be modeled for the i-th and j-th satellites respectively, i, j = 1, 2, ..., n, i and j are not equal, and n represents the number of satellites with the same altitude angle.

4. The method for correcting tropospheric delay taking into account non-isotropy according to claim 3, characterized in that: The error term that cannot be modeled , Ignore it, and because the SPD calculated by the mapping function for the satellites is the same when the altitude angle is the same, it is considered that the residual only contains the error term caused by ignoring the non-isotropic characteristics of the tropospheric delay, that is, the non-isotropic tropospheric delay correction value on the respective paths of the two satellites. The calculation formula is as follows: and As shown; （5） Right now: （6） In the formula, , are the non-isotropic tropospheric delay correction values ​​of the i-th and j-th satellites, respectively. It represents the difference between the non-isotropic tropospheric delay correction values ​​of the i-th and j-th satellites.

5. The method for correcting tropospheric delay taking into account non-isotropy according to claim 4, characterized in that: In step 1, the calculation process of the non-isotropic tropospheric delay correction value is as follows: Firstly, the tropospheric slant delay is calculated by using the ray tracing method and the mapping function method, and the calculated results are recorded as SPD and ; Among them, the SPD calculated by ray tracing method is used as the reference value; Then the non-isotropic tropospheric delay correction value is calculated, which is the reference value SPD of the tropospheric oblique delay at the same altitude angle and The difference between them is shown in formula (7): （7） in, is the non-isotropic tropospheric delay correction value; SPD is the reference value, which is calculated by the ray tracing method; is the tropospheric slant delay calculated by the mapping function method.

6. The method for correcting tropospheric delay taking into account non-isotropy according to claim 4, characterized in that: The step 2 is specifically as follows: The difference in the non-isotropy values ​​between two different satellites in formula (6) is regarded as an equation in which two unknown numbers are subtracted to obtain a known number, that is, the non-isotropy value of each satellite at different time points is regarded as an unknown number; At the same time, ignoring the short-term changes of tropospheric delay, assuming that the non-isotropic value of a certain satellite remains unchanged within six hours, within this time interval, using the satellites observed by the selected station, the satellites that meet the condition that the difference between the elevation angle of all satellites and the elevation angle of the sample satellite within the 30s sampling time is no more than 1° are selected, and the phase residual of their observation data is subtracted from the phase residual of the random sample satellite observation data, and a linear equation system is constructed, such as the formula As shown: （8） In the formula, represents the non-isotropic tropospheric delay correction value of the first satellite, represents the non-isotropic tropospheric delay correction value of the i-th satellite, represents the difference between the non-isotropic tropospheric delay correction values ​​of the i-th and j-th satellites; In formula (8), the azimuth angles of all satellite positions are , rounded to the nearest integer, then and is an integer, and equation (8) is written as In the form of, A represents the coefficient matrix; ； If the linear equations The coefficient matrix A has a full rank characteristic, then the linear equations are solved to obtain the optimal solution of the non-isotropic tropospheric delay correction value corresponding to the azimuth angle. , and the solution is has a unique solution; If the rank is deficient, change the ground station and repeat steps 1-2.

Images