River ecological flow real-time risk early warning method

Abstract

The invention discloses a riverway ecological flow real-time risk early warning method, which comprises the following steps: on the basis of historical dry season data, identifying a delay response relationship meeting physical consistency constraints through water recession baseline stripping and total variation simplex regularization, and generating an evolution parameter set by utilizing adaptive Dirichlet sampling; interval confluence is reversely deduced based on evolution parameters, a two-channel probabilistic runoff forecasting model is constructed, and water recession gating consistency constraints are embedded in training; generating section flow ensemble forecast coupled with double uncertainties through double-layer ensemble sampling of outer layer evolution parameters and inner layer forecast residual errors; and identifying key uncertainty sources based on variance decomposition and outputting risk early warning. According to the method, the problems of non-physical oscillation of evolution parameters in the dry season and lack of physical constraints in forecasting are solved, and accurate attribution of risk sources is realized.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of ecological flow risk assessment and early warning, and particularly relates to a real-time risk early warning method for river ecological flow. BACKGROUND

[0002] The technical difficulty of river ecological flow early warning lies in accurately depicting the downstream response process of the upper stream discharge after a long distance river evolution, and superimposing the prediction component of unknown interval confluence to evaluate the risk of the section flow meeting the standard. Especially in the dry season, the water flow is slow, the river bed boundary resistance is obvious, and the evolution process presents strong nonlinear characteristics; at the same time, the interval confluence is mainly dominated by groundwater base flow supply, and the flow signal is weak and difficult to directly monitor. Building a dynamic model that can couple the evolution lag effect and the randomness of confluence, and systematically quantifying the multi-source uncertainty such as parameter drift and prediction residual, is the technical key to realize high-precision and high-reliability risk early warning.

[0003] In the existing technology, the simplified model (such as the diffusion wave model) of Saint-Venant equation set or the Muskingum method is usually used in river evolution calculation, and the evolution parameters are mostly based on historical typical flood data or assumed as constant values according to the geometric characteristics of the river. In terms of runoff prediction, the new loss method and other lumped hydrological models or deep learning algorithms such as standard long short-term memory network (LSTM) are mostly used for deterministic point prediction. For uncertainty analysis, the existing research mainly uses GLUE (generalized likelihood uncertainty estimation) or Monte Carlo method to generate prediction intervals by randomly disturbing the model input or predetermined parameters.

[0004] However, the above scheme has obvious technical limitations in the application of weak signal environment in the dry season. The traditional evolution parameter identification method based on flood data is difficult to adapt to the characteristics of the dry season, and it is difficult to effectively separate the interference of strong base flow background, resulting in the non-physical numerical oscillation or negative value of the analyzed lag response relationship (unit line), which violates the water balance constraint. The purely data-driven prediction model lacks the hard constraint of physical mechanism, and is easy to output prediction trajectories that violate the law of hydraulic attenuation in the recession period, resulting in the divergence of long-forecasting-period prediction results. In addition, the existing uncertainty quantification framework usually mixes the structural error of evolution parameters and the input residual of prediction model, which is difficult to decouple the independent contributions of evolution model parameter inaccuracy and confluence prediction noise in the topology structure, resulting in difficulty in accurately attributing the risk source. SUMMARY

[0005] The application aims to provide a real-time risk early warning method for river ecological flow, in order to solve at least one of the above problems existing in the prior art.

[0006] The technical scheme provides a real-time risk early warning method for river ecological flow, which comprises the following steps:

[0007] Obtaining historical hydrological observation data of the river channel in the dry season, including upstream discharge and downstream section observation discharge;

[0008] Based on the historical hydrological observation data, a time delay response relationship satisfying the physical consistency constraint is identified, and an evolution parameter set representing the evolution process uncertainty is generated;

[0009] Based on the evolution parameter set, interval confluence characteristics are constructed, and input into a pre-constructed probabilistic runoff prediction model to obtain prediction distribution statistics of the downstream section;

[0010] Based on the evolution parameter set and the prediction distribution statistics, double set sampling is performed to generate a section discharge set prediction coupled with evolution uncertainty and prediction uncertainty;

[0011] Based on the section discharge set prediction, risk probability is calculated and key uncertainty sources are identified, and ecological flow risk warning results are output.

[0012] Beneficial effects, the present application solves the problems of non-physical shock of evolution parameters in the dry season and lack of physical constraints in prediction, and realizes accurate attribution of risk sources. BRIEF DESCRIPTION OF DRAWINGS

[0013] Figure 1 is a river ecological flow real-time risk warning method process schematic diagram provided by the embodiment of the present application.

[0014] Figure 2 is a construction method flow chart of a river ecological flow real-time risk warning model under the coupling uncertainty of dry season evolution and runoff prediction provided by the embodiment of the present application.

[0015] Figure 3 is a dry season evolution uncertainty representation and simulation flow schematic diagram provided by the embodiment of the present application.

[0016] Figure 4 is a runoff prediction uncertainty representation and simulation flow schematic diagram provided by the embodiment of the present application. DETAILED DESCRIPTION

[0017] In order for those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be described in detail below with reference to 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 those skilled in the art without creative labor shall fall within the scope of protection of the present application.

[0018] Embodiment 1, the overall flow of a river ecological flow real-time risk warning method is described in detail, as shown inFigure 1 As shown, by constructing a complete closed-loop system including data cleaning, physical identification, probabilistic prediction and attribution decision, the dual uncertainty problems of nonlinearity of evolution parameters and difficulty in monitoring interval inflow in the ecological flow risk warning of river channel in dry season are solved.

[0019] Step 101, historical hydrological observation data of the river channel in dry season is obtained, and the historical hydrological observation data includes upstream discharge flow and downstream section observed flow.

[0020] The historical hydrological observation data is usually obtained from the basin hydrological monitoring network or the reservoir operation record. Specifically, long sequence daily flow data in the dry season (usually from October to March of the next year) is collected. The upstream discharge flow data reflects the boundary input condition, and the downstream section observed flow data reflects the comprehensive response after the river evolution and interval confluence. In order to eliminate the influence of observation noise on subsequent identification, the original observation data can be preferably preprocessed, for example, high-frequency random fluctuations are removed by using sliding average or filtering algorithm to ensure the stationarity of the data.

[0021] Step 102, based on the historical hydrological observation data, a time-lag response relationship meeting the physical consistency constraint is identified, and an evolution parameter set characterizing the uncertainty of the evolution process is generated.

[0022] Due to the slow flow speed and great influence of river bed boundary in the dry season, the evolution process has strong nonlinearity. Instead of seeking a single deterministic parameter, this embodiment identifies a parameter set meeting the physical law (such as water balance). The time-lag response relationship is specifically characterized by a time-lag kernel weight vector, which reflects the time allocation proportion of upstream flow propagation to downstream. The evolution parameter set is a sampling result of a set of possible time-lag kernel weight vectors, which is used to quantify the evolution uncertainty caused by model structure and parameter selection.

[0023] Step 103, interval confluence features are constructed based on the evolution parameter set, and the interval confluence features are input into a pre-constructed probabilistic runoff prediction model to obtain prediction distribution statistics of the downstream section.

[0024] The probabilistic runoff prediction model combines physical mechanism and data-driven. For each sample in the evolution parameter set, a corresponding interval confluence sequence can be obtained by using the principle of water balance. The probabilistic runoff prediction model is a deep learning model that can output the probability distribution characteristics of the prediction result, for example, the mean and variance (or scale parameter) of the prediction value. The model adopts a double-channel input structure, which receives both the backwater features representing physical laws and the interval confluence features representing data statistical laws, realizing the deep fusion of physical information and data information.

[0025] Step 104, performing double-layer set sampling based on the evolution parameter set and the prediction distribution statistics to generate a cross-section flow set forecast coupled with evolution uncertainty and prediction uncertainty.

[0026] This step realizes the decoupling of uncertainty through a structured double-layer sampling mechanism. The outer loop samples the evolution parameter set to capture the uncertainty of the evolution process, and the inner loop samples the prediction distribution statistics output by the probabilistic runoff forecast model to capture the uncertainty of the forecast model itself and the input data. The final generated cross-section flow set forecast contains multiple possible future flow trajectories, and the divergence degree reflects the overall risk level of the system.

