水工学論文集,第56巻,2012年2月

群知能最適化手法を用いた分布型流出モデルのパラメーター同定
PARAMETER IDENTIFICATION OF DISTRIBUTED RUNOFF MODELUSING THE PARTICLE SWARM OPTIMIZATION METHOD

小槻峻司1・田中賢治2・小尻利治3・浜口俊雄4
Shunji KOTSUKI, Kenji TANAKA, Toshiharu KOJIRI and Toshio HAMAGUCHI
1学生会員

京都大学大学院 工学研究科都市社会工学専攻（〒615-8530 京都市西京区京都大学桂）2正会員 博（工） 京都大学 防災研究所 准教授（〒611-0011 宇治市五ヶ庄）
3正会員 博（工） 京都大学 防災研究所 教授（〒611-0011 宇治市五ヶ庄）
4正会員博（農）

京都大学

防災研究所

助教（〒611-0011 宇治市五ヶ庄）

In this paper, the particle swarm optimization (PSO) is applied into automatic parameter calibration
process of a distributed runoff model. As distributed runoff models require long simulation timecompared with general optimization problems, the number of particles and repeat computation times
should be selected property. We conducted sensitivity experiments for the number of particles and foundthat the PSO has to be applied in following conditions: i) to set the number of particles more than 100 in
the case of calibrating about five parameters, ii) to conduct repeat computations about 25 times. Analyzed
river discharge using identified parameters shows good agreement with the observed one.
Key Words : PSO, Parameter Identification, Distributed Runoff Modelている．多田 4 ）は，群知能最適化手法PSO（ParticleSwarm Optimization）6)を用いてパラメーター同定を行い，SCE-UAを用いた場合と比較して良好な結果を得ている．
本研究では，PSOを用いて分布型流出モデルのパラ現在，治水・利水などを目的に，広く降雨流出モデルメーター同定を行い，計算回数に重きを置いた実験と考が用いられている．その中で，物理的根拠は乏しいが基察を行う．多田4)は，集中型タンクモデルのパラメー底流出の特性を良好に再現する応答モデルとして，概念ター同定に最適化手法を適用したが，分布型モデルに適モデルも良く用いられている．概念モデルは，物理型モ用する場合，一回当たりの計算負荷を考慮する必要があデルと比較して計算が簡易であるがその反面，流出パラる．計算負荷が低いタンクモデルと異なり，膨大な計算メーターを決定するためのパラメーター同定が必要な点が欠点としてあげられる．概念モデルのパラメーターは， を実行不可能なためである．PSO手法では，粒子数や繰り返し計算回数が同定時間に直接影響するが，既往研究経験的な値やトライアンドエラー方式で決定される事がでは，その点については十分に議論されていない．
多い．そこで，客観的指標から自動的に同定する事を目
本論文ではPSOを用いたパラメーター同定手法を説明的として，パラメーターの同定問題を非線形の最適化問すると共に，その結果及び応用性についての議論を行う．
題に置き換え，既存の最適化手法を用いてパラメーターを推定することが広く行われている．
最適化手法には，多点探索によって局所解の問題に対２．パラメーター同定手法処できる大域的探索法がよく用いられている．例えば，遺伝的アルゴリズム（Genetic Algorithm： GA 1））を用いた研究や（田中丸2）など），Duanらによって提案され(1) 粒子群最適化手法の概要
パラメーター同定とは，未知係数を決定するための最たShuffled complex evolution method (SCE-UA法3）)を用い2）4）た研究（田中丸 ；多田 など），進化戦略（Evolution
小化問題（最適化問題）のことであり，キャリブレーションとも呼ばれる．一般に，モデルのパラメーター同Strategy：ES）を用いた研究（藤原5）など）が報告され
１． はじめに定は，その再現精度を目的関数として記述する事で，以下の最適化問題と定義できる．

m in im iz es .t .

f ( )

