obspy去polezero类型仪器仪响应

发表于 更新于 分类于 阅读次数： 本文字数： 1.1k 阅读时长 ≈ 1 分钟

  如果有polezero文件，那可以用attach_paz来读取。然后去仪器响应的脚本如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import numpy as np
import obspy
import matplotlib.pyplot as plt
from obspy import read
from obspy.io.sac import attach_paz
from obspy.signal.invsim import corn_freq_2_paz
from pathlib import Path
fpath="/home/junxie/work/Snet/data/2016_08_16_01_00"
p=Path(fpath)
i=0
for file in p.rglob('*VX*.SAC'): #循环读取文件
    st=read(file, debug_headers=True)
    tr=st[0].copy()
    aa=str(file).split("/")
    pz=str(file)+"_PZ" # 这就是paz文件
    attach_paz(tr,pz)  # 贴到tr中
    paz_1hz = corn_freq_2_paz(1.0, damp=0.707)  # 1Hz instrument
    paz_1hz['sensitivity'] = 1.0
    st.simulate(paz_simulate=paz_1hz) #去仪器响应

  如果没有的话找到零点、极点和放大系数，赋值到paz里面，例如：

1
2
3
4
5
6
7
8
9
10
11
from obspy import read
from obspy.signal.invsim import corn_freq_2_paz
st = read()
paz_sts2 = {'poles': [-0.037004+0.037016j, -0.037004-0.037016j,
                      -251.33+0j,
                      -131.04-467.29j, -131.04+467.29j],
            'zeros': [0j, 0j],
            'gain': 60077000.0,
            'sensitivity': 2516778400.0}
paz_1hz = corn_freq_2_paz(1.0, damp=0.707)
st.simulate(paz_remove=paz_sts2, paz_simulate=paz_1hz)