Nonlin. Processes Geophys., 14, 395–408, 2007
www.nonlin-processes-geophys.net/14/395/2007/© Author(s) 2007. This work is licensed
under a Creative Commons License.

Nonlinear Processes
in Geophysics

Merging particle filter for sequential data assimilation
S. Nakano1,2 , G. Ueno1,2 , and T. Higuchi1,2
1 The Institute of Statistical Mathematics, Research Organization of Information and Systems, Japan
2 Japan Science and Technology Agency, Japan

Received: 8 January 2007 – Revised: 25 April 2007 – Accepted: 5 July 2007 – Published: 16 July 2007

Abstract. A new filtering technique for sequential data assimilation, the merging particle filter (MPF), is proposed.
The MPF is devised to avoid the degeneration problem,
which is inevitable in the particle filter (PF), without prohibitive computational cost. In addition, it is applicable tocases in which a nonlinear relationship exists between a stateand observed data where the application of the ensemble
Kalman filter (EnKF) is not effectual. In the MPF, the filtering procedure is performed based on sampling of a forecast
ensemble as in the PF. However, unlike the PF, each member of a filtered ensemble is generated by merging multiple
samples from the forecast ensemble such that the mean andcovariance of the filtered distribution are approximately preserved. This merging of multiple samples allows the degeneration problem to be avoided. In the present study, the newly
proposed MPF technique is introduced, and its performanceis demonstrated experimentally.

1

Introduction

Data assimilation is performed to obtain the best estimatesof a state of a dynamic system or the evolution of a system
by incorporating observation into a model of the system andis used as an important tool for modeling and prediction ofgeophysical processes. Data assimilation methods are classified into two categories: variational data assimilation andsequential data assimilation. While variational data assimilation is performed by fitting a dynamic model to all of theavailable observations during a period of interest, sequential
data assimilation is an on-line approach that updates the estimation of a state at each observation time. In the present
study, we focus on sequential data assimilation.
Correspondence to: S. Nakano(shiny@ism.ac.jp)

Most sequential data assimilation techniques basicallyconsider a probability density function (PDF) of a state ofa dynamic system. An assimilation process is based on aprior PDF of the current state which is obtained using pastdata and a system model. This prior PDF is then updated toobtain the posterior PDF of the state by incorporating constraints based on observation. The procedure used to obtain
the posterior PDF is called “filtering”. The filtering procedure provides a PDF of the current state considering current
and past observations, which should be a basis for accurate
prediction of future states.
If a PDF of a state is Gaussian and the dynamics of thesystem is linear, then a filtering process can be described by
the algorithm of the Kalman filter. However, since geophysical systems usually contain inherent nonlinearity, it is rarethat the Kalman filter can be applied. The Kalman filter algorithm is sometimes extended by modifying the calculationof covariances of a state by linearizing a system model, andthis extended algorithm is called the extended Kalman filter(EKF). However, for models with high nonlinearity, the EKF
can make errors diverge (e.g., Evensen, 1992). Moreover, fora model with a large number of variables, the EKF requires
a high computational cost. Although the computational cost
could be reduced by using a variant of the EKF, the singular evolutive extended Kalman (SEEK) filter (Pham et al.,1998b), the SEEK filter also requires the linearization of asystem model and it can provide an unstable result for caseswith high nonlinearity.
In order to apply data assimilation to a system with nonlinear dynamics, it is practical to approximate a PDF of astate by an ensemble consisting of many realizations called“particles”. The ensemble Kalman filter (EnKF) (Evensen,
1994; Burgers et al., 1998) is one of such methods, and several variants of this algorithm have also been proposed (e.g.,Anderson, 2001; Whitaker and Hamill, 2002). The EnKF
is applicable to data assimilation of nonlinear systems. In
the EnKF, each particle in an ensemble is updated using a

Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.

396
Kalman gain calculated from the mean and the covariancesof the prior ensemble. However, the EnKF basically assumes
a linear relationship between a state and observed data in calculating a Kalman gain. Therefore, the EnKF does not provide good estimates of a state for cases in which linear approximation of the relationship between a state and observed
data is invalid. In addition, the computational cost of eachfiltering step in the EnKF is large due to repetitive multiplications and additions of matrices. Pham et al. (1998a) haveproposed another ensemble-based filtering method, the singular evolutive interpolated Kalman (SEIK) filter, which isderived as a variant of the SEEK filter. Although the SEIK
filter can work more efficiently than the EnKF (Nerger et al.,2005), it is not applicable to cases with nonlinear observation
as well.
The particle filter (PF) (Gordon et al., 1993; Kitagawa,
1993, 1996; Kitagawa and Gersch, 1996; Higuchi and Kitagawa, 2000; van Leeuwen, 2003), which is sometimes referred to as the sequential importance resampling (SIR) filter, is another method that is based on ensemble approximation of a PDF. In the PF, an estimation of a posterior PDF
is obtained by resampling with replacement from a prior ensemble. As the PF does not require assumptions of linearity
or Gaussianness, it is applicable to general nonlinear problems. In particular, the PF can be applied to cases in which
the relationship between a state and observed data is nonlinear, to which the application of the ensemble Kalman filter(EnKF) is not appropriate. However, the PF often encountersa problem called “degeneration”, which does not occur in the
EnKF. Since resampling procedures are applied recursively,
most of the particles are replaced by particles that fit the observed data better, and the posterior PDF is eventually represented by only a few of the particles among the membersof the initial ensemble. This reduces the validity of ensemble
approximation. This problem could be avoided by increasing
the number of particles in the ensemble. However, in order toincrease the number of particles, a prohibitive computational
cost is often required at each forecast step.
One potential way to avoid the degeneration problem is
to approximate a posterior distribution as a Gaussian distribution. This approach has been proposed by Kotecha andDjurić (2003) under the name of the Gaussian particle filter(GPF), and a similar algorithm was also presented by Anderson and Anderson (1999). In this technique, from an ensemble that represents a filtered posterior distribution, the meanand covariances are calculated to obtain a Gaussian distribution for approximating the filtered distribution. By drawing
random samples from this Gaussian distribution, a filteredensemble is newly generated. In the GPF, although the accuracy of an approximation of a filtered distribution is worse
than in the PF because of the assumption of Gaussianness,no duplicate particles are contained in the ensemble and degeneration does not occur. However, in generating Gaussian
random vectors to make a Gaussian ensemble, we must factorize the covariance matrix, which requires a high computaNonlin. Processes Geophys., 14, 395–408, 2007

