Multi-dimensional quantification method for non-uniform evolution law of flood

Abstract

The application discloses a multi-dimensional quantification method for flood non-uniform evolution law, comprising the following steps: constructing a physically constrained regular path through a multi-channel Bayesian grey weight and a water balance penalty to physically align the measured and natural flood processes; constructing a local deformation tensor field along the trajectory distribution in the reconstructed high-dimensional phase space; decoupling the mixed regulation and storage effect into a stretching tensor component representing the peak clipping effect and a rotating tensor component representing the mispeak effect by using polar decomposition technology to generate a decoupling index set; obtaining a water engineering state vector, constructing a dynamic causal diagram containing a physically constrained regularization term, and quantifying the path-level contribution of the water engineering regulation path to the non-uniform evolution. Through tensor decoupling and physically constrained causal learning, the application realizes accurate extraction and physical attribution of the non-uniform full attribute characteristics of flood, and improves the scientificity of flood regulation and control decision.

Description

TECHNICAL FIELD

[0001] The present application relates to flood forecasting and the like technology, in particular to a multi-dimensional quantification method of flood non-uniform evolution law. BACKGROUND

[0002] As one of the most frequent natural disasters in the world, the accurate quantification of the evolution law of flood is crucial for water conservancy planning and flood control and disaster reduction. Especially under the double driving of climate change and large-scale water conservancy construction, the flood sequence presents obvious non-uniformity characteristics, resulting in a huge difference in morphology and statistical law between the measured flood process and the restored natural process. The traditional frequency analysis method based on the assumption of stationarity has been difficult to meet the needs of modern hydrological risk assessment, therefore, in-depth revealing the nonlinear dynamic variation mechanism of flood under the disturbance of human activities has important scientific value for improving engineering safety and resource utilization efficiency.

[0003] In the prior art, the research on flood non-uniformity mainly focuses on the description of statistical characteristics and the measurement of sequence similarity. The commonly used methods include using phase space reconstruction technology to analyze the chaotic characteristics of flood sequence, or using dynamic time warping (DTW) algorithm to calculate the numerical distance between sequences to evaluate the variation degree. In terms of driving mechanism identification, existing researches mainly identify the mutation points through statistical methods such as Mann-Kendall test (M-K test), and conduct sensitivity analysis by combining Sobol index or Shapley value decomposition to evaluate the contribution degree of rainfall, underlying surface change and other factors to flood elements. The above methods constitute the mainstream technical framework of current non-uniformity quantification.

[0004] However, the existing system faces two major bottlenecks in dealing with flood processes disturbed by strong human activities (such as joint operation of reservoir group), which are difficulty in coupling physical mechanism and unexplainable causal attribution. Specifically, the traditional phase space analysis only measures the Euclidean distance or time sequence deviation between trajectories, ignoring the differential geometric characteristics of trajectory evolution, making it difficult to effectively distinguish the two different physical effects of reservoir regulation, i.e. peak clipping (corresponding to radial contraction / stretching deformation of phase space trajectory) and peak shifting (corresponding to tangential rotation deformation of trajectory), and difficult to realize the fine decoupling of regulation effect. In addition, the existing driving mechanism identification mainly relies on pure data-driven statistical correlation, lacks the physical constraints of water balance equation and engineering topological structure, and is prone to produce pseudo-causal relationships that violate hydrological common sense when dealing with short sample sequences or high noise data, and it is difficult to accurately quantify the path-level contribution of pre-determined regulation path (not a single variable) to non-uniform evolution. SUMMARY

[0005] The present application provides a multi-dimensional quantification method of flood non-uniform evolution law, in order to solve at least one of the above technical problems.

[0006] Technical scheme, according to one aspect of the present application, a multi-dimensional quantification method of flood non-consistent evolution law, comprising:

[0007] Obtain the measured flood hydrograph sequence and the natural flood hydrograph sequence, and determine the physical constraint regular path of the two sequences, which represents the time point matching mapping relationship between the two sequences that meets the physical constraint condition;

[0008] Based on the physical constraint regular path, a local deformation tensor field along the trajectory distribution is established in the reconstructed high-dimensional phase space, and polar decomposition is performed on it to separate out the stretching tensor component representing the peak clipping effect and the rotation tensor component representing the peak error effect, and generate a flood non-consistent decoupling index set accordingly;

[0009] Based on the water project state vector set obtained from the reservoir operation monitoring data, the water project state vector set is taken as the target variable to quantify the path-level contribution of the water project dispatching path to the flood non-consistent evolution.

[0010] Beneficial effects, through the above technical scheme, the present application realizes the accurate extraction and physical attribution of the full attribute characteristics of flood non-consistency, and improves the scientificity of flood control scheduling decision. BRIEF DESCRIPTION OF DRAWINGS

[0011] Figure 1 It is the overall flow process diagram of the multi-dimensional quantification method of flood non-consistent evolution law.

[0012] Figure 2 It is the flow process diagram of determining the physical constraint regular path.

[0013] Figure 3 It is the flow process diagram of determining the expected weight of each channel.

[0014] Figure 4 It is the flow process diagram of generating a flood non-consistent decoupling index set. DETAILED DESCRIPTION

[0015] In order to make the person in the art better understand the present application scheme, the technical scheme in the embodiments of the present application will be described clearly and completely below in combination with the drawings in the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, not all. Based on the embodiments in the present application, all other embodiments obtained by the person skilled in the art without creative labor shall belong to the scope of protection of the present application.

[0016] Embodiment 1, describe the overall flow framework of a multi-dimensional quantification method of flood non-consistent evolution law, as Figure 1The whole process from data input, physical constraint alignment, phase space tensor decoupling to final path level contribution quantification is shown, solving the problems of inaccurate identification of nonlinear dynamic variation of flood sequence under human activity disturbance and unclear physical attribution.

[0017] In step 101, the measured flood hydrograph sequence and the natural flood hydrograph sequence are obtained, and the physical constraint regular path of the two sequences is determined, which represents the time point matching mapping relationship between the two sequences that meets the physical constraint condition.

[0018] In this embodiment, the measured flood hydrograph sequence specifically refers to the cross-section flow time series data affected by reservoir regulation or other water conservancy projects, which is usually obtained through actual observation records of hydrological stations. The natural flood hydrograph sequence specifically refers to the benchmark flow time series data obtained by water balance principle or other hydrological restoration methods, which eliminates the influence of human activities. The two sequences usually have the same sampling frequency and time length, but there are obvious differences in flood peak shape, peak time and total flood volume. The physical constraint regular path is not directly obtained by the traditional dynamic time warping algorithm based on numerical distance, but a certain physical constraint condition is introduced, such as water balance equation or non-negative reservoir capacity constraint, to establish a reasonable time point matching mapping relationship in hydrology. This mapping relationship can avoid matching a certain flow point in the measured process to a time that is physically impossible in the natural process, ensuring that the subsequent shape comparison is carried out on the correct time basis. Specifically, this step constructs a cost function containing a physical penalty term, suppresses the generation of paths that violate physical laws in the dynamic programming optimization process, and finally outputs a regular path composed of a series of time point pairs.

[0019] In step 102, based on the physical constraint regular path, a local deformation tensor field distributed along the trajectory is established in the reconstructed high-dimensional phase space, and polar decomposition is performed on the local deformation tensor field to separate the stretching tensor component representing the peak reduction effect and the rotation tensor component representing the peak shift effect, and generate a set of flood inconsistency decoupling indexes according to the rotation tensor component.

