Koren B., Vuik K. (Editors) Advanced Computational Methods in Science and Engineering

Подождите немного. Документ загружается.

Data Assimilation Algorithms for Numerical Models

The Gaussian measurement noise V

with zero mean and covariance matrix R(k)

represents uncertainties in the measurements. The initial condition X

is assumed to

be Gaussian with mean x

and covariance matrix P

. W

, V

and X

are assumed to

be mutually independent.

Having deﬁned the general stochastic state space representation of the model

(10) and the measurement relation (11) it is desired to combine the measurements

with the information provided by the model to obtain an optimal estimate of the

system state. To solve this ﬁltering problem the probability density of the state X

conditioned on the history of the measurements Z

, Z

, .. ., Z

taken, has to be

determined. Under the assumption described above it can be shown that this con-

ditional density function is Gaussian. Therefore it is completely characterized by

the mean X(k|k) and covariance matrix P(k|k). Moreover, the mean is in this case

equivalent to the least squares estimate or any other meaningful estimate. Recursive

equations to obtain these quantities can be summarized as follows:

Time update:

X(k | k −1) = F(k −1)X(k −1 | k −1) + B(k −1)u

k−1

, (12)

P(k |k−1) = F(k−1)P(k −1 |k−1)F(k −1)

+G(k −1)Q(k −1)G(k−1)

. (13)

Measurement update:

K(k) = P(k | k −1)M(k)



M(k )P(k | k −1)M(k)

+ R(k)



−1

, (14)

X(k | k) = X(k | k −1) + K(k)



−M(k)X(k | k −1)



, (15)

P(k | k) =



I −K(k)M(k)



P(k | k −1), (16)

where K(k) is called the Kalman gain.

The initial conditions for the recursive equations are

X(0 | 0) = x

,P(0 | 0) = P

The Kalman ﬁlter has a predictor-corrector structure. Based on all previous infor-

mation, a prediction of the state vector at time t is made by means of the equations

(12) and (13). Once this prediction is known it is possible to predict the next mea-

surement Z

by:

M(k )X (k | k −1).

When the measurement Z

has become available the difference between this mea-

surement and its predicted value is used to update the estimate of the state vector by

means of the equations (14)-(16). Note that the Kalman ﬁlter gain matrix K(k) does

not depend on the measurements and may be precomputed.

The Kalman ﬁlter does not only produce the optimal estimate of the system state,

but also computes the covariance matrix of the estimation error (13,16). This gives

113

A.W. Heemink et al.

insight into the accuracy of the estimate. If the covariance of the state estimate is

computed with and without the use of measurements, the difference of the two re-

sults indicates the improvement of the state estimate due to the measurement infor-

mation.

The variational approach described in Section 2 is based on minimizing a least

squares criterion. Since the Kalman ﬁlter also produces least squares estimates, the

ﬁlter results will be identical to the results of the variational approach for the same

linear problem formulation. This is, e.g., the case if the initial condition is chosen as

uncertain parameter in the variational approach and in the Kalman ﬁlter the system

noise is set to zero. The Kalman ﬁlter algorithm does not, however, explicitly com-

putes the estimate of the initial condition, but only the optimal estimate of the ﬁnal

state. If necessary the ﬁlter algorithm can be modiﬁed to also produce estimates of

previous states (Kalman smoother). This would produce the same estimate of the

initial condition as obtained with the variational approach. A discussion on optimal

smoothing is, however, out of the scope of this chapter.

For linear systems the Kalman ﬁlter is optimal for Gaussian noise processes.

Linearity of the system results in a Gaussian conditional probability density function

for X

given the measurements available up to and including time k. As a result the

mean X(k|k) and covariance matrix P(k|k) completely determine this probability

density function. Higher order moments do not have to be taken into account to

derive the optimal ﬁlter. This is not the case for non-linear systems. Due to the

non-linearities the density function for X

is non-Gaussian and has to be described

in general by an inﬁnite number of moments. This makes optimal ﬁlters in the non-

linear case computationally infeasible. However, it is possible to derive approximate

ﬁlters for many non-linear systems using linearization techniques.

3.2 Kalman ﬁltering for large-scale systems

3.2.1 Square root ﬁltering

Using the standard Kalman ﬁlter algorithm for systems with a very large dimension

