Transcript Data Assimilation: Data assimilation seeks to characterize

Advanced data assimilation methodsEKF and EnKF
Hong Li and Eugenia Kalnay
University of Maryland
17-22 July 2006
Recall the basic formulation in OI
• OI for a Scalar:
0  w 1
Ta  Tb  w(T0  Tb )
Optimal weight to minimize the analysis error is:
 2f
w 2
 f   o2
• OI for a Vector:
x  Mx
f
i
a
i 1
t=0
t=6h
True atmosphere
xia  xbi  W(yio  Hxbi )
W  BHT [HBHT  R]1
Pa  [I  WH]B
B   xb   xTb
 xb  xb  x t
• B is statistically pre-estimated, and constant with time in its practical
implementations. Is this a good approximation?
OI and Kalman Filter for a scalar
• OI for a scalar:
Tb (ti 1 )  M (Ta (ti );
 b2  (1  a) a2 
 b2
2
2
w 2
;


(1

w)

a
b
 b   o2
1
 a2  const.
1 w
Ta  Tb  w(T0  Tb )
• Kalman Filter for a scalar:
Tb (ti1 )  M(Ta (ti ));  b2 (t)  (L a )(L a )T ; L  dM / dT
 b2
w 2
;  a2  (1  w) b2
2
b  o
Ta  Tb  w(T0  Tb )
• Now the background error variance is forecasted using the linear
tangent model L and its adjoint LT
“Errors of the day” computed with the Lorenz 3-variable
model: compare with rms (constant) error
 z2  (zb  zt )2
=0.08
Not only the amplitude, but also the structure of B is constant in 3D-Var/OI
This is important because analysis increments occur only within the subspace
spanned by B
“Errors of the day” obtained in the reanalysis
(figs 5.6.1 and 5.6.2 in the book)
Note that the mean error went down from 10m to 8m because of
improved observing system, but the “errors of the day” (on a synoptic time
scale) are still large.
In 3D-Var/OI not only the amplitude, but also the structure of B is fixed
with time
Flow independent error covariance
If we observe only Washington, D.C,
we can get estimate for Richmond and
Philadelphia corrections through the
error correlation (off-diagonal term in B).
  12

cov1, 2

B
 ...

cov1,n
cov1, 2
 22
...
cov2,n
... cov1,n 

cov2,n 
... 

...  n2 
In OI(or 3D-Var), the scalar error correlation between two points in the
same horizontal surface is assumed homogeneous and isotropic. (p162 in
the book)
Typical 3D-Var error covariance
  12

cov1, 2

B
 ...

cov1,n
cov1, 2
 22
...
cov2,n
... cov1,n 

cov2,n 
... 

...  n2 
In OI(or 3D-Var), the error correlation between two mass points in the
same horizontal surface is assumed homogeneous and isotropic.(p162
in the book)
Flow-dependence – a simple example (Miyoshi,2004)
There is a cold front in
our area…
What happens in this
case?
This is not appropriate
This does reflects the flow-dependence.
Extended Kalman Filter (EKF)
• Forecast step
xbi  Mxia1
Pib  Li 1Pia1LTi1  Q
* (derive it)
• Analysis step
xia  xbi  Ki (yio  Hxif )
Ki  PibHT [HPib HT  R]1
**
Pia nn  [I  Ki H]Pib  [(Pib )1  HT R1H]1
•
b
Using the flow-dependent Pi , analysis is expected to be improved
significantly
b
However, it is computational expensive. Pi , L i
computing equation
directly is impossible
*
n*n matrix, n~107
Ensemble Kalman Filter (EnKF)
P  L P Li1  Q
b
i
a
i1 i1
T
*
 Although the dimension of Pi f is huge, the rank
( ) <<
Pi nf (dominated by the errors of the day)
1 m f
P   (xk  x t )(xkf  x t )T
m k 1
m 
Ideally
b
i
 Using ensemble method to estimate
Physically,
“errors of day” are the instabilities of
the background flow. Strong instabilities
have a few dominant shapes
(perturbations lie in a low-dimensional
subspace).
 It makes sense to assume that large
errors are in similarly low-dimensional
spaces that can be represented by a
low order EnKF.
1 K f
b
Pi 
(xk  x f )(xkf  x f )T

K  1 k 1
1

X b  X bT
K 1
*
K ensemble members, K<<n
 Problem left: How to update ensemble ?
a
i.e.: How to get xi for each ensemble
member? Using
K times?
**
Ensemble Update: two approaches
1. Perturbed Observations method:
An “ensemble of data assimilations”
 It has been proven that an observational
ensemble is required (e.g., Burgers et al.
Pia nn  [I  Ki H]Pib
1998). Otherwise
is not satisfied.
 Random perturbations are added to the
observations to obtain observations for each
independent cycle
y
o
i (k )
 y  noise
o
i
 However, it introduces a source of sampling
errors when perturbing observations.
a
xbi (k)  Mxi1
(k)
1 K b
P 
(xk  x b )(xkb  x b )T

K  1 k 1
b
i
Ki  PibHT [HPib HT  R]1
xia(k)  xbi (k)  Ki (yio(k)  Hxbi (k) )
Ensemble Update: two approaches
2. Ensemble square root filter
(EnSRF)
 Observations are assimilated to
update only the ensemble mean.
xia  xbi  Ki (yio  H xbi )
 Assume analysis ensemble
perturbations can be formed by
transforming the forecast ensemble
perturbations through a transform
matrix
xbi  Mxia1
1 K b
P 
(xk  x b )(xkb  x b )T

K  1 k 1
b
i
Ki  PibHT [HPib HT  R]1
xbi  xbi  Ki (yio  H xbi )
Xia  Ti Xbi
xia  xia  Xa
T
1
1
X a X a  Pia n  n  [I  Ki H]Pib  [I  Ki H]
X b X bT
k 1
k 1
Xia  Ti Xbi
Several choices of the transform matrix
 EnSRF, Andrews 1968, Whitaker and Hamill, 2002)
