Flood simulation real-time correction method considering system non-Gaussian property

Abstract

The invention discloses a flood simulation real-time correction method considering system non-Gaussian property. The method comprises the following steps: constructing a hydrological model for flood simulation in a basin to be simulated; constructing an ensemble Kalman filter coupled with the hydrological model; training a routing network model based on the hydrological model and an ensemble Kalman filter to process the non-Gaussian property of the hydrological model; and carrying out real-time correction on a hydrological model forecasting result based on the trained routing network model. The method has the advantages of stable and reliable data source, high calculation efficiency, objective and reasonable result and the like, and the simulation forecast precision of the hydrological model on the flood process is improved.

Description

Technical Field

[0001] The present invention relates to hydrology and deep learning technology, and in particular to a real-time correction method for flood simulation considering system non-Gaussianity. Background Art

[0002] The terrestrial hydrological cycle system is undergoing rapid and profound changes, which are reflected in the non-steady-state characteristics of river basin hydrological processes and frequent floods and waterlogging disasters.

[0003] Flood forecasting is one of the primary non-engineering flood prevention and mitigation measures. Generally, flood forecasting schemes are based on historical data. When actual conditions differ from historical data, forecast errors may occur. In such cases, it is necessary to update and correct the original forecast results in real time based on newly acquired observational information. This process of adjusting flood forecasts based on real-time information is known as real-time correction of flood forecasts. State variable correction is a widely used theory in real-time correction. This involves using a sequential assimilation method to gradually assimilate observed river flow cross-sections during the forecast process to update the state variables of the hydrological model, enabling real-time correction of flood forecasts.

[0004] Hydrological models involve the simulation of multiple complex physical processes and are complex, high-dimensional, nonlinear systems whose parameters, simulation outputs, and measurement errors exhibit non-Gaussian properties. Existing real-time correction methods based on sequential assimilation techniques are often constrained to varying degrees by the system's high-dimensional nonlinearity and non-Gaussianity. For example, the Gaussian distribution assumption of the ensemble Kalman filter and the degradation of particle weights in high-dimensional nonlinear systems caused by particle filters hinder accurate flood forecasting based on hydrological models. In contrast, deep learning models excel at extracting complex features from the system and learning non-Gaussian relationships from data. However, they struggle with the non-stationarity of the background error covariance and lack interpretability.

[0005] For example, patent CN117851746A discloses a real-time correction method for time-varying parameters of a hydrological model based on an ensemble Kalman filter. This prior art uses ensemble Kalman filtering technology to update the state of the hydrological model, thereby achieving the purpose of predicting flood processes. However, this scheme does not take into account that the non-Gaussianity in the hydrological model may violate the Gaussian assumption required by the ensemble Kalman filter, which may have an adverse impact on the accuracy of flood forecast results.

[0006] For example, Wang Shuang et al. published "Satellite remote sensing monitoring soil moisture model of Hetao irrigation area based on data call" (Northwest A&F University, 2023.DOI:10.27409 / d.cnki.gxbnu.2023.002245). This technology uses ensemble Kalman filtering technology to assimilate satellite remote sensing soil moisture into the HYDRUS-1D soil moisture model to improve model accuracy. However, this scheme also does not take into account the non-Gaussianity in the soil moisture model, which may violate the Gaussian assumption required by the ensemble Kalman filter, and thus may have an adverse effect on the accuracy of assimilation.

[0007] Furthermore, a single deep learning model and a single input feature pattern have limited learning capabilities for non-Gaussian relationships. For example, "A Basin Runoff Simulation Method Integrating Data Assimilation and Machine Learning, Deng Chao et al., Advances in Water Science, Vol. 34, No. 6, pp. 840-849" discloses a basin runoff simulation method that integrates data assimilation and machine learning. This existing technology uses a single LSTM deep learning model and a single input feature pattern, which can limit its ability to identify non-Gaussian characteristics in hydrological systems, potentially adversely affecting the accuracy of flood forecasting results.

[0008] In view of the above shortcomings, how to combine the existing state variable correction method with the deep learning model to alleviate the constraints of non-Gaussianity in the hydrological model in the real-time correction of flood forecasts and achieve high-precision real-time correction of flood forecasts is exactly the problem that the inventors need to solve. Summary of the Invention

[0009] Purpose of the invention: The purpose of the present invention is to address the deficiencies in the prior art and to provide a real-time correction method for flood simulation that takes into account the non-Gaussianity of the system.

[0010] Technical solution: A real-time correction method for flood simulation considering system non-Gaussianity of the present invention comprises the following steps:

[0011] Step 1: Construct a hydrological model M(.) for flood simulation in the basin to be simulated. The hydrological model M(.) has the function of calculating runoff and confluence; obtain the input observation rainfall and evaporation data U, the optimal parameter set and the initial value of the model state

[0012] Step 2: Construct an ensemble Kalman filter coupled with the hydrological model. The specific method is as follows:

[0013] Step 2.1: Before the hydrological model starts to simulate (t=0), the initial value of the model state Apply Gaussian perturbation and repeat N times to obtain the initial value of the model state after perturbation n∈[1,N]; then set up N hydrological models M that are independent of each other and calculated in parallel in time n (.); Then the corresponding input data U t＝1 and the optimal parameter set and the initial value of the model state after disturbance Input hydrological model M n (.), get the simulated flow rate and X1 and X2 values ​​of N groups at the first moment, recorded as

[0014] Step 2.2: At time t>1, calculate the predicted state of all N hydrological models M(.) and simulated flow ys n,t ;

[0015] Step 2.3: At time t>1, calculate the analytical state quantities of all N hydrological models M(.) at time t

[0016] Step 2.4: The optimal analytical value of the state variable at time t>1 Set as the ensemble average of the analyzed state quantities, the optimal simulated flow of the model at each moment Set to the ensemble mean of the simulated flow;

[0017] Step 3: Build a routing network model that integrates multiple deep learning technologies and soft gating modules. The specific method is as follows:

