Reservoir flow backstepping method based on ensemble Kalman filtering and hydrodynamic simulation

Abstract

The invention discloses a reservoir inflow back-stepping method based on ensemble Kalman filtering and hydrodynamic simulation. The method mainly comprises the following steps: (1) collecting and sorting reservoir upstream river topography, water level measured data and reservoir inflow data; (2) constructing and calibrating a one-dimensional hydrodynamic model according to reservoir area section topographic data and multi-section actually measured water levels; (3) assuming that the flow at the initial moment is in normal distribution and extracting an initial set, constructing a model instance based on the flow and the water level data at the last moment, and generating a state prior set; (4) generating an observation sample set based on the observation water level data and white noise; and (5) updating the storage flow posterior set by using set Kalman filtering, and repeating the steps (3)-(5) by taking the average set as the flow average value at the next moment until the last moment. According to the method, the influence of inaccurate water level measurement and the dynamic reservoir capacity of the reservoir can be reduced, minute-level, smooth and accurate reservoir inflow data can be obtained, and the method can be easily moved to different drainage basins to guide actual scheduling.

Description

Technical Field

[0001] The invention belongs to the field of reservoir regulation and relates to a reservoir inflow back-estimation method based on ensemble Kalman filtering and hydrodynamic simulation. Background Art

[0002] The calculation of inflow is one of the most important basic tasks in reservoir operation. The preparation of reservoir flood forecast and operation plan, the preparation of reservoir operation diagram, the economic evaluation of reservoir operation and the calculation of reservoir flood regulation all require inflow data. At present, the calculation of reservoir inflow mainly uses the static storage capacity method, which lacks consideration of the dynamic storage capacity of the reservoir. At the same time, the traditional method must take a longer calculation period step, and the inflow can only reach the hourly step, which simplifies the flood process and makes it difficult to reflect important information such as the flood peak. Coupled with water level observation errors, dynamic storage capacity influence, characteristic curve errors, and information collection errors, the inversion results of inflow are far from the real process. The unreasonable phenomena such as the jagged "fluctuation" of the reservoir inflow and the negative flow are the defects of the inversion method of inflow flow, which are the comprehensive manifestations of various errors such as large fluctuations in the measured water level data, poor representativeness, and inability to reflect the influence of dynamic storage capacity. The present invention selects a shorter calculation period step to make full use of all observed water level information, and adopts the ensemble Kalman filter method based on the continuity assumption of the inflow flow, and uses the one-dimensional hydrodynamic model as the observation equation of the ensemble Kalman filter to consider the influence of dynamic storage capacity. It is expected to obtain minute-level, smooth, and accurate reservoir inflow data, solve the problem of unreasonable phenomena such as sawtooth fluctuations and systematic errors in the inversion of reservoir inflow flow, and provide basic data support for hydrological model simulation and reservoir scheduling. Summary of the invention

[0003] Aiming at the problem that the reservoir inflow calculated based on the traditional water volume smoothing method has serious jagged fluctuations and low precision, the present invention proposes a reservoir inflow inversion method based on ensemble Kalman filtering and hydrodynamic simulation to improve the reservoir inflow accuracy and eliminate jagged fluctuations. The inflow time obtained by this method has higher accuracy and can make full use of measured data, which is of great significance for guiding reservoir scheduling and hydrological simulation.

[0004] In order to achieve the above purpose, the present method invents a method for inflow back-calculation based on ensemble Kalman filtering and hydrodynamic simulation, which includes the following steps: Step 1: Raw data processing: According to the research needs, collect and organize the measured data of the topography and water level of the river upstream of the reservoir.

[0005] Step 2: Model construction and calibration. A one-dimensional hydrodynamic model is constructed based on the measured cross-sectional data, and the measured water level data and historical runoff data are used as true values ​​to calibrate the model.

[0006] Step 3: Generate a priori set of states. Assume that the flow at the initial moment follows a normal distribution and extract the initial set. Construct multiple hydrodynamic model instances based on the flow and water level data at the previous moment to generate a priori set of states. Step 4: Generate observation samples. Generate an observation sample set based on the observed water level data and white noise.

[0007] Step 5: Ensemble Kalman Filter. Use the ensemble Kalman filter to update the posterior set of incoming traffic, take the average of the posterior set as the mean of the traffic at the next moment, repeat steps (3) to (5), and generate the reversed traffic in a rolling manner.

