做地震数据处理的同行应该都有过这种经历野外采集回来的一条测线因为过沟、过村庄、设备故障或者遇上坏道最终叠前数据里总是缺那么几道。缺道少道看着是小问题但到了偏移成像或者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或者FISTAOMP适合系数非常稀疏且原子间互相干性低的情况地震数据在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必须原样传入算法框架每一次迭代的残差计算都只针对Mask1的位置进行约束。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) % 计算信噪比单位为dBref为无噪参考信号signal为待评价信号 numerator sum(abs(ref(:)).^2); denominator sum(abs(ref(:) - signal(:)).^2); val 10 * log10(numerator / (denominator 1e-12)); end理论上这个ricker_wavelet写法有更简洁的版本但上面这种循环形式可读性更好每列对应一个同相轴的到时。实际跑的时候你会发现第二个同相轴和第三个因为倾角较大重建难度会明显高于近水平的第一个同相轴这也正是后面要讨论的“陡倾角先验不足”问题。4.4 参数怎么调lambdamumaxIter三者各管什么三个关键参数里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谱对比多半能一眼看出问题是出在稀疏性假设不成立还是参数给偏了。当我第一次把连续缺失几十道的中段同相轴干净利落地重建出来时那种满足感是很强的这套东西也因此成了我日常规则化处理流程里保留的一件常用工具。