Basin weather station correlation prediction method and system based on graph neural network modeling

Abstract

This invention relates to a method and system for watershed meteorological station correlation forecasting using graph neural network modeling, and belongs to the field of smart hydrological and meteorological forecasting technology. The method includes acquiring digital elevation model data, multi-source historical meteorological and hydrological data sequences, and meteorological radar volume scan data of the target area; obtaining a directed acyclic graph structure containing K sub-watershed units; constructing a meteorological background field and generating a three-dimensional radar reflectivity factor cube aligned with the grid; generating an initial embedded feature set for each station at all times within the historical time window; constructing a physical prior graph containing topological connection edges and topographically similar edges; generating a multi-step forecast sequence containing future rainfall, future temperature, and future outflow from each sub-watershed for each station; and outputting the trained graph neural network forecasting model. This invention embeds physical constraints such as watershed spatial topology, water balance, and temperature gradient into the model training process, effectively improving the accuracy and physical consistency of hydrological forecasts.

Description

Technical Field

[0001] This invention relates to the field of intelligent hydrological and meteorological forecasting technology, and more specifically, to a method and system for basin meteorological station correlation forecasting based on graph neural network modeling. Background Technology

[0002] With the intensification of global climate change and the frequent occurrence of extreme hydrological events, accurate watershed meteorological and hydrological forecasting has become a key technological support for flood control and disaster reduction, optimized water resource allocation, and efficient utilization of hydropower resources. Currently, watershed forecasting technologies both domestically and internationally have mainly evolved from empirical formula methods and conceptual hydrological models to distributed physical models, and in recent years, machine learning algorithms have been gradually introduced to assist in forecasting.

[0003] However, existing technologies generally suffer from the following problems: Traditional distributed hydrological models (such as SWAT and VIC models) can characterize physical mechanisms, but they are insufficient in representing the spatial heterogeneity of meteorological stations and are difficult to adapt to the characteristics of fragmented terrain and sparse stations in karst landform areas; Conventional machine learning methods (such as support vector machines and random forests) can establish statistical mapping relationships, but they treat each station as an independent sample and ignore the hydraulic connection between upstream and downstream of the watershed and the topographic barrier effect, resulting in significant forecast errors due to the asynchronous propagation of floods between upstream and downstream during rainstorms; Although existing graph neural networks have been introduced into fields such as transportation and social networking, they mostly use static adjacency matrices, which cannot adapt to the rapid evolution of weather systems in the southwest monsoon region and lack explicit constraints on physical laws such as water balance and energy conservation, often resulting in outliers in forecast results that violate physical common sense. In existing technologies, some researchers have attempted to use Kriging interpolation or inverse distance weighting methods for spatial interpolation, or to use long short-term memory networks for single-station time series prediction. However, these methods have failed to effectively integrate the nonlinear correlation between the watershed topology and the spatiotemporal evolution of meteorological elements, and the forecast accuracy and stability are difficult to meet the actual needs of joint operation of cascade hydropower station groups. Summary of the Invention

[0004] The purpose of this invention is to provide a method and system for basin meteorological station correlation forecasting based on graph neural network modeling, in order to improve the aforementioned problems. To achieve the above objective, the technical solution adopted by this invention is as follows: Firstly, this application provides a method for basin meteorological station correlation forecasting based on graph neural network modeling, including: Acquire digital elevation model data, multi-source meteorological and hydrological historical data sequences, and meteorological radar volume scan data of the target area; perform depression filling processing and eight-direction flow accumulation calculation on the digital elevation model data to obtain a directed acyclic graph structure containing K sub-basin units, where each sub-basin unit includes centroid latitude and longitude coordinates, average elevation, and catchment area attributes. Bilinear interpolation is used to map the observation values ​​of each monitoring station in the multi-source meteorological and hydrological historical data sequence to the raster grid of the digital elevation model data to construct the meteorological element background field; after performing ground clutter suppression and range attenuation correction on the meteorological radar volume scan base data, a three-dimensional radar reflectivity factor cube aligned with the raster grid is generated through polar coordinate to rectangular coordinate transformation and vertical linear interpolation. Based on the meteorological background field and the three-dimensional radar reflectivity factor cube, local spatial gradient features and vertical profile features are extracted for each monitoring station. The original observations, local spatial gradient features, vertical profile features, and station altitude, longitude, and latitude in the multi-source meteorological and hydrological historical data sequence are standardized and then vectorized to generate the initial embedded feature set for each station at all times within the historical time window. Based on the directed acyclic graph structure and the spatial coordinates and elevation of each station, a physical prior graph containing topological connection edges and terrain dynamic similarity edges is constructed. The edge weights of the physical prior graph are adaptively updated by dynamic time warping distance and cross-correlation peak coefficient, and new edges are dynamically added according to the temporal similarity threshold to form a dynamic graph structure that evolves over time. The initial embedded feature set and dynamic graph structure are input into a stacked graph-gated recurrent unit, which outputs the spatiotemporal hidden state sequence within the entire historical time window and the final encoder hidden state. Using the final encoder hidden state as the key matrix and value matrix, and the concatenation of the shared query starting embedding vector and the learnable position encoding vector as the query input of the decoder, they are fed into a decoder based on the Transformer architecture to generate a multi-step forecast sequence containing future rainfall, future temperature, and future outflow from each sub-basin. Based on the multi-step forecast sequence and the historical observations at the corresponding time, a hybrid loss function is constructed by combining the watershed water balance equation and the canyon temperature decay rate constraint along the river topology of the directed acyclic graph structure. The model parameters are iteratively updated until convergence, and the trained graph neural network forecast model is output.

[0005] Preferably, based on the meteorological background field and the three-dimensional radar reflectivity factor cube, local spatial gradient features and vertical profile features are extracted for each monitoring station; the original observations, local spatial gradient features, vertical profile features, and station altitude, longitude, and latitude from the multi-source meteorological and hydrological historical data sequence are standardized and then concatenated into vectors to generate an initial embedded feature set for each station at all times within the historical time window, including: Based on the meteorological background field, the central difference scheme of 3×3 grid window is used to calculate the temperature meridional gradient, temperature zonal gradient, pressure meridional gradient and pressure zonal gradient of each station to form a four-dimensional local spatial gradient vector. Among them, the temperature meridional gradient reflects the blocking effect of the east-west topographic barrier on heat transport, and the pressure meridional gradient is used to characterize the pressure gradient force driving mechanism of the canyon wind channel. Based on the three-dimensional radar reflectivity factor cube, reflectivity factor values ​​of 20 standard isobaric surface layers above each station are extracted along the vertical direction to form a vertical profile feature vector. Simultaneously, the convective effective potential energy of the station is calculated using a static approximation formula based on the virtual temperature difference, and the vertical wind shear is calculated based on the horizontal wind vector difference between the 200 hPa and 850 hPa isobaric surfaces. Among them, the convective effective potential energy is used to characterize the unstable energy reserve of atmospheric convection, and the vertical wind shear is used to characterize the modulation effect of the vertical wind shear on the organization of convection. The vertical profile features include the reflectivity factor of each isobaric surface layer and the convective effective potential energy and vertical wind shear derived therefrom. The original observation vectors, local spatial gradient vectors, vertical profile feature vectors, convective effective potential energy, vertical wind shear, and station altitude, longitude, and latitude in the multi-source meteorological and hydrological historical data sequence are Z-score standardized. After vector concatenation, a heterogeneous node initial feature vector with dimension 35 is generated. The initial feature vectors of all stations and all historical moments together constitute the initial embedded feature set.

[0006] Preferably, the physical prior graph, constructed based on the directed acyclic graph structure and the spatial coordinates and elevation of each station, includes topological connecting edges and terrain-dynamically similar edges, comprising: Based on the directed acyclic graph structure, directed topological connection edges are established for stations corresponding to sub-basins with river confluence relationships. The initial weight is set as the inverse of the river length. The direction of the directed topological connection edge is consistent with the river network confluence direction in the directed acyclic graph structure, which is used to characterize the water contribution relationship between upstream sub-basins and downstream sub-basins. Based on the station's altitude, slope, and normalized distance to the nearest watershed, the similarity of each pair of stations is calculated using the topographic dynamic similarity index formula. The topographic dynamic similarity index comprehensively considers the vertical temperature decrease effect caused by altitude differences, the difference in slope runoff velocity affected by slope differences, and the difference in rainfall-runoff response time reflected by the difference in watershed distance. When the similarity index is greater than 0.65, an undirected topographic dynamic similarity edge is established, with the initial weight set to the similarity index value; when the similarity index is less than or equal to 0.65, no topographic dynamic similarity edge is established. Directed topological connecting edges and undirected terrain-dynamic similarity edges are merged into a unified edge set. Combined with the set of station nodes, a physical prior graph containing topological connecting edges and terrain-dynamic similarity edges is constructed.

[0007] Preferably, the step of adaptively updating the edge weights of the physical prior graph through dynamic time-warped distance and cross-correlation peak coefficient, and dynamically adding new edges according to the temporal similarity threshold to form a dynamic graph structure that evolves over time, includes: Using the initial embedded feature set, the dynamic time warp distance of each pair of established edges on its historical feature sequence and the peak cross-correlation coefficient within the past 6-hour time window are calculated. The dynamic time warp distance is solved by dynamic programming algorithm to find the optimal time curvature path, which is used to measure the sequence similarity under nonlinear time distortion. The peak cross-correlation coefficient is determined by sliding window cross-correlation analysis to determine the maximum correlation coefficient and its corresponding lag time, which is used to capture the propagation time lag characteristics of the precipitation system. The physical prior weights of the physical prior graph and the data-driven dynamic similarity are weighted and fused using a dynamic weight fusion formula to obtain the final edge weights at each time step. The balance coefficient in the dynamic weight fusion formula is adaptively adjusted according to the characteristics of the flood and drought seasons in the watershed. During the flood season, data-driven similarity is emphasized to capture sudden associations of rainstorms, while during the drought season, physical prior weights are emphasized to maintain a stable topology. For node pairs without physical prior edges, if their dynamic time-warped distance is less than 0.2 times the historical window length and their cross-correlation peak coefficient is greater than 0.8, new edges are dynamically added, with weights directly obtained from the cross-correlation peak coefficient via Sigmoid mapping. The mechanism of dynamically adding new edges is used to identify strong correlations of meteorological elements temporarily formed during rainstorms and to supplement the strong coupling relationships within rainstorm clouds that are not covered by the fixed topology. All edge weights are organized into an adjacency matrix according to time steps, and combined with the fixed node set, the dynamic graph structure that evolves over time is generated.

[0008] Preferably, the graph-gated recurrent unit, which inputs the initial embedded feature set and the dynamic graph structure into a stacked graph, outputs the spatiotemporal hidden state sequence within the entire historical time window and the final encoder hidden state, including: The feature matrix from the initial embedded feature set and the adjacency matrix with self-loops in the dynamic graph structure are input into the graph convolutional layer. The node features are propagated layer by layer using the degree-normalized convolution formula. After passing through the linear rectified activation function, the node feature matrix of the current graph convolutional layer is output. The adjacency matrix with self-loops is the sum of the edge weight matrix and the identity matrix of the dynamic graph structure, which is used to preserve the node's own feature information. The degree-normalized convolution formula suppresses the feature drowning effect of high-connectivity nodes by symmetrically normalizing the degree matrix, and the elements of the degree diagonal matrix are the sum of the corresponding rows of the node, which is used to adaptively adjust the weight distribution of neighborhood aggregation. The graph convolution operation is embedded in the three types of gating structures of the graph gated recurrent unit: update gate, reset gate, and candidate hidden state. Each gating structure takes the concatenation of the node feature matrix at the current time step and the hidden state matrix at the previous time step as input. After graph convolution, it is activated by the Sigmoid or hyperbolic tangent activation function, and the hidden state matrix at each time step is recursively calculated by the final hidden state update formula. The update gate is used to control the retention ratio of historical hidden states, the reset gate is used to adjust the influence of historical hidden states on candidate hidden states, and the candidate hidden states are used to generate the latent feature representation at the current time step. After traversing all historical time steps, the spatiotemporal hidden state sequence and the final encoder hidden state are output. The spatiotemporal hidden state sequence encodes the comprehensive feature representation of each station after integrating spatial neighborhood association and temporal evolution law.

