1. No category

“What is data assimilation really solving?
What can it solve?
How is the calculation actually done?”
Ichiro Fukumori (JPL)
GODAE Summer School
September, 2004
What are we really solving?
How can we solve it?
OUTLINE
1. What is assimilation really solving?
a) Filtering vs Smoothing
b) Observability & Controllability
c) Data Error and Model Error
2. Examples and Practical Issues; e.g., ECCO
a) A hierarchical approach
b) Practical derivation and implementation
c) Assessment
What is assimilation really solving?
Filtering vs Smoothing
1. Kalman filtering and recursive smoothing are both
least-squares inversions, but solve different parts
of the assimilation problem,
2. Filtered solutions are physically inconsistent,
3. Smoothed solutions can be physically consistent,
4. Estimating model error sources is fundamental.
I.Fukumori/JPL
The Mathematical Problem of Assimilation
Different times
Observation Eqs
M
⎛
⎞ ⎛ M ⎞
Assimilation problem: Given
⎜
⎟ ⎜ y ⎟
Hx t
observations, y, and a model (A,
⎜
⎟ ⎜ t ⎟
⎜
⎟ ⎜ y t +1 ⎟
Hxt +1
G, H), what are the model states
⎜ Eqs
⎟ ⎜
⎟
Model
≈
M
M
(x) and their controls (u)?
⎜
⎟ ⎜
⎟
⎜ xt +1 − Axt − Gut ⎟ ⎜ 0 ⎟
⎜
⎟ ⎜
⎟
wind
⎜ xt + 2 − Axt +1 − Gut +1 ⎟ ⎜ 0 ⎟
⎛ M ⎞
⎜
⎟ ⎜ M ⎟
M
⎜ u
⎟ surface heat flux
⎝
⎠ ⎝
⎠
⎛ M ⎞
⎜ i, j,k ⎟
⎜τ ⎟
⎜ ui +1, j , k ⎟ sea⎛ level M
⎞
i, j ⎟
⎜
⎜
⎟
⎜
⎟
⎜ τ i +1, j ⎟
M
x
y
η
,
x: model state
(
)
⎧
1
1
⎜
⎟
⎜
⎟
⎜
⎟
⎪
⎜ vi , j , k ⎟
⎜ η ( x2 , y2 ) ⎟
M
y: observations
⎟
⎪
x=⎜
⎟, y = ⎜
⎟ , u = ⎜⎜
e.g.,
qi , j ⎟
M
⎪H: observation operator
⎜ M ⎟
⎜
⎟
⎜
⎟
⎨
⎜ Ti , j , k ⎟
⎜ T ( xn , yn , zn ) ⎟
M
A
G
,
:
model
dynamics
⎜
⎟
⎪
⎜
⎟
⎜⎜
⎟⎟
⎜ ωi , j , k ⎟
⎪u: control (e.g., forcing,
M
⎜ M ⎟
⎝
⎠
⎜⎜
⎟⎟
temperature
⎪
⎜
⎟
M ⎠
si , j , k
⎩ and model errors)
⎝
⎜
⎟
interior flux error
⎜ M ⎟
I.Fukumori/JPL
⎝
⎠
The Mathematical Problem of Assimilation
M
⎛
⎞ ⎛ M ⎞
⎜
⎟ ⎜ y ⎟
Hx t
⎜
⎟ ⎜ t ⎟
⎜
⎟ ⎜ y t +1 ⎟
Hxt +1
⎜
⎟ ⎜
⎟
≈
M
M
⎜
⎟ ⎜
⎟
⎜ xt +1 − Axt − Gut ⎟ ⎜ 0 ⎟
⎜
⎟ ⎜
⎟
⎜ xt + 2 − Axt +1 − Gut +1 ⎟ ⎜ 0 ⎟
⎜
⎟ ⎜ M ⎟
M
⎝
⎠ ⎝
⎠
⎛O
⎜L
⎜
⎜L
⎜
⎜L
⎜L
⎜
⎜L
⎜N
⎝
This problem is essentially
an inverse problem
Ea ≈ b
⎛ M ⎞
M
M
M
M
M N⎞ ⎜
⎟ ⎛ M ⎞
x
t ⎟
⎜ y ⎟
H
L⎟ ⎜
⎟ ⎜ xt +1 ⎟ ⎜ t ⎟
H
L⎟ ⎜
⎟ ⎜ y t +1 ⎟
⎟
⎟ ⎜ xt + 2 ⎟ ⎜
≈⎜ M ⎟
L L L L L L⎟
⎜ M ⎟
−A I
−G
L⎟ ⎜
⎟ ⎜ 0 ⎟
⎟
⎟ ⎜ ut ⎟ ⎜
0
−A I
−G L ⎟
⎟
⎜ ut +1 ⎟ ⎜
M
M
M
M
M O ⎟⎠ ⎜⎜
⎟⎟ ⎜⎝ M ⎟⎠
I.Fukumori/JPL
⎝ M ⎠
Least-Squares Solution
Almost all direct assimilation methods can be
recognized as a form of least-squares.
Ea ≈ b
For a linear inverse problem
The least-squares solution is given by
(Gauss-Markov inversion)
aˆ = a 0 + R aa E
T
( ER
R aˆ aˆ = R aa − R aa E
T
E + R bb )
T
aa
( ER
( b − Ea 0 )
E + R bb ) ER aa
T
aa
−1
−1
aˆ : least-squares solution
a 0 , R aa : first guess (prior solution) and its error covariance matrix
R bb : error covariance matrix of b
ˆ
R aa
ˆ ˆ : Error covariance of a
I.Fukumori/JPL
Kalman filter as Least-Squares
Linear inverse problem and its least-squares solution;
aˆ = a 0 + R aa E
T
( ER
E + R bb )
T
aa
−1
( b − Ea 0 )
Ea ≈ b
R aˆ aˆ = R aa − R aa ET ( ER aa ET + R bb ) ER aa
−1
Kalman filter estimate :
a
f
ˆ
ˆ
x
=
A
x
t −1 + Gu t −1
t
“forecast”
“analysis”
xˆ ta = xˆ tf + Pt f HT ( HPt f HT + R t )
−1
(y
t
− Hxˆ tf
)
Pt f = APˆ ta−1 AT + GQG T
ˆP a = P f − P f HT ( HP f HT + R )−1 HP f
t
t
t
t
t
t
There is a direct correspondence between the Kalman filtering
algorithm and least-squares; i.e., Kalman filtering is nothing more
than an inversion of the observation equation/operator H .
Inverse Solution
−1
T
T
aˆ = a 0 + R aa E ( ER aa E + R bb ) ( b − Ea 0 )
is the linear inverse solution of Ea ≈ b for any of following;
T
1. Minimum (expected error) variance: min diag ( a − a )( a − a )
if
T
T
True a
R aa = ( a0 − a )( a0 − a ) , R bb = ( b − Εa )( b − Εa )
)
(
(
T
2. Least-squares: min ( a − a0 ) Wa ( a − a0 ) + ( b − Ea) Wb ( b − Ea)
−1
−1
if Wa = R aa
, Wb = R bb
T
3. Maximum likelihood: max ⎡⎣℘( a ) ⎤⎦
if Gaussian probability distribution
T
−1
℘( a ) = exp − ( a − a0 ) R aa
( a − a0 ) − (b − Ea)T R bb−1 (b − Ea)
(
)
)
Least-squares does not necessarily assume Gaussian statistics.
I.Fukumori/JPL
Kalman filter as Least-Squares
M
⎛
⎞ ⎛ M ⎞
⎜
⎟ ⎜ y ⎟
Hx t
⎜
⎟ ⎜ t ⎟
⎜
⎟ ⎜ y t +1 ⎟
Hxt +1
⎜
⎟ ⎜
⎟
≈
M
M
⎜
⎟ ⎜
⎟
⎜ xt +1 − Axt − Gut ⎟ ⎜ 0 ⎟
⎜
⎟ ⎜
⎟
−
−
x
A
x
G
u
0
t +1
t +1 ⎟
⎜ t +2
⎜
⎟
⎜
⎟ ⎜ M ⎟
M
⎝
⎠ ⎝
⎠
M
⎛
⎞ ⎛ M ⎞
⎜
⎟ ⎜y ⎟
Hxt −1
⎜
⎟ ⎜ t −1 ⎟
⎜
⎟ ⎜ yt ⎟
Hx t
⎜
⎟ ⎜
⎟
≈
A
G
0
x
x
u
−
−
0
0 ⎟
⎜ 1
⎜
⎟
⎜
⎟ ⎜ M ⎟
M
⎜
⎟ ⎜
⎟
A
G
u
0
x
x
−
−
t
t ⎟
⎜ t +1
⎜
⎟
⎜ x − Ax − Gu ⎟ ⎜ 0 ⎟
t −1
t −1 ⎠
⎝ t
⎝
⎠
xˆ = xˆ + Pt H
a
t
f
t
f
T
( HP
t
f
H + Rt )
T
−1
(y
t
− Hxˆ tf
)
• This is an inversion of the
observation equation but NOT an
inversion of the model equations.
y The estimate is the least-squares
solution for the end state, but not
at intermediate times.
y The virtue of Kalman filtering is
that these inverted end states do
not have to be recomputed over
the entire period every time new
observations are obtained.
Kalman filter as Least-Squares
The temporal evolution of filtered estimates are physically
inconsistent (e.g., budgets do not close).
Kalman filter estimate : xˆ tf = Axˆ ta−1 + Gut −1
xˆ = xˆ + Pt H
a
t
f
t
f
T
( HP
t
f
H + Rt )
T
−1
(y
t
− Hxˆ tf
Combining the two steps give,
xˆ = Axˆ
a
t
a
t −1
+ Gut −1 + Pt H
f
T
( HP
t
f
H + Rt )
T
−1
(y
t
− Hxˆ tf
)
Not associated with any particular
physics (e.g., advection, mixing,
external forcing) embodied in A and G.
)
Consistency of Temporal Evolution
xˆ = Axˆ
a
t
a
t −1
+ Gut −1 + Pt H
f
T
( HPt H + R t )
f
T
−1
f
ˆ
y
−
H
x
( t
t
)
Because of the filter’s data increments, the temporal evolution
of the filtered state from xˆ ta−1 to xˆ ta , etc, cannot be physically
accounted for; e.g., budgets cannot be closed.
Filtered
Estimate
xˆ tf = Axˆ ta−1 + Gut −1
xˆ ta = xˆ tf + Pt f HT ( HPt f HT + R t )
xˆ
a
t
−1
(y
t
− Hxˆ tf
)
yt
xˆ ta−1
Smoothed
Estimate
xˆ tf
time
I.Fukumori/JPL
Example: Inconsistencies in
Atmospheric Analyses
24% of mass change in NCEP’s operational
analysis is physically unaccounted for.
1.0
change over 6-hours
assimilation increment
0.8
0.6
0
10
20
mbar
0
2
4
mbar
Smother as Least-Squares
Kalman filter estimate:
xˆ = Axˆ
a
t
a
t −1
+ Gut −1 + Pt H
f
T
( HP
t
f
H + Rt )
T
−1
(y
t
− Hxˆ tf
)
a
ˆ
x
Smoothing problem: Given the updated estimate t invert the
model to correct the prior state x t −1 and control ut −1 that;
s
ˆ
⎛
x
a
s
s
t −1 ⎞
xˆ t = Axˆ t −1 + Guˆ t −1 = ( A G ) ⎜ s ⎟
⎝ uˆ t −1 ⎠
A least-squares solution that exactly satisfies the inverse
problem Ea = b can be obtained by setting R bb = 0 in
the canonical least-squares solution;
By substitution, we have;
⎛ Pa A T ( APa A T + GQ G T ) −1 ⎞
a
s
t −1
t −1
⎛ xˆ t −1 ⎞ ⎛ xˆ t −1 ⎞ ⎜ t −1
⎟ ( xˆ a − Axˆ a − Gu )
=
+
t
t −1
t −1
⎜ s ⎟ ⎜
⎟ ⎜
T
T
T −1 ⎟
a
ˆ
u
u
⎝ t −1 ⎠ ⎝ t −1 ⎠ Qt −1G ( GQt −1G + APt −1A )
⎝
⎠
Smother as Least-Squares
More generally,
−1
T
T
T
a
a
⎛
⎞
⎛ xˆ ts ⎞ ⎛ xˆ ta ⎞ ⎜ Pt A ( APt A + GQt G ) ⎟ s
a
ˆ
ˆ
=
+
x
−
A
x
(
t +1
t − Gu t )
⎜ s⎟ ⎜ ⎟ ⎜
1
−
⎝ uˆ t ⎠ ⎝ ut ⎠ Qt G T ( GQt G T + APta A T ) ⎟
⎝
⎠
⎛ Pa A T ( APa A T + GQ G T ) −1 APa ⎞ ⎛ a s a T ⎞
Lt Pt +1Lt
t
t
t
⎛ P ⎞ ⎛P ⎞ ⎜ t
⎟
⎟
+⎜
⎜ ⎟=⎜ ⎟−⎜
T
1
−
⎝ Q ⎠ ⎝ Qt ⎠ Qt GT ( GQt GT + APta A T ) GQt ⎟ ⎜⎝ M ta Pts+1M ta ⎟⎠
⎝
⎠
a
t
s
t
s
t
where,
( AP A + GQ G )
≡ Q G ( GQ G + AP A )
L ≡P A
a
t
M
a
t
a
t
a
t
T
Least-squares solution;
Ea ≈ b
Extra term due
to errors in xˆ ts+1
t
T
t
T −1
T
T
t
T −1
a
t
aˆ = a 0 + R aa E
T
( ER
E + R bb )
T
aa
−1
( b − Ea 0 )
R aˆ aˆ = R aa − R aa ET ( ER aa ET + R bb ) ER aa
−1
Smoother as Least-Squares
−1
T
T
T
a
a
⎛
⎞
⎛ xˆ ts ⎞ ⎛ xˆ ta ⎞ ⎜ Pt A ( APt A + GQt G ) ⎟ s
a
ˆ
ˆ
=
+
x
−
A
x
(
t +1
t − Gu t )
⎜ s⎟ ⎜ ⎟ ⎜
1
−
⎝ uˆ t ⎠ ⎝ ut ⎠ Qt G T ( GQt G T + APta A T ) ⎟
⎝
⎠
y This is the Rauch-Tung-Striebel (RTS) fixed-interval
smoother and is an inversion of the model equations.
M
⎛
⎞ ⎛ M ⎞
⎜
⎟ ⎜ y ⎟
Hx t
⎜
⎟ ⎜ t ⎟
⎜
⎟ ⎜ y t +1 ⎟
Hxt +1
⎜
⎟ ⎜
⎟
≈
M
M
⎜
⎟ ⎜
⎟
⎜ xt +1 − Axt − Gut ⎟ ⎜ 0 ⎟
⎜
⎟ ⎜
⎟
x
−
A
x
−
G
u
0
t +1
t +1 ⎟
⎜ t +2
⎜
⎟
⎜
⎟
⎜
⎟
M
M
⎝
⎠ ⎝
⎠
y This solves the entire assimilation
problem by sequentially solving
smaller calculations.
y Another virtue is its results’
physical consistency.
Smoother as Least-Squares
⎛ P a A T ( AP a A T + GQ G T )−1 ⎞
t
t
⎛ xˆ ⎞ ⎛ xˆ ⎞ ⎜ t
⎟ xˆ s − Axˆ a − Gu
=
+
( t +1 t
⎜ ⎟ ⎜ ⎟ ⎜
t)
−1 ⎟
a T
T
T
ˆ
u
u
⎝ ⎠ ⎝ t ⎠ ⎜ Qt G ( GQt G + APt A ) ⎟
⎝
⎠
s
t
s
t
a
t
y By construction, this satisfies the model equation,
xˆ ts = Axˆ ts−1 + Guˆ ts−1
making it amenable to process
studies. This should not be confused
with the “strong constraint” solution;
In fact, xˆ ts ≠ Axˆ ts−1 + Gut −1
y Control estimates are of the essence. Smoothed state
estimates can be computed from the control estimates.
y The model for model process noise u t must be physically
based for the result to be physically sensible.
I.Fukumori/JPL
Smoother as Least-Squares
⎛ P a A T ( AP a A T + GQ G T )−1 ⎞
t
t
⎛ xˆ ⎞ ⎛ xˆ ⎞ ⎜ t
⎟ xˆ s − Axˆ a − Gu
=
+
( t +1 t
⎜ ⎟ ⎜ ⎟ ⎜
t)
−1 ⎟
a T
T
T
ˆ
u
u
⎝ ⎠ ⎝ t ⎠ ⎜ Qt G ( GQt G + APt A ) ⎟
⎝
⎠
s
t
s
t
a
t
y The RTS smoother recursively carries information of
formally future observations backwards in time.
y Utilizes all data within a fixed time-interval whereas the
Kalman filter only utilizes data up to that instant.
Data
Kalman filter
time
Smoother
y More accurate than the filter’s estimate.
Other Smoothers
y Fixed-Interval Smoother;
time
Besides the recursive RTS smoother, the adjoint method
and Green’s function optimizations are other commonly
used fixed-interval smoothers.
y Fixed-Point Smoother;
time
y Fixed-Lag Smoother;
time
What are we really solving?
How can we solve it?
OUTLINE
1. What is assimilation really solving?
a) Filtering vs Smoothing
b) Observability & Controllability
c) Data Error and Model Error
2. Examples and Practical Issues; e.g., ECCO
a) A hierarchical approach
b) Practical derivation and implementation
c) Assessment
Observability & Controllability
Kalman Filter is an inversion of observation H.
RTS smoother is an inversion of model evolution A and G.
Identifies the problem the filter & smoother are
solving (inverting).
Observability identifies if the model state that can be determined.
Controllability identifies if controls can be determined.
Identifies the solution that could be determined.
Observability
If there are no errors in the observations and no errors in the
model (except state is unknown), can we determine the
model state uniquely ?
y1 = H1x1 = H1 ( Ax 0 − Gu0 )
y1 + H1Gu0 = H 2 Ax 0
M
)
y t + H t Gut −1 + Ht AGut − 2 + L + Ht A t −1Gu0 ≡ y t = H t A t x 0
)
y
⎛ 1⎞
⎛ H1A ⎞
⎜ M ⎟ = Ox where O ≡ ⎜ M ⎟
1
⎜) ⎟
⎟
⎜
⎜y ⎟
⎜ H At ⎟
⎝ t⎠
⎝ t ⎠
If O is full rank (column), there is a unique solution for
the state, and the system is said to be observable.
Controllability
Can we determine controls that can drive the model to an
arbitrary state?
x1 = Ax 0 + Gu0
x 2 = Ax1 + Gu1 = A 2 x 0 + AGu0 + Gu1
M
x t = Ax t −1 + Gut −1 = A t x 0 + A t −1Gu0 + A t − 2 Gu1 L + Gut −1
⎛ u0 ⎞
⎜ u ⎟
x t − A t x 0 = C ⎜ 1 ⎟ where
⎜ M ⎟
⎜u ⎟
⎝ t −1 ⎠
C ≡ ( G AG L A t −1G )
If C is full rank (row), there is a solution for the control that
satisfies this equation, and the system is said to be controllable.
Observability & Controllability
These conditions are conceptually important. But
the physical state being observable does not mean
any meaningful estimate can be made; i.e.,
a
f
Difference between P and P may be
indistinguishable.
What are we really solving?
How can we solve it?
OUTLINE
1. What is assimilation really solving?
a) Filtering vs Smoothing
b) Observability & Controllability
c) Data Error and Model Error
2. Examples and Practical Issues; e.g., ECCO
a) A hierarchical approach
b) Practical derivation and implementation
c) Assessment
What is assimilation really solving?
“Data Error” and “Model Error”
1. The two are better thought of as observation
constraint error and model constraint error,
respectively, as opposed to errors of the data and
model, per se.
2. “Data Error” includes representation errors, i.e.,
things that the models do not represent,
3. How to specify the two prior errors.
“Data Error” & “Model Error”
M
⎛
⎞ ⎛ M ⎞
⎜
⎟ ⎜ y ⎟
Hx t
⎜
⎟ ⎜ t ⎟
⎜
⎟ ⎜ y t +1 ⎟
Hxt +1
⎜
⎟ ⎜
⎟
≈
M
M
⎜
⎟ ⎜
⎟
⎜ xt +1 − Axt − Gut ⎟ ⎜ 0 ⎟
⎜
⎟ ⎜
⎟
x
−
A
x
−
G
u
0
t +1
t +1 ⎟
⎜ t +2
⎜
⎟
⎜
⎟
⎜
⎟
M
M
⎝
⎠ ⎝
⎠
What are error covariance
matrices R and Q?
These two in effect define the
inverse solution, and therefore,
their understanding and proper
specification are fundamental to
solving the assimilation problem.
Kalman filter estimate :
“forecast”
“analysis”
xˆ tf = Axˆ ta−1 + Gut −1
xˆ ta = xˆ tf + Pt f HT ( HPt f HT + R t )
−1
(y
t
− Hxˆ tf
Pt f = APˆ ta−1 AT + GQG T
ˆP a = P f − P f HT ( HP f HT + R )−1 HP f
t
t
t
t
t
t
)
“Data Error”
R represents the error of the data constraint and
includes model representation error.
The true model state xt is a finite dimensional
representation of the true state of the ocean, w t
xt ≡ Πw t
Π : Function defining the model state.
Observations are related to the true state of the ocean by
E: Function defining observation.
y t = Ew t + ε
ε: Measurement error
In terms of the true model state, y t = Hxt + {Ew t − HΠw t } + ε
(after Cohn, 1997, J. Met. Soc. Jap.)
Error of observation equation y t = Hx t
Ew t − HΠw t : Representation error
“Data Error”
Ew t − HΠw t : Representation error
Error-free
data
Error-free model
equivalent of data
Π : Defines the model state.
E : Defines observation.
H : Defines model equivalent
of data
Examples;
• Meso-scale variability for a non-eddy resolving model,
• External tides for a rigid lid model,
• Skin temperature for model with thick surface layer,
• Micro-structure for most models,
• Barotropic variability for a reduced-gravity model,
• Baroclinic variability for a barotropic model,
• Meteorologist do not force models too close to
observations in numerical weather forecasting,
• Individual float & drifter trajectories for any finite difference
model.
Examples of Drifter Trajectory
(Paduan and Niiler, 1993, JPO)
Example of diverging float
Typical “spaghetti” diagram
trajectories
“Model Error”
Q represents the error of the model constraint.
Model evolution: x t +1 = A ( x t , uˆ )
Ocean evolution: w t +1 = L ( w t , v t )
0
t
Forcing, boundary condition, mixing
parameters, model error sources.
Real forcing & boundary condition.
Functions describing the evolution
The true model evolution can be written as,
xt +1 = Πw t +1 = ΠL ( w t , v t )
{
}
xt +1 = A ( xt , uˆ 0t ) + ΠL ( w t , v t ) − A ( Πw t , uˆ 0t )
Error of the model evolution; x t +1 = A ( x t , uˆ 0t )
ΠL ( w t , v t ) − A ( Πw t , uˆ 0t ) : Model Process Noise
True evolution
Model evolution given
true model state
“Model Error”
ΠL ( w t , v t ) − A ( Πw t , uˆ 0t ) : Model Process Noise
Examples;
1) Differences in uˆ 0t and v (errors in non-state variables);
forcing, boundary condition, model parameters,
2) Errors due to differences in A and L (errors in model
algorithm); finite difference, interaction with scales and
processes ignored
{
}
xt +1 = A ( xt , uˆ 0t ) + ΠL ( w t , v t ) − A ( Πw t , uˆ 0t )
= A ( xt , ut )
Collectively, sources of model process noise are
modeled as elements of the model control vector.
Prior Error Specification
Covariance Matching (Fu et al., 1993, JPO)
Estimate R & Q from model-data comparison.
Data y and its model simulation equivalent m could be written as;
y =s+r
m=s+n
where, s : true signal ( Hx ) , r : data error,
m: simulation of y ( Hx sim ) ,
n : simulation error ( Hx sim − Hx )
Then,
yy T = ssT + sr T + rsT + rr T
mmT = ssT + snT + nsT + nnT
my T = ssT + sr T + nsT + nr T
Prior Error Specification
T
If, sn ≈ 0 etc,
R = rr T = yy T − my T
HPsim HT = nnT = mmT − my T
R estimate.
Implicit Q estimate.
Adjust process noise Q to match the solution of the
“Lyapunov Equation” to that of the <nnT> estimate.
T
T
Solve P as t→∞ Pt = APt −1A + GQG
Then choose Q such that
HPHT ≈ nnT
Prior Error Specification
Error variance estimate (Fukumori et al., 1999, JGR)
R = rrT = yy T − my T
HPsim HT = nnT = mmT − my T
(a) Data Constraint Error
(b) Model Simulation Error
80N
60N
40N
20N
0
20S
40S
60S
80S
80N
60N
40N
20N
0
20S
40S
60S
80S
0
90E
0
180
20
40
60
90W
80
100
0
0
2
(cm )
90E
-5
0
180
5
Q
10
15
90W
20
25
30
0
1st order
match
2
(cm )
HPsim HT
(c) Calibrated Model Process Noise
80N
60N
40N
20N
0
20S
40S
60S
80S
80N
60N
40N
20N
0
20S
40S
60S
80S
0
90E
0
180
1
2
90W
3
4
0
0
2
(dyn/cm2 )
90E
-5
0
180
5
10
15
90W
20
25
30
0
2
(cm )
I.Fukumori/JPL
What are we really solving?
How can we solve it?
OUTLINE
1. What is assimilation really solving?
a) Filtering vs Smoothing
b) Observability & Controllability
c) Data Error and Model Error
2. Examples and Practical Issues; e.g., ECCO
a) A hierarchical approach
b) Practical derivation and implementation
c) Assessment
“Estimating the Circulation
and Climate of the Ocean”
http://www.ecco-group.org
GOAL: Advance data assimilation into an operational
tool to help understand ocean circulation.
• 3 groups: JPL, Scripps Institution of Oceanography,
Massachusetts Institute of Technology.
• Employ as much observations as possible.
• State-of-the-art ocean general circulation model.
• Advanced assimilation (Kalman filter/smoother, adjoint,
Green’s function) characterized by its consistency.
I.Fukumori/JPL
ECCO Near Real-Time Analysis
• Nowcasts every 10-days.
• Highlights @
http://ecco.jpl.nasa.gov/external
• SSH, Temperature profiles & timemean sea level.
• Controls: winds, heat & fresh
water fluxes and mixing
parameters.
• Green’s function, Kalman
filter & RTS Smoother,
Adjoint method)
• LAS server from 1993 to
present: http://www.eccogroup.org/las
I.Fukumori/JPL
ECCO Routine Assimilation
MIT General Circulation Model (MITgcm)
(Marshall et al., 1997, JGR)
• Nonlinear primitive equation model
• Adaptable to massively parallel supercomputers
• Advanced physics (e.g., KPP & GM mixing)
• Global, high resolution
9 80°S~ 80°N
9 0.3°×10m in tropics
warm
Temperature Profile
0
120W 0N
depth (m)
100
200
180E 0N
300
400
cold
500
0
10
20
30
temperature (deg C)
I.Fukumori/JPL
Approximate Kalman filter & Smoother
Kalman filter estimate :
xˆ tf = Axˆ ta−1 + Gut −1
xˆ ta = xˆ tf + Pt f HT ( HPt f HT + R t )
−1
(y
t
− Hxˆ tf
)
Pt f = APˆ ta−1 AT + GQG T
ˆP a = P f − P f HT ( HP f HT + R )−1 HP f
t
t
t
t
t
t
The explicit evaluation of P (“Riccati Equation”) is impractical
because of their large dimension and time-evolution.
ECCO employs three approximations of P;
1. Time-asymptotic
2. Reduced-state
3. Partitioned
ECCO Approximate KFS
1) Time-asymptotic approximation (Fukumori et al., JPO, 1993);
Pt f
error
Pt ≈ P
f
f
Pf
time
P ≡ δ xδ x T ≈ B δ x ′δ x ′T BT ≈ BP′BT
3) Partitioning (Fukumori, MWR, 2002).
space
L
dim (δ x )
L
P ≈ ∑ Bi Pi′ BTi ,
i
dim (δ x i′ )
T
′
′
′
Pi ≡ δ x iδ x i
error
δ x ≈ B1δ x1′ + L + B Lδ x ′L ≈ ∑ Biδ x i′
i
error
2) State reduction (Fukumori and Rizzoli, JGR, 1995).
δ x ≈ Bδ x ′, dim(δ x ) dim(δ x ′)
space
I.Fukumori/JPL
What are we really solving?
How can we solve it?
OUTLINE
1. What is assimilation really solving?
a) Filtering vs Smoothing
b) Observability & Controllability
c) Data Error and Model Error
2. Examples and Practical Issues; e.g., ECCO
a) A hierarchical approach
b) Practical derivation and implementation
c) Assessment
ECCO Approximate KFS
1. Errors due to wind error (process noise) are estimated.
2. Partition global domain into 8 different regions (basins)
• Global barotropic cell
• 7 Regional baroclinic cells
I.
II.
III.
IV.
V.
VI.
VII.
Tropical Indian (>40 ºS)
Tropical Pacific (40 ºS ~ 40 ºN)
Tropical Atlantic (40 ºS ~ 40 ºN)
South Pacific (<20 ºS)
South Atlantic-Indian (<20 ºS)
North Pacific (>20 ºN)
North Atlantic (>20 ºN)
50
0
-50
0
100
200
300
I.Fukumori/JPL
ECCO Approximate KFS
3. State and control reduction: Horizontal velocity and
vertical displacement in terms of vertical dynamic
modes on a coarse horizontal grid. The reduced
control is wind error on the same coarse grid.
Vertical projection Horizontal interpolation
(objective mapping)
(dynamic modes)
Reduced state (amplitude of
dynamic mode on coarse grid)
δ u = Dvel Oau , δ v = Dvel Oav , δ η = Ddisp Oadisp , δ τ = Oaτ
∂T
∂S
δ T = δ η, δ S = δ η, δ h = F (δ η, T, S )
∂z
∂z
3a. Vertical projection: D vel , D disp
d ⎛ 1 dp ⎞ 1
⎜ 2
⎟+ 2 p=0
dz ⎝ N dz ⎠ c
2
d 2η ⎛ N ⎞
+⎜ ⎟ η =0
2
dz
⎝c⎠
I.Fukumori/JPL
ECCO Approximate KFS
3b. Horizontal Interpolation: O
“Objective Mapping” (Bretherton et al., 1976, DSR)
Given measurements ϕ , the underlying field θ is estimated by
A = ϕϕ T
θ = CA −1ϕ where C = θϕ T
This is just another
least-squares
aˆ = a 0 + R aa E
T
( ER
E + R bb )
T
aa
−1
( b − Ea 0 )
Horizontal interpolation, O, is defined as objective mapping
from a coarse grid system to the model grid system.
Example of mapping of a
perturbation at one coarse
grid point;
I.Fukumori/JPL
ECCO Approximate KFS
• Global barotropic cell: global 6º× 6º grid, barotropic UVH
(966 ×3=2898 variables)
• Regional baroclinic cells: regional 5º× 3º grid, baroclinic
UVD (5 gravest modes)
I.
II.
III.
IV.
V.
VI.
VII.
Tropical Indian (>40 ºS ; 308 ×3 ×5=4620)
Tropical Pacific (40 ºS ~ 40 ºN; 787 ×3 ×5=11805)
Tropical Atlantic (40 ºS ~ 40 ºN; 350 ×3 ×5=5250)
South Pacific (<20 ºS ; 633 ×3 ×5=9495)
South Atlantic-Indian (<20 ºS ; 664 ×3 ×5=9960)
North Pacific (>20 ºN; 271 ×3 ×5=4065)
North Atlantic (>20 ºN; 198 ×3 ×5=2970)
50
In comparison, the model’s state
dimension is 8 million (2 million
grid points).
0
-50
0
100
200
300
I.Fukumori/JPL
ECCO Approximate KFS
4. Compute P by formulating and explicitly
deriving partitioned reduced-state models.
4a. Formulate partitioned reduced-state models.
The reduced-state must form a closed system.
A state and control perturbation δx, δu can be written as
δ x = Bx ′ + Nn
: B, N : range & null space of reduced-state
ˆ ′ + Nm
ˆ
δ u = Bu
: The control equivalent of state reduction.
Given the model x t +1 = ℑ ( x t , ut )
the perturbation
satisfies δ x t +1 = ℑ ( x% + δ x t , u% + δ ut ) − ℑ ( x% , u% )
ˆ ′t + Nm
ˆ t − ℑ ( x% , u% )
or Bx ′t +1 + Nnt +1 = ℑ x% + Bx ′t + Nnt , u% + Bu
(
)
ECCO Approximate KFS
(
)
ˆ ′t + Nm
ˆ t − ℑ ( x% , u% )
Bx ′t +1 + Nnt +1 = ℑ x% + Bx ′t + Nnt , u% + Bu
Multiplying by the pseudo inverse of B, B* ( B*B = I , B* N = 0 )
ˆ ′t + Nm
ˆ t − ℑ ( x% , u% )
x ′t + = B* ℑ x% + Bx ′t + Nnt , u% + Bu
1
( (
)
)
For this to be a model for the reduced-state we cannot have
ˆ t ;
dependency on the null space, Nnt , Nm
( (
)
ˆ ′t − ℑ ( x% , u% )
x ′t +1 = B* ℑ x% + Bx ′t , u% + Bu
)
This would be achieved if null (range) space perturbations
remain in the null (range) space, because B* N = 0 ; i.e., state
reduction must form a closed-system for it to be effective.
ECCO Approximate KFS
4b. Derive partitioned reduced-state model in explicit
matrix form to facilitate error evaluation.
( (
)
ˆ ′t − ℑ ( x% , u% )
x ′t +1 = B* ℑ x% + Bx ′t , u% + Bu
)
x ′t +1 = A ′x ′t + G′u′t
A ′ and G′ can be derived as coarse grain Green’s functions;
( A ′)i = A ′ei + G′0 = B* ( ℑ ( x% + Bei , u% ) − ℑ ( x% , u% ) )
ˆ i ) − ℑ ( x% , u% ) )
( G′)i = A ′0 + G′ei = B* ( ℑ ( x% , u% + Be
Similarly, H′ can be derived by,
( H′)i = H′ei = H ( x% + Bei ) − H ( x% )
where H is a function defining model equivalent of data.
ECCO Approximate KFS
( A ′) = A ′e + G′0 = B ( ℑ ( x% + Be , u% ) − ℑ ( x% , u% ) )
*
i
i
i
Implementing B* requires some considerations;
1. O is sparse but O* is not. Thus, we do
O = ( O O ) OT
*
(O O)
T
−1
T
−1
is small and is, therefore, explicitly
computed and stored. OT is implemented as
the transpose (adjoint) of O .
2. η is evaluated by time-integrating vertical
velocity instead of inverting
∂T
∂S
δ T = δ η, δ S = δ η
∂z
∂z
which is singular where stratification is weak.
ECCO Approximate KFS
4c. Derive representative time-asymptotic limit of P.
• Stationary solution to the Riccati Equation requires a stationary
system, i.e., time-invariant A’, G’, H’, Q, R.
• Approximate with a 3-day assimilation cycle;
x ′t +1 = A ′x ′t + G′u′t
If ℑ ( x% , u% ) is a 1-day integration. Then A ′ = C′3 where
( C′)i = B* ( ℑ ( x% + Bei , u% ) − ℑ ( x% , u% ) )
• Approximate all observing operation H’ and R
with representative observations; e.g.,
Assume a full 10-day cycle of altimeter observations is
available every 3-days but with 3.3 times the error variance.
ECCO Approximate KFS
e.g., TOPEX/POSEIDON and Jason have a 3-day subcycle.
ECCO Approximate KFS
5. Implementation
a) Actual assimilation is conducted every 6-hours using all
observations taken within ±3-hours of the assimilation instant.
b) Filtering is performed using an alternate formulation of the
Kalman gain;
K-filter increment ∆xˆ = Pt H
a
t
filter gain
Pt H
f
f
T
T
( HP
t
f
H + Rt )
T
−1
(y
t
− H ( xˆ tf
( HPt H + R t ) = Pta HT R t
f
T
−1
))
−1
This form minimizes computations involving P and
simplifies the matrix inversion that is needed.
ECCO Approximate KFS
K-filter increment
∆xˆ ta = Pta HTt R t
−1
(
y t − H ( xˆ tf
))
c) The K-filter increment is sum of partitioned increments;
∆xˆ ta = ∑ Bi ∆xˆ ′t ,ai
i
where ∆xˆ ′ = Pi′a H′t ,i R t
a
t ,i
Or,
T
−1
(
∆xˆ ′t ,ai = Pi′a BT H t ,i R t
T
y t − H ( xˆ tf
−1
(
))
y t − H ( xˆ tf
))
using the adjoint of H in place of HT.
These calculations are performed from right to left;
There is no explicit computation of product Pi′a H′t ,i T R t −1
because it is computationally inefficient.
ECCO Approximate KFS
d) Anomalies are assimilated for sea level and temperature
profiles because of
i. Unknown marine geoid,
ii. Temporally uncorrelated process noise assumption.
(
(
a
T
−1
∆xˆ ′t a = Pi′ H′t ,i R t ( y t − y ) − H ( xˆ tf ) − m
))
where y and m are time-means of data and model
(simulation) equivalent of data mean, respectively.
This choice of reference amounts to assimilating temporal
anomalies of the data without altering the model time-mean.
e) Integrate original model after state and diagnostic variables
are updated; partitioned reduced-state model is never used
to integrate state.
ECCO Approximate KFS
f)
Smoothing is performed
i. Derive smoothed wind correction by,
∆uˆ = Qi G ( GQi G + AP A
s
i
T
= QG P
T
In partitioned
reduced-state;
a
i
T
f −1
) ( xˆ
T −1
s
i +1
− Axˆ ia − Gui )
s
a
ˆ
ˆ
∆
x
+
∆
x
( i +1 i +1 )
∆uˆ ts = ∑ Bˆ i ∆uˆ ′t ,si
i
where ∆uˆ ′ = Qi′G′i T Pi′ f
s
t ,i
−1
( ∆xˆ ′
s
t +1,i
+ ∆xˆ ′t +a1,i )
ii. Re-run model simulation with smoother corrected wind
because,
xˆ ta+1 + ∆xˆ ts+1 ≠ ℑ ( xˆ ta + ∆xˆ ts , ut + ∆uˆ ts )
due to approximations, but
xˆˆ s = ℑ xˆˆ s , u + ∆uˆ s
t +1
(
t
t
t
)
Skill & Test of Assimilation
Assimilation should make models closer to data,
both assimilated and independent.
Data not explained by model
Data variance explained by models; yy
T
−
( y − Hx )( y − Hx )
T
Filters (forecasts) and Smoothers
Simulations
I.Fukumori/JPL
Skill & Test of Assimilation
Model-data residuals should be comparable to expectation
simulation - forecast residual
1st order match
HPsim HT − HP f HT
forecast - analysis residual
HP f HT − HP a HT
I.Fukumori/JPL
Skill & Test of Assimilation
Comparison to independent data: ECCO assimilation explains
observed polar motion better than ocean simulation does.
Coherence of observed & modeled excitation
(Gross et al., 2003, JGR)
Summary
1. Data assimilation is an inverse
problem of observations & model
with state & model controls as
unknowns;
2. Data constraint error includes
model representation error,
M
⎛
⎞ ⎛ M ⎞
⎜
⎟ ⎜ y ⎟
Hxt
⎜
⎟ ⎜ t ⎟
⎜
⎟ ⎜ y t +1 ⎟
H x t +1
⎜
⎟ ⎜
⎟
≈
M
M
⎜
⎟ ⎜
⎟
⎜ x t +1 − A x t − G u t ⎟ ⎜ 0 ⎟
⎜
⎟ ⎜
⎟
x
A
x
G
u
0
−
−
t +1
t +1 ⎟
⎜ t+2
⎜
⎟
⎜
⎟
⎜
⎟
M
⎝
⎠ ⎝ M ⎠
3. Kalman filters and RTS smoothers are recursive leastsquares inversions,
4. Kalman filters invert observations; RTS smoothers invert
model,
Summary
5. Smoother’s control estimates are fundamental in solving
the assimilation problem,
filter
smoother
6. Temporal evolution of filters is not
physically consistent, but that of
smoothers is,
7. Various approximations & simplifications are possible
that make solving the assimilation problem tractable;
ECCO KFS.
Future Challenge
Modeling different model error sources by identifying
effective corresponding reduced-state approximations B is
a central issue in data assimilation.
Controls need explicit modeling for estimation;
xt +1 = A ( xt , uˆ 0t ) + ΠL ( w t , v t ) − A ( Πw t , uˆ 0t )
{
}
= Axt + Gut
True evolution
Model evolution given
true model state
What are the most effective basis functions B that would
represent effects of various model error sources ?
Effective = Small dimension but closed system.
)
( (
ˆ ′ ) − ℑ ( x% , u% ) )
≈ B ( ℑ ( x% + Bx ′ , u% + Bu
ˆ ′ + Nm
ˆ
% %
x ′t +1 = B* ℑ x% + Bx ′t + Nnt , u% + Bu
t
t − ℑ ( x, u )
*
t
t
)