简介：面向地震波场模拟与高性能计算方向的研究者与工程师，这份资源聚焦二维TTI各向异性介质下的有限差分正演与逆时偏移实现，并借助CUDA将计算密集的波场传播与成像过程迁移到GPU并行加速，适合已具备波动方程与有限差分基础、希望深入理解各向异性RTM算法落地细节的中高级学习者。压缩包共8个文件，约18KB，以4个C源文件、3个CUDA核函数文件及1个Makefile组成，涵盖速度模型读取、拉普拉斯算子、弹性波一阶方程正演、逆时偏移与ADCIGs提取等模块，编译脚本便于直接构建运行。目前已有474人学习下载。读者可从中获得一套结构紧凑、可编译的TTI介质正演与RTM参考实现，对照代码理解各向异性波场延拓、GPU并行策略与角度道集生成的衔接方式，为算法复现与性能优化提供可借鉴的工程骨架。

1. 当各向异性遇上 GPU：TTI 正演与逆时偏移到底难在哪

很多人第一次接触 TTI 介质成像，会以为它只是把各向同性波动方程里的速度换成一个带方向的张量，改几行公式就能跑。真正上手才发现，翻车的往往不是公式，而是算力与访存：TTI 有限差分正演每一步要更新应力、应变多个分量，逆时偏移还要把整个时间步的波场存下来再倒着读，单卡显存和带宽瞬间见底。这套基于 CUDA 的 TTI 有限差分正演与逆时偏移实现，解决的正是“公式能写出来、但跑不动”的断层——它把差分算子、边界吸收、波场存储和成像条件全部搬到 GPU 上，用 ADCIGS（角度域共成像点道集）把偏移结果拆成随入射角变化的道集，让各向异性参数对成像的影响第一次变得可观测、可验证。搞懂它，等于搞懂地震成像里“精度换算力”这笔账是怎么在底层被算平的。

2. TTI 波动方程：为什么不能直接套用各向同性那套

2.1 从刚度张量到一阶速度—应力方程

各向同性介质里，波速是一个标量，P 波和 S 波各走各的；TTI（倾斜横向各向同性）介质里，速度变成一个带对称轴的张量，纵波和横波在传播中会耦合，波前面不再是球面而是椭球甚至更复杂的形状。工程上很少直接解二阶位移方程，因为二阶导数的差分模板更宽、边界处理更麻烦，常见的实现思路是把它降阶成一阶速度—应力方程组：速度分量对时间的一阶导，等于应力张量的空间梯度除以密度；应力分量对时间的一阶导，等于刚度张量与速度梯度的乘积。这样每个时间步只需要做一阶空间差分，模板窄、并行度高，天然适合 GPU 的线程级并行。

关键在于刚度张量的旋转。TTI 的对称轴是倾斜的，需要先把本构关系从局部坐标系旋转到全局坐标系，旋转角通常用倾角和方位角描述。旋转后的刚度矩阵不再是简单的几个常数，而是随空间位置变化的系数场。这意味着差分计算里每一次应力更新都要读取一组随网格点变化的弹性参数，访存压力比各向同性大得多。我一般会这么理解：各向同性正演是“一个速度场打天下”，TTI 正演是“每个网格点都有一套自己的物理规则”，这也是后面所有性能优化的出发点。

2.2 有限差分模板的阶数与稳定性约束

空间差分用几阶，直接决定精度和算力。二阶模板每维只用相邻两个点，误差大但访存少；高阶模板（比如八阶）每维要读四个邻居，波场模拟的频散被压得很低，但每个网格点的计算量和访存量成倍上升。TTI 方程分量多，如果每个分量都用高阶模板，寄存器压力和共享内存占用会迅速逼近上限。常见的取舍是：对波场主分量用较高阶模板保证精度，对辅助分量适当降阶，或者用交错网格把速度和应力错开半个网格存放，用同样的模板阶数换到更高的有效精度。

稳定性是另一道坎。显式时间推进的稳定性条件（CFL 条件）要求时间步长小于波速与网格间距决定的阈值，TTI 介质里最大波速由弹性参数共同决定，倾斜对称轴还会让某些方向的等效速度更高。时间步取大了，波场会指数发散，屏幕上就是一片噪声；取小了，迭代次数暴涨，GPU 利用率反而上不去。实际调参时通常先按最坏情况估一个保守步长，再根据频散误差微调，这个“先稳后准”的顺序不能反。

