地震信號(hào)隨機(jī)噪聲壓制的雙樹(shù)復(fù)小波域雙變量方法

汪金菊，袁力，劉婉如，徐小紅

合肥工業(yè)大學(xué)數(shù)學(xué)學(xué)院，合肥　230009

?

地震信號(hào)隨機(jī)噪聲壓制的雙樹(shù)復(fù)小波域雙變量方法

汪金菊，袁力，劉婉如，徐小紅

合肥工業(yè)大學(xué)數(shù)學(xué)學(xué)院，合肥230009

有效地壓制地震信號(hào)中的噪聲是地震信號(hào)解釋和后續(xù)處理的重要環(huán)節(jié)之一.本文建立兩種雙樹(shù)復(fù)小波域雙變量模型對(duì)地震信號(hào)中的隨機(jī)噪聲進(jìn)行壓制.地震信號(hào)經(jīng)雙樹(shù)復(fù)小波變換后，同一方向?qū)嵅颗c虛部系數(shù)、實(shí)部(或虛部)系數(shù)與對(duì)應(yīng)的模之間存在較強(qiáng)的相關(guān)性.鑒于此，對(duì)同一方向?qū)嵅颗c虛部小波系數(shù)建立雙變量模型，從含噪地震信號(hào)小波系數(shù)中估計(jì)原始信號(hào)的小波系數(shù)，再基于雙樹(shù)復(fù)小波逆變換重構(gòu)得到降噪后的地震信號(hào).進(jìn)一步對(duì)同一方向?qū)嵅?或虛部)系數(shù)與對(duì)應(yīng)的模建立雙變量模型，得到地震信號(hào)隨機(jī)噪聲壓制的第二種雙樹(shù)復(fù)小波域雙變量方法.最后對(duì)合成地震記錄和實(shí)際地震資料中的隨機(jī)噪聲進(jìn)行壓制的實(shí)驗(yàn)結(jié)果證實(shí)本文兩種方法都能夠有效地壓制地震信號(hào)中的隨機(jī)噪聲.

雙樹(shù)復(fù)小波變換;隨機(jī)噪聲;雙變量模型;小波系數(shù)

1　引言

地震信號(hào)降噪的目的在于去除各種干擾，提高信噪比的同時(shí)較好地保留有效信號(hào).隨機(jī)噪聲衰減是地震信號(hào)處理中的一個(gè)重要環(huán)節(jié)，對(duì)噪聲的壓制效果直接影響后續(xù)處理.為了提高信噪比，人們根據(jù)信號(hào)與噪聲的特征差異，設(shè)計(jì)出各種降噪方法.以傅氏變換為基礎(chǔ)的頻率域?yàn)V波方法是地震資料處理中的經(jīng)典方法，該方法根據(jù)有效波和干擾波在頻帶上的差異來(lái)壓制干擾波突出有效波.該方法的不足在于信號(hào)的降噪效果依賴于信號(hào)頻帶和噪聲頻帶分離的程度.若混入噪聲的頻帶和信號(hào)頻帶有重疊時(shí),則降噪效果不是十分理想.在此基礎(chǔ)上，基于傅氏變換的f-x域預(yù)測(cè)地震信號(hào)降噪技術(shù)被提出，該技術(shù)可有效地壓制隨機(jī)噪聲，但不能對(duì)非平穩(wěn)地震信號(hào)進(jìn)行很好的處理.

小波是一種具有多分辨率和時(shí)頻局域性的分析工具，可將信號(hào)的特征表現(xiàn)在不同分辨層次上，信號(hào)的局部特征可以通過(guò)處理變換系數(shù)而改變，因而被廣泛應(yīng)用于信號(hào)降噪.一些學(xué)者利用小波變換這一時(shí)頻分析工具解決非穩(wěn)態(tài)問(wèn)題，相繼提出了二進(jìn)小波變換方法的地震信號(hào)分時(shí)分頻去噪(章珂等，1996)、小波包變換和f-x域預(yù)測(cè)相結(jié)合的降噪算法(宗濤等，1998)、基于小波域貝葉斯方法(張波等，1999)和最小二乘降噪算法(王振國(guó)等，2002).這類降噪方法首先將地震信號(hào)進(jìn)行多尺度分解得到高低頻子帶系數(shù)，然后對(duì)包含隨機(jī)噪聲的子帶系數(shù)進(jìn)行處理，以達(dá)到降噪目的.它們?cè)谝欢ǔ潭壬咸岣吡私翟胄Ч?，但未考慮小波域統(tǒng)計(jì)模型.基于拉普拉斯先驗(yàn)分布的軟閾值法(Berkner and Wells,2002)開(kāi)創(chuàng)性地分析了小波系數(shù)的統(tǒng)計(jì)特征.該降噪算法將拉普拉斯分布作為單個(gè)小波分解系數(shù)的統(tǒng)計(jì)分布，通過(guò)最大后驗(yàn)概率估計(jì)得到軟閾值函數(shù).以單變量概率密度函數(shù)作為先驗(yàn)分布求解閾值的單變量模型算法如軟閾值法，由于只考慮子帶內(nèi)單個(gè)小波系數(shù)的統(tǒng)計(jì)特征，沒(méi)有考慮子帶間小波系數(shù)的相關(guān)性，雖然實(shí)現(xiàn)簡(jiǎn)單，但去噪效果不理想.SURE-LET算法(Luisier et al.，2007)的先驗(yàn)?zāi)Ｐ驮诔叨葍?nèi)較好地模擬小波系數(shù)的概率分布，但忽略了小波系數(shù)尺度間的相關(guān)性，且其先驗(yàn)統(tǒng)計(jì)模型相對(duì)簡(jiǎn)單，不夠準(zhǔn)確.基于貝葉斯估計(jì)的小波域雙變量法(Sendur and Selesnick，2002)利用相鄰尺度間小波系數(shù)的相關(guān)性，取得較好的降噪效果.隨后此方法被應(yīng)用于地震信號(hào)降噪(劉鑫等,2006)，實(shí)驗(yàn)表明該方法具有良好的降噪性能.殷明等(2014)提出非下采樣雙樹(shù)復(fù)小波域的雙變量模型去噪算法，把雙變量統(tǒng)計(jì)模型引入到非下采樣雙樹(shù)復(fù)小波變換同一方向?qū)嵅亢吞摬啃〔ㄏ禂?shù)中，較為有效地去除圖像中的隨機(jī)噪聲.

