A precipitation spatial extrapolation method based on edge-enhanced graph attention network

Abstract

This invention discloses a spatial inference method for precipitable water content based on an edge-enhanced graph attention network, belonging to the field of atmospheric remote sensing and spatial information processing technology. It includes: constructing a graph structure based on the geographic coordinates of GNSS stations, with nodes corresponding to stations and edges representing spatial relationships; generating node features based on the geographic coordinates of the stations, meteorological reanalysis data, topographic elevation, time coding, and observed precipitable water content; generating edge features based on the geographic distance, elevation difference, longitude difference, latitude difference, and normalized distance weights between stations; and inputting the node and edge features into an edge-enhanced graph attention network to aggregate neighbor information and predict precipitable water content through an attention mechanism. This invention also employs a masked self-supervised training strategy and supports adding virtual grid nodes to the graph structure to reconstruct a two-dimensional precipitable water content field. This invention can explicitly model anisotropic spatial dependencies in irregular station networks, achieving high-precision spatial inference of precipitable water content and continuous field reconstruction.

Description

Technical Field

[0001] This invention belongs to the field of atmospheric remote sensing and spatial information processing technology, and particularly relates to a spatial inference method for precipitable water based on edge-enhanced graph attention network. Background Technology

[0002] Precipitable water is a fundamental parameter describing atmospheric radiation transfer, cloud formation, and precipitation processes. Its accurate spatiotemporal characterization plays a crucial supporting role in numerical weather prediction, climate monitoring, and severe convective weather early warning. Among ground-based atmospheric observation methods, the Global Navigation Satellite System (GNSS) network can provide continuous precipitable water estimates around the clock and in all weather conditions, with a temporal resolution ranging from minutes to hours.

[0003] However, GNSS stations exhibit uneven and irregular spatial distribution: they are densely distributed in urban and tectonically active areas, while coverage is sparse in mountainous, arid, and marine regions. This spatial heterogeneity severely restricts the direct use of GNSS precipitable water data in applications requiring gridded products, such as data assimilation and remote sensing verification.

[0004] Existing techniques for extending station-based precipitable water coverage observations into a spatially continuous field mainly follow two routes. The first route employs deterministic interpolation methods, such as inverse distance weighting, ordinary kriging, and spline fitting, to generate grid products from station observations. While these methods are intuitive and theoretically sound, they rely on assumptions of spatial stationarity and isotropy, which are frequently broken in topographically complex regions, where the precipitable water gradient is strongly modulated by topographic and synoptic-scale forcing. The second route fuses GNSS observations with auxiliary information, employing data fusion or regression techniques. While fusion methods can improve spatial coverage, they often inherit systematic biases from auxiliary data, potentially masking the independent information provided by GNSS observations.

[0005] In recent years, machine learning methods have been used for tasks related to precipitable water availability, including station-level inversion, temporal prediction, and spatial modeling. However, fully connected networks and convolutional architectures are designed for regular grids or fixed-length inputs and cannot natively adapt to the variable topology of real GNSS network. Different stations may have different numbers of neighboring stations, and the propagation of information between stations depends not only on distance but also on elevation difference, latitude, and similarity of meteorological conditions—relationships that Euclidean grid methods cannot express.

[0006] Therefore, there is an urgent need for a spatial inference method for precipitable water that can adapt to irregular GNSS station network topology and explicitly model anisotropic spatial dependencies. Summary of the Invention

[0007] To address the aforementioned technical problems, this invention provides a spatial inference method for precipitable water based on edge-enhanced graph attention networks, comprising the following steps: A graph structure is constructed based on the geographic coordinates of multiple Global Navigation Satellite System (GNSS) stations. The graph structure includes nodes and edges, where the nodes correspond to GNSS stations and the edges represent the spatial relationships between the stations. Based on the geographic coordinates, meteorological reanalysis data, topographic elevation, periodic time coding, and precipitable water observations of each global navigation satellite system station, node characteristics of each node are generated. Based on the geographical distance, elevation difference, longitude difference, latitude difference, and normalized distance weight between the two stations, the edge features of the edge connecting the two nodes are generated. The node features and edge features are input into the edge-enhanced graph attention network. The precipitation information of neighboring nodes is aggregated through the attention mechanism to predict the precipitation of the target node. Based on the prediction results of each target node, a spatially continuous precipitation field is generated.