[0020] Specifically, this step firstly maps the observed flood series and the natural flood series into the high-dimensional phase space according to the warping path obtained in the previous step, and reconstructs two corresponding phase space trajectories. For each pair of matching points on the trajectories, the local geometric transformation relationship is calculated, that is, the local deformation tensor field is established. This tensor field can be described mathematically by a Jacobian matrix, which characterizes the stretching, compression and rotation deformation degree of the observed trajectory relative to the natural trajectory in the local neighborhood. Perform polar decomposition on each tensor in the tensor field, and mathematically orthogonally decompose it into a symmetric positive definite matrix and an orthogonal matrix. Among them, the stretching tensor component corresponds to the symmetric positive definite matrix obtained by decomposition, which reflects the scaling of the flood process on the amplitude axis, and physically corresponds to the peak clipping or compensation effect produced by reservoir regulation; the rotation tensor component corresponds to the orthogonal matrix obtained by decomposition, which reflects the rotation of the flood process in the phase space, and physically corresponds to the lag or peak shift effect produced by regulation. Based on the eigenvalues or rotation angles of the two components, the scalarized index set, i.e., the flood non-uniformity decoupling index set, can be further calculated to realize the physical mechanism decoupling and quantitative characterization of the flood regulation effect.

[0021] In some optional embodiments, if the computing resources are limited or only preliminary statistical feature analysis is needed, traditional dynamic time warping distance, hydrograph shape variation rate and flow concentration ratio can be used as alternative indicators. Among them, the hydrograph shape variation rate can be obtained by calculating the ratio of the standard deviation of the observed sequence to the natural sequence, and the flow concentration ratio can be obtained by calculating the ratio of the concentration index of the two. Although the above alternative indicators are difficult to provide deep physical mechanism separation like tensor decoupling, they can still reflect the overall variation degree of flood non-uniformity on a macroscopic level.

[0022] Step 103, based on the reservoir operation monitoring data, obtain the water engineering state vector set, take the flood non-uniformity decoupling index set as the target variable, and quantify the path-level contribution of the water engineering regulation path to the evolution of flood non-uniformity.

[0023] In this embodiment, the water engineering state vector set includes multi-dimensional physical parameters describing the operation state of the reservoir or other water conservancy projects, such as reservoir water level, inflow, outflow, water abandonment, and dispatching rule parameters. The above state vector is synchronized with time and reflects the dynamic process of project scheduling. This step is used to establish the causal relationship between the project scheduling behavior and the flood pattern variation. Specifically, the decoupling index set generated in step 102 is taken as the target variable to be explained, and the water engineering state vector set is taken as the explanation variable. By constructing a causal analysis model, the key scheduling chain leading to the inconsistent change of the flood pattern is identified. The path-level contribution degree is different from the traditional single variable importance score, which quantifies the entire causal chain, such as the specific contribution value of the complete path that the increase of upstream inflow leads to the rise of reservoir water level and then triggers the change of flood discharge rule, finally causing the peak clipping of downstream. Through the above quantification, the key scheduling link leading to the dominant flood variation can be accurately located, providing a scientific basis for optimizing flood control scheduling.

[0024] Embodiment 2: The specific process of data preprocessing and determination of phase space reconstruction parameters is elaborated in detail, which provides basic parameters for reconstruction of phase space trajectory and solves the problem of loss of dynamic characteristics caused by improper selection of delay step and embedding dimension in phase space reconstruction process.

[0025] Step 201: Calculate the time delay parameter of the measured flood hydrograph sequence by using the mutual information method, and determine the delay step of phase space reconstruction.

[0026] Specifically, this step first obtains the measured flood hydrograph sequence q obs (t). In order to reconstruct the phase space trajectory that can maintain the topological structure of the original system, it is necessary to determine the appropriate time delay parameter τ. In this embodiment, the mutual information method is used to determine the parameter. The mutual information method can effectively measure the nonlinear correlation of time series at different times. In specific calculation, for a given delay τ, the mutual information value I(τ) between the sequence q obs (t) and the delayed sequence q obs (t+τ) is calculated. The calculation formula of I(τ) involves the joint probability distribution and the marginal probability distribution of the sequence. Exemplarily:

[0027] I(τ)=ΣP(x(t),x(t+τ))*log2[P(x(t),x(t+τ)) / (P(x(t))*P(x(t+τ)))];

[0028] Wherein, I(τ) is the mutual information value when the time delay is τ; P(x(t)) is the marginal probability distribution of the time series at time t; P(x(t+τ)) is the marginal probability distribution of the time series at time t+τ; P(x(t),x(t+τ)) is the joint probability distribution of the time series at time t and t+τ.

[0029] As the delay τ increases, I(τ) usually presents a decaying trend. The embodiment selects the delay time corresponding to the first time I(τ) reaches a local minimum as the optimal delay step. This selection criterion can ensure that the individual components of the reconstructed vector are neither excessively correlated nor completely uncorrelated, so that the reconstructed phase space trajectory can be fully unfolded and avoid trajectory overlap or redundancy.

[0030] In step 202, the embedding dimension of the measured flood hydrograph sequence is calculated by using the false nearest neighbor method to determine the spatial dimension of the phase space reconstruction.

[0031] After determining the delay step τ, the embedding dimension m, i.e., the dimension of the phase space, also needs to be determined. Too low a dimension will cause the trajectories to be entangled in the phase space, i.e., so-called false nearest neighbors; too high a dimension will introduce unnecessary noise and increase the computational load. In this step, the false nearest neighbor method is used to determine the minimum embedding dimension. Specifically, as the embedding dimension m increases, it is calculated whether the distance between points that are very close in the phase space significantly increases in the m+1-dimensional space. If the distance significantly increases, it means that the proximity of the two points in the low-dimensional space is an illusion caused by projection, i.e., a false nearest neighbor. When the proportion of false nearest neighbors falls below a preset threshold, e.g., less than 10% or 5%, the corresponding m is determined as the optimal embedding dimension. This ensures that the geometric structure of the flood dynamic system is unfolded and not distorted in the high-dimensional phase space.

[0032] In step 203, based on the delay step and the spatial dimension, the measured flood hydrograph sequence and the natural flood hydrograph sequence are mapped to the high-dimensional phase space.

[0033] According to the delay step τ determined in step 201 and the embedding dimension m determined in step 202, a state vector in the phase space can be constructed. For each time i in the time series, an m-dimensional vector X i is constructed. This vector is composed of the observed values of the original sequence at time i and the subsequent delay times, and can be expressed as:

[0034] X i = [x(i), x(i+τ),...,x(i+(m-1)τ)];

[0035] where X i is the reconstructed state vector at the i-th time in the phase space; x(i) is the flow observation value of the original one-dimensional flood hydrograph sequence at the i-th time; τ is the optimal time delay step determined by the mutual information method; and m is the optimal embedding dimension determined by the false nearest neighbor method.

[0036] In the above manner, the one-dimensional measured flood hydrograph sequence and the natural flood hydrograph sequence are respectively converted into a set of m-dimensional vectors, i.e., the measured phase space trajectory and the natural phase space trajectory are reconstructed, the entire nonlinear evolution characteristics of the original flood dynamic system are retained, and a mathematical basis is provided for subsequent calculation of a local deformation tensor field and decoupling of a physical mechanism.

[0037] In some embodiments, in order to eliminate the dimensional influence and accelerate the convergence of subsequent calculation, the measured flood hydrograph sequence and the natural flood hydrograph sequence can also be subjected to standardization processing before phase space reconstruction, for example, the Z-score standardization method is used to convert the data into a sequence with a mean value of 0 and a standard deviation of 1. In addition, for missing values or abnormal values that may exist in the data, linear interpolation or moving average method can be used for pretreatment to ensure the continuity and integrity of the sequence.

[0038] Embodiment 3, a detailed description of a specific implementation process of a dual-process collaborative physical regularization method is as follows: Figure 2 As shown in the figure, by introducing multi-channel grey correlation analysis and Bayesian probability weight, combined with hydrological physical constraints, the problem that the traditional dynamic time warping algorithm lacks physical mechanism constraints and is prone to unreasonable time distortion when processing flood sequences is solved.