Kingsbury(1998)提出了雙樹(shù)復(fù)小波變換(Dual-tree Complex Wavelet Transform，DTCWT),對(duì)傳統(tǒng)小波進(jìn)行了增強(qiáng)，使其具有近似平移不變、多方向選擇和低冗余的優(yōu)良特性，因而被廣泛應(yīng)用于信號(hào)降噪、特征提取和分析等領(lǐng)域(Kingsbury，2001).由于地震信號(hào)具有明顯的多尺度性，相鄰道共深度點(diǎn)有效信號(hào)有較強(qiáng)相關(guān)性，利用雙樹(shù)復(fù)小波變換能較好地檢測(cè)與表示地震信號(hào)的線狀變化特征(如彎曲的同相軸等).雙樹(shù)復(fù)小波變換的提出解決了二維可分離小波不能很好地表示信號(hào)細(xì)節(jié)特征和方向選擇少的不足，且系數(shù)的模具有旋轉(zhuǎn)不變性.這些優(yōu)點(diǎn)使它能更好地表示地震信號(hào)的高維奇異特征.

鑒于此，本文分析并驗(yàn)證了地震信號(hào)的雙樹(shù)復(fù)小波變換系數(shù)中同一方向?qū)嵅颗c虛部小波系數(shù)之間，以及實(shí)部(或虛部)系數(shù)與對(duì)應(yīng)的模之間具有較強(qiáng)的相關(guān)性.從而對(duì)同一方向?qū)嵅颗c虛部系數(shù)，以及實(shí)部(或虛部)系數(shù)與對(duì)應(yīng)的模分別建立雙變量模型.以此從含噪地震信號(hào)小波系數(shù)中估計(jì)恢復(fù)原信號(hào)的小波系數(shù)，得到兩種地震信號(hào)隨機(jī)噪聲壓制的雙樹(shù)復(fù)小波域雙變量方法.通過(guò)理論合成記錄仿真實(shí)驗(yàn)和實(shí)際地震資料處理結(jié)果證實(shí)本文兩種方法都能夠有效地壓制地震信號(hào)中的隨機(jī)噪聲.

圖1　一維雙樹(shù)復(fù)小波分解變換示意圖Fig.1　1D dual-tree complex wavelet transform

2　雙樹(shù)復(fù)小波變換

雙樹(shù)復(fù)小波分解變換包含兩個(gè)平行的分解樹(shù)，需要用到兩個(gè)獨(dú)立的實(shí)小波來(lái)分別計(jì)算實(shí)部與虛部系數(shù).實(shí)際中采用雙樹(shù)算法來(lái)實(shí)現(xiàn)分解變換.

圖1展示了一維信號(hào)的雙樹(shù)復(fù)小波分解變換，其中h0(n)和h1(n)分別是樹(shù)A分解過(guò)程的低通和高通濾波器，它產(chǎn)生實(shí)部的低頻和高頻小波系數(shù)；g0(n)和g1(n)分別是樹(shù)B分解過(guò)程的低通和高通濾波器，它產(chǎn)生虛部的低頻和高頻小波系數(shù)；↓2表示下采樣操作.

雙樹(shù)復(fù)小波所采用的實(shí)部濾波器組函數(shù)ψh(t)與虛部濾波器組函數(shù)ψg(t)構(gòu)成一個(gè)復(fù)小波，即

(1)

