宋肖楠 王華忠 楊 鍇

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

0 引言

在地震數(shù)據(jù)處理中，地震噪聲壓制技術(shù)占據(jù)重要的地位[1]。目前，噪聲壓制方法主要有兩種：一是在Bayes估計理論框架下，以信號滿足線性假設(shè)、噪聲滿足高斯分布為基礎(chǔ)，求解最佳的濾波器系數(shù)，形成最佳預(yù)測濾波器，實現(xiàn)對信號的最佳預(yù)測，如基于AR模型的f-x濾波、Wiener中心預(yù)測濾波等[2-5]；另一種是在已知噪聲統(tǒng)計特征的基礎(chǔ)上，實現(xiàn)對噪聲的精準(zhǔn)預(yù)測，從而間接實現(xiàn)有效信號的恢復(fù)[6-7]。這兩種方法在噪聲壓制中均有較為廣泛的應(yīng)用。在圖像處理中，均值類濾波是壓制高斯白噪的有效途徑，但無法完整保留圖像中的有效信號。

針對該問題，以高斯濾波器為基礎(chǔ)，各向異性濾波器從保護(hù)圖像邊緣和細(xì)節(jié)信息的角度完善濾波器的性能[8-10]。Yaroslavsky[11]以Lee[12]提出的鄰域濾波為基礎(chǔ)，把鄰域內(nèi)各個數(shù)據(jù)點(diǎn)幅值相似程度作為濾波時的加權(quán)系數(shù)。假如濾波中心點(diǎn)即是噪聲點(diǎn)，最好的噪聲壓制方式是將該點(diǎn)的幅值與鄰域內(nèi)所有噪聲點(diǎn)的幅值進(jìn)行加權(quán)平均。該過程需排除有效信息點(diǎn)，僅對噪聲點(diǎn)進(jìn)行加權(quán)平均。

鄰域濾波器以幅值為約束，在一定程度上改善了濾波效果，但仍需發(fā)展更加有效的保護(hù)邊緣信息的方法。Smith等[13]和Tomasi等[14]在鄰域濾波器的基礎(chǔ)上，提出兼顧?quán)徲騼?nèi)數(shù)據(jù)點(diǎn)幅值相似程度和距離因素的雙邊濾波器。Paris等[15]闡述了雙邊濾波器的理論并展示了其應(yīng)用效果。鄰域濾波器與雙邊濾波器的區(qū)別在于處理空間分量的方式不同：對于鄰域濾波器，在鄰域空間內(nèi)的所有數(shù)據(jù)點(diǎn)均會被統(tǒng)一均值處理；而對于雙邊濾波器，更接近濾波中心點(diǎn)的數(shù)據(jù)點(diǎn)被賦予了更大的權(quán)重。

Buades等[16-19]進(jìn)一步拓展了鄰域濾波器，提出非局部均值算法，該算法使用兩個數(shù)據(jù)點(diǎn)對應(yīng)的鄰域間的相似程度評估這兩點(diǎn)間的相似程度。

目前，已有很多保持圖像邊緣信息和結(jié)構(gòu)的濾波方法，比如改進(jìn)的數(shù)字圖像均值濾波算法[20]、改進(jìn)的反擴(kuò)散圖像平滑模型和自適應(yīng)數(shù)值統(tǒng)計濾波算法[21]、改進(jìn)的非局部均值圖像去噪算法[22]、基于統(tǒng)計濾波的自適應(yīng)雙閾值的邊緣檢測算法[23]、自適應(yīng)特征選擇的相關(guān)濾波跟蹤算法[24]和保持圖像邊緣信息的算法[25]等。

