基于3D Curvelet变换的频率域高效地震数据插值方法研究

地震数据采集中, 受障碍物、禁采区、海上拖缆羽状漂移以及经济因素的制约, 观测的地震数据往往在空间方向具有不规则性, 影响了后续自由表面多次波压制(SRME)、全波形反演以及逆时偏移成像精度。另外, 处理阶段对废炮、废道的剔除也是造成数据不规则的重要因素之一。为了提高地震数据的横向连续性, 进而提高后续多道处理以及精细油气藏描述的精度, 有必要进行地震数据的插值重建[1-7]。

地震数据插值重建方法根据实现方式可以分为:基于数学变换的插值方法[4-5, 8-10]、基于预测滤波的插值方法[11]、基于波动方程的插值方法[12-13]以及近几年来发展起来的基于降秩的插值方法[14-17]。其中基于数学变换的插值方法相对简单且容易实现, 随着压缩感知理论的发展得到了广泛应用。常用的稀疏变换包括Fourier变换、Curvelet变换、Dreamlet变换以及Seislet变换等[18-22]。其中Curvelet变换是一种新兴的稀疏变换, 其相对于小波变换增加了角度的信息, 更加适用于曲线奇异性的表征, 在地震数据插值重建以及去噪处理中得到了广泛的应用。基于稀疏变换建立目标泛函, 借助优化方法求解泛函, 最终得到插值重建结果。在众多优化求解方法中, 最流行的是凸集投影方法(projection onto convex sets, POCS), 因为它具有容易理解、方便实现的优点。POCS方法实际上是一种两步法:首先利用迭代硬阈值方法得到变换域的解, 再将其投影到观测平面上[9-10, 23]。由于Curvelet变换的冗余度较大(当最细尺度为Curvelet时, 2D Curvelet变换的冗余度为8~10, 3D Curvelet变换的冗余度为24~32)[22], 计算效率较低, 因此目前仅在2D地震资料处理中得到应用。张华等[24]及王本锋等[1]采用2D Curvelet变换对每一时间切片或频率切片进行插值, 最终实现了3D地震数据的插值重建, 但是该方法没有利用时间或频率的连续性约束。基于3D Curvelet变换的插值重建方法可以充分利用3D数据各个方向的连续性约束, 但是计算量较大, 必须有效提高计算效率。基于地震信号频谱的共轭对称性, 仅对有效频率数据体进行插值重建, 计算效率可提高一倍以上[5]。实际上, 影响插值重建计算效率的因素很多, 例如稀疏变换种类[4, 25-26]、插值重建方法[27]、迭代误差设计[23]以及数据体的合理表征[5]等。本文主要关注数据体的合理表征, 采用规模较小的数据体表征原数据, 通过对规模较小的数据体进行插值重建, 有效提高插值重建的计算效率。

3D Curvelet变换的数学表达式如下[22]:

式中:f为原始信号; ϕ为曲波基; 〈·, ·〉表示内积; 为ϕ的共轭函数; j, l, k分别表示模型域(即Curvelet域)对应的尺度、角度、位置参数; z∈R3表示数据域位置参数, 即二维空间、一维时间(或频率)的函数; x(j, l, k)为三维数据体对应的Curvelet系数。详细的数学理论以及推导参阅文献[22]。

不规则数据体与规则完整数据体的关系可由公式(2)近似表征:

式中:dobs为观测的3D不规则数据体, R为设计的采样算子, d0为规则的完整3D数据体。地震数据插值重建的目的是从观测数据以及采样算子中恢复规则完整数据体。受观测数据频带有限等因素的影响, 方程(2)的求解是不适定的。考虑到完整数据体d0在Curvelet变换域的稀疏性, 基于压缩感知理论, 采用稀疏促进策略, 构建目标泛函Φ(·)如下:

式中:C为3D Curvelet变换, CT为3D Curvelet逆变换; λ为正则化因子, 用以权衡拟合残差以及Curvelet系数稀疏性。公式(3)可由POCS优化方法[4, 5, 8-10]迭代求解:

式中:dk是第k步的更新解; I为与采样算子R同维数的单位矩阵; Tλk是硬阈值函数[28-29], λk是阈值, 可由改进的指数阈值模型[1]确定。

式中:λmax, λmin为设置的最大、最小阈值; N为最大迭代次数; k为当前迭代次数。基于迭代公式(4), 迭代次数达到设定最大迭代次数后, 便可以得到插值重建后的地震数据。

由于观测的地震信号为时间信号, 对时间信号进行一维Fourier变换得到的频谱具有共轭对称性质, 因此负频率成分可以由正频率成分进行精确估计。此外, 由于观测信号的频带有限, 利用有效频带内的频率成分即可对原始信号进行高精度的表征。该策略已在2D地震数据的插值重建中得到应用[5], 本文将该策略推广应用到3D地震数据的插值重建中。

