A method and system for quantifying multi-step prediction uncertainty of watershed water quality based on a two-stream graph state space model

Abstract

The application discloses a kind of based on two-flow chart state space model's watershed water quality multi-step prediction uncertainty quantification method and system, comprising: based on water quality monitoring section and its observation data, construct watershed hydraulic topology chart and space-time tensor;Two-flow chart state space model is constructed, including: two-flow chart convolution space fusion module, time series evolution module based on state space model, quantile regression decoder and global spatiotemporal comprehensive loss function;Based on global spatiotemporal comprehensive loss function, the model is jointly trained and parameter optimization, simultaneously optimize the central tendency prediction and distribution boundary prediction of data in unified computational graph, and carry out uncertainty quantification index evaluation.The application can simultaneously fuse river network physical topology and long time series evolution law, and can output the method with statistical significance of prediction confidence interval, to realize the accurate quantification of future water quality state and its uncertainty.

Description

Technical Field

[0001] This invention belongs to the field of smart water conservancy and environmental monitoring technology, specifically involving a method for processing spatiotemporal sequence data using deep learning technology, and in particular a method and system for quantifying the uncertainty of multi-step water quality forecasting in watersheds based on a dual-flow graph state-space model. Background Technology

[0002] Watershed water quality forecasting is a core technical support for refined water environment management, water pollution prevention and control, and emergency response to sudden water incidents. Accurately forecasting the changing trends of key water quality indicators such as dissolved oxygen, ammonia nitrogen, and total phosphorus over a future period is of great significance for guiding the ecological operation of reservoirs and optimizing the operation of wastewater treatment plants.

[0003] In recent years, data-driven models, represented by Recurrent Neural Networks (RNNs), Long Short-Term Memory Networks (LSTMs), and Transformers, have made significant progress in the field of water quality prediction. However, in actual water environment supervision, existing prediction methods still face the following serious challenges. First, deterministic point predictions cannot reflect decision-making risks. Most existing water quality prediction models fall into the category of "point predictions," that is, they only output a deterministic value (mean) at a future time. However, influenced by factors such as variable meteorological conditions, random sewage discharge behavior, and complex aquatic biochemical reactions, water quality evolution has a high degree of randomness and uncertainty. A single point prediction cannot inform decision-makers of the reliability (confidence interval) of the predicted value, making it difficult for managers to assess the risk of exceeding standards when facing extreme weather or sudden pollution, often resulting in insufficient early warning or overreaction. Second, error accumulation and computational bottlenecks in long-sequence multi-step predictions. To meet the needs of flood control and other requirements, multi-step predictions that last for several days or even weeks are usually required. Traditional RNN / LSTM models use recursive inference, and the prediction error accumulates exponentially with the increase of the step size. While Transformer models can be computed in parallel, their computational complexity increases quadratically when processing extremely long historical sequences (such as capturing seasonal patterns), leading to huge memory consumption and making them difficult to deploy on edge computing devices. Third, spatial topology modeling lacks physical interpretability. Existing graph neural networks (GCNs) typically use undirected graphs or simple distance matrices to define station relationships when processing water quality data. This approach ignores the natural bidirectional asymmetric physical mechanisms of "downstream migration" and "upstream backwater / diffusion" in river networks, resulting in limited accuracy when processing time-lag relationships between upstream and downstream stations. Summary of the Invention

[0004] Purpose of the invention: One objective of this invention is to provide a method for quantifying the uncertainty of multi-step water quality forecasting in watersheds based on a dual-flow graph state-space model. This method can simultaneously integrate the physical topology of the river network and the long-term evolution law, and can output a statistically significant prediction confidence interval, so as to achieve accurate quantification of the future water quality state and its uncertainty.

[0005] Another objective of this invention is to provide a system for quantifying the uncertainty of multi-step watershed water quality forecasting based on a dual-flow graph state-space model.

[0006] Technical solution: The method described in this invention includes the following steps:

[0007] A watershed hydraulic topology map is constructed based on discretely distributed water quality monitoring sections and their multidimensional monitoring data; a dual-flow physical adjacency matrix is ​​constructed to describe the direct transmission effect of upstream water quality status on downstream, as well as the potential impact of downstream hydrological status on upstream; the multidimensional monitoring data is encapsulated into a spatiotemporal tensor containing four-dimensional information of batch, time, space, and features, including an input spatiotemporal tensor and a truth label spatiotemporal tensor.

[0008] A dual-flow graph state-space model is constructed, comprising: a dual-flow graph convolutional spatial fusion module, a temporal evolution module based on the state-space model, a quantile regression decoder, and a global spatiotemporal comprehensive loss function. The dual-flow graph convolutional spatial fusion module utilizes the dual-flow physical adjacency matrix to perform spatial-dimensional graph convolution operations on the input spatiotemporal tensor, extracting a high-dimensional spatial feature tensor containing river network spatial dependencies. The temporal evolution module based on the state-space model is used to perform deep temporal modeling on the high-dimensional spatial feature tensor, outputting a high-order spatiotemporal feature full tensor. The quantile regression decoder receives the high-order spatiotemporal feature full tensor, performs dimensional reshaping and multi-path decoding through a multilayer perceptron, and generates a global probability prediction output tensor containing multiple quantiles.

[0009] The model is jointly trained and its parameters are optimized based on the global spatiotemporal integrated loss function. The central trend prediction and distribution boundary prediction of the data are optimized simultaneously in a unified computation graph, and the uncertainty quantification index is evaluated.

[0010] Furthermore, a watershed hydraulic topology map is constructed based on the discretely distributed water quality monitoring sections and their multidimensional monitoring data; including: defining all water quality monitoring sections within the target watershed as graph nodes, defining the natural channels connecting each section as graph edges, and determining the direction of the edges based on the water flow direction between graph nodes;

[0011] The two-stream physical adjacency matrix includes: the downstream migration matrix and the upstream diffusion matrix. In the watershed hydraulic topology graph, if there is an edge from an upstream node to a downstream node, the elements of the downstream migration matrix are 1; otherwise, they are 0. The upstream diffusion matrix is ​​the transpose of the downstream migration matrix.

[0012] Slices in the input spacetime tensor Representing the In the nth sample, the nth The time step, the first All standardized feature vectors of each node; slices in the truth label spacetime tensor Representing the In the sample, the future number The prediction time step, the first All standardized feature vectors of each node.

