CITES-2019 Bayesian approach to the data assimilation problem based - - PowerPoint PPT Presentation

SLIDE 1

Bayesian approach to the data assimilation problem based

• n the use of ensembles of forecasts and observations

Ekaterina Klimova ICT SB RAS

CITES-2019

SLIDE 2

Introduction

• The task of data assimilation is usually understood as the time-sequential

estimation of an unknown quantity from observational data.

• The purpose of data assimilation - both the preparation of initial fields for

subsequent forecasting and the more general one - the description of the behavior of the studied fields over time, the study of climate, the estimation

• f parameters, etc.

SLIDE 3

Bayesian approach to the data assimilation problem

The time change of the estimated quantity:

1 1, (

)

k k k k k

f

 

  x x η

The observations:

( )

k k k k

h   y x ε

k

η

k

ε

• random errors of forecast and observations

The Bayesian approach consists of applying the Bayes theorem to obtain an

• ptimal estimate from observational data and a forecast:

( | ) ( ) ( | ) ( ) p p p p  y x x x y y

,

SLIDE 4

Bayesian approach to the data assimilation problem

There are various options for assessing the state of data and forecast:

,1

( | ),

l k

p k l  x y

• forecast,

,1

( | )

k k

p x y

• filtration,

,0 ,1

( | )

k k

p x y

• smoothing,

where

,0 1

{ , , , }

k k k

 x x x x

,

,1 1

{ , , }

k k

 y y y

.

SLIDE 5

Bayesian approach to the data assimilation problem: the ensemble Kalman filter

Consider a nonlinear dynamic system

1 1

( )

t t t k k k

f

 

  x x η

An observation equation

( )

t t k k k

h   y x ε

t k

ε

1 t k

η

are Gaussian random variables:

   

1

0,

t t k k

E E



  ε η

   

1 1 1

[ ] , [ ]

T T t t t t t t k k k k k k

E E

  

  ε ε R η η Q

and The ensemble Kalman filter consists of an ensemble of forecasts

,

{ , 1, , }

f n k

n N  x

, , 1 1

( )

f n a n n k k k

f

 

  x x η and an ensemble of analyses

,

{ , 1, , }

a n k

n N  x

, , ,

( ( ))

a n f n n n f n k k k k k k

h     x x K y ε x

SLIDE 6

Bayesian approach to the data assimilation problem: the ensemble Kalman filter

k

K

is a matrix of the form

1

( )

f T f T k k k k k k k 

  K P H H P H R

   

, , 1 1

1 1 , 1 1

N N T T f f n f n n n k k k k k k n n

N N

 

 

 

P dx dx R ε ε

, , ,

{ , 1, , }

f n f n f n k k k

n N    dx x x

, , 1

1

N f n f n k k n

N







x x

• an ensemble of forecast errors

{ , 1, , }

n k n

N  ε

• an ensemble of observation errors

1

{ , 1, , }

n k

n N



 η

• an ensemble of model noise

 

1 1 1

[ ]

T n n k k k

E

  

 η η Q

k

H is the linearized operator

,

( ) ( )

f n f k k k k

h h   x x H ε

The analysis error covariance matrix

 

, , 1

1 1

N T a a n a n k k k n

N



  P dx dx

, , ,

{ , 1, , }

a n a n a n k k k

n N    dx x x

, , 1

1

N a n a n k k n

N







x x

SLIDE 7

Bayesian approach to the data assimilation problem

• 1. Given a large sample of realizations for each of the prior pdfs, the

joint pdfs can be evaluated by integration of each individual realization forward in time using stochastic model equation.

• 2. The prior pdfs do not need to be Gaussian distributed.
• 3. The analysis step of EnKF consists of the updates performed on

each of the model state ensemble members.

• 4. For the case with a linear dynamical model a Gaussian prior pdfs

the variance minimizing analysis equals the maximum likelihood estimate.