[0008] Furthermore, in step 2, original data including river channel topography data, cross-section roughness, historical water level, inflow and the like are collected. It should be noted that the river channel topography data must be underwater topography data, and the measured river channel cross-section data within the range from the reservoir dam site to the upstream backwater area should be collected from the local river management agency, and the distance between adjacent river channels should be as small as possible, preferably 1 to 3 km per cross-section. The river channel cross-section roughness data can be obtained by referring to the relevant literature of the reservoir to obtain the reference value and range of the river channel roughness data. Historical water level data refers to the measured water level data with high time resolution obtained from the water level measuring station in the reservoir area. The present invention is applicable to the situation where there are multiple water level measuring stations upstream of the reservoir area. The historical inflow flow data can adopt the inflow flow data reported historically by the reservoir.

[0009] The data in step 2 are further subjected to one-dimensional hydrodynamic modeling and model calibration. The specific steps include: The Saint-Venant equations are constructed to describe the flood evolution characteristics (cross-section water level, flow), and the river-type reservoir can be regarded as a one-dimensional unsteady flow. The Saint-Venant equations include two partial differential equations, namely the continuity equation and the momentum equation: Formula (1) Formula (2) In the formula, A is the cross-sectional area of ​​water flow; Q is the river flow; q It is the lateral inflow or outflow per unit length of a channel-type reservoir; x and t are independent spatial and temporal variables; g is the acceleration due to gravity; z It is the water depth of the cross-section; is the roughness coefficient; R is the hydraulic radius, ; It is wet week.

[0010] The Presman four-point eccentric implicit scheme is a commonly used numerical method for solving the Saint-Venant equation. It can calculate the flow and water level of all sections of the river in the one-dimensional hydrodynamic model. Figure 1 For simplicity, the differential equations between finite grid points are often used to replace the differential equations of continuous fluids. j and j +1 is the reservoir section j and section adjacent to the river-type reservoir j +1, grid point t and t +1 is the flood season time of river-type reservoirs t With the moment t +1. The time difference equation is shown in equation (3), and the distance difference equation is shown in equation (4). Figure 1 midpoint M Reservoir section j and cross section j +1 middle section j +1 / 2 at time The function of point M The variable value of can be calculated by formula (5). The weighted average of adjacent grid points is used for calculation, and the weight factor Reflects the moment when the variable value is biased t The smaller the weight factor, the more the Presman method tends to the known moment, which is an explicit finite difference method; the larger the weight factor, the more the Presman method tends to the unknown moment, which is an implicit finite difference method.

[0011] Every two adjacent river-type reservoir sections form a sub-river section, and each sub-river section has two equations, including a water level equation and a flow equation, corresponding to four variables. j +1 on t The flow and water level at time +1 can be calculated based on the reservoir section. j and cross section j +1 on t The flow and water level at the time are solved. Simplifying the Saint-Venant equations can obtain the equations in equation (6), which are actually continuous equations for the upper and lower sections. The parameters of the equations can be obtained by knowing the time t The parameters of are simplified as shown in formula (7). If the river-type reservoir is divided into N cross section, then there is N -1 sub-river section, corresponding to 2 N -2 equations, corresponding to 2 Nvariables. Since the number of variables is greater than the number of equations, two boundary conditions need to be added: the water level and flow rate of the initial section or the flow rate (inflow) of the initial section and the flow rate (outflow) of the terminal section. In this study, the flow rates of the initial section and the terminal section are selected as known boundary conditions for solution. After summarizing multiple equations, equation (8) is obtained. By solving it through the pursuit method, the flow rate and water level values ​​of each section of the river-type reservoir can be obtained.

[0012] Formula (3), Formula (4), Formula (5): Formula (6): In the formula, is the average water level of the j-th section at time t+1; is the average flow value of the j-th section at time t+1; Formula (7), In the formula, Calculate the time step for the model; is the section spacing; , are the cross-sectional area and cross-sectional width at point M respectively; , is the flow and water level at point M; Formula (8), In the formula, , are the water level and flow rate of the Nth section respectively; where, , , are the coefficients of the difference format of the continuity equation of the Nth section, , , , is the coefficient of the difference format of the Nth section momentum equation; Further, according to t All cross-section flow and water level at the moment and known boundary conditions ( t +1 time the flow of the initial section and the terminal section, that is, the reservoir t +1 time inflow and outflow) can be obtained t Unknown variables (flow and water level) at each section at time +1.

[0013] For the case of lateral inflow, a virtual river section can be set , then the basic continuity equation is: Formula (9): In the formula, , are the water levels at upstream section i and downstream section i+1 respectively; , are the flow rates at upstream section i and downstream section i+1 respectively; It is the lateral inflow between sections i to i+1.