[0018] Step 3.1: Use three different deep learning models: long short-term memory network (LSTM), gated recurrent unit (GRU), and Transformer model as routing sub-models, denoted as M1_LSTM, M2_GRU, and M3_Transformer;

[0019] Step 3.2, set three different input modes INPUT_M1, INPUT_M2, and INPUT_M3 as the input feature sequences of the three routing sub-models M1_LSTM, M2_GRU, and M3_Transformer respectively;

[0020] Step 3.3: At time t, the input feature vector Z of the soft gating module is t Including flow observation value y t , soil moisture observation value s t , optimal simulated flow provided by the ensemble Kalman filter Optimal analysis value provided by the ensemble Kalman filter The ensemble mean of the predicted state variables provided by the hydrological model

[0021] Step 3.4: At time t, perform gated network calculation based on the two-layer perceptron network method to obtain the original score O of each routing sub-model weight. t ;

[0022] Step 3.5, use the softmax function to t Normalization is performed so that the sum of the weights of each routing sub-model is 1.

[0023] Step 4: Train the routing network model based on the hydrological model M(.) and the ensemble Kalman filter. The training objective is to maximize the sum of the Nash efficiency coefficient (NSE) of the simulated flow and the Nash efficiency coefficient (NSE) of the simulated soil moisture of the routing network model.

[0024] Step 5: Use the routing network model trained in step 4 to perform real-time correction on the forecast results of the hydrological model M(.).

[0025] Furthermore, the detailed steps of step 1 are as follows:

[0026] Step 1.1: Obtain flood data and evaporation data from the hydrological station at the outlet of the basin. Interpolate the flood data into a flood dataset with a time step of 1 hour. Evenly distribute the evaporation of the observation time step over time and interpolate it into an evaporation dataset with a time step of 1 hour.

[0027] Step 1.2: Obtain the coordinates and rainfall data of rain gauges within the basin. Evenly distribute the rainfall of the observation time step over time and interpolate it into a rainfall dataset with a time step of 1 hour. Then, use the Thiessen polygon method to divide the basin into k sub-basins centered on the rain gauge, where k is the number of rain gauges in the basin. The rainfall of each sub-basin is set as the rainfall dataset of the corresponding rain gauge in the basin, and the evaporation of each sub-basin is set as the evaporation dataset of the hydrological station at the basin outlet.

[0028] Step 1.3: Obtain the coordinates and soil moisture data of the soil moisture stations within the watershed. Temporally interpolate the soil moisture data of the observation time step to create a soil moisture dataset with a time step of 1 hour. Set the soil moisture of each sub-basin to the arithmetic mean of all soil moisture stations within the sub-basin. If there is no soil moisture station within the sub-basin, set the soil moisture of the sub-basin to the soil moisture value of the soil moisture station closest to the sub-basin.

[0029] Step 1.4. Select a hydrological model M(.) that has runoff generation and confluence calculation capabilities. The number of sub-basins calculated in the hydrological model M(.) is k. The state quantity set representing soil moisture in the runoff generation calculation is denoted as X1, and the state quantity set representing river water storage in the confluence calculation is denoted as X2. The model state quantity set is denoted as X = [X1, X2].

[0030] Step 1.5: At the first calculation time t = 1, the initial value of X1 is recorded as X1ini, and the value of each element in X1ini is set to half of the maximum soil moisture content; the initial value of X2 is recorded as X2ini, and the value of each element in X2ini is set to the observed flow rate of the outlet hydrological station at t = 1 divided by the number of sub-basins. The sum of the two is recorded as the initial value of the model state quantity

[0031] Step 1.6: Take the observed rainfall and evaporation as the model input data, denoted as U, set the insensitive parameters in the hydrological model M(.) to the recommended values, and use the shuffled complex evolution algorithm SCE-UA to automatically calibrate the parameters of the hydrological model M(.) with the objective function of maximizing the Nash efficiency coefficient between the observed flow and the simulated flow to obtain the optimal parameter set.

[0032] Furthermore, in step 2, during the process of constructing the ensemble Kalman filter coupled with the hydrological model,

[0033] For step 2.1, at time t=1, first initialize the model state Apply a Gaussian perturbation ω with a mean of 0 and a covariance of Q n , for N hydrological models M that are independent of each other and calculated in parallel in time n (.), let n = 1, 2, ..., N, and repeat the following process: input the model at this moment into the data U t＝1 and the optimal parameter set and the initial value of the model state after disturbance Input hydrological model M n (.), get the simulated flow of the nth group at the first moment and the initial X1 and X2 values ​​at the second moment, recorded as

[0034] For step 2.2, let n = 1, 2, ..., N and repeat the process until all N hydrological model simulation processes are completed at time t>1, and we get

[0035] ω n That is, the Gaussian noise with mean 0 and covariance Q added to the model state in step 2.1, and Q is set by expert experience;

[0036] In step 2.3, analyze the state quantity The calculation formula is as follows:

[0037]

[0038]

[0039]

[0040]

[0041] Where y t is the observed flow at time t, ys n,t is the flow simulated by the hydrological model at time t, K t is the Kalman gain matrix at time t, v is Gaussian noise with mean 0 and covariance R, and R is set by expert experience;

[0042] Optimal analysis value in step 2.4 And the optimal simulation flow The calculation formula is as follows:

[0043]

[0044]

[0045] Starting from t=2, repeat steps 2.2 to 2.4 at each moment until all moments are calculated.

[0046] Furthermore, in the process of constructing a routing network model integrating multiple deep learning technologies and soft gating modules in step 3:

[0047] In step 3.1, three deep learning models are selected: long short-term memory (LSTM), gated recurrent unit (GRU), and Transformer as routing sub-models, denoted as M1_LSTM, M2_GRU, and M3_Transformer.

