以前のリビジョンの文書です


tys-sfm_v10

tys-sfm_v10.lzhをダウンロードして解凍すると、「貯留関数マクロ.xls」ができます。 この「貯留関数マクロ.xls」の使用方法を説明します。 ダウンロードについては、水理水文ソフトウェアのHPを参照してください。

マクロの説明

本マクロは、各条件シートに条件を入力するだけで、1流域についての貯留関数による流域流出解析計算が行えます。

  1. 実績流量より、貯留関数係数の回帰を行う。
  2. 任意の貯留関数係数を用いて流出解析を行う。

の2種類の計算を行うことができます。

計算条件とデータの入力

計算条件

計算条件」シートで貯留関数係数の回帰を行うかどうかをチェックしてください。

計算条件シート

計算条件シート

再度計算実行した場合は、計算結果シート(複数)の値がすべて上書きされます。

雨量・流量データ

雨量・流量データシート

雨量・流量データシート


  1. 実績流量より、貯留関数係数の回帰を行う場合 「計算条件」シートの流域面積・基底流量・流路長さ・遅滞時間計算ピッチを入力してください。 「雨量・流量データ」シートの流域平均雨量・実績流量・計算対象範囲指定のデータを入力してください。計算対象範囲指定は入力した実績流量のうち、係数決定に用いるデータの範囲を指定します。開始点を「1」、終了点を「2」と入力してください。
  2. 任意の貯留関数係数を用いて流出解析を行う場合計算条件」シートの流域面積・基底流量・貯留関数係数(K・P値、遅滞時間、流出係数)を入力してください。「雨量・流量データ」シートの流域平均雨量を入力してください。実績流量は、計算には使用しませんが、雨量・流量結果図で計算結果の評価を行うために入力する方が望ましいでしょう。

計算結果

  1. 実績流量より、貯留関数係数の回帰を行う場合
    1. 計算結果」シートに遅滞時間の設定値ごとの貯留量と流出量の関係の計算結果を書き込みます。この中で、貯留量と流出量のループの積算値の最小値を遅滞時間として採用します。尚、遅滞時間の繰り返し計算は、0時間~最大値(流路長を木村の式に適用した場合の遅滞時間+1)を遅滞時間の計算ピッチ分です。ここで求めた遅滞時間より実績流量を用いた最小自乗法による回帰を行い、K・Pを算出します。
    2. 雨量・流量データ」シートの計算流量の列に上記で求めた貯留関数係数を用いて計算した流量を書き込みます。計算結果は「雨量・流量結果図」に反映されます。
    3. 貯留量と流出量の関係」シートに決定した貯留関数係数を用いた貯留量と流出量の関係を出力します。その結果は「貯留量と流出量の関係図」シートに反映されます。
  2. 任意の貯留関数係数を用いて流出解析を行う場合
    1. 雨量・流量データ」シートの計算流量の列に入力した貯留関数係数を用いて計算した流量を書き込みます。計算結果は「雨量・流量結果図」シートに反映されます。
    2. 貯留量と流出量の関係」シートに貯留量と流出量の関係を出力します。その結果は「貯留量と流出量の関係図」シートに反映されます。


計算結果

計算結果

雨量・流量結果図

雨量・流量結果図

貯留量と流出量の関係図

貯留量と流出量の関係図

注意事項

 本マクロは、通常の貯留関数法についての計算用であり、総合貯留関数法には適用していません。通常の貯留関数法と総合貯留関数法では、K・P値はそれぞれ全く異なりますので注意してください。  ここで、通常の貯留関数法とは、貯留関数を貯留量(m3)と流出量(m3/s)の関係で表現したもの、総合貯留関数法とは貯留関数を貯留高(mm)と流出高(mm/hr)の関係で表現したものです。総合貯留関数法は貯留関数法を単位流域面積で表現したものと言い換えることができます。  本や論文によっては、総合貯留関数法を単に「貯留関数法」と表記されている場合がありますので、基礎式による確認が必要です。

VBAソースファイル

 本マクロのソースを公開しますので、参考にしてください。

Option Explicit

Public m As Integer     'ハイドロデータ数
Public m1 As Integer    'ハイドロ区間1
Public m2 As Long       'ハイドロ区間2
Public l As Double      '流域の流路延長(km)
Public a As Double      '流域面積(km2)
Public qi As Double     '基底流量(m3/s)
'
Public Const mnum As Integer = 1000  '配列数
Public r(0 To mnum) As Double       '雨量(mm/h)
Public q(0 To mnum) As Double       '実測流量(m3/s)

Public fr(0 To mnum) As Double      '雨量を流量に換算した値
Public fi(0 To mnum) As Double      '雨量を流量に換算した値の積算値

Public ql(0 To mnum) As Double      '遅滞時間を考慮した流量(基底流量を差し引いたもの)
Public qs(0 To mnum) As Double      '遅滞時間を考慮した流量の積分値(基底流量を差し引いたもの)
Public s(0 To mnum) As Double
Public sa(0 To 10000) As Double     'ループの積算値(時計回りか反時計回りかを判定)
'
Public rt As Double         '対象範囲の累積雨量(mm)
Public f As Double          '流出係数
Public tlk As Double        '遅滞時間(設定値)
Public qlm As Double        '流量の最大値(基底流量を差し引いたもの)
Public tpl As Double        '流量が最大の時間(遅滞時間分差し引いたもの)
Public qtm As Double        '流量の積算(遅滞時間分差し引いたもの)
Public qsm As Double        'Qsの積算
Public fim As Double        'Fiの最大値
Public sam As Double        'ループの最大
Public sas As Double        'ループの最小
Public tld As Double        '遅滞時間(決定分)
Public tl0 As Integer       '遅滞時間(切り上げ)
Public tl1 As Integer       '遅滞時間(整数部)
Public tl2 As Double        '遅滞時間小数部
Public t_Step As Double     '遅滞時間の計算ステップ
Public num_Cal As Integer    '遅滞時間計算回数

