High-Dimensional Nonlinear Data Assimilationwith the Nonlinear Ensemble Transform Filter (NETF)and its Smoother Extension

(1)

High-Dimensional Nonlinear Data Assimilation

with the Nonlinear Ensemble Transform Filter (NETF) and its Smoother Extension

Lars Nerger

,

Paul Kirchgessner,

Alfred Wegener Institute, Bremerhaven, Germany

Julian Tödter, Bodo Ahrens

University of Frankfurt, Frankfurt, Germany NMEFC, Beijing, China, November 9, 2017

(2)

Nonlinear Ensemble Transform Filter & Smoother

Overview

Ø Study new Nonlinear Ensemble Transform Filter – NETF (Tödter & Ahrens, MWR, 2015)

Ø Extend NETF for smoothing

Ø Test filter and smoother in realistic high-dimensional

idealized ocean data assimilation experiments

(3)

Kalman and Nonlinear Filters

(4)

• represent state and its error by ensemble of states

• Forecast:

• Integrate ensemble with numerical model

• Analysis:

• update ensemble mean

• update ensemble perturbations

(both can be combined in a single step)

• Ensemble Kalman filters & NETF: Different definitions of

• weight vector

• Transform matrix

Ensemble filters – ensemble Kalman filters & NETF

X

w ˜

W

x

^a

= x

^f

+ X

⁰^f

w ˜ X

⁰^a

= X

⁰^f

W

N

(5)

a

The Ensemble Kalman Filter (EnKF, Evensen 94)

Ensemble

Analysis step:

Update each ensemble member

Kalman filter

5 EnKF

Init

x

^a₀

⌅ R

ⁿ

, P

^a₀

⌅ R

ⁿ^⇥ⁿ

(41) { x

^a(l)₀

, l = 1, . . . , N } (42) x

^a₀

= 1

N

⇧

N

l=1

x

^a(l)₀

⇥ x

^t₀

⇥

(43)

P ˜

^a₀

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)₀

x

^a₀

⌅⇤

x

^a(l)₀

x

^a₀

⌅

T

⇥ P

^a₀

(44)

P

^a₀

= LL

^T

, L ⌅ R

ⁿ^⇥^q

(45) x

^a(i)₀

= x

^a₀

+ Lb

⁽ⁱ⁾

, b

⁽ⁱ⁾

⌅ R

^q

(46)

⇤ N (0, 1) (47)

Forecast

x

^a(l)_i

= M

_i,i ₁

[x

^a(l)_i ₁

] +

^(l)_i

(48)

Analysis

{ y

_k^o(l)

, l = 1, . . . , N } (49) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

_k^o(l)

H

_k

⌃

x

^f_k^(l)

⌥⌅

(50) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

^o(l)_k

H

_k

x

^f_k^(l)

⌅

(51) K ˜

_k

= P ˜

^f_k

H

^T_k

⇤

H

_k

P ˜

^f_k

H

^T_k

+ R

_k

⌅

1

(52) K

_k

= P

^f_k

H

^T_k

⇤

H

_k

P

^f_k

H

^T_k

+ R

_k

⌅

1

(53)

H

_k

P

^f_k

H

^T_k

+ R

_k

⌅ R

^m^⇥^m

(54) P ˜

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(55)

x

^a_k

:= 1 N

⇧

N

l=1

x

^a(l)_k

(56)

P ˜

^a_k

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)_k

x

^a_k

⌅⇤

x

^a(l)_k

x

^a_k

⌅

T

(57)

5

5 EnKF

Init

x

^a₀

⌅ R

ⁿ

, P

^a₀

⌅ R

ⁿ^⇥ⁿ

(41) { x

^a(l)₀

, l = 1, . . . , N } (42) x

^a₀

= 1

N

⇧

N

l=1

x

^a(l)₀

⇥ x

^t₀

⇥

(43)

P ˜

^a₀

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)₀

x

^a₀

⌅⇤

x

^a(l)₀

x

^a₀

⌅

T

⇥ P

^a₀

(44)

P

^a₀

= LL

^T

, L ⌅ R

ⁿ^⇥^q

(45) x

^a(i)₀

= x

^a₀

