The stochastic EnKF

The stochastic EnKF was first introduced in the geophysics literature as an alternative to the traditional Kalman filter. Geophysical models, such as reservoir models or models of the atmosphere, preclude straightforward implementation of the traditional Kalman filter for two main reasons. Firstly, they are typically of such high dimensions that explicit storage and computation of the full n⇥n Kalman filter covariance matrices become problematic. Secondly, they usually involve a forward model p(xt|xt 1) which is non-linear and/or non-Gaussian. By exploiting the ensemble representations of the forecast and filtering distributions, the EnKF avoids explicit computation of the full n⇥ncovariance matrices and is able to cope with mild features of non-linearity and non-Gaussianity in the forward model.

Like the traditional Kalman filter, each iteration of the EnKF involves a fore-cast step and an update step. The update step can be viewed as an ensemble-based version of the update step of the traditional Kalman filter. Iteration number t starts with the assumption that a forecast ensemble, {x⁽¹⁾_t , . . . , x^(M)_t }, with in-dependent realisations from the forecast distribution p(xt|y1:t 1) is available. In reality, this assumption holds only approximately. The update step is then per-formed by first (implicitly) approximating the forecast model p(xt|y1:t 1) with a Gaussian distributionN(xt; ˆµt,Pˆt) with mean vector ˆµtequal to the sample mean of the forecast ensemble,

ˆ µt= 1

M XM

i=1

x⁽ⁱ⁾_t ,

and covariance matrix ˆPtequal to the sample covariance,

Pˆt=XtX_t^>, (18)

where

Xt= 1

pM 1

⇣

x⁽¹⁾_t µˆt, . . . , x^(M)_t µˆt

⌘

is the so-called forecast ensemble-anomaly matrix. Thereafter, each forecast sam-ple is linearly shifted according to

e

x⁽ⁱ⁾_t =x⁽ⁱ⁾_t + ˆKt

(

yt Htx⁽ⁱ⁾_t +!_t⁽ⁱ⁾

)

, (19) where!_t⁽ⁱ⁾⇠N(yt; 0, Rt),!⁽¹⁾_t , . . . , !_t^(M) are all independent, and

Kˆt=Xt(HtXt)^> HtXt(HtXt)^>+Rt

1 (20)

is the Kalman gain matrix in Eq. (17), onlyPtis replaced with the sample co-variance matrix ˆPt in Eq. (18). The justification of the update in Eq. (19) is that, under the assumption that the Gaussian approximationN(xt; ˆµt,Pˆt) is the correct forecast model, i.e. that x⁽¹⁾_t , . . . , x^(M)_t are independent samples from N(xt; ˆµt,Pˆt), the update yields independent samplesex⁽¹⁾_t , . . . ,ex^(M)_t from the cor-responding correct filtering distribution, which then is a Gaussian N(xt; ˆ˜µt,Pˆ˜t) with mean ˆ˜µtand covariance ˆ˜Ptavailable from the Kalman filter equations in Eqs.

(15) and (16),

µˆ˜t= ˆµt+ ˆKt(ˆµt Htµˆt) (21) and

Pˆ˜t= (In KˆtHt) ˆPt. (22) In geophysical applications, the ensemble size M is typically much smaller than the state dimensionndue to computer demanding forward modelsp(xt|xt 1).

The observation dimension,m, is also smaller thann, though considerably bigger than M. From Eqs. (19) and (20), one should note that explicit storage and computation of the fulln⇥ncovariance matrix ˆPtcan be avoided and that the largest matrices that need to be maintained are of sizen⇥M,n⇥mandm⇥m.

This allows for an efficient numerical implementation compared to the traditional Kalman filter and makes the EnKF computationally feasible also in large-scale applications.

Having generated the filtering ensemble {ex⁽¹⁾_t , . . . ,ex^(M)_t } according to Eq.

(19), the forecast step is performed by generating a new forecast ensemble, {x⁽¹⁾_t+1, . . . , x^(M)_t+1}, by simulating from the Markov forward model,

x⁽ⁱ⁾_t+1|ex⁽ⁱ⁾_t ⇠p

(

xt+1|ex⁽ⁱ⁾_t

)

