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)转载代码或有其他用途请注明出处。