林浩然，邢磊,2*，劉懷山,2，李倩倩，張洪茂

1 中國(guó)海洋大學(xué)海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，青島 266100 2 青島海洋科學(xué)與技術(shù)試點(diǎn)國(guó)家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室，青島 266071

0 引言

面波是地震資料中最為常見(jiàn)的規(guī)則干擾波，其是由體波到達(dá)地表或介質(zhì)分界面時(shí)，在一定條件下相互干涉并疊加產(chǎn)生的次聲波.面波具有低頻、低速、強(qiáng)振幅、能量衰減緩慢等特點(diǎn)，在炮集記錄中呈近似線性分布，面波的存在嚴(yán)重影響中深層有效信號(hào)，顯著降低了地震資料的信噪比(王偉奇等，2022).

在地震資料處理過(guò)程中，研究人員主要依據(jù)面波低頻、低視速度的特點(diǎn)，提出了許多壓制面波的方法，如高通濾波、FK濾波、Radon域?yàn)V波、KL變換域?yàn)V波、維納濾波、SVD濾波、徑向道分解濾波、二維時(shí)間導(dǎo)數(shù)濾波等.這些方法雖具有一定的壓制效果，但是往往會(huì)對(duì)有效信號(hào)造成一定程度的損傷.為了獲得更好的面波壓制效果，提高地震資料信噪比，Deighan和Watts (1997)提出將一維小波變換應(yīng)用于面波壓制工作中.小波變換具有良好時(shí)頻分辨能力，以及信號(hào)奇異性捕獲能力(陳永芮等，2021).近年來(lái)，在小波變換的基礎(chǔ)上，具有更好逼近高維奇異信號(hào)能力的多尺度幾何分析方法，如脊波(Ridgelet)變換(Candès and Donoho,1999;包乾宗等,2007)、曲波(Curvelet)變換(Candès and Donoho,2004)等被引入面波壓制方法研究領(lǐng)域中來(lái).其優(yōu)越之處在于能更好的利用面波與有效波的視速度差異，捕捉地震信號(hào)的奇異特性，同時(shí)能更好的刻畫圖像邊緣信息.然而這些時(shí)頻分析方法都受到海森堡—蓋博測(cè)不準(zhǔn)原理的影響，導(dǎo)致其時(shí)頻分辨能力受到一定限制.為克服這一問(wèn)題并提高它們表示非平穩(wěn)信號(hào)的能力，基于銳化時(shí)頻表示技術(shù)(Kodera et al.,1976)和時(shí)頻重分配思想(Auger and Flandrin,1995)，Daubechies和Maes(1996)、Daubechies等(2011)提出了同步擠壓算法，并證明其可以應(yīng)用于地震數(shù)據(jù)處理領(lǐng)域.劉晗等(2016)提出將同步擠壓小波變換應(yīng)用于面波壓制中.時(shí)頻重分配算法現(xiàn)階段已被推廣至二維(Yang et al.,2015;Lu et al.,2016)和三維空間(Zhu et al.,2021)中，且有多種變體(Clause et al.,2015;Lu and Yang,2018)和改進(jìn)(陳文超等,2009;Oberlin and Meignen,2017;Daubechies et al.,2016).

隨著處理要求的提高和處理技術(shù)的進(jìn)步，近幾年來(lái)提出了一些采用多種手段聯(lián)合壓制面波的方法.包乾宗等(2011)提出了聯(lián)合小波變換和脊波變換的面波壓制方法.馬見(jiàn)青和李慶春(2011)提出了S變換和TT變換聯(lián)合壓制面波.畢云云等(2017)提出了基于離散曲波變換字典和二維局部離散余弦變換字典聯(lián)合壓制面波的方法.徐陽(yáng)等(2018)提出了廣義S變換與二維離散小波變換聯(lián)合壓制面波方法，該方法對(duì)包含面波成分的低頻信號(hào)再次進(jìn)行二維離散小波變換，成功提取出了低頻信息里包含的有效信號(hào)，減少了在有效壓制面波干擾的同時(shí)對(duì)有效反射波的損害.路鵬飛等(2020)提出了曲波變換與傅里葉變換聯(lián)合壓制面波方法.這些多手段聯(lián)合壓制面波的方法，在盡可能保護(hù)有效信號(hào)的原則下，最大限度地將有效信號(hào)和面波分離，進(jìn)一步提升了地震資料的信噪比.

