Water vapor chromatography method, system and equipment based on machine learning and medium

Abstract

The invention provides a water vapor chromatography method, system and device based on machine learning, and a medium, and the method comprises the following steps: constructing a historical training sample set, and employing the training gradient of the historical training sample set to promote a decision tree model, and generating an inclined path wet delay prediction model; establishing a virtual site coordinate in a blank area which does not cover a global navigation satellite system site, constructing a virtual physical feature vector in combination with a digital elevation model and meteorological reanalysis data, and obtaining a virtual oblique path wet delay observation value through the prediction model; and combining a real observation value of a global navigation satellite system station with the virtual oblique path wet delay observation value, establishing an enhanced chromatography equation set, solving by adopting an algebraic reconstruction algorithm, and performing inversion to generate a three-dimensional water vapor density distribution field of the target monitoring area. By implementing the technical scheme provided by the invention, high-precision virtual observation data can be generated in an area with sparse observation stations or uneven signal path distribution, and the stability and spatial resolution of a three-dimensional water vapor inversion result are improved.

Description

Technical Field

[0001] This invention relates to the field of atmospheric remote sensing technology, and in particular to a water vapor tomography method, system, device and medium based on machine learning. Background Technology

[0002] With the rapid development of atmospheric remote sensing technology and the demand for numerical weather prediction, the monitoring accuracy requirements for the three-dimensional distribution of tropospheric water vapor are continuously increasing, and the requirements for the spatiotemporal resolution of data for severe weather warnings are constantly improving. How to use the Global Navigation Satellite System (GNSS) to acquire high-precision three-dimensional water vapor field information around the clock and at low cost has become a core requirement for improving the accuracy of numerical weather prediction and the ability to prevent meteorological disasters.

[0003] In existing technologies, to address the ill-conditioned tomographic inversion problem caused by uneven distribution of GNSS observation signal paths and sparse station locations, mathematical smoothing constraints or the construction of traditional virtual observations are typically introduced. For example, empirical models and mapping functions are used to estimate virtual slant path delays in different azimuths based on the station's zenith wet delay, or horizontal and vertical smoothing constraints are directly applied during equation solving. However, in practical applications, existing technologies are mostly based on the assumption of water vapor isotropy or empirical models, and virtual observation data generated under complex atmospheric conditions may introduce systematic biases. In scenarios with sparse or geometrically uneven observation network distribution, this bias limits the improvement effect on the ill-conditioned nature of the tomographic equations, posing a risk of poor stability and low accuracy in water vapor inversion results. Summary of the Invention

[0004] In view of this, this application provides a water vapor chromatography method, system, device and medium based on machine learning to solve the above problems.

[0005] Firstly, a machine learning-based water vapor chromatography method is provided, which includes:

[0006] Construct a historical training sample set containing historical signal geometric parameters, historical station spatiotemporal parameters, and historical zenith wet delay values, and obtain the ground truth labels of the slant path wet delay corresponding to the historical training sample set;

[0007] The gradient boosting decision tree model is iteratively trained using historical training sample sets and ground truth labels for oblique path wet delay to generate an oblique path wet delay prediction model.

[0008] Based on the spatial distribution data of Global Navigation Satellite System (GNSS) stations within the target monitoring area, establish virtual station coordinates in blank areas not covered by GNSS stations;

[0009] For the virtual station coordinates, a virtual ray path is generated according to a preset azimuth angle sequence and a preset elevation angle sequence to obtain the virtual signal geometric parameters;

[0010] Acquire digital elevation model data and meteorological reanalysis data, extract the spatiotemporal parameters of virtual station coordinates based on digital elevation model data, and calculate the synthetic zenith wet delay value of virtual station coordinates based on meteorological reanalysis data.

[0011] A virtual physical feature vector is constructed based on the virtual signal geometric parameters, the virtual station spatiotemporal parameters, and the synthetic zenith wet delay value. The virtual physical feature vector is then input into the slant path wet delay prediction model to obtain the virtual slant path wet delay observation value.

[0012] An enhanced tomographic equation set was established using real observations from Global Navigation Satellite System (GNSS) stations and virtual slant path wet delay observations.

[0013] An algebraic reconstruction algorithm is used to solve the enhanced tomography equations and invert the three-dimensional water vapor density distribution field of the target monitoring area.

[0014] The above technical solution introduces a machine learning-based oblique path wet delay prediction model. Building upon traditional GNSS water vapor tomography, it trains a predictive model using historical observation data to reflect true atmospheric characteristics, thereby generating high-precision virtual oblique path wet delay observations in areas with sparse stations or uneven signal path distribution. This method does not rely on empirical models or isotropic assumptions; instead, it automatically learns the complex nonlinear relationship between atmospheric water vapor distribution and path delay through a data-driven approach, reducing systematic biases in virtual observations at the source. Consequently, the enhanced tomographic equations are effectively constrained in blank areas, significantly improving the ill-conditioned nature of the tomographic inversion equations and enhancing the stability and accuracy of the three-dimensional water vapor density distribution inversion results.

[0015] Optionally, the preset gradient boosting decision tree model is iteratively trained using historical training sample sets and ground truth labels for oblique path wet delays to generate an oblique path wet delay prediction model, specifically including:

[0016] The 24 hours are divided into multiple consecutive preset time windows according to preset time intervals;

[0017] Based on the time information contained in the spatiotemporal parameters of historical stations, the historical training sample set and the true value labels of the wet delay of the oblique path are filtered to the corresponding preset time window to obtain the time-segmented training data corresponding to each preset time window.

[0018] Multiple preset gradient boosting decision tree sub-models are trained using training data from each time period to obtain multiple trained gradient boosting decision tree models corresponding to each preset time window.

[0019] Establish a calling index relationship between each preset time window and the corresponding trained gradient boosting decision tree sub-model;

[0020] By utilizing the gradient boosting decision tree sub-models after training and the relationship between each call index, a wet delay prediction model for oblique paths is established.

[0021] The above technical solution, by dividing the all-weather time series into multiple preset time windows and training gradient boosting decision tree sub-models for different time periods, can capture the nonlinear temporal characteristics of atmospheric water vapor in the diurnal cycle and weather evolution. This time-segmented modeling mechanism enables the machine learning model to adaptively reflect the water vapor distribution pattern in different time periods, avoiding the prediction distortion problem of a uniform model under time-varying conditions. This further improves the temporal accuracy of virtual oblique path wet delay prediction and provides more stable and reliable virtual observation support for subsequent tomographic inversion.

[0022] Optionally, based on the spatial distribution data of Global Navigation Satellite System (GNSS) stations within the target monitoring area, virtual station coordinates are established in blank areas not covered by GNSS stations, specifically including:

[0023] The target monitoring area is divided into grids according to preset latitude and longitude intervals to obtain multiple geographic grid units;

[0024] Using the spatial distribution data of Global Navigation Satellite System (GNSS) stations within the target monitoring area, the inclusion relationship of each geographic grid cell is determined, and all blank geographic grid cells that do not contain GNSS stations are selected from each geographic grid cell.

[0025] Obtain the longitude and latitude of the geometric center point of each blank geographic grid cell, and establish each set of longitude and latitude as the coordinates of multiple virtual stations, with each set of longitude and latitude corresponding to one virtual station coordinate.

[0026] The above technical solution achieves uniform spatial filling of observation blind spots by dividing the target monitoring area into latitude and longitude grids and establishing virtual station coordinates based on the geometric centers of blank grid cells without GNSS station coverage. This method ensures the balanced and representative distribution of virtual stations within the region, enabling the generated virtual observation paths to effectively supplement the deficiencies of real observation paths in space. This significantly improves the constraints of the tomographic grid at edges and in blank areas, reducing spatial distortion and boundary discontinuities in the inversion results.

[0027] Optionally, the synthetic zenith wet delay value for the virtual station coordinates is calculated based on meteorological reanalysis data, specifically including:

[0028] Using meteorological reanalysis data, four-dimensional interpolation is performed on the longitude, latitude and time of the virtual station coordinates to construct vertical profiles of atmospheric parameters corresponding to multiple atmospheric layers. The vertical profiles of atmospheric parameters include the layer temperature, layer water vapor pressure and layer potential of each atmospheric layer.

[0029] Using a pre-defined empirical formula for wet refractive index, the wet refractive index of each atmospheric layer is calculated based on the stratification temperature and stratified water vapor pressure.

[0030] Calculate the difference between the potentials of adjacent layers in the vertical profile of atmospheric parameters, and convert the difference into the geometric thickness of each atmospheric layer based on gravitational acceleration;

[0031] Using the geographic elevation of the virtual site coordinates as the starting height for integration, the products of each wet refractive index and the corresponding geometric layer thickness are summed to obtain the synthetic zenith wet delay value.

