A fast water vapor tomography method fusing real-time beidou observation and background field constraint

Abstract

The application discloses a kind of fusion real-time Beidou observation and background field constraint's quick water vapor tomography method, comprising: one, based on GNSS inversion real-time atmospheric precipitable water PWV;Two, obtain grid data from NCEP-GFS forecast data and interpolation obtain water vapor density background field;Three, based on atmospheric precipitable water PWV and water vapor density background field obtain SWV residual;Four, based on SWV residual, obtain the water vapor density average correction coefficient corresponding to GNSS signal ray;Five, based on water vapor density average correction coefficient, carry out water vapor density weight correction, obtain the water vapor density correction value of each grid;Six, water vapor density correction value is superimposed to original water vapor density background field, obtains high-precision water vapor density field and realizes quick water vapor tomography.The application no longer introduces horizontal constraint, vertical constraint and the like equation, based on SWV residual, carries out correction to water vapor density, improves precision and solving efficiency.

Description

Technical Field

[0001] This invention belongs to the field of atmospheric water vapor information inversion technology, specifically involving a rapid water vapor tomography method that integrates real-time BeiDou observations with background field constraints. Background Technology

[0002] With the rapid development of Global Navigation Satellite System (GNSS), GNSS water vapor tomography has become an important technique for detecting three-dimensional water vapor due to its advantages such as high precision and high resolution. Water vapor tomography utilizes delayed observations of GNSS signals as they pass through the atmosphere to invert the three-dimensional water vapor density distribution within the troposphere.

[0003] Current tomographic methods typically incorporate tomographic equations with horizontal, vertical, and prior constraints to improve the stability of water vapor density inversion. Regarding horizontal constraints, a Gaussian weighting function is usually used to establish the water vapor density relationship between grids in the horizontal direction. The core parameter of this function is the smoothing factor, whose value is largely determined empirically. In regions with complex lower-level water vapor structures or under severe convective weather conditions, this can lead to significant deviations in the water vapor density inversion results. Regarding vertical constraints, existing tomographic methods rely heavily on empirical settings for water vapor elevation parameters, resulting in vertical biases in the tomographic results.

[0004] Regarding prior constraints, existing tomographic methods typically incorporate radiosonde data, ERA5 reanalysis data, or mid-latitude atmospheric empirical models as prior information to compensate for insufficient observational data. However, this approach not only introduces systematic biases but also, when jointly constraining prior information with GNSS observations, the constraint weights rely on empirical assignments, leading to a significant impact of prior information on the inversion results.

[0005] In addition, due to the uneven distribution of satellite signals, there are grids that are not penetrated by satellite rays. These grids have no effective constraints in the observations, causing ill-posed problems in the tomographic equations.

[0006] Therefore, existing tomography methods generally suffer from problems such as strong empirical dependence of parameters and complex model structures, low computational efficiency, and ill-posedness due to the introduction of multiple constraints. Summary of the Invention

[0007] The technical problem to be solved by this invention is to address the shortcomings of the prior art by providing a rapid water vapor tomography method that integrates real-time BeiDou observations and background field constraints. This method obtains the average correction coefficient of water vapor density based on SWV residuals, corrects the grids traversed by GNSS signal rays using a single weight coefficient or a combination of weight coefficients and observation confidence, and then corrects the grids not traversed by GNSS signal rays. The corrected water vapor density values ​​are then superimposed on the original water vapor density background field to obtain a high-precision water vapor density field, thereby achieving rapid water vapor tomography and improving the accuracy of subsequent weather forecasts.

[0008] To solve the above-mentioned technical problems, the technical method adopted in this invention is: a rapid water vapor tomography method integrating real-time BeiDou observation and background field constraints, which includes the following steps: Step 1: Real-time atmospheric precipitable water volume (PWV) retrieved based on GNSS inversion; Step 2: Obtain grid data from NCEP-GFS forecast data and interpolate to obtain the water vapor density background field; Step 3: Obtain the SWV residual based on the atmospheric precipitable water volume (PWV) and water vapor density background field; Step 4: Obtain the average correction coefficient of water vapor density corresponding to the GNSS signal rays based on the SWV residual; Step 5: Perform water vapor density weight correction based on the average water vapor density correction coefficient to obtain the water vapor density correction value for each grid: Step 501: Determine the number of GNSS signal rays passing through each grid in the water vapor density background field. If the number of GNSS signal rays passing through the grid is 1, proceed to step 502; if the number of GNSS signal rays passing through the grid is greater than 1, proceed to step 503; if the grid is not passed through by GNSS signal rays, proceed to step 504. Step 502: Correct the average water vapor density correction coefficient with weighting coefficient to obtain the water vapor density correction value; Step 503: Correct the average correction coefficient of water vapor density by applying weighting coefficients and observation confidence to obtain the corrected value of water vapor density; Step 504: Correct the water vapor density of the grid that was not pierced by the GNSS signal rays based on the water vapor density correction value of the grid that was pierced by the GNSS signal rays, and obtain the water vapor density correction value of the grid that was not pierced by the GNSS signal rays. Step 6: Superimpose the water vapor density correction value onto the original water vapor density background field to obtain a high-precision water vapor density field and achieve rapid water vapor chromatography.

[0009] The above-mentioned rapid water vapor tomography method integrating real-time BeiDou observations and background field constraints is further refined in the following steps: Step 201: Obtain the three-dimensional grid data of the tomographic region from the NCEP-GFS forecast data; Step 202: Obtain the thermodynamic temperature and vapor pressure of each grid from the three-dimensional grid data of the tomographic region; Step 203, according to The water vapor density of each grid is obtained. ;in, The water vapor pressure of the grid. The thermodynamic temperature of the grid. The constant for water vapor; Step 204: Interpolate the water vapor density of multiple grids using an exponential function in the elevation direction and bilinear interpolation in the latitude and longitude directions to obtain the spatially interpolated water vapor density background field.

