【流出解析】シリーズ 第3回 貯留関数法【後編】
水文・水質・気象
VBA
水文学
流出解析
4.8．別の降雨イベントに対する流出量の推定
前編では，用意した降水量と流出量のデータに対して貯留関数法を適用しました． 入力データ(降水量)から出力データ(流出量)を推定できるように，今回使ったデータに合わせて\( T_{L} \)，\( K \)，\( p \)を調整しました． この作業はいわゆるキャリブレーション(較正)と呼ばれるものです．
※本来はバリデーション(検証)も必要ですが，本記事ではあえて省略しています．
さて，キャリブレーションをすれば流域の特性はわかりますが，これだけでは面白くないので流出量予測をしてみましょう． 実際の運用(防災や災害予測，ダム制御)を考えた時には，入力データ(降水量)のみから流出を予測できたら便利ですし．
はじめに，データを用意します． 今回用意したデータセットは2つの仮想的な降雨イベントを繋げたものです．1つ目の降雨イベントの終わり(\( t=15 \))と2つ目の降雨イベントの始まり(\( t=16 \))は，実際は連続しているわけではないので注意してください．
データを写真で示したワークシートのように入力してください． ワークシート名は「貯留関数法_予測」としてください． 降雨イベントフラグを用意することで，複数の降雨イベントに対する流出量をまとめてに計算できるようにしています．
| 時刻t | 降水量R (mm) | 降雨イベントフラグ 降雨開始：1，降雨終了：2 | 時刻t | 降水量R (mm) | 降雨イベントフラグ 降雨開始：1，降雨終了：2 |
|---|---|---|---|---|---|
| 1 | 1.5 | 1 | 21 | 1.5 | |
| 2 | 2 | 22 | 9.75 | ||
| 3 | 8 | 23 | 6 | ||
| 4 | 7 | 24 | 12.5 | ||
| 5 | 5.5 | 25 | 11 | ||
| 6 | 3 | 26 | 12 | ||
| 7 | 19.5 | 27 | 14 | ||
| 8 | 12 | 28 | 17.5 | ||
| 9 | 25 | 29 | 11.5 | ||
| 10 | 22 | 30 | 7.5 | 2 | |
| 11 | 24 | 31 | 0 | ||
| 12 | 28 | 32 | 0 | ||
| 13 | 35 | 33 | 0 | ||
| 14 | 23 | 34 | 0 | ||
| 15 | 15 | 2 | 35 | 0 | |
| 16 | 0.75 | 1 | 36 | 0 | |
| 17 | 1 | 37 | 0 | ||
| 18 | 4 | 38 | 0 | ||
| 19 | 3.5 | 39 | 0 | ||
| 20 | 2.75 | 40 | 0 | ||
| 1.5 | 41 | 0 |
はじめに，
\[
I=\dfrac{1}{3.6}fAr
\]
とおき，式(1)より
\[
\begin{align}
\dfrac{dS}{dt} &=\dfrac{1}{3.6}fAr-Q _{L}\\
\dfrac{S _{2}-S_{1}}{\Delta t} &= \left( \dfrac{I _{1}+I _{2}}{2}\right)-\left( \dfrac{Q _{L1}+Q _{L2}}{2}\right)
\end{align}
\]
ここで，
\[
R=\dfrac{r _{1}+ r_{2}}{2}
\]
とおくと，
\[
\dfrac{Q _{L2}}{2}+ \dfrac{S _{2}}{ \Delta t} =\dfrac{fA}{3.6}R - \dfrac{Q _{L1}}{2}+\dfrac{S _{1}}{ \Delta t} \\
\]
なお，\( \Delta t \)を1(hr)とします． 添え字1のときを時刻\( t \)，添え字2のときを時刻\( t+1 \)として，この式に式(3)を代入すると，
\[
\begin{align}
\dfrac{Q _{L2}}{2}+ \dfrac{S _{2}}{ \Delta t} &=\dfrac{fA}{3.6}R - \dfrac{Q _{L1}}{2}+\dfrac{S _{1}}{ \Delta t} \\
Q _{L}(t+1)+2K{Q _{L}(t+1)} ^{p}&=\dfrac{2fA}{3.6}R - Q _{L}(t)+2K{Q _{L}(t)} ^{p}\\
\end{align}
\]
\( Q_{L}(t+1) \)に注目するとき，\( Q_{L}(t) \)は1ステップ前の値で既知なので，
\[
\dfrac{2fA}{3.6}R - Q _{L}(t)+2K{Q _{L}(t)} ^{p}=C=const.\\
\]
とおくことができる．したがって，\( X=Q_{L}(t+1) \)とおくと，
\[
G(X)=X+2KX ^{p}-C=0
\]
ここで，\( G \)関数記号である．\( G(X)=0 \)を満たす\( X \)が時刻\( t+1 \)における流出量\( Q_{L}(t+1) \)です． 方程式の解を求めるにはニュートン法を使います．
はじめに，Xにテキトーな値を入れます．これを\( X_{0} \)とします． 今回は\( X_{0}=1(m^{3}/s) \)とします． これが初期条件です．
関数\( G(X) \)上の点(\( X_{0} \)，\( G(X_{0}) \))における接線が(\( X \)軸と接するときの値が第一近似値であり，これを(\( X_{1} \)とします．
\[
X _{1}=X _{0}- \dfrac{G(X _{0})}{G'(X _{0})}
\]
同様に，第二近似値(\( X_{2} \))は
\[
X _{2}=X _{1}- \dfrac{G(X _{1})}{G'(X _{1})}
\]
第\( i+1 \)近似値は
\[
X _{i+1}=X _{i}- \dfrac{G(X _{i})}{G'(X _{i})}
\]
このように計算を繰り返し，(\( X_{i+1} \))と(\( X_{i} \))の差の絶対値が0.01を下回った時，解にたどり着いたと判断し，計算を終了します．
式(3)を使って\( Q_{L} \)を観測される直接流出量\( Q \)に変換し，基底流出量(4.1のデータの\( Q(1)=122 (m^{3}/s) \)とする)を加えて，全流出量を求めます．
※基底流出量の妥当性については記事の最後に考察します．
5．アルゴリズム
第4章の計算プロセスをアルゴリズム(計算手順書)で書き表すと以下のようになります．
6．プログラム
貯留関数法の計算用プログラムです． エクセルのVBA機能を使うと貯留関数法の計算ができます． 2つ目のStorageFunctionMethodプロシージャを実行するだけで一連の計算が自動実行されます．
なお，今回のプログラムを全て理解する必要はないですが，よかったら挑戦してみてください． 一見すると長いプログラムですが，前編の内容を理解できていればプログラムの内容も理解できると思います． ちなみに，エクセルに標準搭載されている関数のみを使っても同じ計算をすることは可能です．
VBAの使い方がわからないという方はこちらのサイトをどうぞ↓
Option Explicit
Function MergeCells(ByVal C1 As Range, ByVal C2 As Range)
'指定した範囲のセルを結合する
Range(Replace(C1.Address, "$", "") & ":" & Replace(C2.Address, "$", "")).Select
With Selection
.HorizontalAlignment = xlCenter
.VerticalAlignment = xlCenter
.WrapText = True
.Orientation = 0
.AddIndent = False
.IndentLevel = 0
.ShrinkToFit = False
.ReadingOrder = xlContext
.MergeCells = False
End With
Selection.merge
End Function
Sub StorageFunctionMethod()
'最終更新日時：20201026
'川村三郎：土木工学プログラム集 水文・水理1，森北出版，1984に準拠
'2.4 貯留関数法
'※注：QLは遅れ時間TLがないと仮定したときの流出量，Qは実際に観測される流出量
'入力データ:
'流域面積(km2)：B2セル
'流路延長(km)：B3セル
'時刻T(h):A列(6行目から入力)
'降雨量R(mm):B列(6行目から入力)
'流出量Q(m3/s):C列(6行目から入力)
'画面の更新を停止
Application.ScreenUpdating = False
'計算開始メッセージの表示
Debug.Print "計算開始 ", Now()
'計算用シートのアクティブ化
Worksheets("貯留関数法_較正").Activate
'変数の定義
Dim i As Integer, Cnt As Integer 'カウント用
Dim A As Double '流域面積(km2)
Dim L As Double '流路延長(km)
Dim Qb As Double '基底流出量(m3/s) (base flow)
Dim Qd As Double '直接流出量(m3/s) (direct flow)
Dim N As Double '流出量のデータ数
Dim t1 As Double 't1，t2はハイドログラフの増水側と減水側で流出量Qがおよそ同じになる時刻(プログラムでは行番号を格納)
Dim t2 As Double 't1，t2はハイドログラフの増水側と減水側で流出量Qがおよそ同じになる時刻(プログラムでは行番号を格納)
Dim diff As Double '時刻t1と時刻t2の流出量Qの差の絶対値(mm)
Dim RT As Double '時刻t1〜t2までの総降雨量(mm)
Dim Q_Max As Double 'ピーク流出量(m3/s)
Dim Q_Max_Row As Double 'ピーク流出量のデータの行番号
Dim TPL As Double '遅れ時間TLがないと仮定したとき，ピーク流出量が観測される時刻(プログラムでは行番号を格納)
Dim TL_Kimura As Double '木村の式による遅れ時間TLの第一近似値(h)
Dim TL As Double '遅れ時間(h)
Dim OptTL As Integer '最適な遅れ時間(h)
Dim f As Double 'ピーク流出係数(-)
Dim delta As Integer '積分における区間の分割幅
Dim QL() As Double '遅れ時間TLがないと仮定したときの流出量
Dim QT As Double '時刻t1〜t2における流出量QLの総量
Dim FI() As Double '時刻t1〜任意時刻tにおける総面積有効雨量
Dim QS() As Double '時刻t1〜任意時刻tにおける流出量QLの総量
Dim S() As Double '任意時刻tにおける貯留量S
Dim SAR As Double 'Q-S曲線において，ループ曲線の上側とx軸(Q)で囲まれる面積
Dim SAD As Double 'Q-S曲線において，ループ曲線の下側とx軸(Q)で囲まれる面積
Dim SA As Double 'Q-S曲線において，ループ曲線の内側の面積
Dim SA_Min As Double 'ループ曲線内部の面積の最小値
'最小二乗法に関する変数
Dim e As Integer '最小二乗法に用いるデータセットの数
Dim p As Double '貯留関数の定数p
Dim K As Double '貯留関数の定数K
Dim Σx As Double 'QLの自然対数の合計
Dim Σy As Double 'Sの自然対数の合計
Dim Σxy As Double 'QLの自然対数とSの自然対数の積和
Dim Σx2 As Double 'QLの自然対数の2乗の合計
'***************Step1：最適な遅れ時間OptTLを求める***************
'変数の初期化
A = Cells(2, 2) '流域面積(km2)
L = Cells(3, 2) '流路延長(km)
N = 0 '流出量のデータ数
Qb = Cells(6, 3) '基底流出量(m3/s)
RT = 0 '時刻t1〜t2における総降雨量(mm)
Q_Max = 0 'ピーク流出量(m3/s)
Q_Max_Row = 0 'ピーク流出量のデータの行番号
delta = 1 '1時間単位で積分計算
e = 0 '最小二乗法に用いるデータセットの数
'ラベル作成
Cells(2, 4) = "直接流出量Qd (m3/s)"
Call MergeCells(Cells(2, 4), Cells(5, 4))
'流出量データの数N，ピーク流出量Q_Max，直接流出量Qdを求める
i = 6
Do
'直接流出量Qd
Cells(i, 4) = Cells(i, 3) - Qb
'ピーク流出量Q_Max
If Q_Max < Cells(i, 3) Then
Q_Max = Cells(i, 3)
Q_Max_Row = i
End If
'流出量データの数N
N = N + 1
i = i + 1
Loop Until Cells(i, 1) = ""
'出力
Debug.Print "データ数は" & N
Debug.Print "ピーク流出量は" & Q_Max & " (m3/s)"
Debug.Print "ピーク流出量のデータは" & Q_Max_Row & "行目"
'配列の再定義・初期化
ReDim QL(1 To N), FI(6 To 6 + (N - 1)), QS(6 To 6 + (N - 1)), S(6 To 6 + (N - 1))
For i = 1 To N
QL(i) = 0 '遅れ時間TLがないと仮定したときの流出量
FI(i + 5) = 0 '時刻t1〜任意時刻tにおける総面積有効雨量
QS(i + 5) = 0 '時刻t1〜任意時刻tにおける流出量QLの総量
S(i + 5) = 0 '任意時刻tにおける貯留量S
Next i
'ハイドログラフの増水側と減水側で同じ流出量になる時刻t1とt2を求める
'このときの流出量はピーク流出量を1〜2割とする
'今回はピーク流出量の1.5割程度のときに増水側と減水側で流出量の差が最小となる時刻とした
i = 6
diff = Cells(Q_Max_Row, 3) '初期化(とりあえず大きな値を放り込みたい)
Do
Cnt = Q_Max_Row + 1
Do
If Abs(Cells(i, 3) - Cells(Cnt, 3)) < diff Then
diff = Abs(Cells(i, 3) - Cells(Cnt, 3))
t1 = i
t2 = Cnt
End If
Cnt = Cnt + 1
Loop Until Cells(Cnt, 1) = ""
i = i + 1
Loop Until Cells(i, 3) >= Q_Max * 0.15
'出力
Debug.Print "時刻t1は" & t1 & "行目(時刻：" & Cells(t1, 1) & ")"
Debug.Print "時刻t2は" & t2 & "行目(時刻：" & Cells(t2, 1) & ")"
Debug.Print "Q(t1)とQ(t2)の差の絶対値は" & diff & " (m3/s)"
'時刻t1〜t2の間の総降雨量を求める
'時刻t1の行の雨量データは直前の1時間の量なので，t1+1行目から合計を求める
For i = t1 + 1 To t2
RT = RT + Cells(i, 2)
Next i
'出力
Debug.Print "時刻t1〜t2の間の総降雨量は" & RT & " (mm)"
'木村の式で遅れ時間TLの第一次近似値を求める
If L < 11.9 Then
TL_Kimura = 0
ElseIf L >= 11.9 Then
TL_Kimura = Int(0.5 + 0.047 * L - 0.56) '四捨五入して整数にする
End If
'出力
Debug.Print "木村の式による遅れ時間TLは" & TL_Kimura & " (h)"
'結果出力用ラベルの作成
Cells(N + 7, 4) = "遅れ時間TL (h)"
Cells(N + 8, 4) = "ピーク流出量Q_Max (m3/s)"
Cells(N + 9, 4) = "遅れ時間TLがないと仮定したとき，ピーク流出量が観測される時刻TPL"
Cells(N + 10, 4) = "時刻t1〜t2に対するQL(t)の積分値QT"
Cells(N + 11, 4) = "ピーク流出係数f"
Cells(N + 12, 4) = "ループ曲線の上側とS=Qで囲まれる面積"
Cells(N + 13, 4) = "ループ曲線の下側とS=Qで囲まれる面積"
Cells(N + 14, 4) = "ループ曲線で囲まれる面積"
Cells(N + 15, 4) = "ループ判定"
Cells(1, 4) = "出力データ"
Cells(2, 5 + TL_Kimura) = "遅れがない場合の直接流出量QL (m3/s)"
Cells(2, 9 + TL_Kimura) = "時刻t1〜任意時刻tにおける総面積有効雨量FI"
Cells(2, 13 + TL_Kimura) = "時刻t1〜任意時刻tにおける流出量QLの総量QS"
Cells(2, 17 + TL_Kimura) = "任意時刻tにおける貯留量S"
Cells(2, 17 + TL_Kimura + 2) = "最適な遅れ時間に基づく貯留関数"
Call MergeCells(Cells(1, 4), Cells(1, 4 + (TL_Kimura + 2) * 4 + 7))
Call MergeCells(Cells(2, 5), Cells(4, 5 + (TL_Kimura + 1)))
Call MergeCells(Cells(2, 9), Cells(4, 9 + (TL_Kimura + 1)))
Call MergeCells(Cells(2, 13), Cells(4, 13 + (TL_Kimura + 1)))
Call MergeCells(Cells(2, 17), Cells(4, 17 + (TL_Kimura + 1)))
Call MergeCells(Cells(2, 17 + (TL_Kimura + 1) + 1), Cells(4, 17 + (TL_Kimura + 1) + 7))
'最適な遅れ時間OptTLを求める
'遅れ時間0〜遅れ時間TL_Kimura+1まで1時間刻みで貯留関数を求める
For TL = 0 To TL_Kimura + 1
'変数の初期化
Q_Max = 0 'ピーク流出量(m3/s)
QT = 0 '時刻t1〜t2に対するQL(t)の積分値QT
SAR = 0 'Q-S曲線において，ループ曲線の上側とx軸(Q)で囲まれる面積
SAD = 0 'Q-S曲線において，ループ曲線の下側とx軸(Q)で囲まれる面積
'実際の流出量は遅れ時間TLの結果生じたものと考え，遅れ時間TLがないと仮定したときの流出量QLを求める(直接流出量をTLだけ戻したデータセットを作る※)
'※TL>0のときは，TL-1のデータも一部使用するため，厳密には直接流出量をTLだけ戻したデータセットというわけではない
'なぜ，TL-1のデータも引き継ぐのかは不明(積分で不都合が生じないようにするため？)
For i = 6 To 6 + (N - 1) - TL
'QLを求める
If i <= 6 + (N - 1) - TL * 2 Then
QL(i - 5) = Cells(i + TL, 4)
End If
Cells(i, 5 + TL) = QL(i - 5)
'直接流出量QLの最大値とその時の行を取得
If Cells(i, 5 + TL) >= Q_Max Then
Q_Max = Cells(i, 5 + TL)
TPL = i
End If
Next i
'TL，Q_Max，TPLを出力
Cells(5, 5 + TL) = "TL=" & TL
Cells(5, 9 + TL) = "TL=" & TL
Cells(5, 13 + TL) = "TL=" & TL
Cells(5, 17 + TL) = "TL=" & TL
Cells(N + 7, 5 + TL) = TL
Cells(N + 8, 5 + TL) = Q_Max
Cells(N + 9, 5 + TL) = TPL - 5
'QTを求める
'for分がt2までではなく，t2-TLで終わっているのは単に計算の省略のため
For i = t1 + 1 To t2 - TL
QT = QT + Cells(i, 5 + TL) * delta
Next i
'ピーク流出係数fを求める
f = QT / ((A / 3.6) * RT)
'QTとfの出力
Cells(N + 10, 5 + TL) = QT
Cells(N + 11, 5 + TL) = f
'時刻t1〜任意時刻tにおける総面積有効雨量
For i = t1 + 1 To t2
FI(i) = FI(i - 1) + A * f * Cells(i, 2) / 3.6
Cells(i, 9 + TL) = FI(i)
Next i
'時刻t1から任意時刻tまでの総流出量QS(i)と時刻tにおける貯留量S(t)を求める
'なぜt1+1から始めてとt2-TL-1までしかS(t)を計算しないかは謎
'↑ループの始点(t1)と終点(t2)が交わるのを避けるための安全策？
For i = t1 + 1 To t2 - 1
'QS(t)を求める
If i <= t2 - TL - 1 Then
QS(i) = QS(i - 1) + (Cells(i - 1, 5 + TL) + Cells(i, 5 + TL)) / 2
End If
Cells(i, 13 + TL) = QS(i)
'任意時刻tにおける貯留量S(t)を求める
If i <= t2 - TL - 1 Then
If FI(i) > QS(i) Then
S(i) = FI(i) - QS(i)
End If
End If
Cells(i, 17 + TL) = S(i)
Next i
'SARを求める
For i = t1 + 1 To TPL - 1
SAR = SAR + (S(i) + S(i + 1)) / 2 * Abs(QL(i - 5 + 1) - QL(i - 5))
Next i
'SADを求める
For i = TPL To t2 - TL - 1
SAD = SAD - (S(i) + S(i + 1)) / 2 * Abs(QL(i - 5) - QL(i - 5 + 1))
Next i
SA = SAR + SAD
'上側面積，下側面積，ループ面積の出力
Cells(N + 12, 5 + TL) = SAR
Cells(N + 13, 5 + TL) = SAD
Cells(N + 14, 5 + TL) = SA
'ループ内の面積が最小となるときの遅れ時間を最適な遅れ時間OptTLとする
If TL = 0 Then
SA_Min = Abs(SA)
OptTL = TL
Else
If SA_Min > Abs(SA) Then
SA_Min = Abs(SA)
OptTL = TL
End If
End If
'ループが時計周りか反時計周りかを判定
If SA > 0 Then
Cells(N + 15, 5 + TL) = "clockwise"
ElseIf SA < 0 Then
Cells(N + 15, 5 + TL) = "counter clockwise"
Else
Cells(N + 15, 5 + TL) = "just fitting"
End If
Next TL
'出力
Debug.Print "最適な遅れ時間TLはTL="; OptTL
'***************Step2：最適な遅れ時間OptTLに対する貯留関数を求める***************
'結果出力用ラベルの作成
Cells(5, 4 * (TL_Kimura + 2) + 5) = "時刻t"
Cells(5, 4 * (TL_Kimura + 2) + 6) = "流出量QL(m3/s)"
Cells(5, 4 * (TL_Kimura + 2) + 7) = "面積有効雨量"
Cells(5, 4 * (TL_Kimura + 2) + 8) = "時刻t1〜任意時刻tにおける総面積有効雨量"
Cells(5, 4 * (TL_Kimura + 2) + 9) = "時刻t1〜任意時刻tにおける流出量QLの総量"
Cells(5, 4 * (TL_Kimura + 2) + 10) = "任意時刻tにおける貯留量S"
Cells(5, 4 * (TL_Kimura + 2) + 11) = "S=K・Q ^ p"
Cells(N + 7, 4 * (TL_Kimura + 2) + 5) = "遅れ時間TL (h)"
Cells(N + 8, 4 * (TL_Kimura + 2) + 5) = "ピーク流出量Q_Max (m3/s)"
Cells(N + 9, 4 * (TL_Kimura + 2) + 5) = "遅れ時間TLがないと仮定したとき，ピーク流出量が観測される時刻TPL"
Cells(N + 10, 4 * (TL_Kimura + 2) + 5) = "時刻t1〜t2における流出量QLの総量"
Cells(N + 11, 4 * (TL_Kimura + 2) + 5) = "ピーク流出係数f"
Cells(N + 12, 4 * (TL_Kimura + 2) + 5) = "貯留関数の定数p"
Cells(N + 13, 4 * (TL_Kimura + 2) + 5) = "貯留関数の定数K"
'変数の初期化
QT = 0
TL = OptTL
For i = 1 To N
FI(i + 5) = 0
QS(i + 5) = 0
S(i + 5) = 0
Next i
'時刻データセットを作成
For i = 6 To 6 + (N - 1)
Cells(i, 4 * (TL_Kimura + 2) + 5) = i - 5
Next i
'遅れ時間TLがないときの流出量QLを求める
For i = 1 To N - TL
If i <= N - TL * 2 Then
QL(i) = Cells((i + 5) + TL, 4)
End If
Cells(i + 5, 4 * (TL_Kimura + 2) + 6) = QL(i)
'直接流出量QLの最大値とその時の行を取得
If QL(i) >= Q_Max Then
Q_Max = QL(i)
TPL = i + 5
End If
Next i
'TL，Q_Max，TPL-5を出力
Cells(N + 7, 4 * (TL_Kimura + 2) + 6) = TL
Cells(N + 8, 4 * (TL_Kimura + 2) + 6) = Q_Max
Cells(N + 9, 4 * (TL_Kimura + 2) + 6) = TPL - 5
'時刻t1〜t2における流出量QLの総量
For i = t1 + 1 To t2 - TL
QT = QT + Cells(i, 5 + TL) * delta
Next i
'ピーク流出係数の計算
f = QT / ((A / 3.6) * RT)
'QTとfの出力
Cells(N + 10, 4 * (TL_Kimura + 2) + 6) = QT
Cells(N + 11, 4 * (TL_Kimura + 2) + 6) = f
'FIを求める
For i = 6 To 6 + (N - 1)
'面積有効雨量
Cells(i, 4 * (TL_Kimura + 2) + 7) = A * f * Cells(i, 2) / 3.6
'時刻t1〜任意時刻tまでの面積有効雨量
If i = 6 Then
FI(i) = Cells(i, 4 * (TL_Kimura + 2) + 7)
Else
FI(i) = FI(i - 1) + Cells(i, 4 * (TL_Kimura + 2) + 7)
End If
Cells(i, 4 * (TL_Kimura + 2) + 8) = FI(i)
Next i
'QS，Sを求める
For i = 7 To 6 + (N - 1) - TL
QS(i) = QS(i - 1) + (QL(i - 5) + QL(i - 5 - 1)) / 2
If FI(i) > QS(i) Then
S(i) = FI(i) - QS(i)
End If
If i <= 6 + (N - 1) - TL - 1 Then
Cells(i, 4 * (TL_Kimura + 2) + 9) = QS(i)
Cells(i, 4 * (TL_Kimura + 2) + 10) = S(i)
End If
Next i
'***************Step3：最小二乗法により貯留関数の定数p，Kを求める***************
'変数の初期化
Σx = 0
Σy = 0
Σxy = 0
Σx2 = 0
e = (N - 1) - TL - 1
'最上二乗法
For i = 7 To 6 + (N - 1) - TL - 1
Σx = Σx + Log(QL(i - 5))
Σy = Σy + Log(S(i))
Σxy = Σxy + Log(QL(i - 5)) * Log(S(i))
Σx2 = Σx2 + Log(QL(i - 5)) ^ 2
Next i
p = (Σx * Σy - e * Σxy) / (Σx * Σx - e * Σx2)
K = Exp((Σy - p * Σx) / e)
'pとKの出力
Cells(N + 12, 4 * (TL_Kimura + 2) + 6) = p
Cells(N + 13, 4 * (TL_Kimura + 2) + 6) = K
Debug.Print "p="; p
Debug.Print "K="; K
'S=KQL ^ pを求める
For i = 6 To 6 + (N - 1) - TL
Cells(i, 4 * (TL_Kimura + 2) + 11) = K * QL(i - 5) ^ p
Next i
'***************Step4：別の降雨イベントの流出量を予測する***************
Call RunoffPrediction(K, p, N, f, A, OptTL, Qb)
'画面を更新
Application.ScreenUpdating = True
'計算終了メッセージの表示
Debug.Print "計算終了 ", Now()
End Sub
Sub RunoffPrediction(ByVal K As Double, ByVal p As Double, ByVal N As Integer, ByVal f As Double, A As Double, ByVal OptTL As Integer, ByVal Qb As Double)
'20201030作成
'貯留関数法を使って降水量データから流出量を予測する
'引数
'K：貯留関数の定数K
'p：貯留関数の定数p
'N：キャリブレーションデータの数(流出過程の開始から終了までの時間をカバーできるある程度大きな値ならなんでもいい)
'f：ピーク流出係数
'A：流域面積(km2)
'OptTL：最適な遅れ時間
'QI：基底流出
'計算開始メッセージの表示
Debug.Print "流出量予測開始 ", Now()
'計算用シートのアクティブ化
Worksheets("貯留関数法_予測").Activate
'変数の定義
Dim i As Integer '汎用カウンター
Dim j As Integer '降水量データを格納する配列R()の要素番号
Dim Cnt As Integer '結果出力用カウンター
Dim e As Integer '入力データ(降水量データ)用カウンター
Dim R() As Double '降水量データ
Dim QL() As Double '遅れ時間TLがないと仮定したときの流出量
Dim Q() As Double '流出量
Dim Tf As Double '流出過程の終了までにかかると想定される時間の最大値
Dim Y As Double '目的関数G(X)
Dim X As Double '目的関数G(X)の説明変数X
Dim C As Double '目的関数G(X)の定数項C
Dim D As Double '関数G(X)の導関数G'(X)
Dim X2 As Double '最適解に対する(ニュートン法で得られた)第n次近似値
'結果出力用ラベルの作成
Cells(1, 5) = "出力データ"
Cells(2, 5) = "No."
Cells(2, 6) = "降水量R (mm)"
Cells(2, 7) = "流出量Q (m3/s)"
Cells(2, 8) = "流出過程フラグ"
If TypeName(Cells(2, 8).Comment) <> "Comment" Then
Cells(2, 8).AddComment
End If
Cells(2, 8).Comment.Visible = True
Cells(2, 8).Comment.Text Text:="1：流出過程開始，2：流出過程終了" & Chr(10) & ""
Call MergeCells(Cells(1, 5), Cells(1, 8))
Call MergeCells(Cells(2, 6), Cells(2, 6))
Call MergeCells(Cells(2, 7), Cells(2, 7))
Call MergeCells(Cells(2, 8), Cells(2, 8))
'変数の初期化
Tf = 1.5 * N
Cnt = 3
'配列の再定義
ReDim QL(1 To Tf), R(1 To Tf), Q(1 To Tf)
'降水量データが尽きるまで計算を繰り返す
e = 3
Do
'変数の初期化
j = 1 '降水量データを格納する配列R()の要素番号
'配列の初期化
For i = 1 To Tf
QL(i) = 0
R(i) = 0
Q(i) = 0
Next i
'降水量データの取得
'降雨イベント開始フラグを取得
Do Until Cells(e, 3) = 1
e = e + 1
'入力データが尽きたら計算終了
If Cells(e, 1) = "" Then
Debug.Print "流出量予測終了"; Now()
Exit Sub
End If
Loop
Debug.Print "降雨イベント開始行：", e
'降雨イベント終了フラグが現れるまで降水量データを取得し続ける
Do
R(j) = Cells(e, 2)
e = e + 1
j = j + 1
Loop Until Cells(e - 1, 3) = 2 Or Cells(e, 1) = ""
Debug.Print "降雨イベント終了行：", e - 1
'ニュートン法により最適解を求める
For i = 1 To Tf
X = 1 '初期条件X=QL(i+1)=1
'目的関数G(X)の定数項C
C = 2 * f * A * R(i) / 3.6 - QL(i) + 2 * K * QL(i) ^ p
'反復法
Do
Y = X + 2 * K * X ^ p - C '目的変数G(X)
D = 1 + 2 * p * K * X ^ (p - 1) '目的変数G(X)
X2 = X - Y / D '最適解に対する(ニュートン法で得られた)第n次近似値
'収束判定
'第n次近似値と第n-1次近似値の差が0.01未満になったら最適解にたどり着いたと判定
If Abs(X2 - X) < 0.01 Then
QL(i + 1) = X
Exit Do
Else
X = X2
End If
Loop
'雨の降り始めから遅れ時間OptTLだけ経過するまでは，流出量への降雨の影響はないと考える
'すなわち，流出量は基底流出に等しい
If i <= OptTL Then
Q(i) = Qb
End If
'流出量を求める
Q(i + OptTL) = QL(i + 1) + Qb
'結果の出力
Cells(Cnt, 5) = Cnt - 2
Cells(Cnt, 6) = R(i)
Cells(Cnt, 7) = Q(i)
If i = 1 Then
'降雨イベントの始まりを表すフラグを立てる
Cells(Cnt, 8) = 1
End If
Cnt = Cnt + 1
'流出量がおおよそ減衰したら計算終了
If i > 6 Then
If Q(i) < 1.1 * Q(5) Then
'降雨イベントの終わりを表すフラグを立てる
Cells(Cnt - 1, 8) = 2
Exit For
End If
End If
Next i
Loop Until Cells(e, 1) = ""
End Sub
7．結果
細かくて申し訳ありませんが，計算結果はこのようになります．
遅れ時間を考慮したハイドログラフ(流出量波形)
ループ曲線
最適な遅れ時間\( OptT_{L} \)において，雨の降り始めから貯留量と直接流出量を計算するとこんな感じになります↓
ちなみに予測結果はこんな感じです．
8．おわりに
最後までお付き合いいただき，ありがとうございました．
さて，実は今回の記事にはいくつか問題があります． 1．\( t_{1} \)，\( t_{2} \)の選び方が恣意的． 2．流出量予測をするときの基底流出量の値として，キャリブレーションの時の値をそのまま使っている． 3．貯留量などの計算において，\( T_{L}=0 \)から\( T_{L-kimura}+1 \)まで順に計算していくが，当該の\( T_{L} \)における計算結果の一部に，一つ前の\( T_{L} \)における計算結果が影響を与える． 4．プログラムは河村先生の書籍(記事の最後に記載)を参考にしたが，プログラムの一部は書籍通りではない． 5．計算過程もエクセル上に表示していること． 6．木村俊晃さんの論文を読んでいない．
1と2については筆者自身も問題に感じていますが，記事の本来の目的とは異なってくるので，今回は触れないこととしました． 将来的には，\( t_{1} \)，\( t_{2} \)の選び方に一定の明確な基準を設ける．流出量予測をするにしても雨の降り始めの流出量を基底流出量とするといった対策が必要だと思っています．
3，4については，なぜ河村先生が書籍に載っていたようなプログラム構成にしたのか私自身まだ十分に理解できていない，あるいは納得のいっていない部分でもあるため，今の段階ではご容赦いただけると幸いです． いずれは解決したいと思っています．
5については，本来はプログラムの計算過程をエクセルに表示することはありませんが，今回はわかりやすさを優先して表示することとしました．
6について．こういった記事を書く場合，原著論文を読むべきですが，古い論文で入手できなかったためお許しください． 今回紹介したこの貯留関数法，実は木村俊晃さんという方が創始したものだそうです． 同じ日本人として誇らしく思いますね．
実は遅れ時間の計算方法は今回示した方法以外でも可能です． その他，貯留関数法についてもっと詳しく知りたい人は下記の角屋先生の講座をご覧ください． 貯留関数法を手計算で解いているみたいなので，参考になるかと思います．
本記事が農業農村工学の発展に少しでも寄与できれば幸いです． 執筆にあたり，書籍の内容について紹介することを快く承諾いただいた森北出版の皆様と，何より書籍の著者である河村三郎先生には心より感謝いたします．
記事の内容に関してお気づきの点やご質問等がありましたら，ご連絡いただけますと幸いです． リクエスト等も受け付けておりますので，ご遠慮なく連絡ください． Twitterアカウント：エビぐんかん@6LxAi9GCOmRigUI メール：nnCreatorCircle@gmail.com
【流出解析】シリーズ 一覧
【流出解析】シリーズ 第1回 流域界の決定と入力データの用意 【流出解析】シリーズ 第2回 貯留関数法【前編】 【流出解析】シリーズ 第3回 貯留関数法【後編】 【流出解析】シリーズ 第4回 流域実蒸発散量の推定―蒸発散比法― 【流出解析】シリーズ 第5回 実蒸発散量の推定―補完法― 【流出解析】シリーズ 第6回 タンクモデル 【流出解析】シリーズ 第7回 積雪・融雪モデル ―降雪量・融雪量の推定― 【流出解析】シリーズ 第8回 タンクモデルの最適化 【流出解析】シリーズ 第9回 単位図法 【流出解析】シリーズ 第10回 流出量・流出高・比流量の違い 【流出解析】シリーズ 第11回 表面流モデル【定常降雨：前編】 【流出解析】シリーズ 第12回 表面流モデル【定常降雨：後編】
検索用キーワード：流出解析，流出モデル，水文学，洪水，治水，利水，数理モデル，貯留関数，表面流，タンクモデル，角屋睦，永井明博，
引用・参考文献
河村三郎：土木工学プログラム集 水文・水理1 実用プログラムとグラフィック表示，森北出版株式会社，1984，pp.55-72 丸山利輔，石川重雄，大槻恭一，高瀬恵次，永井明博，田中丸治哉，駒村正治，赤江剛夫，堀野治彦，三野徹，武田育郎，金木亮一，渡辺紹裕：地域環境水文学，朝倉書店，2014，pp.64-66 角屋睦，永井明博，1980．講座 流出解析手法(その10) 4．貯留法―貯留関数法による洪水流出解析．農業土木学会誌48(10)，747-754．
