基于一种三维低冗余曲波变换和压缩感知理论的地震数据重建

曹静杰, 王尚旭, 李文斌

(1.中国石油大学地球物理与信息工程学院,北京 102249; 2.河北地质大学勘查技术与工程学院,河北石家庄 050031;3.河北地质大学信息工程学院,河北石家庄 050031)

基于一种三维低冗余曲波变换和压缩感知理论的地震数据重建

曹静杰1,2, 王尚旭1, 李文斌3

(1.中国石油大学地球物理与信息工程学院,北京 102249; 2.河北地质大学勘查技术与工程学院,河北石家庄 050031;3.河北地质大学信息工程学院,河北石家庄 050031)

基于稀疏变换的重建方法是地震数据重建的研究热点,其中稀疏变换的性质对重建的效率和质量起关键作用。曲波变换对波场数据有非常稀疏的表示和可靠的数值效果,然而三维曲波变换的冗余度在24～32之间,须消耗很多的内存和计算时间。提出基于一种低冗余度曲波变换的地震数据重建方法,介绍低冗余曲波变换并且分析其优点,给出一个解基于分析的1范数模型的快速迭代阈值方法。数据实验结果表明:该低冗余变换将三维曲波变换的冗余度降低了60%,能够有效提高数据重建的计算效率;基于低冗余度曲波变换的数据重建的计算效率是基于原始曲波变换重建效率的4倍,对于10%的采样比例仍然能够得到较好的重建和去噪效果。

曲波变换; 地震重建; 稀疏优化; 1范数

由于受到地表障碍、坏道、禁区和成本等因素的制约,野外地震数据在空间上经常不满足采样定理,不完备的地震数据会影响波动方程偏移[1]、表面多次波消除[2]和AVO分析[3]等数据处理的效果,因此须重建地震数据以获得完备的数据。基于稀疏变换的重建方法计算效率高,数值效果稳定,因此是重建方法重要的研究方向。该方法利用信号在变换域的稀疏性为约束建立反演模型,并且采用稀疏反演方法求解。稀疏变换的性质是决定这类方法质量的首要因素。地震信号在变换域中的表达越稀疏,则越有利于提高重建质量。在地震信号处理领域经典的稀疏变换有傅里叶变换[4]和Radon类变换[5-6]。最近曲波变换[7]等多尺度变换被用于地震数据重建且展示了很好的数值效果。多尺度变换比单一尺度的变换能够展示更多的信息,因此是信号分析的有力工具。曲波变换打破了传统傅里叶变换和小波变换只能表示单一尺度或者不能稀疏表示多尺度信号的缺点。对于地震数据有非常稀疏的表达,在地震数据重建[8-10]、去噪[11-12]、多次波消除[13-14]、反褶积[15]、偏移[16]和反射数据的增强[17]等问题中都有应用并且展示了很好的数值效果。曲波变换对曲线或曲面形状信号进行多方向和各项异性的表示,为边缘不连续的目标提供了几乎最优的表达[18-21]。虽然曲波变换有多尺度性、多方向性以及无须对数据划分时间窗口等优点,但是曲波变换的冗余度(冗余度指的是变换域系数的维数与原始信号系数维数之比)很大。在二维情况下的原始曲波变换的冗余度约为8,三维数据变换的冗余度为 24～32,需要很多的内存和计算时间,影响大规模数据计算的效率。笔者将一种低冗余三维曲波变换用于三维地震数据重建,这种新的低冗余曲波变换[22]比原始曲波变换的冗余度有很大的下降,因此称之为低冗余度曲波变换,在求解时采用一种解基于分析的1范数正则化模型的迭代软阈值算法,对模拟和真实地震数据的试验证明方法的有效性。

1 三维低冗余度曲波变换

