ECOH/YG 技術論文 No. 10

2009年5月7日版

高波の極値統計解析に対する L-moments 法の適用と分布関数の選択について
Application of L-moments Method to Statistical Analysis of Extreme Wave Heights andDiscussion on Selection of Distribution Functions

合田良実1・久高将信2
Yoshimi GODA and Masanobu Kudaka
1

土木学会名誉会員 工博 横浜国立大学名誉教授
2

土木学会会員(株) エコー（〒110-0014 東京都台東区北上野 2-6-4）(株) エコー（〒110-0014 東京都台東区北上野 2-6-4）

Application of the L-moments method for statistical analysis of extreme wave heights extracted by the
POT method is made for storm wave data at eight stations around Japan, which have the record length of
27 to 38 years. The formulas of the L-moments of the Weibull distribution are newly derived for
estimation of the shape, scale and location parameters from a sample. A new index called TUD (TwentyUp Deviation) is developed for judgment of the degree of goodness of data fitting to a candidate
distribution. Among eight stations, two stations show the best fitting to the Weibull distribution, four to
the General Pareto (GPA) distribution, and two to the General Extreme Value (GEV) distribution. At two
stations fitted to the GPA distribution, the theoretical upper bounds of extreme wave heights are only 7%
to 11% greater than the observed maximum wave height.

Key Words: Extreme wave statistics, GEV, GPA, Weibull distribution, L-moments, NOWPHAS, TUD,

１．まえがき374頁) のアプローチと類似の目的を持っている。