[0009] Preferably, the process uses the final encoder hidden state as the key matrix and value matrix, and the concatenation of the shared query starting embedding vector and the learnable position encoding vector as the query input to the decoder. These are fed into a Transformer-based decoder to generate a multi-step forecast sequence containing future rainfall, future temperature, and future outflow from each sub-basin for each station. This sequence includes: For each future moment in the decoding process, a learnable position encoding vector is initialized and generated using sine and cosine functions to represent the relative positional relationship of the forecast time and enable the model to perceive the difference in forecast difficulty for different lead times. Its even-numbered and odd-numbered dimensions are generated using sine and cosine functions, respectively, and the dimension index is related to the model dimension. The position encoding vector is concatenated with the shared query start embedding vector as the initial query input of the decoder. An 8-head attention mechanism is adopted to divide the final encoder hidden state into 8 subspaces. Each subspace independently learns different field-station association patterns. By introducing a physical constraint mask in the attention weight calculation, specific attention heads can focus on the confluence direction or elevation similarity. At the same time, the final encoder hidden state is linearly transformed and used as the key matrix and value matrix. The query input is linearly transformed to obtain the query matrix. The attention output is calculated by the scaled dot product attention formula. The scaled dot product attention formula prevents gradient vanishing due to excessively large dot product values ​​by dimensional normalization. The outputs of multiple attention heads are concatenated and then linearly transformed, followed by a feedforward network, residual connections, and layer normalization to obtain the decoder hidden state at each future time step. The feedforward network adopts a two-layer fully connected structure with Gaussian error linear unit activation function embedded in the middle to enhance nonlinear fitting ability. The residual connection is used to alleviate the gradient decay problem of deep networks, and the layer normalization is used to stabilize the distribution drift during the training process. The hidden state of the decoder is input into the feedforward network to generate a forecast vector for each future time step. The output layer of the forecast vector uses a piecewise linear activation function. The rainfall output uses ReLU activation to ensure non-negativity, the temperature output uses linear activation to allow for positive and negative deviations, and the flow output uses Softplus activation to ensure positive values ​​and smooth differentiability. The dimension of the forecast vector is equal to twice the sum of the total number of stations and the total number of sub-basins. The first half of the vector contains the future rainfall and temperature forecasts for each station, and the second half contains the future outflow forecasts for each sub-basin. The forecast vectors for all time steps are arranged in chronological order to output a multi-step forecast sequence.

[0010] Preferably, the process involves constructing a hybrid loss function based on a multi-step forecast sequence and corresponding historical observations, combined with the watershed water balance equation and the canyon temperature decay rate constraint along the river channel topology of a directed acyclic graph structure, and iteratively updating the model parameters until convergence. This includes: Based on the predicted outflow of each sub-basin, the predicted rainfall at each station, and the historical observations at the corresponding time in the multi-step forecast sequence, the forecast error term is calculated using the mean square error formula. The mean square error formula uses differentiated weights for each forecast element. The weight of the rainfall forecast error increases with the extension of the forecast period, and the weight of the flow forecast error automatically increases during the flood peak period. The physical consistency loss term was calculated using the water balance constraint formula based on the Thiessen polygon area, runoff coefficient, and baseflow estimate of each station. The Thiessen polygon area was generated by the Voronoi diagram algorithm based on the spatial distribution of the stations to characterize the effective control range of each station. The runoff coefficient was determined by the lookup table method based on the land use type and soil texture to convert rainfall into effective runoff. The baseflow estimate was obtained by segmenting the measured flow using the receding water curve method to characterize the stable recharge of groundwater to river runoff. Based on the canyon temperature lapse rate along the river topology of the directed acyclic graph structure, the temperature gradient error term is calculated using the geodesic gradient consistency formula. The canyon temperature lapse rate is calculated using the environmental lapse rate formula based on the elevation to constrain the spatial continuity of temperature forecast. The geodesic gradient consistency formula calculates the temperature difference between adjacent sub-basins along the shortest path of the river to penalize abnormal gradients that violate topographic and climatic patterns. The total loss is obtained by weighting and summing the forecast error term, water balance constraint loss term, and geodesic gradient consistency loss term. The weight coefficients are dynamically adjusted according to the convergence of the training phase. In the early stage, the forecast error term is emphasized to quickly fit the observed data, and in the later stage, the physical consistency loss term is emphasized to strengthen the satisfaction of physical constraints. An adaptive moment estimator optimizer is used to iteratively update the model parameters until convergence.

[0011] Preferably, the process of filling depressions and calculating eight-directional flow accumulation in the digital elevation model data yields a directed acyclic graph structure containing K sub-basin units, including: Based on the digital elevation model data, depression filling is used to eliminate false depressions. Depression filling is achieved by identifying local minimum elevation points and raising them to overflow elevations to eliminate water traps in the terrain data and ensure the continuity of river network extraction. The flow accumulation amount of each grid cell is calculated using an eight-directional flow accumulation algorithm, with a threshold of 1000 grid cells. The algorithm determines the flow direction based on the steepest slope principle and accumulates the number of grid cells in the upstream contribution area to identify potential river channel locations. The threshold is adaptively adjusted based on the watershed area and resolution: a high threshold is used for large watersheds to extract the main channel, while a low threshold is used for small watersheds to preserve tributary details. Grid cells exceeding the threshold are marked as river channel grid cells, and the river channel grid cells are vectorized to generate a river network, converting grid line features into vector line segments and establishing the topological connection relationship of the river network. The Stellar River Network Classification Algorithm is used to classify the river network. The Stellar River Network Classification Algorithm determines the river level based on the tributary inflow rules. The source river without tributary inflow is classified as Level 1. The level is incremented by one when two rivers of the same level converge. When rivers of different levels converge, the higher level is taken to represent the level of river inflow contribution. The source sub-basin of Level 1 river is used as the smallest unit. Upstream and downstream connections between sub-basins are established through the river network inflow relationship to form a directed acyclic graph structure. Each node corresponds to a sub-basin and includes the centroid latitude and longitude coordinates, average elevation and catchment area of ​​the sub-basin. Each directed edge represents the inflow direction from the upstream sub-basin to the downstream sub-basin.

[0012] Preferably, after outputting the trained graph neural network prediction model, the method further includes a deployment and execution step, which includes: The trained graph neural network prediction model is solidified in a serialized format that supports dynamic graph structures and loaded into the server's graphics processor memory. Then, a pipelined parallel strategy is adopted to divide data acquisition, feature extraction, graph structure update and model inference into four pipeline stages and bind them to independent CPU cores. The feature extraction pipeline stage runs on a general-purpose CPU to perform normalization and vector concatenation operations, while the model inference pipeline stage runs on a high-performance graphics processor to perform graph convolution and attention calculations. A forecasting task is triggered every 15 minutes using a sliding window approach. Each forecast reads rolling historical data from the past 72 hours and performs forward inference to generate hourly forecasts for the next 24 hours. Then, the outflow from each sub-basin in the multi-step forecast sequence is converted into the inflow from each cascade hydropower station through a pre-calibrated river confluence model. The pre-calibrated river confluence model is based on the linear assumption of the relationship between river storage and outflow and is calibrated using parameter estimation methods to determine the storage constant and flow proportion coefficient. The data is overlaid on the watershed digital twin platform in the form of a geographic information system layer. When the forecasted inflow exceeds the safety threshold corresponding to the flood control limit, an orange warning signal is automatically generated. When the forecasted inflow exceeds the design flood level, a red warning signal is automatically generated. When the forecasted inflow does not exceed the safety threshold corresponding to the flood control limit, the normal forecast status is maintained and monitoring continues. When the forecasted inflow is between the flood control limit and the design flood level, a blue attention prompt is generated to remind dispatchers to strengthen tracking and is pushed to the dispatch duty terminal via 5G message.

[0013] Secondly, this application also provides a basin meteorological station correlation forecasting system modeled by graph neural networks, including: Acquisition module: used to acquire digital elevation model data, multi-source meteorological and hydrological historical data sequences, and meteorological radar volume scan data of the target area; to perform depression filling processing and eight-direction flow accumulation calculation on the digital elevation model data, and to obtain a directed acyclic graph structure containing K sub-basin units, where each sub-basin unit includes centroid latitude and longitude coordinates, average elevation, and catchment area attributes; The construction module is used to map the observation values ​​of each monitoring station in the multi-source meteorological and hydrological historical data sequence to the raster grid of the digital elevation model data using bilinear interpolation, and construct the background field of meteorological elements; after performing ground clutter suppression and range attenuation correction on the meteorological radar volume scan base data, it generates a three-dimensional radar reflectivity factor cube aligned with the raster grid through polar coordinate to rectangular coordinate transformation and vertical linear interpolation. The generation module is used to extract local spatial gradient features and vertical profile features for each monitoring station based on the meteorological background field and the three-dimensional radar reflectivity factor cube. After standardizing the original observation values, local spatial gradient features, vertical profile features, and station altitude, longitude, and latitude in the multi-source meteorological and hydrological historical data sequence, the module performs vector concatenation to generate the initial embedded feature set for each station at all times within the historical time window. Update module: Based on the directed acyclic graph structure and the spatial coordinates and elevation of each station, it constructs a physical prior graph containing topological connection edges and terrain dynamic similarity edges; it adaptively updates the edge weights of the physical prior graph through dynamic time warping distance and cross-correlation peak coefficient, and dynamically adds new edges according to the temporal similarity threshold to form a dynamic graph structure that evolves over time. Output module: This module takes the initial embedded feature set and dynamic graph structure as input to a stacked graph-gated recurrent unit, outputting the spatiotemporal hidden state sequence within the entire historical time window and the final encoder hidden state. Using the final encoder hidden state as the key and value matrices, and the concatenation of the shared query starting embedding vector and the learnable positional encoding vector as the decoder's query input, this module feeds into a Transformer-based decoder to generate a multi-step forecast sequence containing future rainfall, future temperature, and future outflow from each sub-basin. Based on the multi-step forecast sequence and corresponding historical observations, a hybrid loss function is constructed by combining the basin water balance equation with the canyon temperature decay rate constraint along the river topology of the directed acyclic graph structure. The model parameters are iteratively updated until convergence, outputting the trained graph neural network forecast model.

[0014] Thirdly, this application also provides a basin meteorological station correlation forecasting device based on graph neural network modeling, comprising: Memory, used to store computer programs; A processor is used to implement the steps of the graph neural network-based watershed meteorological station association forecasting method when executing the computer program.

[0015] Fourthly, this application also provides a readable storage medium storing a computer program, which, when executed by a processor, implements the steps of the above-described method for basin meteorological station association forecasting based on graph neural network modeling.

[0016] The beneficial effects of this invention are as follows: This invention first extracts the river network topology of a watershed based on digital elevation model (DEM) data. It then constructs a directed acyclic graph (DAG) structure with sub-watersheds as nodes using depression filling, an eight-direction flow accumulation algorithm, and the Stradivarius hierarchy algorithm, explicitly encoding the watershed spatial topology information into the model. Based on this, a decoder network based on a transformer architecture is constructed. For each future time step in the decoding process, a learnable location encoding vector is initialized. This location encoding vector is concatenated with a shared query start embedding vector as the initial query input for the decoder. A multi-head attention mechanism is used to divide the encoder's hidden state into multiple subspaces, capturing different types of field-station association patterns. Attention output is calculated through scaled dot product attention, and after passing through a feedforward network, residual connections, and layer normalization, the decoder's hidden state is obtained, ultimately outputting a multi-step forecast sequence. Furthermore, this invention designs a hierarchical composite loss function that weights and sums the forecast error term based on differentiated weights, the physical consistency loss term based on water balance constraints, and the temperature gradient error term based on geodesic gradient consistency. An adaptive moment estimation optimizer is then used to iteratively update the model parameters. While maintaining the high accuracy of the data-driven model, physical constraints such as watershed spatial topology, water balance, and temperature gradient are embedded into the model training process, effectively improving the accuracy and physical consistency of hydrological forecasts.

[0017] Other features and advantages of the invention will be set forth in the following description, and will be apparent in part from the description, or may be learned by practicing embodiments of the invention. The objects and other advantages of the invention may be realized and obtained by means of the structures particularly pointed out in the written description, claims, and drawings. Attached Figure Description

[0018] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation on the scope. For those skilled in the art, other related drawings can be obtained based on these drawings without creative effort.

[0019] Figure 1 This is a schematic diagram of the process of basin meteorological station association forecasting method based on graph neural network modeling as described in this embodiment of the invention; Figure 2 This is a schematic diagram of the structure of the watershed meteorological station association forecasting system modeled by graph neural network as described in this embodiment of the invention; Figure 3 This is a schematic diagram of the structure of the basin meteorological station association forecasting equipment modeled by graph neural network as described in this embodiment of the invention.

[0020] In the diagram: 701, Acquisition module; 702, Construction module; 703, Generation module; 704, Update module; 705, Output module; 800, Watershed meteorological station association forecasting equipment modeled by graph neural network; 801, Processor; 802, Memory; 803, Multimedia component; 804, I / O interface; 805, Communication component. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0022] It should be noted that similar reference numerals and letters in the following figures indicate similar items; therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures. Furthermore, in the description of this invention, terms such as "first," "second," etc., are used only to distinguish descriptions and should not be construed as indicating or implying relative importance.

[0023] Example 1:

[0024] This embodiment provides a method for basin meteorological station correlation forecasting based on graph neural network modeling.

[0025] See Figure 1 The figure shows that the method includes steps S100, S200, S300, S400 and S500.

