小波變換在面波插值中的應(yīng)用

周聰，劉江平，曾祥芝，范承余，曹進(jìn)

(1.中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，湖北 武漢，430074；2.安徽省高速公路控股集團(tuán)有限公司，安徽 合肥，230051)

在二維地震勘探中，因?yàn)槭芤巴馐┕l件(如障礙物、禁區(qū)等)、成本、效益等多方面因素的限制，野外采集道間距往往過大，容易產(chǎn)生空間假頻，從而影響地震資料的多道處理。而利用地震道插值能加密空間采樣率，防止頻散出現(xiàn)和提高信噪比[1]。目前，地震道內(nèi)插的方法有很多。Spitz[2]提出在f-x域中進(jìn)行地震道內(nèi)插，該方法計(jì)算量較大且精度較低。此后，國(guó)九英等[3-5]在此基礎(chǔ)上做了更深入的研究。Wang等[6-7]將小波變換理論引入地震道插值中，并指出小波插值方法是一種能夠消除假頻并擴(kuò)展帶寬的地震道插值方法。Trad等[8-9]采用預(yù)優(yōu)共軛梯度法對(duì)大矩陣求逆，用雙曲Radon變換對(duì)CMP道集數(shù)據(jù)進(jìn)行插值處理，取得較好的插值效果。崔興福等[10]也利用小波變換在時(shí)間和頻率域良好的局部化性質(zhì)，通過一維小波變換重構(gòu)公式實(shí)現(xiàn)地震道內(nèi)插。王維紅等[11]基于帶限正反最小平方拋物Radon變換的Levinson遞推算法，給出利用拋物Radon變換進(jìn)行地震道重建外推的基本原理以及疊前地震數(shù)據(jù)規(guī)則化地處理流程。Zhou等[12]提出在Wavelet-Radon域校正假頻以及地震數(shù)據(jù)的插值，該方法可以非?？煽康鼗謴?fù)高頻信號(hào)，在合成數(shù)據(jù)和實(shí)際數(shù)據(jù)的處理中都取得很好的效果。Herrmann等[13-14]利用Curvelet變換的稀疏特性，成功將其應(yīng)用于地震數(shù)據(jù)的重建。劉國(guó)昌等[15]在Curvelet變換的基礎(chǔ)上，提出基于POCS算法的Curvelet變換震數(shù)據(jù)插值方法。在以往的反射地震勘探中，面波往往被認(rèn)為是一種具有明顯高振幅和頻率特性的規(guī)則干擾波而被壓制[16]。而在多波多分量地震勘探中，由于橫波靜校正的需要，許多學(xué)者都利用瑞雷面波的頻散信息反演淺層橫波速度結(jié)構(gòu)，進(jìn)而求取靜校正量[17-18]。面波的反演主要是基于頻散曲線的反演，而對(duì)面波數(shù)據(jù)進(jìn)行有效的插值可以提高頻率速度譜的信噪比，有利于頻散曲線的精確提取。在大時(shí)空采樣地震記錄中，由于面波的速度很低，其記錄中的高頻成分容易出現(xiàn)超周期現(xiàn)象，降低其頻率速度譜的信噪比，而且容易出現(xiàn)假的高階能量。而對(duì)面波插值，能夠在一定程度上消除假頻，提高其頻率速度譜的信噪比。為此，本文作者利用面波的線性特征，首先對(duì)原始面波數(shù)據(jù)進(jìn)行線性動(dòng)校正(LNMO)處理，然后利用二維小波變換對(duì)其插值，最后反線性動(dòng)校(RLNMO)恢復(fù)信號(hào)。理論和實(shí)際數(shù)據(jù)插值前后頻率速度譜的變化表明：本文方法能有效提高頻率速度譜的信噪比和頻散曲線的提取精度。

1 小波插值的基本原理

