土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10

バイアス補正カルマンフィルタによる
広域分布型流出予測システムのデータ同化
佐山敬洋1 ・立川康人2 ・寶 馨3
1 正会員 京都大学防災研究所助教 社会防災研究部門 (〒 611-0011 宇治市五ヶ庄)

E-mail: sayama@flood.dpri.kyoto-u.ac.jp
2 正会員 京都大学大学院准教授 工学研究科都市環境工学専攻 (〒 615-8540 京都市西京区京都大学桂)
3 フェロー会員 京都大学防災研究所教授 社会防災研究部門 (〒 611-0011 宇治市五ヶ庄)

広域分布型流出予測システムの観測流量によるデータ同化手法を提案する．予測システムは斜面部の流れを表現する流出モデルと河道流れを表現する河道追跡モデルからなる．これらのモデルが持つ全ての状態量を実時間で観測更新することは計算付加が高く，実時間予測システムとしての実行可能性に困難が伴う．そこで，本研究では河道追跡モデルにマスキンガム−クンジ法を用い，河川流量を観測更新すると共に，流出モデルに起因する予測のバイアスを河川流量と同時に逐次推定する方法を提案する．この手法を桂川流域の洪水予測に適用し，斜面部の流出予測バイアスを補正することによって洪水予測精度が向上することを明らかにした．

Key Words : distributed rainfall-runoff model, Muskingum-Cunge Method, flood prediction, data assimilation

1.

を予測している6) ．モデルの初期状態は，実績降雨をはじめに入力した流出シミュレーションによって 1 時間毎に更計画規模に匹敵するような豪雨は，治水整備の不十新している．ただし，河川流量の観測情報は国土交通

分な中小河川流域を中心に，甚大な洪水被害をもたら省の光ネットワークを通じてリアルタイムで入手してす．例えば，平成 16 年の新潟・福島豪雨では，五十嵐いるものの，現行のシステムでは，その観測情報を状

川や刈谷田川など 11 箇所で破堤し，16 名の死者・行方態量の推定には用いていない．今後，このシステムの不明者をもたらした．また，平成 17 年の台風 14 号で予測結果を向上させるためには，流域内の複数地点では，九州地方を中心に総雨量が 1000 mm を超え，五ヶ観測される河川流量を用いて広域分布型流出予測シス

瀬川や大淀川で大規模な越水氾濫が発生した．宮崎県

テムの状態量を同化するアルゴリズムを導入する必要

ではこの災害により 4300 棟を越える住宅が床上・床下がある．

浸水の被害を受けた1) ．こうした中小河川流域では，治流出モデルを観測流量で同化する方法については，

水施設の整備のみに依存した防災は今後も容易ではな

1974 年に日野7) がはじめてカルマンフィルタ8) を流出
モデルに適用して以来，数多くの手法が提案されてきた9),10),11),12) ．例えば，高棹・椎葉13) は単一斜面のキ
ネマティックウェーブモデルを準線形化してカルマン
フィルタを適用し，モデルの状態量を逐次推定する方法を提案している．また，高棹ら14),15),16) の一連の研究は，貯留関数法にカルマンフィルタを導入して実流域に適用するとともに，有色ノイズの応用や，複合流域への応用を示している．また，近年になって，流出予測の不確実性を評価するという観点から，異なるパラメタを持つ複数のフィルタを同時進行させてフィル
タの重みを逐次更新する方法17) やアンサンブルカルマ
ンフィルタ18) を比較的構造の単純な流出モデルに適用
して，パラメタ分布と予測の不確実性を逐次推定する
研究19),20),21),22) が提案されている．

く，詳細で確度の高い洪水予測と，それに基づく的確な住民避難が災害を軽減するための重要な対策となる．
これまでの洪水予測は，流域下流の主要な地点やダ
ム地点を対象に，主に貯留関数法などの集中型モデルを用いて行われてきた2),3),4) ．こうした予測手法は，モデルを同定するための流量が，予測対象地点で得られる場合には有効である．しかし，中小河川を含めて観測流量が得られない流域内部の任意の河道地点を対象に予測するためには，分布型流出モデルを流域全体に適用し，流域内部の河川流量を予測する必要がある．
筆 者 ら は ，流 域 全 体 の 詳 細 な 洪 水 予 測 を 目 指し て ，淀 川 流 域 を 対 象 と し た 広 域 分 布 型 流 出 予
測 シ ス テ ム を 開 発 し て き た5) ．現 行 の シ ス テ ム(http://yodogawa.kuciv.kyoto-u.ac.jp/) で は ，
予測雨量を時々刻々入力し，6 時間先までの河川流量

226

分布型流出モデルのフィルタリングについては，い土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10

まだ確立した方法は無いものの，近年になっていくつか

ンフィルタで計算する分散共分散行列やカルマンゲイの研究成果が報告されている．流域斜面系のフィルタンの次元が二倍に拡大し，計算負荷の大きい非効率なリングとしては，例えば，キネマティックウェーブ型の解法となる．

分布型流出モデルにカルマンフィルタを導入した Kim

デルに変分法のデータ同化を適用した Seo ら26) の研究

Friedland30) は適応カルマンフィルタの計算次元が拡大する問題に対し，行列計算の次元を状態量の数に抑えつつ，未知定数と状態量とを順に求める方法を提案している．また，これによって得られる推定量が，適応カルマンフィルタによる推定量と一致することを理論的に証明している．さらに，Dee and Da Silva31) は，Friedland30) によって提案された方法を拡張し，予測に含まれるバイアスを推定・補正するより実用的な方法を展開している．

がある．ただし，この研究で使用しているモデルの状

そこで本論はこれらの手法を参考とし，広域分布型

態量は全部で 6 個であり，流域全体をグリッドで覆う

流出予測システムの状態量を観測流量で同化する方法

タイプの分布型流出モデルを対象とはしていない．

として，河道流に適用するマスキンガム−クンジモデら

23)の研究や，HBV と呼ばれる分布型流出モデル

24)の状態量を逐次更新する Wohling ら25) の研究がある．
いずれの方法も，分布型流出モデルで計算する状態量すべてをフィルタリングの対象とするのではなく，流域全体の貯留量をひとつの状態量とみなして更新する
点にその特徴がある．また，空間分布する状態量をフィルタリングの対象とする事例としては，SAC-SMA モ

ルにバイアス補正カルマンフィルタを導入する方法を一方，河道追跡モデルにフィルタリングを適用した事
27)

29)

提案する．このシステムを桂川流域に適用し，降雨流は河道網出モデルに起因するバイアスを考慮することの効果を集中型キネマティックウェーブモデルに，それぞれカル分析する．また，下流地点の観測流量を用いてフィル

マンフィルタを適用している．河道追跡モデルの状態

タリングした結果が，上流域の予測流量に及ぼす影響

量をフィルタリングの対象とするこれらの方法は，空を分析する．

例もある．例えば，Fread and Ming

や Shiiba らはダイナミックウェーブモデルに，藤田ら

28)間分布する状態量を複数地点の観測流量で同化できる

本論の構成は以下の通りである．広域分布型流出予点で優れている．ただし，予測のリードタイムが長く測システムについて 2. で概説する．次に，3. では，計なると，河道追跡の予測結果が流出モデルによる斜面

算安定化手法を適用したマスキンガム−クンジモデル

流量の予測結果に依存するため，河川の状態量を更新

を河道網全体に適用して，それを行列形式で表示する

した効果が消えてしまうという欠点がある．つまり，斜方法を述べる．4. では未知定数を含むシステム方程式

面流量の予測モデルに起因する系統的な誤差をフィルのフィルタリング法として分離カルマンフィルタを導タリングで予測・補正しなければ，河道網の状態量を出し，それをもとにバイアス補正カルマンフィルタをフィルタリングする効果は十分に発揮できない．

導出する．5. では具体的な洪水予測のアルゴリズムを従って，河道および斜面の状態量を同時に観測更新

述べる．6. では提案する手法の桂川流域 (833 km2 ) へすることが重要となる．しかし，水系全体を対象とすの適用結果を示し，7. において本論の結論をまとめる．

るような洪水予測モデルを考える場合，これらの状態

2.

量すべてを観測更新の対象とすることは，実行可能性

広域分布型流出予測システムの概要の点で困難さが付きまとう．この問題を解決する一つの方法として，斜面流出モデルが算定する流出量の予本論で用いる広域分布型流出予測システムは，落水

