特集「時空間統計解析：新たなる分野横断的展開」［原著論文］

統計数理（2012）
第 60 巻 第 1 号 73–91
c 2012 統計数理研究所


地域気象観測に基づく波浪の予測手法に関する検討
1

甫喜本 司 ・清水 邦夫

2（受付 2011 年 6 月 29 日；改訂 11 月 14 日；採択 11 月 17 日）

要旨海上で発生する波浪の変化に関する予測方法の開発は，海象のメカニズムを理解することはもとより，海洋に関わる様々な人間活動にとって重要な意義をもっている．しかしながら，海象の変化の予測は現象自体がもつ複雑さに加えて観測データの十分な取得が困難である場合が多く，こうした背景がこの予測問題を困難なものとしてきた．本稿では，この予測問題に対する一つの考え方として，広域の地域気象観測より得られる情報を基にして海上の波浪の変化を予測するための方法について統計的観点より検討を行う．具体的には，気象庁が近年日本全国の広い地域にわたって設置してきた地域気象観測システムであるアメダス（AMeDAS），及び全
国の沿岸域に設置されてきた波浪計からの観測情報を基にして，両者の時間的・空間的な関係を表現する一つの統計モデルを開発する．開発されたモデルの予測における有効性については
予測実験を行って統計的に評価を行った．その結果，広域の観測情報を考慮することによって一つの気象観測局の観測情報のみから予測を行う場合よりも予測精度上の効果があるという可能性が示された．また，このモデルは気圧配置の変化する一年間を通しても概ね波高の予測性能に改善を与える可能性がある点が示された．
キーワード： 海象データ，気象データ，予測，von Mises 過程，時空間モデル．
1.

緒言

海洋波は海洋の災害を引き起こす大きな原因の一つとなっており，その予測は単に現象自体に関する理解にとどまらず，海洋に関わる様々な人間活動にとっても重要となっている．Sverdrup
and Munk（1947）
や Pierson et al.（1960）
をはじめとした論文の中で海洋波の波浪推算法の基礎理論が提案されてから半世紀以上が経った現在，波浪予測に関する技術は近年における現象の計測技術と計算機技術の著しい進歩にも支えられながら大きな進歩を遂げてきた．
海上で発生する波浪の変化は上空で発生している海上風の影響を大きく受けるため，その予測を行うためには波浪自体の観測に加えて海上で発生する気象要因の変化を観測することが望まれる．しかし，伝統的な観測ブイや観測船を用いた計測には設備の問題や機器の故障といったリスクが常に存在する上に，観測を継続的に行うことが困難な場合が多い．気象や海象に関わる物理的要因は，時間的・空間的に相互連関をもちながら変化していることが予想されるが，このことを利用して，陸上で常時観測可能な気象要因の変化に基づいて波浪の変化を予測する
1
2

東京大学大学院 数理科学研究科：〒153–8914 東京都目黒区駒場 3–8–1
慶應義塾大学 理工学部：〒223–8522 神奈川県横浜市港北区日吉 3–14–1

74

統計数理

第 60 巻

第 1 号 2012

ことがある程度可能となるだろうか．日本では，気象庁が常時観測型の地域気象観測システムであるアメダス（Automated Meteorological Data Acquisition System, AMeDAS）の観測情報を提供する地域気象観測局を全国各地に約 1300 箇所設置してきた．また，波浪の観測点についても気象庁や国土交通省港湾局によって波浪を観測する波浪計がおよそ 80 箇所の沿岸域に設置されてきた．本稿では，広域で観測されたアメダスの気象観測記録を活用した波浪の予測問
題について，統計的手法の開発を通して検討を行う．
波浪の変化を検討するための伝統的な考え方の一つは，海表面の変化を確率的な現象と捉えることである．風によって波浪が発達する際の海面変動が従うと考えられる確率過程は風速の状況に応じて異なる．風速が十分弱い下での海面変動は正規過程で近似できるが，風波の発達段階以降における海面変動は正規過程とみなすことができない局面となることが実験的にも知られている．こうした波浪の変化に関する特性を把握するため．海面の上下方向の振動に関する振幅にあたる有義波高（signiﬁcant wave height）の変化に関する分析が行われるようになった．
短期の有義波高の変化は定常性をもつ時系列とみなすことができ，ARMA モデルのクラスで表現することができる（e.g., Spanos, 1983）．一方，波浪の発達現象を含む長期的な変動については定常性が崩れるため，非定常性の観点（e.g., Scheﬀner and Borgman, 1992; Athanassoulis and，あるいは非線形性の観点（e.g., Scotto andStefanakos, 1995; Guedes Soares and Ferreira, 1996）
Guedes Soares, 2000）から統計的なモデル化の検討が行われてきた．また，第 4 章で述べられるように風速や風向といった風の変化に関する観測値に基づくモデル化についても，近年の自然エネルギーの研究の中で積極的に検討されるようになった．しかし，風の変化が上記の有義
波高の変化に与える影響を評価するためのモデルを開発する場合，問題が発生する．その一つ
は風向に関する角度時系列を含んだ多変量時系列データを統計的に取扱うための方法に関する検討が必要となる点である．そしてもう一つの問題は，従来検討されてきた波浪に関するモデ
ルが海上の観測ブイや観測船からの海洋観測で得られた観測データに基づくものが一般的であるため，気象の空間的な変化の特徴を考慮するための方法の開発が必要となる点である．本稿
ではこうした問題に対する一つのモデル化について検討すると共に，広域で観測されるアメダ
スの気象観測情報を利用して，沿岸域で発生する波浪の変化を予測するための方法を開発することを目標とする．また，開発されたモデルに基づく波浪予測の有効性について数値実験より検証する．
本稿の構成は以下の通りである．次節では，本稿で取り扱う陸上の気象観測と波浪観測に関する実際面について述べる．第 3 節ではモデル化を行うための統計的な予備調査として，観測データが有する相関構造に関する検討を行う．第 4 節では上記の結果を手掛かりとして波高変化を予測するためのモデル化について検討する．第 5 節で開発されたモデル構造の有効性について予測実験を通して評価した後，季節毎におけるモデルの実用性を第 6 節で検証する．これらの実験を通したモデル構造に関する評価と問題点は第 7 節でまとめられる．
2.