S. Nakano et al.: Filtering for data assimilation
tional cost if the dimension of a state vector is large. In mostpractical cases, factorization of the covariance matrix withthe dimension of a state vector is not realistic.
There is another way to avoid degeneration which is a variant of the PF referred to as the kernel filter (Hürzeler andKünsch, 1998; Anderson and Anderson, 1999) or the regularized particle filter (Musso et al., 2001). This techniqueapproximates the filtered PDF by a sum of Gaussian functions with small standard deviations centered at the particle
locations, and members of a filtered ensemble is drawn fromthe sum of Gaussian functions. However, in applying thistechnique to high-dimensional models, there is difficulty indesigning a covariance matrix for each of the Gaussian functions. Although a covariance matrix could be made on thebasis of the covariance matrix of an ensemble representing
a prior or posterior PDF, this bring the same problem as the
GPF; that is, the factorization of the covariance matrix is required and the computational cost would become prohibitivein cases that a state vector is high-dimensional.
Thus, there exists no practical method to allow sequential
data assimilation with acceptable computational cost, except
some methods such as the EnKF and the SEIK filter which
also have a disadvantage in that it is not necessarily applicable to cases with nonlinear observations. To overcome
this problem, another technique, the merging particle filter(MPF), is devised. The MPF is an improved algorithm of thePF, in which filtering is performed by merging several particles of a prior ensemble, which is rather similar to the genetic
algorithm (e.g., Goldberg, 1989). This merging procedure allows the degeneration problem to be avoided and requires far
fewer particles than the PF. The primary advantage of the PF
over the EnKF is inherited; that is, the MPF is applicableeven to cases in which the relationship between a state andobserved data is nonlinear. Moreover, since the MPF does
not require the calculation of an inverse matrix, the computational cost at each filtering step is lower than that of the
EnKF. The PF algorithm, which the proposed algorithm isbased on, is reviewed in Sect. 2, and the MPF algorithm isintroduced in Sect. 3. In order to evaluate the performanceof the MPF, the results of a number of experiments are described in Sect. 4. Finally, the effectiveness of the MPF is
discussed and summarized in Sect. 5.

2Particle filter

The following state space model is considered:x k = Fk (x k−1 , v k )y k = Hk (x k ) + wk

(1a)
(1b)

where the vectors x k and y k indicate the state of a system andobserved data at a discrete time T =tk (k=1, . . .), respectively, and the vectors v k and wk denote system noise andobservation noise, respectively. The operator Fk represents
www.nonlin-processes-geophys.net/14/395/2007/

S. Nakano et al.: Filtering for data assimilation397

the temporal evolution of a state from time tk−1 to time tk according to the system model based on the simulation, while
Hk projects the state vector x k to the observation space.
The PF considers a PDF of a state x k , and the PDF isapproximated by an ensemble consisting of a large number of discrete samples called ‘particles’. For example, afiltered distribution at time T =tk−1 , p(x k−1 |y 1:k−1 ), is ap(N )(2)(1)
proximated by particles {x k−1|k−1 , x k−1|k−1 , · · · , x k−1|k−1 }
as
N


1 X
(i)
p(x k−1 |y 1:k−1 ) ≈δ x k−1 − x k−1|k−1
N i=1

(2)

where δ is Dirac’s delta function, and N is the number of particles in the ensemble. Here we expressedp(x k−1 |y 1 , · · · , y k−1 ) as p(x k−1 |y 1:k−1 ). From this ensemble approximation of p(x k−1 |y 1:k−1 ), we obtain an ensemble
approximation of the forecast distribution of the state at thenext observation time T =tk asp(x k |y 1:k−1 ) ≈

N


1 X
(i)
δ x k − x k|k−1 .
N i=1

(3)
(i)

Each particle of the forecast ensemble x k|k−1 is given by
(i)

(i)

p(x k |y 1:k )p(x k |y 1:k−1 ) p(y k |x k )=Rp(x k |y 1:k−1 ) p(y k |x k )dx k
N

 

X
1
(i)
(i)

p y k |x k|k−1 δ x k − x k|k−1
≈P 
(j )
j p y k |x k|k−1 i=1
N
X



(i)
wi δ x k − x k|k−1

i=1

(4)
(i)

(i)

for each x k|k−1 . From Eqs. (4) and (6), we obtain an approximation of p(x k |y 1:k ) using uniformly weighted particles, as
follows:p(x k |y 1:k ) ≈

N
X



(i)
wi δ x k − x k|k−1

i=1

≈

N

X
mi 
(i)
δ x k − x k|k−1
N
i=1

=

N


1 X
(i)
δ x k − x k|k .
N i=1

(7)

Thus, the newly generated ensemble approximates the filtered PDF p(x k |y 1:k ). Equation (7) has the same form as
Eq. (2), which allows us to recursively repeat the above procedure from Eq. (2) to Eq. (7). By repeating the procedure,a sequence of observed data is incorporated into the systemmodel.

(i)

Fk (x k−1|k−1 , v k ) where v k is a realization of the system
noise. This procedure is called the forecast step.
From the forecast distribution p(x k |y 1:k−1 ) and observed
data y k , we obtain a filtered PDF p(x k |y 1:k ) by using Bayes’
theorem, as follows:=(i)contain multiple copies of x k|k−1 belonging to the forecast
ensemble, and the number of copies mi becomesXmi ≈ N wi
mi = N; mi ≥ 0
(6)

(i)

where p(y k |x k|k−1 ) is the likelihood of x k|k−1 given the datay k and the weight wi is defined as(i)p(y k |x k|k−1 )wi = P.
(j )j p(y k |x k|k−1 )(5)This is called the filtering step.
Equation (4) shows that p(x k |y 1:k ) is approximated using
particles weighted by wi . Based on Eq. (4), we obtain a new(1)(N )ensemble {x k|k , · · · , x k|k } which approximates p(x k |y 1:k )
(1)

(N )