以一維小波變換來說明小波插值的基本原理。根據(jù)小波變換的基本原理，信號(hào)在分辨率為2-m-1空間的近似部分可以由分辨率為2-m空間的近似部分和細(xì)節(jié)部分通過小波重構(gòu)來得到[6]。如圖1所示，將原始信號(hào)c0看成是分辨率為20的空間的近似部分，如果該空間的細(xì)節(jié)部分d0可以忽略或估計(jì)，那么分辨率為2-1空間的近似信號(hào)c-1就可以由原始信號(hào)得到。

圖1 小波重構(gòu)示意圖Fig.1 Sketch of wavelet reconstruction

反變換可表示為：

在式(4)中，是已知的，而未知?？梢酝ㄟ^以下2種方法來確定。

(1)方法1。由式(4)可知，若分辨率為20空間的細(xì)節(jié)部分可以忽略，即令，則可得：

在很多情況下，是可以忽略的。這時(shí)，信號(hào)通過插值來增加信息量是一個(gè)有效的處理過程。

(2)方法2。由于忽略相當(dāng)于將信號(hào)進(jìn)行簡(jiǎn)單的去噪處理，剔除該信號(hào)中的高頻部分，因此，對(duì)于需要保留高頻信息的信號(hào)來說，方法1可能達(dá)不到該要求。在圖像插值里面，一般是利用原圖像來確定高頻信息[19-20]。具體做法是分解原圖像得到高頻信息，由于不同分辨率下各個(gè)對(duì)應(yīng)子圖間存在相似性，可以通過一般插值(如3次樣條)得到和原始圖像同一分辨率下的高頻信號(hào)。王真理等[7]指出為得到的信息，可對(duì)原始數(shù)據(jù)進(jìn)行空間方向高通濾波，或采用頻率調(diào)制方法，提高其頻率，并保持波形與原剖面的相似性。

由于利用方法2得到的會(huì)包含一些不真實(shí)的信息以及其實(shí)現(xiàn)的復(fù)雜性，因此，本文嘗試采用方法1。在二維小波重構(gòu)中，需要將3個(gè)細(xì)節(jié)部分置0。通過面波的頻率速度譜來分析插值的效果，特別是其高頻部分頻散能量的變化。本文選取的母小波是db4小波。

2 模擬數(shù)據(jù)試驗(yàn)結(jié)果分析

下面對(duì)數(shù)值模擬的體波和面波記錄做小波插值，分別從體波插值結(jié)果的誤差以及面波插值前后頻率速度譜的變化來說明本文插值方法的有效性。

2.1 體波插值

如圖2所示地震記錄為二維合成記錄，道間距2.5 m，道數(shù)為480，時(shí)間采樣率為0.5 ms，采樣長(zhǎng)度800 ms，雷克子波主頻35 Hz，同相軸的傾角分別為1 ms/道和0.5 ms/道，對(duì)該記錄做3次小波插值試驗(yàn)。將原始數(shù)據(jù)分別均勻抽道成5，10和20 m，然后插值成2.5 m，并與原始合成記錄進(jìn)行比較。由于小波插值是同時(shí)對(duì)時(shí)間和空間進(jìn)行插值，而在地震插值中我們主要關(guān)心道間插值，所以，二維插值實(shí)驗(yàn)中舍棄時(shí)間方向上插值的數(shù)據(jù)，只考慮道間插值。

圖2 二維合成地震記錄Fig.2 2D synthetic seismic data

小波插值實(shí)驗(yàn)結(jié)果見圖3～5。從圖3～5可見：隨著道間距和插值倍數(shù)的增加，小波插值的誤差也越來越大。圖3(b)顯示的插值結(jié)果很好，僅僅是在同相軸相交處稍差，而圖5(b)中傾角大的同相軸的插值結(jié)果要遠(yuǎn)差于傾角小的同相軸。所以，對(duì)于不考慮細(xì)節(jié)部分的小波插值方法來說，雖然其能夠在一定程度上消除假頻，但較嚴(yán)重的假頻反過來又會(huì)影響插值的效果。

2.2 面波插值

