Grid-based runoff forecasting method based on spatio-temporal graph convolution network

Abstract

The application discloses a gridding runoff prediction method based on a space-time graph convolution network, and particularly relates to the technical field of runoff prediction; the application firstly performs grid division on a target basin based on topography and hydrological response units, and constructs a space graph structure; combines historical rainfall, runoff and remote sensing meteorological data to construct a space-time input feature sequence; uses a graph convolution network to extract spatial dependence, combines a long short-term memory network to extract time evolution features, and generates node representation fused with space-time features; outputs runoff prediction results at multiple time steps through a fully connected neural network; for a heavy rainfall event, input features and a space structure can be updated in real time, dynamic completion of missing data and local graph structure reconstruction are realized, and finally a gridding runoff distribution map in a future period is generated; the method has high precision, high timeliness and high adaptability, and is suitable for intelligent runoff prediction under complex terrain and sudden heavy rainfall scenes.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of runoff prediction, in particular to a gridded 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 gridded 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 gridded runoff models have large errors and delayed early warnings in such scenarios. SUMMARY

[0005] The purpose of the present application is to provide a gridded 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 gridded runoff prediction method based on a spatio-temporal graph convolution network, comprising:

[0007] 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 a spatial connection relationship between adjacent grids;

[0008] S2, based on the spatial graph structure, combining historical rainfall-runoff sequence data to construct a spatio-temporal sequence sample data set;

[0009] S3, using a graph convolution network to extract spatial feature dependency relationships between nodes and construct a spatial feature map Fst at each time;

[0010] S4, time characteristics in the sequence are extracted based on a long short-term memory network, time dimensions of the Fst are modeled, and a spatiotemporal feature representation DFst is obtained;

[0011] S5, the spatiotemporal feature representation DFst is input into a fully connected neural network, and future multiple time flow prediction values of the corresponding grid unit are output;

[0012] S6, for the grid unit with an observation gap, a dynamic interpolation completion is performed through an adjacent node feature propagation mechanism in the graph structure;

[0013] S7, if a sudden heavy rainfall event occurs in the target area, the input feature sequence is dynamically adjusted according to real-time radar rainfall data, the spatial graph structure G is updated in real time, and S3 to S6 are re-executed;

[0014] S8, a high-resolution gridded runoff distribution map of the entire target basin in a future set period is output, and a basin-level and city-level runoff process is predicted.

[0015] Preferably, the target basin is divided into grids, and a spatial graph structure G=(V,E) is constructed, comprising:

[0016] The boundary and river network structure of the target basin are extracted, and the minimum spatial scale of gridding is determined in combination with the terrain slope and the division principle of the hydrological response unit;

[0017] The target basin is divided into a plurality of regular or irregular grid units, and each grid unit is assigned a unique node number to construct a node set V;

[0018] Based on the boundary contact relationship and the terrain flow direction of adjacent grid units, a spatial connection rule is defined, an edge set E is generated, and a spatial graph structure G=(V,E) is formed.

[0019] Preferably, based on the spatial graph structure, a spatiotemporal sequence sample data set is constructed in combination with historical rainfall-runoff sequence data, comprising:

[0020] Historical time series rainfall data and measured runoff data corresponding to the node set are obtained;

[0021] The rainfall intensity, evapotranspiration rate, surface temperature and soil moisture of each grid node at the current time and the previous k times are collectively used as a time step input feature vector together with the corresponding static attributes;

[0022] A spatiotemporal sample structure in the form of a three-dimensional tensor is constructed, wherein the first dimension is a time step sequence, the second dimension is a grid node number, and the third dimension is a multi-dimensional meteorological and geographical feature of each node;

[0023] Synchronously build a target runoff output sequence corresponding to the time step, ensuring that each sample corresponds to a grid runoff prediction target for a number of future time steps.