, (23) independently for i= 1, . . . , M. In contrast to the update step, which relies on Gaussian approximations, the forecast step is exact in the sense that under the as-sumption thatex⁽¹⁾_t , . . . ,ex^(M_t ⁾are exact and independent samples from the true fil-tering distributionp(xt|y1:t), Eq. (23) returns forecast samplesx⁽¹⁾_t+1, . . . , x^(M)_t+1 that are exact and independent samples from the true forecast distributionp(xt+1|y1:t).

If the underlying state-space model really is linear-Gaussian, the EnKF is correct in the limit of an infinite ensemble size. The solution it then provides corresponds to that of the Kalman filter. In the more general case of a non-linear, non-Gaussian model, the filter provides a biased solution. However, if there are non-Gaussian features present in the forecast ensemble, it is to some extent possible for the filtering ensemble to capture some of these properties, since it is obtained by linearly shifting the forecast samples.

4.3 EnSRFs

EnSRFs perform a deterministic version of the update in Eq. (19). As the traditional, stochastic EnKF, EnSRFs start by replacing the forecast distribution p(xt|y1:t 1) with a Gaussian approximationN(xt; ˆµt,Pˆt) from which a correspond-ing Gaussian approximationN(xt; ˆ˜µt,Pˆ˜t) to the filtering distributionp(xt|y1:t) fol-lows from Bayes’ rule. While the stochastic EnKF updates the forecast ensemble so that the sample mean and sample covariance of the updated ensemble equal ˆ˜µt

and ˆ˜Ptin expectation, or in the limit of an infinite ensemble size, EnSRFs deter-ministically update the ensemble so that the sample mean and sample covariance equal ˆ˜µtand ˆ˜Ptexactly.

EnSRFs involve computations with matrix square roots. Similarly to the ensemble-anomaly matrixXtof the forecast ensemble, an ensemble-anomaly ma-trixXet of the filtering ensemble can be defined as

Xet= 1

pM 1

⇣ ex⁽¹⁾_t , . . . , ex^(M_t ⁾⌘

, (24)

where

e

x⁽ⁱ⁾_t =ex⁽ⁱ⁾_t 1 M

XM j=1

e

x^(j)_t . (25)

The ensemble-anomaly matricesXtandXetare matrix square roots of the sample covariance matrices of the forecast and filtering ensembles, respectively. The strategy of EnSRFs is to compute an ensemble-anomaly matrixXet by requiring that the sample covariance matrix of the filtering ensemble is exactly equal to ˆ˜Pt. Specifically, this requirement entails that

Pˆ˜t=XetXe_t^>. (26) Using Eqs. (18) and (26), and definingDt=HtXtX_t^>H_t^>+Rt, Eq. (22) can be written as

XetXe_t^>= In XtX_t^>H_t^>D_t¹Ht XtX_t^>. Rearranging terms on the right-hand-side, we obtain

XetXe_t^>=Xt IM X_t^>H_t^>D_t¹HtXt X_t^>.

Thereby, we see that a posterior ensemble-anomaly matrixXet with the desired properties can be obtained as

Xet=XtWtU

whereWt2R^M⇥M is a matrix square root of IM X_t^>H_t^>D_t¹HtXt , i.e.

WtW_t^>= IM X_t^>H_t^>D_t¹HtXt ,

and U is an M ⇥M orthonormal matrix, i.e. U U^> = U^>U = IM. The pos-terior ensemble members ex⁽¹⁾_t , . . . ,ex^(M)_t can thereafter be obtained by inserting

1 M

PM

j=1ex^(j)_t = ˆ˜µt in Eq. (25). Since the matrices Wt and U are not unique, except in the univariate case, a variety of EnSRF schemes can be formulated. As such, several EnSRFs have been proposed in the literature, e.g. Anderson (2001), Bishop et al. (2001), Whitaker and Hamill (2002), and Evensen (2004). J. L.

Anderson (2001)

(a) (b) (c)

Figure 2: Demonstration of covariance localisation. (a) True covariance matrix,P. (b) Sample covariance matrix, ˆP, obtained from M = 20 random samples. (c) Regularised covariance matrix obtained from the Schur product of⇢and ˆP with⇢given by Eq.

(29) usingL= 10.

In document Ensemble updating for a state-space model with categorical variables (sider 22-26)