Computing cross-correlation and spectrogram of two seismic traces (codes included)

LibraryPKU 2022-09-11 发布于北京

展开全文

我想使用互相关计算两个地震轨迹之间的关系。首先，我使用spectrogramObspy 中的方法计算每个轨迹的频谱图。然后，我对两条地震轨迹进行互相关并绘制其时域互相关nf函数和频谱图。

有关如何计算互相关的详细信息，请访问我之前的帖子：

计算并绘制每条迹线的频谱图

我有一个all_stream_HLZ_20071215_080316.mseed包含多个跟踪的 mseed 文件。我在这项研究中使用了前两条迹线。read为了绘制时间序列，我首先使用函数读取数据nfrom Obspy，然后绘制Trace对象可用的“数据”和“时间”方法。

类似帖子

如何在基础上绘制州的边界美国地图

在 bash、c/c++ 和 python 中读取 yaml 输入文件

什么是数据结构中的链表？

基于l-bfgs方法的数值优化

使用 python 创建投影仪幻灯片

from obspy import readfrom obspy import read, Trace, UTCDateTimeimport numpy as npimport pandas as pdimport noise将 matplotlib.pyplot 导入为 pltplt.rcParams['figure.figsize'] = [10,6]plt.rcParams.update({'font.size': 18})plt.style.use('seaborn')filenameZ = ' all_stream_HLZ_20071215_080316.mseed'fignameTrace = 'spectrogram_layout1.png'figxcorr = 'spectrogram_layout2.png'figxcorr2 = 'spectrogram_layout3.png'st = read(filenameZ)tr1 = st[0]tr2 = st[1]sta1 = tr1.stats[' station']fig, axx = plt.subplots(2,2, sharex=True, sharey='row')axx[0, 0].plot(tr1.times(), tr1.data, 'k-', 线宽=0.2, label=tr1.stats['station'])axx[0, 1].plot(tr2.times(), tr2.data, 'k-', linewidth=0.2, label=tr2.stats['station '])axx[0, 0].set_title(f'Trace 1')axx[0, 1].set_title(f'Trace 2')axx[0, 0].set_ylabel('Amplitude')tr1.spectrogram(日志=真，wlen=50，显示=假，轴=axx[1, 0]，c映射='jet', samp_rate=tr1.stats.sampling_rate)tr2.spectrogram(log=True, wlen=50,show=False, axes=axx[1, 1], c map ='jet', samp_rate=tr2.stats .sampling_rate)axx[1, 0].set_title('Spectrogram 1')axx[1, 1].set_title('Spectrogram 2')axx[1, 0].set_xlabel('time (s)')axx[1 , 1].set_xlabel('time (s)')axx[1, 0].set_ylabel('Amplitude')plt.tight_layout()## put legendfor col in axx[0,:]:ll = col.legend( loc=1)plt.setp(ll.get_texts(), color='red') #color legendplt.savefig(fignameTrace, bbox_inches='tight', dpi=300)plt.close('all')

地震轨迹图及其对应的频谱图

使用 Pandas 库计算互相关

为了计算互相关，我使用该crosscorr函数。读者可以参考这篇文章中的这个功能。计算互相关的步骤也与上一篇非常相似。

但是，我使用spectrogramObspy 的方法获得了频谱图。该方法对地震数据进行了很好的优化。Python 中还有其他几个可用的频谱图函数，它们中的大多数都以相同的方式工作。本文获得的频谱图使用 fft 为 50 秒 ( wlen) 的窗口长度，并以对数刻度输出频率。

### Cross Correlation using Pandasseries1, series2 = pd.Series(tr1.data), pd.Series(tr2.data)window = 1maxlag = 500lags = np.arange(-(maxlag), (maxlag), window) # contrainedrs = np.nan_to_num([crosscorr(series1, series2, lag) for lag in lags])traceCCF2 = Trace()traceCCF2.data = rstraceCCF2.times(reftime=tr1.stats.starttime)fig, axx = plt.subplots(2 ,1)axx[0].plot(lags, rs, 'k', linewidth=0.5)axx[0].set_title('Cross Correlation')traceCCF2.spectrogram(log=True, wlen=50,show=False, axes=axx[1], c map ='jet')axx[1].set_title('Spectrogram')plt.tight_layout()plt.savefig(figxcorr2, bbox_inches='tight', dpi=300)plt.close( '全部'）