[0032] The above technical solution reconstructs a vertical profile including stratified temperature, vapor pressure, and geopotential by performing four-dimensional interpolation on meteorological reanalysis data at a virtual site. It then uses the wet refractive index integral formula to calculate the synthetic zenith wet delay value, achieving an accurate physical representation of virtual observations. This method ensures that the zenith wet delay at the virtual site reflects the actual meteorological conditions at that time and place, significantly reducing delay bias caused by empirical model approximations, thereby improving the physical reliability and regional adaptability of virtual observation inputs.

[0033] Optionally, a virtual physical feature vector is constructed based on the virtual signal geometric parameters, the virtual station spatiotemporal parameters, and the synthesized zenith wet delay value, specifically including:

[0034] The satellite azimuth and elevation angles of the virtual ray path are determined as the geometric parameters of the virtual signal;

[0035] The geographic elevation corresponding to the virtual station coordinates is retrieved from the digital elevation model data, and the current normalized time is obtained. The virtual station coordinates, geographic elevation, and current normalized time are then determined as the spatiotemporal parameters of the virtual station.

[0036] The satellite azimuth angle, satellite elevation angle, virtual site coordinates, geographic elevation, current normalized time, and synthetic zenith wet delay value are combined into a virtual physical feature vector.

[0037] The above technical solution combines elements such as the azimuth angle, elevation angle, virtual station coordinates, geographic elevation, time, and synthetic zenith wet delay value of the virtual ray path into a virtual physical feature vector, which fully characterizes the geometric and physical properties of each virtual signal path. This feature vector structure is consistent with the input of real GNSS observation data, enabling the oblique path wet delay prediction model to perform accurate inference in the virtual scene. This ensures that the generated virtual observations are consistent with real observations in terms of distribution characteristics, improving the reliability and usability of virtual observations.

[0038] Optionally, an enhanced tomographic equation set can be established using real observations from Global Navigation Satellite System (GNSS) stations and virtual slant path wet delay observations, specifically including:

[0039] The real slope path wet delay vector corresponding to the real observation value and the virtual slope path wet delay vector corresponding to the virtual slope path wet delay observation value are vertically stacked to generate the enhanced observation vector.

[0040] Based on the spatial geometric relationship between the real signal path and the virtual ray path corresponding to the real observation value and the preset tomographic grid, calculate the real intercept length of the real signal path in the tomographic grid and the virtual intercept length of the virtual ray path in the tomographic grid.

[0041] A real signal path design matrix is ​​constructed based on the real intercept length, and a virtual signal path design matrix is ​​constructed based on the virtual intercept length.

[0042] The real signal path design matrix and the virtual signal path design matrix are vertically stacked to generate an enhanced design matrix, and a linear observation equation is established by combining the enhanced observation vector.

[0043] The above technical solution establishes an enhanced tomographic equation set by stacking real and virtual observations and their corresponding design matrices. This allows for the simultaneous use of real and virtual supplementary constraints during the inversion process. This enhanced equation set significantly increases the number and independence of equations, improves the matrix rank deficiency problem, and ensures sufficient constraints for tomographic inversion even in sparse observation regions. This effectively reduces the ill-posedness of the inversion and enhances the stability of the three-dimensional water vapor distribution estimation.

[0044] Optionally, an algebraic reconstruction algorithm is used to solve the enhanced tomography equations to invert and generate a three-dimensional water vapor density distribution field in the target monitoring area, specifically including:

[0045] Interpolation is performed on each grid node of the tomographic grid based on meteorological reanalysis data to generate a three-dimensional wet refractive index field;

[0046] The discrete values ​​in the three-dimensional wet refractive index field are arranged into a vector to be solved, and the product of the enhancement design matrix and the vector to be solved is calculated to obtain the forward modeling observation vector.

[0047] Calculate the residual vector between the forward modeling observation vector and the enhanced observation vector, and iteratively correct the discrete values ​​in the vector to be solved based on the residual vector until the residual vector satisfies the preset convergence condition;

[0048] Based on meteorological reanalysis data, the temperature parameters of each grid node are determined, and combined with the physical conversion formula between wet refractive index and water vapor density, the converged solution vector is converted into a three-dimensional water vapor density distribution field.

[0049] The above technical solution establishes an initial wet refractive index field based on meteorological reanalysis data and uses an algebraic reconstruction algorithm to iteratively correct the solution vector under residual constraints, thereby obtaining a convergent and stable wet refractive index distribution. Combined with the physical conversion formula between wet refractive index and water vapor density, it achieves accurate inversion from the refractive index field to a three-dimensional water vapor density field. This method not only improves the numerical convergence and physical consistency of the inversion results but also ensures that the obtained water vapor density field outperforms traditional regularized smoothing methods in both spatial resolution and quantitative accuracy.

[0050] Secondly, a machine learning-based water vapor chromatography system is provided, the system comprising:

[0051] The sample construction module is configured to construct a historical training sample set containing historical signal geometric parameters, historical station spatiotemporal parameters, and historical zenith wet delay values, and to obtain the true value labels of the slant path wet delay corresponding to the historical training sample set.

[0052] The model training module is configured to iteratively train a preset gradient boosting decision tree model using historical training sample sets and ground truth labels for oblique path wet delay, thereby generating an oblique path wet delay prediction model.

[0053] The virtual station establishment module is configured to establish virtual station coordinates in blank areas that do not cover global navigation satellite system stations, based on the spatial distribution data of global navigation satellite system stations within the target monitoring area.

[0054] The path generation module is configured to generate virtual ray paths for virtual station coordinates according to a preset azimuth sequence and a preset elevation sequence, thereby obtaining the geometric parameters of the virtual signal.

[0055] The parameter calculation module is configured to acquire digital elevation model data and meteorological reanalysis data, extract the spatiotemporal parameters of the virtual station coordinates based on the digital elevation model data, and calculate the synthetic zenith wet delay value of the virtual station coordinates based on the meteorological reanalysis data.

[0056] The delay prediction module is configured to construct a virtual physical feature vector based on the virtual signal geometric parameters, the virtual station spatiotemporal parameters, and the synthetic zenith wet delay value, and input the virtual physical feature vector into the slant path wet delay prediction model to obtain the virtual slant path wet delay observation value.

[0057] The equation-building module is configured to build an enhanced tomographic equation set using real observations from Global Navigation Satellite System (GNSS) sites and virtual slant path wet delay observations.

[0058] The tomographic inversion module is configured to solve the enhanced tomographic equations using an algebraic reconstruction algorithm, and invert to generate a three-dimensional water vapor density distribution field of the target monitoring area.

[0059] Thirdly, an electronic device is provided, including a processor, a memory, a user interface, and a network interface, wherein the memory is used to store instructions, the user interface and the network interface are both used to communicate with other devices, and the processor is used to execute the instructions stored in the memory to cause the electronic device to perform the method as described in any of the above.

[0060] Fourthly, a computer-readable storage medium is provided, the computer-readable storage medium storing instructions that, when executed, perform the method as described in any of the preceding claims.

[0061] In summary, implementing one or more technical solutions provided in this application has at least the following technical effects or advantages:

[0062] By introducing a data-driven virtual observation generation mechanism into the traditional GNSS water vapor tomography framework, a deep fusion of machine learning models and physical tomography models is achieved. This fusion mechanism enables the tomographic inversion process to form a self-consistent closed loop across three levels: data acquisition, model constraints, and result inversion. It can still generate physically consistent supplementary observation information even when observational data is insufficient; during the tomographic solution stage, enhanced observational constraints can be used to improve equation stability and inversion convergence speed; and at the output level, it can significantly improve the spatial continuity and temporal consistency of the three-dimensional water vapor distribution field. Compared to traditional methods that rely solely on mathematical smoothing or empirical virtual observations, this invention not only improves the quantitative accuracy of the inversion results but also enhances the algorithm's adaptability to complex meteorological conditions and dynamic water vapor structures, thus possessing higher engineering applicability and long-term operational reliability. Attached Figure Description

[0063] Figure 1 This is an exemplary system architecture diagram of a water vapor tomography method or a water vapor tomography system based on machine learning that applies this application;

[0064] Figure 2 This is a flowchart illustrating a machine learning-based water vapor chromatography method disclosed in this application.

[0065] Figure 3 This is a schematic diagram of a machine learning-based water vapor chromatography system disclosed in this application;

[0066] Figure 4 This is a schematic diagram of the structure of an electronic device disclosed in this application.

[0067] Figure reference numerals: 100, System architecture; 101, First terminal device; 102, Second terminal device; 103, Third terminal device; 104, Network; 105, Server; 301, Sample construction module; 302, Model training module; 303, Virtual station establishment module; 304, Path generation module; 305, Parameter calculation module; 306, Delay prediction module; 307, Equation establishment module; 308, Tomographic inversion module; 401, Processor; 402, Communication bus; 403, User interface; 404, Network interface; 405, Memory. Detailed Implementation

[0068] To enable those skilled in the art to better understand the technical solutions in this specification, the technical solutions in the embodiments of this specification will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments.