[0014] By replacing (6) with (9), the recursive relationship can also be obtained.

[0015] A one-dimensional hydrodynamic model can be constructed based on the sorted river section data. The roughness value of each section of the river channel is manually adjusted based on the reference value and the measured water level data, and can also be calibrated with the help of data assimilation methods or other optimization methods.

[0016] Further, the prior set required for the ensemble Kalman filter is generated according to the hydraulic model in step 2. The specific steps are: From the data collection in step 1, it can be known that the observed quantity in the present invention is the measured water level data. When the ensemble Kalman filter model is used to correct the inflow, it is necessary to generate a priori set of system state simulation values ​​based on the priori set of parameter samples. The present invention uses the inflow as a parameter. Based on the continuity assumption of the inflow, that is, the inflow at the initial moment and the final moment in a short period of time is approximately equal.

[0017] Therefore, the following parameter state transfer equation and system state transfer equation can be obtained: Formula (10), Formula (11), In the formula, is the error term, which is assumed to obey a Gaussian distribution with a mean of 0 and a constant variance. Formula (10) shows that the change of inflow in a short period of time is limited. is the system status water level, is the system evolution noise, It is a one-dimensional hydrodynamic model.

[0018] For one-dimensional hydrodynamic models, a warm-up period needs to be set because the one-dimensional hydrodynamic model has many parameters and initial states. There may be water level and flow fluctuations in the early stage of the simulation. This problem can be solved by adding a warm-up period, using historical water level and flow data as input, simulating to a stable moment, and then using the randomly generated inflow at the next moment as input to obtain the water level data at the next moment. The length of the warm-up period needs to be determined based on the number of sections and the initial condition settings of the model, and is generally appropriate for 24 hours. The generated flow set is used as the input of the preheated hydrodynamic model to obtain a simulated water level data set as a priori set of system states, and the produced parameter samples are used as a priori set of parameters for data assimilation.

[0019] Further, observation samples are generated. The observation value in the present invention is the actual water level of the water level station. Since the water level station has measurement errors and the measured water level is a point, it is necessary to consider the error of the actual water level and generate observation samples of the average water level of the section.

[0020] Therefore, the observation sample can be generated by the following formula: Formula (12), In the formula, represents the measured water level at the i-th water level station; Represents the reservoir water level observation error, which follows a Gaussian distribution with a mean of 0. The variance can be adjusted as a parameter, and the generally recommended value is 0.01~0.1m.

[0021] Further, based on the prior set and observation samples generated in steps 3 and 4, the ensemble Kalman filter model is used to roll-update the inflow flow. The key to applying ensemble Kalman filtering is to construct the system state transfer equation and observation equation. The general state transfer equation is: Formula (13), In the formula, is the system state variable; is the system parameter; is the independent white noise of the system, with mean 0 and variance The normal distribution of To simulate the water level, Inbound flow.

[0022] The observation equation is: Formula (14), In the formula, is the observation value vector; Represents the conversion relationship between system state variables and observation values; is the system white noise, with a mean of 0 and a variance of The normal distribution of The actual water level.

[0023] Based on the above system state transfer equation and observation equation, the assimilation process of ensemble Kalman filtering is as follows: Formula (15), Formula (16), Formula (17), In the formula, is the predicted value of the kth set number at time i+1; is the updated value of the kth set number at time i; is the kth set of white noise; is the observed value of the kth set number; is the observation error of the kth set number; is the gain factor, which represents the weight relationship between the predicted value and the observed value. Its calculation formula is as follows: Formula (18), Formula (19), Formula (20), In the formula, is the predicted state variable covariance matrix; is the observed predicted value error covariance matrix; where: , , is the ensemble mean of the predicted state variables, is the ensemble mean of the observed predicted values; N is the number of ensembles; and the superscript T denotes the matrix transpose.

[0024] The state equation and observation equation in the present invention are equations (10) and (11), and the main update parameter Q, the set Kalman gain and the parameter update equation are as follows: Formula (21), Formula (22), Formula (23), Where, Q is the flow parameter; is the gain factor of parameter Q at time i+1; is the covariance matrix of the prediction parameter variables; is the prior value of the kth set number Q; is the posterior value of the kth set number Q; Substitute the prior set generated in step 3 and the observation set generated in step 4 into equations (18) to (22) to calculate the Kalman gain: , , and then substitute the Kalman gain into equations (16) and (23) to obtain the posterior set, calculate the ensemble average of the posterior set to obtain the parameter estimate, and calculate the ensemble average of the system state to obtain the posterior value of the water level of each river section.

