利用波原子分解系數(shù)自適應(yīng)Wiener濾波方法壓制地震數(shù)據(jù)隨機噪聲

劉彥萍 張乃祿 仵 杰 嚴(yán)正國 高建申

(①西安石油大學(xué)電子工程學(xué)院，陜西西安 710065；②陜西省油氣井測控技術(shù)重點實驗室，陜西西安 710065)

0 引言

實際地震勘探中采集的數(shù)據(jù)包含多種干擾波，有效信號受其影響而被扭曲、畸變甚至湮沒。多年來，為了消除干擾噪聲以突出有效信號，人們展開一系列探究并提出許多行之有效的方法。

隨機噪聲是地震數(shù)據(jù)中最常見的背景噪聲，呈高斯性、平穩(wěn)性[1-3]分布。隨機噪聲壓制方法的研究從未停滯，從20世紀(jì)小波濾波、F-X預(yù)測濾波、K-L(Karhunen-Loeve)變換等方法的興起，到21世紀(jì)被廣泛應(yīng)用于地震數(shù)據(jù)處理領(lǐng)域，并被不斷改進，同時又推出許多新方法，且大多取得較顯著成效。

這些方法常用的有：通過估計信號特征設(shè)計的濾波方法，如F-X預(yù)測濾波[4-6]。該方法能很好地提取有效信號特征，但由于在實際應(yīng)用中需假設(shè)一些條件并存在不可預(yù)測的因素，因此該方法所能達(dá)到的濾波效果很有限。對于同相軸為直線型或近似直線型，且相鄰兩地震道的間距不變的地震記錄，該方法具有較好的預(yù)測效果，但對于含有曲線型同相軸的記錄，其去噪效果較差。

有利用相關(guān)性增強有效信號的濾波方法，如K-L變換法[7-8]、SVD(Singular value decomposition)濾波[9]等。K-L變換是基于統(tǒng)計特性的一種變換，其突出優(yōu)點是去相關(guān)性好，在均方誤差意義下具有最佳性能。但使用該方法需首先知道信源的協(xié)方差矩陣并求出其特征值，而這些在維數(shù)較高時不易獲知和求解。然后，通過保留相關(guān)性較強的信號特征值，舍棄相關(guān)性較小的噪聲特征值，再進行重構(gòu)以實現(xiàn)信噪分離。實際應(yīng)用中，該方法對水平或傾角較小的同相軸增強效果明顯，而對于傾角較大的同相軸跟蹤能力不足，因此增強效果很有限。SVD濾波通過對含噪信號分解后所得特征值矩陣中奇異值的合理劃分，將表征有效信號的奇異值保留用以重構(gòu)，得到去噪后的信號。由于奇異值矩陣中的元素值是從大到小排列的，且奇異值減小得很快，一般選取前10%甚至1%的元素即可捕捉有效信號的特征。但該方法分解出的矩陣解釋性不夠強，奇異值取值個數(shù)對計算精度有較大影響。因此，單純采用SVD濾波往往難以得到理想去噪效果。

有基于時頻變換的濾波方法，如基于WVD(Wigner-Ville分布)和PWVD(Pseudo Wigner-Ville分布)的時頻峰值濾波[10-11]。PWVD是WVD的加窗形式，更適用于非線性、非平穩(wěn)信號的處理。在實際應(yīng)用中，時頻峰值濾波法多采用PWVD，以盡可能滿足信號局部線性化的無偏估計條件。因此，也使該方法存在窗長選擇與噪聲壓制之間的矛盾[12-16]。

還有基于多尺度分析的濾波方法，如小波濾波[17-20]、曲波濾波[21-23]、Shearlet濾波[24-25]等。該類方法通過對含噪數(shù)據(jù)進行多尺度分解，根據(jù)有效信號和噪聲分布在不同尺度系數(shù)上的特點，采用閾值劃分等手段實現(xiàn)信噪分離。但這些方法存在無方向性，或不能達(dá)到最優(yōu)的稀疏逼近，或?qū)﹂撝岛瘮?shù)的選擇較苛刻等不足，其濾波效果仍有較大提升空間。

波原子變換是由Demanet等[26]提出的一種新型多尺度幾何分析工具，可看成是二維小波包變換的變體。該變換具有良好的方向特性，每個波包的振動周期和支撐尺寸滿足拋物尺度關(guān)系，即波長約等于支撐尺寸的平方[26-29]，在此意義下，可簡單地將波原子理解為方向小波與Gabor原子的插值。

