Adaptive partitioned water vapor tomography method and system based on multi-scale Bayesian prior

Abstract

The present invention relates to the field of meteorological remote sensing technology, and in particular to a method and system for adaptive partitioning of water vapor tomography based on a multi-scale Bayesian prior. The method comprises the following steps: S1, acquiring and processing multi-source observation data for the area to be analyzed; S2, constructing a multi-scale grid for the tomographic inversion region; S3, establishing a multi-scale representation of water vapor density and constructing a Bayesian prior constraint model; S4, constructing a three-dimensional tomographic inversion equation system; and S5, performing Bayesian inversion solution. The method integrates global navigation satellite system observation data, meteorological station data, and sounding station data, constructs Bayesian prior constraints, and employs a multi-scale representation and adaptive partitioning method to effectively reduce the number of parameters and improve the stability, accuracy, and computational efficiency of the tomographic inversion.

Description

Technical Field

[0001] The present invention relates to the field of meteorological remote sensing technology, and in particular to an adaptive partitioned water vapor tomography method and system based on multi-scale Bayesian prior. Background Art

[0002] As one of the most abundant greenhouse gases in Earth's atmosphere, atmospheric water vapor plays a crucial role in the development of weather systems and climate change. Accurately acquiring three-dimensional information on the distribution of atmospheric water vapor has important scientific and practical significance.

[0003] Traditional methods for obtaining atmospheric water vapor distribution primarily include sounding balloons, microwave radiometers, and meteorological satellites. However, these methods have limitations in terms of temporal and spatial resolution, observation range, and continuity. In recent years, the technology for retrieving atmospheric water vapor distribution has rapidly developed by exploiting the delay effect caused by atmospheric water vapor during the propagation of global navigation satellite system signals.

[0004] However, due to the limited distribution of global navigation satellite system sites and the restricted signal path geometry, the coefficient matrix A of the observation equations is often ill-conditioned. Especially in the vertical direction, because the signal path mainly passes through the atmosphere in an inclined state and lacks a horizontal signal path, the resolution and accuracy in the vertical direction are significantly lower than those in the horizontal direction. In order to improve the spatial resolution, more and smaller voxels need to be divided, resulting in a sharp increase in the number of unknown parameters and a rapid increase in computational complexity, which seriously affects the real-time or quasi-real-time processing capabilities. In particular, in large-scale high-resolution inversion, the computational burden of traditional methods is often unacceptable. The distribution of atmospheric water vapor has significant spatiotemporal variations, especially in severe convective weather systems, where the water vapor distribution may experience drastic local variations. Traditional uniform grid division methods make it difficult to simultaneously take into account both computational efficiency and local fine description capabilities. Summary of the Invention

[0005] (1) Technical issues to be resolved

[0006] The main purpose of the present invention is to provide an adaptive partitioned water vapor tomography method and system based on multi-scale Bayesian prior to solve the above technical problems.

[0007] (2) Technical solution

[0008] In order to achieve the above object, the present invention provides an adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior, comprising the steps of:

[0009] S1, acquiring and processing multi-source observation data of the area to be analyzed, wherein the multi-source observation data includes: satellite data, meteorological data, and sounding data;

[0010] S2, construct a multi-scale grid in the tomographic inversion area:

[0011] S21, constructing a three-dimensional water vapor density field by spatial interpolation based on the sounding data and the meteorological data , calculate the water vapor variation index between adjacent altitude layers ; The observation geometry intensity index calculated based on the path distribution in the satellite data ; and calculate the composite index :

[0012] ;

[0013] in, represents the weight coefficient;

[0014] S22, based on the comprehensive index The value of , divides the area to be analyzed into high-resolution areas and low-resolution areas ; Wherein, in the high resolution area A fine grid is used in the low-resolution area. A coarse grid is used in the

[0015] S3, establish a multi-scale representation of water vapor density and construct a Bayesian prior constraint model:

[0016] S31, based on the high-resolution area and low-resolution areas , establish a multi-scale representation of the three-dimensional water vapor density:

[0017] ;

[0018] in, Indicates position in space Water vapor density at Represents the multi-scale basis function, where the high-resolution region Using fine-scale basis functions, low-resolution areas Using coarse-scale basis functions; represents the basis function expansion coefficient;

[0019] S32, constructing a water vapor density vector to be inverted based on the sounding data and the meteorological data The prior probability distribution of :

[0020] ;

[0021] in, represents the prior water vapor density estimate, represents the prior covariance matrix, , represents the water vapor density vector to be inverted in the i-th grid;

[0022] S4, construct the three-dimensional tomographic inversion equations:

[0023] ;

[0024] in, represents the oblique path water vapor content observation vector, represents the coefficient matrix, represents the water vapor density vector to be inverted, represents the GNSS observation noise vector; represents the water vapor density observation vector at the weather station, represents the coefficient matrix corresponding to the weather station data, represents the observation noise vector of the weather station data; represents the water vapor density observation vector at the sounding station, represents the coefficient matrix corresponding to the sounding data, represents the observation noise vector of the sounding data;

[0025] S5, Bayesian inversion solution, obtains the three-dimensional water vapor density distribution:

[0026] Based on the three-dimensional tomographic inversion equations, an objective function is constructed, and by introducing sparse regularization and optimization algorithms, the optimal water vapor density vector β is iteratively solved to obtain a three-dimensional water vapor density distribution, wherein the objective function includes:

[0027] ;

[0028] in, represents the observation error covariance matrix, represents the regularization parameter, Represents a differential operator matrix based on the grid structure divided in step S2.

[0029] Preferably, the step S1 includes:

[0030] S11, obtaining satellite data obtained by a ground receiver observing the area to be analyzed, obtaining meteorological data obtained by a meteorological station observing the area to be analyzed, and obtaining sounding data obtained by a sounding station observing the area to be analyzed;

[0031] S12, based on the meteorological data and the sounding data, extracting the temperature data, pressure data and humidity data near the ground and at each altitude layer in the area to be analyzed, and constructing a near-surface water vapor density model Vertical profile water vapor density model ;

[0032] S13, based on the satellite data, calculating the total zenith delay through precise single point positioning, and calculating the zenith static delay in combination with the meteorological data, and extracting the slant path water vapor content data.

[0033] Preferably, the step S11 further comprises: performing quality control on the multi-source observation data obtained for the area to be analyzed, and eliminating abnormal values.

[0034] Preferably, the near-surface water vapor density model in step S12 for:

[0035] ;

[0036] Where, e represents the water vapor pressure, T represents the absolute temperature, Characteristic parameters of water vapor;

[0037] The vertical profile water vapor density model include:

[0038] ;

[0039] in, represents the water vapor density at height h, represents the surface water vapor density; represents the ground height, and H represents the water vapor height.

[0040] Preferably, the step S13 includes:

[0041] Based on the satellite data, a precise point positioning method is used to calculate the zenith total delay;

[0042] Combined with the surface pressure data provided by the weather station, the zenith static delay is calculated using the Saastamoinen model;

[0043] The zenith wet delay is obtained by the difference between the zenith total delay and the zenith static delay, and converted into the slant path wet delay;

[0044] The slant path wet delay is converted into slant path water vapor content data based on the mapping function and conversion factor.

[0045] Preferably, the step S22 includes:

[0046] The K-means clustering method is used to cluster the comprehensive index Classify and identify high-resolution areas and low-resolution areas Threshold :

[0047] , ;

[0048] in, The value of is determined by maximizing the ratio of between-class variance to within-class variance.

[0049] Preferably, in step S31, the fine-scale basis function is a wavelet basis function:

[0050] ;

[0051] in, represents the wavelet mother function, represents the grid point coordinates of the high-resolution area, Indicates the grid spacing;

[0052] The coarse-scale basis function is a cubic B-spline basis function:

[0053] ;

[0054] in, represents the cubic B-spline function, Represents the grid point coordinates of the low-resolution region, Indicates the grid spacing.

[0055] Preferably, the prior covariance matrix in step S32 is The construction method includes: based on the sounding data, calculating the matrix elements according to the following function:

[0056] ;

[0057] in, Indicates the The standard deviation of the basis function coefficients, Indicates the The standard deviation of the basis function coefficients, Indicates the relative lengths in the horizontal and vertical directions.

[0058] Preferably, the step S4, constructing a three-dimensional tomography inversion equation package, includes:

[0059] S41, construct the observation equations of the global navigation satellite system: establish the relationship between the water vapor content and water vapor density of the oblique path to form the observation equation group:

[0060] ;

[0061] in, represents the oblique path water vapor content observation vector, represents the coefficient matrix, represents the water vapor density vector to be inverted, represents the GNSS observation noise vector;

[0062] S42, constructing a multi-source data constraint equation group: integrating the meteorological data to construct a near-surface water vapor density constraint equation group:

[0063] ;

[0064] in, represents the water vapor density observation vector at the weather station, represents the coefficient matrix corresponding to the weather station data, represents the observation noise vector of the weather station data;

[0065] S43, integrating the sounding data to construct a vertical profile constraint equation group:

[0066] ;

[0067] in, represents the water vapor density observation vector at the sounding station, represents the coefficient matrix corresponding to the sounding data, represents the observation noise vector of the sounding data;

[0068] S44, combining the equations of steps S41 to S43 to obtain a three-dimensional tomography inversion equation group:

[0069] .

[0070] Preferably, in step S41, establishing the relationship between the water vapor content and the water vapor density of the oblique path includes:

[0071] S411, establishing an integral relationship between the water vapor content and water vapor density of the oblique path: the water vapor content of the oblique path With the three-dimensional water vapor density satisfy:

[0072] ;

[0073] in, represents the kth GNSS signal ray path, is the position parameter on the ray;

[0074] S412, based on the water vapor content of the segmented oblique path being equal to the corresponding water vapor density value multiplied by the segment length, a discretized observation equation is obtained:

[0075] ;

[0076] in, represents the intercept of the k-th ray path in the i-th grid, Represents the water vapor density vector to be inverted in the i-th grid.

[0077] Preferably, the optimization algorithm of step S5 is an alternating direction multiplier method, including the iterative formula:

[0078] ;

[0079] ;

[0080] ;

[0081] in, represents the penalty parameter, represents the soft threshold operator, k represents the number of iterations, and is dynamically adjusted based on the residual change and value.

[0082] The present invention also provides an adaptive partitioned water vapor tomography system based on multi-scale Bayesian prior, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the steps of any one of the above-mentioned adaptive partitioned water vapor tomography methods based on multi-scale Bayesian prior are implemented.

[0083] (3) Beneficial effects

[0084] This application proposes an adaptive partitioned water vapor tomography method and system based on multi-scale Bayesian priors. By constructing a comprehensive index, this method achieves adaptive regional partitioning, rationally allocating computing resources and significantly reducing the overall number of parameters and computational complexity while ensuring high resolution in key areas. Furthermore, the multi-scale basis functions enhance the ability to depict local details while keeping the number of parameters manageable, particularly in weather conditions with drastic water vapor variations. Introducing prior knowledge into the inversion process in the form of a probability distribution provides a theoretical basis for addressing ill-posed problems. Prior constraints within the Bayesian framework not only improve the stability of the inversion but also effectively reduce the uncertainty of the inversion results. By employing the alternating direction method of multipliers (ADMM) in conjunction with a dynamic parameter adjustment strategy, the algorithm achieves efficient and robust convergence performance. By automatically adjusting the regularization and penalty parameters using the residual ratio, the algorithm possesses powerful adaptability, enabling it to adapt to varying observation conditions and water vapor distribution characteristics. BRIEF DESCRIPTION OF THE DRAWINGS