[0048] Step 3.2, the input mode of M1_LSTM is INPUT_M1, with the flow observation value y t , optimal simulation flow As the input feature sequence of the routing sub-model, it focuses on extracting the dynamic error characteristics between the simulated traffic and the observed traffic, capturing the long-term trend and time series change law; the input mode of M2_GRU is INPUT_M2, specifically the traffic observation value y t , soil moisture observation value s t , optimal analysis value Predict the mean of the state quantity set As input, it focuses on the coupling relationship between model state and flow, and enhances the role of state variables in prediction; the input mode of M3_Transformer is INPUT_M3, specifically the flow observation value y t , soil moisture observation value s t , predicted state quantity And analyze the state As input, it focuses on extracting non-Gaussian features in the prediction results of different set members;

[0049] In step 3.3, at time t, the input feature vector Z of the soft gating module is t It consists of the following features: flow observation value y t , soil moisture observation value s t , optimal simulated flow provided by the ensemble Kalman filter Optimal analysis value provided by the ensemble Kalman filter The ensemble mean of the predicted state variables provided by the hydrological model

[0050] For step 3.4, at time t, the gated network calculation is performed based on the two-layer perceptron network method, and the input feature vector Z t First, after the nonlinear transformation of the hidden layer, the original score of each routing sub-model weight is further calculated through the output layer. t :

[0051] h t =σ(W h Z t +b h )

[0052] O t =σ(W o h t +b o )

[0053] Where W h is the hidden layer weight matrix, b h is the hidden layer bias vector, W o is the output layer weight matrix, b o is the output layer bias vector, σ(.) is the nonlinear ReLU activation function;

[0054] In step 3.5, use the softmax function to t Normalize the model so that the sum of the weights of each routing sub-model is 1:

[0055]

[0056] Where, ω k,t is the weight of routing sub-model k at time t, and K is the number of routing sub-models.

[0057] Furthermore, the routing network model training process in step 4 is as follows:

[0058] Step 4.1: The flow observation value y t, soil moisture observation value s t , optimal simulated flow provided by the ensemble Kalman filter Optimal analysis value provided by the ensemble Kalman filter The ensemble mean of the predicted state variables provided by the hydrological model Prediction state And analyze the state The time series of is input into each routing sub-model and soft gating module according to steps 3.2 and 3.3 respectively;

[0059] Step 4.2: Randomly select 80% of the total computation time as training set samples, use the Adam algorithm as the backpropagation optimization algorithm for the routing network model, and use the remaining 20% ​​as validation set samples;

[0060] Step 4.3: Simulated traffic ysl of the routing network model t and state quantity Xl t As the output value, the training objective is set to maximize the sum of the Nash efficiency coefficient (NSE) of the routing network model simulating traffic and the Nash efficiency coefficient (NSE) of the simulated soil moisture. Based on the performance of the validation set, the particle swarm algorithm is used to automatically optimize the hyperparameters (learning rate, number of hidden layers, hidden state / feature dimension, number of attention heads, Dropout ratio, batch size, and number of iterations) to determine the hyperparameters of the routing network model and complete the routing network model training.

[0061] Furthermore, the process of using the trained routing network model to perform real-time correction on the forecast results of the hydrological model M(.) in step 5 is as follows:

[0062] Step 5.1: When using the hydrological model M(.) for flood forecasting, at the forecast start time tf, repeat steps 2.2 to 2.4 to obtain the corresponding flow observation value y tf , soil moisture observation value s tf , optimal simulated flow provided by the ensemble Kalman filter Optimal analysis value provided by the ensemble Kalman filter The ensemble mean of the predicted state variables provided by the hydrological model Prediction state And analyze the state Then, according to steps 3.2 and 3.3, the corresponding features are input into each routing sub-model and soft gating module respectively, and the collective state variable of the prediction start time tf is obtained. N refers to N hydrological models Mn(.) that are independent of each other and calculated in parallel in time;

[0063] The aggregate state variables at time tf The calculation formula is as follows:

[0064]

[0065] Where, Xl tf It is the state quantity output by the routing network model at time tf, is the analytical state quantity at time tf, is the optimal analysis value of the Kalman filter at time tf;

[0066] Step 5.2: Set the forecast period to LT hours and The ensemble average of is the optimal initial state of the hydrological model at the start of the forecast, which is substituted into the hydrological model for flood forecasting.

[0067] Beneficial effects: Compared with the prior art, the present invention has the following advantages:

[0068] 1. This invention introduces deep learning technology into the real-time correction process of flood forecasting based on the ensemble Kalman filter. On the basis of ensuring the objectivity and reliability of the calculation results, it solves the problem of the conflict between the non-Gaussianity in the hydrological model and the Gaussianity assumption of the ensemble Kalman filter.

[0069] 2. The present invention uses observed rainfall, observed evaporation, observed soil moisture and observed flow, and the data source is stable and reliable. Compared with the single deep learning model based on a single input feature sequence, the routing network model based on multiple input feature sequences significantly improves the expression and processing capabilities of the system's non-Gaussianity, and gives full play to the complementary advantages of the three deep models of LSTM, GRU and Transformer in time series capture, and performs feature extraction on different input feature sequences, thereby effectively improving the model's information fusion and generalization capabilities.

[0070] 3. The present invention dynamically adjusts the output weights of each routing sub-model in real time through a soft gating mechanism, enabling the model to adapt to the nonlinear and dynamic change characteristics of different flood processes, significantly improving the robustness of the model and the real-time correction accuracy, and more effectively reducing the uncertainty of the ensemble prediction.

[0071] 4. The functional relationship between the variables in the present invention is clear, which is conducive to the computer automation execution of real-time correction of flood forecasts, is conducive to the use of ensemble Kalman filters to perform real-time correction of flood forecasts in hydrological models with non-Gaussianity, and is conducive to the in-depth development of digital hydrological forecasting research. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 Schematic diagram of the process of the present invention

[0073] Figure 2 The sub-basin distribution map of the Tunxi Basin in the embodiment

[0074] Figure 3 This is a flowchart of the routing network model training in the embodiment

