Troposphere chromatography method, troposphere chromatography system and troposphere chromatography equipment based on GNSS (Global Navigation Satellite System) and medium

Abstract

The invention discloses a troposphere chromatography method, a troposphere chromatography system, troposphere chromatography equipment and a medium based on a GNSS (Global Navigation Satellite System). The troposphere chromatography method comprises the following steps of: dividing a three-dimensional chromatography region into a plurality of continuous voxel blocks in horizontal and vertical three-dimensional directions; obtaining a GNSS error term according to the divided voxel blocks, and discretizing troposphere wet delay in the GNSS error term into an integral of a product of a wet refractive index parameter to be estimated and a signal length; constructing an observation equation, and substituting the obtained integral into a GNSS original pseudo-range and phase observation equation; in consideration of GNSS station network distribution and signal distribution characteristics, a numerical constraint or empirical constraint condition is introduced to compensate matrix rank deficiency; and solving by adopting a PPP technology to obtain a wet refractive index parameter, and converting the wet refractive index parameter into water vapor density for meteorological analysis and forecasting. According to the invention, the precision and reliability of the inversion result of the three-dimensional wet refractive index field under the extreme weather condition are improved.

Description

Technical Field

[0001] The present invention belongs to the technical field of atmospheric sounding, and particularly relates to a GNSS-based tropospheric tomography method, system, device and medium. Background Technique

[0002] Atmospheric water vapor is the hub of the global land-atmosphere water cycle. Its content and distribution changes directly affect global weather and climate change, especially playing a key role in the generation and evolution of extreme weather. Atmospheric water vapor is mainly concentrated in the troposphere. Traditional methods for detecting atmospheric water vapor include radiosondes, water vapor microwave radiometers, spaceborne remote sensing for earth observation, etc. With the development of satellite navigation technology, the method of using the Global Navigation Satellite System (GNSS) to remotely measure water vapor content has gradually emerged due to its advantages such as global coverage, all-weather, high precision, low cost, and high spatio-temporal resolution. When the satellite signal passes through the troposphere, it will produce a delay due to the refraction effect. On the one hand, it is necessary to correct the tropospheric delay error in satellite navigation positioning to improve the accuracy. On the other hand, the GNSS observation can be used to obtain the atmospheric water vapor distribution information through inversion.

[0003] Traditional ground-based GNSS technology relies on ground-based GNSS stations and can only obtain the total water vapor content in the zenith direction of the troposphere above the station or interpolate to obtain the two-dimensional plane. It lacks the water vapor content distribution information in the vertical direction. In the field of geophysics, tomography technology is an important inversion method, aiming to use ray signals and computed tomography (CT) technology to inversely reconstruct the material density information in the study area. The GNSS tropospheric tomography applies CT technology and uses the slant path tropospheric wet delay (SWD) or slant path water vapor content (SWV) suffered by the GNSS signal as the signal to achieve the inversion of high-precision three-dimensional atmospheric wet refractivity or water vapor density distribution information, providing data for synoptic research and weather forecasting.

[0004] The existing tropospheric tomography methods mainly include two steps: 1) First, use the Network Solution or Precise Point Positioning (PPP) method to process the original carrier and pseudorange observations received by the GNSS station network in the study area, obtain the zenith total delay (ZTD) of the stations, and combine the temperature and pressure parameters obtained by the meteorological equipment equipped at the stations to jointly recover the SWD or SWV using the wet projection function; 2) Use the SWD or SWV obtained by each station for tomography to inversely obtain the three-dimensional wet refractivity (Nw) or water vapor density (WVD) distribution in the study area. This tomography method mainly uses the projection function (referring to the wet projection function) to convert the zenith delay to the slant path direction and further establish the relationship between the SWD / SWV observations and the tomography parameters (wet refractivity or water vapor density). However, the modeling process of the empirical function uses long-term data and parameterizes the conversion relationship from the zenith to the slant path direction. This process results in the empirical function being unable to achieve a good conversion under extreme weather conditions, thus introducing errors into the SWD / SWV recovery process. This situation has a greater impact in extreme weather conditions with more significant atmospheric asymmetry and will further affect the accuracy of the tomography results. In addition, the current two-step tomography method requires GNSS data processing first and then tomography operations in terms of process, which is limited in the operationalization of the meteorological field mainly for real-time applications. Therefore, the present invention proposes a GNSS-based tropospheric tomography method to solve the problems that the wet projection function affects the tomography accuracy under extreme weather conditions and the two-step tomography method is not applicable in real-time scenarios. Summary of the Invention