地震數(shù)據(jù)與常規(guī)圖像之間存在很大的差異。地震數(shù)據(jù)本質(zhì)上是由滿足一定時距關(guān)系的地震同相軸與具有一定統(tǒng)計規(guī)律的隨機(jī)噪聲構(gòu)成。地震記錄是由地震子波構(gòu)成的[26]，在地震記錄上，地震子波的波峰和波谷交替出現(xiàn)，幅值有正有負(fù)，使得同相軸具有一定寬度。常規(guī)圖像可由二維函數(shù)f(x,y)定義，圖像中每一個空間位置點(diǎn)(x,y)對應(yīng)一個幅值f，即圖像的灰度或強(qiáng)度[26-28]。在地震數(shù)據(jù)處理中，使用經(jīng)典統(tǒng)計濾波器會影響地震數(shù)據(jù)中的同相軸信息，濾波效果不理想。楊葆軍[29]和Chuai等[30]提出利用相干值掃描方法確定構(gòu)造方向，將保結(jié)構(gòu)的雙邊濾波方法應(yīng)用于地震數(shù)據(jù)處理，這在一定程度上考慮了同相軸信息。考慮到局部的地震數(shù)據(jù)特征，基于地震波相似性，提出了自適應(yīng)濾波器等非局部均值濾波[31-33]。為進(jìn)一步保護(hù)同相軸信息，提出了自適應(yīng)順序統(tǒng)計濾波器和加權(quán)中值濾波器相結(jié)合的方法[34]、最大一致性傾角掃描方法[35]、自適應(yīng)地震保邊界方法[36]、基于非局部中值濾波的地震圖像結(jié)構(gòu)定向濾波方法[37]、改進(jìn)的各向異性擴(kuò)散濾波方法[38-39]、基于迭代啟發(fā)網(wǎng)絡(luò)算法的非平穩(wěn)隨機(jī)噪聲壓制方法[40]等。這些方法均可在壓制噪聲的同時，在一定程度上保留了地震同相軸局部變化信息。

在前人研究的基礎(chǔ)上，本文設(shè)計了一種沿同相軸方向的各向異性濾波器。首先，利用相關(guān)系數(shù)提取同一同相軸相鄰道的相似特征。然后，利用縱、橫向方差比值劃分區(qū)域，使濾波器在不同區(qū)域自適應(yīng)地沿同相軸方向濾波，從而設(shè)計適合地震數(shù)據(jù)規(guī)律的濾波器?；诤铣蓴?shù)據(jù)和實際數(shù)據(jù)，驗證了該方法去噪的有效性并保持原有地震信號同相軸的有效信息。此外，本文對比了所提方法與圖像處理中的各向同性濾波器(均值濾波器、高斯濾波器)、各向異性濾波器(鄰域濾波器、雙邊濾波器、非局部均值濾波器)在地震數(shù)據(jù)處理中的優(yōu)缺點(diǎn)。最后，根據(jù)信噪比、殘差和誤差曲線的定量、定性結(jié)果，分析了統(tǒng)計濾波在地震數(shù)據(jù)處理中應(yīng)用的可能性和局限性。

1 方法原理

1.1 噪聲模型

地震數(shù)據(jù)可看作是確定性信號與滿足各種概率分布的隨機(jī)噪聲的疊加[8-9]

u(x,y)=s(x,y)+η(x,y)

(1)

式中：u(x,y)表示觀測數(shù)據(jù)的幅值，其中(x,y)表示第x道、第y個樣點(diǎn)位置；s(x,y)表示地震信號；η(x,y)表示隨機(jī)噪聲。常見的噪聲模型是高斯白噪聲模型[8]。高斯分布是簡單的、大數(shù)定律定義的有限分布。高斯白噪聲具有的特點(diǎn)是功率譜密度服從均勻分布，幅值服從高斯分布。高斯隨機(jī)變量n的概率密度函數(shù)[4]為

(2)

圖1 均值為標(biāo)準(zhǔn)差為σ的高斯概率密度函數(shù)

1.2 沿同相軸方向的各向異性濾波器

