基于Seislet域分?jǐn)?shù)階閾值去噪算法的地震資料去噪

張入化,黃建平,國運東,雍 鵬,劉定進

(1.中國石油大學(xué)地球科學(xué)與技術(shù)學(xué)院,山東青島266580;2.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

地震資料去噪是地震資料處理中重要環(huán)節(jié)之一,而變換域去噪則是地震資料去噪的主要方法之一[1]。變換域去噪是通過某一數(shù)學(xué)變換方法將地震資料變換至對應(yīng)的變換域中,并根據(jù)有效信號和噪聲在變換域中的特點,去掉噪聲并保留有效信號。頻率域去噪[2]是最常用的去噪方法之一,但該方法無法壓制與有效信號頻帶范圍重合的噪聲。小波變換去噪[3]能夠?qū)⑿盘栠M行多尺度分解后再去噪,但該方法只能沿單一方向進行信號變換,無法適應(yīng)同時在多個方向上變化的信號。脊波變換、曲波變換[4-5]、Shearlet變換[6]、S變換[7]等都是小波變換的一種衍生方式,且都能夠處理在多個方向上變化的信號,因此在地震資料去噪中得到了廣泛應(yīng)用?；跁r頻分析的算法,如結(jié)合經(jīng)驗?zāi)B(tài)分解和獨立分量分析等[8-10],在實際地震資料處理中也取得了不錯的應(yīng)用效果,但這些算法都是基于圖像處理發(fā)展起來的,并不針對地震資料。

為了得到更加適用于地震數(shù)據(jù)的變換方法,FOMEL等[11-12]提出了一種全新的數(shù)學(xué)變換,即Seislet變換。Seislet變換也屬于小波變換類,但Seislet變換是沿著地震資料同相軸的方向進行變換處理,而沿測線方向的小波變換可以看作是零傾角的Seislet變換。小波變換、曲波變換等主要以信號采樣點為基本變換單位,而Seislet變換則以地震道為變換單位,因此Seislet變換能更好地適應(yīng)地震數(shù)據(jù),且能進一步應(yīng)用于地震資料去噪。劉洋等[13]將高階Seislet變換應(yīng)用于地震數(shù)據(jù)隨機噪聲處理,使隨機噪聲得到了有效壓制。勾福巖等[14]利用Seislet變換壓制海洋地震資料的涌浪噪聲,取得了較好的應(yīng)用效果。為獲得更好的去噪效果,劉財?shù)萚15]首次將閾值函數(shù)引入到Seislet變換壓制地震資料噪聲,但閾值函數(shù)本身會影響地震資料去噪效果。

1995年,DONOHO[16]最早提出了小波閾值去噪算法,介紹了硬閾值函數(shù)和軟閾值函數(shù)。硬閾值函數(shù)在閾值處存在斷點,軟閾值函數(shù)則存在恒定的偏差[17]。針對這一缺陷,一些學(xué)者提出了新的閾值函數(shù)。GAO等[18]利用半軟閾值函數(shù)改進軟閾值函數(shù),并取得了一定的效果。許麗群[19]進一步提出了硬軟閾值函數(shù),能夠一定程度地改善硬、軟閾值函數(shù)的缺陷。何柯等[20]在常規(guī)閾值函數(shù)的基礎(chǔ)上提出一種補償閾值函數(shù)算法,對之前受到損傷的與噪聲曲波系數(shù)重疊的弱信號進行補償。余江奇等[21]根據(jù)數(shù)據(jù)噪聲在曲波域的特點,采用了一種新的閾值函數(shù),并將其引入到曲波變換稀疏反褶積中。張華等[22]則在曲波域中的每一個尺度選取不同的閾值因子,達到曲波域中局部閾值去噪的目的。陳浩等[23]針對Shearlet變換閾值去噪方法不能隨尺度和方向變化的不足,提出了隨尺度和方向變化的自適應(yīng)閾值。

前人研究表明,傳統(tǒng)閾值去噪算法不能有效地分離差異較小的有效信號與噪聲,且傳統(tǒng)閾值確定準(zhǔn)則也不能很好地適用于Seislet域。本文在前人研究的基礎(chǔ)上,應(yīng)用Riemann-Liouville分?jǐn)?shù)階積分理論推導(dǎo)出分?jǐn)?shù)階閾值函數(shù);再根據(jù)地震數(shù)據(jù)在Seislet域不同尺度的分布特點,提出基于尺度加權(quán)的閾值確定方法;最后利用合成加噪模型數(shù)據(jù)和實際地震數(shù)據(jù)驗證了本文方法的可行性和有效性。