+ Lb

⁽ⁱ⁾

, b

⁽ⁱ⁾

⌅ R

^q

(46)

⇤ N (0, 1) (47)

Forecast

x

^a(l)_i

= M

_i,i ₁

[x

^a(l)_i ₁

] +

^(l)_i

(48)

Analysis

{ y

_k^o(l)

, l = 1, . . . , N } (49) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

^o(l)_k

H

_k

⌃

x

^f_k^(l)

⌥⌅

(50) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

^o(l)_k

H

_k

x

^f_k^(l)

⌅

(51) x

^a(l)_k

= x

^f_k^(l)

+ K

_k

⇤

y

^(l)_k

H

_k

x

^f_k^(l)

⌅

(52) K ˜

_k

= P ˜

^f_k

H

^T_k

⇤

H

_k

P ˜

^f_k

H

^T_k

+ R

_k

⌅

1

(53) K

_k

= P

^f_k

H

^T_k

⇤

H

_k

P

^f_k

H

^T_k

+ R

_k

⌅

1

(54)

H

_k

P

^f_k

H

^T_k

+ R

_k

⌅ R

^m^⇥^m

(55) P ˜

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(56)

x

^a_k

:= 1 N

⇧

N

l=1

x

^a(l)_k

(57)

P ˜

^a_k

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)_k

x

^a_k

⌅⇤

x

^a(l)_k

x

^a_k

⌅

T

(58)

5 EnKF

Init

x

^a₀

⌅ R

ⁿ

, P

^a₀

⌅ R

ⁿ^⇥ⁿ

(41) { x

^a(l)₀

, l = 1, . . . , N } (42) x

^a₀

= 1

N

⇧

N

l=1

x

^a(l)₀

⇥ x

^t₀

⇥

(43)

P ˜

^a₀

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)₀

x

^a₀

⌅⇤

x

^a(l)₀

x

^a₀

⌅

T

⇥ P

^a₀

(44)

P

^a₀

= LL

^T

, L ⌅ R

ⁿ^⇥^q

(45) x

^a(i)₀

= x

^a₀

+ Lb

⁽ⁱ⁾

, b

⁽ⁱ⁾

⌅ R

^q

(46)

⇤ N (0, 1) (47)

Forecast

x

^a(l)_i

= M

_i,i ₁

[x

^a(l)_i ₁

] +

^(l)_i

(48)

Analysis

{ y

^o(l)_k

, l = 1, . . . , N } (49) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

^o(l)_k

H

_k

⌃

x

^f_k^(l)

⌥⌅

(50) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

^o(l)_k

H

_k

x

^f_k^(l)

⌅

(51) x

^a(l)_k

= x

^f_k^(l)

+ K

_k

⇤

y

_k^(l)

H

_k

x

^f_k^(l)

⌅

(52) K ˜

_k

= P ˜

^f_k

H

^T_k

⇤

H

_k

P ˜

^f_k

H

^T_k

+ R

_k

⌅

1

(53) K

_k

= P

^f_k

H

^T_k

⇤

H

_k

P

^f_k

H

^T_k

+ R

_k

⌅

1

(54) H

_k

P

^f_k

H

^T_k

+ R

_k

⌅ R

^m^⇥^m

(55) P ˜

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(56)

P

^f_k

:= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(57)

x

^a_k

:= 1 N

⇧

N

l=1

x

^a(l)_k

(58)

P ˜

^a_k

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)_k

x

^a_k

⌅⇤

x

^a(l)_k

x

^a_k

⌅

T

(59)

5

5 EnKF

Init

x

^a₀

⌅ R

ⁿ

, P

^a₀

⌅ R

ⁿ^⇥ⁿ

(41) { x

^a(l)₀

, l = 1, . . . , N } (42) x

^a₀

= 1

N

⇧

N

l=1

x

^a(l)₀

⇥ x

^t₀

⇥

(43)

P ˜

^a₀

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)₀

x

^a₀

⌅⇤

x

^a(l)₀

x

^a₀

⌅

T

⇥ P

^a₀

(44)

P

^a₀

= LL

^T

, L ⌅ R

ⁿ^⇥^q

(45) x