[0008] Optionally, constructing the graph structure includes: Based on the geographical coordinates of each global navigation satellite system station, the set of neighboring nodes for each node is determined using a combination of K-nearest neighbor and distance threshold criteria, and edges connecting the node and its neighboring nodes are established.

[0009] Optionally, the generated node features include: Based on the geographic coordinates of each global navigation satellite system station, surface temperature and pressure, topographic elevation, periodic time coding, and precipitable water observation values ​​from meteorological reanalysis data, a multidimensional node feature vector is constructed.

[0010] Optionally, the generated edge features include: A multidimensional edge feature vector is constructed based on the geographical distance, elevation difference, longitude difference, latitude difference, and normalized distance weight between the two stations.

[0011] Optionally, it also includes a mask self-supervised training step: In each training iteration, a portion of the precipitation observations from the Global Navigation Satellite System stations are randomly masked. The precipitation features of the masked stations are set to zero, while the features of other nodes are retained. The training objective is to predict the precipitation of the masked stations.

[0012] Optionally, it also includes: During the training process, the proportion of randomly occluded sites varies within a preset range.

[0013] Optionally, it also includes a two-dimensional precipitation field reconstruction step: In the study area, regularly distributed virtual grid nodes are generated. Each virtual grid node is assigned geographic coordinates, topographic elevation, and surface air temperature and surface air pressure from meteorological reanalysis data. The precipitation characteristics of the virtual grid nodes are set to zero. Virtual grid nodes are added to the graph structure as new nodes and connected to nearby global navigation satellite system station nodes; The trained edge-enhanced graph attention network is applied to the expanded graph structure to predict the precipitation of all virtual grid nodes, generating a spatially continuous precipitation field.

[0014] Optionally, adding virtual grid nodes to the graph structure includes: Based on the geographic coordinates of the virtual grid node and the preset K-nearest neighbor and distance threshold combination criteria, the set of neighbor nodes of the virtual grid node is determined, and the edges connecting the virtual grid node and its neighbor nodes are established.

[0015] Optionally, the edge-enhanced graph attention network includes several stacked graph attention layers, each layer using a multi-head attention mechanism, and calculates attention coefficients by concatenating edge features with the representations of the source and target nodes.

[0016] Optionally, the loss function of the edge-enhanced graph attention network includes a reconstruction loss for occluded nodes and a spatial smoothing regularization term based on edge weights, wherein the edge weights are determined based on the geographical distance between nodes.

[0017] Compared with the prior art, the present invention has the following advantages and technical effects: First, by explicitly modeling the GNSS station network as a graph structure and using an edge-enhanced attention mechanism to aggregate neighbor information, it can adapt to arbitrary topological structures and heterogeneous station densities of irregular station networks, overcome the dependence of traditional interpolation methods on spatial stability and isotropy, and significantly improve the accuracy of precipitable water inference in complex terrain areas.

[0018] Second, by integrating multidimensional side features such as geographical distance, elevation difference, longitude difference, latitude difference, and normalized distance weights, the model can learn anisotropic spatial dependencies dominated by topography and atmospheric gradients, effectively capturing the changes in water vapor gradients from coastal to inland areas and from low to high altitudes.

[0019] Third, a masked self-supervised training strategy is adopted to randomly mask the observations of stations with variable ratios, enabling the model to learn stable spatial interpolation capabilities under different missing data configurations and directly support inductive generalization to unseen locations (including virtual grid nodes). End-to-end reconstruction of two-dimensional precipitable water fields at arbitrary resolution can be achieved without retraining.

[0020] Fourth, through verification by independent radiosonde stations, the precipitable water inferred by this invention maintains a high degree of consistency with the actual radiosonde measurements, proving that the graph structure spatial modeling method can generate a physically meaningful and spatially coherent precipitable water field from sparse and irregularly distributed GNSS observations. Attached Figure Description