測誤差を，河道流追跡モデルのシステム方程式におけ線型の分布型流出モデルと河道網の河道追跡モデルでる予測バイアスとして捉え，逐次その予測バイアスを構成される．淀川流域全体を対象としたシステム5) は，推定・補正する方法が考えられる．つまり，システム方これらのモデルに加えて，琵琶湖の湖沼モデルとダム

程式に未知定数を加え，河道状態量とともにその未知モデルとを統合しているが，本論では，琵琶湖とダム

定数を逐次推定する適応カルマンフィルタを応用するの影響を考慮する必要がない流域および期間を対象とことによって，斜面流量の予測モデルに起因する誤差

して以下の議論を進める．

を逐次推定し，これをもとに河川流量を予測する方法

広域分布型流出予測システムの構成手順は次の通りが考えられる．ただし，この方法を用いたとしても本である．まず，1/25,000 の地形図に記載されている河論で取り扱う広域分布型流出予測システムでは，河道道網を抽出し，それを合流点で分割する．分割した河網を約 3 km 毎の河道区分に分割し，各河道区分への流道が 3 km 以上の場合には，その長さが 3 km 以下にな入量を斜面流量の予測モデルによって計算しているのるように河道をさらに分割する．次に，空間分解能 250

で，適応カルマンフィルタで推定すべき未知乗数の数

m の標高情報をもとに，最急勾配法で落水方向を決定
し，各河道区分に流入する部分流域を抽出する32) ．そは河道状態量の数の二倍になる．つまり，適応カルマ

227土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10

れぞれの部分流域は，矩形斜面の集合であり，各矩形斜面に不飽和・飽和中間流・表面流モデル33) を適用する．
各部分流域に適用した分布型流出モデルによって，河道区分への側方流入量を計算する．河道追跡モデルにttc < Δx /Δt

n

n

Δt

Δt

c > Δx / Δt

Δx / cは，これまでキネマティックウェーブモデルを用いてきたが，本論でカルマンフィルタを導入するにあたって，

c Δt
Δx
j-1

河川流量を状態量にとり，状態量の時間推移が形式上，線形式となるマスキンガム−クンジ法を用いた．
広域分布型流出予測システムの構成とその再現性に

n-1

n-1
j

l

j-1

Δx

j

l

図–1 Multi-step, Multi-reach 法における計算時空間ステッ
プの分割方法5)

ついては文献 に記述しているのでここでは省略し，本論に関連するマスキンガム−クンジ法とその行列表示について 3. でその詳細を述べる．

定する．c を差分区間内で一定と仮定して，式 (3)，(4)3.

マスキンガム−クンジ法による河道追跡(1)

河道区分における計算を式 (2) に代入すれば，
n−1
Qnj = C1 Qnj−1 + C2 Qn−1
+ C4 qjn
j−1 + C3 Qjを得る．ここに C1 から C4 は，c∆t/∆x − 2X
C1 =2(1 − X) + c∆t/∆xc∆t/∆x + 2X
C2 =2(1 − X) + c∆t/∆x
2(1 − X) − c∆t/∆x
C3 =2(1 − X) + c∆t/∆x
2c∆t
C4 =
2(1 − X) + c∆t/∆x

広域分布型流出予測システムの河道追跡計算にはマスキンガム−クンジ法を用いる．上述の通り，本システ
ムでは，河道の合流点および約 3 km 毎の河道地点で河道網を分割しており，その各河道区分に対して，以下に示すマスキンガム−クンジ法を適用する．なお，流域の河道網が p 個の河道区分で構成される場合には，状態量とその時間変化を示す式がそれぞれ p 個ずつ存在することになり，河道網全体の流量を計算するためには，(4) に示す行列計算を解くことになる．

である．c はマニングの抵抗則を用いれば，
Ã √ ! M1 µ ¶ M −1
i
Qo M
c=M
N
B

ここではひとつの河道区分を対象に，時刻 n，n − 1
n−1の河道区分上端での流量 Qn
j−1 ，Qj−1 ，時刻 n − 1 の下端での流量 Qjn−1 ，および単位河川長あたりの側方流
入量 qjn (n − 1 から n の間の平均側方流入量) をもとに，
流量 Qn
j を求める．ここで，Q や q の上添え字と下添
え字はそれぞれ時間ステップと空間ステップを表しており，ここでは時刻 n − 1 の初期条件と空間 j − 1 の境界条件が既に得られていることを前提にして議論を進める．
流量を Q，通水断面積を A とすると連続式は，
∂A ∂Q
+
=q
(1)
∂t
∂x
となる．ここで流れの伝播速度 c を dQ/dA とすれば，

∂Q
+c
∂t

µ

∂Q
−q
∂x

¶
=0

(2)

(6)
(7)
(8)
(9)

(10)

となる．ここに，M : 無次元の流量流積パラメタ ( =5 / 3)，i : 河道の勾配，N : Manning の粗度係数，
B : 河道幅，∆t : 計算の時間ステップ，∆x : 計算の空間ステップである．Qo は参照流量であり，本論ではQo = Qn−1
j−1 とする．また式 (6) から (9) にある X は，
µ
¶
Qo
1
X=
1−
(11)
2
Bci∆x
とすることにより，数値拡散を物理的な拡散として扱うことができる34) ．

(2)

式 (1) は

(5)

Multi-step, Multi-reach 法 35) による数値計
算の安定化

マスキンガム−クンジ法は，係数 C1 ，C2 ，C3 ，C4のように変形できる．マスキンガム−クンジ法ではこを既知の値を用いて設定することにより，形式上，線形の式の各項を次のように差分近似する．

の式で河川流量を計算できる点にその特長がある．ただ

∂Q
=X
∂t

Qnj−1 − Qn−1
j−1∆t

+ (1 − X)n−1

Qnj − Qn−1j∆tし，∆x と ∆t のとり方によっては Ci (i = 1, . . . , 4) や

(3)

X が負の値になって，計算が不安定になることがある．

n−1

Qj − Qj−1
Qnj − Qnj−1
∂Q
=
+
(4)
∂x
2∆x
2∆x
X は空間差分の重み付けを規定するパラメタであり，
マスキンガム−クンジ法では X を式 (11) のように決228

陸ら35) はマスキンガム−クンジ法による河道追跡計
算の安定性を確保する方法として，Multi-step, Multi-

reach 法を提案しており，本論でもこの方法を採用することにした．Multi-step, Multi-reach 法の要点は，

土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10∆x/∆t が流れの伝播速度 c とほぼ等しくなるように ∆x と ∆t を決定することにあって，これによりCi (i = 1, . . . , 4) や X が負の値になるのを防ぐ．本論の適用例では，∆x と ∆t の初期設定値をそれぞれ各
河道区分の長さ (約 3 km)，1 時間としておき，流れの伝播速度 c が ∆x/∆t よりも小さい場合には，計算の空間ステップを c∆t として短くとる．一方，流れの伝播
速度 c が ∆x/∆t よりも大きい場合には，計算の時間ス
テップを ∆x/c として短くとる．

q3
q2
q1

n

n

Q2

n

n

Upstream

Q3n
Q4n

q4

n

q5

n

Q5n

Q1n Downstream
図–2 5 つの河道区分で構成された河道網の例

図–1 に示した例は，黒丸で示した既知の値 Qn−1
j−1 ，
n
Qn−1
，Qn
j−1 から白丸で示した未知の値 Qj を求めて
jいる．図–1 左が c < ∆x/∆t の場合である．Multi-step,Multi-reach 法では，式 (5) の要領で ((j − 1) + c∆t) 地点の流量をまず求め，それを上流の境界条件として，次に ((j − 1) + 2c∆t) 地点の流量を求める．これを繰り返せば，求めるべき地点の流量 Qn
j を得る．なお，計算の空間ステップを分割する場合には，(j − 1) 地点と j
地点の間の初期流量が必要となるが，ここでは，Qn−1
j−1
n−1
と Qj
との流量を線形内挿してこれらの値を求めている．図–1 右が (c > ∆l/∆t) の場合で時間ステップを分割しており，その手順は空間ステップを分割する場
合と同様である．