[0039] Step 301, calculate the multi-channel grey correlation score of the measured flood hydrograph sequence and the natural flood hydrograph sequence at any time point pair, and the multi-channel grey correlation score at least covers the flow amplitude, process morphology and concentration degree dimensions.

[0040] In this embodiment, in order to comprehensively measure the similarity of the measured flood and the natural flood in different dimensions, the system first constructs multi-dimensional feature channels. Specifically, the flow amplitude channel mainly measures the closeness of the flow values at two time points, which can be calculated by using the reciprocal or exponential function of the normalized Euclidean distance; the process morphology channel mainly measures the consistency of the local fluctuation trend, which can be obtained by calculating the difference of the first-order difference or local slope; and the concentration degree dimension channel focuses on the energy distribution of the flood peak area, which can be quantified by calculating the proportion of flow integration in the local time window. For the i-th time point in the measured sequence and the j-th time point in the natural sequence, the grey correlation coefficients of the above K channels are calculated respectively. In the calculation process, grey absolute correlation degree and grey relative correlation degree algorithms can be used respectively, the former focuses on the geometric shape similarity of the values, and the latter focuses on the similarity of the change rate. The grey correlation score s i,j k of each channel k is a value between 0 and 1, and the closer the value is to 1, the more similar the channel characteristics are at the matching point pair.

[0041] Step 302, constructing a position-dependent Dirichlet Bayesian weight field based on the multi-channel grey correlation score, and determining the expected weight of each channel according to the Dirichlet Bayesian weight field.

[0042] This step is used to solve the problem of fixed multi-index weight in traditional methods, and the importance of each feature channel at different matching positions is adaptively adjusted through Bayesian inference. A preset concentration coefficient is introduced to map the multi-channel grey correlation score to a parameter vector of Dirichlet distribution; a Dirichlet distribution describing the competitive relationship of multi-channel weights is constructed based on the parameter vector; the mathematical expectation of the Dirichlet distribution is calculated to obtain the expected weight for weighting different channel mismatch functions, as shown in the following formula. Figure 3

[0043] Specifically, a concentration coefficient η is introduced to map the grey correlation score obtained in the previous step to a parameter vector of Dirichlet distribution, that is:

[0044] α i,j k =η*s i,j k ;

[0045] wherein, α i,j k is the component of the Dirichlet distribution parameter vector corresponding to the kth feature channel at the matching position of the ith point of the measured sequence and the jth point of the natural sequence; η is a concentration coefficient for controlling the concentration degree of the weight distribution, which is usually a positive real number; s i,j k is the grey correlation score calculated by the kth feature channel at the matching position (i,j), and the value range is (0,1].

[0046] The Dirichlet distribution, as the conjugate prior of the multinomial distribution, can naturally describe the probability distribution of the weight vector. For each matching point pair (i,j), its multi-channel weight vector w i,j obeys the Dirichlet distribution with parameter α i,j . In order to use this probability weight in the deterministic algorithm, the mathematical expectation of this distribution is calculated as the final weight value. The specific calculation formula is:

[0047] w i,j k =α i,j k / (Σ m=1 K α i,j m );

[0048] wherein, w​i,j k is the expected weight of the kth feature channel at the matching position (i,j); K is the total number of feature channels; Σ m=1 K α i,j m is the sum of Dirichlet parameters of all channels at this position.

[0049] In the above manner, if a certain point has a much higher grey correlation score on the shape channel than on the amplitude channel, the system will automatically give the shape channel a higher weight, realizing position-dependent adaptive weighting.

[0050] Step 303, constructing a local regularization cost function based on the expected weight and a preset physical constraint penalty term, the physical constraint penalty term being used to suppress time point matching that violates hydrological physical laws.

[0051] In this embodiment, the local regularization cost function d(i,j) is composed of two parts: one part is the multi-channel mismatch distance weighted based on the Bayesian grey weight, and the other part is the physical constraint penalty term. The specific calculation formula can be expressed as:

[0052] d(i,j)=Σ(w i,j k *g k (x i ,y j ))+μ*Ψ(i,j);

[0053] where g k (x i ,y j ) is the difference measure function (such as Euclidean distance or shape difference) of the measured value x i and the natural value y j of the kth channel, μ is the weight coefficient of the physical constraint penalty term, which is usually taken as a large positive value to form a hard constraint, and a typical value is 10 6 .

[0054] The physical constraint penalty term Ψ(i,j) is constructed based on the water balance equation. When the matching relationship of the time point pair leads to a calculated virtual storage change amount that violates the non-negative storage constraint or the water balance principle, the physical constraint penalty term takes a positive value to increase the local regularization cost, otherwise the physical constraint penalty term is zero. Exemplarily:

[0055] Ψ(i,j)=max(0,-ΔV i,j );

[0056] where Ψ(i,j) is the physical constraint penalty value at the matching position (i,j); ΔV i,j is the virtual storage change amount, and the calculation formula is:

[0057] ΔV i,j =(Σ t=1 i q in (t)-Σ t=1 j q out (t))+V init -V min ;

[0058] Where, q in (t) represents the natural inflow sequence (approximate inflow); q out (t) represents the measured outflow sequence (approximate discharge); V init V is the initial storage capacity. min This refers to the dead storage capacity or the minimum allowable water storage capacity.

[0059] Specifically, for any matching hypothesis—that is, assuming the i-th moment of the measured sequence corresponds to the j-th moment of the natural sequence—the system will calculate the virtual water accumulation process under this matching relationship. If the matching causes the calculated virtual reservoir capacity change to violate the non-negative water storage constraint, such as the physically impossible situation where the accumulated outflow is significantly greater than the accumulated inflow and exceeds the initial reservoir capacity, or the time lag law of flood peak propagation is violated, such as the upstream flood peak appearing later than the downstream flood peak, then Ψ(i,j) takes a very large positive value (e.g., infinity); conversely, if the physical law is satisfied, then Ψ(i,j) is zero. This mechanism is equivalent to setting up a physical forbidden zone in the search space of dynamic programming.

[0060] Step 304: Use dynamic programming algorithm to find the globally optimal path that minimizes the cumulative regularization cost, and obtain the physically constrained regularization path.

[0061] Based on the constructed local regularization cost matrix, this step uses a dynamic programming algorithm to recursively calculate the minimum cumulative cost. Let D(i,j) represent the minimum cumulative cost from the starting point (1,1) to point (i,j), and the recursive formula is:

[0062] D(i,j)=d(i,j)+min(D(i-1,j),D(i,j-1),D(i-1,j-1)).

[0063] By backtracking to the path that minimizes D(N,M) (where N and M are the lengths of the two input sequences, respectively), a physically constrained regularized path can be obtained. This path achieves optimal sequence alignment in numerical statistics, eliminating all possible matches that violate hydrophysical common sense, and providing a high-confidence physical benchmark for subsequent tensor decoupling analysis.

[0064] Embodiment 4: detailed description of the specific implementation process of decoupling trajectory deformation in phase space, by introducing the concept of deformation tensor in continuum mechanics, combined with singular value decomposition technology, the physical mechanism of flood regulation effect is decoupled, and specific numerical calculation cases are provided to prove the feasibility of the scheme.

[0065] Step 401: for each pair of matching points on the physically constrained regularized path, calculate the Jacobian matrix of the measured flood mapping sequence relative to the natural flood mapping sequence, which describes the local differential mapping relationship of the trajectory in phase space, wherein the measured flood mapping sequence and the natural flood mapping sequence are the mapping of the measured flood hydrograph sequence and the natural flood hydrograph sequence in high-dimensional phase space.

