潘 輝 印興耀* 李 坤 裴 松

(①中國石油大學(華東)地球科學與技術學院，山東青島266580；②青島海洋科學與技術國家實驗室海洋礦產(chǎn)資源評價與探測技術功能實驗室，山東青島266071)

0 引言

經(jīng)驗模態(tài)分解(EMD)是由Huang等[1]提出的一種自適應處理非線性和非平穩(wěn)信號的方法。EMD將隨機信號分解為固有模式信號，賦予瞬時頻率合理的物理含義和計算方法，創(chuàng)立了以瞬時頻率表征信號變化的基本屬性[2]。隨著圖像處理技術的發(fā)展，EMD在圖像處理領域取得了巨大成功，如圖像去噪、檢測、增強等[3-5]，但在地震勘探領域應用較少。

近年來，匹配追蹤(Matching Pursuit，MP)算法廣泛用于地震勘探領域，如檢測低頻陰影、提取瞬時譜、稀疏反演、去強反射等[6-10]，其中譜分解技術可全頻帶掃描地震數(shù)據(jù)，獲得特定頻率下描述目標層地質(zhì)特征的數(shù)據(jù)體，普遍用于巖性識別及烴類檢測[11-12]。傳統(tǒng)的譜分解技術(如短時傅里葉變換、小波變換和S變換等)可以有效地分析非平穩(wěn)信號，但存在物理意義不明確、時頻譜精度不高等問題。Mallat等[13]提出的MP算法根據(jù)信號的局部結構特征構造超完備字典原子庫，并將信號擴展到一組適應不同時頻特征的原子上，最后用有限個原子稀疏表示原始信號，實現(xiàn)信號的自適應分解。Liu等[14-15]推出了動態(tài)快速MP算法的概念，將信號的瞬時屬性引入MP中，在每次迭代中，信號的瞬時屬性變化引起參與內(nèi)積的時頻原子的變化，進而提高了信號的分解速率。張繁昌等[16-17]充分利用地震信號的瞬時屬性確定時頻原子動態(tài)參數(shù)的掃描范圍，最終采用阻尼最小二乘法確定各個匹配小波的復振幅，進一步減少了掃描參數(shù)，計算效率提高了數(shù)百倍，但獲得的瞬時頻率穩(wěn)定性較差且存在負值。為此，Boashash[18]、劉漢卿等[19]利用連續(xù)相位求取瞬時頻率，但獲得的瞬時頻率易受噪聲影響，并存在頻率異常值。印興耀等[20]引入局部頻率約束時頻原子的搜尋范圍，但得到的局部頻率與實際頻率存在一定誤差。

為此，本文針對瞬時頻率存在負值的情況展開了一系列研究。瞬時頻率是描述單分量信號的頻率隨時間變化的物理屬性，一般的地震信號為多分量信號集合體，在某一個時刻具有諸多頻率成分，該時刻的地震信號的瞬時頻率實際上是該時刻所有頻率的平均值。故對原始信號進行EMD得到單分量信號，引入連續(xù)相位，利用阻尼最小二乘法求出匹配原子的瞬時頻率，由此得到匹配原子主頻；采用改進的動態(tài)MP方法分解地震信號穩(wěn)定性較好，得到的高分辨率單頻譜剖面有助于確定含油氣儲層的范圍。

1 EMD字典

EMD將非平穩(wěn)信號自適應地分解成一系列固有模式信號(IMF)，IMF滿足單分量信號的物理解釋且含有不同的頻率成分，這些IMF包括高頻到低頻成分，能夠表征不同的地質(zhì)、地層的信息。IMF需要滿足兩個條件[1]：①信號零值點與極值點的數(shù)量頂多相差一個；②在任意時刻，信號局部的極大值點和極小值點形成的上、下包絡的平均值等于零。EMD依賴于“篩選”的過程提取IMF，從稀疏分解的角度看，EMD算法的結果非常稀疏，在分解過程中子成分分量只有若干個，遠遠少于其他稀疏時頻分解方法(如MP)。因此，本文分解方法的時頻字典(EMD字典)比其他時頻字典(Morlet小波字典等)更具冗余性。根據(jù)EMD算法獲得的所有可能的IMF的集合定義EMD字典，每一個IMF分量都是平穩(wěn)的。根據(jù)IMF的定義，EMD字典中可以有無數(shù)個時頻原子，Hou等[21]從子空間投影的角度給出了IMF的表達式