[0085] Figure 1 A schematic flow chart of an adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior is provided in this embodiment.

[0086] Figure 2 A schematic diagram of adaptive grid division of an adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior provided in this embodiment.

[0087] Figure 3 A schematic diagram of ray tracing and voxel traversal of a global navigation satellite system signal provided in this embodiment.

[0088] Figure 4 Schematic diagram comparing water vapor density profiles at different altitudes using different methods provided in this embodiment.

[0089] Figure 5 A schematic diagram of the hardware structure of an adaptive partitioned water vapor tomography system based on multi-scale Bayesian priors provided in this embodiment. DETAILED DESCRIPTION

[0090] In order to better explain the present invention and facilitate understanding, the present invention is described in detail below through specific implementation methods in conjunction with the accompanying drawings.

[0091] It should be noted that all directional indications (such as up, down, left, right, front, back, etc.) in the embodiments of the present invention are only used to explain the relative position relationship, movement status, etc. between the various components under a certain specific posture (as shown in the accompanying drawings). If the specific posture changes, the directional indication will also change accordingly.

[0092] In addition, in the present invention, the terms "first," "second," etc. are used for descriptive purposes only and should not be understood as indicating or implying relative importance or implicitly specifying the number of the technical features indicated. Therefore, features specified as "first" or "second" may explicitly or implicitly include at least one of such features. In the description of the present invention, "plurality" means at least two, such as two, three, etc., unless otherwise specifically defined.

[0093] In the present invention, unless otherwise specified or limited, the terms "connection" and "fixation" should be understood in a broad sense. For example, "fixation" can refer to fixed connection, detachable connection, or integration; "connection" can refer to mechanical connection or electrical connection; it can refer to direct connection or indirect connection through an intermediate medium; it can refer to internal communication between two elements or interaction between two elements, unless otherwise specified. For those skilled in the art, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0094] like Figure 1 As shown, this embodiment provides an adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior, which includes steps: S1 to S5.

[0095] S1, acquiring and processing multi-source observation data of the area to be analyzed, wherein the multi-source observation data includes: satellite data, meteorological data, and sounding data; specifically including:

[0096] S11, obtaining satellite data of the area to be analyzed obtained by a ground receiver of a global navigation satellite system (GNSS), meteorological data of the area to be analyzed obtained by a meteorological station, and sounding data of the area to be analyzed obtained by a sounding station. Specifically, the observation data of the area to be analyzed obtained by the global information navigation satellite system receiver are received, the observation data of the area to be analyzed obtained by the meteorological station are obtained by obtaining data of a meteorological station co-located with or adjacent to the global navigation satellite system station, and the observation data of the area to be analyzed obtained by the sounding station are obtained by releasing a sounding balloon in the area to be analyzed.

[0097] S12, based on the meteorological data and the sounding data, extracting the temperature data, pressure data and humidity data of the near-surface and each altitude layer in the area to be analyzed, and constructing a near-surface water vapor density model , providing the near-surface water vapor density at the weather station location to provide a priori constraints on horizontal distribution; and constructing a vertical profile water vapor density model , providing the vertical water vapor distribution at the radiosonde station location to provide a priori constraints on the vertical profile;

[0098] S13, based on the satellite data, calculate the zenith total delay through precise single point positioning, calculate the zenith static delay in combination with the meteorological data, and extract the slant path water vapor content (SWV) data.

[0099] S2: Construct a multi-scale grid for the tomographic inversion region and implement adaptive meshing. This process dynamically divides the study area into grids of different resolutions (high-resolution and low-resolution regions) based on the intensity of water vapor variations and the density of the observation geometry within the study area, optimizing computational efficiency and inversion accuracy. This process specifically includes steps S21 and S22.

[0100] S21, based on the sounding data and the meteorological data, and the near-surface water vapor density model obtained in step S12 Vertical profile water vapor density model , constructing a three-dimensional water vapor density field through spatial interpolation , calculate the water vapor variation index between adjacent altitude layers and the observation geometry strength index calculated based on the national navigation satellite system signal path distribution processed in step S1 , calculate the composite index :

[0101] ;

[0102] in, Represents the weight coefficient (0≤α≤1), the weight coefficient The choice determines whether to place more emphasis on water vapor changes (α close to 1) or observed density (α close to 0); It represents the rate of change calculated from the prior water vapor value, reflecting the spatiotemporal variability of water vapor in the region; Represents the number of GNSS signals passing through a single grid. The denser the rays, the stronger the observation constraints in the area.

[0103] Specifically, the water vapor change index for: ;in, represents the water vapor density gradient operator, || || represents the L2 norm of the vector, Represents the three-dimensional water vapor density field.

[0104] Optionally, the observation geometry intensity index The calculation includes: based on the observation geometry of the global navigation positioning system processed in step S1, determining the three-dimensional voxels that the signal passes through by ray tracing, and calculating the observation geometry intensity index :

[0105] ;

[0106] in, Indicates whether the i-th ray passes through the unit voxel where the point (x, y, z) is located. If it passes, it is 1, otherwise it is 0; Represents the weight factor, which is related to factors such as satellite signal elevation angle and observation time.

