Water vapor chromatography method combined with GNSS / FY-4A

Abstract

The invention relates to a GNSS / FY-4A combined water vapor chromatography method, which comprises the following steps of: processing layered water vapor data of an FY-4A satellite by using GNSS atmospheric precipitable water to obtain corrected FY-4A PWV data; calculating a GNSS / FY-4A SWV signal emitted from the boundary by using the observation data of the GNSS and the corrected FY-4A PWV data; and fusing the FY-4A SWV signal after the boundary signal is recovered into a GNSS water vapor observation equation, constructing a joint GNSS / FY-4A water vapor chromatography model considering the boundary ejection signal, and solving the joint GNSS / FY-4A water vapor chromatography model by using a joint algebraic reconstruction method to obtain a water vapor density result. According to the method, the utilization rate of the signal lines in the chromatography area is greatly improved, the ASIRT is used for resolving the water vapor chromatography model, the real water vapor distribution is better met, and the water vapor chromatography precision is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of three-dimensional water vapor chromatography, and in particular to a water vapor chromatography method combining GNSS / FY-4A. Background Art

[0002] Water vapor is a very active and variable element in the atmosphere, and plays a vital role in atmospheric energy transfer and weather systems. Therefore, water vapor content is considered one of the most important parameters for studying changes in the Earth's atmosphere. In order to obtain the water vapor content and distribution in the atmosphere and provide important data support for weather forecasts and climate research, water vapor detection is often carried out by remote sensing satellite inversion or by using radiosondes, microwave radiometers, sun photometers, etc.

[0003] The Global Navigation Satellite System (GNSS) signal does not arrive at the receiver in a straight line after being sent from the satellite. It will bend due to atmospheric refraction when passing through the troposphere. The resulting delay is called tropospheric delay. When scholars from various countries were studying how to improve positioning accuracy and reduce positioning errors, they unexpectedly discovered that this delay can be used to invert the content of some elements in the atmosphere, and the accuracy of the final result can reach the millimeter level. GNSS has thus begun to intervene in the field of meteorology, bringing new solutions to atmospheric water vapor detection, and gradually developed into an interdisciplinary subject, namely GNSS meteorology. In addition, GNSS / MET (Meteorology) has the characteristics of all-weather, high precision, unaffected by climate, near real-time, wide coverage and low cost, which perfectly meets the many requirements of water vapor monitoring. At present, the GNSS two-dimensional water vapor inversion technology is relatively mature, but the water vapor distribution outside the measuring station can only be obtained by spatial interpolation, and the three-dimensional distribution of water vapor can only be mastered by traditional equipment such as sounding balloons.

[0004] In the field of three-dimensional water vapor tomography, researchers have attempted to achieve and optimize the accuracy of three-dimensional water vapor tomography by adaptively dividing grids, adding empirical constraints, and fusing multivariate data to increase observation signals. However, compared with the number of grids in the entire tomography area, ground GNSS stations are relatively scarce, and it is not realistic to deploy GNSS stations under all grids, which means that the coverage rate of satellite signals in the bottom layer of the study area is low. This leads to problems such as morbidity of the tomography coefficient matrix and poor tomography quality of the bottom grid.

[0005] At present, the accuracy of GNSS two-dimensional water vapor inversion can reach 1 to 2 mm, which can meet the needs of weather forecasting. However, GNSS three-dimensional water vapor tomography is limited by the defects of its signal spatial structure. The underlying grid lacks signal line punctures, resulting in the ill-posed tomography coefficient equation. Its accuracy cannot meet the requirements of water vapor applications. Summary of the invention

[0006] The objective of the present invention is to propose a water vapor tomography method combining GNSS / FY-4A to solve the problems existing in the above-mentioned prior art. The present invention proposes a technology that improves the phenomenon of no signal line piercing in the bottom layer grid caused by the inverted conical structure of the GNSS signal and restores the boundary outgoing signal by introducing the water vapor data of the Fengyun series geostationary meteorological satellite No. 4A (FY-4A) and relying on the remote sensing satellite signal, realizes the complementary advantages, and achieves the purpose of improving the three-dimensional water vapor inversion accuracy.

[0007] To achieve the above objective, the present invention provides the following solutions:

[0008] The water vapor tomography method combining GNSS / FY-4A includes:

