小波域的波阻抗反演方法

1、2000年2月石油地球物理勘探第35卷第1期小波域的波阻抗反演方法王西文刘全新赵应成苏明军刘洪李幼铭摘要(中国石油天然气集团公司西北地质研究所)(中国科学院地球物理研究所)王西文,刘全新,赵应成,苏明军,刘洪,李幼铭.小波域的波阻抗反演方法.石油地球物理勘探,2000,35(1):8996广义线性反演是波阻抗反演中常用的一种方法,该方法反演精度高,但对噪声非常敏感,严重时可使反演出的波阻抗值完全失真。小波变换的方法具有十分明显的分频特性和去噪功能。本文采纳了广义线性反演和小波变换方法的优点,提出了小波域的波阻抗反演方法。通过理论模型试算和实际资料处理,结果表明,该方法在理论上是正确的,实际应用

2、效果也是好的。主题词小波变换波阻抗广义线性反演分频处理去噪ABSTRACTWangXiwen,LiuQuanxin,ZhaoYingcheng,SuMingjun,LiuHongandLiYoum-ing.Amethodforwaveimpedanceinversioninwaveletdomain.OGP,2000,35(1):8996Generalizedlinearinversionisacommonmethodinwaveimpedanceinversion.Themethodbringshighinversionaccuracy;butitisverysensitivetonoise

3、,sothatitwillmakewaveimpedancevaluesdistortedfullyinseverecase.Wavelettrans2formmethodhasobviousfrequencydivisionpropertyandnoiseeliminationfunc2tion.Anewmethodforwaveimpedanceinversioninwaveletdomainisthusdevel2.Theresultsofmodel2oped,whichhasjointadvantagesoftheabovetwomethodsingcomputationandreal

4、dataprocessingshowthatthenewmethodhasrighttheo2rybasisandbringsgoodapplicationeffect.Subjectterms:wavelettransform,waveimpedance,generalizedlinearinversion,frequency2divisionprocessing,noiseelimination引言作为一种数字信号处理手段的小波变换,自80年代问世以来13,在涉及信号分析和处理的邻域4得到迅速发展和应用。在地震资料处理中小波已被用于地震数据压缩、去噪、地震记录道内插、提高地震资料分辨率、数

5、值计算(如波场延拓、波动方程反问题等)、地震资料解WangXiwen,NorthwestGeologicalResearchInstitute,CNPC,LanzhouCity,GansuProvince,730020,China本文于1999年6月18日收到,修改稿于同年9月23日收到。© 1995-2006 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.石油地球物理勘探2000年90释、地震数据特征分析、地震资料采集中小波采样等方面。地震记录在有效波范围内可能含有多种类型的噪声,因此有效地进行信噪分离就显

6、得更加重要。我们知道,反射地震记录可用“变子波模型”来描述。根据这种模型,来自不同反射界面的反射波及多次反射波到达检波器的时间、波形及频谱也都是不同的。地震记录是这些子波的叠加。实际地震记录中还包含面波、随机噪声等干扰。单从时间域或单从频率域要区分或识别它们是很困难的。利用小波变换把地震记录在二维时间频率域(或时间尺度域)进行分解,就有可能利用多次波、面波、随机干扰波及有效信号之间在时间尺度域的差异把它们分开5,6。根据小波变换的特点进行地震资料分频处理和去噪已取得了一定的进展79。地震波阻抗反演是岩性预测的有效手段之一,其中最常用的是广义线性反演10,该方法有较高的精度,但对高频随机噪声也敏

