袁 力

(蕪湖職業(yè)技術(shù)學(xué)院 基礎(chǔ)部, 安徽 蕪湖 241003)

面波是地震信號中一種主要的規(guī)則干擾,降低了信噪比,影響信號處理質(zhì)量.地震記錄中的面波一般呈扇形,有能量高、頻率低等特點.傳統(tǒng)面波壓制方法(如F-K濾波[1]、τ-p變換[2]等)主要利用面波與有效信號在頻帶和視速度方面的不同特征,但壓制面波的同時有效信息也會被壓制.進一步,根據(jù)面波特征在地震記錄中估計面波,再在地震記錄中減去估計面波,方法有維納濾波[3]、K-L變換[4]等.2012年Freire等利用奇異值分解(singular value decomposition,SVD)分離地震記錄剖面的上行波和下行波[5].

20世紀(jì)80年代以來,基于小波變換的面波壓制方法不斷發(fā)展.1996年羅國安等利用面波能量高等特點,在小波域中根據(jù)面波的視速度特點作線性時移使面波逐道相干,再作K-L分解估計面波[6].2004年詹毅等將小波包與SVD結(jié)合,在小波包處理后聯(lián)合自動追蹤同相軸與SVD進行波場分離[7].2013年Boustani等在曲波域利用自適應(yīng)SVD壓制面波,將曲波域按面波能量比例分為不同區(qū)域,不同區(qū)域采用不同面波壓制方法[8].2020年路鵬飛等將含面波地震信號變換到曲波域,對曲波域系數(shù)作傅里葉變換,在此進行面波壓制,最后作傅里葉逆變換、曲波逆變換得到壓制面波后的地震信號[9].2021年孟娟等利用小波譜能量曲線、經(jīng)驗小波變換等方法壓制面波,首先對信號作連續(xù)小波變換,根據(jù)小波譜得到能量曲線,結(jié)合極大值點頻率和鄰域法確定分割邊界,再作經(jīng)驗小波變換得到各固有模態(tài)分量,對根據(jù)各子頻帶的頻率能量確定的面波模態(tài)分量作濾波,從而分離面波和有效波[10].

1998年Kingsbury提出雙樹復(fù)小波變換,增強了傳統(tǒng)小波變換,廣泛應(yīng)用于特征提取和數(shù)據(jù)分析等領(lǐng)域[11].地震信號有多尺度性,有效信號之間有較強相關(guān)性,在面波干擾下,反射波能量相對較弱,影響反射波的連續(xù)性.雙樹復(fù)小波變換能有效提取和表征地震信號的變化和高維奇異特點.作為常用數(shù)據(jù)降維方法,主成分分析(principal component analysis,PCA)對協(xié)方差陣作特征分解,所得特征向量對應(yīng)數(shù)據(jù)的各個主成分,特征值對應(yīng)數(shù)據(jù)在各個主成分的權(quán)重.PCA變換后各主成分彼此不相關(guān),編號越大的主成分所含信息越少.鑒于此,本文采用二維雙樹復(fù)小波變換與主成分分析進行面波壓制.仿真和實際實驗處理結(jié)果表明本文方法能有效壓制面波,保留較多的有效信號信息.

1 雙樹復(fù)小波域主成分分析面波壓制方法

1.1 雙樹復(fù)小波變換

可分離小波信號表示其細(xì)節(jié)特征不好,方向選擇少,雙樹復(fù)小波變換(dual-tree complex wavelet transform,DTCWT)解決了此問題.雙樹復(fù)小波具有近似平移不變性,其系數(shù)的模有旋轉(zhuǎn)不變性[11].雙樹復(fù)小波分解含兩個平行分解樹,用兩個獨立的實小波分別計算實部和虛部系數(shù).實部和虛部濾波器滿足希爾伯特變換,實際分解變換中采用雙樹算法.

圖1為一維雙樹復(fù)小波分解變換示意圖,其中h0(n)和h1(n),g0(n)和g1(n)分別是樹A、樹B分解過程的低通和高通濾波器,分別產(chǎn)生低頻和高頻的實部、虛部小波系數(shù),↓2為下采樣.

圖1 一維雙樹復(fù)小波分解變換示意圖

