别再死记硬背了！用Python+NumPy手把手模拟地震子波合成与分辨率分析
地震勘探作为地球物理学的核心领域之一，其理论体系往往让初学者望而生畏。传统的学习方式强调公式推导和概念记忆，但今天我们将打破常规——通过Python代码实现地震子波合成与分辨率分析的完整流程，让抽象理论变得触手可及。无论你是地质学专业的学生，还是刚进入能源行业的工程师，这套基于NumPy和Matplotlib的实践方案，将帮助你建立直观的物理认知。
1. 环境配置与基础工具链搭建
工欲善其事，必先利其器。我们选择Jupyter Notebook作为交互环境，配合Python科学计算三件套：
# 基础工具包安装 import numpy as np import matplotlib.pyplot as plt from scipy.signal import ricker, convolve plt.style.use('seaborn-whitegrid') # 设置美观的绘图风格提示：建议使用Python 3.8+环境，可通过conda创建专属虚拟环境：
conda create -n seismology python=3.8 numpy scipy matplotlib
地震信号处理离不开几个关键参数，我们先定义基础常量：
# 时间轴参数 sample_rate = 1000 # 采样率(Hz) duration = 1.0 # 信号时长(s) t = np.linspace(0, duration, int(sample_rate * duration), endpoint=False) # 地层参数示例 layer_thickness = [30, 15, 8] # 地层厚度(m) wave_speed = [2000, 2500, 1800] # 纵波速度(m/s)2. 地震子波生成与特征分析
2.1 三种典型子波的代码实现
地震子波是地震记录的"指纹"，不同震源产生的子波形态各异。我们实现三种常见子波：
def ricker_wavelet(freq, t): """雷克子波生成器""" return (1 - 2*(np.pi*freq*t)**2) * np.exp(-(np.pi*freq*t)**2) def ormsby_wavelet(freq_band, t): """奥姆斯比带限子波""" f1, f2, f3, f4 = freq_band numerator = (np.sinc(f4*t)**2)*f4**2 - (np.sinc(f3*t)**2)*f3**2 denominator = f4 - f3 term1 = numerator / denominator numerator = (np.sinc(f2*t)**2)*f2**2 - (np.sinc(f1*t)**2)*f1**2 denominator = f2 - f1 term2 = numerator / denominator return term1 - term2 def minimum_phase_wavelet(freq, t, phase_shift=30): """最小相位子波""" wavelet = ricker_wavelet(freq, t) n = len(wavelet) spectrum = np.fft.fft(wavelet) frequencies = np.fft.fftfreq(n) spectrum *= np.exp(1j * np.deg2rad(phase_shift * frequencies)) return np.real(np.fft.ifft(spectrum))可视化对比（代码示例）：
# 参数设置 central_freq = 30 # 主频(Hz) freq_band = [10, 20, 40, 50] # 频带参数 # 生成子波 ricker = ricker_wavelet(central_freq, t - 0.2) ormsby = ormsby_wavelet(freq_band, t - 0.3) min_phase = minimum_phase_wavelet(central_freq, t - 0.25) # 绘制对比 plt.figure(figsize=(10,6)) plt.plot(t, ricker/np.max(ricker), label='雷克子波') plt.plot(t, ormsby/np.max(ormsby), label='奥姆斯比子波') plt.plot(t, min_phase/np.max(min_phase), label='最小相位子波') plt.title('三种典型地震子波形态对比') plt.xlabel('时间(s)'); plt.ylabel('归一化振幅') plt.legend(); plt.grid(True)2.2 子波频谱特征解析
时域波形只是故事的一半，频域分析同样重要：
def plot_spectrum(wavelet, sample_rate, color, label): n = len(wavelet) freq = np.fft.fftfreq(n, d=1/sample_rate)[:n//2] spectrum = np.abs(np.fft.fft(wavelet))[:n//2] plt.plot(freq, 20*np.log10(spectrum), color, label=label) plt.figure(figsize=(12,5)) plot_spectrum(ricker, sample_rate, 'b', '雷克子波') plot_spectrum(ormsby, sample_rate, 'r', '奥姆斯比子波') plot_spectrum(min_phase, sample_rate, 'g', '最小相位子波') plt.title('子波振幅谱对比(dB)') plt.xlabel('频率(Hz)'); plt.ylabel('振幅(dB)') plt.legend(); plt.grid(True) plt.xlim(0, 100)通过频谱分析可以清晰看到：
- 雷克子波的频带最宽，但存在旁瓣泄漏
- 奥姆斯比子波具有明确的带限特征
- 最小相位子波能量更集中在前端
3. 合成地震记录生成实战
3.1 反射系数模型构建
地层界面反射系数是合成记录的基础：
def build_reflection_model(depths, velocities, densities): """构建反射系数序列""" impedance = velocities * densities rc = np.zeros(len(impedance)-1) for i in range(len(rc)): rc[i] = (impedance[i+1]-impedance[i])/(impedance[i+1]+impedance[i]) return rc # 示例地层参数 depths = np.array([0, 100, 150, 300]) # 地层界面深度(m) velocities = np.array([1800, 2100, 2400, 2000]) # 纵波速度(m/s) densities = np.array([2.1, 2.3, 2.4, 2.2]) # 密度(g/cm³) ref_coeff = build_reflection_model(depths, velocities, densities) print(f"反射系数序列: {ref_coeff}")3.2 褶积运算与合成记录
将子波与反射系数序列进行褶积：
def generate_synthetic(rc, wavelet, sample_rate, depth, velocity): """生成合成地震记录""" # 计算双程走时 two_way_time = 2 * depth / velocity # 创建时间序列 synthetic = np.zeros_like(t) # 执行褶积 for i in range(len(rc)): index = int(two_way_time[i] * sample_rate) if index < len(synthetic): synthetic[index] = rc[i] synthetic = convolve(synthetic, wavelet, mode='same') return synthetic # 生成合成记录 synth_seismic = generate_synthetic(ref_coeff, ricker, sample_rate, depths[1:], velocities[1:]) # 可视化结果 plt.figure(figsize=(12,6)) plt.plot(t, synth_seismic/np.max(synth_seismic), 'b', label='合成记录') plt.stem(depths[1:]/velocities[1:]*2, ref_coeff, 'r', markerfmt='ro', linefmt='r--', basefmt=' ', label='反射系数') plt.title('合成地震记录与反射系数对比') plt.xlabel('双程走时(s)'); plt.ylabel('振幅') plt.legend(); plt.grid(True)4. 分辨率极限的数值实验
4.1 纵向分辨率：薄层识别能力
通过改变地层厚度观察波形变化：
def test_vertical_resolution(wavelet, thickness_range, velocity=2000): """纵向分辨率测试""" plt.figure(figsize=(12,8)) for i, thickness in enumerate(thickness_range): # 创建薄层模型 depths = np.array([0, 100, 100+thickness, 300]) rc = build_reflection_model(depths, velocities, densities) synthetic = generate_synthetic(rc, wavelet, sample_rate, depths[1:], velocities[1:]) plt.subplot(len(thickness_range),1,i+1) plt.plot(t, synthetic, label=f'厚度={thickness}m') plt.legend(); plt.grid(True) plt.suptitle('不同地层厚度下的合成记录响应') # 测试不同厚度(从30m递减到5m) thickness_test = [30, 15, 8, 5] test_vertical_resolution(ricker, thickness_test)实验结果显示：
- 厚度>λ/4时，顶底反射清晰可分
- 厚度≈λ/8时，反射开始叠加
- 厚度<λ/10时，完全无法分辨
4.2 横向分辨率：菲涅尔带模拟
def fresnel_zone(freq, depth, velocity): """计算第一菲涅尔带半径""" wavelength = velocity / freq return np.sqrt(0.5 * wavelength * depth) # 参数敏感性分析 frequencies = np.linspace(10, 100, 10) depths = np.linspace(100, 3000, 10) F, D = np.meshgrid(frequencies, depths) R = fresnel_zone(F, D, 2500) # 可视化 plt.figure(figsize=(10,6)) contour = plt.contourf(F, D, R, levels=20, cmap='viridis') plt.colorbar(contour, label='菲涅尔带半径(m)') plt.title('横向分辨率与频率、深度的关系') plt.xlabel('频率(Hz)'); plt.ylabel('深度(m)')从模拟结果可见：
- 高频信号可显著减小菲涅尔带
- 深层目标的分辨率急剧下降
- 速度变化对分辨率有非线性影响
5. 实际应用中的陷阱与解决方案
5.1 子波提取常见问题
def wavelet_extraction_issues(): """演示子波提取中的典型问题""" # 案例1：噪声干扰 clean_wavelet = ricker_wavelet(30, t-0.2) noisy_wavelet = clean_wavelet + 0.2*np.random.randn(len(t)) # 案例2：相位误差 phase_error_wavelet = minimum_phase_wavelet(30, t-0.2, phase_shift=60) # 可视化对比 plt.figure(figsize=(12,6)) plt.plot(t, clean_wavelet, 'b', label='理想子波') plt.plot(t, noisy_wavelet, 'r--', label='含噪声子波') plt.plot(t, phase_error_wavelet, 'g:', label='相位错误子波') plt.title('子波提取中的典型问题'); plt.grid(True) plt.legend() wavelet_extraction_issues()应对策略：
- 使用多道统计方法提取平均子波
- 结合井资料进行标定校正
- 应用盲反褶积技术优化估计
5.2 分辨率增强技巧
def resolution_enhancement(original): """分辨率增强处理演示""" # 频谱平衡 spectrum = np.fft.fft(original) freq = np.fft.fftfreq(len(original)) enhanced_spec = spectrum * (1 + 0.5*freq**2) balanced = np.real(np.fft.ifft(enhanced_spec)) # 反褶积处理 wavelet = ricker_wavelet(30, t-0.2) deconvolved = deconvolve(original, wavelet) return balanced, deconvolved # 应用示例 balanced, deconvolved = resolution_enhancement(synth_seismic) plt.figure(figsize=(12,6)) plt.plot(t, synth_seismic/np.max(synth_seismic), 'b', label='原始信号') plt.plot(t, balanced/np.max(balanced), 'r--', label='频谱平衡') plt.plot(t, deconvolved/np.max(deconvolved), 'g:', label='反褶积结果') plt.legend(); plt.grid(True)实际项目中，我们发现结合多种处理技术能获得最佳效果。例如先进行反Q滤波补偿高频衰减，再应用谱蓝化技术拓宽频带，最后通过反褶积压缩子波。但要注意避免过度处理引入虚假信息。
