仪器响应#
利用 SAC 零极点文件去除仪器响应#
SAC 零极点文件是一种常见的仪器响应格式。若有 SAC 零极点格式的仪器响应文件，
则可以使用 obspy.io.sac.sacpz.attach_paz() 函数将对应的零极点文件添加到
Trace 的 stats.paz 属性中，再使用
Trace.simulate 函数去除仪器响应。
首先准备一些示例所需的波形数据和 SAC 零极点文件：
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
from obspy.clients.iris.client import Client as irisClient
starttime=UTCDateTime("2022-09-22T06:18:00")
# 下载 IU.ANMO 台站的三分量波形数据
client = Client("IRIS")
st = client.get_waveforms("IU", "ANMO", "00", "BH?", starttime, starttime + 720)
# 下载 IU.ANMO 台站三分量对应的 SAC 零极点文件
irisclient = irisClient()
irisclient.sacpz("IU", "ANMO", "00", "BH1", starttime, filename="IU.ANMO.00.BH1.SACPZ")
irisclient.sacpz("IU", "ANMO", "00", "BH2", starttime, filename="IU.ANMO.00.BH2.SACPZ")
irisclient.sacpz("IU", "ANMO", "00", "BHZ", starttime, filename="IU.ANMO.00.BHZ.SACPZ")
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 1
----> 1 from obspy import UTCDateTime
2 from obspy.clients.fdsn import Client
3 from obspy.clients.iris.client import Client as irisClient
File ~/micromamba/envs/seismo-learn/lib/python3.11/site-packages/obspy/__init__.py:33
30 import sys
32 # don't change order
---> 33 from obspy.core.utcdatetime import UTCDateTime # NOQA
34 from obspy.core.util import _get_version_string
35 __version__ = _get_version_string(abbrev=10)
File ~/micromamba/envs/seismo-learn/lib/python3.11/site-packages/obspy/core/__init__.py:120
12 """
13 obspy.core - Core classes of ObsPy
14 ==================================
(...) 117 .. _NumPy: http://www.numpy.org
118 """
119 # don't change order
--> 120 from obspy.core.utcdatetime import UTCDateTime # NOQA
121 from obspy.core.util.attribdict import AttribDict # NOQA
122 from obspy.core.trace import Stats, Trace # NOQA
File ~/micromamba/envs/seismo-learn/lib/python3.11/site-packages/obspy/core/utcdatetime.py:21
18 import warnings
20 import numpy as np
---> 21 from obspy.core.util.deprecation_helpers import ObsPyDeprecationWarning
24 # based on https://www.myintervals.com/blog/2009/05/20/iso-8601, w/ week 53 fix
25 _ISO8601_REGEX = re.compile(r"""
26 ^
27 ([\+-]?\d{4}(?!\d{2}\b))
(...) 40 $
41 """, re.VERBOSE)
File ~/micromamba/envs/seismo-learn/lib/python3.11/site-packages/obspy/core/util/__init__.py:22
20 # import order matters - NamedTemporaryFile must be one of the first!
21 from obspy.core.util.attribdict import AttribDict
---> 22 from obspy.core.util.base import (ALL_MODULES, DEFAULT_MODULES,
23 NATIVE_BYTEORDER, NETWORK_MODULES,
24 NamedTemporaryFile, _read_from_plugin,
25 create_empty_data_chunk, get_example_file,
26 get_script_dir_name, MATPLOTLIB_VERSION,
27 SCIPY_VERSION, NUMPY_VERSION,
28 CARTOPY_VERSION, CatchAndAssertWarnings)
29 from obspy.core.util.misc import (BAND_CODE, CatchOutput, complexify_string,
30 guess_delta, score_at_percentile,
31 to_int_or_zero, SuppressOutput)
32 from obspy.core.util.obspy_types import (ComplexWithUncertainties, Enum,
33 FloatWithUncertainties)
File ~/micromamba/envs/seismo-learn/lib/python3.11/site-packages/obspy/core/util/base.py:26
23 from pathlib import PurePath
25 import numpy as np
---> 26 import pkg_resources
27 from pkg_resources import get_entry_info, iter_entry_points
29 from obspy.core.util.misc import to_int_or_zero, buffered_load_entry_point
ModuleNotFoundError: No module named 'pkg_resources'
下面以其中一个波形为例，展示如何使用 SAC 零极点文件去除仪器响应：
tr = st[2].copy() # 垂直分量
print(tr)
obspy.io.sac.sacpz.attach_paz() 函数会读取并解析 SAC 零极点文件，并将
解析后的结果附加到 Trace 的 stats.paz 属性中：
from obspy.io.sac.sacpz import attach_paz
attach_paz(tr, paz_file="IU.ANMO.00.BHZ.SACPZ")
Trace.stats.paz 属性是一个词典，包含了 poles、zeros 和 gain 等属性：
from pprint import pprint
pprint(dict(tr.stats.paz))
使用 Trace.simulate 函数去除仪器响应：
tr.simulate(paz_remove=tr.stats.paz)
对于包含多个波形的 Stream 而言，同理。可以对 Stream 中的每个 Trace 进行循环，
为每个 Trace 附加对应的 SAC 零极点文件并去除仪器响应：
for tr in st:
attach_paz(tr, paz_file=f"{tr.id}.SACPZ")
tr.simulate(paz_remove=tr.stats.paz)
也可以在为每个 Trace 附加 SAC 零极点文件后，调用
Stream.simulate 函数批量去除
仪器响应：
for tr in st:
attach_paz(tr, paz_file=f"{tr.id}.SACPZ")
st.simulate(paz_remove="self")
校正到 WWSP 仪器响应#
在使用地核震相研究地核结构时，有时会需要将波形校正到 WWSP
（即 World Wide Standard Seismograph Station short period）仪器响应。
SAC 软件中内置了 WWSP 仪器响应，可以直接使用命令 transfer ... to wwsp 将
波形校正到 WWSP 仪器。而 ObsPy 中没有内置 WWSP 仪器响应的信息，因而需要自行定义
WWSP 仪器响应信息并做校正。
准备示例所需数据：
from obspy import UTCDateTime
from obspy.clients.fdsn import Client
starttime=UTCDateTime("2022-09-22T06:18:00")
client = Client("IRIS")
st = client.get_waveforms("IU", "ANMO", "00", "BH?", starttime, starttime + 720)
定义 WWSP 仪器响应，需要零极点和增益信息。这些信息可以从 SAC 源代码
sac/src/icm/wwsp.c 中获得。定义如下：
paz_wwsp = {
"poles": [
-5.0136607 + 6.4615109j,
-5.0136607 - 6.4615109j,
-8.2981509 + 0.0j,
-8.6940765 - 7.1968661j,
-8.6940765 + 7.1968661j,
],
"zeros": [0j, 0j, 0j],
"gain": 397.54767,
"sensitivity": 1.0,
}
定义好 WWSP 仪器响应后，即可使用 Stream.simulate
函数去除仪器响应：
st.simulate(paz_simulate=paz_wwsp)
该方法也同样适用于其他 SAC 内置仪器响应。