[0027] Step 105, calculating risk probability and identifying key uncertainty sources based on the cross-section flow set forecast, and outputting the ecological flow risk warning result.

[0028] The system calculates the risk probability based on the proportion of samples in the cross-section flow set forecast that are below the preset ecological flow threshold. At the same time, by analyzing the variance composition of the double-layer set samples, the contribution rate of evolution uncertainty and prediction uncertainty can be quantified, and the dominant uncertainty source at the current time can be identified. The ecological flow risk warning result can specifically include risk level, warning lead time, and key uncertainty source identification, providing intuitive and interpretable scientific basis for reservoir dispatching decisions.

[0029] Embodiment 2 describes the data preprocessing process based on recession baseline stripping, and details the lag response identification preparation stage, solving how to extract the net flow signal that can truly reflect the upstream discharge response from the observed data disturbed by background base flow and interval water use.

[0030] Step 201, extracting the stable recession period in the dry season from historical hydrological observation data, calibrating the interval confluence recession curve equation, and generating the recession baseline sequence for the whole period based on the interval confluence recession curve equation.

[0031] In this embodiment, first, the quality control of the obtained original flow sequence is needed. Preferably, the Savitzky-Golay filter is used to smooth the original data, and the filter window length can be set to 5 to 7 days, and the polynomial fitting order is set to 2 or 3 orders, so as to remove measurement noise while preserving waveform characteristics. From the smoothed data, identify the period without rainfall and stable upstream discharge as the recession analysis sample. Use these samples to calibrate the interval confluence recession curve equation, which is usually expressed by an exponential decay model:

[0032] Q t =Q0*e -k't ;

[0033] where, Qt Q0 is the initial flow, k' is the recession coefficient, and e is the natural logarithm base.

[0034] Based on the calibrated parameter k', the recession baseline sequence covering the entire historical period can be recursively generated, which represents the background runoff process under the condition of no new upstream inflow and interval recharge.

[0035] Step 202, from the observed flow at the downstream section, the recession baseline sequence is stripped to obtain the to-be-explained flow sequence representing the interval net runoff response.

[0036] In order to accurately identify the influence of upstream discharge on the downstream section, the background base flow and the interval water consumption must be eliminated. This embodiment obtains the to-be-explained flow sequence through baseline stripping operation. The specific calculation formula is:

[0037] y net (t) = Q obs (t) - Q base (t) + W(t);

[0038] Where y net (t) is the to-be-explained flow sequence at time t, i.e. the net response part that needs to be explained by the evolution model; Q obs (t) is the observed flow at the downstream section; Q base (t) is the recession baseline sequence generated in step 201; W(t) is the interval water consumption flow after reduction calculation.

[0039] This step decomposes the complex observed flow into background and response parts, reducing the difficulty and interference of subsequent identification.

[0040] A specific numerical example is given below to illustrate the above recession baseline stripping process.