[0026] S100: Acquire digital elevation model data, multi-source meteorological and hydrological historical data sequences, and meteorological radar volume scan data of the target area; perform depression filling processing and eight-direction flow accumulation calculation on the digital elevation model data to obtain a directed acyclic graph structure containing K sub-basin units, where each sub-basin unit includes centroid latitude and longitude coordinates, average elevation, and catchment area attributes.

[0027] It is understood that step S100 includes S101, S102, and S103, wherein: S101. Based on the geographic latitude and longitude coordinates and projection parameters of the target area, perform spatial registration and coordinate system unification processing of multi-source heterogeneous data to generate spatiotemporal benchmark aligned digital elevation model data, multi-source meteorological and hydrological historical data sequences, and meteorological radar volume scan base data. The spatial registration adopts an affine transformation and resampling algorithm based on geographic control points to uniformly transform the digital elevation model data under the National 2000 coordinate system, the meteorological and hydrological telemetry station data under the WGS-84 coordinate system, and the radar volume scan base data under the polar coordinate system to the same projection coordinate system, eliminating the station location deviation and terrain obstruction misjudgment caused by the difference in coordinate systems. S102. Based on the spatiotemporal reference aligned digital elevation model data, a depression filling process is performed to eliminate false depressions and an eight-direction flow accumulation algorithm is used to calculate the runoff accumulation of each grid cell, constructing a topographic flow direction matrix containing the river network topology. The depression filling process uses a priority queue-based inundation algorithm to identify local minimum elevation points and raise them to overflow elevations, eliminating the topographic data runoff traps caused by sinkholes and funnels common in karst landform areas. The eight-direction flow accumulation algorithm determines the flow direction of each grid cell based on the steepest slope principle, accumulates the number of grid cells in the upstream contribution area, identifies potential river locations, and generates a runoff accumulation distribution. S103. Based on the topographic flow direction matrix containing the river network topology, the Stradivarius River Network Classification Algorithm is used to classify the river network and extract the source sub-basins as node units, establishing a directed acyclic graph structure containing K sub-basin units, where each sub-basin unit includes centroid latitude and longitude coordinates, average elevation, and catchment area attributes; the Stradivarius River Network Classification Algorithm determines the river level according to the tributary inflow rules, the source river without tributary inflow is level one, the level is incremented by one after two rivers of the same level converge, and the higher level is taken after rivers of different levels converge; by setting a threshold for the cumulative inflow, the main river is extracted and vectorized, and the source area of ​​the level one river is taken as the smallest unit to establish the upstream and downstream topological connection relationship between sub-basins, forming a directed acyclic graph structure to represent the inflow direction of the basin.

[0028] It should be noted that step S103 includes S1031, S1032, and S1033, wherein: S1031. Based on the digital elevation model data, a depression filling process is adopted to eliminate false depressions. The depression filling process identifies local minimum elevation points and raises them to overflow elevations to eliminate water traps in the terrain data and ensure the continuity of river network extraction. It should be noted that the algorithm used in this step starts from the DEM raster boundary, adding boundary pixels to a min-heap priority queue, with the queue key being the pixel elevation. During iteration, the pixel with the lowest elevation is popped from the queue each time, and its unvisited neighboring pixels are checked: if the neighboring pixel's elevation is lower than the current pixel's, it is identified as a depression, its elevation is raised to the current pixel's elevation (i.e., overflow elevation), and it is marked as processed; if the neighboring pixel's elevation is higher than or equal to the current pixel's, it is directly added to the queue. This process continues until all pixels have been processed, outputting a depression-free DEM. In karst watershed applications, this step is particularly important: a DEM without depression filling can cause the eight-direction flow accumulation algorithm to get stuck in local minima, resulting in broken river networks or incorrect confluence paths—for example, misidentifying the entrance of an underground river as the end point of surface confluence, or identifying an isolated sinkhole as a watershed outlet. By raising pseudo-depressions to overflow elevation, the algorithm ensures that surface water can continue to flow downstream along the terrain gradient, laying a continuous topographic foundation for subsequent river network extraction.

[0029] S1032. The flow accumulation of each grid cell is calculated using an eight-direction flow accumulation algorithm, and the flow accumulation threshold is set to 1000 grid cells. The eight-direction flow accumulation algorithm determines the flow direction based on the steepest slope principle and accumulates the number of grid cells in the upstream contribution area to identify potential river channel locations. The threshold is adaptively adjusted according to the watershed area and resolution: a high threshold is used for large watersheds to extract the main channel, and a low threshold is used for small watersheds to preserve tributary details. Grid cells exceeding the threshold are marked as river channel grid cells, and the river channel grid cells are vectorized to generate a river network, so as to convert the grid line features into vector line segments and establish the topological connection relationship of the river network. It should be noted that after calculating the gradient in all directions, if multiple directions have the same maximum gradient, the positive directions of east, south, west, and north are preferred over the diagonal directions to align with the physical law that actual runoff tends to flow along the slope line. The calculation of flow accumulation is either traced backward from the watershed outlet or accumulated forward from the source. The forward accumulation algorithm is as follows: starting from all grids, the flow is transferred along the flow direction to downstream grids, until the accumulation of each grid equals the total number of upstream grids (assuming each grid initially has 1 flow). The runoff accumulation threshold determines which grids are considered channels: grids with accumulation exceeding the threshold are considered permanent channel locations. A threshold of 1000 grid cells is an empirical benchmark, corresponding to the minimum catchment area required to form stable runoff at DEM resolution. The innovation of this step lies in the adaptive threshold adjustment mechanism: for large karst basins with an area exceeding 1000 square kilometers, a high threshold (e.g., 1000 grid cells) is used to extract the main channel, avoiding tributary noise from submerging the main channel morphology; for small karst spring basins with an area less than 100 square kilometers, a low threshold (e.g., 100 grid cells) is used to retain minor tributaries and capture the dispersed confluence characteristics of karst fissure development areas. This adaptive strategy solves the limitations of applying fixed thresholds to basins of different scales—high thresholds in small basins lead to excessively sparse river networks, missing key slope confluence paths; low thresholds in large basins generate a large number of parallel channels, interfering with the identification of the main channel. Grids exceeding the threshold, after being marked as river channel grids, need to be vectorized. This process uses boundary tracing algorithms (e.g., Moore-Neighbor Tracing) to extract the river centerline, converting discrete grids into continuous vector line segments, and establishing topological connections between river segments (start-end coordinates, segment length, upstream and downstream adjacency relationships). Vectorized river networks not only significantly reduce data volume, but also provide a precise topological basis for subsequent Straler classification and sub-basin delineation.

[0030] S1033. The Stellar River Network Classification Algorithm is used to classify the river network. The Stellar River Network Classification Algorithm determines the river level based on the tributary inflow rules. The source river without tributary inflow is classified as Level 1. The level is increased by one after two rivers of the same level converge. After rivers of different levels converge, the higher level is taken to represent the level of river inflow contribution. The source sub-basin of Level 1 river is taken as the smallest unit. The upstream and downstream connections between sub-basins are established through the river network inflow relationship to form a directed acyclic graph structure. Each node corresponds to a sub-basin and contains the centroid latitude and longitude coordinates, average elevation and catchment area of ​​the sub-basin. Each directed edge represents the inflow direction from the upstream sub-basin to the downstream sub-basin.

[0031] It should be noted that the construction of the directed acyclic graph (DAG) structure involves the following key operations: First, the source grids of all first-order rivers are identified and used as seed points to trace the watershed boundaries and delineate the spatial extent of each sub-watershed; second, the geometric attributes of each sub-watershed are calculated—centroid latitude and longitude coordinates (for subsequent spatial interpolation of meteorological data), average elevation (for temperature lapse rate calculation), and catchment area (for water balance constraints); finally, directed edges are established between sub-watersheds based on the river network topology, with the edge direction strictly following the confluence direction (from source to outlet) to ensure that there are no loop paths in the graph (i.e., acyclicity), which perfectly matches the unidirectional confluence characteristics of natural watersheds. Unlike conventional grids or random graphs, this DAG explicitly encodes the hierarchical confluence relationships of karst watersheds: nodes correspond to sub-watershed units with physical boundaries, edges correspond to quantifiable water transfer channels, node attributes contain key parameters of topographic-meteorological coupling, and the directionality of edges constrains the hydraulic logic of information propagation. This structure enables the message passing mechanism of graph neural networks to proceed along the actual confluence direction, rather than through disordered global aggregation, thereby embedding the physical understanding of watershed hydrology into a deep learning architecture.

[0032] S200: Using bilinear interpolation, the observed values ​​of each monitoring station in the multi-source meteorological and hydrological historical data sequence are mapped to the raster grid of the digital elevation model data to construct the meteorological element background field; After performing ground clutter suppression and range attenuation correction on the meteorological radar volume scan base data, a three-dimensional radar reflectivity factor cube aligned with the raster grid is generated through polar coordinate to rectangular coordinate transformation and vertical linear interpolation.

[0033] It is understood that step S200 includes S201, S202, and S203, wherein: S201. Based on the raster grid resolution and spatial range of the digital elevation model data, a bilinear interpolation algorithm is implemented to map the temperature, air pressure, humidity, wind speed and rainfall observations of each monitoring station in the multi-source meteorological and hydrological historical data sequence to the raster grid nodes, generating a spatially continuous meteorological element background field; the bilinear interpolation algorithm performs weighted averaging based on the distance weight between the station and its four adjacent grid nodes, eliminating the element field break caused by uneven spatial distribution of stations, and constructing a meteorological element spatial distribution matrix consistent with the terrain resolution; It should be noted that bilinear interpolation uses a weighted average based on the distance weights between the station and its four adjacent raster nodes, mapping discrete station observations to a continuous grid with the same resolution as the DEM. Considering the significant elevation differences in karst landforms, elevation corrections are applied to elements such as temperature and air pressure before interpolation. The environmental lapse rate is used to extrapolate to the target elevation before horizontal interpolation, ensuring the results conform to topographic and climatic patterns. The generated five-dimensional meteorological background field (temperature, air pressure, humidity, wind speed, and rainfall) is perfectly aligned with the DEM raster, establishing a direct correspondence between meteorological elements and topographic attributes. This step eliminates the element field fragmentation caused by uneven spatial distribution of stations, constructing a meteorological element spatial distribution matrix consistent with the topographic resolution, supporting the extraction of topographic-meteorological coupled features.

[0034] S202. Based on the original polar coordinate format of the meteorological radar volume scan base data, perform ground clutter suppression and range attenuation correction processing to generate quality-controlled radar reflectivity radial data. The ground clutter suppression uses a fuzzy logic algorithm based on reflectivity horizontal texture features to identify and remove fixed ground object echoes. The range attenuation correction uses the radar meteorological equation to compensate for the energy attenuation of the beam with increasing distance, restoring the true value of precipitation reflectivity in distant areas. It should be noted that ground clutter suppression employs a fuzzy logic algorithm based on reflectivity horizontal texture features. This algorithm utilizes multi-dimensional features such as spatial variance, temporal variance, and adjacency correlation to distinguish ground object echoes from precipitation echoes, eliminating false signals from fixed ground objects. Distance attenuation correction uses radar meteorological equations, estimating path integral attenuation based on observations of unsaturated areas at close range, and then performing exponential compensation for distant areas to restore the true value of precipitation reflectivity. This step addresses the observational biases caused by ground clutter and distance attenuation, providing reliable precipitation structure data for the construction of the three-dimensional radar reflectivity factor cube.

[0035] S203. Using quality-controlled radial radar reflectivity data as input, polar coordinate to rectangular coordinate transformation and vertical linear interpolation are performed sequentially to establish a three-dimensional radar reflectivity factor cube aligned with the grid. The polar coordinate to rectangular coordinate transformation uses nearest neighbor interpolation to map the azimuth-range format radar data to a Cartesian grid, and the vertical linear interpolation interpolates the data of irregular height layers to a standard isobaric surface layer, forming a three-dimensional radar reflectivity factor cube with a unified spatial reference, providing data support for subsequent vertical profile feature extraction.

[0036] It should be noted that the polar coordinate to Cartesian coordinate conversion uses nearest-neighbor interpolation to map the azimuth-range format radar data to a Cartesian grid, limiting the maximum conversion distance to control geometric distortion. After converting the radar slant range to altitude using vertical linear interpolation, linear interpolation is performed on the reflectivity of adjacent altitude layers at horizontal grid points to unify it to the standard isobaric surface layer. The generated 3D radar reflectivity factor cube is consistent with the spatial reference of the meteorological background field, forming a dual-field overlay structure to capture precipitation vertical structure information. This achieves spatial reference unification between polar coordinate radar data and the Cartesian grid meteorological field, providing 3D precipitation structure data support for vertical profile feature extraction and station-level convection parameter calculation.

[0037] S300: Based on the meteorological background field and the three-dimensional radar reflectivity factor cube, local spatial gradient features and vertical profile features are extracted for each monitoring station. The original observations, local spatial gradient features, vertical profile features, and station altitude, longitude, and latitude in the multi-source meteorological and hydrological historical data sequence are standardized and then vectorized to generate the initial embedded feature set for each station at all times within the historical time window.