g j ( )  0 ,(1)j  1, ... q(2)ここに，f(θ)：目的関数, θ：パラメーターベクトル，g(θ)：制約条件であり，qは制約条件数である．
最適化問題を解く方法としては，古くは最急降下法から様々な手法が提案されてきている．本研究では，動物や昆虫，魚群の行動にヒントを得て確立された最適化手法である群知能最適化手法を用いる．群知能最適化手法
は，GAに代表される様な従来の最適化手法に比して収図-1 PSO解析手順
束速度が非常に速い点が利点として挙げられる．その代7)
表例としては，ACO（Ant Colony Optimization） やPSO
各粒子が持つ自己最良解を与えるパラメーターと，グ等が挙げられるが，本研究ではPSOを用いる．非連続のループ内で共有されるグループ最良解を与えるパラメー問題しか扱えないACO手法に対してPSOは連続的な問題
ターから，移動ベクトルを算出する．本研究では，坊原に適用可能なため，パラメーター同定に適している．
ら8)により提案されている，移動速度制約機能と局所ラPSOは，粒子群に於いてグループが発見したベストなンダム探索機能を追加して計算を行う．各パラメーター解を共有することにより最適解に収束させる手法である． の制約条件を飛び出した範囲の物については，ベクトル二次元空間，三次元空間の実現象から導出された手法で反射処理を施す事とする．
あるが，多次元空間に於いても適用可能である．
* t(3)
u e i  c 1 r a n d ( i   i )(2) 計算手順
PSOの計算手順を説明する．ここでは， n個の粒子群により，m個のパラメーター同定を行うとして説明を行う．PSOの基本構造は繰り返し計算であり，各期において，全粒子の解析，移動ベクトルの算出，パラメーターベクトルの更新といった手順をたどる(図-1)．以下に具体的な計算手順を説明する．
a) 初期値の決定
制約条件を満たすn個のパラメーターを，乱数を用いてm次元空間に発生させる．初期移動ベクトルについては，収束パラメーターに与える影響が小さい事を予備実験で確認したため，今回の解析では 0とする．
b) 目的関数値の算出
n個の粒子について解析を行い，目的関数を算出する．
なおこれ以後の，各期におけるパラメーター移動量を決定する移動ベクトルの算出時と，移動ベクトルを用いたパラメーターベクトル更新時の計算負荷は低いため，パラメーター同定時間はこの解析時間に依存する．
c) 自己最良解，グループ最良解の算出
n個の粒子それぞれについて，自己最良解(θ*)を算出する．自己最良解とは，それぞれの粒子が試行回数t 回の中で算出した目的関数値の中での最良解である．次に，n個の粒子が持つ各粒子の最良解からグループ最良解(θG*)を算出し，グループ内でその最良解及び最良解を与えるパラメーターベクトルを共有する．θ*は各粒子がそれぞれ記憶し，θG*はグループ内で共有される．
d) 移動ベクトルの算出

*tu g i  c 2 r a n d ( G   i )(4)
u r i  c3V r(5)t t  1 
vi  wvi  uei  ug i  ur itw h erev i  R m ax
 t 1i t t  i  vi

(6)
(7)
(8)

ここで，ue:各粒子の最良解に向かう移動ベクトル成分，ug:グループの最良解に向かう移動ベクトル成分，
ur:，局所ランダム探査ベクトル，Vr: 平均0分散1のn次元正規確率変数ベクトル v:パラメーター移動ベクトル，
Rmax: 探索範囲制約域，rand:(0:1)の乱数である．c, w ,Rmaxは，PSO手法のパラメーターである．PSOと同様に群知能最適化手法であるACOには，得られた目的関数値により各粒子の収束速度が上昇するスキームが組み込まれているが，PSOには組み込まれていない．PSOにおける収束速度は，c,wで支配される．
得られた次期パラメーターベクトルを b)目的関数値の算出 に返す事により繰り返し計算を行う．繰り返し計算を続けパラメーターベクトルを収束させる事で，パラメーターが同定される．
(3) 計算回数
PSOを用いたパラメーター同定を行う際に必要な計算

表-1 同定パラメーター

記号
意味
単位
同定範囲
-1
khbs
(0.0:1.0)
B層の横流出係数
kvbs-1(0.0:1.0)
B層縦流出係数khcs-1(0.0:1.0)
C層の横流出係数-1kvcs(0.0:1.0)
C層縦流出係数khds-1(0.0:1.0)
D層の横流出係数fi(0.3:0.9)
森林,畑地の流出率
ni(0.005:0.3)
森林,畑地の等価粗度 m-1/3 s
nm-1/3 s (0.0005:0.05)
河道の等価粗度
添え字i は，i=f,cであり，fは森林, cは畑地を表す．
時間は，1粒子の計算時間，粒子数，繰り返し計算数の積となる．分布型流出モデルに適用する場合， 集中型のモデルに比して，1粒子の計算時間が非常に大きいた
め，適当な解を同定するための粒子数と繰り返し計算数を設定することが重要である．特に粒子数は，少ないほど計算時間が少なくすむものの，局所解に陥るリスクが大きくなるため，パラメーター数に応じて妥当な粒子数を設定する必要がある．