[0005] To solve the above technical problems, the present invention proposes a GNSS-based tropospheric tomography method, system, device, and medium, which directly discretizes the tropospheric slant path wet delay SWD into the integral along the signal path of the product of the tomography parameter (wet refractivity Nw) and the signal path length (i.e., tomography parameterization), substitutes the parameterized wet refractivity parameter into the GNSS pseudorange / phase observation equation, and uses the precise point positioning PPP method for parameter solution. By directly constructing the functional relationship between the tomography parameter (wet refractivity) and the GNSS original observations, the direct solution of the tomography parameter is realized during the GNSS data processing. This is not only simpler in process and beneficial to real-time meteorological applications, but also avoids the introduction of modeling errors under extreme weather conditions caused by the wet projection function, improving the accuracy and reliability of the inversion results of the three-dimensional wet refractivity field under extreme weather conditions.

[0006] On the one hand, to achieve the above object, the present invention provides a GNSS-based tropospheric tomography method, including:

[0007] Dividing the three-dimensional tomography area into several continuous voxel blocks in the horizontal and vertical three-dimensional directions;

[0008] Obtaining the GNSS error term according to the divided voxel blocks, and discretizing the tropospheric wet delay in the GNSS error term into an integral of the product of the wet refractivity parameter to be estimated and the signal length;

[0009] Constructing an observation equation and substituting the obtained integral into the GNSS raw pseudorange and phase observation equations;

[0010] Considering the GNSS station network distribution and signal distribution characteristics, introducing numerical constraints or empirical constraint conditions to compensate for the matrix rank deficiency;

[0011] Using the PPP technology to solve for the wet refractivity parameter, and converting the wet refractivity parameter into water vapor density for meteorological analysis and forecasting.

[0012] Optionally, dividing the three-dimensional tomography area into several continuous voxel blocks in the horizontal and vertical three-dimensional directions includes: dividing in the horizontal direction by longitude and latitude, and dividing in the vertical direction in a uniform manner or in a lower-dense and upper-sparse manner according to the change of water vapor density.

[0013] Optionally, obtaining the GNSS error term according to the divided voxel blocks includes:

[0014]

[0015] where the coordinates in the Earth-Centered Earth-Fixed coordinate system are respectively (x 1 , y 1 , z 1 ) and (x 2 , y 2 , z 2 ); ds is the GNSS error term.

[0016] Optionally, discretizing the tropospheric wet delay in the GNSS error term into an integral of the product of the wet refractivity parameter to be estimated and the signal length includes:

[0017]

[0018] where and ds i are respectively the wet refractivity within the i-th voxel block and the intercept length of the signal passing through the voxel block; SWD is the integral of the product of the wet refractivity parameter to be estimated and the signal length.

[0019] Optionally, constructing the observation equation includes:

[0020] A PPP model using two types of observations, namely carrier phase and pseudorange, of Beidou / GNSS;

[0021] For the total tropospheric delay along the slant path, it is divided into slant path hydrostatic delay and slant path wet delay according to the generation reasons and characteristics. The slant path hydrostatic delay includes zenith hydrostatic delay and dry mapping function;

[0022] Expand the slant path wet delay into the form of an integral of the product of the wet refractivity parameter to be estimated and the signal length, and complete the construction of the observation equation.

[0023] Optionally, the zenith hydrostatic delay is calculated by using the Saastamoinen, Hopfield, and Black empirical models in combination with barometric pressure and station information.

[0024] Optionally, the numerical constraint is constrained by external observation methods, including using lidar, satellite remote sensing, radiosonde, and ground meteorological stations to constrain the wet refractivity of the voxels at the position to be measured.

