企业官网建设流程全解析

简介：本资源是一套面向地球物理勘探与计算地球科学方向的MATLAB正演模拟工具，专为电子信息工程、数学及计算机专业本科生课程设计、毕业设计等实践环节开发，解决VTI介质中弹性波传播的高精度数值模拟问题。压缩包仅含2个文件（1个核心MATLAB脚本.m + 1张运行结果图.png），总大小16KB，轻量易用，适配MATLAB 2014a至2021a多个版本。已有154人学习下载，说明其在教学实践与算法验证场景中具备较强实用性。用户可直接运行主程序，获得含PML吸收边界条件的高阶交错网格有限差分正演结果；代码采用参数化设计，介质参数、网格步长、时间步进等关键变量均集中定义并附详细中文注释，逻辑清晰、易于修改与拓展，特别适合初学者理解弹性波方程离散化原理及边界处理技术。 前两天从网盘里拉下来一个压缩包，名字叫“基于matlab模拟有限差分VTI介质弹性波方程的高阶交错网格正演+pml吸收边界条件.zip”。我本来觉得这类代码包十个里有八个是跑不起来的课程作业，结果解压之后花了点时间梳理，发现这套程序把地震波正演里的几块硬骨头——VTI各向异性介质、高阶交错网格有限差分、PML吸收边界——串得很完整。这篇文章就把这套程序的原理、实现和调参细节拆开聊聊，给正在做地震波正演、各向异性介质模拟或者打算用MATLAB做数值模拟的同学一个参考。

作为地震正演里最常用的数值手段之一，有限差分法在波动方程求解中的地位不用多讲。但真正上手写一套能用的程序，和上课考试完全是两码事。很多人卡在几个地方：为什么介质要分各向同性和VTI？交错网格到底好在哪？高阶系数怎么算？PML层为什么总在地边界反射？这几点搞不明白，程序跑出来的结果就没法信。这篇文章会把这些都铺开说清楚。

1. 这个项目到底在解决什么问题

1.1 VTI介质是什么，为什么值得专门做正演

VTI（Vertical Transversely Isotropic，垂直横向各向同性）介质，是地震勘探里最常见的各向异性介质模型。你可以把它想象成一本叠起来的书：水平方向书页排列都一样，垂直方向则有明显的层理结构。页岩、薄互层、裂缝发育地层，在地震波传播上都会表现出这种“水平方向性质相同、垂直方向性质不同”的特征。

在VTI介质里，弹性波传播速度和方向有关。纵波速度不再是单独一个数值，而是由垂直方向和水平方向两个基准速度，加上各向异性参数共同决定。如果做正演的时候还拿各向同性公式硬套，炮集记录里的走时会出现系统性偏差，偏移成像位置会错位，这是实际生产中不能接受的误差。

这套程序里用的参数体系，通常不是直接拿刚度矩阵C11、C13、C33、C44硬算，而是用Thomsen提出的各向异性参数。Vp0是垂直方向纵波速度，Vs0是垂直方向横波速度，ε描述纵波各向异性强度，δ控制近垂直方向的速度变化率。理解这几个参数是读懂整个程序的第一步。

1.2 为什么选弹性波而不是声波

很多初学正演的人第一反应是：先做声波方程，简单，跑得快。但声波方程对VTI介质是先天不足的。VTI介质中P波和SV波是耦合的，两者在边界和界面处互相转化，声波方程根本描述不了这种转换现象。如果目标是研究转换波、横波分裂或者各向异性参数反演，就必须上弹性波方程。

代价也很直观：弹性波方程涉及的速度场分量（vx、vz）和应力场分量（τxx、τxz、τzz）加起来五个波场变量，内存和计算量大概比声波方程高一个量级。这个程序的做法是用速度-应力一阶方程组，也就是把弹性波方程的二阶位移形式拆成速度和应力的一阶偏微分方程组，这样方便用交错网格进行时间和空间的同时推进。

二维VTI介质中，平面内只有P-SV波系统，弹性刚度矩阵里C66其实用不到。三维情况下SH波才会牵出C66。所以你看程序里刚度系数更新只操作C11、C13、C33、C44四个量，就是这个道理。这个细节很多人会漏，以为二维程序也要给六个刚度系数。

2. 交错网格与高阶有限差分的核心原理

2.1 交错网格的空间配置