[0009] Processing the stratified water vapor data of the FY-4A satellite by using the GNSS precipitable water vapor to obtain the corrected FY-4A PWV data;

[0010] Calculating the GNSS / FY-4A SWV signal outgoing from the boundary of the water vapor tomography area by using the GNSS observation data and the corrected FY-4A PWV data;

[0011] Integrating the FY-4A SWV signal with the restored boundary outgoing signal into the GNSS water vapor observation equation, constructing a combined GNSS / FY-4A water vapor tomography model, and using the combined algebraic reconstruction method to solve the combined GNSS / FY-4A water vapor tomography model to obtain the water vapor density result.

[0012] Optionally, obtaining the corrected FY-4A PWV data includes:

[0013] Performing matching geographical coordinate processing on the stratified water vapor data of the FY-4A satellite to obtain the FY-4A PWV data;

[0014] Constructing a random forest water vapor correction model integrating elevation;

[0015] Using the random forest water vapor correction model to correct the FY-4A PWV data to obtain the corrected FY-4A PWV data.

[0016] Optionally, calculating the FY-4A SWV signal outgoing from the boundary includes:

[0017] Calculating GNSS SWV by using the GNSS observation data;

[0018] Calculating FY-4A SWV by using the GNSS observation data and the corrected FY-4A PWV data; wherein, the observation data includes: elevation angle, azimuth angle and atmospheric gradient information;

[0019] Calculate the GNSS / FY-4A SWV signal emitted from the boundary using the GNSS SWV and FY-4A SWV mentioned above.

[0020] Optionally, the GNSS SWV is:

[0021]

[0022] where θ and represent the satellite altitude angle and azimuth angle respectively, M wet (θ) represents the wet mapping function, M Δ (θ) represents the atmospheric horizontal gradient mapping function, G NS and G WE represent the gradients in the north-south and east-west directions respectively, R θ represents the observed value phase residual, represents the slant path water vapor content of the calculated GNSS signal, and Π represents the water vapor conversion coefficient.

[0023] Optionally, the FY-4A SWV is:

[0024]

[0025] where SWV FY-4A represents the SWV value corresponding to the signal line from each pixel in the LPW product to the satellite, and PWV FY-4A represents the FY-4A PWV corrected by the random forest model integrating elevation.

[0026] Optionally, the GNSS / FY-4A SWV signal emitted from the boundary is:

[0027]

[0028] where SWV H represents the SWV that passes through the signal line on the side at a certain elevation, and PWV H represents the PWV signal passing through the tomography boundary at a certain elevation.

[0029] Optionally, integrating the FY-4A SWV signal that restores the signal emitted from the boundary into the water vapor observation equation of GNSS to construct a combined GNSS / FY-4A water vapor tomography model includes:

[0030] Integrate the FY-4A SWV signal that restores the signal emitted from the boundary into the water vapor observation equation of GNSS to construct a tomography observation equation;

[0031] Add constraint conditions to the tomography observation equation to construct a combined GNSS / FY-4A water vapor tomography model.

[0032] Optionally, the tomography observation equation is:

[0033]

[0034] where SWV i represents the slant-path water vapor content corresponding to the i-th ray, a i,j represents the intercept of the i-th signal line within the j-th grid, x j represents the water vapor density of the j-th grid;

[0035] The constraint conditions include:

[0036] Horizontal constraint:

[0037] w1x1 + … + w i-1 x i-1 -x i +w i+1 x i+1 +… + w n x n = 0

[0038] where x i represents the water vapor density of the selected grid, w i represents the weighting coefficient corresponding to each grid, vertical constraint:

[0039]

[0040] where x i,j,k+1 is the water vapor density of the upper-layer grid, h k -h k+1 is the vertical spacing between adjacent two-layer grids, h sc is the water vapor scale height, e is the base of the natural logarithm, an irrational number, and its value is approximately 2.71828;

[0041] The combined GNSS / FY-4A water vapor tomography model is:

[0042] y = A·x + ε

[0043] where y is the result vector composed of GNSS SWV, FY-4A SWV and constraints, A is the coefficient matrix composed of the intercepts of the signal lines in the grid, horizontal and vertical constraint weight values, x is the water vapor density of the grid, and ε is the residual vector to be allocated.

[0044] Optionally, the combined algebraic reconstruction method is:

[0045]