^a(i)₀

= x

^a₀

+ Lb

⁽ⁱ⁾

, b

⁽ⁱ⁾

⌅ R

^q

(46)

⇤ N (0, 1) (47)

Forecast

x

^a(l)_i

= M

_i,i ₁

[x

^a(l)_i ₁

] +

^(l)_i

(48)

Analysis

{ y

^o(l)_k

, l = 1, . . . , N } (49) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

_k^o(l)

H

_k

⌃

x

^f_k^(l)

⌥⌅

(50) x

^a(l)_k

= x

^f_k^(l)

+ K ˜

_k

⇤

y

_k^o(l)

H

_k

x

^f_k^(l)

⌅

(51) x

^a(l)_k

= x

^f_k^(l)

+ K

_k

⇤

y

_k^(l)

H

_k

x

^f_k^(l)

⌅

(52) K ˜

_k

= P ˜

^f_k

H

^T_k

⇤

H

_k

P ˜

^f_k

H

^T_k

+ R

_k

⌅

1

(53) K

_k

= P

^f_k

H

^T_k

⇤

H

_k

P

^f_k

H

^T_k

+ R

_k

⌅

1

(54) H

_k

P

^f_k

H

^T_k

+ R

_k

⌅ R

^m^⇥^m

(55) P ˜

^f_k

= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(56)

P

^f_k

:= 1

N 1

⇧

N

l=1

⇤ x

^f_k^(l)

x

^f_k

⌅⇤

x

^f_k^(l)

x

^f_k

⌅

T

(57)

x

^a_k

:= 1 N

⇧

N

l=1

x

^a(l)_k

(58)

P ˜

^a_k

:= 1

N 1

⇧

N

l=1

⇤ x

^a(l)_k

x

^a_k

⌅⇤

x

^a(l)_k

x

^a_k

⌅

T

(59)

5

Ensemble

covariance matrix

Ensemble mean (state estimate)

Expensive to compute

(6)

Efficient use of ensembles

€

Kalman gain

K ˜

_k

= ˜ P

^f_k

H

^T_k

⇣

H

_k

P ˜

^f_k

H

^T_k

+ R

_k

⌘

1

K ˜

_k

= ⇣

P ˜

^f_k

⌘

1

+ H

^T

R

¹

H

1

H

^T

R

¹

Alternative form (Sherman-Morrison-Woodbury matrix identity)

Looks worse: matrices need inversion

n ⇥ n

K ˜

_k

= X

⁰

h

(N 1)I + X

⁰^T

H

^T

R

¹

HX

⁰

i

1

X

⁰^T

H

^T

R

¹

However: with ensemble

Inversion of matrix

(Ensemble perturbation matrix )

P ˜

^f_k

= (N 1)

¹

X

⁰

X

⁰^T

N ⇥ N

X

⁰

= X X ¯

(7)

• Ensemble Transform Kalman filter:

• Transform matrix

• Mean update weight vector

(depends on R and y)

• Transformation of ensemble perturbations

(depends only on R, not y)

ETKF (Bishop et al., 2001)

A

¹

= (m 1)I + (HX

⁰^f

)

^T

R

¹

HX

⁰^f

˜

w = A(HX

⁰^f

)

^T

R

¹

⇣

y Hx

^f

⌘

W = p

(m 1)A

^1/2

⇤

N

(8)

• Avoid changing ensemble members (‘particles’)

• Instead: give particles a weight at change it at the analysis step

• Initial weight: 1/N for all particles

• Weights are given by statistical likelihood of an observation

• Example: With Gaussian observation errors (for each particle i):

• Ensemble mean state computed with weights

• This update does not assume any distribution of the state errors (and is not limited to Gaussian distributations)

Particle filters – fully nonlinear ensemble filters

˜

w

ⁱ

⇠ exp ⇣

0.5(y Hx

^f_i

)

^T

R

¹

(y Hx

^f_i

) ⌘

x

^a

= x

^f

+ X

⁰^f

w ˜ = X

^f

w ˜

(9)

• Ensemble Kalman:

• Transformation according to KF equations

• NETF (Tödter & Ahrens, MWR, 2015)

