有限差分法求解二维声波方程：从理论到波场动画的完整实践

1. 从波动方程到波场动画：一个地球物理从业者的实践笔记

搞地球物理勘探的同行，或者对地震波传播数值模拟感兴趣的朋友，对“波场模拟”这个词应该不陌生。我们每天处理的叠前地震数据，本质上就是地下介质对震源激发波场的响应记录。而想要真正理解这些复杂波形背后的物理机制，或者验证一个新的反演算法，最直接有效的方法就是自己动手“造”一个波场出来。这就是数值模拟的核心价值——它像一个数字沙盘，让我们能在计算机里“看到”波是如何在地下传播、反射、折射和衰减的。

今天我想分享的，就是构建这个数字沙盘最经典、也最基础的一环：使用有限差分法（FDM）求解二维声波波动方程，并生成直观的波场传播动图。别看“各向同性介质”和“二维声波”听起来像是做了很多简化，但这恰恰是理解更复杂模型（如各向异性、粘弹性、三维）的绝佳起点。很多初学者一上来就想跑复杂模型，结果连基本的数值频散都控制不住，波场图看起来像一团乱麻。我从十多年前学生时代开始接触FDM，踩过无数坑，也用它解决过不少实际科研和生产中的正演问题。这篇文章，我就结合自己的经验，把从方程离散化到代码实现，再到最终生成清晰动图的完整流程和核心细节掰开揉碎讲清楚。我们的目标不仅是“跑通”代码，更是要理解每一个参数、每一步操作背后的物理意义和数值考量，最终能生成一张稳定、准确、能真实反映波传播物理过程的动态图像。

2. 理论基石：二维声波方程及其有限差分离散化

在开始写代码之前，我们必须把理论基础打牢。这一步决定了整个模拟的物理正确性和数值稳定性，绝不能含糊。

2.1 物理方程：各向同性介质中的声波近似

我们模拟的物理场景是：在一个二维空间（例如X-Z剖面）中，介质是均匀且各向同性的，即波速在各个方向相同。同时，我们采用声学近似，忽略介质的剪切刚度，只考虑纵波（P波）的传播。在这种情况下，控制波传播的偏微分方程是经典的二阶声波方程：

[ \frac{1}{v^2(x, z)} \frac{\partial^2 p}{\partial t^2} = \frac{\partial^2 p}{\partial x^2} + \frac{\partial^2 p}{\partial z^2} + s(t, x_s, z_s) ]

这里，p(x, z, t)是我们要求解的波场（通常是压力场），v(x, z)是介质的纵波速度，s(t, x_s, z_s)是位于(x_s, z_s)位置的震源项。这个方程描述了压力扰动p在空间中随时间的演化。

注意：这里使用的是压力-速度形式的声波方程，而非位移形式。在油气勘探的地震正演中，压力场是更常被观测和使用的物理量。同时，方程右边包含了震源项s，这通常用一个时间函数（如雷克子波）与空间狄拉克函数的乘积来表示，用于在特定位置激发振动。

2.2 数值离散化：有限差分法的核心思想

有限差分法的精髓，就是用离散网格点上的函数值之差，来近似表示连续函数的导数。我们要在计算机里求解，就必须把连续的时空(x, z, t)离散化。

1. 空间离散：将二维区域划分为均匀的网格。设dx和dz分别为 x 和 z 方向的网格间距，nx和nz为网格点数。那么网格点(i, j)对应的物理位置是(i*dx, j*dz)，该点的波场值记为p[i, j]。

2. 时间离散：将总模拟时间T以时间步长dt进行分割，得到nt = T/dt个时间步。第n个时间步的时刻为n*dt，该时刻的波场记为p^n[i, j]。

3. 导数近似：这是最关键的一步。我们采用最常用的二阶中心差分格式来近似方程中的二阶偏导数。

   • 对时间的二阶导数： [ \frac{\partial^2 p}{\partial t^2} \approx \frac{p^{n+1}[i, j] - 2p^n[i, j] + p^{n-1}[i, j]}{dt^2} ]
   • 对空间的二阶导数（以 x 方向为例）： [ \frac{\partial^2 p}{\partial x^2} \approx \frac{p^n[i+1, j] - 2p^n[i, j] + p^n[i-1, j]}{dx^2} ] z 方向同理。