[0069] In the description of the embodiments of this application, the words "for example" or "for instance" are used to indicate examples, illustrations, or explanations. Any embodiment or design that is described as "for example" or "for instance" in the embodiments of this application should not be construed as being more preferred or advantageous than other embodiments or design options. Rather, the use of the words "for example" or "for instance" is intended to present the relevant concepts in a specific manner.

[0070] In the description of the embodiments of this application, the term "multiple" means two or more. For example, multiple systems means two or more systems, and multiple screen terminals means two or more screen terminals. Furthermore, the terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the indicated technical features. Thus, a feature defined with "first" or "second" may explicitly or implicitly include one or more of that feature. The terms "comprising," "including," "having," and variations thereof all mean "including but not limited to," unless otherwise specifically emphasized.

[0071] Figure 1 An exemplary system architecture diagram is shown, illustrating an embodiment of a machine learning-based water vapor chromatography method or a machine learning-based water vapor chromatography system to which this application can be applied.

[0072] like Figure 1As shown, the system architecture 100 may include a first terminal device 101, a second terminal device 102, a third terminal device 103, a network 104, and a server 105. The network 104 is used as a medium to provide communication links between the terminal devices 101, 102, 103, and the server 105. The network 104 may include various connection types, such as wired or wireless communication links or fiber optic cables, etc.

[0073] Users can use terminal devices 101, 102, and 103 to interact with server 105 via network 104 to receive or send messages, etc. Various communication client applications can be installed on terminal devices 101, 102, and 103, such as model training applications, video recognition applications, web browser applications, social platform software, etc.

[0074] Terminal devices 101, 102, and 103 can be either hardware or software. When terminal devices 101, 102, and 103 are hardware, they can be various electronic devices with displays, including but not limited to smartphones, tablets, e-book readers, MP3 (Moving Picture Experts Group Audio Layer III) players, MP4 (Moving Picture Experts Group Audio Layer IV) players, laptops, and desktop computers, etc. When terminal devices 101, 102, and 103 are software, they can be installed in the aforementioned electronic devices. They can be implemented as multiple software programs or software modules (e.g., multiple software programs or software modules used to provide distributed services) or as a single software program or software module. No specific limitations are imposed here.

[0075] When terminals 101, 102, and 103 are hardware devices, video capture devices can also be installed on them. These video capture devices can be various devices capable of capturing video, such as cameras, sensors, etc. Users can use the video capture devices on terminals 101, 102, and 103 to capture video.

[0076] Server 105 can be a server that provides various services, such as a backend server for processing data displayed on terminal devices 101, 102, and 103. The backend server can analyze and process the received data and can feed back the processing results (such as recognition results) to the terminal devices.

[0077] It should be noted that a server can be either hardware or software. When the server is hardware, it can be implemented as a distributed server cluster consisting of multiple servers, or as a single server. When the server is software, it can be implemented as multiple software programs or software modules (e.g., multiple software programs or software modules used to provide distributed services), or as a single software program or software module. No specific limitations are made here.

[0078] It should be understood that Figure 1 The number of terminal devices, networks, and servers shown is merely illustrative. Depending on implementation needs, any number of terminal devices, networks, and servers can be included. In particular, if the target data does not need to be obtained remotely, the above system architecture may exclude the network and include only terminal devices or servers.

[0079] Figure 2 This is a flowchart illustrating a machine learning-based water vapor chromatography method according to an embodiment of this application. This method can be implemented using a computer program, a microcontroller, or run on a machine learning-based water vapor chromatography system. The computer program can be integrated into an application or run as a standalone utility application. The specific steps of a machine learning-based water vapor chromatography method are described in detail below.

[0080] S201: Construct a historical training sample set containing historical signal geometric parameters, historical station spatiotemporal parameters, and historical zenith wet delay values, and obtain the ground truth labels of the slant path wet delay corresponding to the historical training sample set.

[0081] In this embodiment, the historical training sample set refers to a dataset pre-collected and organized for training a machine learning model. It contains a large amount of input feature data and corresponding target output data, used to guide the model in learning the mapping relationship between input and output. For example, this sample set can be a standardized dataset formed by preprocessing observation data from multiple reference stations in the target monitoring area over the past three years. Each sample data point contains information such as the geometric relationship between a specific time, a specific satellite, and a specific station, the station status, and atmospheric conditions. When constructing the historical training sample set, the 3σ criterion is used to perform quality control on the historical zenith wet delay (ZWD) sequence. The mean μ and standard deviation σ within the sliding window are calculated. If the observed value x at a certain time... i Satisfy |x iIf -μ|>3σ, it is considered a gross error and discarded. For the discarded missing epochs or epochs where the original observation was interrupted, if the number of consecutive missing epochs is less than 3, Lagrange interpolation is used to fill in the missing epochs using data from adjacent epochs; if the number of consecutive missing epochs exceeds 3, the data for that period is marked as invalid and not included in the sample set. Based on the sampling rate of GNSS observations (e.g., 30 seconds or 1 minute), meteorological parameters and ephemeris parameters are aligned to the same timestamp through linear interpolation.

[0082] Specifically, to build a high-precision prediction model, the system needs to retrieve long-period GNSS observation data and corresponding precise ephemeris products from the historical database. For historical signal geometric parameters, the system calculates the relative positional relationship between satellite coordinates and station coordinates at historical moments, extracting the azimuth and elevation angles of the satellite relative to the station. These two parameters describe the specific path direction of the signal through the atmosphere. For historical station spatiotemporal parameters, the system extracts the station's longitude, latitude, ellipsoidal height, and the annualized day or normalized time of the observation time. These parameters reflect the spatial heterogeneity and temporal periodicity of water vapor distribution. For historical zenith wet delay (ZWD) values, the system processes historical observation data using precise point positioning technology or high-precision baseline calculation software, separating the zenith static delay from the calculated total zenith delay to obtain high-precision zenith wet delay data as a key physical feature. While preparing the above input features, it is necessary to obtain the target for model learning, namely, the ground truth label of slant wet delay (SWD). Typically, high-precision GNSS post-processing software is used to project and restore the zenith wet delay onto each signal path direction based on the anisotropic mapping function and gradient parameters. The slant path wet delay value under that path is then calculated and used as the ground truth label. The historical signal geometric parameters, historical station spatiotemporal parameters, and historical zenith wet delay values ​​obtained above are aligned according to timestamps and satellite numbers, combined into a feature vector, and the corresponding slant path wet delay ground truth labels are matched one-to-one with this feature vector to construct a structured and standardized historical training sample set.

[0083] S202: Iteratively train the preset gradient boosting decision tree model using the historical training sample set and the ground truth labels of the oblique path wet delay to generate an oblique path wet delay prediction model.

[0084] For example, the constructed historical training sample set is divided into a training set and a validation set according to a preset ratio. During the actual model training process, to prevent overfitting and improve generalization ability, grid search or Bayesian optimization is used to optimize the key hyperparameters of the gradient boosting decision tree. Specifically, the main hyperparameters adjusted include: the learning rate, used to control the contribution of each tree to the combined model, typically set between 0.01 and 0.1 to balance convergence speed and accuracy; the maximum depth of the tree and the number of leaf nodes, limiting the complexity of the tree and preventing overfitting to historical noisy data; and the feature sampling rate and data sampling rate, randomly selecting some features and samples in each iteration to increase the diversity of the model.

[0085] Furthermore, the training process employs an early stopping strategy: training automatically stops when the loss function (such as root mean square error) on the validation set no longer decreases for N consecutive epochs (e.g., 50 epochs). The loss function L is typically defined as the sum of the root mean square error between the predicted slope wet delay and the ground truth label, plus a regularization term: L = Σ(SWD) pred -SWD true ) 2 +Ω(f). Where, SWD pred This represents the predicted slant path wet delay value output by the model; SWD true Ω(f) represents the ground truth label of the wet delay of the slant path in the historical data; Ω(f) represents the regularization term (such as L1 or L2 regularization), used to penalize overly complex tree structures to prevent overfitting. After training, the system also outputs the importance score of each feature. If a feature is found to have an extremely low gain contribution, it will be removed in subsequent model iterations to achieve feature dimensionality reduction. Finally, a prediction model that can accurately describe the complex nonlinear mapping relationship between multidimensional physical features and wet delay of the slant path is obtained.

[0086] In one possible implementation, a pre-defined gradient boosting decision tree model is iteratively trained using a historical training sample set and ground truth labels for oblique path wet delay to generate an oblique path wet delay prediction model. Specifically, this includes: dividing 24 hours into multiple consecutive pre-defined time windows according to a pre-defined time interval; filtering the historical training sample set and oblique path wet delay ground truth labels to the corresponding pre-defined time windows based on the time information contained in the historical station spatiotemporal parameters, obtaining time-segmented training data for each pre-defined time window; training multiple pre-defined gradient boosting decision tree sub-models using the time-segmented training data for each pre-defined time window, obtaining multiple trained gradient boosting decision tree sub-models corresponding to each pre-defined time window; establishing a call index relationship between each pre-defined time window and the corresponding trained gradient boosting decision tree model; and establishing an oblique path wet delay prediction model using each trained gradient boosting decision tree model and each call index relationship.

