简介：一套面向地震勘探与地球物理研究人员的VTI介质Matlab数值模拟工具包，用于模拟各向异性介质中地震波传播，解决波场正演与可视化问题。VTI介质指纵波速度在垂直与水平方向存在差异的岩石层，这种各向异性普遍存在于实际地壳结构中，会显著影响地震波的传播路径与速度。压缩包共5个文件，含4个.m脚本和1个.mat数据文件，其中.m脚本负责波动方程离散、有限差分求解及完全匹配层吸收边界处理，.mat文件提供自定义色标以渲染波场快照，总大小仅7KB，轻量便携。已有521人学习下载，代码结构清晰，适合具备一定Matlab基础的地球物理专业学生、科研人员或勘探工程师快速上手。运行后可获得不同时刻的波场快照，直观观察地震波在VTI介质中的传播、反射与折射现象，理解吸收边界条件对模拟精度的提升作用，为地震资料解释、地下结构反演与储层评估提供参考。
1. VTI介质数值模拟到底在做什么：一个让波场快照“不圆”的Matlab工程
拿到“VTI numerical stimulation”这个Matlab工程，先别急着解压、跑例程。它背后解决的是一个很具体的地球物理问题：在页岩、薄互层这类垂直横向各向同性（VTI）介质里，各向同性假设会让正演波场的旅行时差几个毫秒，成像点偏移几十米。这个工程的核心产出是数值模拟的波场快照，你能在图上看到qP波波前从圆形变成准椭圆，qSV波出现交叉型的三叉区特征。适合正在做地震正演、逆时偏移调试或岩石物理建模的研究生和工程师，也适合课程设计想找一套能讲清楚原理的数值模拟代码的人。这篇笔记把模型参数、代码骨架和调试过程一次讲透，照着跑完，你至少能亲手输出一张有VTI特征的波场快照。
2. VTI介质的弹性矩阵与Thomsen参数：为什么不能用各向同性代码糊弄
2.1 刚度矩阵里放什么：C11、C13、C33、C44各自的角色
各向同性介质只需要两个弹性参数（拉梅常数λ、μ），因为它的刚度矩阵只有两个独立分量，任意方向波速都一样，波前是标准圆。VTI介质不同：它有水平对称面，垂向是唯一特殊方向，二维x-z平面内的弹性矩阵变成了五个独立分量。对这个标题下的波场模拟来说，更关键的一点是——二维x-z剖面上P-SV波和SH波是解耦的，SH波独立传播，不参与qP和qSV的耦合。这意味着你写Matlab代码时不需要C66，只需要四个分量：C11、C13、C33、C44。
这四个分量的物理意义最好记成“模量”而不是“速度”：C33是垂直方向纵波模量，C11是水平方向纵波模量，C44是垂直方向横波（SV波）模量，C13是最麻烦的耦合项，它把水平应变和垂向应变耦在一起，不直接对应任何单一波速。很多初学者盯着C13发呆，但它正是VTI和非VTI的差别所在：没有C13的耦合作用，波场快照里就不会出现qSV波的典型交叉特征。
VTI刚度矩阵（x-z平面内相关的部分）： | C11 C13 0 | | C13 C33 0 | | 0 0 C44 |2.2 Thomsen参数到刚度系数：一张换算表和一个边界条件
写代码时你不会直接输入C11、C13这种量级在十亿级别的弹性系数，而是用Thomsen参数：垂直P波速度Vp0、垂直S波速度Vs0、密度ρ、各向异性参数ε和δ。ε和δ是VTI介质描述的核心，ε控制水平与垂直P波速度差异，δ控制中等角度（45度附近）的相速度异常。换算公式如下。
| 刚度系数 | 换算公式 | 需要输入的参数 |
|---|---|---|
| C33 | ρ * Vp0² | Vp0、ρ |
| C11 | C33 * (1 + 2ε) | Vp0、ρ、ε |
| C44 | ρ * Vs0² | Vs0、ρ |
| C13 | sqrt((C33 - C44)² + 2δ * C33 * (C33 - C44)) - C44 | Vp0、Vs0、ρ、δ |
这里的C13换算公式是从Thomsen的δ定义反解出来的，有个隐藏坑：根号里面的值必须为正。如果Vp0 > Vs0（正常介质都满足），临界条件约等于δ > -(C33 - C44) / (2C33)。比如Vp0=3000m/s、Vs0=1500m/s时，δ小于-0.375就非法了，程序会在第二十步左右直接炸出NaN。实际地层δ一般处于-0.2到0.3之间，大多数时候安全，但遇到一些特殊夹层还是要检查。
设置参数时，如果你的Matlab变成里用了eps做变量名，会在这一节踩到莫名其妙的坑——eps是Matlab内置的浮点精度函数，直接覆盖它会导致整个会话精度相关的计算全乱。常见做法是用epsilon和delta作为变量名。
2.3 速度-应力方程组：为什么交错网格是这题的默认答案
有了刚度系数，下一步是选控制方程。理论上可以直接对位移二阶波动方程做二阶差分，但VTI里C13会把二阶交叉偏导搅在一起，数值差分对网格位置极其敏感，换个参数就吃不准稳定性。更稳的做法是使用一阶速度-应力方程：
ρ ∂vx/∂t = ∂σxx/∂x + ∂σxz/∂z ρ ∂vz/∂t = ∂σxz/∂x + ∂σzz/∂z ∂σxx/∂t = C11 ∂vx/∂x + C13 ∂vz/∂z ∂σzz/∂t = C13 ∂vx/∂x + C33 ∂vz/∂z ∂σxz/∂t = C44 (∂vx/∂z + ∂vz/∂x)这组方程里没有二阶空间导数，全是交错的一阶偏导，天然适合交错网格有限差分。交错网格的核心思想是让不同的场变量错开半格存放：速度分量落在半网格位置，应力分量落在整网格位置，这样每个一阶导数都可以用两点的中心差分完成，稳定性、频散性质都好很多。下一章的代码就是按这个方案搭的。
3. 用Matlab实现VTI交错网格有限差分：最小可运行骨架与参数表
3.1 交错网格排布：哪个量落在哪个半格
交错网格的实现里，最让新手犯晕的是数组索引和物理位置的对应关系。我一般用统一的数组尺寸，不在Matlab里真的把某个量存成nx+0.5的位置，而是用“差分时错一位索引”来等效半格。具体排布规则如下：vx的物理位置在x方向半格、z方向整格；vz在x方向整格、z方向半格；σxx和σzz在整格；σxz在x和z两个方向都是半格，落在四角点。写差分循环时按这个规则错开索引，效果和真正交错的网格一致。
σxx, σzz 位于 (i, j) vx 位于 (i+1/2, j) vz 位于 (i, j+1/2) σxz 位于 (i+1/2, j+1/2)这个排布表是后面所有循环索引的依据。如果你把σxz放错了位置，比如放在和σxx同一个整格位置，就会出现棋盘状振荡，振幅随步数指数增长，最后发散。后面第5章节会专门说这个现象。
3.2 空间差分系数：二阶起步，四阶进阶
交错网格的空间差分可以用二阶精度，也可以用更高阶。二阶精度实现最简单：对位置x+dx/2处的导数，用x处和x+dx处的两个值相减除以dx。高阶更平滑但代码更长，对于VTI这种各向异性介质，高阶差分能明显减少网格频散。四阶精度交错差分系数是：导数近似等于(-f(x+1.5dx) + 27f(x+0.5dx) - 27f(x-0.5dx) + f(x-1.5dx)) / (24dx)，这套系数用起来比二阶好很多。
模板原则上建议先从二阶跑通流程，确认波场快照形态正确后，再换成四阶。换的时候注意代码里所有导数的计算必须用同一套系数，不能在几个方程里混用二阶和四阶，否则应力更新和速度更新的频散特性不匹配，波前会出现人为的波纹。
3.3 主循环代码：速度更新、应力更新与震源加载
下面这段代码是最小可运行的骨架，模型大小为251×251网格，用二阶交错差分。为了让你看清交错索引的错位规律，我保留了循环写法。Matlab循环慢，但这段代码在小网格下可以在几十秒内跑完，适合先做验证。
% 参数设置 clc; clear; % 介质参数（VTI页岩典型值） Vp0 = 3000; % 垂直P波速度，m/s Vs0 = 1500; % 垂直S波速度，m/s rho = 2200; % 密度，kg/m^3 epsilon = 0.20; % Thomsen epsilon delta = 0.10; % Thomsen delta % Thomsen参数转刚度系数 C33 = rho * Vp0^2; C11 = C33 * (1 + 2*epsilon); C44 = rho * Vs0^2; C13 = sqrt((C33 - C44)^2 + 2*delta*C33*(C33 - C44)) - C44; % 网格参数 nx = 251; nz = 251; % 网格点数 dx = 5; dz = dx; % 空间步长，m dt = 0.0006; % 时间步长，s，满足CFL条件，见3.4 nstep = 2500; % 时间步数 snap_int = 100; % 每100步输出一张波场快照 % 波场数组初始化 vx = zeros(nx, nz); vz = zeros(nx, nz); sxx = zeros(nx, nz); szz = zeros(nx, nz); sxz = zeros(nx, nz); % 震源参数 f0 = 25; % 雷克子波主频，Hz t0 = 1.5 / f0; % 子波峰值延时，s isx = 126; isz = 126; % 震源位置 % 吸收边界衰减带（见4.2节详细解释） L = 40; alpha = 0.006; damp = ones(nx, nz); for k = 1:L w = exp(-alpha * (L - k)^2); damp(k, :) = damp(k, :) * w; damp(nx-k+1, :) = damp(nx-k+1, :) * w; damp(:, k) = damp(:, k) * w; damp(:, nz-k+1) = damp(:, nz-k+1) * w; end % 主循环：交错网格有限差分 for it = 1:nstep t = it * dt; src = (1 - 2*(pi*f0*(t - t0))^2) * exp(-(pi*f0*(t - t0))^2); % 速度更新，注意σxz错位索引 for i = 2:nx-1 for j = 2:nz-1 vx(i,j) = vx(i,j) + dt/rho * ( ... (sxx(i+1,j) - sxx(i,j)) / dx + ... (sxz(i,j) - sxz(i,j-1)) / dz ); vz(i,j) = vz(i,j) + dt/rho * ( ... (sxz(i,j) - sxz(i-1,j)) / dx + ... (szz(i,j+1) - szz(i,j)) / dz ); end end % 应力更新 for i = 2:nx-1 for j = 2:nz-1 dvx_dx = (vx(i,j) - vx(i-1,j)) / dx; dvz_dz = (vz(i,j) - vz(i,j-1)) / dz; sxx(i,j) = sxx(i,j) + dt * (C11*dvx_dx + C13*dvz_dz); szz(i,j) = szz(i,j) + dt * (C13*dvx_dx + C33*dvz_dz); sxz(i,j) = sxz(i,j) + dt * C44 * ( ... (vx(i,j+1) - vx(i,j)) / dz + ... (vz(i+1,j) - vz(i,j)) / dx ); end end % 震源加载：加到法向应力上，近似点膨胀源 sxx(isx,isz) = sxx(isx,isz) + src * C11 * dt; szz(isx,isz) = szz(isx,isz) + src * C33 * dt; % 吸收边界衰减 vx = vx .* damp; vz = vz .* damp; sxx = sxx .* damp; szz = szz .* damp; sxz = sxz .* damp; % 波场快照输出（见4.3节） if mod(it, snap_int) == 0 figure(1); imagesc(vz'); colormap(flipud(gray)); axis image xy; title(sprintf('VTI wavefield, t=%.3f s', it*dt)); drawnow; end end代码里有三个地方需要专门说明交错索引。第一，vx更新用sxz(i,j) - sxz(i,j-1)而不是sxz(i,j+1) - sxz(i,j)，因为vx的物理位置在z方向整格，而sxz在z方向半点，求导必须取整格上下两个半点值才能落在整格。第二，vz更新用sxz(i,j) - sxz(i-1,j)，道理一样，sxz在x方向是半点，取差分后落在vz所在的整格。第三，应力更新时vx和vz的导数也都是错位相减，保证每个导数的计算位置落在被更新变量的物理位置上。这个规律摸透了，改成四阶差分也不容易出错。
3.4 稳定性条件与参数表：先跑小网格再放大
时间步长不是随便取的，交错网格二维弹性波有严格的CFL条件：dt <= dx / (Vmax * sqrt(2))。这里的Vmax取模型中的最大P波速度，在VTI中要按最坏方向算，用sqrt(C11/rho)而不是Vp0。上面代码里sqrt(C11/rho)约等于3556m/s，dx / (Vmax * sqrt(2)) = 0.000994s，折算出安全系数后取dt=0.0006s。如果dt取到0.001s附近，波场会在数百步内发散。
| 参数 | 推荐值 | 说明 |
|---|---|---|
| 网格间距dx | λ_min / 10 | λ_min按最高有效频率2.5*f0计算 |
| 时间步长dt | 0.6 * dx / (Vmax * sqrt(2)) | 再小一点更安全，但会增加计算时间 |
| 雷克子波主频f0 | 10 - 40 Hz | 地震尺度的常用频段 |
| ε、δ | 0.1 - 0.3 | 太小快照看不出各向异性，太大可能出现不稳定 |
| 衰减带宽度L | 30 - 40 格 | 太窄吸收效果差，太宽浪费计算量 |
网格间距的直觉标准是：25Hz主频、Vp0=3000m/s，波长约120m，dx=5m时每个波长有24个网格点，远小于1/10波长，频散可以忽略。如果dx放大到10m甚至15m，快照的高频部分很快出现滞后尾巴。做完小网格验证后再把nx和nz加到501或更高，循环版本会明显变慢，这时可以试试4.4节的矩阵切片优化。
4. 震源、吸收边界与波场快照输出：让第2500步的图还能看
4.1 雷克子波与震源加载方式：三种写法的差别
雷克子波是地震数值模拟的默认震源，公式为(1 - 2(πf0(t-t0))²) exp(-(πf0(t-t0))²)，没有直流分量，频谱是干净的单峰。选f0时要考虑模型尺寸：25Hz对应波长约120m，模型只有1250m宽，勉强能放两三个波长，快照会好看一些。震源加载方式有三种常见写法，效果差异不小。
第一种加载到法向应力上，即sxx += src*C11*dt; szz += src*C33*dt，这是上面代码的方式，近似点膨胀源，在VTI介质中激发的波场以qP波为主，qSV相对较弱。第二种是加载到速度分量vx += src; vz += src，等价于体力源，会在近场引入较强的偶极成分，快照上能明显看到额外的角向辐射花样。第三种是只加载到某一个应力分量上，比如sxx += src，相当于一个水平偶极源，在VTI介质中反倒会激发出更强的耦合qSV波。我的建议是先用第一种跑通流程，因为它的波形干净，后续提取相速度时波前位置容易判读。等你想专门研究qSV波特征时，再去试第三种。
4.2 吸收边界：衰减带参数与“够用不完美”的边界
模型四周如果不处理，波场到达边界后会反射回来，在快照上形成像第二个震源一样的虚假圆弧。3.3节的代码里用的是衰减带方案，也叫乘性阻尼边界：把波场数组在整个模型上乘一个空间衰减系数，让靠近边界的区域的振幅逐步减小。衰减带宽度L取30到40格，系数alpha按0.005到0.008调节。设定时注意用指数函数exp(-alpha*(L-k)^2)，从边界向内部是平滑递增的，不能在边界形成阶跃突变。
必须坦白这个方案的局限：衰减带只是把到达边界的振幅压下去，不是真正的无反射吸收，部分能量会在衰减区内部产生轻微重反射。对小模型、课程设计、参数敏感性分析足够了，但如果你要拿波场做逆时偏移成像，边界反射会污染成像结果，那时候老老实实上PML（完美匹配层），别在衰减带上硬撑。PML实现复杂度至少翻一倍，等第一版跑通了再换不迟。
4.3 波场快照与视频输出：色标、坐标轴和逐帧保存
波场快照是这个标题的核心产物，输出代码值得单独打磨。最基本的imagesc能出图，但要看清VTI的弱振幅特征，要注意三个设置。第一，色标要对称，用红蓝白或黑白对称色标，不要用jet，否则波峰波谷的正负号会被颜色误导。第二，色标范围不能默认铺满最大振幅到最小振幅，建议裁剪到最大振幅的10%到20%，这样qSV波的弱振幅波前面才不会被qP波的大振幅淹没。第三，坐标方向用axis image xy，让z轴向下是正方向，符合地震剖面习惯，免得你从图上读深度时正负搞反。
% 波场快照输出与AVI视频保存 vw = VideoWriter('vti_wavefield.avi'); vw.FrameRate = 20; open(vw); for it = 1:nstep % ... 主循环代码 ... if mod(it, snap_int) == 0 figure(1); imagesc((0:nx-1)*dx, (0:nz-1)*dz, vz'); colormap(flipud(gray)); clim([-1 1] * max(abs(vz(:))) * 0.15); axis image xy; xlabel('x (m)'); ylabel('z (m)'); title(sprintf('VTI wavefield at t = %.3f s', it*dt)); drawnow; writeVideo(vw, getframe(gcf)); end end close(vw);注意clim函数是R2022a之后的名字，老版本要用caxis。如果看到图没反应，先说一句你Matlab版本，再决定改哪一行。快照不仅为了看，还要为后续提取相速度服务：建议把有代表性的快照数据本身用save存成.mat文件，而不只是存图片，不然你要重新跑一遍主循环才能回看某个时刻的波场。
4.4 计算量与优化：循环太慢时改矩阵切片
251×251的网格、2500步，纯循环在Matlab里大概要几分钟，能接受。如果扩到501×501，就要先算一笔账——内层循环次数接近25万，乘以三个更新方程组，再乘2500步，Matlab的循环开销会让单次运行超过半小时。这时我一般改成矩阵切片版本，把内层空间循环一次性向量化。
% 矩阵切片方式替代内层循环（示意） for it = 1:nstep % 速度更新切片写法，边界行最后再处理 vx(2:nx-1, 2:nz-1) = vx(2:nx-1, 2:nz-1) + dt/rho * ( ... (sxx(3:nx, 2:nz-1) - sxx(2:nx-1, 2:nz-1)) / dx + ... (sxz(2:nx-1, 2:nz-1) - sxz(2:nx-1, 1:nz-2)) / dz ); vz(2:nx-1, 2:nz-1) = vz(2:nx-1, 2:nz-1) + dt/rho * ( ... (sxz(2:nx-1, 2:nz-1) - sxz(1:nx-2, 2:nz-1)) / dx + ... (szz(2:nx-1, 3:nz) - szz(2:nx-1, 2:nz-1)) / dz ); % 应力更新同理 % ... end切片写法的风险在于索引错位，每一行都要对照交错网格排布表确认，错一个索引就会在某一步突然冒出NaN。而且边界一行一列在切片里没有被更新，需要单独赋值或依赖衰减带压制，这和循环版本的处理不完全一样。新手如果第一次就跑切片版本，出了bug很难分辨是物理问题还是索引问题，所以我的习惯是先跑通循环版本，确认波场形态正确后，再用切片版本做效率优化。
5. VTI模拟常见问题与排查：五次翻车记录
5.1 波场全是噪点或者Inf：先查CFL和C13正定性
现象：运行到几十步，波场快照突然变成密密麻麻的随机噪点，数值里有NaN或者Inf，整个图是蓝色的死区。
原因：最常见是时间步长超过了CFL限制，导致差分格式数值不稳定，振幅随步数指数爆炸。另一个隐藏原因是C13换算时根号内为负，刚度矩阵不满足正定性，物理上就不该存在这个介质，程序自然翻车。
解决：先检查dt <= 0.6 * dx / (sqrt(C11/rho) * sqrt(2))，不满足就调小dt。再打印C13的值，如果是NaN或者虚数，回去检查Vp0、Vs0、δ的搭配是否符合2.2节的边界条件。这两个因素占了这个标题下80%的发散问题。
5.2 快照里只有P波没有SV波：震源加载方式的锅
现象：波场快照只有一个最外侧的准圆形波前，内部看不到任何qSV波的交叉特征，整张图单调得像各向同性介质。
原因：震源加载在法向应力上且只激发纯膨胀源时，qSV波在某些VTI参数下被激发得很弱，特别是δ较小时，qSV波的振幅本身就比qP波低一个量级。再加上色标范围自动铺满到qP波的最大振幅，qSV波那条弱波前死死地藏在色标阈值以下，肉眼看就是消失。
解决：把clim范围手动设置到最大振幅的10%到15%，或者改用3.1节说的第三类加载方式，把震源只加到σxx分量上，人为制造偶极成分，qSV波的激发效率会明显提升。如果你观察的是固定方向，还可以把vz换成vx分量做快照，不同分量对不同波型的敏感度不一样。
5.3 波前锯齿加高频振荡：网格频散与差分阶数
现象：波场快照的大波前后面拖着一条条细碎的尾部波纹，像锯齿一样，离震源越远越明显。
原因：空间网格相对于最高频率波长太粗。雷克子波虽然主频是f0，但实际频谱延伸到2.5f0附近，这段高频成分在粗网格上传播速度比低频慢，产生网格频散。
解决：把dx缩小到最高有效频率对应波长的1/10以下。比如f0=25Hz时，按2.5f0算最高有效频率约62.5Hz，对应波长48m，dx至少小于4.8m，取5m正好在边界上。预算允许就取dx=2.5m，锯齿基本消失。另一个解法是换四阶差分，它能把频散误差再往下压一到两个量级，但代码改动面广，先试dx。
5.4 边界反射像第二个源：衰减带不是PML
现象：波到达边界后明显折返，快照内部出现一圈向震源方向汇聚的圆弧，振幅甚至和原始波前差不多强。
原因：衰减带宽度太窄或系数太小。特别是衰减区内部的幅值梯度不够平滑，波场经过时发生阻抗突变，一部分能量反射回内部。
解决：加宽L到40格以上，alpha调到0.006附近，并检查衰减剖面的平滑性。衰减带的设计初衷是省事，不要为了省内存把宽度压到20格以下。如果发现加宽了还有强反射，如实说这是衰减带的极限了，这个场景该上PML，PML的衰减函数通常取exp(-2 * (L-k) / L)配合吸收系数，和这里的思路完全不同，别混用。
5.5 图看起来“圆”得可疑：eps、色标裁剪与坐标方向
现象：代码跑得很顺利，快照也干净，但波前是标准圆形，VTI的椭圆特征一点都没有，或者波前有椭圆特征但方位角方向像翻转了。
原因：圆得可疑先查epsilon和delta是不是真的起作用了——检查代码里变量名是否叫eps，Matlab里eps是内置函数，你把它当变量名赋了值，可能在某个运算中被悄悄当成了机器精度常数。各向异性参数如果为0，波前自然变圆。翻转问题则是坐标轴方向：imagesc默认y轴从上到下，如果你不写axis xy，z轴正方向会朝上，VTI的椭圆长轴方向看着就反了。
解决：变量统一用epsilon，参数打印出来核对。坐标轴强制axis image xy。还有最后一个检查项：色标裁剪过狠会把弱各向异性特征完全切掉，ε=0.1、δ=0.05这种小参数本来就很难肉眼分辨，建议先把ε调到0.3、δ调到0.15跑一张确认趋势，再回调到目标值。
6. 从波场快照验证代码对不对：提取qP相速度并与解析解对比
光看波场快照“像”VTI介质还不够，你得有办法量化验证代码没有系统性错误。常用的验证方案是：从快照中提取某个方向的视速度，与Thomsen弱各向异性近似公式对比，差别在5%以内基本说明方程、参数、差分都正确。
具体步骤是：在某一帧快照上，从震源位置出发沿角度θ（以垂直方向z为0°）画一条测线，用interp2插值提取沿线振幅序列，找到最大振幅的位置，距离除以时间就是该方向的群速度近似值。然后用弱各向异性公式Vp(θ) = Vp0 * (1 + δ*sin²θ*cos²θ + ε*sin⁴θ)计算理论相速度，两个速度对比即可。
% 从波场快照提取45度方向视速度 theta = deg2rad(45); Lmax = 1000; % 测线最大长度，m s = linspace(0, Lmax, 2000); xp = x0 + s * sin(theta); zp = z0 + s * cos(theta); % z方向，theta从垂直方向量起 vals = interp2((0:nx-1)*dx, (0:nz-1)*dz, vz', xp, zp, 'cubic'); [~, imax] = max(abs(vals)); V_measured = s(imax) / (it_snap * dt); V_theory = Vp0 * (1 + delta*sin(theta)^2*cos(theta)^2 + epsilon*sin(theta)^4);这里有一个理论细节值得知道：你从快照上拾取的是群速度（能量传播速度），而公式给出的是相速度（波前法向速度），两者在强各向异性介质中角度差可以达到好几度。不过当ε、δ小于0.2时，这个差异在1%到3%以内，用来做代码正确性的粗验证完全够用。实测偏差如果超过10%，别急着怀疑理论公式，先回第5章检查色标裁剪和震源加载方式，这两个是最容易引入系统性读图误差的地方。
另外推荐一个长久受益的小习惯：每次修改参数，都跑同一方向、同一时刻的快照做差分对比。比如想确认ε从0.2改到0.3的效果，用两次数值模拟的快照相减，剩余波场就是各向异性变化造成的响应，比盯着两张彩图猜差异可靠得多。我自己有次把δ设成了负值，快照看着还挺像模像样，就是45度方向的波前速度比垂直方向还慢，对理论公式一查才发现介质参数早就物理非法了。那次之后我再也不敢只听信快照的“形状直觉”，每跑完一组参数，都要做一次定量验证。这套流程希望你也能用起来，希望帮到你。
本文还有配套的精品资源，点击获取
