王 福,王華忠

(波現(xiàn)象與智能反演成像研究組(WPI),同濟大學(xué)海洋與地球科學(xué)學(xué)院,上海200092)

精細(xì)的油藏描述首先需要信噪比較高的地震數(shù)據(jù),然后是基于此數(shù)據(jù)的地震波反演成像?！皟蓪捯桓摺钡卣饠?shù)據(jù)采集代表了當(dāng)前地震數(shù)據(jù)采集技術(shù)的發(fā)展方向,基于“兩寬一高”地震數(shù)據(jù)進行地震波反演成像實現(xiàn)目標(biāo)油藏精細(xì)描述嚴(yán)重受限于地震數(shù)據(jù)的信噪比。

“兩寬一高”地震數(shù)據(jù)采集為去噪提供了更好的數(shù)據(jù)基礎(chǔ),把去噪方法推進到高維空間中能更好地把握數(shù)據(jù)中蘊含的信號特征,從而去除噪聲更為徹底。地震信號分析理論中,實際觀測的地震數(shù)據(jù)被視為確定性信號與各種隨機噪聲疊加形成的所謂隨機信號。隨機信號分析是基于隨機信號的概率分布特征。針對離散實測數(shù)據(jù)的信號分析,基于其統(tǒng)計特征,在常見的一維隨機信號的統(tǒng)計分析中,通常假設(shè)隨機信號是平穩(wěn)的,并滿足高斯分布?？紤]到地震信號是蘊含結(jié)構(gòu)特征的高維信號,結(jié)構(gòu)特征的存在破壞了統(tǒng)計去噪方法的平穩(wěn)性假設(shè),因此基于隨機信號統(tǒng)計特征的去噪方法或濾波器設(shè)計要充分體現(xiàn)地震信號的結(jié)構(gòu)特征。各向異性高斯濾波器和中值濾波器設(shè)計,包括各向異性擴散方程濾波方法都是上述思想的具體體現(xiàn)。另外在貝葉斯框架下用信號模型對實測信號進行預(yù)測,預(yù)測誤差假定是高斯的(尤其假定是高斯白噪聲),信號模型是線性的,優(yōu)化求解得到最佳濾波器,從而實現(xiàn)觀測數(shù)據(jù)中噪聲的壓制,是另一種經(jīng)典的地震數(shù)據(jù)去噪思路。在“兩寬一高”地震數(shù)據(jù)采集已成常規(guī)的情況下,貝葉斯框架下的去噪方法同樣亟需推進到高維空間中。針對高維數(shù)據(jù)去噪,上述兩種去噪方法具有逐漸融合的發(fā)展趨勢,本文重點分析基于隨機信號理論的統(tǒng)計去噪思想與方法。

實際觀測的地震數(shù)據(jù)可視為帶限子波形成的確定性信號與隨機噪聲疊加形成的隨機信號,當(dāng)對地震數(shù)據(jù)進行統(tǒng)計濾波時,濾除的隨機噪聲主要有服從高斯分布的高斯白噪聲和服從拉普拉斯分布的脈沖噪聲。針對高斯白噪聲,主要采用統(tǒng)計均值類濾波器。在地震數(shù)據(jù)去噪方法研究中,關(guān)于統(tǒng)計均值類濾波器的研究目前主要集中在如何構(gòu)建沿著信號結(jié)構(gòu)方向的均值類濾波器和如何依據(jù)地震數(shù)據(jù)變化自適應(yīng)地選擇合適的濾波參數(shù)。LUO等[1]提出了保邊平滑濾波方法,在包含當(dāng)前濾波點且方差最小的數(shù)據(jù)窗內(nèi)進行均值濾波。在保邊濾波的基礎(chǔ)上,AN等[2]設(shè)計了最大一致性傾角掃描算法,準(zhǔn)確地拾取地層傾角信息,沿著同相軸進行濾波。范桃園等[3]通過相干技術(shù)和圖像分割技術(shù)考察了濾波點周圍地震數(shù)據(jù)的結(jié)構(gòu)信息,自適應(yīng)確定均值濾波器的形狀。HALE[4]設(shè)計了沿著構(gòu)造方向的雙邊濾波器,對地震圖像進行增強處理。華崗等[5]在雙邊濾波的基礎(chǔ)上提出一種基于地震數(shù)據(jù)局部結(jié)構(gòu)屬性的三邊濾波算法。BONAR等[6]將非局部均值濾波引入到時空域地震資料處理中。SHANG等[7]提出了自適應(yīng)地震數(shù)據(jù)分塊非局部均值算法(ABNLM),通過估算圖像噪聲方差,選擇合適的濾波參數(shù)。針對脈沖噪聲,主要采用統(tǒng)計排序類濾波器。目前在地震數(shù)據(jù)處理中,關(guān)于統(tǒng)計排序類濾波器的研究主要集中在構(gòu)建沿著構(gòu)造方向的中值濾波器及加權(quán)中值濾波器。LIU等[8]在多個方向上應(yīng)用一維中值濾波器,構(gòu)建了自適應(yīng)二維多級中值濾波器,有效降低了地震數(shù)據(jù)中的高頻隨機噪聲。LIU等[9]提出一維自適應(yīng)非平穩(wěn)時變中值濾波器,依據(jù)地震數(shù)據(jù)改變?yōu)V波器的范圍。劉洋等[10]提出了局部相關(guān)加權(quán)中值濾波技術(shù),通過平面波分解濾波器求出局部地層傾角,在數(shù)據(jù)點傾角方向上應(yīng)用加權(quán)中值濾波器。王偉等[11]假設(shè)地震反射同相軸具有局部線性結(jié)構(gòu),通過結(jié)構(gòu)張量構(gòu)造了依據(jù)地層傾向和地層結(jié)構(gòu)自適應(yīng)變化濾波范圍的中值濾波器。LIU[12]將三維中值濾波引入到地震數(shù)據(jù)的處理中。AL-DOSSARY等[13]基于局部噪聲水平設(shè)計了自適應(yīng)濾波范圍的中值濾波器來壓制三維地震數(shù)據(jù)中的隨機噪聲。尋超等[14]采用矢量中值濾波器,通過相關(guān)分析自適應(yīng)地選擇最佳矢量中值濾波窗對三分量地震數(shù)據(jù)進行濾波。