提案するパラメーター同定手法を，分布型流出モデルのパラメーター同定に適用する．流出モデルはHydroBEAM9) (hydrological model for basin environmentassessment model)に用い，阿武隈川流域におけるモデルのパラメーター同定を行う．
(1) 適用モデルモデルの説明を行うが，紙面の制約から詳細な記述については参考文献に譲る．Hydro-BEAMは分布型の流出モデルであり，表面流出と河道流量はkinematic waveモデルにより計算され，基底流出は線形貯留法により記述される．地表面は森林，畑地，水田，水体，都市の5種の土地利用に区分され，各土地利用別に流出率や斜面の等価粗度が変更可能である．以下では同定パラメーター(表-1)に関係するモデルの数理式を記述する．
斜面系からの直接流出量は，kinematic wave法により記述される．
(9)

ここで，h:水深[m], q:単位幅流量[m2s-1]，r:降雨強度[ms-1], f:流出率であり，土地利用別に変更可能である．
河道内の河川流量も斜面系と同様に，kinematic wave法により記述される．

A Q
 q intx(10)Q A(11)mf

記号wc1c2c3
Rmax

値0.92.02.00.0001
0.5

図-2 1km解像度河道網データ

 sl mn  2 1 m2 1/ 3(12)

ここで，A: 流水面積[m2]，Q: 河川流量[m3 s-1]，qin:
単位幅横流入量[m2 s-1], sl: 河床勾配, n: 河道の等価粗度[m-1/3 s]，m: 河道の法面勾配である．
基底流を表現する線形貯留法は以下の式で記述される．

dS jdt Ij Oj

O j  (kv j  kh j )S j

３．計算モデルと計算条件

h q f  r ( x, t)tx

表-2 PSOパラメーター

(13)
(14)

ここに，S: タンクの貯留高[m], I: タンクへの流入量[m s-1], kh: 横流出係数[s-1]，kv: 縦流出係数[s-1]，sinθ: 斜面勾配，O: タンクからの流出量[m s-1]であり， 添字のjはタンクを表す．B,C,D層の3層のタンクを設置する．
(2) 入力データ
河道網データは，GDBD10)により250mで提供されている日本域の落水方向マップから，1km解像度にアップスケーリングを行うことにより作成した(図-2)．気象強制力データは地域気象観測システム（AMEDAS）及び地上気象観測のデータより作成し，蒸発散量は水文陸面過程モデルSiBUC11)により出力されたデータを用いる．
(3) PSOパラメーターと収束条件
PSOパラメーター(表-2)については，既往研究4) 6) 8)を参照して決定した．いずれも，推奨されている値や，既往研究で良好な結果が得られている値である．パラメー
ターの収束条件については，繰り返し計算5回の間にグループ最適解が更新されなかった場合に，パラメーターが収束していると判断する事とする．
(4) 目的関数
目的関数を変更して予備実験を行い，目的関数の違いを反映可能である事を確認した．同定流量観測点は阿久津，同定期間は2001年6月から10月であり，相対平均二乗誤差(以下CVRMSE)，相対平均誤差（以下MRE）を用いた．定義式は以下である．

1  t  Q o b s ( t )  Q s im ( t )  (15)2nQ
2C VRM SE obs

M RE 

Q o b s ( t )  Q s im ( t )1tn
Q obs (t )

(16)

1Psub
Pbest

0.8
0.6
0.4
0.2
0khb

khcParameter(a) 3パラメーター同定kvb

Normarized Parameter

Normarized Parameter