by resampling the forecast ensemble {x k|k−1 , · · · , x k|k−1 }with a weight of wi for each i. The new ensemble may
www.nonlin-processes-geophys.net/14/395/2007/

3

Merging particle filter

In the PF, a filtered ensemble generated through the resampling procedure contains multiple copies of particles withhigh likelihoods, and particles with low likelihoods are removed from the ensemble. Therefore, after repeating resampling several times, the diversity of the ensemble decreases
and eventually becomes insufficient for validly representing
a PDF. This problem can be avoided by increasing the number of particles. However, due to limited computational resources, it is often impossible to use a sufficient number ofparticles to repeat resampling several times. The MPF, which
we propose in this section, allows us to remake a filtered ensemble while restraining the reduction of its diversity.
The MPF is a modification of the PF. In the MPF, a filteredensemble is constructed based on samples from a forecast ensemble as in the PF. However, each particle of a filtered ensemble is generated as an amalgamation of multiple particles
from the forecast ensemble, which is rather similar to the genetic algorithm. Although this does not ensure that the shapeof the filtered PDF is preserved, the mean and covariance ofthe filtered PDF are approximately preserved (asymptoticallypreserved as the number of particles approaches infinity) in
generating a filtered ensemble.
A filtered ensemble is obtained as follows. When thenumber of particles to be merged is assumed to be n,
we draw n×N samples from the forecast ensemble withweights of wi in Eq. (5), and we thus obtain an ensem(1,1)(n,1)(1,N )(n,N )
ble: {x̂ k|k , · · · , x̂ k|k , · · · , x̂ k|k , · · · , x̂ k|k }. A subset
Nonlin. Processes Geophys., 14, 395–408, 2007

398

S. Nakano et al.: Filtering for data assimilation

(j,1)

(j,N)

{x̂ k|k , · · · , x̂ k|k } from the n × N samples forms an ensemble approximating the filtered PDF, which satisfies

variances given by the new ensemble becomeZp(x k |y 1:k ) ≈

N


1 X
(j,i)
δ x k − x̂ k|k
N i=1

(8)

(i)

n
X

(j,i)

αj x̂ k|k .

N
1 X
(i)
(i)
(x −µk|k )(x k|k −µk|k )T
N i=1 k|k
!
!T
N
n
n
X
X
1 X
(j1 ,i)(j2 ,i)
=
αj x̂
−µk|k
αj2 x̂ k|k −µk|k
N i=1 j =1 1 k|k
j2 =1
1
"
#"
#
N
n
n

 X

 T
X
1 X
(j1 ,i)(j2 ,i)
=
αj x̂
−µk|k
αj2 x̂ k|k −µk|k
N i=1 j =1 1 k|k
j =1

(9)

1

≈

=

j =1

N X
n
1 X

N i=1 j =1n
X≈
In order to ensure that the newly generated ensemble preserves the mean and covariances of the filtered PDF for
N→∞, the merging weights αj are set to satisfy

j =1
n
X

αj2

Z

2

αj2



T

(j,i)
(j,i)
x̂ k|k −µk|k x̂ k|k −µk|k(x k −µk|k )(x k −µk|k )T

j =1

Z

n
X

N


1 X
(i)
δ x k −x k|k dx k
N i=1

=

because it consists of N samples drawn from the forecast
ensemble with weights of wi , as was the case in obtaining
the filtered ensemble in the previous section. Next, we make(1)
(N )a new ensemble consisting of N particles {x k|k , · · · , x k|k } toapproximate p(x k |y 1:k ). Each particle in the new ensemble
is generated as a weighted sum of n samples from the n×N
sample set as:

x k|k =

(x k −µk|k )(x k −µk|k )T

N


1 X
(j,i)
δ x k −x̂ k|k dx k
N i=1

(x k −µk|k )(x k −µk|k )T p(x k |y 1:k ) dx k =6 k|k(12)where 6 k|k is the covariance matrix of p(x k |y 1:k ). Here, weused an approximation as

αj = 1

(10a)

N
1 X
(j ,i)(j ,i)(x̂ 1 − µk|k )(x̂ k|k2 − µk|k )T ≈ 0 (if j1 6 = j2 ),N i=1 k|kαj2 = 1(10b)which is justified because the two sets of samples(j ,1)(j ,N )
(j ,1)(j ,N ){x̂ k|k1 , · · · , x̂ k|k1 } and {x̂ k|k2 , · · · , x̂ k|k2 } are obtained
through independent random sampling and would not correlate with each other. Therefore, the ensemble obtained using
Eq. (9) affords an approximation of p(x k |y 1:k ) preserving
the mean and covariances as

j =1

where each αj is a real number. When the merging weights
satisfy Eq. (10a), the mean of the PDF approximated by the(1)(N )new ensemble {x k|k , · · · , x k|k } becomesp(x k |y 1:k ) ≈
Z
xk

N
X

1
(i)
δ x k − x k|k
N i=1




dx k

N
N X
n
1 X
1 X
(j,i)
(i)
x k|k =
αj x̂ k|k
N i=1
N i=1 j =1
" Z
#
n
N


X
1 X
(j,i)
=
αj x kδ x k − x̂ k|k dx k
N i=1
j =1
Z
n
X
≈
αj x k p(x k |y 1:k ) dx k

=

(11)

j =1

Z
=

x k p(x k |y 1:k ) dx k = µk|k

where µk|k is the mean of the filtered PDF p(x k |y 1:k ). In
addition, if the merging weights αj satisfy Eq. (10b), the coNonlin. Processes Geophys., 14, 395–408, 2007

N


1 X
(i)
δ x k − x k|k .
N i=1

(13)

The number of merged particles n can be chosen almost
arbitrarily. However, in order that the merging proceduremakes sense, n must be equal to or greater than 3. If n=1,
the weight α1 must be 1 in order to satisfy both Eqs. (10a)and (10b), which is obviously equivalent to the normal PF. If
n=2, then one of merging weights must be 1, and the othermust be 0, so as to satisfy both Eqs. (10a) and (10b). Thissetting is also equivalent to the normal PF, which means thatthe merging procedure does not make sense. Although there
is no upper limit for n, it is not necessary to set n to be large.
As shown in the next section, if none of the merging weights
are zero, we would greatly benefit by the merging procedure
even when n is as small as 3.
When n is equal to or greater than 3, there are infinite
allowable sets of the merging weights: {α1 , · · · , αn }. Although there is no definitive way to determine the values ofwww.nonlin-processes-geophys.net/14/395/2007/