[0013] Furthermore, the two-stream graph convolutional spatial fusion module utilizes the two-stream physical adjacency matrix to perform spatial-dimensional graph convolution operations on the input spatiotemporal tensor, extracting a high-dimensional spatial feature tensor containing river network spatial dependencies; including:

[0014] (1) Using the normalized downstream migration matrix, perform downstream graph convolution on the features of each graph node to obtain the downstream feature matrix. : ;in, The learnable weight matrix of the forward convolution kernel. For the hidden layer dimension, For the bias term of the co-current convolution kernel, It is a non-linear activation function. For the normalized downstream migration matrix, For any sample and any time step The corresponding node feature matrix;

[0015] (2) Using the normalized inverse diffusion matrix, perform inverse graph convolution on the features of each graph node to obtain the inverse feature matrix. : ;in, The independent, learnable weight matrix of the inverse convolution kernel. This is the bias term of the countercurrent convolution kernel. This is the normalized countercurrent diffusion matrix;

[0016] (3) The upstream and downstream features are concatenated along the feature dimension, and node-level fusion weight coefficients are generated through a lightweight multilayer perceptron (MLP). ; Utilizing fusion weight coefficients The two sets of features are dynamically weighted and summed to obtain the fused spatial context feature matrix;

[0017] (4) Repeat steps (1) to (3) above for all batches and all time steps of the input spatiotemporal tensor to finally obtain a high-dimensional spatial feature tensor that has fully extracted the spatial dependencies of the river network.

[0018] Furthermore, the temporal evolution module based on the state-space model is used for deep temporal modeling of high-dimensional spatial feature tensors, outputting a high-order spatiotemporal feature full tensor; including:

[0019] (1) Define a continuous linear time-invariant state equation. Decompose the high-dimensional spatial feature tensor output by the dual-flow graph convolutional spatial fusion module along the feature channel. Treat the time series of any node and any feature channel as a continuous dynamic system controlled by the differential equation. The input signal of the state equation is the value of a specific node and a specific feature channel of a certain sample in the high-dimensional spatial feature tensor at a certain time. The evolution of the corresponding latent state vector and output signal follows the continuous linear time-invariant state equation.

[0020] (2) Diagonalize the state transition matrix in the continuous linear time-invariant state equation and initialize the diagonal elements of the state transition matrix with frequency tuning.

[0021] (3) Discretize the continuous linear time-invariant state equations using the zero-order preservation principle to obtain discrete recursive equations;

[0022] (4) Transform the discrete recursive equation into an equivalent global convolution form. Stack the input and output signals of all discrete time steps in time order to construct the input sequence vector and output sequence vector. Use Fast Fourier Transform and the convolution kernel generated by the convolution kernel generator to perform frequency domain convolution calculation on the input and output sequence vectors. Define the frequency domain convolution calculation process from the input sequence vector to the output sequence vector as... Operator;

[0023] (5) Operators are embedded into deep neural network architectures to construct structures containing... A multi-layer stacked deep network, where the input tensor of the 0th layer is a high-dimensional feature tensor, and for the 1st layer of the network... layer, Its input is the output of the previous layer; extract the first... The layer outputs the state vector at the last time step, which forms the high-order spatiotemporal feature full tensor.

[0024] Furthermore, the quantile regression decoder decodes the full tensor of high-order spatiotemporal features into a probability distribution representation of all water quality variables across multiple future time steps; including:

[0025] (1) For the first in the batch The first sample and the first For each graph node, its corresponding feature vector is extracted from the high-order spatiotemporal feature full tensor output by the state-space model-based temporal evolution deduction module. ,Will The input decoder expands the dimension to the total dimension required for prediction through two layers of linear transformation, resulting in a flattened prediction output vector. The flattened prediction output vector is then reshaped to obtain the local prediction tensor of the node.

[0026] (2) Define a containing The set of quantiles of elements, and select These are key quantile anchor points, each corresponding to the last dimension of the local prediction tensor;

[0027] (3) Stack the local prediction tensors of all batches and all nodes in the spatial dimension, and swap the time and node dimensions to construct the final global probability prediction output tensor;

[0028] (4) The confidence interval of any variable can be directly analyzed from the output tensor of global probability prediction.

[0029] Furthermore, the model is jointly trained and its parameters are optimized based on a global spatiotemporal integrated loss function. This simultaneously optimizes the prediction of data central tendency and distribution boundary within a unified computational graph, and evaluates the uncertainty using quantitative indicators; including:

[0030] (1) Define the quantile regression loss function, let For the true value, For probabilistic forecasting models targeting quantiles Calculate the single-point loss based on the predicted value. :

[0031] ;

[0032] (2) Aggregate the losses of all batch steps, nodes, times, variables and quantiles to construct a global spatiotemporal comprehensive loss function;

[0033] (3) The gradient-based optimization algorithm is used to minimize the global spatiotemporal integrated loss function. During the training process, the gradient is calculated using the backpropagation algorithm, and all learnable parameters of the spatial layer, temporal layer and decoding layer are updated at the same time.

[0034] (4) After the model training is completed, the test set data is input into the dual-flow graph state space model to generate the probability prediction results of the test set and verify its uncertainty quantification capability.

[0035] The system corresponding to the method includes:

[0036] The data processing unit is used to construct a watershed hydraulic topology map based on discretely distributed water quality monitoring sections and their multidimensional monitoring data; construct a dual-flow physical adjacency matrix to describe the direct transmission effect of upstream water quality status on downstream, as well as the potential impact of downstream hydrological status on upstream; and encapsulate the multidimensional monitoring data into a spatiotemporal tensor containing four-dimensional information of batch, time, space, and features, including an input spatiotemporal tensor and a truth label spatiotemporal tensor.

[0037] The model building unit is used to construct a dual-flow graph state-space model, including: a dual-flow graph convolutional spatial fusion module, a temporal evolution module based on the state-space model, a quantile regression decoder, and a global spatiotemporal comprehensive loss function. The dual-flow graph convolutional spatial fusion module uses the dual-flow physical adjacency matrix to perform spatial-dimensional graph convolution operations on the input spatiotemporal tensor to extract a high-dimensional spatial feature tensor containing river network spatial dependencies. The temporal evolution module based on the state-space model is used to perform deep temporal modeling on the high-dimensional spatial feature tensor and output a high-order spatiotemporal feature full tensor. The quantile regression decoder receives the high-order spatiotemporal feature full tensor, performs dimensional reshaping and multi-path decoding through a multilayer perceptron, and generates a global probability prediction output tensor containing multiple quantiles.

[0038] The model training and optimization prediction unit is used to jointly train the model and optimize the parameters based on the global spatiotemporal integrated loss function. It simultaneously optimizes the central trend prediction and distribution boundary prediction of the data in a unified computation graph and evaluates the uncertainty quantification index.

[0039] An electronic device for storing and executing the method includes:

[0040] Memory, used to store computer programs;

[0041] A processor for executing the computer program to implement the method.

[0042] A non-volatile storage medium for storing and executing the method, for storing a computer program, wherein the computer program implements the method when executed by a processor.

[0043] The method corresponds to a computer program product, which includes a computer program / instruction that, when executed by a processor, implements the method.

[0044] Beneficial effects: The significant technical effects of this invention include: (1) It realizes the leap from "numerical estimation" to "risk quantification" in water quality prediction; This invention abandons the traditional point prediction mode and innovatively introduces quantile regression decoder and bouncing loss function; The model can not only output the most likely value (median) of future water quality, but also generate a confidence interval with statistical guarantee at the same time; This allows water managers to not only see the "predicted concentration", but also intuitively grasp the "range of uncertainty of prediction": the interval narrows during the stable period, indicating that the prediction is credible; the interval diverges during the sudden change period, indicating that the risk increases. This risk-aware characteristic provides a key basis for scientific decision-making; (2) It breaks through the computational bottleneck of long sequence modeling and improves the accuracy of multi-step prediction. This invention adopts the frequency domain tuned diagonal state space model (S4D-FT) as the core of time evolution, replacing the traditional LSTM or Transformer. By utilizing the characteristic of global convolution in the frequency domain of the S4D model, the time complexity of long sequence modeling is reduced from Reduce to While ensuring extremely low computational resource consumption, it can effectively capture historical dependencies (such as seasonal patterns) over thousands of time steps, significantly reducing the error accumulation problem in multi-step prediction; (3) It integrates hydrodynamic mechanisms, enhancing the physical interpretability of spatial features; This invention constructs a dual-flow graph convolutional network based on river network flow direction. By designing "downstream migration matrix" and "upstream diffusion matrix" respectively, the model can explicitly simulate the process of pollutants migrating downstream with the water flow, as well as the reverse influence of downstream water level on upstream. With the adaptive gating fusion mechanism, the model can automatically The dynamic adjustment of the weight of upstream and downstream information according to the characteristics of the river section (such as rapid flow or slow flow) makes the spatial feature extraction process more in line with the laws of hydrology and physics; (4) It has excellent robustness and generalization ability; through end-to-end joint training, this model can adaptively learn the spatiotemporal evolution patterns of different monitoring stations and different water quality indicators; even if some monitoring data is missing or there is noise interference, the model can still give a robust prediction range by relying on the powerful spatiotemporal context reasoning ability and probability output mechanism, avoiding the problem of severe oscillation false alarms caused by individual abnormal inputs in deterministic models. Attached Figure Description

[0045] Figure 1 This is a flowchart of the method of the present invention;

[0046] Figure 2 This is a schematic diagram of multi-step probability prediction and uncertainty quantification of ammonia nitrogen based on a dual-flow graph state-space model.

[0047] Figure 3 This is a schematic diagram comparing the uncertainty quantification performance of the dual-flow graph state-space model and the classical model under multi-step forecasting. Detailed Implementation

[0048] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments.

[0049] like Figure 1 As shown, the method of the present invention includes the following steps:

[0050] Step S1: Constructing the watershed hydraulic topology map and spatiotemporal tensor. This step aims to abstract the discretely distributed water quality monitoring sections and their multidimensional monitoring data into a computer-processable graph structure and high-dimensional tensor, and explicitly encode the upstream and downstream hydraulic connections of the river network. Specifically, this includes:

[0051] S1-1: Topological abstraction of the watershed monitoring network. All water quality monitoring sections within the target watershed are defined as graph nodes, and the natural channels connecting each section are defined as graph edges, constructing a watershed hydraulic topology graph. Assume there are a total of [number missing] water quality monitoring sections within the watershed. A graph node defines a directed graph. ,in For a set of nodes, They are respectively Each graph node For the edge set, use a Geographic Information System (GIS) or Digital Elevation Model (DEM) to determine the direction of water flow between graph nodes. If the water flow can pass through the node... Direct flow to nodes Then in Establish a line from point to The directed edges, where As an upstream node, It is a downstream node.

[0052] S1-2: Constructing a Two-Stream Physical Adjacency Matrix. To simulate the two physical mechanisms of "downstream migration" and "counterstream diffusion" of water pollutants in the river network, this embodiment constructs a two-stream physical adjacency matrix, including a downstream migration matrix and a counterstream diffusion matrix. First, the downstream migration matrix is ​​constructed. This matrix describes the direct transmission effect of upstream water quality conditions on downstream areas. Downstream migration matrix elements. The definition is as follows:

[0053] ;

[0054] To prevent gradient explosion and ensure numerical stability, Perform row normalization:

[0055] ;

[0056] in, For the normalized downstream migration matrix, for The inverse matrix, For an out-degree diagonal matrix, its diagonal elements .

[0057] Then construct the countercurrent diffusion matrix. This matrix is ​​used to describe the potential impact of downstream hydrological conditions (such as water level backflow and backflow effect) on upstream areas, or for information backpropagation in graph neural networks to enhance information connectivity across the entire basin. Defined as The transpose of is:

[0058] ;

[0059] in, for The transpose of the matrix;

[0060] Similarly to The normalized countercurrent diffusion matrix is ​​obtained by performing row normalization. .pass and This invention can aggregate spatial information from two dimensions: "causal chain" (upstream affects downstream) and "correlation chain" (downstream feeds back to upstream), which can better reflect the physical nature of water flow than traditional undirected graphs.

[0061] S1-3: Multidimensional Monitoring Data Fusion and Spatiotemporal Tensor Construction. Multidimensional monitoring data is encapsulated into a standard tensor containing four dimensions: batch, time, space, and features. For each node... Collect its in Multidimensional monitoring data at all times The multidimensional monitoring data includes endogenous and exogenous features. The feature channel number represents the total number of endogenous and exogenous features. Endogenous features are historical water quality concentrations, including dissolved oxygen and ammonia nitrogen. Exogenous features are meteorological and hydrological driving factors, including rainfall, water level, and flow velocity. The Z-Score standardization method is used to remove dimensions from the multidimensional monitoring data, ensuring a mean of 0 and a variance of 1, resulting in the standardized feature vector. For any sampling time Construct a spacetime tensor, including the input spacetime tensor. And truth label spacetime tensor .in, For training batch size, To predict the step size, The length of the historical backtracking window. The number of nodes. A slice in the input spacetime tensor. Representing the In the nth sample, the nth The time step, the first All standardized feature vectors of nodes ; slices in the spacetime tensor with truth labels Representing the In the sample, the future number The prediction time step, the first All standardized feature vectors of nodes .

[0062] Step S2: Construct a dual-flow graph state-space model. This step aims to build an end-to-end deep learning architecture to extract spatiotemporal evolution patterns and output probability prediction results including uncertainty quantification. The dual-flow graph state-space model includes, in sequence according to the data flow direction, a dual-flow graph convolutional spatial fusion module, a temporal evolution module based on the state-space model (SSM), a quantile regression decoder, and a global spatiotemporal comprehensive loss function; wherein, the dual-flow graph convolutional spatial fusion module receives the input spatiotemporal tensor from step S1. The dual-stream physical adjacency matrix, by aggregating observation information from upstream and downstream stations, outputs a high-dimensional spatial feature tensor that incorporates global spatial correlations. The temporal evolution module based on the state-space model receives the high-dimensional spatial feature tensor. Global convolutional evolution is performed in the time dimension to extract long-sequence historical memories, and finally the high-order spatiotemporal feature full tensor after dimensionality reduction is output. The quantile regression decoder receives the full tensor of the aforementioned high-order spatiotemporal features. By using a multilayer perceptron for dimensionality reshaping and multi-path decoding, a global probability prediction output tensor containing multiple quantiles is generated. The lower bound, median, and upper bound of the prediction interval are constructed. The global spatiotemporal integrated loss function receives the global probability prediction output tensor and compares it with the true value label tensor. By calculating the asymmetric bouncing loss, the total loss value is output, forcing the prediction quantile to approach the true probability distribution, thereby driving the model parameter update in subsequent steps.

[0063] Step S2-1: Construct a two-stream graph convolutional spatial fusion module. This step aims to utilize the two-stream physical adjacency matrix constructed in step S1 to perform spatial fusion on the input spatiotemporal tensor. Spatial-dimensional graph convolution operations are performed. The core of this approach lies in simulating the physical processes of downstream pollutant diffusion and the feedback mechanism from downstream to upstream, thereby extracting feature representations that integrate global spatial topological information of the watershed. Specifically, this includes:

[0064] S2-1-1: Downstream Feature Aggregation. Utilizing the normalized downstream migration matrix. Spatial convolution is performed on the features of each graph node. The physical meaning of this operation is: the current node... It will receive all its direct upstream nodes. The water quality information is used to simulate the process of water flow transporting pollutants from upstream to downstream. Assume that in any sample... and any time step The node feature matrix in the graph is The formula for single-layer convolution of a downstream graph is as follows:

[0065] ;

[0066] in, The characteristic matrix is ​​for downstream flow. The learnable weight matrix of the forward convolution kernel. For the hidden layer dimension, For the bias term of the co-current convolution kernel, It is a non-linear activation function (such as ReLU or GELU).

[0067] S2-1-2: Backflow Feature Aggregation. Utilizing the normalized backflow diffusion matrix. Then, perform spatial convolution similar to that in step S2-1-1, which aggregates downstream features. The physical significance of this operation is to capture the backflow effect of pollutants caused by downstream water level rise, or to invert the upstream state through downstream measuring points in areas without data. The formula for the convolution operation of the upstream map is as follows:

[0068] ;

[0069] in, The independent, learnable weight matrix of the inverse convolution kernel. The reverse flow characteristic matrix, This is the bias term for the inverse convolution kernel.

[0070] S2-1-3: Adaptive Dual-Flow Feature Fusion. Due to the different hydrodynamic conditions of river sections where different monitoring stations are located (e.g., rapid current sections are greatly influenced by downstream flow, while slow-flowing / tidal sections are significantly influenced by upstream flow), the contribution weights of downstream and upstream flow information to prediction are not equal. Therefore, this embodiment designs an adaptive gating fusion module. First, feature stitching and gating weight calculation are performed. The downstream feature matrix... and reverse flow characteristic matrix The components are concatenated along the feature dimension, and node-level fusion weight coefficients are generated using a lightweight multilayer perceptron (MLP). :

[0071] ;

[0072] in, Represents a non-linear activation function. Describes the multilayer perceptron function. This indicates a feature concatenation operation, which merges weight coefficients. It tends to follow the flow characteristics. Then, a weighted fusion is performed.

[0073] Using fusion weight coefficients The two features are dynamically weighted and summed to obtain the fused spatial context feature matrix. :

[0074] ;

[0075] in, This represents element-wise multiplication under the broadcast mechanism.

[0076] S2-1-4: Iterative application of batch and time dimensions. For the input spatiotemporal tensor... All samples and all time steps Repeat steps S2-1-1 to S2-1-3 above. Finally, a high-dimensional spatial feature tensor containing the complete extracted spatial dependencies of the river network is obtained. :

[0077] ;

[0078] After this step, the high-dimensional space feature tensor Each element in the graph (representing the characteristics of a graph node at a certain moment) not only contains its own observation information, but also aggregates the water quality status of its upstream and downstream stations. This "dual-flow mechanism" effectively overcomes the problem of information transmission bias in directed water systems in traditional graph neural networks, laying the foundation for processing sparse data regions (by "borrowing" information from neighboring nodes).

[0079] Step S2-2: Temporal Evolution Module Based on State-Space Model (SSM). This step aims to process the high-dimensional spatial feature tensor output from step S2-1-1. Deep modeling is performed in the time dimension. To address the issues of traditional recurrent neural networks (RNNs) being unable to be trained in parallel and the vanishing gradient problem in long sequences, this embodiment constructs a temporal evolution module based on a state-space model.