曲波[23]变换有两种快速实现方式,分别是基于USFFT(unequally-space fast transform)的快速算法和基于wrapping的快速算法。由于基于wrapping的方法更加简单和容易执行,本文中采用基于wrapping的曲波变换的实现过程。基于wrapping的三维曲波[23]变换包括一个低通逼近子带和一些包含曲波系数的曲波基子带,这些曲波系数由尺度、位置和方向索引[22]。该变换主要由两部分组成:多尺度分解和角度分解。首先,假设输入的三维数据为N=(Nx,Ny,Nz),其中，Nx、Ny、Nz分别表示x、y、z三个方向上的维数,在傅里叶域通过三维Meyer小波将其分解,得到尺度为N,N/2,…,N/2J的以原点为中心的立方体,其中J为尺度的个数。第二步,将这些立方体和角度窗口相乘,分割成各向异性的楔形体,这些楔形体符合抛物线尺度准则,将楔形体包装成立方体后采用三维傅里叶逆变换得到曲波系数。三维低冗余度曲波变换的具体实现方式如下:①输入N=(Nx,Ny,Nz)的三维数据体,定义划分的尺度个数J,在最小尺度的每个面上划分的角度个数Na。 ②多尺度分解。在傅里叶域应用三维Meyer小波变换,得到尺度为N,N/2,…,N/2J的立方体。③角度分解。在每一个尺度的每个方向上,在傅里叶域将小波立方体和角度窗相乘得到楔形体;将这些楔形体包装成中心在原点的平行六面体;对平行六面体应用三维傅里叶逆变换得到曲波系数。

在原始曲波变换中,冗余度主要来源于角度分解。此外,Meyer小波的应用方法产生了额外冗余性,而低冗余曲波变换中Meyer小波没有产生额外冗余性[22]。 新的曲波变换和原始曲波变换在实现上存在许多不同之处,但是其能够降低额外冗余度的主要原因是Meyer小波变换的应用方法不同。

(1)

