Ensemble Kalman filter

The Kalman filter is based on a stochastic formulation of the model equation the discrete state-transition model for the MIKE 1D dynamic system and measurement equation the measurement operator mapping model state to observation space. By including possible errors on the forcing terms the system gets the following stochastic formulation:

\[ { x}_{k~+~1}~=~\Phi({ x}_{k},{ u}_{k}+\varepsilon_{k}) \]

Equation: Stochastic state transition equation with model errors in forcing terms.

where

e_kis a vector containing the model error of the forcing terms.

The stochastic representation of the measurements is given by:

\[ z_k = C_k x_k + n_k \]

Equation: Measurement equation relating observed data to model state through observational operator and measurement error.

where

h_kis a random measurement error vector with zero mean and covariance matrix R_k.

The uncertainty of the model forecast is described by the covariance matrix \(P_k^f\). The Kalman gain to be used in the update equation the kalman filter analysis step combining model forecast with observations via gain is given by

\[ G_k = P_k C_k^T (C_k P_k C_k^T + R_k)^{-1} \]

Equation: Kalman gain matrix balancing model and measurement uncertainties in state updates.

which depends on the model and measurement covariances. The model covariance expresses how the innovation vector should be distributed on the state vector, and the relation between the model and measurement covariances expresses how much weight should be put on the measurement. The smaller measurement variance as compared to the model variance, the more weight is put on the measurement.

A major shortcoming of the above formulation is the assumption that the covariance matrix \(P_k^f\) is known a priori. For non-linear and high-dimensional systems such as MIKE 1D the evaluation of the covariance matrix requires huge computational costs and storage requirements, and hence makes the filtering infeasible for real-time applications.

In the EnKF, the covariance matrix is represented by an ensemble of possible state vectors. Each of these vectors is propagated according to the dynamical system subjected to model errors cf. Eq. the rainfall-driven erosion and detachment rate proportional to rainfall intensity, and the resulting ensemble then provides estimates of the forecast state vector and covariance matrix. In the measurement update, the Kalman gain matrix is applied for each of the forecast state vectors. To account for measurement errors, the measurements are represented by an ensemble of possible measurements. The resulting updated sample provides estimates of the updated state vector and covariance matrix.

The EnKF algorithm is summarised below. It is asumed that measurement errors are uncorrelated i.e. the covariance matrix R_k is a diagonal matrix

\[ R_{k} = \text{diag}\left(\sigma_{1}^{2}, \ldots, \sigma_{p}^{2}\right) \]

Equation: Diagonal measurement error covariance matrix for the Ensemble Kalman Filter.

and hence the sequential updating algorithm in Eq. the sequential measurement update for data assimilation processing independent can be applied.

Let there be given a set of M state vectors \(x_{j,k}^{f}\), where j=1,...,M and k is the time step index.

1. Each member of the ensemble is propagated forward in time according to the stochastic dynamical system given by Eq. the rainfall-driven erosion and detachment rate proportional to rainfall intensity. For each member of the ensemble the model error is drawn from a Gaussian distribution with zero mean and known covariance.

2. The forecast of the state vector is then estimated as the mean of this ensemble:

\[ {{\bar{x}}}_{k}^{f}~=~{\frac{1}{M}}{\sum_{j=1}^{M}}{x}_{j~k}^{f} \]

Equation: Ensemble mean forecast state estimate from multiple ensemble members.

1. An estimate of the covariance matrix is found according to:

(8.12a)         \(P_k^f = \frac{1}{M - 1} S_k^f (S_k^f)^T\),

(8.12b)      \(S_{j,k}^f = x_{j,k}^f - \bar{x}_k^f\),

1. Each of the ensemble state vectors are updated sequentially using Eq. the sequential measurement update for data assimilation processing independent. The gain vector for measurement No. i is given by:

\[ g_k^j = \frac{1}{c_k^j P_{k,i-1}^a (c_k^j)^T + \sigma_i^2} P_{k,i-1}^a (c_k^j)^T \]

Equation: Ensemble Kalman filter gain vector for sequential measurement update.

• In the above sub-/super-script i is a running index indicating elements corresponding to the i´th update using the i´th measurement \(\overline{x}_{k}^{f}\) . The term \(c_k^j\) is the i´th row of the information matrix C_k described earlier.

• The updated state vector and covariance matrix are estimated as:

\[ {{\bar{x}}}_{k}^{a}~=~\frac{1}{M}{\sum_{j~=~1}^{M}}{x}_{j~k}^{a} \quad \text{and} \quad P_k^a = \frac{1}{M - 1} S_k^a (S_k^a)^T \]

Equation: Ensemble mean state and covariance matrix for ensemble Kalman filter analysis.

In the implementation of the EnKF the covariance matrix \(P_{k}\) is never calculated. All calculations are based on \(s_{k}\).

When no data are available for updating, only steps 1-3 are carried out. This is the approach taken in MIKE 1D DA when performing uncertainty assessment and making a stochastic or ensemble forecast.

Model errors described by auto-regressive processes

The model errors are defined by auto-regressive processes of first order. Thus, for each uncertain boundary condition l = 1,2,...,L the model errors are described by

\[ \varepsilon_{i,k+1}^l = \phi_l \varepsilon_{i,k}^l + \delta_{i,k+1}^l \quad \text{and} \quad i = 1, 2, \ldots, M \]

Equation: First-order auto-regressive process for model error at uncertain boundary conditions.

where \(\epsilon_{i, k + 1}^{l}\) is randomly drawn from a Gaussian distribution with zero mean and user defined standard deviation s_d,j. Furthermore, all \(\epsilon_{i, k + 1}^{l}\) are assumed to be uncorrelated. It should be emphasised that to ensure that the errors do not grow without limit the absolute value of the regression coefficients f must be less than unity.

By using an auto-regressive description the model error has some memory that is propageted as part of the model forecast. In this case the model state is augmented with the model errors that are updated as part of the Kalman filter update scheme. In forecast mode, the model errors at the time of forecast are phased out according to the exponential decay defined by the auto-regressive process.