用Python从零实现地震波合成：手把手教你用NumPy和Matplotlib搞定褶积模型

用Python从零实现地震波合成手把手教你用NumPy和Matplotlib搞定褶积模型地震勘探是地球物理研究的重要手段而合成地震记录则是理解地震波传播特性的关键工具。本文将带你用Python从头构建一个完整的地震波合成系统通过代码实现反射系数计算、雷克子波生成和褶积运算最终可视化合成地震记录。无论你是地球物理专业的学生还是对科学计算感兴趣的开发者都能从这篇实战教程中获得实用技能。1. 环境准备与基础概念在开始编码前我们需要明确几个核心概念。地震波合成本质上是通过数学模型模拟地震波在地下介质中的传播过程。对于水平层状模型垂直入射的反射系数R可以通过以下公式计算R (ρ2*v2 - ρ1*v1) / (ρ2*v2 ρ1*v1)其中ρ表示密度v表示波速下标1和2分别代表上下两层介质。我们将使用NumPy进行数值计算Matplotlib进行可视化展示。首先安装必要的库pip install numpy matplotlib反射波的双程旅行时t0即波从地表到界面再返回地表的时间计算公式为t0 2*h / v1这里h是第一层的厚度。理解这些基础物理概念后我们就可以开始构建Python实现。2. 反射系数序列计算让我们首先实现反射系数序列的计算。假设我们有以下模型参数第一层厚度h100m纵波速度v11500m/s密度ρ12000kg/m³第二层纵波速度v22200m/s密度ρ22100kg/m³对应的Python实现如下import numpy as np import matplotlib.pyplot as plt # 模型参数 ρ1 2000 # kg/m³ ρ2 2100 v1 1500 # m/s v2 2200 h 100 # m # 计算反射系数和双程旅行时 R (ρ2*v2 - ρ1*v1) / (ρ2*v2 ρ1*v1) t0 2 * h / v1 print(f反射系数R{R:.4f}, 双程旅行时t0{t0:.4f}s)注意在实际地震数据处理中反射系数通常很小绝对值小于0.3这是因为地下介质参数通常是渐变的。接下来我们需要将反射系数表示为时间序列。在地震数据处理中我们通常使用离散时间序列表示信号# 时间序列参数 dt 0.001 # 采样间隔(s) t_max 0.3 # 最大时间(s) n_samples int(t_max / dt) 1 # 采样点数 # 创建反射系数序列 r np.zeros(n_samples) index int(t0 / dt) # 反射系数出现的位置 r[index] R # 可视化 plt.figure(figsize(8, 4)) plt.stem(np.arange(n_samples)*dt, r, basefmt , use_line_collectionTrue) plt.gca().invert_yaxis() # 地震数据惯例时间向下增加 plt.xlabel(时间(s)) plt.ylabel(反射系数) plt.title(反射系数序列r(t)) plt.grid(True) plt.show()3. 雷克子波生成地震子波是地震勘探中的基本波形雷克子波(Ricker Wavelet)是最常用的零相位子波之一。其数学表达式为w(t) (1 - 2(πft)²) * exp(-(πft)²)其中f是主频。让我们实现一个50Hz雷克子波的生成函数def ricker_wavelet(f, dt, length): 生成雷克子波 参数: f: 主频(Hz) dt: 采样间隔(s) length: 子波长度(s) 返回: (时间数组, 振幅数组) t np.arange(-length/2, length/2, dt) wavelet (1 - 2*(np.pi*f*t)**2) * np.exp(-(np.pi*f*t)**2) return t, wavelet # 生成50Hz雷克子波 f 50 # Hz wavelet_length 0.1 # s t_wavelet, wavelet ricker_wavelet(f, dt, wavelet_length) # 可视化 plt.figure(figsize(8, 4)) plt.plot(t_wavelet, wavelet) plt.xlabel(时间(s)) plt.ylabel(振幅) plt.title(50Hz雷克子波) plt.grid(True) plt.show()雷克子波有几个重要特性值得注意主频决定了子波的胖瘦高频子波更窄低频子波更宽零相位特性意味着最大振幅对应反射界面位置子波长度应足够长以确保振幅衰减到接近零4. 褶积运算实现褶积是地震记录合成的核心运算它将反射系数序列与地震子波结合起来。数学上离散褶积定义为s[n] Σ w[k] * r[n-k]在NumPy中我们可以直接使用np.convolve函数实现褶积。但需要注意相位校正问题# 执行褶积运算 s np.convolve(wavelet, r, modesame) # 由于褶积会使信号长度增加我们需要截取有效部分 valid_length min(len(s), n_samples) s s[:valid_length] # 时间轴 time np.arange(valid_length) * dt # 可视化三个信号 plt.figure(figsize(10, 6)) plt.subplot(311) plt.stem(time, r, basefmt , use_line_collectionTrue) plt.gca().invert_yaxis() plt.title(反射系数序列r(t)) plt.grid(True) plt.subplot(312) plt.plot(t_wavelet, wavelet) plt.title(雷克子波w(t)) plt.grid(True) plt.subplot(313) plt.plot(time, s) plt.gca().invert_yaxis() plt.title(合成地震记录s(t)) plt.xlabel(时间(s)) plt.grid(True) plt.tight_layout() plt.show()关键点零相位子波确保了地震记录中波峰/波谷直接对应反射界面位置这是解释地震数据的重要基础。5. 高级应用与优化基础实现完成后我们可以考虑一些高级应用和优化5.1 多层模型实现真实地下通常有多层介质。我们可以扩展代码处理多层情况# 多层模型参数 layers [ {h: 100, v: 1500, ρ: 2000}, # 第一层 {h: 150, v: 1800, ρ: 2100}, # 第二层 {v: 2200, ρ: 2300} # 半空间 ] # 计算各层反射系数和双程时 interfaces [] t_accum 0 for i in range(len(layers)-1): v1, ρ1 layers[i][v], layers[i][ρ] v2, ρ2 layers[i1][v], layers[i1][ρ] R (ρ2*v2 - ρ1*v1) / (ρ2*v2 ρ1*v1) t_accum 2 * layers[i][h] / layers[i][v] interfaces.append({R: R, t: t_accum}) # 创建反射系数序列 r_multi np.zeros(n_samples) for interface in interfaces: index int(interface[t] / dt) if index n_samples: r_multi[index] interface[R] # 合成地震记录 s_multi np.convolve(wavelet, r_multi, modesame)[:n_samples]5.2 性能优化技巧对于大规模计算我们可以优化褶积运算from scipy.signal import fftconvolve # 使用FFT加速的褶积 s_fast fftconvolve(wavelet, r, modesame)[:n_samples]FFT-based卷积对于长信号效率更高。下表比较了两种方法的性能方法信号长度1000信号长度10000信号长度100000直接卷积0.23ms21.5ms2.1sFFT卷积0.45ms1.2ms12.4ms5.3 添加噪声模拟真实地震数据总包含噪声我们可以添加随机噪声模拟实际情况# 添加高斯白噪声 noise_level 0.1 # 噪声水平 noise np.random.normal(0, noise_level * np.max(np.abs(s)), len(s)) s_noisy s noise # 可视化带噪声的地震记录 plt.figure(figsize(8, 4)) plt.plot(time, s_noisy) plt.gca().invert_yaxis() plt.title(含噪声的合成地震记录) plt.xlabel(时间(s)) plt.grid(True) plt.show()6. 完整代码整合将所有功能整合到一个类中便于复用class SeismicSynthesizer: def __init__(self, dt0.001, t_max0.3): self.dt dt self.t_max t_max self.n_samples int(t_max / dt) 1 self.time np.arange(self.n_samples) * dt def calculate_reflection_coefficients(self, layers): 计算多层模型的反射系数序列 r np.zeros(self.n_samples) t_accum 0 for i in range(len(layers)-1): v1, ρ1 layers[i][v], layers[i][ρ] v2, ρ2 layers[i1][v], layers[i1][ρ] R (ρ2*v2 - ρ1*v1) / (ρ2*v2 ρ1*v1) if h in layers[i]: t_accum 2 * layers[i][h] / layers[i][v] index int(t_accum / self.dt) if index self.n_samples: r[index] R return r def generate_ricker_wavelet(self, f, length0.1): 生成雷克子波 t np.arange(-length/2, length/2, self.dt) wavelet (1 - 2*(np.pi*f*t)**2) * np.exp(-(np.pi*f*t)**2) return wavelet def synthesize(self, layers, f50, noise_level0): 合成地震记录 r self.calculate_reflection_coefficients(layers) wavelet self.generate_ricker_wavelet(f) s fftconvolve(wavelet, r, modesame)[:self.n_samples] if noise_level 0: noise np.random.normal(0, noise_level*np.max(np.abs(s)), self.n_samples) s noise return r, s # 使用示例 synth SeismicSynthesizer() layers [ {h: 100, v: 1500, ρ: 2000}, {h: 150, v: 1800, ρ: 2100}, {v: 2200, ρ: 2300} ] r, s synth.synthesize(layers, f50, noise_level0.05)这个类封装了所有核心功能可以方便地用于不同模型参数的地震记录合成。在实际项目中我通常会先测试简单模型验证代码正确性再逐步增加复杂度。