本文首先在前人研究的基礎(chǔ)上將同步擠壓小波變換推廣至一維空間域，后根據(jù)面波干擾頻率低、視速度低的特點(diǎn)，采用高精度時(shí)頻分析方法同步擠壓小波時(shí)域和空域變換聯(lián)合壓制面波干擾.從時(shí)間、空間、頻率、波數(shù)維度對(duì)地震數(shù)據(jù)進(jìn)行分析，利用時(shí)域同步擠壓小波變換結(jié)合時(shí)頻域?yàn)V波器將面波與有效信號(hào)分離，并對(duì)分離出的面波數(shù)據(jù)體使用空域同步擠壓小波變換進(jìn)行低頻有效信號(hào)的重提取，通過(guò)理論模型測(cè)試驗(yàn)證了所提聯(lián)合方法的有效性，最后通過(guò)實(shí)際資料處理與一些常用方法的處理結(jié)果相比較，證明了所提出的聯(lián)合方法處理結(jié)果信噪比高、保幅性好的優(yōu)點(diǎn).

1 方法原理

1.1 一維空間域同步擠壓小波變換

同步擠壓小波變換(Synchrosqueezing Wavelet Transform, SWT)是由Daubechies等(2011)引入的一種針對(duì)非平穩(wěn)信號(hào)的稀疏表示.SWT是基于小波變換的時(shí)頻重分配，其類似于經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition, EMD)可以將信號(hào)分解為本征模態(tài)函數(shù)，優(yōu)于EMD的是，SWT具有更堅(jiān)實(shí)的理論基礎(chǔ)(Thakur et al.,2013).

同步擠壓算法通過(guò)沿頻率軸壓縮信號(hào)的能量來(lái)聚焦時(shí)頻表示，這種后處理方法的時(shí)頻分辨率極高且具有嚴(yán)格的連續(xù)可逆性.該變換可將時(shí)變信號(hào)分解為如式(1)所示的固有模態(tài)：

(1)

式中，Am(t)為瞬時(shí)振幅，η(t)代表噪聲或誤差，M為模態(tài)數(shù)量，φm(t)為瞬時(shí)相位.每個(gè)模態(tài)的瞬時(shí)頻率是通過(guò)瞬時(shí)相位的偏差(Boashash,2003)來(lái)計(jì)算的，計(jì)算方程為：

(2)

現(xiàn)進(jìn)行空間域同步擠壓小波變換推導(dǎo)，給定一空變信號(hào)s(x)：

s(x)=A(x)cos(kx+φ0)，(3)

式中，x為波動(dòng)距離，φ0為初相，A為振幅，(kx+φ0)為信號(hào)相位，k為空域角頻率也稱為角波數(shù)：

k=2πξ，(4)

式中，ξ為空間頻率.

若需要獲取空變信號(hào)的空間頻率即波數(shù)，通常對(duì)于空間采樣信號(hào)通過(guò)傅里葉變換來(lái)進(jìn)行空間-波數(shù)域變換，得到此空間采樣信號(hào)的波數(shù)譜，其一般定義式為：

(5)

給定母小波函數(shù)φ，對(duì)空變信號(hào)s(x)進(jìn)行連續(xù)小波變換，得到空間-尺度域小波系數(shù)Ws(a,b)(Stéphane,2009)：

(6)

式中a為尺度序列，b為空間偏移.根據(jù)Plancherel定理，式(6)可以寫成波數(shù)域表達(dá)形式：

(7)

(8)

對(duì)Ws(a,b)來(lái)說(shuō)，|Ws(a,b)|表示各空間模態(tài)分量的小波系數(shù)(Li and Liang,2012)：

(9)

和時(shí)間偏移τ一樣，當(dāng)給定一個(gè)空間偏移b，信號(hào)能量主要沿尺度軸擴(kuò)散.忽略沿空間軸的能量擴(kuò)散，對(duì)任意Ws(a,b)不為零的尺度-空間點(diǎn)(a,b)，其瞬時(shí)空間頻率估計(jì)表示為：

(10)

當(dāng)把地震信號(hào)同一時(shí)刻各道的數(shù)據(jù)序列看作是反映一定范圍內(nèi)地層及其地質(zhì)參數(shù)，便構(gòu)成了地震信號(hào)的空間序列.地震信號(hào)的空間序列與時(shí)/深剖面是一一對(duì)應(yīng)的，對(duì)于某一探區(qū)，不同時(shí)刻的空間序列信號(hào)對(duì)應(yīng)于相應(yīng)時(shí)/深剖面的不同深度地層位置(黨建武和黃建國(guó),2003).