交错网格（Staggered Grid）的思路，说穿了就是“你有的我正好没有，我有的你正好没有”。速度分量和应力分量在空间上错开半个网格点放置。以二维模型为例，vx定义在格点(nz, nx)的右边界，vz定义在上边界，法向应力τxx、τzz定义在格点中心，切向应力τxz定义在格点的左上角。

这样安排的最大好处是：每个空间导数都只需要做一次半个网格点的差分，不需要插值。比如∂τxx/∂x，τxx在格点中心，而vx正好在两个τxx的中点，所以这个导数天然就是(vx右边减去vx左边)。差分模板的误差是O(h²)起步，如果用高阶差分系数，还可以进一步压低误差。

对比一下普通网格：如果所有变量都在同一位置，要计算一个变量在另一个变量位置上的值，就不得不做插值。插值本身会引入额外的数值误差，还会扩大差分模板的有效长度，在边界处更麻烦。所以地震波正演现在几乎清一色用交错网格，不是没有原因的。

在时间推进上，交错网格同样采用半时间步：速度场在n+1/2时刻更新，应力场在n+1时刻更新，两者互相利用对方在“中间时刻”的旧值。这样只需要存储相邻时刻的波场，不需要像龙格库塔那样存多级中间量，内存压力小很多。

2.2 高阶差分系数怎么算

普通二阶差分，就是对导数用一个三点模板近似。高阶有限差分则用更多点，让每个点贡献一个加权系数，组合起来使得截断误差达到更小的阶数。程序里最常见的配置是时间二阶、空间十阶，或者时间二阶、空间十二阶。空间阶数提高之后，同样的网格剖分下频散更小，模拟出来的波前面更加光滑。

交错网格的高阶差分系数，形式上是用泰勒展开求出来的。以x方向一阶导数为例，假设用2M个半网格点的加权和来逼近中心点导数，系数a_m满足：

Σ a_m * (m - 0.5) = 1

Σ a_m * (m - 0.5)^(2j-1) = 0，j = 2, 3, ..., M

由于交错网格的采样点只在半整数位置，偶数阶矩自动为零，所以只需要解一个M×M的线性方程组。这个方程组完全可以在MATLAB里用一行代码解出来：

function a = fd_coeff(M) A = zeros(M, M); b = zeros(M, 1); for j = 1:M b(j) = (j == 1); for k = 1:M A(j, k) = (k - 0.5)^(2*j - 1); end end a = A \ b; end

调用一下：fd_coeff(2)得到四阶系数[1.125, -0.0417]，fd_coeff(5)得到十阶系数。把这个函数放到程序里，就不用查表了，想换几阶精度随时改。系数求出来之后，空间导数算子就是把这些系数应用到更新公式里，比如速度对x的偏导：

∂vx/∂x ≈ (1/dh) * Σ_{m=1}^{M} a_m * (vx_at_x_plus - vx_at_x_minus)

注意半网格点的索引关系，写矩阵切片的时候特别容易错位。我在调试这套程序的时候，至少有两次跑出来波场像长了毛的刺猬，最后检查都是索引方向写反了。

2.3 时间推进与稳定性条件

时间上一般用二阶精度，也就是常规的前向/后向差分（蛙跳格式）。每个时间步里，先由旧应力场更新速度场，再由新速度场更新应力场。虽然时间只有二阶，但因为时间步长dt通常取得比空间步长小一个量级，实际数值误差主要来自空间差分。

稳定性条件是正演模拟绕不开的硬约束。交错网格下二维弹性波方程的CFL条件，近似写作：

v_max * dt / dh ≤ 1/√(2)

这里的v_max是模型中的最大波速，dh是空间步长。对于各向异性介质，v_max至少要取最大相速度方向的准纵波速度，比垂直方向的Vp0可能高出不少。

实际操作里我习惯留足余量，把CFL数控制在0.5左右。比如模型最大速度v_max=3000 m/s，空间步长dh=5 m，那么：

dt ≤ 0.5 * dh / v_max = 0.5 * 5 / 3000 ≈ 0.000833 s

取dt=0.5 ms，保证稳定。程序默认参数如果跑发散，第一件事就是检查这个比值，而不是去调整什么边界条件。

3. PML吸收边界：从公式到代码

3.1 拉伸坐标系的思想

人工边界反射是有限差分正演里最让人头痛的问题。模型四周截断了，波传到边界如果不做处理，会原路弹回来，把真实波场盖得什么都看不清。PML（Perfectly Matched Layer，完美匹配层）的原理，是在模型外部加一圈特殊介质，让波进入这个区域之后被指数衰减，并且在PML与内部介质的分界面上尽量不产生反射。