“兩寬一高”地震數(shù)據(jù)是一個典型的五維數(shù)據(jù)體,我們認(rèn)為該地震數(shù)據(jù)是線性和/或非線性的連續(xù)同相軸處在滿足不同統(tǒng)計特征的隨機噪聲中,形成所謂的高維隨機信號,二維以上的地震數(shù)據(jù)都可以認(rèn)為是高維數(shù)據(jù)。這就是基于統(tǒng)計特征進行去噪時所要建立的概念信號模型?；谶@樣的概念信號模型構(gòu)建統(tǒng)計濾波器時,首先需要對數(shù)據(jù)進行加窗處理,接著判別窗內(nèi)數(shù)據(jù)中確定性信號的幾何特征(主要是方向特征),然后基于噪聲的統(tǒng)計特征設(shè)計濾波器(原則上也應(yīng)該在貝葉斯框架下設(shè)計最佳統(tǒng)計濾波器！),最后采用該濾波器進行去噪處理。

本文主要討論統(tǒng)計濾波器的設(shè)計與應(yīng)用。統(tǒng)計濾波的核心思想是根據(jù)觀測數(shù)據(jù)的統(tǒng)計特征設(shè)計合適的統(tǒng)計濾波器。首先討論了地震數(shù)據(jù)中常見隨機噪聲的概率密度函數(shù),重點分析了圖像處理中壓制高斯白噪聲的各向同性高斯濾波器及壓制脈沖噪聲的統(tǒng)計排序類濾波器。各向同性高斯濾波器主要包括均值濾波器和高斯(加權(quán))濾波器,統(tǒng)計排序類濾波器主要有中值濾波器、最大值濾波器、最小值濾波器。指出統(tǒng)計濾波方法是基于信號預(yù)測模型濾波方法的重要補充,在地震數(shù)據(jù)去噪中有著不可替代的作用。由于高維地震數(shù)據(jù)中同相軸的存在破壞了統(tǒng)計濾波的平穩(wěn)性假設(shè),因此,在統(tǒng)計濾波的過程中為了不破壞同相軸(即地震信號)的保真性,必須引入同相軸結(jié)構(gòu)的預(yù)測算子,把各向同性統(tǒng)計濾波器修改成基于同相軸結(jié)構(gòu)信息的各向異性統(tǒng)計濾波器。然后,具體分析了鄰域濾波器、雙邊濾波器、非局部均值濾波器三類各向異性高斯(加權(quán))濾波器的設(shè)計思想,并在非局部均值濾波算法的基礎(chǔ)上,設(shè)計了自適應(yīng)變化搜索窗的非局部均值濾波算法,該方法通過局部數(shù)據(jù)窗的相關(guān)來考察地震數(shù)據(jù)中的結(jié)構(gòu)特征,依據(jù)地震數(shù)據(jù)變化自適應(yīng)地改變非局部均值濾波算法中的搜索窗。最后,用合成地震數(shù)據(jù)驗證了所提方法的有效性。

1 地震數(shù)據(jù)統(tǒng)計濾波器的基本思想

1.1 實測地震數(shù)據(jù)中常見的隨機噪聲

地震信號處理主要包括兩方面的核心內(nèi)容:信號建模和在貝葉斯估計理論下對信號模型中所含參數(shù)的最佳估計。具體的信號分析或處理工作就是對上述兩種內(nèi)容的具體應(yīng)用,譬如去噪、插值、信號恢復(fù)等。比較而言,信號建模更為基礎(chǔ),在希爾伯特空間中構(gòu)建一組(或多組)基函數(shù)形成子空間,在子空間中對信號進行基函數(shù)線性組合的表達(dá)是信號建模的最基本方法或最常用方法。但是實測數(shù)據(jù)通常都是隨機的,針對隨機數(shù)據(jù)(或稱隨機信號)用概率統(tǒng)計的方法進行統(tǒng)計建模,并基于此進行信號分析也是非常普遍的做法。

首先建立如下的概念模型,認(rèn)為實際測量的地震數(shù)據(jù)是確定性信號與滿足各種概率分布的隨機噪聲的疊加,即:

(1)

式中:(x,y,t)表示地震數(shù)據(jù)的空間時間坐標(biāo);u(x,y,t)表示實際測量的地震數(shù)據(jù);s(x,y,t)表示對應(yīng)的確定性信號(在貝葉斯反演框架下,要預(yù)測的信號也是隨機的！);η(x,y,t)是滿足某種概率分布的加性隨機噪聲,并假定各坐標(biāo)點處的噪聲分量為獨立同分布的隨機變量且與該點處的信號s(x,y,t)無關(guān)。在實際地震信號的統(tǒng)計處理中,(1)式的信號模型還可以表達(dá)在常見的五維空間(譬如u(mx,my,hx,hy,t),其中mx,my為中心點坐標(biāo),hx,hy為偏移距坐標(biāo))中,在更高維的空間中進行統(tǒng)計信號分析,統(tǒng)計特征會更符合既定的假設(shè),濾波效果會更好。(1)式定義的信號模型的應(yīng)用非常普遍,當(dāng)前絕大部分地震信號處理方法的構(gòu)建都從此概念模型出發(fā),其核心工作就是實現(xiàn)對信號s(x,y,t)的最佳預(yù)測或最佳(特征)表達(dá),從而達(dá)到最佳去噪、插值、壓縮等目的。

