水文統計

水文量の統計的特性

水文諸量は自然現象である降水に起因しているところから、その発生規模は確率的には小さくても際限なく大きいものと小さいものが起る可能性が常に存在します。一方、水工計画や水理構造物の設計において議論の対象となるのは、数十年に1回というような非常に低い確率で生起するまれな水文事象です。通常の統計解析では平均値近傍の分布特性が最大の関心事であるのに反して、水文統計学ではその両端部の議論が重視されることが大きな特徴です。

したがって、水文統計学は水文資料から有用な確率統計的特性値を引き出し、その情報を水工計画・管理に活用しようとする学問です。当然のことながら、水文統計解析は、基礎資料となる水文観測資料の量的・質的制約を受けます。量的制約は観測期間の長さに依存し、質的制約は資料の精度に関連します。水文観測は自然現象を相手にするため観測のやり直しができず、欠測値の補充ができないことも水文統計解析上の問題点です。また、流域の都市化に伴う流出機構の変化など資料の統計的均質性が疑問視される事柄もありますが、その詳細を短い観測値のうえで実証することは必ずしも簡単ではありません。いずれにせよ、水文統計解析は上述した種々の制約条件のもとで水工計画上の意思決定に必要な情報を供給せざるをえない実状にあります。

水文量の変動特性は大きく分布特性と時系列特性に分類されます。分布特性値として重要な統計量は平均値、標準偏差(あるいは変動係数)およびひずみ係数です。時系列特性とは、水文量が時系列的にみて独立に生起しているかあるいは持続性を示しているかの性質をいいます。この特性値として自己相関係数が重要です。

したがって、これらの分布・時系列の特性値が水工計画上どのような影響を及ぼすかを念頭に入れておいた方がよいです。例えば、治水計画における計画高水流量の決定に際して、確率分布形の非対称性を表現するひずみ係数が無視されるとき(正規分布を仮定)、設計値は過小評価され、危険側の計画となります。一方、利水計画においては貯水池容量の決定問題があります。このとき、ひずみ係数が大きくなるにつれて容量は減少する反面、自己相関係数が大きくなるにつれて容量は増加します。

確率統計水文学の理論体系は1変数分布理論と多変数分布理論に大別されます。前理論の適用にあたっての大きな障害は分布の非正規性と水文量の時間的・空間的相関性です。治水計画ではその性質上、年最大水文量が解析対象となりますが、標本は独立性と等質性を満足しているものと仮定する場合が多いです。近年、治水・水資源システムがますます複雑化し、多変数分布理論を用いて処理しなければならない問題も多くなってきています。

水文量の分布とその解析

確率年と確率水文量

確率変数$X$の確率密度関数を$f(x)$、その分布関数を$F(x)$とすれば、次式が成立します。

$$F(x) = \int_{-\infty}^{x} f(x) dx \tag{4.1}$$

式(4.1)が確率統計解析における基本確率モデルです。

$X$がある特定値$x_T$以下となる確率が$F(x_T)$であり、この$F(x_T)$を非超過確率、$W(x_T) = 1 - F(x_T)$を超過確率といいます。

治水計画においては変数の大きいところが問題となるので、$W(x_T)$が治水安全度の指標とされます。一方、渇水のように小さいときが問題となる場合には$F(x_T)$がその指標とされます。

$X$が$x_T$以上となるようなことが平均的にみて$T$年に1度の割合で生起することが期待されるとき、この$T$をリターンピリオド(再現期間)と呼んでいます。この場合、リターンピリオド$T$と超過確率$W(x_T)$の間には次の関係があります。

$$T = \frac{1}{m W(x_T)} = \frac{1}{m[1 - F(x_T)]} \tag{4.2}$$

ここに、$m$は対象とする水文量の年平均生起回数であり、年最大水文量を取り扱う場合、$m=1$となります。このとき、リターンピリオドは超過確率の逆数となるので、$T$を確率年ともいい、$x_T$を$T$年確率水文量と呼びます。$X$が$x_T$以下となるリターンピリオドも同様に定義されます。

リターンピリオド$T$年の意味は超過確率$1/T$レベルの水文量の生起時間間隔に関する期待値であって、治水計画そのものの危険度が$1/T$でないことに留意する必要があります。

$T$年確率事象が特定の期間$N$年間に少なくとも1回起こる確率$P_r$は、ベルヌイ試行として、次式で与えられます。

