2013, Vol. 56
地震资料噪声压制是地震数据处理中的一个重要流程.针对不同的地震干扰类型, 出现了不同的去噪方法, 但最大的挑战仍然是如何做到保幅去噪.能否做到保幅去噪主要取决于对混杂在一起的信号和干扰进行有效分离的程度.调和分析方法通过诸如旋转、伸缩、平移等变换处理, 将数据从原始的记录空间投影到不同的正交子空间, 为我们提供了另一种分析和处理地震数据的视角, 因此, 如何针对地震资料的特点设计有针对性的特征函数或特征记录成为一个重要的研究方向.
地震记录中反映地下结构的一次反射信息和伴随的干扰具有不同的记录模式, 主要体现在出现周期上的不同、出现位置上的不同(记录时间和记录检波器位置两个方面), 以及能量上的区别等.具体来讲, 连续的反射地层会对应具有一定规律的、在时间空间位置上相关性较强的曲线模式, 而随机噪声会无规律的出现, 多次波是一次反射波的重复回放, 面波则以长周期、强能量的形式出现在近地表.
傅里叶变换通过三角函数将信号从记录空间投影到频率空间, 根据出现的周期差异将信息分开.傅里叶变换在信号分析中具有划时代的意义, 其在各个领域都有广泛的应用, 后期演变出的各种变换方法大都是建立在傅里叶变换的基本理念之上.小波变换相对于傅里叶变换在地震数据处理上具有明显优势, 原因在于, 小波变换可以高效地检测出信号中的奇异点, 而傅里叶变换需要对不连续点进行全局表达, 但小波变换的局限性在于不能对地震记录中常见的曲线进行稀疏表达.Curvelet[1]变换综合了傅里叶变换和小波分析的优势, 可以采用少量系数表达奇异曲线特征, 因而Curvelet变换又被称为高维小波变换.
作为一种新的调和分析方法, Curvelet变换同时具有多尺度、多方向及时间空间定位的优点, 因而在分析诸如图像边界等[2]各向异性结构中具有显著的优势.自1999年第一代Curvelet变换[3]提出至今, 其在图像处理方面得到广泛应用, 并被推广应用于勘探地球物理中, 主要体现在地震数据噪声消除、地震数据插值、地震偏移、地震反褶积等方面[4-11].
此外, 经典的Curvelet变换是一种二进的分析方式, 主要体现在基函数的伸缩和旋转两个方面.具体来讲, 一旦确定分析采用的尺度数和低频分量的方向分辨率, 整个记录将根据二进模式进行自动分解.该方法优点在于简单、易于操作(只需两个参数), 缺点在于分析模式单一化、统一化.这种变换模式对于那些在频率和方向变化上满足二次幂或拟二次幂的信号或干扰, 可以有效分解, 但这种要求对于地震数据来讲过于苛刻.统一化在地震资料处理方面造成的问题主要体现在:1.不能将信号和干扰有效地投影到不同的正交空间; 2.最低频分量方向分辨率较低.为此, 本文设计了一套能够较有效反映不同地震记录特征的自适应基函数.为区别于经典的“二次幂”Curvelet变换, 我们称之为“非二次幂” Curvelet变换.
本文介绍了非二次幂Curvelet变换的基本原理, 并使用简单模型对比了经典Curvelet变换和非二次幂Curvelet变换在分解方式上的区别, 为验证本文方法在地震资料信噪分离上的有效性, 分别采用了模型算例和一个含有面波干扰的实际地震记录进行测试, 获得了较好的测试结果.
2 非二次幂Curvelet变换 2.1 Curvelet变换:各向异性结构的最优稀疏表达在多元分析中, 各向异性结构成为研究中的主要特征, 例如二维勘探地震记录以曲线或类曲线特征为主导, 地震波前也具有类曲波结构, 为了研究这类多元各向异性特征, 1999年Candès和Donoho最早提出了完备的、多方向的Curvelet构架[3].与各向同性的Wavelet体系不同, Curvelet是一种调和分析中各向异性的构造体系.
假设存在一组完备的各向异性Curvelet特征函数φj, l, kd, 其中j=(0, 1, …, N)代表尺度, N为最大尺度, l=(1, …, L)代表方位, kd代表d维空间中的位置(d≥2), 则d维空间信号fd可投影到基函数张成的空间中, 映射系数称之为Curvelet系数[2]:
|
(1) |
根据Curvelet变换的正交性和线性性质, 信号fd可表示为特征基函数φj, l, kd的线性加权和, 权系数可由(1)式唯一确定:
|
(2) |
Curvelet变换的过程和系数的排列方式如图 1所示.(图 2a, 图 2b)给出了5个典型的Curvelet基函数(使用不同的数字标记)在时间空间域和频率域的显示, 从中可以看出, Curvelet1, 3, 4, 5具有明显不同的频率和方向特性, 而Curvelet1, 2具有相同的频率和方向特征, 但具有不同的时间空间域位置特征.此外, Curvelet特征函数满足抛物尺度关系.以Curvelet 4为例, 其时间空间和频率域的有效支撑区域均满足length≈width2.
|
图 1 Curvelet变换过程及相应的Curvelet系数 (a)时空域的二维地震记录; (b)该地震记录对应的F-K频谱; 对该记录进行4个尺度的频率分解和第二尺度8个方向的角度分解, 可获得如图 1c所示的Curvelet系数(由于对称性, 只显示一半的系数), 横坐标代表角度, 从左向右-90°~90°变化, 纵坐标代表尺度, 从上往下尺度依次增大. Fig. 1 Curvelet transform and visualization of Curvelet coefficients (a)2D seismic record with x axis indicating trace number, and y axis notating recording time in second; (b) Shows the corresponding F-K spectrum; through applying four maximum scale analysis for the record, as well as 8 directions in angular decomposition at the second scale, we obtain the Curvelet coefficients displayed in(c) (due to the conjugate symmetric property, only half of the coefficients are visualized). The horizontal axis indicates angle change from -90°~90° and vertical axis represents scale(increase from top to bottom). |
|
图 2 Curvelet基函数在时间空间域和频率波数域的显示 (a)展示了时间空间域5个典型的Curvelet基函数(使用不同的数字标记具有不同频率、方向和位置特征的Curvelet基函数); (b)对应的频率波数域频谱, 表明Curvelets具有多尺度、多方向(由F-K域的楔形体限定, 且Curvelets在T-K域的方向垂直于F-K楔形体的主方向)、空间频率定位特性, 以及满足抛物尺度关系. Fig. 2 Curvelets in spatial and frequency domain (a) Shows five classical Curvelets in spatial domain(different digital numbers indicate Curvelets with distinct frequencies, orientations or locations), the corresponding functions in frequency domain are displayed in 2b, which indicates that Curvelets are multi-scale, multi-directional(confined in the wedge of F-K domain and the direction of curvelets in T-K space is perpendicular to the principal direction of F-K wedge), and localized both in space and frequency spaces, satisfying parabolic scaling relationship. |
Curvelet变换提供了一种对C2奇异曲线最优的稀疏表达方式.在N项逼近的误差分析中, Curvelet变换可以达到N-2的收敛速率, 理论证明, 这是C2奇异曲线的最优近似手段.与此相比, 小波变换的重构误差仅能以N-1的速率收敛.图 3对比了小波变换对地震波前的各向同性表达方式和Curvelet变换稀疏的各向异性近似.
|
图 3 对比小波变换和Curvelet变换对各向异性结构的非线性逼近方式 (a)小波变换对某一地震同相轴的各向同性表示; (b)Curvelet变换的各向异性表示. Fig. 3 Comparison of non-linear approximation styles of anisotropic structures by wavelets and curvelets (a)Shows the isotropic representation of a seismic hyperbolic event by wavelets; (b)Visualizes the anisotropic approximation by Curvelets. |
数值计算中, 特征基函数φj, l, kd可以通过对应的频率域函数Uj, l(w)来定义, 从图 2b可以看出,
Uj, l(w)的支撑区域是由Wj(w)(径向窗)和Vl(w)
(角度窗)联合限定的楔形区域, 其中径向窗函数定义为:Wj(w)=W(2-jw); 角度窗函数定义为:
首先, 不同于二次幂的、全局的尺度变换, 非二次幂Curvelet变换针对处理目标和频率分辨率设计尺度分解参数.假设在频率域有N个重要的研究区域, 以其中一个在频率域支撑区间为(w1, w2)∈ Ej(j=1, …, N)的目标为例, 可引入:
|
(3) |
其中w12, w11为目标在K1方向上的上限和下限;
w22, w21为目标在K2方向上的上限和下限.若假设该尺度区间在K1方向上可细分为Nj, 1个子带, 在K2方向上可细分为Nj, 2个子带, 则子尺度窗函数的频率分辨率为:
|
(4) |
在目标能量集中、需要整体分析的情况下, 取Nj, 1=Nj, 2=1, 则(4)式定义的尺度函数可简化为
|
(5) |
在经典Curvelet变换中, 给定某一尺度分析精度, 整个频率空间即可自动进行二进方式尺度剖分.以图 4为例, 图 4a为二次幂的Curvelet尺度剖分方式, 其尺度2、3和4的频谱范围如图 5a所示, 可以看出, 面波能量(红色椭圆内)被分散到三个尺度中, 各尺度内的精细角度分解会导致面波能量的进一步分散, 若要消除面波就需要对所有包含面波能量的Curvelet系数进行处理, 这将导致较低的计算效率, 而且误差较大.本文提出的非二次幂尺度分解, 可根据处理目标设计合适的参数和分解方式, 如图 4b所示, 其尺度2的频谱范围如图 5b所示, 可以看出这种分解方式使得面波能量更加集中、特征更加清晰, 便于对面波成分进行精细的后续处理.
|
图 4 经典二次幂Curvelet和本文提出的非二次幂Curvelet尺度分解方式对比 (a)经典Curvelet二次幂的尺度分解方式, 两个白框之间为一个尺度; 红色椭圆中的面波能量被分散到三个尺度中; (b)中的非二次幂尺度分解将面波成分集中到一个尺度中. Fig. 4 Comparison of the classical dyadic Curvelet and the proposed non-dyadic Curvelet based scale decomposition methods (a)Shows the dyadic scale decomposition pattern(a scale is defined by two white rectangular circles), energies of surface waves are dispersed over three scales; (b)Shows the non-dyadic scale decomposition style, in which surface waves are concentrated in one scale. |
|
图 5 经典二次幂Curvelet和本文提出的非二次幂Curvelet尺度分解结果对比 (a)经典Curvelet二次幂的尺度分解结果, 可以看出, 红色椭圆中的面波能量被分散到三个尺度中, 破坏了其原有特征; (b)非二次幂尺度分解结果, 面波能量较为集中. Fig. 5 Comparison of the classical dyadic Curvelet and the proposed non-dyadic Curvelet based scale decomposition results (a)Shows the dyadic scale decomposition resultss, in which surface waves are spread over three scales, failing to preserve the original features; (b) Shows the non-dyadic scale decomposition result, in which surface waves are concentrated in a single scale as expect. |
经典Curvelet变换的角度窗函数定义为
|
(6) |
由于对称性, 以
|
(7) |
根据定义的非二次幂尺度和角度窗函数, 频率域窗函数可表示为:Uj, l:=Wj×Vlj.
以一个包含4条具有不同方向的直线同相轴的简单模型, 对比常规二次幂角度分解方式和本文提出的非二次幂角度分解方式的区别.图 6a为常规Curvelet多尺度多方向分解结果, 当第二尺度被均匀细分为8个角度区间时(采用图 1c中的Curvelet系数排列方式), 二次幂分解模式决定了第三和第四尺度的角度分解个数为16, 可以明显看出, 这种分解模式下精细尺度均得到了恰当分解, 但是第二尺度角度分解不足, 不同方向的同相轴没有有效分开.若增加低尺度的角度分解个数使其充分分解, 将会导致高频部分角度分解过量, 给后续处理带来困难.非二次幂的角度分解方式, 使得各个尺度的角度分解均是独立的, 因此每个尺度内的角度分解参数可根据该频带信息的分析精度要求确定, 故而可以较有效地反映地震资料的主要方向特征, 使得各个频带(暂不考虑最低频带)均得到恰当分析, 不仅可以提高分析效果, 而且可以提高运算效率.基于本文提出的非二次幂Curvelet分解结果如图 6b所示, 可以明显看出, 这种自适应的分解方式可以将各个尺度下的方向信息有效提取出来, 清晰地反映出资料的主要方向特性.
|
图 6 经典Curvelet的分解结果(a)及本文提出的非二次幂变换分解结果(b) 从(a)可以看出, 第三、四尺度得到恰当分解, 但第二尺度分辨率不足, 未实现彻底分解; (b)中根据各个尺度的角度分辨率设计自适应的参数, 使得各个成分均得到适当分解. Fig. 6 (a)Shows the directional analysis results via the classical Curvelet transform, from which we can see clearly that the third and fourth scale are decomposed properly while the second scale is far from proper decomposition with lower angular resolution; (b) displays the result using the proposed non-dyadic decomposition, in which all angular analysis parameters are designed according to the corresponding resolutions, rendering all components properly unraveled. |
当j=0时, 尺度窗函数为W0(w), 传统Curvelet变换对最内尺度进行无方向性的各向同性分析.为了提高低频成分的分析精度, 非二次幂Curvelet变换对最内尺度进行方向性分析.定义θ10, θ20为角度域的上、下限, 则方向分辨率为
|
(8) |
则最内尺度的频率域窗函数可表示为U0, l:=W0× Vl0, l=1, 2, 3, …, L0.
以图 7为例对比常规Curvelet变换和本文提出的方法对低频成分处理方式的不同.图 7(a-c)展示了对最内尺度进行常规Curvelet各向同性分析的结果, 其中图 7a为最内尺度对应的Curvelet系数, 图 7b为仅保留最内尺度Curvelet系数重构的记录, 图 7c为该记录所对应的F-K频谱.可以看出, 这种各向同性的处理方式混合了各个方向的信息成分, 不利于进一步的精细分析和处理.非二次幂Curvelet变换对最内尺度引入各向异性分解, 提高了低频成分的分析精度.图 7d即是对最内尺度多方向分解的结果, 可以明显看出这种各向异性的处理方法是非常有效的, 低频部分不同方向的信息被提取出来, 为后续的处理、解释等提供更精细的低频信息.
|
图 7 经典Curvelet变换和本文提出的非二次幂Curvelet变换对低频成分的处理方式对比 (a)最内尺度的Curvelet系数，(b、c)为其对应的时间空间域和频率波数域记录，可以看出，经典的Curvelet变换对最内尺度不进行多方向分析; d)本文提出的非二次幂Curvelet变换对最内尺度进行各向异性分析得到的变换系数，结果表明非二次幂Curvelet变换可以有效实现低频成分的精细处理. Fig. 7 Comparison of the analysis of low-frequency components via the classical Curvelet transform and the proposed non-dyadic Curvelet transform (a)Displays the Curvelet coefficients at the lowermost scale; (b) and(c) show the corresponding data in the time-spatial and frequency domain, from which we can see obviously that the classical Curvelets transform treats the lowermost scale isotropically, thus the corresponding records contain information from multi directions; (d) shows us the non-dyadic Curvelet coefficients through the proposed anisotropic analysis, which demonstrates that the introduced non-dyadic Curvelet transform realizes higher-resolution analysis for the low-frequency component. |
为了验证非二次幂Curvelet变换方法在地震数据处理中的效果, 采用二维合成和实际地震资料对其进行信噪分离测试.图 8a所示的模型数据包含两个反射同相轴和线性干扰, 首先通过分析数据本身的主要频率和方向特征, 确定多尺度多方向分解参数; 然后对其进行非二次幂Curvelet变换, 将线性干扰与有效信号在非二次幂Curvelet域中分开, 并采用阈值法去除线性干扰对应的系数; 最后将变换系数重构回原始时空域内.信噪分离结果如图 8b所示, 图 8c为二者的差值, 图 9(a-c)分别为图 8的三个记录所对应的F-K频谱.本算例表明了非二次幂Curvelet变换可以将不同方向的信息进行有效的分离.
|
图 8 采用非二次幂Curvelet变换对合成记录进行信噪分离测试 (a)合成数据：记录中有200个地震道，记录长度为3s，记录中有两个有效同相轴和两条线性干扰; (b)采用本文提出的非二次幂Curvelet变换消除线性干扰后的记录; (c)为二者的差值剖面. Fig. 8 A test of the proposed non-dyadic Curvelet transform in signal-noise-separation for a synthetic record (a)Synthetic data with 200 traces and 3 s recording length for each receiver, which contains two desirable curves and two noisy lines; (b)Shows that the record after denoising using the proposed non-dyadic Curvelet transform; (c)Displays their differences. |
|
图 9 频率波数域显示 (a-c)分别为去噪前、去噪后和去除的噪声对应的频谱. Fig. 9 Frequency-wavenumber representation of seismic records in Fig. 8, before denoising(a), after denoising(b)and noise removed(c)respectively. |
图 10a为某地区二维数据体中的单炮地震记录, 可以看出记录中面波能量较强, 分布范围广, 与有效信息混叠严重, 破坏了有效信息的连续性和完整性.基于该记录面波和有效信息在频率、方向上的差异, 我们对其进行非二次幂Curvelet分析, 将面波成分和有效反射信息成分投影到不同的非二次幂Curvelet空间, 通过消除面波对应的系数来保留有效信号成分.压制面波干扰后的地震记录如图 10b所示, 图 10c为去噪前后的差值记录.去噪结果表明, 由于非二次幂Curvelet变换继承了Curvelet变换多尺度多方向以及最优稀疏化表示地震信息的优势, 并且充分考虑到了待处理资料的主要频率和方向特性, 因而可以较好地将面波和有效信息分开, 有效信号的保真度较高.
|
图 10 采用非二次幂Curvelet变换对实际二维地震记录进行面波分离测试 (a)去噪前的地震记录; (b)去噪后的地震记录; (c)差值记录. Fig. 10 A test of the proposed non-dyadic Curvelet transform in 2D seismic data surface wave separation (a)Shows the original seismic data; (b)Displays the seismic record after denoising; (c)The difference record. |
本文根据地震数据的方向和频率特点, 设计了一套有针对性的、高效的特征函数, 实现了尺度域、角度域的非二次幂自适应分解, 以及对最低频成分的各向异性分解, 从而提升了Curvelet变换对地震数据有效频带、方向和低频信息的精细分析和处理能力.非二次幂Curvelet变换在地震资料面波干扰压制中的较好应用, 进一步验证了本文方法的有效性.
致谢作者感谢两位匿名评委提出的宝贵意见.
| [1] | Candès E J, Donoho D L. Curvelets-A surprisingly effective nonadaptive representation for objects with edges.//Cohen A, Rabut C, Schumaker L eds. Curves and Surfaces. Nashville:Vanderbilt University Press, 2000:105-120. |
| [2] | Candès E J, Demanet L, Donoho D, et al. Fast discrete curvelet transforms. Multiscale Modeling and Simulation , 2006, 5(3): 861-899. DOI:10.1137/05064182X |
| [3] | Candès E J, Donoho D L. Curvelets:A Surprisingly Effective Nonadaptive Representation for Objects with Edges. Nashville: Vanderbilt University Press, 1999 . |
| [4] | Herrmann F J, Verschuur E. Curvelet-domain multiple elimination with sparseness constraints. Society of Exploration Geophysicists, Expanded Abstracts , 2004: 1333-1336. |
| [5] | Neelamani R, Baumstein A I, Gillard D G, et al. Coherent and random noise attenuation using the curvelet transform. The Leading Edge , 2008, 27(2): 240-248. DOI:10.1190/1.2840373 |
| [6] | Chauris H, Nguyen T. Seismic demigration/migration in the curvelet domain. Geophysics , 2008, 73(2): S35-S46. DOI:10.1190/1.2831933 |
| [7] | Wang Y B, Dong S Q, Xue Y W. Surface waves suppression using interferometric prediction and curvelet domain hybrid L1/L2 norm subtraction. SEG Expanded Abstracts 28 , 2009: 3292-3296. |
| [8] | 包乾宗, 高静怀, 陈文超. Curvelet域垂直地震剖面波场分离. 西安交通大学学报 , 2007, 41(6): 650–654. Bao Q Z, Gao J H, Chen W C. Wave field separation of vertical seismic profiling in Curvelet domain. Journal of Xi'an Jiaotong University (in Chinese) , 2007, 41(6): 650-654. |
| [9] | Verschuur D J, Berkhout A J, Wapenaar C P A. Adaptive surface-related multiple elimination. Geophysics , 1992, 57(9): 1166-1177. DOI:10.1190/1.1443330 |
| [10] | Wang D L, Tong Z F, Tang C, et al. An iterative curvelet thresholding algorithm for seismic random noise attenuation. Applied Geophysics , 2010, 7(4): 315-324. DOI:10.1007/s11770-010-0259-8 |
| [11] | Yuan Y H, Wang Y B, Liu Y K, et al. Curvelet domain adaptive least-squares subtraction of internal multiples. Journal of Seismic Exploration , 2011, 20(3): 273-288. |