简介：面向地震勘探资料处理学习者与地球物理专业师生的MATLAB偏移成像代码合集，聚焦校正反射点位置、消除波传播造成的成像畸变，适合具备一定MATLAB基础、希望从脚本层面理解偏移算法实现的中高级读者。压缩包共24个文件，全部为.m源码，体积约50KB，内容涵盖频率-波数域偏移、有限差分偏移、分裂步偏移、基尔霍夫偏移以及时间-深度转换、锥形窗、数据预处理等辅助模块，结构紧凑，便于按算法逐项研读。已有2015人学习下载，说明其在课程实验与小规模数据验证中具有一定参考价值。读者可据此对照波动方程偏移的基本流程，理解频域变换、波场延拓与成像条件在代码中的落地方式，并借助不同脚本比较各类方法的适用条件，为后续使用实际地震数据开展偏移实验、调试参数与可视化分析提供可复用的起点。

1. 地震波偏移成像在 MATLAB 里跑通，多数人卡在一个很具体的画面：叠加剖面上还算连续的反射同相轴，一偏移就变成一圈圈弧形假象，深部糊成一片灰。偏移成像要干的事很朴素——按速度模型把被绕射和倾斜地层打散的能量，重新放回它真正的空间位置。手上已有观测记录或正演数据、会写 MATLAB 循环但没系统做过波动方程成像的人，可以照着搭出一套二维声波逆时偏移的最小闭环，顺带知道哪几个参数一动就翻车。

2. 有限差分逆时偏移的最小可跑闭环：从速度模型到成像剖面

在 MATLAB 里做地震波偏移成像，核心从来不是「会不会写 FFT」，而是能不能把四个环节串成一个自洽的循环：正演震源波场、记录地表数据、把数据倒着注回模型、在相同时刻把两个波场相乘累加。这一段把可运行的最小版本直接写出来，网格约定、量纲、存储策略、成像条件全部定死，后面几章再往上加精度和存储优化。

2.1 网格约定与速度模型：先把坐标和量纲确定下来

MATLAB 的矩阵索引是「行、列」，地震波场则习惯「深度、水平距离」。两者能对上：行索引对应深度 z 方向（向下为正），列索引对应水平 x 方向。一旦这个约定定下来就别再改，我自己血的教训是中途换了一次行列约定，反传波场和正传波场对不上，成像结果靠近对称位置，排查了大半天才找到。

量纲统一用 SI：速度 m/s，网格间距 m，时间 s，频率 Hz。别混用 km 和 m，也别让速度写 2.0 而网格写 5，这类错误不会报错，只会让成像结果整体偏移。

速度模型对偏移来说必须是平滑的。偏移只负责把能量归位，不负责反演速度；把速度模型建得过于精细反而容易把绕射能量打散，成像剖面反而更乱。常见做法是用一个背景速度加几个高斯型异常体，既可控又能验证算法：

% 网格与速度模型：行=深度，列=距离，单位统一 m nx = 401; nz = 201; dx = 5; dz = 5; % 网格间距 5 m x = (0:nx-1)*dx; z = (0:nz-1)*dz; [Z,X] = ndgrid(z,x); % 生成 nz x nx 的坐标网格 % 背景 1800 m/s，叠加一个高速高斯异常体 v = 1800 + 700*exp(-((X-850).^2/(2*150^2) + (Z-420).^2/(2*90^2))); v(1:24,:) = 1500; % 顶部低速层，模拟近地表

速度场用高斯函数，而不是矩形块，原因很实际：矩形块的四个角是速度的阶跃面，正演时会在角点产生强烈的绕射，成像剖面上会出现四条从角点发散出去的弧，干扰判断。用平滑的异常体，绕射能量集中在异常体本身，看起来更接近真实地质体的响应。150和90分别是横向、纵向的高斯半径，控制异常体的胖瘦，调这两个数就能做出不同尺度的目标体。

2.2 二维声波方程的 4 阶差分算子与单步更新函数

二维常密度声波方程写成u_tt = v²(u_xx + u_zz)，时间方向用 2 阶中心差分，空间方向用 4 阶中心差分。空间阶数选 4 阶，而不是 2 阶，是精度和成本之间的折中：2 阶差分在同样的网格间距下网格频散非常明显，剖面上会出现一串沿传播方向拖尾的「波纹」；4 阶差分的算子系数是(-1/12, 4/3, -5/2, 4/3, -1/12)，同样网格下频散误差降低一个数量级以上，代价只是每个点多重两个点，在 MATLAB 里就是多写两行索引。

function un = step2d(u, up, v, dt, dx, dz, damp) % 二维声波方程单步推进：时间 2 阶、空间 4 阶 % u : 当前时刻波场 [nz x nx] % up : 上一时刻波场 [nz x nx] % v : 速度模型 [nz x nx]，m/s % damp : 每步乘的吸收系数 [nz x nx]，内部为 1 c = (v*dt).^2; L = zeros(size(u)); L(3:end-2,3:end-2) = ... (-1/12*(u(3:end-2,1:end-4) + u(3:end-2,5:end)) ... + 4/3 *(u(3:end-2,2:end-3) + u(3:end-2,4:end-1)) ... - 5/2 * u(3:end-2,3:end-2)) / dx^2 ... + (-1/12*(u(1:end-4,3:end-2) + u(5:end,3:end-2)) ... + 4/3 *(u(2:end-3,3:end-2) + u(4:end-1,3:end-2)) ... - 5/2 * u(3:end-2,3:end-2)) / dz^2; un = (2*u - up + c.*L) .* damp;