[0066] In this embodiment, first, based on the physically constrained regularized path obtained in the previous embodiment, the correspondence between the points X obs (i) and X nat (j) on the measured phase space trajectory is determined. In order to describe the evolution from the natural state to the measured state, it is regarded as a nonlinear transformation mapping in high-dimensional space. For a point on the trajectory, the Jacobian matrix J approximately represents the local linear mapping relationship in the neighborhood of the point. Assuming that the dimension of the phase space is m, then J is an m x m matrix, and the element J uv represents the partial derivative of the u-th component of the measured vector with respect to the v-th component of the natural vector. In specific implementation, due to discrete time series data, it is difficult to directly derive, and central difference method or local neighborhood least squares fitting method can be used for approximation. For example, in the k-neighborhood of the natural trajectory point X nat (j), find the optimal matrix J that satisfies the equation X obs -X obs_center = J*(X nat -X nat_center ), where center represents the neighborhood center point. This matrix completely records the stretching, compression, rotation and shear deformation information of the flood dynamic system locally.

[0067] Step 402: approximate the Jacobian matrix by using the central difference method or the local neighborhood least squares fitting method to obtain the local deformation tensor field describing the local stretching, compression and rotation characteristics of the trajectory in each dimension.

[0068] In this step, the set of Jacobian matrices calculated for each matching point pair on the entire trajectory is defined as the local deformation tensor field.

[0069] Exemplarily, the Jacobian matrix is approximated by the central difference method, and the calculation formula is:

[0070] Let the measured flood mapping sequence be X r ={x r(1),x r (2),...,x r (N)}, the natural flood mapping sequence is X n ={x n (1),x n (2),...,x n (N)}, where each point x(i)∈R m Let F be an m-dimensional vector. For the i-th corresponding point pair, construct the local deformation tensor F(i):

[0071] ;

[0072] in It represents the outer product.

[0073] For example, the least-squares approximation formula for calculating the local deformation tensor (Jacobi matrix) is as follows:

[0074] ;

[0075] Among them, F i Let be the local deformation tensor (Jacobi matrix) at the i-th corresponding point pair; arg min F This indicates the search for the matrix F that minimizes the objective function; x r (j), x r (i) are the phase space vectors of the measured flood mapping sequence at times j and i, respectively; x n (j), x n (i) are the phase space vectors of the natural flood mapping sequence at times j and i, respectively; k is the radius of the local neighborhood window (e.g., 5).

[0076] Step 403: Perform singular value decomposition on the Jacobian matrix at each time point, decomposing the Jacobian matrix into the product of an orthogonal rotation matrix and a symmetric positive definite stretching matrix.

[0077] This step utilizes the extreme value decomposition theorem to physically separate the tensors mentioned above. Mathematically, any non-singular matrix F can be uniquely decomposed into the form F = R * U, where R is an orthogonal matrix representing a pure rotation transformation, and U is a symmetric positive definite matrix representing a pure stretching transformation. In practical calculations, singular value decomposition (SVD) is usually performed on F first, i.e., F = V * Σ * W. T Based on the results of SVD, the rotation matrix R = V * W can be constructed. T And the stretching matrix U=W*Σ*W T Continuing with the numerical example above, after decomposing the tensor F, we can obtain the rotation matrix R and the stretching matrix U. Here, matrix R describes the phase rotation angle of the flood process in phase space at that moment, while matrix U describes the stretching ratio along each principal axis.

[0078] Step 404, determining a rotation matrix as a rotation tensor component, the rotation tensor component representing a phase shift effect of the flood process on a time axis; and determining a symmetric positive definite stretch matrix as a stretch tensor component, the stretch tensor component representing a scaling effect of the flood process on a magnitude axis.

[0079] In the embodiment, the mathematical decomposition has clear hydrological physical meaning. The rotation of the phase space trajectory corresponds to the phase shift of the waveform in the time domain, i.e. the lag or advance of the occurrence time of the flood peak, i.e. the effect of peak shifting. The stretching or compression of the phase space trajectory corresponds to the increase or decrease of the waveform amplitude in the time domain, i.e. the reduction or flattening of the flood peak, i.e. the effect of peak clipping. Through the decomposition, the storage and regulation effects mixed together are separated into two orthogonal components, solving the technical problem that the traditional indicators are difficult to distinguish the peak clipping and peak shifting.

[0080] Step 405, calculating the maximum eigenvalue of the stretch tensor component, constructing a peak clipping effect index SPI based on the deviation of the maximum eigenvalue from a preset non-storage and regulation reference value; analyzing a rotation angle from the rotation tensor component, constructing a peak shifting effect index TPI based on the rotation angle; calculating the determinant of the local deformation tensor field, constructing a comprehensive deformation index CDI based on the degree of deviation of the determinant from the volume-preserving transformation; and combining the peak clipping effect index, the peak shifting effect index and the comprehensive deformation index to form a flood non-uniformity decoupling index set, as shown in the following table. Figure 4

[0081] Specifically, for the stretch matrix U, the eigenvalue sequence is calculated. The maximum eigenvalue λ1 reflects the maximum deformation ratio at the moment. The peak clipping effect index SPI can be defined as: 1 minus the ratio of λ1 to the reference value (usually 1). In the above numerical case, if the eigenvalues of U are {0.71, 0.86, 0.92, 0.96, 0.99}, the maximum stretch ratio is 0.99, or the minimum ratio in the main reduction direction is 0.71. According to the specific definition, if the reduction is concerned, SPI=1-0.71=0.29 can be taken, indicating that about 29% is reduced. For the rotation matrix R, the corresponding rotation angle θ can be calculated. The peak shifting effect index TPI can be defined as the ratio of θ to π. In the above case, if the calculated rotation angle is about 0.18 radian, then TPI reflects a phase shift of about 10.3 degrees. The comprehensive deformation index CDI is obtained by calculating the determinant det(F), representing the compression or expansion rate of the local phase space volume. Finally, the non-uniformity of the entire flood process is quantified as three independent time series indexes of SPI, TPI and CDI.

[0082] Optionally, the rotation angle analysis (for peak shifting effect) θ i =arccos((trace(R i )-1) / 2).​

[0083] where θi is the rotation angle at the i-th time point (unit: radian); arccos is the inverse cosine function; trace(Ri) is the trace of the rotation matrix Ri (i.e., the sum of the main diagonal elements); Note: this formula is applicable to the principal rotation angle approximation of 3-dimensional and above phase spaces. i i ) is the trace of the rotation matrix Ri (i.e., the sum of the main diagonal elements); Note: this formula is applicable to the principal rotation angle approximation of 3-dimensional and above phase spaces. i

[0084] SPI = (1 / N) *∑(1-λ1(i) / λ i=1 N ) for i = 1, 2, …, N. ref

[0085] where SPI is the average peak clipping effect index of the entire sequence; N is the sequence length; λ1(i) is the maximum eigenvalue of the stretching matrix Ui (representing the maximum stretching direction) at the i-th time point; λ is the reference eigenvalue under no storage regulation, usually taking a value of 1. i ref

[0086] TPI = (1 / N) *∑(θ i=1 N / π) for i = 1, 2, …, N. i

[0087] where TPI is the average peak clipping effect index of the entire sequence, taking a value in the range [0, 1], and the larger the value, the more obvious the time phase shift; θi is the rotation angle at the i-th time point. i

[0088] Comprehensive deformation index CDI = det(F) ;

[0089] CDI represents the overall deformation degree of the trajectory, and det(F) = 1 indicates a volume-preserving transformation, and deviation from 1 indicates compression or expansion.

[0090] As an alternative to the above embodiments, in scenarios where computing resources are limited or deep machine learning coupling is not required, traditional statistical-based methods can also be used to generate the index set. For example, directly calculating the Euclidean distance between the measured sequence and the natural sequence as the overall difference index; calculating the ratio of the standard deviations of the two as the process line shape variation rate FSVR, if the ratio is less than 1, it indicates that the fluctuation is weakened; calculating the ratio of the concentration indices of the two as the flow concentration ratio FCR. Although the physical meaning of the above index is not as clear as tensor decoupling, it can also reflect the macroscopic features of non-uniformity to some extent.

