Gridding runoff forecasting method based on space-time diagram convolutional network

Abstract

The invention discloses a grid runoff forecasting method based on a space-time diagram convolutional network, and particularly relates to the technical field of runoff forecasting. The method comprises the following steps: firstly, performing grid division on a target drainage basin based on a landform and hydrological response unit, and constructing a space map structure; combining historical rainfall, runoff and remote sensing meteorological data to construct a space-time input feature sequence; extracting a spatial dependency relationship by using a graph convolutional network, extracting a time evolution feature by combining a long-short term memory network, and generating a node representation fused with the spatial-temporal feature; outputting a runoff prediction result of multiple time steps through a full-connection neural network; for a heavy rainfall event, input characteristics and a space structure can be updated in real time, dynamic completion and local graph structure reconstruction of missing data are realized, and finally a gridding runoff distribution graph in a future time period is generated. The method has high precision, high timeliness and high adaptability, and is suitable for runoff intelligent prediction in complex terrains and sudden heavy rainfall scenes.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of runoff prediction, in particular to a meshed runoff prediction method based on a spatio-temporal graph convolution network. BACKGROUND

[0002] With the increasing frequency of extreme weather events and the acceleration of urbanization, the frequency and intensity of flood disasters have increased significantly, which puts higher requirements on urban drainage systems and water resource management. As a key component of the hydrological prediction system, the accuracy of meshed runoff prediction directly affects the real-time and reliability of key decisions such as flood control scheduling, urban drainage, and reservoir regulation.

[0003] Currently, the mainstream runoff prediction method mainly relies on traditional physical models or statistical models. Such methods require a large amount of basic data such as topography, soil parameters, and land use types, have long modeling periods, are sensitive to parameters, and have poor generalization ability. Especially in complex terrain areas or under meteorological mutation conditions, the model performance decreases significantly. At the same time, existing methods often ignore the dynamic coupling relationship of hydrological variables in space and time, making it difficult to fully exploit the spatio-temporal correlation information contained in high-frequency remote sensing and radar rainfall data.

[0004] In urban small watersheds and mountainous rapid-response watersheds, the runoff formation process has the characteristics of "rapid rise, short duration, and peak burst". The hydrological response behavior of such watersheds is highly nonlinear and complex with a long lag, requiring high prediction accuracy. However, due to the sparse distribution of rain gauges and the existence of missing or abnormal observation data, traditional meshed runoff models have large errors and delayed early warnings in such scenarios. SUMMARY

[0005] The purpose of the present application is to provide a meshed runoff prediction method based on a spatio-temporal graph convolution network to solve the problems in the background art.

[0006] To achieve the above purpose, the present application provides the following technical solution: a meshed runoff prediction method based on a spatio-temporal graph convolution network, comprising: S1, obtaining rainfall data, historical runoff observation data, and slope and soil type data of a target watershed, performing grid division on the target watershed, and constructing a spatial graph structure G=(V,E), wherein V represents a set of grid nodes, and E represents the spatial connection relationship between adjacent grids; S2, based on the spatial graph structure, combining historical rainfall-runoff sequence data to construct a spatio-temporal sequence sample data set; S3, using a graph convolution network to extract the spatial feature dependency relationship between nodes and construct a spatial feature map Fst at each time; S4, based on a long short-term memory network to extract the time features in the sequence, modeling Fst in the time dimension to obtain a fused spatio-temporal feature representation DFst; S5, inputting the space-time feature representation DFst into a fully connected neural network, and outputting future multiple time step runoff prediction values corresponding to the grid cells; S6, for the grid cells with missing observations, performing dynamic interpolation completion through an adjacent node feature propagation mechanism in the graph structure; S7, if a sudden heavy rainfall event occurs in the target area, dynamically adjusting the input feature sequence according to real-time radar rainfall data, and real-time updating the spatial graph structure G, and re-executing S3 to S6; S8, outputting a high-resolution gridded runoff distribution map of the entire target basin in a future set period, and predicting the runoff process at the basin level and the city level.

[0007] Preferably, the grid division of the target basin is performed, and a spatial graph structure G=(V,E) is constructed, comprising: extracting the boundary and river network structure of the target basin, combining the terrain slope and the hydrological response unit division principle, and determining the minimum spatial scale of gridding; dividing the target basin into a plurality of regular or irregular grid cells, and assigning a unique node number to each grid cell to construct a node set V; based on the boundary contact relationship and the terrain flow direction of adjacent grid cells, defining a spatial connection rule, generating an edge set E, and forming a spatial graph structure G=(V,E).

[0008] Preferably, based on the spatial graph structure, a time-space sequence sample data set is constructed in combination with historical rainfall-runoff sequence data, comprising: obtaining historical time series rainfall data and measured runoff data corresponding to the node set; the rainfall intensity, evapotranspiration rate, surface temperature and soil moisture of each grid node at the current time and the previous k time steps are used as time step input feature vectors together with their corresponding static attributes; constructing a three-dimensional tensor form of space-time sample structure, wherein the first dimension is the time step sequence, the second dimension is the grid node number, and the third dimension is the multi-dimensional meteorological and geographical features of each node; synchronously constructing the target runoff output sequence corresponding to the time step, ensuring that each sample corresponds to the gridded runoff prediction target of several future time steps.