其實(shí)際情況下的一般離散型式為：

(11)

式中，i為時(shí)序采樣點(diǎn)，q為地震記錄的道號(hào).

(12)

離散形式的瞬時(shí)頻率估計(jì)為：

(13)

需要注意的是，式(13)中得到的瞬時(shí)空間頻率估計(jì)往往并不是瞬時(shí)空間頻率本身(只有在輸入信號(hào)為理想簡(jiǎn)諧信號(hào)時(shí)才為瞬時(shí)頻率)，他是我們用來(lái)計(jì)算實(shí)際空間頻率的校正量.

基于上述概念推導(dǎo)出空間同步擠壓小波變換(Spatial Domain Synchrosqueezing Wavelet Transform, SSWT)的離散形式為：

(14)

這里把a(bǔ)和ω看成連續(xù)函數(shù)，空間同步擠壓小波變換的連續(xù)形式為：

(15)

式中，A(x)={a;Ws(a,x)≠0}.

式(14)、(15)將時(shí)頻重分配思想推廣至空間-波數(shù)域，且同步擠壓只發(fā)生在空間頻率軸.

空間域同步擠壓小波變換得到的同步擠壓空間波數(shù)表示，可以使用類似于時(shí)域同步擠壓小波變換的重構(gòu)方式，重構(gòu)得到原始信號(hào)：

(16)

如前文所述，野外采集得到的地震信號(hào)是時(shí)空維度的，通過(guò)相速度我們可以建立起地震時(shí)域序列與空域序列的聯(lián)系：

(17)

式中，k為角波數(shù)，Ω為角頻率，Cp為相速度.由于地震橫波如面波具有很強(qiáng)的頻散現(xiàn)象其傳播速度取相速度，地震體波的相速度和群速度基本相差不大，因此我們可以將式(17)中的相速度近似規(guī)定為地震波傳播速度，且地震波波速越大，其角波數(shù)越小.對(duì)于任意t時(shí)刻，有：

(18)

式中，Ω、φt為常數(shù).

信號(hào)s(x)的連續(xù)小波變換(Continuous wavelet transform, CWT)為：

(19)

瞬時(shí)空間頻率為：

(20)

信號(hào)s(x)的SSWT為：

(21)

通過(guò)式(19)、(21)建立起了CWT、SSWT與地震波波速之間的聯(lián)系，并將其展示在空間波數(shù)譜中.

由于面波具有視速度低的特性，其在炮集記錄中呈現(xiàn)近乎垂直的圓錐狀，與有效信號(hào)差異明顯，因此在對(duì)地震數(shù)據(jù)進(jìn)行空間域同步擠壓小波變換時(shí)，面波干擾將主要出現(xiàn)在高波數(shù)區(qū)，據(jù)此可選合適濾波方式對(duì)低頻段面波干擾和有效信號(hào)分離，恢復(fù)部分低頻有效信號(hào).

1.2 時(shí)域同步擠壓小波變換與空域同步擠壓小波變換聯(lián)合面波壓制方法

綜合上述方法原理，針對(duì)面波與有效信號(hào)的振幅、頻率、視速度差異，采用時(shí)域同步擠壓小波變換與空域同步擠壓小波變換聯(lián)合方法壓制面波干擾，具體實(shí)現(xiàn)步驟如下：

(1)對(duì)目標(biāo)單炮記錄進(jìn)行頻譜分析，確定有效信號(hào)和面波干擾的頻帶.選取最優(yōu)母小波函數(shù)，按道對(duì)輸入的單炮記錄進(jìn)行時(shí)域同步擠壓小波正變換(Time Domain Synchrosqueezing Wavelet Transform, TSWT)，得到每道信號(hào)對(duì)應(yīng)的時(shí)頻譜.

(2)根據(jù)面波和有效信號(hào)在時(shí)頻譜上能量顯示差異，面波能量位于低頻段沿時(shí)間軸聚集，并呈條帶狀分布，有效信號(hào)能量沿頻率軸聚集且沿頻率軸正方向收斂，同時(shí)結(jié)合頻譜分析得到的面波與有效信號(hào)各自的頻帶，在時(shí)頻譜上確定面波分布區(qū).

(3)選擇合適時(shí)頻窗口將面波區(qū)域內(nèi)的數(shù)據(jù)充零.后進(jìn)行時(shí)域同步擠壓小波反變換，得到去除面波干擾的單道時(shí)域信號(hào)，按道依次重復(fù)步驟1～3，即可得到去除面波干擾后的地震記錄.同時(shí)，將原始單炮記錄與此記錄匹配相減，得到一次濾波濾除的面波數(shù)據(jù)體.