ここで，Qsim:解析流量[m3 s-1]，Qobs:観測流量[m3 s-1]で(a) 6月1日～10月31日あり，バーは対象期間の平均値を表す．MREとCVRMSEを比較した場合，誤差の二乗項を含むCVRMSEはピーク流量の再現性を重視し，MREは低水期の再現性を重視した指標となる．共に低い値が精度の高い解析結果となる．
この実験では，両目的関数に対して3パラメーター（khb, kvb, khc）の同定を行い，同定されたパラメーターによる流量解析値の比較を行った（図-3）．黒丸が観測(b) 8月1日～10月1日
流量であり，青線がMREに対して同定されたパラメー
図-3 目的関数変更実験による流量解析結果
ターによる解析流量，赤線がCVRMSEに対して同定され（青：MRE同定，赤：CVRMSE同定，黒丸：観測）た解析流量である．粒子数を100とし，それぞれ3組の初最適化問題としてみた場合，Pbestのパラメーターが解期値から収束させたため，MRE（青）,CVRMSE（赤）析精度の良い結果であり，Psubのパラメーターは精度が共に3本の線が引かれている．特に図-3(b)からは，悪い結果となる．2つの群に別れた理由は，Pbestの収束
CVRMSE同定による解析流量が，MRE同定による解析流パラメーターであるkh量に比して，ピーク流量の再現性が高い事が見て取れる．
b, kvb,の値が制約条件の同定範囲下限0に非常に近い事が原因であると考えられる．khb, kvb,
この予備実験より，PSOを用いたパラメーター同定によが0の近傍に最適解（Pbest）があるものの，Psubで同定り目的関数の違いを表現可能である事が確認された．こされた範囲に局所最適解（Psub）があり，Pbestの近くにれ以後，本論文中はCVRMSEを目的関数とする．
初期パラメーターベクトルが存在しなかった場合には，Psubに収束したと考えられる．これは，最適解が同定制
４．適用結果と考察約の上限や下限近くに存在し，かつ探索範囲の中央付近に局所解が存在する場合には，最適パラメーターを見逃す可能性があると言う事である．
(1) 収束パラメーター
タンクモデルのパラメーターとしてみた場合は，Pbest
同定パラメーター数が3つの場合（khb, kvb, khc）と, 6つのパラメーターはB層以下の基底流出が非常に少ないこの場合（khb, kvb, khc, ff, fc, n）の実験を行う(表-1)．粒子数ととを意味しており，非現実的である．この様に，一見非同定する観測地点と同定期間は，3-(4)目的関数変更実験
現実な解であっても，数学的に最適な解を見つけてしまと同じ条件である．初期パラメーターを10組発生させてう危険性もPSO手法は有していると言える．得られる情収束させた．同定された計10組のパラメーターのばらつ
報から物理的に妥当な範囲に探索範囲を設定する事が重きを図-4に示す．縦軸は同定されたパラメーターの同定要である．Pbestに収束したパラメーターとPsubに収束し範囲での正規化値を表し，同一の線でつながれたものがたパラメーターによる流量解析結果を図-5に示す．2001
1組のパラメーターベクトルである．赤は収束パラメー
年は同定に用いた年であるが，同定されたパラメーターターによる目的関数CVRMSEが0.5以上（以下Psub収束のを2002年，2003年に適用した場合も流量が再現されていパラメーター），青は0.5未満のパラメーター（以下る事が分かる．
Pbest収束のパラメーター）である．図-4より，タンクの図-4(b)での同定されたパラメーターのばらつきを見
同定パラメーター（khb, kvb, khc）には，大きく2つの群がると，他のパラメーターに比して，（ff, fc,）のばらつきあるのが分かる．
1

Psub
Pbest

0.8
0.6
0.4
0.2
0khb

khc

kvb

fffc

nParameter(b) 5パラメーター同定

