Journal of Japan Society of Civil Engineers, Ser.
B1of(HydraulicVol. 67,
No. 4,2011,
I_1-I_6,
2011.
Annual Journal
Hydraulic Engineering),
Engineering, JSCE,
Vol.55,February

DUAL STATE-PARAMETER UPDATING SCHEME
ON A CONCEPTUAL HYDROLOGIC MODELUSING SEQUENTIAL MONTE CARLO FILTERS
Seong Jin NOH1, Yasuto TACHIKAWA2, Michiharu SHIIBA3 and Sunmin KIM4
1Student Member of JSCE, Graduate Student, Dept. of Urban and Env. Eng., Kyoto University(Nishikyo-ku, Kyoto 615-8540, Japan) E-mail:seongjin.noh@gmail.com
2Member of JSCE, Dr. Eng., Associate Professor, Dept. of Civil and Earth Resources Eng., Kyoto University(Nishikyo-ku, Kyoto 615-8540, Japan)
3Member of JSCE, Dr. Eng., Professor, Dept. of Civil and Earth Resources Eng., Kyoto University(Nishikyo-ku, Kyoto 615-8540, Japan)
4Member of JSCE, Dr. Eng., Lecturer, Dept. of Civil and Earth Resources Eng., Kyoto University(Nishikyo-ku, Kyoto 615-8540, Japan)

Applications of data assimilation techniques have been widely used to improve upon thepredictability of hydrologic modeling. Among various data assimilation techniques, sequential Monte
Carlo (SMC) filters, known as “particle filters”, provide the capability to handle non-linear andnon-Gaussian state-space models. This paper proposes a dual state-parameter updating scheme (DUS)
based on SMC methods to estimate both state and parameter variables of a hydrologic model. We
introduce a kernel smoothing method for the robust estimation of uncertain model parameters in the DUS.
The applicability of the dual updating scheme is illustrated using the implementation of the storage
function model on a middle-sized Japanese catchment. We also compare performance results of DUS
combined with various SMC methods, such as SIR, ASIR and RPF.

Key Words: Data assimilation, sequential Monte Carlo, particle filter, storage function modelin various areas, including pattern recognition,
target tracking, financial analysis, and robotics.
Only in recent years has the application of thesemethods been included in hydrology research2),3),4).
In the practical use of hydrologic models,estimated states are highly sensitive to theuncertainty of model parameters. Furthermore,
there is no guarantee that parameters calibrated
from previous data are the optimum in the current
prediction. Therefore, updating state variables
based on inappropriate parameters will likely
increase uncertainty in the forecasting of hydrologic
models. In this respect, sequential estimates of theparameters and state variables are needed to enable
the model to generate accurate forecasts.
In this paper, we propose a dual state-parameter
updating scheme (DUS) based on the SMC filters
for the estimation of both the state and parametervariables of a hydrologic model. A kernelsmoothing method is introduced for the robust
estimation of uncertain model parameters in the

1. INTRODUCTION
Identification and minimization of uncertainty
are key issues in the hydrologic prediction. Dataassimilation is a way to integrate information froma variety of sources to improve model accuracy,
considering the uncertainty in both themeasurement and modeling system. Among data
assimilation techniques, the sequential Monte Carlo(SMC) methods are a Bayesian learning process inwhich the propagation of all uncertainties is carriedout by a suitable selection of randomly generated
particles without any assumptions about the natureof the distributions. Unlike the various Kalmanfilter-based methods that are basically limited to thelinear system equation and the assumption ofGaussian distribution errors, SMC filters have theadvantage of being applicable to non-linear,
non-Gaussian, state-space models. Since theirintroduction in 19931), the application of thesepowerful and versatile methods has been increasing

I_1

Propagate statewith uniform weight

