Differential evolution algorithm-based multiple-factor optimization methods for data assimilation

Abstract

The methods of searching for optimized parameters have substantial effects on the forecast accuracy of ensemble data assimilation systems. The selection of these factors is usually performed using trial-and-error methods, and poor parameterizations may lead to filter divergence. Combined with the local ensemble transform Kalman filtering method (LETKF), a technique for an automated search of the best configuration (parameters) of a data assimilation system is proposed. To obtain better assimilation, a differential evolution (DE) algorithm-based multiple-factor parameterization method results in the corresponding circumstances. By combining with fast-searching DE algorithms, we may retrieve the most ideal parameter combinations. Several numerical experiments performed with the Lorenz-96 model show that new methods performed better than the original one-parameter optimization methods. As the basis of DE methods, the best combinations of the local radius and the covariance inflation parameter, which can guarantee the best DA performances in the corresponding circumstances, are retrieved. It is found that the new method is capable of outperforming previous search algorithms under both perfect and imperfect model scenarios, and the calculation cost in Lorenz-96 model is lower. However, how to apply the new proposed method to more complex atmospheric or land surface models requires further verification.

Keywords

Differential evolution local ensemble transform Kalman filter Lorenz-96 chaos system multiple-factor optimization

1. Introduction

Data assimilation (DA) provides an objective methodology to combine prior information, observational data, numerical model predictions, and their corresponding error statistics in order to produce a better estimate of the state of a geophysical system [1]. To reduce sampling errors and remove spurious long-distance correlations in the ensemble-based error covariance, inflation [2] and localization [3, 4] have been introduced into ensemble Kalman filters (EnKF) to offset the underestimation of the true variance. Moreover, DA involves two key control parameters, i.e., the inflation factor and localization length scale, which may significantly influence the optimization performance of the DA. Therefore, the optimal values of those tuning parameters, which may depend on different cases (e.g., models and ensemble sizes), are traditionally determined using trial-and-error methods. Meanwhile, it is difficult to select appropriate factors in real data assimilation due to the lack of real model state, limited observational information and the missing model and observational errors [1, 5].

A series of techniques have proposed within the DA literature to deal with these tuning factors. For example, a relationship between localization lengths and ensemble sizes in domain localization and observation localization has been proposed by Kirchgessner et al. [6]. In addition, an efficient serial scheme was introduced in the EnKF analysis step for the purpose of mitigating the effects of observations’ perturbations [7]. Moreover, Raanes et al. [8] studied multiplicative inflation with the complementary scaling of the state covariance in the EnKF. In practice, Gharamti [9] proposed an enhanced adaptive inflation algorithm for ensemble filters to combat the loss of variance caused by various models and sampling errors in the forecast. Additionally, Attia and Constantinescu [10] presented an optimal adaptive tuning experimental design framework to minimize the uncertainty during the DA processes. However, there was not a theoretical method presented to select the best parameters, especially for the multiple factors that influence the whole DA process.

Inspired by Genetic Algorithm (GA) theories, a new data assimilation method coupled with the crossover principles of vector algorithms based on ensemble transform Kalman filters (ETKFs) was proposed to solve the difficult problem of error adjustment factor searches [11]. Subsequent to this work, a new parameterization method was developed using evolutionary strategies that were both relatively easy to program and computationally efficient [12]. However, only one parameter could be selected within those frameworks. In this study, inspired by the previous evidence of error handling methods, we extend the previous studies to propose a multiple-factor optimization method with DA to choose better values within the DA framework. The search algorithm belongs to the class of derivative-free evolutionary optimization methods designed to mimic biological (Darwin’s) evolution. In the data analysis part, two or more factors could be selected using DE algorithm-based DA methods, which provide the best assimilation results in certain circumstances.

The remainder of this paper is organized as follows. In Section 2, the classical local ensemble transform Kalman filter (LETKF) is reviewed as the theoretical basis of this study. Then, a DE-based parameter optimization system is introduced in Section 3. Lorenz-96 model is used to verify the new method in Section 4. A summary of the main results followed by a general discussion concludes the paper in Section 5.

2. Related works and theoretical basis