[0107] S22, based on the comprehensive index The value of , divides the area to be analyzed into high-resolution areas and low-resolution areas ; Wherein, in the high resolution area A fine grid is used in the low-resolution area. A coarse grid is used.

[0108] The high-resolution area This is the composite index These areas have either drastic water vapor changes or dense observation rays (or both), and require a fine grid for detailed characterization.

[0109] The low-resolution area This is the composite index In these areas, the water vapor is relatively stable and / or the observation rays are sparse. A coarse grid can meet the requirements, greatly reducing the number of parameters to be estimated and improving computational efficiency.

[0110] In other embodiments, the step S22 divides the area to be analyzed into high-resolution areas by setting a threshold. and low-resolution areas .

[0111] S3, establishing a multi-scale representation of water vapor density and constructing a Bayesian prior constraint model: Using the background information from the sounding data and the meteorological data, an initial estimate of the spatial distribution of water vapor density is constructed. Multi-scale basis functions are used to represent the water vapor field, providing important constraints for the subsequent Bayesian inversion. Because the spatial distribution of water vapor is complex and multi-scale (e.g., varying dramatically in some areas and more gently in others), the use of a multi-scale representation method can better adapt to this characteristic and improve the accuracy of the subsequent water vapor tomography inversion. This specifically includes steps S31 and S32.

[0112] S31, based on the high-resolution area and low-resolution areas , establish a multi-scale representation of the three-dimensional water vapor density:

[0113] ;

[0114] in, Indicates position in space Water vapor density at represents the spatial position vector, ; Represents a multi-scale basis function, which is used to describe the spatial variation characteristics of water vapor density at different scales; among them, the high-resolution area Using fine-scale basis functions, low-resolution areas Using coarse-scale basis functions; Represents the basis function expansion coefficient, which is the coefficient to be determined, and N represents the total number of basis functions.

[0115] S32, constructing a water vapor density vector based on the sounding data and meteorological data The prior probability distribution of :

[0116] ;

[0117] in, represents the estimated value of the prior water vapor density obtained based on step S1, represents the prior covariance matrix based on the data statistics of step S1, , Represents the water vapor density vector to be inverted in the i-th grid.

[0118] Specifically, the Based on the sounding data and the meteorological data in the basis function The water vapor density prior estimate obtained by projection is used to calculate the initial value of the water vapor density vector β using the existing, relatively sparse but reliable sounding station and meteorological station data. .

[0119] S4: Constructing a three-dimensional tomographic inversion equation set, establishing a mathematical relationship between the multi-source observation data and the water vapor density vector β, integrating all available data into a unified model, and constructing a three-dimensional tomographic inversion equation set to provide a basis for subsequent inversion calculations. This specifically includes steps S41 to S43.

[0120] S41, constructing the global navigation satellite system observation equations: Based on the grid structure divided in step S2 and the multi-scale basis functions defined in step S3, the relationship between the oblique path water vapor content and the water vapor density is established to form the observation equation group:

[0121] ;

[0122] in, represents the oblique path water vapor content observation vector, represents the coefficient matrix of GNSS observations, represents the water vapor density vector to be inverted as defined in step S3, It represents the GNSS observation noise vector, which represents the various errors in the observation process.

[0123] Specifically, establishing the relationship between the water vapor content and the water vapor density of the oblique path includes:

[0124] S411, establishing an integral relationship between the water vapor content and water vapor density of the oblique path: the water vapor content of the oblique path With the three-dimensional water vapor density satisfy:

[0125] ;

[0126] in, represents the kth GNSS signal ray path, is the position parameter on the ray, Indicates water vapor density .

[0127] S412, based on the water vapor content SWV of the oblique path of each segment being equal to the corresponding water vapor density value multiplied by the segment length, a discretized observation equation is obtained:

[0128] ;

[0129] in, represents the intercept of the k-th ray path in the i-th grid, Represents the water vapor density vector to be inverted in the i-th grid.

[0130] S42, constructing a multi-source data constraint equation group: Based on the step S1, the meteorological data is integrated to construct a near-surface water vapor density constraint equation group:

[0131] ;

[0132] in, represents the water vapor density observation vector at the weather station, represents the coefficient matrix corresponding to the weather station data, represents the observation noise vector of the weather station data;

[0133] S43, based on step S1, the sounding data are integrated to construct a set of vertical profile constraint equations:

[0134] ;

[0135] in, represents the water vapor density observation vector at the sounding station, represents the coefficient matrix corresponding to the sounding data, represents the observation noise vector of the sounding data.

[0136] S44, combining the equations of steps S41 to S43 to obtain a three-dimensional tomography inversion equation group:

[0137] .

[0138] S5, Bayesian inversion solution, to obtain the three-dimensional water vapor density distribution: Based on the three-dimensional tomographic inversion equations, combined with the prior water vapor density estimate and the prior covariance matrix Based on the differential operator matrix of the grid structure divided in step S2, an objective function is constructed, and by introducing sparse regularization and optimization algorithms, the optimal water vapor density vector β is iteratively solved to obtain a three-dimensional water vapor density distribution, wherein the objective function includes:

[0139] ;

[0140] in, represents the observation error covariance matrix, represents the regularization parameter, The difference operator matrix representing the grid structure divided based on the step S2 acts on the vector β representing the water vapor density at the discrete grid points to calculate its spatial difference, which is used to impose spatial smoothness or structural regularization constraints.

[0141] Furthermore, the near-surface water vapor density model in step S12 include:

[0142] ;