波原子變換對于紋理模型具有最優(yōu)的稀疏表示。一條地震記錄可看作為一幅紋理圖像，利用波原子變換可以更好地捕捉其特征[27，29]。在此基礎(chǔ)上，需采用適宜的劃分方法使波原子分解系數(shù)中的信號成分與噪聲成分有效分離，舍棄噪聲系數(shù)，然后通過波原子重構(gòu)得到去噪后的地震記錄。

本文采用對波原子分解所得系數(shù)進行自適應(yīng)Wiener濾波[30-36]的方案，設(shè)計出性能更優(yōu)且易于實現(xiàn)的濾波方法壓制地震勘探隨機噪聲以突出有效信號。通過實驗驗證，該方法較之于閾值分離方法可達(dá)到更優(yōu)的效果。

1 方法原理

1.1 波原子變換原理

相比于小波變換、曲波變換等，波原子變換對具有振蕩性的函數(shù)或具有豐富方向紋理特征的信號能呈現(xiàn)最優(yōu)稀疏表示[26，28]。定義波原子為φμ(x)，μ=(j,m,n)，其中m=(m1,m2)，n=(n1,n2)。五個參量j、m1、m2、n1、n2是整數(shù)，且有j≥0，m≥0,n∈Z。定義相空間中一點(xμ,ωμ)為

xμ=2-jnωμ=π2jm

(1)

式中：C1、C2是兩個正常數(shù)，根據(jù)實際應(yīng)用情況進行合理選擇，可都取為1；位置向量xμ和波向量ωμ分別表征φμ(x)的空域和頻域中心。

波包{φμ}的框架元稱為波原子，波原子圍繞相空間點(xμ,ωμ)遵循局部化條件。即對于任意M>0，有

(2)

(3)

對于任意f(x)∈L2(R2)，在2-j尺度上的空間域一維波原子系數(shù)為

(4)

(5)

Hilbert變換記為H，在二維情形下，定義正交基及對偶正交基分別為

(6)

(7)

(8)

那么，二維波原子變換系數(shù)為

(9)

該系數(shù)參量是一個關(guān)于j、m、n三參數(shù)的三維矩陣。j為該矩陣包含的系數(shù)矩陣組數(shù)，這里j分別取1和2，得到兩組系數(shù)矩陣，每組系數(shù)矩陣都是二維矩陣，其維度分別為m1×n1、m2×n2。在后續(xù)濾波中，需對這兩組系數(shù)分別進行處理；在最后的重構(gòu)過程中，還需將處理過的兩組系數(shù)進行組合。

對多道地震信號做二維濾波能很好地保持有效信號之間的相關(guān)性，這也是其優(yōu)于一維濾波之處。本文選擇二維波原子變換正是充分考慮地震記錄中各道有效信號之間的相關(guān)性，使有效信號得到更大程度的凸顯。

1.2 自適應(yīng)Wiener濾波

維納濾波器[30-31]是Wiener在二十世紀(jì)四十年代提出的一種基于最小均方誤差準(zhǔn)則的自適應(yīng)最佳線性濾波器。該方法被應(yīng)用于眾多領(lǐng)域，如地震勘探[32-33]、生物醫(yī)學(xué)[34]、遙感圖像[35]以及語音信號處理[36]等。

方法中需處理二維波原子分解系數(shù)，因此需采用二維自適應(yīng)Wiener濾波。假設(shè)一個不含噪信號為ξ(l,k)，加性隨機噪聲為η(l,k)，那么含噪信號為

s(l,k)=ξ(l,k)+η(l,k)

(10)

式中l(wèi)、k分別為二維信號的行、列元素序號。

(11)

式中L、K分別為二維信號的行、列元素個數(shù)。二維濾波時，需選取濾波局域窗(掩模)。那么，在局域窗中采樣點的局部均值和方差可分別表示為

(12)

(13)

式中：r1×r2為局域窗維度；P、Q分別為局部區(qū)域內(nèi)行、列元素個數(shù)，p、q為相應(yīng)元素序號。對于該區(qū)域內(nèi)的采樣點，其二維自適應(yīng)Wiener濾波可表示為

(14)

式中ν2為噪聲方差。

如果噪聲方差未知，則可用所有局部估計方差的均值代替。根據(jù)式(14)求出各局部區(qū)域的信號估計值，再將這些估計值合并為一個整體矩陣即為最終估計值。Wiener濾波的此種實現(xiàn)方法有別于其他傳統(tǒng)實現(xiàn)方法。傳統(tǒng)實現(xiàn)方法中，需求解Wiener濾波器的響應(yīng)函數(shù)，然后用輸入信號卷積(時域?qū)崿F(xiàn))或乘積(頻域?qū)崿F(xiàn))此響應(yīng)函數(shù)而得到輸出信號。以最小均方誤差準(zhǔn)則(式(11))為濾波誤差約束條件。

