Method for analyzing driving mechanism of reservoir group dispatching on non-consistent evolution of flood

Abstract

The application discloses a kind of reservoir group scheduling to the driving mechanism analysis method of non-consistent evolution of flood, comprising: using topological coupling physical manifold constraint denoising model to obtain smooth reservoir flood sequence, and the dynamic topological feature representing propagation time lag and sensitivity is constructed by graph neural ordinary differential equation based on hydraulic propagation mechanism;Structural causal model containing scheduling ability index, exogenous hydrology driving and topological feature is constructed, and nonlinear dependence relationship is fitted using causal generalized additive model;Perform counterfactual inference, calculate local causal driving index and global causal cumulative effect index by intervention operator;Based on causal index, the optimization scheduling strategy of inhibiting variation is generated.The application can realize the systematic analysis from data physics reduction to causal mechanism decoupling, accurately quantify the driving contribution of scheduling behavior to non-consistency of flood, and provide decision support for scientific flood control of river basin.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of hydrological analysis, and particularly relates to a driving mechanism analysis method for non-uniform evolution of reservoir group regulation on flood. BACKGROUND

[0002] Non-uniform evolution of flood is an important challenge currently faced by basin flood control safety and water resources management. Driven by global climate change and large-scale water conservancy construction, the flood process in the basin presents obvious non-stationary characteristics in the spatial and temporal distribution. Revealing the specific driving mechanism of reservoir group regulation behavior in flood variation has important scientific significance and engineering value for stripping natural variation and human intervention influence and formulating a scientific adaptive regulation strategy.

[0003] Existing researches mainly adopt statistical analysis or hydrological simulation methods for attribution analysis. Traditional methods are mostly based on Mann-Kendall test, coefficient of variation method to identify mutation points of flood sequences, or utilize soil and water resources assessment model SWAT, hydrological engineering center hydrological simulation system HEC-HMS and other hydrological models to restore natural runoff, and infer human activity influence by comparing the difference between the restored sequence and the measured sequence. In recent years, deep learning models based on long short-term memory network LSTM and graph convolution network GCN have also been widely applied to flood process simulation, trying to capture the evolution law of hydrological elements through data-driven methods.

[0004] However, the existing technology still has systematic defects in dealing with complex reservoir group driving mechanism, mainly manifested as lack of physical cause decoupling ability and dynamic structure capturing ability. Traditional statistical correlation analysis cannot distinguish the mixed effects of natural hydrological fluctuations and human regulation, and it is difficult to accurately strip out the independent causal contribution of the regulation strategy to the non-uniform evolution through counterfactual inference in mathematics; the existing data preprocessing methods lack physical manifold constraints, and the denoising process easily destroys the water balance and waveform characteristics, resulting in distorted basic data; the modeling method based on static graph structure ignores the physical propagation time delay and flow sensitivity evolution of water flow in the river network, and cannot accurately depict the dynamic time-varying characteristics of the topological relationship of the reservoir group. SUMMARY

[0005] The present application relates to the technical field of hydrological analysis, and particularly relates to a driving mechanism analysis method for non-uniform evolution of reservoir group regulation on flood.

[0006] Technical scheme, according to one aspect of the present application, a driving mechanism analysis method for non-uniform evolution of reservoir group regulation on flood, comprising:

[0007] obtaining a smooth inflow flood sequence denoised by topological coupling and physical manifold constraint, and a dynamic topological feature vector representing the upstream and downstream hydraulic propagation time delay relationship of the reservoir group;

[0008] a structural causal model containing the scheduling ability index sequence and the flood inconsistency index is constructed based on the smoothed reservoir inflow flood sequence, the dynamic topological feature vector, and the real-time collected exogenous hydrological driving data;

[0009] The variable relationship in the structural causal model is nonlinearly fitted to determine the causal dependence relationship of the scheduling ability index sequence on the flood inconsistency index;

[0010] Based on the causal dependence relationship, counterfactual inference is performed, the scheduling ability index sequence is forcibly set to a reference value through an intervention operator, and the conditional expectation of the flood inconsistency index under the counterfactual scenario is calculated;

[0011] Based on the difference between the conditional expectation and the observed expectation, a local causal driving index representing the instantaneous impact of reservoir scheduling on flood variation is calculated, and the local causal driving index is integrated in the time dimension to obtain a global causal cumulative effect index representing the long-term cumulative contribution;

[0012] An optimal scheduling strategy for inhibiting the evolution of flood inconsistency is generated based on the global causal cumulative effect index.

[0013] According to another aspect of the present application, a system for analyzing the driving mechanism of reservoir group scheduling on the evolution of flood inconsistency includes:

[0014] a memory for storing a computer program;

[0015] a processor for executing the computer program to implement the steps of the method for analyzing the driving mechanism of reservoir group scheduling on the evolution of flood inconsistency according to any one of the above.

[0016] The above technical solutions can realize systematic analysis from data physical restoration to causal mechanism decoupling, accurately quantify the driving contribution of scheduling behavior to flood inconsistency, and provide decision support for scientific flood control in a basin. BRIEF DESCRIPTION OF DRAWINGS

[0017] Figure 1 The figure is a schematic diagram of the overall flow of the method for analyzing the driving mechanism of reservoir group scheduling on the evolution of flood inconsistency.

[0018] Figure 2 The figure is a schematic diagram of the process for obtaining a smoothed reservoir inflow flood sequence constrained by topologically coupled physical flow.

[0019] Figure 3 The figure is a schematic diagram of the process for obtaining a dynamic topological feature vector representing the hydraulic propagation time lag relationship between the upstream and downstream of the reservoir group.

[0020] Figure 4 The figure is a schematic diagram of the parallel subsystem coordination architecture for constructing a deep reinforcement learning model. Detailed Implementation

[0021] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0022] Example 1 details the overall process of the analytical method for the driving mechanism of reservoir group scheduling on the non-uniform evolution of floods, such as... Figure 1 As shown, by integrating physical manifold denoising, structural causal inference, and counterfactual analysis techniques, a complete closed loop from data cleaning to mechanism analysis and strategy generation is achieved, solving the technical problem that traditional correlation analysis cannot separate the causal contributions of natural variation and human scheduling.

[0023] Step 101: Obtain the smooth inflow flood sequence after denoising by physical manifold constraints with topological coupling, and the dynamic topological feature vector characterizing the time delay relationship of hydraulic propagation between the upstream and downstream reservoir groups.

[0024] In this step, the raw hydrological data must first be acquired and preprocessed. A smoothed inflow flood sequence refers to a reservoir inflow time series that has undergone specific denoising processing to remove high-frequency observation noise and conforms to the laws of physical conservation. Specifically, the system collects raw operational records of the reservoir group in real time or reads them from a historical database. These raw operational records typically contain information such as water level and outflow, and often contain non-physical fluctuations caused by observation errors or environmental interference, i.e., sawtooth noise. To obtain a high-quality smooth sequence, this embodiment employs a topologically coupled physical manifold constraint denoising method, which will be described in detail in subsequent embodiments. This method utilizes a conditionally invertible normalized flow model to map the data to a latent space and combines this with water conservation constraints for diffusion denoising, resulting in a flow sequence that is both smooth and conforms to hydraulic laws.

[0025] The dynamic topological feature vector is a characteristic representation used to quantitatively describe the complex hydraulic connection within the reservoir group. In a large river basin, the outflow of the upstream reservoir propagates to the downstream reservoir through the river channel, and there is a significant time delay effect in this process, and the mutual influence strength between each reservoir dynamically changes with the size of the flow. In this step, a graph neural differential equation model based on the hydraulic propagation mechanism is constructed to model this dynamic relationship. Specifically, the system uses the basic attribute data of the reservoirs (such as geographical location, reservoir capacity curve) and the above smoothed inflow flood sequence to solve the topological evolution equation by numerical integration, generating an evolution trajectory containing time-varying node states and time-varying edge weights. The evolution trajectory constitutes the dynamic topological feature vector. This feature vector encodes the static river network structure, captures the time delay characteristics and flow sensitivity of water flow propagation, and provides a structural basis for subsequent causal analysis.

[0026] In step 102, based on the smoothed inflow flood sequence, the dynamic topological feature vector, and the real-time collected exogenous hydrological driving data, a structural causal model containing a sequence of scheduling capability indices and a flood inconsistency index is constructed.

[0027] On this basis, this step establishes a mathematical model describing the causal dependence relationship between key variables. Exogenous hydrological driving data refers to natural meteorological and hydrological factors that are not controlled by reservoir scheduling behavior, such as natural rainfall intensity in the basin, inter-branch inflow, etc. It is usually obtained through meteorological satellites, rain gauge networks, or hydrological simulators to represent the natural background of flood evolution. The sequence of scheduling capability indices is a comprehensive physical index that quantifies the reservoir's potential for regulation and storage at different times. Specifically, this index can be calculated based on the effective flood control capacity of the reservoir, the current occupied reservoir capacity, and the safety management factor of the scheduling strategy. For example, the calculation formula of the improved reservoir scheduling index MRI can be expressed as:

[0028] M i,t =(V total,i -V current,i (t))*(1-S j ) / W A ;

[0029] Where V total,i is the total flood control capacity of the i-th reservoir, V current,i (t) is the used reservoir capacity at time t, S j is a safety management factor reflecting the conservatism of the scheduling strategy, with a value range of 0 to 1, and W A is the reference total flood volume. The higher the index, the stronger the reservoir's ability to regulate floods.

