简介：针对地震逆时偏移（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(t+1) = 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 存储与策略：源波场快照的三种通道

反传过程中，检波波场每个时刻都要和同时刻的源波场做互相关，但源波场是正传时的，时间顺序相反。所以要么先把源波场存下来，要么在反传时重新算一遍。三种策略我都在工程里用过：

全快照是最容易实现的方式。举例：一个 120×80 网格、700 步的模型，double 类型快照需要 80×120×700×8 ≈ 54MB，MATLAB 直接存没问题。但如果换成一个 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 最常用的震源子波，延迟 tw=1/f0 这一行非常关键，它让子波从接近零的幅度开始增长，避免波场在时间轴上出现半周期错位。若省略延迟，成像深度会系统性偏移。

参数说明：dx、dz 取 10 米，模型横向 1.2km、深度 800m，属于小模型。dt=0.001 秒需要满足稳定性条件，后文会专门讲怎么算。f0=30Hz 对应的波长在背景速度下约 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-ii+1, :) = pml(end-ii+1, :) .* att; pml(:, ii) = pml(:, ii) .* att; pml(:, end-ii+1) = pml(:, end-ii+1) .* 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 在这个网格尺寸下效果可以，吸收层内部波场衰减到接近零。震源和接收刻意放在 npml+1 行，因为如果放在更浅的位置，记录会包含吸收层的衰减痕迹，反传时会产生人为的振幅不连续。

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、子波主频、成像条件四个参数。我把它们列成速查表，每一列对应一种常见病。

4.1 dt 与 dx 的稳定性约束和频散控制

时间步长不是随便取的。对时间二阶、空间二阶的声波方程，稳定性条件是 dt ≤ dx / (v_max·√2)。举例：v_max=3400 m/s、dx=10m 时，dt 上限约 2.08ms，我取 1ms 留了余量。如果 dx 加密到 5m，dt 上限降到 1.04ms，得取 0.5ms 量级。违反的后果很直接：波场在计算中指数增长，几百步内变 NaN，图像一片雪花。

频散和空间步长绑定。一个波长内网格数太少，不同频率成分的传播速度不一致，波前面会出现"龟裂"状震荡。我一般保证每个波长至少 5 个网格，稳妥取 8~10 个。用一个关系式倒推：dx ≤ v_min / f0 / 8，其中 v_min 是最小速度、f0 是子波主频。若 f0=30Hz、v_min=2200m/s，dx 应不大于 9m，10m 勉强及格，但也能跑。

4.2 PML吸收边界的两个关键参数

PML 参数出问题，成像剖面最外圈会出现一圈向内的弧形伪影，看起来像括弧。两个参数起决定性作用：吸收层厚度 npml 和理论反射系数 R。业界通用公式是 α_max = -3·v·log(R) / (2·L)，L 是吸收层物理厚度（npml·dx），R 取 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 时记住一个换算：f0=30Hz、v=2500m/s 对应波长约 83 米，能分辨约 20 米厚的薄层；如果你要分辨的目标是 50 米以上的构造，完全够用。

比主频更容易踩坑的是子波相位。我在 3.1 里特意让子波延迟一个主频周期 tw=1/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 原型跑得更顺，也希望你在调参时多留存几组对比剖面，它们比任何文档都更能说明问题。

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