Transcript Ensemble methods in data assimilation: an introduction

Introduction to
Data Assimilation
Peter Jan van Leeuwen
IMAU
Basic estimation theory
E{e0} = 0
E{em} = 0
E{e02} = s02
E{em2} = sm2
E{e0em} = 0
T0 = T + e0
Tm = T + em
Assume a linear best estimate: Tn = a T0 + b Tm
with Tn = T + en
Find a and b such that:
E{en} = 0
b=1-a
E{en2} minimal
2
s
m
a = __________
s02 + sm2
Basic estimation theory
2
2
s
s
m
0
Solution: Tn = _______ T0 + _______ Tm
s02 + sm2
s02 + sm2
1 = ___
1 + ___
1
___
and
sn 2
s02 sm2
Note: sn smaller than s0 and sm !
Best Linear Unbiased Estimate BLUE
Can we generalize this?
• More dimensions
• Nonlinear estimates (why linear?)
• Observations that are not directly
modeled
• Biases
The basics: probability density functions
P(u)
0.5
1.0
u (m/s)
The model pdf
P[u(x1),u(x2),T(x3),..
u(x1)
u(x2)
T(x3)
Observations
• In situ observations: e.g.
sparse hydrographic
observations, irregular in
space and time
• Satellite observations: e.g. of
the sea-surface
QuickTime™ and a
TIFF (Uncompressed) decompressor
are needed to see this picture.
Data assimilation: general formulation
NO INVERSION !!!
Bayes’ Theorem
Conditional pdf:
Similarly:
Combine:
Even better:
Filters and smoothers
Filter: solve 3D
problem several
times
Smoother: solve 4D
problem once
Time
Note: the model is (highly) nonlinear!
Pdf evolution in time
Model equation:
Pdf evolution: Kolmogorov’s equation
(Fokker-Planck equation)
(Ensemble) Kalman Filter
Only consider mean and covariance
At observation times:
- Assume p(d|) and p() are Gaussian, and use Bayes
-The mean of the product of 2 Gaussians is equal to linear
combination of the 2 means: E{|d} = a E{} + b E{d|}
- But we have seen this before in the first example !
(Ensemble) Kalman Filter II
2
2
s
s
m
0
Old solution: Tn = _______ T0 + _______ Tm
s02 + sm2
s02 + sm2
But now for covariance matrices:
mnew = R (P+R)-1 mold + P (P+R)-1d
Kalman Filter notation:
with Kalman gain
mnew = mold + K (d - H mold)
K = PHT (HPHT + R)-1
The error covariance:
tells us how model variables co-vary
For example SSH at point x with SSH at point y:
PSSH SSH(x,y) = E{ (SSH(x) - E{SSH(x)}) (SSH(y) - E{SSH(y)}) }
Or SSH at point x and SST at point y:
PSSH SST(x,y) = E{ (SSH(x) - E{SSH(x)}) (SST(y) - E{SST(y)}) }
Spatial
correlation of
SSH
and SST in
the Indian
Ocean
x
x
Haugen and Evensen, 2002
Covariances
between
model
variables
Haugen and Evensen, 2002
Summary on Kalman filters:
• Gaussian pdf’s for model and observations
• Propagation of error covariance P
If N operations for state vector evolution,
then N2 operations for P evolution…
Problems:
• Nonlinear dynamics, so non-Gaussian statistics
• Evolution equation for P not closed
• Size of P (> 1,000,000,000,000) ….
Propagation of pdf in time:
ensemble or particle methods
Example of Ensemble Kalman Filter (EnKF)
MICOM model with
1.3 million model variables
Observations:
Altimetry, infra-red
Validated with
hydrographic
observations
SST (-2K to +2K)
SSH (-10 cm to +10 cm)
RMS difference with XBT-data
Spurious covariances
?
Localization in EnKF-like methods
EnKF:
with Kalman gain
Local updating: restrict update using only
local covariances:
Schurproduct, or direct cut-off
Ensemble Kalman Smoother (EnKS)
Basic idea: use covariances over time.
Efficient implementation:
1) run EnKF, store ensemble at observation times
2) add influence of data back in time using
covariances at different times
Nonlinear filters
10
Probability density
function of layer
thickness
of first layer
at day 41
during
data-assimilation
No Kalman filter
No variational
methods
8
6
4
2
0
408 412 416 420 424 428 432 436 440 444 448 452 456
The particle filter
(Sequential Importance Resampling SIR)
Ensemble
with
Particle filter
SIR-results for a quasi-geostrophic ocean
model around South Africa with 512 members
Smoothers: formulation
Model error
Initial error
Observation error
Boundary errors etc. etc.
Smoothers: prior pdf
Smoothers: posterior pdf
Assume all errors are Gaussian:
model
initial
observation
Smoothers in practice:
Variational methods
Assume Gaussian pdf for model errors and observations:
in which
model dynamics
initial condition
Find min J from variational derivative:
J is costfunction or penalty function
model-obs misfit
Gradient descent methods
J
1’
3
4
6 5
model variable
2
1
The Euler-Lagrange equations
Forward integrations
Backward integrations
Nonlinear two-point boundary value problem
solved by linearization and iteration
4D-VAR strong constraint
Assume model errors negligible:
In practice only a few linear and one or two nonlinear
iterations are done….
No error estimate (Hessian too expensive and unwanted…)
Example 4D-VAR: GECCO
•
1952 through 2001 on a 1º global grid with 23 layers in the vertical,
using the ECCO/MIT adjoint technology.
•
Model started from Levitus and NCEP forcing and uses state of the
art physics modules (GM, KPP).
•
Control parameters: initial temperature and salinity fields
and the time varying surface forcing,
The Mean Ocean Circulation, global
Residual values can reveal inconsistencies in
data sets (here geoid).
MOC at 25N
Bryden et al.
(2005)
Error estimates
J
X
Local curvature from
second derivative of J,
the Hessian
Other smoothers
Representers, PSAS, Ensemble Kalman smoother, ….
Simulated annealing (Metropolis Hastings), …
Relations between model
variables
• Covariance gives linear correlations between variables
• Adjoint gives linear correlation between variables along a
nonlinear model run (linear sensitivity)
• Pdf gives full nonlinear relation between variables (nonlinear
sensitivity)
Parameter estimation
Bayes:
Looks simple, but we don’t observe model parameters….
We observe model fields, so:
in which Hhas to be found from model integrations
Example: ecosystem modeling
29 parameters
of which
15 were estimated
and 14 were kept
Fixed.
Estimated
parameters
from particle
filter (SIR)
All other methods that
were tried, including
4D-VAR and EnKF failed.
Losa et al, 2001
Estimate size of model error
Brasseur et al, 2006
Why data assimilation?
• Forecasts
• Process studies
• Model improvements
- model parameters
- parameterizations
• ‘Intelligent monitoring’
Conclusions
• Evolution of pdf with time is essential ingredient
• Filters: dominated by Kalman-like methods, but
moving towards nonlinear methods (SIR etc.)
• Smoothers: dominated by 4D-VAR,
New ideas needed!