7、感。高频噪声在广义线性反演中对解病态方程组影响很大,即使加了正则化因子11也无法完全消除高频噪声的影响。因此,基于小波变换分频特性和去噪功能及广义线性反演精度高的特点,本文提出了小波域的波阻抗反演方法。方法原理设有地震信号S(t),用褶积模型表示为S(t)=+-w(t)R(t-)d(1)式中:R(t-)为反射系数;w(t)为地震子波。式(1)的小波变换表示为10S(b,a)=(b-w,a)=w(t-aa-+-+S(t)gdt=aw(b-+,a)R()d(2)(3)gdtt,bR;a>0a式中:a是小波函数的尺度因子;b是时间因子;g是g的共轭。式(3)构造出的新aa(b-,a)是地震子波

8、w(t-)与小波函数g子波w的互相关。小波反变换公式为12S(b)=Cga02+-+-S(t)gRdtdaa(4)式中:Cg为常数;gR为g的实部(因为我们仅讨论小波反变换后的地震信号,而不aa是研究瞬时特征12,所以只取实部)。由式(1)(4)进一步可得S(b)=令式(4)可写成wca(b-(b-)=(b-CaCgR()0,a)dad,a)da(5)(6)0gS(b)=w-+c(b-)R()d(7)式(3)和式(6)就是重新构造地震子波的表达式。式(7)又可写成离散形式© 1995-2006 Tsinghua Tongfang Optical Disc Co., Ltd. All

9、rights reserved.第35卷第1期王西文等:小波域的波阻抗反演方法91NcS(b)=w(b)3j=1(b-Zj+1+ZjTj)(8)式中:S(b)为小波域分频处理后的地震子波记录;wc(b)为小波域构造出的新子波;Zj为波阻抗序列;(b-Tj)为狄拉克函数。(b-,a)和wc(b-)的构造。图中横坐标序号0是50Hz的我们用图1理论模型说明w14是小波变换展开尺度因子对应于频率5150Hz。Ricker子波,横坐标序号1(b-,a)子波的构造图1w(b-,a)。图2是根据式(6)在尺度因子根据式(3)可得到对应不同尺度因子a的子波w对应于3040Hz频率段(图中横坐标为-0.51.

10、0)构造出的子波wc(b-)。子波wc(b-)的构造充分体现了小波域的分频特征。当反演分辨率较高的波阻抗时,在小波域选择高频段来构造子波wc(b-);反之,当反演分辨率较低的波阻抗时,在小波域选择低频段来构造子波wc(b-)。而且构造子波wc(b-)所选择的频段要与式(8)中S(b)(小波分频处理的地震记录)的频段保持一致,这样才能取得较好的效果。在式(8)中今b=it,并对式(8)进行离散取样,可以导出偏导数矩阵A的矩阵元素(9)aij=5Zj式(9)的扩展表达式为aij=图2子波wc(b-)在3040Hz频率下的构造c2(it)3i-(Zj+Zj-1)(j-1)t-ci-jt2(it)3(

11、Zj+Zj+1)(10)(11)最终,地震道反演波阻抗的迭代格式为T-1TZk=Zk+1+(Ak-1Ak-1+I)Ak-1(S-Sk-1)© 1995-2006 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.石油地球物理勘探2000年92式中:下标k为迭代次数;Sk-1为由Zk-1按式(7)合成的地震记录;S为小波分频处理过的地震记录;为正则化因子;为控制迭代收敛速度的加速收敛因子(0<1)。式(1)(11)刻画了小波域的波阻抗反演方法。反演质量控制为了控制反演迭代速度和精度,采用统计参数和地震记录拟合

12、方差来控制反演质量。线性相关系数对于数据对(xi,yi),i=1,2,N,线性相关系数R定义为N(xR=Ni=1i)(yi-y)-x(12)N(xi)2-x(yi=1i)2-y、式中:xy分别是xi和yi的平均值;R-1,1。R=1为完全正相关;R=-1为完全负相关。当式(11)中Sk-1和地震记录S之间拟合的相关系数R0.95时,我们认为反演正确。学生分布概率函数13学生分布概率用A(t v)表示,用于检验两个观察分布是否有相同的值。A(t v)为有v重自由度的概率分布,若两均值相等,则统计量t(用于度量均值差)小于观测值。如果A(t v)趋于0,则两个均值差很小,两个观察分布趋于相同均值。

13、t定义为t=R1-R2t(13)式中:R为式(12)定义的相关系数;N为测量点数。t的学生概率分布定义为A(t v)=1 2vB,22-t1+v2-2dx(14)(15)限定值A(0 v)=0A( v)=1v+t2A(t v)=1-I2,2(16)式(16)为A(t v)与不完全B函数Ix(a,b)的关系式。可用式(14)计算学生分布概率(Prob),其值越小相关性就越好。z分布当相关性非常显著时,为了观察这种显著性差异,用下面变换评价相关性z=ln21-R(17)式中R为式(12)定义的相关系数。上述三个统计量用于控制迭代反演。地震记录拟合方差地震记录拟合方差定义为© 1995-2

14、006 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.第35卷第1期王西文等:小波域的波阻抗反演方法93NS=(Si-Si,k-1)2N(18)式中:Si为小波域处理后的地震记录元素;Si,k-1为k-1次迭代合成的地震记录元素;N为采样点数。该方差主要用于判断地震记录的拟合程度。模型试算图3是采用本文提出的小波域波阻抗反演的波阻抗模型。在小波域进行分频去噪处理后得到了地震记录(图中横坐标序号2)。由图3看出,迭代反演10次(图中横坐标序号为312)就得到了非常好的反演结果。而笔者按照广义线性反演的方法在信噪比为5的情

15、况下反演迭代20次后仍然可以看出,噪声对波阻抗反演影响很大,以致无法得到正确的波阻抗值。图3小波域波阻抗反演方法的波阻抗模型在广义线性反演中为使方程的解具有稳定性,通常在式(11)中设计一个正则化因子来压制小特征值对误差的放大作用,但同时也压制了一些有用的高频信息,所以在图3中没有反演出理论模型的“角点”。尽管如此,在无约束条件下还是反演出了理论模型的基本特征。应用实例鄯善油田位于吐鲁番地区,鄯善构造处在吐鲁番拗陷台北凹陷鄯善弧形构造带的近弧顶偏西部位,其为一个完整的穹窿背斜构造,含油面积近27km2。该油藏属构造岩性油气藏,生、储、盖、圈闭条件齐全。上侏罗统齐古组为区域性盖层,中侏罗统七克台

16、组(J2q)、三间房组(J2s)和西山窑组(J2x)为储层,其中三间房组(J2s)油藏为主力油藏。图4是一条经过精细处理的鄯善构造的地震测线(图中CDP号和时间都是相对的)。图中© 1995-2006 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.石油地球物理勘探2000年94图4Y227S地震剖面图5小波域波阻抗反演剖面© 1995-2006 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.第35卷第1期王西文等:小波域的波

17、阻抗反演方法95显示出J2q和J1b分别是二套煤系地层的强反射界面。J2s和J2x反射同相轴连续性不是太好,能量也相对较弱。图5是在无约束条件下小波域波阻抗反演剖面,该剖面的波阻抗值与岩性对应较好。从J2q到J1b波阻抗同相轴连续性增强,清楚地反映出J2q、J2s、J2x、J1b之间的层间岩性结构。图6是反演波阻抗(图5)每个地震道误差统计分布和地震记录拟合方差,可以看出:几乎所有地震道相关系数(R)趋于1;t学生分布概率(Prob)趋于0;z分布范围为15,均值约为图6反演误差统计分布和方差3;地震记录拟合方差S的范围为1560,均值约为35(方差单位与地震记录幅值单位相同,这里仅用此量作为

18、参考量来判断拟合程度)。统计分布和方差表明,反演精度是高的。结束语小波域的波阻抗反演方法利用小波变换中分频特性和去噪功能,并结合了广义线性反演方法精度高的特点,将两种方法有机地结合起来。并且可根据地质目标选择适当频段(尺度区间)进行分频反演波阻抗值。该方法特别适用于无约束条件下的波阻抗反演。© 1995-2006 Tsinghua Tongfang Optical Disc Co., Ltd. All rights reserved.石油地球物理勘探2000年96参考文献1GrossmannAandMorletJ.Decompositionoffunctionintowavelets

19、ofconstantshape,andrelatedtrans2.in.ed.StreitL,eds,LecturesonRecentResults,Singapore:WorldSci2forms,InMathematics+Physics165entifical,1985,1352GrossmannA,Kronland2MartinetR,MorletJ.Readingandunderstandingcontinuouswavelettrans2forms,In:Wavelets,Time2FrequencyMethodsandPhaseSpace,1stInt.WaveletsConf.Marseille,1989,2203DaubechI.Orthonormalbaseofwaveletswithfinitesupport2connectionwithdiscretefilters,In:141Wavelets,Time2FrequencyMethodsandPhaseSpace,1stInt.WaveletsConf.Marseille,1989,1264刘贵忠,邸双亮.小波分析及其应用,西安电子科技大学出版社,1