保姆级教程：用Python+NumPy从零模拟地震波传播（附代码）

用Python+NumPy构建地震波传播模拟器：从理论到可视化实战

地震波模拟是地球物理研究和工程勘探的核心技术之一。想象一下，当你运行几行代码就能看到地震波在地壳中传播的动态过程，这种将抽象理论转化为可视化结果的体验，正是本文要带给你的。我们将完全从零开始，使用Python和NumPy构建一个功能完整的地震波模拟器，避开繁琐的数学推导，直接进入代码实现的快车道。

1. 环境准备与基础概念

在开始编码前，我们需要配置合适的开发环境。推荐使用Anaconda创建专属的Python 3.8+环境：

conda create -n seismic python=3.8 numpy matplotlib scipy jupyter conda activate seismic

地震波主要分为三种类型，它们在模拟中需要区别对待：

提示：在实际模拟中，P波总是最先到达，这个特性可以用来验证代码的正确性。

理解波动方程是模拟的基础。我们采用简化的二维声波方程作为起点：

import numpy as np def acoustic_2d(dt, dx, dz, velocity): """ 二维声波方程有限差分核心 dt: 时间步长 dx, dz: 空间步长 velocity: 速度模型矩阵 """ # 初始化波场 next_wave = np.zeros_like(velocity) # 计算空间二阶差分 next_wave[1:-1,1:-1] = (velocity[1:-1,1:-1]*dt)**2 * ( (wave[2:,1:-1] - 2*wave[1:-1,1:-1] + wave[:-2,1:-1])/dx**2 + (wave[1:-1,2:] - 2*wave[1:-1,1:-1] + wave[1:-1,:-2])/dz**2) return next_wave

2. 构建地震子波发生器

地震子波是模拟的"震源"，雷克子波(Ricker Wavelet)因其数学简洁和物理意义明确而成为首选。下面实现一个可调节主频和相位的子波生成器：

def ricker_wavelet(freq, dt, length, phase=0): """ 生成雷克子波 freq: 主频(Hz) dt: 时间采样间隔(s) length: 子波长度(点数) phase: 相位调整(0-360) """ t = np.linspace(-length/2, length/2, length) t = t * dt wavelet = (1 - 2*(np.pi*freq*t)**2) * np.exp(-(np.pi*freq*t)**2) # 相位调整 if phase != 0: spectrum = np.fft.fft(wavelet) angles = np.angle(spectrum) + np.deg2rad(phase) wavelet = np.real(np.fft.ifft(np.abs(spectrum)*np.exp(1j*angles))) return wavelet

不同子波特性对比：

• 主频影响：高频子波分辨率高但穿透力弱
• 相位特性：
  • 零相位：能量对称分布（理论模型）
  • 最小相位：能量前部集中（实际勘探）
  • 混合相位：能量中部集中

可视化子波（Jupyter中使用）：

import matplotlib.pyplot as plt freqs = [10, 20, 30] wavelets = [ricker_wavelet(f, 0.001, 200) for f in freqs] plt.figure(figsize=(10,6)) for w, f in zip(wavelets, freqs): plt.plot(w, label=f'{f}Hz') plt.legend(); plt.title('不同主频雷克子波对比'); plt.show()

3. 介质建模与波场传播

真实地下介质非常复杂，我们先从简单的分层模型开始构建速度场：

def build_velocity_model(size, layers): """ 构建分层速度模型 size: (width, depth) 模型尺寸 layers: [(depth_top, vp, vs), ...] 各层参数 """ vp_model = np.ones(size) * layers[0][1] # P波速度 vs_model = np.ones(size) * layers[0][2] # S波速度 for i in range(1, len(layers)): depth, vp, vs = layers[i] vp_model[:, depth:] = vp vs_model[:, depth:] = vs return vp_model, vs_model

典型岩石速度参考值：

实现完整的波场更新循环：

def seismic_simulation(model_size, velocity, source_pos, receiver_pos, dt=0.001, steps=1000): """ 完整的地震波模拟流程 """ # 初始化波场和记录器 wavefield = np.zeros(model_size) records = np.zeros((steps, len(receiver_pos))) # 生成震源子波 source = ricker_wavelet(20, dt, 100) for step in range(steps): # 注入震源 if step < len(source): wavefield[source_pos] += source[step] # 更新波场 wavefield = acoustic_2d(dt, 10, 10, velocity) # 边界吸收（防止反射） wavefield = apply_absorbing_bc(wavefield) # 记录接收器数据 for i, pos in enumerate(receiver_pos): records[step, i] = wavefield[pos] return records