n of the state vector would impose an unacceptable great computational burden. The

most time consuming part of the ﬁlter algorithm for large-scale systems is the ﬁrst

term in the time update of the covariance matrix(13):

F(k −1)P(k −1 | k −1)F(k −1)

For large-scale systems it is not practical to store the matrix F(k −1) and to de-

termine its transpose. F(k −1) merely represents a sequence of linear operations.

Therefore the time update equation of the covariance matrix is rewritten as:

F(k −1)



F(k −1)P(k −1 | k −1)



From this equation we can see that for every time step k the operator F(k −1) has to

114

Data Assimilation Algorithms for Numerical Models

be applied 2n times: ﬁrst for updating the n columns of the covariance matrix, then

for updating the n rows of the resulting matrix. As a result the computational effort

required for the algorithm is approximately equivalent with 2n model simulations.

In order to obtain a computationally feasible ﬁlter, simpliﬁcations have to be

introduced. However, there are serious problems with the more obvious simpliﬁca-

tions that one may consider. E.g., modifying the ﬁlter algorithm by forcing correla-

tions at large distances to be zero usually results in a covariance matrix that is not

semi positive deﬁnite. The last condition is physically not possible and often results

in a total failure of the ﬁlter algorithm. To avoid this problem a square root algorithm

can be used. This algorithm is based on the fact that a covariance matrix P can be

factorized as:

P(k −1 | k −1) = L(k −1 | k −1)L(k −1 | k −1)

where L is a square root of P. The factorization is not unique and there exist many

square roots of P. By rewriting the ﬁlter algorithm in terms of L instead of P a

more robust algorithm is obtained because approximations of L will never result in

a covariance matrix that is not semi positive deﬁnite. Most ﬁlter algorithms that have

been developed for large-scale systems are based on a square root representation of

the covariance matrix P where the square root of P is approximated by a matrix with

a reduced number of columns r. Hence, we have:

F(k −1)P(k −1 | k −1)F(k −1)

F(k −1)L(k −1 | k −1)



F(k −1)L(k −1 | k −1)



As a result only the r columns of L(k −1 | k −1) have to be updated yielding a

computational effort that is approximately equivalent with only r model simulations.

3.2.2 Classical Ensemble Kalman ﬁlter

The Ensemble Kalman ﬁlter (EnKF) introduced by Evensen (1994) is based on a

representation of the probability density of the state estimate by a ﬁnite number N

of randomly generated system states

(k −1 | k −1), i = 1,...,N.

The optimal estimate and the square root of the covariance matrix of the estimation

error are now given by:

X(k −1 | k −1) =

∑

i=1

(k −1 | k −1), (17)

L(k −1 | k −1) =

[ξ

(k −1 | k −1) −X(k −1 | k −1) ··· ξ

(k −1 | k −1) −X(k −1 | k −1)].

(18)

115

A.W. Heemink et al.

The square root L(k −1 | k −1) deﬁnes an approximation of the covariance matrix

P(k −1 | k −1) with rank N:

(k −1 | k −1) =

N −1

L(k −1 | k −1)L(k −1 | k −1)

. (19)

(k −1 | k −1) is however never actually computed. Using the algorithm, ﬁrst the

initial ensemble of state vectors is generated with mean x

and covariance matrix P

Then for updating the ensemble, realizations of the system noise and measurement

noise processes are generated too. The Ensemble Kalman ﬁlter algorithm can be

summarized as follows:

Time update:

(k | k −1) = f



(k −1 | k −1), k −1



+ G(k)w

, (20)

X(k | k −1) =

∑

i=1

(k | k −1), (21)

L(k | k −1) = [ξ

(k | k −1) −X (k | k −1) ··· ξ

(k | k −1) −X (k | k −1)]

(22)

Measurement update:

K(k) =

N−1

L(k | k −1)L(k | k −1)

M(k )



N−1

M(k )L(k | k −1)L(k | k −1)

M(k )

+ R(k)



−1

(23)

(k | k) = ξ

(k | k −1) + K(k)



−M(k)ξ

(k | k −1) + v



. (24)

The time update equation is computationally dominant. As a result the computa-

tional effort required for the Ensemble Kalman ﬁlter is approximately N model

