教師なし機械学習で地震波をクラスタリングする

Pythonの機械学習ライブラリ、scikit-learnの教師なし機械学習で、地震波をクラスタリングしてみました。

地震波は、K-NETから2011年3月11日の地震波をダウンロードして、予めランダムに選択した地震波の応答スペクトルを計算しておき、これをデータにして5つのクラスタに分けます。

以下がプログラムです。

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()

データを準備したり、結果を描画するための処理が多いのですが、機械学習自体は以下のたった1行で完結します。

model = KMeans(n_clusters=5, init='random').fit(feature)

このプログラムでグループ分けされた応答スペクトルを描画すると以下のようになります。

応答スペクトルのクラスタリング(1)

少し分かりづらいですが、greenのグループが応答スペクトルが大きいグループ、blueのグループが応答スペクトルが小さいグループで、blue < cyan < red < magenta < green の順で大きくなっています。

これらの地震波を観測点ごとに、上図と同じ色分けで日本地図に描画すると以下のようになります。
(地図データは国土地理院のデータを使用しています。)

観測点プロット(1)

容易に想像できますが、震源に近い観測点で大きな応答スペクトルになっていて、遠いところで小さな応答スペクトルになっています。

さて、地震波はその特徴として周期特性がありますが、この周期特性を元にしてクラスタリングする為に、応答スペクトルを地震波の最大値で基準化してグループ分けしてみます。

地震波の時間刻みが0.01秒であれば、解像度は50Hzで、応答スペクトルの最短周期の0.02秒が地震波の最大値になるので、この値で応答スペクトル全体を徐算します。

プログラムは以下の1行のコメントを外して入れ替えます。

# feature = np.append(feature, df['acc'])
feature = np.append(feature, df['acc']/df['acc'][0])

この結果を描画したものが以下です。

応答スペクトルのクラスタリング(2)

redが長周期成分が卓越しているグループ、greenが短周期成分が卓越しているグループで、green < cyan < magenta < blue < red の順で卓越している周期成分が長周期寄りになります。

最初のグループ分けと同様にこれらを日本地図に描画したものが以下です。

観測点プロット(2)

最初のケースと比較すると、各グループの際立った地域特性は発見できません。これは、観測波の周期特性はその点の地盤の特性が大きく影響するためだと思います。

北海道や関西に長周期寄りの観測波が多いのは何か理由があるのかもしれません。

scikit-learnの公式サイトには豊富なサンプルがあるので、直ぐに機械学習を試すことができます。

※2021/10/30追記
長周期成分が卓越している記録波が遠方に多いのは、長周期成分は遠方まで到達しやすいからではないかとのコメントをいただきました。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

CAPTCHA