[0009] Preferably, the spatial feature dependency relationship between nodes is extracted by using the graph convolution network, and a spatial feature graph Fst at each time step is constructed, comprising: slicing the constructed space-time input tensor in the time dimension according to the time step, extracting the node feature matrix corresponding to the current time step as the input of the graph convolution network; The node features are aggregated with the adjacency node features in the spatial graph structure by using a graph convolution operation, and a multiplication operation based on a normalized adjacency matrix and a trainable weight matrix is used to extract the spatial correlation features; A nonlinear activation function is applied to each layer of graph convolution output. The multiple graph convolution layers are repeatedly executed to gradually fuse the multi-level spatial features, and finally obtain the spatial feature map Fst of each node at the current time step.

[0010] Preferably, the long short-term memory network is used to extract the time features in the sequence, and the time dimension of Fst is modeled to obtain the fused spatio-temporal feature representation DFst, including: The spatial feature maps generated at the continuous multiple time steps are arranged in time sequence to construct a three-dimensional input sequence; The three-dimensional input sequence is input into the bidirectional long short-term memory network model, and the short-term dynamic change and long-term dependence features of the node in the time dimension are extracted by forward and backward time sequence learning; The hidden state output by the long short-term memory network is connected and transformed to obtain a representation vector of each node in the current time period. The node output vector after fusing the time features is connected with the corresponding graph convolution result in sequence to obtain the final fused spatio-temporal feature representation DFst of each node.

[0011] Preferably, the spatio-temporal feature representation DFst is input into a fully connected neural network, and the future multiple time step runoff prediction values of the corresponding grid cells are output, including: The fused spatio-temporal feature representation of each node is combined into a feature matrix in the order of node number, and the feature matrix is normalized; The normalized fused spatio-temporal feature matrix is input into a fully connected neural network containing at least two hidden layers, wherein each hidden layer realizes feature mapping through a trainable weight matrix and a bias vector; An output node consistent with the number of target prediction steps is set at the output layer of the fully connected neural network, and the hidden layer output is mapped into the future multiple time step runoff prediction values through linear transformation; The multiple time step runoff prediction values of each node are reconstructed into a gridded runoff prediction sequence in the order of nodes, which is used to generate the gridded runoff spatial distribution result in the future time period.

[0012] Preferably, for the grid cells with observation gaps, the dynamic interpolation completion is performed through the adjacency node feature propagation mechanism in the graph structure, including: A set of nodes with observation gaps is identified, a mask matrix is constructed to identify the missing positions, and the feature values of the missing nodes in the input tensor are set to zero vectors or average feature values. The feature information of adjacent nodes is propagated to the missing node by weighting through the normalized adjacency matrix using the forward propagation capability of the trained graph convolution network; The updated feature vector of the missing node is reconstructed, and the output features of the adjacent nodes are integrated by weighted summation to realize local spatial continuity completion.

[0013] Preferably, when a sudden heavy rainfall event occurs in the target area, the input feature sequence is dynamically adjusted according to real-time radar rainfall data, and the spatial graph structure G is updated in real time, including: Receiving and analyzing real-time radar rainfall data of the target area to generate a dynamic rainfall intensity distribution map; If the radar rainfall intensity of any grid cell exceeds the set threshold, it is marked as a heavy rainfall cell; The input feature sequence of the node corresponding to the marked area is updated by insertion, and the newly added high-frequency rainfall data and the latest soil moisture inversion value are combined to form an extended input sequence; At the same time, the terrain flow direction and connection relationship between the nodes in the heavy rainfall area are recalculated, the spatial graph structure subgraph in the region is dynamically reconstructed and the corresponding subgraph area in the original graph is replaced to form the updated spatial graph structure.

[0014] Preferably, the output is a high-resolution gridded runoff distribution map of the entire target watershed in a future set period, including: Arranging the future multiple time step runoff prediction values of each node in time sequence to generate a runoff prediction matrix with dimensions of node number multiplied by prediction step length; Based on the mapping relationship between the nodes and the original geographic coordinates, the runoff value at each prediction time step is restored to the corresponding grid space position to construct a high-resolution runoff distribution map.

[0015] In the above technical solution, the technical effects and advantages provided by the present application are: 1. The gridded runoff prediction method fusing graph convolution network and long short-term memory network is used to fully mine the high-order correlation features of hydrological variables in spatial topological structure and time dynamic evolution, and significantly improve the modeling accuracy of runoff response process under complex terrain conditions and sudden heavy rainfall events.

[0016] 2. The dynamic input updating mechanism and spatial graph structure subgraph level reconstruction method proposed by the present application support rapid response to regional hydrological state changes triggered by real-time radar rainfall, have strong online adaptability and spatial local sensitivity, and significantly enhance the fine prediction ability of the model for high-risk areas such as urban waterlogging points and mountain torrent confluence areas. BRIEF DESCRIPTION OF DRAWINGS

[0017] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0018] Figure 1 This is a flowchart of the method of the present invention. Detailed Implementation

[0019] 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 embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0020] For examples, please refer to Figure 1 As shown in this embodiment, the gridded runoff forecasting method based on spatiotemporal graph convolutional networks includes: S1. Obtain rainfall data, historical runoff observation data, and slope and soil type data for the target watershed. Divide the target watershed into grids and construct a spatial graph structure G=(V,E), where V represents the set of grid nodes and E represents the spatial connection relationship between adjacent grids.