[0038] It is understood that step S300 includes S301, S302, and S303, wherein: S301. Based on the meteorological background field, the temperature meridional gradient, temperature zonal gradient, pressure meridional gradient, and pressure zonal gradient of each station are calculated using the central difference scheme of a 3×3 grid window, forming a four-dimensional local spatial gradient vector. Among them, the temperature meridional gradient reflects the blocking effect of the east-west topographic barrier on heat transport, and the pressure meridional gradient is used to characterize the pressure gradient force driving mechanism of the canyon wind channel. It should be noted that this step is based on the meteorological background field and uses a central difference scheme with a 3×3 grid window to calculate the spatial gradients of temperature and pressure around each station. The core of the central difference scheme is to take the target station as the center, take the difference of elements between adjacent grids in the east-west (meridian) and north-south (zonal) directions, divide it by twice the grid spacing, and obtain the rate of change of the station in the four directions. Specifically, the meridional temperature gradient is the ratio of the temperature difference between adjacent grids in the east and west directions to twice the meridional distance, the zonal temperature gradient is the ratio of the temperature difference between adjacent grids in the north and south directions to twice the zonal distance, and the pressure gradient is calculated similarly. The meridional temperature gradient reflects the blocking effect of east-west topographic barriers on heat transport—when warm and humid airflow encounters north-south trending mountains, the windward slope is forced to rise and cool, while the leeward slope sinks and warms, forming a significant east-west temperature difference. This gradient quantifies the intensity of this topographic thermal effect. The meridional pressure gradient is used to characterize the pressure gradient force driving mechanism of canyon wind channels—the pressure difference between the deep canyons and plateau surfaces in karst areas drives the formation of canyon winds that alternate between day and night (such as the "narrowing effect" of the Wujiang Canyon). This gradient is directly related to the intensity and direction of local circulation.

[0039] The four-dimensional local spatial gradient vector (temperature meridional, temperature latitudinal, air pressure meridional, air pressure latitudinal) transforms the station from an isolated observation point into a field point with spatial context information, enabling the subsequent graph neural network to perceive the spatial variation trend of meteorological elements, rather than relying solely on single-point values.

[0040] S302. Based on the three-dimensional radar reflectivity factor cube, the reflectivity factor values ​​of 20 standard isobaric surface layers above each station are extracted along the vertical direction to form a vertical profile feature vector. At the same time, the convective effective potential energy of the station is calculated using the static approximation formula with virtual temperature difference, and the vertical wind shear is calculated based on the horizontal wind vector difference between the 200 hPa and 850 hPa isobaric surfaces. Among them, the convective effective potential energy is used to characterize the unstable energy reserve of atmospheric convection, and the vertical wind shear is used to characterize the modulation effect of the vertical wind shear on the organization of convection. The vertical profile features include the reflectivity factor of each isobaric surface layer and the convective effective potential energy and vertical wind shear derived therefrom. S303. The original observation vectors, local spatial gradient vectors, vertical profile feature vectors, convective available potential energy, vertical wind shear, and station elevation, longitude, and latitude from the multi-source meteorological and hydrological historical data sequence are Z-score standardized. After vector concatenation, a heterogeneous node initial feature vector with dimension 35 is generated. The initial feature vectors of all stations at all historical moments together constitute the initial embedded feature set. The formula for calculating convective available potential energy is as follows: In the formula, For convective effective potential energy, It is the acceleration due to gravity. To balance the height, For free convection height, The body is weak and cold. For the ambient temperature, Let the integral variable represent altitude. Altitude The differential variable.

[0041] It is important to note that the vertical profile characteristics not only include the original reflectivity value, but more importantly, two key convective parameters derived from it: convective available potential energy (CAPE) and vertical wind shear. CAPE is calculated using a static approximation formula based on the virtual temperature difference, integrating from the free convection height (LFC) to the equilibrium height (EL). It represents the ratio of the difference between the virtual temperature of the air parcel and the virtual temperature of the environment, characterizing the atmospheric convective instability energy reserve—the higher the CAPE value, the greater the buoyant acceleration experienced by the air parcel, and the stronger the convection development. In karst regions, the spatial distribution of CAPE often exhibits a characteristic of high values ​​in valleys and low values ​​in plateaus, closely related to the orographic lifting triggering mechanism. Strong vertical wind shear favors the tilted development of convective storms, separating rising and descending airflows and prolonging the storm's lifespan; under weak shear conditions, convection tends to develop vertically, easily leading to localized heavy precipitation. This parameter is of significant indicative importance for distinguishing between the short-duration heavy precipitation and persistent rainstorms common in karst regions. The three-dimensional radar reflectivity factor cube is compressed into a station-level vertical profile feature vector. Two key convective parameters, convective effective potential energy and vertical wind shear, are extracted to quantify the modulating effect of atmospheric convective instability energy reserves and wind field vertical shear on convective organization.

[0042] Furthermore, Z-score standardization is performed independently for each feature component, subtracting the mean of that component across all stations and all historical moments, and dividing by the standard deviation. This eliminates numerical scale differences caused by different physical units (temperature in °C, air pressure in hPa, reflectivity in dBZ, etc.), making each feature component comparable under a uniform scale of zero mean and unit variance. Vector concatenation connects the five standardized features in a fixed order, generating a 35-dimensional heterogeneous node initial feature vector (5-dimensional original observation + 4-dimensional spatial gradient + 20-dimensional vertical profile + 1-dimensional CAPE + 1-dimensional vertical wind shear + 3-dimensional static attributes). This dimension is fixed and independent of station location, ensuring that feature vectors from different stations and at different times can be directly input into the same neural network structure. By eliminating dimensional differences in multi-source heterogeneous features through Z-score standardization and generating fixed-dimensional heterogeneous node initial feature vectors through vector concatenation, a comprehensive node representation including ground observations, spatial gradients, upper-air structures, physical parameters, and static attributes is constructed, providing complete and standardized node attribute input for the graph neural network.

[0043] S400: Based on the directed acyclic graph structure and the spatial coordinates and altitude of each station, a physical prior graph containing topological connection edges and terrain dynamic similarity edges is constructed. The edge weights of the physical prior graph are adaptively updated by dynamic time warping distance and cross-correlation peak coefficient, and new edges are dynamically added according to the temporal similarity threshold to form a dynamic graph structure that evolves over time.

[0044] It is understood that in this step, S400 includes S401, S402, and S403, wherein: S401. Based on the directed acyclic graph structure, establish directed topological connection edges for stations corresponding to sub-basins with river confluence relationships. The initial weight is set as the inverse of the river length. The direction of the directed topological connection edge is consistent with the river network confluence direction in the directed acyclic graph structure, which is used to characterize the water contribution relationship between upstream sub-basins and downstream sub-basins. Specifically, each directed edge in the directed acyclic graph (DAG) structure is traversed to identify the corresponding upstream and downstream sub-basins, thereby locating the meteorological and hydrological stations belonging to these two sub-basins and establishing directed topological connections between these stations. The direction of the edges strictly follows the river network confluence direction, i.e., from the upstream station to the downstream station, consistent with the confluence direction in the DAG structure. Through a reciprocal transformation, the river length is converted into a weight value between 0 and 1; the shorter the length, the greater the weight, representing a close hydraulic connection; the longer the length, the smaller the weight, representing a looser hydraulic connection. This weight initialization method transforms the geographical distance metric into a graph structure association strength metric, providing a physically reasonable initial state for subsequent dynamic weight updates. By converting the river network confluence relationship into directed edges in the graph structure and initializing edge weights through the reciprocal of the river length, a topological constraint representing the direction of water contribution is established, providing the graph neural network with a message passing path consistent with hydrological principles.

[0045] S402. Based on the station's altitude, slope, and normalized distance to the nearest watershed, the similarity of each pair of stations is calculated using the topographic dynamic similarity index formula. The topographic dynamic similarity index comprehensively considers the vertical temperature decrease effect caused by altitude differences, the difference in slope runoff velocity affected by slope differences, and the difference in rainfall runoff response time reflected by the difference in watershed distance. When the similarity index is greater than 0.65, an undirected topographic dynamic similarity edge is established, with the initial weight set to the similarity index value. When the similarity index is less than or equal to 0.65, no topographic dynamic similarity edge is established. The three attributes correspond to different hydrological and meteorological processes: altitude differences lead to a vertical decrease in temperature, affecting precipitation phase and evaporation capacity; slope differences affect runoff velocity, with faster runoff on steep slopes and slower runoff on gentle slopes; and watershed distance differences reflect the timeliness of rainfall-runoff response, with faster response near the watershed and slower response far from the watershed. The topographic dynamic similarity index formula uses a negative exponential form, summing the standardized differences of the three attributes and then taking the exponent to output a similarity value between 0 and 1. The standardization process uses the watershed statistical standard deviation as a scale to eliminate the influence of different attribute dimensions, making altitude differences (meters), slope differences (degrees), and distance differences (dimensionless) comparable. The exponential transformation ensures that similarity monotonically decreases as differences increase, and the decay rate is naturally adjusted by the exponential function characteristics—similarity approaches 1 for small differences and rapidly approaches 0 for large differences. The similarity threshold is set at 0.65, a value based on empirical values ​​from watershed hydrological similarity studies. When the similarity index is greater than 0.65, the two stations are considered to have sufficient topographic dynamic similarity, and an undirected topographic dynamic similarity edge is established, with the initial weight set to the similarity index value. When the similarity index is less than or equal to 0.65, the similarity is considered insufficient, and no topographic dynamic similarity edge is established to avoid weakly associated edges interfering with message aggregation in the graph structure. This threshold screening mechanism ensures that the establishment of topographic dynamic similarity edges has a clear physical basis, rather than being a simple full connection. In other words, undirected topographic dynamic similarity edges and directed topological connection edges complement each other: the former, based on the similarity of spatial attributes, captures the synchronous response relationship between stations within the same climate zone (e.g., mountaintop stations with similar altitudes may be affected by clouds and fog simultaneously); the latter, based on the connectivity of river network topology, captures the hydraulic transmission relationship between upstream and downstream. The fusion of the two constructs a composite graph structure that includes both physical confluence constraints and spatial similarity associations.

[0046] S403. Merge the directed topological connection edges and the undirected terrain dynamic similarity edges into a unified edge set. Combine this with the site node set to construct a physical prior graph containing both topological connection edges and terrain dynamic similarity edges. The formula for calculating the terrain dynamic similarity index is as follows: In the formula, The terrain dynamic similarity index between station p and station q. , The elevations of stations p and q are respectively. , The slopes for stations p and q are respectively. and Let be the normalized distances from stations p and q to the nearest watershed, respectively. , , These represent the statistical standard deviations of altitude, slope, and watershed distance within the watershed.

[0047] It should be noted that the physical prior graph serves as the baseline framework for the subsequent dynamic graph structure evolution. The dynamic update process adjusts edge weights and adds or removes dynamic edges based on real-time observation data, while the topological skeleton (directed edges) and similar skeleton (undirected edges) remain stable, ensuring that the evolution of the graph structure is always constrained within a physically reasonable range.

[0048] It is understandable that step S400 also includes S404, S405, and S406, where: S404. Using the initial embedded feature set, calculate the dynamic time warp distance of each pair of established edges on its historical feature sequence and the peak cross-correlation coefficient within the past 6-hour time window; wherein, the dynamic time warp distance is solved by dynamic programming algorithm to find the optimal time curvature path, which is used to measure the sequence similarity under nonlinear time distortion; wherein, the peak cross-correlation coefficient is determined by sliding window cross-correlation analysis to determine the maximum correlation coefficient and its corresponding lag time, which is used to capture the propagation time lag characteristics of the precipitation system; It should be noted that in the correlation of meteorological stations in a watershed, the propagation of precipitation systems from upstream to downstream is often accompanied by time lag and velocity variations—for example, the difference in the advance velocity of cold air fronts in different topographic regions. DTW can capture the sequence similarity under this nonlinear time distortion, while traditional distance metrics will underestimate the correlation strength due to phase differences. The peak cross-correlation coefficient is determined through sliding window cross-correlation analysis. For each pair of nodes, with a fixed past 6-hour time window, the cross-correlation coefficient of one sequence relative to another is calculated at different time lags (-3 hours to +3 hours), and the maximum value is taken as the peak cross-correlation coefficient, with the corresponding time lag being the optimal lag time. This coefficient directly quantifies the linear correlation strength of meteorological elements between the two stations, while the optimal lag time reveals the propagation direction and velocity of the precipitation system—a positive lag indicates that the upstream leads the downstream, and a negative lag indicates that the downstream leads the upstream (possibly reflecting topographic waves or feedback mechanisms).

