Using Analysis State to Construct a Forecast Error Covariance Matrix in Ensemble Kalman Filter Assimilation

Data assimilation is a method to produce an optimal description of the system state based on model prediction and observations (Talagrand, 1997). Accurately estimating the forecast error covariance matrix is one of the most important steps in data assimilation. If it is correctly known, the analysis states can be generated by minimizing the objective function, which is a fairly technical aspect and can be accomplished with existing engineering solutions (Reichle, 2008).

Ensemble Kalman filter (EnKF) is a popular data assimilation scheme, which has been widely used in atmospheric and oceanic fields since it was first introduced by Evensen (1994a, b). In EnKF, the forecast error covariance matrix is estimated as the sampling covariance matrix of the ensemble forecast states. However, past works on EnKF have found that the sampling errors in such estimation, resulting from the limited ensemble size, as well as the poor initial perturbations and model error, can generally lead to an underestimate of the forecast error covariance matrix and eventually result in filter divergence (e.g. Anderson and Anderson, 1999; Constantinescu et al., 2007).

To compensate for this defect, using forecast error covariance matrix inflation to increase ensemble variance is increasingly important. One of the inflation techniques is additive inflation, in which a noise is added to the ensemble forecast states that samples the probability distribution of model error (Hamill and Whitaker, 2005). Another widely used forecast error covariance matrix inflation technique is multiplicative inflation; that is, to multiply the matrix by an appropriate factor.

At the early state, the inflation factor is estimated empirically. (Wang and Bishop, 2003) developed an approach to optimize inflation factor. Their scheme was further improved by (Li et al., 2009). All these schemes are based on first moment estimation of squared observation-minus-forecast residual, which was first introduced to data assimilation by (Dee, 1995). (Zheng, 2009) and (Liang et al., 2012) proposed another approach to optimize the inflation factor, but using maximum likelihood estimation (MLE), which was introduced by Dee and colleagues (Dee and da Silva, 1999; Dee et al., 1999). Anderson (2007, 2009) also proposed an inflation method based on a Bayesian approach, and (Miyoshi, 2011) further simplified Anderson's inflation approach by making a number of additional simplifying assumptions. However, their study was limited to the case that observational error is spatially independent.

In this paper, analysis states are used to construct the forecast error covariance matrix (which is different from the sampling covariance matrix of ensemble forecast states used in conventional EnKF). The methodology documented by (Liang et al., 2012) is further developed by introducing an adaptive procedure for estimating the proposed forecast error covariance matrix. By definition, the ensemble forecast error should be represented by ensemble forecast states minus true states rather than by ensemble forecast states minus the ensemble mean of the forecast states (Evensen, 2003). For the model with large error, the ensemble mean of forecast states can be far from the true states. Therefore, the sampling covariance matrix of the ensemble forecast states can be far from the true forecast error covariance matrix. As a result, the estimated analysis state cannot be accurate enough. In reality, it is impossible to know the true state exactly, but the analysis state is a better estimation of the true state than the forecast state. Therefore, in this paper, we propose to use the information feedback from analysis state to update the forecast error covariance matrix. In fact, our proposed forecast error covariance matrix is a combination of multiplicative and additive inflation. Bai and Li (2011) also used the feedback from the analysis state to improve assimilation, but in a different way.

The remainder of the paper is organized as follows. Section 2 summarizes the EnKF scheme with the MLE inflation scheme of the forecast error covariance matrix (Zheng, 2009; Liang et al., 2012), as well as our proposed structure of the forecast error covariance matrix and its adaptive estimation procedure. Sections 3 and 4 respectively present the assimilation results from the low-order Lorenz model (Lorenz, 1996) and the more realistic 2D Shallow Water Equation (SWE) model with a larger correlated observational system. In particular, our approach is compared with two published forecast error covariance matrix inflation techniques: the MLE inflation scheme (Liang et al., 2012), and the first moment inflation scheme (Wang and Bishop, 2003). Finally, a discussion and conclusions are given in section 5.