µ¶
om+C3 1 −
; m = 1, . . . , nl , A2,0 = 0nlnm−1
A3 = A3,nl ; A3,m = C1 A3,m−1 + C2no lm+C3 ; m = 1, . . . , nl , A3,0 = 0nl
1 − C1nl
A4 = C4
1 − C1

(14)

(15)
(16)

となる．ここに，nl (= ∆l/c∆t) は Multi-step, Multi-

行列の次元や時間ステップの幅が変化することになり，

reach 法における空間分割数であり，A2 , A3 は漸化式を解くことにより得られる．同様に，時間ステップを
nt (= c∆t/∆l) 分割する場合の係数は，
nm−1
A1 = A1,nt ; A1,m = C3 A1,m−1 + C2no tm+C1 ; m = 1, . . . nt , A1,0 = 0)
(17)
nt
¶
µ
n
m−1
A2 = A2,nt ; A2,m = C3 A2,m−1 + C2 1 −
nt
µ
¶
o
m
+C1 1 −
; m = 1, . . . , nt , A2,0 = 0(18)
nt
A3,m = C3nt
(19)nt
1 − C3A4 = C4(20)
1 − C3

そのままではカルマンフィルタを適用することができとなる．A1 ，A2 は漸化式を解いて得られる．これによなくなる．従って，次のような工夫が必要となる．

り，Multi-step，Multi-reach 法を適用した各河道区分(3)

時空間ステップを変更する場合の係数について

Multi-step, Multi-reach 法は，計算の時空間ステップを変化させることにより，河道追跡計算の安定性を確保する．計算の時空間ステップが変化することは，上流から下流に向かって順次計算を進める場合には問題ではない．しかし，本論のようにカルマンフィルタを適用する場合には，状態量の時間変化を行列形式で表現する必要がある．Multi-step, Multi-reach 法を適用すると，のマスキンガム−クンジ式をの線形の式 (12) で表すこ河道区分を空間的に分割する図–1 の (c < ∆l/∆t) のとができる．

ケースをまず考えてみる．Multi-step, Multi-reach 法
では，上述の通り c∆t ずつ下流に移動しながら流量を求めていく．従って，最終的に求まる下端の流量 Qn
j は

(4)

式 (5) を満たさなくなってしまう．

マスキンガム−クンジ法の河道網全体への適用

上ではひとつの河道区分を対象にした場合のマスキ

ただし，式 (5) を解析的に解けば，未知量である Qnj

ンガム−クンジ法による河道追跡計算法を述べた．次と既知量であるその他の流量との関係を，式 (5) と同じに複数の河道区分が存在する河道網を対象にした場合に，どのように式 (12) を連立させ，河道網全体の流量線形の式 (12) で記述することができる．

を計算するかについて説明する．ここでは簡単のため

n−1
Qnj = A1 Qnj−1 + A2 Qn−1
+ A4 qjn (12)j−1 + A3 Qj

図–2 に示す河道網を対象にして以下の議論を進める．

そのときの式 (12) の係数 Ai (i = 1, . . . , 4) は，

図–2 は 5 つの河道区分が河道網を構成している例である．それぞれの河道区分には式 (12) を適用する．例

A1 = C1nl
(13)
µ
¶
n
m−1
A2 = A2,nl ; A2,m = C1 A2,m−1 + C2 1 −nl

えば，時刻 n の河道区分 3 の下端流量を Qn
3 ，河道区
分 4 の下端流量を Qn
4 とすれば，時刻 n の河道区分 2229

土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10の下端流量 Qn
2 は，フィルタでは ξ n の白色正規性を仮定して，その期待値¡¢
Qn2 = A21 (Qn3 + Qn4 ) + A22 Qn−1
+ Qn−1
3
4+A23 Qn−1+ A24 q2n2と分散共分散行列を次のように設定する．

(21)となる．なお，A の肩につけた 2 という数字は，それE{ξ n } = 0(25)E{ξ n ξ Tm } = Qn δnm(26)が河道区分 2 を対象にして求めた係数であることを意

ここに，E{·} は期待値操作を表す．δnm はクロネッカー

味している．同様の式を各河道区分で構成すれば，河のデルタであり，Qn はシステムノイズの大きさを表す

道網全体の河川流量を次のように行列表示できる．
 n 
1 −A11000Q1 n 
22
 0
1−A1 −A1
0   Q2 



 0
 n 
0
1
0
0 

  Q3 



0
0
1
−A41   Qn4  0(p × p) 次元の分散共分散行列である．

0


0

0

0

1



一方，観測ベクトル y n (r ×1 次元) は，係数行列 H n(r × p 次元) を用いて以下のように表す．

y n = H n xn + η n(27)

ここに，η n は観測ノイズであり，システムノイズと同Qn5様，白色正規性を仮定して次のように設定する．

Qn−1
A13 A12 0
0
0
1



 0 A23 A22 A22 0   Qn−1  2n−1 

3
 Q

=
0
0
A
0
0
3
 3



n−1 
4
4 
0
0 A3 A2   Q4
 0

n−1
5
0
0
0
0 A3
Q5
 1
 n 
A4 0
0
0
0
q1



 0 A24 00
0   q2n 



 n 
+
0 A34 00 
 0
  q3  (22)
00 A44 0   q4n 
 0
0
0
0
0 A54
q5n

E(η n ) = 0(28)

E(η n η Tm ) = Rn δnm(29)

ここに，Rn は観測ノイズの大きさを示す (r × r) 次元の分散共分散行列である．
式 (24) をシステム方程式，式 (27) を観測方程式として，観測量をもとに状態量を逐次推定するカルマンフィ
ルタのアルゴリズムは Kalman8) によって提案され，それ以降，降雨流出予測の分野にも数多く適用されてきた．ただし，式 (24) における定数項 D n dn を河道区
間ごとに分布型流出モデルで算定し，その結果を既知

式 (22) に示したように，Qn
j を下流から上流に向かっの定数として与える本論の適用例では，D n dn の推定て順にならべて行列表示すれば，左辺の係数行列は上によって生ずる誤差バイアスを考慮にいれた上で，状三角行列になるので，逆行列を容易に求めることがで態量 (河川流量) を推定・予測することが望ましい．そきる．その逆行列を両辺に左側からかけた式を，行列こで本論では，式 (24) と式 (27) に未知のバイアス項 β

表示すれば，(p × 1 次元) を加えた式 (30) と式 (31) をシステム方程式と観測方程式とし，状態量と未知のバイアスを併せて推定するカルマンフィルタを導入することを考えた．

xn = An xn−1 + D n dn

(23)となる．ここに x は各河道区分の流量を要素に持つベクトルであり，d は降雨流出モデルによって計算される各河道区分の側方流入量を要素に持つベクトルであ

xn = An xn−1 + D n dn + B n β + ξ n(30)

y n = Hxn + C n β + η n(31)

る．これをシステム方程式として，次章で述べるカル

ここに，B n は (p × p) 次元の係数行列である．基礎式

マンフィルタの一連の式を適用する．

に一般性を持たせるため，観測方程式にも C n β (C :r × p 次元の係数行列) の項を加えて以下の導出を行う．

4.

バイアス補正カルマンフィルタの導出(1)問題の定式化(2)

適応カルマンフィルタ

未知定数 β を含むシステム方程式に対しては，以下に示す適応カルマンフィルタを利用するのが一般的な

マスキンガム−クンジ法の状態ベクトル xn (p × 1 次解法である36) ．これは，推定すべき状態量 xn と未知定元) は，時刻 n とともに線形的に変化する量であり，斜面部の流出モデルで計算する側方流入量を dn (p × 1 次

数 β n を併せて新たに拡大したベクトル z n = [xTn β Tn ]T元) として，モデルの不確かさを表現するシステムノイ(2p × 1 次元) を定義し，この z n をカルマンフィルタで次のように推定する．

ズ項を ξ n を加えれば，

xn = An xn−1 + D n dn + ξ n

0 0
ẑ −
n = Φn ẑ n−1 + D n dn

(24)

(32)

0
T
P−
(33)
n (z) = Φn P n−1 (z)Φn + Qn
h
i−1
−
T
T
K(z) = P −(34)n (z)Ln Ln P n (z)Ln + Rnのように表すことができる．ここに，An と D n は (p×p)の次元を持つ係数行列である．また，通常のカルマン

230

土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10£−¤ẑ n = ẑ −n + K(z) y n − Ln ẑ n

(35)