(4)對(duì)分離出來(lái)的低頻數(shù)據(jù)體按時(shí)序進(jìn)行SSWT，選擇合適空頻窗口將面波區(qū)域內(nèi)的系數(shù)充零，反變換為各時(shí)刻點(diǎn)對(duì)應(yīng)的空變信號(hào)，即提取出在時(shí)頻域?yàn)V波中損傷的部分有效信號(hào).將其疊加到第3步得到的一次濾波后的地震記錄中，完成了對(duì)面波干擾的聯(lián)合壓制，有效保護(hù)了低頻段反射波信號(hào).

本文基于Gauss、Bump、Cmhat、Morlet四種母小波函數(shù)進(jìn)行了信號(hào)重構(gòu)(重構(gòu)誤差見(jiàn)圖1).結(jié)果表明Morlet小波的重構(gòu)誤差遠(yuǎn)小于其他小波，且能較好地表征信號(hào)時(shí)頻特征(圖2).因此本文選用Morlet小波作為母小波函數(shù).

圖1 四種母小波函數(shù)重構(gòu)誤差比較Fig.1 Comparison of reconstruction errors of four mother wavelet functions

圖2 四種母小波函數(shù)的TSWT時(shí)頻譜Fig.2 The time frequency spectrum of TSWT of four mother wavelet functions

2 理論模型測(cè)試

為驗(yàn)證所提方法的有效性，首先對(duì)正演得到的理論模型進(jìn)行處理，模型參數(shù)見(jiàn)表1.

表1 地層模型參數(shù)Table 1 Formation model parameters

圖3 理論模型合成記錄Fig.3 Theoretical model synthesis record

圖3為采用彈性波波場(chǎng)模擬方式生成的理論模型合成記錄，單震源單邊激發(fā)，震源位于地表，二維均勻各項(xiàng)同性介質(zhì)，共設(shè)計(jì)了2層水平反射層，采樣間隔為1 ms，共有750個(gè)采樣點(diǎn)，總記錄長(zhǎng)度為0.75 s，80道接收，最小偏移距0 m，道間距6.25 m，最大偏移距493.75 m.面波主頻為20 Hz，反射波主頻為50 Hz，整個(gè)合成記錄不含隨機(jī)噪聲.

圖4為理論模型第35道數(shù)據(jù)進(jìn)行TSWT后得到的時(shí)頻譜.時(shí)頻分析后可得，根據(jù)面波與有效信號(hào)能量聚集差異和頻率差異可以將其區(qū)分，但面波與有效信號(hào)在部分頻段存在重合.因此，當(dāng)進(jìn)行面波區(qū)充零處理時(shí)將損傷部分低頻有效信號(hào)(徐陽(yáng)等,2018).

圖4 理論模型合成記錄第35道的時(shí)頻譜Fig.4 The time frequency spectrum of channel 35 of theoretical model synthesis record

對(duì)理論模型合成記錄進(jìn)行TSWT，確定面波區(qū)域并將其充零，得到不含面波干擾的模型記錄，如圖5a所示.其與原始理論模型合成記錄匹配相減，得到一次濾波后濾出的面波干擾，如圖5b所示.圖5a與圖3對(duì)比發(fā)現(xiàn)，其反射波的同相軸振幅發(fā)生了改變，分析可得是由于充零過(guò)程導(dǎo)致低頻有效信號(hào)收到了損傷，這部分有效信號(hào)在圖5b中也得到了體現(xiàn).

圖5 (a) 經(jīng)TSWT壓制面波后的理論模型記錄； (b) TSWT分離出的面波記錄； (c) 經(jīng)SSWT重構(gòu)出的有效信號(hào)； (d) 經(jīng)聯(lián)合方法處理所得到的最終記錄Fig.5 (a) Record of theoretical model after suppressing ground roll by TSWT； (b) Ground roll separated by TSWT； (c) Valid signal reconstructed by SSWT； (d) Final record of theoretical model after suppressing ground roll by combination method

為恢復(fù)被損傷的低頻有效信號(hào)，對(duì)經(jīng)TSWT后分離出的面波數(shù)據(jù)體再次按時(shí)序?qū)臻g序列進(jìn)行SSWT，得到空間波數(shù)譜.如圖6所示為0.32 s時(shí)刻的空間波數(shù)譜對(duì)比.