Using the notations of (Ide et al., 1997), a nonlinear discrete-time forecast and linear observational system is written as:

x_t,i= M_i-1(x_a,i-1)+η_i, (1)

Y_o,i=H_iX_t,i+ε_i, (2)

where i is the time step index; x_t,i ={ x_t,i (1), x_t,i (2),…, x_t,i (n)}^T is the n-dimensional true state vector at the time step i; x_a,,i-1 ={ x_a,i-1 (1), x_a,i-1 (2),…, x_a,i-1 (n)}^T is the n-dimensional analysis state vector, which is an estimate of x_t,i-1; M_i is the nonlinear forecast operator, such as a weather forecast model; y_o,I is the observational vector with dimension P_i; H_i (linear in this paper) is the measurement matrix of dimension P_i×n that maps the model states to the observed variables; η_i andε_i are the forecast error and observational error vectors, which are assumed to be statistically independent of each other, time-uncorrelated, and have zero mean and covariance matrices i and R_i respectively. The goal of EnKF assimilation is to find a series of analysis states x_a,i that is sufficiently close to the true state x_t,i, using information provided by M_i and y_a,i.

It is well known that any EnKF assimilation scheme should include a forecast error inflation scheme; otherwise, the EnKF may diverge (Anderson and Anderson, 1999). In this paper, a procedure for estimating the multiplicative inflation factor of P_i is proposed based on the MLE principle (e.g. Dee and da Silva, 1999; Zheng, 2009; Liang et al., 2012). The basic filter algorithm in this paper uses perturbed observations (Burgers et al., 1998) without localization (Houtekamer and Mitchell, 2001). The estimation steps of this algorithm equipped with MLE inflation (Zheng, 2009; Liang et al., 2012) are as follows.

Step 1: Calculate the perturbed forecast states

X_f,I,j = M_i-1(x_a,i-1,j) (3)

where x_a,i-1,j is the perturbed analysis states derived at the previous time step (1jm and m is the number of ensemble members). Define the forecast state x_f,I to be the ensemble mean of x_f,I,j.

and then inflate it to the form of . The inflation factorλ_i is estimated by minimizing the objective function

where is the observation-minus-forecast residual and L_j(λ) is the 2log-likelihood of d_j. The relatively efficient calculations for the determinant and inverse in Eq. (5) are documented in Appendix A, as well as in (Liang et al., 2012).

Step 3: Compute the perturbed analysis states

where is the estimated inflation factor andε’_i,j is normally distributed with mean zero and covariance matrix R_i. Then, define the estimate of the analysis state x_a,i to be the ensemble mean of x_a,I,j. Finally, set i=i+1 and return to Step 1.

By Eqs. (1) and (3), the ensemble forecast error is x_f,i,j - x_t,i. x_f,i is an estimate of x_t,i without knowing observations. The ensemble forecast error is initially estimated as x_f,i,j - x_t,i, which is used to construct the forecast error covariance matrix in section 2.1. However, owing to limited ensemble size and model error, x_t,i can be biased. Therefore, x_f,i,j - x_t,i can be a biased estimate of x_f,i,j - x_t,i.

In this paper, we propose to use observations for improving the estimation of ensemble forecast error. The idea is as follows: After analysis state x_a,i is derived, it is generally a better estimate of x_t,i than the forecast state x_t,i. So x_t,i in Eq. (4) is substituted by x_a,i for updating the forecast error covariance matrix. This procedure is repeated iteratively until the corresponding objective function [Eq. (5)] converges. For the computational details, Step 2 in section 2.1 is modified to the following adaptive procedure:Step 2a: Use Step 2 in section 2.1 to inflate the initial forecast error covariance matrix to . Then, use Step 3 in section 2.1 to estimate the initial analysis state and set .

