Reservoir group-lake system flood control risk adaptive analysis method and system

Abstract

The invention discloses a reservoir group-lake system flood control risk adaptive analysis method and system. The method comprises the following steps: acquiring a basin hydrological sequence and a pre-configured physical causal constraint matrix, and constructing a VineCopula joint distribution and initial Bayesian network conforming to physical topology according to the basin hydrological sequence and the pre-configured physical causal constraint matrix; risk sensitive weights of the sample points are calculated based on the Bayesian network, and weighted sampling density inclining towards the high-risk area is constructed; generating a self-adaptive scene set based on the density, solving an optimal scheduling scheme, feeding back a scheduling result to the Bayesian network to update parameters, carrying out closed-loop iteration of perception-sampling-scheduling-learning until risk assessment is converged, and determining a comprehensive flood control risk assessment result of the system based on the converged risk factor Bayesian network. According to the method, the problems of poor physical interpretability and scarcity of extreme risk samples of a traditional model are solved, and accurate capture and quantification of extreme risks of a complex flood control system are realized.

Description

Technical Field

[0001] This invention relates to the field of flood control safety and risk management in water conservancy projects, and in particular to an adaptive analysis method and system for flood control risks in reservoir-lake systems. Background Technology

[0002] With the increasing frequency of extreme hydrological events triggered by climate change, complex flood control systems composed of reservoirs, lakes, and rivers face severe challenges. Accurately identifying the correlation structure of multiple risk factors and quantifying the comprehensive flood control risk of the system under complex scheduling is of great research significance for formulating scientific watershed flood control plans and emergency response plans.

[0003] Currently, risk analysis for reservoir-lake systems mainly employs a combination of joint distribution modeling based on Copula functions and Monte Carlo stochastic simulation. Existing research typically utilizes VineCopula (a high-dimensional joint distribution modeling method) to construct the joint distribution of multidimensional hydrological variables and describes risk transmission paths through Bayesian networks. In scenario simulation, Latin hypercube sampling (LHS) or importance sampling is often used to generate flood combination scenarios, which are then input into reservoir operation models to calculate risk indicators after system regulation.

[0004] However, existing technologies still face significant technical bottlenecks when dealing with high-dimensional complex systems. These bottlenecks are mainly manifested in the lack of physical causality due to purely data-driven approaches and the underestimation of extreme risks caused by static sampling. Specifically, traditional VineCopula structure learning relies solely on selecting the strongest statistical correlations, ignoring the objective spatiotemporal causality of water flow evolution. This often results in statistically strongly correlated but physically meaningless erroneous connections in the model, and the generated Bayesian network exhibits uncertain orientation and lacks interpretability. Furthermore, existing stochastic simulations employ one-time static sampling, failing to perceive the risk distribution after scheduling. This leads to extremely sparse samples in the most vulnerable high-risk areas of the system (such as severe flood combinations encountering unfavorable scheduling boundaries), resulting in insufficient revelation of system failure mechanisms under extreme conditions and distorted risk assessment results. Summary of the Invention

[0005] The purpose of this invention is to provide an adaptive analysis method and system for flood control risks in reservoir-lake systems, in order to solve one of the aforementioned problems in the existing technology.

[0006] According to one aspect of this application, a method for adaptive analysis of flood control risk in a reservoir-lake system includes:

[0007] Obtain the watershed hydrological sequence and pre-configured physical causal constraint matrix, and construct the VineCopula joint distribution model and initial risk factor Bayesian network that satisfy the physical causal topological constraints.

[0008] Based on the initial risk factor Bayesian network, the posterior probability of each sample point causing the system risk index to exceed the preset risk threshold is calculated, and the risk-sensitive weight function is obtained.

[0009] The risk-sensitive weighting function is coupled to the VineCopula joint distribution model to construct the risk-weighted sampling density, and a risk-adaptive scenario set is generated based on this density.

[0010] The optimal scheduling scheme for the reservoir group is solved for the risk adaptive scenario set and fed back to the risk factor Bayesian network to update the network parameters. The steps of calculating the risk-sensitive weight function and updating the network parameters are repeated until the risk assessment results converge.

[0011] Based on the converged risk factor Bayesian network, the comprehensive flood control risk assessment results of the system are determined.

[0012] According to another aspect of this application, an adaptive analysis system for flood control risk in a reservoir-lake system includes:

[0013] processor;

[0014] Memory, used to store computer programs;

[0015] The processor executes the steps described in any one of the methods of the present invention when executing a computer program.

[0016] Beneficial effects: This invention solves the problems of poor physical interpretability and scarcity of extreme risk samples in traditional models, and achieves accurate capture and quantification of extreme risks in complex flood control systems. Attached Figure Description

[0017] Figure 1 A schematic diagram of the overall process of the adaptive analysis method for flood control risk of reservoir group-lake system provided in the embodiments of this application.

[0018] Figure 2 This is a schematic diagram of the first-layer tree structure process for obtaining the VineCopula joint distribution model, provided in an embodiment of this application.

[0019] Figure 3 This is a schematic diagram of the Bayesian network process for obtaining the initial risk factor provided in the embodiments of this application.

[0020] Figure 4 This is a schematic diagram illustrating the process of feeding the optimal scheduling scheme back to the risk factor Bayesian network to update network parameters, as provided in an embodiment of this application. 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: A comprehensive implementation scheme for an adaptive analysis method of flood control risk in a reservoir-lake system under physical constraints is provided, such as... Figure 1 As shown, by introducing a physical causal constraint matrix and a risk adaptive iteration mechanism, the technical problems of lack of physical interpretability of statistical structures, scarcity of extreme risk samples, and low simulation efficiency in traditional pure data-driven methods are solved, thereby achieving accurate quantification and dynamic assessment of the risks of complex flood control systems.

[0023] Step 101: Obtain the watershed hydrological sequence and the pre-configured physical causal constraint matrix, and based on the watershed hydrological sequence and the physical causal constraint matrix, construct the VineCopula joint distribution model and the initial risk factor Bayesian network that satisfy the physical causal topological constraints.

[0024] In this embodiment, obtaining the watershed hydrological sequence specifically refers to collecting a long series of historical measured data from reservoirs, lakes, and river sections within the study area, which forms the basis for subsequent probabilistic modeling and risk analysis. Specifically, the data to be collected includes, but is not limited to: river channel calculation parameters (such as wave velocity and diffusion coefficient), reservoir storage processes, measured reservoir outflow processes, lake storage processes, measured lake outflow processes, measured starting water levels of reservoirs and lakes, and external river water level processes. Simultaneously, engineering characteristic data also need to be collected, such as water level-storage capacity curves, water level-discharge capacity curves, and characteristic water levels (such as flood control limits and high flood control levels) for reservoirs and lakes. To obtain a high-quality flow sequence, the inflow to the reservoir is calculated using a back-calculation formula based on the water balance principle. The specific calculation formula is as follows:

[0025] Q i (t)=(V i (t+1)-V i (t)) / Δt+q i (t);

[0026] Among them, Q i (t) represents the inflow rate of reservoir i at time t, V i (t+1) and V i (t) represents the water storage capacity of the reservoir at times t+1 and t, respectively. i(t) represents the outflow from reservoir i at time t, and Δt is the time step. After calculation, the moving average method is usually used to eliminate the sawtooth effect in the inflow process caused by errors in the reservoir capacity curve. For inflow within a lake section, the Muskingan method must first be used to calculate the outflow from the upstream reservoir to the inflow section. The Muskingan method's flow evolution formula is:

[0027] Q out (t)=C0Q in (t)+C1Q in (t-1)+C2Q out (t-1);

[0028] Among them, the coefficients C0, C1, and C2 are determined by the channel parameters K and X, and Q out (t) and Q out (t-1) represent the outflow discharge at the downstream section of the river at time t and the previous time t-1, respectively. in (t) and Q in (t-1) represents the inflow discharge at the upstream section of the river at time t and the previous time t-1, respectively, where K is the storage constant and X is the weighting factor. After obtaining the response discharge of the reservoir outflow at the lake inflow section, the inflow process of the uncontrolled watershed in the lake area is again derived by combining the measured water level changes and lake outflow using the water balance equation. Based on the compiled data on flood peaks, flood volumes, and initial regulation water levels for each flood event, a multidimensional risk factor sample set is constructed. In some embodiments, the initial regulation water level can be the real-time water level during the regulation period.

[0029] While acquiring hydrological sequences, a physical causal constraint matrix needs to be defined. This matrix is ​​a pre-configured logical matrix used to digitally represent the objective physical relationships between various risk factors within the watershed, such as upstream flow inevitably affecting downstream water level, and inflow inevitably affecting outflow. The process of constructing the VineCopula joint distribution model and the initial risk factor Bayesian network using this matrix and the hydrological sequences is achieved by integrating physical prior knowledge into statistical learning algorithms. Specifically, firstly, marginal distributions are fitted to each univariate risk factor. Optional distribution functions include normal distribution, log-normal distribution, Pearson type III distribution, and Gumbel distribution, and the optimal distribution is selected through the KS test and the AIC criterion (Akaike Information Criterion). When constructing the multivariate joint distribution, it does not rely entirely on statistical correlation of the data, but rather mandates that the connection structure between variables must conform to the topological relationship defined by the physical causal constraint matrix, resulting in a physically interpretable first-level tree structure of VineCopula, which is further transformed into the initial risk factor Bayesian network. This process ensures that the subsequent risk inference logic conforms to the basic laws of hydraulics.

