A GNSS water vapor tomography method based on dynamic variable scale voxels

Abstract

The present application relates to the technical field of GNSS meteorology, and particularly relates to a GNSS water vapor tomography method based on dynamic variable scale voxels. The present application comprises the following steps: S10: collecting GNSS data and preprocessing the same to obtain a scanning area of GNSS water vapor signals and GNSS signals; S20: determining a dynamic variable scale tomography boundary corresponding to the GNSS signals by using a boundary optimization algorithm; S30: determining a dynamic variable scale voxel block according to a water vapor factor invariable rule based on water vapor vertical profile information provided by sounding data; S40: constructing a GNSS water vapor tomography observation equation and a constraint equation according to a spatial position relationship between the GNSS water vapor signals and the dynamic variable scale voxel block; S50: solving the tomography equation set by using an optimization algebra reconstruction algorithm, and performing accuracy evaluation on the solving result. Through the above method, the structural stability of the GNSS water vapor tomography model is effectively improved, and the inversion accuracy of a three-dimensional atmospheric water vapor density real-time field is improved.

Description

Technical Field

[0001] This invention relates to the field of GNSS meteorology technology, specifically to a GNSS water vapor tomography method based on dynamic variable-scale voxels. Background Technology

[0002] Atmospheric water vapor is one of the most abundant and important greenhouse gases in the Earth system. Its rapid spatial and temporal variations play a crucial role in the global water cycle and energy balance, and it is a key meteorological element in scientific research such as extreme weather forecasting and climate change. In recent years, meteorological natural disasters have shown a severe trend of increasing frequency, intensity, and duration, causing adverse effects on people's lives and property, socio-economic development, and the ecological environment. There is an urgent need to develop high-precision, high-spatial-temporal, and all-weather refined atmospheric water vapor detection technologies to meet the major needs of my country's construction of an integrated air-ground-human coordinated atmospheric detection technology system, enhance the ability to prevent extreme weather events, and serve important areas such as national meteorological early warning and disaster prevention and mitigation.

[0003] GNSS water vapor tomography, as an emerging atmospheric water vapor detection technology, has rapidly become a highly promising high-spatiotemporal resolution atmospheric water vapor detection method due to its advantages such as high observation accuracy, high spatiotemporal resolution, all-weather monitoring, and low operating costs. However, the mismatch between the "box-shaped" tomographic region and the "inverted cone-shaped" spatial geometry of the GNSS signal cluster in the GNSS water vapor tomography model leads to a large number of voxels in the three-dimensional discretized tomographic structure being unable to be penetrated by GNSS signal lines. This causes ill-posedness in the water vapor tomography equations, severely limiting the inversion accuracy of the tomographic results, which is the key challenge of GNSS water vapor tomography. Many researchers have adopted different methods to address this key problem: on the one hand, by fusing multi-source water vapor data such as radiosonde data, standard atmospheric models, radio occultation data, atmospheric reanalysis data, and atmospheric infrared detectors, the quality of the initial values ​​for water vapor tomography solutions can be improved, thereby enhancing the inversion accuracy of GNSS water vapor tomography. Furthermore, some scholars have addressed the structural defects of GNSS signals by adding multi-system GNSS signals, side observation signals, virtual observation signals, InSAR water vapor signals, and remote sensing water vapor signals to GNSS water vapor tomography models. This significantly increases the number of penetrating voxels in the three-dimensional tomography model, thereby improving the quality of water vapor tomography results. However, the inconsistency between the spatiotemporal resolution of multi-type water vapor data and the spatiotemporal resolution of GNSS data remains a bottleneck limiting the practical application of these studies. On the other hand, some scholars have optimized the spatial discretization method of GNSS water vapor tomography models. For example, they have designed the optimal voxel horizontal resolution, constructed a dual-resolution discretization method, and appropriately increased the size of upper-level voxels. By adjusting the position and size of voxels, they aim to ensure that as many voxels as possible are penetrating by GNSS signals, thereby alleviating the ill-posedness of GNSS water vapor tomography models. However, the voxels in these methods are all of fixed scale, and the fixed boundaries of the three-dimensional tomographic region are determined empirically and discretized accordingly, lacking rigorous scientific discretization criteria. More importantly, atmospheric water vapor is not uniformly distributed in physical space, but exhibits a significant vertical decrease. Existing technologies do not consider the important impact of this variation on the voxels of the tomographic model, and cannot meet the requirements for inverting a high-precision three-dimensional water vapor density field. Summary of the Invention