[0143] Where, e represents the water vapor pressure (hPa), T represents the absolute temperature (K), Indicates the characteristic parameters of water vapor, specifically the water vapor ratio gas constant, which does not change with the environment. ;

[0144] The vertical profile water vapor density model include:

[0145] ;

[0146] in, represents the water vapor density at height h, represents the surface water vapor density; represents the ground height, and H represents the water vapor height.

[0147] Specifically, the surface water vapor density It is the reference parameter in the sounding model, indicating the height Water vapor density at the surface, the surface water vapor density The data source can be set according to the actual situation, for example, it can be extracted from the ground layer observation value of the sounding data, or when the sounding station and the meteorological station are located at the same time, the ground water vapor density The data source can come from the near-surface water vapor density calculated by meteorological stations.

[0148] Preferably, the step S13 comprises: calculating the Zenith Total Delay (ZTD) by using the Precise Point Positioning (PPP) method based on the observation data of the global navigation satellite system for the area to be analyzed; and calculating the Zenith Hydrostatic Delay (ZHD) by using the Saastamoinen model in combination with the ground pressure data provided by the meteorological station:

[0149] ;

[0150] in, is the latitude of the measuring station (°), H is the geodetic height of the measuring station (km), and P is the surface air pressure of the measuring station (hPa). The Saastamoinen model is a model widely used in geodesy and meteorology to calculate the zenith total delay (ZTD).

[0151] The zenith wet delay (ZWD) is obtained from the difference between the zenith total delay and the zenith static delay: ;

[0152] Convert the Zenith Wet Delay (ZWD) to Slant Path Wet Delay (SWD):

[0153] ;

[0154] in, The mapping function of wet delay ZWD is represented by is the atmospheric horizontal gradient mapping function, Represents the gradient of the satellite signal in the north-south direction, It represents the gradient of satellite signal in the east-west direction, e is the satellite altitude angle, is the satellite azimuth.

[0155] The slant path wet delay is converted to slant path water vapor content (SWV) data based on the mapping function and conversion factor:

[0156] ;

[0157] in, Represents the conversion factor, which is calculated as follows: , is the water vapor gas constant, is the atmospheric refraction constant.

[0158] Specifically, step S11 further includes: performing quality control on the multi-source observation data obtained for the area to be analyzed to remove outliers, and also includes data checking and missing value interpolation on the observation data of the area to be analyzed by the meteorological station, which is conducive to improving the inversion accuracy.

[0159] Furthermore, the step S22 includes:

[0160] The K-means clustering method is used to cluster the comprehensive index Classify and identify high-resolution areas and low-resolution areas Threshold :

[0161] , ;

[0162] in, The value of is determined by maximizing the ratio of between-class variance to within-class variance.

[0163] Optionally, in step S31, the fine-scale basis function is a wavelet basis function:

[0164] ;

[0165] in, represents the wavelet mother function, represents the grid point coordinates of the high-resolution area, Indicates the grid spacing.

[0166] The coarse-scale basis function is a cubic B-spline basis function:

[0167] ;

[0168] in, represents the cubic B-spline function, Represents the grid point coordinates of the low-resolution region, Indicates the grid spacing.

[0169] Furthermore, the prior covariance matrix in step S32 is constructed by calculating the matrix elements according to the following function based on the sounding station observation data obtained in step S1:

[0170] ;

[0171] in, Indicates the The standard deviation of the basis function coefficients, Indicates the The standard deviation of the basis function coefficients, Indicates the relative lengths in the horizontal and vertical directions.

[0172] As a preferred embodiment of the present invention, the optimization algorithm of step S5 is an alternating direction multiplier method, including the iterative formula:

[0173] ;

[0174] ;

[0175] ;

[0176] in, represents the penalty parameter, represents the soft threshold operator, k represents the number of iterations, and is dynamically adjusted based on the residual change and value; it reduces the amount of calculation while improving the inversion accuracy of complex water vapor structures.

[0177] As a preferred embodiment of the present invention, step S5 further includes: setting convergence conditions and convergence judgment: introducing convergence conditions in the iterative solution process: when the relative improvement of M consecutive iterations is less than the threshold Stop iteration when , and this condition is met M times continuously or when the number of iterations reaches the preset maximum value The iteration is stopped when the threshold For example, it can be set to arrive The number of iterations M can be set to 3 to 5, for example, with a preset maximum value of Typically set to 50 to 100.

[0178] Furthermore, the observation error covariance matrix The following adaptive estimation method is used: Initialization: Calculate the initial error covariance of each type of observation data: ;

[0179] in, represents the initial error covariance of the global navigation positioning system observation data, represents the initial error covariance of the weather station observation data, represents the initial error covariance of the observation data at the radiosonde station;

[0180] Iterative update: During the iteration process, the residual error between the slant path water vapor content of the kth inversion and the slant path water vapor content observation is updated:

[0181] ;

[0182] in, is the residual value of the kth iteration; is the smoothing factor, and its value range is [0,1].

[0183] The following describes an adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior in this embodiment through a specific operation process.

[0184] Step 1: Observation data processing and quality control: Process the observation data of CORS stations in a certain area: Use the precise single point positioning method to calculate the zenith total delay (ZTD) of each station using satellite signal data.

[0185] Combined with weather station data, the Saastamoinen model is used to calculate the zenith static delay (ZHD) and thus the zenith wet delay (ZWD).

[0186] Utilization of sounding data: In the region, the sounding station releases sounding balloons twice a day (at 00UTC and 12UTC) to obtain data such as air pressure, temperature and relative humidity at different altitudes.

[0187] The water vapor density at each height of the sounding station is calculated using the formula:

[0188] ;

