地球物理学报  2015, Vol. 58 Issue (10): 3783-3790   PDF    
基于复曲波变换的一次波和多次波分离方法
董烈乾, 李培明, 张奎, 汪长辉, 祝杨, 王泽    
中国石油集团东方地球物理公司国际勘探事业部, 河北涿州 072751
摘要:多次波压制方法的研究一直都是地震数据处理中非常重要的一个课题.由于常用的多次波匹配方法主要针对多次波模型和实际多次波存在的振幅或相位的差异进行匹配校正,而无法直接校正多次波模型和实际多次波存在时移误差.本文构建了一种复曲波变换的算法,利用复曲波变换的时移不变性质,通过调整复曲波系数的振幅和相位实现对多次波模型振幅和时移误差的校正.为了更好地保护有效信号,在一次波和多次波分离前,引入一个非线性屏蔽滤波器,可以事先分离出大部分有效波,然后再将剩余部分数据作为输入数据,在复曲波域进行剩余一次波和多次波分离.最后通过模型试算和实际资料处理验证了本文提出的一次波和多次波分离方法的有效性.
关键词多次波压制     复曲波变换     时移误差     屏蔽滤波器    
Primary and multiple separation method based on complex curvelet transform
DONG Lie-Qian, LI Pei-Ming, ZHANG Kui, WANG Chang-Hui, ZHU Yang, WANG Ze    
BGP international, CNPC, Hebei Zhuozhou 072751, China
Abstract: The study of de-multiple methods is a very important in seismic data processing. With increasing difficulties of oil exploration and the enhancement of seismic imaging accuracy, the existence of multiples seriously degrades the signal to noise ratio of seismic data, and interferes the identification of primary waves. Meanwhile, it may also cause some false geological features in seismic sections. Therefore, how to effectively suppress multiple is always a focused topic in seismic exploration.
The current multiple suppression methods can be classified into two categories: based on filtering methods and based on prediction-and-subtraction methods. When using wave equation and considering the properties of wave-field propagation, the prediction-and-subtraction methods are suitable for complex seismic data. However, standard matching or subtraction methods can only correct the amplitude or phase errors between the multiple model and actual multiples in the stage of subtraction. As for the misalignment errors, these methods cannot achieve better results. So, a new curvelet transform named complex curvelet transform (CCT) is proposed. Taking advantage of shift invariance property of CCT, we can use the phase and amplitude of the data's and multiple model's CCT coefficients to correct misalignment and amplitude errors between the multiple model and actual multiples. In addition, for the purpose of protecting primary waves further, a non-linear masking filter is applied in advance, which can separate most of primary waves firstly, then recover the remaining primary waves using the CCT-based separation method.
To demonstrate the validity of the CCT-based separation method, the shot gathers of a concave model are simulated. Firstly, the non-linear masking filter is applied to separate the primary and multiple with partial primary. After doing a simple L2 norm matching to the multiple, the CCT-based separation method is adopted to obtain the residual primary. Combining the residual primary and the primary separated by the non-linear masking filter, the final de-multiple data is obtained. Comparing to the standard matching methods, such as the pseudo multi-channel matching method, the CCT-based separation method can directly correct the amplitude and phase of the multiple model, and protect the primary greatly while suppressing the multiple.
The test on the field data also shows the proposed method is applicable and effective. The tests on synthetic and field data show that in the case of misalignment error existence, taking advantage of the properties of shift invariance and amplitude changing with coefficient of CCT, the phase and amplitude of the multiple model can be corrected directly to fit the actual multiple. The application of a non-linear masking filter can protect the primary better while suppressing the multiple. And the new method is proved to be applicable and effective.
Key words: Multiple suppression     Complex curvelet transform     Misalignment error     Masking filter    
1 引言