[0010] The aforementioned rapid water vapor tomography method integrating real-time BeiDou observations and background field constraints is further refined in the following steps: Step 2041: Record the grid in the latitude and longitude direction as the horizontal grid, and the grid in the elevation direction as the vertical grid; Step 2042: Denote the multiple vertical grids with the same latitude and longitude as the i-th horizontal grid as the i-th group of vertical grids; where i is a positive integer; Step 2043: Obtain the water vapor density of the i-th vertical grid from the water vapor density of multiple grids; Step 2044: Construct the i-th exponential function ρ i(Z) =ρ i0 ×exp(-Z / H); where ρ i(Z) To determine the water vapor density at a location with the same latitude and longitude as the i-th horizontal grid and an elevation of Z, ρ i0 For the water vapor density at the same latitude and longitude as the i-th horizontal grid and at zero elevation, exp(·) is a power operation with the natural constant e as the base, Z is the elevation, and H is the water vapor elevation parameter; Step 2045: Based on the water vapor density of the i-th set of vertical grids, in ρ i(Z) Given Z, the least squares method is used to fit the i-th exponential function to obtain the water vapor density ρ at the same latitude and longitude as the i-th horizontal grid and at zero elevation. i0 Given the water vapor elevation parameter H, the i-th exponential function after fitting is obtained; Step 2046: Based on the fitted i-th exponential function, at the latitude and longitude of the i-th horizontal grid, interpolate along the elevation direction by inputting different elevations to obtain the water vapor density of different height layers. Step 2047: For the water vapor density at the same altitude layer, use bilinear interpolation in the latitude and longitude directions to obtain the spatially interpolated water vapor density background field.

[0011] The above-mentioned rapid water vapor tomography method integrating real-time BeiDou observations and background field constraints is further refined in the following steps: Step 301: Denote the grids through which any GNSS signal ray passes as the 1st grid, ..., the nth grid, ..., the Nth grid; where n and N are both positive integers, 1 ≤ n ≤ N; Step 302: Obtain the water vapor density of the nth grid from the water vapor density background field. and according to The oblique path water vapor content of the GNSS signal ray based on the water vapor density background field was obtained. ;in, The length of the GNSS signal ray passing through the nth grid; Step 303: Obtain the oblique path water vapor content (SWV) based on the atmospheric precipitable water vapor (PWV) inversion. Step 304, according to =SWV- SWV residuals are obtained. ; Step 305: Following the methods in steps 301 to 304, obtain the SWV residuals corresponding to multiple GNSS signal rays.

[0012] The above-mentioned rapid water vapor tomography method integrating real-time BeiDou observations and background field constraints is further refined in the following steps: Step 401, according to The weighting coefficients of the GNSS signal ray passing through the nth grid are obtained. ; Step 402: Construct the SWV residual equation The average correction coefficient for water vapor density corresponding to the GNSS signal rays was obtained. ; Step 403: Following the methods in steps 401 to 402, obtain the average correction coefficient of water vapor density corresponding to multiple GNSS signal rays.

[0013] The aforementioned rapid water vapor tomography method integrating real-time BeiDou observations and background field constraints is further refined in the following steps: Step 5021: Following the method in step 401, obtain the weight coefficients of the grid. The average water vapor density correction coefficients for the GNSS signal rays passing through the grid are obtained from the average water vapor density correction coefficients corresponding to multiple GNSS signal rays and denoted as [missing value]. ; Step 5022, according to The water vapor density correction value of the grid is obtained. .

[0014] The aforementioned rapid water vapor tomography method integrating real-time BeiDou observations and background field constraints is further refined in the following steps: Step 5031: If the number of GNSS signal rays passing through the grid is greater than 1, then the GNSS signal rays passing through the grid are respectively denoted as the 1st GNSS signal ray, ..., the mth GNSS signal ray, ..., the Mth GNSS signal ray; where M is a positive integer greater than or equal to 2, and 1≤m≤M; Step 5032, according to The elevation angle uncertainty component corresponding to the m-th GNSS signal ray is obtained. ;in, Let m be the elevation angle of the m-th GNSS signal ray; Step 5033, according to The carrier-to-noise ratio uncertainty component corresponding to the m-th GNSS signal ray is obtained. ;in, The signal-to-carrier-to-noise ratio of the m-th GNSS signal ray; Step 5034, according to The residual uncertainty component corresponding to the m-th GNSS signal ray is obtained. ;in, Let m be the SWV residual corresponding to the m-th GNSS signal ray. The median of the SWV residuals. The median of the absolute values ​​of the SWV residuals. It is the absolute value; Step 5035, according to The overall uncertainty corresponding to the m-th GNSS signal ray is obtained. ; Step 5036, according to The observation confidence level corresponding to the m-th GNSS signal ray is obtained. ;in, For adjustment coefficients; Step 5037, according to The water vapor density correction value obtained when the m-th GNSS signal ray passes through the grid is obtained. ;in, This is the average correction factor for the water vapor density of the grid as it is traversed by the m-th GNSS signal ray; The weighting coefficient corresponding to the m-th GNSS signal ray passing through the grid. Step 5038: If the observation confidence level corresponding to all GNSS signal rays is less than the preset threshold, then abandon the observation confidence level correction, and proceed according to... The water vapor density correction value of the grid is obtained. Otherwise, according to The water vapor density correction value of the grid is obtained. .

[0015] The aforementioned rapid water vapor tomography method integrating real-time BeiDou observations and background field constraints is further refined in the following steps: Step 5041: Obtain the water vapor density correction value corresponding to the grid that is traversed by the GNSS signal ray in the height layer where the grid is located as a sample point. If the number of sample points is greater than or equal to 6, proceed to steps 5042 to 5044; if the number of sample points is greater than or equal to 3 but less than 6, proceed to step 5045; if the number of sample points is less than 3, the water vapor density correction value is zero. Step 5042: Construct a bivariate quadratic polynomial regression equation ρ=Ax 2 +By 2 +Cxy+Dx+Ey+F; where A is the first fitting coefficient, B is the second fitting coefficient, C is the third fitting coefficient, D is the fourth fitting coefficient, E is the fifth fitting coefficient, and F is the sixth fitting coefficient; x is the longitude independent variable, y is the latitude independent variable; and ρ is the water vapor density correction value dependent variable. Step 5043: Using the least squares method, based on the sample points in step 5041, fit and solve the bivariate quadratic polynomial regression equation to obtain the first to sixth fitting coefficients. Step 5044: Substitute the longitude and latitude of the grid that were not traversed by the GNSS signal rays in the same altitude layer into the fitted bivariate quadratic polynomial regression equation to obtain the water vapor density correction value of the grid. Step 5045: Construct a bivariate linear polynomial regression equation ρ=A´x+B´y+C´, and obtain the water vapor density correction value of the grid according to the methods in steps 5043 and 5044; where A´ is the first fitting parameter, B´ is the second fitting parameter, and C´ is the third fitting parameter.