如圖6(a)所示為根據(jù)3層介質(zhì)模擬出的面波記錄，其道間距為2 m，時(shí)間采樣率為0.5 ms，雷克子波主頻為20 Hz。將其均勻抽道得到道間距4 m的試驗(yàn)信號(hào)，如圖6(b)所示。利用小波變換直接對(duì)該記錄插值，得到結(jié)果如圖7(a)所示。由圖7(a)可見：插值效果不理想，插值出來的中高頻部分幾乎是拷貝其左邊道的同相軸信息。若同相軸斜率較大，則其高頻部分容易出現(xiàn)假頻信息，會(huì)影響插值效果?？梢酝ㄟ^2種數(shù)據(jù)處理的方法來處理：一是通過坐標(biāo)旋轉(zhuǎn)，將數(shù)據(jù)轉(zhuǎn)動(dòng)某一角度使得大傾角同相軸的斜率變小，但此時(shí)在新坐標(biāo)系中的縱橫向采樣率發(fā)生變化，難以處理且可行性有待驗(yàn)證；另一種是對(duì)于線性信號(hào)來說(如面波)，可以先對(duì)原始數(shù)據(jù)做類似線性動(dòng)校正的處理，然后進(jìn)行小波插值，最后反線性動(dòng)校正恢復(fù)信號(hào)即可。后者即相當(dāng)于單坐標(biāo)軸的旋轉(zhuǎn)。

圖3 小波插值試驗(yàn)(1倍插值)Fig.3 Wavelet interpolation experiment (1× interpolation)

圖4 小波插值試驗(yàn)(2倍插值)Fig.4 Wavelet interpolation experiment (2× interpolation)

圖5 小波插值試驗(yàn)(3倍插值)Fig.5 Wavelet interpolation experiment (3× interpolation)

圖6 理論面波記錄Fig.6 Synthetic rayleigh waves

圖7 小波變換對(duì)瑞雷面波插值Fig.7 Rayleigh wave interpolation using wavelet transform

下面對(duì)圖6(b)所示的面波記錄進(jìn)行一定程度的線性動(dòng)校正處理。由于要往上“拉伸”面波，因此，需要在面波記錄的零時(shí)刻前端補(bǔ)零，以防止面波被拉出邊界。經(jīng)線性動(dòng)校處理后的結(jié)果如圖7(b)所示，再對(duì)該記錄做小波插值，然后做反線性動(dòng)校處理恢復(fù)信號(hào)，如圖7(c)所示。將圖6(a)，7(a)和7(c)中第36道抽出來進(jìn)行對(duì)比，結(jié)果見圖7(d)。由圖7(d)可見：經(jīng)過線性動(dòng)校處理的插值效果要優(yōu)于直接利用小波插值，特別是在面波記錄的后半部，如0.4～0.6 s之間。從整個(gè)記錄剖面來看，除了記錄2端的幾道，結(jié)合線性動(dòng)校處理的小波插值能夠很好地恢復(fù)面波信號(hào)。對(duì)小波插值方法來說，記錄2端道的插值效果較差是其缺陷，本文作者認(rèn)為在小波插值過程中，默認(rèn)記錄以外的部分為零值，所以，插值出來的邊界道能量要弱。該處理方法通過單坐標(biāo)軸的旋轉(zhuǎn)，可減小線性同相軸的斜率，避開直接對(duì)大時(shí)空采樣地震記錄進(jìn)行小波插值，對(duì)于同相軸為近似線性的面波記錄的道內(nèi)插有較好的效果。該處理過程的重點(diǎn)部分是線性動(dòng)校正量的估算。以面波為例，面波在地震記錄上呈“掃帚”狀，可以大致畫出“掃帚”的上、下界限，即可以估算出面波的最大速度和最小速度；然后可以取平均速度作為參考速度；最后，利用道間距計(jì)算出校正量即可。

2.3 面波的頻率速度譜的變化

在轉(zhuǎn)換波靜校正中，地表橫波速度模型的獲取尤其關(guān)鍵，而面波速度非常趨近于橫波，因此，通過面波反演來獲取橫波速度是可行的。面波的反演主要是基于頻散曲線的反演，而對(duì)面波數(shù)據(jù)進(jìn)行插值可以提高頻率速度譜的信噪比，有利于頻散曲線的提取。

