基于先驗(yàn)信息約束的核磁共振數(shù)據(jù)反演新方法

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

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

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

1 基于先驗(yàn)信息約束的NMR反演新方法

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

(1)

式中,b(t)為采集的NMR回波數(shù)據(jù);exp(-t/T2)為核函數(shù);T2為NMR橫向弛豫時(shí)間;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為回波個(gè)數(shù);q為T2布點(diǎn)數(shù)。

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

(3)

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

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

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

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

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

(4)

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

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

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

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

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

(5)

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

(6)

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

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

公式(6)離散形式為

(7)

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

(8)

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

(9)

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

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

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

(10)

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

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

(11)

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

其中

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

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

圖3 指數(shù)雙曲正弦變換在時(shí)域和T2域?qū)?yīng)的函數(shù)形式Fig.3 Exponential hyperbolic sine transform in time domain andT2domain

1.2 梅林變換求取

(12)

(13)

(14)

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

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

(15)

(16)

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

(17)

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

(18)

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

模擬兩種具有雙峰結(jié)構(gòu)的T2分布地層模型,給定地層模型的孔隙度為10%,如圖4所示,T2布點(diǎn)數(shù)為64,其中圖4(a)為大孔占優(yōu)T2分布地層模型,圖4(b)為小孔占優(yōu)T2分布地層模型,縱坐標(biāo)φ為孔隙度。根據(jù)此地層模型正演合成回波間隔為0.2 ms,回波個(gè)數(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ù)的反演速度,通常需要對(duì)回波數(shù)據(jù)進(jìn)行壓縮,本文中采用截?cái)嗥娈愔捣纸?truncated singular value decomposition)方法進(jìn)行數(shù)據(jù)壓縮[20]。經(jīng)過壓縮后,圖5所示的每組回波均被壓縮至20個(gè)回波。

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

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

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

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

3 測(cè)井資料處理

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

本文中選取一段致密砂巖儲(chǔ)層NMR測(cè)井資料,其中NMR測(cè)井回波間隔為0.2 ms,回波個(gè)數(shù)為3 000,分別采用ILT、ILT+S、ILT+T和ILT+ST方法進(jìn)行處理,處理結(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道為各方法計(jì)算的孔隙度(實(shí)線為ILT+ST孔隙度,虛線為ILT孔隙度,黑點(diǎn)為巖心孔隙度)。由第4道和第5道反演結(jié)果對(duì)比可知,ILT+S方法主要提高了T2譜短弛豫部分的幅值精度;由第4道和第6道反演結(jié)果對(duì)比可知,ILT+T方法主要提高了短弛豫部分譜峰的聚焦性;由第4、5和第6道反演結(jié)果對(duì)比可知,ILT方法在T2譜短弛豫部分反演精度不高,不僅幅值偏低,而且譜峰聚焦性差,甚至反演不出短弛豫組分,ILT+ST方法則綜合了ILT+S和ILT+T方法的優(yōu)點(diǎn),同時(shí)提高了T2譜短弛豫部分的幅值精度和譜峰聚焦性;第8道孔隙度對(duì)比則驗(yàn)證了ILT+ST方法的準(zhǔn)確性,ILT方法由于對(duì)T2譜短弛豫部分反演精度不高,導(dǎo)致NMR孔隙度與巖心分析孔隙度相比普遍偏低,ILT+ST方法則由于對(duì)T2譜短弛豫部分反演精度的提高得到了更準(zhǔn)確的孔隙度。實(shí)際井資料處理結(jié)果驗(yàn)證了ILT+方法的有效性和實(shí)用性。

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

4 結(jié) 論

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

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

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