XバンドMPレーダ情報を用いた大気場の推定と降雨流出過程の不確実性
A study on estimation of atmospheric fields with X-band MP radarand uncertainty of rainfall-runoff process
山田 朋人1
Tomohito YAMADA1

北海道大学大学院 工学研究院 准教授

要旨

本研究では，複数のXバンドMPレーダによる観測データに対して変分法の一種であるMUSCAT法を適用することによって豪雨時の3次元風速場を推定した．これによって，2014年9月に道央地域を襲った線状降水帯の大気場の特徴を明らかにした．
降雨流出過程に関しては，確率微分方程式による定式化を行い，同過程に内在する不確実性の主要な要素として考えられる降雨強度および初期土壌水分量がピーク流量に与える不確実性を議論可能とし，複数の山地流域を対象に不確実性の大きさを明らかにした．
《キーワード：XバンドMPレーダ；変分法；風速場；降雨流出過程；不確実性；
Fokker-Planck方程式；確率微分方程式》― 73 ―

1. はじめに近年，北海道では線状降水帯による集中豪雨の発生が増加傾向にある（Yamada et al. 2012）．線状降水帯とは複数の積乱雲により組織化された線状に伸びる降水帯のことであり，その規模は長さ100km以上に及ぶ．数時間に渡り同じ場所に留まるその定在性は局所的な集中豪雨をもたらし大規模な災害に繋がる．また，その形成機構はバックビルディング現象と呼ばれる．バックビルディング現象について説明する．i）雨を降らす積乱雲下の地表面には降雨による冷気外出流が生じ周囲が冷やされる．ii）そこへ暖湿塊が流入するとその空気は冷たい空気とぶつかり，上空へ持ち上げられ新たな積乱雲が形成される．iii）成長した積乱雲は風下へ移動するが，暖湿流の供給が継続することで新たな積乱雲が連続的に形成され続け，複数の積乱雲が線状をなすことで線状降水帯が形成される（小倉 1999）．線状降水帯は前線や温帯低気圧等の大規模な気象擾乱に伴い発生する場合もあるが，それらから離れた地域でも発生することがあるため，その発生位置の予測は困難である（吉崎, 加藤 2007）．
北海道における線状降水帯の代表事例として2014年9月11日に発生した北海道道央圏の大雨が挙げられる．この事例は11日午後19時頃に室蘭市付近から雨雲が発生し始め，11日午前0時頃から同日の朝にかけて道央およびその周辺に強い雨が降り，総雨量が300mmを超える豪雨となった．この事例によって北海道には初の大雨特別警報が発表され，札幌市は道路冠水98件，土砂崩れ9件，避難者数479名であったと報告している（札幌市の報告書より）．
9月11日午前0時の北海道の大気場として，対流圏中層における気温と風速ベクトルを図1，対流圏下層における気温と水蒸気フラックスを図2に示す．対流圏中層では-16℃以下の寒気が流れ込んでおり，対流圏下層では南東から9℃以上の暖かい湿った空気が流入していたことが確認できる．したがって，大気は不安定な状態であったことがわかる．また，北海道道央において降雨が強かった9月11日午前0時から6時の積算降雨量を図3に示す．積算雨量は多い場所で200mmを超え，雨域が線状に伸びていることがわかる．

2. 使用データ
本研究では気象庁の札幌レーダと新千歳空港レーダの2つのドップラーレーダによる観測値を用い，線状降水帯の雨域や風の場の解析を行う．2つのレーダの位置を図4に示す．ドップラーレーダはレーダの送信周波数と受信周波数の差である周波数偏移を検出することで，受信信号の位相の変化を検出することができ，ドップラー効果を利用して測定対象の動きを推定するものである（深尾 2004）．
取得した観測データは反射強度とドップラー速度である．反射強度とは雲の中に含まれる降水粒子から反射して戻ってくる電波の強さのことであり，降水粒子の数が多いほど，また粒子が大きいほど強くなる（単位は［dBZ］）．また，反射強度は降水強度として捉えることができ，経験式を用いることで降水強度を推定することが可能となる．ドップラー速度とは降水粒子のレーダビーム方向の速度であり，レーダビームの送信周波数と降水粒子から反射された受信周波数の違いから算出される（単位は［m/s］）（気象庁資料より）．

