obspy初步学习

秃然的凉意

已于 2024-10-09 22:40:18 修改

阅读量1.1k

点赞数 11

文章标签：笔记

于 2024-10-09 22:39:38 首次发布

本文链接：https://blog.youkuaiyun.com/qq_74825666/article/details/142798328

版权

obspy初步学习

学习目标：

读取数据（mseed和SAC）并查看信息
如何对简单的波形图像进行绘制以及matplotlib连用
如何进行滤波等处理

学习内容：

读取数据并查看信息：

调用obspy.read方法进行对数据的读取：

import obspy
st = obspy.read("http://examples.obspy.org/RJOB_061005_072159.ehz.new")

其中读出的数据处于一个数据流（stream）中，stream可以容纳多个地震数据。为了保证方便批量处理，所以类似于滤波等方法全部位于stream的类中。例如：st.filter("lowpass", freq=0.1, corners=2) #这里的例子我们采用低通;截至频率0.1hz;二阶滤波器

2.查看信息：

我们这里的信息主要分两大类：1. 是地震文件的表头信息：包括震源位置，震源深度，台站位置，采样频率、周期等表头信息;2. 是地震文件正文部分信息：即地震数据。查看的方法是先从数据流中导出我们要查看的某一个数据，然后进行单独查看。

导出数据：

st = obspy.read("http://examples.obspy.org/RJOB_061005_072159.ehz.new")
#这里我们在再一个数据加入数据流，模拟我们真实的多个台站同时处理的情况
st += obspy.read(r"./RJOB_061005_072159.ehz.new")
tr = st[0]
print(tr.stats)#获取表头信息
print(tr.data)#获取存储数据信息
print(tr)#获取简要的表头信息（包括starttime-endtime；采样率等信息）
print(tr.times("matplotlib"))#可以导出各个采样点的时间t
#==注意==：如果不加"matplotlib"就会导致输出的只是一个object类
#无法显示其具体的时间t

调用数据：

tr.stats.(starttime/npts/sampling_rate/……)
#只要tr.stats出现的都属于是stats下面可调用的参量

对图像进行绘图：

分为两大分支：1.利用obspy内部集成的plot直接绘制；2.将数据导出，采用matplotlib进行绘制。

st.plot() #对于地震数据流进行画图（由于类似SAC里面绘图，我们可能要同时绘制多个台站接受到的数据，所以其数据流类型的内部集成了一个plot绘图方便绘制多台站）：

dt = st[0].stats.starttime
#这里我们也可以直接st.plot()直接绘制无修饰图像
#这里我们将绘图进行一定的调控，使得绘制图像为红色
#横轴坐标进行旋转以及时分的格式输出
#绘制的时间是开始后推3600s到3600+120s结束的两分钟地震图像
st.plot(color="red",tick_rotation=5,tick_format='%I:%M %p'
,starttime=dt + 60*60,endtime=dt + 60*60 + 120)

可以根据设置starttime和endtime来决定我们的图像画的是哪个时间区间的图像。可以通过设置tick_format设置横轴的变量格式。

#dayplot格式图，将数据按一定的时间间隔平铺到图像中。
#type选择输出图片为dayplot图，interval表示图片的时间间隔为60min
#right_vertical_labels表示有标轴为False不显示
#vertical_scaling_range表示幅值增加5000倍
#one_tick_per_line表示每一个line一个纵轴信息
#color表示设置的颜色变化按照四个一组变化
#show_y_UTC_label关闭纵轴UTC_label的时间戳
#events={’min_magnitude’: 6.5}识别6.5级以上的地震
st.plot(type="dayplot", interval=60, right_vertical_labels=False,
... vertical_scaling_range=5e3, one_tick_per_line=True,
... color=[’k’, ’r’, ’b’, ’g’], show_y_UTC_label=False,
... events={’min_magnitude’: 6.5})

matplotlib绘图：

import obspy
import matplotlib.pyplot as plt
st = obspy.read("http://examples.obspy.org/RJOB_061005_072159.ehz.new")
t = st[0].time("matplotlib")#数据转为matplotlib可识别格式
am = st[0].data
fig = plt.figure(1)
ax = fig.add_subplot(1,1,1)
ax.plot(t,am,"b-")
#自动调整date格式，从原先的天转变为时分秒格式
ax.xaxis_date()
#自动调整横轴的时间标准格式，auto表示自动
#fmt表示format。autofmt_xdate自动横轴时间排版调整
fig.autofmt_xdate()

上面的tr.times(“matplotlib”)时间获取是计算机诞生的一个总体时间，而我们采用np.arange(0,tr.stats.npts/tr.stats.sampling_rate,
tr.stats.delta)，实际上是计算采样总时间，将其轴归到0处所绘制的图像。

滤波等预处理：

滤波：