[0025] The parameter estimation value is the inflow required by the present invention, and the water level posterior value is the state data with a smaller variance obtained by combining simulation and measured data. The updated water level data and inflow are input into the one-dimensional model as the data of the previous moment, and the inflow at the next moment is regenerated as a parameter. Steps 3 to 5 are repeated for rolling solution to obtain the inflow data at continuous moments.

[0026] Compared with the prior art, the advantages and beneficial effects of the present invention are as follows: The inflow inversion method based on ensemble Kalman filtering and one-dimensional hydrodynamic simulation provided by the present invention fully considers the influence of dynamic storage capacity caused by uneven water surface line during flood evolution, and improves the calculation accuracy of inflow flow; uses minute-level measured water level data for data assimilation, makes full use of measured data information, and can well restore the flood peak process; the present invention adopts the flow continuity assumption, and the inflow flow process derived has the characteristics of small fluctuation and good continuity, which can better guide reservoir scheduling. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] Figure 1 Schematic diagram of Pressman's method; Figure 2 Schematic diagram of side inflow; Figure 3 The results of the calibration of the Shibaozhai section; Figure 4 Yunyang section calibration results; Figure 5 Wushan section calibration results; Figure 6 The results of the Zigui section calibration; Figure 7 Comparison of incoming traffic Figure 8 Flow chart of the method of the present invention. DETAILED DESCRIPTION

[0028] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and examples. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0029] like Figure 1 As shown, the present invention includes the following process: Step 1: Raw data processing: According to the research needs, collect and organize the measured data of the topography and water level of the river upstream of the reservoir.

[0030] Step 2: Model construction and calibration. A one-dimensional hydrodynamic model is constructed based on the measured cross-sectional data, and the measured water level data and historical runoff data are used as true values ​​to calibrate the model.

[0031] Step 3: Generate a priori set of states. Assume that the flow at the initial moment follows a normal distribution and extract the initial set, build multiple hydrodynamic model instances based on the flow and water level data at the previous moment, and generate a priori set of states.

[0032] Step 4: Generate observation samples. Generate an observation sample set based on the observed water level data and white noise.

[0033] Step 5: Ensemble Kalman Filter. Use the ensemble Kalman filter to update the posterior set of incoming traffic, take the average of the posterior set as the mean of the traffic at the next moment, repeat steps (3) to (5), and generate the reversed traffic in a rolling manner.

[0034] Example 1 The present invention takes the Three Gorges Dam as an example. The Three Gorges Reservoir is a typical river-type reservoir with a huge dynamic storage capacity. Its inflow cannot be obtained through actual measurement or traditional static storage capacity estimation methods. The inflow inversion method based on ensemble Kalman filtering and hydrodynamic simulation provided by the present invention is used.

[0035] In step 1, according to the research needs, collect and organize the measured data of river topography, inflow and water level upstream of the reservoir. First, collect a total of 287 sections from the Three Gorges Dam site to the upstream Cuntan section, and the distance between adjacent sections is calculated by the center coordinates of the river section. By consulting the Three Gorges Dispatching Center, collect 2-hour inflow data in August 2020 and October 2022, 15-minute interval measured water level data of 8 water level stations, and the recommended value of river section roughness of 0.03~0.08. If there are individual data missing or incomplete information, they are interpolated and extended by linear interpolation.

[0036] In the present invention, in step 2, a one-dimensional hydrodynamic model is constructed based on the Saint-Venant equation and manually calibrated.

[0037] The specific steps are as follows: first, construct the Saint-Venant equations to describe the flood evolution characteristics (cross-section water level, flow), and the river-type reservoir can be regarded as a one-dimensional unsteady flow. The Saint-Venant equations include two partial differential equations, namely the continuity equation and the momentum equation: Formula (1), Formula (2): In the formula, A is the cross-sectional area of ​​water flow; Qis the river flow; q It is the lateral inflow or outflow per unit length of a channel-type reservoir; x and t are independent spatial and temporal variables; g is the acceleration due to gravity; z It is the water depth of the cross-section; is the roughness coefficient; R is the hydraulic radius, ; It is wet week.

