Shearlet域自適應閾值地震數據隨機噪聲壓制

薛 林 程 浩* 鞏恩普 陳毅軍

(①東北大學深部金屬礦山安全開采教育部重點實驗室，遼寧沈陽 110004； ②東北大學資源與土木工程學院，遼寧沈陽 110004)

0 引言

地震勘探在地質調查和資源探測方面具有重要作用[1]。實際地震數據往往混雜隨機噪聲，嚴重影響數據的分析和處理效果。所以有效去除噪聲、獲得高信噪比的地震數據是地震勘探首先要解決的問題[2-4]。

隨機噪聲壓制主要是基于有效信號與隨機噪聲在頻率或視速度方面的差異性，常用的方法有中值濾波[5-7]、奇異值分解[8-10]、f-x域預測濾波[11-13]等。小波變換[14-16]作為一種多尺度時頻分析手段，能夠多尺度分析地震數據，但在處理多維地震數據時效果不明顯。而Curvelet變換[17-19]和Contourlet變換[20-22]能夠更好地描述地震信號的局部特征，具有突出的稀疏表示能力，可以實現很好的去噪效果。

Guo等[23-24]進一步研究了小波理論，結合多尺度幾何分析和膨脹仿射系統，提出了Shearlet變換理論。Shearlet變換繼承了小波變換的局部性和多尺度特性。由于Shearlet變換具有很強的方向性，而且其頻域支撐是各向異性的，因此Shearlet變換可以用最少的稀疏系數逼近奇異曲線，實現信號的最優稀疏表示。同時Shearlet變換可多尺度描述高維信號的幾何特性，且方向的數目隨著尺度的增加而增加，因而可以克服小波變換對多維地震數據處理效果不佳的缺陷。

由于Shearlet變換獨特的優勢，其在多個領域得到了廣泛的使用。鄭晶等[25]首先對地震數據進行奇異值分解(SVD)去除直達波，然后通過Shearlet變換改進閾值壓制隨機噪聲；李娟等[26]對地震數據進行Shearlet分解，通過增大信號分布集中的方向系數來增強信號與噪聲的差異，然后利用峰度特性拾取地震初至波；鄭升等[27]對大量含隨機噪聲的沙漠地震數據進行Shearlet變換，將得到的系數作為輸入標簽用于訓練深度殘差卷積神經網絡模型，得到預測噪聲和有效信號的變換系數；王憲楠等[28]利用探地雷達數據在Shearlet域的差異，采用硬閾值壓制探地雷達數據中的隨機噪聲；劉昕等[29]將非下采樣與Shearlet變換結合以壓制隨機噪聲，消除傳統Shearlet變換去噪中存在的吉布斯現象；吳一全等[30]基于非下采樣Shearlet變換對圖像的低頻和高頻分量分別進行奇異值分解和各向異性擴散，以達到去噪目的，很好地保留了圖像的細節信息。

傳統的Shearlet變換閾值噪聲壓制方法，主要考慮信號在Shearlet域尺度層面的特性，針對不同的尺度層設置自適應閾值，卻忽略了地震數據在方向上的變化特性。本文通過分析信號在Shearlet域方向上的分布特點，提出閾值隨尺度和方向同時自適應變化的噪聲壓制方法。分別對一組理論模型和兩組實際數據進行噪聲壓制試驗，并與傳統閾值去噪結果對比，驗證本文的尺度和方向自適應閾值法能有效壓制噪聲，最大程度上保留有效信號，提高地震數據的信噪比。

1 Shearlet變換基本理論

1.1 連續Shearlet變換

Guo等[31]基于合成小波理論，結合仿射系統和多尺度幾何分析，提出了帶有合成膨脹系統的幾何變換方法。合成連續膨脹系統在二維空間表達式為

ψAB(φ)=φj,l,k(x)=|detA|j/2×

φ(BlAjx-k∶j,l∈Z,k∈Z2)

(1)

若ψAB(φ)滿足Parseval框架[32-33]，則稱為合成小波。Aj控制尺度變換，Bl控制幾何變換。