可見(jiàn)對(duì)0.32 s的空變信號(hào)序列進(jìn)行同步擠壓小波變換后得到的空間波數(shù)譜較連續(xù)小波變換后能量沿波數(shù)軸有了明顯的集中，同步擠壓小波變換后的空頻分辨率遠(yuǎn)高于連續(xù)小波變換.如圖6c所示，偏移距在0～100 m能量較強(qiáng)且基本對(duì)應(yīng)同一波數(shù)，結(jié)合如圖5的輸入數(shù)據(jù)不難得到偏移距為0～100 m時(shí)有效信號(hào)同相軸接近水平，說(shuō)明其視速度相差不大，因此對(duì)應(yīng)在空間波數(shù)譜上的波數(shù)接近一致，而面波視速度較低且在同相軸不水平連續(xù)，因此在空間波數(shù)譜上位于高波數(shù)區(qū)域且能量軸傾斜發(fā)散.這正是在空間波數(shù)譜上區(qū)分面波和有效信號(hào)的依據(jù).對(duì)分離的面波數(shù)據(jù)體按時(shí)序進(jìn)行SSWT，抽取同樣時(shí)刻數(shù)據(jù)得到的空間波數(shù)譜如圖6f所示.在時(shí)頻域被損傷的有效信號(hào)即模型中的兩套水平層位的反射波部分能量對(duì)應(yīng)一波數(shù)軸刻度集中呈長(zhǎng)條狀分布，而面波干擾能量團(tuán)在空間波數(shù)譜中高波數(shù)區(qū)呈傾斜“堅(jiān)果”狀分布，如圖6e、f中紅框內(nèi)所示.圖6充分說(shuō)明了SSWT的進(jìn)行空頻分析的可行性和其較連續(xù)小波變換更卓越的空頻分辨能力.

將面波區(qū)充零，并重構(gòu)即可提取低頻段部分被損傷的有效信號(hào).據(jù)圖6所示，面波干擾得到有效壓制同時(shí)有效信號(hào)的同相軸連續(xù)性和形態(tài)得到極大改善，與理論模型合成記錄接近.同時(shí)，本文還使用TSWT對(duì)各方法的濾波結(jié)果進(jìn)行了單頻地震數(shù)據(jù)提取，提取出的四組數(shù)據(jù)從低頻至高頻涵蓋了面波與有效信號(hào)的優(yōu)勢(shì)頻帶，并計(jì)算了各目標(biāo)頻段以及整個(gè)炮集記錄的的信噪比(具體算法請(qǐng)見(jiàn)附錄).如圖7所示，對(duì)比原始記錄、TSWT濾波后記錄和使用聯(lián)合方法所得各有效頻段信噪比結(jié)果可得，經(jīng)聯(lián)合方法處理后得到的地震記錄的分頻記錄和整體記錄的信噪比均高于前者，尤其是在低頻段.以上結(jié)果證明了SSWT有效恢復(fù)在TSWT濾波過(guò)程中被損傷的部分反射波，TSWT和SSWT的聯(lián)合方法提升了地震記錄的信噪比.

圖6 0.32 s時(shí)刻信號(hào)空間波數(shù)譜對(duì)比(a) 面波分離后合成記錄的波數(shù)-歸一化振幅譜； (b) 面波分離后合成記錄經(jīng)空間連續(xù)小波變換后得到的空間波數(shù)譜； (c) 面波分離后合成記錄經(jīng)SSWT后得到的空間波數(shù)譜； (d) 分離出的面波記錄的波數(shù)-歸一化振幅譜； (e) 分離出的面波記錄經(jīng)空間連續(xù)小波變換后得到的空間波數(shù)譜； (f) 分離出的面波記錄經(jīng)SSWT后得到的空間波數(shù)譜.Fig.6 Comparison of spatial wavenumber spectrum of signal at 0.32 s(a) The wavenumber-normalized amplitude spectrum of the synthetic record after ground roll separation； (b) The spatial wavenumber spectrum of the synthetic record after ground roll separation by Spatial Domain Continuous Wavelet Transform (SCWT)； (c) The spatial wavenumber spectrum of the synthetic record after ground roll separation by SSWT； (d) The wavenumber-normalized amplitude spectrum of the ground roll； (e) The spatial wavenumber spectrum of the ground roll by SCWT； (f) The spatial wavenumber spectrum of the ground roll by SSWT.