[0038] The Presman four-point eccentric implicit scheme is a commonly used numerical method for solving the Saint-Venant equation. It can calculate the flow and water level of all sections of the river in the one-dimensional hydrodynamic model. Figure 1 For simplicity, the differential equations between finite grid points are often used to replace the differential equations of continuous fluids. j and j +1 is the reservoir section j and section adjacent to the river-type reservoir j +1, grid point t and t +1 is the flood season time of river-type reservoirs t With the moment t +1. The time difference equation is shown in equation (3), and the distance difference equation is shown in equation (4). Figure 1 midpoint M Reservoir section j and cross section j +1 middle section j +1 / 2 at time The function of point M The variable value of can be calculated by formula (5). The weighted average of adjacent grid points is used for calculation, and the weight factor Reflects the moment when the variable value is biased t The smaller the weight factor, the more the Presman method tends to the known moment, which is an explicit finite difference method; the larger the weight factor, the more the Presman method tends to the unknown moment, which is an implicit finite difference method.

[0039] Every two adjacent river-type reservoir sections form a sub-river section, and each sub-river section has two equations, including a water level equation and a flow equation, corresponding to four variables. j +1 on t The flow and water level at time +1 can be calculated based on the reservoir section. j and cross section j +1 on tThe flow and water level at the time are solved. Simplifying the Saint-Venant equations can obtain the equations in equation (6), which are actually continuous equations for the upper and lower sections. The parameters of the equations can be obtained by knowing the time t The parameters of are simplified as shown in formula (7). If the river-type reservoir is divided into N cross section, then there is N -1 sub-river section, corresponding to 2 N -2 equations, corresponding to 2 N variables. Since the number of variables is greater than the number of equations, two boundary conditions need to be added: the water level and flow rate of the initial section or the flow rate (inflow) of the initial section and the flow rate (outflow) of the terminal section. In this study, the flow rates of the initial section and the terminal section are selected as known boundary conditions for solution. After summarizing multiple equations, equation (8) is obtained. By solving it through the pursuit method, the flow rate and water level values ​​of each section of the river-type reservoir can be obtained.

[0040] Formula (3), Formula (4), Formula (5): Formula (6): In the formula, is the average water level of the j-th section at time t+1; is the average flow value of the j-th section at time t+1; Formula (7), In the formula, Calculate the time step for the model; is the section spacing; , are the cross-sectional area and cross-sectional width at point M respectively; , is the flow rate and water level at point M.

[0041] Formula (8), In the formula, , are the water level and flow rate of the Nth section respectively; where, , , are the coefficients of the difference format of the continuity equation of the Nth section, , , , are the coefficients of the difference format of the Nth section momentum equation.

[0042] Further, according to tAll cross-section flow and water level at the moment and known boundary conditions ( t +1 time the flow of the initial section and the terminal section, that is, the reservoir t +1 time inflow and outflow) can be obtained t Unknown variables (flow and water level) at each section at time +1.

[0043] For the case of side inflow, such as Figure 2 As shown, a virtual river section can be set , then the basic continuity equation is: Formula (9): In the formula, , are the water levels at upstream section i and downstream section i+1 respectively; , are the flow rates at upstream section i and downstream section i+1 respectively; It is the lateral inflow between sections i to i+1.

[0044] By replacing (6) with (9), the recursive relationship can also be obtained.

[0045] A one-dimensional hydrodynamic model can be constructed based on the sorted river section data. The roughness value of each section of the river channel is manually adjusted based on the reference value and the measured water level data, and can also be calibrated with the help of data assimilation methods or other optimization methods.

[0046] The roughness calibration results are shown in Table 1.

[0047] Table 1 River channel roughness coefficient flood Baishatuo Shibaozhai Wanxian Yunyang Fengjie Wushan Padang Zigui 2020.8 0.0294 0.0404 0.0400 0.0297 0.0710 0.0771 0.0631 0.0598 The roughness coefficients of other sections along the Three Gorges Reservoir are determined by distance linear interpolation.

[0048] The water level determination results of some sections are shown in Figures (3~6).

[0049] In the present invention, in step 3, it is assumed that the flow at the initial moment obeys a normal distribution and an initial set is extracted, and multiple hydrodynamic model instances are constructed based on the flow and water level data at the previous moment to generate a state priori set.

[0050] Step 3 specifically includes: Based on the continuity assumption of the inflow flow, that is, the inflow flow at the initial moment and the final moment in a short period of time is approximately equal.

[0051] Therefore, the following parameter state transfer equation and system state transfer equation can be obtained: Formula (10), Formula (11), In the formula, is the error term, which is assumed to obey a Gaussian distribution with a mean of 0 and a constant variance. Formula (10) shows that the change of inflow in a short period of time is limited. is the system status water level, is the system evolution noise, It is a one-dimensional hydrodynamic model.

