用Python动态模拟地震波：从代码实现理解P波与S波本质

当第一次接触地震波理论时，那些关于P波、S波的抽象描述总让人困惑——为什么纵波能在液体中传播而横波不能？为什么面波破坏力更强？传统教材的静态图示难以展现波动的动态过程，而这正是编程可视化能够突破的认知边界。本文将用不到100行Python代码，带您构建一个可交互的地震波模拟器，在动态图形中直观理解波的本质差异。

1. 环境配置与基础波动模型

1.1 搭建Python科学计算环境

推荐使用Anaconda创建专属虚拟环境，避免依赖冲突：

conda create -n seismic python=3.9 conda activate seismic pip install numpy matplotlib ipywidgets

核心工具链选择：

• NumPy：处理大规模矩阵运算
• Matplotlib：生成高质量动态可视化
• IPywidgets：创建交互式控制面板

1.2 一维波动方程离散化

采用有限差分法模拟波在介质中的传播，核心方程如下：

def wave_propagation(u_prev, u_current, c, dt, dx): """一维波动方程有限差分解算""" u_next = 2*u_current - u_prev + (c*dt/dx)**2 * ( np.roll(u_current,1) - 2*u_current + np.roll(u_current,-1)) return u_next

参数说明：

注意：Courant条件要求c*dt/dx ≤ 1以保证数值稳定性

2. P波与S波的动态模拟对比

2.1 纵波（P波）特性实现

P波模拟关键点在于质点振动与传播方向平行。以下代码创建压缩-膨胀交替的波动效果：

def simulate_p_wave(steps=500): # 初始化位移场 u = np.zeros(NX) u_prev = np.hanning(NX) * 0.5 # 初始扰动 for _ in range(steps): u_next = wave_propagation(u_prev, u, c_p, DT, DX) # 边界吸收条件 u_next[0] = u_next[-1] = 0 yield u_next u_prev, u = u, u_next

可视化技巧：用箭头图标注质点运动方向：

quiv = ax.quiver(x, np.zeros_like(x), u, np.zeros_like(u), scale=20, color='r')

2.2 横波（S波）的特殊约束

S波的介质刚性要求体现在代码中需要：

1. 设置剪切模量参数
2. 添加固体介质检测逻辑：

if medium == 'fluid': raise ValueError("横波不能在流体中传播！")

典型参数对比：

实验：尝试修改介质密度ρ，观察波速变化规律

3. 复杂波场现象模拟

3.1 莫霍面波速跃变模拟

通过设置速度突变边界模拟地壳-地幔界面：

# 创建速度模型 c = np.where(x < NX//2, 6.8, 8.1) # 莫霍面位置

关键现象观察：

• 波到达界面时部分能量反射
• 透射波发生折射（Snell定律）
• 使用能量衰减公式模拟实际传播：

u_next *= np.exp(-0.01*x) # 指数衰减

3.2 面波生成与频散效应

瑞利波模拟需要：

1. 添加自由表面边界条件
2. 实现椭圆极化运动：

def rayleigh_wave(): u_vertical = simulate_p_wave() * 0.7 u_horizontal = simulate_s_wave() * 0.3 return u_vertical + u_horizontal * 1j # 复合振动

频散现象实现技巧：

freqs = np.linspace(0.1, 0.5, 5) # 多频率成分 waves = [np.sin(2*np.pi*f*x) for f in freqs] composite = sum(w * np.exp(-0.1*i) for i,w in enumerate(waves))

4. 地震勘探应用实战

4.1 合成地震记录生成

完整工作流代码框架：

def synthetic_seismogram(): # 1. 加载测井数据 velocity, density = load_well_log() # 2. 计算反射系数 rc = (velocity[1:]*density[1:] - velocity[:-1]*density[:-1]) / \ (velocity[1:]*density[1:] + velocity[:-1]*density[:-1]) # 3. 雷克子波生成 wavelet = ricker(freq=25, length=0.1) # 4. 褶积运算 return np.convolve(rc, wavelet, mode='same')

4.2 波场快照技术

实现多时刻波场状态保存与回放：

from matplotlib.animation import FuncAnimation def update(frame): line.set_ydata(snapshots[frame]) return line, ani = FuncAnimation(fig, update, frames=len(snapshots), interval=50, blit=True)

交互功能增强：

• 添加滑块控制波速参数
• 实现介质类型实时切换
• 支持波前标记与传播计时

from ipywidgets import interact @interact(v_p=(2.0, 8.0, 0.1), v_s=(1.0, 4.5, 0.1)) def update_waves(v_p=6.8, v_s=3.9): # 重新初始化模拟参数 simulate(v_p, v_s)

通过调整参数可以直观发现：当横波速度降为零时（模拟流体情况），S波立即消失，而P波仍能传播——这正是海上勘探主要利用P波的原因。