Embedding Kalman Filter into a Distributed Hydrological Model

Kim

Embedding Kalman filter into a Distributed Hydrological Model
Sunmin KIM1), Yasuto TACHIKAWA2), and Kaoru TAKARA2)
1) Graduate school of Urban and Environment Engineering, Kyoto University
Yoshida Honmachi, Sakyo-ku, Kyoto, 606-8501, JAPANE-mail: sunmin@rdp.dpri.kyoto-u.ac.jp
2) Disaster Prevention Research Institute, Kyoto University
Uji, Kyoto 611-0011, JAPAN
Abstract
This paper describes a methodology for embedding a Kalman filter into a distributed hydrologicalprediction system for real-time flood prediction. In most cases it is complicated to formulate a Kalmanfilter algorithm in the system structures of distributed models. However, the Kalman filter theory is
successfully coupled with a distributed hydrological model to update spatially distributed statevariables by using several techniques proposed here. To acquire the total water storages of a basinfrom discharge observations at the outlet, a Q-S curve is used as an observation equation. The ratio
method is introduced to update the distributed storage amount of a basin, maintaining the spatiallydistributed water storage pattern of the basin by multiplying the distributed state variables in the modelby a specified ratio. A Monte Carlo simulation is adopted to predict state variables and error
covariance propagations. The methodology for deciding system noise is also discussed. A distributedmodel coupled with the Kalman filter theory gives updated simulation results with improved forecasting
accuracy.

1. INTRODUCTION
To obtain an accurate prediction in a real-time rainfall-runoff simulation it is essential to set model
parameters effectively and to give properly assessed initial state variables. Modern real-timeobservation systems make it possible to access various hydrological data while real-time simulation iscarried out. Prediction accuracy should be improved when the parameters or state variables are
updated with real time observed data. Also, the uncertainty caused by improper model parameters,
initial state variables and input data can be reduced if filtering theory, such as the Kalman filter, isincorporated in a hydrological model.
R.E. Kalman (1960) published his famous paper describing a recursive solution, which was later
named as the Kalman filter, to discrete data linear filtering problems. Having potential for broader use,
Kalman filter has been enhanced as Extended Kalman filter for nonlinear systems. The Kalman filter isan optimal recursive data processing algorithm to estimate the state variables for minimizing the error
statistically. It combines all available observation data, plus prior knowledge about the system andmeasuring devices, to produce an estimate of the desired variables in such a manner that the error isminimized statistically (Maybeck, 1979). A more detailed description about Kalman filter theory can befound at Jazwinski (1970) and a good discussion of the filter with several application cases to thehydrological system is given by Bras and Rodriguez-Itulbe (1985).
Since Hino (1974) initially adapted the Kalman filter theory to a hydrological system, numerous studies
have been carried out to use the filter theory in the field of hydrology. Takaso et al. (1989) describes
real-time flood forecasting based on a stochastic state-space formulation of rainfall-runoff systems
coupled with the Kalman filtering-prediction theory and its application. In the research, the storage
function method was used to couple with the filter taking the storage amounts of sub-basins as thestate vector. Lee and Singh (1999) showed upgraded simulation results of the tank model when thestate vector of the Kalman filter is composed of the model’s parameters. The storage function method
and the tank model are often used lumped models in Japan, Korea and many other countries for flood
forecasting and watershed modeling. While the Kalman filter has been applied to many lumped
models for better simulation or more accurate forecasting, it has hardly ever been applied to
distributed hydrological models. One of the main reasons is that unlike lumped models, in most cases
it is complicated to formulate the Kalman filter algorithm in the system structures of distributed models.

WSUD2004 Kim, Tachikawa & Takara

278

Embedding Kalman Filter into a Distributed Hydrological Model

Kim

