Dual state-parameter optimal estimation of one-dimensional open channel model using ensemble Kalman filter*

LAI Rui-xun (赖瑞勋), FANG Hong-wei (方红卫), HE Guo-jian (何国建)

Department of Hydraulic Engineering and State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China

Yellow River Institute of Hydraulic Research, Zhengzhou 450003, China, E-mail: lairuixun@163.com

YU Xin (余欣), YANG Ming (杨明), WANG Ming (王明)

Yellow River Institute of Hydraulic Research, Zhengzhou 450003, China

Dual state-parameter optimal estimation of one-dimensional open channel model using ensemble Kalman filter*

LAI Rui-xun (赖瑞勋), FANG Hong-wei (方红卫), HE Guo-jian (何国建)

Department of Hydraulic Engineering and State Key Laboratory of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China

Yellow River Institute of Hydraulic Research, Zhengzhou 450003, China, E-mail: lairuixun@163.com

YU Xin (余欣), YANG Ming (杨明), WANG Ming (王明)

Yellow River Institute of Hydraulic Research, Zhengzhou 450003, China

(Received August 25, 2012, Revised November 2, 2012)

In this paper, both state variables and parameters of one-dimensional open channel model are estimated using a framework of the Ensemble Kalman Filter (EnKF). Compared with observation, the predicted accuracy of water level and discharge are improved while the parameters of the model are identified simultaneously. With the principles of the EnKF, a state-space description of the Saint-Venant equation is constructed by perturbing the measurements with Gaussian error distribution. At the same time, the roughness, one of the key parameters in one-dimensional open channel, is also considered as a state variable to identify its value dynamically. The updated state variables and the parameters are then used as the initial values of the next time step to continue the assimilation process. The usefulness and the capability of the dual EnKF are demonstrated in the lower Yellow River during the water-sediment regulation in 2009. In the optimization process, the errors between the prediction and the observation are analyzed, and the rationale of inverse roughness is discussed. It is believed that (1) the flexible approach of the dual EnKF can improve the accuracy of predicting water level and discharge, (2) it provides a probabilistic way to identify the model error which is feasible to implement but hard to handle in other filter systems, and (3) it is practicable for river engineering and management.

Ensemble Kalman Filter (EnKF), lower Yellow River, water-sediment regulation, inverse problem

Introduction

The state variables and the parameters are two essential factors to determine the accuracy of a hydrodynamic model. The state variables which are in a state of flux are propagated over time controlled by the numerical model while the parameters describe the dynamic characteristics of a river, such as bed roughness. The numerical models, however, with some limitations or hypotheses or both, are imperfect to describe all aspects of a natural river as its complexity can be controlled not only with hydraulics but also with geological theories or sedimentation[1]. In recent decades, the improvements in prediction accuracy of water levels and discharge, for example, can be classified in three ways: one is to estimate the parameters using sets of algorithms[2-6], another is to construct an optimal scheme by coupling the one- and two-dimensional models or by coupling the hydrodynamic model with the hydrology model[7,8], and the third is to train the optimal algorithms using historical data sets[9].There are many sources of uncertainties associated with these applications which are essential to affect the model output. Although parts of these studies contain the effects of parameter uncertainties on generating accurate forecasts, other sources of errors, such as input, output, model structure, and observations, have not yet been considered. In other words,the uncertainties of the model structure are added to parameters’ error.

Fig.1 The flowchart of dual state-parameter estimation using EnKF

The Kalman filter, which, in a linear system, takes into account all sources of the uncertainties from both the model structure and observations, is one of the widely-used sequential assimilation techniques. Although much research using the Kalman filter has obtained optimal roughness parameters[10], their drawbacks can be listed as follows: (1) some studies considered the parameters to be of static values or the assimilation is applied only to state variables, (2) the Kalman filter is limited to the linear system, and the filter process may result in computational instability as the hydrodynamic equations is a nonlinear system. Furthermore, the accuracy of the predicted value would decrease if either the observations are scarce or the measure error is large, (3) it is difficult to define the hydrodynamic model error. To deal with the drawbacks of the Kalman filter, a Monte Carlo based Ensemble Kalman Filter (EnKF) was introduced. Compared to the traditional Kalman filter, one of the advantages of the EnKF is that it need not define the model error directly.

