三维曲波变换在地震资料去噪处理中的应用研究

张之涵,孙成禹,姚永强,肖广锐

(1.中国石油大学(华东)地球科学与技术学院,山东青岛266580;2.中海石油(中国)有限公司天津分公司渤海石油研究院,天津300452)

随着油气勘探进程的不断推进,地震勘探逐渐从浅层、平原地带向深层以及地表复杂地域转移,复杂工区的噪声干扰现象越来越严重。噪声压制的效果将直接影响后续的地震资料处理、解释工作,因此,地球物理工作者越来越重视去噪方法的研究。

压制地震噪声的方法有很多,如多项式拟合、中值滤波、KL变换、F-K滤波和奇异值分解等[1]。其中有一种思路是将地震数据进行稀疏表示,在变换域内进行阈值处理,这就要求稀疏表示所采用的基函数能够捕捉到地震波前,用较少的系数就可以表示地震波的主要特征。同时,复杂的地震数据局部特征变化明显,因此,也要求所采用的稀疏表示基具有良好的局部识别能力[2]。可用来进行地震数据稀疏表示的方法有很多,比较有代表性的包括傅里叶变换、小波变换、Surfacelet变换、曲波(Curvelet)变换等[3]。不过,地震波在地下传播时其波前面是一曲面,傅里叶变换和小波变换都不能很好地对其进行稀疏表示。近些年兴起的曲波变换是一种新的信号分析和图像处理方法,它由曲线状的基元构成,在保持传统小波多尺度性的同时,还具有多方向性和各向异性,更加适合表示地震波前特征,因此,可以对地震数据进行理想的稀疏表示。1999年,Candès等[4]在其博士论文中提出在脊波(Ridgelet)变换的基础上实现连续曲波变换,也就是第一代曲波变换。2002年,Candès等[5]提出了新的Curvelet框架体系,称之为第二代曲波变换,实现过程比第一代曲波变换简单、快速。2004年,Herrmann等[6]最先将曲波变换引入地震数据处理领域,成功实现了曲波域多次波衰减。国内外许多学者在利用曲波变换进行地震资料的噪声压制方面做了很多工作[7-12],将二维曲波变换分别应用于叠前地震资料、微地震资料的噪声压制以及面波压制中等。2005年,Candès等[13]首先提出了可以将二维曲波变换拓展到多维;随后,Ying等[14]在Candès等[13]的基础上提出了三维曲波变换的实现方法,并将其应用于图像处理当中。Dubois等[15]将三维曲波变换应用到动态图像处理中,取得了很好的效果。

目前,地球物理工作者一般都利用二维曲波变换进行地震资料去噪处理,但对于三维地震资料而言,二维曲波变换并没有考虑数据中第三维的信息,对整体的去噪效果造成影响。为此,我们将三维曲波变换引入到地震数据去噪处理中。对三维地震数据进行三维曲波正变换,在曲波域根据有效信号和噪声的特征差异,采用改进的非线性阈值方法进行处理,再进行三维曲波反变换,最终完成整个处理流程。

1 曲波变换基本理论

曲波变换结合了小波变换的多尺度特点和脊波变换的各向异性特点,能够对信号进行理想的稀疏表示,在系数项个数相同的情况下,曲波变换的误差最小。因此,曲波变换更适合用于信号分析。

1.1 二维曲波变换

(1)

式中:Rθ表示以θ为弧度的旋转。曲波系数可以简单地表示为两个元素f∈L2(R2)和φj,l,k的内积,因为连续曲波变换在频率域进行,可以得到频率域的曲波变换为

(2)

二维曲波变换对沿着光滑曲线上有奇异值的数据有着高效的表达,但是在处理三维数据时,却没有充分利用三维数据体内所包含的丰富信息。三维曲波变换针对沿着二维曲面上有边缘奇异值的数据,给出了最佳的描述方法,具有敏感的方向特性以及各向异性,能更好地捕捉二维变换所缺失的整体信息。为此,我们引入三维曲波变换。

1.2 三维曲波变换

三维曲波变换保持了二维曲波变换的性质[13-14],它以一个拥有抛物线尺度性质的楔形体作为其频率域支撑。在算法上,三维曲波变换与二维曲波变换相似,可以表示为两个元素f∈L2(R3)和φj,l,k的内积,即

(3)

(4)

(5)

其中,