将上述近似代入连续的波动方程，我们就能得到关于离散波场值p^{n+1}[i, j]的递推公式（也称为更新公式）：

[ p^{n+1}[i, j] = 2p^n[i, j] - p^{n-1}[i, j] + \frac{v[i, j]^2 dt^2}{dx^2} \left( p^n[i+1, j] - 4p^n[i, j] + p^n[i-1, j] + p^n[i, j+1] + p^n[i, j-1] \right) + s^n[i, j] dt^2 ]

这里我假设了dx = dz，所以分母统一用了dx^2。这个公式就是我们代码迭代的核心：已知当前时刻n和前一时刻n-1的整个空间波场，我们可以直接计算出下一时刻n+1每个网格点的波场值。这种显式的时间推进格式计算效率非常高。

2.3 稳定性条件：CFL准则

显式格式的“阿喀琉斯之踵”是稳定性。如果时间步长dt选得太大，计算会迅速发散，结果毫无意义。稳定性由著名的CFL（Courant-Friedrichs-Lewy）条件约束。对于我们的二维声波方程，使用二阶差分格式时，稳定性要求为：

[ v_{max} \cdot dt \cdot \sqrt{\frac{1}{dx^2} + \frac{1}{dz^2}} \leq C_{max} ]

其中v_max是模型中的最大波速，C_max是一个常数，对于这种中心差分格式，通常取C_max = 1 / \sqrt{2} \approx 0.707。当dx = dz时，条件简化为：

[ dt \leq \frac{dx}{v_{max} \cdot \sqrt{2}} ]

实操心得：在实际编程中，我通常会取一个安全系数，比如dt = 0.8 * dx / (v_max * sqrt(2))。这为数值误差留出了余量，确保长时间模拟也不会失稳。记住，v_max是你的模型全局最大速度，如果模型中有高速层（如盐体），必须用它来计算dt。

3. 从公式到代码：Python实现的关键步骤与技巧

理论清晰后，我们就可以用代码将其实现。我习惯用Python（NumPy）来做原型开发和教学，因为它语法简洁，可视化方便。下面我将分步拆解代码实现，并穿插我积累的一些关键技巧。

3.1 环境与模型参数设置

首先，定义模拟的全局参数。这部分虽然基础，但参数设置不合理会直接导致模拟失败或结果失真。

import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation # ========== 模拟参数 ========== nx, nz = 401, 201 # 网格点数 (x, z方向) dx, dz = 10.0, 10.0 # 网格间距 (米) dt = 0.001 # 时间步长 (秒) nt = 1000 # 总时间步数 T = nt * dt # 总模拟时间 (秒) # ========== 速度模型 ========== # 创建一个简单的层状速度模型作为例子 v = np.ones((nz, nx)) * 2000.0 # 背景速度 2000 m/s v[100:, :] = 2500.0 # 在深度100格以下，速度变为2500 m/s (模拟一个界面) # 可以在此处添加更复杂的构造，如透镜体、断层等 # v[50:80, 150:250] = 1800.0 # 例如，一个低速透镜体 # ========== 震源参数 ========== f0 = 20.0 # 主频 (Hz) src_x, src_z = nx//2, 10 # 震源位置 (网格索引，置于浅层中心) src_type = 'ricker' # 震源子波类型 # ========== 稳定性检查 ========== v_max = np.max(v) cfl = v_max * dt * np.sqrt(1/dx**2 + 1/dz**2) print(f"CFL数: {cfl:.3f}") if cfl > 0.707: print(f"警告：CFL数 {cfl:.3f} > 0.707，模拟可能不稳定！建议减小dt.") # 可以自动调整dt # dt = 0.9 * 0.707 / (v_max * np.sqrt(1/dx**2 + 1/dz**2))

关键技巧1：模型初始化。对于初学者，强烈建议从一个均匀速度模型开始（v = np.ones((nz, nx)) * 2000.0）。先确保波在均匀介质中能产生完美的同心圆扩散，然后再引入速度界面或异常体。这能帮你快速判断是算法问题还是模型问题。

关键技巧2：网格与波长。一个经验法则是，每个最小波长内至少需要8-10个网格点，才能较好地抑制数值频散。最小波长λ_min = v_min / f_max，其中f_max是震源子波的最高有效频率（对于雷克子波，约为主频f0的2.5倍）。检查一下dx和dz是否满足dx < λ_min / 10。