[0021] The accompanying drawings, which form part of this application, are used to provide a further understanding of this application. The illustrative embodiments and descriptions of this application are used to explain this application and do not constitute an undue limitation of this application. In the drawings: Figure 1 This is a schematic diagram of the method flow according to an embodiment of the present invention; Figure 2 This is a schematic diagram of the California GNSS station network structure according to an embodiment of the present invention; Figure 3 This is a convergence curve of the training and validation loss of the edge-enhanced graph attention network according to an embodiment of the present invention. Figure 4 This is a scatter plot of predicted and observed values ​​from the edge-enhanced graph attention network in an embodiment of the present invention. Figure 5 This is a spatial distribution diagram of the root mean square error of the California GNSS station network according to an embodiment of the present invention. Figure 6 This is a scatter plot of the inferred precipitable water volume and the measured precipitable water volume by radiosonde, according to an embodiment of the present invention. Figure 7 This is a bar chart showing the station-by-station verification error of the six radiosonde stations in this embodiment of the invention; Figure 8 This is a comparison chart of the time series of precipitable water inferred by the edge-enhanced graph attention network and the precipitable water obtained by radiosonde, according to an embodiment of the present invention. Figure 9 This is a comparison diagram of the two-dimensional precipitable water field reconstructed by the edge-enhanced graph attention network in an embodiment of the present invention and the water vapor in the ERA5 total column. Detailed Implementation

[0022] It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0023] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0024] Example 1 like Figure 1 As shown, this embodiment provides a spatial inference method for precipitable water based on edge-enhanced graph attention networks, including the following steps: A graph structure is constructed based on the geographic coordinates of multiple Global Navigation Satellite System (GNSS) stations. The graph structure includes nodes and edges. Nodes correspond to GNSS stations, and edges represent the spatial relationships between stations. Based on the geographic coordinates, meteorological reanalysis data, topographic elevation, periodic time coding, and precipitable water observations of each global navigation satellite system station, node characteristics of each node are generated. Based on the geographical distance, elevation difference, longitude difference, latitude difference, and normalized distance weight between the two stations, the edge features of the edge connecting the two nodes are generated. The node features and edge features are input into the edge-enhanced graph attention network. The precipitation information of neighboring nodes is aggregated through the attention mechanism to predict the precipitation of the target node. Based on the prediction results of each target node, a spatially continuous precipitation field is generated.

[0025] In this embodiment, California is used as an example. California is an ideal test area for PWV spatial inference due to its diverse topographical conditions (from coastal plains at sea level to the Sierra Nevada Mountains with elevations exceeding 4,000 meters) and the significant climatic gradient between the humid Pacific coast, arid inland valleys, and mountainous environments. The distribution of GNSS stations inherits this geographical complexity, with dense station coverage in urban corridors (Los Angeles, San Francisco Bay Area) and sparse coverage in the eastern desert and mountainous inland areas.

[0026] The hourly PWV estimates were obtained from the Nevada Geodetic Laboratory (NGL), which uses precise point positioning to routinely process global GNSS observation data, providing tropospheric products including total zenith delay (TWD) and integrated PWV. This embodiment implemented a four-level quality control process: (1) physical range filtering to remove values ​​below 0 or above 100 mm; (2) sliding window 3σ filtering to remove transient spikes; (3) cross-validation with ERA5 Total Column Water Vapor (TCWV) to identify systematic anomalies; and (4) retaining only the correlation coefficient R with ERA5. 2 Stations with a score >0.90. After quality control, 193 stations were retained, with an average data integrity rate exceeding 98%.