$$P_r = 1 - \left(1 - \frac{1}{T}\right)^N \tag{4.3}$$

例えば、年最大日雨量資料が独立でかつ等質な時系列と仮定します。このとき、100年確率雨量以上の雨が10年間に少なくとも1回起る確率は約10%であり、決して小さい確率とはいえません。

すなわち、異常降雨がいつどこで起っても不思議ではありません。また、$T$年確率事象が少なくとも1回起るか、あるいは全く起らない半々の確率($P_r = 0.5$)となる期間は式(4.3)より$N \fallingdotseq 0.7T$となります。

式(4.1)の確率密度関数$f(x)$が既知であれば、特定値$x_T$に対する確率年$T$、あるいは特定の確率年$T$に対する水文量$x_T$は理論的には求まります。

しかしながら、実際の適用にあたってはいくつかの困難な問題が生じます。

計算の簡便性

実用的観点からは、式(4.1)の分布関数が変数$x$の陽形式で表現できる確率密度関数$f(x)$が選定できれば解析に有利です。その分布形としては、指数分布、一般化極値分布(グンベル分布を含む)、平方根指数型最大値分布等があげられます。

一方、水文統計解析で常用されている正規分布、対数正規分布、ガンマ型分布等については、確率密度関数を解析的には積分はできません。このため、数表や数値積分、ないしは近似式の利用によって確率水文量等を算定せざるをえません。なお、Hastingsによって提案された正規分布に関する近似式は精度が高く、また、ガンマ型分布に関する近似式として、後述するWilson-Hilferty変換があります。

確率分布モデルの選択

従来から数多くの確率分布形が水文量の頻度分布解析に用いられてきました。しかしながら、ある水文量が特定の確率分布モデルに従うという物理的根拠は何もありません。

現在主用されている確率密度関数は母数の数によって、2母数モデルと3母数モデルに大別されます。母数の数が多くなるほど観測資料への適合度が高くなる場合が多い反面、積率法による母数推定にみられるようにその推定誤差が大きくなる問題もあります。従来、20〜30年間の観測資料からは統計的適合度検定だけでは確率分布形選択の優劣をつけがたいという経緯もあり、分布形状の多様性がその結果として生じたとする意見もあります。

わが国においても、気象・水文観測体制が広域でかなり整備されてきており、観測資料も60〜100年にわたってデータベース化が進んでいます。現在、この資料集積と相まって、分布形選択に関する系統的調査・研究の機運が高まっています。

母数推定法

確率密度関数に含まれる母数を観測資料から推定する方法も、現在非常に多岐にわたっています。

従来、計算の簡便さから確率紙を用いた図式推定法や積率法が多用されてきましたが、近年の電子計算機の発達により複雑な最尤法による解の導出も可能となってきています。以上のほかに、クオンタイル法、セクスタイル法、最大エントロピー法、後述するPWM法があげられます。

推定量の重要な性質として不偏性(母数推定値の期待値が母数に一致する)と有効性(不偏推定値が最小分散を有する)があります。この特性は標本数、分布特性(変動特性とひずみ特性)および確率密度関数によっても異なるため、最適な母数推定法を選定することは、一般的には困難です。積率法においては、不偏推定値を得るための小標本によるひずみ係数の補正が必要です。また、水文資料の中に極端に小さい値や大きい値が含まれているときには、標本ひずみ係数は見かけ上大きくなる場合があります。したがって、それらの値を異常値とみなしてよいかどうかを検定する方法も提案されています。

適合度検討

あてはめた分布の適合度を視覚的にとらえる手段として、従来から確率紙を利用する方法、観測値のヒストグラムと理論曲線とを比較する方法が広範に用いられています。

水文解析では通常極端に大きい値と小さい値の確率評価が問題となります。したがって、分布モデルと母数推定の優劣の比較においては資料全体に対する全般的適合度のほかに、分布の両端部における形状や適合度が最重視されます。以上の観点に立って、実用上も容易に行える適合度の数値化による評価法が提案されています。

確率水文量の推定

治水計画においては超過確率が小さい部分での$T$年確率水文量の推定が必要となりますが、この問題は基本的に外挿補間法と同等です。

したがって、確率水文量は用いた標本数、分布形および推定法によっても異なった値をとります。

特に、水文資料の蓄積が進んでいる現在、標本数の違いによって確率水文量が大きく変動しないような確率分布モデルと母数推定法を選定することが望ましいです。

主要確率分布モデル

