简介逆时偏移RTMMatlab代码包面向地震勘探与成像算法学习者实现叠后/叠前逆时偏移的基本流程采用波场正传和反传互相关成像并支持对简单高角度速度模型进行测试。压缩包共7个文件体积约15KB包含3个m源码脚本、2个mat数据文件、1个md说明文档和1个licence许可证源码与演示数据划分清晰便于直接运行和阅读修改。描述中提到通过求解程函方程获取激发时间成像条件是弹性波逆时偏移的难点而这份代码以较直观的方式展示了RTM核心思路配合拉普拉斯滤波可进一步增强成像对比度适合理解逆时偏移原理并快速上手实操。目前已有678人学习下载可作为地震数据成像、Matlab算法练习及课程设计的参考资源。1. 逆时偏移 MATLAB 代码先搞清楚它在解决什么再谈怎么跑逆时偏移圈内习惯叫 RTM是地震成像里少数能把陡倾角构造和速度剧烈变化都照顾到的算法。拿到这份 MATLAB 代码目标不是跑出一个彩色剖面就完事而是理解它怎样把震源波场正传、把接收波场反传再在零延迟处做互相关成像。适合刚接触逆时偏移的学生、需要处理二维模型验证算法的工程师或者想从正演过渡到成像的从业者。代码包能直接用但里面有不少需要自己动手调的参数否则容易翻车。下面我会按正演核心、成像主循环、常见问题排查的顺序把这套资源讲清楚。2. 逆时偏移的正演核波动方程离散与三个关键参数RTM 的物理核心是正演模拟成像质量很大程度上由正演算得准不准决定。常见做法是用时间域有限差分法在规则网格上离散波动方程原因很简单MATLAB 的矩阵切片操作天然适合做二维波场更新代码短调试直观。选一阶速度-应力方程而不是二阶声波方程是为了方便加 PML 吸收边界这个概念在 2.2 会专门说。2.1 一阶应力-速度方程为什么代码里要拆成两个分量一阶速度-应力方程把二阶声波方程拆成压力的时间导数和质点速度的时间导数两组式子。二维声波情况下速度场又分成水平分量 vx 和垂直分量 vz所以内存里至少维护三个二维矩阵一个压力场 p、两个速度分量场。代码包里的初始化一般是这样% 定义模型空间与网格参数 nx 401; nz 201; % 网格点数横向401点纵向201点 dx 10; dz 10; % 网格间距单位米 nt 2000; % 时间步数视数据长度而定 dt 0.001; % 时间步长单位秒2.2节讲解取值依据 % 纵波速度模型单位 m/s vp 2000 * ones(nz, nx); % 背景速度 vp(100:end, :) 2500; % 第二层速度用于构造两层模型这段代码的逻辑是先铺一个均匀背景速度再把第 100 行以下替换成高速层构成一个简单的两层介质。参数说明里有两处要留意dx 和 dz 的取值不是越小越好常用经验是一个最短波长内至少有 10 个网格点nt 与 dt 的乘积要覆盖地震记录的总长度别只看道数。速度场与压力场的交错更新是正演的主干简化版如下% 交错网格更新注意这里暂时去掉 PML 项方便看清核心 for it 1:nt % 先更新压力场散度项 p(2:end, :) p(2:end, :) ... - vp(2:end, :).^2 .* dt ./ dx .* (vx(2:end, :) - vx(1:end-1, :)); % 再更新水平速度分量 vx(:, 1:end-1) vx(:, 1:end-1) ... - dt ./ dx .* (p(:, 2:end) - p(:, 1:end-1)); % 最后更新垂直速度分量 vz(1:end-1, :) vz(1:end-1, :) ... - dt ./ dz .* (p(2:end, :) - p(1:end-1, :)); end这个更新顺序不是随意写的。压力场依赖速度场的空间导数而速度场更新又依赖压力场交错排列能保证二阶精度这是有限差分里的经典做法。逻辑说明里我一般会强调一点p 和 vx、vz 的索引边界要留出至少一个点的余量否则更新时会访问到矩阵外的位置。参数含义中vp.^2 是介质弹性参数换成密度变化模型时这里要额外除以密度但那份代码包给的是常密度声波假设。2.2 网格步长与时间步长从 CFL 条件到参数表网格步长和时间步长直接决定稳定性这也是新手最容易翻车的地方。常见做法是在给定速度模型的条件下先算最大速度 vmax再用 CFL 条件反推时间步长上限。代码包里已经有现成脚本% 自动计算时间步长上限 vmax max(vp(:)); % 模型中最大速度 cfl 0.6; % CFL 系数一般取 0.6 上下 dt_max cfl * min(dx, dz) / (vmax * sqrt(2)); dt dt_max; % 实际使用值逻辑说明CFL 条件是显式有限差分稳定性的底线大于这个值波场会指数发散表现为成像剖面里出现越来越大的随机噪声。参数上cfl 不建议超过 0.7因为还要给 PML 吸收层留余量dx 和 dz 通常保持一致避免空间各向异性误差。下面这张表是项目里常用的推荐值参数推荐范围说明dx/dz510 m按主频 3040 Hz 与 2500 m/s 估算一个波长至少 10 个点cfl0.50.7大于 0.7 容易出现低幅高速噪声这类噪声很难靠滤波去掉dt由上式自动计算不要手动设一个整数除非你清楚自己在做什么nt记录时长 / dt时间总长不够深部界面就没法成像2.3 PML 边界吸收实现与参数敏感性没有吸收边界正演波场会在模型边缘反弹反向传回来时叠加在真实反射上让成像结果出现一圈高亮弧形。这份代码用的是 PML 边界实现上是在原方程右侧加一个衰减项衰减系数随距离呈抛物线变化。给一段典型的 PML 设置片段% PML 层参数设置 npml 20; % 吸收层厚度单位网格数 pml_a 1.5; % 衰减系数调节因子 sigma pml_a * vmax / (dx * 2); % 理论上限参考值 % 向右边界注入衰减项 for i nx-npml1:nx sigma_x(i) sigma * ((i - (nx-npml)) / npml)^2; end逻辑说明衰减项不是越大越好过大会形成人为反射过小又吸收不干净。经验做法是先按 vmax 估算上限再在 0.4 到 2.0 倍之间扫几个值观察边界处波场幅度是否迅速下降。参数说明里npml20 对应 200 米左右的吸收层对大多数二维模型够用如果模型横向极大pml 厚度可以保持但衰减系数要微调。这段是排查边界伪影的关键后面第 4 章会再讲。3. 用这份代码跑通一次二维 RTM从合成记录到成像剖面这一章是完整的落地流程建议把代码包解开后对照着看。我用一个两层模型做例子目标是输出一张干净的零延迟互相关成像剖面。3.1 准备速度模型与震源子波先给速度模型做平滑。真实代码包里通常会有一个 smooth 函数避免突变界面产生过强的人为绕射。之后生成雷克子波作为震源。% 平滑速度模型 vp0 vp; % 原始阶梯模型 sm 5; % 平滑半径单位网格点 vp imgaussfilt(vp0, sm); % 高斯平滑避免数值绕射 % 生成雷克子波主频 25 Hz dt 0.001; f0 25; t 0:dt:0.2; w (1 - 2 * (pi * f0 * (t - 1/f0)).^2) .* exp(-(pi * f0 * (t - 1/f0)).^2);逻辑说明平滑半径 sm 不宜过大否则会把薄层抹掉雷克子波的起跳时间设为 1/f0是为了保证零时刻子波近似为零减少正演初值引入的伪影。参数说明里主频 f0 与模型规模相关dx10m 时可取 2535 Hz再高就需要把 dx 加密到 5m否则空间采样不足。3.2 主循环一正演传播与炮记录存取正演阶段要为多个炮点重复运行。每炮先在震源位置注入子波然后更新波场把模型边界处的一条或几条记录保存下来作为合成地震记录。代码包里的核心循环是% 单炮正演并保存接收记录 recv zeros(nr, nt); % 接收记录nr 为接收道数 p zeros(nz, nx); vx p; vz p; % 初始化波场 for it 1:nt p(src_z, src_x) p(src_z, src_x) w(it); % 震源加载 % 更新压力场与速度场省略 2.3 的 PML 项 upd_wavefield(); recv(:, it) p(recvz, recvx); % 抽取接收点值 end save(shot_1.mat, recv, vp, dt, dx); % 保存供逆时步使用这里有个细节震源加载不是简单赋值而是在每个时间步累加 w(it)物理含义是把震源看成时间序列的点源。逻辑说明里提醒一句代码包经常会把upd_wavefield()单独的脚本你检查循环时重点看它是否在每个时间步内调用两次漏调会直接导致记录里少一半能量。内存方面我一般会用 single 存储 recv可以省一半内存精度影响不大recv single(recv);参数说明nr 是接收道数对应野外采集时的道间距如果模型横向是 401 点接收道通常间隔 1 道就行别过度抽取。3.3 主循环二逆时外推与零延迟互相关成像逆时偏移最核心的步骤是先把接收记录作为边界或震源项逆时加载同时把震源波场从正演阶段的末端倒播回去在每个时间步对源波场和接收波场做互相关并累加。% 逆时外推与成像条件 image zeros(nz, nx); % 成像累加器 for it nt:-1:1 % 逆时加载接收记录把第 it 时刻的记录反向注入边界 p(recvz, recvx) recv(:, it); % 用同样的差分算子更新波场 upd_wavefield(); % 零延迟互相关成像条件 image image p_src .* p_recv; % 前提同时保存源波场 end但这里会有一个经典问题直接保存每一时刻的源波场 p_src 需要存储 nt 个二维矩阵模型稍大就爆内存。代码包里的常用做法是分两步走——先正演一遍并把每个时间步的 p_src 写盘再在逆时时读回。如果磁盘空间有限就得用更高级的随机边界存储方案但那份代码没做到那一步。参数说明里我建议把 update 函数单独封装正演和逆时调用同一份这样能避免两边差分不一致导致的相位误差。这个坑在 4.3 还会遇到。3.4 成像输出归一化与去噪处理零延迟互相关得到的原始 image 噪声大、能量不均匀直接显示基本没法看。代码包的尾部通常有一段简单后处理% 对成像结果做幅度归一化 image image ./ max(abs(image(:))); % 拉普拉斯滤波器去除低频背景噪声 lap [0 -1 0; -1 4 -1; 0 -1 0]; image_clean imfilter(image, lap, replicate); image_clean image_clean ./ max(abs(image_clean(:)));逻辑说明拉普拉斯滤波本质是取二阶差分能压制成像剖面里横向不变的低频条带同时保留陡倾角构造。参数说明里’replicate’ 是边界填充方式如果换成默认的零填充会在模型上下边界产生两条假亮线。如果深部能量太弱还可以加一个 AGC 窗口窗口长度取主波长的一半代码包里一般都有现成函数这里就不展开。4. 逆时偏移常见问题排查三个反复踩过的坑4.1 成像剖面出现平行横条纹现象剖面里出现一组间距稳定的横纹理跟真实构造走向无关呈全图均匀分布。原因最常见的是震源子波长度没有取完整或者子波起跳时间太短导致正演初值引入直流分量。解决办法检查雷克子波首尾是否都归零。经验上把 t 范围从 0.2 秒扩大一倍再去掉前面的零填充条纹就会消失。如果还在就把第三节里的image image ./ max(abs(image(:)))换成image image ./ max(abs(mean(image,2)))先压横向均值。4.2 模型边界处出现强反射弧现象成像剖面四周有连续高亮往内部衰减像相框一样包住整个模型。原因PML 吸收不彻底边界反射的波场参与互相关加上逆时过程中边界加载记录也引入了旁瓣。解决办法把 npml 从 20 加到 40重新调整吸收系数同时检查 2.3 代码里衰减项是不是只加在了压力场上速度场也要各加一段对称的衰减项。我见过某个工程里漏写了 vx 的衰减项结果边界伪影比真实反射还强调了半天才发现。4.3 内存溢出或计算中途卡死现象程序运行几十秒后弹 out of memory或整体变慢到无法接受。原因波场矩阵默认是 double 精度三个二维矩阵加一个图像累加器本身不算大但保存源波场时用了 double很快就会撑爆内存。解决办法统一在初始化时加 single 转换p single(p); vx single(vx); vz single(vz);此外源波场分块保存也是个办法每次只保存 50 个时间步逆时读回时同步更新对应段的成像结果。这样内存消耗从 nt 个矩阵降到 50 个矩阵代价只是循环多跑几遍。5. 一个进阶技巧角度域成像条件与强制检查主频的习惯原始代码里的零延迟互相关成像条件得到的是每个点的能量叠加深度越深照明越弱尤其构造两侧会出现能量空洞。我常用的一个改进是把它改成简易的角度域成像条件在互相关前对源波场和接收波场分别做希尔伯特变换然后根据相位相关性给强角度贡献打折。代码只需要改一小段% 角度域成像条件简化版 src_h imag(hilbert(p_src)); % 源波场相位旋转90度 rec_h imag(hilbert(p_recv)); % 接收波场相位旋转90度 image image p_src .* p_recv src_h .* rec_h;逻辑说明这个近似把互相关从纯能量叠加变成包含相位加权高角度入射的道对成像贡献被压制等效于照明补偿。实际验证时成像剖面在断层下盘的噪声会明显减少。参数上hilbert 每次调用都计算整场跑大数据时比较耗时我一般只在前几十个时间步用这个条件后面步数走零延迟几乎不影响效果。从那以后我每次跑逆时偏移前都强制检查三件事子波主频是否匹配网格间距、时间步长是否满足 CFL、PML 层厚度是否足够。这三项出了问题后面再调成像条件也救不回来。希望帮到你。本文还有配套的精品资源点击获取