从地震波到AI预测：用Python模拟P波与S波传播（附完整代码）

用Python模拟地震波P波与S波的动态可视化与物理机制解析地震波是地球物理勘探和地震学研究的基础理解P波纵波和S波横波的传播特性对于预测地震影响、分析地下结构至关重要。本文将带你用Python构建一个完整的波传播模拟系统从基础物理公式推导到交互式动画生成最后探讨这些模拟数据如何训练机器学习模型。不同于教科书式的理论讲解我们会通过代码让抽象的地震波物理变得直观可见。1. 地震波基础与数值模拟原理地震波本质上是在弹性介质中传播的机械振动。P波和S波作为体波Body Wave的两种基本类型其区别不仅体现在传播方向上更反映在介质响应机制中。P波通过介质的压缩和膨胀传递能量而S波则依赖剪切形变。这种本质差异导致了两者在传播速度、穿透能力等方面的显著区别。波动方程的离散化是数值模拟的核心。对于一维P波传播我们可以从简化的声波方程出发# 一维声波方程离散形式 def p_wave_1d_simulation(c, dx, dt, nt): c: 波速 (m/s) dx: 空间步长 (m) dt: 时间步长 (s) nt: 时间步数 nx int(1.0 / dx) 1 u np.zeros(nx) # 当前时刻波场 u_prev np.zeros(nx) # 前一时刻波场 u_next np.zeros(nx) # 下一时刻波场 # 初始化脉冲源 u[nx//2] 1.0 u_prev[nx//2] 1.0 # 有限差分迭代 for it in range(1, nt): u_next[1:-1] 2*u[1:-1] - u_prev[1:-1] \ (c*dt/dx)**2 * (u[2:] - 2*u[1:-1] u[:-2]) u_prev, u u, u_next yield u对于S波模拟需要考虑横向位移和剪切模量的影响。各向同性介质中的S波速度计算公式为$$ v_s \sqrt{\frac{\mu}{\rho}} $$其中$\mu$为剪切模量$\rho$为介质密度。相比之下P波速度公式为$$ v_p \sqrt{\frac{K \frac{4}{3}\mu}{\rho}} $$$K$表示体积模量。这两个公式直接解释了为什么在相同介质中$v_p v_s$以及为什么S波无法在流体中传播流体$\mu0$。表常见岩石的P波和S波速度对比岩石类型密度 (g/cm³)VP (km/s)VS (km/s)VP/VS 比值花岗岩2.6-2.75.0-6.03.0-3.51.65-1.75玄武岩2.8-3.06.0-7.03.5-4.01.70-1.80砂岩2.2-2.63.0-5.01.8-3.01.60-1.70页岩2.0-2.42.5-4.01.5-2.51.65-1.802. 二维波场模拟与Matplotlib动态可视化将一维模拟扩展到二维空间可以更真实地展现波前传播和不同介质边界的反射/折射现象。我们使用交错网格有限差分法来实现这一扩展def wave_propagation_2d(vp_model, vs_model, rho_model, dx, dz, dt, nt): vp_model: 2D P波速度模型 (m/s) vs_model: 2D S波速度模型 (m/s) rho_model: 2D 密度模型 (kg/m³) dx, dz: 空间采样间隔 (m) dt: 时间采样间隔 (s) nt: 时间步数 nx, nz vp_model.shape mu rho_model * vs_model**2 # 剪切模量 lam rho_model * vp_model**2 - 2*mu # 拉梅常数 # 初始化波场变量 vx np.zeros((nx, nz)) # x方向质点速度 vz np.zeros((nx, nz)) # z方向质点速度 sxx np.zeros((nx, nz)) # 正应力xx分量 szz np.zeros((nx, nz)) # 正应力zz分量 sxz np.zeros((nx, nz)) # 剪应力xz分量 # 震源函数 - Ricker小波 def ricker(freq, t): arg (np.pi * freq * (t - 1.0/freq))**2 return (1.0 - 2.0*arg) * np.exp(-arg) # 时间迭代 for it in range(nt): # 更新应力分量 sxx[1:,1:] dt * (lam[1:,1:]2*mu[1:,1:]) * (vx[1:,1:]-vx[:-1,1:])/dx \ dt * lam[1:,1:] * (vz[1:,1:]-vz[1:,:-1])/dz szz[1:,1:] dt * lam[1:,1:] * (vx[1:,1:]-vx[:-1,1:])/dx \ dt * (lam[1:,1:]2*mu[1:,1:]) * (vz[1:,1:]-vz[1:,:-1])/dz sxz[:-1,:-1] dt * mu[:-1,:-1] * ((vx[:-1,1:]-vx[:-1,:-1])/dz \ (vz[1:,:-1]-vz[:-1,:-1])/dx) # 更新速度分量 t it * dt source 1e8 * ricker(25, t) # 25Hz主频震源 vx[10,10] source # 点源位置 vx[:-1,1:-1] dt/rho_model[:-1,1:-1] * \ ((sxx[1:,1:-1]-sxx[:-1,1:-1])/dx \ (sxz[:-1,1:-1]-sxz[:-1,:-2])/dz) vz[1:-1,:-1] dt/rho_model[1:-1,:-1] * \ ((sxz[1:-1,:-1]-sxz[:-2,:-1])/dx \ (szz[1:-1,1:]-szz[1:-1,:-1])/dz) yield vx, vz, sxxszz # 返回速度场和标量波场注意实际模拟时需要满足CFL稳定性条件dt ≤ min(dx,dz)/√2/max(vp)。对于典型的地壳速度5km/s和10m网格间距时间步长应小于1.4ms。利用Matplotlib的FuncAnimation可以创建生动的传播过程动画。以下代码片段展示如何生成交互式可视化def create_wave_animation(simulator, frames100): fig, (ax1, ax2) plt.subplots(1, 2, figsize(12,5)) # 初始化图像 img1 ax1.imshow(np.zeros((100,100)), cmapseismic, vmin-1e-5, vmax1e-5) ax1.set_title(P波波场) img2 ax2.imshow(np.zeros((100,100)), cmapseismic, vmin-1e-6, vmax1e-6) ax2.set_title(S波波场) def update(frame): vx, vz, p next(simulator) img1.set_array(p.T) img2.set_array(vx.T) return img1, img2 ani animation.FuncAnimation(fig, update, framesframes, interval50, blitTrue) plt.close() return HTML(ani.to_jshtml())图1展示的模拟结果中可以清晰观察到P波波前呈圆形对称扩散传播速度明显快于S波在流体层如设置vp1500m/svs0m/s中S波能量完全消失遇到高速层时波前发生明显折射现象自由表面产生复杂的反射波和转换波3. 复杂地质模型构建与波场特征分析真实地下结构远非均匀介质我们的模拟需要包含断层、不整合面、盐丘等复杂构造。通过定义速度模型函数可以灵活构建各种地质场景def create_geologic_model(nx300, nz150, dx10): 构建包含断层、沉积层和盐丘的速度模型 vp np.ones((nx,nz)) * 3000 # 背景速度3km/s vs np.ones((nx,nz)) * 1800 rho np.ones((nx,nz)) * 2.4 # 沉积层序列 layers [ {z0:30, dz:20, vp:2500, vs:1500, rho:2.2}, {z0:50, dz:30, vp:3200, vs:1900, rho:2.5}, {z0:80, dz:40, vp:3800, vs:2300, rho:2.7} ] for layer in layers: z_start int(layer[z0]/dx) thickness int(layer[dz]/dx) vp[:,z_start:z_startthickness] layer[vp] vs[:,z_start:z_startthickness] layer[vs] rho[:,z_start:z_startthickness] layer[rho] # 添加倾斜断层 for i in range(nx): offset int(0.2*i) if 100offset nz: vp[i,100offset:] 500 vs[i,100offset:] 300 rho[i,100offset:] 0.3 # 添加盐丘构造 salt_x, salt_z 180, 60 salt_r 40 for i in range(nx): for j in range(nz): if (i-salt_x)**2 (j-salt_z)**2 salt_r**2: vp[i,j] 4500 # 盐岩典型速度 vs[i,j] 2600 rho[i,j] 2.2 return vp, vs, rho典型构造的波场响应特征断层效应断层面产生清晰的反射波断层下盘出现波场阴影区断面两侧波至时间出现跳跃盐丘响应高速盐体导致波前畸变盐体边界产生强反射盐下区域出现速度下拉假象薄层调谐效应当层厚度接近1/4波长时产生相长干涉调谐厚度计算公式$h_{tuning} \frac{v}{4f}$对于30Hz主频和3000m/s速度调谐厚度约25m表不同地质构造的典型地震响应特征构造类型P波响应特征S波响应特征识别要点断层同相轴错断、断面反射更明显的波场中断断面两侧波组关系背斜同相轴向上弯曲类似P波但分辨率低顶部是否有振幅异常河道下切谷反射、内部杂乱更好的砂泥岩区分平面形态分析盐丘顶部强反射、内部透明能量快速衰减速度异常分析4. 从数值模拟到机器学习地震波属性预测地震波模拟生成的海量数据为机器学习提供了理想的训练样本。我们可以构建一个端到端的流程将波动方程模拟与深度学习相结合数据生成流程随机生成包含不同构造的速度模型运行波动方程模拟获取波场快照提取振幅、相位、频率等属性标注对应的岩性、流体等信息def generate_training_data(num_samples1000): 生成训练数据集 features [] labels [] for _ in range(num_samples): # 随机生成速度模型 vp, vs, rho random_velocity_model() # 运行模拟获取波场 simulator wave_propagation_2d(vp, vs, rho, dx10, dz10, dt0.001, nt500) wavefield np.array([p for _,_,p in simulator]) # 提取特征 amp np.max(np.abs(wavefield), axis0) freq dominant_frequency(wavefield) phase instantaneous_phase(wavefield) # 生成标签岩性分类 lithology classify_lithology(vp, vs, rho) features.append(np.stack([amp, freq, phase], axis-1)) labels.append(lithology) return np.array(features), np.array(labels)深度学习模型架构示例PyTorchclass WaveNet(nn.Module): def __init__(self, in_channels3, num_classes5): super().__init__() self.encoder nn.Sequential( nn.Conv2d(in_channels, 32, 3, padding1), nn.ReLU(), nn.MaxPool2d(2), nn.Conv2d(32, 64, 3, padding1), nn.ReLU(), nn.MaxPool2d(2), nn.Conv2d(64, 128, 3, padding1), nn.ReLU(), nn.MaxPool2d(2) ) self.decoder nn.Sequential( nn.ConvTranspose2d(128, 64, 3, stride2), nn.ReLU(), nn.ConvTranspose2d(64, 32, 3, stride2), nn.ReLU(), nn.ConvTranspose2d(32, num_classes, 3, stride2, output_padding1) ) def forward(self, x): x self.encoder(x) x self.decoder(x) return x模型训练关键技巧使用弹性变形(Elastic Deformation)增强训练数据采用Dice损失函数处理类别不平衡引入注意力机制聚焦关键波场特征使用迁移学习从合成数据过渡到实际数据在实际项目中这种融合物理模拟和机器学习的方法已经展现出巨大潜力。例如某研究团队使用20000个合成模型训练的神经网络将地震波速预测准确率提高了40%特别是在复杂构造区域的预测效果显著优于传统反演方法。