企业官网建设流程全解析

简介：弹性波方程是地球物理学中描述地震波传播的核心方程，正演模拟则借助数值求解方法，预测地震波在介质中的传播路径、波形振幅与到达时刻。这份面向地震波数值模拟学习者与科研人员的MATLAB算法资源，利用MATLAB的矩阵运算特性，采用10阶差分精度对弹性波方程进行离散求解，相比低阶格式能更细腻地刻画高频波场细节，适合用来学习高阶有限差分方法的实现思路。压缩包内仅含1个MATLAB脚本文件，包体约2KB，代码短小精悍，便于逐行研读和修改测试；但由于未加入吸收边界条件，初次运行时更偏重展示无边界情形下的波场基本行为。目前已有133人学习，资源可作为高阶地震波模拟的入门样例，在此基础上自行添加PML或人工边界条件，即可进一步优化模拟结果在实际地质构造中的物理合理性。

1. 弹性波方程正演模拟为什么值得做：声波近似掩盖了太多真实信息

接手过实际工区正演任务的人大概都经历过这种别扭：声波方程正演跑得飞快，模拟出来单分量地震记录也像模像样，可一旦把振幅、相位谱拿去和野外观测数据对比，立刻就能看出不对劲——横波信息完全空白，AVO 响应不对，转换波更是无从谈起。问题不在代码写错了，而在物理模型本身。声波方程假设介质只传播纵波，地层内没有剪切应力，这在储层描述、裂缝预测和流体识别这些真正需要弹性参数约束的场景里，等于主动丢掉了一大半有效信息。

弹性波方程的地震波场正演模拟，简单说就是在给定速度模型、密度模型和衰减参数的条件下，用数值方法同时求解纵波和横波在介质中的传播过程，输出多分量波场快照和合成地震记录。它解决的问题很具体：帮你在没有野外采集数据之前就能验证观测系统能不能收到目标层的转换波，帮全波形反演准备敏感度核函数，帮微地震事件定位判断震源机制产生的波场特征。适合谁？做地震采集设计、VSP 资料处理、FWI 和微地震解释的工程师和研究生，只要你的工作对象是弹性波而不是纯声波，这篇文章里网格参数怎么定、边界怎么吸收、震源怎么加载这些事，都会直接落在你的日活里。

2. 弹性波方程的理论基础：从双曲型方程组到交错网格有限差分

2.1 一阶速度-应力方程：为什么数值实现几乎都选它

弹性波正演可选的方程形式有两种：二阶位移方程和一阶速度-应力方程。教科书里推导方便，多用位移表示，但真正落到数值计算，主流实现几乎都选择一阶速度-应力方程。

原因在离散误差上。二阶位移方程直接对空间求二阶导，需要至少三个网格点才能构造中心差分，而一阶速度-应力方程组里空间导数和时间导数都只到一阶，配合交错网格可以让空间差分模板只跨半个网格步长，相同网格密度下相位误差更小，各向异性伪影更少。对于含自由表面的模型，一阶形式处理应力自由边界条件也更直接。

在二维各向同性介质中，一阶速度-应力方程的标准形式是：

∂vx/∂t = (1/ρ) * (∂σxx/∂x + ∂σxz/∂z) ∂vz/∂t = (1/ρ) * (∂σxz/∂x + ∂σzz/∂z) ∂σxx/∂t = (λ+2μ) * ∂vx/∂x + λ * ∂vz/∂z ∂σzz/∂t = λ * ∂vx/∂x + (λ+2μ) * ∂vz/∂z ∂σxz/∂t = μ * (∂vx/∂z + ∂vz/∂x)

这里 vx、vz 是质点振动速度分量，σxx、σzz 是正应力，σxz 是剪切应力，ρ 是密度，λ 和 μ 是拉梅常数。实际程序里一般不直接存 λ 和 μ，而是换算成纵波速度 vp、横波速度 vs 和密度 ρ，空间任意点的 λ = ρ(vp² - 2vs²)，μ = ρvs²。

给参数做换算的时候注意一个常见误用：有人为了省事直接在空间上线性插值 vp 和 vs，这种做法会在介质分界面处产生虚假的薄层反射，因为 vp² 的线性平均并不等于平均介质中的等效速度。我一般先对 vp、vs、ρ 分别做插值，再算 λ 和 μ；如果模型本身是逐点给定的块状模型，插值顺序带来的误差在低频段不明显，但涉及高角度界面时还是严格一点好。

2.2 交错网格上的空间差分离散：精度与实现的平衡