[0087] In the embodiments of this application, the gradient boosting decision tree sub-model refers to an independent functional component constituting the overall prediction system. Each sub-model is a complete and independently trained machine learning model used to establish a nonlinear mapping relationship between input features and oblique path wet delay within a specific time range. For example, if the day is divided into 24 time periods, 24 sub-models are generated accordingly. The 8th sub-model is specifically responsible for predicting the atmospheric wet delay conditions between 07:00 and 08:00 in the morning, thereby capturing the subtle features of water vapor changes over time.

[0088] Specifically, considering the significant diurnal variation periodicity of atmospheric water vapor distribution, the system divides the 24 hours of a day into multiple consecutive preset time windows according to preset time intervals (e.g., every 1 hour or every 2 hours) to improve the model's ability to capture water vapor changes at different times. It iterates through the historical station spatiotemporal parameters of each sample data point in the historical training sample set, reads the timestamp information, and determines the time period to which each sample belongs based on this time information. The historical training sample set and its corresponding oblique path wet delay ground truth labels are then precisely filtered and distributed to their respective preset time windows, thus obtaining independent time-segmented training data corresponding to each preset time window. Multiple preset gradient boosting decision tree sub-models are initialized, and supervised learning iterative training is performed on the corresponding time-segment sub-models using the time-segmented training data from each group. By continuously optimizing the loss function to minimize the prediction error, multiple trained gradient boosting decision tree sub-models corresponding one-to-one with each preset time window are obtained. To ensure accurate model matching based on input time in subsequent application phases, a call index relationship is established between each preset time window (as the index key) and the corresponding trained gradient boosting decision tree sub-model (as the index value). This relationship essentially defines a time-based model lookup table. The trained gradient boosting decision tree sub-models and the aforementioned call index relationship are logically encapsulated and integrated to establish a sloped path wet delay prediction model with all-weather, high-precision prediction capabilities.

[0089] S203: Based on the spatial distribution data of Global Navigation Satellite System (GNSS) stations within the target monitoring area, establish virtual station coordinates in blank areas not covered by GNSS stations.

[0090] For example, the system reads the geographical extent of the target monitoring area and the latitude and longitude positions of all existing Global Navigation Satellite System (GNSS) stations within the area, analyzing the sparsity and coverage blind spots of the current physical observation network. To meet the spatial density requirements of high-resolution water vapor tomography, the system plans several supplementary observation points in these blank areas without deployed physical receivers, based on a preset uniform distribution principle, and calculates their corresponding geographical coordinates. These coordinates are established as virtual stations, thereby logically expanding the discrete and irregular physical observation network into a uniformly covered observation network, providing a spatial reference for the subsequent generation of virtual observation data.

[0091] In one possible implementation, based on the spatial distribution data of Global Navigation Satellite System (GNSS) stations within the target monitoring area, virtual station coordinates are established in blank areas not covered by GNSS stations. Specifically, this includes: dividing the target monitoring area into grids according to preset latitude and longitude intervals to obtain multiple geographic grid units; using the spatial distribution data of GNSS stations within the target monitoring area to determine the inclusion relationship of each geographic grid unit, and filtering out all blank geographic grid units that do not contain GNSS stations from each geographic grid unit; obtaining the longitude and latitude of the geometric center point of each blank geographic grid unit, and establishing each set of longitude and latitude as multiple virtual station coordinates, with one set of longitude and latitude corresponding to one virtual station coordinate.

[0092] In this embodiment, virtual station coordinates refer to the location information of simulated observation points set by algorithmic logic in areas lacking physical observation equipment, used as a spatial reference for generating simulated atmospheric wet delay data. For example, in mountainous areas with complex terrain or vast lake areas, since it is impossible to deploy real global navigation satellite system receivers, the system will establish latitude and longitude points at specific locations (such as grid centers) in these areas, and these points are the virtual station coordinates.

[0093] Specifically, to effectively supplement observation data and improve the spatial resolution of the tomographic model, the system acquires the geographical boundaries of the target monitoring area and divides the target monitoring area into regular grids according to preset latitude and longitude intervals (e.g., 0.1 degrees or 0.25 degrees each in the longitude and latitude directions), thereby constructing multiple geographical grid units covering the entire monitoring range. The system imports precise coordinate data of Global Navigation Satellite System (GNSS) stations, i.e., spatial distribution data, and traverses each geographical grid unit. By comparing the numerical relationship between the station coordinates and the grid boundaries, it determines whether there are actual GNSS stations within the latitude and longitude range of the unit. Based on this inclusion relationship determination result, grids containing at least one actual station are marked as covered grids and removed, thus filtering out all blank geographical grid units that do not contain GNSS stations. For each filtered blank geographical grid unit, the system calculates the average value of its longitude boundary and the average value of its latitude boundary to obtain the longitude and latitude of the geometric center point of each blank geographical grid unit. These calculated geometric center points are used as supplementary observation points, and each set of longitude and latitude is established as multiple virtual station coordinates to ensure that one set of longitude and latitude corresponds to one virtual station coordinate, thereby constructing a uniform observation network that combines virtual and real data.

[0094] S204: For the virtual station coordinates, generate a virtual ray path according to the preset azimuth angle sequence and the preset elevation angle sequence to obtain the virtual signal geometric parameters.

[0095] In this embodiment, a virtual ray path refers to a signal propagation trajectory from a virtual ground station to a specific direction in the sky, simulated mathematically in the absence of real satellite signal coverage. This trajectory is used as an integration path in subsequent steps to detect atmospheric water vapor distribution. For example, if the system assumes a satellite is located 45 degrees northeast of the virtual station at an elevation angle of 30 degrees, then the straight line segment connecting the virtual station and this assumed location constitutes a virtual ray path.

[0096] Specifically, to ensure the acquisition of multi-angle coverage atmospheric sounding information even in areas lacking real observation data, the system uses the established virtual station coordinates as the geometric origin and loads pre-configured scanning parameters, including a preset azimuth sequence covering all directions (e.g., starting at 0 degrees and setting an azimuth every 30 degrees) and a preset elevation sequence covering different elevation angles (e.g., starting at a 10-degree cutoff elevation and setting an elevation every 15 degrees). For each virtual station coordinate, the program iterates and combines each angle value in the preset azimuth sequence with each angle value in the preset elevation sequence to construct a series of spatial vector lines pointing from the ground to the tropopause, thereby generating a virtual ray path covering the virtual station. After generating the virtual ray path, the system introduces a terrain occlusion verification mechanism. Based on the Digital Elevation Model (DEM), it checks whether the generated ray path passes through terrain obstacles. For any sampling point P on the ray path... k (x) k y k , z k ), query the corresponding coordinates (x) in the DEM k y k ) surface elevation H dem If there exists any point z satisfying k <H dem If the ray path is blocked by terrain (i.e., the signal is unreachable), the system will automatically discard the invalid path. The azimuth and elevation angle values ​​corresponding to these paths are extracted and recorded as key data describing the spatial morphology of the path, thus obtaining the virtual signal geometric parameters.

[0097] S205: Acquire digital elevation model data and meteorological reanalysis data, extract the spatiotemporal parameters of the virtual station coordinates based on the digital elevation model data, and calculate the synthetic zenith wet delay value of the virtual station coordinates based on the meteorological reanalysis data.

[0098] For example, the system indexes the digital elevation model based on the established latitude and longitude coordinates of the virtual station, accurately matches and extracts the geographic elevation information of that location, thereby completing the three-dimensional spatial and temporal attributes of the virtual station. Simultaneously, using meteorological reanalysis data as a large-scale meteorological background field, the system extracts the atmospheric temperature, humidity, and pressure parameters at the location and time of the virtual station, and calculates the total water vapor effect in the vertical direction at that location using a meteorological physical model, converting it into an equivalent synthetic zenith wet delay value. This value serves as the vertical physical reference for subsequently constructing virtual tilt path observations.

[0099] In one possible implementation, the synthetic zenith wet delay value of the virtual station coordinates is calculated based on meteorological reanalysis data. Specifically, this includes: using meteorological reanalysis data to perform four-dimensional interpolation on the longitude, latitude, and time of the virtual station coordinates to construct vertical profiles of atmospheric parameters corresponding to multiple atmospheric layers. The vertical profiles of atmospheric parameters include the layer temperature, layer vapor pressure, and layer potential of each atmospheric layer; using a preset empirical formula for wet refractive index, the wet refractive index of each atmospheric layer is calculated based on the layer temperature and layer vapor pressure; the difference between adjacent layer potentials in the vertical profiles of atmospheric parameters is calculated, and the difference is converted into the geometric layer thickness of each atmospheric layer based on gravitational acceleration; using the geographical elevation of the virtual station coordinates as the starting height for integration, the products of each wet refractive index and the corresponding geometric layer thickness are summed to obtain the synthetic zenith wet delay value.

