EMD濾波方法在NMR測井中的應(yīng)用

程 誠，程晶晶，王光偉，薛志波

(1.華中科技大學(xué)控制科學(xué)與工程系，湖北武漢430074;2.中海油田服務(wù)股份有限公司，北京101149)

核磁共振(NMR)測井因能直接測量任意巖性儲層自由流體的滲流體積特性，在油田勘探中已得到廣泛應(yīng)用［1－2］。現(xiàn)代核磁共振測井采用交替相位CPMG脈沖序列測量NMR回波衰減信號，通過對其反演獲取地層流體橫向弛豫時間T2分布。這些分布數(shù)據(jù)包含地層流體物性與孔隙幾何形態(tài)信息，可定量評價束縛水飽和度、滲透率、流體類型等地層屬性［3］。T2譜反演過程需求解一組非適定性的線性方程組，而混雜在信號里的噪聲較易造成解的偏離。在原始信號信噪比較低時，反演結(jié)果偏離尤為嚴(yán)重。因此，對信號進行高效去噪是提高核磁共振測井解釋精度的關(guān)鍵。

經(jīng)驗?zāi)B(tài)分解［4－5］(empiricalmode decomposition，EMD)方法是HUANG等于1998提出的一種信號分析方法。它依據(jù)數(shù)據(jù)自身的時間尺度特征來進行信號分解，與建立在先驗性基函數(shù)的小波分解方法有本質(zhì)區(qū)別［6－7］，在處理非平穩(wěn)、非線性數(shù)據(jù)上具有非常明顯的優(yōu)勢。EMD將復(fù)雜信號分解為有限個本征模態(tài)函數(shù)(intrinsic mode function，IMF)，所分解出來的各個IMF分量包含原信號的不同時間尺度的局部特征信號，可以將其看作是時空濾波過程，利用這個性質(zhì)可以對信號進行濾波分析和降噪處理。

筆者將EMD算法應(yīng)用于核磁共振測井?dāng)?shù)據(jù)消噪處理中，改善原始回波數(shù)據(jù)信噪比，并在此基礎(chǔ)上進行T2譜反演，避免因信噪比低造成反演結(jié)果的偏離，為最終準(zhǔn)確獲取地層孔隙流體T2分布奠定基礎(chǔ)。

1 經(jīng)驗?zāi)B(tài)分解

經(jīng)驗?zāi)B(tài)分解，又叫Huang變換，是一種全新的處理非平穩(wěn)數(shù)據(jù)序列的方法，它將時間信號分解成一系列具有不同特征尺度的固有模態(tài)函數(shù)(IMF)，每個固有模態(tài)函數(shù)必須滿足以下兩個條件［8－9］:

(1)從全局特性上看，局部極值點與過零點的數(shù)目必須相等，或者最多相差一個;

(2)在任意時刻點，局部極大值點的包絡(luò)線(上包絡(luò)線)和局部極小值點的包絡(luò)線(下包絡(luò)線)的平均值必須為零。

EMD分解方法基于以下假設(shè)條件:①數(shù)據(jù)至少有兩個極值，一個極大值和一個極小值;②數(shù)據(jù)局部時域特性是由極值點間的時間尺度唯一確定;③如果數(shù)據(jù)沒有極值點但有拐點，則可以通過對數(shù)據(jù)微分一次或多次求得極值，再通過積分來獲得分解結(jié)果。

經(jīng)驗?zāi)B(tài)分解的一般步驟為:

(1)初始化ri(t)=x(t)，i=1。x(t)為要分析的時間序列;

(2)循環(huán)操作，獲取IMF分量。令h(t)=r1(t)，計算h(t)的極大值點和極小值點;將離散的極值點用三次樣條函數(shù)插值到整個時間段上，得到由所有局部極大值點構(gòu)成的上包絡(luò)線和所有局部極小值點構(gòu)成的下包絡(luò)線，分別記為u(t)和v(t)，并計算上下包絡(luò)的均值為:

計算信號與上、下包絡(luò)線的均值差為:

判斷h(t)是否趨向于零，若是，則滿足IMF的第二個條件，那么h(t)為第一個IMF，該循環(huán)終止，進入步驟(3);否則，以h(t)為輸入信號繼續(xù)循環(huán)，直至得到一個IMF，記為IMF1(t);

(3)令r1(t)=x(t)－IMF1(t);

(4)如果r1(t)仍然有至少2個極值點，回到步驟(2)繼續(xù)循環(huán)，直到所有的IMF分量均已被取出，分解過程結(jié)束，這時得到的余項rn(t)為趨勢分量，是一個單調(diào)信號或其值小于某個預(yù)先給定的閾值。

經(jīng)過上述分解過程后，得到n個IMF分量IMF1(t)，IMF2(t)，…，IMFn(t)及余項 rn(t)，于是原始信號可表示為:

EMD分解得到的IMF序列包含原信號的不同時間尺度的局部特征，是多通帶濾波結(jié)構(gòu)，與小波分解方法相比，這種方法是直觀的、直接的，經(jīng)EMD分解的各階本征模態(tài)函數(shù)能夠完全重構(gòu)，幾乎沒有能量損失。信號經(jīng)EMD分解后，噪聲主要集中在前幾個特定的IMF，對這幾個IMF進行處理，即直接利用時空濾波器進行消噪;對于混雜的信號和噪聲，直接利用閾值去噪的方法進行消噪，即可設(shè)置硬(或軟)門限去噪，筆者在此選擇軟門限去噪。

軟閾值方式將大于閾值的信號值與閾值作差，小于閾值的信號值直接置零，表達式如下:

式中:xi為原信號值;λcut為閾值參數(shù);x*i為閾值去噪后的信號值。

2 核磁共振信號EMD去噪處理

圖1 刻度筒原始回波串信號衰減曲線

如圖1所示，核磁共振測井的原始數(shù)據(jù)是幅度隨時間衰減的回波信號，其數(shù)學(xué)模型為［10］:

上述方程的求解往往是非適定的，為了求解該問題的數(shù)值解，通常先將積分方程離散化，式(5)可用下列離散方程組描述:

式中:m，n分別為測量過程中的回波個數(shù)和弛豫分量個數(shù);ti為采集時間(通常是回波間隔(TE)的整數(shù)倍);T2j為預(yù)先選擇的在對數(shù)尺度上等間隔分布的弛豫時間。

對于橫向弛豫信號，T2j為第j個弛豫分量的橫向弛豫時間，fj為第j個弛豫分量對總的橫向弛豫信號的貢獻。

信噪比是評價信號質(zhì)量的標(biāo)準(zhǔn)，對上述離散化回波模型信號，若將信號定義為x(t)，將對其進行最小二乘法擬合后的信號定義為x'(t)，則信號信噪比計算方法如下:

其中，norm(x(t))為x(t)的歐幾里得長度，信噪比的單位為dB。

2.1 實驗室刻度筒數(shù)據(jù)處理

圖1(a)為核磁共振測井儀在實驗室刻度筒中獲取的原始回波串?dāng)?shù)據(jù)，信號經(jīng)過paps對、相位旋轉(zhuǎn)［11］等處理后，得到的信號如圖1(b)所示。

如圖1(a)所示，原始回波信號噪聲很大，按照式(6)計算，信噪比為22 dB，經(jīng)過paps、相位旋轉(zhuǎn)等預(yù)處理之后，抑制了部分噪聲，信噪比提高到25 dB，但是信噪比仍不能滿足反演的要求，利用EMD方法對經(jīng)過預(yù)處理之后的信號進行消噪處理，如圖2所示，信號經(jīng)過EMD分解后自適應(yīng)地得到6階模態(tài)函數(shù)(如圖2(b)～圖2(g)所示)和殘差resigna l(如圖2(h)所示)。從圖2中可以看到每一個IMF分量都有不同的振幅和頻率，模態(tài)分量imf1～imf3包含了信號較高頻的成分，噪聲主要包含在這些分量中。對EMD分解后的imf1、imf2、imf3分別進行消噪處理，從圖2中可以看出，imf1可以作為噪聲直接置零濾去，imf2、imf3分別以0.5和0.3為閾值，按式(7)進行軟閾值去噪處理，最后重構(gòu)得到消噪后的回波串信號。

圖2 回波串信號EMD分解

圖3 (a)為核磁共振測井信號EMD消噪后的信號。與經(jīng)過預(yù)處理信號相比，噪聲幾乎完全得到抑制，如圖3(c)所示，噪聲道信號幅度基本保持在±100 mV以內(nèi)，其分布也接近符合正態(tài)分布的白噪聲，EMD濾波后信噪比由22 dB提高到46 dB。圖3(b)為經(jīng)過平滑濾波后的回波信號，其信噪比為35 dB，濾波效果沒有EMD理想。

圖3 去噪后的回波串信號

如圖4所示，把EMD去噪前后的回波串信號作在同一坐標(biāo)下，細線為原始信號，粗線為經(jīng)EMD去噪的信號，去噪后的信號基本位于原始信號的中間，表明EMD濾波在去噪的同時原始信號的特征得到了較好的保留。

圖4 核磁共振回波串信號去噪前后對比

2.2 現(xiàn)場測井?dāng)?shù)據(jù)處理

如圖5所示，該圖中的數(shù)據(jù)是標(biāo)準(zhǔn)油井的測井回波串信號，圖5(a)為原始信號經(jīng)過預(yù)處理后的回波串信號，圖5(b)為經(jīng)過EMD濾波處理之后的信號，圖5(c)為經(jīng)過平滑濾波處理之后的信號。相對于實驗室刻度筒數(shù)據(jù)，核磁共振現(xiàn)場測井?dāng)?shù)據(jù)信噪比更低，噪聲也更大，經(jīng)過預(yù)處理后信號信噪比只有11 dB，EMD濾波后信噪比提高到了31 dB，平滑濾波后信號信噪比只有17 dB。

圖5 標(biāo)準(zhǔn)井中EMD濾波前后回波串信號對比

3 EMD去噪后回波信號的反演處理

3.1 實驗室數(shù)據(jù)反演處理

如圖6所示，該圖是實驗室數(shù)據(jù)經(jīng)過EMD濾波處理之后進行反演的T2分布，經(jīng)EMD濾波的信號與原始回波信號反演結(jié)果存在較大的區(qū)別。原始回波信號信噪比低，只有22 dB，部分有用數(shù)據(jù)被刪除，某些奇異點反演后不穩(wěn)定，造成部分反演數(shù)據(jù)的丟失，因而無法顯示完整的T2譜信息。從圖6中可知，原始回波信號T2譜形狀畸變較大，而經(jīng)過EMD濾波處理之后，T2譜形狀畸變得到了明顯的改善，譜較光滑，對稱性良好，同時，從圖6中可以看出，T2譜的時間范圍為45～180 ms，符合自由水信號的弛豫時間范圍。

圖6 回波串信號EMD濾波反演結(jié)果對比圖

對數(shù)據(jù)進行平滑濾波，且在同樣的條件下進行T2譜反演，反演結(jié)果對比圖如圖7所示，經(jīng)過平滑濾波后得到的反演結(jié)果，與只經(jīng)過簡單的信號預(yù)處理后的T2譜接近，EMD濾波后高信噪比的信號反演結(jié)果更好。

圖7 反演結(jié)果對比圖

3.2 現(xiàn)場測井?dāng)?shù)據(jù)反演處理

對現(xiàn)場數(shù)據(jù)進行反演處理，該信號經(jīng)過疊加后進行EMD濾波處理，反演后T2分布如圖8所示，圖8中顯示的是雙峰，前一個譜峰時間為0.1～3 ms，為束縛水信號，后一個峰譜時間為9～20 ms，為自由水信號。

圖8 標(biāo)準(zhǔn)油井中EMD濾波后回波串信號T2譜

經(jīng)過EMD濾波處理后的數(shù)據(jù)具有良好的重復(fù)性，圖9為重復(fù)測量反演結(jié)果，對T2譜進行積分，即得孔隙度，所得值為624.2和621.9，該結(jié)果表明，經(jīng)過經(jīng)驗?zāi)B(tài)分解方法濾波后的現(xiàn)場信號信噪比得到明顯改善，孔隙度波動很小，基于該反演結(jié)果為后續(xù)解釋和評價提供了較好的基礎(chǔ)。

圖9 標(biāo)準(zhǔn)油井中重復(fù)測量反演T2譜對比圖

4 結(jié)論

受測井環(huán)境與工作機理限制，核磁共振測井原始信號信噪比低，反演解釋較為困難。筆者提出并實現(xiàn)將EMD濾波方法應(yīng)用到核磁共振測井回波串信號的消噪處理中。實驗表明，分解的模態(tài)函數(shù)能反映原始信號的固有特性，抑制信號的大部分噪聲，實驗室回波數(shù)據(jù)信噪比由22 dB提高到46 dB，標(biāo)準(zhǔn)油井回波數(shù)據(jù)信噪比由11 dB提高到31 dB，基于該濾波后的數(shù)據(jù)反演T2分布時間范圍更接近理論值，譜畸變得到明顯改善，多次重復(fù)測量中孔隙度結(jié)果基本一致。

［1］ GEORGE C，肖立志，MANFRED P.核磁共振測井原理與應(yīng)用［M］.北京:石油工業(yè)出版社，2007:2－75.

［2］王為民，李培，葉朝輝.核磁共振弛豫信號的多指數(shù)反演［J］.中國科學(xué)(A 輯)，2001，31(8):730 －736.

［3］HUANG N E.The empiricalmode decomposition and the Hilbert spectrum for nonlinear and non－stationary time series analysis［J］.Procedures of the Royal Society of London，1998，454(3):903 －995.

［4］HUANG N E，SHEN Z，LONG S R.A new view of nonlinear water waves:the Hilbert spectrum ［J］.Annual Review of Fluid Mechanics，1999(31):417 －457.

［5］ 孫延奎.小波分析及其應(yīng)用［M］.北京:機械工業(yè)出版社，2005:233－244.

［6］高靜，李朝偉，董云峰，等.空空導(dǎo)彈導(dǎo)引頭小波降噪?yún)?shù)優(yōu)選仿真研究［J］.航空兵器，2010，47(5):48－54.

［7］徐曉剛，徐冠雷，王孝通，等.經(jīng)驗?zāi)Ｊ椒纸?EMD)及其應(yīng)用［J］.電子學(xué)報，2009，37(3):581 －585.

［8］孫偉峰，彭玉華，許建華.基于EMD的激光超聲信號去噪方法［J］.山東大學(xué)學(xué)報:工學(xué)版，2008，38(5):121－126.

［9］原宏壯，陸大衛(wèi)，張辛耘，等.測井技術(shù)新進展綜述［J］.地球物理學(xué)進展，2005，20(3):786 －795.

［10］徐立強，何宗斌，郭書生.回波信號生成關(guān)鍵技術(shù)［J］.工程地球物理學(xué)報，2008，5(1):42 －47.

［11］王忠東，肖立志，劉堂宴.核磁共振弛豫信號多指數(shù)反演新方法及其應(yīng)用［J］.中國科學(xué)(G輯)，2003，33(4):323 －333.