[0024] Preferably, the spatial feature inter-node dependence relationship is extracted by using the graph convolution network to build a spatial feature graph Fst at each time, including:

[0025] The constructed spatio-temporal input tensor is sliced in the time dimension according to the time step to extract the node feature matrix corresponding to the current time step as the input of the graph convolution network;

[0026] The node features and their adjacent node features are aggregated in the spatial graph structure by using the graph convolution operation, and the multiplication operation based on the normalized adjacency matrix and the trainable weight matrix is used to extract the spatial correlation features;

[0027] A nonlinear activation function is applied to each layer of graph convolution output;

[0028] Multiple graph convolution layers are repeatedly executed to gradually fuse multi-level spatial features, and finally obtain the spatial feature graph Fst of each node at the current time step.

[0029] Preferably, the time features in the sequence are extracted based on the long short-term memory network to model the time dimension of Fst to obtain a fused spatio-temporal feature representation DFst, including:

[0030] The spatial feature graphs generated at consecutive multiple time steps are arranged in time sequence to build a three-dimensional input sequence;

[0031] The three-dimensional input sequence is input into a bidirectional long short-term memory network model to simultaneously extract short-term dynamic changes and long-term dependence features of nodes in the time dimension through forward and backward time sequence learning;

[0032] 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 that fuses the time features;

[0033] The node output vector after fusing the time features is connected with the corresponding graph convolution result of the previous sequence to obtain the final fused spatio-temporal feature representation DFst of each node.

[0034] Preferably, the spatio-temporal feature representation DFst is input into a fully connected neural network to output runoff prediction values of a plurality of future time steps corresponding to the grid cells, including:

[0035] The obtained fused spatio-temporal feature representation of each node is combined into a feature matrix in the order of node numbers, and the feature matrix is normalized;

[0036] The normalized fused spatio-temporal feature matrix is input into a fully connected neural network comprising at least two hidden layers, wherein each hidden layer realizes feature mapping through a trainable weight matrix and a bias vector;

[0037] 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 to runoff prediction values at multiple future time steps through linear transformation;

[0038] The runoff prediction values at multiple time steps of each node are reconstructed into a gridded runoff prediction sequence in node order, which is used to generate a gridded runoff spatial distribution result in a future period.

[0039] Preferably, for the grid cells with missing observations, a dynamic interpolation completion is performed through an adjacent node feature propagation mechanism in the graph structure, including:

[0040] A set of nodes with missing observations 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;

[0041] The feature information of adjacent nodes is propagated to the missing nodes through a normalized adjacency matrix weighting by using the forward propagation capability of the trained graph convolutional network;

[0042] The updated feature vectors of the missing nodes are reconstructed, and the output features of the adjacent nodes are integrated through weighted summation to realize local spatial continuity completion.

[0043] Preferably, if 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:

[0044] Real-time radar rainfall data of the target area is received and analyzed to generate a dynamic rainfall intensity distribution map;

[0045] If the radar rainfall intensity of any grid cell exceeds a set threshold, the grid cell is marked as a heavy rainfall cell by comparing with a historical statistical threshold;

[0046] The input feature sequence of the nodes 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;

[0047] Meanwhile, the terrain flow direction and connection relationship between the nodes in the heavy rainfall area are recalculated, the spatial graph structure subgraph in the area is dynamically reconstructed and the corresponding subgraph area in the original graph is replaced to form an updated spatial graph structure.

[0048] Preferably, the output is a high-resolution gridded runoff distribution map of the entire target basin in a future set period, including:

[0049] arranging the runoff prediction values of each node corresponding to a plurality of future time steps in time sequence to generate a runoff prediction matrix with a dimension of the number of nodes multiplied by the prediction step length;

[0050] 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.

[0051] In the above technical solution, the technical effects and advantages provided by the present application are:

[0052] 1. The grid runoff prediction method fusing graph convolution network and long short-term memory network provided by the present application fully excavates the high-order correlation characteristics of hydrological variables in spatial topological structure and time dynamic evolution, significantly improving the modeling accuracy of runoff response process under complex topographic conditions and sudden heavy rainfall events.

[0053] 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

[0054] In order to more clearly illustrate the technical solutions in the embodiments of the present application or the prior art, the drawings needed in the embodiments will be briefly introduced. Obviously, the drawings described below are only some embodiments described in the present application, and other drawings can also be obtained by those skilled in the art based on these drawings.

[0055] Figure 1 The method flowchart of the present application. DETAILED DESCRIPTION

[0056] In order to make the purpose, technical scheme and advantages of the embodiments of the present application clearer, the technical scheme in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are part of the embodiments of the present application, not all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the present application.

[0057] Embodiment, please refer to Figure 1 As shown in the figure, the grid runoff prediction method based on spatio-temporal graph convolution network described in the present embodiment comprises:

[0058] S1, acquire rainfall data, historical runoff observation data, and slope, soil type data of a target watershed, divide the target watershed into grids, and construct a spatial graph structure G=(V, E), wherein V represents a set of grid nodes, and E represents a spatial connection relationship between adjacent grids.

[0059] First, a digital elevation model (DEM) is used as the basic terrain input data source. The DEM data has a spatial resolution of no less than 30 meters and covers the complete geographical range of the target watershed. By using flow direction analysis and catchment area calculation methods, the hydrological boundary and main river network structure of the watershed can be accurately determined. The specific steps are as follows:

[0060] The D8 (single flow direction) algorithm is used to calculate the flow direction of each grid cell from the DEM to determine the flow path from high to low;

[0061] Based on the flow direction results, the catchment area of each grid cell is further calculated, i.e., the total area of the flow accumulation region above the cell;

[0062] Set a catchment area threshold, for example, 1 square kilometer, if the catchment area of a 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;

[0063] According to the river network flow location and the change of the natural slope of the watershed, the watershed boundary is identified, and the vector boundary surface data is generated.

[0064] After the above operations, according to the extracted river network density, slope change amplitude, and hydrological response unit distribution, the minimum spatial scale of the grid processing is determined comprehensively. Specifically, priority is given to ensuring that each grid can cover a single hydrological response unit and have hydrological independence. The final selected grid size depends on the average distribution length of the maximum river network density area and the minimum slope change area. For example, in small and medium-sized mountainous watersheds, a 100m x 100m square grid is recommended; while in urban plain watersheds, to improve response accuracy, a finer scale such as a 50m x 50m regular grid or an irregular polygonal unit generated automatically according to the terrain can be selected.

[0065] After determining the grid scale, the target watershed is divided into several spatial units according to the above standard, and each grid unit represents a node in the spatial graph structure. Each grid unit needs to have complete geographic coordinate attributes and can be uniquely located in space. To construct the node set V, the specific operations are as follows:

[0066] According to the vector boundary of the watershed and the grid distribution of the DEM data, the entire watershed spatial range is discretized into a regular or irregular grid set;

[0067] A unique number is assigned to each grid cell, starting from 1 and increasing, forming a grid identification ID, which is used for subsequent node attribute attribution and update;

[0068] Extract the static attribute features of each grid cell, including rainfall data, historical runoff observation data, and slope, soil type data; where the slope is obtained by averaging the elevation gradient within the cell, and the land type data is obtained by superimposing the raster classification map;

[0069] All numbered grid cell sets are defined as node set V, that is, 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.

[0070] To describe the mutual influence and runoff propagation path between grid cells, the connection relationship between nodes needs to be constructed to form the edge set E in the spatial graph structure. The construction of edge set E includes:

[0071] First, determine the adjacent relationship. For any two grid cells, if they share a boundary or vertex, they are defined as adjacent cells, denoted as adjacency relationship. Through the construction of the adjacency matrix, all adjacent node pairs are coded to form the initial connection edge pair;