雙樹(shù)復(fù)小波變換的優(yōu)良性能來(lái)源于ψC(t)的解析性，為使ψC(t)滿足解析性，則需要ψh(t)與ψg(t)構(gòu)成一組希爾伯特變換對(duì)，即

(2)

使兩個(gè)正交小波函數(shù)組成希爾伯特變換對(duì)的充要條件是：兩低通濾波器滿足半采樣延遲條件(Coifman and Donoho，1995).則對(duì)正交小波基，有

(3)

(4)

(5)

式(5)與式(2)等價(jià)，稱式(5)為半幀移條件.由于半幀移相當(dāng)于在每個(gè)尺度對(duì)低通信號(hào)的采樣提高到原來(lái)的2倍，這就極大地避免了由于二元下抽樣導(dǎo)致的走樣現(xiàn)象產(chǎn)生.雙樹(shù)結(jié)構(gòu)和特殊濾波器的選擇使變換具有近似平移不變性和頻率無(wú)偏性.

二維雙樹(shù)復(fù)小波變換可通過(guò)一維雙樹(shù)復(fù)小波的張量積得到，分解時(shí)相當(dāng)于對(duì)行和列分別進(jìn)行一維雙樹(shù)復(fù)小波變換.因?yàn)殡p樹(shù)復(fù)小波變換是冗余的，對(duì)于一維信號(hào)，產(chǎn)生2∶1的冗余度，而對(duì)于n維的信號(hào)，將產(chǎn)生2n∶1的冗余度.二維雙樹(shù)復(fù)小波變換過(guò)程中，各尺度生成2個(gè)低頻子帶和方向分別為±15°、±45°和±75°的6個(gè)不同高頻子帶，且每個(gè)方向子帶有實(shí)部與虛部?jī)蓚€(gè)部分.對(duì)上一尺度的低頻子帶再進(jìn)行分解操作可達(dá)到多尺度分解的目的，得到不同尺度不同方向的小波分解系數(shù).

雙樹(shù)復(fù)小波變換后同一尺度內(nèi)的各子帶維數(shù)相同，同一位置的小波系數(shù)在同一尺度內(nèi)的各子帶位置固定，這樣有利于分析同一位置小波系數(shù)間的統(tǒng)計(jì)特性.用合適的分布擬合小波系數(shù)分布，并通過(guò)分析系數(shù)間的相關(guān)性特征，建立合適的統(tǒng)計(jì)模型，可提高信號(hào)降噪方法的能力(Hill et al.，2012).

3　雙樹(shù)復(fù)小波子帶系數(shù)相關(guān)性分析

基于小波域的降噪算法中，如何利用小波系數(shù)的統(tǒng)計(jì)分布特征一直是地震信號(hào)降噪算法領(lǐng)域的一個(gè)研究熱點(diǎn).小波系數(shù)具有零均值和重尾性，它的邊緣分布近似服從拉普拉斯分布、廣義高斯分布或混合高斯分布等非高斯統(tǒng)計(jì)分布.但是僅僅利用這些邊緣分布特征描述小波系數(shù)特征是不夠準(zhǔn)確的，因?yàn)檫@類邊緣統(tǒng)計(jì)分布特征忽略了小波系數(shù)之間的相關(guān)性.因此，如何建立較為準(zhǔn)確的數(shù)學(xué)模型來(lái)描述系數(shù)之間的相關(guān)關(guān)系，從而提高降噪效果是研究的難點(diǎn)(Yin and Liu，2012).

小波系數(shù)具有傳遞和聚集的特性，即幅值較大系數(shù)的鄰域系數(shù)幅值也較大，且在其他尺度同一位置系數(shù)的幅值也較大(Shi and Selesnick，2007).Sendur和Selesnick(2002)提出的雙變量閾值降噪法考慮了小波系數(shù)間的相關(guān)性，把小波系數(shù)的統(tǒng)計(jì)特征引入模型以提高降噪效果.由于雙樹(shù)復(fù)小波的實(shí)部和虛部濾波器構(gòu)成希爾伯特變換對(duì)，因此信號(hào)經(jīng)雙樹(shù)復(fù)小波變換后，同一方向?qū)嵅颗c虛部系數(shù)具有相關(guān)性.另外雙樹(shù)復(fù)小波系數(shù)的模具有旋轉(zhuǎn)不變性，為了將這一性質(zhì)應(yīng)用于雙變量模型，本文還考慮了實(shí)部(或虛部)系數(shù)與模的相關(guān)關(guān)系.為了衡量實(shí)部與虛部系數(shù)、實(shí)部(或虛部)系數(shù)和模之間的相關(guān)性，本文采用互信息(Liu and Moulin，2001)作為度量標(biāo)準(zhǔn).對(duì)于兩個(gè)變量z1和z2之間的互信息I(z1;z2)可表示為

(6)