[0049] The two metrics complement each other: DTW captures overall morphological similarity and is robust to time distortion; cross-correlation captures the strength of linear correlation and reveals time-lag characteristics. Together, they constitute a data-driven dynamic similarity representation, providing input for subsequent weight fusion.

[0050] S405. The physical prior weights of the physical prior graph and the data-driven dynamic similarity are weighted and fused using a dynamic weight fusion formula to obtain the final edge weights at each time step. The balance coefficient in the dynamic weight fusion formula is adaptively adjusted according to the characteristics of the flood and drought seasons in the watershed. During the flood season, data-driven similarity is emphasized to capture sudden associations of rainstorms, while during the drought season, physical prior weights are emphasized to maintain a stable topology. It should be noted that the dynamic weight fusion formula adopts a linear weighting form, and the balance coefficient is adaptively adjusted according to the seasonal characteristics of flood and drought in the basin. During the flood season, rainfall is frequent and weather systems are active, and the association patterns between stations change rapidly. At this time, the weight of dynamic similarity is increased (e.g., 0.7) so that the model can capture temporary strong associations formed by sudden rainstorms. During the dry season, the weather is stable and the association patterns are fixed. At this time, the weight of physical prior weights is increased (e.g., 0.7) to maintain a stable topological structure and avoid false associations caused by data noise. The Sigmoid function maps the similarity scores before fusion to the interval between 0 and 1, realizing the normalization of weights. The peak cross-correlation coefficient of the numerator directly reflects the association strength, while the dynamic time warp distance of the denominator serves as an adjustment factor—the larger the DTW, the larger the denominator and the smaller the overall score, reflecting the penalty mechanism for time distortion. By using an adaptive balance coefficient to weight and fuse physical prior weights with data-driven similarity, the edge weights are dynamically adjusted according to the seasonal characteristics of flood and drought in the basin. During the flood season, it captures sudden associations caused by rainstorms, and during the dry season, it maintains a stable topology, realizing the adaptive response of the graph structure to the evolution of hydrological and meteorological conditions.

[0051] S406. For node pairs without physical prior edges, if their dynamic time-warped distance is less than 0.2 times the historical window length and their cross-correlation peak coefficient is greater than 0.8, then new edges are dynamically added, and the weights are directly obtained from the cross-correlation peak coefficients via Sigmoid mapping. The mechanism of dynamically adding new edges is used to identify strong correlations of meteorological elements temporarily formed during rainstorms and to supplement the strong coupling relationships within rainstorm clouds that are not covered by the fixed topology. All edge weights are organized into an adjacency matrix according to time steps, and combined with the fixed node set, a dynamic graph structure that evolves over time is generated.

[0052] It should be noted that the triggering conditions for dynamic edge addition are: the dynamic time-normalized distance is less than 0.2 times the historical window length (highly similar sequence morphology), and the peak cross-correlation coefficient is greater than 0.8 (extremely strong linear correlation). These two thresholds jointly filter out station pairs that are not covered by the physical prior map but whose data clearly indicate a strong correlation—typically two stations within a rainstorm cloud cluster. These stations may be located in different sub-basins, without direct river connections, but are influenced by the same convective system and exhibit highly synchronized changes in meteorological elements. The physical prior map is constructed based on topography and river networks, suitable for stable climate zoning and confluence relationships. Strong convective systems during rainstorms may cross topographic boundaries, forming temporary, non-topographically constrained strong correlations. Dynamic edge addition enables the model to capture these special correlation patterns under extreme weather conditions. Finally, all edge weights are organized into an adjacency matrix according to time steps, combined with a fixed node set, to generate a dynamic graph structure that evolves over time. This structure keeps the node set unchanged, but the edge weight matrix and edge existence change over time, forming a time-varying graph sequence. The graph structure at each time step is a snapshot of the associated state of the meteorological stations in the watershed at that moment, supporting the graph neural network to perform message passing and state updates along the time axis.

[0053] In other words, by dynamically adding temporary strong correlation edges during the rainstorm period through dual threshold screening, the internal coupling relationships of rainstorm clouds that are not covered by the fixed topology are supplemented. The time-varying edge weights are organized into an adjacency matrix sequence to generate a dynamic graph structure that evolves over time, providing time-varying topological input for the spatiotemporal graph neural network.

[0054] S500: The initial embedded feature set and dynamic graph structure are input into a stacked graph-gated recurrent unit, which outputs the spatiotemporal hidden state sequence within the entire historical time window and the final encoder hidden state. Using the final encoder hidden state as the key matrix and value matrix, and the concatenation of the shared query starting embedding vector and the learnable position encoding vector as the query input of the decoder, they are fed into the decoder based on the Transformer architecture to generate a multi-step forecast sequence containing future rainfall, future temperature, and future outflow from each sub-basin. Based on the multi-step forecast sequence and the historical observations at the corresponding time, a hybrid loss function is constructed by combining the watershed water balance equation and the canyon temperature decay rate constraint along the river topology of the directed acyclic graph structure. The model parameters are iteratively updated until convergence, and the trained graph neural network forecast model is output.

[0055] It is understood that in this step, S500 includes S501, S502, S503, S504, S505, S506, and S507, wherein: S501. Input the feature matrix from the initial embedded feature set and the self-loop adjacency matrix of the dynamic graph structure into the graph convolutional layer. Use the degree-normalized convolution formula to propagate the node features layer by layer. After passing through the linear rectified activation function, output the node feature matrix of the current graph convolutional layer. The self-loop adjacency matrix is ​​the sum of the edge weight matrix and the identity matrix of the dynamic graph structure, which is used to preserve the node's own feature information. The degree-normalized convolution formula suppresses the feature drowning effect of high-connectivity nodes by symmetrically normalizing the degree matrix. The elements of the degree diagonal matrix are the sum of the corresponding rows of the node, which is used to adaptively adjust the weight distribution of neighborhood aggregation. S502. The graph convolution operation is embedded into the three types of gating structures of the graph gated recurrent unit: update gate, reset gate, and candidate hidden state. Each gating structure takes the concatenation of the node feature matrix at the current time step and the hidden state matrix at the previous time step as input. After graph convolution, it is activated by the Sigmoid or hyperbolic tangent activation function, and the hidden state matrix at each time step is recursively calculated by the final hidden state update formula. The update gate is used to control the retention ratio of historical hidden states, the reset gate is used to adjust the influence of historical hidden states on candidate hidden states, and the candidate hidden states are used to generate the potential feature representation at the current time step. S503. After traversing all historical time steps, output the spatiotemporal hidden state sequence and the final encoder hidden state. The spatiotemporal hidden state sequence encodes the comprehensive feature representation of each station after integrating spatial neighborhood association and temporal evolution law.

[0056] It should be noted that this step inputs the feature matrix from the initial embedded feature set and the self-loop adjacency matrix of the dynamic graph structure into the graph convolutional layer. Layer-by-layer propagation and aggregation of node features are achieved through a degree-normalized convolution formula. The self-loop adjacency matrix is ​​the sum of the edge weight matrix and the identity matrix of the dynamic graph structure, used to preserve the node's own feature information. The degree-normalized convolution formula suppresses the feature overwhelming effect of high-connectivity nodes by symmetrically normalizing the degree matrix. The elements of the degree diagonal matrix are the sum of the corresponding rows of the self-loop adjacency matrix, used to adaptively adjust the weight distribution of neighborhood aggregation. After a nonlinear transformation introduced by a linear rectified activation function, the node feature matrix of the current graph convolutional layer is output. Subsequently, the graph convolution operation is embedded into three types of gating structures in the graph gating recurrent unit: update gate, reset gate, and candidate hidden state. The gating structure takes the concatenation of the current node feature matrix and the previous hidden state matrix as input. After graph convolution, it is activated by either the Sigmoid or hyperbolic tangent activation function. The hidden state matrix at each time step is recursively calculated using the final hidden state update formula. The update gate controls the retention ratio of historical hidden states, and the reset gate adjusts the influence of historical hidden states on candidate hidden states. The candidate hidden states are used to generate the potential feature representation of the current time step. After traversing all historical time steps, the spatiotemporal hidden state sequence and the final encoder hidden state are output. The spatiotemporal hidden state sequence encodes the comprehensive feature representation of each station after fusing spatial neighborhood association and temporal evolution. The final encoder hidden state compresses the information of the entire historical time window, providing hierarchical spatiotemporal representation input for the subsequent decoder.

[0057] It should be noted that a deployment execution step is included after step S500, which includes: The trained graph neural network prediction model is solidified in a serialized format that supports dynamic graph structures and loaded into the server's graphics processor memory. Then, a pipelined parallel strategy is adopted to divide data acquisition, feature extraction, graph structure update and model inference into four pipeline stages and bind them to independent CPU cores. The feature extraction pipeline stage runs on a general-purpose CPU to perform normalization and vector concatenation operations, while the model inference pipeline stage runs on a high-performance graphics processor to perform graph convolution and attention calculations. A forecasting task is triggered every 15 minutes using a sliding window approach. Each forecast reads rolling historical data from the past 72 hours and performs forward inference to generate hourly forecasts for the next 24 hours. Then, the outflow from each sub-basin in the multi-step forecast sequence is converted into the inflow from each cascade hydropower station through a pre-calibrated river confluence model. The pre-calibrated river confluence model is based on the linear assumption of the relationship between river storage and outflow and is calibrated using parameter estimation methods to determine the storage constant and flow proportion coefficient. The data is overlaid on the watershed digital twin platform in the form of a geographic information system layer. When the forecasted inflow exceeds the safety threshold corresponding to the flood control limit, an orange warning signal is automatically generated. When the forecasted inflow exceeds the design flood level, a red warning signal is automatically generated. When the forecasted inflow does not exceed the safety threshold corresponding to the flood control limit, the normal forecast status is maintained and monitoring continues. When the forecasted inflow is between the flood control limit and the design flood level, a blue attention prompt is generated to remind dispatchers to strengthen tracking and is pushed to the dispatch duty terminal via 5G message.

[0058] It should be noted that this step solidifies the trained graph neural network forecasting model into a serialized format that supports dynamic graph structures and loads it into the server's graphics processing unit (GPU) memory. Then, a pipelined parallel strategy is employed to divide data acquisition, feature extraction, graph structure update, and model inference into four pipeline stages, each bound to an independent CPU core. The feature extraction pipeline stage runs on a general-purpose CPU to perform normalization and vector concatenation operations, while the model inference pipeline stage runs on a high-performance GPU to perform graph convolution and attention calculations. This heterogeneous hardware collaboration achieves efficient offloading of computational tasks. A forecasting task is triggered every 15 minutes using a sliding window approach. Each forecast reads rolling historical data from the past 72 hours and performs forward inference to generate hourly forecasts for the next 24 hours. Subsequently, the sub-forecasts in the multi-step forecast sequence are... The outflow from the basin is converted into the inflow of each cascade hydropower station through a pre-calibrated river confluence model. The pre-calibrated river confluence model is based on the linear assumption of the relationship between river storage and outflow, and uses parameter estimation methods to calibrate the storage constant and flow proportion coefficient to realize the flood evolution calculation from the sub-basin outlet to the hydropower station inflow. Finally, it is displayed on the basin digital twin platform in the form of a geographic information system layer, and a graded early warning response is implemented based on the comparison results of the predicted inflow and multi-level safety thresholds: an orange warning signal is generated when the flood limit water level is exceeded, a red warning signal is generated when the design flood level is exceeded, normal monitoring is maintained when the flood limit water level is not exceeded, and a blue attention prompt is generated when the water level is between the two. At the same time, the data is pushed to the dispatch duty terminal via 5G message, forming a fully automated business chain from model inference, flood evolution to early warning issuance.

[0059] S504. For each future moment in the decoding process, a learnable position encoding vector is initialized and generated using sine and cosine functions to represent the relative positional relationship of the forecast time and enable the model to perceive the difference in forecast difficulty for different lead times. Its even-numbered and odd-numbered dimensions are generated using sine and cosine functions respectively, and the dimension index is related to the model dimension. The position encoding vector is concatenated with the shared query start embedding vector as the initial query input of the decoder. It should be noted that the location encoding vector is generated using sine and cosine functions, with its even-numbered and odd-numbered dimensions generated using sine and cosine functions respectively. The dimension index is related to the model dimension. In the watershed meteorological forecasting scenario, the forecasting difficulty varies significantly across different lead times: 0-6 hour nowcasts mainly rely on extrapolation of radar echoes, resulting in low uncertainty; 6-24 hour short-term forecasts are significantly more uncertain due to model initialization errors and boundary conditions. Location encoding, by assigning a unique location identifier to each future moment, enables the decoder to distinguish inputs from different lead times and adaptively adjust attention weights and output confidence. The location encoding vector is concatenated with a shared query start embedding vector, serving as the initial query input for the decoder. The query start embedding vector is a learnable parameter shared across all future moments, providing the initial semantic anchor for the forecasting task; the location encoding vector changes with time, providing relative information about temporal location. The concatenation of the two forms a complete query representation that includes both task semantics and temporal location.