3.2 震源子波与波场初始化

震源是模拟的“发动机”，它的设计直接影响波场特征。

def ricker_wavelet(t, f0, t0): """ 生成雷克子波。 参数： t: 时间序列 f0: 主频 (Hz) t0: 时间延迟，用于控制子波峰值出现的时间 返回： 子波振幅序列 """ # 标准的雷克子波公式 tau = np.pi * f0 * (t - t0) return (1.0 - 2.0 * tau**2) * np.exp(-tau**2) # 生成震源时间函数 t = np.arange(nt) * dt t0 = 1.0 / f0 # 通常将峰值延迟约一个主周期 src_time_func = ricker_wavelet(t, f0, t0) # 初始化波场数组 # 我们通常需要三个时间层：过去(p0)，现在(p1)，未来(p2) p0 = np.zeros((nz, nx)) # p^{n-1} p1 = np.zeros((nz, nx)) # p^{n} p2 = np.zeros((nz, nx)) # p^{n+1}

为什么用雷克子波？雷克子波是零相位子波，频谱明确，能量集中，且是地震数据处理中常用的子波模型。它比简单的正弦波或高斯脉冲更接近实际震源。t0的引入是为了让子波在t=0时从零开始，更符合物理实际。

波场数组的“三明治”结构：这是实现时间递推的经典方法。在每一个时间步，我们用p0（n-1）和p1（n）计算p2（n+1），然后滚动更新：p0, p1, p2 = p1, p2, p0。这样只需要三个数组在内存中循环，节省了大量空间，尤其对于大规模三维模拟至关重要。

3.3 核心迭代循环与边界处理

这是整个模拟的“心脏”。我们需要在循环中完成波场更新、震源注入，并处理边界。

# 预计算系数矩阵，避免在循环中重复计算，提升效率 c = (v**2) * (dt**2) / (dx**2) # 注意：这里假设dx=dz # 用于存储每一帧波场，用于后续生成动图 (每隔一定步数存储一次) snapshot_interval = 5 # 每5个时间步存一帧 snapshots = [] for it in range(nt): # --- 1. 应用核心有限差分更新公式 (内部区域) --- # 使用数组切片操作，避免低效的Python循环。这是性能关键！ # 更新内部网格点 (i=1:nx-2, j=1:nz-2) p2[1:-1, 1:-1] = (2 * p1[1:-1, 1:-1] - p0[1:-1, 1:-1] + c[1:-1, 1:-1] * (p1[1:-2, 1:-1] + p1[2:, 1:-1] + p1[1:-1, 1:-2] + p1[1:-1, 2:] - 4 * p1[1:-1, 1:-1])) # --- 2. 注入震源 --- # 在震源点位置添加震源项。注意乘以dt^2已在系数c中体现，这里只需加子波幅值。 src_val = src_time_func[it] p2[src_z, src_x] += src_val # 简单点源注入 # 更精确的做法是考虑震源的空间分布，例如使用高斯分布平滑注入到周围几个点 # --- 3. 边界条件处理 --- # 最简单的：Dirichlet边界（固定边界），直接将边界点设为零。 # 但这会产生强烈的虚假反射。下面介绍吸收边界条件(ABC)的一种简单实现。 # 我们这里先使用最简单的固定边界，后续再讨论吸收边界。 p2[0, :] = 0.0 # 上边界 (地表) p2[-1, :] = 0.0 # 下边界 p2[:, 0] = 0.0 # 左边界 p2[:, -1] = 0.0 # 右边界 # --- 4. 存储快照 (用于动画) --- if it % snapshot_interval == 0: # 注意：存储副本，而非引用 snapshots.append(p2.copy()) # --- 5. 滚动更新时间层 --- p0, p1, p2 = p1, p2, p0 print("时间迭代完成！")

性能关键：向量化操作。注意，波场更新公式p2[1:-1, 1:-1] = ...完全使用了NumPy的数组切片和广播机制，没有使用Python的for循环。这是将计算从慢速的Python解释器转移到快速的C/Fortran底层库的关键，性能可能有数百倍的提升。初学者最容易犯的错误就是用嵌套循环去更新每个(i, j)点。