[0100] In this embodiment, the atmospheric parameter vertical profile refers to a data sequence describing the distribution of atmospheric physical properties with altitude above a target location, used to finely reflect the atmospheric thermodynamic state at different levels from the ground to the upper atmosphere. For example, it includes a list of numerical values ​​for meteorological elements such as temperature, air pressure, and humidity at each potential altitude above a virtual station, and serves as the basic physical model for calculating the atmospheric refraction delay effect.

[0101] Specifically, to obtain a high-precision meteorological background field at the virtual station, the system first accesses the global meteorological model database to acquire high spatiotemporal resolution meteorological reanalysis data (such as ERA5 or NCEP data). Since this data is typically based on standard grid points, the system needs to resample the grid data in the spatiotemporal dimension using four-dimensional linear interpolation or spline interpolation algorithms, based on the precise longitude and latitude of the virtual station coordinates and the current observation time. This resampling constructs a vertical profile of atmospheric parameters corresponding to multiple atmospheric layers (such as isobaric layers from the bottom to the top of the troposphere) above the virtual station. This vertical profile details the stratified temperature, stratified vapor pressure, and stratified geopotential height for each atmospheric layer. Based on the theory of electromagnetic wave propagation in inhomogeneous media, and using a pre-defined empirical formula for wet refractive index (such as the Thayer or Smith-Weintraub formula), the stratified temperature and vapor pressure of each layer are substituted to calculate the wet refractive index N of each atmospheric layer. wet The formula is as follows: , where N wet k1 represents the wet refractive index of a given atmosphere (unit: dimensionless); e represents the stratification vapor pressure of that atmosphere (unit: hPa); T represents the stratification absolute temperature of that atmosphere (unit: K); k2' is the atmospheric refractive index constant, typically taken as approximately 22.1 K / hPa; k3 is the atmospheric refractive index constant, typically taken as approximately 373900 K. 2 / hPa.

[0102] Furthermore, to determine the physical thickness of each atmospheric layer, the difference between the potentials of adjacent layers in the vertical profile of atmospheric parameters is calculated, and local gravity parameters are introduced. Based on gravitational acceleration, the difference is converted into the geometric thickness of each atmospheric layer. The process of simulating the vertical traversal of the signal through the atmosphere is simulated. Using the geographic elevation of the virtual station coordinates as the starting height for integration, the product of each moist refractive index and its corresponding geometric thickness is accumulated and summed layer by layer (i.e., discrete numerical integration is performed), thereby obtaining the synthetic zenith wet delay (ZWD) value for the virtual station. , where N wet，i Δh represents the wet refractive index of the i-th atmospheric layer; i The geometric thickness of the i-th atmospheric layer (in meters) is represented by the potential difference between adjacent layers, converted based on gravitational acceleration.

[0103] S206: Construct a virtual physical feature vector based on the virtual signal geometric parameters, virtual station spatiotemporal parameters, and synthetic zenith wet delay value, and input the virtual physical feature vector into the slant path wet delay prediction model to obtain the virtual slant path wet delay observation value.

[0104] For example, the system standardizes and assembles the acquired multidimensional physical parameters (including spatial geometric information, station spatiotemporal state, and atmospheric background values) according to a preset feature format to form an input vector that can characterize the virtual observation scene. The system then calls a pre-trained oblique path wet delay prediction model to perform inference calculations on this input vector. Utilizing the nonlinear mapping relationship learned by the model, the vertical zenith wet delay is mapped and extended to the oblique path wet delay along a specific oblique path direction, thereby obtaining highly realistic virtual observation values.

[0105] In one possible implementation, a virtual physical feature vector is constructed based on the virtual signal geometric parameters, the virtual station spatiotemporal parameters, and the synthetic zenith wet delay value. Specifically, this includes: determining the satellite azimuth and satellite elevation angle of the virtual ray path as the virtual signal geometric parameters; retrieving the geographic elevation corresponding to the virtual station coordinates from the digital elevation model data and obtaining the current normalized time; determining the virtual station coordinates, geographic elevation, and current normalized time as the virtual station spatiotemporal parameters; and combining the satellite azimuth, satellite elevation angle, virtual station coordinates, geographic elevation, current normalized time, and synthetic zenith wet delay value into the virtual physical feature vector.

[0106] In this embodiment, the virtual physical feature vector refers to a multi-dimensional data combination constructed to meet the input requirements of a machine learning model. It encapsulates geometric features, geographical environment features, and meteorological background features along the simulated ray path, and is used to drive the model to predict the oblique path wet delay. For example, it is a standardized numerical array containing satellite azimuth, elevation angle, station latitude and longitude, elevation, time, and synthetic zenith wet delay value, which can comprehensively describe the physical state under the virtual observation scenario.

[0107] Specifically, to construct feature data that conforms to the preset model input dimensions, the system first parses the virtual ray path data generated in the previous steps, extracting the satellite azimuth and elevation angles describing spatial pointing, and defining these two angle values ​​as the virtual signal geometric parameters. Based on the virtual station's planar location, the system indexes and queries the digital elevation model data to retrieve the precise geographic elevation corresponding to the virtual station coordinates. Simultaneously, it reads the current system time and normalizes it (e.g., converting it to an annual day or a value between 0 and 1) to obtain the current normalized time. Then, it integrates the virtual station coordinates (longitude and latitude), the retrieved geographic elevation, and the calculated current normalized time to determine the virtual station's spatiotemporal parameters. Following the data arrangement order defined during model training, the system sequentially splices and stacks the satellite azimuth, satellite elevation, virtual station coordinates, geographic elevation, current normalized time, and the synthetic zenith wet delay value calculated in the previous step, combining them into a high-dimensional virtual physical feature vector. This vector fully describes the physical state under virtual observation conditions and can be directly input into the prediction model.

[0108] S207: Establish an enhanced tomographic equation set using real observations from Global Navigation Satellite System (GNSS) stations and virtual slant path wet delay observations.

[0109] For example, the system fuses sparsely distributed real observation data from the Global Navigation Satellite System with dense virtual observation data generated by machine learning predictions, thereby significantly expanding the effective observation sample size for tomographic inversion. By calculating the geometric intercepts of all signal propagation paths (including real signal paths and virtual ray paths) in the discretized three-dimensional grid model, a coefficient matrix describing the linear relationship between observations and grid medium parameters is constructed. This expands the initial set of equations, which originally exhibited ill-conditioned or rank-deficient characteristics due to insufficient ray coverage, into an enhanced linear set of equations with sufficient constraints and higher numerical stability.

[0110] In one possible implementation, an enhanced tomographic equation set is established using real observations and virtual slant path wet delay observations from Global Navigation Satellite System (GNSS) stations. Specifically, this includes: vertically stacking the real slant path wet delay vector corresponding to the real observations and the virtual slant path wet delay vector corresponding to the virtual slant path wet delay observations to generate an enhanced observation vector; calculating the real intercept length of the real signal path within the tomographic grid and the virtual intercept length of the virtual ray path within the tomographic grid based on the spatial geometric relationship between the real signal path and the virtual ray path corresponding to the real observations and a preset tomographic grid; constructing a real signal path design matrix based on the real intercept length and a virtual signal path design matrix based on the virtual intercept length; vertically stacking the real signal path design matrix and the virtual signal path design matrix to generate an enhanced design matrix, and establishing a linear observation equation in conjunction with the enhanced observation vector.

[0111] In this embodiment, the enhanced tomography equations refer to a linear mathematical model established by integrating real global navigation satellite system observation data with virtual observation data generated by machine learning. This model describes the physical summation between the total wet delay along the signal path and the water vapor density (or wet refractive index) within each grid the path passes through. For example, it is typically expressed in the form y=Ax, where y is an enhancement vector containing both real and virtual observations, A is an enhancement design matrix describing the geometric information of all rays passing through the grid, and x is the grid parameter to be inverted. By introducing virtual equations, this equation set effectively improves the ill-conditioned problems of equations caused by insufficient observed rays in traditional tomography.

