从地震波到地下油藏：手把手教你用Python模拟P波与S波传播（附代码）

从地震波到地下油藏手把手教你用Python模拟P波与S波传播附代码当地质工程师面对复杂的地下构造时地震波就像一把打开地球内部奥秘的钥匙。想象一下当你点击运行按钮屏幕上瞬间呈现出P波在地壳中快速穿行的动态画面而S波则以稍慢的速度紧随其后——这种直观的展示方式比任何教科书上的静态图表都更能帮助我们理解地震波的特性。本文将带你用Python构建一个完整的地震波模拟系统从基础理论到代码实现一步步揭开地震波传播的神秘面纱。1. 地震波基础与Python环境搭建地震波主要分为体波和面波两大类其中体波又包含纵波P波和横波S波。P波的传播速度通常比S波快1.73倍左右这个比值与介质的泊松比直接相关。在实际勘探中P波总是最先到达检波器而S波虽然迟到却携带了更多关于岩石剪切性质的信息。P波与S波的关键区别传播介质P波可在固体和流体中传播S波仅限固体质点振动P波为压缩-膨胀运动S波为剪切运动速度差异Vp/Vs≈√[2(1-ν)/(1-2ν)]ν为泊松比要模拟这些现象我们需要配置以下Python环境# 安装必要库 pip install numpy matplotlib scipy ipywidgets推荐使用Jupyter Notebook进行交互式开发方便实时观察波形变化。核心计算库NumPy将处理矩阵运算Matplotlib负责可视化而SciPy提供必要的科学计算函数。提示对于大规模3D模拟建议考虑使用PyOpenCL或CuPy加速计算普通2D案例用纯NumPy即可满足需求。2. 构建一维波动方程求解器地震波传播遵循弹性波动方程其简化的一维形式可表示为∂²u/∂t² v² ∂²u/∂x²其中u(x,t)为位移场v为波速。我们用有限差分法进行离散化求解import numpy as np import matplotlib.pyplot as plt def simulate_1d_wave(vp2000, dx5, dt0.001, duration1.0): nx 500 # 空间点数 nt int(duration/dt) # 时间步数 # 初始化波场 u np.zeros(nx) u_prev np.copy(u) u_next np.copy(u) # 设置震源Ricker子波 t np.arange(nt)*dt fc 25 # 主频(Hz) src (1-2*(np.pi*fc*t)**2)*np.exp(-(np.pi*fc*t)**2) # 有限差分系数 c (vp*dt/dx)**2 # 时间迭代 for it in range(nt): # 注入震源 u[10] src[it] # 更新波场 u_next[1:-1] 2*u[1:-1] - u_prev[1:-1] c*(u[2:] - 2*u[1:-1] u[:-2]) # 边界条件简单吸收 u_next[0] u[1] u_next[-1] u[-2] # 更新状态 u_prev[:] u u[:] u_next # 实时显示 if it % 10 0: plt.clf() plt.plot(u) plt.ylim(-1, 1) plt.title(fTime {it*dt:.3f}s) plt.pause(0.01) return u这段代码实现了Ricker子波震源生成二阶空间导数的中心差分近似时间上的显式迭代求解简单边界处理参数敏感性分析参数物理意义典型值范围影响效果vpP波速度1500-6000 m/s速度越高波长越长dx空间步长1-10 m越小精度越高但计算量增大dt时间步长0.1-1 ms需满足CFL稳定性条件fc主频10-50 Hz高频成分分辨率更高3. 二维弹性波场模拟进阶扩展到二维情况后我们需要同时考虑P波和S波的耦合作用。各向同性介质中的弹性波方程可表示为ρ ∂²u/∂t² (λ2μ)∇(∇·u) - μ∇×(∇×u)其中u(ux,uz)为位移矢量场λ和μ为拉梅常数。对应的Python实现def elastic_2d_simulation(): # 模型参数 nx, nz 200, 100 dx, dz 10, 10 # 网格间距(m) dt 0.001 # 时间步长(s) nt 1000 # 时间步数 # 介质参数两层模型 vp np.ones((nz, nx))*2000 vs np.ones((nz, nx))*1200 rho np.ones((nz, nx))*2000 vp[nz//2:,:] 2500 # 下层速度更高 vs[nz//2:,:] 1500 rho[nz//2:,:] 2200 # 计算模量 mu rho*vs**2 lam rho*vp**2 - 2*mu # 初始化波场 ux np.zeros((nz, nx)) uz np.zeros((nz, nx)) vx np.zeros((nz, nx)) vz np.zeros((nz, nx)) sxx np.zeros((nz, nx)) szz np.zeros((nz, nx)) sxz np.zeros((nz, nx)) # 震源位置和参数 src_x, src_z nx//2, 10 fc 25 # 主频(Hz) for it in range(nt): t it*dt # 注入震源垂直力源 src_val np.exp(-(np.pi*fc*(t-0.1))**2) vz[src_z, src_x] src_val * dt/rho[src_z, src_x] # 更新速度场 vx[:,:] dt/rho * (np.diff(sxx, axis1, append0)/dx np.diff(sxz, axis0, append0)/dz) vz[:,:] dt/rho * (np.diff(sxz, axis1, append0)/dx np.diff(szz, axis0, append0)/dz) # 更新应力场 dux_dx np.diff(ux, axis1, prepend0)/dx duz_dz np.diff(uz, axis0, prepend0)/dz dux_dz np.diff(ux, axis0, prepend0)/dz duz_dx np.diff(uz, axis1, prepend0)/dx sxx[:,:] dt * ((lam2*mu)*dux_dx lam*duz_dz) szz[:,:] dt * (lam*dux_dx (lam2*mu)*duz_dz) sxz[:,:] dt * mu * (dux_dz duz_dx) # 更新位移场 ux vx*dt uz vz*dt # 可视化 if it % 10 0: plt.figure(figsize(12,5)) plt.subplot(121) plt.imshow(ux, cmapseismic, vmin-1e-6, vmax1e-6) plt.title(fHorizontal displacement at t{t:.3f}s) plt.colorbar() plt.subplot(122) plt.imshow(uz, cmapseismic, vmin-1e-6, vmax1e-6) plt.title(Vertical displacement) plt.colorbar() plt.tight_layout() plt.pause(0.01) plt.close()这段代码展示了速度-应力格式的弹性波方程求解交错网格有限差分方法点震源激发机制双层介质模型设置典型地质界面模拟结果对比界面类型P波反射特征S波反射特征转换波产生硬-软界面反射系数正反射系数负强PS转换软-硬界面反射系数负反射系数正强SP转换流体-固体全反射无透射S波显著波型转换4. 合成地震记录制作实战有了波场模拟基础我们可以进一步制作合成地震记录——这是连接理论模拟与实际勘探数据的桥梁。基本流程包括构建地下模型def create_layer_model(): # 三层模型页岩-砂岩-灰岩 vp np.array([2500, 3000, 4500]) # P波速度(m/s) vs np.array([1200, 1500, 2500]) # S波速度 rho np.array([2300, 2400, 2600]) # 密度(kg/m³) thickness [500, 800] # 层厚度(m) return vp, vs, rho, thickness计算反射系数def calculate_reflection_coefficients(vp, vs, rho): # 垂直入射P波反射系数 z vp * rho # 波阻抗 rc (z[1:]-z[:-1])/(z[1:]z[:-1]) return rc子波生成与褶积def generate_synthetic_seismogram(): vp, vs, rho, thickness create_layer_model() rc calculate_reflection_coefficients(vp, vs, rho) # 生成Ricker子波 dt 0.002 # 采样间隔(s) t np.arange(-0.1, 0.1, dt) fc 30 # 主频(Hz) wavelet (1-2*(np.pi*fc*t)**2)*np.exp(-(np.pi*fc*t)**2) # 构建反射系数序列 sample_per_layer [int(t/dt) for t in thickness] rc_series np.zeros(sum(sample_per_layer)) rc_series[0] rc[0] rc_series[sample_per_layer[0]] rc[1] # 褶积生成合成记录 synth np.convolve(rc_series, wavelet, modesame) # 添加噪声 noise 0.05*np.random.randn(len(synth)) return synth noise合成记录质量影响因素子波参数选择主频过高 → 分辨率提升但信噪比下降主频过低 → 层间干涉严重相位误差 → 同相轴位置偏移模型参数准确性速度误差5% → 时间偏移可达数十ms密度估计偏差 → 振幅信息失真薄层厚度误差 → 调谐效应变化在实际项目中我们通常会使用更专业的库如bruges来处理复杂的地震正演问题import bruges as bg # 使用bruges生成更专业的合成记录 vp [1500, 2000, 2500] rho [1800, 2100, 2300] thickness [100, 150] rc bg.reflection.acoustic_reflectivity(vp, rho) wavelet bg.filters.ricker(duration0.1, dt0.001, f30) synth bg.models.convolve(rc, wavelet)5. 实际应用案例油气储层响应分析通过前面的技术积累我们现在可以模拟一个完整的油气藏地震响应。假设我们有一个砂泥岩互层储层其中砂岩含气def gas_sand_simulation(): # 五层模型上覆泥岩-盖层-含气砂岩-含水砂岩-基底 vp [2200, 2400, 1800, 2100, 2800] # m/s vs [1100, 1200, 900, 1100, 1500] # m/s rho [2100, 2250, 1900, 2200, 2500] # kg/m³ thickness [300, 100, 50, 150] # m # 计算反射系数 rc_p bg.reflection.reflectivity(vp, vs, rho, theta0) rc_s bg.reflection.reflectivity(vp, vs, rho, theta0, waveS) # 生成子波 wavelet_p bg.filters.ricker(duration0.1, dt0.001, f35) wavelet_s bg.filters.ricker(duration0.15, dt0.001, f25) # 制作合成记录 synth_p bg.models.convolve(rc_p, wavelet_p) synth_s bg.models.convolve(rc_s, wavelet_s) # 含气砂岩特征分析 gas_sand_amplitude synth_p[thickness[0]thickness[1]] brine_sand_amplitude synth_p[sum(thickness[:3])] amplitude_ratio gas_sand_amplitude / brine_sand_amplitude plt.figure(figsize(10,6)) plt.plot(synth_p, labelP波记录) plt.plot(synth_s, labelS波记录) plt.axvline(thickness[0]thickness[1], colorr, linestyle--) plt.axvline(sum(thickness[:3]), colorb, linestyle--) plt.legend() plt.title(f含气砂岩振幅比: {amplitude_ratio:.2f}) plt.show()油气检测关键指标振幅异常含气砂岩通常表现为强负反射亮点振幅随偏移距变化AVO特征速度比异常Vp/Vs比值明显降低通常1.5泊松比显著下降频率衰减高频成分被强烈吸收主频向低频方向移动储层参数敏感性表参数变化P波速度影响S波速度影响密度影响综合地震响应孔隙度↑显著↓显著↓轻微↓振幅↑频率↓含气饱和度↑显著↓轻微↓轻微↓强振幅异常泥质含量↑↓↓↓↑振幅↓频率↑压力↑↑↑-振幅↓速度↑在完成这些模拟后建议将结果与实际野外数据进行对比分析。一个常见的实践是使用公开数据集如SEG盐体模型或Marmousi模型进行基准测试评估模拟算法的可靠性。