由于经过SRME多维褶积得到的多次波模型在振幅、时间等动力学特征上与实际多次波存在差异(Verschuur D J,1992),因此需要对存在的差异性进行匹配或校正后,才能够使多次波模型与实际多次波进行直接相减或分离,从而达到消除多次波的目的.目前常用的匹配相减或一次波和多次波分离的方法主要包括:基于误差能量最小的滤波方法(Wang,2003; 李学聪,2010; 李振春,2011; 董烈乾,2013a);基于独立变量分析(陆文凯,2004)和聚束滤波的方法(胡天跃,2002);基于曲波域的一次波和多次波分离方法(Herrman et al.,2008a2008b);模式匹配滤波方 法(Spitz S,1999);预测误差滤波方法(Guitton,2003; Fomel,2009; 董烈乾,2013b)等. 基于模式识别匹配滤波方法和预测误差滤波方法可以对预测的多次波模型进行较好匹配,而且容忍度也比较好,但前提是多次波同相轴具有线性和可预测性,而且当一次波和多次波重叠时,匹配效果不佳;基于误差能量最小的滤波方法是在时间域或频率域设计一个维纳滤波器,对多次波模型振幅和相位进行校正;基于独立分量的方法在一次波和多次波非正交情况下,也可以取得不错的效果;基于实曲波变换的一次波和多次波分离方法则是根据一次波和多次波在曲波域被稀疏分解到不同的尺度和方向上,通过一种非线性最优化的方法将一次波从原始数据中恢复出来.但上述方法都是针对一次波和多次波模型存在的振幅或相位的差异进行匹配校正,但对于存在时移误差的情况,上述方法获得的一次波效果不佳.因此本文构建了复曲波变换的算法,利用复曲波变换的时移振幅不变和改变复曲波系数相位可以改变同相轴相对位置的性质设计了基于复曲波变换的一次波和多次波分离方法.

2 复曲波变换的构建

对于m维地震数据,实曲波变换(RCT)可以表示为

其中,x1,…,xm表示m维的地震数据,表示尺度,R(x1,…,xm)表示实曲波变换函数,Γ,originalreal表示实曲波系数.

复曲波变换(CCT)算法和实曲波算法近似相同(Neelamani et al.,2010).首先构造信号的实曲波函数,一对具有相同倾角、位置和频率,但相互之间存在90°相移的实曲波函数构成了复曲波函数的实部和虚部,如图 1所示.则根据式(1),m维地震数据的复曲波变换可以表示为:

图 1 组成复曲波函数的两个存在90°相移的实曲波函数Fig. 1 Two real Curvelet functions with 90° phase shift that compose complex Curvelet function
其中,C(x1,…,xm)表示复曲波变换函数,Γ,originalreal表示复曲波系数.复曲波变换函数C(x1,…,xm)的实部和虚部是由两个不同实曲波函数R(x1,…,xm)构成.

3 复曲波变换的性质

通过复曲波变换的构建过程看出,复曲波变换和实曲波变换同样具有信噪分离的特性.但更为重要的一点是,复曲波变换的特殊构造克服了实曲波变换时移敏感性的缺陷.图 2是对一水平同相轴进行实曲波变换和复曲波变换的示意图.图 2a为位置在64点的一水平同相轴,2b为向下时移3个位置点(即位置在67)的同相轴,2c和2d分别为2a和2b所对应的实曲波变换系数,可以看出,曲波系数发生了很大的变化(中间部分横线颜色的深浅表示曲波系数的大小),2e和2f分别为2a和2b所对应的复曲波变换系数,可以看出,曲波系数基本没有什么变化.图 3a是对图 2a水平同相轴上下时移后,同一尺度、方向和位置的实曲波系数和复曲波系数随时移误差的变化.实线为复曲波系数随时移误差的变化曲线,下三角虚线和虚线分别表示的是实曲波系数和具有相同倾角、位置但存在90°相移的实曲波系数随时移误差的变化曲线.可以看出即使很小的时间移动,实曲波系数也会发生很大的变化,而复曲波系数基本保持一致.图 3b为地震数据复曲波系数的相位随时移误差的变化曲线,该曲线近似是线性的,可以利用这种特性,改变相位可以对同相轴进行相应的时移.

图 2 实曲波变换和复曲波变换时移敏感性对比图 (a) 原始图形; (b) 下移3个位置点后的图形; (c) 图形(a)对应的RCT系数; (d) 图形(b)对应的RCT系数; (e)图形(a)对应的CCT系数; (f) 图形(b)对应的CCT系数. Fig. 2 Time shift sensitivity contrast of RCT and CCT (a) Original data; (b) Data after 3 points downward shift; (c) RCT coefficient of figure (a); (d) RCT coefficient of figure (b); (e) CCT coefficient of figure (a); (f) CCT coefficient of figure (b).