The LETKF was proposed as a more computationally efficient version of the ETKF [3, 4]. In this scheme, the number of ensemble members is $k$ , in which the ensemble perturbation matrix $X$ is defined as follow.

$\displaystyle X=\{{x^{(1)}}-\bar{x},x^{(2)}-\bar{x},\ldots,x^{(k)}-\bar{x}\}$ (1)

where $x^{({i})}$ and $\bar{x}$ are the $i^{\text{th}}$ ensemble state vector and the mean of state estimate, respectively. The analysis ensemble mean $\bar{X}_{p}^{a}$ , is given by the LETKF update [4]:

$\displaystyle\bar{X}_{p}^{a}=\bar{X}_{p}^{b}+X_{p}^{b}\bar{W}^{a}$ (2)

Refer to Eq. (2), $\bar{X}_{p}^{b}$ is the $m_{p}$ vector of the background ensemble mean, where $m_{p}$ represents the number of the estimated parameters; $X_{p}^{b}$ is the $m_{p}\times k$ background ensemble perturbation matrix. The optimal linear-combination weights $({\bar{W}^{a}})$ are given by the following expression:

$\displaystyle\bar{W}^{a}=\widetilde{P}^{a}({Y^{b}})^{T}R^{-1}({y^{o}-\bar{y}^{% b}})$ (3)

where $R$ is the $l\times l$ observation error covariance matrix, $y^{o}$ is the $l$ -dimensional observation vector, $\widetilde{P}^{a}$ is the $k\times k$ analysis error covariance matrix, tilde means that the matrix is in the ensemble subspace, and $Y^{b}$ is an $l\times k$ matrix composed of the perturbations of the background ensemble in the observation space. $\bar{y}^{b}$ is the background ensemble mean in the observation space, i.e., $\bar{y}^{b}=\frac{1}{k}\sum_{i=1}^{k}{H({x^{b(i)}})}$ , where $x^{b(i)}$ denotes the $i^{\text{th}}$ member of the background ensemble, and $H$ is the observation operator. To calculate the columns of the $Y^{b}$ matrix, we can apply $H$ to each background ensemble member $x^{b(i)}$ and from each member $H({x^{b(i)}})$ subtract the ensemble mean $\bar{y}^{b}$ to obtain the background ensemble perturbations in the observation space. It is worth noting that these weights contain the information about state variables in the observation space, as well as error correlation between the parameters and state variables.

The estimated parameter disturbance is obtained by Eq. (4)

$\displaystyle X_{p}^{a}=X_{p}^{b}W^{a}$ (4)

Where $W^{a}$ is the transformation matrix, which is defined in Eq. (5)

$\displaystyle W^{a}=[{({k-1})\widetilde{P}^{a}}]^{\frac{1}{2}}$ (5)

The analysis covariance matrix $\widetilde{P}^{a}$ of parameters is computed by the observation error covariance matrix which has no localization, as we will be discuss further later. After each data assimilation step, we can use the optimal parameters and the analyzed state for the next step.

3. A DE-based parameter optimized DA system

In the LETKF-based data assimilation system, many factors affect the assimilation results. The performance of DA systems is largely based on the optimal choice of the key parameters of the DA methods, such as the covariance inflation factors, the localized radius and the ensemble numbers. Unfortunately, the optimization methods that have been used to solve static data assimilation problems face difficulties in direct dynamic DA applications. The main reason is that the traditional GAs usually converge to a single optimum, while in dynamic DA problems, the solution landscape continues to change, which means that there are several local optimal points that will be revealed during the solving process. One of the ways to overcome the shortcomings of traditional genetic algorithm-based methods is to apply a differential evolution algorithm. Figure 1 illustrates the simultaneous estimation cycle.

Figure 1.

Schematic representation of the optimal inflation factors and localization radius searching DA cycles using the Differential Evolution (DE) algorithm.

