简介基于CUDA的二维TTI介质有限差分正演与逆时偏移实现是一份面向地震勘探与高性能计算开发者的代码资源。内容涵盖TTI具有垂直对称轴的横向各向同性介质中的波动方程离散、GPU并行加速策略、逆时偏移RTM成像以及ADCIGs角度域共成像点道集输出等关键技术有助于理解各向异性介质中地震波传播模拟与成像算法的工程落地。压缩包共8个文件以C源程序、CUDA内核文件和Makefile构建脚本为主整体仅18KB结构紧凑。已有470人浏览学习适合具备一定地震成像与CUDA编程基础的读者参考其并行实现思路用于算法验证或二次开发。 做地震勘探数据成像的人应该都听过TTI这一串字母。TTI介质、有限差分正演、逆时偏移RTM组合在一起是当前各向异性复杂构造成像绕不开的一套组合拳。很多盆地的页岩、裂缝型储层或者盐下构造波场传播路径和速度都会随方向变化如果硬按各向同性处理走时和振幅都会算偏成像位置自然也对不上。这个项目要做的就是基于CUDA把TTI介质里的有限差分波场模拟和逆时偏移整条链路跑通让正演和偏移都能在GPU上加速。适合正在做地震波场数值模拟的研究生、处理复杂构造数据的地球物理工程师以及想认真把科学计算程序迁移到CUDA后端的开发者。如果你是刚开始接触这块我的建议是先把“为什么非要用TTIRTM”想清楚。各向同性RTM只能处理速度随空间变化的简单介质遇到TTI介质波前面的形态和传播方向都变了正演模拟的震源波场不准偏移成像自然也会散。TTI介质里的有限差分正演就是把倾斜对称轴带来的方向各向异性通过参数体放进差分方程里再结合RTM双向波场延拓和互相关成像条件把复杂构造的反射界面“照”出来。整个过程计算量非常大所以GPU并行几乎是必然选择。1. 项目整体设计为什么TTI正演和RTM要一起做1.1 场景需求各向异性模型下的双程波成像在实际地震资料里页岩地层往往表现出很强的VTI特性也就是对称轴垂直的横向各向同性但当构造倾斜、地层褶皱之后对称轴通常不再垂直而是随空间变化这时候就要用TTI来描述。TTI介质的核心参数不只有速度还要有Thomsen各向异性参数epsilon、delta以及对称轴的倾角和方位角。正演模拟时这些参数一起参与波动方程求解所以计算模板比各向同性情况要复杂得多。RTM用的是双程波方程理论上没有倾角限制可以处理回转波、棱柱波这类强复杂波场。但双程波意味着震源波场要沿时间正向传播、检波点波场要沿时间反向传播再把同一时刻的波场做互相关成像。这里最大的拦路虎是存储和计算三维模型动辄上亿网格点一个炮集要跑几千个时间步震源波场的每一时刻都可能需要参与成像。所以你会发现正演和RTM天然绑定在一起正演不光是用来合成地震记录更是RTM内部最核心的计算模块。这个项目把二者放在一条流水线里而不是当成两个独立程序正好省掉了大量重复建模和数据传递工作。1.2 算法选型有限差分与逆时偏移的组合逻辑为什么正演用有限差分不用有限元或者谱方法我的理解是有限差分赢在简单、规则、好并行。TTI介质虽然比各向同性多了一些交叉偏导数项但差分模板仍然是局部操作每个网格点的更新只依赖周围有限个邻居点。这种局部性特别适合GPU每个线程算一个点数据访问的访存模式相对规整性能容易做高。逆时偏移的存储压力也是算法选型时不能忽略的因素。如果每步快照都写盘I/O会直接把GPU省下来的时间吃掉。我用的是有限差分正演生成波场快照配合checkpointing做中间状态保存在反向延拓时再逐步重构波场。这样既保留了RTM需要的全时刻波场信息又不需要把整个四维波场塞进显存或者硬盘。整体思路可以概括为正演负责产生可靠的波场RTM负责把波场转化为成像结果GPU并行负责让这个组合在可接受的时间内跑完。1.3 CUDA并行策略怎么铺开CUDA端的并行策略我采用的是最直接也最稳定的“一个线程负责一个网格点”。空间维度用二维或三维block组织每个线程根据全局坐标找到自己对应的模型点读取当前时刻和前一时刻的波场计算下一时刻的波场并写回。时间步循环放在host端一个时间步调一次内核内核对整个网格做一次更新。这里有两个关键点。第一空间差分模板是局部操作但每个线程都要访问周围多个邻居点如果全从全局内存读访存开销会非常大。我习惯用shared memory做tile缓存先把当前block需要的子区域连同halo边界一起加载到共享内存再计算空间导数这比反复访问全局内存快得多。第二波场用双缓冲切换当前时刻波场和下一时刻波场交替使用避免在同一时刻覆盖数据。这个方案看着朴素但调试起来非常舒服后面加PML边界、加TTI参数体、加RTM成像条件时出问题都能快速定位。2. 核心算法细节与CUDA落地要点2.1 TTI伪声学方程与差分模板的选择TTI实现里最常用的是伪声学qP波方程不是完整的弹性波方程。完整弹性波方程当然更接近物理真实但计算量太大且自由表面边界和横波伪影处理都很复杂。伪声学方程通过VTI模型旋转坐标得到强行让波场以qP波为主在大多数地震成像场景里精度足够。差分格式上我选的是二次时间精度、八阶空间精度。时间二阶是因为时间步长在CFL条件下通常取得比较小降低时间误差比提升空间精度更划算空间高阶则是为了保证波数采样充足减少数值频散。八阶空间模板意味着每个网格点要访问左右各四个邻居点加上TTI方程里的交叉偏导项模板系数比各向同性多不少。也有人推荐Abbott-Lonescu六点隐式有限差分格式时间精度和稳定性指标都更漂亮但隐式格式意味着每步要解大型线性方程组GPU上要额外配稀疏求解器工程复杂度不是一般地高。我建议大部分场景先用显式格式把流程跑通再评估是否值得上隐式。TTI伪声学方程还有一个特别容易踩的坑epsilon或者delta参数不合理时会出现非物理的伪SV波能量甚至直接数值发散。所以在模型准备阶段我一般会先对各向异性参数做平滑和约束保证系数矩阵在正定区间内。速度模型也一样强横向突变的参数体在差分计算里容易激发数值噪声光滑过渡的模型反而能得到更稳定的成像结果。2.2 核函数设计、共享内存和显存布局核函数设计上常见做法是把当前波场、前一时刻波场、下一时刻波场三个数组都放在全局内存里每次调用更新内核后交换指针。为了减少全局内存访问我会为每个block开一块shared memory tile。以八阶空间精度的二维计算为例block内部计算区域可以取16×16再额外加载上下左右各4个点的halo区域也就是共享内存实际尺寸为24×24。加载时注意边界判断内部点直接从全局内存批量拷贝边界点单独处理。共享内存有个容易被忽视的问题bank conflict。二维数组按行存储时同一行的不同列会让相邻线程访问同一bank的不同地址轻则降低访存吞吐重则让我们辛苦优化的性能打回原形。我在共享内存行尾加一个元素的padding让行首地址错开就能明显减少bank冲突。代码结构大致是这样__global__ void wave_update_kernel( const float* cur, const float* prev, float* next, const float* vp2, const float* theta, int nx, int nz, float dt2, float inv_dz2) { int ix blockIdx.x * blockDim.x threadIdx.x; int iz blockIdx.y * blockDim.y threadIdx.y; if (ix nx || iz nz) return; int idx iz * nx ix; // 这里以二阶空间精度示意实际代码是八阶模板 float lap -4.f * cur[idx] cur[idx - 1] cur[idx 1] cur[idx - nz] cur[idx nz]; next[idx] 2.f * cur[idx] - prev[idx] dt2 * vp2[idx] * lap * inv_dz2; }这段代码只是最小示意TTI完整实现时还要把epsilon、delta、倾角旋转项都加进laplacian计算里但内核的组织方式是一样的。显存布局上我坚持用float精度避免double精度带来的显存翻倍和带宽翻倍。对于一个2D模型几个波场数组加参数体已经能轻松占掉几个GB3D模型更应该精打细算。2.3 逆时偏移的存储与成像条件RTM与正演最大的区别在于它需要“同时”使用正向和反向传播的波场。最朴素的做法是把每个时刻的震源波场都存下来反向延拓时直接读但这对显存或者磁盘的要求高到不现实。我采用的是checkpointing加边界重构的组合方案每隔N步保存一次完整波场快照同时保存这N步之间的吸收边界数据。反向时先从最近的checkpoint恢复波场再用保存的边界数据重新正演把中间时刻的震源波场重建出来。成像条件我用的是零延迟互相关。这个条件实现简单在每个时间步把正向震源波场和反向检波点波场相乘并累加最后得到成像剖面。但直接互相关得到的剖面往往存在浅层强振幅、深层弱照明的问题最好在成像累加时用震源波场的能量做归一化也就是做照明补偿。成像后的低频噪声也几乎不可避免我会再做空间方向的拉普拉斯滤波以及自动化增益控制让深层弱反射信息不至于被掩盖。这些后处理看起来只是锦上添花实际上决定了一张剖面能不能直接用。3. 从零搭建工程环境、代码与执行流程3.1 开发环境与CUDA多版本切换开发环境我常用Ubuntu 22.04搭配NVIDIA驱动和CUDA Toolkit。用nvidia-smi可以看到驱动版本用nvcc --version看的是CUDA编译器的版本这两个经常不是同一个数字不要一见不一致就以为装错了。驱动版本决定能支持的最高CUDA runtime而项目里可以通过多版本CUDA并存来兼容不同依赖。很多人在装CUDA时踩过的坑是电脑里已经有一个PyTorch编译好的CUDA版本又想装另一个CUDA Toolkit结果因为覆盖软链导致原来能跑的深度学习代码炸了。我现在的习惯是安装时保持/usr/local/cuda作为软链实际目录用/usr/local/cuda-11.8、/usr/local/cuda-12.3这样区分。需要切版本时只改PATH和LD_LIBRARY_PATH或者直接用update-alternatives管理。另一个常见报错是CUDA error: no kernel image is available for execution这种八成是GPU架构和编译目标不匹配。RTX 3060算力是8.6RTX 4060 Ti是8.9H系列是9.0编译时不能用单个-archsm_86到处跑最好用-gencode同时生成多个架构的cubin或者让编译器嵌入PTX用JIT在运行时做兼容。项目结构我会分得很干净tti_rtm/ ├── include/ │ ├── model.h │ └── rtm_common.h ├── src/ │ ├── main.cpp │ ├── io.cpp │ ├── tti_kernel.cu │ └── pml.cu ├── CMakeLists.txt └── params.jsonCMake里要显式设置CUDA架构列表不要完全依赖CMAKE_CUDA_ARCHITECTURES的默认值同时把NDEBUG在release模式下打开否则边界检查和断言会拖慢内核。3.2 检查设备与CUDA环境是否正常正式开始前先确认机器能看到GPU。最简单的命令是nvidia-smi能显示驱动、显存和进程信息。接着编译并运行CUDA自带的deviceQuery例子确认你的代码能用与GPU匹配的compute capability。找不到cuda samples也是常见问题安装时选择完整toolkit才会带上samples。如果你不想重装也可以写一个两三行的cudaGetDeviceProperties小程序自己查效果一样。如果机器上同时存在多个CUDA版本我在bashrc里只保留一个默认版本的环境变量不把多个路径同时放进LD_LIBRARY_PATH。因为不同版本的libcudart混进同一个进程轻则报版本不匹配重则cudaFree直接crash。使用PyTorch的时候用torch.version.cuda查看它对应的CUDA版本如果RuntimeError: CUDA error: cublas status execution failed大概率是显存不够或者PyTorch与驱动/编译架构不匹配。这里有个实用技巧先用torch.cuda.is_available()和torch.zeros(1).cuda()做最小验证能少猜很多问题。3.3 正演到RTM的完整工程流水线整个工程流水线可以拆成六个步骤。第一步是建模和参数平滑读入速度场、epsilon、delta、倾角、方位角生成GPU端参数体。第二步是波场初始化分配三个波场数组和PML边界数组加载震源子波。第三步是正演主循环每个时间步调用一次波场更新内核同时把每N步的checkpoint和吸收边界数据写盘。第四步是反向重建震源波场从最后一个checkpoint开始用保存的边界数据把中间时刻波场逐步恢复出来。第五步是反向延拓检波点波场并在同一时刻计算互相关成像累积。第六步是后处理对成像结果做拉普拉斯滤波、照明归一化和振幅增益。以我的测试模型为例二维网格取1200×800点空间步长10米时间步长0.8毫秒一共跑3000个时间步。在RTX 3060上一个正演大约一到两分钟量级同样计算放到纯CPU串行程序上起码是小时级。这里数值和具体实现、编译器优化都有关系但差距足以说明GPU方案的价值。RTM因为要做反向重建和互相关耗时大约是正演的三到四倍整体仍然在可接受范围内。4. 踩坑清单CUDA环境、显存与成像质量排查4.1 环境类问题速查我把实际工程里最常遇到的问题整理成一张表排查时先对症状再动手。现象可能原因排查动作运行时报no kernel image available编译架构与GPU算力不匹配用nvidia-smi查算力重新编译时指定对应arch或保留PTXcuBLAS execution failed显存不足或PyTorch/CUDA环境混用减少batch检查LD_LIBRARY_PATH用最小例子验证程序启动就退出但无输出驱动与toolkit版本不匹配对比nvidia-smi和nvcc --version必要时升级驱动deviceQuery不到GPU权限问题或容器未透传显卡加--gpus all确认用户有权限访问设备多版本CUDA切换后编译错乱头文件和库文件混用了不同版本清理缓存只保留一个版本的环境变量这里我想多说一句不要在系统层面反复卸载重装CUDA。切换版本用目录隔离和环境变量组合远比“装一个卸载一个”安全。很多深度学习框架会捆绑自己的CUDA runtime那个装在你自己的conda环境里和系统toolkit并不冲突除非你手动把libcudart.so随便软链到系统目录。4.2 计算与成像类问题排查波场模拟最容易遇到的是NaN发散。排查顺序我固定为先检查时间步长是否满足CFL条件再检查速度模型里有没有异常大值最后检查epsilon和delta参数是否在稳定区间。如果只是局部发散把模型参数平滑一遍基本能解决。PML吸收边界参数不当时边界区域也会出现缓慢增长这时候把PML厚度从20个网格点增加到40个衰减系数相应调大一般能压住反射。RTM剖面出现横轴、条带或者浅层异常亮斑本质上不是GPU的问题而是成像条件噪声。低频噪声用拉普拉斯滤波可以压制但要注意滤波窗口不要太大否则深层有效信号也没了。照明不均匀的问题用震源能量归一化效果很明显。还有一个容易被忽视的点如果速度模型存在强横向突变差分模板在突变处会产生人为散射建议在正演前把模型做空间平滑但不是无脑平均而是保持界面位置大致不变。成像剖面看起来“糊”先别急着换算法。检查一下波场输出精度是不是被float限制得厉害2D小模型可以用double对比一次能直观看到精度损失在哪。很多情况下是checkpoint间隔取得太大中间波场重建精度不够导致互相关能量分散这时减小checkpoint间隔比换更高阶差分格式更立竿见影。5. 实测体会先2D再3D先稳定再提速5.1 性能分析与渐进式验证我把最常见的劝告再说一遍不要一上来就啃3D TTI RTM。先做2D用最简单的水平层状模型正演结果和解析解对比确认波场基本正确再逐步加倾角、加各向异性参数、加复杂构造最后才考虑3D。每加一个模块单独验证一次。性能分析我用NVIDIA Nsight Systems看整体时间分布用Nsight Compute看具体kernel的访存和计算效率。优化顺序是先消除明显浪费比如不必要的cudaMemcpy、没用的同步、过度复杂的if判断再调block尺寸和共享内存大小。实测下来block从64提升到256后性能可能翻一倍但继续加到1024反而可能因为寄存器溢出和调度压力变慢。每个模型的最优block规模不一定相同值得花一晚上做一组参数扫描。5.2 最后一点个人经验真正让这套程序跑得顺的不是我写了多花哨的CUDA代码而是数据结构足够规整、每个内核足够小、每步结果足够可验证。数据访问连续、参数体放在显存里只读、主循环里不掺杂任何临时分支这三条做到位性能自然差不到哪里去。如果要说一句最值得记住的体会我会说先让整个流程在CPU上小规模可复现再上GPU优化。这个顺序能省掉你最多的调试时间。波场模拟涉及大量中间数据你只有知道正确结果长什么样才能判断GPU并行到底对不对。等把这套流程吃透了再看那些高阶有限差分格式、更复杂的各向异性系数、更大规模的三维偏移思路都会清晰很多。本文还有配套的精品资源点击获取