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

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

0 引言

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

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

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

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

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

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

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

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

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

1 方法原理

1.1 波原子變換原理

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

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

(1)

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

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

(2)

(3)

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

(4)

(5)

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

(6)

(7)

(8)

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

(9)

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

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

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

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

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

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

(10)

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

(11)

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

(12)

(13)

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

(14)

式中ν2為噪聲方差。

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

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

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

2.1 處理過(guò)程概述

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

2.2 模擬地震記錄處理

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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