水文統計解析上有用と考えられる代表的な確率分布モデルと実際に適用する際の留意点は次のとおりです。

対数正規分布

対数正規分布が古くから水文統計学の分野に導入された最大の理由は、水文量の対数変換標本値が近似的に正規分布に従う場合が多いことが経験的に知られていたことによります。

現在多用されている対数正規分布の表現形式は次の2種類に大別されます。

1. $x$を対数正規変量としたとき、関係式$x = a + \exp(y)$中の変量$y$は平均値$\mu_y$、分散$\sigma_y^2$の正規分布$N(\mu_y, \sigma_y^2)$に従う。ここに、$a$は分布下限値。この形式は欧米でよく用いられています。
2. 関係式$x = a + c\exp(kz)$における変量$z$は標準正規分布$N(0, 1)$に従う。ここに、$a, c, k$は分布母数。この形式はわが国で主に用いられているいわゆるSlade型分布に相当します。

上記2型式の母数の間には$c = \exp(\mu_y), k = \sigma_y$の互換性があるので、どちらの形式を用いても実際上問題とはなりません。

ただし、正規分布の近似度を測る指標としては、対数変換標本値$y$のひずみ係数が十分な精度で零に近いかどうかを検証すればよいわけで、その意味では2の形式の方が実用上便利です。

対数正規変量$x$の平均値$\mu_x$、分散$\sigma_x^2$、ひずみ係数$\gamma_x$と分布母数の間には次の関係式が成り立ちます。

$\mu_x = a + c\phi^{1/2} \qquad (4.5)$ $\sigma_x^2 = c^2[\phi(\phi - 1)] \qquad (4.6)$ $\gamma_x = (\phi - 1)^{1/2} (\phi + 2) \qquad (4.7)$

ここに、$\phi = \exp(k^2)$。積率法による母数推定には、$\mu_x, \sigma_x^2, \gamma_x$にそれぞれ標本積率値を代入して、式(4.5)〜(4.7)を$a, c, k$について解けばよいです。

式(4.7)は$\phi$に関して1実根を有し、その根は次式で得られます。

$\phi = [s + (s^2 - 1)^{1/2}]^{1/3} + [s - (s^2 - 1)^{1/2}]^{1/3} - 1 \qquad (4.8)$

ここに、$s = 1 + (\gamma_x / 2)$。すなわち、$\gamma_x$が既知であれば、$\phi$は陽形式で求められ、母数$k$は$k = (\ln \phi)^{1/2}$で計算されます。他の母数$a$と$c$は式(4.5)、式(4.6)より容易に算定可能です。

対数正規分布の他の母数推定法として、岩井法で代表されるクオンタイル法があります。この推定法は、特に分布両端部での適合度を高めることに主眼をおいた有用な方法であり、積率法と同様に簡便であるため、従来からわが国でも多用されてきています。

クオンタイル法では、まず分布下限値$a$を次式で推定します。

$a = \frac{x_{(1)} - x_{(N)} z_m^{1/\{1 + [(z_m - z_1)/(z_N - z_m)]^2\}}}{1 - z_m^{1/\{1 + [(z_m - z_1)/(z_N - z_m)]^2\}}} \qquad (4.9)$

ここに、$x_{(1)}, x_{(N)}, x_m$はそれぞれ標本最小値、最大値、メディアンです。他の母数$c, k$は最尤法を併用して、式(4.10)(4.11)で計算されます。

$\ln c = \frac{1}{N} \sum_{i=1}^N \ln (x_i - a) \qquad (4.10)$ $k^2 = \frac{1}{N} \sum_{i=1}^N [\ln (x_i - a) - \ln c]^2 \qquad (4.11)$

ここに、$x_i$は水文観測値、$N$は標本数です。

ガンマ型分布

水文統計で使用されるガンマ分布族としては、ピアソン3型分布、対数ピアソン3型分布等があげられます。ここでは、ピアソン3型分布を代表例としてとりあげます。

$y$をピアソン3型分布に従う変量とすると、その密度関数は次式で与えられます。

$f(y) = \frac{1}{\Gamma(b)} \left(\frac{y - c}{a}\right)^{b-1} \exp\left(-\frac{y - c}{a}\right) \qquad (4.12)$