[0004] To address the aforementioned issues, this invention provides a GNSS water vapor tomography method based on dynamically variable-scale voxels. This method fully considers the spatial dynamic variations and inverted conical geometry of GNSS signals. Based on the dynamic scanning area of ​​the GNSS signal, a boundary optimization algorithm is used to determine the optimal dynamically variable-scale tomography boundary, resolving the spatial geometric mismatch between the GNSS signal and the fixed tomography region. Furthermore, by combining the vertical decrease in atmospheric water vapor content with the water vapor factor invariance criterion, the optimal size of the dynamically variable-scale voxels within different altitude layers is scientifically determined. This ensures a high degree of matching between the GNSS water vapor tomography model framework and the actual distribution of atmospheric water vapor, thereby achieving high-precision inversion of the three-dimensional water vapor density field.

[0005] To achieve the above-mentioned objectives, the present invention adopts the following technical solution:

[0006] A GNSS water vapor chromatography method based on dynamic variable-scale voxels includes the following steps:

[0007] S10: Collect observation data from GNSS stations within the study area, process the GNSS data using high-precision data processing software, and then determine the oblique path water vapor observation value of the GNSS water vapor signal and the scanning area of ​​the GNSS signal.

[0008] S20: The tropospheric study area is discretized into multiple vertical stratified intervals in the vertical direction. Based on the scanning area of ​​the GNSS signal in each vertical stratified interval, the dynamic variable-scale tomographic boundary of each height layer is determined by the boundary optimization algorithm.

[0009] S30: Based on the sounding water vapor profile data of the first three days of the tomographic epoch, the exponential interpolation algorithm is used to vertically interpolate it to the height of each layer interval to obtain the average water vapor density information of different layer intervals. Then, the variable scale voxel block information of each layer interval is determined by using the water vapor factor invariance criterion.

[0010] S40: Based on the intercept information of GNSS signals in variable-scale voxel blocks, the tomographic observation equations of GNSS signals are constructed. At the same time, based on the spatial positional relationship between dynamic variable-scale voxels within the entire tomographic model framework, horizontal constraint equations and vertical constraint equations are established, thereby forming a set of water vapor tomographic equations based on dynamic variable-scale voxels.

[0011] S50: The optimized algebraic reconstruction algorithm is used to solve the water vapor tomography equations based on dynamic variable-scale voxels to obtain the water vapor tomography results. The accuracy of the tomography results is evaluated using synchronous sounding water vapor profile data. If the accuracy index meets the requirements, the water vapor tomography results can express the actual three-dimensional water vapor density field in the study area. Otherwise, steps S20 to S50 are repeated to recalculate the water vapor tomography results.

[0012] Furthermore, as a preferred embodiment of the present invention, step S10 includes the following steps:

[0013] S11: Process the raw GNSS observation data using the professional software GAMIT to obtain tropospheric delay information for each GNSS station. This tropospheric delay information includes the zenith tropospheric delay (ZTD) of the GNSS station and the east-west atmospheric horizontal gradient (G). EW and the north-south atmospheric horizontal gradient G NS ;

[0014] S12: Remove the tropospheric dry delay ZHD at the GNSS station from the ZTD to obtain the tropospheric wet delay ZWD. The tropospheric dry delay ZHD can be calculated using formula (1).

[0015]

[0016] Among them, P S This indicates the surface air pressure of the GNSS station. H and H represent the latitude and elevation of the GNSS station, respectively;

[0017] S13: Calculate the east-west horizontal wet gradient of the GNSS station. and the horizontal moisture gradient in the north-south direction From G respectively EW and G NS Extract it, and then use formula (2) to calculate the slant path water vapor content (SWV) of the GNSS signal;

[0018]

[0019] Where Π represents the conversion factor, mf w and mf g Let e ​​and α represent the wet projection function and the gradient projection function, respectively, and let e and α represent the elevation angle and azimuth angle of the GNSS signal, respectively.

[0020] S14: Determine the scanning area of ​​the GNSS signal in the three-dimensional troposphere based on the elevation and azimuth information of the GNSS signal.