注意：实际应用中需要加入稳定性条件检查，确保CFL条件满足：max(v)*dt/dx < 0.7

4. 高级技巧与实战优化

经过基础模拟后，我们需要解决几个关键问题才能获得逼真的结果：

边界反射处理使用完美匹配层(PML)吸收边界：

def apply_pml(wavefield, pml_width=20, damping=0.1): """ 应用PML吸收边界 """ for i in range(pml_width): # 左边界 wavefield[:, i] *= np.exp(-damping*(pml_width-i)) # 右边界 wavefield[:, -i-1] *= np.exp(-damping*(pml_width-i)) # 上边界 wavefield[i, :] *= np.exp(-damping*(pml_width-i)) # 下边界 wavefield[-i-1, :] *= np.exp(-damping*(pml_width-i)) return wavefield

各向异性处理修改波动方程加入各向异性参数：

def anisotropic_update(wave, vp, vs, epsilon, delta, dt): """ 考虑各向异性的波场更新 epsilon, delta: Thomsen各向异性参数 """ # 计算空间导数 dxx = np.gradient(np.gradient(wave, axis=0), axis=0) dzz = np.gradient(np.gradient(wave, axis=1), axis=1) dxz = np.gradient(np.gradient(wave, axis=0), axis=1) # 各向异性修正项 term1 = (1 + 2*epsilon) * dxx term2 = dzz term3 = np.sqrt(1 + 2*delta) * dxz return wave + (vp**2 * dt**2) * (term1 + term2 + term3)

性能优化技巧使用Numba加速计算：

from numba import jit @jit(nopython=True) def wave_propagation_numba(wave, velocity, dt, dx): next_wave = np.zeros_like(wave) for i in range(1, wave.shape[0]-1): for j in range(1, wave.shape[1]-1): lap = (wave[i+1,j] - 2*wave[i,j] + wave[i-1,j])/dx**2 \ + (wave[i,j+1] - 2*wave[i,j] + wave[i,j-1])/dx**2 next_wave[i,j] = 2*wave[i,j] - next_wave[i,j] + (velocity[i,j]*dt)**2 * lap return next_wave

典型性能对比：

5. 结果可视化与分析

模拟结果的直观展示是理解波传播的关键。我们实现多维度可视化：

波场快照动画

import matplotlib.animation as animation def animate_wavefield(wavefields): fig = plt.figure(figsize=(10,8)) im = plt.imshow(wavefields[0], cmap='seismic', vmin=-0.1, vmax=0.1) def update(frame): im.set_array(wavefields[frame]) return [im] ani = animation.FuncAnimation(fig, update, frames=len(wavefields), interval=50, blit=True) plt.close() return ani

地震记录分析

def plot_seismogram(records, dt): plt.figure(figsize=(12,6)) plt.imshow(records.T, aspect='auto', cmap='seismic', extent=[0, len(records)*dt, len(records[0]), 0]) plt.colorbar(label='振幅') plt.xlabel('时间(s)'); plt.ylabel('接收器编号') plt.title('共炮点地震记录')

实际应用中的关键观察点

1. 初至识别：第一个到达的波永远是P波
2. 波型转换：当P波遇到界面时会产生转换S波
3. 振幅衰减：高频成分随距离快速衰减
4. 界面反射：振幅极性变化指示波阻抗增减

通过调整参数观察现象变化：

# 改变主频观察分辨率差异 for freq in [10, 20, 30]: source = ricker_wavelet(freq, 0.001, 100) records = seismic_simulation(..., source) plot_seismogram(records) # 改变速度模型观察反射特征 models = [ [(0, 2.0, 1.0), (100, 2.5, 1.2)], # 弱反射 [(0, 2.0, 1.0), (100, 3.0, 1.7)] # 强反射 ] for vp, vs in models: velocity = build_velocity_model(..., [(0, vp, vs)]) records = seismic_simulation(..., velocity)

在完成基础模拟后，尝试将这些技术应用到实际场景中，比如设计一个包含断层、背斜等地质构造的复杂模型，观察地震波在这些结构中的传播特征。