简介：这份资源面向地震勘探与地球物理方向的学习者和研究人员，提供利用MATLAB实现地震波偏移成像的完整代码集合，帮助校正反射点位置、消除波传播导致的图像失真，从而更准确地揭示地下地质构造。压缩包共24个文件，全部为m脚本，整体约50KB，涵盖傅里叶变换、逆时迁移、频率-空间域偏移、数据预处理与成像可视化等环节，并包含速度模型构建、时间-深度转换及多种偏移算法实现，便于读者结合自身地震数据开展实验与算法对比研究。目前已有2010人学习下载，适合具备一定MATLAB编程基础、希望深入理解偏移成像原理并动手实践的中高级用户参考。

1. 地震波偏移成像：从叠后数据到地下构造的MATLAB实现

拿到一份叠后地震数据，想看看地下构造长什么样，绕不开的一步就是偏移成像。简单说，偏移就是把地震记录中因绕射、倾斜反射等原因偏离真实空间位置的能量，重新归位到它该在的地方。没有偏移，你看到的剖面就像一张被揉皱的地图，断层、尖灭、盐丘边界全糊在一起。这套MATLAB代码包干的就是这件事：输入叠后时间剖面，输出偏移归位后的深度/时间域图像。适合谁？地球物理专业的学生做课程设计、科研院所的初级处理人员快速验证算法、以及想用MATLAB复现经典偏移方法的开发者。它不依赖商业软件，纯代码实现，改起来比黑匣子灵活得多。

2. 技术选型与目录结构：为什么用MATLAB做偏移

2.1 选MATLAB而不是C/Python的实在理由

地震波偏移成像的核心是波场延拓和成像条件，计算量大，但算法验证阶段更看重矩阵操作和可视化效率。MATLAB在这两点上几乎是“开箱即用”：fft2、ifft2、稀疏矩阵、imagesc直接出图，省去大量底层胶水代码。这套代码包大概率围绕几种经典方法展开——相移偏移（Phase-Shift）、分步傅里叶偏移（SSF）、或者有限差分法。相移法在频率-波数域做延拓，适合横向速度变化不大的介质；分步傅里叶在相移基础上加了空间域的相位校正，能处理中等横向变速。代码里通常会有对应的主脚本和函数文件，而不是一个巨型.m文件从头写到尾。

常见做法是：主脚本负责读数据、设参数、调函数、画图；函数文件按偏移方法拆开，比如phase_shift_migration.m、ssf_migration.m、fd_migration.m。这样你换方法时只改一行调用，不用动数据加载和显示部分。

2.2 典型目录结构长什么样

虽然没有拿到实际压缩包，但按这类项目的惯例，目录层级不会太深，大致如下：

seismic_migration_matlab/ ├── main.m % 主入口：读数据、选方法、跑偏移、出图 ├── data/ │ └── stacked_section.mat % 叠后地震数据，通常是二维矩阵 ├── src/ │ ├── phase_shift.m % 相移偏移核心函数 │ ├── ssf.m % 分步傅里叶偏移 │ ├── fd_migration.m % 有限差分偏移 │ ├── extrapolate.m % 波场延拓公共模块 │ └── imaging_condition.m % 成像条件（零时刻/互相关） ├── utils/ │ ├── read_segy.m % 读SEGY格式（如果支持） │ ├── taper.m % 边界吸收/衰减 │ └── plot_section.m % 统一绘图函数 └── results/ └── migrated.mat % 偏移结果保存

main.m里通常有清晰的参数区，比如dx、dt、fmax、vel、method。你拿到手第一件事就是打开main.m，看它默认读哪个数据、用什么速度、跑哪种方法。如果数据是.mat，直接load；如果是SEGY，得确认read_segy.m能不能解析你的文件头。

2.3 数据输入输出的约定

叠后数据一般是二维矩阵，行对应时间采样，列对应道号。时间采样间隔dt单位通常是秒或毫秒，空间采样间隔dx单位是米。速度模型可能是一个常数，也可能是一个与深度相关的层状模型，甚至是一个二维速度场。代码里如果支持变速，速度矩阵的尺寸必须和数据矩阵匹配，否则延拓时索引会越界。

输出通常是偏移后的图像矩阵，以及配套的坐标轴信息。plot_section.m一般会做归一化、加色标、标注时间/深度轴。有些实现还会把偏移前后的剖面对比画在一起，方便你直观判断绕射是否收敛、断层是否清晰。

注意：如果代码包里没有data/目录，说明数据需要你自己准备。常见做法是用MATLAB内置的seismic示例数据，或者用简单的合成模型（如点绕射、倾斜层）先跑通流程。

3. 核心模块拆解：相移偏移与分步傅里叶的实现细节

3.1 相移偏移的频域延拓逻辑

相移偏移的核心思想：在频率-波数域，波场延拓的相位因子是解析的。对二维情况，延拓步长dz对应的相位为：

phase = exp(1i * dz * sqrt((omega/v)^2 - kx^2))

