基于經(jīng)驗(yàn)?zāi)B(tài)分解字典的自適應(yīng)匹配追蹤譜分解方法及其在油氣檢測(cè)中的應(yīng)用

潘 輝 印興耀* 李 坤 裴 松

(①中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島266580；②青島海洋科學(xué)與技術(shù)國(guó)家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室，山東青島266071)

0 引言

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

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

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

1 EMD字典

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

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

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

(1)

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

(2)

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

(3)

一個(gè)二維地震信號(hào)s(t)經(jīng)過(guò)EMD后，表示為

(4)

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

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

2 基于EMD字典的自適應(yīng)MP算法

2.1 快速M(fèi)P基本原理

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

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

(5)

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

S=〈S,gy1〉gy1+R1S

(6)

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

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

(7)

其中g(shù)yn+1滿足

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

(8)

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

(9)

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

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

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

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

連續(xù)信號(hào)x(t)的復(fù)地震道為

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

(10)

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

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

θ=P(φ)=φ+rπ

(11)

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

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

(12)

Q為求主值相位的算子，即P是Q的逆運(yùn)算

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

(13)

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

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

(14)

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

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

(15)

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

(16)

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

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

(17)

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

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

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

圖1 由不同方法獲得的合成信號(hào)瞬時(shí)頻率

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

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

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

2.3 改進(jìn)MP算法

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

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

本文搜索方法的每次迭代過(guò)程可簡(jiǎn)單表示為

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

i=1，2，…，M

(18)

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

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

(19)

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

a=[GTG+εI]-1gTc

(20)

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

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

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

圖5 無(wú)噪聲條件下理論信號(hào)的迭代分解過(guò)程

圖6 不同時(shí)頻分析方法得到的時(shí)頻譜

3 模型試算

3.1 一維模型試算

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

圖7 多尺度地震信號(hào)

圖8 圖7的重構(gòu)信號(hào)殘差

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

圖10 圖9的重構(gòu)信號(hào)殘差

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

圖12 圖11的重構(gòu)信號(hào)殘差

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

圖14 圖13的重構(gòu)信號(hào)殘差

3.2 二維模型試算

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

圖15 由改進(jìn)MP算法重構(gòu)的M區(qū)地震剖面

圖16 圖15的重構(gòu)剖面殘差

4 實(shí)際資料處理

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

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

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

圖18 基于EMD改進(jìn)MP追蹤與傳統(tǒng)三參數(shù)的動(dòng)態(tài)快速M(fèi)P方法計(jì)算時(shí)間對(duì)比

5 結(jié)論

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

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

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

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