統(tǒng)計濾波的噪聲壓制思想是基于實測數(shù)據(jù)的統(tǒng)計特征，利用施加高斯白噪后的數(shù)據(jù)統(tǒng)計量信息，設(shè)計合適的統(tǒng)計濾波器。對含噪圖像施加由濾波窗和參數(shù)確定的濾波算子，設(shè)計濾波器并濾波，使地震信號接近理論信號，并使隨機(jī)噪聲的統(tǒng)計量(期望、方差等)等于或接近零。統(tǒng)計量是隨機(jī)數(shù)據(jù)的某種加權(quán)疊加，其結(jié)果依賴于施加高斯白噪后的數(shù)據(jù)的統(tǒng)計特征。高斯分布是一種大數(shù)定理定義的極限分布。

本文假設(shè)數(shù)據(jù)中包含的噪聲是高斯白噪。地震數(shù)據(jù)是由帶限的振蕩子波組合而成，直接對其進(jìn)行簡單地加權(quán)求和處理會導(dǎo)致子波畸變。為此，設(shè)計沿同相軸方向的各向異性濾波器，在同相軸的方向進(jìn)行加權(quán)疊加，盡可能地保留原始有效信號。該濾波器具有以下幾點(diǎn)特征：沿著同相軸方向濾波； 濾波窗是扁平的； 鄰域窗的形狀和大小需根據(jù)地震信號特征變化進(jìn)行設(shè)計。因此，在理論上若能沿著同相軸方向?qū)Φ卣鹦盘栠M(jìn)行濾波，濾波算子應(yīng)有如下的性質(zhì)

(3)

式中：F是濾波算子；E(·)表示求數(shù)學(xué)期望。因此，針對不同的結(jié)構(gòu)特征，采取不同的濾波策略。理論上，只要沿著同相軸方向濾波，濾波后可以保護(hù)原始的同相軸信息。

本文通過計算鄰道的相關(guān)系數(shù)求取局部的同相軸信息。兩道的相關(guān)系數(shù)計算公式為

r(x,y,x′,y′)=

(4)

式中：(x′,y′)表示當(dāng)前鄰域窗內(nèi)的某一點(diǎn)；a∈[-m,m]，其中m是相關(guān)窗長度的一半。

獲得相關(guān)系數(shù)后，通過鄰域窗振幅和、相似系數(shù)范圍、相似點(diǎn)總數(shù)三個參數(shù)表征數(shù)據(jù)之間的相似程度。其中，振幅和閾值用來排除不含同相軸信息的濾波點(diǎn)，相似系數(shù)范圍和相似點(diǎn)總數(shù)用來提取相鄰?fù)噍S的相似特征。最后，根據(jù)同一同相軸相鄰道的相似程度確定同相軸的方向。

由于地震子波具有一定的延續(xù)長度，鄰域濾波窗應(yīng)該設(shè)計成扁平狀。子波寬度決定濾波器寬度，同相軸的平緩程度決定濾波器的長度。由于地震剖面往往包含噪聲、弱反射、斷層、不整合面等，因此有必要針對不同地質(zhì)情況，構(gòu)建不同尺寸的濾波器。本文通過縱、橫向方差確定濾波器橫向長度。方差縱橫比定義為

(5)

選用幾個不同尺寸的鄰域窗，分別計算縱、橫向方差比，根據(jù)比值粗略地劃分不同的區(qū)域。最后，自適應(yīng)地在不同區(qū)域采用不同尺寸的鄰域窗進(jìn)行濾波。