地上気象と波浪に関する観測の実際

本稿では波浪予測について，北海道松前郡松前町におけるケーススタディを取り扱う．松前町は漁業を主要な産業の一つとし，海に関わる活動に携わる人口が少なくない．その一方で，付近の海域（松前沖）は短時間で時化の状態になりやすく，天候の急変に伴い海難事故が発生しやすいことでも知られている．松前沖の沿岸域には気象庁によって超音波式沿岸波浪計が設置され，1979 年より波浪の常時観測が続けられてきた．波浪計は松前沖の沖合（42◦ 24 38 N，
140◦ 05 50 E），水深約 50m の海底に観測センサーである超音波送受波器が設置されており，超音波を用いて海面の水位変動を計測した後に，監視局である函館海洋気象台へデータが送られ

地域気象観測に基づく波浪の予測手法に関する検討

図 1.

75

波浪計と地域気象観測局の位置．

る．一方，観測海域の周辺にある地域気象観測局は最寄りとなる松前町をはじめとして江差町，奥尻町など周辺の諸地域に設置されている．各観測局では地上における気温，降水量，日照時
間，風向，風速をはじめとした多くの気象要因が観測されている．
図 1 は波浪計（wave recorder）の設置されている位置（松前沖）とその周辺地域に設置されてい，奥尻町（Okushiri），江差町（Esashi）, 森町（Mori），る地域気象観測局の中から松前町（Matsumae）
函館市（Hakodate），及び大間町（Ohma）の 6 地域の観測局の位置を示したものである．本稿の分析ではこの 6 地域の観測記録を使用する．図 2 は 2010 年 4 月 1 日から 5 月 31 日までの 2 か
月間に波浪計より得られた 1/3 有義波高（m）の変化，及び同時点において松前町，江差町，大間町の気象観測局で得られた風速
（m/s）
と風向
（rad．
）
の毎正時における観測記録を示している．
風速と風向は毎正時の直近 10 分間の観測データより平均値を求めたものであることに注意す
る．1/3 有義波高は毎正時の 25 分前から 5 分前までの 20 分間に収集された海面の水位変動の観測値よりゼロアップクロス法に基づいて海面振幅（波高）の系列を算出した後，得られた海面
振幅の最大値から上位 3 分の 1 の波高に関する平均値をとったものであり，目視観測による波高と近いという実際面からの理由より海洋学をはじめとする分野で広く用いられている．風向の観測値については，アメダスから出力された 16 方位の観測データを仰角に変換してラジアン単位で表示したものであり，真東を原点として時計回りを正にとっていることに注意する．
図 2 は春季の観測記録にあたる．この季節は気圧配置の移行期にあたるために不安定になりやすく，一年の中で風速と風向が非常に変化しやすい時期となる．
松前沖の波浪計の設置場所と松前町の地域気象観測局とは近い場所にあるため，この観測局
における風の観測値を海上風の観測値の近似値と考えることは不自然ではない．しかし，この観測局における風の変化に基づいて波高を予測すれば十分かという点については疑問が生じる．その背景の一つとして風向の変化の特徴が挙げられる．図 2 で観察されるように，風向が一定の方向を示している時間はそれほど長くはない．このことから，海上風を推測するためには様々な方位から風が吹く可能性をモデルに考慮することが効果的であると考えることができる．風速の変化についても観測点に最も近い松前町の風速の変化が波高の変化と必ずしも同期しているとは限らないと予想される局面があることが観察される．例えば，観測開始後およそ700 時間後から 1500 時間後までの期間において発生した波高の変化に注目すると，最寄りにあ

統計数理

76

図 2.

第 60 巻

第 1 号 2012

有義波高（松前沖）と風速，及び風向に関する観測記録の例（上から順に松前町，江差町，大間町）．

る松前町で観測された風速の変化よりも，やや遠方にある江差町で観測されたそれとの方がむしろ同期しているようにも観察される．この現象の物理的背景として地形的な特徴が風の挙動に及ぼす影響があると予想される．松前町は北海道の南西部，渡島半島の海岸線に位置しているが，半島は山の多い陸地が広がっている．また，その他の方位には津軽海峡，日本海といった開かれた海に広く取り囲まれている．このため，海側から吹く風に対しては風波の生成過程となり風送時間と風送距離に依存しながら風波が発達する．しかし，例えば北東の方位から吹く風に対しては山に風の流れが乱される可能性がある．従って，広域的な観測で得られた風速や風向の観測情報を考慮することによって効果的な波高の予測が得られるのではないかということが予想される．
3.

観測記録に基づく相関構造の予備分析

この節では，観測データを基に統計的なモデル化を行うための予備分析として相関構造の特徴について調査を行う．以下では WH t を時点 t における 1/3 有義波高とし，WS t ，WD t をそれぞれ地上観測局で観察された風速，及び風向を表すものとする．最初に，風向の角度時系列
データと風速や波高に関する線形時系列データが同時に観測された際の取り扱いについて検討

地域気象観測に基づく波浪の予測手法に関する検討

図 3.

77

{∇WS t }，{∇ cos(WD t )}，{∇(WS t cos(WD t ))} の変化と相互相関係数（（a）–（c））．