下面本文從插值前后面波頻率速度譜的變化來分析插值的結(jié)果。提取頻散曲線可以通過2個(gè)線性變換實(shí)現(xiàn)[21-22]：首先對(duì)原始數(shù)據(jù)做傾斜疊加(即線性Radon變換)，將數(shù)據(jù)變?yōu)槁?截距域(p-τ)；其次，沿截距τ方向做一維傅里葉變換得到慢度-頻率域(p-f)；最后，將慢度轉(zhuǎn)換成速度即可獲得速度-頻率域(v-f)。

如圖8(a)～(c)所示分別是圖6(a)，6(b)和7(b)所示面波通過線性Radon變換得到的頻率速度譜。對(duì)比圖8(a)和8(b)可見：均勻抽道后，在能譜的高頻高速區(qū)域出現(xiàn)一個(gè)較強(qiáng)的假高階能量。這可能是由于抽道后，面波的高頻部分產(chǎn)生超周期現(xiàn)象，其高頻信息的同相軸斜率由小變大，因此，其高頻信息的速度增大，出現(xiàn)了如圖8(b)中箭頭A所示的假高階能量。在箭頭B所示的低速區(qū)域，也產(chǎn)生干擾信息。對(duì)面波插值后，如圖8(c)所示，箭頭A處的假高階能量團(tuán)明顯減弱，而箭頭B處的噪音干擾也得到改善。從插值前后頻率速度譜的變化情況說明本文的插值方法可以對(duì)面波的高頻信息進(jìn)行有效插值，并能夠提高其頻率速度譜的信噪比。

3 實(shí)例分析

下面將本文的小波插值方法應(yīng)用于實(shí)際數(shù)據(jù)。同樣，根據(jù)記錄插值前后面波頻率速度譜的變化來驗(yàn)證插值方法的有效性。

3.1 野外地震記錄插值

圖8 理論面波數(shù)據(jù)的頻率速度譜Fig.8 Dispersive images in f-v domain of synthetic data

如圖9(a)所示為野外采集的原始地震記錄，時(shí)間采樣率為2 ms，記錄長(zhǎng)度為7 s，道間距為20 m，記錄道數(shù)36道(起始道為零道)。圖9(a)中上部黑色部分為體波。圖9(b)所示為將原始地震記錄均勻抽道后得到的地震記錄，道間距為40 m。利用本文的插值方法對(duì)該記錄進(jìn)行插值，得到圖9(c)所示的插值結(jié)果。對(duì)比圖9(a)和9(c)可以看到：面波插值效果較好，得到的面波記錄與原始記錄中的面波基本一致，沒有帶入背景噪聲。將圖9(a)和9(c)中第16道抽出來對(duì)比，得到圖9(d)?？梢钥吹矫娌ㄖ黧w部分的波形匹配較好，只是在面波的尾部如時(shí)間軸4.3 s附近重構(gòu)效果稍差。從圖9(a)，(c)和(d)均可以發(fā)現(xiàn)體波部分發(fā)生較大變化。這是由于在插值前進(jìn)行線性動(dòng)校處理，體波同相軸多是非線性的，在“拉平”面波的同時(shí)，將體波的同相軸破壞。因此，本文的插值方法無法同時(shí)對(duì)體波和面波進(jìn)行有效插值。由于體波速度較高，其假頻現(xiàn)象要比面波記錄中的少。當(dāng)需要對(duì)體波進(jìn)行有效插值時(shí)，可以直接利用小波變換對(duì)體波插值。

圖9 實(shí)際地震數(shù)據(jù)Fig.9 Real seismic data

3.2 插值前后頻散曲線的變化