圖8 聯(lián)合方法面波壓制效果對(duì)比圖(a) Alaska地區(qū)第20炮原始記錄； (b) 時(shí)域同步擠壓小波變換壓制面波干擾后的記錄； (c) 使用聯(lián)合方法壓制面波干擾后的記錄.Fig.8 Comparison of ground roll suppression effects of combined method(a) The original record of the 20th source in Alaska；(b) The record after ground roll suppression by TSWT；(c) The record after ground roll suppression by combined method.

圖9 SSWT的輸入與輸出(a) TSWT后濾除的低頻記錄； (b) SSWT重構(gòu)結(jié)果.Fig.9 Input and output of SSWT(a) Low frequency record filtered byTSWT；(b) Reconstruction result of SSWT.

圖10 時(shí)空域同步擠壓小波變換切片(a) TSWT濾波前15 Hz頻率切片； (b) TSWT濾波后15 Hz頻率切片； (c) SSWT濾波前0.001 rad·m-1波數(shù)切片； (d) SSWT濾波后0.001 rad·m-1波數(shù)切片.Fig.10 Slices of SSWT and TSWT(a) 15 Hz frequency slice before TSWT； (b) 15 Hz frequency slice after TSWT； (c) 0.001 rad·m-1 wavenumber slice before SSWT； (d) 0.001 rad·m-1 wavenumber slice after SSWT.

圖11 原始單炮記錄第30道TSWT時(shí)頻譜和分離出的低頻數(shù)據(jù)體第320 ms SSWT空間波數(shù)譜(a)和(c) 面波區(qū)充零前； (b)和(d) 面波區(qū)充零后.Fig.11 The time frequency spectrum of channel 30 of TSWT and spatial wavenumber spectrum at 320 ms of SSWT of the separated low-frequency data(a) and (c) Before zero filling in the ground roll area； (b) and (d) After zero filling in the ground roll area.

圖12 3組數(shù)據(jù)的頻譜對(duì)比Fig.12 Spectrum comparison of three datasets

3 實(shí)際資料處理

本文選取美國(guó)Alaska地區(qū)某一測(cè)線實(shí)際數(shù)據(jù)進(jìn)行處理，震源采用炸藥震源，總道數(shù)5568道(曾祥堃，2015).抽取第20炮作為輸入單炮記錄，每炮共96道，道間距33.53 m，最小偏移距55 m，最大偏移距5525 m，3000個(gè)采樣點(diǎn)，采樣間隔2 ms，記錄長(zhǎng)度5 s，如圖8a所示，面波呈“圓錐狀”存在，近偏移距的淺層反射波和遠(yuǎn)偏移距中、深層反射波與面波干擾重疊，有效信號(hào)難以識(shí)別，信噪比嚴(yán)重降低.按道依次進(jìn)行TSWT，得到各道的時(shí)頻譜，結(jié)合數(shù)據(jù)分析結(jié)果，確定面波分布區(qū)和時(shí)頻濾波窗口將其充零.整個(gè)記錄的面波分離完成后，將得到的不含面波干擾的地震記錄(圖8b)，與原始記錄匹配相減，即可得到經(jīng)TSWT后濾除的低頻數(shù)據(jù)體，如圖9a所示.低頻數(shù)據(jù)體中含有若干條反射波同相軸，結(jié)合圖10a、b的TSWT濾波前后的15Hz頻率切片來(lái)看，分析該現(xiàn)象是由于低頻段有效信號(hào)與面波存在重疊，且為徹底壓制面波干擾而擴(kuò)大使用時(shí)頻濾波窗口等綜合原因?qū)е?

為減少對(duì)有效信號(hào)的損害，將時(shí)域同步擠壓小波變換分離出的低頻數(shù)據(jù)體按時(shí)序進(jìn)行SSWT，得到空間波數(shù)譜，此時(shí)空域的地震數(shù)據(jù)將轉(zhuǎn)換至空間、波數(shù)域中.在空間波數(shù)譜中確定面波區(qū)將其充零，如圖11所示.圖10c、d為SSWT濾波前后空頻濾波窗口內(nèi)0.001 rad·m-1波數(shù)切片，濾波前對(duì)應(yīng)的重構(gòu)結(jié)果未混有有效信號(hào)，濾波后該波數(shù)對(duì)應(yīng)的面波干擾基本被完全壓制，證實(shí)了有效信號(hào)與面波干擾視速度差異的SSWT方法的優(yōu)越性.