する．図 2 で示された海象に関する物理要因の変化は非定常性の強い時系列とみなすことができるため，以下では Box and Jenkins（1976, Chapter 4）で提案されている取扱いに従って，各変量について 1 階の階差変動{∇WH t },∇WH t = WH t − WH t−1などに着目する．図 3 は {∇WS t }，{∇ cos(WD t )} と {∇WH t } に関するそれぞれの変動，及び（a）{∇WS t } と {∇WH t }（，b）{∇ cos(WD t )} と {∇WH t }，及び（c）{∇(WS t cos(WD t ))} と {∇WH t }の 3 つのケースに関して，相互相関係数を
N−τ
1
t=1 (∇WH t − ∇W H)(∇WS t+τ − ∇W S)N ρ̂(τ )=  , τ = 0, 1, . . . ,
N
N
1
1
2
2
(∇WH
−
∇WH
)
(∇WS
−
∇WS
)
t
t
t=1
t=1
N
N
∇W H =NN
1 
1 ∇WH t , ∇W S =
∇WS t
N t=1
N t=1

などによって推定した結果である（点線は Bartlett の上下限を示す）．∇(WS t cos(WD t )) を一つの変量として定義した場合， ∇WH t との相互相関は ∇WS t や ∇ cos(WD t ) を単独で取り扱った場合に比べて高くなる傾向を示している．
次に，上記の相互相関係数が空間的にどのような変化をするかという点について検討する．
図 4 は図 2 で挙げた松前町，江差町，大間町の各地域気象観測局で得られた観測値を基に 3

統計数理

78

図 4.

第 60 巻

第 1 号 2012{∇(WS t cos(WD t ))} と {∇WH t } との相互相関係数に関する時間的変化（上段から順に [150–250]， [300–400]，[450–500] の各期間で推定）．

つの期間，[150–250]，[300–400]，及び [450–550] における {∇WS t cos(WD t )} と {∇WH t } との相互相関係数の変化を調べたものである．江差町の地域観測局は松前町の北部に，そして大間町のそれは東部に位置する．最初の時期において，波高の変化は観測海域の最寄りに設置されている松前町で観測される風と最も相関が強くなっていたが，時間の経過と共に徐々に小さくなっている．一方，時間の経過と共に相互相関が高くなる地域が北部にある江差町，さらに大間町へ移行していく傾向が観察される．この結果は，広域的に観測された風の情報を考慮する
ことが波高の予測にとってより効果的となる可能性を示唆している．
4.

地上気象観測に基づく波高の予測方法

この節では前節の予備分析で得られた風速と風向の特徴を考慮に入れて，1/3 有義波高の変化を予測するための一つのモデル化を検討する．風速の変化については ARMA モデルに基づく方法（e.g., Philippopoulos and Deligiorgi, 2009）
や GARCH モデルに基づくモデル化（e.g., Tol,が検討されてきた．一方，風向については ARMA モデルに基づくモデル
1997; Liu et al., 2011）
化が試みられているが（e.g., Erdem and Shi, 2011），観測域における気象の変化や地形的特性に強く影響を受けるため，風速とは異なる特徴をもつ変化が観測されることが一般的である．このように線形的な変化と方向の変化を同時に観測したデータに対して統計モデルをあてはめる
場合，方向のサイン変換とコサイン変換を説明変数とした多変量自己回帰モデルが検討されて

地域気象観測に基づく波浪の予測手法に関する検討

79

きた（Hokimoto and Shimizu, 2008）．以下では，このようなモデルの構造に風向の方向時系列
としての特性，及び風速と風向の空間的な特性を考慮することによって，予測の精度を改善できる可能性があるかという点に関心がある．4.1 項では，風向の時間変化が円周上の値をとる確
率変数に基づく一つの確率過程に従うことを仮定した上で，この特性を考慮に入れた多変量時
系列モデルを検討する．また 4.2 項では，開発されたモデルを広域の観測局で観測された風の観測データに対しても適用するための方法について提案する．以下では，観測データ {WH t }，{WS t }，{WD t }（t = 1, . . . , T ）に基づいて WH T +l（l = 1, . . . , L）の変化を予測する局面を考える．
4.1 単独の地域観測局に基づくモデル化
地域観測局が一つの場合のモデル化から検討をはじめることにする．最初に風向の変化が従う確率過程を定義する．本稿では風向の変化 {WD t }（−π ≤ WD t ≤ π ）が Breckling（1989）によって検討されている 1 次の von Mises 過程に従うと仮定する．いま，1 時点前の観測値に基づく WD t の条件付分布が以下を集中度ベクトルとしたある von Mises 分布に従うものと仮定
する．


 


cos(µ(WD),t )cos(WD t−1 )1ν(WD),t, k0 > 0, k1 > 0
= k1
+ k0
0
sin(µ(WD),t )sin(WD t−1 )ここで，µ(WD ),t は von Mises 分布の平均風向（−π ≤ µ(WD ),t ≤ π ），ν(WD),t はこの分布のばらつ
きを表す集中度（ν(WD),t > 0）を意味する．また，k1 ，k0 は正値をとる未知係数で，集中度の平均風向に関するサイン成分とコサイン成分を過去の風向によって説明できる成分とそれ以外の成分に分解した際の影響の度合をそれぞれ意味している．この下で， WD t の条件付確率密度関数は以下の形で書くことができる．
1
exp{k1 cos(WD t − WD t−1 ) + k0 cos WD t }(4.1)
f (WD t |WD t−1 ) =2πI0 (ν(WD),t )ここで，関数 I0 (·) は 0 次の第 1 種変形ベッセル関数を表す．
（4.1）は平均風向が
k1 sin(WD t−1 )µ(WD ),t = tan−1(4.2)
k1 cos(WD t−1 ) + k0
かつ集中度が(4.3)

