Abstract

In this paper, a Moving Horizon Estimator (MHE) and an Unscented Kalman Filter (UKF) are applied and compared for state estimation in flood forecasting. The investigations are based on a conceptual rainfall-runoff model proposed by Lorent/Gevers for streamflow forecasting. Data for the investigations was collected from the region Trusetal in Germany. Streamflow prediction, especially for watersheds with fast response to intense rain, require the knowledge of the current state of the system (e.g., soil moisture content). Firstly, a Moving Horizon Estimator (MHE) was applied for the state estimation, due to our good experience with it in other applications, its ability to deal with non-Gaussian disturbances and the fact that the hydrologic model is nonlinear, and its states satisfy equality and inequality constraints. Due to computational intensity of the MHE, an UKF was also implemented for comparison. Even though theory and most literature conclude the superiority of MHE to UKF, in this application example the results show that the UKF and the MHE produce almost similar results with UKF slightly better, which might be due to several reasons such as problems with the initialization of the hessian matrix, choice of prediction horizon and existence of local optima in MHE. Therefore, comprehensive investigations were performed in this respect.

1. Introduction

During the last decade several flood events with partially disastrous extent occurred in Central Europa leading to losses of human life as well as to serious economic damages ¹. By now most of the large rivers are therefore equipped with operation warning systems. However, more frequent heavy thunderstorms with high rainfall intensities have also led to an increasing occurrence probability of flash floods in smaller streams. Consequently, there is a need to provide comprehensive forecast modules for such events to provide in time warning to endangered communities ^{2, 3, 4}. Now, a pilot flash flood warning system for one municipality in the federal state Thuringia (Germany) is under development.

Despite the progress in high resolution data sources for hydrologic models (e. g. radar-based precipitation measurements or satellite-based soil characteristics and land use data), which force the application of physical based distributed model approaches, simple (lumped) conceptual models or even statistical “black-box” models still outperform the more involved model concepts regarding online forecasts of river flow ^{3, 4}.

A successful application of rainfall-runoff models for flood forecasting purposes depends on the availability of the current model state. Fundamentally different sources of uncertainty, such as (a) error in input data to the model; (b) error in measurement of output from nature;(c) error in model structure; (d) non-optimal values of model parameters, exist in rainfall- runoff modelling. Hence an adjustment with respect to random or systematic errors in the input variables (precipitation) is of practical relevance, as uncertainties in the input variables are to be considered as substantial source of errors ⁵. The Kalman filtering technique and the state space theory has proven to be powerful mathematical tools for treating various uncertainties in mathematical modelling ⁶. Therefore, the (extended) Kalman filter has been the most common technique for state estimation of conceptual rainfall-runoff models ^{7, 8, 9, 10}. Recently, the ensemble Kalman filter has gained attraction for data assimilation of hydrologic models because of its improved robustness compared to other filtering techniques ¹¹ and its ability to treat nonlinear distributed models without the need of an explicit model linearization ^{12, 13}. Most physical systems, exhibit besides nonlinear dynamics also states subject to hard constraints, such as nonnegative soil moisture, concentrations, or pressures. Therefore, we propose the application of MHE techniques, which can also take constraints on the noise and the state space, as well as asymmetric probability distributions, into account. Another state estimation approach, which is from control theoretical point of view like moving horizon state estimation, is termed variational data assimilation, and is discussed in ¹⁴.

The flash flood warning system under development is based on a conceptual rainfall-runoff model, which is used for approximately 30 years for operational flow forecasting in Thuringia ¹⁵. This paper presents a comprehensive comparison of UKF with MHE in this hydrological state estimation application to provide a reliable starting point for the forecast procedure.

The paper is organized as follows. First the application area is introduced briefly. The rainfall-runoff model is presented in section III together with results for model calibration and validation. Section IV describes and give comparison results of the application of UKF and MHE algorithms for state estimation in the rainfall-runoff model. In Kalman filtering the process noise is the uncertainty in the model description assumed to originate both from uncertainty in the precipitation input and from parameter uncertainty. The moving horizon state estimator is extended to include uncertainty in the inputs. Results of including input uncertainties into the MHE state estimation routine are given in section V.

2. Materials and Methods

The catchment area of the Truse stream is situated in the low mountain range of the Thuringia Forest and has an area of approximately 28 km² (see Figure 1). It is covered mainly by coniferous woodland (60 %). The average annual precipitation ranges from 1288 mm at the highest spot (Großer Inselsberg) to 928 mm in the lower area of the catchment. The flow of the Truse stream is recorded since 1961. The mean annual outflow is 0.55 m³/s and the highest flow peak has been observed at 1981/08/14 with 14.6 m³/s. The time span for the change-over from normal flow conditions to the peak flow is very short (30 . . . 60 min).