[0041] Suppose the observed data of a river on the tth day is as follows: the observed flow at the downstream section Q obs (t) = 85 m 3 / s, the recession baseline sequence Q base (t) = 72 m 3 / s (calculated from the recession coefficient k' = 0.05 / day of calibration), and the reduced interval water consumption flow W(t) = 3 m 3 / s. According to the formula of step 202, we have:

[0042] y net (t) = 85 - 72 + 3 = 16 m 3 / s.

[0043] The to-be-explained flow sequence value is 16 m 3 / s represents the net response flow at the downstream section, which is generated by the upstream discharge propagating downstream through the river channel after the influence of background recession flow and inter-basin flow is removed. This value will be used as the target variable for subsequent lag response identification.

[0044] Step 203, based on the to-be-interpreted flow sequence and the upstream discharge, perform lag response identification.

[0045] On this basis, use the processed pure signal (to-be-interpreted flow sequence and upstream discharge) as input to enter the next stage of evolution parameter identification process.

[0046] This step-by-step processing strategy ensures that the identification process focuses on the water flow evolution itself, rather than being masked by the base flow fluctuation.

[0047] Embodiment 3, details how to construct and solve the optimization model with multiple physical constraints to obtain physically consistent lag kernel weight vector.

[0048] Step 301, based on the preset maximum lag extension length, use the upstream discharge to construct an extended lag response matrix.

[0049] In order to capture the possible long lag effect in dry season, a larger maximum lag extension length L is set, for example, 10 days or 15 days. Based on this length, an extended lag response matrix X is constructed. The matrix is usually designed as a Toeplitz structure (diagonal constant matrix), and the column vectors correspond to the upstream daily discharge at different lags. Specifically, the matrix element X t,k represents the upstream input flow I(t-k) at time t with a lag of k time steps. This matrix structure converts the convolution operation into matrix multiplication in linear algebra, which is convenient for subsequent optimization solution.

[0050] Step 302, to minimize the error between the to-be-interpreted flow sequence and the reconstructed flow, a least squares optimization model is established, which includes a sparse smoothing regularization term and a physically consistent constraint term.

[0051] Traditional least squares method is prone to oscillation and negative solution when dealing with hydrological data with multiple collinearity, which violates the physical common sense. This embodiment establishes a regularized least squares optimization model, and the objective function J(h) is specifically as follows:

[0052] J(h)=||Xh-y net ||2 2 +λ tv *||Dh||1;

[0053] Wherein, J(h) is the target function value to be minimized; h is the lag kernel weight vector to be identified; X is the extended lag response matrix; y netis the flow series to be explained; ||...||2 2 denotes the square of the Euclidean norm, which is used to measure the data fitting error; λ tv is a regularization parameter, which is preferably in the range of 0.01 to 1.0, and controls the balance between data fitting and curve smoothing: a smaller λ tv value focuses on data fitting, which may cause local oscillation of the lag kernel curve; a larger λ tv value focuses on smoothness, which may over-smooth the peak features. In practical applications, the optimal λ tv value can be selected by cross-validation or L-curve criterion, and is preferably 0.1; D is a first-order difference matrix; ||...||1 denotes the L1 norm.

[0054] The second term is a sparse smoothing regularizer, which is a total variation (TV) regularizer. This term effectively suppresses non-physical oscillation by penalizing the gradient of the weight vector (i.e., the difference between adjacent lag weights), and forces the lag kernel curve to be smooth.

[0055] Step 303: Solve the least squares optimization model to obtain a physically consistent lag kernel weight vector.

[0056] Step 304: The sparse smoothing regularizer is a total variation regularizer, which is used to suppress non-physical oscillation of the lag kernel weight vector between adjacent lag steps.

[0057] As in step 302, the total variation regularizer not only smooths the curve, but also to some extent preserves the edge features, which is suitable for describing unit hydrograph processes with single or double peak morphology. Compared with L2 regularizer (ridge regression), TV regularizer can better avoid the phenomenon of over-smoothing the peak value for the sake of smoothing. In some alternative embodiments, ridge regression can also be used as an alternative, but its suppression effect on oscillation is usually weaker than that of TV regularizer.

[0058] Step 305: The physically consistent constraint term is a simplex constraint, which is used to ensure that each element of the lag kernel weight vector is non-negative and the sum of the elements is equal to 1.

[0059] In order to ensure the physical balance of water, when solving the above objective function J(h), a strict hard constraint, i.e., the physically consistent constraint term, must be applied. The constraint is specifically expressed as a simplex constraint:

[0060] h i >=0 and ∑(h i )=1;

[0061] where h i is the i-th element of the lag kernel weight vector h.

[0062] The non-negative constraint ensures that the water flow will not flow backward, and the constraint of 1 ensures that the water released upstream will eventually pass through the downstream section (without considering transmission losses). This constraint makes the optimization problem a constrained convex optimization problem, which can be solved by the projection gradient descent method or the interior point method.

[0063] In step 306, the time-lag kernel weight vector is determined by repeatedly solving the least squares optimization model after dividing the historical hydrological observation data into multiple sub-sample windows, and selecting and retaining the effective time-lag set based on a stability threshold.

[0064] In order to improve the robustness of the identification result, the embodiment adopts a resampling strategy. The historical hydrological observation data is divided into B sub-sample windows, and the above optimization model is independently solved on each window. The frequency of each time-lag step k being identified as a non-zero weight is counted. The calculation formula of the frequency is:

[0065] freq k = (∑ i=1 B I(h k,i > δ)) / B;

[0066] Wherein, freq k is the frequency of the time-lag step k being identified as an effective time-lag, B is the total number of sub-sample windows, I is an indicator function, which takes a value of 1 when the condition is met, otherwise 0, h k,i is the weight value of the identified time-lag step k in the i-th sub-window, and δ is the threshold for determining the weight, which is preferably set to 0.01. Based on the preset stability threshold F s , which is preferably set to 0.8, the effective time-lag set is screened. Only when freq k > F s , the time-lag step is considered to be stable and effective. The final time-lag kernel weight vector can take the average or median of all sub-window results, obtaining a set of evolution parameters that satisfy both physical constraints and statistical stability.

[0067] Embodiment 4: Details the quantification process of evolution parameter uncertainty, solves the problem of fixed prior parameters in traditional methods that are difficult to reflect the real-time identification effect of the model, realizes the adaptive characteristics that the greater the error, the more scattered the distribution, and the smaller the error, the more concentrated the distribution.

[0068] In step 401, the identification residual sequence and its residual variance of the least squares optimization model are calculated.

[0069] In this embodiment, the optimal time-lag kernel weight vector h opt, combined with the extended response matrix X, to calculate the reconstructed flow sequence. The reconstructed flow is compared with the actual flow sequence y net to obtain the identification residual sequence. The variance of the residual sequence is calculated, denoted as σ res 2 . The statistical quantity directly reflects the explanation ability of the current linear evolution model for the observation data. If the residual variance is small, it means that the linear evolution assumption is established, and the identification result is highly reliable; if the residual variance is large, it means that there is nonlinear interference or model deviation, and the uncertainty of the identification result is large.

[0070] In step 402, a posterior probability distribution satisfying the simplex constraint is constructed with the lag kernel weight vector as the distribution center.

[0071] In order to quantify the uncertainty of the parameters, the embodiment does not regard the evolution parameters as a fixed vector, but regards them as random variables subject to a certain probability distribution. Considering that the evolution coefficients must satisfy the physical constraint of non-negativity and sum of 1 (simplex constraint), a Dirichlet distribution is preferably used as the posterior distribution. The Dirichlet distribution is a multivariate probability distribution defined on a simplex, which naturally satisfies the physical constraint and avoids the destruction of statistical properties caused by secondary normalization after sampling using a Gaussian distribution. When constructing, the lag kernel weight vector h opt is set as the expectation (center) of the distribution, so as to ensure that the sampling results fluctuate around the physically consistent optimal solution.

[0072] In step 403, the concentration coefficient of the posterior probability distribution is dynamically adjusted based on the residual variance, so that the dispersion degree of the posterior probability distribution is positively correlated with the residual variance.

[0073] When the posterior probability distribution is a Dirichlet distribution, the concentration coefficient of the posterior probability distribution is dynamically adjusted based on the residual variance, specifically by using an inverse mapping function to determine the parameter vector of the Dirichlet distribution.

[0074] The inverse mapping function is: κ=κ0 / max(σ res 2 ,ε);

[0075] Wherein, κ is the concentration coefficient, κ0 is a preset reference scale constant, σ res 2 is the residual variance, and ε is a preset truncation threshold to prevent the denominator from being zero.

[0076] The parameter vector of the Dirichlet distribution is determined by the product of the lag kernel weight vector and the concentration coefficient.

[0077] The above steps define the mathematical mechanism of adaptive adjustment in detail. The parameter vector a of Dirichlet distribution can be expressed as a = k * h opt . Where the concentration coefficient k determines the degree of peak of the distribution. In order to achieve adaptation, the above inverse mapping function is designed in this embodiment. Specifically, the pre-set reference scale constant k0 can be set to a value between 10 and 100 according to experience, for example, 50; the pre-set truncation threshold e is a very small positive number, for example, 10 -6 . When the model fitting effect is good (s res 2 mall), the calculated k value becomes large, making the Dirichlet distribution highly concentrated in the center h opt , and the sampling samples are small in difference; on the contrary, when the fitting effect is poor, the k value becomes small, the distribution tends to be flat, and the sampling samples are divergent, which reflects the greater parameter uncertainty.

[0078] Step 404, sampling from the adjusted posterior probability distribution multiple times to generate an evolution parameter set.

[0079] Using the determined parameter vector a, M times of random sampling are performed from the Dirichlet distribution, for example, M = 100, to obtain a set containing M time-lag weight vectors, that is, an evolution parameter set. In addition, in the rolling update process, the posterior parameter can also be updated using the exponential forgetting mechanism, and the update formula is:

[0080] a new = p * a old + (1-p) * a obs ;

[0081] Where a new is the updated Dirichlet distribution parameter vector (i.e. new posterior parameter), a old is the Dirichlet parameter vector estimated at the last time (or the last time window), represents the accumulation of historical information, p is the forgetting factor, typically 0.2, and a obs is the instantaneous parameter estimate calculated based on the latest day data.

[0082] This mechanism ensures that the parameter set can adapt to the slow changes of the river characteristics over time.

[0083] Embodiment 5, details the specific architecture of the probabilistic runoff prediction model, solves the problems of poor physical interpretability of traditional black box models and difficulty in quantifying prediction uncertainty by designing a physical + data dual-channel input structure and a distribution output layer.

[0084] At step 501, based on the evolution parameter samples in the evolution parameter set, the interval confluence sequence is backstepped using the water balance principle, and the sliding statistical characteristics of the interval confluence sequence are calculated.

[0085] In this embodiment, for each sample in the evolution parameter set, the interval confluence sequence is backstepped using the water balance equation. In order to enhance the model's ability to capture local trends, a sliding window is applied to the interval confluence sequence, for example, a window length of 7 days, and statistical characteristics within the window are calculated, including mean, standard deviation, and coefficient of variation. The statistical characteristics constitute the enhanced input of the data-driven part.

[0086] At step 502, a probabilistic runoff prediction model with a dual-channel input structure is constructed, in which the physical channel inputs the recession base sequence extracted from historical hydrological observation data, and the data channel inputs the interval confluence sequence and its sliding statistical characteristics.

[0087] The long short-term memory network (LSTM) used in this embodiment is used as the model skeleton, but the input layer is designed specifically. Specifically, the model input layer is divided into two independent channels:

[0088] Physical channel: directly input the recession base sequence Q base and its first-order difference. This part of information represents the determined physical decay law, which constitutes the physical prior basis of the model.

[0089] Data channel: input the interval confluence sequence and its sliding statistical characteristics. This part of information represents random fluctuations affected by rainfall, evaporation, and human activities, and the neural network is responsible for mining its nonlinear law.

[0090] In some optional embodiments, in addition to LSTM, the probabilistic runoff prediction model can also use gated recurrent unit (GRU) or Transformer architecture as the basic feature extractor, as long as the above dual-channel input structure is maintained.

[0091] In a specific embodiment, the structure of the LSTM network is configured as follows: the physical channel contains 1 layer of LSTM unit, and the hidden layer dimension is 32; the data channel contains 2 layers of stacked LSTM units, and each layer has a hidden layer dimension of 64. The outputs of the two channels are concatenated and then passed through a fully connected layer with a dimension of 128 for feature fusion, and then sent to the mean head and scale head, respectively. The sequence input length is set to 30 days, corresponding to a month of historical window. During training, the Adam optimizer is used, the initial learning rate is set to 0.001, the batch size is 32, and the number of training periods is dynamically determined according to the convergence of the validation set loss, usually 100 to 200 periods. To prevent overfitting, a Dropout layer is introduced after the fully connected layer, with a dropout probability of 0.2.

[0092] Step 503, output the conditional mean and conditional scale of the cross-section flow prediction distribution in the future prediction period using the probabilistic runoff prediction model.

[0093] In order to realize probabilistic prediction, the output layer of the model does not directly output a single flow value, but outputs the parameters of the assumed probability distribution. In this embodiment, it is assumed that the prediction error follows a Gaussian distribution (normal distribution), and the output layer of the model is designed as two parallel fully connected layers (Head):

[0094] Mean Head: output the conditional mean μ t , representing the most likely flow prediction value;

[0095] Scale Head: output the conditional scale σ t , i.e. standard deviation, representing the uncertainty range of the prediction. In order to ensure that the scale parameter σ t is always positive, a Softplus activation function (an activation function that ensures the output is positive) is applied at the output end of the scale head:

[0096] ;

[0097] where σ t is the conditional scale (standard deviation) of the model output at time t, which is guaranteed to be non-negative, and z t is the original output value of the neural network scale head (Scale Head) linear layer.

[0098] In other embodiments, in view of the characteristics of the skew distribution of hydrological data, it can also be assumed that the prediction value follows a Gamma distribution, in which case the model outputs the shape parameter k t and the scale parameter θ t .

[0099] Embodiment 6: By introducing a physical constraint term controlled by a gating variable into the loss function, the problem of model performance degradation caused by the imposition of physical constraints by traditional PINN (Physical Information Neural Network) in non-physical dominant periods (such as reservoir operation, heavy rainfall recharge) is solved.

[0100] Step 601, construct a loss function containing a negative log-likelihood term and a conditional physical consistency constraint term.

[0101] In this embodiment, the total loss function L total is composed of two parts:

[0102] L total =L nll +L phy ;

[0103] where L nllThis is the negative log-likelihood term, used to optimize the probability distribution parameters. For the Gaussian distribution assumption, its formula is:

[0104] ;

[0105] Among them, L nll This is the negative log-likelihood loss term. T is the total number of time steps in the forecast period. σ t y represents the conditional scale of the model prediction at time t. t Let μ be the true value of the flow rate observation at time t. t Let t be the conditional mean of the model prediction at time t.

[0106] This factor drives the mean of the model output to approximate the observed value y. t At the same time, adjust the scale σ t To reflect the magnitude of the error.

[0107] L phy This is a conditional physical consistency constraint term, used to penalize predictions that violate physical laws under certain conditions.

[0108] Step 602: During the training iteration process, monitor in real time whether the physical dominance condition is met in the current time period.

[0109] Step 603: When the physical dominance condition is met, the physical consistency constraint term is activated to force the conditional mean to follow the physical decay law.

[0110] The activation state of the conditional physical consistency constraint is controlled by a gating variable, and the calculation logic of the gating variable is as follows:

[0111] Gated variables are constructed based on the rate of change of upstream outflow and preset interval replenishment indicators.

[0112] The gating variables are constructed based on the rate of change of upstream downstream flow and the preset interval replenishment indication, and specifically include the introduction of interval water use plan change criterion.

[0113] The activation logic of the gated variable is as follows: the gated variable is activated only when the rate of change of the upstream outflow is lower than the preset fluctuation threshold, the interval replenishment indicator shows no obvious rainfall replenishment, and the pre-acquired interval water use plan remains stable.

[0114] The above steps describe the logic of the gating mechanism in detail. To determine when to activate the physical constraint, a binary gating variable g is constructed. t g t The calculation integrates three dimensions of criteria:

[0115] Upstream boundary criterion: Calculate the rate of change of upstream downstream discharge |dQ upIf / dt| is less than the preset fluctuation threshold, such as 5%, it is considered to be boundary stable.

[0116] Inter-regional replenishment criteria: Check the interval rainfall data or soil moisture data. If it is lower than the replenishment threshold, it is considered that there is no significant replenishment.

[0117] Human interference criterion: Check the water use plan of the interval. If the change in water use in adjacent time periods is less than the sudden change threshold, such as 10%, it is considered that the water use is stable.

[0118] Only when all three conditions above are met simultaneously is it determined that the current state is dominated by pure receding water, and g is set. t =1; otherwise set g t =0. The stringent activation logic ensures that physical constraints will not incorrectly suppress the model's dynamic response capability during reservoir peak shaving or heavy rainfall.

[0119] Step 604: Minimize the loss function to optimize the model parameters.

[0120] When the upstream outflow is stable and there is no strong inter-regional replenishment, the gating variable is set to the active state; the specific form of the conditional physical consistency constraint is:

[0121] ;

[0122] Among them, L phy For the conditional physical consistency constraint term, λ phy As a penalty weight, g t For gated variables, and These are the conditional mean values ​​of the current time and the previous time, respectively, e -k' This is the physical attenuation factor determined based on the interval confluence drainage curve equation.

[0123] During training, when g t When =1, L phy The term is now active. This formula uses ReLU logic (i.e., max(0,...)) to check if the predicted mean at adjacent time points conforms to an exponential decay law. Physical decay factor e -k' The k' value can be directly adopted from the drainage coefficient calibrated in the aforementioned embodiments. If the predicted value... Greater than the theoretical attenuation value This indicates that the model's prediction violates the receding water pattern, in which case a positive loss value is generated as a penalty; otherwise, the loss is 0.

[0124] Penalty weight λ phy It is usually set to a large positive number, such as 1.0 or 10.0, to ensure the priority of physical constraints. Specifically, λ phy The selection of λ needs to balance the accuracy of data fitting and the degree to which physical constraints are satisfied: when λphy When λ is too small, the penalty effect of physical constraints is insufficient, and the model may output predictions that violate the decay law during the recession phase; when λ... phy When the value is too large, the physical constraints become too strong, potentially suppressing the model's ability to respond to changes in the actual confluence. In practical applications, it is recommended to compare different values ​​on a validation set. phy Given the root mean square error of prediction and the number of physical violations, we select a λ value that achieves Pareto optimality for both. phy Value. Based on experimental experience, for typical dry season scenarios, λ phy A value of 1.0 usually achieves a good balance.

[0125] By minimizing the total loss function, the model learns to fit the data while also incorporating the physical logic of the dry season.

[0126] Example 7: By constructing a structured two-layer sampling system, the traditional mixed uncertainty is explicitly separated into evolution parameter uncertainty and prediction residual uncertainty, and the computable attribution of key uncertainty sources is realized based on the variance decomposition principle.

[0127] Step 701: Execute the outer loop: Iterate through each evolution parameter sample in the evolution parameter set and use the probabilistic runoff forecasting model to calculate the corresponding conditional mean and conditional scale.

[0128] In this embodiment, the architecture of the two-layer set sampling is designed as a nested loop structure. The outer loop is mainly responsible for capturing the uncertainty of the evolution process. Assume that the evolution parameter set contains M samples, for example, M=50, denoted as {h1,h2,...,h...}. M For each sample h i First, the corresponding inter-regional runoff sequence is inversely derived using the water balance principle and then input into the probabilistic runoff forecasting model. The model outputs the predicted distribution parameters for the next T forecast periods under these evolution parameters: the conditional mean sequence μ. i (t) and conditional scaling sequence σ i (t), where t=1,...,T. This process generates M sets of center prediction trajectories and their corresponding uncertainty bandwidths, each representing the best prediction result under a predetermined evolutionary assumption.

[0129] Step 702, execute the inner loop: for each evolution parameter sample corresponding to the conditional mean and conditional scale, generate multiple sets of forecast residual samples that satisfy the distribution characteristics by random sampling, and superimpose them to obtain the secondary flow prediction value.

[0130] The inner loop is primarily responsible for capturing forecast uncertainties caused by data noise and model structure errors. For each pair of distribution parameters {μ} obtained from the outer loop... i (t),σi (t)}, N times random sampling, for example, N = 20. The specific sampling formula is:

[0131] q i,j (t) = μ i (t) + σ i (t) * z j ;

[0132] Where z j is the jth random number drawn from the standard normal distribution N(0, 1). For each time t, the formula generates N specific flow values q i (t) that fluctuate around the mean μ i,j (t). z j represents the normalized forecast residual sample. By superimposition, the abstract distribution parameters are converted into specific flow value samples.

[0133] Step 703, combine all the secondary flow prediction values generated by the outer loop and the inner loop to form a double-layer set prediction sample.

[0134] After the above two loops, a large set containing M*N flow trajectories is obtained, for example, 50*20=1000 trajectories. These trajectories form a double-layer set prediction sample, which can be mathematically represented as a three-dimensional matrix P with dimensions M*N*T. Each trajectory corresponds to an evolution parameter-random noise combination. The structured sample set not only provides common probability prediction intervals such as 90% confidence intervals through statistical quantiles, but also provides necessary data basis for subsequent variance decomposition.

[0135] Step 704, decompose the total variance of the double-layer set prediction sample using the variance decomposition principle.

[0136] To find the risk source, the embodiment uses the variance decomposition technique (ANOVA). According to the total variance formula (Law of Total Variance), the total variance Var(Y) can be decomposed into the variance of the conditional mean and the mean of the conditional variance. In the context of double-layer set, the two components correspond to evolution uncertainty and prediction uncertainty respectively. This decomposition directly links abstract statistical concepts to physical processes.