ν(WD),t =(k1 cos(WD t−1 ) + k0 )2 + (k1 sin(WD t−1 ))2である von Mises 分布の確率密度関数
f (WD t |WD t−1 ) =1exp(ν(WD),t cos(WD t − µ(WD),t ))2πI0 (ν(WD),t )で書き換えることができる．これは von Mises 分布の位置パラメータである µ(WD),t と集中度
ν(WD),t が共に過去の観測値である WD t−1 と（k0 , k1 ）に依存しながら変化することを意味している．k0 と k1 が共に十分小さい場合には ν(WD),t も十分小さくなるため，この分布は一様分布で近似できる．一方，k0 ，k1 が大きくなるにつれて，ν(WD),t の値が大きくなるため µ(WD),tの周辺に集中した形状の分布に変化する．
の各変化に対する f (WD t |WD t−1 ) の変化を示
上記の傾向をやや具体的にみるため，（k0 , k1 ）した例を図 5 に示す．以下の計算は全て WD t−1 = 1.0 として行った．最初の例は k0 = 0.1 の下で k1 が（0.05, 0.3, 1.0, 3.0）の 4 つのケースにおける f (WD t |WD t−1 ) の変化を，また 2 番目の例では k1 = 0.1 とした下で k0 が（0.05, 0.3, 1.0, 3.0）の 4 つのケースにおける変化が示されている．ここで，それぞれのケースが実線，点線，破線，及び太い破線に対応していることに注意する．k0 が小さい場合，f (WD t |WD t−1 ) は k1 の値が大きくなるほど WD t−1（=1.0）の周辺

統計数理

80

第 60 巻

第 1 号 2012

図 5. （k0 , k1 ）の変化に対する f (WD t |WD t−1 ) の変化．

に集中した分布へ変化する傾向をもつ．これに対して，k1 の影響が小さい場合には k0 が大きくなると共に原点の周辺に集中した分布に変化する傾向がある．一方，3 番目と 4 番目の例は
k0 と k1 がそれぞれ 0.8 の場合における f (WD t |WD t−1 ) の変化について同様に調べたものである．平均については上記と同じ変化の傾向があること，及び平均の周辺に集中する度合が上記に比べて増すことが観察される．
上記の時間的な変化の過程を推定するため，以下では ki = exp(ci )（i = 0, 1）と表し（c0 , c1 ）をパラメータとしたモデルを検討する．いま，WD t に関する条件付密度関数が任意の時点 t で
f (WD t |WD t−1 , . . . , WD 1 ) = f (WD t |WD t−1 )を満たすものと仮定する．このとき風向の観測値に基づく尤度関数 f (WD 1 , . . . , WD T ) は，初期時点における確率密度関数を f (WD 1 ) とするとき
T1
exp exp(c1 ) cos(WD j − WD j−1 ) + exp(c0 ) cos(WD j ) · f (WD 1 )
2πI(ν0 (WD),j )j=2と書けることから，未知パラメータ（c0 , c1 ）に関する最尤推定を行うことができる．これらの推定量をそれぞれ ĉ0 , ĉ1 とするとき，（4.2）と（4.3）に関する推定量をそれぞれ以下で構成する．

exp(ĉ1 ) sin(WD t−1 )(4.4)µ̂(WD),t = tan−1
exp(ĉ1 ) cos(WD t−1 ) + exp(ĉ0 )(4.5)

ν̂(WD),t =

(exp(ĉ1 ) cos(WD t−1 ) + exp(ĉ0 ))2 + (exp(ĉ1 ) sin(WD t−1 ))2

地域気象観測に基づく波浪の予測手法に関する検討

81

図 6 は上記で得られた {µ(WD ),t } の推定値の挙動とその自己相関係数，及び {ν(WD),t } の推定値に関する変化をそれぞれ示した例である（ν(WD),t については底を 2 とした対数で表示していることに注意する）．前者の変化は時間の経過と共に比較的緩やかな変化を示す傾向がある
一方で，後者のそれについては平均や分散が時間と共に変化する傾向が観察される．このことから，{µ(WD ),t } については定常性をもつ変動，{ν(WD),t } については非定常な変動とみなすことが可能である．
次に，風向が上記の過程に従っている場合に風速や有義波高の変化へ与える影響を表現するためのモデル化を行う．以下では，Hokimoto and Shimizu（2008）
で検討されたモデルの構造を参考にして，上記で検討された時点 t における {µ(WD),t } と {ν(WD),t } の変化が波高，風速，及び風向の成分と同期する過程をモデル化する．具体的には，∇WH t に対して以下のように記述する．
∇WH t =p(1)

αi ∇WH t−i +

i=1

+

p
K



(1)

βi,k ∇(ν(WD),t−i WS t−i cos(kWD t−i ))i=1 k=1

p
K



(1)

γi,k ∇(ν(WD),t−i WS t−i sin(kWD t−i )) +i=1 k=1

+

p


(1)

p

(1)
δi cos(µ(WD),t−i )
i=1

(1)

(1)

2
ωi sin(µ(WD),t−i ) + εt , εt ∼ WN (0, σ1,h
)

i=1

図 6.

{µ(WD ),t } の推定値の変化と自己相関係数，及び {ν(WD),t } の推定値の変化．

統計数理

82

第 60 巻

第 1 号 20122
ここで p，K は次数，（α, β, γ, δ, ω ）は未知係数，そして ε(1)は平均 0，分散 σ1,hのホワイトノイtズ過程に従う確率変数とする．同様に ∇(ν(WD),t WS t cos(hWD t )) と ∇(ν(WD),t WS t sin(hWD t ))（h = 1, . . . , K ）に対して∇(ν(WD),t WS t cos(hWD t ))
=

p


(2)

αi ∇WH t−i +

i=1

+

p
K



(2)

βi,k ∇(ν(WD),t−i WS t−i cos(kWD t−i ))i=1 k=1

p
K



(2)

γi,k ∇(ν(WD),t−i WS t−i sin(kWD t−i )) +i=1 k=1

+

p


p

(2)
δi cos(µ(WD),t−i )
i=1

(2)

(2)

(2)

2
ωi sin(µ(WD),t−i ) + εt , εt ∼ WN (0, σ2,h)i=1の形でモデル化する．さらに，cos(µ(WD),t ) と sin(µ(WD ),t ) の変動についても
cos(µ(WD),t ) =

p


(3)