[0075] Figure 4 The flood forecast results of Tunxi River Basin "TX-20150809082010041008" in the example are as follows: Figure 4 a is a real-time correction method using only ensemble Kalman filtering; Figure 4 b is the real-time correction effect of the scheme that captures non-Gaussianity using only a single LSTM model and a single input feature pattern; Figure 4 c is the real-time correction effect using the solution of the present invention;

[0076] Figure 5 This is the flood forecast result of Tunxi River Basin "TX-20110609082017062308" in the example. Figure 5 a is a real-time correction method using only ensemble Kalman filtering; Figure 5 b is the real-time correction effect of the scheme that captures non-Gaussianity using only a single LSTM model and a single input feature pattern; Figure 5 c is the real-time correction effect using the solution of the present invention. DETAILED DESCRIPTION

[0077] The technical solution of the present invention is described in detail below, but the protection scope of the present invention is not limited to the embodiments.

[0078] like Figure 1 As shown, the real-time correction method for flood simulation considering the non-Gaussianity of the system of the present invention comprises the following steps:

[0079] Step 1: Construct a hydrological model M(.) for flood simulation in the basin to be simulated. The hydrological model M(.) has the function of calculating runoff and confluence; obtain the input observation rainfall and evaporation data U, the optimal parameter set and the initial value of the model state

[0080] Step 2: Construct an ensemble Kalman filter coupled with the hydrological model. The specific method is as follows:

[0081] Step 2.1: Before the hydrological model starts to simulate (t=0), the initial value of the model state Apply Gaussian perturbation and repeat N times to obtain the initial value of the model state after perturbation n∈[1,N]; then set up N hydrological models M that are independent of each other and calculated in parallel in time n (.); Then the corresponding input data U t＝1 and the optimal parameter set and the initial value of the model state after disturbance Input hydrological model M n (.), get the simulated flow rate and X1 and X2 values ​​of N groups at the first moment, recorded as

[0082] Step 2.2: At time t>1, calculate the predicted state of all N hydrological models M(.) and simulated flow ys n,t ;

[0083] Step 2.3: At time t>1, calculate the analytical state quantities of all N hydrological models M(.) at time t

[0084] Step 2.4: The optimal analytical value of the state variable at time t>1 Set as the ensemble average of the analyzed state quantities, the optimal simulated flow of the model at each moment Set to the ensemble mean of the simulated flow;

[0085] Step 3: Build a routing network model that integrates multiple deep learning technologies and soft gating modules. The specific method is as follows:

[0086] Step 3.1: Use three different deep learning models: long short-term memory network (LSTM), gated recurrent unit (GRU), and Transformer model as routing sub-models, denoted as M1_LSTM, M2_GRU, and M3_Transformer;

[0087] Step 3.2, set three different input feature sequence modes INPUT_M1, INPUT_M2, and INPUT_M3 as the input feature sequences of the three routing sub-models M1_LSTM, M2_GRU, and M3_Transformer respectively. The input mode of M1_LSTM is INPUT_M1, specifically the traffic observation value y t , optimal simulation flow As the input feature sequence of the routing sub-model, it focuses on extracting the dynamic error characteristics between the simulated traffic and the observed traffic, capturing the long-term trend and time series change law; the input mode of M2_GRU is INPUT_M2, specifically the traffic observation value y t , soil moisture observation value s t , optimal analysis value Predict the mean of the state quantity set As input, it focuses on the coupling relationship between model state and flow, and enhances the role of state variables in prediction; the input mode of M3_Transformer is INPUT_M3, specifically the flow observation value y t, soil moisture observation value s t , predicted state quantity And analyze the state As input, it focuses on extracting non-Gaussian features in the prediction results of different set members;

[0088] Step 3.3: At time t, the input feature vector Z of the soft gating module is t It consists of the following features: flow observation value y t , soil moisture observation value s t , optimal simulated flow provided by the ensemble Kalman filter Optimal analysis value provided by the ensemble Kalman filter The ensemble mean of the predicted state variables provided by the hydrological model

[0089] Step 3.4: At time t, perform gated network calculation based on the two-layer perceptron network method to obtain the original score O of each routing sub-model weight. t ;

[0090] Step 3.5, use the softmax function to t Normalization is performed so that the sum of the weights of each routing sub-model is 1.

[0091] The routing network model in the present invention can effectively deal with the non-Gaussian problem of the system. By simultaneously introducing LSTM, GRU and Transformer models as routing sub-models, different input feature patterns are extracted and fused respectively, thereby enhancing the model's ability to capture the nonlinear characteristics of complex flood processes; at the same time, the soft gating mechanism is used to dynamically adjust the contribution weights of the sub-models in real time, so that the model has real-time adaptive adjustment capabilities, more accurately adapts to the complex time series characteristics and non-Gaussianity in the real-time flood forecasting process, and significantly improves the accuracy of real-time correction of flood simulation.

[0092] Step 4: Train the routing network model based on the hydrological model M(.) and the ensemble Kalman filter. The training objective is to maximize the sum of the Nash efficiency coefficient (NSE) of the simulated flow and the Nash efficiency coefficient (NSE) of the simulated soil moisture of the routing network model.

[0093] Step 5: Use the routing network model trained in step 4 to perform real-time correction on the forecast results of the hydrological model M(.). The specific method is as follows:

[0094] Step 5.1: When using the hydrological model M(.) for flood forecasting, at the forecast start time tf, repeat steps 2.2 to 2.4 to obtain the corresponding flow observation value y tf , soil moisture observation value s tf , optimal simulated flow provided by the ensemble Kalman filter Optimal analysis value provided by the ensemble Kalman filter The ensemble mean of the predicted state variables provided by the hydrological model Prediction state And analyze the state Then, according to steps 3.2 and 3.3, the corresponding features are input into each routing sub-model and soft gating module respectively, and the collective state variable of the prediction start time tf is obtained.

