A modified iterative ensemble Kalman filter data assimilation method

High nonlinearity is a typical characteristic associated with data assimilation systems. Additionally, iterative ensemble based methods have attracted a large amount of research attention, which has been focused on dealing with nonlinearity problems. To solve the local convergence problem of the iterative ensemble Kalman filter, a modified iterative ensemble Kalman filter algorithm was put forward, which was based on a global convergence strategy from the perspective of a Gauss-Newton iteration. Through self-adaption, the step factor was adjusted to enable every iteration to approach expected values during the process of the data assimilation. A sensitivity experiment was carried out in a low dimensional Lorenz-63 chaotic system, as well as a Lorenz-96 model. The new method was tested via ensemble size, observation variance, and inflation factor changes, along with other aspects. Meanwhile, comparative research was conducted with both a traditional ensemble Kalman filter and an iterative ensemble Kalman filter. The results showed that the modified iterative ensemble Kalman filter algorithm was a data assimilation method that was able to effectively estimate a strongly nonlinear system state.


Introduction
As in highly nonlinear dynamical and observational models, introducing iterations to a Kalman filter has proven to be beneficial for estimation problems [7] . To improve the performance of the traditional ensemble Kalman filter (EnKF) [5] , the iterative ensemble Kalman filter (IEnKF) was recently proposed [13] . At the cost of additional iterations and propagations of the ensemble, this IEnKF was shown to significantly outperform the EnKF, especially in cases with large time intervals between updates [3] . The IEnKF should be considered for use in cases where there is a well-constrained system with nonlinear growth during the assimilation cycle (the propagation was initially nonlinear but became linear through repeated iteration). Later on, Bocquet and Sakov (2012) proved that the IEnKF can be taken as a lag-one smoother and be further extended to a fixed-lag smoother: the iterative ensemble Kalman smoother (IEnKS). The IEnKS is an ensemble variational method that does not require the use of the tangent linear of the evolution or the adjoint of observation models. By using a nonlinear least square formulation, Gu and Oliver [6] were able to introduce an ensemble randomized maximum likelihood (EnRML) method. In an EnRML, the sensitivities are calculated using an ensemble based method. However, for high-dimensional problems, the ensemble approximation of the sensitivity matrix is often poor [4] . An iterated EnKF [9] was also proposed with a running in place scheme, which was referenced by Sakov [13] and applied by Penny [12] in a realistic Ocean General Circulation Model (OGCM) reanalysis. In further study of all iteration-based Ensemble Kalman filters, we found that one downside was that only when certain local conditions are satisfied can the Because the Newton method, which was implicitly adopted [13] , is only one of many minimization schemes to minimize the cost function of data assimilation systems, there is still a large degree of freedom in choosing an iterative scheme, such as the Levenberg-Marquardt algorithm [3] . In this study, an IEnKF with a quick local convergence is combined with a global convergence strategy for the purpose of achieving a modified iterative ensemble Kalman filter algorithm (modified IEnKF: MIEnKF) that leads to every iteration in the process of assimilation to tend towards the optimal value. Thus, even if the initial value was not accurately estimated, the higher-accuracy solution of the iterative convergence can still be obtained. Experimental tests were conducted on this algorithm with both ensemble Kalman and iterative ensemble Kalman filters in some typical, corresponding chaotic systems to confirm the algorithm.

Gauss-Newton derivation of the IEnKF
A Gauss-Newton iteration algorithm is an approximate application of a Newton iteration algorithm, which is an iterative method type for nonlinear equations. In this study, the updated problem of the EnKF observation was transformed to solve the minimum cost function, and the observation iterative equation was deduced from the angle of a Gauss-Newton iteration algorithm, which thereby proved that an IEnKF is a type of Gauss-Newton iteration process.
The IEnKF was introduced [13] , in which the prior probability density function is (1) where l x indicates the l-th member forecast ensemble. It need not match the mean b x , and the error covariance matrix B of the prior probability density function, where  is the normalized constant and y indicates the observation vector at a certain update step.   | p y x is the likelihood probability and satisfies: where h() is the observation operator and R is the observation error covariance matrix. Therefore, while solving the maximum posteriori estimation of x,  can be ignored, and the cost function of the state variable can be shown as follows: Moreover, when

SY
, then the following can be obtained: Jl  xx (9) According to Eq. (9), the IEnKF update can be transferred to the solution of a nonlinear least squares method, and the Gauss-Newton iterative algorithm can be used to obtain the solution of Eq. (9) as follows: I is the unit matrix. Then, when substituted into Eq. (10), the following can be obtained: When the gain is set as (12) This equation is the updated iterative formula of the iterative ensemble Kalman filter.
It was determined that the observed iteration of a traditional iterative ensemble Kalman filter algorithm was a Gauss-Newton iterative process.
Due to the large spread of initial value deviations, as well as other factors during the process of data assimilation, the most recent iteration value was chosen as the final value in every iteration of the original IEnKF, and the new iteration's result was considered to be closer to the expected value when compared to the former iteration. However, this method may reduce the convergent performance in certain applications [13] .