αi ∇WH t−i +

i=1

p
K



(3)

βi,k ∇(ν(WD),t−i WS t−i cos(kWD t−i ))i=1 k=1

p
p
K



(3)
(3)
γi,k ∇(ν(WD),t WS t−i sin(kWD t−i )) +δi cos(µ(WD ),t−i )+i=1 k=1

+

p


i=1

(3)

(3)

(3)

2
ωi sin(µ(WD),t−i ) + εt−i , εt ∼ WN (0, σ3,h)i=1の形でそれぞれ記述する．いま，時点 t における状態を（2K + 3）次のベクトル(K)yt≡ (∇WH t , ∇(ν(WD),t WS t cos (WD t )), ∇(ν(WD),t WS t sin (WD t )), . . . ,∇(ν(WD),t WS t cos (KWD t )), ∇(ν(WD),t WS t sin (KWD t )), cos(µ(WD),t ), sin(µ(WD),t ))で定義すると上記の表現が
(K)

yt

(4.6)

(K)

= A1

(K)

(K)

(K)

y t−1 + · · · + A(K)p y t−p + δ t ,

(K)

δt

∼ WN (0, Σ(K) )と多変量自己回帰モデルで表される．ここで，A(K)は未知の係数行列である．
i （i = 1, . . . , p）
（4.6）
を観測データに基づいて同定する際には以下のような方法で行う．まず {WD t } に基づいて，先に述べた最尤法により c0 ，c1 を推定した後，{µ̂(WD ),t } と {ν̂(WD ),t }（t = 1, . . . , T ）を（4.4），（4.5）を用いてそれぞれ推定する．そして，(K)yt≡ (∇WH t , ∇(ν̂(WD),t WS t cos (WD t )), ∇(ν̂(WD),t WS t sin (WD t )), . . . ,∇(ν̂(WD),t WS t cos (KWD t )), ∇(ν̂(WD),t WS t sin (KWD t )), cos(µ̂(WD),t ), sin(µ̂(WD),t ))によって y(K)を構成した後，これに（4.6）
をあてはめる．
（4.6）
は多変量自己回帰モデルであり，t未知の係数ベクトル A(K)は最小 2 乗法で推定することができる（例えば，Brockwelli （i = 1, . . . , p）
．時点 T から l 期先の線形予測量は次のように構成できる．
and Davis, 1996）
(4.7)

(K)

(K) (K)(K) (K)(K)z T +l−1 + Â2 z T +l−2 + · · · + Â(K)p z T +l−p ,
(K)
(K)
(K)
(K)
z T +l−m = ŷ T +l−p (l > p)z T +l−m = yT +l−p (l ≤ p),ŷ T +l = Â1

ここで Â(K)は A(K)の最小 2 乗推定量である．こうして WHT からの l 期先予測値，WHT +lii（l = 1, . . . , L）は ŷ (K)の予測量を基に得ることができる．
T +l

地域気象観測に基づく波浪の予測手法に関する検討

83

4.2 広域の気象観測情報を考慮するためのモデル化
次に，前項で検討されたモデルを複数の地域観測局で得られた観測データに対して適用する
ための方法について考える．以下では，実際面の観点も考慮し K = 1 として検討を続ける．第2 節で行った観察より，風向の変化は時間の経過と共に速い速度で変化するという構造をモデ
ルに反映する方が自然と考えられる．そこで，時間によって異なる地域観測局の方向から海上(K)に風が吹き付ける構造を考慮に入れたモデルを考える．提案するモデルは，（4.6）
において y tを以下で定義したものとする．
(s∗ )yt|s∗ ≡ (∇WH t , ∇(νt

(s∗ )

WS t

(s∗ )

cos (WD t(s∗ ))), ∇(νt

(s∗ )

WS t

(s∗ )

sin (WD t

)),


(s∗ )
(s∗ )
cos(µ(W D),t ), sin(µ(W D),t ))∗∗∗ここで，WS t(s ) ，WD t(s ) ，νt(s ) はそれぞれ，S 箇所の地域観測局 s（s = 1, . . . , S ）の中で最も波浪の発達に与える影響が大きいと考えられる観測局 s∗ において観測された WS t ，WD t ，νtの値とする．
s∗ の値は観測値に基づいて統計的に推定する．具体的には，地域観測局 s の中で 1 期先予測における平均 2 乗予測誤差が最も小さくなるときの s を与えることにする．各観測局 s におい(s)て得られた風速と風向の観測データをそれぞれ {WS (s)t }，{WD t } とするとき，
(s)

yt

(s)

(s)

(s)

(s)

(s)

(s)

= (∇WH t , ∇(νt WS t cos (WD t )), ∇(νt WS t sin (WD t )),
(s)

(s)

cos(µ(WD),t ), sin(µ(WD ),t ))

の変化に（4.6）をあてはめることにより WH t の 1 期先予測値を得る．各時点 t において得られ(s)
た波高の 1 期先予測値を WH t （t = 1, . . . , T ）とするとき，局所的な平均予測誤差平方和
W (s, τ (s)) =

1
τ (s)

T


(s)

(WH t − WH t )2

t=T −τ (s)+1∗が最も小さくなるときの s を s として選択する．ここで τ (s) は時間平均をとる際の時間幅であるが，WS (s)や WD (s)tt の変化は観測局毎に異なる非定常性を示す可能性があることから s に依存することに注意する．こうして選択された s∗ に基づいて（4.7）から ŷ T +l|s∗ を得ることによって WH T +l の予測値が求められる．
5.

モデルの予測性能に関する検討

