Note

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

GW 検出器におけるシューマン共鳴解析

シューマン共鳴は、地球と電離層のキャビティ内で落雷によって励起されるグローバルな電磁共鳴です。 約 7.83、14.3、20.8、26.4 Hz に鋭いスペクトルピークとして現れ、KAGRA などの第 2 世代重力波検出器の 感度帯域に直接重なります。

このノートブックでは gwexpy を用いたエンドツーエンドの特性評価ワークフローを示します:

1. ブートストラップ PSD (bootstrap_spectrogram) — 非対称信頼区間付きのロバストなスペクトル推定

2. ローレンツ Q ファクターフィット (fit_series + lorentzian_q) — 各モードの共鳴周波数・Q 値・振幅を計測

3. 共分散構造 (BifrequencyMap) — ブートストラップが返す周波数ビン間相関を可視化し、GLS フィットに利用

4. 時間追跡 — 観測ウィンドウ内でのモード振幅の時間変化を監視

前提知識:

• スペクトログラム基礎

• フィッティング上級編

• ブートストラップ PSD & GLS フィット

セットアップ

# ruff: noqa: I001
import matplotlib.pyplot as plt
import numpy as np
from astropy import units as u

from gwexpy.fitting import fit_series
from gwexpy.fitting.models import lorentzian_q
from gwexpy.frequencyseries import BifrequencyMap  # noqa: F401 (参照用)
from gwexpy.spectral import bootstrap_spectrogram
from gwexpy.timeseries import TimeSeries

1. モックデータ: シューマン共鳴を含む磁力計データ

地球の背景電磁場を測定する 1 軸磁力計 (ELF 帯) をシミュレートします。 4 つのシューマン共鳴を ローレンツ型スペクトル形状の狭帯域ノイズ として生成し、 広帯域ノイズフロアに重ね合わせます。

モード	周波数	Q 値	ピーク ASD
SR1	7.83 Hz	5.0	4.0 nT/√Hz
SR2	14.3 Hz	4.5	2.5 nT/√Hz
SR3	20.8 Hz	4.0	1.5 nT/√Hz
SR4	26.4 Hz	3.5	1.0 nT/√Hz

rng = np.random.default_rng(42)
fs = 512        # サンプリングレート [Hz]
T  = 300.0      # 継続時間 [s]
n  = int(fs * T)

# シューマン共鳴パラメータ (地球-電離層キャビティモード)
SR_FREQS = [7.83, 14.3, 20.8, 26.4]   # Hz — 第 1〜4 モード
SR_Q     = [5.0,  4.5,  4.0,  3.5]    # Q 値
SR_AMP   = [4.0,  2.5,  1.5,  1.0]    # ピーク ASD [nT/√Hz]
NOISE_FLOOR = 0.3                       # 広帯域フロア [nT/√Hz]

# ── 周波数領域合成 ───────────────────────────────────────────────────────────
f = np.fft.rfftfreq(n, d=1.0 / fs)   # 片側周波数軸 [Hz]

# 目標 ASD をローレンツ和 + フロアで構築
asd_target = np.full_like(f, NOISE_FLOOR)
for f0, q, A in zip(SR_FREQS, SR_Q, SR_AMP):
    gamma = f0 / (2.0 * q)            # 半値半幅 [Hz]
    asd_target += A * gamma / np.sqrt((f - f0) ** 2 + gamma ** 2)

# rfft 振幅に変換: Welch PSD ≈ asd_target² となるよう正規化
amp = asd_target * np.sqrt(n * fs / 2.0)
amp[0] = 0.0  # DC 成分をゼロに

Z = (rng.standard_normal(len(f)) + 1j * rng.standard_normal(len(f))) / np.sqrt(2)
x = np.fft.irfft(amp * Z, n=n)

mag = TimeSeries(x, dt=1.0 / fs, unit=u.nT,
                 name="K1:PEM-MAG_EXV_EAST_X_DQ", t0=0)
print(f"継続時間: {T:.0f} s | fs: {fs} Hz | N: {n:,}")

# クイック確認: ASD でシューマン共鳴ピークを確認
fig, ax = plt.subplots(figsize=(10, 4))
asd_raw = mag.asd(fftlength=16.0, overlap=8.0)
ax.semilogy(asd_raw.frequencies.value, asd_raw.value, lw=0.8,
            color='steelblue', label='単一 ASD 推定 (16 s FFT)')