假設(shè)噪聲是高斯白噪，通過綜合衡量地震數(shù)據(jù)的相似程度，構(gòu)建符合地震數(shù)據(jù)特征的濾波器，利用相關(guān)系數(shù)提取同一同相軸相鄰道的相似特征，使得濾波器在不同區(qū)域自適應(yīng)地沿同相軸方向濾波。本文提出的方法可以考慮地震數(shù)據(jù)中信號的變化特征，因此，可獲得較好的濾波效果。圖2為本文所提統(tǒng)計濾波方法處理流程，可概括為：首先，計算地震數(shù)據(jù)的縱、橫向方差比，劃分不同的區(qū)域。然后，在某個區(qū)域的某個濾波點(diǎn)對應(yīng)的鄰域內(nèi)，如果鄰域內(nèi)幅值和的絕對值小于閾值1，那么認(rèn)為該濾波點(diǎn)是純噪聲點(diǎn)，進(jìn)行均值濾波； 反之，認(rèn)為該鄰域內(nèi)含有同相軸信息。對于含有同相軸信息的鄰域，如果相似系數(shù)大于閾值2且相似點(diǎn)總數(shù)大于閾值3，則沿根據(jù)相似程度確定的同相軸方向進(jìn)行均值濾波； 反之，認(rèn)為該鄰域內(nèi)的大部分點(diǎn)為噪聲點(diǎn)，直接進(jìn)行均值濾波。

圖2 統(tǒng)計濾波方法處理流程

2 數(shù)值實驗

2.1 水平同相軸

采用理論合成數(shù)據(jù)從定性和定量角度分別驗證沿同相軸方向的各向異性濾波器壓制噪聲的可行性和有效性。圖3a為采用主頻為30Hz的雷克子波合成的地震記錄，該數(shù)據(jù)橫向空間樣點(diǎn)數(shù)為40，縱向時間樣點(diǎn)數(shù)為300，時間采樣間隔和空間采樣間隔分別為2ms和1m。選取信噪比為6.94dB的水平同相軸數(shù)據(jù)，沿同相軸方向分別使用均值濾波器、高斯濾波器和鄰域濾波器濾波，結(jié)果和相應(yīng)殘差如圖3c～圖3e所示，可見，殘差中幾乎不含有效信號，三種濾波方式均能有效去除噪聲并保留有效信號。進(jìn)一步計算三種方法處理結(jié)果的信噪比，分別為24.54、24.54、24.56dB。結(jié)果表明，本文方法在理論上是可行的； 只要滿足濾波器沿同相軸方向，分別和均值、高斯、鄰域濾波器相結(jié)合均可獲得很好的濾波效果，因此，本文方法僅采用沿同相軸方向均值濾波器。

圖3 對水平同相軸數(shù)據(jù)三種方法的濾波結(jié)果及殘差對比

對于信噪比為6.94、-2.07dB的水平同相軸數(shù)據(jù)沿同相軸方向做均值濾波，計算每個樣點(diǎn)的相對誤差，部分結(jié)果如圖4所示。圖4從定量角度表明，要使濾波結(jié)果相對誤差降低到1%以下，對于信噪比為6.94dB數(shù)據(jù)，濾波器包含的樣點(diǎn)長度至少達(dá)到80個； 對于信噪比為-2.07dB的數(shù)據(jù)，濾波器包含的樣點(diǎn)長度至少達(dá)到600個。因此信噪比越高，濾波器可以越短，對點(diǎn)數(shù)的限制越小。

圖4 對不同信噪比的合成水平同相軸數(shù)據(jù)均值濾波器結(jié)果的相對誤差曲線

統(tǒng)計濾波器可以達(dá)到理想濾波效果，即徹底壓制地震數(shù)據(jù)中的隨機(jī)噪聲，但應(yīng)用時需要盡量滿足假設(shè)條件。而實際地震數(shù)據(jù)中的同相軸基本上不呈簡單的直線型，數(shù)據(jù)更為復(fù)雜。因此，在實際數(shù)據(jù)應(yīng)用中存在很大困難和局限性。

2.2 雙曲型同相軸