PML在数学上的解释有几种，最具操作性的是复坐标拉伸。简单说，把偏微分方程中的空间坐标x换成它的复函数：

x̃ = x - (i/ω) ∫ d(x) dx

i是虚数单位，ω是角频率，d(x)是随位置变化的衰减系数。这样做的效果是：进入PML区域的波场在频域上自然带上一个指数衰减因子，波还没到达截断边界，能量已经基本耗尽。

在时域实现里，这个拉伸会引入辅助变量。程序里最常见的版本是CPML（卷积PML），它对不同方向的导数分别做记忆变量卷积更新。比如x方向的导数项：

∂/∂x 的空间差分结果，叠加一个记忆变量ψ

ψ的更新公式为：

ψ^{n+1} = b_x * ψ^n + a_x * (∂/∂x的差分结果)

b_x = exp(-(d_x + α_x) * dt)

a_x = d_x * (b_x - 1) / (d_x + α_x)

其中d_x是衰减函数，α_x是一个很小的正数，用来保证极低频情况下的稳定性。每个方向各自维护一组记忆变量，最后应力更新的导数项等于原始差分加上记忆变量，PML就实现了。

3.2 衰减函数怎么选

衰减函数d(x)通常不是常数，而是在PML层内从分界面处的0逐渐增大到外边界处的最大值d_max。程序里最常用的是二次渐变和三次渐变，二次渐变已经够用，三次渐变更稳一点。

d(x) = d_max * (x / L_pml)^2

这里x是当前点到PML内边界的距离，L_pml是PML层的物理厚度。d_max的经验公式：

d_max = - (M + 1) * v_max * ln(R) / (2 * L_pml)

M是渐变幂次，R是理论反射系数，通常取0.001。举个例子：v_max=3000 m/s，L_pml=100 m（20个5 m网格），M=2，R=0.001：

d_max = 3 * 3000 * ln(1000) / 200 ≈ 3 * 3000 * 6.9 / 200 ≈ 310.5

这个值大概在几百量级是正常的。如果d_max给得太大，PML区和内部区域的分界面会有虚假反射；给得太小，边界吸收不干净，波会从截断边界反弹回来。

3.3 PML更新实现要点

实现PML的时候，最容易踩的坑是角点处两个方向的衰减同时生效。在PML的角部区域，x方向和z方向都要施加阻尼，所以那里需要同时维护两组记忆变量，不能只做单方向处理。有些简化程序为了省事把角点直接置零，波到了角落还是会被反射，效果很差。

程序里通常配置的PML厚度是10到30个网格点。厚度太小，吸收效果不理想；厚度太大，计算量白白增加。我一般用20个点左右，乘以空间步长就是实际的吸收层宽度。检查PML是否有效，最简单的方法是跑一个均匀介质模型，观察波场快照的边界区域有没有明显的二次弧状反射波；或者把边界处的道集振幅调出来，看它们是不是比内部道逐级衰减到很小。

PML衰减函数要从零开始渐变，千万不要在PML入口处直接给一个大的d值。这种突变等同于人为设置了一个反射界面，边界反射比不设置PML还严重。

4. 完整正演流程拆解

4.1 模型参数与网格设计

一套典型参数的二维VTI模型可以这样定：模型大小400×400网格点，dh=5 m，dt=0.5 ms，总时间2 s，中心频率f0=25 Hz。介质参数用Thomsen参数表示时，先给一组基准值：Vp0=2500 m/s，Vs0=1400 m/s，ε=0.2，δ=0.05，ρ=2200 kg/m³。

这组参数需要转换成波场更新用的刚度系数。在MATLAB里可以这样算：

Vp0 = 2500; Vs0 = 1400; epsilon = 0.2; delta = 0.05; rho = 2200; C33 = rho * Vp0^2; C44 = rho * Vs0^2; C11 = C33 * (1 + 2 * epsilon); C13 = sqrt((C33 - C44)^2 + 2 * delta * C33 * (C33 - C44)) - C44;

算出来的值大致在以下量级：C33约1.375×10^10 Pa，C44约4.312×10^9 Pa，C11约1.925×10^10 Pa，C13约5.79×10^9 Pa。给刚度系数赋值的顺序直接影响程序的运行效率，MATLAB里把这些系数先构建成矩阵（每个网格点一份），再参与向量化运算，比在循环里反复计算要快得多。