图 3 实曲波变换和复曲波变换误差分析(a) 实曲波变换和复曲波变换系数随时移误差的变化曲线; (b) 复曲波变换相位随时移误差的变化曲线.Fig. 3 Error analysis of RCT and CCT(a) Variance curves of RCT′s and CCT′s coefficient with time shift error; (b) Variance curve of CCT′s phase with time shift error

复曲波变换实部和虚部90°的相移关系,类似于傅里叶基函数(exp(-jωt))的实部和虚部(cos(ωt)和sin(ωt))(Ramesh,2010s).傅里叶变化算法中,对某一信号进行Δt时移,这不会对信号的振幅产生任何影响,但相位会变化ωΔt.因此,假如保持信号振幅不变,相位变化ωΔt,则信号在时间域会有Δt的时移.和傅里叶变换理论类似,反射同相轴沿垂直于倾角方向小的时移会引起该反射同相轴所对应的复曲波系数相位发生变化,改变复曲波系数的相位同样也会使反射同相轴发生相应的时移.精确来说,如果复曲波系数相位改变ε,则该系数所表征的反射轴会沿振荡的方向产生ε/ω的时移量,ωC(x1,…,xm)的中心频率.图 4为复曲波系数相位改变与地震数据相对应的复曲波函数之间的变化关系.可以看出,随着相位的变化,与之相对应的地震反射沿着垂直于倾角方向相应的移动.

图 4 改变复曲波系数相位与地震数据相对应的复曲波函数的变化关系(a) π/2; (b) π/4; (c) 0; (d) -π/4; (e) -π/2.Fig. 4 Relationship between CCT coefficient phase and corresponding CCT function
4基于复曲波变换的一次波和多次波分离算法

在经过全局的振幅匹配校正后,令原始数据和模型数据复曲波系数相减后剩余误差最小化:

其中,γ*θ*分别表示最优的振幅校正因子和最优相位校正因子;Γ,originalComplexΓ,modelComplex分别表示原始数据尺度为的复曲波系数和模型数据尺度为的复曲波系数;γ表示控制模型数据复曲波系数幅值的大小;θ表示控制模型数据复曲波系数相位的大小;γ,minγγ,maxθ,minθθ,max,{γ,minγ,max}和{θ,minθ,max}分别表示给定的求取γ*θ*的变化范围.振幅和相位变换范围的选取是很重要的,直接影响最终模型数据校正后的精确度和迭代速度.通常情况下,模型数据复曲波系数的振幅误差和相位误差是可以大致预测的.在经过全局整体校正和子波匹配后,多次波模型与实际多次波时移一般在-2 ms和2 ms之间;经过简单振幅匹配后,多次波模型振幅和实际多次波振幅误差一般在30%左右.因此,根据实际情况,择优选取振幅和相位的变化范围对提高计算效率是比较重要的.

假设选取一次波模型数据进行校正,通过最小平方最优化方法求取γ*θ*后,则获得校正后的一次波数据为:

以多次波模型为输入模型数据为例,则基于复曲波变换的一次波和多次波分离步骤为:

1)对输入的原始数据和多次波模型数据进行全局的振幅匹配校正(通常进行较为粗糙的L2模匹配滤波),目的是使原始数据中多次波振幅和模型数据多次波振幅近似在一个数量级上,再对两者进行复曲波变换,得到相应的复曲波系数Γ,originalComplexΓ,modelComplex

2)在给定的{γ,minγ,max}和{θ,minθ,max}范围内,最小化式(3),通过解析法求解最优的γ*θ*,即令一未知变量γ已知,求出θ,然后再根据求解的θ,求解γ

3)将求解得到的γ*θ*反代回式(4),得到动力学特征校正后的数据.

为更有效的保护有效信号,本文引入了一种非线性屏蔽滤波器Φ,该屏蔽滤波器可以事先分离出大部分有效信号,其表达式如下:

其中,Am,表示预测多次波模型的振幅,Ao表示原始输入数据的振幅,ε表示加权因子,n表示控制滤波器平滑程度的参数.该滤波器可以将原始输入数据分为两部分:包含部分一次波能量的多次波能量Φp和完全只包含一次波能量的(1-Φ)p.在实际处理时,只是选取Φp进行处理,最后将分离的一次波和(1-Φ)p相加即得到最终多次波压制的结果,该策略可以更好的保持有效波的能量.

5 模型和实际资料试算

