陳偉，楊浪，陳名德，連辰浩，江源

1.油氣資源與勘探技術教育部重點實驗室(長江大學)，湖北 武漢 430100 2.長江大學地球物理與石油資源學院，湖北 武漢 430100

隨著全球原油需求日益增長，油田勘探程度不斷提升，促使勘探重點轉向隱蔽油氣藏，該類油氣藏受埋藏深度影響，多由致密砂巖組成，具有超高壓、超深、低孔、低滲等特點，這對勘探技術帶來新的挑戰(zhàn)，而瞬時譜分析技術給油氣勘探向更深處發(fā)展帶來希望。地震信號是非線性的、非平穩(wěn)的信號，它們的能量分布可以通過時頻分析獲得，因此可以在時域和頻域，即聯(lián)合時頻分布中完全建立信號能量或強度的分布。時頻分析方法可以觀察時域和頻域信號的演化，并顯示信號的局部時頻特征。因此，時頻分析方法已被廣泛應用于地震波分析中。瞬時譜分析技術應用到地球物理領域中檢測油氣藏的存在現(xiàn)已成為熱點。瞬時譜分析技術主要依賴于時頻分析方法的優(yōu)劣，如今廣泛應用的時頻分析方法主要分為兩類，一類為線性時頻分析方法，如Gabor變換，小波變換，S變換等；另一類為非線性時頻分析方法，如Hilbert-Huang變換等。

Fourier變換只能將時間域的信號整體變換到頻率域，將信號的頻率成分和分布表示出來，雖然使信號擁有了頻率分辨率，卻丟失了時間分辨率，不能反映出信號的頻率特征隨時間的變化情況。為了解決該問題，GABOR[1]于1946年首次提出在時頻域內(nèi)分解一維信號的方法，即短時Fourier變換。短時Fourier變換中最具代表的為Gabor變換。Gabor變換選取高斯函數(shù)作為窗函數(shù)，在信號的每一小段時間間隔[2]內(nèi)進行Fourier變換，該窗口在時間軸上進行平移，對每一段信號進行局部Fourier分析后，再將變換得到的頻率信息放到時間-頻率的二維平面上，從而得到時頻譜。Gabor變換選取高斯函數(shù)作為窗函數(shù)的優(yōu)點在于高斯函數(shù)的Fourier變換還是高斯函數(shù)，不會對變換結果產(chǎn)生不良影響，可以同時提供時域和頻域局部化的信息。另外，根據(jù)Heisenberg測不準原理，變換得到的結果不能同時擁有極高的時間分辨率和頻率分辨率，極大地影響了分析結果的好壞。而高斯函數(shù)窗口面積已達到Heisenberg測不準原理的下界，能兼顧時間軸和頻率軸的分辨率。Gabor變換的基函數(shù)不能成為正交基[3]，為保證不丟失信息，只能采用非正交基參與信號分析，使得計算復雜度更高。Gabor變換的窗口固定不變，無法針對高低頻信號給出可以調(diào)節(jié)的時頻窗，割斷了頻率與窗口長度的內(nèi)在聯(lián)系[4]；對短時突發(fā)信號進行分析時，可能產(chǎn)生時間或頻率上的模糊現(xiàn)象[5]。MORLET[6]等人從短時Fourier變換中推出小波變換，能對信號同時在時間域和頻率域內(nèi)進行局部化分析。與Fourier變換不同，小波變換有兩個變量：尺度和平移量。尺度控制小波函數(shù)的伸縮，平移量控制小波函數(shù)的平移。尺度就對應于頻率(反比)，平移量就對應于時間。尺度增大，則時窗伸展，頻寬收縮，帶寬變窄，中心頻率降低，時間分辨率降低而頻率分辨率增高；尺度減小，則相反[7]。這符合實際問題中高頻信號持續(xù)時間短、低頻信號持續(xù)時間長的自然規(guī)律。小波變換能準確反映待測信號的幅頻特性。在小波變換中，時頻窗口的長度和寬度隨著參數(shù)尺度的變化而變化[8]。雖然小波變換克服了短時Fourier變換單一分辨率分析的不足，但是由于引入的尺度因子與頻率沒有直接的聯(lián)系，且在小波變換中沒有明顯表現(xiàn)出來，因此小波變換的結果不是一種真正的時頻譜。小波分析的另一個問題是其具有的自適應特點，一旦基本小波被選定，就必須用它來分析所有待分析的數(shù)據(jù)。STOCKWELL等[9]結合短時Fourier變換和小波變換提出S變換，S變換相對于小波變換的不同在于其公式中加入了相位項，用于進行相位校正。它與小波變換一樣存在基本小波，其基本小波由簡諧波函數(shù)與高斯函數(shù)相乘組成，其中，簡諧波函數(shù)在時間域內(nèi)作伸縮變換，而高斯函數(shù)則作伸縮和平移變換[10]。S變換既克服了短時Fourier變換和小波變換的不足，也繼承了它們的優(yōu)點[11-12]。但是S變換的窗函數(shù)隨頻率以固定的趨勢改變，基本小波一經(jīng)選定則無法更改。PINNEGAR等[13]、高靜懷等[14]在S變換的基礎上提出廣義S變換，通過各種可調(diào)的屬性參數(shù)控制窗函數(shù)，使其具有更高的適應性。