定義伸縮矩陣H由各向異性膨脹矩陣A和剪切矩陣B構成

(2)

則可以得到如下連續仿射系統

ψH(φ)=ψa,b,k(φ)

a∈R+,b∈R,k∈R2}

(3)

式中：t為平移參數；ψa,b,k(φ)為連續剪切波；ψa,b,k(x)為剪切母函數。剪切波ψa,b,k(φ)在各個尺度層對應的頻域劃分如圖1所示。

基于剪切波函數，可以定義連續Shearlet變換為

C=SH{f}=〈f,φi,l,k〉

(4)

式中：C為Shearlet系數矩陣； SH{·}為Shearlet變換； 〈·，·〉表示內積；f為含噪信號。

同樣可以對Shearlet系數矩陣C進行反變換重構原始函數

f=SH-1{C}

(5)

式中SH-1表示Shearlet逆變換。

圖1 Shearlet頻域剖分示意圖

1.2 離散Shearlet變換

假設尺度參數a=2-j,剪切參數b=-l(j,l∈Z)，用點k∈Z2代替連續Shearlet變換中的平移參數t∈R2，得到Shearlet變換基元素的表達式為

(6)

對于任意f∈L2(R2)，離散Shearlet變換需滿足

(7)

2 Shearlet變換閾值去噪

2.1 Shearlet變換閾值去噪原理

基于Shearlet變換最優稀疏表示的特點，地震數據經Shearlet變換后得到各分解尺度和方向的Shearlet系數。若Shearlet基函數的方向與地震信號的方向大致相同，則對應較大的Shearlet系數；當兩方向相差較大時，Shearlet系數偏小。在實際地震數據中，有效信號分布具有方向性，而隨機噪聲是隨機分布的，所以地震數據經Shearlet變換后，有效信號相比隨機噪聲會對應較大的Shearlet系數。這樣就可以設置一個閾值，利用閾值函數去除隨機噪聲對應的較小的系數，達到去除隨機噪聲的目的。然后對處理后的Shearlet系數進行逆變換，得到去噪后的地震信號。

假設待處理的含噪地震數據模型為

f(t)=u(t)+n(t)

(8)

式中：u(t)為有效信號；n(t)為隨機噪聲。Shearlet變換閾值去噪算法的實質就是從f(t)中得到u(t)。該過程可以分以下三個步驟。

(1)對含噪地震數據f(t)進行Shearlet變換，獲取Shearlet系數

C(j,l,k)=SH{f(t)}=〈f(t),φj,l,k〉

(9)

(2)對得到的Shearlet系數進行閾值化

(10)

(11)

(12)

式中：Cnew(j,l,k)表示閾值化處理后的系數；N表示地震數據的總采樣點數； median(·)表示對矩陣所有元素取中值。

(3)對閾值化處理后的系數Cnew進行逆變換，得到去噪后的有效信號

(13)

2.2 自適應閾值函數

傳統閾值函數在壓制隨機噪聲時，僅對不同分解尺度設置不同的閾值，實現了閾值隨分解尺度的自適應。但因為Shealret變換具有多尺度、多方向的特點，所以僅隨尺度自適應閾值不能充分發揮Shearlet變換的優勢，不能達到滿意的去噪效果。雖然地震數據的分布與Shearlet變換的尺度和方向相關性很強，但是在確定的尺度內，卻有相似的分布規律，即在一些方向上分布大量的地震數據，而總存在一個方向，該方向地震數據分布最少。而在含有隨機噪聲的地震數據中，因為隨機噪聲的分布具有隨機性，所以在有效信號分布最少的方向上就會存在有大量的隨機噪聲，用此噪聲來預估實際噪聲，并從地震數據中去除噪聲，分離出有效地震信號。為此，本文分析了地震數據在Shearlet域不同尺度和不同方向上分布的差異，提出一種衡量地震數據噪聲的方法，即將有效信號分布最少的尺度和方向上分布的噪聲作為該尺度的噪聲，并將其剔除，得到分布在各尺度和方向上的有效地震數據。