[0072] Second, determine the flow direction priority relationship. Based on the obtained DEM data and flow direction matrix, for any pair of adjacent cells, calculate their slope and water flow path. 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 set as a directed edge, indicating the direction of the directed edge from node v1 to node v2, that is, the directionality of runoff transmission in space;

[0073] Third, assign a weight to each edge. The weight can be calculated by combining factors such as elevation difference, slope coefficient, and land roughness influence factor. For example, let the elevation difference between two adjacent nodes vi and vj be Δh, and the slope angle be θ, then the edge weight wij can be defined as: wij=α·Δh+β·tan(θ)+γ·f(LU), where α, β, γ are adjustable coefficients, and f(LU) is a land use factor function that can set a resistance coefficient according to different land surface types. Finally, all edge pairs that meet the connection relationship and flow direction rules are summarized to construct the edge set E: E= (i,j=1,2,…,N); where N is the total number of nodes, denotes a directed edge from node vi to node vj, is an edge element in the graph structure, used to describe the spatial connection or runoff flow direction between grid cells. The existence of this edge represents the spatial correlation between node vi and node vj, or the existence of a hydrological flow path. When there is a flow or adjacency relationship between vi and vj, and the weight of the corresponding edge wij>0, let =1; otherwise, let =0.

[0074] After the construction of the node set V and the edge set E is completed, the spatial graph structure G is finally formed, which is defined as: G=(V,E), where V represents the set of all grid nodes of the target watershed, and E represents the set of directed connection relationships between nodes.

[0075] S2, based on the spatial graph structure, combined with historical rainfall-runoff sequence data, construct a spatio-temporal sequence sample data set.

[0076] First, extract the historical rainfall data and measured runoff data within the coverage of the target watershed. Rainfall data should be in the form of grid or raster (such as longitude-latitude grid), and each time step of rainfall data should be distributed to each grid cell in the node set through bilinear interpolation or inverse distance weighting method. If the measured runoff data is the observation value of the confluence point, then through the watershed response model or the geomorphology weighting method, it is reversely mapped to each node in space to realize the runoff estimation at the grid scale.

[0077] For any node in the node set, the dynamic meteorological features at the current time and the previous k time steps need to be extracted at each time step, and its static attributes are combined to construct the input feature vector. Set the length of the historical time window as k, then the input feature vector at each time step t is composed of the following parts:

[0078] Dynamic meteorological features: including current and previous k time steps of rainfall intensity, evapotranspiration rate, surface temperature and soil moisture, etc. Each type of feature is a time series with length k+1;

[0079] Static geographic features: including average elevation, slope, land use type, soil type, land cover index, etc. Each type of feature is a fixed value or a type-coded, with fixed length;

[0080] The final input feature vector length is 4 times (k+1) plus the number of static features, and the input vector is unified as a one-dimensional float number sequence for subsequent model processing.

[0081] Specifically, for time step t and node vn, its input feature is denoted as is a real value vector with length d, where d is the total length of all features mentioned above. This input vector will be used as the input of the graph neural network and time series model.

[0082] Based on the time step input feature vector and the node set number mentioned above, a complete spatio-temporal input sample tensor can be constructed. This tensor has a three-dimensional structure, and its dimensions are defined as follows:

[0083] The first dimension is the time dimension, representing T consecutive time steps, and the value of T is determined by the model training window setting, usually between 100 and 1000;

[0084] The second dimension is the spatial dimension, representing all grid cells in the node set, denoted as N nodes;

[0085] The third dimension is the feature dimension, representing the input feature dimension of each node at each time step, which is d.

[0086] The final input tensor is denoted as X, with a shape of (T, N, d), i.e. for each time step, each node corresponds to a multi-dimensional feature vector containing meteorological dynamic and geographical static information. This tensor can be directly input into the spatio-temporal graph convolutional network model to realize cross-time and cross-space feature learning.