[0030] Step 102: Based on the initial risk factor Bayesian network, calculate the posterior probability of each sample point causing the system risk index to exceed the preset risk threshold, and obtain the risk-sensitive weight function.

[0031] This step is the perception stage of the risk adaptive iteration. The risk factor Bayesian network is a directed acyclic graph model that not only encodes the dependencies between variables but also stores the conditional probability table (CPT). The system risk indicator is usually selected as the highest water level Z of key flood control engineering nodes such as lakes during the scheduling period. max Preset risk threshold R th It can be set according to flood control standards, such as guaranteed water level or warning water level. Risk-sensitive weighting function w k (x) represents the occurrence of a high-risk event (i.e., Z) in the system under the current combination of risk factors x. max >R th The probability of (R) is calculated. Specifically, the forward inference function of a Bayesian network is used, with the value of the risk factor vector x as the evidence variable. The value of the target node (system risk indicator) is calculated to be greater than R. th The probability of occurrence. This probability value is the risk weight of that sample point. By traversing the risk factor space or calculating over a large number of sample points, a risk-sensitive weight function covering the entire space can be obtained. This function can identify a region of factor combinations that, although their probability of occurrence may be low (i.e., statistically long-tail), would lead to serious consequences (high risk in a risk sense) once they occur.

[0032] Step 103: Couple the risk-sensitive weight function to the VineCopula joint distribution model to construct the risk-weighted sampling density, and generate a risk-adaptive scenario set based on the risk-weighted sampling density.

[0033] This step is the execution phase of the risk-adaptive iteration. Traditional Copula sampling is usually based on the original joint probability density, which tends to generate a large number of samples in safe regions with high probability density, while the samples are sparse in high-risk regions. To solve this problem, this embodiment constructs a distorted sampling density, namely the risk-weighted sampling density f. k (x). The density function is derived from the original VineCopula density c(u) and the risk-sensitive weight w. kThe product of (x) constitutes the risk factor combination. It is understandable that the higher the probability that a sample point leads to a high risk in the system, the higher its probability density in the new sampling distribution. Based on a non-standard probability density function, Markov chain Monte Carlo (MCMC) sampling or improved Latin hypercube sampling (LHS) techniques are used to generate a new set of risk factor combination samples, i.e., a risk-adaptive scenario set. This scenario set has obvious adverse flood control characteristics, containing a large number of extreme or severe combinations that threaten the flood control system, such as a scenario where a major upstream flood encounters a high initial water level, thus improving the efficiency and relevance of risk analysis.

[0034] Step 104: Solve for the optimal scheduling scheme of the reservoir group for the risk adaptive scenario set, and feed the optimal scheduling scheme of the reservoir group back to the risk factor Bayesian network to update the network parameters. Repeat the steps of calculating the risk-sensitive weight function and updating the network parameters until the risk assessment results converge.

[0035] This step is the closed-loop feedback loop of the risk-adaptive iteration. For each generated risk-adaptive scenario, i.e., a specific combination of flood inflow and initial water level, the reservoir-lake joint optimization scheduling model needs to be run to solve for the optimal discharge strategy. The scheduling model aims to minimize the average remaining excess flood volume of the lakes, considering physical constraints such as water balance, water level limits, and discharge capacity limitations. The solution process can employ intelligent optimization algorithms such as the Sparrow Search Algorithm (SSA) to obtain the optimal scheduling scheme and the corresponding highest lake water level under that scenario. The newly generated scenario-scheduling-result data pairs are then added to the historical training sample set, and the parameters (conditional probability table) of the Bayesian network are retrained using the Expectation-Maximization (EM) algorithm.

[0036] As the proportion of high-risk scenarios in the training samples increases, the updated Bayesian network's ability to predict extreme risks will be enhanced. At this point, the system enters the next iteration: the risk-sensitive weights are recalculated based on the updated network to generate more targeted new samples. The perception-sampling-scheduling-learning loop continues until convergence conditions are met, such as the relative error of the system failure probability calculated in the two rounds being less than a preset threshold (e.g., 0.001), or the preset maximum number of iterations is reached.

[0037] Step 105: Based on the converged risk factor Bayesian network, determine the comprehensive flood control risk assessment results of the system.

[0038] Once the iterations converge, the risk factor Bayesian network has fully learned the full spectrum of risk characteristics from normal to extreme scenarios. At this point, the network can be used for the final risk assessment. The assessment results include, but are not limited to: the system risk probability under any given combination of risk factors (forward prediction), and the most likely combination of each causative factor when the system experiences a certain level of risk (reverse diagnosis). Furthermore, the overall failure probability of the system can be calculated, i.e., the weighted average risk considering all possible combinations of risk factors and their natural occurrence probabilities. The quantitative assessment results provide a scientific and reliable basis for decision-making in flood control planning, dispatching rule formulation, and emergency response plan development.

[0039] Example 2 describes the specific implementation process of constructing a VineCopula joint distribution model and an initial risk factor Bayesian network that satisfy physical causal topological constraints. By defining the physical causal constraint matrix, screening the causal compatible edge set, and applying the physical constraint maximum weight spanning tree algorithm, the constructed statistical model is guaranteed to conform to hydrophysical laws.

[0040] Step 201: Based on the watershed hydrological sequence, determine multiple risk factors and define the elements A of the physical causal constraint matrix. ij If risk factor i is a direct cause of risk factor j, then A ij =1, if risk factor j is the direct cause of risk factor i, then A ij =-1, otherwise A ij =0.

[0041] In this embodiment, the physical causal constraint matrix A is a d×d square matrix, where d is the total number of risk factors. This matrix is ​​used to explicitly encode expert knowledge judgments about causal relationships between variables. Matrix element A ij The rules for determining the value of strictly follow the physical mechanism: when a change in variable i necessarily leads to a change in variable j, and there is a clear physical transmission path (such as water flow evolution), i is determined to be the direct cause of j, denoted as A. ij =1; conversely, if j is the direct cause of i, let A be the cause of i. ij =-1; If there is no direct physical causal relationship between the two, or the causal relationship is unclear, then denote A as -1. ij =0.

[0042] The specific causal determination follows the following four priority rules: First, spatial causality, where upstream variables are the cause of downstream variables; second, process causality, where input variables at the same engineering node are the cause of output variables; third, state causality, where initial state variables (such as the initial water level) are the cause of process variables (such as the water level during the scheduling period); and fourth, external causality, where external boundary conditions (such as the water level of an external river) are the cause of internal variables of the system.

[0043] Step 202: Calculate the Kendall rank correlation coefficients between each pair of risk factors based on the watershed hydrological sequence, and construct the Kendall correlation coefficient matrix.

[0044] After determining the physical constraints, it is necessary to quantify the statistical correlation between the data. The Kendall rank correlation coefficient τ is a nonparametric statistic used to measure the strength of the dependency between two variables. It does not depend on the marginal distribution of the variables and is suitable for hydrological nonlinear variables. For any two risk factors X... i and X j The Kendall correlation coefficient is calculated using the following formula:

[0045] τ ij =(CD) / sqrt((C+D+n1)(C+D+n2));

[0046] Where C is the number of consistent pairs, D is the number of inconsistent pairs, and n1 and n2 are the adjustment amounts for the parallel groups. By iterating through all variable pairs, a d×d symmetric matrix T is calculated, where the elements T... ij =|τ ij | represents the absolute value of the correlation coefficient, characterizing the connection strength.

[0047] Step 203: Construct a causal compatible edge set. For any two risk factors i and j, if their corresponding A ij Not 0, or A ij If the Kendall rank correlation coefficient is 0 and the absolute value of the correlation coefficient exceeds the preset correlation threshold, then the edge connecting i and j will be added to the causal compatible edge set.

[0048] This step is crucial for fusing physical constraints with statistical data characteristics. Causal compatible edge set E c It is the candidate set of all possible edges. The selection rule reflects the principle of prioritizing physical properties while also considering statistical factors. If there is a clear causal relationship between the physical causal constraint matrix indicator variables i and j, i.e., A... ij If the edge (i,j) is equal to 1 or -1, regardless of its statistical relevance, it must be included in the candidate set to prevent the loss of important physical connections due to sampling errors in the data samples. If the physical matrix does not indicate a direct causal relationship, i.e., A... ij If |τ = 0, then we examine its statistical correlation. ij |Exceeded the preset significance threshold τ th (For example, using 0.1 or determining based on confidence tests) indicates that there may be an indirect connection or a correlation not captured by the physical matrix, thus adding the edge to the candidate set. Conversely, if there is neither physical causality nor a very weak statistical correlation, the edge is directly removed. This pre-screening mechanism narrows the space for subsequent structure search and eliminates physically meaningless noisy connections.

[0049] Step 204: Search for the maximum weight spanning tree within the causal compatible edge set to obtain the first-level tree structure of the VineCopula joint distribution model. For example... Figure 2 As shown, the specific process includes the following:

[0050] We construct a weighted graph on the causal compatible edge set using the absolute value of the Kendall rank correlation coefficient as the edge weight.

[0051] The initial tree node set contains any risk factor, and the tree edge set is empty;

[0052] Iteratively execute the following steps until the tree node set contains all risk factors: among the candidate edges connecting the tree node set and the non-tree node set, select the edge with the largest weight that belongs to the causal compatible edge set, and add the edge and the non-tree node it connects to to the tree edge set and the tree node set respectively.