[0025] Optionally, the empirical constraint refers to a constraint mode established according to the distribution characteristics of water vapor density or wet refractivity. The empirical constraint includes horizontal constraint, vertical constraint, and layer top constraint.

[0026] On the other hand, to achieve the above object, the present invention also provides a GNSS-based tropospheric tomography system, including:

[0027] Voxel block division module, discretization integration module, observation equation construction module, constraint introduction module, and wet refractivity parameter conversion module;

[0028] The voxel block division module is used to divide the three-dimensional tomography area into several continuous voxel blocks in the horizontal and vertical three-dimensional directions;

[0029] The discretization integration module is used to obtain the GNSS error term according to the divided voxel blocks, and discretize the tropospheric wet delay in the GNSS error term into the integral of the product of the wet refractivity parameter to be estimated and the signal length;

[0030] The observation equation construction module is used to construct the observation equation and substitute the obtained integral into the GNSS original pseudorange and phase observation equations;

[0031] The constraint introduction module is used to consider the rank deficiency of the compensation matrix due to the GNSS station network distribution and signal distribution characteristics, and introduce numerical constraints or empirical constraint conditions;

[0032] The wet refractivity parameter conversion module is used to solve the wet refractivity parameter by using PPP technology and convert the wet refractivity parameter into water vapor density for meteorological analysis and forecasting.

[0033] A computer device, comprising: at least one processor; and a memory storing a computer program that can run on the processor, wherein the processor, when executing the program, performs the steps of the method.

[0034] A computer-readable storage medium storing a computer program, wherein the computer program, when executed by a processor, performs the steps of the method.

[0035] Technical effects of the present invention: The present invention discloses a GNSS-based tropospheric tomography method, system, device and medium, which directly constructs a functional relationship between tomography parameters and GNSS raw observations and solves them during GNSS data processing. The present invention simplifies the data processing flow and has the ability to perform real-time processing using filtering, making it more suitable for real-time meteorological monitoring services. The present invention does not rely on wet projection functions and directly establishes a functional relationship between tomography parameters and GNSS raw observations, avoiding the inapplicability of wet projection functions under extreme weather conditions and thus the influence of errors introduced during the conversion of zenith delay to the slant path direction. On the one hand, the present invention improves the accuracy and reliability of the inversion results, and on the other hand, it better adapts to different meteorological conditions and improves the all-weather ability of GNSS tropospheric tomography technology. Description of the Drawings

[0036] The drawings forming a part of this application are used to provide a further understanding of this application. The schematic embodiments and descriptions thereof of this application are used to explain this application and do not constitute an improper limitation of this application. In the drawings:

[0037] Figure 1 It is a schematic flow chart of a GNSS-based tropospheric tomography method according to an embodiment of the present invention;

[0038] Figure 2 It is a schematic structural diagram of a GNSS-based tropospheric tomography system according to an embodiment of the present invention. Detailed Embodiments

[0039] It should be noted that, without conflict, the embodiments in this application and the features in the embodiments can be combined with each other. The following will refer to the drawings and combine the embodiments to detail this application.

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

[0041] Such asFigure 1 As shown in the figure, in this embodiment, a GNSS-based tropospheric tomography method is provided, including:

[0042] Dividing the three-dimensional tomography area into several continuous voxel blocks in the horizontal and vertical three-dimensional directions;

[0043] Obtaining GNSS error terms according to the divided voxel blocks, and discretizing the tropospheric wet delay in the GNSS error terms into an integral of the product of the wet refractivity parameter to be estimated and the signal length;

[0044] Constructing an observation equation and substituting the obtained integral into the GNSS raw pseudorange and phase observation equations;

[0045] Considering the distribution characteristics of the GNSS station network and signals, introducing numerical constraints or empirical constraint conditions to compensate for the matrix rank deficiency;

[0046] Using the PPP technology to solve for the wet refractivity parameter, and converting the wet refractivity parameter into water vapor density for meteorological analysis and forecasting.

