Approximating Kalman filtering in large dynamical systems Harri - PowerPoint PPT Presentation

Approximating Kalman filtering in large dynamical systems Harri Auvinen, Heikki Haario and Tuomo Kauranne Department of Applied Mathematics Lappeenranta University of Technology tuomo.kauranne@lut.fi 4 January 2006 – Typeset by Foil T EX –

Outline � Introduction to data assimilation � Kalman filter � Extended Kalman filter � 3D-Var � 4D-Var � Variational Kalman filter � Simulated assimilation results � Conclusions – Typeset by Foil T EX – 1

Introduction to data assimilation The formulation of the general data assimilation (state estimation) problem for discrete time steps t = 1 , 2 , 3 , ..., n contains an evolution or prediction equation and an observation equation: x ( t ) = M t ( x ( t − 1)) + E t (1) y ( t ) = K t ( x ( t )) + e t , (2) where M t is a model of the evolution operator, x t is the state of the process at time t and E t is a vector valued stochastic process. The second equation connects the state estimate x ( t ) to the measurements y ( t ) , with their associated observation errors e t , by way of an observation operator K t . – Typeset by Foil T EX – 2

Popular assimilation methods: � Kalman filter – Extended Kalman filter (EKF) – Fast Kalman filter (FKF) – Reduced Rank Kalman filter (RRKF) – Ensemble Kalman filter (EnKF) � 3D-Var � 4D-Var (used in operational weather forecasting) – Typeset by Foil T EX – 3

Basic linear Kalman filter The basic Kalman filter operates on a linear version of the general assimilation problem, equations (1) and (2): x ( t ) = M t x ( t − 1) + E t (3) y ( t ) = K t x ( t ) + e t , (4) where M t is the linear evolution operator and similarly K t is the linear observation operator. The Kalman filter operates sequentially in time. – Typeset by Foil T EX – 4

Kalman filter algorithm (1/2) Let x est ( t − 1) be an estimate of state x ( t − 1) and S est ( t − 1) be the corresponding error covariance matrix of the estimate. At time t the evolution operator is used to produce an a priori estimate x a ( t ) and its covariance S a ( t ) : x a ( t ) = M t x est ( t − 1) (5) M t S est ( t − 1) M T S a ( t ) = t + SE t , (6) where SE t is the covariance of the prediction error E t . – Typeset by Foil T EX – 5

Kalman filter algorithm (2/2) The next step is to combine x a ( t ) with the observations y ( t ) made at time t to construct an updated estimate of the state and its covariance: S a ( t ) K T t ( K t S a ( t ) K T t + Se t ) − 1 = (7) G t x est ( t ) = x a ( t ) + G t ( y ( t ) − K t x a ( t )) (8) S est ( t ) = S a ( t ) − G t K t S a ( t ) , (9) where G t is the Kalman gain matrix, which is functionally identical to the maximum a posteriori estimator. In a more general case, when the evolution model and/or the observation model is non-linear, the Extended Kalman Filter (EKF) is required. – Typeset by Foil T EX – 6

Extended Kalman filter algorithm (1/3) The filter uses the full non-linear evolution model equation (1) to produce an a priori estimate: x a ( t ) = M t ( x est ( t − 1)) . In order to obtain the corresponding covariance S a ( t ) of the a priori information, the prediction model is linearized about x est ( t − 1) : ∂ M t ( x est ( t − 1)) = (10) M t ∂ x M t S est ( t − 1) M T S a ( t ) = t + SE t . (11) The linearization in equation (10) is produced by taking finite differences between evolved small state vector perturbations, which is a computationally very expensive operation for large models, unless an – Typeset by Foil T EX – 7

adjoint code to the model M has been produced. In both cases, integrating all the columns in M forward in time in equation (11) is expensive. – Typeset by Foil T EX – 8

Extended Kalman filter algorithm (2/3) The observation operator is linearized at the time of the observation about the a priori estimate x a ( t ) in order to obtain K t , which is then used for calculating the gain matrix: ∂ K t ( x a ( t )) = (12) K t ∂ x S a ( t ) K T t ( K t S a ( t ) K T t + Se t ) − 1 . = (13) G t – Typeset by Foil T EX – 9

Extended Kalman filter algorithm (3/3) After this, the full non-linear observation operator is used to update x a ( t ) and this is then used to produce a current state estimate and the corresponding error estimate: x est ( t ) = x a ( t ) + G t ( y ( t ) − K t ( x a ( t ))) . (14) S est ( t ) = S a ( t ) − G t K t S a ( t ) . (15) If the linearization of the observation operator at x a ( t ) is not good enough to construct x est ( t ) , it is necessary to carry out some iterations of the last four equations. – Typeset by Foil T EX – 10