Dimf(t)=a(t)cosθ(t)

a(t)、θ′(t)∈V(θ,λ)θ′(t)≥0

(1)

式中：Dimf(t)為經(jīng)EMD的IMF函數(shù)；a(t)為瞬時振幅，t為時間；θ(t)為瞬時相位(遞增函數(shù))，其導數(shù)θ′(t)可以表示瞬時頻率。由于IMF的前提假設是平穩(wěn)的，故信號的能量和瞬時頻率不隨時間發(fā)生劇烈變化，因此a(t)和θ′(t)屬于由θ(t)定義的諧波列向量所組成的子空間V(θ,λ)

(2)

通過使用上述IMF作為超完備字典中的原子，即可定義EMD字典

(3)

一個二維地震信號s(t)經(jīng)過EMD后，表示為

(4)

式中：IMFm(t)為第m個IMF；Rm(t)為殘余量。

EMD過程是極具啟發(fā)性的，雖然是基于經(jīng)驗的分解，不同于傳統(tǒng)的稀疏時頻分解方式，但仍可以視作基于EMD字典逐步搜索最優(yōu)原子的過程。

2 基于EMD字典的自適應MP算法

2.1 快速MP基本原理

Mallat等[13]提出的MP算法的核心是：根據(jù)信號的特性選定母函數(shù)創(chuàng)建一個冗余字典，將信號投影到原子庫的所有原子上，選擇與信號相關系數(shù)最大的原子作為最優(yōu)匹配原子，并從原始信號中減去；對剩余信號繼續(xù)上述操作，記錄每次迭代產(chǎn)生的最優(yōu)匹配子波，繼續(xù)迭代下去，直到剩余信號的能量小于閾值或達到預先設置的最大迭代次數(shù)，則信號分解完成，原始信號可以由分解得到的一系列最優(yōu)匹配原子表征。

|〈S,gy1〉|=supi∈(1,…,l)|〈S,gyi〉|

(5)

式中：〈S,gyi〉為最優(yōu)時頻原子與信號的內(nèi)積，yi是字典矩陣的列索引，i為D的列索引；sup表示取|〈S,gyi〉|的上確界。這樣，S由第一次迭代搜尋獲得的最優(yōu)匹配原子gy1的投影分量及其信號殘差R1S兩部分組成，即

S=〈S,gy1〉gy1+R1S

(6)

將R1S作為新信號，重新搜索原子庫尋找最優(yōu)匹配原子，并從原始信號中減去，得到剩余信號，繼續(xù)迭代下去。到第n+1次迭代時，經(jīng)過匹配分解后殘存信號為Rn+1S，搜索得到的最優(yōu)匹配原子為gyn+1，則第n次迭代獲得的殘存信號RnS為

RnS=〈RnS,gyn+1〉gyn+1+Rn+1S

(7)

其中gyn+1滿足

|〈RnS,gyn+1〉|=supi∈(1,2,…，l)|〈RnS,gyi〉|

(8)

若經(jīng)過K步迭代、分解后，匹配子波的次數(shù)達到預先設置的次數(shù)或迭代、分解后信號殘存能量遠低于能量閾值時，就完成對信號的分解。原始信號可以表示為K個最優(yōu)匹配原子的線性組合與殘差RK+1S之和

(9)

2.2 基于連續(xù)相位阻尼最小二乘法求取瞬時頻率

傳統(tǒng)的動態(tài)MP方法以信號的瞬時屬性作為先驗信息構建動態(tài)子波庫，顯著提高了每次迭代的內(nèi)積計算速度，極大地加快了分解效率。傳統(tǒng)MP的最優(yōu)子波的搜索鄰域由瞬時屬性確定，瞬時頻率計算結果存在“負頻率”的問題，且受噪聲干擾較大，得到的結果極不穩(wěn)定。當獲得的頻率不合理時，只能全局搜尋時頻原子的頻率，這并不是真正意義上的動態(tài)搜尋。為此，引入連續(xù)相位概念計算瞬時頻率[22-23]。根據(jù)正則化理論，利用正則化算子得到的局部域值求取數(shù)據(jù)點處的頻率信息，進而求得準確的信號瞬時頻率屬性信息，并且計算結果穩(wěn)定，抗噪聲能力強。因此，本文提出的瞬時相位計算方法可以直接找到最優(yōu)時頻原子。

