薛 林 程 浩* 鞏恩普 陳毅軍
(①東北大學(xué)深部金屬礦山安全開采教育部重點(diǎn)實(shí)驗(yàn)室,遼寧沈陽 110004; ②東北大學(xué)資源與土木工程學(xué)院,遼寧沈陽 110004)
地震勘探在地質(zhì)調(diào)查和資源探測方面具有重要作用[1]。實(shí)際地震數(shù)據(jù)往往混雜隨機(jī)噪聲,嚴(yán)重影響數(shù)據(jù)的分析和處理效果。所以有效去除噪聲、獲得高信噪比的地震數(shù)據(jù)是地震勘探首先要解決的問題[2-4]。
隨機(jī)噪聲壓制主要是基于有效信號與隨機(jī)噪聲在頻率或視速度方面的差異性,常用的方法有中值濾波[5-7]、奇異值分解[8-10]、f-x域預(yù)測濾波[11-13]等。小波變換[14-16]作為一種多尺度時(shí)頻分析手段,能夠多尺度分析地震數(shù)據(jù),但在處理多維地震數(shù)據(jù)時(shí)效果不明顯。而Curvelet變換[17-19]和Contourlet變換[20-22]能夠更好地描述地震信號的局部特征,具有突出的稀疏表示能力,可以實(shí)現(xiàn)很好的去噪效果。
Guo等[23-24]進(jìn)一步研究了小波理論,結(jié)合多尺度幾何分析和膨脹仿射系統(tǒng),提出了Shearlet變換理論。Shearlet變換繼承了小波變換的局部性和多尺度特性。由于Shearlet變換具有很強(qiáng)的方向性,而且其頻域支撐是各向異性的,因此Shearlet變換可以用最少的稀疏系數(shù)逼近奇異曲線,實(shí)現(xiàn)信號的最優(yōu)稀疏表示。同時(shí)Shearlet變換可多尺度描述高維信號的幾何特性,且方向的數(shù)目隨著尺度的增加而增加,因而可以克服小波變換對多維地震數(shù)據(jù)處理效果不佳的缺陷。
由于Shearlet變換獨(dú)特的優(yōu)勢,其在多個(gè)領(lǐng)域得到了廣泛的使用。鄭晶等[25]首先對地震數(shù)據(jù)進(jìn)行奇異值分解(SVD)去除直達(dá)波,然后通過Shearlet變換改進(jìn)閾值壓制隨機(jī)噪聲;李娟等[26]對地震數(shù)據(jù)進(jìn)行Shearlet分解,通過增大信號分布集中的方向系數(shù)來增強(qiáng)信號與噪聲的差異,然后利用峰度特性拾取地震初至波;鄭升等[27]對大量含隨機(jī)噪聲的沙漠地震數(shù)據(jù)進(jìn)行Shearlet變換,將得到的系數(shù)作為輸入標(biāo)簽用于訓(xùn)練深度殘差卷積神經(jīng)網(wǎng)絡(luò)模型,得到預(yù)測噪聲和有效信號的變換系數(shù);王憲楠等[28]利用探地雷達(dá)數(shù)據(jù)在Shearlet域的差異,采用硬閾值壓制探地雷達(dá)數(shù)據(jù)中的隨機(jī)噪聲;劉昕等[29]將非下采樣與Shearlet變換結(jié)合以壓制隨機(jī)噪聲,消除傳統(tǒng)Shearlet變換去噪中存在的吉布斯現(xiàn)象;吳一全等[30]基于非下采樣Shearlet變換對圖像的低頻和高頻分量分別進(jìn)行奇異值分解和各向異性擴(kuò)散,以達(dá)到去噪目的,很好地保留了圖像的細(xì)節(jié)信息。
傳統(tǒng)的Shearlet變換閾值噪聲壓制方法,主要考慮信號在Shearlet域尺度層面的特性,針對不同的尺度層設(shè)置自適應(yīng)閾值,卻忽略了地震數(shù)據(jù)在方向上的變化特性。本文通過分析信號在Shearlet域方向上的分布特點(diǎn),提出閾值隨尺度和方向同時(shí)自適應(yīng)變化的噪聲壓制方法。分別對一組理論模型和兩組實(shí)際數(shù)據(jù)進(jìn)行噪聲壓制試驗(yàn),并與傳統(tǒng)閾值去噪結(jié)果對比,驗(yàn)證本文的尺度和方向自適應(yīng)閾值法能有效壓制噪聲,最大程度上保留有效信號,提高地震數(shù)據(jù)的信噪比。
Guo等[31]基于合成小波理論,結(jié)合仿射系統(tǒng)和多尺度幾何分析,提出了帶有合成膨脹系統(tǒng)的幾何變換方法。合成連續(xù)膨脹系統(tǒng)在二維空間表達(dá)式為
ψ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構(gòu)成
(2)
則可以得到如下連續(xù)仿射系統(tǒng)
ψH(φ)=ψa,b,k(φ)
a∈R+,b∈R,k∈R2}
(3)
式中:t為平移參數(shù);ψa,b,k(φ)為連續(xù)剪切波;ψa,b,k(x)為剪切母函數(shù)。剪切波ψa,b,k(φ)在各個(gè)尺度層對應(yīng)的頻域劃分如圖1所示。
基于剪切波函數(shù),可以定義連續(xù)Shearlet變換為
C=SH{f}=〈f,φi,l,k〉
(4)
式中:C為Shearlet系數(shù)矩陣; SH{·}為Shearlet變換; 〈·,·〉表示內(nèi)積;f為含噪信號。
同樣可以對Shearlet系數(shù)矩陣C進(jìn)行反變換重構(gòu)原始函數(shù)
f=SH-1{C}
(5)
式中SH-1表示Shearlet逆變換。
圖1 Shearlet頻域剖分示意圖
假設(shè)尺度參數(shù)a=2-j,剪切參數(shù)b=-l(j,l∈Z),用點(diǎn)k∈Z2代替連續(xù)Shearlet變換中的平移參數(shù)t∈R2,得到Shearlet變換基元素的表達(dá)式為
(6)
對于任意f∈L2(R2),離散Shearlet變換需滿足
(7)
基于Shearlet變換最優(yōu)稀疏表示的特點(diǎn),地震數(shù)據(jù)經(jīng)Shearlet變換后得到各分解尺度和方向的Shearlet系數(shù)。若Shearlet基函數(shù)的方向與地震信號的方向大致相同,則對應(yīng)較大的Shearlet系數(shù);當(dāng)兩方向相差較大時(shí),Shearlet系數(shù)偏小。在實(shí)際地震數(shù)據(jù)中,有效信號分布具有方向性,而隨機(jī)噪聲是隨機(jī)分布的,所以地震數(shù)據(jù)經(jīng)Shearlet變換后,有效信號相比隨機(jī)噪聲會(huì)對應(yīng)較大的Shearlet系數(shù)。這樣就可以設(shè)置一個(gè)閾值,利用閾值函數(shù)去除隨機(jī)噪聲對應(yīng)的較小的系數(shù),達(dá)到去除隨機(jī)噪聲的目的。然后對處理后的Shearlet系數(shù)進(jìn)行逆變換,得到去噪后的地震信號。
假設(shè)待處理的含噪地震數(shù)據(jù)模型為
f(t)=u(t)+n(t)
(8)
式中:u(t)為有效信號;n(t)為隨機(jī)噪聲。Shearlet變換閾值去噪算法的實(shí)質(zhì)就是從f(t)中得到u(t)。該過程可以分以下三個(gè)步驟。
(1)對含噪地震數(shù)據(jù)f(t)進(jìn)行Shearlet變換,獲取Shearlet系數(shù)
C(j,l,k)=SH{f(t)}=〈f(t),φj,l,k〉
(9)
(2)對得到的Shearlet系數(shù)進(jìn)行閾值化
(10)
(11)
(12)
式中:Cnew(j,l,k)表示閾值化處理后的系數(shù);N表示地震數(shù)據(jù)的總采樣點(diǎn)數(shù); median(·)表示對矩陣所有元素取中值。
(3)對閾值化處理后的系數(shù)Cnew進(jìn)行逆變換,得到去噪后的有效信號
(13)
傳統(tǒng)閾值函數(shù)在壓制隨機(jī)噪聲時(shí),僅對不同分解尺度設(shè)置不同的閾值,實(shí)現(xiàn)了閾值隨分解尺度的自適應(yīng)。但因?yàn)镾healret變換具有多尺度、多方向的特點(diǎn),所以僅隨尺度自適應(yīng)閾值不能充分發(fā)揮Shearlet變換的優(yōu)勢,不能達(dá)到滿意的去噪效果。雖然地震數(shù)據(jù)的分布與Shearlet變換的尺度和方向相關(guān)性很強(qiáng),但是在確定的尺度內(nèi),卻有相似的分布規(guī)律,即在一些方向上分布大量的地震數(shù)據(jù),而總存在一個(gè)方向,該方向地震數(shù)據(jù)分布最少。而在含有隨機(jī)噪聲的地震數(shù)據(jù)中,因?yàn)殡S機(jī)噪聲的分布具有隨機(jī)性,所以在有效信號分布最少的方向上就會(huì)存在有大量的隨機(jī)噪聲,用此噪聲來預(yù)估實(shí)際噪聲,并從地震數(shù)據(jù)中去除噪聲,分離出有效地震信號。為此,本文分析了地震數(shù)據(jù)在Shearlet域不同尺度和不同方向上分布的差異,提出一種衡量地震數(shù)據(jù)噪聲的方法,即將有效信號分布最少的尺度和方向上分布的噪聲作為該尺度的噪聲,并將其剔除,得到分布在各尺度和方向上的有效地震數(shù)據(jù)。
類比式(11)閾值隨尺度變化的規(guī)律設(shè)置方向自適應(yīng)項(xiàng),結(jié)合尺度自適應(yīng)性,提出一種新的隨尺度和方向自適應(yīng)變化的閾值函數(shù)
(14)
式中:ej,l為尺度j、方向l上的Shearlet變換系數(shù)矩陣的L2范數(shù),j=1,2,…,J,J為分解總尺度;σ為噪聲標(biāo)準(zhǔn)差;η為1~10的整數(shù),微調(diào)去噪效果; min(·)為最小值函數(shù)。
以一合成理論單炮地震記錄(圖2)為例說明本文方法的實(shí)現(xiàn)過程。該地震記錄共包含250道,每道有1000個(gè)采樣點(diǎn),采樣頻率為1ms。對該數(shù)據(jù)進(jìn)行Shearlet變換,用于分析地震數(shù)據(jù)在Shearlet域不同分解尺度和方向上的分布特征。圖3為求取各分解尺度和各方向的Shearlet系數(shù)矩陣的L2范數(shù),發(fā)現(xiàn)在不同尺度上地震數(shù)據(jù)的分布存在差異,并且在確定尺度內(nèi)不同方向上的地震數(shù)據(jù)也存在著更明顯差異。對圖3進(jìn)行排序得到地震數(shù)據(jù)經(jīng)Shearlet變換后在不同尺度、不同方向上的分布規(guī)律(圖4)。
圖2 合成理論地震數(shù)據(jù)
為了驗(yàn)證本文自適應(yīng)閾值去噪方法的優(yōu)越性,分別對一組理論地震數(shù)據(jù)和疊前、疊后實(shí)際地震數(shù)據(jù)應(yīng)用本文方法和傳統(tǒng)閾值法進(jìn)行噪聲壓制試驗(yàn),比較二者的去噪效果,并統(tǒng)計(jì)兩種方法去噪后的信噪比、峰值信噪比和均方差。
圖3 不同尺度和方向Shearlet變換系數(shù)矩陣L2范數(shù)柱狀圖
圖4 不同尺度和方向Shearlet變換系數(shù)矩陣L2范數(shù)折線圖
利用理論單炮地震記錄進(jìn)行噪聲壓制試驗(yàn)。試驗(yàn)選擇正交小波作為分解母函數(shù),設(shè)置4個(gè)分解層,角度參數(shù)為[1,1,2,2]。理論單炮記錄如圖2所示。在理論地震數(shù)據(jù)中加入白噪聲后,地震數(shù)據(jù)信噪比為5.1539dB(圖5a)。在單炮記錄中分布著大量的隨機(jī)噪聲,能量較弱的同向軸被隨機(jī)噪聲淹沒,無法有效識(shí)別。分別利用基于Shearlet變換的傳統(tǒng)硬閾值和自適應(yīng)閾值兩種去噪方法對合成的單炮記錄進(jìn)行隨機(jī)噪聲處理,結(jié)果如圖5b和圖5c所示。
分析圖5可知,自適應(yīng)閾值方法能有效壓制隨機(jī)噪聲,處理后的整個(gè)地震記錄很干凈,能夠清楚識(shí)別出有效信號,同向軸的連續(xù)性也很好。圖6為圖2和圖5紅框區(qū)域局部放大圖,可以發(fā)現(xiàn),模擬理論地震數(shù)據(jù)在700~1000ms存在能量較弱的有效信號,添加噪聲后能量較弱的同向軸被噪聲淹沒。硬閾值去噪剖面(圖6c)上,箭頭所指同向軸幾乎不可分辨; 而利用本文的自適應(yīng)閾值進(jìn)行噪聲壓制(圖6d),同向軸連續(xù)性較好,有效信號得到了很好的重構(gòu),說明該方法比硬閾值能更好地保護(hù)有效信號、尤其是弱信號的能力。
分析去噪差剖面(圖7),硬閾值法雖然去除了隨機(jī)噪聲,但不能識(shí)別能量較弱的有效信號,造成有效信號損失嚴(yán)重,而且差剖面上也有較多的同向軸出現(xiàn)(圖7a),表明硬閾值不能準(zhǔn)確區(qū)分信號與噪聲,去噪的同時(shí)對有效信號造成損傷;用本文方法得到的去噪差剖面(圖7b),幾乎沒有同向軸出現(xiàn),說明該方法能夠有效區(qū)分信號與噪聲,在去除隨機(jī)噪聲的同時(shí)保留更多的有效地震信息,獲得高信噪比的地震數(shù)據(jù)。
圖5 兩種閾值去噪結(jié)果對比
圖6 圖2和圖5中紅色方框區(qū)域的局部放大
圖7 兩種閾值方法所去除的噪聲道集
表1給出了兩種閾值去噪方法去除的理論含噪數(shù)據(jù)的峰值信噪比(PSNR)、信噪比(SNR)和均方誤差(MSE)??梢姡赃m應(yīng)閾值法相比硬閾值法,去噪后的剖面具有較大的SNR和PSNR,但MSE卻更小,說明自適應(yīng)閾值去噪效果優(yōu)于硬閾值法。綜上所述,本文提出的自適應(yīng)閾值法相比于硬閾值法能夠更有效去除理論地震數(shù)據(jù)中的噪聲。
表1 兩種閾值去噪效果指標(biāo)
圖8a為一褶皺構(gòu)造的實(shí)際地震疊后剖面,數(shù)據(jù)含有大量的噪聲,同向軸的連續(xù)性較差,部分能量較弱的有效信號被噪聲淹沒。
分別利用基于Shearlet變換傳統(tǒng)硬閾值和自適應(yīng)閾值對上述實(shí)際地震數(shù)據(jù)進(jìn)行隨機(jī)噪聲壓制,結(jié)果分別如圖8b和圖8c所示。與圖8a對比可以發(fā)現(xiàn),圖8b中地震剖面干凈平整,說明隨機(jī)噪聲基本被剔除,但是同向軸的連續(xù)性很差;圖8c中,同向軸的連續(xù)性更好,且邊界信息更加清晰,凸顯了被隨機(jī)噪聲淹沒的能量較弱的同向軸(圖中箭頭所指位置)。圖9為兩種方法去噪的差值剖面,相比于硬閾值去噪差剖面(圖9a),自適應(yīng)閾值剖面(圖9b)中幾乎無有效地震信號分布,說明自適應(yīng)閾值方法具有更好的保留地震信號的能力,能夠有效提高去噪后地震數(shù)據(jù)的信噪比。
圖8 實(shí)際疊后地震數(shù)據(jù)去噪效果分析
實(shí)際采集的疊前地震記錄如圖10a所示,分別應(yīng)用硬閾值法和自適應(yīng)閾值法對其進(jìn)行去噪處理,結(jié)果見圖10b和圖10c。圖11 是對圖10中紅色線框區(qū)域的放大圖。圖12是不同方法所剔除的噪聲。對比可以發(fā)現(xiàn),經(jīng)自適應(yīng)閾值法去噪后,隨機(jī)噪聲得到很好的壓制,地震道集更平整,同相軸也更清晰,且能夠更好地分辨出深部能量較弱的有效地震信號,豐富了地震數(shù)據(jù)的細(xì)節(jié)信息。經(jīng)自適應(yīng)閾值法處理后得到的差值相比傳統(tǒng)硬閾值法明顯地具有更少的同相軸。表明本文方法去除了隨機(jī)噪聲的同時(shí)能在較大程度上保護(hù)地震數(shù)據(jù)的完整性,實(shí)現(xiàn)比硬閾值更好的去噪效果。
圖9 實(shí)際疊后地震數(shù)據(jù)不同方法去除的噪聲剖面
圖10 實(shí)際疊前地震數(shù)據(jù)去噪效果分析
圖11 圖10紅色方框部分的局部放大
圖12a 實(shí)際疊前地震數(shù)據(jù)硬閾值法去除的噪聲
圖12b 實(shí)際疊前地震數(shù)據(jù)自適應(yīng)閾值法去除的噪聲
傳統(tǒng)閾值法壓制隨機(jī)噪聲只考慮地震信號的分布在尺度上存在差異,因而只能去除部分噪聲,閾值設(shè)置不合理。為了提高去噪效果,本文統(tǒng)計(jì)并分析了理論地震數(shù)據(jù)在Shearlet域尺度和方向的分布規(guī)律,提出一種估計(jì)隨機(jī)去噪方法,并且設(shè)置尺度和方向自適應(yīng)閾值壓制隨機(jī)噪聲?;诶碚摂?shù)據(jù)和實(shí)際地震數(shù)據(jù)的計(jì)算結(jié)果,得到以下結(jié)論。
(1)Shearlet變換具有多尺度、多方向特性,地震數(shù)據(jù)在Shearlet域不同尺度和不同方向存在差異。對不同方向上地震數(shù)據(jù)的統(tǒng)計(jì)分析表明,在每一尺度內(nèi),總存在某一方向上的地震數(shù)據(jù)分布最少,而噪聲在該方向的分布相對較多。
(2)因?yàn)殡S機(jī)噪聲的分布具有隨機(jī)性,所以將噪聲分布相對較多的方向(即該方向上有效地震信號最少)的噪聲強(qiáng)度作為該尺度內(nèi)隨機(jī)噪聲的強(qiáng)度。這樣可最大程度上分離噪聲,得到較為純凈的有效信號。
(3)基于地震數(shù)據(jù)在Shearlet域各尺度和方向的分布差異,設(shè)置尺度和方向自適應(yīng)閾值,去除隨機(jī)噪聲。參考尺度對閾值影響的處理方式,構(gòu)造方向自適應(yīng)項(xiàng),提出一種新的隨尺度和方向自適應(yīng)變化的閾值,解決了傳統(tǒng)閾值去噪過程中損失部分有效信號的問題。
(4)利用兩種不同的閾值法對理論地震數(shù)據(jù)和實(shí)際地震數(shù)據(jù)進(jìn)行隨機(jī)噪聲壓制試驗(yàn),結(jié)果表明本文提出的基于Shearlet變換的自適應(yīng)閾值估算法相比于傳統(tǒng)的硬閾值法能更有效去除地震數(shù)據(jù)中的隨機(jī)噪聲,保留更多的有效信號。