A large number of state variables based on a fine grid cell hydrologic system also makes it harder toapply the Kalman filter.
In this research, to avoid the computational burden for updating each state variable, several
techniques are introduced for applying the Kalman filter to a distributed hydrological model. The Q-Scurve which is determined under steady state assumption on a study basin is used for the observation
equation to update the simulated total storage amount of the basin with discharge observations. To
consider the spatial pattern of the updated storage amount in every grid-cell of a distributedhydrological system, the ratio method is adopted. The ratio method is used for updating spatiallydistributed storage amount in the model by multiplying by a ratio calculated from the updated total
storage amount and the simulated storage amount. For the prediction algorithm, stochastic analysis isadopted to predict state variables and error covariance of the next updating step. Monte Carlosimulation is an effective technique to analyze the effect of error covariance propagation. The
methodology for deciding system error variance is also discussed.

2. COUPLING OF CDRMV3 WITH KALMAN FILTER
The study mainly focuses on a coupling method of the Kalman filter to the Cell based Distributed Runoff
Model Version 3 (CDRMV3, http://fmd.dpri.kyoto-u.ac.jp/~flood/product/cellModel/cellModel.html).
The state variable to be updated is the total storage amount and its spatial distribution is calculatedusing water depth at all computational nodes in the model. The water depths are easily converted tostorage amount by multiplying the cell area. The parameters of the CDRMV3 are calibrated before
applying the Kalman filter and do not change when state variables are updated. Uncertainties causedby systems and observations are considered in the error covariance of the filter, though uncertainty
caused by rainfall forecasting is not accounted for in this study. Radar observed rainfall data which iscalibrated by ground gauges are used as forecast rainfall data.

2.1

Brief Model Description of CDRMV3

CDRMV3 is a one dimensional physically based distributed hydrologic model developed at Flood
Disaster Research Laboratory of Disaster Prevention Research Institute, Kyoto University. The model
solves the Kinematic wave equation using Lax Wendroff scheme on every node in a cell (Kojima et al.,2003). Discharge and water depth propagate to the next cell according to a predefined routine order
determined in accordance with DEM data. An advantage of the CDRMV3 is that the stage-discharge
relationship of each cell reflects the topographic and physical characteristics of its own cell. Specified
stage-discharge relationship, which incorporates saturated-unsaturated flow mechanism, is included ineach cell (Tachikawa et al., 2004). Because of the variable slope and roughness coefficient, each cell
has its own relationship:vC d C (h / d C )β , (0 ≤ h < d C )q (h) = vC d C + v a (h − d C ), (d C ≤ h < d S )mvC d C + v a (h − d C ) + α (h − d S ) , (d S ≤ h)
where

(1)

vC = k C i , v a = k a i , α = i / n .
kc : hydraulic conductivity of unsaturated layer
ka : hydraulic conductivity of saturated layern : roughness coefficienti : slope of grid

The stage-discharge relationship is expressed by three equations corresponding to the water levels
divided into three layers (see Equation 1). When the water depth h is lower than the depth ofunsaturated layer (0≤ h<dc), flow is described by Darcy’s law with a degree of saturation, (h/dc)β, and ahydraulic conductivity kc. If the h increases, flow from the saturated layer is considered with a different
hydraulic conductivity ka of saturated layer. After the water depth is greater than the soil layer ds,
overland flow is added by using the Manning’s equation. According to this mechanism, the equations

WSUD2004 Kim, Tachikawa & Takara

279

Embedding Kalman Filter into a Distributed Hydrological Model

Kim

between discharge per unit width q and water depth h are formulated. Figure 1 shows the graphicalrelationship between q and h. More detail on the specified state-discharge relationship and the model
structure can be found in Tachikawa et al. (2004).
αβFigure 1. Relationship between unit width discharge and water depth in the CDRMV3
The model is applied to the Kamishiiba basin (211km2) of Kyushu area. Four different flood types ofthe basin are selected for this study. Flood period and maximum discharge of the four events areshown in Table 1.
Table 1. Flood events used in the study
EVENT
Event 979
Event 996
Event 998
Event 999

2.2

Flood Term
Sep 15~19, 1997
Jun 24~Jul 3, 1999
Aug 1~7, 1999
Sep 22~27, 1999

Max Discharge
1203.0 m3/s
210.0 m3/s
489.0 m3/s
644.0 m3/s

Measurement Update Algorithm

In the measurement update algorithm of the Kalman filter, an observation equation which specifies a
relation between observed data and state values is necessary. The observation vector yk can be
described as a linear combination of a state vector xk as Equation 2. The observations are affected bywhite noise wk which has a covariance matrix Rk. The m×n matrix H relates the state vector to theobservation. In the measurement update algorithm, the state vector x(klk-1) and error covarianceP(klk-1) which are estimated through system at time step k-1, are updated at time step k with use of
the n×m matrix K. The matrix K which is called Kalman gain is chosen to minimize the updated error
covariance P(klk). In the algorithm, ^ indicates estimated value and T indicates the transpose of amatrix.
Observation equation

yk = H k xk + wk ; wk ~ N (0, Rk )(2)

Measurement update algorithmxˆ (k k ) = xˆ (k k − 1) + K k ( y k − H k xˆ (k k − 1))(3)

P (k k ) = P (k k − 1) − K k H k P(k k − 1)
K k = P(k k − 1) H ( H k P(k k − 1) H + Rk )
T
k

T
k

(4)
−1

(5)

The difference, yk-Hkx(k│k-1), which is called the residual or innovation reflects the discrepancy
between the estimated observation Hkx(k│k-1) and the actual observation yk. If the total storage
amount to be updated is measured directly, the residuals are easily obtained. However observed
quantities are discharge or river stage rather than distributed storage amount. When we think of thestorage amount as a state variable, we can easily come up with a nonlinear relationship between
discharge and total water storage in a subject basin as shown in Equation 6. In the storage function

WSUD2004 Kim, Tachikawa & Takara

280

Embedding Kalman Filter into a Distributed Hydrological Model

Kim

method, when it is coupled with the Kalman filter, the observation equation is adopted after
linearization of Equation 6.
(6)
S (t ) = KQ (t ) P
where S(t) : total water storage in a basin
Q(t) : discharge at the outlet
K, P : constants.
Unlike the storage function method the CDRMV3 does not have general relation between observed
discharge and the storage amount. The relationship between the discharge at the outlet and the total
amount of storage represents a loop shape as shown in Figure 2, whose shapes are different for eachflood event. However, it is still possible to get a relation in a specific case like a steady state condition.
After reaching the steady state condition with a given constant rainfall on the subject basin, the total
storage amount that corresponds to the given rainfall intensity can be acquired by multiplying cell areawith water depths of every cell and sum up these entire amounts. The cell size in this study is
250m×250m. Applying various rainfall intensities, the Q-S curve that is used for the observation
equation is obtained as shown in Figure 2. Even though states during a runoff simulation are notsteady, the difference between the two curves from the steady state and unsteady state seems
acceptable.
The reason for needing the observation equation is to get the matrix H in the measurement update
algorithm. The conversion matrix H relates total storage amount and observed discharge at time step k.
More specifically, it stands for the gradient of the Q-S curve in accordance with simulation results at
updating time step. As only one observation is available in the Kamishiiba basin, H is a scalar value inthis study.
Q-S curve
Event 979

runoff simulation(unsteady state)
two different Hs(steady state)

Figure 2. Two different curves between steady and unsteady state
In processing the measurement update algorithm, a couple of problems exist to be considered. These
are basically caused by the steady state assumption when making use of the Q-S curve. Because therelationship is determined under the steady state assumption, there is always some departure towards
the unsteady state, even if it is small.
At first, two different H values could be given at each time step as seen in the Figure 2. At each time
step, discharge and total storage amount are obtained from simulation results. Those two values fromunsteady state simulation may not match on the Q-S curve. This mismatch often gives two different H;
one is from the simulated discharge and the other is from the simulated storage amount. However,through several tests, it is checked that those two different H values do not make a recognizable
difference to the filtered results. For this reason, an arithmetic average of those two H values is usedin the application of the CDRMV3.
Another problem occurs while getting residuals in the measurement update algorithm. According to theconventional equations of the Kalman filter, the residuals are calculated by the use of the observation
equation. As shown in Equation 7, converting state variables needs the H value which contains thesteady state error described already. On the other hand, residuals can be directly calculated from thedifference between observed discharge and simulated discharge as in Equation 8. Filtered results

WSUD2004 Kim, Tachikawa & Takara

281

Embedding Kalman Filter into a Distributed Hydrological Model

Kim

from those two residuals show significant difference. In conclusion, the residuals from directly using
simulated discharge gives much better filtered results. More detail discussion is provided in Section 3.
(7)( yk − H k xˆ (k k − 1)) ⇒ (OQk − H k × SS k + Ak )( yk − H k xˆ (k k − 1)) ⇒ (OQk − SQk )(8)

where OQk : observed discharge
SQk : simulated result of discharge
SSk : simulated result of total storage amount
Hk: conversion matrix, gradient of the Q-S curve
Ak: optional input for the linearization
After updating the total storage amount through the measurement update algorithm, the updated
storage amount should be distributed to each cell in a subject basin. One efficient way to update eachcell’s storage amount is using a specific ratio calculated from the updated total storage amount andthe simulated storage amount. The calculated ratio is applied to all water depths of each cell in themodel, which has the same spatial distribution pattern with the simulation result before updating as
shown in Figure 3. For example, if the simulated storage amount at a specific time step is 1.03476E+8
m3 and the updated value is 0.97292E+8 m3, all water depths of each cell at this time are multiplied by
the ration 0.9402, and then the simulation starts again using each cell’s updated state variable. Thismethod which is named as the ratio method, offers efficient and effective updating skill of statevariables considering its spatial distribution pattern (Kim et al., 2004).

State variablesof every cell

Updated storage amount
Simulated result one

Figure 3. Concept of the ratio method for updating spatially distributed state variables

2.3

Time Update Algorithm

The Kalman filter is an algorithm to optimize the state vector x of a discrete time controlled process
which is governed by a linear difference equation. The n×n matrix F in the system equation relates thestate at previous time step k to the state at current step k+1. The systems are continuously affected bywhite noise, vk, with covariance matrix Qk respectively. The matrix Bk relates optional control input tothe state x. The time update algorithm is for projecting forward the current state and error covarianceto obtain the estimation for the next time step. The estimated error covariance P is a n×n matrix.
System equation

xk +1 = Fk xk + Bk + vk ; vk ~ N (0, Qk )(9)xˆ (k + 1 k ) = Fk xˆ (k k ) + Bk

(10)

P(k + 1 k ) = Fk P(k k ) F + Qk

(11)Time update algorithmTk

In the CDRMV3, a complicated relation exists between the present state variables and the next statevariables; in this case the storage amount and its spatial distribution. Each cell responds
interdependently to the next step’s state variable according to its present state variable and other input
data such as rainfall. So it is impractical to formulate the system matrix Fk, which is essential to update
the error covariance P(k+1│k), even though it is possible theoretically. Rather than use the conventional
concept of the Kalman filter theory as shown in the schematic drawing (a) in Figure 4, the Monte Carlosimulation technique (drawing (b) in Figure 4) is useful to solve this kind of problem.

WSUD2004 Kim, Tachikawa & Takara

282

Embedding Kalman Filter into a Distributed Hydrological Model

Kim

By the concept of the Monte Carlo simulation, many sets of random variables are generated at time
step k; hundreds of storage values were generated in this study. Each random variable is defined byany possible storage amount within a range of probability distribution, N(Xk, σk); where, σk=Pk0.5.

Qkxˆ (k k )Hydrologic

P(k k )SYSTEM

Fk xˆ ( k k ) + Bk

Fk P ( k k ) FkT + Qk(a) Estimation using Kalman filter

P(k + 1 k ) = (σ k +1 ) 2 + Qk

Hydrologic

X k, σ k

SYSTEM

Xk+1, σ k + 1(b) Estimation using Monte Carlo simulation(σ k +1 ) 2 ≈ Fk P(k k ) FkT(σ k ) 2 = P(k k )Figure 4. Concepts of time update algorithm
The ratio method is used again at this point to reset the water stage at each cell by multiplying by theratio of generated storage amount to the mean storage amount Xk. After a simulation repeatedly
calculates multiple input sets until the next update time step k+1, another probability distribution,
N(Xk+1,σk+1), is calculated from the simulated results. Now, the estimated state Xk+1 is the mean valueof the probability distribution and the estimation error covariance FkPkFkT is regarded as (σk+1)2. Adding
the system error covariance Qk completes the error covariance Pk+1 at time step k+1. Estimation error
covariance FkPkFkT means propagation of the error covariance Pk through the simulation, and Qk
stands for a new generated or added system covariance during simulation from time step k to the nexttime step k+1. The newly added system covariance is caused by system structure or new input datasuch as rainfall. The methodology to determine the system error covariance, Qk, is discussed in thefollowing section.

Figure 5. Propagations of probability density of storage amount and discharge (Event 979)
Figure 5 shows an example of the probability distribution propagation during three hours from Monte
Carlo simulation results. The three lines on the bottom represent the variation of mean value andstandard deviation from the distribution at each calculation time step. While the distributions of storage
amount follow the normal distribution pattern, the distributions of discharge which are composed byeach simulation results set from the storage amount variables do not always stick to the normal
distribution. This phenomenon is caused by the nonlinear relationship between storage and discharge.
The different distribution pattern of discharge and storage amount will also affect the form of the error.

2.4

Setting the Observation Noise and System NoiseOne of the difficulties in applying the Kalman filter is determining the error covariance of the systemand the observation. Although the Kalman filter provides an algorithm for better forecasting by
updating the state estimates, its success depends on determination of the error covariance which

WSUD2004 Kim, Tachikawa & Takara

283

Embedding Kalman Filter into a Distributed Hydrological Model

Kim

requires proper judgment by the hydrologist. Because the hydrological system is a natural system
which varies in time and space, it is impossible to get the true value. This means that it is also
impossible to get the error covariance which is based on noises to the true value of hydrologicalvariables. For this reason, error covariance has been assumed when the Kalman filter is applied tohydrological models. The only action we can take is to try to get a reasonable error covariance with theleast assumptions.

2.4.1 Observation NoiseThe basic assumption of the Kalman filter is that system and observation noises are white andGaussian. Thus, at any point in time, the noise value is not correlated in time and the probability
density curve of noise takes on the shape of a normal bell-shape. This assumption can be justified
physically by the fact that a system or observation noise is typically caused by a number of small
sources (Maybeck, 1979). It is reasonable to see the observation noise from this point of view. Usuallyobserved values are acquired through a conversion from stage to discharge by using a stagedischarge relation or rating curve. Observation noises are mainly caused by misreading of gauge
stage, interpolation during conversion from stage to discharge, and improper stage-discharge
relationship. Also it is not difficult to accept that the noise is normally distributed to the true value. The
observed data used in this study are acquired from the Kamishiiba dam inflow data. The inflow data
are calculated mainly by converting the dam reservoir stage into discharge with considering therelease from the dam for various purposes. Since more research is needed to determine theobservation noise reasonably, an assumed observation error variance is used in this paper. The
assumed error variance is mentioned in each filtered result.

2.4.2 System NoiseIf observed data are assumed to be true values, differences between simulation results and theobserved data could be regarded as system noise. First, it needs to be examined whether thevariances of simulation results to the observed data have a normal distribution. When the offline
simulation results from the CDRMV3 are examined, the variances have a normal probability
distribution to the observed data as sown in the Figure 6. Under this examination, the RMSE given by
Equation 12 could be regarded as a standard deviation of system noise wk. The mean value and
RMSE from the offline simulation results of four different events are shown at Table 2. We can figureout that the mean values of variances are around zero and the RMSEs are around 30m3/s. Following
this analysis, the standard deviation of system noise in terms of discharge is set as 30m3/s in thisstudy.

(a) Offline simulation results and observed data
65.9
50
34.1
18.2
2.3
-13.6
-29.4
-45.3
-61.2
-77.1
0

(b) Variances of simulation results to the observed data

50

100

150

200

250

300

350

Frequency(c) Probability distribution of the variances

Figure 6. Variances of system error and its probability distribution (Event 998)

WSUD2004 Kim, Tachikawa & Takara

284

Embedding Kalman Filter into a Distributed Hydrological Modeln

RMSE =∑ (Q − Q )
i =1

S ,i

Kim

2

O ,i

(12)

n

where RMSE : Root Mean Square Error
QS : simulated discharge
QO : observed dischargen : number of QS or QO values
Table 2. The RMSE and mean values of each event
EVENT
Event 979
Event 996
Event 998
Event 999

RMSE
33.64 m3/s
21.15 m3/s
20.85 m3/s
23.73 m3/s

MEAN
-2.22 m3/s
-10.80 m3/s
3.32 m3/s-1.02 m3/s

Then, the discharge RMSE is converted to the error variance of the total storage amount. System
noise in terms of discharge can be translated to the noise in terms of storage amount by using the Q-Scurve as seen in the Figure 7. The term ‘sd’ in the figure means the standard deviation of discharge
variances, which was previously determined to be 30m3/s. Three discharges from simulation results at
a specific time step will match with three different storage amounts on the Q-S curve. Because of
nonlinearity of the relation, the differences between the upper value and lower value to the Sk will bedifferent. Using those two differences of storage amount, Supk–Sk and Sk-Sdnk, the system error
variance Qk can be calculated as shown in Equation 13.

Qk = ( Supk − S k ) × ( S k − Sdnk )(13)

where Qk : system error variance at time step kSup, S, Sdn : converted storage amounts

SQk+sd
SQk
SQk-sd

Sdnk
Sk
Supk

Figure 7. Conversion of noise from discharge to storage amount form
There is one important checking point about the Gaussian assumption of the Kalman filter for anonlinear system. Because the relationship between discharge and storage amount is nonlinear, if theprobability distribution of storage amount is a Gaussian distribution, the distribution of discharge will
not follow the normal distribution, and vice versa. The distributions of the variables are no longer
normal after undergoing their respective nonlinear transformations. However, this nonlinear effect on
probability distribution is not significant in this study. As shown in Figure 5, the probability distributionof discharges, which is calculated from the simulation results according to each storage amounthaving normal distribution, can be roughly regarded as a normal distribution.

WSUD2004 Kim, Tachikawa & Takara

285

Embedding Kalman Filter into a Distributed Hydrological Model

Kim

3. ANALYSIS OF RESULTS
3.1

Results from two different residual types

The Kalman filter, a recursive data processing algorithm, is successfully coupled to the distributedhydrological model, CDRMV3. To check the filtering results effectively, the observed data are
assumed to have no error. If there is no observation error, which means that the observed data is thetrue value, the filtered results and observed data should match exactly. When the results are checked,
there are discrepancies at the beginning and after the peak of the hydrograph compared to the otherpart as seen in Figure 8.
Three reasons are considered to explain these discrepancies. The first one is because of the steadystate assumption of the Q-S curve. When the steady state assumption is made, discharge and storageis expressed as a single-valued function (see Figure 9 (a)). On the other hand, different storage values
occurs in unsteady state conditions, even though the discharge is the same as shown in Figure 9 (b).
The storage amount at the beginning of the hydrograph is different from the amount at the falling limb
or after the peak of the hydrograph. The differences in storage amount, “Storage B-A” at the beginningof the runoff and “Storage C-A” after the peak, cause the under estimation of discharge at the risinglimb and the over estimation at the falling limb of the hydrograph.

(a) Event 979(b) Event 996

Figure 8. Updated results comparisons from two different residuals
at the beginning of runoff

at the ending of runoff

Storage B

Storage A

Discharge Q(a) Steady state

Storage C

Discharge Q

Discharge Q(b) Unsteady state

Figure 9. Conceptual storage amount distributions according to the state differences
Secondly, every observation update is given at every an hour while the calculation time step is ten
minute. Five simulation results between nearest observation update would make its own hydrograph.
When observation updates are carried at every calculation time step, the discrepancies are decreased
a lot. The last reason is that the filtered results presented as the discharge hydrograph are one
calculation time step ahead prediction results. The updated storage amounts are used to update eachcells water depth, and then the next step’s discharges are calculated using the updated water depthsof the previous step. Even thought it is small, it results in some discrepancies in the hydrograph.
Two different updated results are presented in Figure 8; the hydrograph labeled “using equation 7” isthe result when Equation 7 is used to calculate the residuals in the measurement update algorithm,
and the other one labeled “using equation 8” is using Equation 8 in the algorithm. Those two kinds ofresults are obtained under the setting that the observation error variance is zero and the standard

WSUD2004 Kim, Tachikawa & Takara

286

Embedding Kalman Filter into a Distributed Hydrological Model

Kim

deviation of system noise is 30m3/s. Direct use of the simulated discharge to obtain residuals gives
much better filtered results than using the H value from the Q-S curve. When these two residuals at
each updating step are compared, around 10 to 30m3/s of difference is observed. In the case of “using
Equation 8”, the discrepancies to the observed data are not recognizable as much as in the “using
Equation 7” case. Every result from now on is obtained using Equation 8.

3.2

Comparison the updated results from different error variance

Figure 10 shows the various updated results by setting a different error variance. The label “oe30”
means 30m3/s of standard deviation as an observation noise and “se30” means 30m3/s of standard
deviation as a system noise, “se0” means no system noise, and “oe0” stands for no observation noise
which means the observed data are regarded as true values.
When the hydrographs from the case “se0.0:oe30” are examined, the filtered results are exactly thesame as the results from offline simulation because the system is regarded as a perfect one to do asimulation. On the other hand, the case “se30:oe0.0” shows that the filtered results trace the observed
data. The case “se30:oe30” show that the filtered results are located between the offline simulation
results and the observed data. Even though the noises are set to the same value as 30m3/s, thefiltered results are closer to the observed data than offline simulation results. Several reasons are
considered to explain this phenomenon such as nonlinear observation equation, different form of error
variance and initial error variance. While the observation error variance is in terms of discharge, thesystem error variance is in terms of storage amount which is transformed from the discharge noise.
The initial error variance of the Kalman filter also affects the filtered results; the initial error variance is
set as 20m3/s of discharge noise in this study. Further research is needed to determine how thesefactors affect the filtered results.

(a) Event 979(b) Event 996

Figure 10. Updated results from variant error variance

3.3

Another noise form: Coefficient of Variation

It is sometimes unreasonable to set the standard deviation of noise as a constant value, no matter
how much discharge there is. When discharge is comparably small, for example 50 or 100m3/s,
30m3/s of noise is a large value to set as a standard deviation. On the other hand, if discharge is over
1000m3/s, that amount of the noise could be negligible. For this reason, it is more reasonable to set
the noise not as a deterministic value but as a ratio to discharge, such as a form of Coefficient of
Variation (C.V.) which is calculated as a ratio of standard deviation to mean value. The 30m3/s ofdeterministic system noise can be converted to the C.V. by Equation 14. The converted C.V. value
from the four flood events gives a 10% value for noise. Every result here after comes from the C.V.
error form. The results shown in Figure 11 are obtained when both observation noise and system
noise are set at 10% of observed discharge and simulated discharge. Figure 12 shows the variation of
the standard deviation of discharge results using two different noise forms. It can be easily recognizedthat the variations from C.V. form show reasonable results in accordance with discharge variations.

WSUD2004 Kim, Tachikawa & Takara

287

Embedding Kalman Filter into a Distributed Hydrological ModelKim∑ ( y × cv) = ∑ ( x − y )
2

i

i

n

2

i

n

(14)

where cv : coefficient of variation
yi : simulated discharge
xi : observed dischargen : number of discharge values(a) Event 979(b) Event 996

Figure 11 Filtered results when the noises are C.V. form

Figure 12. Standard deviation of discharge

3.4

Prediction Efficiency

To check the prediction accuracy after coupling with the Kalman filter, 1hr, 6hr and 12hr prediction
results are compared. Figure 13 shows the prediction results when the system noise and observationnoise are equally set as 10% of C.V. Table 3 shows the RMSE from the prediction results. The RMSE
is calculated by Equation 12. As expected, prediction for short lead times shows higher accuracy. It isinteresting that even prediction of 12hr ahead gives quite good accuracy compared to the short leadtime forecasting. One main reason for this is the use of recorded rainfall data. Simulation andprediction are carried out under an assumption that we know exactly the expected rainfall. Because
this paper focuses on coupling the Kalman filter to a distributed hydrological model, it needs todecrease the factors which impact on the filtered results. If the uncertainty of the rainfall forecast isprovided, it should be added in the system error variance Qk.

(a) Event 979(b) Event 996

Figure 13. Prediction results when the C.V. is 10%

WSUD2004 Kim, Tachikawa & Takara

288

Embedding Kalman Filter into a Distributed Hydrological Model

Kim

Table 3. The RMSE of predicted results (unit: m3/s)
EVENT
Event 979
Event 996
Event 998
Event 999

1 hr predict
37.18
11.42
16.45
28.00

6 hr predict
39.64
17.24
21.93
34.72

12 hr predict
37.28
20.14
22.17
27.85

4. CONCLUSION
The Kalman filter was successfully coupled with the distributed hydrological model, CDRMV3, to
update the state variables. Rather than formulate an impractical algorithm to apply the filter, several
techniques, such as the Q-S curve, ratio method and Monte Carlo simulation are used. Total storage
amount from the Q-S curve is used as an observed storage amount. The ratio method is used for
updating each water stage of every cell in the model by multiplying by the ratio calculated fromupdated total storage amount and the simulated storage amount. For the prediction algorithm, Monte
Carlo simulation is adopted to predict state variables and error covariance at the next step. Monte
Carlo simulation is an effective technique to analyze the propagation of error covariance.
The CDRMV3 using Kalman filter yields better results than the CDRMV3 without the filter in terms of
RMSE and computed hydrographs. Prediction efficiency after coupling with the Kalman filter, for 1hr,
6hr and 12hr prediction results, shows quite good accuracy when compared with the observed datawhen prediction are carried out under an assumption that the expected rainfall is exactly known.
Research to overcome the steady state assumption on the Q-S curve is needed to improve the filteredresults. A methodology to include the uncertainty of rainfall forecasting in the system error covarianceis a further research issue.

5. ACKNOWLEDGEMENTS
A part of this study is supported by "System Modelling Approaches for Assessment of Interaction
Between Social Changes and Water Cycle (PI: Prof. Kaoru Takara, Disaster Prevention ResearchInstitute, Kyoto University)", which is conducted under the CREST (Core Research for Evolutional
Science and Technology) Program by the Japan Science and Technology Agency.