以下では，前節で検討されたモデルがどの程度有効な予測を与えるかという点について統計的観点より評価を行う．
5.1 予測実験
予測実験では，観測データに基づいて外挿による 1/3 有義波高の長期予測の値を計算した
後，実測値との予測誤差について統計的に評価する．具体的な手順は以下の通りである．はじめに，サンプル数 100（100 時間に相当）の観測データにモデルをあてはめて未知パラメータを推定し，予測値を 5 期先
（5 時間先に相当）まで求める．観測データの範囲で任意に予測の開始時点を選んでこの計算を繰り返す．次に，得られた予測値と実測値に基づいて平均 2 乗予測
誤差（MSE），及び相関係数（COR）を評価する．また，伝統的に用いられる時系列モデルを導入して同様の評価を行い比較する．モデルをあてはめる際の次数については，赤池情報量規準（AIC）を用いて選択した．

統計数理

84

第 60 巻

第 1 号 2012

5.2 直近の地域観測局における観測記録に基づく予測性能
最初に，地域気象観測局が一つのケースからみていくことにする．第 2 節で示した春季における波高の変化について，観測海域の最寄りに位置する松前町の観測局で得られた風速と風向の観測情報に基づいて予測する実験を行った．この実験で比較対象としたモデルは以下の通りである．
（i） WH t = pi=1 αi WH t−i + ε1,t , ε1,t ∼ W N (0, σ12 )（ii） ∇WH t = pi=1 βi ∇WH t−i + ε2,t , ε2,t ∼ W N (0, σ22 )（iii） y t = A1 y t−1 + · · · + Ap y t−p + δ t , δ t ∼ WN (0, Σ)，
y t = (∇WH t , ∇WS t )（iv） y t = A1 y t−1 + · · · + Ap y t−p + δ t , δ t ∼ WN (0, Σ)，
y t = (∇WH t , ∇(WS t cos(WD t )))（v） y t = A1 y t−1 + · · · + Ap y t−p + δ t , δ t ∼ WN (0, Σ)，
y t = (∇WH t , ∇(νt WS t sin(WD t )), sin(µ(WD),t ))（vi） y t = A1 y t−1 + · · · + Ap y t−p + δ t , δ t ∼ WN (0, Σ)，
y t = (∇WH t , ∇(νt WS t cos(WD t )), cos(µ(WD),t ))

（i）
と
（ii）
は波高 {WH t } の観測値のみに基づいて予測を行うもので，前者は 1 変量の定常 AR
モデル，
（ii）
は非定常時系列の予測を念頭においた 1 変量の ARIMA（p, 1, 0）
モデルである．
（iii）
は風速の変化のみを共変量として考慮したモデル，（iv）から（vi）は風速と風向の変化を共変量として考慮したモデルである（第 4 節で検討されたモデルは（v）と（vi）にあたる）．
表 1 は上記の各モデルを用いて外挿による予測を行った際の精度の比較を行った結果であり，（i）から（vi）を用いて 5 期先までを予測した際の MSE，及び COR の値を示している．
（i）と
（ii）
の比較より，1 階の階差変動に着目した ARIMA（p, 1, 0）
モデルの構造には定常 AR モデルを用いた場合に比べて予測精度の効果があることが認められ，非定常なモデルが有効となることが確認できる．次に，（ii）と（iii）及び（ii）と（iv）とをそれぞれ比較すると，いずれも後者の方が予測ステップ数が増す毎に予測性能が良くなることが認められる．すなわち，風の変化を共変量として考慮することによって，波高の変化のみに基づく予測よりもより有効な予測を行う可能性があることが認められる．さらに
（iv）
と
（v）
，
（iv）
と
（vi）
をそれぞれ比較すると，von Mises
過程に基づく風向の変化の構造を考慮することによって予測精度がさらに改善される可能性が示唆されている．風向の変化に仮定された von Mises 分布のパラメータの変化が観測データの変化と同期しやすく，この性質が予測精度の改善に結びついたものと考えることができる．
5.3 広域の気象観測情報を考慮した場合の効果
次に複数の地域気象観測局を考慮した場合の予測効率について検討を行う．表 2 は第 2 節で
表 1.

直近の地域気象観測情報に基づく予測性能の比較（実験回数：130）．

地域気象観測に基づく波浪の予測手法に関する検討
表 2.

85

広域の地域気象観測情報に基づく予測性能の比較（実験回数：130）．

とりあげた 6 つの観測局における風速と風向の観測データを基に， 5.2 項と同様の予測実験を行った場合の MSE と COR を示している．ここで比較対象としたモデルは以下の通りである．
vii） y t = A1 y t−1 + · · · + Ap y t−p + δt , δ t ∼ WN (0, Σ)，

(1)
(1)
(6)
(6)
y t = (∇WH t , ∇(WS t cos(WD t )), . . . , ∇(WS t cos(WD t )))viii） y t = A1 y t−1 + · · · + Ap y t−p + δt , δ t ∼ WN (0, Σ)，

(s∗ )
(s∗ )
y t = (∇WH t , ∇(WS t cos(WD t )))ix） y t = A1 y t−1 + · · · + Ap y t−p + δt , δ t ∼ WN (0, Σ)，

(s∗ )
(s∗ )
(s∗ )
(s∗ )
y t = (∇WH t , ∇(νt WS t sin(WD t )), sin(µ(W D),t ))x） y t = A1 y t−1 + · · · + Ap y t−p + δt , δ t ∼ WN (0, Σ)，

(s∗ )
(s∗ )
(s∗ )
(s∗ )
y t = (∇WH t , ∇(νt WS t cos(WD t )), cos(µ(W D),t ))（vii）は最も基本的なモデルとして，6 つの地域観測局における風の観測情報を全て同時点における状態ベクトルとして考慮したものである．表 2 の結果は，（vii）を用いた場合には 5.2 項で用いた 1 地域を基にしたモデルよりも予測精度が低下してしまう傾向があることを示してい
る．
（iii）
から
（vi）
における状態ベクトルの次元は最大 3 であったのに対して（vii）
では 7 となるため，（4.6）
をあてはめる際に推定すべき係数行列の要素は後者の方が多くなる．このことはモデ
ルのあてはめに伴う推定誤差は（vii）の方が大きくなり，結果として予測誤差が大きくなってしまうことが考えられる．
（viii）は未知パラメータの数が（vii）に比べて少なく，前節で行った予測実験の結果と比較すると予測効率の改善に結びついていることが認められる．さらに（viii）の構造に角度分布の集中度の変化を考慮した
（ix）
と
（x）
は，
（viii）
による予測効率をさらに改善する可能性があることがわかる．
6.