Step 2b: Update the forecast error covariance matrix as

Then, adjust the forecast error covariance matrices to , where is estimated by minimizing the following objective function

If , where δ is a predetermined threshold to control the convergence of Eq. (8), then estimate the kth updated analysis state as

set k=k+1 and return to Eq. (7); otherwise, accept as the estimated forecast error covariance matrix at the ith time step, and go to Step 3 in section 2.1.

A flowchart of our proposed assimilation scheme is shown in Fig. 1. Moreover, our proposed forecast error covariance matrix [Eq. (7)] can be expressed as

which is a multiplicatively inflated sampling error covariance matrix plus an additive inflation matrix (see Appendix B for the proof).

Figure 1.  Flowchart of the proposed assimilation scheme.

For ideal models, the "true" state x_t,i can be calculated. So, we can use the RMSE of the analysis state to evaluate the accuracy of the assimilation results. The RMSE at the ith step is defined as

where ||.|| represents the Euclidean norm and n is the dimension of the state vector. In principle, a smaller RMSE means a better performance of the assimilation scheme.

In this section we apply our proposed data assimilation schemes to a nonlinear dynamical system with properties relevant to realistic forecast problems: the Lorenz-96 model (Lorenz, 1996) with model error and linear observational system. We evaluate the performances of the covariance inflation methods described in section 2 through the following experiments.

The Lorenz-96 model (Lorenz, 1996) is a strongly nonlinear dynamical system with quadratic nonlinearity, governed by the equation

where k=1,2,…,40. To make Eq. (12) meaningful for all values of k, we define r_-1= r_39, r₀= r₄₀, r₄₁= r₁. The dynamics of Eq. (12) are "atmosphere-like" in that the three terms on the right-hand side consist of a nonlinear advection-like term, a damping term, and an external forcing term respectively. These terms can be thought of as some atmospheric quantity (e.g. zonal wind speed) distributed on a latitude circle.

In our assimilation schemes, we set F=8. For this parameter setting, the leading Lyapunov exponent implies an error-doubling time of about eight time steps, and the fractal dimension of the attractor is 27.1 (Lorenz and Emanuel, 1998). The initial condition is chosen to be r_k=F when k≠20 and r₂₀ =1.001F. The initial ensemble perturbation states are generated by adding a white noise with standard deviation F/40 to the initial value. We solve Eq. (12) using a fourth-order Runge-Kutta time integration scheme (Butcher, 2003) with a time step of 0.05 nondimensional units to derive the true state. This is roughly equivalent to six hours in real time, assuming that the characteristic time scale of the dissipation in the atmosphere is five days (Lorenz, 1996).

In this study, we assume the synthetic observations are generated at every model grid point by adding random noises which are multivariate-normally distributed with zero mean and covariance matrix R_i to the true state. In our experiments, the observational errors are spatially correlated, which may potentially be the case for remote sensing observations and radiances data. The leading-diagonal elements of R_i areσ²₀ =1 and the off-diagonal elements at site pair (j,k) are

We would like to introduce model error in the Lorenz-96 model, because model error is inevitable in real dynamic systems. Thus, we chose different values of F in our assimilation schemes, while retaining F=8 when generating the "true" state. We assimilate observations every four time steps for 2000 steps and select the ensemble size as 30. The predetermined thresholdδto control the convergence of Eq. (8) is set to be one.

In this section we further test our proposed adaptive assimilation scheme using a larger model, the 2D SWE model (see Liang et al., 2012).

: The SWEs used in the present study are the equations of motion that can be used to describe the horizontal structure of an atmosphere. They describe the evolution of an incompressible fluid in response to gravitational and rotational accelerations. The barotropic nonlinear SWE consists of the following partial differential equations (Lei and Stauffer, 2009):

