obspy

ObsPy - 地震学数据处理

快速参考

from obspy import read, UTCDateTime
from obspy.clients.fdsn import Client

# Read local file (MiniSEED, SAC, etc.)
st = read("data.mseed")
tr = st[0]                          # First trace
print(tr.stats)                     # Metadata

# Fetch from FDSN
client = Client("IRIS")
t = UTCDateTime("2023-02-06T01:17:00")
st = client.get_waveforms("IU", "ANMO", "00", "LHZ", t, t + 3600)
st.plot()

主要类

类	用途
Stream	多个 Trace 对象的容器
Trace	包含数据和元数据的单条波形
UTCDateTime	精确时间处理
Inventory	台站／通道元数据
Catalog	地震事件信息

基本操作

读取和检查

st = read("data.mseed")              # Auto-detect format
tr = st[0]
print(tr.stats.station, tr.stats.channel, tr.stats.sampling_rate)

滤波和处理

st.detrend("demean")                 # Remove mean
st.detrend("linear")                 # Remove trend
st.taper(max_percentage=0.05)        # Taper edges
st.filter("bandpass", freqmin=0.1, freqmax=10.0)

从 FDSN 获取波形

client = Client("IRIS")
t1 = UTCDateTime("2023-02-06T01:17:00")
st = client.get_waveforms("IU", "ANMO", "00", "LHZ", t1, t1 + 3600)
st.write("output.mseed", format="MSEED")

搜索地震

client = Client("USGS")
cat = client.get_events(
    starttime=UTCDateTime("2023-01-01"),
    endtime=UTCDateTime("2023-12-31"),
    minmagnitude=7.0
)
event = cat[0]
print(f"M{event.magnitudes[0].mag} at {event.origins[0].latitude}")

获取台站元数据

inv = client.get_stations(
    network="IU", station="ANMO",
    level="response"                 # Required for response removal
)

移除仪器响应

st = client.get_waveforms("IU", "ANMO", "00", "LHZ", t, t + 3600)
inv = client.get_stations(network="IU", station="ANMO", level="response")
st.remove_response(inventory=inv, output="VEL")  # VEL, DISP, or ACC

裁剪和选择

st.trim(UTCDateTime("2023-01-01"), UTCDateTime("2023-01-01T01:00:00"))
st_z = st.select(channel="*Z")       # Vertical only
st_bh = st.select(channel="BH*")     # BH channels

合并和处理间隙

st.print_gaps()                      # Check for gaps
st.merge(method=1, fill_value="interpolate")

写入数据

st.write("output.mseed", format="MSEED")
st.write("output.sac", format="SAC")

错误处理

from obspy.clients.fdsn.header import FDSNNoDataException

try:
    st = client.get_waveforms("IU", "ANMO", "00", "LHZ", t1, t2)
except FDSNNoDataException:
    print("No data available")

常见提示

1. 滤波前务必去趋势并加窗，以避免伪影
2. 获取台站信息以进行仪器校正时，使用 level="response"
3. 处理前使用 st.print_gaps() 检查间隙
4. 查询支持通配符：station="A*"、channel="BH?"
5. UTCDateTime 接受 ISO 字符串、时间戳和 datetime 对象

何时使用及替代方案

工具	最适用场景
obspy	完整的地震学工作流、FDSN 数据访问、波形处理
segyio	用于反射地震数据的快速、底层 SEG-Y 文件 I/O
scipy.signal	不涉及地震学专属功能的通用信号处理

需要地震学背景信息时使用 obspy：例如 FDSN 客户端、移除仪器 响应、地震目录或标准地震数据格式。

仅处理 SEG-Y 反射地震文件时改用 segyio。对于支持 inline/crossline 访问的大型 3D 数据体，segyio 速度更快。

只需要通用滤波或频谱分析，且不需要地震学专属元数据或数据访问时，改用 scipy.signal。

常见工作流

获取并处理地震波形

- [ ] Initialize FDSN client: `Client("IRIS")`
- [ ] Search events with `client.get_events()` for target magnitude/region
- [ ] Fetch waveforms with `client.get_waveforms()` for desired stations
- [ ] Get station metadata with `level="response"` for instrument correction
- [ ] Preprocess: detrend, taper, remove instrument response
- [ ] Filter to frequency band of interest
- [ ] Trim to analysis window and export or plot

参考资料

• FDSN 数据中心 - 可用数据中心及其功能
• 通道代码 - 通道命名规范详情
• 故障排查 - 常见问题及解决方案

脚本

• scripts/fetch_earthquake.py - 获取地震波形
• scripts/batch_process.py - 批量处理地震文件