基于EMD時(shí)頻分析的聲波測(cè)井油層的信息提取

張 婷，聶春燕

(1.長(zhǎng)春理工大學(xué) 電子信息工程學(xué)院，長(zhǎng)春 130022;2.長(zhǎng)春大學(xué) 電子信息工程學(xué)院，長(zhǎng)春 130022)

基于EMD時(shí)頻分析的聲波測(cè)井油層的信息提取

張 婷1，聶春燕2

(1.長(zhǎng)春理工大學(xué) 電子信息工程學(xué)院，長(zhǎng)春 130022;2.長(zhǎng)春大學(xué) 電子信息工程學(xué)院，長(zhǎng)春 130022)

陣列聲波測(cè)井信號(hào)是典型的非線性、非平穩(wěn)信號(hào)。文中采用EMD(經(jīng)驗(yàn)?zāi)B(tài)分解)的時(shí)頻分析方法，對(duì)油層聲波信息提取儲(chǔ)集層性質(zhì)。首先對(duì)信號(hào)進(jìn)行EMD分解，得到有限個(gè)固有模態(tài)函數(shù)(IMF)，再次對(duì)每個(gè)IMF做Hilbert變換，求得信號(hào)的Hilbert譜的三維分布以及Hilbert邊際譜和瞬時(shí)能量譜，仿真結(jié)果表明，油層中縱波能量最小，橫波能量較大，一般大于斯通利波的能量。斯通利波能量衰減嚴(yán)重。

陣列聲波;經(jīng)驗(yàn)?zāi)B(tài)分解;時(shí)頻分析;Hilbert變換

0 引言

聲波測(cè)井技術(shù)經(jīng)歷了幾十年的發(fā)展，已經(jīng)成為地球物理測(cè)井學(xué)科的重要領(lǐng)域，是發(fā)展最快和應(yīng)用最廣泛的現(xiàn)代測(cè)井方法之一。隨著聲波測(cè)井技術(shù)的發(fā)展，對(duì)聲波測(cè)井處理的解釋的要求也越來(lái)越高。傳統(tǒng)方法的局限性和不適應(yīng)性暴露的越來(lái)越明顯，如傅立葉變換不能作局部分析，無(wú)法由函數(shù)的傅里葉變換知道函數(shù)在任一點(diǎn)處的狀態(tài)［1］。長(zhǎng)期以來(lái)信號(hào)處理的對(duì)象局限于確定性信號(hào)或是統(tǒng)計(jì)量不隨時(shí)間變化的平穩(wěn)信號(hào)，但非平穩(wěn)信號(hào)的廣泛存在是不爭(zhēng)的事實(shí)，尤其陣列聲波測(cè)井信號(hào)都是典型的非平穩(wěn)信號(hào)。由于受到信號(hào)處理理論的發(fā)展限制，對(duì)非平穩(wěn)信號(hào)的分析過(guò)去人們一直是沿用平穩(wěn)信號(hào)的處理方法來(lái)作近似，效果不夠理想。時(shí)頻分析將非平穩(wěn)的陣列聲波測(cè)井信號(hào)映射到時(shí)頻的二維、三維空間，可以很好地表示出信號(hào)的頻率成分隨時(shí)間的變化規(guī)律。

1 陣列聲波測(cè)井信號(hào)的時(shí)頻分析機(jī)理

時(shí)頻分析(JTFA)即時(shí)頻聯(lián)合域分析(Joint Time-Frequency Analysis)的簡(jiǎn)稱，作為分析時(shí)變非平穩(wěn)信號(hào)的有力工具，成為現(xiàn)代信號(hào)處理研究的一個(gè)熱點(diǎn)。

時(shí)頻分析旨在構(gòu)造一種時(shí)間和頻率的聯(lián)合密度函數(shù)，以揭示信號(hào)中所包含的頻率分量及其演化特性，最終目的是要建立一種分布，以便能在時(shí)間和頻率上同時(shí)表示信號(hào)的能量或者強(qiáng)度，得到這種分布后，就可以對(duì)信號(hào)進(jìn)行分析、處理，提取信號(hào)中所包含的特征信息，或者綜合得到具有期望的時(shí)頻分布特征信號(hào)。

陣列聲波測(cè)井信號(hào)與一般的信號(hào)一樣，是有層次結(jié)構(gòu)的，層次結(jié)構(gòu)越大越復(fù)雜，其間的結(jié)合能亦越大，要分離開(kāi)這些層次結(jié)構(gòu)，就需要更復(fù)雜的算法和更大的計(jì)算量?；贓MD的時(shí)頻分析為這提供了有利的技術(shù)手段。許多情形下，如井眼不規(guī)則、層界面、裂縫以及非均質(zhì)各向異性等導(dǎo)致的波形畸變和異常時(shí)，通過(guò)計(jì)算機(jī)自動(dòng)識(shí)別并計(jì)算這些參數(shù)會(huì)有許多困難。另外縱波與橫波的相對(duì)幅度比、低頻斯通利波與高頻偽瑞利波的相對(duì)幅度比、縱波主頻、斯通利波主頻、斯通利波反射系數(shù)等有待進(jìn)一步開(kāi)發(fā)利用。在某些情況下，這些新參數(shù)比聲速或衰減等常規(guī)參數(shù)更有效。在實(shí)際應(yīng)用中當(dāng)?shù)貙訔l件比較復(fù)雜導(dǎo)致波形變化異常時(shí)，利用陣列聲波波形計(jì)算幅度衰減會(huì)有較大困難［2］。本文采用EMD的時(shí)頻分析方法，計(jì)算同一深度點(diǎn)測(cè)量得到的波形在不同時(shí)間、不同頻率的相對(duì)幅度變化情況，充分利用陣列聲波測(cè)井所獲取的波列信息，有效提取信號(hào)中的特征。

2 基于EMD的信號(hào)處理方法

為把信號(hào)的頻譜分析精確到每一個(gè)時(shí)間點(diǎn)，頻率點(diǎn)上，美國(guó)學(xué)者N.E.Huang等提出了經(jīng)驗(yàn)?zāi)B(tài)分解法(EMD)?；诮?jīng)驗(yàn)的模態(tài)分解方法，將一個(gè)時(shí)間序列信號(hào)分解成有限個(gè)不同時(shí)間尺度的固有模態(tài)函數(shù)IMF(intrinsic mode function)，然后將每個(gè)IMF進(jìn)行Hilbert Huang變換，得到時(shí)頻平面上的能量分布譜圖，用來(lái)分析信號(hào)時(shí)頻譜特征的信號(hào)分析方法［3］。

EMD算法［4－6］通過(guò)不斷剔除信號(hào)的極大值和極小值連接的上下包絡(luò)曲線的均值，將原始信號(hào)X(t)分解成如下形式［4］:

(1)式中，Ci為每一個(gè)IMF，rn為殘余函數(shù)，代表信號(hào)的平均趨勢(shì)。EMD的分解過(guò)程其實(shí)就是一個(gè)“篩分”過(guò)程，在這個(gè)過(guò)程中，不僅消除了模態(tài)波形的疊加，同時(shí)使波形輪廓更加對(duì)稱。

由于分解過(guò)程中多次應(yīng)用三次樣條插值法，造成信號(hào)兩端點(diǎn)附近出現(xiàn)失真，嚴(yán)重時(shí)可能使后續(xù)分解失去繼續(xù)研究的意義，所以在運(yùn)用FMD方法對(duì)非線性信號(hào)進(jìn)行分解時(shí)，必須進(jìn)行端點(diǎn)抑制［7－8］。

對(duì)于任意一個(gè)時(shí)間序列X(t)，都能得到它的希爾伯特變換結(jié)果Y(t)，［4］即:

(2)式中P為柯西主值。

X(t)與Y(t)構(gòu)成的復(fù)數(shù)信號(hào)［4］:

定義瞬時(shí)頻率 ω(t)［4］:

把(2)到(4)式所表示的變換用于每個(gè)固有模態(tài)函數(shù)序列，可以表示為［4］:

這里省略了殘余項(xiàng)rn，Re表示取實(shí)部。(5)式稱為Hilbert譜，記為［5］:

進(jìn)一步可得到Hilbert邊際譜為［5］:

其中T為數(shù)據(jù)序列長(zhǎng)度。

作為Hilbert邊際譜的附加結(jié)果，可以定義Hilbert瞬時(shí)能量譜如下［6］:

即瞬時(shí)能量譜表示了信號(hào)能量隨時(shí)間的分布情況。

此外，在對(duì)儲(chǔ)集層陣列聲波測(cè)井獲取的原始波列數(shù)據(jù)進(jìn)行EMD分解前要經(jīng)過(guò)預(yù)處理，包括環(huán)境校正、均衡化、去增益、去噪等，以保證波列數(shù)據(jù)的可靠、準(zhǔn)確。

3 某井油層數(shù)據(jù)的處理及分析

運(yùn)用Hilbert-Huang變換處理陣列聲波數(shù)據(jù)可以得到時(shí)間－頻率－能量譜、Hilbert邊際譜、Hilbert瞬時(shí)能量譜;邊際譜表示的是能量(振幅)在時(shí)間上的分布，瞬時(shí)能量譜是相同頻率能量在時(shí)間軸上的分布，能夠更直觀的體現(xiàn)出了能量和波至的關(guān)系［9］。

以某井2446.125米處的原始陣列波形記錄為例。常規(guī)測(cè)井對(duì)該點(diǎn)所處的地層解釋為油層。信號(hào)通過(guò)EMD分解，可以分解為數(shù)個(gè)有一定物理意義的IMF分量，每個(gè)IMF表征信號(hào)在某一特征尺度參數(shù)上的模態(tài)。油層原始波列信號(hào)的EMD分解示意圖如圖1所示，該信號(hào)自適應(yīng)地分解為7階固有模態(tài)函數(shù)(IMF)，每一個(gè)IMF分量都有不同的振幅和頻率，分解順序是按頻率從高至低進(jìn)行的。IMF1分量為信號(hào)的高頻，IMF2分量、IMF3分量、IMF4分量分別為信號(hào)的優(yōu)勢(shì)頻率子頻帶，體現(xiàn)了原始信號(hào)中顯著的信息。IMF5分量～I(xiàn)MF7分量為低頻分量，Res為信號(hào)本身微弱的趨勢(shì)或儀器的漂零。EMD分解實(shí)際上反映了該方法的自適應(yīng)濾波能力，獲得的IMF時(shí)序曲線大體上為單一諧波，是Hilbert變換進(jìn)行時(shí)頻估計(jì)的基礎(chǔ)。

圖1 油層陣列聲波測(cè)井的原始波列信號(hào)及其各IMF分量

基于EMD油層陣列聲波信號(hào)的Hilbert邊際譜如圖2所示，直觀地顯示了油層陣列聲波的大部分能量集中在5kHz以下。某一頻率出現(xiàn)的時(shí)刻可以在Hilbert時(shí)頻譜中找到?；贓MD油層陣列聲波信號(hào)的Hilbert幅度譜的三維示意圖，即時(shí)間－頻率－能量譜如圖3所示。顯示了縱波、橫波、斯通利波等不同振型波在不同時(shí)間、不同頻率出現(xiàn)的分布特征。基于EMD油層陣列聲波信號(hào)的Hilbert瞬時(shí)能量譜如圖4所示，顯示了波能分布是隨時(shí)間變化的。

圖2 基于EMD油層陣列聲波的Hilbert邊際譜

圖3 基于EMD油層陣列聲波的三維Hilbert時(shí)頻譜

綜合圖1－圖4，可知，油層測(cè)井信號(hào)的縱波初至?xí)r間為1000μs左右，主頻在8kHz左右。橫波時(shí)間在1500μs左右出現(xiàn)，主頻在7kHz左右。由于偽瑞利波的速度永遠(yuǎn)小于橫波速度，所以在橫波波至之后到達(dá)，和橫波疊加在一起。斯通利波的波至?xí)r間應(yīng)該在2500μs以后，頻率在0kHz～5kHz。能量譜出現(xiàn)的兩個(gè)峰值一個(gè)出現(xiàn)在低頻帶(0kHz～5kHz)上，一個(gè)出現(xiàn)在較高的頻帶(9kHz～13kHz)上。由此可見(jiàn)，低頻帶主要表示了斯通利波的能量，高頻帶表示了縱波、橫波和偽瑞利波能量的疊加。

圖4 基于EMD油層陣列聲波的Hilbert瞬時(shí)能量譜

干層測(cè)井信號(hào)的各種波能量都比較均勻，斯通利波能量最大，橫波分布較寬與斯通利波連在一起，橫波出現(xiàn)時(shí)間較早，多分成兩段。其明顯特征是偽瑞利波在橫波波至?xí)r間處能量較小，在縱波波至?xí)r間(3000μs)處較發(fā)育［10］。與干層測(cè)井信號(hào)各組分波分布特征相比，油層測(cè)井信號(hào)受到滲透率油水飽和度的影響，斯通利波能量衰減較嚴(yán)重。油層中縱波能量最小，橫波能量較大，一般大于斯通利波的能量。

4 結(jié)語(yǔ)

EMD依據(jù)信號(hào)本身的固有特性進(jìn)行分解，保證了陣列聲波測(cè)井信號(hào)分解后的非平穩(wěn)特性。分解后得到的IMF序列經(jīng)過(guò)Hilbert變換計(jì)算得到三維時(shí)頻譜圖，可以更好地表示出信號(hào)的頻率成分隨時(shí)間的變化規(guī)律，以及各組分波能量的衰減程度。

結(jié)果表明，與干層測(cè)井信號(hào)相比，油層測(cè)井信號(hào)由于受到滲透率油水飽和度的影響，斯通利波能量衰減較嚴(yán)重。油層中橫波能量較大，一般大于斯通利波的能量。

［1］ 王建華.聲波測(cè)井技術(shù)綜述［J］.工程地球物理學(xué)報(bào)，2006，3(5):395－400.

［2］ 李曉春.HHT在陣列聲波信號(hào)處理中的應(yīng)用［D］.長(zhǎng)春:吉林大學(xué)碩士學(xué)位論文，2009.

［3］ 付瑤等.EMD信號(hào)處理方法在LabVIEW和MATLAB中的實(shí)現(xiàn)［J］.北京機(jī)械工業(yè)學(xué)院學(xué)報(bào)，2008，23(2):23－27.

［4］ 王祝文，劉菁華，聶春燕.基于Hilbert-Huang變換的陣列聲波測(cè)井信號(hào)時(shí)頻分析［J］.地球科學(xué)－中國(guó)地質(zhì)大學(xué)學(xué)報(bào)，2008，33(3):387－392.

［5］ 鐘佑明.希爾伯特－黃變換瞬時(shí)信號(hào)分析理論的研究［D］.重慶:重慶大學(xué)博士學(xué)位論文，2002.

［6］ 郭喜平，王立東.經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)新算法及應(yīng)用［J］.內(nèi)蒙古科技大學(xué)學(xué)報(bào)，2008(5):70－72.

［7］M.G.A min，Spatialti me-frequency distributions for direction finding and blinds ource sep araiton［J］.in P roc.SP IE，1999，37(23):62－67.

［8］ 陳忠，鄭時(shí)雄.EMD信號(hào)分析方法邊緣效應(yīng)的分析［J］.數(shù)據(jù)采集與處理，2003，20(1):165－169.

［9］ 張學(xué)濤，王祝文，原鏡海.利用時(shí)頻分析方法在陣列聲波測(cè)井中區(qū)分油水層［J］.巖性油氣藏，2008，20(1):101－104.

［10］ 聶春燕，王祝文，杜丙錄.時(shí)頻分析在測(cè)井聲波信息提取中應(yīng)用［J］.長(zhǎng)春大學(xué)學(xué)報(bào)，2009，19(12):38－41.

Information Extraction of Acoustic Logging for Oil Layers Based on EMD Time-frequency Analysis

ZHANG Ting1，NIE Chun-yan2
(1.College of Electronics and Information Engineering，Changchun University of Science and Technology，Changchun 130022，China;2.College of Electronic Information Engineering，Changchun University，Changchun 130022，China)

Array acoustic logging signal is a typical non-linear and non-stationary signal.This paper extracts acoustic information reservoir properties in oil layers by adopting time-frequency analysis of experience mode decomposition(EMD).First of all，the signals are decomposed into a finite number of intrinsic modal functions(IMF)by EMD，and then 3-D distribution of Hilbert spectrum，Hilbert marginal spectrum and instantaneous energy spectrum are obtained by Hilbert transform.Simulation results show that the energy of P-wave is least and that of S-wave is stronger in oil layers，generally speaking，it is stronger than that of Stoneley Wave.The energy attenuation of Stoneley Wave is serious.

array acoustic;experience mode decomposition(EMD);time-frequency analysis;Hilbert transformation

P319.2

A

1009－3907(2011)08－0015－04

2011-05-22

張婷(1986-)，女，山東濟(jì)寧人，碩士研究生，主要從事信號(hào)檢測(cè)與信號(hào)處理方面研究。聶春艷(1966-)，女，吉林長(zhǎng)春人，教授，博士，主要從事信號(hào)檢測(cè)與處理、智能儀器等方面的研究。

責(zé)任編輯:吳旭云