[0047] Furthermore, dividing the three-dimensional tomography area into several continuous voxel blocks in the horizontal and vertical three-dimensional directions includes: dividing in the horizontal direction along the longitude and latitude directions, and dividing in the vertical direction in a uniform manner or in a way of being denser at the bottom and sparser at the top according to the change of water vapor density.

[0048] Furthermore, obtaining GNSS error terms according to the divided voxel blocks includes:

[0049]

[0050] Among them, the coordinates in the Earth-centered Earth-fixed coordinate system are respectively (x 1 , y 1 , z 1 ) and (x 2 , y 2 , z 2 ); ds is the GNSS error term.

[0051] Furthermore, discretizing the tropospheric wet delay in the GNSS error terms into an integral of the product of the wet refractivity parameter to be estimated and the signal length includes:

[0052]

[0053] Among them, and ds i are respectively the wet refractivity within the i-th voxel block and the intercept length of the signal passing through the voxel block; SWD is the integral of the product of the wet refractivity parameter to be estimated and the signal length.

[0054] Furthermore, constructing an observation equation includes:

[0055] A PPP model using two types of observations, namely carrier phase and pseudorange, of Beidou / GNSS;

[0056] For the total tropospheric delay along the slant path, it is divided into slant path hydrostatic delay and slant path wet delay according to the generation causes and characteristics, where the slant path hydrostatic delay includes zenith hydrostatic delay and dry projection function;

[0057] Expand the slant path wet delay into the form of an integral of the product of the wet refractivity parameter to be estimated and the signal length to complete the construction of the observation equation.

[0058] Furthermore, the zenith hydrostatic delay is calculated by combining the Saastamoinen, Hopfield, and Black empirical models with barometric pressure and station information.

[0059] Furthermore, the numerical constraint is constrained by external observation methods, including using lidar, satellite remote sensing, radiosonde, and ground meteorological stations to constrain the wet refractivity of the voxels at the position to be measured.

[0060] Furthermore, the empirical constraint refers to the constraint mode established according to the distribution characteristics of water vapor density or wet refractivity, and the empirical constraint includes horizontal constraint, vertical constraint, and layer top constraint.

[0061] As Figure 1 shown, the specific implementation steps of the present invention are as follows:

[0062] S1. Divide the three-dimensional tomography area into several continuous voxel blocks in the horizontal and vertical three-dimensional directions. The sizes of the voxel blocks can be equal or unequal. Optionally, in the horizontal direction, it is usually divided in the longitude and latitude directions, and in the vertical direction, it is generally divided in a uniform manner or in a lower-dense and upper-sparse manner according to the change of water vapor density.

[0063] S2. Under the condition of ignoring the influence of the earth's curvature, the signal path from the GNSS satellite to the antenna position of the station can be regarded as a straight line in space. Establish a straight line equation in space and find the coordinates of the intersection points of the straight line and each horizontal and vertical stratified plane, that is, the coordinates of the piercing points.

[0064] S3. After finding the coordinates of the piercing points on all stratified planes, calculate the intercept length ds of the GNSS signal path in each voxel block; taking a voxel block on a non-bottommost plane as an example, a certain GNSS signal passing through this voxel block will generate 2 piercing points, and their coordinates in the earth-centered earth-fixed coordinate system are (x 1 , y 1 , z 1 ) and (x 2 , y 2 , z 2) Then the intercept length ds of the signal in this voxel block is expressed as:

[0065]

[0066] Among them, the coordinates in the Earth-centered Earth-fixed coordinate system are (x 1 , y 1 , z 1 ) and (x 2 , y 2 , z 2 ); ds is the GNSS error term.

[0067] S4. Discretize the tropospheric wet delay SWD into the integral of the product of the wet refractivity parameter to be estimated and the signal length.

[0068] S4-1. There are mainly two existing discretization methods, namely the classical model proposed by Flores et al. and the interpolation integral model based on the Newton-Cotes formula proposed by Perler et al.

