2013, Vol. 56
Curvelet变换于1999年由Candès[1]等人首次提出，此后被广泛应用于图像处理领域[2-5].在2004年的SEG年会上，Herrmann等做了关于Curvelet域多次波衰减、基于Curvelet变换和小波变换的AVO反演等报告[6-7]，自此Curvelet变换开始被广泛应用于地震数据处理中[8-18]，其中信噪分离是Curvelet变换一个重要的应用方向[13-18].一般而言，信号和噪声可以在诸如波长、能量、位置、方向因素的一个或几个方面加以区分，这些因素在Curvelet分解中都是非常重要的参数，经Curvelet变换后，各个因素均会在Curvelet域得到明显体现.Curvelet变换在信噪分离上的独特优势，为预测匹配滤波方法提供了强有力的保证.
Curvelet变换根据最大分解尺度及低频成分的方向分辨率等参数进行自动二进分解，具体来讲，每增加两个分解尺度，角度方向网格加密一倍.该分解模式的优点在于简单、易于理解和实现.对于能量在各个频率区间和方向上分布比较均匀的图像或数据而言，Curvelet变换是一种较为理想的分析和处理工具.然而，在分析那些对于特定频率范围或方向上有明显偏向的数据时，经典二进Curvelet变换遇到了瓶颈，地震数据就是个比较典型的例子.震源性质、地应力方向、地层对不同频率和方向地震波的衰减作用、检波器在接收频段和方向上的设计等因素，导致我们难以得到在各个频率和方向上能量分布均衡的地震数据.经典Curvelet变换可以高效率地对具有一般特征的数据或图像进行自动二进分解，但难以做到针对地震数据进行自适应分析.针对这个问题，本文采用了一种改进的Curvelet变换，即自适应的非二次幂Curvelet变换[14].该变换的主要思想在于根据数据本身的特征和先验信息，设计灵活的、有针对性的滤波窗函数，并根据各区间方向分辨率进行自动剖分.该方法继承了经典Curvelet变换在高维数据分析上的优势，并弥补了其在处理地震资料上的局限性.就效果而言，非二次幂变换方法灵活性高、针对性强、分解结果分辨率高；就效率而言，相对于全局的二次幂分解，局部的非二次幂算法可以节省对信息量稀疏区域的繁琐处理时间，当然，效率提高的程度取决于数据本身的特征.
基于非二次幂Curvelet变换，本文提出了一种改进的最小二乘匹配滤波算法（Non-Dyadic Curvelet based Least-Squares Method，NDCLSM）.该方法首先将输入数据和预测的噪声模型变换到非二次幂Curvelet域，然后在噪声能量集中的非二次幂Curvelet子记录上进行最小二乘匹配滤波.该方法具有非二次幂Curvelet变换多尺度、多方向、根据数据特征进行自适应参数选择等特点，适用性强，可广泛应用于诸如面波[18]、层间多次波和自由表面多次波[19-25]等的自适应匹配相减.本文对非二次幂Curvelet变换和NDCLSM算法的原理进行了介绍，并使用三维实际地震资料验证了本文方法的有效性.
2 非二次幂Curvelet变换Curvelet变换是一种多尺度、多角度的时频分析方法，通常在频率域内实现，Curvelet正变换可表示为：
|
(1) |
其中Uj, l（w）为尺度j角度l定义下的Curvelet基函数对应的频率域窗函数，通常Uj, l（w）的支撑区域是由Wj（w）（径向窗）和Vl（w）（角度窗）的支撑区间联合限定的楔形区域.xj, lk=Rθl-1（k12-j，k22-j/2）为Curvelet域坐标为k=（k1，k2）的Curvelet系数在时间空间域内对应的位置，其中
Curvelet变换基于二进的尺度角度分解，尺度窗函数定义为：Wj（w）=W（2-jw）；角度窗函数定义为：
图 1显示了Curvelet变换的多分辨率、多角度和空间定位等特性，与小波变换等变换方法相比，Curvelet变换更适于高维信号处理，特别是当信号具有交叉和重叠等现象的时候.
|
图 1 几个Curvelet函数在时间空间域（左）和频率波数域（右）的显示，表明Curvelet具有多尺度、多方向、空间频率定位特性，以及满足抛物尺度关系. Fig. 1 Several Curvelets in time-spatial domain (left) and frequency-wavenumber domain (right) indicative of Curvelet's multi-scale, multi-directional and localized in T-X and F-K domains, satisfying parabolic scaling relationship. |
非二次幂Curvelet变换基于目标特征进行自适应的尺度角度分解.假设目标区域在频率域可定义为：（w1，w2）∈Ej，则可引入
|
(2) |
其中w12，w11为目标在k1方向上的上限和下限；w22，w21为目标在k2方向上的上限和下限.尺度窗函数可定义为：
|
(3) |
引入
|
(4) |
对应于目标区域Ej的方向范围，
|
(5) |
则频率域窗函数为：Uj, l:=Wj×Vlj.对非二次幂Curvelet变换的细节描述可参考文献[14].
经典Curvelet变换的分解方式未能充分考虑数据本身特点，从而影响了变换的效果和效率，非二次幂Curvelet变换继承了Curvelet变换的优点，并克服了其局限性和非自适应性.这里我们分别使用两种方法对一个“米”字形模型进行分解，如图 2a所示，该模型包含四个方向的分量：水平、垂直和两个对角.图 2b为二次幂Curvelet变换结果，从分解结果可以看出，第三、第四尺度均得到合适的分解，但第二尺度方向分解不足，未将交叉分量区分开.而增加第二尺度角度分解总数则会导致高频分量分解过量而模糊化.图 3为非二次幂Curvelet变换的分解结果，各个尺度、各个方向均得到恰当剖分.
|
图 2 Curvelet变换分解方式 （a）时间空间域的合成模型；（b）对该“米”字形模型进行Curvelet多角度多尺度分解，以第二尺度角度分解个数为8得到的Curvelet系数. Fig. 2 Decomposition by Curvelet transform (a)Synthetic data in time-spatial domain; (b)Curvelet coefficients by dyadic multi-scale and multi-angular Curvelet decomposition (number of angles in the second scale equals 8). |
|
图 3 （a）为图 2a中的模型数据；（b）为基于数据的非二次幂Curvelet多尺度多角度分解结果.可以看出，该分解方式不同于传统Curvelet变换二次幂的剖分模式，且第一尺度进行了各向异性分解，该剖分方式较好地体现了数据本身的特征，可以将模型中“米”字形的各个分量区分开来. Fig. 3 (a)Synthetic data shown in Fig.2a; (b)Non-dyadic decomposition result.Clearly, this decomposition differs from the conventional Curvelet transform in dyadic partitioning, and the innermost scale is decomposed anisotropically.Using such decomposition, main features in the model are well reflected, and components in the "star" model areseparated properly. |
自适应匹配滤波的目的是使去噪后的残差最小化，而使残差的L2范数最小的自适应滤波称为最小二乘匹配滤波（LSM）.最小二乘自适应匹配滤波过程可以表示为：
|
(6) |
其中d为含噪地震记录，m为预测的噪声干扰模型，G为自适应滤波因子.LSM在时间空间域得到了广泛的应用，但是该算法要求有效信号与噪声互不相关，而实际地震数据很难符合这个要求.以面波为例，一般情况下面波干扰分布范围广，能量强，与有效信号重叠严重，因此如果直接在时间空间域进行最小二乘匹配，会损伤有效信号并导致面波能量残留.而非二次幂Curvelet变换的高度稀疏性和区分信噪的能力使得基于这种变换的最小二乘匹配算法能在充分压制噪声的同时，有效保留反射信号，并大幅提高算法的稳定性和准确性.
NDCLSM的算法流程可以归纳为以下几点：通过分析输入数据和噪声干扰的特征，对记录及预测的噪声干扰进行非二次幂Curvelet变换；然后分析噪声干扰能量主要集中的子域，并实施最小二乘意义下的匹配算法，最小化噪声干扰和实际资料中噪声干扰所对应的非二次幂Curvelet系数的残差能量.
在非二次幂Curvelet域的最小二乘匹配滤波可以表示为：
|
(7) |
其中dc为实际地震记录对应的非二次幂Curvelet系数，mc为预测的噪声干扰记录对应的非二次幂Curvelet系数，Gc为将预测噪声干扰对应的非二次幂Curvelet系数进行旋转、尺度化和平移化的变换或滤波矩阵.求解（7）式残差的L2范数最小可以得到最小二乘意义下的最佳滤波因子Gc.为了保证算法的稳定性，我们引入阻尼因子ε2，并同时使能量残差和非二次幂Curvelet系数达到最小，即
|
(8) |
上述目标函数可以使用共轭梯度等最优化方法求解.在噪声干扰能量集中的非二次幂Curvelet子记录上使用最小二乘算法得到对应的滤波因子
在实际地震资料噪声处理中，NDCLSM算法可广泛用于诸如面波、表面多次波和层间多次波等的自适应匹配相减，本文将该方法应用于面波提取.面波是陆上地震勘探中最常见的噪声干扰之一，特别是在中国西部山区、高原地区、沙漠地区等广泛存在.面波能量强，大范围地覆盖在有效信号之上，通常存在频散现象，严重影响着后续的处理和解释工作.传统的去噪方法主要利用面波频率低、视速度小的特点，很难做到在消除面波的同时保护有效信号，并且难以保存对后续地质构造研究和解释极为重要的低频数据.此外，在天然地震学研究中，有效提取面波成分对反演S波速度、解释地壳和上地幔结构具有重要意义.
为了测试本文方法在面波提取中的效果，我们对实际三维地震资料进行了处理.图 4a为某地区三维地震数据中的一个单炮记录，可以看出该数据中面波能量非常强，覆盖在有效反射信号同相轴之上，使得反射信号中断、模糊甚至畸变，这种数据无法直接用于后续的处理和解释.我们提取数据中低频低视速度的成分作为初始面波模型，如图 4b所示.图 4c为本文方法压制面波干扰后的结果.可以看出面波成分压制得比较彻底，保证了反射信号的完整性，使隐藏在面波能量下的反射信号清晰可见，并与周围的同相轴较为自然地连为一体.
|
图 4 对一个共炮点道集记录的处理结果 （a）一个带有面波干扰的炮集；（b）预测的面波模型；（c）应用本文方法得到的处理结果. Fig. 4 Processing result of a common shot record (a)Shoot gather with surface waves interferences; (b)Predicted ground-roll model; (c)Result using the proposed filtering method. |
由于NDCLSM算法效率较高，我们可以实现对大量数据的快速处理.图 5a为一个三维数据体的部分显示，我们对该数据体的单炮记录重复图 4中的去噪过程.可以看出，处理前面波几乎覆盖了整个数据体，有效信号难以识别，处理后的结果如图 5b所示，面波能量得到了有效去除，有效反射信号变得清晰、自然而且连续，表明了本文方法在处理大数据量实际资料中的有效性.由于面波能量随着距离的增大而迅速衰减，因此其能量主要集中在近偏移距道集.为了更清楚的显示去噪效果，我们将数据按共偏移距道集排列，并显示近偏移距部分，如图 6a所示，图 6b为图 6a对应的处理结果.图 6c-6d为图 6a中矩形区域处理前和处理后的显示结果.
|
图 5 三维实际资料的处理结果 （a）处理前的记录；（b）应用本文方法处理后的结果. Fig. 5 Processing results of a 3-D field data (a)Records before filtering; (b)Result after filtering by the proposed method. |
|
图 6 对三维数据体共偏移距道集近偏移距部分的处理结果 (a)处理前记录; (b)应用本文方法处理后的结果; (c)(a)中矩形区域的放大显示; (d)对(c)进行去噪处理的结果. Fig. 6 Processing results of near offset parts of 3D common offset gather (a) Data before filtering; (b) Result filtering by the proposed method; (c) Zoom view of the rectangular part in (a); (d) Result of de-noise processing to (c). |
通过分析地震数据的频率和方向特征，本文对地震数据进行自适应非二次幂Curvelet分解，最大程度地将噪声干扰和有效反射信号分离开来，在此基础上，结合最小二乘匹配滤波算法，实现了自适应信噪分离，提高了常规最小二乘匹配算法的稳定性和匹配结果的准确度.三维实际数据的处理结果验证了本文方法的可行性和有效性.
致谢感谢两位匿名评审人提出的宝贵意见.
| [1] | Candès E J, Donoho D L. Curvelets:A Surprisingly Effective Nonadaptive Representation for Objects with Edges, Curves and Surfaces. TN: Vanderbilt University Press, 1999 . |
| [2] | Candès E J, Guo F. New multiscale transforms, minimum total variation synthesis:Applications to edge-preserving image reconstruction. Signal Processing , 2002, 82(11): 1519-1543. DOI:10.1016/S0165-1684(02)00300-6 |
| [3] | Starck J, Murtagh F, Candes E, Donoho D. Gray and color image contrast enhancement by the Curvelet transform. IEEE Transactions on Image Processing , 12: 706-717. DOI:10.1109/TIP.2003.813140 |
| [4] | Candès E J, Donoho D L. Curvelets, multiresolution representation, and scaling laws. Wavelet Applications in Signal and Image Processing VⅢ, A. Aldroubi, A. F. Laine, M. A. Unser eds., Proc. SPIE 4119. |
| [5] | Starck J L, Candès E J, Donoho D L. The Curvelet transform for image denoising. IEEE Trans. Image Processing , 2002, 11(6): 670-684. DOI:10.1109/TIP.2002.1014998 |
| [6] | Herrmann F J, Verschuur E. Curvelet-domain multiple elimination with sparseness constraints. Society of Exploration Geophysicists, Expanded Abstracts , 2004, 23(1): 1333-1336. |
| [7] | Hennenfent G, Herrmann F J. Three-term amplitude-versus-offset (AVO) inverse on revisited by Curvelet and wavelet transforms. Society of Exploration Geophysicists, Expanded Abstracts , 2004, 23: 211-214. |
| [8] | Chauris H, Nguyen T. Seismic demigration/migration in the Curvelet domain. Geophysics , 2008, 73(2): S35-S46. DOI:10.1190/1.2831933 |
| [9] | Douma H, de Hoop M. Leading-order seismic imaging using Curvelets. Geophysics , 2007, 72(6): S231-S248. DOI:10.1190/1.2785047 |
| [10] | Hermann F J, Wang D, Hennenfent G, et al. Curvelet-based seismic data processing:A multiscale and nonlinear approach. Geophysics , 2008, 73(1): 706-716. |
| [11] | Herrmann F G, Hennenfent G. Non-parametric seismic data recovery with Curvelet frames. Geophys. J. Int. , 2008, 173(1): 233-248. DOI:10.1111/gji.2008.173.issue-1 |
| [12] | Naghizadeh M, Sacchi M D. Beyond alias hierarchical scale Curvelet interpolation of regularly and irregularly sampled seismic data. Geophysics , 2010, 75(6): WB189-WB202. DOI:10.1190/1.3509468 |
| [13] | 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 |
| [14] | 袁艳华, 王一博, 刘伊克, 常旭. 非二次幂Curvelet变换及其在地震噪声压制中的应用. 地球物理学报 , 2013, 56(3): 1023–1032. Yuan Y H, Wang Y B, Liu Y K, Chang X. Non-dyadic Curvelet transform and its application in seismic noise elimination. Chinese J. Geophys. (in Chinese) , 2013, 56(3): 1023-1032. |
| [15] | Wang Y B, Yuan Y H, Liu Y K, Chang X. Geophysical data oriented Curvelet-like transform. SEG Abstract , 2011: 3663-3667. |
| [16] | Yuan Y H, Wang Y B, Liu Y K, Chang X. Curvelet domain adaptive least-squares subtraction of internal multiple. Journal of Seismic Exploration , 2011, 20(3): 273-288. |
| [17] | 包乾宗, 高静怀, 陈文超. 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) (in Chinese) , 2007, 41(6): 650-654. |
| [18] | Wang Y B, Xue Y W, Dong S Q. Surface waves suppression using interferometric prediction and Curvelet domain hybrid L1/L2 norm subtraction. SEG Expanded Abstracts , 2009, 28: 3292. |
| [19] | 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 |
| [20] | Spitz S. Pattern recognition, spatial predictability and subtraction of multiple events. The Leading Edge , 1999, 18(1): 55-58. DOI:10.1190/1.1438154 |
| [21] | Guitton A. Multiple attenuation with multidimensional prediction-error filter. 73th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts , 2003: 57-74. |
| [22] | 金德刚, 常旭, 刘伊克. 逆散射级数法预测层间多次波的算法改进及其策略. 地球物理学报 , 2008, 51(4): 1209–1217. Jin D G, Chang X, Liu Y K. Algorithm improvement and strategy of internal multiples prediction based on inverse scattering series method. Chinese J. Geophys. (in Chinese) , 2008, 51(4): 1209-1217. |
| [23] | Lin D, Young D, Huang Y, et al. 3-D SRME application in the Gulf of Mexico. 74th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded abstracts , 2004, 23: 1257-1260. |
| [24] | Berkhout A J, Verschuur D J. Estimation of multiple scattering by iterative inversion, part I:theoretical considerations. Geophysics , 1997, 62(5): 1586-1595. DOI:10.1190/1.1444261 |
| [25] | Weglein A B, Carvalho F A, Stolt P M. An inverse scattering series method for attenuating multiples in seismic reflection data. Geophysics , 1997, 62(6): 1975-1989. DOI:10.1190/1.1444298 |