边界条件：模拟的“隐形杀手”。固定边界（p=0）会像一堵坚硬的墙，将传播到边界的波完全反射回模型内部，严重干扰有效波场。对于波场模拟动图，这种反射是致命的，会让画面充满杂乱的回波。因此，吸收边界条件（Absorbing Boundary Condition, ABC）或完美匹配层（PML）是必备的。上面代码中我故意用了固定边界，是为了让大家先看到问题。下面我们马上来改进它。

3.4 实现简单的吸收边界条件

一个简单有效的吸收边界是衰减边界。在边界附近的一定层数内，对波场施加一个衰减系数，使波逐渐减弱至零。

# 在初始化参数部分增加 absorb_width = 30 # 吸收边界宽度（网格点数） absorb_coeff = 0.99 # 每时间步的衰减系数（小于1） # 在核心迭代循环中，替换掉原来的固定边界代码块： # --- 3. 吸收边界条件 (衰减型) --- # 上边界吸收层 for i_abs in range(absorb_width): coeff = absorb_coeff ** (absorb_width - i_abs) # 越靠近边界衰减越强 p2[i_abs, :] *= coeff # 下边界 for i_abs in range(absorb_width): coeff = absorb_coeff ** (absorb_width - i_abs) p2[-(i_abs+1), :] *= coeff # 左边界 for i_abs in range(absorb_width): coeff = absorb_coeff ** (absorb_width - i_abs) p2[:, i_abs] *= coeff # 右边界 for i_abs in range(absorb_width): coeff = absorb_coeff ** (absorb_width - i_abs) p2[:, -(i_abs+1)] *= coeff

这种衰减边界实现简单，对于非垂直入射的波有一定效果，但对垂直入射或掠入射的波吸收效果不佳，且可能会引起数值反射。对于高质量的动图和生产级的模拟，我强烈建议实现PML。PML通过在边界区域引入复数坐标拉伸，使波在进入该区域后指数衰减，理论上可以实现近乎完美的吸收。虽然PML实现更复杂（需要分裂波场并引入额外的记忆变量），但网上有许多开源实现（如devito框架中的PML）可以参考。对于本文的入门目标，衰减边界在模型足够大、边界反射尚未到达主要观测区域时，已经可以生成不错的动图了。

4. 波场可视化：生成清晰动图的科学与艺术

模拟出的数据是三维数组（两个空间维，一个时间维），将其转化为直观的动图，是交流和展示成果的关键。这里有很多细节决定了动图是“专业”还是“业余”。

4.1 静态快照与动态图生成

我们先绘制几个关键时刻的静态波场快照，检查模拟是否正常。

# 绘制几个时间步的波场快照 fig, axes = plt.subplots(2, 3, figsize=(15, 8)) time_indices = [0, nt//4, nt//2, 3*nt//4, nt-1] plot_snapshots = [snapshots[i] for i in [0, len(snapshots)//4, len(snapshots)//2, 3*len(snapshots)//4, -1]] for idx, (ax, snapshot) in enumerate(zip(axes.flat, plot_snapshots)): im = ax.imshow(snapshot, cmap='seismic', aspect='auto', extent=[0, nx*dx/1000, nz*dz/1000, 0], # 转换为公里，深度向下为正 vmin=-np.max(np.abs(snapshot))*0.1, # 动态调整色标范围，突出波前 vmax=np.max(np.abs(snapshot))*0.1) ax.scatter(src_x*dx/1000, src_z*dz/1000, c='yellow', s=50, marker='*', label='Source') ax.set_xlabel('Distance (km)') ax.set_ylabel('Depth (km)') ax.set_title(f'Time Step ~{time_indices[idx]*dt:.2f}s') ax.legend() plt.colorbar(im, ax=ax, label='Pressure') plt.tight_layout() plt.show()

色标（Colormap）的选择：seismic是地震数据可视化的标准色标，中间白色代表零值，两端的红色和蓝色分别代表正负振幅，非常符合人的直觉。vmin和vmax的设置很重要，我通常设为全局最大振幅的一个比例（如10%），这样可以压制强振幅，让微弱的波前和反射波更清晰。如果直接用snapshot的绝对最大最小值，强震源附近的振幅会淹没所有细节。

接下来，是生成动图的核心。