自Taner 等[24]初次將瞬時頻率引入地震勘探領域之后，瞬時頻率成為基本地震屬性之一。人們又提出了求取瞬時頻率的優(yōu)化算法。高靜懷等[25]提出基于小波變換域的瞬時頻率計算方式；尹繼堯等[26]基于TK能量的最大幅值計算瞬時頻率；劉漢卿等[19]采用連續(xù)相位求取瞬時頻率；印興耀等[20]利用局部頻率替代瞬時頻率去除強反射，獲得了較好效果。

針對傳統(tǒng)MP方法計算瞬時頻率出現(xiàn)負值及計算過程魯棒性差的問題，本文引入連續(xù)相位，借助阻尼最小二乘法，加入整形正則化算子[27]對數(shù)據(jù)平滑處理，得到了穩(wěn)定的瞬時頻率，并獲得匹配原子主頻。采用改進的動態(tài)MP算法對地震信號譜分解，具有良好的穩(wěn)定性，獲得的高分辨率單頻譜剖面有助于確定含油氣儲層的范圍。

連續(xù)信號x(t)的復地震道為

Z(t)=x(t)+jh(t)=A(t)ejθ(t)

(10)

式中：h(t)為x(t)的Hilbert變換；A(t)為地震道包絡；θ(t)為信號的瞬時相位。

設θ和φ分別為連續(xù)相位和主值相位，P、Q分別為計算θ、φ的算子，則連續(xù)相位定義為

θ=P(φ)=φ+rπ

(11)

式中r為正整數(shù)。由式(11)可得

Δθ(t)=θ(t+1)-θ(t)∈(0，π]

(12)

Q為求主值相位的算子，即P是Q的逆運算

φ=Q(θ)=P-1(θ)

(13)

Q(θ)位于區(qū)間(-π，π]，則

Δθ(j)=P{Q[Δφ(j)]}

(14)

將式(14)代入式(12)可得連續(xù)相位的計算公式

θ(t+1)=P{Q[Δφ(t)]}+θ(t)

(15)

傳統(tǒng)的瞬時頻率f(t)是θ(t)的變化率，即

(16)

用連續(xù)相位替代瞬時相位φ(t)，引入阻尼最小二乘法，并加入整形正則化算子求取瞬時頻率

f(t)=[e2I+(MB)T(MB)]-1(MB)Tl

(17)

式中：e為權系數(shù)，一般取B中元素最大值的1%～5%；M為整形正則化算子；l和B分別為l(t)和b(t)組成的向量(矩陣)；I為單位陣。

本文基于連續(xù)相位求解的瞬時頻率有效避免了“負頻率”現(xiàn)象，相較于直接求逆過程，利用整形正則化平滑算子的阻尼最小二乘法不會出現(xiàn)頻率異常值，得到的瞬時頻率曲線更平滑、真實，且更突顯高頻成分。

為了驗證瞬時頻率計算結果的合理性，設計了合成信號(圖1a)及加噪合成信號(圖2a、圖3a、圖4a)，并由不同方法計算瞬時頻率(圖1～圖4)。由圖1可見：①瞬時頻率求導計算結果出現(xiàn)很多“負頻率”現(xiàn)象(圖1b)，連續(xù)相位求導計算結果雖然避免了“負頻率”現(xiàn)象，但在0.4～0.5s出現(xiàn)瞬時頻率異常(圖1c)。②連續(xù)相位常規(guī)阻尼最小二乘法(圖1d)、連續(xù)相位改進阻尼最小二乘法(圖1f)計算結果沒有出現(xiàn)瞬時頻率異常，且后者的計算結果更平滑，并且突顯了頻率異常點。③局部頻率(圖1e)與連續(xù)相位改進阻尼最小二乘法(圖1f)計算結果接近真實頻率，并且突顯了高頻成分，但后者更接近真實頻率。由圖2～圖4可見，基于連續(xù)相位改進阻尼最小二乘法(圖2f、圖3f、圖4f)能很好地處理加噪地震數(shù)據(jù)，計算結果相對于局部頻率曲線(圖2e、圖3e、圖4e)更光滑，且突顯了頻率異常點，證明本文方法具有很好的抗噪性。

圖1 由不同方法獲得的合成信號瞬時頻率

圖2 由不同方法獲得的加噪合成信號(信噪比為1)瞬時頻率

圖3 由不同方法獲得的加噪合成信號(信噪比為2)瞬時頻率