[0189] Where e is the water vapor pressure (hPa) and T is the absolute temperature (K).

[0190] All observation data are then quality-checked to remove obvious outliers and ensure data reliability. This includes: for global navigation satellite system observation data, satellite signals with elevation angles below 10° are removed; for weather station data, temperature, air pressure, and relative humidity are checked to see if they are within a reasonable range; and for sounding data, abnormal temperature and humidity values ​​are removed.

[0191] Step 2: Figure 2 As shown in Figure 2, the study area is divided into high-resolution areas according to the water vapor distribution characteristics and observation geometry intensity. and low-resolution areas .

[0192] Construct the initial water vapor density field by interpolating sounding data and meteorological data , calculate the water vapor change index in the study area , this index reflects the degree of change in water vapor distribution, and the calculation formula is .

[0193] Calculate the observation geometry intensity index : Counts the number of GNSS signal rays that pass through each grid (i.e., the number of times each grid is traversed by a ray; the greater the number of traversals, the stronger the observation constraint). This index reflects the ability of the observation data to constrain each grid.

[0194] like Figure 3 As shown in Figure 2, the GNSS signal propagates from the satellite to the receiver, forming a ray path that passes through different grid cells, providing observation constraints for each grid. Then, the composite index is calculated. , where the weight coefficient The threshold is set to 0.6, which has achieved good results in this embodiment. K-means clustering is used to divide the study area into high-resolution areas and low-resolution areas, and the threshold is determined by maximizing the ratio of the between-class variance to the within-class variance. .

[0195] Step 3: Establish a multi-scale representation of water vapor density and a Bayesian prior constraint model. In this step, a multi-scale representation of water vapor density is established and a Bayesian prior constraint model is constructed. In the multi-scale representation, the water vapor density distribution is expressed as .

[0196] For low-resolution areas , using cubic B-spline basis function; for high-resolution areas , using wavelet basis functions to better characterize local details. Based on the historical data of the sounding station, the prior probability distribution is constructed:

[0197] ;

[0198] in is the prior water vapor density estimate, calculated from the sounding data. The prior covariance matrix It reflects the correlation of water vapor density at different locations, and its elements are calculated as

[0199] .

[0200] In this example, the parameters of the prior distribution, horizontal correlation length, and = =50 km, vertical correlation length =1 km, these parameters can better describe the spatial distribution characteristics of water vapor in the area.

[0201] Step 4: Establish a 3D water vapor tomography inversion model. This step establishes the integral relationship between the slant path water vapor content (SWV) and the water vapor density, discretizing the continuous integral relationship into a system of linear equations. First, establish the GNSS observation equations. Ray tracing is used to calculate the length of each signal path through each grid cell and construct the coefficient matrix A.

[0202] Ray tracing uses an accurate three-dimensional ray bending model, taking into account the influence of atmospheric refractive index changes with altitude. Then, by fusing weather station data and sounding data, near-surface water vapor density constraints and vertical profile constraints are constructed to form a complete three-dimensional tomographic inversion equation set:

[0203] .

[0204] Step 5: Adaptive iterative solution and parameter update: Under the Bayesian framework, the objective function is constructed and solved using the alternating direction multiplier method (ADMM). The objective function includes:

[0205] .

[0206] Initial regularization parameter Set to 0.1, the penalty parameter Set to 1. During the iteration process, these parameters are dynamically adjusted based on the residual change. When the relative improvement of 3 consecutive iterations is less than Or stop the iteration when the number of iterations reaches 100.

[0207] like Figure 4 Figure 2 shows the water vapor density profile at 00 UTC on August 26, 2019, at the radiosonde station in the region, compared with the inversion results using the traditional method and the proposed method. The inversion results using the proposed method (solid red line) are closer to the measured values ​​at the radiosonde station (solid black line), especially at altitudes of 2–3 km and 6–8 km, where the inversion accuracy is significantly better than that of the traditional method (dashed blue line).

[0208] like Figure 5 As shown, the present application also provides an adaptive partitioned water vapor tomography system based on multi-scale Bayesian prior, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the steps of an adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior as described in any one of the above items are implemented.

[0209] Figure 5 FIG. 1 is a schematic diagram of the hardware structure for running an adaptive partitioned water vapor tomography method based on multi-scale Bayesian priors provided by an embodiment of the present invention. Figure 5 As shown, this embodiment / computer 6 includes: a processor 60, a memory 61, and a computer program 62 stored in the memory 61 and executable on the processor 60, such as a program for executing an adaptive partitioned water vapor tomography method based on a multi-scale Bayesian prior. When the processor 60 executes the computer program 62, it implements the steps described above in the embodiments for executing an adaptive partitioned water vapor tomography method based on a multi-scale Bayesian prior. Alternatively, when the processor 60 executes the computer program 62, it implements the functions of the modules / units described above in the various device embodiments.

[0210] For example, the computer program 62 may be divided into one or more modules / units, which are stored in the memory 61 and executed by the processor 60 to implement the present invention. The one or more modules / units may be a series of computer program instruction segments capable of implementing specific functions, and the instruction segments are used to describe the execution process of the computer program 62 in the computer 6.

[0211] The computer 6 can be a computing device such as a desktop computer, a notebook, a PDA, or a cloud server. The computer 6 device can include, but is not limited to, a processor 60 and a memory 61. It can be understood by those skilled in the art that Figure 5 This is only an example of computer 6 and does not constitute a limitation on computer 6. It may include more or fewer components than shown in the figure, or a combination of certain components, or different components. For example, the computer 6 may also include input and output devices, network access devices, buses, etc.

