IGU-BD3C-5地震仪SAC原始数据频率分析

IGU-BD3C-5地震仪是教学中常用的地震仪，其输出的文件名形如590001488.0001.2023.07.18.07.05.52.000.Z.sac。

本文针对该仪器（以及该输出文件名格式），结合shell语言、SAC（结合仪器厂家给出的pz1文件）、Python3，在前人基础上完善该程序：

#! /bin/bash
#! /usr/bin/python3
#Optimized by ljw in 2023/8/9, China University of Geoscience(Wuhan)

sac << EOF
r *.?.sac
rglitches
rtr
rmean
taper 0.2
transfer from polezero subtype pz1 to none freq 0.1 0.2 100 200
mul 1.0e9
w over
q
EOF
rm -rf Z
rm -rf E
rm -rf N
for file in *.sac
do
rm -rf name.txt
echo $file >> name.txt
python3 << EOF
import struct
import pylab
import numpy as np
from scipy.fftpack import fft,ifft 
class sacfile_wave:
    def read(self, sFile):
        f = open(sFile, 'rb')
        hdrBin = f.read(632)
        sfmt = 'f' * 70 + 'I ' * 40 + '8s ' * 22 + '16s';
        hdrFmt = struct.Struct(sfmt)
        self.m_header = hdrFmt.unpack(hdrBin)
        self.dt = self.m_header[0]
        npts = int(self.m_header[79])
        fmt_data = 'f' * npts
        dataFmt = struct.Struct(fmt_data)
        dataBin = f.read(4 * npts)
        f.close()
        self.m_data = dataFmt.unpack(dataBin)
    def draw(self, sImageFile):
        npts = len(self.m_data)
        xd = range(1, npts + 1)
        dt = self.dt
        NFFT = 250
        Fs = int(1.0 / dt)
        x0 = np.arange(0, npts/Fs, np.round(dt,2))
        pylab.subplot(311)
        pylab.plot(x0, self.m_data, linewidth=.6)
        # pylab.xlabel('Time / second')
        # pylab.xlabel('时间 / 秒', fontproperties = 'FangSong', fontsize = 20)
        pylab.subplot(312)
        Pxx, freqs, bins, imm = pylab.specgram(self.m_data, NFFT=NFFT, Fs=Fs,
                                                noverlap=NFFT * 19 / 20,cmap=pylab.get_cmap('jet'))
        pylab.subplot(313)
        n = len(self.m_data)
        k = np.arange(n)
        T = n / Fs
        frq = k / T
        frq1 = frq[range(int(n / 2))]
        YY = np.fft.fft(self.m_data)
        Y = np.fft.fft(self.m_data) / n
        Y1 = Y[range(int(n / 2))]
        Y1[0] = Y1[1]
        pylab.plot(frq1, abs(Y1))
        pylab.savefig(sImageFile, dpi = 1200)
    def exportAsc(self, sAscFile):
        f2 = open(sAscFile, "wt")
        sdataAsc = [str(x) for x in self.m_data]
        sDataAsc = '\n'.join(sdataAsc)
        f2.writelines(sDataAsc)
        f2.close()
if __name__ == "__main__":
    with open("name.txt","r",encoding='utf-8') as f:
	    z1 = f.read()
    print(z1)
    z2 = z1[5:9]
    z3 = z1[15:19]
    z4 = z1[20:22]
    z5 = z1[23:25]
    z6 = z1[39:40]
    sacfile = z1.strip()
    sac = sacfile_wave()
    sac.read(sacfile)
    sac.draw(z2+'台站'+z3+'年'+z4+'月'+z5+'日'+z6+'方向.png')
EOF
rm -rf name.txt
done
mkdir E
mkdir N
mkdir Z
mv *E方向.png E
mv *N方向.png N
mv *Z方向.png Z

说明：

（1）该代码的python3部分主要来自博客园某位用户，现在搜不到原文（respect！！！！），本人稍作修整改了程序接口使之更为灵活；SAC部分去除仪器响应需要结合厂家的pz1文件，否则SAC部分将会报错，无法去除仪器响应；

（2）该代码运行在Linux环境下，在Linux(Ubuntu23.04)中正确安装python3+SAC以及相关依赖库；

（3）IGU-BD3C-5地震仪默认输出的标准命名的地震数据(形如
590001488.0001.2023.07.18.07.05.52.000.Z.sac)作为数据输入；

（4）将（批量的）地震数据放在名为sac的文件夹中，将代码放在sac文件夹中，sac文件夹中无其他无关数据；

（5）使用bash运行该代码；

（6）该代码运行较慢，2GB的SAC数据需要运行约1小时；

（7）运行结果节选（从上至下依次为波形图、时频图、频谱图）：

（8）在运行前请自主备份数据，数据丢失本人概不负责；

（9）转载代码或有其他用途请注明出处。