【流出解析】貯留関数法のプログラムを作ってみる

はじめに
基礎方程式の導出
数値解法
ニュートン法
入力データ
プログラムと計算結果
引用文献

はじめに

もともと流出解析はやっていたので，貯留関数法は知っていた（聞いたことがある程度だ）が，実際にその計算をしたことはなかった．そこで，少し時間と気持ちに余裕があるので，貯留関数法のプログラムを作ることにした．

貯留関数法とは，流域および河道を一つの貯水池と考えて，そこでの貯留量と流出量の間に一義的な関係性を仮定し，貯留量を媒介変数として，降雨量（流入量）から流出量を求める手法のことである．

基礎方程式の導出

貯留関数法では，河川流域に降った雨が流域に貯留された後に流出する過程を，流域を下図のような貯留タンクに見立てて考える．

f:id:okome29taberu:20200426202628p:plain:w400

流域（タンク）に降った降雨 $r$ mm/sはタンクに一時的に貯留され，貯留量 $S$ mmに応じた流量 $q$ mm/sがタンクから流出する．この時の降雨というのは，流域に降った雨のうち，樹木による遮断や地中による貯留の効果によって流域からの流出とならない損失分を除いたもの（有効雨量）である．ここで，貯留量 $S$ と流出量 $q$ の間には以下の式の関係があると仮定する．

$S=kq^p$ 　　（運動方程式）

ここで，k，pはパラメータである．次に，Δt時間におけるタンクの貯留量変化ΔSは，

$ΔS=rΔt-qΔt$

と表せられるので，以下の式が得られる．

$\frac{dS}{dt}=r-q$ 　　（連続式）

これら貯留関数法の運動方程式と連続式を連立して解くと，流域からの流出量を得られる．

数値解法

はじめに，貯留関数法の連続式を前進差分法で近似すると，

$\frac{S_{i+1}-S_{i}}{Δt}=r_{i+1}-\frac{q_{i+1}+q_{i}}{2}$

となり， $S=kq^{p}$ であるから，

$\frac{kq_{i+1}^{p}-kq_{i}^{p}}{Δt}=r_{i+1}-\frac{q_{i+1}+q_{i}}{2}$

となる．さらに， $q_{i+1}$ について整理すると，

$\frac{2k}{Δt }q_{i+1}^{p}+q_{i+1}-2(r_{i+1}-\frac{q_{i}}{2}+\frac{k}{Δt}q_{i}^{p})= 0$

となる． $q_{i+1}$ を求めるためには，この式を数値的に解けばよいので，本稿ではニュートン法で求める．

ニュートン法

ニュートン法は，主として多項式の実根を求めるために使用される．例えば， $x$ に関する方程式 $f(x)=0$ の根の近似値を $x_{1}$ とし，これを最初の試行値とすると，次の試行値 $x_{2}$ は $x_{1}$ における関数 $f(x_{1})=0$ の接線を求めればよい．そして，接線とx軸の交点が次の試行値 $x_{2}$ となる．接線の勾配を $m$ とすると，

$m=\frac{f(x_{1})}{x_{1}-x_{2}}=f'(x_{1})$

となるので，

$x_{2}=x_{1}-\frac{f(x_{1} )}{m}= x_{1}-\frac{f(x_{1})}{f'(x_{1})}$

以下同様に，

$x_{n+1}= x_{n}-\frac{f(x_{n})}{f'(x_{n})}$

として（ $x_{1}，x_{2}，…，x_{n+1}$ ）を求めていけば，最終的に真値に収束する．この式を用いて $q_{i+1}$ を求めるとすると， $q_{i+1}=x$ ， $r_{i+1}-\frac{q_{i}}{2}+\frac{k}{Δt}q_{i}^{p}=C$ とおいて，

$f(x)=\frac{2k}{Δt}x^{p}+x-2C$

$f'(x)=\frac{2kp}{Δt}x^{p-1}+1$

をそれぞれ使えばよい．

入力データ

貯留関数法による流出量の計算に必要なデータは，降雨量 $r$ およびパラメータ $k$ ， $p$ である．これらの値は解析対象とする洪水によって降雨量が異なるのはもちろんのこと，パラメータも流域に固有であるため各流域で異なる．本稿では，降雨量データの入手やパラメータの設定の手間を省くため，財）北海道河川防災研究センター・研究所（2006）に記載されている湧別川・丸瀬布地点における2001年9月10日～15日に発生した洪水のデータを利用する．そこでは，パラメータは $k=22.7602$ ， $p=0.6$ と設定されている．

プログラムと計算結果

本稿では，以下のようなプログラムをC++で書いた．ただC++は勉強し始めたばかりであるため，適切な書き方でないかもしれないが，ゆるしてほしい．

#include <iostream>
#include <iomanip>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

using namespace std;

// Parameter setting
double k=22.7602;
double p=0.6;
double dt=1.;       // time step [h]
double lim=1.0e-6;  // 収束判断値

// f(x)
double f (double x, double r, double xold){
  return 2.*k/dt*pow(x,p) + x -2.*(r-xold/2.+k/dt*pow(xold,p));
}

// f'(x)
double df (double x){
  return 2.*k*p/dt*pow(x,p-1.)+1.;
}

// Newton Raphson method
double NewtonRaphson (double x, double r){

  int i=0;
  int imax=10000;   // 最大反復回数
  double xn;        // 近似値＝次の試行値
  double xold=x;    // 時間tにおける流出量[mm/h]
  
  // 反復開始
  while(i<=imax){
    i++;
    xn=x-f(x,r,xold)/df(x); 
    if(fabs(xn-x)<lim) break;   // 値が収束していた場合は反復から抜ける
    x=xn;                       // 近似値を次の試行値に置き換える
  }
  x=xn;
  return x;
}

int main (){

    char fileName[256];       // File name
    string ss;                // Dammy line
  double q;             // 時間tにおける流出量[mm/h]
  double qn;            // 時間t+1における流出量[mm/h]
  double qobs;          // 実測流出量[mm/h]
  double r;             // 時間tからt+1における有降雨量[mm/h]
  double tmax=94;       // 計算時間[h]

  // 入力データファイルの設定
  ifstream ifs("./data.txt");
  if(!ifs){
    cout<<"cannot open file"<<endl;
    exit(1);
  }
  getline(ifs,ss);

  // 出力データファイルの設定
  ofstream ofs("./out.txt");
  ofs<<"time[h]\train[mm/h]\tQsim[mm/h]\tQobs[mm/h]"<<endl;


  //流出量の初期値
  q=0.0001;
  
  // 計算開始
  for(int t=1;t<tmax; t++){
    ifs>>ss>>r>>qobs;       // 時間tからt+1における有降雨量の読み込み
    qn=NewtonRaphson(q,r);
    q=qn;                   // 流出量の値の更新

    // 計算結果の出力
    ofs<<std::fixed<<std::setprecision(4)<<t<<"\t"<<r<<"\t"<<q<<"\t"<<qobs<<endl;
    cout<<t<<endl;
  }

  return 0;
}

上記のプログラムを用いた計算結果は以下の通りである．定量的な評価を本稿では実施していないので，強くは言えないが，計算値は，大まかではあるものの実測値を再現できているように見える．この計算結果以上に再現性を上昇させるためには，パラメータを最適化するなどの手続きが必要であろう．

f:id:okome29taberu:20200426211816p:plain:w500

引用文献

（財）北海道河川防災研究センター・研究所（2006）実践流出解析ゼミQ&A編

雑記

日々の記録になります．