L-moments 法を用いた地域頻度解析についてはこれを解説した書籍を Hosking・Wallis (1997) が出版しており，
ま た 数 理 統 計 プ ロ グ ラ ム の デ ー タ ベ ー ス で ある
http://lib.stat.cmu.edu の中の lmoments から計算を実行するための FORTRAN プログラムをダウンロードできるようになっている。水文分野以外では外山・水野(2002) が全国アメダス地点の確率降水量の解析に用いており，高波については van Gelder ほか (2000) がオラ
ンダ沿岸の波浪データをこの手法を使って解析している。
本論文では L-moments 法を用いた極値統計解析の方法について解析し，NOWPHAS のわが国沿岸の波浪データベースから抽出した高波の資料について適用した事例を紹介する。ただし，Hosking・Wallis (1997) が推奨している L-moments 法には，高波の極値統計で多用さ極値統計解析にはいろいろな手法が開発されてきており，極値データの標本への当てはめ方法についてもProbability Weighted Moments (PWM，確率重み付き積率法) が Greenwood ほか (1979) によって提案され，米国における河川流量解析等に使用されてきた。さらに，
Hosking (1990) はこの PWM 法を発展させた L-moments法を発表している。米国においては河川流量の解析では単一河川のみを扱うのではなく，特性をほぼ同じくする流域内の複数の河川が類似の極値分布関数に従うと見なして解析する手法が主流であり，その際には L-moments法を用いた研究が多い。地域ごとに共通の極値分布関数を当てはめて解析する方式は地域頻度解析 (Regional
Frequency Analysis) と名付けられており，広域の海岸ごとに高波の母分布関数を見いだそうとする合田 (2008,-1-

れてきたワイブル分布が取り上げられていない。これまの順序統計量から次のような r 次平均量 br を計算する。
1 n(10)
b0 = ∑ x j
n j =1での各港湾における高波の極値解析資料との連続性から
みて，L-moments 法を適用するに当たってもワイブル分
布を含める必要がある。そこで，本論文ではワイブル分布に対する L-moments の理論値を示し，その適用方法を紹介する。

２．L-moments の定義とその計算

b1 =1 n ( j − 1)xj∑
n j = 2 (n − 1)

(11)

b2 =

1 n ( j − 1)( j − 2)xj∑
n j =3 (n − 1)(n − 2)

(12)

br =

1 n ( j − 1)( j − 2) K ( j − r )xj∑
n j =3 (n − 1)(n − 2) K (n − r )

(13)

一般に，

L-moments というのは，順序統計量の線形結合(linear combinations of order statistics) として Hosking(1990) が名付けたものである。いま，確率変量 X の分布関数を F (X ) とし，X がある閾値を超えない確率が Pであるような閾値を x (P ) と表示する。X (P ) は X のquantile function と呼ばれ，波高の極値統計でいうところの再現確率統計量である。これを数式で表示すれば
F[ x( P)] = P, x = F −1 ( P)

参考として，lmoments のデータファイルの中の年最
大風速の一つの標本（フロリダ州タンパ）について計算した結果を表−1 に示す。左から第 2 欄が各年の最大風
速，第 3 欄が順序統計量に並べ替えたものであり，第 4〜6 欄が br を計算するための行である。br の計算結果は最下行に示されている。

(1)

この確率統計量すなわち quantile function を用いて，ま表−1 n = 10 の場合の r 次平均量 br の計算表

ず次の確率重み付き積率 βr を定義する。
1

β r = ∫ x ( P ) P r dP

j
1
2
3
4
5
6
7
8
9
10

(2)

0

L-moments はこの積率 βr を用いて次式のように定義される。

λ1 = β0(3)λ2 = 2β1 – β0(4)

λ3 = 6β2 – 6 β1 + β0(5)λ4 = 20β3 – 30 β2 + 12 β1– β0(6)

さらに，上記のL-moments から相互の比率を次の用に定b0 = xj

b1

b2

b3

45
47
65
50
56
55
37
53
44
42

37
42
44
45
47
50
53
55
56
65
49.4

4.666667
9.777778
15.00000
20.88889
27.77778
35.33333
42.77778
49.77778
65.00000
27.1

1.222222
3.750000
7.833333
13.88889
22.08333
32.08333
43.55556
65.00000
18.94167

0.535714
2.238095
5.952381
12.61905
22.91667
37.33333
65.00000
14.65952

この r 次平均量 br を積率 βr の普遍推定値であり，標本
L-moments の値は，br を用いて次のように計算する。

義する。

L-CV (Coefficient of L-variation)： τ = λ 2 / λ1x(7)l 1 = b0(14)

l 2 = 2b1 − b0(15)

1 次の L-moment λ1 は母集団の平均値を与え，2 次の L-

l 3 = 6b 2 −6b1 + b0

(16)

moment λ2 は標準偏差に関するものの偏差そのものでは

l 4 = 20b3 − 30b 2 +12b1 − b0

(17)

L-skewness (L-ひずみ度)：
L-kurtosis (L-尖鋭度)：

τ 3 = λ3 / λ 2 (8)τ 4 = λ 4 / λ 2 (9)0.5

なく，母集団の標準偏差を σ で表すと 0.5σ ≤ λ2 ≤ σ /3

そして，この標本 L-moments を用いて次のような標の範囲の値をとる (Hosking・Wallis (1997), p. 35 より推

本 L-moments 比を算定する。

定)。
式 (2)〜(9) は母集団に関するものである。いま，高波の観測・推算資料からからある閾値を超える極値時系列
データ (POT) が n 個，すなわち大きさが n の標本が得られたとすると，これを昇順に並べ替えて順序統計量に変

L-CV (Coefficient of L-variation)： t = l 2 / l 1(18)

L-skewness (L-ひずみ度)：

t3 = l 3 / l 2(19)

L-kurtosis (L-尖鋭度)：

t 4 = l 4 / l 2 (20)表−1の例では，標本L-momentsと標本L-moments比が換する。すなわち， x1 ≤ x 2 ≤ K ≤ x n と並べ替える。こ2

(C) 一 般 化 極 値 分 布

次のように算定される。
l 1 = 49.4 , l 2 = 4.8, l 3 = 0.450, l 4 = 0.7405,(Generalized Extreme-ValueDistribution: GEV)
分布関数：

t = 0.0972, t 3 = 0.0938, t 4 = 0.1543⎧exp[−{1 − k ( x − B ) / A}1 / k ], k > 0 : − ∞ ≤ x < B + A/k⎪
F ( x) = ⎨k < 0 : B − A/k < x < ∞⎪exp[− exp{−( x − B ) / A}],k = 0: − ∞ < x < ∞⎩なお，標本に対して l1, l2, …, t, t3, t4, …などの記号を用いるのは，母分布関数に対する推定値であることを明示するためである。

(29)
確率統計量 (quantile)：

３．極値分布関数と L-moments 比の関係

⎧ B + A{1 − (− ln P ) k } / k , k ≠ 0
x = F −1 ( P ) = ⎨k =0⎩ B − A ln(− ln P ),Greenwood ほか (1979) による PWM 法は，標本について確率を重みとした積率を計算することによって適合候
補の分布関数の形状・尺度・位置母数を速やかに推定することを目的としていた。L-moment 法においても，各種の極値分布関数について L-moments の理論値を計算
し，分布関数の母数との関係を求めておくことによって分布関数の母数推定が可能となる。
Hosking (1990) は式 (2)〜(9) を用いて各種の極値分布関数に対するλ1, λ2, τ3, および τ4 の理論値ならびに母数の推定式を計算した。この結果は Hosking・Wallis (1997)の解説書にも付録としてまとめられている。そのうちから，高波の極値統計解析に使用する分布関数に関わるものを以下に再録する。ただし，ワイブル分布についてはGreenwood ほか (1979) を参照して今回新しく導いたも
のである。

(30)

L-moments：
λ1 = B + A{1 − Γ(1 + k )} / k , λ2 = A(1 − 2− k ) Γ(1 + k ) / k ,τ 3 = 2 (1 − 3− k ) /(1 − 2− k ) − 3,
τ 4 = {5 (1 − 4− k ) − 10 (1 − 3− k ) + 6 (1 − 2− k )} /(1 − 2− k )(31)
母数推定式：k ≈ 7.8590c + 2.9554c 2 , c =A=

λ2 k(1 − 2 − k ) Γ(1 + k )

log 2
2
−
3 + τ 3 log 3(32), B = λ1 − A {1 − Γ(1 + k )} / k

一般化極値分布は，式 (29) の表現では k < 0 のときに
極値 II 型分布，k > 0 のときに極値 III 型分布となり，前者では確率統計量が下限値を持ち，後者では上限値を持(A) 指数型分布

つ。また，k = 0 のときは極値 I 型分布である。ただし，分布関数：
F ( x ) = 1 − exp{−( x − B ) / A},

高波の極値統計で使われている極値 II 型は式 (29) の形

B≤x<∞

(21)

状母数を k’ = –1/k と変換したものに等しい。

確率統計量 (quantile)：−1

x = F ( P) = B − A ln(1 − P)
L-moments：

なお，式 (32) における形状母数 k は近似推定式であり，

Hosking (1990) が導いたものである。

(22)(D) 一 般 化 パ レ ー ト 分 布Distribution: GPA)
分布関数：