首先选取洼陷模型正演产生的炮记录.图 5a为原始单炮记录,图 5b为预测得到的多次波模型.应用屏蔽滤波器首先保存大部分一次波能量,如图 5c所示.再将剩余的波场作为输入的原始波场,与多次波模型进行一次简单的L2模的匹配相减处理,主要目的是使输入的原始数据与模型数据在振幅量级上相匹配;再利用复曲波变换进行剩余一次波和多次波分离,本文是直接对一次波模型进行动力学校正(原始波场减去多次波模型后的波场作为模型数据输入);振幅校正因子的范围选取的是从0.5到1.5,相位误差的校正范围为±3 ms,校正分离后的剩余一次波波场如图 5d所示,将图 5c图 5d相加,得到分离后的一次波波场,图 5f所示.与伪多道匹配方法(图 5e)相比较可以看出,基于复曲波变换的一次波分离方法能够很好地将一次波和多次波分离,同时能够保持一次波基本不受损害,而伪多道的方法仍存在较明显的多次波残余.

图 5 洼陷模型数据一次波恢复对比图(a) 原始炮记录; (b) 输入的多次波模型; (c) 应用屏蔽滤波器获得的一次波; (d) 基于复曲波变换分离的剩余一次波; (e) 伪多道方法获得的一次波; (f) 基于复曲波变换最终获得的一次波.Fig. 5 Contrast diagrams of primary recovery of sunken model data(a) Original shot gather; (b) Multiple model; (c) Preserved primary; (d) Residual primary obtained by CCT-based separation method; (e) Primary obtained by pseudo multi-channel matching method; (f) Primary obtained by CCT-based separation method.

选取某海上二维测线数据进行试算.图 6a为选取的单炮数据,可以看出数据中广泛发育了能量较强的多次波,严重地干扰了有效反射波,多次波问题突出,严重影响了后续成像及油气藏预测的可靠性.利用本文设计的一次波和多次波分离方法进行处理,该测试中振幅校正因子的范围选取的是从0到1.5,相位误差的校正范围为±3 ms.图 6e图 6f分别为应用伪多道匹配方法和本文基于复曲波变换的方法得到的一次波结果对比图.可以看出,伪多道匹配方法压制多次波后仍存在多次波残余,而基于复曲波变换的一次波分离方法能够很好地将一次波和多次波分离.图 7为抽取的近偏移距剖面,可以看出本文方法可以更好地压制中深层的多次波能量,而且可以很好地保持较微弱的有效同相轴,可以为为后续的处理提供更好的数据.

图 6 某海上数据一次波恢复对比图(a) 原始炮记录; (b) 输入的多次波模型; (c) 应用屏蔽滤波器获得的一次波; (d) 基于复曲波变换分离的剩余一次波; (e) 伪多道方法获得的一次波; (f) 基于复曲波变换最终获得的一次波. Fig. 6 Contrast diagrams of primary recovery of ocean data (a) Original shot gather; (b) Multiple model; (c) Preserved primary; (d) Residual primary obtained by CCT-based separation method; (e) Primary obtained by pseudo multi-channel matching method; (f) Primary obtained by CCT-based separation method.

图 7 某海上数据共偏移距剖面对比图 (a) 原始近偏移距数据; (b) 多次波模型数据; (c) 伪多道方法获得的一次波数据; (d) 本文方法获得的一次波数据.Fig. 7 Contrast diagrams of common-offset section of ocean data(a) Original near-offset data; (b) Multiple model; (c) Primary obtained by pseudo multi-channel matching method; (d) Primary obtained by CCT-based separation method.
6 结论

由于常规多次波匹配或分离方法仅仅是针对多次波模型和实际多次波存在的振幅或相位的差异进行校正,而对于存在时移误差的情况,多次波匹配或分离效果不佳.本文构建了一种复曲波变换的算法,并利用该算法的时移振幅不变和改变复曲波系数相位可以改变同相轴相对位置的性质,设计了基于复曲波变换的一次波和多次波分离方法.该方法通过控制模型数据复曲波系数的相位对模型数据存在的时移误差进行校正,通过调整复曲波系数的幅值实现对模型数据的振幅误差的校正.为了更好的保持有效信号,在一次波和多次波进行分离前,本文引入了一个非线性屏蔽滤波器,该滤波器可以将事先保存大部分有效波能量,然后再将剩余的部分数据作为输入数据,在复曲波域进行剩余一次波和多次波分离.模型测试和实际资料处理验证本文方法较常规匹配方法更具优越性.