for f0, lab in zip(SR_FREQS, ['SR1', 'SR2', 'SR3', 'SR4']):
    ax.axvline(f0, color='red', ls='--', lw=0.8, alpha=0.7)
    ax.text(f0 + 0.15, 0.35, lab, color='red', fontsize=8)

ax.set_xlim(4, 35)
ax.set_xlabel('周波数 [Hz]')
ax.set_ylabel('ASD [nT/√Hz]')
ax.set_title('磁力計 ASD — シューマン共鳴の 4 ピークが確認できる')
ax.legend()
plt.tight_layout()
plt.show()

2. ブートストラップ PSD 推定

bootstrap_spectrogram はスペクトログラムの時間カラムをリサンプリングし、 ロバストな PSD 推定値 (中央値または平均値) と非対称信頼区間を提供します。 return_map=True にすると、周波数ビン間の相関を定量化する 共分散 ``BifrequencyMap`` cov_map(f1, f2) も返します。これは次のステップの GLS フィットに使用されます。

# スペクトログラム計算: 16 s FFT、50 % ハン窓オーバーラップ
spec = mag.spectrogram2(fftlength=16.0, overlap=8.0, window='hann')
print(f"スペクトログラム shape: {spec.shape}  (n_times × n_freqs)")

# ブートストラップ — (PSD, 共分散 BifrequencyMap) を返す
psd_boot, cov_map = bootstrap_spectrogram(
    spec,
    n_boot=500,
    method='median',
    ci=0.68,
    fftlength=16.0,
    overlap=8.0,
    return_map=True,
    ignore_nan=True,
)
print(f"ブートストラップ PSD shape : {psd_boot.shape}")
print(f"共分散マップ shape: {cov_map.shape}  ← BifrequencyMap(f1, f2)")

# ブートストラップ PSD と ±1σ 信頼区間をプロット
fig, ax = plt.subplots(figsize=(10, 4))

f_psd = psd_boot.frequencies.value
y_psd = psd_boot.value

diag = cov_map.diagonal(method='mean')
y_var = np.interp(f_psd, diag.frequencies.value, np.abs(diag.value))
y_lo = np.sqrt(np.maximum(y_psd - np.sqrt(y_var), 1e-12))
y_hi = np.sqrt(y_psd + np.sqrt(y_var))

ax.semilogy(f_psd, np.sqrt(y_psd), lw=1.5, color='steelblue',
            label='ブートストラップ中央値 ASD')
ax.fill_between(f_psd, y_lo, y_hi, alpha=0.25, color='steelblue', label='±1σ')

for f0 in SR_FREQS:
    ax.axvline(f0, color='red', ls='--', lw=0.7, alpha=0.6)

ax.set_xlim(4, 35)
ax.set_xlabel('周波数 [Hz]')
ax.set_ylabel('ASD [nT/√Hz]')
ax.set_title('ブートストラップ PSD — 中央値 ± 1σ 信頼区間')
ax.legend()
plt.tight_layout()
plt.show()

3. ローレンツ Q ファクターフィット

各シューマンモードを Q ファクター型ローレンツ分布 でモデル化します:

\[S(f) = \frac{A\,\gamma^2}{(f - f_0)^2 + \gamma^2}, \quad \gamma = \frac{f_0}{2Q}\]

fit_series に cov=cov_map を渡すことで、オーバーラップ FFT セグメント間の 相関を考慮した 一般化最小二乗法 (GLS) でフィットを行い、正しいパラメータ誤差が得られます。

# 各シューマンモードを専用の周波数ウィンドウで個別にフィット
fit_ranges = [(6.0, 10.5), (11.5, 17.5), (17.5, 24.5), (23.0, 30.5)]