Φj(ω1,ω2,ω3)=φ(2-jω1)·φ(2-jω2)·φ(2-jω3)

(6)

图1 三维曲波时频域表示[13]

2 三维曲波变换去噪原理

2.1 基于三维曲波变换进行地震资料去噪处理的可行性分析

前文提到,在地震资料去噪处理中,可以将地震数据进行稀疏表示,在变换域内进行阈值处理。在实际工作中,我们一般不直接在物理空间域进行去噪处理,而是先对原始数据进行某种变换,再对其在该变换域内的系数做阈值处理。主要原理是利用有效信号和随机噪声在变换域内,因特征差异而被分为不同大小的系数,再将小于某个设定阈值的所有系数做置零处理以去除噪声。因此,要达到理想的去噪效果,所使用的变换方法要能较稀疏地表示地震数据,即只保留很少的系数就可以理想地表示原始数据的主要特征,而大量的小系数可以直接滤除,且这些系数对主要特征影响不大[2]。

稀疏表示理论被广泛地应用于信号和图像处理等领域,其思想是用信号的稀疏性进行数学建模,通过少量非零的表示系数反映信号的内在结构和本质属性,得到信号在某个变换域上最简洁的一种表达方式[16]。它认为信号可以以低于奈奎斯特采样频率的频率进行采样,且仍能精确重建原始信号。地震信号作为信号的一种,由于其本身的一些特性使得一些变换能对其进行稀疏表达,这就为在地震数据处理中引入稀疏表示理论打下了坚实的基础。在地球物理领域,稀疏表示的理论经历了傅里叶变换、小波变换、脊波变换和曲波变换等阶段[17]。

曲波变换属于稀疏函数的表示范畴,具有多尺度以及良好的方位特性,它能够对数据边缘部分的直线或曲线提供更加稀疏的表示,从而降低了运算量并更好地保护数据边缘的特征。跟小波变换等不同,曲波系数不仅跟位置(空间域)和尺度(频率域)有关,还和方向有关。在地震勘探领域,曲波变换能将地震信号的能量集中到曲波系数上,而随机噪声的能量变换后仍为噪声,并且幅度相同,这样我们就能更加清楚地分辨出有效信号和随机噪声。此外,曲波变换不会使数据的边缘产生模糊,很适合处理地震数据。因此,在地震数据处理中,选择合适的阈值,对地震数据的曲波系数进行阈值处理,就可以达到去噪的目的。

2.2 三维曲波变换阈值法去噪原理

目前业内一般都利用二维曲波变换进行地震资料去噪处理,但对于三维地震资料而言,二维曲波变换并没有考虑第三维的信息,这会对横向上同相轴的连续性及整体的去噪效果造成影响。从理论上讲,处理三维地震资料应采用三维曲波变换,但长久以来,制约三维曲波变换应用的主要原因是其巨大的计算机内存占用量。随着计算机技术的发展,三维曲波变换处理三维数据变得可行。因此在本文中,我们引入三维曲波变换对三维地震资料进行处理。

设某三维地震信号S(x,y,z)=F(x,y,z)+N(x,y,z),其中，F(x,y,z)表示有效信号,N(x,y,z)表示随机噪声。对其进行三维曲波变换,得到不同尺度、不同方向的曲波系数,记为C。在三维曲波域,有效信号沿不同角度具有较强的相关性,其能量主要集中在有限的尺度和方向分量上;而随机噪声不具有相关性,其能量分布在各个尺度和方向分量上。我们利用曲波系数的相关性来区分有效信号与随机噪声。对两个不同尺度方向的曲波系数进行相关运算,得到相关系数R,并进行归一化处理。

(7)

式中:C1(j,l,k)和C2(j,l,k)代表两个不同尺度方向的曲波系数;R′为归一化处理后的相关系数。比较R′与C2(j,l,k),若R′≥C2(j,l,k),则认为该系数代表有效信号;否则,则认为代表随机噪声。

利用阈值法对这些曲波系数进行处理,就可以达到去噪的目的。阈值法去噪处理的关键在于阈值和阈值函数的选取,这将直接影响地震资料的去噪效果以及重构后地震资料的保真程度。

1995年,Donoho[18]提出了一种固定阈值估计方法,即

(8)