DUS. We illustrate its applicability for hydrologic
forecasting on a middle-sized Japanese catchment
using a conceptual hydrologic model.
The paper is organized in the following way.
Section 2 outlines the Bayesian filtering theory; thesequential Monte Carlo filters, known as particle
filters, which are based on the sequential importancesampling (SIS); and parameter inference
approaches in SMC. In Section 3, the case study
demonstrating the applicability of the SMC filters ispresented. The SMC filters are applied for real-time
forecasting of river discharge of the Katsura River
catchment using the storage function (SF) model.
Sequential data assimilation is performed by two
different schemes via the SMC filters: state only
updating and dual state-parameter updating.
Comparisons of the performance results of various
SMC filters are presented. Section 4 summarizes
the methodology and the analysis results.

{xt 1|t 1 , 1 / n}

Estimate likelihoodusing measurement

Update weights{xt|t 1 , wt }Resample

Predict the next step
{xt |t , 1 / n}Fig. 1 A single cycle of SMC.

p( xt | y1:t ) n w  (x  x )
i
t

t

i
t

(3)

i 1

where xti and wti denote the ith posterior state(“particle”) and its weight, respectively, and  ()denotes the Dirac delta function.
Since it is usually impossible to sample from thetrue posterior PDF, an alternative is to sample froma proposal distribution, also called importancedensity, denoted by q( xt | yt ) . After the several
steps of computation, the recursive weight updating
can be derived as follows:

2. METHODOLOGY(1) Bayesian filtering theory
To define the problem of the Bayesian filtering,consider a generic dynamic state-space model
which can be described as follows:(1)xt  f ( xt 1, , ut )  t t ~ N (0,Wt )(2)
yt  h( xt , )   t
 t ~ N (0,Vt )wti  wti1p( yt | xti ) p( xti | xti1 )q( xti | xti1, yt )(4)

The choice of importance density is one of themost critical issues in the design of SMC filters6).
The most popular choice is the transitional prior as(5)q( xti | xti1 , yt )  p( xti | xti1 )By substituting Eq. (5) into Eq. (4), the weight
updating becomes(6)
wti  wti1 p( yt | xti )With these particles and associated weights, theestimated state vector x̂t is the weighted mean ofparticles as:

where xt  n x is the nx dimensional vector denoting
the system state at time t. The operator
f : nx  nx expresses the system transition in
response to the forcing data ut, and parameters θ.
h : n x  nx expresses the measurement function.
ωt and νt represent the model and the measurement
error, and Wt and Vt represent the covariance of theerror, respectively. In particular, we seek filteredestimates of xt based on the set of all available
measurements y1:t = {yi, i = 1, …, t}. In the Bayesianrecursive estimation, if the system and measurement
models are non-linear and non-Gaussian, it is notpossible to construct the posterior probability density
function (PDF) of the current state xt given all themeasurement analytically. When the analytic
solution is intractable, an optimal solution can beapproximated by the SMC filters explained in thenext section.

xˆt 

n

w x

i i
t t |t 1

(7)

i 1

The SIS algorithm is a Monte Carlo method that
forms the basis for most SMC filters. A common
problem with the SIS algorithm is the degeneracy
phenomenon, where after a few iterations, all but
one particle will have negligible weight. The
degeneracy phenomenon can be reduced by
performing the resampling step whenever asignificant degeneracy is observed. A graphicalrepresentation of SMC is illustrated in Fig. 1. At the
top we start with a uniformly weighted randommeasure. Then we use the received measure yt tocompute its importance weight of each particle. If
necessary, a resampling step is executed to selectimportant particles with a uniform weight. If thenumber of particles is n, the weight is 1/n. The last
step is a prediction introducing process noise.

(2) Sequential importance sampling (SIS)
Sequential Monte Carlo (SMC) filters are a set ofsimulation-based methods that provide a flexible
approach to computing the posterior distribution
without any assumptions about the nature of the
distributions. The key idea of SMC is based on
point mass (“particle”) representations ofprobability densities with associated weights as5):

I_2

