企业官网建设流程全解析

简介：一套基于 MATLAB 实现的逆时偏移（RTM）示例代码，适合地球物理、地震勘探方向的初学者或研究者。逆时偏移是地震偏移成像中的常用技术，涵盖叠后与叠前成像思路；实现中采用激发时间成像条件，通过程函方程求取旅行时，完成波场正传与反传的互相关成像。整套代码包含 7 个文件，其中有 3 个 MATLAB 脚本（主程序、正演建模函数、反演建模函数）、2 个数据文件（速度模型与单炮记录）、1 份说明文档和 1 份许可证文件，整体压缩包仅 15KB，结构简洁、轻量易部署。当前已有 678 人学习，适合快速上手逆时偏移流程。借助数据文件中的简单高角度速度模型，可直接运行观察高角度构造的成像效果；配合拉普拉斯滤波，还能进一步突出边界、提升剖面清晰度，是理解 RTM 原理与参数调试的实用入门资源。

1. 逆时偏移matlab代码：不是论文复现，是你明天就能跑通的最小骨架

搜“逆时偏移matlab代码”的人，多半是想给复杂构造做叠前深度偏移，但一翻参考书就被波动方程、成像条件劝退。我想说，逆时偏移（Reverse Time Migration，RTM）没有想象中那么玄，它的核心只是用同一个波动方程把震源波场正向传播一遍、把检波器记录反向传播一遍，再把两者在相同时刻的波场做乘积累加。对做勘探地震初值的工程师和研究者来说，在Matlab里先写一个二维声波RTM骨架，往往比直接调商业软件更能看清参数与成像结果之间的关系。这套代码适合验证算法、跑理论模型、给三维方案做原型，不适合直接把生产数据往里面灌。

2. RTM先把原理立住：两种波场、三个步骤、一个成像条件

逆时偏移能处理复杂构造成像，靠的不是某个复杂技巧，而是正、反两个波场在时间和空间上的“相遇”。这一章把原理拆成可执行的三段式流程，后面第三章的代码就按这个流程组织。

2.1 正向源波场与反向接收波场：两场相遇的地方就是界面

地震采集时，人工震源激发的波在地下传播，遇到波阻抗差异产生反射，反射波回到地表被检波器记录。RTM的思路是重新利用同一份速度模型：让源波场从震源点按时间正向传播，同时让检波器记录到的时间序列按时间反向“灌”回地下。当源波场与反向接收波场在某个空间点、某个时间步出现相位对齐时，这个点极有可能就是反射界面的像点。

这个原理可以写成最简洁的零延迟互相关成像条件：

I(x,z) = Σ_t S(x,z,t) R(x,z,t)

其中 S 是正向传播的源波场，R 是反向传播的接收波场。因为在反射界面上，两个波场的走时相等，乘积在时间上累加后自然会增强；而在非反射点，两个波场通常不会同步，累加值相对较小。这就是“两场相遇”的物理含义。

为什么要叫“逆时”？不是把时间倒转回去破坏物理规律，而是把接收记录当成边界驱动，让随时间循环从最后一炮道走向第一炮道。接收波场相当于在每个检波器位置向地下“重新回放”一段波，凡是能与源波场走时对齐的能量，最终都会在成像剖面上显示出来。这也解释了RTM为什么能成像回转波、棱柱波这类射线法难以描述的复杂波场。

2.2 为什么选二阶声波方程而不是一阶应力速度方程组

实现声波RTM有两条常用路线：二阶声波方程 p_tt = v²∇²p，和一阶速度-应力方程组。我一般会先在Matlab中用二阶声波方程做二维原型，原因很直接：变量少，内存占用小，向量化方便，调试时能一眼看懂更新式。

二阶方程只需要压力或位移一个标量场，时间递推用当前和前一个时刻两个场就够了，核心更新式写出来就是：

p_new = 2*p_now - p_old + v²·dt²·Laplace(p_now)