Where μ is the fluid velocity in the first direction (or zonal velocity); νis the fluid velocity in the second direction (or meridional velocity); h is the fluid height (or depth); g is the acceleration due to gravity, which is set to 9.8 m s^-2; f is the Coriolis coefficient associated with the Coriolis force, defined as the constant 10^-4 s^-1 using the f -plane approximation; and the diffusion coefficientκis given as 10⁴ m² s^-1. Further, L and D are the length and width of the rectangular domain of integration, which are set to 500 km and 300 km respectively. The rectangular domain is divided into 50x31 dimensions using a uniform grid spacing of 10 km in the two directions. The time step is 30 seconds and the Lax-Wendroff method is used to integrate the model with respect to time. The periodic boundary condition is used at the first direction boundaries, and a free-slip rigid wall boundary condition (where μ and h are defined from the values one point inside the boundary) is used at the second direction boundaries. See (Liang et al., 2012) for a more detailed explanation.

where H₀, H₁ and H₂ are set to 50.0 m, 5.5 m and 3.325 m respectively. The initial ensemble height values are created by adding a white noise with unit variance to the initial height. The initial velocity field and the ensemble members are derived from the initial height fields with the geostrophic relation.

The observational errors are specified asσ²_oh =0.5 m² s^-2, σ²_oh =0.5 m² s^-2, σ²_oh=1.0 m², and are spatially correlated with each other. The correlation coefficient between the grid points z_k and z_k’ is ρ^d(zk,zk’), where d(z_k,z_k’) is the distance between the two grid points and the parameter ρ is set to 0.5. Furthermore, the observational errors of different variables are assumed to be independent of each other. We use the diffusion coefficient k=5x10^-4m² s^-1 to generate model error in our assimilation process. The ensemble size is 100. The predetermined threshold δ is also set to be one.

Figure 4.  The same as Fig. 3, but for the SWE model.

We integrate SWE for 72 hours and the results are considered as the "true" states and have synthetic observations of μ, ν and h at randomly located 310 grid points for every 3 hours in a 48-hour assimilation period. A 24-hour forecast is taken after the 48-hour assimilation period. The time series of analysis RMSE of componentsμ, ν and h of the three assimilation schemes in the 72-hour period are plotted in Fig. 4. Since the first 24-hour period is treated as spin-up and the last 24-hour period is for validation of prediction, we only calculate the time-mean values between 25-hour and 48-hour. The results are shown in Table 2.

Figure 4 and Table 2 show that the time-mean RMSE of EnKF with our proposed structure of forecast error covariance matrix is smaller than the EnKF with the MLE inflation scheme and with the first moment inflation scheme, especially for the component ν.

In all data assimilation schemes, the analysis state is estimated as a weighted average of model forecasts and observation. The weights are related to the corresponding error covariance matrices. Therefore, the performance of data assimilation depends crucially on how accurately the forecast and observational error covariance matrices can be estimated. An assimilation scheme that can estimate them more accurately should have a better assimilation result.

In fact, the true forecast error should be represented as the ensemble forecast states minus the true state. Since in real problems the true state is not available, we use the ensemble mean of the forecast states instead. Consequently, the forecast error covariance matrix is represented as the sampling covariance matrix of the ensemble forecast states. However, if the model error is large, the ensemble mean of the forecast states may be far from the true state. In this case, the estimated forecast error covariance matrix will remain far from the truth, no matter which inflation technique is used.

For verifying the aforementioned point, EnKF with error covariance matrix inflation was applied to the Lorenz-96 model and SWE model, but the forecast state x_f,i in the forecast error covariance matrix [Eq. (4)] was substituted by the true state x_t,i. The corresponding RMSEs are also shown in Figs. 2-4 and Tables 1-2. All the figures and tables show that the RMSE was very significantly reduced compared with the MLE inflation scheme and the first moment inflation scheme.