# 生成波场传播动图 fig, ax = plt.subplots(figsize=(10, 6)) # 初始化图像对象。使用第一个快照来确定全局色标范围，保持动图颜色一致。 vmax = np.max(np.abs(snapshots[0])) * 0.2 # 设置一个合适的固定范围 im = ax.imshow(snapshots[0], cmap='seismic', aspect='auto', extent=[0, nx*dx/1000, nz*dz/1000, 0], vmin=-vmax, vmax=vmax) ax.scatter(src_x*dx/1000, src_z*dz/1000, c='yellow', s=100, marker='*', edgecolors='black', label='震源') ax.set_xlabel('水平距离 (km)') ax.set_ylabel('深度 (km)') ax.set_title('二维声波波场传播模拟 (各向同性介质)') plt.colorbar(im, ax=ax, label='压力场振幅') ax.legend(loc='upper right') ax.grid(True, linestyle='--', alpha=0.3) # 动态更新函数 def update(frame): im.set_array(snapshots[frame]) ax.set_title(f'二维声波波场传播模拟 | 时间: {frame * snapshot_interval * dt:.2f} s') return [im] # 创建动画对象 ani = FuncAnimation(fig, update, frames=len(snapshots), interval=50, blit=True) # interval控制帧间隔(ms) # 保存为GIF或MP4文件 print("正在生成动画，这可能需要一些时间...") ani.save('wavefield_propagation.gif', writer='pillow', fps=20, dpi=150) # 保存为GIF # 如需更高清，可保存为MP4（需要安装ffmpeg） # ani.save('wavefield_propagation.mp4', writer='ffmpeg', fps=20, dpi=150) print("动画已保存为 'wavefield_propagation.gif'") plt.close(fig) # 关闭图形，避免重复显示

动图参数调优：

• interval=50：控制动画播放时每帧的间隔（毫秒）。50ms对应约20帧/秒（FPS），比较流畅。
• fps=20：保存文件时的帧率。GIF一般20fps足够，MP4可以更高。
• dpi=150：输出分辨率。对于博客或演示，150dpi在清晰度和文件大小间取得平衡。
• 保持色标一致：动图中所有帧必须使用相同的vmin/vmax，否则颜色会闪烁，干扰观察。这里用第一帧的振幅来设定全局范围。

4.2 高级可视化技巧：突出物理现象

一张好的波场动图，不仅要“能动”，更要能清晰地展示物理过程。以下是一些进阶技巧：

1. 叠加速度模型轮廓：在波场图上以半透明等高线或颜色填充的方式叠加速度模型，可以一目了然地看到波前在速度界面处的变化（反射、折射）。

   # 在创建im后添加 ax.contour(v, levels=[2100], colors='gray', linewidths=1, alpha=0.7, extent=[0, nx*dx/1000, nz*dz/1000, 0]) # 画出速度2100 m/s的轮廓线

2. 绘制射线路径或波前标记：对于简单的层状模型，可以计算并绘制理论射线路径或波前时刻图，与数值结果对比，验证模拟的准确性。

   # 计算并绘制从震源出发的直达波理论走时曲线（均匀介质） # 这里只是一个示意，实际需要根据速度模型计算 # theta = np.linspace(0, 2*np.pi, 100) # r = v[src_z, src_x] * current_time # x_ray = src_x*dx/1000 + r*np.cos(theta) # z_ray = src_z*dz/1000 + r*np.sin(theta) # ax.plot(x_ray, z_ray, 'k--', linewidth=0.5, alpha=0.5)

3. 多视图对比：创建子图，同时显示波场快照、对应的速度模型、以及某一测线（如地表）的地震记录（单道或多道），信息量更丰富。

   fig, axes = plt.subplots(1, 3, figsize=(18, 5)) # 左图：速度模型 im1 = axes[0].imshow(v, cmap='viridis', aspect='auto', extent=...) # 中图：波场快照 im2 = axes[1].imshow(snapshot, cmap='seismic', aspect='auto', extent=...) # 右图：地表接收记录（需要事先在循环中记录地表各点的波场时间序列） # axes[2].imshow(seismogram, aspect='auto', extent=...， cmap='seismic') # seismogram是 (nt, nx) 的数组

5. 常见问题排查与模型设计进阶

即使代码逻辑正确，第一次运行也很可能得不到理想的波场图。下面是我总结的几个最常见的问题及其解决方法。

5.1 数值频散：波场图中的“锯齿”与“毛刺”

现象：波前本应是光滑的圆弧，但在模拟中出现了锯齿状、网格状的图案，或者高频成分传播速度变慢，导致波包散开。原因：网格不够精细，无法分辨波的最小波长。这是有限差分法固有的误差，源于用有限精度的差分近似导数。解决方案：