[0080] S2-2-1: Define the continuous linear time-invariant (LTI) state equation. The high-dimensional feature tensor output from step S2-1... Decompose along the characteristic channels. Consider the time series of any node and any characteristic channel as a continuous dynamic system governed by differential equations. Let the input signal... for A specific node and a specific feature channel of a certain sample at time [time] The value of is denoted as and its corresponding potential state vector is denoted as . The output signal is denoted as The system's evolution follows the following continuous linear time-invariant (LTI) state equations:

[0081] ;

[0082] ;

[0083] in, for The derivative with respect to time, The state transition matrix represents the evolution dynamics of the control system. The preset state space order (e.g., 64); For input control matrix; To output the observation matrix; This is a direct join item.

[0084] S2-2-2: Diagonalization of the State Transition Matrix and Frequency Tuning Initialization. To reduce computational complexity and improve numerical stability, this embodiment employs a diagonalized state space (S4D) strategy. The first step is to perform diagonalization constraints. Define the state transition matrix. It is a diagonal matrix, that is , Let be the system's eigenvalue. This allows the above continuous linear time-invariant (LTI) state equation to be decomposed into: Each independent scalar differential equation is solved in parallel. Then, frequency tuning initialization is performed. Considering the long-term and short-term dependencies of water quality data, [further details are needed]. The diagonal elements are initialized with a specific distribution. Specifically, real part Initialize to negative values ​​to ensure system stability (state decay), imaginary part It is initialized with frequency components corresponding to different time scales, thereby covering a variety of frequency features from hourly abrupt changes to seasonal trends.

[0085] S2-2-3: Discretization Derivation Based on Zero-Order Preservation. Since water quality monitoring data is sampled discretely at fixed time intervals (e.g., 1 hour), the continuous state equation must be converted into a discrete recursive form. A learnable time step parameter is introduced. (Represents the physical scaling factor of the sampling interval). Discretization is performed using the zeroth-order hold principle, assuming that at discrete time steps... arrive Between (time step is) ), input signal Keep constant, that is For the state equation In the interval Integrating above:

[0086] ;

[0087] in, for The latent state vector at time t, for The latent state vector at time t, The input signal is within the integration interval. Let it be the integration variable; Recorded as , Recorded as Substituting the zeroth-order preservation assumption Simplifying, we get:

[0088] ;

[0089] in, The discretized state transition matrix is ​​calculated using the following formula: ; The formula for calculating the discretized input matrix is ​​as follows: , express The inverse matrix, Represents the identity matrix. When the diagonalization strategy in step S3-2 is used to make... In this case, the above matrix operations can be simplified to operating on each element of the vector. This allows for scalar operations, thus avoiding the need to invert high-dimensional matrices. The diagonal elements and The The individual components are calculated as follows:

[0090] ;

[0091] ;

[0092] in, express The One diagonal element, Indicates the first A system characteristic value, express The One portion, For matrix The Each element. Using the above formula, the dual-flow graph state-space model avoids complex matrix inversion operations, thus significantly reducing computational complexity. At sampling time... Reading the output, the output equation can be directly discretized as follows:

[0093] ;

[0094] in, For a moment The output signal, , .