(3) Variant of SMC filters
Several variants of SMC filters have been
proposed in the literature to overcome thedegeneracy and sample impoverishment and to
improve selection of importance density. The
sample importance resampling (SIR) filter isderived from the SIS algorithm by performing theresampling step at every time index. The auxiliary
SIR (ASIR) filter performs the resampling step at
the previous time step, attempting to mimic the
optimal importance density. The regularizedparticle filter (RPF) was suggested as a method toimprove the sample diversity. It is worth noting thatthese filters can be (and often are) combined6).

Start
Initialize the system by sampling the model state andparameter vectors
Assign the particle weights uniformly
wti  1 / n i  1,, nSample the parameters from ith normal componentof the kernel density
ti ~ N (. | ati1  (1  a)t 1 , h2Vt1 ) i  1,, nPropagate the n model states forward in timethrough model operator f(.)
xti  f ( xti1 ,ti , ut )  t t ~ N (0,Wt ) i  1,, nPredict the system output through the operator h(.)
yti  h( xti ,ti )  t  t ~ N (0,Vt ) i  1,, n

Estimate the likelihoodp[ yt  h( xti ,ti ) | Vt ]Forecast based on
updated particles

i
t 1

t

i
t 1

2

t 1

i
t 1

i
t

t

i
t

i  1,, nt

Implement the resampling for states andparameterst= t +1

t<T
Stop

Fig. 2 Flowchart of the dual state-parameter updating schemewith kernel smoothing via the SIR particle filter.

Japan, and covers an area of 1,100 km2 (887 km2 at
the Katsura station). There are 13 rainfall
observation stations and 6 river flow observation
stations. The Hiyoshi dam is located upstream,
and the outflow record from that reservoir has been
considered to be input data in a hydrologic model.
(2) Hydrological model and simulation condition
The storage function (SF) model9) is one of themost commonly used conceptual hydrologic models
for flood prediction due to its simple numerical
procedure and its proper regeneration of nonlinear
characteristics of flood runoff. The state-space formof the SF model adapted in this catchment is asfollows:n w N ( | m , h V  )n w p[ y  h( x , ) | V ]i 1(4) Parameter inference
Identification of parameter uncertainty isessential to obtain unbiased data assimilation. To
handle inference of the unknown parameters, theconcept of “artificial evolution” can be applied.
That means that the parameter vector  isfluctuated at each time step, adding an independent,
zero-mean normal increment as follows:(8)t  t 1   t
 t ~ N (0, s 2Vt1 )
where  t is random noise, Vt 1 is the variance ofparameter particles at time t-1 before resampling,and s is a small tuning parameter. The drawbackof this approach is that estimated posterior
distribution of parameters becomes more diffuse
compared to the actual ones7). Kernel smoothing8) is
one remedy for this problem and is accomplishedby determining the covariance of parameters basedon particles from previous time points. The smoothkernel density can be a mixture of Gaussian
densities as follows:P( t | y1 ,, yt 1 ) ~Update the particle weights
wti1 p[ yt  h( xti , ti ) | Vt ]

wti 

(9)

i 1

1/ p

ds(t )
 s(t ) 
 re (t  TL )  

dt
 k 

where h is the variance reduction parameter. The
kernel locations mti1 are specified by a shrinkage
rule forcing the particles to be closer to their mean:(10)
mti1  ati1  (1  a)t 1 with a  1  h 2
It can be verified that the mixture probability in
Eq. (9) has a covariance matrix Vt1 and that it does
not increase over time8). A dual state-parameter
updating scheme with kernel smoothing via the SIR
particle filter can be summarized in Fig. 2.

 t

1/ p

q sim (t ) Adown  s(t )   qdam (t  Tdam ) (12)3.6  k 
(13)
qt  qsim (t )   t

3. IMPLEMENTATION(1) Study area
The SMC filters were applied to the Katsura
River catchment (Fig. 3) to improve the river flow
forecasting. This catchment is located in Kyoto,Fig. 3 The Katsura River catchment.

I_3

(11)