図-4 パラメーターの収束値 （赤は目的関数CVRMSEが0.5以上，青はCVRMSEが0.5未満のパラメーターベクトル）(2) 同定に必要な計算回数
タンクモデルに関する同定パラメーターが3つの場合（khb, kvb, khc）と5つの場合（khb, kvb, khc, kvc, khd），畑地・森林・河川の流出率と粗度に関する同定パラメーターが3つの場合（n, ff, fc,）と5つの場合（n, ff, fc,,nf,nc）の4実験を行う(表-1)．同定対象パラメーターを変えることで，模擬的に目的関数への感度が異なる2つの流出モデルのパラメーター同定を想定した．同定地点を白河，期間を2001年8月から9月とする．一回の計算負荷を減らし，4(1)より多くの粒子数での実験を行うためである．粒子数を10, 25, 50, 100, 250, 500とし，それぞれ初期パラメーターを4組発生させて収束させた．
タンクモデル，5つのパラメーターに関して同定された計24組のパラメーターのばらつきを図-6に示す．図の見方は，図-4と同様である．ただし，粒子数が50以下の12組と，100以上の12組を別々に図示し，収束しているパラメーター（赤線：Pcov）からばらついているものを青線（Pdiv）で図示している．粒子数100以上の場合に比して，粒子数50以下の収束パラメーターのばらつきが大きい事が分かる．寄与度が同程度の5つまでのパラメーターを同定する場合，粒子数を100以上に設定すれば，ほぼばらつきなく同定可能であると言える．しかし，粒子数100以上の場合でも青線のベクトル（Pdiv）の様に，異なる解に収束する可能性があるため，複数の計算を行うことが必要であろう．一方，タンクモデルの3つのパラメーターに対して同定を行った場合は，ほぼばらつきなく同定された．畑地・森林・河川の流出率と粗度に関して同定を行った場合にも同様の結果が得られた．
各粒子実験の際の繰り返し計算数を示す（図-7）．各粒子数4組の実験の平均値と，バーは最大値・最小値を示している．図-7からは，各粒子実験において収束に必要な繰り返し回数が概ね変化しないことが分かる．
ところで，最適化において最も重要なのは局所解に陥らないことであり，局所解に陥るリスクが少ない，最低限の粒子数を見積もる必要がある．局所解を「同定されたパラメーターによる目的関数値が全24回の実験平均値(a) 2001年，CVRMSE（Psub：0.63417，Pbest：0.32377）(b) 2002年，CVRMSE（Psub：0.52966，Pbest：0.35616）(c) 2003年，CVRMSE（Psub：0.49298，Pbest：0.40642）
図-5 同定されたパラメーターによる流量解析結果（同定：2001年，検証：2002，2003年）

が大きい事が分かる．この様にばらつきの大きいパラメーターには，二つの原因があると考えられる．一点目
は，目的関数値に与える影響が低いため，6次元空間に於いて影響が高いパラメーター（khb, kvb, kvc,n）が優先して同定されたためだと考えられる．二点目は，水収支のみを考えた場合，森林・畑地の流出率の大小によりバランスが可能であるためだと考えられる．厳密には森林・畑地の粗度が異なるため流出にかかる時間が異なるが，粗度の違いが結果に与える影響は小さいため，それ以外のパラメーターも同時に同定した場合にばらつくと考えられる．（ff, fc,）の様に，目的関数への寄与度が低い，若しくは水収支のバランスが可能なパラメーターについては，再度同定する事が良いと考えられる．実際に（khb, kvb, kvc,n）を固定し，（ff, fc,）のみを収束させる追試実験を行い，両パラメーターがばらつきなく同定される事を確認している．
0.8Pdiv
Pcov

0.60.40.20khb

khc

kvb
Parameter

khd

kvc

Normarized Parameter

Normarized Parameter

RMSE MORE THAN 100BIRD
110.8

Pcov
Pdiv

0.6
0.4
0.20khb

khc

kvb
Parameter

khd

kvc

3param5param102550100250

Number of Particle(a) タンクパラメーター同定

500

Number of Repeat

Number of Repeat

図-6 5つのタンクパラメーターの収束値 （左：粒子数10,25,50の収束値，右：粒子数100,250,500の収束値）
3025201510503025201510503param
5param102550100250

Number of Particle(b) 流出率・粗度パラメーター同定

図-7 収束に必要な繰り返し計算数 （赤線：3パラメーター同定，青線：5パラメーター同定）500

表-3 粒子数別の局所解に陥る確率[%]

粒子数102550100250500

タンク3param. 5 param.
507575502525002525250

流出率・粗度
3 param. 5 param.
5025050025000000

高いタンクパラメーターの同定に要する計算回数を最大値として見積もる必要がある．本研究から，5つ程度のパラメーターの同定に，粒子数100以上，繰り返し計算数25回程度が必要との指標を得た．ただし，局所解に収束する可能性があり，複数の計算を行う事が必要である．
本研究では，観測値を再現する様なパラメーターの同定を行ったが，パラメーター同定は観測情報と，同定パラメーターの関係法則を見つける上でも重要な作業である．例えば，今回同定したタンクのパラメーターは，土壌物性に支配されるはずである．観測点のある場所で複数の実験を行い，グローバルプロダクトから得られる土壌物理データと同定パラメーターとの関係を見つけ，水文貧観測地域に適用する，といった応用を検討している．

