Obspy / Scipy 带通滤波器 feat. SAC

2020.09.14

Bandpass 带通滤波器是信号处理是常见的滤波工具。本文介绍 Python 中分别使用 scipy 和 obspy 实现的零相位偏移的的两种带通滤波器实现，并类比 SAC 中的滤波操作的结果。

测试信号生成

使用 CPS 合成理论地震记录并保存为sac文件：

文件 create_seismology.sh

#!/bin/sh
cat > dfile << EOF
0.5,0.01,512,0,0
EOF

cat > soil.mod << EOF
MODEL
simple soil model
ISOTROPIC
KGS
FLAT EARTH
1-D
CONSTANT VELOCITY
LINE08
LINE09
LINE10
LINE11
HR      VP      VS  RHO QP  QS  ETAP ETAS FREFP FREFS
0.010000 1.417382 0.250000 1.898545  0.0 0.0 0.0 0.0 1.0 1.0
0.010000 1.584648 0.350000 1.952237  0.0 0.0 0.0 0.0 1.0 1.0
0.010000 1.740763 0.450000 1.998638  0.0 0.0 0.0 0.0 1.0 1.0
00. 1.815069 0.500000 2.019633  0.0 0.0 0.0 0.0 1.0 1.0
EOF

sprep96 -M soil.mod -d dfile  -R -L -NMOD 3
sdisp96 >sdisp96.out
scomb96 -R -FREQ -FMIN 0.1 -FMAX 100 -CMIN 0.1 -CMAX 3
scomb96 -L -FREQ -FMIN 0.1 -FMAX 100 -CMIN 0.1 -CMAX 3
sregn96 -HS 0 -HR 0 -T >sregn96.out
slegn96 -HS 0 -HR 0 -T >slegn96.out
spulse96 -d dfile -EX -D -i | fmech96 -fx 1 -fy 0 -fz 0 -A 45 -B 225 -ROT| f96tosac

生成了三个 sac 文件：B00101Z00.sac、B00102R00.sac、 B00103T00.sac

以B00101Z00.sac文件进行测试：进入 sac 命令行，读入并使用4阶带通滤波器分别进行一次和正反两次滤波后保存：

r B00101Z00.sac
bp n 4 p 1 c 3 10
w n4p1.sac
r B00101Z00.sac
bp n 4 p 2 c 3 10
w n4p2.sac

其中p参数的取值可以参考：

p 取1时，对波形做一次带通滤波，由于滤波器存在相位延迟，因而导致波形的峰值出现了时间延迟，因而会影响到震相的最大峰值的拾取，但对震相的初至却没有影响。 p 取2时，对波形做正反两次带通滤波，此时不存在相位延迟，因而不会影响到最大峰值的拾取，但震相的初至则存在时间上的提前。 ——SAC参考手册：滤波

Obspy实现

from obspy.signal.filter import bandpass
import obspy
import matplotlib.pyplot as plt
import numpy as np

tvec = np.linspace(0, 0.01, 512)
st = obspy.read('B00101Z00.sac')

## METHOD 1
p2 = st.copy()
p2[0].filter("bandpass", freqmin=3, freqmax=10, corners=4, zerophase=True)

## METHOD 2
data_p1 = bandpass(st[0].data, freqmin=3, freqmax=10, df=100, corners=4, zerophase=True)
n4p2 = obspy.read('n4p2.sac')

plt.plot(tvec, n4p2[0].data, label='SAC', linewidth=5)
plt.plot(tvec, data_p1, label='obspy METHOD 1')
plt.plot(tvec, p2[0].data, label='obspy METHOD 2', c='w')
plt.legend()
plt.show()

可以看到输出绘图完全一致。通过zerophase=True参数指定进行了正反两次带通滤波，从而消除了相位延迟。

以上代码分别展示了两种带通滤波器的使用：

METHOD 1 为调用obspy.core.trace.Trace.filter
METHOD 2 为调用obspy.signal.filter.bandpass

需要注意的是METHOD 2中的采样率df为必要参数，METHOD 1则无需填写。此外，METHOD 1 的本质也是调用METHOD 2。

Scipy实现

import obspy
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal

tvec = np.linspace(0, 0.01, 512)
st = obspy.read('B00101Z00.sac')

band = np.array([3, 10])
Wn = band * 2 / 100
b, a = signal.butter(4, Wn, 'bandpass') 

## METHOD 1

pp2 = signal.filtfilt(b, a, st[0].data)

## METHOD 2
p1 = signal.lfilter(b, a, st[0].data) 
p2 = signal.lfilter(b, a, p1[::-1])
p2 = p2[::-1]

n4p2 = obspy.read('n4p2.sac')

plt.plot(tvec, n4p2[0].data, label='SAC', linewidth=5)
plt.plot(tvec, pp2, label='scipy METHOD 1')
plt.plot(tvec, p2, label='scipy METHOD 2', c='w')
plt.legend()
plt.show()