[0046] where, represents the water vapor density value of the k + 1-th iteration of the j-th element in the grid, λ represents the iteration relaxation factor, represents the value of the previous iteration, m represents the number of signal lines, i represents the i-th ray, n represents the number of grid cells, j represents the j-th pixel, a ij represents the intercept of the i-th ray in the j-th pixel, P i represents the weight information of the i-th GNSS signal, SWV i represents the slant-path water vapor content corresponding to the i-th ray.

[0047] The beneficial effects of the present invention are as follows:

[0048] The present invention integrates FY-4A data to effectively constrain the underlying grid, solves the problem of signal structure defects in traditional water vapor tomography, improves the overall grid piercing rate of tomography, and alleviates the ill-posedness of the tomography equation system. Restoring the GNSS / FY-4A SWV signal emitted from the boundary greatly improves the utilization rate of signal lines in the tomography area. Using ASIRT to solve the water vapor tomography model is more consistent with the actual water vapor distribution and improves the accuracy of water vapor tomography. BRIEF DESCRIPTION OF THE DRAWINGS

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

[0050] Figure 1 is the water vapor correction flowchart of the embodiment of the present invention;

[0051] Figure 2 is the schematic flowchart of the water vapor tomography method combining GNSS / FY-4A of the embodiment of the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0052] The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the drawings in the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, rather than all embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts belong to the scope of protection of the present invention.

[0053] To make the above objects, features, and advantages of the present invention more obvious and understandable, the present invention will be further described in detail below with reference to the drawings and specific embodiments.

[0054] In this embodiment, aiming at problems such as the limited coverage of water vapor inversion stations, low three-dimensional tomography accuracy, and ill-conditioned tomography coefficient matrix, the problem of insufficient inversion station density is solved by interpolation and fusion of multi-source data. In the field of three-dimensional water vapor tomography, this embodiment solves the contradiction of low station coverage rate by means of adaptive grid division, adding empirical constraints, and fusing multi-source data to increase observation signals. And the FY-4A data is fused for water vapor tomography correction, effectively improving the defects of the GNSS signal spatial structure and enhancing the water vapor accuracy. At the same time, due to the differences in the accuracy of different data, when jointly using GNSS / FY-4A data for water vapor tomography, this technology also solves difficult problems such as FY-4A water vapor correction, FY-4A SWV estimation, and solution of tomography observation equations respectively.

[0055] As Figure 2 shown, this embodiment proposes a water vapor tomography method combining GNSS / FY-4A, including:

[0056] Processing the stratified water vapor data of the FY-4A satellite using the GNSS precipitable water vapor to obtain the corrected FY-4A PWV (Precipitable Water Vapor content observed by the Fengyun-4A geostationary meteorological satellite) data;

[0057] Calculating the GNSS / FY-4A SWV (slant path water vapor content of the Global Navigation Satellite System / Fengyun-4A geostationary meteorological satellite) signals emitted from the boundary of the water vapor tomography area using the GNSS observation data and the corrected FY-4A PWV data;

[0058] Integrating the FY-4A SWV signals that restore the signals emitted from the boundary into the GNSS water vapor observation equation, constructing a combined GNSS / FY-4A water vapor tomography model, and using the combined algebraic reconstruction method to solve the combined GNSS / FY-4A water vapor tomography model to obtain the water vapor density result.

[0059] Specifically, this embodiment involves three aspects: 1. Calculating the airborne PWV (Precipitable Water Vapor) data; 2. Separately calculating the GNSS SWV (slant path water vapor content) data, FY-4A SWV data, and the SWV data of the signals penetrating the boundary; 3. Fusing the GNSS / FY-4A data that restores the signals emitted from the boundary to construct a tomography equation and solve the water vapor equation.

[0060] Furthermore, obtaining the corrected FY-4A PWV data includes:

[0061] Performing geographical coordinate matching processing on the stratified water vapor data of the FY-4A satellite to obtain the FY-4A PWV data;

[0062] Construct a random forest water vapor correction model integrating elevation;

[0063] Use the random forest water vapor correction model to correct the FY-4A PWV data and obtain the corrected FY-4A PWV data.