一阶速度-应力方程组在二维声波情形下至少需要压力和两个速度分量三个场，写PML吸收边界时虽然会更自然，但代码量更大，每个场的索引和边界条件都要单独处理，原型迭代周期会被拖长。

二阶方程的代价是：严格的PML需要把空间导数分裂成不同方向的分量，或者引入辅助场，代码复杂度明显上升。对于只想验证算法、跑理论模型的人来说，用指数衰减层做吸收边界已经能满足需求。这个取舍在后面第四章还会继续展开。

2.3 一套可执行的三段式流程

在贴代码之前，先把RTM的落地路径定下来。我建议所有第一次写的开发者都按下面三步走，不要一上来就并行做多炮：

1. 建立速度模型和震源子波，设定空间网格间距、时间步长和记录长度；
2. 正向传播震源波场，同时记录检波器位置的波场值作为单炮地震道，并把每个时间步的源波场快照保存下来；
3. 从最后一炮道开始，把地震记录反向注入波场，边逆推边与保存的源波场做零延迟互相关累加，最后做拉普拉斯滤波。

这套流程最容易被忽略的是第二步的波场保存策略。如果每个时间步都保存全波场，内存会很快耗尽；如果只保存边界，逆时重建时又需要额外写一套重放逻辑。我给出的第三章代码先用“保存全波场快照”的方式降低理解门槛，到了第五章再讨论如何在内存受限时改造。

3. 把RTM拆成能直接跑的Matlab代码：雷克子波、正演循环、逆时成像

这一章的代码按可读性优先，不追求性能最优。你可以把三个小节里的代码片段拼进同一个.m脚本，也可以按函数结构单独封装。运行前先确认Matlab当前目录可写，否则波场快照数组可能会把内存占满。

3.1 雷克子波与速度模型：最小可跑参数

RTM需要一个人工震源子波，最常用的是雷克子波。它的特点是主瓣尖锐、旁瓣有规律，很接近实际震源的频带特征。先写一个独立的子波函数：

function w = ricker(f0, t, delay) % f0 子波主频，单位Hz % t 时间向量，单位s % delay 子波峰值出现的延迟时间，单位s ta = t - delay; w = (1 - 2*(pi*f0*ta).^2) .* exp(-(pi*f0*ta).^2); end

雷克子波的关键参数是主频 f0。主频越高，时间上压缩越紧，对速度模型的细节分辨能力越强，但网格频散也越明显。delay 参数用来保证子波在正时间轴上有完整的波形，避免 t=0 时刻出现截断。实际使用中，我会让 delay 至少为 1/f0，这样子波起跳时不会打断波动方程的传播。

接着定义速度模型和观测系统参数：

% 速度模型与网格参数 dx = 10; dz = 10; % 空间步长，单位m nx = 201; nz = 201; % 横向/纵向网格数 x = (0:nx-1)*dx; z = (0:nz-1)*dz; % 坐标向量 v = 1500*ones(nz, nx); % 均匀背景速度 1500 m/s v(60:80, :) = 2500; % 中浅层水平高速体，模拟强反射面 dt = 0.001; tmax = 0.6; % 时间步长与记录长度 nt = round(tmax/dt); t = (0:nt-1)*dt; f0 = 25; delay = 0.1; % 雷克子波主频和延迟 w = ricker(f0, t, delay); sx = 100; sz = 30; % 震源位置（x网格号，z网格号） rx = [40 70 100 130 160]; % 5个检波器的x网格号 rz = 30*ones(1,5); % 检波器与震源同深度

这里要注意矩阵维度的约定：速度模型 v 采用“行对应 z、列对应 x”的二维矩阵，所以 v(60:80,:)=2500 表示 z 从 600m 到 800m 的一层高速体。震源放在 sz=30，也就是 300m 深度，位于高速体上方；检波器与震源同深度，模拟一道等偏移距接收的近似零炮检距记录。

