企业官网建设流程全解析

简介：面向地震数据重建与压缩感知算法学习者，这是一份轻量的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优化无法区分；随机欠采样则把混叠能量打散成近似均匀的随机噪声背景，重建算法在寻找稀疏解时可以自然抑制这些伪影。由此可以对比三种采样方式的表现。

纯随机采样在局部会出现大范围空道，导致中浅层覆盖次数波动大；抖动采样（Jittered Sampling）在规则网格上叠加有界随机位移，既保留整体均匀性又引入随机性。抖动幅度一般取网格距的0.5至1倍，过小退化为规则抽稀，过大接近纯随机采样，都有损重建效果。

2.3 L1范数优化：欠定方程的正确解结构

观测过程写为 y=Ax+n，A包含欠采样掩码与稀疏变换，n是观测噪声。由于丢道，方程欠定，最小二乘解会把能量分配到所有系数上，得到充满空间假频的波场。L1范数惩罚项 min ||x||₁，约束 ||Ax-y||₂≤ε，强制解稀疏，只保留下少数大系数，从而恢复出干净波场。目标函数是凸的，存在全局最优解且不依赖初值，这对批量处理地震数据非常重要。

ε的取值与噪声水平直接相关。残差阈值小于噪声水平时，算法会拟合噪声，重建结果出现颗粒状伪影；阈值过松则会丢掉弱信号。实际处理中通过交叉验证曲线标定ε，具体做法放到第四章参数校准部分展开。

3. 地震压缩感知的稀疏变换选型与观测矩阵搭建

3.1 四种变换的适用场景与参数选择

不同稀疏基对地震数据的表达能力差异很大，选错重建质量会明显下降。项目中常见的选择是FK、小波、Curvelet、Seislet四种变换。它们各自的特征对比如下。

选型经验是：信噪比高、构造平缓的资料优先FK变换，速度快且稳定；构造复杂或信噪比低的资料切到Curvelet；手头有可靠速度模型时试Seislet，它从预测属性出发，压缩率通常最高，但对前期处理质量最敏感。实际项目中为了控制计算成本，经常先用FK重建快速验证参数是否正确，再用Curvelet跑正式结果。

3.2 抖动采样掩码与复合观测算子的代码实现

观测矩阵的工程实现核心是丢道逻辑与随机扰动。下面的代码生成抖动采样掩码，适用于二维炮集数据。

import numpy as np def jitter_mask(nt, nx, ratio=0.4, jitter=3): """生成地震道抖动采样掩码 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, size=n_keep) idx = np.clip(base + offset, 0, nx - 1) idx = np.unique(idx) mask = np.zeros((nt, nx), dtype=bool) mask[:, idx] = True return mask # 使用示例 mask = jitter_mask(nt=512, nx=200, ratio=0.3, jitter=2) 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, dtype=np.complex128) data[self.mask] = y return self.transform.fwd(data)

逻辑说明：forward对应 y=RF⁻¹x，adjoint对应 x=F(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, lam=0.01, n_iter=200, tol=1e-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), dtype=complex) 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个点，分别重建并计算重建结果与原始数据的SNR，SNR曲线的拐点处就是合适的量级。这个扫描在单炮规模上跑，几分钟时间就能完成，再把选定的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%。

这种级联方式有两个直接好处：一是每层的正交变换更匹配窄频带信号，稀疏性更好，重建精度反而高于全频带一次解；二是每层数据体变小，内存和耗时下降，适合大炮集批处理。

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