λ1 = B + A, λ 2 = (1 / 2) A, τ 3 = 1 / 3, τ 4 = 1 / 6 (23)母数推定式：
A = 2λ 2 ,B = λ1 − A(24)(Generalized Pareto⎧1 − [1 − k ( x − B) / A]1 / k 1 ≥ k > 0 : B ≤ x < B + A/k⎪
F ( x) = ⎨k < 0:
B≤x<∞⎪1 − exp{−( x − ξ ) / α }, k = 0 :B≤x<∞⎩(B) ガンベル分布（極値Ⅰ型分布）
分布関数：
F ( x ) = exp[ − exp{−( x − B ) / A}], B ≤ x < ∞ (25)(33)

確率統計量 (quantile)：−1

x = F ( P) = B − A log(1 − ln P)

確率統計量 (quantile)：(26)

L-moments：
λ1 = B + Aγ , λ 2 = A ln 2, τ 3 ≅ 0.1699, τ 4 ≅ 0.1504 (27)

⎧ B + A{1 − (1 − P) k } / k , k ≠ 0
x = F −1 ( P) = ⎨k =0⎩ B − A ln(1 − P),

母数推定式：
A = λ2 / ln 2,

B = λ1 − Aγ

(34)

L-moments：
λ1 = B + A /(1 + k ), λ 2 = A /{(1 + k )(2 + k )},τ 3 = (1 − k ) /(3 + k ), τ 4 = (1 − k ) ( 2 − k ) /{(3 + k ) (4 + k )}(28)

ここに γ はオイラーの定数, 0.5772....である。

(35)

3

母数推定式：k = (1 − 3τ 3 ) /(1 + τ 3 ),

3.5

(36)

B = λ1 − ( 2 + k ) λ 2

Weibull Distribution
3.0
Shape Parameter, k

A = (1 + k ) ( 2 + k ) λ 2 ,一般化パレート分布は一般化極値分布のような極大・極小値に対する漸近関数ではないけれども，ポアソ
ン分布に従う確率変量のうち，ある閾値を超えるものに着目する場合などを表示するのに適切な分布関数といわれる。北野ほか (2002) は Kodiak 沖の高波データについ

k-cal
k-est

2.5
2.0
1.5
1.0

てこの分布の適用性を検討している。なお，この分布関数の特性等については Coles (2001) などを参照されたい。

0.5
0

0.1

0.2

(E) ワイブル分布

F ( x) = 1 − exp[−{( x − B) / A}k ], B < x < ∞が数値計算結果から導いた近似式である。すなわち，形(38)

状母数 k をいろいろ変えて 3[次の積率 λ3 を計算し，そのL-moments：

結果に対して 6次の多項式を当てはめたものである。近λ1 = B + AΓ(1 + 1 / k ), λ 2 = A(1 − 2 −1 / k ) Γ(1 + 1 / k ),−1 / k−1 / k) /(1 − 2

−1 / k) − 10 (1 − 3

似精度は ±0.2% である。

),

−1 / k図−2 はワイブル分布における母数 k と L-skewness の) + 6 (1 − 2−1 / k)} /(1 − 2−1 / k)関係を示している。なお，この図では数値計算の結果と式 (40) による推定値の両方を示しているが，図上では(39)

両者の違いを見分けることがむずかしい。

母数推定式：
k = 285.3 λ36 – 658.6 λ35 + 622.8 λ34 – 317.2 λ33+ 98.52 λ32 – 21.256 λ3 + 3.5160

λ2(1 − 2−1 / k)Γ(1 + 1 / k )４．L-moments 比に基づく極値分布関数の選別(40)

前節に紹介した極値分布関数のうち，指数分布と極値(41), B = λ1 − A Γ(1 + 1 / k )I 型分布は 2 母数型分布であって， L-skewness と Lkurtosis は一定値である。他の分布関数は 3 母数型分布であるため，形状母数 k の値によってL-skewness とLkurtosis の値が連続的に変化し，両者の関係は分布関数に特有な一つの曲線で表される。ここでは紹介していな
いが，4 母数型分布分の場合には，両者がある面内で変化する関係を示す。
図−3 は，前節の分布関数に加えて Pearson III 型分布
も含め，二つの L-skewness 比の関係をプロットしたものである。2 母数型分布の指数型分布とガンベル分布はそれぞれ一つの点で表される。指数型分布の点（四角にクロス）から左下へ延びる曲線は一般化パレート分布(GPA) で形状母数が k > 0 の場合を表し，同じく右上へ上昇する曲線は k < 0 の場合である。なお，k > 0 の場合には前述のように確率統計量 x の値に上限がある。式(33) の関数形は Hosking・Wallis (1997) によるもので，一般の数理統計学では形状母数の符号が逆であるのが普通である。
ガンベル分布 (FT-I) の点（白十字）の左へ緩やかに

極値分布関数の特性を決めるのは形状母数 k である。

L-moments 法では L-skewness の値によって母数 k の値が一義的に定まる。図−1 は GPA と GEV について母数
k と L-skewness の関係を示したものである。
1
GEV k
GPA k

Shape Parameter, k

0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
0

0.6

ワイブル分布の形状母数 k を推定する式 (40) は著者ら

x = F −1 ( P) = B + A[{− log(1 − P)}1 / k ]

A=

0.5

(37)