P n (z) = [I − K(z)Ln ] P −
n (z)

(36)

£
¤
−
x̃n = x̃−n + K̃(x) y n − H n x̃nhi −P̃ n (x) = I − K̃(x)H n P̃ n (x)ここに，ẑ n (2p × 1 次元) は z n の推定量であり，K(z)

(44)
(45)

ここに，x̃n はバイアスを無視して推定した状態量ベク(2p × r 次元) はカルマンゲイン，P −n (z) は z の事前分
散共分散行列である．上式の – の記号は z n や P n
を前時間ステップ (n − 1) の情報をもとに推定したことを意味する．また，Φn (2p × 2p)，Ln (r × 2p)，D 0n d0n(2p × 1)，Qn (2p × 2p) は，それぞれ状態推移行列，状態量から観測量への変換行列，既知の定数ベクトル，システムノイズの分散共分散行列であり，以下のように
定義する．
"
#
h
i
An B nΦn =Ln = H n C n
(37)

I
"
#
"
#
D n dnQn 
0 0
0
D n dn =
Qn =
(38)

 

トルであり，P̃ n (x) はバイアスを無視して推定した x̃nの分散共分散行列である．
次にバイアス β が既知である場合を想定する．バイアスが既知であると仮定した場合の状態ベクトルの事−

前推定値 X̂ n と事後推定値 X̂ n は，−X̂ n = An X̂ n−1 + D n dn + B n β
(46)
h
i
−
−
X̂ n = X̂ n + K̃(x) y n − H n X̂ n − C n β (47)のように表される．バイアス項 β は既知の定数として扱うため，事前および事後の分散共分散行列とカルマ
ンゲインは式 (42)，(45)，(43) と同様である．

b)

パラメタ U n , V n , S n の定義と計算法

バイアスを無視した場合のフィルタリング前後の推定

適応カルマンフィルタは，式 (30)，(31) で表される
状態空間モデルにおいて，最適推定量を得るための一量 x̃−
n ，x̃n と，バイアスが既知であると仮定した場合

般的な解法である．しかし，状態量ベクトルとそのバのフィルタリング前後の推定量 X̂ n ，X̂ n との関係を，

イアスベクトルがそれぞれ (p × 1) の次元を持ち，かつX̂ n = x̃−n + U nβ(48)

p が大きくなる本論の問題においては，計算する行列の次元が (2p × 2p) に拡大するため，計算負荷が大きくなるという欠点を持つ．

X̂ n = x̃n + V n β

(49)

−

−

と仮定し，係数行列 U n と V n (p × p 次元) を求めることを考える．
−(3)

X̂ n と x̃−n との関係は，式 (41)，(46)，(48)，(49) から，

分離カルマンフィルタ

この問題を克服するため，Friedland30) は行列の次元−X̂ n − x̃−n = U nβ³´= An X̂ n−1 − x̃n−1 + B n β

を (2p × 2p) に拡大することなく，状態量ベクトル xnとそのバイアスベクトル β を (p × p) の次元のままで順次求める方法を提案している．この方法は，分離カル= (An V n−1 + B n ) β

マンフィルタと呼ばれ，Friedland の提案以降，より一となるので，U n と V n−1 との関係は，般化した問題への展開37),38) や，Friedland とは異なっ

U n = An V n−1 + B n

た方法でのアルゴリズムの導出が進められてきた．以30)

下では，Friedland
した Igagni

39)に比べて，簡潔的な導出方法を示(51)と表される．

の方法に従って，分離カルマンフィルタ

同様に，X̂ n と x̃n との関係は，式 (44)，(47)，(48)，の導出過程を追う．

a)

(50)

(49) から，

バイアス項 β を無視した場合と既知の場合とのカ

X̂ n − x̃n = V n β

i
h
−
−
−
= X̂ n − x̃−
n + K̃(x) −H n (X̂ n − x̃n ) − C n βhi= U n − K̃(x) (H n U n + C n ) β(52)

ルマンフィルタの比較
分離カルマンフィルタを導出するにあたり，まず式(30) における未知定数の項 (バイアス項) を無視した場
合のシステム方程式と観測方程式を想定する．
xn = An xn−1 + D n dn + ξ n(39)

y n = H n xn + η n(40)となる．いま S n (r × p 次元) を，

Sn = H nU n + C n(53)と定義すれば，V n と U n との関係は，この問題に対し，バイアス項を無視した通常のカルマ

V n = U n − K̃(x)S n

ンフィルタは次のようになる．

(54)

(41)

となる．式 (51)，(53)，(54) をみれば，V n−1 をもとにP̃ n (x) = An P̃ n−1 (x)ATn + Qn
(42)
h
i−1
−
−
K̃(x) = P̃ n (x)H Tn H n P̃ n (x)H Tn + Rn
(43)

U n , S n , V n を求める式になっており，この三つの式を再帰的に解くことにより，式 (48)，(49) における U nと V n を求めることができる．

x̃−n = An x̃n−1 + D n dn−

231

土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10c)

β を推定するためのカルマンフィルタ

以上をまとめると，分離カルマンフィルタは，式 (41)から (45) でバイアスを無視した推定量を求め，式 (55)未知の定数項 β を逐次推定する方程式を以下に導出から (64) でバイアスを推定し，式 (65) から (68) でバする．まず，β の伝達方程式は，−

β̂ n = β̂ n−1

イアスを補正する一連のアルゴリズムであった．なお，(55)

分離カルマンフィルタの推定量は，適応カルマンフィルとする．β の観測方程式については，以下のように r nタによる推定量と理論的に一致することが Friedland30)(r × 1 次元) を定義し，によって証明されている．

r n = y n − H n x̃−hn −i= y n − H n X̂ n − U n β

(4)

(56)

分離カルマンフィルタは，行列の次元を拡大するこ式 (56) の右辺第 1 項と第 2 項にそれぞれ −C n β と C n βとなく，状態量と未知定数を順次推定できる．ただし，を加えれば，

実際の問題で合理的な推定量を得るためには，P (x) や³´−
r n = y n − H n X̂ n − C n β(57)+ (H n U n + C n ) β(58)= S n β + vn

P (β) の初期値を適切に設定する必要がある．これは，上述のアルゴリズムを工学的な問題に適用するうえで十分に留意すべき点である．

(59)となるので，この式 (59) を β の観測方程式と見立てる．
−

ここに，v n は (y n − H n X̂ n − C n β) で定義する r 次
元のベクトルであり，v n の分散共分散行列は次のように表される．
−E{v n v Tn } = H n P̃ x (n)H Tn + Rn

(60)

これより β を逐次推定するためのカルマンフィルタ
は次のようになる．

P−n (β) = P n−1 (β)(61)TK(β) = P −n (β)S n ×
i−1
h
−
T
(62)
H n P̃ n (x)H Tn + S n P −n (β)S n + Rnhi−−β̂ n = β̂ n + K(β) r n − S n β̂ n
(63)

P n (β) = [I − K(β)S n ] P −n (β)

バイアス補正カルマンフィルタ

Dee and Da Silva31) は，推定結果が P (β) の初期値に大きく依存するという分離カルマンフィルタの問題点を指摘したうえで，バイアスの分散共分散行列 P (β)
を P (x) から簡単な方法で求め，より安定したフィル
タリングのアルゴリズムを提案している．そのアルゴ
リズムは，分離カルマンフィルタと同様に行列の次元を抑えて計算し，かつより少ない手順で状態量とバイアスを推定することができる．本論では，このアルゴ
リズムをバイアス補正カルマンフィルタとよび，分離カルマンフィルタの一連の式からバイアス補正カルマ
ンフィルタを導出する．
a)バイアス b の定義

まず，ここからの議論で用いるバイアス (p × 1 次元)
を，

(64)

−
b−
n ≡ E[x̃n − xn ]−上式の中で，P̃ n (x) と r n はバイアスを無視したカルマ(69)

ンフィルタで求めることができる．また，S n は式 (51)，と定義する．bn は x の推定量 x̃ に含まれるバイアスで(53)，(54) で求めることができる．つまり，バイアスを無視したカルマンフィルタで x̃ を推定し，バイアス推定のためのカルマンフィルタで β を推定し，それらの値を式 (48) と (49) に代入すれば，バイアスを補正した
状態の事前推定値 x̂−n と事後推定値 x̂n を得る．