+σ/2より大きい場合」と定義して，各粒子数の局所解に陥る確率を算出した（表-3）．σ2は分散である．タンクパラメーター同定における局所解に陥る可能性が高いのは，パラメーター感度が高く，多くの局所解が存在するためだと考えられる．図-7において，タンクパラメーター同定に要する繰り返し計算数が，流出率・粗度
参考文献
パラメーター同定に要する計算数より多かったのも，局
1) Melanie Mitchell, Stephanie Forrest and J. H. Holland. The royal
所解から抜けるのに数回の試行が必要であった可能性が
road for genetic algorithms, Proceedings of the First Europeanある．一方で，同じ粒子数では同定パラメーター数が多Conference on Artificial Life, Vol1001-1, pp.245-254, 1992.
いほど局所解に陥る可能性が高いと考えられるが，タン2) 田中丸治哉. タンクモデル定数の大域的探索，農業土木学
クパラメーター同定では異なっている．この点は，実験会論文集，Vol.178，pp103-112, 1995.
数を多くして確認する必要がある．
3) Duan, Q. Sorooshian, S. and Gupta V.K. Effective and efficient
以上の議論を踏まえて，PSO手法によるパラメーター
global optimization for conceptual rainfall-runoff models, Water
同定を他の分布型流出モデルに適用する場合を考える．
Resources Letters, Vol.28, pp1015-1031, 1992.
その際には，一般に感度が最も高いと考えられる，タン
4) 多田毅．PSO アルゴリズムによる流出モデルパラメータ
クパラメーターの同定に要する計算回数を最大値としての最適化, 水文・水資源学会誌, Vok.20-5, pp.450-461, 2007.
見積もれば良いといえる．本実験では，タンクパラメー
5) 藤原洋一，田中丸治哉，畑武志，多田明夫．進化戦略にターの同定には粒子数100以上，繰り返し計算数が25回
よる流出モデル定数の最適同定，農業土木学会論文集，程度であり，総試行回数は2500回という結果が得られた．
Vol.227，pp.119-129, 2003.
一回の試行時間を3分とすると，2500回の計算に要する
6) James Kennedy and Russsell Eberhart. Particle Swarm
時間は5日程度であり，十分に実行可能である．タンク
Optimization, Proceedings of IEEE International Conference on
パラメーターの感度はモデル・流域・同定期間によってNeural Networks, pp1942-1948, 1995.
異なるが，オーダーとしては変わらないと考えられるた
7) Marco Dorigo and Gianni Di Caro. Ant Algorithm for Discreteめ，今回の結果は一つの指標となりうる．また，粗度やOptimization, Artificial Life, Vol.5-2, pp.137-172, 1999.
流出率など，感度の低いパラメーターの同定では，要す
8) 坊原尚記，作田健，香月智，嶋丈示．PSO を用いた鋼製る繰り返し計算数が少なくなると考えられる．
防えん堤の最適断面設計，土木学会論文集，A-62, pp.681692, 2006.

５．結論

9)

小尻利治，東海明宏，木内陽一. シュミレーションモデル
での流域環境評価手順の開発，京都大学防災研究所年報，

本研究では粒子群最適化手法PSOを用いた分布型流出モデルのパラメーター同定を行った．得られた知見をまとめると以下の様になる．
・PSOを用いたパラメーター手法の定式化を行った．目的関数変更実験では，目的関数に応じたパラメーター同定が可能である事を示した．
・目的関数への寄与度の高いパラメーターが優先されて同定される事が分かった．そのため，寄与度の低いパラメーターは，再度同定計算を行う方が望ましい．
・PSOによるパラメーター同定を行う場合，感度の最も第 41 号Ｂ-2，pp119-134,1998.
10) Yuji Masutomi, Yusuke Inui, Kiyoshi Takahashi and Yuzuru
Matsuoka. Development of highly accurate global polygonaldrainage basin data, Hydrological Processes, Vol:7186-23,pp.572-584, 2009.
11) Kenji Tanaka. Development of the new land surface scheme
SiBUC commonly applicable to basin water management andnumerical weather prediction model, doctoral dissertation, KyotoUniversity, 2004.
（2011.9.30受付）