网格参数设计要同时考虑稳定性和频散。最大波速是各向异性介质中沿水平方向传播的准P波速度，取Vp_horiz ≈ Vp0 * sqrt(1 + 2ε) = 2500 * sqrt(1.4) ≈ 2958 m/s。CFL条件检查下来，0.5 ms的dt是安全的。

4.2 震源设置

震源加载方式直接影响波场的激励类型。程序里一般用雷克子波（Ricker wavelet）作为时间函数：

s(t) = (1 - 2π² f0² (t - t0)²) * exp(-π² f0² (t - t0)²)

t0取1.2/f0左右，保证子波起始段接近零。中心频率f0=25 Hz时，t0=0.048 s，对应时间步在约96步时开始发力。

在交错网格里，震源的加载位置也要注意分量对应关系。如果是压力型震源，就加在应力分量τxx、τzz上；如果是力源，就加在速度分量vx或vz上；如果做的是爆炸源模拟，一般是在σxx和σzz上同时加一个负的震源时间函数。程序里常见的是加载在vz分量上模拟垂直力源，或者加载在τxx、τzz上模拟爆炸源，具体看你要对比什么资料。

加载格式上要注意：震源项加到波场更新公式的右端时，单位要换算对。速度-应力方程中，力的单位是加速度，所以震源项需要除以密度再乘以dt。如果直接拿子波值加进去，振幅比例全乱了，后期做波场快照会看不出有效信息。模拟中我用到的做法是：

src = (1 - 2*pi^2*f0^2*(t - t0)^2) * exp(-pi^2*f0^2*(t - t0)^2); vz(isrc, jsrc) = vz(isrc, jsrc) + src * dt / rho(isrc, jsrc);

这样源项的量纲才正确。

4.3 主循环实现

主循环是程序的灵魂，结构上是“先速度后应力，先内部网格后边界PML”。更新每个变量时，先用内部区域的交错差分公式在二维矩阵切片上做向量化运算，再单独处理PML区域。下面是一个简化的速度场vx更新代码片段（省略PML部分）：

% 更新 vx % Dx_tauxx: x方向 tauxx 导数 Dx_tauxx = (1/dh) * ( sum_a * (tauxx(nz, 2:nx+1) - tauxx(nz, 1:nx)) ); % 这里 sum_a 表示用高阶系数做加权叠加，需要展开成多阶循环或矩阵乘法 % Dz_tauxz: z方向 tauxz 导数 vx(2:nz+1, 2:nx) = vx(2:nz+1, 2:nx) ... + dt/rho(2:nz+1, 2:nx) .* (Dx_tauxx + Dz_tauxz);

实际写的时候，高阶差分系数展开成M个相邻点的加权求和，代码会比较长。一个实用的做法是写一个通用函数diff_staggered(field, axis, coeff)，内部用circshift或者矩阵切片完成半网格点差分，但需要注意边界范围。用circshift有个麻烦，它会循环填充边界，必须在循环之后把PML区域正确覆盖掉。

应力更新公式与之对称。比如τxx的更新需要vx对x的导数和vz对z的导数：

∂τxx/∂t = C11 ∂vx/∂x + C13 ∂vz/∂z

注意C11和C13在VTI介质里是不同的刚度值，各向同性情况下C11=C33=λ+2μ，C13=λ，程序就退化到各向同性了。这也是验证程序正确性的一个重要手段。

PML区域的更新建议单独用一个循环块处理。因为PML内部的记忆变量更新和主区域的差分格式不同，如果把两者混在一个向量化表达式里，代码可读性差，也容易出错。我的做法是主区内层用向量化，PML边界用循环或者部分向量化，两者分得很清楚。

4.4 波场快照与炮集输出

正演的结果一般有两种输出方式：一是动态波场快照，即固定时刻观察整个模型区域的波场分布；二是固定检波器位置记录时间序列，形成单炮记录。程序里这两者都要做。

动态波场快照在MATLAB里很简单，每若干时间步用imagesc画一帧，稍作停顿再更新。这样你就能亲眼看到P波和SV波在界面的转换、波前在PML区域的衰减过程，非常直观。我习惯每20个时间步存一帧，最后合成GIF或者AVI，方便反复观察。