The series of Kalman filters, including the EnKF, was originally developed for dynamic systems while in this paper a strategy of dual state-parameter assimilation is introduced for hydrodynamic systems. With observation, the updated states and parameters are considered as the initial values of the next time step to continue the assimilation process. The usefulness and the capability of the dual ensemble Kalman filter are demonstrated in studying the lower Yellow River during the water-sediment regulation in 2009.

1. Data assimilation and equations

1.1 Ensemble Kalman Filter (EnKF)

Apart from the traditional Kalman filter, the EnKF represents the model uncertainty by perturbing the forcing data and parameters:

where Uit+1is the ith ensemble member of forcing data at time t+1 and θit+1is the ith ensemble of para- meter at time t+1. The model error is displayed by ξit+1and N(0,Ri+1) refers to the normal distributiont with a mean of 0 and a covariance of R.

Driven by the forcing data, the fluxes of both state variables and parameters can be described by the theory of nonlinear stochastic-dynamic processes. The model forecast of each ensemble member is as follows

where X represents the state variable, Xi-t+1is the ith forecast ensemble at present time t+1, Xti+is the ith optimal ensemble at the prior time t and ωit+1is the error of thi ensemble at time +1t.

Similar to the forcing data, the ensemble of the observation is generated by perturbing the measurements

where Yit+1is the ith trajectory of the observation by adding the error ε with the covariance of S. Generally, the noise of the observation can be considered as a normal distribution.

Then the prior error covariance of the model can be approximately equal to

where Kt+1is the Kalman gain matrix and H is the transitional operator weighting the value between observation and model.

1.2 State-space description of equations

The Saint-Venant equation considering the lateral inflow can be written as

where B is the width of the cross section, Q total volume discharge, Z the water level, lq the side flow discharge per unit channel length, R the hydraulic radius and C the Chézy coefficient which can be calculated as

where n is the roughness. Using the Preissmann scheme, the state-space description of the Saint-Venant equation can be constructed as follows:

where ωtQZis the model error, νt+1QZthe observation error of water level and discharge, Ht+1QZthe observation operator and its value is equal to 1 once there is observational data, UtQZthe control matrix and Φt+1/tQZthe transport matrix.

Table 1 Error distributions of observations and initial values of parameters

The purpose of the roughness inverse problem is not only for accurate prediction, but also one of the ways to deepen understanding of the system characteristics[11]. For the roughness inverse problem, observation is needed. However, the bed roughness can not be measured directly. This paper constructs the observation of roughness using the difference of the water level between two cross sections[12]. The roughness parameter which is only included in the momentum equation and its finite difference equation can be yielded from Eqs.(10) and (11):

where tΔ is the time step and xΔ is the distance between cross sections. Then the difference of the water level can be written as

Using the Taylor series expansion and neglecting the second-order term, the observation for roughness in Eq.(14) can be rewritten as

So the state-space form of the roughness can be given as

where ρ refers to the initial difference of roughness from the previous time to the present time.

Fig.2 Plot of the lower Yellow River with tributaries and hydrological stations

Fig.3 Root Mean Square Errors (RMSEs) with and without data assimilation (DA in the legend refers to data assimilation)

2. Assimilation process of dual state-parameter