Most (lumped) conceptual catchment area models are based on one or more (nonlinear) storages, which account for the basic hydrologic processes (soil moisture content, vegetation interception) and transform the precipitated water into the water delivered to the river. The flow along the river is described based on simplified routing approaches (e. g. kinematic wave or Muskingum-Cunge).

A Lorent-Gevers type storage in combination with two parallel dead-time/delay elements for the routing of the surface and the base flow, respectively is applied to the considered catchment area (see Figure 2) ¹⁶. This model is used successfully since 1980 for operational flow forecasting of the Werra river, which is the receiving river of the Truse stream.

The inputs of the model are area averaged values of the gross precipitation and the potential evapotranspiration which are derived from different numerical forecast products of the German Weather Service (DWD) and data of the land use/vegetation cover. As the relevant forcing input data are provided with a time resolution of, the time step for the model was chosen accordingly to

The rainfall lowered by the directly evaporating part is decomposed into the effective rainfall which constitutes the direct surface runoff, and the inflow to the storage, which accounts for the long-term behavior of the catchment area:

The water storage of the catchment area is described by a linear balance equation:

The infiltration depends linearly on the current storage content:

and forms the base flow from the catchment area. The evapotranspiration terms and are derived from the potential evapotranspiration and additional requirements in order to fulfill the overall water balance of the model:

The inflow to the water storage is determined by the following relation

which was proposed in ¹⁶. This relation ensures an increase in the direct surface flow with the amount of stored water and precipitation which corresponds to the natural behavior of catchment area. accounts for the limited ability of the soil to hold water and is a specific runoff parameter ().

The effective rainfall is split up in to two parallel delay elements by the parameter to distinguish a fast and slow response of the surface runoff :

The same applies to the infiltration from the storage, which forms the base runoff . Based on a cross correlation analysis for the catchment outflow and the area averaged precipitation a dead time of was selected for all flow routing elements.

The knowledge of the current state of the catchment area model is a prerequisite for the forecast of the catchment outflow. As the state variables of the considered conceptual rainfall-runoff model have no clear physical meaning, it is impossible to provide measurements for these states. Therefore, a state estimation procedure must be applied to gain the necessary state information. However, this is also required for mechanistic (distributed) catchment area models because the effort to retrieve the full state information is not feasible.

The EKF remains the most widely used state estimation algorithm for nonlinear systems. However, the extended Kalman filter, or EKF, is at best an ad hoc solution to a difficult problem, and hence there exist many barriers to the practical implementation of EKFs (see, e.g., ¹⁷). Some of these problems include the inability to accurately incorporate physical state constraints and poor use of the nonlinear model. Therefore, to overcome these problems we apply and compare the UKF and the MHE algorithms in this paper.

The Due to the flaws of the EKF (EKF is difficult to tune, the Jacobian can be hard to derive, and it can only handle limited amount of nonlinearity), we apply the UKF in this paper, which was first proposed by Julier et al ^{18, 19}. The UKF differs from the EKF in the manner in which the Gaussian random variables are represented for propagating through the system dynamics. In the EKF the state distribution, which is approximated by the Gaussian random variables is propagated analytically through the first order linearization of the nonlinear system. This procedure can cause errors in the posterior mean and covariance of the transformed Gaussian random variables leading to sub optimal performance of the filter. In the UKF a deterministic sampling approach is utilized. As with the EKF the state distribution is approximated by the Gaussian random variables, but this time is now represented by a minimal set of carefully chosen points, which completely capture the posterior mean and covariance accurately to the second order (Taylor series expansion) of any nonlinearity ¹⁸, whereas the EKF only achieves the first order accuracy. No explicit Jacobian calculations are necessary for the UKF which is very import in many hydrologic models.

The basic framework for the EKF involves estimation of the state of a discrete-time nonlinear dynamic system,

(1)

(2)

where represent the unobserved state of the system and is the only observed signal. The inputs and the process noise drives the dynamic system, and the observation noise is given by Note that we are not assuming additivity of the noise sources. The system dynamic model and are assumed known.

The UKF progresses in 4 stages, namely initialization, unscented transformation, prediction, and correction.