有效信號(hào)提取結(jié)果如圖9b所示，可見(jiàn)明顯同相軸.將重構(gòu)結(jié)果疊加到TSWT壓制面波干擾后的記錄中，得到聯(lián)合面波壓制方法的最終結(jié)果，如圖8c所示.圖8c中反射波能量尤其是面波覆蓋區(qū)和中深層段，顯著增強(qiáng).通過(guò)圖12中3組數(shù)據(jù)頻譜對(duì)比，可以更直觀地看出使用聯(lián)合方法后，低頻段能量被有效恢復(fù).根據(jù)圖13信噪比對(duì)比結(jié)果來(lái)看，經(jīng)聯(lián)合方法處理后的低頻段數(shù)據(jù)以及整體地震數(shù)據(jù)的信噪比均高于原始數(shù)據(jù)和單一使用TSWT處理后的結(jié)果.從實(shí)際資料的處理結(jié)果可以看出本文所述面波壓制聯(lián)合方法在有效壓制面波的同時(shí)能極大程度地保留低頻段有效信號(hào).

4 討論

在這一部分中，主要討論時(shí)域和空域同步擠壓小波變換聯(lián)合方法的優(yōu)勢(shì)和通過(guò)實(shí)際結(jié)果呈現(xiàn)出的該方法尚未解決的缺點(diǎn)以及未來(lái)的研究方向.

4.1 聯(lián)合方法實(shí)際應(yīng)用優(yōu)勢(shì)

時(shí)域和空域同步擠壓小波變換聯(lián)合面波壓制方法是高精度時(shí)頻分析方法在地震資料處理面波壓制領(lǐng)域的綜合應(yīng)用.如前文分析，理論上使用TSWT一次濾波后可以將面波徹底分離，再對(duì)提取出的面波進(jìn)行SSWT，重提取部分低頻有效信號(hào).對(duì)于實(shí)際數(shù)據(jù)時(shí)頻分析來(lái)說(shuō)，精確定位面波區(qū)，徹底分離面波干擾，避免過(guò)多損傷非目標(biāo)頻段的有效信號(hào)才是主要目的.

時(shí)域同步擠壓算法可以有效解決傳統(tǒng)時(shí)頻分析方法能量沿頻率軸發(fā)散的問(wèn)題.雖然未使用同步擠壓算法的時(shí)頻分析方法得到的時(shí)頻譜看起來(lái)能量成團(tuán)分布，但是相鄰能量出現(xiàn)多處重疊，這導(dǎo)致以往的時(shí)頻域?yàn)V波方法為了徹底分離面波，存在擴(kuò)大濾波窗口(系數(shù)充零區(qū))的現(xiàn)象，會(huì)損傷非目標(biāo)頻段的有效信號(hào)，而TSWT會(huì)有效改善這一問(wèn)題，但仍未徹底解決.因此時(shí)頻窗口選擇范圍越精確，信號(hào)的保幅性越好，雖然這會(huì)犧牲一定的計(jì)算效率.

本文所提出的SSWT重提取在一次濾波中被濾除的有效信號(hào)，充分利用了反射波和面波的視速度差異.結(jié)合圖5c和圖9b的重構(gòu)結(jié)果，是聯(lián)合方法中SSWT可以恢復(fù)部分被損傷的低頻有效信號(hào)的直接證據(jù)，這一優(yōu)勢(shì)是單一方法和其他聯(lián)合面波壓制方法不具備的.圖14的四種方法處理結(jié)果對(duì)比再次證實(shí)了上述觀點(diǎn).圖14a中可見(jiàn)基于視速度差異的FK濾波對(duì)面波壓制效果較為徹底但是會(huì)損傷淺層有效信號(hào)如直達(dá)波，同時(shí)會(huì)產(chǎn)生大量空間假頻；圖14b和圖14c對(duì)比證明了TSWT的面波壓制能力和保幅能力優(yōu)于CWT但還是丟失了部分中深層有效信號(hào)；圖14d使用聯(lián)合方法得到的處理結(jié)果面波壓制力和保幅性兼優(yōu)，壓制效果最為理想.

4.2 聯(lián)合方法的缺點(diǎn)及未來(lái)研究方向

如上述分析在使用同步擠壓小波變換后對(duì)時(shí)頻、空頻濾波窗口內(nèi)系數(shù)充零時(shí)沒(méi)有尋找到最優(yōu)處理方式.本研究下一步可以改進(jìn)濾波方式，在閾值迭代去噪、自適應(yīng)時(shí)窗、深度學(xué)習(xí)等領(lǐng)域進(jìn)行探索.