其中p(z1)和p(z2)為z1和z2的概率密度函數(shù)，p(z1,z2)為z1和z2的聯(lián)合概率密度函數(shù).并利用Peng等(2005)所述方法計(jì)算互信息值I(z1;z2).互信息值越大，表示變量之間的相關(guān)性越強(qiáng).為了便于實(shí)驗(yàn)，對(duì)采樣點(diǎn)數(shù)為512，道數(shù)為512的4個(gè)合成地震記錄進(jìn)行2層雙樹(shù)復(fù)小波變換，對(duì)分解得到的子帶系數(shù)進(jìn)行互信息計(jì)算，分別計(jì)算父子系數(shù)、相鄰系數(shù)、實(shí)部與虛部系數(shù)、實(shí)部(或虛部)系數(shù)與模之間的互信息值，列出部分子帶間互信息值如下表1所示.

表1　雙樹(shù)復(fù)小波變換子帶系數(shù)的互信息

在表1中r、Pr和Nr分別表示實(shí)部系數(shù)、其對(duì)應(yīng)的父系數(shù)和鄰域系數(shù)；i、Pi和Ni分別表示與r同一方向虛部系數(shù)、其對(duì)應(yīng)的父系數(shù)和鄰域系數(shù)，m是實(shí)部與虛部系數(shù)對(duì)應(yīng)的模，其中

(7)

當(dāng)式(6)中的變量z1和z2分別為r、Pr、Nr、i、Pi、Ni、m中的兩者時(shí)，I(z1;z2)表示這兩個(gè)子帶系數(shù)之間的互信息值.

從表1中看出，同一方向的實(shí)部與虛部系數(shù)、實(shí)部(或虛部)系數(shù)與模存在著較強(qiáng)的相關(guān)性.相鄰尺度間系數(shù)即某一子帶系數(shù)與其父系數(shù)子帶大小不同，系數(shù)位置不固定，子帶系數(shù)與其鄰域系數(shù)雖子帶大小相同，在互信息度量時(shí)，考慮了取其鄰域的平均值，這些都影響了相關(guān)性的計(jì)算與運(yùn)用.由雙樹(shù)復(fù)小波分解變換特點(diǎn)知，同一方向的實(shí)部與虛部系數(shù)、實(shí)部(或虛部)系數(shù)與模子帶大小一致，系數(shù)位置固定，便于度量與應(yīng)用.所以本文考慮運(yùn)用實(shí)部與虛部系數(shù)、實(shí)部(或虛部)系數(shù)與模之間的相關(guān)性建立統(tǒng)計(jì)模型.

4　雙樹(shù)復(fù)小波域雙變量模型

4.1實(shí)部與虛部系數(shù)的雙變量模型

(8)

在給定含噪地震信號(hào)s的前提下，地震信號(hào)降噪的目的是根據(jù)一定的準(zhǔn)則盡可能地恢復(fù)原始地震信號(hào)x.在雙樹(shù)復(fù)小波域中，含噪模型可以表示為

(9)

其中，y為含噪小波系數(shù)，w為原始地震信號(hào)小波系數(shù)，n為噪聲的小波系數(shù).

考慮同一方向?qū)嵅颗c虛部小波系數(shù)之間相關(guān)性，構(gòu)造雙變量模型.假設(shè)

(10)

其中y1和y2是含噪地震信號(hào)同一方向的實(shí)部和虛部小波系數(shù)，w1和w2是對(duì)應(yīng)的原信號(hào)小波系數(shù)的實(shí)部和虛部，n1和n2為對(duì)應(yīng)的噪聲小波系數(shù)的實(shí)部和虛部.式(10)可以改寫(xiě)成向量的形式：

(11)

其中w=(w1,w2)T，y=(y1,y2)T，n=(n1,n2)T.

(12)

利用條件概率公式進(jìn)行推導(dǎo),式(12)可變?yōu)?/p>

(13)

為了計(jì)算方便，進(jìn)行對(duì)數(shù)變換后式(13)等價(jià)于

(14)

假設(shè)噪聲是獨(dú)立同分布的高斯白噪聲，這種噪聲在經(jīng)過(guò)雙樹(shù)復(fù)小波變換后仍服從高斯分布，即概率密度函數(shù)為

(15)

設(shè)原地震信號(hào)同一方向?qū)嵅颗c虛部小波系數(shù)的雙變量聯(lián)合概率密度函數(shù)為

(16)

其中a是待定參數(shù)，σ2表示被估計(jì)信號(hào)小波系數(shù)的方差.

(17)

(18)

4.2實(shí)部(或虛部)系數(shù)與模的雙變量模型

基于文中第3節(jié)的相關(guān)性分析，實(shí)部(或虛部)系數(shù)與模之間具有較強(qiáng)的相關(guān)性.故將實(shí)部與虛部系數(shù)的雙變量模型推廣建立實(shí)部(或虛部)系數(shù)與對(duì)應(yīng)模的雙變量模型.假設(shè)：

(19)

假設(shè)原地震信號(hào)實(shí)部系數(shù)與模的雙變量聯(lián)合概率密度函數(shù)為

(20)

