2. 清华大学自动化系, 智能技术与系统国家重点实验室, 北京 100084;
3. 中国石油大学(北京)油气资源与探测国家重点实验室, 北京 102249
2. State Key Laboratory of Intelligent Technology, Department of Automation, Tsinghua University, Beijing 100084, China;
3. State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249, China
地震数据采集中, 受障碍物、禁采区、海上拖缆羽状漂移以及经济因素的制约, 观测的地震数据往往在空间方向具有不规则性, 影响了后续自由表面多次波压制(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]等。本文主要关注数据体的合理表征, 采用规模较小的数据体表征原数据, 通过对规模较小的数据体进行插值重建, 有效提高插值重建的计算效率。
1 方法原理 1.1 3D Curvelet变换3D Curvelet变换的数学表达式如下[22]:
$ x\left( {j,l,k} \right) = \left\langle {f,{\phi _{j,l,k}}} \right\rangle = \int_{{{\bf{R}}^3}} {\mathit{\boldsymbol{f}}\left( z \right){{\bar \phi }_{j,l,k}}\left( z \right){\rm{d}}z} $ | (1) |
式中:f为原始信号; ϕ为曲波基; 〈·, ·〉表示内积; 为ϕ的共轭函数; j, l, k分别表示模型域(即Curvelet域)对应的尺度、角度、位置参数; z∈R3表示数据域位置参数, 即二维空间、一维时间(或频率)的函数; x(j, l, k)为三维数据体对应的Curvelet系数。详细的数学理论以及推导参阅文献[22]。
1.2 3D地震数据插值重建方法不规则数据体与规则完整数据体的关系可由公式(2)近似表征:
$ {\mathit{\boldsymbol{d}}_{{\rm{obs}}}} = \mathit{\boldsymbol{R}}{\mathit{\boldsymbol{d}}_0} $ | (2) |
式中:dobs为观测的3D不规则数据体, R为设计的采样算子, d0为规则的完整3D数据体。地震数据插值重建的目的是从观测数据以及采样算子中恢复规则完整数据体。受观测数据频带有限等因素的影响, 方程(2)的求解是不适定的。考虑到完整数据体d0在Curvelet变换域的稀疏性, 基于压缩感知理论, 采用稀疏促进策略, 构建目标泛函Φ(·)如下:
$ \mathit{\Phi }(x) = \left\| {{\mathit{\boldsymbol{d}}_{{\rm{obs}}}} - \mathit{\boldsymbol{R}}{\mathit{\boldsymbol{C}}^{\rm{T}}}\mathit{\boldsymbol{x}}} \right\|\;_2^2 + \lambda \left\| \mathit{\boldsymbol{x}} \right\|{\;_0} $ | (3) |
式中:C为3D Curvelet变换, CT为3D Curvelet逆变换; λ为正则化因子, 用以权衡拟合残差以及Curvelet系数稀疏性。公式(3)可由POCS优化方法[4, 5, 8-10]迭代求解:
$ {\mathit{\boldsymbol{d}}_k} = {\mathit{\boldsymbol{d}}_{{\rm{obs}}}} + \left( {\mathit{\boldsymbol{I - R}}} \right){\mathit{\boldsymbol{C}}^{\rm{T}}}{\mathit{\boldsymbol{T}}_{{\lambda _k}}}(\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{d}}_{k - 1}}) $ | (4) |
式中:dk是第k步的更新解; I为与采样算子R同维数的单位矩阵; Tλk是硬阈值函数[28-29], λk是阈值, 可由改进的指数阈值模型[1]确定。
$ \begin{array}{l} {\lambda _k} = {\lambda _{{\rm{max}}}}{{\rm{e}}^{c\sqrt {\left( {k - 1} \right)/(N - 1)} }}\\ c = {\rm{ln}}({\lambda _{{\rm{min}}}}/{\lambda _{{\rm{max}}}})\;\;\;\;k = 1, 2, \ldots, N \end{array} $ | (5) |
式中:λ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)), 因此对规模减半的数据体进行插值重建, 计算效率可以提高一倍。在有效频率域, 观测数据与完整数据之间的关系为:
$ \mathit{\boldsymbol{d}}_{{\rm{obs}}}^f = {\mathit{\boldsymbol{R}}^f}\mathit{\boldsymbol{d}}_0^f $ | (6) |
式中:dfobs, df0为观测数据以及完整数据体dobs和d0对应的有效频率数据体; Rf是与R对应的有效频率域采样算子。将3D Curvelet变换作用于频率域数据体, 假设其Curvelet系数也是稀疏的[5], 根据公式(6), 基于3D Curvelet变换的POCS迭代插值重建公式为:
$ \mathit{\boldsymbol{d}}_k^f = \mathit{\boldsymbol{d}}_{{\rm{obs}}}^f + ({\mathit{\boldsymbol{I}}^f} - {\mathit{\boldsymbol{R}}^f}){\mathit{\boldsymbol{C}}^{\rm{T}}}{\mathit{\boldsymbol{T}}_{{\lambda _k}}}(\mathit{\boldsymbol{Cd}}_{k - 1}^f) $ | (7) |
式中:dkf是第k步的更新解, If为与采样算子Rf同维数的单位矩阵。迭代到最大迭代次数后, 便可得到有效频率域的完整数据体, 借助Fourier逆变换, 便可以得到时间域完整的3D数据体。
2 数值算例 2.1 模拟数据用模拟得到的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变换的计算效率较低, 因此高效的稀疏变换以及稀疏字典的研究是我们另一个研究方向。
[1] |
王本锋, 陈小宏, 李景叶, 等. POCS联合改进的Jitter采样理论曲波域地震数据重建[J].
石油地球物理勘探, 2015, 50(1): 20-28 WANG B F, CHEN X H, LI J Y, et al. Seismic data reconstruction based on POCS and improved Jittered sampling in the Curvelet domain[J]. Oil Geophysical Prospecting, 2015, 50(1): 20-28 |
[2] |
王本锋. 基于反演理论的地震信号处理方法研究[D]. 北京: 中国石油大学(北京), 2015
WANG B F. Research on seismic data processing technology by inverse theory[D]. Beijing: China University of Petroleum, 2015 |
[3] |
曹静杰, 王本锋. 基于一种改进凸集投影方法的地震数据同时插值和去噪[J].
地球物理学报, 2015, 58(8): 2935-2947 CAO J J, WANG B F. An improved projection onto convex sets method for simultaneous interpolation and denoising[J]. Chinese Journal of Geophysics, 2015, 58(8): 2935-2947 DOI:10.6038/cjg20150826 |
[4] | WANG B F, WU R S, GENG Y, et al. Dreamlet-based interpolation using POCS method[J]. Journal of Applied Geophysics, 2014, 109(10): 256-265 |
[5] | WANG B F. An efficient POCS interpolation method in the frequency-space domain[J]. IEEE Geoscience and Remote Sensing Letters, 2016, 13(9): 1384-1387 DOI:10.1109/LGRS.2016.2589260 |
[6] | HERRMANN F J, HENNENFENT G. Non-parametric seismic data recovery with Curvelet frames[J]. Geophysical Journal International, 2008, 173(1): 233-248 DOI:10.1111/gji.2008.173.issue-1 |
[7] | HENNENFENT G, HERRMANN F J. Simply denoise:wavefield reconstruction via Jittered undersampling[J]. Geophysics, 2008, 73(3): V19-V28 DOI:10.1190/1.2841038 |
[8] | WANG B F, CHEN X H, LI J Y, et al. An improved weighted projection onto convex sets method for seismic data interpolation and denoising[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2016, 9(1): 228-235 DOI:10.1109/JSTARS.2015.2496374 |
[9] | WANG B F, LI J Y, CHEN X H. A novel method for simultaneous seismic data interpolation and noise removal based on the L0 norm constraint[J]. Journal of Seismic Exploration, 2015, 24(2): 187-204 |
[10] | WANG B F, WU R S, CHEN X H, et al. Simultaneous seismic data interpolation and denoising with a new adaptive method based on dreamlet transform[J]. Geophysical Journal International, 2015, 201(2): 1180-1192 |
[11] | SPITZ S. Seismic trace interpolation in the F-X domain[J]. Geophysics, 1991, 56(6): 785-794 DOI:10.1190/1.1443096 |
[12] | RONEN J. Wave-equation trace interpolation[J]. Geophysics, 1987, 52(7): 973-984 DOI:10.1190/1.1442366 |
[13] | FOMEL S. Seismic reflection data interpolation with differential offset and shot continuation[J]. Geophysics, 2003, 68(2): 733-744 DOI:10.1190/1.1567243 |
[14] | GAO J J, Sacchi D M, Chen X H. A fast reduced-rank interpolation method for prestack seismic volumes that depend on four spatial dimensions[J]. Geophysics, 2013, 78(1): V21-V30 DOI:10.1190/geo2012-0038.1 |
[15] | LÓPEZ O, KUMAR R, YILMAZ Ö, et al. Off-the-grid low-rank matrix recovery and seismic data reconstruction[J]. IEEE Journal of Selected Topics in Signal Processing, 2016, 10(4): 658-671 DOI:10.1109/JSTSP.2016.2555482 |
[16] | MA J W. Three-dimensional irregular seismic data reconstruction via low-rank matrix completion[J]. Geophysics, 2013, 78(5): V181-V192 DOI:10.1190/geo2012-0465.1 |
[17] |
马继涛, 王建花, 刘国昌. 基于频率域奇异值分解的地震数据插值去噪方法研究[J].
石油物探, 2016, 55(2): 205-213 MA J T, WANG J H, LIU G C. Seismic data noise attenuation and interpolation using singular value decomposition in frequency domain[J]. Geophysical Prospecting for Petroleum, 2016, 55(2): 205-213 |
[18] |
高建军, 陈小宏, 李景叶, 等. 基于POCS方法指数阈值模型的不规则地震数据重建[J].
应用地球物理, 2010, 7(3): 229-238 GAO J J, CHEN X H, LI J Y, et al. Irregular seismic data reconstruction based on exponential threshold model of POCS method[J]. Applied Geophysics, 2010, 7(3): 229-238 |
[19] | GAN S W, WANG S D, CHEN Y K, et al. Compressive sensing for seismic data reconstruction via fast projection onto convex sets based on seislet transform[J]. Journal of Applied Geophysics, 2016, 130(7): 194-208 |
[20] | GAN S W, WANG S D, CHEN Y K, et al. Dealiased seismic data interpolation using seislet transform with low-frequency constraint[J]. IEEE Geoscience and Remote Sensing Letters, 2015, 12(10): 2150-2154 DOI:10.1109/LGRS.2015.2453119 |
[21] | LIU W, CAO S Y, GAN S W, et al. One-step slope estimation for dealiased seismic data reconstruction via iterative seislet thresholding[J]. IEEE Geoscience and Remote Sensing Letters, 2016, 13(10): 1462-1466 DOI:10.1109/LGRS.2016.2591939 |
[22] | CANDÈS E, DEMANET L, DONOHO D, et al. Fast discrete Curvelet transforms[J]. Multiscale Modeling & Simulation, 2006, 5(3): 861-899 |
[23] | WANG B F, CHEN X H, WANG T. A novel reconstruction error definition for seismic data interpolation[J]. Expanded Abstracts of 85th Annual Internat SEG Mtg, 2015: 3875-3879 |
[24] |
张华, 陈小宏. 基于Jitter采样和曲波变换的三维地震数据重建[J].
地球物理学报, 2013, 56(5): 1637-1649 ZHANG H, CHEN X H. Seismic data reconstruction based on Jittered sampling and Curvelet transform[J]. Chinese Journal of Geophysics, 2013, 56(5): 1637-1649 DOI:10.6038/cjg20130521 |
[25] | CAO J J, ZHAO J T. Simultaneous seismic interpolation and denoising based on sparse inversion with a 3D low redundancy Curvelet transform[J]. Exploration Geophysics, 2017, 48(4): 422-429 |
[26] | CAO J J, ZHAO J T. 3D seismic interpolation with a low redundancy, fast Curvelet transform[J]. Journal of Seismic Exploration, 2015, 24(2): 121-134 |
[27] | ZU S H, ZHOU H, CHEN Y K, et al. Interpolating big gaps using inversion with slope constraint[J]. IEEE Geoscience and Remote Sensing Letters, 2016, 13(9): 1369-1373 DOI:10.1109/LGRS.2016.2587301 |
[28] | BLUMENSATH T, DAVIES M E. Iterative hard thresholding for compressed sensing[J]. Applied and Computational Harmonic Analysis, 2009, 27(3): 265-274 DOI:10.1016/j.acha.2009.04.002 |
[29] | BLUMENSATH T, DAVIES M E. Iterative thresholding for sparse approximations[J]. Journal of Fourier Analysis and Applications, 2008, 14(5): 629-654 |