企业官网建设流程全解析

简介：面向地震勘探与地球物理研究人员，该资料包提供了一套基于Matlab的VTI（垂直横向各向同性）介质地震波场数值模拟程序。程序实现了含PML吸收边界的波动方程有限差分求解，可输出任意时刻波场快照，适用于各向异性介质中纵波传播特征的分析与教学演示。压缩包共5个文件，其中4个为.m脚本（含主程序、差分系数计算及波场推进函数），另有1个.mat调色板数据文件，整体仅7KB，非常轻量易读，便于使用者二次修改与调试运行。已有521人学习下载，适合正在学习地震数值模拟或从事VTI介质正演研究的初学者参考。通过实际运行代码，可直观观察波前传播、反射与折射现象，并理解PML边界在压制人工反射中的作用，为后续地震成像与参数反演打下基础。

1. 把 VTI 介质 MATLAB 数值模拟跑通：关键是输出一张可信的波场快照

VTI 介质（垂直横向各向同性）的 MATLAB 数值模拟，大多数人拿到代码后第一件事，就是跑出一张能在指定时刻看出纵波、横波和各向异性特征的波场快照。这份压缩包对应的是 VTI 介质的二维波场模拟，核心是用有限差分裂解弹性波方程，输出不同时刻的应力或速度快照。对地震勘探、岩石物理建模和刚入波动方程数值模拟的从业者来说，它的价值在于省去从零搭框架的时间：改参数、看快照、再改参数，比对着教材推公式直观得多。适合已经会用 MATLAB 读脚本的人，也适合想在项目中先验证 ε、δ 会不会对波前产生可观测影响的人。如果你只是想把代码存进网盘，那就没必要往下读了；想在半小时内看到第一张波场图，继续。

2. VTI 介质波动方程的参数体系：先弄清 ε、δ 和刚度矩阵再动手

2.1 从各向同性到 VTI：为什么模型里要同时给 QP/QSv 速度和各向异性参数

各向同性模拟只需要一个纵波速度和一个横波速度，或者对应的拉梅常数，VTI 介质多了两个关键的各向异性参数：ε 和 δ。严格说还有 γ，但 γ 只影响纯横波 SH 波的分裂，二维 P-SV 波场模拟里用不到，所以大部分代码里看不到它的影子。真正影响纵波和 QSv 波的是 ε 和 δ：ε 描述水平方向纵波速度相对垂直方向纵波速度的差异，δ 描述近垂直方向的速度各向异性梯度。很多刚接触的人只给速度不给 ε、δ，跑出来的波场快照其实跟各向同性没有本质区别，这不是代码的问题，而是模型没有进入 VTI 状态。

场景上，如果目标是模拟页岩层或薄互层，Vp0 可以取 2800~3200 m/s，Vs0 取 1200~1400 m/s，ε 在 0.1 到 0.3 之间，δ 通常比 ε 小，取 0.05 到 0.15。这个组合来自实验室岩芯测量和井间观测，虽然不是唯一答案，但它是让波场快照“看起来真实”的常见起点。速度和各向异性参数不能混着给，比如有代码会把 ε 直接加到速度上，结果水平速度变成 Vp0×(1+ε)，这相当于把 Thomsen 参数当成了速度放大系数，后期走时曲线会对不上参考解。

2.2 由 Thomsen 参数计算刚度矩阵：切到 MATLAB 前的关键一步

VTI 介质的二维 P-SV 波场由弹性刚度矩阵中的四个非零分量决定：C11、C13、C33、C44。用 MATLAB 实现时，常见做法是先由 Vp0、Vs0、ρ 算出背景刚度，再用 ε、δ 修正水平和交叉分量。具体公式如下：

% 输入基本模型参数 Vp0 = 3000; % 垂直方向纵波速度, m/s Vs0 = 1400; % 垂直方向横波速度, m/s rho = 2200; % 密度, kg/m^3 epsilon = 0.2; % Thomsen epsilon delta = 0.1; % Thomsen delta C33 = rho * Vp0^2; C44 = rho * Vs0^2; C11 = C33 * (1 + 2 * epsilon); % 由 delta 的定义反解 C13 tmp = 2 * delta * C33 * (C33 - C44) + (C33 - C44)^2; C13 = sqrt(tmp) - C44; % 打印出来检查刚度矩阵数值范围 fprintf('C11=%.1e C13=%.1e C33=%.1e C44=%.1e\n', C11, C13, C33, C44);

