用Python+NumPy模拟地震波传播：从代码实现到可视化分析

地震波传播的数值模拟一直是地球物理学和地震工程领域的核心课题。传统教学中，学生往往被大量公式推导所困扰，难以直观理解波场传播的物理过程。本文将带你用Python的科学计算库NumPy和可视化工具Matplotlib，构建一个完整的地震波传播数值模型。通过代码实现，你将看到反射波、折射波如何在界面处产生，以及不同地层参数对波场的影响。

1. 环境准备与基础概念

在开始编码前，我们需要明确几个关键概念。地震波在介质中的传播遵循弹性波方程，而数值模拟的核心就是离散化这个偏微分方程。对于二维情况，我们可以使用声波近似，简化方程为：

# 二维声波方程离散形式
def wave_propagation(u, u_prev, c, dt, dx, dy):
    """
    u: 当前时刻波场
    u_prev: 前一时刻波场
    c: 介质速度场
    dt: 时间步长
    dx, dy: 空间步长
    """
    laplacian = (np.roll(u, 1, axis=0) + np.roll(u, -1, axis=0) - 2*u)/(dx**2) + \
                (np.roll(u, 1, axis=1) + np.roll(u, -1, axis=1) - 2*u)/(dy**2)
    u_next = 2*u - u_prev + (c*dt)**2 * laplacian
    return u_next

安装必要的Python库：

pip install numpy matplotlib scipy ipywidgets

关键参数说明：

参数	描述	典型值
dx, dy	空间离散间隔	5-20 m
dt	时间步长	0.001-0.01 s
nt	时间步数	1000-5000
nx, ny	空间网格数	100-500

注意：时间步长dt必须满足CFL稳定性条件：dt ≤ min(dx,dy)/√2 * cmax，其中cmax是介质中的最大波速

2. 构建地层模型与波场初始化

地震波模拟的第一步是定义地下介质模型。我们创建一个包含速度突变的简单两层模型：

