基于噪聲水平估計(jì)的加權(quán)核范數(shù)最小化噪聲壓制方法研究

唐 杰,張文征,戚瑞軒,李 聰

(中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島266580)

因各種干擾的存在,地震資料中會(huì)存在隨機(jī)噪聲,有些衰減噪聲的濾波處理方法會(huì)造成地震反射中的不連續(xù)性信息被平滑而變得模糊[1-2]。為了克服過度平滑的問題,可以根據(jù)地震數(shù)據(jù)中的噪聲水平采用自適應(yīng)的濾波處理方法。噪聲水平估計(jì)的常見方法有:基于小波分解的方法、基于主成分分析的方法和基于尺度不變性的方法[3]。基于小波分解的噪聲估計(jì)方法根據(jù)小波變換的特點(diǎn)來估計(jì)噪聲的標(biāo)準(zhǔn)方差,小波變換后數(shù)據(jù)的能量主要集中在尺度大的子帶,而尺度小的高頻子帶其系數(shù)較小、能量較低。當(dāng)噪聲比例較高時(shí),可將最高頻率子帶的系數(shù)全部看成噪聲,由此來估計(jì)噪聲的標(biāo)準(zhǔn)方差。基于主成分分析的噪聲估計(jì)方法利用主成分分析(PCA)提取數(shù)據(jù)中特征值最小的區(qū)域,估計(jì)噪聲水平[4-5]。這兩種方法得到的估計(jì)結(jié)果低于真實(shí)噪聲水平。而基于尺度不變性的估計(jì)方法利用峰度值不隨尺度改變、只由噪聲引起其系統(tǒng)變化的特性進(jìn)行噪聲估計(jì),效果要明顯優(yōu)于前兩種方法[6],尤其對(duì)邊緣信息豐富的低噪聲地震數(shù)據(jù)進(jìn)行估計(jì)時(shí)效果更好[7-8]。

傳統(tǒng)壓制隨機(jī)噪聲的方法可分為空間域和變換域兩種?？臻g域方法包括:均值濾波、中值濾波及各向異性擴(kuò)散濾波等[9-10]。變換域方法包括:傅里葉變換濾波、小波變換濾波、曲波變換閾值去噪等[11-12]。均值濾波是濾除高斯噪聲的常規(guī)算法,均值濾波算法將窗口內(nèi)所有數(shù)據(jù)點(diǎn)的均值作為中心數(shù)據(jù)的值,對(duì)所有數(shù)據(jù)都采用無差別濾波,這樣會(huì)使得非噪聲數(shù)據(jù)點(diǎn)的值也被改變,容易破壞邊緣細(xì)節(jié)信息,造成數(shù)據(jù)模糊。獨(dú)立分量分析技術(shù)(ICA)是實(shí)現(xiàn)盲源分離的有效手段,它考慮數(shù)據(jù)的高階統(tǒng)計(jì)特性,從而更加全面地表達(dá)數(shù)據(jù)的本質(zhì)特征[13]。K次迭代奇異值分解(K-SVD)數(shù)據(jù)降噪算法是一種基于K-means算法和奇異值分解(SVD)的方法[14],該方法降噪計(jì)算時(shí)間較長,且對(duì)高比例的噪聲數(shù)據(jù)降噪效果略顯不足[15]。為了解決這些問題,DONG等[16]提出了非局部集中稀疏表示(NCSR)數(shù)據(jù)修復(fù)算法,對(duì)信噪比較低的地震數(shù)據(jù)的降噪效果有了較大提升,但是它構(gòu)建用于稀疏表示所需字典時(shí)的計(jì)算量比較大,此外還可能過度地平滑數(shù)據(jù),造成數(shù)據(jù)失真。為了克服這種方法的缺點(diǎn),GU等[17]提出了加權(quán)核范數(shù)最小化(WNNM)降噪算法,該算法能根據(jù)矩陣奇異值刻畫數(shù)據(jù)差異,給定不同的權(quán)值,突顯數(shù)據(jù)中重要的信息,具有良好的去噪效果。在此基礎(chǔ)上發(fā)展起來的多道加權(quán)核范數(shù)最小化噪聲壓制算法考慮多道之間的噪聲差異,利用各道之間的冗余信息,同時(shí)考慮區(qū)分不同噪聲的統(tǒng)計(jì)性信息,去噪效果較為理想[18-19]。

