劉佳賓，王四巍，王艷艷，丁 聰，陳鈞龍

(1.華北水利水電大學(xué)巖土工程與水工結(jié)構(gòu)研究院, 鄭州 450046;2.河南省巖土力學(xué)與結(jié)構(gòu)工程重點(diǎn)實(shí)驗(yàn)室，鄭州 450046；3.河南省交通規(guī)劃設(shè)計(jì)研究院股份有限公司，鄭州 451450)

0 引言

噪聲壓制是地震勘探資料成像處理的重要環(huán)節(jié)，噪聲壓制效果是決定最終成像效果的關(guān)鍵因素之一。最早被人們使用的地震信號(hào)數(shù)字濾波方法是傅里葉變換[1]，該方法通過(guò)將地震信號(hào)由時(shí)間域轉(zhuǎn)換至頻率域，建立合適的濾波器，濾去干擾波的頻譜成分，然后反變換至?xí)r間域來(lái)實(shí)現(xiàn)去噪目標(biāo)。1946年，Gabor.D[2]提出了短時(shí)傅里葉變換方法，目前該方法仍被廣泛使用。但是傅里葉變換在頻率域的分辨率與時(shí)間域的分辨率二者互相制約：時(shí)間分辨率與頻率分辨率二者不能同時(shí)任意提高。傅里葉變換是針對(duì)平穩(wěn)信號(hào)進(jìn)行的時(shí)頻分析，對(duì)于地震信號(hào)這樣的非平穩(wěn)信號(hào)，傅里葉變換只能給出整體信號(hào)中包含的某一頻率分量的平均值，不能完整把握信號(hào)在某一時(shí)刻的本質(zhì)特征，因此嚴(yán)重制約了該方法精度。信號(hào)在時(shí)頻聯(lián)合平面上具有多個(gè)不同的能量聚集區(qū)域，這種信號(hào)稱之為多分量信號(hào)，即時(shí)頻面上的每一個(gè)分量對(duì)應(yīng)一個(gè)信號(hào)分量。1948年，J.Ville[3]把Wigner分布引入非平穩(wěn)信號(hào)時(shí)頻分析領(lǐng)域中，該方法的優(yōu)點(diǎn)在于不需要受到窗函數(shù)的制約，但是在處理多分量信號(hào)時(shí)，會(huì)出現(xiàn)交叉項(xiàng)，并不適用于復(fù)雜地震信號(hào)[4]。1982年，Morlet J.[5]提出了一種新的時(shí)頻分析方法即小波變換。該方法通過(guò)選擇尺度因子與平移因子，根據(jù)信號(hào)的特點(diǎn)對(duì)時(shí)窗進(jìn)行調(diào)整，可以在信號(hào)的不同位置提供不同的分辨率。該方法在濾波過(guò)程中可以對(duì)信號(hào)的局部特征進(jìn)行更好的識(shí)別，更準(zhǔn)確地區(qū)分有效信號(hào)和噪聲，取得更好的濾波效果。段春節(jié)等人[6]將小波變換方法與相干體技術(shù)相結(jié)合，提出了一種基于小波變換的相干體計(jì)算的新方法，得到突出特定頻帶的相干體。吳雅娟等人[7]通過(guò)改進(jìn)的小波閾值法在有效濾除測(cè)井信號(hào)中的噪聲的同時(shí)，又可以保留曲線的細(xì)節(jié)信息。

小波變換方法自20世紀(jì)90年代，C.S.Luchele等人[8]利用小波分析刻畫地質(zhì)體的不均勻性；高靜懷[9]采用解析小波的方法求取地震信號(hào)的瞬時(shí)參數(shù)，深入地分析了不同尺度下地震記錄小波變換結(jié)果及其對(duì)應(yīng)的瞬時(shí)參數(shù)含義。王姣等人[10]針對(duì)互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解(CEEMD)舍棄高頻分量的去噪方法和小波閾值去噪方法存在的不足，利用小波閾值去噪方法對(duì)地震信號(hào)進(jìn)行濾波。王江濤[11]基于小波多尺度的分析技術(shù)，提出了基于小波分析的動(dòng)態(tài)變形信息提取多尺度分析方法,并應(yīng)用于強(qiáng)噪聲背景下提取弱信號(hào)及信號(hào)的奇異性檢驗(yàn)，有效提取了動(dòng)態(tài)變形信息。阮清青[12]通過(guò)對(duì)修正S變換與常規(guī)時(shí)頻分析方法的對(duì)比，根據(jù)合成信號(hào)以及實(shí)測(cè)地震記錄的時(shí)頻分析得知修正S變換較常規(guī)的時(shí)頻分析方法具有更好的時(shí)頻分辨率和能量聚集性。