[0060] S505 employs an 8-head attention mechanism to divide the final encoder hidden state into 8 subspaces. Each subspace independently learns different field-station association patterns. By introducing a physical constraint mask in the attention weight calculation, specific attention heads can focus on confluence direction or elevation similarity. Simultaneously, the final encoder hidden state is linearly transformed into the key matrix and value matrix. The query input is linearly transformed to obtain the query matrix, and the attention output is calculated using a scaled dot product attention formula. The scaled dot product attention formula uses dimension normalization to prevent gradient vanishing caused by excessively large dot product values. S506. The outputs of multiple attention heads are concatenated and then linearly transformed, followed by a feedforward network, residual connections, and layer normalization to obtain the decoder hidden state at each future time step. The feedforward network adopts a two-layer fully connected structure with Gaussian error linear unit activation function embedded in the middle to enhance nonlinear fitting ability. The residual connection is used to alleviate the gradient decay problem of deep networks, and the layer normalization is used to stabilize the distribution drift during the training process. It should be noted that the division of labor among attention heads reflects the diversity of meteorological correlations in the watershed: some attention heads focus on capturing the confluence-driven relationship between upstream and downstream stations—this type of correlation has a clear direction and lag, and the attention weights should decrease along the confluence direction; some attention heads focus on capturing the topographic wave propagation effect between stations at the same altitude—this type of correlation has spatial proximity and synchronicity, and the attention weights should be concentrated in areas of similar altitude. Through multi-head parallelism, the model can simultaneously learn correlation patterns driven by multiple physical mechanisms, avoiding the averaging tendency of single-head attention. Finally, the encoder hidden state, after linear transformation, serves as the key and value matrices, and the query input, after linear transformation, serves as the query matrix. The linear transformation projects the encoder hidden state and query input to different subspaces, giving keys, values, and queries independent expressive power. The scaled dot product attention formula calculates the dot product of the query and key, and normalizes it by dividing by the square root of the key vector dimension to prevent the Softmax gradient vanishing problem caused by excessively large dot product values. The attention output is a weighted sum of values, with weights obtained by Softmax normalization of the similarity between the query and key. By learning multiple field-station association modes such as confluence drive and topographic wave propagation in parallel with 8 heads, and calculating the similarity between query and key by scaling dot product attention and normalizing the weighted sum, the encoder's hidden state can be used to adaptively extract information from the decoder.

[0061] Understandably, the concatenation of multi-head outputs combines the attention results from eight subspaces into a complete vector, and the linear transformation projects it into a space consistent with the model's dimensions. The feedforward network employs a two-layer fully connected structure, embedding a Gaussian error linear unit activation function in between to enhance nonlinear fitting capabilities. Compared to ReLU, the Gaussian error linear unit has a smooth non-zero gradient in the negative region, alleviating the neuron death problem and making deep network training more stable. Residual connections add the input and output of the feedforward network, forming skip connections, mitigating the gradient decay problem in deep networks. During backpropagation, residual connections provide shortcuts for gradient flow, enabling efficient updates of underlying parameters. Layer normalization normalizes the feature dimensions of each sample, stabilizing distribution drift during training and reducing dependence on learning rate and initialization sensitivity. By enhancing nonlinear fitting capabilities through a two-layer feedforward network, mitigating gradient decay through residual connections, stabilizing the training distribution through layer normalization, and generating the decoder's hidden state for each future time step, the network achieves this.

[0062] S507. Input the hidden state of the decoder into the feedforward network to generate a forecast vector for each future time step. The output layer of the forecast vector uses a piecewise linear activation function. The rainfall output uses ReLU activation to ensure non-negativity, the temperature output uses linear activation to allow for positive and negative deviations, and the flow output uses Softplus activation to ensure positive values ​​and smooth differentiability. The dimension of the forecast vector is equal to twice the sum of the total number of stations and the total number of sub-basins. The first half of the vector contains the future rainfall forecast and the future temperature forecast for each station, and the second half contains the future outflow forecast for each sub-basin. Arrange the forecast vectors for all time steps in chronological order to output a multi-step forecast sequence.

[0063] It should be noted that the output layer of the forecast vector employs piecewise linear activation functions, selecting different activation methods for different physical constraints: the rainfall output uses ReLU activation to ensure non-negativity, consistent with the physical nature of precipitation; the temperature output uses linear activation to allow for positive and negative deviations, adapting to bidirectional temperature fluctuations; the flow output uses Softplus activation to ensure positive values ​​and smooth differentiability, avoiding the non-differentiability problem of ReLU at zero point, while maintaining the non-negativity constraint of flow. The dimension of the forecast vector is equal to twice the sum of the total number of stations and the total number of sub-basins. The first half contains the future rainfall and temperature forecasts for each station, and the second half contains the future outflow forecasts for each sub-basin. This output structure covers the core elements of basin meteorological and hydrological forecasting—meteorological elements at the station scale and hydrological responses at the sub-basin scale—supporting a complete forecast chain from precipitation to runoff. The forecast vectors for all times are arranged in chronological order, outputting a multi-step forecast sequence. This sequence is the model's complete prediction of the basin's state for the next 24 hours, including hourly temporal resolution, and can be directly used for scheduling decisions and early warning issuance.

[0064] By using piecewise linear activation functions to satisfy the physical constraints of different elements, a multi-step forecast sequence covering meteorological elements at stations and hydrological responses in sub-basins is generated, forming a complete 24-hour forecast of the basin's state at hourly resolution.

[0065] This step also includes: Based on the predicted outflow of each sub-basin, the predicted rainfall at each station, and the historical observations at the corresponding time in the multi-step forecast sequence, the forecast error term is calculated using the mean square error formula. The mean square error formula uses differentiated weights for each forecast element. The weight of the rainfall forecast error increases with the extension of the forecast period, and the weight of the flow forecast error automatically increases during the flood peak period. The physical consistency loss term was calculated using the water balance constraint formula based on the Thiessen polygon area, runoff coefficient, and baseflow estimate of each station. The Thiessen polygon area was generated by the Voronoi diagram algorithm based on the spatial distribution of the stations to characterize the effective control range of each station. The runoff coefficient was determined by the lookup table method based on the land use type and soil texture to convert rainfall into effective runoff. The baseflow estimate was obtained by segmenting the measured flow using the receding water curve method to characterize the stable recharge of groundwater to river runoff. Based on the canyon temperature lapse rate along the river topology of the directed acyclic graph structure, the temperature gradient error term is calculated using the geodesic gradient consistency formula. The canyon temperature lapse rate is calculated using the environmental lapse rate formula based on the elevation to constrain the spatial continuity of temperature forecast. The geodesic gradient consistency formula calculates the temperature difference between adjacent sub-basins along the shortest path of the river to penalize abnormal gradients that violate topographic and climatic patterns. The total loss is obtained by weighted summing of the forecast error term, the water balance constraint loss term, and the geodesic gradient consistency loss term. The weighting coefficients are dynamically adjusted based on the convergence progress during the training phase. Initially, the forecast error term is emphasized to quickly fit the observed data; later, the physical consistency loss term is emphasized to strengthen the satisfaction of physical constraints. An adaptive moment estimation optimizer is used to iteratively update the model parameters until convergence. The formula for calculating the water balance constraint loss term is as follows: In the formula, The loss term is the water balance constraint term, and F is the forecast duration. Indexing future moments For the total number of sub-basins, Let be the predicted outflow from the i-th sub-basin at time τ. M represents the time step, and M represents the total number of stations. Let be the predicted rainfall value for the j-th station at time τ. Let j represent the area of ​​the Thiessen polygon. Let be the runoff coefficient of the area where the j-th station is located. This is the estimated base current at time τ.

[0066] It should be noted that this step, based on the predicted outflow from each sub-basin, the predicted rainfall at each station, and the measured values ​​in the multi-step forecast sequence, uses the mean square error formula to calculate the forecast error term. The mean square error formula employs differentiated weights for each forecast element; the rainfall forecast error weight increases with the lead time to strengthen the penalty for uncertainties in long-term precipitation, while the flow forecast error weight automatically increases during peak flood periods to improve the fitting accuracy for critical flood control periods. Simultaneously, based on the Thiessen polygon area, runoff coefficient, and baseflow estimate at each station, the water balance constraint formula is used to calculate the physical consistency loss term. The Thiessen polygon area is generated using the Voronoi diagram algorithm based on the spatial distribution of the stations. To characterize the effective control range of each station, the runoff coefficient is determined by looking up tables based on land use type and soil texture to convert rainfall into effective runoff. The baseflow estimate is obtained by segmenting the measured flow using the receding water curve method to characterize the stable recharge of groundwater to river runoff. Based on the canyon temperature lapse rate along the river topology of the directed acyclic graph structure, the geodesic gradient consistency formula is used to calculate the temperature gradient error term. The canyon temperature lapse rate is calculated based on the elevation using the environmental lapse rate formula to constrain the spatial continuity of temperature forecasts. The geodesic gradient consistency formula calculates the temperature difference between adjacent sub-basins along the shortest path of the river to penalize abnormal gradients that violate topographic and climatic patterns.

[0067] The forecast error term, water balance constraint loss term, and geodesic gradient consistency loss term are then weighted and summed to obtain the total loss. The weight coefficients are dynamically adjusted according to the convergence of the training phase. In the early stage, the focus is on the forecast error term to quickly fit the observed data, and in the later stage, the focus is on the physical consistency loss term to strengthen the satisfaction of physical constraints. An adaptive moment estimation optimizer is used to iteratively update the model parameters until convergence, so that the model achieves a balance between data fitting and physical rationality, and outputs a trained graph neural network forecast model that matches the true observation value and conforms to the hydrological laws of the watershed.

[0068] In summary, this invention constructs a physical prior graph that integrates the similarity between the river network topology and topographic dynamics of a digital elevation model. It introduces dynamic time-warped distance and cross-correlation peak coefficients to achieve adaptive updates of edge weights, forming a dynamic graph structure that evolves with the weather system. A graph-gated recurrent unit with embedded graph convolution operations is used to extract spatiotemporal latent states, and a Transformer decoder is combined to generate multi-step forecast sequences. Furthermore, a hybrid loss function incorporating water balance constraints and geodesic gradient consistency constraints is constructed to embed the physical laws of the watershed into the model training process. This method effectively solves the technical challenges of sparse and non-uniform spatial distribution of meteorological stations in karst landform areas and asynchronous propagation of upstream and downstream floods due to spatiotemporal drift of rainstorm centers. It achieves a deep integration of physical interpretability and data-driven capabilities, significantly improving the accuracy and reliability of watershed meteorological and hydrological forecasts.

[0069] Example 2:

[0070] like Figure 2 As shown, this embodiment provides a watershed meteorological station correlation forecasting system modeled by graph neural networks. See [link to relevant documentation]. Figure 2 The system includes: Acquisition module 701: used to acquire digital elevation model data, multi-source meteorological and hydrological historical data sequences and meteorological radar volume scan data of the target area; to perform depression filling processing and eight-direction flow accumulation calculation on the digital elevation model data to obtain a directed acyclic graph structure containing K sub-basin units, wherein each sub-basin unit includes centroid latitude and longitude coordinates, average elevation and catchment area attributes. Module 702: Used to map the observation values ​​of each monitoring station in the multi-source meteorological and hydrological historical data sequence to the raster grid of the digital elevation model data using bilinear interpolation, and construct the meteorological element background field; after performing ground clutter suppression and range attenuation correction on the meteorological radar volume scan base data, it generates a three-dimensional radar reflectivity factor cube aligned with the raster grid through polar coordinate to rectangular coordinate transformation and vertical linear interpolation; Generation module 703: Based on the meteorological background field and the three-dimensional radar reflectivity factor cube, it extracts local spatial gradient features and vertical profile features for each monitoring station; after standardizing the original observation values, local spatial gradient features, vertical profile features, and station altitude, longitude, and latitude in the multi-source meteorological and hydrological historical data sequence, it performs vector concatenation to generate the initial embedded feature set for each station at all times within the historical time window. Update module 704: Based on the directed acyclic graph structure and the spatial coordinates and elevation of each station, it constructs a physical prior graph containing topological connection edges and terrain dynamic similarity edges; it adaptively updates the edge weights of the physical prior graph through dynamic time warping distance and cross-correlation peak coefficient, and dynamically adds new edges according to the temporal similarity threshold to form a dynamic graph structure that evolves over time. Output module 705: This module takes the initial embedded feature set and the dynamic graph structure as inputs to a stacked graph-gated recurrent unit, outputting the spatiotemporal hidden state sequence and the final encoder hidden state within the entire historical time window. Using the final encoder hidden state as the key matrix and value matrix, and the concatenation of the shared query starting embedding vector and the learnable position encoding vector as the query input to the decoder, this module feeds into a Transformer-based decoder to generate a multi-step forecast sequence containing future rainfall, future temperature, and future outflow from each sub-basin. Based on the multi-step forecast sequence and the corresponding historical observations, a hybrid loss function is constructed by combining the watershed water balance equation and the canyon temperature decay rate constraint along the river topology of the directed acyclic graph structure. The model parameters are iteratively updated until convergence, and the trained graph neural network forecast model is output.

