何 峰，周亞同，趙翔宇，石超君

(河北工業(yè)大學(xué) 電子信息工程學(xué)院 天津市電子材料與器件重點實驗室，天津 300401)

地震信號由震源的振動產(chǎn)生，沿地層傳播后被檢波器采集接收[1]。在采集過程中地震信號往往含有高斯隨機(jī)噪聲。去噪時噪聲水平未知，只能進(jìn)行盲去噪。目前信號盲去噪算法主要分三類：①噪聲估計與非盲去噪結(jié)合的去噪[2]；②信號評估與無參非盲去噪[3]；③同時進(jìn)行噪聲估計并去噪[4]。

本文擬采用第一類盲去噪算法，首先對信號進(jìn)行噪聲估計，然后利用噪聲標(biāo)準(zhǔn)差作為非盲去噪輸入?yún)?shù)。目前噪聲估計主要為三類：①基于塊算法[5]，將信號規(guī)則或非規(guī)則分塊，再找出強(qiáng)度變換(噪聲引起)最小塊，此塊即標(biāo)準(zhǔn)差最小的塊；②基于高通濾波算法[6]，將信號經(jīng)過高通濾波，再計算含噪信號與濾波后信號的差值，得到噪聲標(biāo)準(zhǔn)差。但對復(fù)雜紋理信號，此算法有時不準(zhǔn)確，會發(fā)生過估計；③基于統(tǒng)計算法：Zoran等[7]提出一種基于統(tǒng)計的估計算法，分析DCT(Discrete Cosine Transform)濾波的信號后，發(fā)現(xiàn)噪聲的存在導(dǎo)致了峰值變化，利用該影響準(zhǔn)確估計出噪聲標(biāo)準(zhǔn)差。本文所用估計算法為基于塊的算法[8]。

常見變換域非盲去噪算法主要有：小波去噪[9]、脊波去噪[10]、輪廓波去噪[11]、曲波去噪[12]等。波原子變換是一種新型的多尺度幾何分析工具，相比小波、曲波等，對震蕩型紋理型信號具有“最優(yōu)”稀疏表示，且能同時表征沿振蕩與穿過該振蕩方向的信息。已有研究學(xué)者將其應(yīng)用于信號去噪領(lǐng)域，如Haddad等[13]將波原子應(yīng)用至典型紋理型信號(指紋)壓縮，表現(xiàn)出波原子優(yōu)勢。Zhang等[14]提出一種波原子結(jié)合循環(huán)平移的算法，較好地抑制了信號邊緣易出現(xiàn)的偽吉布斯現(xiàn)象。楊寧等[15]將波原子應(yīng)用至疊前地震信號的信噪分離，利用各尺度間系數(shù)相關(guān)性去噪。這些研究表明波原子對振蕩紋理型信號的處理優(yōu)勢。而地震同相軸的排列組合可視作紋理，故波原子變換亦適用于處理地震信號。

本文提出一種基于波原子變換的三維地震信號盲去噪算法，利用噪聲估計得到信號的噪聲標(biāo)準(zhǔn)差，再對信號進(jìn)行循環(huán)平移并在波原子變換域處理，按不同尺度分層選取閾值并進(jìn)行修正，然后利用改進(jìn)閾值函數(shù)處理波原子系數(shù)，再波原子反變換與反循環(huán)平移得到去噪后信號。合成與實際地震信號實驗表明，本文算法能夠有效去除地震信號噪聲。

1 基于塊的地震信號噪聲估計

針對紋理復(fù)雜場景易出現(xiàn)過估計問題，從含噪地震信號中選出不含高頻分量的低秩信號塊，再從利用主成分分析(Principal Component Analysis，PCA)估計噪聲。

1.1 基于PCA的噪聲估計

設(shè)含噪信號塊的生成模型為

yi=zi+ni,i=1,2,…,M

(1)

式中：M為分塊數(shù)量；zi為無噪信號塊；ni為均值為0，方差為σ2的高斯白噪聲。將數(shù)據(jù)方差投影到某個軸，則投影向量的方差為

(2)

式中：u為單位向量；V(uTyi)為zi在u方向的方差。可利用PCA求解，則