上述方法都是在Fourier變換的基礎上發(fā)展起來的，但是Fourier變換存在時間域和頻率域的相互制約，無法兼顧時間軸和頻率軸的分辨率，并且在處理非線性非平穩(wěn)信號時效果較差。1998年，HUANG等[15]首次提出希爾伯特-黃變換，這是一種用于處理非線性、非平穩(wěn)信號的時頻分析新方法，該方法由經(jīng)驗模態(tài)分解[16]和Hilbert譜分析兩部分組成。經(jīng)驗模態(tài)分解就是將非平穩(wěn)信號按頻率由高到低分解為一系列固有模態(tài)函數(shù)，這一系列固有模態(tài)分量能反映出原始信號在不同頻率上的局部特征。對每個固有模態(tài)分量進行Hilbert變換，得到每一個固有模態(tài)分量隨時間變換的瞬時頻率和瞬時幅值，求得振幅-頻率-時間的三維譜分布[17]。把所有固有模態(tài)分量的Hilbert譜進行綜合之后即可求出原始信號的Hilbert譜，以此來得到信號的時頻屬性。希爾伯特-黃變換很適合處理非線性非平穩(wěn)信號，在地震資料處理和瞬時譜分析中有著廣泛應用。然而經(jīng)驗模態(tài)分解存在模態(tài)混疊現(xiàn)象。模態(tài)混疊主要是指同一個固有模態(tài)分量當中出現(xiàn)了不同尺度或頻率的信號，或者同一尺度或頻率的信號被分解到多個不同的固有模態(tài)分量當中[18]，嚴重影響了數(shù)據(jù)分析的結果。為了克服該問題，學者提出了總體經(jīng)驗模態(tài)分解算法[19]。總體經(jīng)驗模態(tài)分解方法是一種噪音分析法，即在每次信號分解的過程中對原始信號添加高斯白噪聲，該方法能夠很好地解決經(jīng)驗模態(tài)分解算法在分解過程中出現(xiàn)的端點效應[20]和模態(tài)混疊現(xiàn)象?？傮w經(jīng)驗模態(tài)分解方法在原始信號中加入了高斯白噪聲，所以分解后的結果很難重構得到原始信號，分解不具備很好的完備性。其分解依賴添加白噪聲幅值和集成次數(shù)，如果參數(shù)選擇不合適，不僅不能抑制模態(tài)混淆，而且會出現(xiàn)偽分量；且也無法保證分解得到的分量滿足固有模態(tài)分量的定義條件[21]。總體經(jīng)驗模態(tài)分解方法需要的計算量非常大，計算效率較低，不利于實際地震資料的處理。TORRES等[22]在2010年提出完備總體經(jīng)驗模態(tài)分解，與總體經(jīng)驗模態(tài)分解不同的是，該方法向原始信號中加入了正負成對形式的白噪音，這樣重構信號中的殘余輔助信號被有效消除，從而簡化計算量。完備總體經(jīng)驗模態(tài)分解既繼承了經(jīng)驗模態(tài)分解的正交性、自適應性以及能夠處理非平穩(wěn)信號的特點，又在總體經(jīng)驗模態(tài)分解的基礎上簡化了計算量，幾乎可以重構完整的原始信號。