(1) Select the initial estimation For systems with deterministic state equation the diagonal elements of covariance matrix can be interpreted as the uncertainty to the selected initial values .

(2) Compute the sigma points of dimension from the given mean and covariance of the filter density :

whereby correspond to the i^th row of the root matrix. The parameter determines the width of the distribution of the sigma points about the mean

(3) The sigma points are subjected to the dynamics of the system. In other words, and are computed for all sigma points Their mean and covariance can be estimated from the resulting sample for a-priori-density and likelihood.

(3)

(4)

(5)

(6)

where

and The weights are defined by and For better comparability and to avoid non-local effects, the estimation should be independent of dimension of the state space. This can be achieved by setting

(4) The predicted moments are updated with the data by the Kalman correction equations.

(7)

(8)

with and

(9)

The steps (2) to (4) are to be continued, until either all available data points are utilized in the estimation, or the estimations has converged. The detailed procedure of the UKF can be found in ¹⁸.

The Together with the development of model predictive controllers the application of optimization-based methods for state estimation has gained attention during the last years. Compared to other techniques optimization-based methods, which are summarized as moving horizon state estimation ^{20, 21, 22}, have the advantage of a direct utilization of the nonlinear process model. This method can be adopted easily to problem specific needs, like e. g. the additional estimation of model parameters or exogenous inputs and can handle constraints on the estimation parameters in a straightforward manner (e, g. non-negativity of the storage content of the rainfall-runoff model). A drawback of this method is the significant increased computing time, which has no negative implications for the considered field of application because of low real-time requirements.

The moving horizon state estimation at the point in time is stated as follows:

(10a)

subject to

(10b)

(10c)

(10d)

(10e)

The process equations (10b) comprise the rainfall-runoff model described in section 3 with the state variables , the exogenous inputs and additional parameters to be estimated. The real measurements are denoted by It is assumed, that the measurement noise follows a zero-mean Gaussian process and the process noise is small compared to the uncertainties of the exogenous input variables and can be neglected. The matrices and as well as the estimation horizon are design parameters of the state estimation algorithm.

3. Results and Discussion

The remaining model parameters as well as the time constants of the four delay elements are determined using a least squares method. Based on experiences from a model of the Werra catchment area constraints are formulated to push the model parameters to reasonable ranges and to separate the time constants of the delay elements sufficiently. Figure 3 shows results for the calibration period from 2019/11/01 to 2021/10/31. The model has a good performance with a Nash-Sutcliffe efficiency of 0.74 and a bias of -8.46%. The optimized model parameters stay within the predefined bounds with exception of the time constant for the fast delay element of the base flow. Removing the corresponding constraint from the least squares problem the best fit of the measurements is reached for the minimum reasonable time constant of 1h. It seems, that the approach of Lorent/Gevers for the inflow to the storage directs especially in cases with high precipitation intensity too less water to the direct surface runoff, leading to an underestimation of the flow peak and a too slow decay of the outflow (see also Figure 4 for the flood event at day 103 during the validation period).

The model has been validated for the period from 2010/03/31 to 2010/08/11 resulting in a slightly decreased Nash-Sutcliffe of 0.68 and a bias of 4.0%. The rest difference between measured and simulated outflow after the second rainfall event (around day 50) may result from errors in the radar-based precipitation data, as there is almost no rain recorded during this event.

Figure 5 shows the effective rainfall in dependence of the total rainfall and the current storage content As can be seen, only for an almost saturated storage the total rainfall take effect immediately as surface runoff. Otherwise, the main portion of the total rainfall is directed to the storage. This behavior can be changed slightly by increasing the specific runoff parameter () leading to a reduced overall model performance. However, an accurate reproduction of the flash flood dynamics may require further investigations with respect to a more suitable model structure.

The performance of the two estimators (UKF, MHE) was investigated on the estimation of the states of the rainfall-runoff model for the validation period of the process model as discussed in section 3. For the MHE the investigations were carried out under the assumption of a perfect knowledge of the forcing input variables (total precipitation and potential evapotranspiration) whereas the UKF can take additive uncertainties of the inputs as part of the process model (see eq 1). The simulated model output was assumed to be the measured output. Under these (nominal) conditions the state of the rainfall-runoff model is reconstructed perfectly by the moving horizon state estimator (see. Figure 6-Figure 9 (nominal case)).