(3)

(4)

1.2 低秩信號塊選擇

信號結(jié)構(gòu)的復(fù)雜度可利用梯度協(xié)方差計算。信號塊yi的梯度矩陣為Gyi=[DhyiDvyi]，其中Dh與Dv為水平與垂直導(dǎo)數(shù)算子。梯度協(xié)方差矩陣為

(5)

式中：T為轉(zhuǎn)置運(yùn)算符；V為Cyi的特征向量；s1與s2為特征值。協(xié)方差矩陣的跡能夠反映出信號塊的紋理強(qiáng)度。則定義紋理強(qiáng)度ξi為：ξi=trace(Cyi)，其中trace()為求跡運(yùn)算符。ξi越小，則信號越平滑，紋理越弱。梯度矩陣對噪聲極敏感。此時需考慮一種低秩塊的極端情況：理想的無噪信號塊zf的梯度為0，即Gzf=[DhzfDvzf]=[0 0]，則得yf的梯度矩陣為

Gyf=[Dh(zf+n)Dv(zf+n)]=
[DhnDvn]=Gn

(6)

則含噪塊yf的紋理強(qiáng)度為

(7)

為簡化問題求解，用Gamma分布來近似ξn分布，其概率密度函數(shù)為

(8)

式中：α為形狀參數(shù)；β為尺度參數(shù)。

假設(shè)給出的含噪塊為理想塊(即無噪信號塊的梯度為0)，置信區(qū)間為：P(0<ξ(n)<τ)=1-δ，其中τ為閾值，δ為顯著性水平。若信號塊的ξi小于閾值則認(rèn)為其為弱紋理塊(低秩塊)。閾值τ可由δ與σn定義

(9)

式中：Gamma-1(δ,α,β)為Gamma反分布;N2為信號塊中的像素數(shù)量。

1.3 噪聲估計迭代過程

選出低秩塊就可準(zhǔn)確估計噪聲標(biāo)準(zhǔn)差。但選擇弱紋理塊的閾值τ應(yīng)可變。用迭代來選擇需要的塊并估計噪聲。具體過程如下：

步驟3用τk+1從含噪信號中選出弱紋理塊Wk+1；

步驟4然后用Wk+1與τk+1估計噪聲標(biāo)準(zhǔn)差；

步驟5重復(fù)步驟2～步驟4，直至噪聲標(biāo)準(zhǔn)差不變。

2 波原子變換

波原子是一種新型的多尺度幾何分析工具[16]，是一種遵循拋物比例尺度關(guān)系的二維小波包的變體。對于給定精度，以一個N×N大小的二維信號為例，僅需o(N)波原子系數(shù)就可稀疏表示該信號，若要表示相同精度，則需曲波系數(shù)o(N3/2)，小波和Gabor系數(shù)o(N2)。波原子能以最少的系數(shù)最優(yōu)稀疏表示信號，因此本文選擇波原子作為變換基。

2.1 波原子定義

若將波原子表示為wμ(f)，其中μ=(j,m,n)=(j,m1,m2,n1,n2),j,m1,m2,n1,n2∈Z，設(shè)(fμ,ωμ)為相空間中任意的一點，且滿足

(10)

式中：C1,C2為正數(shù)；j,m,n分別為尺度、頻率與空間系數(shù)；fμ與±ωμ為空間域wμ(f)與頻率域wμ(ω)中心。

若wμ(f)滿足空域局部條件

|wμ(f)|≤CM2j(1+2j(ω-ωμ))-M

(11)

且滿足頻域局部條件

(12)

則稱波包{wμ}為波原子。

2.2 波原子構(gòu)造及系數(shù)求解

(13)

(14)

離散波原子變換的系數(shù)為

coeffj,m,n=

(15)

(16)

(17)

3 基于波原子變換的三維地震信號盲去噪過程

3.1 總體去噪流程

本文基于波原子變換的三維地震信號盲去噪的總體去噪流程為：

步驟1利用噪聲估計算法估計信號噪聲標(biāo)準(zhǔn)差；

步驟2對地震信號進(jìn)行循環(huán)平移運(yùn)算；

步驟3循環(huán)平移后信號做波原子變換得到系數(shù)；