圖4 由不同方法獲得的加噪合成信號(信噪比為5)瞬時頻率

2.3 改進MP算法

傳統(tǒng)動態(tài)MP方法將信號轉換到Hilbert空間，并在信號包絡幅值最大值處展開、分解信號。首先對信號在先驗信息確定的搜尋鄰域內(nèi)的所有原子進行匹配，找出相關系數(shù)最大的時頻原子作為最優(yōu)原子，這種分解策略造成每一次迭代歷經(jīng)搜尋范圍內(nèi)所有的原子，僅獲得一個最優(yōu)匹配的時頻原子，影響了信號的收斂速度[28-29]。為了加快MP算法的收斂速度，采用多原子分解方法，即不僅僅在包絡幅值最大值處進行，而是找到所有的信號能量極大值位置，并設置特定的搜索條件，在滿足搜索條件的時間點處進行匹配[30]，通過一次迭代得到多個最優(yōu)匹配子波，并利用最小二乘法求出每個子波的幅值參數(shù)。

多原子分解的動態(tài)MP算法需要在每一次迭代過程中搜索滿足一定條件的極大值附近的所有原子，從而獲得最優(yōu)匹配的時頻原子，搜索條件的限制嚴重降低了信號的分解速度。故本文引入EMD思想，將每個地震波形自適應地分解為多個不同波形的組合，每一次迭代只在所有信號包絡幅值極大值處分解，當信號的殘存能量與原始信號總能量的比值非常小時停止分解。在實際的地震信號分析中，當殘存能量遠小于原始信號總能量時，可以忽略殘存量。故當殘留能量小于原始信號總能量的千分之一時，文中MP停止迭代，能在一次迭代中快速生成多個匹配子波，極大提高了分解效率和分解精度。

本文搜索方法的每次迭代過程可簡單表示為

gyi={T0,f∈U[f(T0)],θ∈U[θ(T0)]}

i=1，2，…，M

(18)

式中：M表示每一次迭代搜尋確定的最優(yōu)匹配子波的個數(shù)；T0為每個最優(yōu)匹配子波在搜索位置處的時間；U[f(T0)]和U[θ(T0)]分別為在搜索位置處的頻率與相位的搜索鄰域集合，f(T0)和θ(T0)分別為搜索位置處信號的瞬時頻率信息和相位信息。

經(jīng)過n次、n+1次匹配分解后，分別得到信號的殘差RnS、Rn+1S，假定第n+1次迭代搜索到M個符合條件的匹配原子gyi(i=1，2，…，M)，則

(19)

式中ai(i=1，2，…，M)為第n+1次迭代搜索獲得的每個匹配子波的幅值，即

a=[GTG+εI]-1gTc

(20)

式中：a=(a1,a2,…,aM)T為采用多原子搜尋策略通過一次迭代確定的M個時頻原子的振幅構成的列向量；c代表第n+1次迭代后的Rn+1S構成的N維列向量，N≥M為信號列向量的元素個數(shù)；G=[gy1,gy2,…,gyM]為最優(yōu)原子組成的N×M階原子庫矩陣，其中每個元素表征一個最優(yōu)時頻原子；I為M階單位矩陣；ε為阻尼因子。

利用11個Morlet小波合成理論信號，采用基于EMD字典的改進MP算法分解理論信號。圖5為無噪聲條件下理論信號的迭代分解過程。由圖可見，每一次迭代都篩選出信號的所有極大值點，并逐個匹配原子，經(jīng)過11次迭代后，重構信號與理論信號波形基本吻合，重構殘差幾乎為零。

為進一步驗證聯(lián)合EMD與連續(xù)相位求解瞬時頻率方法的實用性和穩(wěn)定性，對圖5的合成理論信號利用基于EMD字典的匹配追蹤Wigner-Ville分布(EMP-WVD)進行時頻表征，并對比、分析短時傅里葉變換(STFT)、連續(xù)小波變換(CWT)、S變換(ST)以及EMP-WVD的時頻譜計算結果(圖6)?？梢姡珽MP-WVD時頻譜(圖6d)的時頻分辨率遠高于常規(guī)時頻分析方法(圖6a～圖6c)，時頻能量聚集能力更強，能清晰定位原子時頻信息，為儲層預測奠定了基礎。

圖5 無噪聲條件下理論信號的迭代分解過程

圖6 不同時頻分析方法得到的時頻譜

3 模型試算

3.1 一維模型試算