[0137] Step 705, calculate the variance of the conditional mean on the evolution parameter set as the evolution uncertainty contribution.

[0138] In specific calculation, first, for each prediction period t, extract M conditional means μ1(t),..., μ M (t) generated by the outer loop. Calculate the variance of this group of means, which is the evolution uncertainty contribution V routing(t). The formula is:

[0139] ;

[0140] Among them, V routing (t) represents the variance contribution at time t attributed to the uncertainty of the evolution parameters; M is the total number of samples in the outer loop (evolution parameter set); μ i (t) is the conditional mean of the prediction corresponding to the i-th evolution parameter sample; This is the average of the predicted conditional means for all M evolution parameter samples. This metric measures how much the predicted center value has drifted due to changes in the evolution parameter h. A large value indicates that the uncertainty of the evolution model has dominated the divergence in the prediction.

[0141] Step 706: Calculate the mean of the conditional scale squared over the set of evolution parameters as a contribution to the forecast uncertainty.

[0142] Extract the M conditional scales σ1(t),...,σ generated by the outer loop. M (t). Calculate the average of the squares (i.e., variances) of this set of scales, which is the contribution V to the forecast uncertainty. forecast (t). The formula is:

[0143] V forecast (t)=(1 / M)*Σ(σ i (t)) 2 ;

[0144] Among them, V forecast (t) represents the variance contribution of the uncertainty at time t attributed to the prediction model (and input noise), M is the total number of samples in the outer loop (evolution parameter set), and σ i (t) represents the prediction condition scale corresponding to the i-th evolution parameter sample.

[0145] This metric measures the average level of uncertainty of the forecast model itself, given a fixed evolution parameter. A large value indicates that no matter how accurately the evolution parameter is selected, noise in the input data or the model's own accuracy limitations can still lead to significant prediction errors.

[0146] Step 707: Compare the contributions of evolution uncertainty and forecast uncertainty, and identify the one with the larger contribution as the key source of uncertainty at the current warning time.

[0147] For each forecast period t, compare V routing (t) and V forecast The magnitude of (t). If V routing (t)>V forecastIf the ratio of the evolution contribution to the total variance is greater than 0.8 (t), it is determined that the current risk mainly comes from the uncertainty of the evolution process, and it is suggested to strengthen the monitoring or calibration of the river evolution characteristics; otherwise, it is determined that the risk comes from the insufficient prediction accuracy, and it is suggested to optimize the interval confluence prediction model or improve the quality of input data.

[0148] The calculable attribution results provide direct guidance for improving the pertinence of the early warning system. For example, in a specific numerical case, if the total variance is calculated to be 12 at a certain time, the evolution contribution is 10, and the prediction contribution is 2, the system will clearly indicate that the evolution uncertainty dominates, and the decision maker should focus on the rationality of the evolution parameters.

[0149] Embodiment 8: The post-processing and operation mechanism of the early warning system from model output to final decision are described in detail. The probability bias problem is solved through online calibration, the timeliness of the system is ensured through rolling update, and the practical value of the early warning is improved through the calculation of the early warning lead time.

[0150] Step 801: Calculate the proportion of samples below the preset ecological flow threshold in the double-layer ensemble prediction sample to obtain the uncalibrated risk probability.

[0151] After obtaining the double-layer ensemble prediction sample, for any prediction period t, the number of samples with flow values less than the preset ecological flow threshold Q eco is counted low . The uncalibrated risk probability P raw (t) is calculated as:

[0152] P raw (t) = Count low / (M*N).

[0153] For example, if 200 out of 1000 samples are below the ecological threshold, the uncalibrated risk probability is 20%. Due to the deviation of the model assumption (such as normal distribution) from the actual situation, the original probability often has the problem of overconfidence or underestimation, and cannot be directly used for issuing early warnings.

[0154] Step 802: Obtain the paired data of historical prediction probabilities and historical actual observation events within the sliding window, and construct a probability calibration mapping function.

[0155] In order to correct the above deviation, the system maintains a sliding window, such as the past 90 days, to collect historical model prediction probabilities P hist and corresponding actual observation results (1 if actual water shortage occurs, otherwise 0) of the same period. Using this set of paired data, a probability calibration mapping function f calPreferably, an isotonic regression is used as the calibration algorithm, which can ensure that the calibrated probability monotonically increases with the original probability and does not require the assumption of a specific function form. The calibration function is essentially a lookup table or a piecewise linear function that maps the probability considered by the model to the frequency actually occurring historically.

[0156] At step 803, the uncalibrated risk probability is corrected using the probability calibration mapping function to obtain a calibrated risk probability, and a warning level is determined based on the calibrated risk probability.

[0157] The P raw (t) calculated in step 801 is input into the calibration function to obtain the calibrated risk probability P cal (t)=f cal (P raw (t)). According to the preset risk classification standard, the warning level is determined. For example, the following can be set:

[0158] IV (blue warning): 0.2 <= P cal <0.4;

[0159] III (yellow warning): 0.4 <= P cal <0.6;

[0160] II (orange warning): 0.6 <= P cal <0.8;

[0161] I (red warning): P cal >=0.8.

[0162] The level directly determines the urgency of the response measures.

[0163] A sequence of calibrated risk probabilities corresponding to each future prediction period is obtained; according to the chronological order, the first prediction period time when the first calibrated risk probability exceeds the corresponding risk level threshold is searched; the time difference between the prediction period time and the current time is determined as the warning lead time, and is associated with the risk level and output.

[0164] This process defines the calculation logic of the warning lead time T lead . Assuming that the current day is day 0, the model outputs a sequence of calibrated risk probabilities for the next 7 days. If the system determines that a level II warning is triggered (threshold value is 0.6), the system will start scanning from day 1 backward to find the first time t cal that satisfies P alert (t) >=0.6. If the threshold value is exceeded for the first time on day 3, the warning lead time T lead =3 days. This indicates that the manager has a 3-day time window to perform reservoir water replenishment scheduling to avoid an ecological water shortage event.

[0165] The latest upstream discharge and downstream observed discharge are obtained by advancing in fixed time steps.

[0166] The posterior distribution parameters in the identification of hysteresis response are updated using an exponential forgetting mechanism, giving higher weight to the latest observation data.

[0167] The updated posterior distribution parameters are used to regenerate the evolution parameter set, and the risk warning at the next time is performed.

[0168] The system is not one-off, but daily rolling. When the time advances to the next day, the posterior distribution parameters α of the evolution model need to be updated online after obtaining the new measured data. The update formula uses an exponential forgetting mechanism:

[0169] α new =ρ*α old +(1-ρ)*α obs ;

[0170] Where α new is the updated Dirichlet distribution parameter vector (i.e. new posterior parameter), α old is the Dirichlet parameter vector estimated at the last time (or last time window), α obs is the instantaneous parameter estimate calculated based on the latest daily data, and ρ is the forgetting factor, for example 0.2. This mechanism makes the weight of old information decay exponentially, so that the evolution parameter set always reflects the latest hydraulic characteristics of the river channel, realizing the adaptive evolution of the model.