[0091] ​​​​​​​Further, in order to fuse features of different scales, the application can also include a feature coordination mapping step. This step constructs the SPI and TPI indexes obtained by decoupling the tensor and traditional statistical indexes such as flood peak change rate and flood volume change rate into a multi-dimensional feature vector, and performs dimension reduction fusion through manifold learning or feature projection algorithm to generate a feature set comprehensively representing the overall picture of the non-uniform evolution of the flood.

[0092] Embodiment 5: details the specific construction and attribution analysis process of the physical constraint driven mechanism dynamic graph PC-DCG, and solves the problems of false causality and lack of engineering interpretability of pure data driven methods by introducing physical prior constraints in causal structure learning.

[0093] Step 501: based on the water balance equation and the preset water conservancy engineering topology structure, a physical constraint adjacency matrix describing the physical causal relationship between the components of the water conservancy state vector is constructed. The construction of the physical constraint adjacency matrix satisfies the following conditions: when there is a direct hydrological connection or a dispatching topology connection between two state vector components, the adjacency element at the corresponding position is set to exist connection, otherwise it is set to no connection; the hydrological connection at least includes the mass conservation relationship composed of reservoir capacity change, inflow, outflow and loss term.

[0094] In this embodiment, first, the water conservancy state vector z t at each moment is defined t = [S t , Q in , Q out ,...] T , where S t represents the reservoir capacity, Q in represents the inflow, and Q out represents the outflow. In order to construct the physical constraint adjacency matrix A phys , the explicit hydrological equation and engineering topology need to be determined. For example, according to the water balance equation S t+1 = S t + Q in - Q out - Loss, it can be determined that S t , Q in and Q out are all direct physical causes of S t+1 , where Loss is the system water loss term. Therefore, in the adjacency matrix A phys , from S t to S t+1 , from Q in to S t+1 , and from Q out to S t+1The corresponding elements of the matrix are set to 1 (or a predetermined value representing a connection), while the elements between nodes that are physically impossible to have a direct causal relationship are set to 0. The matrix represents the skeleton of the physical prior knowledge of the system.

[0095] In step 502, the optimization objective function containing the physical constraint regularization term is constructed with the set of flood non-consistent decoupling indicators as the target variables, and the data-driven adjacency matrix reflecting the data-driven causal structure is learned. The optimization objective function containing the physical constraint regularization term is composed of a data fitting term and a physical constraint deviation penalty term. The physical constraint deviation penalty term adopts the L1 norm of the difference between the adjacency matrix to be learned and the physical constraint adjacency matrix, which is used to constrain the learned causal structure to approximate the physical prior structure.

[0096] This step combines the observation data to modify the causal structure. A prediction model is constructed to predict the current non-consistent decoupling indicator y t from the historical state z t , i.e. the SPI or TPI obtained in the previous embodiment. While training the model to optimize the data fitting accuracy, a physical constraint regularization term is added to the loss function. The regularization term usually adopts the form of L1 norm, i.e. λ*||A-A phys ||1, where A is the adjacency matrix to be learned, and λ is the weight coefficient. The role of the regularization term is to punish the degree of deviation of the learned matrix A from the physical prior matrix A phys .

[0097] Exemplarily, the calculation formula of the physical constraint structure learning objective function is:

[0098] L=Σ t ||y t -f(z t ;A)|| 2 +λ*||A-A phys ||1;

[0099] Where L is the total loss function; y t is the observed value of the target variable (such as the SPI indicator) at time t; f(z t ;A) is a prediction model based on the state vector z t and the adjacency matrix A; λ is the penalty coefficient (weight coefficient) of the physical constraint regularization term; A is the data-driven adjacency matrix to be learned; A phys is the physical prior adjacency matrix constructed based on the hydraulic topology; ||.||1 is the L1 norm (sum of absolute values of elements) of the matrix.

[0100] By minimizing the total objective function containing the data fitting error and the physical constraint penalty, the data-driven adjacency matrix A dataThe learned structure can explain the statistical law of the observed data and will not deviate from the basic physical common sense.

[0101] Step 503, the physical constraint adjacency matrix is weighted and fused with the data-driven adjacency matrix to generate a final dynamic causal graph structure.

[0102] In order to balance the difference between the theoretical model and the actual data, the two matrices are fused in this step. A weighted formula is used:

[0103] A final =γ*A phys +(1-γ)*A data ;

[0104] Wherein, A final is the final generated dynamic causal graph structure matrix; γ is the fusion weight coefficient of the physical prior, and the value range is [0, 1]; A data is the data-driven adjacency matrix obtained by optimizing the learning objective function through the physical constraint structure.

[0105] The fused matrix determines the directed edge structure of the final dynamic causal graph.

[0106] Step 504, in the dynamic causal graph structure, identify the directed path from the driving factor node to the target variable node, and calculate the path level contribution degree by integrating the edge level contribution degree along the directed path. The calculation of the path level contribution degree includes: Shapley value decomposition is performed on each edge in the dynamic causal graph structure to obtain the edge level contribution degree; the edge level contribution degrees of all edges constituting the directed path are accumulated or integrated to obtain the path level contribution degree representing the influence degree of the complete regulation chain on the non-consistency evolution.

[0107] After the dynamic causal graph is constructed, this step performs path level attribution analysis. All directed paths from the source node (such as changes in warehouse inflow or dispatching rules) to the target node (non-consistency index) are identified by using a graph search algorithm (such as depth-first search). For example, a typical path may be: increase in warehouse inflow -> increase in water level -> trigger flood control dispatch -> reduce discharge flow -> increase SPI index in downstream. For each path, its contribution needs to be quantified. In this embodiment, the Shapley value decomposition method is used to calculate the edge level contribution degree of each edge in the path, and then the contribution degrees of all edges in the path are accumulated or integrated to obtain the path level contribution degree of the complete path. The quantification method can answer the engineering problem of which specific dispatching chain leads to the non-consistency change of the flood.

[0108] Exemplarily, the path level contribution degree calculation formula is:

[0109] C(π')=Σ (u->v)∈π'φ(u->v);

[0110] where C(π') is the total contribution degree of path π'; π' is a directed path from the driver node to the target variable node; (u->v) represents a directed edge on path π'; φ(u->v) is the edge-level Shapley value contribution degree of edge u->v.

[0111] In some optional embodiments, as an alternative or supplement to the above attribution analysis, if it is difficult to construct an accurate physical graph, a contribution quantification method based on mutation stage adaptive weighting can also be used. This method first divides the time series into several evolution stages based on the mutation point detection results. For each stage, the contribution rate of the principal component is calculated. Further, an adaptive weight function is constructed, which includes a time decay factor and a contribution jump factor. The time decay factor gives higher weight to recent scheduling behavior, and the contribution jump factor emphasizes the importance of factors whose contribution rates change dramatically before and after the mutation. By weighting and aggregating the contribution rates of each stage, and combining the block bootstrap sampling technique to evaluate the confidence interval of the results, the robust quantification of non-consistent driving factors can also be achieved. The above method focuses on the time-varying analysis of statistical characteristics, and can be a useful supplement to the physical causal graph analysis.

[0112] Embodiment 6, provides a preferred alternative or supplement to the attribution analysis. In certain practical engineering scenarios, if accurate reservoir scheduling rules or hydraulic topology data are lacking, it is difficult to construct an accurate physical constraint adjacency matrix, or when it is necessary to focus on evaluating the non-consistent contribution of different evolution stages, this embodiment provides a principal component contribution quantification method based on mutation stage adaptive weighting. By introducing time decay and contribution jump mechanisms, and combining the block bootstrap sampling technique, the problem of ignoring stage differences and time sequence correlation in traditional methods in long sequence non-consistent quantification is solved.

[0113] Step 601, divide the principal component sequence into multiple evolution stages based on the mutation point set, and calculate the intra-stage contribution rate and inter-stage contribution difference rate for each stage.

[0114] In this embodiment, first, based on the mutation point set T' = {t1, t2,..., tK} identified in the previous step, the principal component sequence of the entire flood non-consistent index is divided into K+1 evolution stages. K

[0115] ;

[0116] where N is the total number of sequence samples (sequence length).

[0117] ​For example, if 1985 and 2003 are identified as mutation points, the sequence is divided into a benchmark stage before 1985, a transition stage from 1985 to 2003, and a recent stage after 2003.

[0118] For each stage k and each principal component p, its intra-stage contribution rate is calculated :

[0119]

[0120] where λ p (t) is the eigenvalue of the pth principal component at time t, λ j (t) is the eigenvalue of the jth principal component at time t, and P is the total number of principal components.

[0121] This index reflects the degree of explanation of the pth principal component to the overall non-consistent variation in this period, which is usually calculated by the average proportion of eigenvalues in this stage. At the same time, in order to capture the mutation characteristics in the evolution process, the inter-stage contribution difference rate ΔCR is calculated. This index quantifies the jump amplitude of principal component contribution across mutation points, and the calculation formula can be expressed as the absolute value of the difference between the contributions of the two adjacent stages divided by the maximum value of the two.

[0122]

[0123] ΔCR ranges from 0 to 1, and the larger the value of ΔCR, the more dramatic the change in the mechanism of the principal component before and after the mutation, which is a key factor driving the evolution of non-consistency.

[0124] Step 602, construct a time decay weighting function, dynamically adjust the weight according to the time distance of each stage from the current time point and the inter-stage contribution difference rate, and form a stage-adaptive weighted contribution value.

[0125] This step solves the problem of dilution of recent influences caused by using global uniform weights in traditional methods. This embodiment constructs an adaptive weight function that integrates time-dimension features of near-large and far-small and physical-dimension features of dramatic changes. Specifically, for the kth stage, its comprehensive weight w k is composed of two parts.

[0126] The first part is the time decay factor w time (k), which adopts an exponential decay form, and the formula is:

[0127]

[0128] where w time (k) is the time decay weight of the kth stage; γ' is the decay coefficient, which is determined by cross-validation, and the initial value is 0.5; tcurrent t is the current analysis time, usually the last year of the sequence; t mid (k) is the midpoint of the kth stage; T total is the time span of the entire sequence.

[0129] This factor ensures that the closer to the current time, the higher the weight of the stage, which is in line with the principle of attaching importance to recent rules in hydrological forecasting.

[0130] The second part is the contribution jump factor w jump (k), the formula is:

[0131]

[0132] where w jump (k) is the contribution jump weight of the kth stage; β is the jump enhancement coefficient, the initial value is 1.0, and the missing term of the boundary stage is 0; is the contribution difference rate from the k-1th stage to the kth stage; is the contribution difference rate from the kth stage to the k+1th stage.

[0133] This factor gives higher weight to the stage with a dramatic mutation, because the mutation often contains key driving mechanism information.

[0134] The final comprehensive weight w k is: w k =w time (k)*(w jump (k) / (∑ j=0 K w time (j)*w jump (j)))).

[0135] Based on this weight, the contribution rate of each stage is weighted and aggregated to obtain the adaptive weighted contribution value ACR p of each principal component.

[0136]

[0137] To illustrate the process more intuitively, a specific numerical calculation case is provided below: Assume that based on the data of a certain watershed from 1970 to 2020 for a total of 50 years, two mutation points in 1985 and 2003 are identified, dividing the sequence into three stages: Stage0 (1970-1984), Stage1 (1985-2002), and Stage2 (2003-2020). After calculation, the contribution rate of a certain principal component PC1 in the three stages is 0.52, 0.41, and 0.35, respectively. Calculate the difference rate between stages: the change from Stage0 to Stage1 is |0.52-0.41| / 0.52=0.21; the change from Stage1 to Stage2 is |0.41-0.35| / 0.41=0.15. Set the current time as 2020, the decay coefficient γ' as 0.5, and the jump enhancement coefficient β as 1.0.

[0138] Calculate the weight of each stage: for Stage2 (the recent stage), the time distance is the smallest, the decay factor is the largest (close to 1), and the contribution change is stable, with a moderate jump factor. For Stage0 (the early stage), the time distance is the largest, and the decay factor is the smallest. After normalization calculation, the final weight distribution may be: w0=0.18, w1=0.32, w2=0.50. The final weighted contribution value of PC1 is: ACR1=0.18*0.52+0.32*0.41+0.50*0.35=0.40. This result is more accurate than the simple arithmetic mean value (0.43) in reflecting the actual driving strength of the principal component under the current inconsistent state.

[0139] Step 603, construct a time-related block bootstrap sampling strategy based on the weighted contribution value, evaluate the uncertainty interval of the contribution value matrix through block bootstrap method, and generate a final contribution value matrix with time robustness.

[0140] Considering the obvious time correlation of flood sequence, traditional independent and identically distributed sampling will destroy this structure, leading to underestimation of uncertainty. This step uses block bootstrap sampling method. According to the autocorrelation coefficient of the sequence, the optimal block length b is determined. For example:

[0141]

[0142] Where b is the optimal block length; ceil is the ceiling function; N is the total number of sequence samples (sequence length); ρ1 is the lag-1 autocorrelation coefficient.

[0143] The original sequence is divided into several consecutive blocks with length b. Instead of sampling a single time point, the consecutive data blocks are randomly sampled and spliced into a new sample sequence when resampling. Repeat this process B times (e.g. 1000 times), and repeat the calculation of steps 601 and 602 for each bootstrap sample to obtain B sets of weighted contribution values.

[0144] Calculate the statistical characteristics of the B sets of values, and output the final contribution value matrix containing the point estimate (such as the median) and the uncertainty interval (such as the 2.5% and 97.5% quantiles).

[0145] The above method not only quantifies who drives it, but also gives a reliable range of driving strength, enhancing the scientificity of the results.

[0146] Embodiment 7, which shows the preprocessing link for mutation point identification, describes a Bayesian network mutation point identification method combined with adversarial training and attention mechanism. This embodiment is usually used as a preprocessing step for constructing a physical constraint driven mechanism dynamic graph PC-DCG and attribution analysis, providing an accurate set of mutation points as input.

[0147] Step 701, using a principal component analysis method combined with an attention mechanism on the flood inconsistency index sequence to extract mutually independent principal components.

[0148] In this embodiment, the input data is the flood inconsistency index sequence generated by the previous embodiment, including the tensor decoupling indicators SPI and TPI, and the traditional statistical indicators. In order to eliminate the redundancy between indicators and reduce the dimension, a multi-dimensional input matrix is first constructed. Unlike traditional PCA, this embodiment introduces a multi-head attention mechanism to capture the correlation characteristics of the indicator sequence in long time series. Specifically, the initially extracted principal component sequence is taken as the query Query, the key Key and the value Value, the attention score matrix is calculated, the principal components are reconstructed by weighting, and the dimension reduction effect is optimized by minimizing the reconstruction error (such as KL divergence). Finally, by adaptively adjusting the cumulative contribution rate threshold (such as 85%), the optimized principal component sequence that retains most of the variation information and is mutually independent is selected.

[0149] Step 702, constructing a Bayesian network model on the principal component sequence, and inputting the statistics of M-K test and Pettitt test as prior probability.

[0150] To overcome the defect that a single statistical test method is sensitive to noise, a probability inference model is constructed in this step. A Bayesian network is used to model the principal component sequence, and the nodes represent the mutation states at different times. In order to introduce prior knowledge, Mann-Kendall trend test and Pettitt mutation point test are performed on the original sequence respectively to obtain test statistics such as U statistics and K statistics. The statistics are normalized and mapped to the parameters of the Beta distribution as the prior probability P(Mutation) of each node in the Bayesian network. The above fusion strategy uses the computational convenience of traditional statistical methods to provide a reasonable initial guide for Bayesian reasoning.

[0151] Step 703, performing adversarial training optimization on the network based on the attention enhanced features, generating an optimized network model and identifying a set of mutation points.

[0152] In view of the random noise and outliers widely existing in the flood sequence, the embodiment introduces an adversarial training mechanism in the training process of the Bayesian network. Specifically, a small amount of adversarial perturbation generated by gradient calculation, such as Gaussian noise or perturbation in the opposite direction of the gradient, is superimposed on the input principal component features, forcing the network to make consistent mutation judgments even when facing contaminated data. The original classification error and the classification error of the adversarial sample are included in the optimization objective function. Through Min-Max game training, the robustness of the model under complex hydrological conditions is improved. After training, the time points with posterior probability exceeding a certain threshold (such as 0.95) output by the network are identified as the set of mutation points T'. This set will be used as key time node information for stage division or path attribution.

[0153] Embodiment 8 briefly describes some auxiliary or preferred feature fusion processing methods involved in the present application, further improving the comprehensiveness of the quantization results.

[0154] In some preferred embodiments, in order to balance micro-mechanism and macro-statistical characteristics, the present application adopts a feature collaborative mapping method. Specifically, the physical indicators (SPI, TPI) obtained based on tensor decoupling are aligned with the indicators (such as flood peak modulus coefficient, period flood volume deviation rate) calculated based on traditional statistical methods. Use t-SNE (t-distributed stochastic neighbor embedding) or UMAP (uniform manifold approximation and projection) and other manifold learning algorithms to map heterogeneous indicators to a unified low-dimensional embedding space. In the mapping process, a mutual information maximization objective function is constructed to preserve the nonlinear dependency relationships originally existing between different indicators. The finally generated full-attribute feature indicators contain deep physical mechanism information and intuitive statistical order information, providing more abundant information input for subsequent mutation identification and attribution analysis.

[0155] Embodiment 9, describe the subsequent processing steps about the spatiotemporal revelation of the evolution law of non-consistency, after the identification and quantification of driving factors are completed, further evolution simulation is carried out by using dynamic Bayesian network, and the transition law of the non-consistency state of flood is revealed by combining Markov chain analysis, so that the problem that the future evolution trend cannot be predicted by simple attribution analysis is solved.

[0156] Step 901, according to the identified dominant influence factors and driving intensity results, dynamic Bayesian network evolution simulation is adopted to generate a non-consistent mutation time series of flood.

[0157] In this embodiment, first, based on the dominant influence factors (i.e. principal components or physical variables with contribution degree exceeding a set threshold, such as reservoir capacity change and inflow) identified in the previous step and the determined dynamic causal diagram structure, an evolution model of dynamic Bayesian network DBN is constructed. The model discretizes the time axis into several time slices (t = 1, 2,..., T), and the node relationship in each time slice is described by a static Bayesian network, and the state transition between time slices is described by a transition probability distribution P(X t |X t-1 ). The expectation maximization EM algorithm is used to learn the model parameters. The driving factor data under different scenarios are input, for example, the rainfall is increased by 10% or the reservoir dispatching rule is adjusted in the next 10 years, and the inference algorithm (such as particle filtering) is used to perform forward recursion simulation between time slices. The simulation output is a generated non-consistent mutation time series of flood, which shows the dynamic evolution trajectory of the non-consistency intensity of flood over time under the action of certain driving mechanism.

[0158] Step 902, the transition probability matrix is calculated by Markov chain path dependence analysis to reveal the spatiotemporal law of non-consistency evolution.

[0159] Based on the simulation time series generated in step 901, this step discretizes it into several typical non-consistency state levels, such as low non-consistency, medium non-consistency and high non-consistency. The transition characteristics between states are analyzed by using Markov chain theory. Specifically, the frequency n ij of transition from state i to state j is counted, the transition probability p ij is calculated n ij / Σ j n ij , and the state transition probability matrix P is constructed. The matrix directly reflects the evolution law of the system. For example, if the matrix shows that the probability of transition from low state to high state increases significantly, it reveals that under the current climate and dispatching conditions, the flood system has a trend of evolving to high non-consistency (i.e. more difficult to predict, more variable). The path dependence-based spatiotemporal law analysis provides forward-looking decision support for formulating long-term flood control and disaster reduction plans and revising reservoir dispatching rules.

[0160] According to one aspect of the present application, in the phase space trajectory deformation tensor and regulation effect decoupling method, a specific numerical case is shown as follows:

[0161] A typical flood event in a certain basin in 2010 is selected. The measured flood hydrograph is regulated by a large reservoir, and the natural flood hydrograph is restored by water balance. Both hydrographs contain 168 hours of points. The parameters τ=3 and m=5 are used to map to a 5-dimensional phase space.

[0162] The deformation tensor is calculated for 168 corresponding points. The calculation result of the 84th point (near the flood peak) is:

[0163] The deformation tensor F(84) is:

[0164]

[0165] The polar decomposition result is:

[0166] The rotation angle θ(84) of the rotation matrix R(84) is 0.18 rad ≈ 10.3°;

[0167] The eigenvalues λ of the stretching matrix U(84) are {0.71, 0.86, 0.92, 0.96, 0.99}.

[0168] Based on the full sequence of 168 points, the three-dimensional decoupling index is calculated:

[0169] The peak clipping effect index SPI is 0.23 (indicating an average peak clipping of about 23%);

[0170] The time shift effect index TPI is 0.12 (indicating an average phase shift of about 12% × π ≈ 22°);

[0171] The comprehensive deformation index CDI is 0.31 (indicating a moderate overall deformation).

[0172] According to one aspect of the present application, in the dynamic quantification of the contribution value of the abrupt change stage adaptive weighting, a specific numerical case is shown as follows:

[0173] Based on 50 years of data from 1970 to 2020, the set of mutation points is identified as T'={1985, 2003}, and the sequence is divided into three stages:

[0174] Stage0:1970-1984 (15 years);

[0175] Stage1:1985-2002 (18 years);

[0176] Stage2:2003-2020 (17 years);

[0177] There are 3 principal components, and the contribution rate of each stage is:

[0178] Stage0: PC1 contribution rate is 0.52, PC2 contribution rate is 0.31, and PC3 contribution rate is 0.17;

[0179] Stage1: PC1 contribution rate is 0.41, PC2 contribution rate is 0.38, and PC3 contribution rate is 0.21;

[0180] Stage2: PC1 contribution rate is 0.35, PC2 contribution rate is 0.42, and PC3 contribution rate is 0.23.

[0181] Inter-stage contribution difference rate:

[0182] ΔCR 0->1 ={0.21,0.18,0.19};

[0183] ΔCR 1->2 ={0.15,0.10,0.09}。

[0184] Taking 2020 as the current time point, the attenuation coefficient γ'=0.5, and the jump enhancement coefficient β=1.0, the adaptive weight is calculated:

[0185] w0=0.18;

[0186] w1=0.32;

[0187] w2=0.50。

[0188] Weighted contribution value:

[0189] ACR1=0.18×0.52+0.32×0.41+0.50×0.35=0.40;

[0190] ACR2=0.18×0.31+0.32×0.38+0.50×0.42=0.39;

[0191] ACR3=0.18×0.17+0.32×0.21+0.50×0.23=0.21。

[0192] The 95% confidence interval obtained by the block bootstrap method (b=5, B=1000):

[0193] ACR1:[0.35,0.45];

[0194] ACR2:[0.34,0.44];

[0195] ACR3:[0.17,0.25]。

[0196] The application adopts phase space local deformation tensor field construction and polar decomposition technology, captures the differential mapping relationship between the measured and natural trajectories through the Jacobian matrix, and orthogally decomposes the tensor into a symmetric positive definite stretching matrix and an orthogonal rotation matrix. The amplitude scaling (peak clipping) and phase offset (peak shifting) of the flood process are completely decoupled from the mathematical bottom, the peak clipping effect index SPI and the peak shifting effect index TPI with clear physical meaning can be output, and the physical and quantitative characterization of the regulation and storage effect is realized. The problem that the existing method is difficult to distinguish the peak clipping and peak shifting physical effects is solved.

[0197] The application introduces a dynamic atlas of physical constraint driving mechanism, encodes the water balance equation and the reservoir topological structure into a physical prior matrix, and uses an L1 regular term to constrain causal structure learning, so as to eliminate false edges that violate physical common sense. On this basis, by integrating the edge level contribution along the directed path, not only the key factors are identified, but also the driving strength of specific scheduling chain to non-consistent evolution is accurately quantified, and the engineering interpretability and robustness of the analysis result are enhanced. The problem of lacking physical constraints in driving mechanism identification and being easy to produce false causality is solved.

[0198] In addition, the application also adopts double-process collaborative physical regularization, uses a water balance penalty term to suppress unreasonable matching, and ensures that the subsequent tensor analysis is established on the correct time benchmark. The problem of matching that may violate physical laws in time alignment is solved.

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

Claims

1. A multi-dimensional quantification method of flood non-uniform evolution law, characterized in that, The method comprises the following steps: obtaining a measured flood hydrograph sequence and a natural flood hydrograph sequence, and determining a physically constrained regular path of the two sequences, the physically constrained regular path representing a time point matching mapping relationship between the two sequences that satisfies a physical constraint condition; based on the physically constrained regular path, establishing a local deformation tensor field along a trajectory distribution in a reconstructed high-dimensional phase space, performing polar decomposition on the local deformation tensor field, separating out a stretching tensor component representing a peak clipping effect and a rotation tensor component representing a peak misplacement effect, and generating a set of flood non-uniformity decoupling indexes based on the stretching tensor component and the rotation tensor component; obtaining a set of water project state vectors based on reservoir operation monitoring data, and taking the set of water project state vectors as target variables to quantify a path-level contribution degree of a water project dispatching path to flood non-uniformity evolution; wherein the quantification of the path-level contribution degree of the water project dispatching path to the flood non-uniformity evolution comprises: based on a water balance equation and a preset water project topological structure, constructing a physically constrained adjacency matrix describing a physical causal relationship between components of the water project state vector; taking the set of flood non-uniformity decoupling indexes as target variables, constructing an optimization objective function including a physically constrained regularization term, and learning a data-driven adjacency matrix reflecting a data-driven causal structure; performing weighted fusion on the physically constrained adjacency matrix and the data-driven adjacency matrix to generate a final dynamic causal graph structure; in the dynamic causal graph structure, identifying a directed path from a driving factor node to a target variable node, and calculating the path-level contribution degree by integrating edge-level contribution degrees along the directed path.

2. The method of claim 1, wherein, obtaining a measured flood hydrograph sequence and a natural flood hydrograph sequence, and determining a physically constrained regular path of the two sequences, comprises: calculating multi-channel grey correlation scores of the measured flood hydrograph sequence and the natural flood hydrograph sequence at any time point, the multi-channel grey correlation scores including flow amplitude, process morphology and concentration degree dimensions; constructing a position-dependent Dirichlet Bayesian weight field based on the multi-channel grey correlation scores, and determining expected weights of each channel based on the Dirichlet Bayesian weight field; constructing a local regularization cost function based on the expected weights and a preset physical constraint penalty term, the physical constraint penalty term being used to suppress time point matching that violates hydrological physical laws; solving a global optimal path that minimizes the cumulative regularization cost by using a dynamic programming algorithm to obtain the physically constrained regular path.

3. The method of claim 2, wherein, constructing a position-dependent Dirichlet Bayesian weight field based on the multi-channel grey correlation scores, and determining expected weights of each channel based on the Dirichlet Bayesian weight field comprises: introducing a preset concentration coefficient to map the multi-channel grey correlation scores into a parameter vector of Dirichlet distribution; constructing a Dirichlet distribution describing the competitive relationship of multi-channel weights based on the parameter vector; calculating the mathematical expectation of the Dirichlet distribution to obtain the expected weights for weighting different channel mismatch functions.

4. The method of claim 2, wherein: the physical constraint penalty term is constructed based on a water balance equation; when the matching relationship of the time point pair causes a calculated virtual reservoir storage variation to violate the non-negative storage constraint or the water balance principle, the physical constraint penalty term takes a positive value to increase the local regularization cost, otherwise the physical constraint penalty term is zero.

5. The method of claim 1, wherein, establishing a local deformation tensor field along the trajectory distribution in the reconstructed high-dimensional phase space, including: for each pair of matched points on the regularized path, calculating the Jacobian matrix of the measured flood mapping sequence relative to the natural flood mapping sequence, which describes the local differential mapping relationship of the trajectory in the phase space, wherein the measured flood mapping sequence and the natural flood mapping sequence are the mappings of the measured flood hydrograph sequence and the natural flood hydrograph sequence in the high-dimensional phase space, respectively; using the central difference method or the local neighborhood least squares fitting method to approximate the Jacobian matrix, and obtaining the local deformation tensor field describing the local stretching, compression and rotation characteristics of the trajectory in each dimension.

6. The method of claim 5, wherein, performing polar decomposition on the local deformation tensor field to separate the stretching tensor component representing the peak clipping effect and the rotation tensor component representing the peak misplacement effect, including: performing singular value decomposition on the Jacobian matrix at each time point to decompose it into the product of an orthogonal rotation matrix and a symmetric positive definite stretching matrix; determining the orthogonal rotation matrix as the rotation tensor component, which represents the phase shift effect of the flood process on the time axis; determining the symmetric positive definite stretching matrix as the stretching tensor component, which represents the scaling effect of the flood process on the amplitude axis.

7. The method of claim 6, wherein, generating a set of flood non-uniformity decoupling indicators, including: calculating the maximum eigenvalue of the stretching tensor component, and constructing a peak clipping effect index based on the deviation of the maximum eigenvalue from a preset non-storage reference value; analyzing the rotation angle from the rotation tensor component, and constructing a peak misplacement effect index based on the rotation angle; calculating the determinant of the local deformation tensor field, and constructing a comprehensive deformation index based on the degree of deviation of the determinant from the volume-preserving transformation; combining the peak clipping effect index, the peak misplacement effect index and the comprehensive deformation index to form the set of flood non-uniformity decoupling indicators.

8. The method of claim 1, wherein, The construction of the physical constraint adjacency matrix satisfies the following conditions: when there is a direct hydrodynamic connection or a dispatching topological connection between two state vector components, the adjacency element at the corresponding position is set to exist connection, otherwise it is set to no connection; the hydrodynamic connection at least includes the mass conservation relationship composed of storage capacity change, inflow, outflow and loss term.

9. The method of claim 1, wherein: the optimization objective function including the physical constraint regularization term is composed of a data fitting term and a physical constraint deviation penalty term, the physical constraint deviation penalty term adopts the L1 norm of the difference between the learned adjacency matrix and the physical constraint adjacency matrix, and is used to constrain the learned causal structure to approximate the physical prior structure; the calculation of the path-level contribution degree includes: performing Shapley value decomposition on each edge in the dynamic causal graph structure to obtain an edge-level contribution degree; accumulating or integrating the edge-level contribution degrees of all edges constituting the directed path to obtain a path-level contribution degree representing the influence degree of the complete regulation chain on the non-uniformity evolution.

Images