def create_velocity_model(nx, ny, dx, dy):
    """创建两层速度模型"""
    v = np.ones((ny, nx)) * 1500  # 上层速度1500 m/s (水或疏松沉积层)
    v[ny//2:, :] = 3000  # 下层速度3000 m/s (致密岩石)
    return v

# 模型参数
nx, ny = 300, 300
dx, dy = 10, 10  # 10米网格
v = create_velocity_model(nx, ny, dx, dy)

# 计算稳定时间步长
dt = 0.9 * min(dx, dy) / (np.sqrt(2) * v.max())

波场初始化需要三个数组：当前时刻波场(u)、前一时刻波场(u_prev)和下一时刻波场(u_next)。我们还需要定义震源函数，通常使用Ricker子波：

def ricker_wavelet(t, f0):
    """Ricker子波(墨西哥帽小波)"""
    return (1 - 2*(np.pi*f0*t)**2) * np.exp(-(np.pi*f0*t)**2)

# 震源参数
f0 = 15  # 主频15Hz
t0 = 1.5 / f0  # 延迟时间
source_loc = (50, 150)  # 震源位置

3. 波场传播模拟实现

有了上述准备，我们可以实现完整的时间步进循环。下面是核心模拟代码：

def run_simulation(nt, nx, ny, dx, dy, dt, v, source_loc, f0, t0):
    # 初始化波场
    u = np.zeros((ny, nx))
    u_prev = np.zeros_like(u)
    u_next = np.zeros_like(u)
    
    # 存储地震记录
    seismogram = np.zeros(nt)
    receiver_loc = (250, 150)  # 检波器位置
    
    for it in range(nt):
        # 计算当前时间
        t = it * dt
        
        # 添加震源
        if it < len(seismogram):
            u[source_loc] += ricker_wavelet(t - t0, f0) * (dt**2)
        
        # 波场传播
        u_next = wave_propagation(u, u_prev, v, dt, dx, dy)
        
        # 边界条件(简单吸收边界)
        u_next[0,:] = u_next[-1,:] = u_next[:,0] = u_next[:,-1] = 0
        
        # 记录地震道
        seismogram[it] = u[receiver_loc]
        
        # 更新波场
        u_prev, u = u, u_next
        
    return seismogram, u

模拟过程中几个关键点：

• 边界处理 ：简单设置为0会导致强烈反射，实际应用中应使用吸收边界条件(PML)
• 数值频散 ：空间和时间离散不足会导致高频成分出现数值频散
• 计算效率 ：纯Python循环较慢，实际项目可考虑使用Numba加速或转向C++实现

4. 结果可视化与分析

模拟完成后，我们需要直观展示波场传播过程和地震记录特征。以下是几个关键可视化：

def plot_wavefield(u, v, cmap='seismic'):
    """绘制波场快照"""
    plt.figure(figsize=(10,8))
    plt.imshow(u, cmap=cmap, extent=[0, nx*dx, ny*dy, 0])
    plt.colorbar(label='振幅')
    plt.contour(v, levels=[2250], colors='k')  # 标记速度界面
    plt.scatter(*source_loc[::-1], c='r', marker='*', s=200, label='震源')
    plt.scatter(*receiver_loc[::-1], c='g', marker='v', s=100, label='检波器')
    plt.xlabel('距离(m)')
    plt.ylabel('深度(m)')
    plt.legend()
    plt.title('波场快照')

def plot_seismogram(seismogram, dt):
    """绘制地震记录"""
    plt.figure(figsize=(12,4))
    t = np.arange(len(seismogram)) * dt
    plt.plot(t, seismogram)
    plt.xlabel('时间(s)')
    plt.ylabel('振幅')
    plt.title('单道地震记录')
    plt.grid(True)

通过分析模拟结果，我们可以观察到：

1. 直达波 ：从震源直接传播到检波器的波
2. 反射波 ：在速度界面反射后到达检波器的波
3. 折射波 ：当下层速度大于上层时，在临界角后出现的折射波

典型地震记录特征：

波型	到达时间	振幅特征	频率特征
直达波	最早	通常最大	与震源一致
反射波	中等	取决于反射系数	可能高频衰减
折射波	最晚	通常较小	低频成分较多

5. 高级扩展与参数研究

基础模型运行成功后，我们可以进行多种扩展研究：

不同地层结构的波场响应

# 复杂速度模型示例
def create_complex_model(nx, ny):
    v = np.ones((ny, nx)) * 1800
    v[100:200, 100:200] = 2500  # 高速透镜体
    v[200:, :] = 3200  # 基底
    v[150:180, 50:250] = 2200  # 低速层
    return v

各向异性介质模拟

在横向各向同性介质中，波速随传播方向变化，需要修改波场更新公式：

def anisotropic_propagation(u, u_prev, v_h, v_v, dt, dx, dy):
    # v_h: 水平方向波速
    # v_v: 垂直方向波速
    laplacian_h = (np.roll(u, 1, axis=1) + np.roll(u, -1, axis=1) - 2*u)/(dx**2)
    laplacian_v = (np.roll(u, 1, axis=0) + np.roll(u, -1, axis=0) - 2*u)/(dy**2)
    u_next = 2*u - u_prev + dt**2 * (v_h**2 * laplacian_h + v_v**2 * laplacian_v)
    return u_next

吸收边界条件实现

# 简单的海绵吸收边界
def apply_sponge_boundary(u, sponge_width=20):
    ny, nx = u.shape
    for i in range(sponge_width):
        # 左边界
        u[:, i] *= (i / sponge_width)**2
        # 右边界
        u[:, nx-1-i] *= (i / sponge_width)**2
        # 上边界
        u[i, :] *= (i / sponge_width)**2
        # 下边界
        u[ny-1-i, :] *= (i / sponge_width)**2
    return u

6. 交互式可视化与教学应用

为了让模拟过程更加直观，我们可以创建交互式可视化工具：

from ipywidgets import interact, IntSlider

def interactive_wavefield(step):
    """交互式波场查看"""
    plt.figure(figsize=(10,8))
    plt.imshow(wavefields[step], cmap='seismic', vmin=-0.1, vmax=0.1)
    plt.colorbar()
    plt.title(f'波场传播 - 时间步 {step}')
    
interact(interactive_wavefield, step=IntSlider(min=0, max=nt-1, step=10, value=0))

教学应用建议：

1. 波前传播可视化 ：展示惠更斯原理如何在实际波场中体现
2. 反射系数实验 ：修改界面两侧波阻抗，观察反射振幅变化
3. 临界角演示 ：调整入射角度，展示折射波产生条件
4. 频散分析 ：比较不同空间离散下的数值频散现象

7. 性能优化与大规模模拟

当模型规模增大时，纯Python实现会变得缓慢。以下是一些优化策略：

使用Numba加速

from numba import jit

@jit(nopython=True)
def wave_propagation_numba(u, u_prev, c, dt, dx, dy):
    # Numba加速版本
    u_next = np.zeros_like(u)
    for i in range(1, u.shape[0]-1):
        for j in range(1, u.shape[1]-1):
            laplacian = (u[i+1,j] + u[i-1,j] - 2*u[i,j])/(dx**2) + \
                        (u[i,j+1] + u[i,j-1] - 2*u[i,j])/(dy**2)
            u_next[i,j] = 2*u[i,j] - u_prev[i,j] + (c[i,j]*dt)**2 * laplacian
    return u_next

并行计算策略

对于三维模型或长时间模拟，可以考虑：

• 使用Dask进行分布式计算
• 使用CUDA实现GPU加速
• 采用域分解方法进行MPI并行

内存优化技巧

对于超大模型，可以：

• 使用内存映射文件处理大型数组
• 采用checkpointing技术只保存关键时间步
• 降低输出采样率减少存储需求