where s is catchment storage (mm), t is time (hr),Adown is the downstream area from the dam (km2),
qsim is simulated river discharge (m3/s), qt and qdam
are observed discharge at the Katsura gauging
station and at the Hiyoshi dam (m3/s), TL and Tdam
are the lag time parameters of catchment andoutflow from the dam reservoir (hr), and k and p aremodel parameters. ωt and νt are the state and themeasurement error, respectively. Effective rainfall
re is estimated as follows:1
fraccum  Rsa
raccum  Rsa(14)

Discharge (m3/s) Discharge (m3/s) Storage S (mm) Rainfall (mm/h)

where r is rainfall (mm/hr), f is the runoff
coefficient, raccum is the accumulated rainfall amount(mm), and Rsa is the saturation amount (mm). Arealmean values of hourly observed rainfall from the 13
gauging stations were used as model input. Six
model parameters, including k, p, TL, Tdam, f and Rsa,
have been estimated from the events of 2004. In thestate only updating scheme, pre-calibrated
parameter values were used. On the other hand, thedual state-parameter updating scheme has been
performed on five model parameters, excluding
Tdam, which showed stable values compared toothers. Both simulations were performed by the SIR
particle filter with 3,000 particles. Covariance ofthe error of system (Wt) and measurement (Vt) wereassumed to be 4 mm and 10% of the current
observed discharge, respectively.

(4) Dual state-parameter updating scheme
In the dual state-parameter updating scheme,
initial values of each parameter have been set touniform distribution with widths that coverdeviations of pre-calibrated parameter distributions.
In other words, true static values of parameters are
assumed to be located within these initial
distributions. Inference of five parameters (e.g., k,p, TL, f and Rsa) was performed by the kernelsmoothing method in DUS. The value of kernelsmoothing parameter a in Eq. (9) was set as 0.95.

0
10

(a)

20
40

20

(b)
0

0
10

(a)

20
40

20

(b)

0

400

400

200

200

(c)
0

(c)
0

400

400
200
0
0

Discharge (m3/s) Discharge (m3/s) Storage S (mm) Rainfall (mm/h)

re  r  f  with f   (3) State only updating scheme
Fig. 4 shows the simulation results of the state
only updating via the SIR particle filter compared toobservations and a deterministic prediction. In thisscheme, particles are resampled in each observation
time step, and catchment storage (s) is fluctuatedaccording to the system noise.
While updated river discharge using a state only
updating scheme shows good conformity between
observation and simulation (Fig. 4(c)), a forecast
based on the same particles does not reproduce theriver flow properly compared to a deterministic
prediction (Fig. 4(d)). To compare off-line optimal
parameters with those calibrated from the past event(Table 1), several parameters show quite different
values. In this respect, it can be inferred that state
updating based on inappropriate parameters may be
one of the causes misleading the forecast.

200

(d)
50

100

150

Time (hour)

0
0

(d)
50

100

150

Time (hour)Fig. 4 Results of the state only updating via the SIR particle
filter from 11 to 16 July 2007. (a) Hourly precipitation.
(b) Catchment storage. (c) Updated river discharge. (d)3-hour-lead forecasted river discharge. Black dots
represent observed discharge. Blue line and area
represent mean value and 95% confidential interval,
respectively. Dashed line represents a deterministic
modeling case.

Fig. 5 Results of the dual state-parameter updating via the SIR
particle filter from 11 to 16 July 2007. (a) Hourly
precipitation. (b) Catchment storage. (c) Updated riverdischarge. (d) 3-hour-lead forecasted river discharge.
Black dots represent observed discharge. Blue line andarea represent mean value and 95% confidential interval,
respectively. Dashed line represents a deterministic
modeling case.

I_4

Table 1. Parameter information.

400

Model
Pre-calibrated
Off-line optimum
Initial range
parameters from 2004 events from 2007 events for dual updating

17.0
0.6
3.8
0.33
82.0
4.0

30.0
0.66
6.0
0.65
105.0
4.0

