简介：这份资源面向地震学与地球物理数值模拟方向的学习者和研究者，提供基于弹性波方程的地震波场正演模拟MATLAB实现。弹性波方程描述地震波在地壳中的传播行为，正演模拟通过数值方法预测波在不同地质结构中的传播路径、振幅变化与到达时间，而该脚本采用10阶精度差分算法离散方程，能更细腻地刻画高频波场细节，提升模拟精度。压缩包为rar格式，内含1个m文件，整体约2KB，属于轻量级脚本资源，可直接在MATLAB环境中运行与二次修改。目前已有133人学习下载，适合作为高阶有限差分数值模拟的入门参考或算法验证素材。读者可从中获取高阶精度离散思路、波场正演流程框架，并据此进一步引入完美匹配层等边界条件，优化计算效率与物理合理性。

1. 弹性波方程正演模拟：一个 10 阶精度 MATLAB 脚本能跑出什么

如果你手头只有一份jiaocuo2_10jie.m，没有说明书、没有示例数据、没有依赖清单，第一反应大概率是「这玩意儿能直接跑吗」。我拿到这个包的时候也是同样的疑问。它做的事情其实很明确：用 MATLAB 对弹性波方程做正演模拟，空间差分精度拉到 10 阶，输出地震波场在时间上的演化。弹性波方程描述的是应力与质点速度之间的耦合关系，比声波方程多出横波分量，因此波场里能同时看到 P 波和 S 波，这对岩性识别和裂缝检测场景更贴近实际。10 阶差分意味着在同样的网格密度下，高频成分的数值频散被压得更低，波形更干净。代价是每步计算要访问更多邻点，内存和耗时都会上去。这份资源适合做地震波传播算法验证、教学演示，或者作为更复杂逆时偏移流程的正演引擎原型。不适合直接拿去做三维大规模生产，因为脚本本身没有边界吸收，能量会在网格边缘反弹回来。

2. 弹性波方程的一阶速度-应力形式与 10 阶差分实现

2.1 为什么用一阶方程组而不是二阶位移方程

弹性波方程常见写法有两种：二阶位移形式直接解质点位移，物理直观但边界处理麻烦；一阶速度-应力形式把速度分量和应力分量作为独立变量，用一阶偏微分方程组联立求解。后者在数值实现上更受欢迎，原因是它天然适合交错网格（staggered grid），速度和应力定义在不同半网格点上，差分模板更紧凑，自由表面和内部界面的处理也更自然。jiaocuo2_10jie.m走的就是这条路。脚本里你会看到vx、vz、sxx、szz、sxz这几组变量，分别对应水平速度、垂直速度、正应力和剪应力。它们按时间步交替更新，先算应力再算速度，或者反过来，取决于脚本采用的蛙跳格式。

交错网格的另一个好处是压制棋盘格振荡。如果速度和应力定义在同一套网格点上，中心差分会出现奇偶失耦，波场里冒出高频噪声。交错半个网格后，这种振荡被自然消除。代价是初始化震源和接收器时需要做空间插值，脚本里通常用简单平均或者直接取最近点，精度要求高的话得用高阶插值算子。

2.2 10 阶差分系数怎么来的

有限差分系数来自泰勒展开。对于一阶导数，2 阶中心差分只用左右各一个点，误差项是二阶；4 阶用左右各两个点，误差降到四阶；10 阶需要左右各五个点，共十个邻点参与计算。系数通过解一个线性方程组得到，方程组由泰勒展开的矩条件构成。MATLAB 里可以用A\b直接解，也可以查表硬编码。jiaocuo2_10jie.m大概率是硬编码系数，因为这样运行时更快，也避免每次启动都解方程。

10 阶差分在波数域的有效带宽比 2 阶宽得多。举个例子，同样每个波长采样 5 个点，2 阶差分已经严重频散，波形后面拖尾巴；10 阶差分还能保持波形基本不失真。这就是高阶精度的价值。但要注意，阶数越高，边界处理越棘手。内部点用 10 阶模板没问题，靠近网格边缘的点没有足够邻点，必须降阶或者用单侧差分。脚本如果没有做降阶处理，边缘附近会出现异常反射，这也是后面避坑章节要展开的。

2.3 脚本核心循环的代码骨架

下面这段代码是根据jiaocuo2_10jie.m的典型结构还原的骨架，不是原文件逐行复制，但变量命名和更新顺序与常见实现一致。你可以对照自己的脚本看结构是否吻合。