從統(tǒng)計意義下建立一個模型實現(xiàn)對信號的預(yù)測,就是建立(1)式中噪聲滿足的概率分布模型,而聯(lián)合概率密度函數(shù)是刻畫隨機信號特征的最基本形式。在地震勘探中,最基本的隨機噪聲類型主要有:滿足高斯分布的高斯白噪聲、滿足拉普拉斯分布的脈沖噪聲以及滿足均勻分布的隨機噪聲。滿足均勻分布的隨機噪聲具有重要的數(shù)學(xué)意義,但物理世界中應(yīng)該很少有如此類型的噪聲,尤其在勘探地震的實際數(shù)據(jù)中。

高斯白噪聲滿足的概率密度函數(shù)為:

(2)

式中:σ2為η(x,y,t)的方差。實際地震數(shù)據(jù)包含的噪聲往往并非高斯白噪聲,更多的是一般高斯噪聲,甚至是高斯有色噪聲,設(shè)計統(tǒng)計濾波器時這是非常重要的考慮因素。

脈沖噪聲(或稱椒鹽噪聲)的概率密度函數(shù)為:

(3)

若Pa或Pb為零,則脈沖噪聲稱為單極脈沖噪聲;若Pa和Pb都不為零,且近似相等時,則脈沖噪聲類似于在圖像上隨機分布的胡椒和鹽粉微粒,稱為雙極脈沖噪聲,也稱為椒鹽噪聲[15]。此類噪聲在實際地震數(shù)據(jù)中經(jīng)常出現(xiàn):譬如未知源的野值,同時源激發(fā)數(shù)據(jù)按主炮偽解碼后副炮信號表現(xiàn)為脈沖噪聲等。

均勻分布的隨機噪聲滿足的概率密度函數(shù)為:

(4)

式中:Ω代表隨機噪聲η(x,y,t)的取值范圍。服從均勻分布的隨機噪聲在進行理論(數(shù)值)分析時非常有用。

另外,基于隨機噪聲的概率分布特征進行噪聲η(x,y,t)的壓制時,對信號s(x,y,t)是有假定條件的。即假設(shè)信號是局部緩變的,最好是不變的,設(shè)計任何統(tǒng)計濾波器時必須充分注意這個假設(shè)條件。

1.2 基于高斯分布的隨機噪聲壓制的濾波器設(shè)計思想

1.2.1 均值濾波器的設(shè)計思想

(5)

很顯然,該問題的解為:

(6)

由于實際地表觀測的地震信號都是離散的,在離散情形下,最簡單的估計方程為:

(7)

1.2.2 線性高斯加權(quán)濾波器的設(shè)計思想

理論上,實測(隨機)數(shù)據(jù)u(x,y,t)的統(tǒng)計期望值的計算式為:

(8)

當(dāng)概率密度函數(shù)為高斯函數(shù)時,(8)式稱為(各向同性)高斯(平滑)加權(quán)濾波器。針對實測的離散隨機數(shù)據(jù),(8)式可改寫為:

(9)

(9)式含義為:在給定的空間加權(quán)模板上,信號(或圖像)中的每一個離散點(x,y,t)的濾波結(jié)果為空間加權(quán)模板與以該離散點為中心的模板所覆蓋區(qū)域內(nèi)的觀測值的加權(quán)平均。式中:w(i,j,k)為定義的空間加權(quán)模板,其大小為(2a+1)×(2b+1)×(2c+1),當(dāng)w(i,j,k)中所有系數(shù)都相等時,(9)式就退化成(7)式。

一般地,設(shè)計線性高斯加權(quán)濾波器時,理論上應(yīng)該用高斯概率密度函數(shù)來定義加權(quán)函數(shù)w(i,j,k),即:

(10)

式中:μ為期望,被認(rèn)為等于0;σ2為方差,其大小由噪聲水平?jīng)Q定;C為常數(shù)。由于(10)式定義的概率密度函數(shù)是鐘形函數(shù),其積分總和等于1,在離散情形下,(9)式中的加權(quán)系數(shù)w(i,j,k)可以定義為:

(11)

很顯然,加權(quán)系數(shù)w(i,j,k)是“鐘形”的,h控制它的“緊度”。當(dāng)h趨于0時,w(i,j,k)接近于脈沖函數(shù),高斯加權(quán)模板w(i,j,k)的中心點系數(shù)趨于1,其它點趨于0,對當(dāng)前點無濾波作用;h越大,線性高斯加權(quán)濾波器對信號(圖像)的平滑作用越明顯。相比于(7)式定義的均值濾波器,(9)式定義的線性高斯加權(quán)濾波器考慮了濾波模板內(nèi)不同數(shù)據(jù)點的加權(quán)作用。一般地,應(yīng)按照偏離濾波模板中心點的距離在方差的控制下決定加權(quán)系數(shù),更本質(zhì)地,還應(yīng)該按照信號s(x,y,t)的變化特征來構(gòu)造更復(fù)雜的加權(quán)函數(shù),高維數(shù)據(jù)空間中各向異性高斯加權(quán)濾波器據(jù)此產(chǎn)生。這是當(dāng)前統(tǒng)計去噪領(lǐng)域中重要的研究內(nèi)容,也是本文的核心關(guān)注點。