[0069] S4-2. The classical model assumes that the wet refractivity within the divided voxel block is uniform and invariant. The slant-path wet delay within a certain voxel block can be directly expressed as the product of the wet refractivity N w and the intercept ds; while the interpolation integral model takes into account the non-uniform distribution characteristic of water vapor in space. Through the wet refractivities at the 8 node positions of this voxel block, the horizontal bilinear interpolation method and the Newton-Cotes formula are used to characterize the slant-path wet delay within this voxel block. Compared with the interpolation integral model, the classical model involves relatively fewer parameters and the calculation operation is simpler and more efficient; however, the interpolation model takes into account the non-uniform distribution characteristic of water vapor and is closer to the real situation.

[0070] S4-3. The present invention takes the classical model as an example for description. Specifically, it is assumed that the wet refractivity N w within a certain voxel block is uniformly distributed and invariant during the tomography period. The slant-path wet delay can be expressed as the product of the wet refractivity N w and the intercept length ds; and for the entire signal path, the total slant-path wet delay SWD of this signal can be expressed as the integral form shown in the following formula:

[0071]

[0072] Among them, and ds i are respectively the wet refractivity within the i-th voxel block and the intercept length of the signal passing through the voxel block; SWD is the integral of the product of the wet refractivity parameter to be estimated and the signal length.

[0073] S5. Construct the observation equation:

[0074] S5-1. For the PPP model using the carrier phase and pseudorange observations of Beidou / GNSS, the pseudorange \(P\) and carrier phase \(L\) observation equations are generally expressed as follows:

[0075]

[0076] where the subscripts \(r\) and \(j\) represent the receiver and frequency number respectively, and the superscripts \(s\) and \(Q\) represent the satellite and satellite system respectively; represents the geometric distance between the satellite antenna and the receiver antenna, in meters; \(c\) represents the speed of light in vacuum, in meters per second; \(t\) r and \(t\) s,Q represent the receiver and satellite clock biases respectively, in seconds; \(T\) r is the total tropospheric delay (Slant Total Delay, STD) on the signal path between the receiver and the satellite, in meters; \(\lambda\) j represents the wavelength of frequency \(j\); represents the ambiguity of frequency \(j\); \(\gamma\) j is the ionospheric delay amplification factor for the L1 frequency observation value; is the slant ionospheric delay of the first frequency pseudorange observation value, in meters; and represent the pseudorange hardware delays of frequency \(j\) at the receiver end and satellite end respectively, in meters; and represent the carrier phase hardware delays of frequency \(j\) at the receiver end and satellite end respectively, in meters; and are the sums of the observation noise, multipath effects and other unmodeled errors of the pseudorange and carrier phase respectively, in meters.

[0077] S5-2. For the total slant path tropospheric delay STD, it can be divided into the slant hydrostatic delay (SHD) and the slant wet delay SWD according to the generation reasons and characteristics, as follows:

[0078] \(T\) r = SHD + SWD;

[0079] where the slant hydrostatic delay SHD can be expressed by the zenith hydrostatic delay (ZHD) and the dry projection function,

[0080] SHD = mf d (E)×ZHD;

[0081] where \(mf\) d represents the dry projection function, and its value is related to the satellite elevation angle \(E\). ​

[0082] S5-3. S5-3, ZHD can be calculated by combining the Saastamoinen, Hopfield, and Black empirical models with barometric pressure and station information.

[0083] S5-4. S5-4. Expand the slant-path wet delay SWD into the integral form described in S4, then the observation equation is rewritten as:

[0084]

[0085] The above observation equation does not contain the wet projection function, avoiding the errors introduced by its modeling process and improving the accuracy of the tomography results.

[0086] S6. S6. Use the constraint equation to compensate for the rank deficiency of the matrix. Due to the GNSS station network distribution and the characteristics of the GNSS signal geometric configuration, there are some low-level voxels without signal penetration. Empirical constraint conditions or numerical constraint conditions can be added to ensure the accuracy and stability of the overall tomography parameter solution.

[0087] S6-1. S6-1. Empirical constraints refer to the constraint mode established based on the distribution characteristics of water vapor density or wet refractivity, mainly including horizontal constraints, vertical constraints, and top-of-layer constraints.