Similar to the DA systems [11], Fig. 2 shows the schematic diagram of the DA systems coupled with the DE algorithms. The idea behind this approach and the connection between the DA models and the DE is shown in Fig. 1. The DA systems in Fig. 2 will be the same as those in traditional DA systems. However, the fitness minimum mRMSE function defined later will be selected based on the real DA systems. The DE part of Fig. 2 will handle to the select the multiple optimal factors, such as the combination (referred to as the “gene” in the DE algorithm) of the local radius parameter and the covariance inflation parameter. In terms of the coupling strategy, our method is similar to those presented in previous studies [14]. The choice of the fitness function can refer to those in the literature for different DA systems.

Figure 2.

Flowchart of a DE-based parameter optimized DA system. The DE search part is connected to the normal DA part serially, which forms a kind of feedback mechanism of the DA systems.

The DE algorithm is a population-based, heuristic global search algorithm [15]. Like other evolutionary algorithms, the technique proposed in this paper defines a “population” of parameter vectors (“individuals”). “Individuals” in the “population” give “birth” to a new “generation” whose “individuals” (“children”) are computed by a number of simple rules called “mutation”, “crossover”, etc. From the “parents”, “individuals” in the new “generation” are evaluated by a “fitness function” (the time and space averaged forecast RMSE). “Individuals” with the smallest fitness function “survive” and “give birth” to the next “generation”.

3.1 Initialization

The population for which the data assimilation method parameters are considered as vectors is initialized. The original population is defined as follow:

$\displaystyle G_{i,t}(i=1,\ldots,NP;t=1,\ldots,\text{maxgen})$ (6)

Here, $N P$ is the population size and “maxgen” means the maximum generation. In this scheme, employing a trial-and-error scheme to search for the most suitable population parameter obtains a $N P$ value of 20. In the following experiment, we will choose the local radius as one vector and the covariance inflation parameter as the other. Therefore, the number of vectors is two, which is the component of “gene” in the DE. Therefore, the data assimilation method parameter vectors in $G_{i,t}$ are $g_{i,t}^{j}$ ( $j\in[{1,2}]$ ; $i=1,\ldots,NP$ ). The vectors of the individuals $G_{i,t}$ are represented by Eqs (7) and (8):

$\displaystyle G_{i,t}=\{{g_{i,t}^{j}}\}$ (7) $\displaystyle g_{i,t}^{j}=({g_{\min,t}^{j}+\textit{rand}({0,1})*({g_{\max,t}^{% j}-g_{\min,t}^{j}})})$ (8)

where rand (0, 1) denotes the uniformly distributed random variables within the range (0, 1). The prescribed minimum and maximum parameter bounds are:

$\displaystyle G_{\min}=\{{g_{\min,1}^{j},\ldots,g_{\min,\text{maxgen}}^{j}}\},% j\in[{1,2}]\ \ \ G_{\max}=\{{g_{\max,1}^{j},\ldots,g_{\max,\text{maxgen}}^{j}}% \},j\in[{1,2}].$

3.2 Mutation operation

The difference between the differential evolution algorithm and the ordinary GA is that the former employs a difference strategy to achieve the variation. The main idea of the difference strategy is to randomly select two different individuals in the population and take the adjustment value of the individuals’ difference with a third individual to obtain the trial vector according to Eq. (9):

$\displaystyle W_{i,t}=G_{r_{1},t}+B*(G_{r_{3},t}-G_{r_{2},t}),r_{1},r_{2},r_{3% }\in[{1,\ldots,NP}]\text{ and }r_{1}\neq r_{2}\neq r_{3}$ (9)

where $G_{r_{1},t}$ denotes the $r_{1}^{\text{th}}$ individual in the $t^{\text{th}}$ generation. The scaling factor $B$ represents a positive control parameter for scaling the difference vector within the range (0, 2). $r_{1}$ , $r_{2}$ and $r_{3}$ are randomly generated from the matrix $[1,2,\ldots,NP]$ and are different from each other. The $i^{\text{th}}$ population $G_{i,t}(i=1,\ldots,NP)$ generates the trial vector $W_{i,t}(i=1,\ldots,NP)$ via this mutation.

3.3 Crossover operation

The crossover operation is used to reduce the overly fast convergence of a vector perturbation. The population $V_{i,t}(i=1,\ldots,NP)$ is the result of the crossover of the $t^{\text{th}}$ generation $G_{i,t}(i=1,\ldots,NP)$ and the trial vector $W_{i,t}(i=1,\ldots,NP)$ . In the basic version, DE employs the binomial (uniform) crossover, which is defined as follows:

$\displaystyle V_{i,t}^{j}=\left\{{\begin{array}[]{ll}W_{i,t}^{j}&(\text{if }% \textit{rand}(0,1)\leqslant CR\text{ or }{j=j_{\textit{rand}}})\\ G_{i,t}^{j}&\text{otherwise}\\ \end{array}}\right.$ (10)

In Eq. (10), the crossover rate $C R$ is a constant within the range [0, 1] specified by the user ,which controls the fraction of parameter values copied from the mutant vector. $G_{i,t}^{j}$ denotes the $j^{\text{th}}$ vector value of the $i^{\text{th}}$ individual of the $t^{\text{th}}$ generation. $j_{\textit{rand}}$ is a random number that is chosen in the vector $[1,2,\ldots,NP]$ . The binomial crossover operation can hold at least one parameter from $W_{i,t}^{j}$ to the corresponding element in the trial vector $V_{i,t}^{j}$ , if $\textit{rand}(0,1)\leqslant CR$ or $j=j_{\textit{rand}}$ .

3.4 Selection operation

The selection is the decisive operation in the DE algorithm that changes the population composition. The standard used to determine whether the generation population individual $V_{i,t}(i=1,\ldots,NP)$ is an individual of the next generation population is based on which individual has the smaller mRMSE, defined by Eq. (11).

The mRMSE of time is used to evaluate the performance of the DA system. The mRMSE is computed according to Eq. (11):

$\displaystyle\textit{mRMSE}=\frac{1}{t_{\max}}\sum_{k=1}^{t_{\max}}{\sqrt{% \frac{1}{n}\sum_{i=1}^{n}{(x_{k,i}-\bar{x}_{k,i})^{2}}}}$ (11)

Where $t_{\max}$ is the total time. $x_{k,i}$ is the state estimate value of the $i^{\text{th}}$ dimension. $\bar{x}_{k,i}$ denotes the mean value of the $i^{\text{th}}$ state estimate, and $n$ is the ensemble numbers.

The selection standard is defined by Eq. (12):

$\displaystyle G_{i,t+1}^{j}=\left\{\begin{array}[]{l}V_{i,t}(\text{if }\textit% {mRMSE}(V_{i,t})\leqslant\textit{mRMSE}(G_{i,t}))\\ G_{i,t}(\text{if }\textit{mRMSE}(V_{i,t})>\textit{mRMSE}(G_{i,t}))\\ \end{array}\right.$ (12)

The $t+1$ generation population is obtained. Along with the increase in the generation, the optimal parameter combination in the parameter space of the DA system method can been determined.

In addition, because of the bounds of the parameter space of the DA system, there may be a cross-border situation in the mutation process, resulting in meaningless system parameters. In this paper, the evolutionary optimization process is continued by using a random value in the parameter space instead of the crossed value.

3.5 The algorithm of DE with the choice of B and CR

The performance of the DE algorithm depends mainly on two aspects, namely the selected test vector generation strategy and the relevant parameter values used. In the past decade, many empirical guidelines have been proposed for choosing trial vector generation strategies and their associated control parameter settings [15, 16, 17]. In our case, according to Storn [15], we choose a reasonable value for the crossover operation $C R$ of 0.9, and set the good initial choice of $B$ to be 0.85. The algorithmic description of DE is summarized in Algorithm 1.

Algorithm 1: Algorithmic description of DE
1: Generate a uniformly distributed random initial population of $N P$ solutions according to Eqs (7) and (8)
2: while terminal condition is not satisfied
3: for $i=1 \textbf{to} NP$
4: Generate three random indexes $r_{1}$ , $r_{2}$ and $r_{3}$ with $r_{1}\neq r_{2}\neq r_{3}$ //mutation
5: $W_{i,t}=G_{r_{1},t}+B*(G_{r_{2},t}-G_{r_{3},t})$ //end mutation
6: $j_{\textit{rand}}=\textit{randind}({1,2})$ //crossover
7: for $i=1\textbf{ to }NP$
8: if $\textit{rand}({({0,1})\leqslant CR \text{or} j=j_{\textit{rand}}})$
9: $V_{i,t}^{j}=W_{i,t}^{j}$
10: else
11: $V_{i,t}^{j}=G_{i,t}^{j}$
12: end if
13: end for //end crossover
14: if $\textit{mRMSE}(V_{i,t})\leqslant\textit{mRMSE}(G_{i,t})$ //selection (The selection is based on which individual has the smaller
mRMSE)
15: $G_{i,t+1}^{j}=V_{i,t}$
16: else
17: $G_{i,t+1}^{j}=G_{i,t}$
18: end if //end selection
19: end for
20: end while

4. Numerical experiments

The Lorenz-96 model is a strongly nonlinear dynamical system and is used to verify the validity of the DA algorithm [18]. The state dimension is chosen as 40 and the governing equations are given in the following:

$\displaystyle\frac{dx_{i}}{dt}=-x_{i-2}x_{i-1}+x_{i-1}x_{i+1}-x_{i}+F$ (13)

The main components of the experimental system are a high-dimensional nonlinear Lorenz-96 model, the LETKF and the DE algorithm. In this experiment, the optimal parameter combination of the local radius (irad) parameter and the covariance inflation (infla) parameter are sought, and the corresponding mRMSE is the minimum value in the parameter space. In the following experiments, for a normal DA system parameter setting, the ensemble number chosen is 10. The total number of iteration steps is 10000. To avoid the transient effect, the first 1000 steps are omitted. The observation error covariance is defined as 2. Usually atmospheric data are dense in the Northern Hemisphere (over United States, Europe, …) and very sparse in places like the Pacific. A good strategy would be to observe the first 20 variables and see how the optimization algorithm work. The observation numbers are set to 20 and are the fraction of one point at which observations are taken.

4.1 The characteristics of two factors in the DA process

To check the trend of the mRMSE with changing parameters in the DA process, the mRMSE for LETKF in the parameter space of the local radius and covariance inflation is calculated. In the experiment, we set the range of covariance error inflation and localized radius as follows:

$\displaystyle\textit{infla}\in[1.00,\ldots,1.20]$ $\displaystyle\textit{irad}\in[1,\ldots,15]$

From Fig. 3, we can draw several conclusions. (1) With the increase in the local radius (irad), the mRMSE of LETKF and the growth rate of mRMSE become larger. (2) The experiments with severe model errors (Fig. 3a–d) produce very similar results to those of $F=$ 8. However, more irregular changes are present in the presence of moderate model errors ( $F=$ 3; $F=$ 9; $F=$ 13) than in the normal case ( $F=$ 8). (3) Due to the irregular changes of mRMSE in the parameter space of the local radius and covariance inflation, the minimal mRMSE and the optimal parameter combination cannot be obtained by direct judgments. Using a DE algorithm with a global search capability is a new way to find the optimal parameter combination for the minimal mRMSE.

Figure 3.

Time mRMSE of inflation and local radius vectors. (a)–(d) with severe model error and (b) without model error.

4.2 The DE optimization process

Using the local radius and the covariance inflation parameter as the gene in the DE framework, the possibility of searching for the optimal parameter combination is investigated. The individual distribution of the DA method parameter population with the population update is studied. In the experiment, the number of populations is set to be 15 and the total number of generations is set to be 20. The forcing factor $F=$ 3 is considered as an experiment with severe model error.

Figure 4.

The individual distributions with the evolutions with severe model error ( $F=$ 3) and the populations of the 1 ${}^{\text{st}}$ , 4 ${}^{\text{th}}$ , 8 ${}^{\text{th}}$ , 12 ${}^{\text{th}}$ , 16 ${}^{\text{th}}$ and 20 ${}^{\text{th}}$ iterations. The red points represent the individual “genes” in the figure. The final selection shown in the last figure is the best individual, and the blue point is the random choice for poor selection.

Figure 4 shows the individual distributions of the initial population and the populations of the 1st, 4th, 8th, 12th, 16th and 20th generations. The red points in Fig. 4 represent the individuals “genes”. The individual distribution of the initial population is random. (1) With the increase in the generation of the DE algorithm, the individuals with more minimal mRMSE values are preserved, such that the individuals in Fig. 4 gradually gather to the optimal individual. (2) The higher the population generation is, the more frequently the same vectors appear. This results in an increase in the number of repetitive individuals appearing in Fig. 4 and a decrease in the number of red points in Fig. 4. (3) In the 12th generation, the optimal individual appears. In the 20th generation, the vectors of all the individuals are the same, verifying that this is the optimal individual in the population. The vectors of the optimal individual, which is the optimal value of the local radius ( $\textit{irad}=$ 2) and the covariance inflation parameter ( $\textit{infla}=$ 1.14), are the optimal parameter combinations of this DA method. Further experiments will prove this, as shown in the following section.

Figure 5.

The individual distributions with the evolutions with severe model error ( $F=$ 9) and the populations of the 1 ${}^{\text{st}}$ , 4 ${}^{\text{th}}$ , 8 ${}^{\text{th}}$ , 12 ${}^{\text{th}}$ , 16 ${}^{\text{th}}$ and 20 ${}^{\text{th}}$ iterations. Other illustrations are the same as those in Fig. 3. After the 16 ${}^{\text{th}}$ iteration, the best individuals are more clearly displayed.

In Fig. 5, the experimental conditions are set; the number of the population is 15, and the number of total generations is set to 20. The observation numbers are set to 20. The forcing factor $F=$ 9 is considered as an experiment with severe model error. In the 12th generation, the vectors of all the individuals are the same, verifying that this is the optimal individual in the population. The vectors of the optimal individual, which are the optimal value of the local radius ( $\textit{irad}=$ 2) and the covariance inflation parameter ( $\textit{infla}=$ 1.06), are the optimal parameter combinations of this DA method.

4.3 The convergence of the new algorithms

Furthermore, the convergence of the minimal mRMSE values in the population with increasing generations is studied. The number of the population is set to 15 and the number of the total generations to 20.

Figure 6.

The convergence of the DE-based parameter-optimized DA system. (a) without model error ( $F=$ 8) and (b) with severe model error ( $F=$ 9).

As shown in Fig. 6, with the increase in the population generation, the optimal individual in the current population becomes gradually closer to the optimal individual in the global parameter space. At the same time, when the number of generations reaches a certain value, the optimal individual is selected and will not fluctuate again, which indicates that the method in this paper has better convergence. Although each experiment has different waveforms, all can show a good convergence. Notably, in the presence of severe model error, the new DE algorithms perform moderately better than the same cases without model error.

4.4 The differential values of the vectors

In the DE algorithm, the vector difference values between the vectors reflect the differences between the population individuals. The number of populations is set to 15 and the number of the total generations to 20.

In the process of the optimization of the parameters of the DA system method using the DE algorithm, the differential values of the vectors reflect the degrees of individual differences. Figure 7 shows that with the increase in the population generations, the maximum values of the difference are gradually reduced, and the global distribution is closer to the zero value in Fig. 7. This verifies that the search using the method of the global region shifts to the local region. Finally, the exact combination of the specific vectors corresponding to the optimal individual (minimal mRMSE) is sought, and in the subsequent population generations, the vector differential distributions no longer change.

Figure 7.

The vectors difference value distributions of the initial population, the 5 ${}^{\text{th}}$ population, 15 ${}^{\text{th}}$ population and 20 ${}^{\text{th}}$ population.

Table 1

Influence of generation steps on the optimization result

Generation steps	8	12	16	20	24
Case 1: 40 Obs.Nums without the best parameters	550.6 s	896.5 s	1321.9 s	1652.8 s	1962.6 s
Case 2: 20 Obs.Nums used the best combination of parameters	652.7 s	969.5 s	1323.4 s	1643.5 s	1953.7 s

4.5 The algorithm of DE with the choice of F

After the searching part of the experiments, the optimal vector combination ( $\textit{irad}=$ 3, $\textit{infla}=$ 1.00) in the parameter vector space is selected to confirm the data assimilation performance. As shown in Fig. 8, a comparison of the mRMSEs of the optimal vectors, suboptimal and poor selections was given. We concluded that the mRMSE of the optimal vector combination is the smallest of all vector combinations.

Figure 8.

A comparison of the mRMSEs of the optimal vectors, suboptimal and poor selections with different strength forcings $F$ .

To simulate that the atmospheric data are dense in the Northern Hemisphere (over the United States, Europe …) and very sparse in places such as the Pacific, a good strategy is taken to observe the first 20 variables in Lorenz-96 models and to see how the proposed algorithms work. As shown in Fig. 9, comparing the mRMSE of the first 20 observation used the best combination parameters with the 40 observations without the best parameters setting, we concluded that the former is better than the latter. Meanwhile, the computational complexity, as shown in Table 1 (The time unit is seconds), is increased with the increase of generation size. The running time for both of them increases gradually. However, when the generations steps exceed 16 generations, the running time of the two cases is almost the same. Therefore, combined with the LETKF, a DE algorithm-based multiple-factor method was developed to obtain better assimilation results.

Figure 9.

A comparison of the mRMSEs of the first 20 observation used the best combination of parameters with the 40 observations without the best parameter settings for $F=$ 8.

5. Discussions

Evolutionary (and related) optimization techniques can be interesting to meteorologists, but there are quite a few papers in meteorology and other applied sciences that discuss this approach. Significant effort has been made to seeking to enhance the state of DA technology by developing new coupled DA methods aimed at maximize the strengths of both forms of data assimilation while offsetting the weaknesses of each. However, Tian and Xie [19] suggested that several issues, such as how to choose the appropriate localization radius and the relaxation coefficient, still need to be addressed. In this study, in considering the optimization of DA parameters, a new scheme was developed to improve error parameterization methods [12]. By applying the DE, the best combination of the two key factors of LETKF was generated through a typical procedure using the DE operation. The impact of the new searching method on the robustness of the DA system under selected special conditions, such as different force parameters, was then investigated. We use the Lorenz-96 model to verify the proposed optimization method. The numerical experiments show that for the Lorenz-96 model with 40 variables, two parameters of the LETKF analysis can be successfully optimized with the “population” size equal to 15 and the number of “generations” equal to 20. The results indicate that the proposed DE method has the strong robustness when the system environment changes. Additionally, when compared with the traditional filtering methods, the proposed robust optimization method can improve the assimilation performance. The main differences of the technique presented in the paper are twofold. (a) Two parameters are optimized here instead of one parameter [11]. (b) The “mutation” process is somewhat different here.

6. Conclusions

There is a well-established area of parameter estimation within DA [3, 20]. With that approach, the state and parameters can be estimated simultaneously, implying, potentially, much greater efficiency than with “black-box” schemes, such as the one described in this paper. In particular, the estimation of parameters is embedded in the DA system and thus is achieved “on-line”, without the need to examine multiple systems’ configurations (which is costly). The optimization of the parameters in the DA method in an efficient manner is a priority for optimizing the performance index of the mRMSE. In this paper, the parameter optimization scheme of the DA system based on the DE algorithm is proposed due to the non-enumeration of multiparameter vector space and the discontinuity of the mRMSE index function. The effectiveness of this scheme is verified by the change in the distribution of the parameter population, the convergence of the optimal combination, and the distribution of the vector difference values. Simultaneously, the heuristic convergence of the DE method for an optimization problem is presented.

In an ensemble DA system, the state dimension of the models is typically very high, and in most cases small ensemble can be used. The localization technique is used to calculate the local analysis of each grid point to solve the problems of false correlation and lack of freedom. Additionally, different localization methods were proposed, and extensive tuning of the localization parameters is necessary to achieve the optimal results, including adaptive localization methods [21, 22]. Normally, in order to obtain most of the relevant covariances, a sufficiently large localization radius should be selected. Therefore, the relationship between ensemble numbers and localization radii should be further investigated in the new scheme in the future. Meanwhile, multi-objective dynamic DE methods can be adopted to achieve better assimilation results.

Footnotes

Acknowledgments

This work is supported by the NSFC (National Natural Science Foundation of China) project (grant number: 41861047, 41461078) and the Northwest Normal University young teachers’ scientific research capability upgrading program (NWNU-LKQN-17-6).

References

Anderson

J.L.

and Anderson

S.L.

, A Monte Carlo implementation of the non-linear filtering problem to produce ensemble assimilations and forecasts, Monthly Weather Review (1999), 2741–2758.

Houtekamer

P.L.

and Mitchell

H.L.

, Data assimilation using an ensemble Kalman filter technique, Monthly Weather Review (1998), 796–811.

Ott

et al., A local ensemble Kalman filter for atmospheric data assimilation, Tellus A: Dynamic Meteorology and Oceanography (2004), 415–428.

Hunt

B.R.

Kostelich

E.J.

and Szunyogh

, Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter, Physica D: Nonlinear Phenomena (2007), 112–126.

Hamill

T.M.

and Snyder

, Using improved background-error covariances from an ensemble Kalman filter for adaptive observations, Monthly Weather Review (2002), 1552–1572.

Kirchgessner

Nerger

and Bunse-Gerstner

, On the choice of an optimal localization radius in ensemble Kalman filter methods, Monthly Weather Review (2014), 2165–2175.

Hoteit

et al., Mitigating observation perturbation sampling errors in the stochastic EnKF, Monthly Weather Review (2015), 2918–2936.

Raanes

P.N.

Marc

, and Alberto

, Adaptive covariance inflation in the ensemble Kalman filter by Gaussian scale mixtures, Quarterly Journal of the Royal Meteorological Society (2019), 53–75.

Gharamti

, Enhanced adaptive inflation algorithm for ensemble filters, Monthly Weather Review (2017), 623–640.

10.

Attia

and Constantinescu

, An Optimal Experimental Design Framework for Adaptive Inflation and Covariance Localization for Ensemble Filters, arXiv preprint arXiv:1806.10655, 2018.

11.

Bai

Y.L.

and Li

, Evolutionary algorithm-based error parameterization methods for data assimilation, Monthly Weather Review (2011), 2668–2685.

12.

Bai

Y.L.

and Huang

, Handling error propagation in sequential data assimilation using an evolutionary strategy, Advances in Atmospheric Sciences (2013), 1096–1105.

13.

Ruiz

J.J.

Pulido

and Miyoshi

, Estimating model parameters with ensemble-based data assimilation: A review, Journal of the Meteorological Society of Japan (2013), 79–99.

14.

Wang

and Cai

, Robust data assimilation in hydrological modeling – A comparison of Kalman and H-infinity filters, Advances in Water Resources (2008), 455–472.

15.

Storn

and Price

, Differential evolution – A simple and efficient heuristic for global optimization over continuous spaces, Journal of Global Optimization (1997), 341–359.

16.

Rönkkönen

Kukkonen

and Price

K.V.

, Real-parameter optimization with differential evolution, IEEE Xploe Conference: Evolutionary Computation. [Online]. Available: doi: 10.1109/CEC. 2005.

17.

Omran

M.G.H.

Salman

and Engelbrecht

A.P.

, Self-adaptive differential evolution: in Lecture Notes in Artificial Intelligence, in: International Conference on Computational Intelligence & Security Springer-Verlag, 2005, pp. 192–199.

18.

Lorenz

E.N.

and Emanuel

K.A.

, Optimal sites for supplementary weather observations: Simulation with a small model, Journal of the Atmospheric Sciences (1998), 399–414.

19.

Tian

Xie

and Sun

, A POD-based ensemble four-dimensional variational assimilation method, Tellus A: Dynamic Meteorology and Oceanography (2011), 805–816.

20.

Frei

and Künsch

H.R.

, Sequential state and observation noise covariance estimation using combined ensemble Kalman and particle filters, Monthly Weather Review (2012), 1476–1495.

21.

Anderson

J.L.

, Localization and sampling error correction in ensemble kalman filter data assimilation, Monthly Weather Review (2012), 2359–2371.

22.

Bishop

and Hodyss

, Ensemble covariances adaptively localized with ECO-RAP. Part 1: Tests on simple error models, Tellus A: Dynamic Meteorology and Oceanography (2009), 84–96.