在高斯白噪聲的假設(shè)下,線性高斯加權(quán)濾波器較均值濾波器效果更好,這是統(tǒng)計濾波的理論要求,(8)式清楚地說明了這一點。實際計算中的問題是很難得到噪聲的實際概率密度函數(shù),(10)式實際上是人為給出的,(11)式定義的加權(quán)函數(shù)是在(10)式的啟發(fā)下給出的。針對實測隨機數(shù)據(jù)的統(tǒng)計去噪,本質(zhì)上需要確定噪聲所滿足的概率密度函數(shù),然后通過(5)式定義的均方誤差最小的逼近問題來估計概率密度函數(shù)中的系數(shù),最后利用估計的最佳濾波器進行加權(quán)濾波處理,這才是理論上完善的統(tǒng)計濾波方法,但是數(shù)值計算不易實現(xiàn)。

1.3 基于拉普拉斯分布的隨機噪聲壓制的濾波器設(shè)計思想

在統(tǒng)計濾波器中,除了壓制高斯白噪聲的均值濾波器和線性高斯加權(quán)濾波器,針對服從拉普拉斯分布的脈沖噪聲,主要采用統(tǒng)計排序類濾波器,譬如中值濾波器以及它的變種(最大值濾波器及最小值濾波器)。

拉普拉斯分布的基本形式為:

(12)

式中:α和β分別為拉普拉斯分布的位置參數(shù)和尺度參數(shù)。根據(jù)拉普拉斯分布可以從估計理論的角度來說明統(tǒng)計中值濾波器的設(shè)計思想。

假設(shè)局部信號緩變,設(shè)um(x,y,t)是以(x,y,t)為中心點的局部鄰域Sxyt內(nèi)實測(隨機)數(shù)據(jù){u(xi,yj,tk)|(xi,yj,tk)∈Sxyt}的中位數(shù),同樣地,在某種程度上,它是對其中蘊含的確定性信號s(x,y,t)的估計結(jié)果。在離散數(shù)據(jù)情形下,根據(jù)噪聲滿足拉普拉斯分布的假設(shè),可將上述估計問題表示為絕對值誤差最小意義下的最優(yōu)化問題:

(13)

求解上式,有:

(14)

即:

(15)

式中:sign為符號函數(shù)。當(dāng)um(x,y,t)u(xi,yj,tk)時,sign為正;當(dāng)um(x,y,t)

(16)

其含義為在一個以(x,y,t)為中心的鄰域Sxyt中取N=NxNyNt個數(shù)據(jù)的中位數(shù)來替代點(x,y,t)的值,作為該點的濾波結(jié)果。很明顯,窗口Sxyt中N=NxNyNt個數(shù)據(jù)的中位數(shù)作為濾波結(jié)果并不代表這是最佳的濾波結(jié)果,是否最佳取決于噪聲真實的概率密度函數(shù)。然而實際上很難知道此概率密度函數(shù),根據(jù)經(jīng)驗,可以取中值的百分?jǐn)?shù)作為濾波結(jié)果,可能會得到更好的濾波結(jié)果,不同的百分比選擇可以構(gòu)建不同的變種中值濾波器,譬如最大值濾波器、最小值濾波器等。

很明顯,中值濾波是一種非線性濾波,在平穩(wěn)信號及加性高斯白噪聲的情況下,此時觀測值的分布為正態(tài)分布,觀測值的中位數(shù)與平均值一致,可以代表對一組觀測值中信號的最佳估計。在這種理論假設(shè)下,中值濾波可以壓制高斯白噪聲,但是,這種假設(shè)太苛刻,實際上還是盡量不要寄希望于用這樣的方法壓制高斯白噪聲。當(dāng)噪聲的概率分布符合脈沖分布時,觀測數(shù)據(jù)中存在少數(shù)特別大和/或特別小的觀測值,中位數(shù)可以不受它們的影響,能得到更好的估計結(jié)果um(x,y,t),因此濾除椒鹽噪聲用中值濾波器及其相應(yīng)的變種濾波器是比較好的選擇。

相對于線性高斯加權(quán)濾波器,中值濾波器及其相應(yīng)的變種濾波器對信號突變的保持效果更好,對信號的光滑作用小。但是,即便如此,此類統(tǒng)計排序濾波器還是假定實測數(shù)據(jù)中蘊含的信號盡可能平穩(wěn),否則對信號幅值的改變也是不可接受的。地震數(shù)據(jù)的特征是反映各種波現(xiàn)象的同相軸分布在滿足不同概率分布假設(shè)的隨機噪聲中,因此不能認(rèn)為實測數(shù)據(jù)中的信號是滿足平穩(wěn)性假設(shè)的,尤其在小窗口內(nèi)更是如此。在高維數(shù)據(jù)空間中,首先預(yù)測高維地震數(shù)據(jù)蘊含的結(jié)構(gòu)信息,然后沿著信號的結(jié)構(gòu)方向設(shè)計合理的中值濾波器進行濾波是合理的邏輯選擇。

2 高維數(shù)據(jù)空間中高斯加權(quán)濾波器設(shè)計思想

前已述及,無論是統(tǒng)計均值濾波器、線性高斯加權(quán)濾波器還是中值濾波器的理論設(shè)計中,都假設(shè)實測隨機數(shù)據(jù)中蘊含的信號是平穩(wěn)的、變化緩慢的,最好是不變的。實際上信號的變化可能相當(dāng)劇烈,這就會導(dǎo)致濾波后的信號受到傷害,統(tǒng)計均值濾波器、線性高斯加權(quán)濾波器會對信號進行光滑化處理,中值濾波器會改變信號的幅值。