% 网格与模型参数 nx = 400; nz = 400; % 网格点数 dx = 10; dz = 10; % 空间步长(m) dt = 0.001; nt = 1500; % 时间步长(s)与总步数 npml = 0; % 本脚本未加PML，边界为自由表面或刚性边界 % 差分系数(10阶一阶导数中心差分) c = [1.211242675782, -0.089721679688, 0.013842773438, ... -0.001765659878, 0.000118679471]; c = [fliplr(c), 0, c]; % 对称扩展为10点模板 % 波场变量初始化 vx = zeros(nz, nx); vz = zeros(nz, nx); sxx = zeros(nz, nx); szz = zeros(nz, nx); sxz = zeros(nz, nx); % 震源 isrc = round(nx/2); jsrc = round(nz/2); f0 = 25; % 主频(Hz) t = (0:nt-1) * dt; src = (1 - 2*(pi*f0*(t-1/f0)).^2) .* exp(-(pi*f0*(t-1/f0)).^2); % 时间循环 for it = 1:nt % 应力更新：使用速度的空间导数 for i = 6:nz-5 for j = 6:nx-5 dvx_dx = sum(c .* vx(i, j-5:j+5)); dvz_dz = sum(c .* vz(i-5:i+5, j)); % 本构关系更新 sxx, szz, sxz % ...(略去具体弹性参数乘法) end end % 速度更新：使用应力的空间导数 for i = 6:nz-5 for j = 6:nx-5 dsxx_dx = sum(c .* sxx(i, j-5:j+5)); dsxz_dz = sum(c .* sxz(i-5:i+5, j)); % 运动方程更新 vx, vz % ... end end % 加震源 sxx(jsrc, isrc) = sxx(jsrc, isrc) + src(it); % 边界处理：本脚本未做吸收，边缘直接截断 end

这段骨架里最关键的是c数组，它就是 10 阶一阶导数中心差分系数。sum(c .* vx(i, j-5:j+5))这一行完成了沿 x 方向的 10 点加权求和。注意循环范围从 6 到 nz-5，因为边缘 5 个点没有完整模板，脚本直接跳过了。跳过意味着边缘不更新，波传过去会被硬截断，产生强烈反射。这是无边界脚本的典型行为。

参数方面，dx和dz决定空间分辨率，一般要求每最短波长至少采样 5 到 10 个点。dt受 CFL 条件约束，10 阶差分下 CFL 上限比 2 阶略小，通常取dt <= 0.5 * dx / vmax比较安全。f0是震源主频，越高频越考验差分精度，10 阶能撑到每波长 4 个点左右还不明显频散。nt要保证波场传到感兴趣区域之外，否则截断效应会污染记录。

2.4 弹性参数与稳定性条件

弹性波正演需要纵波速度vp、横波速度vs和密度rho。脚本里通常先定义vp和vs模型，再换算拉梅常数lambda = rho*(vp^2 - 2*vs^2)和mu = rho*vs^2。本构关系更新时用lambda和mu乘应力散度项。如果脚本只给了vp和rho，那它可能退化成了声波近似，横波信息就丢了。拿到脚本后先搜vs或mu，确认是不是真弹性。

稳定性方面，10 阶差分的 CFL 条件比 2 阶更紧。经验公式是dt * vmax * sqrt(1/dx^2 + 1/dz^2) <= C，其中C随阶数升高而减小。2 阶时C约 1，10 阶时可能降到 0.6 左右。如果跑几步就出现数值爆炸，先把dt减半试试。另一个稳定性隐患是弹性参数突变界面，lambda和mu在界面两侧差异大时，交错网格需要做参数平均，否则会引入虚假绕射。常见做法是对mu取调和平均，对lambda取算术平均，脚本里如果直接取节点值，界面精度会掉。

3. 在 MATLAB 里跑通脚本：从参数检查到波场快照

3.1 运行前的依赖与版本确认

jiaocuo2_10jie.m是纯 MATLAB 脚本，不依赖额外工具箱的概率较大，但保险起见还是检查一下有没有用到signal、image这类工具箱函数。在命令行输入which jiaocuo2_10jie确认脚本在路径里，然后edit jiaocuo2_10jie打开看一眼开头有没有addpath或者load语句。如果有load，说明它依赖外部模型文件，压缩包里应该还有.mat或.dat，别只把.m拷出来就跑。