[0016] Compared with the prior art, the present invention has the following advantages: 1. This invention solves the problems of existing tomographic methods, such as strong parameter dependence on constraint equations, high complexity, low computational efficiency, and ill-posedness. This invention improves the accuracy and solution efficiency of tomographic results, enabling rapid water vapor tomography, and can correct the water vapor density of grids that have not been penetrated by satellite rays.

[0017] 2. This invention is based on GNSS inversion of real-time atmospheric precipitable water volume (PWV), and inverts oblique path water vapor content (SWV) based on real-time PWV data, providing SWV observation data; The purpose of obtaining gridded data from NCEP-GFS forecast data and interpolating it to obtain the water vapor density background field is to indirectly obtain the slant path water vapor content (SWV) based on the water vapor density background field. Then, the difference between the observed SWV data and the indirectly obtained slant path water vapor content is calculated to obtain the SWV residual.

[0018] 3. This invention connects GNSS observation information with the water vapor density background field through SWV residuals as a key variable for correction. When the water vapor density background field can accurately describe the atmospheric water vapor density distribution, the SWV calculated based on the water vapor density background field should be basically consistent with the SWV of GNSS observation. In this way, the SWV residuals can not only quantify the error of the water vapor density background field, but also serve as an observation constraint, effectively introducing GNSS observation information into the water vapor density background field correction process so as to obtain the average correction coefficient of water vapor density in the future.

[0019] 4. This invention obtains the SWV residual, and since the SWV residual originates from the deviation of each grid on the GNSS signal propagation path based on the background field water vapor density value, the average water vapor density correction coefficient can be obtained through the constructed SWV residual equation.

[0020] 5. When the present invention performs water vapor density correction based on the average water vapor density correction coefficient, it can perform correction based on different situations. When the number of GNSS signal rays passing through the grid is 1, only the average water vapor density correction coefficient is corrected by the weighting coefficient; when the number of GNSS signal rays passing through the grid is greater than 1, it is corrected by the weighting coefficient and the observation confidence level to obtain the water vapor density correction value.

[0021] 6. When the number of GNSS signal rays passing through the grid is greater than 1, the present invention, after correction by weighting coefficient and observation confidence, not only considers the weighting coefficient correction of the water vapor density of each grid relative to the SWV observation data, preserving the difference in the contribution of different grids on the path to the water vapor density deviation, but also considers that the same grid may be passed through by multiple different rays, thus resulting in multiple water vapor density correction values. Through observation confidence, high-quality rays are fully utilized and low-quality rays are suppressed, thereby obtaining more stable results.

[0022] 7. When correcting the water vapor density of grids not traversed by GNSS signal rays, this invention cannot directly correct the data using the methods described above because the grids are not traversed by any GNSS signal rays. However, since the water vapor density at the same altitude layer has a continuous distribution characteristic, the invention fully utilizes the water vapor density correction values ​​corresponding to grids traversed by GNSS signal rays at the same altitude layer. Based on the number of grids traversed by rays at the same horizontal plane, a binary quadratic polynomial or a binary linear polynomial regression equation is fitted and solved. Finally, the longitude and latitude of the grids not traversed by GNSS signal rays at the same altitude layer are substituted into the above polynomial regression equation to obtain the water vapor density correction values ​​of the grids not traversed by GNSS signal rays. This spatial fitting solution fully utilizes the continuous variation characteristic of water vapor density in the horizontal direction, effectively compensating for the information loss problem caused by insufficient satellite observation coverage.

[0023] 8. This invention eliminates the need for equations for horizontal constraints, vertical constraints, and prior constraints, thus reducing the impact of various constraints on the tomographic results. Furthermore, this invention significantly improves both the accuracy and efficiency of water vapor density inversion, providing a simple, effective, and rapid method for GNSS-based water vapor tomography.

[0024] In summary, this invention obtains the average correction coefficient for water vapor density based on SWV residuals, corrects the grids traversed by GNSS signal rays using a single weight coefficient or a combination of weight coefficients and observation confidence, and then corrects the grids not traversed by GNSS signal rays. The corrected water vapor density values ​​are then superimposed onto the original water vapor density background field to obtain a high-precision water vapor density field, thereby achieving rapid water vapor tomography and improving the accuracy of subsequent weather forecasts.

[0025] The technical method of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0026] Figure 1 This is a flowchart of the method of the present invention.

[0027] Figure 2 This is a comparison chart of water vapor density obtained by the method of the present invention and existing methods.

[0028] Figure 3 This is a comparison chart of the RMSE of water vapor density obtained during rainfall periods using the method of this invention and existing methods, and the reference value.

[0029] Figure 4 This is a comparison chart of the RMSE of water vapor density obtained by the method of the present invention and existing methods during non-rainfall periods with reference values.

[0030] Figure 5 This is a comparison chart of the computation time of the method of the present invention and existing methods. Detailed Implementation