其中t是待定參數(shù)，σ2表示被估計(jì)信號(hào)小波系數(shù)的方差.

(21)

(22)

5　雙樹(shù)復(fù)小波域雙變量模型降噪方法

5.1參數(shù)估計(jì)

(23)

其中yi為高頻子帶系數(shù).

(24)

(25)

其中N是鄰域U的大小.利用式(24)和(25)可得σ2的估計(jì)值:

(26)

對(duì)地震數(shù)據(jù)降噪效果的評(píng)價(jià)可以選用信噪比(SNR)進(jìn)行定量評(píng)價(jià).SNR的公式定義如下：

(27)

為了研究參數(shù)a和t對(duì)降噪效果的影響，本文選取5個(gè)合成地震信號(hào)數(shù)據(jù)，并對(duì)其分別加入標(biāo)準(zhǔn)差為10、20、30的高斯白噪聲進(jìn)行實(shí)驗(yàn).下面的圖2展示了當(dāng)參數(shù)a的取值范圍為[0,5]時(shí)采用雙樹(shù)復(fù)小波域?qū)嵅颗c虛部系數(shù)組成的雙變量模型降噪方法對(duì)SNR的影響，以及參數(shù)t取值范圍為[0,8]時(shí)實(shí)部(或虛部)系數(shù)與模組成的雙變量模型降噪方法對(duì)SNR的影響.圖2(a、b、c)中的5條曲線分別表示了5個(gè)不同地震數(shù)據(jù)降噪后SNR隨著參數(shù)a取值變化而變化.觀察得知在一定范圍內(nèi)，SNR值先是隨著參數(shù)a取值的增大而逐漸上升，在[1.5，3]間達(dá)到最大，即在此處取得較好的降噪效果.當(dāng)參數(shù)a取值繼續(xù)增大時(shí)SNR平緩地減小，對(duì)于加入了不同大小方差噪聲的合成地震數(shù)據(jù)皆有這一特征.圖2(d、e、f)中的5條曲線分別表示了5個(gè)不同地震數(shù)據(jù)降噪后SNR隨著參數(shù)t取值變化而變化，得到參數(shù)t最優(yōu)取值范圍為[2,4].因此本文取參數(shù)a經(jīng)驗(yàn)值為3，參數(shù)t經(jīng)驗(yàn)值為4.

圖2　參數(shù)a及t對(duì)SNR的影響(a)、(b)和(c)噪聲標(biāo)準(zhǔn)差分別為10，20，30時(shí)參數(shù)a對(duì)SNR的影響；(d)、(e)和(f)噪聲標(biāo)準(zhǔn)差分別為10，20，30時(shí)參數(shù)t對(duì)SNR的影響.Fig.2　The influence of parameter a and t on SNR(a)、(b)and (c)The influence of parameter a on SNR with standard deviation equals to 10,20,30；(d)、(e)and (f)The influence of parameter t on SNR with standard deviation equals to 10,20,30.

5.2降噪方法的算法步驟

本文中的雙樹(shù)復(fù)小波域?qū)嵅颗c虛部系數(shù)組成的雙變量降噪方法的算法步驟如下：

步驟1：對(duì)含噪地震信號(hào)進(jìn)行雙樹(shù)復(fù)小波分解變換，得到每一尺度2個(gè)低頻子帶和6個(gè)不同方向的高頻子帶，且每個(gè)子帶有實(shí)部與虛部?jī)蓚€(gè)部分；

步驟3：低頻系數(shù)保持不變,而對(duì)高頻同一方向?qū)嵅颗c虛部系數(shù)，運(yùn)用式(17)，(18)估計(jì)出不含噪地震信號(hào)的實(shí)部與虛部高頻小波系數(shù).從而得到各尺度各方向的小波分解系數(shù)的估計(jì)值；

步驟4：對(duì)小波分解系數(shù)的估計(jì)值進(jìn)行雙樹(shù)復(fù)小波逆變換重構(gòu)得到降噪后的地震信號(hào).

本文提出的雙樹(shù)復(fù)小波域?qū)嵅?或虛部)系數(shù)與模組成的雙變量降噪方法的算法步驟與上面步驟類似，不同之處在于步驟3利用式(21)和(22)估計(jì)出不含噪地震信號(hào)高頻小波系數(shù)的實(shí)部與虛部.

6　仿真實(shí)驗(yàn)與分析

為了驗(yàn)證本文降噪方法壓制地震信號(hào)隨機(jī)噪聲的有效性，利用MATLAB2014a進(jìn)行算法的仿真實(shí)驗(yàn).首先生成合成地震信號(hào)，并對(duì)其加入標(biāo)準(zhǔn)差為30的高斯白噪聲，對(duì)添加噪聲后的地震信號(hào)進(jìn)行一層二維雙樹(shù)復(fù)小波分解變換，所得到的小波系數(shù)如圖3所示.