3.2 正演循环：二阶有限差分与吸收边界

正演循环是整个RTM里最重的一段代码。它承担两个任务：一是产生单炮地震记录 seis，二是保存源波场快照 src_store 供后面成像使用。我的实现里用了一个简单的指数衰减边界，而不是完整PML，目的是先让代码逻辑清晰起来。

% 简单衰减边界：四周20个网格点做指数衰减 nb = 20; coeff = 0.03; amp = exp(-(coeff*(nb:-1:1)).^2); % 一维衰减系数 az = ones(nz,1); ax = ones(1,nx); for ii = 1:nb az(ii) = amp(ii); az(nz-ii+1) = amp(ii); ax(ii) = amp(ii); ax(nx-ii+1) = amp(ii); end alpha = az * ax; % 二维衰减剖面，nz×nx p0 = zeros(nz,nx); p1 = zeros(nz,nx); src_store = zeros(nz,nx,nt); % 保存源波场快照，教学用 seis = zeros(nt, numel(rx)); % 单炮地震记录 for it = 1:nt lap = zeros(nz,nx); lap(2:end-1,2:end-1) = p1(2:end-1,1:end-2) + p1(2:end-1,3:end) ... + p1(1:end-2,2:end-1) + p1(3:end,2:end-1) ... - 4*p1(2:end-1,2:end-1); lap = lap/(dx*dx); % dx=dz p2 = 2*p1 - p0 + (dt*dt)*(v.*v).*lap; p2(sz,sx) = p2(sz,sx) + (dt*dt)*w(it); p2 = p2 .* alpha; for kr = 1:numel(rx) seis(it,kr) = p2(rz(kr),rx(kr)); end src_store(:,:,it) = p2; p0 = p1; p1 = p2; % 时间推进 end

这里有几个参数必须解释清楚。

lap 使用的是五点拉普拉斯差分模板。索引写成 p1(2:end-1,1:end-2) 而不是直接对整矩阵操作，是为了避免把边界外的数据卷进来。边界一圈的 lap 保留为 0，也就是边界处不进入差分更新，再加上 alpha 衰减层，人工边界反射能被压到一个可接受水平。

时间递推式 p2 = 2*p1 - p0 + dt²·v²·lap 是二阶声波方程的标准显式格式。dt² 项不是随便加的，它和速度平方、拉普拉斯算子一起决定波场能否在网格上正确传播。如果去掉 dt²，波场会在几个时间步内直接发散成 NaN。

震源项在 p2 上叠加 dt²·w(it)，这样做是为了一阶方程的量纲对齐。实际实现中也可以不乘 dt²，只影响振幅，不影响同相轴位置，但我建议保留这个系数，便于之后做震源归一化时保持一致性。

alpha 衰减层的效果取决于 nb 和 coeff。nb 越大吸收带越宽，coeff 越大靠近边界的衰减越快。但衰减太快会把波“掰断”产生反射，所以 coeff 通常取 0.02 到 0.05 之间。更多细节在第四章继续展开。

这个循环结束后，seis 就是我们这一炮的合成记录，src_store 则是后续成像所需的源波场时间序列。如果读者手头有野外采集的实际道集，可以直接用真实数据替换 seis，但要注意真实数据里没有包含震源子波的正演过程，逆推前通常还要做带通滤波。

3.3 逆时外推与零延迟互相关成像

成像循环把时间从最后一刻往前扫描，每一步都把对应时刻的地震记录在检波器位置重新注入，同时将反向波场与提前保存的源波场做互相关叠加。代码结构与正演循环很接近，只是在末端多了一行累加式。