經(jīng)驗小波變換作為一種全新的自適應性算法在信號處理中得到應用。經(jīng)驗小波變換能夠更好地分解出原始信號中固有的本征信號，因此其較經(jīng)驗模態(tài)分解和總體經(jīng)驗模態(tài)分解算法有更高的自適應性。經(jīng)驗小波變換以成熟的小波理論為基礎，其數(shù)學理論基礎非常充分，自身具有較高的計算效率?？紤]到經(jīng)驗小波變換是基于經(jīng)驗模態(tài)分解的最新自適應算法，本文先介紹了經(jīng)驗模態(tài)分解和經(jīng)驗小波變換的基本原理，然后提出基于經(jīng)驗小波變換的瞬時譜分析技術，最后將其應用到致密砂巖氣藏的檢測中。本文首次將經(jīng)驗小波變換應用到瞬時譜分析技術中檢測致密砂巖氣藏，得到了較好的效果。

1 基本原理

1.1 經(jīng)驗模態(tài)分解

HUANG等[15]人認為，只有固有模態(tài)分量的瞬時頻率具有物理意義。經(jīng)驗模態(tài)分解基于簡單的假設，即任何數(shù)據(jù)都包含不同的簡單固有模式。這些振蕩模式中的每一種都由固有模態(tài)分量表示，具有以下定義：①在整個信號段中，極值的數(shù)量和過零的數(shù)量必須相等或相差最多一個；②在任何時候，由局部最大值定義的包絡的平均值和由局部最小值定義的包絡為零。固有模態(tài)函數(shù)定義的第一個條件類似于靜態(tài)高斯過程的窄帶要求。第二個條件使用局部要求而非全局要求，這使得由對稱波形引起的不需要的波動不會出現(xiàn)在瞬時頻率中。

固有模態(tài)函數(shù)的定義表明，數(shù)據(jù)的內(nèi)部振動模式由固有模態(tài)分量表征。每個固有模態(tài)分量只涉及一種振蕩模式。振蕩也將相對于“局部均值”對稱。線性或非線性的固有模態(tài)分量可以在簡諧分量中具有恒定的幅度和頻率，也可以在時間函數(shù)中具有可變幅度和頻率。經(jīng)驗模態(tài)分解通過將多組分信號簡化為單組分函數(shù)的集合[23]，解決了從多組分信號中計算有意義瞬時頻率的難題。最終的復雜信號是由固有模態(tài)分量的重疊形成的。經(jīng)驗模態(tài)分解的目的是獲得固有模態(tài)分量。通過經(jīng)驗模態(tài)分解方法分解的固有模態(tài)分量，可以通過篩選過程根據(jù)信號本身的相鄰極值點之間的延遲來定義和區(qū)分。

HUANG等[24]表明，為了獲得關于經(jīng)驗模態(tài)分解詳細而經(jīng)驗性的統(tǒng)計知識，經(jīng)驗模態(tài)分解基本上充當了與小波分解相關的二元濾波器組[25]，經(jīng)驗模態(tài)分解方法將信號分解為固有模態(tài)分量的步驟如下：

1)識別時間序列信號X(t)的所有極值點，使用X(t)的所有極大點確定上包絡u(t)，并用所有的極小點確定下包絡v(t)，X(t)滿足：

v(t)≤X(t)≤u(t)

(1)

然后，上包絡線和下包絡線的平均曲線m(t)是：

(2)

設h1(t)=X(t)-m(t)，那么h1(t)就是固有模態(tài)分量。

2)由于過沖和俯沖，包絡樣條將近似產(chǎn)生一個新的極值，它會影響原始極值點的位置和大小。因此，h1(t)不完全滿足固有模態(tài)分量條件。為了獲得所需的h1(t)，讓h1(t)代替X(t)。對應于h1(t)，上包絡線為u1(t)，下包絡線為v1(t)，然后重復這個過程：

(3)

h2(t)=h1(t)-m1(t)

(4)

?