如圖10(a)～(c)所示分別是圖9(a)～(c)通過線性Radon變換得到的頻率-速度譜。對(duì)比圖10(a)和(b)可以看到：抽道前后其頻率速度譜發(fā)生明顯變化，抽道后地震記錄的頻率速度譜的信噪比降低很多，在高頻高速的區(qū)域出現(xiàn)了3個(gè)假能量團(tuán)，而且2個(gè)能量特別強(qiáng)，這有可能是面波高頻率信息超周期的現(xiàn)象造成的；在高頻低速的區(qū)域也出現(xiàn)了較多的噪音干擾；其二階的頻散能量發(fā)生明顯彎曲，尾部的能量團(tuán)也產(chǎn)生斷裂，對(duì)面波插值后，高頻高速區(qū)域的3個(gè)假能量團(tuán)均已消失，說明面波中高頻部分的假頻信息得到改善；高頻低速區(qū)域的噪音干擾也已消除，二階的頻散能量團(tuán)較插值前更圓滑，而且其尾部斷裂的能量團(tuán)也得到改善。對(duì)比圖10(a)和(c)可見：除了三階頻散能量團(tuán)沒有較好地恢復(fù)外，其他的與原始記錄的頻率速度譜基本相同。這可能是由于抽道后該部分信息已經(jīng)缺失，而插值無法恢復(fù)缺失的信息，只能突出已有的有效信息，壓制干擾信息。因此，實(shí)際資料的頻率速度譜的變化也說明本文插值方法能夠較好地應(yīng)用于實(shí)際，提高頻散-速度譜的信噪比，有利于精確提取面波頻散曲線。

圖10 實(shí)際地震資料的頻率速度譜Fig.10 Dispersive images in f-v domain of real seismic data

4 結(jié)論

(1)基于小波變換的基本原理，采用db4母小波，實(shí)現(xiàn)了對(duì)瑞雷面波的時(shí)空插值。

(2)根據(jù)瑞雷面波低速和線性特點(diǎn)，將線性動(dòng)校正引入小波插值中，有效實(shí)現(xiàn)了對(duì)瑞雷面波的插值，改善了瑞雷面波的插值效果。

(3)在大時(shí)空采樣情況下，對(duì)瑞雷面波進(jìn)行插值，可有效提高頻率速度譜的信噪比，有利于瑞雷面波頻散曲線的精確提取。

[1]羅丹.f-k域地震道插值方法研究[D].成都: 成都理工大學(xué)地球物理學(xué)院, 2009: 1-3.LUO Dan.The research of methods to seismic trace interpolation in f-k domain[D].Chengdu: Chengdu University of Technology.College of Geophysics, 2009: 1-3.

[2]Spitz S.Seismic trace interpolation in the f-x domain[J].Geophysics, 1991, 56(6): 785-794.

[3]國(guó)九英, 周興元, 俞壽朋.f-x域等道距道內(nèi)插[J].石油地球物理勘探, 1996, 31(1): 28-34.GUO Jiu-ying, ZHOU Xing-yuan, YU Shou-peng.Iso-interval trace interpolation in f-x domain[J].Oil Geophysical Prospecting,1996, 31(1): 28-34.

[4]Porsani M J.Seismic trace interpolation using half-step prediction filters[J].Geophysics, 1999, 64(5): 1461-1467.

[5]WANG Yang-hua.Seismic trace interpolation in the f-x-y domain[J].Geophysics, 2002, 67(4): 1232-1239.

[6]WANG Zhen-li, LI Yan-da.Trace interpolation using wavelet transform[C]//64th SEG Annual International Meeting, Los Angles, USA: SEG, 1994, 13: 729-730.

[7]王真理, 李衍達(dá).基于小波的信號(hào)部分重構(gòu)與插值[J].清華大學(xué)學(xué)報(bào): 自然科學(xué)版, 1995, 35(5): 67-72.WANG Zhen-li, LI Yan-da.Signal partial reconstruction and interpolation based on wavelet transform[J].Journal of Tsinghua University: Science and Technology, 1995, 35(5): 67-72.

[8]Trad D O, Ulrych T J, Sacchi M D.Accurate interpolation with high-resolution time-variant Radon transforms[J].Geophysics,2002, 67(2): 644-656.