地震數(shù)據(jù)一定是高維的,尤其在“兩寬一高”地震勘探中,地震去噪一定要在高維空間中進行。高維空間中的地震信號有顯著的結(jié)構(gòu)特征,沿著信號的結(jié)構(gòu)方向設(shè)計統(tǒng)計濾波器,更容易滿足統(tǒng)計濾波器要求信號必須是平穩(wěn)變化的假設(shè),各向異性高斯加權(quán)濾波器和自適應(yīng)信號結(jié)構(gòu)特征的中值濾波器就是基于這樣的想法提出的。

在前述統(tǒng)計濾波器設(shè)計思想的基礎(chǔ)上,高維數(shù)據(jù)空間中自適應(yīng)數(shù)據(jù)結(jié)構(gòu)特征的濾波器設(shè)計的重點是如何檢測出隨機數(shù)據(jù)中蘊含的信號結(jié)構(gòu)特征,目前在地震數(shù)據(jù)中常用的測量結(jié)構(gòu)特征的方法主要包括:①基于空間數(shù)據(jù)相關(guān)的傾角掃描方法[2-3,10,14,18];②結(jié)構(gòu)張量方法[5,11,19-21];③局部平面波分解[10,22-25]。第一種方法的抗噪性強,測量結(jié)果比較穩(wěn)健,相關(guān)方法選擇不好時精度較低;第二種方法精度高、局部性好,但不穩(wěn)健;第三種方法最好,采用局部平面波匹配的方法,穩(wěn)健性和計算效率都有保證。

2.1 保持信號結(jié)構(gòu)的非線性高斯加權(quán)濾波器

通過對數(shù)據(jù)的加權(quán)平均來壓制隨機噪聲是一種重要的去噪方法,但其假設(shè)是數(shù)據(jù)中蘊含的信號(幅值)是緩變的,否則一般的線性高斯加權(quán)類濾波器就會對信號進行較為嚴(yán)重的平滑處理。因此,在高維數(shù)據(jù)空間中,設(shè)計保持信號結(jié)構(gòu)的高斯加權(quán)濾波器變得非常必要,其基本思想是自適應(yīng)信號幅值的變化來調(diào)整高斯加權(quán)濾波器,而不是僅僅按前述的“距離”關(guān)系定義高斯加權(quán)濾波系數(shù),更進一步地,應(yīng)該根據(jù)高維信號(圖像)中蘊含的結(jié)構(gòu)模式來設(shè)計加權(quán)濾波器,當(dāng)然,中值濾波器也應(yīng)該根據(jù)高維信號(圖像)中蘊含的結(jié)構(gòu)模式合理設(shè)計。

2.1.1 鄰域濾波器

鄰域濾波器[16-17]的設(shè)計思想是通過計算鄰域內(nèi)各點數(shù)據(jù)幅值(或圖像灰度值)與濾波中心點數(shù)據(jù)幅值的相似程度確定濾波模板系數(shù)。

針對高維數(shù)據(jù)空間中的濾波中心點x(x=(x,y,t)),其空間鄰域定義為:

(17)

式中:y為鄰域內(nèi)任意一點;d為鄰域半徑。在(17)式定義的鄰域內(nèi)考慮離散信號在幅值上的相似性,設(shè)計出如下鄰域濾波器:

(18)

式中:uNF(x)為點x處的鄰域濾波結(jié)果;wr(y)為依據(jù)幅值相似性定義的高斯加權(quán)系數(shù);C(x)為歸一化參數(shù);h為濾波參數(shù);u(y)為鄰域內(nèi)任一點y的含噪聲數(shù)據(jù)。顯然,鄰域濾波器不僅考慮了空間上的相鄰關(guān)系‖x-y‖

參數(shù)h依據(jù)幅值相似性決定高斯加權(quán)的權(quán)重大小,在當(dāng)前濾波點的鄰域內(nèi),鄰域濾波器的加權(quán)系數(shù)根據(jù)像素點幅值的相似性來決定,對于幅值差異較小的點給予較大的權(quán)重,對于幅值差異較大的點給予較小的權(quán)重,從而自適應(yīng)數(shù)據(jù)幅值的變化來調(diào)整加權(quán)系數(shù),降低了對信號的光滑化程度,能更好地保持信號的結(jié)構(gòu)特征。

2.1.2 雙邊濾波器

事實上,任何空間相鄰的物理信號大多都有更緊密的關(guān)系,即幅值上具有相似性,在鄰域濾波器的基礎(chǔ)上,TOMASI等[26]于1998年提出了雙邊濾波算法,同時考慮了鄰域內(nèi)任一點與濾波中心點之間幅值的相似程度和二者的距離關(guān)系來定義高斯加權(quán)系數(shù)。雙邊濾波器的定義如下:

(19)

式中:uBF(x)為點x處的雙邊濾波結(jié)果;wd(y)為依據(jù)鄰域Ω內(nèi)點y與濾波中心點x之間的距離定義的高斯加權(quán)系數(shù);wr(y)為依據(jù)幅值相似性定義的高斯加權(quán)系數(shù);ρ為空間域濾波參數(shù);h為與幅值有關(guān)的濾波參數(shù);C(x)為歸一化參數(shù)。

