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

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

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

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

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

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

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

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

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

1.1 基于PCA的噪聲估計

設含噪信號塊的生成模型為

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

(1)

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

(2)

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

(3)

(4)

1.2 低秩信號塊選擇

信號結構的復雜度可利用梯度協方差計算。信號塊yi的梯度矩陣為Gyi=[DhyiDvyi]，其中Dh與Dv為水平與垂直導數算子。梯度協方差矩陣為

(5)

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

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

(6)

則含噪塊yf的紋理強度為

(7)

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

(8)

式中：α為形狀參數；β為尺度參數。

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

(9)

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

1.3 噪聲估計迭代過程

選出低秩塊就可準確估計噪聲標準差。但選擇弱紋理塊的閾值τ應可變。用迭代來選擇需要的塊并估計噪聲。具體過程如下：

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

步驟4然后用Wk+1與τk+1估計噪聲標準差；

步驟5重復步驟2～步驟4，直至噪聲標準差不變。

2 波原子變換

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

2.1 波原子定義

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

(10)

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

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

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

(11)

且滿足頻域局部條件

(12)

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

2.2 波原子構造及系數求解

(13)

(14)

離散波原子變換的系數為

coeffj,m,n=

(15)

(16)

(17)

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

3.1 總體去噪流程

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

步驟1利用噪聲估計算法估計信號噪聲標準差；

步驟2對地震信號進行循環平移運算；

步驟3循環平移后信號做波原子變換得到系數；

步驟4對變換系數進行閾值選取和修正；

步驟5構造閾值函數，基于該函數對波原子變換系數進行篩選處理；

步驟6對處理后系數進行波原子反變換；

步驟7對反變換結果進行反循環平移，得到當前循環平移的結果；

步驟8重復步驟2～步驟7，至所有循環平移完成；

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

3.2 地震信號循環平移

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

(18)

3.3 閾值選取及修正

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

(19)

式中：K為可調參數；σn為噪聲標準差；N為波原子系數長度，但該閾值為全局閾值。通常噪聲的波原子系數在低尺度時較大，在高尺度時較小，因此閾值參數也應該隨尺度的變化而變化，即閾值應與波原子變換尺度成反比關系，但變化范圍不宜過大。且全局閾值會導致地震信號細節的丟失，因此本文采用分層閾值

(20)

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

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

(21)

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

3.4 閾值函數構造

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

(22)

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

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

(23)

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

(24)

式中：T介于0～1。盡管效果有所改善，但T是固定值，針對不同尺度的系數，很多學者提出了改進算法[20]。考慮到軟硬閾值的優缺點，要克服硬閾值函數在閾值點處不連續的缺點，還有解決軟閾值函數處理前后存在的固有偏差。因此需找到一種函數，在閾值點處連續，且應為高階可導的。綜合軟硬閾值不同特點，本文構造了一種新的閾值函數

(25)

4 實驗結果與分析

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

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

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

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

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

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

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

4.2 合成地震信號去噪

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

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

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

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

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

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

4.3 實際地震信號去噪

本實驗信號源自北海F3數據體，F3是北海位于荷蘭的一個區塊，已采用OpendTect進行了傾角導向濾波。信號范圍為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 實際地震信號去噪評價指標Tab.3 Denoising evaluation index of real seismic signal

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

5 結 論

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

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