[0064] Specifically, in this embodiment, GNSS-inverted water vapor data is used to extract the FY-4A PWV data at the GNSS stations. The high-precision global navigation satellite system precipitable water vapor is used to process the layered water vapor product (LPW, Layered Precipitable Water Vapor) of the FY-4A satellite, and a random forest water vapor correction model integrating elevation is proposed to obtain the corrected FY-4A PWV product. The specific process is as Figure 1 shown.

[0065] Furthermore, calculating the FY-4A SWV signal emitted from the boundary includes:

[0066] Calculate GNSS SWV using the GNSS observation data;

[0067] Calculate FY-4A SWV using the GNSS observation data and the corrected FY-4A PWV data; where the observation data includes: elevation angle, azimuth angle, and atmospheric gradient information;

[0068] Calculate the GNSS / FY-4A SWV signal emitted from the boundary using GNSS SWV and FY-4A SWV.

[0069] Specifically, in this embodiment, GNSS SWV and FY-4A SWV are respectively calculated using the elevation angle, azimuth angle, and atmospheric gradient information of the observation signal, and then a water vapor proportion factor is constructed to calculate the GNSS / FY-4A SWV signal emitted from the boundary. It includes:

[0070] A. GNSS SWV calculation

[0071] GNSS SWV represents the total water vapor amount per unit volume of the atmosphere along the signal path. With the help of the mapping function, the PWV above the GNSS station can be projected onto the GNSS signal path to obtain GNSS SWV. The relationship formula for calculating SWV from PWV considering the atmospheric horizontal gradient and the observed value phase residual above the station is shown in Equation (1):

[0072]

[0073] In Equation (1), θ and are the satellite elevation angle and azimuth angle, M wet (θ) is the wet mapping function, M Δ(θ) is the atmospheric horizontal gradient mapping function, G NS and G WE Represent the north-south and east-west gradients, R θ is the observed phase residual.

[0074] B. FY-4A SWV calculation

[0075] When the FY-4A satellite scans and images, it continuously receives water vapor signals reflected by ground pixels. Similar to GNSS signals, these signals are processed to obtain the atmospheric precipitable water value PWV in the zenith direction of the pixel. In order to restore the SWV value corresponding to the pixel to the satellite signal line, it is also necessary to obtain the atmospheric horizontal gradient value, projection function value, and altitude and azimuth information of each signal at the pixel position. Referring to formula (1), the calculation method of FY-4A SWV is shown in formula (2):

[0076]

[0077] In formula (2), M wet (θ), Π, M Δ The meanings of parameters such as (θ) and SWV are the same as those in formula (1). FY-4A Represents the SWV value corresponding to the signal line from each pixel to the satellite in the LPW product, PWV FY-4A Represents the FY-4APWV corrected by the random forest model fused with elevation.

[0078] C. Calculation of boundary crossing signal SWV

[0079] When the traditional tomography method establishes the coefficient matrix, it is necessary to eliminate the signal line that passes through the side of the tomography area, because the distance that the signal line passes through the tomography area cannot represent all the water vapor content it contains, which leads to a large number of side-penetrating signals not being effectively utilized and the grids near the tomography boundary lack effective constraints. If these signals are restored, the number of observation equations can be greatly increased. Starting from the SWV calculation formula, the calculation formula for the side-penetrating signal line SWV can be derived:

[0080]

[0081] The meanings of the parameters in formula (3) are the same as those in formula (2).

[0082] The PWV is obtained by integrating the water vapor density along the vertical direction, as shown in formula (4):

[0083]

[0084] Among them, ρ w represents water vapor density, h represents the current plane height, and H represents the cutoff height.

[0085] The antiderivative of Equation (4) is obtained by the method of indefinite integration, and then the Newton-Leibniz formula is used to substitute the upper and lower limits into the antiderivative to find the difference, and Equation (5) can be obtained:

[0086]

[0087] After the parameters in Equation (5) are obtained, the PWV values at different heights can be calculated. Considering that the PWV will change on the day of the tomography experiment, but the ratio between the vertical water vapor densities will not change much in a short time. Therefore, the ratio between PWV H and PWV 总 will also remain unchanged.

[0088] The PWV corresponding to different heights can be calculated according to Equation (5) H , and the water vapor proportion factor is calculated according to Equation (6). Substituting the height H at which the signal line penetrates out of the side wall into Equation (7), the PWV H value at the specified height on the day can be calculated inversely. For the FY-4A pixel points outside the tomography area, whose signal lines penetrate into the side wall and penetrate out of the top layer, the PWV outside the tomography area can also be calculated by referring to the above method H , and then subtracting PWV H from the total PWV of the pixel can obtain the PWV of the signal line within the tomography area.