% 逆时偏移：把地震记录倒着注入，与源波场互相关 p0 = zeros(nz,nx); p1 = zeros(nz,nx); image = zeros(nz,nx); for it = nt:-1:1 lap = zeros(nz,nx); lap(2:end-1,2:end-1) = p1(2:end-1,1:end-2) + p1(2:end-1,3:end) ... + p1(1:end-2,2:end-1) + p1(3:end,2:end-1) ... - 4*p1(2:end-1,2:end-1); lap = lap/(dx*dx); p2 = 2*p1 - p0 + (dt*dt)*(v.*v).*lap; for kr = 1:numel(rx) p2(rz(kr),rx(kr)) = p2(rz(kr),rx(kr)) + (dt*dt)*seis(it,kr); end p2 = p2 .* alpha; image = image + src_store(:,:,it) .* p2; % 零延迟互相关累加 p0 = p1; p1 = p2; end

为什么在逆时循环里仍然使用和正演一样的差分格式？因为声波方程时间对称，把地震记录按时间倒序作为源项注入时，扩散方向就相当于波形从检波点反向往地下传播。它不是严格意义上的波动方程反演，而是同一个算子配合逆序驱动，这也是“反向回放”这个说法的由来。

seis(it,kr)在这里被乘上 dt²，与正演时的源加载保持同一尺度。如果你在正演时没有乘 dt²，这里也不要乘，否则振幅基准不一致。互相关累加式image = image + src_store(:,:,it) .* p2就是公式 I(x,z)=Σ_t S·R 的离散实现。

跑完这段代码后，image 里通常会出现很强的低频背景。不要急着把它当最终结果，先用 imagesc(x,z,image) 看一眼，你会发现反射同相轴基本能被辨认，但画面上有很多大面积低频斑块。这些低频能量来自零延迟互相关在波传播路径上的累积，处理方法是后文会讲的拉普拉斯滤波和震源归一化。

4. 参数边界与成本控制：CFL、衰减层取舍和网格间距怎么定

很多人拿到RTM代码后第一件事是改速度模型、跑数据，结果画面里全是噪声或干脆发散。问题往往不在成像条件，而在最基础的物理参数。这一章把最容易出问题的三个参数讲清楚。

4.1 衰减层参数不是越大越好：nb与coeff的关系

我在第三章代码里写的衰减层是一个按指数规律向模型内部衰减的乘性系数。边界最外侧系数接近 0，越往内越接近 1，形成一个不反射的“软墙”。很多初学者以为把nb加到50、coeff加到0.1就能把边界吸收得干干净净，实际上衰减系数过陡时，波场在衰减层内会产生新的反射，反而比衰减不够更难看。

我一般先用几组参数做对比试验。下面这组参数可以作为起点：

代码中生成alpha的部分可以单独提出来，方便调参：

nb = 30; coeff = 0.03; amp = exp(-(coeff*(nb:-1:1)).^2);

注意：alpha是矩阵乘到每个时间步波场值上的，所以它也会影响边界内部的波场幅值。如果模型很大且衰减层很宽，靠近边界的真实成像振幅会被压低。因此正式计算后要把边界部分裁掉，不要用边界附近的像点做岩性判断。

4.2 CFL条件：跑崩了先查这个

二阶显式有限差分对时间步长有硬性限制，通常叫CFL条件。Matlab里一旦出现NaN，最先检查的一定是dt是否过大，而不是成像条件写错。二维声波方程在dx=dz时，时间步长上限是：

dt_max = dx / (vmax · sqrt(2))

实际取值要留余量，我习惯取dt_max的0.9倍。可以在脚本开头加一个自动检查：

vmax = max(v(:)); dt_max = dx/(vmax*sqrt(2)); if dt >= dt_max error('当前dt太大，请把dt降到 %.6f s 以下', 0.9*dt_max); end

如果dx不等于dz，上限公式要改成：

dt_max = 1 / (vmax * sqrt(1/dx² + 1/dz²))

这个检查必须在正演循环之前执行，而不是等波场发散后再去查。还有一个容易踩的坑：速度模型里如果存在极大速度值，比如把某些网格设成了5000甚至更高，CFL上限会按你给定格式自动计算，但如果你在建模时不注意网格速度单位，也会让dt超过上限。