[0071] It should be noted that the specific methods by which each module performs operations in the system described in the above embodiments have been described in detail in the embodiments related to the method, and will not be elaborated here.

[0072] Example 3:

[0073] Corresponding to the above method embodiments, this embodiment also provides a basin meteorological station association forecasting device modeled by graph neural network. The basin meteorological station association forecasting device modeled by graph neural network described below and the basin meteorological station association forecasting method modeled by graph neural network described above can be referred to in correspondence.

[0074] Figure 3 This is a block diagram illustrating a watershed meteorological station correlation forecasting device 800 based on graph neural network modeling, according to an exemplary embodiment. Figure 3 As shown, the basin meteorological station association forecasting device 800 modeled by the graph neural network includes a processor 801 and a memory 802. The basin meteorological station association forecasting device 800 also includes one or more of the following: a multimedia component 803, an I / O interface 804, and a communication component 805.

[0075] The processor 801 controls the overall operation of the graph neural network-modeled watershed meteorological station association forecasting device 800 to complete all or part of the steps in the graph neural network-modeled watershed meteorological station association forecasting method described above. The memory 802 stores various types of data to support the operation of the graph neural network-modeled watershed meteorological station association forecasting device 800. This data may include, for example, instructions for any application or method operating on the graph neural network-modeled watershed meteorological station association forecasting device 800, as well as application-related data such as contact data, sent and received messages, images, audio, video, etc. The memory 802 can be implemented using any type of volatile or non-volatile storage device or a combination thereof, such as Static Random Access Memory (SRAM), Electrically Erasable Programmable Read-Only Memory (EEPROM), Erasable Programmable Read-Only Memory (EPROM), Programmable Read-Only Memory (PROM), Read-Only Memory (ROM), magnetic storage, flash memory, magnetic disk, or optical disk. The multimedia component 803 may include a screen and an audio component. The screen may be, for example, a touchscreen, and the audio component is used to output and / or input audio signals. For example, the audio component may include a microphone for receiving external audio signals. The received audio signals may be further stored in the memory 802 or transmitted via the communication component 805. The audio component also includes at least one speaker for outputting audio signals. I / O interface 804 provides an interface between processor 801 and other interface modules, such as keyboards, mice, or buttons. These buttons can be virtual or physical. Communication component 805 is used for wired or wireless communication between the watershed meteorological station correlation forecasting device 800 (modeled by the graph neural network) and other devices. Wireless communication includes, for example, Wi-Fi, Bluetooth, Near Field Communication (NFC), 2G, 3G, or 4G, or a combination thereof. Therefore, the corresponding communication component 805 may include a Wi-Fi module, a Bluetooth module, or an NFC module.

[0076] In an exemplary embodiment, the basin meteorological station association forecasting device 800 modeled by graph neural network can be implemented by one or more application-specific integrated circuits (ASICs), digital signal processors (DSPs), digital signal processing devices (DSPDs), programmable logic devices (PLDs), field programmable gate arrays (FPGAs), controllers, microcontrollers, microprocessors, or other electronic components to perform the above-described basin meteorological station association forecasting method modeled by graph neural network.

[0077] In another exemplary embodiment, a computer-readable storage medium including program instructions is also provided. When executed by a processor, these program instructions implement the steps of the graph neural network-modeled watershed meteorological station association forecasting method described above. For example, the computer-readable storage medium may be the memory 802 including the program instructions described above. These program instructions may be executed by the processor 801 of the graph neural network-modeled watershed meteorological station association forecasting device 800 to complete the graph neural network-modeled watershed meteorological station association forecasting method described above.

[0078] Example 4:

[0079] Corresponding to the above method embodiments, this embodiment also provides a readable storage medium. The readable storage medium described below can be referred to in conjunction with the graph neural network modeling method for watershed meteorological station association forecasting described above.

[0080] A computer program is stored on a readable storage medium, and when the computer program is executed by a processor, it implements the steps of the watershed meteorological station association forecasting method modeled by graph neural network in the above method embodiment.

[0081] Specifically, the readable storage medium can be a USB flash drive, a portable hard drive, a read-only memory (ROM), a random access memory (RAM), a magnetic disk, or an optical disk, or any other readable storage medium capable of storing program code.

[0082] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A watershed meteorological station correlation forecasting method based on graph neural network modeling, characterized in that, include: Acquire digital elevation model data, multi-source meteorological and hydrological historical data sequences, and meteorological radar volume scan data of the target area; perform depression filling processing and eight-direction flow accumulation calculation on the digital elevation model data to obtain a directed acyclic graph structure containing K sub-basin units, where each sub-basin unit includes centroid latitude and longitude coordinates, average elevation, and catchment area attributes. Bilinear interpolation is used to map the observation values ​​of each monitoring station in the multi-source meteorological and hydrological historical data sequence to the raster grid of the digital elevation model data to construct the meteorological element background field; after performing ground clutter suppression and range attenuation correction on the meteorological radar volume scan base data, a three-dimensional radar reflectivity factor cube aligned with the raster grid is generated through polar coordinate to rectangular coordinate transformation and vertical linear interpolation. Based on the meteorological background field and the three-dimensional radar reflectivity factor cube, local spatial gradient features and vertical profile features are extracted for each monitoring station. The original observations, local spatial gradient features, vertical profile features, and station altitude, longitude, and latitude in the multi-source meteorological and hydrological historical data sequence are standardized and then vectorized to generate the initial embedded feature set for each station at all times within the historical time window. Based on the directed acyclic graph structure and the spatial coordinates and elevation of each station, a physical prior graph containing topological connection edges and terrain dynamic similarity edges is constructed. The edge weights of the physical prior graph are adaptively updated by dynamic time warping distance and cross-correlation peak coefficient, and new edges are dynamically added according to the temporal similarity threshold to form a dynamic graph structure that evolves over time. The initial embedded feature set and dynamic graph structure are input into a stacked graph-gated recurrent unit, which outputs the spatiotemporal hidden state sequence within the entire historical time window and the final encoder hidden state. Using the final encoder hidden state as the key matrix and value matrix, and the concatenation of the shared query starting embedding vector and the learnable position encoding vector as the query input of the decoder, both are fed into the decoder based on the Transformer architecture to generate a multi-step forecast sequence containing future rainfall, future temperature, and future outflow from each sub-basin. Based on the multi-step forecast sequence and the historical observations at the corresponding times, a hybrid loss function is constructed by combining the watershed water balance equation and the canyon temperature decay rate constraint along the river topology of the directed acyclic graph structure. The model parameters are iteratively updated until convergence, and the trained graph neural network forecast model is output.

2. The basin meteorological station correlation forecasting method based on graph neural network modeling according to claim 1, characterized in that, Based on the meteorological background field and the three-dimensional radar reflectivity factor cube, local spatial gradient features and vertical profile features are extracted for each monitoring station. The original observations, local spatial gradient features, vertical profile features, and station altitude, longitude, and latitude from the multi-source meteorological and hydrological historical data sequence are standardized and then concatenated into vectors to generate an initial embedded feature set for each station at all times within the historical time window, including: Based on the meteorological background field, the central difference scheme of 3×3 grid window is used to calculate the temperature meridional gradient, temperature zonal gradient, pressure meridional gradient and pressure zonal gradient of each station to form a four-dimensional local spatial gradient vector. Among them, the temperature meridional gradient reflects the blocking effect of the east-west topographic barrier on heat transport, and the pressure meridional gradient is used to characterize the pressure gradient force driving mechanism of the canyon wind channel. Based on the three-dimensional radar reflectivity factor cube, reflectivity factor values ​​of 20 standard isobaric surface layers above each station are extracted along the vertical direction to form a vertical profile feature vector. Simultaneously, the convective effective potential energy of the station is calculated using a static approximation formula based on the virtual temperature difference, and the vertical wind shear is calculated based on the horizontal wind vector difference between the 200 hPa and 850 hPa isobaric surfaces. Among them, the convective effective potential energy is used to characterize the unstable energy reserve of atmospheric convection, and the vertical wind shear is used to characterize the modulation effect of the vertical wind shear on the organization of convection. The vertical profile features include the reflectivity factor of each isobaric surface layer and the convective effective potential energy and vertical wind shear derived therefrom. The original observation vectors, local spatial gradient vectors, vertical profile feature vectors, convective available potential energy, vertical wind shear, and station elevation, longitude, and latitude from the multi-source meteorological and hydrological historical data sequences are Z-score standardized. After vector concatenation, a heterogeneous node initial feature vector with dimension 35 is generated. The initial feature vectors from all stations at all historical moments constitute the initial embedded feature set. The formula for calculating convective available potential energy is as follows: In the formula, For convective effective potential energy, It is the acceleration due to gravity. To balance the height, For free convection height, The body is weak and cold. For the ambient temperature, Let the integral variable represent altitude. Altitude The differential variable.

3. The basin meteorological station correlation forecasting method based on graph neural network modeling according to claim 1, characterized in that, Based on the directed acyclic graph structure and the spatial coordinates and altitude of each station, a physical prior graph is constructed, including topological connecting edges and terrain-dynamically similar edges, which includes: Based on the directed acyclic graph structure, directed topological connection edges are established for stations corresponding to sub-basins with river confluence relationships. The initial weight is set as the inverse of the river length. The direction of the directed topological connection edge is consistent with the river network confluence direction in the directed acyclic graph structure, which is used to characterize the water contribution relationship between upstream sub-basins and downstream sub-basins. Based on the station's altitude, slope, and normalized distance to the nearest watershed, the similarity of each pair of stations is calculated using the topographic dynamic similarity index formula. The topographic dynamic similarity index comprehensively considers the vertical temperature decrease effect caused by altitude differences, the difference in slope runoff velocity affected by slope differences, and the difference in rainfall runoff response time reflected by the difference in watershed distance. When the similarity index is greater than a preset first threshold, an undirected topographic dynamic similar edge is established, with the initial weight set to the similarity index value; when the similarity index is less than or equal to the preset first threshold, no topographic dynamic similar edge is established. Directed topological connection edges and undirected terrain dynamic similarity edges are merged into a unified edge set. Combined with the set of station nodes, a physical prior graph containing both topological connection edges and terrain dynamic similarity edges is constructed. The formula for calculating the terrain dynamic similarity index is as follows: In the formula, The terrain dynamic similarity index between station p and station q. , The elevations of stations p and q are respectively. , The slopes for stations p and q are respectively. and Let be the normalized distances from stations p and q to the nearest watershed, respectively. , , These represent the statistical standard deviations of altitude, slope, and watershed distance within the watershed.

4. The watershed meteorological station correlation forecasting method based on graph neural network modeling according to claim 1, characterized in that, The method involves adaptively updating the edge weights of the physical prior graph using dynamic time-warped distance and cross-correlation peak coefficient, and dynamically adding new edges based on a temporal similarity threshold to form a dynamic graph structure that evolves over time. This includes: Using the initial embedded feature set, the dynamic time warp distance of each pair of established edges on its historical feature sequence and the peak cross-correlation coefficient within the past 6-hour time window are calculated. The dynamic time warp distance is solved by dynamic programming algorithm to find the optimal time curvature path, which is used to measure the sequence similarity under nonlinear time distortion. The peak cross-correlation coefficient is determined by sliding window cross-correlation analysis to determine the maximum correlation coefficient and its corresponding lag time, which is used to capture the propagation time lag characteristics of the precipitation system. The physical prior weights of the physical prior graph and the data-driven dynamic similarity are weighted and fused using a dynamic weight fusion formula to obtain the final edge weights at each time step. The balance coefficient in the dynamic weight fusion formula is adaptively adjusted according to the characteristics of the flood and drought seasons in the watershed. During the flood season, data-driven similarity is emphasized to capture sudden associations of rainstorms, while during the drought season, physical prior weights are emphasized to maintain a stable topology. For node pairs without physical prior edges, if their dynamic time-warped distance is less than a preset proportion of the historical window length and their cross-correlation peak coefficient is greater than a preset similarity threshold, new edges are dynamically added. The weights are directly obtained from the cross-correlation peak coefficient via a Sigmoid mapping. This mechanism of dynamically adding new edges is used to identify strong correlations among meteorological elements temporarily formed during heavy rain periods, supplementing strong coupling relationships within heavy rain cloud clusters that are not covered by the fixed topology. All edge weights are organized into an adjacency matrix according to time steps, and combined with a fixed node set to generate a dynamic graph structure that evolves over time. The dynamic weight fusion formula is as follows: In the formula, Let be the edge weight between node p and node q at time t. For balance coefficient, For physical prior initial weights, For the Sigmoid function, Scaling factor This represents the peak cross-correlation coefficient between node p and node q over the past 6-hour window. As a regulating factor, Let be the dynamic time-normalized distance between the historical feature sequences of node p and node q.