% b, K
%  K
Xa  (I  KH)X
 EAKF (Anderson 2001)
X a  AX b
 ETKF (Bishop et al. 2001)
X a  Xb T
 LETKF (Hunt, 2005)
Based on ETKF but perform analysis
simultaneously in a local volume
surrounding each grid point
one observation
at a time
An Example of the analysis corrections from
3D-Var (Kalnay, 2004)
An Example of the analysis corrections from
EnKF (Kalnay, 2004)
Summary of LETKF (Hunt, 2005)
Forecast step (done globally): advance the ensemble 6 hours
xb(i)
 M(xa(i)
i  1,..., k
n
n1 )
Analysis step (done locally). Local model dimension m, locally used s obs
X  xb(1) | ... | xb(k) 
Xb  X  x b
Matrix of ensembles (mxk)
Matrix of ensemble perturbations (mxk)
Yb  H(X)  H(xb )  HXb
Pb  (k 1)1 XbT Xb
Matrix of ens. obs. perturbations (sxk)
b
in model space, so P
(Pa )1  (Pb )1  HT R1H
 (k  1)1 I
in ensemble space (kxk)
in model space, so that
%a )1  (k 1)I  (HXb )T R1 (HXb )
(P
in ensemble space (kxk)
Summary of LETKF (cont.)
Forecast step (done globally): advance the ensemble 6 hours
xb(i)
 M(xa(i)
i  1,..., k
n
n1 )
Analysis step (done locally). Local model dimension m, locally used obs s
Xb  X  x b
Yb  H(X)  H(xb )  HXb
1
a
b T
1
b
%
(P )  (k  1)I  (HX ) R (HX ) 
in ensemble space (kxk)
%a Xb
Pa  X aT Xa  XbT P
%a )1/2
Xa  Xb (P
Ensemble analysis perturbations in model space (mxk)
%a YbT R1 (yo  yb )
wan  P
n
n
xna  xna  Xb wan
in model space (mxm)
Analysis increments in ensemble space (kx1)
Analysis in model space (mx1)
We finally gather all the analysis and analysis perturbations from each grid point
and construct the new global analysis ensemble (nxk)
Summary steps of LETKF
1) Global 6 hr ensemble forecast
starting from the analysis ensemble
Obs.
operator
2) Choose the observations used for
each grid point
3) Compute the matrices of forecast
perturbations in ensemble space Xb
Ensemble forecast
at obs. locations
4) Compute the matrices of forecast
perturbations in observation space Yb
Pb
5) Compute
in ensemble space
space and its symmetric square root
6) Compute wa, the k vector of
perturbation weights
7) Compute the local grid point
analysis and analysis perturbations.
8) Gather the new global analysis
ensemble.
LETKF
Ensemble
forecast
Ensemble analysis
GCM
Observations
EnKF vs 4D-Var
 EnKF is simple and model
independent, while 4D-Var requires the
development and maintenance of the
adjoint model (model dependent)
Obs.
operator
Ensemble forecast
at obs. locations
4D-Var can assimilate asynchronous
observations, while EnKF assimilate
observations at the synoptic time.
 Using the weights wa at any time 4DLETKF can assimilate asynchronous
observations and move them forward or
backward to the analysis time
Disadvantage of EnKF:
Low dimensionality of the ensemble in
EnKF introduces sampling errors in the
estimation of Pb . ‘Covariance
localization’ can solve this problem.
LETKF
Ensemble
forecast
Ensemble analysis
GCM
Observations