式中:σ为噪声标准差,大小为曲波系数中值与经验值0.6745的比值;N为地震信号的长度。此阈值估计方法会损伤曲波系数。2013年,张博[17]提出了一种新的阈值估计方法,并通过实验证实了其有效性。我们在文献[17]的基础上,结合三维曲波变换的性质,提出一种改进的阈值估计方法,即

(9)

式中:e代表与尺度方向有关的经验参数;σ代表噪声标准差,与传统的固定阈值不同,这里σ是最高尺度层曲波系数标准差的估计值;σc代表曲波域内所用曲波系数的标准差估计值。

在阈值函数的选取上,除采用常用的软阈值法和硬阈值法外,还选择了1995年Gao等[19]提出的一种半软阈值函数,表达式为

(10)

式中:λ1和λ2分别为估计的两个阈值。为了尽可能地保护有效信号不受损害,设定λ1=λ=eσσc,λ2=1.5λ=1.5eσσc。

2.3 三维曲波变换阈值法去噪处理的实现步骤

综上所述,提出的基于三维曲波变换的地震资料去噪方法的实现步骤如下。

1) 选择合适的尺度和方向,对三维地震记录进行三维曲波正变换,得到相应的曲波系数c(i,j,k)。

2) 利用公式(7)判别曲波系数,若R′≥C2(j,l,k),认为该系数代表有效信号,否则认为该系数代表随机噪声。

3) 利用公式(9)进行阈值λ的估计。

4) 在曲波域,选择合适的阈值函数对曲波系数进行阈值处理。将小于阈值的曲波系数看作是对噪声的表示,置零;大于阈值的曲波系数看作是对有效信号的表示,利用阈值函数进行处理。

5) 对经过阈值处理后的曲波系数进行重构,即进行三维曲波逆变换,得到重构后的三维地震数据即为去噪后的地震记录。

3 模型试算

通过模型试算验证方法的正确性。建立一个水平层状模型,如图2所示,通过正演得到三维地震模拟记录,炮点置于地面中心处。在模型记录中加入不同强度的噪声,使得地震记录的信噪比分别大于1,等于1和小于1。利用提出的三维曲波变换阈值法进行去噪处理,使用了硬阈值、软阈值和半软阈值3种函数并比较它们的效果。

图2 水平层状模型

当模型记录的信噪比大于1时,去噪效果如图3 所示。同时,为了进一步验证三维曲波变换的去噪效果,抽取炮点所在处的二维剖面,如图4所示。

经硬阈值去噪后记录的信噪比为3.08,经软阈值去噪后记录的信噪比为4.50,经半软阈值去噪后记录的信噪比为4.30。从三维图像以及二维剖面中可以看出,曲波阈值法在信噪比较高的情况下,能较好地去除噪声,效果较好,并且有效信号的损失很小。在曲率较大的地方,有效保留了有效信号。在同相轴交叉的地方,同相轴形态没有发生改变。对比可见,半软阈值和软阈值方法的去噪效果要好于硬阈值去噪方法。

图3 信噪比大于1时三维记录的去噪效果a 加噪记录; b 硬阈值去噪记录; c 软阈值去噪记录; d 半软阈值去噪记录

图4 信噪比大于1时二维剖面的去噪效果

改变模拟记录的信噪比,当信噪比等于1时,三维记录的去噪效果如图5所示。同时,为了进一步验证三维曲波变换的去噪效果,抽取炮点所在处的二维剖面,如图6所示。

经硬阈值去噪后记录的信噪比为2.01,经软阈值去噪后记录的信噪比为3.02,经半软阈值去噪后记录的信噪比为2.50。从三维图像以及二维剖面中可以看出,曲波阈值法在信噪比等于1的情况下,仍能较好地去除噪声,并且有效信号的损失很小。在同相轴弯曲以及同相轴相交的部位,较好地保留了有效信号。3种阈值方法相比,软阈值处理效果好于硬阈值,而半软阈值方法在提高信噪比、保留有效信号等方面,均能取得良好的效果。

图5 信噪比等于1时三维记录的去噪效果a 加噪记录; b 硬阈值去噪记录; c 软阈值去噪记录; d 半软阈值去噪记录

图6 信噪比等于1时二维剖面的去噪效果

为检验该方法对低信噪比地震记录的去噪效果,降低模拟记录的信噪比。当信噪比为0.50时,去噪效果如图7所示。同时,抽取炮点所在处的二维剖面,如图8所示。