― 74 ―

図1 : 2014年9月11日午前0時（JST）の対流圏中層における気温［℃］および風速ベクトル［m/s］

図2 : 2014年9月11日午前0時（JST）の対流圏下層における気温［℃］および水蒸気フラックス［kg/kg m/s］

― 75 ―

図3 : 2014年9月11日の午前0時から午前6時における積算降雨量［mm］（使用データはXRAIN）

図4 : 解析に使用したレーダの位置図（Google Earthより）― 76 ―

3. 次元風速場の推定法（MUSCAT法）
2台のドップラーレーダで別々に同一の3次元空間を観測し，風速3成分（u，v，w）を得る方式をデュアルドップラー法という（深尾, 2004）．本研究ではその中でも変分法を用いた3次元風速場を推定する手法であるMUSCAT法を用いる（Bousquet and Chong，1998）．変分法とは汎関数（関数を含む関数）が最小となるような最適な関数の形を求める方法である．最も簡単な例として2点間の最短経路を考える問題がある．何も制約条件を与えなければ，単純に直線が解を与えるが，制約条件を与えれば解は単純ではない上に複数存在し得る．こういった制約条件下における関数形の最適化問題を変分法という．
MUSCAT法は(1)式（下にまとめて記述）の汎関数 F が最小となるような風速3成分（u，v，w）を各高度平面において求める方法である．
u, v, wは最終的に疎行列方程式となり共役勾配法によって解くことで定められる．最小化にあたっては，(2)式を満たす必要がある．(3)式，(4)式，(5)式はそれぞれ汎関数 F における関数A，B，Cである．
Aはデータフィット項といい，解析された風速3成分（u，v，w）から計算されるドップラー速度と観測されたドップラー速度の差を表す．Bは連続式による束縛条件を意味する．Cはフィルター項といい，解析される風速3成分（u，v，w）の2次微分とその微分の二乗和を表しており，微小な風の変動を抑える平滑化項である（Bousquet, O. and M. Chong, 1998; Yamada, T. 2013）．

F (u, v, w) = ∫ [ A(u, v, w) + B(u, v, w) + C (u, v, w)]dxdys∂F= 0,∂u

∂F= 0,∂vn∂F=0∂w（1）（2）n ( p)1 p qA(u, v, w) = ∑ ∑ ωq [α q u + β q v + γ q ( w + vT ) − Vq ]2
N p =1 q =1
 ∂u ∂v 1 ∂ρw B(u, v, w) = µ1  + + ∂x ∂y ρ ∂z 

（3）

2（4）