这里最关键的是C13的反解。有些简化版代码会把 C13 直接写成C13 = sqrt(C11 * C33)，这是各向同性假设下的做法，用于 VTI 会让 δ 明显偏掉，波场快照的椭圆特征不真实。反解公式里带有C33 - C44，当 Vs0 接近 0 时会出现数值不稳定，所以流体层或软沉积层里要单独设置横波速度下限，常见的下限是 200~300 m/s。

参数说明：epsilon控制水平与垂直纵波速度的相对差，差值越大，波前面在水平方向拉得越长；delta控制近垂直方向的各向异性梯度，也是决定纵波在 45° 附近是否出现三叉戟趋势的主要参数。初次跑通时建议用epsilon=0.2, delta=0.1，先验证代码逻辑，再按目标地质体修改。用上面推荐参数算出来的一组参考刚度矩阵为：C11 约 2.77×10^10、C13 约 1.30×10^10、C33 约 1.98×10^10、C44 约 4.31×10^9，单位 Pa。如果你算出来的 C13 和 C33 差不多大，说明 δ 没进到刚度矩阵里。

2.3 差分格式与稳定性：网格步长、时间步长怎么定

这份资源大概率采用一阶速度-应力交错网格有限差分，因为它处理 VTI 介质比二阶位移方程更直接：应力分量和速度分量分布在半个网格间距上，更新公式不需要高阶导数，天然适合处理 C11、C13、C33、C44 随空间变化的情况。空间上一般用二阶精度，时间上也用二阶精度。要保证波场快照不出频散和数值爆炸，网格步长和时间步长必须同时满足两个约束。

% 建议参数设置示例 f0 = 25; % Ricker 子波主频, Hz fmax = 2.5 * f0; % 有效最高频率估为 2.5 倍主频 vmax = 3400; % 模型最大波速, m/s vmin = 1400; % 模型最小波速, m/s % 每波长网格数经验值是 12，偏少会出频散 dx = vmin / (fmax * 12); dz = dx; % 二维模拟保持空间步长各向同性 % 时间步长取库朗数 0.6 左右 dt = 0.6 * dx / (vmax * sqrt(2)); fprintf('dx=%.2f m, dt=%.3f ms\n', dx, dt * 1000);

逻辑说明：dx由最小速度和有效最高频率共同决定，而不是用主频。如果用主频去算，低频成分没问题，高频成分会在快照上出现频散，波前面附近拖着一串高频尾巴。dz = dx是二维模拟的默认选择，长宽比不一致会让波前面变成几何椭圆，这种假椭圆和介质各向异性造成的真椭圆非常难分辨。sqrt(2)来自二阶二维差分算子对时间稳定性的约束，具体系数会因差分格式不同而浮动，稳妥做法是先把库朗数压到 0.5。

拿到资源后，我一般先把dx按上面的公式重新算一遍，再和压缩包里自带参数对比。如果代码内置f0=20 Hz、dx=5 m，而你想把主频改到 40 Hz，那dx必须同步缩小，否则快照上的波前面会像蜈蚣一样伸出无数条高频毛刺。这个动作能提前挡掉 80% 的“快照看起来脏”问题。

3. 把压缩包里的 MATLAB 代码跑通：主脚本、震源与波场快照输出

3.1 拿到压缩包后的第一个动作：把路径理顺

这份资源解压后，典型结构是一个入口脚本加几个函数文件。入口脚本负责定义网格模型、调用初始化和时间循环；函数文件负责震源加载、边界吸收、波场更新。很多新手会直接双击某个函数文件然后点“运行”，结果报“未定义函数或变量”，这不是代码坏了，而是入口没找对。我拿到任何下载源码包都会先把 MATLAB 当前目录切到解压目录，再在命令行里输入which加入口名，确认入口脚本在搜索路径里。

另外注意 MATLAB 安装路径和工程路径不要带中文，也不要在脚本里暴露绝对路径。常见做法是在主脚本开头写一段自动路径处理：

% 自动获取当前脚本所在目录并加入路径 scriptDir = fileparts(mfilename('fullpath')); cd(scriptDir); addpath(genpath(scriptDir));