[0095] The aggregate state variables at time tf The calculation formula is as follows:

[0096]

[0097] Among them, Xl tf It is the state quantity output by the routing network model at time tf, is the analytical state quantity at time tf, is the optimal analysis value of the Kalman filter at time tf;

[0098] Step 5.2: Set the forecast period to LT hours and The ensemble average of is the optimal initial state of the hydrological model at the start of the forecast, which is substituted into the hydrological model for flood forecasting.

[0099] Example

[0100] This embodiment demonstrates real-time calibration of the Tunxi hydrological station at the outlet of the Tunxi Basin and the Xin'an River hydrological model. The specific method is as follows:

[0101] Step 1: Construct a hydrological model for flood simulation in the basin to be simulated:

[0102] Step 1.1: Obtain flood and evaporation data from the Tunxi hydrological station at the outlet of the Tunxi River Basin. Interpolate the flood data into a flood dataset with a time step of 1 hour. Evenly distribute the evaporation of the observation time step over time and interpolate it into an evaporation dataset with a time step of 1 hour.

[0103] Step 1.2: Obtain the coordinates and rainfall data of 9 rain gauges in the Tunxi Basin, evenly distribute the rainfall of the observation time step in time, interpolate it into a rainfall dataset with a time step of 1 hour, and use the Thiessen polygon method to divide the basin into 9 sub-basins with the rain gauge as the center (see Figure 2 ), the rainfall of each sub-basin is set to the rainfall dataset of the corresponding rain gauge station in the basin, and the evaporation of each sub-basin is set to the evaporation dataset of the hydrological station at the basin outlet;

[0104] Step 1.3: Obtain the coordinates and soil moisture data of the nine soil moisture stations within the Tunxi watershed. Interpolate the soil moisture data of the observation time step in time to obtain a soil moisture dataset with a time step of 1 hour. Set the soil moisture of each sub-basin as the arithmetic mean of all soil moisture stations within the sub-basin. If there is no soil moisture station within the sub-basin, set the soil moisture of the sub-basin to the soil moisture value of the soil moisture station closest to the sub-basin. The unit of the soil moisture observation value of the soil moisture station (cm 3 / cm 3 ) is inconsistent with the unit (mm) of the state quantity representing soil moisture in the Xin'anjiang model. The unit conversion of the soil moisture observation value at the soil moisture station is carried out using the following formula:

[0105]

[0106] Where s t is the soil moisture observation value at the soil moisture station after unit conversion at time t (mm), θ t is the soil moisture observation value at the soil moisture station at time t (cm 3 / cm 3 ), WM is the average tension water capacity of the basin, which is set to WM = 160 mm according to the Xinanjiang model manual, θ wp and θ fld are the wilting water content and field capacity of the soil, respectively, which are obtained from Table 2 in the literature "Anderson RM, Koren VI, Reed SM. Using SSURGO data to improve Sacramento Model a priori parameter estimates[J]. Journal of Hydrology, 2006, 320(1-2):103-116." according to the soil type at the soil moisture station.

[0107] Step 1.4: Let the Xin'an River hydrological model be denoted by M(.), the number of calculated sub-basins be 9, the tensile water content W be used as the state variable representing soil moisture, and the state variable set X1 = [W1, W2, ..., W9]; the outflow P at the outlet of each sub-basin be used as the state variable representing river water storage, and the state variable set be labeled X2 = [P1, P2, ..., P9]; the model state variable set is denoted by X = [X1, X2];

[0108] Step 1.5: At the first calculation time t = 1, the initial value of X1 is recorded as X1ini, and the value of each element in X1ini is set to half of the maximum soil moisture content, that is, 80 mm; the initial value of X2 is recorded as X2ini, and the value of each element in X2ini is set to the observed flow rate at the outlet hydrological station at t = 1 divided by the number of sub-basins, that is, 1.67 m3 / s, and the two are recorded together as the initial value of the model state quantity

[0109] Step 1.6: Use observed rainfall and evaporation as model input data, denoted as U.

[0110] According to the Xin'anjiang model manual, the insensitive parameters of the Xin'anjiang model are set as follows: upper layer tension water capacity WU = 20 mm, lower layer tension water capacity WL = 80 mm, deep evapotranspiration conversion coefficient C = 0.17, tension water storage capacity curve power B = 0.3, free water storage capacity curve power EX = 1.5, soil flow recession coefficient CI = 0.6, and Muskingum method calculation parameter XE = 0.4.

[0111] The sensitive parameters of the Xin'anjiang model include the evapotranspiration conversion coefficient K, the free water storage capacity SM, the free water outflow coefficient KI to the soil flow, the free water outflow coefficient KG to the groundwater runoff, the groundwater recession coefficient CG, the river network recession coefficient CS and the hysteresis L. The objective function is to maximize the Nash efficiency coefficient between the observed flow and the simulated flow. The shuffled complex evolution algorithm (abbreviated as SCE-UA algorithm) is used to automatically calibrate the sensitive parameters of the Xin'anjiang model. The optimal parameter set of the sensitive parameters of the Xin'anjiang model is obtained.

[0112] Step 2: Construct an ensemble Kalman filter coupled with the hydrological model:

[0113] Step 2.1: At t = 1, apply a Gaussian perturbation with a mean of 0 and a covariance of 5.0 to the initial value X1ini of the model state X1 of the Xinanjiang model, and apply a Gaussian perturbation with a mean of 0 and a covariance of 10.0 to the initial value X2ini of the model state X2. Repeat this 50 times to obtain the initial value of the model state after the perturbation. Set up 50 independent Xinanjiang models M that are calculated in parallel in time. n (.), let n = 1, 2, ..., 50 respectively, and repeat the following process:

[0114] The model input data U at this moment t＝1 and the optimal parameter set and the initial value of the model state after disturbance Input hydrological model M n (.), get the simulated flow of the nth group at the first moment and the initial X1 and X2 values ​​at the second moment, recorded as

