water

Article
Immune Evolution Particle Filter for Soil Moisture
Data Assimilation
Feng Ju * , Ru An * and Yaxing Sun
School of Earth Sciences and Engineering, Hohai University, Nanjing 211100, China; yaxingsun@njust.edu.cn
* Correspondence: jufenghhu@126.com (F.J.); anrunj@163.com (R.A.)




Received: 15 December 2018; Accepted: 14 January 2019; Published: 26 January 2019


Abstract: Data assimilation (DA) has been widely used in land surface models (LSM) to improve
model state estimates. Among various DA methods, the particle filter (PF) with Markov chain
Monte Carlo (MCMC) has become increasingly popular for estimating the states of the nonlinear and
non-Gaussian LSMs. However, the standard PF always suffers from the particle impoverishment
problem, characterized by loss of particle diversity. To solve this problem, an immune evolution
particle filter with MCMC simulation inspired by the biological immune system, entitled IEPFM,
is proposed for DA in this paper. The merit of this approach is in imitating the antibody diversity
preservation mechanism to further improve particle diversity, thus increasing the accuracy of
estimates. Furthermore, the immune memory function refers to promise particle evolution process
towards optimal estimates. Effectiveness of the proposed approach is demonstrated by the numerical
simulation experiment using a highly nonlinear atmospheric model. Finally, IEPFM is applied to a
soil moisture (SM) assimilation experiment, which assimilates in situ observations into the Variable
Infiltration Capacity (VIC) model to estimate SM in the MaQu network region of the Tibetan Plateau.
Both synthetic and real case experiments demonstrate that IEPFM mitigates particle impoverishment
and provides more accurate assimilation results compared with other popular DA algorithms.

Keywords: immune evolution algorithm; particle filter; Markov chain Monte Carlo; soil moisture;
data assimilation; Variable Infiltration Capacity

1. Introduction
Soil moisture (SM) plays a key role in the interactions between the hydrosphere, the biosphere,
and the atmosphere by governing the partitioning of mass and energy fluxes between the land and the
atmosphere [1]. As such, understanding SM is pivotal in various relevant fields, such as water resource
management, drought warning, flood and landslide modelling and prediction, irrigation management,
and even economic and policy analysis [2,3]. Generally, SM can be obtained through observation
and modeling. At the local scale, in situ observation techniques provide accurate measurements of
SM at different depths, but they are unable to characterize the SM at large spatial scales. However,
at large spatial scales, microwave remote sensing data provide a technique of estimating only the near
surface SM, limiting the microwave penetration depth. In addition, another fundamentally different
way to obtain SM information is by the application of a physical-based spatially distributed Land
Surface Model (LSM) [4]. Different methods have their respective advantages and drawbacks. The data
assimilation (DA), incorporating observed data with model simulations and keeping physical states of
the models up-to-date, is a widely-used technique in hydrology and meteorology [5–7].
Several DA algorithms have been used in hydrologic data assimilation. The Kalman filter (KF) and
its variants are widely known DA methods, but the fatal disadvantages of these methods arise from the
assumptions of the linear dynamic system with Gaussian noise. To extend KF to the nonlinear systems,
much effort has been devoted to this aspect, including extended KF (EKF). However, EKF might

Water 2019, 11, 211; doi:10.3390/w11020211 www.mdpi.com/journal/water

Water 2019, 11, 211 2 of 18

generate instabilities, which leads to divergence and inaccurate results in strong nonlinear models [8].
To deal with the drawbacks of EKF, the ensemble Kalman filter (EnKF) was proposed, which used
ensembles to approximate the probability distribution function of state variables [9]. However, the
common limitation of all the filtering techniques based on KF is filter evolution characterizing by
the Gaussian error assumption. Also, when EnKF is applied in the strongly nonlinear LSM, the filter
process is greatly simplified and restricted to the linear updating rule [6]. To overcome the mentioned
limitations, a non-parametric Monte Carlo sampling-based method in the form of particle filtering (PF)
was proposed [10].
The PF is a viable alternative DA technique. It handles the propagation of non-Gaussian
distributions through nonlinear models, unlike EnKF restrictively assumes the error distributions are
normal [11,12]. Thus, the PF has been successfully applied to estimate state variables or parameters
in nonlinear LSMs [13–15]. However, the general PF continues to suffer from some serious problems.
One of them is particle impoverishment characterizing, where most particles share a few distinct values
and the posterior distribution is insufficiently approximated, which leads to the misleading estimation
results. In order to combat the impoverishment problem, a variety of methods have been proposed.
Elsheikh et al. combined the nested sampling algorithm with the standard PF to mitigate the particle
impoverishment problem [16]. Guingla et al. proposed improving strategy based on the Markov chain
Monte Carlo (MCMC) resample move step and redefining the importance density function within the
Gaussian PF structure [17]. Bi et al. used similar ideas to assimilate brightness temperatures for SM
assimilation [18]. In this paper, a new candidate particle generating strategy for the MCMC move step
is proposed, which belongs to artificial intelligence (AI) algorithms in nature and is different from
previous studies. With the development of AI, modified PFs embedded with intelligent and heuristic
algorithms have been successfully applied in different areas, including hydrologic data assimilation,
streamflow forecasting, temperature downscaling, fault detection, and so on [15,19–21]. Recently,
Qu et al. established a machine learning model combining remote sensing data with experimental
data to reveal the SM for the construction and management of sponge cities [22]. Abbaszadeh et al.
proposed a genetic evolutionary PF algorithm combining a genetic algorithm (GA) with MCMC to
enhance hydrologic prediction [23]. Despite lots of attention being paid to AI algorithms in hydrologic
DA, little research on intelligent PFs is inspired by the biological immune system.
In this paper, a new method, named the immune evolution particle filter with MCMC (IEPFM), is
attempted to alleviate the particle impoverishment. Compared with the other methods, IEPFM owns
immune memory and diversity preservation strategies. In order to evaluate the effectiveness and
efficiency of this new DA method, two case studies are implemented. The first one uses a numerical
experiment with the Lorenz96 atmospheric model, characterizing 40 dimensions and having high
nonlinearity, and the second one is a real case of SM data assimilation with a distributed hydrological
model. Three different DA algorithms, including IEPFM, EnKF, and the differential evolution particle
filter with MCMC (DEPFM), are applied in both experiments at the same time. It is noted that
DEPFM follows the method of Vrugt et al. to enhance sample diversity, because it outperforms the
standard PF and demonstrates its effectiveness in many practical applications [15,24,25]. Comparing
the performance of IEPFM with EnKF and DEPFM, it is found that IEPFM has the best data
assimilation effect.
The main contributions of this paper are as follows: (1) Inspired by the biological immune system,
the immune evolution conception is integrated into PF with MCMC to increase particle diversity.
(2) The proposed algorithm demonstrates usefulness and effectiveness, both in the high dimensional
and nonlinearity model or in the hydrological model for SM assimilation. (3) The key parameters
of immune evolution are explored during experiments and the reference values are provided for
practical applications.
The rest of this paper is organized as follows: Section 2 outlines the theory of Bayesian filtering,
the general PF and its drawbacks, and then introduces IEPFM. Section 3 illustrates the detail process
of the numerical experiment and the SM assimilation experiment basing on the aforementioned theory.
Water 2019, 11, 211 3 of 18

The performances of the proposed approaches are evaluated and discussed in this Section. Section 4
concludes the whole paper.

2. Materials and Methods

2.1. Bayesian Filtering

In a sequential DA process, the evolution of the simulated system states can be represented
as follows:
xt = f (xt −1 , θ ) + ωt , (1)

yt = h (xt ) + vt , (2)

where the subscript t denotes the time step, f(·) is a nonlinear function expressing the system transition
(i.e., LSM) from time t − 1 to t; h(·) represents the measurement function producing forecasted
observations; xt ∈ Rn denotes a vector of the state variables at time step t, and θ ∈ Rd is a vector of the
model parameters; yt ∈ Rm is a vector of the observations, ωt is considered as process noise (i.e., model
error), and vt is the measurement noise. In most cases, ωt and vt are assumed to be independent and
as white noise, with mean zero and covariance represented by Qt and Rt , respectively.
Under the framework of Bayesian filtering, the purpose of the state estimation is to seek the
posterior probability density function (pdf) p(xt |y1:t ) of state xt based on the set of all available
observations y1:t . Such a filter consists of essentially two stages: prediction and update [26]. These two
steps are formulated as:

Z 

p xt y1:t−1 = p(xt |xt−1 ) p xt−1 y1:t−1 dxt−1 , (3)


p(y1:t |xt ) p(xt ) p(yt |xt ) p xt y1:t−1
p(xt |y1:t ) = = R 
, (4)
p(y1:t ) p(yt |xt ) p xt y1:t−1 dxt


where p xt−1 y1:t−1 denotes the posterior pdf at last time step t − 1, p(xt |xt−1 ) is the transition pdf of


the state variables. The prior pdf p xt y1:t−1 can be obtained through the prediction step (3) at time
step t; p(yt |xt ) is the likelihood function for time step t. When the new observation yt is available,
the posterior pdf p(xt |y1:t ) can be calculated through the update step (4). In the update stage, the
new observation yt is used to modify the prior density to obtain the required posterior density of the
current state.
The optimal Bayesian solution is difficult to determine since the evaluation of the integral
might be intractable [17]. Therefore, some approximate solutions, such as EnKF and PF, are used in
practical applications.

2.2. Particle Filter Method

The PF approach is an implementation of Bayesian filtering by Sequential Monte Carlo (SMC)
sampling. The main idea is to express the required posterior pdf of state variables by a set of random
particles with associated weights. As the number of particles increases, the true posterior pdf of
N
state variables is equivalently approached [26]. In detail, consider a mass of particles xit , wit i=1 are


sampled from the posterior pdf p(xt |y1:t ), where N is the particle number, xit represents the random
state variables, and wit denotes the associated weight. The posterior pdf of the state variables is
approximated as a discrete function:

N  
p(xt |y1:t ) ≈ ∑ t
w i
δ x t − x i
t , (5)
i =1

where δ(·) is the Dirac delta function.

Water 2019, 11, 211 4 of 18

For the most nonlinear systems, the analytic solution for the posterior pdf p(xt |y1:t ) is often
difficult to obtain. To deal with the difficulty in sampling indirectly from the posterior pdf,
the sequential importance sampling (SIS) method, which is the basic PF form, is often adopted.
The particles are generated from a known importance density, also called proposal distribution, which
is easily sampled and represented by q(xt |y1:t ). The importance weight for each particle is updated in
a recursive form, as follows:
p yt xit p xit xit−1



i i
w t = w t −1
, (6)
q xi  xi , y

t t −1 t

Equation (6) provides the mechanism to sequentially update the importance weights, given
an appropriate choice of the proposal distribution q xit xit−1 , yt . Consequently, the choice of the


proposal distribution is one of the most critical issues in the design of PF. The optimal proposal density
makes the variance of importance weights minimal and significantly improves the efficiency of PF.
An appropriate choice for the proposal distribution is usually expressed for simplicity as follows:
     
q xit xit−1 , yt = p xit xit−1 , (7)
 

When the transitional prior pdf is chosen as the proposal distribution, the importance weights
only depend on their past values and the likelihood p yt xit . Then, the weights updating formula in


Equation (6) is reduced to   

wit = wit−1 p yt xit , (8)


The normalized weight is calculated by

wit
eit =
w
2 , (9)
∑iN=1 wit

2.3. Problem Formulation

The degeneracy phenomenon is one of the major drawbacks of the basic PF scheme. In fact,
after a few iterations, all but one particle have negligible weights. These particles’ contribution to the
approximation of the posterior pdf is close to zero. This implies that a large computational resource
is wasted to update these particles. The degeneracy of PF is measured by the effective sample size
Ne f f [26], which is defined by
1
Ne f f =
2 , (10)
N
∑ i =1 weit

where weit is the normalized weight. According to the Ne f f definition, the smaller Ne f f indicates the
more severe degeneracy. To solve the particle degeneracy problem, a resampling procedure is usually
adopted in the PF. This process eliminates the particles with small important weights and replaces
them by multiplying particles with large weights. The process is illustrated in Figure 1. In the past
decades, various resampling methods have been proposed in literature, and the most commonly used
are multinomial resampling, stratified resampling, and residual resampling [13]. After the resampling
process, all particles are the same weight. Thus, the posterior pdf is approximated by

1 N  
p(xt |y1:t ) ≈ ∑
N i =1
δ xt − xit , (11)
procedure may lead to a dramatic loss of diversity of the offspring particles, characterized by the new
particles sharing a few distinct values gathering at a few select places. In other words, after several
recursive calculations, the new particle set is difficult to reflect the true posterior pdf. Therefore, how
to make appropriate changes to the multiplicative particles according to certain rules to obtain
diversity
Water and
2019, 11, 211reduce impoverishment is the core issue. This is the original idea of the modified 5 of 18
particle filter algorithm based on the biological immune system.

Resampling

Impoverishment

Figure 1. The particle impoverishment phenomenon caused by resampling procedure.

However,
2.4. Immune the other
Evolution serious
Particle Filterproblem
with MCMCis particle impoverishment sparked by the resampling
method. In Figure 1, the solid circles represent the parent particles and the hollow circles denote
2.4.1.
the Immuneparticles.
offspring EvolutionThe
Algorithm (IEA)
size of the circles represents the corresponding weights of the particles.
The blue solid circles are selected particles and the grey solid circles are eliminated particles. The above
The immune system is a multifunction biological defense system, which has evolved to protect
curve is the prior pdf, and the curve shown below presents the true posterior pdf. As described
bodies against infections from viruses, bacteria, fungi, and also worms and tumor cells [27]. Once
in Figure 1, in the resampling procedure, only the particles with large weights are selected as the
pathogens invade the body, they will become antigens and provoke the immune response. When the
parent particles, others are abandoned. The offspring particles are replicas of the selected particles.
immune system detects the foreign invasive antigens, the immune cells will propagate a mass of
This procedure may lead to a dramatic loss of diversity of the offspring particles, characterized by
cloned antibodies. The cloned antibodies are selected from the top candidate antibodies listed in
the new particles sharing a few distinct values gathering at a few select places. In other words, after
descending order by selection probability. Through the division and differentiation of cells, immune
several recursive calculations, the new particle set is difficult to reflect the true posterior pdf. Therefore,
cells generate diverse antibodies to destroy the antigens. During the immune study process, the cells
how to make appropriate changes to the multiplicative particles according to certain rules to obtain
propagating the antibodies will be saved as memory cells. When the same or similar antigens appear
diversity and reduce impoverishment is the core issue. This is the original idea of the modified particle
again, the memory cells will quickly react and produce better antibodies [28]. The immune memory
filter algorithm based on the biological immune system.
is the distinct feature of the immune system.
The artificial
2.4. Immune immune
Evolution evolution
Particle algorithm
Filter with MCMC is a type of algorithm inspired by biological immune
system theories. It simulates the adaptive capacity of the immune system to propagate antibodies
2.4.1. Immune Evolution Algorithm (IEA)
The immune system is a multifunction biological defense system, which has evolved to
protect bodies against infections from viruses, bacteria, fungi, and also worms and tumor cells [27].
Once pathogens invade the body, they will become antigens and provoke the immune response.
When the immune system detects the foreign invasive antigens, the immune cells will propagate a
mass of cloned antibodies. The cloned antibodies are selected from the top candidate antibodies listed
in descending order by selection probability. Through the division and differentiation of cells, immune
cells generate diverse antibodies to destroy the antigens. During the immune study process, the cells
propagating the antibodies will be saved as memory cells. When the same or similar antigens appear
again, the memory cells will quickly react and produce better antibodies [28]. The immune memory is
the distinct feature of the immune system.
The artificial immune evolution algorithm is a type of algorithm inspired by biological immune
system theories. It simulates the adaptive capacity of the immune system to propagate antibodies
against the foreign invasive antigens. It implements the recognition of invasive antigens, and the
generation of diverse antibodies, immune memory, self-regulation, and other functions [29]. There have
been many successful applications of IEA, such as optimization, data mining, video target tracking,
industrial applications, and Internet of Things services [29–32].

2.4.2. Immune Evolution Particle Filter Proposed

