小波變換在面波插值中的應用

周聰，劉江平，曾祥芝，范承余，曹進

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

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

1 小波插值的基本原理

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

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

反變換可表示為：

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

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

在很多情況下，是可以忽略的。這時，信號通過插值來增加信息量是一個有效的處理過程。

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

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

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

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

2.1 體波插值

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

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

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

2.2 面波插值

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

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

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

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

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

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

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

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

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

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

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

3 實例分析

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

3.1 野外地震記錄插值

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

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

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

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

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

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

4 結(jié)論

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

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

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

[1]羅丹.f-k域地震道插值方法研究[D].成都: 成都理工大學地球物理學院, 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]國九英, 周興元, 俞壽朋.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]王真理, 李衍達.基于小波的信號部分重構(gòu)與插值[J].清華大學學報: 自然科學版, 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)泉.小波變換實現(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變換法近偏移距波場外推[J].地球物理學進展, 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]劉國昌, 陳小宏, 郭志峰, 等.基于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].中南大學學報: 自然科學版, 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].武漢理工大學學報, 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].武漢理工大學學報, 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.