[0030] The flood inconsistency index is a statistical quantity used to characterize the degree of variation of the flood waveform relative to the natural state. Specifically, the index can include the relative deviation of the flood peak flow, the deviation of the total flood volume, and the dynamic time warping (DTW) distance between flood hydrographs. For example, the inconsistency index Y t can be defined as a vector:

[0031] Y t = [ΔQ peak , D DTW ] T ;

[0032] where ΔQ peak is the difference ratio of the flood peak after regulation and the natural flood peak, and D DTW is the shape distance of the two hydrographs. The structural causal model is a directed acyclic graph whose nodes include the above-mentioned regulation capacity index sequence, exogenous hydrological driving data, dynamic topological feature vector, and flood inconsistency index. In this model, the regulation capacity, exogenous driving, and topological features are regarded as causal parent nodes, which collectively point to the flood inconsistency index as a child node, explicitly establishing the driving logic relationship and avoiding the risk of causal inversion in traditional correlation analysis.

[0033] In step 103, the nonlinear fitting of the variable relationship in the structural causal model is performed to determine the causal dependence relationship of the regulation capacity index sequence on the flood inconsistency index.

[0034] After the graph structure is constructed, it is necessary to quantitatively describe the specific influence function of the parent nodes on the child nodes. Since the impact of reservoir regulation on flood waveforms usually exhibits complex nonlinear characteristics, for example, when the reservoir capacity is sufficient, the regulation impact is small, while near full reservoir, the impact increases dramatically, and simple linear regression cannot accurately depict it. This embodiment uses the causal generalized additive model (cGAM) to fit the nonlinear dependence relationship.

[0035] Specifically, it is assumed that the flood inconsistency index follows a certain probability distribution, such as a normal distribution or a gamma distribution, and the location parameter (such as the mean), the scale parameter (such as the variance), and the shape parameter of this distribution are all modeled as the sum of smooth functions of each parent node. For example, the location parameter μ can be represented as:

[0036] μ = f M (M t ) + f R (R t ) + f T (T t );

[0037] where f M , f R , f TThe spline smoothing functions are respectively for scheduling ability, exogenous driving and topological characteristics. Through maximum likelihood estimation or Bayesian inference method, the specific form of each smoothing function can be determined by training the model with historical data, and the causal dependence relationship can be obtained.

[0038] In step 104, counterfactual inference is performed based on the causal dependence relationship, the scheduling ability index sequence is forcibly set to the reference value by the intervention operator, and the conditional expectation of the flood inconsistency index in the counterfactual scenario is calculated.