確率統計量 (quantile)：τ 4 = {5 (1 − 4

0.4

図−2 ワイブル分布における形状母数 k とL-skewness τ3 との関係

分布関数：

τ 3 = 3 − 2 (1 − 3

0.3

L-skewness, τ3

0.1

0.2
0.3
0.4
L-skewness, τ

0.5

0.6

3

図−1 GPA および GEV 分布における形状母数 k とL-skewness τ3 との関係

4

下がる線は一般化極値分布 (GEV) のうちで確率変数 x

いる。すなわち，それぞれの標本に含まれているデータの値に上限がある場合に相当し，極値 III 型分布 (FT-III)

数の少なさを複数地点全体のデータ数で補う方式である。

である。右へ緩やかに上がる線は極値 II 型分布 (FT-II)

地域頻度解析は年最大値資料を対象としており，高波である。ワイブル分布は極値 III 型分布と上限値があるのように極値時系列法 (POT) で抽出されたデータに対一般化パレート分布の間に位置するが，ワイブル分布はしてどのように地域頻度解析を行うかは今後の課題であ

上限値がなく，右裾がどこまでも延びるのが特徴的である。本論文ではまず単一地点の極値資料に対してあらかる。

じめ分布関数の候補を選定し，L-moments を用いて形
0.5

0.4

L-Kurtosis, τ4状・尺度・位置の 3 母数を推定する方式について検討を

GEV (k = 0)
GEV (k < 0)
GEV (k < 0)
GEV (k > 0)
Weibull
Weibull
Exponential
GPA (k > 0)
GPA (k < 0)Pearson III

0.3

進める。以下においては，この方式を LMN 法と略称する。
k = - 0.4

k = - 0.2

0.2

５．使用する高波資料の概要

k = - 0.3

先に筆頭著者は NOWPHAS の1989年までのデータ

k = 0.75

k = - 0.1

ベースを使用し，日本沿岸の高波の極値分布関数についk = 0.85て考察した（合田 2008, pp. 374-375）。それから既に10

k = 1.0
k = 1.2

0.1

k = 2.0

年を経過し，波浪データベースも拡充されたところから，

k = 1.7 k = 1.4

長期間にわたる波浪データが蓄積されている 8 地点を選
0

0

0.1

0.2

0.3

0.4

0.5

び，POT 法で抽出した高波資料の解析を行った。これ

0.6

らの高波資料の概要を表−1 に示す。

L-Skewness, τ3

図−3 各種分布関数における L-skewness と L-kurtosis の関係

表−1 検討に用いた高波資料の概要

図中のワイブル分布と極値 II 型分布については，形状

デー

波高閾値

タ数

(m)

(m)

38

923

3.0 (0.5) 5.0

10.65

地点名

観測期間

年数

酒田

1970-2008

母数に特定の値を与えたときのL-skewness とL-kurtosisの対応点を示してある。また，ワイブル分布と Pearson

III 型分布は L-moments 比の値が隣接している。すなわ
ち，後者はワイブル分布で代替させることが可能であり，極値統計の標本当てはめの候補に入れなくても良いと考えられる。
ワイブル分布と一般化パレート分布 (GPA) は指数分布を共有している。指数分布の点の左下側では，両者のL-kurtosis の値に差異が見られる。すなわち，ワイブル
分布では L-kurtosis の値が余り変化しないのに対し，一般化パレート分布では L-skewness の減少につれて Lkurtosis も減少する。指数分布の点の右上側ではワイブル分布と一般化パレート分布の違いが小さいが，後者のほうが L-kurtosis がやや大きめの値を示す。
このように，L-skewness と L-kurtosis の二つのパラメータを用いることによって，極値分布関数の差異を判別できる可能性が考えられる。図−3 に示したのは母集団についての結果であり，データ数に限りのある標本の場合には統計的変動性の影響によって母集団ほどの差別化はむずかしいであろう。
地域頻度解析の場合には，複数地点の極値資料を収集し，これらが同じ母集団に属することを検定した上で，複数地点全体に共通な分布関数を見いだす方式を用いて5(H1/3)max

金沢

1970-2008

37722

3.0 (0.5) 5.0

8.14

むつ小川原

1971-2008

37551

2.5 (0.5) 5.0

9.56

鹿島

1972-2008

36598

2.5 (0.5) 5.0

7.50

波浮

1973-2008

35640

3.0 (0.5) 5.0

8.48

志布志

1980-2008

28482

1.5 (0.5) 4.0

10.30

中城湾

1975-2008

27599

2.0 (0.5) 5.0

12.08

那覇

1973-2008

35656

2.5 (0.5) 4.5

9.24

観測年数は1年単位の欠測期間を除外したものである
が，月単位の欠測期間の調整は行っていない。また，合田・小長谷 (1998) が行ったような，高波のピーク時の部分欠測も補填していない。また，気象原因別（台風系と低気圧系など）のデータ区分もおこなっていない。したがって，この資料だけでは確率波高の推定作業を進めることはできない。
表中のデータ数は，波高閾値が最も低いときのものであり，それから 0.5 m ごとに閾値を高めていったので，解析を行ったデータ数は閾値ごとに異なる。

６．LMN 法における最適合分布の判定
極値統計解析では，与えられた標本に対して複数の分布関数を当てはめ，そのうちで適合度の最も高い関数を仮に母集団分布と見なす方式が一般的である。最小 2 乗

異なる最適合分布が得られてしまう。

法の場合には，確率変量とその基準化変量との間の相関

たとえば，図−6 は中城湾の高波データについて波高係数が高いほど適合度が良いとみなす（これを若干修正の閾値を Hc = 2.0 (1.0) 5.0 m と 4 通りに変えたときの元したのが合田が提案した MIR 基準である）。最尤法では尤

の波高データと，確率統計量として推定された波高値と度を判定基準に使うけれども，最適関数判定の感度がそを比較した図である。データは波高の大きさの順に並べれほど良いものではない。

られておりこうした表示法は Quantile – Quantile (Q – Q)

LMN 法の場合には，まず L-skewness の値によって分
布関数ごとの形状母数が算定される。そうすると，その分布関数の L-kurtosis の理論値が計算されるので，その値と標本から求めた値との差を判定基準とする方式が一つ考えられる。図−5 は，沿岸 8 地点のうちの 4 地点について波高閾値を変えることによって，L-skewness とL-kurtosis の関係がどのように変わるかも示したものである。

プロットと呼ばれる。ここでは当てはめ分布関数をワイブル分布に固定しているが，波高閾値によって形状母数が変化している。
16

Predicted Wave Height, Hpred (m)L-kurtosis, τ4

0.2

0
0.1

GEV
Weibull
GPA
SakataKanazawa
Nakagusuku
Naha
0.2

0.3

0.4

Hpred = Hobs

14

0.3

0.1

Naka-gusuku Wan
Hc = 2.0 m (k = 0.78)
Hc = 3.0 m (k = 1.00)
Hc = 4.0 m (k = 1.23)
Hc = 5.0 m (k = 1.40)

12

10

8

6

0.5

L-Skewness, τ34

図−5 波高閾値を変えたときの沿岸波浪の L-skewness

6

8

10

12

Observed Wave Height, Hobs (m)とL-kurtosis の関係

図−6 中城湾のデータに波高閾値を変えてワイブル分布を

酒田港では Hc = 3.0 m のデータが左下に位置し，閾値当てはめたときの順序統計量と確率統計量の比較(Q–Q プロット)を上げるにつれて表示点がワイブル分布の曲線に沿って右上に移動する。金沢港では Hc = 5.0 m のデータが一つ

LMN 法では，式 (1) 〜(4) の r 次平均量 br を計算する

だけ左上に離れており，他の閾値のデータは表示点がGPA の曲線上の一カ所に固まっている。中城湾の資料
では，Hc = 2.0 m のデータが右上に位置し，閾値を上げるにつれて表示点がワイブル分布の曲線の下方を左下へ向けて大きく移動し，Hc = 5.0 m のときには GEV の関係
曲線の左上へ飛び跳ねる。那覇港の場合には，Hc = 3.0
m のデータが左下に位置し，閾値を上げるにつれて表示点が GEV の曲線へ向かって大きく右上に移っていく。
波高の閾値によって適合分布関数が変化することはあり得るけれども，図−5 に示す中城湾や那覇港のケースは変化が極めて大きく，単一の極値資料に対して Lskewness と L-kurtosis の関係だけに基づいて最適合の分布関数を選ぶ方式には疑問を抱かせる。実際に表−1 の高波資料について，標本から得られた L-kurtosis の値とその理論値との差を分布関数ごとに比較し，最低値を示
す分布関数を調べてみると，波高の閾値を上げるたびに

6

際に，データを小さい方から昇順に並べ替えた順序統計量を使う。閾値を上げたときにはその閾値よりも小さな
波高データは切り捨てられ，それよりも大きなデータの順位が切り下げられる。n 個のデータの中で小さい方から j 番目，あるいは大きい方から m 番目のデータに対す
る確率統計量を求めるときには，その非超過確率 Pj‫׃‬n を与える必要があり，このために次のプロッティング・ポ
ジション公式を使う。

j − 0.35m − 0.65(42), Pm:n = 1 −nn
式 (42) の非超過確率 Pj::n は GPA 分布に対してHosking・Wallis (1997) が推奨 (33頁) しているもので，他の分布関数については問題があるかもしれないと述べている。しかし，他に適切な公式が提示されていないので，本論文では Q – Q プロットを行うときの非超過確率
Pj:n =

として式 (42) を使うことにする。なお，Pm:n は n = j +布関数ごとの適合度の順位は変わらなかった。

m – 1 の関係から導かれたものである。

中城湾の高波極値データについて TUD による適合度

最小 2 乗法の場合には，プロッティング・ポジションを計算すると図−7 のようになり，高波抽出の閾値を上の取り方によって確率統計量に bias が生じることがあげるに従って適合度が高まることが分かった。図−7 でるため，公式の選定には十分な注意が必要であるけれどは，当てはめ候補の極値分布関数として3 母数型分布でも，LMN 法では分布関数の当てはめ自体にはプロッ

ある GVE, GPA, およびワイブル分布の 3 種類を取り上ティング・ポジション公式の選択は影響を与えない。

げている。３節では 5 種類の関数を紹介したけれども，そのうちの指数型分布は GPA 分布およびワイブル分布

図−6 では，あらかじめ分布関数としてワイブル分布を固定した条件で，極値波高データを抽出する閾値の影に含まれているので，上記の 3 分布関数を用いれば良い

響を調べたものである。図から明らかなように，閾値のと考えたものである。先に図−5 に示した中城湾のL-設定によって確率統計量が観測された波高の順序統計量

skewness とL-kurtosis の関係の変化を図−7 の TUD の値の変化から読みとることはむずかしい。
高波抽出の閾値によって適合度が変化するのは，中城湾の場合には台風系と低気圧系という母集団の異なる高波が混在しているためである。合田・小長谷 (1978)の場合には抽出された高波の全てについて気象原因を確認し，台風系とそれ以外の高波に区分して解析したけれ
ども，今回はそうした吟味を行っていない。波高の閾値が低い場合には波高が比較的に小さい非台風系のデータが数多く含まれ，閾値を上げるに従って台風系の高波のみが解析対象となったものであろう。図−7 によれば，
閾値を Hc ≥ 4.0 m とすることによって，TUD の値が安定するので，そうしたデータを対象として中城湾の高波の確率波高を計算すればよいと思われる。

と一致する度合いが相当に変化する。
設計波高を選定するために極値統計解析を行う場合には，主として波高の上位のデータに対する極値分布関数の適合度に関心が持たれる。そこで，まず極値データの観測値の順序統計量と確率統計量の相対偏差 (%) を次のように定義する。

⎛ H pred, m − H obs, m ⎞⎟ ×100
Em = ⎜⎜⎟
H obs,m⎝⎠(43)

ここに，Hpred, m は上から第 m 番目の確率統計量，Hobs, m
は観測波高の第 m 番目の順序統計量である。そして，上位 20 個の 2 乗平均偏差を次のように定義し，それが最小となるものを最適合とすることを提案する。

TUD =1 20 2∑ Em
20 m =1(44)

７．各地点の高波の極値分布特性
本節では，表−1 に示した 8 地点の高波資料に対して記号の TUD は，Twenty-Up Deviation の頭文字を取った

LMN 法で極値分布関数を当てはめた結果について述べ略号である。20 個を取り上げたことに特別の理論的根拠はない。各種分布関数を当てはめたときの適合度を比14較し，最もよく適合するものを選ぶ判断基準として使う

Naka-gusuku Wan

ためのものである。上位 10 個についても試行してみた

Hpred = Hobs
Predicted Wave Height, Hpred (m)Twenty Up Deviation, TUD (%)

35

Nakagusuku Wan

30

Weibull
GPA
GEV

25

Weibull (k = 1.36)
GPA (k = 0.262)
GEV (k = - 0.086)

12

ところ，TUD の絶対値はやや大きくなったものの，分

20
15

10

8

6

10
5
0

4

1.5

2.0

2.5

3.0

3.5

4.0

4.5

5.0

6

8

10

12

Observed Wave Height, Hobs (m)5.5Threshold Wave Height, Hc (m)

図−8 中城湾の高波 (Hc = 4.5 m) に対する極値分布関数の当てはめ結果

図−7 中城湾の高波データに対する TUD 値

7る。図−8 は中城湾の高波データであって，図−7 の結りもわずかながら適合度に劣るけれども，Hc = 5.0 m の果を参照して高波抽出の閾値を Hc = 4.5 m と設定した場

場合の図−10 では両者の違いはわずかである。GEV 分合を示している。この高波データでは観測値の第 1 位が布は上位の順序統計量に対して過大評価，中位の順序統12.08 m，第 2 位が 11.93 m と近接しており，このため確率統計量の推定結果が前者に対しては過大，後者に対しては過小な結果となっている。 3 種類の分布関数は適合度に大差がなく，どの関数も適用可能である。GEV
分布は順序統計量の第 1 位に対してやや過大な値を与えるけれども，第 2 位以下に対しては他の二つの分布関数よりも適合度がよい。

計量に対して過小評価となり，酒田港には不向きである。
金沢港については，極値分布関数の適合度指標のTUD 値を図−11，閾値を Hc = 5.0 m と設定した場合の順序統計量と確率統計量の比較を図−12 に示す。

Twenty Up Deviation, TUD (%)

20

Twenty Up Deviation, TUD (%)

14

Sakata12

Weibull
GPA
GEV

10
8

Kanazawa
Weibull
GPA
GEV

15

10

5

6

0
2.5

4

3.0

2

4.0

4.5

5.0

5.5

図−11 金沢港の高波データに対する TUD 値

0
2.5

3.0

3.5

4.0

4.5

5.0

5.5

10

Threshold Wave Height, Hc (m)Kanazawa
Hpred = Hobs
Predicted Wave Height, Hpred (m)

図−9 酒田港の高波データに対する TUD 値
13

SakataPredicted Wave Height, Hpred (m)3.5

Threshold Wave Height, Hc (m)12

Hpred = Hobs

11

Weibull (k = 1.15)
GPA (k = 0.119)
GEV (k = -0.169)

10

Weibull (k = 1.47)
GPA (k = 0.329)
GEV (k = - 0.049)

9

8

7

9
6
6

8

7

8

9

Observed Wave Height, Hobs (m)図−12 金沢港の高波 (Hc = 5.0 m) に対する極値分布関数の当てはめ結果

7

6
6

7

8

9

10

11

Observed Wave Height, Hobs (m)金沢港の高波資料では，閾値を変えても GPA 分布の優位性が変わらない。ワイブル分布や GEV 分布を当て

図−10 酒田港の高波 (Hc = 5.0 m) に対する極値分布関数の当てはめ結果

はめると，上位の順序統計量に対して過大評価となってしまう。なお，高波抽出の波高閾値が Hc = 3.5 m のとき次に，酒田港の高波データに対する極値分布関数の適に適合度指標の TUD 値が増大する理由は解明できてい合度指標の TUD 値を図−9，閾値を Hc = 5.0 m と設定しない。ただし，閾値レベルによる適合度の変化が少ない

た場合の順序統計量と確率統計量の比較を図−10 に示ところから，金沢港の高波資料も同じ母集団に属するとす。酒田港の場合にはワイブル分布が最適合であり，閾

みなされよう。

値を変えても適合度はそれほど変化しない。すなわち，

次にむつ小川原港についての検討結果を図−13, 14 に酒田港の波浪データはそのほとんどが同じ高波の母集団示す。図−13 は適合度指標の TUD 値が波高閾値によっに所属すると推測される。GPA 分布はワイブル分布よ8て変化する状況を示している。ワイブル分布と GPA 分

分布は全般的に適合度が低く，鹿島港の高波に対する極布は波高閾値の影響をほとんど受けず，どちらかといえ値分布としては不適当である。他の地点では波高閾値をばワイブル分布のほうが適合度がよい。一方，GEV 分

上げるに従って TUD 値が小さくなる，すなわち波高

布は波高閾値が低いときには適合度が極めて低いけれど

データに当てはめる分布関数の適合度が高まるのであるも，波高閾値が Hc = 4.0 m 以上ではワイブル分布やが，鹿島港の場合には波高閾値が Hc = 3.5 m のときに

GPA 分布よりも高波のデータによく適合するようになる。図−14 は波高閾値が Hc = 5.0 m のときの順序統計量と確率統計量の比較を示したものである。第 1 位のデータについては GEV 分布がやや高めの推定値を与えるものの，波高が Hobs = 6.0〜7.0 m の範囲では GEV 分
布がデータによく適合しており，TUD 値が最小となる
ことを裏付けている。

TUD 値が最小となっている。
この波高閾値が Hc = 3.5 m のときの順序統計量と確率統計量の比較を示したのが図−16 である。図で明らかなように，GPA 分布が波高の観測データからの乖離が最も少なく，他の分布よりも TUD 値が小さいことに対応している。
20

KashimaMutsu-Ogawara

12

Twenty Up Deviation, TUD (%)

Twenty Up Deviation, TUD (%)

14
Weibull
GPA
GEV

10
8
6
4
2

10

5

0

0
2.0

2.5

3.0

3.5

4.0

4.5

5.0

Weibull
GPA
GEV

15

2.0

2.5

3.0

3.5

4.0

4.5

5.0

5.5

Threshold Wave Height, Hc (m)5.5Threshold Wave Height, Hc (m)

図−15 鹿島港の高波データに対する TUD 値

図−13 むつ小川原港の高波データに対する TUD 値
Kashima11

9

Mutsu-Ogawara
Predicted Wave Height, Hpred (m)Predicted Wave Height, Hpred (m)H pred = H obs

10

Weibull (k = 1.08)
GPA (k = 0.064)
GEV (k = - 0.201)

9

8

7

6

5

H pred = H obs
Weibull (k = 1.16)
GPA (k = 0.126)
GEV (k = - 0.164)

8

7

6

5

5

6

7

8

9

5

10

Observed Wave Height, Hobs (m)67
Observed Wave Height, Hobs (m)8

図−14 むつ小川原港の高波 (Hc = 5.0 m) に対する極値分布
関数の当てはめ結果

図−16 鹿島港の高波 (Hc = 3.5 m) に対する極値分布
関数の当てはめ結果

鹿島港についての検討結果は，図−15, 16 に示す通り

その次は波浮港であり，このデータに対する検討結果である。図−15 は適合度指標の TUD 値が波高閾値には図−17, 18 に示す通りである。図−17 は適合度指標のよって変化する状況を示しており，ワイブル分布とTUD 値が波高閾値によって変化する状況を示しており，
閾値レベルが上るにつれて TUD 値が急速に減少している。中城湾の場合と同様に，波高レベルの低い高波のグ

GPA 分布は波高閾値の影響をあまり受けず，どちらかといえば GPA 分布のほうが適合度がよい。一方，GEV

9

GEV の適合度が劣り，波高閾値が低い間はワイブル分布が優位であり，波高閾値が高くなると GPA 分布がやや優位に立つ。

ループと高いグループとで，高波を発生させた気象原因が異なる可能性が強い。極値分布関数として GEV 分布
は全閾値レベルで適合度が低い。ワイブル分布と GPA
分布とでは，波高閾値が低いときは前者の適合度が高く，

35

Twenty Up Deviation, TUD (%)

波高閾値が高まると後者のほうが適合度が高くなる。
図−18 は波高閾値が Hc = 5.0 m のときの順序統計量と確率統計量の比較を示したものである。GPA 分布が高波のデータによく適合していることが認められる。
30

Twenty Up Deviation, TUD (%)

Habu
25

Weibull
GPA
GEV

20

Shibushi30

Weibull
GPA
GEV

25
20
15
10
5
0
1.0

1.5

2.0

2.5

3.0

3.5

4.0

4.5

Threshold Wave Height, Hc (m)15

図−19 志布志港の高波データに対する TUD 値

10
5

12

Shibushi
3.0

3.5

4.0

4.5

5.0

Hpred = Hobs

5.5

Predicted Wave Height, Hpred (m)02.5

Threshold Wave Height, Hc (m)

図−17 波浮港の高波データに対する TUD 値

Habu

Predicted Wave Height, Hpred (m)10

Hpred = Hobs
Weibull (k = 1.39)
GPA (k = 0.281)
GEV (k = - 0.075)

9

Weibull (k = 1.54)
GPA (k = 0.370)
GEV (k = - 0.026)
9

6

3
8

4

6

8

10

12

Observed Wave Height, Hobs (m)図−20 志布志港の高波 (Hc = 4.0 m) に対する極値分布
関数の当てはめ結果

7

6
8

Naha
5

6

7

8

Twenty Up Deviation, TUD (%)

5
9Observed Wave Height, Hobs (m)図−18 波浮港の高波 (Hc = 5.0 m) に対する極値分布
関数の当てはめ結果

さらに図−19, 20 には志布志港の高波データに対する検討結果であり，図−19 は適合度指標の TUD 値の変化，
図−20 は波高閾値が Hc = 4.0 m のときの順序統計量と確Weibull
GPA
GEV

6

4

2

0
2.0

2.5

3.0

3.5

4.0

4.5

5.0

Threshold Wave Height, Hc (m)率統計量の比較を示している。

図−21 那覇港の高波データに対する TUD 値

波高閾値による極値分布関数の適合度の変化は波浮港のデータと類似しており，閾値レベルがあがるにつれて最後に，那覇港の高波データについての検討結果を図

適合度が一様に向上する。3 種類の分布関数の中では−21, 22 に示す。図−21 は適合度指標の TUD 値の変化

10

4 地点のうち，金沢港では波高が Hupper = 8.73 m で頭打ちになると予測され，既往最大観測値の 8.14 m との差が僅か (7%) しかない。また，志布志港でも (Hobs)max =10.30 m に対して Hupper = 11.48 m とその差が小さい(11%)。また，両地点における高波の発生特性を考えると，発生波高がこれらの上限値で抑えられるということは考えにくい。
極値分布関数の当てはめの適合度だけで判定すると一般化パレート (GPA) 分布が最適合となっても，設計のための確率波高を選定する際には，上限値を持たないワ
イブル分布を当てはめることを優先すべきであると考えられる。この点については今後さらに検討したい。

10

Naha
Predicted Wave Height, Hpred (m)Hpred = Hobs
Weibull (k = 0.98)
GPA (k = - 0.015)
GEV (k = - 0.250)
8

6

4
4

6

8

10

Observed Wave Height, Hobs (m)図−22 那覇港の高波 (Hc = 4.5 m) に対する極値分布
関数の当てはめ結果

９．地域頻度解析のための今後の課題についてであり，他の諸港とは異なって GEV 分布の適合度がワ

今回の検討結果では，地点ごとに最適合と判定される

イブル分布や GPA 分布よりも高いケースが多く，特に極値分布関数が異なり，また POT 法で高波データを抽

波高閾値が Hc = 4.5 m のときには適合度が最も高くなっ出するときの閾値によっても最適合の分布関数が逆転すている。この状況は，図−22 の順序統計量と確率統計ることも少なくなかった。日本沿岸の設計波高を合理的量の比較によく現れている。

に設定するためには，高波の母集団についての検討をさらに高度なものとする必要がある。
最初に述べたように，米国の河川流量解析では母集団

８．GPA 分布における上限値の存在についてが同じと推定される複数の河川の流量データを総合的に今回の検討結果では，一般化パレート (GPA) 分布が解析する地域頻度解析の手法が用いられている。この方適合する地点がかなり多かった。波高閾値の選定にもよ法は毎年最大値を取り扱っており，各河川の流量データ

るが，適合度指標の TUD 値が最小であったのが 8 地点

をそれぞれの平均値で除して正規化した無次元流量を対では 4 地点を数えた。しかし，GPA 分布は式 (33) の表象として LMN 法を適用する。複数河川の無次元流量の式で形状母数が正の時に極値が上限値 xupper = B + A/k でデータに対して共通の極値分布関数が得られたならば，頭打ちになるという特性がある。そこで，この理論的上

河川ごとの平均流量を乗じることによって流量の極値分限値が観測最大波高に対してどの程度離れているかを調布を復元する。
これに対して高波の極値統計解析では，データ数を多べたのが表−2 である。

くして確率統計量の信頼区間を狭めるために POT 法を用いるのが標準であり，単純に地点ごとの高波データの表−2 GPA 分布を適用したときの波高の上限値
地点

Hc
(m)

N

k

平均値で正規化することには疑問が残る。POT 法にふA

B

Hupper(Hobs)max

(m)

(m)

(m)

(m)

15.18

10.65

また，今回の検討では波高の閾値レベルを上げること

酒田

5.0

244

0.1187

1.212

4.97

さわしい地域頻度解析法を開発する必要がある。

金沢

5.0

162

0.3294

1.236

4.98

8.73

8.14

によって分布関数の適合度が向上したけれども，それに

むつ小川原

5.0

63

0.0644

0.977

4.99

20.15

9.56

鹿島

3.5

214

0.1265

1.110

3.47

12.24

7.50

付随して利用できるデータの個数が減少している。今回

波浮

5.0

66

0.2811

1.547

4.87

10.38

8.48

志布志

4.0

42

0.3704

2.868

3.74

11.48

10.30

の幅はデータ個数の平方根に逆比例するので，データは

中城湾

4.5

88

0.2621

2.440

4.45

13.76

12.08

できるだけ多く残してデータ個数を多くしておきたい。

那覇

4.5

51

–0.0150.8331

4.42

∞

9.24

は信頼区間について検討していないけれども，信頼区間分布関数の適合度と信頼区間の幅の兼ね合いをどのように取るのかが今後の課題である。

表−2 の地点のうち，那覇港のデータに対しては形状

Hosking・Wallis (1997) が提唱している地域頻度解析では，最適分布関数が定まった後でモンテカルロ法の試行によって信頼区間幅を推定するようになっている。彼ら

母数が負値であるので上限値はない。また，酒田港と中城湾についてはワイブル分布が最適合，むつ小川原港と那覇港については GEV 分布が最適合である。その他の11

量が上限値で頭打ちとなる危険性が例証された。

の方式では，特定の地域内の複数の標本に適合する分布関数の候補として 4 母数あるいは 5 母数型含めているの

謝 辞で，モンテカルロ法以外に信頼区間幅を推定することは困難であろう。しかし，地域共通の分布関数の候補を3本論文で使用した高波のデータは港湾空港技術研究所

母数型に限定するならば，あらかじめモンテカルロ法をで所管している NOWPHAS のデータベースに基づくも駆使して信頼区間幅を推定するための経験式を導くことのであり，データ使用を承認してくださった河合弘泰

も可能であろう。

海象情報研究領域長のご厚意に深謝する次第である。

いずれにしても，まずは高波のデータを気象原因別に分類することから作業を進める必要があろう。

参考文献

１０．まとめ北野利一・間瀬 肇・喜岡 渉・矢野陽一郎 (2002)：一般化
パレート分布による極値波浪解析−拡張形状母数の推定−，
海岸工学論文集，第49巻， pp. 161-165.

本論文は，高波の極値統計解析に L-moments 法を適

合田良実 (2008)：「耐波工学 – 港湾・海岸構造物の耐波設計–」，鹿島出版会，430p.

用するための一つの試論であり，今後検討を進めるべき課題が多く残されている。ここまでの段階で明らかに

合田良実・小長谷 修・永井紀彦 (1998)：極値波浪統計の母分布関数に関する実証的研究，海岸工学論文集，第45巻，pp.211-215.

なった事項を列挙すると以下のようになろう。
1)

L-moments 法でこれまで考慮されていなかったワイ

外山奈央子・水野 量 (2002)：L-moments を用いた地域頻度解
析による全国アメダス地点における確率降水量の推定，気象庁研究時報，54巻 5−6号合併号，pp. 55-100.

ブル分布を当てはめ候補の極値分布関数にふくめることができた。
2)

標本から計算される L-skewness に基づいてワイブ

Coles, S. (2001): An Introduction to Statistical Modeling ofExtreme Values, Springer, 208p.

ル分布，一般化極値 (GEV) 分布，および一般化パGreenwood, J A., Landwehr, J. M., Matalas, N. C., and J. R.
Wallis (1978): Probability weighted moments: Definition andrelation to parameters of several distributions expressable ininverse form, Water Resources Res., Vol. 15, No. 5, pp. 10491064.

レート (GPA) 分布の形状・尺度・位置母数の 3 母数を求める方法が明確化された。
3)