[0089]

[0090] PWV H =λ H *PWV 总 (7)

[0091] After calculating the PWV H value and substituting it into Equation (3), the corresponding SWV H of the signal line can be restored.

[0092] Furthermore, integrating the FY-4A SWV signal that restores the signal emitted from the boundary into the water vapor observation equation of GNSS, the joint GNSS / FY-4A water vapor tomography model is constructed, including:

[0093] Integrating the FY-4A SWV signal that restores the signal emitted from the boundary into the water vapor observation equation of GNSS to construct a tomography observation equation;

[0094] Adding constraint conditions to the tomography observation equation to construct a joint GNSS / FY-4A water vapor tomography model.

[0095] Specifically, in this embodiment, the FY-4A SWV that restores the boundary emission signal is incorporated into the GNSS water vapor observation equation, constraint conditions are added, a tomographic observation equation for the combined GNSS / FY-4A is constructed, and the water vapor tomography model is solved using the ASIRT algorithm. This includes:

[0096] The propagation speed of the GNSS signal is the same as the speed of light. When a single signal line penetrates, the water vapor within the tomographic area remains constant. However, tomography requires a large number of signal lines, which means that signals at different times need to be incorporated into one equation. Assuming that the water vapor density within the tomographic area remains unchanged during this time window and the water vapor density is equal everywhere within each grid. Then, the sum of the products of the water vapor density within the grid and the intercept of the signal line of that grid is equal to the SWV of that signal line. Calculate the intercepts and SWVs of all signal lines in a loop and combine them into an observation equation:

[0097] SWV = A·X (8)

[0098] In Equation (8), A represents the coefficient matrix, X represents the water vapor density matrix of all grids, and SWV represents the corresponding observed value matrix. Solving this equation can calculate the water vapor density within each grid, and the entire process of construction and solution is called water vapor tomography.

[0099] The fitness function for the tomographic observation equation is expressed as:

[0100] minf(x) = (SWV - A GNSS ·x) T ·P·(SWV - A GNSS ·x), x ∈ R + (9)

[0101] In Equation (9), P is the identity matrix.

[0102] Based on Equation (9), establish the relationship equation between the SWV of each signal line and the intercept and water vapor density parameters within the grid. Combining all the equations can obtain the tomographic observation matrix, as shown in Equation (10):

[0103]

[0104] Among them, a i,j represents the intercept of the i-th signal line within the j-th grid, and x j represents the water vapor density of the j-th grid.

[0105] To alleviate the ill-posedness of the tomography equations, the present invention chooses to introduce FY-4A data so that each underlying grid contains pixel points. If the quality of the LPW product is extremely good during a certain period, it can ensure that all grids have signal lines passing through. However, in reality, the LPW product is always affected by the environment, and there is no PWV value for pixels in some areas, so it cannot be guaranteed that all grids have signal lines passing through. At this time, it is necessary to add certain constraint conditions to make the tomography model reach the solvable condition. Usually, the added constraint conditions are divided into two categories: horizontal and vertical.

[0106] (1) Horizontal constraint. Considering the actual distribution law of water vapor, it can be considered that the change of water vapor density in the horizontal direction is continuous and smooth. The water vapor density of the discretized grid should also be like this. Based on this, the horizontal constraint equation (11) can be established to add constraint conditions for those grids without tomography signals passing through:

[0107] w1x1+…+w i-1 x i-1 -x i +w i+1 x i+1 +…+w n x n =0 (11)

[0108] Among them, x i represents the water vapor density of the selected grid, and w i represents the weighting coefficient corresponding to each grid. This coefficient is related to the horizontal distance between grids and can be calculated using the Gaussian weighting function, as shown in equation (12).

[0109]

[0110] Among them, i, j, and k respectively represent the row, column, and height numbers of the selected grid, in, jn, and k represent the positions of other grids horizontally, d in,kn,k represents the horizontal distance from the center point of other grids to the center point of the selected grid, σ is the smoothing factor, and the larger its value is set, the larger the stable range is and the better the smoothing degree is. Generally, it is set to 1-2 times the horizontal grid spacing. In particular, when i = in and j = jn, the value is 1.