Ø Mean update from Particle Filter weights: for all particles i

Nonlinear ensemble transform filter - NETF

Ø Ensemble update

• Transform ensemble to fulfill analysis covariance (like KF, but not assuming Gaussianity)

• Derivation gives

( : mean-preserving random matrix; useful for stability) (Almost same formulation: Xiong et al., Tellus, 2006)

⇤

W = p

m ⇥

diag( ˜ w) w ˜ w ˜

^T

⇤

^1/2

⇤

˜

w

ⁱ

⇠ exp ⇣

0.5(y Hx

^f_i

)

^T

R

¹

(y Hx

^f_i

) ⌘

p N

(10)

• Mean state update

• Analysis covariance matrix

with

Derivation of NETF

x

^a

= x

^f

+ X

⁰^f

w ˜ = X

^f

w ˜

P

^a

= X

i=1,m

˜

w

_i

(x

^f_i

x

^a

)(x

^f_i

x

^a

)

^T

W = p

m ⇥

diag(w) w ˜ w ˜

^T

⇤

^1/2

⇤ P

^a

= 1

m X

^f

W

²

(X

^f

)

^T

X

i=1,N

N

p N

(11)

• ETKF parameterizes ensemble distribution by a Gaussian distribution

• NETF uses particle filter weights to ensure correct update of ensemble mean and covariance

• Filter update:

• in ETKF is linear in observations

• in NETF is nonlinear in observations

Difference of ETKF and NETF

˜

w = A(HX

⁰^f

)

^T

R

¹

⇣

y Hx

^f

⌘

˜

w

ⁱ

⇠ exp ⇣

0.5(y Hx

^f_i

)

^T

R

¹

(y Hx

^f_i

) ⌘

(12)

• Smoother: Update past ensemble with future observations

• Rewrite ensemble update as

• Filter:

Ensemble Smoothers – ETKS & NETS

X

^a_k_|_k

= X

^f_k_|_k ₁

W ˆ

_k

analysis time Observations used up to time

• Smoother at time

Ø works likewise for ETKS and NETS Ø also possible for localized filters

X

^a_i_|_k

= X

^f_i_|_k ₁

W ˆ

_k

i < k

See, e.g., Nerger, Schulte & Bunse-Gerstner, QJRMS 140 (2014) 2249–2259

(13)

Experiments with small Lorenz-96 model

(14)

Configuration of Lorenz-96 model experiments

Lorenz-96:

• 1-dimensional period wave

• Chaotic dynamics

Configuration for assimilation experiments

• State dimension: 80

• Observed: 40 grid points

• Time steps between analysis steps: 8

• Double-exponential observation errors (for even stronger nonlinearity)

• Experiment length: 5000 time steps

• Observation error standard deviation: 1

➜ this is a difficult case for the assimilation

www.data-assimilation.net

(15)

• Performance for small model (Lorenz-96)

• NETF beats ETKF for ensemble size larger 30

Performance of NETF – Lorenz-96

20 30 40 50 60 70

ensemble size 1.1

1.2 1.3 1.4 1.5 1.6 1.7 1.8

MRMSE

EKTF filter NETF filter

(16)

• Time period over which smoothing is performed: smoother lag

Typical behavior with nonlinear models

• Fast reduction of error short lag

• Error increase for large lag (caused by nonlinarity)

➜ There is an optimal lag with minimum error

Appliction of smoother

0 50 100 150 200

Lag (time steps) 1

1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45

MRMSE

LETKS LNETS

L. Nerger, S. Schulte, A. Bunse-Gerstner (2014) QJR. Meteorol. Soc. 140: 2249

(17)

• Performance for small model (Lorenz-96)

• Blue: Smoother

• NETS beats ETKS for ensemble size 40 and larger

• Smoother slightly stronger for ETKS

• NETS better than ETKF filter for N=70

Performance of NETF – Lorenz-96

20 30 40 50 60 70

ensemble size 0.9

1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8

MRMSE

EKTF filter ETKS smoother NETF filter NETS smoother

(18)

NETF/NETS with

high-dimensional ocean model

(19)

Assimilation into NEMO

European ocean circulation model

Model configuration

• box-configuration “SEABASS”

