1.はじめに
計画規模に匹敵する，あるいはそれを上回る極めて大きな洪水がしばしば発生し、特に県管理の中小河川流域において多くの被害が出ている．こうした河川流域は今後とも高い治水整備を実現することは容易ではなく，被害を軽減するためには，確度の高い洪水予測情報を得て，それに基づく的確な避難体制を構築することが洪水災害を軽減する基本と考えられる.中小河川流域を洪水予測の対象とする場合，観測情報の蓄積が不十分なことが多く，流出予測モデルを構築することが困難であることが多い．この困難を克服するためには，特定の中小河川を洪水予測の対象とするのではなく，水系全体を対象とする実時間流出予測システムを構築する必要がある．また，複数観測地点で得られる観測値を用いてモデル状態量をフィルタリングにより同化し、それを時々刻々，予測の初期値として用いれば，観測情報のない地点においても観測による初期状態量更新の効果が及ぶので，観測情報のない地点を含めて予測精度を向上させることができる.従来，集中型流出モデルとカルマンフィルタとを組み合わせた流出予測手法が開発され，実用的にも用いられてきた．佐山らﾘはこれを発展させ，分布型流出予測システムにカルマンフィルターを導入した実時間流出予測手法を提案した．分布型流出モデルにカルマンフィルタを適用する場合，実時間計算を可能とするために効率的な状態量の取り方と計算アルゴリズムの開発が重要となる. そのためこの手法ではDee and Da Silva²〉の手法をもとに，流域斜面の流出モデルに起因する誤差バイアスを，状態量とともに時々刻々フィルタリングを行い推定しながら，状態量の予測計算を行うアルゴリズムを提案した．しかし，この手法は洪水追跡モデルの状態方程式としてマスキンガム-クンジ法による線形式を用いており，状態方程式の線形化の扱いに不十分な点があった．そこで本研究では，佐山の提案した手法で用いた状態方程式の線形化の方法を見直し，より精度の高いアルゴリズムを提示するとともに，実流域に適用してその有効性を確認する.
2.カルマンフィルタの洪水追跡モデルへの導入
非線形の微分方程式を離散時間の線形方程式に変換したのち、カルマンフィルタ理論を適用する．本研究では，洪水追跡モデルとして，キネマティックウェーブモデルの基礎式を逐次線形化した方程式を新たに用いた．以下にその基礎式を導出する． (1)キネマティックウェーブモデルの線形化本研究で用いる広域分布型流出予測システム³は，落水線型の分布型斜面流出モデルと河道網の洪水追跡モデルで構成される．淀川の流出予測システムでは河道の合流点および約3kmごとの河道地点で河道網を分割しており，その各河道区分に対して，洪水追跡モデルを適用している.時刻をt，計算区間上端からの距離をxとして，断面平均流量をQ，通水断面積をA，単位幅あたりの側方流入量をqとすると，連続式は
となる．一方，運動方程式は，抵抗則としてManning公式を用いると, Q = αAmと書くことができる. αは，河道断面が広幅矩形であると想定した場合，mを無次元の流量流積パラメタ(= 5/3), Iを河道勾配，NをManningの粗度係数，Bを河道幅としてα = √I/(NBm-¹)とし，洪水期間中は変化しないとして一定の値を与える.流れの伝播速度をc = f(Q)とすると,
となり，連続式(1)は
と変形することができる。左辺各項を，空間的に(1–X) : X，時間的に1：1の重みを付けて次のように差分近似する.
とする．非線形関数g (Q',)をあらかじめ既知なQ', = Qoのまわりに，Taylor級数展開を用いて線形近似すると，次式が得られる.
とすると，式(2)は以下のようになる.
うことができる．また，Qoは直前に推定された流量を用いて, Qo
としている. マスキンガム-クンジ法との違いは，各係数C 1からC 5における流れの伝播速度に，時刻th-1とt,でそれぞれ得られる伝播速度を用いていることである。なお，流量や勾配などの条件に応じてΔtとΔxを自動的に調整する陸ら4のMulti-step,  Multi-reach法を用い、数値計算の安定性を確保している.
(2) 洪水追跡モデルの河道網全体への適用
これまでは一つの河道区分を対象にした場合の洪水追跡計算法を述べた．次に複数の河道区分が存在する河道網を対象にした場合に，どのように式(5)を連立させ，河道網全体の流量を計算するかについて説明する.とこでは簡単のため、図-1に示す河道網を対象にして以下の議論を進める．図-1は5つの河道区分が河道網を構成している例である．それぞれの河道区分に式(5)を適用すると，以下のようになる．
式の簡略化のため，測方流入量は時刻t,時のもののみを記している. なお，Cの肩につけた数字は，その番号の河道区分を対象にして求めた係数であることを意味している．上式を整理すると，以下のように行列表示できる.
以下行列計算の説明を簡単にするため，式(6)の左辺の係数行列をXn，右辺第一項，二項の係数行列をそれぞれYn,  Znとすると，式(6)は次式で表される.
式(6)に示したように，Q';を下流から上流に向かって順に並べて表示すれば，X,は上三角行列になるので，逆行列X,1を比較的容易に求めることができる．その逆行列を両辺にかけた式をベクトル表示すれば，
である．式(7)を状態量の伝達方程式として，次章で述ベるカルマンフィルタ理論の一連の式に適用する.
3.バイアス補正カルマンフィルタ
洪水追跡モデルの状態量をフィルタリングの対象とする手法は，空間分布する状態量を複数地点の観測流量で同化できる点で優れている．しかし，予測のリードタイムが長くなると，洪水追跡の予測結果が斜面流出の予測結果に依存するため、河川の状態量を更新した結果が消えてしまうという欠点がある。つまり，斜面流出の予測モデルに起因する系統的な誤差をフィルタリングで予測・補正しなければ，河道網の状態量をフィルタリングする効果は十分に発揮できない.従って，河道および斜面の状態量を同時に観測更新することが重要となる．しかし，水系全体を対象とするような洪水予測モデルを考える場合，これらの状態量すべてを観測更新の対象とすることは，実行可能性の点で困難さが付きまとう．そこで本研究では，斜面流出モデルが算定する流出量の予測誤差を，洪水追跡モデルのシステム方程式における予測バイアスとして捉え，逐次その予測バイアスを推定・補正する方法を用いた．以下にその導出方法を簡単に示す．
(1)問題の定式化
カルマンフィルタ理論で扱うシステム方程式は，前節の式(7)に，モデルの不確かさをあらわすシステムノイズ項を加えて,
のように表すことができる．対象流域の河道区分数をpとすると,  Qn-1とqn,  ξnは(p × 1)次元のベクトルとなり、AnとDnは(p × p)次元の係数行列となる.一方，η,を観測ノイズとして，観測方程式を以下のように装す.
データ同化に用いる観測流量がr個，すなわち観測地点がr地点ある場合，観測ベクトルQnobsは， (r × 1)の次元になる．また，H nは，(r × p)次元の係数行列であり，状態空間を観測空間に線形写像する役割を担っている.バイアス補正カルマンフィルタでは，斜面流出モデルで得られる斜面からの流出量に誤差が発生すると考え，河道における連続式(1)は単位幅当たりの未知のバイアス項βを導入し、以下のように表される.
となる．ここで,
親測方程式は以下のように表される.
ここで，B,は(p × p)次元の係数行列である．基礎式の一般性を持たせるため、観測方程式にもC n βの項を加えている.
(2)バイアス補正カルマンフィルタ
式(11)と式(12)から状態量とともに斜面流出モデルに起因するバイアスも逐次推定する適応カルマンフィルタかは，計算の次元が拡大するため，計算負荷が大きくなり，リアルタイム予測に用いるには難点がある.そこでDee and Da Silva²)は，行列計算の次元を状態量の数の次元に抑え，バイアスと状態量を順次求める方法を提案した．彼らは，時刻ｎにおけるバイアス補正前後の流量の推定値をそれぞれŽņ,  Ôn, バイアスをb,とすると，流量のバイアス補正を，