2.3 边界吸收：PML 为什么比简单衰减更可靠

有限计算区域必然有边界，波打到边界会反射回来，污染整个波场。最简单的做法是在边界附近乘一个衰减系数，但衰减太弱反射压不住，太强又会在边界内侧产生虚假反射。PML（完美匹配层）的思路是在边界外扩一层特殊介质，让波进入后不反射、快速衰减，理论上对任意入射角都能吸收。实现上通常把空间导数做复坐标拉伸，在差分代码里体现为对波场分量额外乘一个衰减因子并累加辅助变量。

PML 的代价是每个边界网格点要多存几个辅助波场，显存占用上升，而且拉伸参数选得不好会在层内产生不稳定。血泪经验是：PML 层数不能太薄，通常十几个网格起步；衰减剖面的最大值要和主频、网格间距匹配，拍脑袋设一个常数往往在长时模拟后爆掉。这套实现把 PML 参数做成可配置项，正是为了让使用者能在精度和稳定性之间反复试。

3. CUDA 映射：差分算子怎么落到线程和共享内存上

3.1 三维网格到线程块的映射策略

GPU 上的有限差分，核心问题是怎么把三维网格切成线程块。最直观的映射是每个线程负责一个网格点，线程块取一个立方体区域（比如 8×8×8 或 16×16×4）。块内线程读邻居时，如果邻居落在块边界外，就得走全局内存，延迟高、带宽浪费。所以常见做法是先把当前块需要的波场区域（含 halo 区）整体搬进共享内存，块内所有差分读取都命中共享内存，只有 halo 区从全局内存加载一次。

这里有个容易被忽略的细节：TTI 方程有多个波场分量，每个分量都要做差分，如果每个分量单独开一次 kernel，共享内存的 halo 数据就要重复加载多次。更高效的做法是把相关分量打包进同一个 kernel，一次加载、多次复用。代价是寄存器压力上升，occupancy 下降，需要在“少读几次全局内存”和“多驻留几个线程块”之间找平衡。我一般会先看 profiler 里的共享内存 bank 冲突和全局内存吞吐，再决定块尺寸，而不是一上来就抄一个“看起来很美”的配置。

3.2 共享内存分块与 halo 交换的代码骨架

下面这段是典型的共享内存分块骨架，用伪 CUDA C 写出核心逻辑，重点看 halo 加载和同步的位置：

// 块尺寸 BLOCK_X × BLOCK_Y × BLOCK_Z，halo 半径为 R（由差分阶数决定） __global__ void tti_stress_update(float *vx, float *vy, float *vz, float *sxx, float *syy, float *szz, float *sxy, float *sxz, float *syz, const float *c11, const float *c13, /*...*/ int nx, int ny, int nz) { // 共享内存 tile 含 halo：每维多 2R 个点 __shared__ float tile_vx[(BLOCK_X + 2*R) * (BLOCK_Y + 2*R) * (BLOCK_Z + 2*R)]; // ... 其余分量同理 int tx = threadIdx.x, ty = threadIdx.y, tz = threadIdx.z; int gx = blockIdx.x * BLOCK_X + tx - R; int gy = blockIdx.y * BLOCK_Y + ty - R; int gz = blockIdx.z * BLOCK_Z + tz - R; // 每个线程负责加载 tile 中一个或多个点，含 halo 区 for (int idx = threadIdx.x + threadIdx.y * blockDim.x + threadIdx.z * blockDim.x * blockDim.y; idx < (BLOCK_X + 2*R) * (BLOCK_Y + 2*R) * (BLOCK_Z + 2*R); idx += blockDim.x * blockDim.y * blockDim.z) { // 将 idx 解码为 tile 内坐标，再映射到全局坐标做边界判断 // 越界点填零或做镜像，避免读脏数据 tile_vx[idx] = load_or_zero(vx, gx, gy, gz, nx, ny, nz); } __syncthreads(); // 必须等 halo 全部就位 // 用共享内存做高阶差分，更新应力分量 // 每个线程只算自己负责的中心点，邻居从 tile 里取 float dvx_dx = fd8(tile_vx, tx, ty, tz, /*方向 x*/); // ... 组合刚度系数，写回全局内存 }