交错网格的核心思想是把不同物理量放在半个网格步长的错位位置上。以二维笛卡尔网格为例，法向应力 σxx、σzz 定义在整网格点上，速度分量 vx 定义在 x 方向半网格偏移点，vz 定义在 z 方向半网格偏移点，剪切应力 σxz 定义在 x 和 z 同时偏移半网格的位置。这样每个空间导数在相邻半个网格点之间用中心差分构造，只需要两点模板就能达到二阶精度。

空间差分阶数是精度和计算量的博弈。四阶空间差分（每阶导数需要 4 个网格点）是性价比最高的选择，能在每波长 5～6 个网格点的情况下把数值频散压到肉眼不可见；要模拟长传播距离或做高精度全波形反演，可以升到八阶甚至十阶，但需要额外处理边界附近的单向差分。时间推进上多数实现用二阶龙格-库塔或蛙跳格式，时间步长受 CFL 条件限制，二维情况下大约要满足 vmax * dt / dx < 0.6～0.7（四阶空间差分实际经验值更紧，通常取 0.4～0.5 更稳妥）。

一个常见折衷：空间上做高阶，时间上仍保持二阶。因为时间阶数不够会在波形尾部造成轻微的非对称误差，而空间阶数不够会造成波长尺度内的振荡尾巴，两种误差的特征不同，前者更容易靠吸收边界压住，后者会在全剖面到处漏。

这一节落到代码里，核心就是构造差分系数。四阶中心差分在均匀网格上的系数是 [1/12, -2/3, 0, 2/3, -1/12]/dx，交错网格上的半网格差分系数会不一样，要按物理量的空间位置分两套系数处理。你要是直接拿整网格系数套到交错网格上，波场会立刻出现棋盘状高频振荡。

3. 跑通第一个弹性波正演模拟：模型、参数与最小代码骨架

3.1 设计一个两层介质模型：密度、纵波速度、横波速度怎么给

第一步先不追求复杂构造，用最简单的水平两层模型验证程序正确性。假定模型尺寸 2000 m × 2000 m，网格间距 dx = dz = 10 m，即 200 × 200 个网格点。上层介质 vp = 2500 m/s，vs = 1400 m/s，ρ = 2200 kg/m³，代表风化层；下层介质 vp = 4000 m/s，vs = 2300 m/s，ρ = 2400 kg/m³，代表致密砂岩。分界面放在深度 1000 m 处。

参数给定时注意泊松比的约束。在各向同性弹性介质中，vp / vs 的下限是 √2，低于这个值意味着泊松比为负，能量上不稳定；上限没有硬约束，但 vp / vs 超过 3.5 在沉积岩里是极端的，数值上会让 λ 远大于 μ，导致正应力更新方程中两项大数相消，浮点误差被放大。

模型用 NumPy 数组初始化，每个物理量一个数组，不压缩存储：