ここに、$a > 0$のとき$c \leq y < \infty$、$a < 0$のとき$-\infty < y \leq c$、$a$は尺度母数、$b$は形状母数($b > 0$)、$c$は位置母数、$\Gamma(b)$はガンマ関数です。表-4.2の分布形状は$a > 0$として、母数$b$の選択により、左側から$0 < b < 1$、$b = 1$(指数分布)、$b > 1$、$b >> 1$としたガンマ型分布で統一的に表現できるのがこの分布が多用される理由でもあります。

式(4.12)の母統計量は次式で与えられます。

$\mu_y = c + ab, \quad \sigma_y^2 = ab^2, \quad \gamma_y = \frac{2}{\sqrt{b}} \qquad (4.13)$

ここに、$\mu_y$は平均値、$\sigma_y^2$は分散、$\gamma_y$はひずみ係数です。

$x$を原水文量として、対数変換値$y = \ln x$の変量$y$が式(4.12)の密度関数で表わされるとき、$x$を対数ピアソン3型分布に従う変量といいます。したがって、対数ピアソン3型分布の形状は母数$a$に依存します。すなわち、$\gamma_y > 0$のとき$a > 0$、$\exp(c) \leq x < \infty$、$\gamma_y < 0$のとき$a < 0$、$0 \leq x \leq \exp(c)$。

水文量の対数変換標本値は、負のひずみ係数をもつ場合が多いです。したがって、式(4.13)による対数ピアソン3型分布の母数推定法(積率法)では、あてはめた分布形が上限値をもつ頻度が大きくなります。このため、この上限値をこえる水文量が起りえないという物理的不都合が生ずるので、適用にあたっては注意を要します。

ガンマ型分布の実際問題への適用にあたっては、標準ガンマ分布が基本となります。標準ガンマ分布は式(4.12)で$w = (y - c)/a$の変数変換を行った1母数ガンマ分布です。すなわち、

$f(w) = \frac{1}{\Gamma(b)} w^{b-1} \exp(-w), \quad 0 \leq w < \infty \qquad (4.14)$

したがって、標準ガンマ分布の母統計量は式(4.13)で$c = 0$、$a = 1$とすればよいです。また、標準ガンマ分布の分布関数$F(w)$は以下のように表わせます。

$p = F(w) = \begin{cases} \frac{1}{\Gamma(b)} \int_0^w t^{b-1} \exp(-t) dt, & a > 0 \\ 1 - \frac{1}{\Gamma(b)} \int_{-\infty}^w t^{b-1} \exp(-t) dt, & a < 0 \end{cases} \qquad (4.15)$

ガンマ型分布を用いた確率水文量推定問題では、式(4.15)において非超過確率$p$に対する標準ガンマ変量$w$を求める必要があります。もし、$p$に対する$w_p$が算定できれば、ピアソン3型変量$y_p$は$y_p = c + aw_p$で、また対数ピアソン3型変量$x_p$は$x_p = \exp(c + aw_p)$で容易に求まります。

ガンマ型分布は母数$b$の選択によってかなりの範囲の形状に対処できる反面、式(4.15)に不完全ガンマ関数比の解法を含むのが実際上の適用にあたっての難点です。このため式(4.15)の近似解として導出されたのがWilson-Hilferty変換です。この変換は次の理論的根拠に基づいています。

すなわち、$(w/b)^{1/3}$は平均値$1 - (1/9b)$、分散$1/9b$をもつ正規分布で近似できます。したがって、次式が成立します。

$w = b\left(1 - \frac{1}{9b} + \frac{t}{9\sqrt{b}}\right)^3 \qquad (4.16)$

ここに、$w$は式(4.14)に従う標準ガンマ変量、$t$は$N(0, 1)$に従う標準正規変量です。

なお、形状母数$b$が0.44より小さくなる(ひずみ係数$\gamma_y > 3.0$)と、式(4.16)の精度が落ちるので、このときは修正Wilson-Hilferty変換を用いた方がよいです。

一般化極値分布

グンベル分布として知られている最大値分布の第1形式は、年最大洪水量の分布関数としての有用性が指摘されて以来、水文統計解析で数多く適用されてきています。Jenkinsonは3種の極値極限分布を1つの式形に統一して、一般化極値分布(Generalized Extreme Value分布：GEV分布)の導入を図りました。

GEV分布の分布関数$F(x)$は次式で定義されます。

$F(x) = \begin{cases} \exp\left[-\left\{1 - k\frac{x - c}{a}\right\}^{1/k}\right], & k \neq 0 \\ \exp\left[-\exp\left\{-\frac{x - c}{a}\right\}\right], & k = 0 \end{cases} \qquad (4.17)$