[0088] S6-2. S6-2. For the numerical constraint method, external observation methods can be used for constraint, such as using lidar, satellite remote sensing, radiosonde, and ground meteorological stations to constrain the wet refractivity of specific-position voxels to enhance the well-posedness and solution accuracy of the overall parameter equation.

[0089] S7. S7. After correcting the ionospheric delay, antenna error, tidal error, and hydrostatic tropospheric delay error terms in the GNSS observations, use the precise point positioning PPP method to estimate and solve the tomography parameters (wet refractivity) and other GNSS parameters (receiver clock error, receiver coordinates, carrier phase ambiguity, inter-system bias).

[0090] S8. S8. According to actual needs, convert the wet refractivity into water vapor density for meteorological analysis and forecasting.

[0091] As Figure 2 shown, the present invention also provides a GNSS-based tropospheric tomography system, including: a voxel block division module, a discretized integration module, an observation equation construction module, a constraint introduction module, and a wet refractivity parameter conversion module;

[0092] The voxel block division module is used to divide the three-dimensional tomography area into a number of continuous voxel blocks in the horizontal and vertical three-dimensional directions;

[0093] The discretized integration module is used to obtain GNSS error terms based on the divided voxel blocks, and discretize the tropospheric wet delay in the GNSS error terms into an integral of the product of the wet refractivity parameter to be estimated and the signal length;

[0094] The observation equation construction module is used to construct an observation equation and substitute the obtained integral into the GNSS raw pseudorange and phase observation equations;

[0095] The constraint introduction module is used to introduce numerical constraints or empirical constraint conditions to compensate for the rank deficiency of the matrix in consideration of the GNSS network distribution and the distribution characteristics of signals;

[0096] The wet refractivity parameter conversion module is used to solve the wet refractivity parameter by using PPP technology and convert the wet refractivity parameter into water vapor density for meteorological analysis and forecasting.

[0097] A computer device includes: at least one processor; and a memory storing a computer program that can run on the processor, wherein the processor executes the steps of the method when executing the program.

[0098] A computer-readable storage medium stores a computer program, wherein the computer program executes the steps of the method when executed by a processor.

[0099] The present invention discloses a GNSS-based tropospheric tomography method, system, device and medium, which directly constructs a functional relationship between tomography parameters and GNSS raw observations and solves them during GNSS data processing. The present invention simplifies the data processing flow and has the ability to perform real-time processing using filtering, making it more suitable for real-time meteorological monitoring services. The present invention does not rely on wet projection functions and directly establishes a functional relationship between tomography parameters and GNSS raw observations, avoiding the influence of errors introduced during the conversion of zenith delay to the slant path direction due to the inapplicability of wet projection functions under extreme weather conditions. On the one hand, the present invention can improve the accuracy and reliability of the inversion results, and on the other hand, it can better adapt to different meteorological conditions and improve the all-weather ability of GNSS tropospheric tomography technology.

[0100] In view of the problems of the traditional two-step GNSS tropospheric tomography method, such as relying on empirical wet projection functions and complex calculation processes, the present invention proposes a one-step tropospheric tomography method with strong applicability, high precision, and being simpler and more effective. This method directly constructs the functional relationship between tomography parameters and GNSS raw observations, and directly solves the tomography parameters during the GNSS data processing, avoiding the possibility of introducing errors in the process of converting the zenith delay to the slant path direction by the wet projection function. On the one hand, this method solves the problem of large tomography errors caused by the inapplicability of the wet projection function in extreme weather conditions, improving the applicability, precision, and reliability of the model; on the other hand, the data processing flow is more simplified, with the ability to perform real-time processing using filtering, and it is more suitable for real-time meteorological monitoring services.

[0101] The above are only the preferred specific embodiments of this application, but the protection scope of this application is not limited thereto. Any changes or substitutions that can be easily thought of by those skilled in the art within the technical scope disclosed in this application should be covered by the protection scope of this application. Therefore, the protection scope of this application should be subject to the protection scope of the claims.

Claims

