从地震波到合成记录：手把手教你用Python模拟地震勘探核心流程
当地质工程师面对一片未知的地下区域时，地震勘探就像给地球做CT扫描。本文将带你用Python代码重现这个神奇的过程——从基础波动理论到生成可解释的合成地震记录。不需要昂贵的勘探设备，只需一台电脑和几行代码，你就能在屏幕上看到地下结构的"回声"。
1. 搭建地震波模拟环境
在开始前，确保你的Python环境已安装以下核心库：
import numpy as np import matplotlib.pyplot as plt from scipy.signal import convolve from scipy.fft import fft, ifft关键工具链配置：
- NumPy：处理大型矩阵运算
- SciPy：提供专业科学计算函数
- Matplotlib：可视化波形与记录
提示：推荐使用Jupyter Notebook进行交互式开发，方便实时查看每个步骤的输出结果
地震波速度模型是模拟的基础。我们可以用字典构建一个简化的地层速度剖面：
velocity_model = { '表层': {'vp': 1500, 'vs': 500, 'density': 1.8}, # 单位：m/s, g/cm³ '砂岩层': {'vp': 2800, 'vs': 1500, 'density': 2.2}, '页岩层': {'vp': 3200, 'vs': 1800, 'density': 2.4} }2. 构建地震子波发生器
雷克子波(Ricker Wavelet)是地震模拟中最常用的零相位子波，其数学表达式为：
$$ w(t) = (1-2\pi^2 f^2 t^2)e^{-\pi^2 f^2 t^2} $$
对应的Python实现：
def ricker_wavelet(freq, duration, dt): """ 生成雷克子波 :param freq: 主频(Hz) :param duration: 持续时间(s) :param dt: 时间采样间隔(s) """ t = np.arange(-duration/2, duration/2, dt) return (1 - 2*(np.pi**2)*(freq**2)*(t**2)) * np.exp(-(np.pi**2)*(freq**2)*(t**2))子波参数对比实验：
| 主频(Hz) | 波长(m) | 分辨率 | 适用场景 |
|---|---|---|---|
| 30 | 50 | 较低 | 深层勘探 |
| 60 | 25 | 中等 | 中层勘探 |
| 100 | 12.5 | 较高 | 浅层勘探 |
绘制不同主频子波的效果：
plt.figure(figsize=(10,4)) for freq in [30, 60, 100]: wavelet = ricker_wavelet(freq, 0.1, 0.001) plt.plot(wavelet, label=f'{freq}Hz') plt.legend(); plt.title('不同主频雷克子波对比')3. 计算地层反射系数序列
反射系数是连接地层属性与地震响应的桥梁，计算公式为：
$$ R = \frac{\rho_2 v_2 - \rho_1 v_1}{\rho_2 v_2 + \rho_1 v_1} $$
Python实现步骤：
- 定义地层界面深度
- 计算各层波阻抗(密度×速度)
- 求取界面反射系数
def calculate_reflection_coefficients(model): layers = list(model.keys()) rc = [] for i in range(1, len(layers)): z1 = model[layers[i-1]]['density'] * model[layers[i-1]]['vp'] z2 = model[layers[i]]['density'] * model[layers[i]]['vp'] rc.append((z2 - z1)/(z2 + z1)) return np.array(rc)实际应用示例：
ref_coeff = calculate_reflection_coefficients(velocity_model) print(f"反射系数序列: {ref_coeff}")4. 合成地震记录生成
将子波与反射系数序列进行褶积运算，即得到合成地震记录：
def generate_synthetic_seismogram(wavelet, rc_series, time_window): """ :param wavelet: 地震子波 :param rc_series: 反射系数序列 :param time_window: 时间窗口(s) """ # 创建时间轴 t = np.linspace(0, time_window, len(rc_series)) # 褶积运算 synthetic = convolve(rc_series, wavelet, mode='same') return t, synthetic可视化对比分析：
wavelet = ricker_wavelet(40, 0.2, 0.001) t, syn_seis = generate_synthetic_seismogram(wavelet, ref_coeff, 2) plt.figure(figsize=(12,6)) plt.subplot(211) plt.stem(t[:len(ref_coeff)], ref_coeff, 'r') plt.title('反射系数序列') plt.subplot(212) plt.plot(t, syn_seis) plt.title('合成地震记录') plt.tight_layout()5. 高级模拟技巧与实践
5.1 添加随机噪声模拟真实环境
def add_noise(signal, snr_db): """ 添加高斯白噪声 :param snr_db: 信噪比(dB) """ signal_power = np.mean(signal**2) noise_power = signal_power / (10**(snr_db/10)) noise = np.random.normal(0, np.sqrt(noise_power), len(signal)) return signal + noise5.2 时变子波模拟地层吸收效应
def time_variant_wavelet(freq, t, q_factor): """ 模拟地层吸收导致的子波变化 :param q_factor: 品质因子 """ return ricker_wavelet(freq, 0.1, 0.001) * np.exp(-np.pi*freq*t/q_factor)5.3 多道地震剖面生成
def generate_seismic_profile(model, traces=10, spacing=50): profile = [] for i in range(traces): # 添加随机微调模拟实际勘探 adjusted_model = adjust_velocity_model(model, i*0.01) rc = calculate_reflection_coefficients(adjusted_model) _, seis = generate_synthetic_seismogram(wavelet, rc, 2) profile.append(seis) return np.array(profile).T plt.imshow(generate_seismic_profile(velocity_model), aspect='auto', cmap='seismic')6. 结果解释与地质意义
通过调整参数观察地震记录变化：
- 主频影响：高主频子波能分辨更薄地层，但穿透深度有限
- 噪声水平：信噪比低于15dB时，有效信号开始被淹没
- 速度变化：10%的速度差异会导致同相轴明显偏移
典型陷阱与解决方案：
子波旁瓣干扰：
- 现象：主波峰两侧出现虚假反射
- 对策：应用反褶积处理
调谐效应：
- 现象：薄层产生增强/减弱干涉
- 对策：采用谱分解技术
多次波污染：
- 现象：重复出现的相似波形
- 对策：设计Radon变换滤波器
# 薄层调谐效应演示 thin_layer_model = { '上覆层': {'vp': 3000, 'density': 2.3}, '薄层(5m)': {'vp': 2800, 'density': 2.1}, '下伏层': {'vp': 3200, 'density': 2.4} }在地质解释实践中，我经常发现初学者容易过度依赖自动追踪结果。实际上，合成记录与真实数据的匹配度达到70%就已经是很好的结果，剩余差异往往包含重要的地质信息。