[0053] Output the tree edge set as the first-level tree structure.

[0054] This step employs an improved Prim algorithm on the causal compatible edge set E. c Construct a maximum weighted spanning tree (MST). The specific process is as follows: Initialize the tree node set V. T Includes any one risk factor, such as variable 1, and tree edge set E. T Empty. Proceed to the iteration process: in all connected tree node sets V T AND non-tree node set V\V T Among the candidate edges, find the weights (i.e., |τ) ij |) The largest edge, the selected edge must belong to the causally compatible edge set E. c After finding the optimal edge (u,v) that satisfies the conditions, add the new node v to V. T Add edge (u,v) to E T Repeat this step until V. T It contains all d risk factors. The final output edge set E T The first-level tree structure T1 constitutes VineCopula.

[0055] Optionally, solving for the maximum weight spanning tree of the physical constraints includes the following process:

[0056] In the candidate edge set (causally compatible edge set E) c Construct a weighted graph G=(V,E) c ,W), where the node set V represents all risk factors, and the edge weight W ij =|τ ij |

[0057] An improved Prim algorithm is used to solve for the maximum weight spanning tree. The algorithm steps are as follows:

[0058] (1) Initialization: Select any node v0 and add it to the tree node set V T ={v0}, tree edge set E T ={}, where {} is an empty set;

[0059] (2) Iteration: on all connections V T With V\V T Among the candidate edges, select the edge (u,v) with the largest weight, where u∈V. T ,v∈V\V T Add v to V T Add (u,v) to E T The node set V represents all risk factors.

[0060] (3) Termination: When V T Stop when =V, output spanning tree T1=(V,E) T ).

[0061] The spanning tree T1 is the first-level tree structure of VineCopula, and also the skeleton structure of the Bayesian network.

[0062] This tree not only statistically maximizes the sum of correlations, but all its connections are also verified by the physical causal matrix. For each edge in the selected tree structure, the optimal PairCopula function is selected according to the AIC criterion, and its parameters are estimated using the maximum likelihood method.

[0063] Optionally, the selection and parameter estimation of the PairCopula function includes the following procedures:

[0064] Based on the first-level tree structure, for each edge (i,j)∈E T The optimal PairCopula function is selected. The candidate Copula function family includes: GaussianCopula, Student-tCopula, ClaytonCopula, GumbelCopula, FrankCopula, and their rotational forms.

[0065] The parameter θ of each candidate Copula is estimated using the maximum likelihood method:

[0066] ;

[0067] In the formula, c(·) is the Copula density function; F i (·) and F j (·) represents the marginal distribution functions of variables i and j, respectively; n is the number of samples.

[0068] Selecting the optimal Copula using the AIC criterion:

[0069] ;

[0070] In the formula, k is the number of parameters; The maximum likelihood value is used. The Copula function with the smallest AIC value is selected as the PairCopula for that edge.

[0071] Step 205: Determine the Bayesian network skeleton based on the first-layer tree structure, and calculate the causal score for each risk factor. The causal score is the number of times the risk factor is a cause in the physical causal constraint matrix minus the number of times it is a result. Orient each edge in the Bayesian network skeleton: if the A of the edge connecting risk factors i and j... ij =1, then the orientation is i->j; if A ij =-1, then the orientation is j->i; if A ij =0, then the risk factor with the higher causality score points to the risk factor with the lower causality score; the directed network is then subjected to directed acyclic detection and correction to obtain the initial risk factor Bayesian network, such as... Figure 3 As shown.

[0072] The first layer of VineCopula's tree structure is an undirected graph, while Bayesian networks are directed acyclic graphs (DAGs), thus the edge orientation problem must be solved. This step proposes a deterministic orientation strategy, calculating the causal score S for each variable. i The score is based on the physical causal constraint matrix A, calculated using the formula S. i =∑(A ik ==1)-∑(A ki ==1), which is the number of times variable i acts as a direct cause minus the number of times it acts as a direct result. A higher score indicates that the variable is more likely to be upstream (parent node) in the causal chain. Then, orient each edge (i,j) in the tree structure: if matrix A already explicitly specifies (A... ij If A ≠ 0, then it directly follows the physical rules, such as A ij =1, then the orientation is i pointing to j. If matrix A does not specify (A ij If the score is 0, then inference is made based on the causal score, with the direction defined as pointing from the node with the higher score to the node with the lower score (S). i >S j Then i->j). This rule avoids the arbitrariness of orientation in traditional methods. After initial orientation, a topology sorting algorithm is used to detect whether there are loops in the network. If a loop exists, the edge with the smallest statistical weight (Kendall coefficient) in the loop is identified, and its direction is reversed or removed from the Bayesian network structure (but retained in the Vine computation) until the DAG constraint is satisfied. The final network structure retains the strong correlations captured by VineCopula and has a clear physical causal flow direction.

[0073] Optionally, edge orientation based on the causal constraint matrix includes the following steps:

[0074] For each undirected edge (i,j)∈E in the first-level spanning tree T1 T Orient the data and generate a set of directed edges. :

[0075] ;

[0076] The calculation rules for the orientation function D(i,j) are as follows:

[0077] ;

[0078] In the formula, S i The causal score for variable i is defined as:

[0079] ;

[0080] This is calculated by subtracting the number of times a variable acts as a cause from the number of times it acts as a result. Variables with higher causal scores are more likely to act as parent nodes.

[0081] When D(i,j)=1, the directed edge is i->j; when D(i,j)=-1, the directed edge is j->i.

[0082] Optionally, the verification and correction of directed acyclicity includes the following steps:

[0083] For the oriented graph G D =(V, Perform directed acyclicity (DAG) verification:

[0084] (1) Using the topological sorting algorithm to detect graph G D Does it contain directed loops?

[0085] (2) If a loop C={v1->v2->...->v...} exists, k ->v1}, calculate the weight W of each edge in the loop. {vi,vi+1} ;

[0086] (3) Remove the edge with the smallest weight (v) in the loop. m ,v m+1 ): ;

[0087] (4) Repeat steps (1)-(3) until there are no loops in the diagram;

[0088] (5) Output the risk factor Bayesian network skeleton structure G that satisfies the DAG constraint BN =(V, ).

[0089] Optionally, the iterative construction of the high-level Vine structure includes the following steps:

[0090] Based on the PairCopula function of the first-level tree, the second to (d-1)th level trees are constructed sequentially according to the R-Vine construction rules:

[0091] (1) For a tree of level l (l≥2), the node set is the edge set of the tree of level (l-1);

[0092] (2) Two nodes are adjacent in the l-th level tree if and only if they share a node in the (l-1)-th level tree;

[0093] (3) Model each edge using conditional Copula, and use the same parameter estimation method as step - PairCopula function selection and parameter estimation;

[0094] (4) Iterate until the construction of the (d-1)th level tree is completed.

[0095] Finally, the complete R-VineCopula model and the Bayesian network skeleton structure of risk factors were obtained.

[0096] By introducing a physical causal constraint matrix and a constrained maximum spanning tree algorithm, this invention achieves the following: First, the obtained first-layer tree structure of VineCopula conforms to the physical causal laws of the reservoir-lake system, improving the interpretability of the model. Second, the introduction of physical constraints reduces the structure learning error of pure data-driven methods under small sample conditions, improving the robustness of the model. Third, the pre-screening of the candidate edge set reduces the size of the structure search space, improving the computational efficiency of the algorithm. Fourth, it lays a physical foundation for the edge orientation of the Bayesian network, ensuring consistency between the VineCopula structure and the Bayesian network structure.

[0097] By introducing a deterministic edge orientation method based on physical causal constraints, this invention achieves the following: First, it solves the problem of unclear rules for transforming undirected graphs into directed graphs, providing a complete and operable edge orientation algorithm. Second, the orientation of all edges has clear physical basis or mathematical rules, eliminating the arbitrariness and uncertainty of orientation. Third, the orientation results conform to the physical causal laws of the reservoir-lake system, giving the Bayesian network clear physical interpretability. Fourth, through loop detection and correction mechanisms, it ensures that the output results strictly satisfy the directed acyclic graph constraints, and can be directly used for probabilistic inference in Bayesian networks. Fifth, it achieves the collaborative construction of the VineCopula structure and the Bayesian network structure, both based on the same set of edge orientation rules, ensuring structural consistency.

[0098] Example 3 describes the adaptive risk sampling and scheduling-risk closed-loop iterative process, aiming to solve how to effectively capture black swan scenarios with extremely low probability of occurrence but huge destructive power in an infinite space of risk factor combinations, and to transform static risk analysis into a dynamic exploration-utilization process by constructing mathematical risk weights and closed-loop feedback mechanisms.

[0099] Step 301: Calculate the posterior probability that each sample point will cause the system risk index to exceed the preset risk threshold, and obtain the risk-sensitive weight function, which is calculated using the following formula:

[0100] w k (x)=P(R>R th |X=x;BN k );

[0101] Among them, w k (x) is the risk-sensitive weight function in the k-th iteration, x is the risk factor vector, and R is the system risk index. th To preset risk thresholds, BN k This is the Bayesian network of risk factors for the current round.

[0102] In this embodiment, the risk-sensitive weighting function w k (x) serves as a bridge connecting the statistical space and the risk space. Specifically, the risk factor vector x is represented as (x1, x2, ..., x...). d ), where d is the dimension, and the components include, but are not limited to, the peak flow and flood volume of each reservoir, the inflow between lakes, and the initial regulating water level of each project. The system risk index R is specifically selected as the highest water level Z of the lake during the entire scheduling period T. max The preset high-risk threshold R th It can be set to the lake's guaranteed water level (e.g., 12.5m) or the highest disaster water level in history.