[0112] Specifically, to address the instability of tomographic inversion results caused by relying solely on sparse real sites, the system first vectorizes and integrates the observation data. The real oblique path wet delay vector, composed of all real observations, is placed on top, while the virtual oblique path wet delay vector, composed of virtual oblique path wet delay observations obtained using the prediction model, is placed below. These two vectors are then vertically stacked to generate a higher-dimensional, more information-rich enhanced observation vector. A ray tracing algorithm is used to process the geometric projection relationship. Based on the real signal path corresponding to the real observations and the simulated virtual ray path, the trajectories of these paths traversing a pre-defined tomographic grid (i.e., a three-dimensional voxel set) are calculated. The distance of each line segment passing through each independent grid cell is precisely measured, thereby calculating the real intercept length of the real signal path within the tomographic grid and the virtual intercept length of the virtual ray path within the tomographic grid. In calculating the intercept length, this embodiment specifically employs the Siddon algorithm or a ray tracing algorithm based on planar intersection. The ray equation is assumed to be P(t) = P start +t·(P end -P start ), where t∈[0,1], Pstart P represents the coordinate vector of the starting point of the ray (signal propagation path) in three-dimensional space. end Let represent the coordinate vector of the endpoint of the ray (signal propagation path) in three-dimensional space. Calculate the intersection parameters t of the ray with all X, Y, and Z planes of the mesh; sort all intersection parameters to obtain the sequence t0, t1, ..., t m ; Calculate the intercept length L between adjacent intersection points voxel :L voxel =(t n+1 -t n )·||P end -P start ||, where L voxel t represents the intercept length of a ray within a specific grid cell (unit: meters); n+1 , t n This represents the parameter values ​​corresponding to the entry and exit points of the ray when it crosses the mesh boundary; ||P end -P start || represents the total physical length from the ray's starting point (station) to its ending point (satellite or tropopause).

[0113] Furthermore, these intercept lengths are used as coefficients to fill the corresponding rows and columns of the matrix. A true signal path design matrix describing the true observation geometry is constructed based on the true intercept lengths, and a virtual signal path design matrix describing the virtual observation geometry is constructed based on the virtual intercept lengths. To maintain strict consistency with the dimensional structure of the observation vectors, the true signal path design matrix and the virtual signal path design matrix are vertically stacked, combining to generate an enhancement design matrix with a significantly increased number of rows. Based on the fundamental principles of tomography, this enhancement design matrix (coefficient matrix) is combined with the aforementioned enhanced observation vectors (constant vectors) to establish a linear observation equation with more sufficient constraints (i.e., the enhanced tomographic equation system), laying the foundation for subsequent high-precision solutions.

[0114] S208: An algebraic reconstruction algorithm is used to solve the enhanced tomography equations and invert the three-dimensional water vapor density distribution field of the target monitoring area.

[0115] For example, the system selects meteorological reanalysis data as the initial value for iteration and uses an algebraic reconstruction algorithm to successively approximate the solution of the above-mentioned enhanced tomography equations. By iteratively calculating the difference between the forward modeling observations and the enhanced observations, and correcting the wet refractive index within the grid accordingly, the residuals of the equations converge to the optimal solution. A thermodynamic conversion relationship is introduced to transform the mathematically solved wet refractive index field into a physically meaningful three-dimensional water vapor density distribution field, thereby achieving a refined reconstruction of the vertical and horizontal structure of atmospheric water vapor.

[0116] In one possible implementation, an algebraic reconstruction algorithm is used to solve the enhanced tomography equations to invert and generate a three-dimensional water vapor density distribution field of the target monitoring area. Specifically, this includes: interpolating each grid node of the tomography grid based on meteorological reanalysis data to generate a three-dimensional wet refractive index field; arranging the discrete values ​​in the three-dimensional wet refractive index field into a solution vector, and calculating the product of the enhanced design matrix and the solution vector to obtain the forward modeling observation vector; calculating the residual vector between the forward modeling observation vector and the enhanced observation vector, and iteratively correcting the discrete values ​​in the solution vector according to the residual vector until the residual vector meets the preset convergence condition; determining the temperature parameters of each grid node based on meteorological reanalysis data, and combining the physical conversion formula between wet refractive index and water vapor density to convert the converged solution vector into a three-dimensional water vapor density distribution field.

[0117] In this embodiment, the algebraic reconstruction algorithm refers to a numerical solution method based on row iteration, mainly used to handle large and sparse linear equation systems in computed tomography. It minimizes the residuals of the equation system by repeatedly correcting the estimates of unknown variables, thereby reconstructing the internal structure of the object. For example, in this water vapor tomography scenario, the algorithm does not attempt... Figure 1 Instead of calculating the inverse matrix all at once, we start from an initial atmospheric model and continuously adjust the wet refractive index values ​​within the grid using the observed signal delay until the delay calculated by the model is highly consistent with the actual observed delay.

[0118] Specifically, to obtain high-precision inversion results and accelerate computation, the system introduces prior information to construct initial values ​​for iteration. Based on meteorological reanalysis data, spatial interpolation is performed on each grid node of the tomographic grid to estimate the initial wet refractive index at each node, thereby generating an initial three-dimensional wet refractive index field. According to the grid index order, the discrete values ​​in the three-dimensional wet refractive index field are vectorized and rearranged to form the solution vector required for algorithm iteration (i.e., the initial solution vector). The product of the enhancement design matrix and this solution vector is calculated to simulate the current observation values, obtaining the forward calculation observation vector. During the iteration process, each correction to the unknowns follows the following algebraic reconstruction algorithm formula:

[0119]

[0120] Where, x k+1 x represents the vector to be solved after the (k+1)th iteration (i.e., the updated mesh wet refractive index value); k y represents the vector to be solved in the k-th iteration; λ represents the relaxation factor, which typically ranges from (0, 2) and is used to control the convergence speed and stability; i This represents the i-th observation (SWD value) in the augmented observation vector; a ijThis represents the element in the i-th row and j-th column of the design matrix (i.e., the intercept length of the i-th ray in the j-th grid). This represents the forward modeling observation at the k-th iteration. This represents the transpose of the row vector (i.e., the column vector).

[0121] Furthermore, by comparing the simulated and actual values, the residual vector between the forward-modeled observation vector and the enhanced observation vector is calculated. Using the Kaczmarz method or its variant, this residual is projected back onto the grid along the ray path. The discrete values ​​in the solution vector are iteratively corrected based on the residual vector. This process is repeated until the residual vector meets the preset convergence condition (e.g., the sum of squared residuals is less than a threshold). Since the direct solution yields the wet refractive index, it is necessary to determine the Kelvin temperature and other air temperature parameters for each grid node based on meteorological reanalysis data. Then, using the physical conversion formula between wet refractive index and water vapor density (which describes the thermodynamic relationship between water vapor density, temperature, and wet refractive index), the wet refractive index values ​​in the converged solution vector are converted one by one into a meteorologically significant three-dimensional water vapor density distribution field. The physical conversion formula is as follows:

[0122] Where, ρ v This represents the water vapor density obtained from the inversion (unit: g / m³). 3 ); N wet R represents the wet refractive index of the mesh after iterative convergence. v The specific gas constant representing water vapor is taken as 461.5 J·kg⁻¹. -1 ·K -1 T represents the air temperature (in K) of the grid node, determined by meteorological reanalysis data; k2' and k3 are the same as the atmospheric refractive index constants in the aforementioned empirical formula for wet refractive index.

[0123] Figure 3 This is a schematic diagram of a machine learning-based water vapor chromatography system according to an embodiment of this application. This system can be implemented through software, hardware, or a combination of both, forming all or part of the overall system. Figure 3 As shown, the system includes:

[0124] The sample construction module 301 is configured to construct a historical training sample set containing historical signal geometric parameters, historical station spatiotemporal parameters and historical zenith wet delay values, and to obtain the true value label of the slant path wet delay corresponding to the historical training sample set.

[0125] The model training module 302 is configured to iteratively train a preset gradient boosting decision tree model using historical training sample sets and ground truth labels of oblique path wet delay, thereby generating an oblique path wet delay prediction model.

[0126] The virtual station establishment module 303 is configured to establish virtual station coordinates in blank areas that do not cover global navigation satellite system stations based on the spatial distribution data of global navigation satellite system stations in the target monitoring area;

[0127] The path generation module 304 is configured to generate a virtual ray path for virtual station coordinates according to a preset azimuth sequence and a preset elevation sequence, thereby obtaining the virtual signal geometric parameters.

[0128] The parameter calculation module 305 is configured to acquire digital elevation model data and meteorological reanalysis data, extract the virtual station spatiotemporal parameters of the virtual station coordinates based on the digital elevation model data, and calculate the synthetic zenith wet delay value of the virtual station coordinates based on the meteorological reanalysis data.

[0129] The delay prediction module 306 is configured to construct a virtual physical feature vector based on the virtual signal geometric parameters, the virtual station spatiotemporal parameters and the synthetic zenith wet delay value, and input the virtual physical feature vector into the oblique path wet delay prediction model to obtain the virtual oblique path wet delay observation value.

[0130] Equation building module 307 is configured to build an enhanced tomographic equation set using real observations from Global Navigation Satellite System (GNSS) stations and virtual slant path wet delay observations.

[0131] The tomographic inversion module 308 is configured to solve the enhanced tomographic equations using an algebraic reconstruction algorithm to invert and generate a three-dimensional water vapor density distribution field of the target monitoring area.