[0052] For the one-dimensional hydrodynamic model, a warm-up period needs to be set, because the one-dimensional hydrodynamic model has many parameters and initial states. There may be water level and flow fluctuations in the early stage of the simulation. This problem can be solved by increasing the warm-up period, using historical water level and flow data as input, simulating to a stable moment, and then using the randomly generated inflow flow at the next moment as input to obtain the water level data at the next moment. The length of the warm-up period needs to be determined based on the number of sections and the initial condition settings of the model. In the embodiment of the present invention, it is set to 24 hours. The generated flow set is used as the input of the hydrodynamic model after preheating, and a simulated water level data set is obtained as a system state prior set. The produced parameter samples are used as parameter prior sets for data assimilation.

[0053] In the present invention, in step 4, an observation sample is further generated. The observation value in the present invention is the measured water level of the water level station. Since the water level station has measurement errors and the measured water level is a point, it is necessary to consider the error of the measured water level and generate an observation sample of the average water level of the section.

[0054] Therefore, the observation sample can be generated by the following formula: Formula (12), In the formula, represents the measured water level at the i-th water level station; Represents the reservoir water level observation error, which follows a Gaussian distribution with a mean of 0. The variance can be adjusted as a parameter, and the generally recommended value is 0.01~0.1m.

[0055] In the present invention, in step 5, the ensemble Kalman filter is used to update the posterior set of incoming traffic.

[0056] The key to applying ensemble Kalman filtering is to construct the system state transfer equation and observation equation. The general state transfer equation is: Formula (13), In the formula, is the system state variable; is the system parameter; is the independent white noise of the system, with mean 0 and variance The normal distribution of To simulate the water level, Inbound flow.

[0057] The observation equation is: Formula (14), In the formula, is the observation value vector; Represents the conversion relationship between system state variables and observation values; is the system white noise, with a mean of 0 and a variance of The normal distribution of The actual water level.

[0058] Based on the above system state transfer equation and observation equation, the assimilation process of ensemble Kalman filtering is as follows: Formula (15), Formula (16), Formula (17), In the formula, is the predicted value of the kth set number at time i+1; is the updated value of the kth set number at time i; is the kth set of white noise; is the observed value of the kth set number; is the observation error of the kth set number; is the gain factor, which represents the weight relationship between the predicted value and the observed value. Its calculation formula is as follows: Formula (18), Formula (19), Formula (20), In the formula, is the predicted state variable covariance matrix; is the observed predicted value error covariance matrix; where: , , is the ensemble mean of the predicted state variables, is the ensemble mean of the observed predicted values; N is the number of ensembles; and the superscript T denotes the matrix transpose.

[0059] The state equation and observation equation in the present invention are equations (10) and (11), and the main update parameter Q, the set Kalman gain and the parameter update equation are as follows: Formula (21), Formula (22), Formula (23), Where, Q is the flow parameter; is the gain factor of parameter Q at time i+1; is the covariance matrix of the prediction parameter variables; is the prior value of the kth set number Q; is the posterior value of the kth set number Q; Substitute the prior set generated in step 3 and the observation set generated in step 4 into equations (18) to (22) to calculate the Kalman gain: , , and then substitute the Kalman gain into equations (16) and (23) to obtain the posterior set, calculate the ensemble average of the posterior set to obtain the parameter estimate, and calculate the ensemble average of the system state to obtain the posterior value of the water level of each river section.

[0060] The parameter estimation value is the inflow required by the present invention, and the water level posterior value is the state data with a smaller variance obtained by combining simulation and measured data. The updated water level data and inflow are input into the one-dimensional model as the data of the previous moment, and the inflow at the next moment is regenerated as a parameter. Steps 3 to 5 are repeated for rolling solution to obtain the inflow data at continuous moments.

[0061] The results of the inflow flow calculated by the present invention are compared with those calculated by the traditional water balance method. Figure 7 shown.

[0062] Based on the above conclusions, the inflow inversion method based on ensemble Kalman filtering and one-dimensional hydrodynamic simulation provided by the present invention fully considers the influence of dynamic storage capacity caused by uneven water surface line during flood evolution, and improves the calculation accuracy of inflow. The present invention uses minute-level measured water level data for data assimilation, makes full use of measured data information, and can restore the flood peak process well. The inferred inflow process has the characteristics of small fluctuation, good continuity, and fine time scale, which can better guide reservoir scheduling.