[0031] like Figure 1 The rapid water vapor tomography method shown, which integrates real-time BeiDou observations with background field constraints, includes the following steps: Step 1: Real-time atmospheric precipitable water volume (PWV) retrieved based on GNSS inversion; Step 2: Obtain grid data from NCEP-GFS forecast data and interpolate to obtain the water vapor density background field; Step 3: Obtain the SWV residual based on the atmospheric precipitable water volume (PWV) and water vapor density background field; Step 4: Obtain the average correction coefficient of water vapor density corresponding to the GNSS signal rays based on the SWV residual; Step 5: Perform water vapor density weight correction based on the average water vapor density correction coefficient to obtain the water vapor density correction value for each grid: Step 501: Determine the number of GNSS signal rays passing through each grid in the water vapor density background field. If the number of GNSS signal rays passing through the grid is 1, proceed to step 502; if the number of GNSS signal rays passing through the grid is greater than 1, proceed to step 503; if the grid is not passed through by GNSS signal rays, proceed to step 504. Step 502: Correct the average water vapor density correction coefficient with weighting coefficient to obtain the water vapor density correction value; Step 503: Correct the average correction coefficient of water vapor density by applying weighting coefficients and observation confidence to obtain the corrected value of water vapor density; Step 504: Correct the water vapor density of the grid that was not pierced by the GNSS signal rays based on the water vapor density correction value of the grid that was pierced by the GNSS signal rays, and obtain the water vapor density correction value of the grid that was not pierced by the GNSS signal rays. Step 6: Superimpose the water vapor density correction value onto the original water vapor density background field to obtain a high-precision water vapor density field and achieve rapid water vapor chromatography.

[0032] In this embodiment, the specific process of step two is as follows: Step 201: Obtain the three-dimensional grid data of the tomographic region from the NCEP-GFS forecast data; Step 202: Obtain the thermodynamic temperature and vapor pressure of each grid from the three-dimensional grid data of the tomographic region; Step 203, according to The water vapor density of each grid is obtained. ;in, The water vapor pressure of the grid. The thermodynamic temperature of the grid. The constant for water vapor; Step 204: Interpolate the water vapor density of multiple grids using an exponential function in the elevation direction and bilinear interpolation in the latitude and longitude directions to obtain the spatially interpolated water vapor density background field.

[0033] In this embodiment, the specific process of step 204 is as follows: Step 2041: Record the grid in the latitude and longitude direction as the horizontal grid, and the grid in the elevation direction as the vertical grid; Step 2042: Denote the multiple vertical grids with the same latitude and longitude as the i-th horizontal grid as the i-th group of vertical grids; where i is a positive integer; Step 2043: Obtain the water vapor density of the i-th vertical grid from the water vapor density of multiple grids; Step 2044: Construct the i-th exponential function ρ i(Z) =ρ i0 ×exp(-Z / H); where ρ i(Z) To determine the water vapor density at a location with the same latitude and longitude as the i-th horizontal grid and an elevation of Z, ρ i0 For the water vapor density at the same latitude and longitude as the i-th horizontal grid and at zero elevation, exp(·) is a power operation with the natural constant e as the base, Z is the elevation, and H is the water vapor elevation parameter; Step 2045: Based on the water vapor density of the i-th set of vertical grids, in ρ i(Z) Given Z, the least squares method is used to fit the i-th exponential function to obtain the water vapor density ρ at the same latitude and longitude as the i-th horizontal grid and at zero elevation. i0 Given the water vapor elevation parameter H, the i-th exponential function after fitting is obtained; Step 2046: Based on the fitted i-th exponential function, at the latitude and longitude of the i-th horizontal grid, interpolate along the elevation direction by inputting different elevations to obtain the water vapor density of different height layers. Step 2047: For the water vapor density at the same altitude layer, use bilinear interpolation in the latitude and longitude directions to obtain the spatially interpolated water vapor density background field.

[0034] In this embodiment, the specific process of step three is as follows: Step 301: Denote the grids through which any GNSS signal ray passes as the 1st grid, ..., the nth grid, ..., the Nth grid; where n and N are both positive integers, 1 ≤ n ≤ N; Step 302: Obtain the water vapor density of the nth grid from the water vapor density background field. and according to The oblique path water vapor content of the GNSS signal ray based on the water vapor density background field was obtained. ;in, The length of the GNSS signal ray passing through the nth grid; Step 303: Obtain the oblique path water vapor content (SWV) based on the atmospheric precipitable water vapor (PWV) inversion. Step 304, according to =SWV- SWV residuals are obtained. ; Step 305: Following the methods in steps 301 to 304, obtain the SWV residuals corresponding to multiple GNSS signal rays.

[0035] In this embodiment, the specific process of step four is as follows: Step 401, according to The weighting coefficients of the GNSS signal ray passing through the nth grid are obtained. ; Step 402: Construct the SWV residual equation The average correction coefficient for water vapor density corresponding to the GNSS signal rays was obtained. ; Step 403: Following the methods in steps 401 to 402, obtain the average correction coefficient of water vapor density corresponding to multiple GNSS signal rays.

[0036] In this embodiment, the specific process of step 502 is as follows: Step 5021: Following the method in step 401, obtain the weight coefficients of the grid. The average water vapor density correction coefficients for the GNSS signal rays passing through the grid are obtained from the average water vapor density correction coefficients corresponding to multiple GNSS signal rays and denoted as [missing value]. ; Step 5022, according to The water vapor density correction value of the grid is obtained. .

[0037] In this embodiment, the specific process of step 503 is as follows: Step 5031: If the number of GNSS signal rays passing through the grid is greater than 1, then the GNSS signal rays passing through the grid are respectively denoted as the 1st GNSS signal ray, ..., the mth GNSS signal ray, ..., the Mth GNSS signal ray; where M is a positive integer greater than or equal to 2, and 1≤m≤M; Step 5032, according to The elevation angle uncertainty component corresponding to the m-th GNSS signal ray is obtained. ;in, Let m be the elevation angle of the m-th GNSS signal ray; Step 5033, according to The carrier-to-noise ratio uncertainty component corresponding to the m-th GNSS signal ray is obtained. ;in, The signal-to-carrier-to-noise ratio of the m-th GNSS signal ray; Step 5034, according to The residual uncertainty component corresponding to the m-th GNSS signal ray is obtained. ;in, Let m be the SWV residual corresponding to the m-th GNSS signal ray. The median of the SWV residuals. The median of the absolute values ​​of the SWV residuals. It is the absolute value; Step 5035, according to The overall uncertainty corresponding to the m-th GNSS signal ray is obtained. ; Step 5036, according to The observation confidence level corresponding to the m-th GNSS signal ray is obtained. ;in, For adjustment coefficients; Step 5037, according to The water vapor density correction value obtained when the m-th GNSS signal ray passes through the grid is obtained. ;in, This is the average correction factor for the water vapor density of the grid as it is traversed by the m-th GNSS signal ray; The weighting coefficient corresponding to the m-th GNSS signal ray passing through the grid. Step 5038: If the observation confidence level corresponding to all GNSS signal rays is less than the preset threshold, then abandon the observation confidence level correction, and proceed according to... The water vapor density correction value of the grid is obtained. Otherwise, according to The water vapor density correction value of the grid is obtained. .