这段代码里，R由差分阶数决定，八阶模板R=4，tile 每维要多出 8 个点，共享内存占用随块尺寸立方增长，块开太大直接编译失败或 occupancy 掉到个位数。__syncthreads()的位置不能省也不能提前，halo 没加载完就做差分，读到的是上一轮的旧值，波场会出现块状伪影——这种 bug 在单卡小模型上偶尔看不出来，放大网格后才暴露，排查起来非常费劲。边界判断里“填零或镜像”的选择也有讲究，填零相当于在边界引入一个强反射面，必须配合 PML 才能压住。

3.3 多分量更新的寄存器与访存权衡

TTI 一阶方程里，速度三个分量、应力六个分量，每个时间步都要更新。如果全部塞进一个 kernel，寄存器用量很容易超过每线程上限，编译器开始 spill 到本地内存，性能断崖式下跌。常见的拆分方式是：速度更新一个 kernel，应力更新一个 kernel，交替执行。这样每个 kernel 的寄存器压力可控，但波场要在全局内存里多读写一轮。

另一种思路是把不参与当前计算的中间量放进共享内存而不是寄存器，用共享内存换寄存器。这招在分量多、模板窄的场景下有效，但共享内存本身也是稀缺资源，块尺寸就得相应缩小。实际调优时我习惯先固定一个中等块尺寸，用 profiler 看寄存器数、共享内存占用和 occupancy 三个指标，再决定往哪个方向挪。没有放之四海皆准的配置，只有针对当前方程形式和硬件代际反复试出来的那一组。

4. 逆时偏移与 ADCIGS：波场存储和成像条件的底层账

4.1 逆时偏移为什么要“存波场再倒放”

逆时偏移的基本图像是：震源波场从时间零正推到最大时刻，检波点波场从最大时刻反推到零，两者在同一时刻、同一位置相遇时做互相关，就得到成像值。问题在于，互相关要求两个波场在同一时刻对齐，而检波点波场是倒着算的，正推的震源波场必须“等”它。最直接的办法是把正推过程中每个时间步的震源波场全部存下来，倒推时按时间逆序读出来。这就是逆时偏移最吃显存的地方：波场大小乘以时间步数，三维模型下轻松到几十 GB。

常见的缓解手段是 checkpoint 技术：只存若干关键时刻的波场快照，倒推时从最近的 checkpoint 重新正推一小段，用计算换显存。checkpoint 间隔越小，重算量越大但显存越省；间隔越大则相反。这套实现把存储策略做成可调参数，本质上是在显存和算力之间画一条可移动的线。踩坑点在于 checkpoint 的重算必须和原始正推在数值上完全一致，浮点累加顺序稍有不同，重算出的波场就和存档对不上，互相关成像会出现细微但系统性的偏差。

4.2 互相关成像条件与角度道集的生成逻辑

普通逆时偏移的成像条件是震源波场与检波点波场的零延迟互相关，得到的是叠加成像。ADCIGS 要更进一步：它不只记录“这里有没有反射”，还要记录“反射来自哪个入射角”。实现上通常借助波场的方向信息——比如用波场的空间梯度估计传播方向，或者用 Poynting 矢量判断能量流向，再把震源侧和检波点侧的方向合成入射角，把互相关值累加到对应的角度 bin 里。

角度分 bin 的粒度是个取舍。bin 太粗，角度信息糊成一团，各向异性参数的影响看不出来；bin 太细，每个 bin 分到的能量少，道集信噪比低，还得额外做平滑。常见的做法是分几十个 bin，配合角度域的加权和插值。这里有个容易翻车的地方：方向估计在波场复杂区域（比如绕射点、强横向变速区）本身就不稳，Poynting 矢量在振幅接近零的地方方向会乱跳，直接分 bin 会引入噪声。稳妥的做法是加一个振幅阈值，低于阈值的点不参与角度累加，或者对方向做局部平滑。

4.3 显存与带宽：checkpoint 与波场压缩的取舍

把账算细一点：一个三维波场分量在单精度下占 4 字节，九个分量、512×512×512 网格就是约 4.8 GB 一个时间步，存几百步根本不可能。所以实际系统几乎一定要做分量裁剪（只存参与互相关的必要分量）和空间降采样（checkpoint 存粗网格，重算时插值回细网格）。降采样会引入高频损失，对成像分辨率有影响，但对角度道集的宏观趋势影响有限，这是很多实现愿意接受的近似。