[0212] The processor 60 may be a central processing unit (CPU), other general-purpose processors, digital signal processors (DSP), application-specific integrated circuits (ASIC), field-programmable gate arrays (FPGA), other programmable logic devices, discrete gate or transistor logic devices, discrete hardware components, etc. The general-purpose processor may be a microprocessor or any conventional processor.

[0213] The memory 61 can be an internal storage unit of the computer 6, such as the computer's hard drive or memory. Alternatively, the memory 61 can be an external storage device of the computer 6, such as a plug-in hard drive, a Smart Media Card (SMC), a Secure Digital (SD) card, a flash memory card, etc., equipped on the terminal device. Furthermore, the memory 61 can include both the internal storage unit of the computer 6 and an external storage device. The memory 61 is used to store the computer program and other programs and data required by the terminal device. The memory 61 can also be used to temporarily store data that has been output or is about to be output.

[0214] Those skilled in the art can clearly understand that, for the convenience and brevity of description, only the division of the above-mentioned functional units and modules is used as an example for illustration. In actual applications, the above-mentioned functions can be distributed and completed by different functional units and modules as needed, that is, the internal structure of the device can be divided into different functional units or modules to complete all or part of the functions described above. The functional units and modules in the embodiment can be integrated into one processing unit, or each unit can exist physically alone, or two or more units can be integrated into one unit. The above-mentioned integrated unit can be implemented in the form of hardware or in the form of software functional units. In addition, the specific names of the functional units and modules are only for the convenience of distinguishing each other, and are not used to limit the scope of protection of this application. The specific working process of the units and modules in the above-mentioned system can refer to the corresponding process in the aforementioned method embodiment, and will not be repeated here.

[0215] In the above embodiments, the description of each embodiment has its own focus. For parts that are not described or recorded in detail in a certain embodiment, reference can be made to the relevant description of other embodiments.

[0216] Those skilled in the art will appreciate that the units and algorithm steps of each example described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are performed in hardware or software depends on the specific application and design constraints of the technical solution. Professionals and technicians can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of the present invention.

[0217] In the embodiments provided by the present invention, it should be understood that the disclosed devices / terminal equipment and methods can be implemented in other ways. For example, the device / terminal equipment embodiments described above are merely illustrative. For example, the division of the modules or units is merely a logical function division. In actual implementation, there may be other division methods, such as multiple units or components can be combined or integrated into another system, or some features can be ignored or not executed. Another point is that the mutual coupling or direct coupling or communication connection shown or discussed can be through some interfaces, indirect coupling or communication connection of devices or units, which can be electrical, mechanical or other forms.

[0218] The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed across multiple network units. Some or all of these units may be selected to achieve the purpose of this embodiment according to actual needs.

[0219] In addition, the functional units in the various embodiments of the present invention may be integrated into a single processing unit, each unit may exist physically separately, or two or more units may be integrated into a single unit. The aforementioned integrated units may be implemented in the form of hardware or software functional units.

[0220] If the integrated module / unit is implemented as a software functional unit and sold or used as a standalone product, it can be stored in a computer-readable storage medium. Based on this understanding, the present invention can implement all or part of the process steps in the above-mentioned method embodiments by using a computer program to instruct the relevant hardware. The computer program can be stored in a computer-readable storage medium. When executed by a processor, the computer program can implement the steps of each of the above-mentioned method embodiments. The computer program includes computer program code, which can be in source code form, object code form, executable file, or some intermediate form. The computer-readable medium can include: any entity or device capable of carrying the computer program code, recording medium, USB flash drive, mobile hard drive, magnetic disk, optical disk, computer memory, read-only memory (ROM), random access memory (RAM), electric carrier signal, telecommunication signal, and software distribution medium. It should be noted that the content of the computer-readable medium can be appropriately increased or decreased based on the requirements of legislation and patent practice in a jurisdiction. For example, in some jurisdictions, based on legislation and patent practice, computer-readable media does not include electric carrier signals and telecommunication signals.

[0221] The above are merely specific application examples of the present invention and do not constitute any limitation on the scope of protection of the present invention. In addition to the above embodiments, the present invention may also have other implementations. Any technical solution formed by equivalent replacement or equivalent transformation falls within the scope of protection claimed by the present invention.

Claims

1. An adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior, characterized by: The following steps are involved: S1, acquiring and processing multi-source observation data of the area to be analyzed, wherein the multi-source observation data includes: satellite data, meteorological data, and sounding data; S2, construct a multi-scale grid in the tomographic inversion area: S21, constructing a three-dimensional water vapor density field by spatial interpolation based on the sounding data and the meteorological data , calculate the water vapor variation index between adjacent altitude layers ; The observation geometry intensity index calculated based on the path distribution in the satellite data ; and calculate the composite index : ； in, represents the weight coefficient; S22, based on the comprehensive index The value of , divides the area to be analyzed into high-resolution areas and low-resolution areas ; Wherein, in the high resolution area A fine grid is used in the low-resolution area. A coarse grid is used in the S3, establish a multi-scale representation of water vapor density and construct a Bayesian prior constraint model: S31, based on the high-resolution area and low-resolution areas , establish a multi-scale representation of the three-dimensional water vapor density: ； in, Indicates position in space Water vapor density at Represents the multi-scale basis function, where the high-resolution region Using fine-scale basis functions, low-resolution areas Using coarse-scale basis functions; represents the basis function expansion coefficient; S32, constructing a water vapor density vector to be inverted based on the sounding data and the meteorological data The prior probability distribution of : ； in, represents the prior water vapor density estimate, represents the prior covariance matrix, , represents the water vapor density vector to be inverted in the i-th grid; S4, construct the three-dimensional tomographic inversion equations: ； in, represents the oblique path water vapor content observation vector, represents the coefficient matrix, represents the water vapor density vector to be inverted, represents the GNSS observation noise vector; represents the water vapor density observation vector at the weather station, represents the coefficient matrix corresponding to the weather station data, represents the observation noise vector of the weather station data; represents the water vapor density observation vector at the sounding station, represents the coefficient matrix corresponding to the sounding data, represents the observation noise vector of the sounding data; S5, Bayesian inversion solution, obtains the three-dimensional water vapor density distribution: Based on the three-dimensional tomographic inversion equations, an objective function is constructed, and by introducing sparse regularization and optimization algorithms, the optimal water vapor density vector β is iteratively solved to obtain a three-dimensional water vapor density distribution, wherein the objective function includes: ； in, represents the observation error covariance matrix, represents the regularization parameter, Represents a differential operator matrix based on the grid structure divided in step S2.

