做地震数据处理的同行应该都有过这种经历：野外采集回来的一条测线，因为过沟、过村庄、设备故障或者遇上坏道，最终叠前数据里总是缺那么几道。缺道少道看着是小问题，但到了偏移成像或者AVO分析阶段，空道带来的采空效应和空间假频会让人非常头疼。传统做法是用线性插值、f-k插值或者反漏频（anti-leakage）方法去补，这些方法在缺口小、倾角缓的时候还凑合，一旦出现连续几十道缺失，插值结果基本就是一片模糊，反射同相轴被抹平，波场特征也丢了。
压缩感知（Compressive Sensing, CS）的思路完全不同，它换个角度问问题：与其在“补齐”之后再做去噪，不如直接在重建过程中利用“地震数据本身可以被稀疏表示”这个先验。只要数据在某个变换域（比如f-k域、curvelet域）足够稀疏，那么即使在空间方向严重欠采样，也有可能通过非线性优化把完整波场恢复出来。我去年在一个二维地震工区的规则化处理里完整走了一遍这套流程，Matlab实现下来效果很稳，这里把原理、代码、调参经验和踩过的坑一次性分享出来。
这篇内容适合三类人：正在做地震数据规则化或者缺失道重建的同行、对压缩感知算法感兴趣但不知道如何落地到实际数据的同学、以及已经有Matlab基础想找一份能直接跑起来的核心代码作为起点的人。
1. 从缺道到重建：为什么传统插值在面对地震数据时不够用
1.1 缺道问题的真实形态：不只是“少了一道”
野外实际数据里的缺道和我们在教科书上画的“规则采样后挖掉几个点”差别很大。一条二维测线拿到手上，常见的情况是：
- 个别检波点因为埋置条件差，记录完全作废，表现为孤立坏道；
- 炮点附近有障碍物或者穿越公路，一连十几炮的排列无法铺开，形成缺口段；
- 3D观测系统的边缘区域本身就不完整，规则化后会发现大面积的不规则网格。
这些缺缝在时间-空间域未必很显眼，但变换到f-k域看，能量会沿波数轴涂抹开，形成明显的假频泄露。尤其是连续缺道的时候，缺失宽度在空间上超过一个空间采样间隔的几倍以上，重建难度会指数上升。
传统线性插值的问题在于，它只利用了数据和邻道之间的局部相关性，本质上是在做“平滑补齐”。当地震数据含有多个不同视速度的同相轴时，线性插值无法区分交叉能量，结果就是同相轴之间互相污染。f-k插值虽然把数据变到频率-波数域做预测，但它的前提是数据在f-k域呈现可预测性，一旦空间采样不规则，这类方法的效果会大打折扣。
1.2 压缩感知提供了什么不同视角
压缩感知的基本假设是：一段信号本身虽然看起来维度很高，但在某个线性变换下可以用很少的系数近似表达。对地震数据来说，这个假设非常自然——波场是由有限个反射同相轴构成的，每个同相轴在f-k域就是一条沿特定方向的能量线，在curvelet域则对应少量尺度-方向系数。
基于这个假设，重建问题就从“插值”变成了“稀疏系数恢复”：知道观测算子（即哪些道缺失，哪些道保留），找一个稀疏系数向量，让它在经过逆变换之后与观测数据吻合。数学上写成优化问题就是：
[ \min_{\mathbf{x}} |\mathbf{x}|_1 \quad \text{s.t.} \quad |\mathbf{A}\mathbf{x} - \mathbf{b}|_2 \le \epsilon ]
这里的A是观测矩阵，由稀疏变换和采样掩码复合而成；b是实际观测到的残缺数据。只要稀疏变换选择恰当，这个看似简单的优化问题能够解决传统方法无能为力的连续缺道场景。这也是为什么最近几年规则化处理和一发及多发的混采分离都在往CS框架上靠。
2. 重建算法落地的数学地基：稀疏表示、观测矩阵、迭代收缩
2.1 选哪个稀疏变换域：f-k 域与 curvelet 域
稀疏变换的选择直接决定重建质量的上限。如果你面对的地震数据主要包含线性或近似线性的反射同相轴，f-k域（二维傅里叶变换）就足够了。但是如果数据里有绕射波、断层附近的复杂波场，或者你处理的是叠后数据且地层倾角变化很快，f-k域的稀疏性会明显下降，这时候curvelet域的表现更好，因为curvelet同时具备方向性和多尺度特性，对曲线状同相轴也能给出稀疏表达。
我自己的习惯是：第一步先看数据集的f-k谱，如果能量集中在少数几条线上，直接用二维FFT做稀疏变换，代码简单、迭代快。如果f-k谱是弥散的，切换curvelet工具包比较稳妥。需要注意的是，curvelet变换的Matlab实现有多个版本，接口不一，用的时候务必确认正向变换和逆变换是否严格互逆，否则重建结果里会出现系统性的幅度畸变。
为了降低门槛，文章里的完整代码用归一化二维FFT作为稀疏变换，这样读者不需要额外安装工具包就能跑通流程。理解了整个框架之后，把变换算子替换成curvelet非常容易，只需要改两个函数句柄。
2.2 欠采样方式：随机掩码和规则掩码的差别
理论上压缩感知要求观测矩阵和稀疏表示基满足低互相干性，这意味着采样方式最好是随机的。地震数据重建里，“随机”体现在空间方向上随机地保留一部分地震道、丢掉其余道。一组掩码矩阵Mask，尺寸与完整数据相同，1表示该道该时刻的采样点被观测到，0表示缺失。
有一点常被忽略：时间方向不要做随机欠采样，只沿空间方向（道方向）欠采样。因为地震数据的时间方向始终是完整采样的，压缩感知利用的是空间维度的稀疏性来完成道间插值。如果对时间轴也随机抽，重建问题变得病态并且完全不必要。
实际操作中比纯均匀随机更好的是jittered采样——大致保持等间距的基础上叠加随机抖动。因为纯随机分布的炮检距产生的空道簇容易聚集，jittered方式让缺失道的空间分布更均匀，重建更稳定，也接近可控采样采集设计的实际做法。
2.3 迭代收缩阈值算法（ISTA）为什么适合这个场景
重建优化问题里的L1范数项让整体目标函数不可微，不能用简单梯度下降直接求解。迭代收缩阈值算法（ISTA）是最容易理解且稳定的一种求解方式：每一轮先对当前系数计算保真项梯度，做一个梯度下降更新，再执行一次软阈值操作来促进稀疏性。软阈值操作就是：
[ S_\tau(x) = \operatorname{sign}(x) \cdot \max(|x| - \tau, 0) ]
用生活化的类比，梯度下降确保重建结果不偏离实际观测，软阈值则负责压制掉那些不重要的微弱系数，让能量集中在少量大系数上。
为什么不提最常被推荐的OMP或者FISTA？OMP适合系数非常稀疏且原子间互相干性低的情况，地震数据在f-k域虽然相对稀疏，但并没有稀疏到只有几个系数，OMP的贪婪策略性能不稳定。FISTA虽然收敛更快，但引入了动量项，刚开始跑通代码时如果参数没配好，误差曲线会出现震荡，反而不利于理解。先用ISTA把基线版做出来，后续优化成FISTA只是多几行代码的事。
3. 完整处理链路：从带缺口的炮集数据到重建后的完整波场
3.1 数据准备：如何从SEG-Y或者文本矩阵进入算法
假设你已经从SEG-Y里读出了一炮共炮点道集，或者一条二维测线的cmp道集，数据组织是一个二维矩阵d，维度是时间采样点数nt乘以道数nx。下一步就是生成采样掩码：
% 假设 d 是 nt x nx 的二维地震记录 % bad_flag 可以是一个向量，标记哪些道缺失，1表示缺失，0表示保留 bad_flag = zeros(1, nx); bad_flag([3, 4, 5, 18:30, 101:105]) = 1; % 示意：孤立坏道+连续缺道 Mask = ones(size(d)); Mask(:, bad_flag == 1) = 0; d_obs = d .* Mask; % 缺失道全部置零一个常见误区是直接对缺失道的数值做零填充然后进算法。如果Mask置零后不记录哪些位置是缺失的，优化过程就会把零值当成真实观测去拟合，重建结果会被拉向零幅度。因此Mask必须原样传入算法框架，每一次迭代的残差计算都只针对Mask==1的位置进行约束。
3.2 归一化傅里叶变换：正则性比变换本身更重要的一段细节
Matlab原生的fft2和ifft2实际上是互逆关系，但它们并不是等距算子——fft2不会自动做归一化，直接拿fft2作为正交基用会导致梯度步长估计困难。标准的做法是做一个归一化处理，构造一对严格等距的正交变换算子：
N = numel(d); F = @(x) fft2(x) / sqrt(N); % 正变换 Finv = @(X) sqrt(N) * ifft2(X); % 逆变换这样定义后，F和Finv满足 Finv(F(x)) = x，且算子范数为1，迭代算法的步长可以安心地取1附近。这个小细节是代码能否稳定收敛的关键。我第一次实现的时候没做归一化，直接用fft2，步长无论怎么调都容易发散，后来把归一化补上，问题立刻消失。
3.3 重建主流程的伪代码视角
完整迭代流程可以概括为四步循环：
- 计算当前重建数据的保真残差，只在Mask为1的位置考虑；
- 把残差从数据域变换回稀疏系数域，得到梯度方向；
- 沿负梯度方向更新稀疏系数；
- 对更新后的系数做软阈值收缩，完成一次稀疏性约束。
注意这里的变量域转换：矩阵d是时空域数据，矩阵X是变换域系数。重建过程完全在系数域运作，每次迭代只通过一次正变换和一次逆变换在时空域和稀疏域之间来回跳跃。
初始化时可以把X设为零矩阵。迭代若干次后，由于保真项的作用，被观测道的位置上重建结果会越来越接近真实数据，缺失道则靠着空间方向上的稀疏约束一点点被“长出来”。
4. Matlab核心代码逐段拆解与调参经验
4.1 主函数代码：seismic_recon_ista
下面这个函数是完整可运行的ISTA压缩感知重建实现，把输入、输出、参数、迭代收敛判据全部封装在一起。直接粘贴到Matlab里保存为seismic_recon_ista.m即可调用。
function [d_rec, X, hist] = seismic_recon_ista(d_obs, Mask, lambda, mu, maxIter, tol) % seismic_recon_ista 基于压缩感知的二维地震数据缺道重建（ISTA版） % 输入： % d_obs : nt x nx 带缺失道的二维地震记录（缺失道请置零） % Mask : nt x nx 逻辑/数值矩阵，1表示观测道，0表示缺失道 % lambda : L1约束系数，控制稀疏性与保真度的平衡（可选） % mu : 步长，建议使用归一化FFT时取 0.9~1.0（可选） % maxIter : 最大迭代次数（可选，默认300） % tol : 收敛阈值，用相对变化量判断（可选，默认1e-6） % 输出： % d_rec : 重建后的完整数据 % X : 最终的f-k域系数，用于诊断稀疏性 % hist : 结构体，包含每轮迭代的残差与相对误差 [nt, nx] = size(d_obs); N = nt * nx; % 归一化傅里叶变换算子，保证正交等距 F = @(x) fft2(x) / sqrt(N); Finv = @(X) sqrt(N) * ifft2(X); if nargin < 3 || isempty(lambda) % 经验初值：取观测数据在稀疏域峰值幅度的1/100左右 temp = abs(F(d_obs)); lambda = 0.05 * max(temp(:)); end if nargin < 4 || isempty(mu) mu = 1.0; end if nargin < 5 || isempty(maxIter) maxIter = 300; end if nargin < 6 || isempty(tol) tol = 1e-6; end % 初始化：稀疏系数置零，对应时空域数据也置零 X = zeros(size(d_obs)); d_rec = zeros(size(d_obs)); hist.residual = zeros(1, maxIter); hist.relError = zeros(1, maxIter); for iter = 1:maxIter % 1. 保真项残差：只在观测位置计算 residual = Mask .* d_rec - d_obs; % 2. 梯度：残差变换回稀疏系数域 grad = F(residual); % 3. 系数域梯度下降更新 X_new = X - mu * grad; % 4. 软阈值收缩，推进稀疏性 thresh = mu * lambda; X_new = sign(X_new) .* max(abs(X_new) - thresh, 0); % 5. 逆变换回到数据域 d_rec_new = real(Finv(X_new)); % 6. 收敛判断 relChange = norm(d_rec_new - d_rec, 'fro') / (norm(d_rec, 'fro') + 1e-12); res = norm(Mask .* (d_rec_new - d_rec), 'fro') / sqrt(nnz(Mask)); X = X_new; d_rec = d_rec_new; hist.residual(iter) = res; hist.relError(iter) = relChange; if relChange < tol hist.residual = hist.residual(1:iter); hist.relError = hist.relError(1:iter); break; end end end4.2 合成数据快速验证脚本
光说不练假把式，给一个可以直接跑通的验证脚本。我用一个简单的三层层状模型正演合成炮集，然后随机抽走约30%的地震道，再用上面的函数重建，跑完就能看到缺失道被恢复出来。
% 生成一个简易合成地震记录：水平层状模型，三组反射同相轴 nt = 512; dt = 0.004; nx = 128; dx = 10; t = (0:nt-1) * dt; x = (0:nx-1) * dx; d = zeros(nt, nx); % 第一组同相轴：近水平 t1 = 0.5 + 0.0005 * (x / dx); d = d + ricker_wavelet(t, t1, 35); % 第二组同相轴：带倾角 t2 = 0.8 + 0.002 * (x / dx); d = d + 0.7 * ricker_wavelet(t, t2, 30); % 第三组同相轴：反向倾斜 t3 = 1.1 - 0.0015 * (x / dx); d = d + 0.4 * ricker_wavelet(t, t3, 40); % 随机抽道：大约抽掉30%的道 rng(2024); bad_flag = rand(1, nx) < 0.3; Mask = ones(size(d)); Mask(:, bad_flag) = 0; d_obs = d .* Mask; % 调用重建函数 [d_rec, X, hist] = seismic_recon_ista(d_obs, Mask, [], [], 400, 1e-7); % 重建质量定量评估 SNR_obs = snr_db(d, d_obs); SNR_rec = snr_db(d, d_rec); fprintf('观测数据SNR = %.2f dB\n', SNR_obs); fprintf('重建数据SNR = %.2f dB\n', SNR_rec);4.3 配套辅助函数
上面用到了ricker_wavelet和snr_db两个自编函数，也一并给出：
function w = ricker_wavelet(t, tc, f0) % 从中心时间tc附近截取Ricker子波，并映射到对应采样点 for k = 1:length(tc) for it = 1:length(t) dt_ = t(it) - tc(k); w(it, k) = (1 - 2 * (pi * f0 * dt_).^2) * exp(-(pi * f0 * dt_).^2); end end endfunction val = snr_db(ref, signal) % 计算信噪比，单位为dB；ref为无噪参考信号，signal为待评价信号 numerator = sum(abs(ref(:)).^2); denominator = sum(abs(ref(:) - signal(:)).^2); val = 10 * log10(numerator / (denominator + 1e-12)); end理论上这个ricker_wavelet写法有更简洁的版本，但上面这种循环形式可读性更好，每列对应一个同相轴的到时。实际跑的时候你会发现第二个同相轴和第三个因为倾角较大，重建难度会明显高于近水平的第一个同相轴，这也正是后面要讨论的“陡倾角先验不足”问题。
4.4 参数怎么调：lambda，mu，maxIter三者各管什么
三个关键参数里，lambda是对重建结果影响最大也最容易调错的一个。lambda太小，稀疏约束基本不起作用，缺失道的能量会被随机分配到大量系数上，重建结果看起来“脏”，像是没有去干净的噪声；lambda太大，所有系数都被压得过狠，重建结果会变得过平滑，弱反射同相轴直接消失，振幅信息失真。
一个比较省心的做法是先用稀疏域系数的峰值幅度做参照系。上文代码里lambda默认取的是观测数据稀疏系数的5%，这个取值在多数二维叠后工区上都能给出不错的结果。如果你处理的资料信噪比很低，建议把lambda增大一倍，因为强噪声本身在稀疏域也是分布的，需要更强的收缩才能压制。
mu的取值在归一化FFT框架下可以放松，取0.9至1.0都行。千万别取2以上，迭代会震荡甚至直接发散。maxIter不用设太大，300次对于二维数据已经足够；如果你把算子换成更复杂的curvelet，收敛速度会慢不少，此时需要适当放开迭代上限到800至1000次。
5. 重建结果怎么评价，以及几个必须绕开的坑
5.1 定量评价与定性评价结合
定量评价最常用的是SNR，这个指标在合成数据验证阶段是万能的——因为你知道真实地下的完整波场是什么。以我那个三层模型为例，抽掉30%道之后，观测数据SNR大约只有5至6dB，重建之后能提升到16至20dB，这个提升幅度在主观剖面上一眼就能看出来。但SNR并不是越高越好，有些极端参数设置会让重建结果过度偏向噪声模型，SNR反而虚高，所以还要结合f-k谱做评价。
定性看三样东西：重建剖面上缺失道位置的同相轴是否连续；断层边界和弱振幅同相轴是否被保留；f-k谱上是否还残留沿波数方向的零乱能量。第一样过关说明整体重建成功，第二样决定这个方法能不能用到AVO等保幅处理里，第三样判断参数是否还有潜力可挖。
5.2 坑一：陡倾角同相轴重建质量明显变差
这是压缩感知地震重建最典型的痛点。深层反射或者大倾角绕射波在f-k域的能量紧邻Nyquist边界，采样率不足的时候，这部分能量即使稀疏也容易被欠采样混淆。如果你在处理的时候发现近中倾角的同相轴恢复得不错，但陡倾角同相轴出现了锯齿状断层，那不是代码的问题，是信息量本身不够。
缓解方式有两种：一是把稀疏变换从f-k域换到curvelet域或者剪切波域，这类多尺度几何变换对方向性的刻画更精细；二是有条件的话做多道联合重建，把相邻几个炮集放在一个矩阵里同时重建，共享空间结构信息。切到curvelet之后陡倾角会好一些，但计算时间至少翻三到五倍。
5.3 坑二：收敛判据只看残差容易被骗
如果你只盯着保真残差下降就停下来，很可能会得到一个“观测位置拟合得不错、缺失位置还是空的”伪收敛结果。原因是ISTA的保真项只约束观测道，缺失道的系数只要有微微一点能量就会被收缩掉，两者互相拉扯，到后期残差几乎不再变化，但缺失道能量还在缓慢增长。
我一般两个指标同时看：一个是保真残差，一个是重建数据相邻两次迭代的相对变化量。后者下降到1e-6以下时才能说真正收敛。上面代码里用relChange作为收敛条件，配置文件里也建议把tol保持在这个量级。
5.4 坑三：边界效应
二维FFT隐含着周期性延拓假设，数据边界不连续时，重建结果会在左右边界出现明显的伪同相轴。处理脚本里直接对全工区做重建时边界效应最明显，对策有两个：一是重建前对道集做边缘衰减，用余弦斜坡把边界道幅度逐渐降到零，重建后再恢复边界幅度；二是在每一维都填充一定数量的零道，重建完成后裁掉，这个办法简单但会增加计算量。
5.5 坑四：纯随机掩码其实不是最优解
实际工区的缺失道分布很少是完全随机，通常是障碍物导致的连续缺口。连续缺口对重建是恶意场景，因为缺失道簇内部的相邻关系全部丢失，稀疏恢复难度大增。预处理阶段可以用jittered过采样思路重新组织数据：先对完整采集网格做轻度规则化，把连续大缺口拆分到多个重建子块里处理，每个子块内缺失比例尽量均匀。这个步骤虽小，对最终剖面质量的影响却非常显著。
6. 往实际工区数据扩展时要改哪几处，以及我个人的体会和建议
把这份代码从合成数据接回到真实工区数据时，有四处必须改动。第一处是数据读入方式，SEG-Y文件需要先从segy工具包或者fread逐道读取，得到nt*nx的矩阵后再进重建流程。第二处是数据规模，真实工区的单炮道集可能达到几千道，矩阵尺寸大了之后每次迭代的fft2计算量都不会小，建议先截取测线的一个扇段做测试，跑通参数后再铺到全道长跑。第三处是Mask的构造方式，真实数据里坏道分布要从数据质量标识里提取，不能像我演示代码那样直接随机赋值。第四处是振幅保真度要求，如果用这套结果做AVO属性分析，lambda要往小调，确保弱反射振幅不被过度收缩。
关于数值效率，还有一个很实用的小技巧：在跑大规模数据之前，先用降采样的低分辨率版本定参数。具体做法是把每条道隔几个采样点抽稀，快速跑一轮迭代，把lambda的合理范围摸出来，再切回全分辨率正式重建。这个技巧能省下大量试错时间。
我个人的感觉是，压缩感知重建和很多地球物理反演问题一样，90%的工作在数据准备和参数诊断，真正跑算法只占很小一部分。如果你第一次跑出来的结果不理想，先别急着改算法结构，画出重建前后数据的f-k谱对比，多半能一眼看出问题是出在稀疏性假设不成立，还是参数给偏了。当我第一次把连续缺失几十道的中段同相轴干净利落地重建出来时，那种满足感是很强的，这套东西也因此成了我日常规则化处理流程里保留的一件常用工具。