4.3 网格间距：一个波长至少要多少个点

RTM用网格波场近似连续介质，网格太粗会让波场产生数值频散。经验规则是：一个波长至少要有5个网格点，稳妥一点取10个。判断方法并不复杂：

vmin = min(v(:)); fmax = 2.5*f0; % 主频25Hz，高频端约50Hz dx_req = vmin/(fmax*10); if dx > dx_req fprintf('建议把dx降到 %.1f m，当前模型会有网格频散\n', dx_req); end

如果成像剖面上沿着强反射界面出现一连串高频颤动，像波浪尾巴一样贴在界面旁边，基本就是网格太粗导致的高频波场被离散化扭曲。处理方法有两种：一是把dx减小到满足一个波长10个点，二是降低震源主频。主频降下来会牺牲分辨率，但能快速判断问题方向。

我自己的习惯是：先用现有网格直接成像，如果频散不严重就不加密；如果严重，先降f0试算一炮，确认频散消失后再决定是否加密网格。这比一上来就把网格翻倍要省时间。

5. 逆时偏移避坑实录：低频噪音、内存翻车和边界反射

跑RTM第一次出图时，几乎每个人都会遇到下面几个问题。我把它们按“现象 → 原因 → 解决”的方式记录下来，遇到哪条就对哪条处理。

5.1 低频噪音铺满全图

现象：成像图上除了反射同相轴，还有大面积的黑白斑块，低频背景把所有浅中层信息压住，整张图像蒙了一层雾。

原因：零延迟互相关成像条件不仅反射界面处有能量，在源波场和反向波场的传播路径上也会产生噪声。尤其浅层，两个波场都走过，互相关值被反复累加，形成低波数背景。这是RTM最常见的现象，不是代码bug。

解决：先用拉普拉斯滤波把低频背景压下去。简单做法是直接对image取负拉普拉斯：

lap_img = zeros(size(image)); lap_img(2:end-1,2:end-1) = image(2:end-1,1:end-2) + image(2:end-1,3:end) ... + image(1:end-2,2:end-1) + image(3:end,2:end-1) ... - 4*image(2:end-1,2:end-1); image = -lap_img;

这个操作会增强高频噪声，所以在滤波之后最好补一个中值滤波或带通滤波。很多商业软件里的“Laplacian filter”选项，本质就是这一行。

5.2 内存直接溢出

现象：脚本跑到正演循环时，Matlab提示内存不足，尤其是当我把nx、nz调到200以上，nt调到1000以上后，连保存src_store这一步都过不去。

原因：src_store(:,:,:)是一个三维数组，尺寸是nz×nx×nt，double类型下每个数字占8字节。200×200×1000就有32亿个元素，约25.6GB内存，一般工作站根本扛不住。

解决：教学场景下先砍模型规模，比如把nx、nz降到100，nt降到400，跑通后再放大。生产级做法是只保存边界场，而不是全波场。伪代码如下：

% 正演时只保存边界 edge_store(1,1,:) = p2(1,:); % 上边界 edge_store(1,2,:) = p2(nz,:); % 下边界 edge_store(2,1,:) = p2(:,1); % 左边界 edge_store(2,2,:) = p2(:,nx); % 右边界 % 逆推时用边界值重建源波场，再与反向波场互相关

边界重建需要对差分算子做时间反向重放，代码量会比全波场保存大得多。如果你只是验证算法，建议先接受全波场保存法的内存消耗，把模型调小，而不是一开始就写边界重建。

5.3 边界反射像鬼影

现象：成像图四边出现弧形强干扰，像用圆规画出的弧线，高速层越浅，弧线越明显。

原因：衰减层没有把边界吸收干净，波遇到人工边界后反射回模型，反射波在互相关累加里被当成有效能量成像。