[0021] First, Digital Elevation Model (DEM) data is used as the basic topographic input data source. The DEM data has a spatial resolution of at least 30 meters and covers the entire geographical area of ​​the target watershed. By utilizing flow direction analysis and catchment area calculation methods, the hydrological boundaries and main river network structure of the watershed can be accurately determined. The specific steps are as follows: The D8 (i.e., unidirectional flow) algorithm is used to calculate the flow direction of each grid cell in the DEM to determine the confluence path of water flow from high to low. Based on the flow direction results, the catchment area of ​​each grid cell is further calculated, which is the total area of ​​the catchment area above the cell. Set a catchment area threshold, such as 1 square kilometer. If the catchment area of ​​a certain grid cell is greater than the threshold, it is identified as a river network node, thereby extracting the main river network structure in the target watershed. Based on the confluence of river networks and the changes in the natural slope of the watershed, the watershed boundary is identified, and vector boundary surface data is generated.

[0022] After completing the above operations, the minimum spatial scale for gridding is determined based on the extracted river network density, slope variation range, and hydrological response unit distribution. Specifically, priority is given to ensuring that each grid can cover a single hydrological response unit and has hydrological independence. The final selected grid size depends on the average distribution length of the area with the highest river network density and the area with the lowest slope variation. For example, in small to medium-sized mountainous watersheds, a 100m × 100m square grid is recommended; while in urban plain watersheds, to improve response accuracy, a finer scale can be selected, such as a 50m × 50m regular grid, or irregular polygonal units can be automatically generated based on the terrain.

[0023] After determining the grid scale, the target watershed is divided into several spatial units according to the above criteria. Each grid unit represents a node in the spatial map structure. Each grid unit must have complete geographic coordinate attributes and be spatially unique. The specific steps for constructing the node set V are as follows: Based on the watershed vector boundary and the raster distribution of DEM data, the entire watershed spatial range is discretized into a set of regular or irregular grids; Each grid cell is assigned a unique number, starting from 1 and incrementing to form a grid identifier ID, which is used for the subsequent attribution and updating of node attributes; Static attribute features of each grid cell are extracted, including rainfall data, historical runoff observation data, and slope and soil type data; the slope is obtained by averaging the elevation gradient within the cell, and the land type data is obtained by overlaying a raster classification map. Define the set of all numbered grid cells as the node set V, i.e., V={v1,v2,...,vn}; where vi represents the i-th grid node numbered i, and n is the total number of grids in the watershed.

[0024] To describe the interactions between grid cells and runoff propagation paths, it is necessary to construct the connections between nodes, forming an edge set E in the spatial graph structure. The construction of edge set E includes: First, determine adjacency relationships. For any two mesh cells, if they share a boundary or vertex, they are defined as adjacent cells, denoted as an adjacency relationship. Using the adjacency matrix construction method, all adjacent node pairs are encoded to form initial connection edge pairs; Secondly, the priority relationship of flow direction is determined. Based on the obtained DEM data and flow direction matrix, the slope aspect and flow path are calculated for any pair of adjacent cells. If the elevation of node v1 is greater than that of node v2, and v2 is in the main flow direction of v1, then the edge is... Setting it as a directed edge indicates a directed edge from node v1 to node v2, which represents the directionality of runoff propagation in space; Next, assign weights to each edge. Weights can be calculated using a weighted combination of factors such as elevation difference, slope coefficient, and land roughness influence factor. For example, if the elevation difference between two adjacent nodes vi and vj is Δh, and the slope angle is θ, then the edge weight wij can be defined as: wij = α·Δh + β·tan(θ) + γ·f(LU), where α, β, and γ are adjustable coefficients, and f(LU) is a land use factor function, which can be set according to different land surface types. Finally, summarize all edge pairs that satisfy the connection relationship and flow direction rules to construct the edge set E: E = (i,j=1,2,…,N); where N is the total number of nodes. This represents a directed edge from node vi to node vj, an edge element in a graph structure used to describe the spatial connection between grid cells or the direction of runoff. The existence of this edge indicates a spatial structural association between nodes vi and vj, or the existence of a hydrologically significant flow path. When there is a flow or adjacency relationship between vi and vj, and the weight wij of the corresponding edge is greater than 0, let... =1; otherwise, let =0.

[0025] After the node set V and edge set E are constructed, the spatial graph structure G is finally formed, which is defined as: G=(V,E), where V represents the set of all grid nodes in the target watershed and E represents the set of directed connections between nodes.

[0026] S2. Based on the spatial graph structure, and combined with historical rainfall-runoff sequence data, construct a spatiotemporal sequence sample dataset.

[0027] First, historical rainfall data and measured runoff data are extracted from the target watershed coverage area. Rainfall data should be in gridded or rasterized form (e.g., longitude-latitude grid), and the rainfall data for each time step is distributed to each grid cell in the node set using bilinear interpolation or inverse distance weighting. If the measured runoff data are observations at the confluence point, they are spatially back-mapped to each node using a watershed response model or geomorphic weighting method to achieve runoff estimation at the grid scale.

