把測站波形丟進 ObsPy 前，先確認取樣率

上個月幫朋友檢查一份 2023 年的強震儀紀錄，儀器是 Basalt，取樣率設 200 Hz，事件是規模 5 左右、震央距離約 40 公里的近場紀錄。他用 ObsPy 讀進來畫頻譜，結果在 25 Hz 附近出現一個很平的高台，跑濾波之後整段波形都在抖。檔案本身沒壞，問題出在讀檔時沒做 detrend，也沒注意儀器內建的 anti-alias 行為。

我先把最小可用的檢查流程貼出來，這段可以直接放進 script：

from obspy import read
st = read("basalt_2023.mseed")
print(st)
for tr in st:
    print(tr.id, tr.stats.sampling_rate, tr.stats.npts, tr.data.dtype)
st.detrend("linear")
st.taper(0.02)
st.filter("highpass", freq=0.5, corners=4, zerophase=True)

關鍵是 st.detrend 與 taper 要在濾波前做。線性去趨勢把感測器零點漂移拿掉，taper 則避免濾波時兩端邊界反射。zerophase=True 是為了不讓相位延遲影響到 P 波到時的判讀，代價是前後各會多一點邊緣效應，所以 taper 不能省。

為什麼會出現 25 Hz 的高台：200 Hz 取樣率的 Nyquist 是 100 Hz，但這台儀器在 25 Hz 以上有明顯的數位濾波滾降，頻譜上看到的是儀器響應，不是地動。要驗證很簡單，把 tr.stats.response 印出來，或直接對照儀器手冊的響應曲線。如果沒有響應檔，至少要知道這條曲線的存在，不然把 25 Hz 當成訊號去解釋，結論會完全走偏。

失效的地方有幾個。第一，取樣率如果不是整數倍，taper 的比例要重算，固定 0.02 在短紀錄上會吃掉太多資料。第二，zerophase=True 搭配 lowpass 在高頻端容易出現 ringing，尤其是資料長度不到 10 秒的餘震紀錄。第三，如果檔案是多通道合併的 miniSEED，st.detrend 會逐 trace 處理，但 st.filter 的預設是整段 stream 一起做，通道數一多記憶體會爆。我自己踩過一次，8 通道、取樣率 500 Hz、長度 30 分鐘，濾波時直接吃掉快 2 GB。

後來我的習慣是先把每個 trace 的 stats.sampling_rate 與 npts 印出來，再決定要不要分段濾波。這一步很笨，但省掉後面重跑的時間。

你們在處理近場強震紀錄時，會先看儀器響應曲線，還是直接憑取樣率判斷可用頻帶？