あり，分離カルマンフィルタの説明で用いた β とは異

−

d)

−
x̂−
n = x̃n + U n β̂ n(65)x̂n = x̃n + V n β̂ n(66)

なっていることに注意する．式 (65) と (69) から b−n とβ との間には，−b−n ≡ E[x̂n − xn ] − E[U n β̂ n−1 ]= −U n β

定量は，−

バイアス補正後の状態量ベクトルとバイアスベクト
ルの分散共分散行列は以下の式で計算できる39) ．

(67)

P̂ n (x) = E{(x̂n − xn )(x̂n − xn )}= P̃ n (x) + V n P n (β)V Tn−b̂n = −U n β̂ n(72)b̂n = −U n β̂ n(73)となる．また，バイアス推定量の分散共分散行列は，−−
TP̂ n (x) = E{(x̂−n − xn )(x̂n − xn ) }−(71)の関係があるので，フィルタリング前後のバイアス推

バイアス補正後の分散共分散行列

T= P̃ n (x) + U n P −n (β)U n

(70)

(68)

となる．

232

−
T
P−
n (b) = U n P n (β)U n(74)

P n (b) = U n P n (β)U Tn(75)土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10b)

状態量とバイアスの更新の関係がある．これは，バイアスを無視して推定した

式 (53) と (54) から V n は，状態量の分散共分散行列 P̃ n (x) に，バイアスの分散共V n = U n − K̃(x)[H n U n + C n ]−分散行列 P −n (b) を加えた値が，バイアス補正後の状態

(76)

−

量の分散共分散行列 P̂ n (x) となることを意味する．

となり，観測にはバイアスがないと仮定 ( C n = 0 ) す

分離カルマンフィルタは，この P −n (b) に相当するバ

れば，

イアスの分散共分散行列を逐次更新するアルゴリズム

V n = [I − K̃(x)H n ]U n

(77)となっているが，バイアス補正カルマンフィルタでは，−P−n (b) が P n (x) と同じ空間構造で生ずると仮定し，式(83) を用いて P −n (b) を簡単に推定する．

となる．この V n を式 (66) に代入することで，バイアス補正後の状態量 x̂n の推定式を以下のように得ることができる．

−

P−n (b) = γ P̂ n (x)x̂n = x̃n + V n β̂ n
= x̃n + [I − K̃(x)H n ]U n β̂ n
= x̃n − [I − K̃(x)H n ]b̂n量のパラメタであり 0 から 1 の値をとる．γ が 1 に近−= x̃−n + K̃(x)[y n − H n x̃n ]ければ，バイアス推定のカルマンゲインが相対的に大きくなり，観測量と予測量との差の情報は主にバイア−[I − K̃(x)H n ]b̂n

スの更新に用いることになる．一方，γ が 0 に近けれ= x̃−n − b̂n+K̃(x)[y n − H n (x̃−n − b̂n )]

(83)

−
−
P̃ n (x) = (1 − γ)P̂ n (x)(84)−
ここで，γ は P −n (b) と P̃ n (x) の割合を決めるスカラ

ば，状態量のカルマンゲインが相対的に大きくなり，観(78)

測量と予測量との差の情報は主に状態量の更新に用い

ここで，x̃−n はバイアス補正前の予測状態量である．前ることになる．なお，γ が 0 のときはバイアスは更新

述の分離カルマンフィルタでは x̃−n を式 (41) で求めた．

されないので，バイアス補正カルマンフィルタの推定

つまり，バイアスを無視したカルマンフィルタを実行

結果はバイアスを考慮しない通常のカルマンフィルタ

することにより x̃n−1 を求め，それを初期値として時になる．

刻 n の状態量を予測した．一方，バイアス補正カルマ
ンフィルタでは，x̃−n を求める際に以下の式を用いる．

x̃−n = An x̂n−1 + D n dn

(5)

バイアス補正カルマンフィルタのまとめ

以下にバイアス補正カルマンフィルタのアルゴリズ(79)

ムをまとめる．

これは，バイアス補正後の状態量 x̂n−1 を初期値として，x̃−n を予測することを意味する．これにより，最新• 状態量とバイアスの予測のバイアス推定量を状態量の伝達式にフィードバック
するとともに，バイアスを無視したカルマンフィルタのアルゴリズムを省略できるので，分離カルマンフィ

−

(87)

• 分散共分散行列の予測−P̃ n (x) = An P̃ n−1 (x)ATn + Qx (n)−1−P̃ n (x)P̂ n (x) =
1−γ
γ
−
−
P̃ n (x)P n (b) =1−γ• バイアスの更新

T
K(b) = P −n (b)H n ×−−1[H n P̃ n (x)H n + H n P −(80)
n (b)H n + Rn ]として求めることができる．同様に，式 (63) の両辺に−U n を掛けることで，b̂n の更新方程式は次のように
なる．

(88)
(89)
(90)

T
−
T
K(b) = P −n (b)H n [H n P n (b)H n

(81)

−

以上がバイアス補正カルマンフィルタにおける，状態量 x̂n とバイアス b̂n の推定法である．

+H n P̃ n (x)H Tn + Rn ]−1

(91)

−
−
b̂n = b̂n − K b [y n − H n (x̃−n − b̂n )](92)• 状態量の更新

状態量とバイアスの分散共分散行列−

K̃(x) = P̃ n (x)H Tn

バイアスを補正した状態量の分散共分散行列は，式−×[H n P̃ n (x)H Tn + Rn ]−1(67) と (74) から，−−P̂ n (x) = P̃ n (x) + P −
n (b)

(86)

−
x̂−
n = x̃n − b̂nbn を推定するためのカルマンゲイン K(b) は，式(62) の両辺に −U n を掛けて，K(b) を U n K(β) と定
義すれば，

c)

(85)

−
b̂n = b̂n−1• 予測結果のバイアス補正

ルタに比べて計算の手順が少なくなる．

−−
b̂n = b̂n − K(b)(y n − H n x̃−
n + H n bn )x̃−n = An x̂n−1 + D n dnx̂n = (x̃−n − b̂n )

(82)

233

(93)

土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10+K̃(x)[y n − H n (x̃−n − b̂n )] (94)

Katsura River Basin(833 km2)• 分散共分散行列の更新−P̃ n (x) = [I − K̃(x)H n ]P̃ n (x)

5.

(95)

洪水予測のアルゴリズム

Lake
Biwa

マスキンガム−クンジ法による河道追跡計算に，バイアス補正カルマンフィルタを導入する．以下に，そのアルゴリズムを組み込んだ広域分布型流出予測シス
テムによる流量予測の具体的な手順を述べる．現在時

Tenryuji

刻を n とし，時刻 n − 1 のフィルタリング後の状態量

Kameokaが得られている状況を想定する．

Step 1 時刻 n − 1 から n までの観測降雨を入手する．
Step 2 各部分流域に適用した分布型流出モデルを実行し，各河道区分への側方流入量 dn を計算する．
Step 3 マスキンガム−クンジモデルを実行し，時刻 nの河川流量 x̃−n を計算する．
−
Step 4 時刻 n − 1 で求めたバイアス b̂n−1 (= b̂n ) を−−もとに，状態量を補正する (x̂−n = x̃n − b̂n )．
−−Step 5 状態量の分散共分散行列 P̃ n (x)，P̂ n (x) と，バイアスの分散共分散行列 P −n (b) を計算する．
Step 6 時刻 n の観測流量 y n を入手する．
Step 7 式 (91) から式 (95) を実行し，バイアス，状態
量，分散共分散行列を更新する．
Step 8 Step 2 で得た流出モデルの状態量と，Step 7で更新した河道モデルの状態量を初期条件として，
時刻 n + 1 までの予測流量を計算を実行し，Step 6で推定したバイアスをもとに予測結果を補正する．
Step 9 Step 7 で得たバイアス補正後の状態量を初期値として，n + 2 までの予測計算を実行し，同様に予測結果をバイアス補正する．この予測計算を必要とするリードタイムまで繰り返す．

KatsuraGauging Sta.

Yodo River Basin(7,281 km2)
10 km

図–3 桂川流域 (833 km2 ) における流量観測所の位置

本論の適用例では，1992 年 8 月 18 日 0 時から 23 日