為驗證基于EMD和連續(xù)相位求取瞬時頻率方法的抗噪性，利用M區(qū)A井的測井資料建立了一維模型。圖7為多尺度地震信號，圖8為實際信號與重構信號的殘差。由圖可見：EMD從地震信號最高頻率分量開始，依次分解出頻率范圍逐漸降低的IMF成分，可以有效獲得低頻和高頻成分(圖7)；殘差基本為0(圖8)。圖9～圖14分別為多尺度加噪地震信號及其重構信號殘差。由圖可見，殘差基本為0(圖10、圖12、圖14)，表明本文方法能夠完全重構原始地震信號，具有良好的抗噪性。

圖7 多尺度地震信號

圖8 圖7的重構信號殘差

圖9 多尺度加噪地震信號(信噪比為1)

圖10 圖9的重構信號殘差

圖11 多尺度加噪地震信號(信噪比為2)

圖12 圖11的重構信號殘差

圖13 多尺度加噪地震信號(信噪比為5)

圖14 圖13的重構信號殘差

3.2 二維模型試算

針對復雜薄互層二維模型開展基于EMD改進MP算法抗噪性測試。圖15為由改進MP算法重構的M區(qū)地震剖面，圖16為圖15的重構剖面殘差。由圖可見：基于EMD字典的改進MP方法的重構剖面(圖15b)與實際地震剖面(圖15a)極相似，但前者分辨率得較高，殘差(圖16a)只剩隨機噪聲；對地震記錄加噪后，殘差基本為隨機噪聲(圖16b～圖16d)，且重構剖面的分辨率均有提升(圖15d、圖15f、圖15h)，進一步驗證了本文方法的可行性。

圖15 由改進MP算法重構的M區(qū)地震剖面

圖16 圖15的重構剖面殘差

4 實際資料處理

在對一維及二維模型測試的基礎上，為進一步驗證文中方法對實際資料的應用效果，利用中國N區(qū)地震資料進行測試。圖17為N區(qū)多尺度瞬時譜剖面。由圖可見：①EMD是從信號最高頻率成分開始，依次分解出頻率范圍逐漸降低的IMF分量(圖17b～圖17d)，因此可以有效獲得地震信號的低頻和高頻成分。②當頻率較低時，儲層底部能量很強，當頻率較高時，儲層底部出現(xiàn)能量減弱現(xiàn)象，當頻率再次增大時，儲層下方能量基本消失，能夠觀測到明顯的低頻陰影現(xiàn)象。因此，聯(lián)合EMD和改進MP譜分解技術檢測低頻陰影的效果更好，結合測井資料可更好地指示油氣儲層。

圖17 N區(qū)多尺度瞬時譜剖面

圖18為基于EMD改進MP追蹤與傳統(tǒng)三參數(shù)的動態(tài)快速MP方法計算時間對比。由圖可見，本文提出的基于EMD的自適應MP方法顯著提升了MP分解效率。

圖18 基于EMD改進MP追蹤與傳統(tǒng)三參數(shù)的動態(tài)快速MP方法計算時間對比

5 結論

本文結合EMD稀疏分解與MP算法，提出了一種基于EMD字典的稀疏時頻分解算法。此外，在連續(xù)相位的阻尼最小二乘反演求解瞬時頻率的方法中加入整形平滑算子約束，有效地避免了奇異值，同時進一步提高了MP計算效率。經(jīng)過理論分析與實驗測試得到以下認識：

(1)基于整形平滑算子約束的連續(xù)相位求解瞬時頻率方法的計算過程并不依賴于每個數(shù)據(jù)點的瞬時值，而是基于數(shù)據(jù)點的局部鄰域值求取的。因此即使在部分信號微弱或缺失的情況下，該方法也能從相鄰時窗中提取有效信息，從而得到合理的頻率值，有效降低了噪聲敏感度。

(2)基于EMD字典的MP算法，在確定動態(tài)快速MP先驗信息的搜索范圍時，以連續(xù)相位替換瞬時相位求解瞬時頻率，提高了運算效率，但繼承了EMD的端點效應，在重構地震數(shù)據(jù)時，未能較好地重構端點數(shù)據(jù)。

(3)本文將基于EMD字典的快速MP方法應用于儲層含油氣性預測，對比不同尺度瞬時譜剖面可以清楚地看到明顯的低頻陰影現(xiàn)象，且提高了分解效率，從而進一步驗證了方法可行性。