经硬阈值去噪后记录的信噪比为1.40,经软阈值去噪后记录的信噪比为2.36,经半软阈值去噪后记录的信噪比为1.87。从三维图像以及二维剖面中可以看出,在信噪比小于1的情况下,很难从地震剖面上识别出连续同相轴,但是采用曲波去噪,能较好地去除噪声。在同相轴弯曲的地方,有效信号有一定程度的损失。3种阈值方法相比,硬阈值在曲率比较大的地方能更好地保留有效信号;软阈值处理能更多地去除随机噪声,但会对有效信号造成损害;半软阈值方法既能很好地去除噪声、又能保证同相轴特别是弯曲同相轴的信息不受损害。

通过模型试算可以看到,经过三维曲波变换处理后,随机噪声得到了压制,同时没有损害有效同相轴的信息,在同相轴曲率较大的部分能达到较好的效果。当地震信息中含有弱信号时,即使地震记录的信噪比较低,还是可以提取出有效弱信号。

在去噪过程中,阈值的选取十分重要,文中使用的3种阈值方法各有特点。硬阈值法保留弱信号的能力强,在曲率较大的地方去噪效果好;在相同阈值情况下,软阈值法去除噪声的效果好于硬阈值,但是会损害弯曲同相轴的信息;半软阈值法能较好地保留有效同相轴的信息,同时达到很好的去噪效果。因此,针对不同类型的数据可以选择不同的阈值进行处理:对于水平同相轴,可以选取半软阈值法或软阈值法进行处理;当资料中含有弱的弯曲同相轴时,采用半软阈值法或硬阈值法进行处理，以更好地保留弱同相轴信息。

图7 信噪比小于1时三维记录的去噪效果a 加噪记录; b 硬阈值去噪记录; c 软阈值去噪记录; d 半软阈值去噪记录

图8 信噪比小于1时二维剖面的去噪效果

4 实际资料测试

为进一步测试三维曲波变换在实际地震资料中的应用效果,将其应用于实际三维地震资料去噪处理。资料的in-line方向为250道,cross-line方向为200道,每道600个采样点。我们选取半软阈值函数进行处理。

图9为实际三维地震资料去噪前、后效果对比。从图9中可以看出,去噪处理后的同相轴变得连续、清晰,一些湮没在噪声里的同相轴显现出来,浅层的强同相轴没受影响,地震记录的信噪比提高。为了更为直观的显示,抽取二维剖面如图10所示。从图10中可以看出,浅层的同相轴明显变得清晰,0.2～0.4s的同相轴从噪声中恢复,同时有效同相轴保存得较为完好,资料的信噪比得到了明显的提高。为了定量显示三维曲波变换的去噪效果,我们计算测线上每个二维剖面的信噪比并将其连成曲线(图11)。从图11可以看出,经过三维曲波变换后的实际资料,信噪比得到了明显提高。

图9 实际三维地震资料去噪前(a)、后(b)效果对比

图10 实际三维地震资料抽取二维剖面去噪前(a)、后(b)效果对比

图11 实际三维地震资料去噪前、后各测线信噪比

5 结论与认识

我们引入三维曲波变换,结合改进的非线性阈值法,提出了一种基于三维曲波变换的地震资料去噪方法。通过理论模型和实际资料的去噪处理试验,得出以下结论与认识。

1) 在三维曲波域中,通过相关计算的方法判别曲波系数,使用改进的非线性阈值法对曲波系数进行处理,取得了良好的去噪效果。该方法不仅能够很好地压制地震资料中的随机噪声,而且能够很好地保护有效信号。模型试算和实际资料测试验证了此方法的正确性和可行性。

2) 3种阈值方法的试算结果表明,硬阈值法保留弱信号的能力强,在曲率较大处能取得较好的去噪效果;在相同阈值情况下,软阈值法压制噪声的能力优于硬阈值法,但会损害同相轴的信息;改进的半软阈值法能够在保护有效同相轴的基础上有效压制随机噪声。

3) 基于三维曲波变换的地震资料去噪方法对计算机的存储能力和算法的计算效率提出了更高的要求,接下来需要进行进一步的研究工作。

三维曲波变换具有良好的空间-频率-角度特性,能够很好地对地震数据进行稀疏表示,使得其在地震数据处理中具有更大的应用前景。

参 考 文 献

[1] 张军华.地震资料去噪方法——原理、算法、编程及应用[M].东营:中国石油大学出版社,2011:1-355