[0171] In actual operation, it also includes data quality control and abnormal processing mechanism:

[0172] For input data anomalies: when detecting negative values, abnormal peak values exceeding the historical extreme value range, or continuous missing data for more than three days in upstream discharge or downstream observed discharge, the system performs the following processing strategies:

[0173] (1) For single-point outliers, linear interpolation of adjacent time is used for replacement;

[0174] (2) For continuous missing data (missing days not more than 3 days), extrapolation based on the recession law is used for filling;

[0175] (3) For large-scale data anomalies or long-term missing, the system suspends the warning output and issues a data quality alarm.

[0176] For model prediction anomalies: when the conditional scale σ t output by the probabilistic runoff prediction model exceeds the conditional mean μ tWhen the prediction confidence is lower than a preset proportion threshold (e.g., 200%), it is determined that the prediction confidence is too low, and the system will mark a warning information of too large prediction uncertainty in the early warning result, prompting the decision maker to refer carefully.

[0177] Embodiment 9, provide alternative technical implementation.

[0178] Regarding the alternative solution of hysteresis response identification, in the foregoing embodiment, total variation (TV) regularization is used to constrain the smoothness of the hysteresis kernel. As an alternative, ridge regression (L2 regularization) can be used instead of the TV regularization term. At this time, the regularization term in the objective function becomes λ tv *||h||2 2 Although the ridge regression has weaker suppression ability on oscillation than the TV regularization under some abrupt signals, it is more convenient to calculate, and acceptable results can also be obtained in a smooth river channel. In addition, for simplex constraints, in addition to using the projection gradient method for solving, a quadratic programming (QP) standard solver can also be used for processing, as long as the non-negative and normalization conditions are met.

[0179] Regarding the alternative solution of the output distribution of the probabilistic prediction model, in the foregoing embodiment, it is assumed that the flow prediction error follows a Gaussian distribution (normal distribution). However, in some small and medium-sized rivers, the flow in the dry season may present a significant skewed distribution. As a preferred alternative, the output distribution can be set as a Gamma distribution. At this time, the output shape parameter k t and the scale parameter θ t of the probabilistic runoff prediction model will be adjusted. The negative log-likelihood term in the loss function will also be replaced by the likelihood function form of the Gamma distribution. This alternative solution is more suitable for hydrological data with strict non-negativity and skewness characteristics.

