别再死记硬背了用PythonNumPy手把手模拟地震子波合成与分辨率分析地震勘探作为地球物理学的核心领域之一其理论体系往往让初学者望而生畏。传统的学习方式强调公式推导和概念记忆但今天我们将打破常规——通过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 python3.8 numpy scipy matplotlib地震信号处理离不开几个关键参数我们先定义基础常量# 时间轴参数 sample_rate 1000 # 采样率(Hz) duration 1.0 # 信号时长(s) t np.linspace(0, duration, int(sample_rate * duration), endpointFalse) # 地层参数示例 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_shift30): 最小相位子波 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, d1/sample_rate)[:n//2] spectrum np.abs(np.fft.fft(wavelet))[:n//2] plt.plot(freq, 20*np.log10(spectrum), color, labellabel) 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[i1]-impedance[i])/(impedance[i1]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, modesame) 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, markerfmtro, linefmtr--, basefmt , label反射系数) plt.title(合成地震记录与反射系数对比) plt.xlabel(双程走时(s)); plt.ylabel(振幅) plt.legend(); plt.grid(True)4. 分辨率极限的数值实验4.1 纵向分辨率薄层识别能力通过改变地层厚度观察波形变化def test_vertical_resolution(wavelet, thickness_range, velocity2000): 纵向分辨率测试 plt.figure(figsize(12,8)) for i, thickness in enumerate(thickness_range): # 创建薄层模型 depths np.array([0, 100, 100thickness, 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,i1) plt.plot(t, synthetic, labelf厚度{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, levels20, cmapviridis) 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_shift60) # 可视化对比 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滤波补偿高频衰减再应用谱蓝化技术拓宽频带最后通过反褶积压缩子波。但要注意避免过度处理引入虚假信息。