6. REFERENCES
Bras, R. L. and Rodriguex-Itulbe, I. (1985) Random function and hydrology. Chapter 8 and 9, AddisonWesley, Reading, Mass., 425-520.
Hino, M. (1974) Kalman fiter for the prediction of hydrologic runoff system - application of theory, Journalof JSCE, Volume 221, 39-47.
Jazwinski, A. H. (1970) Stochastic Processes and Filtering Theory. Academic Press, New York.
Kalman, R. E. (1960) A New Approach to Linear Filtering and Prediction Problems, Transaction of the
ASME – Journal of Basic Engineering, 83, 35-45.
Kim, S., Tachikawa, Y., and Takara, K., (2004) Real-time Updating of State Variables in a Distributed
Hydrological Model, Annuals of Disaster Prevention Research Institute, Kyoto Univ., No. 47B.
Kojima, T. and Takara, K. (2003) A grid-cell based distributed flood runoff model and its performance,
Weather radar information and distributed hydrological modeling (Proceedings of symposium HS03 held
during IUG2003 at Sapporo, July 2003), IAHS Publ. No. 282, 234-240.
Lee, Y. H. and Singh, V. P. (1999) Tank model using Kalman filter, Journal of Hydrologic Engineering,
ASCE, Volume 4, No. 4, 344-349.
Maybeck, Peter S. (1979) Stochastic Models, Estimation, and Control, Volume 1, Academic Press, NewYork.

WSUD2004 Kim, Tachikawa & Takara

289

Embedding Kalman Filter into a Distributed Hydrological Model

Kim

Tachikawa, Y., Nagatani, G., and Takara, K. (2004) Development of stage-discharge relationshipequation incorporating saturated- unsaturated flow mechanism, Annual journal of hydraulic engineering
Vol. 48, 7-12.
Takasao, T., Shiiba, M., and Takara, K., (1989) Stochastic state- space techniques for flood runoff
forecasting, Pacific international seminar on water resources systems, Tomamu, 117-132.

WSUD2004 Kim, Tachikawa & Takara

290