Curvelet域面波衰减方法研究
张恒磊1, 2，刘天佑1, 2，李红巧3
(1. 中国地质大学(武汉) 地球物理与空间信息学院，湖北 武汉，430074；
2. 中国地质大学(武汉) 构造与油气资源教育部重点实验室，湖北 武汉，430074；
3. 江西科技师范学院 光电子与通信重点实验室，江西 南昌，330013)
摘要：在面波压制方法中，常规面波消除方法都是单一的利用面波某种特性，如Fourier变换域方法利用波组的频率差异等，难以有效地压制面波，并且根据单一的波组特征差异进行面波消除很容易损伤有效波信息。Curvelet变换可以对时空信号进行最稀疏表达，能够获得最优的非线性逼近。分析地震信号面波与有效反射波在Curvelet域(j，θ，k)三维空间的特征差异，将空间域波组方向与Curvelet域角度变量联系起来，指出可以在Curvelet域利用波组的频率、角度和空间位置差异实现波场分离，并设计非线性阈值函数对面波系数进行衰减，避免传统软、硬阈值函数的不足。研究结果表明：基于Curvelet域的面波压制方法受随机噪声的影响较小，它可以有效地压制面波干扰，同时对有效波信息的保真度高。
关键词：
中图分类号：P631.4 文献标志码：A 文章编号：1672-7207(2011)08-2372-07
Attenuation of surface wave in Curvelet domain
ZHANG Heng-lei1, 2, LIU Tian-you1, 2, LI Hong-qiao3
(1. Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China;
2. Key Laboratory of Tectonics and Petroleum Resources of Ministry of Education,
China University of Geosciences, Wuhan 430074, China;
3. Key Laboratory of Photoelectron & Communication, Jiangxi Science & Technology Normal University,
Nanchang 330013, China)
Abstract: The traditional methods can suppressing the surface wave to some extent, but they only use single characteristic. For instance, Fourier transform methods use the frequency characteristic between surface wave and effective wave. When these methods are employed to attenuate surface wave, the signal components may be harmed. Curvelet transform is characterized by optimum sparseness constraint condition that it can deal with line-like phenomena in high dimension. The signal characteristic in (j,θ,k) of Curvelet domain between surface wave and desirable signal was analyzed, and the Curvelet angle was linked with the direction of wave group, so surface wave will be separated from effective wave based on the differences in frequency, angle and position characteristics. And besides, a non-linear thresholding was designed to suppress the surface wave coefficients which can avoid the discontinuity caused by the hard or soft thresholding. The results show its feasibility and effectiveness in attenuating surface wave, and it can provide superior surface wave attenuation with minimal impact on the desirable signal components.
Key words: surface wave; Curvelet domain; non-linear thresholding; attenuation; separation of wave field
油气勘探工作的日益复杂化需要不断提高地震勘探目标，要求野外采集的地震资料有较高的信噪比和分辨率，才能满足精细构造解释和储层预测工作要求。在反射地震勘探中面波是一种具有明显的高振幅和频散特性的规则干扰波，在油气勘探地震数据处理中，面波的压制工作是关键的一步：面波压制的效果直接影响地震数据后续处理的效果，最终会影响地质构造解释的结果。传统消除面波的方法主要是利用面波与反射波的频带差异及视速度差异，如高通滤波、切除法、F-K滤波等方法[1-3]，这些方法曾起到一定的作用，但都存在相应的局限性，如：高通滤波法会损失中深层的低频有效波信息；切除法在切除面波成分的同时会将包括在面波中的有效波切除；F-K滤波只适应地层倾角较缓的地区，对复杂地质条件下的面波去除效果不佳，且该方法由于采用二维傅里叶变换，处理结果存在炕席现象[4]，容易产生新的相干干扰。20世纪80年代兴起的小波变换由于其时频局部分析能力而被广泛应用于地震资料处理领域[5-8]。但是二维小波变换只能提供水平、垂直以及对角3个方向的分量信息，角度分辨率较低[9],不能对复杂多变的有效波和面波特征进行有效的表示和检测，导致面波压制的效果不佳。为了克服小波在高维线奇异的局限性，Candès 等[10-11]在小波理论的基础上，提出了一种特别适合于表示具有线状变化特征信号的分析方法——Ridgelet变换。包乾宗等[12]研究了在Ridgelet域压制面波的方法，得到了较好的处理效果。该方法将地震记录拓展到(a，τ，p)三维空间，同时利用地震记录的视速度、时间和尺度域特性差别，实现有效信号与面波的分离。但是，Ridgelet变换也有它的局限性，Candès 等[10, 13-14]指出：Ridgelet变换只能针对直线状的信号特征进行有效处理。自1999年以来Candès等[15-17]提出了Curvelet变换，其思路是基于曲线分割的思想，即在无限小的尺度上将每一小段曲线近似为直线，对每一段直线应用Ridgelet变换。以微线段而非点作为基本的信号表达单位，因此，能够像小波很好地处理点奇异性一样，Curvelet变换能很好地处理高维信号的线奇异特征。近年来，该方法被迅速应用于地球物理领域[18-22]，但主要集中在噪声压制和成像等方面，张恒磊等[23]论述了Curvelet变换的性质及在噪声压制方面的应用，目前还没有见到国内关于Curvelet域压制面波的方法报道。在此，本文作者分析Curvelet变换的数学原理，指出其尺度与方向参数可以表示复杂的地震波组特征，可以将时空域波组方向与(j，θ，k)三维空间的尺度与角度联系起来。从物理角度分析在Curvelet域面波和有效信号的分布特征，提出利用Curvelet变换压制面波的方法，简单有效，对大量地震资料处理具有一定的实践意义。
1 方法原理
面波干扰会降低地震资料的信噪比，影响资料处理的质量，特别对叠前资料处理影响更为严重。传统的根据单一的波组特征差异进行面波消除的方法很容易破坏有效波信息，不能实现信号的保真处理。Curvelet变换可以提供多尺度多方向的细节信息，面波与有效反射波在Curvelet域(j，θ，k)三维空间有明确的物理意义，因此，从频率、角度和空间位置差异实现有效反射波和面波的分离将有很好的应用效果。
1.1 Curvelet分析与旋转角度剖分
在二维空间情况下，设空间域变量为x，频率域变量为w，频率域的极坐标为r和θ。假设φj(x)为Curvelet“母”函数，则j尺度下的各Curvelet都可以由该“母”函数φj(x)经过旋转平移得到。取角度剖分量A，定义旋转角度
，
m=1，2，…，A?2[j/2]； 0＜θj,m≤2π (1)
其中：
表示对j/2取整。
平移参数k=(k1，k2)
Z2，则尺度j，方向θj,m，位置
处的Curvelet变换[15]为：
(2)
式中：Rθ=
，为θ的旋转。
在构造Curvelet的过程中，每个尺度下的角度分量由尺度因子和角度剖分量共同决定(图1所示为反应的是角度剖分量A取8时的Curvelet频率域数字化块)，突出了Curvelet的多方向特征。由式(1)和图1可以看出：Curvelet的方向分量是随着分解尺度的增大而增大的，这种性质表示它能够有效地描述二维图像中的方向细节。
对于
，其Curvelet变换可表示成与Curvelet函数
的内积[15-16]：
(3)
图1 Curvelet频率域数字化块
Fig.1 Curvelet tiling of frequency
通过以上分析，给定一个Curvelet，经过伸缩、平移和旋转可以生成L2(R2)(平方可积函数空间)的紧标架，这就意味着它具有重构公式[15-16]：
(4)
式中：φj,m,k表示由指标(j，m，k)确定的Curvelet族；M表示指标集。
Curvelet变换是小波变换的二维各项异性扩展，它和小波变换相似，同样可以伸缩平移。只是由于Curvelet是在二维情况下分解，因此，平移由2个参数决定；此外Curvelet变换和小波变换的主要不同在于Curvelet变换具有旋转性，可以对信号在多方向上进行有效分析。
1.2 波场Curvelet域物理意义及面波压制方法
由式(1)和(2)可知：每一个Curvelet都有特定的频率、方向和位置，而地震资料中不同成分的信息都可以从这3个方面来区分。实际地震资料中有效反射波是相对较长、较平缓的波组特征，在Curvelet域其能量主要集中在反映小角度的分量中；面波是相对较长、较陡的具有频散特征的低频波组，在Curvelet域其能量主要集中在反映相对大角度的分量中。Curvelet变换可以将时空域信号f(x)变换到(j，θ，k)域，因此，在Curvelet域可以很容易的将不同信号进行分离。
图2所示为不同波组信号经过Curvelet分离的结果。可以看出：Curvelet域(j，θ，k)三维空间可以将时空域不同频率、方向以及空间位置的波组很好地分离开来：尺度信息体现波组频率差异，θ参数反映波组的方向差异，即主要体现视速度差异。这就是本文基于Curvelet域压制面波方法的物理基础。另一方面，在Curvelet域，有效反射波在相对小尺度段呈现的是小区域、强幅值的Curvelet系数；面波在相对大尺度段呈现的是小区域、强幅值的Curvelet系数；而随机噪声由于分布在全尺度、所有方向，呈现出大区域、低幅值的Curvelet系数，因此，可以利用k参数体现出的空间位置差异及幅值差异对波场进行进一步分离。
设地震剖面记录f(x，t)，其中x表示CDP号，t表示时间。设g(x，t)为面波，r(x，t)为有效反射波，n(x，t)为随机噪声，则实际记录可以表示为：
根据Curvelet线性变换的性质，在Curvelet域，有：
式中：Cf(j，θ，k)、Cr(j，θ，k)、Cg(j，θ，k)和Cn(j，θ，k)分别为f(x，t)，r(x，t)，g(x，t)和n(x，t)的Curvelet变换系数。在原始信号中，有效反射波与面波及随机干扰叠加在一起无法分离，通过Curvelet 变换把地震信号拓展到Curvelet域(j，θ，k)三维空间，可以实现不同波组的分离。为了防止有效信号损失，本文设计了非线性阈值衰减对面波系数进行处理，即认为在特定方向分量上(如图2所示)，大于阈值的系数为面波引起，处理法则如下：
(5)
式中：Cfj,θ表示有面波存在的j尺度θ方向的Curvelet分解系数；
表示面波衰减后的估计系数。
图2 Curvelet波组特征分解
Fig.2 Curvelet representations
该方法主要包括以下4个步骤：
(1) 对信号作频谱分析和速度分析，合理选择变换尺度和角度剖分量A；
(2) 对信号作Curvelet分解，得到(j，θ，k)三维空间的向量组；
(3) 分析(j，θ，k)空间各向量组特征，选择有面波存在的j，θ参数，并通过能量确定阈值t，利用式(5)法则对Cfj,θ进行衰减；
(4) 通过Curvelet反变换重构信号得到压制面波后的地震剖面。
在处理过程中，变换尺度j、角度剖分量A以及对Cfj,θ的处理法则是关系到分离效果的3个重要因素。变换尺度一般选择3阶或4阶，角度剖分量一般取8时即可，具体可通过步骤1中的分析，以所得到θj,m参数能够反映波组方向为基本原则(如图2所示)。另外对系数Cfj,θ的处理法则，本文采用非线性阈值函数克服了传统硬阈值存在间断点的现象，同时也避免了软阈值减去一个阈值而引起重构失真的问题[23]。
2 模型实验
为对本文方法进行验证，进行了模型数值实验分析。图3所示为在合成的记录剖面(2条弯曲同相轴中心频率为30 Hz，幅值为1，模拟反射波；1条中心频率为10 Hz，幅值为2，模拟面波)上加入随机噪声所得的信号剖面(信噪比为0.5)及其处理结果。图3(b)所示为采用二维小波分析，利用传统阈值收缩方法分解到三阶得到的处理结果。经分析可知：小波变换受其角度分辨率的限制，不能很好地刻画二维信号中的曲线变化特征，在低信噪比情况下利用阈值收缩压制低频低波数面波的同时，会衰减一部分有效信息。从图3(b)和4(c)可以看出：剖面中的面波被较好地消除，但同时反射波能量变弱，尤其损伤了大倾角、弯曲同相轴部分的有效波信息。
图3 合成数据处理结果
Fig.3 Synthetic model data
图4 剖面信号第30道波形对应的频谱
Fig.4 Spectrum analysis of 30th trace
Curvelet变换可以提供(j，θ，k)三维空间多种信息。分析图2可知：不同特征的波组在Curvelet域可以通过j和θ参数很容易区分。基于这种Curvelet变换的物理基础，采用在Curvelet域对面波及有效波进行分离。将分解尺度参数设为3，对面波系数按照式(5)进行衰减，同时对随机噪声利用阈值收缩进行压制，得到图3(c)所示的处理效果。可以看出：面波和随机噪声干扰都被较好地分离出去，且有效波信息基本不受破坏，信号保真度高(图4(d))。图3(d)所示数据为图3(a)与图3(c)所示数据之差，从图3可以看出：该方法可以很好地分离出面波干扰，同时对随机噪声干扰也有较好的压制效果。
3 实例应用
对某工区单炮记录(图5)进行分析处理。采用小波分析及本文方法分别进行处理，以作对比研究。图6所示为小波分析处理的结果。可以看出：在小波域虽然可以衰减面波能量，但由于小波变换只能刻画一维信号的点奇异特征，不能对二维信号的线奇异特征进行有效描述，利用小波阈值收缩后重构信号会使波组发生畸变。此外二维小波变换的角度分辨率低，对于具有弯曲同相轴等复杂特征的地震资料，很难实现波场有效分离，且在弯曲同相轴部分很容易丢失信号，如剖面A和B部分。
图5 原始地震剖面
Fig.5 Original seismic data
在本文Curvelet研究应用中，通过分析，对该段地震剖面取剖分量A=8时，认为面波干扰主要集中在j=3，θ为0，π/8，π以及9π/8的方向分量中，利用阈值衰减加以去除，得到面波压制后的结果如图7所示。对比分析可以看出利用Curvelet域方法压制面波效果要明显优于小波域传统方法的压制面波效果，面波能量被很好地衰减，且在面波区域的有效波信息也较好的突出出来；同时，随机噪声干扰也得到相应的压制(图8)。处理结果重构精度较高，剖面视觉效果好，同相轴清晰、连续性好。
图6 小波阈值方法处理结果
Fig.6 Result by wavelet thresholding
图7 本文方法处理结果
Fig.7 Result by this work
图8 本文方法去除的面波
Fig.8 Difference between Fig.5 and Fig.7
4 结论
(1) Curvelet变换的多尺度多方向特征可以实现对地震信号的最稀疏表达，即在Curvelet域可以实现不同波组间较好地分离。因此，在消除干扰波的同时，可以保持有效波信息，同时受随机噪声的影响较小，对高精度地震勘探资料的保真处理具有重要意义。
(2) 文中系数处理策略采用的是阈值收缩法，所以算法的主要时间集中在Curvelet正变换与逆变换。对于文中实际数据处理，在硬件配置 CPU为AMD3000+，内存512 M，Matlab计算平台下正反变换所占用的CPU时间分别为0.985 s和1.594 s。
(3) 在Curvelet域实现对地震资料面波分离，方法简单有效。若能研究应用更为先进的系数处理法则，在更细致的层次上识别处理面波及其他干扰波，将会取得更好的处理效果。此外，Curvelet可以处理3D数据体，对三维地震勘探中的大数据量地震资料处理具有重要的实践意义。
致谢
加州理工学院应用和计算数学中心Candès E J教授以及CurveLab实验室的所有成员提供了Curvelet源代码，在此表示衷心感谢。
参考文献：
[1] Galbraith J N, Wiggins R A. Characteristics of optimum multi-channel stacking filters[J]. Geophysics, 1968, 33(1): 36-48.
[2] 刘法启, 张关泉. 小波变换与F-K算法在滤波中的应用[J]. 石油地球物理勘探, 1996, 31(6): 782-791.
Liu Fa-qi, Zhang Guan-quan. Application of wavelet transform and F-K algorithm in filtering[J]. Oil Geophysical Prospecting, 1996, 31(6): 782-791.
[3] Herrmann R B, Russell D R. Ground roll: Rejection using adaptive phase-matched filters[J]. Geophysics, 1990, 55(6): 776-781.
[4] 夏洪瑞, 周开明. 地震资料处理中相干干扰消除方法分析[J]. 石油物探, 2003, 42(4): 526-528.
Xia Hong-run, Zhou Kai-ming. Elimination of coherent interference in seismic data processing[J]. Geophysical Prospecting for Petroleum, 2003, 42(4): 526-528.
[5] 高静怀, 汪文秉, 朱光明, 等. 地震资料处理中小波函数的选取研究[J]. 地球物理学报, 1996, 39(3): 389-397.
Gao Jing-huai, Wang Wen-bing, Zhu Guang-ming, et al. On the choice of wavelet functions for seismic data processing[J]. Chinese J Geophys, 1996, 39(3): 389-397.
[6] 高静怀, 毛剑, 满蔚仕, 等. 叠前地震资料噪声衰减的小波域方法研究[J]. 地球物理学报, 2006, 49(4): 1155-1163.
Gao Jing-huai, Mao Jian, Man Wei-shi, et al. On the denoising method of prestack seismic data in wavelet domain[J]. Chinese J Geophys, 2006, 49(4):1155-1163.
[7] Deighan A J, Watts D R. Ground-roll suppressing using the wavelet transform[J]. Geophysics, 1997, 62(6): 1896-1903.
[8] 张三宗, 徐义贤. 地震记录小波域高阶相关叠加技术[J]. 地球物理学报, 2006, 49(2): 554-560.
Zhang San-zong, Xu Yi-xian. Higher-order correlative stacking for seismic data in the wavelet domain[J]. Chinese J Geophys, 2006, 49(2): 554-560.
[9] Mallat S. A theory for multiresolution signal decomposition: The wavelet representation[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1989, 11(7): 674-693.
[10] Candès E J, Donoho D L. Ridgelets: A key to high-dimensional intermittency[J]. Phil Trans Roy Soc London A, 1999, 357(1760): 2495-2509.
[11] Donoho D L. Orthonormal Ridgelets and linear singularities[J]. SIAM J Math Anal, 1998, 31(5): 1030-1061.
[12] 包乾宗, 高静怀, 陈文超. 面波压制的Ridgelet域方法[J]. 地球物理学报, 2007, 50(4): 1210-1215.
Bao Qian-zong, Gao Jing-huai, Chen Wen-chao. Ridgelet domain method of ground-roll suppression[J]. Chinese J Geophys, 2007, 50(4): 1210-1215.
[13] Arivazhagan S, Ganesan L, Subash Kumar T G. Texture classification using Ridgelet transform[J]. Pattern Recognition Letters, 2006, 27(16): 1875-1883.
[14] Zhang Heng-lei, Liu Tian-you, Zhang Yun-cui. Denoising of seismic data via multi-scale Ridgelet transform[J]. Earthquake Science, 2009, 22(5): 493-498.
[15] Candès E J, Donoho D L. Curvelets: A surprisingly effective nonadattive representation for objects with edges[C]//Cohen A, Rabut C, Schumaker L L. Curves and Surface Fitting. Nashville: Vanderbilt University Press, 2000: 105-120.
[16] Candès E J, Donoho D L. Continuous curvelet transform[J]. Appl Comput Harmon Anal, 2005, 19(2): 198-222.
[17] Starck J L, Fadili M J. Curvelets and ridgelets [EB/OL]. [2007-10-24]. http://jstarck.free.fr/curvencyclop09.pdf.
[18] Starck J L, Candès E J, Donoho D L. The curvelet transform for image denoising [J]. IEEE Transactions on Image Processing, 2002, 11(6): 67-684.
[19] Neelamani R, Baumstein A I, Gillard D G, et al. Coherent and random noise attenuation using the Curvelet transform[J]. The Leading Edge, 2008, 27(2): 240-248.
[20] Herrmann F J, Hennenfent G, Non-parametric seismic data recovery with Curvelet frames[J]. Geophysics J Int, 2008, 173(1): 233-248.
[21] Herrmann F J, Wang D L, Hennenfent G, et al. Curvelet-based seismic data processing: A multiscale and nonlinear approach[J]. Geophysics, 2008, 73(1): A1-A5.
[22] Douma H, de Hoop M. Leading-order seismic imaging using Curvelets[J]. Geophysics, 2007, 72(6): S231-S248.
[23] 张恒磊, 张云翠, 宋双, 等. 基于Curvelet域的叠前地震资料去噪方法[J]. 石油地球物理勘探, 2008, 43(5): 508-513.
Zhang Heng-lei, Zhang Yun-cui, Song Shuang, et al. Curvelet domain-based prestack seismic data denoise method[J]. Oil Geophysical Prospecting, 2008, 43(5): 508-513.
(编辑 杨幼平)
收稿日期：2010-09-03；修回日期：2010-11-05
基金项目：教育部高等学校博士学科点专项科研基金资助项目(20050491504)；国家自然科学基金资助项目(40374041)
通信作者：张恒磊(1984-)，男，安徽明光人，博士，从事地球物理弱信号处理及综合地球物理方法研究；电话：15927458203；E-mail：sky0610@163.com