'
Public e As Double          '最小自乗法の対象となるデータの個数
Public sm As Double         '貯留高さの最大値
Public sar As Double        '上昇ループの積算値
Public sad As Double        '下降ループの積算値
'
Public k As Double          'K値
Public p As Double          'P値

Public Sub Main()
'貯留関数法1流域のみ、係数を回帰で決定
    Dim ws1 As Worksheet
    Dim ws2 As Worksheet
    Dim ws3 As Worksheet
    Dim ws4 As Worksheet
    Dim i As Integer
    Dim tl As Integer
    Dim t As Integer
    Set ws1 = Sheets("計算条件")
    Set ws2 = Sheets("雨量・流量データ")
    Set ws3 = Sheets("計算結果")
    Set ws4 = Sheets("貯留量と流出量の関係")
    'データ読み込み
    DataRead1 ws1
    DataRead2 ws2
    
    If Sheet1.CheckBox1.Value = False Then GoTo CalSkip '係数回帰をするかどうかのチェック
    
    '遅滞時間の初期定義(木村の式)
    If l < 11.9 Then
        tl = 0
    Else
        tl = Int(0.5 + 0.047 * l - 0.56)   '切り上げ
    End If
    '適当な遅滞時間の判別
    '**********************
        num_Cal = CInt(CDbl(tl + 1) / t_Step)
    '**********************
    For t = 0 To num_Cal
        tlk = t_Step * CDbl(t)  '遅滞時間設定値
        tl1 = Int(tlk)          '遅滞時間整数部
        tl2 = tlk - tl1         '遅滞時間少数部
        If tl2 = 0 Then
            tl0 = tl1
        Else
            tl0 = tl1 + 1
        End If
    
        CalQl
        
        'DataWrite0 T   'デバック
        
        CalQtm
        f = qtm / ((a / 3.6) * rt)
        CalFr
        CalSa t
        '
        DataWrite1 ws3, t   '計算経過出力
        '
    Next t
    'ループの最大・最小判定
    sam = 0
    For t = 0 To num_Cal
        If Abs(sa(t)) >= sam Then sam = Abs(sa(t))
    Next t
    sas = sam
    For t = 0 To num_Cal
        If Abs(sa(t)) <= sas Then
            sas = Abs(sa(t))
            tld = t_Step * CDbl(t)  '遅滞時間設定値
        End If
    Next t
    tl1 = Int(tld)          '遅滞時間整数部
    tl2 = tld - tl1         '遅滞時間少数部
    If tl2 = 0 Then
        tl0 = tl1
    Else
        tl0 = tl1 + 1
    End If
    
    CalQl
    
    CalQtm
    
    f = qtm / ((a / 3.6) * rt)  '流出係数
    
    fim = 0
    For t = 1 To m
        fr(t) = a * f * r(t) / 3.6
        fim = fim + fr(t)
        fi(t) = fim
    Next t
    
    qsm = 0#
    For t = 2 To m - tl0
        qsm = qsm + (ql(t - 1) + ql(t)) / 2#
        qs(t) = qsm
        If fi(t) - qs(t) >= 0 Then s(t) = fi(t) - qs(t)
    Next t
    
    sm = 0
    For t = 2 To m - tl0 - 1
        If s(t) >= sm Then sm = s(t)
    Next t
    
    Cal_KP  '最小自乗法によるK・Pの計算
    
    DataWrite2 ws1      '貯留関数式書き込み
    DataWrite3 ws4      '貯留量と流出高書き込み
    Graph_Refresh       'グラフ修正
    
CalSkip:    '流量計算
    
    DataRead3 ws1
        
    Cal_Qt ws2          '求めた係数で流量を計算
        
    Graph_Refresh2       'グラフ修正
    
End Sub

Private Sub DataRead1(ws As Worksheet)
'データ読み込み1(流域諸元)
    ws.Select
    a = ws.Cells(3, "D")        '流域面積(km2)
    qi = ws.Cells(4, "D")       '基底流量(m3/s)
    l = ws.Cells(5, "D")        '流路長(km)
    t_Step = ws.Cells(6, "D")   '遅滞時間計算ピッチ(hr)
    If t_Step = 0 Then t_Step = 1
End Sub

Private Sub DataRead2(ws As Worksheet)
'データ読み込み1(ハイドロ・ハイエトデータ)
    Dim row_num As Integer
    Dim i As Integer
    '①ハイドロ・ハイエト読み込み
    row_num = 3
    i = 0
    Do Until ws.Cells(row_num, "C") = ""    'データがあるまでループ
        i = i + 1
        r(i) = ws.Cells(row_num, "C")
        q(i) = ws.Cells(row_num, "D")
        row_num = row_num + 1
    Loop
    m = i       'ハイドロ・ハイエトデータ数定義
    '計算対象範囲指定
    For i = 1 To m
        If ws.Cells(i + 2, "E") = 1 Then m1 = i         '開始点
        If ws.Cells(i + 2, "E") = 2 Then m2 = i         '終了点
    Next i
    '期間の累積雨量
    rt = 0
    For i = m1 To m2
        rt = rt + r(i)
    Next i