C (u, v, w) = µ 2 [ J 2 (u ) + J 2' (u ) + J 2 ( v ) + J 2' ( v )+ J 2 ( w) + J 2' ( w)]（5）各パラメータの詳細は以下の通りである．u，v，w : x, y, z方向の風速成分［m/s］，p : レーダ番号，q :p番のレーダにおける観測番号，N : 全領域におけるqの総和，ωq : cressmanの重み係数，αq，βq，γq :方向余弦，νT : 雨粒の終端落下速度［m/s］，Vq : 観測された動径方向成分［m/s］，μ1，μ2 : 正規化重み付きパラメータ，ρ: 大気の密度［kg/m3］，J2 : 微分演算子
MUSCAT法では，ある高度面で求まる風速3成分は1つ上の高度の下面境界条件として使用される．
まず，下面の境界条件としてw=0［m/s］を与え，風速分布を最低高度から順々に各高度で求める．
MUSCAT法では各高度で水平風速u，vの誤差を可能な限り小さくすることで高高度における鉛直風wの誤差が大きくなることを防いでいる（清水ら 2007）．この手法では互いのレーダのビームが重なる範囲でのみ風の場が算出される．

― 77 ―

4. 確率情報付きの流出量推定手法
降雨流出過程では得られる流出量の精度は流域固有のパラメータのほか，入力情報である降雨強度や対象とする降雨直前の流域貯水量に依存する．しかしながら，これらの入力情報は本来不確実性を有するものである．降雨強度を算定するために使用される地上雨量計に着目すると，我が国では全国に約1300ヶ所展開されている気象庁の地域気象観測システム（通称アメダス）や国土交通省水文水質データベースなど諸外国と比較すると空間的に密な観測網を有するものの，これらは約数km～十数kmごとの分布である．そのため地上観測による不確実性は存在する．一方，レーダは時空間的に詳細な降雨情報が得られるものの，それは上空での値であり，またレーダからの距離によっても精度は大きく依存する．
降雨流出過程に着目すれば，たとえ完全に雨量を観測できたとしても，ピーク流量に寄与する雨量を推定できなければ流量の予測には不確実が生じる．既存の降雨流出過程は決定論的手法が採用されているが，降雨分布の不確実性（もしくは不確定性）を考慮した場合の流量の不確実性を定量化することは将来予測結果に含まれる不確実性の相対的な大きさを議論する上で不可欠であろう．そこで本研究では既存の降雨流出過程に確率微分方程式を適用し，入力値である降雨強度の不確実性に伴う流量の不確実性を明らかにすることを目的とする．一連の研究成果は吉見ら（2015）やSupraba and Yamada（2015）において詳細な説明がなされており，ここではその概要を紹介する．
本研究で使用する流出計算手法は山田（2003）に準拠している．同研究は単一斜面における一般化した降雨流出の基礎式を提案しており，その概要は以下の通りである．単一斜面に対して幅広矩形断面を想定すると，連続式は式(6)として表される．式(7)に示すように運動則を断面平均流速は水深の冪乗に比例するという形式で表現することで多様な流出形態に対応可能である．一般的なKinematic Wave理論
における式(7)の表現では，単位幅流量 q は水深の冪乗に比例する形式で記述することが多いため，冪数 m の値が重要な意味を有する．

∂h ∂q+= r (t )∂t ∂xv = αh m ,q = vh = αh m + 1（6）（7）ここに，v：断面平均流速［mm/h］，h：湛水深［mm］，q：単位幅流量［mm2/h］，α，mは流域特性を表すパラメータである．また，r(t)は有効降雨強度［mm/h］であり，流出計算の入力値となる．以下，入力降雨強度とする．直接流出は流出寄与域（河道及び河道近傍の湿潤領域）のみからの流出と考えると，流出寄与斜面長は実地形上の斜面長より十分短いものと考慮することができ，式(8)に示す変数分離形の近似式が仮定できる．

q( x, t ) ≅ xq* (t )

（8）ここに，q*：流出高［mm/h］である．また，流出寄与斜面長 L の末端であるx=Lとして式(8)を用いると，式(6)と式(7)を連立させて得られる偏微分方程式は流出高に関する常微分方程式に変形でき，さらに一般化を施すことによって，単一斜面からの降雨流出の基礎式（式(9)）が得られる．

― 78 ―

dq *dtβ= a0 q* ( r ( t ) − q* )（9）上式のαとmは土壌・地形特性を示す流出パラメータである．
1次元空間において，確率的な挙動を示すある時系列 x(t)の微小変化量 dx の時間発展は，以下に示す
伊藤過程によって x(t)の動きを表現することが可能である．

dx ( t ) = y ( x ( t ) , t ) dt + z ( x ( t ) , t ) dw ( t )（10）

ここにy(x(t), t)は時系列 x(t)と時間 t の任意の関数であり，dw(t)はWiener過程 w(t)の微小時間変化量である．なお，dw(t)は平均が0，標準偏差√t の正規分布に従っている．右辺第1項は決定論的な項（ドリフト項）であり，第2項は確率論的に決定される項である．上式は伊藤の確率微分方程式と呼ばれ，これと数学的に等価であるFokker-Planck方程式は式(11)によって記述される．

∂p( x(t ), t )∂y ( x(t ), t ) p( x(t ), t ) 1 ∂ 2 z 2 ( x(t ), t ) p ( x(t ), t )
=−+∂t
∂x2
∂x 2



（11）

降雨強度が各時刻において平均値とその周囲に分布したものであるとすれば，r ( t ) = r ( t ) + r’ ( t )（12）として表すことができる．ここで確率解析学の概念を適用し，入力降雨強度の分散がσ^2で表されるとすれば，式(11)はββdq * = a 0 q* ( r ( t ) − q* ) dt + a 0 q* σ 2 √TL dw ( t )（13）となる．この決定論的な項と確率項から構成される流出高に関する式をFokker-Planck方程式として記述
すると，β(β21 ∂ aq σ √TL∂ p ( q * ) ∂ aq * { r − q* }p( q* )+=∂t∂ q*2∂ q*2) p( q )
2*



（14）

として表され，この式によって降雨の不確実性に関する流出高の不確実性を求めることができる．また流出高に関するFokker-Planck方程式の定常解は― 79 ―

p (q * ) = p022

(a q β σ √TL )exp*( (2aσ 2 TL

1−βq*2 − βr q*−1−β 2−β))（15）となる．この式で得られる流出高の確率分布関数から水位の確率分布関数も求めることができる．

5. 3次元風速場
札幌レーダと新千歳空港レーダの観測値を組み合わせ，MUSCAT法により解析した結果について述べていく．結果としては高度ごとの水平断面，そして任意位置における鉛直断面を得ることができ，ともに反射強度の分布と風速ベクトルが描画されている（図5-図7）．これらの図は，今回の事例において降雨が強くなり始める時間帯に着目しており，左からそれぞれ日本時間で2014年9月11日の午前2時52分，午前3時2分，午前3時12分と10分間隔で推移していく．図5と図6はそれぞれ標高1.4km，3.5kmの位置における水平断面である．白丸が札幌レーダ，黒丸が新千歳空港レーダの位置である．新千歳レーダを原点とし，x 軸を東西距離，y 軸を南北距離とする．図7は図5，6中の線状降水帯に沿った赤線における鉛直断面である．縦軸は標高，横軸は水平距離を意味する．切断位置を図5，6に示すが切断位置は時間によらず同位置である．また，それぞれの図には反射強度の分布を示すカラーバーと風向風速の基準となるベクトルの大きさと基準値が示されている．図5，図6に示される水平断面について考察する．標高1.4kmを雲底高度として下層，3.5kmを中層と考える．両図の反射強度を見ると雨域が南西から北東へ伸びていることがわかり，線状降水帯の形状を有している．左から右へ時間を追っていくと反射強度の強い領域が北東へ移動していることが見て取れる．また，反射強度は降水強度として捉えるができる
ため，31dBZ以上の反射強度が強い部分（図では黄色）では積乱雲が発達していると解釈できる．
風の場を見てみると，下層ではどの時間でも全体として10～18m/s程の南東風が吹いているが，特に反射強度が強い部分では線状降水帯の進行方向に対して側方から強い風が収束している．2014年9月11日の道央には南東から暖かく湿った空気が流入しており，この南東風は降水帯に対し水蒸気を供給し続け積乱雲の発達，線状降水帯の維持に寄与していると考えられる．中層では約15m/sの南よりの風が吹いており，線状降水帯の進行方向に沿うような向きである．また，中層の風は図5のそれぞれの断面の中央付近における降水帯の終端において蛇行しており南西風に遷移している．
図7に示す鉛直断面について考察する．それぞれの断面は右端が降水帯の進行方向であり方角は北東となっている．反射強度の分布を見ると，雲の高さは約7~8kmである．34dBZ以上の反射強度の強い部分に着目すると積乱雲の位置がよくわかり，時間を追うと積乱雲は線状降水帯の進行方向へ移動していることが確認できる．鉛直方向の風の場を見てみると，上昇流と一部において下降流が確認できる．強い鉛直風の領域は線状降水帯の進行方向へ約17kmの周期で伝播しており反射強度が強い箇所ほど，また高度が高いほど大きい傾向にある．

― 80 ―

図5 : 標高1.4［km］
における反射強度［dBZ］と風速ベクトル［m/s］○:札幌レーダ, ●:新千歳空港レーダ
左から2014年9月11日午前2時52分, 3時2分, 3時12分（JST）

図6 : 標高3.5［km］
における反射強度［dBZ］と風速ベクトル［m/s］○:札幌レーダ, ●:新千歳空港レーダ
左から2014年9月11日午前2時52分, 3時2分, 3時12分（JST）

図7 : 鉛直断面における反射強度［dBZ］と風速ベクトル［m/s］
左から2014年9月11日午前2時52分, 3時2分, 3時12分（JST）― 81 ―

6. 確率微分方程式による流出量
図8は a)決定論的に流出計算（従来の流出計算）した結果，b)流出高についての確率微分方程式を用いた計算結果，c)流出高についてのFokker-Planck方程式から得られる計算結果である．いずれの計算においても，入力降雨強度は，ピーク降雨強度が50mm/hの1周期2時間のSin波形の降雨を与えており，
σ =5mm/hとして一定値，TL=1/6hを与えている．σ は降雨強度の変化に伴って，時間ないし降雨強度の関数として表せるはずであるが，本論では確率過程論に基づく，流出過程における不確実性評価の枠組
みの提示を主題としているため，この点は参考文献に譲り，最も単純な場合による結果を示している．
図8b）の確率微分方程式を解いた計算の試行回数は1000回であり，図8d）の確率微分方程式から得た結果のヒストグラムとFokker-Planck方程式から得た確率密度関数を，流出高ハイドログラフのピーク付近5.5hour時点において比較している．また，図8d）のヒストグラムと確率密度関数の比較からも確率微分方程式とFokker-Planck方程式が同一の現象を記述していることがわかる．常微分方程式型の貯留型の流出計算手法を確率微分方程式として解いた計算結果と，理論的に導き出された確率微分方程式に対応するFokker−Planck方程式を解いた解が一致しており，流出過程を確率過程と考えた場合に容易に任意時間の確率密度関数を求められる事がわかる．
本研究が使用する降雨流出過程には有効降雨強度が与えられる．有効降雨強度は降雨直前の流域の飽和度を損失雨量という値で代表することができ，損失雨量は保水能分布を通して有効降雨強度に影響を与えるものである．ここでは国土交通省水文水質データベースにおいてデータが公開されている全国106ヶ所の山地流域から6つの流域の計18の降雨イベントを対象に，降雨強度の不確実性と保水能の不確実性の双方が流出量の不確実性に与える影響を議論する．なお，本検討では流出量の不確実性に与える影響は2つの独立した要素(降雨強度と保水能)であるため，これらの確率分布は足し合わせが可能となる．そのため，Fokker-Plack方程式として記述される降雨流出過程の拡散項は2つ存在すると扱うことができる．保水能に起因する不確実性は過去に発生した強降雨イベントから定量化し，一方で降雨強度の不確実性は10, 15, 20%とした場合の流出量の不確実性を議論する．
図9は6つの山地流域の計18の降雨イベントにおける降雨強度と保水能の双方の不確実性が流出量の不確実性に与える影響を示す（左から降雨強度の不確実性を10, 15, 20%とした場合）．この図より流出量の不確実性に与える降雨強度と保水能の不確実性は概ね同程度であることがわかる．