(5)

hk(t)=hk-1-mk-1(t)

(6)

重復該過程直到每個hk(t)滿足固有模態(tài)分量條件。然后，獲得第一個固有模態(tài)分量C1(t)和信號r1(t)的剩余部分：

C1(t)=hk(t)

(7)

r1(t)=X(t)-C1(t)

(8)

在信號的剩余部分繼續(xù)使用經(jīng)驗模態(tài)分解。分解繼續(xù)進行，直到信號的剩余部分單調(diào)或其值小于預定值。分解后得到的固有模態(tài)分量和余量為：

r2(t)=r1(t)-C2(t)

?

rn(t)=rn-1(t)-Cn(t)

(9)

X(t)可以表示為固有模態(tài)分量和余量之和：

X(t)=C1(t)+C2(t)+…+Cn(t)+rn(t)

(10)

1.2 經(jīng)驗小波變換

經(jīng)驗小波變換的目標是通過構建自適應小波來提取不同的單分量。該方法的步驟如下：

1)將快速Fourier變換應用于信號f(t)，從而獲得頻譜X(w)，其中f(t)是離散信號，t={ti}i=1，2，…，M(M表示樣本數(shù))。在Fourier頻譜中找出最大值M={Mi}1，2，…，N(N表示最大值的數(shù)量)，并推導它們的相應頻率W={Wi}i=1，2，…，N。

2)獲得Fourier頻譜和邊界集的適當分段。將每個段的邊界Ωi定義為兩個連續(xù)最大值的中心：

(11)

式中：Wi和Wi+1為兩個頻率；邊界集是Ω={Ωi}i=1，2，…，N-1。

3)定義一組由一個低通濾波器和基于邊界的N-1個帶通濾波器組成N個小波濾波器。尺度函數(shù)φ1(W)和經(jīng)驗小波ψi(W)的Fourier變換的表達式如下：

(12)

(13)

(14)

求取尺度函數(shù)和小波函數(shù)以提取不同單分量。近似系數(shù)由分析信號f(t)的內(nèi)積與經(jīng)驗尺度函數(shù)表示：

(15)

類似地，細節(jié)系數(shù)是通過經(jīng)驗小波的分析信號f(t)的內(nèi)積得到的：

(16)

式中：Wf(i，t)表示第t個時間點的第i個濾波器組的細節(jié)系數(shù)。

怎樣將Fourier譜進行分段在經(jīng)驗小波變換中至關重要，其直接關系到對原始信號進行分解后的自適應程度。經(jīng)驗小波變換將原始信號進行不同的分割，比如對某個中心頻率的緊支撐部分進行分割。假設斷點數(shù)目為N，意味著需要N+1個邊界。除了起點0和終點σ以外，還需要N-1個邊界。為了找到這些邊界，首先對信號頻譜的局部極大值點進行降序排列(起點0和終點σ包括在內(nèi))。假設找到了M個極大值點，下面兩種情況將會出現(xiàn)：①M≥N：算法發(fā)現(xiàn)了足夠的極值點以便于分割原始信號，但只取前N-1個極大值點；②M≤N：信號沒有預期的那么多模態(tài)，將這M個極值點保留，并添加一些近似值直到極值點達到N個。

(17)

(18)

這樣原始信號可以通過下面的式子來重構：

(19)

為簡單起見，經(jīng)驗小波變換所蘊含的經(jīng)驗模態(tài)函數(shù)可以定義為：

(20)

(21)

1.3 基于經(jīng)驗小波變換的瞬時譜分析技術

瞬時譜分析(instantaneous spectral analysis)能夠對地震信號進行連續(xù)時頻分析，觀察任一時間點處振幅變化情況及頻率分布。

基于經(jīng)驗小波變換的瞬時譜分析技術實現(xiàn)步驟為：①利用經(jīng)驗小波變換將某一單道地震信號分解成按頻率高低分布的固有模態(tài)分量；②求取各個固有模態(tài)分量中的瞬時頻率分布；③將各個固有模態(tài)分量的頻率分布依據(jù)時間疊加得到該地震道的時頻剖面；④依次選取二維地震記錄中每個單道記錄，重復進行步驟①～③；⑤從各時頻剖面中提取特定頻率分量，根據(jù)對應的地震道位置依次疊加，最終得到特定頻率的瞬時譜。