[0027] This embodiment uses three types of auxiliary data to construct the node and edge features of the graph: (1) ERA5 reanalysis data. Hourly 2m air temperature (T) for each station location was extracted from the ERA5 global reanalysis. 2m ) and surface air pressure (P s The field (bilinear interpolation) uses these meteorological variables as surrogate features to capture the synoptic-scale atmospheric conditions that control PWV variations. The ERA5 TCWV is also used as an independent reference for two-dimensional field assessments.

[0028] (2) Digital Elevation Model (DEM). The station orthometric height is obtained from metadata and converted to geodetic height using the EGM2008 geoid model to ensure consistency with the GNSS coordinate system. The virtual grid node elevations used for 2D reconstruction are taken from the ERA5 topographic field (approximately 9 km resolution) or the Copernicus GLO-30 DEM (30 m resolution).

[0029] (3) Radiosonde observations. Six radiosonde stations within or near the California network (World Meteorological Organization (WMO) numbers: 72293, 72388, 72489, 72493, 72582, and 72393) provided independent PWV estimates derived from the vertical integral of the humidity profile. These observations were used for independent validation only and were not used for model training.

[0030] Specifically, the dataset is divided by time to prevent data leakage: observations from 2019–2020 constitute the training set, and data from 2021 constitute the validation set. This time division ensures that the model is evaluated on future observations not seen during training, reflecting actual business deployment conditions.

[0031] Furthermore, the construction of the graph structure includes: determining the set of neighboring nodes for each node based on the geographical coordinates of each global navigation satellite system station, using a combination of K-nearest neighbor and distance threshold criteria, and establishing edges connecting the node with its neighboring nodes.

[0032] The generated node features include: constructing a multi-dimensional node feature vector based on the geographic coordinates of each global navigation satellite system station, surface air temperature and surface air pressure, topographic elevation, periodic time coding, and precipitable water observation values ​​from meteorological reanalysis data.

[0033] The generation of edge features includes: constructing a multidimensional edge feature vector based on the geographical distance, elevation difference, longitude difference, latitude difference, and normalized distance weights between two stations.

[0034] Specifically, the GNSS station network is represented as a directed graph G=(V,E), where each node v i ∈V corresponds to a GNSS station, and each edge e ij ∈E represents station v i and v jThe spatial relationships between them are determined. Edges are constructed using a combination of k-Nearest Neighbors (KNN) and distance thresholds: each station is connected to its k nearest neighbors (k=8), and additional edges are added for any station within a 200km distance threshold. This hybrid strategy strikes a balance between local connectivity in dense urban areas and long-distance connectivity in sparse areas. Figure 2 A graph structure was constructed on the California GNSS network based on a combination of KNN (k=8) and distance threshold (200km) criteria. The transparency of the edges reflects the distance between stations, with denser connections in coastal and urban areas.

[0035] Each node carries an 11-dimensional feature vector x. i ∈R 11 This includes PWV values ​​(observations from unshaded stations; values ​​from shaded stations are set to zero), geographic coordinates (longitude, latitude), geodetic height, and ERA5 meteorological variables (T). 2m P s This includes periodic time encoding (sine and cosine values ​​for months and hours). Periodic encoding preserves the periodicity of daily variations and seasonal signals while avoiding discontinuities at boundaries.

[0036] Each edge carries a five-dimensional feature vector e. ij ∈R 5 , encoding geographical distance d ij Elevation difference ∆h ij Longitude difference ∆λ ij Latitude difference ∆φ ij and normalized distance weights ( (This represents the average station spacing). These edge features explicitly encode the anisotropic spatial structure of the station network: elevation differences capture the strong correlation between PWV and altitude, directional components (longitude and latitude differences) encode the water vapor gradient from the coast to the inland and zonal, and distance weights provide soft distance attenuation.

[0037] Furthermore, the edge-enhanced graph attention network includes several stacked graph attention layers. Each layer uses a multi-head attention mechanism and calculates attention coefficients by concatenating edge features with the representations of the source and target nodes. In this embodiment, it consists of three sequential components: a feature embedding layer, a stacked edge-enhanced GAT module, and an output prediction head.

[0038] (1) Feature embedding. The 11-dimensional original node features are transformed by linear transformation. The image is projected into a 64-dimensional latent space and then injected with high-frequency spatial information through Fourier position encoding.

[0039] (2) Edge-enhancing GAT module. The model contains three stacked GAT modules. In each module, multi-head attention is computed using four parallel attention heads, and the hidden layer dimension is 64. Nodes in the ℓth layer i and j The attention coefficient between them is: Among them W and W e Let be the learnable weight matrix, 'a' be the attention vector, '∥' denote concatenation, and N( i ) is a node i The neighbor set. Edge feature vector e ij Attention scores are calculated after concatenating the source and target node representations, enabling the network to adjust information flow based on inter-site geographical relationships. Each GAT module contains residual skip connections and layer normalization to stabilize training.

[0040] Node representation is updated by aggregating neighbor messages: Where H=4 represents the number of attention heads, and the average between heads is used for aggregation. Dropout (p=0.1) is applied to the attention coefficients during training.

[0041] (3) Prediction Head. Final Node Representation The predicted PWV value of each node is output through a two-layer MLP (64→32→1, ELU activation). .

[0042] Furthermore, this method also includes a mask self-supervised training step: In each training iteration, a portion of the precipitable water observations from Global Navigation Satellite System (GNSS) stations are randomly obscured. The precipitable water features of the obscured stations are set to zero, while the features of other nodes are retained. The training objective is to predict the precipitable water of the obscured stations. During training, the proportion of randomly obscured stations varies within a preset range.

[0043] The loss function of the edge-enhanced graph attention network includes the reconstruction loss of occluded nodes and a spatial smoothing regularization term based on edge weights, wherein the edge weights are determined based on the geographical distance between nodes.

[0044] Specifically, PWVGAT employs a mask self-supervised training strategy inspired by a mask autoencoder. In each training iteration, m stations are randomly selected as the mask set M, where m are uniformly sampled to ensure the mask ratio |M| / |V| is within the range [0.2, 0.5]. For masked stations, the PWV features in the input vector are set to zero, while the remaining 10 auxiliary features are retained. The model is trained to predict the PWV values ​​of masked stations under conditions of complete observations of unmasked stations and auxiliary features of all stations.

[0045] The training objective combines the reconstruction loss of occluded sites with spatial smoothing regularization: in and These are the actual and predicted PWV values, respectively. For distance-based edge weights, =0.01 controls the strength of the smoothing penalty. The smoothing term encourages spatially coherent predictions by penalizing the prediction differences between neighboring nodes (weighted by distance).

[0046] The model uses the AdamW optimizer (learning rate 0.001, weight decay 10). −4 Optimization is performed using cosine annealing scheduling (maximum 200 epochs). An early stopping strategy based on validation loss (tolerance of 20 epochs) prevents overfitting. All input features are standardized to zero mean and unit variance using training set statistics. The random seed is fixed at 42 to ensure repeatability.

[0047] The variable masking ratio during training serves a dual purpose: first, to prevent the model from relying on a fixed neighborhood density, forcing it to learn stable spatial interpolation under different missing data configurations; second, to directly support inductive inference, i.e., at deployment time, the model can predict the PWV of completely new locations (virtual grid nodes), treating these locations as masked nodes in the extended graph.

[0048] Furthermore, this method also includes a two-dimensional precipitation field reconstruction step: In the study area, regularly distributed virtual grid nodes are generated. Each virtual grid node is assigned geographic coordinates, topographic elevation, and surface air temperature and surface air pressure from meteorological reanalysis data. The precipitation characteristics of the virtual grid nodes are set to zero. Virtual grid nodes are added to the graph structure as new nodes and connected to nearby global navigation satellite system station nodes; The trained edge-enhanced graph attention network is applied to the expanded graph structure to predict the precipitation of all virtual grid nodes, generating a spatially continuous precipitation field.

[0049] Including the addition of virtual grid nodes to the graph structure includes: Based on the geographic coordinates of the virtual grid node and the preset K-nearest neighbor and distance threshold combination criteria, the set of neighbor nodes of the virtual grid node is determined, and the edges connecting the virtual grid node and its neighbor nodes are established.

[0050] Specifically, to reconstruct a continuous PWV field, a regular virtual node grid with the required spatial resolution was used to cover the study area (0.25° in this study, generating 1763 grid nodes). Each virtual node was assigned geographic coordinates and elevation from an external DEM, as well as meteorological features (T). 2m P s Extracted from ERA5 at grid coordinates. Virtual nodes are connected to the graph using the same KNN and distance threshold criteria as station nodes, constructing an extended graph. Since virtual nodes have no PWV observations, their PWV input features are set to zero, similar to occluded nodes during training. The trained PWVGAT model performs a forward propagation on the extended graph, simultaneously predicting the PWV of all virtual grid nodes, generating a two-dimensional PWV field. Table 1 summarizes the complete model configuration and training hyperparameters.

[0051] Table 1 The following verifications were also performed in this embodiment: Training convergence: Figure 3 The training and validation loss convergence curves for PWVGAT are shown. The vertical dashed line marks the optimal epoch (epoch 67), after which early stopping is triggered. The model converges after approximately 80 epochs, with the optimal validation loss occurring in epoch 67. The training and validation loss curves show a stable decrease with no significant separation, indicating good generalization and no overfitting. The final validation metrics (normalized feature space) are: RMSE = 0.2194, MAE = 0.1550, bias = −0.0017, R0 2 =0.9736.

[0052] This embodiment also underwent an accuracy evaluation compared to a comparative method: PWVGAT was evaluated against two comparative methods using the same validation set and the same masked site configuration: (1) Inverse distance weighted interpolation (IDW): a classic deterministic interpolation method that estimates the value of the target location as a distance-weighted average of the neighboring observation stations.

[0053] (2) (MLP): Using the same 10-dimensional auxiliary features as PWVGAT but without graph structure, it represents a general machine learning method without spatial adjacency information.

[0054] Table 2 presents the comparison results. The RMSE of PWVGAT is significantly lower than that of the two comparison methods, decreasing by 46.6% compared to IDW and 78.3% compared to MLP. R²=0.9736 indicates that PWVGAT explains more than 97% of the PWV variance, while IDW explains 90.8% and MLP only 43.7%. The poor performance of MLP demonstrates that auxiliary features alone, without explicit spatial information from neighboring sites, are insufficient for accurate PWV inference; spatial structure is a necessary condition.

[0055] Table 2 All indicators in the table are calculated in the normalized feature space. Optimal values ​​are bolded, and arrows indicate directions for improvement (↓ lower is better, ↑ higher is better).

[0056] Figure 4 The scatter plot shows the PWVGAT predicted values ​​and observed values. The data points are closely clustered around the 1:1 line in the entire PWV range, with no systematic bias at the low and high value ends. Figure 5 The spatial distribution of RMSE at each station is shown. Overall, the performance is uniform within the station network, with only a slight increase in RMSE at a few stations in areas with complex terrain.

[0057] To evaluate PWVGAT performance in physical units and compare it with completely independent reference data, the model predictions were compared with the integral PWV of humidity profiles from six radiosonde stations in or near the study area. Figure 6 The sounding data was not used for model training or hyperparameter tuning, providing truly independent validation.

[0058] Table 3 lists the station-by-station and overall validation statistics. In 3933 time-matched samples, the overall RMSE of PWVGAT was 2.12 mm, the MAE was 1.62 mm, the bias was −1.40 mm, and the correlation coefficient was 0.979. The slight negative bias indicates a slight tendency for PWVGAT to underestimate PWV relative to radiosonde, possibly reflecting a systematic difference between the GNSS and radiosonde atmospheric column integration methods.

[0059] Table 3 Where n is the number of time-matched samples.

[0060] The performance of each station exhibits a physically interpretable spatial gradient. Figure 7The lowest RMSE (1.35 mm, r = 0.993) was observed at station 72293, located in a mild coastal environment with dense surrounding GNSS coverage, while the highest RMSE (3.49 mm, r = 0.939) was observed at station 72393, located in an inland desert with sparse surrounding station density. This gradient reflects that PWVGAT performance is controlled by the local density and representativeness of surrounding GNSS observations.

[0061] Figure 8 The time series comparison of station 72293 throughout 2021 is presented. PWVGAT closely tracks the diurnal variation and synoptic-scale fluctuations of radiosonde PWV across the full dynamic range, accurately capturing summer moisture peaks and winter dry periods.

[0062] The two-dimensional PWV field reconstruction is as follows: The inductive inference capability of PWVGAT was validated by reconstructing a two-dimensional PWV field covering California using a 0.25° regular grid. The reconstructed field, generated by 184 active sites and 1763 virtual grid nodes on July 15, 2021, at 12:00 UTC, exhibited the expected spatial pattern: lower PWV along the Pacific coast, higher PWV in the warmer Central Valley, and lower PWV in the high-altitude regions of the Sierra Nevada and Cascade Mountains.

[0063] Figure 9 The three-panel comparison between the reconstruction field and the ERA5 TCWV is shown. Figure 9 (a) in the figure represents the PWV field inferred by PWVGAT. Figure 9 (b) in the image represents the ERA5 TCWV reference field. Figure 9 In the diagram, (c) represents the difference between the two fields (PWVGAT−ERA5). The two fields exhibit a highly consistent spatial pattern. The difference field reveals that PWVGAT is slightly underestimated relative to ERA5 in some high-altitude areas, while producing comparable or slightly higher values ​​in coastal areas. Table 4 summarizes the comparative results of the three virtual node elevation assignment schemes. Using ERA5 terrain as the elevation source yields the best consistency (RMSE=2.89mm, correlation coefficient 0.970), followed by Copernicus GLO-30 DEM (RMSE=3.03mm, correlation coefficient 0.967) and IDW interpolated station elevations (RMSE=3.54mm, correlation coefficient 0.957). The best performance of ERA5 terrain reflects self-consistency because the reference field (ERA5 TCWV) shares the same atmospheric background field with ERA5 terrain. Although the RMSE of GLO-30 is slightly larger, it represents a more realistic description of the terrain and may be superior in applications requiring independence from reanalysis assumptions.

[0064] Table 4 Its assessment is limited to land grids ( n =1268).

[0065] Experimental results show that explicit graph structure modeling has a clear advantage over classical interpolation and general machine learning methods in PWV spatial inference. The 46.6% reduction in RMSE compared to IDW is particularly significant, as IDW, although simple, is widely used in geodesy for GNSS PWV gridding and is considered a competitive comparative method. The improvement stems from two aspects. First, the attention mechanism learns data-driven, station-specific aggregation weights, rather than relying on a pre-defined distance decay function. This flexibility allows the model to adapt to the heterogeneous station geometries of real networks. Second, multi-dimensional edge features encode physical relationships (elevation difference, direction components) completely ignored by IDW, making it possible for the model to handle the strong anisotropy of the PWV field in mountainous terrain. The low performance of MLP (R²=0.437) reveals another important insight: auxiliary features alone (meteorological variables, coordinates, time encoding) are insufficient for accurate PWV inference. In the absence of PWV observations from neighboring stations, MLP is essentially a regression model that maps environmental conditions to PWV, equivalent to a simplified version of the empirical PWV model. The graph structure provides crucial missing information, namely the actual PWV status of neighboring sites, validating the hypothesis of the necessity of site network spatial information.

[0066] The radiosonde validation results (RMSE = 2.12 mm, r = 0.979) place PWVGAT's performance within the accuracy range of existing GNSS PWV products and radiosonde validation reports. PWVGAT infers PWV in locations lacking direct GNSS observations, yet achieves accuracy comparable to co-station direct GNSS measurements. The -1.40 mm system negative bias is consistent with known differences in GNSS and radiosonde PWV, stemming from variations in sampling geometry, representativeness, and dry bias correction of the radiosonde humidity sensor. The station-by-station error gradient (from 1.35 mm for 72293 coastal stations to 3.49 mm for 72393 inland stations) is physically plausible. Coastal stations benefit from dense surrounding GNSS coverage and relatively uniform oceanic air mass conditions, while inland stations in complex terrain face sharper water vapor gradients, making sparse observations more difficult to distinguish. This spatial error pattern suggests that the primary limiting factor for the current method is station density rather than model capacity.

[0067] The reconstructed 2D field exhibits spatial characteristics consistent with known California climatology: a strong coastal-inland PWV gradient, moisture accumulation in the central valley, and a topographic dehumidification effect from the Sierra Nevada Mountains. Comparison with the ERA5 TCWV (RMSE=2.89mm, r=0.970) shows that PWVGAT generates a spatially coherent field from an irregular network of stations without requiring reanalysis of the explicit physical constraints (mass conservation, hydrostatic equilibrium) employed by the system. The sensitivity of the elevation assignment scheme reveals the importance of consistent elevation representation. The ERA5 topography achieves the best statistical consistency because it shares the same topographic surface with the ERA5 TCWV reference field; while the Copernicus GLO-30 DEM provides higher resolution topographic information, it introduces scale mismatch when compared with the 0.25° ERA5 product. For applications requiring higher resolution reconstruction (e.g., 0.05°), the GLO-30 DEM is more suitable, but validation based on independent high-resolution reference data (e.g., dense radiosonde profiles or GNSS occultation profiles) is necessary.

[0068] The PWVGAT proposed in this embodiment is the first to apply graph attention networks to PWV spatial inference of GNSS station networks. By representing the GNSS station network as a graph with multidimensional node and edge features, PWVGAT learns anisotropic spatial attention weights, adapting to the irregular geometry, heterogeneous station density, and complex terrain of real-world networks. The masked self-supervised training strategy achieves inductive generalization to unseen locations, directly supporting two-dimensional field reconstruction without retraining.

[0069] The above are merely preferred embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A spatial inference method for precipitable water based on edge-enhanced graph attention network, characterized in that, Includes the following steps: A graph structure is constructed based on the geographic coordinates of multiple Global Navigation Satellite System (GNSS) stations. The graph structure includes nodes and edges, where the nodes correspond to GNSS stations and the edges represent the spatial relationships between the stations. Based on the geographic coordinates, meteorological reanalysis data, topographic elevation, periodic time coding, and precipitable water observations of each global navigation satellite system station, node characteristics of each node are generated. Based on the geographical distance, elevation difference, longitude difference, latitude difference, and normalized distance weight between the two stations, the edge features of the edge connecting the two nodes are generated. The node features and edge features are input into the edge-enhanced graph attention network. The precipitation information of neighboring nodes is aggregated through the attention mechanism to predict the precipitation of the target node. Based on the prediction results of each target node, a spatially continuous precipitation field is generated.

2. The method according to claim 1, characterized in that, The construction of the graph structure includes: Based on the geographical coordinates of each global navigation satellite system station, the set of neighboring nodes for each node is determined using a combination of K-nearest neighbor and distance threshold criteria, and edges connecting the node and its neighboring nodes are established.

3. The method according to claim 1, characterized in that, The generated node features include: Based on the geographic coordinates of each global navigation satellite system station, surface temperature and pressure, topographic elevation, periodic time coding, and precipitable water observation values ​​from meteorological reanalysis data, a multidimensional node feature vector is constructed.

4. The method according to claim 1, characterized in that, The features of generated edges include: A multidimensional edge feature vector is constructed based on the geographical distance, elevation difference, longitude difference, latitude difference, and normalized distance weight between the two stations.

5. The method according to claim 1, characterized in that, It also includes the mask self-supervised training step: In each training iteration, a portion of the precipitation observations from the Global Navigation Satellite System stations are randomly masked. The precipitation features of the masked stations are set to zero, while the features of other nodes are retained. The training objective is to predict the precipitation of the masked stations.

6. The method according to claim 5, characterized in that, Also includes: During the training process, the proportion of randomly occluded sites varies within a preset range.

7. The method according to claim 1, characterized in that, It also includes the steps for reconstructing the two-dimensional precipitation field: In the study area, regularly distributed virtual grid nodes are generated. Each virtual grid node is assigned geographic coordinates, topographic elevation, and surface air temperature and surface air pressure from meteorological reanalysis data. The precipitation characteristics of the virtual grid nodes are set to zero. Virtual grid nodes are added to the graph structure as new nodes and connected to nearby global navigation satellite system station nodes; The trained edge-enhanced graph attention network is applied to the expanded graph structure to predict the precipitation of all virtual grid nodes, generating a spatially continuous precipitation field.

8. The method according to claim 7, characterized in that, Adding virtual grid nodes to the graph structure includes: Based on the geographic coordinates of the virtual grid node and the preset K-nearest neighbor and distance threshold combination criteria, the set of neighbor nodes of the virtual grid node is determined, and the edges connecting the virtual grid node and its neighbor nodes are established.

9. The method according to claim 1, characterized in that, The edge-enhanced graph attention network comprises several stacked graph attention layers. Each layer uses a multi-head attention mechanism and calculates attention coefficients by concatenating edge features with the representations of the source and target nodes.

10. The method according to claim 1, characterized in that, The loss function of the edge-enhanced graph attention network includes the reconstruction loss of occluded nodes and a spatial smoothing regularization term based on edge weights, which are determined based on the geographical distance between nodes.

Images