[0028] For any node in the node set, dynamic meteorological features of the current time and several previous time steps need to be extracted at each time step, and combined with its static attributes to construct the input feature vector. Setting the historical time window length to k, the input feature vector at each time step t consists of the following parts: Dynamic meteorological characteristics include four categories: rainfall intensity, evapotranspiration rate, surface temperature, and soil moisture at the current time and k times ago. Each category of characteristics is a time series of length k+1. Static geographic features include average elevation, slope, land use type, soil type, land cover index, etc. Each feature is a fixed value or categorized code with a fixed length. The final input feature vector length is 4 times (k+1) plus the number of static features. The input vector is uniformly a one-dimensional floating-point sequence for subsequent model processing.

[0029] Specifically, for time step t and node vn, the input features are denoted as: , where is a real-valued vector of length d, and d is the total length of all the features mentioned above. This input vector will be used as the input to the graph neural network and time series model.

[0030] Based on the aforementioned time-step input feature vectors and node set numbers, a complete spatiotemporal input sample tensor can be constructed. This tensor has a three-dimensional structure, and its dimensions are defined as follows: The first dimension is the time dimension, representing T consecutive time steps. The value of T depends on the model training window setting, and is usually between 100 and 1000. The second dimension is the spatial dimension, representing all grid cells in the node set, denoted as N nodes; The third dimension is the feature dimension, which represents the input feature dimension of each node at each time step, and is denoted as d.

[0031] The final input tensor is denoted as X, with a shape of (T, N, d), meaning that at each time step and at each node, there is a multi-dimensional feature vector containing both dynamic meteorological and static geographical information. This tensor can be directly fed into the spatiotemporal graph convolutional network model to achieve feature learning across time and space.

[0032] All features are normalized so that their value range is uniformly between 0 and 1.

[0033] To achieve supervised training of the runoff prediction model, a corresponding target output sequence needs to be configured for each input sequence while constructing the input sample tensor. The specific method is as follows: The target prediction step size h is selected, which means that the runoff value of each input sequence needs to be predicted for the next h consecutive time steps; For the input sequence at time step t, the measured runoff values ​​from time step t+1 to t+h are extracted and used as the prediction target sequence for the current sample. The target output is a three-dimensional tensor Y with dimensions (T, N, h), which is the predicted runoff change value for each node at each input time step, h steps in the future, in cubic meters per second. By aligning the input tensor X with the output tensor Y and pairing them with samples, a training sample set (X,Y) is formed, which can be used for supervised training of deep learning models.

[0034] S3. Use graph convolutional networks to extract spatial feature dependencies between nodes and construct the spatial feature map Fst at each time step.

[0035] First, the constructed three-dimensional spatiotemporal input tensor needs to be sliced ​​along the time dimension. Let this tensor be X, with dimensions T×N×s, where T represents the number of time steps, N represents the total number of nodes, and s represents the input feature dimension of each node.

[0036] In graph convolution operations, the network can only process two-dimensional feature matrices. Therefore, for any time step t, a time slicing operation is performed to extract the two-dimensional feature matrix Xt of the corresponding node, with a dimension of N×s, representing all the input features of each node at that time step. This matrix will serve as the input to the graph convolutional network at the current time step. This process can be formally described as: In the tensor X, with the time dimension fixed at t, extract the two-dimensional matrix Xt∈ , which serves as the input features for graph convolution operations.

[0037] The core of graph convolution lies in weighted aggregation of each node's own features and the features of its neighboring nodes to extract spatial dependencies. This process relies on the adjacency matrix A ∈ [the graph structure]. , where Aij represents whether there is a connection between node i and node j (1 if a connection exists, 0 otherwise), and this matrix has been constructed in the spatial graph structure G=(V,E). To ensure the stability of feature propagation and the controllability of the numerical range, a normalized adjacency matrix  is used, and its calculation method is as follows: ;in, This means adding the identity matrix I to the adjacency matrix (indicating the introduction of a self-join), i.e. Then calculate the degree matrix D, where Dii is the degree of node i, which is the sum of the elements in the i-th row of A'.

[0038] The formula for the first layer of graph convolution is: ,in The first layer of graph convolution outputs the feature matrix. Let be the trainable weight matrix, σ represent the non-linear activation function, and h be the dimension of the output features of this layer. This operation means that for each node, its output feature is the weighted sum of its own features and the features of all its neighboring nodes, multiplied by the weight matrix, and then subjected to a non-linear transformation.

[0039] This step enables spatial aggregation of node features, mapping the original node attributes to a new spatial feature space to capture the hydrological response correlations inherent in the spatial structure.

[0040] To enhance the model's expressive power and ability to fit complex spatial relationships, a non-linear activation function needs to be introduced after the output of each graph convolution layer. The Rectified Linear Unit (ReLU) function is recommended, and its definition is: For any input x, ReLU(x) = max(0,x).

[0041] In practical applications, the ReLU function remains unchanged for positive values ​​and sets them to zero for negative values, which can effectively alleviate the gradient vanishing problem and enhance the sparsity of the model.

[0042] After performing the ReLU activation operation, the feature matrix becomes: This non-linear operation enhances the ability of graph convolutional layers to fit non-linear spatial relationships and is an indispensable part of building deep graph neural networks.

[0043] A single-layer graph convolution can only perceive the first-order neighbor information of a node. To capture the spatial dependencies between nodes at greater distances, the graph convolution operation needs to be stacked in multiple layers, with the input of each layer being the output of the previous layer, thereby aggregating higher-order neighbor information layer by layer.