fit_results = []
for i, (f0, q0, A0, (flo, fhi)) in enumerate(
        zip(SR_FREQS, SR_Q, SR_AMP, fit_ranges)):
    result = fit_series(
        psd_boot,
        lorentzian_q,
        x_range=(flo, fhi),
        cov=cov_map,
        p0={'A': A0 ** 2, 'x0': f0, 'Q': q0},
        limits={'A': (0, 500), 'x0': (flo, fhi), 'Q': (1.0, 100.0)},
    )
    fit_results.append(result)
    p, e = result.params, result.errors
    print(
        f"SR{i + 1}: f0 = {p['x0']:.3f} ± {e['x0']:.3f} Hz  |"
        f"  Q = {p['Q']:.2f} ± {e['Q']:.2f}  |"
        f"  A = {p['A']:.4f} ± {e['A']:.4f} nT²/Hz"
    )

# フィット済みローレンツ曲線をブートストラップ PSD に重ねてプロット
fig, axes = plt.subplots(1, 4, figsize=(15, 4))
colors = ['steelblue', 'darkorange', 'seagreen', 'crimson']
plot_ranges = [(5.5, 11.0), (11.5, 17.5), (17.5, 24.5), (22.5, 30.5)]
labels = ['SR1 (7.83 Hz)', 'SR2 (14.3 Hz)', 'SR3 (20.8 Hz)', 'SR4 (26.4 Hz)']

for ax, result, fr, color, label in zip(axes, fit_results, plot_ranges, colors, labels):
    psd_crop = psd_boot.crop(*fr)
    f_crop = psd_crop.frequencies.value
    y_crop = psd_crop.value

    ax.semilogy(f_crop, y_crop, '.', color=color, ms=2.5, label='ブートストラップ PSD')

    f_fine = np.linspace(*fr, 300)
    ax.semilogy(f_fine, lorentzian_q(f_fine, **result.params),
                'k-', lw=1.5, label='ローレンツフィット')

    p = result.params
    ax.set_title(f"{label}\nf₀={p['x0']:.3f} Hz, Q={p['Q']:.2f}", fontsize=9)
    ax.set_xlabel('周波数 [Hz]')
    ax.legend(fontsize=7)

axes[0].set_ylabel('PSD [nT²/Hz]')
plt.suptitle('シューマン共鳴モードへの GLS ローレンツフィット (cov_map 使用)', y=1.01)
plt.tight_layout()
plt.show()

4. BifrequencyMap による共分散構造の解析

bootstrap_spectrogram が返す共分散マップ cov_map(f1, f2) は、 異なる周波数のスペクトル推定値がどの程度相関しているか を表します。 狭帯域共鳴を持つ磁力計データでは以下が期待されます:

• 共鳴周波数における対角線沿いの高共分散 (f1 ≈ f2) — オーバーラップ FFT ウィンドウから生じる隣接ビン間の相関

• 共鳴付近のオフ対角構造 — ピークが 1 ビンを超えて広がることを反映

ここで使用する BifrequencyMap メソッド:

• .diagonal(method='mean') — 分散プロファイルを返す (行列の対角線)

• .get_slice(at=f0, axis='f1') — f1 固定での 1 次元共分散スライスを返す

from matplotlib.colors import LogNorm

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# ── 左: 2-D 共分散マップ (4〜35 Hz) ──────────────────────────────────────
f_lo, f_hi = 4.0, 35.0
f_vals = cov_map.frequency1.value
mask = (f_vals >= f_lo) & (f_vals <= f_hi)
cov_sub = np.abs(cov_map.value[np.ix_(mask, mask)])
f_sub   = f_vals[mask]

im = axes[0].pcolormesh(f_sub, f_sub, cov_sub,
                         norm=LogNorm(vmin=cov_sub[cov_sub > 0].min(), vmax=cov_sub.max()),
                         cmap='inferno', shading='auto')
fig.colorbar(im, ax=axes[0], label='|共分散| [nT⁴/Hz²]')

for f0 in SR_FREQS:
    axes[0].axvline(f0, color='cyan', lw=0.6, ls='--')
    axes[0].axhline(f0, color='cyan', lw=0.6, ls='--')

axes[0].set_xlabel('f₁ [Hz]')
axes[0].set_ylabel('f₂ [Hz]')
axes[0].set_title('ブートストラップ共分散マップ (4〜35 Hz)')

# ── 右: 対角線 = ビン毎の分散 ────────────────────────────────────────────
diag = cov_map.diagonal(method='mean')
f_d = diag.frequencies.value
m_d = (f_d >= f_lo) & (f_d <= f_hi)
axes[1].semilogy(f_d[m_d], np.abs(diag.value[m_d]), lw=1.2, color='steelblue')

