基于Python库ObsPy的STA/LTA的事件检测和挑选

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')