[0111] (2) Vertical constraint. As mentioned above, water vapor decreases exponentially with the increase of elevation in the vertical direction. Although the water vapor density of the discretized grid is not continuous, the ratio change of the water vapor density values between the upper and lower layers also follows this law. Therefore, the vertical constraint can be established using the exponential function, as shown in equation (13).

[0112]

[0113] Equation (13) is in the form of the ratio of the water vapor density of the upper and lower layer grids, where x i,j,k+1is the water vapor density of the upper - layer grid, h k -h k+1 is the vertical spacing between adjacent two - layer grids, h sc is the water vapor elevation.

[0114] Convert (13) into the form of subtracting the water vapor densities of two layers, as shown in Equation (14):

[0115]

[0116] The parameter meanings in Equation (14) are the same as those in Equation (13).

[0117] Expand the formula to each layer of the grid, that is, form the vertical constraint equation, as shown in Equation (15):

[0118] 0 = H v X (15)

[0119] In Equation (15), H v represents the coefficient matrix of the vertical constraint equation, and this coefficient can be calculated according to the local sounding data.

[0120] Combining the above - mentioned observation equation and constraint conditions, the final expression of the water vapor tomography observation equation for the combined GNSS / FY - 4A data can be formed, as shown in Equation (16):

[0121]

[0122] Equation (16) can be abbreviated as:

[0123] y = A·x+ε (17)

[0124] Among them, y is the result vector composed of GNSS SWV, FY - 4A SWV and constraints, A is the coefficient matrix composed of the intercepts of signal lines, horizontal and vertical constraint weight values in the grid, x is the water vapor density of the grid, the dimension number of this vector is the same as the number of grids, and ε is the residual vector to be allocated, which is reflected into x according to the weight during the iterative solution.

[0125] Considering that the SIRT method is not affected by the observation order during the solution and can take into account all observation signals when iteratively allocating residuals, the present invention selects the ASIRT algorithm as the solution method for the combined GNSS / FY - 4A water vapor tomography model. The iterative formula of the ASIRT algorithm is as shown in Equation (18):

[0126]

[0127] This embodiment proposes a three-dimensional water vapor tomography scheme that integrates GNSS and FY-4A observation data. The global pressure and temperature model (Global Pressure and Temperature 3, GPT3) is used to calculate the atmospheric horizontal gradient value at the corresponding position of the pixel, solve the SWV of the FY-4A signal, construct a water vapor proportion factor model, restore the GNSS / FY-4A SWV signal emitted from the boundary, and then integrate the FY-4A SWV of the restored boundary-emitted signal into the GNSS water vapor observation equation to construct a combined GNSS / FY-4A water vapor tomography model. The combined algebraic reconstruction method is used to solve the tomography equation, significantly improving the water vapor tomography accuracy. The two achieve complementary advantages, aiming to improve the accuracy of three-dimensional water vapor inversion.

[0128] The integration of FY-4A data in this embodiment can effectively constrain the underlying grid, solve the problem of signal structure defects in traditional water vapor tomography, improve the overall grid piercing rate of tomography, and alleviate the ill-posedness of the tomography equation system. Restoring the GNSS / FY-4A SWV signal emitted from the boundary greatly improves the utilization rate of the signal lines in the tomography area. Using the ASIRT to solve the water vapor tomography model is more consistent with the actual water vapor distribution, improving the accuracy of water vapor tomography.

[0129] The embodiments described above are only descriptions of the preferred embodiments of the present invention, and do not limit the scope of the present invention. Without departing from the design spirit of the present invention, various deformations and improvements made by those of ordinary skill in the art to the technical solutions of the present invention shall fall within the protection scope determined by the claims of the present invention.

Claims

1. A water vapor tomography method combining GNSS / FY-4A, characterized in that Including: Processing the stratified water vapor data of FY-4A satellite using GNSS precipitable water vapor to obtain the corrected FY-4A PWV data; Calculating the GNSS / FY-4A SWV signal emitted from the boundary of the water vapor tomography region using the GNSS observation data and the corrected FY-4A PWV data; Incorporating the FY-4A SWV signal after restoring the boundary signal into the GNSS water vapor observation equation, constructing a combined GNSS / FY-4A water vapor tomography model, and solving the combined GNSS / FY-4A water vapor tomography model using the combined algebraic reconstruction method to obtain the water vapor density result.