步驟4對變換系數(shù)進(jìn)行閾值選取和修正；

步驟5構(gòu)造閾值函數(shù)，基于該函數(shù)對波原子變換系數(shù)進(jìn)行篩選處理；

步驟6對處理后系數(shù)進(jìn)行波原子反變換；

步驟7對反變換結(jié)果進(jìn)行反循環(huán)平移，得到當(dāng)前循環(huán)平移的結(jié)果；

步驟8重復(fù)步驟2～步驟7，至所有循環(huán)平移完成；

步驟9對每次循環(huán)平移得到的去噪結(jié)果求平均值，得到最終的去噪地震信號。

3.2 地震信號循環(huán)平移

波原子變換缺乏平移不變性，會導(dǎo)致在邊緣或紋理等不連續(xù)點臨域內(nèi)產(chǎn)生偽吉布斯現(xiàn)象，造成去噪信號的失真[17]。為了防止該現(xiàn)象，Coifman等[18]提出一種循環(huán)平移算法，即將信號做循環(huán)平移，然后閾值去噪，再做反循環(huán)平移。每次偽吉布斯現(xiàn)象會在不同位置出現(xiàn)。三維地震信號循環(huán)平移表達(dá)式為

(18)

3.3 閾值選取及修正

記波原子變換系數(shù)的閾值為λ，作為波原子去噪算法中的重要參數(shù)，會直接影響去噪效果。Donoho等[19]提出一種VisuShrink閾值法

(19)

式中：K為可調(diào)參數(shù)；σn為噪聲標(biāo)準(zhǔn)差；N為波原子系數(shù)長度，但該閾值為全局閾值。通常噪聲的波原子系數(shù)在低尺度時較大，在高尺度時較小，因此閾值參數(shù)也應(yīng)該隨尺度的變化而變化，即閾值應(yīng)與波原子變換尺度成反比關(guān)系，但變化范圍不宜過大。且全局閾值會導(dǎo)致地震信號細(xì)節(jié)的丟失，因此本文采用分層閾值

(20)

式中：j為波原子變換的尺度；Nj為第j層波原子變換系數(shù)的長度?？紤]到地震信號的特點，波原子變換的系數(shù)偏大，因此閾值也應(yīng)相應(yīng)變大，但不宜過大，因此考慮引入指數(shù)函數(shù)，且指數(shù)參數(shù)不宜過大，應(yīng)為一接近零的值。經(jīng)由實驗發(fā)現(xiàn)，隨地震信號的含噪量的增大，閾值也應(yīng)隨之稍有增大。因此在多尺度分層閾值選取的基礎(chǔ)上，引入噪聲參數(shù)，提出了隨信號含噪量自適應(yīng)變化的閾值修正

λnew=λ·e1/(m·n·d)

(21)

式中：m為波原子分解最大尺度；n為一常數(shù)，本文取n=4；d與噪聲大小有關(guān)d=2/σn。

3.4 閾值函數(shù)構(gòu)造

硬閾值函數(shù)保留了信號的局部特性，但由于其不連續(xù)性，導(dǎo)致處理信號時存在一定的波動與振蕩。硬閾值函數(shù)為

(22)

式中：cλ為閾值函數(shù)處理后得到的系數(shù)；c為波原子系數(shù)；λ為閾值。

軟閾值函數(shù)在閾值處連續(xù)，處理系數(shù)時會更平滑，但會損失一部分高頻信息，且存在固有的偏差，會使信號模糊。軟閾值函數(shù)為

(23)

軟閾值函數(shù)的連續(xù)性比硬閾值好很多，但為克服其存在的固有偏差，又提出了軟硬閾值折衷算法

(24)

式中：T介于0～1。盡管效果有所改善，但T是固定值，針對不同尺度的系數(shù)，很多學(xué)者提出了改進(jìn)算法[20]?？紤]到軟硬閾值的優(yōu)缺點，要克服硬閾值函數(shù)在閾值點處不連續(xù)的缺點，還有解決軟閾值函數(shù)處理前后存在的固有偏差。因此需找到一種函數(shù)，在閾值點處連續(xù)，且應(yīng)為高階可導(dǎo)的。綜合軟硬閾值不同特點，本文構(gòu)造了一種新的閾值函數(shù)

