把測站波形接進自家管線：一個小踩坑紀錄

上個月幫朋友整理一份山區小網格的強震紀錄，環境很普通：一台跑 Ubuntu 22.04 的舊筆電、Python 3.10，資料是四個測站連續三天的加速度檔，取樣率 100 Hz，換算下來約 2.6 GB 的 miniseed。目標不是做研究，只是想離線看 P 波到時跟幾個觸發窗，順便標一下可疑的雜訊段。

做法很土。先用 ObsPy 把檔案讀進來，切出要看的時間窗，再自己算一個簡單的 STA/LTA 當觸發器：

from obspy import read
from obspy.signal.trigger import classic_sta_lta, trigger_onset

st = read("sta01.mseed")
tr = st.select(component="Z")[0]
tr.detrend("linear")
tr.filter("bandpass", freqmin=1.0, freqmax=20.0, corners=2)

sta, lta = 0.5, 15.0
ratio = classic_sta_lta(tr.data, int(sta * tr.stats.sampling_rate),
                        int(lta * tr.stats.sampling_rate))
on = trigger_onset(ratio, 3.0, 1.5)
print(len(on), on[:5])

這段能直接跑，但真正花時間的不是程式，是參數。1–20 Hz 的帶通對近距離小震還行，對遠一點的體波就偏窄；STA 0.5 秒在 100 Hz 下只有 50 個樣點，短窗太短，單一尖峰就會誤觸。我把 STA 拉到 1 秒、門檻提到 3.5，誤觸才降下來，代價是尾端一些小事件被吃掉。

為什麼要自己算而不用現成套件？因為我需要看每個觸發窗前後的原始波形，套件吐出來的清單對不上我要標記的時間軸。自己算的好處是中間變數都留著，ratio 陣列可以直接疊在圖上。缺點也很明顯：這不是偵測器，只是門檻法，遇到連續背景振動或儀器階躍就會失效。我踩過的坑是某站換電池那段有一步直流跳動，filter 之後仍留下低頻擺動，ratio 一路偏高，最後是手動把那 40 秒標成無效。

另一個邊界是檔案格式。miniseed 裡同一個 channel 可能被切成很多段，read 進來是一個 Stream，直接取 [0] 只拿到第一段。正確做法是先 merge，缺資料的地方會補成 masked array，後面的 filter 就要處理填充值，不然觸發結果會出現假的突起。

我目前的做法是先 merge、檢查 gaps、再逐段處理，寧可慢一點也不要讓填充值混進 ratio。

你們在處理多段 miniseed 時，是習慣先轉成單一連續檔，還是保留原始分段、在流程裡自己管時間軸？