Note

このページは Jupyter Notebook から生成されました。 ノートブックをダウンロード (.ipynb)

ObsPy 連携による地震データ解析

このチュートリアルでは MiniSEED / SAC / FDSN 形式の地震データを読み込み、gwexpy で解析する方法と、obspy ↔ gwexpy 間の双方向変換を解説します。

重力波検出器科学での用途

• KAGRA・LIGO・Virgo 環境モニターの地震計データ読み込み

• 地震ノイズ解析と検出器感度へのカップリング

• 地震チャンネルと重力波データのクロス相関

レガシーコード移行メモ:

obspy.read() で読み込んで手動で numpy 処理していたコードは、

TimeSeries.read(format='miniseed') + gwexpy の高レベル API に置き換えられます。

API_MAPPING.md Category 4 の interop リファレンスも参照してください。

import numpy as np
import matplotlib.pyplot as plt
from astropy import units as u

from gwexpy.timeseries import TimeSeries, TimeSeriesDict

try:
    import obspy
    OBSPY_AVAILABLE = True
    print(f"ObsPy バージョン: {obspy.__version__}")
except ImportError:
    OBSPY_AVAILABLE = False
    print("ObsPy がインストールされていません。pip install obspy でインストールしてください。")
    print("このチュートリアルでは API の紹介のみ行い、obspy が必要なセルはスキップします。")

1. ルート 1: TimeSeries.read(format='miniseed')

gwexpy の I/O レイヤーは MiniSEED ファイルを直接読み込み、正しい単位と時間軸を持つ TimeSeries を返します。

# --- ルート 1: gwexpy I/O 経由で MiniSEED を直接読み込み ---
# 実際のファイルパスに置き換えてください:
#
# ts_seismic = TimeSeries.read(
#     "path/to/seismic.mseed",
#     format='miniseed',
# )
#
# SAC 形式:
# ts_seismic = TimeSeries.read("path/to/seismic.sac", format='sac')
#
# FDSN ウェブサービス（ネットワーク接続が必要）:
# ts_seismic = TimeSeries.read(
#     "IU.ANMO.00.BHZ",
#     format='miniseed',
#     start=1234567890,
#     end=1234568090,
# )

print("ルート 1: TimeSeries.read(format='miniseed') または format='sac'")
print("GPS 時間軸・単位・チャンネル名を持つ TimeSeries が返されます。")

2. ルート 2: from_obspy_trace() – obspy オブジェクトから変換

すでに obspy.Trace / obspy.Stream がある場合は from_obspy_trace() で gwexpy に変換します。

if OBSPY_AVAILABLE:
    rng = np.random.default_rng(0)
    fs_seis = 100.0
    n_seis  = 6000
    t_seis  = np.arange(n_seis) / fs_seis

    # 地震ノイズのシミュレーション
    seis_data = (
        1e-7 * rng.normal(0, 1, n_seis)
        + 5e-7 * np.sin(2 * np.pi * 0.15 * t_seis)  # 0.15 Hz マイクロセイズム
        + 2e-7 * np.sin(2 * np.pi * 1.0  * t_seis)  # 1 Hz ノイズ
    )

    stats = obspy.core.Stats()
    stats.network  = "K1"
    stats.station  = "KAGRA"
    stats.location = "00"
    stats.channel  = "BHZ"
    stats.sampling_rate = fs_seis
    stats.starttime = obspy.UTCDateTime("2023-01-01T00:00:00")
    stats.npts = n_seis

    tr = obspy.Trace(data=seis_data, header=stats)
    print("obspy Trace:", tr)

    # gwexpy TimeSeries に変換
    ts_seismic = TimeSeries.from_obspy_trace(tr, unit=u.m / u.s)
    print(f"\nt0={ts_seismic.t0:.3f}, dt={ts_seismic.dt}, N={len(ts_seismic.value)}, unit={ts_seismic.unit}")