(25)

4 實驗結(jié)果與分析

為驗證本文算法對三維地震信號的去噪效果，以合成地震信號為例，驗證閾值修正、改進(jìn)閾值函數(shù)與循環(huán)平移對傳統(tǒng)波原子去噪效果的提升。然后以合成與實際地震信號為例，將小波變換(Wavelet Transform，WT)、雙樹復(fù)小波變換(Dual Tree-Complex Wavelet Transform,CPLX_DT_WT)、曲波變換(Curvelet Transform,CT)以及傳統(tǒng)波原子變換(Wave Atoms Transform,WAT)與本文算法對比。采用輸入信噪比(SNR_IN)作為輸入地震信號的含噪大小衡量指標(biāo)，以輸出信噪比(SNR_OUT)、均方誤差(MSE)以及峰值信噪比(PSNR)作為去噪質(zhì)量評價指標(biāo)。

4.1 去噪處理過程前后性能對比

本實用采用合成地震信號，并歸一化至[0,1]。信號共41×41=1 681道，每道85個采樣點。為定量分析去噪效果，分別添加信號最大幅值的5%，10%，15%，20%，25%與30%的噪聲。圖1為歸一化后信號。

圖1 三維合成地震信號Fig.1 3D synthetic seismic signal

表1為各去噪手段處理前后的性能對比，其中“未修正”表示多尺度閾值與硬閾值函數(shù)；“閾值修正”表示分層多尺度閾值后進(jìn)行閾值修正，再硬閾值函數(shù)處理系數(shù)；“改進(jìn)閾值函數(shù)”表示分層多尺度閾值并進(jìn)行閾值修正，然后采用改進(jìn)閾值函數(shù)；“循環(huán)平移”表示在“改進(jìn)閾值函數(shù)”前對地震信號進(jìn)行循環(huán)平移。其中參數(shù)設(shè)置為p=0.9，q=10(σ≥25%)或15(σ<25%)。對比指標(biāo)采用PSNR。

表1 去噪處理過程性能對比Tab.1 Performance comparison of denoising process

從表1中列與列兩兩對比可知，閾值修正前后PSNR有了一定提升，噪聲較小時提升不明顯，但隨信號含噪量增大，去噪效果提升越明顯，表明閾值修正對去噪效果的提升；改進(jìn)閾值函數(shù)比硬閾值函數(shù)處理系數(shù)能獲得更優(yōu)秀去噪效果；循環(huán)平移處理前后去噪效果提升最明顯，即使在含噪量僅為5%時，也將PSNR從34.089 3提升至34.769 2。但循環(huán)平移會造成運(yùn)算時間增加。從整體分析可看出，數(shù)據(jù)從左至右，PSNR從小變大，且隨信號含噪量增大，PSNR值也越大。表明閾值修正、改進(jìn)閾值函數(shù)與循環(huán)平移處理的有效性。

4.2 合成地震信號去噪

本實驗針對合成三維地震信號去噪，所用地震信號同第“4.1”節(jié)，添加噪聲大小分別為信號最大幅值的5%，10%，15%，20%，25%及30%。對比算法采用全局經(jīng)驗閾值。因為經(jīng)驗閾值往往為能夠取得各算法最好去噪效果的閾值。圖2為對圖1合成地震信號添加噪聲后信號。

圖2 含噪的合成三維地震信號Fig.2 3D synthetic seismic signal with noise

表2 合成地震信號去噪評價指標(biāo)Tab.2 Denoising evaluation index of synthetic seismic signal

從表2可知，同一噪聲水平下，針對各重建算法，由于MSE較小，因此去噪效果差距不明顯，但從SNR_OUT與PSNR可以看出，WT去噪質(zhì)量稍差，而CPLX_DT_WT與CT，WAT的去噪效果優(yōu)于WT。本文算法在任何噪聲水平下均取得最佳的去噪效果，明顯優(yōu)于其他對比算法。圖3為添加噪聲幅度為5%時，各算法去噪效果圖。

圖3 各算法去噪結(jié)果Fig.3 Denoising results of each algorithm