從(19)式可以看出,雙邊濾波的加權(quán)系數(shù)同時取決于空間域濾波參數(shù)ρ和幅值相似參數(shù)h。鄰域濾波器和雙邊濾波器在灰度域內(nèi)都進行高斯濾波,不同的是,對于鄰域濾波器,在鄰域內(nèi)的所有像素點空間加權(quán)系數(shù)都為1,即空間域進行均值濾波,對于雙邊濾波器,靠近濾波中心點的像素點空間加權(quán)系數(shù)wd(y)大,即在空間域進行高斯濾波,因此相較鄰域濾波,雙邊濾波器具有更好的保持邊緣特征(或結(jié)構(gòu)特征)的能力。

2.1.3 非局部均值濾波器

無論是鄰域濾波器或是雙邊濾波器的設(shè)計都僅僅是在一個鄰域內(nèi)基于單個像素點距離和/或幅值相似性進行的。當(dāng)前信號分析的原則性思想是利用更高維空間的數(shù)據(jù)中蘊含的信號的結(jié)構(gòu)信息對信號進行更徹底的分析,即使在基于高維數(shù)據(jù)定義的局部鄰域內(nèi),基于單個像素點相似性的濾波器也難以準(zhǔn)確地感知到高維信號的結(jié)構(gòu)(或紋理)特征,有必要基于單個像素點的局部鄰域的相似性設(shè)計加權(quán)濾波器。

依據(jù)上述思想,BUADES等[27]設(shè)計了一種非局部均值濾波器,核心思想是在一個較大搜索窗I(甚至整幅圖像中)對上述單個像素點x的局部鄰域Ω搜索相似的局部鄰域,可認(rèn)為在進行“升維”處理(我們認(rèn)為這是一種升維處理),即將I內(nèi)任一點y都視為一個新的中心點,據(jù)此定義一個新的鄰域Ω′,通過考察鄰域Ω和鄰域Ω′的模式相似性(用圖像處理的語言,是兩個鄰域的模式相似性),定義非局部均值濾波器。

非局部均值濾波器定義為:

(20)

式中:uNLM(x)為點x處的非局部均值濾波結(jié)果;wp(y)為依據(jù)鄰域Ω和鄰域Ω′的模式相似性定義的加權(quán)系數(shù);Z(|uΩ(x)-uΩ′(y)|2)是模式相似性度量函數(shù);h為濾波參數(shù);C(x)為歸一化參數(shù)。

度量兩個鄰域內(nèi)信號(圖像)結(jié)構(gòu)模式相似性的方法有很多,直接對信號結(jié)構(gòu)模式進行檢測就是合乎邏輯的方法,譬如前述的局部傾角掃描法、結(jié)構(gòu)張量方法和平面波分解方法。當(dāng)然也可以根據(jù)鄰域Ω和鄰域Ω′內(nèi)信號(圖像)幅值差的某種范數(shù)來判斷其中包含的信號(圖像)模式是否一致,并根據(jù)一致性的定量評判結(jié)果來設(shè)計非局部均值濾波器。

Z(|uΩ(x)-uΩ′(y)|2)可以定義為:

(21)

(21)式再次引入方差為σ2的高斯函數(shù)Gσ加權(quán)來考察兩個鄰域中模式的相似度,這樣要優(yōu)于不考慮高斯加權(quán)。

相比于鄰域濾波器和雙邊濾波器,非局部均值濾波器以鄰域Ω和鄰域Ω′中的模式相似性替代用孤立的像素點來考慮像素點之間的相似性,具有較強的抗噪能力,可以更加準(zhǔn)確地檢測出信號的結(jié)構(gòu)信息,能夠更好地實現(xiàn)保持信號結(jié)構(gòu)特征的去噪。

2.2 高維數(shù)據(jù)中信號結(jié)構(gòu)導(dǎo)引濾波器設(shè)計的進一步評述

針對高維數(shù)據(jù)空間中的概念信號模型,u(x,y,t)=s(x,y,t)+η(x,y,t),濾除噪聲的基本思想是設(shè)計最佳的信號預(yù)測器壓制噪聲或者設(shè)計最佳的噪聲統(tǒng)計預(yù)測模型從而壓制噪聲。前者是選擇合適的基函數(shù)族(或基函數(shù)字典)進行線性組合并在貝葉斯估計理論下得到最佳的信號預(yù)測器。后者是前面討論的在信號的結(jié)構(gòu)導(dǎo)引下設(shè)計高斯加權(quán)濾波器。實際上,在線性(此處線性的含義是:信號可線性預(yù)測或信號線性變化或信號平緩變化)最小二乘(隱含地指噪聲是高斯分布的)意義下,二者的差異非常小。由于二者的實現(xiàn)方式和控制參數(shù)不同,有必要沿著這兩種方式發(fā)展高維數(shù)據(jù)空間中的濾波器,處理“兩寬一高”數(shù)據(jù)中不同類型的噪聲,尤其是針對偏離高斯分布的噪聲兩種方法都有各自廣闊的發(fā)展空間。

基于擴散方程進行圖像(信號)的濾波也是一類壓制隨機噪聲的方法,擴散方程本質(zhì)上也是一種濾波器,通過局部鄰域加權(quán)平均的思想,即可實現(xiàn)基于擴散方程的去噪。BUADES等[28]和SCHERZER[29]論證了鄰域濾波器與Perona-Malik擴散方程是漸近一致的。利用各向異性擴散方程濾波更容易實現(xiàn)保持結(jié)構(gòu)特征的隨機噪聲壓制,尤其是對于存在結(jié)構(gòu)連續(xù)性的圖像(信號)。各向異性高斯加權(quán)濾波器與各向異性擴散方程濾波本質(zhì)上相同,因為實現(xiàn)方式和參數(shù)控制上的不同使得兩類方法并存。