End Sub

Private Sub CalFr()
    Dim t As Integer
    fim = 0
    For t = m1 To m2
        fr(t) = a * f * r(t) / 3.6
        fim = fim + fr(t)
        fi(t) = fim
    Next t
End Sub

Private Sub CalQl()
'遅滞時間分戻した流量の定義(基底流量含まず)
    Dim t As Integer
    qlm = 0
    For t = 1 To m2
        If t - tl1 >= 0 Then
            ql(t - tl1) = q(t) + (q(t + 1) - q(t)) * tl2 - qi  '遅滞時間を考慮した流量(基底流量を差し引いたもの)
            If ql(t - tl1) >= qlm Then
                qlm = ql(t - tl1)    '最大流量
                tpl = t - tlk        '最大流量の発生時間
            End If
        End If
    Next t
End Sub

Private Sub CalQtm()
'遅滞時間を考慮した流量の積算(計算区間内)
    Dim t As Integer
    qtm = 0
    For t = m1 To m2 - tl0
        qtm = qtm + ql(t)
    Next t
End Sub

Private Sub CalSa(tl As Integer)
'貯留高さと流出高をプロットした場合のループ形状
    Dim t As Integer
    qsm = 0: sar = 0: sad = 0
    
    For t = m1 To m2 - tl0 - 1
        qsm = qsm + (ql(t - 1) + ql(t)) / 2#
        qs(t) = qsm
        If fi(t) - qs(t) < 0 Then
            s(t) = 0#
        Else
            s(t) = fi(t) - qs(t)
        End If
    Next t
    '上昇ループの積算値計算
    For t = m1 To tpl - 1
        sar = sar + (s(t) + s(t + 1)) * Abs(ql(t + 1) - ql(t)) / 2#
    Next t
    '下降ループの積算値計算
    For t = tpl To m2 - tl0 - 1
        sad = sad - (s(t) + s(t + 1)) * Abs(ql(t) - ql(t + 1)) / 2#
    Next t
    'ループ全体の積算値計算
    sa(tl) = sar + sad
End Sub