― 82 ―

図8 : 決定論及び確率論的に求めたハイドログラフと流出量の確率密度関数の時間発展，及び確率微分方程式（SDE）から求めたヒストグラムとFokker-Planck方程式（FPE）で求めた確率密度関数の比較（σ=5mm/h，TL=1/6h）である．
a）はSin波形の降雨強度を入力としている．b）はその入力降雨強度に微小な乱れ r'（t）としてホワイトノイズが加わった降雨強度を入力としてSDEを解いている（SDEの計算試行回数N=1000）．c）は初期条件として極めて分散の小さい正規分布を与え，降雨の標準偏差σ，時定数TLは一定値としてFPEを解いている．d）はt=5.5hour
時点のq*のヒストグラムで，SDEとFPEの計算結果から得られる平均と標準偏差はそれぞれ，SDE：平均30.14mm/h，標準偏差1.53mm/h，FPE：平均30.12mm/h，標準偏差1.52mm/hであり，両者の差は極めて小さい．
同図は吉見ら（受理）より．

― 83 ―

図9 : 降雨強度と保水能の不確実性がピーク流出量の不確実性に与える影響
図の左側から降雨強度の不確実性が10, 15, 20%における結果を示す．

謝辞：本研究で使用したXバンドMPレーダデータは，国土交通省北海道開発局よりご提供頂いた．ここに記して深謝の意を表する．MUSCAT法による3次元風速場の推定に関しては，気象研究所第一研究室の山田芳則氏からご助言を頂戴した．
参考文献
1） 吉見和紘, Cho-Wen Wang, 山田正, 山田朋人: 確率過程論に基づいた降雨流出過程における不確実性評価の理論的枠組の提示. 土木学会論文集B1(水工学), 受理.
2） 北野慈和, 山田朋人: 北日本を通過する爆弾低気圧と太平洋ブロッキングとの関係. 土木学会論文集
B1(水工学), 受理.
3） Supraba, I. and T. J. Yamada, 2015: Uncertainty of peak runoff height associated with uncertainty of water
holding capacity and rainfall pattern. Annual Journal of Global Environment Engineering, JSCE, 71, 5,
I_17-I_22.
4） 鈴木章弘, 山田朋人, 2015: XバンドMPレーダと地上雨量計の降雨検知時間差を利用した降雨量の再推定. 土木学会論文集G(環境)，V0l.71，No.5，I_257-I_262.
5） 山田朋人, 和智光貴, 鈴木英一, 2015: 石狩川流域泥炭地における地盤沈下量に基づく地下水位変化の推定．土木学会論文集G(環境), Vol.71, No.5, I_55-I_60.
6） Inatsu, M., T. Sato, T. J. Yamada, R. Kuno, S. Sugimoto, M. A. Farukh, Y. N.Pokhrel, and S. Kure, 2015:
Multi-GCM by multi-RAM experiments for dynamical downscaling on summertime climate change in
Hokkaido. Atmos. Sci. Lett., 16,297-304.
7） 山田朋人, 2015: 将来の短時間豪雨の推定と降雨の不確実性―北海道を対象とする総合的ダウンスケール手法の開発と適用―. 日本環境衛生センター 月刊誌生活と環境, 39-44.
8） Supraba, I. and T. J. Yamada, 2015: Potential water storage capacity of mountainous catchments based oncatchment characteristics. Journal of Japan Society of Civil Engineers, Ser. B1 (Hydraulic Engineering),
Vol.71, No. 4, I_151-I_156.
9） 吉見和紘，山田正，山田朋人, 2015：確率微分方程式の導入による降雨流出過程における降雨の不確実性の評価，土木学会論文集B1（水工学），Vol.71，No.4，I_259-I_264.
10）Supraba, I. and T. Yamada, 2014: Catchment storage estimation based on total rainfall - total loss rainfall
relationship for 47 catchments in Japan. Journal of Japan Society of Civil Engineers, Ser.B1 (Hydraulic― 84 ―Engineering), 70, I_169-I_174.
11） Yamada, T.J., M.A. Farukh, T. Fukushima, M. Inatsu, T. Sato, Y. N. Pokhrel, and T. Oki, 2014: Extreme
precipitation intensity in future climates associated with the Clausius-Clapeyron-like relationship,
Hydrological Research Letters, Vol. 8, No. 4, p.108-113, DOI: 10.3178/hrl.8.108.

― 85 ―