[0087] All features are normalized to have a value range of 0 to 1.

[0088] In order to realize the supervised training of the runoff prediction model, for each input sequence, the corresponding target output sequence needs to be configured at the same time as the input sample tensor is constructed. The specific method is as follows:

[0089] Select the target prediction step h, which represents the future continuous h time steps of runoff values that need to be predicted for each input sequence;

[0090] For the input sequence at time step t, the runoff measured values from time step t+1 to t+h are intercepted as the prediction target sequence corresponding to the current sample;

[0091] The target output is a three-dimensional tensor Y with dimensions (T, N, h), i.e. for each input time step, for each node, the runoff change value in the future h steps is predicted, with the unit of cubic meters per second;

[0092] Align the input tensor X and the output tensor Y and pair the samples to form the training sample set (X, Y), which can be used for supervised training of the deep learning model.

[0093] S3, use the graph convolutional network to extract the spatial feature dependency relationship between nodes and construct the spatial feature map Fst at each time step.

[0094] 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.

[0095] 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:

[0096] 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.

[0097] 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'.

[0098] 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.

[0099] 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.

[0100] 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:

[0101] For any input x, ReLU(x) = max(0, x).

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

[0103] After performing the ReLU activation operation, the feature matrix becomes: This nonlinear operation improves the fitting ability of the graph convolution layer for nonlinear spatial relationships and is an indispensable part of building deep graph neural networks.

[0104] Single-layer graph convolution can only perceive first-order neighbor information. To capture spatial dependencies between nodes at a longer distance, multiple layers of graph convolution operations need to be stacked, with the input of each layer being the output of the previous layer, thereby aggregating high-order neighbor information layer by layer.

[0105] Let the number of graph convolution stacking layers be L, and the weight matrix of each layer be (l = 1,..., L), then the calculation method of the output feature matrix of the l-th layer is: , where .

[0106] After L layers of stacking, the final output spatial feature map is denoted as Fst, i.e., Fst= , with a dimension of N x h', where h' is the output dimension of the last layer of graph convolution. This spatial feature map Fst represents the spatial response characteristics formed by each node after considering the spatial influence of its neighbors and high-order neighbors at the current time step t, and will serve as the input basis for subsequent time series modeling (such as recurrent neural networks) to complete the joint learning of spatio-temporal coupling features.

[0107] S4, based on the long short-term memory network, extract the time features in the sequence, model the time dimension of Fst, and obtain the fused spatio-temporal feature representation DFst.

[0108] Stack the spatial feature maps obtained at each time step in the previous stage in chronological order. Let the spatial feature map at each time step t be Fst, with a dimension of N x h, where N is the total number of nodes and h is the feature dimension output by the graph convolution network. Concatenate Fst at consecutive T time steps to form a three-dimensional input sequence X_LSTM with a dimension of T x N x h. The sequence representation is:

[0109] X_LSTM = {Fst_1, Fst_2,..., Fst_T} ∈ .

[0110] The input sequence is processed by a bidirectional long short-term memory network. The network structure is composed of two parallel LSTM sub-networks, which process the input data in the forward and reverse directions along time, and output hidden states in both directions at each time step.

[0111] The three-dimensional input sequence X_LSTM is sliced in the node dimension, and for each node vn, a spatial feature sequence Xn = {Fst_1[n], Fst_2[n],..., Fst_T[n]} ∈ is extracted at T time steps. The sequence is input into the LSTM network and processed through time steps to generate hidden states at each time step. Assuming the hidden state dimension is d, the forward LSTM output is → ∈ , and the backward LSTM output is ← ∈ The output of the bidirectional LSTM can be represented as: bi = concat( →, ←) ∈ where concat represents the concatenation operation in the feature dimension, and 2s is the dimension of the fused hidden state.