本文研究了基于噪聲水平估計(jì)的加權(quán)核范數(shù)最小化噪聲壓制方法,在基于尺度不變性的噪聲方差估計(jì)方法的基礎(chǔ)上,采用加權(quán)核范數(shù)最小化方法去除地震數(shù)據(jù)中的隨機(jī)噪聲,采用理論模型數(shù)據(jù)和實(shí)際地震資料驗(yàn)證了方法的有效性。

1 方法原理

基于噪聲水平估計(jì)的加權(quán)核范數(shù)最小化噪聲壓制方法采用分塊處理的思路。首先利用塊匹配的思想對(duì)地震數(shù)據(jù)進(jìn)行分塊處理,然后利用基于尺度不變性的噪聲水平估計(jì)方法來計(jì)算分塊數(shù)據(jù)中的隨機(jī)噪聲方差,利用所得到的方差來歸一化加權(quán)核范數(shù)最小化算法中F范數(shù)的保真項(xiàng),保證加權(quán)核范數(shù)最小化算法能在壓制隨機(jī)噪聲的同時(shí)保護(hù)有效信號(hào),從而實(shí)現(xiàn)根據(jù)地震數(shù)據(jù)本身噪聲水平自適應(yīng)壓制隨機(jī)噪聲的目的。

1.1 基于尺度不變性的噪聲水平估計(jì)

為了實(shí)現(xiàn)自適應(yīng)的去噪處理,需要知道地震數(shù)據(jù)中隨機(jī)噪聲的分布特征,本文采用噪聲方差估計(jì)來衡量數(shù)據(jù)中的噪聲水平[20-21]。

用y表示某個(gè)含噪數(shù)據(jù),x表示與y對(duì)應(yīng)的不含噪聲的原始數(shù)據(jù),η表示噪聲,則:y=x+η。假定噪聲高斯分布,x可以用廣義高斯分布擬合,x的廣義高斯分布的峰度值κx(α)直接取決于其形狀參數(shù)α:

式中:Γ是標(biāo)準(zhǔn)伽馬函數(shù);κx(α)表示數(shù)據(jù)x的峰度值?？梢钥闯?峰度值κ與形狀參數(shù)α成反比。在實(shí)際地震剖面的噪聲估計(jì)中,峰度值κ為用方差平方歸一化處理的四階中心矩:

式中:μ4表示四階中心距;σ2為方差。由于噪聲的獨(dú)立性,含噪數(shù)據(jù)y的方差滿足:

y的四階中心矩滿足:

利用前面計(jì)算的方差的平方進(jìn)行歸一化處理:

根據(jù)上述結(jié)果,可以在給定無噪數(shù)據(jù)的情況下,從含噪數(shù)據(jù)中獲得邊緣濾波響應(yīng)分布的峰度值。但是,在實(shí)際處理中,無噪數(shù)據(jù)是未知的,因此假設(shè)原始無噪數(shù)據(jù)具有尺度不變特性,其邊緣濾波器響應(yīng)數(shù)據(jù)分布的峰度值是未知常數(shù),而在數(shù)據(jù)中加入噪聲會(huì)導(dǎo)致整個(gè)數(shù)據(jù)的峰度值發(fā)生變化。用以離散余弦變換(DCT)為基的濾波器對(duì)圖像作卷積處理,基函數(shù)為:

具體步驟如下:

1) 根據(jù)輸入的含噪數(shù)據(jù)y和規(guī)定的分塊大小N,由公式(6)得到N×N個(gè)以離散余弦變換為基的濾波器矩陣b;

2) 用每個(gè)濾波器矩陣b對(duì)圖像進(jìn)行卷積,生成響應(yīng)數(shù)據(jù)yi;

圖1給出了采用不同方法得到的隨機(jī)噪聲水平估計(jì)結(jié)果,噪聲為白噪。從圖1可以看出,基于尺度不變性的噪聲估計(jì)方法的估計(jì)水平要優(yōu)于小波分解估計(jì)和主成分分析方法,后二者估計(jì)結(jié)果與實(shí)際結(jié)果相比偏低?？梢悦黠@看出,在添加中、低比例噪聲時(shí),基于尺度不變性的噪聲估計(jì)方法具有較好且穩(wěn)定的估計(jì)結(jié)果,雖然添加高比例噪聲時(shí)其估計(jì)效果不如低比例噪聲時(shí)理想,但是估計(jì)結(jié)果仍然優(yōu)于其它兩種算法。因此基于尺度不變性的噪聲估計(jì)方法在估計(jì)復(fù)雜低比例噪聲時(shí),效果較好。