其中omega是角频率，kx是横向波数，v是介质速度。代码里通常先对数据做二维FFT，然后对每个频率切片乘以相位因子，再反变换回空间域。下面是一个简化版的相移延拓函数：

function wavefield = phase_shift_extrap(wavefield, dx, dt, dz, v) % 输入：wavefield 当前波场（时间-空间域） % 输出：延拓dz后的波场 [nt, nx] = size(wavefield); % 转到频率-波数域 FK = fftshift(fft2(wavefield)); % 频率轴和波数轴 omega = 2*pi*[(0:nt/2-1), (-nt/2:-1)]' / (nt*dt); kx = 2*pi*[(0:nx/2-1), (-nx/2:-1)] / (nx*dx); [KX, OMEGA] = meshgrid(kx, omega); % 相移因子，只对传播波（sqrt内为正）做延拓 kz2 = (OMEGA/v).^2 - KX.^2; kz = sqrt(max(kz2, 0)); % 倏逝波置零 phase = exp(1i * dz * kz); % 延拓并反变换 FK_extrap = FK .* phase; wavefield = real(ifft2(ifftshift(FK_extrap))); end

逻辑说明：fft2把波场转到频率-波数域，fftshift把零频移到中心。omega和kx的构造要保证正负对称，否则相位因子会错。kz2为负时对应倏逝波，物理上不传播，直接置零避免数值爆炸。real取实部是因为原始数据是实信号，反变换后虚部由数值误差引入。

参数怎么改：dz是延拓步长，越小精度越高但计算越慢，一般取v/(4*fmax)到v/(2*fmax)之间。v如果是常数，直接填数值；如果是变速，这个函数就不够用了，得换分步傅里叶或有限差分。

3.2 分步傅里叶偏移的补偿项

分步傅里叶偏移在相移基础上，把速度分解为背景速度v0和扰动dv。背景部分在频率-波数域做相移，扰动部分在空间域做时移。核心代码片段：

function wavefield = ssf_extrap(wavefield, dx, dt, dz, v0, dv) % 背景相移 wavefield = phase_shift_extrap(wavefield, dx, dt, dz, v0); % 扰动时移：在空间域乘以相位 [nt, nx] = size(wavefield); t = (0:nt-1)' * dt; shift_phase = exp(1i * 2 * pi * dv(:)' * t / v0^2); % 近似时移 wavefield = wavefield .* shift_phase; end

逻辑说明：dv是速度扰动，通常是一个与空间位置有关的向量。shift_phase的推导基于Born近似，把速度扰动对波场的影响近似为时移。这里dv(:)'转成行向量，t是列向量，乘出来是nt×nx的矩阵，逐元素乘到波场上。

参数怎么改：v0一般取速度模型的平均或平滑背景，dv = v - v0。如果dv太大，这个近似会失效，偏移结果会出现假象。常见做法是限制dv/v0不超过0.2。

3.3 成像条件与边界处理

成像条件通常用零时刻成像：波场延拓到每个深度后，取t=0时刻的波场值作为该深度的成像结果。代码里就是每延拓一步，把wavefield(1,:)存到成像矩阵的对应行。

边界处理是偏移里的“玄学”环节。不做边界吸收，波场会在边缘反射，偏移剖面上出现弧形假象。常见做法是加一个 taper 窗，或者用吸收边界条件。下面是一个简单的 taper 函数：