索引3:end-2是内点范围，边界两圈留给吸收层，不做差分计算。un计算出来之后马上乘以damp，把边界附近的能量逐步衰减掉——这一步放在最后，是为了让差分算子始终作用在未经衰减的场值上，避免吸收系数污染内点。

这个函数每调用一次推进一个时间步，是整个程序里最热的代码。MATLAB 里做这类循环，有几个提升速度的做法值得早点用上：把u、up、v、damp全部预先转成single；L不要每次zeros重新分配，改成一个在外层声明好的持久变量；把c = (v*dt).^2提到循环外面算好，因为它不随时间变。

2.3 正传记录与反传成像：把记录倒着注回去

逆时偏移的物理图像很直白：正传时震源波场从震源向外扩散，遇到界面产生反射，反射被检波器记录下来；反传时把记录当成震源，从检波点位置按时间倒序注入，波场会沿着原来的路径反向传播回去，在反射界面处「聚焦」。把正传波场和反传波场在同一时刻、同一位置相乘再累加，聚焦处乘积大，其他位置乘积小，成像条件就成立。

观测系统设置如下。把震源和检波点放在同一条水平线上，相当于一条海底电缆，这样能避开自由表面产生的强烈多次波干扰，把注意力集中在成像本身。

正传阶段把每一步的波场抽稀存储，抽稀间隔dec = 4。这一步是后面第 4 章存储优化的引子，先按最省事的方式实现，把成像闭环跑通。

% ---- 正传：记录数据并抽稀保存波场 ---- dec = 4; % 每 4 步存一次正传波场 ns = floor(nt/dec); Sf = zeros(nz, nx, ns, 'single'); rec = zeros(nt, nrec, 'single'); u1 = zeros(nz,nx,'single'); u2 = zeros(nz,nx,'single'); for it = 1:nt un = step2d(u1, u2, v, dt, dx, dz, damp); u2 = u1; u1 = un; u1(isz,isx) = u1(isz,isx) + src(it); % 点源注入 rec(it,:) = u1(rz, rxs); if mod(it,dec) == 0 Sf(:,:,it/dec) = u1; end end

u2 = u1; u1 = un;这两行的顺序不能反。step2d需要「当前」和「上一时刻」两个场，推进完之后旧的当前场变成上一时刻场，新算出来的场成为当前场。写反了程序不报错，但波场会原地不动或者发散，属于典型的静默失败。

源注入直接加src(it)，不乘dt²/(dx·dz)。严格来说这个常数因子是量纲需要的，但成像条件和最终的归一化会把它消掉，为了代码可读性先略去。如果要做振幅保真的反演类工作，这个因子就必须补回去。

反传阶段把记录倒序注入，同时在抽稀的时间点上做互相关累加：

% ---- 反传：记录倒序注入，遇抽稀点做互相关 ---- ur1 = zeros(nz,nx,'single'); ur2 = zeros(nz,nx,'single'); img = zeros(nz,nx); for it = nt:-1:1 unr = step2d(ur1, ur2, v, dt, dx, dz, damp); ur2 = ur1; ur1 = unr; ur1(rz, rxs) = ur1(rz, rxs) + rec(it,:); if mod(it,dec) == 0 img = img + Sf(:,:,it/dec) .* ur1; % 零延迟互相关 end end

mod(it,dec)==0在倒序循环里同样成立，it/dec从nt/dec递减到 1，正好和正传时的存储顺序一一对应。这里只用了 1/4 的时刻做互相关，成像项少了四分之三，但因为互相关函数在时间上是平滑的，少几个采样点带来的误差远小于视觉可辨的程度。真要较真，把dec改成 1 再跑一遍对比，会发现差异主要出现在深部弱能量区。

2.4 照明补偿与 Laplace 滤波：让深部能量能看

直接算出来的img几乎不能看。两个原因：一是几何扩散让深部能量比浅部弱好几个数量级，画出来只有顶部一条亮线；二是零延迟互相关成像条件本身会产生强烈的低频背景，表现为剖面上大片的灰雾和界面处的「画弧」。

照明补偿解决第一个问题。把正传波场的能量在时间上累加，得到每个点的照明强度，再用它去除成像结果：

illum = sum(Sf.^2, 3); % 正传照明强度 img = img ./ (illum + 0.02*max(illum(:))); % 加小量防止除零

分母里的0.02*max(illum(:))是稳定项，取值太小吃不消深部极弱的照明，除法会把噪声放大成雪花；取得太大又退化成不做补偿。0.02是经验值，速度横向变化剧烈时可以往上调到0.05。注意照明补偿只能压平能量，不能凭空造出信息，照明为零的地方补偿之后还是一片噪声，这是方法的边界。

第二个问题用 Laplace 滤波处理，本质是高通：把成像结果减去它自己的大尺度平滑版本，留下界面附近的尖锐部分。

k = ones(15)/225; % 15x15 盒式平滑核 imgl = conv2(conv2(img, k, 'same'), k', 'same'); imgf = img - imgl; % 高通后的成像结果 imgf = imgf / max <p> <a href="https://download.csdn.net/download/u014775170/8849083" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>