致谢 感谢东方地球物理公司中青年科技创新基金项目(11-01-2015)的资助,感谢中国石油大学(华东)地震波传播与成像课题组对本文提供的帮助,感谢本文匿名评审专家提出的意见和建议.

参考文献
[1] Dong L Q, Li Z C, Yang S C, et al. 2013a. Non-causal matching filter multiple elimination method based on correlation and iteration. Chinese J. Geophys. (in Chinese), 56(10): 3542-3551, doi: 10.6038/cjg20131029.
[2] Dong L Q, Li Z C, Yang S C, et al. 2013b. An improved method of surface-related multiple elimination. Progress in Geophys. (in Chinese), 28(6): 3148-3152, doi: 10.6038/pg20130641.
[3] Fomel S. 2009. Adaptive multiple subtraction using regularized nonstationary regression. Geophysics, 74(1): V25-V33.
[4] Guitton A. 2003. Multiple attenuation with multidimensional prediction-error filters.//73th Annul International Meeting, SEG Expanded Abstracts, 1945-1948.
[5] Herrmann F J, Moghaddam P, Stolk C C. 2008a. Sparsity- and continuity-promoting seismic image recovery with curvelet frames. Applied and Computational Harmonic Analysis, 24(2): 150-173.
[6] Herrmann F J, Wang D L, Hennenfent G, et al. 2008b. Curvelet-based seismic data processing: a multiscale and nonlinear approach. Geophysics, 73(1): A1-A5.
[7] Hu T Y, Wang R Q, Wen S L. 2002. Multiple attenuation of seismic data from South China Sea by using beam-forming filtering method. OGP (in Chinese), 37(1): 18-23.
[8] Li X C, Liu Y K, Chang X, et al. 2010. The adaptive subtraction of multiple using the equipoise multichannel L1-norm matching. Chinese J. Geophys. (in Chinese), 53(4): 963-973, doi: 10.3969/j.issn.0001-5733.2010.04.021.
[9] Li Z C, Liu J H, Guo C B, et al. 2011. Amplitude-preserved multiple suppression based on expanded pseudo-milti-channel matching. OGP (in Chinese), 46(2): 207-210, 231.
[10] Lu W K, Luo Y, Zhao B, et al. 2004. Adaptive multiple wave subtraction using independent component analysis. Chinese J. Geophys. (in Chinese), 47(5): 886-891, doi: 10.3321/j.issn:0001-5733.2004.05.021.
[11] Neelamani R, Baumstein A, Ross W S. 2010. Adaptive subtraction using complex-valued curvelet transforms. Geophysics, 75(4): V51-V60.
[12] Spitz S. 1999. Pattern recognition, spatial predictability, and subtraction of multiple events. The Leading Edge, 18(1): 55-58.
[13] Verschuur D J, Berkhout A J, Wapenaar C P A. 1992. Adaptive surface-related multiple elimination. Geophysics, 57(9): 1166-1177.
[14] Wang Y H. 2003. Multiple subtraction using an expanded multichannel matching filter. Geophysics, 68(1): 346-354.
[15] 董烈乾, 李振春, 杨少春等. 2013a. 基于相关迭代的非因果匹配滤波器多次波压制方法. 地球物理学报, 56(10): 3542-3551, doi: 10.6038/cjg20131029.
[16] 董烈乾, 李振春, 杨少春等. 2013b. 一种改进表层多次波压制方法. 地球物理学进展, 28(6): 3148-3152, doi: 10.6038/pg20130641.
[17] 胡天跃, 王润秋, 温书亮. 2002. 聚束滤波法消除海上地震资料的多次波. 石油地球物理勘探, 37(1): 18-23.
[18] 李学聪, 刘伊克, 常旭等. 2010. 均衡多道1范数匹配多次波衰减的方法与应用研究. 地球物理学报, 53(4): 963-973, doi: 10.3969/j.issn.0001-5733.2010.04.021.
[19] 李振春, 刘建辉, 郭朝斌等. 2011. 基于扩展伪多道匹配的保幅型多次波压制方法. 石油地球物理勘探, 46(2): 207-210, 231.
[20] 陆文凯, 骆毅, 赵波等. 2004. 基于独立分量分析的多次波自适应相减技术. 地球物理学报, 47(5): 886-891, doi: 10.3321/j.issn:0001-5733.2004.05.021.