[9]Trad D O.Interpolation and multiple attenuation with migration operators[J].Geophysics, 2003, 68(6): 2043-2054.

[10]崔興福, 劉東奇, 張關(guān)泉.小波變換實(shí)現(xiàn)地震道插值[J].石油地球物理勘探, 2003, 38(增刊): 93-97.CUI Xing-fu, LIU Dong-qi, ZHANG Guan-quan.Seismic trace interpolation using wavelet transform[J].Oil Geophysical Prospecting, 2003, 38(Suppl): 93-97.

[11]王維紅, 劉洪.拋物Radon變換法近偏移距波場(chǎng)外推[J].地球物理學(xué)進(jìn)展, 2005, 20(2): 289-293.WANG Wei-hong, LIU Hong.Near offset wavefield extrapolation based on parabolic radon transform[J].Progress in Geophysics, 2005, 20(2): 289-293.

[12]ZHOU Yu, Ferguson J, McMechan G, et al.Wavelet-Radon domain dealiasing and interpolation of seismicdata[J].Geophysics, 2007, 72(2): 41-49.

[13]Herrmann F J, Hennenfent G.Non-parametric seismic data recovery with curvelet frames[J].Geophys J Int, 2008, 148(1):233-248.

[14]Naghizadeh M, Sacchi D.Beyond alias hierarchical scale curvelet interpolation of regularly and irregularly sampled seismic data[J].Geophysics, 2010, 75(6): WB189-WB202.

[15]劉國(guó)昌, 陳小宏, 郭志峰, 等.基于Curvelet變換的缺失地震數(shù)據(jù)插值方法[J].石油地球物理勘探, 2011, 46(2): 237-246.LIU Guo-chang, CHEN Xiao-hong, GUO Zhi-feng, et al.Missing seismic data rebuilding by interpolation based on Curvelet transform[J].Oil Geophysical Prospecting, 2011, 46(2):237-246.

[16]張恒磊, 劉天佑, 李紅巧.Curvelet域面波衰減方法研究[J].中南大學(xué)學(xué)報(bào): 自然科學(xué)版, 2011, 42(8): 2372-2378.ZHANG Heng-lei, LIU Tian-you, LI Hong-qiao.Attenuation of surface wave in Curvelet domain[J].Journal of Central South University: Science and Technology, 2011, 42(8): 2372-2378.

[17]孟小紅, 郭良輝.利用地震瑞利波速度反演求取P-SV波橫波靜校正量[J].石油地球物理勘探, 2008, 42(4): 448-453.MENG Xiao-hong, GUO Liang-hui.Using velocity inversion of seismic Rayleigh wave to compute S-wave statics of P-SV wave[J].Oil Geophysical Prospecting, 2007, 42(4): 448-453.

[18]Bansal R, Ross W, Lee S, et al.A novel approach to estimating near-surface S-wave velocity and converted-wave receiver statics[C]//79th SEG Annual International Meeting, Houston,USA: SEG, 2009, 28: 1192-1196.

[19]盧玨.基于小波的圖像插值研究[J].武漢理工大學(xué)學(xué)報(bào), 2003,25(1): 81-83.LU Jue.Research of image interpolation based on wavelet[J].Journal of Wuhan University of Technology, 2003, 25(1): 81-83.

[20]陶洪久, 柳健, 田金文.基于小波變換和插值的超分辨率圖像處理算法[J].武漢理工大學(xué)學(xué)報(bào), 2002, 24(8): 63-66.TAO Hong-jiu, LIU Jian, TIAN Jin-wen.Surper resolution image processing algorithm based on wavelet transform and bilinear interpolation[J].Journal of Wuhan University of Technology, 2002, 24(8): 63-66.

[21]George A M, Mathew J Y.Analysis of dispersive waves by wave field transformation[J].Geophysics, 1981, 46(6): 869-874.

[22]LUO Yin-he, XIA Jiang-hai, Richard D M, et al.Rayleigh-wave dispersive energy imaging using a high-resolution linear radon transform[J].Pure and Applied Geophysics, 2008, 165(5):903-922.