[0115] Step 2.2: At time t (t>1), calculate the predicted state quantity through the hydrological model and simulated flow ysn,t , repeat each set (i.e., let n = 1, 2, ..., 50) until all 50 hydrological model simulation processes are calculated:

[0116]

[0117] Where, ω n To add Gaussian noise with mean 0 and covariance Q to the model state, all off-diagonal elements of Q are 0. For X1, all diagonal elements of Q are 5.0, and for X2, all diagonal elements of Q are 10.0.

[0118] Step 2.3, at time t (t>1), repeat each set (i.e., let n = 1, 2, ..., 50), and calculate the state quantity by the following formula

[0119]

[0120]

[0121]

[0122]

[0123] Where y t is the observed flow at time t, ys n,t is the simulated flow of the hydrological model at time t, K t is the Kalman gain matrix at time t, v is Gaussian noise with mean 0 and covariance R, since there is only one outlet of the Tunxi Basin, the matrix y t 、ys n,t 、 v and R are both 1 row and 1 column, and R is set to 10.0;

[0124] Step 2.4: Optimal analytical value of the state variable at time t (t>1) Set to analyze the state quantity to get the ensemble average value, the optimal simulation flow of the model at each moment Set to the ensemble average of the simulated flow, the expression is as follows:

[0125]

[0126]

[0127] Step 2.5: Starting from t=2, repeat steps 2.2 to 2.4 moment by moment until all moments are calculated.

[0128] Step 3: Build a routing network model that integrates multiple deep learning technologies and soft gating modules:

[0129] Step 3.1: Use three different deep learning models: long short-term memory network (LSTM), gated recurrent unit (GRU), and Transformer model as routing sub-models, denoted as M1_LSTM, M2_GRU, and M3_Transformer;

[0130] Step 3.2, set three different input feature sequence modes INPUT_M1, INPUT_M2, and INPUT_M3 as the input feature sequences of the three routing sub-models M1_LSTM, M2_GRU, and M3_Transformer respectively. The input mode of M1_LSTM is INPUT_M1, specifically the traffic observation value y t , optimal simulation flow As the input feature sequence of the routing sub-model, it focuses on extracting the dynamic error characteristics between the simulated traffic and the observed traffic, capturing the long-term trend and time series change law; the input mode of M2_GRU is INPUT_M2, specifically the traffic observation value y t , soil moisture observation value s t , optimal analysis value Predict the mean of the state quantity set As input, it focuses on the coupling relationship between model state and flow, and enhances the role of state variables in prediction; the input mode of M3_Transformer is INPUT_M3, specifically the flow observation value y t , soil moisture observation value s t , predicted state quantity And analyze the state As input, it focuses on extracting non-Gaussian features in the prediction results of different set members;

[0131] Step 3.3: At time t, the input feature vector Z of the soft gating module is t It consists of the following features: flow observation value y t , soil moisture observation value s t , optimal simulated flow provided by the ensemble Kalman filter Optimal analysis value provided by the ensemble Kalman filter The ensemble mean of the predicted state variables provided by the hydrological model

[0132] Step 3.4: At time t, perform gated network calculation based on the two-layer perceptron network method to obtain the original score O of each routing sub-model weight. t :

[0133] ht =σ(W h Z t +b h )

[0134] O t =σ(W o h t +b o )

[0135] Where W h is the hidden layer weight matrix, b h is the hidden layer bias vector, W o is the output layer weight matrix, b o is the output layer bias vector, σ(.) is the nonlinear ReLU activation function;

[0136] Step 3.5, use the softmax function to t Normalize the model so that the sum of the weights of each routing sub-model is 1:

[0137]

[0138] Where, ω k,t is the weight of routing sub-model k at time t, and K is the number of routing sub-models, K=3.

[0139] Step 4: Train the routing network model based on the hydrological model M(.) and the ensemble Kalman filter. The training process is shown in Figure 3 :

[0140] Step 4.1: The flow observation value y t , soil moisture observation value s t , optimal simulated flow provided by the ensemble Kalman filter Optimal analysis value provided by the ensemble Kalman filter The ensemble mean of the predicted state variables provided by the hydrological model Prediction state And analyze the state The time series of is input into each routing sub-model and soft gating module according to steps 3.2 and 3.3 respectively;

[0141] Step 4.2: Randomly select 80% of the total computation time as training set samples, use the Adam algorithm as the backpropagation optimization algorithm for the routing network model, and use the remaining 20% ​​as validation set samples;

[0142] Step 4.3: Simulated traffic ysl of the routing network model t and state quantity Xl tAs the output value, the state quantity Xl t Contains soil moisture X1l t The training objective is set to maximize the sum of the Nash efficiency coefficient (NSE) of the routing network model simulating traffic and the Nash efficiency coefficient (NSE) of the simulated soil moisture. Based on the performance of the validation set, the particle swarm algorithm is used to automatically optimize the hyperparameters (learning rate, number of hidden layers, hidden state / feature dimension, number of attention heads, Dropout ratio, batch size, and number of iterations) to determine the hyperparameters of the routing network model and complete the training of the routing network model. The final hyperparameter values ​​used are as follows: learning rate is 0.003, number of hidden layers is 3, hidden state / feature dimension is 128, number of attention heads is 4, Dropout ratio is 0.2, batch size is 64, and number of iterations is 500.

[0143] Step 5: Perform real-time correction on the hydrological model forecast results based on the trained routing network model:

[0144] Step 5.1: When using the hydrological model for flood forecasting, at the start of the forecast, denoted as time tf, repeat steps 2.2 to 2.4 to obtain the corresponding flow observation value y tf , soil moisture observation value s tf , optimal simulated flow provided by the ensemble Kalman filter Optimal analysis value provided by the ensemble Kalman filter The ensemble mean of the predicted state variables provided by the hydrological model Prediction state And analyze the state Then, according to steps 3.2 and 3.3, the corresponding features are input into each routing sub-model and soft gating module respectively, and the collective state variable of the prediction start time tf is obtained.