S. Nakano et al.: Filtering for data assimilation

399

Forecast ensemble (N particles)

Resampling (N particles)

Forecast ensemble (N particles)
State xFig. 1. PF scheme. The value of a state x is on the horizontal axis assuming that the state x is scalar.

the weights, it would be preferable to set them such thatno two weights are equal to each other and that none of
the weights become zero in order to reinforce the diversityof the filtered ensemble. Under this setting, two duplicate(1)(N )particles in the filtered ensemble {x k|k , · · · , x k|k } can begenerated only from two identical sets of n merged parti(1)(N )cles drawn from the forecast ensemble {x k|k−1 , · · · , x k|k−1 },if duplicate particles are not contained in the forecast en(i)semble. When the probability that particle x k|k−1 is drawn
from the forecast ensemble is wi (0≤wi <1), the probability
(i1 )
(in )
that a sequence of n particles {x k|k−1, · · · , x k|k−1} is drawn
Qn
Qnn
is j =1 wij . Since j =1 wij ≤(max wi ) , the number of
duplicate particles contained in the filtered ensemble is, at
most, approximately N×(max wi )n for the MPF, while it is
N× max wi for the PF.
Figures 1 and 2 show schematically the respective procedures of the PF and the MPF when the number of merging
particles is set to be 3. In the PF, a filtered ensemble is simply obtained by resampling. In the MPF with 3 merging particles, after 3N particles are sampled from the forecast ensemble, the 3N particles are divided into N combinationsof 3 particles, and the 3 particles in each combination are
merged to obtain a new particle. Even from combinations of
the same 3 particles, different particles can be made with different sets of weights. Thus, the filtered ensemble obtainedwith the MPF contains diverse particles in comparison withthat obtained with the PF.

4
4.1

Numerical experiments
Lorenz 63 model

We performed a numerical experiment to test the MPF. Although this method is actually devised for data assimilation
for high-dimensional models, we first used a simple model,
the Lorenz 63 model (Lorenz, 1963), to investigate the behaviors of the method. The Lorenz 63 model is described by
the following equations:dx= −s(x − y)dtdy
= rx − y − xzdtdz
= xy − bz.
dt

(14b)
(14c)

In the conventional parameter setting, the three parameters
are set as follows: s=10, r=28, and b=8/3. One time stepin integrating the system equation was set to be 0.01.
Initially, we ran this model to generate a sequence of measurement data for this test. The data were generated every 20
time step with errors of a standard deviation of 2.0. It was assumed that all of the components of the state vector, x, y, andz, could be observed. In this situation, the observation vector
at each observation time resides in the same vector space as
the state vector.
The generated data were assimilated into the model using
the PF and the MPF. In this and the following experiments,
we assume additive system noise, and thus Eqs. (1a) and (1b)are rewritten as follows.
x k = F (x k−1 ) + v k
y k = H (x k ) + wk

www.nonlin-processes-geophys.net/14/395/2007/

(14a)

(15a)
(15b)

Nonlin. Processes Geophys., 14, 395–408, 2007

400

S. Nakano et al.: Filtering for data assimilation

Forecast ensemble (N particles)

Resampling (3N particles)

Merging

Forecast ensemble (N particles)
State xFig. 2. Scheme of the MPF, in which the number of merging particles is set to be 3. The value of a state x is on the horizontal axis assumingthat the state x is scalar.

where the subscript k in Fk and Hk is omitted because thesystem and observation models considered here are timeindependent. In applying the MPF, the number of merged
particles was set to n=3, and the weights αj were set as follows:3
4√
13 + 1α2 =8√13 − 1α3 = −
8

α1 =

(16a)
(16b)
(16c)

which satisfies Eqs. (10a) and (10b). In both the PF and the
MPF, we need to calculate the likelihood p(y k |x k ) wherey k is the observation vector (xko , yko , zko ), and x k is the statevector (xk , yk , zk ) at time T =tk . Assuming that observationnoise wk obeys a Gaussian distribution with zero mean anda diagonal covariance as diag(σ 2 , σ 2 , σ 2 ), the likelihood becomes"#1||y k − x k ||2p(y k |x k ) = √
exp −
(17)
2σ 2
2π σ
where we set σ =3. The system noise was assumed to be a
Gaussian noise with zero mean and a diagonal covariance as
Nonlin. Processes Geophys., 14, 395–408, 2007

diag(0.01, 0.01, 0.01). Particles of the forecast ensemble atthe initial time step (T =t1 ) were generated from a Gaussian
distribution where the mean was given by the value of thedata at the same time step and the standard deviation was 4.0
for each component.
Figure 3 shows the x-component of the state vector x k as
estimated by the MPF, where the number of particles was set
to N=64, and Fig. 4 shows that estimated by the PF, where
the number of particles was also set to N=64. Here, the estimate was given by the average over the ensemble members.
In each figure, the black squares indicate the test data thatwere assimilated into the model, and the red line indicates
the true trajectory of the state. Finally, the blue line indicates the state estimated through data assimilation. As seenin these figures, the MPF successfully estimated the state,while the estimate by the PF largely deviated from the true
state after around time step 6360. Figures 5 and 6 show thesame data as shown in Figs. 3 and 4, respectively, but are focused on the period from time step 6000 to time step 7000.
While the true state began to decrease after time step 6360,
the estimate by the PF began to increase, and the PF failed
to trace the true trajectory thereafter. Estimates by the MPF
also increased after time step 6360. However, this result was
improved by the filtering at time step 6380, after which the
MPF again successfully traced again the true state.
www.nonlin-processes-geophys.net/14/395/2007/

S. Nakano et al.: Filtering for data assimilation
30

401
30

State
MPF
Data

20

20

10

10

0

0

-10

-10

-20

-20

-30

0

2000

4000

6000

8000

State
MPF
Data

10000

-30
6000

6200

6400

6600

6800

7000

Step

Step

Fig. 3. Result of the experiment of data assimilation by the MPF
for the Lorenz 63 model. The number of particles was set to N=64.
The black squares indicate the test data that were assimilated into
the model. The red line indicates the true state of x, and the blueline indicates the estimation of x as a result of the data assimilation.

