王 偉,高靜懷,陳文超*,朱振宇
1 西安交通大學電子與信息工程學院波動與信息研究所,西安 710049
2 中海石油研究總院,北京 100027
基于結(jié)構(gòu)自適應(yīng)中值濾波器的隨機噪聲衰減方法
王 偉1,高靜懷1,陳文超1*,朱振宇2
1 西安交通大學電子與信息工程學院波動與信息研究所,西安 710049
2 中海石油研究總院,北京 100027
本文提出一種保護斷層、裂縫等地層邊緣特征的結(jié)構(gòu)自適應(yīng)中值濾波器,用于衰減地震資料中的隨機噪聲.基于地震反射同相軸局部呈線型結(jié)構(gòu)的假設(shè),采用梯度結(jié)構(gòu)張量估計地層傾向,分析地層結(jié)構(gòu)的規(guī)則程度,在此基礎(chǔ)上引入地震剖面中線型和橫向不連續(xù)性兩種結(jié)構(gòu)特征的置信度量.結(jié)構(gòu)自適應(yīng)中值濾波器根據(jù)這兩種置信度量調(diào)整濾波器窗函數(shù)的尺度和形狀,根據(jù)地層傾角調(diào)整濾波器窗函數(shù)的方向,從而使得濾波操作窗能夠最佳匹配信號的局部結(jié)構(gòu)特征.將本文方法用于合成和實際數(shù)據(jù)的處理,并與兩種常用中值濾波方法進行對比,結(jié)果表明,該方法能夠更好地解決地震剖面的隨機噪聲衰減和有效信號保真的問題,在增強反射同相軸的橫向一致性的同時有效保持了剖面內(nèi)的地層邊緣和細節(jié)特征,顯著改善了地震資料的品質(zhì).
中值濾波,自適應(yīng)濾波,噪聲衰減,構(gòu)造保護,梯度結(jié)構(gòu)張量
在石油天然氣勘探開發(fā)中,斷層、通道和裂縫等結(jié)構(gòu)特征以及河道砂體等沉積特征是發(fā)現(xiàn)和描述油藏的基礎(chǔ),分辨和分析地震資料中這些結(jié)構(gòu)特征具有重要的意義.然而,由于這些區(qū)域地質(zhì)形態(tài)的復雜性,使得地震數(shù)據(jù)中這些區(qū)域易受噪聲的干擾,無論采用人工解釋還是相干和邊緣檢測等自動檢測算法來描述和刻畫這些區(qū)域都變得較為困難.傳統(tǒng)的地震資料去噪方法,如F-K域傾向濾波、f-x域預測濾波和多道相干濾波等,多是利用地震反射同相軸的橫向一致性,易造成對傾斜同相軸振幅的壓制,也會模糊小的斷層和裂縫等不連續(xù)結(jié)構(gòu),甚至會引起大的斷層兩側(cè)同相軸的錯誤連接.如何能夠較好地平衡地震噪聲衰減和有效信號保護的關(guān)系受到了廣泛的關(guān)注[1-6],而能夠保護和增強斷裂等邊緣結(jié)構(gòu)的疊后噪聲衰減技術(shù)尤其受到地震解釋工作者的重視[7-9].
Luo等[10]基于次序統(tǒng)計的思想給出一種地震數(shù)據(jù)保邊濾波(EPS,Edge-Preserving Smoothing)算法,該方法采用Kuwahara多窗分析技術(shù),計算當前分析點周圍各個子窗內(nèi)數(shù)據(jù)的均值和方差,將最小方差對應(yīng)的子窗的均值作為當前點的濾波輸出.AlBinHassan等[11]發(fā)展了這種保邊濾波思想,給出了三維情況下的子窗剖分方法,將EPS算法用于三維地震數(shù)據(jù)的保邊濾波處理.Lu等[12]基于多項式擬合的思想推廣了一維的EPS算法,該方法通過計算包含當前分析點的一系列平移子窗內(nèi)信號對給定階次的多項式擬合誤差,選擇出最小擬合誤差對應(yīng)的子窗的信號估計值作為當前點的濾波輸出.楊培杰等[13]也基于這種多窗分析思想給出一種方向性邊界保持斷層增強技術(shù),使得相干體中的斷層圖像得到更加清晰顯示.Liu等[14-15]應(yīng)用圖像處理中的二維多級中值濾波器消除地震隨機噪聲,并分析了濾波器尺度參數(shù)選擇對噪聲衰減和信號細節(jié)保護的影響.此外,F(xiàn)ehmers等[16]引入微分方程濾波法用于地震資料的保邊濾波處理,該方法采用梯度結(jié)構(gòu)張量估計地層方向,度量地層的連續(xù)性,以此約束擴散濾波方程的濾波方向和濾波程度.Lavialle等[17]結(jié)合線型和面型擴散濾波模型給出專門針對保護和增強斷層構(gòu)造的斷層保持擴散濾波器.孫夕平等[18]和王旭松等[19]通過在擴散方程中引入二階導數(shù)項、張爾華等[20]通過修正擴散方程中的保邊約束項進一步改善非線性各向異性擴散濾波器的邊緣保護性能.然而,基于各向異性擴散方程的濾波方法以原始地震數(shù)據(jù)作為初始條件,通過擴散時間的迭代推移形成濾波作用,由于不好估計出獲得最佳濾波效果的擴散時間,導致不能準確預測其濾波性態(tài)[21].
經(jīng)典的二維中值濾波法具有優(yōu)異的噪聲衰減性能,尤其適用于消除非平穩(wěn)信號中的峰值噪聲.但是,由于在濾波過程中其濾波窗口保持恒定,且濾波操作沒有方向特性,使得采用中值濾波法處理具有層狀結(jié)構(gòu)的地震信號時,必然會引起有用信號的損傷和不連續(xù)結(jié)構(gòu)特征的模糊.一些改進的中值濾波方法如二維多級中值濾波器[14-15],雖然有較強的信號細節(jié)保護能力,但是其對噪聲的衰減能力不夠.本文結(jié)合微分方程濾波法中利用梯度結(jié)構(gòu)張量刻畫地震反射層局部結(jié)構(gòu)特征的方法,定義了帶有能量變化強度信息的橫向不連續(xù)性置信度量,結(jié)合線型置信度量和橫向不連續(xù)性置信度量,給出一種二維結(jié)構(gòu)自適應(yīng)中值濾波方法,使中值濾波操作窗函數(shù)根據(jù)信號局部結(jié)構(gòu)特征自適應(yīng)地作方向延伸和尺度伸縮,從而既能最大限度地消除隨機噪聲,又能有效保護有用信號和地層邊緣等細節(jié)特征.
有多種方法可用于地層結(jié)構(gòu)的規(guī)則性分析[16,22-25],梯度結(jié)構(gòu)張量法是其中一種簡單而又有效的途徑.記二維地震剖面為u(x,t),則該二維剖面對應(yīng)的梯度結(jié)構(gòu)張量由下式描述
其中,和T分別為梯度算子和轉(zhuǎn)置算子,?表示卷積運算,Gρ表示尺度為ρ的二維高斯函數(shù)Gρ(x,t)=exp(-(x2+t2)/2ρ2).在梯度結(jié)構(gòu)張量的定義中,由梯度向量Δu及其轉(zhuǎn)置向量形成的張量矩陣與尺度為ρ的高斯低通濾波器作用,確定了梯度結(jié)構(gòu)張量可度量的信號特征的最大尺度.將該對稱半正定矩陣Sρ(Δu)作矩陣特征分解
其中,特征值μ1≥μ2≥0反應(yīng)了信號在位置(x,t)的高斯鄰域Gρ內(nèi)沿特征方向的振幅變化強度;特征向量v1和v2給出局部正交方位,v1對應(yīng)于局部振幅變化最大的方向,即信號梯度方向,而v2對應(yīng)于局部振幅變化最小的方向,即反射同相軸的取向.根據(jù)梯度結(jié)構(gòu)張量兩個特征值的取值情況,可分辨不同的信號結(jié)構(gòu)特征,文獻[26]中引入線型信號結(jié)構(gòu)的置信度量:
該置信度量在區(qū)間[0,1]之間取值,衡量局部信號特征與線型結(jié)構(gòu)的相似程度,同時亦表征局部方位估計的準確程度.若CL→1,則表示局部信號特征趨于線型結(jié)構(gòu),估計的地層傾角的準確度越高;反之,若CL→0,則局部信號結(jié)構(gòu)背離了線型假設(shè),這可能由干擾噪聲和斷層、裂縫等橫向不連續(xù)性造成.根據(jù)線型結(jié)構(gòu)置信度量可較好地識別局部呈線型結(jié)構(gòu)的反射同相軸;但是,由于該置信度量值判斷的是信號能量的相對變化程度,容易受到噪聲的干擾,并且在無明顯信號結(jié)構(gòu)的區(qū)域內(nèi)取值沒有意義,從而很難從干擾背景中明確地識別斷層、裂縫等橫向不連續(xù)結(jié)構(gòu).
基于上述分析,本文在線型結(jié)構(gòu)置信度量的基礎(chǔ)上,給出一種含有信號橫向能量變化強度信息的置信度量值,用于識別斷層、裂縫等橫向不連續(xù)結(jié)構(gòu).
(4)式中(1-CL)項表示信號局部結(jié)構(gòu)特征相對線型結(jié)構(gòu)的背離程度,μ2表征在給定度量鄰域內(nèi)在均方誤差最小意義下信號沿局部一致方向的能量變化強度[24].在反射同相軸呈線型結(jié)構(gòu)的區(qū)域,由于CL→1,即(1-CL)→0,此時沿地層傾向的信號變化相對較小,即μ2→0,從而使得CI→0.在含有斷層、裂縫等橫向不連續(xù)性的區(qū)域,由于CL→0,即(1-CL)→1;且在橫向不連續(xù)區(qū)域,由于信號結(jié)構(gòu)背離了局部線型結(jié)構(gòu)假設(shè),此處信號的一致性取向不明確,因而沿局部一致性方向的信號變化強度也相對較大,即μ2?0,從而使得CI?0.因此,根據(jù)橫向不連續(xù)性置信度量CI可以很好地識別斷層和裂縫等橫向不連續(xù)結(jié)構(gòu).由于信號橫向能量變化項μ2的存在,對于噪聲能量弱于信號能量的情況,CI具有良好的抗噪能力;此外,CI也能夠判別無明顯信號結(jié)構(gòu)的區(qū)域(此時μ2→0,即CI→0),這對采用CI指導濾波過程尤為重要,可以避免一些濾波假象的產(chǎn)生.
對于二維地震剖面u(x),記x=(x,t)為采樣點的空間、時間位置,定義位置x處的結(jié)構(gòu)自適應(yīng)中值濾波操作為:
式中,M(x,y)表示位置x處的結(jié)構(gòu)自適應(yīng)濾波窗口(見圖1a),y為當前濾波窗口包含的采樣點的空間、時間位置,(x)表示位置x處的結(jié)構(gòu)自適應(yīng)中值濾波的輸出結(jié)果.結(jié)構(gòu)自適應(yīng)濾波窗M(x,y)是濾波位置x的函數(shù),其大小和形狀取決于x處分析鄰域內(nèi)的信號結(jié)構(gòu)特征:
其中,·表示內(nèi)積運算;n(x)和n⊥(x)分別對應(yīng)于x處的地層傾向及梯度方向,取為梯度結(jié)構(gòu)張量的兩個特征向量,即n(x)=v2(x),n⊥(x)=v1(x);σ1(x)和σ2(x)分別表示橢圓形濾波窗口的長軸和短軸,它們共同決定了結(jié)構(gòu)自適應(yīng)中值濾波器的濾波尺度以及濾波操作的方向選擇性.濾波器尺度參數(shù)的設(shè)計是影響結(jié)構(gòu)自適應(yīng)中值濾波器性態(tài)的關(guān)鍵,根據(jù)結(jié)構(gòu)自適應(yīng)濾波操作的要求,本文結(jié)合線型結(jié)構(gòu)置信度量CL和橫向不連續(xù)性置信度量CI,給出如下的濾波器尺度參數(shù)自適應(yīng)選擇策略:
式中Rmax規(guī)定了橢圓形濾波窗口的最大尺寸,g(·)為關(guān)于CI(x)的單調(diào)減函數(shù),其取值范圍限定為(0,1],文中取g(·)為指數(shù)函數(shù)
其中,β為針對CI(x)的閾值參數(shù),控制指數(shù)函數(shù)的衰減速度,β取值越小,指數(shù)函數(shù)的下降越快;反之,指數(shù)函數(shù)的衰減變緩.
結(jié)構(gòu)自適應(yīng)中值濾波器根據(jù)信號的局部結(jié)構(gòu)特征具有靈活的濾波尺度調(diào)節(jié)功能和方向選擇特性,控制最大濾波尺度的σ1(x)由橫向不連續(xù)性置信度量CI(x)決定,而調(diào)整方向選擇特性的σ2(x)由線型結(jié)構(gòu)置信度量CL(x)決定.在斷層、裂縫等橫向不連續(xù)區(qū)域由于CI(x)?0,使得橢圓型窗函數(shù)的長軸σ1(x)趨短;而在其他區(qū)域隨CI(x)→0,有σ1(x)→Rmax.在反射同相軸呈線型結(jié)構(gòu)區(qū)域,由于CL(x)→1,使得橢圓型窗函數(shù)的短軸σ2(x)萎縮;而在其他區(qū)域隨CL(x)→0,有σ2(x)→σ1(x).因此,如圖1b所示,在反射同相軸局部呈線型延伸的區(qū)域,濾波窗口沿同相軸取向拉伸為狹長的各向異性窗口(位置x1區(qū)域);在斷層等橫向不連續(xù)區(qū)域,濾波窗口自適應(yīng)地萎縮為小窗以匹配此處的不連續(xù)結(jié)構(gòu)特征(位置x2區(qū)域);而在無明顯信號結(jié)構(gòu)的噪聲干擾區(qū)域,濾波窗口延展為大的近似各向同性窗口(位置x3區(qū)域).
圖1 結(jié)構(gòu)自適應(yīng)中值濾波器窗口的構(gòu)造(a)濾波器窗口函數(shù);(b)信號結(jié)構(gòu)自適應(yīng)濾波策略.Fig.1 Filter window construction of structure-adaptive median filter(a)Filter window function;(b)Signal structure-consistent filtering mechanism.
聯(lián)合式(7)和(8)可見,結(jié)構(gòu)自適應(yīng)中值濾波器的濾波性能受到參數(shù)Rmax和β共同約束,最大尺寸參數(shù)Rmax控制濾波器的平滑性能,閾值參數(shù)β則調(diào)整濾波器對斷層等邊緣結(jié)構(gòu)的保持程度.考慮到文中采用梯度結(jié)構(gòu)張量估計地層的傾角和規(guī)則程度,梯度結(jié)構(gòu)張量的尺度參數(shù)ρ決定了其對局部信號結(jié)構(gòu)的度量孔徑,因此為了確?;趦煞N置信度量約束的濾波過程可信,要求濾波窗口的最大尺寸Rmax不超過2ρ,即取Rmax≤2ρ.對于閾值參數(shù)β,由圖2可見,當β取值偏小時,稍大的橫向不連續(xù)性置信度量CI(x)即引起濾波窗口的尺寸減小,因而會使一些弱的斷層、裂縫等邊緣結(jié)構(gòu)得到有效保護,但同時也會減弱濾波器在橫向不連續(xù)性區(qū)域的濾波效果;當β取值偏大時,只有很大的CI(x)才會導致濾波窗口萎縮到足夠小,因而在這種情況下濾波器的噪聲衰減能力得到增強,但在一定程度上會造成弱邊緣結(jié)構(gòu)的模糊和壓制.
圖2 結(jié)構(gòu)自適應(yīng)中值濾波窗的長軸控制函數(shù)Fig.2 Controlling functions of the main axis of the structure-adaptive median filter kernel
結(jié)構(gòu)自適應(yīng)中值濾波器的邊緣結(jié)構(gòu)保持性能受閾值參數(shù)β控制,選擇合適的β是保證濾波過程中能夠有效保持斷層、裂縫等地層邊緣的關(guān)鍵.然而,地震資料通常從淺層到深層,即使在同一層段內(nèi)信號橫向能量變化的動態(tài)范圍波動也很大,很難統(tǒng)一地選擇一個全局閾值參數(shù)β.本文給出一種區(qū)域分割的閾值參數(shù)自適應(yīng)選擇策略,即根據(jù)待處理剖面內(nèi)橫向能量變化在時間和空間上的分布情況,將剖面Ω分割為若干具有一定邊界重疊的子區(qū)域Ω=∪Ωi,在每個子區(qū)域Ωi內(nèi)取橫向不連續(xù)性置信度量CI(x)的極大值,即
從而通過下式得到該區(qū)域內(nèi)的自適應(yīng)閾值參數(shù)β為
式中α為根據(jù)保邊緣濾波要求而全局確定的百分比因子,thr是由剖面的整體噪聲干擾水平所確定的基底噪聲閾值.
通過將待處理剖面分割為若干子區(qū)域,由(9)和(10)式將(8)式中對閾值參數(shù)β的絕對取值轉(zhuǎn)化為對百分比因子α的相對取值,從而能夠統(tǒng)一調(diào)整結(jié)構(gòu)自適應(yīng)中值濾波器的邊緣保持性態(tài).在實際資料處理中,并不需要精確地確定子區(qū)域的大小,一種典型的選擇是每個子區(qū)域內(nèi)包含100道,每道包含150個采樣點,即可滿足濾波過程保持斷層、裂縫等邊緣特征的要求.
為了驗證結(jié)構(gòu)自適應(yīng)中值濾波器的保護邊緣去噪的有效性,本文選取主頻為20Hz的Ricker子波,采用射線追蹤方法得到一個含有斷層的復雜地層結(jié)構(gòu)模型的合成地震數(shù)據(jù)(見圖3a),在合成數(shù)據(jù)中加入5dB的帶限高斯隨機噪聲,得到圖3b的含噪數(shù)據(jù).
圖3 含有斷層的復雜地層結(jié)構(gòu)的合成地震數(shù)據(jù)(a)原始數(shù)據(jù);(b)加噪數(shù)據(jù).Fig.3 Synthetic seismic data of complex geology model with a fault(a)Clean data;(b)Noisy data.
采用結(jié)構(gòu)自適應(yīng)中值濾波器處理加噪的合成數(shù)據(jù),考慮到合成數(shù)據(jù)的信噪比很低,在梯度結(jié)構(gòu)張量的計算中,二維高斯低通濾波函數(shù)取較大的尺度參數(shù)ρ=4.0,結(jié)構(gòu)自適應(yīng)中值濾波器的最大濾波尺寸取為Rmax=ρ=4.0,即對應(yīng)最大尺度為9點的濾波窗口.針對圖3b中黑色長方形框所選定的典型區(qū)域,分析閾值百分比因子α對濾波結(jié)果的影響.調(diào)整α從10%開始以10%的間隔變化到200%,記錄每次濾波結(jié)果的信噪比,得到該區(qū)域的濾波結(jié)果的信噪比隨閾值百分比因子α的變化曲線,如圖4所示.從圖中可以看出,信噪比曲線隨α的增大而升高,且曲線在α=90%后升高逐漸變緩.這是由于隨著α的增加,濾波器的保邊緣性能減弱,當α>90%后,濾波過程已造成斷面處波形的損傷,在一定程度上抵消了濾波性能增強帶來的信噪比增加.因而在處理該加噪的合成數(shù)據(jù)時,選擇閾值百分比因子α=90%,此時對應(yīng)的結(jié)構(gòu)自適應(yīng)中值濾波器的各向異性窗函數(shù)的長軸σ1和短軸σ2的取值分布如圖5所示.對比圖5a和圖5b可見,結(jié)構(gòu)自適應(yīng)中值濾波器的尺度參數(shù)能夠很好地匹配合成數(shù)據(jù)的信號結(jié)構(gòu)特征,在局部呈線型結(jié)構(gòu)的同相軸區(qū)域橢圓形濾波窗口的長軸拉伸、短軸萎縮,形成有強方向選擇特性的各向異性窗口;在斷層區(qū)域長軸和短軸同時萎縮,形成小尺度分析窗口;而在無信號特征的區(qū)域,長軸和短軸均被拉伸,形成大尺度的各向同性窗口.
圖4 不同α對應(yīng)的結(jié)構(gòu)自適應(yīng)中值濾波結(jié)果的信噪比Fig.4 SNR of filtered results by structure-adaptive median filter using differentα
圖5 針對加噪合成數(shù)據(jù)的結(jié)構(gòu)自適應(yīng)中值濾波器的窗口參數(shù)(a)橢圓形濾波窗口的長軸σ1;(b)橢圓形濾波窗口的短軸σ2.Fig.5 Window parameters of structure-adaptive median filter corresponding to the noisy synthetic data(a)Main axesσ1of the elliptic filter window;(b)The second axesσ2of the elliptic filter window.
加噪合成數(shù)據(jù)的結(jié)構(gòu)自適應(yīng)中值濾波結(jié)果及其得到的噪聲剖面分別如圖6a和圖6d所示.為了對比本文方法的處理效果,使用經(jīng)典的二維中值濾波器和具有信號細節(jié)保護能力的二維多級中值濾波器處理加噪合成數(shù)據(jù).取兩種中值濾波器的窗口尺寸與結(jié)構(gòu)自適應(yīng)中值濾波器的最大窗口尺寸相同,即二維中值濾波器選擇9×9的濾波窗口,二維多級中值濾波器取9點的濾波長度,兩者對應(yīng)的濾波結(jié)果以及得到的噪聲剖面分別如圖6b-c和圖6e-f所示.
從圖6a可見,經(jīng)過結(jié)構(gòu)自適應(yīng)中值濾波后,隨機噪聲得到很大程度地衰減,反射同相軸的一致性得到增強,剖面整體變得非常干凈;從圖6d可見,結(jié)構(gòu)自適應(yīng)中值濾波器得到的噪聲剖面中只有隨機分布的干擾噪聲,無明顯的有效信號結(jié)構(gòu),斷層結(jié)構(gòu)在濾波過程中得到了有效保持.對比圖6b-c中兩種中值濾波方法的濾波結(jié)果,9×9的二維中值濾波器消除了大部分的隨機噪聲,但也造成有效信號的損傷,從其得到的噪聲剖面圖6e可見,濾波過程引起斷層信息的模糊,也造成反射同相軸的損傷,尤其對傾斜反射結(jié)構(gòu)的損傷更為明顯;9點濾波長度的二維多級中值濾波器雖然有較強的信號細節(jié)保護能力,但是從其得到的噪聲剖面圖6f可見,二維多級中值濾波器對隨機噪聲衰減程度有限.
從上述分析、比較可以看出:在取相同濾波窗口尺度的條件下,相比二維中值濾波器和二維多級中值濾波器,結(jié)構(gòu)自適應(yīng)中值濾波器能夠根據(jù)地震剖面的信號結(jié)構(gòu)特征,自適應(yīng)地構(gòu)造與信號結(jié)構(gòu)相吻合的濾波窗口,在消除隨機噪聲和保護斷層等邊緣結(jié)構(gòu)之間找到了很好折中;結(jié)構(gòu)自適應(yīng)中值濾波器的噪聲衰減性能接近9×9的二維中值濾波器,而其對斷層及有效信號保護能力又優(yōu)于9點濾波長度的二維多級中值濾波器.
圖6 加噪合成數(shù)據(jù)的濾波結(jié)果比較(a)結(jié)構(gòu)自適應(yīng)中值濾波結(jié)果;(b)二維中值濾波結(jié)果;(c)二維多級中值濾波結(jié)果;(d)結(jié)構(gòu)自適應(yīng)中值濾波器得到的噪聲剖面;(e)二維中值濾波器得到的噪聲剖面;(f)二維多級中值濾波器得到的噪聲剖面.Fig.6 Comparison of filtering results of the noisy synthetic data(a)Result obtained by structure-adaptive median filter;(b)Result obtained by 2Dmedian filter;(c)Result obtained by 2Dmultistage median filter;(d)Noise removed by structure-adaptive median filter;(e)Noise removed by 2Dmedian filter;(f)Noise removed by 2D multistage median filter.
圖7 XX油田實際地震剖面Fig.7 Real seismic profile from XX oil field
基于理論分析結(jié)果,進一步將結(jié)構(gòu)自適應(yīng)中值濾波器用于疊后地震剖面的濾波處理.圖7中給出某油田的疊后純波地震資料,剖面中含有豐富的不連續(xù)結(jié)構(gòu)特征,除了兩條主要的斷層外,在中央部位有一處雜亂反射區(qū)域,在右側(cè)偏下部有一處振幅異常區(qū)域.由于受到隨機噪聲和采集腳印的干擾,同相軸的連續(xù)性受到一定影響,不利于層位的橫向追蹤、不連續(xù)結(jié)構(gòu)的自動檢測和拾取等后繼的解釋和處理.運用本文的結(jié)構(gòu)自適應(yīng)中值濾波器處理該剖面,選擇100道、每道150個采樣點的子區(qū)域分割;在梯度結(jié)構(gòu)張量的計算中取高斯濾波器的尺度參數(shù)ρ=3.0;結(jié)構(gòu)自適應(yīng)中值濾波器的最大濾波尺寸取為Rmax=4.0,即對應(yīng)最大為9點的濾波窗口;為了達到較好的邊緣保持性能,取閾值百分比因子α=50%,對應(yīng)的濾波結(jié)果和得到的噪聲剖面分別如圖8a和圖8d所示.作為對比,同樣采用兩種中值濾波方法處理該剖面.由于二維中值濾波器的濾波窗口選擇得過大,對有效信號損傷嚴重,而選擇得過小,會導致噪聲衰減不夠,故此處取二維中值濾波器的窗口尺寸為5×5,而二維多級中值濾波器采用9點濾波長度,即與結(jié)構(gòu)自適應(yīng)中值濾波器的最大濾波尺寸相當,對應(yīng)的濾波結(jié)果和它們得到的噪聲剖面分別如圖8b-c、e-f所示.
圖8 實際地震資料濾波結(jié)果比較(a)結(jié)構(gòu)自適應(yīng)中值濾波結(jié)果;(b)二維中值濾波結(jié)果;(c)二維多級中值濾波結(jié)果;(d)結(jié)構(gòu)自適應(yīng)中值濾波器得到的噪聲剖面;(e)二維中值濾波器得到的噪聲剖面;(f)二維多級中值濾波器得到的噪聲剖面.Fig.8 Comparison of filtering results of real seismic profile(a)Result obtained by structure-adaptive median filter;(b)Result obtained by 2Dmedian filter;(c)Result obtained by 2Dmultistage median filter;(d)Noise removed by structure-adaptive median filter;(e)Noise removed by 2Dmedian filter;(f)Noise removed by 2Dmultistage median filter.
對比結(jié)構(gòu)自適應(yīng)中值濾波前后的剖面圖7和圖8a可見,結(jié)構(gòu)自適應(yīng)中值濾波過程消除了大部分的干擾噪聲,改善了反射同相軸的一致性,斷層構(gòu)造得到增強,振幅異常區(qū)域變得清晰.分析結(jié)構(gòu)自適應(yīng)中值濾波得到的噪聲剖面(圖8d)可見,淺層的干擾噪聲非常強,中深層的干擾噪聲相對較弱,剖面中除隨機噪聲外還有一些大角度的相干噪聲;對比原始地震剖面圖7,噪聲剖面中看不到明顯的有效信號結(jié)構(gòu).分析圖8b中二維中值濾波器的濾波結(jié)果,從視覺上看噪聲基本被壓制,信噪比得到很大的提高,但是信號的邊緣結(jié)構(gòu)尤其斷層、雜亂反射等被嚴重模糊,這一點可以從其得到的噪聲剖面(圖8e)得到印證.從圖8c和圖8f可見,二維多級中值濾波器僅僅沿同相軸走向有微弱的噪聲衰減作用,濾波能力非常有限.以上結(jié)果表明,結(jié)構(gòu)自適應(yīng)中值濾波器對淺層的強噪聲和深層的弱噪聲均能夠有效壓制,不僅能夠衰減隨機噪聲,也能夠消除部分相干噪聲,并且很好地保持了有效信號和斷層、振幅異常等地層邊緣和細節(jié)特征;結(jié)構(gòu)自適應(yīng)中值濾波器不僅具有二維中值濾波器優(yōu)異的噪聲衰減能力,而且又不會因保護地層邊緣和細節(jié)特征的需求而像二維多級中值濾波器那樣引起濾波能力的明顯降低.
對原始地震剖面和三種濾波方法的處理結(jié)果進行頻譜分析,圖9給出相應(yīng)剖面的平均振幅譜.對比處理前后的頻譜發(fā)現(xiàn),三種方法處理后的高頻段隨機噪聲幅度均要低于處理前,且頻譜的主要結(jié)構(gòu)特點沒有被破壞;然而結(jié)構(gòu)自適應(yīng)中值濾波器和二維多級中值濾波器均較好地保持了主頻和有用信息頻段的振幅譜,而二維中值濾波器對主頻和有用信息頻段的振幅譜的損傷較為嚴重.以上結(jié)果說明結(jié)構(gòu)自適應(yīng)中值濾波器能夠有效壓制地震隨機噪聲,并且?guī)缀醪粨p傷有用信息頻段的頻譜.
圖9 實際地震資料濾波結(jié)果平均振幅譜對比(a)圖7的平均振幅譜;(b)圖8(a)的平均振幅譜;(c)圖8(b)的平均振幅譜;(d)圖8(c)的平均振幅譜.Fig.9 Comparison of average spectrum of filtering results of real seismic profile(a)Average spectrum of Fig.7;(b)Average spectrum of Fig.8(a);(c)Average spectrum of Fig.8(b);(d)Average spectrum of Fig.8(c).
在實際地震資料的處理中,中值濾波器的應(yīng)用已經(jīng)比較成熟.本文針對現(xiàn)有的一些中值濾波方法在疊后地震資料隨機噪聲衰減方面存在的缺陷,提出了一種結(jié)構(gòu)自適應(yīng)中值濾波器.該方法依據(jù)梯度結(jié)構(gòu)張量對分析鄰域內(nèi)地層結(jié)構(gòu)規(guī)則程度的估計情況,引入了地震信號的兩種結(jié)構(gòu)特征的置信度量,利用它們調(diào)控結(jié)構(gòu)自適應(yīng)中值濾波器的濾波窗口的尺度和形狀,使得濾波窗口能夠最佳地匹配處理鄰域內(nèi)的地層結(jié)構(gòu)特征.通過合成模型和實際地震資料的處理結(jié)果表明,相比經(jīng)典的二維中值濾波器及二維多級中值濾波器,本文方法不僅能夠有效地衰減隨機噪聲、提高剖面的信噪比,而且能夠較好地保持有效信號和斷層、裂縫等地層邊緣和細節(jié)特征.本文方法可作為一個預處理過程,濾波結(jié)果進一步用于層位的自動追蹤、對噪聲敏感屬性的提取等后繼的解釋和處理流程.
(References)
[1] 李月,楊寶俊,趙雪平等.檢測地震勘探微弱同相軸的混沌振子算法.地球物理學報,2005,48(6):1428-1433.Li Y,Yang B J,Zhao X P,et al.An algorithm of chaotic vibrator to detect weak events in seismic prospecting records.Chinese J.Geophys.(in Chinese),2005,48(6):1428-1433.
[2] 高靜懷,毛劍,滿蔚仕等.疊前地震資料噪聲衰減的小波域方法研究.地球物理學報,2006,49(4):1155-1163.Gao J H,Mao J,Man W S,et al.On the denoising method of prestack seismic data in wavelet domain.Chinese J.Geophys.(in Chinese),2006,49(4):1155-1163.
[3] Neelamani R,Baumstein A I,Gillard D G,et al.Coherent and random noise attenuation using the curvelet transform.The Leading Edge,2008,27(2):240-248.
[4] 劉洋,F(xiàn)omel S,劉財?shù)?高階seislet變換及其在隨機噪聲消除中的應(yīng)用.地球物理學報,2009,52(8):2142-2151.Liu Y,F(xiàn)omel S,Liu C,et al.High-order seislet transform and its application of random noise attenuation.Chinese J.Geophys.(in Chinese),2009,52(8):2142-2151.
[5] 陸文凱.基于離散余弦變換的地震隨機噪聲壓制技術(shù).石油地球物理勘探,2011,46(2):202-206.Lu W K.Seismic random noise suppression based on the discrete cosine transform.Oil Geophysical Prospecting(in Chinese),2011,46(2):202-206.
[6] 林紅波,李月,徐學純.壓制地震勘探隨機噪聲的分段時頻峰值濾波方法.地球物理學報,2011,54(5):1358-1366.Lin H B,Li Y,Xu X C.Segmenting time-frequency peak filtering method to attenuation of seismic random noise.Chinese J.Geophys.(in Chinese),2011,54(5):1358-1366.
[7] Chopra S,Marfurt K J.Seismic attributes–a historical perspective.Geophysics,2005,70(5):3SO-28SO.
[8] Chopra S,Marfurt K J.Emerging and future trends in seismic attributes.The Leading Edge,2008,27(3):298-318.
[9] Chopra S,Misra S,Marfurt K J.Coherence and curvature attributes on preconditioned seismic data.The Leading Edge,2011,30(4)386-393.
[10] Luo Y,Marhoon M,Al-Dossary S,et al.Acquistion Processing—Edge-preserving smoothing and applications.The Leading Edge,2002,21(2):136-158.
[11] AlBinHassan N M,Luo Y,Al-Faraj M N.3Dedgepreserving smoothing and applications.Geophysics,2006,71(4):P5-P11.
[12] Lu Y H,Lu W K.Edge-preserving polynomial fitting method to suppress random seismic noise.Geophysics,2009,74(4):V69-V73.
[13] 楊培杰,穆星,張景濤.方向性邊界保持斷層增強技術(shù).地球物理學報,2010,53(12):2992-2997.Yang P J,Mu X,Zhang J T.Orientational edge preserving fault enhance.Chinese J.Gephys.(in Chinese),2010,53(12):2992-2997.
[14] 劉財,王典,劉洋等.二維多級中值濾波技術(shù)在隨機噪聲消除中的應(yīng)用初探.石油地球物理勘探,2005,40(2):163-167.Liu C,Wang D,Liu Y,et al.Preliminary study of using 2D multi-level median filtering technique to eliminate random noises.Oil Geophysical Prospecting(in Chinese),2005,40(2):163-167.
[15] Liu C,Liu Y,Yang B J,et al.A 2Dmultistage median filter to reduce random seismic noise.Geophysics,2006,71(5):V105-V110.
[16] Fehmers G C,H?cker C F W.Fast structural interpretation with structure-oriented filtering.Geophysics,2003,68(4):1286-1293.
[17] Lavialle O,Pop S,Germain C,et al.Seismic fault preserving diffusion.Journal of Applied Geophysics,2007,61(2):132-141.
[18] 孫夕平,杜世通,湯磊.相干增強各向異性擴散濾波技術(shù).石油地球物理勘探,2004,39(6):651-655.Sun X P,Du S T,Tang L.Coherent-enhancing anisotropic diffusion filtering technique.Oil Geophysical Prospecting(in Chinese),2004,39(6):651-655.
[19] 王旭松,楊長春.對地震圖像進行保邊濾波的非線性各向異性擴散算法.地球物理學進展,2006,21(2):452-457.Wang X S,Yang C C.An edge-preserving smoothing algorithm of seismic image using nonlinear anisotropic diffusion equation.Progress in Geophys.(in Chinese),2006,21(2):452-457.
[20] 張爾華,王偉,高靜懷等.非線性各向異性擴散濾波器用于三維地震資料噪聲衰減與結(jié)構(gòu)特征增強.地球物理學進展,2010,25(3):866-870.Zhang E H,Wang W,Gao J H,et al.Non-linear anisotropic diffusion filtering for 3Dseismic noise removal and structure enhancement.Progress in Geophys.(in Chinese),2010,25(3):866-870.
[21] Gómez G.Local smoothness in terms of variance.Proceedings of the BMVC,2000:815-824.
[22] Marfurt K J,Sudhaker V,Gersztenkorn A,et al.Coherency calculations in the presence of structural dip.Geophysics,1999,64(1):104-111.
[23] Randen T,Monsen E,Signer C,et al.Three-dimensional texture attributes for seismic data analysis.70th Annual International Meeting,SEG Expanded Abstracts,2000:668-671.
[24] Wang Y E,Luo Y,Al-faraj M N.Inverse-vector approach for smoothing dips and azimuths.Geophysics,2008,73(6):P9-P14.
[25] Gao D L.Latest developments in seismic texture analysis for subsurface structure,facies,and reservoir characterization:A review.Geophysics,2011,76(2):W1-W13.
[26] Bakker P.Image structure analysis for seismic interpretation[Ph.D′s thesis].Netherlands:Technische Universiteit Delft,2002.
Random seismic noise suppression via structure-adaptive median filter
WANG Wei1,GAO Jing-Huai1,CHEN Wen-Chao1*,ZHU Zhen-Yu2
1 School of Electronic and Information Engineering,Xi′an Jiaotong University,Xi′an 710049,China
2 CNOOC Research Institute,Beijing100027,China
This paper presents a structure-adaptive median filter for reducing random seismic noise which preserves seismic edges and details such as faults and fractures.By considering seismic horizons as line-like structures,this filtering framework relies on the computation of the Gradient Structure Tensors,which provides dips of seismic events and the regularity of local seismic structures,based on which two confidence measures are defined for different seismic structures.The structure-adaptive median filter adjusts the shape of the filter kernel according to the two confidence measures and aligns the filtering kernel along local dips to make the filtering kernel optimally matched with different local geological features.The proposed filter has been applied to both the synthetic and real data,and compared with two widely used median filtering schemes.The processing and comparison results show that the proposed structure-adaptive median filter is more suitable to balance suppressing random seismic noise and preserving signals which preserves seismic edges and details while enhancing the coherence of seismic horizons,and improves the quality of seismic images significantly.
Median filter,Adaptive filtering,Noise suppression,Structure preserving,Gradient structure tensor
10.6038/j.issn.0001-5733.2012.05.030
P631
2011-04-06,2011-06-28收修定稿
國家自然科學基金重點項目(40730424)、國家自然科學基金面上項目(40674064)、國家油氣重大專項(2011ZX05023-005-009)資助.
王偉,男,1983年生,現(xiàn)為西安交通大學在讀博士生,主要從事地震信號處理和地震屬性分析方面的研究.E-mail:feitianliuhuo@gmail.com
*通訊作者陳文超,男,1970年生,副教授,2000年于西安交通大學獲博士學位,主要從事小波分析及盲信號分離在地震信號處理、解釋中的應(yīng)用研究.E-mail:wencchen@m(xù)ail.xjtu.edu.cn
王偉,高靜懷,陳文超等.基于結(jié)構(gòu)自適應(yīng)中值濾波器的隨機噪聲衰減方法.地球物理學報,2012,55(5):1732-1741,
10.6038/j.issn.0001-5733.2012.05.030.
Wang W,Gao J H,Chen W C,et al.Random seismic noise suppression via structure-adaptive median filter.Chinese J.Geophys.(in Chinese),2012,55(5):1732-1741,doi:10.6038/j.issn.0001-5733.2012.05.030.
(本文編輯 劉少華)