压缩感知地震重建：从稀疏变换到FISTA的工程实践

简介面向地震数据重建与压缩感知算法学习者这是一份轻量的MATLAB示例包聚焦利用稀疏性与L1范数优化完成缺失地震道重建。压缩感知理论允许在低于奈奎斯特采样率下采样只要地震信号在傅里叶域或小波域满足稀疏条件即可通过求解优化问题恢复完整地震波形包中MATLAB脚本演示了从稀疏基构造、观测矩阵选取到正则化重建的核心过程运行后可直观查看压缩感知重建效果并可作为进一步修改实验的起点。整个压缩包仅含1个m文件总大小约2KB代码量精简依赖少适合研究生、工程师快速阅读和调试。目前已有213人学习/下载对刚接触地震数据压缩感知、希望从代码层面理解重建流程的读者具有直接参考价值也能为后续降低存储成本、提高成像质量的研究应用打下基础。1. 压缩感知地震重建用1/4数据换取完整波场地震勘探成本的大头压在外业采集上实际布设明显受地形、障碍物制约很难按理想网格采满道。高密度采集成本压力下用25%~40%的道数完成同等覆盖非常常见这时常规插值误差急剧放大。压缩感知地震重建要替代的正是这个环节只要波场在某个变换域可稀疏表示且观测过程与稀疏基满足非相干性重建就变成有唯一解的稀疏优化问题而不是对空道做局部拟合。零道距、坏道剔除、采集脚印消除、规则化重排放进CS框架是同一类问题。从 chongbianxie 这类压缩感知地震重建码包常见的落地过程看稀疏变换选型、欠采样设计、L1求解器配置、质量评价构成完整主线下文按这条主线逐步展开。2. 压缩感知的数学前提与地震波场的映射2.1 稀疏性判断用系数衰减测试做变换域取舍压缩感知成立的第一前提是信号在某一变换域具有稀疏性也就是大部分变换系数的模接近零。地震数据在时空域看起来信息量大转到频率-波数域FK域后则会显著变稀疏能量集中在与视速度对应的扇形区域内其余位置接近于零。对于含陡倾角、断层、盐丘的地区FK域的稀疏程度不够需要引入多尺度、多方向的Curvelet变换或Seislet变换才能用少数系数表达复杂波前。工程上判断变换是否够稀疏不靠主观印象跑一次系数衰减测试就有量化结果。以二维中值滤波后的叠前数据切片为例做正变换后把系数按模从大到小排序再统计累积能量曲线。前10%系数解释总能量85%以上时后续重建的稳定性和保幅性都有保障低于70%则说明基函数与数据形态不匹配需要换域。测试成本就是一次正变换加一次排序建议在项目启动时把候选变换逐一跑一遍用数据说话。2.2 非相干性随机采样为什么优于规则抽稀第二个条件是观测矩阵与稀疏基之间满足非相干性。直观来看规则抽稀在FK域产生规则的空间假频假频能量与真实信号在确定位置重叠L1优化无法区分随机欠采样则把混叠能量打散成近似均匀的随机噪声背景重建算法在寻找稀疏解时可以自然抑制这些伪影。由此可以对比三种采样方式的表现。采样方式变换域伪影重建难度施工复杂度规则抽稀规则假频与真实信号重叠高常规插值和CS都吃力低纯随机采样均匀随机噪声背景低CS重建效果好中空道概率大抖动采样弱规则性加随机扰动低工程上最常用低施工易实现纯随机采样在局部会出现大范围空道导致中浅层覆盖次数波动大抖动采样Jittered Sampling在规则网格上叠加有界随机位移既保留整体均匀性又引入随机性。抖动幅度一般取网格距的0.5至1倍过小退化为规则抽稀过大接近纯随机采样都有损重建效果。2.3 L1范数优化欠定方程的正确解结构观测过程写为 yAxnA包含欠采样掩码与稀疏变换n是观测噪声。由于丢道方程欠定最小二乘解会把能量分配到所有系数上得到充满空间假频的波场。L1范数惩罚项 min ||x||₁约束 ||Ax-y||₂≤ε强制解稀疏只保留下少数大系数从而恢复出干净波场。目标函数是凸的存在全局最优解且不依赖初值这对批量处理地震数据非常重要。ε的取值与噪声水平直接相关。残差阈值小于噪声水平时算法会拟合噪声重建结果出现颗粒状伪影阈值过松则会丢掉弱信号。实际处理中通过交叉验证曲线标定ε具体做法放到第四章参数校准部分展开。3. 地震压缩感知的稀疏变换选型与观测矩阵搭建3.1 四种变换的适用场景与参数选择不同稀疏基对地震数据的表达能力差异很大选错重建质量会明显下降。项目中常见的选择是FK、小波、Curvelet、Seislet四种变换。它们各自的特征对比如下。变换域稀疏表达特征适用数据形态计算复杂度主要限制FK变换能量集中在视速度锥水平缓倾角、高信噪比低对构造突变、空道敏感小波变换时频局部化强、方向性弱纵向突变、横向平缓低倾斜同相轴表达弱Curvelet多尺度多方向断层、陡倾角、复杂构造高内存占用大、参数多Seislet沿同相轴预测高信噪比、速度稳定中依赖叠加速度场质量选型经验是信噪比高、构造平缓的资料优先FK变换速度快且稳定构造复杂或信噪比低的资料切到Curvelet手头有可靠速度模型时试Seislet它从预测属性出发压缩率通常最高但对前期处理质量最敏感。实际项目中为了控制计算成本经常先用FK重建快速验证参数是否正确再用Curvelet跑正式结果。3.2 抖动采样掩码与复合观测算子的代码实现观测矩阵的工程实现核心是丢道逻辑与随机扰动。下面的代码生成抖动采样掩码适用于二维炮集数据。import numpy as np def jitter_mask(nt, nx, ratio0.4, jitter3): 生成地震道抖动采样掩码 Parameters ---------- nt : int 时间采样点数 nx : int 空间道数 ratio : float 保留道比例0.4 表示保留 40% 的道 jitter : int 随机抖动范围单位为道间距 Returns ------- mask : numpy.ndarray, (nt, nx), bool 保留位置为 True n_keep int(np.floor(nx * ratio)) base np.linspace(0, nx - 1, n_keep).astype(int) offset np.random.randint(-jitter, jitter 1, sizen_keep) idx np.clip(base offset, 0, nx - 1) idx np.unique(idx) mask np.zeros((nt, nx), dtypebool) mask[:, idx] True return mask # 使用示例 mask jitter_mask(nt512, nx200, ratio0.3, jitter2) obs full_data * mask # full_data 为完整观测数据逻辑说明base在空间方向均匀分布offset引入最大±jitter道间距的随机位移np.unique去除抖动后落在同一道号上的重复项。这样得到的掩码在大尺度上覆盖均匀在小尺度上具备随机性兼顾施工可行性与重建质量。真实项目中应按炮集的道头坐标生成掩码而不是按道序号操作否则观测几何与野外布设不匹配。获得掩码后还需要把稀疏变换和采样算子合成一个整体。二维数据的完整观测矩阵尺寸可能是百万乘百万量级显式存储不现实因此代码里统一用函数封装正向与伴随过程。class CSOperator: 压缩感知观测算子 A R o F^{-1} forward: x - 反变换到数据域 - 按掩码抽取 adjoint: y - 掩码填零 - 正变换到变换域 def __init__(self, transform, mask): self.transform transform # 提供 fwd / inv 接口 self.mask mask def forward(self, x): data self.transform.inv(x) return data[self.mask] def adjoint(self, y): data np.zeros_like(self.mask, dtypenp.complex128) data[self.mask] y return self.transform.fwd(data)逻辑说明forward对应 yRF⁻¹xadjoint对应 xF(Rᵀy)两个过程的调用必须配对梯度方向才不会反。transform接口只需提供fwd和inv两个方法FK、小波、Curvelet的实现都可以适配这样在对比不同稀疏域的重建效果时切换成本极低。4. 地震数据重建实现FISTA求解器与参数校准4.1 ISTA与FISTA的收敛特性差异L1求解器有一阶方法和内点法等选择地震数据规模大内存占用低的一阶方法更现实。ISTA每轮迭代做一次梯度下降加一次软阈值收缩收敛速度偏慢。FISTA在ISTA基础上引入Nesterov动量用前两轮解的差做外插把收敛率从O(1/k)提升到O(1/k²)同样迭代次数下残差下降幅度更明显。两种方法每轮计算成本几乎相同都只需一次正向算子和一次伴随算子因此工程默认优先用FISTA。个别场景下FISTA会因动量项产生振荡重建结果出现横向条纹。这时先检查梯度步长是否偏大其次检查lam是否过小若两者都正常仍振荡退回ISTA即可。从实际数据表现看FISTA在200轮内达到的精度ISTA通常要跑四五百轮以上才能接近。4.2 FISTA重建的完整可运行代码下面代码以FK域作为稀疏基封装一个可直接运行的FISTA重建流程。替换变换接口后可直接切换到Curvelet等变换。import numpy as np class FKTransform: FK域稀疏基正变换为二维FFT逆变换取实部 def fwd(self, data): return np.fft.fft2(data) def inv(self, coeff): return np.real(np.fft.ifft2(coeff)) def soft_threshold(x, thr): 软阈值算子对每个系数做收缩 return np.sign(x) * np.maximum(np.abs(x) - thr, 0) def fista_cs(y, mask, lam0.01, n_iter200, tol1e-6): FISTA 求解 min 0.5||Ax - y||^2 lam||x||_1 y : 欠采样观测数据 (nt, nx) mask : 采样掩码 (nt, nx) lam : 稀疏正则化系数 op CSOperator(FKTransform(), mask) x np.zeros_like(np.fft.fft2(y), dtypecomplex) z x.copy() t 1.0 L 2.0 # Lipschitz常数上界估计 for it in range(n_iter): grad op.adjoint(op.forward(z) - y) x_prev x x soft_threshold(z - grad / L, lam / L) t_new (1 np.sqrt(1 4 * t**2)) / 2 z x ((t - 1) / t_new) * (x - x_prev) t t_new rel np.linalg.norm(x - x_prev) / (np.linalg.norm(x_prev) 1e-12) if rel tol: break return op.transform.inv(x)参数说明lam控制稀疏惩罚强度和数据拟合项的平衡调大则结果干净但弱信号被削平调小则保幅更好但伪影增多。L是梯度Lipschitz常数取算子A的谱范数平方的上界使用归一化FFT时谱范数为1加观测掩码后上界不超过2。若迭代中间目标函数不降反升把L乘以2再继续。逻辑说明y是观测数据op.forward(z)把当前变换域系数先反变换到数据域再抽取减去真实观测得到残差op.adjoint把残差填零并正变换得到梯度方向的更新量。soft_threshold实现的是近端梯度步阈值lam/L对应L1项的近端映射。t_new按Nesterov公式更新让外插量逐步加大。4.3 正则化系数与停止准则的标定方法实际处理包里比如 chongbianxie.zip 这类压缩感知地震重建代码包lam的标定极少靠经验拍板推荐用L曲线扫描完成。取一块有代表性的叠前切片lam从1e-4到1e-1按对数均匀取8到10个点分别重建并计算重建结果与原始数据的SNRSNR曲线的拐点处就是合适的量级。这个扫描在单炮规模上跑几分钟时间就能完成再把选定的lam推广到全区批处理。对于信噪比偏低的数据拐点会变得不明显这时建议取拐点偏大一侧优先保证构造主导的波场被还原。停止准则方面除了代码中的相对变化量还可以记录目标函数值F(x)0.5||Ax-y||²lam||x||₁。FISTA的目标函数不是单调递减每20轮记录一次以相邻两次记录变化率小于1e-4作为收敛信号。实际项目中重建结果的视觉变化往往远早于数值收敛就趋于静止继续迭代主要是改善振幅精度对构造形态影响有限。因此建议设n_iter上限为300到500轮防止异常数据把单次运行时间拖到不可接受。提示FISTA目标函数不单调判断收敛时至少间隔20轮比较目标函数值单独盯相邻两次迭代的数值波动容易误判为已经收敛。5. 重建质量验证与级联重建的实用技巧5.1 用盲道SNR和频谱曲线看重建是否真实质量评价不能只看参与重建道的拟合误差否则过拟合会被误判为高精度。项目中的做法是先随机抽出5%的道充当盲测道不参与重建完成后对比这些道上的参考数据与重建结果。def snr_db(ref, rec): 计算信噪比单位 dB noise ref - rec return 10 * np.log10(np.sum(ref ** 2) / (np.sum(noise ** 2) 1e-12)) blind_snr snr_db(blind_ref, blind_recon)盲道SNR能说明算法把未观测信息恢复到了什么程度。对于信噪比大于20dB的资料重建后盲道SNR一般不应低于15dB低于10dB则要重点检查稀疏域选择或lam设定是否合理。同时对比重建前后的频谱曲线重点关注低频端和高频端是否有整体抬高或凹陷。频谱异常通常指向观测矩阵与稀疏基失配而不是迭代次数不够。5.2 复杂构造区重建的边界条件与失效信号断层、盐丘边界、陡倾角区是压缩感知重建最容易失守的地方。原因在于这些区域能量在变换域中扩散到多个尺度和方向稀疏性变差。常见应对是把曲线型同相轴所在的频带单独切出来重建并适当增大正则化系数以抑制变换域的旁瓣。如果构造过于复杂单靠CS重建难以还原断裂细节这时把速度模型信息作为约束加入目标函数比盲目调lam更有效。5.3 频率分层级联重建把单次CS求解拆成多个小尺度的地震重建问题全频带一次性重建在大工区上耗时很长且低频与高频的最优lam往往不一样。常用技巧是按频带分三层先重建低频部分约束大尺度构造形态再把低频重建结果作为初始值依次推高中频、高频部分。每层的lam按该频带能量动态设定通常取当前频带最大振幅的1%到5%。这种级联方式有两个直接好处一是每层的正交变换更匹配窄频带信号稀疏性更好重建精度反而高于全频带一次解二是每层数据体变小内存和耗时下降适合大炮集批处理。本文还有配套的精品资源点击获取