The aim of this study is solving the particle impoverishment problem. Based on the same demand
of improving individual diversity and quality, PF and immune evolution can be linked together.
Corresponding to the biological immune system, the foreign invasion antigens indicate the state
Water 2019, 11, 211 6 of 18

or parameters of optimal problems to be solved; the antibodies generated by the immune system
represent the candidate particles for the variable estimation. The distinct features, which are useful for
improving diversity, are discussed below.
Similar to the basic GA, the immune system produces offspring antibodies by genetic operators
such as selection, crossover, and mutation. However, the selection process is not only according to
the fitness of individuals, but also considering the diversity of the selected ones. According to the
crossover rate, parent antibodies exchange part of their information to generate the new offspring
antibodies. The crossover operator expands the search space of candidate solutions and finally achieves
the global search purpose. At the same time, the mutation may also happen in the offspring antibodies,
according to the mutation rate. This process produces new genes which are not included in the initial
population, or recovered genes lost during the selection. The mutation operator maintains population
diversity and prevents premature convergence. This antibody diversity preservation mechanism is
worth learning to solve the impoverishment problem of standard PF. Until reaching the termination
condition, the iteration is finished.
In practice, as the first step of this mechanism implementation, the particle selection is
implemented to achieve the goal of promotion and suppression of particles based on the fitness
and the concentration. The term fitness measures the degree of matching between observations and
candidate particles by fitness function. The term concentration indicates the similarity degree of
the particle with other particles by concentration function. The particles with higher fitness will be
promoted and the ones with higher concentration will be suppressed. In this way, the diversity of
particles is further ensured. Thus, the problem of super individuals is mitigated and the solution space
is extended. In addition, inspired by the immune system, the immune memory is designed to store the
optimal particles and generate better offspring particles for each iteration.
The abovementioned diversity preservation idea, immune memory, and genetic operators are
adopted in the proposed algorithm IEPFM. These functions need to be improved to fit the particle filter
in DA. Figure 2 displays the algorithm process of IEPFM. The idea of immune evolution is described
in the red dotted frame. Note that the aim of immune evolution process is creating new candidate
particles for MCMC simulation, which happens after residual resampling.
The detailed assimilation process includes the following steps:
N
Step (1). At the assimilation initial time step t = 0 initializes particle set x0i , w0i i=1 from the state


variable prior pdf p(x0 ), where w0i = 1/N.

N
Step (2). The particles set state propagates according to the model and particle set xit−1 , wit−1 i=1 .

N
The new particles xit , wit i=1 are forecasted through model, where wit = 1/N.


Step (3). If the new observation is available at current time step t, use the new observation to
calculate the selection probability of particles. Otherwise, skip to step (13) and the final estimated state
is obtained.
Step (4). At each observation time, the memory cell is constructed. The length of memory cell is
M, half of the total number of particles (i.e., M = N/2).
Step (5). In this study, the fitness is proportional to the likelihood. However, a common choice of
the likelihood density function is the Gaussian distribution that describes the misfit between the state
predictions xit and the observations yt , scaled by the (usually a priori defined) standard deviation of
the observation error σ [17]. Thus, the fitness can be estimated as follows:
m
T
!
yt − xit yt − xit

1
f itnessit = √ exp − , (12)
2πσ 2σ2

N
Step (6). Calculate the normalized weights of particles xit , w
eit i=1 according to (8) and (9). Then,


the effective sample size Ne f f can be calculated according to (10). If Ne f f ≤ Nr , go to the next step,
N
otherwise, forecast particles xit , wit i=1 through the model and skip to step (13).

 fitnessti = 
1  t t t t 
 exp  − , (12)
 2πσ   2σ 2 
 

{x }
N
Step (6). Calculate the normalized weights of particles i
t
, w ti according to (8) and (9). Then,
i =1

the effective sample size N eff can be calculated according to (10). If N eff ≤ Nr , go to the next step,
Water 2019, 11, 211 7 of 18
{x }
N
i i
otherwise, forecast particles t
,w t
through the model and skip to step (13).
i =1

Start

Initialize particles set x t −1

N
New observation?

Y
Calculate the weights of particles

N
Particle degeneracy?

Y
Residual resampling Forecast particles x t through Model

Calculate fitness Memory

Forecast particles x t through Model
and concentration Replace cell

Compute selection
probability Update

Crossover and p
Forecast particles x t through Model
mutation operation

N
Evolution terminated? Calculate the Metropolis acceptance probability α

Y
N
Proposal particles set x t −1
p U[0,1] ≤ α
Y

Keep the particles xt Accept the proposal particles x tp

Optimal estimation state

t=t+1 N
Is it time terminated?

Y
end

Figure 2.2.The
Figure Theoverview of IEPFM.
overview of IEPFM.

Step (7). Calculate the concentration of particles as follow:

nsimi
Cti = , (13)
N
where nsimi is the number of particles with high affinity of particle i. Particles with high concentration
should reduce their chosen probability, and vice versa.
Step (8). The selection probability of each particle is obtained as:

f itnessit Ci
pchoose = ω · pfit + (1 − ω ) · pcon = ω · + (1 − ω ) · , (14)
N j β
∑ f itnesst
j =1

where ω and β are adjustment constants. In IEPFM, the selection probability of each particle is
determined by its fitness and concentration. Then, the particle set will be sorted by selection probability.
The optimal particle, which has the highest selection probability, will be stored in the memory cell at
every iteration. If the memory cell is full, the worst particle of them will be replaced.
Water 2019, 11, 211 8 of 18