但是，該誤差公式用到信號的期望值，而該期望值在實際中很難得到。式(14)所示的Wiener濾波實現(xiàn)過程不依賴于期望信號，主要是利用局部均值與方差得到有效信號的估計值，是一種更實用的自適應(yīng)濾波方法。

2 在地震數(shù)據(jù)處理中的應(yīng)用

2.1 處理過程概述

前文述及，對波原子分解系數(shù)進行濾波處理，需去除不必要的系數(shù)，即表征噪聲的系數(shù)，保留表征有效信號的系數(shù)，然后對這些有效信號系數(shù)進行波原子重構(gòu)即可得到濾波后的信號。本文采用自適應(yīng)Wiener濾波對分解系數(shù)進行處理，是利用其在最小均方誤差意義下可達(dá)到最佳線性濾波的特點。

2.2 模擬地震記錄處理

首先針對含有線性同相軸的模擬地震記錄進行實驗。該記錄為2ms采樣，道間距為10m，共有50道。其中包含四個線性同相軸，每個同相軸由主頻為30Hz的Ricker子波構(gòu)成。對該不含噪模擬地震記錄(圖1a)加入高斯白噪聲，使其信噪比約為-5dB(圖1b)。

SVD閾值(方法1)、二維多尺度小波分解閾值(方法2)和二維波原子分解閾值(方法3)三種濾波方法中，閾值的選取須綜合考慮隨機噪聲壓制與有效信號保幅，并歷經(jīng)多次實驗測試。三種方法濾波后重構(gòu)的地震記錄依次如圖1c～圖1e所示。

對含噪記錄進行二維波原子分解系數(shù)自適應(yīng)Wiener濾波(本文方法)時，考慮記錄中有效同相軸分布情況，選取的局域窗為長寬不等的矩形窗，即r1=5、r2=187。當(dāng)噪聲強度改變時，可適當(dāng)調(diào)整r1和r2值。由此，本文方法濾波后重構(gòu)的地震記錄如圖1f所示。

從濾波結(jié)果可見：方法1對背景噪聲的去除及對有效信號的保留效果均較差(圖1c中矩形框)；方法2對有效同相軸的損傷較嚴(yán)重且在同相軸邊緣產(chǎn)生陰影和模糊現(xiàn)象(圖1d中橢圓和矩形框)；方法3濾波效果較前二者有所提高，但模擬記錄兩邊地震道有效信號損失較嚴(yán)重(圖1e中矩形框)；本文方法處理效果最好，在隨機噪聲壓制和有效信號保持方面均表現(xiàn)較優(yōu)(圖1f)。據(jù)圖中右側(cè)能量值色標(biāo)也能看出圖1f中有效信號能量顯著強于圖1c～圖1e。經(jīng)定量計算，得到四種方法濾波后記錄的信噪比分別為0.1686、7.6282、7.8338和8.7075dB。

對上述純凈記錄加入不同強度隨機噪聲，得到不同信噪比含噪記錄，分別采用四種方法做濾波實驗，濾波前后地震記錄信噪比數(shù)據(jù)如表1所示。

從表1可見：方法1所能達(dá)到的信噪比很低，對原含噪記錄信噪比的提高較少；方法2、方法3的濾波后記錄能達(dá)到較高信噪比，其中方法3信噪比提升幅度稍高于方法2，但二者都存在同相軸局部損失較多的問題，且方法2還會導(dǎo)致有效同相軸邊緣不清晰；本文方法相對于其他三種方法優(yōu)勢更明顯，對有效同相軸的增強及邊緣保留效果最好，對含噪記錄信噪比的提升能力最強。

本次實驗采用一臺高性能便攜式工作站。在相同計算機條件下，上述四種濾波方法處理圖1b所示含噪記錄的時間依次為：0.0071、0.7044、3.4747、3.9913s?？梢姺椒?雖用時最短，但效果不佳；方法2用時次短，但效果也欠理想；方法3和本文方法的用時雖較前二者長，但本文方法處理效果更好。

再對含有雙曲同相軸的地震記錄(圖2)做濾波實驗。該數(shù)據(jù)采樣間隔為1ms，道間距為20m，采用主頻為35Hz的Ricker子波，共有76道。層速度分別為2300、2500、2900m/s。對不含噪記錄(圖2a)加入高斯白噪聲使其信噪比為-5dB(圖2b)，三種方處理后的地震記錄分別如圖2c～圖2e所示。

由于圖2a同相軸的形態(tài)和分布與圖1a不同，因此需調(diào)整Wiener濾波窗函數(shù)。調(diào)整時須兼顧隨機噪聲壓制與有效信號保持。此時，選取r1=5、r2=117，即可得到本文方法處理后的地震記錄(圖2f)。