1 方法原理

1.1 Seislet變換基本原理

Seislet變換主要是利用小波提升算法[24](見附錄A),并應(yīng)用平面波解構(gòu)濾波器(PWD)[25],沿著地震同相軸的方向?qū)Φ卣饠?shù)據(jù)進行分解。

PWD-Seislet變換定義的預(yù)測算子P和更新算子U如下:

(1)

(2)

Seislet變換不再像小波變換沿著單一角度進行分解,而是根據(jù)地震同相軸的方向沿著整個地震數(shù)據(jù)空間進行分解,能更為精細地刻畫地下地質(zhì)特征。因此有效信號在Seislet域中能被很好地壓縮,其Seislet系數(shù)集中在較低的尺度。

1.2 分?jǐn)?shù)階閾值算法

閾值去噪算法最早應(yīng)用于小波變換去噪,現(xiàn)在廣泛應(yīng)用于曲波變換去噪、Shearlet變換去噪等方法中。稀疏變換將信號投影到對應(yīng)的變換域,有效信號在合適的變換域下具有稀疏特征,即有效信號的能量集中在有限的變換域系數(shù)范圍,而噪聲則遍布整個變換空間。一般認為代表有效信號的變換系數(shù)幅值大于代表噪聲的變換系數(shù)幅值,在確定了閾值之后,大于閾值的變換系數(shù)被認為是由有效信號變換而來的,小于閾值的變換系數(shù)則被認為是由噪聲變換而來的,這時可以采用閾值函數(shù)對變換系數(shù)進行處理,最后將處理后的系數(shù)進行反變換,達到壓制噪聲、保留有效信號的目的。

傳統(tǒng)的閾值函數(shù)主要有硬閾值函數(shù)和軟閾值函數(shù)。

硬閾值函數(shù)是將系數(shù)絕對值小于閾值的系數(shù)直接置為0,將系數(shù)絕對值大于等于閾值的系數(shù)保持原值,但在閾值處存在斷點,如圖1a所示。硬閾值函數(shù)如下:

(3)

式中:s為Seislet稀疏變換后的系數(shù);λ為閾值。

軟閾值函數(shù)將系數(shù)絕對值大于閾值的系數(shù)等于閾值與該系數(shù)的差值,這樣就能確保閾值函數(shù)的連續(xù)性,但軟閾值函數(shù)處理后的系數(shù)存在較大的偏差,重構(gòu)信號與原信號差異較大,如圖1b所示。軟閾值函

數(shù)如下:

(4)

式中:sgn(·)為符號函數(shù)。

硬閾值函數(shù)在閾值附近存在斷點,重構(gòu)信號會產(chǎn)生偽吉布斯現(xiàn)象,一般較少采用。常用的主要是軟閾值函數(shù),但軟閾值函數(shù)處理后的系數(shù)與原系數(shù)相比具有恒定的偏差。近幾年來,分?jǐn)?shù)階微積分理論在科學(xué)和工程學(xué)的很多領(lǐng)域得到了應(yīng)用。本文引入Riemann-Liouville分?jǐn)?shù)階積分理論[25-27],其公式如下:

(5)

將分?jǐn)?shù)階積分理論應(yīng)用到閾值函數(shù)中,推導(dǎo)得到分?jǐn)?shù)階閾值函數(shù):

圖1 不同閾值函數(shù)處理前、后系數(shù)對比a 硬閾值; b 軟閾值; c 分?jǐn)?shù)階閾值(設(shè)閾值為2)

(6)

式中:f(s)為某一載體函數(shù),此處用硬閾值函數(shù)作為載體函數(shù);s為Seislet變換系數(shù)。經(jīng)過多次測試,本文選用的分?jǐn)?shù)階階數(shù)α為0.5階。被該分?jǐn)?shù)階積分公式積分后的閾值函數(shù)T(s)既能增強函數(shù)本身的連續(xù)性,又能保留原函數(shù)的特征。公式(6)的數(shù)學(xué)意義相當(dāng)于對硬閾值函數(shù)做0.5階的分?jǐn)?shù)階積分處理,在保留硬閾值函數(shù)原有特征的基礎(chǔ)上,進一步增強其在閾值處的連續(xù)性,進而得到分?jǐn)?shù)階閾值函數(shù)。與傳統(tǒng)的硬閾值函數(shù)和軟閾值函數(shù)相比,分?jǐn)?shù)階閾值函數(shù)很好地控制了系數(shù)偏差,在端點處也得到了一定的平滑,如圖1c所示。