[0180] Regarding the alternative solution of the gating variable criterion, in the foregoing embodiment, the gating variable g t is a binary variable based on a hard threshold (0 or 1). In some more refined implementations, g t can be designed as a continuous variable between [0, 1] and softened by using a Sigmoid function or fuzzy logic. For example, g t =Sigmoid(-α'*(|dQ / dt|-threshold)), where α' is a tunable scaling parameter (also known as steepness parameter or gain) that controls the transition steepness of the Sigmoid function, dQ / dt is the rate of change, and threshold is a preset threshold. This soft gating mechanism allows the physical constraints to intervene in a gradual manner, avoiding the gradient discontinuity caused by the hard switch, and is more conducive to the stability of model training.

[0181] Embodiment 10 provides an optional scheme of a river ecological flow real-time risk early warning model construction method under the coupling uncertainty of low water evolution and runoff prediction, comprising the following steps:

[0182] S1, physical consistent evolution parameter identification and uncertainty quantification under water recession constraint: based on the historical data of the smooth section of the reservoir discharge flow in the dry season, the interval confluence recession curve equation is calibrated; the maximum lag time extended response matrix is constructed, the recession baseline stripping, total variation regularization and simplex constraint are introduced, and the physically consistent lag time kernel weight is identified; the evolution parameter set is generated through adaptive posterior simplex sampling, and the uncertainty of the dry season evolution process is quantified;

[0183] S2, embedding of distributed output physical information into prediction and double-layer set generation: based on the evolution parameter set, the interval confluence sequence is backstepped, a double-channel LSTM model of output condition mean and scale is constructed, and recession gate consistency constraint is introduced in training; the probability set prediction of the cross-section flow is generated through the double-layer set method of outer evolution parameter sampling and inner prediction residual sampling;

[0184] S3, online calibration and computable attribution rolling risk early warning: based on the double-layer set sample, the uncalibrated risk probability is calculated, and the online probability calibration is carried out through the sliding window historical paired data; the variance decomposition is carried out on the set prediction sample, and the key uncertainty source is attributed; combined with the calibrated probability and the preset threshold, the early warning level and the early warning lead time are determined, and the rolling update mechanism is established to realize continuous early warning.

[0185] According to one aspect of the present application, the physical consistent evolution parameter identification and uncertainty quantification under the water recession constraint further comprises:

[0186] The daily discharge process of the upstream reservoir in the dry season, the daily observed flow process of the downstream ecological section and the interval water use plan are selected as the historical smooth section data, the interval confluence recession curve equation is calibrated, and the recession baseline sequence is generated.

[0187] The observed flow of the section is deducted from the recession baseline and corrected for water use, to obtain the to-be-interpreted flow; the maximum lag time extended response matrix is constructed, the least square model with total variation regularization and simplex constraint is established, the lag time kernel weight vector is identified, and the effective lag time set is determined through the stable selection method;

[0188] Taking the identified lag time kernel weight as the center, the residual variance is adaptively controlled, the Dirichlet posterior distribution is constructed and sampled, and the evolution parameter set is generated; the posterior parameters are updated in an exponential forgetting manner.

[0189] The determining the effective lag time set further comprises:

[0190] The observed section flow is stripped using the recession curve sequence to obtain a to-be-interpreted flow sequence;

[0191] An extended response matrix is constructed with the maximum lag time L, and the column vectors of the extended response matrix correspond to daily upstream discharge at different lag times;

[0192] An optimization model with a total variation regularization term is established, the constraint condition is that the lag kernel weight is non-negative and the sum is 1, and a physically consistent evolution coefficient vector is obtained by solving;

[0193] Based on the solving results of multiple sub-sample windows, the selection frequency of each lag time is calculated, and the lag time exceeding the stability threshold is determined as the stable and effective lag time set.

[0194] The generating evolution parameter set and the rolling updating posterior parameter further comprises:

[0195] The lag kernel weight vector is used to construct a parameter vector of a Dirichlet distribution;

[0196] The residual variance is used to adaptively adjust the concentration coefficient in the distribution, so that the distribution is dispersed when the error is large, and the distribution is concentrated when the error is small;

[0197] Samples of the evolution parameter set satisfying the simplex constraint are generated by sampling from the Dirichlet posterior distribution;

[0198] At each rolling time, the Dirichlet distribution parameters are updated in an exponential forgetting manner to fuse new and old information.

[0199] According to one aspect of the application, the physical information embedded in the prediction and the double-layer set generation output by the distribution further comprises:

[0200] Based on each sample in the evolution parameter set, the corresponding interval concentration sequence is obtained by water balance backstepping, and quality control and feature extraction are performed;

[0201] A double-channel LSTM model is constructed, wherein the recession sequence is input into the physical channel, and the interval concentration sequence and its sliding statistical features are input into the data channel; the model output is the conditional mean and the conditional scale of the prediction period;

[0202] A recession gate consistency constraint term is introduced into the training loss function, and the predicted mean is forced to follow the recession curve decay law only in the recession dominant period;

[0203] For each outer evolution parameter sample, corresponding prediction distribution parameters are generated, and two-level set prediction samples are formed by inner residual sampling.

[0204] The double-channel LSTM model is further constructed, and comprises:

[0205] The model output layer is designed with two heads, respectively outputting the conditional mean μ(t,τ) and the conditional scale s(t,τ) of each forecast period τ, where s(t,τ)>0;

[0206] The training loss function is composed of a negative log-likelihood term and a recession consistency constraint term, where the recession consistency term controls the activation timing through a gating variable;

[0207] The gating variable is determined by the upstream depletion rate, the interval supply indication, and the water use plan mutation criterion.

[0208] According to one aspect of the present application, the online calibration and computable attributable rolling risk warning further comprises:

[0209] Based on the double-layer set sample, the sample proportion of the section flow below the ecological threshold in each forecast period is counted to obtain the uncalibrated risk probability;

[0210] The historical prediction probability and actual threshold event label in the sliding window are used to fit a monotone segmented calibration function to map the uncalibrated probability online to obtain the calibrated risk probability;

[0211] The double-layer set sample is subjected to variance decomposition to calculate the evolution uncertainty contribution, the prediction uncertainty contribution, and the interaction term contribution, and the maximum contribution term is determined as the key uncertainty source;

[0212] According to whether the calibrated risk probability exceeds the preset risk classification threshold, the risk level of each forecast period is determined, and the earliest forecast period reaching the predetermined warning level is calculated as the warning lead time;

[0213] With fixed time steps, the latest observation and scheduling information are integrated, and the physical consistency evolution parameter identification and uncertainty quantification under the recession constraint are repeatedly performed until the process of calculating the earliest forecast period reaching the predetermined warning level as the warning lead time is achieved, realizing dynamic updating of risk warning.

[0214] According to one aspect of the present application, a method for constructing a real-time risk warning model of river ecological flow under the coupling uncertainty of low-flow evolution and runoff prediction is provided, and the overall process is as shown in Figure 2 The method is trained based on historical data, and through the fusion of physical consistency evolution parameter identification, distribution output runoff prediction, double-layer set probability calculation, online calibration, and attributable rolling warning, dynamic and probabilistic monitoring and warning of ecological flow risk are realized.

[0215] Step S01, low-flow evolution uncertainty representation and simulation are performed, as shown in Figure 3 .

[0216] Data from a stable period during the dry season (usually from September to February of the following year) with no effective rainfall, no sudden changes in reservoir scheduling, and stable water use plans within the region were selected, including the daily outflow Q from the upstream reservoir. rel (t), Daily observed flow rate at downstream ecological section Q obs (t) and the interval daily water use plan (interval water use flow after restoration calculation) W(t).

[0217] Savitzky-Golay smoothing was applied to the original flow sequence to suppress short-term fluctuations. Based on the principle of water balance, multiple sets of interval confluence sequences were derived through back-calculation. Using the period of continuous stable upstream discharge as a marker, the recession segment of the interval confluence was extracted from the back-calculation results. An exponential or linear model was used to fit the recession segment, and the parameters θ of the recession curve equation were calibrated. rec And generate the full-time recession baseline sequence Q. base (t).

[0218] Furthermore, physical consistency time delay kernel identification under the constraint of water discharge is performed.

[0219] Calculate the unexplained flow y net (t):

[0220] y net (t)=Q obs (t)-Q base (t)+W(t);

[0221] Where W(t)>0 indicates that the water intake in the interval is added back after stripping to restore the natural runoff state.

[0222] Construct the extended response matrix X of the maximum delay L, whose elements are:

[0223] X(t,k)=Q rel (tk), k=0,1,…,L;

[0224] Establish an optimization model with total variation regularization and simplex constraints:

[0225] ;

[0226] Constraints:

[0227] ;

[0228] Where w is the evolution coefficient (weight vector), w k λ is the weighting coefficient corresponding to the time delay k; tv This is the total variation regularization coefficient, used to enhance the continuity of the time-delay kernel.

[0229] Identifying stable and effective time delay:

[0230] The historical data is divided into B sub-sample windows, and the above optimization is solved respectively to obtain the weight vector w (b) The selection frequency π k of each lag k is calculated

[0231] ;

[0232] A stable and effective lag set Ω is defined: Ω = {k | π k ≥ ρ};

[0233] In the formula, ε w is the weight judgment threshold, for example 0.01, is the stability threshold, for example 0.8.

[0234] Further, the adaptive posterior simplex sampling generates an evolution parameter set.

[0235] The Dirichlet posterior distribution parameter is constructed:

[0236] ;

[0237] Wherein, is the weight of physical consistent lag kernel identification, and α0 is the prior baseline parameter, which can be 1.0.

[0238] Based on the residual variance σ res 2 The concentration coefficient κ is adaptively set:

[0239] ;

[0240] Wherein, κ0 is a scale constant, and ε is a small positive number to avoid division by zero.

[0241] Sampling from the Dirichlet distribution generates an evolution parameter set:

[0242] ;

[0243] In the rolling update, the posterior parameter is updated in an exponential forgetting manner:

[0244] ;

[0245] Wherein, η is the forgetting factor, for example 0.2, is a full 1 vector.

[0246] Step S02, runoff prediction uncertainty representation and simulation are performed, as shown in Figure 4 .

[0247] For each sample w (m) in the evolution parameter set, the interval concentration sequence is backstepped through water balance.

[0248] ;

[0249] Quality control is performed on the backstepping sequence (eliminate negative values and abnormal mutations), and the mean, standard deviation and coefficient of variation in the 7-day sliding window are calculated as data channel input features.

[0250] Further, the two-channel LSTM model for distribution output is constructed and trained.

[0251] Model input:

[0252] Physical channel: Baseline sequence Q base (t) and its first-order difference.

[0253] Data channel: Interval confluence sequence Q lateral (t) and its sliding statistical features (mean, standard deviation, and coefficient of variation).

[0254] Model output:

[0255] The model sets two output heads, respectively outputting the conditional mean μ(t,τ) and the conditional scale s(t,τ) of each prediction period τ, where s(t,τ)>0 is guaranteed by the softplus activation function.

[0256] Design the loss function, the total loss function L total is composed of the negative log-likelihood term L nll and the consistency constraint term L rec of the backwater:

[0257] L total =L nll +β rec ·L rec ;

[0258] Where:

[0259] β rec is the backwater constraint weight coefficient;

[0260] ;

[0261] ;

[0262] H is the set of prediction time ranges, Q true (t+τ) is the true observed flow, f rec (μ(t,τ),θ rec ) is the theoretical backwater function, θ rec is the physical parameter of the backwater function (backwater curve equation parameter);

[0263] The gated variable g(t)∈{0,1} is determined by a combination of criteria such as the upstream discharge rate of change and the abrupt change in water use in the interval, and is activated (set to 1) only during the period when water recedes.

[0264] Furthermore, a two-layer ensemble forecast is generated.

[0265] Outer loop (evolutionary uncertainty): for each evolution parameter sample w (m) Repeat the above steps to obtain the corresponding forecast distribution parameter μ. (m) (t,τ) and s (m) (t,τ).

[0266] Inner loop (forecast uncertainty): For each outer sample m, sample N standard normal random numbers z. (n) ~N(0,1), generate secondary flow samples:

[0267] ;

[0268] The final result is a two-layer ensemble forecast sample of size M×N.

[0269] Step S03: Construct and continuously update the real-time risk early warning model for river ecological flow.

[0270] Perform online probability calibration.

[0271] Collect sliding windows Historical data (e.g., the last 90 days) ,in, Labels for actual subthreshold events. Q represents the current uncalibrated risk probability. eco This is a preset ecological flow threshold.

[0272] The monotonic calibration function C is fitted using the equal-frequency binning method or ordinal-preserving regression. τ (·)

[0273] Current uncalibrated risk probability Perform calibration:

[0274] .

[0275] Furthermore, attribution of key sources of uncertainty is conducted:

[0276] Calculate the outer conditional mean:

[0277] ;

[0278] Perform variance decomposition:

[0279] ;

[0280] where, U evo is the contribution of evolution uncertainty; U for is the contribution of prediction uncertainty; U int is the contribution of interaction term, the remaining part after the sum of the contributions of the first two terms is subtracted from the total variance.

[0281] Determine the critical source source(t,τ):

[0282] ;

[0283] On this basis, the early warning output and rolling update are carried out.

[0284] Set the risk classification threshold γ r , such as γ1=0.5, γ2=0.8, determine the risk level according to the calibrated probability:

[0285] ;

[0286] Calculate the early warning lead time of reaching the predetermined early warning level r * :

[0287] ;

[0288] Roll forward with fixed time steps (usually 1 day):

[0289] 1) Include the latest observed flow Q obs (t), the outflow Q rel (t) and the updated scheduling plan;