最后,基于拉普拉斯分布下的L1范數(shù)最小導(dǎo)出的中值濾波及變種方法同樣需要拓展到高維數(shù)據(jù),也同樣有必要在信號結(jié)構(gòu)特征的導(dǎo)引下設(shè)計出更好的濾波器。

3 自適應(yīng)搜索窗的非局部均值濾波器

3.1 自適應(yīng)搜索窗的非局部均值濾波器原理

在地震勘探中,由于地震圖像與常規(guī)圖像存在明顯不同,HALE[4]給出了地震圖像中邊緣產(chǎn)生的兩種原因:①地震數(shù)據(jù)對應(yīng)于波阻抗變化引起的反射波,由于子波帶限,地震數(shù)據(jù)表現(xiàn)為波峰和波谷交替出現(xiàn)的正弦特征,這與常規(guī)圖像邊緣存在明顯差異;②地震反射橫向不連續(xù)性對應(yīng)的邊緣(如斷層等)。對地震數(shù)據(jù)中的高斯白噪聲進行濾波處理時,直接進行均值濾波或線性高斯加權(quán)濾波在平滑噪聲的同時會傷害有效地震信號。

考慮到地震數(shù)據(jù)中蘊含的結(jié)構(gòu)特征,在將上述非局部均值濾波器用于壓制地震記錄中的高斯白噪聲時,存在兩個關(guān)鍵的問題:首先是如何精確地檢測出地震數(shù)據(jù)中蘊含的信號的結(jié)構(gòu)特征,使得濾波器沿著同相軸方向進行濾波;其次是濾波參數(shù)的選取,小濾波參數(shù)不能較好地濾除噪聲,大濾波參數(shù)在濾除噪聲的同時會傷害有效信號,需要考慮如何設(shè)計合適的濾波參數(shù),在濾除隨機噪聲的同時盡可能地保護有效地震信號。

對于如何精確地檢測出地震數(shù)據(jù)中蘊含的信號的結(jié)構(gòu)特征這一問題,主要有基于空間數(shù)據(jù)相關(guān)的傾角掃描方法、結(jié)構(gòu)張量方法和局部平面波分解等方法?？紤]到相關(guān)類方法抗噪性強,測量結(jié)果比較穩(wěn)健,本文選擇采用局部數(shù)據(jù)窗的相關(guān)算法來考察數(shù)據(jù)點之間的結(jié)構(gòu)特征,設(shè)計了自適應(yīng)地震數(shù)據(jù)變化搜索窗的非局部均值濾波器,將(20)式中的固定搜索窗I改成對于數(shù)據(jù)點u(x)自適應(yīng)搜索窗Ix,即:

(22)

式中:uAWNLM為自適應(yīng)搜索窗的非局部均值濾波結(jié)果。圖1為局部窗滑動相關(guān)示意圖,在濾波點處開窗(圖1中藍(lán)色窗),與相鄰地震道的滑動數(shù)據(jù)窗(圖1中綠色窗)進行相關(guān)計算。

相關(guān)系數(shù)計算公式為:

r=

(23)

式中:(x,t)表示當(dāng)前濾波點位于第x道第t個采樣時間;(x′,t′)表示滑動窗口的中心位于第x′道第t′個采樣時間;r表示兩個數(shù)據(jù)窗之間的相關(guān)系數(shù)。進行相關(guān)計算的數(shù)據(jù)窗大小為(2m+1)×(2n+1)。當(dāng)相關(guān)系數(shù)最大時,將此窗(圖1中紅色窗)的中心點以及同地震道的上、下兩點加到自適應(yīng)搜索窗,繼續(xù)進行下一道搜索,當(dāng)最大相關(guān)系數(shù)小于閾值時,停止修改自適應(yīng)搜索窗。

圖1 局部開窗滑動相關(guān)示意

在濾波點(x,t)處開窗對下一道地震數(shù)據(jù)進行相關(guān)計算,尋找結(jié)構(gòu)相似點的具體流程如表1所示,完成該點右邊地震道相關(guān)搜索以后還需進行左邊地震道的計算。

對于當(dāng)前濾波點(x,t),完成自適應(yīng)搜索窗的計算后,當(dāng)搜索窗內(nèi)的點數(shù)大于一定閾值時,可以認(rèn)為該濾波點周圍包含較多相關(guān)信息,計算得到的自適應(yīng)搜索窗是沿著地震數(shù)據(jù)同相軸截取的數(shù)據(jù)窗,在自適應(yīng)搜索窗內(nèi)對該點進行非局部均值濾波,反之,則認(rèn)為該點無結(jié)構(gòu)信息,在該點處取大小固定的搜索窗進行非局部均值濾波。

3.2 數(shù)值算例

為了測試自適應(yīng)搜索窗的非局部均值濾波算法,