• 5. For a nonlinear dynamical model the pdfs for the model evolution

will become non-Gaussian. In this case analysis will provide only an approximate solution.

SLIDE 8

Approaches to the implementation of ensemble Kalman filter

Ensemble of forecasts

( ) ( ) 0( ) ( )

( )

a i f i i f i k k k k k k

x x K y H x   

Ensemble of analysis

( ) ( ) ( ) ( )

,

a i a a i a i a i T a

x x dx dx dx P   

Estimation error (skill)

( ) ( ) ( ) ( )

{ }, { }

a a i a i f f i f i

X x x X x x    

Ensemble spread

( ) ( )

{ }, { }

a t a i f t f i

x x x x      

“True value”

1

( )

t t t k k k k

x f x 

 



«Stochastic filter» (EnKF) «Deterministic filter» (ESRF, ETKF, LETKF)

( ) ( ) ( ) 1

( )

f i a i i k k k k

x f x 

 



( ) ( ) a i f i k k k

dx A dx 

The transformation of ensembles of forecast

SLIDE 9

Practical implementation of ensemble algorithms

• 1. Algorithms with the transformation of forecast ensembles.
• 2. Local algorithms.
• 3. Methods to increase ensemble spread.

SLIDE 10

Local Ensemble Transform Kalman Filter- LETKF (Hunt et al, 2007)

( ) ( ) ( ) ( )

{ } { }

a a i a i f f i f i

X x x X x x    

• the ensembles of analysis and

forecasts

 

1 1 1/2 1 ( ) ( ) ( )

1) ( 1) / ( ) , 2) ( 1) 3) ( ), ( ) 4) , 5) 6)

a f T f f f a a a a f f T

• a i

f f a i a i T a f f a a f a f

P k I Y R Y Y HX W k P w P C y Hx C Y R x x X w w i th row of W x x X w P X P X 

  

                       

LETKF:

• “inflation

factor”.



SLIDE 11

Ensemble π-algorithm

The ensemble π-algorithm is a stochastic filter in which the analysis step is performed

• nly for the ensemble mean.

The ensemble of analysis errors D is a matrix the columns of which are vectors

{ , 1, , }

n k n

N  dx

T T 1 T

( ) ,



  D I Π B

1 T 2

( 0,25 ) 0,5 ,    Π C I I B is a matrix with columns {

, 1, , }

n k n

N  b

, , n f n f n k k k

  b x x

Ε

T T 1 1 2

1 ( ) . 1 N



     C F H R HB Ε C C

is a matrix with columns

n k

ε

• the ensemble of observation errors

Klimova E. A suboptimal data assimilation algorithm based on the ensemble Kalman

• filter. Quarterly Journal of the Royal Meteorological Society. 2012.

DOI:10.1002/qj.1941. Klimova E.G. The Kalman stochastic ensemble filter with transformation of perturbation ensemble. Siberian. J. Nun, Math. /Sib. Branch of Russ. Acad. of Sci. – Novosibirsk, 2019. Vol. 22 N 1. P.27-40.

SLIDE 12

Classical particle filter

Ensemble

( ) l

x

• f states representing the prior probability distribution

b k

p at time k t . The analysis step at time k t : calculation of new weights

( ) ( ) ( ) ( , )

( ) ( | ) ( )

a l l b l k l k k k

w p x cp y x p x  

SLIDE 13

Gaussian particle filter

The Gaussian particle filter treats each particle as a Gaussian probability distribution

( , ) 1 1 1

( | ) ( , )

L b l k k k l

p x Y N x B

  



( , ) 1

( | ) ( , )

L a l k k k l

p x Y N x B



 The analysis ensembles

( , ) a l

x are calculated by treating each particle as an individual Gaussian distribution:

( , ) ( , ) ( , ) 1 1 1

( ) ( ) ( )

a l b l b l k k k T k k k k k k k k k k