[0044] Let the number of graph convolution stacked layers be L, and the weight matrix of each layer be... (l=1,...,L), then the output feature matrix of the l-th layer The calculation method is as follows: ,in .

[0045] After L layers of stacking, the final output spatial feature map is denoted as Fst, i.e., Fst = The spatial feature map Fst has a dimension of N×h', where h' is the output dimension of the last layer of graph convolution. This spatial feature map Fst represents the spatial response features formed by each node after integrating the spatial influences of its neighbors and higher-order neighbors at the current time step t. It will serve as the input basis for subsequent time series modeling (such as recurrent neural networks) to complete the joint learning of spatiotemporal coupling features.

[0046] S4. Based on the Long Short-Term Memory network, extract the temporal features in the sequence, model the temporal dimension of Fst, and obtain the fused spatiotemporal feature representation DFst.

[0047] The spatial feature maps obtained in the previous stage are stacked in chronological order. Let the spatial feature map at each time step t be Fst, with a dimension of N×h, where N is the total number of nodes and h is the feature dimension output by the graph convolutional network. The Fst maps from T consecutive time steps are concatenated to form a 3D input sequence X_LSTM with a dimension of T×N×h. This sequence is represented as: X_LSTM={Fst_1,Fst_2,...,Fst_T}∈ .

[0048] The above input sequence is processed using a bidirectional long short-term memory network. This network structure consists of two parallel LSTM subnetworks, which process the input data in forward and reverse time order, respectively, and output the hidden states in two directions at each time step.

[0049] The 3D input sequence X_LSTM is sliced ​​along the node dimension. For each node vn, its spatial feature sequence at T time steps is extracted: Xn={Fst_1[n],Fst_2[n],...,Fst_T[n]}∈ This sequence serves as the input to the LSTM network, and is processed progressively through time steps to generate the hidden state at each time step. Assuming the hidden state dimension is d, the output of the forward LSTM is... →∈ Output to LSTM ←∈ The output of a bidirectional LSTM can be expressed as: bi = concat( →, ←)∈ , where concat represents the connection operation on the feature dimension, and 2s is the dimension of the hidden state after fusion.

[0050] To construct a temporal representation of each node within the current time period, global features for the entire time window need to be extracted from the bidirectional LSTM output sequence. A common approach is to use average pooling in the temporal dimension to extract the average of the hidden states of each node over T time steps, which serves as its temporal feature vector. For node vn, its temporal feature vector is denoted as: ;in, This represents the fusion hidden state of node n at time step t, where 2s is the time feature dimension.

[0051] Finally, the temporal feature vector Tn of each node is concatenated with the graph convolution output feature Fst_T[n] of its corresponding time step to construct the fused spatiotemporal feature representation DFst_n: DFst_n=concat(Fst_T[n],Tn)∈ Where Fst_T[n] represents the spatial features of node n at the last time step, and Tn represents the temporal features within its corresponding time window. These two features are fused into a single spatiotemporal feature vector through a concatenation operation. The above operation is performed on all nodes to obtain the fused feature matrix DFst∈ , where N is the total number of nodes, and (h+2s) is the fusion feature dimension of each node.

[0052] S5. Input the spatiotemporal feature representation DFst into the fully connected neural network and output the future runoff forecast values ​​for the corresponding grid cells at multiple time points.

[0053] The fused spatiotemporal feature representation of each grid node (denoted as DFstn) is arranged sequentially in ascending order of node number to form the feature matrix DFst∈ , where N represents the total number of nodes, and d' represents the spatiotemporal feature dimension of each node after fusion (usually the sum of the graph convolution output dimension and the temporal feature dimension). To improve the numerical stability and training convergence efficiency of the neural network, the feature matrix needs to be normalized. Using the min-max normalization method, the following operations are performed on each feature column: ;in This represents the original value of the i-th node in the j-th feature dimension. For the minimum original value, For the maximum original value, This is the normalized value. The normalized feature matrix DFst retains the original dimensions, but its numerical range is compressed to 0 to 1.

[0054] The normalized feature matrix DFst∈ The input is fed into a fully connected neural network with a multi-layer structure. This network contains at least two hidden layers, each of which performs an affine transformation from input to output through a trainable weight matrix and bias vector, and enhances feature representation through a non-linear activation function.

[0055] Let the weights of the first hidden layer be W1∈ Bias is Then the output of the first layer is: ; where σ represents the activation function, preferably the rectified linear unit function, which is defined as: σ(x)=max(0,x); Similarly, the weights of the second hidden layer are Bias is Its output is: The number of hidden layers and the number of neurons per layer can be adjusted according to the actual data scale and computing power, but it is necessary to ensure that the network has sufficient nonlinear fitting ability to express complex hydrological response characteristics.

[0056] The output layer has the same number of output nodes as the target prediction step size h, with each node outputting a runoff prediction value for one future time step. Let the output layer weights be... Bias is Then the final output of each node is: Where Rn∈ This represents the sequence of predicted runoff values ​​for node n over the next h time steps. This output is a continuous variable, with units of cubic meters per second, representing the surface runoff intensity within the corresponding grid cell. The network is trained by minimizing the mean squared error loss function between the predicted and measured values. In the formula, The predicted values ​​Rn output by all nodes are taken as measured values. These predicted values ​​are then reconstructed into an N×h matrix R, representing the runoff prediction values ​​for N grid cells over the next h time steps. For each time step t∈{1,2,...,h}, the t-th column of matrix R is extracted to generate a complete gridded runoff prediction layer, corresponding to the spatial runoff distribution at the future time step t.