圖3　含噪合成地震信號(hào)及其一層雙樹(shù)復(fù)小波分解系數(shù)(a)含噪地震信號(hào)；(b)低頻系數(shù)；(c)6個(gè)角度方向的高頻系數(shù).Fig.3　The noisy synthetic seismic signal and the one layer dual tree complex wavelet transform coefficients(a)The noisy seismic signal；(b)The low frequency coefficients；(c)The high frequency coefficients of the six angles.

圖3顯示了含噪合成地震信號(hào)經(jīng)過(guò)一層雙樹(shù)復(fù)小波分解后的低頻與高頻小波系數(shù).觀察圖3(b，c)發(fā)現(xiàn)隨機(jī)噪聲主要集中于高頻系數(shù)中.直觀上可看出同一方向?qū)嵅颗c虛部系數(shù)有較為明顯的相關(guān)性.因此本文方法只對(duì)高頻系數(shù)進(jìn)行處理，而保持低頻系數(shù)不變.最后對(duì)壓制噪聲后的小波系數(shù)進(jìn)行雙樹(shù)復(fù)小波逆變換，得到壓制噪聲后的地震信號(hào).

為了進(jìn)一步說(shuō)明本文方法的有效性，對(duì)合成地震數(shù)據(jù)加入不同大小方差的高斯白噪聲，將本文方法與基于父子系數(shù)的離散小波域雙變量閾值(DWT Bishrink)降噪方法(Abdelnour and Selesnick，2001)和基于父子系數(shù)的雙樹(shù)復(fù)小波域雙變量閾值(DTCWT Bishrink)降噪方法(Sendur and Selesnick，2002，劉鑫等,2006)進(jìn)行比較，其中各小波分解層數(shù)均為4層.本文中雙樹(shù)復(fù)小波域?qū)嵅颗c虛部系數(shù)組成的雙變量降噪方法記為方法1，雙樹(shù)復(fù)小波域?qū)嵅?或虛部)系數(shù)與模組成的雙變量降噪方法記為方法2.以SNR與均方誤差(MSE)作為評(píng)價(jià)標(biāo)準(zhǔn)，MSE的公式定義如下：

(28)

從表2中可以看出，本文兩種方法得到的SNR值明顯大于其他兩種方法得到的值，MSE值明顯小于其他兩種方法得到的數(shù)值.因此本文兩種方法的降噪效果更好，且多數(shù)情況下方法2的結(jié)果稍優(yōu)于方法1的結(jié)果.

為了能夠更好地體現(xiàn)實(shí)驗(yàn)效果，在仿真實(shí)驗(yàn)中取加噪地震信號(hào)其中的某一道地震數(shù)據(jù)進(jìn)行觀察.選擇降噪效果較好的DTCWT Bishrink方法與本文方法2進(jìn)行降噪效果比較.表2中降噪前信噪比為-4.50dB地震信號(hào)的第11道地震數(shù)據(jù)的降噪結(jié)果局部放大比對(duì)圖及對(duì)應(yīng)的振幅譜比較圖如圖4所示.

圖4　單道地震信號(hào)降噪結(jié)果及對(duì)應(yīng)的振幅譜(a)原始地震信號(hào)、加噪地震信號(hào)、DTCWT Bishrink方法和方法2降噪結(jié)果局部放大比較圖；(b)地震信號(hào)及其降噪后信號(hào)對(duì)應(yīng)的振幅譜；(c)振幅譜的局部放大圖.Fig.4　The denoised results of the single channel seismic signal and their corresponding amplitude spectrums(a)The local enlarged figure of the ideal seismic signal、the noisy seismic signal、the denoised seismic signals via the DTCWT Bishrink and the method 2；(b)The amplitude spectrums of the seismic signal and the denoised seismic signals via the DTCWT Bishrink and the method 2；(c)The local enlarged figure of the amplitude spectrums.

降噪前DWTBishrinkDTCWTBishrink方法1方法2SNRSNRMSESNRMSESNRMSESNRMSE7.5719.0428.539820.1022.341021.8514.942421.9214.7066-0.3913.22108.940014.0390.290015.9358.229516.0057.3045-4.5010.09223.681510.80189.973712.67123.530912.75121.3709

圖5　實(shí)際地震信號(hào)降噪結(jié)果(a)實(shí)際地震信號(hào)；(b)本文方法1降噪后的地震信號(hào)；(c)本文方法2降噪后的地震信號(hào);(d)DWT Bishrink降噪后的地震信號(hào)；(e)DTCWT Bishrink降噪后的地震信號(hào).Fig.5　The denoised results of the field seismic signal(a)The field seismic signal；(b)The denoised seismic signal via the method 1；(c)The denoised seismic signal via the method 2；(d)The denoised seismic signal via the DWT Bishrink；(e)The denoised seismic signal via the DTCWT Bishrink.

由圖4中的(a)—(c)子圖可以看出，本文方法能夠較好地壓制隨機(jī)噪聲.從圖4(b、c)振幅譜圖中可看出，高頻噪聲得到了有效壓制，且低頻中混有的噪聲也有所減弱.