Fig. 5. Result of the experiment of data assimilation by the MPF
for the Lorenz 63 model from time step 6000 to time step 7000 for
the x-component.

30
30

State
PF
Data

20

20
10

10

0

0

-10

-10

-20

-20
-30

State
PF
Data

-30
6000
0

2000

4000

6000

8000

10000

Step

Fig. 4. Result of the experiment of data assimilation by the PF for
the Lorenz 63 model. The number of particles was set to N=64.
The black squares indicate the test data that were assimilated into
the model. The red line indicates the true state of x and the blue line
indicates the estimation of x as a result of the data assimilation.

In order to clarify why the PF failed to trace the true trajectory, histograms of the ensemble for x around time step
6360 are shown in Fig. 7. At time step 6340, a gap appeared
around −1<x<0 in the filtered ensemble in the result by the
PF. This gap expanded remarkably at the next forecast step,
resulting in a large gap in the forecast ensemble at time step
6360. While the x value of the true state was −0.125 at thistime step, as indicated by the dashed line in each panel, nomembers of the ensemble were distributed around the true
state. In contrast, no distinct gap appeared in the filtered ensemble at time step 6340 in the result obtained by the MPF.
Thus, there were only small gaps in the forecast ensemble atthe next time step.
www.nonlin-processes-geophys.net/14/395/2007/

6200

6400

6600

6800

7000

Step

Fig. 6. Result of the experiment of data assimilation by the PF for
the Lorenz 63 model from time step 6000 to time step 7000 for the
x-component.

We conducted experiments with various numbers of particles using the PF and the MPF. Table 1 shows the root-meansquare of deviations from the true state over 50 000 time steps
for all of the components for each experiment. In this table,
the results obtained using the EnKF (Evensen, 1994; Burgers et al., 1998), which is widely used for data assimilation,
are also displayed for reference. Even if the number of particles was increased into 128, the PF provided a worse estimation than the MPF. When the number of particles was
set to N=256, the estimates by the PF became as good as
those by the MPF. Although Kivman (2003) has pointed outthat the PF tends to provide better estimations than the EnKF,
this table shows that the EnKF yields lower errors when thenumber of particles is small. However, even in such cases,
the MPF gives better estimations than the EnKF.
Nonlin. Processes Geophys., 14, 395–408, 2007

402

S. Nakano et al.: Filtering for data assimilation
PF filtered (t = 6340)

MPF filtered (t = 6340)

50

50

40

40

30

30

20

20

10

10

0

0
-20

-15

-10

-5

0

5

10

15

20

-20

-15

-10

-5

x
PF forecast (t = 6360)

50

0

5

10

15

20

15

20

15

20

x
MPF forecast (t = 6360)

50

40

40

30

30

20

20

10

10

0

0
-20

-15

-10

-5

0

5

10

15

20

-20

-15

-10

-5

x
PF filtered (t = 6360)

50

0

5

10

x
MPF filtered (t = 6360)

50

40

40

30

30

20

20

10

10

0

0
-20

-15

-10

-5

0

5

10

15

20

x

-20

-15

-10

-5

0

5

10

x

Fig. 7. Histograms of the distribution of x in the ensemble around time step 6360. The left-hand panels show the distributions for the results
obtained by the PF, and the right-hand panels show the distributions for the results obtained by the MPF. The upper panels show the filtered
distributions at time step 6340. The middle panels show the forecast distributions at time step 6360. The lower panels show the filtered
distribution at time step 6360. In the middle and lower panels, for reference, the true state of x is indicated with dashed lines.

4.2

Lorenz 96 model

In order to evaluate the performance of the MPF for models
on higher dimension, we performed another experiment usNonlin. Processes Geophys., 14, 395–408, 2007

ing the Lorenz 96 model (Lorenz and Emanuel, 1998), which
is described by the following equations:dxj= (xj +1 − xj −2 )xj −1 − xj + f
dt

(18)

www.nonlin-processes-geophys.net/14/395/2007/

S. Nakano et al.: Filtering for data assimilation

403

Table 1. Root-mean-square deviations from the true state over
50 000 time steps for an experiment using the Lorenz 63 model.
3000

MPF

EnKF

4.55
3.87
0.87
0.86

1.00
0.91
0.92
0.91

1.34
1.29
1.29
1.29

4000

5000

Table 2. Root-mean-square deviations from the true state from time
step 3000 to time step 20 000 for an experiment using the Lorenz
96 model. Since the result has converged to the limit, we omitted
to calculate the deviations for N >8192 for the EnKF and those for
N>65 536 for the MPF.

N =128N =256
N =512
N=1024
N=2048
N=4096
N=8192
N=16 384
N=32 768
N=65 536
N =131 072
N =262 144

PF

MPF

EnKF

3.47
3.10
2.94
2.26
1.60
1.29
1.08
0.96
0.84
0.83
0.79
0.77

1.74
1.03
0.90
0.84
0.83
0.81
0.81
0.80
0.80
0.80
–
–

0.91
0.88
0.87
0.87
0.86
0.86
0.86
–
–
–
–
–

(for j 6 = 20)(for j = 20).

7000

9000

10000

0

10

20

30

40

j

(19a)
(19b)

After we iterated the model through 2000 time steps to allow fluctuations in the system to develop sufficiently, thedata were generated every 10 time steps with errors having a standard deviation of 1.5. It was assumed that we
can observe xj if j is an even number (j =2, . . . , 40); thatis, if half of the state variables are observed. In assimilating these test data, the system noise was assumed to be a
Gaussian noise with zero mean and a diagonal covariance asdiag(0.25, . . . , 0.25). Particles of the forecast ensemble atthe initial time step (T =t1 ) were generated from a Gaussian
distribution with mean 2.0 and variance 2.0 for each component. Again, in applying the MPF, the number of merged
particles was set to n=3, and the weights αj were set according to Eq. (16). The likelihood was calculated as follows:
www.nonlin-processes-geophys.net/14/395/2007/

6000

8000