ここに、$k > 0$のとき$-\infty < x \leq c + (a/k)$であり、$k < 0$のとき$c + (a/k) \leq x < \infty$となります。$c$と$a$は、それぞれ位置母数と尺度母数です。また、$k$は極値分布がどの形式に属するかを決める形状母数です。特に、$k = 0$のとき、GEV分布はよく知られているグンベル分布に一致します。

式(4.17)の逆分布関数$x(F)$は次式で与えられます。

$x(F) = \begin{cases} c + \frac{a}{k}(1 - (-\ln F)^k), & k \neq 0 \\ c - a\ln(-\ln F), & k = 0 \end{cases} \qquad (4.18)$

GEV分布は、対数正規分布やガンマ型分布などと異なり、分布母数が既知の場合、非超過確率$F(x)$に対応する確率水文量$x$が式(4.18)により容易に求められます。この性質は実務面での適用にあたってもきわめて重要な特徴といえます。

また、一般化極値分布は、実際上$-1/3 \leq k \leq 1/3$の母数に対して広範囲の分布形状を表現できます。

一般化極値分布の母数推定法として、Sextile法、最尤法、確率重みつき積率法(Probability Weighted Moments: PWM法)が試みられています。このうち、実用上の簡便さを重視すれば、PWM法が有用でしょう。PWM法は特に、分布関数$F_0 = F(x)$とその逆分布関数$x \equiv x(F)$が陽形式で表現できる確率分布形の母数推定に威力を発揮します。

確率重みつき積率を以下のように定義します。

$M_{p,r,s} = E[X^p\{F(X)\}^r(1 - F(X))^s] = \int_0^1 \{x(F)\}^p F^r(1 - F)^s dF \qquad (4.19)$

ここに、$E$は期待値演算子であり、$p, r, s$は実数です。$M_{p,0,0}$は確率変数$X$の原点のまわりの$p$次の積率です。水文量のように、観測資料数が少ない場合には、高次の積率は標本誤差が大きくなるので、実用的には$M_{1,r,s}$を用いることが望ましいです。$r$と$s$の決定には各々の確率分布形について式(4.19)の積分が可能なように選べばよいです。

GEV分布については$\beta_r \equiv M_{1,r,0} = E[X\{F(X)\}^r]$($r = 0, 1, 2$)の積分が可能です。すなわち、

$\beta_r = c + \frac{a}{r + 1} [1 - (r + 1)^{-k}\Gamma(1 + k)], \quad k > -1 \qquad (4.20)$

ここに、$\Gamma(\cdot)$はガンマ関数です。したがって、式(4.20)より次の関係式が得られます。

$a = \frac{\beta_0 - 2\beta_1}{1 - 2^{-k}\Gamma(1 + k)} \qquad (4.21)$ $c = \beta_0 - a[1 - \Gamma(1 + k)] \qquad (4.22)$ $\frac{2\beta_1 - \beta_0}{3\beta_2 - \beta_0} = \frac{1 - 2^{-k}}{1 - 3^{-k}} \qquad (4.23)$

式(4.23)の実用推定式として次式が提案されています。

$k = 7.8590d + 2.9554d^2 \qquad (4.24)$

ここに、

$d = \frac{2\beta_1 - \beta_0}{3\beta_2 - \beta_0} - \frac{\ln 2}{\ln 3} \qquad (4.25)$

他の母数$a$と$c$は式(4.21)、式(4.22)により容易に算出されます。PWM法では$\beta_r$($r = 0, 1, 2$)の標本推定値が必要であり、次式を用いて計算します。

$\beta_0 = \frac{1}{N} \sum_{j=1}^N X_{(j)} \qquad (4.26)$

$\beta_1 = \frac{1}{N(N - 1)} \sum_{j=1}^N (j - 1)X_{(j)} \qquad (4.27)$

$\beta_2 = \frac{1}{N(N - 1)(N - 2)} \sum_{j=1}^N (j - 1)(j - 2)X_{(j)} \qquad (4.28)$

ここに、$X_{(j)}$は$N$個の標本値を大きさの順序に並べ換えたときの小さい方より$j$番目の値。

以下は、文章内容を変えずに、マークダウン形式で整形し、数式をLaTeX形式に変換したものです。

水文量の時系列特性

水文量時系列は表-4.1に示したように、傾向成分、周期成分および偶然変動成分からなります。

傾向成分は、いわゆる水文量の経年的変化が増加傾向にあるのか減少傾向にあるのかといった長期的傾向(トレンド)を意味します。