[0290] 2) Re-execute the online probability calibration to the early warning output and rolling update step, update the evolution parameter posterior, the prediction distribution, the calibrated probability and the attribution result;

[0291] 3) Output the updated risk level, early warning lead time and critical uncertainty source, realize continuous and adaptive ecological flow risk early warning.

[0292] The application adopts the backwater baseline stripping technology to purify the net confluence signal, and jointly applies total variation (TV) regularization and simplex constraint, effectively suppresses the oscillation in numerical calculation, forces the identified lag kernel weight to meet the non-negative, normalized and smoothness requirements, and ensures that the evolution model strictly follows the water balance principle. Solve the problem of non-physical oscillation and negative value in the evolution parameter identification under weak signal in dry season.

[0293] The application constructs a physical and data dual-channel model, and embeds a backwater gate consistency constraint in a training loss function. The mechanism automatically activates a physical penalty term in a backwater dominant period by using a gate variable, forces the model output to follow an exponential decay law, and eliminates the phenomenon of divergent or physically incorrect prediction results in a long prediction period. The problem of a pure data-driven prediction model violating the physical decay law in the backwater stage is solved.

[0294] The application establishes a double-layer set sampling architecture of evolution parameters and prediction residuals, realizes separation and independent quantification of two types of uncertainties in a topological structure by using a nested structure of outer-layer Dirichlet adaptive sampling and inner-layer random sampling, and combines a variance decomposition technique, thereby solving the problem of difficult accurate positioning of risk sources. The problem of difficult decoupling and attribution of evolution parameter errors and prediction residuals is solved.

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

Claims

1. A method for real-time risk early warning of river ecological flow, characterized in that, include: Obtain historical hydrological observation data of the river during the dry season, including upstream discharge and downstream cross-sectional observation flow; Based on historical hydrological observation data, the time-delay response relationship that satisfies the physical consistency constraint is identified, and a set of evolution parameters representing the uncertainty of the evolution process is generated. Based on the evolution parameter set, the inter-regional confluence characteristics are constructed and input into the pre-constructed probabilistic runoff forecasting model to obtain the predicted distribution statistics of the downstream section. Two-layer ensemble sampling is performed based on the evolution parameter set and the prediction distribution statistics to generate a cross-sectional flow ensemble forecast that couples evolution uncertainty and forecast uncertainty; Based on cross-sectional flow ensemble forecasts, risk probabilities are calculated and key sources of uncertainty are identified, resulting in an ecological flow risk warning.

2. The method according to claim 1, characterized in that, Based on historical hydrological observation data, the time-delay response relationships that satisfy physical consistency constraints are identified, including: Extract the stable water receding section during the dry season from historical hydrological observation data, calibrate the interval confluence water receding curve equation, and generate the water receding baseline sequence for the entire time period. By stripping the recession baseline sequence from the downstream cross-sectional flow observations, an unexplained flow sequence characterizing the inter-regional net confluence response is obtained; Based on the flow sequence to be explained and the upstream outflow flow, a time-delayed response identification is performed; Execution delay response identification includes: Based on the preset maximum delay extension length, an extended delay response matrix is ​​constructed using the upstream outflow. To minimize the error between the traffic sequence to be explained and the reconstructed traffic, a least squares optimization model containing sparse smooth regularization terms and physical consistency constraints is established. Solve the least squares optimization model to obtain a physically consistent time-delay kernel weight vector.

3. The method according to claim 2, characterized in that, The least squares optimization model satisfies the following conditions: The sparse smoothing regularization term is the total variation regularization term, used to suppress non-physical oscillations of the lagging kernel weight vector between adjacent lagging steps; The physical consistency constraint is a simplex constraint, used to ensure that each element of the time-delay kernel weight vector is non-negative and the sum of the elements is equal to 1; The time delay kernel weight vector is determined by repeatedly solving the least squares optimization model after dividing historical hydrological observation data into multiple sub-sample windows, and then filtering and retaining the effective time delay set based on the stability threshold.

4. The method according to claim 2, characterized in that, Generate a set of evolution parameters that characterize the uncertainties in the evolution process, including: Calculate the identification residual sequence and its residual variance of the least squares optimization model; Using the time-delay kernel weight vector as the distribution center, construct a posterior probability distribution that satisfies the simplex constraint; The concentration coefficient of the posterior probability distribution is dynamically adjusted based on the residual variance, so that the dispersion of the posterior probability distribution is positively correlated with the residual variance. Multiple samplings are performed from the adjusted posterior probability distribution to generate a set of evolution parameters.

5. The method according to claim 4, characterized in that, The posterior probability distribution is a Dirichlet distribution; the concentration coefficient of the posterior probability distribution is dynamically adjusted based on the residual variance, which is a parameter vector of the Dirichlet distribution determined by the inverse proportional mapping function; The inverse proportional mapping function is: κ = κ0 / max(σ) res 2 ,ε); Where κ is the concentration coefficient, κ0 is the preset baseline scale constant, and σ res 2 ε is the residual variance, and ε is the preset truncation threshold to prevent the denominator from being zero. The parameter vector of the Dirichlet distribution is determined by the product of the time-delay kernel weight vector and the concentration coefficient.

6. The method according to claim 2, characterized in that, Inter-regional runoff characteristics are constructed based on the evolution parameter set and input into a pre-constructed probabilistic runoff forecasting model, including: Based on the evolution parameter samples in the evolution parameter set, the interval confluence sequence is inferred by using the water balance principle, and the sliding statistical characteristics of the interval confluence sequence are calculated. A probabilistic runoff forecasting model with a dual-channel input structure was constructed, in which the physical channel inputs the recession baseline sequence extracted from historical hydrological observation data, and the data channel inputs the interval confluence sequence and its sliding statistical characteristics. The conditional mean and conditional scale of the predicted cross-sectional flow distribution within the future forecast period are output using a probabilistic runoff forecast model.

7. The method according to claim 6, characterized in that, The training process for a probabilistic runoff forecasting model includes: Construct a loss function that includes a negative log-likelihood term and a conditional physical consistency constraint term; During the training iteration process, it is monitored in real time whether the physical dominance condition is met in the current time period; When the physical dominance condition is met, the physical consistency constraint term is activated to force the conditional mean to follow the physical decay law. Minimize the loss function to optimize the model parameters.

8. The method according to claim 7, characterized in that, The activation state of the conditional physical consistency constraint is controlled by a gating variable, and the calculation logic of the gating variable is as follows: A gated variable is constructed based on the rate of change of upstream outflow and a preset interval replenishment indicator; when the upstream outflow is stable and there is no strong interval replenishment, the gated variable is set to the active state. The specific form of the conditional physical consistency constraint is as follows: ； Among them, L phy For the conditional physical consistency constraint term, λ phy As a penalty weight, g t For gated variables, and These are the conditional mean values ​​of the current time and the previous time, respectively, e -k' This is the physical attenuation factor determined based on the interval confluence drainage curve equation.

9. The method according to claim 6, characterized in that, Generate a ensemble forecast of cross-sectional flows that combines evolutionary and forecast uncertainties, including: Execute the outer loop: Iterate through each evolution parameter sample in the evolution parameter set and use the probabilistic runoff prediction model to calculate the corresponding conditional mean and conditional scale; Execute the inner loop: For each evolution parameter sample, the conditional mean and conditional scale are used to generate multiple sets of forecast residual samples that satisfy the distribution characteristics through random sampling, and then superimposed to obtain the secondary flow prediction value; The second-level flow predictions generated by all outer and inner loops are combined to form a two-layer ensemble forecast sample.

10. The method according to claim 9, characterized in that, Risk probabilities are calculated and key sources of uncertainty are identified based on cross-sectional flow ensemble forecasts, including: The total variance of the two-level ensemble forecast sample is decomposed using the principle of variance decomposition. The variance of the conditional mean over the set of evolution parameters is calculated as a contribution to the evolution uncertainty. The mean of the conditional scale squared over the set of evolution parameters is calculated as a contribution to the forecast uncertainty. By comparing the contributions of evolution uncertainty and forecast uncertainty, the one with the larger contribution is identified as the key source of uncertainty at the current warning time.

Images