Three-dimensional variational method (3D-Var) The idea of 3D-Var method is to avoid the computation of the gain completely by looking for the analysis as an approximate solution to the equivalent minimization problem defined by the cost function: ( x a ( t ) − x t ) T ( S a ( t )) − 1 ( x a ( t ) − x t ) J ( x t ) = (16) ( y t − K ( x t )) T Se − 1 + t ( y t − K ( x t )) , (17) where K is the nonlinear observation operator. Opposite to the EKF method, the solution is retrieved iteratively by performing several evaluations of the cost function and of its gradient using a suitable descent algorithm. The approximation lies in the fact that only a small number of iterations are performed. – Typeset by Foil T EX – 11

Illustration of the variational cost-function minimization – Typeset by Foil T EX – 12

The four-dimensional variational method (4D-Var) 4D-Var is a generalization of the less expensive 3D-Var method and it also supersedes earlier statistical interpolation methods, which do not take into account model dynamics. In its general form, it is defined as the minimization of the following cost function: J ( x 0 ) = J b + J o (18) ( x b − x 0 ) T S a ( t 0 ) − 1 ( x b − x 0 ) = (19) T � ( y ( t ) − K t ( M t ( x 0 ))) T Se − 1 + t ( y ( t ) − K t ( M t ( x 0 ))) . (20) t =0 4D-Var method assumes that the model is perfect over the assimilation period T . This leads to the so-called strong constraint formalism. – Typeset by Foil T EX – 13

Illustration of the 4D-Var assimilation – Typeset by Foil T EX – 14

Approximating EKF � (Strong constraint) 4DVAR : � S est ( t ) is a static background error matrix S a ( t 0 ) and there is no model error S est ( t ) ≡ S a ( t 0 ) – Typeset by Foil T EX – 15

Approximating EKF � Reduced Rank Kalman Filtering RRKF : � RRKF behaves like the Extended Kalman Filter in a low dimensional subspace determined at the beginning of the assimilation window � ˆ � F T S S est ( t ) = M t L T LM T t F I � L is the orthogonal matrix that transforms model variables into control variables. ˆ S is the model error covariance matrix in the chosen k -dimensional subspace of the control space and F is its cross-correlation matrix with the complement subspace – Typeset by Foil T EX – 16

Approximating EKF � Fast Kalman Filtering FKF : � FKF produces results numerically identical to the full Kalman filter EKF but has one order of magnitude lower complexity � FKF complexity is still superlinear, but it suits medium sized systems well – Typeset by Foil T EX – 17

Approximating EKF � Ensemble Kalman Filtering EnKF : � A low dimensional approximation En − S est to S a ( t 0 ) is propagated by the full nonlinear model and there is no model error � ( M ( s i )( t ) − M ( s i )( t ))( M ( s i )( t ) − M ( s i )( t )) T En − S est ( t ) = s i ∈ V � V is a sample of vectors, such as a multinormal sample modulated by a low rank approximation to the analysis error covariance matrix S a ( t 0 ) at initial time – Typeset by Foil T EX – 18

Variational Kalman filter � The idea of the VKF method is to replace the most time consuming part of the Extended Kalman Filter (EKF) method, the updating of the prediction error covariance matrix S est ( t − 1) equation (11) in order to obtain S a ( t ) , by a low rank update that only requires a few dozen model integrations. � The Variational Kalman Filter can be viewed as a predictor corrector scheme in the Cartesian sum of the state and covariance spaces. In the scheme, a 3D-Var method is used to find both the state estimate and the error covariance matrix. � Dominant Hessian singular vectors at time t + 1 are approximated by evolving the search directions from time t . � The rank of the error covariance approximation can be regulated. – Typeset by Foil T EX – 19

Variational Kalman filter algorithm (1/3) The VKF method uses the full non-linear prediction equation (1) to construct an a priori estimate from the previous state estimate: x a ( t ) = M t ( x est ( t − 1)) . (21) The corresponding approximated covariance S a ( t ) of the a priori information is available from previous time step of VKF method. In order to avoid the computation of the Kalman gain we perform a 3D-Var optimization with equivalent cost function. As a result of the optimization we get the state estimate x est ( t ) for present time t . – Typeset by Foil T EX – 20