这段代码的作用是把压缩包解压后的所有子目录一次性加进搜索路径。fileparts(mfilename('fullpath'))在脚本以普通.m方式执行时能可靠返回文件所在目录；如果你用实时脚本.mlx跑，这个方法可能返回空值，所以建议启动后先切到脚本目录再运行入口文件。

3.2 Ricker 子波与加载方式：震源别只放在一个应力分量上

VTI 介质模拟常用的震源是 Ricker 子波，因为它频谱干净、带宽可控。子波公式如下：

function src = ricker_wavelet(t, t0, f0) % 标准 Ricker 子波 % t: 时间序列, s % t0: 延迟时间, s % f0: 主频, Hz src = (1 - 2 * (pi * f0 * (t - t0)).^2) .* exp(-(pi * f0 * (t - t0)).^2); end

主频f0决定子波的脉冲宽度，延迟t0一般取1.2 / f0，保证子波在 t=0 附近接近零，避免加载时产生突变。主频的选取要和网格步长联动，前面 2.3 节的dx公式就是围绕它来的。使用时的常见做法是把时间序列先算好，再按索引加载：

t = 0:dt:(nt-1) * dt; src = ricker_wavelet(t, 1.2 / f0, f0); % 震源位置，单位：网格点 isx = round(Nx / 2); isz = round(Nz / 3); % 放在偏上位置，避免边界影响 % 加载到正应力分量，注意乘以 dt sxx(isx, isz) = sxx(isx, isz) + dt * src(it); szz(isx, isz) = szz(isx, isz) + dt * src(it);

加载方式对 VTI 快照的影响很大。二维 P-SV 模拟里，爆炸源同时加在正应力分量 σxx 和 σzz 上，这样激发的波场以纵波为主，Sv 波对称出现在两侧；如果只在 σzz 上加载，波场会出现明显的横纵波混杂，快照看起来像两个错开的震源。把子波加载到应力分量上并乘dt，是速度-应力方程离散化的标准做法。如果你发现快照里的子波被“多积了一次分”，看上去低频漂移明显，多半就是没有乘dt，或者在速度分量而不是应力分量上加载了同一个源。

3.3 速度-应力交错网格主循环：五分钟看懂的更新式

有了模型参数、震源和边界，核心的时间循环就是反复更新三个应力分量和两个速度分量。VTI 介质和各项同性不同，应力-应变关系要用前面算出的 C11、C13、C33、C44。主循环里的关键更新式通常是这个形态：

% 速度更新，使用交错网格空间差分 vx(2:end-1, 2:end-1) = vx(2:end-1, 2:end-1) + ... dt / rho(2:end-1, 2:end-1) .* ... ((sxx(2:end, 2:end-1) - sxx(1:end-1, 2:end-1)) / dx + ... (sxz(2:end-1, 2:end) - sxz(2:end-1, 1:end-1)) / dz); vz(2:end-1, 2:end-1) = vz(2:end-1, 2:end-1) + ... dt / rho(2:end-1, 2:end-1) .* ... ((sxz(2:end, 2:end-1) - sxz(1:end-1, 2:end-1)) / dx + ... (szz(2:end-1, 2:end) - szz(2:end-1, 1:end-1)) / dz); % 应力更新，代入 VTI 刚度矩阵 sxx(2:end-1, 2:end-1) = sxx(2:end-1, 2:end-1) + dt * ... (C11(2:end-1, 2:end-1) .* (vx(2:end, 2:end-1) - vx(1:end-1, 2:end-1)) / dx + ... C13(2:end-1, 2:end-1) .* (vz(2:end-1, 2:end) - vz(2:end-1, 1:end-1)) / dz); szz(2:end-1, 2:end-1) = szz(2:end-1, 2:end-1) + dt * ... (C13(2:end-1, 2:end-1) .* (vx(2:end, 2:end-1) - vx(1:end-1, 2:end-1)) / dx + ... C33(2:end-1, 2:end-1) .* (vz(2:end-1, 2:end) - vz(2:end-1, 1:end-1)) / dz); sxz(2:end-1, 2:end-1) = sxz(2:end-1, 2:end-1) + dt * ... C44(2:end-1, 2:end-1) .* ... ((vz(2:end, 2:end-1) - vz(1:end-1, 2:end-1)) / dx + ... (vx(2:end-1, 2:end) - vx(2:end-1, 1:end-1)) / dz);