顯而易見，小波變換已廣泛應(yīng)用于地震數(shù)據(jù)壓縮、地震數(shù)據(jù)特征的分析、以及多尺度反演等領(lǐng)域，且其處理效果主要依賴于小波基函數(shù)的選擇。本文通過(guò)實(shí)例對(duì)地震信號(hào)不同小波基函數(shù)的濾波效果進(jìn)行比對(duì)，討論小波基函數(shù)的選取，并采用二維小波包濾波分析方法對(duì)地震剖面進(jìn)行處理，驗(yàn)證濾波效果。

1 小波變換濾波方法

1.1 小波變換原理

小波變換是采用待處理信號(hào)的一組基本小波函數(shù)在空間上的投影對(duì)信號(hào)進(jìn)行描述。由不同尺度因子(a)與時(shí)移因子(b)構(gòu)成的小波基表達(dá)式如下：

(1)

式中Ψ(t)表示某類小波基函數(shù)，上式須滿足a≠0。

使用連續(xù)小波變換(CT)把信號(hào)x(t)映射到時(shí)間—尺度平面：

(2)

式中Ψ*(t)表示某類小波基函數(shù)的共軛表達(dá)式。

(2)式的逆變換為：

(3)

式中，CΨ是小波的容許性條件，即當(dāng)CΨ滿足下式(4)時(shí)，可以由小波變換的結(jié)果重構(gòu)原始信號(hào)。

(4)

其中Ψ(w)為Ψ(t)的傅里葉變換。

通過(guò)改變小波基的尺度因子，可以改變其對(duì)所處理信號(hào)的分辨率。當(dāng)遇到低頻信號(hào)時(shí)，増大小波函數(shù)在時(shí)間上的尺度，降低被處理信號(hào)的時(shí)間分辨率，從而達(dá)到提高其頻率分辨率的目的；遇到高頻信號(hào)時(shí)，縮小小波函數(shù)的時(shí)間尺度，使之具有更高的時(shí)間分辨率，更好的表述信號(hào)的高頻部分。本文采用二尺度分解的方法，首先將待處理的地震信號(hào)S分解為高頻(d1)與低頻(a1)兩部分,然后對(duì)低頻(a1)部分再次進(jìn)行分解，得到高頻(d2)與低頻(a2)分量。此時(shí)原始信號(hào)S等于a2、d2與d1之和，對(duì)轉(zhuǎn)換至小波域的地震信號(hào)進(jìn)行閾值處理即可達(dá)到濾波的效果。小波域中較小的小波系數(shù)主要包含的是隨機(jī)噪聲，所以對(duì)較小的小波系數(shù)進(jìn)行壓制后，即可重構(gòu)得到壓制噪聲之后的信號(hào)。

1.2 小波變換濾波原理

閾值去噪簡(jiǎn)而言之就是對(duì)信號(hào)進(jìn)行分解，然后對(duì)分解后的系數(shù)進(jìn)行閾值處理，最后重構(gòu)得到去噪信號(hào)。對(duì)原始信號(hào)進(jìn)行小波變換后，得到各尺度系數(shù)。隨后確定閾值λ，若小波系數(shù)小于λ，則認(rèn)為該系數(shù)主要由噪聲引起，需要對(duì)該部分系數(shù)進(jìn)行去除。若小波系數(shù)大于λ，則認(rèn)為該系數(shù)主要是由信號(hào)引起，對(duì)該部分系數(shù)加以保留。對(duì)處理后的小波系數(shù)進(jìn)行小波反變換，即可得到濾波之后的信號(hào)。

常見的閾值函數(shù)有：硬閾值函數(shù)(小波系數(shù)的絕對(duì)值低于閾值的置零，高于的保留不變)；軟閾值函數(shù)(小波系數(shù)的絕對(duì)值低于閾值的置零，高于的系數(shù)shrinkage處理)。本文采用軟閾值函數(shù)進(jìn)行去噪處理。軟閾值函數(shù)公式為：

(5)

式中：w表示小波系數(shù)，T為給定閾值，sign(*)符號(hào)函數(shù)。

1.3 小波變換濾波效果評(píng)價(jià)

本文采用求取均方差與相似度的方法對(duì)地震記錄不同小波基的濾波效果進(jìn)行評(píng)價(jià)。

均方差σ是總體各單位標(biāo)準(zhǔn)值與其平均數(shù)離差平方的算術(shù)平均數(shù)的平方根。它反映組內(nèi)個(gè)體間的離散程度。定義式如下：

(6)