類比式(11)閾值隨尺度變化的規律設置方向自適應項，結合尺度自適應性，提出一種新的隨尺度和方向自適應變化的閾值函數

(14)

式中：ej,l為尺度j、方向l上的Shearlet變換系數矩陣的L2范數，j=1,2,…,J，J為分解總尺度；σ為噪聲標準差；η為1～10的整數，微調去噪效果； min(·)為最小值函數。

以一合成理論單炮地震記錄(圖2)為例說明本文方法的實現過程。該地震記錄共包含250道，每道有1000個采樣點，采樣頻率為1ms。對該數據進行Shearlet變換，用于分析地震數據在Shearlet域不同分解尺度和方向上的分布特征。圖3為求取各分解尺度和各方向的Shearlet系數矩陣的L2范數，發現在不同尺度上地震數據的分布存在差異，并且在確定尺度內不同方向上的地震數據也存在著更明顯差異。對圖3進行排序得到地震數據經Shearlet變換后在不同尺度、不同方向上的分布規律(圖4)。

圖2 合成理論地震數據

為了驗證本文自適應閾值去噪方法的優越性，分別對一組理論地震數據和疊前、疊后實際地震數據應用本文方法和傳統閾值法進行噪聲壓制試驗，比較二者的去噪效果，并統計兩種方法去噪后的信噪比、峰值信噪比和均方差。

圖3 不同尺度和方向Shearlet變換系數矩陣L2范數柱狀圖

圖4 不同尺度和方向Shearlet變換系數矩陣L2范數折線圖

3 理論數據測試

利用理論單炮地震記錄進行噪聲壓制試驗。試驗選擇正交小波作為分解母函數，設置4個分解層，角度參數為[1,1,2,2]。理論單炮記錄如圖2所示。在理論地震數據中加入白噪聲后，地震數據信噪比為5.1539dB(圖5a)。在單炮記錄中分布著大量的隨機噪聲，能量較弱的同向軸被隨機噪聲淹沒，無法有效識別。分別利用基于Shearlet變換的傳統硬閾值和自適應閾值兩種去噪方法對合成的單炮記錄進行隨機噪聲處理，結果如圖5b和圖5c所示。

分析圖5可知，自適應閾值方法能有效壓制隨機噪聲，處理后的整個地震記錄很干凈，能夠清楚識別出有效信號，同向軸的連續性也很好。圖6為圖2和圖5紅框區域局部放大圖，可以發現，模擬理論地震數據在700～1000ms存在能量較弱的有效信號，添加噪聲后能量較弱的同向軸被噪聲淹沒。硬閾值去噪剖面(圖6c)上，箭頭所指同向軸幾乎不可分辨； 而利用本文的自適應閾值進行噪聲壓制(圖6d)，同向軸連續性較好，有效信號得到了很好的重構，說明該方法比硬閾值能更好地保護有效信號、尤其是弱信號的能力。

分析去噪差剖面(圖7),硬閾值法雖然去除了隨機噪聲，但不能識別能量較弱的有效信號，造成有效信號損失嚴重，而且差剖面上也有較多的同向軸出現(圖7a)，表明硬閾值不能準確區分信號與噪聲，去噪的同時對有效信號造成損傷；用本文方法得到的去噪差剖面(圖7b)，幾乎沒有同向軸出現，說明該方法能夠有效區分信號與噪聲，在去除隨機噪聲的同時保留更多的有效地震信息，獲得高信噪比的地震數據。

圖5 兩種閾值去噪結果對比

圖6 圖2和圖5中紅色方框區域的局部放大

圖7 兩種閾值方法所去除的噪聲道集

表1給出了兩種閾值去噪方法去除的理論含噪數據的峰值信噪比(PSNR)、信噪比(SNR)和均方誤差(MSE)。可見，自適應閾值法相比硬閾值法，去噪后的剖面具有較大的SNR和PSNR，但MSE卻更小，說明自適應閾值去噪效果優于硬閾值法。綜上所述，本文提出的自適應閾值法相比于硬閾值法能夠更有效去除理論地震數據中的噪聲。