However, since the true state x_t,i is not known, we used the analysis state x_a,i to substitute the forecast state x_f,i, because x_a,i is closer to x_t,i than x_f,i. To achieve this goal, a proposed structure of the forecast error covariance matrix and an adaptive estimation procedure for estimating the new structure have been introduced here to constantly improve the assimilation result. As shown in sections 3 and 4, the RMSEs of corresponding analysis states are indeed smaller than those of the other two schemes.

The -2log-likelihood function of observation-minus-forecast seems a good objective function to quantify the goodness-of-fit of the forecast error covariance matrix. In most cases in this study, the number of iterations was only 3 or 4 and the objective function decreased sharply. In fact, as shown in Tables 1-2, the smaller -2log-likelihood always corresponded to a smaller RMSE of the analysis state. On the other hand, the improved forecast error covariance matrix indeed leads to the improvement of analysis state. All the figures and tables also show that analysis RMSE of EnKF with the MLE inflation scheme is smaller than that with the first moment inflation scheme. This may not be generally true, but could be case dependent.

Similar to other EnKF assimilation with single parameter inflation schemes, this study also assumed the inflation factor is constant in space. Apparently, this is not the case in many practical applications, especially when the observations are unevenly distributed. Persistently applying the same inflation values that are reasonably large to address problems in densely observed areas to all state variables can systematically overinflate the ensemble variances in sparsely observed areas (Hamill,2005; Anderson,2009). Even if the adaptive procedure for estimation of the forecast error covariance matrix is applied, the problem may still exist to some extent. In the two case studies conducted in this paper, the observational systems were relatively evenly distributed.

From the assimilation results shown in Fig. 2, the improvement is minor when there is no model error or model error is small (e.g. F=8,7,9). The improvement is more significant when model error is larger (e.g. F=8,12). Since the model error in numeric weather prediction (NWP) models is generally not very large (otherwise, the 4DVar scheme will not work), we suggest that the proposed scheme may not be able to significantly improve NWP data assimilation. However, we expect the proposed scheme may be more useful to improve hydrological data assimilation and land surface data assimilation, since the model errors in hydrological models and land surface models are generally larger.

In future work, we will investigate how to modify the adaptive procedure to suit a system with unevenly distributed observations. We also plan to use our proposed scheme to improve the estimation of the forecast error covariance matrix in the data assimilation schemes developed by our colleagues (e.g. Liu et al., 2008; Tian et al., 2011; Bai and Li, 2011). In this way, we hope that data assimilation capabilities in China will be considerably improved.

Efficient calculations of and

Provided that R_i^-1 is inexpensive to obtain, : can be calculated using the Sherman-Morrison-Woodbury identity (Golub and van Loan, 1996), as is proposed in (Tippett et al., 2003):

where I is the p_i×p_j identity matrix, and the n×m matrix X_f,i is the square root of :

For ensemble filters, X_f,i can be obtained through the deviations of ensemble forecast states to their mean:

In this way, only m×m matrices need to be inverted and thus the computational cost is significantly reduced.

can be calculated efficiently using the singular value decomposition of the matrix :

where R_i^-1/2 is the square root of R_i^-1, U_i and V_i are orthogonal matrices; D_i is a p_i×m matrix whose leading-diagonal elements are the singular values of R_i^-1/2H_iX_f,i. Then, we have:

where d_i,j is the jth diagonal element of D_i.

In fact,

(B1)

Since x_f,i is the ensemble mean forecast, we have

(B2)

and similarly

(B3)

So the last two terms of Eq. (B1) vanish. Therefore

(B4)

In this sense, what we proposed is a scheme combining a multiplicative inflation and an additive inflation.

Suppose that the forecast errors and observational errors are uncorrelated, the following relationship is satisfied:

where <> represents the mathematical expectation operator. The vector d_i in this operator is regarded as a random vector.

As is stated in (Wang and Bishop, 2003), for the global observational networks the number of independent scalar elements within R_i^-1/2d_i is large, and thus distributes closely around its mean value . With this assumption, Eq. (C1) becomes: