【流出解析】シリーズ 第9回 単位図法<h1>【流出解析】シリーズ 第9回 単位図法</h1>
NN広場
農業農村工学のための総合サイト
NNテックブログ
農業農村工学向け解説記事
NN広場 農業農村工学のための総合サイト NNテックブログ 農業農村工学向け解説記事
【流出解析】シリーズ 第9回 単位図法
水文・水質・気象
Excel
VBA
モデル
水文学
流出解析

1.はじめに

今回は『単位図法』について紹介したいと思います. 表面流モデルに関する記事を現在執筆中ですが,思ったよりも時間がかかりそうなので繋ぎの記事として単位図法に登場してもらいました笑

単位図法というと「変な数列の式が出てきてよくわからん」という方もいるかもしれませんが,実際にはものすごくシンプルなモデルで計算も簡単です. ぜひ,皆さんも計算にチャレンジしてみてください!

今回使用した計算用のワークシートは下記URLからダウンロードできます. ご自由にお使いください. 単位図.xlsm

【流出解析】シリーズ 一覧

【流出解析】シリーズ 第1回 流域界の決定と入力データの用意 【流出解析】シリーズ 第2回 貯留関数法【前編】 【流出解析】シリーズ 第3回 貯留関数法【後編】 【流出解析】シリーズ 第4回 流域実蒸発散量の推定―蒸発散比法― 【流出解析】シリーズ 第5回 実蒸発散量の推定―補完法― 【流出解析】シリーズ 第6回 タンクモデル 【流出解析】シリーズ 第7回 積雪・融雪モデル ―降雪量・融雪量の推定― 【流出解析】シリーズ 第8回 タンクモデルの最適化 【流出解析】シリーズ 第9回 単位図法 【流出解析】シリーズ 第10回 流出量・流出高・比流量の違い 【流出解析】シリーズ 第11回 表面流モデル【定常降雨:前編】 【流出解析】シリーズ 第12回 表面流モデル【定常降雨:後編】

2.単位図法とは

単位図法の概要については下記の書籍の言葉を借りて説明したいと思います.

丸山利輔,石川重雄,大槻恭一,高瀬恵次,永井明博,田中丸治哉,駒村正治,赤江剛夫,堀野治彦,三野徹,武田育郎,金木亮一,渡辺紹裕:地域環境水文学,朝倉書店,2014,pp.63-64,77.

単位図法は,単位図というものを使って流出量を計算する流出解析手法で,Shermanという人によって1932年に提案されました. 単位図法は海外では今でも使われてるそうですが,日本ではあまり用いられていないそうです. 単位図法は線形性を仮定した手法ですが,日本の河川の流出特性は非線形性が強いため単位図では流出現象を十分に説明できないそうです. 実際,うちの某先生も「日本にはあんまり合わない」って言ってました.

単位図は以下の1つの定義と3つの仮定から構成されます. 1.単位図の定義 単位時間\( t_{0} \)に有効降雨強度\( I_{e} \)があると,それにより直接流出が生じる.いま\( R=I_{e}t_{0}=1 \)の単位有効雨量を想定し,この単位有効雨量から生じる直接流出ハイドログラフを単位図とする.

要は,1という量の雨が降った時の流出ハイドログラフ(流出量の時系列グラフ)を単位図と呼ぶということです.

2.不変性の仮定 流域が開発(ゴルフ場化,宅地化など)されるなど流域に変化が起きない限りは単位図の形状は不変と仮定する.

3.比例仮定 雨の量がR倍になれば,直接流出量もR倍になる.

4.重ね合わせの仮定 ある連続した時間の降雨がある場合,単位時間ごとの有効雨量は独立であると仮定する. このとき全流出量(直接流出+基底流出)ハイドログラフは各々の有効降雨から生じるハイドログラフの和として求まると仮定する.

言葉による説明だけだとピンとこないので,図で見てみましょう! 単位図.png

いま1という量の雨(単位有効雨量\( R_{e}=1 \))があるとします.なお,単位時間t0=1とします. このような単位雨量があったときに,例えば,最初の1時間(t=0〜1)では流出量が2,次の1時間(t=1〜2)では流出量が3,最後の1時間(t=2〜3)では流出量が1とします.この図のことを単位図と呼び,各時間の流出量の値を単位図成分や単位図流量と呼びます. そして,流域に変化がない限りこの波形(2,3,1)は不変です. これが,単位図の定義不変性の仮定です.

ちなみに単位図成分\( u(t) \)は以下のようになります.

\[ u(1)=2\\ u(2)=3\\ u(3)=1 \]

※ 1の雨が6の流出に変わっていることに違和感を覚えた方もいるかもしれませんが,これは単位有効雨量であることをイメージしやすくするためにあえてこうしています.どうしても違和感がある方は流域面積を6と考え,流域面積と有効雨量の積である流出量が6になるとでも考えておいてください.

比例仮定は,割と簡単にイメージできると思うので特に補足説明はありません.

重ね合わせの仮定は,ざっくり言えば各時間の雨に対して先ほどの単位図を書き,それを重ね合わせたものが直接流出量になるということです. 流出量計算ハイドログラフ_例.png

いま,5時間分の雨があったとします. 各時間の有効雨量Re(t)は次の通りとします.

\[ Re(1)=2,Re(2)=3,Re(3)=4,Re(4)=3,Re(5)=1 \]

まず,1時間目の有効雨量Re(1)は2なので,比例仮定より単位図の値を2倍にしたグラフを書きます(図中の赤色のカラム). 1時間目の有効雨量によって生じる各時間tの直接流出量Q1(t)は次式のようになります.

\[ \begin{align} Q1(1)&=Re(1) \times u(1)=2 \times 2=4\\ Q1(2)&=Re(1) \times u(2)=2 \times 3=6\\ Q1(3)&=Re(1) \times u(3)=2 \times 1=2\\ \end{align} \]

直接流出量をQとすると最初の1時間で流出するのはQ1(1)のみなので,

\[ Q(1)=Q1(1)=4 \]

次に2時間目の有効雨量Re(2)は3なので,比例仮定より単位図の値を3倍にしたグラフを書きます(図中の黄色のカラム). 2時間目の有効雨量によって生じる各時間tの直接流出量Q2(t)は次式のようになります.

\[ \begin{align} Q2(1)&=Re(2) \times u(1)=3 \times 2=6\\ Q2(2)&=Re(2) \times u(2)=3 \times 3=9\\ Q2(3)&=Re(2) \times u(3)=3 \times 1=3\\ \end{align} \]

2時間目に流出するのはQ1(2)とQ2(1)なので,

\[ Q(2)=Q1(2)+Q2(1)=6+6=12 \]

このように各時間の雨を独立と考えて計算し,足し合わせたものを全流出量とするのが重ね合わせの仮定です.

同様に,3時間目の有効雨量Re(3)は4なので,比例仮定より単位図の値を4倍にしたグラフを書きます(図中の黄緑色のカラム). 3時間目の有効雨量によって生じる各時間tの直接流出量Q3(t)は次式のようになります.

\[ \begin{align} Q3(1)&=Re(3) \times u(1)=4 \times 2=8\\ Q3(2)&=Re(3) \times u(2)=4 \times 3=12\\ Q3(3)&=Re(3) \times u(3)=4 \times 1=4\\ \end{align} \]

2時間目に流出するのはQ1(3)とQ2(2)とQ3(1)なので,

\[ Q(3)=Q1(3)+Q2(2)+Q3(1)=2+9+8=19 \]

これをまとめると

\[ \begin{align} Q(1)&=Re(1) \times u(1)\\ Q(2)&=Re(2) \times u(1)+Re(1) \times u(2)\\ Q(3)&=Re(3) \times u(1)+Re(2) \times u(2)+Re(1) \times u(3)\\ : \end{align} \]

これを一般化すると,時間ステップnの直接流出量Q(n)は次式で表されます.

\[ \begin{align} Q(n)&=Re(n) \times u(1)+Re(n-1) \times u(2)+・・・+Re(2) \times u(n-1)+Re(1) \times u(n)\\ &=\sum _{j=1}^{n} Re(n-(j-1)) \times u(j) \end{align} \]

こうして見ると,一見よくわからない単位図の数列式も,ブロック(各時間の雨に対する流出量)の積み上げを表しているに過ぎないことがわかると思います.

3.単位図法の実際

3.1.単位図の作成方法

ここまでは,単位図成分の値や数がわかっている前提で話を進めてきました. しかし,現実には単位図が初めからわかっている場合はほぼありません. これまで単位有効雨量(Re=1)について話を進めてきましたが,実際の流域できれいな単位雨量が観測されるとも限りません.

そこで,流量データのみを使って単位図を求める方法を紹介します.

まず,単位時間t0=2(h),単位有効雨量Re0=10(mm)とします. 流域全体に一様に降雨が見られた時の流量データを用意します. 例えば,以下のような流量データが得られたとします.この流量Qは全流出量に相当します. なお,流域面積Aをここでは500(km2)とします.

時刻 t (h) 流量 Q(t) (m3/s)
0 0.8
2 13.0
4 43.8
5 62.0
8 67.5
10 46.0
12 33.8
14 22.0
16 14.5
18 9.3
20 4.6
22 1.8

下図のように,洪水流出の始まりと終わりの流量を結んだ赤色の直線の値を各時刻tの流出量の値を基底流出量Q_base(t)とします. 基底流出量.png

直接流出量Q_directおよび基底流出量Q_baseは次の通りです.

時刻 t (h) 流量 Q(t) (m3/s) 直接流出量 Q_direct(t) (m3/s) 基底流出量 Q_base(t) (m3/s)
0 0.8 0.8 0.0
2 13.0 0.9 12.1
4 43.8 1.0 42.8
6 62.0 1.1 60.9
8 67.5 1.2 66.3
10 46.0 1.3 44.7
12 33.8 1.3 32.5
14 22.0 1.4 20.6
16 14.5 1.5 13.0
18 9.3 1.6 7.7
20 4.6 1.7 2.9
22 1.8 1.8 0.0

直接流出量の総和ΣQ_directを求めます.

\[ \Sigma Q \_ direct = \sum Q\_direct(t)=303.5 \]

流域面積A(km2)の流域で単位時間t0(h)に単位有効雨量Re0(mm)の降雨があったとき,流域からの直接流出量の総和ΣQ_unit(m3/s)は次式で表されます.

\[ \Sigma Q \_ unit = \left( \dfrac{Re0}{1000} \right) \times \left( \dfrac{1}{3600 \times t0} \right) \times \left(A \times 10 ^{6} \right) \]

今回は例として流域面積Aを500 (km2),単位時間t0を2 (h),単位有効雨量Re0を10 (mm)としているので,

\[ \begin{align} \Sigma Q \_ unit &= \left( \dfrac{10}{1000} \right) \times \left( \dfrac{1}{3600 \times 2} \right) \times \left(500 \times 10 ^{6} \right)\\ &\fallingdotseq 694 \end{align} \]

単位図成分u(t),単位有効雨量に対する直接流出量の総和ΣQ_unit,観測された直接流出量Q_direct(t),観測された直接流出量の総和ΣQ_directの関係は次式のようになります.

\[ u(t) : Q\_direct(t)= \Sigma Q \_ unit : \Sigma Q \_ direct \]

したがって,単位図成分u(t)は

\[ u(t) = \dfrac{\Sigma Q \_ unit} {\Sigma Q \_ direct} \times Q\_direct(t) \]
時刻 t (h) 単位図成分 u(t) (m3/s)
0 0.00
2 27.69
4 97.93
6 139.35
8 151.70
10 102.28
12 74.36
14 47.14
16 29.75
18 17.62
20 6.64
22 0.00

3.2.単位図を用いた流出量の計算

いま,上記の単位図を持つ流域に下表に示す降雨イベントがあったとします.

時刻 t (h) 流出係数f 降雨強度R (mm/h)
2 0.40 56.2
4 0.50 70.5
6 0.80 40.4

単位図法を使って計算すると,この雨に対する全流出量は下図のようになります. ただし,基底流出量は0.8 (m3/s)としています.

流出量計算ハイドログラフ.png

ちなみに計算用ワークシートはこんな感じです↓ 計算開始ボタンを押したら計算が開始されます. 単位図法.xlsm - Excel 2021_01_27 10_31_22.png

4.おわりに

最後までお付き合い頂き,ありがとうございました.

筆者は単位時間は1時間とばかり思っていましたが,今回参考にした書籍では2時間としていました. なぜ2時間にしているのかご存じの方は教えてください.

今回の記事の内容はいかがだったでしょうか? 本記事が水文学やかんがい排水など勉強を支える一助になれば幸いです.

本記事の執筆にあたり,書籍の内容について紹介することを快く承諾いただいた森北出版の皆様と,何より書籍の著者である河村三郎先生には心より感謝いたします.

それでは,また次回の記事でお会いしましょう.

今後も農業農村工学(水文学,かんがい排水,土壌物理,水理学)を中心に記事を執筆していきたいと思います. リクエスト等も受け付けておりますので,ご遠慮なく連絡ください. Twitterアカウント:エビぐんかん@6LxAi9GCOmRigUI メール:nnCreatorCircle@gmail.com

検索用Key Words: 単位図法,単位図成分,単位図,流出解析,Unit Hydrograph,Sherman,1932,降雨,降水,水文事象,水文学,計算方法,解き方,プログラム,マクロ,やり方,使い方

引用・参考文献

(1) 河村三郎:土木工学プログラム集1 水文・水理1 実用プログラムとグラフィック表示,森北出版,1984,pp.28-37. (2) 丸山利輔,石川重雄,大槻恭一,高瀬恵次,永井明博,田中丸治哉,駒村正治,赤江剛夫,堀野治彦,三野徹,武田育郎,金木亮一,渡辺紹裕:地域環境水文学,朝倉書店,2014,pp.63-64,77.

(注意)グラフの横軸は体裁をよくするためテキストボックスを上から重ねています.入力データを変えるとズレるので適宜編集してご利用ください.

今回使用した計算用のワークシートは下記URLからダウンロードできます. ご自由にお使いください. 単位図.xlsm

【流出解析】シリーズ 一覧

【流出解析】シリーズ 第1回 流域界の決定と入力データの用意 【流出解析】シリーズ 第2回 貯留関数法【前編】 【流出解析】シリーズ 第3回 貯留関数法【後編】 【流出解析】シリーズ 第4回 流域実蒸発散量の推定―蒸発散比法― 【流出解析】シリーズ 第5回 実蒸発散量の推定―補完法― 【流出解析】シリーズ 第6回 タンクモデル 【流出解析】シリーズ 第7回 積雪・融雪モデル ―降雪量・融雪量の推定― 【流出解析】シリーズ 第8回 タンクモデルの最適化 【流出解析】シリーズ 第9回 単位図法 【流出解析】シリーズ 第10回 流出量・流出高・比流量の違い 【流出解析】シリーズ 第11回 表面流モデル【定常降雨:前編】 【流出解析】シリーズ 第12回 表面流モデル【定常降雨:後編】

5.計算用プログラム

計算用プログラムです.ご自由にお使いください. ちなみに,わざわざプログラムを書く必要はなくエクセル関数でも単位図法の実装は可能です.

Option Explicit
Sub UnitHydrographMethod()
'プログラムの内容:単位図法によって流出量の求める
'作成年月日:20210124
'最終更新日時:20210125

'入力データ:降水量P,流出係数f,流出量Q

'画面更新の停止
Application.ScreenUpdating = False
'計算開始メッセージ
Debug.Print "計算開始 "; Now()
'ワークシートのアクティブ化
Worksheets("単位図法").Activate

'エラー防止
On Error Goto ErrorFlag

