本文简单介绍如何利用 ObsPy 下载理论地震图。

# 可将本段代码全部复制到 obs_test.py 文件中, 在终端进行操作
#
# $ source activate obspy   # 在当前终端激活 obspy 运行环境 
# $ python obs_test.py        # 执行脚本
# $ source deactivate        # 退出 obspy 运行环境
#
# 第一步, 设置绘图方式
import matplotlib.pyplot as plt
plt.style.use("ggplot")

# 第二步, 利用 obspy 的 fdsn 客户端连接 IRIS 数据中心
import obspy
from obspy.fdsn import Client
# 在新版本中应该使用 obspy.clients.fdsn 代替 obspy.fdsn
c = Client("IRIS")

# 第三步, 获取地震信息(可选),见图 获取地震信息
cat = c.get_events(starttime=obspy.UTCDateTime(2010, 5, 24, 16),
                   endtime=obspy.UTCDateTime(2010, 5, 24, 17),
                   minmagnitude=6.4)
print("Event depth in km:", cat[0].origins[0].depth / 1000.0)
cat.plot(projection="ortho");

获取地震信息

# 第四步, 获取地震波形,见图 获取地震波形
st = c.get_waveforms(network="IU", station="ANMO", location="00", channel="BHZ",
                     starttime=obspy.UTCDateTime(2010, 5, 24, 16, 18, 28),
                     endtime=obspy.UTCDateTime(2010, 5, 24, 17, 18, 28),
                     attach_response=True)
st.plot()

获取地震波形

# 第五步, 滤波, 去除仪器响应,见图 去除仪器响应
st.remove_response(output="DISP", pre_filt=(0.005, 0.01, 4, 8))
st.plot()

去除仪器响应

# 第六步, 下载理论地震图,见图 下载理论地震图
st_synth = obspy.read("http://service.iris.edu/irisws/syngine/1/query?"
                      "network=IU&station=ANMO&components=Z&"
                      "eventid=GCMT:C201005241618A&dt=0.05&"
                      "units=displacement&"
                      "model=iasp91_2s&label=Tutorial&format=miniseed")
st_synth.plot()

下载理论地震图

# 第七步, 计算理论到时
from obspy.taup import TauPyModel

m = TauPyModel("IASP91")
arr = m.get_travel_times(
    source_depth_in_km=582.1, distance_in_degree=54.16,
    phase_list=["P"])
print(arr)

# 第八步, 对比P波到时附近观测波形和理论波形,见图 观测波形和理论波形对比图
# Combine into a single Stream object.
st_all = st.copy() + st_synth.copy()

# Filter to a period band from 1 to 10 seconds to highlight
# the P phase and force the same band limit on both.
st_all.detrend("demean")
st_all.taper(0.05)
st_all.filter("bandpass", freqmin=0.1, freqmax=1.0, corners=3)

# Cut window around the P arrival
starttime = st_synth[0].stats.starttime + 470
endtime = st_synth[0].stats.starttime + 550
st_all.trim(starttime, endtime)

st_all.plot()

观测波形和理论波形对比图

参考文献