改進的統計濾波方法在地震數據處理中的應用及局限性

宋肖楠 王華忠 楊 鍇

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

0 引言

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

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

鄰域濾波器以幅值為約束，在一定程度上改善了濾波效果，但仍需發展更加有效的保護邊緣信息的方法。Smith等[13]和Tomasi等[14]在鄰域濾波器的基礎上，提出兼顧鄰域內數據點幅值相似程度和距離因素的雙邊濾波器。Paris等[15]闡述了雙邊濾波器的理論并展示了其應用效果。鄰域濾波器與雙邊濾波器的區別在于處理空間分量的方式不同：對于鄰域濾波器，在鄰域空間內的所有數據點均會被統一均值處理；而對于雙邊濾波器，更接近濾波中心點的數據點被賦予了更大的權重。

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

目前，已有很多保持圖像邊緣信息和結構的濾波方法，比如改進的數字圖像均值濾波算法[20]、改進的反擴散圖像平滑模型和自適應數值統計濾波算法[21]、改進的非局部均值圖像去噪算法[22]、基于統計濾波的自適應雙閾值的邊緣檢測算法[23]、自適應特征選擇的相關濾波跟蹤算法[24]和保持圖像邊緣信息的算法[25]等。

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

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

1 方法原理

1.1 噪聲模型

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

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

(1)

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

(2)

圖1 均值為標準差為σ的高斯概率密度函數

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

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

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

(3)

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

本文通過計算鄰道的相關系數求取局部的同相軸信息。兩道的相關系數計算公式為

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

(4)

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

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

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

(5)

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

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

圖2 統計濾波方法處理流程

2 數值實驗

2.1 水平同相軸

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

圖3 對水平同相軸數據三種方法的濾波結果及殘差對比

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

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

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

2.2 雙曲型同相軸

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

圖5 不同濾波器的濾波結果比較

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

2.3 Sigsbee 2A模型

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

圖7 部分Sigsbee 2A模型數據的濾波結果

2.4 實際數據

采用某地區實際數據測試本文方法的應用效果。偏移成像剖面如圖8a所示。在方差比分選出的不同區域中分別濾波得到成像剖面(圖8b)，數據殘差剖面(圖8c)中無明顯的連續同相軸，在壓制噪聲的同時較好地保留了有效信號。可見，采用本文方法可較好地壓制噪聲，在保留原始同相軸和斷層信息的同時，可提高同相軸的連續性。

圖8 實際數據濾波結果

3 結束語

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

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