function w = taper2d(w, ntaper) % 对二维波场四周做余弦衰减 [nt, nx] = size(w); taper_t = ones(nt, 1); taper_t(1:ntaper) = 0.5*(1 - cos(pi*(0:ntaper-1)'/ntaper)); taper_t(end-ntaper+1:end) = flipud(taper_t(1:ntaper)); taper_x = ones(1, nx); taper_x(1:ntaper) = 0.5*(1 - cos(pi*(0:ntaper-1)/ntaper)); taper_x(end-ntaper+1:end) = fliplr(taper_x(1:ntaper)); w = w .* (taper_t * taper_x); end

逻辑说明：taper_t和taper_x分别是时间和空间方向的余弦窗，外积后乘到波场上。ntaper一般取20到50个采样点。这个操作在每次延拓前做一次，能显著压制边缘反射。

提示：如果偏移结果边缘有强能量弧形，先检查taper有没有加，再加检查taper宽度够不够。

4. 跑通流程：从数据加载到偏移出图

4.1 准备数据与参数文件

第一步是确认数据格式。如果代码包自带.mat文件，直接load即可。如果是SEGY，用read_segy.m读，注意检查采样间隔和道间距是否和代码里的dt、dx一致。不一致的话，要么改代码参数，要么对数据重采样。

参数文件通常是一个结构体或脚本里的变量区。关键参数包括：

dz如果取太大，相位因子会混叠，偏移结果出现假频。一个血泪经验：先拿一个点绕射模型试，偏移后点绕射应该收敛成一个点，如果散成十字或弧形，多半是dz或fmax设错了。

4.2 主脚本调用与运行

main.m的典型结构：

%% 加载数据 load('data/stacked_section.mat'); % 变量名假设为 data [nt, nx] = size(data); %% 参数设置 dt = 0.002; dx = 20; fmax = 40; v = 2500; dz = v/(4*fmax); ntaper = 30; %% 边界衰减 data = taper2d(data, ntaper); %% 偏移 method = 'phase_shift'; % 可选 'ssf', 'fd' switch method case 'phase_shift' image = phase_shift_migration(data, dx, dt, dz, v); case 'ssf' image = ssf_migration(data, dx, dt, dz, v0, dv); case 'fd' image = fd_migration(data, dx, dt, dz, v); end %% 显示 figure; subplot(1,2,1); imagesc(data); title('偏移前'); colormap(gray); axis tight; subplot(1,2,2); imagesc(image); title('偏移后'); colormap(gray); axis tight;

运行直接敲main。如果报错“未定义函数”，检查src/和utils/有没有加到MATLAB路径里，用addpath(genpath('src'))和addpath(genpath('utils'))。

4.3 结果判读与常见现象

偏移后的剖面，绕射波应该收敛，倾斜反射归位到正确倾角，断层更清晰。如果偏移后还不如偏移前，先别怀疑算法，按这个顺序查：速度对不对、dz是不是太大、taper够不够、数据有没有做过静校正。速度给错是最常见的翻车原因——速度偏大，同相轴过偏；速度偏小，同相轴欠偏。

注意：偏移不是万能的。如果原始数据信噪比太低，偏移会把噪声也归位，结果可能更乱。先做去噪再偏移。

5. 避坑与排查：偏移成像里最容易翻车的几个点

5.1 偏移后同相轴位置不对

现象：偏移后反射轴偏离了地质认识上的位置，或者倾斜层倾角变了。原因：速度模型不准。偏移对速度非常敏感，速度误差10%就能让同相轴偏移几十米。解决：用速度谱或层析反演重新拾取速度，或者先用常速扫描，看哪个速度能让绕射收敛最好。

5.2 边缘出现强弧形能量

现象：剖面左右边缘有规律的弧形同相轴。原因：边界反射。波场延拓到边界时没有吸收，能量反射回来。解决：加taper，或者用PML吸收边界。taper宽度至少20个点，太窄了吸收不干净。

5.3 偏移结果出现假频

现象：剖面上有细密的条纹或棋盘格。原因：延拓步长dz太大，或者fmax设得过高，导致相位因子混叠。解决：dz满足dz <= v/(4*fmax)，fmax取数据有效频带的上限，不要盲目设高。

5.4 内存溢出或计算太慢

现象：跑大剖面时MATLAB卡死或报“Out of Memory”。原因：二维FFT对大数据矩阵很吃内存，尤其是三维数据。解决：分块处理，把数据按道分成若干块，逐块偏移后拼接。或者降采样，但降采样会损失高频。

5.5 成像条件取错时刻

现象：偏移结果整体模糊，同相轴不聚焦。原因：成像条件没有取零时刻，或者取错了索引。解决：检查成像条件代码，确认每延拓一步取的是wavefield(1,:)，而不是中间某个时刻。

6. 二次开发：换速度模型、换偏移方法、加并行

6.1 从常速扩展到变速

常速相移偏移的代码改造成变速，最直接的路子是换分步傅里叶。把v拆成v0和dv，v0取平滑背景，dv做扰动补偿。代码改动集中在延拓函数，主脚本里多传一个dv参数。如果横向速度变化剧烈，分步傅里叶也不够，得上有限差分或逆时偏移，但计算量会大一个量级。

6.2 加一个自己的偏移方法

代码包的switch结构就是给你留的口子。新建一个my_migration.m，输入输出接口和现有函数保持一致，然后在main.m里加一个case 'my_method'。接口约定：输入data, dx, dt, dz, v，输出image矩阵，尺寸和data一致。这样你可以在不改动其他部分的情况下，对比不同方法的效果。

6.3 用parfor加速频率循环

相移偏移里对每个频率切片的处理是独立的，可以用parfor并行。把频率循环从for改成parfor，前提是循环内没有跨迭代的依赖。改之前先确认MATLAB Parallel Computing Toolbox可用。加速比大概等于核数，但要注意内存开销也会翻倍。

% 串行版 for iw = 1:nt FK(iw,:) = FK(iw,:) .* phase(iw,:); end % 并行版 parfor iw = 1:nt FK(iw,:) = FK(iw,:) .* phase(iw,:); end

逻辑说明：parfor把频率循环分到多个worker上，每个worker处理一部分频率。FK和phase是切片变量，MATLAB会自动处理数据分发。注意parfor里不能有break、return，也不能动态扩展数组。

6.4 验证改造是否正确的笨办法

改完代码别急着跑实际数据，先拿合成模型试。点绕射模型：一个点绕射偏移后应该收敛成一个点。倾斜层模型：偏移后倾角应该和模型一致。这两个测试过了，再上实际数据。我一般会在results/里存一份基准结果，每次改完代码跑一遍合成模型，和基准对比，误差超过5%就回退。从那以后我每次改偏移代码都强制走一遍合成模型测试，省得在真实数据上浪费时间。希望帮到你。

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