[0063] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the methods and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for inflow inference based on ensemble Kalman filtering and hydrodynamic simulation, characterized by: The method comprises the following steps: Step 1: Raw data collection: Collect and organize river topography data, historical water levels, and inflow data; Step 2: Model construction and calibration: A one-dimensional hydrodynamic model is constructed based on the measured cross-section data, and the measured water level data and historical runoff data are used as true values ​​to calibrate the model; Step 3: Generate a priori set of states: Assume that the flow at the initial moment obeys a normal distribution and extract the initial set, build multiple hydrodynamic model instances based on the flow and water level data at the previous moment, and generate a priori set of states; Step 4: Generate observation samples, generate observation sample sets based on observed water level data and white noise; Step 5: Ensemble Kalman Filter: Use the ensemble Kalman filter to update the posterior set of incoming traffic, take the average of the posterior set as the mean of the traffic at the next moment, repeat steps (3) to (5), and generate the reversed traffic in a rolling manner.

2. The method for inflow inference based on ensemble Kalman filtering and hydrodynamic simulation according to claim 1 is characterized by: In step 1, The river channel topography data is underwater topographic data, which collects measured river channel cross-section data from the reservoir dam site to the upstream backwater area. The distance between adjacent river channels is 1~3km as one cross-section. Obtain reference values ​​and ranges of river channel roughness data; Historical water level data refers to the measured water level data with high time resolution obtained from water level measuring stations in the reservoir area; The inflow data adopts the inflow data reported historically by the reservoir.

3. The method for inflow inference based on ensemble Kalman filtering and hydrodynamic simulation according to claim 1 is characterized by: In step 2, the specific process of model construction and calibration is as follows: 2.1 The evolution characteristics of river-type reservoir floods are described by the Saint-Venant equations. They are the continuity equation and the momentum equation: Formula (1); Formula (2); In the formula, A is the cross-sectional area of ​​water flow; Q is the river flow; q It is the lateral inflow or outflow per unit length of a channel-type reservoir; x and t are independent spatial and temporal variables; g is the acceleration due to gravity; z It is the water depth of the cross-section; is the roughness coefficient; R is the hydraulic radius, ; It is wet week; 2.2 Pressman four-point eccentric implicit method is used to calculate the flow and water level of all sections of the river in the one-dimensional hydrodynamic model. j and j +1 is the reservoir section j and section adjacent to the river-type reservoir j +1, grid point t and t +1 is the flood season time of river-type reservoirs t With the moment t +1, the time difference equation is shown in equation (3), the distance difference equation is shown in equation (4), point M Reservoir section j and cross section j +1 middle section j +1 / 2 at time The function of point M The variable value of is calculated by formula (5), using the weighted average of adjacent grid points for calculation, and the weight factor Indicates the time when the variable value is biased t the extent of; 2.3 Every two adjacent river-type reservoir sections form a sub-river section. Each sub-river section has two equations, including a water level equation and a flow equation, corresponding to four variables, reservoir section j and section j +1 on t The flow and water level at time +1 are based on the reservoir section j and cross section j +1 on t The flow and water level at the time are solved to obtain the equation group in equation (6). The equation group is the continuous equation of the upper and lower sections. The parameters of the equation group are obtained by the known time t The parameters of are simplified as shown in formula (7); Formula (3), Formula (4), Formula (5): Formula (6): In the formula, is the average water level of the j-th section at time t+1; is the average flow value of the j-th section at time t+1; Formula (7), In the formula, Calculate the time step for the model; is the section spacing; , are the cross-sectional area and cross-sectional width at point M respectively; , is the flow rate and water level at point M.

4. The method for inflow inversion based on ensemble Kalman filtering and hydrodynamic simulation according to claim 3 is characterized by: 2.4 River-type reservoirs are divided into N cross section, then there is N -1 sub-river section, corresponding to 2 N -2 equations, corresponding to 2 N variables, and two boundary conditions are added, namely, the water level and flow of the initial section or the flow-inflow of the initial section and the flow-outflow of the terminal section; the flow of the initial section and the terminal section is selected as the known boundary conditions for solution, and the equations are summarized to obtain equation (8), which is solved by the pursuit method to obtain the flow and water level values ​​of each section of the river-type reservoir; Formula (8), In the formula, , are the water level and flow rate of the Nth section respectively; where, , , are the coefficients of the difference format of the continuity equation of the Nth section, , , , is the coefficient of the difference format of the Nth section momentum equation; 2.5 According to t All cross-section flow and water level at the moment and known boundary conditions ( t +1 time the flow of the initial section and the terminal section, that is, the reservoir t +1 time inflow and outflow) to obtain t +1 The unknown variables of each section are flow and water level, For the case of lateral inflow, a virtual river section is set , the basic continuity equation is: Formula (9): In the formula, , are the water levels at upstream section i and downstream section i+1 respectively; , are the flow rates at upstream section i and downstream section i+1 respectively; It is the lateral inflow between sections i to i+1.