Traditionally, in the optimal process of stochastic systems, the parameters are considered to be time invariant. Within the framework of the dual state-parameter system, both the state variables and parameters are calculated or predicted simultaneously. One method for such combined state variables and parameters is provided by joint estimation in that both stateand parameter vectors are concatenated to form a single vector. The drawback of this method is that it will result in unstable estimation when the number of the states or parameters increases, especially in the nonlinear system[13]. Another way to cope with these limitations is to estimate the state variable and the parameter separately, in which two interactive filters are designed either to estimate the optimal states or to solve the inverse problem of parameters[14].

Figure 1 illustrates the process of dual state-parameter estimation using the EnKF with its application in one-dimensional hydrodynamic model. The total number of the ensemble is estimated first at time t, and then the observations such as the water level, discharge and roughness parameters are perturbed. The normal distribution of observations and the rational range of parameters are listed in Table 1. According to the regulation of river measurement, the measurement uncertainty of water level is less than three centimeters. The stochastic error of the discharge is restricted to normal distribution with the confidence level of 95%. So the error of discharge accounts for about 5% of its total volume. The initial value of roughness is from 0.03 in the input and generally decreases to 0.01 in the output.

Using Eq.(3), the ensemble of the water level or discharge is propagated and each member represents one realization of the model replicates. Then, the prior covariance of the model is calculated using Eq.(5) and, once the observation is available, the optimal states can be yielded with Eq.(7). Next, using Eq.(15), such optimal water level and discharge are adopted to calculate observed value of roughness. After that, inverse roughness, which is considered to be the optimal value at this time, is evolved using the Kalman filter algorithm. Finally, the updated states and roughness values are taken as the initial values in the next time step to continue the process of assimilation.

Fig.4 Simulated discharge with and without assimilation against observations at the six hydrological stations from Huayuankou to Luokou

3. Results and analysis

One of the goals of regulating water and sediment processes is to flush out the sediment and to enhance flow conveyance capacity of the main channel.Figure 2 shows the study area, which is located in the lower Yellow River, from Xiaolangdi to Lijin crossing Henan and Shandong provinces, about 650 km, with eight hydrological stations along the river: Xiaolangdi, Huayuankou, Jiahetan, Gaocun, Sunkou, Aishan, Luokou and Lijin, while Xiaolangdi and Lijin acting as the input and output boundaries respectively.

During the 7th regulation in 2009, the data assimilation scheme described above was adopted from 8 am of June 15 to 8 am of July 13, a total of 672 h. To make the simulated close to the true process of a natural flood, lateral inflow, irrigation, evaporation and infiltration was calculated in the numerical model. During the regulation, the total inflow of 51 m3/s from Qinhe and Yiluohe rivers, the outflow for irrigation of 388 m3/s, and the lose due to evaporation and infiltration of 120 m3/s were taken into account.

3.1 Error analysis for water level and discharge

Water level and discharge are two most important state variables in one-dimensional hydrodynamic model. The RMSE is used to measure the accuracy of predicted state variables. Figure 3(a) plots the RMSEs of discharge with and without data assimilation. Without data assimilation, the discharge RMSE decreased slightly from 195.6 in the Huayuankou Station to 168.1 in the Jiahetan Station and increased greatly to about 270 in the Lijin Station. However, with data assimilation, the RMSE dropped obviously from 106.5 in the Huayuankou Station to 76.4 in the Lijin Station.

Figure 3(b) shows the result of RMSEs of water levels with and without data assimilation. Without the assimilation method, the RMSE decreased gradually from 0.44 in the Huayuankou Station to 0.24 in the Sunkou station, and it remained almost stable at about 0.4 from the Aishan Station to the Luokou Station. Sources of error were attributed primarily to the unchanged roughness value and, secondly, to the absence of sediment calculation. However, with data assimilation method, the RMSE dropped dramatically from 0.04 in the Huayuankou station to 0.09 in the Lijin Station.

3.2 Analysis of predicted discharge

During the regulation of flow and sediment processes in 2009, the simulated peak flow of the Huayuankou station was 4 220 m3/s, almost equal to the bankfull discharge, and maintained for a short time. For this reason, flood was controlled in the main channel.