[0057] S6. For grid cells with observation gaps, dynamic interpolation is performed to fill in the gaps through the feature propagation mechanism of adjacent nodes in the graph structure.

[0058] First, the feature dimensions of each node in the input tensor are scanned at each time step. If any key dynamic feature of a node (such as rainfall intensity, evapotranspiration rate, soil moisture) or all features are missing values ​​(which can be represented as NaN or null values), then the node is marked as an "observation missing node".

[0059] Construct the mask matrix M∈ Where N is the total number of nodes and d is the feature dimension. The element Mi,j in the mask matrix is ​​defined as follows: if the j-th feature of the i-th node has an observation, then Mi,j=1; if it is missing, then Mi,j=0. In subsequent feature propagation operations, this mask matrix is ​​used to identify missing positions and initialize the initial feature vector of the missing nodes, using one of two methods: Zero-vector imputation: Set all feature vectors Xi of the missing nodes to 0; Average feature imputation: The feature of missing nodes is set to the mean of the feature of non-missing nodes in the corresponding feature dimension.

[0060] In the constructed graph structure G=(V,E), the connection relationships between each node are determined by the normalized adjacency matrix. According to the propagation rules of graph convolution, the node feature update formula is as follows: Where: X represents the initial node feature matrix (including missing features); W is the weight of the trained graph convolutional layer; H is the output representation after node feature propagation.

[0061] During forward propagation, for missing nodes, since their initial features are set to zero or the average value, their updated values ​​are entirely derived from neighboring nodes through weighted propagation based on adjacency relationships, reflecting the spatial dependency characteristic of the graph structure. Normalized adjacency matrix. The construction method is as follows: .

[0062] In the graph convolutional forward propagation output H, the updated feature vector of each node contains information from its neighboring nodes. For nodes with missing features, their updated feature vectors... It is actually a weighted linear combination of the features of adjacent nodes.

[0063] This update result can be viewed as the interpolated complete vector for node i. In terms of the missing feature dimension, this output is written into the original feature matrix as the new valid feature of the node: Where ⊙ represents element-wise multiplication, Let be the mask vector for node i. This represents the completed feature vector of the i-th node after graph convolution feature propagation and interpolation. This operation ensures that non-missing values ​​are preserved, and missing values ​​are replaced with the graph convolution output.

[0064] Finally, a complete feature matrix is ​​formed, which can be used in subsequent time series modeling steps to ensure the integrity of the input to the subsequent LSTM network and avoid training or prediction interruptions due to null values.

[0065] S7. If a sudden heavy rainfall event occurs in the target area, dynamically adjust the input feature sequence based on real-time radar rainfall data, update the spatial map structure G in real time, and re-execute S3 to S6.

[0066] First, real-time radar reflectivity data of the target area is received via a meteorological data interface, with a temporal resolution of no more than 5 minutes and a spatial resolution of no less than 1 kilometer. The reflectivity-rainfall intensity relationship is then used for conversion. Where R is the rainfall intensity (unit: millimeters per hour); Z is the radar reflectivity (unit: millimeters to the power of 6 per cubic meter), and commonly used empirical coefficients are a=0.0365, b=0.625.

[0067] The converted rainfall intensity values ​​are mapped in raster form to the corresponding spatial range of the original grid node set, forming a radar rainfall intensity distribution map with the same resolution as the node set. This map is used to identify local high-intensity rainfall areas to support subsequent dynamic response operations.

[0068] Based on historical climate statistics of the target area, a heavy rainfall identification threshold Rth is set, for example, 50 mm / hour in urban watersheds, and can be appropriately reduced to 30 mm / hour in mountainous or small-to-medium watersheds.

[0069] For each node vi, if the current radar rainfall intensity Ri ≥ Rth within its mapping region, then mark it as a "heavy rainfall unit node" and add it to the set. This labeling process is performed in real time, ensuring that the model can quickly identify and respond to sudden extreme rainfall events.

[0070] For all nodes marked as heavy rainfall units, their original input feature sequences consist of dynamic meteorological and static geographical features at a fixed time step. Now, it is necessary to introduce higher-frequency observation features along the time dimension. The specific steps are as follows: Extract the newly added radar rainfall intensity sequence of this node within the last L time steps (e.g., the last 30 minutes, one value every 5 minutes, for a total of 6 data points), and denote it as . ; Simultaneously, the latest soil moisture inversion results (such as those obtained from remote sensing inversion or ground observation) are acquired and populated into each new time step; The newly added data sequence is concatenated chronologically to the end of the corresponding time dimension of the original input tensor to construct the extended input sequence. Where T is the original time step and d is the feature dimension. This extended sequence will serve as the input basis for this node in the current forecast period, ensuring that the model captures the immediate impact of abrupt rainfall changes on the runoff response.

[0071] Since heavy rainfall can significantly alter local runoff paths, it is necessary to recalculate the connectivity between grid nodes within the heavy rainfall area to form a new spatial subgraph structure. The specific steps are as follows: Extract the grid cells corresponding to the set of heavy rainfall nodes Vrain; Calculate the main direction and slope value of each node in the area based on real-time elevation data or DEM data; Determine the potential confluence path of each node under heavy rainfall conditions, reconstruct the adjacency relationship based on the maximum slope principle, construct a new edge set Erain, and then form a local subgraph Grain=(Vrain,Erain); In the original graph structure G=(V,E), the edge set Eold of the corresponding node subset is replaced with Erain to generate the updated graph structure G′=(V,E′), where .

