用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. 反射系数序列计算

让我们首先实现反射系数序列的计算。假设我们有以下模型参数：

• 第一层：厚度h=100m，纵波速度v1=1500m/s，密度ρ1=2000kg/m³
• 第二层：纵波速度v2=2200m/s，密度ρ2=2100kg/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_collection=True) 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, mode='same') # 由于褶积会使信号长度增加，我们需要截取有效部分 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_collection=True) 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[i+1]["v"], layers[i+1]["ρ"] 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, mode='same')[:n_samples]

5.2 性能优化技巧

对于大规模计算，我们可以优化褶积运算：

from scipy.signal import fftconvolve # 使用FFT加速的褶积 s_fast = fftconvolve(wavelet, r, mode='same')[:n_samples]

FFT-based卷积对于长信号效率更高。下表比较了两种方法的性能：

5.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, dt=0.001, t_max=0.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[i+1]["v"], layers[i+1]["ρ"] 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, length=0.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, f=50, noise_level=0): """合成地震记录""" r = self.calculate_reflection_coefficients(layers) wavelet = self.generate_ricker_wavelet(f) s = fftconvolve(wavelet, r, mode='same')[: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, f=50, noise_level=0.05)

这个类封装了所有核心功能，可以方便地用于不同模型参数的地震记录合成。在实际项目中，我通常会先测试简单模型验证代码正确性，再逐步增加复杂度。