[0038] In this embodiment, the specific process of step 504 is as follows: Step 5041: Obtain the water vapor density correction value corresponding to the grid that is traversed by the GNSS signal ray in the height layer where the grid is located as a sample point. If the number of sample points is greater than or equal to 6, proceed to steps 5042 to 5044; if the number of sample points is greater than or equal to 3 but less than 6, proceed to step 5045; if the number of sample points is less than 3, the water vapor density correction value is zero. Step 5042: Construct a bivariate quadratic polynomial regression equation ρ=Ax 2 +By 2+Cxy+Dx+Ey+F; where A is the first fitting coefficient, B is the second fitting coefficient, C is the third fitting coefficient, D is the fourth fitting coefficient, E is the fifth fitting coefficient, and F is the sixth fitting coefficient; x is the longitude independent variable, y is the latitude independent variable; and ρ is the water vapor density correction value dependent variable. Step 5043: Using the least squares method, based on the sample points in step 5041, fit and solve the bivariate quadratic polynomial regression equation to obtain the first to sixth fitting coefficients. Step 5044: Substitute the longitude and latitude of the grid that were not traversed by the GNSS signal rays in the same altitude layer into the fitted bivariate quadratic polynomial regression equation to obtain the water vapor density correction value of the grid. Step 5045: Construct a bivariate linear polynomial regression equation ρ=A´x+B´y+C´, and obtain the water vapor density correction value of the grid according to the methods in steps 5043 and 5044; where A´ is the first fitting parameter, B´ is the second fitting parameter, and C´ is the third fitting parameter.

[0039] In this embodiment, it should be noted that the atmospheric precipitable water volume (PWV) retrieved based on GNSS in step one can refer to conventional methods in the field, specifically the GNSS water vapor retrieval method in patent CN202411799199.9, "An Assimilation Method Based on GNSS Multidimensional Water Vapor Cooperative WRF Meteorological Model".

[0040] In this embodiment, it should be noted that the water vapor gas constant... The value is 461.495 J / (kg·K).

[0041] In this embodiment, it should be noted that the NCEP-GFS forecast data is forecast data from the U.S. National Center for Environmental Prediction's Global Forecast System.

[0042] In this embodiment, it should be noted that the oblique path water vapor content (SWV) obtained by inverting atmospheric precipitable water vapor (PWV) can refer to conventional methods in the field, and the wet delay mapping function (Mw) in the inversion process can adopt the conventional global mapping function (GMF) in the field.

[0043] In this embodiment, it should be noted that the number of sample points in step 5041 is generally greater than 3. According to statistics, the average ray coverage of the low-level grid is 50%, the average ray coverage of the middle-level grid is 75.3%, and the average ray coverage of the high-level grid is 86.8%. The number of grids that are traversed by rays is much greater than 3. In actual use, if the number of sample points is less than 3, water vapor density correction is not performed.

[0044] In this embodiment, it should be noted that the adjustment coefficient... The value is 0.01.

[0045] In this embodiment, the step two involves obtaining grid data and interpolating it from the NCEP-GFS forecast data because the NCEP-GFS forecast data divides the continuous atmosphere into a three-dimensional grid, allowing the acquisition of thermodynamic temperature and water vapor pressure for each grid. However, the NCEP-GFS forecast data has a latitude and longitude spatial resolution of 0.25° × 0.25°, and originally consists of 64 mixed σ-vertical model layers in the height direction. After UPP processing, the publicly released GRIB2 isobaric surface product is obtained. The GRIB2 isobaric surface product provides meteorological elements for water vapor chromatography, such as thermodynamic temperature, relative humidity, and air pressure, at 1000 hPa, 975 hPa, 950 hPa, 925 hPa, 900 hPa, 850 hPa, 800 hPa, 750 hPa, and 700 hPa. A total of 26 standard isobaric layers were established, including 650hPa, 600hPa, 550hPa, 500hPa, 450hPa, 400hPa, 350hPa, 300hPa, 250hPa, 200hPa, 150hPa, 100hPa, 70hPa, 50hPa, 30hPa, 20hPa, and 10hPa. The pressure layers were converted into upper terrestrial layers using conventional methods found in the research and application of key technologies for ground-based GNSS water vapor inversion. Therefore, the NCEP-GFS forecast data used 26 upper terrestrial layers for subsequent construction of the tomographic background field. The sampling times for the NCEP-GFS forecast data were daily at UTC 00:00, UTC 06:00, UTC 12:00, and UTC 18:00.

[0046] Since the spatial resolution of the grid division of the NCEP-GFS forecast data is inconsistent with that of the tomographic region grid division, it is necessary to perform spatial interpolation on the water vapor density background field. Exponential interpolation is used in the vertical direction and bilinear interpolation is used in the horizontal direction to obtain the water vapor density background field after spatial interpolation of the tomographic region.

[0047] In this embodiment, the exponential interpolation in the vertical direction in step 204 is based on the fact that the distribution of water vapor density in the vertical direction follows a near-exponential decay law. This allows for the accurate capture of detailed features of near-surface water vapor accumulation and rapid decrease in the upper troposphere. It also has a certain smoothing and filtering effect on local observation noise, avoiding excessive interference of individual outliers on the overall profile, thereby ensuring the stability of the numerical solution of subsequent tomographic inversion. The computational complexity is low, making it easy to analyze the distribution of water vapor density in the vertical direction.