對一維雙樹復(fù)小波作張量積得二維雙樹復(fù)小波變換,即對行和列分別進行一維雙樹復(fù)小波變換.一維雙樹復(fù)小波變換產(chǎn)生2∶1冗余度,而對n維信號,產(chǎn)生2n∶1冗余度.二維雙樹復(fù)小波分解時,各尺度產(chǎn)生2個低頻子帶和方向分別為±15°、±45°和±75°的6個不同高頻子帶,均有實部和虛部兩部分,對上一尺度的低頻子帶再進行分解可得多尺度分解,產(chǎn)生不同尺度、不同方向的小波系數(shù)[12].對含面波的仿真地震信號作二維雙樹復(fù)小波變換,分解層數(shù)為1.分解得到的所有子帶系數(shù)見圖2.

圖2 雙樹復(fù)小波變換后仿真信號的各子帶系數(shù)

圖2a為含面波的仿真地震信號,圖2b為經(jīng)雙樹復(fù)小波變換后仿真地震信號的低頻子帶實、虛部小波系數(shù),圖2c為±15°、±45°、±75°各方向高頻實、虛部子帶小波系數(shù).

圖2b顯示低頻子帶中含大量的反射波信息,但同時也含大量的面波信息,因此在該子帶應(yīng)采用合適的方法進行面波壓制.在圖2c中,±75°實、虛部子帶中反射波信息少,而面波成分多,考慮在該子帶應(yīng)用閾值方法壓制面波;±45°實、虛部子帶中同時含有反射波和面波信息,信息量均較弱,為了更好地壓制面波,在該子帶考慮進行面波壓制;±15°實、虛部子帶主要包含反射波信息,面波信息很弱,所以在這一子帶考慮保持不變.

1.2 主成分分析(PCA)

主成分分析是一種常用的統(tǒng)計分析方法[13],是對標(biāo)準(zhǔn)化后數(shù)據(jù)的協(xié)方差陣作特征分解,分解后的特征向量即為數(shù)據(jù)的主成分,各主成分分量之間不相關(guān),特征值越大的主成分所含信息越多,選取特征值較大的幾個特征向量,其余元素置零,并取該矩陣的轉(zhuǎn)置作為一個乘法因子,將該乘法因子乘以標(biāo)準(zhǔn)化后數(shù)據(jù)矩陣就可以得到重構(gòu)后的新數(shù)據(jù).

對地震信號進行雙樹復(fù)小波變換后的小波系數(shù)作PCA,由于其中編號大的主成分中含反射波的大部分特征,反射波信息主要包含在前幾個主成分分量中,其余主成分主要含噪聲信息,利用PCA壓制其他信號能量,較好地保留了有效信號信息.

為了說明PCA應(yīng)用于壓制小波系數(shù)面波成分的有效性,我們對低頻子帶系數(shù)進行PCA處理,利用PCA處理前、后經(jīng)雙樹復(fù)小波域變換后合成地震信號的低頻子帶系數(shù)和經(jīng)雙樹復(fù)小波域變換后反射波的低頻子帶系數(shù)計算信噪比,且選擇保留的主成分個數(shù)分別為2,6,10,12,16,結(jié)果見表1.

表1 低頻子帶系數(shù)PCA處理前、后的信噪比 dB

觀察表1數(shù)據(jù)可知,經(jīng)PCA處理后雙樹復(fù)小波域低頻子帶系數(shù)信噪比大多情況下得到大幅提升.PCA處理時,后面的成分仍含有用信息,為了減少重構(gòu)后的信息損失,選擇合適的主成分個數(shù)不可忽視.根據(jù)實驗結(jié)果我們發(fā)現(xiàn)保留主成分個數(shù)為10左右時,信噪比有峰值產(chǎn)生.

為了更直觀地顯示,我們將圖2中低頻子帶系數(shù)和高頻±45°實、虛部子帶系數(shù)進行PCA處理,選擇前10個主成分重構(gòu),結(jié)果見圖3.

圖3 低頻和高頻±45°實虛部子帶系數(shù)PCA處理結(jié)果圖

從圖3可以看出,經(jīng)PCA處理后、小波子帶系數(shù)中的面波成分得到較好的壓制,反射波信息得以保留.本文主要利用PCA變換消除小波系數(shù)間的相關(guān)性,去除冗余信息,保留有效信息.