2. The adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior according to claim 1 is characterized in that: The step S1 comprises: S11, obtaining satellite data obtained by a ground receiver observing the area to be analyzed, obtaining meteorological data obtained by a meteorological station observing the area to be analyzed, and obtaining sounding data obtained by a sounding station observing the area to be analyzed; S12, based on the meteorological data and the sounding data, extracting the temperature data, pressure data and humidity data near the ground and at each altitude layer in the area to be analyzed, and constructing a near-surface water vapor density model Vertical profile water vapor density model ; S13, based on the satellite data, calculating the total zenith delay through precise single point positioning, and calculating the zenith static delay in combination with the meteorological data, and extracting the slant path water vapor content data.

3. The adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior according to claim 2 is characterized in that: The step S11 further includes: performing quality control on the acquired multi-source observation data of the area to be analyzed to eliminate outliers; The near-surface water vapor density model in step S12 for: ； Where, e represents the water vapor pressure, T represents the absolute temperature, Characteristic parameters of water vapor; The vertical profile water vapor density model include: ； in, represents the water vapor density at height h, represents the surface water vapor density; Indicates the ground height, H indicates the water vapor height; The step S13 includes: Based on the satellite data, a precise point positioning method is used to calculate the zenith total delay; Combined with the surface pressure data provided by the weather station, the zenith static delay is calculated using the Saastamoinen model; The zenith wet delay is obtained by the difference between the zenith total delay and the zenith static delay, and converted into the slant path wet delay; The slant path wet delay is converted into slant path water vapor content data based on the mapping function and conversion factor.

4. The adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior according to claim 1 is characterized in that: The step S22 includes: The K-means clustering method is used to cluster the comprehensive index Classify and identify high-resolution areas and low-resolution areas Threshold : ， ； in, The value of is determined by maximizing the ratio of between-class variance to within-class variance.

5. The adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior according to claim 1 is characterized in that: In step S31, the fine-scale basis function is a wavelet basis function: ； in, represents the wavelet mother function, represents the grid point coordinates of the high-resolution area, Indicates the grid spacing; The coarse-scale basis function is a cubic B-spline basis function: ； in, represents the cubic B-spline function, Represents the grid point coordinates of the low-resolution region, Indicates the grid spacing.

6. The adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior according to claim 5, characterized in that: The prior covariance matrix in step S32 The construction method includes: based on the sounding data, calculating the matrix elements according to the following function: ； in, Indicates the The standard deviation of the basis function coefficients, Indicates the The standard deviation of the basis function coefficients, Indicates the relative lengths in the horizontal and vertical directions.

7. The adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior according to claim 1 is characterized in that: The step S4, constructing a three-dimensional tomography inversion equation package, includes: S41, construct the observation equations of the global navigation satellite system: establish the relationship between the water vapor content and water vapor density of the oblique path to form the observation equation group: ； in, represents the oblique path water vapor content observation vector, represents the coefficient matrix, represents the water vapor density vector to be inverted, represents the GNSS observation noise vector; S42, constructing a multi-source data constraint equation group: integrating the meteorological data to construct a near-surface water vapor density constraint equation group: ； in, represents the water vapor density observation vector at the weather station, represents the coefficient matrix corresponding to the weather station data, represents the observation noise vector of the weather station data; S43, integrating the sounding data to construct a vertical profile constraint equation group: ； in, represents the water vapor density observation vector at the sounding station, represents the coefficient matrix corresponding to the sounding data, represents the observation noise vector of the sounding data; S44, combining the equations of steps S41 to S43 to obtain a three-dimensional tomography inversion equation group: 。 8. The adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior according to claim 7 is characterized in that: In step S41, establishing the relationship between the water vapor content and the water vapor density of the oblique path includes: S411, establishing an integral relationship between the water vapor content and water vapor density of the oblique path: the water vapor content of the oblique path With the three-dimensional water vapor density satisfy: ； in, represents the kth GNSS signal ray path, is the position parameter on the ray; S412, based on the water vapor content of the segmented oblique path being equal to the corresponding water vapor density value multiplied by the segment length, a discretized observation equation is obtained: ； in, Represents the intercept of the k-th ray path in the i-th grid.

9. The adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior according to claim 1, characterized in that: The optimization algorithm of step S5 is the alternating direction multiplier method, which includes the iterative formula: ； ； ； in, represents the penalty parameter, represents the soft threshold operator, k represents the number of iterations, and is dynamically adjusted based on the residual change and value.

10. An adaptive partitioned water vapor tomography system based on multi-scale Bayesian prior, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, the steps of the adaptive partitioned water vapor tomography method based on multi-scale Bayesian prior are implemented as described in any one of claims 1 to 9.

Images