基于信号子空间分解的三维地震资料高分辨率处理方法

陆文凯,丁文龙,张善文等.基于信号子空间分解的三维地震资料高分辨率处理方法.地球物理学报,2005,48(4):896～901Lu W K,Ding W L ,Zhang S W ,et al.A high 2res olution processing technique for 32D seismic data based on signal sub 2space decom position.Chinese J .G eophys .(in Chinese ),2005,48(4):896～901基于信号子空间分解的三维地震资料高分辨率处理方法陆文凯1,丁文龙1,张善文2,肖焕钦2,赵铭海21清华大学自动化系信息处理所,智能技术和系统国家重点实验室,北京　1000842胜利油田有限公司,山东东营　257100摘　要　提出了一种新的基于信号子空间分解技术的32D 地震资料高分辨率处理方法.利用信号子空间分解技术,不仅可以分离32D 地震资料中的信号和噪声,而且可以进一步根据地震同相轴的倾角不同,将混合信号进行分解得到单个同相轴信号.通过丢弃噪声子空间,只对不同信号子空间重构的信号利用谱白化技术进行高分辨率处理,然后累加所有处理结果,从而达到既提高地震资料的分辨率,又提高地震资料信噪比的目的.对合成资料和实际地震资料的处理结果表明,此法具有好的应用前景.关键词　信号子空间分解,奇异值分解,谱白化,分辨率,信噪比文章编号　0001-5733(2005)04-0896-06 中图分类号　P631收稿日期　2004-12-31,2005-04-28收修定稿基金项目　国家科技攻关(2003BA613A 202),国家自然科学基金(40474040)和中国石油天然气集团公司创新基金资助.作者简介　陆文凯,男,1969年生,副教授,1991年毕业于清华大学自动化系,获学士学位,1996毕业于石油大学(北京)地球科学系,获地球物理博士学位,研究方向是信号处理及其应用.E 2mail :lwkm f @A high 2resolution processing technique for 3-D seismic databased on signal sub 2space decompositionLU Wen-K ai 1,DI NG Wen-Long 1,ZHANG Shan-Wen 2,XI AO Huan-Qin 2,ZHAO Ming -Hai 21Dept .o f Automation ,T singhua Univer sity ,K ey State Lab o f Intelligent Technology and System ,Beijing 100084,China 2Shengli Oilfield Limited Company ,Shandong Dongying 257100,ChinaAbstract This paper proposes a new signal sub 2space decom position (SS D )based technique to im prove the res olution of 32D seismic data processing.By exploiting the SS D technique ,we can decom pose not only the signal and noise in 32D seismic data ,but als o the com posite event into single events according to their different dips.A fter rem oval of the noise space ,we integrate all the reconstructions corresponding to the signal sub 2spaces with im proved res olution using the spectral whitening technique.Since we only im prove the res olution of the signal and rem ove the noise at the same time ,we can im prove the res olution and the signal-to-noise ratio (S NR )of the seismic data simultaneously.The application results of the synthetic and real data show our method has a g ood perspective of applications.K eyw ords Signal sub-space decom position ,Singular value decom position ,S pectral whitening ,Res olution ,Signal-to-noise ratio第48卷第4期2005年7月地　球　物　理　学　报CHI NESE JOURNA L OF GE OPHY SICSV ol.48,N o.4July ,20051　引　言在地震资料处理中,信噪比和分辨率是评价地震资料品质的两个重要指标.具有高信噪比和高分辨率的地震资料,可以用来进行有效的反演和地质解释,为实现精细油藏描述奠定坚实的基础.因此,如何有效提高地震资料的分辨率和信噪比,是地震资料处理中的一个关键问题.众所周知,高信噪比是提高地震资料分辨率的基础,地震资料的有效频宽制约着地震资料的分辨率的提高[1].即地震资料的分辨率和信噪比是相互矛盾的,一般情况下,提高地震资料的分辨率会导致其信噪比的降低.为了解决这个矛盾,人们往往采用先提高地震资料的分辨率再提其高信噪比,或者先提高地震资料的信噪比再提高其分辨率的策略[2].另外一类方法是利用某种变换,将地震资料进行分解,然后利用高信噪比部分约束低信噪比部分的分辨率提高,小波变换就被用来实现多分辨率地震信号反褶积[3].信号子空间分解技术被广泛地应用于信号压缩,噪声压制等领域.奇异值分解是实现信号子空间分解的一个重要手段.奇异值分解技术常常被用来提高地震资料的信噪比.与图像处理中的纹理方向检测技术相结合,奇异值分解技术可以有效地提高地震资料的信噪比[4,5],并被用来实现地震信号的插值[6].由于通过奇异值分解,原始地震资料被分解为具有不同信噪比的信号子空间,所以可以根据每个信号子空间的信噪比进行自适应的高分辨率处理,达到了既提高地震资料的分辨率,又提高其信噪比的目的[7].上述基于奇异值分解的地震信号处理方法都是二维时空域信号处理技术.随着32D地震资料的广泛应用,32D地震资料处理起的作用越来越大.利用奇异值分解对32D地震资料在f-x-y域进行分解,可以有效地提高32D地震资料的信噪比[8].本文提出的一种新的基于信号子空间分解技术的32D地震资料高分辨率处理方法,就是利用奇异值分解先在f-x-y域将32D地震资料中的信号和噪声进行分离,并将不同倾角同相轴重叠形成的复合同相轴按倾角进行分解.由于地震信号通过奇异值分解后,每个信号子空间对应的奇异值大小代表了该子空间的信噪比高低,奇异值大的子空间,其信噪比高.新方法利用谱白化技术根据每个子空间的信噪比的高低,对各个子空间重构的信号分别进行高分辨率处理,自适应地加强信号的高低频成分.最后通过丢弃噪声子空间,只累加所有信号子空间的处理结果重建整个32D地震资料,从而达到既提高地震资料的分辨率,又提高地震资料的信噪比的目的.对人工合成和实际的32D地震资料的处理结果表明,此法具有好的应用前景.2　理论与算法211　32D地震信号的奇异值分解[8]将32D地震信号d(x,y,t)的每一道进行傅里叶变换,得到对应的f-x-y域32D数据体D(x,y,ω).32D地震信号的奇异值分解是通过对D(x,y,ω)的每一个频率切片进行奇异值分解来实现的.假设32D数据体的纵测线和横测线数目分别为N和M,则对应某一个频率ω可以得到一个22D切片,我们用一个N×M的矩阵Z(ω)表示,其元素为zij(ω) =D(i,j,ω)(i为纵测线号,j为横测线号).根据奇异值分解理论,Z(ω)可写成Z(ω)=UΣV T=∑rk=1σku k v T k,(1)式中,r为Z(ω)的秩,U,V分别是Z(ω)[Z(ω)]T的特征向量ui(i=1,…,N),vj(j=1,…,M)组成的矩阵,又被称为Z(ω)的奇异向量矩阵,Σ是Z(ω) [Z(ω)]T的特征值按递减顺序组成的对角矩阵,又被称为Z(ω)的奇异值矩阵,这些矩阵分别表示成Σ=diag(σ1,σ2,…,σr),(2)U=[u1,u2,…,u N],(3)V=[v1,v2,…,v M]T.(4) 利用奇异值分解来提高地震信号的信噪比是通过去掉一些对应小奇异值的噪声子空间,即选取一个滤波参数l(l<r),重构一个矩阵 Z(ω),来逼近原始矩阵Z(ω),Z(ω)=∑lk=1σku k v T k.(5) 处理完所有频率切片后,可以得到一个提高信噪比后的32D数据体 D(x,y,ω),通过反傅里叶变换,就可以得到滤波后的32D地震数据体 d(i,j,ω).2.2　自适应提高分辨率利用每一个奇异值及其相应的奇异向量,我们可以重构出一个子空间(信号或噪声).也就是说,频率切片Z(ω)可以分解为r个子空间:798　4期陆文凯等:基于信号子空间分解的三维地震资料高分辨率处理方法Z (ω)=∑rk =1z k (ω),(6)z k (ω)=σk u k v Tk ,k =1,2,…,r .(7) 将所有频率切片都分解为r 个子空间后,则对应相同下标的信号子空间可以组合到一起,形成和32D 数据体D (i ,j ,ω)对应的子数据体c k (i ,j ,ω).也就是说,32D 数据体D (i ,j ,ω)可以分解为r 个子数据体:D (i ,j ,ω)=∑rk =1ck(i ,j ,ω),(8)其中,子数据体c k (i ,j ,ω)在整个数据体D (i ,j ,ω)所占的能量比例为λk =σk∑rk =1σk.(9) 由于σ1≥σ2≥…≥σr ,所以λ1≥λ2≥…≥λr .考虑到信号都集中在对应大奇异值的子空间中,而噪声则一般都集中在对应小奇异值的子空间中[4,5,7,8],所以对应大奇异值的子空间,其信噪比高,有效频宽相对原始地震数据来说得到拓宽[7].对第k 个子空间,可以根据其对应的能量比例确定其为噪声子空间还是信号子空间[4,5].如果为噪声子空间,则不参加信号重构;如果是信号子空间,则根据λk 对其进行自适应地提高分辨率处理,然后参加进一步的信号重构.谱白化是一种展宽频谱的基本方法,是一种“纯振幅”的滤波过程.谱白化处理既可以在频率域中实现,也可以在时间域中实现.考虑到32D 地震信号经过奇异值分解后,再重构得到单个子数据体c k (i ,j ,ω)是频率域的数据,我们选择了频率域的谱白化技术.在频率域中将每一道地震信号的振幅谱的各个极大值尖峰连成没有极小值的外包线,将其进行平滑滤波后,得到地震信号的振幅谱包络e k (i ,j ,ω).则实现谱白化的滤波器设计为f k (i ,j ,ω)=v (i ,j )e k (i ,j ,ω)+ε(λk )v (i ,j ),(10)式中,v (i ,j )为e k (i ,j ,ω)的最大值,ε(λk )为“白噪成分”,是用来折衷信噪比和分辨率的.一般来说,“白噪成分”越大,分辨率降低,信噪比增强,而“白噪成分”越小,则分辨率提高,信噪比降低.在我们提出的方法中,“白噪成分”ε(λk )是随λk 而自适应地设定,总的来说,ε(λk )应该随λk 的增大而变小,即对信噪比高的子空间,我们展宽其频谱宽一些,而对信噪比低的子空间,我们展宽其频谱小一些,甚至不展宽或者缩窄其频谱.我们的方法采用了如下的自适应调节策略:ε(λk )=ε0+ε1(1-λk ),(11)式中,ε0,ε1为常数.假设前l 个子空间为信号子空间,经过高分辨率处理后,重构得到的地震信号为D (i ,j ,ω)=∑l k =1ck(i ,j ,ω)f k (i ,j ,ω).(12) 经过反傅里叶变换,得到既提高分辨率,又提高信噪比的32D 地震数据.本方法的处理流程如图1所示.图1　本方法流程图Fig.1　Flow 2chart of the proposed method3　资料处理为了检验本文方法的有效性,我们将其用来处理人工合成的32D 信号和实际地震数据.人工合成的三维数据共16条测线,每条测线40道,每道地震信号800ms ,时间采样率为4ms.图2a 显示了第一条测线剖面,剖面是由一个尖灭模型和一个Ricker 子波通过褶积形成的,剖面中共有两条同相轴,一条为水平同相轴,一条为倾斜同相轴.第一条测线的尖灭点在第11道的400ms 处,而第i 条测线的尖灭点位于第11+i 道的400ms 处.整个合成数据体加上-12dB 的高斯白噪声,图2b 显示了加噪后的第一条测线的剖面.对加噪后的32D 数据体进行奇异值分解后,用第一个奇异值对应的子空间重构的数据体的第一条测线剖面显示在图2c 中,用第二和第三两个奇异值对应的子空间重构的数据体的第一条测线剖面显示在图2d 中,用剩余所有奇异值对应的子空898地球物理学报(Chinese J.G eophys.)48卷　图2　合成地震资料的处理(a)无噪地震剖面;(b)加噪地震剖面;(c)第一个奇异值对应子空间的重构结果;(d)第二和第三个奇异值对应子空间的重构结果;(e)其余奇异值对应子空间的重构结果;(f)前三个奇异值对应子空间的重构结果(即滤波结果);(g)常规谱白化技术处理结果;(h)本文方法处理结果.Fig2Processing of synthetic seismic data(a)N oisy2free seism ic profile;(b)N oisy seism ic profile;(c)Profile reconstructed by the first singular value;(d)Profile reconstructedby the2th and3th singular values;(e)Profile reconstructed by the other singular values;(f)Filtered result using the first three singular values;(g)Result obtained by the spectra whitening technique;(h)Result obtained by our method.998　4期陆文凯等:基于信号子空间分解的三维地震资料高分辨率处理方法图3　实际地震资料处理(a )原始测线剖面;(b )本方法处理结果;(c )谱白化处理结果.Fig.3　Application of the real seismic data(a )Original profile ;(b )Result obtained by our method ;(c )Result obtained by the spectra whitening technique.间重构的数据体的第一条测线剖面显示在图2e 中.可以看出,利用第一个奇异值对应的子空间进行信号重构就可以分离出水平同相轴,而利用第二和第三个奇异值对应的子空间进行信号重构就可以分离出倾斜同相轴,剩余所有奇异值对应的子空间重构的数据体基本上为随机噪声.显然,利用前三个奇异值对应的子空间重构,就可以得到噪声压制后的数据体,其第一条测线剖面显示在图2f 中.常规谱白化的结果显示在图2g 中,所处理的频率范围0～45H z ,谱白化“白噪成分”为振幅谱包络的最大值的01001.利用我们的方法进行高分辨率处理后得到的数据体的第一条测线剖面显示在图2h 中,在这次试验中,所处理的频率范围为0～45H z ,l =3,且ε0=01001,ε1=011.可以看出,我们的方法既能有效地提高地震资料的分辨率,又能提高其信噪比.为了进一步对处理结果进行量化评价,我们将反射系数的0～45H z 低通滤波信号作为期望结果,计算它和两种处理结果的归一化互相关系数.理论上讲,某种方法的处理结果和期望结果的归一化互相关系数越接近1,则说明该方法的性能好.在上述试验中,我们方法得到的结果和期望结果的归一化互相关系数为016158,而谱白化的结果和期望结果的归一化互相关系数为015257.我们处理了一个实际的32D 数据,共有157条测线,每条测线1037道,每道4000点,时间采样率为1ms.在实际数据处理时,我们的方法采用了分块的策略,且块与块之间存在重叠以消除边界效应,所以可以适应空变和时变的情况.在这次试验中,分块尺寸为16×40×128,即16条测线,每条测线40道,每道128点,且3个方向都是50%重叠.图3a 显示了第120条测线的一部分,图3b 和3c 分别给出了我们方法和常规谱白化的结果.在本次试验中,我们方法所处理的频率范围为0～145H z ,l =3,且ε0=01001,ε1=0101.而谱白化技术所处理的频率范围0～145H z ,谱白化“白噪成分”为振幅谱包络的最大值的01001.图4给出了图3地震资料相应的振幅009地球物理学报(Chinese J.G eophys.)48卷　谱.从图3右上角部分的尖灭可以看出,本文提出的方法相对常规的谱白化技术可以更有效地提高地震资料的分辨率,从而给出更多的地质沉积现象.图4所给出的振幅谱也可以得出同样的结论.图4　和图3中数据对应的振幅谱Fig.4　Am plitude spectrum of the data shown in Fig.34　结　论32D 地震资料经过奇异值分解以后,利用不同子空间进行重构,不仅可以分离信号和噪声,而且可以根据地震同相轴的倾角差异,分离重叠在一起的复合地震信号.本文给出的方法,根据每个信号子空间的能量比重,自适应地设计频谱展宽滤波器,从而达到既提高地震资料的分辨率,又提高地震资料的信噪比的目的.初步的人工合成和实际的32D 地震资料处理给出了令人满意的结果.参考文献(References )[1]李庆忠.走向精确勘探的道路.北京:石油工业出版社,1994Li Q Z.The W ay to Obtain a Better Res olution in Seism ic Prospecting (in Chinese ).Beijing :Petroleum Industry Press ,1994[2]李鲲鹏,刘业新,李衍达等.小波变换的过零点特性与地震勘探信号的信噪比和分辨率.地球物理学报,1997,40(4):561～569Li Q P ,Liu Y X ,Li Y D ,et al.Z ero crossing property of wavelet trans form and im provement of both seism ic signal noise ratio and res olution.Chinese J .G eophys .(in Chinese ),1997,40(4):561～569[3]章　珂,李衍达,刘贵忠等.多分辨率地震信号反褶积.地球物理学报,1999,42(4):529～535Zhang K,Li Y D ,Liu G Z ,et al.Multires olution seism ic signal deconv olution.Chinese J .G eophys .(in Chinese ),1999,42(4):529～535[4]陆文凯,牟永光.一种改进的S VD 滤波器.石油地球物理勘探,1996,31(5):736～741Lu W K,Mu Y G.An im proved S VD filter.Oil G eophysicalProspecting (in Chinese ),1996,31(5):736～741[5]Lu W.Adaptive noise attenuation of seism ic image using singular value decom position and texture direction detection.ICIP2002,2002[6]陆文凯,李衍达.S VD 任意道间距内插.石油地球物理勘探,1997,32(4):582～588Lu W K,Li Y D.Any -interval trace interpolation using S VD method.Oil G eophysical Prospecting (in Chinese ),1997,32(4):582～588[7]陆文凯,李衍达.S VD 分解提高地震资料的分辨率.石油地球物理勘探,1998,33(增刊):145～149Lu W K,Li YD.Im provement of seism ic signal S NR and res olution using S VD method.Oil G eophysical Prospecting (in Chinese ),1998,33(Suppl.):145～149[8]S tewart R.T rickett.F -x -y eigenimage noise suppression.G eophysics ,2003,68(2):751～759109　4期陆文凯等:基于信号子空间分解的三维地震资料高分辨率处理方法。