王通,劉建勛,王興宇,李廣才,田密
(1.國家現(xiàn)代地質(zhì)勘查工程技術(shù)研究中心,河北 廊坊 065000;2.中國地質(zhì)科學(xué)院 地球物理地球化學(xué)勘查研究所,河北 廊坊 065000)
在深地探測領(lǐng)域中,深反射地震技術(shù)被認(rèn)為是探測巖石圈、解決深部地質(zhì)問題最為準(zhǔn)確有效的方法之一[1],深反射地震探測方法采用大藥量激發(fā),長時(shí)間觀測,地震波能夠帶來莫霍面及以下的深部信息,深反射地震數(shù)據(jù)能夠?qū)ι畈康刭|(zhì)結(jié)構(gòu)實(shí)現(xiàn)有效的成像。
深反射地震接收來自地層深部的反射信號(hào),由于大地濾波作用使得反射信號(hào)能量很弱,受背景隨機(jī)噪聲干擾嚴(yán)重,一些微弱的構(gòu)造響應(yīng)會(huì)被隨機(jī)噪聲淹沒,對深反射地震數(shù)據(jù)處理和解釋帶來很大影響。深反射地震數(shù)據(jù)的處理結(jié)果直接關(guān)系到地震剖面的整體成像質(zhì)量,高分辨率地震數(shù)據(jù)處理技術(shù)對深反射地震探測至關(guān)重要[2-3]。朱小三等[4]指出怎樣除去深部信號(hào)中噪聲的同時(shí)又能最大限度地保留有效信號(hào)的能量是深反射地震資料去噪的關(guān)鍵。
通過對深反射地震數(shù)據(jù)的處理發(fā)現(xiàn),常規(guī)處理(異常振幅噪聲衰減、局部異常振幅壓制、相干噪聲和非相干噪聲壓制、球面擴(kuò)散補(bǔ)償、地表一致性補(bǔ)償、地震子波反褶積等)后深反射地震數(shù)據(jù)中大部分噪聲能夠得到很好的壓制,但是來自莫霍面及以下構(gòu)造的弱反射信號(hào)不能有效地從背景噪聲中呈現(xiàn)出來。傳統(tǒng)的隨機(jī)噪聲去除方法(頻率域?yàn)V波、f-x域去噪、Radon域去噪等)不能滿足深部弱信號(hào)高精度成像的需求。
20世紀(jì)80年代,小波變換(Wavelet)[5-6]作為一種多尺度變換工具,在去噪和圖形處理領(lǐng)域得到了快速的發(fā)展,利用小波變換方法分尺度對地震數(shù)據(jù)中隨機(jī)噪聲進(jìn)行壓制取得了很好的應(yīng)用效果。程浩等[7]采用基于小波變換的自適應(yīng)閾值,利用分層自適應(yīng)因子,對礦山微震信號(hào)中所包含的隨機(jī)噪聲進(jìn)行壓制。羅紅梅等[8]采用基于小波變換的經(jīng)驗(yàn)?zāi)B(tài)EMD分解的地震記錄重建方法,較好地提高了弱信號(hào)的能量信息。陳亮等[9]提出一種基于神經(jīng)網(wǎng)絡(luò)改進(jìn)小波的地震數(shù)據(jù)隨機(jī)噪聲去除方法,采用級(jí)聯(lián)BP神經(jīng)網(wǎng)絡(luò)模型提取出多類別隨機(jī)噪聲信號(hào),實(shí)現(xiàn)地震數(shù)據(jù)的隨機(jī)信號(hào)壓制?;谛〔ㄗ儞Q的隨機(jī)噪聲壓制方法實(shí)現(xiàn)了噪聲閾值的尺度自適應(yīng),但是還不滿足地震數(shù)據(jù)多方向性的需求[10-12]。
在Wavelet多尺度變換的基礎(chǔ)上,多尺度多方向的Curvelet變換應(yīng)運(yùn)而生,有效的應(yīng)用到地震數(shù)據(jù)隨機(jī)噪聲壓制中。劉磊等[13]、金丹等[14]、何柯等[15]利用Curvelet域的多種閾值方法對二維地震數(shù)據(jù)實(shí)現(xiàn)了弱信號(hào)提取和隨機(jī)噪聲壓制,張華等[16]采用曲波變換對地震數(shù)據(jù)抽取時(shí)間切片,實(shí)現(xiàn)了三維地震數(shù)據(jù)隨機(jī)噪聲壓制。孫苗苗等[17]采用基于曲波域稀疏約束的OVT域地震數(shù)據(jù)去噪方法,有效提取了OVT域地震數(shù)據(jù)的中、深反射層的弱有效信號(hào),在壓制強(qiáng)隨機(jī)噪聲的同時(shí)減少了弱有效信號(hào)的損失。但是,Curvelet變換的數(shù)學(xué)結(jié)構(gòu)復(fù)雜,且在變換的某一尺度角度域中系數(shù)矩陣尺度與原始數(shù)據(jù)不一致。
相較于Curvelet變換,Shearlet變換[18]數(shù)學(xué)結(jié)構(gòu)更為簡單,同樣具有多尺度、多方向特性,利用更少的系數(shù)逼近曲線。對于一個(gè)N×N階的圖像,在Shearlet域中各尺度方向的數(shù)據(jù)體仍是N×N階,且各方向系數(shù)仍然保留時(shí)間的維度。程浩等[19]采用尺度自適應(yīng)三維Shearlet變換對地震數(shù)據(jù)中隨機(jī)噪聲進(jìn)行壓制。童思友等[20]在Shearlet域局部閾值的基礎(chǔ)上改進(jìn)貝葉斯閾值,將信噪比與閾值函數(shù)有機(jī)關(guān)聯(lián),并將信噪比作為閾值設(shè)定的因素,實(shí)現(xiàn)了地震資料中隨機(jī)噪聲的壓制。薛琳等[21]結(jié)合Shearlet變換多尺度、多方向特性,提出一種隨尺度和方向同時(shí)自適應(yīng)變化的閾值實(shí)現(xiàn)地震數(shù)據(jù)隨機(jī)噪聲的壓制。以上對Shearlet域內(nèi)隨機(jī)噪聲估計(jì)的參數(shù)過于單一,對于深層弱反射信號(hào)下隨機(jī)噪聲壓制效果欠佳。
本文應(yīng)用結(jié)合信噪比(SNR)、L2范數(shù)、隨機(jī)噪聲殘差的Shearlet域尺度角度自適應(yīng)閾值方法對深反射地震數(shù)據(jù)中隨機(jī)噪聲進(jìn)行壓制,將地震數(shù)據(jù)變換到Shearlet域中,由于地震數(shù)據(jù)具有方向性,而隨機(jī)噪聲是無規(guī)律分布的,所以同一尺度下不同角度上有效信號(hào)與隨機(jī)噪聲比存在微弱差異,在此基礎(chǔ)上利用信噪比(SNR)、隨機(jī)噪聲殘差和Shearlet系數(shù)L2范數(shù)對每個(gè)角度域數(shù)據(jù)進(jìn)行自適應(yīng)閾值求取,實(shí)現(xiàn)更為精準(zhǔn)的隨尺度角度自適應(yīng)變化的隨機(jī)噪聲壓制方法。
Shearlet變換是基于合成小波理論,結(jié)合仿射系統(tǒng)和多尺度幾何分析,帶有合成膨脹系統(tǒng)的幾何變換方法,是一種多尺度幾何分析工具,具有最優(yōu)的非線性誤差逼近性能,數(shù)學(xué)結(jié)構(gòu)嚴(yán)謹(jǐn),計(jì)算復(fù)雜度低。連續(xù)的仿射系統(tǒng)在L2(Rn)具有以下形式[18]:
{TtDMψ,M∈G,t∈Rn}
(1)
式中ψ∈L2(Rn),Tt是平移算子,并且定義:Ttf(x)=f(x-t),DM表示伸縮算子,定義:DMf(x)=|detM|-1/2f(M-1x),G?GLn(R)是一組矩陣,當(dāng)n=2,G是含有2個(gè)參數(shù)的擴(kuò)張組:
(2)
矩陣Mas可以表示為Mas=SsAa,式中S為剪切矩陣,A是尺度矩陣:
(3)
式中:s為剪切參數(shù),a為尺度參數(shù)。令i,j,k分別代表尺度、方向和系數(shù)位置序號(hào),則由i,j,k定義的實(shí)數(shù)域連續(xù)可積函數(shù)為:
φi,j,k(x)=|detA|i/2φ(SjAkx-k) ,
(4)
式中:x為自變量。若φ滿足Parserval框架,則對于任何一個(gè)平方可積的函數(shù)f的Shearlet變換表示為:
(5)
式中:C為Shearlet系數(shù)矩陣;SH{·}為Shearlet變換;〈·,·〉表示內(nèi)積;f為含噪信號(hào)。
同樣可以對Shearlet系數(shù)矩陣C進(jìn)行反變換重構(gòu)原始函數(shù):
f=SH-1{C} ,
(6)
式中,SH-1表示Shearlet逆變換。
假設(shè)地震記錄中有效信號(hào)為P0(t),隨機(jī)噪聲為n(t),則觀測到的包含隨機(jī)噪聲的地震數(shù)據(jù)為:
P(t)=P0(t)+n(t) 。
(7)
傳統(tǒng)閾值法壓制隨機(jī)噪聲是通過計(jì)算地震數(shù)據(jù),篩選出一個(gè)隨機(jī)噪聲強(qiáng)度相應(yīng)的閾值參數(shù)Th,這里Th是一個(gè)常數(shù)。利用閾值參數(shù)Th對地震數(shù)據(jù)分選,大于Th的數(shù)據(jù)被視為有效信號(hào)進(jìn)行保留,小于Th的數(shù)據(jù)被視為隨機(jī)噪聲進(jìn)行壓制,最終獲得去除隨機(jī)噪聲后的地震數(shù)據(jù)。
Shearlet是一種稀疏變換方法,地震數(shù)據(jù)經(jīng)過Shearlet變換后得到分尺度、分方向表示稀疏的Shearlet系數(shù)。Shearlet變換具有多尺度多角度特性,在特定的尺度角度中,當(dāng)?shù)卣饠?shù)據(jù)同相軸走向與Shearlet變換基函數(shù)的方向大致相同時(shí),地震信號(hào)就對應(yīng)著較大的Shearlet系數(shù);當(dāng)?shù)卣饠?shù)據(jù)同相軸走向與Shearlet變換基函數(shù)的方向相差較大時(shí),相應(yīng)的Shearlet系數(shù)較小。在實(shí)際地震數(shù)據(jù)中,有效的地震信號(hào)通常具有方向性,對應(yīng)著較大的Shearlet系數(shù),而隨機(jī)噪聲并無規(guī)律而言,通常對應(yīng)著較小的Shearlet系數(shù)。對含噪聲的地震數(shù)據(jù)P(t)進(jìn)行Shearlet稀疏變換,獲得對應(yīng)的Shearlet系數(shù),如下:
Cp(j,l,k)=SH{P(t)}=〈P(t),φj,l,k〉 。
(8)
利用稀疏變換后的Shearlet系數(shù)估計(jì)隨機(jī)噪聲的方差,再通過閾值估計(jì)算子獲得對應(yīng)的閾值參數(shù):
(9)
(10)
式中:N表示地震數(shù)據(jù)的總采樣點(diǎn)數(shù);median(·)表示對數(shù)據(jù)矩陣中所有元素取中值。
對稀疏變換后的Shearlet系數(shù)進(jìn)行閾值處理
(11)
式中:CTh(j,l,k)表示閾值處理后的Shearlet系數(shù)矩陣,保留較大的Shearlet系數(shù),去除包含隨機(jī)噪聲較小的Shearlet系數(shù),就實(shí)現(xiàn)了Shearlet變換閾值法去除隨機(jī)噪聲的目的[21]。
式(10)利用隨機(jī)噪聲的殘差及尺度參量對閾值大小進(jìn)行估計(jì),所以隨機(jī)噪聲的閾值僅隨尺度的不同而發(fā)生改變,實(shí)現(xiàn)了閾值隨尺度變化的自適應(yīng)。
對一單炮地震數(shù)據(jù)進(jìn)行Shearlet變換,分別求取不同尺度角度下地震數(shù)據(jù)的信噪比和L2范數(shù),如圖1和2所示,分別展示了第1、2、3、4尺度下各個(gè)角度的信噪比和L2范數(shù)柱狀圖。通過對柱狀圖分析可知,在不同的尺度角度域中,數(shù)據(jù)的信噪比、L2范數(shù)相差很大。
單一的利用尺度參量進(jìn)行隨機(jī)噪聲閾值估計(jì)不能最大程度地利用Shearlet的稀疏特性,同時(shí)也不能達(dá)到最優(yōu)的去噪效果。地震數(shù)據(jù)同相軸具有方向性,而隨機(jī)噪聲比較均一,在Shearlet變換的某些方向上,地震數(shù)據(jù)的分布將會(huì)很少,這時(shí)Shearlet系數(shù)絕大部分表示隨機(jī)噪聲,用此部分Shearlet系數(shù)就能估計(jì)出某一尺度下不同角度隨機(jī)噪聲的能值范圍,即將有效信號(hào)分布最小的尺度和方向上分布的噪聲作為該尺度的噪聲,并將其剔除,得到分布在各尺度和方向上的有效地震數(shù)據(jù)。隨尺度和方向自適應(yīng)變化的閾值函數(shù)如下:
圖2 Shearlet域不同尺度角度下地震數(shù)據(jù)的L2范數(shù)Fig.2 L2 norm of seismic data at different scale angles in Shearlet domain
(12)
式中:ej,l為尺度j、方向l上的稀疏變換系數(shù)矩陣的L2范數(shù),j=1,2,…,J,J為分解總尺度;σ為噪聲標(biāo)準(zhǔn)差; min(·)為最小值函數(shù);λj與信噪比相關(guān),
(13)
SNR(j)為不同尺度Shearlet系數(shù)單獨(dú)構(gòu)成的地震信噪比,由相鄰地震道互相關(guān)求得[20]:
式中:Qi,i+1為第i道和第i+1道互相關(guān)函數(shù)最大值;Qi,i(0)為第i道自相關(guān)函數(shù)最大值;N為道數(shù)。
自適應(yīng)閾值函數(shù)將信噪比、隨機(jī)噪聲殘差、Shearlet系數(shù)L2范數(shù)作為閾值設(shè)定的因素,即不同尺度、角度的閾值估計(jì)權(quán)值系數(shù)不同,可以自適應(yīng)求取相應(yīng)的閾值,以達(dá)到最優(yōu)的隨機(jī)噪聲壓制效果,避免有效信號(hào)的損傷。
為驗(yàn)證基于Shearlet稀疏變換自適應(yīng)閾值法對隨機(jī)噪聲壓制的有效性,截取部分Sigsbee速度模型(圖3a)采用有線差分進(jìn)行正演模擬,采用道距20 m、炮距20 m,采樣間隔1 ms,8 s接收,加入隨機(jī)噪聲的單炮記錄如圖3b所示。隨機(jī)噪聲對深層弱反射信號(hào)干擾很大,部分微弱的反射信號(hào)在隨機(jī)噪聲的背景下已經(jīng)很難用肉眼識(shí)別。采用傳統(tǒng)Shearlet域閾值法和Shearlet域尺度角度自適應(yīng)閾值法對此單炮記錄進(jìn)行隨機(jī)噪聲壓制,結(jié)果如圖4。對比圖3和圖4可以看出,基于Shearlet變換的傳統(tǒng)閾值方法和本文提出的尺度角度自適應(yīng)閾值方法都能夠起到壓制單炮記錄中隨機(jī)噪聲的作用,主要的反射同相軸成像效果得到了提升。但是,圖4a的局部弱反射同相軸并沒有得到提升,有些本該存在的同相軸受到了損傷,相應(yīng)的圖4c中隨機(jī)噪聲波場含有有效反射波同相軸信息。圖4b的隨機(jī)噪聲去除效果明顯優(yōu)于圖4a,并且通過觀察去除的隨機(jī)噪聲波場(圖4d)未發(fā)現(xiàn)有效反射同相軸信息,壓制隨機(jī)噪聲的同時(shí),有效保護(hù)了中深層的弱反射信息。圖4b中構(gòu)造邊緣的弱反射同相軸成像效果得到了提升,單炮記錄上對構(gòu)造細(xì)節(jié)的呈現(xiàn)優(yōu)于圖4a,表明本文提出的Shearlet域尺度角度自適應(yīng)閾值方法在隨機(jī)噪聲壓制和弱信號(hào)保護(hù)上是有效的。圖4e、f分別為傳統(tǒng)閾值方法去噪的第三、第四尺度結(jié)果,圖4g、h分別為本文方法去噪的第三、第四尺度結(jié)果,在單獨(dú)的尺度域中,本文提出的Shearlet域尺度角度自適應(yīng)閾值方法在隨機(jī)噪聲壓制和弱信號(hào)保護(hù)上也是優(yōu)于傳統(tǒng)Shearlet域閾值方法的。
圖3 Sigsbee模型(a)及包含隨機(jī)噪聲的單炮記錄(b)Fig.3 Sigsbee model(a) and single shot record containing random noise(b)
a—傳統(tǒng)Shearlet域閾值方法;b—Shearlet域尺度角度自適應(yīng)閾值方法;c—傳統(tǒng)Shearlet域閾值方法去除噪聲;d—Shearlet域尺度角度自適應(yīng)閾值方法去除的噪聲;e—傳統(tǒng)Shearlet域閾值方法第三尺度結(jié)果;f—傳統(tǒng)Shearlet域閾值方法第四尺度結(jié)果;g—Shearlet域尺度角度自適應(yīng)閾值方法第三尺度結(jié)果;h—Shearlet域尺度角度自適應(yīng)閾值方法第四尺度結(jié)果a—traditional Shearlet domain threshold method;b—Shearlet domain scale angle adaptive threshold method;c—random noise by traditional Shearlet domain threshold method;d—random noise by Shearlet domain scale angle adaptive threshold method;e—the 3rd scale result of traditional Shearlet domain threshold method;f—the 4th scale result of traditional Shearlet domain threshold method;g—the 3rd scale result of Shearlet domain scale angle adaptive threshold method;h—the 4th scale result of Shearlet domain scale angle adaptive threshold method圖4 單炮記錄隨機(jī)噪聲去除效果Fig.4 Shearlet domain threshold method for the removal of random noise
為驗(yàn)證本文提出的方法對疊后數(shù)據(jù)的適用性,抽取含隨機(jī)噪聲模擬數(shù)據(jù)的零偏移距數(shù)據(jù),如圖5,理論上零偏移距數(shù)據(jù)相當(dāng)于一次疊加的疊后剖面。分別采用傳統(tǒng)Shearlet域閾值法和本文提出的Shearlet域尺度角度自適應(yīng)閾值法對此零偏移距數(shù)據(jù)進(jìn)行隨機(jī)噪聲壓制,結(jié)果如圖6所示。對比圖6a和圖6b,兩種方法對隨機(jī)噪聲都起到了很好的去除效果,但是,在交互薄層、構(gòu)造邊緣及深部構(gòu)造上,圖6b保留的有效信息更全,更能適應(yīng)高精度高保真度的地震數(shù)據(jù)處理要求。通過在零偏移距地震數(shù)據(jù)上測試可知,Shearlet稀疏變換自適應(yīng)閾值法的隨機(jī)噪聲壓制方法同樣適用于疊后地震數(shù)據(jù),能夠取得理想的處理效果。
圖5 包含隨機(jī)噪聲的零偏移距數(shù)據(jù)Fig.5 Zero-offset record containing random noise
a—傳統(tǒng)Shearlet域閾值方法;b—Shearlet域尺度角度自適應(yīng)閾值方法;c—傳統(tǒng)Shearlet域閾值方法去除的噪聲;d—Shearlet域尺度角度自適應(yīng)閾值方法去除的噪聲a—traditional Shearlet domain threshold method;b—Shearlet domain scale angle adaptive threshold method;c—Random noise by traditional Shearlet domain threshold method;d—Random noise by Shearlet domain scale angle adaptive threshold method圖6 零偏移距數(shù)據(jù)隨機(jī)噪聲去除效果Fig.6 Shearlet domain threshold method for the removal of random noise
為探明松遼盆地基底構(gòu)造結(jié)構(gòu)及盆山構(gòu)造演化過程,本單位近幾年在松遼盆地部署了多條深反射地震測線。圖7a為松遼盆地西部的一條深反射疊加剖面,觀測系統(tǒng)為道距25 m、炮距100 m,960道接收,采樣間隔1 ms,25 s記錄。經(jīng)過常規(guī)地震數(shù)據(jù)處理流程處理后(包含靜校正、異常振幅壓制、面波壓制、先行干擾波壓制、補(bǔ)償、反褶積等),剖面能夠呈現(xiàn)出主要的地質(zhì)構(gòu)造走勢。但是受隨機(jī)噪聲的影響,剖面整體信噪比并不高,中淺層構(gòu)造細(xì)節(jié)刻畫不清晰,中深層莫霍面上下反射同相軸不連續(xù),給后期的解釋工作帶來一定干擾。
將本文提出的Shearlet域尺度角度自適應(yīng)閾值隨機(jī)噪聲壓制方法應(yīng)用到松遼盆地深反射地震數(shù)據(jù)中。圖7b為經(jīng)過Shearlet域自適應(yīng)閾值隨機(jī)噪聲壓制后的疊加剖面,對比圖7a和圖7b可以發(fā)現(xiàn),經(jīng)過隨機(jī)噪聲壓制后,剖面的成像質(zhì)量明顯提升,中深層的構(gòu)造細(xì)節(jié)和莫霍面的走勢都得到清晰呈現(xiàn)。圖7c為去除的隨機(jī)噪聲,可以發(fā)現(xiàn)本文方法在實(shí)際數(shù)據(jù)應(yīng)用中只是去除了背景場的隨機(jī)噪聲部分,未對中深層的弱有效反射信號(hào)帶來損傷。圖8a、c對應(yīng)圖7a中的黃色和藍(lán)色方框圈定的位置,圖8b、d對應(yīng)圖7b中的黃色和藍(lán)色方框圈定的位置,通過局部放大進(jìn)行對比可以看出,經(jīng)過Shearlet域尺度角度自適應(yīng)閾值隨機(jī)噪聲壓制方法處理后,微弱的地震同相軸和構(gòu)造邊緣位置都得到了清晰的成像,剖面的精細(xì)成像質(zhì)量得到了很大的提升。此外,通過圖9去噪前后頻譜對比分析發(fā)現(xiàn),去噪后10~65 Hz主頻帶能量沒有損耗,證明了該方法對地震剖面實(shí)現(xiàn)了保幅處理。
a—松遼盆地深反射地震數(shù)據(jù)疊加剖面;b—Shearlet域尺度角度自適應(yīng)閾值去噪后剖面;c—去除的隨機(jī)噪聲a—superimposed profile of deep reflection seismic data in Songliao Basin;b—Shearlet domain scale angle adaptive threshold denoised profiles;c—random noise圖7 松遼盆地地震數(shù)據(jù)隨機(jī)噪聲去除效果Fig.7 Random noise removal effect of field seismic data in Songliao Basin
a—圖7a中黃色方框放大;b—圖7b中黃色方框放大;c—圖7a中藍(lán)色方框放大;d—圖7b中藍(lán)色方框放大a—enlarge the yellow box area in figure 7a;b—enlarge the yellow box area in figure 7b;c—enlarge the blue box area in figure 7a;d—enlarge the blue box area in figure 7b圖8 圖7局部放大對比Fig.8 Zoom in of Fig.7
圖9 去噪前、后剖面頻譜對比Fig.9 Spectrum comparison of profile before and after denoising
本文提出的結(jié)合Shearlet尺度域角度域地震數(shù)據(jù)的信噪比、隨機(jī)噪聲殘差及L2范數(shù)作為閾值設(shè)定因素的Shearlet域自適應(yīng)閾值方法,充分考慮了不同尺度—角度內(nèi)地震數(shù)據(jù)分布具有方向性而隨機(jī)噪聲分布無規(guī)律性的特征,實(shí)現(xiàn)了最大限度地保護(hù)地震數(shù)據(jù)有效信號(hào),去除背景場的隨機(jī)噪聲。
深反射地震數(shù)據(jù)中深部反射信號(hào)能量較弱,受背景場的隨機(jī)噪聲干擾嚴(yán)重,利用本文提出方法有效去除隨機(jī)噪聲的同時(shí),對深層弱信號(hào)起到了保護(hù)作用,提高了深反射地震剖面信噪比,提升了剖面整體成像質(zhì)量。由此可以證明,Shearlet域自適應(yīng)閾值方法能夠有效地應(yīng)用到深反射地震勘探領(lǐng)域,具有廣泛的應(yīng)用前景。