A Kernel Extension of the Ensemble Transform Kalman Filter

Abstract

Data assimilation methods are mainly based on the Bayesian formulation of the estimation problem. For cost and feasibility reasons, this formulation is usually approximated by Gaussian assumptions on the distribution of model variables, observations and errors. However, when these assumptions are not valid, this can lead to non-convergence or instability of data assimilation methods. The work presented here introduces the use of kernel methods in data assimilation to model uncertainties in the data in a more flexible way than with Gaussian assumptions. The use of kernel functions allows to describe non-linear relationships between variables. The aim is to extend the assimilation methods to problems where they are currently unefficient. The Ensemble Transform Kalman Filter (ETKF) formulation of the assimilation problem is reformulated using kernels and show the equivalence of the two formulations for the linear kernel. Numerical results on the toy model Lorenz 63 are provided for the linear and hyperbolic tangent kernels and compared to the results obtained by the ETKF showing the potential for improvement.

A Construction of \(P^a\) When K is Invertible

We consider the expression of $\mathbf {P^a}$ given by (17):

$$\begin{aligned} \mathbf {P^a} = \mathbf { K P^{\alpha } K} \end{aligned}$$

and substitute for $\mathbf {P^{\alpha }}$ its approximation by the hessian (18):

$$\begin{aligned} \mathbf {P^a} = \textbf{K}[(N-1) \mathbf {K + K \Pi _H \Pi _H K}]^{-1}\textbf{K} \end{aligned}$$

(26)

Since $\textbf{K}$ is invertible,

$$\begin{aligned} \Leftrightarrow \mathbf {P^a} = \frac{1}{N-1} (\textbf{K}^{-1} + \frac{1}{N-1} \mathbf {\Pi _H \Pi _H})^{-1} \end{aligned}$$

(27)

We the apply Woodbury identity:

$$\begin{aligned} \Leftrightarrow \mathbf {P^a} = \frac{1}{N-1} (\mathbf {K - K \Pi _H} ((N-1) \mathbf {I_{n+p} + \Pi _H K \Pi _H})^{-1} \mathbf {\Pi _H K}) \end{aligned}$$

(28)

$$\begin{aligned} \Leftrightarrow \mathbf {P^a} = \frac{1}{N-1} (\textbf{K} - \begin{bmatrix} \mathbf {0_n} &{} \mathbf {K_{HX}} \\ \mathbf {0_{pn}} &{} \mathbf {K_H} \end{bmatrix} \begin{bmatrix} (N-1) \mathbf {I_n} &{} \mathbf {0_{np}} \\ \mathbf {0_{pn}} &{} (N-1) \mathbf {I_p + K_H} \end{bmatrix}^{-1} \begin{bmatrix} \mathbf {0_n} &{} \mathbf {0_{np}} \\ \mathbf {K_HX}^{\top } &{} \mathbf {K_H} \end{bmatrix}) \end{aligned}$$

(29)

Let $\mathbf {K_H}$ decompose as $\mathbf {K_H = U_H \Sigma _H U_H}^{\top }$ with $\mathbf {\Sigma _H} = diag([\lambda _i]_{1 \le i \le p})$ and $[\lambda _i]_{1 \le i \le p}$ the eigenvalues of $\mathbf {K_H}$. We finally obtain:

$$\begin{aligned} \mathbf {P^a} = \frac{1}{N-1} \begin{bmatrix} \mathbf {K_X - K_{HX} U_H} diag(\frac{1}{(N-1)+ \lambda _i}) \mathbf {U_H}^{\top } \mathbf {K_{HX}}^{\top } &{} \mathbf {K_{HX} - K_{HX} U_H} diag(\frac{\lambda _i}{(N-1)+ \lambda _i}) \mathbf {U_H}^{\top } \\ \mathbf {K_{HX}}^{\top } - \mathbf {U_H} diag(\frac{\lambda _i}{(N-1)+ \lambda _i}) \mathbf {U_H}^{\top } \mathbf {K_{HX}}^{\top } &{} \mathbf {U_H} diag(\frac{\lambda _i (N-1)}{(N-1)+ \lambda _i}) \mathbf {U_H}^{\top } \end{bmatrix} \end{aligned}$$

(30)

This gives an explicit expression for $\mathbf {P^a}$ which does not require numerical inversion.