在SSWT重構(gòu)有效信號(hào)環(huán)節(jié)，主要依據(jù)的是面波和有效信號(hào)視速度的差異即在地震剖面中同相軸傾斜程度不同.然而由于面波的頻散效應(yīng)部分相鄰道往往會(huì)形成低傾角同相軸，這些假層位同相軸振幅往往較小，但會(huì)不可避免的出現(xiàn)在重構(gòu)結(jié)果中，如圖5c所示.如何去除這些假層位低傾角同相軸是未來(lái)研究的重點(diǎn)問(wèn)題.

在圖13中45 Hz的信噪比對(duì)比中聯(lián)合方法處理后的地震數(shù)據(jù)信噪比略低于僅使用TSWT處理.我們抽取了三組45 Hz地震數(shù)據(jù)及其時(shí)頻能量切片，發(fā)現(xiàn)不存在明顯的面波干擾，如圖15b、e所示，因此判斷是由于重構(gòu)誤差或是信噪比算法誤差所致.聯(lián)合方法處理得到的整炮地震記錄信噪比是最高的未與理論相背.

5 結(jié)論

本文將同步擠壓小波變換推廣至空間域，并對(duì)其實(shí)現(xiàn)過(guò)程進(jìn)行了公式推導(dǎo)與實(shí)際驗(yàn)證.并通過(guò)時(shí)域與空域同步擠壓小波變換聯(lián)合方法壓制面波干擾，并對(duì)理論和實(shí)際資料進(jìn)行處理，調(diào)整最優(yōu)參數(shù)，取得了較好的處理效果.該聯(lián)合方法的優(yōu)點(diǎn)在于：

(1)使用高精度時(shí)頻分析方法時(shí)域同步擠壓小波變換，結(jié)合時(shí)頻域?yàn)V波與空頻域?yàn)V波，對(duì)面波壓制徹底.

(2)對(duì)一次濾波后分離的面波干擾再次進(jìn)行空域同步擠壓小波變換，從空間、波數(shù)域中判別有效信息與面波干擾的地震屬性差異，能有效提取部分低頻反射波，切實(shí)保護(hù)了低頻段有效信息，提高了地震數(shù)據(jù)的信噪比.

(3)本文方法的實(shí)現(xiàn)，證明了在面波壓制中進(jìn)行時(shí)空域?yàn)V波，從而恢復(fù)低頻段有效信號(hào)的方法的優(yōu)越性.

圖13 野外單炮記錄各有效頻段信噪比對(duì)比Fig.13 Comparison of signal-to-noise ratio of each valid frequency band recorded by actualdata

圖14 四種面波壓制方法實(shí)際資料處理結(jié)果對(duì)比(a) FK濾波； (b) CWT； (c) TSWT； (d) TSWT & SSWT.Fig.14 Comparison of actual data processing results of four ground roll suppression methods(a) FK filtering； (b) CWT； (c) TSWT； (d) TSWT & SSWT.

圖15 實(shí)際記錄45 Hz切片(a) TSWT 45 Hz 切片； (b) SSWT重構(gòu)結(jié)果45 Hz切片； (c) 聯(lián)合方法最終記錄45 Hz切片； (d) TSWT 45 Hz時(shí)頻譜； (e) SSWT重構(gòu)結(jié)果 45 Hz時(shí)頻譜； (f) 聯(lián)合方法最終記錄45 Hz時(shí)頻譜.Fig.15 45 Hz slices of actual record(a) TSWT； (b) Reconstruction result of SSWT； (c) Final record with combination method； (d) The time frequency spectrum of TSWT at 45 Hz； (e) The time frequency spectrum of reconstruction result of SSWT at 45 Hz； (f) The time frequency spectrum of final record with combination method at 45 Hz.

附錄

本文計(jì)算信噪比的方法為相關(guān)時(shí)間移法及其改進(jìn)算法(牛聰?shù)?2006;王紅玲,2007).

Qi,i+1(ti,i+1)=max{Qi,i+1(τ)}，(A1)

其中相鄰道間最大的互相關(guān)值Qi,i+1(ti,i+1)作為信號(hào)的功率譜：

(A2)

(A3)

N為地震道數(shù)，使用各道的自相關(guān)值作為地震記錄的總能量，則噪聲的能量為：

(A4)

信噪比為：

(A5)

當(dāng)對(duì)含有面波干擾的原始地震記錄進(jìn)行信噪比估算時(shí)，需要將使用聯(lián)合方法濾波后的最終記錄作為有效信號(hào)能量Es的輸入記錄.總能量E則用原始地震記錄通過(guò)相關(guān)時(shí)移法計(jì)算.