圖1 不同方法的噪聲水平估計(jì)結(jié)果

1.2 多道加權(quán)核范數(shù)最小化

對(duì)于疊后含噪數(shù)據(jù)y的局部分塊yj,可以采用塊匹配的方法在數(shù)據(jù)中找到它的非局部分塊[22]。將地震數(shù)據(jù)分成指定大小的數(shù)據(jù)塊,選取yj作為參考?jí)K,其它為候選塊,從候選塊中找出與參考?jí)K相似的塊,將這些非局部相似分塊組成矩陣Yj,其中兩個(gè)數(shù)據(jù)塊之間的相似性采用(8)式表示:

假設(shè)Xj和Nj分別是無噪數(shù)據(jù)分塊矩陣和噪聲數(shù)據(jù)分塊矩陣,則Yj=Xj+Nj,Xj為低秩矩陣,低秩矩陣逼近方式可以用于由Yj估計(jì)Xj。聚集所有經(jīng)過降噪的分塊,則可以估計(jì)整個(gè)數(shù)據(jù)[23]。

式中:‖·‖w,*表示加權(quán)矩陣奇異值的和,該方法的主要問題是確定權(quán)重向量w。權(quán)重向量與對(duì)應(yīng)的奇異值成反比,權(quán)重公式為:

式中:λi(Xj)是數(shù)據(jù)塊組Xj對(duì)應(yīng)的第i個(gè)奇異值;c為常數(shù);m是Yj中相似塊的數(shù)量。為了避免分母為0,令ε=10-16。λi(Xj)的奇異值不能直接獲得,可以近似為:

其中,λi(Yj)是Yj的第i個(gè)奇異值。

根據(jù)實(shí)際數(shù)據(jù)選擇迭代次數(shù)L和迭代權(quán)重參數(shù)δ,L的設(shè)置原則為處理結(jié)果趨于穩(wěn)定,根據(jù)試驗(yàn)結(jié)果分析,當(dāng)L=6,δ=0.1時(shí),去噪效果較好。

2 理論模型試算

2.1 模型數(shù)據(jù)去噪效果分析

建立如圖2a所示的合成記錄剖面,圖2b為含噪數(shù)據(jù),圖2c為采用本文方法得到的去噪結(jié)果,圖2d為濾除的噪聲剖面。圖3為圖2黑色框區(qū)域的放大圖像。由圖2和圖3可以看出,去除的噪聲中基本不包含有效的結(jié)構(gòu)信息,說明本文去噪方法能夠保持剖面的邊界特征;去噪剖面中基本沒有隨機(jī)噪聲的殘留,說明本文方法在未傷害有效信號(hào)的前提下,能夠有效地壓制加入的噪聲。

圖2 模型數(shù)據(jù)采用本文方法去噪前、后的結(jié)果a 合成數(shù)據(jù); b 含噪數(shù)據(jù); c 去噪結(jié)果; d 去除的噪聲

圖3 圖2的局部放大結(jié)果a 合成數(shù)據(jù); b 含噪數(shù)據(jù); c 去噪結(jié)果

對(duì)圖2b所示含噪數(shù)據(jù)分別采用SVD方法、正則化非平穩(wěn)自回歸(RNA)方法和本文方法進(jìn)行去噪處理,得到的結(jié)果如圖4所示。采用SVD方法得到的去噪結(jié)果(圖4a)噪聲殘留嚴(yán)重,去噪效果較差,去噪剖面上存在偽影現(xiàn)象;RNA方法的去噪結(jié)果較SVD方法要好(圖4b),但是去噪剖面上仍然殘留噪聲;與前兩種方法相比,本文方法去噪效果最好(圖4c),噪聲得到了壓制。對(duì)比圖4中箭頭所指區(qū)域,本文方法對(duì)剖面中不連續(xù)性信息保留效果較好,說明本文方法能夠在保護(hù)有效信號(hào)的同時(shí),壓制數(shù)據(jù)中的隨機(jī)噪聲。

2.2 斷面保持效果分析

為了驗(yàn)證本文方法的邊界保持效果,建立含斷層的模型,如圖5a所示;對(duì)圖5a所示模型加入噪聲得到的結(jié)果如圖5b所示;采用本文方法對(duì)圖5b 進(jìn)行去噪處理后的結(jié)果如圖5c所示;圖5d為去除的噪聲。由圖5d可以明顯看出,本文方法能夠有效壓制原始剖面中的隨機(jī)噪聲,而且能較好地保護(hù)邊緣處的不連續(xù)性信息,噪聲剖面中無有效信號(hào)。

