金國文, 謝然紅, 徐紅軍, 郭江峰, 高 倫

(1.中國石油大學(北京)油氣資源與探測國家重點實驗室,北京 102249; 2.中國石油大學(北京)地球探測與 信息技術(shù)北京市重點實驗室,北京 102249; 3.中國石油勘探開發(fā)研究院,北京 100083)

核磁共振(NMR)測井直接測量地層孔隙流體信息,在流體識別、孔隙度計算和束縛水評價等方面具有獨特優(yōu)勢[1-3]。NMR回波數(shù)據(jù)反演是NMR測井應用的基礎(chǔ),國內(nèi)外學者提出了多種反演方法[4-8]。截斷奇異值分解法[9-10]和Tikhonov正則化方法[11-13]是兩種常用的核磁共振反演方法,其中Tikhonov正則化方法又包括模平滑、斜率平滑和曲率平滑方法。但總體而言,現(xiàn)有的反演方法對于NMRT2譜短弛豫組分的反演精度還有待提高。致密砂巖、頁巖等非常規(guī)儲層孔徑小,在T2譜上反映為短弛豫組分居多,而反演結(jié)果的不準確則制約了核磁共振測井在非常規(guī)儲層中的應用。積分變換方法的發(fā)展與應用為提高NMRT2譜短弛豫組分的反演精度提供了有效途徑[14-17]。筆者將積分變換方法應用于核磁共振數(shù)據(jù)反演中,提出一種基于先驗信息約束的NMR反演新方法(ILT+),并通過數(shù)值模擬和實際井資料處理驗證新方法的有效性。

1 基于先驗信息約束的NMR反演新方法

在完全極化的條件下,NMR測井采集的CPMG回波串幅度可用如下第一類Fredholm積分方程來表示:

(1)

式中,b(t)為采集的NMR回波數(shù)據(jù);exp(-t/T2)為核函數(shù);T2為NMR橫向弛豫時間;f(T2)為T2譜幅度;ε(t)為高斯白噪聲。

公式(1)的離散矩陣形式為

Ap×qfq×1=bp×1.

(2)

其中

Aij=exp(-ti/T2,j);

fq×1=[f(T2,1),…,f(T2,q)]T;

bp×1=[b(t1)-ε(t1),…,b(tp)-ε(tp)]T;

i=1,…,p,j=1,…,q.

式中,p為回波個數(shù);q為T2布點數(shù)。

方程(1)的求解通常采用Laplace逆變換(ILT),但核磁共振數(shù)據(jù)ILT是一個嚴重病態(tài)問題,噪聲項導致解的擾動[18]。為了解決這一問題,通常采用正則化方法,求解如下目標函數(shù):

(3)

式中,W為權(quán)重矩陣;b為測量數(shù)據(jù)矩陣;L為系數(shù)矩陣;f為待反演T2譜的離散形式;α為正則化參數(shù)。

基于先驗信息約束的NMR回波數(shù)據(jù)反演(ILT+),通過積分變換方法直接從NMR原始回波數(shù)據(jù)中提取先驗信息作為約束條件來構(gòu)建新的反演目標函數(shù),具體為采用指數(shù)雙曲正弦變換(ESHT)求取錐形區(qū)域面積(S),采用梅林變換(MT)求取T2的ω次冪的均值(T),最后求解目標函數(shù)得到更為準確的NMRT2譜。具體步驟如圖1所示。

圖1 ILT+反演方法流程Fig.1 Workflow of ILT+

1.1 指數(shù)雙曲正弦變換求取錐形區(qū)域面積

傳統(tǒng)方法中,用T2截止值(單一截止值或漸進截止值)在反演得到的T2譜上劃分小孔隙和大孔隙,小于T2截止值的錐形區(qū)域面積對應束縛流體體積,大于T2截止值的錐形區(qū)域面積對應自由流體體積[19]。錐形區(qū)域面積可用如下公式計算:

(4)

式中,S為錐形區(qū)域面積;TC為T2截止值(可通過經(jīng)驗或?qū)嶒灚@得);K(T2,TC)為截止值函數(shù),由于本文中采用漸進截止值,可用漸進階躍函數(shù)近似。K(T2,TC)需要滿足如下特定值條件:

K(T2,TC)|T2=0=0,

K(T2,TC)|T2=TC=n,

式中,n、m為設定值。如圖2所示,實線為地層T2分布模型,虛線代表T2截止值為TC時對應的漸進階躍函數(shù),陰影部分為所求錐形區(qū)域面積。

假設K(T2,TC)對應的拉普拉斯逆變換函數(shù)k(t,TC)存在,即滿足關(guān)系式

(5)