参数说明：这里的sxx(2:end,2:end-1) - sxx(1:end-1,2:end-1)是对 x 方向偏导的中心差分，但在交错网格上两个数组错开半个网格步长，导数中心正好落在另一个变量的位置上。C11、C13、C33、C44 既可以是常数，也可以是随空间变化的矩阵；如果是变化模型，要用逐点乘法.*，写成C11 * (vx(...))会被 MATLAB 当成矩阵乘法直接报维度错。rho同理，速度更新时要在同一位置取密度值。

这个循环是资源里最容易出问题的地方。如果代码里把 C13 和 C33 的位置搞混，快照的波前面会有明显的“菱形变形”；如果全部换成各向同性公式，快照又会在长偏移距处偏离真实 VTI 波场。我第一次跑通时把 C11、C13 同时用成了 C33，快照上横波速度明显偏快，排查了半天才发现少给了一个epsilon。

3.4 波场快照提取与画图：imagesc 加合适的色标

波场快照本质上是某一时刻某个分量的空间分布。资源里常见的输出方式是每 N 步把vz或sxx存到二维数组里，再绘制出来。画图时建议用imagesc，因为它按像素显示，没有等高线的平滑误差：

figure; imagesc(x, z, squeeze(vz(:,:,it))); axis ij; % 让 z 轴方向向下，符合地震剖面习惯 axis equal tight; colormap(gray); caxis([-0.02 0.02]); % 动态范围固定，避免每帧自动缩放 xlabel('Distance (m)'); ylabel('Depth (m)'); title(sprintf('Wavefield snapshot at t=%.2f ms', t(it) * 1000));

这里axis ij很关键。默认imagesc的图像纵坐标是第一维，倒过来会让快照看起来上下颠倒，很多人第一次出图发现“波往空气里跑”就是这个原因。caxis建议固定为对称区间，否则每帧自动定标会让弱纵横波忽隐忽现。最近几个 MATLAB 版本把caxis改成了推荐clim，如果你的版本报caxis警告，直接改成clim([-0.02 0.02])即可。

如果要批量保存快照，不要把图一帧帧截图，建议直接把波场数组写成.mat文件，后面统一出图：

save(sprintf('snapshot_%04d.mat', it), 'vz', 'x', 'z', 't');

这样做的好处是后期换色标、换动态范围、换分量的时候不用重跑模拟。

4. 边界条件与网格参数：PML 吸收层和网格精度怎么搭配

4.1 三种常见边界处理：反射、海绵、PML 的取舍

VTI 介质模拟的边界处理决定了快照里有没有假的边界反射。最简单的处理是刚性边界，也就是不处理，波打到边界后反射回来，快照在晚时刻会看到一组以边界为轴的二次圆弧，看起来很饱满，其实全是虚假信号。如果只是验证短时间内的波前面形状，可以把模拟时间压缩到波还没到边界之前，但这种做法大大浪费模型空间，也限制了你观察远处的各向异性差异。

海绵边界是中间方案，在边界附近把波场按权重吸收，但吸收系数分布不当会在边界产生新的反射。更可靠的是 PML，它通过在边界区域修改波动方程，在频域上实现无反射吸收。PML 参数再多也不离谱，VTI 介质的 PML 与各向同性 PML 在实现上基本一致，因为 PML 作用在波场分量上，而不是刚度矩阵上。也就是说，如果你已经在各向同性代码里写好了 PML，直接迁移到 VTI 时不需要改吸收层公式，只需要保证入射到 PML 的分量是完整的。

4.2 PML 衰减系数和厚度的经验值

PML 最常翻车的两个参数是厚度和衰减系数。厚度太薄，低频波会直接穿过，看起来像边界没有被吸收；衰减系数太大，PML 内部的反射会变成波场快照里的“亮边”。我常用的经验值是：PML 厚度取 15 到 20 个网格点，衰减系数按下面这个公式随位置渐变：

