1、Obspy安装
官方建议利用Anaconda安装 ObsPy,国内用户可能需要更换镜像源。
# 创建obspy运行环境
$ conda create -n obspy python=3.8
# 查看所有环境名称
$ conda info -e
# 激活obspy环境
$ conda activate obspy
# 安装obspy预编译包
(obspy) $ conda install obspy -c conda-forge
2、源代码
2-1、导入库
import numpy as np
import matplotlib.pyplot as plt
import obspy
from obspy.signal.trigger import classic_sta_lta, plot_trigger, trigger_onset
2-2、STA / LTA参数设置
windows = 5 # 截取的时窗长度
THRESH_ON = 8 # STA/LTA触发条件(通常设置为长、短时窗比值的0.8倍)
THRESH_OFF = 5 # STA/LTA关闭条件
SHORT_WINDOW = 0.1 # 短窗长度
LONG_WINDOW = 1 # 长窗长度
2-3、读入数据信息
read_input = obspy.read('./seis.sac')[0]
sampling_interval = read_input.stats.sac.delta
sampling_frequency = int(1/sampling_interval)
total_point = read_input.stats.npts # 输入数据的总点数
total_time = total_point*sampling_interval # 输入数据的总时间长度
total_t = np.arange(0, total_time, sampling_interval)
2-4、STA / LTA触发计算及画图
# obspy.signal.trigger.classic_sta_lta 计算标准STA/LTA
cft = classic_sta_lta(read_input.data, SHORT_WINDOW/sampling_interval, LONG_WINDOW/sampling_interval)
# 画图
plt.figure(figsize=(12, 6))
plt.subplot(211)
plt.plot(total_t, read_input.data, 'k', linewidth=1)
plt.ylabel('Amplitude')
plt.subplot(212)
plt.plot(total_t, cft, 'y', linewidth=1)
# 添加触发和终止信息
cft1 = trigger_onset(cft, THRESH_ON, THRESH_OFF) # STA/LTA触发器
for i in range(len(cft1)):
plt.vlines(cft1[i, 0]*sampling_interval, min(cft), max(cft), colors='r', linestyles='dashed')
plt.vlines(cft1[i, 1]*sampling_interval, min(cft), max(cft), colors='b', linestyles='dashed')
plt.xlabel('Time (s)')
plt.ylabel('STA/LTA')
plt.savefig('sla_lta.png', dpi=1000, bbox_inches='tight')
2-5、计算一个长信号中触发及终止的次数
for i in range(len(cft1)):
tmp = (cft1[i, 1] - cft1[i, 0])*sampling_interval # 每个事件的时间长度
# 取事件窗口前1/3的时间和事件窗口后2/3的时间
tmp = int((windows - tmp)/3)
# 裁减数据
input_copy = read_input.copy()
input_copy.trim(read_input.stats.starttime + int(cft1[i, 0]*sampling_interval) - tmp, \
read_input.stats.starttime + int(cft1[i, 0]*sampling_interval) - tmp + windows)
# if (len(input_copy.data) % sampling_frequency) != 0:
# input_copy.data = input_copy.data[: -1]
# if len(input_copy.data) != sampling_frequency*windows:
# temp_endtime = input_copy.stats.endtime
# input_copy = input.copy()
# input_copy.trim(temp_endtime - windows, temp_endtime)
# print(len(input_copy.data))
# input_copy.data = input_copy.data[: -1]
input_copy.write('./seis_new_' + str(i+1) + '.sac', format='SAC')