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

王本锋,陆文凯,陈小宏,王志凯

(1.同济大学海洋与地球科学学院,上海200092;2.清华大学自动化系,智能技术与系统国家重点实验室,北京100084;3.中国石油大学(北京)油气资源与探测国家重点实验室,北京102249)

地震数据采集中,受障碍物、禁采区、海上拖缆羽状漂移以及经济因素的制约,观测的地震数据往往在空间方向具有不规则性,影响了后续自由表面多次波压制(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]:

(1)

1.2 3D地震数据插值重建方法

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

(2)

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

(3)

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

(4)

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

(5)

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

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

(6)

(7)

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,计算效率提高近一倍,验证了理论分析的正确性。

图1 模拟数据a 完整数据体; b 不规则缺失数据体(缺失50%)

2.2 实际地震数据

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

图2 模拟数据测试结果a,b 有效频率域POCS方法重建结果及重构误差; c,d 时间域POCS方法重建结果及重构误差

图3 模拟算例重构信噪比曲线

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

图4 实际资料a 完整数据体; b 不规则缺失数据体(缺失50%)

图5 实际资料重建结果a 有效频率域POCS方法; b 时间域POCS方法

图6 实际资料重构信噪比曲线

3 结论与认识

本文基于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

[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

[6] HERRMANN F J,HENNENFENT G.Non-parametric seismic data recovery with Curvelet frames[J].Geophysical Journal International,2008,173(1):233-248

[7] HENNENFENT G,HERRMANN F J.Simply denoise:wavefield reconstruction via Jittered undersampling[J].Geophysics,2008,73(3):V19-V28

[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

[9] WANG B F,LI J Y,CHEN X H.A novel method for simultaneous seismic data interpolation and noise removal based on the L0norm 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

[12] RONEN J.Wave-equation trace interpolation[J].Geophysics,1987,52(7):973-984

[13] FOMEL S.Seismic reflection data interpolation with differential offset and shot continuation[J].Geophysics,2003,68(2):733-744

[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

[16] MA J W.Three-dimensional irregular seismic data reconstruction via low-rank matrix completion[J].Geophysics,2013,78(5):V181-V192

[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

[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

[23] WANG B F,CHEN X H,WANG T.A novel reconstruction error definition for seismic data interpolation[J].Expanded Abstracts of 85thAnnual 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

[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

[28] BLUMENSATH T,DAVIES M E.Iterative hard thresholding for compressed sensing[J].Applied and Computational Harmonic Analysis,2009,27(3):265-274

[29] BLUMENSATH T,DAVIES M E.Iterative thresholding for sparse approximations[J].Journal of Fourier Analysis and Applications,2008,14(5):629-654