2. The water vapor tomography method combining GNSS / FY-4A according to claim 1, characterized in that, Obtaining the corrected FY-4A PWV data includes: Performing geospatial coordinate matching processing on the stratified water vapor data of FY-4A satellite to obtain FY-4A PWV data; Constructing a random forest water vapor correction model integrating elevation; Using the random forest water vapor correction model to correct the FY-4A PWV data to obtain the corrected FY-4A PWV data.

3. The GNSS / FY-4A combined water vapor tomography method according to claim 1, characterized in that Calculating the FY-4A SWV signal emitted from the boundary includes: Calculating GNSS SWV using the GNSS observation data; Calculating FY-4A SWV using the GNSS observation data and the corrected FY-4A PWV data; wherein, the observation data includes: elevation angle, azimuth angle, and atmospheric gradient information; Calculating the GNSS / FY-4A SWV signal emitted from the boundary using the GNSS SWV and FY-4A SWV.

4. The water vapor tomography method combining GNSS / FY-4A according to claim 3, wherein The GNSS SWV is: where, θ and represent the satellite elevation angle and azimuth angle respectively, M wet (θ) represents the wet mapping function, M Δ (θ) represents the atmospheric horizontal gradient mapping function, G NS and G WE represent the north-south and east-west gradients respectively, R θ represents the observed value phase residual, represents the slant-path water vapor content of the calculated GNSS signal, and Π represents the water vapor conversion coefficient.

5. The water vapor tomography method combining GNSS / FY-4A according to claim 3, wherein The FY-4A SWV is: Among them, SWV FY-4A represents the SWV value corresponding to the signal line from each pixel in the LPW product to the satellite, and PWV FY-4A represents the FY-4A PWV corrected by the random forest model of the fused elevation.

6. The water vapor tomography method combining GNSS / FY-4A according to claim 3, characterized in that The GNSS / FY-4A SWV signal emitted from the boundary is: Among them, SWV H indicates that the signal line SWV passes through the side at a certain elevation, and PWV H indicates the PWV signal passing through the tomographic boundary at a certain elevation.

7. The water vapor tomography method combining GNSS / FY-4A according to claim 1, wherein Incorporating the FY-4A SWV signal after restoring the boundary emitted signal into the GNSS water vapor observation equation to construct a combined GNSS / FY-4A water vapor tomography model includes: Incorporating the FY-4A SWV signal after restoring the boundary emitted signal into the GNSS water vapor observation equation to construct a tomography observation equation; Adding constraint conditions to the tomography observation equation to construct a combined GNSS / FY-4A water vapor tomography model.

8. The GNSS / FY-4A combined water vapor tomography method according to claim 7, characterized in that The tomography observation equation is: Among them, SWV i represents the slant-path water vapor content corresponding to the i-th ray, a i,j represents the intercept of the i-th signal line within the j-th grid, x j represents the water vapor density of the j-th grid; The constraint conditions include: Horizontal constraint: w1x1 + … + w i-1 x i-1 -x i +w i+1 x i+1 +… + w n x n =0 where x i represents the water vapor density of the selected grid, and w i represents the weighting coefficient corresponding to each grid. Vertical constraint: where x i,j,k+1 is the water vapor density of the upper grid, h k -h k+1 is the vertical spacing between adjacent grids, h sc is the water vapor elevation, and e is the base of the natural logarithm; The combined GNSS / FY-4A water vapor tomography model is: y = A·x + ε wherein, y is the result vector composed of GNSS SWV, FY-4A SWV, and constraints, A is the coefficient matrix composed of the signal line intercepts, horizontal, and vertical constraint weight values in the grid, x is the grid water vapor density, and ε is the residual vector to be allocated.

9. The water vapor tomography method combining GNSS / FY-4A according to claim 7, wherein The combined algebraic reconstruction method is: Among them, represents the water vapor density value of the j-th element in the grid at the (k + 1)-th iteration, λ represents the iteration relaxation factor, represents the previous iteration value, m represents the number of signal lines, i represents the i-th ray, n represents the number of grid cells, j represents the j-th pixel, a ij represents the intercept of the i-th ray in the j-th pixel, P i represents the weight information of the i-th GNSS signal, SWV i represents the slant path water vapor content corresponding to the i-th ray.

Images