• ¼^o resolution

• 121x81 grid points, 11 layers (state vector ~300,000)

• wind-driven double gyre (a nonlinear jet and eddies)

• medium size SANGOMA benchmark

True sea surface height at 1st analysis time

Longitude (degree)

Latitide (degree)

−60 −55 −50 −45 −40 −35 −30

24 28 32 36 40

44 −0.6

−0.4

−0.2 0 0.2 0.4 0.6

True sea surface height at last analysis time

Longitude (degree)

Latitide (degree)

−60 −55 −50 −45 −40 −35 −30

24 28 32 36 40

44 −0.6

−0.4

−0.2 0 0.2 0.4 0.6

(20)

PDAF: A tool for data assimilation

PDAF - Parallel Data Assimilation Framework

§ a program library for ensemble data assimilation

§ provide support for parallel ensemble forecasts

§ provide fully-implemented & parallelized filters and smoothers (EnKF, LETKF, NETF, EWPF … easy to add more)

§ easily useable with (probably) any numerical model

(applied with NEMO, MITgcm, FESOM, HBM, TerrSysMP, …)

§ run from laptops to supercomputers (Fortran, MPI & OpenMP)

§ first public release in 2004; continued development

§ ~250 registered users; community contributions Open source:

Code, documentation & tutorials at http://pdaf.awi.de

L. Nerger, W. Hiller, Computers & Geosciences 55 (2013) 110-118

(21)

single program

Indirect exchange (module/common) Explicit interface

state time

state

observations

mesh data

Model

initialization time integration post processing

Ensemble Filter

Initialization analysis

ensemble transformation

Observations

quality control obs. vector obs. operator

obs. error

Core of PDAF

Logical separation of assimilation system

modify parallelization

Nerger, L., Hiller, W. Software for Ensemble-based DA Systems –

(22)

Extending a Model for Data Assimilation

Extension for data assimilation

revised parallelization enables ensemble forecast

plus:

Possible model-specific

adaption:

for NEMO:

handle leapfrog time

stepping

Start

Stop Do i=1, nsteps

Initialize Model

Initialize coupler Initialize grid & fields

Time stepper

in-compartment step coupling

Post-processing

Model

single or multiple executables coupler might be separate program

Initialize parallel.

Aaaaaaaa

Aaaaaaaa aaaaaaaaa

Stop

Initialize Model

Initialize coupler Initialize grid & fields

Time stepper

in-compartment step coupling

Post-processing Init_parallel_PDAF

Do i=1, nsteps Init_PDAF

Assimilate_PDAF Start

Initialize parallel.

(23)

Features of online-coupled DA program

• minimal changes to model code when combining model with filter algorithm

• model not required to be a subroutine

• no change to model numerics!

• model-sided control of assimilation program (user-supplied routines in model context)

• observation handling in model-context

• filter method encapsulated in subroutine

• complete parallelism in model, filter, and ensemble integrations

Aaaaaaaa Aaaaaaaa aaaaaaaaa

Start

Stop

Initialize Model

generate mesh Initialize fields

Time stepper

consider BC Consider forcing

Post-processing init_parallel_pdaf

Do i=1, nsteps init_pdaf

assimilate_pdaf

(24)

Online coupling: Minimal changes to NEMO

Add to mynode (lin_mpp.F90) just before init of myrank

#ifdef key_USE_PDAF

CALL init_parallel_pdaf(0, 1, mpi_comm_opa)

#endif

Add to nemo_init (nemogcm.F90) at end of routine

#ifdef key_USE_PDAF

CALL init_pdaf()

#endif

Add to stp (step.F90) at end of routine

#ifdef key_USE_PDAF

CALL assimilate_pdaf()

#endif

Modify dyn_nxt (dynnxt.F90)

#ifdef key_USE_PDAF

IF((neuler==0 .AND. kt==nit000).OR.assimilate)

#else

Aaaaaaaa Aaaaaaaa aaaaaaaaa

Start

Stop

Initialize Model

generate mesh Initialize fields

Time stepper

consider BC Consider forcing

Post-processing init_parallel_pdaf

Do i=1, nsteps init_pdaf

assimilate_pdaf