[0112] To construct the time dimension representation of each node in the current time period, the global features of the entire time window are extracted from the bidirectional LSTM output sequence. The commonly used method is average pooling in the time dimension, which extracts the average value of the hidden state of each node at T time steps as its time feature vector. For node vn, its time feature vector is denoted as: ; where represents the fused hidden state of node n at time step t, and 2s is the time feature dimension.

[0113] Finally, the time feature vector Tn of each node is concatenated with its corresponding graph convolution output feature Fst_T[n] to construct the fused spatio-temporal feature representation DFst_n: DFst_n = concat(Fst_T[n], Tn) ∈ ; where Fst_T[n] represents the spatial feature of node n at the last time step, and Tn represents the time feature in its corresponding time window, which are fused into a single spatio-temporal feature vector through the 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 fused feature dimension of each node.

[0114] S5, input the spatio-temporal feature representation DFst into a fully connected neural network to output the future multiple time step runoff prediction values of the corresponding grid cell.

[0115] The fusion spatio-temporal features of each grid node (denoted as DFstn) are arranged in ascending order of node number to form a feature matrix DFst∈ , where N represents the total number of nodes, and d' represents the dimension of the fusion spatio-temporal features of each node (usually the sum of the output dimension of the graph convolution and the dimension of the time feature). To improve the numerical stability and training convergence efficiency of the neural network, the feature matrix needs to be normalized. The minimum-maximum normalization method is used to perform the following operation on each column of features: ; where represents the original value of the i-th node in the j-th feature dimension, is the minimum original value, is the maximum original value, is the normalized value. The normalized feature matrix DFst maintains the original dimension and the numerical range is compressed to 0 to 1.

[0116] The normalized feature matrix DFst∈ is input into a fully connected neural network with a multi-layer structure. The network contains at least two hidden layers, each implementing an affine transformation from input to output through a trainable weight matrix and bias vector, and enhancing feature expression ability through a nonlinear activation function.

[0117] Let the first hidden layer weight be W1∈ , and the bias be , then the first layer output is: ; where σ represents the activation function, preferably the rectified linear unit function, defined as: σ(x) = max(0, x);

[0118] Similarly, the second hidden layer weight is , the bias is , and the output is: ; The number of hidden layers and the number of neurons in each layer can be adjusted according to the actual data size and computing power, but it needs to ensure that the network has sufficient nonlinear fitting ability to express complex hydrological response characteristics.

[0119] The output layer is set with an output node corresponding to the target prediction step h, and each node outputs the runoff prediction value of the corresponding future time step. Let the output layer weight be , and the bias be , then the final output of each node is: ; where Rn∈ represents the runoff prediction value sequence of node n in the future h time steps. The output is a continuous variable, with units of cubic meters per second, representing the intensity of surface runoff within the grid cell corresponding to the node. The network is trained by minimizing the mean squared error loss function between the predicted value and the measured value: ; where, is the measured value; the predicted values Rn of all nodes are reconstructed in node number order into a matrix R with dimensions N x h, representing the runoff prediction values of N grid cells in the future h time steps. For each time step t e {1, 2,..., h}, the t-th column of the matrix R is extracted to generate a complete gridded runoff prediction layer, corresponding to the spatial runoff distribution at the future t-th time step.

[0120] S6, for the grid cells with observation gaps, the dynamic interpolation completion is performed through the adjacent node feature propagation mechanism in the graph structure.

[0121] First, scan the feature dimension of each node in the input tensor at each time step, if any one of the key dynamic features (such as rainfall intensity, evapotranspiration rate, soil moisture) or all features of a node are missing values (represented as NaN or null), mark the node as an "observation gap node".

[0122] A mask matrix M e is constructed, where N is the total number of nodes and d is the feature dimension. The elements Mi,j in the mask matrix are defined as follows: if the i-th node has an observation value for the j-th feature, then Mi,j = 1; if there is a missing value, then Mi,j = 0. In the subsequent feature propagation operation, the missing position is identified using the mask matrix, and the initial feature vector of the missing node is initialized using one of the following two methods:

[0123] Zero vector filling method: set all the feature vectors Xi of the missing nodes to 0;

[0124] Average feature filling method: set the missing node features to the average of the non-missing nodes in the corresponding feature dimension.

[0125] In the graph structure G = (V, E) that has been constructed, the connection relationship between each node is represented by the normalized adjacency matrix . According to the propagation rule of graph convolution, the node feature update formula is as follows: ; where: X represents the initial node feature matrix (including missing values); W is the trained graph convolution layer weight; H is the output representation after node feature propagation.

[0126] During the forward propagation process, for the missing nodes, since their initial features are set to zero or average values, their updated values are completely propagated from the adjacent nodes through the weighted adjacency relationship, reflecting the spatial dependence features in the graph structure. The normalized adjacency matrix The configuration mode is as follows: .

[0127] In the graph convolution forward propagation output H, the updated feature vector of each node contains information from adjacent nodes. For the node with missing values, its updated feature vector is actually a weighted linear combination of the features of adjacent nodes.

[0128] This update result can be regarded as the interpolation completion vector of node i. In the missing feature dimension, write this output as the new valid feature of the node to the original feature matrix: ; wherein ⊙ represents element-wise multiplication, is the mask vector of node i, represents the completion feature vector of the i-th node after the graph convolution feature propagation and interpolation completion process. This operation ensures that non-missing values are retained and missing values are replaced by the graph convolution output.

[0129] Finally, a complete completion feature matrix is formed, which can be used for subsequent time series modeling steps, ensuring the integrity of the input of the subsequent LSTM network and avoiding interruptions in training or prediction caused by null values.

[0130] S7, if a sudden heavy rainfall event occurs in the target area, dynamically adjust the input feature sequence according to the real-time radar rainfall data, and update the spatial graph structure G in real time, and re-execute S3 to S6.

[0131] First, receive the real-time radar reflectivity data of the target area through the meteorological data interface, with a time resolution of not higher than 5 minutes and a spatial resolution of not lower than 1 kilometer. Convert using the reflectivity-rainfall intensity relationship: ; wherein R is the rainfall intensity (unit: millimeters per hour); Z is the radar reflectivity (unit: millimeters per cubic meter); and the commonly used empirical coefficients are a=0.0365, b=0.625.

[0132] Map the converted rainfall intensity value to the original grid node set corresponding to the spatial range in the form of a grid, 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 and support subsequent dynamic response operations.

[0133] According to the historical climate statistical data of the target area, set the heavy rainfall identification threshold Rth, for example, set it to 50 millimeters per hour in urban watersheds, and appropriately reduce it to 30 millimeters per hour in mountainous or small watersheds.

[0134] For each node vi, if the current radar rainfall intensity Ri in its mapping area is greater than or equal to Rth, mark it as a "heavy rainfall unit node" and add it to the set The marking process is performed in real time, ensuring that the model can quickly identify and respond to sudden extreme rainfall events.

[0135] For all nodes marked as heavy rainfall cells, their original input feature sequence is the dynamic meteorological and static geographical features at fixed time steps. Now, higher-frequency observation features need to be introduced in the time dimension. The specific steps are as follows:

[0136] Extract the newly added radar rainfall intensity sequence of the node within the last L time steps (e.g., 30 minutes, with 5-minute values, a total of 6 data), denoted as ;

[0137] At the same time, obtain the latest soil moisture inversion results (such as obtained by remote sensing inversion or ground observation) and fill them into each newly added time step;

[0138] Concatenate the above newly added data sequence to the end of the corresponding node time dimension of the original input tensor in chronological order 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 the node in the current prediction period, ensuring that the model captures the immediate impact of rainfall mutations on runoff response.