則將公式(5)代入公式(4)中,結(jié)合公式(1)推導可得

(6)

式中,I{b(t)}表示對b(t)作積分變換。由公式(6)可知,若可求得核函數(shù)k(t,TC)滿足公式(5),則可直接通過對NMR原始數(shù)據(jù)作積分變換求取錐形區(qū)域面積。

圖2 地層T2分布模型和T2截止值 為TC時對應的漸進階躍函數(shù)Fig.2 T2distributions and tapered heaviside function withT2cutoff value equalingTC

公式(6)離散形式為

(7)

式中,tE為回波間隔;N為回波個數(shù)。進一步可得積分變換方法所求錐形區(qū)域面積對應的方差:

(8)

式中,σS為錐形區(qū)域面積標準差;σε為噪聲標準差。由公式(8)可知,只有當核函數(shù)k(t,TC)滿足平方可積,即能量E有限時:

(9)

所求錐形區(qū)域面積對應方差才是有限的,此時通過積分變換求取錐形區(qū)域面積的方法才可行。

本文中選取的核函數(shù)為指數(shù)雙曲正弦函數(shù),由于原始雙曲正弦函數(shù)不滿足公式(9)能量有限原則,故增加指數(shù)衰減項,如下式所示:

k(t,TC)=λe-βtsh(αt).

(10)

式中,sh(αt)為雙曲正弦函數(shù);λ、β、α為未知參數(shù);t為時間。

對公式(10)作拉普拉斯變換可得

(11)

根據(jù)K(T2,TC)滿足的特定值條件,計算得到各未知參數(shù)分別為

其中

通過數(shù)值模擬確定經(jīng)驗參數(shù)n=0.5,m=0.3。

由于所選核函數(shù)為指數(shù)雙曲正弦函數(shù),這里稱通過該積分變換求取錐形區(qū)域面積的方法為指數(shù)雙曲正弦變換法,簡記為ESHT。圖3顯示了TC=33 ms時對應的拉普拉斯變換函數(shù)對。

圖3 指數(shù)雙曲正弦變換在時域和T2域?qū)暮瘮?shù)形式Fig.3 Exponential hyperbolic sine transform in time domain andT2domain

1.2 梅林變換求取

(12)

(13)

(14)

1.3 反演目標函數(shù)構(gòu)建

將通過ESHT方法和MT方法從原始回波數(shù)據(jù)中提取的先驗信息作為約束條件,添加到原反演方法目標函數(shù)公式(3)中,構(gòu)建新反演方法的目標函數(shù):

(15)

(16)

式中,Na為所選TC個數(shù);Nm為所選ω個數(shù)。

(17)

式中,T2,min和T2,max分別為離散化的T2中的最小值和最大值。

(18)

2 數(shù)值模擬結(jié)果分析

模擬兩種具有雙峰結(jié)構(gòu)的T2分布地層模型,給定地層模型的孔隙度為10%,如圖4所示,T2布點數(shù)為64,其中圖4(a)為大孔占優(yōu)T2分布地層模型,圖4(b)為小孔占優(yōu)T2分布地層模型,縱坐標φ為孔隙度。根據(jù)此地層模型正演合成回波間隔為0.2 ms,回波個數(shù)為5 000的回波串,并加入噪聲水平為0.5 pu的高斯白噪聲,正演結(jié)果如圖5所示。

圖4 大孔占優(yōu)和小孔占優(yōu)T2分布模型Fig.4 T2distributions with more macropore and more micropore

圖5 大孔占優(yōu)和小孔占優(yōu)T2分布模型正演所得回波串Fig.5 Echo trains corresponding toT2distributions with more macropore and more micropore

為了提高NMR回波數(shù)據(jù)的反演速度,通常需要對回波數(shù)據(jù)進行壓縮,本文中采用截斷奇異值分解(truncated singular value decomposition)方法進行數(shù)據(jù)壓縮[20]。經(jīng)過壓縮后,圖5所示的每組回波均被壓縮至20個回波。

由圖6(a)和圖7(a)通過兩種模型反演結(jié)果進行定性分析可知,ILT方法和ILT+方法在T2譜長弛豫部分反演結(jié)果基本相同,但在短弛豫部分,ILT+方法相比ILT方法反演結(jié)果精度更高,尤其是對于小孔占優(yōu)模型,ILT+方法對反演結(jié)果精度的提高效果更明顯。ILT+方法反演效果受所加先驗信息種類及個數(shù)的影響,由圖7(a)可知:ILT+S方法主要提高T2譜短弛豫部分幅值精度,對譜峰聚焦性提高較小;ILT+T方法提高了T2譜幅值精度和譜峰聚焦性,但對T2譜幅值精度的改善效果弱于ILT+T方法;ILT+ST方法則綜合了ILT+S方法和ILT+T方法的優(yōu)點,同時提高了T2譜短弛豫組分幅值精度和譜峰聚焦性,對幅值精度改善效果優(yōu)于ILT+S,對譜峰聚焦性改善效果優(yōu)于ILT+T。