[0021] As a further preferred embodiment of the present invention, the plurality of vertical stratification intervals in S20 are 15 stratification intervals, and the heights of each stratification interval from the Earth's surface to the top of the troposphere are 0.4km, 0.4km, 0.4km, 0.4km, 0.4km, 0.4km, 0.4km, 0.4km, 1.0km, 1.0km, 1.5km, 1.5km and 2.0km respectively.

[0022] As a further preferred technical solution of the present invention, the boundary optimization algorithm in step S20 is as follows: First, determine the puncture point coordinate information of each GNSS signal in the vertical layering interval, and then obtain the smallest polygon containing all puncture points. Second, using each side of the smallest polygon as the side length, construct rectangles containing the smallest polygon in sequence to form multiple rectangles. Finally, compare the area information of these rectangles and select the longest deformation with the smallest area as the optimal tomographic region boundary of the layering interval. Repeat this algorithm for different layering intervals in sequence to determine the dynamic variable-scale tomographic boundary of the entire three-dimensional tomographic model.

[0023] Furthermore, as a preferred embodiment of the present invention, in step S30, the exponential interpolation algorithm is as shown in formula (3):

[0024]

[0025] Where ρ(P0) represents the water vapor density at the interpolation point P0, ρ(P1) and ρ(P2) represent the water vapor density at the observation points P1 above and P2 below P0, respectively, corresponding to h1 and h2 representing the height difference between the interpolation point P0 and the two observation points P1 and P2, respectively, and Δh representing the vertical distance between the two observation points. scale Indicates the water vapor elevation value;

[0026] The constant water vapor factor criterion is shown in formula (4):

[0027]

[0028] Among them, s i and s j These represent the horizontal side lengths of the voxel blocks in the i-th and j-th layers, respectively. Correspondingly, hei i andhei j Then, ρ(radiosonde) represents the height of these two vertically layered intervals, respectively. i ) and ρ(radiosonde j ) represents the sounding water vapor density value at these two heights.

[0029] Furthermore, as a preferred embodiment of the present invention, in step S40, the tomographic observation equation of the GNSS signal is as shown in formula (5):

[0030]

[0031] Among them, SWV q This represents the oblique path water vapor content observation value of the q-th GNSS signal. This represents the puncture length of the signal in the variable-scale voxel block at the k-th layer, i-th row, and j-th column. If the signal does not puncture the voxel block, then... =0, ρ(voxel ijk ) represents the water vapor density value of the voxel block, and n, m, and t represent the total number of rows, columns, and layers of the GNSS water vapor tomography model based on dynamic variable-scale voxels.

[0032] The horizontal constraint equation is shown in equation (6):

[0033]

[0034] Wherein, ρ(voxel center ) represents the water vapor density value of the voxel block to be constrained within the vertical stratification interval, ρ(surround) g ω represents the water vapor density value of the g-th voxel surrounding the given voxel within the same stratified interval. g The weight information corresponding to the g-th voxel block can be calculated using formula (7), where l g σ represents the distance from the g-th voxel block to the voxel block to be constrained, and σ is the smoothing factor;

[0035]

[0036] The vertical constraint equation is shown in equation (8):

[0037]

[0038] Wherein, ρ(voxel k ) and ρ(voxel k+1 ) represent the water vapor density values ​​of the k-th and k+1-th voxel blocks, respectively, h k and h k+1 These represent the heights of the two layered intervals, respectively;

[0039] The water vapor chromatography equation set based on dynamic variable-scale voxels is shown in equation (9):

[0040]

[0041] Among them, A GNSS A H A V Let X represent the tomographic model coefficient matrices corresponding to the GNSS signal, horizontal constraint, and vertical constraint, respectively, and let X represent the parameter vector to be determined, which consists of the unknown water vapor density of all dynamically variable scale voxels.

[0042] Furthermore, as a preferred embodiment of the present invention, in step S50, the optimized algebraic reconstruction algorithm is as shown in formula (10):

[0043]

[0044] in, Let λ represent the water vapor density value of the j-th voxel during the k-th iteration, and let λ represent the relaxation factor, which plays a role in regulating convergence during the iteration process. ij Y represents the element in the i-th row and j-th column of the coefficient matrix of the tomographic model. i G represents the element in the i-th row of the column vector of the tomographic model. i The weight of the i-th observation signal can be calculated using formula (11), where e i Indicates the elevation angle of the i-th GNSS signal;

