简介针对地震逆时偏移RTM原理学习与代码实现需求这份Matlab代码以简洁的工程化方式提供了可运行的叠前/叠后逆时偏移示例。代码以简单速度模型为对象通过波场正传与反传的互相关完成成像较适合高角度构造偏移试验配合说明文档可帮助理解叠后爆炸反射面成像、叠前激发时间成像条件以及程函方程求取单程旅行时等核心难点。资源共7个文件主要为3个M函数正演模拟、反演模拟与主程序RTM、2个MAT数据文件速度模型与地震记录、1份Markdown说明及许可证压缩包仅15KB轻量易用便于直接运行和修改试验参数。目前已有678人学习下载适合地球物理专业学生、地震数据处理初学者或需要快速搭建RTM测试环境的科研人员也可作为后续开展拉普拉斯滤波、改善高角度成像效果的参考代码。1. 逆时偏移是什么为什么复杂构造成像绕不开它做过逆时偏移RTM的人都知道速度模型一旦出现横向突变或构造翻转克希霍夫偏移和单向波动方程偏移的成像结果往往直接翻车剖面缺块、回转波归错位、角度构造被抹平。RTM的思路是把源波场从震源处正传、检波波场从记录处反传两个波场在同一时刻逐点互相关把所有能被波场到达的反射点归位到正确深度。它的实现基础并不玄学一个双程波动方程时间步进加一个互相关累加就够了。用MATLAB写一个可运行的最小RTM半天时间足够且非常适合做速度模型验证、成像条件对比和论文公式的落地复现。适合三类人处理解释人员验证速度模型、研究生复现文献算法、想摸清成像底层的开发者。2. 从原理到代码RTM的四个核心模块RTM不是黑匣子。把这四个模块拆出来你会发现它和普通的声波正演模拟只差一个反传成像的尾巴波场外推、成像条件、存储策略、整体流程。逐个说清楚后面改代码就不慌。2.1 波场传播双程波动方程的有限差分离散RTM用的是双程声波波动方程∂²p / ∂t² v² · (∂²p / ∂x² ∂²p / ∂z²)这个方程描述压力场在介质中的完整传播既包含下行波也包含上行波所以能自动处理回转波和棱柱波。离散时最常用的是时间二阶、空间二阶中心差分格式p(t1) 2·p(t) - p(t-1) v²·dt²·∇²p(t)这个格式简单、易实现、占内存小是MATLAB原型和教学代码的首选。空间二阶差分在MATLAB里用矩阵切片就能写好不需要循环我有很长时间都在用下面这个函数function lap laplacian2d(p, dx, dz) % 计算压力场的拉普拉斯算子空间二阶中心差分 % 输入 p: nz×nx 的波场矩阵 % 输出 lap: 与 p 同尺寸的拉普拉斯场 [nz, nx] size(p); lap zeros(nz, nx); % x 方向二阶差分 lap(:, 2:end-1) (p(:, 3:end) - 2*p(:, 2:end-1) p(:, 1:end-2)) / dx^2; % z 方向二阶差分 lap(2:end-1, :) lap(2:end-1, :) ... (p(3:end, :) - 2*p(2:end-1, :) p(1:end-2, :)) / dz^2; % 边界不参与差分由边界条件负责 end逻辑说明这个函数每次返回一个完整的拉普拉斯场时间循环里直接调用即可。dx 和 dz 必须带单位米而且要一致否则差分算子会把各向异性误差带进波场。边界行和边界列不计算差分置零靠PML吸收边界吸收出射波这是常规做法。参数说明p 是双精度矩阵nz 是行数深度、nx 是列数横向。dx 和 dz 的单位通常取 5~20 米。如果模型是 400×1000 网格这个函数每步要执行约 5 毫秒对单炮 RTM 的总耗时影响不大瓶颈主要在快照存储和反传循环次数。2.2 成像条件零延迟互相关为什么能成图成像条件的物理直觉是如果在某个深度存在反射界面那么当源波场到达该界面时会激发一个反射波这个反射波到达地表被记录在逆时偏移中我们把记录到的波场反传回地下当反传波场到达那个界面的同一时刻源波场也恰好到达同一个位置。两个波场在该点相遇、相位一致乘积为正且最大。对每个网格点沿着整个时间轴累加就得到该点的成像值I(x,z) Σ_t S(x,z,t) · R(x,z,t)S 是源波场R 是反传的检波波场。这个零延迟互相关最大的优点是稳健不需要事先知道反射系数和入射角对所有构造类型统一处理。但它有两个明显毛病。第一浅部照明强、深部照明弱直接互相关的成像振幅随深度衰减明显第二整个成像结果里会带着宽低频背景界面轮廓被淹掉一半。针对第一个问题工业界普遍改用归一化成像条件或激发振幅成像条件本质是用源波场能量做除法把照明不均匀校正掉针对第二个问题最通用的补救是拉普拉斯滤波把低频背景当作空间二阶导数滤除掉。这两招我在第4章会给出具体代码和参数。2.3 存储与策略源波场快照的三种通道反传过程中检波波场每个时刻都要和同时刻的源波场做互相关但源波场是正传时的时间顺序相反。所以要么先把源波场存下来要么在反传时重新算一遍。三种策略我都在工程里用过存储策略内存/磁盘开销适用场景主要缺点全快照nz×nx×nt×8 字节double小模型、教学、算法验证大模型直接爆内存检查点每 K 步存一帧源波场中等模型反传时需重算中间步耗时增加随机边界重建不存源波场快照大模型、生产级需要设计随机边界实现复杂全快照是最容易实现的方式。举例一个 120×80 网格、700 步的模型double 类型快照需要 80×120×700×8 ≈ 54MBMATLAB 直接存没问题。但如果换成一个 500×1000 网格、3000 步的模型快照就是 500×1000×3000×8 12GB内存直接告急。真到了那一步先把 double 改成 single内存减半还不够就用检查点策略。2.4 RTM代码的整体流程与最小框架把上面三块拼接起来单炮 RTM 的流程如下读入速度模型、震源子波和炮记录正传源波场逐时刻保存快照或做检查点把炮记录作为边界条件做反传从最后一个时刻往前推每推进一步就把当前源波场快照和检波波场相乘累加最后做一次拉普拉斯滤波输出。我用伪代码归纳一下% RTM 单炮最小流程伪代码 % 1 正传源波场 for it 1:nt pf wave_step(pf, pf_old, vp, dx, dz, dt, pml); snap(:, :, it) pf; % 保存源波场快照 end % 2 反传检波波场 成像 for it nt:-1:1 pr wave_step(pr, pr_old, vp, dx, dz, dt, pml); pr inject_record(pr, rec_data(it, :), rec_z); % 在表层注入炮记录 image image snap(:, :, it) .* pr; % 零延迟互相关成像 end % 3 拉普拉斯滤波去低频 image laplacian2d(image, dx, dz);逻辑说明正传和反传可以用同一个波场步进函数差别只在时间方向和加载项。正传时加载的是震源子波反传时加载的是炮记录序列。把这一步抽象成函数后替换PML、成像条件、速度模型都只改一处不会牵一发动全身。这个模块化结构是整套代码好维护的关键。3. 最快的落地路径改写一份可跑的逆时偏移MATLAB代码纸上谈兵结束直接给一份能落地的单炮 RTM 最小实现。我会把代码拆成三段模型准备、正传和炮记录合成、反传成像。你把三段存成同一目录下的脚本再把 2.1 的 laplacian2d 函数保存成 laplacian2d.m 就能跑。这里用合成数据验证不需要真实采集数据。3.1 准备速度模型与震源子波先定义网格参数、速度模型和雷克子波。我用一个简单的高速透镜体作为目标体正传产生的反射波将来自透镜边界速度界面清晰方便判断成像深度对不对。% RTM_demo.m —— 单炮逆时偏移最小实现 clear; clc; % 网格与时间参数 nx 120; nz 80; % 横向/纵向网格数 dx 10; dz 10; % 空间采样单位米 dt 0.001; % 时间采样单位秒 nt 700; % 时间步数 % 速度模型背景 2200 m/s嵌入一个高速透镜体 vp ones(nz, nx) * 2200; xg ((1:nx) - 0.5) * dx; % 网格中心横向坐标 zg ((1:nz) - 0.5) * dz; % 网格中心深度坐标 [XX, ZZ] meshgrid(xg, zg); lens (XX - 600).^2 / 180^2 (ZZ - 400).^2 / 120^2 1; vp(lens) 3400; % 雷克子波延迟一个主频周期近似零相位 f0 30; % 主频 30 Hz t (0:nt-1) * dt; tw 1 / f0; src (1 - 2 * (pi * f0 * (t - tw)).^2) .* exp(-(pi * f0 * (t - tw)).^2);逻辑说明速度模型用矩阵方式直接构造透镜区域内速度从 2200 突增到 3400 m/s这是我们要成像的反射界面。雷克子波是 RTM 最常用的震源子波延迟 tw1/f0 这一行非常关键它让子波从接近零的幅度开始增长避免波场在时间轴上出现半周期错位。若省略延迟成像深度会系统性偏移。参数说明dx、dz 取 10 米模型横向 1.2km、深度 800m属于小模型。dt0.001 秒需要满足稳定性条件后文会专门讲怎么算。f030Hz 对应的波长在背景速度下约 73 米每个波长约 7.3 个网格满足频散控制要求。速度模型里的透镜中心在横向 600m、深度 400m 处对应网格坐标 60、40 附近。3.2 正传源波场并同时合成炮记录正传这一步做的事把震源子波加载到源点按双程波动方程向前推 700 步每一步保存源波场快照同时把表层一整排网格点的波场记录成炮记录。这里我加了一个简化 PML 吸收边界防止边界反射在成像阶段捣乱。% —— 正传源波场保存快照并在浅层合成炮记录 —— npml 15; % PML 吸收层厚度单位网格 pml ones(nz, nx); pml_w 0.002; % 衰减强度越大吸收越强 for ii 1:npml att 1 - pml_w * (npml - ii 1) / npml; pml(ii, :) pml(ii, :) .* att; pml(end-ii1, :) pml(end-ii1, :) .* att; pml(:, ii) pml(:, ii) .* att; pml(:, end-ii1) pml(:, end-ii1) .* att; end % 震源和接收位置放在吸收层内侧避免记录被 PML 衰减 sz npml 1; rec_z npml 1; rec_x 1:nx; % 全排列接收 % 波场与存储变量 pf zeros(nz, nx); pf_old pf; snap zeros(nz, nx, nt); % 源波场快照 rec_data zeros(nt, nx); % 炮记录 for it 1:nt lap laplacian2d(pf, dx, dz); pf_new 2 * pf - pf_old vp.^2 .* lap * dt^2; pf_new pf_new .* pml; % 边界吸收 pf_new(sz, floor(nx/2)) pf_new(sz, floor(nx/2)) src(it) * dt^2; pf_old pf; pf pf_new; snap(:, :, it) pf; rec_data(it, :) pf(rec_z, rec_x); % 记录浅层波场作为后续反传输入 end逻辑说明时间循环内部的顺序是固定的先算拉普拉斯再做时间递推再乘吸收系数最后加震源。震源加载用了 src(it)*dt^2这是压力源在速度-压力格式下的标准离散方式。snap 和rec_data 都是在每个时间步完成后立即赋值保证快照与记录时间对齐。这里 rec_data 是从同一个正传过程里合成出来的炮记录真实资料处理时直接替换成野外采集记录就行。参数说明npml 取 15 是中等厚度对 80 行深度的模型来说约占 19% 的网格偏保守但安全。pml_w 决定衰减梯度强度0.002 在这个网格尺寸下效果可以吸收层内部波场衰减到接近零。震源和接收刻意放在 npml1 行因为如果放在更浅的位置记录会包含吸收层的衰减痕迹反传时会产生人为的振幅不连续。3.3 反传检波波场与成像输出反传从最后一个时刻开始把炮记录从浅层注入计算域向过去方向外推。每推进一步用当前检波波场乘以同时刻的源波场快照并累加得到互相关成像剖面。最后用拉普拉斯滤波压掉低频背景再显示。% —— 反向传播检波波场零延迟互相关成像 —— pr zeros(nz, nx); pr_old pr; image zeros(nz, nx); for it nt:-1:1 lap laplacian2d(pr, dx, dz); pr_new 2 * pr - pr_old vp.^2 .* lap * dt^2; pr_new pr_new .* pml; % 在浅层注入当前时刻的炮记录 pr_new(rec_z, rec_x) pr_new(rec_z, rec_x) rec_data(it, :) * dt^2; pr_old pr; pr pr_new; % 零延迟互相关成像条件 image image snap(:, :, it) .* pr; end % 拉普拉斯滤波去低频背景噪声 image_f laplacian2d(image, dx, dz); figure; imagesc(xg, zg, image_f); axis xy; colormap(gray); xlabel(Horizontal (m)); ylabel(Depth (m)); title(RTM image);逻辑说明反传循环的方向是反的但步进公式和正传一模一样。关键区别在注入项正传加载的是震源反传加载的是炮记录序列且逐时刻往回收。成像累加放在循环体尾部用第 it 时刻的检波波场去乘第 it 时刻的源波场快照这一步的时间对齐必须严格否则成像深度会乱。参数说明反传用的速度模型、PML 系数必须和正传完全一致任何一处改了都要重新正传才能保持一致性。如果换成真实炮记录要记得先对记录做直达波切除或时窗衰减否则剖面浅部会出现强水平条带我在第5章具体讲这个坑。4. 参数怎么调从成像效果倒推四个关键参数代码跑通只是第一步成像质量不理想时先别怀疑算法按顺序检查 dx/dt、PML、子波主频、成像条件四个参数。我把它们列成速查表每一列对应一种常见病。参数作用常规取值调大影响调小影响dx / dz空间采样每波长 5~10 格频散更小、耗时增加频散明显、成像模糊dt时间采样≤ dx / (v_max·√2)数值不稳定、直接 NaN计算时间变长npml吸收层厚度10~30 网格边界反射更少但耗时增加弧形伪影出现f0子波主频与目标分辨率和 v 匹配分辨率高、穿透变浅界面模糊、穿透加深4.1 dt 与 dx 的稳定性约束和频散控制时间步长不是随便取的。对时间二阶、空间二阶的声波方程稳定性条件是 dt ≤ dx / (v_max·√2)。举例v_max3400 m/s、dx10m 时dt 上限约 2.08ms我取 1ms 留了余量。如果 dx 加密到 5mdt 上限降到 1.04ms得取 0.5ms 量级。违反的后果很直接波场在计算中指数增长几百步内变 NaN图像一片雪花。频散和空间步长绑定。一个波长内网格数太少不同频率成分的传播速度不一致波前面会出现龟裂状震荡。我一般保证每个波长至少 5 个网格稳妥取 8~10 个。用一个关系式倒推dx ≤ v_min / f0 / 8其中 v_min 是最小速度、f0 是子波主频。若 f030Hz、v_min2200m/sdx 应不大于 9m10m 勉强及格但也能跑。4.2 PML吸收边界的两个关键参数PML 参数出问题成像剖面最外圈会出现一圈向内的弧形伪影看起来像括弧。两个参数起决定性作用吸收层厚度 npml 和理论反射系数 R。业界通用公式是 α_max -3·v·log(R) / (2·L)L 是吸收层物理厚度npml·dxR 取 0.001~0.0001。R 越小衰减越强但衰减曲线过于陡峭会在吸收层内部产生数值反射反而引入新噪声。我的经验npml 少于 10 基本没用多于 30 对计算量拖累明显。中等模型取 15~20 最合适。如果发现最外圈伪影先把 npml 加到 20~25再调 R不要一上来就猛增衰减强度。MATLAB 教学代码里常见的乘衰减系数方案是简化版边界反射能压到-20dB左右够学习但不够生产。真的要压到-40dB以下建议换 C-PML 或分裂场 PML实现会复杂一些但排错逻辑不变。4.3 子波主频与相位主频 f0 决定分辨率和穿透深度。f0 太高波长太短深部能量衰减快而且更容易频散f0 太低界面成像模糊陡倾角构造的分辨率跟不上。调 f0 时记住一个换算f030Hz、v2500m/s 对应波长约 83 米能分辨约 20 米厚的薄层如果你要分辨的目标是 50 米以上的构造完全够用。比主频更容易踩坑的是子波相位。我在 3.1 里特意让子波延迟一个主频周期 tw1/f0就是为了让子波零相位化。如果省略这一步成像界面会整体偏移约半个波长换算成深度约 v·(1/f0)/4这个系统误差很难从成像图上看出来但和已知模型一对比就露馅。真实资料里如果不知道子波的准确相位常规做法是用自相关或井数据估计零相位子波。4.4 成像条件与照明补偿零延迟互相关的天生问题是浅部能量强、深部能量弱因为浅部源波场和检波波场都强深部都弱。要让深部构造显形需要做照明补偿。最常用的改进是激发振幅成像条件% 在成像累加处对源波场能量做归一化 image image snap(:, :, it) .* pr ./ (snap(:, :, it).^2 eps);逻辑说明多了一次源波场能量除法浅部因为能量大被自动压权重深部因为能量小被抬权重。eps 是防止除零的稳定项建议取值比波场振幅平方小几个数量级比如波场振幅约 1e-2 时eps 取 1e-6 合适。注意这个除法会放大噪声所以它在使用时要和拉普拉斯滤波配成组合只做除法不做滤波剖面会非常脏。参数说明eps 太小会除零太大会把补偿效果吃掉。更稳妥的做法是先把能量图单独存下来成像累加完成后统一除一次这样能直观看到照明分布我一般用后者。5. 逆时偏移MATLAB落地避坑五个常见问题和排查思路跑 RTM 的过程里我踩过的坑大多不是算法推导错误而是参数和边界处理不当。下面按现象 → 原因 → 解决写五个高频问题照着排查能省大半天时间。5.1 现象成像剖面有横向连续的水平条带浅部尤其明显浅部出现连续强振幅水平同相轴且和地下界面位置无关通常是直达波在作祟。合成炮记录里直达波能量比反射波强一个数量级反传时它会作为强边界源被注入计算域和源波场做互相关后形成宽水平条带。解决分两步第一步在反传前对炮记录做直达波切除或时窗衰减时窗起点按直达波走时 t dist/v_near 来定第二步保留拉普拉斯滤波把残余低频压掉。做完这两步浅部条带会明显干净。5.2 现象剖面最外围出现一圈弧形伪影像括弧外圈括弧几乎都是边界反射。简化 PML 的衰减强度不足或吸收层太薄时出射波到达计算边界后被反射回计算域与源波场形成伪相关。排查方法造一个完全均匀的速度模型不做任何反射体只放一个震源正传看边界处的波场幅值。如果边界处还能看到明显回传能量就是 PML 参数问题。先把 npml 从 15 加到 25再调大衰减强度通常能压制到背景水平。5.3 现象深部完全没有成像或只在浅部有一团模糊深部没图像多数不是算法问题而是记录长度或接收孔径不够。深部反射波需要足够长的记录时间才能到达接收点计算一下目标深度 z_max 对应的双程走时 t_max ≈ 2·z_max / v_avg把 nt 取到 t_max/dt 的 1.2 倍以上。接收排列如果只在震源附近一小段深部反射能量会被偏移出排列覆盖范围。对策是扩大 rec_x 的接收范围或者多炮叠加。这一步在做合成数据验证时可以先量好理论值排错会快得多。5.4 现象程序运行到一半内存耗尽或者卡死内存爆炸是 RTM 在 MATLAB 里最常见的翻车点。全快照的占用是 nz×nx×nt×8 字节double 类型下 200×400 网格、2000 步就有 1.28GB已经肉眼可见地卡了。第一招是把所有波场数组改成 single内存直接减半第二招是改用检查点策略每 K 步存一帧源波场反传时从最近检查点重新正传给后续时刻。检查点会显著增加计算时间但这是 MATLAB 里做大模型的必经之路没有捷径。5.5 现象成像界面位置和已知模型明显对不上深一度或浅一度成像深度系统性偏移优先怀疑子波延迟和速度模型尺度。子波没延迟时源波场整个时间轴提前约 1/f0互相关成像的深度也相应偏移这个我在 3.1 中已经用 tw 处理了。速度模型如果整体偏快或偏慢成像深度会按比例错位。验证方法是单独造一个水平界面模型80 网格深、速度 2000m/s正演合成记录后做 RTM看成像界面是否落回 80 网格。能通过这个测试说明代码正确问题在目标模型本身。6. 进阶用MATLAB并行与GPU让RTM从单炮走向多炮6.1 parfor 多炮并行最省事的提速手段多炮 RTM 的每炮完全独立天然适合并行。在 Parallel Computing Toolbox 的环境下把单炮 RTM 包成一个函数用 parfor 分发给多个 worker 即可。我通常把速度模型作为只读变量传入每炮只更改炮记录和震源位置% 多炮并行每炮独立成像最终叠加 images zeros(nz, nx, nshot); parfor shot 1:nshot % 每炮读取自己的炮记录速度模型 vp 只读共享 image_shot run_single_rtm(vp, rec_data_set(:, :, shot), ... src, dx, dz, dt, nt, pml); images(:, :, shot) image_shot; end final_image sum(images, 3) / nshot; % 多炮叠加提升信噪比逻辑说明parfor 里每炮独立持有 wavefield 数组不会互相污染。叠加是 RTM 生产中提高信噪比的标准手段单炮成像噪声大、多炮叠加后构造轮廓清晰稳定。如果机器是多核 CPU这个改动能把 16 炮的总耗时从串行的 16 倍压到接近 nshot/核数 的比例。显存充足时还能把核心函数改成 gpuArray 版本主要瓶颈是每炮的快照存储用 single 可以缓解。参数说明nshot 不能开太大同时并行太多 worker 会让内存峰值叠加。我一般让 parfor 并行数不超过物理核数的一半剩下留给系统和其他 MATLAB 窗口。多炮叠加时各炮成像振幅可能有量级差异最好先做归一化再叠加否则能量大的炮会主导最终剖面。6.2 验证RTM代码正确性的一个技巧水平界面单炮检验我现在的习惯是任何新模型上手前先造一个最简单的地层模型——两层速度中间一个水平界面界面深度放在网格 40 附近。正演合成记录再跑一遍 RTM看成像深度是否精确落在界面位置、振幅极性是否正确。这一步通过说明波场外推、成像条件、快照存储三个环节自洽如果连这个简单模型都错位那就没必要去跑复杂构造模型先回头查子波延迟和 PML 参数。这个最小验证的价值在于把代码问题和模型问题分开。复杂模型成像不好可能是速度模型不准、采集孔径不足、子波相位未知这些和 RTM 代码本身没关系。而水平界面测试只检验代码链路所以我会把它当作每次修改代码后的回归用例。做 RTM 这些年我最大的体会是别怕慢先在简单模型上把每一步想明白再上复杂数据才是真正省时间的路子。希望这些踩坑记录能帮你把 RTM 的 MATLAB 原型跑得更顺也希望你在调参时多留存几组对比剖面它们比任何文档都更能说明问题。本文还有配套的精品资源点击获取