綜合表2和圖4可以看出，本文方法在信噪比上較一些經(jīng)典的降噪方法有一定的提高，而且降噪得到的地震信號(hào)與原地震信號(hào)具有更小的差異.壓制隨機(jī)噪聲的同時(shí)較好地保留了地震信號(hào)的大部分細(xì)節(jié)信息.

7　實(shí)際地震資料處理

為了進(jìn)一步測(cè)試本文方法對(duì)地震信號(hào)降噪的效果，采用DWT Bishrink方法、DTCWT Bishrink方法、本文方法1和本文方法2對(duì)實(shí)際采集的地震信號(hào)進(jìn)行降噪處理，各小波分解層數(shù)均為4層.圖5顯示了實(shí)際地震信號(hào)降噪效果.

觀察圖5可知，本文提出的兩種方法在壓制隨機(jī)噪聲的同時(shí)也保留了更多的有效信號(hào)，具有較高的保真度.

8　結(jié)論

本文分析了地震信號(hào)經(jīng)雙樹(shù)復(fù)小波變換后同一方向?qū)嵅颗c虛部小波系數(shù)、實(shí)部(或虛部)系數(shù)與模的相關(guān)性.通過(guò)計(jì)算互信息進(jìn)行驗(yàn)證，證實(shí)了實(shí)部與虛部系數(shù)、實(shí)部(或虛部)系數(shù)與模之間存在較強(qiáng)的相關(guān)關(guān)系.因此對(duì)小波系數(shù)實(shí)部與虛部、實(shí)部(或虛部)系數(shù)與模分別建立雙變量模型，得到了兩種地震信號(hào)隨機(jī)噪聲壓制的雙變量方法.本文方法不僅在信噪比上較一些經(jīng)典的降噪算法有一定的提高，而且較好地保留了地震信號(hào)中的有效信息.當(dāng)然本文方法還有需要進(jìn)一步改進(jìn)的地方，例如小波系數(shù)之間的相關(guān)性可以考慮采用其他統(tǒng)計(jì)模型進(jìn)行分析，以獲得更好的降噪效果.

致謝作者非常感謝兩位匿名評(píng)審專家提出的寶貴意見(jiàn).

References

Abdelnour A F,Selesnick I W.2001.Nearly symmetric orthogonal wavelet bases.//Proceedings of the IEEE International Conference on Acoustics,Speech,Signal Processing (ICASSP).IEEE.Berkner K,Wells R O Jr.2002.Smoothness estimates for soft-threshold denoising via translation-invariant wavelet transforms.Applied and Computational Harmonic Analysis,12(1):1-24,doi:10.1006/acha.2001.0366.

Coifman P R,Donoho D L.1995.Translation-invariant de-noising.//Wavelets and Statistics.New York:Springer,125-150.

Donoho D L,Johnstone J M.1994.Ideal spatial adaptation by wavelet shrinkage.Biometrika,81(3):425-455.

Hill P,Achim A,Bull D.2012.The Undecimated dual Tree Complex Wavelet Transform and its application to bivariate image denoising using a Cauchy model.//2012 19th IEEE International Conference on Image Processing (ICIP).Orlando,FL:IEEE,1205-1208.

Kingsbury N.1998.The dual-tree complex wavelet transform:A new technique for shift invariance and directional filters.//Proceedings of European signal processing Conference.Rhodes,319-322.Kingsbury N.2001.Complex wavelets for shift invariant analysis and filtering of signals.Applied and Computational Harmonic Analysis,10(3):234-253,doi:10.1006/acha.2000.0343.

Liu J,Moulin P.2001.Information-theoretic analysis of interscale and intrascale dependencies between image wavelet coefficients.IEEE Transactions on Image Processing,10(11):1647-1658,doi:10.1109/83.967393.

Liu X,He Z H,Huang D J.2006.Seismic data denoising research based on wavelet transform.Petroleum Geophysics (in Chinese),4(4):15-18.Luisier F,Blu T,Unser M.2007.A new SURE approach to image denoising:Interscale orthonormal wavelet thresholding.IEEE Transactions on Image Processing,16(3):593-606,doi:10.1109/TIP.2007.891064.

Peng H C,Long F H,Ding C.2005.Feature selection based on mutual information criteria of max-dependency,max-relevance,and min-redundancy.IEEE Transactions on Pattern Analysis and Machine Intelligence,27(8):1226-1238,doi:10.1109/TPAMI.2005.159.

Sendur L,Selesnick I W.2002.Bivariate shrinkage with local variance estimation.IEEE Signal Processing Letters,9(12):438-441,doi:10.1109/LSP.2002.806054

Shi F,Selesnick I W.2007.An elliptically contoured exponential mixture model for wavelet based image denoising.Applied and Computational Harmonic Analysis,23(1):131-151,doi:10.1016/j.acha.2007.01.007.