A Modified IEnKF scheme
The Gauss-Newton iteration algorithm was found to have stronger local convergence. When the initial value was close enough to the optimal value, the convergence of the iterative sequences was ensured. During this process, the Newton step is always used for examining whether it is closer to the expected value by iterations. However, when the initial estimate is quite different from the expected value, the iteration of the Newton step may lead to a greater deviation from the expected value. In this study, through adding a step factor,  , the value of the step factor will adaptively adjust every iteration to ensure that each step of the iteration tends towards the expected. It has been proven that a Gauss-Newton method with a step factor  presented a general convergence [14] .
To determine the step factor, the cost function J (x) was minimized with reference to search step factor i  along the Newton direction, enabling J (x) to be minimized along the Newton direction, and so on, as follows: Through solving Eq. (10), it was found that the multiple step factor  met the conditions, so it could be used to calculate ( ) , which corresponded to each  many times, and the  value was then taken to enable the ( x minimum through comparison. Due to the multiple calculations for the status value, a large amount of calculating was consumed by this process in each iteration. When the data assimilation system showed higher requirements for real time, it would determine the step factor  according to the J (x) gradient descent direction.
In regards to the use of the Newton step, once it approached the expected value, a super-linear convergence rate was achieved. In each iteration, it was necessary to test whether the current step reduced the error or not. If the answer was negative, back-tracking was conducted along the Newton direction, until a step meeting the conditions was located. This step did not necessarily support the J(x) minimum in such a direction, but was the needed step to be able to make it decline. Because the Newton direction was downward for the J(x) direction, an acceptable step needed to be found in each of the back-tracking processes. Then, the following was obtained: In an actual situation, 4 10    should be realized. Meanwhile, according to the specific situation, the step factor should be as large as possible. Through the back-tracking step, the easiest way was to halve the step factor in each of the back-tracking processes. Only when  was approaching 1, could the super-linear convergence rate of the Newton iterative method be obtained.
Based on the above global convergence strategy, it was determined that the Newton iterative algorithm may cause a larger norm of the system output value at the point of the new iteration value, under the local condition of not satisfying Newton's lemma. For example, the new iteration value deviated more from the optimal value, which sometimes resulted in such divergence as Therefore, a Newton iterative algorithm with a global convergence strategy was used to solve this problem. If is taken as the new iteration value 1 i  x . Until the assimilation norm of the new iteration value was less than () i J x , the value satisfying the conditions was assigned to 1 i x to satisfy the local conditions of Newton's lemma in the next assimilation. However, multiple calculations were required for this method, which aggravated the burden of the calculation and affected the convergence speed so that the weighted average method could be used to improve the speed of the convergence. Table 1 illustrates the algorithm: a modified iterative ensemble Kalman filter.  Under certain local conditions, a Gauss-Newton iterative algorithm descends along the Newton direction. In regards to an IEnKF, if the norm of the new iteration assimilation result is greater than that of the previous assimilation result, it has been determined that both assimilation results become distributed on both sides of the optimal value absolutely. For example, when x , the weighted average algorithm can be used. The new iteration value is made to return along the old iteration value direction, where the optimal solution is closer to the true value and can be obtained as the new iteration value. The comparison with IEnKF algorithm is as shown in Figure 1.  Figure 1, the x-axis shows the number of iterations, while the y-axis indicates the state changes. The red dotted box indicates the specific process of using the MIEnKF for reducing the errors at the i+1-th iteration, when the i+1-th iteration value greatly deviates from the i-th iteration value.
According to the gradient descent direction of the cost function, the step  is adaptively adjusted to make assimilation results eventually converge, and so on. So, the MIEnKF can better adapt to a strongly nonlinear system.

Numerical experiments
With different parameterizations and physical representations determining the sensitivity of the procedure of different models, observing system simulation experiments (OSSE) are designed to enable examination of the performance of data assimilation procedures [12] . In this section, the MIEnKF is compared with a standard IEnKF method, as well as a traditional EnKF method for research purposes.

Low-dimensional chaotic system Lorenz-63
A Lorenz-63 model is a three-variable chaotic dynamical system, of which the control equation of a state evolution can be given by the following ordinary differential equation: where x and true x are the estimate and true value of the model state, respectively, and n is the dimension of the model state vector.

Experimental results and analyses
In this section, the assimilation performances of the IEnKF and MIEnKF are examined. For the same setup, the modified effects were found to be different.

Influences of ensemble size.
In the experiment, we changed ensemble sizes from N=5 to N=80; the model step T = 25 was selected, the inflation factor was set as infl = 1.05, and the observation error  analysis-rmse The research showed that: (1) with an increase of ensemble size, the performance of the three algorithms improved in agreement with typical results from the literature [1] . While the EnKF seems to improve with increased ensemble size, the IEnKF and MIEnKF show very little improvement. The larger the ensemble size, the better the assimilation effect. However, the assimilation time was correspondingly longer. (2) Among the three types of algorithms, the performance of the EnKF was found to be the least effective. In cases with a small ensemble size, the rmse was larger. The IEnKF was determined to have a slightly better performance than the EnKF. However, the modified MIEnKF was able to meet the demand of the actual assimilation more effectively, and the rmse was minimal.

Influences of the inflation factor.
An inflation factor is a type of error processing method for improving assimilation effects. If the inflation factor is changed, the performance of the different algorithms will change greatly.
In this study, the results showed the following: (1) When the inflation factor was increased, three types of algorithms had diverse results, wherein the performance of the EnKF improved the most, the improvement of the IEnKF was not obvious, and no influence was made on the MIEnKF.
(2) When the step of the inflation factor was equally reduced, it was found that the EnKF had the worst performance, and the rmse of the IEnKF decreased from 0.237 to 0.236, while no changes were made on the MIEnKF. Therefore, when combined with the performance curve diagram of the three algorithms, as shown in Figure 3, it can be confirmed that no matter what type of situation, the performance of the EnKF was the worst. For example, the rmse value of the EnKF presented multiple peak-trough phenomena as the equal interval of the inflation factor increased or decreased. Therefore, it was particularly important to select the inflation factor. To focus on this problem, a genetic algorithm was used to select the optimal inflation factor from previous studies [1] and to obtain the optimal data assimilation results under limited computing resources.
(3) When compared with the IEnKF, this algorithm was able to adapt well to the data assimilation of a strongly nonlinear system. In addition, when combining the algorithm performance variations caused by ensemble size changes, it could be seen that the rmse of the MIEnKF was smaller, and the same was found for the inflation factor. Therefore, the MIEnKF had a more robust performance than the IEnKF.   In this section, to further detect the performance of the modified iterative ensemble Kalman filter, a Lorenz-96 model that had 40 state variables [11] was chosen. A Lorenz-96 system is a more complex nonlinear dynamic system, and it is used to simulate the time evolution of atmospheric variables. The equation is as follows:

High-dimensional chaotic system
The discrete numerical experiments adopted a fourth-order Runge-Kutta method, where 1, 2, , jM  represent the loop label, forcing parameters of 8  F , and variables of 40  M . The time step dt was taken as a 0.05 dimensionless unit fixed value.
The performances were compared among the three methods (EnKF, IEnKF, and MIEnKF) by changing the observation error and the inflation factors. The ensemble size was set as N = 30, the observation variance as obs_var = 1, and the model step as T = 25. The inflation factor infl changed from 1.01 to 1.20 and was set at intervals of 0.01. Its performance in data assimilation was then observed. Figure 4 illustrates the observation error influence on the performance of the three algorithms' assimilation. In the experiment, the ensemble size was set as N = 30, the model step as T = 25, and inflation factor as infl = 1.05. The observation error, obs_var, ranged between 0.5 and 8, and the step size was 0.5. The results confirmed that the modified algorithm could still finish with an accurate estimation, even when the interval between the observation errors was larger by contrast. In the early stages, the EnKF could basically meet the true state. However, great deviations appeared in later stages. The floating range of the IEnKF was relatively small, and therefore the assimilation performance was relatively poor because the analogue system error was larger. When combined with the previous analysis, the new method was found to have a greatly improved performance compared to the first two algorithms.
In addition, in terms of computational costs, the average time for the EnKF was 18.25 s, while the average time of the IEnKF and MIEnKF were approximately 35.25 s and 44.5 s, respectively. We found that the MIEnKF needed a larger amount of calculation, which took longer; therefore, reducing the computation time while maintaining the premise of satisfying the estimation precision is still a subject for further study.

Conclusions and further discussions
In terms of a linear Gaussian system, the EnKF was the optimal algorithm for data assimilation. However, as the operational-scale models increase in resolution and complexity, nonlinearities pose greater problems to today's most widely used data assimilation methods. These nonlinearities are particularly critical for determining the optimal assimilation method for a strongly nonlinear system. A MIEnKF was proposed in this study for the purpose of exploring narrating an IEnKF from the perspective of a Newton iteration and to determine the optimal iteration value by judging global convergence assimilation results. This iterative version of an EnKF is used to minimize the local (i.e., limited to on assimilation cycle) cost function and investigate the performance of such a scheme in a number of experiments withtwo common nonlinear chaotic models. The only difference between the MIEnKF method and the IEnKF method is the method of minimization: we used a Gauss-Newton derivation of the IEnKF. This is the reason the results of some experiments with the MIEnKF method are almost a straight line. Due to the OSSE setup for such synthetic experiments, "truth" is defined as when the model is integrated (run) for a set of initial conditions. "Observations" are generated from the "truth" by adding some type of noise. Assuming the error from the observations and forecast ensemble are ideal because of the applied minimization method, we can guarantee the assimilation performance is better than other methods. Meanwhile, when compared with the EnKF, the computing burden of the MIEnKF is increased with greater iteration numbers. Therefore, to satisfy the premise of estimation precision, a solution for reducing the calculation burden of the MIEnKF is an important direction for future study.