[0103] In the k-th iteration, the Bayesian network BN k Having already learned parameters using data from the first k-1 rounds, it possesses a preliminary risk reasoning ability. Calculate w. k The process of (x) is also the forward reasoning process of a Bayesian network: inputting specific risk factor sample values ​​x into the network as evidence, and using variable elimination or cluster tree propagation algorithms, deriving that the values ​​of leaf nodes (i.e., risk indicator nodes) fall within the interval (R). th The probability of (+∞). For example, for a combination x0 (high water level upstream, high water level downstream), if BN k The probability that it will cause the lake water level to exceed the standard is predicted to be 0.9, then w k (x0) = 0.9; For a safe combination x1 (dry year), the predicted probability of exceeding the standard is 0.001, then wk (x1) = 0.001. Through the above method, the risk-sensitive weight function w k (x) assigns extremely high weights to high-risk areas while suppressing the weights of safe areas.

[0104] Step 302: Construct the risk-weighted sampling density using the following formula:

[0105] f k (x)∝w k (x)·c(u1,...,u d ;θ k )·∏ i=1 d f i (x i );

[0106] Among them, f k (x) represents the risk-weighted sampling density of the k-th iteration, c(·) is the density function of the VineCopula joint distribution model, and θ k For the current Copula parameters, u i Risk factor x i The marginal cumulative distribution function value, where d is the dimension of the risk factor, f i (x i ) represents risk factor x i The marginal probability density function; generating a risk-adaptive scenario set, including: based on the risk-weighted sampling density f k (x) A multidimensional risk factor sample is generated using an improved Latin hypercube sampling method.

[0107] In this embodiment, the risk-weighted sampling density f k (x) is an artificially distorted probability distribution. The formula contains c(·)·∏f i (x i The risk factor represents the joint probability density of the risk factors that objectively exist in nature, i.e., the prior distribution, while the first term, the risk-sensitive weight function w, represents the risk factor. k (x) is a correction factor, i.e., risk appetite. The product of the two makes f k (x) preserves the physical correlation structure between variables (controlled by Copula) while concentrating the probability quality in the high-risk region.

[0108] To efficiently sample from non-standard distributions, this embodiment employs an improved Latin hypercube sampling (LHS) method. The specific implementation process includes: based on the risk-weighted sampling density f... k The marginal distribution characteristics of (x) divide the value range of each dimension into M equally probable intervals, due to the introduction of a risk-sensitive weight function w. k(x) These intervals are not of equal width in the original physical space, with higher-risk areas having denser intervals. A value is randomly selected from each interval in each dimension, and simulated annealing is used to pair (shuffle) the sampled values ​​across dimensions to minimize the correlation error between sampled points, making it as close as possible to the risk-weighted sampling density f. k The multidimensional correlation structure defined by (x). The final risk adaptive scenario set X is generated. new It contains M samples, for example, M=1000, and the samples mathematically follow the formula f. k The (x) distribution physically represents the most concerning potential disaster-causing flood combinations based on current understanding.

[0109] In some alternative implementations, to avoid overly concentrated sampling leading to a decrease in exploration capability, a risk-weighted sampling density f can be constructed. k When (x) is expressed, a smoothing coefficient α is introduced, and the formula is modified to f. k (x)∝(w k (x)) α c(...), where 0 < α < 1. When α is small, the sampling distribution is close to the prior distribution, focusing on global search; when α is close to 1, the sampling distribution is highly concentrated in high-risk areas, focusing on local utilization.

[0110] Step 303: Incorporate the samples from the risk adaptive scenario set and their corresponding optimal scheduling scheme results into the historical training sample set; based on the updated training sample set, re-estimate the conditional probability table parameters of the risk factor Bayesian network using the EM algorithm to obtain the updated risk factor Bayesian network; calculate the difference between the risk estimates in the high-risk interval between the current iteration and the previous iteration. If the difference is less than the preset convergence threshold or the number of iterations reaches the preset upper limit, the risk assessment result is considered converged. Figure 4 As shown.

[0111] This step details the parameter update mechanism of the closed-loop feedback. For each risk adaptive scenario x generated in step 302... j (j=1...M), calculate the highest water level Z of the corresponding lake using the reservoir group scheduling model. maxj To form a complete data pair (x j Z maxjNew data is appended to the sample library from the previous round, gradually increasing the proportion of high-risk samples. The Expectation-Maximization (EM) algorithm is used to relearn the Bayesian network parameters. The EM algorithm consists of two steps: the expectation step (E-step), which fills in any missing data or calculates the expectation of latent variables using the current parameter estimates; and the maximization step (M-step), which updates the conditional probability table (CPT) of each node in the network by maximizing the log-likelihood function. For example, for the node representing lake water level, the EM algorithm adjusts its conditional probability distribution based on the new samples when the parent node (e.g., inflow rate) takes different values, allowing the network to more accurately reflect the water level response characteristics under extreme conditions. To determine whether the iteration has terminated, the convergence index ε = ||P is calculated. k (R>R th )-P k-1 (R>R th ε|| represents the difference between the overall system failure probability calculated in the two iterations. When ε is less than a preset threshold (e.g., 10), the probability is considered to be 10. -4 When the network's perception of risk has stabilized, the iteration ends.

[0112] Example 4 provides an alternative implementation scheme for constructing a VineCopula joint distribution model and an initial risk factor Bayesian network that satisfy physical causal topological constraints. Unlike the rule-based construction method in Example 2, this example adopts a score-search-based optimization method, which is suitable for scenarios with complex physical constraints or soft constraints.

[0113] Step 401: Construct a VineCopula joint distribution model and an initial risk factor Bayesian network that satisfy physical causal topological constraints. A score-search-based structure learning method is used, with the objective function being:

[0114] max G L(G)-λ1Φ(G,G vine )-λ2Ψ(G,C phy );

[0115] Where G is the Bayesian network structure to be searched, L(G) is the log-likelihood function, and Φ(G,G) is the maximum likelihood function. vine () represents the network structure G and the first-level tree structure G of VineCopula. vine Inconsistency penalty term, Ψ(G,C) phy The network structure G violates the physical causality constraint matrix C. phy The physical topology penalty term is used, and λ1 and λ2 are weight coefficients. During the search process, the directed acyclic graph that maximizes the objective function is selected as the initial risk factor Bayesian network.

[0116] In this embodiment, the determination of the Bayesian network structure is transformed into a mathematical optimization problem. The objective function consists of three parts: the first part, L(G), is the Bayesian Information Criterion (BIC) or log-likelihood score, used to measure the degree of fit of the network structure G to the observed data; a larger value indicates a better statistical fit; the second part, Φ(G,G...)... vine ) is a structural consistency constraint used to penalize G for deviating from the first-level tree structure G of VineCopula. vine The difference between them. Specifically, if G vine If an edge (i,j) exists in G but does not exist in G (regardless of direction), then Φ increases by 1; this term is to ensure that the Bayesian network inherits the strongly correlated skeleton captured by the Copula model. The third part Ψ(G,C) phy ) is a physical topology penalty term used to penalize violations of the physical causality constraint matrix C. phy The behavior. The specific calculation rule is: traverse each directed edge i->j in G, if the physical matrix specifies A... ji =1, meaning j is the cause of i, the direction should be j->i, or A ij If Ψ = 0 and physical connections are prohibited, then Ψ will be increased by a larger penalty value, such as 100. λ1 is used to control the weight of statistical consistency, and the value range can be selected from 0.5 to 2.0; λ2 is used to control the weight of physical compliance. Since physical constraints have higher priority, the value range can be selected from 5.0 to 20.0. The weight coefficients λ1 and λ2 (for example, values ​​of 1.0 and 10.0) are used to balance the relationship between data fitting, statistical consistency and physical compliance.

[0117] To solve this optimization problem, hill climbing or simulated annealing algorithms can be used. The search process begins with a random structure or an empty graph. In each iteration, it attempts to add, remove, or transform edges on G, calculating the change in the score of the new structure. If the score increases, the operation is accepted; otherwise, it is rejected with a certain probability. During the search process, a tabu list is used to avoid loop searches, and the DAG property is always checked. The final converged structure G is the optimal initial risk factor Bayesian network, which fits the data to the greatest extent possible while taking into account the correlation structure of VineCopula and the causal logic of the physical system.

[0118] Example 5 describes the solution of the optimal scheduling scheme for a reservoir group for a risk-adaptive scenario set.

[0119] In this embodiment, for each generated risk adaptive scenario, the inflow process Q of each reservoir is included. i (t) and the starting water level Z i (0) A mathematical model for joint optimal scheduling of reservoir groups and lakes is established. The objective function of this model is set as minimizing the average remaining excess flood volume of lakes during the scheduling period, and its expression is:

[0120] MinF=(1 / T)∑ t=1 T [max(0,V lake (t)-V safe )];

[0121] Where T is the total number of scheduling periods, and V lake (t) represents the water volume of the lake at time t, V safe This represents the safe water storage capacity of the lake (corresponding to the guaranteed water level). The physical meaning of this objective function is to utilize the impounding effect of the upstream reservoir group as much as possible to reduce the time and volume of water that the downstream lake maintains at a high water level.

[0122] The constraints of the model include:

[0123] (1) Water balance equation: V i (t+1)=V i (t)+(Q in_i (t)-Q out_i (t))Δt, this equation holds for all reservoirs and lakes;

[0124] (2) Water level constraint: Z min_i <=Z i (t)<=Z max_i Ensure that the water level is between the minimum and maximum allowable water levels;

[0125] (3) Discharge capacity constraint: Q out_i (t)<=Q capacity (Z i (t)), that is, the discharge flow rate cannot exceed the maximum discharge capacity at the current water level;

[0126] (4) Amplitude constraint: |Q out_i (t)-Q out_i (t-1)|<=ΔQ limit This is to ensure the feasibility and operability of the scheduling plan.

[0127] Among them, Z min_i Z represents the minimum allowable water level of the i-th reservoir. i Z(t) represents the water level of the i-th reservoir / lake at time t. max_i Let Q be the maximum allowable water level of the i-th reservoir. capacity (Z i (t) represents the current water level Z. i The maximum discharge capacity under (t), ΔQ limit This represents the maximum allowable flow rate variation.

[0128] To solve this non - linear and multi - constraint optimization model, the sparrow search algorithm (SSA) is adopted in this embodiment. This algorithm simulates the foraging and anti - predation behaviors of a sparrow population and features a fast convergence speed and strong robustness. During the solution process, the downstream discharge of each reservoir in T time periods is used as the optimization variable (individual position).

[0129] In the algorithm iteration, the position update formula for the discoverer (Producer) is as follows:

[0130] X i,j t+1 =X i,j t exp(-i / (α*iter max )) (when R2 < ST);

[0131] Among them, X i,j t represents the position of the i - th sparrow in the j - th dimension at the t - th iteration, α is a random number, iter max is the maximum number of iterations, R2 is the warning value, and ST is the safety threshold.

[0132] The joiner (Scrounger) then follows the discoverer for position update. Through multiple iterations, such as 500 times, the discharge sequence represented by the optimal individual in the population is the optimal scheduling scheme under this scenario. This scheme minimizes the flood control risk of the downstream lake and provides a real sample reflecting the optimal flood - fighting ability of the system for subsequent Bayesian network parameter update.

[0133] Embodiment 6 describes how to perform structural expansion based on the converged risk - factor Bayesian network and achieve a quantitative assessment of the system's comprehensive flood control risk by combining probability inference techniques, solving the problem that relying solely on risk factors cannot directly map to system consequences, and establishing a complete causal chain from hazard - causing factors to disaster consequences by introducing intermediate process variable nodes.

[0134] Step 601: Define intermediate process variables, which at least include the maximum downstream discharge during the reservoir scheduling period, the highest water level during the reservoir scheduling period, and the total inflow flood volume into the lake.

[0135] In this embodiment, to more accurately describe the risk transfer process, intermediate process variables need to be introduced based on the original risk - factor Bayesian network. These variables are not the original inputs of the system (i.e., non - risk factors), but process state variables with clear physical meanings calculated by the scheduling model. Specifically, the following three types of intermediate variables are defined: The first type is the reservoir state variable. For example, the maximum downstream discharge q max_i and the highest water level Z max_iThe first category directly reflects the reservoir's regulation and storage behavior and its own risks; the second category is transitive variables, such as the total peak flow Q evolving to the lake inlet section. _in_total The third category is convergent variables, such as the total flood volume W flowing into the lake during the scheduling period. h Intermediate variables act as a bridge connecting rainfall / inflow (source) and lake level (consequence), allowing the network structure to explicitly express the physical mechanism by which a major flood leads to increased reservoir discharge, resulting in a rise in lake level.

[0136] Step 602: Based on the converged risk factor Bayesian network, add directed edges to expand the network structure according to the following rules: If the risk factor is the input of the scheduling model and directly affects an intermediate process variable, add a directed edge from the risk factor to the intermediate process variable; if there is a hydraulic transmission relationship between two intermediate process variables, add a directed edge from the upstream variable to the downstream variable; if the intermediate process variable directly affects the system risk index, add a directed edge from the intermediate process variable to the system risk index.

[0137] In this embodiment, the network expansion process no longer relies on simple data learning, but is based on deterministic physical and scheduling logic. The specific expansion rules are implemented as follows: For the scheduling input edge, due to the reservoir inflow peak Q... pi (Risk factors) directly determine the flood control scheduling decisions of the reservoir and affect the maximum outflow q of the reservoir. max_i (Intermediate variable), therefore add a directed edge Q. pi ->q max_i Regarding the hydraulic transmission edge, considering that the discharge from the upstream reservoir needs to go through the river channel to reach the downstream lake, a discharge from the maximum discharge q of the reservoir is added. max_i The total flood volume W flowing into the lake h The directed edge. Regarding the risk convergence edge, due to the highest water level Z of the lake... max (System risk indicators) mainly depend on the total flood volume W flowing into the lake. h And the lake's water level adjustment Z 0h (Risk factor), therefore a directed edge W is added. h ->Z max and Z 0h ->Z max Through hierarchical expansion, a complete Bayesian network comprising a risk factor layer, an intermediate process layer, and a risk indicator layer was constructed. Using the risk-adaptive scenario set and its corresponding scheduling results generated in the aforementioned embodiments as training data, the conditional probability table (CPT) of newly added nodes was learned using either maximum likelihood estimation (MLE) or expectation-maximization (EM) algorithms, thus completing the parameterization of the expanded network.

[0138] Step 603: Perform probabilistic reasoning based on the extended Bayesian network and output the comprehensive flood control risk assessment results of the system.

[0139] In this embodiment, the comprehensive risk assessment includes three dimensions: positive prediction, reverse diagnosis, and overall risk calculation.

[0140] Predict risk probability using forward inference. Given a set of observed risk factor evidence e, such as the inflow peak of Reservoir 1 falling within the interval [5000, 6000]m. 3 / s, input the evidence into the Bayesian network to calculate the target node Z. max Posterior probability distribution P(Z) in different water level ranges max |e) can help decision-makers predict the flood control situation under the current water conditions.

[0141] The reverse diagnostic function is used to identify the causes of the disaster. Evidence e' is set as the highest water level Z of the lake. max If the water level exceeds the guaranteed level, the posterior probability distribution of each risk factor X is inferred backwards using the Bayesian formula P(X|e')=P(e'|X)P(X) / P(e'). For example, the calculation results may show that when the lake exceeds the guaranteed level, the probability of the reservoir's water level being too high increases by 30%, identifying this factor as the main cause of this risk.