由定義可知，濾波后信號(hào)與原始地震記錄之間的均方差越小，二者相似程度越高，濾波結(jié)果越接近原始地震記錄。相似度ρ的定義式如下：

(7)

2 一維信號(hào)的濾波結(jié)果分析

以鄂爾多斯盆地北部地震資料為例，該區(qū)主要目的層氣藏多以巖性控制為主，少部分為構(gòu)造—巖性型油氣藏。原始地震剖面如圖1所示，在主要目的層(1 700～1 900ms)山西組、太原組、石盒子組，主要發(fā)育有T9b+c、T9d、T9e、T9f四組反射波，形成一個(gè)強(qiáng)振幅、連續(xù)性好的密集反射段。測(cè)線西部局部地段反射雜亂、連續(xù)性差，出現(xiàn)反射同相軸的合并、分叉等現(xiàn)象。提取該地震剖面中的第308道地震記錄，對(duì)其進(jìn)行小波分析。小波變換的尺度越大，越容易獲取信號(hào)的低頻特征；反之，尺度越小，獲取的主要是信號(hào)的高頻特性。地震信號(hào)的小波變換是一種先驗(yàn)的時(shí)頻分析方法，其小波分析結(jié)果不僅取決于信號(hào)本身，而且也與所采用的分析小波有關(guān)。在進(jìn)行地震信號(hào)的薄互層高分辨率分析時(shí)，恰當(dāng)?shù)剡x擇小波函數(shù)是至關(guān)重要的。

在常用的小波基函數(shù)中，考慮正交性、雙正交性、緊支撐等性質(zhì)，選取Daubechies[13]、Coiflets、Symlets三種函數(shù)(縮寫分別為db、coif、sym，表示形式為：dbN、coifN、symN；N為消失矩)。

圖1 待處理地震剖面Figure 1 Seismic section to be processed

首先選用db3小波，該小波系是一系列二進(jìn)制小波的總稱，是由法國(guó)學(xué)者Daubechies提出的，故稱為db小波。Daubechies小波不僅是連續(xù)和正交的，而且是支集最小的，這種小波的濾波器個(gè)數(shù)少，在地震信號(hào)的分解和重構(gòu)中計(jì)算量小。該小波沒有明確的解析表達(dá)式，其有效支撐長(zhǎng)度為2N，消失矩為N。實(shí)際地震信號(hào)中的噪音通常表現(xiàn)為高頻成分，因此對(duì)信號(hào)進(jìn)行二尺度一維分解，所得結(jié)果如圖2所示。隨后采用軟閾值處理方法進(jìn)行小波閾值濾波，其結(jié)果如圖3所示。

然后采用sym2小波對(duì)信號(hào)進(jìn)行分析。Symlets小波與dbN小波相比，在連續(xù)性、支集長(zhǎng)度、濾波器長(zhǎng)度等方面與dbN小波一致，但symN小波具有更好的對(duì)稱性，即一定程度上能夠減少對(duì)信號(hào)進(jìn)行分析和重構(gòu)時(shí)的相位失真。Symlets小波也是對(duì)一系列小波的總稱，這類正交小波的支撐長(zhǎng)度為2N-1，濾波器的長(zhǎng)度為2N，消失矩為N，具有相似的對(duì)稱性。此次采用尺度函數(shù)與db3相近的sym2小波對(duì)信號(hào)進(jìn)行二尺度一維分解，進(jìn)行軟閾值濾波后，得到的結(jié)果如圖4所示。

s：原始信號(hào)；a2：二階低頻分量；d1：一階高頻分量；d2：二階高頻分量圖2 db3小波變換結(jié)果Figure 2 Transformed results of db3 wavelet

圖3 db3小波閾值去噪結(jié)果Figure 3 Denoised results of db3 wavelet threshold

圖4 sym2小波閾值去噪結(jié)果Figure 4 Denoised results of sym2 wavelet threshold

隨后采用coif1小波對(duì)原始信號(hào)進(jìn)行變換。Coiflet的小波函數(shù)的2N階矩為零，具有比dbN更好的對(duì)稱性。Coiflet小波系屬于雙正交的小波，該類小波的支撐長(zhǎng)度為6N-1，濾波器的長(zhǎng)度為6N，消失矩為2N，其對(duì)稱性要好于db小波。選用尺度函數(shù)與db3小波接近的coif1小波對(duì)原始信號(hào)進(jìn)行二尺度一維分解，進(jìn)行軟閾值濾波后所得結(jié)果如圖5所示。

圖5 coif1小波閾值去噪結(jié)果Figure 5 Denoised results of coif1 wavelet threshold