Zhang J H.Seismic data denoising method:principle,algorithm,programming and application[M].Dongying:China University of Petroleum Press,2011:1-355

[2] 唐刚.基于压缩感知和稀疏表示的地震数据重建与去噪[D].北京:清华大学,2010

Tang G.Seismic data reconstruction and denoising based on compressive sensing and sparse representation[D].Beijing:Tsinghua University,2010

[3] Shan H,Ma J W,Yang H Z.Comparisons of wavelets,contourlets and curvelets in seismic denoising[J].Journal of Applied Geophysics,2009,69(2):103-115

[4] Candès E J,Donoho D L.Curvelets:a surprisingly effective nonadaptive representation for objects with edges[M].Nashville:Vanderbilt University Press,1999:1-16

[5] Candès E J,Donoho D L.New tight frames of curevelets and optimal representations of objects with smooth singularities[R].USA:Department of Statistics,Stanford University,2002

[6] Herrmann F J,Verschuur E.Curvelet-domain multiple elimination with sparseness constraints[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004,1333-1336

[7] Herrmann F J,Wang D,Hennenfent G,et al.Curvelet-based seismic data processing:a multiscale and nonlinear approach[J].Geophysics,2007,73(1):A1-A5

[8] 张恒磊,张云翠,宋双,等.基于Curvelet域的叠前地震资料去噪方法[J].石油地球物理勘探,2008,43(5):508-513

Zhang H L,Zhang Y C,Song S,et al.Curvelet domain-based prestack seismic data denoise method[J].Oil Geophysical Prospecting,2008,43(5):508-513

[9] 仝中飞,王德利,刘冰.基于Curvelet变换阈值法的地震数据去噪方法[J].吉林大学学报(地球科学版),2008,38(S1):48-52

Tong Z F,Wang D L,Liu B.Seismic data denoise based on curvelet transform with the threshold method[J].Journal of Jilin University (Earth Science Edition),2008,38(S1):48-52

[10] 吴爱弟,赵秀玲.基于曲波变换的地震信号去噪方法[J].中国石油大学学报(自然科学版),2010,34(3):30-33

Wu A D,Zhao X L.Seismic signal denoising method based on curvelet transform[J].Journal of China University of Petroleum(Edition of Natural Science),2010,34(3):30-33

[11] 董烈乾,李振春,王德营,等.第二代Curvelet变换压制面波方法[J].石油地球物理勘探,2012,46(6):897-904

Dong L Q,Li Z C,Wang D Y,et al.Ground-roll suppression based on the second generation Curvelet transform[J].Oil Geophysical Prospecting,2011,46(6):897-904

[12] 姜宇东,杨勤勇,何柯,等.基于曲波变换的地面微地震资料去噪方法研究[J].石油物探,2012,51(6):620-624

Jiang Y D,Yang Q Y,He K,et al.Surface microseismic data denoising method based on curevelet transform[J].Geophysical Prospecting for Petroleum,2012,51(6):620-624

[13] Candès E,Demanet L,Donoho D,et al.Fast discrete curvelet transforms[J].Multiscale Modeling & Simulation,2006,5(3):861-899

[14] Ying L,Demanet L,Candès E.3D discrete curvelet transform[J].International Society for Optics and Photonics (Optics & Photonics 2005),2005,SPIE591413

[15] Dubois S,Péteri R,Ménard M.A 3D discrete curvelet based method for segmenting dynamic textures[C]∥Institute of Electrical and Electronics Engineers.Image Processing(ICIP),2009 16thIEEE International Conference on.USA:Institute of Electrical and Electronics Engineers,2009:1373-1376

[16] 邓承志.图像稀疏表示理论及其应用研究[D].武汉:华中科技大学,2008

Deng C Z.Research on image sparse representation theory and its applications[D].Wuhan:Huazhong University of Science and Technology,2008

[17] 张博.基于稀疏变换的地震数据去噪方法研究[D].杭州:浙江大学,2013

Zhang B.Research of seismic data denoising methods based on sparse transform[D].Hangzhou:Zhejiang University,2013

[18] Donoho D L.De-noising by soft-thresholding[J].IEEE Transactions on Information Theory,1995,41(3):613-627

[19] Gao H Y,Bruce A G.WaveShrink and semisoft shrinkage[J].Statistical Science Research Report,1995(39):1-25