import numpy as np import matplotlib.pyplot as plt from obspy import read # 读取面波数据并画图。 st = read('MASW_DATA/Sample_Data/*.SAC') dt = st[0].stats.delta data = [] scale = 0.05 dx = 2 plt.figure(figsize=(8, 6)) for i, tr in enumerate(st): d = tr.data data.append(d) t = np.arange(len(d)) * dt plt.plot(t, d*scale+(i+1)*dx, lw=1, color='b') plt.xlabel('Time (s)') plt.ylabel('Offset (m)') plt.tight_layout() plt.savefig('Surface_wave.png') plt.show() # 二维FFT。 d = np.array(data) n = len(d[0]) # m为空间方向的采样点数，m增大可以让FK谱光滑一点，以达到插值效果。 m = len(d[:, 0]) * 5 D = np.zeros((m, n)) D[:len(d[:, 0])] = d # 时间采样率。 fs = 1 / dt # 空间采样率 xs = 1 / dx # 频率 (赫兹)。 f = np.arange(-n//2, n//2) * fs / (n-1) # 波数 (每米)。 k = 2 * np.pi * np.arange(-m//2, m//2) * xs / (m-1) # 二维FFT。 fk = np.fft.fft2(D) # 作图。 pmin = -10 P = abs(np.fft.fftshift(fk)); P /= P.max(); P = 10 * np.log10(P) P2 = abs(fk); P2 /= P2.max(); P2 = 10 * np.log10(P2) plt.figure(figsize=(11, 8)) plt.subplot(221) plt.pcolormesh(f, k, P2, cmap='magma', vmin=pmin, vmax=0) plt.xlabel('Frequency (s$^{-1}$)') plt.ylabel('Wave number (2$\pi$m$^{-1}$)') plt.subplot(222) plt.pcolormesh(f, k, P, cmap='magma', vmin=pmin, vmax=0) plt.plot([f[n//2], f[-1], f[-1], f[n//2], f[n//2]], [k[0], k[0], k[m//2], k[m//2], k[0]], lw=2, ls='--', color='r') plt.xlabel('Frequency (s$^{-1}$)') plt.ylabel('Wave number (m$^{-1}$)') plt.subplot(223) plt.pcolormesh(f[n//2:], k[:m//2], P[:m//2, n//2:], cmap='magma', vmin=pmin, vmax=0) plt.xlabel('Frequency (s$^{-1}$)') plt.ylabel('Wave number (m$^{-1}$)') plt.subplot(224) plt.pcolormesh(f[n//2:], abs(k[:m//2][::-1]), P[:m//2, n//2:][::-1], cmap='magma', vmin=pmin, vmax=0) cbar = plt.colorbar() cbar.set_label(r'FK spectra (dB)') plt.xlim(0, 100) plt.xlabel('Frequency (s$^{-1}$)') plt.ylabel('Wave number (m$^{-1}$)') plt.tight_layout() plt.show()

这段代码主要实现了面波数据的读取和处理，包括画出原始数据的波形图和进行二维FFT后绘制FK谱图。具体步骤如下：

1. 读取面波数据：使用obspy中的read函数读取指定路径下的所有.SAC文件，并将各个通道的数据保存到data列表中。

2. 画出原始数据的波形图：将每个通道的数据与对应的时间轴绘制在同一张图上，用不同的颜色区分不同通道。

3. 进行二维FFT：将原始数据进行二维FFT，得到FK谱，用于分析面波的频率和波数分布情况。

4. 绘制FK谱图：将FK谱进行可视化，绘制四个子图。其中，第一个子图展示了整个FK谱的情况，第二个子图在第一个子图的基础上加入了一个红色方框，表示主要分析区域。第三个子图展示了主要分析区域的上半部分，第四个子图展示了主要分析区域的下半部分。

总的来说，这段代码实现了对面波数据的初步处理和分析，为后续的面波分析提供了基础。