觀察四種濾波方法所得濾波記錄(圖2)可知：方法1對背景噪聲的消除能力很有限且對有效同相軸損失較多(圖2c中矩形和橢圓框)；方法2和方法3能較有效地壓制隨機噪聲，但對有效信號的損失不容忽視，前者會造成同相軸邊緣模糊(圖2d中矩形和橢圓框)，后者會使雙曲同相軸拱起部分損失較嚴(yán)重(圖2e中橢圓框)；而本文方法能在隨機噪聲壓制和有效信號保持方面做到很好權(quán)衡(圖2f)。四種方法處理后記錄的信噪比分別為0.3301、7.1616、7.4222、9.1073dB。

圖2 含有雙曲同相軸的模擬地震記錄及其處理結(jié)果

同樣，對上述純凈記錄加入不同強度隨機噪聲，得到不同信噪比含噪記錄，分別采用四種方法進行濾波實驗，濾波前、后的信噪比數(shù)據(jù)如表2所示。分析對比該表數(shù)據(jù)，也能得到與表1相似的結(jié)論。

表2 含噪的雙曲同相軸記錄濾波前、后信噪比(單位：dB)

2.3 實際地震數(shù)據(jù)處理

選取實際共炮點記錄做濾波對比測試。該共炮點記錄共有168道，采樣間隔為1ms，僅截取記錄的上部(0～2048ms)。首先分別采用方法2、方法3和本文方法進行濾波處理(圖3)。此次處理選取的Wiener濾波窗為正方形，即r1=r2=95。

從濾波結(jié)果(圖3)可見：方法2(圖3b)和方法3(圖3c)濾波效果明顯不理想，對有效同相軸造成較大程度的損失及畸變；而本文方法使有效同相軸更清晰、連續(xù)，能量得到顯著增強(圖3d中梯形和矩形框)。

圖3 實際地震記錄濾波效果對比(一)

對該實際記錄還應(yīng)用方法1和一維時頻峰值濾波法進行處理。方法1處理時，若選取過少特征值，就會導(dǎo)致有效同相軸損失太多；若選取較多特征值，就幾乎沒有去噪效果。對于此次實際記錄，選取分解后奇異值矩陣中的前60個特征值進行重構(gòu)較適宜(圖4b)。但該方法對背景噪聲的壓制效果非常有限，且對有效同相軸造成了一定程度的畸變。

進行時頻峰值濾波時，濾波窗長不宜太大或太小，若太大雖能較好地壓制隨機噪聲，但有效信號損失較大；若太小則不能很好地壓制隨機噪聲。對于該實際記錄，選取濾波窗長為13點(圖4c)。與本文方法所得結(jié)果(圖4d)對比，可見一維時頻峰值濾波處理結(jié)果雖然背景噪聲有所壓制，但同相軸的清晰度和連續(xù)性提升效果欠理想，而本文方法處理后的記錄中同相軸的清晰度更高、連續(xù)性更好。

圖4 實際地震記錄濾波效果對比(二)

3 結(jié)論

本文提出二維波原子分解系數(shù)自適應(yīng)Wiener濾波方法，對地震資料進行隨機噪聲消減處理以增強有效同相軸。其實現(xiàn)過程是，先對地震記錄進行二維波原子分解，再對所得系數(shù)進行二維自適應(yīng)Wiener濾波，最后對濾波后的系數(shù)進行波原子重構(gòu)得到同相軸增強的記錄。通過對人工合成記錄與實際記錄的濾波處理，得到以下認(rèn)識和結(jié)論：

(1)通過對模擬地震記錄的濾波實驗，并與SVD閾值濾波、二維多尺度小波分解閾值濾波及二維波原子分解閾值濾波結(jié)果作比較，驗證本文所提濾波方法在地震隨機噪聲壓制及有效信號保持方面具有優(yōu)越性。

(2)采用二維多尺度小波分解閾值濾波、二維波原子分解閾值濾波、SVD閾值濾波、一維時頻峰值濾波及二維波原子分解系數(shù)自適應(yīng)Wiener濾波五種方法對實際共炮點記錄進行濾波處理，并將處理后的結(jié)果作對比，發(fā)現(xiàn)本文本文方法能夠較徹底地壓制記錄中的隨機噪聲，同時使有效同相軸更清晰、連續(xù)性更好。

因此，應(yīng)用本文方法對地震數(shù)據(jù)進行處理，可為后續(xù)的地震成像及綜合解釋提供可靠資料，對提高地震資料解釋的準(zhǔn)確性具有重要意義。