10.0~40.0
0.4~0.9
3.0~7.0
0.1~0.8
50.0~150.0
4.0

Simulation (m3/s)

k (-)
p (-)
TL (hr)f (-)
Rsa (mm)Tdam (hr)Dual updating
State only updating

Fig. 5 illustrates the simulation results of the dual
state-parameter updating. Compared with the state
only updating case, a forecast by the dual updating
scheme shows better conformity with observations(Fig. 5(d)). Furthermore, the unexpected drawdownof hydrograph in the rising part (Fig. 4(d)) is notshown in the dual updating case. Traces of thecatchment storage s present different patterns in
Fig. 4(b) and Fig. 5(b), whereas updated dischargehydrographs show similar traces in both cases.
Fig. 6 presents the traces of parameter
distribution. One can observe a significant
reduction of parameter uncertainty for allparameters after the first flood peak. In comparison
with the off-line optimum (Table 1), estimatedparameters show similar ranges, especially inparameter k, TL and f.

100

100

200

300

400

Observation (m3/s)

Fig. 7 Scatter diagram of simulation results. Cross dots
represent results of state only updating. Circle dots
represent results of dual state-parameter updating.
Table 2. Statistics on model accuracy.
State only
Deterministic
Dual updating
updating
RMSE (m3/s)

44.6

36.4

20.4

Nash
-Sutcliffe

0.73

0.82

0.94

It is worth noting that when the artificial
evolution is applied for parameter inference insteadof kernel smoothing in the dual updating scheme,
estimated parameters present more diffusive
distributions and unstable inference is produced
resulting in different posterior distributions at eachsimulation. However, inference from the kernelsmoothing presents relatively consistent results
because there is less uncertainty of parameters.
In the scatter diagram shown in Fig. 7, the dual
state-parameter updating scheme presents enhanced
simulation results in the overall flow regime fromhigh flow to low flow. Additionally, the model
accuracy criteria shown in Table 2 confirm that the
DUS is superior to other simulations.

40
30

k

200

0
0

50

20
10

1
0.8

p

300

0.6
0.4
0.2

8

TL

6
4

f

2
1

0.5

Rsa

0
150

100
50
0
0

50

100

150

Time(hour)Fig. 6 Traces of parameter k, P, TL, f, Rsa of the SF model using
dual state-parameter updating of the SIR particle filter
from 11 to 16 July 2007. Black lines represent medianvalue, and gray area represents 95% confidential
interval.

I_5(5) Comparison of various SMC filters
Several different versions of the SMC filters,such as SIR, ASIR, and RPF with the MCMC move
step, were implemented under the same simulation
conditions. The Markov chain Monte Carlo(MCMC) move step of RPF, which is used for
improving sample diversity in the resampling step,is based on the Metropolis-Hastings algorithm10).
The dual state-parameter updating scheme has been
adapted in all the cases with 3,000 particles. A
comparison of the simulated discharge hydrographis illustrated in Fig. 8. There is no significant
difference in the estimated 3-hour-ahead forecasting
via three SMC filters.
Although three SMC filters reproduce riverdischarge properly in the first flood peak (1~30hour) and the recession part, all the SMC methods

Discharge (m3/s)

scheme, parameter inference was performed by the
kernel smoothing method. A significant reductionof parameter uncertainty was observed for allparameters after the first flood peak, and estimatedparameter distributions showed good conformity
with off-line optimum. Performance results of SIR
and RPF showed similar forecasting accuracy,
while ASIR resulted in a slightly higher number of
errors than others. However, RMSE statistics ofthree SMC filters presented stable results when thenumber of particles was over 1,000.
The SMC filters are applicable to more complexhydrologic models, such as process-based andspatially distributed hydrologic models, in which it
is difficult to use the conventional data assimilation
methods. We will examine the performance of the
SMC filters on a distributed hydrologic model.

400
Obs.
SIR
ASIR
RPF

300
200
100
0
0

50

100

150