[0072] S8. Output a high-resolution gridded runoff distribution map of the entire target watershed for a future set time period, and forecast runoff processes at the watershed and city levels.

[0073] First, organize the runoff prediction results for each node in the future prediction period. Assuming the total number of nodes is N and the target prediction step size is H (for example, if the time resolution is 1 hour for the next 6 hours, then H=6), then each node corresponds to a prediction vector of length H, representing its surface runoff (in cubic meters per second) at each future time step.

[0074] Arrange the predicted vectors of all nodes in numerical order to form a two-dimensional matrix R with dimensions N rows × H columns: ;in This represents the predicted runoff value for the Nth node at prediction time step H. The runoff value for each node in the prediction matrix is ​​then mapped back to its corresponding geographical location. During the aforementioned graph structure construction phase, a unique spatial coordinate mapping relationship has been established for each node within the target watershed grid, typically represented by the grid center point coordinates or the vector boundary center.

[0075] Let the spatial grid center corresponding to node i be (xi, yi). Then, at any prediction time step t, the following will be performed: Mapping to two-dimensional spatial coordinates (xi, yi). This operation transforms node data into a spatially distributed form.

[0076] The two-dimensional runoff distribution map obtained by down-interpolation at each prediction time step is taken as a layer, and a complete gridded runoff layer sequence is generated in time step order: R={R1(x,y),R2(x,y),…,RH(x,y)}; where each Rt(x,y) is a two-dimensional spatial matrix with a resolution consistent with the input grid (e.g., 100m×100m per grid), and the value represents the surface runoff intensity within the unit of the region.

[0077] This layer sequence can be used to visualize runoff evolution, showcasing dynamic characteristics such as runoff diffusion, convergence, and peak shift through consecutive image frames. Administrative boundaries, river networks, and infrastructure layers can be further overlaid to achieve geospatial correlation.

[0078] To meet the needs of urban flood control, reservoir scheduling, and emergency response, it is necessary to extract runoff process information for key areas from the above runoff layer sequence.

[0079] First, define the set of key monitoring units. It can be determined in the following way: For urban areas: Select corresponding nodes in areas such as historical waterlogging points, the vicinity of drainage pumping stations, and underpasses; For mountainous watersheds: select the outlet node of the catchment area or the node with the lowest terrain; For any monitoring node vs∈S, extract its predicted sequence in the R matrix: Qs={rs,1,rs,2,…,rs,H}; plot the curve of Qs over time to obtain key indicators such as peak time, peak size, runoff duration, and rising and falling rates, which can be used to achieve early perception of the development trend of local flooding or the evolution of flood peaks.

[0080] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A gridded runoff forecasting method based on spatiotemporal graph convolutional networks, characterized in that: include: S1. Obtain rainfall data, historical runoff observation data, and slope and soil type data for the target watershed. Divide the target watershed into grids and construct a spatial graph structure G=(V,E), where V represents the set of grid nodes and E represents the spatial connection relationship between adjacent grids. S2. Based on the aforementioned spatial map structure, and combined with historical rainfall-runoff sequence data, construct a spatiotemporal sequence sample dataset; S3. Use graph convolutional networks to extract spatial feature dependencies between nodes and construct the spatial feature map Fst at each time step; S4. Extract the temporal features from the sequence based on the Long Short-Term Memory Network, model the temporal dimension of Fst, and obtain the fused spatiotemporal feature representation DFst. S5. Input the spatiotemporal feature representation DFst into a fully connected neural network and output the future runoff forecast values ​​for the corresponding grid cells at multiple time points. S6. For grid cells with observation gaps, dynamic interpolation is performed to fill in the gaps through the feature propagation mechanism of adjacent nodes in the graph structure. S7. If a sudden heavy rainfall event occurs in the target area, dynamically adjust the input feature sequence based on real-time radar rainfall data, update the spatial map structure G in real time, and re-execute S3 to S6. S8. Output a high-resolution gridded runoff distribution map of the entire target watershed for a future set time period, and forecast runoff processes at the watershed and city levels.

2. The gridded runoff forecasting method based on spatiotemporal graph convolutional networks according to claim 1, characterized in that: The step of dividing the target watershed into grids and constructing a spatial graph structure G=(V,E) includes: Extract the boundaries and river network structure of the target watershed, and determine the minimum spatial scale of the grid by combining the principles of topographic slope and hydrological response unit division; The target watershed is divided into several regular or irregular grid units, and a unique node number is assigned to each grid unit to construct a node set V; Based on the boundary contact relationship between adjacent grid cells and the terrain flow direction, spatial connection rules are defined to generate an edge set E, forming a spatial graph structure G=(V,E).