Private Sub Cal_KP()
'最小自乗法による貯留関数係数の計算
    Dim M3 As Integer
    Dim t As Integer
    Dim Xi As Double
    Dim Xx As Double
    Dim Yi As Double
    Dim XY As Double
    Dim P0 As Double
    Dim K0 As Double
    Dim K1 As Double
    'M3 = Int(0.8 * M)
    M3 = m2 - tl0
    e = 0
    'For T = 2 To M3
    For t = m1 To M3
        If s(t) <> 0 Then e = e + 1
    Next t
    Xi = 0: Xx = 0: Yi = 0: XY = 0
    'For T = 2 To M3
    For t = m1 To M3
        Xi = Xi + Log(ql(t))
        Xx = Xx + Log(ql(t)) * Log(ql(t))
        Yi = Yi + Log(s(t))
        XY = XY + Log(ql(t)) * Log(s(t))
    Next t
    P0 = (Xi * Yi - e * XY) / (Xi * Xi - e * Xx)
    K0 = (Yi - P0 * Xi) / e
    K1 = Exp(K0)
    k = K1
    p = P0
    k = Int(0.5 + K1 * 100#) / 100#     '桁合わせ
    p = Int(0.5 + P0 * 1000#) / 1000#   '桁合わせ
End Sub
Private Sub DataWrite1(ws As Worksheet, t As Integer)
'結果書き込み1(遅滞時間計算経過)
    Dim Sa_State As String
    If t = 0 Then
        ws.Cells(2, "A") = "遅滞時間(hr)"
        ws.Cells(2, "B") = "最大流量(m3/s)"
        ws.Cells(2, "C") = "最大流量時間(hr)"
        ws.Cells(2, "D") = "流量積算値(m3/s)"
        ws.Cells(2, "E") = "流出係数"
        ws.Cells(2, "F") = "上昇ループの積算値"
        ws.Cells(2, "G") = "下降ループの積算値"
        ws.Cells(2, "H") = "ループ全体の積算値"
        ws.Cells(2, "I") = "ループの状態"
    End If
    '
    Sa_State = ""
    If sa(t) > 0 Then Sa_State = "時計回り"
    If sa(t) < 0 Then Sa_State = "反時計回り"
    If sa(t) = 0 Then Sa_State = "ぴったりOK"
    '
    ws.Cells(t + 3, "A") = tlk
    ws.Cells(t + 3, "B") = qlm
    ws.Cells(t + 3, "C") = tpl
    ws.Cells(t + 3, "D") = qtm
    ws.Cells(t + 3, "E") = f
    ws.Cells(t + 3, "F") = sar
    ws.Cells(t + 3, "G") = sad
    ws.Cells(t + 3, "H") = sa(t)
    ws.Cells(t + 3, "I") = Sa_State
    
End Sub

Private Sub DataWrite2(ws As Worksheet)
'結果書き込み2(貯留関数式出力)
    ws.Select
    ws.Cells(10, "B") = k
    ws.Cells(11, "B") = p
    ws.Cells(12, "B") = tld
    ws.Cells(13, "B") = f
End Sub

Private Sub DataWrite3(ws As Worksheet)
'貯留量と流出高書き込み
    Dim t As Integer
    ws.Select
    ws.Cells(3, "D") = tld
    For t = 1 To m - tl0
        ws.Cells(t + 2, "A") = ql(t)
        ws.Cells(t + 2, "B") = s(t)
    Next t
End Sub

Private Sub DataRead3(ws As Worksheet)
'データ読み込み3(貯留関数式元)
    ws.Select
    k = ws.Cells(10, "B")
    p = ws.Cells(11, "B")
    tld = ws.Cells(12, "B")
    f = ws.Cells(13, "B")
End Sub

Private Sub Cal_Qt(ws As Worksheet)
'求めた貯留関数を用いて流量を再度計算
    Dim t As Integer
    Dim tm As Integer
    Dim mr As Integer   '雨量・流量データ数
    Dim qql(0 To mnum) As Double
    Dim qqt(0 To mnum) As Double
    Dim j As Integer
    Dim c As Double
    Dim d As Double
    Dim x As Double
    Dim x2 As Double
    Dim y As Double
    mr = m
    tl1 = Int(tld)          '遅滞時間整数部
    tl2 = tld - tl1         '遅滞時間少数部
    'If Tl2 <> 0 Then Tl1 = Tl1 + 1
     'ニュートン法による流量計算
    qql(1) = 0
    For t = 1 To mr
        c = 2# * f * a * r(t) / 3.6 - qql(t) + 2# * k * qql(t) ^ p
        x = 1
        j = j + 1
Step1:  y = x + (2# * k * x ^ p) - c
        d = 1# + 2# * p * k * x ^ (p - 1#)
        x2 = x - y / d
        If Abs(x2 - x) <= 0.0001 Then GoTo Step2
        x = x2
        GoTo Step1
Step2:  qql(t + 1) = x
        If t <= tl1 Then qqt(t) = qi
        qqt(t + tl1) = qi + qql(t + 1) + (qql(t + 1) - qql(t)) * tl2
        If j <= 100 Then GoTo Step3
        If qqt(t) < 1.1 * qqt(99) Then GoTo Skip1
Step3:
    Next t
Skip1:
    For t = 1 To mr
        ws.Cells(t + 2, "G") = qqt(t)
    Next t
End Sub

Private Sub Graph_Refresh()
'グラフの修正
    Dim g_range As String
    Sheets("貯留量と流出量の関係図").Select
    g_range = "A3:B" & CStr(m - tl0 + 2)
'
    ActiveChart.PlotArea.Select
    ActiveChart.SetSourceData Source:=Sheets("貯留量と流出量の関係").Range(g_range), PlotBy:=xlColumns

End Sub

Private Sub Graph_Refresh2()
'グラフの修正
    Dim g_range As String
    Sheets("雨量・流量結果図").Select
    g_range = "B2:D" & CStr(m + 2) & ",F2:G" & CStr(m + 2)
'
    ActiveChart.PlotArea.Select
    ActiveChart.SetSourceData Source:=Sheets("雨量・流量データ").Range(g_range), PlotBy:=xlColumns
End Sub

Private Sub DataWrite0(t As Integer)
'結果書き込み(デバック)
    Dim ws As Worksheet
    Dim i As Integer
    Set ws = Sheets("Ql")
    ws.Cells(2, t + 2) = t
    ws.Cells(3, t + 2) = tld
    ws.Cells(4, t + 2) = tl1
    ws.Cells(5, t + 2) = tl2
    For i = 1 To m2 - tl0
        'If I - Tl1 >= 0 Then
        '    ws.Cells(I + 5, T + 2) = Ql(I - Tl1)
        'End If
        ws.Cells(i + 5, t + 2) = ql(i)
    Next i
End Sub

tys-sfm_v_2_0_1

tys-sfm_v_2_0_1.lzhをダウンロードして解凍すると、「tys-sfm_v_2_0_1.xls」ができます。 この「tys-sfm_v_2_0_1.xls」の使用方法を説明します。 ダウンロードについては、水理水文ソフトウェアのHPを参照してください。

マクロの説明

本マクロは、各条件シートに条件を入力するだけで、1流域についての総合貯留関数法と貯留関数法の2種類の流域流出解析計算が行えます。

  1. 実績流量より、貯留関数係数の回帰を行う。
  2. 任意の貯留関数係数を用いて流出解析を行う。

の2種類の計算を行うことができます。

計算条件とデータの入力

計算条件

計算条件」シートで貯留関数係数の回帰を行うかどうかをチェックしてください。

計算条件シート

計算条件シート

雨量・流量データ

雨量・流量データシート

雨量・流量データシート


  1. 実績流量より、貯留関数係数の回帰を行う場合 「計算条件」シートの流域面積・基底流量・流路長さ・遅滞時間計算ピッチを入力してください。 「雨量・流量データ」シートの流域平均雨量・実績流量・計算対象範囲指定のデータを入力してください。計算対象範囲指定は入力した実績流量のうち、係数決定に用いるデータの範囲を指定します。開始点を「1」、終了点を「2」と入力してください。
  2. 任意の貯留関数係数を用いて流出解析を行う場合計算条件」シートの流域面積・基底流量・貯留関数係数(K・P値、遅滞時間、流出係数)を入力してください。「雨量・流量データ」シートの流域平均雨量を入力してください。実績流量は、計算には使用しませんが、雨量・流量結果図で計算結果の評価を行うために入力する方が望ましいでしょう。

計算結果

  1. 実績流量より、貯留関数係数の回帰を行う場合
    1. 計算結果」シートに遅滞時間の設定値ごとの貯留量と流出量の関係の計算結果を書き込みます。この中で、貯留量と流出量のループの積算値の最小値を遅滞時間として採用します。尚、遅滞時間の繰り返し計算は、0時間~最大値(流路長を木村の式に適用した場合の遅滞時間+1)を遅滞時間の計算ピッチ分です。ここで求めた遅滞時間より実績流量を用いた最小自乗法による回帰を行い、K・Pを算出します。
    2. 雨量・流量データ」シートの計算流量の列に上記で求めた貯留関数係数を用いて計算した流量を書き込みます。計算結果は「雨量・流量結果図」に反映されます。
    3. 貯留量と流出量の関係」シートに決定した貯留関数係数を用いた貯留量と流出量の関係を出力します。その結果は「貯留量と流出量の関係図」シートに反映されます。
  2. 任意の貯留関数係数を用いて流出解析を行う場合
    1. 雨量・流量データ」シートの計算流量の列に入力した貯留関数係数を用いて計算した流量を書き込みます。計算結果は「雨量・流量結果図」シートに反映されます。
    2. 貯留量と流出量の関係」シートに貯留量と流出量の関係を出力します。その結果は「貯留量と流出量の関係図」シートに反映されます。


計算結果

計算結果

雨量・流量結果図

雨量・流量結果図

貯留量と流出量の関係図

貯留量と流出量の関係

貯留量と流出量の関係図

貯留量と流出量の関係図

注意事項

 本プログラムは、総合貯留関数法と通常の貯留関数法の両方が行えますが、総合貯留関数法と通常の貯留関数法では、K・P値はそれぞれ全く異なりますので注意してください。  ここで、総合貯留関数法とは貯留関数を貯留高(mm)と流出高(mm/hr)の関係で表現したもの、通常の貯留関数法とは、貯留関数を貯留量(m3)と流出量(m3/s)の関係で表現したものです。総合貯留関数法は貯留関数法を単位流域面積で表現したものと言い換えることができます。  本や論文によっては、総合貯留関数法を単に「貯留関数法」と表記されている場合があり、特殊な手法ではなく、かなり一般的に利用されていますので、基礎式による確認が必要です。

VBAソースファイル

 本マクロのソースを公開しますので、参考にしてください。

Option Explicit

Public m As Integer     'ハイドロデータ数
Public m1 As Integer    'ハイドロ区間1
Public m2 As Long       'ハイドロ区間2
Public l As Double      '流域の流路延長(km)
Public a As Double      '流域面積(km2)
Public qi As Double     '基底流量(m3/s)
'
Public Const mnum As Integer = 1000 '配列数
Public r(0 To mnum) As Double       '雨量(mm/h)
Public q(0 To mnum) As Double       '実測流量(m3/s)
Public fr(0 To mnum) As Double      '雨量を流量に換算した値
Public fi(0 To mnum) As Double      '雨量を流量に換算した値の積算値
Public ql(0 To mnum) As Double      '遅滞時間を考慮した流量(基底流量を差し引いたもの)
Public qs(0 To mnum) As Double      '遅滞時間を考慮した流量の積分値(基底流量を差し引いたもの)
Public s(0 To mnum) As Double
Public sa(0 To 10000) As Double     'ループの積算値(時計回りか反時計回りかを判定)
'
Public rt As Double         '対象範囲の累積雨量(mm)
Public f As Double          '流出係数
Public tlk As Double        '遅滞時間(設定値)
Public qlm As Double        '流量の最大値(基底流量を差し引いたもの)
Public tpl As Double        '流量が最大の時間(遅滞時間分差し引いたもの)
Public qtm As Double        '流量の積算(遅滞時間分差し引いたもの)
Public qsm As Double        'Qsの積算
Public fim As Double        'Fiの最大値
Public sam As Double        'ループの最大
Public sas As Double        'ループの最小
Public tld As Double        '遅滞時間(決定分)
Public tl0 As Integer       '遅滞時間(切り上げ)
Public tl1 As Integer       '遅滞時間(整数部)
Public tl2 As Double        '遅滞時間小数部
Public t_step As Double     '遅滞時間の計算ステップ
Public num_cal As Integer    '遅滞時間計算回数
'
Public e As Double          '最小自乗法の対象となるデータの個数
Public sm As Double         '貯留高さの最大値
Public sar As Double        '上昇ループの積算値
Public sad As Double        '下降ループの積算値
Public k As Double          'K値
Public p As Double          'P値
'
Public cal_flag As Boolean  'Trueであれば総合貯留関数法

Public Sub Main()
'貯留関数法1流域のみ、係数を回帰で決定
    Dim ws1 As Worksheet
    Dim ws2 As Worksheet
    Dim ws3 As Worksheet
    Dim i As Integer
    Dim tl As Integer
    Dim t As Integer
    
    'シートのセット
    Set ws1 = Sheets("計算条件")
    Set ws2 = Sheets("雨量・流量データ")
    Set ws3 = Sheets("計算結果")
    
    'データ読み込み
    data_read1 ws1
    data_read2 ws2
    
    If Sheet1.CheckBox1.Value = True Then
        '係数回帰をするかどうかのチェック
        If Sheet1.OptionButton1.Value = True Then
            '総合貯留関数法
            cal_flag = True
        Else
            '貯留関数法
            cal_flag = False
        End If
    
        '遅滞時間の初期定義(木村の式)
        If l < 11.9 Then
            tl = 0
        Else
            tl = Int(0.5 + 0.047 * l - 0.56)   '切り上げ
        End If
    
        '適当な遅滞時間の判別
        '**********************
        num_cal = CInt(CDbl(tl + 1) / t_step)
        '**********************
        For t = 0 To num_cal
            tlk = t_step * CDbl(t)  '遅滞時間設定値
            tl1 = Int(tlk)          '遅滞時間整数部
            tl2 = tlk - tl1         '遅滞時間少数部
            If tl2 = 0 Then
                tl0 = tl1
            Else
                tl0 = tl1 + 1
            End If
    
            cal_ql
            cal_qtm
        
            'data_write0 t   'デバック
            
            If cal_flag = True Then
                '総合貯留関数法
                f = qtm / rt
            Else
                '貯留関数法
                f = qtm / ((a / 3.6) * rt)
            End If
        
            cal_fr
        
            cal_sa t
        
            data_write1 ws3, t   '計算経過出力
        '
        Next t
        
        'ループの最大・最小判定
        sam = 0
        For t = 0 To num_cal
            If Abs(sa(t)) >= sam Then sam = Abs(sa(t))
        Next t
        sas = sam
        For t = 0 To num_cal
            If Abs(sa(t)) <= sas Then
                sas = Abs(sa(t))
                tld = t_step * CDbl(t)  '遅滞時間設定値
            End If
        Next t
        tl1 = Int(tld)          '遅滞時間整数部
        tl2 = tld - tl1         '遅滞時間少数部
        If tl2 = 0 Then
            tl0 = tl1
        Else
            tl0 = tl1 + 1
        End If
    
        cal_ql
    
        cal_qtm
    
        If cal_flag = True Then
            '総合貯留関数法
            f = qtm / rt                '流出係数
        Else
            '貯留関数法
            f = qtm / ((a / 3.6) * rt)  '流出係数
        End If
    
        fim = 0
        For t = 1 To m
            If cal_flag = True Then
                '総合貯留関数法
                fr(t) = f * r(t)
            Else
                '貯留関数法
                fr(t) = a * f * r(t) / 3.6
            End If
            fim = fim + fr(t)
            fi(t) = fim
        Next t
    
        qsm = 0#
        For t = 2 To m - tl0
            qsm = qsm + (ql(t - 1) + ql(t)) / 2#
            qs(t) = qsm
            If fi(t) - qs(t) >= 0 Then s(t) = fi(t) - qs(t)
        Next t
    
        sm = 0
        For t = 2 To m - tl0 - 1
            If s(t) >= sm Then sm = s(t)
        Next t
    
        cal_kp  '最小自乗法によるK・Pの計算
    
        data_write2 ws1      '貯留関数式書き込み
        data_write3          '貯留高と流出高書き込み
        graph_refresh        'グラフ修正
    End If
    
    '流量計算
    
    data_read3 ws1
        
    cal_qt ws2          '求めた係数で流量を計算
        
    graph_refresh2       'グラフ修正
    
End Sub

Private Sub data_read1(ws As Worksheet)
'データ読み込み1(流域諸元)
    ws.Select
    a = ws.Cells(3, "D")        '流域面積(km2)
    qi = ws.Cells(4, "D")       '基底流量(m3/s)
    l = ws.Cells(5, "D")        '流路長(km)
    t_step = ws.Cells(6, "D")   '遅滞時間計算ピッチ(hr)
    If t_step = 0 Then t_step = 1
End Sub

Private Sub data_read2(ws As Worksheet)
'データ読み込み1(雨量・流量データ)
    Dim row_num As Integer
    Dim i As Integer
    
    '雨量・流量読み込み
    row_num = 3
    i = 0
    Do Until ws.Cells(row_num, "C") = ""    'データがあるまでループ
        i = i + 1
        r(i) = ws.Cells(row_num, "C")
        q(i) = ws.Cells(row_num, "D")
        row_num = row_num + 1
    Loop
    m = i       '雨量・流量データ数定義
    '計算対象範囲指定
    For i = 1 To m
        If ws.Cells(i + 2, "E") = 1 Then m1 = i         '開始点
        If ws.Cells(i + 2, "E") = 2 Then m2 = i         '終了点
    Next i
    '期間の累積雨量
    rt = 0
    For i = m1 To m2
        rt = rt + r(i)
    Next i
End Sub

Private Sub cal_fr()
'frの決定
    Dim t As Integer
    fim = 0
    For t = m1 To m2
        If cal_flag = True Then
            '総合貯留関数法
            fr(t) = f * r(t)
        Else
            '貯留関数法
            fr(t) = a * f * r(t) / 3.6
        End If
        '
        fim = fim + fr(t)
        fi(t) = fim
    Next t
End Sub

Private Sub cal_ql()
'遅滞時間分戻した流出高の定義(基底流量含まず)
    Dim t As Integer
    qlm = 0
    For t = 1 To m2
        If t - tl1 >= 0 Then
            ql(t - tl1) = q(t) + (q(t + 1) - q(t)) * tl2 - qi  '遅滞時間を考慮した流量(基底流量を差し引いたもの)
            
            '総合貯留関数法
            If cal_flag = True Then
                ql(t - tl1) = ql(t - tl1) / a * 3.6 '流量を流出高へ変換
            End If
            
            If ql(t - tl1) >= qlm Then
                qlm = ql(t - tl1)    '最大流量
                tpl = t - tlk        '最大流量の発生時間
            End If
        End If
    Next t
End Sub

Private Sub cal_qtm()
'遅滞時間を考慮した流量の積算(計算区間内)
    Dim t As Integer
    qtm = 0
    For t = m1 To m2 - tl0
        qtm = qtm + ql(t)
    Next t
End Sub

Private Sub cal_sa(tl As Integer)
'貯留高さと流出高をプロットした場合のループ形状
    Dim t As Integer
    qsm = 0: sar = 0: sad = 0
    
    For t = m1 To m2 - tl0 - 1
        qsm = qsm + (ql(t - 1) + ql(t)) / 2#
        qs(t) = qsm
        If fi(t) - qs(t) < 0 Then
            s(t) = 0#
        Else
            s(t) = fi(t) - qs(t)
        End If
    Next t
    '上昇ループの積算値計算
    For t = m1 To tpl - 1
        sar = sar + (s(t) + s(t + 1)) * Abs(ql(t + 1) - ql(t)) / 2#
    Next t
    '下降ループの積算値計算
    For t = tpl To m2 - tl0 - 1
        sad = sad - (s(t) + s(t + 1)) * Abs(ql(t) - ql(t + 1)) / 2#
    Next t
    'ループ全体の積算値計算
    sa(tl) = sar + sad
End Sub

Private Sub cal_kp()
'最小自乗法による貯留関数係数の計算
    Dim m3 As Integer
    Dim t As Integer
    Dim xi As Double
    Dim xx As Double
    Dim yi As Double
    Dim xy As Double
    Dim p0 As Double
    Dim k0 As Double
    Dim k1 As Double
    
    m3 = m2 - tl0
    e = 0
    For t = m1 To m3
        If s(t) <> 0 Then e = e + 1
    Next t
    xi = 0: xx = 0: yi = 0: xy = 0
    For t = m1 To m3
        xi = xi + Log(ql(t))
        xx = xx + Log(ql(t)) * Log(ql(t))
        yi = yi + Log(s(t))
        xy = xy + Log(ql(t)) * Log(s(t))
    Next t
    p0 = (xi * yi - e * xy) / (xi * xi - e * xx)
    k0 = (yi - p0 * xi) / e
    k1 = Exp(k0)
    k = k1
    p = p0
    k = Int(0.5 + k1 * 100#) / 100#     '桁合わせ
    p = Int(0.5 + p0 * 1000#) / 1000#   '桁合わせ
End Sub

Private Sub data_write1(ws As Worksheet, t As Integer)
'結果書き込み1(遅滞時間計算経過)
    Dim sa_state As String
    If t = 0 Then
        ws.Select
        Cells.Select
        Selection.ClearContents 'セルの値をクリア
        ws.Cells(2, "A") = "遅滞時間(hr)"
        ws.Cells(2, "B") = "最大流量(m3/s)"
        ws.Cells(2, "C") = "最大流量時間(hr)"
        ws.Cells(2, "D") = "流量積算値(m3/s)"
        ws.Cells(2, "E") = "流出係数"
        ws.Cells(2, "F") = "上昇ループの積算値"
        ws.Cells(2, "G") = "下降ループの積算値"
        ws.Cells(2, "H") = "ループ全体の積算値"
        ws.Cells(2, "I") = "ループの状態"
    End If
    '
    sa_state = ""
    If sa(t) > 0 Then sa_state = "時計回り"
    If sa(t) < 0 Then sa_state = "反時計回り"
    If sa(t) = 0 Then sa_state = "ぴったりOK"
    '
    ws.Cells(t + 3, "A") = tlk
    ws.Cells(t + 3, "B") = qlm
    ws.Cells(t + 3, "C") = tpl
    ws.Cells(t + 3, "D") = qtm
    ws.Cells(t + 3, "E") = f
    ws.Cells(t + 3, "F") = sar
    ws.Cells(t + 3, "G") = sad
    ws.Cells(t + 3, "H") = sa(t)
    ws.Cells(t + 3, "I") = sa_state
End Sub

Private Sub data_write2(ws As Worksheet)
'結果書き込み2(貯留関数式出力)
    ws.Select
    ws.Cells(10, "B") = k
    ws.Cells(11, "B") = p
    ws.Cells(12, "B") = tld
    ws.Cells(13, "B") = f
End Sub

Private Sub data_write3()
'貯留高と流出高書き込み
    Dim t As Integer
    Sheet8.Select
    If cal_flag = True Then
        '総合貯留関数法
        Sheet8.NAME = "貯留高と流出高の関係"
        Sheet8.Cells(1, "A") = "流出高と貯留高の関係"
        Sheet8.Cells(2, "A") = "流出高(mm/hr)"
        Sheet8.Cells(2, "B") = "貯留高(mm)"
    Else
        '貯留関数法
        Sheet8.NAME = "貯留量と流出量の関係"
        Sheet8.Cells(1, "A") = "流出量と貯留量の関係"
        Sheet8.Cells(2, "A") = "流出高(m3/s)"
        Sheet8.Cells(2, "B") = "貯留高(m3)"
    End If
    Sheet8.Cells(3, "D") = tld
    For t = 1 To m - tl0
        Sheet8.Cells(t + 2, "A") = ql(t)
        Sheet8.Cells(t + 2, "B") = s(t)
    Next t
End Sub

Private Sub data_read3(ws As Worksheet)
'データ読み込み3(貯留関数式元)
    ws.Select
    k = ws.Cells(10, "B")
    p = ws.Cells(11, "B")
    tld = ws.Cells(12, "B")
    f = ws.Cells(13, "B")
End Sub

Private Sub cal_qt(ws As Worksheet)
'求めた貯留関数を用いて流量を再度計算
    Dim t As Integer
    Dim tm As Integer
    Dim mr As Integer   '雨量・流量数
    Dim qql(0 To mnum) As Double
    Dim qqt(0 To mnum) As Double
    Dim qi2 As Double
    Dim j As Integer
    Dim c As Double
    Dim d As Double
    Dim x As Double
    Dim x2 As Double
    Dim y As Double
    
    mr = m
    tl1 = Int(tld)          '遅滞時間整数部
    tl2 = tld - tl1         '遅滞時間少数部
    
    If cal_flag = True Then
        '総合貯留関数法
        qi2 = qi / a * 3.6
    Else
        '貯留関数法
        qi2 = qi
    End If
    
     'ニュートン法による流量計算
    qql(1) = 0
    For t = 1 To mr
        If cal_flag = True Then
            '総合貯留関数法
            c = 2# * f * r(t) - qql(t) + 2# * (k * qql(t) ^ p)
        Else
            '貯留関数法
            c = 2# * f * a * r(t) / 3.6 - qql(t) + 2# * k * qql(t) ^ p
        End If
        '
        x = 0.0001
        j = j + 1
Step1:  y = x + (2# * k * x ^ p) - c
        d = 1# + 2# * p * k * x ^ (p - 1#)
        x2 = x - y / d
        If Abs(x2 - x) <= 0.0001 Then GoTo Step2
        x = x2
        GoTo Step1
Step2:  qql(t + 1) = x
        If t <= tl1 Then qqt(t) = qi2
        qqt(t + tl1) = qi2 + qql(t + 1) + (qql(t + 1) - qql(t)) * tl2
        If j <= 100 Then GoTo Step3
        If qqt(t) < 1.1 * qqt(99) Then GoTo Skip1
Step3:
    Next t
Skip1:
    For t = 1 To mr
        '総合貯留関数法
        If cal_flag = True Then
            qqt(t) = qqt(t) * a / 3.6
        End If
        ws.Cells(t + 2, "G") = qqt(t)
    Next t
End Sub

Private Sub graph_refresh()
'グラフの修正
    Dim g_range As String
    グラフ1.Select
    If cal_flag = True Then
        '総合貯留関数法
        グラフ1.NAME = "貯留高と流出高の関係図"
        ActiveChart.ChartTitle.Select
        Selection.Characters.Text = "貯留高と流出高の関係図"
        ActiveChart.Axes(xlCategory).AxisTitle.Select
        Selection.Characters.Text = "流出高(mm/hr)"
        ActiveChart.Axes(xlValue).AxisTitle.Select
        Selection.Characters.Text = "貯留高(mm)"
    Else
        '貯留関数法
        グラフ1.NAME = "貯留量と流出量の関係図"
        ActiveChart.ChartTitle.Select
        Selection.Characters.Text = "貯留量と流出量の関係図"
        ActiveChart.Axes(xlCategory).AxisTitle.Select
        Selection.Characters.Text = "流出量(m3/s)"
        ActiveChart.Axes(xlValue).AxisTitle.Select
        Selection.Characters.Text = "貯留量(m3)"
    End If
    g_range = "A3:B" & CStr(m - tl0 + 2)
'
    ActiveChart.PlotArea.Select
    ActiveChart.SetSourceData Source:=Sheet8.Range(g_range), PlotBy:=xlColumns

End Sub

Private Sub graph_refresh2()
'グラフの修正
    Dim g_range As String
    Sheets("雨量・流量図").Select
    g_range = "B2:D" & CStr(m + 2) & ",F2:G" & CStr(m + 2)
'
    ActiveChart.PlotArea.Select
    ActiveChart.SetSourceData Source:=Sheets("雨量・流量データ").Range(g_range), PlotBy:=xlColumns
End Sub

Private Sub data_write0(t As Integer)
'結果書き込み(デバック)
    Dim ws As Worksheet
    Dim i As Integer
    Set ws = Sheets("デバッグQl")
    ws.Cells(2, t + 2) = t
    ws.Cells(3, t + 2) = tld
    ws.Cells(4, t + 2) = tl1
    ws.Cells(5, t + 2) = tl2
    For i = 1 To m2 - tl0
        'If i - tl1 >= 0 Then
        '    ws.Cells(i + 5, t + 2) = ql(i - tl1)
        'End If
        ws.Cells(i + 5, t + 2) = ql(i)
    Next i
End Sub

tys-sfb_v1_0_0

tys-sfb_v1_0_0.lzhをダウンロードして解凍すると、「tys-sfb_v1_0_0.xls」ができます。 この「tys-sfb_v1_0_0.xls」の使用方法を説明します。 ダウンロードについては、水理水文ソフトウェアのHPを参照してください。

マクロの説明

本マクロは、各条件シートに条件を入力するだけで、複数流域の総合貯留関数による流出解析計算が行えます。

  1. 複数流域を一括して流出解析を行う。

計算条件とデータの入力

計算条件

計算条件」シートで貯留関数係数の回帰を行うかどうかをチェックしてください。

計算条件シート

計算条件シート

 
流域流出計算/貯留関数法.1344491215.txt.gz · 最終更新: 2012/08/09 14:46 by tys
[unknown button type]
 
特に明示されていない限り、本Wikiの内容は次のライセンスに従います: CC Attribution-Share Alike 3.0 Unported
Recent changes RSS feed Donate Powered by PHP Valid XHTML 1.0 Valid CSS Driven by DokuWiki