## Xa Xb aXX Vyy Xb136

The variance of this estimate is

In spite of its simplicity, this example shows some of the key features of advanced linear data assimilation methods. First, note that the optimal estimate in Eq. (13.6)

is a linear combination of the background xb and the residual y - xb. The term

is called the analysis increment. Note the limits ox 0 (perfect background) and °y ^ 0 (perfect data), which yield xa=xb and xa =y, respectively. Furthermore, the estimated variance of the optimum (13.7) is less than the variance of either the background or the data, separately; combining information from the background and the observations has reduced uncertainty.

13.3.2.1 Example 2: Estimation of a Vector (Optimal Interpolation)

It is customary and instructive to generalize the above example to the estimation of a vector, assuming all errors Gaussian. This is Gauss-Markov smoothing, which forms the basis for many estimation algorithms. When the unknown vector represents values on a regular spatial grid, this procedure is known as optimal interpolation (Bretherton et al. 1976). The notation used here follows Ide et al. (1997).

Assume one wishes to estimate a vector x e RN, given a background xb, and a vector of observations y e RM. Optimal interpolation (also called objective analysis) is typically performed in order to smoothly interpolate sparse observations onto a regularly-spaced spatial grid, in which case x might represent, say, the value of sea-surface height at the grid points. Assume that each element of the observation vector y = {yi }f=x may be represented as the action of a linear operator on x(ti),

where h; e RlxN. For example, h. might extract the value of x at spatial coordinate (\$p X) in geographic latitude-longitude coordinates. Next, define the matrix h e RM xN by collecting the measurement operators together so that y=Hx. More generally, the observation operator H maps the model variables x to an equivalent of the observation vector, in terms of the locations and characteristics (e.g., units) measured.

To apply Bayes' Theorem, it is necessary to state the probability densities of x - xb and the measurement error s = y - Hx. These shall both be assumed to be multivariate Gaussian with zero means; the covariance of B e RNxN, i.e.,

and the covariance of e is denoted r e rmxM. Applying Eq. (13.2), one finds the pdf of x conditioned on y is proportional to exp ( - 1J(x)), where

J(x) = (x - xb)TB-1(x - xb) + (y - Hx)TR-1(y - Hx), (13.11)

which is the objective function that forms the basis for so-called variational data assimilation. The objective function, also called the cost function or penalty function, is to be minimized with respect to vector x.

With a bit of linear algebra, one finds the value x=xa which minimizes J is xa = xb + K(y - Hxb), (13.12)

where the analysis increment is Sxa=K(y - Hxb), and K takes the form

The full expression for the analysis error, the error covariance of xa, is denoted pa e rNxN

Remarks:

1. Equation (13.12) can be derived as the Best Linear Unbiased Estimator (BLUE), which minimizes expected error (eTK(xa - x))2, where ek e Rw is the basis vector pointing in direction k. Similarly, the estimator also minimizes the expected mean square error Tr((xa - x)(xa - x)T)/N. These facts are the basis for the equivalence of the BLUE, variational-, and Kalman Filter-based state estimates when the model dyanamics and measurement operators are linear.

2. When Optimal Interpolation is used in practice, the above formulas are often simplified so that the analysis at each grid point is computed from just the data within some nearby radius of influence Lorenc (1986); Daley (1991).

3. J(x) written above is also called the penalty function or cost function. The analysis field xa is its minimizer. Also, 1J is the negative log-likelihood function when the errors are Gaussian.

4. Because H is linear, J is convex and possesses a unique minimum. When H represents a nonlinear operator, there may be multiple minima.

5. Additional constraints may be added to the objective function, say, to suppress certain dynamics. These can complicate the solution procedure considerably and may obscure the failure of B or R to properly account for the covariance structure of the background and measurement errors.

6. The condition number of the objective function refers to the ratio of the largest to smallest eigen-values of the Hessian matrix of second derivatives of J. The eigenvalue spectrum of the Hessian can be interpreted as the curvature of principle axes of the iso-surfaces of J.

7. When the assumptions regarding the Gaussian errors are correct, the Hessian matrix H = d2 J/dx2 and the analysis error covariance Pa are related by Pa = 1H-1.

8. The above formalism can be applied to continuous fields, rather than vectors in RN. In that case x is generally a vector function and J is then a penalty functional. The stationarity condition for the minimum V J = 0 must be derived using the calculus of variations; the result is the Euler-Lagrange equation. There are close ties to the theory of smoothing splines, with B-1 being a positive-definite-symmetric differential operator (Wahba 1990). ### 2044 60 80 100 20 40 60 80 100

Fig. 13.4 Example 2: Optimal Interpolation. The left panel shows the sea-surface height ^(x, y) to be estimated from observations of n along idealized satellite ground tracks (open black circles filled with color indicating measured value), and from three observations of the surface current (filled black dots with arrows). The right panel shows the interpolated field (black contours) overlaid on the true sea-surface height. Parameters in the upper-right-hand corner are defined in the text. One can see that the large-scale flow features, such as the large anti-cyclonic eddy, have been reconstructed in the interpolated fields

### 2044 60 80 100 20 40 60 80 100

Fig. 13.4 Example 2: Optimal Interpolation. The left panel shows the sea-surface height ^(x, y) to be estimated from observations of n along idealized satellite ground tracks (open black circles filled with color indicating measured value), and from three observations of the surface current (filled black dots with arrows). The right panel shows the interpolated field (black contours) overlaid on the true sea-surface height. Parameters in the upper-right-hand corner are defined in the text. One can see that the large-scale flow features, such as the large anti-cyclonic eddy, have been reconstructed in the interpolated fields

9. The interpretation in terms of continuous fields is essential to properly understanding the conditioning of the objective function as model resolution is increased to the continuum limit (Bennett and Budgell 1987; Bennett 1992). The language of linear algebra is not adequate for analyzing spatial regularity (differentiability) of the analysis increments.

An example of optimal interpolation is shown in Fig. 13.4. The vector x to be estimated is the sea-surface height on a regular 1 km grid within a box of dimensions 100 km by 100 km. Denote the sea-surface height field by ^(x, y), with elements of x being ^(xj, yj) on the grid. In this idealized setup there are 30 observations of sea-surface height

yi = e x for i = 1,..., 30 along three satellite ground tracks (color-filled dots). It is assumed that the currents are in geostrophic balance with the surface pressure gradient, and there are 3 observations of near-surface current (black-dots with arrows). The u-component is given by 