[0145]

[0146] Among them, Xl tf It is the state quantity output by the routing network model at time tf, is the analytical state quantity at time tf, is the optimal analysis value of the Kalman filter at time tf;

[0147] Step 5.2: If the forecast period is set to 6 hours, The ensemble mean of is the optimal initial state of the hydrological model at the start of the forecast. It is substituted into the hydrological model for flood forecasting. The forecast results of two floods in the Tunxi Basin with a forecast period of 6 hours are shown in Figure 4 and Figure 5 .

[0148] In the Tunxi River Basin "TX-2015080908" flood forecast of this embodiment, compared with the observed value, the flood peak relative error of the present invention is 12%, which is 39% lower than when no real-time correction is performed, 15% lower than when only the ensemble Kalman filter is used for real-time correction, and 8% lower than when only the LSTM model is used for real-time correction to capture non-Gaussianity; in the Tunxi River Basin "TX-2017062308" flood forecast, compared with the observed value, the flood peak relative error of the present invention is 1%, which is 32% lower than when no real-time correction is performed, 17% lower than when only the ensemble Kalman filter is used for real-time correction, and 7% lower than when only the LSTM model is used for real-time correction to capture non-Gaussianity, proving that the present invention is more beneficial to improving the accuracy of flood forecasting than the real-time correction method of flood forecasting using only the ensemble Kalman filter and the real-time correction method of capturing non-Gaussianity using only a single LSTM model and a single input feature pattern.

[0149] In summary, the present invention combines the update results of the ensemble Kalman filter with the routing network model, introduces a variety of different deep learning models and different input feature sequence patterns, and fully utilizes the advantages of the deep learning model to further process the non-Gaussianity that may exist in the hydrological model, thereby having a positive impact on the accuracy of flood forecast results.

Claims

1. A real-time correction method for flood simulation considering system non-Gaussianity, characterized in that: The following steps are involved: Step 1: Construct a hydrological model M(.) for flood simulation in the basin to be simulated. The hydrological model M(.) has the function of calculating runoff and confluence; obtain the input observation rainfall and evaporation data U, the optimal parameter set and the initial value of the model state Step 2: Construct an ensemble Kalman filter coupled with the hydrological model M(.). The specific method is as follows: Step 2.1: Before the hydrological model starts to simulate, the initial value of the model state Apply Gaussian perturbation and repeat N times to obtain the initial value of the model state after perturbation Then set up N hydrological models M that are independent of each other and calculated in parallel in time n (.); Then the corresponding input data U t＝1 and the optimal parameter set and the initial value of the model state after disturbance Input hydrological model M n (.), we get N groups of simulated flow, X1 and X2 values ​​at the first moment, and the vector after X1 and X2 are combined is recorded as Among them, X1 is the state quantity set of soil moisture in runoff calculation, and X2 is the state quantity set of river water storage in runoff calculation; Step 2.2: At time t>1, the observation data Ut and the optimal parameter group and Input hydrological model Mn(.) and calculate the predicted state quantity of all N hydrological models M(.) and simulated flow ys n,t ; Step 2.3: At time t>1, calculate the analytical state quantities of all N hydrological models M(.) at time t Step 2.4: The optimal analytical value of the state variable at time t>1 Set as the ensemble average of the analyzed state quantities, the optimal simulated flow of the model at each moment Set to the ensemble mean of the simulated flow; Step 3: Build a routing network model. The specific method is as follows: Step 3.1: Use the long short-term memory network LSTM, gated recurrent unit GRU, and Transformer models as routing sub-models, denoted as M1_LSTM, M2_GRU, and M3_Transformer respectively; Step 3.2, set three different input feature sequence modes INPUT_M1, INPUT_M2, and INPUT_M3 as the input feature sequences of the three routing sub-models M1_LSTM, M2_GRU, and M3_Transformer respectively; Step 3.3: At time t, the input feature vector Z of the soft gating module is t Includes the following features: flow observation value y t , soil moisture observation value s t , optimal simulated flow provided by the ensemble Kalman filter Optimal analysis value provided by the ensemble Kalman filter The ensemble mean of the predicted state variables provided by the hydrological model Step 3.4: At time t, perform gated network calculation based on the two-layer perceptron network method to obtain the original score O of each routing sub-model weight. t ; Step 3.5, use the softmax function to t Perform normalization so that the sum of the weights of each routing sub-model is 1; Step 4: Train the routing network model based on the hydrological model M(.) and the ensemble Kalman filter. The training objective is to maximize the sum of the Nash efficiency coefficient NSE of the routing network model simulating flow and the Nash efficiency coefficient NSE of the simulating soil moisture. Step 5: Use the routing network model trained in step 4 to perform real-time correction on the forecast results of the hydrological model M(.). The specific method is as follows: Step 5.1: When using the hydrological model M(.) for flood forecasting, at the forecast start time tf, repeat steps 2.2 to 2.4 to obtain the corresponding flow observation value y tf , soil moisture observation value s tf , optimal simulated flow provided by the ensemble Kalman filter Optimal analysis value provided by the ensemble Kalman filter The ensemble mean of the predicted state variables provided by the hydrological model Prediction state And analyze the state Then the corresponding features are input into each routing sub-model and soft gating module respectively, and the collective state variables at the prediction start time tf are obtained. The aggregate state variables at time tf The calculation formula is as follows: Among them, Xl tf It is the state quantity output by the routing network model at time tf, is the analytical state quantity at time tf, is the optimal analysis value of the Kalman filter at time tf; Step 5.2: Set the forecast period to LT hours and The ensemble average of is the optimal initial state of the hydrological model at the start of the forecast, which is substituted into the hydrological model for flood forecasting.