通過(guò)選用三種不同的小波基函數(shù)對(duì)第308道數(shù)據(jù)進(jìn)行了閾值濾波后，采用求取濾波后信號(hào)與原始地震信號(hào)的均方差與相似度的方法進(jìn)行比較，選取最優(yōu)化的濾波結(jié)果。詳情見表1。

表1 濾波結(jié)果的均方差與相似度

由上表可知，從均方差來(lái)看，sym2去噪結(jié)果均方差最小，說(shuō)明其對(duì)噪聲消除效果最好，db3去噪結(jié)果次之，coif1去噪結(jié)果均方差最大；從相似度來(lái)看，db3去噪結(jié)果與原始信號(hào)的相似程度最高，coif1去噪結(jié)果與原始信號(hào)相似程度次之，而sym2去噪結(jié)果與原始信號(hào)相似程度最差。綜合考慮來(lái)看，對(duì)該工區(qū)地震記錄進(jìn)行小波閾值濾波時(shí)選用db3小波基函數(shù)最為合適。

3 小波包對(duì)二維信號(hào)的濾波結(jié)果分析

經(jīng)過(guò)上一節(jié)的討論可以得出結(jié)論，選用db3小波基函數(shù)對(duì)該工區(qū)地震記錄進(jìn)行濾波處理效果最好，因此采用db3小波對(duì)二維地震剖面進(jìn)行小波包濾波。

小波包分析是小波分析的延伸，該方法對(duì)信號(hào)的時(shí)頻分析更為細(xì)致。其原理如下：小波分解實(shí)際是將信號(hào)分成高頻與低頻兩個(gè)部分，篩選出高頻部分之后對(duì)剩余的低頻部分再一次做分解，再次篩選出高頻與低頻部分，直到達(dá)到規(guī)定的尺度。小波包分析與之不同之處在于，針對(duì)每一次篩選出來(lái)的高頻部分繼續(xù)做分解，也就是說(shuō)對(duì)每一次分解得到的所有結(jié)果繼續(xù)進(jìn)行分解。由于小波變換不具備自適應(yīng)的特征，所以高頻成分中可能也包含有低頻信息，因此采用小波包分析能夠?qū)⒏哳l成分中的低頻信息也分解出來(lái)，其時(shí)頻分析能力會(huì)得到進(jìn)一步的提升。

(a)原始信號(hào)

(b)濾波信號(hào)圖6 二維地震數(shù)據(jù)的濾波效果對(duì)比Figure 6 Filtering effects comparison of 2D seismic data (left: original; right: denoised)

圖7 A、B區(qū)域?qū)Ρ确糯髨DFigure 7 Comparison of A and B regions enlarged images

對(duì)二維地震數(shù)據(jù)進(jìn)行濾波處理后，選取地震剖面中的幾處細(xì)節(jié)進(jìn)行放大對(duì)比，對(duì)濾波效果進(jìn)行驗(yàn)證，如圖6、7所示。從圖中可以看出，原始信號(hào)剖面圖中存在許多“毛刺”現(xiàn)象，而進(jìn)行小波閾值濾波之后，濾波后的圖中高頻干擾明顯降低，同向軸的形態(tài)更為清晰圓滑。圖7為濾波效果局部放大圖，效果更佳明顯?？梢姴捎胐b3小波進(jìn)行濾波之后，所得剖面中的隨機(jī)噪音得到了很好的壓制，有效信號(hào)更加突出，信噪比得到了提高。

4 結(jié)論

本文通過(guò)對(duì)實(shí)際工區(qū)資料進(jìn)行小波變換濾波處理，并對(duì)結(jié)果進(jìn)行分析總結(jié)后，可以得出一下幾點(diǎn)結(jié)論。

①不同的小波基函數(shù)選取會(huì)對(duì)濾波效果產(chǎn)生影響，需要結(jié)合質(zhì)特點(diǎn)和地震響應(yīng)特征，多嘗試幾種不同的小波基函數(shù)，并從中找到效果最佳的小波基函數(shù)進(jìn)行分析。

②針對(duì)該工區(qū)地震資料，db3小波閾值濾波效果最好。處理后的反射層特征變得更清晰，同相軸連續(xù)性增強(qiáng)，資料頻帶向高頻方向拓寬，高頻成分明顯加強(qiáng)，低頻部分有所保持；隨機(jī)噪音得到了很好的壓制，有效信號(hào)更加突出，信噪比得到了提高。

③小波包濾波能夠有效降低地震剖面中的高頻干擾，其結(jié)果可用于精細(xì)解譯地震波信號(hào)，更好地服務(wù)于氣田地震信號(hào)的解譯工作。