for f0 in SR_FREQS:
    axes[1].axvline(f0, color='red', ls='--', lw=0.7)

axes[1].set_xlim(f_lo, f_hi)
axes[1].set_xlabel('周波数 [Hz]')
axes[1].set_ylabel('対角分散 [nT⁴/Hz²]')
axes[1].set_title('共分散マップの対角線 (≈ PSD² / N_bootstrap)')

plt.tight_layout()
plt.show()

# SR1 (7.83 Hz) における共分散スライス — 共鳴の周波数幅を反映
cov_slice = cov_map.get_slice(at=7.83, axis='f1')
f_sl = cov_slice.frequencies.value
m_sl = (f_sl >= f_lo) & (f_sl <= f_hi)

fig, ax = plt.subplots(figsize=(10, 4))
ax.semilogy(f_sl[m_sl], np.abs(cov_slice.value[m_sl]), lw=1.2, color='darkorange')
ax.axvline(7.83, color='red', ls='--', lw=1, label='SR1 (7.83 Hz)')
for f0 in SR_FREQS[1:]:
    ax.axvline(f0, color='gray', ls=':', lw=0.7)

ax.set_xlim(f_lo, f_hi)
ax.set_xlabel('周波数 f₂ [Hz]')
ax.set_ylabel('|Cov(f₁=7.83 Hz, f₂)| [nT⁴/Hz²]')
ax.set_title('SR1 における共分散スライス — ピーク幅が共鳴の帯域幅を反映')
ax.legend()
plt.tight_layout()
plt.show()

5. 振幅の時間追跡

シューマン共鳴の振幅は世界的な落雷活動に応じて変動します (日変化・季節変化)。 スペクトログラムから各モード周辺の帯域内平均パワーを抽出し、 ASD 時系列に変換して時間変化を監視します。

fig, ax = plt.subplots(figsize=(10, 4))
times_s = spec.times.value
colors_t = ['steelblue', 'darkorange', 'seagreen', 'crimson']

for f0, color, label in zip(SR_FREQS, colors_t,
                             ['SR1 7.83 Hz', 'SR2 14.3 Hz',
                              'SR3 20.8 Hz', 'SR4 26.4 Hz']):
    # ±0.5 Hz 帯域内の狭帯域パワー
    spec_band = spec.crop_frequencies(f0 - 0.5, f0 + 0.5)
    amp_t = np.sqrt(spec_band.value.mean(axis=1))
    ax.plot(times_s, amp_t, lw=0.8, color=color, label=label)

ax.set_xlabel('時間 [s]')
ax.set_ylabel('帯域内平均 ASD [nT/√Hz]')
ax.set_title('シューマン共鳴振幅の時間変化 (300 s モック観測)')
ax.legend(ncol=2, fontsize=9)
plt.tight_layout()
plt.show()

まとめ

ステップ	ツール	出力
ロバスト PSD	bootstrap_spectrogram(return_map=True)	中央値 PSD + 1σ 帯
不確かさモデル	BifrequencyMap (ブートストラップから)	周波数ビン間共分散行列
ピーク特性評価	fit_series(lorentzian_q, cov=cov_map)	各モードの f₀、Q、振幅
時間監視	Spectrogram.crop_frequencies()	振幅時系列

ポイント

• ブートストラップリサンプリングは、単一セグメント推定がノイジーな狭帯域ピークに対して 信頼できる不確かさ推定を与えます。

• BifrequencyMap 共分散は GLS フィットに不可欠です。fit_series に cov=cov_map を 渡すことで、周波数ビン間の相関が自動的に考慮されます。

• lorentzian_q パラメタライゼーションにより、各シューマンモードの物理的に意味のある Q ファクターが直接得られます。

• このワークフローは実際の KAGRA PEM 磁力計データにそのまま適用できます。

次のステップ

• 実データの読み込み: TimeSeries.read('K1:PEM-MAG_EXV_EAST_X_DQ', start, end).

• BifrequencyMap.propagate() と組み合わせて、DARM ノイズへの磁気結合寄与を推定する。

• GPS 時間セグメントをループして、シューマン共鳴パラメータの日変化・季節変化を追跡する。