[0142] Calculate the overall flood control risk of the system. Since the Copula model describes the joint encounter probability of risk factors, while Bayesian networks describe the system response probability under a given combination of factors, the overall failure probability P of the system is calculated. fail Both need to be combined. The specific calculation formula is as follows:

[0143] P fail =∑ j=1 M P copula (x j ∈Ω j )P BN (R>R th |x j ∈Ω j );

[0144] The entire risk factor space is discretized into M hypercube intervals Ω. j P copula (·) represents the combination of factors calculated using the VineCopula model that falls within the interval Ω. j The joint probability (encounter probability) within, P BN (·) is calculated using an extended Bayesian network when the factor combination is at Ω. j The conditional probability of system risk indicators exceeding the limit is given by this formula. This formula achieves the mathematical coupling of hydrological uncertainty (characterized by Copula) and hydraulic / dispatch uncertainty (characterized by BN), resulting in a more comprehensive and accurate assessment than traditional single methods.

[0145] Example 7 provides a data-driven VineCopula construction and non-adaptive risk analysis method as an alternative to Examples 2 and 3. This example uses a purely data-driven statistical learning method to construct the model and employs static sampling for risk analysis. It is suitable for application scenarios where the physical mechanism is not yet clear, it is difficult to construct a physical causal constraint matrix, or where high computational efficiency is required while the accuracy requirement for extreme risks is relatively low.

[0146] Step 701: Construct the VineCopula model structure using a purely data-driven maximum spanning tree algorithm.

[0147] In this embodiment, instead of introducing a physical causal constraint matrix A, the model is directly constructed based on the statistical correlation between risk factors. Specifically, the Kendall rank correlation coefficients between each pair of risk factors are calculated, and the absolute values ​​of the correlation coefficients are used as edge weights to construct a complete graph. The standard Prim algorithm or Kruskal algorithm is directly applied to search for the maximum spanning tree (MST) on the complete graph, which serves as the first-level tree structure of VineCopula.

[0148] After determining the tree structure, the AIC criterion is used to select the PairCopula function. For constructing Bayesian networks, structure learning algorithms such as the K2 algorithm or PC algorithm can be used, based on data mining of dependencies between variables, or experts can directly specify the network structure based on experience without imposing physical rules. This method can be automated, eliminating the need for pre-planning complex physical topologies, but with a small sample size, it may learn spurious correlations that violate physical common sense.

[0149] Step 702: Generate a risk scenario set based on static Latin hypercube sampling.

[0150] In this embodiment, a risk-sensitive weighting function w is not introduced. k (x) and closed-loop iteration mechanism. After constructing the VineCopula joint distribution model, the conventional Latin hypercube sampling (LHS) method is directly used to generate N risk factor samples. Specifically, the marginal distribution intervals of each variable are divided into N equal parts with equal probability. Random sampling is performed within each sub-interval, and the samples are arranged and combined according to the correlation structure of the Copula. The generated sample set distribution conforms to historical statistical patterns, i.e., it has a high probability of generating common small and medium-sized flood scenarios and a low probability of generating extreme flood scenarios.