for j =1, . . . , J . Here, x−1 =xJ −1 , x0 =xJ , and xJ +1 =x1 . In
this study, J was set to be 40; that is, the dimension of a statevector is 40. The forcing term f was set to be 8. One time
step was set to be 0.005. In order to generate data for theexperiment, we ran this model from the initial condition as
xj = 8.0
xj = 8.008

Time step

N=64
N=128
N=256
N=512

PF

Fig. 8. Result of the experiment of data assimilation by the MPF for
the Lorenz 96 model for every 1000 times step from 3000 to 10 000.
In this experiment, the number of particles was set to N=512. The
red and blue lines indicate the true state and the estimate by the
MPF, respectively.

"||y k − H x k ||2p(y k |x k ) = √
exp −
2σ 2
2π σ
1

#
(20)

where y k is the observation vector (y1,k , . . . , y20,k ) andσ was set to be 3. The operator H extracts the observable components from the state vector x k . Since we assume that we can observe xj for an even number of j ,H x k =(x2,k x4,k . . . x40,k )T .
Figures 8 and 9 show the estimation by the MPF and thatby the PF, respectively. In the experiments shown in thesefigures, the number of particles was set to N=512. The abscissa indicates j , and the value of xj for each j for every
1000 time step from 3000 to 10 000 is shown in these figures. As shown in Fig. 9, the PF often deviates from the true
state (e.g., at time step 9000). On the other hand, the MPF
successfully estimates the state over the period shown here.
Table 2 shows the root-mean-square of the deviations from
Nonlin. Processes Geophys., 14, 395–408, 2007

404

S. Nakano et al.: Filtering for data assimilation
result using the PF is gradually improved, and the root-meansquare of the deviations for the PF seems to converge to aslightly better value than that for the MPF, probably because
the MPF does not preserve the shape of the PDF while thePF can faithfully preserve the shape of the filtered PDF with
abundant particles. For cases that N is larger than 262 144,
we did not perform experiments because they need too muchcomputational resources, and we could not confirm the value
which the root-mean-square deviation for the PF converged
to. Thus, the result of the PF with a further large ensemble
size possibly converges to a further good value than that for
N=262 144. However, the use of such an enormous numberof particles is not realistic, and it seems to provide only minor
improvement of the estimation accuracy even if it were possible. For practical applications to high-dimensional systems,
the use of the MPF or the EnKF with much fewer particles
would be effectual.

3000

4000

Time step

5000

6000

7000

8000

4.3
9000

10000

0

10

20

30

40

j
Fig. 9. Result of the experiment of data assimilation by the PF for
the Lorenz 96 model for every 1000 time steps from 3000 to 10 000.
In this experiment, the number of particles was set to N=512. The
red and blue lines indicate the true state and the estimate by the PF,
respectively.

the true state from time step 3000 to time step 20 000 forvarious numbers of particles. Again, for reference purposes,
the results using the EnKF are also shown in this table. We
omitted the calculation of the deviations for N>8192 for the
EnKF and those for N>65 536 for the MPF which requires
much computational resources and cost, because the valueof the root-mean-square deviation has converged to the limit
and the estimate would not be improved any more even if N
increased.
When N is small, the MPF fails to estimate the state, while
the EnKF achieves a robust estimation of the state. However,
the estimation accuracy of the MPF is remarkably improved
when N=256, and it becomes better than that of the EnKF
when N ≥1024. In comparison with the PF, the MPF provides good estimates without requiring a large number of particles. In this experiment, the MPF requires only 1024 particles to obtain as good accuracy as the PF with 32 768 particles. As the number of ensemble members N increases, theNonlin. Processes Geophys., 14, 395–408, 2007

Lorenz 96 model with nonlinear observation

Another experiment was performed to examine whether the
MPF works for the Lorenz 96 model with nonlinear observations. In this experiment, we assumed that we can observe only an absolute value |xj | if j is an even number(j =2, . . . , 40). The data were generated every 10 time steps
by taking the absolute values of xj containing errors with astandard deviation of 1.5. As in the previous experiment, thesystem noise was assumed to be a Gaussian noise with zeromean and a diagonal covariance as diag(0.25, . . . , 0.25), andparticles of the forecast ensemble at the initial time step
(T =t1 ) were generated from a Gaussian distribution withmean 2.0 and variance 2.0 for each component. The number and the weights of the merged particles in applying the
MPF were also the same as in the previous experiment. The
likelihood was calculated as follows:#"||y k − H (x k )||2
1
(21)
exp −
p(y k |x k ) = √
2σ 2
2π σ
where H (x k )=(|x2,k | |x4,k | . . . |x40,k |)T and σ =3.
Figures 10 and 11 show the estimation by the MPF andthat by the PF, respectively. In the experiments shown inthese figures, the number of particles was set to N=1024.
The abscissa indicates j , and, again, the value of xj for eachj for every 1000 time step from 3000 to 10 000 is shown. As
shown in Fig. 10, the MPF successfully estimates the state.
On the other hand, the PF failed to estimate the state. Table 3
shows the root-mean-square of deviations from the true statefrom time step 3000 to time step 20 000 for various numbersof particles. The results using the EnKF are also shown inthis table again. Here, it should be noted that the algorithm of
the EnKF must be modified to apply to cases with nonlinear
observations because the EnKF basically assumes a linear
relationship between a state and observed data. In applying
the EnKF to this particular experiment, according to Evensen
www.nonlin-processes-geophys.net/14/395/2007/

405

3000

3000

4000

4000

5000

5000

6000

6000

Time stepTime step

S. Nakano et al.: Filtering for data assimilation

7000

7000

8000

8000

9000

9000

10000

10000

0

10

20

30

40

Fig. 10. Result of the experiment of data assimilation by the MPF
for the Lorenz 96 model with nonlinear observations for every 1000
time steps from 3000 to 10 000. In this experiment, the number ofparticles was set to N=1024. The red and blue lines indicate thetrue state and the estimate by the MPF, respectively.

(2003), we define a new state vector x 0k =[x Tk , (H (x k ))T ]Tsuch that the observation model becomes linear, and the state
space model in Eqs. (15a) and (15b) is accordingly rewritten
into a new state space model as follows:
x 0k = F 0 (x 0k−1 , v k )y k = H 0 x 0k + wk .
Here the operators F 0 and H 0 are defined as:
F 0 (x 0k−1 , v k ) = F 0 x k−1 , H (x k−1 ), v kF (x k−1 ) + v k =
H F (x k−1 ) + v k
H 0 = (Odimx k Idimy k )