MATLAB 版本方面，脚本如果用了隐式扩展或者string类型，需要 R2016b 以上。老版本可能报Undefined function错误。我一般先在 R2020a 上试，没问题再往下兼容。如果脚本里全是for循环，那版本兼容性很好，R2010 都能跑，只是慢。

3.2 关键参数的手动调整

脚本开头通常有一块参数区，你需要根据自己机器和模型改几个值。下面这张表列出了最常改的参数和推荐范围。

改dt的时候注意 CFL。一个快速估算：dt_max = 0.5 * min(dx,dz) / max(vp(:))。如果vp最大 4000 m/s，dx=10，那dt_max约 0.00125 s。脚本里如果设了 0.002，大概率会炸。先按这个公式算一遍再跑。

3.3 运行与波场可视化

直接run jiaocuo2_10jie或者在编辑器里按 F5。如果脚本没有内置绘图，跑完只会得到一堆变量在 workspace 里。你可以在最后加一段快照代码：

% 在时间循环内每隔200步存一张快照 if mod(it, 200) == 0 figure(1); imagesc(vx); colorbar; title(['vx at step ', num2str(it)]); axis image; drawnow; end

imagesc把vx矩阵按颜色映射画出来，axis image保持纵横比，drawnow强制刷新。这样你能实时看到波前从震源向外扩散。P 波和 S 波因为速度不同会逐渐分离，P 波在前，S 波在后，波前面形状取决于模型是否各向同性。如果看到波场里出现网格状高频噪声，多半是差分系数用错或者 CFL 超了。

跑完后可以用save('wavefield.mat', 'vx', 'vz', 'sxx', 'szz', 'sxz')把波场存下来，方便后续做偏移或者分析。如果只想看某个接收点的波形，在循环里加rec(it) = vx(jsrc, isrc+50)这样的语句，把时间序列记下来再plot。

3.4 性能瓶颈与向量化改造

纯for循环的 MATLAB 脚本跑 400×400 网格、1500 步，可能要几分钟到十几分钟。瓶颈在双重循环里的sum(c .* ...)，每次都要分配临时数组。向量化改造的思路是把内层j循环去掉，用矩阵切片一次性算完所有点。比如沿 x 方向的差分可以写成：

% 向量化沿x方向的10阶差分 dvx_dx = zeros(nz, nx); for k = 1:10 dvx_dx(:, 6:nx-5) = dvx_dx(:, 6:nx-5) + ... c(k) * (vx(:, 6:nx-5+k-5) - vx(:, 6:nx-5-k+5)); end

这样只保留外层k循环，内层全部矩阵化。MATLAB 对矩阵运算做了优化，速度能提升几倍。再进一步可以用gpuArray把波场搬到 GPU，但需要 Parallel Computing Toolbox，而且显存要够。400×400 单精度大概占几 MB，普通显卡都能吃下，加速比可能到 10 倍以上。不过脚本原本没做 GPU 适配，你得自己改数据类型和函数调用。

4. 无边界正演模拟的避坑与排查

4.1 波场在边缘反弹形成假反射

现象：波传到网格边缘后，过一段时间从边缘往回传，和真实波场混在一起，记录上出现无法解释的同相轴。原因：脚本没有加吸收边界，边缘点被直接截断，相当于刚性边界，反射系数接近 1。解决：在网格四周加 PML 层，层数一般 10 到 20 个点，在 PML 内引入阻尼因子逐步衰减波场。如果不想改脚本结构，至少把模型区域放大，让感兴趣的区域远离边缘，等反射回来时已经超出记录时间。常见做法是模型四周各留 50 个点的缓冲。

4.2 高阶差分边缘降阶导致的异常

现象：边缘附近出现高频振荡或者波形畸变，和内部点结果对不上。原因：10 阶模板需要左右各 5 个邻点，边缘 5 个点没有完整模板，脚本直接跳过不更新，相当于那里波场被冻结。解决：对边缘点做降阶处理，比如最外面 5 层用 2 阶差分，往里逐层升阶。或者用单侧差分模板，但系数要重新推导。偷懒的办法是把有效波场区域缩小，边缘 5 个点不参与成像。

4.3 CFL 超限导致数值爆炸

现象：跑几十步后波场值变成NaN或Inf，图像全白或全黑。原因：dt太大，不满足 CFL 条件，蛙跳格式发散。解决：按dt <= 0.5 * min(dx,dz) / max(vp(:))重算dt，先减半跑通再逐步调大。另外检查vp模型里有没有异常大值，比如 10000 m/s 这种，会把 CFL 上限拉得很低。