圖6(b)和圖7(b)給出了反演結(jié)果與模型之間的絕對誤差,其中黑色實線為ILT反演方法誤差,綠色虛線為ILT+S反演方法誤差,粉紅色虛線為ILT+T反演方法誤差,藍色虛線為ILT+ST反演方法誤差?？梢钥闯?ILT方法和ILT+方法在T2譜長弛豫部分反演結(jié)果誤差基本相同,但在短弛豫部分,ILT+方法相比ILT方法反演結(jié)果誤差更小,反演精度更高,并且ILT+ST方法反演結(jié)果誤差最小,ILT+S方法和ILT+T次之,ILT方法反演結(jié)果誤差最大。

圖6 大孔占優(yōu)模型反演結(jié)果及誤差對比Fig.6 Comparison of inversion results and errors of model with more macropore

圖7 小孔占優(yōu)模型反演結(jié)果及誤差對比Fig.7 Comparison of inversion results and error of model with more micropore

3 測井資料處理

數(shù)值模擬結(jié)果顯示ILT+方法相比于ILT方法反演結(jié)果精度有所提高,尤其是對于T2譜短弛豫部分,效果更明顯。非常規(guī)儲層孔徑較小,小孔部分占較大比重,ILT+方法對T2譜短弛豫部分反演精度的提高使得該方法在小孔占優(yōu)的非常規(guī)儲層中具有廣闊的應用前景。

本文中選取一段致密砂巖儲層NMR測井資料,其中NMR測井回波間隔為0.2 ms,回波個數(shù)為3 000,分別采用ILT、ILT+S、ILT+T和ILT+ST方法進行處理,處理結(jié)果如圖8所示。其中,第1道為深度曲線(DEPTH);第2道為井徑曲線(CAL)和自然伽馬曲線(GR);第3道為原始NMR回波數(shù)據(jù)(ECHO);第4道為ILT方法反演結(jié)果;第5道為ILT+S方法反演結(jié)果;第6道為ILT+T方法反演結(jié)果;第7道為ILT+ST方法反演結(jié)果;第8道為各方法計算的孔隙度(實線為ILT+ST孔隙度,虛線為ILT孔隙度,黑點為巖心孔隙度)。由第4道和第5道反演結(jié)果對比可知,ILT+S方法主要提高了T2譜短弛豫部分的幅值精度;由第4道和第6道反演結(jié)果對比可知,ILT+T方法主要提高了短弛豫部分譜峰的聚焦性;由第4、5和第6道反演結(jié)果對比可知,ILT方法在T2譜短弛豫部分反演精度不高,不僅幅值偏低,而且譜峰聚焦性差,甚至反演不出短弛豫組分,ILT+ST方法則綜合了ILT+S和ILT+T方法的優(yōu)點,同時提高了T2譜短弛豫部分的幅值精度和譜峰聚焦性;第8道孔隙度對比則驗證了ILT+ST方法的準確性,ILT方法由于對T2譜短弛豫部分反演精度不高,導致NMR孔隙度與巖心分析孔隙度相比普遍偏低,ILT+ST方法則由于對T2譜短弛豫部分反演精度的提高得到了更準確的孔隙度。實際井資料處理結(jié)果驗證了ILT+方法的有效性和實用性。

圖8 ILT和ILT+方法反演結(jié)果對比Fig.8 Comparison of inversion results of ILT and ILT+

4 結(jié) 論

(1)針對NMR測井中T2譜短弛豫組分反演不精確的問題,提出基于先驗信息約束的ILT+反演方法。該新方法提高了T2譜反演精度,不僅T2譜幅值更準確,而且譜峰更聚焦,尤其是對于T2譜短弛豫組分反演結(jié)果的改善效果更明顯。

(2)施加不同先驗信息作為約束條件對反演結(jié)果的改進效果不同。ILT+S方法主要提高T2譜短弛豫組分幅值精度,ILT+T方法主要提高T2譜短弛豫組分譜峰聚焦性,ILT+ST則綜合了上述兩種方法的優(yōu)點,同時提高了T2譜短弛豫組分幅值精度和譜峰聚焦性。

(3)新方法對于小孔占優(yōu)的致密砂巖、頁巖等非常規(guī)儲層的核磁共振回波數(shù)據(jù)反演具有廣闊的應用前景。