[0045] G i =sin 2 (e i (11)

[0046] The accuracy index is the root mean square error (RMSE), which can be calculated using formula (12):

[0047]

[0048] in, This represents the water vapor density value of the i-th voxel calculated by the chromatography model. This represents the water vapor density value of the same voxel provided by the radiosonde station, and r represents the total number of voxels.

[0049] Furthermore, as a preferred embodiment of the present invention, for the horizontal and vertical constraint observation equations, the corresponding weight information G i All are 1.

[0050] Furthermore, as a preferred embodiment of the present invention, when the root mean square error (RMSE) is less than 1.5 g / m 3 At that time, the accuracy of the water vapor chromatography results met the requirements.

[0051] The GNSS water vapor chromatography method based on dynamically variable scale voxels described in this invention has the following technical advantages compared with existing technologies:

[0052] This invention provides a GNSS water vapor tomography method based on dynamically variable-scale voxels. It fully considers the spatial dynamic variation characteristics and inverted conical geometry of GNSS signals. Based on the dynamic scanning area of ​​the GNSS signal, a boundary optimization algorithm is used to determine the optimal dynamically variable-scale tomography boundary, solving the problem of spatial geometric mismatch between the GNSS signal and the fixed tomography region. Combining the vertical decrease law of atmospheric water vapor content, the optimal size of the dynamically variable-scale voxels in different altitude layers is scientifically determined using the water vapor factor invariance criterion. This ensures that the GNSS water vapor tomography model framework highly matches the actual distribution of atmospheric water vapor, thereby achieving high-precision inversion of the three-dimensional water vapor density field. Attached Figure Description

[0053] Figure 1 The diagram shown is a flowchart of a GNSS water vapor chromatography method based on dynamically variable scale voxels according to an embodiment of the present invention.

[0054] Figure 2 The diagram shown is a simplified flowchart of a GNSS water vapor chromatography method based on dynamically variable scale voxels according to an embodiment of the present invention.

[0055] Figure 3 The diagram shown is a schematic diagram illustrating the principle of constant water vapor factor according to an embodiment of the present invention.

[0056] Figure 4 The figure shown is a schematic diagram of a GNSS water vapor tomography model based on dynamically scaled voxels according to an embodiment of the present invention. Detailed Implementation

[0057] The present invention will be further explained in detail below with reference to the accompanying drawings, so that those skilled in the art can better understand and implement the present invention. However, the following examples are only used to explain the present invention and are not intended to limit the present invention.

[0058] This embodiment provides a GNSS water vapor chromatography method based on dynamically scaled voxels, such as... Figures 1-2 As shown, it includes the following steps:

[0059] S10: Collect observation data from GNSS stations within the study area, process the GNSS data using high-precision data processing software, and then determine the oblique path water vapor observation value of the GNSS water vapor signal and the scanning area of ​​the GNSS signal.

[0060] S20: The tropospheric study area is discretized vertically into multiple vertical stratified intervals. Based on the scanning area of ​​the GNSS signal in each vertical stratified interval, a boundary optimization algorithm is used to determine the dynamic variable-scale tomographic boundary for each height layer. The multiple vertical stratified intervals consist of 15 stratified intervals, with heights from the Earth's surface to the top of the troposphere of each interval being 0.4 km, 0.4 km, 0.4 km, 0.4 km, 0.4 km, 0.4 km, 0.4 km, 1.0 km, 1.0 km, 1.5 km, 1.5 km, and 2.0 km.

[0061] S30: Based on the sounding water vapor profile data of the first three days of the tomographic epoch, the exponential interpolation algorithm is used to vertically interpolate it to the height of each layer interval to obtain the average water vapor density information of different layer intervals. Then, the variable scale voxel block information of each layer interval is determined by using the water vapor factor invariance criterion.

[0062] S40: Based on the intercept information of GNSS signals in variable-scale voxel blocks, the tomographic observation equations of GNSS signals are constructed. At the same time, based on the spatial positional relationship between dynamic variable-scale voxels within the entire tomographic model framework, horizontal constraint equations and vertical constraint equations are established, thereby forming a set of water vapor tomographic equations based on dynamic variable-scale voxels.

[0063] S50: The optimized algebraic reconstruction algorithm is used to solve the water vapor tomography equations based on dynamic variable-scale voxels to obtain the water vapor tomography results. The accuracy of the tomography results is evaluated using synchronous radiosonde water vapor profile data. If the accuracy index meets the requirements, the water vapor tomography results can express the actual three-dimensional water vapor density field in the study area. Otherwise, steps S20, S30, S40, and S50 are repeated to recalculate the water vapor tomography results.

[0064] Step S10 includes the following steps: S11: Process the raw GNSS observation data using the professional software GAMIT to obtain the tropospheric delay information of each GNSS station, wherein the tropospheric delay information includes the zenith tropospheric delay ZTD of the GNSS station and the east-west atmospheric horizontal gradient G. EW and the north-south atmospheric horizontal gradient G NS ;

[0065] S12: Remove the tropospheric dry delay ZHD at the GNSS station from the ZTD to obtain the tropospheric wet delay ZWD. The tropospheric dry delay ZHD can be calculated using formula (1).

[0066]

[0067] Among them, P S This indicates the surface air pressure of the GNSS station. H and H represent the latitude and elevation of the GNSS station, respectively;

[0068] S13: Calculate the east-west horizontal wet gradient of the GNSS station. and the horizontal moisture gradient in the north-south direction From G respectively EW and G NS Extract it, and then use formula (2) to calculate the slant path water vapor content (SWV) of the GNSS signal;

[0069]

[0070] Where Π represents the conversion factor, mf w and mf g Let e ​​and α represent the wet projection function and the gradient projection function, respectively, and let e and α represent the elevation angle and azimuth angle of the GNSS signal, respectively.

[0071] S14: Determine the scanning area of ​​the GNSS signal in the three-dimensional troposphere based on the elevation and azimuth information of the GNSS signal.

[0072] In step S20, the boundary optimization algorithm is as follows: First, determine the coordinate information of the puncture point of each GNSS signal in the vertical layering interval, and then obtain the smallest polygon containing all puncture points. Second, using each side of the smallest polygon as the side length, construct rectangles containing the smallest polygon in sequence to form multiple rectangles. Finally, compare the area information of these rectangles and select the longest deformation with the smallest area as the optimal tomographic region boundary of the layering interval. Repeat this algorithm for different layering intervals in sequence to determine the dynamic variable-scale tomographic boundary of the entire three-dimensional tomographic model.

[0073] like Figure 3 The diagram illustrates the principle of constant water vapor factor. Based on the vertical profile data of the radiosonde, an exponential interpolation algorithm is used to determine the average water vapor density at different altitudes. This principle is then used to accurately determine the size information of variable-scale voxels within different strata. In step S30, the exponential interpolation algorithm is shown in formula (3):

[0074]

[0075] Where ρ(P0) represents the water vapor density at the interpolation point P0, ρ(P1) and ρ(P2) represent the water vapor density at the observation points P1 above and P2 below P0, respectively, corresponding to h1 and h2 representing the height difference between the interpolation point P0 and the two observation points P1 and P2, respectively, and Δh representing the vertical distance between the two observation points. scale Indicates the water vapor elevation value;

[0076] The constant water vapor factor criterion is shown in formula (4):

[0077]

[0078] Among them, s i and s j These represent the horizontal side lengths of the voxel blocks in the i-th and j-th layers, respectively. Correspondingly, hei i andhei j Then, ρ(radiosonde) represents the height of these two vertically layered intervals, respectively. i ) and ρ(radiosonde j ) represents the sounding water vapor density value at these two heights.