Figure 4 shows the measured discharge in comparison with the predicted result. Without data assimilation, the predicted processes of peak flow and its attenuation were correspondent with the observations. As for the simulated value at each hydrological station, the predicted discharge values at Huayuankou and Jiahetan was in accordance with the observations since the lateral inflow and outfow were taken into consideration. For the stations from Gaocun to Lijin, however, the predicted discharges were much higher than the observations, and those discrepancies can be ascribed to (1) the inconsistence of the irrigation lose calculated in the model with the real process and (2) the errors from numerical model.

It is proved that, with assimilation system accompanied by observed discharge, a prediction of high accuracy has been gained.

Fig.5 Results of sediment erosion and deposition after the 7th water-sediment regulation (the negative number means erosion while the positive one means deposition E and D means erosion and deposition, X-H means river from Xiaolangdi to Huayuankou, H-J means Huayuankou to Jiahetan, J-G means Jiahetan to Gaocun, G-S means Gaocun to Sunkou, S-A means Sunkou to Aishan, A-L means Aishan to Luokou and L-L means Luokou to Lijin)

3.3 Analysis of water level and inverse roughness

Since the sediment was not calculated in the current model, it is necessary to understand the results of erosion and deposition before analyzing the predicted water level. Figure 5 illustrates the results of erosion and deposition after the 7th water-sediment regulation. It shows that, along the lower Yellow River, except for a small amount of sediment deposition from Aisha to Luokou stations, large amount of erosion happened and the number went even up to 9.96 ×106t from Gaocun to Sunkou stations.

Figure 6 shows the predicted water level compared with the inverse roughness. Without data assimilation, the simulated water level was a little higher than the observations from Huayuankou to Aishan while it was lower at the Luokou station. One reason for the discrepancies was unreasonable initial value of roughness, and another was that there is no sediment calculation adopted to adjust the change of bedform. Compared with the results of sediment transport plotted in Fig.5, it is believed that, if the sediment calculation is taken into the model, there will be a better simulated result in which the predicted water level will decrease from the Huayunkou to Aishan stations in the course of eroding while increase at the Luokou Station in the course of deposition.

It is illustrated that, with assimilation system accompanied by observed water level, a more accurateprediction has been gained.

The inverse roughness, which determines the accuracy of the calculated water level and depends on a number of factors, was estimated using Kalman filter algorithm. Three conclusions about inverse roughness could be reached. Firstly, as the discharge increases, the inverse roughness drops simultaneously. This phenomenon, which is identical to the empirical relations between roughness and discharge, is obviously illustrated at the four stations including those at Huayuankou, Gaocun, Sunkou and Aishan. Secondly, without the assimilation, the roughness parameter is a unique value which is 0.025 from Xiaolangdi to Jiahetan, 0.015 from Gaocun to Sunkou and 0.01 from Aishan to Lijin. But the inverse roughness values, which varies for each cross section and changes with water level and discharge, are gained using assimilation method. Such inverse roughness values for each cross section at each time are beneficial to getting a better predicted result. Thirdly, the inverse roughness values drop gradually from Xiaolandi to Lijin, which is similar to the changes of the roughness values without data assimilation.

Fig.6 Roughness using inverse method at the six hydrological stations

Observations are essential to the EnKF. Unfortunately, the data obtained by current measurements always lags far behind the process of assimilation. Therefore, to achieve an ideal system of real time estimation, much research is still needed.

4. Conclusion

In this paper, a data assimilation framework using the dual EnKF for one-dimensional hydrodynamic model has been constructed, and its usefulness and applicability have been demonstrated in the prediction of water level and discharge in the downstream areas of the Yellow River. During the water-sediment regulation in 2009, the errors both with water level and discharge, and of the difference between observation and prediction, are analyzed. The rational value of the inverse roughness is discussed considering the whole process of flooding. It is believed that (1) the flexible approach of the dual EnKF can improve the accuracy of predicted water level and discharge, (2) it provides a probabilistic way to identify the model error which is feasible to implement but hard tohandle in other filter systems, (3) the inverse problem of roughness can be used in order to understand its activity in such a dynamic system in greater depth, and (4) the algorithm need not store all the past information in order to be practicable for river engineering and management.