[0151] Step 703: Conduct a risk assessment based on a single scheduling simulation. Input the generated N sets of static samples into the reservoir group scheduling model to calculate N sets of lake water level results. Directly based on the N sets of input and output data, calculate the frequency of lake water levels exceeding the threshold, using this as an estimate of the system failure probability: P fail ≈(N fail / N); where N fail This represents the number of samples that caused the failure. Alternatively, a static Bayesian network can be trained using N sets of data for subsequent inference analysis. This alternative approach is simple and has a one-shot computational cost, but when dealing with black swan events, the variance in the estimation of tail risk is large because the N samples may only contain a very few extreme scenarios, resulting in lower accuracy than the adaptive iterative method.

[0152] The data-driven method in this embodiment can serve as a rapid pre-evaluation module within the system of this invention. In practical applications, the rapid mode of this embodiment can be run first to obtain a risk overview. If a high potential risk is found, the system can then be switched to the physical constraint adaptive mode for refined analysis.

[0153] Example 8: Taking a typical two-reservoir-one-lake flood control system in the Yangtze River Basin as an example, the system includes two large reservoirs connected in parallel upstream and a large lake connected to the Yangtze River downstream to receive incoming water.

[0154] The flood control system described in this embodiment consists of Reservoir A, Reservoir B, and Lake C. Reservoir A is located on a northwestern tributary of the basin, controlling a drainage area of ​​approximately 15,000 square kilometers. Reservoir B is located on a northeastern tributary of the basin, controlling a drainage area of ​​approximately 12,000 square kilometers. The outflows from both reservoirs flow into Lake C after evolving through their respective downstream channels. Lake C also receives water from the uncontrolled area between the two reservoirs. The lake's outlet is connected to the Yangtze River main stream via a sluice gate, and its outflow capacity is affected by the backwater effect of the Yangtze River's main stream water level.

[0155] Based on the above topology, seven risk factors are defined. Risk factor 1 is the peak inflow discharge into Reservoir A, denoted as Q. p1 The unit is m 3 / s. Risk factor number 2 is the peak inflow discharge of Reservoir B, denoted as Q. p2 The unit is m 3 / s. Risk factor number 3 is the peak flow of the lake in section C, denoted as Q. ph The unit is m 3 / s. Risk factor number 4 is the starting water level of reservoir A, denoted as Z. 01 The unit is meters. Risk factor number 5 is the starting water level of reservoir B, denoted as Z. 02 The unit is meters. Risk factor number 6 is the water level of Lake C, denoted as Z. 0hThe unit is meters. Risk factor number 7 is the highest water level during the Yangtze River main stream's scheduling period, denoted as Z. w The unit is meters.

[0156] The construction of the physical causal constraint matrix A follows four priority rules, which are in descending order of priority: mandatory causal rules, spatial transit rules, state influence rules, and external boundary rules.

[0157] The first type is the mandatory causal rule, applicable to variable pairs with direct hydraulic connections and undisputed causal directions. When a change in the value of variable i necessarily and directly leads to a change in the value of variable j, and this transitive relationship is governed by the fundamental equations of hydraulics, let A... ij =1. For example, the peak inflow of water into reservoir A is Q. p1 The scheduling behavior of Reservoir A is directly determined by the reservoir water balance equation, which in turn affects the outflow. However, it should be noted that this matrix only describes the relationship between the original risk factors; the reservoir outflow is an intermediate process variable and will be handled in subsequent network expansion.

[0158] The second type is the spatial transfer rule, applicable to upstream and downstream variable pairs with river hydraulic connections. When the outflow from an upstream node becomes the inflow from a downstream node after evolving through the river channel, the upstream variable is the cause of the downstream variable. In this embodiment, the outflow from Reservoir A evolves into Lake C, and the outflow from Reservoir B also evolves into Lake C. Since the reservoir outflow is driven by the inflow, from the perspective of the original factors, Q... p1 and Q p2 All are related to the inflow from Lake C. However, the peak flood Q in the lake area... ph The runoff mainly originates from rainfall in the uncontrolled area between the two reservoirs, and is related to Q. p1 and Q p2 The relationship needs to be analyzed specifically. If the three flood peaks are caused by the same rainfall event, then there is a common driving factor rather than a direct causal relationship. In this case, A ij Set to 0. If the flood peak in the interval is mainly replenished by the lateral overflow of the upstream reservoir in the river channel, then a causal relationship exists, A. ij Set to 1. This embodiment assumes that the interval watershed generates runoff independently, therefore A 13 and A 23 Set all to 0.

[0159] The third category is the state influence rule, which applies to the relationship between initial state variables and process variables. The initial water level, as the initial condition for scheduling calculations, affects the water level change process throughout the entire scheduling period. In this embodiment, the initial water level Z of reservoir A... 01 The impact is on the scheduling process of Reservoir A itself, not on other risk factors. Since this matrix describes the relationships between the original risk factors, there is no causal relationship between the reservoir's initial water level and its inflow peak; therefore, A... 14 and A41 All are set to 0. Lake initial water level Z 0h The same rule applies.

[0160] The fourth category is external boundary rules, which apply to the influence of external boundary conditions on internal variables of the system. (Example: Yangtze River main stream water level Z) w As an external boundary condition, Z affects the outflow capacity of lake C through its backwater effect, thereby influencing the water level process during the lake's regulation period. From a causal perspective, Z... w The lake level is the result, while the initial water level of the lake is the cause. 0h It precedes the scheduling process in time, and Z w The description refers to the highest water level during the scheduling period; the two overlap in time series and require further analysis. If Z... w If the predicted future water level process does not affect the already determined starting water level, A 76 Set Z to 0. w If the real-time water level represents the water level synchronized with the current scheduling period, then there is an interaction between the two rather than a one-way causal relationship. This situation will be handled according to the dispute resolution rules.

[0161] The rules for handling disputes are now outlined. In the actual matrix construction process, the causal relationships between some variable pairs may be disputed or uncertain, and these need to be handled according to the following rules.

[0162] Dispute Resolution Rule 1: Handling Ambiguous Timing. When the temporal relationship between two variables is unclear, analyze the dominant direction of their physical processes. If the change in variable i precedes the response of variable j in most scenarios, then i is determined to be the cause of j. Let A... ij =1. If the changes of the two occur almost simultaneously and influence each other, then it is determined that there is no direct causal relationship. Let A be an example. ij =0. In this embodiment, the lake's initial water level Z is... 0h With the Yangtze River water level Z w All are variables related to the scheduling period, but the initial water level is defined as the state at the start of the scheduling, while Z... w Defined as the highest value during the scheduling period, with the start time of the scheduling occurring first; therefore, if there is any influence, it should be Z. 0h Influences subsequent processes on Z w The response, but this goes beyond the direct causality at the level of the original factors, therefore let A... 67 and A 76 All are 0.

[0163] Dispute Resolution Rule Two: Distinguishing Between Indirect and Direct Causation. If variable i affects variable j through intermediate variable m, i.e., the causal chain is i to m to j, then there is no direct causal relationship between i and j. Let A... ij=0. The intermediate variable m will be processed during the Bayesian network expansion stage. In this embodiment, the peak inflow flood Q into the reservoir is... p1 The reservoir's outflow is affected by its scheduling, and this outflow, in turn, influences the lake's inflow through river channel evolution. Here, both the reservoir outflow and the river's inflow are intermediate variables; therefore, at the original factor level, Q... p1 With Z 0h Assuming there is no direct causal relationship between A and B, 16 =0.

[0164] Dispute Resolution Rule Three: Handling of Common Driving Factors. If two variables are simultaneously driven by a third factor and exhibit statistical correlation, but there is no direct physical transmission relationship between them, then let A... ij =0. For example, the peak inflow of floodwater into reservoir A is Q. p1 The flood peak Q entering Reservoir B p2 They may be strongly correlated due to the same rainstorm, but the water from reservoir A will not physically flow to reservoir B, therefore A... 12 and A 21 Both are set to 0. Similarly, the starting water level Z of the two reservoirs is... 01 and Z 02 They may exhibit correlation due to unified scheduling procedures, but they are not hydraulically connected, therefore A 45 and A 54 Set all to 0.

[0165] Rule Four for Dispute Resolution: Handling Weak Causal Relationships. When a causal relationship exists physically but is weak or only manifests under extreme conditions, its inclusion in the matrix can be chosen based on the research objective. If the objective is to capture extreme risks, it is recommended to include it and set A. ij =1. If the objective is routine risk analysis, this can be ignored and set A. ij =0. This embodiment aims at extreme risk analysis and tends to incorporate causal constraints for boundary cases.

[0166] Rule Five for Dispute Resolution: Handling Disagreements Among Experts. When multiple experts in the field disagree on the causal relationship of a pair of variables, the majority vote principle applies. If more than half of the experts support the existence of a causal relationship, then let A... ij Let A be either 1 or -1. If the opinions are similar, then let A be... ij =0, leaving the association of the variable pair to data-driven statistical learning for judgment.

[0167] Based on the above rules, the complete physical causal constraint matrix A of this embodiment is now given. The matrix is ​​7 rows and 7 columns, with rows and columns arranged according to risk factor numbers 1 to 7.