Wang Z G,Zhou X X,Li J.2002.Noise attenuation using minimum smooth filter based on wavelet transform.Oil Geophysical Prospecting (in Chinese),37(6):594-600.

Yin M,Liu W.2012.Image denoising using mixed statistical model in nonsubsampled Contourlet transform domain.Acta Photonica Sinica,41(6):751-756,doi:10.3788/gzxb20124106.0751.Yin M,Bai R F,Xing Y,et al.2014.Denoising algorithm by Nonsubsampled Dual-tree Complex Wavelet domain bivariate model.Acta Photonica Sinica (in Chinese),43(10):1010004,doi:10.3788/gzxb20144310.1010004.

Zhang B,Yin X Y,Wu G C,et al.1999.Random noise elimination using wavelet analysis and inversion.Oil Geophysical Prospecting (in Chinese),34(6):635-641.Zhang K,Liu G Z,Zou D W,et al.1996.Seismic data time-frequency Domain Denoising based on dyadic wavelet transform.Acta Geophysica Sinica (in Chinese),39(2):265-271.Zong T,Meng H Y,Jia Y L,et al.1998.Denoising method based on wavelet packet and forward-backward prediction in F-X domain.Acta Geophysica Sinica (in Chinese),41(S1):337-346.

附中文參考文獻(xiàn)

劉鑫,賀振華,黃德濟(jì).2006.基于小波變換的地震資料去噪處理研究.油氣地球物理,4(4):15-18.

王振國(guó),周熙襄,李晶.2002.基于小波變換的最小光滑濾波去噪.石油地球物理勘探,37(6):594-600.

殷明,白瑞峰,邢燕等.2014.基于非下采樣雙樹(shù)復(fù)小波域的雙變量模型去噪算法.光子學(xué)報(bào),43(10):1010004,doi:10.3788/gzxb20144310.1010004.

張波,印興耀,吳國(guó)忱等.1999.用小波分析和反演方法聯(lián)合去除隨機(jī)噪聲的研究.石油地球物理勘探,34(6):635-641.

章珂,劉貴忠,鄒大文等.1996.二進(jìn)小波變換方法的地震信號(hào)分時(shí)分頻去噪處理.地球物理學(xué)報(bào),39(2):265-271.

宗濤,孟鴻鷹,賈玉蘭等.1998.小波包F-X域前后向預(yù)測(cè)去噪.地球物理學(xué)報(bào),41(S1):337-346.

(本文編輯胡素芳)

Dual-tree complex wavelet domain bivariate method for seismic signal random noise attenuation

WANG Jin-Ju,YUAN Li，LIU Wan-Ru,XU Xiao-Hong

School of Mathematics，Hefei University of Technology，Hefei 230009，China

Seismic signal noise attenuation is important in processing and interpreting seismic signal subsequently.Two dual-tree complex wavelet domain bivariate methods for seismic signal random noise attenuation are proposed.After the dual-tree complex discrete wavelet transform,the real and imaginary parts of the wavelet coefficients have dependency in the same direction.The real or imaginary parts and the corresponding magnitudes of the wavelet coefficients have dependency in the same direction.So we construct the bivariate model for the real and imaginary parts in the same direction.Using the model the wavelet coefficients of the original seismic signal are estimated.The denoised seismic signal is achieved based on the dual-tree complex wavelet inverse transform.The proposed algorithm is also extended to the real or imaginary parts and the corresponding magnitudes of dual-tree complex wavelet coefficients in the same direction.Through experiments on the synthetic seismic record and the field seismic data,the results demonstrate that the proposed two methods can attenuate random noise effectively.

Dual-tree complex wavelet transform;Random noise;Bivariate model;Wavelet coefficients

汪金菊，袁力，劉婉如等.2016.地震信號(hào)隨機(jī)噪聲壓制的雙樹(shù)復(fù)小波域雙變量方法.地球物理學(xué)報(bào)，59(8):3046-3055,

10.6038/cjg20160827.

Wang J J,Yuan L,Liu W R,et al.Dual-tree complex wavelet domain bivariate method for seismic signal random noise attenuation.Chinese J.Geophys.(in Chinese),59(8):3046-3055,doi:10.6038/cjg20160827.

國(guó)家重大科研裝備研制項(xiàng)目“深部資源探測(cè)核心裝備研發(fā)”(ZDYZ2012-1)—06子項(xiàng)目“金屬礦地震探測(cè)系統(tǒng)”—05課題“系統(tǒng)集成野外試驗(yàn)與處理軟件研發(fā)”，中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(J2014HGXJ0072,2015HGZX0018)，國(guó)家級(jí)大學(xué)生創(chuàng)新創(chuàng)業(yè)訓(xùn)練計(jì)劃項(xiàng)目(201410359075)共同資助.

汪金菊，女，1978年生，博士，講師，主要從事信號(hào)處理等工作.E-mail:wjj@hfut.edu.cn

10.6038/cjg20160827

P631

2015-12-01，2016-01-19收修定稿