x x K y H x K B H H B H R B I K H B

  

      

SLIDE 14

Nonlinear ensemble filter (T.Bengtsson, C.Snyder, D.Nychka J. of Geoph. Res. V. 108 No D24 2003)

Suppose that

1 , , , 1

( | ) ( , )

L f f f k k k l k l k l l

p x Y N x P 

 

 

, , , 1

( | ) ( , )

L a a a k k k l k l k l l

p x Y N x P 





, , , , , , , ,

( ) ( ) ( )

a f f k l k l k l k k k l f f T k k l k k k l k a f k l k k k l

x x K y H x K P H H P H R P I K H P       

, , , 1 f k l l a k l L f k j j j

w w   







1/2 1 , , , ,

( ) exp[ 1/ 2( ) ( ) ( )]

f T f T f T f l k k l k k k k l k k l k k k k l

w H P H R y H x H P H R y H x

 

     

Ensemble of analysis

, , , ,

( )

a f f k j k j k j k j k k j

x x K y H x      Sampling according to

( | )

k k

p x Y

SLIDE 15

The behavior of ensemble spread in the ensemble Kalman filter (stochastic filter)

The stochastic ensemble Kalman filter can be written in the following form:

1 1

( ) ( ) ( )

n n n n n k k k k k k k k

f

 

         x I K H x η K y ε

The optimal estimate in the ensemble Kalman filter is the ensemble mean value

n k

x

Deviation from the mean (spread) simulates the estimate error

n n n k k k

  dx x x

1 1 1

( )( ( ) ( ) )

n n n n n k k k k k k k k

f f

  

     dx I K H x x η K ε

A ‘theoretical’ estimation error (skill)

t t n k k k

  dx x x

1 1 1

( )( ( ) ( ) )

t t n t t k k k k k k k k

f f

  

     dx I K H x x η K ε

SLIDE 16

The behavior of ensemble spread in the ensemble Kalman filter (deterministic filter)

The deterministic ensemble Kalman filter (analysis step) consists of the equation for the mean

, , ,

( ( ))

a n f n n f n k k k k k

h    x x K y x

and an estimate of the ensemble of analysis errors such that the corresponding covariance matrix satisfies the Kalman filter equation

( )

a f k k k k

  P I K H P

The transformation of an ensemble of forecast errors into analysis errors in a deterministic filter can be represented in the form of left multiplication

, , a n f n k k k

 dx A dx

f T a k k k k

 A P A P

where The ensemble of analyses of the deterministic filter:

, , , ,

( )

a n a n f n f n k k k k k

   x x A x x

SLIDE 17

The equation for ensemble spread in the ensemble Kalman filter (stochastic filter)

The equation for ensemble spread in the stochastic Kalman filter

1 1

( )( )

n n n n k k k k k k k k  

    dx I K H F dx η K ε

Instead of the nonlinear model operator f we take the linearized operator

k

F

1 1 1

( ,0) ( , ) ( , )( )

k k n n n n k i i i i i i i

k k i k i

  

   

 

dx Ψ dx Ψ K ε Ψ I K H η

1

( , ) ( )

k j j j j i

k i

 

 



Ψ I K H F k i 

( , ) k k  Ψ I

SLIDE 18

The estimation error

The estimation error (the deviation of the mean from the ‘true’ value)

1 1

( )( )

t t t t k k k k k k k k  

    dx I K H F dx η K ε

The simulated estimation error tends to the theoretical error if the random vectors

• f observational errors and model noise being simulated have the same covariance

matrices as the true ones.

SLIDE 19

The behavior of ensemble spread in the ensemble Kalman filter (deterministic filter)

Writing the formula for the analysis perturbations in terms of ‘left multiplication’, we obtain the following equation for ensemble spread in the deterministic Kalman filter:

1 1

( )

n n n k k k k k  

  dx A F dx η

det det 1 1

( ,0) ( , )