[0132] Based on the above embodiments, as an optional embodiment, the model training module 302 is specifically used for: dividing 24 hours into multiple consecutive preset time windows according to preset time intervals; filtering historical training sample sets and true values ​​of wet delay along the oblique path to the corresponding preset time windows according to the time information contained in the spatiotemporal parameters of historical stations, obtaining time-segmented training data corresponding to each preset time window; training multiple preset gradient boosting decision tree sub-models using the training data of each time segment, obtaining multiple trained gradient boosting decision tree sub-models corresponding to each preset time window; establishing a calling index relationship between each preset time window and the corresponding trained gradient boosting decision tree model; and establishing an oblique path wet delay prediction model using each trained gradient boosting decision tree model and each calling index relationship.

[0133] Based on the above embodiments, as an optional embodiment, the virtual station establishment module 303 is specifically used for: dividing the target monitoring area into grids according to a preset latitude and longitude interval to obtain multiple geographic grid units; using the spatial distribution data of the Global Navigation Satellite System (GNSS) stations within the target monitoring area to determine the inclusion relationship of each geographic grid unit, and filtering out all blank geographic grid units that do not contain GNSS stations from each geographic grid unit; obtaining the longitude and latitude of the geometric center point of each blank geographic grid unit, and establishing each set of longitude and latitude as multiple virtual station coordinates, with one set of longitude and latitude corresponding to one virtual station coordinate.

[0134] Based on the above embodiments, as an optional embodiment, the parameter calculation module 305 is specifically used to: use meteorological reanalysis data to perform four-dimensional interpolation on the longitude, latitude, and time of the virtual station coordinates to construct atmospheric parameter vertical profiles corresponding to multiple atmospheric layers. The atmospheric parameter vertical profiles include the layer temperature, layer vapor pressure, and layer potential of each atmospheric layer; use a preset empirical formula for wet refractive index to calculate the wet refractive index of each atmospheric layer based on the layer temperature and layer vapor pressure; calculate the difference between adjacent layer potentials in the atmospheric parameter vertical profiles and convert the difference into the geometric layer thickness of each atmospheric layer based on gravitational acceleration; and sum the products of each wet refractive index and the corresponding geometric layer thickness using the geographical elevation of the virtual station coordinates as the starting height for integration to obtain the synthetic zenith wet delay value.

[0135] Based on the above embodiments, as an optional embodiment, the delay prediction module 306 is specifically used to: determine the satellite azimuth and satellite elevation angle of the virtual ray path as virtual signal geometric parameters; retrieve the geographic elevation corresponding to the virtual station coordinates from the digital elevation model data, and obtain the current normalized time; determine the virtual station coordinates, geographic elevation, and current normalized time as virtual station spatiotemporal parameters; and combine the satellite azimuth, satellite elevation angle, virtual station coordinates, geographic elevation, current normalized time, and synthetic zenith wet delay value into a virtual physical feature vector.

[0136] Based on the above embodiments, as an optional embodiment, the equation establishment module 307 is specifically used to: vertically stack the real oblique path wet delay vector corresponding to the real observation value and the virtual oblique path wet delay vector corresponding to the virtual oblique path wet delay observation value to generate an enhanced observation vector; calculate the real intercept length of the real signal path in the tomographic grid and the virtual intercept length of the virtual ray path in the tomographic grid according to the spatial geometric relationship between the real signal path and the virtual ray path corresponding to the real observation value and the preset tomographic grid; construct a real signal path design matrix based on the real intercept length and a virtual signal path design matrix based on the virtual intercept length; vertically stack the real signal path design matrix and the virtual signal path design matrix to generate an enhanced design matrix, and establish a linear observation equation in combination with the enhanced observation vector.

[0137] Based on the above embodiments, as an optional embodiment, the tomographic inversion module 308 is specifically used for: interpolating each grid node of the tomographic grid based on meteorological reanalysis data to generate a three-dimensional wet refractive index field; arranging the discrete values ​​in the three-dimensional wet refractive index field into a solution vector, and calculating the product of the enhancement design matrix and the solution vector to obtain the forward modeling observation vector; calculating the residual vector between the forward modeling observation vector and the enhancement observation vector, and iteratively correcting the discrete values ​​in the solution vector according to the residual vector until the residual vector meets the preset convergence condition; determining the temperature parameters of each grid node based on meteorological reanalysis data, and combining the physical conversion formula between wet refractive index and water vapor density to convert the converged solution vector into a three-dimensional water vapor density distribution field.

[0138] It should be noted that the system provided in the above embodiments is only illustrated by the division of the above functional modules. In actual applications, the above functions can be assigned to different functional modules as needed, that is, the internal structure of the device can be divided into different functional modules to complete all or part of the functions described above. In addition, the system and method embodiments provided in the above embodiments belong to the same concept, and the specific implementation process can be found in the method embodiments, which will not be repeated here.

[0139] This embodiment also discloses an electronic device, as shown in the reference. Figure 4 The electronic device may include: at least one processor 401, at least one communication bus 402, user interface 403, network interface 404, and at least one memory 405.

[0140] The communication bus 402 is used to enable communication between these components.

[0141] The user interface 403 may include a display screen and a camera. Optionally, the user interface 403 may also include a standard wired interface and a wireless interface.

[0142] The network interface 404 may optionally include a standard wired interface or a wireless interface (such as a Wi-Fi interface).

[0143] The processor 401 may include one or more processing cores. The processor 401 connects to various parts of the server using various interfaces and lines, and performs various server functions and processes data by running or executing instructions, programs, code sets, or instruction sets stored in memory 405, and by calling data stored in memory 405. Optionally, the processor 401 may be implemented using at least one hardware form of Digital Signal Processing (DSP), Field-Programmable Gate Array (FPGA), or Programmable Logic Array (PLA). The processor 401 may integrate one or a combination of several of the following: Central Processing Unit (CPU), Graphics Processing Unit (GPU), and modem. The CPU primarily handles the operating system, user interface, and applications; the GPU is responsible for rendering and drawing the content required for display; and the modem handles wireless communication. It is understood that the modem may also be implemented as a separate chip without being integrated into the processor 401.

[0144] The memory 405 may include random access memory (RAM) or read-only memory. Optionally, the memory 405 may include a non-transitory computer-readable storage medium. The memory 405 may be used to store instructions, programs, code, code sets, or instruction sets. The memory 405 may include a program storage area and a data storage area, wherein the program storage area may store instructions for implementing an operating system, instructions for at least one function (such as touch function, sound playback function, image playback function, etc.), instructions for implementing the above-described method embodiments, etc.; the data storage area may store data involved in the above-described method embodiments, etc. Optionally, the memory 405 may also be at least one storage device located remotely from the aforementioned processor 401. Figure 4As shown, the memory 405, which serves as a computer storage medium, may include an operating system, a network communication module, a user interface module, and an application program for a water vapor chromatography method based on machine learning.

[0145] exist Figure 4 In the electronic device shown, the user interface 403 is mainly used to provide an input interface for the user and to obtain the user input data; while the processor 401 can be used to call an application program stored in the memory 405 for a water vapor chromatography method based on machine learning. When executed by one or more processors 401, the electronic device executes one or more methods as described in the above embodiments.

[0146] It should be noted that, for the sake of simplicity, the foregoing method embodiments are all described as a series of actions. However, those skilled in the art should understand that this application is not limited to the described order of actions, as some steps may be performed in other orders or simultaneously according to this application. Furthermore, those skilled in the art should also understand that the embodiments described in the specification are preferred embodiments, and the actions and modules involved are not necessarily essential to this application.

[0147] In the above embodiments, the descriptions of each embodiment have different focuses. For parts not described in detail in a certain embodiment, please refer to the relevant descriptions of other embodiments.

[0148] In the several embodiments provided in this application, it should be understood that the disclosed apparatus can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of units is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the shown or discussed mutual couplings or direct couplings or communication connections may be through some service interfaces; indirect couplings or communication connections between apparatuses or units may be electrical or other forms.

[0149] The units described as separate components may or may not be physically separate. 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 the units can be selected to achieve the purpose of this embodiment, depending on actual needs.

[0150] Furthermore, the functional units in the various embodiments of this application can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit. The integrated unit can be implemented in hardware or as a software functional unit.

[0151] If the integrated unit is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage device (CMD). Based on this understanding, the technical solution of this application, in essence, or the part that contributes to the prior art, or all or part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a memory 405 and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods of the various embodiments of this application. The aforementioned memory 405 includes various media capable of storing program code, such as a USB flash drive, external hard drive, magnetic disk, or optical disk.

[0152] The foregoing description is merely an exemplary embodiment of this disclosure and should not be construed as limiting the scope of this disclosure. Any equivalent changes and modifications made in accordance with the teachings of this disclosure shall still fall within the scope of this disclosure. Those skilled in the art will readily conceive of other embodiments of this disclosure upon considering the disclosure in this specification. This application is intended to cover any variations, uses, or adaptations of this disclosure that follow the general principles of this disclosure and include common knowledge or customary techniques in the art not described in this disclosure. The specification and embodiments are considered exemplary only, and the scope of this application is defined by the claims.

Claims