4.4 震源加载方式不对导致波场不对称

现象：波场不是从震源向外均匀扩散，而是偏向一侧或者出现方向性。原因：震源只加在sxx上，没有同时加在szz上，或者加的位置在交错网格上偏移了半个网格。解决：爆炸震源应该同时给sxx和szz加相同的时间函数，且加在应力节点上。如果加在速度节点上，需要做空间插值。检查脚本里震源是加在sxx(jsrc,isrc)还是vx(jsrc,isrc)，前者是应力源，后者是速度源，波形不一样。

4.5 弹性参数界面处理粗糙引入虚假绕射

现象：在两层模型界面附近出现不该有的绕射波，能量不强但污染记录。原因：lambda和mu在界面两侧直接取节点值，交错网格上半网格处的参数没有做平均。解决：对mu取调和平均，lambda取算术平均，或者用更精确的等效参数方法。如果脚本里参数是直接赋值给整个矩阵，没有平均步骤，界面精度会掉一个档次。

5. 从无边界到 PML：给 10 阶正演加吸收边界的实操

无边界脚本跑通之后，下一步自然是加 PML。PML 的核心思想是在边界区域引入一个阻尼项，让波场在传播过程中按指数衰减，到达网格边缘时已经衰减到接近零，反射自然就没了。实现上分两步：定义阻尼剖面，然后在波场更新方程里乘上衰减因子。

阻尼剖面常用多项式形式：d(x) = d0 * (x/L)^n，其中L是 PML 层厚度，n通常取 2 到 4，d0和最大速度、主频有关，经验公式是d0 = -3 * vmax * log(0.001) / (2 * L)。在 MATLAB 里可以这样生成一维剖面再扩展到二维：

npml = 20; % PML层数 L = npml * dx; % PML物理厚度 vmax = max(vp(:)); d0 = -3 * vmax * log(0.001) / (2 * L); n = 3; % 多项式阶数 % 一维阻尼剖面，从内到外递增 x_pml = (1:npml) * dx; d_x = d0 * (x_pml / L).^n; % 扩展到二维：左右和上下各一层 d_left = zeros(nz, nx); d_left(:, 1:npml) = repmat(fliplr(d_x), nz, 1); d_right = zeros(nz, nx); d_right(:, end-npml+1:end) = repmat(d_x, nz, 1); d_top = zeros(nz, nx); d_top(1:npml, :) = repmat(fliplr(d_x)', 1, nx); d_bottom = zeros(nz, nx); d_bottom(end-npml+1:end, :) = repmat(d_x', 1, nx); % 总阻尼，角落区域叠加 d_total = d_left + d_right + d_top + d_bottom;

这段代码生成一个和波场同尺寸的阻尼矩阵d_total，内部区域为零，PML 区域逐渐增大。然后在时间循环里，每次更新完vx、vz、sxx、szz、sxz之后，乘上衰减因子：

% 在时间循环内，波场更新之后 vx = vx .* exp(-d_total * dt); vz = vz .* exp(-d_total * dt); sxx = sxx .* exp(-d_total * dt); szz = szz .* exp(-d_total * dt); sxz = sxz .* exp(-d_total * dt);

exp(-d_total * dt)是衰减因子，d_total越大衰减越快。PML 层内波场被逐步压低，到达边缘时已经接近零，反射自然消失。注意d_total在角落区域是左右和上下阻尼之和，衰减更强，这是对的，因为角落波场来自两个方向。

加完 PML 后重新跑一遍，对比边缘反射是否消失。一个简单的验证方法：把震源放在中心，记录边缘点的波形，加 PML 前会看到明显的二次到达，加 PML 后二次到达应该被压到很小。如果还有残留，增大npml或者调大d0。但d0太大会导致 PML 内部本身产生反射，因为阻尼变化太剧烈。一般npml取 10 到 20，n取 2 到 4，d0按公式算，微调范围不超过 20%。

还有一个细节：PML 只对特定方向的波有效，对于掠射波吸收效果会变差。如果模型里波几乎平行于边界传播，PML 可能需要加厚或者改用卷积 PML。不过对于教学和中小规模模拟，标准 PML 已经够用。

从那以后我每次拿到无边界正演脚本，都先跑一遍看边缘反射，再决定加 PML 还是直接放大模型。这个习惯帮我省了很多事后排查的时间。希望帮到你。

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