[0039] This step answers the counterfactual question of how the flood variation will change if the scheduling ability changes. The intervention operator is the do operator in causal inference, denoted as do(M t =m). Its mathematical meaning is to cut off all incoming edges of the node M t pointing to the scheduling ability in the structural causal graph, and forcibly fix its value to a specific reference value m, while keeping the distribution of other parent nodes (exogenous driving R t and topological characteristics T t unchanged. This operation is mathematically different from the traditional conditional probability calculation P(Y|M=m), which blocks the backdoor path of confounding factors and can purely separate the causal effect of scheduling behavior. In specific implementation, the system substitutes the reference value m into the causal generalized additive model fitted in step 103, keeps the values of f R (R t ) and f T (T t ) unchanged, and calculates the expected value of the model output, i.e. the conditional expectation E[Y t |do(M t =m)]. This expected value represents the theoretical value of the flood inconsistency index under the assumed scheduling scenario.

[0040] In step 105, based on the difference between the conditional expectation and the observed expectation, the local causal driving index representing the instantaneous influence of reservoir scheduling on flood variation is calculated, and the global causal cumulative effect index representing the long-term cumulative contribution is obtained by integrating the local causal driving index in the time dimension, and the optimal scheduling strategy for inhibiting the evolution of non-consistency is generated based on the global causal cumulative effect index.

[0041] By quantifying the changes before and after counterfactual intervention, the specific driving index is obtained. The local causal driving index LCE is used to measure the marginal contribution of a small change in scheduling ability to flood inconsistency at a specific time. Specifically, it can be calculated by finite difference method:

[0042] LCE t =(E[Y t |do(M t =m+Δm)]-E[Y t |do(Mt = m + (m - m)) / Am;

[0043] where Am is a small perturbation. The larger the absolute value of the index, the higher the sensitivity of the scheduling behavior at that moment to the flood variation. The global causal cumulative effect index GCE is the integral of the local index over the entire flood process or a specific scheduling period. Specifically, the trapezoidal rule can be used to numerically integrate the LCE t sequence:

[0044] GCE = ∑(LCE t + LCE t+1 )*Δt / 2. The index reflects the overall contribution of reservoir scheduling to the non-uniform evolution of the flood in the entire event. Through the two indexes, it can be clearly identified which period and which reservoir is the dominant factor leading to the variation of the flood waveform, providing a quantitative basis for subsequent scheduling optimization.

[0045] Embodiment two, the specific process of obtaining the smooth incoming flood sequence is described in detail, as Figure 2 shown, a denoising method based on topological coupling conditional reversible normalization flow and diffusion model is proposed to solve the common sawtooth noise and water imbalance problem in actual reservoir operation data. By applying physical constraints in the latent space, high-fidelity data restoration is achieved.

[0046] Step 201, collect the historical operation records of the reservoir group, and based on the water balance principle, obtain the sawtooth incoming flood sequence containing non-physical fluctuations.

[0047] In this step, first, the basic data is backstepped. The incoming flow of the reservoir is usually difficult to measure directly, but is calculated indirectly through the change of reservoir capacity and outgoing flow. The system collects the historical water level Z t and outgoing flow Q out,t sequence of the reservoir. Using the water level-storage capacity curve V = f(Z) of the reservoir, the water level sequence is converted to the storage capacity sequence V t . According to the discrete time water balance equation: Q in,t = (V t+1 - V t ) / Δt + Q out,t , the incoming flow sequence Q in,t is calculated. Since the calculation of storage capacity V t is sensitive to water level observation errors, and Δt is usually small (such as 1 hour), small water level fluctuations will be amplified by the difference operation, resulting in the calculated incoming flow sequence Q in,t showing sharp, sawtooth fluctuations that do not conform to hydrological laws. Therefore, this sequence is called the sawtooth incoming flood sequence and needs further denoising.

[0048] Step 202, constructing a conditional invertible normalizing flow model conditioned on the dynamic topological feature vector, and using the conditional invertible normalizing flow model to establish a bidirectional reversible mapping between the jagged reservoir inflow flood sequence and the Gaussian latent variable in the potential space.

[0049] To effectively remove complex non-Gaussian noise, the embodiment introduces a generative model technique. The conditional invertible normalizing flow model is a deep neural network composed of a series of reversible transformation layers, which can map data of complex distribution to a simple known distribution, such as standard normal distribution.

[0050] Specifically, let the jagged sequence be a vector q, and the dynamic topological feature vector be a condition τ, the model defines a reversible mapping z = T θ (q; τ), where T θ is a conditional invertible transformation function parameterized by parameters θ. By maximizing the log-likelihood function, the model learns to transform the flood sequence q containing noise into a variable z in the latent space. Ideally, if q is a pure flood sequence, z should obey the standard Gaussian distribution; and due to the inclusion of noise in q, the distribution of z will deviate. This mapping not only preserves all the information of the data (due to reversibility), but also encodes complex physical and topological constraints into the transformation function T θ .

[0051] Step 203, constructing a denoising diffusion probability model in the latent space, performing forward noise addition on the Gaussian latent variable to learn the noise distribution, and removing noise by reverse diffusion sampling to obtain a denoised latent manifold representation.

[0052] After mapping the data to the latent space, a diffusion model is used for denoising. Unlike direct denoising in the original data space, diffusion in the latent space can significantly reduce computational complexity and capture higher-level semantic features. Specifically, in the forward process, the system gradually adds Gaussian noise to the latent variable z until it becomes pure noise, learning the noise distribution pattern. In the reverse process, a denoising network is trained to predict and remove the noise added at each step. The reverse sampling process is not blind, but is guided by physical manifold regularization (see Step 205 for details). Through multi-step reverse sampling, the system recovers the pure latent variable from pure noise, i.e., the denoised latent manifold representation z clean . clean represents the optimal representation in the latent space that conforms to the laws of flood physics.

[0053] Step 204, using the inverse transformation of the conditional invertible normalizing flow model to map the denoised latent manifold representation back to the physical data space, obtaining a smooth reservoir inflow flood sequence that conforms to the conservation of hydraulics constraints.

[0054] Utilizing the invertible property of the normalized flow model, a smoothed inflow flood sequence q is obtained through inverse transformation. clean :

[0055] q clean =T θ -1 (z clean ;τ), mapping the denoised latent variables back to the original physical data space. Because T θ Topological features have already been fused during training, and z clean The generated smooth inflow flood sequence q is denoised using physical constraints. clean Not only is it smooth and even, eliminating sawtooth-like fluctuations, but it also strictly conforms to the hydraulic laws of the reservoir group in terms of overall trend and waveform characteristics.

[0056] Step 205, the training process of the conditionally invertible normalized flow model and the denoised diffusion probability model includes physical manifold regularization constraints, which include:

[0057] The formula for calculating the water conservation constraint is as follows:

[0058] ;

[0059] Among them, L mass Let λ be the loss value of the water conservation constraint term. mass Here, K represents the weighting coefficient for the water conservation constraint, L represents the total number of flood events, and q represents the total number of time steps for a single flood event. l k Let q be the inflow rate of the k-th flood at the l-th time step generated by the model. out,l k Let be the outflow from the reservoir at time step l for the k-th flood.

[0060] The cross-event manifold regularization term is calculated as follows:

[0061] ;

[0062] Among them, L manifold The loss value for the cross-event manifold regularization term, z k and z k' Let |k| and |k'| be the latent manifold representation vectors of the k-th and k'-th floods in the latent space, respectively. 2 w represents the square of the Euclidean distance. k,k' The similarity weight between the events is calculated based on the distance between the exogenous hydrological driving data of the two floods.

[0063] a water mass conservation constraint term to penalize the deviation between the time integral of the generated flood hydrograph and the theoretical total mass calculated from the outflow and reservoir storage variation; an inter-event manifold regularization term to compute event weights based on the similarity of the exogenous hydrological driving data of different flood events, and to penalize the Euclidean distance between the latent manifold representations of similar flood events in the latent space according to the event weights.

[0064] To ensure the physical plausibility of the generated results, the loss function terms in the training process are defined in detail in this step. The water mass conservation constraint term L mass is used to enforce the overall volume of the generated flood hydrograph. According to the principle of reservoir operation, the reservoir storage variation is usually zero or a constant value within a complete flood event period, so the total inflow mass should be approximately equal to the total outflow mass (ignoring minor factors such as evaporation and seepage). This constraint term enforces the model to generate q clean that is consistent with the water mass conservation principle by penalizing the deviation between the integral of the inflow and the integral of the outflow. mass is a positive hyperparameter used to adjust the strength of this constraint.

[0065] The inter-event manifold regularization term L manifold utilizes the similarity between different flood events. In hydrology, floods with similar rainfall causes and basin backgrounds should be close to each other in the latent space, forming a low-dimensional physical manifold. This constraint term calculates the Euclidean distance between the latent space representations z k and z k' of any two floods k and k', and uses the weight w k,k' for weighting. The weight w k,k' reflects the similarity of the external driving conditions, which can be defined as:

[0066] w k,k' = exp(-dist(R k ,R k' ) / σ);

[0067] where dist(R k ,R k' ) represents the distance measure between the exogenous hydrological driving data R k and R k' , and σ is a scale parameter or bandwidth parameter that controls the speed of similarity decay. This formula indicates that the more similar the exogenous driving is, the greater the weight is. By minimizing this loss, the model is encouraged to map physically similar floods to close positions in the latent space, making the generated latent manifold more continuous and compact, and improving the model's generalization ability and denoising stability in the case of scarce samples or extreme flood scenarios.

[0068] In addition to the above scheme, after obtaining the smoothed reservoir inflow sequence, the natural flood can be further reconstructed in combination with the interval inflow. Specifically, the smoothed upstream reservoir inflow is calculated to the downstream section using the Muskingum flood routing algorithm, and the inflow of the interval tributaries is superimposed to obtain the natural flood hydrograph without the influence of the reservoir. This step provides baseline data for subsequent calculation of non-consistency indicators (such as comparison of regulated flood and natural flood).

[0069] The embodiment effectively solves the noise problem in reservoir inflow backstepping by combining reversible flow model, diffusion model and explicit physical regularization constraint, and provides high-quality data basis for subsequent causal analysis.

[0070] Embodiment three, the specific implementation process of obtaining dynamic topological feature vector is described in detail, as shown in Figure 3 The introduction of explicit physical propagation mechanism (time delay and attenuation) and flow sensitivity dynamic weight solves the technical problem that traditional static graph neural network cannot depict the time-varying hydraulic connection of reservoir group.

[0071] Before entering the ordinary differential equation modeling, as a preferred embodiment, the multi-source heterogeneous data is first preprocessed by using the graph convolution network. Specifically, the system collects the static attributes (such as geographic location, design reservoir capacity) and dynamic attributes (such as water level, flow sequence) of the reservoir, and constructs the initial graph structure. The node features are aggregated by using multi-layer graph convolution layer, and the formula can be expressed as:

[0072] ;

[0073] Where σ'(.) is a nonlinear activation function such as ReLU, A is an adjacency matrix, D is a degree matrix, H is a feature matrix, H l is the node feature matrix of the lth layer, H l+1 is the updated node feature matrix of the l+1th layer, W l is the learnable weight matrix of the lth layer.

[0074] The node embedding vector obtained after preprocessing is used as the initial state x(0) of the subsequent differential equation, which provides the starting information containing spatial correlation for the evolution process.

[0075] Step 301, based on the smoothed reservoir inflow sequence and the basic attribute data of the reservoir group, a graph neural differential equation model is constructed, and the reservoir is defined as a graph node, and the river connection between the reservoirs is defined as a graph edge.

[0076] In this step, the system maps the reservoir group in the physical world to the topological structure in the graph neural differential equation. The graph node not only represents the reservoir entity, but also its state vector x iAlso carries the real-time physical state of the reservoir. The edge represents the entity river, whose weight w ji Not a fixed connection relationship, but represents the real-time hydraulic influence strength of the upstream reservoir j on the downstream reservoir i. Unlike traditional discrete-time recurrent neural networks, the ordinary differential equation model describes the evolution of the state as a continuous-time flow field, which is highly consistent with the physical nature of the continuous propagation of water flow in the river channel, and can more accurately handle non-uniformly sampled data and long-short term dependencies.

[0077] Step 302, establish a topological evolution equation describing the evolution of the state of the reservoir node with time, the topological evolution equation includes a local dynamics term describing the operation characteristics of a single reservoir, and a propagation coupling term describing the coupling influence of upstream nodes and parallel nodes on the current node.

[0078] This step establishes the dynamics of the system. The topological evolution equation is specifically represented as:

[0079] ;

[0080] The left side of the equation dx i / dt represents the state vector x i of the i-th reservoir node with respect to time t, i.e. the rate of change of the state.

[0081] The first term f local on the right side of the equation is a local dynamics function describing the operation characteristics of a single reservoir, which describes the self-evolution law of the reservoir without external inflow interference, such as the change of reservoir capacity with discharge. The local dynamics term includes a reservoir capacity change function based on the water balance principle of the reservoir, and a flow regulation function based on the reservoir water level-storage curve and discharge capacity curve.

[0082] The second term is the upstream propagation coupling term, which includes an upstream propagation function that explicitly introduces a water flow propagation time delay parameter, which is used to map the outflow of the upstream reservoir node at the lag time to the driving contribution of the inflow of the current reservoir node, U i is the set of upstream reservoirs of the i-th reservoir node, w ji (t) is the edge weight from upstream node j to current node i at time t, g up is a coupling function describing the upstream propagation influence, τ ji is the water flow propagation time delay from node j to node i.

[0083] The third term is the parallel coupling term, P i is the set of reservoirs parallel to the i-th reservoir node, i.e. the reservoirs at the same basin level, w ki (t) is the edge weight from parallel node k to current node i, g para is a function describing the parallel coupling influence.

[0084] Through the above decomposition, the model clearly separates self-evolution from interaction.

[0085] Step 303, the state vector is defined as:

[0086] x i =[V i ,Q in,i ,Q out,i ,Z i ] T ;

[0087] where V i , Q in,i , Q out,i , Z i are the storage, inflow, outflow, and water level of the i-th reservoir, respectively.

[0088] The local dynamics function f local is specifically represented as:

[0089] ;

[0090] where κ1 is the inflow recovery coefficient, Q in,i eq is the equilibrium inflow, φ(V i , Z i ) is the theoretical discharge capacity function based on reservoir characteristics, and ψ(V i ) is the water level calculation function based on the storage curve.

[0091] The coupling function g up describing the upstream propagation effect is specifically represented as:

[0092] ;

[0093] where η ji is the propagation attenuation coefficient from the upstream node j to the current node i, and Q out,j (t-τ ji ) represents the outflow of the upstream node j at the time lagged by the time delay τ ji .

[0094] This step finely defines the physical equation, ensuring the physical interpretability of the model. For the local dynamics function f local , the first component Q in,i -Q out,i directly corresponds to the water balance equation in physics, i.e., the change rate of storage is equal to the inflow minus the outflow. The second component introduces the recovery coefficient κ1 to simulate the tendency of the inflow to the equilibrium value Q in,i eqregression process of the base flow. The third and fourth components constrain the outflow and water level to satisfy the physical characteristics of hydraulic structures, φ(V i , i ) is usually determined by the orifice flow equation Q = C * A * (2gh) 0.5 , ψ(V i ) is determined by the measured water level-storage curve. For the upstream propagation function g up , the focus is on the second component η ji *Q out,j (t-τ ji ). This explicitly introduces a time delay τ ji , the outflow of the upstream reservoir at time t-τ ji , propagates and attenuates (described by the coefficient η ji , 0 < η ji < 1) in the river channel to form the inflow increment of the downstream reservoir at the current time t. Compared with the traditional synchronous correlation modeling, the above time-lag modeling method can more truly reflect the physical evolution process of the flood wave.

[0095] Step 304, the topological evolution equation further includes an edge weight evolution equation, which is defined as:

[0096] ;

[0097] wherein dw ij / dt represents the derivative of the edge weight w ij with respect to time t, λ w is a weight recovery coefficient, w ij base is a basic edge weight determined based on the river network distance, μ w is a sensitivity response coefficient, σ(.) is a Sigmoid activation function, is the partial derivative of the inflow of the current node i with respect to the outflow of the upstream node j, used to represent the real-time sensitivity of the flow change.

[0098] The edge weight evolution equation defines the functional relationship between the change rate of the graph edge weight and the flow sensitivity response; the flow sensitivity response is determined by calculating the partial derivative of the inflow of the current reservoir node with respect to the outflow of the upstream reservoir node, and is used to dynamically enhance the weight of the graph edge with strong hydraulic connection.

[0099] This step further defines the dynamics of the network structure. The edge weight is no longer a static value, but changes in real time with the closeness of the hydraulic connection. The first term of the equation -λ w *(w ij -w ij base) is the regression term, which ensures the weights tend to the base weights w determined by geographical distance in the absence of perturbations ij base , maintaining the structural stability of the topology. The second term of the equation is the driving term, where the partial derivative measures the sensitivity of downstream inflow to upstream outflow. During flood propagation, this partial derivative will increase significantly when the flood peak arrives, indicating that the control effect of the upstream on the downstream is enhanced. Through the Sigmoid function σ, the sensitivity is mapped to the interval (0, 1) and multiplied by the response coefficient μ w , so that the edge weight w ij increases accordingly. The above mechanism enables the model to automatically focus on the upstream reservoir that plays a dominant role at the current time, dynamically reconstructing the topology structure.

[0100] Step 305, use a numerical integrator to solve the topology evolution equation to obtain the evolution trajectory containing the time-varying node state vector and time-varying edge weight matrix, and take the evolution trajectory as the dynamic topology feature vector.

[0101] Since the above equation usually cannot obtain an analytical solution, the numerical integration method is used for solving. Preferably, the fourth-order Runge-Kutta (RK4) solver or the Dormand-Prince solver with adaptive step size (also known as the Dormand-Prince method in Chinese) is used. Starting from the initial time t = 0, gradually integrating according to the time step h, the node state x(t) and the edge weight w(t) at all times within the entire flood regulation period can be obtained. This complete spatiotemporal evolution trajectory is extracted, flattened or encoded to form a dynamic topology feature vector T t , which is the key input of the subsequent causal analysis module, providing the physical structure information behind the flood evolution.

[0102] Embodiment four, the specific implementation of structural causal analysis and counterfactual inference is elaborated in detail, by introducing the improved reservoir index MRI and the quantitative definition of non-uniformity index, combining the cGAM model and the do operator, the accurate quantification of the reservoir dispatching driving mechanism is realized.

[0103] Step 401, the nonlinear dependence relationship between variables in the structural causal model is described by the causal generalized additive model, where the flood non-uniformity index Y

[0104] In this step, the input and output variables of the model need to be first clarified. The flood non-uniformity index Y tis calculated based on the post-regulation flood and the natural flood. For example, the dynamic time warping (DTW) distance is used to measure the waveform difference, combined with the peak value deviation, defined as follows:

[0105] Y t =α*D DTW +β*|Q peak,adj -Q peak,nat | / Q peak,nat ;

[0106] where D DTW is the dynamic time warping distance, Q peak,adj is the post-regulation flood peak flow, Q peak,nat is the natural flood peak flow, and α and β are weighting coefficients.

[0107] The scheduling capacity index sequence M i,t needs to be calculated explicitly. In this embodiment, the following formula is preferably used:

[0108] ;

[0109] where H j R and H A are the catchment areas of the reservoir j and the common section A, V total,j -V curr,j (t) is the available flood control storage, S j is the safety management factor, 0 represents aggressive scheduling, and 1 represents conservative scheduling, and W A is the total flood volume of the section. The physical meaning of the formula is that the greater the remaining flood control storage, the higher the catchment area ratio, and the more aggressive the scheduling strategy, the stronger the reservoir's reshaping ability of the flood. Based on the above variables, a causal generalized additive model (cGAM) is constructed, and it is assumed that Y t obeys the distribution D(μ,σ,ν), and the distribution parameters are modeled as functions of the parent nodes.

[0110] The functional expression of the position parameter is as follows:

[0111] ;

[0112] where μ j,t is the position parameter of the jth flood non-uniformity index at time t, α j,0 is a constant intercept term, N is the total number of reservoirs, M i,t is the scheduling capacity index of the ith reservoir at time t, s μ,j i (.) is a spline smoothing function for the scheduling capacity index, R t is exogenous hydrological driving data, s μ,j R(.) is a spline smoothing function for exogenous drivers, T t is the dynamic topological feature vector, s μ,j T (.) is a spline smoothing function for topological features.

[0113] This step describes the important construction of cGAM in detail. The location parameter μ usually corresponds to the conditional expectation E[Y]. The s(.) in the formula is preferably expanded using B-spline or thin-plate regression spline basis functions. For example:

[0114] ;

[0115] where b k is the basis function, and β k is the coefficient to be estimated. The additive structure lies in its interpretability: assuming that the effects of each factor are additive, the contribution of the scheduling ability can be isolated from the total change. Similarly, the scale parameter σ and the shape parameter ν can also be modeled in an additive form to capture the changes in heteroscedasticity and higher moments.

[0116] Step 403, the local causal driving index is defined as the sensitivity or partial derivative of the conditional expectation of the flood inconsistency indicator to the change in the scheduling ability index in the counterfactual scenario when the single-reservoir scheduling ability index is forcibly changed; the global causal cumulative effect index is defined as the time integral of the local causal driving index within a preset analysis period, used to represent the total amount of cumulative contribution of reservoir scheduling behavior to the evolution process of flood inconsistency.

[0117] This step defines the physical meaning of the two important quantitative indicators. The local causal driving index LCE is a counterfactual sensitivity indicator, which does not measure the correlation in the observed data, but answers how much the inconsistency indicator Y is expected to change if the scheduling ability of reservoir i is artificially increased by one unit at time t. This definition excludes the confounding effects of natural rainfall R and topological structure T. The global causal cumulative effect index GCE is the time accumulation of LCE, which answers how much inconsistency change is caused by the total amount of scheduling behavior of reservoir i in the entire flood season or a certain flood process.

[0118] Step 404, the calculation formula of the local causal driving index is:

[0119] ;

[0120] where LCE i,j,t is the local causal driving index of the i-th reservoir to the j-th flood inconsistency indicator at time t, denotes the partial derivative operator with respect to the scheduling ability index, E[·|do(·)] denotes the counterfactual conditional expectation after performing the intervention operator, M 0i,t is the reference value under the current actual dispatching scheme;

[0121] The calculation formula of the global causal cumulative effect index is:

[0122]

[0123] wherein, GCE i,j is the global causal cumulative effect index of the i th reservoir on the j th flood non-uniformity index, t start and t end are the starting and ending time of the analysis period respectively, and Δt is the discrete time step.

[0124] This step gives the specific calculation and numerical implementation method of the index. For the calculation of LCE, since the cGAM model has been fitted, the counterfactual expectation E[Y|do(M=m)] is equivalent to fixing R t and T t are the observed values, and the output obtained by replacing M input with m. The partial derivative can be calculated by analytical method (derivation on spline basis function) or finite difference method. Preferably, the central difference method is adopted:

[0125] LCE≈(E[Y|do(M=m+ε)]-E[Y|do(M=m-ε)]) / (2ε);

[0126] wherein, ε is a small increment. For the calculation of GCE, the trapezoidal integral formula is used to sum the discrete LCE time series. By calculating the GCE values of all reservoirs and sorting, the control reservoirs that play a leading role in the evolution of flood non-uniformity can be identified.

[0127] Example five, the specific implementation of the causal feedback reinforcement learning dispatch is elaborated in detail, the LCE index obtained by the foregoing causal analysis is converted into a control signal, and the intelligent dispatching that takes into account flood control safety and non-uniformity control is realized by improving the reward function and the policy network architecture.

[0128] Step 501, generating an optimized dispatching strategy for suppressing the evolution of non-uniformity includes: constructing a deep reinforcement learning model, and defining an immediate reward function containing multiple constraint targets; the immediate reward function is composed of three parts: a flood control safety penalty term for punishing the cross-section flow exceeding the limit, a non-uniformity level penalty term for punishing the variation degree of flood wave shape, and a causal driving intensity penalty term for punishing the strong causal contribution of the dispatching behavior to the variation.

[0129] In this step, an agent is constructed, whose state space S t includes the reservoir water level, inflow, LCE index and time coding at the current time, and the action space A​t To make the discharge decision for each reservoir. The traditional scheduling only focuses on flood control safety, which can easily lead to excessive damage to the natural flood wave shape. This embodiment introduces a causal driving strength penalty term. The principle is: if a certain scheduling operation leads to a sharp increase in LCE value, i.e. the operation strongly drives the inconsistency evolution, a larger negative reward is given to force the agent to explore a mild scheduling strategy that has less impact on inconsistency, and to suppress the generation of variation at the source.

[0130] Step 502, the calculation formula of the instant reward function is:

[0131] ;

[0132] Wherein, R t is the instant reward value at time t, λ1, λ2, λ3 are non-negative weight coefficients of flood control safety, inconsistency level and causal driving strength respectively, Q max,t is the maximum discharge of the section, Q safe is the safety discharge threshold, NCI t is the comprehensive inconsistency level index, N is the number of reservoirs, J is the number of inconsistency indexes, |LCE i,j,t | is the absolute value of the local causal driving index.

[0133] This step quantifies the multi-objective trade-off mechanism. The first term -λ1*Q max,t / Q safe is a softened form of hard constraint, which ensures that the section discharge does not exceed the safety threshold and guarantees the flood control bottom line. The second term -λ2*NCI t focuses on the results and punishes the deviation of the final flood wave shape from the natural wave shape. The third term -λ3*∑|LCE| focuses on the process and mechanism. By punishing the absolute value of the local causal driving index LCE, the model will tend to choose scheduling actions with LCE close to zero, i.e. choose strategies that statistically mimic nature and contribute the least to variation. The setting of weight coefficient λ usually follows the principle of λ1>>λ3>λ2, i.e. priority is given to safety, further inhibition of causal driving, and finally consideration of wave shape similarity.

[0134] As a preferred embodiment, this embodiment improves the structure of the policy network of reinforcement learning and adopts a structure-improved Transformer. In view of the long dependence characteristics of the flood sequence, a dynamic position perception mechanism with gating adjustment is introduced in the position encoding of the encoder:

[0135] PE(pos,t)=sin(pos / 10000 2i' / d )+g t *cos(t / 10000 2i' / d );

[0136] where PE(pos, t) is the dynamic position encoding, pos is the sequence position index, t is the physical timestamp indicating the real-world time corresponding to the time step, i' is the dimension index of the position encoding vector, d is the model dimension, g t is the adaptive coefficient generated by the gating network, used to dynamically adjust the weight of the time information according to the water regime change.

[0137] In addition, at the decoder end, the quantum-inspired encoding technology is combined to use the Hadamard quantum gate to superimpose and entangle the feature vectors, enhancing the network's high-dimensional expression ability for the coupled state space of multiple reservoirs. This architecture design improves the decision stability and convergence speed of the agent in a complex hydraulic environment.

[0138] Embodiment six, another preferred embodiment or alternative for obtaining a smooth reservoir inflow sequence is provided.

[0139] Unlike the method based on latent space manifold mapping, this embodiment adopts a hierarchical constraint diffusion model based on hydrological physical consistency. This method directly operates on the reservoir inflow sequence in the data space, constructs a three-layer constraint system containing global constraints, local constraints and time sequence constraints, and combines a dynamic weight adjustment mechanism, without complex manifold mapping calculations, to achieve high-quality denoising that meets physical conservation. It is more suitable for application scenarios where computing resources are limited or real-time requirements are high.

[0140] A multi-scale physical constraint system for the flood process is constructed, including three levels.

[0141] The first layer is the global water balance constraint, which is used to ensure the total water balance of the entire flood process. Its calculation formula can be expressed as:

[0142] ;

[0143] where C global is the global water imbalance, Q in,t and Q out,t are the inflow and outflow at the t period, V start and V end are the initial and final reservoir capacities, ε g is the global tolerance threshold, and Δt is the discrete time step, used to convert the flow to volume.

[0144] The second layer is the local water balance constraint, which divides the entire period into K sub-periods, and establishes a constraint for each sub-period k:

[0145] ;

[0146] where C local kΩkis the water imbalance for the kth local time window (or sub-process) k ΔV is the set of time for the sub-period k ε is the change of reservoir capacity in this period l ε is the change of reservoir capacity in this period

[0147] The third layer is the temporal continuity constraint, which limits the rate of change of flow between adjacent periods based on hydrological prior knowledge. Specifically, for the rising period, we require:

[0148] |Q in,t+1 -Q in,t |≤α rise *Q in,t ;

[0149] For the falling period, we require:

[0150] Q in,t+1 ≥Q in,t *exp(-β recession );

[0151] where α rise is the maximum increasing coefficient and β recession is the falling constant. This layer of constraint effectively eliminates the non-physical abrupt jump.

[0152] Inverse diffusion sampling strategy with dynamically adjusted constraint strength. In the inverse denoising process of the diffusion model, the noise level gradually decreases with the diffusion step s. This embodiment introduces a constraint strength adjustment function ω(s) related to the diffusion step. The function is defined as:

[0153] ω(s)=ω min +(ω max -ω min )*(1-s / S) γ ;

[0154] where s decreases from S (pure noise) to 0 (pure data), ω min and ω max are the lower and upper bounds of the constraint weight respectively, and γ is the adjustment index.

[0155] This mechanism ensures that the constraint weight is small in the early stage of inverse diffusion (s close to S), allowing the model to fully explore the data distribution; while in the later stage (s close to 0), the constraint weight increases significantly, forcing the generated results to strictly converge to the physical constraint surface. At each sampling step, the gradient of the physical constraint loss function is used to correct the sampling result:

[0156] x s-1 * =x s-1 -ω(s)*▽x L physics (x s-1 );

[0157] where x s-1 is the original sample (unmodified) generated by the diffusion model (or other generative model) at the s-1 step of the reverse denoising process, x s-1 * is the sample after physical constraint gradient modification, closer to the feasible solution that satisfies the conservation law, ∇ x represents the gradient operator with respect to the vector variable x, L physics is the weighted sum of the above three constraint violation degrees.

[0158] In addition, in order to further enhance the hydrological feature preservation ability of the denoising result, the embodiment also introduces a conditional diffusion guidance mechanism based on hydrological features. Key feature vectors are extracted from the original jagged sequence:

[0159] h=[Q peak ,T peak ,W total ,T rise ,T recession ,Q base ] T , including peak flow, peak time, total flood volume, rising duration, falling duration and base flow. The feature vector h is integrated into the denoising network U-Net through the cross-attention mechanism, so that the denoising process is not only smoothing the fluctuations, but also reconstructing and repairing the key hydrological elements. Through the physical consistency verification link, the deviation of the generated sequence and the original sequence in the total amount, peak and shape is calculated, and if the deviation exceeds the threshold, the gradient descent method is used to iteratively modify the result until all physical indicators are met.

[0160] Embodiment seven, the specific structure and hardware implementation of the driving mechanism analysis system of reservoir group dispatching on the non-consistent evolution of flood are described in detail. The system is the physical carrier for implementing the above-mentioned method embodiments, and can be deployed in the dispatching center of the basin management institution or the cloud server cluster.

[0161] A driving mechanism analysis system of reservoir group dispatching on the non-consistent evolution of flood, comprising: a memory for storing a computer program; a processor for executing the computer program to realize the steps of the driving mechanism analysis method of reservoir group dispatching on the non-consistent evolution of flood according to any one of the embodiments of the present application.

[0162] The system mainly comprises a memory and a processor. The memory is used to store computer programs and related data, including but not limited to a historical hydrological database, a static topology attribute library, pre-trained denoising model parameters, graph neural differential equation model parameters, and causal analysis model parameters. The processor is used to execute the computer programs to implement each step in the above method embodiments. Specifically, in terms of functional module division, the system comprises a data acquisition and processing module, a topology evolution modeling module, a causal driving analysis module, and a scheduling strategy generation module.

[0163] The data acquisition and processing module is used to execute the process of obtaining a smooth reservoir inflow flood sequence and a dynamic topology feature vector, and a denoising process. The module receives water level and rainfall data from a water regime telemetry station in real time through a communication interface, calculates jagged reservoir inflow using a built-in water balance back-calculation algorithm, and calls a conditional reversible normalization flow model or a hierarchical constraint diffusion model loaded in the processor to output a smooth reservoir inflow flood sequence. At the same time, the module is also responsible for collecting meteorological forecast data as an exogenous hydrological driving input.

[0164] The topology evolution modeling module is used to execute graph neural differential equation solving. The module maintains a dynamic graph structure object in memory, uses a numerical integral solver (such as an RK4 solving unit) to calculate the topology evolution equation in real time, updates the node state vector and edge weight matrix, and generates a dynamic topology feature vector. The module supports parallel computing and can quickly process high-dimensional state evolution of large reservoir groups.

[0165] The causal driving analysis module is an important analysis unit of the system, used to perform structural causal modeling and counterfactual inference. The module internally integrates a causal generalized additive model cGAM fitting engine and a causal inference engine. The fitting engine uses historical data to regularly update the parameters of cGAM, such as the spline function coefficients of location, scale, and shape parameters. The causal inference engine responds to user analysis requests, performs do operator intervention, calculates counterfactual expectations under different scheduling scenarios, and outputs local causal driving index LCE and global causal cumulative effect index GCE in real time, and displays them on the user terminal in the form of visual charts, intuitively presenting a contribution heat map of each reservoir to flood variation.

[0166] The scheduling strategy generation module is used to execute the reinforcement learning scheduling in embodiment five. The module is deployed with a structure-improved Transformer policy network, which receives LCE index from the causal analysis module as part of the reward signal. In offline training mode, the module uses historical flood scenarios for a large number of simulation iterations to optimize the policy network parameters; in online operation mode, the module quickly generates optimized scheduling instructions containing non-consistency suppression targets based on the current water regime state, for decision makers to refer to or directly send to the reservoir on-site control unit.

[0167] At the hardware level, the processor can adopt one or more combinations of central processing unit (CPU), graphics processing unit (GPU), tensor processing unit (TPU), or field programmable gate array (FPGA). Considering that the manifold denoising model and the Transformer network involve a large number of matrix operations, a high-performance GPU cluster is preferably configured for acceleration. The memory can include high-speed random access memory (RAM) and non-volatile memory (such as a hard disk, flash memory) to meet the read-write requirements of massive hydrological data. The system also includes a display and an input device for human-computer interaction. Through the above-mentioned cooperation of software and hardware, the system can efficiently and stably realize deep analysis and accurate regulation of the non-consistent driving mechanism of floods.

[0168] Embodiment eight, a preferred architecture for building a deep reinforcement learning model is described in detail - a parallel subsystem coordination architecture, which further solves the problem of dimension disaster caused by high state space dimension in large reservoir group dispatching. By decomposing the giant system into several subsystems and introducing a coordinator mechanism, the unification of local decision and global optimization is realized, as shown in Figure 4 .

[0169] Step 801, based on the reservoir group dispatching model of each stage, the system is divided into several subsystems, the reservoir constraints of the subsystems are mapped into state vectors, and the target discharge is mapped into action vectors, and a high-dimensional action-state space is constructed.

[0170] In this step, first, the decoupling and space reconstruction of the system are performed. For a complex network composed of dozens of reservoirs in a basin, direct use of centralized control will cause the state space to explode exponentially. Therefore, according to the closeness of the hydraulic connection, such as the dynamic edge weight calculated in the foregoing embodiments, the entire reservoir group is divided into K relatively independent subsystems, such as the upstream cascade group and the midstream cascade group, using the spectral clustering algorithm. For the kth subsystem, the internal reservoir operation constraints (such as reservoir capacity V t , inflow Q in,t , maximum water level limit Z max ) are mapped into local state vectors s k,t . At the same time, the target discharge of all reservoirs in the subsystem is mapped into a local action vector a k,t . In this way, the original global space with a dimension of M*N is divided into K subspaces with smaller dimensions, reducing the learning difficulty of the individual model.

[0171] Step 802, configure an independent sub-policy network for the subsystem, introduce a coordinator network to converge the hidden layer states of the sub-policy network to generate a global coordination signal, and establish a policy improvement type deep deterministic policy gradient of the parallel sub-policy-coordinator architecture.

[0172] This step constructs a hierarchical and collaborative control architecture. Specifically, an independent sub-policy network π is configured for each subsystem k. sub,k (s k ;θ k ), where π sub,k (.) represents the local policy network of the k-th subsystem, s k Let θ be the local state of subsystem k. k For sub-policy network π sub,k The parameters are specified. This network typically employs a multilayer perceptron (MLP) structure, with the local state s as the input. k It outputs local action suggestions. To prevent suboptimal solutions from arising due to a lack of communication between subsystems—for example, an upstream unit cutting off water flow to generate its own electricity, causing downstream units to dry up—this embodiment introduces a global coordinator network π. coord The coordinator receives the hidden feature representations h of all sub-policy networks. k A global coordination signal c is generated through a convergence function (such as mean aggregation or attention aggregation). t =σ'(Agg(h1,...,h K ), where σ'(.) is the nonlinear activation function, Agg(.) is the convergence function, and K is the total number of subsystems. Furthermore, the coordination signal c t Feedback is given to each sub-policy network to correct its final output: a k,t =π sub,k (s k ,c t During training, an improved DDPG algorithm (Deep Deterministic Policy Gradient Algorithm) is used, and a Critic network (evaluation network) evaluates the global state value including coordination signals, guiding each subsystem to pursue local optima while taking into account the global flood control objective.

[0173] Step 803: Dynamically reduce the dimensionality of the high-dimensional state-action space based on the policy-improved deep deterministic policy gradient, extract the low-dimensional feature space, train the policy network in this space, inject adversarial noise to enhance the robustness of the policy network, and obtain the subsystem coordination model.

[0174] To further improve the training efficiency, the embodiment also introduces a self-encoder for dynamic dimension reduction. An encoder Enc(s, a) (an encoder function or neural network module that encodes the input state-action pair (s, a) to achieve dynamic dimension reduction) is constructed to compress the high-dimensional joint state-action pair into a low-dimensional feature vector z, and search for a policy on the low-dimensional manifold. To enhance the robustness of the model to hydrological uncertainties (such as prediction errors), an adversarial noise z' = z + ε' is injected into the low-dimensional feature z during the training phase, where the small random disturbance ε' follows a truncated normal distribution. By forcing the policy network to still output stable scheduling instructions under noise interference, the anti-interference ability of the model in actual application is improved.

[0175] Embodiment Nine, details the specific internal implementation details of the structure improved Transformer as a policy network. This network structure is designed for the long-term and short-term dependencies and non-stationary characteristics of hydrological sequences, integrates gate position perception and quantum-inspired encoding technology, and is a computing engine for achieving high-precision scheduling.

[0176] Step 901, a dynamic position perception mechanism with gating adjustment is introduced in the encoder position encoding to capture non-stationary time series dependencies, while the subsystem feature coupling weight is added in the multi-head self-attention, and a quantum-inspired decoding technology is combined to establish a structure improved Transformer.

[0177] In this step, the position encoding of the standard Transformer is reformed for hydrology. The traditional sinusoidal position encoding is static and cannot perceive the dynamic phase changes of the flood waveform, such as early flood or late flood. This embodiment proposes a dynamic gating position perception mechanism, whose calculation formula is:

[0178] PE(pos, t) = sin(pos / 10000 2i' / d )+g t *cos(t / 10000 2i' / d );

[0179] where pos is the sequence position index, t is the physical time step, i' is the dimension index, and d is the model dimension. g t is an adaptive coefficient generated by a gating network (a simple fully connected layer) according to the current water regime characteristics. When the water regime is stable, g t tends to 0, degenerating into standard encoding; when a mutation occurs, g t increases, introducing explicit constraints on physical time to capture non-stationary time series dependencies.

[0180] In addition, in the multi-head self-attention mechanism, the physical coupling information between subsystems is explicitly injected. The standard attention score calculation formula is Softmax(QK Tsqrt(d k )),the embodiment corrects it to:

[0181] Attention(Q,K,V) = Softmax(QK T / sqrt(d k )+ gamma * C sub )*V;

[0182] wherein (Q,K,V) is a query (Query), key (Key), and value (Value) matrix, Attention(Q,K,V) is the corrected attention output, which fuses data-driven similarity and physical structure constraints, d k is the dimension of the key / query vector of each attention head, C_sub is the subsystem coupling matrix constructed based on the dynamic topological feature T t , and gamma is a learnable coupling weight coefficient. This improvement ensures that the attention mechanism can focus on the upstream or downstream nodes with strong connections in the water system, rather than blindly calculating global dependencies.

[0183] Step 902, map the multivariate input feature sequence into a high-dimensional embedding vector, input the high-dimensional embedding vector into the multi-head self-attention layer and feedforward neural network of the structure-improved Transformer encoder, extract the coupling time series features between subsystems, use the decoder to generate the outflow sequence of each subsystem, and output the flood process after the reservoir group regulation in each stage.

[0184] At the decoder end, the embodiment introduces quantum-inspired encoding technology to enhance the diversity and search ability of the generated sequence. Specifically, the Hadamard quantum gate operator H is used to transform the hidden layer state |psi in > of the decoder:

[0185] |psi out > = H |psi in >;

[0186] wherein |psi in > is the initial quantum state (input state) of the system, and |psi out > is the final state (output state) of the system after time evolution or operation.

[0187] Mathematically, it is equivalent to orthogonal rotation and superposition processing of the feature vector, so that the model can explore multiple potential scheduling trajectories in the state space simultaneously (similar to quantum superposition state), effectively jumping out of the local optimal solution. |psi out > is mapped to the outflow sequence Q out of each subsystem through a linear projection layer, completing the generation of scheduling decisions.

[0188] According to one aspect of this application, the construction of structural causal graphs and structural equations in the process of quantifying the driving effect of phased scheduling of reservoir groups on non-uniform flood evolution is as follows:

[0189] Define the key variables at time t:

[0190] M i,t Indicates the first The reservoir is at all times The scheduling capacity index (MRI value) is given by i, which is the reservoir number index, with values ​​of 1, 2, ..., N, where N is the number of reservoirs.

[0191] M t Let represent the index vector of the scheduling capacity of all reservoirs at time t, i.e.:

[0192] M t =[M 1,t M 2,t ,…,M N,t ] T The superscript T indicates the transpose operation.

[0193] R t The exogenous hydrological driving feature vector at time t represents factors such as rainfall intensity field within the watershed, upstream water flow, and initial soil moisture content.

[0194] T t This represents the topological feature vector output by the dynamic topology model, used to characterize the upstream and downstream relationships, parallel relationships, and topological connectivity of a multi-reservoir system at time t.

[0195] Y t The vector representing the flood inconsistency index at time t is, for example, the principal component score, which is composed of original indices such as DTW distance and peak deviation.

[0196] ε t M ε represents the exogenous noise vector in the scheduling capacity exponential equation. t Y This represents the exogenous noise vector in the non-consistent index equation.

[0197] The above variables constitute the node set: V={M t ,R t ,T t ,Y t}

[0198] A structural cause-effect graph can be represented as a directed graph as follows:

[0199] G=(V,E), where G represents the structural causal graph object, and E represents the set of directed edges in the graph, describing the causal directions. For example:

[0200] Mt =f M (R t ,T t ,ε t M );

[0201] Y t =f Y (M t ,R t ,T t ,ε t Y );

[0202] In the formula, the function f M (.) represents the generating function of the dispatch capacity index, with the input being the exogenous hydrological driving vector R. t Topological feature vector T t and noise vector ε t M The output is the scheduling capability vector M. t .

[0203] function f Y (.) represents the generating function for the flood inconsistency index, with the input being the dispatch capacity vector M. t Exogenous hydrological driving vector R t Topological feature vector T t and noise vector ε t Y The output is a non-consistent index vector Y. t .

[0204] According to one aspect of this application, the construction and counterfactual computation of a causal generalized additive model (cGAM) are carried out in the process of quantifying the driving effect of phased scheduling of reservoir groups on the non-uniform evolution of floods, specifically as follows:

[0205] Based on the structural equation, the non-consistent index component Y j,t (For the j-th index) Select a suitable probability distribution family. Assume Y j,t Obtained by parameter μ j,t σ j,t ν j,t Distribution:

[0206] Y j,t ~D(μ j,t ,σ j,t ,ν j,t );

[0207] Wherein, the symbol Y j,t Represents vector Y tthe jth component of the jth index, j is the index number index; the symbol ~ means the distribution subject to …; the symbol D(.) means a given distribution family, such as normal distribution, Gamma distribution or other distributions suitable for non-uniform indexes; μ j,t represents the mean or location parameter; σ j,t represents the scale parameter; v j,t represents the shape parameter.

[0208] In the causal generalized additive model, the relationship between the parameter and the causal parent node is represented by an additive smooth function, for example:

[0209] ;

[0210] ;

[0211] In the formula, α j,0 represents the constant term of the jth index on the location parameter, β j,0 represents the constant term of the jth index on the scale parameter; the function s μ,j i (.) and s σ,j i (.) represent the smooth contribution function of the MRI component M i,t ; the commonly used spline function is used to realize it; the function s μ,j R (.) and s σ,j R (.) represent the smooth contribution of the exogenous hydrological driving vector R t ; the function s μ,j T (.) and s σ,j T (.) represent the smooth contribution of the topological feature vector T t .

[0212] Causal intervention is represented by the do operator. When the scheduling ability vector M t is intervened, it is forcibly set to a given vector , which is denoted as . Under this intervention, the conditional expectation of the jth non-uniform index can be written as:

[0213] ;

[0214] In the formula, the symbol E[.] represents the mathematical expectation operator, the symbol | represents under the condition of …, and the symbol do(.) represents the causal intervention operator, , which represents the scheduling ability vector after intervention.

[0215] and the traditional conditional expectation ​Unlike the above, the expectation is to directly replace the independent variable M in the structural equation. t This allows for the estimation of inconsistencies in counterfactual scenarios while keeping exogenous noise and other structures constant.

[0216] According to one aspect of this application, in the topologically coupled reversible physical manifold diffusion-flow denoising method, the conditionally reversible normalized flow constructs the physical-topological manifold, specifically as follows:

[0217] Suppose there are K historical flood events, numbered k = 1, 2, ..., K. For the k-th flood event, let the discrete time series of the sawtooth inflow flood be:

[0218] q k =[q1 k ,q2 k ,…,q L k ] T ;

[0219] In the formula, q k Indicates the first The inflow vector of the flood; q l k This represents the inbound flow value of the event at time step l; l is the time step index, with values ​​of 1, 2, ..., L; L represents the total number of discrete time steps.

[0220] Let the dynamic topological eigenvector corresponding to the k-th flood be:

[0221] τ k =[τ1 k ,τ2 k ,…,τ P k ] T ;

[0222] In the formula, τ k τ represents the topological eigenvector of the k-th flood; P k This represents the p-th topological feature component, such as certain upstream and downstream connectivity indicators or series-parallel structure indicators; p is the feature index, with values ​​of 1, 2, ..., P; P represents the dimension of the topological feature.

[0223] Define a parameter as Conditional invertible transformation T θ :q k ,τ k ->z k ;

[0224] That is: z k =T θ (q k ,τk );

[0225] where T θ (.) denotes the invertible normalization flow function; z k denotes the representation vector of the kth flood in the latent space; arrow -> denotes the mapping relationship; parameter θ denotes the set of all trainable parameters in the transformation.

[0226] Assume that the prior distribution of z k in the latent space is a multivariate standard normal distribution:

[0227] p Z (z k )=N(z k ;0,I);

[0228] where p Z (.) denotes the probability density function of the latent variable; N(.;0,I) denotes a multivariate normal distribution with mean vector zero and covariance matrix I; 0 denotes the all-zero vector; I denotes the identity matrix.

[0229] According to the transformation formula of the invertible flow, the probability density of the kth flood in the observation space is:

[0230] ;

[0231] where p (q k |τ k ) denotes the conditional probability density of the outflow vector q k under the condition of topological characteristics τ flow ; symbol det(.) denotes the matrix determinant; the numerator in the fractional form denotes the Jacobian matrix of the transformation function with respect to the input vector, and the denominator denotes the derivative with respect to the input vector q mass .

[0232] In the training stage, a loss function based on the log-likelihood is defined, and a water conservation constraint is added:

[0233] ;

[0234] where L out,l denotes the total loss function of the invertible flow; log(.) denotes the natural logarithm; λ k is the weight coefficient of the water conservation constraint; q l=1 L denotes the outflow of the kth flood at the lth time step; ∑ l k denotes the time summation of the inflow sequence of the event; ∑ l=1L q out,l k denotes the time summation of the outflow sequence.

[0235] After training, the reversible flow transformation T θ The inflow sequence satisfying the water conservation constraint is compressed into a nearly Gaussian representation in the latent space, and the dynamic topological features are explicitly encoded.

[0236] According to an aspect of the present application, in the topologically coupled reversible physical manifold diffusion-flow denoising method, the latent space diffusion denoising and the local physical constraint are specific to:

[0237] In the latent space, the latent representation z k A discrete-time diffusion process is constructed. Let the diffusion time step be t = 1, 2, …, T, and define the forward diffusion process as:

[0238] ;

[0239] where z t k denotes the latent vector of the kth flood event at diffusion time step t; α t is a scaling coefficient located in the interval (0, 1), which is usually determined by a set of β t coefficients through the relationship α t = 1-β t ; ε t k denotes a Gaussian noise vector at time step t, which follows a standard normal distribution; the subscript t-1 denotes the variable of the previous diffusion time step.

[0240] In the reverse denoising process, a noise estimation network ε ϕ (.) with parameter ϕ is constructed, and the input is the noisy latent vector z t k , topological features τ k , and diffusion time step t. The training objective is to minimize the following loss:

[0241] ;

[0242] where L diff denotes the total loss of the diffusion model; ||.||2 denotes the Euclidean norm; λ phys is the weight coefficient of the local physical constraint loss; L phys denotes the physical constraint loss term, such as monotonicity constraint, slope constraint, etc.

[0243] To introduce local physical constraints into the latent space generation result, a physical loss needs to be defined after generating the smooth inflow hydrograph sequence. Let the smooth inflow sequence obtained through the reversible mapping be:

[0244]

[0245] Then the slope constraint loss can be defined as:

[0246]

[0247] wherein, represents the smooth inflow flow of the kth flood at the lth time step; represents the smooth inflow flow of the previous time step; max represents the maximum slope threshold value allowed; max(0,.) represents taking the larger value in the parentheses and zero, used to define the penalty only when the threshold is exceeded; phys L slope and other physical constraint terms, such as peak monotonicity constraint, wave front timing constraint, etc.

[0248] According to an aspect of the present application, a method for analyzing the driving mechanism of reservoir group scheduling on non-uniform evolution of flood can further include the following steps:

[0249] Collecting multi-source heterogeneous hydrological data of reservoir geographic coordinates, operation time, and river network connection information, obtaining an optimized time-varying graph structure through dynamic topological mapping and graph convolution network (GCN) spatiotemporal feature extraction, driving node state and edge weight update in an ordinary differential equation (ODE) model by a topological evolution equation, optimizing the topological evolution trajectory through a graph attention mechanism and an adaptive solving strategy, and establishing a dynamic reservoir group topological model;

[0250] Based on the dynamic reservoir group topological model, gradually incorporating the scheduling system according to the reservoir operation sequence to obtain the topological model at each stage, collecting historical flood and operation records, combining the water balance equation to back-propagate the jagged inflow hydrograph, using the forward diffusion of the denoising diffusion probability model to add noise, and in the inverse diffusion sampling, combining the water conservation knowledge guidance mechanism to generate a smooth inflow hydrograph, and based on the Muskingum method to calculate the inflow to the downstream superposition area to obtain the natural flood without reservoirs;

[0251] ​​A dispatching model framework is established with section safety as the target, each stage topological model is incorporated into the establishment of each stage dispatching model, the dispatching model is mapped into a high-dimensional state-action space, a strategy improvement type DDPG (deep deterministic policy gradient algorithm) is used for dynamic dimension reduction and training to generate a dispatching decision sequence, a structure improvement type Transformer embedded quantum heuristic coding is used to extract time sequence features, and a post-regulation flood process of each stage is obtained, and a non-consistency index sequence of the flood is calculated by comparing the natural flood;

[0252] Based on the effective flood control storage of the reservoir, the dispatching strategy and the uncertainty factor of subsequent rainfall, an improved reservoir index MRI sequence is established, a generalized additive model is used to fit the MRI and the non-consistency index sequence to obtain a nonlinear model structure, a Bootstrap sampling (bootstrap sampling method) is used to generate and screen the nonlinear mapping function with the best accuracy, and based on the first-order difference and trapezoidal integral of the mapping function, an MRI influence intensity index and a cumulative influence index are defined to quantify the driving effect of the phased dispatching of the reservoir group on the non-consistency evolution of the flood.

[0253] The application replaces the traditional statistical correlation analysis by constructing a structural causal model and introducing a counterfactual intervention (do operator). This method mathematically cuts off the confounding path, successfully separates the influence of reservoir dispatching on flood non-consistency from natural rainfall and topological changes, accurately quantifies the driving mechanism, and solves the causal decoupling problem.

[0254] The application uses a topologically coupled conditional reversible normalization flow and diffusion model, and introduces water conservation and manifold regularization constraints in the latent space. It solves the water imbalance and waveform distortion problems caused by traditional filtering methods, provides high-fidelity data that meets the laws of hydraulics for analysis, and solves the problem of data physical consistency.

[0255] The application establishes a graph neural ordinary differential equation model that introduces a propagation time lag term and a flow sensitivity weight evolution equation. It overcomes the limitations of static graph models that cannot describe the delay of flood wave evolution and the dynamic changes of structure, accurately describes the time-varying hydraulic connection of the reservoir group, feeds the causal driving index back to the reinforcement learning reward function to realize the active optimization of the dispatching strategy, and solves the problem of dynamic topology capture.

[0256] The above detailed description of the preferred embodiments of the application, but the application is not limited to the specific details in the above embodiments, within the technical concept of the application, various equivalent transformations of the technical solutions of the application can be made, and these equivalent transformations all belong to the protection scope of the application.

Claims

1. A method for analyzing the driving mechanism of reservoir group operation on the non-uniform evolution of floods, characterized in that, The method comprises the following steps: obtaining a smooth reservoir inflow flood sequence denoised by topological coupling physical manifold constraint, and a dynamic topological feature vector representing the hydraulic propagation time lag relationship between upstream and downstream reservoir groups; based on the smooth reservoir inflow flood sequence, the dynamic topological feature vector and the real-time collected exogenous hydrological driving data, a structural causal model containing a scheduling capability index sequence and a flood inconsistency index is constructed; nonlinear fitting is performed on the variable relationship in the structural causal model to determine the causal dependence relationship of the scheduling capability index sequence on the flood inconsistency index; based on the causal dependence relationship, counterfactual inference is performed, the scheduling capability index sequence is forcibly set to a reference value through an intervention operator, and the conditional expectation of the flood inconsistency index in the counterfactual scenario is calculated; based on the difference between the conditional expectation and the observed expectation, a local causal driving index representing the instantaneous influence of reservoir scheduling on flood variation is calculated, and the local causal driving index is integrated in the time dimension to obtain a global causal cumulative effect index representing the long-term cumulative contribution; based on the global causal cumulative effect index, an optimized scheduling strategy for inhibiting the evolution of inconsistency is generated; wherein, the dynamic topological feature vector representing the hydraulic propagation time lag relationship between upstream and downstream reservoir groups is obtained by: based on the smooth reservoir inflow flood sequence and the basic attribute data of the reservoir group, a graph neural differential equation model is constructed, a topological evolution equation describing the evolution of the state of the reservoir node with time is established, and the topological evolution equation is expressed as: ； where dx i is the state vector of the ith reservoir node i is the derivative with respect to time t, f local is the local dynamics function describing the behavior of a single reservoir, U i is the set of upstream reservoirs of the ith reservoir node, w ji is the edge weight from upstream node j to the current node i at time t, g up is the coupling function describing the upstream propagation influence, τ ji is the water flow propagation time lag from node j to node i, P i is the set of reservoirs in parallel with the ith reservoir node, w ki is the edge weight from parallel node k to the current node i at time t, g para is the function describing the parallel coupling influence; the topological evolution equation is solved by using a numerical integral solver to obtain an evolution trajectory containing a time-varying node state vector and a time-varying edge weight matrix, and the evolution trajectory is taken as the dynamic topological feature vector.

2. The method of claim 1, wherein, The smooth reservoir inflow flood sequence denoised by the topological coupling physical manifold constraint is obtained by: collecting the historical operation records of the reservoir group, and inversely deducing a jagged reservoir inflow flood sequence containing non-physical fluctuations based on the water balance principle; a conditional reversible normalizing flow model is constructed based on the dynamic topological feature vector, and a bidirectional reversible mapping between the jagged reservoir inflow flood sequence and the latent space of the Gaussian latent variable is established by using the conditional reversible normalizing flow model; a denoising diffusion probability model is constructed in the latent space, forward noise is added to the Gaussian latent variable to learn the noise distribution, and noise is removed by reverse diffusion sampling to obtain a denoised latent manifold representation; the denoised latent manifold representation is mapped back to the physical data space by using the inverse transformation of the conditional reversible normalizing flow model to obtain a smooth reservoir inflow flood sequence conforming to the hydrodynamic conservation constraint.

3. The method of claim 2, wherein, The training process of the conditional reversible normalizing flow model and the denoising diffusion probability model contains a physical manifold regularization constraint, and the physical manifold regularization constraint comprises: a water conservation constraint term, the calculation formula is: ； wherein L mass is the loss value of the water conservation constraint term, λ mass is the weight coefficient of the water conservation constraint, K is the total number of flood events, L is the total number of time steps of a single flood event, q l k is the inflow of the kth flood event at the lth time step generated by the model, q out,l k is the outflow of the kth flood event at the lth time step. a cross-event manifold regularization term, the calculation formula is: ； where L manifold is the loss value of the cross-event manifold regularization term, z k and z k' are the latent manifold representation vectors of the kth event and k'th event flood in latent space, respectively, ||·||2 2 denotes the squared Euclidean distance, w k,k' is the similarity weight between events, which is calculated based on the distance of the exogenous hydrological driving data of the two floods.

4. The method of claim 1, wherein, The state vector is defined as x i = [V i , Q in,i , Q out,i , Z i ] T , where V i , Q in,i , Q out,i , Z i are the storage, inflow, outflow, and water level of the ith reservoir, respectively. Local dynamics function f local Specifically represented as: ； wherein, κ1 is the reservoir inflow recovery coefficient, Q in,i eq To balance the reservoir inflow, is a theoretical discharge capacity function based on reservoir characteristics, ψ(V i is a water level calculation function based on the reservoir capacity curve. Coupling function g describing the upstream propagation influence up Specifically represented as: ； Where, η ji Let Q be the propagation attenuation coefficient from upstream node j to current node i. out,j (t-τ ji This indicates that upstream node j is lagging behind by a time delay τ. ji The outbound flow rate at any given time.

5. The method of claim 1, wherein, The topological evolution equation also includes an edge weight evolution equation, which is defined as: ； where dw ij denotes the edge weight w ij denotes the derivative of time t, λ w is the weight recovery coefficient, w ij base is the basic edge weight determined based on river network distance, μ w is the sensitivity response coefficient, σ(·) is the Sigmoid activation function, is the partial derivative of the current node i's inflow to the upstream node j's outflow, which is used to represent the real-time sensitivity of flow change.

6. The method of claim 1, wherein, The nonlinear dependence relationship between variables in the structural causal model is described by a causal generalized additive model, wherein the flood inconsistency index is subject to a preset probability distribution, and the location parameter, scale parameter and shape parameter of the probability distribution are all constructed as additive smooth functions of the causal parent nodes; The additive smoothing function is linearly superimposed by a smoothing component for the scheduling capability index sequence, a smoothing component for the exogenous hydrological driving data, and a smoothing component for the dynamic topology feature vector; The functional expression of the position parameter is: ； where μ j,t is the position parameter of the jth flood non-uniformity index at time t, a j,0 is the constant intercept term, N is the total number of reservoirs, M i,t is the scheduling ability index of the ith reservoir at time t, s μ,j i (.) is the spline smoothing function for scheduling ability index, R t is the exogenous hydrological driving data, s μ,j R (.) is the spline smoothing function for exogenous driving, T t is the dynamic topological feature vector, s μ,j T (.) is the spline smoothing function for topological features.

7. The method of claim 1, wherein, The calculation formula of the local causal driving index is: ； where LCE i,j,t is the local causal driving index of the i-th reservoir to the non-uniformity index of the j-th flood at time t, denotes the partial derivative operator of the dispatching capacity index, E[·|do(·)] denotes the counterfactual conditional expectation after performing the intervention operator, M 0 i,t is the reference value under the current actual dispatching scheme; The calculation formula of the global causal cumulative effect index is: ； where GCE i,j is the global causal cumulative effect index of the ith reservoir on the jth flood inconsistency index, t start and t end are the starting and ending time of the analysis period, respectively, t is the discrete time step.

8. The method of claim 1, wherein, The optimization scheduling strategy for generating inhibition of inconsistent evolution includes: constructing a deep reinforcement learning model, and defining an immediate reward function containing multiple constraint targets; The calculation formula of the immediate reward function is: ； wherein R t is the instant reward value at time t, λ1, λ2, λ3 are non-negative weight coefficients of flood control safety, non-uniformity level and causal driving strength respectively, Q max,t is the maximum discharge of the cross section, Q safe is the safety discharge threshold, NCI t is the comprehensive non-uniformity level index, N is the number of reservoirs, J is the number of non-uniformity indexes, |LCE i,j,t is the absolute value of the local causal driving index.

9. A system for analyzing driving mechanism of non-uniform evolution of flood for reservoir group operation, characterized in that, The method comprises the following steps: A memory is configured to store a computer program. A processor is configured to execute the computer program to implement the steps of the reservoir group scheduling driving mechanism analysis method for flood inconsistent evolution according to any one of claims 1 to 8.

Images