3. The gridded runoff forecasting method based on spatiotemporal graph convolutional networks according to claim 1, characterized in that: Based on the aforementioned spatial map structure, and combined with historical rainfall-runoff sequence data, a spatiotemporal sequence sample dataset is constructed, including: Obtain historical time-series rainfall data and measured runoff data corresponding to the node set; The rainfall intensity, evapotranspiration rate, surface temperature, and soil moisture of each grid node at the current time and the k previous times, together with their corresponding static attributes, are used as the time step input feature vector. Construct a spatiotemporal sample structure in the form of a three-dimensional tensor, where the first dimension is the time step sequence, the second dimension is the grid node number, and the third dimension is the multidimensional meteorological and geographical features of each node. Simultaneously construct the target runoff output sequence for the corresponding time step to ensure that each sample corresponds to the gridded runoff prediction target for several future time steps.

4. The gridded runoff forecasting method based on spatiotemporal graph convolutional networks according to claim 1, characterized in that: The step of using a graph convolutional network to extract spatial feature dependencies between nodes and constructing a spatial feature map Fst at each time step includes: The constructed spatiotemporal input tensor is sliced ​​along the time dimension by time step, and the node feature matrix corresponding to the current time step is extracted and used as the input of the graph convolutional network. The node features and their neighboring node features are aggregated in the spatial graph structure using graph convolution operations. Spatial correlation features are extracted by multiplication operations based on the normalized adjacency matrix and the trainable weight matrix. A non-linear activation function is applied to the output of each graph convolution layer; By repeatedly executing multiple graph convolutional layers, multi-level spatial features are gradually fused together, ultimately yielding the spatial feature map Fst of each node at the current time step.

5. The gridded runoff forecasting method based on spatiotemporal graph convolutional networks according to claim 1, characterized in that: The method of extracting temporal features from the sequence based on a Long Short-Term Memory network, modeling the temporal dimension of Fst, and obtaining the fused spatiotemporal feature representation DFst includes: The spatial feature maps generated at multiple consecutive time steps are arranged in chronological order to construct a three-dimensional input sequence. The three-dimensional input sequence is fed into a bidirectional long short-term memory network model. Through forward and backward time series learning, the short-term dynamic changes and long-term dependency features of nodes in the time dimension are extracted simultaneously. The hidden states output by the Long Short-Term Memory network are connected and transformed to obtain the representation vector of each node in the current time period, which incorporates the temporal features. The node output vector after fusing temporal features is concatenated with the corresponding convolution result of the preceding graph to obtain the final fused spatiotemporal feature representation DFst for each node.

6. The gridded runoff forecasting method based on spatiotemporal graph convolutional networks according to claim 1, characterized in that: The process of inputting the spatiotemporal feature representation DFst into a fully connected neural network and outputting future runoff forecast values ​​for corresponding grid cells includes: The obtained spatiotemporal feature representations of each node are combined into a feature matrix according to the node numbering order, and the feature matrix is ​​normalized. The normalized fused spatiotemporal feature matrix is ​​input into a fully connected neural network containing at least two hidden layers, where each hidden layer achieves feature mapping through a trainable weight matrix and a bias vector. In the output layer of the fully connected neural network, the number of output nodes is set to be the same as the number of target prediction steps. The output of the hidden layer is then mapped to the runoff forecast values ​​for multiple future time steps through a linear transformation. The multi-time runoff forecast values ​​of each node are reconstructed into a gridded runoff prediction sequence according to the node order, which is used to generate gridded runoff spatial distribution results for future periods.

7. The gridded runoff forecasting method based on spatiotemporal graph convolutional networks according to claim 1, characterized in that: The method for dynamically interpolating and filling in observation gaps in grid cells using the adjacent node feature propagation mechanism in the graph structure includes: Identify the set of nodes with observation gaps, construct a mask matrix to identify the missing locations, and set the eigenvalues ​​of the missing nodes in the input tensor to the zero vector or the average eigenvalue; By leveraging the forward propagation capability of a trained graph convolutional network, the feature information of adjacent nodes is propagated to the missing nodes through a weighted normalized adjacency matrix. The updated feature vectors of missing nodes are reconstructed, and the output features of adjacent nodes are integrated by weighted summation to achieve local spatial continuity completion.

8. The gridded runoff forecasting method based on spatiotemporal graph convolutional networks according to claim 1, characterized in that: If a sudden heavy rainfall event occurs in the target area, the input feature sequence is dynamically adjusted based on real-time radar rainfall data, and the spatial map structure G is updated in real time, including: Receive and analyze real-time radar rainfall data of the target area to generate a dynamic rainfall intensity distribution map; If the radar rainfall intensity of any grid cell exceeds the set threshold by comparing it with historical statistical thresholds, it is marked as a heavy rainfall cell. The input feature sequence of the corresponding node in the labeled area is inserted and updated, and the newly added high-frequency rainfall data is combined with the latest soil moisture inversion value to form an extended input sequence; Simultaneously, the terrain flow direction and connection relationship between nodes in the heavy rainfall area are recalculated, the spatial graph structure subgraphs within the area are dynamically reconstructed, and the corresponding subgraph regions in the original graph are replaced to form an updated spatial graph structure.

9. The gridded runoff forecasting method based on spatiotemporal graph convolutional networks according to claim 1, characterized in that: The output includes a high-resolution gridded runoff distribution map of the entire target watershed over a future set time period, including: Arrange the future runoff forecast values ​​for multiple time steps corresponding to each node in chronological order to generate a runoff prediction matrix with the dimension of the number of nodes multiplied by the prediction step length. Based on the mapping relationship between nodes and original geographic coordinates, the runoff value at each predicted time step is restored to the corresponding grid spatial location to construct a high-resolution runoff distribution map.

Images