本文采用合成地震數(shù)據(jù)進行去噪試驗(圖2)。合成地震記錄的子波主頻為40Hz,時間采樣間隔為2ms,共41道,如圖2a所示。圖2b為加入-1dB高斯白噪聲的數(shù)據(jù)。圖2c為采用21×21的固定搜索窗,7×7的鄰域窗,濾波參數(shù)h=2.5,σ為1的非局部均值濾波結(jié)果,可以看到,非局部均值濾波雖然可以剔除一定噪聲,但也損失了部分有效信號。圖2d為采用自適應(yīng)搜索窗的非局部均值濾波結(jié)果,在該算法(表1)中,局部相關(guān)的數(shù)據(jù)窗大小為11×3,相關(guān)系數(shù)閾值設(shè)為0.35,搜索窗內(nèi)的點數(shù)閾值設(shè)為9,固定搜索窗的大小為21×21,鄰域窗的大小為7×7,濾波參數(shù)h=2.5,σ為1。從圖2d中可以看到采用自適應(yīng)搜索窗的非局部均值濾波算法在濾除一定噪聲的同時也能保護有效地震信號。取濾波參數(shù)h=5.0,其余參數(shù)與圖2d的濾波參數(shù)相同,可得到圖2e所示的自適應(yīng)搜索窗的非局部均值濾波結(jié)果。圖2f為圖2c 所示的濾波結(jié)果與合成地震數(shù)據(jù)之間的差剖面,從中可以明顯看到殘留的同相軸信息。圖2g為圖2d 所示的濾波結(jié)果與合成地震數(shù)據(jù)之間的差剖面,與圖2f相比,其中殘留的同相軸信息極少。圖2h 為圖2e所示的濾波結(jié)果與合成地震數(shù)據(jù)之間的差剖面。對比圖2d和圖2e 可知,隨著濾波參數(shù)的增大,自適應(yīng)搜索窗的非局部均值濾波算法能在壓制更多隨機噪聲的同時較好地保護有效地震信號。

圖2 合成地震數(shù)據(jù)去噪效果a 原始記錄; b 含-1dB高斯白噪地震數(shù)據(jù); c 濾波參數(shù)h=2.5的非局部均值濾波結(jié)果; d 濾波參數(shù)h=2.5的自適應(yīng)搜索窗非局部均值濾波結(jié)果; e 濾波參數(shù)h=5.0的自適應(yīng)搜索窗非局部均值濾波結(jié)果;f,g,h分別為與c,d,e對應(yīng)的濾波結(jié)果與合成地震數(shù)據(jù)之間的差剖面

上述數(shù)值試驗結(jié)果表明,自適應(yīng)搜索窗的非局部均值濾波算法在壓制隨機噪聲和保護有效地震信號上的表現(xiàn)優(yōu)于固定搜索窗的非局部均值濾波算法。

4 結(jié)論與討論

“兩寬一高”地震數(shù)據(jù)采集越來越普及,但地震數(shù)據(jù)預(yù)處理和地震波反演成像的方法技術(shù)滯后于數(shù)據(jù)采集技術(shù)的發(fā)展。海量數(shù)據(jù)中蘊含的與地下地質(zhì)結(jié)構(gòu)和巖性參數(shù)相關(guān)的信息不能被充分提取出來,信噪比低是關(guān)鍵的制約因素。與“兩寬一高”地震數(shù)據(jù)采集技術(shù)一同進步的首先應(yīng)該是在高維空間中在新的理論框架下提出的各種噪聲壓制技術(shù),盡可能滿足后續(xù)反演成像方法對噪聲的假設(shè)及對信噪比的要求。

高維地震數(shù)據(jù)去噪的核心在于對表達(dá)在不同域中的地震數(shù)據(jù)所蘊含的信號特征有透徹的理解,并建立合適的表達(dá)、預(yù)測及描述方法。對隨機信號而言,隨機信號的統(tǒng)計特征是設(shè)計合適的濾波器的基礎(chǔ)。在高維地震數(shù)據(jù)的去噪過程中,由于都需要在已知信號結(jié)構(gòu)特征的基礎(chǔ)上進行最佳濾波器的設(shè)計,因此基于線性信號預(yù)測器與貝葉斯反演理論的最佳濾波(包括插值)方法與基于高維隨機信號統(tǒng)計特征的濾波方法正在走向融合。

本文系統(tǒng)地分析了在隨機噪聲服從高斯分布和拉普拉斯分布下統(tǒng)計濾波器設(shè)計的基本思想,并重點討論了高維數(shù)據(jù)情形下,保持信號結(jié)構(gòu)特征的各種各向異性高斯加權(quán)濾波器的設(shè)計思想。考慮到高維地震數(shù)據(jù)的特點,設(shè)計了自適應(yīng)搜索窗的非局部均值濾波器,并通過數(shù)值試驗驗證了該方法在濾除高維地震數(shù)據(jù)中高斯白噪聲的有效性。針對高維地震數(shù)據(jù)中復(fù)雜的信號變化和多變的噪聲類型,該類方法還有較大的發(fā)展空間。

目前,統(tǒng)計濾波器已廣泛用于壓制地震數(shù)據(jù)中的隨機噪聲,與各種地震信號結(jié)構(gòu)提取方法結(jié)合實現(xiàn)了沿著信號結(jié)構(gòu)方向的濾波。但是,在高維數(shù)據(jù)空間中進行各向異性統(tǒng)計濾波器的最佳設(shè)計目前并沒有有效的方法,也缺乏真正保持邊界特征的最佳濾波方法,缺乏自動適應(yīng)隨機信號統(tǒng)計特征變化的濾波器設(shè)計方法。筆者認(rèn)為基于隨機信號統(tǒng)計特征的高維濾波方法在壓制隨機噪聲(包括不相干噪聲)方面有明顯的優(yōu)點,在“兩寬一高”地震勘探中需要深入研究此類方法,進一步提高地震數(shù)據(jù)的信噪比。

致謝:感謝中國石油天然氣股份有限公司勘探開發(fā)研究院、西北分院,中海石油(中國)有限公司北京研究中心、湛江分公司以及中國石油化工股份有限公司石油物探技術(shù)研究院、勝利油田分公司對波現(xiàn)象與智能反演成像研究組(WPI)研究工作的資助與支持。