别再死记硬背了！用Python模拟地震波传播，5分钟搞懂反射波与折射波

weixin_30292745

582人浏览 · 2026-05-26 11:48:51

weixin_30292745 · 2026-05-26 11:48:51 发布

用Python模拟地震波：5行代码可视化反射与折射原理

地震勘探听起来像是地质学家的专利？其实只需几行Python代码，你就能在屏幕上看到地震波如何在地下舞蹈。作为曾经被"斯奈尔定律"折磨过的程序员，我发现用代码重现波前传播后，那些抽象概念突然变得触手可及。本文将带你用NumPy和Matplotlib，从零构建一个微型地震模拟器——不需要任何地球物理背景，只要会写 import numpy 就能开始这场视觉盛宴。

1. 准备工作：搭建数字实验室

在开始模拟前，我们需要配置一个轻量级计算环境。推荐使用Jupyter Notebook交互式开发，每一步都能即时看到波形变化。

基础工具包安装 ：

pip install numpy matplotlib ipywidgets

核心库的作用很简单：

NumPy ：处理波场数据的矩阵运算
Matplotlib ：将数字矩阵转化为动态图像
ipywidgets （可选）：创建交互控件实时调整参数

定义介质属性时，地球物理学家用 波阻抗 （密度与波速的乘积）描述地层特性。我们可以用字典模拟两种介质的交界：

medium = {
    'upper': {'velocity': 1500, 'density': 2.0},  # 水层典型参数
    'lower': {'velocity': 3000, 'density': 2.5}   # 砂岩典型参数
}

注意：实际勘探中速度单位是m/s，密度单位是g/cm³。模拟时保持单位统一即可，比例关系更重要。

2. 生成地震子波：数字震源诞生记

真实地震勘探使用炸药或振动车产生震波，我们则用数学函数合成 雷克子波 （Ricker Wavelet）——最接近实际爆炸产生的地震脉冲。

def ricker_wavelet(freq, duration, dt):
    """生成雷克子波"""
    t = np.arange(-duration/2, duration/2, dt)
    return (1 - 2*(np.pi*freq*t)**2) * np.exp(-(np.pi*freq*t)**2)

# 示例：主频30Hz的子波
wavelet = ricker_wavelet(freq=30, duration=0.1, dt=0.001)
plt.plot(wavelet); plt.title("雷克子波波形"); plt.show()

关键参数对比表：

参数	物理意义	典型值范围	模拟影响
freq	主频	10-100Hz	高频分辨率高但衰减快
duration	持续时间	0.05-0.2秒	影响波形的周期数
dt	时间间隔	0.0001-0.001秒	数值稳定性关键

3. 波场传播：惠更斯原理的代码实现

惠更斯原理说"波前每点都是新震源"，用代码表达就是卷积运算。我们构建一个2D波场矩阵，每个时间步根据相邻点计算新状态：

def propagate(wavefield, velocity, dt, dx):
    """二维波场传播核心算法"""
    next_wave = 2*wavefield[1] - wavefield[0] + 
               (velocity*dt/dx)**2 * np.diff(wavefield[1], 2)
    return np.vstack([wavefield[1], next_wave])

模拟反射/折射现象时，需要在介质交界处处理边界条件。 斯奈尔定律 告诉我们波在界面处的"分裂"规律：

def snells_law(theta_i, v1, v2):
    """计算反射角和折射角"""
    theta_r = theta_i  # 反射角等于入射角
    theta_t = np.arcsin(v2/v1 * np.sin(theta_i))  # 折射角
    return theta_r, theta_t

运行完整模拟后，你会看到这样的波前演化过程：

初始圆形波前从震源向外扩展
遇到界面时部分能量反射回上层介质
剩余能量以新角度折射进入下层
折射波速度突变导致波前弯折

4. 可视化技巧：让波动物理现形

静态图像难以展现波动过程，我们使用Matplotlib的动画功能：

from matplotlib.animation import FuncAnimation

def update(frame):
    """动画帧更新函数"""
    global wavefield
    wavefield = propagate(wavefield, velocity, dt, dx)
    im.set_array(wavefield[1])
    return [im]

ani = FuncAnimation(fig, update, frames=100, blit=True)
HTML(ani.to_jshtml())  # 在Jupyter中嵌入动画

进阶技巧：添加交互控件实时调整参数

from ipywidgets import interact

@interact(v1=(1000, 5000), v2=(1000, 5000), angle=(0, 90))
def interactive_simulation(v1=1500, v2=3000, angle=30):
    """交互式参数调整"""
    theta_r, theta_t = snells_law(np.radians(angle), v1, v2)
    print(f"反射角：{np.degrees(theta_r):.1f}°  折射角：{np.degrees(theta_t):.1f}°")
    # 重新运行模拟...

5. 从模拟到实践：解读地震剖面

将多个接收点的波形按位置排列，就得到地震勘探的核心数据—— 地震剖面 。用Python可以模拟检波器阵列的响应：

receivers = np.linspace(0, 1000, 20)  # 20个均匀分布的检波器
traces = []
for x in receivers:
    trace = simulate_receiver(x, source_pos=500)
    traces.append(trace)

# 绘制地震道集
plt.figure(figsize=(10,6))
plt.imshow(np.array(traces).T, aspect='auto', 
          cmap='seismic', vmin=-0.5, vmax=0.5)
plt.colorbar(label='振幅'); plt.title("合成地震记录");

常见波型识别特征：