從圖3去噪結(jié)果可知，WT與CPLX_DT_WT的去噪結(jié)果圖3(a)與圖3(b)表面仍稍有較明顯的噪聲痕跡；CT與WAT的去噪效果如圖3(c)與圖3(d)，只殘余輕微的噪聲痕跡；而本文提出的算法，圖3(e)可知，本文算法的去噪效果最好，基本將噪聲去除。本實驗說明，相對其他對比算法，本文算法能以較大優(yōu)勢去除合成地震信號的噪聲。

4.3 實際地震信號去噪

本實驗信號源自北海F3數(shù)據(jù)體，F(xiàn)3是北海位于荷蘭的一個區(qū)塊，已采用OpendTect進(jìn)行了傾角導(dǎo)向濾波。信號范圍為In-line：400-527；Cross-line：700-827；Time：480-732。In-line與Cross-line方向空間采樣率為1，時間采樣率為4 ms，故信號大小為128×128×64。將信號歸一化至[0，1]。圖4(a)、圖4(b)分別為原始實際三維地震信號的三維與橫切片圖。圖5為添加信號后的信號。

圖4 三維實際地震信號Fig.4 3D real seismic signal

圖5 含噪的實際地震信號Fig.5 Real seismic signal with noise

表3 實際地震信號去噪評價指標(biāo)Tab.3 Denoising evaluation index of real seismic signal

首先從圖6(a)和圖6(b)可知，WT與CPLX_DT_WT去噪算法的去噪效果較差，仍可在表面觀察到噪點的存在，未完全去除噪聲，這是由于WT算法運(yùn)算簡單，算法復(fù)雜度較低，對復(fù)雜的地震信號處理存在一定劣勢，但運(yùn)算速度較快。從圖6(c)、圖6(d)、圖6(e)可知，其余三種去噪算法較好地完成了去噪，表面幾乎觀察不到噪聲的存在。三維圖的表面以及去噪結(jié)果切片圖中都不再含有噪點，可觀察到地震信號清晰的紋理與結(jié)構(gòu)。其次從表5中也可以定理分析出，除本文算法之外的其它幾種對比算法，去噪效果最好的是WAT算法，其次是曲波算法、雙樹復(fù)小波算法以及小波算法，證明了波原子相對曲波、雙樹復(fù)小波以及小波對地震信號的處理優(yōu)勢。去噪效果最佳的是本文算法，證明了本文算法中用到的分層閾值選取及閾值修正、改進(jìn)閾值函數(shù)以及循環(huán)平移等處理手段的有效性，表明本文算法的準(zhǔn)確性。

5 結(jié) 論

(1)波原子作為一種新型多尺度幾何分析工具，能夠比小波、曲波、Gabor原子等工具更稀疏的表示紋理豐富的地震信號，本文提出的基于波原子變換的三維地震信號盲去噪算法，直接將含噪三維地震信號進(jìn)行去噪，且無需給出噪聲的標(biāo)準(zhǔn)差，提高了去噪效率。該算法的創(chuàng)新點主要為：①結(jié)合基于塊的噪聲估計算法對含噪地震信號進(jìn)行噪聲估計，用于盲去噪；②采用分層多尺度閾值，并進(jìn)行閾值修正；③對閾值函數(shù)進(jìn)行改進(jìn)，改善硬閾值、軟閾值的缺陷；④利用循環(huán)平移解決偽吉布斯現(xiàn)象，并提升地震信號的去噪效果。

(2)本文提出的閾值修正與改進(jìn)閾值處理函數(shù)，結(jié)合噪聲估計與循環(huán)平移算法，較好地完成了合成三維地震信號與實際三維地震信號的去噪，在SNR_OUT，MSE，MAE及PSNR等效果指標(biāo)上要優(yōu)于小波變換、雙樹復(fù)小波變換、曲波變換以及傳統(tǒng)波原子變換去噪算法，表明了本文算法的有效性與準(zhǔn)確性，為地震信號的變換域去噪提供了一種新思路。下一步工作將繼續(xù)研究閾值選取與修正、改進(jìn)閾值函數(shù)等，使其更加具有自適應(yīng)性。