評判瞬時譜分析剖面時遵循以下原則：①時頻分布沿著頻率的振幅疊加值與信號的瞬時振幅值近似相等；②時頻分布沿著時間的振幅疊加值與信號的瞬時頻率值近似相等；③地震剖面上明顯的構造，在瞬時譜剖面上也必須表現(xiàn)出來；④地震剖面上明顯的構造異常旁瓣，在瞬時譜剖面上不能單獨出現(xiàn)；⑤一個單獨的同相軸在瞬時譜剖面上是連續(xù)的。

2 合成數(shù)據(jù)

為驗證本文所提方法的有效性，現(xiàn)構造一個典型的非平穩(wěn)信號，分別利用經(jīng)驗小波變換和其他時頻分析方法對其進行對比，如短時Fourier變換，S變換，經(jīng)驗模態(tài)分解。

該信號具體形態(tài)(見圖1)為20 Hz背景余弦波，在300 ms處疊加100 Hz Morlet子波，在1 070 ms和1 100 ms處疊加了兩個30 Hz Richer子波。在1 300～1 700 ms之間，信號具有3種不同的頻率組分：7，30和40 Hz，其中7 Hz的頻率是不連續(xù)的，并且包含少于一個周期的部分，出現(xiàn)在1.37，1.51和1.65 s。該信號具有顯著特征，為典型的非平穩(wěn)信號，適合用于檢驗時頻分析方法的優(yōu)劣。

圖1 測試信號Fig.1 Test signal

2.1 信號分解對比

經(jīng)驗模態(tài)分解將原始信號分解成了7個固有模態(tài)分量(見圖2(a))，在300 ms處信號頻率大幅突變，導致固有模態(tài)分量1(IMF1)中雖提取出了Morlet子波的高頻成分，但混雜了20 Hz背景余弦波、1 100 ms處的Ricker子波等低頻成分。同樣，固有模態(tài)分量2(IMF2)、固有模態(tài)分量3(IMF3)、固有模態(tài)分量4(IMF4)中均是各種組分的高頻部分和低頻部分混雜在同一固有模態(tài)分量中，難以識別出信號包含的具體頻率組分，使信號分析復雜化，這屬于典型的模態(tài)混疊現(xiàn)象。圖2(b)是經(jīng)驗小波變換分解原始信號的結果，IMF1表示1 300～1 700 ms之間的低頻組分(7 Hz)。IMF2中包含20 Hz背景余弦波、1 300～1 700 ms處低頻組分及1 070 ms和1 100 ms處兩個Ricker子波。而IMF4僅提取了300 ms處的Morlet子波，沒有其他組分混雜。IMF3則為少量Morlet子波剩余組分和Ricker子波剩余組分。經(jīng)驗小波變換分解效果明顯優(yōu)于經(jīng)驗模態(tài)分解分解效果。

注：IMF為固有模態(tài)分量。圖2 不同方法的信號分解結果對比Fig.2 Comparison of signal decomposition results of different methods

2.2 時頻譜對比

圖3(a)是經(jīng)過128ms時間窗的短時Fourier變換得到的時頻譜，30 Hz背景余弦波、300 ms處的100 Hz Morlet子波、1 300～1 700 ms處的低頻部分的輪廓已凸顯出來，但時頻分辨率過低，能量顯示較為分散，只能觀察到各頻率組分的大致分布情況，受Heisenberg測不準原理制約，缺乏實際研究價值。

圖3 不同方法處理原始信號得到的時頻譜對比Fig.3 Comparison of time spectrum obtained by processing raw signals with different methods

圖3(b)是通過S變換得到的時頻譜，在300 ms處的頻率分辨率很低，無法識別出100 Hz Morlet子波，1 300～1 700 ms處的時頻分辨率較短時Fourier變換相比有明顯提升，能量顯示更為集中。