假设观测数据dobs的维数为nx×ny×nt, 沿着时间轴对其进行Fourier变换, 得到有效频率数据体dfobs, 其维数为nx×ny×nf, 其中nf≤nt/2。对于n×n×n规模的数据体, 3D Curvelet变换的计算量[22]为O(n3lg(n)), 因此对规模减半的数据体进行插值重建, 计算效率可以提高一倍。在有效频率域, 观测数据与完整数据之间的关系为:

式中:dfobs, df0为观测数据以及完整数据体dobs和d0对应的有效频率数据体; Rf是与R对应的有效频率域采样算子。将3D Curvelet变换作用于频率域数据体, 假设其Curvelet系数也是稀疏的[5], 根据公式(6), 基于3D Curvelet变换的POCS迭代插值重建公式为:

式中:dkf是第k步的更新解, If为与采样算子Rf同维数的单位矩阵。迭代到最大迭代次数后, 便可得到有效频率域的完整数据体, 借助Fourier逆变换, 便可以得到时间域完整的3D数据体。

用模拟得到的3D数据对上述方法进行了测试, 该数据共有128炮, 每炮128道, 每道128个时间采样点; 炮间距以及道间距均为24m, 时间采样率为8ms。图 1a为某一固定排列的三维数据体的平面展示(图中①为等时间切片, ②为共炮点道集, ③为共检波点道集), 其随机缺失50%的地震道集如图 1b所示, 可以看出, 地震数据的横向连续性遭到破坏, 这将影响后续多道处理的精度。

对不规则缺失数据体进行频谱分析后, 选取最大有效频率为奈奎斯特频率, 即fmax=fN=62.5Hz, 利用本文提出的有效频率域POCS方法以及常用的时间域POCS方法对不规则数据体进行插值重建。设最大迭代次数为50, 重建结果以及重构误差如图 2所示。由图 2a和图 2c可以看出, 两种方法均可以对缺失地震数据进行较好的重建, 且重构误差(图 2b, 图 2d)在允许的误差范围内。通过重构信噪比进一步比较了两种方法的异同。图 3为两种方法重构信噪比曲线, 可见随着迭代次数的增加, 两种方法均趋于收敛, 收敛时的信噪比分别为24.0dB和26.3dB, 有效频率域POCS方法精度略低。将重建后的有效频率数据体转换到时间域, 最终的重构信噪比为22.1dB。有效频率域POCS方法的精度相对时间域POCS方法的精度稍低的原因之一是时间采样率较大, 当时间采样率较小时, 有效频率域POCS方法与时间域POCS方法的精度相当, 甚至略好[5], 具体的影响因素及其量化分析留作后续研究。在计算效率上, 时间域POCS方法所用时间为2918s, 有效频率域POCS方法所用时间为1512s, 计算效率提高近一倍, 验证了理论分析的正确性。

从实际海洋拖缆数据中截取一小块3D数据体对本文方法的有效性进行了测试。该数据体共120炮, 每炮120道, 每道151个时间采样点; 炮间距及道间距均为25m, 时间采样率为8ms。图 4a为该数据体平面展示图(图中①为等时间切片, ②为共炮点道集, ③为共检波点道集), 其随机缺失50%的地震道集如图 4b所示, 可以看出, 缺失后地震数据的横向连续性变差。

对实际数据进行频谱分析, 根据有效频带分布范围, 选择fmax=fN=62.5Hz, 利用本文提出的有效频率域POCS方法及时间域POCS方法对图 4b中不规则数据体进行插值重建, 最大迭代次数根据经验设置为50。图 5a为本文方法重建结果, 图 5b为时间域POCS方法重建结果, 两种方法重建结果与完整数据体的一致性均较好。图 6展示了两种方法重构信噪比曲线, 可见随着迭代次数的增加, 两种方法均趋于收敛。将重建后的有效频率数据体转换到时间域, 最终的重构信噪比为17.0dB, 与时间域方法的重构信噪比17.3dB相近。但是, 有效频率域POCS方法所用时间为1719s, 时间域POCS方法所用时间为3396s, 即有效频率域POCS方法的计算效率相对时间域POCS方法提高了近一倍, 与理论分析结果一致, 进一步验证了本文方法的有效性。当减小最大有效频率时, 得到的有效频率数据体规模将减小, 可以提高插值重建的计算效率, 但是以牺牲插值重建精度为代价, 具体的量化分析将留作后续研究。

本文基于3D Curvelet稀疏变换, 利用频率域地震数据的共轭对称性质, 得到规模减半的有效频率数据体, 应用凸集投影(POCS)方法, 研究了高效的3D地震数据插值重建方法。研究结果表明:有效频率域POCS方法是一种有效的插值重建手段; 对规模减半的有效频率数据体进行插值重建, 能在保证插值精度的同时, 提高约一倍的计算效率。

本文仅仅围绕地震数据插值重建问题展开研究, 关于插值去噪一体化以及考虑衰减信息的插值重建, 将是我们下一步的研究计划。由于Curvelet变换的计算效率较低, 因此高效的稀疏变换以及稀疏字典的研究是我们另一个研究方向。