k n n n k i i i

k k i

 

  dx Ψ dx Ψ A η

det 1

( , )

k j j j i

k i

 



Ψ A F

The formula for the deterministic filter lacks the term with

n i i

K ε

which simulates ensemble spread as a function of observational data distribution and

• bservational error covariances.

SLIDE 20

Equation for estimation error (particle filter)

1

ˆ

N a a a j j j

x x 



 ( )

a f

• j

j

I KH K      

1

( )

L a f

• t

j j j

I KH K    



  



Given N independent samples from a density p , an estimator of p can be obtained as a mixture of N Gaussian densities. In that case

1,

,

N

x x

SLIDE 21

Methods of improving convergence in the ensemble Kalman filter

Some of the most frequently used methods of improving convergence of the ensemble Kalman filter are multiplicative inflation and additive inflation. Let us consider the ensemble spread modification in general form. In the case of analysis step it has the form

n n n k k k k

   dx dx β

n k

β

is a random vector with a specified covariance matrix in the case of forecast step, the form

1 1

( )( )

n n n n n k k k k k k k k k k k

 

 

     dx I K H F dx η β K ε

SLIDE 22

Methods of improving convergence in the ensemble Kalman filter

In the case of analysis step the formula for a stochastic filter

1 1 1 1 1

( ,0) ( , ) ( , ) ( ) ( , )

k k k n n n n n k i i i i i i i i i i i

k k i k i k i  

    

    

  

dx Ψ dx Ψ K ε Ψ I K H η Ψ β

in the case of forecast step

1 1 1 1 1

( ,0) ( , ) ( , )( ) ( , )( )

k k k n n n n n k i i i i i i i i i i i i

k k i k i k i 

    

     

  

dx Ψ dx Ψ K ε Ψ I K H η Ψ I K H β

1

( , ) ( )

k j j j j j i

k i 

 





Ψ I K H F

j

K

is calculated using the modified covariance matrices

SLIDE 23

Methods of improving convergence in the ensemble Kalman filter

For a deterministic filter in the case of analysis step

det det 1 det 1 1

( ,0) ( , ) ( , )

k k n n n n k i i i i i i

k k i k i 

  

  

 

dx Ψ dx Ψ A η Ψ β in the case of forecast step

det det 1 det 1 1

( ,0) ( , ) ( , )

k k n n n n k i i i i i i i

k k i k i 

  

  

 

dx Ψ dx Ψ A η Ψ A β

det 1

( , )

k j j j j i

k i 

 

  Ψ A F

The equation for the error when modifying the ensemble spread:

1 1 1

( ,0) ( , ) ( , )( )

k k t n t t k t t i i t i i i i i

k k i k i

  

   

 

dx Ψ Δ Ψ K ε Ψ I K H η

1

( , ) ( )

k t j j j j i

k i

 





Ψ I K H F

SLIDE 24

Methods of improving convergence in the ensemble Kalman filter

• 1. The perturbation ensembles of deterministic and stochastic filters with the

thus modified ensemble spread do not correspond to the error ensemble.

• 2. For we obtain a version of multiplicative inflation.
• 3. For we obtain a version of additive inflation. In this case, additive

inflation can be specified so that the covariance matrix coincides with the matrix obtained when using multiplicative inflation:

n k 

β 1

k

 

n n k k

 β δQ ξ

2

( 1)

k k k

   δQ P

SLIDE 25

Numerical experiments. Lorenz-96 model

The equations of the model

1 2 1 1 1 1 1

d ( ) , 1, , dt , x ,

j j j j j J J

x x x x x F j J x x x

     

      

1,

,

J

x x

(J=40) are the variables being forecasted; a fourth-order finite-difference Runge-Kutta scheme;

0.05 Δt 

corresponds to 6 hours (t=1 is taken for five days);

F =8.

To simulate ‘true values’ in the numerical data assimilation experiments,

t