5. The basin meteorological station correlation forecasting method based on graph neural network modeling according to claim 1, characterized in that, The graph-gated recurrent unit, which inputs the initial embedded feature set and the dynamic graph structure into a stacked array, outputs the spatiotemporal hidden state sequence within the entire historical time window and the final encoder hidden state, including: The feature matrix from the initial embedded feature set and the adjacency matrix with self-loops in the dynamic graph structure are input into the graph convolutional layer. The node features are propagated layer by layer using the degree-normalized convolution formula. After passing through the linear rectified activation function, the node feature matrix of the current graph convolutional layer is output. The adjacency matrix with self-loops is the sum of the edge weight matrix and the identity matrix of the dynamic graph structure, which is used to preserve the node's own feature information. The degree-normalized convolution formula suppresses the feature drowning effect of high-connectivity nodes by symmetrically normalizing the degree matrix, and the elements of the degree diagonal matrix are the sum of the corresponding rows of the node, which is used to adaptively adjust the weight distribution of neighborhood aggregation. The graph convolution operation is embedded in the three types of gating structures of the graph gated recurrent unit: update gate, reset gate, and candidate hidden state. Each gating structure takes the concatenation of the node feature matrix at the current time step and the hidden state matrix at the previous time step as input. After graph convolution, it is activated by the Sigmoid or hyperbolic tangent activation function, and the hidden state matrix at each time step is recursively calculated by the final hidden state update formula. The update gate is used to control the retention ratio of historical hidden states, the reset gate is used to adjust the influence of historical hidden states on candidate hidden states, and the candidate hidden states are used to generate the latent feature representation at the current time step. After traversing all historical time steps, the spatiotemporal hidden state sequence and the final encoder hidden state are output. The spatiotemporal hidden state sequence encodes the comprehensive feature representation of each station after integrating spatial neighborhood association and temporal evolution law.

6. The basin meteorological station correlation forecasting method based on graph neural network modeling according to claim 1, characterized in that, The process uses the final encoder hidden state as the key matrix and value matrix, and the concatenation of the shared query starting embedding vector and the learnable positional encoding vector as the query input to the decoder. These are fed into a Transformer-based decoder to generate a multi-step forecast sequence containing future rainfall, future temperature, and future outflow from each sub-basin for each station. This sequence includes: For each future moment in the decoding process, a learnable position encoding vector is initialized and generated using sine and cosine functions to represent the relative positional relationship of the forecast time and enable the model to perceive the difference in forecast difficulty for different lead times. Its even-numbered and odd-numbered dimensions are generated using sine and cosine functions, respectively, and the dimension index is related to the model dimension. The position encoding vector is concatenated with the shared query start embedding vector as the initial query input of the decoder. An 8-head attention mechanism is adopted to divide the final encoder hidden state into 8 subspaces. Each subspace independently learns different field-station association patterns. By introducing a physical constraint mask in the attention weight calculation, specific attention heads can focus on the confluence direction or elevation similarity. At the same time, the final encoder hidden state is linearly transformed and used as the key matrix and value matrix. The query input is linearly transformed to obtain the query matrix. The attention output is calculated by the scaled dot product attention formula. The scaled dot product attention formula prevents gradient vanishing due to excessively large dot product values ​​by dimensional normalization. The outputs of multiple attention heads are concatenated and then linearly transformed, followed by a feedforward network, residual connections, and layer normalization to obtain the decoder hidden state at each future time step. The feedforward network adopts a two-layer fully connected structure with Gaussian error linear unit activation function embedded in the middle to enhance nonlinear fitting ability. The residual connection is used to alleviate the gradient decay problem of deep networks, and the layer normalization is used to stabilize the distribution drift during the training process. The hidden state of the decoder is input into the feedforward network to generate a forecast vector for each future time step. The output layer of the forecast vector uses a piecewise linear activation function. The rainfall output uses ReLU activation to ensure non-negativity, the temperature output uses linear activation to allow for positive and negative deviations, and the flow output uses Softplus activation to ensure positive values ​​and smooth differentiability. The dimension of the forecast vector is equal to twice the sum of the total number of stations and the total number of sub-basins. The first half of the vector contains the future rainfall and temperature forecasts for each station, and the second half contains the future outflow forecasts for each sub-basin. The forecast vectors for all time steps are arranged in chronological order to output a multi-step forecast sequence.

7. The basin meteorological station correlation forecasting method based on graph neural network modeling according to claim 1, characterized in that, The method involves constructing a hybrid loss function based on multi-step forecast sequences and corresponding historical observations, combined with the watershed water balance equation and the canyon temperature decay rate constraint along the river topology of a directed acyclic graph structure. The model parameters are iteratively updated until convergence, including: Based on the predicted outflow of each sub-basin, the predicted rainfall at each station, and the historical observations at the corresponding time in the multi-step forecast sequence, the forecast error term is calculated using the mean square error formula. The mean square error formula uses differentiated weights for each forecast element. The weight of the rainfall forecast error increases with the extension of the forecast period, and the weight of the flow forecast error automatically increases during the flood peak period. The physical consistency loss term was calculated using the water balance constraint formula based on the Thiessen polygon area, runoff coefficient, and baseflow estimate of each station. The Thiessen polygon area was generated by the Voronoi diagram algorithm based on the spatial distribution of the stations to characterize the effective control range of each station. The runoff coefficient was determined by the lookup table method based on the land use type and soil texture to convert rainfall into effective runoff. The baseflow estimate was obtained by segmenting the measured flow using the receding water curve method to characterize the stable recharge of groundwater to river runoff. Based on the canyon temperature lapse rate along the river topology of the directed acyclic graph structure, the temperature gradient error term is calculated using the geodesic gradient consistency formula. The canyon temperature lapse rate is calculated using the environmental lapse rate formula based on the elevation to constrain the spatial continuity of temperature forecast. The geodesic gradient consistency formula calculates the temperature difference between adjacent sub-basins along the shortest path of the river to penalize abnormal gradients that violate topographic and climatic patterns. After normalizing each loss term, the total loss is obtained by weighted summing of the prediction error term, water balance constraint loss term, and geodesic gradient consistency loss term. The weighting coefficients are dynamically adjusted based on the convergence progress during training. Initially, emphasis is placed on the prediction error term to quickly fit the observed data; later, emphasis is placed on the physical consistency loss term to enhance the satisfaction of physical constraints. An adaptive moment estimation optimizer is used to iteratively update the model parameters until convergence. The formula for calculating the water balance constraint loss term is as follows: In the formula, The loss term is the water balance constraint term, and F is the forecast duration. Indexing future moments For the total number of sub-basins, Let be the predicted outflow from the i-th sub-basin at time τ. M represents the time step, and M represents the total number of stations. Let be the predicted rainfall value for the j-th station at time τ. Let j represent the area of ​​the Thiessen polygon. Let be the runoff coefficient of the area where the j-th station is located. This is the estimated base current at time τ.

8. The watershed meteorological station correlation forecasting method based on graph neural network modeling according to claim 1, characterized in that, The process of filling depressions and calculating eight-directional flow accumulation in the digital elevation model data yields a directed acyclic graph structure containing K sub-basin units, including: Based on the digital elevation model data, depression filling is used to eliminate false depressions. Depression filling is achieved by identifying local minimum elevation points and raising them to overflow elevations to eliminate water traps in the terrain data and ensure the continuity of river network extraction. The flow accumulation amount of each grid cell is calculated using an eight-directional flow accumulation algorithm, with a threshold of 1000 grid cells. The algorithm determines the flow direction based on the steepest slope principle and accumulates the number of grid cells in the upstream contribution area to identify potential river channel locations. The threshold is adaptively adjusted based on the watershed area and resolution: a high threshold is used for large watersheds to extract the main channel, while a low threshold is used for small watersheds to preserve tributary details. Grid cells exceeding the threshold are marked as river channel grid cells, and the river channel grid cells are vectorized to generate a river network, converting grid line features into vector line segments and establishing the topological connection relationship of the river network. The Stellar River Network Classification Algorithm is used to classify the river network. The Stellar River Network Classification Algorithm determines the river level based on the tributary inflow rules. The source river without tributary inflow is classified as Level 1. The level is incremented by one when two rivers of the same level converge. When rivers of different levels converge, the higher level is taken to represent the level of river inflow contribution. The source sub-basin of Level 1 river is used as the smallest unit. Upstream and downstream connections between sub-basins are established through the river network inflow relationship to form a directed acyclic graph structure. Each node corresponds to a sub-basin and includes the centroid latitude and longitude coordinates, average elevation and catchment area of ​​the sub-basin. Each directed edge represents the inflow direction from the upstream sub-basin to the downstream sub-basin.

9. The basin meteorological station correlation forecasting method based on graph neural network modeling according to claim 1, characterized in that, The output of the trained graph neural network prediction model also includes a deployment and execution step, which includes: The trained graph neural network prediction model is solidified in a serialized format that supports dynamic graph structures and loaded into the server's graphics processor memory. Then, a pipelined parallel strategy is adopted to divide data acquisition, feature extraction, graph structure update and model inference into four pipeline stages and bind them to independent CPU cores. The feature extraction pipeline stage runs on a general-purpose CPU to perform normalization and vector concatenation operations, while the model inference pipeline stage runs on a high-performance graphics processor to perform graph convolution and attention calculations. A forecasting task is triggered every 15 minutes using a sliding window approach. Each forecast reads rolling historical data from the past 72 hours and performs forward inference to generate hourly forecasts for the next 24 hours. Then, the outflow from each sub-basin in the multi-step forecast sequence is converted into the inflow from each cascade hydropower station through a pre-calibrated river confluence model. The pre-calibrated river confluence model is based on the linear assumption of the relationship between river storage and outflow and is calibrated using parameter estimation methods to determine the storage constant and flow proportion coefficient. The data is overlaid on the watershed digital twin platform in the form of a geographic information system layer. When the forecasted inflow exceeds the safety threshold corresponding to the flood control limit, an orange warning signal is automatically generated. When the forecasted inflow exceeds the design flood level, a red warning signal is automatically generated. When the forecasted inflow does not exceed the safety threshold corresponding to the flood control limit, the normal forecast status is maintained and monitoring continues. When the forecasted inflow is between the flood control limit and the design flood level, a blue attention prompt is generated to remind dispatchers to strengthen tracking and is pushed to the dispatch duty terminal via 5G message.

10. A basin meteorological station association forecasting system based on graph neural network modeling, wherein the basin meteorological station association forecasting method based on graph neural network modeling as described in claim 1 is characterized in that, include: Acquisition module: used to acquire digital elevation model data, multi-source meteorological and hydrological historical data sequences, and meteorological radar volume scan data of the target area; to perform depression filling processing and eight-direction flow accumulation calculation on the digital elevation model data, and to obtain a directed acyclic graph structure containing K sub-basin units, where each sub-basin unit includes centroid latitude and longitude coordinates, average elevation, and catchment area attributes; The construction module is used to map the observation values ​​of each monitoring station in the multi-source meteorological and hydrological historical data sequence to the raster grid of the digital elevation model data using bilinear interpolation, and construct the background field of meteorological elements; after performing ground clutter suppression and range attenuation correction on the meteorological radar volume scan base data, it generates a three-dimensional radar reflectivity factor cube aligned with the raster grid through polar coordinate to rectangular coordinate transformation and vertical linear interpolation. The generation module is used to extract local spatial gradient features and vertical profile features for each monitoring station based on the meteorological background field and the three-dimensional radar reflectivity factor cube. After standardizing the original observation values, local spatial gradient features, vertical profile features, and station altitude, longitude, and latitude in the multi-source meteorological and hydrological historical data sequence, the module performs vector concatenation to generate the initial embedded feature set for each station at all times within the historical time window. Update module: Based on the directed acyclic graph structure and the spatial coordinates and elevation of each station, it constructs a physical prior graph containing topological connection edges and terrain dynamic similarity edges; it adaptively updates the edge weights of the physical prior graph through dynamic time warping distance and cross-correlation peak coefficient, and dynamically adds new edges according to the temporal similarity threshold to form a dynamic graph structure that evolves over time. Output module: Used to input the initial embedded feature set and dynamic graph structure into the stacked graph gating recurrent unit, and output the spatiotemporal hidden state sequence and the final encoder hidden state within the entire historical time window; Using the final encoder hidden state as the key matrix and value matrix, and the concatenation of the shared query starting embedding vector and the learnable position encoding vector as the query input of the decoder, both are fed into the decoder based on the Transformer architecture to generate a multi-step forecast sequence containing future rainfall, future temperature, and future outflow from each sub-basin. Based on the multi-step forecast sequence and the historical observations at the corresponding times, a hybrid loss function is constructed by combining the watershed water balance equation and the canyon temperature decay rate constraint along the river topology of the directed acyclic graph structure. The model parameters are iteratively updated until convergence, and the trained graph neural network forecast model is output.

Images