Time (hour)Fig. 8 Forecasted river discharge (3 hour ahead) by three SMC
filters from 11 to 16 July 2007. Black dots represent
observed discharge.
(a)

(b)

Fig. 9 Sensitivity analysis of the effects of particle numbers on
the prediction accuracy. (a) Updated river discharge. (b)Forecasted river discharge.

REFERENCES
1) Gordon, N. J., Salmond, D. J. and Smith, A. F. M., Novelapproach to nonlinear/non-Gaussian Bayesian state
estimation, Proc. Inst. Electr. Eng., Vol.140(2),pp.107-113, 1993.
2) Moradkhani, H., Hsu, K. L., Gupta, H. and Sorooshian, S.,
Uncertainty assessment of hydrologic model states andparameters: Sequential data assimilation using the particle
filter, Water Resour. Res., Vol.41, W05012, doi:10.1029/
2004WR003604, 2005.
3) Smith, P. J., Beven, K. J. and Tawn, J. A., Detection of
structural inadequacy in process-based hydrologicalmodels: A particle-filtering approach. Water Resour. Res.,
Vol.44, W01410, doi:10.1029/2006WR005205, 2008.
4) Salaman, P. and Feyen, L., Assessing parameter,
precipitation, and predictive uncertainty in a distributedhydrological model using sequential data assimilation withthe particle filter, J. Hydrol., Vol. 376, pp.428-442, 2009.
5) Arulampalam, M. S., Maskell, S., Gordon, N. and Clapp, T.,
A tutorial on particle filters for online nonlinear/
non-Gaussian Bayesian tracking, IEEE Trans. Signal
Proces., Vol.50(2), pp.174-188, 2002.
6) Ristic, B., Arulampalam, S. and Gordon, N.: Beyond theKalman filter: Particle filters for tracking applications,
Artech House, 2004.
7) Moradkhani, H., Sorooshian, S., Gupta, H. V. and Houser, P.
R., Dual state-parameter estimation of hydrological models
using ensemble Kalman filter, Adv. Water Resour., Vol.28,
pp.135-147, 2005.
8) Liu, J. and West, M., Combined parameter and state
estimation in simulation-based filtering, Sequential Monte
Carlo in Practice, Doucet, A., de Freitas, N. and Gordon,
N. eds., Springer-Verlag, pp.197-223, 2001.
9) Kimura, T.: The flood runoff analysis method by the storage
function model, The Public Works Research Institute,
Ministry of Construction, Japan, 1961.
10) Robert, C. P. and Casella, G.: Monte Carlo statisticalmethods, Springer, New York, 1999.

overestimate the discharge during 65-80 hours.
Uncertainty of forcing data (rainfall) and noconsideration of spatial heterogeneity in the SFmodel are plausible reasons.
Sensitivity analysis was performed concerning
the effects of particle numbers on the prediction
accuracy (Fig. 9). RMSE statistics of simulated
discharge show stabilized accuracy in both updating
and forecasting via three SMC filters when thenumber of particles exceeds 1,000. In terms offorecasting accuracy, SIR and RPF show similar
RMSE statistics, while ASIR presents a slightly
higher number of errors than others. Researchers
also stated in a previous study6) that if theimportance density of ASIR does not characterize
the transitional prior p( xt | xti1 ) for some reason(e.g., process noise is large), the use of ASIR caneven degrade the performance. The simulation timefor 1,000 particles is less than 2 min in three SMC
filters, which is short enough to be applied forreal-time forecasting.

4. CONCLUSIONS
Sequential Monte Carlo (SMC) filters wereapplied to a conceptual hydrologic model, thestorage function model, using state only updating
and the dual state-parameter updating scheme. The
river discharge forecast via the SMC filters wascompared with observations. The forecast provided
by the dual state-parameter updating scheme wassuperior to that of state only updating anddeterministic modeling in terms of the model
accuracy criteria, a scatter diagram, and simulatedhydrographs. In the dual state-parameter updating(Received September 30, 2010)

I_6