季節別にみたモデルの有効性

前節の実験結果は春季の観測データに適用したケースであるが，季節毎に異なる気圧配置の特徴を有する日本の場合には，開発されたモデルが効果的な予測を与えるか否かという点についても季節毎に検討する必要がある．図 7 は松前沖と松前町において 2010 年の春季
（4 月 1日–5 月 31 日），夏季（7 月 1 日–8 月 31 日），秋季（10 月 1 日–11 月 31 日），及び冬季（1 月 1 日–2
月 28 日）に観測された 1/3 有義波高，風速，及び風向の観測記録を示したものである．一見して，季節毎にこれらの時間的変化の特徴は少し異なっているように観察される．特に冬季は北西の季節風が比較的安定して吹く傾向があるため，他の季節に比べて波高，風向に関する時間的変化の特徴が少し異なっているように観察される．
6.1 予測性能の特徴
予測実験は，松前町の観測局のみの情報に基づくモデルとして 5.2 項で用いた
（i）から
（iv）ま
での各モデル，及び 6 地域の観測情報に基づくモデルとして（vii），（viii），（x）の各モデルを用い

統計数理

86

図 7.

第 60 巻

第 1 号 2012

季節毎の有義波高，風速，風向に関する観測記録（松前沖，松前町）．上から順に有義波
高（m），風速（m/s），風向（rad.）．

て前節と同様な予測精度の実験を行った．表 3 は季節別の MSE の結果を示している．モデル
は春季，夏季，秋季については他の予測性能を改善する構造をもっていると評価される．こ（x）のことは 6 地域の風の変化を考慮に入れることで 1 地域の観測情報に基づく予測の精度を改善できる可能性を示唆している．冬季に関しては他のモデルと同程度の予測結果にとどまり，より効果的な結果を得るには至らなかった．これは風向の変化の分布に関する集中度が他の季節に比べて高く，時間と共に風向変化の構造が変化するというモデルの構造が不利になった結果

地域気象観測に基づく波浪の予測手法に関する検討

表 3.

季節毎の MSE の比較（実験回数：130）．

87

統計数理

88

第 60 巻

第 1 号 2012と考えられる．なお，同じように 6 地域の風の変化を考慮に入れたモデル（vii）
は未知パラメー
タの数が増えることが原因となって予測効率の改善に十分な効果を与えることは難しいことに注意する．
6.2 波高の変化に影響を与える風向の推測
波高の変化に影響を与える風の流れが時間と共にどのように変化しているかという点を調査することは現象を知る上で重要となる．以下では本稿で検討されたモデルがこの点をどの程度説明できるかということについて検討してみよう．図 8 は 6 地域における風向のヒストグラム
について季節毎に示したものである．風向は真北を 1 として時計回りに 16 方位で示されていることに注意する．一方，図 9 はモデル（x）に基づいて観測値から推定された地域観測局 s∗ に関して季節毎の分布を示したものである．ここで，横軸にある MA，OK，ES，MO，HA，OH
はそれぞれ松前町，奥尻町，江差町，森町，函館市，及び大間町の地域気象観測局であることに注意する．
物理的には地域気象観測局の上空で観測される風向と海上の波浪との間には以下の傾向があることが予想される．春季にあらわれやすい東風の場合，大間町や函館市の周辺域の上空で観測される風の流れが松前町を経由して松前沖の観測点上空に達する可能性が高く，これらの地域の観測局と波浪の観測点の上空における風の変化には比較的高い相関があることが予想される．夏季にあらわれやすい南西の風や秋季にあらわれやすい西風の場合，西側が広く海と接し十分な風送距離を有する松前町や江差町，及び奥尻町で観測される風の流れが波浪の観測点上空にも達するものと予想される．冬季になると顕著になる西北西方向からの季節風は奥尻町，江差町，松前町で観測される風の流れが波浪の観測点上空にも達すると共に，大間町の上空付
近も同方向からの風の流れが通過するものと予想される．
一方，図 9 をみると春季の波高の変化を予測するためには函館市の上空における風の影響が他の地域に比べ高くなっており，次いで松前町，江差町の上空で発生する風の影響が高くなっている．これは，図 8 より春季においては東風が最も高い頻度で発生し，次いで南西から西より吹く風が高い頻度で発生していることが原因として考えられる．夏季になると，図 9 で推定された分布は松前町，奥尻町で観測される風の頻度が相対的に高くなる．これは夏型の気圧配置への移行に伴い，図 8 で示されるように風向の分布が南西風を中心としたものに変わっているからである．秋季に入ると松前町，江差町で観測される風の影響が相対的に高い頻度で推定される．この時期は冬型の気圧配置へ移行していく時期にあたり，冬季に頻度の高い西北西の風の頻度が高くなっていることがこのような原因を生む原因と考えられる．冬季になると大間町，松前町，奥尻町，江差町の頻度が相対的に高く推定される．この季節は全域が冬型の気圧配置に覆われ，西北西の季節風が高い頻度で発生していることが上記の 4 地域の頻度を相対的に高くした原因と考えられる．
7.

結論と今後の展望

本稿では，アメダスの地域観測局で得られた陸上における広域的な風の変化に基づいて海上で発生する波高の変化を予測する統計的な方法について検討を行った．海象の変化には対象とする海域付近における固有な気象や海洋の物理的環境，及び周辺の地形的な特徴などが複雑に影響するため，観測データに一般的な統計モデルをあてはめても意味のある結果が得にくい場合が多い．本稿では，広域的な風の変化と風向の変化が従うプロセスを考慮に入れながら時間的・空間的な変化に関する一つのモデル化を試みた．その結果，複数の地域観測局から得られた風の観測情報が効果的に反映されることによって，一般的な時系列モデルを用いた場合の波
高の予測性能に改善を与える可能性があることが示された．本稿では風向の変化が von Mises