0 時までの 5 日間の洪水イベント (1992 年イベント, 流
域平均雨量: 140 mm) と，1993 年 8 月 14 日 0 時から
19 日 0 時までの 5 日間の洪水イベント (1993 年イベン
ト, 流域平均雨量: 175 mm) の二つを計算対象洪水とした．1992 年イベントは一山洪水であり，1993 年イベ
ントは複数のピークが連続する洪水である．
入力降雨については，流域内 13ヶ所で観測された地
上雨量 (時間分解能: 1 時間) を最近隣法で空間内挿したものを用いる．フィルタリング手法の性能を明らかにすることが本論の目的なので，降雨の予測が完全に6.

桂川流域への適用(1)

対象流域と計算条件

なされた状況を想定して，予測降雨についても実績降雨と同じ降雨分布を用いた．
状態量の更新に用いる河川流量は，特に断らない限

バイアスカルマフィルタを導入した広域分布型流出り，桂地点一地点の観測流量 (時間分解能: 1 時間) と予測システムを，桂川流域 (桂地点上流: 833 km2 ) にし，(2) b) で亀岡地点の流量をフィルタリングに用い適用する．計算対象期間の 1992 年と 1993 年当時は上た例を示す．

流域に主要なダムが存在しないので，自然流況で予測性能を検証できる．また，図–3 に示すように，主河道

バイアス補正カルマンフィルタで加えるノイズにつに沿って複数の流量観測所が存在するので，流域下端でいては，試行錯誤的に決定し，システムノイズの分散をフィルタリングした効果の上流に及ぼす影響や，複数

10 [m6 /s2 ]，観測ノイズの分散を 100 [m6 /s2 ] とした．
バイアスの分散共分散行列の大きさを規定するパラメ
タ γ はまず 0.5 に設定し，その値を変更した場合の予測結果への影響を (2) c) で検討した．広域分布型流出予測システムのパラメタは既往の研究5) で同定したパラ
メタと同じ値を用いた．

地点の流量を更新に用いた場合の効果などを分析できる．なお，広域分布型流出予測システムは桂川流域の河道網を 193 の河道区分に分割しており，バイアス補正
カルマンフィルタの主な行列計算の次元は (193 × 193)
になる．

234

土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10

(2)

結果と考察

a)

バイアス補正の効果

1412

図–4 に 1992 年の洪水イベントを対象にしてバイア

Rainfall [mm]10ス補正カルマンフィルタを適用した結果を示す．リード
タイムをそれぞれ 0, 1, 2, 3 時間として予測したときの

8
6桂地点の流量を示している．また，バイアス補正カル

4マンフィルタの結果 (Bias filter) と比較するために，バ2

イアスを考慮しない通常のカルマンフィルタ (Kalman0

まず，フィルタリング直後の解析結果に着目すると，バイアス補正カルマンフィルタと通常のカルマンフィル
タとの違いはほとんど見られず，いずれも観測流量に

0

800

20

40

60
Time [hr]

80

100

120

'92, Katsura, Updated estimate
Bias filter
with observation Kalman filter
Observed
Simulated
Bias

700
600
Discharge [m3/s]

filter) と，オフラインシミュレーション (Simulated) の結果を併せて示す．リードタイム 0 時間の予測，すな
わちフィルタリング直後の解析結果を示すパネルには，推定した桂地点のバイアス (Bias) も併記している．

500
400
300
200

近い値をとっていることがわかる．フィルタリングを実100行した後の河川流量は，システムノイズと観測ノイズ

00の大きさをいかにとるかに依存しており，この計算例

20

40

60

80

100

120

Time [hr]800

したバイアスの時系列は，概ねオフラインシミュレー

700

ションと観測流量の差の時系列に類似している．バイ

600

アスの大きさは，例えばピークに近い計算開始後 40 時
間の結果で見れば，バイアスが約 70 m3 /s，オフライ
ン計算値と観測流量との差が約 100 m3 /s であるから，

Discharge [m3/s]

では解析結果が観測流量に近い値になっている．推定'92, Katsura, 1-hr ahead prediction

Bias filter
Kalman filter
Observed

500
400
300
200

オフライン計算値と観測流量との差の約 7 割程度で推100移していることになる．また降雨終了後では，計算流

00量と観測流量との差の時間的変化が小さいために，推定されたバイアスの値も概ね同じ値で推移することが

800

分かる．

700

次に予測結果に着目する．リードタイムが 1 時間の

20

40

60
Time [hr]'92, Katsura, 2-hr ahead prediction

80

100

120

Bias filter
Kalman filter
Observed

Discharge [m3/s]

600

場合，バイアスを考慮することによって，洪水の立ち上がり部や逓減部の予測が観測流量に近づく．また，ピー
ク直前もバイアスを考慮することによって予測結果が

500
400
300
200

下方修正され，観測流量に近づいている．

100

予測のリードタイムが 3 時間の場合に着目すると，バ0

イアスを考慮しない通常のカルマンフィルタの予測は

0800

ルタリングの効果が見られないことがわかる．これは，700

観測情報をもとに河川流量を更新したとしても，3 時間

600
Discharge [m3/s]

オフラインシミュレーションの結果とほぼ一致し，フィが経過すれば，河川の流量が斜面部の流出モデルに大きく依存するためである．それに対して，バイアス補正
カルマンフィルタは流出モデルに起因するバイアスを

20

40

60
Time [hr]'92, Katsura, 3-hr ahead prediction

80

100

120

Bias filter
Kalman filter
Observed

500
400
300

逐次推定して，河川流量の予測を補正するので，リー200

ドタイムが 3 時間の場合でもフィルタリングの効果が

100
0

現れる．

0

同様の傾向は図–5 に示した 1993 年イベントの結果でも見られる．特にオフラインシミュレーションの結果

20

40

60
Time [hr]

80

100

120

図–4 1992 年イベントを対象にしたフィルタリング後の流量と予測流量 (リードタイム: 1-hr, 2-hr, 3-hr)

が継続的に観測流量よりも大きくなっている計算開始か

235

14

800

12

700

10

600
Discharge [m3/s]

Rainfall [mm]土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10

86
4

Bias filter'92, Tenryuji, 3-hr ahead prediction
Bias filter (two stations)Observed

500
400
300
200

2

100
0
0

20

40

60

80

100

120

0
0

Time [hr]'93, Katsura, Updated estimate
with observation

450

Discharge [m3/s]

400Bias filter
Kalman filter
Observed
Simulated

300

80

100

120

250
200
150
Bias

ら 40 時間後以降は，バイアスを考慮する効果がはっき50りと現れる．ただし，30 時間後のピーク前後では，洪

0
0
500

20

40

60
Time [hr]'93, Katsura, 1-hr ahead prediction

450

Discharge [m3/s]

60

図–6 桂地点の流量のみをフィルタリングに用いた場合と，桂地点・亀岡地点の流量をフィルタリングに用いた場合の，天竜寺地点における 3 時間先の予測流量の比較

350

100

80

100

120

水の立ち上がり部で推定したバイアスが正の値をとっBias filter
Kalman filter
Observedているのに対し，ピーク時のオフラインシミュレーショ

400

ンが観測結果に近いので，予測結果を小さく見積もる

350

結果になっている．

300
250

b)

200
150

上流地点におけるフィルタリングの効果

図–6 は，リードタイムを 3 時間とした場合の天竜寺

10050地点における予測結果である．桂地点の約 5 km 上流に0

天竜寺地点は位置する．桂地点だけで流量を更新する

0

20

40

60

80

100

120

場合と天竜寺から約 10 km 上流の亀岡地点でも流量をTime [hr]500'93, Katsura, 2-hr ahead prediction

450

更新する場合について，バイアス補正カルマンフィル

Bias filter
Kalman filter
Observed

タによる予測結果を示している．

400
Discharge [m3/s]

40Time [hr]

500

350

桂地点だけで流量を更新する場合にまず着目すると，

300
250

バイアス補正カルマンフィルタの結果は図–4 に示した

200

オフラインシミュレーションのそれより観測流量に近

150
100

づいていることが分かる．これは，下流桂地点で流出

50

モデルのバイアスを検知したことが，上流地点の予測

0
0
500

20

40

60
Time [hr]'93, Katsura, 3-hr ahead prediction

450

80

100

にも効果を及ぼすことを意味する．また，上流の亀岡

120地点にも観測情報がある場合には，さらに予測結果がBias filter
Kalman filter
Observed