Step (9). Crossover can operate between two cross parents, which are selected randomly according
to the crossover probability pc . Considering the state of double-coding and the simplicity, parent
particles produce offspring by linear arithmetic crossover. The crossover operator is defined as:
(
x0 1 = kx1 + (1 − k)x2
, (15)
x0 2 = (1 − k )x1 + kx2

where k is a random number generated from the interval [0, 1].

Step (10). Mutation strategy is simulated to promote the diversity of particles. Mutation can
operate between two mutation parents, which are selected randomly according to the mutation
probability pm . The mutation operator is expressed as:
0
x = 2xH − xL + ε, (16)

where xH and xL represent the particles with high and low selection probability, respectively; ε is a
perturbation vector, which is drawn from the uniform distribution on U(−b, b) with b, which is small
compared to the width of the target distribution.
Step (11). Cycle from step (5) to step (10) until the maximum
n evolution generation is achieved.
N  N
Otherwise, use MCMC simulation to obtain optimal particle set x̂it , ŵit i=1 , where ŵit i=1 = 1/N.
Further details are described below.
Step (12). Calculate the optimal estimate state as

N
xt =
e ∑ x̂ti × ŵit , (17)
i =1

Step (13). If time step is terminated, the whole process stops. Otherwise, make t = t + 1 and
return to step (2), and cycle the above procedures.

2.4.3. Markov chain Monte Carlo (MCMC) Simulation

The MCMC and sequential Monte Carlo methods have emerged as the two main tools to sample
high dimensional probability distributions [33]. The key idea of MCMC simulation is based on using
Markov chain to explore the target distribution. The immune evaluation particle filter procedure
p
generates candidate particles xt−1 . Then, each candidate particle is re-evaluated from time t − 1 to t by
Equation (1), and the Metropolis acceptance probability α is used to determine whether to move from
p
xt−1 to xt−1 or not; α is calculated as follows

  p 
p(yt xt )
α = min1, , (18)
p (yt |xt )

Only when µ ≤ α, where µ ∼ U[0, 1], the candidate particle is accepted and replaces the old
particle. After the IEPFM process, the new particle set will have a distribution closer to the posterior
pdf, and will also have more diversity, thus reducing the potential of particle impoverishment.

3. Results and Discussion

In this section, firstly, the Lorenz96 model is used as the numerical simulation experiment to
compare the performance of IEPFM with EnKF and DEPFM. Then, IEPFM is applied for soil moisture
data assimilation. The performance of each method is evaluated against the true values using the
average root mean square error (RMSE) over all state variables.
Water 2019, 11, 211 9 of 18

3.1. Lorenz96 Model Experiment

The Loren96 model [34] is a high dimensional and strongly nonlinear model often used to validate
the new DA algorithms. This model simulates the time evolution of an atmospheric quantity and is
described by the ordinary differential equations:

dx j

= x j+1 − x j−2 x j−1 − x j + F, j = 1, 2 · · · , J, (19)
dt
The cyclic boundary conditions are defined by x−1 = x J −1 , x0 = x J and x1 = x J +1 ; F is a constant
external forcing term, which is usually chosen as F = 8 for high chaotic behavior.
In this paper, the 40 dimensions Lorenz96 model is chosen for the experiment. The fourth-order
Runge-Kutta scheme with a step of ∆t = 0.05, corresponding to six hours in real life, is used to
integrate the model. The model is initialized by choosing F = 8, except for F = 8.01 at variable x20 ,
and running for 1000 time steps. The end states of that run are considered as the initial condition for
the DA experiment. The true state variables are obtained by integrating the model. All of the model
state variables are observed every five time steps and added to the Gaussian observational noise with
2 = 0.2. The initial particles or ensembles are generated by adding the Gaussian
error covariance σobs
2 = 0.1. The model error covariance is σ2
noise with initial covariance σinit model = 0.1.

3.1.1. Sensitivity Analysis

To further understand the influence of key parameter settings of IEPFM, a set of experiments are
conducted for sensitivity analysis. The average RMSEs are used to evaluate the performance of each
setting average for 30 independent runs. The crossover probability pc is varied from 0.1 to 0.9 by steps
of 0.1. The experimental results are summarized in Figure 3. For the same mutation probability, the
value of RMSE is near 1.2 when pc ≥ 0.5. The lowest RMSE is obtained when pc = 0.8. Thus, pc = 0.8 is
Water 2019,for
selected 11 FOR PEER REVIEW
IEPFM. 10

Figure 3.
Figure Boxplots
3. Box plotsof
ofthe
theRMSE
RMSE with
with different
different crossover
crossover probability.
probability.

The same as the crossover probability, the mutation probability ranges from 0.1 to 0.9 steps by 0.1,
The same as the crossover probability, the mutation probability ranges from 0.1 to 0.9 steps by
and the corresponding RMSEs are presented by box plots in Figure 4. It is clearly observed that the
0.1, and the corresponding RMSEs are presented by box plots in Figure 4. It is clearly observed that
results are almost unanimously obtained when pm ≥ 0.3. Therefore, the parameter setting of crossover
the results are almost unanimously obtained when pm ≥ 0.3 . Therefore, the parameter setting of
and mutation probability has no significant impact on IEPFM as a result of MCMC simulation. This
crossover
conclusionand mutation probability hasby no significantet impact onthis
IEPFM
study,aspca=result of pMCMC
is consistent with that drawn Abbaszadeh al. [23]. In 0.8 and m = 0.3
simulation. This conclusion is consistent with that drawn by
are adopted, since RMSE resulting from this setting is slightly better. Abbaszadeh et al. [23]. In this study,
As and
p c =0.8 =0.3 ofare
pm size
for the adopted,
the particlesince RMSE 5resulting
set, Figure gives barfrom
plotsthis
ofsetting is slightly
the RMSEs better.
and computing time
for different ensemble sizes (i.e., 25, 50, 100, 150, and 200) for 30 independent runs. As shown in
Figure 5, the RMSE gradually decreases and computing time sharply increases, along with ensemble
 Figure 3. Box plots of the RMSE with different crossover probability.

The same as the crossover probability, the mutation probability ranges from 0.1 to 0.9 steps by
0.1, and the corresponding RMSEs are presented by box plots in Figure 4. It is clearly observed that
Water 2019, 11, 211 10 of 18
the results are almost unanimously obtained when pm ≥ 0.3 . Therefore, the parameter setting of
crossover and mutation probability has no significant impact on IEPFM as a result of MCMC
size becomingThis
simulation. larger. It can is
conclusion beconsistent
found that thethat
with performance of IEPFM improves
drawn by Abbaszadeh et al. [23].as
Inthe
thisensemble
study,
sizep increases.
=0.8 and But
p the
=0.3 RMSE
are shows
adopted, almost
since no
RMSE obvious
resulting change
from since
this N
setting >
is 100. After
slightly comprehensive
better.
c m
consideration of RMSE and computing time, ensemble size is set to 100 in this study.

Water 2019, 11 FOR PEER REVIEW 11

Figure 4. Box plots of the RMSE with different mutation probability.
Figure 4. Box plots of the RMSE with different mutation probability.

As for the size of the particle set, Figure 5 gives bar plots of the RMSEs and computing time for
different ensemble sizes (i.e., 25, 50, 100, 150, and 200) for 30 independent runs. As shown in Figure
5, the RMSE gradually decreases and computing time sharply increases, along with ensemble size
becoming larger. It can be found that the performance of IEPFM improves as the ensemble size
increases. But the RMSE shows almost no obvious change since N>100. After comprehensive
consideration of RMSE and computing time, ensemble size is set to 100 in this study.

Figure 5.5.Bar
Figure Barplots
plotsofofthe
theRMSE
RMSE and
and computing timewith
computing time withthe
thesame
sameconditions.
conditions.

3.1.2. Performance
3.1.2. PerformanceComparisons
Comparisons
Three
Three DA algorithmswith
DA algorithms withdifferent
different particles
particles or ensemblemembers
or ensemble membersare aresimulated
simulated forfor 200
200 time
time
steps
stepsand thethe
and state
statevariables
variablesareareestimated
estimated by by assimilating
assimilating the theobservations
observationsinto intothethe Lorenz96
Lorenz96 model.
model.
AllAll
40 40
state
state variables are selected for comparison. The average RMSEs of 30 independent runs areare
variables are selected for comparison. The average RMSEs of 30 independent runs
used to evaluate
used to evaluate thetheperformance
performanceof ofDA
DAalgorithms.
algorithms. In In addition,
addition,for forIEPFM,
IEPFM,the the crossover
crossover probability
probability
pc and the the
p c and mutation
mutation probabilitypmp are
probability m
setset
are toto0.80.8and
and0.3,
0.3,respectively.
respectively. TheTheparameters
parametersofof DEPFM
DEPFM
areare
referred to Vrugt et al. [25]. Table 1 shows the experiment results of EnKF,
referred to Vrugt et al. [25]. Table 1 shows the experiment results of EnKF, DEPFM, and IEPFM DEPFM, and IEPFM
with
withdifferent
differentensemble
ensemblesize. size. From
From Table
Table 1, 1, it
it can
can bebe observed
observedthat thatEnKF
EnKFperforms
performspoorly
poorlyin in
thethe
highly
highlynonlinear
nonlineardynamics
dynamics model.
model. However,
However, thetheparticle filters
particle with
filters MCMC
with MCMC simulation improve
simulation improvefilter
performance significantly, and the REMSs diminish with the increment of
filter performance significantly, and the REMSs diminish with the increment of the ensemble size.the ensemble size. The best
state
Theestimation
best stateresults belong
estimation to IEPFM,
results belongwhich offerswhich
to IEPFM, the lowest
offersRMSE at all different
the lowest RMSE atensemble sizes,
all different
as ensemble
shown in sizes, as shown in the table.
the table.
To intuitively present the effectiveness of IEPFM, the average RMSE (N = 100) of every time step
is displayed in Figure 6. AsTable shown 1. Comparison
in Figure of 6, REMS for the
the state Lorenz96
values model.
are not correctly estimated by EnKF,
except for the beginning of
Ensemble Sizethe assimilation.
EnKF This is most likely
DEPFM because normal
IEPFM distribution is assumed
to derive from EnKF analysisN = 25 scheme, limiting 4.6244 its performance 3.2236for nonlinear2.4721 and non-Gaussian models.
IEPFM gives much more N = accurate
50 estimation
4.6222 results than1.8439 EnKF and DEPFM. 1.8256
N = 100 4.6128 1.5260 1.3661
N = 150 4.6710 1.2242 1.1160
N = 200 4.6685 1.1864 1.1038
p c and the mutation probability p m are set to 0.8 and 0.3, respectively. The parameters of DEPFM
are referred to Vrugt et al. [25]. Table 1 shows the experiment results of EnKF, DEPFM, and IEPFM
with different ensemble size. From Table 1, it can be observed that EnKF performs poorly in the
highly nonlinear dynamics model. However, the particle filters with MCMC simulation improve
filter2019,
Water performance
11, 211 significantly, and the REMSs diminish with the increment of the ensemble 11 size.
of 18
The best state estimation results belong to IEPFM, which offers the lowest RMSE at all different
ensemble sizes, as shown in the table.
Table 1. Comparison of REMS for the Lorenz96 model.
Table 1. Comparison of REMS for the Lorenz96 model.
Ensemble Size EnKF DEPFM IEPFM
Ensemble Size
N = 25 EnKF
4.6244 DEPFM
3.2236 IEPFM
2.4721
NN= =2550 4.6222
4.6244 1.8439
3.2236 1.8256
2.4721
NN==50100 4.6128
4.6222 1.5260
1.8439 1.3661
1.8256
N = 150 4.6710 1.2242 1.1160
NN= =
100
200 4.6128
4.6685 1.5260
1.1864 1.3661
1.1038
N = 150 4.6710 1.2242 1.1160
N = 200 4.6685 1.1864 1.1038

Water 2019, 11 FOR PEER REVIEW 12

assumed to derive from EnKF analysis scheme, limiting its performance for nonlinear and non-
6. Comparison
Figure IEPFM
Gaussian models. of the assimilation
gives much results of average RMSEthan
at every time step.
Figure 6. Comparison of the more accurate
assimilation estimation
results results
of average EnKF
RMSE at every and
time DEPFM.
step.
3.1.3. Comparison of Computation Time
3.1.3.To
Comparison
intuitively of Computation
present Time
the effectiveness of IEPFM, the average RMSE (N=100) of every time step
To compare
is displayed in the
Figurecomplexity
6. As shown ofinthe three6,different
Figure the state algorithms,
values thecorrectly
are not average computation time
To compare the complexity of the three different algorithms, the averageestimated by EnKF,
computation time
corresponding
except for the toto Table 1
beginning is provided in Figure 7. It is clearly seen that the computation time increases
corresponding Table 1 isof the assimilation.
provided in Figure 7.This is mostseen
It is clearly likely
thatbecause normal distribution
the computation is
time increases
with the increments of ensemble size. The obvious difference is that EnKF increases asymptotically
with the increments of ensemble size. The obvious difference is that EnKF increases asymptotically
and approximately linearly. However, the computation time of other methods increases rapidly when
and approximately linearly. However, the computation time of other methods increases rapidly
the ensemble size is greater than 100. This finding is of no surprise, since the IEPFM and DEPFM
when the ensemble size is greater than 100. This finding is of no surprise, since the IEPFM and
are combining PF with MCMC simulation, which improve particles diversity through differential
DEPFM are combining PF with MCMC simulation, which improve particles diversity through
evolution or immune evolution during the filtering process. From Figure 7, EnKF is the most efficient,
differential evolution or immune evolution during the filtering process. From figure 7, EnKF is the
computationally, and IEPFM the least. It is not hard to see from Table 1 that EnKF exhibits rather
most efficient, computationally, and IEPFM the least. It is not hard to see from Table 1 that EnKF
poor precision. In addition, the gap between IEPFM and DEPFM grows from ensemble size N = 150,
exhibits rather poor precision. In addition, the gap between IEPFM and DEPFM grows from ensemble
since IEPFM searches solution space based on crossover and mutation, and spends time on sorting
size N=150, since IEPFM searches solution space based on crossover and mutation, and spends time
particles. Considering the RMSE index and computation time, IEPFM is suitable for assimilation of
on sorting particles. Considering the RMSE index and computation time, IEPFM is suitable for
high dimensional and nonlinearity problems for much smaller ensemble size.
assimilation of high dimensional and nonlinearity problems for much smaller ensemble size.

Figure
Figure 7. Comparison of
7. Comparison of the
the average
average computation
computation time
time at
at different
different ensemble
ensemble sizes.
sizes.

3.2. Soil Moisture Data Assimilation Experiment

3.2.1. Study Area

Water 2019, 11, 211 12 of 18

3.2. Soil Moisture Data Assimilation Experiment

3.2.1. Study Area

To understand land-atmosphere interactions of the Tibetan Plateau and their influence on the
climate of South-East Asia, the SM monitoring of this region is most important. The Tibetan Plateau
observatory of plateau scale soil moisture and temperature consists of three regional scale in situ
reference networks, including the Maqu network, the Naqu network, and the Ngari network. The
Maqu network was installed on the north-eastern fringe of the Tibetan Plateau (33◦ 300 –34◦ 150 N,
101◦ 380 –102◦ 450 E) in July 2008. The network is located in the Yellow River Source Region and in the
south of Maqu County in the Gansu province, China (Figure 8). It is also helpful for water management
sustainability research, such as reducing river sedimentation, soil water retention, and flood mitigation
of the Yellow River [35]. It covers the large valley of the Yellow River and the hills around it, which are
characterized by a uniform land cover of short grassland grazed by sheep and yaks. In this area, the
elevations of stations range from 3430 m to 3750 m above mean sea level. According to the Koeppen
Classification System, this site has a cold and wet climate, with rainy summers and dry winters due to
monsoons [36].
Water 2019, 11 FOR PEER REVIEW 13

Figure 8. (a)
Figure8. (a) Location
Location of
of the
the study
study area
area (black
(black rectangle);
rectangle); (b)
(b) SRTM
SRTM DEM
DEM of
of Maqu
Maqu monitoring
monitoring network
network
and location of the monitoring sites of the network.
and location of the monitoring sites of the network.
3.2.2. The Land Surface Model and Data Preparation
3.2.2. The Land Surface Model and Data Preparation
The LSM used in this study is the Variable Infiltration Capacity (VIC) model. The VIC model
The LSM used in this study is the Variable Infiltration Capacity (VIC) model. The VIC model is
is a semi-distributed LSM accounting for the balances of water and energy, originally developed
a semi-distributed LSM accounting for the balances of water and energy, originally developed jointly
jointly at the University of Washington and Princeton University. It is a grid cell based macroscale
at the University of Washington and Princeton University. It is a grid cell based macroscale hydrology
model used to simulate various hydrologic variables, such as SM, evapotranspiration, snow
accumulation and melt, runoff, baseflow, and various heat fluxes [37]. The surface runoff and
baseflow are subsequently routed through the river network based on grids, and simulate streamflow
at selected points within the basin. The model can run either at a sub-daily time step with a full energy
Water 2019, 11, 211 13 of 18

hydrology model used to simulate various hydrologic variables, such as SM, evapotranspiration, snow
accumulation and melt, runoff, baseflow, and various heat fluxes [37]. The surface runoff and baseflow
are subsequently routed through the river network based on grids, and simulate streamflow at selected
points within the basin. The model can run either at a sub-daily time step with a full energy balance,
or at daily time step in water balance mode.
This study performs the VIC model version 4.2.c. In order to drive the model, VIC requires
meteorological forcing files, vegetation data, and soil information for each grid cell. The meteorological
forcing files, including at least precipitation, maximum and minimum temperature, and the average
wind speed, are extracted from the China ground climate daily dataset (version 3.0) developed by
the China Meteorological Information Center. The vegetation information is obtained from Maryland
land cover dataset and the soil information is gathered from the China Soil Map Based Harmonized
World Soil Database (version 1.1) provided by Cold and Arid Regions Sciences Data Center at Lanzhou
(http://westdc.westgis.ac.cn).
To obtain the VIC model parameters of the study area, simulations were performed over the
Source Region of the upper Yellow River in China at a spatial resolution of 0.25◦ , which made 250
model grid cells. The VIC model was calibrated by forcing the model and adjusting parameters
that govern infiltration and baseflow recession to match simulated streamflow, with naturalized
observations obtained from the selected hydrologic station in the same recording period of 2007–2013.
The assumption made is that the adjusted parameters are the optimal values, and keep constant in the
following experiment.

3.2.3. Observations Preparation

The Maqu network spans an area of approximately 40 km × 80 km, covering more than one
resolution grid cell. In the study area, the observed soil moisture of participating assimilation grid
cells is calculated from monitoring sites by inverse distance weighted interpolation method. Since
the depth of top soil layer of VIC model is 10 cm and all record data of the Maqu network is at
least at 5 and 10 cm, the only soil moisture at 10 cm is updated by the assimilation of observed soil
moisture. When assimilating the weighted mean of in situ SM data, the biases between the single point
measurement-based and model-based SM cannot be avoided. The bias correction is necessary and
required before assimilation experiment. In this paper, the cumulative distribution function (CDF)
matching technique is used to rescale the in situ soil moisture into the VIC model space. Finally, the
VIC model begins to run from March 1, 2011, to eliminate the spin-up period, but only the data after
March 31, 2012, is it used for the data assimilation experiments at daily time step.

3.2.4. Results and Analysis

The SM assimilation experiments are performed to assess whether IEPFM can improve the
VIC model simulations better than other DA methods in the forecast period. The focus is on the
effectiveness of IEPFM in solving particle impoverishment. The classical assimilation algorithms EnKF
and DEPFM are used to assimilate observed data into the VIC model to estimate surface SM on a daily
time step scale. Then, the forecast results of three DA methods are compared with the observations.
The ensemble size is set to 100 and the initial particles are sampled from a Gaussian distribution with
mean zero and standard deviation σmodel = 0.05, obtained from the model error. According to the
specific calibration result of the Maqu network, the RMSE is reduced from 0.06 to 0.02 m3 /m3 , which
can be considered as the absolute accuracy of each station of this network [36]. One option is to use this
accuracy value to characterize the observation error, which is same across the three methods. The other
parameters of evolution algorithms are followed to the numerical simulation experiment.
1. Soil moisture analysis
The DA improvements of EnKF, DEPFM, and IEPFM are assessed by comparing the simulations
to the in situ soil moisture observations every two days, shown in Figure 9. The daily precipitation
3.2.4. Results and Analysis
The SM assimilation experiments are performed to assess whether IEPFM can improve the VIC
model simulations better than other DA methods in the forecast period. The focus is on the
Water 2019, 11, 211
effectiveness of IEPFM in solving particle impoverishment. The classical assimilation algorithms 14 of 18

EnKF and DEPFM are used to assimilate observed data into the VIC model to estimate surface SM
on a is
data daily
also time step scale.
displayed during Then, the forecast
the simulation resultsrepresented
period, of three DA by methods
stem plotsare
in compared
the figure. with
Figurethe9
observations. The ensemble size is set to 100 and the initial particles are sampled from
indicates the following points: Firstly, the results of the three assimilation methods can reflect the a Gaussian
distribution
overall change with mean
trend zero
of the and standard
surface SM. It is not hard toσsee
deviation m odthat
el
=0.05 , obtained
strong fromexists
consistency the model
betweenerror.
the
According
change to of
trend theSM
specific calibration
and the resultdata.
precipitation of the Maqu network,
Especially with thethe RMSE of
increase is precipitation,
reduced from the0.06SM
to
0.02 m3/m3increases
obviously , which can be flood
in the considered
seasonas in the absolute
June. accuracy
In that period, ofgap
the each station the
between of this network
model [36].
simulation
One option
result is toand
with DA use observations
this accuracy has
value to characterize
increased. the observation
According error,
to the forecast SM which is same
curves of theacross
three
the three it
methods, methods.
could beThe other
clearly parameters
seen that IEPFM of generates
evolutionthealgorithms are followed
most accurate to theby
SM, followed numerical
DEPFM,
simulation
and experiment.
EnKF produces the worst results farthest from the observations.

Figure 9. Comparison of the assimilation results of EnKF, DEPFM, and IEPFM with VIC model.
Figure 9. Comparison of the assimilation results of EnKF, DEPFM, and IEPFM with VIC model.

2.
1. Performance
Soil moisture metrics
analysis
The DA improvements
The evaluation of the three of EnKF, DEPFM, and
DA methods’ IEPFM areisassessed
performance based onby comparing
the the simulations
three metrics: the RMSE,
to the in situ soil moisture observations every two days, shown in Figure 9. The
the mean absolute bias (MAB), and the correlation coefficient (R). The lower RMSE or MAB, and daily precipitation
datagreater
the is also R,
displayed
indicateduring
better the simulationofperiod,
performance represented
forecast SM. Table by stem plots
2 presents theinRMSE
the figure.
(m3 /mFigure 9
3 ), the
indicates
MAB (m3 the
/m3following
), and the points:
R (−) ofFirstly, the DA
the three results of the three
methods. assimilation
The conclusion canmethods
be drawn canthat
reflect
IEPFM the
overall change trend of the surface SM. It is not hard to see that strong consistency exists
algorithm leads to the best statistic scores for the assimilation experiment. IEPFM shows a strong between the
correspondence with the observation, as indicated by the highest correlation (0.9502) and the lowest
RMSE (0.0295 m3 /m3 ) and MBA (0.0225 m3 /m3 ). As revealed by Table 2, the forecast SM obtained
from IEPFM is closer to the observed data than DEPFM and EnKF.

Table 2. Comparison of the performance of three DA algorithms.

DA Algorithms RMSE MAB R

EnKF 0.0388 0.0303 0.9103
DEPFM 0.0372 0.0291 0.9181
IEPFM 0.0295 0.0225 0.9502

3. Assimilation frequency analysis

In order to further understand the effects of different assimilation frequencies on prediction, four
experiments were designed. When other assimilation parameters were consistent, only the assimilation
cycles of four group experiments changed following every two days, four days, eight days, and sixteen
days. It meant that only one observation was available and could be assimilated into model per cycle.
The results of experiments are illustrated in Figure 10. Overall, similar trends can be observed. With
the extension of assimilation cycle, the prediction error increases. As shown in the figure, IEPFM
performs best with the lowest RMSE among the three algorithms, no matter how the assimilation cycle
changes, and DEPFM produces more accurate predictions than EnKF.
 assimilation cycles of four group experiments changed following every two days, four days, eight
days, and sixteen days. It meant that only one observation was available and could be assimilated
into model per cycle. The results of experiments are illustrated in Figure 10. Overall, similar trends
can be observed. With the extension of assimilation cycle, the prediction error increases. As shown in
the2019,
Water figure, IEPFM performs best with the lowest RMSE among the three algorithms, no matter how
11, 211 15 of 18
the assimilation cycle changes, and DEPFM produces more accurate predictions than EnKF.

Figure Bar
10.10.
Figure plots
Bar plotsofofthe
theRMSEs
RMSEsobtained
obtained by
by three methods with
three methods withthe
thesame
sameparameters
parameters and
and different
different
assimilation frequencies.
assimilation frequencies.

4. Particles diversity analysis

4. Particles diversity analysis
The particle diversity based on IEPFM is illustrated in Figure 11. Two subfigures describe the
The particle
particles diversity
and their basednormalized
associated on IEPFM weights
is illustrated in Figure experiment
of assimilation 11. Two subfigures describe
using DEPFM andthe
particles
IEPFMand when their
theassociated normalized
assimilation time stepsweights ofand
are t = 45 assimilation experiment
t = 70. Although using DEPFM
it is difficult andjudge
to directly IEPFM
Water 2019, 11 FOR PEER REVIEW 16
when the assimilation time steps are t = 45 and t = 70. Although it is difficult to directly judge the
accuracy of the
the accuracy oftwo
the algorithms from the
two algorithms figure,
from it is clear
the figure, that
it is particles
clear generated
that particles by the two
generated bymethods
the two
are able to
methods represent
are the true posterior
able to represent the true distributions. In each figure,
posterior distributions. In each the greenthe
figure, andgreen
blue andparticles
blue
concentrate the surrounding
particles concentrate regions of
the surrounding observation,
regions which also
of observation, whichindicates the convergence
also indicates to the
the convergence
estimates. Overall,
to the estimates. the bluethe
Overall, particle’s aggregation
blue particle’s curve ofcurve
aggregation IEPFMofisIEPFM
below, is and covered
below, andwith the green
covered with
curve of DEPFM in Figure 11. This means that blue particles tend to gather
the green curve of DEPFM in Figure 11. This means that blue particles tend to gather to a more to a more limiting and
actual target distribution, and IEPFM generates more diversity particles with large
limiting and actual target distribution, and IEPFM generates more diversity particles with large weights. In addition,
the MCMC
weights. In move stepthe
addition, is used
MCMCto reject
movethe particles
step is usedwhich move
to reject theoutside
particlesthe true posterior
which distributions
move outside the true
and
posterior distributions and accept the particles with excellent performance. Through themost
accept the particles with excellent performance. Through the process of IEPFM, of the
process of
particles with normalized weights less than 0.04 are eliminated, and the large-weight
IEPFM, most of the particles with normalized weights less than 0.04 are eliminated, and the large- particles are
preserved and evolved.
weight particles are preserved and evolved.

(a) (b)

Figure 11. Comparison of the particle distribution and normalized weights produced by DEPFM and
Figure 11.
IEPFM on Comparison of the particle
different assimilation days.distribution
(a) t = 45, and
(b) tnormalized weights
= 70. The blue and produced by DEPFM
green points present and
the
IEPFM on
particles differentby
generated assimilation
IEPFM anddays. (a) t respectively.
DEPFM, = 45, (b) t = 70.
TheThe blue and green
corresponding points
SM of present
red line the
denotes
particles
the generated by IEPFM and DEPFM, respectively. The corresponding SM of red line denotes
observation.
the observation.

The scatter diagram Figure 12 is selected for illustration of the particle’s diversity performance
in another way. From the figure, it is hard to find the phenomenon of particle impoverishment
characterized by many particles sharing the same values, clustering together at a few places, like erect
Water 2019, 11, 211 16 of 18

The scatter diagram Figure 12 is selected for illustration of the particle’s diversity performance
in another way. From the figure, it is hard to find the phenomenon of particle impoverishment
characterized by many particles sharing the same values, clustering together at a few places, like
erect short line segments. Because of crossover and mutation, particles could maintain the diversity
and avoid collapsing local optimums. Figure 12 also exhibits that blue particles scatter around the
Water 2019, 11 FOR PEER REVIEW 17
observations and their location center is closer to the observation than green particles. In summary,
after genetic operation and MCMC moves in IEPFM procedure, the diversity of particles are promoted.

(a) (b)

Figure 12. Scatter diagrams of the particles generated by IEPFM and DEPFM on different assimilation
Figure(a)
days. 12.t Scatter diagrams
= 45, (b) of the
t = 70. The particles
blue generated
and green circlesby IEPFMthe
present and DEPFMgenerated
particles on different
by assimilation
IEPFM and
days. (a) t = 45, (b) t = 70. The blue and green circles present the particles
DEPFM, respectively. The corresponding SM of red line denotes the observation. generated by IEPFM and
DEPFM, respectively. The corresponding SM of red line denotes the observation.
4. Conclusions
4. Conclusions
An IEPFM algorithm inspired by the biological immune system has been proposed in this paper to
alleviate the particle
An IEPFM impoverishment
algorithm inspired byproblem of theimmune
the biological general PF. The distinguishing
system has been proposed feature of IEPFM
in this paper
is
to aalleviate
special strategy
the particle of generating
impoverishment candidate particles
problem for general
of the MCMC PF. move Thestep based on thefeature
distinguishing immune of
evolution
IEPFM is concept.
a specialThe performance
strategy of the proposed
of generating candidateIEPFM in estimating
particles for MCMC themove
posterior
stepstate
based variables
on the
is validated
immune through
evolution the two
concept. case
The studies, including
performance the high
of the proposed dimensional
IEPFM and strongly
in estimating nonlinear
the posterior state
atmospheric model andthrough
variables is validated the VICthe hydrological model. This
two case studies, study obtains
including the highthe following findings.
dimensional and strongly
Both assimilation
nonlinear atmospheric results model demonstrate
and the VIC that IEPFM consistently
hydrological model. This receives the bestthe
study obtains performance
following
compared
findings. with EnKF and DEPFM. The candidates with more prior information may lead to more
accurate
Bothposterior
assimilation distribution of state variables
results demonstrate duringconsistently
that IEPFM particle evolution
receivesprocession. The ability
the best performance
of MCMC with
compared simulation
EnKF also improvesThe
and DEPFM. adequate
candidatesparticle
withdiversity,
more prior thusinformation
the particlemay impoverishment
lead to more
problem
accurate is mitigated.
posterior distribution of state variables during particle evolution procession. The ability of
MCMC An important
simulationand alsointeresting
improvesfindingadequateis that the performance
particle diversity, of IEPFM
thus the appears
particle almost unaffected
impoverishment
by the setting of
problem is mitigated. crossover probability and mutation probability, which are the key parameters of the
immune An evolution
importantalgorithm. Based on
and interesting this finding,
finding is thatcrossover and mutation
the performance suggested
of IEPFM probabilities
appears almost
are set to 0.8by
unaffected andthe0.3,setting
respectively. The benefit
of crossover of this approach
probability is in reducing
and mutation the influence
probability, which of aresubjective
the key
setting
parametersof keyofparameters
the immune on algorithm
evolution performance.
algorithm. Based on this finding, crossover and mutation
Additionally,
suggested the ensemble
probabilities are set to size of particles
0.8 and set to 100The
0.3, respectively. receives
benefitsimilar
of thisRMSE as the
approach larger
is in size.
reducing
This finding isofuseful
the influence for decreasing
subjective setting of computational
key parameterstime of IEPFM,performance.
on algorithm since the computational demand
increases rapidly when
Additionally, particle number
the ensemble exceeds this
size of particles value.
set to 100 receives similar RMSE as the larger size.
This Finally,
finding isituseful
is notfordifficult
decreasing to see that EnKFtime
computational behaves differently
of IEPFM, since the in computational
the two experiments.
demand
The performance is not improved significantly
increases rapidly when particle number exceeds this value. with the increase of ensemble size. The results indicate
that EnKF
Finally,hasit limitations in dealing
is not difficult to seewiththatstrongly nonlineardifferently
EnKF behaves systems. in the two experiments. The
performance is not improved significantly with the increase of ensemble size. The results indicate
that EnKF has limitations in dealing with strongly nonlinear systems.
However, IEPFM will be validated further with other observations, such as brightness
temperatures. Furthermore, some improvements could be made on IEPFM and VIC model. The
parameters of VIC model also impact the efficiency of assimilation. Thus, the optimization problem
of these model parameters will probably be investigated in the future.
Water 2019, 11, 211 17 of 18

However, IEPFM will be validated further with other observations, such as brightness
temperatures. Furthermore, some improvements could be made on IEPFM and VIC model.
The parameters of VIC model also impact the efficiency of assimilation. Thus, the optimization
problem of these model parameters will probably be investigated in the future.

Author Contributions: Conceptualization & methodology, F.J. and R.A.; validation, F.J. and Y.S.; investigation,
F.J. and Y.S.; software, F.J; writing—original draft preparation, F.J.; writing—review and editing, R.A. and Y.S.;
funding acquisition, R.A.
Funding: This research is funded jointly by the Key Project in the National Science & Technology Pillar Program
during the Twelfth Five-Year Plan Period, grant number 2013BAC03B04; the National Natural Science Foundation
of China, grant number 41871326 and 41271361.
Acknowledgments: The authors would like to thank Prof. Su and Dr. Zeng for kindly providing the in situ data.
The authors also thank the anonymous reviews for their critical reviews and constructive comments.
Conflicts of Interest: The authors declare no conflict of interest.

References
1. Parrens, M.; Wigneron, J.P.; Richaume, P.; Mialon, A.; Al Bitar, A.; Fernandez-Moran, R.; Kerr, Y.H.
Global-scale surface roughness effects at L-band as estimated from SMOS observations. Remote Sens.
Environ. 2016, 181, 122–136. [CrossRef]
2. Brocca, L.; Ciabatta, L.; Massari, C.; Camici, S.; Tarpanelli, A. Soil Moisture for Hydrological Applications:
Open Questions and New Opportunities. Water 2017, 9, 140. [CrossRef]
3. Kirby, J.M.; Mainuddin, M.; Ahmad, M.D.; Gao, L. Simplified Monthly Hydrology and Irrigation Water
Use Model to Explore Sustainable Water Management Options in the Murray-Darling Basin. Water Resour.
Manag. 2013, 27, 4083–4097. [CrossRef]
4. Lievens, H. SMOS soil moisture assimilation for improved hydrologic simulation in the Murray Darling
Basin, Australia. Remote Sens. Environ. 2015, 168, 146–162. [CrossRef]
5. Huang, C.; Xin, L.; Ling, L. Retrieving soil temperature profile by assimilating MODIS LST products with
ensemble Kalman filter. Remote Sens. Environ. 2008, 112, 1320–1336. [CrossRef]
6. Moradkhani, H. Hydrologic remote sensing and land surface data assimilation. Sensors 2008, 8, 2986–3004.
[CrossRef] [PubMed]
7. Dumedah, G.; Coulibaly, P. Evolutionary assimilation of streamflow in distributed hydrologic modeling
using in-situ soil moisture data. Adv. Water Resour. 2013, 53, 231–241. [CrossRef]
8. Weerts, A.H.; El Serafy, G.Y.H. Particle Filtering and Ensemble Kalman Filtering for state updating with
hydrological conceptual rainfall runoff models. Water Resour. Res. 2006, 42, 123–154. [CrossRef]
9. Evensen, G. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo
methods to forecast error statistics. J. Geophys. Res. Ocean. 1994, 99, 10143–10162. [CrossRef]
10. Gordon, N.J.; Salmond, D.J.; Smith, A.F.M. Novel approach to nonlinear/non-Gaussian Bayesian state
estimation. IEE Proc. F Radar Signal Process. 1993, 140, 107–113. [CrossRef]
11. Montzka, C.; Moradkhani, H.; Weihermüller, L.; Franssen, H.J.H.; Canty, M.; Vereecken, H. Hydraulic
parameter estimation by remotely-sensed top soil moisture observations with the particle filter. J. Hydrol.
2011, 399, 410–421. [CrossRef]
12. Mattern, J.P.; Dowd, M.; Fennel, K. Particle filter-based data assimilation for a three-dimensional biological
ocean model and satellite observations. J. Geophys Res. Ocean. 2013, 118, 2746–2760. [CrossRef]
13. Fan, Y.R.; Huang, G.H.; Baetz, B.W.; Li, Y.P.; Huang, K.; Li, Z.; Xiong, L.H. Parameter uncertainty and
temporal dynamics of sensitivity for hydrologic models: A hybrid sequential data assimilation and
probabilistic collocation method. Environ. Model. Softw. 2016, 86, 30–49. [CrossRef]
14. Hartanto, I.M.; Van Der Kwast, J.; Alexandridis, T.K.; Almeida, W.; Song, Y.; van Andel, S.J.; Solomatine, D.P.
Data assimilation of satellite-based actual evapotranspiration in a distributed hydrological model of a
controlled water system. Int. J. Appl. Earth Obs. Geoinf. 2017, 57, 123–135. [CrossRef]
15. Vrugt, J.A.; ter Braak, C.J.; Diks, C.G.; Schoups, G. Hydrologic data assimilation using particle Markov chain
Monte Carlo simulation: Theory, concepts and applications. Adv. Water Resour. 2013, 51, 457–478. [CrossRef]
16. Elsheikh, A.H.; Hoteit, I.; Wheeler, M.F. A nested sampling particle filter for nonlinear data assimilation. Q. J.
R. Meteorol. Soc. 2014, 140, 1640–1653. [CrossRef]
Water 2019, 11, 211 18 of 18

17. Plaza Guingla, D.A.; Pauwels, V.R.; De Lannoy, G.J.; Matgen, P.; Giustarini, L.; De Keyser, R. Improving
particle filters in rainfall-runoff models: Application of the resample-move step and the ensemble Gaussian
particle filter. Water Resour. Res. 2013, 49, 4005–4021. [CrossRef]
18. Bi, H.; Ma, J.; Wang, F. An Improved Particle Filter Algorithm Based on Ensemble Kalman Filter and Markov
Chain Monte Carlo Method. IEEE Trans. Geosci. Remote Sens. 2015, 53, 6134–6147. [CrossRef]
19. Dumedah, G.; Coulibaly, P. Integration of an evolutionary algorithm into the ensemble Kalman filter and the
particle filter for hydrologic data assimilation. J. Hydroinf. 2014, 16, 74–94. [CrossRef]
20. Mechri, R.; Ottlé, C.; Pannekoucke, O.; Kallel, A. Genetic particle filter application to land surface temperature
downscaling. J. Geophys. Res. Atmos. 2014, 119, 2131–2146. [CrossRef]
21. Yin, S.; Zhu, X. Intelligent Particle Filter and Its Application to Fault Detection of Nonlinear System.
IEEE Trans. Ind. Electron. 2015, 62, 3852–3861. [CrossRef]
22. Qu, Y.; Qian, X.; Song, H.; Xing, Y.; Li, Z.; Tan, J. Soil Moisture Investigation Utilizing Machine Learning
Approach Based Experimental Data and Landsat5-TM Images: A Case Study in the Mega City Beijing. Water
2018, 10, 423. [CrossRef]
23. Abbaszadeh, P.; Moradkhani, H.; Yan, H.X. Enhancing hydrologic data assimilation by evolutionary Particle
Filter and Markov Chain Monte Carlo. Adv. Water Resour. 2018, 111, 192–204. [CrossRef]
24. Zhang, H.J. State and parameter estimation of two land surface models using the ensemble Kalman filter
and the particle filter. Hydrol. Earth Syst. Sci. 2017, 21, 4927–4958. [CrossRef]
25. Vrugt, J.A. Markov chain Monte Carlo simulation using the DREAM software package: Theory, concepts,
and MATLAB implementation. Environ. Model. Softw. 2016, 75, 273–316. [CrossRef]
26. Arulampalam, M.S.; Maskell, S.; Gordon, N.; Clapp, T. A tutorial on particle filters for online
nonlinear/non-Gaussian Bayesian tracking. IEEE Trans. Signal Process. 2002, 50, 174–188. [CrossRef]
27. Dudek, G. An Artificial Immune System for Classification with Local Feature Selection. IEEE Trans. Evolut.
Comput. 2012, 16, 847–860. [CrossRef]
28. De Castro, L.N.; Timmis, J.I. Artificial immune systems as a novel soft computing paradigm. Soft Comput.
2003, 7, 526–544. [CrossRef]
29. Dasgupta, D.; Yu, S.H.; Nino, F. Recent Advances in Artificial Immune Systems: Models and Applications.
Appl. Soft Comput. 2011, 11, 1574–1587. [CrossRef]
30. Han, H.; Ding, Y.S.; Hao, K.R.; Liang, X. An evolutionary particle filter with the immune genetic algorithm
for intelligent video target tracking. Comput. Math. Appl. 2011, 62, 2685–2695. [CrossRef]
31. Wang, M.; Feng, S.; He, C.; Li, Z.; Xue, Y. An Artificial Immune System Algorithm with Social Learning and
Its Application in Industrial PID Controller Design. Math. Probl. Eng. 2017. [CrossRef]
32. Yang, Z.; Ding, Y.; Jin, Y.; Hao, K. Immune-Endocrine System Inspired Hierarchical Coevolutionary
Multiobjective Optimization Algorithm for IoT Service. IEEE Trans. Cybern. 2018. [CrossRef] [PubMed]
33. Andrieu, C.; Doucet, A.; Holenstein, R. Particle Markov chain Monte Carlo methods. J. R. Stat. Soc. 2010, 72,
269–342. [CrossRef]
34. Lorenz, E.N. Predictability: A problem partly solved. Proc. Semin. Predict. 1995, 1, 1–18.
35. Bryan, B.A.; Gao, L.; Ye, Y.; Sun, X.; Connor, J.D.; Crossman, N.D.; Liu, Z. China’s response to a national
land-system sustainability emergency. Nature 2018, 559, 193–204. [CrossRef] [PubMed]
36. Su, Z.; Wen, J.; Dente, L.; Velde, R.; Wang, L.; Ma, Y.; Hu, Z. The Tibetan Plateau observatory of plateau scale
soil moisture and soil temperature (Tibet-Obs) for quantifying uncertainties in coarse resolution satellite and
model products. Hydrol. Earth Syst. Sci. 2011, 15, 2303–2316. [CrossRef]
37. Liang, X.; Lettenmaier, D.P.; Wood, E.F.; Burges, S.J. A simple hydrologically based model of land surface water
and energy fluxes for general circulation models. J. Geophys. Res. Atmos. 1994, 99, 14415–14428. [CrossRef]

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access
article distributed under the terms and conditions of the Creative Commons Attribution
(CC BY) license (http://creativecommons.org/licenses/by/4.0/).