表1 兩種閾值去噪效果指標

4 實際地震數據去噪效果分析

4.1 疊后地震數據

圖8a為一褶皺構造的實際地震疊后剖面，數據含有大量的噪聲，同向軸的連續性較差，部分能量較弱的有效信號被噪聲淹沒。

分別利用基于Shearlet變換傳統硬閾值和自適應閾值對上述實際地震數據進行隨機噪聲壓制，結果分別如圖8b和圖8c所示。與圖8a對比可以發現，圖8b中地震剖面干凈平整，說明隨機噪聲基本被剔除，但是同向軸的連續性很差；圖8c中，同向軸的連續性更好，且邊界信息更加清晰，凸顯了被隨機噪聲淹沒的能量較弱的同向軸(圖中箭頭所指位置)。圖9為兩種方法去噪的差值剖面，相比于硬閾值去噪差剖面(圖9a)，自適應閾值剖面(圖9b)中幾乎無有效地震信號分布，說明自適應閾值方法具有更好的保留地震信號的能力，能夠有效提高去噪后地震數據的信噪比。

圖8 實際疊后地震數據去噪效果分析

4.2 疊前地震數據

實際采集的疊前地震記錄如圖10a所示，分別應用硬閾值法和自適應閾值法對其進行去噪處理，結果見圖10b和圖10c。圖11 是對圖10中紅色線框區域的放大圖。圖12是不同方法所剔除的噪聲。對比可以發現，經自適應閾值法去噪后，隨機噪聲得到很好的壓制，地震道集更平整，同相軸也更清晰，且能夠更好地分辨出深部能量較弱的有效地震信號，豐富了地震數據的細節信息。經自適應閾值法處理后得到的差值相比傳統硬閾值法明顯地具有更少的同相軸。表明本文方法去除了隨機噪聲的同時能在較大程度上保護地震數據的完整性，實現比硬閾值更好的去噪效果。

圖9 實際疊后地震數據不同方法去除的噪聲剖面

圖10 實際疊前地震數據去噪效果分析

圖11 圖10紅色方框部分的局部放大

圖12a 實際疊前地震數據硬閾值法去除的噪聲

圖12b 實際疊前地震數據自適應閾值法去除的噪聲

5 結論與認識

傳統閾值法壓制隨機噪聲只考慮地震信號的分布在尺度上存在差異，因而只能去除部分噪聲，閾值設置不合理。為了提高去噪效果，本文統計并分析了理論地震數據在Shearlet域尺度和方向的分布規律，提出一種估計隨機去噪方法，并且設置尺度和方向自適應閾值壓制隨機噪聲。基于理論數據和實際地震數據的計算結果，得到以下結論。

(1)Shearlet變換具有多尺度、多方向特性，地震數據在Shearlet域不同尺度和不同方向存在差異。對不同方向上地震數據的統計分析表明，在每一尺度內，總存在某一方向上的地震數據分布最少，而噪聲在該方向的分布相對較多。

(2)因為隨機噪聲的分布具有隨機性，所以將噪聲分布相對較多的方向(即該方向上有效地震信號最少)的噪聲強度作為該尺度內隨機噪聲的強度。這樣可最大程度上分離噪聲，得到較為純凈的有效信號。

(3)基于地震數據在Shearlet域各尺度和方向的分布差異，設置尺度和方向自適應閾值，去除隨機噪聲。參考尺度對閾值影響的處理方式，構造方向自適應項，提出一種新的隨尺度和方向自適應變化的閾值，解決了傳統閾值去噪過程中損失部分有效信號的問題。

(4)利用兩種不同的閾值法對理論地震數據和實際地震數據進行隨機噪聲壓制試驗，結果表明本文提出的基于Shearlet變換的自適應閾值估算法相比于傳統的硬閾值法能更有效去除地震數據中的隨機噪聲，保留更多的有效信號。