2. The real-time correction method for flood simulation considering system non-Gaussianity according to claim 1 is characterized in that: The detailed steps of step 1 are as follows: Step 1.1: Obtain flood data and evaporation data from the hydrological station at the outlet of the basin. Interpolate the flood data into a flood dataset with a time step of 1 hour. Evenly distribute the evaporation of the observation time step over time and interpolate it into an evaporation dataset with a time step of 1 hour. Step 1.2: Obtain the coordinates and rainfall data of rain gauges within the basin. Evenly distribute the rainfall of the observation time step over time and interpolate it into a rainfall dataset with a time step of 1 hour. Then, use the Thiessen polygon method to divide the basin into k sub-basins centered on the rain gauge, where k is the number of rain gauges in the basin. The rainfall of each sub-basin is set as the rainfall dataset of the corresponding rain gauge in the basin, and the evaporation of each sub-basin is set as the evaporation dataset of the hydrological station at the basin outlet. Step 1.3: Obtain the coordinates and soil moisture data of the soil moisture stations within the watershed. Temporally interpolate the soil moisture data of the observation time step to create a soil moisture dataset with a time step of 1 hour. Set the soil moisture of each sub-basin to the arithmetic mean of all soil moisture stations within the sub-basin. If there is no soil moisture station within the sub-basin, set the soil moisture of the sub-basin to the soil moisture value of the soil moisture station closest to the sub-basin. Step 1.

4. Select a hydrological model M(.) that has runoff generation and confluence calculation capabilities. The number of sub-basins calculated in the hydrological model M(.) is k. The state quantity set representing soil moisture in the runoff generation calculation is denoted as X1, and the state quantity set representing river water storage in the confluence calculation is denoted as X2. The model state quantity set is denoted as X = [X1, X2]. Step 1.5: At the first calculation time t = 1, the initial value of X1 is recorded as X1ini, and the value of each element in X1ini is set to half of the maximum soil moisture content; the initial value of X2 is recorded as X2ini, and the value of each element in X2ini is set to the observed flow rate of the outlet hydrological station at t = 1 divided by the number of sub-basins. The sum of the two is recorded as the initial value of the model state quantity Step 1.6: Take the observed rainfall and evaporation as the model input data, denoted as U, set the insensitive parameters in the hydrological model M(.) to the recommended values, and use the shuffled complex evolution algorithm SCE-UA to automatically calibrate the parameters of the hydrological model M(.) with the objective function of maximizing the Nash efficiency coefficient between the observed flow and the simulated flow to obtain the optimal parameter set.

3. The real-time correction method for flood simulation considering system non-Gaussianity according to claim 1 is characterized in that: In step 2, the process of constructing the ensemble Kalman filter coupled with the hydrological model is as follows: For step 2.1, at time t=1, first initialize the model state Apply a Gaussian perturbation ω with a mean of 0 and a covariance of Q n , for N hydrological models M that are independent of each other and calculated in parallel in time n (.), let n = 1, 2, ..., N, and repeat the following process: input the model at this moment into the data U t＝1 and the optimal parameter set and the initial value of the model state after disturbance Input hydrological model M n (.), get the simulated flow of the nth group at the first moment and the initial X1 and X2 values ​​at the second moment, recorded as For step 2.2, let n = 1, 2, ..., N and repeat the process until all N hydrological model simulation processes are completed at time t>1, and we get In step 2.3, analyze the state quantity The calculation formula is as follows: Where, K t is the Kalman gain matrix at time t, v is the observed data Gaussian noise with mean 0 and covariance R; Optimal analysis value in step 2.4 And the optimal simulation flow The calculation formula is as follows: Starting from t=2, repeat steps 2.2 to 2.4 at each moment until all moments are calculated.

4. The real-time correction method for flood simulation considering system non-Gaussianity according to claim 1 is characterized in that: The specific construction contents of the routing network model are as follows: Step 3.1, determine the input pattern and the corresponding input feature sequence; The input mode of M1_LSTM is INPUT_M1, with the flow observation value y t , optimal simulation flow As the input feature sequence of the routing sub-model; The input mode of M2_GRU is INPUT_M2, with the flow observation value y t , soil moisture observation value s t , optimal analysis value Predict the mean of the state quantity set As input; The input mode of M3_Transformer is INPUT_M3, with the flow observation value y t , soil moisture observation value s t , predicted state quantity And analyze the state is the input; Step 3.2: Determine the input feature vector Z of the soft gating module at time t t ; Step 3.3: The original score O of each routing sub-model weight at time t t The calculation formula is as follows: h t =σ(W h Z t +b h ) The t =σ(W o h t +b o ) Where W h is the hidden layer weight matrix, b h is the hidden layer bias vector, W o is the output layer weight matrix, b o is the output layer bias vector, σ(.) is the nonlinear ReLU activation function; Step 3.4, use the softmax function to t Normalize the weights of each routing sub-model so that the sum is 1. The calculation formula is as follows: Where, ω k,t is the weight of routing sub-model k at time t, and K is the number of routing sub-models.

5. The real-time correction method for flood simulation considering system non-Gaussianity according to claim 1 is characterized in that: The routing network model training process in step 4 is as follows: Step 4.1: The flow observation value y t , soil moisture observation value s t , optimal simulated flow provided by the ensemble Kalman filter Optimal analysis value provided by the ensemble Kalman filter The ensemble mean of the predicted state variables provided by the hydrological model Prediction state And analyze the state The time series of are input into each routing sub-model and soft gating module respectively; Step 4.2: Randomly select 80% of the total computation time as training set samples, use the Adam algorithm as the backpropagation optimization algorithm for the routing network model, and use the remaining 20% ​​as validation set samples; Step 4.3: Simulated traffic ysl of the routing network model t and state quantity Xl t As the output value, the state quantity Xl t Contains soil moisture X1l t The training objective is set to maximize the sum of the Nash efficiency coefficient NSE of the routing network model simulating traffic and the Nash efficiency coefficient NSE of simulating soil moisture. The particle swarm algorithm is used to automatically optimize the hyperparameters, determine the routing network model hyperparameters, and complete the routing network model training.

Images