采用理論雙曲型同相軸合成數(shù)據(jù)(圖5a)，驗證沿同相軸方向的各向異性濾波器在壓制噪聲和保持地震數(shù)據(jù)同相軸中有效信息的優(yōu)勢。圖5b是添加高斯白噪后的地震記錄，信噪比為6.66dB。采用均值濾波器(圖5c)、高斯濾波器(圖5d)、鄰域濾波器(圖5e)、雙邊濾波器(圖5f)、非局部鄰域濾波器(圖5g)進(jìn)行濾波，濾波結(jié)果的信噪比分別為11.81、8.86、12.09、12.10、4.91dB； 采用鄰域窗(37×37)和相關(guān)窗(17×1)濾波后地震記錄(圖5h)的信噪比為20.97dB ?？梢?，采用本文方法可有效去除噪聲并保留有效信號，剖面的信噪比得到有效提升。圖6為圖5各濾波結(jié)果的殘差，本文方法濾波結(jié)果的殘差中幾乎不含有效信號(圖6f)，證明了該方法在保持同相軸中的有效信息里具有優(yōu)勢。

圖5 不同濾波器的濾波結(jié)果比較

圖6 不同濾波器濾波后的殘差對比

2.3 Sigsbee 2A模型

采用較復(fù)雜的理論模型Sigsbee 2A測試本文方法的有效性和自適應(yīng)分區(qū)域窗的必要性。圖7a為信噪比為6.63dB的部分Sigsbee 2A模型的成像剖面。模型的橫向和縱向的采樣點(diǎn)數(shù)均為541×261，縱、橫向采樣間隔均為10m。采用本文所提計算不同大小窗的方法，根據(jù)縱橫向方差比劃分的區(qū)域如圖7b所示。圖7c為采用固定窗的濾波結(jié)果，圖7d 為圖7c與原始成像剖面的殘差，圖7e為自適應(yīng)分區(qū)域窗的濾波結(jié)果，圖7f為圖7e與原始成像剖面的殘差。對比圖7c和圖7e發(fā)現(xiàn)，自適應(yīng)地在不同區(qū)域采取不同的鄰域窗參數(shù)進(jìn)行濾波，可以考慮成像剖面中信號的變化特征，提高有效信號的保真度，獲得較好的濾波結(jié)果。若不考慮信號特征，那么固定窗濾波能夠有效地壓制噪音，但殘差中仍有較多的有效信息。

圖7 部分Sigsbee 2A模型數(shù)據(jù)的濾波結(jié)果

2.4 實際數(shù)據(jù)

采用某地區(qū)實際數(shù)據(jù)測試本文方法的應(yīng)用效果。偏移成像剖面如圖8a所示。在方差比分選出的不同區(qū)域中分別濾波得到成像剖面(圖8b)，數(shù)據(jù)殘差剖面(圖8c)中無明顯的連續(xù)同相軸，在壓制噪聲的同時較好地保留了有效信號?？梢姡捎帽疚姆椒奢^好地壓制噪聲，在保留原始同相軸和斷層信息的同時，可提高同相軸的連續(xù)性。

圖8 實際數(shù)據(jù)濾波結(jié)果

3 結(jié)束語

本文提出沿同相軸方向的各向異性濾波器，使濾波器沿同相軸方向在不同區(qū)域自適應(yīng)地濾波。首先，濾波需從統(tǒng)計方面考慮； 其次，統(tǒng)計濾波需跟隨地震同相軸，自適應(yīng)地設(shè)計與不同區(qū)域地質(zhì)特征相匹配的扁平的窗函數(shù)。通過理論數(shù)據(jù)和實際數(shù)據(jù)測試，驗證了本文所提方法可較好地壓制噪聲，保留的同相軸具有較好的相干性和連續(xù)性； 最后，與經(jīng)典的圖像濾波方法相比，本文所提方法在地震數(shù)據(jù)濾波中適用性更強(qiáng)。但在實際的數(shù)據(jù)處理中，統(tǒng)計濾波器的長度和噪音的分布情況需要進(jìn)一步分析和測試。

感謝中石油勘探開發(fā)研究院及西北分院、中海油研究院和湛江分公司、中石化物探技術(shù)研究院和勝利油田分公司對波現(xiàn)象與智能反演成像研究組(WPI)研究工作的資助與支持。