[0048] Bilinear interpolation is used in the horizontal direction of latitude and longitude because the distribution of atmospheric water vapor in the horizontal direction is usually continuous and gradual. As a classic two-dimensional interpolation method, bilinear interpolation is very suitable for such smoothly changing fields and can ensure that the interpolated horizontal surface is continuous, avoiding unnatural jumps in the data and conforming to the characteristics of the horizontal distribution of water vapor. It also has low computational cost and is very suitable for processing large-scale grid data.

[0049] In this embodiment, it should be noted that the oblique path water vapor content (SWV) is the core observation of GNSS water vapor tomography, and its physical essence is the line integral of water vapor density along the satellite signal propagation path. Because the continuous atmosphere is divided into a three-dimensional grid, the tomographic region is usually discretized into several independent three-dimensional grids, thus converting the line integral of water vapor density along the satellite signal propagation path into a discrete form of summation of water vapor density along the satellite signal propagation path.

[0050] In this embodiment, the SWV residual is obtained in step three based on the water vapor density background field and atmospheric precipitable water volume (PWV). This is because the SWV residual is a key variable connecting GNSS observation information with the water vapor density background field correction. When the water vapor density background field can accurately describe the atmospheric water vapor density distribution, the SWV calculated based on the water vapor density background field should be basically consistent with the SWV of GNSS observation. In this way, the SWV residual can not only quantify the water vapor density background field error, but also serve as an observation constraint, effectively introducing GNSS observation information into the water vapor density background field correction process so as to obtain the average correction coefficient of water vapor density in the future.

[0051] In this embodiment, the SWV residual equation constructed in step 402 takes into account that the SWV residual comes from the deviation of each grid on the signal propagation path based on the water vapor density background field, and the water vapor density correction coefficient of the same ray passing through each grid is the same, so it is the average water vapor density correction coefficient, which represents the average deviation of the water vapor density on the GNSS signal ray path relative to the background field.

[0052] In this embodiment, the number of GNSS signal rays passing through the grid is judged in step five so that the number of GNSS signal rays passing through the grid is 1, and only the weight coefficient correction in step 502 is performed; when the number of GNSS signal rays passing through the grid is greater than 1, not only the weight coefficient correction but also the observation confidence correction is performed in step 503, which improves the correction accuracy.

[0053] In this embodiment, when correcting the water vapor density of the grid that was not traversed by GNSS signal rays in step 504, it cannot be directly corrected using the data obtained by the above method because it was not traversed by any GNSS signal rays. However, since the water vapor density at the same altitude layer has a continuous distribution characteristic, the water vapor density correction value corresponding to the grid that was traversed by GNSS signal rays in the same altitude layer is fully utilized. Based on the number of grids traversed by rays on the same horizontal plane, a binary quadratic polynomial or a binary linear polynomial regression equation is fitted and solved. Finally, the longitude and latitude of the grids that were not traversed by GNSS signal rays in the same altitude layer are substituted into the above polynomial regression equation to obtain the water vapor density correction value of the grid that was not traversed by GNSS signal rays. This spatial fitting solution fully utilizes the continuous variation characteristic of water vapor density in the horizontal direction, effectively compensating for the information loss problem caused by insufficient satellite observation coverage.

[0054] In this embodiment, the weighting coefficient correction in steps 5021 and 5037 first considers reasonably distributing the average water vapor density correction coefficient of each GNSS signal ray to each grid through which the ray passes, thereby realizing the reasonable distribution of the overall scale of the ray to each grid scale, while preserving the differences in the contribution of different grids on the path to the water vapor density deviation.

[0055] In this embodiment, steps 5032 to 5034 consider the elevation angle uncertainty component, carrier-to-noise ratio uncertainty component, and residual uncertainty component, respectively. This is because in actual observations, the same grid may be traversed by multiple different rays, resulting in multiple water vapor density correction values. Due to differences in the observation quality of different satellite rays, rays with low elevation angles are more affected by multipath effects, while rays with high elevation angles have relatively stable observation quality. Furthermore, factors such as the carrier-to-noise ratio and elevation angle also affect the SWV inversion accuracy. Therefore, the elevation angle uncertainty component, carrier-to-noise ratio uncertainty component, and residual uncertainty component are multiplied to obtain the comprehensive uncertainty, which in turn yields the observation confidence. This is used to calculate the observation confidence for each ray. Weighted fusion of observation confidence is performed during grid correction, fully utilizing high-quality rays and suppressing low-quality rays, thus obtaining a more stable result.

[0056] In this embodiment, a preset threshold is set with a value of 0.4. If a grid is traversed by multiple GNSS signal rays, but the confidence level of all rays is lower than the preset threshold, then the observation correction is abandoned and the background field value is retained to avoid low-quality observations polluting the background field.

[0057] In this embodiment, in order to verify the effectiveness of the method of the present invention, 13 GNSS stations and 1 radiosonde station in the satellite positioning reference station network (SatRef) of Hong Kong, China were selected for the experiment, with the time periods being July 21-26, 2015 (rainy season) and August 2-6, 2015 (non-rainy season).

[0058] In this embodiment, the tomographic region (i.e., the Hong Kong region of China) ranges from 113.87°E to 114.35°E, from 22.2°N to 22.55°N, and from 0 to 10 km in terms of latitude and longitude, respectively. The spatial resolution of the tomographic region grid division, i.e., the spatial resolution of the interpolated grid, is 0.06°E×8 (8 grids in the longitude direction) and 0.05°N×7 (7 grids in the latitude direction). The height layers are 0.2km×10, 0.3km×8, 0.4km×6, 0.6km×4, ​​and 0.8km×1, resulting in 29 height layers in the elevation direction.

[0059] To compare the advantages of the method of this invention, the traditional tomography method is referred to as the T method; the water vapor density calculated based on NCEP-GFS forecast data is referred to as the G1 method; and the method of this invention is referred to as the G2 method. The NCEP-GFS forecast data used (with a time resolution of 6 hours) is the forecast field of the previous moment, with a forecast step size of 6 hours, and the GNSS data is all data of the previous 10-minute time window.