向上する．

400
Discharge [m3/s]

20

本論で行ったシミュレーションはフィルタリングの

350
300

検証を目的としているので，あえて用いる観測流量の250

地点数を制限して計算した．しかし，洪水予測の実務

200においては，可能な限り多地点の信頼できる観測流量

150
100

を用いてフィルタリングすればよい．フィルタリング

50の効果は流量観測が無い上下流地点にも及ぶので，従

0
0

20

40

60

80

100

来洪水予測の難しかった中小河川においても，バイア

120Time [hr]ス補正カルマンフィルタを導入した分布型流出モデルを用いることにより，確度の高い洪水予測を実現でき

図–5 1993 年イベントを対象にしたフィルタリング後の流量と予測流量 (リードタイム: 1-hr, 2-hr, 3-hr)

る可能性がある．

236

土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10500'93, Katsura, 3-hr ahead
prediction

450

400

350
300
250
200
150
Bias (γ = 0.25)
Bias (γ = 0.75)

100

350
300
250
200
150
100
50

50
0

0
0

20

40

60

80

100

120

0

20

Time [hr]

40

60

80

100

120

Time [hr]図–7 パラメタ γ による 3 時間先の予測流量とバイアス推定の違い (桂地点, 1993 年イベント)

図–8 適応カルマンフィルタによる 3 時間先の予測流量 (桂地点, 1993 年イベント)

d)
c)

0-hr
1-hr
2-hr
3-hr
Observed
Simulated'93, Katsura, Adaptive Filter

450

Discharge [m3/s]

Discharge [m3/s]

400

500

Bias filter (γ = 0.25)
Bias filter (γ = 0.75)
Observed

パラメタ γ の影響について

適応カルマンフィルタとの比較

バイアス補正カルマンフィルタは，4. に示したように，
適応カルマンフィルタから派生した方法であり，両者バイアス補正カルマンフィルタには，バイアスの分の予測結果の違いを比較しておく必要がある．図–8 は，散共分散行列の大きさを決める γ というパラメタがある．4. (4) で前述の通り，γ は観測量と予測量との差をバイアスの補正に用いるか，状態量の補正に用いるかを決定するパラメタである．γが 0 に近ければ状態量のカルマンゲインが，1 に近ければバイアスのカル
マンゲインが大きくなる．分布型流出モデルと河道モ
デルの不確実性が定量的に明らかであれば，γ の値を理論的に決定することができると考えられる．しかし，両者のモデルの不確実性評価に関する研究は現状では発展途上にあり，上述のシミュレーションでは，分布型
流出モデルと河道モデルによる相対的な誤差が概ね同じ大きさで出現することを仮定して γ の値を 0.5 に設定することにした．ここでは，この γ の値を変えて予測計算を行った場合に，それが予測結果に及ぼす影響を考察する．
図–7 は 1993 年の洪水を対象にしたリードタイム 3
時間の予測結果である．γ を 0.25 と 0.75 に変えて計算している．γ が大きくなればバイアス推定量の分散共分散行列である P (b) が大きくなり，バイアスの推定値
は予測と観測との差に応じてより鋭敏に変化することになる．例えば，30 hr 付近のバイアス推定値に着目す
ると，流出ピークに達するまでは予測流量の方が観測

1993 年の洪水イベントを対象に適応カルマンフィルタ
による桂地点の予測結果を示している．計算開始後 30
時間後から 40 時間後に現れる洪水のピーク付近は，バ
イアス補正カルマンフィルタのように過小評価をすることもなく，適応カルマンフィルタの予測結果は良好である．また，図–4 に示した通常のカルマンフィルタ(Kalman filter) に比べて，各リードタイムで適応カル
マンフィルタの予測が観測流量と近づいている．つまり，未知定数をシステムに加えたことによる予測結果の向上を確認できる．しかし，リードタイムが長くなるほど，予測結果がオフラインシミュレーションの結果に近づき，リードタイムが 3 時間の場合にはフィル
タリングの効果が現れない．この傾向は，通常のカル
マンフィルタと同様である．両者の手法の違いは，主に未知定数項（バイアス項）の分散共分散行列の取り扱いにあり，適応カルマンフィルタはそれを逐次更新していくのに対し，バイアス補正カルマンフィルタでは状態量の分散共分散行列とバイアス項のそれが同じ空間構造を持つと仮定して推定している．これによりバイアス補正カルマンフィルタでは行列計算の次元が状態量の数以上に拡大することがないので，より大きな流域への展開にも適している．

流量より大きく，γ = 0.75 の場合は γ = 0.25 の場合に比べてバイアス推定値が大きくなっている．その結果，

7.

おわりに

流出ピークの付近では γ = 0.75 の予測流量は観測流量に比べて小さくなっている．一方，バイアスの出現傾

広域分布型流出予測システムの観測流量データ同化

向が安定する後半部では γ = 0.75 の方が γ = 0.25 に手法として，河道網に適用したマスキンガム−クンジ

比べて予測流量が観測流量に近い．なお γ の予測流量モデルのフィルタリング法を提案した．以下に，本論に対する鋭敏性はそれほど大きなものではないことがで示した手法の概要と予測シミュレーションの結論を

図–7 の分析結果から確認できる．

まとめる．

237

土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10(a) 計算安定化手法を導入したマスキンガム−クンジ
法をもとに，空間分布する河川流量の時間推移を行列形式で表現する方法を示した．
(b) マスキンガム−クンジ法によって推定する河川流量をフィルタリングの対象とした．また，流出モ
デルに起因するバイアスをシステム方程式に加えて，河川流量とともに推定すべき対象とした．
(c) 状態量とともに未知定数を推定するフィルタリン
グ法として，適応カルマンフィルタを示したうえ
で，その代替法として，行列の次元を拡大せずに状態量と未知定数を推定する分離カルマンフィル
タの導出方法を示した．さらに，バイアスの推定問題に特化して，より計算効率を高めた方法として，バイアス補正カルマンフィルタの導出方法を説明した．
(d) バイアス補正カルマンフィルタをマスキンガム−
クンジ法に適用して，桂川流域の洪水を予測した．
その結果，通常のカルマンフィルタではフィルタリングの効果をほとんど発揮しない 3 時間先の予測についても，バイアス補正カルマンフィルタではその効果を発揮し，予測結果が全般的に観測流量に近づくことを示した．
(e) 流量観測が得られない上流の地点においても，バイアス補正カルマンフィルタの方法では，フィルタリングの効果が及ぶことを明らかにした．ただし，観測地点から離れればその効果は小さくなるので，
洪水予測の実務ではできるだけ多くの流量情報をフィルタリングに用いるべきである．提案した方法は，複数地点の流量観測情報を同じアルゴリズ
ムで同化することができて，流量情報を増やすことによって予測結果が向上する事例を示した．
以上のように，本論で提案した方法は，流出モデルに起因する予測のバイアスを，カルマンフィルタのア
ルゴリズムで状態量とともに逐次推定し，予測結果の補正を行うものである．今後はこの手法を淀川全流域体に適用するとともに，稼働中のリアルタイム予測シ
ステムに導入する予定である．
謝辞：本研究は，科学研究費補助金・若手研究 (B)18760373 (代表: 佐山敬洋)，および (財) 河川環境管理財団・河川整備基金助成事業の補助を得ました．こ
こに記して謝意を申し上げます．
参考文献
1) 国土交通省 : 水害レポート 2004, 2005,
http://www.mlit.go.jp/river/pamphlet_jirei/.
2) 宝 馨，高棹琢馬，椎葉充晴 : 洪水流出の確率予測における実際的手法, 第 28 回水理講演会論文集, pp. 415 –
422, 1984.
3) 橋本識秀，兪 朝夫，星 清 : 洪水流出予測における実238