else:
    # obspy なしのフォールバック
    rng = np.random.default_rng(0)
    fs_seis = 100.0
    n_seis  = 6000
    t_seis  = np.arange(n_seis) / fs_seis
    seis_data = (
        1e-7 * rng.normal(0, 1, n_seis)
        + 5e-7 * np.sin(2 * np.pi * 0.15 * t_seis)
        + 2e-7 * np.sin(2 * np.pi * 1.0  * t_seis)
    )
    ts_seismic = TimeSeries(
        seis_data, dt=(1/fs_seis)*u.s, t0=0*u.s,
        unit=u.m/u.s, name="K1:KAGRA.BHZ",
    )
    print("合成 TimeSeries を作成しました（obspy なし）")

3. 地震データの信号処理

ts_seismic.plot(xscale='seconds');

# 帯域フィルタリング
ts_lf = ts_seismic.bandpass(0.05, 0.5)   # 0.05–0.5 Hz: 地球の固有振動 + マイクロセイズム
ts_hf = ts_seismic.bandpass(1.0,  10.0)  # 1–10 Hz: 人工ノイズ帯

fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True)
axes[0].plot(np.arange(len(ts_seismic.value)) / fs_seis, ts_seismic.value, lw=0.5, label="Raw")
axes[1].plot(np.arange(len(ts_lf.value))      / fs_seis, ts_lf.value,      lw=0.8, label="0.05–0.5 Hz", color="C1")
axes[2].plot(np.arange(len(ts_hf.value))      / fs_seis, ts_hf.value,      lw=0.8, label="1–10 Hz",    color="C2")

for ax in axes:
    ax.set_ylabel(f"[{ts_seismic.unit}]")
    ax.legend(loc="upper right")
    ax.grid(True, alpha=0.3)
axes[2].set_xlabel("時間 [s]")
fig.suptitle("地震データ: 帯域フィルタリング")
plt.tight_layout()

# 振幅スペクトル密度（ASD）
asd = ts_seismic.asd(fftlength=10.0, overlap=0.5)

fig, ax = plt.subplots(figsize=(8, 4))
ax.loglog(asd.frequencies.value, asd.value)
ax.set_xlabel("周波数 [Hz]")
ax.set_ylabel(f"ASD [{asd.unit}]")
ax.set_title("地震 ASD")
ax.grid(True, which="both", alpha=0.3)
ax.axvline(0.15, color="C1", linestyle="--", label="0.15 Hz マイクロセイズム")
ax.axvline(1.0,  color="C2", linestyle="--", label="1.0 Hz ノイズ")
ax.legend()
plt.tight_layout()

4. obspy に戻す変換

if OBSPY_AVAILABLE:
    tr_filtered = ts_lf.to_obspy_trace()
    print("obspy Trace に変換:", tr_filtered)

    from gwexpy.interop import to_obspy
    tr_generic = to_obspy(ts_lf)
    print("to_obspy() 経由:", tr_generic)
else:
    print("obspy が利用できないため、ラウンドトリップ変換をスキップします")

5. 移行ガイド: obspy + numpy → gwexpy

移行前（obspy + scipy）:

import obspy
from scipy import signal

st = obspy.read("seismic.mseed")
tr = st[0]

tr_bp = tr.copy().filter('bandpass', freqmin=0.1, freqmax=1.0)
f, psd = signal.welch(tr.data, fs=tr.stats.sampling_rate, nperseg=1024)
asd = np.sqrt(psd)   # 単位なし

移行後（gwexpy）:

from gwexpy.timeseries import TimeSeries

ts = TimeSeries.read("seismic.mseed", format='miniseed')
# または: ts = TimeSeries.from_obspy_trace(obspy.read("seismic.mseed")[0])

ts_bp = ts.bandpass(0.1, 1.0)
asd   = ts.asd(fftlength=10.0, overlap=0.5)  # 単位・軸情報付き
asd.plot()  # そのままプロット可能

まとめ

操作	gwexpy API
MiniSEED 読み込み	TimeSeries.read(path, format='miniseed')
SAC 読み込み	TimeSeries.read(path, format='sac')
obspy Trace から変換	TimeSeries.from_obspy_trace(tr, unit=...)
obspy Trace に変換	ts.to_obspy_trace() または to_obspy(ts)
帯域フィルタ	ts.bandpass(fmin, fmax)
ASD	ts.asd(fftlength=..., overlap=...)
スペクトログラム	ts.spectrogram(stride=..., fftlength=...)

関連チュートリアル:

• 相互運用性の基礎

• カップリング関数解析