'************変数の定義************
Dim i As Long, j As Long, Cnt As Long     'カウンタ
Dim A As Double                 '流域面積(km2)
Dim CntQ As Long                '流出量データの数-1を格納する
Dim CntR As Long                '降水量データの数-1を格納する
Dim t() As Integer              '時刻を格納する配列(要素番号は0からスタート)
Dim Q() As Double               '流出量を格納する配列(要素番号は0からスタート)
Dim f() As Double               '流出係数を格納する配列(要素番号は0からスタート)
Dim R() As Double               '降雨量を格納する配列(要素番号は0からスタート)
Dim Re() As Double              '有効雨量を格納する配列(要素番号は0からスタート)
Dim t0 As Integer               '単位時間t0(h)
Dim R0 As Double                '単位雨量R0(mm)
Dim Re_u() As Double            '流出量の成分Re(j)×Q_unit(i).多次元配列(行:単位図成分数,列:降雨量データ数)
Dim t_runoff_obs As Integer     '単位図作成のために用意した流出量データの総時間(h) t_runoff_obs=CntQ*t0
Dim t_runoff_cal As Integer     '単位図と雨量データから計算される流出量データの総時間(h) t_runoff_cal=(m+k-1)*t0
Dim Q_unit() As Double          '単位図成分(m3/s)
Dim Q_base() As Double          '基底流出量(m3/s)
Dim Q_direct() As Double        '直接流出量(m3/s)
Dim sigma_Q_unit As Double      '単位図成分の総和(m3)
Dim sigma_Q_direct As Double    '直接流出量の総和(m3)
Dim Ratio_Unit_Obs As Double    '単位図成分の総和と直接流出量の総和の比


'************変数の初期化************
A = 500
t0 = 2
R0 = 10
sigma_Q_unit = 0
sigma_Q_direct = 0


'************結果出力用ラベルの作成************
Cells(1, 6) = "時間 t (h)"
Cells(1, 7) = "基底流出量Q_base (m3/s)"
Cells(1, 8) = "直接流出量Q_direct (m3/s)"
Cells(1, 9) = "単位図成分Q_unit (m3/s)"
Rows("1:1").WrapText = True


'************入力データの取得************
'流出量データおよび降水量データの数を取得
i = 2
Do
    If Cells(i, 2) <> "" And Cells(i + 1, 2) = "" Then
        CntQ = i
    End If

    '降水量データは(時間ステップ)2時間目以降からカウント開始
    If i >= 3 And Cells(i, 4) <> "" And Cells(i + 1, 4) = "" Then
        CntR = i
    End If
    i = i + 1
Loop Until Cells(i, 1) = ""

'配列の要素数を決定
CntQ = CntQ - 2
CntR = CntR - 2

'単位図作成のために用意した流出量データの総時間を計算
t_runoff_obs = CntQ * t0

'配列の再定義
ReDim Q(0 To CntQ), Q_direct(0 To CntQ), Q_base(0 To CntQ), t(0 To CntQ), Q_unit(0 To CntQ)
ReDim f(0 To CntR), R(0 To CntR), Re(0 To CntR)
ReDim Re_u(0 To CntQ, 0 To CntR)

'入力データを配列に格納
i = 2
Do
    t(i - 2) = Cells(i, 1)
    If i - 2 <= CntQ Then
        Q(i - 2) = Cells(i, 2)
    End If

    '降水量データは(時間ステップ)2時間目以降からカウント開始
    If i - 2 <= CntR Then
        f(i - 2) = Cells(i, 3)
        R(i - 2) = Cells(i, 4)
    End If
    i = i + 1
Loop Until Cells(i, 1) = ""


'************************第1部 単位図成分を求める************************
'************基底流出量,直接流出量,直接流出量の総和を求める************
'流出の始まりQ(0)と終わりQ(CntQ)を直線で結んだ直線の各時刻tの値を基底流出量Q_base(t)とする
For i = 0 To CntQ
    Q_base(i) = Q(0) + Int((Q(CntQ) - Q(0)) / t_runoff_obs * t(i) * 10 + 0.5) / 10  '基底流出量の変化量の小数点第二位を四捨五入
    Q_direct(i) = Abs(Q(i) - Q_base(i))                                             '直接流出量は非負とする
    sigma_Q_direct = sigma_Q_direct + Q_direct(i)
Next i


'************単位図成分の総和と単位図成分を求める************
'単位図成分の総和
sigma_Q_unit = (R0 / 1000) / (3600 * t0) * (A * 10 ^ 6)

'単位図成分の総和と直接流出量の総和の比を求める
Ratio_Unit_Obs = sigma_Q_unit / sigma_Q_direct

'各時刻tの単位図成分Q_unit(t)を求める
For i = 0 To CntQ
    Q_unit(i) = Int(Q_direct(i) * Ratio_Unit_Obs * 100 + 0.5) / 100     '基底流出量の変化量の小数点第三位を四捨五入
Next i