際的課題とその解決法, 水工学論文集, 第 36 巻, pp. 567– 572, 1992.
4) (財) 北海道河川防災研究センター・研究所 : 実践流出解
析ゼミ [講義テキスト編], 2006.
5) 佐山敬洋，立川康人，寶 馨，市川 温 : 広域分布型流出予測システムの開発とダム群治水効果の評価, 土木学会
論文集, No. 803 / II - 73, pp. 13 – 27, 2005.
6) 立川康人，佐山敬洋，宝 馨，松浦秀起，山崎友也，山路昭彦，道広有理 : 広域分布型物理水文モデルを用いた実時間流出予測システムの開発と淀川流域への適用, 自然
災害科学, 2007.
7) 日野幹雄 : 水文流出系予測へのカルマン・フィルター理論の適用, 土木学会論文報告集, 第 221 号, pp. 39 – 47,
1974.
8) Kalman, R. E. : A new approach to linear filtering andprediction problems, J. of Basic Engineering, Vol. 82,
pp. 35 – 45, 1960.
9) Wood, E. F. and Szollosi-Nagy, A. : An adaptive algorithm for analyzing short-term structual and parameter changes in hydrologic prediction models, Water
Resour. Res., Vol. 14, No. 4, pp. 577 – 581, 1978.
10) Kitandis, P. K. and Bras, R. L. : Real-time forecasting with a conceptual hydrologic model 1. analysis of
uncertainty, Water Resour. Res., Vol. 16, No. 6, pp.
1025 – 1033, 1980.
11) O’Connell, P. E. and Clarke, R. T. : Adaptive hydrological forecasting - a review, Hydrol. Sci. Bull., Vol.
26, No. 2, pp. 179 – 205, 1981.
12) Refsgaard, J. C. : Validation and intercomparison ofdifferent updating procedures for real-time forecasting, Nord. Hydro., Vol. 28, pp. 65 – 84, 1997.
13) 高棹琢馬，椎葉充晴 : 状態空間法による流出予測 - kinematic wave 法を中心として -, 京都大学防災研究所年報,
第 23 号, B-2, pp. 211 – 226, 1980.
14) 高棹琢馬，椎葉充晴，宝 馨 : 集中型流出モデルの構成と流出予測手法, 京都大学防災研究所年報, 第 25 号, B-2,
pp. 221 – 243, 1982.
15) 高棹琢馬，椎葉充晴，宝 馨 : 貯留モデルによる実時間
流出予測に関する基礎的研究, 京都大学防災研究所年報，
第 25 号, B-2, pp. 245 – 267, 1982.
16) 高棹琢馬，椎葉充晴，宝 馨 : 複合流域における洪水流
出の確率予測手法, 京都大学防災研究所年報, 第 26 号,
B-2, pp. 181 – 196, 1983.
17) 立川康人，椎葉充晴，市川 温 : 貯留関数法のモデルパラ
メータの不確定性を考慮した実用的な実時間予測手法, 水
文・水資源学会誌, Vol. 10, No. 6, pp. 617 – 626, 1997.
18) Evansen, G. : Sequential data assimilation with anonlinear quasi-geostrophic model using Monte Carlo
methods to forecast error statistics, J. Geophys. Res.,
Vol. 99, pp. 10,143 – 10,162, 1994.
19) Moradkhani, H. M., Hsu, K. L., Gupta, H. andSorooshian, S. : Uncertainty assessment of hydrologic
model states and parameters: sequential data assimilation using the particle filter, Water Resour. Res.,
Vol. 41, W05012, 2005.
20) Weerts, A. H. and El Serafy, G. Y. H. : Particle filtering and ensemble Kalman filtering for state updatingwith hydrological conceptual rainfall-runoff models,
Water Resour. Res., Vol. 42, W09403, 2006.
21) Romanowicz, R. J., Young, P. C. and Beven, K. J.
: Data assimilation and adaptive forecasting of waterlevels in the river Severn catchment, United Kingdom,
Water Resour. Res., Vol. 42, W06407, 2006.
22) Vrugt, J. A. and Robinson B. A. : Treatment of uncertainty using ensemble methods: comparison of se-

土木学会論文集Ｂ Vol.64 No.4, 226-239, 2008. 10
quential data assimilation and Bayesian model averaging, Water Resour. Res., Vol. 43, W01411, 2007.
23) Kim, S., Tachikawa, Y. and Takara, K. : Applying arecursive update algorithm to a distributed hydrologic
model, J. of Hydrol. Eng., ASCE, Vol. 12, No. 3, pp.
336 – 344, 2007.
24) Lindstrom, G., Johansson, B., Persson, M., Gardelin,
M. and Bergstrom, S. : Development and test of the
distributed HBV-96 hydrological model, J. of Hydrol.,Vol. 201, pp. 272 – 288, 1997.
25) Wohling, T. and Lennartz, F. and Zappa M. : Technical note: updating procedure for flood forecasting with conceptual HBV-type models, Hydrol. Earth
Syst. Sci., Vol. 10, pp. 783 – 788, 2006.
26) Seo, D. J., Koren, V. and Cajina, N. : Real-timeassimilation of radar-based precipitation data andstreamflow observations into a distributed hydrological model, Proc. of symposium HS03 held duringIUGG2003 at Sapporo, July 2003, IAHS Publ., no.
282, pp. 138 – 142, 2003.
27) Fread, D. L. and Ming, J. : Real-time dynamic floodrouting with NWS FLDWAV model using Kalman
Filter updating, Engineering Hydrology: Proc. of the
Symposium, San Francisco, July 1993, ASCE, pp. 946– 951, 1993.
28) 藤田 暁，大東秀光，上坂 薫，椎葉充晴，立川康人，市川温 : 分布型流出モデルに基づくダム流入量実時間予測モデルについて, 水工学論文集, 第 45 巻, pp. 115 – 120,
2001.
29) Shiiba, M., Laurenson, X. and Tachikawa, Y. : Realtime stage and discharge estimation by a stochasticdynamic flood routing model, Hydrol. Process., Vol.
14, pp. 481 – 495, 2000.

30) Friedland, B. : Treatment of bias in recursive filtering,IEEE Trans. on Auto. Contl., Vol. 14, No. 4, pp. 359– 367, 1969.
31) Dee, D.P. and Da Silva, A.M. : Data assimilation in
the precense of forecast bias, Q.J.R. Meteorol. Soc.,
Vol. 124, pp. 269 – 295, 1998.
32) 市川 温，村上将道，立川康人，椎葉充晴 : 流域地形の新たな数理表現形式に基づく流域流出系シミュレーショ
ンシステムの開発, 土木学会論文集, No. 691 / II - 57,
pp. 42 – 52, 2001.
33) 立川康人，永谷 言，寶 馨 : 飽和・不飽和流れの機構を導入した流量流積関係式の開発, 水工学論文集, 第 48 巻,
pp. 7 – 12, 2004.
34) Cunge, J.A. : On the subject of a flood propagation
computation method (Muskingum Method), J. of Hydraulic Res., Vol. 7, No. 2, pp. 205 – 230, 1969.
35) 陸 旻皎，小池俊雄，早川典生 : Multi-step, Multi-reachMuskingum-Cunge 法を用いた分布型水文モデルの開発,
水文・水資源学会，Vol. 12, No. 5, pp. 384 – 390, 1999.
36) Jazwinski, A. H. : Stochastic Processes and Filtering
Theory, Academic Press, New York, 1970.
37) Hsieh, C. and Chen, F. : Optimal solution of thetwo-stage Kalman estimator, IEEE Trans. on Auto.
Contl., Vol. 44, No. 1, pp. 194 – 199, 1999.
38) Ignagni, M. : Optimal and suboptimal separate-bias
Kalman estimators for a stochastic bias, IEEE Trans.
on Auto. Contl., Vol. 45, No. 3, pp. 547 – 551, 2000.
39) Ignagni, M. : An alternate derivation and extension ofFriedland’s two-stage Kalman estimator, IEEE Trans.
on Auto. Contl., Vol. 26, No. 3, pp. 746 – 750, 1981.
(2008. 2. 21 受付)

DATA ASSIMILATION OF A DISTRIBUTED RAINFALL-RUNOFF
PREDICTION SYSTEM BY KALMAN FILTER WITH BIAS CORRECTION
Takahiro SAYAMA, Yasuto TACHIKAWA and Kaoru TAKARAThis study proposes a data assimilation method for a distributed rainfall-runoff prediction system. The
system is composed by a rainfall-runoff model and a river routing model. Since it is computationallyinefficient for updating all the model variables on the real-time basis, the proposed filtering method takes
river discharges, which are simulated by Muskingum-Cunge method, as the state variables. It also sequentially estimates and collects prediction biases induced by the rainfall-runoff model. The application to the
Katsura river basin shows that the filtering with bias estimation and correction improves the accuracy offlood predictions.

239