(22a)
(22b)

(23a)

(23b)

where Odimx k is a zero matrix whose dimension is the sameas x k and Idimy k is an identity matrix whose dimension iswww.nonlin-processes-geophys.net/14/395/2007/

0

10

20

30

40

j

j

Fig. 11. Result of the experiment of data assimilation by the PF
for the Lorenz 96 model with nonlinear observations for every 1000
time steps from 3000 to 10 000. In this experiment, the number ofparticles was set to N=1024. The red and blue lines indicate thetrue state and the estimate by the PF, respectively.

the same as y k , and thus H 0 extracts H (x k ) from the vector
x 0k . The EnKF is then applied to this new state space model.
Since the results have converged to the limit, we omitted to
calculate the deviations for N >8192 for the EnKF and thosefor N >65 536 for the MPF.
As shown in this table, if an enormous number of particles
are not allowed, the MPF provides much better results than
the PF. The MPF requires only 1024 particles to achieve better accuracy than the PF with 32 768 particles, as well as theprevious experiment. With an enormous number of particles,
the PF apparently provides better results than the MPF. In
this experiment in which only absolute values are allowed tobe observed, the filtered PDF may often have multiple modes
and moments of higher order than the second moment could
then be significant. This situation would limit the accuracy of
the MPF, which preserves only the first two moments, even
if infinite ensemble members are used. However, the PF
requires more than at least 65 536 particles to obtain better
Nonlin. Processes Geophys., 14, 395–408, 2007

406

S. Nakano et al.: Filtering for data assimilation

Table 3. Root-mean-square deviations from the true state from time
step 3000 to time step 20 000 for an experiment using the Lorenz 96
model with nonlinear observation. Since the result has converged to
the limit, we omitted to calculate the deviations for N >8192 for the
EnKF and those for N >65 536 for the MPF.
MPF

EnKF

4.17
4.01
3.66
3.70
3.15
2.65
2.07
1.80
1.23
1.19
1.04
1.00

3.56
2.47
1.50
1.20
1.19
1.14
1.14
1.13
1.13
1.13
–
–

1.75
1.94
1.93
1.98
1.99
1.99
1.99
–
–
–
–
–

4000

5000

Time step

N=128
N=256
N=512
N=1024
N=2048
N=4096
N=8192
N =16 384
N =32 768
N =65 536
N=131 072
N=262 144

PF

3000

6000

7000

8000

9000

accuracy than the MPF. Thus, as far as the number of particles is not allowed to increased to more than at least 65 536,
the degeneration problem of the PF is more serious than the
problem concerning high order moments of the MPF in thisexperiment. In comparison between the MPF and the EnKF,
when N is small, the EnKF provides better estimations again,
although estimations by the EnKF are not so good. WhenN ≥512, the estimation accuracy of the MPF is remarkably
improved to be much better than that of the EnKF. Actually,
the EnKF does not effectively work in this experiment. Figure 12 shows the estimation by the EnKF, where the numberof particles was set to N=1024. It is indicated that the estimates by the EnKF often significantly deviate from the true
state which means that the EnKF fails to capture the variation of the true state. Thus, for this experiment, the use of the
MPF would be the most effectual.
5

Summary and discussion

We proposed a new algorithm, the MPF, for realizing practical sequential data assimilation. The MPF provides anensemble-based approximation of the filtered PDF such thatthe mean and covariance are approximately preserved. The
MPF allows the problem of degeneration, which occurs in
the PF, to be avoided. It must be noted that the MPF does
not preserve the shape of the filtered PDF while the PF canfaithfully preserve the shape of the filtered PDF with abundant particles. Therefore, if a sufficient number of particles
is used, the PF should provide a better estimation than the
MPF. In particular, in cases that the filtered PDF is significantly non-Gaussian, the MPF possibly provides a rather bad
estimate. In application to a high-dimensional system, however, it is not realistic to use a sufficient particles to avoid deNonlin. Processes Geophys., 14, 395–408, 2007

10000

0

10

20

30

40

j

Fig. 12. Result of the experiment of data assimilation by the EnKF
for the Lorenz 96 model with nonlinear observations for every 1000
time steps from 3000 to 10 000. In this experiment, the number ofparticles was set to N=1024. The red and blue lines indicate thetrue state and the estimate by the EnKF, respectively.

generation, and therefore the PF should fail to approximate
the filtered PDF. Indeed, as illustrated in Sect. 4.2, the PF
provides a worse estimation of the state than the MPF for the
Lorenz 96 model, until the number of particles in the ensemble was increased to at least 65 536. Since usual geophysicalmodels are of much higher dimension than the Lorenz 96
model, although they could be less nonlinear, a hopelessly
large number of particles would be required in order to use
the PF. The MPF requires far fewer particles than the PF andthus would be a more effectual algorithm.
In addition, the MPF is applicable to cases in which the relationship between a state and observed data is nonlinear. For
cases with nonlinear observations, the EnKF does not necessarily provide a good estimation of the state. As illustrated
in Sect. 4.3, even if the number of particles is increased, theestimation by the EnKF is not improved, whereas that by the
MPF is remarkably improved. Therefore, the MPF would be
the best method of sequential data assimilation with nonlinear observations.
www.nonlin-processes-geophys.net/14/395/2007/

S. Nakano et al.: Filtering for data assimilation

407

Table 4. Comparison among the algorithms for sequential data assimilation with a high-dimensional nonlinear system.

Nonlinear observation
Necessary number of particles
Cost of filtering

PF

MPF

EnKF

OK
Exceedingly many
Low

OK
MediumLow

Ineffectual for some cases
Relatively few
High