1. 加密网格：这是最根本的方法。确保dx和dz小于v_min / (G * f_max)，其中G是每个波长所需的网格点数，对于二阶差分，G至少取10-15。对于主频f0=20Hz，v_min=2000m/s，f_max≈50Hz，则dx < 2000/(10*50) = 4米。我们之前设的dx=10米可能就偏大了。
2. 使用高阶差分格式：将空间二阶差分（使用3个点）升级到四阶（使用5个点）或更高阶。高阶格式在相同网格下能更精确地近似导数，显著抑制频散。代价是计算量稍增，边界处理更复杂。
3. 降低震源主频：如果研究目标允许，使用更低主频的震源（如f0=10Hz），可以增大最小波长，从而在相同网格下满足采样要求。

5.2 边界反射干扰有效波场

现象：在模拟中后期，模型边界出现明显的同心圆状波纹，并向内传播，与真实的反射波混杂。原因：边界条件吸收效果不佳。解决方案：

1. 增大吸收层宽度和优化衰减系数：将absorb_width增加到50甚至100，并精细调整absorb_coeff。可以尝试使用非均匀衰减系数，例如使用二次或指数衰减函数。
2. 实现PML：如前所述，这是工业标准和学术研究的首选。虽然编码复杂，但一旦实现，可以一劳永逸地解决绝大多数边界反射问题。建议寻找成熟的代码模块进行集成。
3. 扩大模型尺寸：在主要研究区域和物理边界之间设置足够大的“缓冲区域”。让边界反射需要很长时间才能传播到关注区域，这样在感兴趣的模拟时间内，边界反射尚未到达。这是最省事但最耗内存的方法。

5.3 震源注入引起的数值噪声

现象：在震源点附近出现高频的“噪声环”，或者整个波场出现不期望的对称模式。原因：将点源近似为一个网格点上的狄拉克函数，会引入高频成分，这些高频成分更容易产生数值频散。解决方案：

1. 震源平滑：不将能量注入单个点，而是注入到一个小的空间区域（如3x3网格），并赋予其一个空间分布（如高斯分布）。这相当于对震源进行了空间低通滤波。

   src_radius = 2 for iz in range(src_z-src_radius, src_z+src_radius+1): for ix in range(src_x-src_radius, src_x+src_radius+1): dist = np.sqrt((iz-src_z)**2 + (ix-src_x)**2) if dist <= src_radius: weight = np.exp(-(dist**2)/(2*(src_radius/2)**2)) # 高斯权重 p2[iz, ix] += src_val * weight

2. 使用更光滑的震源子波：雷克子波本身是光滑的。避免使用方波、尖脉冲等包含丰富高频成分的子波。

5.4 设计有意义的模型：从均匀介质到复杂构造

当基础代码稳定后，就可以通过设计不同的速度模型v[x, z]来研究各种地质现象。

1. 水平层状模型：如上文示例，研究波在界面上的反射和透射。可以计算反射系数，与Zoeppritz方程的理论解对比。
2. 倾斜界面/断层模型：

   v = np.ones((nz, nx)) * 2000.0 # 创建一个倾斜的断层/界面 for iz in range(nz): ix_boundary = int(0.3*nx + 0.2*iz) # 倾斜的界面 v[iz, ix_boundary:] = 3000.0 # 界面右侧速度更高

   观察波在倾斜界面的反射波、透射波以及可能产生的绕射波。
3. 高速透镜体/低速异常体：

   # 在背景速度中嵌入一个高速透镜体 v = np.ones((nz, nx)) * 2500.0 v[80:120, 150:250] = 3500.0 # 高速体 # 或者一个低速空洞 v[80:120, 150:250] = 1800.0 # 低速体

   观察波的聚焦（高速体）或散射（低速体）现象。
4. 起伏地表模型：将模型上边界（iz=0）设置为非水平，并相应地调整边界条件，可以模拟地形对波传播的影响。

通过这些模型实验，你可以直观地“看到”地震勘探中遇到的各种波现象，这对于理解地震数据剖面、验证偏移成像算法、甚至向非专业人士解释地球物理概念都极具价值。生成这些不同模型的动图并对比，本身就是一份极好的研究笔记或教学材料。