[0168] The first row of the matrix corresponds to risk factor Q. p1 The elements are as follows: A 11=0 indicates that A has no causal relationship with itself. 12 =0 indicates Q p1 With Q p2 In cases where there is no direct causal relationship but a shared driving force, A 13 =0 indicates Q p1 With Q ph A is an independent runoff region with no direct causal relationship. 14 =0 indicates Q p1 With Z 01 No causal relationship, A 15 =0 indicates Q p1 With Z 02 Irrelevant, A 16 =0 indicates Q p1 With Z 0h There is no direct causal relationship; the intermediate variable must be passed on. A 17 =0 indicates Q p1 With Z w There is no direct relationship.

[0169] The second row of the matrix corresponds to the risk factor Q. p2 The elements are as follows: A 21 =0, A 22 =0, A 23 =0, A 24 =0, A 25 =0, A 26 =0, A 27 =0. The reason for the value of the element in this row is similar to that of the first row.

[0170] The third row of the matrix corresponds to risk factor Q. ph The elements are as follows: A 31 =0, A 32 =0, A 33 =0, A 34 =0, A 35 =0, A 36 =0 indicates Q ph With Z 0h There is no direct causal relationship; the influence of inter-regional flow on lakes requires an evolutionary process. (A) 37 =0.

[0171] The fourth row of the matrix corresponds to risk factor Z. 01 The elements are as follows: A 41 =0, A 42 =0, A 43 =0, A 44 =0, A 45 =0, A 46 =0, A 47 =0.

[0172] The fifth row of the matrix corresponds to risk factor Z.02 The elements are as follows: A 51 =0, A 52 =0, A 53 =0, A 54 =0, A 55 =0, A 56 =0, A 57 =0.

[0173] The sixth row of the matrix corresponds to risk factor Z. 0h The elements are as follows: A 61 =0, A 62 =0, A 63 =0, A 64 =0, A 65 =0, A 66 =0, A 67 =0.

[0174] The seventh row of the matrix corresponds to risk factor Z. w The elements are as follows: A 71 =0, A 72 =0, A 73 =0, A 74 =0, A 75 =0, A 76 =1 indicates Z w As an external boundary condition, it affects the lake's water level response during the scheduling period, thus impacting risk analysis in relation to Z. 0h There is interaction, A 77 =0.

[0175] It should be noted that the matrix exhibits a large number of zero elements because this embodiment strictly distinguishes between direct causal relationships at the level of the original risk factors, excluding a large number of indirect relationships transmitted through intermediate process variables. These indirect relationships will be characterized by adding intermediate nodes during the Bayesian network expansion stage.

[0176] In certain applications, if researchers wish to directly characterize end-to-end causal chains at the level of raw factors, a loosely defined causal matrix can be used. For example, if it is believed that the inflow peak of an upstream reservoir will ultimately affect the water level of a downstream lake, then A can be directly defined as... 16 =1 and A 26 =1. The above processing method simplifies the network structure, but sacrifices the ability to finely characterize intermediate processes. The specific definition method to be adopted should be determined comprehensively based on the research accuracy requirements and data availability.

[0177] Furthermore, the physical causal constraint matrix is ​​not static. When changes occur in the river basin system, such as the commissioning of a new reservoir or the completion of river straightening projects, or when extreme events reveal previously unidentified risk transmission paths, the matrix elements should be updated promptly. It is recommended that the river basin management agency organize expert consultations before the start of each flood season to review and revise the matrix as necessary.

[0178] For example, the numerical values ​​of intermediate results such as the Kendall coefficient matrix, maximum spanning tree construction, and weight calculation are shown below:

[0179] Assuming we have obtained the following Kendall rank correlation coefficient matrix T based on historical flood data (only a portion is shown):

[0180] τ(Q p1 Q p2 ) = 0.65;

[0181] τ(Q p1 Q ph )=0.42;

[0182] τ(Q p1 Z 0h ) = 0.35; ...

[0184] Based on this matrix and the physical causal constraint matrix, a causal compatible edge set is constructed, containing edges: (1,2), (1,3), (2,3), (6,7)....

[0185] Applying Prim's algorithm, initially selecting node 1, the iterative process is as follows:

[0186] First iteration: Select edge (1,2), weight 0.65, and add node 2;

[0187] Second iteration: Select edge (1,3), weight 0.42, and add node 3; ...

[0189] The first-level tree structure is obtained as follows: {(1,2),(1,3),(3,6),(6,7),(4,1),(5,2)}.

[0190] According to one aspect of this application, the position update of the Producer has been described in Embodiment 5. To fully implement the SSA algorithm, update rules for the Scrounger and Alerter / Vigilant also need to be defined.

[0191] For the joiner Scrounger, the position update formula is:

[0192] When i > n / 2, Xij t+1 =Q·exp((X worst t -X ij t ) / i 2 );

[0193] When i <= n / 2, X ij t+1 =X P t+1 +|X ij t -X P t+1 |A + L.

[0194] Among them, X ij t+1 X is the updated position of the i-th individual in the (t+1)-th generation and the j-th dimension, Q is a random number following a normal distribution, and X... worst t X is the position of the worst individual (with the worst fitness) in generation t. ij t X is the position of the i-th joiner in the t-th generation on the j-th dimension. P t+1 It is the position of the discoverer (Producer) in generation t+1 (usually the current best individual, i.e., the global best solution), A + =A T (AA T ) -1 A is a 1×d' matrix with elements of 1 or -1, and L is a matrix of all 1s. This indicates that when a newcomer is in the latter half of the population (i>n / 2), its fitness is poor, and it needs to fly to other places to forage; otherwise, it will search around the vicinity of the discoverer.

[0195] For the Alerter, its position update formula is:

[0196] ;

[0197] in, It represents the updated position of the i-th individual in the (t+1)-th generation and the j-th dimension. It is the position of the globally optimal individual in the t-th generation in the j-th dimension. β is the position component of the i-th individual in the t-th generation and j-th dimension, and β is the step size control parameter (a random number following a normal distribution). f is the position of the worst individual in the t-th generation in the j-th dimension, K is a random number in [-1, 1], and f i f g f wThese are the fitness values ​​for the current individual, the globally optimal, and the globally worst, respectively, with ε being a minimal constant. This formula simulates the behavior of a sparrow moving to a safe area (optimal position) or making an emergency adjustment near its current position when it senses danger.

[0198] According to one aspect of this application, based on a risk factor Bayesian network, and considering the hydraulic connections of the reservoir-lake system, the input-output structure of the scheduling model, and system risk indicators, the network structure is directionally extended starting from risk factor variables until the overall flood control system model is completed. Specifically, it can also be as follows:

[0199] Based on the flood control objectives of the reservoir-lake system, system risk index variables are defined. The highest water level Z during the lake's operation period is used as the indicator. max As a systemic risk indicator:

[0200] ;

[0201] In the formula, Z h (t) represents the lake's water level at time t; T represents the number of time periods in the scheduling period.

[0202] Define risk level classification criteria: When Z max Exceeding the warning water level Z w When Z is in a state of alert, the system is in an alert state; when Z is in a state of alert, the system is in an alert state. max Exceeding the guaranteed water level Z g At that time, the system was in a risky state.

[0203] Based on the scheduling process of the reservoir-lake system, the following intermediate process variables are defined:

[0204] (1) The maximum discharge q of reservoir i during the scheduling period max_i : ;

[0205] (2) The highest water level Z during the scheduling period of reservoir i max_i : ;

[0206] (3) Total flood volume W flowing into the lake h : ;

[0207] In the formula, q i Z(t) represents the outflow from reservoir i at time t. i Q(t) represents the water level of the i-th reservoir / lake at time t. h (t) represents the inflow to the lake at time t, which includes the flow from the reservoir discharge to the inflow section and the inflow from the lake area.

[0208] Based on the risk factor Bayesian network skeleton, directed edges are added from risk factors to intermediate process variables and from intermediate process variables to risk indicator variables according to the input-output relationship and hydraulic connection of the scheduling model.

[0209] Edges are added according to the following rules:

[0210] Rule 1 (Scheduling Input Edges): If risk factor X i Y is the input variable of the scheduling model and the intermediate variable. j X is the output variable of the scheduling model, and X i Directly affects Y j The calculation then involves adding a directed edge X. i ->Y j .

[0211] Rule 2 (Hydraulic Transmission Edge): If the intermediate variable Y i With Y j There is a hydraulic transfer relationship between them. For example, if the outflow from the reservoir evolves to the inflow section, then a directed edge Y is added. i ->Y j .

[0212] Rule 3 (Risk Convergence Edge): If the intermediate variable Y i If the value of the risk indicator R is directly affected, then a directed edge Y is added. i ->R.

[0213] For the extended Bayesian network G' BN =(V',E') performs a directed acyclicity verification, where V' contains risk factors, intermediate variables, and risk indicators, and E' contains skeleton edges and newly added edges.

[0214] If loops exist, they should be corrected according to the method used in constructing Bayesian network structures to ensure that the network satisfies the DAG constraints.

[0215] Output the complete Bayesian network G' from the risk factor variable through intermediate process variables to the risk indicator variable. BN The network structure reflects the propagation path of flood risk in the reservoir-lake system.

[0216] According to another aspect of this application, an adaptive analysis system for flood control risk of a reservoir-lake system is provided, comprising:

[0217] processor;

[0218] Memory, used to store computer programs;

[0219] The method described in any one of the procedures of this invention is implemented when the processor executes a computer program.