simulations. The standard deviation of the statistical errors in the state estimate

converges very slowly with the sample size (≈ 1/N). This is one of the very few

drawbacks of this Monte Carlo type approach. Note that for the time update only

simulations with the original non-linear model are used. The tangent linear model

is not required.

3.2.3 Reduced Rank Square Root Kalman ﬁlter

The reduced rank square root ﬁlter algorithm (Verlaan and Heemink 1997) is based

on a factorization of the covariance matrix of the state estimate according to:

(k −1 | k −1) = L(k −1 | k −1)L(k −1 | k −1)

where L(k −1 | k −1) is a matrix with the q leading eigenvectors l

(k −1 | k −1)

116

Data Assimilation Algorithms for Numerical Models

of P(k −1 | k −1) as columns. These eigenvectors are scaled with the square root

of the corresponding eigenvalue. The computationally dominating part of the ﬁlter

algorithm is to determine:

F(k −1)L(k −1 | k −1).

Here the tangent linear model has to be applied to every column l

(k −1 | k −1) of

L(k −1 | k −1):

F(k −1)l

(k −1 | k −1), i = 1, ···, q.

To avoid the use of the tangent linear model we use the approximation:

F(k −1)l

(k −1 | k −1) ≈



f (

X(k −1 | k −1) + εl

(k −1 | k −1), (k −1) − f (

X(k −1 | k −1), k −1)



where ε represents a small perturbation. This approximation results in a time up-

date where only simulations with the original nonlinear model are required. The

algorithm can now be summarized as follows:

Time update:

X(K | k −1) = f



X(k −1 | k −1), k −1



+ B(k −1)u

k−1

, (25)

(k |k−1)=



f (

X(k−1|k−1)+εl

(k−1|k−1),(k−1)−f (

X(k−1 |k−1), k−1)



(26)

L(k |k−1) =

(k |k−1)

(k |k−1) ···

(k |k−1) G(k −1)Q

(k−1)

(27)

By computing:

L(k | k −1)

F(k −1)L(k −1 | k −1)



F(k −1)L(k −1 | k −1)



G(k)Q(k)G(k)

F(k −1)L(k −1 | k −1)L(k −1 | k −1)

F(k −1)

G(k)Q(k)G(k)

(28)

it is easy to show that equation (28) is equivalent to the time update (16) of the

original Kalman ﬁlter algorithm and that

L(k | k −1) is a square root approximation

of the covariance matrix P(k | k −1). However, a serious problem with this square

root is that its dimension has increased with respect to the dimension of the square

root L(k −1 | k −1) of P(k −1 | k −1). This would increase the computational

effort of the algorithm in every time step. This can be avoided by computing a new

approximate square root of P(k | k −1) with reduced dimension:

L(k | k −1) = Π (k | k −1)

L(k | k −1), (29)

117

A.W. Heemink et al.

where Π(k | k −1) is a projection onto the q leading eigenvectors of the matrix

L(k | k −1)

Using this reduction step we obtain again an approximate square root L(k | k −1) of

P(k | k −1) with only q columns.

Measurement update (for scalar measurements only):

β(k) =



N −1

M(k )L(k | k −1)L(k | k −1)

M(k )

+ R(k)



−1

, (30)

K(k) = L(k | k −1)L(k | k −1)

M(k )

β(k), (31)

X(k | k) =

X(k −1 | k −1) + K(k)



−

X(k −1 | k −1)



, (32)

L(k | k) = L(k | k −1) −

k(k)M(k)L(k | k −1)

1 + β(k)R(k)

. (33)

It is easy to show that for scalar measurements this measurement update is equiv-

alent to the measurement update (14)-(16) of the original Kalman ﬁlter algorithm.

Independent measurements (R(k) is diagonal) can be processed one by one using

the update equations (30)-(33). In case measurements are not independent of each

other they have to be transformed.

3.2.4 Discussion

One of the main practical advantages of the Kalman ﬁlter based algorithms is that

they are completely separate from the model implementation. The underlying model

f (X

,k) is used a number of times for different system states in order to compute

approximations of the mean and covariance of the state estimate. Neither a tangent

linear model nor an adjoint model implementation is required. As a result the coding

of the algorithm is very simple. Only an interface between the model code and the

ﬁlter code has to be developed.

For both the Ensemble Kalman ﬁlter as the Reduced Rank approach the choice

of the ensemble size is very essential for the performance of the algorithms. The

computational effort is approximately proportional with the ensemble size, so the

ensemble size has to be chosen as small as possible. However, if the ensemble size

is chosen too small the algorithm can become very inaccurate or unstable. Usually a

simple ”trial and error” procedure is the only way to get insight into an appropriate

value for the ensemble size.

Recently, many new variants of the ensemble type ﬁlter algorithms have appeared

in literature. Heemink et al. (2001) have been examining different approaches which

combine ideas from RRSQRT ﬁltering and the EnKF to derive computationally

118

Data Assimilation Algorithms for Numerical Models

more efﬁcient methods. The new algorithm, COFFEE (Complementary Orthogo-

nal subspace Filter For Efﬁcient Ensembles) is a RRSQRT algorithm where the

truncated modes are not neglected, but taken into account by a number of random

ensemble members (see also 6.8). Bishop et al. (2001) used an implementation of

the EnKF in an observation system simulation experiment. Ensemble-predicted er-

ror statistics were used to determine the optimal conﬁguration of future targeted

observations. The methodology was named Ensemble Transform Kalman Filter.

The EnKF can also be related to some other sequential ﬁlters such as the Singu-

lar Evolutive Extended Kalman (SEEK) ﬁlter by Pham et al. (1998) and the Error

Subspace Statistical Estimation (ESSE) by Lermusiaux and Robinsin (1999a,b) and

Lermusiaux (2001), which can be interpreted as an EnKF where the analysis is com-

puted in the space spanned by the EOFs of the ensemble. Anderson (2001) proposed

a method denoted the Ensemble Adjustment Kalman Filter, where the analysis is

computed without adding perturbations to the observations. Whitaker and Hamill

(2002) proposed another version of the EnKF where the perturbation of observa-

tions are avoided. The scheme provides a better estimate of the analysis variance

by avoiding the sampling errors of the observation perturbations. This is essentially

a Monte Carlo implementation of the square root ﬁlter and was named EnSRF. A

summary of the square root ﬁlters by Bishop et al. (2001), Anderson (2001), and

Whitaker and Hamill (2002) has been given by Tippet et al. (2003).

4 A software environment for data assimilation: COSTA

Data assimilation and calibration techniques are widely used in various modeling ar-

eas such us meteorology, oceanography and chemistry. Unfortunately it is very hard

to reuse the existing software implementing these techniques because they are in

general very model speciﬁc. Because existing software cannot be reused it is neces-

sary to program them from scratch in order to extend models with data assimilation

and model calibration techniques. COSTA

(Velzen, 2006 and Velzen and Verlaan,

2007 ) is a modular framework for data assimilation and model calibration. Within

the COSTA framework it is possible to combine models with the available data as-

similation and model calibration methods without the need of additional software

development. The usage of COSTA signiﬁcantly reduces the development costs for

implementing data assimilation for simulation models, especially when a model has

to be combined with various methods.

When new data assimilation methods are developed in COSTA they can be tested

on various models. Therefore their performance can be better determined and they

can be compared to other methods. Finally, data assimilation methods can be reused

and do not need to be programmed. COSTA simpliﬁes the application of data as-

similation methods such that it becomes available to a wider group of users and

application areas. Data assimilation and model calibration software that is specially

http://costapse.org

119

A.W. Heemink et al.

developed for a single model on the other hand are in general computationally ef-

ﬁcient. COSTA is set up in order to be as computationally efﬁcient as possible,

without losing its generic properties. The aim is that applications developed with

COSTA have a computational performance comparable to tailor made implemen-

tations. COSTA is set up such that it can be used for a range of different models

varying from small to large scale models and in sequential as well as parallel envi-

ronments.

4.1 COSTA components

COSTA deﬁnes a number of building blocks called components. A data assimilation

application is created by combining these components. The components can be split

into two groups; the high level components; the model, the observations and the

data assimilation method, and the low level components including among others the

vector, matrix, time and tree-vector.

A data assimilation application is created by combining the high level compo-

nents as illustrated in Figure 1. In most cases one or more of these components are

created using existing software. For that purpose some interfacing between the ex-

isting code and COSTA needs to be provided. This is illustrated by the red parts of

the components in Figure 1. The low level components provide a generic way for

data transfer between the high level components.

The principle of COSTA is that the programmer can choose whether to use a

generic COSTA implementation of a component or program a specialized imple-

mentation of a component. The COSTA development strategy is to quickly try out

alternative methods or models using the generic implementations of the components

and in a later stage improve performance or functionality by replacing some of the

components by specialized implementations, when necessary.

The components in COSTA are similar to classes in object-oriented languages.

This means that there can be multiple instances of a component and the data inside a

component is not directly accessible from outside. The data of a component can only

be manipulated by a set of functions, together called the interface of the component.

Different components with the same interface can therefore be interchanged and

reused for other applications. An additional advantage of shielding the data inside

a component is that the data does not need to be directly accessible. The data can

therefore e.g. be stored on ﬁle or in memory on another computer. It is therefore

possible to use the same interface for a black box model, where all interaction takes

place by reading box model output and writing COSTA input ﬁles as models that

use parallel computing.

120

Data Assimilation Algorithms for Numerical Models

Fig. 1 Creating COSTA components and combine them into an application. In this example exist-

ing model code and observation handling code is transformed into COSTA components. For this

purpose an interface layer (red) needs to be created. Thereafter it becomes possible to combine the

model with various data assimilation and calibration methods without any additional adjustments.

4.2 The COSTA software

COSTA is an open source software package developed at Delft University of Tech-

nology. The software can be freely downloaded at sourceforge.net, project costapse

and it is released under the LGPL license. Under this license it is possible to freely

use COSTA in combination with open source models and also in combination with

proprietary models. The COSTA project exists of two parts; the deﬁnition of com-

ponents and their interfaces and a software library, containing various data assimila-

tion methods, calibration methods, default implementation of various COSTA com-

ponents and tools for creating new COSTA components. The design of COSTA is

object oriented. This does not imply that COSTA can only be used in an object ori-

ented language. COSTA is developed in C and Fortran to ensure high performance

and currently it can be used from C, Fortran and Java.

The COSTA software environment is built up in different layers. The bottom

layer is not speciﬁc for data assimilation but provides the tools for programming

the COSTA components and contains software written by third parties. This layer

contains basic components like vectors, matrices and time but also robust and ef-

ﬁcient numerical methods, input/output facilities and data storage. On top of this

layer, the data assimilation and calibration speciﬁc components are constructed like

121

A.W. Heemink et al.

tree-vector, model, observations, and component builders. The top layer contains

complete (utility) programs.This is illustrated in Figure 2.

BLAS Functions LAPACK XML HDF SQLITE

Vector/MatrixStochastic utilities

Relation Tables Interpolations

Observation descriptions

Treevector

ModelStochastic observations

Model combiner Black−box model builder

Parallel and Domain Decomposition model builder

Data assimilation method

DA workbench

Third party implementations

Not data assimilation

specific

Data assimilation

and tools

specific components

Data assimilation

and utility programs

Level of

abstraction

Visualization

Fig. 2 An overview of the layers in COSTA. The bottom layer contains non data assimilation and

calibration speciﬁc software. The middle layer contains data assimilation components and the top

layer contains utility programs

A generic implementation is available for most COSTA components. The pro-

grammer can always choose to create his own implementation of a component. This

can be necessary for increasing the performance or the addition of new models and

methods. The interface of a component can consist of a large number of methods.

Fortunately this does not mean that all these methods need to be implemented when

programming a new version of a component. Some of the methods are redundant

and can be derived from others. The interface of the vector for example contains a

method for setting all values at once and a method for setting a single value in the

vector. Setting the whole vector can, although less efﬁciently, be implemented using

the method for setting a single value. COSTA will, when possible, automatically

implement missing methods in the case an implementation of a component does not

provide all the methods, The interfaces of the components contain this redundancy

because they simplify the usage of the component saving programming time and

yielding better software. The extra methods also provide ﬂexibility for optimizing

the implementation of the component.

We will illustrate why redundant methods can be very useful using the following

example. At the analysis step of a Kalman ﬁlter, the state is adjusted towards the

observations. A possible implementation, in pseudo code is:

122