HTML
--> --> -->However, due to the limited sample size and imprecise dynamics model evolution function, it may be easy for the forecast variance to be underestimated (Anderson and Anderson, 1999; Evensen, 2009). Therefore, over time, the observations may have a small impact on the estimation process, and further lead to the phenomenon of filter divergence.
Therefore, dealing with the problem of forecast error variance underestimation becomes a necessary way to prevent filter divergence in this case, and the variance inflation technique is an effective procedure to address this problem.
In early research, the variance inflation parameter was determined by repeated experiments or prior knowledge (Anderson and Anderson, 1999). (Wang and Bishop, 2003) proposed an online estimation of the inflation parameter by a sequence of innovation statistics. Based on their work, (Li et al., 2009) further developed an algorithm that can simultaneously estimate the inflation parameter and observation errors. (Wu et al., 2013) proposed a method to estimate the inflation parameter with a second-order least-squares method. Although these methods can effectively address the phenomenon of forecast variance underestimation, sometimes the distance between the actual observation and its estimate given the analysis state of the above methods may be too large. In this situation, the estimate of the state cannot be regarded as a reasonable output of the corresponding assimilation method (Anderson and Anderson, 1999). Besides, the true initial condition is usually unknown in most practical applications, and a data assimilation method is frequently initialized by guessing or prior experience, and thus the estimation process may take a long time to become stable if the initial state is settled far away from the true value. Hence, reducing this time is also very important in practical applications, particularly when the estimate is needed as soon as possible (Linderoth et al., 2011; Wang et al., 2014).
To relieve the phenomenon of filter degeneracy and reduce the convergence time, here we judge the output of the data assimilation method in real time based on a confidence region constructed with the observation information. This idea is quite similar to the quality control method. Quality control is an online procedure to monitor and control an ongoing production process (MacGregor and Kourti, 1995; Montgomery, 2007). If the outputs of the process fall outside a pre-specified limit, the process is regarded as out of control and an adjustment measure should be taken. Based on this idea, we propose a method using the idea of quality control to address the issues of EnKF stated above, and we call this method EnCR, in which "CR" stands for "Confidence Region". In the new method, the filter process is regarded as a production process and should be adjusted when the predicted observation given the analysis value falls outside of the pre-specified limits around the actual observation. This is the basic idea of the EnCR method. In this way, the predicted observation given the analysis state will not differ too much from the actual observation, and thus the estimate (output) will be more reasonable.
The rest of this paper is organized as follows: In section 2, we summarize the popular EnKF method and give a brief description of the phenomenon of forecast error variance underestimation. In section 3, we propose a new method to address the problem of inflation parameter estimation in the EnKF method. Various simulation results in the Lorenz-63 and Lorenz-96 models are presented in section 4. A summary of the new method is detailed in section 5. The technical proofs are given in the Appendix.
Consider a nonlinear discrete time state space model below: \begin{equation} \label{eq1} \left\{ \begin{array}{rcl} {x}_t&=&{M}_t({x}_{t-1})+{\eta}_t ,\\ {y}_t&=&{H}_t{x}_t+{\varepsilon}_t , \end{array} \right. \ \ (1)\end{equation} where t is the time index, \(x_t\in\Re^r\) is an r-dimensional state vector at time t, Mt(·) is a known nonlinear forecast operator at time t, and ηt is a white, Gaussian distribution with mean 0 and variance Qt and is often referred to as the model error. Moreover, \(y_t\in\Re^n\) is an n-dimensional observation vector at time t, Ht is a known linear observation operator at time t, and εt is a white, Gaussian distribution with mean 0 and variance Rt and is often referred to as the observation error. The error terms are assumed to be statistically independent and time-uncorrelated.
When the forecast error operator Mt(·) is linear and Qt is known, model (2) can be solved by the classical Kalman filter. Furthermore, when Mt(·) is nonlinear and Qt is known, the extend Kalman filter, EnKF or the particle filter (PF) (Gordon et al., 1993; Salmond and Gordon, 2005) can be used to solve this model. However, when Mt(·) is nonlinear and Qt is unknown, model (2) can be solved by the EnKF or PF, in which the forecast ensemble is obtained by the evolution of dynamic model without involving the error term. Meanwhile, the corresponding ensemble can be used to estimate the forecast variance. Before proceeding with the EnKF procedure, we denote the ensemble at time t as {xt,i,i=1,2,…,m}, where m represents the sample size. In this paper, we denote the forecast ensemble with subscript f and the analysis ensemble with subscript a. Thus, the EnKF procedure can be stated as follows:
First, the initial distribution is usually assumed to be known and generated from a Gaussian distribution: $$ {x}_{0,i,{\rm a}}\sim N(\mu_{0},Q_0),1\le i\le m , $$ where μ0 and Q0 are assumed to be known and usually specified by experience.
The forecast ensemble and its error variance can be estimated by \begin{eqnarray} \label{eq2} {x}_{t,i,{\rm f}}&=&{M}_t({x}_{t-1,i,{\rm a}}),{x}_{t,{\rm f}}=\frac{1}{m}\sum_{i=1}^m{x}_{t,i,{\rm f}} ,\ \ (2)\\ \label{eq3} \hat{P}_{t\vert t-1}&=&\frac{1}{m-1}\sum_{i=1}^m({x}_{t,i,{\rm f}}-{x}_{t,{\rm f}})({x}_{t,i,{\rm f}}-{x}_{t,{\rm f}})' .\ \ (3) \end{eqnarray}
Then, the analysis ensemble can be obtained by \begin{eqnarray} \label{eq4} {x}_{t,i,{\rm a}}&=&{x}_{t,i,{\rm f}}+\hat{P}_{t\vert t-1}{H}'_t({H}_t\hat{P}_{t\vert t-1}{H}'_t+R_t)^{-1}({y}_t+{\gamma}_{t,i}-{H}_t{x}_{t,i,{\rm f}}) ,\nonumber\\ &&1\le i\le m ,\ \ (4) \end{eqnarray} where γt,i~ N(0,Rt),i=1,2,… m. Here, yt+γt,i can be viewed as a perturbed observation of the observation at time t. Therefore, the analysis state can be estimated by \begin{equation} \label{eq5} {x}_{t,{\rm a}}=\frac{1}{m}\sum_{i=1}^m{x}_{t,i,{\rm a}} . \ \ (5)\end{equation}
From the procedure described above, we can ascertain that the estimated forecast error variance \(\hat{P}_t\vert t-1\) in Eq. (5) is underestimated because the true forecast error variance can be represented by \begin{eqnarray} \label{eq6} P_{t\vert t-1}&=&E(({x}_t-{x}_{t,i,{\rm f}})({x}_t-{x}_{t,i,{\rm f}})')\nonumber\\ &=&\mathop{\lim}\limits_{m\to\infty}(\hat{P}_{t\vert t-1}+E(({x}_t-{x}_{t,{\rm f}})({x}_t-{x}_{t,{\rm f}})'))\nonumber\\ &=&\mathop{\lim}\limits_{m\to\infty}(\hat{P}_{t\vert t-1})+Q_t . \ \ (6)\end{eqnarray}
Hence, over time, the forecast error variance may be significantly underestimated and the phenomenon of filter degeneracy may occur.
Moreover, even though the model error Qt is known, the problem of forecast error variance underestimation may also exist because of the existence of spurious correlations over long spatial distances or between variables known to be uncorrelated (Evensen, 2009). Hence, the filter may also become degenerated over time, as depicted in the simulation example of the perfect model case in subsection 4.2.2 or the explanations of Evensen (2009, Chapter 15) about this issue.
Therefore, the forecast error variance is typically multiplied by a parameter that is bigger than one to relieve the problem of forecast error variance underestimation; that is, \begin{equation} P_{t\vert t-1}=\lambda_t\hat{P}_{t\vert t-1} ,\ \ (7) \end{equation} where Λt is the adjustment (inflation) parameter. Then, the inflated forecast error variance \(\lambda_t\hat{P}_t\vert t-1\) is used to replace the original \(\hat{P}_t\vert t-1\) in Eq. (5). The adjustment parameter can be estimated by many methods, e.g., the W-B (Wang and Bishop, 2003) or the SLS [Second-order Least Squares, (Wu et al., 2013)] methods. However, these methods do not consider whether the forecast observation given the analysis state differs too much from the newest observation. Next, we present the estimation of Λt with the constraints of confidence region.
From Eq. (8), we can see that if the relative probability of yt given xt is too small——that is, the prediction observation given the analysis state diverges too far from the actual observation——the final conditional PDF of xt will also be very small. This means that xt is unlikely to be the true state (Anderson and Anderson, 1999)——that is, the analysis state xt, a cannot be regarded as a reasonable estimate if the value of \(p(y_t\vert x_t,\rm a)\) is too small. Therefore, we can evaluate the outcome using this perspective. From this point of view, we construct a confidence region based on the observation to calculate the inflation parameter in the EnKF.
2
3.1. EnCR method
Based on the descriptions above, we introduce a confidence region to testify the feasibility of the analysis state. This idea is similar to the quality control method. Quality control is a method that aims at monitoring the extent to which products meet specifications. By monitoring the quality characteristics, the abnormality of a production procedure can be detected quickly and a measure can be taken in time to guarantee that the process is under control (MacGregor and Kourti, 1995; Montgomery, 2007). Statistical principles are widely used in quality control and, here, we briefly introduce its procedure based on the χ2 statistic.Suppose the quality characteristic is x, \(\hat\Sigma\) is the corresponding estimated in-control covariance, and μ is a pre-specified value. The control region can be stated as $$ D_L=\{{x}:({x}-{\mu})'\Sigma^{-1}({x}-{\mu})<\chi_\alpha^2\} , $$ where χα is the α quantile of the χ2 distribution. If x falls in the region DL, the production process is considered to be normal; otherwise, the process is regarded as out of control and some adjustment procedure needs to be taken. Usually, α is set as 0.99; that is $$ P(({x}-{\mu})'\Sigma^{-1}({x}-{\mu})\le L)=0.99 . $$
Similarly, we can view the EnKF method with an unknown adjustment parameter as a production process, and the analysis state can be viewed as the corresponding product of the process. Hence, we can determine the parameter based on the idea of quality control and construct the corresponding control region as follows: First, the product (analysis state) xt, a can be expressed as \begin{equation} \label{eq8} {x}_{t,{\rm a}}(\lambda_t)={x}_{t,{\rm f}}+\lambda_t\hat{P}_{t\vert t-1}{H}'_t(\lambda_t{H}_t\hat{P}_{t\vert t-1}{H}'_t+R_t)^{-1} ({y}_t-{H}_t{x}_{t,{\rm f}}) , \ \ (9)\end{equation} which is controlled by parameter Λt, and xt, f represents the prior forecast of the state. The observation can be viewed as the pre-specified value, and the control statistic can be constructed by \begin{equation} \label{eq9} u_t(\lambda_t)=({y}_t-{H}_t{x}_{t,{\rm a}}(\lambda_t))'V_t^{-1}({y}_t-{H}_t{x}_{t,{\rm a}} (\lambda_t)) , \ \ (10)\end{equation} where Vt is the in-control covariance matrix that can be calculated by $$ V_t=E(({y}_t-{H}_t{x}_{t,{\rm a}}(\lambda_t))({y}_t-{H}_t{x}_{t,{\rm a}}(\lambda_t))') . $$ The parameter Λt needs to be adjusted when xt, a falls outside the control region (confidence region): \begin{equation} D=\{\lambda_t:u_t(\lambda_t)<L\} , \ \ (10)\end{equation} where L is the α quantile of the χ2(n) distribution and is usually set between 0.90 and 0.99, and n is the dimension of yt.
Theorem 3.1: Under assumption (7), we have \begin{equation} \label{eq10} V_t=R_t(\lambda_t{H}_t\hat{P}_{t|{t-1}}{H}'_t+R_t)^{-1}R_t \ \ (12)\end{equation}
Theorem 3.2: Under assumption (7), the analysis state variance of EnCR can be expressed as \begin{equation} \label{eq11} E(({x}_t-{x}_{t,i,{\rm a}})({x}_t-{x}_{t,i,{\rm a}})')=(\lambda_t^{-1}\hat{P}_{t|{t-1}}^{-1}+{H}'_t R_t^{-1}{H}_t)^{-1} . \ \ (13)\end{equation}
Theorem 3.3: The control statistic can be expressed as \begin{equation} \label{eq12} u_t(\lambda_t)=({y}_t-{H}_t{x}_{t,{\rm f}})'(\lambda_t{H}\hat{P}_{t\vert t-1}{H}'+R)^{-1} ({y}_t-{H}_t{x}_{t,{\rm f}}) , \ \ (14)\end{equation}
Based on the results of theorem 3.3, we can ascertain that Λt is a monotonically decreasing function of ut(Λt); therefore, the adjustment parameter can be constructed by Eqs. (11) and (13). That is, if the initial value \(\lambda_t\notin D\), we increase Λt until it falls in the region D. In this paper, the parameter Λt that first falls in region D is used as its estimate. We also note that the calculation of Λt in Eq. (13) only involves the dimension of the observation, which is time-saving, especially in some high-dimensional situations, if the dimension of the observation is much smaller than the state.
Next, the algorithm procedure of the EnCR is presented.
2
3.2. EnCR procedure
Based on the descriptions stated above, the EnCR procedure can be expressed as follows:(1) Set the initial values x0 and Q0; and define L to be the α quantile of χ2(n). Here, we set α=0.99 in the simulation examples.
(2) Generate the initial ensemble x0,1, a,x0,2, a,?,x0,m, a from N(u0,Q0) and set t=1.
(3) Set the initial parameter Λt=1 and calculate the forecast state and its variance by \begin{eqnarray*} {x}_{t,i,{\rm f}}&=&{M}_t({x}_{t-1,i,{\rm a}}) ,{x}_{t,{\rm f}}=\frac{1}{m}\sum_{i=1}^m{x}_{t,i,{\rm f}} ,\\ \hat{P}_{t|{t-1}}&=&\frac{1}{m-1}\sum_{i=1}^m({x}_{t,i,{\rm f}}-{x}_{t,{\rm f}})({x}_{t,i,{\rm f}}-{x}_{t,{\rm f}})' . \end{eqnarray*}
(4) If the inequality below holds, go to step (6); otherwise, go to step (7): $$ ({y}_t-{H}_t{x}_{t,{\rm f}})'(\lambda_t{H}\hat{P}_{t\vert t-1}{H}'+R)^{-1}({y}_t-{H}_t{x}_{t,{\rm f}})\ge L , $$
(5) If Λt<K, we increase the value until it enters the feasible region of step (5), where K represents the upper bound, and we set it to 100 in the simulations; then go to step (5).
(6) Generate the observation samples γ1,γ2,…,γm from N(0,Rt) and evaluate the analysis ensemble, analysis state and its variance by \begin{eqnarray*} {x}_{t,i,{\rm a}}&=&{x}_{t,{\rm f}}+\lambda_t\hat{P}_{t|{t-1}}{H}'_t(\lambda_t{H}_t\hat{P}_{t|{t-1}}{H}'_t+R_t)^{-1} ({y}_t+{\gamma}_i-{H}_t{x}_{t,i,{\rm f}}) ,\\ && 1\le i\le m ,\\ {x}_{t,{\rm a}}&=&\frac{1}{m}\sum_{i=1}^m{x}_{t,i,{\rm a}} .\\ P_{t,{\rm a}}&=&(\lambda_t^{-1}\hat{P}_{t|{t-1}}^{-1}+{H}'_t R_t^{-1}{H}_t)^{-1} . \end{eqnarray*}
(7) If t reaches the last time point, stop the algorithm; otherwise, go to step (4).
Then, the analysis state of all the observation time can be obtained.
2
4.1. Simulation results of the Lorenz-63 model
The Lorenz-63 model is a set of nonlinear differential equations with three variables (Lorenz, 1963). The solution of this model has chaotic behavior and is very sensitive to the initial condition settings. Moreover, this model has been examined by various data assimilation methods for their potential applications with other strongly nonlinear and chaotic models, such as oceanic and atmospheric models (Palmer, 1993; Anderson and Anderson, 1999; Evensen, 2009; Sheng and Li, 2015).4.1.1 Forecast model
First, the Lorenz-63 model can be expressed as \begin{equation} M(x)=\left\{ \begin{array}{l} \dfrac{dx_1}{dt}=ax_2-ax_1\\ \dfrac{dx_2}{dt}=(b-x_3)x_1-x_2\\ \dfrac{dx_3}{dt}=x_1x_2-cx_3 \end{array} \right., \ \ (15)\end{equation} where x=(x1,x2,x3) is the state and a,b,c is the model parameter. Here, we adopt the parameter settings of (Evensen, 2009) and set a=10, b=28, c=8/3.
We use a forecast model with a model error added on model (15): \begin{equation} \label{eq13} {x}_t={M}({x}_{t-1})+{\eta}_t , \ \ (16)\end{equation} where {ηt} is a sequence of independent and identical random variables that obey a normal distribution of N(0,Qt).
4.1.2. Observation model
Here, we consider an observation operator whose dimension is less than three: \begin{equation} \label{eq14} {y}_t={H}_t{x}_t+{\varepsilon}_t , \ \ (17)\end{equation} where {εt} is a sequence of independent and identical random variables that obey a normal distribution of N(0,Rt). The observation matrix H is set as $$ H_t=\left( \begin{array}{c@{\quad}c@{\quad}c} 1 & 2 & 3\\ 1 & 1 & 1 \end{array} \right) , $$ which is independent of time.
4.1.3. Description of the setup
In this simulation experiment, we set the time step as 0.05 and solve the Lorenz model (15) with a fourth-order Runge-Kutta integration scheme, and the observations are obtained every four time steps. The initial condition is given by x0=(1,2,3)', and the true state series are generated by Eq. (14) and denoted by x0,x0.2,x0.4,?,x30.
The observations are generated by Eq. (17) and denoted by y0.2,y0.4,?,y30. We repeat this procedure 200 times and at each time we can obtain an estimate sequence. Hence, the root-mean-square error (RMSE) of the state in the kth dimension of state at time t can be calculated by \begin{equation} \label{eq15} {\rm RMSE}_k(t)=\sqrt{\frac{1}{n}\sum_{i=1}^n(x_{k,t}-x_{k,t,i,{\rm f}})^2} , \ \ (18)\end{equation} where xk,t is the kth dimension of state at time t; and xk,t,i, f is its ith forecast ensample estimate, where i=1,…,200. For simplicity, we denote RMSEk(t) simply as "RMSE". It is obvious that a smaller RMSE usually indicates a better performance of the corresponding assimilation method. Moreover, we also use the time-averaged RMSE value as a criterion to compare the performance of different assimilation methods.
4.1.4. Initial condition settings
Since in practical situations the true initial distribution is usually not available, we set an inaccurate initial condition for all the assimilation methods in the simulation experiments to study the robustness of the new method. That is, in the numerical experiments of subsection 4.1.5 and subsection 4.1.6, the mean value of the initial distribution used to estimate the state diverges a distance of 10 from the true initial value used to generate the state. The initial distribution is set with a normal distribution with mean (11,12,13) and variance 0.52I3. The ensemble size is 30 in all cases. Additionally, in subsection 4.1.7, the influence of different initial condition settings is further studied via numerical experiments.
4.1.5. Assimilation results of one case
First, we compare the performance of EnCR with the W-B and SLS methods in a situation when the model error variance Qt=0.012I3 and the observation error variance Rt=I2.
![](http://www.iapjournals.ac.cn/fileDQKXJZ/journal/article/dqkxjz/2018/3/PIC/0256-1530-35-3-334/170045fig1_mini.jpg)
Figure 1 displays the results of the RMSE of the three dimensions of the state based on the W-B, SLS and EnCR methods, separately. From this figure, we can see that the convergence rate of the EnCR method is slightly faster than that of the other two methods, and the RMSE values of the SLS and EnCR methods are clearly smaller than that of the W-B method. To further compare the performances of the EnCR and SLS methods, we remove the RMSE line of the W-B method and obtain Fig. 2. From this figure, it is clear that the EnCR method obtains the smallest RMSE overall.
![](http://www.iapjournals.ac.cn/fileDQKXJZ/journal/article/dqkxjz/2018/3/PIC/0256-1530-35-3-334/170045fig2_mini.jpg)
Moreover, Table 1 demonstrates the time-averaged RMSE values of the EnKF, W-B, SLS and EnCR methods. This table indicates that the W-B, SLS and EnCR methods significantly reduce the time-averaged RMSE of the EnKF method, and the EnCR method achieves the smallest time-averaged RMSE of all the methods. Furthermore, to study the implications of different variance settings, the four methods are compared under various variance settings in the next subsection.
4.1.6. Implications of the model and observation error variance settings
In this subsection, we study the estimated accuracy and stability of the EnCR method under different error variance settings.
First, the observation error variance is set as Rt=0.12I3, and the value of model error variance is set as Qt=σ12I3, where σ1 is a varying parameter. The results are demonstrated in Table 2. From this table, we can see that as σ1 increases, the estimated precision of all the methods decreases, which is consistent with our intuition. When Qt≤ 0.12I2, the time-averaged RMSEs of the SLS and EnCR methods are quite close; and when Qt>0.12I2, the EnCR method obtains the smallest time-averaged RMSE. Moreover, when the model error variance is large, i.e. the last several columns of Table 2, the time-averaged RMSE of the W-B, SLS and EnCR methods are much smaller than that of the EnKF method, and the EnCR method maintains the smallest time-averaged RMSE value, which suggests that the estimated accuracy of the EnCR method is better than that of the other methods.
Using the information provided in Table 3, we study the influence of the observation error variance when the model error variance is fixed. That is, we set the model error variance Qt=0.12I3 and the observation error variance Rt=σ22I3 with different σ2 values. From this table, we can see that, overall, the estimated precision of the EnCR method is better than that of the other methods. Moreover, it is also notable that, even when the observation variance is relatively large, the new method performs better than the other three methods. Overall, Table 3 shows that the performance of the new method is significantly superior to the other three methods.
Based on the above results, the estimation results of the EnCR method are superior to the other three methods, particularly when the model error variance is large. Furthermore, as the error variance increases, the estimation precision of all the methods decreases. In short, the simulation results indicate that the new method has a higher estimation precision and is more stable than the other methods.
4.1.7. Implications of the initial condition settings
In this section, the implications of the initial condition settings are examined for different assimilation methods. Here, we set the model error variance Qt to 0.012I3 and the observation error variance to Rt=I2. For description simplicity, the variance of the initial condition is set to 0.5I3, and we set the same deviation for the simulation initial value of the three state components to the true initial value and denote it by Diff——that is, \(\rm Diff=\hat{x}_0,i-x_0,i,i=1,2,3\), where \(\hat{x}_0,i\) and x0,i are the simulation initial state and the true state of x0 in the ith component, respectively.
Table 4 gives the estimation results of all the methods under different deviation Diff settings. The time-averaged RMSE of the EnKF method increases significantly as the deviation increases. However, the time-averaged RMSE of the W-B, SLS and EnCR methods grows very slowly as the bias of the model initial condition increases, which indicates that the three methods are not extremely sensitive to the initial condition settings. The estimated precision of the EnCR method is evidently superior to the other methods in this experiment.
It is also interesting to study the estimated inflation values of the three methods with different initial condition settings. As we know, the initial state is biased and the variance inflation technique should be adopted in the first several steps. Here, we consider the mean value of the inflation parameter estimation in the first six steps. Table 5 gives the mean value of the inflation parameter with an increase in the bias of the initial state settings. From this table, we can see that the estimation of the inflation parameter of the new method is more reasonable than the other two methods. That is, when the deviation of the initial state is zero, a smaller inflation should be adopted, and with an increase in the deviation, the inflation parameter should also increase. The inflation estimation of the W-B method is always very large, and the variation of the SLS method is not as reasonable as the new method.
2
4.2. Simulation results of the Lorenz-96 model
The Lorenz-96 model is a 40-dimension differential equations model, and it is widely used in various data assimilation studies. In this section, we use this model to study the performance of the EnCR method.4.2.1. Model description and parameter settings
The Lorenz-96 model can be expressed by \begin{equation} \dfrac{dX_k}{dt}=(X_{k+1}-X_{k+2})X_{k-1}-X_k+F , \ \ (19)\end{equation} where k=1,2,…,K(K=40) and the boundary conditions are assumed to be cyclic; that is, X-1=XK-1, X0=XK, XK+1=X1. Here, F is the external forcing term and we set it to F=8 to generate the true state values in this experiment. The solution of the Lorenz-96 model has chaotic behavior and mimics the temporal evolution of a scalar meteorological quantity on a circle of latitude, and the three terms of the right-hand side of Eq. (19) can be viewed as an advection-like term, a damping term and an external forcing term, respectively (Wu et al., 2013). Besides, similar to the Lorenz63 model, we use the fourth-order Runge-Kutta time integration scheme to solve the state model (19), and set a time step of 0.05 non-dimensional units to drive the true state. Besides, assuming the characteristic time scale of the dissipation in the atmosphere is five days, the time step here is roughly equivalent to six hours in real time (Lorenz, 1996).
![](http://www.iapjournals.ac.cn/fileDQKXJZ/journal/article/dqkxjz/2018/3/PIC/0256-1530-35-3-334/170045fig3_mini.jpg)
![](http://www.iapjournals.ac.cn/fileDQKXJZ/journal/article/dqkxjz/2018/3/PIC/0256-1530-35-3-334/170045fig4_mini.jpg)
In this model, we adopt the parameter settings of (Wu et al., 2013) and set X0,k=F, \(k\neq 20\), X0.20=1.01F, with the time step set as 0.2. To achieve stationary estimation results, we obtain the observations every four time steps over a duration of 100 000.
In this study, the observations are generated by an identical observation operator with a Gaussian distribution error: \begin{equation} y_t=x_t+\varepsilon_t,\varepsilon_t\sim N(0,R_t), \ \ (20)\end{equation} where the observation error variance Rt is spatially correlated with \begin{equation} \label{eq16} R_t(i,j)=\sigma_0^20.5^{\min\{|i-j|,40-|i-j|\}} , \ \ (21)\end{equation} in which we set σ0=1.
In this study, we use the average state estimate error to compare the performance of different methods, and denote it as A-RMSE: \begin{equation} \label{eq17} \mbox{A-RMSE}(t)=\sqrt{\frac{1}{K}\sum_{k=1}^K(\hat{x}_{k,t,{\rm f}}-x_{k,t})^2} , \ \ (22)\end{equation} where K is the dimension of state, xk,t is the true state of the kth dimension at time t, and \(\hat{x}_k,t,\rm f\) is the corresponding forecast ensemble estimate. Similar to the RMSE criterion, a smaller A-RMSE usually indicates a better performance of the corresponding assimilation method.
4.2.2. Results when F is correctly specified
First, we consider a situation when the forcing term F is correctly specified——that is, there has no model error when estimating the state. Here, we use a normal distribution to generate the initial ensemble: \begin{equation} x_{0,j}=x_0+\gamma_j,\gamma_j\sim N(0.05^2I_{40}),1\leq j\leq 30. \ \ (23)\end{equation}
Figure 3 shows the results of the A-RMSE value of the four methods when the sample size is 20. To make the results clearer, only the A-RMSE results of the first 100 s (500 steps) are displayed in the figure (the trend of A-RMSE after 100 s is similar to the trend around the 100 s). When the sample size is small, the EnCR method obtains the best estimation results, followed by the W-B method. However, the results of the SLS and EnKF methods to a certain degree show up the phenomenon of filter degeneracy in this case.
Figure 4 presents the results of the A-RMSE value when the sample size is 80. In this case, the EnKF method relieves the phenomenon of filter degeneracy in the initial stage. However, over time, the EnKF method begins to degenerate at approximately 15 s. Overall, the other three methods prevent the filter degeneration phenomenon quite well, and the estimated accuracy of the SLS and EnCR methods is slightly better than that of the W-B method.
Additionally, Table 6 provides the results of the time-averaged A-RMSE over 100 000 time steps for all the methods under different sample size settings. Overall, the A-RMSE of the EnCR method is smaller and more stable than that of the other methods. Moreover, the estimation results of the EnKF method are the poorest, which is coincident with our intuition.
4.2.3. Results when F is incorrectly specified
Since the true model is usually not available, a model used in practical situations is always an approximation of the real one. Therefore, studying the estimation performance of a data assimilation method in the situation when the model is incorrectly specified is meaningful. Here, we also use the A-RMSE criterion to examine the estimation robustness of different models.
We also adopt the model parameter settings in subsection 4.2.2 and set F=8 to generate the true states; however, we use different values of F when we estimate the states, such as F=4,5,…,12.
Figure 5 shows the results of the time-averaged A-RMSE over 100 000 time steps of the four methods with different forcing term values when the sample size is 20. When the sample size is small, the time-averaged A-RMSE of the EnCR method is smallest among all methods. Moreover, the time-averaged A-RMSE values of the SLS and EnKF methods are significantly larger than those of the other two methods. It is also very interesting that the minimum value of the time-averaged A-RMSE of the SLS and EnKF methods is not achieved at the point when F=8, in which the value of F is same with that of the model when we generate the data. This is a little inconsistent with our intuition and is possibly because when the sample size is very small, both methods are degenerated (as shown in Fig. 3), and thus the observation information contributes very little to the state estimation process. Next, we verify this conjecture by increasing the sample size in the following simulations.
![](http://www.iapjournals.ac.cn/fileDQKXJZ/journal/article/dqkxjz/2018/3/PIC/0256-1530-35-3-334/170045fig5_mini.jpg)
Figures 6-7 show the results when the sample size is 80 and 100, respectively. First, we can see that when the sample size is large, the EnKF and SLS methods achieve their minimum time-averaged A-RMSE when F=8, which verifies the above conjecture. When F is around the true value, there has very little difference in estimation among the SLS, W-B and EnCR methods. However, when F is much smaller than 8, the W-B method achieves the smallest time-averaged A-RMSE, but the results of the EnCR and W-B methods are almost indistinguishable. When F is much larger than 8, the EnCR method achieves the best precision among all the methods. Based on the above results, we can conclude that the EnCR method demonstrates preferable performance and robust estimation results compared with the other methods.
![](http://www.iapjournals.ac.cn/fileDQKXJZ/journal/article/dqkxjz/2018/3/PIC/0256-1530-35-3-334/170045fig6_mini.jpg)
![](http://www.iapjournals.ac.cn/fileDQKXJZ/journal/article/dqkxjz/2018/3/PIC/0256-1530-35-3-334/170045fig7_mini.jpg)
2
4.3. Summary of the experimental results
From the experiments using the Lorenz-63 model, we can draw the following conclusions:(1) Overall, the EnCR method performs best. The SLS method ranks second and, when the observation error variance is small, the precision values of the SLS and EnCR methods are very close. However, when the observation error variance is large, the performance of the EnCR method is clearly better than the other methods.
(2) With the model error variance increases, the precision of the EnCR and SLS methods tends to be similar and better than that of the W-B and EnKF methods.
(3) The EnCR, W-B and SLS methods are not very sensitive to the initial condition settings, and the estimation results of the EnCR method are better than those of the other two methods under different initial deviation settings, which suggests the new method produces a more robust and highly precise estimation result.
From the experiments using the Lorenz-96 model, we can draw the following conclusions:
(1) When the model forcing parameter is correctly specified and the sample size is small, the EnCR method is better than the other methods. With the sample size increases, the precision values of the W-B, SLS and EnCR methods become very close to one another.
(2) When the model forcing parameter is incorrectly specified, the EnCR method is better than the other methods when the sample size is small. When F used in the model is smaller than the true value and the sample size is large, the W-B method obtains the best results, and the results of the EnCR method are very close to the results of the W-B method in that case. Moreover, when F is larger than the true value, the EnCR method achieves the best results and is clearly better than the other methods.
Using the idea of a confidence region to estimate the parameter of the EnKF method can also be extended to other unknown parameter estimation problems within data assimilation methods. In our future studies, the idea of a confidence region will be further studied in a similar direction.