經(jīng)驗(yàn)?zāi)B(tài)分解法在地震資料處理中的應(yīng)用

李世濤 張勇 范存輝 梁亮

(1.中海油能源發(fā)展股份有限公司上海采油技術(shù)服務(wù)分公司，上海 200030;2.西南石油大學(xué)油氣藏地質(zhì)及開(kāi)發(fā)工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，成都 610500)

小波變換在石油勘探中已被廣泛應(yīng)用，它本質(zhì)上是一種自適應(yīng)加窗富氏變換，要求小波窗內(nèi)的信號(hào)短時(shí)平穩(wěn)，沒(méi)有擺脫富氏變換的局限性;其次，小波基的多樣性及零相位不適應(yīng)非零相位子波生成的信號(hào)，即小波基的有限長(zhǎng)度還會(huì)造成信號(hào)能量的泄漏:因此，基于小波變換的時(shí)頻分析很難適應(yīng)非平穩(wěn)時(shí)間序列的時(shí)頻分析。為此，本次研究應(yīng)用Huang等人(1996)提出的經(jīng)驗(yàn)?zāi)B(tài)分解法(Empirical Mode Decomposition，簡(jiǎn)稱 EMD)中的固有模態(tài)函數(shù)(Intrinsic Mode Function，簡(jiǎn)稱IMF)作為基函數(shù)，對(duì)非平穩(wěn)時(shí)間序列的時(shí)頻分析方法進(jìn)行深入研究，把該研究成果應(yīng)用于地震和測(cè)井資料聯(lián)合處理解釋之中。

1 經(jīng)驗(yàn)?zāi)B(tài)分解法時(shí)頻分析

經(jīng)驗(yàn)?zāi)B(tài)分解是對(duì)一種非線性、非平穩(wěn)信號(hào)的新的時(shí)頻分析方法。它能夠?qū)?shù)據(jù)進(jìn)行線性化、平穩(wěn)化處理，并在處理過(guò)程中保留信號(hào)自身特性;它利用非平穩(wěn)信號(hào)的統(tǒng)計(jì)特征尺度隨時(shí)間信號(hào)變化的局部時(shí)變特征，從原信號(hào)中提取若干個(gè)固有模態(tài)函數(shù)(IMF)和一個(gè)殘余量，用多個(gè)固有模態(tài)函數(shù)表征信號(hào)的局部特征，用殘余分量表征信號(hào)緩慢變化的趨勢(shì)。固有模態(tài)函數(shù)需滿足2個(gè)條件:(1)整個(gè)數(shù)據(jù)集合中，極值點(diǎn)數(shù)目和過(guò)0點(diǎn)數(shù)目必須相等或最多相差一個(gè);(2)在任何一點(diǎn)局部區(qū)間，由局部最大值和最小值形成的包絡(luò)均值為0。下面介紹非平穩(wěn)信號(hào)經(jīng)驗(yàn)?zāi)B(tài)分解算法。

(1)EMD分解步驟，分解過(guò)程采用以下“篩分”算法。

第1步，對(duì)信號(hào)x(t)求取所有的局部極值點(diǎn)，用樣條函數(shù)分別擬合所有的局部極大值點(diǎn)和極小值點(diǎn)，形成上下包絡(luò)線bm(t)和bn(t)。

第2步，計(jì)算上下包絡(luò)線的均值:

并將原始序列x(t)減去包絡(luò)均值:

此時(shí)判別h1，1(t)是否滿足IMF的2個(gè)條件。如果不滿足，將h1，1(t)作為原始數(shù)據(jù)重復(fù)上述處理過(guò)程，直至序列

滿足 IMF的條件。這樣，就得到 IMF的第一個(gè)分量:

第3步，由信號(hào)x(t)分離出第一個(gè)分量c1(t)，得剩余量:

將r1(t)作為一個(gè)新的原始序列，重復(fù)上述步驟，依次求出第2，3，…，n個(gè) IMF ci(t)，當(dāng)剩余量rn(t)成為一個(gè)單調(diào)函數(shù)或小于預(yù)定值時(shí)分解結(jié)束。

這樣，原始信號(hào)x(t)可表示為:

(2)篩分終止條件及樣條函數(shù)的選擇。在實(shí)際IMF求取時(shí)，IMF的第2個(gè)條件很難滿足。為此，應(yīng)用Huang等人提出的收斂準(zhǔn)則，定義:

當(dāng)Is≤ξ時(shí)，則終止篩分。對(duì)局部求取的上下極值點(diǎn)用樣條插值函數(shù)插值組成包絡(luò)。插值函數(shù)的邊界端點(diǎn)處理亦是EMD分解的一個(gè)難題。本次研究運(yùn)用自適應(yīng)外推法進(jìn)行邊界處理。

2 希爾伯特變換和瞬時(shí)譜計(jì)算

由式(6)可知，信號(hào)x(t)經(jīng)EMD分解得到的每一個(gè)IMF的ci(t)均為單分量?，F(xiàn)對(duì)每一個(gè)單分量進(jìn)行富氏變換:

式中:P表示柯西積分主值。由此可得解析信號(hào):

Am(t)、θ(t)為瞬時(shí)振幅和相位，其計(jì)算公式為:

由瞬時(shí)相位可求得瞬時(shí)頻率:

將振幅值用時(shí)間和頻率表示，可得到希爾伯特幅值譜:

將Hm(ω，t)對(duì)時(shí)間積分，得邊間譜:

由希爾伯特幅值譜和邊間譜，可給出信號(hào)序列的平穩(wěn)度函數(shù):

對(duì)應(yīng)時(shí)間域有:

式中:

3 應(yīng)用實(shí)例

3.1 地震資料的EMD分解

圖1展示了地震道xi(t)，EMD分解全過(guò)程包括IMF單分量ck(t)及分解時(shí)產(chǎn)生的逐級(jí)上、下極值包絡(luò)bkn(t)。

圖1 地震道EMD分解全過(guò)程

圖2 圖1中EMD分解后的IMF k

圖2 展示了圖1中EMD分解后的IMFk包絡(luò)bkn和單分量ck結(jié)果 。圖3顯示的是地震道x(t)，通過(guò)EMD分解得到的IMF分量ck(t)及IMF的第2個(gè)條件篩分準(zhǔn)則中的近似包絡(luò)值bkn(t)。IMF分量ck(t)的頻率隨著k的增大逐漸降低，波長(zhǎng)逐漸增大。各個(gè)分量ck(t)包含有不同的時(shí)間尺度及分辨率的信號(hào)特征，這種特征與信號(hào)自身的特點(diǎn)是自適應(yīng)的。其EMD分解后有6個(gè)IMF單分量。以下逐級(jí)IMF分量ck(t)的值及包絡(luò)平均值bkn(t)逐次依0.5的幾何級(jí)數(shù)衰減，為方便展示，全部作歸一化處理。圖3表明，具有高頻成分的主分量c1(t)的平穩(wěn)性較差，具有低頻成分的平穩(wěn)性較好。

圖3 圖1中EMD分解后局部放大的IMF k

3.2 地震資料EMD分解IMF主分量分析剖面效果

圖4 是某測(cè)線偏移剖面圖，圖5展示了該地震剖面EMD分解后得到的主分量、次分量及第3分量剖面。圖6是某測(cè)線井段附近剖面的EMD分解主分量剖面與小波分解剖面圖，展示了主分量剖面與分解的1級(jí)包絡(luò)剖面的關(guān)系:在包絡(luò)剖面上絕對(duì)幅度越大表示該處的平穩(wěn)性越差、地層的非白噪性越強(qiáng);在主分量剖面上，幾處低頻高幅處是儲(chǔ)層及含油氣層，而小波分解剖面上沒(méi)有這種特征。

圖4 某測(cè)線偏移剖面

圖5 某測(cè)線EMD單分量IMF剖面

圖6 井段剖面EMD分解與小波分解剖面

4 地震資料EMD分解后單分量IMF k的時(shí)頻分析

地震資料x(chóng)(t)在經(jīng)過(guò)EMD分解后獲得近似平穩(wěn)的固有模態(tài)函數(shù)IMFk，個(gè)數(shù)比較少只有6～7個(gè)，很難直接應(yīng)用這么少量的IMF作時(shí)頻分析。因此采用基于小波基的時(shí)頻分析方法對(duì)逐個(gè)IMFk作時(shí)頻分析的有益嘗試，取得了一定成效。

圖7是地震道EMD及連續(xù)小波變換CWT分解后所作的瞬時(shí)Hilbert幅值譜對(duì)比圖。兩者對(duì)比表明EMD幅值譜較CWT幅值譜能量強(qiáng)的特點(diǎn)。圖8是地震道EMD分解后計(jì)算得到的時(shí)域和頻域的平穩(wěn)度曲線。

由圖7—圖9可知，EMD分解后所作的瞬時(shí)Hilbert幅值譜比連續(xù)小波變換CWT分解后的譜的分辨率高，并與地震道波形匹配自然。

圖7 瞬時(shí)Hilbert幅值譜對(duì)比

圖8 地震道EMD分解的平穩(wěn)度曲線

圖9 地震剖面對(duì)應(yīng)的EMD分解Hilbert 譜

5 結(jié)語(yǔ)

本次研究對(duì)經(jīng)驗(yàn)?zāi)B(tài)分解法在地震數(shù)據(jù)處理領(lǐng)域中的應(yīng)用作了探索論證，并取得一定的應(yīng)用效果。運(yùn)用經(jīng)驗(yàn)?zāi)B(tài)分解法和小波變換2種方法進(jìn)行試驗(yàn)比較，得到4個(gè)方面的結(jié)論。

(1)在信號(hào)分頻處理中，基于小波變換的時(shí)頻分析應(yīng)用方便靈活，具有廣泛的應(yīng)用性，優(yōu)于常規(guī)富氏變換的分頻處理技術(shù);信號(hào)的小波分解是一個(gè)褶積過(guò)程，但不利于非平穩(wěn)數(shù)據(jù)的時(shí)頻分析。

(2)在地震數(shù)據(jù)x(t)的EMD分解中，由于分解得到的IMFk一般很少，而且只有前幾個(gè)分量在地震數(shù)據(jù)時(shí)頻分析的頻段內(nèi)，因此一般較難用于時(shí)頻分析。而EMD的第一個(gè)分量c1(t)具有較寬的頻率成分，因此，與地震解釋配合，EMD方法的主分量可用于地震屬性分析。

(3)由于測(cè)井資料與地震數(shù)據(jù)的相似特性，可在重構(gòu)測(cè)井曲線，儲(chǔ)層參數(shù)反演方面有良好的應(yīng)用前景。

(4)在對(duì)地震資料的EMD分解作時(shí)頻分析時(shí)，插值函數(shù)的邊界端點(diǎn)處理是EMD分解的一個(gè)難題。本次研究采用自適應(yīng)外推法作為邊界條件，但需要進(jìn)一步研究模態(tài)函數(shù)IMFk之間的分形插值方法，使少量的IMFk插成較密的IMFk簇，提高地震資料EMD分解法的時(shí)頻分析效果。

［1］ Huang E N，Wu C M，Long R S，et al.A Confidence Limit for the Empirical Mode Decomposition and Hilbert Spectral Analysis［C］//Proceedings Royal Society of London.2003:2317-2345.

［2］楊世錫，胡勁松，吳昭同，等.基于高次樣條插值的經(jīng)驗(yàn)?zāi)B(tài)分解方法研究［J］.浙江大學(xué)學(xué)報(bào)(工學(xué)版)，2004，3(3):267-270.

［3］程磊，瞿偉廉.基于Hilbert-Huang變換理論的非平穩(wěn)數(shù)據(jù)處理［J］.建筑科學(xué)與工程學(xué)報(bào)，2007，24(1):26-30.

［4］王祝文，劉菁華，聶春燕，等.Hilbert-Huang變換在提取陣列聲波信號(hào)動(dòng)力特性中的應(yīng)用［J］.地球物理學(xué)進(jìn)展，2008，23(2):450-455.

［5］ 俞壽朋.寬帶 Ricker子波［J］.石油地球物理勘探，1996，31(5):605-615.

［6］田芳，高靜懷，陳文超，等.廣義S變換與薄互層地震響應(yīng)［J］.地球物理學(xué)報(bào)，2003，46(4):526-532.

［7］蓋廣洪.經(jīng)驗(yàn)?zāi)B(tài)分解的一種改進(jìn)算法［J］.西安交通大學(xué)學(xué)報(bào)，2004，38(11):1199-1202.