import numpy as np nx, nz = 200, 200 dx = dz = 10.0 vp = np.ones((nz, nx)) * 2500.0 vs = np.ones((nz, nx)) * 1400.0 rho = np.ones((nz, nx)) * 2200.0 vp[nz//2:, :] = 4000.0 vs[nz//2:, :] = 2300.0 rho[nz//2:, :] = 2400.0 mu = rho * vs**2 lam = rho * vp**2 - 2.0 * mu

这段代码里nz//2是数组行方向的中线，对应深度 1000 m。λ 和 μ 用vp**2而不是vp直接参与计算，目的是保持拉梅常数与速度之间的平方关系，后续更新应力时不要写成(lam + 2*mu) * dvx_dx + lam * dvz_dz以外的东西，拉梅常数算错最常见就是把2*mu的系数漏掉或者把lam乘到剪切项上。

3.2 稳定性条件与网格尺寸：为什么 CFL 数不能随便取

网格尺寸和时步的决策直接决定一个正演算例从「能跑」变成「跑得对」。

网格间距的约束来自频散控制。通常要求最短波长内至少有 5～6 个网格点，而最短波长取决于横波速度和震源最高频率：λ_min = vs_min / f_max。如果 vs_min 是 1400 m/s，震源最高频率 30 Hz，λ_min ≈ 47 m，网格间距取 10 m 大约 4.7 个点，勉强达到四阶差分的及格线；想要更好的波形保真度，要么网格加密到 7～8 m，要么把震源主频降到 20 Hz。

时间步长的约束来自 CFL 条件。对二维四阶空间差分，经验上限在 0.45 附近，这里 CFL 数的定义是 vmax * dt / dx。vmax 取整个模型的最大纵波速度 4000 m/s，dx 取 10 m 时：

vmax = np.max(vp) cfl = 0.45 dt = cfl * dx / vmax print(dt)

计算得到 dt ≈ 0.001125 s，取整为 0.001 s 安全余量更大。注意这里必须用最大纵波速度而不是平均速度或最大横波速度，因为弹性波正演里纵波总是跑得比横波快，如果按 vs_max 去定 dt，模拟会在纵波经过高速度区时直接数值爆炸。时间步数等于模拟时长除以 dt，模拟 1 s 的波场记录需要 1000 步，这个数量级适合在笔记本上调试。

有些实现里会建议用二阶时间精度的 CFL 上限 0.7 甚至 0.9，那是针对声波标量方程和低阶差分格式的结论，弹性波方程里横波与纵波耦合、自由表面边界条件又会让局部不稳定性提前出现，我自己的习惯是把 CFL 控制在 0.35～0.45，慢是慢一点，但不用半夜爬起来看程序有没有发散。

3.3 最小可运行代码骨架：时间循环里到底在算什么

用 Python 写一个纯 NumPy 实现，不考究性能，目的是把每个计算环节摊开看清楚。主循环里只做五件事：更新切应力对速度的贡献、更新正应力对速度的贡献、更新正应力、更新剪切应力、加载震源项。

nt = 1000 vx = np.zeros((nz, nx)) vz = np.zeros((nz, nx)) sxx = np.zeros((nz, nx)) szz = np.zeros((nz, nx)) sxz = np.zeros((nz, nx)) # 四阶交错网格差分系数，dx 方向 c1 = 9.0 / 8.0 / dx c2 = -1.0 / 24.0 / dx for it in range(nt): # 由应力更新速度分量 dvx_dx = c1 * (sxx[2:-2, 3:-1] - sxx[2:-2, 1:-3]) + \ c2 * (sxx[2:-2, 5:-1] - sxx[2:-2, 0:-4]) # 注意这里垂直方向用 dz 对应系数 vx[2:-2, 2:-2] += dt / rho[2:-2, 2:-2] * dvx_dx # 同理更新 vz，由 szz 和 sxz 的组合 # ... 此处省略 vz、sxx、szz、sxz 的同类更新 # 震源项以集中力形式加载到 vz 分量 vz[src_z, src_x] += dt * source_wavelet[it] / rho[src_z, src_x]

上面代码里最值得解释的是索引偏移。交错网格差分时，速度点和应力点不在同一套网格索引上，卷积模板要在前后各错开两个点（四阶精度用 4 个点），所以边界上各留两个网格点不做差分更新，这些点的值需要另外通过吸收边界或单向差分补充。很多初版实现跑出来的图边界上有一圈高振幅亮斑，就是没处理这层索引内外之差导致的。

实际跑这种规模的模型，纯 NumPy 双层循环慢得离谱，正确做法是把所有运算写成数组切片的广播形式，利用 NumPy 在 C 层面做循环。上面代码里刻意写了半截，因为完整实现需要把五个分量方程全部展开，建议自己动手补齐 vx 在 z 方向的应力贡献、vz 的 x/z 两个方向贡献以及三个应力分量的速度散度组合，补完这一步，你对弹性波正演的理解能超过大半只跑过现成程序的人。

4. 吸收边界与震源加载：让模拟结果真正可用的两个环节

4.1 完全匹配层（PML）：怎么加才能压住人为边界反射

模拟区域的四条边如果不做处理，波传到边界会被反射回来，叠加在有效波场上形成假同相轴。简单吸收带的思路是把边界区域内的波场振幅乘以衰减系数，但这会带来边界反射和频散；工程上长期使用且效果可靠的是完全匹配层。

PML 的核心思想是把波动方程在吸收层内的解替换为一种在数学上与内部解匹配、但在吸收层内指数衰减的变换形式。常见的实现方法是复坐标伸缩法：对吸收层内的空间导数算子做一个拉普拉斯域变换，引入阻尼剖面 σ(x)。实现时通常需要把速度和应力在 PML 区域内拆分成两个方向的分量，比如 vx 的 x 方向贡献和 z 方向贡献分别用不同的衰减系数更新，这就让代码量一下翻倍。

实际建模时长宽各加 20 个网格点的 PML 层，总厚度大约 200 m。衰减剖面从内到外按二次或三次函数从 0 增加到最大值，最大值通常取 σ_max = -3 * vp_max * log(0.001) / (2 * L_pml)，其中 L_pml 是 PML 层物理厚度。注意如果 σ_max 取大了，PML 层本身会变成一个引射反射界面；取小了又压不住低频反射。判断 PML 是否调好的标准很简单：看波场快照在边界附近有没有圆弧状弱能量——有则说明衰减还不够或剖面过渡不光滑。

一个更隐蔽的坑在介质参数进入 PML 时发生突变。如果 PML 区域的 vp、vs、ρ 直接沿用到内部边缘的赋值，而内部模型在边缘处本身就有横向不均匀性，PML 和内部之间的阻抗差会生成虚假反射。稳妥做法是把 PML 区域内的材料参数设成内部边缘值的平滑延拓。

4.2 震源子波与加载方式：爆炸源、集中力源的区别

震源类型决定了波场中 P 波和 S 波能量的比例关系，这是弹性波正演和声波正演一个明显的操作差异所在。常见两种源：

爆炸源只在正应力分量中加载，不加载剪切应力，理论上只辐射 P 波，在均匀各向同性介质中不应有 S 波。垂直集中力源同时扰动正应力和剪切应力，天然产生 P 波和 SV 波，更接近可控震源或微地震事件的辐射特性。

加载时建议用主频 25 Hz 的雷克子波：

def ricker(freq, dt, nt): t = np.arange(nt) * dt t0 = 1.0 / freq factor = (np.pi * freq)**2 return (1.0 - 2.0 * factor * (t - t0)**2) * \ np.exp(-factor * (t - t0)**2)

子波零延迟时间 t0 取一个主频周期（1/freq），确保 t=0 时振幅接近零，避免零时刻直接加载一个巨大冲量，否则波场中会混入直流分量，后面的走时分析全是错的。

加载位置与方式：震源网格点下标 (iz_src, ix_src) 通常放在自由表面以下 5～10 个网格点，避免自由表面反射和直达波在起跳点纠缠。把子波按时间步逐点加到对应分量的网格上，使用如下形式：

source_wavelet = ricker(25.0, dt, nt) for it in range(nt): # 在震源点添加垂直集中力源 vz[src_z, src_x] += dt * source_wavelet[it] / rho[src_z, src_x] # 如果采用爆炸源，则改为： # sxx[src_z, src_x] += dt * source_wavelet[it] # szz[src_z, src_x] += dt * source_wavelet[it]

注意分母上的密度。集中力源加载的是力密度（单位体积力），物理量纲上是加速度，所以加速度 = 力 / 密度；而爆炸源加载的是应力直接叠加，不除密度。这个细节搞反的后果是振幅尺度差一个密度量级，且 P/S 波相对振幅全乱。

还有一类源是纯剪应力源，用于模拟横波辐射，但实际数据处理中用得少。如果你做微地震正演，震源要换成矩张量加载，规则是给九个应力分量同时加不同权重的子波。初次练习不必碰这个，先用垂直集中力源跑通主流程最重要。

5. 弹性波正演模拟避坑：5 个高频翻车现场与排查思路

5.1 模拟爆炸：波场在几百步内出现高频振荡

现象：程序运行到 200～400 步时，整个波场快照布满细密的棋盘状纹路，振幅指数级增长，边界反射还没走到就变成噪声场了。

原因：时间步长超过 CFL 限值，或者局部区域 vp 异常大（比如模型里混入了 6000+ m/s 的高速体但 dt 按平均值算）。另一种可能是不小心把整网格差分系数用到了交错网格上，造成高频截断不满足。

解决：重新按 vmax = max(vp) 计算 dt，CFL 取 0.4 以下重跑。如果高频振荡只出现在震源附近并随着距离增长而扩散，检查空间差分系数 c1、c2 的正负号——四阶交错网格系数的符号顺序写错一个，波场会在一个波长的距离内彻底散掉。

5.2 近震源处有拖尾振荡，像梳子齿一样

现象：波场快照中震源周围有一圈周期性的细密波纹，不是正常环形波前。

原因：空间网格太粗，最短波长内只有两三个网格点，四阶差分无法分辨该波长成分，数值频散把能量在网格尺度上散开。

解决：加密网格到最短波长 6 个点以上，或降低震源主频。如果不想同时增加计算量，可以只把震源附近区域加密，但交界面处要插值衔接，麻烦，不如先整体加密看结果趋势再决定。另一种有效手段是把空间差分阶数从四阶提到八阶，但对应的 PML 层差分公式要同步改，不能只换系数。

5.3 模拟记录里背景有一条水平直线状强振幅

现象：整个模拟时间窗内，固定某一深度处的所有检波器都收到一个水平延伸的强能量条带，走时与震源距离无关。

原因：这是经典的「零频漂移」，常见于采用爆炸源但未让子波均值归零。雷克子波如果截断时间窗太短（比如只取主频半个周期），子波本身带有直流分量，加载到应力分量后形成静态偏移，波场快照上就是一条水平的直流亮线。检查子波均值是否为零，可加一步：

source_wavelet -= np.mean(source_wavelet)

我见过不少程序在子波上的这个减法比吸收边界调参更能改善记录质量，先做这个再碰 PML。

5.4 转换波看得到但极性反了

现象：合成三分量记录中，PS 波（转换横波）同相轴的极性符号与理论相反，或者与纵波记录的相对极性关系不对。

原因：应力与速度分量之间的符号约定不一致。速度-应力方程写成∂vx/∂t = (1/ρ) * ∂σxx/∂x还是= (1/ρ) * (∂σxx/∂x + ...)，以及剪切应力更新方程中∂vx/∂z与∂vz/∂x的相对正负号，不同教材有不同约定。换个符号体系，波场极性就整体翻转。解决方式是锚定一个参考解来校验：在均匀半空间模型中计算垂直分量上 P 波初动方向，P 波到达时质点运动方向应与传播方向同向（即向下传播时向下运动），按这个基准去校正符号，别相信代码注释里写的极性说明。

5.5 边界反射压不住，PML 失效

现象：PML 区域内部波场仍有弧形反射回传到内部区域，反射波振幅比理论值高一个量级，而且反射看起来来自 PML 层内某个位置而不是外部边界。

原因：最常见是 PML 参数剖面设计错了。衰减系数 σ 从内边界处直接跳到最大值，产生了阻抗突变面；或 PML 区域的空间差分阶数低于内部区域，导致入射波在 PML 入口处发生频散；还有一种是把 PML 内部的速度参数设置与内部不连续，形成人工反射界面。

解决：先检查 σ 从 0 平滑过渡到 σ_max 的形态，用二次或三次函数；再确认 PML 区域差分方式与内部一致；最后将 PML 厚度加大到 30 个网格点。调好后做一个空模型测试：均匀半空间，震源放在中心，看边界反射振幅是否低于主波振幅的 1%。

6. 用波场快照与合成记录验证模拟：走时、振幅与偏振的综合检验

跑通一个正演程序只是第一步，真正让模拟结果具备说服力的是它能不能通过三层验证：解析解对比、走时一致性、偏振特性。我的习惯是每换一种新介质模型，先做这三件事再谈后续应用。

第一层，解析解对比。在各向同性均匀半空间中，兰姆问题（Lamb's problem）的解析解给出了垂向集中力源产生的波场精确表达式，这是「金标准」。不要拿自己离散模型的数值解直接跟解析解逐点对比，因为网格频散和震源离散化总会带来细微差异。我一般取相距震源数十个波长的检波点，对比 P 波和 S 波到达时刻的均方根误差，误差小于一个时间步长即可认为程序正确。

第二层，走时检验。在两层模型中由射线理论手工计算 P 波和 PS 转换波的理论走时，与合成记录中拾取的走时对比。P 波反射走时公式 t_P = sqrt((2h)² + x²) / vp，PS 转换波走时要分段计算下行 P 波和上行 S 波路径，没有解析闭式解，可以用简单的二分法求 Snell 定律条件下的最小走时。偏差超过两个时间步长时，优先怀疑震源子波零时刻定义和检波器基准面高程。

第三层，偏振检验。弹性波的三分量合成记录里，P 波质点振动方向平行于传播方向，S 波垂直于传播方向。取直达 P 波的时窗，对水平分量和垂直分量作质点运动图（hodogram），P 波应该是一条过原点的直线，斜率等于入射角的正切；如果看到椭圆形状，说明三分量增益不匹配或者模拟输出坐标系定义错了。

这三层验证做完，才能放心把模拟结果用于 AVA 响应分析或者观测系统参数论证。过去我在一个 VSP 正演项目里吃过亏：当时跨过了解析解对比直接做观测系统设计，结果合成波形里 P 波初动极性没问题，但后续的转换波组参与野外数据的对比始终差半拍，最后排查发现是 Vz、Vx 分量的坐标轴方向定义与采集记录约定不一致，白做了一个月的观测方案。后来就把偏振检验作为每次正演任务结束前必须跑一遍的例行程序，再没出过这类错。

弹性波正演到这里骨架就齐了：方程、网格、CFL、PML、震源、验证，一层一层搭起来，每一步都有对应的检查办法。调参遇到玄学问题别急着改代码，先回去看子波均值和 CFL 数，这两个地方占八成翻车原因；剩下的二成用参考解对比定位，比反复试参数更接近问题本质。希望这篇笔记能帮你少走几段弯路。

本文还有配套的精品资源，点击获取