4. 流出計算のアルゴリズム
バイアス補正カルマンフィルタを用いて流出計算を行う場合の計算アルゴリズムを述べる。現在時刻はt,,で，時刻tn-1のフィルタリング後にバイアスを補正した状態量Q,,-1とバイアスδ,-1が得られている状況を想定する、本研究では、計算の時間ステップを１時間とし，３時間先までの予測流量を得ることを目標とする.図-2に実時間流出計算アルゴリズムの概略図を示す.なお、バイアスは予測計算では変化しないため、
が成り立つ. ここで，記号’~’は，個数分だけ前の時刻の観測流量をもとに推定された状態量であることを表している．以下に，バイアス補正カルマンフィルタを適用した場合の流出計算のアルゴリズムを示す. 1)時刻tn-1からtnの間の現況雨量データを入手し,現在時刻t,まで現況流出計算を実行し，状態量Ž,を得る. 2)時刻tn-1で求めたバイアスδn-1(= δn)をもとに,バイアスを補正した状態量Qを計算する. 3)時刻t,の観測流量Ωn,obsを入手できたら，カルマンフィルタによりバイアスδ 。の更新計算を行い, b,を得る. 4) 更新されたバイアスδηをもとにQ をバイアス補正するとともに，観測流量Ωn,obsを用いてカルマンフィルタにより更新し、Q,を求める. 5)状態量がフィルタリングにより更新されると，以後予測計算を行う．更新された状態量Q,を初期値とし、1時間先までの予測降雨データを用い、Õ n +1を求める. 予測計算によって得られたÕ n +1にバイアスδπ+1を用いて補正された流量Q n+1を1時間先の予測流量とする. 6)  2時間先，3時間先の予測流量Ô¬¬¬  ô¬¬ も同様に，予測計算によって得られた流量をバイアス補正した値とする. 7)  l時間時刻が経過し現在時刻がtn+1となったら、一連の計算を同様に行い，予測流量を得る．
5.桂川流域への適用
淀川流域のうち、桂川流域(流域面積1100km² :図-3)の桂地点上流域(流域面積833km²)を対象に，本手法を適用して流出予測計算を行った．計算対象洪水として， 1992年8月18日から8月22日までの洪水イベントを用いた．この時期には，上流域に主要なダムが存在しないため，自然流況で予測性能を検証できる.図-4に，バイアス補正を行った場合と行わない場合の，桂地点における3時間先予測計算結果を示す．カルマンフィルタ理論を適用していないオフラインでのシミュレーション結果と観測流量も併せて示している.システムノイズξπ、観測ノイズη,の分散はともに100 m6/s²としている．バイアス補正を導入する場合，バイアス補正に関するパラメータγ（状態量の誤差分散から定めるバイアス補正に関する誤差分散の寄与率）は0.25としたﾘ. このパラメタは0から1の値をとり，1に近いほど観測値の情報は主にバイアスの更新に用いることになる．一方，０に近ければその情報は主に状態量の更新に用いることになる．なお，γが０のときはバイアスは更新されないので，バイアス補正を考慮しない通常のカルマンフィルタによる結果が得られる．図-4より，バイアス補正を行わない場合，オフラインでのシミュレーション結果に近づき，フィルタリングの効果が十分に反映されていないことがわかる．この傾向は，予測のリードタイムが長くなるにつれ，顕著に現れることがわかった．一方バイアス補正を行った場合では，斜面流出に起因する誤差バイアスを逐次推定して河川流量の予測を補正しているため，リードタイムが長くなってもフィルタリングの効果が現れている.次に，更新計算に用いる観測流量の個数の違いが及ぼす効果を図-5,  6,  7,  8に示す. 図-5には，桂地点のみ，桂地点と亀岡地点，および亀岡地点のみの観測流量をそれぞれ更新計算に用いた場合の，天竜寺地点における3時間先予測計算結果を示す. 図-3に示すように,桂地点の約5km上流に天竜寺地点が，さらに約10km上流に亀岡地点が位置している．図-5から天竜寺地点では， 1地点のみの観測流量を更新計算に用いるより、桂・亀岡両地点の観測流量を用いた方が予測の精度は上がっていることがわかる.図-6, 7, 8には，ハイドログラフの立ち上がり始めの８月１９日０時，立ち上がる途中の同日９時，ピーク時の同日18時におけるバイアス補正量(式(13)のbn)を示している．これらの図はそれぞれ更新計算に桂地点のみ、亀岡地点のみ，桂・亀岡両地点での観測流量を用いて計算されたバイアスを示している．図-6では，桂地点に近い下流で集中的に比較的大きいバイアス補正がなされ，上流域ではバイアス補正による効果があまり及んでいないことがわかる．図-7では，亀岡地点より上流でバイアス補正の効果が表れ，亀岡地点より下流ではバイアス補正の効果があまり及んでいない．しかし，桂・亀岡両地点で得られた観測流量を更新計算に用いて計算を行った図-8では，下流から中流にまで比較的小さいバイアス補正がなされ，より広い地点にバイアス補正の効果が及んでいることがわかる．つまり，複数の観測地点から得られる観測情報を更新計算に用い予測計算後にバイアス補正を行うことで、それ以外の観測情報のない地点にまで更新計算の効果が及び状態量が修正され予測精度が向上することになる.
6.結論
本研究では，広域分布型流出予測システムの洪水追跡モデルにキネマティックウェーブモデルの基礎式を線形化した方程式を新たに用いることで，計算の精度を向上させるとともに，理論的にも確かな方法でカルマンフィルタを用いた予測計算を行うことが可能になった．バイアス補正が流出計算に及ぼす効果については、状態量更新計算に関しては，バイアスを考慮していない場合と比べてあまり違いが見てとれなかったが，予測計算に関しては，バイアスを考慮した場合のほうが、考慮していない場合より明らかに観測流量に近づき，予測精度が向上する結果が得られた．また，観測流量が得られない地点では，他の上流及び下流に位置する観測地の観測流量を用いて状態量の更新を行うと，予測精度が向上するという結果も得られた．この結果から，観測情報の警積が不十分な中小河川流域においても洪水予測の高精度化が見込まれることがわかった.謝辞:本研究を進めるにあたり，(財)河川環境管理財団の河川整備基金助成事業による補助をいただいた．記して謝意を表する. (a) 8/19 0:00 (b) 8/19 9:00図-6 桂地点のみの観測流量を更新計算に用いた場合のバイアスの変動図-7 亀岡地点のみの観測流量を更新計算に用いた場合のバイアスの変動図-8 桂・亀岡地点の観測流量を更新計算に用いた場合のバイアスの変動