多次波是地震资料处理中的一个重要环节.目前压制多次波的方法可分为两类[1]:一类是基于有效波和多次波之间差异的滤波方法,例如Radon变换[2-4]和F-K 滤波等.另一类是基于波动理论的预测相减法[5-8],有波场外推法、预测减去法和逆散射级数法.抛物Radon变换根据一次波与多次波的剩余时差差异滤除多次波[2],该方法简单易实现,在实际工作中得到广泛应用.
抛物Radon变换的实质是不同曲率方向各地震道的振幅叠加,由于这一变换不是正交变换,因此反变换后不能完全重构原数据.为了减少重构误差,Hampson[2]建议首先定义抛物Radon 变换的反变换,然后应用最小二乘方法求得地震资料的抛物Radon变换.但是由于实际地震资料是离散的并且是有限孔径的,抛物Radon变换并不能将实际资料的同相轴聚焦为一点,而产生剪刀状的发散,降低了同相轴的分辨率,滤波后重构的数据发生失真,而失真的位置恰恰发生在近偏移距处,不利于AVO 的分析[9].Thorson等[10]提出了零阶正则化时域高分辨率Randon 变换;Sacchi和Ulrych[11]提出一种基于最小熵的非线性稀疏反演,这些方法提高了同相轴Radon域分辨率.Trad[12]实现了时间域的高分辨率Radon 变换,较频率域实现分辨率更高.但是当同相轴较近和AVO 较明显时,稀疏反演不能取得较好的效果[11, 13-14].上述高分辨率抛物Radon变换在一定程度上提高了一次波和多次波的分辨率和保幅特性.但是如果稀疏反演的约束过强,就会对AVO 特性产生较大影响.例如,如果一个严格抛物路径的AVO 曲线映射到抛物Radon变换的一个曲率点上,可以得到最高分辨率,但是却丢失了AVO的变化特性[15],因此抛物Radon 变换分辨率和AVO 特性是一对矛盾,难以同时兼顾.同时稀疏反演在实际应用中仍存在一些问题[16],例如反演参数难以设定,计算量太大,会引入一些假象等等.与上述稀疏反演获得高分辨率的思路不同,基于Gauss-Seidel方法的Radon 变换[17] 和局部Radon 变换[18-19]通过迭代方法,逐步获取最大能量的剩余时差参数,从而获得高分辨率的Radon 变换.但是该方法在AVO 复杂的情况下,分辨率会下降.
基于动校正CDP 道集的地震信号在横向上保持一定连续性的前提下,李庆忠提出了剔除拟合法求取P波正入射剖面[20],该方法克服了多次波并且能够保留AVO 现象,获得AVO 参数;Johansen等[21]利用单位正交多项式拟合CDP 道集的AVO特性,表明AVO 可以用较少的低阶正交多项式系数表征,滤除高阶系数可以有效压制随机噪声;压制多次波的正交多项式谱减法根据一次波和多次波正交多项式谱的分布特征,可以有效分离与一次波混叠的多次波,并保留AVO 特性[22].上述三种方法均是基于振幅随炮检距的横向连续性,通过滤除AVO振幅高频变化的成分,达到压制多次波的目的.
抛物Radon变换压制多次波利用了一次波和多次波的速度差异,其变换仅仅是某一方向的简单叠加,不利于保留振幅特性;正交多项式变换仅仅提取了AVO 的多项式特征,没有考虑速度因素.本文提出一种将一次波、多次波速度差异特性和AVO横向连续特性结合在一起的压制多次波方法.该方法以正交多项式变换和抛物Radon变换为基础,对CDP道集实现多方向的正交多项式变换,不同方向同相轴用一个曲率参数表征其速度特性,提高一次波和多次波的速度分辨率;利用正交多项式谱表征其AVO,保留其振幅特性.该变换将二维的CDP 数据映射到三维空间,导致了信息冗余,根据同相轴的能量分布特征设计相应的滤波器可以将冗余信息滤除;最后通过正交多项式反变换和Radon反变换获得压制多次波的CDP道集.
2 多方向正交多项式变换 2.1 正交多项式变换通常AVO 曲线可以用偶数阶的多项式表征为
(1) |
这里,d(x)表示地震信号振幅随炮检距x的变化,{cj,j=1,2,…}是多项式的各阶系数.通过对地震数据的拟合,可以获得各阶拟合系数{cj,j=1,2,…}.由于采用常规代数多项式拟合,各阶系数的计算相互依赖,一旦拟合阶次发生变化,所有的系数必须重新计算,因此该拟合方法的计算较为复杂.
Johansen等[21]提出利用单位正交多项式表征AVO,即
(2) |
其中{Pj(x),j=0,1,…,N}是在炮检距坐标x上构造的单位正交多项式集,C(j)是地震数据d(x)的第j阶正交多项式的分解系数,由于拟合空间采用的正交多项式空间,各个基函数相互独立,因此C(j)仅依赖于该阶正交多项式,并有
(3) |
C(j)称为数据d(x)的正交多项式系数谱(Spectrum of orthogonal polynomial coefficients).
因为{Pj(x),j=0,1,…,N}是单位正交多项式集,该变换属于正交变换.根据Pasval定理,空间域的数据能量等于正交多项式域能量,即正交多项式系数谱可以表征数据d(x)的能量,满足
(4) |
Johansen研究表明AVO 特性用较少的几个低阶正交多项式系数即可表征,滤除高阶系数可以压制随机噪声.对于动校正CDP道集,一次波被校平,其振幅特性可以用几个低阶正交多项式系数表征,而多次波由于仍存在剩余时差,横向提取时表现为叠加在一次波上的一些大的振幅点子[16],因此其正交多项式谱表征向高阶延拓,滤除高阶系数就可以部分压制多次波[18].
2.2 抛物Radon变换的基本原理抛物Radon变换定义为
(5) |
其实质是将不同曲率q方向的振幅叠加,反映了同相轴的零阶特性.τ 为截距.由于采集系统的孔径效应,实际同相轴不能聚焦在一个曲率参数上,而是发生剪刀状的发散,降低了各同相轴的分辨率,其中近偏移距影响水平方向的扩散,远偏移距影响对角线的扩散.
反变换定义为
(6) |
由于Radon 变换分辨率降低,一次波和多次波在Radon域相互重叠,切除过多可以有效压制多次波,但会损伤一次波;切除太少就会残留多次波,因此在一次波和多次波在Radon 域的分辨率影响到多次波的压制效果.为提高Radon变换分辨率和保幅的效果,一些学者提出高分辨率Radon 变换,即采用稀疏反演方法提高一次波和多次波在Radon 域的分辨率,但是Radon域分辨率和AVO 保留不能得到同时兼顾.
2.3 多方向正交多项式变换正交多项式变换反映了AVO 曲线的各阶正交多项式谱特性,描述了曲线的振幅变化特点,对于二维地震剖面没有考虑其方向因素;Radon 变换对不同方向的同相轴实现振幅叠加,仅得到了不同方向各道振幅的零阶特性,没有对高阶多项式特性进行研究,因此将正交多项式变换的振幅特性和Radon变换的方向特性结合可得到多方向的正交多项式变换.
由式(3)可知,正交多项式变换是AVO 振幅与各阶正交多项式的内积运算,同理,通过不同方向各道振幅与不同阶正交多项式的内积运算可以得到多方向正交多项式变换,即
(7) |
该式表示CDP道集在某一截距τ,沿曲率q的路径上各道振幅不同阶的正交多项式分布特征,即多方向正交多项式变换,该方法综合考虑了正交多项式变换的振幅描述与Radon 变换的方向特性.相对于Radon变换,不仅研究了AVO 的零阶特性,还考虑了AVO 的高阶多项式特性.
正交多项式变换空间是完备的正交空间,因此某一同相轴仅用其正交多项式谱就可以完全重构,即一个同相轴只需一个曲率参数和其正交多项式谱就可表征;对于不同剩余时差的同相轴,可以根据他们的q参数加以区别.该方法将同相轴分辨率提高到了最大,有利于不同同相轴的识别.
3 同相轴的剩余时差参数q的提取多方向正交多项式变换将二维的CDP 数据d(x,t)映射到了三维空间m(τ,q,j),仅用一个曲率参数可以表征同相轴,大大提高了速度分辨能力.同Radon 变换一样,由于采集系统的有限性,一个同相轴的各阶多项式谱并不仅仅映射到一个曲率参数,仍然存在剪刀状发散,因此如何有效提取同相轴的剩余时差参数是重构一次波的重要因素.
在这里采用能量谱识别同相轴.由式(4)和(7)得到地震道集在不同曲率上的能量分布为
(8) |
当采样的剩余时差与同相轴剩余时差一致时,各道的波形没有相位差,能量的叠加达到最强;如果剩余时差采样与同相轴剩余时差不一致,则各道波形存在相位差,叠加能量削弱,因此利用Radon 域的能量谱局部极大值可以识别同相轴.假设识别同相轴的有效剩余时差子空间是{q′i,i=1,2,…,M},与正变换同样道理,由式(2)和式(7)可以得到多方向正交多项式反变换:
(9) |
为验证算法的有效性,利用褶积模型生成一个复杂的合成记录如图 1a所示,包含12个一次波(图 1b)和8个多次波(图 1c).动校正后,一次波校正过量,其曲率参数为负,同相轴上扬,而多次波校正不足,其曲率参数为正;假设一次波和多次波的剩余时差满足严格的抛物线规律.该合成记录共有100道,道间距10m,最大偏移距1000m.
为了研究该方法对AVO 的影响,各个一次波和多次波的振幅随炮检距的变化用二阶多项式描述:
其中a0,a1,a2 均为实系数,分别表示AVO 函数的截距、梯度和曲率特性.12个一次波的AVO 曲线按时间顺序如图 3 蓝线所示.这些AVO 曲线表现了较强的振幅变化,有些出现了过零点.合成记录中,一次波和多次波相互交错,使得同相轴AVO 出现畸变,这将影响AVO 的参数反演.
合成记录的多方向正交多项式变换的零阶系数谱、一阶系数谱和二阶系数谱分别如图 1(d,e,f)所示.一次波和多次波由于剩余时差的不同分别分布在曲率大于零和小于零的区域.由于采集系统的有限孔径效应,各阶系数谱均出现了剪刀状的发散.近偏移距影响水平方向的发散,远偏移距影响对角线的发散,由于正交多项式变换的一阶和二阶系数计算受偏移距加权影响,因此可以看到一阶和二阶系数谱的对角线扩散尤为强烈.为了识别不同曲率的同相轴,根据式(8)利用零阶、一阶和二阶系数谱得到τ-q域的能量分布如图 2所示,当在能量谱中增加三阶系数谱能量时,能量的变化非常小,因此确定表征各同相轴AVO 特性的多项式阶次为2.由图 2可以明显地看到,每一个同相轴对应着一个能量的峰值,找到每一个极大值就获得了每一个同相轴的曲率参数,结合其正交多项式谱即可重构该同相轴.
图 1(g,h,i)分别是重构的合成记录和一次波、多次波,与理论合成记录的误差如图 1(j,k,l)所示,可以看到一次波和多次波基本得到重构,仅有较小的误差存在;其中误差主要存在于一次波和多次波近偏移距混叠部分,因此当一次波和多次波混叠时对重构的影响较大.
图 3是12条一次波的理论AVO 曲线(蓝线),合成记录的AVO 曲线(绿线)以及重构的AVO 曲线(红线),结果表明利用正交多项式谱重构的AVO曲线基本恢复了真实的AVO,尤其是一次波和多次波相互不重叠的情况.若一次波和多次波相互重叠,则拟合的结果有一定的误差.
为了验证该方法的抗噪性能,对上述合成数据叠加了随机噪声,如图 4(a,b,c)所示,信噪比-8.3dB.利用本文方法得到的合成记录,一次波和多次波剖面及其误差剖面分别如图 4(d,e,f)和(g,h,i)所示,有效数据得到了很好的恢复,信噪比提高到17.9dB,误差剖面上主要是随机噪声的能量,因此该方法具有很好的抗噪性能.
图 5a是某一动校正后的CMP 道集,多次波严重,且多次波的振幅特性较为稳定,为了更好地保留一次波振幅,在这里采用二阶多项式直接拟合多次波,而后从原始数据中减掉拟合的多次波.由于实际多次波的振幅变化并不是严格地服从多项式特征,第一次相减后会有部分多次波残留,在残留剖面上按照本文方法对多次波再次拟合、相减,经过三次迭代后,压制多次波结果如图 5b所示,可以看到大部分的多次波已经被去掉,图 5c是常规抛物Radon变换去多次波的结果.为了分析对一次波的影响,对压制多次波后的剖面与原始数据相减得到多次波的剖面如图 6所示.图 6a是本文方法得到的差剖面,仅有多次波的能量,因此利用本文方法压制多次波后对一次波损伤较小,有利于AVO 的分析;图 6b是原始数据与常规Radon 变换压制多次波后的差剖面,可以看到差剖面当中不仅含有多次波,在浅层中含有大量的一次波能量,对一次波的影响较大.
本文提出一种将一次波、多次波速度差异特性和AVO 横向连续性结合在一起的新的压制多次波方法.该方法以正交多项式变换和Radon变换为基础,对CDP道集不同方向的同相轴进行正交多项式变换,获得振幅随炮检距变换的正交多项式表征.相对于Radon变换,在AVO 零阶多项式特性的基础上进一步考虑了AVO 的高阶多项式特性,有利于保留一次波的真振幅.该变换将二维的CMP 数据映射到三维空间,导致了信息冗余,根据Radon域同相轴能量最大的特点,可以有效提取不同同相轴的剩余时差参数,并用一个曲率参数表征同相轴的剩余时差,大大提高了速度识别能力.对实际数据的处理表明,该方法可以较好地拟合多次波,同时对一次波的损伤较小,有利于进一步的AVO 分析.
[1] | Weglein A. Multiple attenuation: an overview of recent advances and the road ahead. The Leading Edge , 1999, 18(1): 40-44. DOI:10.1190/1.1438150 |
[2] | Hampson D. Inverse velocity stacking for multiple elimination. Journal of Canada Social Exploration Geophysics , 1986, 22(1): 44-55. |
[3] | Foster D J, Mosher C C. Suppression of multiple reflections using the Radon transform. Geophysics , 1992, 57(3): 386-395. DOI:10.1190/1.1443253 |
[4] | Yilmaz O, Taner T M. Discrete plane-wave decomposition by least-mean-square-error method. Geophysics , 1994, 59(6): 973-982. DOI:10.1190/1.1443657 |
[5] | Weglein A, Hsu S Y, Terenghi P, et al. Multiple attenuation: recent advances and the road ahead. The Leading Edge , 2011, 30(8): 864-875. DOI:10.1190/1.3626494 |
[6] | Morley L, Claerbout J. Predictive deconvolution in shot-receiver space. Geophysics , 1983, 48(5): 515-531. DOI:10.1190/1.1441483 |
[7] | Weglein A B, Gasparotto F A, Carvalho P M, et al. An inverse-scattering series method for attenuating multiples in seismic reflection data. Geophysics , 1997, 62(6): 1975-1989. DOI:10.1190/1.1444298 |
[8] | 黄新武, 孙春岩, 牛滨华, 等. 基于数据一致性预测与压制自由表面多次波――理论研究与试处理. 地球物理学报 , 2005, 48(1): 173–180. Huang X W, Sun C Y, Niu B H, et al. Surface-related multiple prediction and suppression based on data-consistence: A theoretical study and test. Chinese J. Geophys. (in Chinese) , 2005, 48(1): 173-180. |
[9] | Kabir N M M, Marfurt K J. Toward true amplitude multiple removal. The Leading Edge , 1999, 18(1): 66-73. DOI:10.1190/1.1438158 |
[10] | Thorson J R, Claerbout J F. Velocity-stack and slant-stack stochastic inversion. Geophysics , 1985, 50(12): 2727-2741. DOI:10.1190/1.1441893 |
[11] | Sacchi M, Ulrych T J. High-resolution velocity gathers and offset space reconstruction. Geophysics , 1995, 60(4): 1169-1177. DOI:10.1190/1.1443845 |
[12] | Trad D. Implementations and Applications of the Sparse Radon Transform. Columbia: University of British Columbia, 2001 . |
[13] | Schonewille M, Aaron P. Application of time-domain high-resolution Radon demultiple. SEG Expanded Abstract , 2007(1): 2565-2567. |
[14] | Hargreaves N, Cooper N. High-resolution Radon demultiple. SEG Expanded Abstract , 2001: 1325-1328. |
[15] | Schonewille M. High-resolution transform and amplitude preservation. SEG Expanded Abstract, 2002. |
[16] | Trad D, Ulrych T, Sacchi M. Latest views of the sparse Radon transform. Geophysics , 2003, 68(1): 386-399. DOI:10.1190/1.1543224 |
[17] | Ng M, Perz M. High resolution Radon transform in the t-x domain using "intelligent" prioritization of the Gauss-Seidel estimation sequence. SEG Expanded Abstracts , 2004: 2160-2163. |
[18] | Wang J, Ng M, Perz M. Fast high-resolution Radon transforms by greedy least-squares method. SEG Expanded Abstracts , 2009: 3128-3132. |
[19] | Wang J F, Ng M, Perz M. Seismic data interpolation by greedy local Radon transform. Geophysics , 2010, 75(6): WB225-WB234. DOI:10.1190/1.3484195 |
[20] | 李庆忠. 用剔除拟和法求纵波正入射剖面——一种取代水平叠加的处理技术. 石油地球物理勘探 , 1995, 30(2): 147–167. Li Q Z. Producing a normal-incident P-wave section by deletion fitting method: a processing technique for replacing horizontal stack. Oil Geophysical Prospecting (in Chinese) , 1995, 30(2): 147-167. |
[21] | Johansen T A, Brulan L, Lutro J. Tracking the AVO by using orthogonal polynomials. Geophysical Prospecting , 1995, 43(2): 245-261. DOI:10.1111/gpr.1995.43.issue-2 |
[22] | 薛亚茹, 陈小宏, 陆文凯. 压制多次波的正交多项式谱减法. 地球物理学报 , 2009, 52(3): 817–823. Xue Y R, Chen X H, Lu W K. Orthogonal polynomial spectrum subtraction for multiple attenuation. Chinese J. Geophys. (in Chinese) , 2009, 52(3): 817-823. |