The estimation horizon has been selected to 240h, because of the long-time constants of the base flow delay elements. The situation changes if the simulated output is augmented with a zero-mean Gaussian noise process (see. Figure 6-Figure 9 (with disturbances)). While the estimated output is still in a very good coincidence with the measurements, the state of the storage as well as of the two base flow delay elements is not estimated accurately especially in periods with almost no rain. However, a jump of the state variables of the storage and the fast base flow delay element in one direction is accompanied by a jump of the state of slow base flow delay element in the other direction resulting in an overall smooth base outflow (see also Figure 9). The state of the two delay elements for the surface outflow is reconstructed with a good accuracy because it depends directly on the (known) rainfall. Note, that in periods with a high rainfall intensity (e. g. between day 40 and 60) the storage content can be reconstructed also by the state estimation procedure as the transformation of the total rainfall into the effective rainfall depends then strongly on the storage content.

The problems in adjusting the states to the right values by the state estimation procedure may result from the model structure with the four parallel pathways for the outflow composition. However, this may also indicate a general problem of deriving the correct state information for conceptual rainfall-runoff models (like e. g. tank or HBV model) solely from the catchment outflow as their model structures do not differ strongly.

Figure 10 shows the composition of the total outflow from the four flow fractions for the validation period of the catchment area model. The comparatively low outflow error for the first two rainfall events results from an estimated full storage (). Therefore, the rainfall during this time takes effect completely as surface runoff. This certainly conflicts with amount and decreasing behavior of the outflow in the previous period.

The results of the comparison between the MHE and the UKF state estimators are shown in Figure 11- Figure 17. It can be seen that both estimators produce similar results with the UKF slightly better in estimating all components of the flow. The errors of the MHE increases if the measured signal is back propagated instead of the simulated signal as can be seen when comparing Figure 11 and Figure 17.

The uncertainty with respect to precipitation (intensity, spatial and temporal variability) is one major source of errors in rainfall-runoff modelling. The moving horizon state estimator can be easily extended to determine uncertain input parameters. Assuming a correct identification of the temporal and spatial rainfall pattern by radar, a multiplicative approach can be used to eliminate magnitude errors. However, especially for small catchments in mountainous areas this might be not sufficient. Therefore, the following approach including an additive term for the last six hours of the estimation horizon is used to adjust the rainfall:

The parameters and are optimization variables, which can be subject to bounds derived from the expected uncertainty. The weighting matrices and in equation 10a are assigned as diagonal matrices, where the output error is penalized much stronger than the parameters for precipitation adjustment

Figure 18b shows the estimated outflow for the validation period of the catchment area model. Compared to the pure state estimation problem, the peak flow of the flash flood is reproduced better by estimating an additional rainfall of about over the last six hours prior to the rainfall event (see Figure 18a). In contrast to this offset, the scaling parameter is reduced to about 0.7 to force the fast declination of the catchment outflow after the peak. In this way the additional parameters also compensate for structural model errors.

4. Conclusions

A conceptual rainfall runoff model is proposed as a core element of a flash flood warning system for a small catchment area in Thuringia (Germany). The model is calibrated based on a time series for two years resulting in a sufficient overall model performance.

A good approximation of the current model state is a prerequisite for a successful online application of the rainfall-runoff model. Two state estimators are used and compared for this purpose. Given the same estimator tuning, model, and measurements as the UKF, simulation results show that MHE does not provide improved state estimation in our case as will be expected due to its ability to incorporates physical state constraints into an optimization and optimizes over a trajectory of states and measurements. It seems as if MHE fail to converge to the true state when the system model and measurement are such that multiple states constellation satisfy the steady-state measurement. The correct identification of the model states is even under idealized conditions (assumption of correct model and input variables) difficult with the MHE. Especially the state variables determining the long-term behavior (base flow) are hard to distinguish regarding their effect on the measurable output variable (catchment outflow). The introduction of additional parameters to adjust the precipitation as uncertain input variable led to a reduced approximation error for the measured output variable based on real world data.

Looking at the results of the MHE, issues of global versus local optimization, initialization of the Hessen Matrix and the selection of an arrival cost seem to have substantial impact on the behavior of MHE. These problems are very difficult to address but would be partly solved by implementing a global optimization strategy (problem: usually not realizable in real time), approximating the arrival cost with a uniform prior and making the estimation horizon reasonably long.

Even though according to the results of the MHE, the system may yield multiple optima corresponding to both physically realizable and unrealizable states if not constrained, UKF did not generate physically unrealizable states.

Statement of Competing Interests

There are no relevant financial or non-financial competing interests to the paper.

Acknowledgements

Thanks to the colleagues who worked in the data preparation and proof reading.

References