For cases in which the relationship between a state in thesystem and observed data is linear, the EnKF basically provides a good estimation without a large number of particles.
However, the EnKF tends to require a higher computational
cost at each filtering step in applying to a high-dimensionalmodel, because it involves many multiplications and additions between matrices. In addition, even if the numberof particles is taken to be small, estimates using the EnKF
can be affected by spurious correlations between distant locations, and thus localization on the covariance matrix (Ottet al., 2004) might be required to avoid this problem. On theother hand, a computational cost at each filtering step is notserious in the MPF, because neither iterative calculations of
inverse matrices nor numerous multiplications between matrices are required. Therefore, for cases in which a systemmodel does not require a great deal of computational time,
the MPF may perform better than the EnKF.
Table 4 summarizes the characteristics of the algorithmsof the PF, the MPF, and the EnKF. In the cases of a nonlinear relationship between a state and observed data, the EnKF
does not necessarily work, whereas the PF or the MPF canbe applied. The PF requires an exceedingly large numberof particles, which imposes prohibitive computational cost at
each forecast step. The MPF requires far fewer particles than
the PF, although the EnKF requires fewer particles than the
MPF. As for the computational cost at each filtering step, the
EnKF requires a larger computational cost than the PF andthe MPF. The high computational cost at each filtering stepwould become serious in the case that the number of assimilated data is large. On the other hand, the increase in thenumber of particles causes a high computational cost at each
forecasting step, which becomes serious for the case in whicha system model requires a great deal of computational time.
Therefore, for the case in which only linear observations areused, the choice between the MPF and the EnKF should bemade based on the considerations of the dimension of the observation vector and the computational cost required by thesystem model.
Acknowledgements. This study was supported by the Japan
Science and Technology Agency (JST) under the Core Research
for Evolutional Science and Technology (CREST) program, andpartially supported by the Transdisciplinary Research Integration
Center, Research Organization of Information and Systems(ROIS/TRIC) as a Function and Induction Research Project.

www.nonlin-processes-geophys.net/14/395/2007/

Edited by: O. TalagrandReviewed by: P. J. van Leeuwen and two other anonymous referees

References
Anderson, J. L.: A ensemble adjustment Kalman filter for data assimilation, Mon. Wea. Rev., 129, 2884–2903, 2001.
Anderson, J. L. and Anderson, S. L.: A Monte Carlo implementation of the nonlinear filtering problem to produce ensemble
assimilations and forecasts, Mon. Wea. Rev., 127, 2741–2758,
1999.
Burgers, G., van Leeuwen, P. J., and Evensen, G.: Analysis scheme
in the ensemble Kalman filter, Mon. Wea. Rev., 126, 1719–1724,
1998.
Evensen, G.: Using the extended Kalman filter with a multilayer
quasi-geostrophic model, J. Geophys. Res., 97(C11), 17 905–
17 924, 1992.
Evensen, G.: Sequential data assimilation with a nonlinear quasigeostrophic model using Monte Carlo methods to forecast error
statistics, J. Geophys. Res., 99(C5), 10 143–10 162, 1994.
Evensen, G.: The ensemble Kalman filter: theoretical formulation and practical implementation, Ocean Dynam., 53, 343–367,
doi:10.1007/s10236-003-0036-9, 2003.
Goldberg, D. E.: Genetic algorithms in search, optimization andmachine learning, Addison-Wesley, Reading, 1989.
Gordon, N. J., Salmond, D. J., and Smith, A. F. M.: Novel approach
to nonlinear/non-Gaussian Bayesian state estimation, IEE Proceedings F, 140, 107–113, 1993.
Higuchi, T. and Kitagawa, G.: Knowledge discovery and selforganizing state space model, IEICE Transactions on Information and Systems, E83-D, 36–43, 2000.
Hürzeler, M. and Künsch, H. R.: Monte Carlo approximations forgeneral state space models, J. Comp. Graph. Statist., 7, 175–191,
1998.
Kitagawa, G.: Monte Carlo filtering and smoothing method for nonGaussian nonlinear state space model, Inst. Statist. Math. Res.
Memo., 1993.
Kitagawa, G.: Monte Carlo filter and smoother for non-Gaussiannonlinear state space models, J. Comp. Graph. Statist., 5, 1–25,
1996.
Kitagawa, G. and Gersch, W.: Smoothness priors analysis of timeseries, chap. 6, Springer-Verlag, New York, 1996.
Kivman, G. A.: Sequential parameter estimation for stochastic systems, Nonlin. Process. Geophys., 10, 253–259, 2003.
Kotecha, J. H. and Djurić, P. M.: Gaussian particle filtering, IEEE
Trans. Signal Processing, 51, 2592–2601, 2003.
Lorenz, E. N.: Deterministic nonperiodic flow, J. Atmos. Sci., 20,
130–141, 1963.

Nonlin. Processes Geophys., 14, 395–408, 2007

408
Lorenz, E. N. and Emanuel, K. A.: Optimal sites for supplementary
weather observations: Simulations with a small model, J. Atmos.
Sci., 55, 399–414, 1998.
Musso, C., Oudjane, N., and Le Gland, F.: Improving regularizedparticle filters, in: Sequential Monte Carlo methods in practice,edited by Doucet, A., de Freitas, N., and Gordon, N., chap. 12,
p. 247, Springer-Verlag, New York, 2001.
Nerger, L., Hiller, W., and Schröter, J.: A comparison of error subspace Kalman filters, Tellus, 57A, 715–735, 2005.
Ott, E., Hunt, B. R., Szunyogh, I., Zimin, A. V., Kostelich, E. J.,
Corazza, M., Kalnay, E., Patil, D. J., and Yorke, J. A.: A local
ensemble Kalman filter for atmospheric data assimilation, Tellus,
56A, 415–428, 2004.

Nonlin. Processes Geophys., 14, 395–408, 2007

S. Nakano et al.: Filtering for data assimilation
Pham, D. T., Verron, J., and Gourdeau, L.: Singular evolutive
Kalman filters for data assimilation in oceanography, C. R. Acad.
Sci. Ser. II, 326, 255–260, 1998a.
Pham, D. T., Verron, J., and Roubaud, M. C.: A singular evolutive
extended Kalman filter for data assimilation in oceanography, J.
Mar. Syst., 16, 323–340, 1998b.
van Leeuwen, P. J.: A variance-minimizing filter for large-scale applications, Mon. Wea. Rev., 131, 2071–2084, 2003.
Whitaker, J. S. and Hamill, T. M.: Ensemble data assimilation without perturbed observations, Mon. Wea. Rev., 130, 1913–1924,
2002.

www.nonlin-processes-geophys.net/14/395/2007/