また、季節水文量(例えば、月降水量や月流量)時系列には1年の卓越周期が認められます。

以上2つの確定変動成分の検出・同定のために移動平均法や(直交)多項式による方法が用いられます^94。

観測時系列から確定変動成分を差し引けば、偶然(確率的)変動成分が得られます。なお、周期成分の検出と除去の程度を確認するために自己相関係数やスペクトルを調べればよいです。

確率的変動成分は統計的性質が時間的に変化しない定常確率過程とそうでない非定常確率過程に分類されます。実際面では定常時系列理論の応用がほとんどであり、2次の積率(共分散)までを考える弱定常過程を扱う場合が多いです。

なお、Box-Jenkins(ボックス・ジェンキンズ)によって開発された非定常時系列理論は水文学でも多くの適用例がみられます^94, 111。この理論の特徴は、「均質な非定常性」という概念で、時系列値そのものは定常でないが、時系列の差分は定常性を示すという考えです。

線形定常過程

線形定常過程を表現する一般形は次式で与えられます。

$z_t=\sum_{i=1}^p \phi_i z_{t-i} - \sum_{j=1}^q \theta_j \eta_{t-j} + \eta_t$ (4.29)

ここに、$Z_t$は標準正規分布($N(0,1)$)に従う時系列、$\eta_t$は白色雑音過程ともいわれる$N(0,\sigma_{\eta}^2)$に従う独立正規変量、$\phi_i$、$\theta_j$は係数です。

式(4.29)のモデルを自己回帰・移動平均(autoregressive-moving average)過程とし、ARMA($p,q$)と表わすことが多いです。

さらに、ARMA($p,q$)過程には2つの特殊例があります。

1つは自己回帰(autoregressive)過程で式(4.29)の右辺第1項と第3項からなり、AR($p$)と表現します。

もう1つは移動平均(moving average)過程で、第2項と第3項からなり、MA($q$)で表わします。

特に、AR($p$)過程は線形流出解析の単位図法の構造と基本的に同等であり、AR過程の段通次数決定法および$\phi_i$の係数同定法も提案されています^47。

水文量模擬発生モデル

水工計画、特に水資源計画にあたっては十分な長さの流量観測資料が存在していることはまれであり、このことが計画の策定をより困難にしている場合が多いです。この種の困難を克服するために用いられる手法が、モンテカルロ法による水文量資料の模擬発生(data generation)です。ここでは、年流量と月流量(季節流量)発生手法を中心に述べます。

(1) AR(1)過程による年流量時系列発生

AR(1)過程は1次マルコフ過程とも呼ばれます。AR係数$\phi_1$と1次自己相関係数$\rho_1$の間には$\phi_1=\rho_1$が成立します。また、自己相関係数$\rho_k$ ($k \geq 2$)は$\rho_k=\rho_1^k$ ($k$はラグタイムを示す)で与えられます。さらに、白色雑音過程$\eta_t$の分散は$\sigma^2_{\eta}=1-\rho_1^2$となります。したがって、式(4.29)のAR(1)過程表現は次式のようにも書き換えることができます。

$Z_t = \rho_1 Z_{t-1} + \sqrt{1-\rho_1^2} \varepsilon_t$ (4.30)

ここに、$\varepsilon_t$は$N(0,1)$に従う白色雑音過程です。

なお、年流量時系列の確率分布が非対称性を示すときには、4.2.3で述べた対数正規分布を用いると解析が容易となります。すなわち

$x_t = a + c \exp(k z_t)$ (4.31)

ここに、$x_t$は原年流量時系列です。

年流量時系列の模擬発生計算手順は次のようになります。

①年流量標本値から分布母数$a,c,k$を推定する。

②対数標本値$z_t$の1次自己相関係数$\rho_1$を計算する。

③式(4.30)の$Z_0$と$\eta_1$に標準正規乱数を発生させて$Z_1$を得る。同様に$\eta_t$に標準正規乱数を逐次発生させて、式(4.30)を適用していけば任意の長さの$Z_t$時系列が得られる。

④$Z_t$系列を式(4.31)に代入すれば、年流量時系列$x_t$が模擬発生されたことになります。

(2) ARMA(l,l) 過程による年流量時系列発生

式(4.29)よりARMA(1,1) 過程は次式で定義されます。

$Z_t = \phi Z_{t-1} + \eta_t - \theta \eta_{t-1}$ (4.32)