解决：三个方向同时调。第一，把nb从20加宽到40；第二，把coeff从0.03调小到0.02；第三，给模型四周扩边后再跑，成像后裁掉扩边部分。扩边代码很简单：

v_pad = padarray(v, [30 30], 'replicate'); % 四周扩30个网格 % 跑完RTM后裁回 image = image_pad(31:end-30, 31:end-30);

注意扩边时要用'replicate'，也就是用最外层速度值向外复制，不能用0填充。如果填充0，波到速度突变处照样会产生强反射，等于白扩。

5.4 成像位置对，但像点旁瓣多

现象：反射界面在主像点上下出现一串等间距同相轴，看起来像多个界面叠在一起。

原因：雷克子波本身有正负旁瓣，零延迟互相关会把子波的旁瓣也保留在成像值里。另外，单炮成像的振幅差异很大，直接把多炮像叠加，强炮旁瓣会盖住弱炮主瓣。

解决：用震源归一化替代简单互相关。分母用源波场能量做补偿：

denom = sum(src_store.^2, 3) + eps; image = image ./ denom;

这样做能明显压旁瓣和浅层低频，但也会放大低信噪比区域的噪声。更稳妥的做法是在分母里加一个正则化项，比如取0.001倍的最大能量：

denom = sum(src_store.^2, 3); image = image ./ (denom + 0.001*max(denom(:)));

正则化参数看着像玄学，实际上它控制的是“对弱能量区域放大的程度”。我会先用不归一化的像判断同相轴位置，再用归一化像做最终显示，两者结合着看。

6. 进阶验证技巧：拉普拉斯滤波、震源归一化和单点模型自检

前五章把RTM打通了，但通不代表对。最后一章讲三个我每次跑新模型都会做的动作，它们能帮你快速判断成像结果是否可信。

6.1 拉普拉斯滤波的正确用法

拉普拉斯滤波不是锐化工具，它的本质是去掉成像剖面里的低频背景。使用前先看一眼image直方图，如果低频斑块占据很大面积，再决定滤波强度。滤波后一定要用同样的色标显示前后对比，不能只看滤波后的图，否则容易被“高对比度假清晰”欺骗。

imagesc(x, z, image); axis image; colormap gray

如果滤波后的真实反射界面变碎、变细，说明滤波过度，可以改用带通滤波代替。带通滤波保留中频反射主体，只压零频和最低频成分，效果通常会比纯拉普拉斯更自然。

6.2 用单散射点模型校验定位

在跑复杂构造之前，先在均匀背景里放一个小扰动体，比如在第三章的速度模型里加一个点状高速块：

v(100:103, 90:93) = 1700; % 在1500背景中放点散射体

跑完RTM后，成像亮点应当落在散射体上，深度误差不超过一个网格。如果亮点整体偏移半个波长以上，第一检查子波delay是否扣除；第二检查震源位置和检波器位置是否标错。这步是RTM里最重要的自检，它能直接暴露坐标系和子波定义的问题。

6.3 多炮叠加前的归一化

多炮RTM的结果是每炮成像的叠加。单炮像之间振幅差异很大，直接叠加往往让浅层强反射主导整张图。我的习惯是每炮先除以自身能量，再叠加：

img_shot = img_shot / (norm(img_shot(:)) + 1e-6); final_img = final_img + img_shot;

这个做法能让深层弱反射也有机会被看到，但代价是噪声也被同步放大。所以叠加后我会再做一次拉普拉斯滤波，把被放大的高频噪声压回去。

RTM这个方向，代码只是起点，真正花时间的是把边界、频散、成像条件三个环节逐一验证。我自己每次跑新模型都坚持先单点、再多层、后复杂构造的顺序；每改一个参数只跑单炮，不跑叠加，只有确认单炮像的第一个亮点位置正确后，才会继续加炮。这样排错时，你不会在“震源子波错了还是边界泄漏”之间反复猜。希望帮到你。

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