% 构造 PML 区域的衰减系数，d0 是最大衰减 L = 20; % PML 厚度，网格数 frac = (1:L) / L; % 从边界到内部线性增长 d0 = -3 * vmax / (2 * L * dx) * log(1e-3); % 目标反射系数 1e-3 d = d0 * frac.^2; % 平方渐变，减弱 PML 内部反射

参数说明：vmax是介质最大波速，dx是空间步长。frac.^2让衰减系数从边界向内部平滑过渡，而不是突然从 0 跳到 d0，这样能大幅减少 PML 内边界处的反射。log(1e-3)里的反射系数 R=1e-3 是常见目标，想更强就改成 1e-4，但注意d0变大后时间步长要相应缩小，否则会不稳定。

在波场快照里判断 PML 好坏，不用看振幅值，而是看边界附近是否出现“亮环”。如果出现，先不要调衰减系数，把 PML 厚度加厚 5 个网格点再看。如果波场在 PML 内部整体震荡，说明d0过大或库朗数偏大，把时间步长压到原来的一半通常能压下去。

4.3 每波长网格数：快照不出现频散的最小要求

网格步长和 PML 是同一条链上的两个环节。如果dx太大，即便 PML 完美，快照也会在高频位置出现频散：波前面后面跟着一串同相位的“毛刺”。在 VTI 介质里，Sv 波速度比 P 波低不少，频散往往先出现在 Sv 波上，因为同频率下 S 波波长更短。

经验值是每最短波长至少 10 到 15 个网格点。所谓最短波长，用vmin / fmax，其中vmin是模型里所有波速的最小值，包括横波。很多人只考虑 P 波速度，结果 P 波快照看着干净，Sv 波像乱码。VTI 模拟里 Sv 波是不可忽视的物理信号，在页岩模型中 Sv 波往往携带各向异性信息，所以建议把最小横波速度也放进dx计算：

vmin_s = min(Vs0(:)); % 最小横波速度 dx = vmin_s / (2.5 * f0 * 15); % 每波长 15 个网格，留出余量

如果你发现快照里只有 P 波干净、Sv 波毛糙，第一步不是加滤波，而是减小dx重新算。滤波会把物理上的 Sv 波和高频噪声一起抹掉，违背了用数值模拟看波场的初衷。网格总数也随之上升，二维模型从 500×500 变成 800×800，MATLAB 内存占用大约翻三倍，但换来的波前平滑度值得这个代价。

5. VTI 仿真避坑指南：从波场快照反推的五个常见问题

5.1 快照满屏棋盘状干扰：网格上变量位置没配对

现象：波场快照上同时出现间隔一个网格点的黑白相间干扰，整个波前像蒙了一层棋盘格。原因：交错网格有限差分对变量的位置有严格要求。标准做法是vx位于整网格，vz位于半网格，sxx、szz与vz对齐，而sxz落在网格中心。如果代码里三个应力分量全定义在整点上，速度更新时会出现响应为零的寄生模式，也就是棋盘状噪声。解决：检查主循环里的差分表达式。正确的交错网格导数(u(i+1)-u(i))/dx对应导数位置在i+1/2；如果你看到(u(i+1)-u(i-1))/(2*dx)这类中心差分出现在速度-应力更新里，那多半不是交错网格，需要整套调整变量布局。这个坑我在从位移方程转速度-应力方程时踩过，改完数组错位后快照立刻干净了。

5.2 波前面方向对、速度不对：ε 和 δ 的权重被忽略

现象：波前面沿水平和垂直方向的传播距离不同，说明各向异性特征有了，但具体走时和入射角对不上参考解，尤其大偏移距处偏快或偏慢。原因：VTI 介质的相速度和群速度由 ε、δ 共同控制。如果代码把 δ 忽略了，只保留 ε，得到的是椭圆各向异性近似；对弱各向异性尚可，但对页岩这种 ε=0.2、δ=0.1 的介质，会在 45° 附近相位速度明显偏低，波前面出现错误的内凹。解决：回到刚度矩阵公式，用 2.2 节的C13反解公式重新算一遍。先做个 sanity check：把 ε、δ 都设成 0，快照应退化到各向同性圆形波前；再把 δ 设成 0、ε 保持 0.2，快照应出现椭圆波前；最后同时给 ε、δ，看 45° 附近是否有三叉戟趋势。这三步能快速定位是公式还是参数写错。