閾值是有效信號與噪聲存在區(qū)別的臨界點。閾值越小,噪聲壓制的程度越小,閾值越大,噪聲壓制的程度越大,所以閾值的選擇會直接影響閾值去噪算法的效果[16]。噪聲的Seislet系數(shù)與有效信號的Seislet系數(shù)在Seislet域的不同尺度不同,所以需要對每一個尺度的噪聲Seislet系數(shù)和有效信號Seislet系數(shù)進行分析,以此計算閾值。本文根據(jù)地震數(shù)據(jù)在Seislet域的分布特點,并在DONOHO給出的全局閾值[28]基礎(chǔ)上提出了一種尺度加權(quán)閾值,尺度越小,權(quán)重越大,尺度越大,權(quán)重越小,且得到的閾值在每一個尺度上都不同,公式如下:

(7)

式中:σ是噪聲方差估計公式[29];max(·)表示取最大值;max(j)-j為權(quán)重;Nj表示尺度的層數(shù);Ni表示每一層系數(shù)的個數(shù);j為尺度;i為尺度上每個系數(shù)的序號;w為調(diào)節(jié)參數(shù),最大不超過分解的最大尺度,w越大對噪聲壓制的程度越大;sj,i為不同尺度上的Seislet系數(shù);λj為不同尺度上的閾值;median(·)表示取中值。

1.3 Seislet分?jǐn)?shù)階閾值算法流程

Seislet分?jǐn)?shù)階閾值算法流程見圖2。該流程使用平面波解構(gòu)濾波器(PWD)來估計地震數(shù)據(jù)的局部傾角,也可以使用線性Randon變換、預(yù)測誤差濾波器、有限差分法等方法來估計地震數(shù)據(jù)的局部傾角[30]。在該流程中,尺度加權(quán)閾值的設(shè)置非常重要,第一次設(shè)置的閾值w為最大尺度的一半,再根據(jù)去噪效果的好壞進行調(diào)節(jié),直到達到最好的去噪效果。尺度加權(quán)閾值利用Seislet域低尺度中有效信號分量遠多于高尺度中有效信號分量的特點,在Seislet域中設(shè)置更為合理的閾值,也可以根據(jù)去噪效果人為調(diào)控參數(shù)w,得到更好的去噪效果。去噪效果的好壞由人為抽取的單道振幅頻率、去噪結(jié)果剖面以及壓制的噪聲進行綜合分析得出。

圖2 Seislet分?jǐn)?shù)階閾值算法流程

2 數(shù)值測試

2.1 反射波雙曲線模型測試

為了驗證本文提出的Seislet分?jǐn)?shù)階閾值算法的可行性和有效性,首先利用共炮點反射波時距曲線公式直接模擬得到反射波雙曲線模型。地震反射波時距曲線公式如下:

(8)

式中:x是炮檢距;v為地下介質(zhì)速度;h為深度。

該模型共有2條反射曲線,如圖3a所示,道間距為10m,共計512道,時間采樣間隔為2ms,采樣時間為1.024s。

利用公式(9)計算信噪比,將添加噪聲的反射波雙曲線模型(圖3b)信噪比RSN控制在0.25dB。

(9)

式中:Si,j表示有效信號第i道第j個采樣點的振幅值;Ni,j表示添加的隨機噪聲第i道第j個采樣點的振幅值,計算結(jié)果的單位為dB。

進行Seislet變換之前,需要獲得局部傾角估計數(shù)據(jù)。本文使用平面波解構(gòu)濾波器(PWD)對含噪模型(圖3b)進行局部傾角估計,得到的局部傾角估計值如圖3c所示,由圖3可知,因為噪聲的幅值較大,對局部傾角的估計造成一定的影響,能量團的分布較為散亂,但在局部傾角估計圖中也能夠識別出反射波曲線。圖3d為含噪模型的Seislet系數(shù)分布,在尺度小于50的范圍內(nèi)系數(shù)值比較集中,主要為有效信號的Seislet系數(shù),但在整個尺度范圍都有Seislet系數(shù)的分布,說明隨機噪聲的Seislet系數(shù)在Seislet域中仍然分布雜亂,難以集中。在模型測試中,分?jǐn)?shù)階閾值算法和傳統(tǒng)硬、軟閾值算法都采用全局閾值。

