2. University of Alberta, Edmonton T6G 2G7
2. University of Alberta, Edmonton, T6G 2G7, Canada
1 引言
多次波的压制是地震数据处理人员面临的一个很重要的问题,抛物拉东变换是一种应用广泛的有效衰减多次波的基本方法(Yilmaz et al., 1994;Sacchi et al., 1995;熊登等,2009;郭梦秋等,2012).Thorson(1985),Hampson(1986)和Kostov(1990)论述了用最小平方方法进行拉东变换的求解,并指出如何减轻由有限孔径引起的拖尾效应.为避免反演大的矩阵,拉东变换常在频率-空间(f-x)域进行求解.另外,拉东变换也可在时间域进行求解,Thorson(1985)和Yilmaz(1994)提出了在时间域进行的高分辨率最小平方倾斜叠加因子.最小平方拉东变换方法可获得能量最小化的解,但最小平方解的分辨率还不足以区分动校正量比较小的同相轴.为克服这一缺点,有不少专家学者基于抛物拉东变换陆续提出了几种时间域和频率域的高分辨率拉东变换反演方法(Sacchi et al., 1995;Cary,1998;Herrmann et al., 2000;Hargreaves et al., 2001;Sahonewill et al., 2002;Ng et al., 2004;Wang et al., 2009),其中,Sacchi(1995)使用稀疏先验对同相轴的曲率方向进行了约束,从而在拉东域获得与输入地震数据相对应的稀疏的数据体显示.与经典最小平方拉东变换相比,高分辨率拉东变换使变换所得的拉东域数据体的能量更加集中,在同样动校正量差的情况下,高分辨率拉东变换可更好地区分有 效波和多次波.由于高分辨率反演算子是Hermitian 算子,不具有Toeplitz结构,这使得求解很耗时.针对这一问题,Sacchi(1999)给出了快速求解高分辨率抛物拉东变换的方法,利用循环矩阵和共轭梯度算法来反演Hermitian算子,与Cholesky分解相比,该方法大大提高了运算效率.这一改进使得高分辨率拉东变换的应用更加普遍,而且有利于处理大的地震数据体.
尽管不少学者提出了提高拉东域转换数据体分辨率的方法,但是,当有效波和多次波间的动校正量差较小或需处理的地震数据体具有振幅随偏移距变化(Amplitude Versus Offset,AVO)现象时,用拉东变换去除多次波时会遇到有效波振幅失真的问题,尤其是在近、远偏移距处,振幅扰动更明显(Trad et al., 2003;Sahonewill et al., 2007;Kabir et al., 1999,2007).Kabir(1999)指出由于抛物拉东变换不是正交变换,在正变换过程中会出现假象,反变换不能完全重构原始地震数据,尤其会使近偏移距处的振幅失真,进而导致AVO分析出现误差.对于抛物拉东变换在对实际资料进行多次波衰减时是否会损坏AVO振幅信息这一问题,2007年,Recorder杂志专门邀请了该领域的两位著名专家BP的 Kabir和Delft的Verschuur进行了分析讨论.Kabir(2007)指出由于拉东变换中的正、反拉东变换并不完全是彼此的逆过程,因此,在用最小平方进行拉东变换时会在近偏移距处产生假象.即使是比较先进的加权最小平方或高分辨率拉东变换也不能完全消除这些假象.Verschuur(2007)指出,与最小平方拉东变换相比,高分辨率拉东变换虽然有了明显的提高,但是,在拉东域仍然存在拖尾效应,拉东域的多次波和有效波仍存在重叠现象.虽然高分辨率拉东变换大大减小了近偏移距处的振幅泄漏,但反演所恢复的有效波振幅的可靠性仍待商榷.
为消除Radon变换对有效波振福的损失,保护原始数据的AVO现象,近年来不少学者进行了相关的研究,Wu(2008,2009)提出了基于模式识别的多次波识别和衰减技术,用以保护有效波的AVO信息.Nowak(2006)采用加权解削减近、远偏移距的假 象进而达到保护AVO信息的目的.薛亚茹等(2012)提出一种基于Radon变换和正交多项式变换的多方向正交多项式变换压制多次波的方法.从本质上来说,常规抛物拉东变换假设剩余动校正量能表示为比较好的抛物线形式,而且是在没有考虑振幅随偏移距的变化的条件下进行的多次波压制,但是实际地震资料不可能符合以上情况.本文基于常规抛物拉东变换提出了一种新的保幅的抛物拉东变换表达式,并形成基于保幅拉东变换的多次波衰减方法.该方法同时兼顾反射波形状和振幅特性的转换,可在去除多次波的同时有效保护地震数据的AVO现象.
2 基本理论拉东变换是一种应用广泛的有效衰减多次波方法,它可在频率域或时间域进行,由Cary(1998)及Bancroft和Cao(2004)的研究可知,频率域的拉东变换速度比较快,而时间域的方法则更加稳定,能更好地分离相干信号和噪音.拉东变换一般分为线性拉东变换、双曲拉东变换和抛物拉东变换三种形式.Hampson(1986)提出了时间域的经典抛物拉东变换,在该变换中,时间域的地震道集可以在拉东变换域表示为常振幅抛物轴叠加的形式.其正演算子的离散形式可表示为
其中,d(h,t)是叠前CMP道集,h是偏移距,t是在偏移距h处的时间,τ是在零偏移距处的截距时 间,p通常指远偏移距处的剩余校正时间量.m(τ,p)表示拉东域的转换数据体.
方程(1)可把拉东域的一个点转换成时间-偏移距域的一条抛物线,它只是形状上的转换,没有考虑振幅的变化,因此,方程(1)不能模拟振幅随偏移距的变化,即地震数据的AVO特性,而AVO现象是进行地下目的储层准确描述的一个重要特征.
考虑到地震数据的AVO现象,本文给出了保幅拉东变换(AVO-Preserving Radon Transform,AVO-P RT)的正演表达式:
其中,函数F(h)表示振幅随偏移距的近似变化规律.在上式中,时间域的地震数据通过拉东变换可用转换域的两个数据体mA(τ,p)和mB(τ,p)来表示.在本文中,采用了F(h)=(h/hmax)2,hmax指最大偏移距.
利用(2)式,通过保幅拉东变换可以从地震道集中同时估计mA(τ,p)和mB(τ,p).与(1)式相比,由于(2)式不仅考虑了同相轴的形状,也考虑了同相轴幅度随偏移距的变化,因此,利用该式可以模拟地震数据的AVO信息.
3 算法求解方程(2)可以简化为下面的矩阵表达形式:
从表达式可以看出,保幅拉东变换与常规拉东 变换的形式类似,不同之处在于常规拉东变换只 使用其中的 L A和m A.以(2)式为基础,本文采用Sacchi(1995)提出的高分辨率拉东变换方法对方程进行求解,进而得到高分辨率的保幅拉东变换方法.由于考虑了地震数据的AVO现象,与常规拉东变换相比,保幅拉东变换的正演和共轭算子均包含两项.
下式给出了求解拉东域转换数据体的目标函数J,使目标函数最小化可获得拉东域的数据体mA(τ,q)和mB(τ,q).
其中,规则化项R是对转换数据体mA(τ,p)和mB(τ,p)的稀疏约束,W(t,h)是数据空间的加权因子.如果只是反演mA(τ,p),该算法则变为常规时间域的高分辨率抛物拉东变换.平衡参数μ用于平滑数据匹配与稀疏约束间的关系.另外,该算法还用到两个标定参数sA和sB.用于稀疏约束的规则化项可以表示为R(x,y)=log(1+x2+y2)的形式,该形式与Sacchi(1995)提出的Cauchy稀疏约束类似,但此处对应的是稀疏约束组,也就是把mA(τ,p)和mB(τ,p)看作一组,同时对二者进行稀疏约束.最后可对方程(3)采用叠代重加权最小平方(IRLS)进行最小化求解. 4 模型处理与分析
本文用Ostr and er(1984)提出的第三类AVO模型合成包含有效波和多次波的地震道集,其模型数据见表 1.利用常规拉东变换和保幅拉东变换分别对在远、近偏移距处具有不同动校正差的地震道集进行多次波衰减,Δτhmax和Δτ0分别表示最大和最 小偏移距处多次波(M)和有效波(R)间的动校正量差,其参数示意图如图 1所示.
为产生多次波,假设在Ostr and er模型的气层之上有多套地层,多次波来自于其中的一个反射界面,有效波来自Ostr and er模型气层的顶界面,产生的多次波和有效波彼此相交.
基于以上假设,叠前地震道集由表 1中所示的经典AVO模型产生的地层反射系数与30 Hz的地震子波褶积产生,最小和最大偏移距分别为100 m和1000 m.图 2给出了Δτhmax和Δτ0均为0.03 s,有效波和多次波间的动校正量差比较大的合成记录.图 3给出了Δτhmax和Δτ0均为0.01 s,有效波和多次波间的动校正量差比较小的合成记录.从图中可以看出,有效波振幅随着偏移距的增加,幅度增加,具有明显的第三类AVO特征,而多次波的AVO现象则不明显,与实际的地质现象相吻合,这是因为多次波的产生常常需要经过多次反射,地层的吸收使得多次波的能量减小,振幅随偏移距变化的特性也受到削弱.
下面分别采用图 2和图 3所示的Δτ0、Δτhmax均为0.03 s和Δτ0、Δτhmax均为0.01 s,即不同动校正差的合成CMP道集对算法进行详细的分析.
对于Δτ0和Δτhmax均为0.03 s的合成道集数据,在图 2所示的初始合成记录的基础上按信号与噪音的最大能量之比添加信噪比为3:1的部分随机噪音.用常规和保幅拉东变换方法分别对含噪数据进行高分辨率拉东变换,所得拉东域数据体如图 4所示.图 4a是常规拉东变换结果,图 4b和4c是保幅拉东变换结果,其中,4b是mA,4c是mB.
对图 4a来说,虽然是高分辨率拉东变换结果,但仍存在拖尾效应,在拉东域所得的多次波和有效波的能量有部分重叠.这意味着由于时间域抛物同相轴的AVO现象,在拉东域有某一特定的“区域”而不是特定的“点”与时间域的同相轴相对应.因此,对具有AVO现象的地震同相轴,用常规拉东变换进行多次波和有效波的分离有它本身的局限性.如果拖尾效应的影响超过多次波和有效波间的分离界限,在估计的有效波上会存在假象,进而影响数据的AVO特征.而对于图 4b和4c中的保幅拉东变换结果,由于考虑了数据的振幅变化,地震同相轴的形状和幅度分别被转换到拉东域的不同数据体中,从而减小了转换数据的拖尾效应,使多次波和有效波在拉东域得到更好的分离.
由于是合成数据,可计算得到估计的有效波与 实际的有效波间的差异.为更好地对比两种变换方 法的效果,图 5给出了两种方法估计的有效波与实际有效波间的差,图 5a是常规方法所得有效波与实际有效波的差,图 5b是保幅拉东变换结果与实际有效波的差.从图中可以看出,与图 5b相比,图 5a的误差比较大,尤其在远偏移距处.这表明,常规拉东变换所估计的有效波不仅没得到完全恢复,而且还包含残留的部分多次波,而保幅拉东变换则比常规拉东变换能获得更好的振幅恢复.
虽然从图 5中可以明显看出两种方法所得结果的差异.但是需要说明的是,由于多次波和有效波间的动校正量差比较大,在拉东域进行多次波与有效波的分离时,常规拉东变换引起的拖尾效应比较小,从而对有效波的影响不大.
在Ostrander模型气层顶部界面的反射轴附近取一小时窗,计算在该时窗内,通过两种方法估计所得有效波振幅与实际有效波振幅间的均方根误差(Mean Square Error,MSE)MSE=E(PTrue-P)2,其结果如图 6所示.虚线显示的是常规高分辨率方法的结果(High Resolution Radon Transform,HR-RT),实线显示的是保幅高分辨率方法的 结果(AVO-Preserving High Resolution Radon Transform,AP-HR-RT).
从图中可以看出,在近偏移距处二者的差异比较小,但在远偏移距处二者的差异比较大.通过常规方法获得的远偏移距处振幅的能量被削弱,保幅方法很好地保护了远偏移距处的振幅信息.这是因为在远偏移距处振幅随偏移距的变化比较明显,常规方法因为没考虑振幅的变化,从而不能很好地恢复远偏移距处的有效波振幅,而保幅变换考虑了数据 的AVO效应,从而减小了远偏移距处振幅的能量泄漏,使有效波得到很好的恢复.
4.2 动校正量差较小的合成CMP道集分析对Δτ0和Δτhmax均为0.01 s的地震道集数据,其合成记录如图 3所示,从图中可清楚看出,由于多次波和有效波间的动校正量之差比较小,二者非常接近,这将使二者的分离比较困难.同样在该初始合成记录上按信号与噪音的最大能量之比添加信噪比为3 ∶ 1的部分随机噪音.两种拉东变换方法对加噪地震数据所得拉东域数据体如图 7所示.图 7a是常规拉东变换的结果,图 7b和7c是保幅拉东变换的结果,其中,7b是mA,7c是mB.从图上可清楚看出保幅拉东变换方法相对常规方法的改进之处.常规拉东变换没有考虑地震数据的AVO效应,由于拖尾效应,在拉东域所得的多次波和有效波的能量相对比较分散,而且有部分重叠.保幅拉东变换考虑了数据的振幅变化,地震同相轴的形状和幅度分别被转换到拉东域的不同数据体中,更好地满足了抛物拉东变换的条件,从而减小了转换数据的拖尾效应,使多次波和有效波在拉东域得到更好的分离.
图 8显示了通过两种方法恢复的有效波与实际有效波间的差异.图 8a是常规方法所得结果,图 8b是保幅拉东变换所得结果.可以看出,与图 8b相比,图 8a误差比较大,尤其在近、远偏移距处.这是由于有效波和多次波间的动校正量差比较小,振幅变化引起的拖尾效应使得多次波和有效波在拉东域不能得到较好的分离,从而使所恢复的有效波的振幅误差比较大.相对来说,由于保幅拉东变换考虑了数据的AVO特性,减小了拉东域数据的拖尾效应,尽管动校正量差比较小,仍能使多次波和有效波进行分离,从而比常规拉东变换获得了更好的振幅恢复.
同样的,对小动校正量差的情况,也在Ostr and er 模型气层顶部界面的反射轴附近取一小时窗,然后在该时窗内计算两种方法估计所得有效波振幅与实际有效波振幅间的均方根误差,其结果如图 9所示.虚线显示的是常规高分辨率方法的结果,实线显示的是保幅高分辨率方法的结果.从图中可以看出,保幅拉东变换的MSE误差比常规方法的误差小,特别在近偏移距和远偏移距,常规方法的误差尤其大.
对小的动校正量差的情况,通过常规拉东变换可衰减多次波的能量,但在近偏移距处,存在振幅能量的减小,而在远偏移距处,恢复的有效波又包含部分多次波,即存在多次波的剩余.这是因为在近偏移距,由于有效波和多次波在拉东域比较难分离,使得部分有效波的能量映射到多次波区域而被衰减,使得恢复的有效波存在能量泄漏.在远偏移距处,部分多次波的能量被映射到有效波区域,使得对多次波进行衰减后仍存在多次波能量的剩余.相对来说,由于保幅拉东变换考虑了地震数据的AVO效应,无 论在近偏移距还是在远偏移距,估计所得的有效波 的振幅都相对比较准确.
总之,以上模型测试结果表明,通过常规拉东变换可衰减多次波,但伴随有效波振幅能量的泄漏或多次波能量的剩余;与常规拉东变换相比,保幅拉东变换减小了转换数据体在拉东域的拖尾效应,可在不破坏有效波振幅的情况下衰减多次波,从而可更好地恢复地震数据的有效波信息.
5 实例分析为验证保幅拉东变换的实际应用效果,本文还对某实际工区动校正后的CMP道集数据进行了应用分析.图 10显示了色标范围一致的地震数据和通过两种拉东变换方法所恢复的有效波.其中,图 10a为叠前地震道集数据,选取的时间范围为1.6 s到4.8 s,目的储层在2.7 s附近,图 10b和图 10c分别为通过高分辨率常规拉东变换和高分辨率保幅拉东变换反演所得的有效波数据.从图 10a中可以看出,有效波基本上被校平,多次波有明显的剩余校正量,而且多次波并不是振幅为常数的完美抛物轴.从图 10b和图 10c所示的有效波数据上可以看出,多次波基本被消除,有效波能量得到了较好的恢复,但在目的储层所在的2.7 s附近,即在图中红色箭头指示位置处,二者有比较明显的差异.
为更好地展示二者间的差异,本文对图 10中红色虚线矩形框内的数据进行了放大显示,其结果如图 11所示.图 11a为放大后的地震道集、图 11b和11c分别为常规拉东变换和保幅拉东变换反演所得的有效波.从图 11b中可看到,在红色箭头所示区域,常规拉东变换所得有效波的同相轴发生了明显的扭曲,破坏了地震数据的AVO特征,进而会影响目的储层的识别,而且在蓝色箭头所示区域有明显的多次波剩余;在图 11c所示的保幅拉东变换所得结果上,红色箭头所示区域的有效波得到很好的恢复,有效保护了地震数据的AVO现象,为目的储层的识别奠定了基础,而且蓝色箭头所示区域的多次波也得到了很好的去除.
图 12展示了色标范围一致的、通过两种拉东变换方法估算得到的多次波.图 12a为高分辨率常规拉东变换所得多次波,图 12b为高分辨率保幅拉东变换所得结果.从图中可以看出,与常规拉东变换相比,保幅拉东变换可有效去除更多的多次波信息,从而可更加准确地恢复有效波的能量.
本文提出了一种通过拉东变换模拟具有AVO特征地震反射波数据的新算法.该算法是在常规拉东变换基础上进行的修改,它把地震道集转换成拉东域的两部分,其中的一部分与常规拉东变换类似,模拟了零偏移距处的反射振幅,另一部分模拟了反射振幅AVO特征.通过IRLS算法可对新的拉东算子进行反演,当然也可采用其他的稀疏反演算法.与常规拉东变换相比,无论是合成地震道集还是实际地震数据,保幅拉东变换均有比较明显的改进.当然,对于用拉东变换进行地震数据的AVO分析问题,仍需要进行进一步的研究.总之,保幅拉东变换因考虑了地震数据的AVO现象,有助于保护有效波的振幅信息,可使恢复的有效波数据体更准确地反映地下地质体信息,为后续的反演解释工作奠定有利的基础.
[1] | Bancroft J, Cao Z H.2004.Multiple attenuation using the space-time radon transform and equivalent offset gathers.74th SEG Extended Abstracts. |
[2] | Cary P W.1998.The simplest discrete radon transform.68th.Ann.Internat.Mtg., Soc.Expl.Geophys., SEG.Expanded Abstracts, 1999-2002. |
[3] | Guo M Q, Zhao Y L, Zuo S J, et al.2012.Combined multiple attenuation in marine seismic data processing.Oil Geophysical Prospecting (in Chinese), 47(4): 537-544. |
[4] | Hampson D.1986.Inverse velocity stacking for multiple elimination.Journal of the CSEG, 22(1): 44-55. |
[5] | Hargreaves N, Cooper N, Whiting P.2001.High-resolution Radon demultiple.70th Ann.Internat Mtg., Soc.Expi.Geophys..Expanded Abstracts, 1325-1328. |
[6] | Herrmann P, Mojesky T, Magesan M, et al.2000.De-aliased, High-Resolution Radon transforms.70th Ann.Internat Mtg., Soc.Expi.Geophys..Expanded Abstracts. |
[7] | Kabir N M M, Marfurt K J.1999.Toward true amplitude multiple removal.The Leading Edge, 18(1): 66-73. |
[8] | Kabir N M M, Verschuur E.2007.Does parabolic Radon transform multiple removal hurt amplitudes for AVO analysis? CSEG Recorder, 10-14. |
[9] | Kostov C.1990.Toeplitz structure in slant-stack inversion.60th Ann.Internat Mtg., Soc.Expi.Geophys..Expanded Abstracts, 1618-1621. |
[10] | Ng M, Perz M.2004.High resolution Radon transform in the t-x domain using "intelligent" prioritization of the Gauss-Seidel estimation sequence.74th Ann.Internat Mtg., Soc.Expi.Geophys..Expanded Abstracts, 2160-2163. |
[11] | Nowak J E, Imhof G M.2006.Amplitude preservation of Radon-based multiple-removal filters.Geophysics, 71(5): V123-V126. |
[12] | Ostrander W J.1984.Plane-wave reflection coefficients for gas sands at non-normal angles of incidence.Geophysics, 49(10): 1637-1648. |
[13] | Sacchi M, Ulrych T.1995.High-resolution velocity gathers and offset space reconstruction.Geophysics, 60(4): 1169-1177. |
[14] | Sacchi M.1999.Fast high resolution parabolic Radon transform.69th Ann.Internat Mtg., Soc.Expi.Geophys..Expanded Abstracts, |
[15] | Sahonewill M, Zwartjes P.2002.High-resolution transforms and amplitude preservation.EAGE 64th Conference and Exhibition. |
[16] | Sahonewill M, Aaron P.2007.Applications of time-domain high-resolution radon demultiple.77th Ann.Internat Mtg., Soc.Expi.Geophys..Expanded Abstracts, 2565-2569. |
[17] | Thorson J R, Claerbout J F.1985.Velocity-stack and slant stochastic inversion.Geophysics, 50(12): 2727-2741. |
[18] | Trad D, Ulrych T, Sacchi M.2003.Latest views of the sparse radon transform.Geophysics, 68(1): 386-399. |
[19] | Wang J F, Ng M.2009.Greedy Least-Squares and its application in Radon Transforms.CSEG National Convention, 5-8. |
[20] | Wu W J.2008.A new method of multiple attenuation: multiple identification and subtraction.CSEG National Convention, 1-4. |
[21] | Wu W J.2009.A new method of multiple attenuation: multiple identification and subtraction (II).CSEG National Convention, 21-24. |
[22] | Xiong D, Zhao W, Zhang J F.2009.Hybrid-domain high-resolution parabolic Radon transform and its application to demultiple.Chinese J.Geophys. (in Chinese), 52(4): 1068-1077. |
[23] | Xue Y R, Chen X H, Ma J T.2012.Multiples attenuation based on directional orthogonal polynomial transform.Chinese J.Geophys. (in Chinese), 55(10): 3450-3458. |
[24] | Yilmaz O, Taner M T.1994.Discrete plane-wave decomposition by least-mean-square-error method.Geophysics, 59(6): 973-982. |
[25] | 郭梦秋, 赵彦良, 左胜杰等.2012.海上地震资料处理中的组合压制多次波技术.石油地球物理勘探, 47(4): 537-544. |
[26] | 熊登, 赵伟, 张剑锋等.2009.混合域高分辨率抛物Radon变换及在衰减多次波中的应用.地球物理学报, 52(4): 1068-1077. |
[27] | 薛亚茹, 陈小宏, 马继涛.2012.多方向正交多项式变换压制多次波.地球物理学报, 55(10): 3450-3458. |