可以对于数据流进行滤波，也可以对单独拿出来的trace台站数据进行处理。（滤波类型：lowpass，bandpass，highpass）（注意：带通的核心就是freqmin和freqmax的规定[这也是和高通低通的区别所在]）

#==下面展示的是低通滤波==
from obspy import read 
# 读取数据 st = read("path/to/your/seismogram.mseed") 
# 获取第一个 Trace 
tr = st[0] 
# 创建 Trace 的副本，以避免修改原始数据 
tr_lowpass = tr.copy() 
# 应用低通滤波器 
tr_lowpass.filter('lowpass', freq=1.0, corners=2, zerophase=True)

#==下面展示的是带通滤波==
from obspy import read 
# 创建带通滤波的副本 
tr_bandpass = tr.copy() 
# 应用带通滤波器 
#==注意带通的核心就是freqmin和freqmax的规定==
tr_bandpass.filter('bandpass', freqmin=0.1, freqmax=1.0, corners=2, zerophase=True)

#==下面展示的是高通滤波==
from obspy import read 
# 创建另一个 Trace 的副本 
tr_highpass = tr.copy() 
# 应用高通滤波器 
tr_highpass.filter('highpass', freq=0.1, corners=2, zerophase=True)

下采样（降频处理）：

地震原始数据为了保证更多的细节，往往采用较高的采样频率。但这样会导致我们数据处理效率低下（因为数据处理量大）。在保证原数据不丢失的情况下，降低采样率是我们所追求的。由于奈奎斯特采样频率，至少要保证2f才不会堆叠。所以我们下采样之前，要先滤波，取出高频噪声，保证下采样可以顺利进行。

import matplotlib.pyplot as plt
import obspy
import numpy as np
st = obspy.read("https://examples.obspy.org/RJOB_061005_072159.ehz.new")
tr = st[0]
tr0 = tr.copy()	#tr1只进行滤波，不进行下采样
#==进行滤波和下采样==：
#factor表示间隔为原先的4倍;strict_length表示严格按照factor长度
#<注释>：当为True时截断在最后间隔，False会自动进行补充
#例如我们要以256为间隔，对1000个数据点下采样。
#True为744截断，False会通过补零变成1024.
#这个decimate()方法是默认会先进行低通滤波
#f截默认为0.4*(tr.stats.sampling)/4，是因为奈奎斯特采样频率要保证其不过半
tr.decimate(factor=4,strict_length=False)
#这是完整格式，no_filter=False表示优先进行滤波。
#tr.decimate(factor=4,strict_length=False,no_filter=False)
t = np.arange(0,tr.stats.npts/tr.stats.sampling_rate,
tr.stats.delta)
#只进行滤波
tr0.filter("lowpass",freq=0.4 * tr.stats.sampling_rate / 4.0)
t0 = np.arange(0,tr0.stats.npts/tr0.stats.sampling_rate,
tr0.stats.delta)
plt.plot(t0,st[0].data,"k",label="raw")
plt.plot(t0,tr0.data,"r",label="lowpass")
plt.plot(t,tr.data,"b",label="lowpass/new_sampling")
plt.xlabel('Time [s]')
#只显示82s-83.5s图像
plt.xlim(82, 83.5)
plt.legend()
plt.show()

实际中我们发现数据发生了相移，可以采用tr0.filter(“lowpass”,freq=0.4 * tr.stats.sampling_rate / 4.0,zerophase=True)去除掉滤波时的相移。但是降采样时的相移还没办法解决。下面我们采用直接利用<步长>降采样获取：

import matplotlib.pyplot as plt
import obspy
import numpy as np
st = obspy.read("https://examples.obspy.org/RJOB_061005_072159.ehz.new")
tr = st[0]
#建立新的数据
tr_new = tr.copy()
#回去新的采样率，导出步长，采取strict_length=True的选项
origin_sampling = tr.stats.sampling_rate
new_sampling = origin_sampling/4
#预处理
tr.filter("lowpass",freq=0.4*new_sampling,zerophase=True)
#步长：
factor = int(origin_sampling/new_sampling)
#降采样(这里该数据和采样率后，采样点数也降低了)：
tr_new.data = tr.data[::factor]
tr_new.stats.sampling_rate = new_sampling
t_new = np.arange(0,tr_new.stats.npts/tr_new.stats.sampling_rate,tr_new.stats.delta)
#绘图：
plt.plot(t_new,tr_new.data,"b",label="lowpass/new_sampling")
plt.xlabel('Time [s]')
#只显示82s-83.5s图像
plt.xlim(82, 83.5)
plt.legend()
plt.show()

SAC中lh o设置起始时间：

思路核心：先将数据读入数据流stream中;然后按时间排序获得事件最早时间;可以利用matplotlib中subplot获得并列图像