其中ν从0到1逐渐变化并且满足ν((x)+ν(1-x)=1。Meyer 尺度函数定义为

(2)

最小尺度的小波在傅里叶域的支撑为[-2/3,-1/6]∪[1/6,2/3], 因此超出了Shannon 带。原始曲波变换隐含的假设周期边界条件,计算一个周期信号的小波变换等价于将信号在周期小波基上分解。最小尺度小波超出Shannon 带的部分被其以ξ=1/2为中心沿着竖轴做镜像代替,数据的支撑是原始数据的 4/3倍,对于三维数据来说,产生的额外的冗余度为 (4/3)3。在低冗余度曲波变换中,尺度函数和小波函数的支撑首先被缩小为原来的3/4,然后为了保持单位划分的一致性,对最小尺度的小波进行修改来压制其下降的尾部,小波在[-1/2,-1/4]∪[1/4,1/2]之间变为一个固定值。这样对于三维Meyer小波来说,在傅里叶域内没有增加额外的冗余度。

对于原始曲波变换,三维数据的冗余度为 24～32,而低冗余度曲波变换的冗余度约为10[22]。除了在最小尺度上保持方向选择性和低冗余性,该低冗余曲波变换保持了等距性,并且能够快速精确重建,即该低冗余度曲波变换是一个Parseval 紧框架,C*C=I,其中C为曲波分解算子,C*为其伴随。C*也是逆算子,因此该变换具有快速的重构算法。关于低冗余度曲波变换的详细实现见参考文献[22]。

2 基于压缩感知理论的地震数据重建

地震采集过程是地震信号的采样过程,因此地震重建问题可以根据压缩感知理论求解[7]。压缩感知理论证明了可以根据少量的测量数据高精度地重建原始数据,从而降低地震勘探的成本。地震采集过程可以表示为如下数学模型:

Rx+ε=d.

(3)

式中,R为采样矩阵;x为完整的地震数据;ε为可加性噪声;d为采样的数据。由于数据d是不完备且不满足采样定理的,理论上存在无穷多x满足式(3)。根据压缩感知理论,假设地震数据在某个变换域中是稀疏的,即s=Ψx是稀疏的或者可压缩的,其中Ψ是一个稀疏变换算子,令Ψ为上述的低冗余曲波变换,由于其是紧框架,其共轭转置也是其逆,即Ψ*=Ψ-1,则公式(3)变为

RΨ*s+ε=d.

(4)

其中Ψ*为Ψ的共轭转置。基于s的稀疏性假设,可以通过解以下优化问题求出s:

(5)

其中σ≥0是对噪声能量的估计,该问题称为基追踪问题[24]。 当σ=0时,上述模型针对高信噪比数据或者不含噪声数据的重建;当σ>0时,上述模型可以求解含随机噪声数据的重建,即可以实现同时重建和去噪,因此上述模型既可以求解不含噪声数据的重建,也可以求解含随机噪声数据的重建。当s被解出后,对于低冗余曲波变换来说,x可以由x=Ψ*s求出。由于问题(5)是一个约束优化问题,对于该问题的求解方法很少[10],一般求解其等价的无约束优化问题

(6)

其中λ是一个正则参数,用来平衡残差的2范数和解的1范数,当λ选择合适时,问题(5)和(6)的解相等。为了使问题(6)的解与数据模型相匹配,对于不含噪声数据的重建,只须选择λ很小,接近于0即可;对于含随机噪声的重建,若σ已知,则可以通过偏差原则,广义偏差原则[25]等方法求解与σ匹配的正则参数λ,若σ已知,则可以采用类似于广义交叉检验原则(GCV)或者L曲线准则的方法求解合适的λ。由于广义交叉检验原则(GCV)或者L曲线准则不能保证一定获得可靠的λ[26],因此λ常常通过人工经验选取来获得。

地震数据重建属于大规模计算问题,因此需要快速稳定的计算方法。对于上述1范数正则化的问题(6)的求解,不能采用经典的方法如共轭梯度法或者牛顿法[27-28],传统的内点法等由于用到矩阵求逆,因此计算速度无法与经典方法相比,须研究解1范数正则化模型的快速算法。阈值类方法[29]是针对1范数正则化模型求解的一类重要方法,该方法根据1范数可分的性质得到,很多解问题(6)的方法,如不动点法、逼近点法等都含有软阈值运算。最近研究发现采用如下的模型的数值效果比基于模型(6)的效果要好:

(7)

该模型称为基于分析的模型,以时间空间域的数据x为变量,问题(6)被称为基于合成的模型,以变换域中的系数s为变量。虽然两者之间可以相互推导,但是基于分析的模型(7)的数值效果要稍好于前者。给出一个解问题(7)的简单易行的阈值算法,该算法以x为变量进行迭代求解。算法如下:

(1)输入采样的数据d,采样矩阵R,稀疏变换Ψ,最大迭代步数N,控制停机的参数ω和最小阈值τmin,令k=0。

(5)输出最终解xfinal=xk。

如果将算法中软阈值运算变为硬阈值或者其他改进的阈值运算,上述算法也能得到稀疏解。本文中步长κk固定为1,同时还有更加先进的步长选择方法如BB法、回溯法等方法,步长的选择对算法有较大的影响。阈值的下降可以是线性的,也可以是指数下降,本文中采用指数下降τk+1=τ0ec(k-1)/(N-1),其中c=ln(τ0/τmin),实际上阈值的下降方式对解的改善不大,最小阈值τmin的大小决定着解的质量,因为τmin与噪声的能量相关。当噪声很小时,选择τmin很小,当噪声较大时τmin应适当调节。由于噪声的能量实际上也是很难估计的,本文中采用尝试法选择τmin。本文中算法与固定阈值的阈值类算法不同,实际上是带有连续化策略的阈值法,其中令阈值逐渐下降就是优化算法中的连续化策略,这比固定阈值参数的算法需要更少的迭代次数和更好的数值效果,因此本算法实际上由两层组成,外层控制阈值的下降,内层只用1次迭代,已经证明本算法比固定最小阈值的方法数值表现要改进很多[26]。

3 数值模拟

图1 三维合成数据重建实验Fig.1 3D synthetic data reconstruction experiment

参数采样率50%CandesWoiselle采样率20%CandesWoiselle采样率10%CandesWoiselle采样率5%CandesWoiselle计算时间t/s10661767561907591761082271信噪比/dB203365202013135222171638876601036314116576232

用一个海洋地震数据验证基于低冗余曲波变换的真实数据的重建质量。该海洋数据为512×128×128,即时间轴上的采样个数为512,时间采样率2 ms,在inline和crossline 方向上的采样个数为128,相邻地震道间隔为25 m。图2给出海洋数据的一个切片,图3是随机采样的数据,采样的道数为总道数的50%。图4是采用低冗余曲波变换和迭代软阈值法的重建结果,其信噪比为 16.559 6 dB。为了更加清晰的显示重建的效果,将其中一道原始数据和其相应的重建道进行对比,如图5所示,重建数据和原始数据吻合得很好。

图2 真实海洋数据Fig.2 Field marine data

图3 随机采样道数为50%的海洋数据Fig.3 Randomly sampled marine data with 50% traces sampled

图4 基于低冗余曲波变换的海洋数据重建结果Fig.4 Reconstructed marine data based on low-redundancy curvelet transform

图5 一道原始海洋数据与重建数据对比Fig.5 Comparison of one original trace and reconstructed trace of marine data

为了验证不同采样比例下该低冗余曲波变换的重建效果,对一个含噪声的三维海洋数据分别采样20%和10%的地震道进行重建,该试验同时也验证了低冗余曲波变换对含噪声数据重建的可靠性。由于计算机内存的限制,截取一小部分数据进行计算,数据为200×200×200,即在时间上有200个采样点,空间两个方向上分别有200个采样,时间采样率为4 ms,道间距为12.5 m。 图6(a)为原始数据的一个切片显示,图6(b)为随机采样 20%地震道的显示,图6(c)为图6(b)采用低冗余曲波变换得到的重建结果。由图可见,不仅未采样的地震道被重建出来,而且噪声的能量得到了压制。 图6(d)为原始数据和重建结果的误差,视觉上有效信号很少。为了测试更低采样率的数据能否得到较好的重建,将采样率降至10%,即随机选择10%的地震道进行重建,如图7(b)所示。图7(c)为基于低冗余曲波变换的重建结果,图7(d)为原始数据和重建数据的差。由此可见,图7(c)的结果是可以接受的,并且图7(d)中有效信号很少。为了在频率中比较,图8(a)为原始数据图7(a)的一个空间切片的频谱,图8(b)为重建数据(图7(c))相同空间切片的频谱,图8(c)为图7(d)相同空间切片的频谱,可以看出基于低冗余曲波变换的重建可以滤掉高波数分量和一些高频分量。本试验说明:①低冗余曲波变换对于非常低的采样比例的数据也能够得到稳定的结果,当采样比例为5%时,数值效果会变差。②通过选择合适的参数,采用的模型和算法能够实现含噪声数据的重建,即问题(6)和(7)的模型可以解决含随机噪声的重建问题。③相比于二维变换,三维变换需要的采样比例更少。之前的二维数据重建研究表明,当以曲波变换为稀疏变换,并且采用阈值类算法如凸集投影法[26]求解时,二维数据的采样比例在50%左右才能得到较好的结果,一旦低于这个比例,则数值效果会变差,三维变换由于将两个空间维数的信息整合起来,具有更稀疏的表达,因此能够适合更低的采样比例。同样,五维数据重建也适合更低的采样比例。

图6 三维真实数据随机采样20%的重建试验Fig.6 3D field data reconstruction with 20% traces sampled randomly

图7 三维真实数据随机采样10%的重建Fig.7 3D field data reconstruction with 10% traces sampled randomly

图8 真实数据和重建数据的频谱比较Fig.8 Spectrum comparation of original data and reconstructed data

4 结束语

将一种低冗余度曲波变换用于加速地震数据重建,给出了一种快速的迭代阈值法实现地震数据重建和去噪,提高了基于曲波变换的地震数据重建的计算效率和重建质量。数据试验表明，基于低冗余度曲波变换的重建时间约为基于原始曲波变换的1/4,且能够达到可靠的重建结果,可以适合更低采样比例的数据重建。基于多尺度变换的重建存在一定的局限性,由于曲波变换是局部的、多尺度变换,如果没有对小尺度构造充分的采样,则基于这种多尺度变换的重建不能得到满意的效果。

[1] LIU B, SACCHI M. Minimum weighted norm interpolation of seismic records [J]. Geophysics,2004,69(6):1560-1568.

[2] NAGHIZADEH M. Parametric reconstruction of multidimensional seismic records [D]. Edmonton:University of Alberta,2009.

[3] SACCHI M, LIU B. Minimum weighted norm wavefield reconstruction for AVA imaging [J]. Geophysical Prospecting, 2005,53(6):787-801.

[4] ZWARTJES P, GISOLF A. Fourier reconstruction of marine-streamer data in four spatial coordinates [J]. Geophysics,2006,71(6):V171-V186.

[5] SACCHI M, VERSCHUUR D, ZWARTJES P. Data reconstruction by generalized deconvolution [C/OL]//2004 SEG Annual Meeting. Denver,Colorado, October 10-15,2004:SEG Technical Program Expanded Abstracts, 2004:1989-1992 [2009-03-05].http://library.seg.org/doi/abs/10.1190/1.1843303.

[6] 薛亚茹,杨静,钱步仁. 基于高阶稀疏Radon变换的预测多次波自适应相减方法[J]. 中国石油大学学报(自然科学版),2015,39(1):43-49.

XUE Yaru, YANG Jing, QIAN Buren. Multiples prediction and adaptive subtraction with high-order sparse Radon transform [J]. Journal of China University of Petroleum (Edition of Natural Science),2015,39(1):43-49.

[7] HERRMANN F, MOGHADDAM P, STOLK C. Sparsity-and continuity-promoting seismic image recovery with curvelet frames[J]. Applied and Computational Harmonic Analysis,2008,24:50-173.

[8] 曹静杰,王彦飞,杨长春. 地震数据压缩重构的正则化与零范数稀疏最优化方法[J]. 地球物理学报,2012,55(2):596-607.

CAO Jingjie, WANG Yanfei, YANG Changchun. Seismic data restoration based on compressive sensing using the regularization and zero-norm sparse optimization[J]. Chinese Journal of Geophysics (in Chinese),2012,55(2):596-607.

[9] YANG P, GAO J, CHEN W. Curvelet-based POCS interpolation of nonuniformly sampled seismic records[J]. Journal of Applied Geophysics,2012,79(1):90-99.

[10] 曹静杰,王本锋. 基于一种改进凸集投影方法的地震数据同时插值和去噪[J]. 地球物理学报,2015,58(8):2935-2947.

CAO Jingjie, WANG Benfeng. An improved projection onto convex sets method for simultaneous interpolation and denoising[J]. Chinese Journal of Geophysics (in Chinese),2015,58(8):2935-2947.

[11] HENNENFENT G, HERRMANN F. Seismic denoising with non-uniformly sampled curvelets [J]. Computing in Science & Engineering,2006,8(3):16-25.

[12] NEELAMANI R, BAUMSTEIN A, GILLARD D, et al. Coherent and random noise attenuation using the curvelet transform[J]. The Leading Edge,2008,27(2):240-248.

[13] HERRMANN F, WANG D, HENNENFENT G, et al. Curvelet-based seismic data processing:a multiscale and nonlinear approach[J]. Geophysics, 2008,73(1): A1-A5.

[14] LIN T, HERRMANN F. Robust estimation of primaries by sparse inversion via one-norm minimization[J]. Geophysics,2013,78(3):R133-R150.

[15] KUMAR V, HERRMANN F. Deconvolution with curvelet-domain sparsity[C/OL]//2008 SEG Annual Meeting. Las Vegas,Nevada, November 9-14,2008:SEG Technical Program Expanded Abstracts, c2008:1996-2000 [2011-05-07].http://library.seg.org/doi/abs/10.1190/1.3059287.

[16] CHAURIS H, NGUYEN T. Seismic demigration/migration in the curvelet domain[J]. Geophysics, 2008,73(2):S35-S46.

[17] KUMAR V, OUEITY J,CLOWES R M, et al. Enhancing crustal reflection data through curvelet denoising[J]. Tectonophysics,2011,508(1):106-116.

[18] CANDES E, DONOHO D. Curvelets: a surprisingly effective nonadaptive representation for objects with edges[R].Stanford Univ Ca Dept of Statistics, 2000.

[19] STARKC J, CANDES E, DONOHO D. The curvelet transform for image denoising[J]. IEEE Transaction on Image Processing,2002,11(6):670-684.

[20] CANDES E, DONOHO D. New tight frames of curvelets and optimal representation of objects with piecewise C2 singularities[J]. Communications on Pure and Applied Mathematics,2004,57(2):219-266.

[21] CANDES E, DEMANET L, DONOHO, et al. Fast discrete curvelet transforms[J]. Multiscale Modeling & Simulation,2006,5(3):861-899.

[22] WOISELLE A, STARCK J, FADILI J. 3D Data denoising and inpainting with the low-redundancy fast curvelet transform[J]. Journal of Mathematical Imaging and Vision,2011,39(2):121-139.

[23] YING L, DEMANET L, CANDES E. curvelet transform[C/OL].//Proc. SPIE wavelets XI, San Diego. 2005, vol. 591413 [2008-9-4]. https://www.spiedigitallibrary.org/conference-proceedings-of-spie/5914/1/3D-discrete-curvelet-transform/10.1117/12.616205.short?SSO=1.

[24] CHEN S, DONOHO D, SAUNDERS M. Atomic decomposition by basis pursuit[J]. SIAM Journal on Scientific Computing,1998,20(1):33-61.

[25] 肖庭延,于慎根,王彦飞. 反问题的数值解法[M]. 北京:科学出版社,2003.

[26] YAGOLA A, LEONOV A, TITARENKO V. Data errors and an error estimation for ill-posed problems[J]. Inverse Problems in Science and Engineering,2002,10(2):117-129.

[27] CAO J, WANG Y, WANG B. Accelerating seismic interpolation with a gradient projection method based on tight frame property of curvelet[J]. Exploration Geophysics,2015,46(3):253-260.

[28] 王彦飞. 反问题的计算方法及其应用[M]. 北京:高等教育出版社,2007.

[29] CAO J, WANG Y, ZHAO J, et al. A review on restoration of seismic wavefields based on regularization and compressive sensing [J]. Inverse Problems in Science and Engineering,2011,19(5):679-704.

(编辑 修荣荣)

Seismicreconstructionwith3Dlow-redundancycurvelettransformandcompressedsensingtheory

CAO Jingjie1,2, WANG Shangxu1, LI Wenbin3

(1.CollegeofGeophysicsandInformationEngineeringinChinaUniversityofPetroleum,Beijing102249,China; 2.CollegeofExplorationTechnologyandEngineering,HebeiGEOUniversity,Shijiazhuang050031,China; 3.CollegeofInformationEngineering,HebeiGEOUniversity,Shijiazhuang050031,China)

Sparse-transform-based seismic data reconstruction is a hot topic in seismic reconstruction, where properties of the sparse transform may influence the results of reconstruction greatly. Curvelet transform is a multi-scale, multi-directional, and local transform which has nearly the sparsest expression for seismic data. However, this transform is a highly redundant transform with redundancy about 24-32 for three dimensional data. To improve the efficiency of curvelet based reconstruction, this paper proposed a low-redundancy curvelet-transform based seismic reconstruction. The new transform was introduced first and its merits for seismic signal processing were analyzed, followed by an iterative thresholding method for analysis-based L1-norm regularized models. Numerical experiments illustrate that the low-redundancy transform can reduce 60% redundancy of the original 3D curvelet transform, thus improves greatly the computational efficiency. The reconstruction computational efficiency based on the low redundancy transform is 4 times of the original curvelet based reconstruction, for example, even for 10% sampling ratio, this low-redundancy curvelet can get acceptable results.

curvelet transform; seismic reconstruction; sparse optimization; one-norm

2016-05-23

国家自然科学基金项目(41674114);河北省自然科学基金项目(D2017403027);河北省高校百名优秀创新人才支持计划Ⅲ(SLRC2017024);中国博士后科学基金资助项目(2016M600171，2017T100137)

曹静杰(1982-)，男，副教授，博士，硕士生导师，研究方向为地震信号处理和地球物理反演算法。E-mail: cao18601861@163.com。

1673-5005(2017)05-0061-08

10.3969/j.issn.1673-5005.2017.05.007

P 631.4

A

曹静杰,王尚旭,李文斌.基于一种三维低冗余曲波变换和压缩感知理论的地震数据重建[J].中国石油大学学报(自然科学版),2017,41(5):61-68.

CAO Jingjie, WANG Shangxu, LI Wenbin. Seismic reconstruction with 3D low-redundancy curvelet transform and compressed sensing theory [J]. Journal of China University of Petroleum(Edition of Natural Science), 2017,41(5):61-68.