圖4 對(duì)模型數(shù)據(jù)采用不同方法得到的去噪結(jié)果a SVD方法; b RNA方法; c 本文方法

圖5 含斷層模型采用本文方法去噪前、后的結(jié)果a 合成數(shù)據(jù); b 含噪數(shù)據(jù); c 去噪結(jié)果; d 去除的噪聲

3 實(shí)際數(shù)據(jù)處理

對(duì)某工區(qū)實(shí)際地震資料采用本文方法進(jìn)行去噪處理。圖6a為過該工區(qū)內(nèi)的一條疊加地震剖面,該剖面含有較高比例的隨機(jī)噪聲,同相軸模糊、邊界及斷點(diǎn)不夠清晰、地質(zhì)接觸關(guān)系模糊,給后續(xù)處理解釋工作帶來極大困難。采用本文方法處理后的地震剖面如圖6b所示。對(duì)比圖6a和圖6b可以看出,去噪后的地震剖面中同相軸連續(xù)性得到加強(qiáng),斷點(diǎn)清晰、斷層展布特征清楚。圖6c和圖6d分別為典型區(qū)域(圖6a和圖6b中黑框所示區(qū)域)放大顯示結(jié)果。從圖6c和圖6d可以看出,原始數(shù)據(jù)中的噪聲得到有效衰減,處理后的地震數(shù)據(jù)中反射波同相軸的連續(xù)性得到了增強(qiáng),且斷層特征保持良好,說明本文方法能為層位追蹤、斷層識(shí)別和構(gòu)造解釋等后續(xù)處理提供信噪比較高的地震數(shù)據(jù)。圖7a 和圖7b分別給出了對(duì)實(shí)際地震資料采用本文方法進(jìn)行去噪處理前、后的瞬時(shí)振幅屬性。對(duì)比圖7a 和圖7b可以看出,去噪前瞬時(shí)振幅屬性較為雜亂,難以識(shí)別斷點(diǎn)和巖性突變點(diǎn);采用本文方法去噪處理后改善了地震資料的品質(zhì),瞬時(shí)振幅屬性變得清晰和容易識(shí)別,不連續(xù)性信息得到很好的突出,增強(qiáng)了地質(zhì)體邊緣的清晰度(如圖中箭頭所指區(qū)域)。

圖6 實(shí)際地震資料采用本文方法去噪處理前、后的地震剖面a 去噪前的數(shù)據(jù); b 去噪后的數(shù)據(jù); c 圖6a局部放大結(jié)果; d 圖6b局部放大結(jié)果

圖7 對(duì)實(shí)際地震資料采用本文方法進(jìn)行去噪處理前(a)、后(b)的瞬時(shí)振幅屬性

4 結(jié)論

地震資料去噪處理需要既能夠衰減噪聲,提高信噪比,同時(shí)又能保留和增強(qiáng)地震反射信息中有效的不連續(xù)性信息。本文研究了基于噪聲水平估計(jì)的加權(quán)核范數(shù)最小化噪聲壓制方法,首先采用基于尺度不變性噪聲估計(jì)方法,對(duì)塊匹配后的地震數(shù)據(jù)進(jìn)行噪聲估計(jì),根據(jù)噪聲方差估計(jì)歸一化加權(quán)核范數(shù)最小化噪聲壓制算法中的保真項(xiàng),實(shí)現(xiàn)對(duì)地震數(shù)據(jù)的自適應(yīng)去噪處理。主要結(jié)論包括:

1) 用本文方法實(shí)現(xiàn)了對(duì)地震數(shù)據(jù)的去噪處理,并保護(hù)有效信號(hào),尤其在同相軸不連續(xù)區(qū)域,避免了采用簡單閾值處理方法產(chǎn)生的偽吉布斯震蕩現(xiàn)象;

2) 數(shù)值試驗(yàn)和實(shí)際數(shù)據(jù)的應(yīng)用結(jié)果表明,本文方法能夠在不破壞有效信息的同時(shí),自適應(yīng)地衰減地震數(shù)據(jù)中的隨機(jī)噪聲,處理后的地震數(shù)據(jù)信噪比得到提高,有利于后期的構(gòu)造解釋、斷層和斷點(diǎn)識(shí)別、層位追蹤、幾何屬性提取等。