[0060] In this embodiment, using radiosonde data RS as the true value, the water vapor densities obtained by the T method, G1 method, and G2 method at Doy 203 UTC 00:00 (the time of rainfall) and Doy 215 UTC 00:00 (the time of non-rainfall) are compared. Figure 2 As shown. By Figure 2 It can be seen that the water vapor density results obtained by method T show strong oscillations, while method G1 shows a large deviation in water vapor density values ​​at the middle and lower levels. However, the corrected method G2 shows strong consistency with radiosonde data. Statistical analysis shows that, compared with radiosonde data, the average root mean square error of methods T, G1, and G2 is 1.60 g / m². 3 2.00g / m 3 and 0.35g / m 3 Therefore, this invention can obtain high-precision three-dimensional water vapor density.

[0061] In this embodiment, to verify the accuracy of the chromatography results, the water vapor density obtained from ERA5 data was used as a reference value. The root mean square error (RMSE) of the water vapor density compared to the reference value was calculated using the T method, G1 method, and G2 method. Figure 3 and Figure 4 As shown, the T method yields lower accuracy in obtaining water vapor density at lower levels, while the G2 method shows better improvement in lower-level water vapor density. Statistical analysis shows that during periods of rainfall, the RMSE of the G2 method is 1.02 g / m³. 3 Compared to methods T and G1, which have 1.99 g / m³ and 1.45 g / m³ respectively. 3 All results showed improvement. During non-rainy periods, the RMSE of the G2 method was 1.05 g / m³. 3 Compared to 2.20 g / m² in methods T and G1 3 and 1.57g / m 3 All of these improvements are evident. This indicates that the method of the present invention has high applicability under different weather conditions.

[0062] In this embodiment, to further verify the superiority of the proposed method, the computation time of each method is compared. All methods are computed under the same conditions. Figure 5 A comparison of the computation time of the method of this invention and existing methods shows that the computation time of the T method, G1 method, and G2 method at 0:00 and 12:00 on a single day is approximately 108s, 9s, and 16s, respectively. The G2 method shows a significant improvement in computation efficiency, which further confirms the superiority of the method of this invention.

[0063] In summary, the method of this invention obtains the average correction coefficient of water vapor density based on SWV residuals, corrects the grids that are traversed by GNSS signal rays using a single weight coefficient or a combination of weight coefficient and observation confidence, and then corrects the grids that are not traversed by GNSS signal rays. The corrected water vapor density values ​​are then superimposed on the original water vapor density background field to obtain a high-precision water vapor density field, thereby achieving rapid water vapor tomography and improving the accuracy of subsequent weather forecasts.

[0064] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Any simple modifications, alterations, or equivalent structural changes made to the above embodiments based on the technical essence of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A rapid water vapor tomography method integrating real-time BeiDou observations and background field constraints, characterized in that, The method includes the following steps: Step 1: Real-time atmospheric precipitable water volume (PWV) retrieved based on GNSS inversion; Step 2: Obtain grid data from NCEP-GFS forecast data and interpolate to obtain the water vapor density background field; Step 3: Obtain the SWV residual based on the atmospheric precipitable water volume (PWV) and water vapor density background field; Step 4: Obtain the average correction coefficient of water vapor density corresponding to the GNSS signal rays based on the SWV residual; Step 5: Perform water vapor density weight correction based on the average water vapor density correction coefficient to obtain the water vapor density correction value for each grid: Step 501: Determine the number of GNSS signal rays passing through each grid in the water vapor density background field. If the number of GNSS signal rays passing through the grid is 1, proceed to step 502; if the number of GNSS signal rays passing through the grid is greater than 1, proceed to step 503; if the grid is not passed through by GNSS signal rays, proceed to step 504. Step 502: Correct the average water vapor density correction coefficient with weighting coefficient to obtain the water vapor density correction value; Step 503: Correct the average correction coefficient of water vapor density by applying weighting coefficients and observation confidence to obtain the corrected value of water vapor density; Step 504: Correct the water vapor density of the grid that was not pierced by the GNSS signal rays based on the water vapor density correction value of the grid that was pierced by the GNSS signal rays, and obtain the water vapor density correction value of the grid that was not pierced by the GNSS signal rays. Step 6: Superimpose the water vapor density correction value onto the original water vapor density background field to obtain a high-precision water vapor density field and achieve rapid water vapor chromatography.

2. The rapid water vapor tomography method according to claim 1, which integrates real-time BeiDou observations and background field constraints, is characterized in that: The specific process of step two is as follows: Step 201: Obtain the three-dimensional grid data of the tomographic region from the NCEP-GFS forecast data; Step 202: Obtain the thermodynamic temperature and vapor pressure of each grid from the three-dimensional grid data of the tomographic region; Step 203, according to The water vapor density of each grid is obtained. ;in, The water vapor pressure of the grid. The thermodynamic temperature of the grid. The constant for water vapor; Step 204: Interpolate the water vapor density of multiple grids using an exponential function in the elevation direction and bilinear interpolation in the latitude and longitude directions to obtain the spatially interpolated water vapor density background field.

3. A rapid water vapor tomography method according to claim 2, which integrates real-time BeiDou observations and background field constraints, is characterized in that: The specific process of step 204 is as follows: Step 2041: Record the grid in the latitude and longitude direction as the horizontal grid, and the grid in the elevation direction as the vertical grid; Step 2042: Denote the multiple vertical grids with the same latitude and longitude as the i-th horizontal grid as the i-th group of vertical grids; where i is a positive integer; Step 2043: Obtain the water vapor density of the i-th vertical grid from the water vapor density of multiple grids; Step 2044: Construct the i-th exponential function ρ i(Z) =ρ i0 ×exp(-Z / H); where ρ i(Z) To determine the water vapor density at a location with the same latitude and longitude as the i-th horizontal grid and an elevation of Z, ρ i0 For the water vapor density at the same latitude and longitude as the i-th horizontal grid and at zero elevation, exp(·) is a power operation with the natural constant e as the base, Z is the elevation, and H is the water vapor elevation parameter; Step 2045: Based on the water vapor density of the i-th set of vertical grids, in ρ i(Z) Given Z, the least squares method is used to fit the i-th exponential function to obtain the water vapor density ρ at the same latitude and longitude as the i-th horizontal grid and at zero elevation. i0 Given the water vapor elevation parameter H, the i-th exponential function after fitting is obtained; Step 2046: Based on the fitted i-th exponential function, at the latitude and longitude of the i-th horizontal grid, interpolate along the elevation direction by inputting different elevations to obtain the water vapor density of different height layers. Step 2047: For the water vapor density at the same altitude layer, use bilinear interpolation in the latitude and longitude directions to obtain the spatially interpolated water vapor density background field.