带宽方面，逆时偏移的访存模式是“写一遍、倒着读一遍”，全局内存吞吐是瓶颈之一。用纹理内存或只读缓存加载存档波场，能利用空间局部性减少重复读取；把互相关累加放在共享内存里做局部归约，再写回全局，也能省不少带宽。这些优化单独看都是常规手段，组合起来才决定一套逆时偏移在单卡上能不能跑完一个像样的模型。

5. 避坑与排查：TTI 正演偏移里最容易翻车的几件事

5.1 波场发散：先查 CFL 再查 PML

现象是迭代若干步后波场数值指数增长，画面变成一片高值噪声。原因通常有两个：时间步长超过 CFL 稳定上限，或者 PML 衰减参数过强导致层内不稳定。解决顺序是先按最大等效速度重算 CFL 上限，把时间步降下来试跑；如果仍发散，再检查 PML 层数和衰减剖面，把最大衰减系数调低、层数加厚。不要一上来就怀疑差分公式写错，稳定性问题占了这类故障的大半。

5.2 块状伪影：共享内存同步和 halo 越界

现象是成像结果上出现规则的方块状条纹，周期和线程块尺寸一致。原因多半是 halo 区加载不完整，或者__syncthreads()位置不对，导致部分线程读到未初始化的共享内存。解决方法是逐块打印 halo 加载的全局坐标范围，确认边界点被正确填零或镜像；再检查同步点是否在所有写共享内存的操作之后。这类 bug 在小网格上可能被边界条件掩盖，务必用能整除块尺寸和非整除两种网格各测一遍。

5.3 角度道集发虚：方向估计在低振幅区失效

现象是 ADCIGS 在某些深度段角度分布弥散、看不出清晰趋势。原因是 Poynting 矢量或梯度方向在振幅接近零的区域数值不稳定，噪声被分进了各个角度 bin。解决办法是加振幅阈值门限，低能量点不参与角度累加；或者对方向场做空间平滑后再分 bin。阈值不能设太高，否则有效的小角度反射也被滤掉，需要在道集上对比几组阈值再定。

5.4 重算不一致：checkpoint 的浮点累加顺序

现象是开启 checkpoint 后成像结果与全存波场版本有系统性偏差，且偏差随重算段数增加。原因是重算时的浮点运算顺序与原始正推不完全一致，累加误差被互相关放大。解决办法是保证重算路径和正推路径走同一份 kernel、同一套归约顺序，避免为了省事在重算时换用不同实现。浮点一致性在逆时偏移里不是洁癖，是正确性要求。

5.5 显存溢出：先裁分量再降 checkpoint 密度

现象是模型稍大就报显存不足。原因往往是九个分量全存、checkpoint 间隔又设得很密。解决顺序是先确认哪些分量真正参与互相关，裁掉冗余分量；再逐步拉大 checkpoint 间隔，观察成像质量下降是否可接受；最后才考虑空间降采样。顺序反了，容易在精度已经受损的情况下还在纠结显存，白走弯路。

6. 精度、算力与显存的三方博弈：这套实现的设计取舍

把前面几章串起来看，TTI 有限差分正演与逆时偏移在 GPU 上的每一个设计选择，本质都是精度、算力、显存三者之间的博弈，而且这三者不是线性关系。差分阶数提高，精度上去了，算力和访存同时上涨；checkpoint 间隔拉大，显存省了，重算算力翻倍；角度 bin 分细，ADCIGS 分辨率提高，每个 bin 的能量和信噪比下降。没有哪一项可以单独调到最优，只能针对目标模型和硬件代际找一组“够用且稳定”的组合。

验证一套实现是否可信，我习惯做三件事：用均匀 TTI 模型对比解析解或半解析解的波前形状，确认各向异性参数被正确响应；用简单层状模型看 ADCIGS 的角度趋势是否与理论入射角一致；再关掉 checkpoint、关掉降采样跑一个小模型，作为“金标准”去校准优化版本的偏差。这三步走完，才敢把参数往大模型上推。

从那以后我每次调这类正演偏移的 kernel，都强制先跑一遍小网格的稳定性测试和 checkpoint 一致性对比，再动大模型——省下的不是调试时间，是重新算一遍的算力。希望帮到你。

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