1. A water vapor chromatography method based on machine learning, characterized in that, The method includes: Construct a historical training sample set containing historical signal geometric parameters, historical station spatiotemporal parameters, and historical zenith wet delay values, and obtain the true value labels of the slant path wet delay corresponding to the historical training sample set; The preset gradient boosting decision tree model is iteratively trained using the historical training sample set and the ground truth labels of the oblique path wet delay to generate an oblique path wet delay prediction model. Based on the spatial distribution data of Global Navigation Satellite System (GNSS) stations within the target monitoring area, virtual station coordinates are established in blank areas that do not cover the GNSS stations. For the virtual station coordinates, a virtual ray path is generated according to a preset azimuth angle sequence and a preset elevation angle sequence to obtain the virtual signal geometric parameters; Acquire digital elevation model data and meteorological reanalysis data, extract the spatiotemporal parameters of the virtual station coordinates based on the digital elevation model data, and calculate the synthetic zenith wet delay value of the virtual station coordinates based on the meteorological reanalysis data. A virtual physical feature vector is constructed based on the virtual signal geometric parameters, the virtual station spatiotemporal parameters, and the synthetic zenith wet delay value. The virtual physical feature vector is then input into the slant path wet delay prediction model to obtain the virtual slant path wet delay observation value. An enhanced tomographic equation set was established using the actual observations from the Global Navigation Satellite System (GNSS) stations and the virtual oblique path wet delay observations. The enhanced tomography equations are solved using an algebraic reconstruction algorithm to generate a three-dimensional water vapor density distribution field in the target monitoring area.

2. The method according to claim 1, characterized in that, The step of iteratively training a preset gradient boosting decision tree model using the historical training sample set and the ground truth labels of the oblique path wet delay to generate an oblique path wet delay prediction model specifically includes: The 24 hours are divided into multiple consecutive preset time windows according to preset time intervals; Based on the time information contained in the spatiotemporal parameters of the historical stations, the historical training sample set and the true labels of the oblique path wet delay are filtered to the corresponding preset time windows to obtain the time-segmented training data corresponding to each preset time window. Multiple preset gradient boosting decision tree sub-models are trained using the training data from each of the aforementioned time periods to obtain multiple trained gradient boosting decision tree models corresponding to each of the aforementioned preset time windows. Establish a calling index relationship between each preset time window and the corresponding trained gradient boosting decision tree sub-model; The oblique path wet delay prediction model is established by utilizing the trained gradient boosting decision tree sub-models and the calling index relationships.

3. The method according to claim 1, characterized in that, The step of establishing virtual station coordinates in blank areas not covered by the Global Navigation Satellite System (GNSS) stations within the target monitoring area, based on the spatial distribution data of GNSS stations, specifically includes: The target monitoring area is divided into grids according to preset latitude and longitude intervals to obtain multiple geographic grid units; Using the spatial distribution data of the Global Navigation Satellite System (GNSS) stations within the target monitoring area, the inclusion relationship of each geographic grid cell is determined, and all blank geographic grid cells that do not contain the GNSS stations are selected from each geographic grid cell. Obtain the longitude and latitude of the geometric center point of each blank geographic grid cell, and establish each set of longitude and latitude as multiple virtual station coordinates, with each set of longitude and latitude corresponding to one virtual station coordinate.

4. The method according to claim 1, characterized in that, The calculation of the synthetic zenith wet delay value of the virtual station coordinates based on the meteorological reanalysis data specifically includes: Using the meteorological reanalysis data, four-dimensional interpolation is performed on the longitude, latitude and time of the virtual station coordinates to construct vertical profiles of atmospheric parameters corresponding to multiple atmospheric layers. The vertical profiles of atmospheric parameters include the layer temperature, layer water vapor pressure and layer potential of each atmospheric layer. Using a preset empirical formula for wet refractive index, the wet refractive index of each atmospheric layer is calculated based on the stratification temperature and the stratification vapor pressure. Calculate the difference between the potential of adjacent layers in the vertical profile of the atmospheric parameters, and convert the difference into the geometric thickness of each atmospheric layer based on the gravitational acceleration. Using the geographical elevation of the virtual site coordinates as the starting height for integration, the products of each wet refractive index and the corresponding geometric layer thickness are summed to obtain the synthetic zenith wet delay value.

5. The method according to claim 4, characterized in that, The construction of the virtual physical feature vector based on the virtual signal geometric parameters, the virtual station spatiotemporal parameters, and the synthetic zenith wet delay value specifically includes: The satellite azimuth and satellite elevation angle of the virtual ray path are determined as the geometric parameters of the virtual signal; The geographic elevation corresponding to the virtual station coordinates is retrieved from the digital elevation model data, and the current normalized time is obtained. The virtual station coordinates, the geographic elevation, and the current normalized time are determined as the spatiotemporal parameters of the virtual station. The satellite azimuth angle, the satellite elevation angle, the virtual station coordinates, the geographic elevation, the current normalized time, and the synthetic zenith wet delay value are combined to form the virtual physical feature vector.

6. The method according to claim 1, characterized in that, The process of establishing an enhanced tomographic equation set using the actual observations from the Global Navigation Satellite System (GNSS) stations and the virtual slant path wet delay observations specifically includes: The real slant path wet delay vector corresponding to the real observation value and the virtual slant path wet delay vector corresponding to the virtual slant path wet delay observation value are vertically stacked to generate an enhanced observation vector. Based on the actual signal path corresponding to the actual observation value and the spatial geometric relationship between the virtual ray path and the preset tomographic grid, calculate the actual intercept length of the actual signal path within the tomographic grid and the virtual intercept length of the virtual ray path within the tomographic grid; A real signal path design matrix is ​​constructed based on the real intercept length, and a virtual signal path design matrix is ​​constructed based on the virtual intercept length. The real signal path design matrix and the virtual signal path design matrix are vertically stacked to generate an enhanced design matrix, and a linear observation equation is established by combining the enhanced observation vector.

7. The method according to claim 6, characterized in that, The process of solving the enhanced tomographic equations using an algebraic reconstruction algorithm to invert and generate the three-dimensional water vapor density distribution field of the target monitoring area specifically includes: Based on the meteorological reanalysis data, interpolation is performed on each grid node of the tomographic grid to generate a three-dimensional wet refractive index field; The discrete values ​​in the three-dimensional wet refractive index field are arranged into a vector to be solved, and the product of the enhancement design matrix and the vector to be solved is calculated to obtain the forward modeling observation vector. Calculate the residual vector between the forward modeling observation vector and the enhanced observation vector, and iteratively correct the discrete values ​​in the vector to be solved based on the residual vector until the residual vector satisfies the preset convergence condition; Based on the meteorological reanalysis data, the temperature parameters of each grid node are determined, and combined with the physical conversion formula between wet refractive index and water vapor density, the converged solution vector is converted into the three-dimensional water vapor density distribution field.

8. A water vapor chromatography system based on machine learning, characterized in that, The system includes: The sample construction module is configured to construct a historical training sample set containing historical signal geometric parameters, historical station spatiotemporal parameters, and historical zenith wet delay values, and to obtain the true value label of the slant path wet delay corresponding to the historical training sample set. The model training module is configured to iteratively train a preset gradient boosting decision tree model using the historical training sample set and the ground truth labels of the oblique path wet delay, thereby generating an oblique path wet delay prediction model. The virtual station establishment module is configured to establish virtual station coordinates in blank areas that do not cover the Global Navigation Satellite System (GNSS) stations, based on the spatial distribution data of GNSS stations within the target monitoring area. The path generation module is configured to generate a virtual ray path for the virtual station coordinates according to a preset azimuth sequence and a preset elevation sequence, thereby obtaining virtual signal geometric parameters. The parameter calculation module is configured to acquire digital elevation model data and meteorological reanalysis data, extract the virtual station spatiotemporal parameters of the virtual station coordinates based on the digital elevation model data, and calculate the synthetic zenith wet delay value of the virtual station coordinates based on the meteorological reanalysis data. The delay prediction module is configured to construct a virtual physical feature vector based on the geometric parameters of the virtual signal, the spatiotemporal parameters of the virtual station, and the synthetic zenith wet delay value, and input the virtual physical feature vector into the oblique path wet delay prediction model to obtain the virtual oblique path wet delay observation value. The equation establishment module is configured to establish an enhanced tomographic equation set using the real observations from the Global Navigation Satellite System (GNSS) stations and the virtual oblique path wet delay observations. The tomographic inversion module is configured to solve the enhanced tomographic equations using an algebraic reconstruction algorithm to invert and generate a three-dimensional water vapor density distribution field of the target monitoring area.

9. An electronic device, characterized in that, The device includes a processor, a memory, a user interface, and a network interface. The memory is used to store instructions. The user interface and the network interface are both used to communicate with other devices. The processor is used to execute the instructions stored in the memory to cause the electronic device to perform the method as described in any one of claims 1-7.

10. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores instructions that, when executed, perform the method as described in any one of claims 1-7.

Images