この記事はBrainPad AdventCalendar 201710日目の記事です。
波浪の統計量をPythonで解析した事例が少なそうなので、この記事ではまず手始めに、波の観測データから波高と周期をPythonを使って求めてみました。
また、平均でも最大でもない、波浪解析で特有の有義という統計量を紹介します。
はじめに
実際の海面では、いろいろな波高と周期の波の集まりで、それらが不規則に重なり合い水面が変動しています。このような波は一般に不規則波と呼ばれます。
では、サンプルとして、ある港の実際の波浪観測地点の20分間の観測データを見てみましょう。
import csv
import matplotlib.pyplot as plt
time = []
water_level = []
# 時間(time)と水位(water_level)の2列のcsvを読込
with open('wave_sample_20m.csv', 'r') as f:
reader = csv.reader(f)
header = next(reader)
for x, y in reader:
time.append(float(x))
water_level.append(float(y))
plt.plot(time, water_level)
plt.title("wave sample data (20min)")
plt.xlabel("time (sec)")
plt.ylabel("water lavel (m)")
plt.show()
＜サンプルデータ（20分間）＞
サンプルデータは0.5秒間隔で20分間の1200個のデータをプロットした時系列グラフになります。不規則に変動してますね。では、この波の波高、周期はどのくらいでしょうか？後ほど説明しますが、波高は波の山から谷までの高さになります。
グラフを見た感じ、波高は1m弱？周期は10秒強？
波浪の解析方法
上記の不規則波においても、波高・周期は一義的に表されます。天気予報とかにある波高１ｍも不規則波を一義的に表しています。
では、どのように解析するか？不規則波の解析方法には大きく下記２通りの方法があります。
- 波別解析法：観測した水位変動から、波を個別に定義し、１つ１つの波の諸元（波高・周期）を統計的に処理する。
（例）ゼロアップクロス法 - スペクトル解析：時系列データや空間変動の周期性や類似性を定量的に解析する方法。不規則波をいろいろな周波数の成分波の重ね合わせとし、それぞれの周波数とエネルギーの関係を解析する。
（例）フーリエ解析
噛み砕いて言うと、１つ１つの波に分割し統計処理する方法が波別解析法。定常確率過程としてスペクトラム、確率密度関数で処理する方法がスペクトル解析となります。
前者の方が簡単そうですね。今回は手始めに波別解析法のゼロアップクロス法を用いて、上で示したサンプルデータの波高、周期を解析してみます。
ゼロアップクロス法
まず、１つ１つの波の波高と周期を下図のようにイメージします。
観測データから得られた時系列の波形の、平均値から波形が負(-)から正(+)へ横切った時点から、水面が上昇して下降し次に同じ状態になるまでを１つの波として扱います。ゼロからアップにクロスするということで、ゼロアップクロス法と言います。
＜ゼロアップクロス法での波高と周期＞
ゼロアップクロス法により、20分間の波浪データを分割すると、通常は100波以上の波形に分割されます。では、波浪データの波高は、これら100波以上の波を平均した波高でしょうか？それとも最大の波高でしょうか？
平均でもなく、最大でもない、有義って何？
波の波高・周期は、平均波でも最大波でもなく、有義波と言うもので表します。
有義波は分割された波形の波高の高い方から順に全体の1/3の波で選ばれ、これらの波の波高、周期を平均したものを有義波高、有儀周期と呼びます。（例えば全体で100波あれば、波高の上位33個の波を平均する。）
（気象庁HP 波浪の知識）
また、この有義波は、（漁師や船員など）の熟練の観測者が目視で観測する波高や周期に近いと言われています。
コード例
前置きが長くなりましたね。それでは、先ほどのサンプルデータからゼロアップクロス法で有義波を求めてみましょう。
from statistics import mean
ave = mean(water_level)
WL = [x - ave for x in water_level] #平均水位からの変動
x1 = 0
waves = []
for x in range(1,len(time)):
if WL[x-1] < 0 and WL[x] > 0: #0(平均)をクロスする時点
height = max(WL[x1:x]) - min(WL[x1:x]) #個々の波高
period = (x - x1) * 0.5 #個々の周期
waves.append([height, period])
x1 = x
waves.sort(key=lambda x:x[0]) #個々の波高を小さい順にsort(*reverse()ではkeyが使えなかった...)
wave_height = [x[0] for x in waves]
wave_period = [x[1] for x in waves]
left = [x for x in range(len(waves))]
plt.bar(left,wave_height)
plt.ylabel("wave height (m)")
plt.show()
全部で134波ありました。これから上位1/3の44波を平均し、有義波高、周期を求めます。
# 有義波を算出
n = int(len(waves)/3) #上位44波
H13 = mean(wave_height[-n:]) #0.81
P13 = mean(wave_period[-n:]) #10.9
# ついでに、全体の平均波と最大波も算出
Hmean = mean(wave_height[:]) #0.52
Pmean = mean(wave_period[:]) #8.9
Hmax = mean(wave_height[-1:]) #1.43
Pmax = mean(wave_period[-1:]) #11.0
以上から、有義波を算出し、今回のサンプルデータは、波高0.81m、周期10.9秒と算出されました。
＜サンプルデータ（再掲）＞
まとめ
波の波高・周期は、波高の上位1/3の波の平均である有義波という定義で表すことができ、熟練の観測者が目視で観測する波高や周期に近い数値になると言われています。
今回は、Pythonを使い、ゼロアップクロス法を用いて波浪の統計量を算出し、平均波、有義波、最大波の関係が以下のようになりました。
| 波高 | 周期 | 定義 | |
|---|---|---|---|
| 平均波 | 0.52m | 8.9秒 | 全体134波の平均波高、平均周期 |
| 有義波 | 0.81m | 10.9秒 | 波高上位44波の平均波高、平均周期 |
| 最大波 | 1.43m | 11.0秒 | 波高が最大となる波の波高、周期 |
また機会がございましたら、次回はFFT（高速フーリエ変換）によるスペクトル解析を用いて、波浪解析をPythonでやってみたいと思います。