1. A tropospheric tomography method based on GNSS, characterized in that: include: Divide the three-dimensional tomographic region into a number of continuous voxel blocks in the horizontal and vertical three-dimensional directions; Obtaining a GNSS error term according to the divided voxel blocks, and discretizing the tropospheric wet delay in the GNSS error term into an integral of the product of the wet refractive index parameter to be estimated and the signal length; Construct the observation equation and bring the obtained integral into the GNSS original pseudorange and phase observation equation; Considering the distribution characteristics of GNSS station network and signals, numerical constraints or empirical constraints are introduced to compensate for the matrix rank deficiency. The wet refractive index parameter is obtained by using the PPP technology, and the wet refractive index parameter is converted into water vapor density for meteorological analysis and forecasting.

2. The GNSS-based tropospheric tomography method according to claim 1, characterized in that: The three-dimensional tomography area is divided into a number of continuous voxel blocks in the horizontal and vertical three-dimensional directions, including: the horizontal direction is divided in the longitude and latitude directions, and the vertical direction is divided in a uniform manner or in a sparse manner according to the change of water vapor density.

3. The GNSS-based tropospheric tomography method according to claim 1, characterized in that: The GNSS error terms obtained according to the divided voxel blocks include: Among them, the coordinates in the Earth-centered Earth-fixed coordinate system are (x1, y1, z1) and (x2, y2, z2) respectively; ds is the GNSS error term.

4. The GNSS-based tropospheric tomography method according to claim 1, characterized in that: The process of discretizing the tropospheric wet delay in the GNSS error term into the product of the wet refractive index parameter to be estimated and the signal length includes: in, and ds i are the wet refractive index within the i-th voxel block and the intercept length of the signal passing through the voxel block; SWD is the integral of the product of the wet refractive index parameter to be estimated and the signal length.

5. The GNSS-based tropospheric tomography method according to claim 1, characterized in that: Constructing the observation equation includes: The PPP model uses two observation quantities: carrier phase and pseudorange of Beidou / GNSS; The total tropospheric delay of the slant path is divided into the slant path static delay and the slant path wet delay according to the causes and characteristics. The slant path static delay includes the zenith static delay and the dry projection function. The oblique path wet delay is expanded into the form of an integral of the product of the wet refractive index parameter to be estimated and the signal length, thereby completing the construction of the observation equation.

6. The GNSS-based tropospheric tomography method according to claim 1, characterized in that: The numerical constraints are constrained by external observation methods, including using laser radar, satellite remote sensing, sounding, and ground meteorological stations to constrain the wet refractive index of the voxels at the location to be measured.

7. The GNSS-based tropospheric tomography method according to claim 1, characterized in that: The empirical constraints refer to constraint patterns established according to water vapor density or wet refractive index distribution characteristics, and the empirical constraints include horizontal constraints, vertical constraints, and layer top constraints.

8. A system for the GNSS-based tropospheric tomography method according to any one of claims 1 to 7, characterized in that: include: Voxel block partitioning module, discretization integration module, observation equation construction module, constraint introduction module and wet refractive index parameter conversion module; The voxel block division module is used to divide the three-dimensional tomography area into a plurality of continuous voxel blocks in the horizontal and vertical three-dimensional directions; The discretization integration module is used to obtain the GNSS error term according to the divided voxel blocks, and discretize the tropospheric wet delay in the GNSS error term into the integral of the product of the wet refractive index parameter to be estimated and the signal length; The observation equation construction module is used to construct the observation equation and bring the obtained integral into the GNSS original pseudorange and phase observation equation; The constraint introduction module is used to introduce numerical constraints or empirical constraints to compensate for the matrix rank deficiency in consideration of the distribution characteristics of the GNSS station network and the signal; The wet refractive index parameter conversion module is used to adopt the PPP technology to solve the wet refractive index parameter, and convert the wet refractive index parameter into water vapor density for meteorological analysis and forecasting.

9. A computer device comprising: at least one processor; and a memory storing a computer program executable on the processor, wherein the processor executes the steps of the method according to any one of claims 1 to 7 when executing the program.

10. A computer-readable storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the steps of the method according to any one of claims 1 to 7 are performed.

Images