對含噪模型進行Seislet變換后,分別采用硬閾值算法、軟閾值算法、分?jǐn)?shù)階閾值算法去噪,得到不同的去噪結(jié)果(圖4)。圖4a為經(jīng)過硬閾值去噪后的結(jié)果,噪聲得到部分衰減,但還存在部分噪聲。經(jīng)過軟閾值算法去噪后(圖4b),大部分噪聲得到了壓制,整體信噪比高于使用硬閾值算法去噪的結(jié)果。采用分?jǐn)?shù)階閾值算法的去噪結(jié)果中(圖4c),噪聲幾乎全被壓制,去噪效果明顯優(yōu)于前2種算法。因為局部傾角估計算法受到噪聲影響導(dǎo)致估計結(jié)果不準(zhǔn),所以圖4c的去噪結(jié)果存在部分扭曲現(xiàn)象。采用公式(9)計算去噪后數(shù)據(jù)的信噪比,并以其衡量3種算法的去噪效果發(fā)現(xiàn),硬閾值和軟閾值算法去噪結(jié)果的信噪比分別為1.72dB和2.13dB,而分?jǐn)?shù)階閾值算法去噪結(jié)果的信噪比為3.82dB,去噪后的信噪比明顯高于傳統(tǒng)硬、軟閾值算法。

為了進一步驗證方法的有效性,對3種算法去噪后的結(jié)果、加噪模型數(shù)據(jù)以及未加噪模型數(shù)據(jù)分別抽取第256道進行頻譜分析,結(jié)果如圖5所示,可以看出,未加噪數(shù)據(jù)(圖5a)的主頻在50Hz左右,幅值由中間向兩邊降低,呈現(xiàn)近似的正態(tài)分布;加噪數(shù)據(jù)(圖5b)的頻率在整個頻帶范圍都呈現(xiàn)出一種無規(guī)律的分布,高頻率分量較多;雖然圖5c、圖5d、圖5e的高頻分量都有所減少,但圖5e的高頻分量最少,且其主頻率成分與圖5a的最為接近。綜合分析圖4和圖5可知,傳統(tǒng)硬、軟閾值算法對隨機噪聲的壓制效果不甚理想。

2.2 實際地震資料測試

使用某陸上工區(qū)的實際單炮地震記錄(圖6a)進行測試。該單炮記錄的道間距為10m,共128道,采樣間隔為1ms,采樣時間1.4s。由圖6a可知,該單炮記錄淺部能量和初至能量較強,底部能量較弱但存在異常振幅,整個單炮記錄上存在一定的隨機噪聲和線性噪聲。圖6b為該單炮記錄的局部傾角估計,由圖可見,能量團的分布層次分明,如果某塊區(qū)域能量團顏色較為一致,則說明這塊區(qū)域的局部傾角值大致相同。利用圖6b所示局部傾角估計結(jié)果進行Seislet變換得到Seislet系數(shù)(圖6c)。在實際地震資料測試中,首先在分?jǐn)?shù)階閾值算法、硬閾值算法、軟閾值算法中采用尺度加權(quán)閾值來進行地震資料去噪,然后再在分?jǐn)?shù)階閾值算法中分別采用尺度加權(quán)閾值和全局閾值進行去噪效果對比。

圖4 不同算法的去噪結(jié)果對比a 硬閾值算法; b 軟閾值算法; c 分?jǐn)?shù)階閾值算法

圖5 加噪前、后模型數(shù)據(jù)及其采用不同閾值法去噪后單道記錄頻譜對比a 未加噪模型; b 加噪模型; c 硬閾值算法; d 軟閾值算法; e 分?jǐn)?shù)階閾值算法

對實際單炮記錄進行Seislet變換后,再進行硬閾值算法、軟閾值算法、分?jǐn)?shù)階閾值算法進行去噪,結(jié)果見圖7。經(jīng)過硬閾值算法去噪后的結(jié)果(圖7a)能夠在一定程度上壓制隨機噪聲,但整體去噪程度不如圖7b和圖7c。軟閾值算法(圖7b)和分?jǐn)?shù)階閾值算法(圖7c)都能較好地壓制噪聲,在局部區(qū)域的同相軸連續(xù)性和噪聲壓制方面,分?jǐn)?shù)階閾值去噪算法的效果優(yōu)于軟閾值函數(shù)。將未去噪的單炮數(shù)據(jù)分別和3種不同方法去噪后的結(jié)果相減,得到的噪聲記錄如圖8 所示。硬閾值算法去除的噪聲相對較少,且未能壓制一些振幅較強的噪聲(黃色箭頭);軟閾值算法壓制的噪聲中卻去掉了一定的有效信號(紅色圈和紅色箭頭),且軟閾值算法去掉了部分初至信息;分?jǐn)?shù)階閾值算法的去噪效果比前兩種閾值算法更好,對有效信號的傷害最輕。在分?jǐn)?shù)階閾值算法去噪中分別采用尺度加權(quán)閾值和全局閾值,并對去噪結(jié)果和壓制的噪聲作對比分析(圖9)發(fā)現(xiàn),采用全局閾值的分?jǐn)?shù)階閾值去噪算法未能壓制底部的異常振幅(黃色箭頭),且同相軸附近仍存在部分噪聲(紅色箭頭),而采用尺度加權(quán)閾值的分?jǐn)?shù)階閾值去噪算法能更有效地壓制噪聲。