N(F /4;F /2)  x

and =1000 time steps are made.

t

N

Initial data for forecasting:

d t

(0) 0) , ( N(0,s )   x x δ δ

SLIDE 26

Numerical experiments

The following parameters are specified for the numerical experiments: an ensemble of initial fields:

d

(0) (0) ,

n n n

N(0,s ),n=1, ,N   x x δ δ

• bservations:

t (0)

( ) , N 0,ε   y x δ δ

an ensemble of observations with perturbations:

( )

n n n

, N 0,ε ,n=1, ,N   y y δ δ

model noise:

n 

η

in simulating the ‘truth’: 0.01

t 

η The observations are available at each of the J=40 model grid points. The experimental period has a length 3000 time steps, with assimilation being done at each time step or at every four time steps.

1 s =ε 

N=20

2

  R I

in all experiments:

SLIDE 27

Numerical experiments

The numerical experiments are performed for ten versions of the ‘truth’, and all estimates were calculated as average values over these ten versions. The following estimates were considered:

12 2 , , 1 1

1 1 ( )

K L t k i k i k i

rms x x K L

 

       

 

• the root-mean-square error averaged over K=10

versions of calculations (k is the number of a version and i is the number of a grid node)

12 2 , , 1 1 1

1 1 ( ) ( 1)

K L N n n k i k i k i n

sp x x K L N

  

        

 

• the mean value of the covariance matrix trace

calculated over K=10 versions of calculations (n is the number of an ensemble member).

SLIDE 28

Numerical experiments

Two series of experiments were performed. In the first series, observational data were simulated at one time step intervals, and in the second series, at four time step intervals. The following experiments were performed in each of the series: Experiment 1. In this experiment, multiplicative inflation was used. =1.1 in the first series and =1.2 in the second series). Experiment 2. In this experiment, additive inflation was used so that the change in the covariance matrix coincided with the change made in experiment 1.

 



SLIDE 29

Numerical experiments

2800 2820 2840 2860 2880 2900 2920 2940 2960 2980 3000 0.15 0.16 0.17 0.18 0.19 0.2 0.21 time step rms-sp rms: mult.infl. rms: additive infl. sp: mult.infl. sp: additive infl.

The results of the first series of experiments for the ensemble π-algorithm

SLIDE 30

Numerical experiments

2800 2820 2840 2860 2880 2900 2920 2940 2960 2980 3000 0.15 0.16 0.17 0.18 0.19 0.2 0.21 0.22 0.23 0.24 time step rms-sp rms: mult.infl. rms: additive infl. sp: mult.infl. sp: additive infl.

The results of the first series of experiments for the LETKF algorithm

SLIDE 31

Numerical experiments

2800 2820 2840 2860 2880 2900 2920 2940 2960 2980 3000 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 time step rms-sp rms: mult.infl. rms: additive infl. sp: mult.infl. sp: additive infl.

The results of the second series of experiments for the ensemble π-algorithm

SLIDE 32

Numerical experiments

2800 2820 2840 2860 2880 2900 2920 2940 2960 2980 3000 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 time step rms-sp rms: mult.infl. rms: additive infl. sp: mult.infl. sp: additive infl.

The results of the second series of experiments for the LETKF algorithm.

SLIDE 33

Conclusion

• In the ensemble approach, at the analysis step it is important to specify

the ensembles corresponding to the density of the analysis error

• distribution. It is necessary to take into account the ensemble of errors of
• bservation.
• To regulate the convergence of ensemble algorithms, it is preferable to

use additive inflation.

• The results of the investigations show that ensemble spread in stochastic

filters rather than in deterministic filters is closer to the theoretical estimation error.

• Multiplicative inflation and additive inflation change the general formula

for ensemble spread.

• The formula for ensemble spread in a stochastic filter with additive

inflation rather than with multiplicative inflation is closer to that for the estimation error.

SLIDE 34