[0079] like Figure 4 The diagram shows a schematic of a GNSS water vapor tomography model based on dynamically variable-scale voxels. Based on the spatial relationship between the GNSS signal and the variable-scale voxels, a tomographic observation equation for the GNSS signal is established. Simultaneously, horizontal and vertical constraint equations between the dynamically variable-scale voxels are established, thus forming a set of water vapor tomography equations based on dynamically variable-scale voxels. In step S40, the tomographic observation equation for the GNSS signal is shown in formula (5):

[0080]

[0081] Among them, SWV q This represents the oblique path water vapor content observation value of the q-th GNSS signal. This represents the puncture length of the signal in the variable-scale voxel block at the k-th layer, i-th row, and j-th column. If the signal does not puncture the voxel block, then... =0, ρ(voxel ijk ) represents the water vapor density value of the voxel block, and n, m, and t represent the total number of rows, columns, and layers of the GNSS water vapor tomography model based on dynamic variable-scale voxels.

[0082] The horizontal constraint equation is shown in equation (6):

[0083]

[0084] Wherein, ρ(voxel center ) represents the water vapor density value of the voxel block to be constrained within the vertical stratification interval, ρ(surround) g ω represents the water vapor density value of the g-th voxel surrounding the given voxel within the same stratified interval. g The weight information corresponding to the g-th voxel block can be calculated using formula (7), where l gσ represents the distance from the g-th voxel block to the voxel block to be constrained, and σ is the smoothing factor;