import numpy as np
import matplotlib.pyplot as plt
import obspy
st = obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE")
st += obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE.1")
st += obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE.2")
#st要用于时间先后排序从而确定最早的时间，所以会自带排序工具  
st.sort(['starttime'])  
#这里采用timestamp是为了获取一个方便画图的1970年起的单位为s的时间  
#st[0].stats.starttime获得的时间是UTCDateTime 类型的对象  
#方便进行时间上的运算，s的时间只能方便数值计算和绘图。 
dt = st[-1].stats.starttime.timestamp  
num = len(st)#获取要画的子图数量  
fig,axe = plt.subplots(4,1,sharex=True)  
fig.suptitle("multi-wave")  
for i in range(num):  
    t = np.linspace(st[i].stats.starttime.timestamp - dt,  
                    st[i].stats.endtime.timestamp - dt,  
                    st[i].stats.npts)  
    axe[i].plot(t, st[i].data,linewidth=0.5)  
    axe[i].grid(True)  # 确保这一行在循环内部，并且缩进正确  
axe[2].set_xticks([0,1000,2000,3000,4000])
plt.xlim([0,5000])#就可以实现完全的地震o归零

#在上述的基础上，我们把三个一样的数据合并
#==注意==：数据是一致的但是由于仪器等因素，可能会缺失;
#这时候合并保证了获得连续的数据。
st.merge(method=1)  
#此时合并到了第一张图上:
#(==注意==：此时数据长度由于合并发生了变化，为了顺利绘图，要重新获得t_new)  
t_new = np.linspace(st[0].stats.starttime.timestamp - dt,  
                st[0].stats.endtime.timestamp - dt,  
                st[0].stats.npts)  
axe[3].plot(t_new,st[0].data,"r",linewidth=0.5,label="merge_picture")  
axe[3].legend()  
plt.show()

这里的merge(method=0，1，2)

method=0：简单地连接相邻的 Trace 片段，即使它们之间有间隙。
method=1：仅当相邻的 Trace 片段之间的差距小于等于 1 秒时才合并。
method=2：仅当相邻的 Trace 片段之间的差距小于等于 0.05 秒时才合并。

去除台站响应：

思路核心：1.给数据流中的每一个文件导入其对应的台站响应;2.利用simulate方法除去台站响应参数

import matplotlib.pyplt as plt
import obspy
#AttribDict用于给trace加入台站信息
from obspy.core.util import AttribDict
st = obspy.read("https://examples.obspy.org/agfa.mseed")
#添加台站参数paz,只能装字典类型的
st[0].stats.paz = AttribDict({'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],  
                              'zeros': [0j, 0j],  
                              'sensitivity': 205479446.68601453,  
                              'gain': 1.0})  
#poles极点，zeros零点，sensitivity灵敏度，gain仪器放大系数
st[0].stats.coordinates = AttribDict({'latitude': 48.108589,  
                                      'elevation': 0.450000,  
                                      'longitude': 11.582967})
#利用simulate方法去除仪器响应
st.simulate(paz_remove='self')

绘制图像包络线：

思路核心是使用obspy.signal.filter的子库。这里我们模拟一个真实的地震包络线处理：1.缺失数据填充;2.去除台站响应;3.滤波;4.绘制包络线。

import numpy as np
import matplotlib.pyplot as plt
import obspy
from obspy.signal.filter import envelope
st = obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE")
st += obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE.1")
st += obspy.read("https://examples.obspy.org/dis.G.SCZ.__.BHE.2")
#缺失数据合并  
st.merge(method=1)  
tr = st[0]  
data1 = tr.data  
#去除台站响应,导入相应参数  
tr.stats.paz = AttribDict({'poles': [(-0.03736 - 0.03617j), (-0.03736 + 0.03617j)],  
                              'zeros': [0j, 0j],  
                              'sensitivity': 205479446.68601453,  
                              'gain': 1.0})  
tr.simulate(paz_remove='self')  
data2 = tr.data  
#滤波  
tr.filter('lowpass', freq=0.5, corners=2, zerophase=True)  
data3 = tr.data  
#包络线只要是提供一组数据就可以绘制。  
envelope_tr_data = envelope(data3)  
#绘图  
#这里为了获得每次处理后的结果我们选择绘制4张图  
t = np.arange(0,tr.stats.npts/tr.stats.sampling_rate,tr.stats.delta)  
fig,axe = plt.subplots(3,1,sharex=True)#共x轴  
data_processing = [data1,data2,data3]  
label = ['Merge data', 'remove station response', 'Filter data/envelope']  
for i in range(3):  
    axe[i].plot(t,data_processing[i],linewidth=0.5,label=label[i])  
    axe[i].grid(True)  
    axe[i].legend()  
    axe[i].set_ylabel('Amplitude')  
axe[2].plot(t,-1*envelope_tr_data,"k:")  
axe[2].set_xlabel('Time [s]')  
plt.xlim(800, 1200)  
plt.show()