[1] SCHUMM S. A. River variability and complexity[M]. Cambridge, UK: Cambridge University Press, 2005, 4-8.

[2] TANG Hong-wu, XIN Xiao-kang and DAI Wen-hong et al. Parameter identification for modeling river network using a genetic algorithm[J]. Journal of Hydrodynamics, 2010, 22(2): 246-253.

[3] De DONCKER L., TROCH P. and VERHOEVEN R. et al. Determination of the Manning roughness coefficient influenced by vegetation in the river Aa and Biebrza river[J]. Environmental Fluid Mechanics, 2009, 9(5): 549-567.

[4] BAO Wei-min, ZHANG Xiao-qin and QU Si-min. Dynamic correction of roughness in the hydrodynamic model[J]. Journal of Hydrodynamics, 2009, 21(2): 255-263.

[5] FORZIERI G., DEGETTO M. and RIGHETTI M. et al. Satellite multispectral data for improved floodplain roughness modeling[J]. Journal of Hydrology, 2011, 407(1-4): 41-57.

[6] ARICÒ C., NASELLO C. and TUCCIARELLI T. Using unsteady-state water level data to estimate channel roughness and discharge hydrograph[J]. Advances in Water Resources, 2009, 32(8): 1223-1240.

[7] KIM J., WARNOCK A. and IVANOV V. Y. et al. Coupled modeling of hydrologic and hydrodynamic processes including overland and channel flow[J]. Advances in Water Resources, 2012, 37(3): 104-126.

[8] JOWETT I. G., DUNCAN M. J. Effectiveness of 1D and 2D hydraulic models for instream habitat analysis in a braided river[J]. Ecological Engineering, 2012, 48(11): 92-100.

[9] XU S., HUANG W. Estimating extreme water levels with long-term data by GEV distribution at Wusong Station near Shanghai City in Yangtze Estuary[J]. Ocean Engineering, 2011, 38(2): 468-478.

[10] WU Xiao-ling, WANG Chuan-hai and CHEN Xi et al. Kalman filtering correction in real-time forecasting with hydrodynamic model[J]. Journal of Hydrodynamics, 2008, 20(3): 391-397.

[11] RAOL J. R., GIRIJA G. and SINGH J. Modelling and parameter estimation of dynamic systems[M]. London, UK: The Institution of Engineering and Technology, 2004, 1-2.

[12] LI Da-yong, DONG Zeng-chuan and LIU Ling et al. Real-time alternative updating model coupling 1-D unsteady flow computation with Kalman filter[J]. Journal of Hydraulic Engineering, 2007, 38(3): 330-336(in Chinese).

[13] MORADKHANI H., SOROOSHIAN S. and GUPTA H. V. et al. Dual state-parameter estimation of hydrological models using ensemble Kalman filter[J]. Advances in Water Resources, 2005, 28(2): 135-147.

[14] LU H., YU Z. and ZHU Y. et al. Dual state-parameter estimation of root zone soil moisture by optimal parameter estimation and extended Kalman filter data assimilation[J]. Advances in Water Resources, 2011, 34(3): 395-406.

10.1016/S1001-6058(11)60397-2

* Project supported by the National Basic Research and Development Program of China (973 Program, Grant No. 2011CB403306), the Ministry of Water Resources’ Special Funds for Scientific Research on Public Causes (Grant No. 200901023), and the Central Scientific Institutes Foundation for Public Service (Grant No. HKY-JBYW-2012-5).

Biography: LAI Rui-xun (1981-), Male, Ph. D. Candidate Senior Engineer