ここに、$\phi_1 = \phi$、$\theta_1 = \theta$と略記してあります。

ARMA(1,1)過程の自己相関係数と白色雑音過程$\eta_t$の分散は次式で与えられます。

$\rho_1 = \frac{(1-\phi \theta)(\phi - \theta)}{1-2\phi \theta + \theta^2}$ (4.33)

$\rho_k = \phi^{k-1} \rho_1, \quad (k \geq 2)$ (4.34)

$\sigma_{\eta}^2 = \frac{1-\theta^2}{1-2\phi \theta + \theta^2}$ (4.35)

AR(1)過程の自己相関係数は、ラグタイム$k$が増加するのに伴い急速に減衰するのに対し、ARMA(1,1)過程のそれは定数$\phi$と$\theta$の選択によっては緩やかに減衰します。この特性は、両過程の大きな相違点です。したがって、AR(1)過程は短期持続性、ARMA(1,1)過程は長期持続性モデルを代表します。

特に、ARMA(1,1)モデルはHurst(ハースト)現象と呼ばれる年流量時系列の長期持続性を保存できる最も簡単なモデルです^128。

(3) AR(l)過程による月流量時系列発生

式(4.30)は月流量時系列模擬発生モデルに容易に拡張できます。すなわち

$z_{t,j} = \rho_{j,j-1} z_{t,j-1} + \sqrt{1-\rho_{j,j-1}^2} \varepsilon_{t,j}$ (4.36)

ここに、$z_{t,j}$は第$t$年第$j$月の$N(0,1)$に従う正規変量、$\rho_{j,j-1}$は$(j-1)$月と$j$月の間の相関係数、$\varepsilon_{t,j}$は$N(0,1)$に従う白色雑音成分です。

月流量が対数正規分布に従うとした場合、式(4.31)に準じて月流量の分布母数を推定したのち、$\rho_{j,j-1}$を計算します。月流量発生の手順も年流量の方法に準拠します。

この方法の最大の欠点は、すべての月流量間の相関係数が保存されないために、年流量の1次自己相関係数と分散が実測年流量のそれらに比べて、一般に小さくなることです。したがって、必要貯水池容量は過少評価されます。

この問題を解決するための方策として、次のモデルが提案されています。

(4) Disaggregation(分割)モデルによる月流量発生(01),129)

この方法では、まずAR(1)過程かARMA(1,1)過程で時系列特性が再現できるように年流量を発生させます。次に、すべての月流量の確率分布・相関特性が保存できるように発生年流量時系列を月流量時系列に分割します。

確率降雨

流量資料は降雨資料に比べて、観測期間の長さと資料の統計的均質性という点で劣ります。このため、洪水防御計画や市街地排水計画は、一般に降雨量を基礎として策定され、確率降雨と降雨波形の決定が中心課題となります^125。計画降雨が設定されれば、なんらかの流出モデルを媒介としてピーク流量や計画ハイドログラフが決定されます。

ところで、降雨量は表4.2に示したように、解析時間単位によって確率・時系列特性がまったく異なったものとなります。特に、時間単位が小さくなるにつれて持続性が無視できなくなります。また、ハイエトグラフを規定する変量として、総雨量、継続時間およびピーク雨量があげられますが、この変量間には高い相関性がみられます。

従来の計画降雨量の策定では、変量間の相関性を無視した独立事象を仮定し、単変数理論を基礎とした経験公式が主流を占めてきました。しかしながら、最近計画基本量としての降雨量の機率統計的変動特性がますます重要視され、理論的検討も急速に進歩してきています。たとえば、降雨量の時系列解析手法^107、一雨降雨の確率モデルの開発^130、時間単位と水文量の分布特性を理論的に説明しようとする研究成果等^107,108があげられます。

確率降雨の決定においては、降雨継続時間$t$(一般には10分~2日)と降雨強度$I$の関係、すなわち、降雨強度曲線が必要となります。降雨強度曲線式としては、次の3形式がよく使われます。

Talbot型 $I = \frac{a}{t+b}$ (4.39) Sherman型 $I = \frac{a}{(t+b)^n}$ (4.40) 石黒型 $I = \frac{a}{b+t}$ (4.41)

ただし、$a,b,n$は地域によって異なる定数とされます。これらの式形は理論的に誘導されたものではなく、定数も経験的に決定されているといって過言ではありません。

