用Python动态模拟地震波从代码实现理解P波与S波本质当第一次接触地震波理论时那些关于P波、S波的抽象描述总让人困惑——为什么纵波能在液体中传播而横波不能为什么面波破坏力更强传统教材的静态图示难以展现波动的动态过程而这正是编程可视化能够突破的认知边界。本文将用不到100行Python代码带您构建一个可交互的地震波模拟器在动态图形中直观理解波的本质差异。1. 环境配置与基础波动模型1.1 搭建Python科学计算环境推荐使用Anaconda创建专属虚拟环境避免依赖冲突conda create -n seismic python3.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参数说明变量物理意义典型取值c波速P波:5-8km/sdt时间步长0.001sdx空间步长0.1m注意Courant条件要求c*dt/dx ≤ 1以保证数值稳定性2. P波与S波的动态模拟对比2.1 纵波P波特性实现P波模拟关键点在于质点振动与传播方向平行。以下代码创建压缩-膨胀交替的波动效果def simulate_p_wave(steps500): # 初始化位移场 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), scale20, colorr)2.2 横波S波的特殊约束S波的介质刚性要求体现在代码中需要设置剪切模量参数添加固体介质检测逻辑if medium fluid: raise ValueError(横波不能在流体中传播)典型参数对比波型传播介质速度公式典型速度P波固/液/气√((K4G/3)/ρ)6.8km/sS波仅固体√(G/ρ)3.9km/s实验尝试修改介质密度ρ观察波速变化规律3. 复杂波场现象模拟3.1 莫霍面波速跃变模拟通过设置速度突变边界模拟地壳-地幔界面# 创建速度模型 c np.where(x NX//2, 6.8, 8.1) # 莫霍面位置关键现象观察波到达界面时部分能量反射透射波发生折射Snell定律使用能量衰减公式模拟实际传播u_next * np.exp(-0.01*x) # 指数衰减3.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(freq25, length0.1) # 4. 褶积运算 return np.convolve(rc, wavelet, modesame)4.2 波场快照技术实现多时刻波场状态保存与回放from matplotlib.animation import FuncAnimation def update(frame): line.set_ydata(snapshots[frame]) return line, ani FuncAnimation(fig, update, frameslen(snapshots), interval50, blitTrue)交互功能增强添加滑块控制波速参数实现介质类型实时切换支持波前标记与传播计时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_p6.8, v_s3.9): # 重新初始化模拟参数 simulate(v_p, v_s)通过调整参数可以直观发现当横波速度降为零时模拟流体情况S波立即消失而P波仍能传播——这正是海上勘探主要利用P波的原因。
别再死记硬背P波S波了！用Python模拟地震波传播，直观理解勘探原理
用Python动态模拟地震波从代码实现理解P波与S波本质当第一次接触地震波理论时那些关于P波、S波的抽象描述总让人困惑——为什么纵波能在液体中传播而横波不能为什么面波破坏力更强传统教材的静态图示难以展现波动的动态过程而这正是编程可视化能够突破的认知边界。本文将用不到100行Python代码带您构建一个可交互的地震波模拟器在动态图形中直观理解波的本质差异。1. 环境配置与基础波动模型1.1 搭建Python科学计算环境推荐使用Anaconda创建专属虚拟环境避免依赖冲突conda create -n seismic python3.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参数说明变量物理意义典型取值c波速P波:5-8km/sdt时间步长0.001sdx空间步长0.1m注意Courant条件要求c*dt/dx ≤ 1以保证数值稳定性2. P波与S波的动态模拟对比2.1 纵波P波特性实现P波模拟关键点在于质点振动与传播方向平行。以下代码创建压缩-膨胀交替的波动效果def simulate_p_wave(steps500): # 初始化位移场 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), scale20, colorr)2.2 横波S波的特殊约束S波的介质刚性要求体现在代码中需要设置剪切模量参数添加固体介质检测逻辑if medium fluid: raise ValueError(横波不能在流体中传播)典型参数对比波型传播介质速度公式典型速度P波固/液/气√((K4G/3)/ρ)6.8km/sS波仅固体√(G/ρ)3.9km/s实验尝试修改介质密度ρ观察波速变化规律3. 复杂波场现象模拟3.1 莫霍面波速跃变模拟通过设置速度突变边界模拟地壳-地幔界面# 创建速度模型 c np.where(x NX//2, 6.8, 8.1) # 莫霍面位置关键现象观察波到达界面时部分能量反射透射波发生折射Snell定律使用能量衰减公式模拟实际传播u_next * np.exp(-0.01*x) # 指数衰减3.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(freq25, length0.1) # 4. 褶积运算 return np.convolve(rc, wavelet, modesame)4.2 波场快照技术实现多时刻波场状态保存与回放from matplotlib.animation import FuncAnimation def update(frame): line.set_ydata(snapshots[frame]) return line, ani FuncAnimation(fig, update, frameslen(snapshots), interval50, blitTrue)交互功能增强添加滑块控制波速参数实现介质类型实时切换支持波前标记与传播计时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_p6.8, v_s3.9): # 重新初始化模拟参数 simulate(v_p, v_s)通过调整参数可以直观发现当横波速度降为零时模拟流体情况S波立即消失而P波仍能传播——这正是海上勘探主要利用P波的原因。