標本に対する極値分布関数の適合度判定の指標として，順序統計量の上位 20 個の 2 乗平均偏差を用い

Hosking, J. R. M. (1990): L-moments: Analysis and estimation of
distributions using linear combinations of order statistics, J.
Roy. Statistical Soc., Series B, 52, pp. 105-24.

る TUD (Twenty-Up Deviation) を新たに考案した。
4)

日本沿岸で長期間にわたる波浪観測データが利用できる 8 地点について L-moments 法による極値統計解

Hosking, J. R. M. and J. R. Wallis (1997): “Regional Frequency
Analysis,” Cambridge Univ. Press, 224p.

析を試みたところ，波高の閾値レベルによる TUD

van Gelder, P. H. A. J. M. (2000): Statistical Methods for the RiskBased Design of Civil Structures, Ph.D Thesis Delft
University of Technology, 249p.

変化から高波の母集団が単一あるいは複数混在を判別する可能性が示された。

8 地点中でワイブル分布を最適とするものが 2 地点，
van Gelder, P. H. A. J. M, J. De Ronde, and N. W. Neykov (2000):
GPA 分布を最適とするものが 4 地点，GEV 分布を
Regional frequency analysis of extreme wave heights: trading
space for time, Coastal Engineering 2000 (Proc. 26th ICCE,最適とするものが 2 地点であった。
Sydney), ASCE, pp. 1099-1112.
6) GPA 分布は適合する地点が多いものの，確率統計

5)12