長尾^114は確率分布モデルとして、式(4.12)のガンマ型分布($c=0$)を採用し、さらに降雨量時系列にAR(1)過程を導入することによって、Sherman形式の降雨強度曲線に理論的裏づけを行っています。例えば、日降雨量の自己相関係数が零と仮定して、確率日降雨量の時間配分問題を考えてみます。このとき24時間雨量$R_{24}$の任意の$T$時間雨量$R_T$への配分率は次式の理論的推定式として誘導されています。

$\frac{R_T}{R_{24}} = \left(\frac{T}{24}\right)^{\alpha}$ (4.42)

式(4.42)中の指数$\alpha$は所要確率年(あるいは超過確率)と式(4.12)のガンマ分布の形状母数$b$のみの関数です。確率統計学的考察から得られた指数$\alpha$は、わが国で経験的に知られている$\alpha=1/2$~$1/3$とほぼ対応しています。降雨量の時系列特性を考慮した場合にも、降雨量配分率と降雨継続時間の関係がSherman形式と類似の式で表現できた点は注目に値します。

計画規模の降雨量が決まれば、降雨波形つまり降雨量の時間分布を設定する必要があります。その一つの方法は特定確率年の降雨強度-継続時間$I-T$曲線からピークをはさむ各継続時間の平均降雨強度が等確率年になるようにハイエトグラフを作成する方法です^131。この方法においては、降雨強度曲線が各継続時間の年最大値標本値から独立に確率評価されているために、一雨期間内の雨量相互間の相関関係を再現できない問題が起きます。しかも、上記の方法で設定した計画降雨波形の確率年は、実際の降雨波形の確率年よりも大きい値となることがあるので注意を要します。

他の一つは既往の豪雨パターンに基づく引き伸ばし法^121です。この方法は実績降雨波形を重視しているものの、引き伸ばし率の限界について明確な根拠がありません。

最近、端野^132は2変数指数型確率分布に基づいて、所要確率年の総雨量に対するピーク時間雨量を決定する、いわゆる条件付確率降雨強度曲線と条件付確率波形を提案しています。この方法によれば1時間単位の雨量資料から5分程度までの短時間降雨強度を合理的に推定でき、しかも、実績降雨引き伸ばし限界の確率論的評価にも利用可能です。

多変数統計理論の適用

水工計画では、例えば多地点間の水文量の関係、本支川問題など、多変数の確率的手法の導入が要求されるような問題が非常に多くなってきています。多変数正規分布は、現在の確率統計水文学で重要な位置を占めていることは否定できません。また、前節までに紹介してきたように、実際面の要求から2変数ガンマ分布理論の応用研究もかなり進展してきています。さらに、多変数最大エントロピー分布の提案もされています^133。

水工計画において、多変数統計理論がどのように実務面に適用可能かを以下に大別して示します^111。

(1) 多変数からなる結合確率の評価 計画基本量としての水文量が複数個存在し、その生起の仕方が計画上安全であるか危険であるかの状態が判別できる場合、各状態に対する水文量の同時生起の可能性を確率評価するような事例です。その適用例をいくつかあげます。

① 長尾^134は洪水の河道ピーク水位と総流量が2変数対数正規分布に従うとして1河道、1貯留施設で構成されるモデル水系の洪水処理計画のあり方を考察しています。

② 栗田・江藤^135は、2支川のピーク流量の結合確率密度関数が2変数指数分布に従うとして、両者の相関が完全独立の場合と完全従属の場合について、支川の疎通能力と本川の氾濫危険度の関係を非常に簡単な式形の理論解で表現しています。

③ 室田ら^136は、貯留施設と排水施設間の治水機能分担のあり方を評価するために、等危険度線という概念を提案しています。この概念はピーク降雨強度と総雨量、またはピーク降雨と降雨継続時間の結合確率分布を基礎として構築されています。

(2) 指定変量による条件付変量の推定 ある変量の値が既知あるいは指定された場合に、他の変量の大きさや出現の可能性を確率的に推定する問題です。回帰分析による水文現象の時空間構造の推定、あるいは水文量の時空間相関特性を考慮した多地点水文量の模擬発生手法がその応用例です^111。なお、この問題では多変数正規分布が多用されます。

(3) 多変量群から類似小数群の判別 例えば、流域内の多数の雨量観測所の観測値を整理・統合するといった問題に相当します。成分分析、因子分析、判別解析等の多変量解析手法は、近年の電子計算機の発達・普及により広く利用されるようになってきています。