[0085]

[0086] The vertical constraint equation is shown in equation (8):

[0087]

[0088] Wherein, ρ(voxel k ) and ρ(voxel k+1 ) represent the water vapor density values ​​of the k-th and k+1-th voxel blocks, respectively, h k and h k+1 These represent the heights of the two layered intervals, respectively;

[0089] The water vapor chromatography equation set based on dynamic variable-scale voxels is shown in equation (9):

[0090]

[0091] Among them, A GNSS A H A V Let X represent the tomographic model coefficient matrices corresponding to the GNSS signal, horizontal constraint, and vertical constraint, respectively, and let X represent the parameter vector to be determined, which consists of the unknown water vapor density of all dynamically variable scale voxels.

[0092] In step S50, the optimized algebraic reconstruction algorithm is as shown in formula (10):

[0093]

[0094] in, Let λ represent the water vapor density value of the j-th voxel during the k-th iteration, and let λ represent the relaxation factor, which plays a role in regulating convergence during the iteration process. ij Y represents the element in the i-th row and j-th column of the coefficient matrix of the tomographic model. i G represents the element in the i-th row of the column vector of the tomographic model. i The weight of the i-th observation signal can be calculated using formula (11), where e i Let G represent the elevation angle of the i-th GNSS signal. For the horizontal and vertical constrained observation equations, the corresponding weight information G... i All are 1;