单炮记录则是把每个时间步在预设检波器位置上的vz或者vx分量抽出来，存成一个nt × nrec的矩阵。对于二维VTI模型，检波器一般布设在模型表层，间距与网格间距一致。输出前要注意把有效分量选择好——如果震源是垂直力源，vz记录里的P波和SV波到时信息比较清晰，vx则可以观察横波偏振特征。

5. 常见问题与调参经验

5.1 波场发散、频散这些老问题

正演跑着跑着波场突然变成满屏噪声，通常是CFL条件被打破，dt超出稳定上限。此时把dt调小一半，再观察是否恢复稳定。另一个可能原因是刚度系数赋值错误，比如C13算出来是负数，或者各向异性参数过大，导致相速度出现非物理数值。检查方法是在均匀介质中跑一次小的模型，看波前是否光滑，波形低频是否正常。

网格频散表现为波场尾部拖着一条条振铃状尾巴，看起来像“梳子齿”。解决思路是提高空间差分阶数，或者加密空间网格。如果程序本来就是十阶差分，还是明显频散，那大概率是空间步长太大——一个波长至少要覆盖5到10个网格点。观测最小波长λ_min = v_min / f_max，f_max按震源主频的两三倍估算，然后反推dh。

5.2 PML参数没调好的表现

PML调参问题比稳定性问题更隐蔽。边界反射的典型表现是：在波场快照的四周能看到与主波前同心但晚到的圆弧状波场，振幅逐渐向内扩散。如果出现这种情况，优先检查PML层厚度和d_max。PML厚度太小（比如只有5个点）时，即使d_max再大，也不容易压住低频段的反射。

另一个常见问题是PML区和内部区域分界处参数不连续，在分界线上产生一条明显的亮线。这是d值一开始就给得太大造成的。把渐变函数改为从0开始的三次幂函数，情况会好很多。

5.3 MATLAB性能优化

MATLAB的循环是出了名的慢，主循环如果写成双层for循环，400×400网格、4000个时间步可能要跑到天荒地老。实测经验是：先把模型网格double类型改成single类型，内存占用减半，计算速度也有提升。然后把内部区域的波场更新全都改成矩阵切片运算，MATLAB对这类向量化操作优化得非常好。

如果要做多炮正演，可以在外部对每一炮使用parfor并行，此时每炮内部的循环仍用普通for或者向量化。但要注意parfor下随机数生成和文件写入的顺序问题，最好每炮独立保存，最后再合并。在MATLAB R2022b上如果运行报错，可以先检查是不是并行池占用内存过高，关掉并行池再跑，很多时候“Error 9”一类的报错都跟资源不足有关。

虚拟机上跑MATLAB的性能损失真的很大，一个物理机半分钟完成的模型，虚拟机能拖到好几分钟。如果这只是偶尔跑一次，忍一忍也就算了；如果要做参数扫描和批量模拟，建议直接装原生Linux或者Windows版本，不要再叠一层虚拟化。

5.4 验证程序正确性

程序跑通只是第一步，跑得对不对才是关键。我的调试顺序是这样的：先把各向异性参数全部设为0（ε=δ=0），让程序退化成各向同性弹性波正演。在均匀各向同性介质中，P波波前是一个完美的圆，SV波波前是另一个小圆，两个圆同心嵌套。如果快照上看到这两个圆，说明应力-速度更新和震源加载基本正确。

接着只在模型左半部分设置一个垂直界面，右侧改成另一种介质，看反射波和透射波能不能正确产生，界面处有没有数值伪影。最后再加入VTI各向异性参数，观察准P波波前是否从圆形变成椭圆形状。椭圆长轴方向和振幅大小随ε、δ的变化趋势，要和各向异性介质理论吻合。

在做这些验证的时候，记录每一组实验的固定时间步波场快照，累积成一组“基准图”。以后程序改动引入新bug，拿新的快照和基准图对比，能很快定位问题。这套方法论，比跑到一半对着数值发呆高效得多。

根据我个人实际跑这套程序的体会，最值得改进的地方是PML区域的代码结构。程序里把PML更新和主区更新混在一起，读起来费劲，调试时也不够灵活。如果你只是拿这套程序做基础研究或者课程作业，建议先保持原样跑通，再按自己的需求把PML部分拆出来单独维护。这个压缩包里的代码给了一个很好的起步框架，剩下的可以按你的模型和地质场景去扩展——比如加上自由地表边界、加入地形起伏，替换成TTI介质的倾斜对称轴，一步步把正演工具往前推。

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