地域気象観測に基づく波浪の予測手法に関する検討

図 8. 季節別の風向の観測値に関するヒストグラム．

図 9.

季節別に推定された s∗ の分布．

89

90

統計数理

第 60 巻

第 1 号 2012

過程に従うことを仮定してモデル化を行ったが，より有効となる確率過程が存在する可能性も考えられる．こうした検討をさらに進めることによって広域的な風の変化に対する理解や，これが海上の波浪の発達や減衰の変化へ及ぼす影響について，より精密に分析することが可能になるものと期待される．

参 考 文 献Athanassoulis, G. A. and Stefanakos, C. N.（1995）. A nonstationary stochastic model for long-term（C8）, 16149–16162.
time series of signiﬁcant wave height, Journal of Geophysical Research, 100
Box, G. E. P. and Jenkins, G. M.（1976）. Time Series Analysis, Forecasting and Control（revised
edition）, Holden-Day, San Francisco.
. The Analysis of Directional Time Series: Applications to Wind Speed andBreckling, J.（1989）
Direction, Lecture Notes in Statistics, 61, Springer-Verlag, Berlin.
Brockwell, P. J. and Davis, R. A.（1996）. Introduction to Time Series and Forecasting, SpringerVerlag, New York.
. ARMA based approaches for forecasting the tuple of wind speed andErdem, E. and Shi, J.（2011）（4）, 1405–1414.
direction, Applied Energy, 88
Guedes Soares, C. and Ferreira, A. M.（1996）. Representation of non-stationary time series of signiﬁcant wave height with autoregressive models, Probabilistic Engineering Mechanics, 11,
139–148.
Hokimoto, T. and Shimizu, K.（2008）. An angular-linear time series model for waveheight prediction,
Annals of the Institute of Statistical Mathematics, 60, 781–800.
. Comprehensive evaluation of ARMA-GARCH（-M）approaches
Liu, H., Erdem, E. and Shi, J.（2011）
for modeling the mean and volatility of wind speed, Applied Energy, 88, 724–732.
. Statistical simulation of wind speed in Athens, Greece
Philippopoulos, K. and Deligiorgi, D.（2009）
based on Weibull and ARMA models, International Journal of energy and environment, 3（4）,
151–158.
. Practical Methods for Observing and ForePierson, W. J., Neumann, G. and James, R. W.（1960）
casting Ocean Waves by Means of Wave Spectra and Statistics, Reprint edition, U.S. Navy
Hydrographic Oﬃce, Washington, D. C.
Scheﬀner, N. W. and Borgman, L. E.（1992）. Stochastic time-series representation of wave data,
Journal of Waterway, Port, Coastal and Ocean Engineering, 118（4）, 337–351.
. Modelling the long-term time series of signiﬁcant wave
Scotto, M. G. and Guedes Soares, C.（2000）
height with non-linear threshold models, Coastal Engineering, 40, 313–327.
Spanos, P. D.（1983）. ARMA algorithms for ocean wave modelling, Journal of Energy Resources
Technology, 105, 300–309.
. Wind, Sea and Swell : Theory of Relation for Forecasting,
Sverdrup, H. U. and Munk, W. H.（1947）
No. 601, U.S. Navy Hydrographic Oﬃce, Washington, D. C.
. Autoregressive conditional heteroscedasticity in daily wind speed measurements,
Tol, R. S. J.（1997）
Theoretical and Applied Climatology, 56, 113–122.

Proceedings of the Institute of Statistical Mathematics Vol. 60, No. 1, 73–91 (2012)91

Predicting Waveheight Based on Ground-based Monitoring of Wind
Tsukasa Hokimoto1 and Kunio Shimizu2
1 Graduate School of Mathematical Sciences, The University of Tokyo
2 Department of Mathematics, Keio University

The ocean wave is one of the physical factors that cause serious sea disasters, and its
prediction provides information for various human activities related to the sea. However,the prediction has been a diﬃcult problem until now, because it is usually diﬃcult tocarry out constant monitoring of changes in various meteorological factors relating to a
sea area of interest. To solve this problem, in the present article, we develop a statisticalpredictor for waveheight based on changes in meteorological factors obtained by groundbased monitoring.
In our country, the Japan Meteorological Agency has set up many regional stations
for ground-based meteorological monitoring from an Automated Meteorological Data Acquisition System (AMeDAS), and sensors for observing wave recorders in many coastal
areas. In this article, we use measured data on wind speed and wind direction from multiple AMeDAS stations and data on waveheight obtained from a wave recorder, measured
at Matsumae-oki, Hokkaido, Japan. Some preliminary statistical analyses suggest thatthe correlation structure obtained from the regional AMeDAS station, which gives the
maximum impact on the cross correlation between the wind and the waveheight, changes
over time. Therefore, we developed a space-time model for predicting the waveheight.
More precisely, assuming that the change of wind direction follows von Mises process, wedeveloped a nonstationary time series model based on wind speed, wind direction as wellas waveheight, which takes into account the possibility that the wind direction changes
within a short time.
In order to examine the applicability of the developed model, we carried out prediction experiments. The results suggested that the developed model can improve prediction
accuracy by using the predictors based on the observation at the nearest AMeDAS stationfrom the observation point of the waveheight. It also suggested that the above eﬀectiveness on the prediction tends to be robust throughout a year, although there is room for
improving its accuracy in winter. The developed model gave a time series structure fortaking into account the von Mises process. Further consideration on the directional process which is more ﬂexible for expressing the change of wind direction remains as a futurework.

Key words: Sea state data, meteorological data, prediction, von Mises process, space-time model.