[0095] G i =sin 2 (e i (11)

[0096] The accuracy index is the root mean square error (RMSE), which can be calculated using formula (12):

[0097]

[0098] in, This represents the water vapor density value of the i-th voxel calculated by the chromatography model. This represents the water vapor density value of the same voxel provided by the radiosonde station, where r represents the total number of voxels. When the root mean square error (RMSE) is less than 1.5 g / m³... 3 At that time, the accuracy of the water vapor chromatography results met the requirements.

[0099] The specific implementation schemes described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific implementation schemes of the present invention and are not intended to limit the scope of the present invention. Any equivalent changes and modifications made by those skilled in the art without departing from the concept and principles of the present invention should fall within the scope of protection of the present invention.

Claims

1. A GNSS water vapor chromatography method based on dynamically scaled voxels, characterized in that, Includes the following steps: S10: Collect observation data from GNSS stations within the tropospheric study area, process the observation data using high-precision data processing software, and then determine the oblique path water vapor observation value of the GNSS water vapor signal and the scanning area of ​​the GNSS signal; S20: The tropospheric study area is discretized into multiple vertical stratified intervals in the vertical direction. Based on the scanning area of ​​the GNSS signal in each vertical stratified interval, the dynamic variable-scale tomographic boundary of each vertical stratified interval is determined by the boundary optimization algorithm. S30: Based on the sounding water vapor profile data of the first three days of the tomographic epoch, the exponential interpolation algorithm is used to vertically interpolate it to the height of each layer interval to obtain the average water vapor density information of different layer intervals. Then, the variable scale voxel block information of each layer interval is determined by using the water vapor factor invariance criterion. S40: Based on the intercept information of GNSS signals in variable-scale voxel block information, construct the tomographic observation equation of GNSS signals. At the same time, establish horizontal constraint equations and vertical constraint equations according to the spatial positional relationship between dynamic variable-scale voxels in the entire tomographic model framework, and then build a water vapor tomographic equation set based on dynamic variable-scale voxels. S50: The optimized algebraic reconstruction algorithm is used to solve the water vapor tomography equations based on dynamic variable-scale voxels to obtain the water vapor tomography results. The accuracy of the tomography results is evaluated using synchronous sounding water vapor profile data. If the accuracy index meets the requirements, the water vapor tomography results can express the actual three-dimensional water vapor density field in the study area. Otherwise, steps S20 to S50 are repeated to recalculate the water vapor tomography results.

2. The GNSS water vapor chromatography method based on dynamically variable-scale voxels according to claim 1, characterized in that, Step S10 includes the following steps: S11: Process the raw GNSS observation data using the professional software GAMIT to obtain tropospheric delay information for each GNSS station. This tropospheric delay information includes the zenith tropospheric delay (ZTD) of the GNSS station and the east-west atmospheric horizontal gradient. and the horizontal atmospheric gradient in the north-south direction ; S12: Remove the tropospheric dry delay ZHD at the GNSS station from the zenith tropospheric delay ZTD to obtain the tropospheric wet delay ZWD. The tropospheric dry delay ZHD is calculated using formula (1). （1） in, This indicates the surface air pressure of the GNSS station. and These represent the latitude and elevation of the GNSS station, respectively. S13: Calculate the east-west horizontal wet gradient of the GNSS station. and the horizontal moisture gradient in the north-south direction From respectively and Extract it, and then use formula (2) to calculate the slant path water vapor content (SWV) of the GNSS signal; （2） in, Indicates the conversion factor. and Let these represent the wet projection function and the gradient projection function, respectively. and These represent the elevation angle and azimuth angle of the GNSS signal, respectively. S14: Determine the scanning area of ​​the GNSS signal in the three-dimensional troposphere based on the elevation and azimuth information of the GNSS signal.

3. The GNSS water vapor chromatography method based on dynamically variable-scale voxels according to claim 1, characterized in that, The S20 consists of 15 vertical stratification intervals, with heights of 0.4 km, 0.4 km, 0.4 km, 0.4 km, 0.4 km, 0.4 km, 0.4 km, 0.4 km, 0.4 km, 0.4 km, 1.0 km, 1.0 km, 1.5 km, 1.5 km and 2.0 km from the Earth's surface to the top of the troposphere.

4. The GNSS water vapor chromatography method based on dynamically variable-scale voxels according to claim 1, characterized in that, The boundary optimization algorithm in step S20 is as follows: First, determine the puncture point coordinates of each GNSS signal within the vertical layering interval, and then obtain the smallest polygon containing all puncture points. Next, using each side of the smallest polygon as the side length, construct rectangles containing the smallest polygon in sequence to form multiple rectangles. Finally, compare the area information of these rectangles and select the rectangle with the smallest area as the optimal tomographic region boundary of the layering interval. Repeat this algorithm for different layering intervals in sequence to determine the dynamic variable-scale tomographic boundary of the entire three-dimensional tomographic model.

5. The GNSS water vapor chromatography method based on dynamically variable-scale voxels according to claim 1, characterized in that, In step S30, the exponential interpolation algorithm is as shown in formula (3): （3） in, Indicates the point to be interpolated The water vapor density value at that location, and Then they respectively represent The observation point above and the observation point below The water vapor density value, correspondingly, and These represent the points to be interpolated. To two observation points and The height difference, This represents the vertical distance between two observation points. Indicates the water vapor elevation value; The constant water vapor factor criterion is shown in formula (4): （4） in, and These represent the horizontal side lengths of the voxel blocks in the i-th and j-th layers, respectively. and These represent the heights of the two vertically separated sections, and This represents the sounding water vapor density values ​​at these two altitudes.

6. The GNSS water vapor chromatography method based on dynamically variable-scale voxels according to claim 1, characterized in that, In step S50, the optimized algebraic reconstruction algorithm is as shown in formula (10): （10） in, This represents the water vapor density value of the j-th voxel during the (k+1)-th iteration. This represents the water vapor density value of the j-th voxel during the k-th iteration. This represents the relaxation factor, which plays a role in regulating convergence during the iteration process. This represents the element in the i-th row and j-th column of the coefficient matrix of the tomography model. This represents the element in the i-th row of the column vector of the tomography model. Let represent the weight of the i-th observed signal, calculated using formula (11), where Indicates the elevation angle of the i-th GNSS signal; （11） The accuracy index is the root mean square error (RMSE), which is calculated using formula (12): （12） in, This represents the water vapor density value of the i-th voxel calculated by the chromatography model. This represents the water vapor density value of the same voxel provided by the radiosonde station, and r represents the total number of voxels.

7. The GNSS water vapor chromatography method based on dynamically variable-scale voxels according to claim 6, characterized in that, For the horizontal and vertical constraint observation equations, their corresponding weight information All are 1.

8. The GNSS water vapor chromatography method based on dynamically variable-scale voxels according to claim 6, characterized in that, When the root mean square error (RMSE) is less than 1.5 g / m 3 At that time, the accuracy of the water vapor chromatography results met the requirements.

Images