圖3(c)顯示通過經(jīng)驗模態(tài)分解得到的瞬時譜，圖中30 Hz背景余弦波，1 300～1 700 ms處低頻組分以及1 100 ms附近的兩個Ricker子波均清晰可見，時頻分辨率較前兩種方法有質的飛躍。在信號兩端可見明顯的發(fā)散現(xiàn)象，屬于典型的端點效應，導致200～1 800 ms間有大量低頻成分混雜，特別是在300 ms處100 Hz Morlet子波部分，其上下均有能量泄露，信號失真十分嚴重。

經(jīng)驗小波變換處理后的瞬時譜(見圖3(d))首次將1 300～1 700 ms處的7 Hz頻率特征準確地刻畫出來。30 Hz背景余弦波、Ricker子波以及1 300～1 700 ms處其他頻率組分均以最高時頻分辨率呈現(xiàn)出來，且300 ms處高頻突變部分標識準確，能量非常集中，信號兩端未出現(xiàn)端點效應。對比上述4種時頻分析方法處理結果，可以發(fā)現(xiàn)經(jīng)驗小波變換明顯優(yōu)于其他3種方法。

3 實際地震資料處理

為了驗證本文方法在含氣檢測中的有效性，針對中國某地區(qū)的實際資料進行瞬時譜分析。該區(qū)儲集層物性總體較差，屬致密-超致密砂巖區(qū)塊。該區(qū)天然氣的富集、產(chǎn)出和高產(chǎn)與裂縫的發(fā)育程度關系極為密切?？紤]到裂縫的存在以及天然氣的影響使地震波速度明顯降低，造成物性參數(shù)轉換上具有不同程度的突變特征，因此該地區(qū)的地震資料具有非線性非平穩(wěn)的特性，適合利用本文的經(jīng)驗小波變換對該區(qū)儲集層的含氣特征進行研究。

取實際地震資料的部分疊后剖面(1 550至2 000道中1 500～2 000 ms，見圖4)進行應用。地震資料經(jīng)過瞬時譜分析技術處理后常出現(xiàn)含氣區(qū)下方振幅異常的低頻陰影現(xiàn)象，低頻陰影現(xiàn)象是含氣檢測的主要依據(jù)。用基于經(jīng)驗小波變換的瞬時譜分析技術處理地震資料，提取3～7 Hz低頻、13～17 Hz中頻以及23～27 Hz高頻的瞬時譜剖面，結果如圖5所示。圖5(a)3～7 Hz低頻瞬時譜剖面，在含氣儲層中信號能量迅速衰減，在其下方(紅色橢圓內(nèi))存在強大的能量分布，出現(xiàn)低頻陰影現(xiàn)象；圖5(c)23～27 Hz高頻瞬時譜剖面，當頻率增大時，含氣儲層中能量增大，其下方的低頻陰影現(xiàn)象消失。結合低頻和高頻瞬時譜剖面，對比能量分布差異可以很好地確定含氣區(qū)域。

圖4 中國某地區(qū)實際資料疊后地震剖面Fig.4 Post-stack seismic profile of actual data in a region of China

注：黃色橢圓表示儲層位置，紅色橢圓表示低頻陰影處，紅色豎線表示A井所在位置。圖5 基于經(jīng)驗小波變換的瞬時譜分析技術處理得到的瞬時譜剖面Fig.5 Instantaneous spectrum section processed by instantaneous spectrum analysis technique based on empirical wavelet transform

4 結束語

經(jīng)驗小波變換是一種全新的自適應分析算法，具有很高的自適應性，其時頻分析解決了模態(tài)混疊問題，能夠有效地去除線性噪音，也不易丟失數(shù)據(jù)。本文將經(jīng)驗小波變換應用到瞬時譜分析技術中，得到了一種新的瞬時譜分析技術，與基于傳統(tǒng)時頻分析方法的瞬時譜分析技術相比解決了Heisenberg測不準原理的制約等難題，具有更高的時間和頻率分辨率。本文將基于經(jīng)驗小波變換的瞬時譜分析技術應用于實際致密砂巖氣藏資料的含氣檢測中，一維地震道分析效果優(yōu)于其他方法，二維地震道分析中低頻陰影更加明顯，低頻與高頻的瞬時譜剖面結合能夠更好地反映地層細節(jié)，更有效地識別致密砂巖氣藏的含氣巖層。