[0139] Since heavy rainfall can significantly change the local runoff path, it is necessary to recalculate the connection relationship between grid nodes in the heavy rainfall area to form a new spatial subgraph structure. The specific operation is as follows:

[0140] Extract the grid cells corresponding to the heavy rainfall node set Vrain;

[0141] Based on real-time elevation data or DEM data, calculate the main flow direction and slope value of each node in the region;

[0142] Determine the potential flow path of each node under heavy rainfall conditions, reconstruct the adjacency relationship based on the maximum slope value principle, construct a new edge set Erain, and then form a local subgraph Grain=(Vrain,Erain);

[0143] In the original graph structure G=(V,E), replace the edge set Eold corresponding to the node subset with Erain to generate the updated graph structure G′=(V,E′), where .

[0144] S8, output the high-resolution gridded runoff distribution map of the entire target watershed within a future set period, and predict the runoff process at the watershed and city levels.

[0145] Firstly, the runoff prediction results of each node in the future prediction period are sorted. Assuming that the total number of nodes is N, the target prediction step is H (for example, 6 hours in the future, and the time resolution is 1 hour, then H = 6), each node corresponds to a prediction vector with a length of H, which represents the surface runoff (unit: cubic meters per second) at each time step in the future.

[0146] The prediction vectors of all nodes are arranged in order according to the number to form a two-dimensional matrix R with dimensions N rows by H columns: ; wherein represents the runoff prediction value of the Nth node at the prediction time step H. The runoff value of each node in the prediction matrix is restored to its corresponding geographical location. In the aforementioned graph structure construction stage, a unique spatial coordinate mapping relationship for each node in the target watershed grid has been established, usually using grid center point coordinates or vector boundary center.

[0147] Let the spatial grid center corresponding to node i be (xi, yi), then at any prediction time step t, map to the two-dimensional spatial coordinates (xi, yi). This operation converts node data into a spatial distribution form.

[0148] The two-dimensional runoff distribution map obtained by interpolation at each prediction time step is taken as a layer, and a complete sequence of gridded runoff layers is generated in sequence of time steps: R = {R1(x, y), R2(x, y), …, RH(x, y)}; wherein each Rt(x, y) is a two-dimensional spatial matrix with the same resolution as the input grid (e.g. 100 meters x 100 meters per grid), and the numerical value represents the surface runoff intensity in the area unit.

[0149] This layer sequence can be used to visualize the runoff evolution process, showing the dynamic characteristics of runoff diffusion, convergence, peak movement, etc. through consecutive image frames. Administrative boundaries, river networks, infrastructure layers can be further superimposed to achieve geographical spatial correlation.

[0150] To meet the needs of urban flood control, reservoir scheduling and emergency response, the runoff process information of key areas needs to be extracted from the above runoff layer sequence.

[0151] Firstly, set the key monitoring unit set which can be determined as follows:

[0152] For urban areas: select nodes corresponding to historical waterlogging points, surrounding areas of drainage pumping stations, sunken interchanges, etc.

[0153] For mountainous watersheds: select nodes at the outlet of the catchment area or the lowest point of the terrain.

[0154] For any monitoring node vs∈S, extract its prediction sequence in the R matrix: Qs={rs,1,rs,2,…,rs,H}; draw the curve of Qs with time, obtain its peak time, peak size, runoff duration, rising and falling rate and other key indicators, which are used to realize the early perception of the development trend of local waterlogging or the evolution of flood peak.

[0155] The above merely provides a specific implementation of the present application, but the protection scope of the present application is not limited thereto, any person skilled in the art can easily think of changes or replacements within the technical range disclosed by the present application, which shall be covered within the protection scope of the present 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, specifically 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 areas in the original graph are replaced to form an updated spatial graph structure. 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 step of inputting the spatiotemporal feature representation DFst into a fully connected neural network and outputting the future runoff forecast values ​​for the 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 target prediction step size. The output of the hidden layer is mapped to the runoff forecast value 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: 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