[0220] The reservoir-lake system flood control risk adaptive analysis system is specifically a computer system, including a central processing unit (CPU) that can perform various appropriate actions and processes based on programs stored in read-only memory (ROM) or loaded from the storage section into random access memory (RAM). For example, it can perform operations such as constructing physical causal constraint matrices, calculating risk-sensitive weights, estimating VineCopula parameters, scheduling and solving using the sparrow search algorithm, and performing Bayesian network inference. The RAM also stores various programs and data required for system operation, such as historical hydrological sequences, physical causal constraint matrices, and pre-trained Bayesian network parameters. The CPU, ROM, and RAM are interconnected via a bus.

[0221] Input / output (I / O) interfaces are connected to the bus. Multiple components in a computer system are connected to the I / O interface, including: input components such as keyboards, mice, and touchscreens, used for users to input physical causal constraint rules or set risk thresholds; output components such as cathode ray tube (CRT) displays, liquid crystal displays (LCDs), and speakers, used to display risk assessment results, scheduling plan charts, or issue flood warnings; storage components such as hard drives and solid-state drives (SSDs), used for persistent storage of hydrological data and model parameters; and communication components such as local area network (LAN) cards and modems, which communicate with external devices (such as hydrological and meteorological telemetry stations and reservoir scheduling center servers) via networks (such as the Internet and LANs).

[0222] This invention defines a physical causal constraint matrix and a causal compatible edge set, which forcibly eliminates statistical noise connections that do not conform to hydraulic connections during the structure learning stage; and uses causal scoring rules to determine the edge directions of the Bayesian network, so that the final constructed model not only fits the statistical characteristics of the data, but also strictly follows the physical laws from upstream to downstream and from input to output, eliminating the uncertainty of the model structure, improving the interpretability of the model, and solving the problem of missing physical causality.

[0223] The present invention establishes a risk-adaptive iterative mechanism. By calculating the posterior probability (risk-sensitive weight) of high-risk events using a Bayesian network, it dynamically distorts the Copula sampling density, forcing the sampling process to tend towards generating unfavorable flood control combinations of system risk factors. The optimized scheduling results are then fed back to the network for parameter updates. This closed-loop mechanism solves the problem of sparse samples in long-tail regions in traditional static sampling, enabling accurate capture of extreme flood combinations. It can effectively reveal the true risk level of the system under complex scheduling and solve the problem of underestimating extreme risks.

[0224] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all fall within the protection scope of the present invention.

Claims

1. A method for adaptive analysis of flood risk for a reservoir-lake system, characterized in that, The application relates to a method for comprehensive flood risk assessment of a system, comprising the following steps: obtaining a hydrological sequence of a river basin and a preconfigured physical causal constraint matrix, and constructing a VineCopula joint distribution model and an initial risk factor Bayesian network meeting the physical causal topological constraint based on the hydrological sequence and the preconfigured physical causal constraint matrix; calculating the posterior probability of each sample point causing the system risk index to exceed a preset risk threshold based on the initial risk factor Bayesian network, and obtaining a risk sensitivity weight function; coupling the risk sensitivity weight function to the VineCopula joint distribution model to construct a risk weighted sampling density, and sampling to generate a risk adaptive scenario set based on the risk weighted sampling density; solving an optimal scheduling scheme of a reservoir group for the risk adaptive scenario set, and feeding back the optimal scheduling scheme to the risk factor Bayesian network to update network parameters, and cyclically executing the steps of calculating the risk sensitivity weight function and updating the network parameters until the risk assessment result converges; determining a system comprehensive flood risk assessment result based on the converged risk factor Bayesian network.

2. The method of claim 1, wherein, The method for constructing the VineCopula joint distribution model and the initial risk factor Bayesian network meeting the physical causal topological constraint comprises the following steps: determining a plurality of risk factors based on a hydrological sequence of a catchment area, defining elements A of a physical causal constraint matrix ij = 1 if the risk factor i is a direct cause of the risk factor j, A ij = -1 if the risk factor j is a direct cause of the risk factor i, otherwise A ij = 0 ij ; calculating Kendall rank correlation coefficients between two risk factors based on the hydrological sequence of the river basin, and constructing a Kendall correlation coefficient matrix; The causal compatible edge set is constructed. For any two risk factors i and j, if A ij is not 0, or A ij is 0 and the absolute value of the Kendall rank correlation coefficient exceeds a preset correlation threshold, an edge connecting i and j is added to the causal compatible edge set. searching for a maximum weight spanning tree in a causal compatible edge set to obtain a first layer tree structure of the VineCopula joint distribution model.

3. The method of claim 2, wherein, The method for searching for the maximum weight spanning tree in the causal compatible edge set to obtain the first layer tree structure of the VineCopula joint distribution model comprises the following steps: taking the absolute value of the Kendall rank correlation coefficient as an edge weight, and constructing a weighted graph on the causal compatible edge set; initializing a tree node set to contain any one risk factor, and initializing a tree edge set to be empty; iteratively executing the following steps until the tree node set contains all the risk factors: in candidate edges connecting the tree node set and a non-tree node set, selecting an edge with the maximum weight and belonging to the causal compatible edge set, and adding the edge and the non-tree node connected by the edge into the tree edge set and the tree node set respectively; outputting the tree edge set as the first layer tree structure.

4. The method of claim 3, wherein, The method for constructing the VineCopula joint distribution model and the initial risk factor Bayesian network meeting the physical causal topological constraint further comprises the following steps: determining a Bayesian network skeleton based on the first layer tree structure, and calculating a causal score of each risk factor, wherein the causal score is the number of times that the risk factor acts as a cause minus the number of times that the risk factor acts as a result in the physical causal constraint matrix; Direction of each edge in the Bayesian network skeleton: If A ij =1 for the edge connecting risk factors i and j, then the direction is i->j; if A ij =-1, then the direction is j->i; if A ij =0, then the risk factor with higher causal score points to the risk factor with lower causal score. performing directed acyclic detection and correction on the directed network to obtain the initial risk factor Bayesian network.

5. The method of claim 1, wherein, The method for calculating the posterior probability of each sample point causing the system risk index to exceed the preset risk threshold to obtain the risk sensitivity weight function adopts the following formula: ； where w k (x) is the risk-sensitive weight function of the kth iteration, x is the risk factor vector, R is the systematic risk indicator, R th is the preset risk threshold, BN k is the risk factor Bayesian network of the current round.

6. The method of claim 5, wherein, The method for constructing the risk weighted sampling density adopts the following formula: ； where f k (x) is the risk-weighted sampling density of the kth iteration, c(·) is the density function of the Vine Copula joint distribution model, θ k is the current Copula parameter, u i is the marginal cumulative distribution function value of the risk factor x i , d is the dimension of the risk factor, f i (x i ) is the marginal probability density function of the risk factor x i .

7. The method of claim 1, wherein, The method for feeding back the optimal scheduling scheme to the risk factor Bayesian network to update the network parameters comprises the following steps: integrating the samples in the risk adaptive scenario set and the corresponding optimal scheduling scheme results into a historical training sample set; re-estimating the conditional probability table parameters of the risk factor Bayesian network based on the updated training sample set by using an EM algorithm to obtain an updated risk factor Bayesian network; and The difference between the risk estimation value of the high-risk interval in the current iteration and the last iteration is calculated, and if the difference is less than a preset convergence threshold or the number of iterations reaches a preset upper limit, it is determined that the risk assessment result converges.

8. The method of claim 1, wherein, A VineCopula joint distribution model and an initial risk factor Bayesian network satisfying the physical causal topological constraint are constructed, and a structure learning method based on score search is adopted, and the objective function is: max G L(G) - λ1Φ(G, G vine ) - λ2Ψ(G, C phy ); wherein G is the Bayesian network structure to be searched, L(G) is a log-likelihood function, Φ(G, G vine ) is a network structure G inconsistency penalty term with the VineCopula first layer tree structure G vine , Ψ(G, C phy ) is a physical topology penalty term of the network structure G violating the physical causal constraint matrix C phy , and λ1 and λ2 are weight coefficients. wherein G is the Bayesian network structure to be searched, L(G) is a log-likelihood function, Φ(G, G vine ) is a network structure G inconsistency penalty term with the VineCopula first layer tree structure G vine , Ψ(G, C phy ) is a physical topology penalty term of the network structure G violating the physical causal constraint matrix C phy , and λ1 and λ2 are weight coefficients. In the search process, the directed acyclic graph that maximizes the objective function is selected as the initial risk factor Bayesian network.

9. The method of claim 1, wherein, Based on the converged risk factor Bayesian network, the system comprehensive flood control risk assessment result is determined, including: Intermediate process variables are defined, including at least the maximum discharge during the reservoir operation period, the highest water level during the reservoir operation period, and the total inflow volume of the lake; Based on the converged risk factor Bayesian network, the directed edges are added to expand the network structure according to the following rules: If the risk factor is the input of the operation model and directly affects a certain intermediate process variable, a directed edge from the risk factor to the intermediate process variable is added; If there is a hydraulic transmission relationship between two intermediate process variables, a directed edge from the upstream variable to the downstream variable is added; If the intermediate process variable directly affects the system risk index, a directed edge from the intermediate process variable to the system risk index is added; Based on the expanded Bayesian network, probability reasoning is performed to output the system comprehensive flood control risk assessment result.

10. A reservoir group-lake system flood risk self-adaptive analysis system, characterized in that, It includes: A processor; A memory for storing a computer program; The processor executes the computer program to implement the method of any one of claims 1 to 9.

Images