'************計算結果の出力************
For i = 0 To CntQ
    Cells(i + 2, 6) = t(i)
    Cells(i + 2, 7) = Q_base(i)
    Cells(i + 2, 8) = Q_direct(i)
    Cells(i + 2, 9) = Q_unit(i)
Next i



'************************第2部 降雨量と流出係数のデータから流出量ハイドログラフを求める************************
'************有効雨量を求める************
For i = 0 To CntR
    Re(i) = R(i) * f(i)     '小数点第三位を四捨五入
Next i

'************流出量を求める************
t_runoff_cal = (CntQ + 1) + CntR - 1
ReDim Q_direct(0 To t_runoff_cal)       '配列Q_directは使い回し

'配列の初期化
For i = 0 To t_runoff_cal
    Q_direct(i) = 0
Next i

For j = 1 To CntR
    For i = 0 To CntQ
        Q_direct(i + j - 1) = Q_direct(i + j - 1) + Int(Re(j) / R0 * Q_unit(i) * 100 + 0.5) / 100 '小数点第三位を四捨五入
        Re_u(i, j) = Int(Re(j) / R0 * Q_unit(i) * 100 + 0.5) / 100
    Next i
Next j

'************計算結果の出力************
'結果出力用ラベルの作成
Cells(1, 10) = "時刻 t(h)"
Cells(1, 11) = "有効降雨強度Re(mm)"
Cells(1, 12) = "直接流出量Q_direct(m3/s)"
Cells(1, CntR + 13) = "時刻 t(h)"
Cells(1, CntR + 14) = "有効降雨強度Re(mm)"
Cells(1, CntR + 15) = "単位図成分Q_unit(m3/s)"
Cells(1, CntR + 16) = "全流出量Q(m3/s)"

'直接流出量
For i = 0 To t_runoff_cal
    Cells(i + 2, 10) = i * 2
    Cells(i + 2, 12) = Q_direct(i)
Next i

'有効雨量と直接流出量
For j = 1 To CntR
    'ラベル作成
    Cells(1, j + 12) = "Re_" & j & "×u(t)"

    Cells(j + 2, 11) = Re(j)
    For i = 0 To CntQ
        Cells(i + j + 1, j + 12) = Re_u(i, j)
    Next i
Next j

'************時刻t=0の流出量を基底流出量として,全流出量を求める************
'線形補間して1時間単位の全流出量を求める
Cnt = 2
For i = 0 To t_runoff_cal * t0
    Cells(Cnt, CntR + 13) = Cnt - 2
    '時刻が偶数の時
    If i Mod 2 = 0 Then
        '有効雨量
        If i / 2 <= CntR Then
            Cells(Cnt, CntR + 14) = Re(i / 2)
        Else
            Cells(Cnt, CntR + 14) = 0
        End If
        '単位図成分
        If i / 2 <= CntQ Then
            Cells(Cnt, CntR + 15) = Q_unit(i / 2)
        Else
            Cells(Cnt, CntR + 15) = 0
        End If
        '単位図成分
        Cells(Cnt, CntR + 16) = Q_direct(i / 2) + Q(0)


    '時刻が奇数の時
    ElseIf i Mod 2 <> 0 Then
        '有効雨量
        If i / 2 <= CntR Then
            Cells(Cnt, CntR + 14) = (Re((i - 1) / 2) + Re((i + 1) / 2)) / 2
        Else
            Cells(Cnt, CntR + 14) = 0
        End If
        '単位図成分
        If i / 2 <= CntQ Then
            Cells(Cnt, CntR + 15) = (Q_unit((i - 1) / 2) + Q_unit((i + 1) / 2)) / 2
        Else
            Cells(Cnt, CntR + 15) = 0
        End If
        '単位図成分
        Cells(Cnt, CntR + 16) = (Q_direct((i - 1) / 2) + Q_direct((i + 1) / 2)) / 2 + Q(0)
    End If
    Cnt = Cnt + 1
Next i

'画面更新の再開
Application.ScreenUpdating = True
'計算終了メッセージ
Debug.Print "計算終了 "; Now()
MsgBox ("計算が終了しました.")
Exit Sub

'エラー防止
ErrorFlag:
    MsgBox ("エラーが発生しました")
End Sub