1.3 雙樹復(fù)小波域PCA面波壓制方法

雙樹復(fù)小波變換具有表示信號高維奇異特征的PCA能力,而PCA是一種較好的特征提取方法,所以本文將兩者結(jié)合,基于雙樹復(fù)小波,對不同子帶區(qū)域采用不同方法壓制面波.首先對含面波信號作雙樹復(fù)小波變換.由圖2知,面波信號主要集中于低頻子帶和±45°、±75°方向的高頻子帶中.因此,對低頻和±45°高頻子帶實、虛部系數(shù)作PCA處理壓制其中的面波成分;對±75°高頻子帶實、虛部系數(shù)作閾值處理;對±15°高頻子帶實、虛部系數(shù)不作處理,旨在壓制面波的同時較多地保留有效信號.

閾值方法是經(jīng)典降噪方法,閾值函數(shù)各有優(yōu)、缺點[14].基于分位數(shù)的閾值公式[15]如下:

P(x≤T)=α,

其中:P(·)表示概率,α的取值范圍為[0.1,0.3].

基于雙樹復(fù)小波PCA的面波壓制方法步驟如下:

步驟1 對含面波的地震信號進行二維雙樹復(fù)小波變換,得到低頻和不同尺度方向的高頻實、虛部小波系數(shù);

步驟2 對分解得到的低頻實、虛部系數(shù)均作PCA處理,對各尺度±45°高頻實、虛部系數(shù)作PCA、±75°高頻實、虛部系數(shù)作基于分位數(shù)的閾值處理,±15°高頻實、虛部系數(shù)保持不變,得到壓制面波后的小波系數(shù);

步驟3 對處理后的小波系數(shù)作二維雙樹復(fù)小波逆變換,得到面波壓制后的地震信號.

2 實驗

2.1 仿真實驗

為了測試本文方法的效果,用Ricker子波模擬二維地震信號,在理想地震信號中加入面波信號,得到含面波的合成仿真地震信號.因?qū)嶋H地震信號往往含隨機噪聲,故仿真實驗中的模擬地震信號含面波、反射波和隨機噪聲.用零均值加性高斯白噪聲模擬隨機噪聲.本文實驗均采用MATLAB軟件進行編程處理.

為了驗證本文方法(DTCWT_PCA_TH)的有效性,將本文方法與PCA、DWT_PCA、DTCWT_SVD_TH、DTCWT_PCA方法進行對比實驗:

(ⅰ)將對合成地震信號直接進行主成分分析的面波壓制方法記為PCA.

(ⅱ)將對低頻和水平方向高頻子帶作PCA,其余系數(shù)不變的普通小波面波壓制方法記為DWT_PCA.

(ⅲ)將低頻系數(shù)均用奇異值分解(SVD),各尺度分解得到的±45°高頻實、虛部系數(shù)作SVD,±75°高頻實、虛部系數(shù)作基于分位數(shù)的閾值處理,±15°方向的高頻實、虛部系數(shù)保持不變的面波壓制方法記為DTCWT_SVD_TH.

(ⅳ)將低頻系數(shù)均作PCA,各尺度分解得到的±45°和±75°高頻實、虛部系數(shù)作PCA,±15°高頻實、虛部系數(shù)保持不變的面波壓制方法記為DTCWT_PCA.

仿真實驗中,普通小波基函數(shù)為Symlets小波,分解層數(shù)為2,雙樹復(fù)小波分解層數(shù)為1.利用信噪比(SNR,dB)衡量各面波壓制方法的效果,實驗結(jié)果見圖4.

a 壓制面波前SNR=-0.801 6 dB

b 壓制面波前SNR=0.556 4 dB圖4 保留不同成分個數(shù)時各方法壓制面波后信噪比對比

圖4a、圖4b分別為加入標(biāo)準(zhǔn)差為15和10的高斯白噪聲的各類方法壓制面波后的信噪比折線圖.橫坐標(biāo)為PCA和SVD重構(gòu)時保留的主成分個數(shù),縱坐標(biāo)為信噪比,其中圖4a中壓制面波前信噪比為-0.801 6 dB,閾值公式中α=0.1.圖4b中壓制面波前信噪比為0.556 4 dB,閾值公式中α=0.2.由圖4可知:

(ⅰ)無閾值的DTCWT_PCA方法和本文方法在兩種不同的隨機噪聲強度下,在一定程度上優(yōu)于直接對原始數(shù)據(jù)進行處理的PCA方法和普通小波域PCA方法.

(ⅱ)本文方法和直接對原始數(shù)據(jù)進行處理的PCA方法其保留主成分個數(shù)需取合適大小,實驗中保留個數(shù)取為[8,12].

(ⅲ)無閾值的DTCWT_PCA方法保留主成分個數(shù)過多信噪比反而下降,原因應(yīng)在于面波壓制過少.

(ⅳ)本文方法優(yōu)于DTCWT_SVD_TH方法,說明雙樹復(fù)小波域的主成分分析閾值方法優(yōu)于雙樹復(fù)小波域的奇異值分解閾值方法.

地震學(xué)里常采用頻率波數(shù)(FK)方法提取頻散曲線,將信號從時間空間域變換到頻率波數(shù)域,本文采用二維傅里葉變換.再將低頻部分移到中心,與頻率和波數(shù)對應(yīng)起來,由于離散傅里葉變換譜的對稱性,只需取其一半圖像,再將負(fù)的波數(shù)軸取絕對值,得到FK譜.為了更好地顯示本文方法的有效性,將PCA、DWT_ PCA、DTCWT_SVD_TH、DTCWT_PCA方法與本文方法壓制含面波的合成地震信號的效果進行展示(圖5),其中降噪前信噪比為-0.801 6 dB,PCA與SVD重構(gòu)中保留的主成分均為前10個.

圖5 合成地震信號面波壓制結(jié)果圖及對應(yīng)FK譜圖

圖5為各方法對合成地震信號的面波壓制效果及對應(yīng)FK譜,其中圖5a為含隨機噪聲的合成地震信號及對應(yīng)FK譜,圖5b為用本文方法壓制面波的效果及其對應(yīng)FK譜圖,圖5c～圖5f分別為用DWT_PCA,DTCWT_PCA,DTCWT_SVD_TH和PCA方法壓制面波的效果及其對應(yīng)FK譜圖.從圖5c、圖5d可以看出DWT_PCA,DTCWT_PCA兩種方法的面波壓制效果不明顯,仍有大量面波信息被保留.從圖5e、圖5f可以看出DTCWT_SVD_TH,PCA兩種方法的面波壓制效果明顯,但反射波也被壓制,有效信號損失過多.從圖5b可以看出本文方法較好地壓制了面波,且反射波較完整清晰.綜合圖4、圖5,說明本文方法能有效壓制面波,保留有效信號的較多能量和特征.

2.2 實際地震信號面波壓制實驗

在實際地震信號中,面波具有低頻率、低視速度、高振幅等特征和頻散特性.為了進一步說明本文方法壓制面波的性能較好,用DWT_PCA,DTCWT_SVD_TH和本文方法對含面波的實際地震信號壓制面波,壓制面波后的效果見圖6.

圖6 實際地震信號面波壓制效果圖

圖6a為原始實際地震信號,圖6b為用DWT_PCA方法壓制面波的效果,圖6c為用DTCWT_SVD_TH方法壓制面波的效果,圖6d為用本文方法壓制面波的效果.圖6d信號比圖6b、圖6c更清晰,面波壓制效果更好,說明本文方法面波壓制效果較好,大多數(shù)面波被濾除,面波壓制后有效信號的能量與質(zhì)量增強,波形變得清晰,該方法在壓制面波的同時還去除了較多的隨機噪聲.

3 結(jié)論

雙樹復(fù)小波域的主成分分析方法可以有效壓制面波,而閾值方法作為經(jīng)典降噪方法雖在閾值選取上有其局限性,但本文采用分位數(shù)作為閾值選取標(biāo)準(zhǔn),取得了一定效果.本文在雙樹復(fù)小波域,結(jié)合主成分分析方法,基于分位數(shù)的閾值降噪方法,有針對性地在不同子帶運用兩種方法,壓制面波信號的同時較好地保留了有效信號.本文還有待改進之處,例如,在雙樹復(fù)小波域可以考慮采用一些更有效的估計有效信號數(shù)學(xué)模型的方法,如擬合方法等,以期更好地壓制面波,恢復(fù)有效信號.