圖6 實際單炮記錄(a)、局部傾角估計(b)及Seislet系數(shù)分布(c)

圖7 采用不同閾值算法去噪后的單炮記錄a 硬閾值算法; b 軟閾值算法; c 分?jǐn)?shù)階閾值算法

圖8 對實際單炮記錄采用不同閾值方法去除的噪聲對比a 硬閾值; b 軟閾值; c 分?jǐn)?shù)階閾值

圖9 采用不同閾值確定準(zhǔn)則去噪后的單炮記錄(左)對比及去除的噪聲(右)對比a 采用尺度加權(quán)閾值; b 全局閾值

為進一步驗證分?jǐn)?shù)階閾值去噪算法的有效性,分別對實際單炮記錄、硬閾值算法、軟閾值算法、分?jǐn)?shù)階閾值算法的去噪結(jié)果抽取第64道數(shù)據(jù)進行頻譜分析,結(jié)果如圖10所示。圖10a為實際單炮記錄的第64道頻譜,可見存在一定的高頻分量,但幅值不高,頻率抖動次數(shù)較多,這些高頻分量為對應(yīng)的背景噪聲,主頻應(yīng)在50～80Hz。硬閾值算法處理后(圖10b),頻譜的高頻分量(大于50Hz)并沒有太過明顯的改觀。相較于硬閾值算法,軟閾值算法處理結(jié)果的單道頻譜在50Hz和80Hz處的幅值均有所下降,且趨于平緩(圖10c)。結(jié)合硬、軟閾值算法去噪的結(jié)果(圖7a 和圖7b)可見,軟閾值算法的去噪效果優(yōu)于硬閾值算法。經(jīng)過分?jǐn)?shù)階閾值算法去噪后的頻譜(圖10d) 的高頻分量幅值低于軟閾值算法去噪后的結(jié)果,且頻率變化也較為緩慢(見紅色箭頭),再結(jié)合不同閾值算法的去噪結(jié)果(圖7c)可見,分?jǐn)?shù)階閾值去噪算法能夠有效地壓制背景噪聲,且效果優(yōu)于硬閾值去噪算法和軟閾值去噪算法。

圖10 實際單炮記錄第64道頻譜分析a 未去噪單炮; b 硬閾值算法去噪; c 軟閾值算法去噪; d 分?jǐn)?shù)階閾值算法去噪

3 結(jié)論

本文在Riemann-Liouville分?jǐn)?shù)階積分公式的基礎(chǔ)上得到分?jǐn)?shù)階閾值函數(shù),增強了閾值處的連續(xù)性,進而提出一種用于Seislet域去噪的尺度加權(quán)閾值函數(shù),實現(xiàn)了閾值計算的局部精細化,應(yīng)用了尺度加權(quán)閾值的分?jǐn)?shù)階閾值算法能有效壓制噪聲并降低有效信號的損失。分別采用硬閾值算法、軟閾值算法、分?jǐn)?shù)階閾值算法對理論模型數(shù)據(jù)和實際地震資料進行處理,通過對去噪效果、信噪比數(shù)據(jù)以及頻譜對比分析,結(jié)果表明:傳統(tǒng)硬、軟閾值算法噪聲壓制不徹底,且不能有效分離差異較小的有效信號和噪聲;相較于傳統(tǒng)閾值算法,分?jǐn)?shù)階閾值算法能更好地壓制噪聲,減少對有效信號的損傷;與全局閾值相比,本文提出的尺度加權(quán)閾值在Seislet域中具有更好的適應(yīng)性,在閾值去噪中效果更佳。

雖然Seislet變換能夠?qū)^為復(fù)雜的地震數(shù)據(jù)進行變換處理,但Seislet變換較為依賴局部傾角的估計結(jié)果,一旦信噪比過低,便會影響Seislet變換的準(zhǔn)確性。如何提高局部傾角的估計精度或者減少Seislet變換對局部傾角數(shù)據(jù)的依賴是今后的改進方向。