4. A rapid water vapor tomography method according to claim 1, which integrates real-time BeiDou observations and background field constraints, is characterized in that: The specific process of step three is as follows: Step 301: Denote the grids through which any GNSS signal ray passes as the 1st grid, ..., the nth grid, ..., the Nth grid; where n and N are both positive integers, 1 ≤ n ≤ N; Step 302: Obtain the water vapor density of the nth grid from the water vapor density background field. and according to The oblique path water vapor content of the GNSS signal ray based on the water vapor density background field was obtained. ;in, The length of the GNSS signal ray passing through the nth grid; Step 303: Obtain the oblique path water vapor content (SWV) based on the atmospheric precipitable water vapor (PWV) inversion. Step 304, according to =SWV- SWV residuals are obtained. ; Step 305: Following the methods in steps 301 to 304, obtain the SWV residuals corresponding to multiple GNSS signal rays.

5. A rapid water vapor tomography method according to claim 4, which integrates real-time BeiDou observations and background field constraints, is characterized in that: The specific process of step four is as follows: Step 401, according to The weighting coefficients of the GNSS signal ray passing through the nth grid are obtained. ; Step 402: Construct the SWV residual equation The average correction coefficient for water vapor density corresponding to the GNSS signal rays was obtained. ; Step 403: Following the methods in steps 401 to 402, obtain the average correction coefficient of water vapor density corresponding to multiple GNSS signal rays.

6. A rapid water vapor tomography method according to claim 5, which integrates real-time BeiDou observation and background field constraints, is characterized in that: The specific process of step 502 is as follows: Step 5021: Following the method in step 401, obtain the weight coefficients of the grid. The average water vapor density correction coefficients for the GNSS signal rays passing through the grid are obtained from the average water vapor density correction coefficients corresponding to multiple GNSS signal rays and denoted as [missing value]. ; Step 5022, according to The water vapor density correction value of the grid is obtained. .

7. A rapid water vapor tomography method according to claim 6, which integrates real-time BeiDou observation and background field constraints, is characterized in that: The specific process of step 503 is as follows: Step 5031: If the number of GNSS signal rays passing through the grid is greater than 1, then the GNSS signal rays passing through the grid are respectively denoted as the 1st GNSS signal ray, ..., the mth GNSS signal ray, ..., the Mth GNSS signal ray; where M is a positive integer greater than or equal to 2, and 1≤m≤M; Step 5032, according to The elevation angle uncertainty component corresponding to the m-th GNSS signal ray is obtained. ;in, Let m be the elevation angle of the m-th GNSS signal ray; Step 5033, according to The carrier-to-noise ratio uncertainty component corresponding to the m-th GNSS signal ray is obtained. ;in, The signal-to-carrier-to-noise ratio of the m-th GNSS signal ray; Step 5034, according to The residual uncertainty component corresponding to the m-th GNSS signal ray is obtained. ;in, Let m be the SWV residual corresponding to the m-th GNSS signal ray. The median of the SWV residuals. The median of the absolute values ​​of the SWV residuals. It is the absolute value; Step 5035, according to The overall uncertainty corresponding to the m-th GNSS signal ray is obtained. ; Step 5036, according to The observation confidence level corresponding to the m-th GNSS signal ray is obtained. ;in, For adjustment coefficients; Step 5037, according to The water vapor density correction value obtained when the m-th GNSS signal ray passes through the grid is obtained. ;in, This is the average correction factor for the water vapor density of the grid as it is traversed by the m-th GNSS signal ray; The weighting coefficient corresponding to the m-th GNSS signal ray passing through the grid. Step 5038: If the observation confidence level corresponding to all GNSS signal rays is less than the preset threshold, then abandon the observation confidence level correction, and proceed according to... The water vapor density correction value of the grid is obtained. Otherwise, according to The water vapor density correction value of the grid is obtained. .

8. A rapid water vapor tomography method according to claim 7, which integrates real-time BeiDou observation and background field constraints, is characterized in that: The specific process of step 504 is as follows: Step 5041: Obtain the water vapor density correction value of the grid that is traversed by the GNSS signal ray in the height layer where the grid is located as a sample point. If the number of sample points is greater than or equal to 6, proceed to steps 5042 to 5044. If the number of sample points is greater than or equal to 3 but less than 6, proceed to step 5045; if the number of sample points is less than 3, the water vapor density correction value is zero. Step 5042: Construct a bivariate quadratic polynomial regression equation ρ=Ax 2 +By 2 +Cxy+Dx+Ey+F; where A is the first fitting coefficient, B is the second fitting coefficient, C is the third fitting coefficient, D is the fourth fitting coefficient, E is the fifth fitting coefficient, and F is the sixth fitting coefficient; x is the longitude independent variable, y is the latitude independent variable; and ρ is the water vapor density correction value dependent variable. Step 5043: Using the least squares method, based on the sample points in step 5041, fit and solve the bivariate quadratic polynomial regression equation to obtain the first to sixth fitting coefficients. Step 5044: Substitute the longitude and latitude of the grid that were not traversed by the GNSS signal rays in the same altitude layer into the fitted bivariate quadratic polynomial regression equation to obtain the water vapor density correction value of the grid. Step 5045: Construct a bivariate linear polynomial regression equation ρ=A´x+B´y+C´, and obtain the water vapor density correction value of the grid according to the methods in steps 5043 and 5044; where A´ is the first fitting parameter, B´ is the second fitting parameter, and C´ is the third fitting parameter.

Images