[0095] S2-2-4: Parallel computation view based on full sequence convolution. To solve the discrete recursive equation in step S2-2-3 ( When processing long sequences, it is necessary to proceed by time step. Addressing the issues of sequential computation and the inability to train in parallel, this step transforms it into an equivalent global convolutional form. First, all discrete time steps defined in step S2-2-3 are... input signal and output signal Stack them in chronological order. Construct the input sequence vector and output sequence vector :

[0096] ;

[0097] ;

[0098] At this point, the entire output sequence vector is solved. The process can be equivalently transformed into an input sequence vector. Convolution operation with a global system impulse response (convolution kernel). Assume the initial state... The zero vector (if non-zero initial states need to be considered, it can be regarded as a superposition of bias terms), according to the properties of linear time-invariant systems, output sequence vector. It can be represented as:

[0099] ;

[0100] in, This represents a one-dimensional causal convolution operation. This represents element-wise multiplication of directly joined terms. Convolution kernel. The calculation formula is:

[0101] ;

[0102] in, This represents the convolution kernel generation operator, the first... element Represents the moment The input signal evolves through the internal state of the system. After the step, for the current moment The contribution weights of the output signal. To reduce the high computational cost of directly calculating time-domain convolution, the calculation is performed in the frequency domain using the convolution theorem and Fast Fourier Transform (FFT):

[0103] ;

[0104] in, Indicates fill to Fast Fourier Transform after length to avoid circular convolution artifacts; This represents element-wise multiplication in the frequency domain; Indicates the inverse Fourier transform; Indicates the truncation before There are one effective time step. For ease of subsequent description, the above will be derived from the input sequence vector. To output sequence vector The frequency domain convolution calculation process is defined as follows: Operator:

[0105] .

[0106] S2-2-5: Construction and Feature Output of Multi-Layer Stacked Deep Networks. (The above...) Operators are embedded into deep neural network architectures to construct structures containing... A multi-layer stacked deep network. Define the input tensor of the 0th layer of the multi-layer stacked deep network. The high-dimensional space feature tensor output in step S2-1 ,Right now:

[0107] ;

[0108] For the first multi-layer stacked deep network layer Its input is the output of the previous layer. The output is The calculation process employs residual connectivity and layer normalization. The calculation formula is as follows:

[0109] ;

[0110] ;

[0111] in, This is an intermediate transition tensor, representing the feature representation after state-space evolution and mixing but before being processed by an activation function. For layer normalization operation, The activation function for the Gaussian error linear unit;

[0112] In the formula In the calculation, the system will use tensors Each feature channel in (i.e., for each sample) Each node Each feature dimension time series These are respectively used as the input sequence vectors in step S2-2-4. Frequency domain convolution operations are performed independently and in parallel to obtain the corresponding output sequence vectors. After combination, it forms To map historical information to the future, the first layer of a multi-layer stacked deep network is extracted. Layer output The state vector at the last time step. This state vector has passed through... Long-distance memory compression aggregates information from the entire historical window. These state vectors form a high-order spatiotemporal feature full tensor. :

[0113] ;

[0114] in, For multi-layer stacked deep networks The deep feature tensor output by the layer, This indicates that the last index is taken in the time dimension.

[0115] Step S2-3: Construct a quantile regression decoder. This step aims to convert the high-order spatiotemporal features output from step S2-2 into a full tensor. Decoding the Future Each time step, all Characterization of the probability distribution of each water quality variable.

[0116] S2-3-1: Multidimensional Decoding Mapping of Spatiotemporal Features. Establish a mapping relationship from the latent feature space to the physical prediction space. Utilize a multilayer perceptron with shared weights as the decoder to perform parallel decoding of the high-order features of each graph node. For the first node in the batch... The first sample and the first Each graph node, from Extract its corresponding vector :

[0117] ;

[0118] in, This represents the input feature vector of the decoder. Tensor No. The sample, the first A slice of a graph node;

[0119] Will The input decoder reduces the dimension from [previous dimension] through two layers of linear transformation. Expanding to the total dimensions required for prediction:

[0120] ;

[0121] in, and These are the parameters for the decoder's intermediate layer; It is the decoder output layer weight matrix. Indicates the decoder bias; The predicted output vector is flattened. The number of quantiles. The flattened prediction output vector. Perform dimensional reshaping to obtain the local prediction tensor of the node. Elements in the local prediction tensor Representing the The sample, the first The site, in the future The prediction time step, for the first The first water quality variable, the first Predicted values ​​for each quantile.

[0122] S2-3-2: Definition and Physical Meaning of Quantile Anchor Points. To clarify the physical statistical meaning of the output values ​​of the dual-flow graph state-space model, a quantile anchor point needs to be predefined. quantile set of elements ,in and In this embodiment, the following is set: For local prediction tensors The last dimension, its physical meaning is defined as follows: Corresponding lower bound anchor point ( This indicates that, given historical information, the probability that the actual value is lower than the predicted value is 10%. Corresponding median anchor point ( This indicates that the probability of the true value being lower or higher than this value is 50%. Corresponding upper bound anchor point ( This indicates that the probability of the actual value being lower than the predicted value is 90% (i.e., the probability of it being higher than the predicted value is only 10%).

[0123] S2-3-3: Construction of the global probability prediction output tensor. (This applies to all batches.) and all nodes Local prediction tensor Stacking the elements in the spatial dimension and swapping the time and node dimensions, we construct the final global probability prediction output tensor, denoted as . .

[0124] S2-3-4: Generation and analysis of prediction confidence intervals. Output tensor based on global probability prediction. This allows direct parsing of the confidence interval for any variable. For example, for the first... The sample, the first The site, the The 80% confidence interval for each variable at the h-th step in the future. Represented as:

[0125] ;

[0126] in, Tensor The middle position is ( ) elements, Tensor The middle position is ( () elements.

[0127] Step S2-4: Joint Model Training and Uncertainty Assessment Based on Bouncing Loss. This step aims to construct an asymmetric loss function system to simultaneously optimize the central tendency prediction and distribution boundary prediction of the data in a unified computational graph.

[0128] S2-4-1: Definition of the quantile regression loss function. This is used to train a dual-flow graph state-space model to output specific quantiles. A quantile regression loss function is constructed using the bouncing ball loss function. The bouncing ball loss function is an asymmetric absolute value error that applies different weights to the prediction bias based on its direction (overestimation or underestimation), thereby forcing the predicted value to approximate a specific quantile of the true distribution. Let... For the true value, For the state-space model of the two-flow graph, the quantiles The predicted value, single-point loss The calculation formula is as follows:

[0129] ;

[0130] S2-4-2: Construction of the Global Spatiotemporal Comprehensive Loss Function. To achieve end-to-end training, the losses from all batch steps, nodes, times, variables, and quantiles are aggregated to construct the global spatiotemporal comprehensive loss function. Based on the truth label spatiotemporal tensor constructed in step S1 The global probability prediction output tensor output from step S2-4 is Construct a global spatiotemporal integrated loss function :

[0131] ;

[0132] in, For the first The preset value of the quantile. Tensor The middle position is ( ) elements, Tensor The middle position is ( () elements.

[0133] Step S3: Joint training and parameter optimization of the model based on the global spatiotemporal integrated loss function. Gradient-based optimization algorithms (such as the AdamW optimizer) are used to minimize... During training, the gradient is calculated using the backpropagation algorithm, and all learnable parameters of the spatial layer, temporal layer, and decoding layer are updated simultaneously. The spatial layer parameters include the two-stream convolutional kernel weights from step S2-1. and etc.; the time-level parameters include the state-space model discretization parameters in step S2-2. , , , , The decoding layer parameters include the multilayer perceptron weights from steps S2-3. , And bias.

[0134] Step S4: Uncertainty Quantification Index Evaluation. After the dual-flow graph state-space model is trained, the test set data is input into the dual-flow graph state-space model to generate the probability prediction results for the test set and verify its uncertainty quantification capability. Based on the definition in step S2-3-2, tensor slicing... This is the lower bound prediction value. This is the upper bound of the predicted value. This embodiment uses two core metrics to measure the effectiveness of uncertainty quantification: Predicted Interval Coverage (PICP) and Average Predicted Interval Width (MPIW). Predicted Interval Coverage is used to assess the reliability of the confidence interval (i.e., whether it truly contains the true value). For a given confidence level (e.g., 80%), the closer the Predicted Interval Coverage (PICP) is to 80% (or higher), the more reliable the prediction.

[0135] ;

[0136] in, This is the indicator function (1 if the condition is met, 0 otherwise). The average prediction interval width (MPIW) is used to evaluate the sharpness of the confidence interval (i.e., whether the interval is narrow enough). Under the premise of meeting coverage requirements, the smaller the MPIW, the lower the uncertainty and the higher the prediction accuracy.

[0137] ;

[0138] Through step S4, this invention not only achieves accurate point prediction of future water quality concentrations (through median optimization), but more importantly, it generates a statistically guaranteed "dynamic risk corridor." During periods of low water levels, data fluctuations are small, and the dual-flow diagram state space model automatically converges to a smaller MPIW (narrow interval), providing a high-confidence forecast. During periods of dramatic hydrological changes (such as rainstorms and sewage discharge), the dual-flow diagram state space model automatically outputs a larger MPIW (wide interval), accurately reflecting the current forecast uncertainty and indicating to managers that the forecast results at this time have a greater risk of fluctuation and need to be judged in conjunction with human experience.

[0139] Another embodiment of the present invention also provides a system corresponding to the method, comprising:

[0140] The data processing unit is used to construct a watershed hydraulic topology map based on discretely distributed water quality monitoring sections and their multidimensional monitoring data; construct a dual-flow physical adjacency matrix to describe the direct transmission effect of upstream water quality status on downstream, as well as the potential impact of downstream hydrological status on upstream; and encapsulate the multidimensional monitoring data into a spatiotemporal tensor containing four-dimensional information of batch, time, space, and features, including an input spatiotemporal tensor and a truth label spatiotemporal tensor.

[0141] The model building unit is used to construct a dual-flow graph state-space model, including: a dual-flow graph convolutional spatial fusion module, a temporal evolution module based on the state-space model, a quantile regression decoder, and a global spatiotemporal comprehensive loss function. The dual-flow graph convolutional spatial fusion module uses the dual-flow physical adjacency matrix to perform spatial-dimensional graph convolution operations on the input spatiotemporal tensor to extract a high-dimensional spatial feature tensor containing river network spatial dependencies. The temporal evolution module based on the state-space model is used to perform deep temporal modeling on the high-dimensional spatial feature tensor and output a high-order spatiotemporal feature full tensor. The quantile regression decoder receives the high-order spatiotemporal feature full tensor, performs dimensional reshaping and multi-path decoding through a multilayer perceptron, and generates a global probability prediction output tensor containing multiple quantiles.

[0142] The model training and optimization prediction unit is used to jointly train the model and optimize the parameters based on the global spatiotemporal integrated loss function. It simultaneously optimizes the central trend prediction and distribution boundary prediction of the data in a unified computation graph and evaluates the uncertainty quantification index.

[0143] Another embodiment of the present invention also provides an electronic device, comprising:

[0144] Memory, used to store computer programs;

[0145] A processor for executing the computer program to implement the method.

[0146] Another embodiment of the present invention provides a non-volatile storage medium for storing a computer program, wherein the computer program implements the method described when executed by a processor.

[0147] Another embodiment of the present invention provides a computer program product, including a computer program / instructions that, when executed by a processor, implement the method described.

[0148] like Figure 2 As shown, this invention not only outputs the median predicted value of future ammonia nitrogen concentration, but also simultaneously outputs the corresponding prediction confidence interval. When the predicted object enters a rapid rise or fall phase, the prediction interval given by the model widens significantly, indicating that this invention can sense the increased uncertainty during abrupt changes; in phases with smaller fluctuations, the prediction interval converges relatively, indicating that the model has a high confidence level for stable processes. Figure 3 As shown, under the statistical significance of the test set, with the increase of the prediction step size, the prediction interval coverage of both the present invention and the comparison model shows a decreasing trend, while the average prediction interval width shows an increasing trend. This is a normal phenomenon caused by the accumulation of errors in multi-step forecasting. Compared with the comparison model, the present invention generally exhibits higher prediction interval coverage and smaller average prediction interval width at each prediction step size, indicating that the present invention can maintain high prediction accuracy while ensuring interval reliability.

Claims

1. A method for quantifying the uncertainty of multi-step watershed water quality forecasting based on a dual-flow graph state-space model, characterized in that, Includes the following steps: A watershed hydraulic topology map is constructed based on discretely distributed water quality monitoring sections and their multidimensional monitoring data; a dual-flow physical adjacency matrix is ​​constructed to describe the direct transmission effect of upstream water quality status on downstream, as well as the potential impact of downstream hydrological status on upstream; the multidimensional monitoring data is encapsulated into a spatiotemporal tensor containing four-dimensional information of batch, time, space, and features, including an input spatiotemporal tensor and a truth label spatiotemporal tensor. A dual-flow graph state-space model is constructed, comprising: a dual-flow graph convolutional spatial fusion module, a temporal evolution module based on the state-space model, a quantile regression decoder, and a global spatiotemporal comprehensive loss function. The dual-flow graph convolutional spatial fusion module utilizes the dual-flow physical adjacency matrix to perform spatial-dimensional graph convolution operations on the input spatiotemporal tensor, extracting a high-dimensional spatial feature tensor containing river network spatial dependencies. The temporal evolution module based on the state-space model is used to perform deep temporal modeling on the high-dimensional spatial feature tensor, outputting a high-order spatiotemporal feature full tensor. The quantile regression decoder receives the high-order spatiotemporal feature full tensor, performs dimensional reshaping and multi-path decoding through a multilayer perceptron, and generates a global probability prediction output tensor containing multiple quantiles. The quantile regression decoder decodes the full tensor of high-order spatiotemporal features into a probability distribution representation of all water quality variables across multiple future time steps; including: (1) For the first in the batch The first sample and the first For each graph node, its corresponding feature vector is extracted from the high-order spatiotemporal feature full tensor output by the state-space model-based temporal evolution deduction module. ,Will The input decoder expands the dimension to the total dimension required for prediction through two layers of linear transformation, resulting in a flattened prediction output vector. The flattened prediction output vector is then reshaped to obtain the local prediction tensor of the node. (2) Define a containing The set of quantiles of elements, and select These are key quantile anchor points, each corresponding to the last dimension of the local prediction tensor; (3) Stack the local prediction tensors of all batches and all nodes in the spatial dimension, and swap the time and node dimensions to construct the final global probability prediction output tensor; (4) The confidence interval of any variable can be directly analyzed from the output tensor of the global probability prediction. The model is jointly trained and its parameters are optimized based on a global spatiotemporal integrated loss function. This simultaneously optimizes the prediction of data central tendency and distribution boundary within a unified computational graph, and evaluates uncertainty using quantitative indicators. This includes: (1) Define the quantile regression loss function, let For the true value, For probabilistic forecasting models targeting quantiles Calculate the single-point loss based on the predicted value. : ； (2) Aggregate the losses of all batch steps, nodes, times, variables and quantiles to construct a global spatiotemporal comprehensive loss function; (3) The gradient-based optimization algorithm is used to minimize the global spatiotemporal integrated loss function. During the training process, the gradient is calculated using the backpropagation algorithm, and all learnable parameters of the spatial layer, temporal layer and decoding layer are updated at the same time. (4) After the model training is completed, the test set data is input into the dual-flow graph state space model to generate the probability prediction results of the test set and verify its uncertainty quantification capability.

2. The method according to claim 1, characterized in that, Constructing a watershed hydraulic topology map based on discretely distributed water quality monitoring sections and their multidimensional monitoring data; including: defining all water quality monitoring sections within the target watershed as graph nodes, defining the natural channels connecting each section as graph edges, and determining the direction of the edges based on the water flow direction between graph nodes; The two-stream physical adjacency matrix includes: the downstream migration matrix and the upstream diffusion matrix. In the watershed hydraulic topology graph, if there is an edge from an upstream node to a downstream node, the elements of the downstream migration matrix are 1; otherwise, they are 0. The upstream diffusion matrix is ​​the transpose of the downstream migration matrix. Slices in the input spacetime tensor Representing the In the nth sample, the nth The time step, the first All standardized feature vectors of each node; slices in the truth label spacetime tensor Representing the In the sample, the future number The prediction time step, the first All standardized feature vectors of each node.

3. The method according to claim 1, characterized in that, The two-stream graph convolutional spatial fusion module utilizes the two-stream physical adjacency matrix to perform spatial-dimensional graph convolution operations on the input spatiotemporal tensor, extracting a high-dimensional spatial feature tensor containing river network spatial dependencies; including: (1) Using the normalized downstream migration matrix, perform downstream graph convolution on the features of each graph node to obtain the downstream feature matrix. : ;in, The learnable weight matrix of the forward convolution kernel. For the hidden layer dimension, For the bias term of the co-current convolution kernel, It is a non-linear activation function. For the normalized downstream migration matrix, For any sample and any time step The corresponding node feature matrix; (2) Using the normalized inverse diffusion matrix, perform inverse graph convolution on the features of each graph node to obtain the inverse feature matrix. : ;in, The independent, learnable weight matrix of the inverse convolution kernel. This is the bias term of the countercurrent convolution kernel. This is the normalized countercurrent diffusion matrix; (3) The upstream and downstream features are concatenated along the feature dimension, and node-level fusion weight coefficients are generated through a lightweight multilayer perceptron (MLP). ; Utilizing fusion weight coefficients The two sets of features are dynamically weighted and summed to obtain the fused spatial context feature matrix. (4) Repeat steps (1) to (3) above for all batches and all time steps of the input spatiotemporal tensor to finally obtain a high-dimensional spatial feature tensor that has fully extracted the spatial dependencies of the river network.

4. The method according to claim 1, characterized in that, The temporal evolution module based on the state-space model is used for deep temporal modeling of high-dimensional spatial feature tensors, outputting a full tensor of high-order spatiotemporal features; including: (1) Define a continuous linear time-invariant state equation. Decompose the high-dimensional spatial feature tensor output by the dual-flow graph convolutional spatial fusion module along the feature channel. Treat the time series of any node and any feature channel as a continuous dynamic system controlled by the differential equation. The input signal of the state equation is the value of a specific node and a specific feature channel of a certain sample in the high-dimensional spatial feature tensor at a certain time. The evolution of the corresponding latent state vector and output signal follows the continuous linear time-invariant state equation. (2) Diagonalize the state transition matrix in the continuous linear time-invariant state equation and initialize the diagonal elements of the state transition matrix with frequency tuning. (3) Discretize the continuous linear time-invariant state equations using the zero-order preservation principle to obtain discrete recursive equations; (4) Transform the discrete recursive equation into an equivalent global convolution form. Stack the input and output signals of all discrete time steps in time order to construct the input sequence vector and output sequence vector. Use Fast Fourier Transform and the convolution kernel generated by the convolution kernel generator to perform frequency domain convolution calculation on the input and output sequence vectors. Define the frequency domain convolution calculation process from the input sequence vector to the output sequence vector as... Operator; (5) Operators are embedded into deep neural network architectures to construct structures containing... A multi-layer stacked deep network, where the input tensor of the 0th layer is a high-dimensional feature tensor, and for the 1st layer of the network... layer, Its input is the output of the previous layer; extract the first... The layer outputs the state vector at the last time step, which forms the high-order spatiotemporal feature full tensor.

5. A system for quantifying the uncertainty of multi-step watershed water quality forecasting based on a dual-flow graph state-space model, characterized in that, include: The data processing unit is used to construct a watershed hydraulic topology map based on discretely distributed water quality monitoring sections and their multidimensional monitoring data. A dual-stream physical adjacency matrix is ​​constructed to describe the direct transmission effect of upstream water quality status on downstream, as well as the potential impact of downstream hydrological status on upstream. Multidimensional monitoring data is encapsulated into a spatiotemporal tensor containing four-dimensional information of batch, time, space, and features, including an input spatiotemporal tensor and a truth label spatiotemporal tensor. The model building unit is used to construct a dual-flow graph state-space model, including: a dual-flow graph convolutional spatial fusion module, a temporal evolution module based on the state-space model, a quantile regression decoder, and a global spatiotemporal comprehensive loss function. The dual-flow graph convolutional spatial fusion module uses the dual-flow physical adjacency matrix to perform spatial-dimensional graph convolution operations on the input spatiotemporal tensor to extract a high-dimensional spatial feature tensor containing river network spatial dependencies. The temporal evolution module based on the state-space model is used to perform deep temporal modeling on the high-dimensional spatial feature tensor and output a high-order spatiotemporal feature full tensor. The quantile regression decoder receives the high-order spatiotemporal feature full tensor, performs dimensional reshaping and multi-path decoding through a multilayer perceptron, and generates a global probability prediction output tensor containing multiple quantiles. The model training and optimization prediction unit is used to jointly train the model and optimize the parameters based on the global spatiotemporal integrated loss function. It simultaneously optimizes the central trend prediction and distribution boundary prediction of the data in a unified computation graph and evaluates the uncertainty quantification index.

6. An electronic device, characterized in that, include: Memory, used to store computer programs; A processor for executing the computer program to implement the method as described in any one of claims 1-4.

7. A non-volatile storage medium, characterized in that, Used to store a computer program, wherein the computer program, when executed by a processor, implements the method as described in any one of claims 1-4.

8. A computer program product comprising a computer program / instructions, characterized in that, When the computer program / instructions are executed by the processor, they implement the method described in any one of claims 1-4.

Images