5.3 边界附近振幅炸掉：PML 衰减系数取值太大

现象：PML 内部不是安静吸收，而是出现随时间增长的大振幅震荡，波场快照在边界处像着火一样。原因：PML 衰减系数d0取值过大，导致 PML 区域的计算波数变成复数，反射异常放大入射波。另一个原因是使用 PML 时时间步长没按 PML 区域的最大速度重新核算库朗数。解决：把d0降一半，把库朗数压到 0.4 以下。PML 厚度从 20 个网格点开始，不要一上来就用 10 个点追求省事。如果资源里的 PML 是复频率移形式，还要检查频率参数是否匹配子波主频；不匹配会在低频段出现慢速拖尾反射。

5.4 快照看起来和各向同性没区别：刚度矩阵用成了椭圆近似

现象：无论怎么改 ε、δ，波场快照都只有速度大小变化，形状始终保持圆形。原因：这是最隐蔽的一种。代码里可能用了C13 = C11 - 2*C44或C13 = C11 - 2*C66这类基于各向同性的关系式，这类公式在 VTI 介质里不满足 Christoffel 方程，导致刚度矩阵退化为各向同性，只有速度缩放，没有方向性。解决：把刚度矩阵的四个分量全部打印出来，检查C11 / C33是否等于1 + 2*epsilon，以及delta是否影响C13。如果C13对delta完全不敏感，代码里的刚度矩阵就是假的 VTI。换成标准 Thomsen 参数反解后，快照立刻有方向特征。

5.5 波场里 Sv 波极性混乱：震源的分量加载方式不对

现象：近震源处的 Sv 波图案不对称，或者一侧波前面极性反转，和理论波形对不上。原因：爆炸源应在两个正应力分量上同时加载等幅子波，才能激发对称的纵波和 Sv 波。如果资源里为了省事只在szz上加载，会产生类似单力源的辐射图案，Sv 波在震源两侧的极性完全不对称。解决：把震源加载改成sxx(isx,isz) += dt*src; szz(isx,isz) += dt*src;，并确认加载位置在网格整点上。快照检查时看 Sv 波在两个方向上的相位是否相反、振幅是否对称，这是一个非常灵敏的验证指标。

6. 进阶用法：用波形记录和参数扫描验证你的 VTI 模型

波场快照只能看形状，走时曲线才能验证速度。常见做法是在模型右侧布置一排检波器，记录vz分量在不同时刻到达的波峰时间，画出偏移距与走时的关系。VTI 介质中，近偏移距走时主要受垂直速度控制，远偏移距走时受 ε 影响。如果快照显示波前面早到了，但走时曲线对不上，优先检查刚度矩阵而不是网格参数。走时提取可以这样写在主循环后面：

% 在时间循环中收集某检波器位置的 vz 波形 rep = 1:Nx-200; % 右侧检波器网格范围 seis(:, ir) = vz(rep, isz); % 保存为二维记录，列对应偏移距

循环结束后用findpeaks找到第一波峰对应的时刻，和理论走时对比。这一步能很快看出 ε 和 δ 的差异：δ=0.1 会让近距离走时相对垂直速度偏差 0.5% 到 1%，这个量在快照上看不清，在走时曲线上一清二楚。

进阶用法是把资源里的主脚本改写成函数，输入是一组 ε/δ 组合，输出是快照或走时。用parfor可以同时跑 8 到 12 组参数。扫描时建议先固定 δ，把 ε 从 0 到 0.4 按 0.05 步进扫一遍，得到一组有递进关系的快照，比单独看一张更能理解各向异性波场的行为。做参数扫描时保存文件名一定要带参数，我习惯写snapshot_eps020_delta010_it300.mat，否则连续跑几十组后，你根本分不清哪张图是哪个模型。另外，每组参数的子波主频保持一致，否则对比会因频散差异失去意义。

其实这一段教训来自我第一次用这份资源时：上来就调 ε、δ，快照看着“差不多像样”，结果走时曲线和同样参数的理论值差了 3%。从那以后，我每次跑 VTI 模拟都会先保留一套各向同性参考模型，把参考波场跑出来，再叠加各向异性模型，快照和走时一起对比。这个过程虽然折腾，但能拦下大量隐蔽的参数错误。希望帮到你。

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