樂友喜 楊 濤 曾賢德
(①中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島 266580;②黃河勘測規(guī)劃設(shè)計(jì)研究院有限公司,河南鄭州 450003;③中國能源建設(shè)集團(tuán)新疆電力設(shè)計(jì)院有限公司,新疆烏魯木齊 830001)
地震信號(hào)中混入隨機(jī)噪聲時(shí),會(huì)大大降低地震資料的信噪比,嚴(yán)重影響對有效信號(hào)的檢測和準(zhǔn)確估計(jì)。研究表明,由于在各個(gè)時(shí)刻具有不確定性,決定了隨機(jī)噪聲呈現(xiàn)在零點(diǎn)處自相關(guān)函數(shù)值最大,而在其他點(diǎn)處自相關(guān)函數(shù)值迅速衰減的特點(diǎn)。對于理想高斯白噪聲,其自相關(guān)函數(shù)值在零點(diǎn)處應(yīng)表現(xiàn)為一個(gè)尖脈沖;而對于其他一般信號(hào),由于存在固有的關(guān)聯(lián)度,其自相關(guān)函數(shù)值的變化規(guī)律明顯不同于隨機(jī)噪聲[1]。
現(xiàn)有的地震信號(hào)去噪方法已相當(dāng)成熟,主要包括Wiener濾波[2]、譜減法[3]、小波變換[4-5]等傳統(tǒng)的(去噪)算法及其改進(jìn)方法。
KSVD(K Singular Value Decomposition)字典稀疏去噪法屬于基本壓縮感知理論[6-7]的改進(jìn)算法,該方法通過大量樣本信號(hào)訓(xùn)練獲得過完備冗余字典,能更自適應(yīng)地以稀疏基表征[8-9]信號(hào),從而更精確地重構(gòu)信號(hào),實(shí)現(xiàn)去噪目的。盡管KSVD字典[10-11]稀疏去噪可使信號(hào)中大部分隨機(jī)噪聲得到有效抑制或消除,但在信號(hào)的各個(gè)波峰和波谷處總會(huì)殘留一些毛刺[12]。
經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposi- tion,EMD)可對信號(hào)局部時(shí)變特征做自適應(yīng)時(shí)頻分解[13-14],非常適用于非平穩(wěn)、非線性信號(hào)分析,但它存在模態(tài)混疊和端點(diǎn)效應(yīng)問題[15]。
完備總體經(jīng)驗(yàn)?zāi)B(tài)分解(Complete Ensemble Empirical Mode Decomposition,CEEMD)[16]彌補(bǔ)了EMD存在的缺陷,在分解信號(hào)的每一階段添加一種特定白噪聲,同時(shí)計(jì)算唯一的殘差以獲取每個(gè)符合定義的固有模態(tài)分量(IMF),這是一個(gè)完備的分解過程。CEEMD提供的模態(tài)分離效果明顯好于EMD,同時(shí)對原始信號(hào)的重構(gòu)更為精確,且計(jì)算效率相對更高[17]。
本文針對傳統(tǒng)地震信號(hào)的隨機(jī)噪聲去除方法的不足及KSVD字典稀疏去噪后信號(hào)中總會(huì)殘留毛刺的問題,提出一種CEEMD與KSVD字典訓(xùn)練相結(jié)合的去噪方法,以期能夠在實(shí)際地震資料中更好地消除隨機(jī)噪聲的影響。
總體經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition,EEMD)在一定程度上解決了EMD的模態(tài)混疊效應(yīng),提高了EMD的分解質(zhì)量;同時(shí)也造成了新問題,如加入不同高斯白噪聲會(huì)導(dǎo)致分解出的IMF不是唯一解,與IMF的原始定義不能很好吻合,這樣就失去它所代表的物理含義。EEMD算法的完備性不強(qiáng),對處理后結(jié)果做累加難以完美重構(gòu)原來信號(hào);同時(shí)受到加總平均[12]的影響,在很大程度上,噪聲添加越多分解結(jié)果越好,但其計(jì)算成本也會(huì)相應(yīng)地增加。
Torres等[18]于2011年提出了一種噪聲自動(dòng)適應(yīng)的完備總體經(jīng)驗(yàn)?zāi)B(tài)分解(CEEMD),該算法與EEMD相同點(diǎn)在于也利用了噪聲進(jìn)行輔助分析,但具體實(shí)現(xiàn)過程與EEMD區(qū)別很大。此算法不僅解決了模態(tài)混疊問題,同時(shí)提高了原信號(hào)重構(gòu)的精度、提升了計(jì)算效率也降低了成本,總體上明顯優(yōu)于EMD、EEMD方法。其具體步驟描述如下。
(1)將不同的白噪聲分別與原來信號(hào)合成T個(gè)混合信號(hào),應(yīng)用EMD對其處理,計(jì)算集合平均并將其作為原來信號(hào)的第一個(gè)固有模態(tài)函數(shù)IMF1
(1)
式中:Fj(·)為EMD處理后獲取的第j階模態(tài);ωi為i個(gè)高斯白噪聲;εk為每一階段白噪聲加入的比例大??;x(t)為初始信號(hào)。
(2)令r0(t)=x(t),對k=1、…、K計(jì)算第k階殘差rk(t)
rk(t)=rk-1(t)-IMFk(t)
(2)
(3)對rk(t)+εkFk[ωi(t)]做EMD處理,獲取對應(yīng)的IMF1;計(jì)算總體平均并將其作為IMFk+1
(3)
(4)重復(fù)步驟(2)和步驟(3),直到殘差信號(hào)不能分解為止,得到最終殘差
(4)
KSVD算法[19]主要包括字典更新和稀疏編碼兩個(gè)步驟,獲取字典的同時(shí)也獲得稀疏系數(shù)矩陣。
在稀疏編碼過程中,常用的稀疏分解方法為正交匹配追蹤算法(Orthogonal Matching Pursuit,OMP)。字典更新的常用方法有MOD(Method of Optimal Directions)和KSVD,兩者差異之處在于字典更新原子時(shí)所用更新機(jī)制不同: MOD是對字典中所有原子進(jìn)行更新,而KSVD只針對字典中某個(gè)原子進(jìn)行。
在字典[20]更新過程中,KSVD學(xué)習(xí)字典使用了奇異值分解更新字典原子,與固定基字典相比,訓(xùn)練的字典具有更能反映信號(hào)特征的優(yōu)點(diǎn)。KSVD算法具體流程(圖1)如下。
圖1 KSVD稀疏去噪方法框圖
參數(shù)初始化: 給定訓(xùn)練樣本X=[x1,x2,…,xN], 初始超完備字典D, 正則化參數(shù)λ和迭代次數(shù)J。
(5)
(3)達(dá)到迭代次數(shù)J即結(jié)束,否則返回步驟(1)。
經(jīng)過上述步驟(1)~步驟(3),即可得到學(xué)習(xí)后的KSVD字典,在整個(gè)KSVD運(yùn)算中,每次字典更新時(shí)都能保證總誤差單調(diào)下降,因此該算法具有較好收斂特性。此外,KSVD算法還具有一定的去噪功能。
CEEMD去噪方法盲目舍棄高頻IMF模態(tài)分量有可能導(dǎo)致高頻有效信號(hào)損失或高頻噪聲壓制不完全。KSVD字典稀疏去噪盡管去除了大部分噪聲,但弱隨機(jī)噪聲尚存在于局部信號(hào)中,尤其在信號(hào)各個(gè)波峰和波谷處尚殘留一些較明顯毛刺。鑒于此,本文將CEEMD與KSVD相結(jié)合,通過互補(bǔ)方式彌補(bǔ)上述方法存在的兩點(diǎn)不足。
CEEMD-KSVD聯(lián)合去噪具體實(shí)現(xiàn)步驟(圖2,其中“過渡模態(tài)”是指介于噪聲主導(dǎo)模態(tài)與信號(hào)主導(dǎo)模態(tài)之間的模態(tài))如下:
(1)采用CEEMD對觀測含噪信號(hào)y模態(tài)分解獲得本征模態(tài)函數(shù)集。
(2)求取每個(gè)IMF分量的峰度,識(shí)別噪聲主導(dǎo)模態(tài),過渡模態(tài)和信號(hào)主導(dǎo)模態(tài)。
(3)噪聲主導(dǎo)模態(tài)分量:直接舍棄。
(4)過渡模態(tài)分量: ①對過渡模態(tài)分量累加合成新信號(hào)x,并對其二次CEEMD分解;②利用各個(gè)模態(tài)分量求取的峰度,舍棄噪聲主導(dǎo)模態(tài);③對步驟②中的剩余IMF分量累加合成新信號(hào)x1,并與觀測信號(hào)y作為訓(xùn)練樣本,獲得KSVD學(xué)習(xí)字典;④利用上步字典及OMP算法對信號(hào)x1稀疏去噪獲得去噪結(jié)果y1。
(5)信號(hào)主導(dǎo)模態(tài)分量: ①各個(gè)信號(hào)主導(dǎo)模態(tài)累加重構(gòu)得到新信號(hào)z;②將信號(hào)z和y作為訓(xùn)練樣本,獲取KSVD學(xué)習(xí)字典;③利用學(xué)習(xí)字典及OMP算法對信號(hào)z稀疏去噪得到y(tǒng)2。
圖2 本文去噪方法流程框圖
(6)將步驟(4)和步驟(5)中的去噪結(jié)果y1和y2重構(gòu)合并,得到最終去噪結(jié)果。
峰度又稱峰態(tài)系數(shù),是用于表征數(shù)據(jù)的概率密度分布曲線在平均值處峰值高低的特征數(shù)[21]。對于服從高斯分布的隨機(jī)噪聲信號(hào),其峰度F接近0;對于服從超高斯分布的隨機(jī)噪聲信號(hào),其峰度大于0;反之, 服從亞高斯分布的隨機(jī)噪聲信號(hào),其峰度小于0。
下文中峰度評判標(biāo)準(zhǔn): 噪聲主導(dǎo)模態(tài)是|F|≤0.1, 過渡模態(tài)是0.1<|F|≤1.0, 信號(hào)主導(dǎo)模態(tài)是|F|>1.0。
設(shè)計(jì)如圖3所示的褶積模型,該模型共有三層,上下兩層為水平層,中間層為傾斜層,并與上覆層斜交。橫向共有25道,縱向上有300個(gè)采樣點(diǎn),采樣頻率為500Hz,雷克子波主頻為25Hz,模型中加入兩種不同比例高斯白噪聲,分別采用五種方法對含噪模擬地震數(shù)據(jù)進(jìn)行去噪,通過分析驗(yàn)證五種方法去噪效果。
圖3 無噪模擬地震記錄
特別說明: 本文所使用的剖面數(shù)據(jù)的字塊為8×8,字典原子的維數(shù)是64×256;由于過渡模態(tài)部分合成信號(hào)的隨機(jī)噪聲較強(qiáng),使用的迭代次數(shù)為20,信號(hào)主導(dǎo)模態(tài)部分合成信號(hào)的隨機(jī)噪聲相對弱小,為提高計(jì)算效率,使用的迭代次數(shù)為10;同時(shí)分割子塊尺寸和字典維數(shù)不變時(shí),迭代次數(shù)越多,字典就會(huì)變得越稀疏,計(jì)算效率也就越高;反之迭代次數(shù)不變時(shí),子塊尺寸和字典維數(shù)越大,計(jì)算效率越低。本文加入的高斯白噪聲頻譜范圍分布廣泛,包含地震頻帶(圖4)。
加入較高比例的高斯白噪聲構(gòu)成信噪比為0.88的加噪模擬地震記錄(圖5a), 圖5b~圖5f為五種方法的去噪結(jié)果。從圖中可看到:F-X域去噪和CEEMD分頻去噪效果最差;小波軟閾值去噪效果也不太理想;KSVD稀疏去噪效果較好,但尚存一些毛刺現(xiàn)象;本文方法基本保存了原始信號(hào)的有效信息,去噪效果最好。表1為五種方法去噪結(jié)果評價(jià)指標(biāo)對比。從各方法去噪結(jié)果及其評價(jià)指標(biāo)的對比得知,五種方法中后兩種去噪效果較好,尤其本文方法相比其他方法去噪后信噪比最高,有效信息損失最少,是一種有效保幅的去噪方法。信噪比計(jì)算公式如下
圖4 高斯白噪聲(a)及其頻譜特征(b)
表1 五種方法去噪結(jié)果評價(jià)指標(biāo)對比
圖5 含較高比例噪聲模擬地震記錄五種方法去噪效果對比
(6)
式中:S1為去噪后的數(shù)據(jù);S為無噪數(shù)據(jù);RS/N為去噪后信噪比。
加入較低比例的高斯白噪聲構(gòu)成信噪比為3.13的加噪模擬地震記錄(圖6a),圖6b~圖6f為五種方法的去噪結(jié)果。從圖中可知: 五種方法都能將大部分噪聲壓制,但F-X域去噪和CEEMD分頻去噪有明顯的噪聲殘余;小波軟閾值去噪沒有毛刺現(xiàn)象,但在局部改變了信號(hào)的幅值;KSVD稀疏去噪效果較好,與本文方法去噪效果相差無幾,但是仔細(xì)觀察仍存在輕微的毛刺現(xiàn)象。表2為五種方法去噪結(jié)果評價(jià)指標(biāo)對比。因此,綜合對比五種方法的去噪結(jié)果及其評價(jià)指標(biāo)可看出:對高信噪比資料,后三種方法去噪效果均較好,尤其KSVD稀疏去噪與本文方法去噪效果只有些許差別。
圖6 含較低比例噪聲模擬地震記錄五種方法去噪效果對比
表2 五種方法去噪結(jié)果評價(jià)指標(biāo)對比
通過上述分析可得到以下認(rèn)識(shí): F-X域去噪和CEEMD分頻去噪不管是在高信噪比還是低信噪比資料中盡管能壓制一定程度的噪聲,但都未獲取理想去噪結(jié)果;小波軟閾值去噪和KSVD稀疏去噪效果質(zhì)量與信噪比呈線性關(guān)系,信噪比越高,效果越好;本文方法實(shí)用性強(qiáng),在不同信噪比資料中都能獲取較好的去噪效果。
為了進(jìn)一步測試本文去噪方法對實(shí)際地震資料的處理,將其應(yīng)用于實(shí)際地震資料。
圖7a為某一工區(qū)的疊后地震記錄,該地震記錄的CDP數(shù)目為250,采樣點(diǎn)個(gè)數(shù)為500,采樣時(shí)間為2ms,從圖中可見地震剖面上由于存在大量的隨機(jī)噪聲致使剖面同相軸不連續(xù)或錯(cuò)斷,剖面顯示不清晰,對地質(zhì)解釋帶來諸多不便。為了消除隨機(jī)噪聲的影響,分別利用上述五種方法進(jìn)行去噪處理,對比這五種方法去噪后的成果剖面(圖7b~圖7f)可知: F-X域去噪與CEEMD分頻去噪效果最差,去噪后還存在許多噪聲殘留;小波軟閾值去噪方法除掉了大量隨機(jī)噪聲,仍存有少許隨機(jī)噪聲,去噪效果一般;KSVD稀疏去噪與本文方法去噪效果相比前三種較好,但KSVD稀疏去噪在一些細(xì)節(jié)(弱幅度有效信號(hào))上不如本文方法,由圖7d和圖7e可知,本文方法充分壓制了隨機(jī)噪聲,同相軸相對更加連續(xù),有效信息得以保留,地震剖面清晰度提高。圖8為各種去噪方法獲得的殘差剖面圖。由殘差剖面可知,F(xiàn)-X域、KSVD稀疏、小波等去噪法都或多或少地去除了部分有效信號(hào),尤其F-X域去噪損失有效信號(hào)最嚴(yán)重;CEEMD分頻去噪,當(dāng)舍棄模態(tài)IMF1時(shí),去除的基本上是隨機(jī)噪聲,但是去除IMF1+IMF2時(shí),導(dǎo)致一部分有效信號(hào)的丟失,而本文方法去噪不僅隨機(jī)噪聲得以完全壓制,同時(shí)又不傷害有用信號(hào),取得了很好的去噪效果。
圖7 針對實(shí)際資料五種方法去噪效果對比
圖8 五類去噪方法對應(yīng)的殘差剖面
(1)本文研究表明: F-X域去噪方法能去除部分隨機(jī)噪聲,但會(huì)導(dǎo)致假同相軸、有效波波形畸變;由于閾值很難精準(zhǔn)選取,小波去噪方法在去除大量噪聲的同時(shí)也丟失了部分有效信號(hào);CEEMD去噪方法因盲目舍棄IMF分量可能造成高頻有效信息缺失;KSVD稀疏去噪總有殘留毛刺,且當(dāng)資料信噪比較低時(shí),去噪效果也隨之變差;本文提出的CEEMD-KSVD聯(lián)合去噪方法的去噪效果優(yōu)于其他四種方法,能顯著提高地震資料品質(zhì)。
(2)隨機(jī)噪聲的頻率大多較高,CEEMD分頻去噪方法通過舍棄前若干個(gè)模態(tài),主要用于去除高頻部分的隨機(jī)噪聲,KSVD稀疏去噪方法是以隨機(jī)噪聲與信號(hào)在整個(gè)定義域上稀疏分解系數(shù)分布差異較大為基礎(chǔ),兩種方法都立足于隨機(jī)噪聲的特性。因此,本文去噪方法更適用于隨機(jī)噪聲的壓制,對于其他類型噪聲,其適用性有待測試。