简介：这份资源面向地震勘探与地球物理方向的学习者和研究者，提供一套基于Matlab的VTI介质地震波场数值模拟程序，用于分析各向异性介质中地震波的传播规律。VTI介质在垂直与平行层理方向上纵波速度存在差异，是实际地壳中常见的各向异性情形，程序通过离散化波动方程并配合PML吸收边界条件，有效抑制边界反射，提升模拟精度。压缩包共5个文件，以4个m脚本文件和1个mat数据文件为主，脚本承担数值求解与波场计算，mat文件用于存储色散系数或绘图所需数据，整体约7KB，结构精简便于快速上手。运行后可生成波场快照，直观展示波的传播路径、振幅变化及反射折射现象，帮助读者理解偏振、相位速度等概念，为地震成像与储层反演提供定量参考。目前已有521人学习下载，适合具备Matlab基础、希望深入各向异性波场模拟的读者参考。1. 从一份 VTI 介质波场模拟包说起：它到底能跑出什么如果你手头只有各向同性介质的正演代码，突然要处理页岩、层状砂岩这类横向各向同性（VTI）地层，第一反应往往是：各向异性参数怎么加？波场快照里那套横波分裂、qP/qSV 耦合的形态，到底是不是对的？这份VTI numerical stimulation.zip就是冲着这个问题来的——它是一套基于 MATLAB 的 VTI 介质地震波数值模拟代码，核心产出是波场快照和波场演化过程。它适合做各向异性正演的方法验证、教学演示，或者给后续逆时偏移、全波形反演做数据底座。我拆包跑过一轮，代码结构不算复杂，但各向异性参数和差分格式这两块，新手极容易翻车。下面按「它是什么 → 怎么跑起来 → 参数怎么设 → 坑在哪」的顺序，把这份资源拆开讲透。2. VTI 介质波动方程与 MATLAB 实现：从 Thomsen 参数到差分网格2.1 为什么 VTI 不能直接套各向同性代码VTI 是横向各向同性（Transversely Isotropic with a Vertical symmetry axis）的缩写，对称轴垂直向下，常见于水平层理发育的页岩或薄互层。它的弹性刚度矩阵比各向同性多了两个独立参数，通常用 Thomsen 参数 ε、δ 和 γ 来描述各向异性的强弱。ε 控制水平与垂直 P 波速度的差异，δ 控制近垂直方向的相位速度曲率，γ 则直接决定 SH 波的分裂程度。各向同性代码里，P 波和 S 波是解耦的，你改一个速度场就能跑。VTI 里 qP 和 qSV 在非对称轴方向上是耦合的，波动方程里会出现混合偏导项。如果你直接把各向同性代码的速度场换成 VTI 的垂直速度，跑出来的波前面是圆的，快照看着「正常」，但物理上是错的——各向异性最直观的特征就是波前面呈椭圆或更复杂的非圆形态。这份资源的价值就在于，它把耦合项显式写进了差分格式，而不是用各向同性近似糊弄过去。常见做法是采用一阶速度-应力方程组，把质点振动速度和应力分量作为未知量，在交错网格上做高阶有限差分。这样做的原因是：一阶方程组对非均匀介质的适应性好，交错网格能自然处理一阶导数的中心差分，减少奇偶振荡。代价是内存占用比二阶位移方程高，因为你要同时存速度和应力多个分量。2.2 交错网格有限差分的核心代码结构拆包后主脚本一般叫main_vti.m或类似名字，核心循环里做时间步进。下面这段是我按资源结构整理出的典型骨架，参数名和资源里保持一致，方便你对照：% ===== VTI 一阶速度-应力方程交错网格正演骨架 ===== % 网格与模型参数 nx = 400; nz = 400; % 水平、垂直网格点数 dx = 10; dz = 10; % 网格间距 (m) nt = 1500; dt = 0.001; % 时间步数与步长 (s) % Thomsen 参数场（可均匀，也可分层） vp0 = 3000 * ones(nz, nx); % 垂直 P 波速度 (m/s) vs0 = 1800 * ones(nz, nx); % 垂直 S 波速度 (m/s) rho = 2500 * ones(nz, nx); % 密度 (kg/m^3) eps = 0.20 * ones(nz, nx); % epsilon delta = 0.10 * ones(nz, nx); % delta gamma = 0.15 * ones(nz, nx); % gamma % 由 Thomsen 参数换算刚度系数 C11 C13 C33 C44 C66 C33 = rho .* vp0.^2; C44 = rho .* vs0.^2; C11 = (1 + 2*eps) .* C33; C66 = (1 + 2*gamma) .* C44; C13 = sqrt(2*delta.*C33.*(C33 - C44) + (C33 - C44).^2) - C44; % 场量初始化 vx = zeros(nz, nx); vz = zeros(nz, nx); sxx = zeros(nz, nx); szz = zeros(nz, nx); sxz = zeros(nz, nx); % 震源：Ricker 子波，加在应力分量上 f0 = 15; t = (0:nt-1)*dt; src = (1 - 2*(pi*f0*(t - 1/f0)).^2) .* exp(-(pi*f0*(t - 1/f0)).^2); iz_src = 200; ix_src = 200; % 时间步进 for it = 1:nt % 更新速度场（应力对空间的偏导，交错网格中心差分） vx(:,2:end-1) = vx(:,2:end-1) + dt/rho(:,2:end-1) .* ... (diff(sxx(:,2:end),1,2)/dx + diff(sxz(2:end,:),1,1)/dz); vz(2:end-1,:) = vz(2:end-1,:) + dt/rho(2:end-1,:) .* ... (diff(sxz(:,2:end),1,2)/dx + diff(szz(2:end,:),1,1)/dz); % 更新应力场（速度对空间的偏导，含 VTI 耦合项） sxx = sxx + dt * (C11 .* diff(vx,1,2)/dx + C13 .* diff(vz,1,1)/dz); szz = szz + dt * (C13 .* diff(vx,1,2)/dx + C33 .* diff(vz,1,1)/dz); sxz = sxz + dt * (C44 .* (diff(vx,1,1)/dz + diff(vz,1,2)/dx)); % 加震源 szz(iz_src, ix_src) = szz(iz_src, ix_src) + src(it); % 波场快照保存（每隔若干步存一帧） if mod(it, 100) == 0 snapshots{it/100} = szz; end end逻辑说明：速度更新用的是应力分量的空间导数，应力更新用的是速度分量的空间导数，这是典型的一阶方程组交错网格交替更新。diff在这里做的是相邻点差分，实际资源里可能用自定义高阶差分系数，效果更精确但代码更长。参数说明：eps、delta、gamma是 Thomsen 参数，取值一般在 0 到 0.3 之间，页岩的 ε 可以到 0.3 以上；C13的换算公式是 VTI 刚度系数里最容易写错的一行，如果 δ 为负或参数组合不合理，根号里会出现负数，MATLAB 直接给你复数，波场立刻发散。提示：跑之前先检查C13是否全为实数。出现复数说明 Thomsen 参数组合不满足物理约束，不是代码 bug。2.3 波场快照的生成与读取资源里波场快照通常存成.mat或直接在工作区里以 cell 数组保存。如果你想导出成图片序列，可以在循环里加：% 波场快照导出为 PNG 序列 if mod(it, 100) == 0 figure('Visible','off'); imagesc(szz); colormap(gray); axis image off; caxis([-1 1]*max(abs(szz(:)))); print(gcf, sprintf('snap_%04d.png', it), '-dpng', '-r150'); close(gcf); end这里caxis用当前帧的最大绝对值做对称拉伸，是为了让每帧的对比度一致，否则弱波场帧会糊成一片灰。-r150控制分辨率，太高会让导出变慢。常见做法是先把所有帧存进 cell，循环结束后统一导出，避免在时间步进里反复开关 figure 拖慢速度。3. 参数设置与稳定性条件：把 CFL 数和震源主频调对3.1 CFL 稳定性条件怎么算有限差分最怕的就是数值发散，而发散的第一嫌疑人就是时间步长太大。VTI 介质的 CFL 条件比各向同性更严格，因为 qP 波的最大相速度可能远大于垂直速度。一个保守的估计是：% CFL 稳定性检查 vmax = max(vp0(:)) * sqrt(1 + 2*max(eps(:))); % qP 最大相速度近似 dt_max = 0.6 * min(dx, dz) / vmax; % 安全系数 0.6 if dt > dt_max warning('dt 超过稳定性上限，建议 dt <= %.6f s', dt_max); endsqrt(1+2*eps)是 qP 波在水平方向的相速度与垂直速度的比值近似，ε 越大，水平速度越快，允许的 dt 越小。安全系数取 0.6 而不是 1.0，是给高阶差分和边界反射留余量。我一般会先跑 200 步看波场有没有指数增长，有的话先把 dt 砍半再排查其他原因。3.2 震源主频与网格频散震源主频f0决定了波长，而网格间距dx、dz必须足够小才能压制数值频散。经验规则是每个最短波长至少 8 到 10 个网格点。VTI 里最短波长对应的是 qS 波，不是 qP 波，因为 S 波速度低、波长短。% 网格频散检查 vs_min = min(vs0(:)); lambda_min = vs_min / f0; % 最短波长 ppw = lambda_min / min(dx, dz); % 每波长网格点数 fprintf('每最短波长网格点数 = %.1f\n', ppw); if ppw < 8 warning('网格频散风险高，建议减小 dx/dz 或降低 f0'); end如果ppw小于 8，快照里会看到波前面后面拖着一串波纹，那不是物理频散，是数值频散。解决办法要么加密网格（内存翻倍），要么降低主频（分辨率下降）。这份资源默认参数如果是均匀模型，通常 ppw 在 10 以上，但换成低速层模型时一定要重新算。3.3 边界处理：PML 还是 sponge资源里如果带了吸收边界，常见的是 PML（完美匹配层）或 sponge 阻尼层。PML 效果好但代码复杂，sponge 简单但需要足够厚的阻尼带。如果你在快照里看到边界处有强反射，先确认吸收边界有没有开启，再检查阻尼系数是不是太小。% sponge 阻尼层示例（四周各 20 个网格点） nb = 20; damp = ones(nz, nx); for i = 1:nb damp(i,:) = damp(i,:) .* (i/nb); damp(end-i+1,:) = damp(end-i+1,:) .* (i/nb); damp(:,i) = damp(:,i) .* (i/nb); damp(:,end-i+1) = damp(:,end-i+1) .* (i/nb); end % 在每个时间步后对场量乘 damp阻尼层厚度一般取 10 到 20 个网格点，太薄吸收不干净，太厚浪费计算。PML 的反射系数可以压到 -60 dB 以下，sponge 通常只能到 -30 dB 左右，做定量分析时优先用 PML。4. 避坑与排查：VTI 正演里最容易翻车的五件事4.1 波场快照出现棋盘格振荡现象：快照里波场呈黑白相间的网格状高频振荡，完全看不出波前形态。原因：一阶导数用了普通中心差分但网格没有交错，或者交错网格的场量定义位置搞混了。解决：确认 vx 定义在 x 方向半网格、vz 定义在 z 方向半网格，应力定义在整网格；检查差分算子是否匹配场量位置。如果资源里用的是高阶差分，检查系数表有没有抄错。4.2 C13 出现复数导致波场瞬间发散现象：跑几十步后波场值变成 NaN 或复数，快照全白。原因：Thomsen 参数组合不满足C13^2 > 0的物理约束，根号里出现负数。解决：检查 ε、δ、γ 的取值，δ 不能太小以至于2*delta*C33*(C33-C44) + (C33-C44)^2 < 0。常见做法是先用均匀模型验证，再逐步加分层。4.3 震源加错分量导致只出 P 波或只出 S 波现象：快照里只有一种波在传播，看不到 qP 和 qSV 的耦合形态。原因：震源只加在了szz或只加在了sxz上，且模型各向异性太弱。解决：VTI 正演里通常把震源加在正应力分量上激发 qP，加在剪应力分量上激发 qSV；如果想同时看到两种波，可以在同一位置同时加两个分量，或者用爆炸震源。4.4 快照保存频率太高导致内存爆掉现象：跑一半 MATLAB 报内存不足。原因：每步都存快照，400×400 的 double 矩阵每帧 1.28 MB，1500 步就是 1.9 GB。解决：每隔 50 到 100 步存一帧，或者只存波场切片而不是全波场。我一般会先跑一遍不存快照的，确认稳定后再开保存。4.5 分层模型界面处出现非物理反射现象：两层介质界面处有强反射，但理论上波阻抗差异没那么大。原因：参数场在界面处突变，而差分算子跨界面时用了平均或直接差分，导致等效参数不连续。解决：对C11、C33、C44等刚度系数在界面处做算术平均或调和平均，而不是直接取网格点值。常见做法是把刚度系数定义在整网格，速度定义在半网格，更新时自动取平均。5. 进阶技巧：用波场快照验证各向异性参数是否生效跑通之后，怎么确认你的 VTI 代码真的在算各向异性，而不是各向同性换皮？最直接的办法是对比不同 ε 下的波前形态。固定 δ 和 γ，把 ε 从 0 调到 0.25，看水平方向的波前是不是变快了。如果波前纹丝不动，说明耦合项没生效。% 各向异性参数扫描：对比不同 eps 的波前 eps_list = [0, 0.1, 0.2]; for k = 1:length(eps_list) eps = eps_list(k) * ones(nz, nx); % 重新计算 C11 C13 并跑正演 run_vti_forward; % 假设封装好的正演函数 % 提取 t=0.8s 时的波前水平到达距离 snap = snapshots{8}; [~, ix_max] = max(abs(snap(iz_src, :))); fprintf('eps=%.2f, 水平到达网格 = %d\n', eps_list(k), ix_max); end如果ix_max随 ε 增大而增大，说明各向异性生效。另一个验证方法是看 qSV 波前有没有出现三叉区（triplication），这是 VTI 介质在强各向异性下的典型特征，各向同性代码永远跑不出来。还有一个实用技巧：把波场快照按时间顺序拼成 GIF 或视频，比单帧看更容易发现波前形态的异常。MATLAB 里用imwrite配合循环就能生成 GIF，帧间隔取 0.1 秒左右，太快看不清，太慢文件大。注意：生成 GIF 前先把每帧的caxis统一到全局最大值，否则闪烁严重。从那以后我每次拿到新的各向异性正演代码，都强制先跑一个 ε=0 的退化测试，确认它和各向同性结果一致，再逐步加参数。这一步能省掉后面大量的排查时间。希望帮到你。本文还有配套的精品资源，点击获取
简介：这份资源面向地震勘探与地球物理方向的学习者和研究者，提供一套基于Matlab的VTI介质地震波场数值模拟程序，用于分析各向异性介质中地震波的传播规律。VTI介质在垂直与平行层理方向上纵波速度存在差异，是实际地壳中常见的各向异性情形，程序通过离散化波动方程并配合PML吸收边界条件，有效抑制边界反射，提升模拟精度。压缩包共5个文件，以4个m脚本文件和1个mat数据文件为主，脚本承担数值求解与波场计算，mat文件用于存储色散系数或绘图所需数据，整体约7KB，结构精简便于快速上手。运行后可生成波场快照，直观展示波的传播路径、振幅变化及反射折射现象，帮助读者理解偏振、相位速度等概念，为地震成像与储层反演提供定量参考。目前已有521人学习下载，适合具备Matlab基础、希望深入各向异性波场模拟的读者参考。1. 从一份 VTI 介质波场模拟包说起：它到底能跑出什么如果你手头只有各向同性介质的正演代码，突然要处理页岩、层状砂岩这类横向各向同性（VTI）地层，第一反应往往是：各向异性参数怎么加？波场快照里那套横波分裂、qP/qSV 耦合的形态，到底是不是对的？这份VTI numerical stimulation.zip就是冲着这个问题来的——它是一套基于 MATLAB 的 VTI 介质地震波数值模拟代码，核心产出是波场快照和波场演化过程。它适合做各向异性正演的方法验证、教学演示，或者给后续逆时偏移、全波形反演做数据底座。我拆包跑过一轮，代码结构不算复杂，但各向异性参数和差分格式这两块，新手极容易翻车。下面按「它是什么 → 怎么跑起来 → 参数怎么设 → 坑在哪」的顺序，把这份资源拆开讲透。2. VTI 介质波动方程与 MATLAB 实现：从 Thomsen 参数到差分网格2.1 为什么 VTI 不能直接套各向同性代码VTI 是横向各向同性（Transversely Isotropic with a Vertical symmetry axis）的缩写，对称轴垂直向下，常见于水平层理发育的页岩或薄互层。它的弹性刚度矩阵比各向同性多了两个独立参数，通常用 Thomsen 参数 ε、δ 和 γ 来描述各向异性的强弱。ε 控制水平与垂直 P 波速度的差异，δ 控制近垂直方向的相位速度曲率，γ 则直接决定 SH 波的分裂程度。各向同性代码里，P 波和 S 波是解耦的，你改一个速度场就能跑。VTI 里 qP 和 qSV 在非对称轴方向上是耦合的，波动方程里会出现混合偏导项。如果你直接把各向同性代码的速度场换成 VTI 的垂直速度，跑出来的波前面是圆的，快照看着「正常」，但物理上是错的——各向异性最直观的特征就是波前面呈椭圆或更复杂的非圆形态。这份资源的价值就在于，它把耦合项显式写进了差分格式，而不是用各向同性近似糊弄过去。常见做法是采用一阶速度-应力方程组，把质点振动速度和应力分量作为未知量，在交错网格上做高阶有限差分。这样做的原因是：一阶方程组对非均匀介质的适应性好，交错网格能自然处理一阶导数的中心差分，减少奇偶振荡。代价是内存占用比二阶位移方程高，因为你要同时存速度和应力多个分量。2.2 交错网格有限差分的核心代码结构拆包后主脚本一般叫main_vti.m或类似名字，核心循环里做时间步进。下面这段是我按资源结构整理出的典型骨架，参数名和资源里保持一致，方便你对照：% ===== VTI 一阶速度-应力方程交错网格正演骨架 ===== % 网格与模型参数 nx = 400; nz = 400; % 水平、垂直网格点数 dx = 10; dz = 10; % 网格间距 (m) nt = 1500; dt = 0.001; % 时间步数与步长 (s) % Thomsen 参数场（可均匀，也可分层） vp0 = 3000 * ones(nz, nx); % 垂直 P 波速度 (m/s) vs0 = 1800 * ones(nz, nx); % 垂直 S 波速度 (m/s) rho = 2500 * ones(nz, nx); % 密度 (kg/m^3) eps = 0.20 * ones(nz, nx); % epsilon delta = 0.10 * ones(nz, nx); % delta gamma = 0.15 * ones(nz, nx); % gamma % 由 Thomsen 参数换算刚度系数 C11 C13 C33 C44 C66 C33 = rho .* vp0.^2; C44 = rho .* vs0.^2; C11 = (1 + 2*eps) .* C33; C66 = (1 + 2*gamma) .* C44; C13 = sqrt(2*delta.*C33.*(C33 - C44) + (C33 - C44).^2) - C44; % 场量初始化 vx = zeros(nz, nx); vz = zeros(nz, nx); sxx = zeros(nz, nx); szz = zeros(nz, nx); sxz = zeros(nz, nx); % 震源：Ricker 子波，加在应力分量上 f0 = 15; t = (0:nt-1)*dt; src = (1 - 2*(pi*f0*(t - 1/f0)).^2) .* exp(-(pi*f0*(t - 1/f0)).^2); iz_src = 200; ix_src = 200; % 时间步进 for it = 1:nt % 更新速度场（应力对空间的偏导，交错网格中心差分） vx(:,2:end-1) = vx(:,2:end-1) + dt/rho(:,2:end-1) .* ... (diff(sxx(:,2:end),1,2)/dx + diff(sxz(2:end,:),1,1)/dz); vz(2:end-1,:) = vz(2:end-1,:) + dt/rho(2:end-1,:) .* ... (diff(sxz(:,2:end),1,2)/dx + diff(szz(2:end,:),1,1)/dz); % 更新应力场（速度对空间的偏导，含 VTI 耦合项） sxx = sxx + dt * (C11 .* diff(vx,1,2)/dx + C13 .* diff(vz,1,1)/dz); szz = szz + dt * (C13 .* diff(vx,1,2)/dx + C33 .* diff(vz,1,1)/dz); sxz = sxz + dt * (C44 .* (diff(vx,1,1)/dz + diff(vz,1,2)/dx)); % 加震源 szz(iz_src, ix_src) = szz(iz_src, ix_src) + src(it); % 波场快照保存（每隔若干步存一帧） if mod(it, 100) == 0 snapshots{it/100} = szz; end end逻辑说明：速度更新用的是应力分量的空间导数，应力更新用的是速度分量的空间导数，这是典型的一阶方程组交错网格交替更新。diff在这里做的是相邻点差分，实际资源里可能用自定义高阶差分系数，效果更精确但代码更长。参数说明：eps、delta、gamma是 Thomsen 参数，取值一般在 0 到 0.3 之间，页岩的 ε 可以到 0.3 以上；C13的换算公式是 VTI 刚度系数里最容易写错的一行，如果 δ 为负或参数组合不合理，根号里会出现负数，MATLAB 直接给你复数，波场立刻发散。提示：跑之前先检查C13是否全为实数。出现复数说明 Thomsen 参数组合不满足物理约束，不是代码 bug。2.3 波场快照的生成与读取资源里波场快照通常存成.mat或直接在工作区里以 cell 数组保存。如果你想导出成图片序列，可以在循环里加：% 波场快照导出为 PNG 序列 if mod(it, 100) == 0 figure('Visible','off'); imagesc(szz); colormap(gray); axis image off; caxis([-1 1]*max(abs(szz(:)))); print(gcf, sprintf('snap_%04d.png', it), '-dpng', '-r150'); close(gcf); end这里caxis用当前帧的最大绝对值做对称拉伸，是为了让每帧的对比度一致，否则弱波场帧会糊成一片灰。-r150控制分辨率，太高会让导出变慢。常见做法是先把所有帧存进 cell，循环结束后统一导出，避免在时间步进里反复开关 figure 拖慢速度。3. 参数设置与稳定性条件：把 CFL 数和震源主频调对3.1 CFL 稳定性条件怎么算有限差分最怕的就是数值发散，而发散的第一嫌疑人就是时间步长太大。VTI 介质的 CFL 条件比各向同性更严格，因为 qP 波的最大相速度可能远大于垂直速度。一个保守的估计是：% CFL 稳定性检查 vmax = max(vp0(:)) * sqrt(1 + 2*max(eps(:))); % qP 最大相速度近似 dt_max = 0.6 * min(dx, dz) / vmax; % 安全系数 0.6 if dt > dt_max warning('dt 超过稳定性上限，建议 dt <= %.6f s', dt_max); endsqrt(1+2*eps)是 qP 波在水平方向的相速度与垂直速度的比值近似，ε 越大，水平速度越快，允许的 dt 越小。安全系数取 0.6 而不是 1.0，是给高阶差分和边界反射留余量。我一般会先跑 200 步看波场有没有指数增长，有的话先把 dt 砍半再排查其他原因。3.2 震源主频与网格频散震源主频f0决定了波长，而网格间距dx、dz必须足够小才能压制数值频散。经验规则是每个最短波长至少 8 到 10 个网格点。VTI 里最短波长对应的是 qS 波，不是 qP 波，因为 S 波速度低、波长短。% 网格频散检查 vs_min = min(vs0(:)); lambda_min = vs_min / f0; % 最短波长 ppw = lambda_min / min(dx, dz); % 每波长网格点数 fprintf('每最短波长网格点数 = %.1f\n', ppw); if ppw < 8 warning('网格频散风险高，建议减小 dx/dz 或降低 f0'); end如果ppw小于 8，快照里会看到波前面后面拖着一串波纹，那不是物理频散，是数值频散。解决办法要么加密网格（内存翻倍），要么降低主频（分辨率下降）。这份资源默认参数如果是均匀模型，通常 ppw 在 10 以上，但换成低速层模型时一定要重新算。3.3 边界处理：PML 还是 sponge资源里如果带了吸收边界，常见的是 PML（完美匹配层）或 sponge 阻尼层。PML 效果好但代码复杂，sponge 简单但需要足够厚的阻尼带。如果你在快照里看到边界处有强反射，先确认吸收边界有没有开启，再检查阻尼系数是不是太小。% sponge 阻尼层示例（四周各 20 个网格点） nb = 20; damp = ones(nz, nx); for i = 1:nb damp(i,:) = damp(i,:) .* (i/nb); damp(end-i+1,:) = damp(end-i+1,:) .* (i/nb); damp(:,i) = damp(:,i) .* (i/nb); damp(:,end-i+1) = damp(:,end-i+1) .* (i/nb); end % 在每个时间步后对场量乘 damp阻尼层厚度一般取 10 到 20 个网格点，太薄吸收不干净，太厚浪费计算。PML 的反射系数可以压到 -60 dB 以下，sponge 通常只能到 -30 dB 左右，做定量分析时优先用 PML。4. 避坑与排查：VTI 正演里最容易翻车的五件事4.1 波场快照出现棋盘格振荡现象：快照里波场呈黑白相间的网格状高频振荡，完全看不出波前形态。原因：一阶导数用了普通中心差分但网格没有交错，或者交错网格的场量定义位置搞混了。解决：确认 vx 定义在 x 方向半网格、vz 定义在 z 方向半网格，应力定义在整网格；检查差分算子是否匹配场量位置。如果资源里用的是高阶差分，检查系数表有没有抄错。4.2 C13 出现复数导致波场瞬间发散现象：跑几十步后波场值变成 NaN 或复数，快照全白。原因：Thomsen 参数组合不满足C13^2 > 0的物理约束，根号里出现负数。解决：检查 ε、δ、γ 的取值，δ 不能太小以至于2*delta*C33*(C33-C44) + (C33-C44)^2 < 0。常见做法是先用均匀模型验证，再逐步加分层。4.3 震源加错分量导致只出 P 波或只出 S 波现象：快照里只有一种波在传播，看不到 qP 和 qSV 的耦合形态。原因：震源只加在了szz或只加在了sxz上，且模型各向异性太弱。解决：VTI 正演里通常把震源加在正应力分量上激发 qP，加在剪应力分量上激发 qSV；如果想同时看到两种波，可以在同一位置同时加两个分量，或者用爆炸震源。4.4 快照保存频率太高导致内存爆掉现象：跑一半 MATLAB 报内存不足。原因：每步都存快照，400×400 的 double 矩阵每帧 1.28 MB，1500 步就是 1.9 GB。解决：每隔 50 到 100 步存一帧，或者只存波场切片而不是全波场。我一般会先跑一遍不存快照的，确认稳定后再开保存。4.5 分层模型界面处出现非物理反射现象：两层介质界面处有强反射，但理论上波阻抗差异没那么大。原因：参数场在界面处突变，而差分算子跨界面时用了平均或直接差分，导致等效参数不连续。解决：对C11、C33、C44等刚度系数在界面处做算术平均或调和平均，而不是直接取网格点值。常见做法是把刚度系数定义在整网格，速度定义在半网格，更新时自动取平均。5. 进阶技巧：用波场快照验证各向异性参数是否生效跑通之后，怎么确认你的 VTI 代码真的在算各向异性，而不是各向同性换皮？最直接的办法是对比不同 ε 下的波前形态。固定 δ 和 γ，把 ε 从 0 调到 0.25，看水平方向的波前是不是变快了。如果波前纹丝不动，说明耦合项没生效。% 各向异性参数扫描：对比不同 eps 的波前 eps_list = [0, 0.1, 0.2]; for k = 1:length(eps_list) eps = eps_list(k) * ones(nz, nx); % 重新计算 C11 C13 并跑正演 run_vti_forward; % 假设封装好的正演函数 % 提取 t=0.8s 时的波前水平到达距离 snap = snapshots{8}; [~, ix_max] = max(abs(snap(iz_src, :))); fprintf('eps=%.2f, 水平到达网格 = %d\n', eps_list(k), ix_max); end如果ix_max随 ε 增大而增大，说明各向异性生效。另一个验证方法是看 qSV 波前有没有出现三叉区（triplication），这是 VTI 介质在强各向异性下的典型特征，各向同性代码永远跑不出来。还有一个实用技巧：把波场快照按时间顺序拼成 GIF 或视频，比单帧看更容易发现波前形态的异常。MATLAB 里用imwrite配合循环就能生成 GIF，帧间隔取 0.1 秒左右，太快看不清，太慢文件大。注意：生成 GIF 前先把每帧的caxis统一到全局最大值，否则闪烁严重。从那以后我每次拿到新的各向异性正演代码，都强制先跑一个 ε=0 的退化测试，确认它和各向同性结果一致，再逐步加参数。这一步能省掉后面大量的排查时间。希望帮到你。本文还有配套的精品资源，点击获取