5. The method for inflow inference based on ensemble Kalman filtering and hydrodynamic simulation according to claim 4 is characterized by: 2.6 A one-dimensional hydrodynamic model is constructed based on the sorted river section data. The roughness value of each section is manually adjusted based on the reference value and the measured water level data, or calibrated with the help of data assimilation methods or other optimization methods.

6. The method for inflow inference based on ensemble Kalman filtering and hydrodynamic simulation according to claim 1 is characterized by: In the step S3, Taking the inflow flow as a parameter, the inflow flow at the initial moment and the final moment in a short period of time is equal, and the following parameter state transfer equation and system state transfer equation are obtained: Formula (10), Formula (11), In the formula, is the error term, which is assumed to obey a Gaussian distribution with a mean of 0 and a constant variance; is the system status water level, Evolution noise for the system; The generated flow set is used as the input of the preheated hydrodynamic model to obtain the simulated water level data set as the system state prior set, and the produced parameter samples are used as the parameter prior set for data assimilation.

7. The method for inflow inference based on ensemble Kalman filtering and hydrodynamic simulation according to claim 1 is characterized by: In the step S4, The observation samples are generated by the following formula: Formula (12), In the formula, represents the measured water level at the i-th water level station; Represents the reservoir water level observation error, obeys a Gaussian distribution with a mean of 0, and the variance is adjusted as a parameter. The recommended value is 0.01~0.1m.

8. The method for inflow inference based on ensemble Kalman filtering and hydrodynamic simulation according to claim 1 is characterized by: In the specific operation of the ensemble Kalman filter described in step 5: The state transfer equation is: Formula (13), In the formula, To simulate the water level, is the inbound flow; is the independent white noise of the system, with mean 0 and variance Normal distribution of The observation equation is: Formula (14), In the formula, is the measured water level; Represents the conversion relationship between system state variables and observation values; is the system white noise, with a mean of 0 and a variance of The normal distribution of .

9. The method for inflow inference based on ensemble Kalman filtering and hydrodynamic simulation according to claim 8 is characterized by: Based on the above system state transfer equation and observation equation, the assimilation process of ensemble Kalman filtering is as follows: Formula (15), Formula (16), Formula (17), In the formula, is the predicted value of the kth set number at time i+1; is the updated value of the kth set number at time i; is the kth set of white noise; is the observed value of the kth set number; is the observation error of the kth set number; is the gain factor, which indicates the weight relationship between the predicted value and the observed value. Its calculation formula is as follows: Formula (18), Formula (19), Formula (20), In the formula, is the predicted state variable covariance matrix; is the observed predicted value error covariance matrix; where: , , is the ensemble mean of the predicted state variables, is the ensemble mean of the observed predicted values; N is the number of ensembles; the superscript T indicates the matrix transpose; The main update parameter Q, the set Kalman gain and the parameter update equation are as follows: Formula (21), Formula (22), Formula (23), Where, Q is the flow parameter; is the gain factor of parameter Q at time i+1; is the covariance matrix of the prediction parameter variables; is the prior value of the kth set number Q; is the posterior value of the kth set number Q; Substitute the prior set generated in step 3 and the observation set generated in step 4 into equations (18) to (22) to calculate the Kalman gain: , , and then substitute the Kalman gain into equations (16) and (23) to obtain the posterior set, calculate the ensemble average of the posterior set to obtain the parameter estimate, and calculate the ensemble average of the system state to obtain the posterior value of the water level of each river section.

10. The method for inflow inversion based on ensemble Kalman filtering and hydrodynamic simulation according to claim 9, characterized in that: The parameter estimation value is the required inflow, and the water level posterior value is the state data with smaller variance obtained by combining simulation and measured data. The updated water level data and inflow are input into the one-dimensional model as the data of the previous moment, and the inflow at the next moment is regenerated as a parameter. Steps 3 to 5 are repeated for rolling solution to obtain the inflow data at consecutive moments.

Images