Pythonの機械学習ライブラリ、scikit-learnの教師なし機械学習で、地震波をクラスタリングしてみました。
地震波は、K-NETから2011年3月11日の地震波をダウンロードして、予めランダムに選択した地震波の応答スペクトルを計算しておき、これをデータにして５つのクラスタに分けます。
以下がプログラムです。
import glob
import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt
import geopandas as gpd
if __name__ == '__main__':
# 予め計算しておいた各記録波の応答スペクトルを読込む
files = glob.glob('*.spec.csv')
feature = np.array([])
for file in files:
df = pd.read_csv(file, index_col=None)
feature = np.append(feature, df['acc'])
# feature = np.append(feature, df['acc']/df['acc'][0])
# 1次元numpy配列を2次元に変換する（応答スペクトルの計算ポイント数300）
feature = feature.reshape(-1, 300)
# 5種類のグループにクラスタリングする
model = KMeans(n_clusters=5, init='random').fit(feature)
# 学習結果のラベルを取得する
labels = model.labels_
# 応答スペクトルをクラスタ毎に色分けしてプロットする
colors = ["red", "green", "blue", "magenta", "cyan"]
plt.grid()
plt.xscale('log')
plt.xlabel('Period(sec)')
plt.ylabel('Response Acceleration Spectrum')
for label, v, file in zip(labels, feature, files):
plt.plot(df['period(sec)'], v, color=colors[label])
plt.show()
# 各観測点の情報を読込む
geo = pd.read_csv('geo.csv', index_col='name')
# クラスタ毎に色分けして地図上にプロットする
jpnShp = gpd.read_file('./gm-jpn-all_u_2_2/coastl_jpn.shp')
ax = jpnShp.plot(figsize=(10, 20), color='black', linewidth=.5)
for label, file in zip(labels, files):
plt.scatter(geo.at[file[:6],'Long'], geo.at[file[:6],'Lat'], marker='o', s=10, color=colors[label], alpha=1.0, linewidths = 1)
plt.show()データを準備したり、結果を描画するための処理が多いのですが、機械学習自体は以下のたった１行で完結します。
model = KMeans(n_clusters=5, init='random').fit(feature)このプログラムでグループ分けされた応答スペクトルを描画すると以下のようになります。
少し分かりづらいですが、greenのグループが応答スペクトルが大きいグループ、blueのグループが応答スペクトルが小さいグループで、blue < cyan < red < magenta < green の順で大きくなっています。
これらの地震波を観測点ごとに、上図と同じ色分けで日本地図に描画すると以下のようになります。
（地図データは国土地理院のデータを使用しています。）
容易に想像できますが、震源に近い観測点で大きな応答スペクトルになっていて、遠いところで小さな応答スペクトルになっています。
さて、地震波はその特徴として周期特性がありますが、この周期特性を元にしてクラスタリングする為に、応答スペクトルを地震波の最大値で基準化してグループ分けしてみます。
地震波の時間刻みが0.01秒であれば、解像度は50Hzで、応答スペクトルの最短周期の0.02秒が地震波の最大値になるので、この値で応答スペクトル全体を徐算します。
プログラムは以下の１行のコメントを外して入れ替えます。
# feature = np.append(feature, df['acc'])
feature = np.append(feature, df['acc']/df['acc'][0])この結果を描画したものが以下です。
redが長周期成分が卓越しているグループ、greenが短周期成分が卓越しているグループで、green < cyan < magenta < blue < red の順で卓越している周期成分が長周期寄りになります。
最初のグループ分けと同様にこれらを日本地図に描画したものが以下です。
最初のケースと比較すると、各グループの際立った地域特性は発見できません。これは、観測波の周期特性はその点の地盤の特性が大きく影響するためだと思います。
北海道や関西に長周期寄りの観測波が多いのは何か理由があるのかもしれません。
scikit-learnの公式サイトには豊富なサンプルがあるので、直ぐに機械学習を試すことができます。
※2021/10/30追記
長周期成分が卓越している記録波が遠方に多いのは、長周期成分は遠方まで到達しやすいからではないかとのコメントをいただきました。