王飛飛, 章海寧, 湯天知, 劉堂宴

(1.海洋地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室(同濟(jì)大學(xué)), 上海 200092; 2.中國(guó)石油集團(tuán)測(cè)井有限公司, 陜西 西安 710077)

0 引 言

核磁共振T2譜提供的巖石物性和流體特性的信息如巖石孔隙大小分布、有效孔隙度、滲透率、可動(dòng)流體與不可動(dòng)流體飽和度等為地質(zhì)勘探和油田開發(fā)過(guò)程中的儲(chǔ)層評(píng)價(jià)、產(chǎn)量預(yù)測(cè)和開發(fā)方法的制定提供了重要依據(jù)[1-2]。核磁共振測(cè)井原始回波串?dāng)?shù)據(jù)通過(guò)多指數(shù)反演獲得用于地層評(píng)價(jià)的T2譜。Butler等[3]討論了罰函數(shù)法(BRD),Dunn等[4]討論了非負(fù)最小二乘曲率平滑法,王為民等[5-6]討論了奇異值分解(SVD)法,翁愛(ài)華等[7]研究了核磁共振測(cè)井?dāng)?shù)據(jù)高分辨率反演方法,王忠東等[8]討論了聯(lián)合迭代重建算法(SIRT),陳華等[9]將差分進(jìn)化算法(DE)引入到核磁共振弛豫信號(hào)反演。其中,以SVD算法為一般應(yīng)用算法。本文考慮孔隙結(jié)構(gòu),在球管模型的基礎(chǔ)上,優(yōu)化常規(guī)的對(duì)數(shù)均勻布點(diǎn)方式,進(jìn)而優(yōu)化了核矩陣,并基于Tikhonov正則化思想,實(shí)現(xiàn)了SVD反演算法和共軛梯度CG算法的改進(jìn)。

1 核磁共振多指數(shù)弛豫模型

核磁共振擴(kuò)散理論表明,單個(gè)孔隙內(nèi)磁化強(qiáng)化信號(hào)的橫向弛豫具有單指數(shù)衰減規(guī)律。巖石內(nèi)部孔隙介質(zhì)是一系列大小不等的孔隙群組成,所測(cè)量的總磁化強(qiáng)度信號(hào)是一系列指數(shù)衰減信號(hào)的迭加,其計(jì)算模型為

d1=f1e-t1/T2,1+f2e-t1/T2,2+…+fne-t1/T2,n

d2=f1e-t2/T2,1+f2e-t2/T2,2+…+fne-t2/T2,n

? ? ? ?

dm=f1e-tm/T2,1+f2e-tm/T2,2+…+fne-tm/T2,n

(1)

式中,m,n分別為回波個(gè)數(shù)和弛豫分量個(gè)數(shù)(布點(diǎn)個(gè)數(shù)),m>n;ti(i=1,2,…,m)為第i個(gè)回波采集時(shí)間,是回波間隔(TE)的整數(shù)倍;fj(j=1,2,…,n)為第j個(gè)弛豫分量的初始幅度,經(jīng)過(guò)刻度,可以表示為第j種孔隙占總孔隙的百分比;di為第i個(gè)回波的幅度;T2j為橫向弛豫時(shí)間常數(shù)。矩陣形式

Ym×1=Am×nfn×1

(2)

通常,m>n,而且核矩陣A的條件數(shù)非常大,因此,該方程是個(gè)病態(tài)問(wèn)題。另外,物理意義上fj≥0即方程(2)還應(yīng)該滿足非負(fù)條件。

圖1 x和y與ε的關(guān)系曲線(a)及典型的L曲線圖(b)

1.1 Tikhonov截?cái)嗥娈愔捣?TTSVD)

Tikhonov正則化方法是解決病態(tài)問(wèn)題應(yīng)用最為普遍的方法[10-11],其估算準(zhǔn)則為

OF=‖Af-d‖2+ε‖Lf‖2=min

(3)

式中,ε為正則化因子,控制殘差和解的約束大小的權(quán)重;這里L(fēng)取單位矩陣,對(duì)應(yīng)零階正則化。該目標(biāo)函數(shù)兼有數(shù)據(jù)方差項(xiàng)和模型長(zhǎng)度項(xiàng)2項(xiàng)內(nèi)容,為了使得目標(biāo)函數(shù)達(dá)到最小,對(duì)f求偏導(dǎo)數(shù),并令其為0,得到正則化解

freg=(ATA+εI)-1ATd

(4)

大多數(shù)地球物理反演問(wèn)題,ε為0或足夠大都難以取得模型的最優(yōu)解,有必要尋找一種折中方案。其中,L曲線法是一種很好的確定正則化因子的方法。

L曲線法的基本思想是對(duì)于不同的正則化因子在雙對(duì)數(shù)坐標(biāo)下數(shù)據(jù)方差項(xiàng)和模型長(zhǎng)度項(xiàng)之間的關(guān)系曲線呈現(xiàn)出“L”形,其拐角處所對(duì)應(yīng)的正則化因子即是所要選擇的最佳正則化因子。

對(duì)于連續(xù)的L曲線,可以通過(guò)求取L曲線的曲率,曲率最大點(diǎn)所對(duì)應(yīng)的正則化因子即為所求。

(5)

顯然,x和y都是正則化因子的函數(shù)。通常x與ε呈現(xiàn)遞減關(guān)系,y與ε呈現(xiàn)遞增關(guān)系[見(jiàn)圖1(a)]。曲率k可以表示為

(6)

實(shí)際過(guò)程中,L曲線并不連續(xù)。一般做法是通過(guò)預(yù)先給定一系列正則化因子ε,可以得到對(duì)應(yīng)的x和y,在雙對(duì)數(shù)坐標(biāo)下繪出兩者的關(guān)系曲線。這種情況下可以構(gòu)造一系列的三角形,最大面積(也可以理解為L(zhǎng)曲線上的點(diǎn)到起始點(diǎn)連線段距離最大)所對(duì)應(yīng)的正則化因子即為最佳的正則化因子[見(jiàn)圖1(b)]。

不考慮正則化因子,則正則化解為最小二乘(LS)解

fLS=(ATA)-1ATd

(7)

對(duì)核矩陣A進(jìn)行SVD分解

(8)

式中,U和V滿足

UTU=VTV=In

(9)

則最小二乘解的奇異值分解形式為

(10)

奇異值λi減小得非常快,當(dāng)i值很大時(shí),對(duì)應(yīng)的奇異值λi很小,此時(shí)即是觀測(cè)數(shù)據(jù)中含有較小的誤差都將使得最小二乘解較大地偏離真實(shí)值。實(shí)際測(cè)量的回波數(shù)據(jù)中包含有噪音,為了避免上述現(xiàn)象,通過(guò)設(shè)置一個(gè)閾值,將奇異值中小于閾值的賦值為0,即可消除小的奇異值對(duì)最小二乘解得影響,達(dá)到消除方程病態(tài)性的目的。應(yīng)該注意到,去除的奇異值越多,損失的有用信息也越多。因此,有必要尋找一個(gè)合理的閾值,既能保證病態(tài)方程得意消除,又能保留足夠多的有用信息。

關(guān)于奇異值截?cái)?王為民等[5]提出利用信噪比進(jìn)行截?cái)嗟姆椒?/p>

(11)

林峰等[12]通過(guò)分析最佳奇異值保留個(gè)數(shù)與信噪比的關(guān)系,提出了改進(jìn)的奇異值截?cái)嗨惴?并且提出,在回波個(gè)數(shù)m=500、布點(diǎn)個(gè)數(shù)n=30的條件下可以取a=1.45,b=16,有

(12)

這里,信噪比SNR的定義方式為

(13)

非負(fù)約束方面,姜瑞忠等[13]從向量空間的角度對(duì)SVD算法進(jìn)行了分析,提出了一種新的實(shí)現(xiàn)非負(fù)約束的迭代方案。

由此,TTSVD反演算法可以概述為2步:L曲線法確定正則化因子,得到正則化解作為初始模型;TSVD迭代,加入非負(fù)約束條件。

1.2 共軛梯度法(CG)

CG法是由Hesteness和Stiefel于1952年為求解線性方程組而提出的,是解決大型線性和非線性方程組的最優(yōu)化有效的算法之一,其基本思想是把共軛性與最速下降法相結(jié)合,利用已知點(diǎn)處的梯度構(gòu)造一組共軛方向,并沿著此組方向進(jìn)行搜索,求出目標(biāo)函數(shù)的極小點(diǎn)。CG的算法流程如下。

(1) 初始化:初值x0∈Rn,計(jì)算殘差r0=b-Ax0和d0=ATr0,給定最大迭代次數(shù)和算法終止常數(shù)err>0;置k=0,開始迭代;

(2) 若‖rk‖

(4) 計(jì)算:βk=‖ATrk‖2/‖Ark-1‖2,dk=ATrk+βkdk-1;

k=k+1,回到(2)。

1.3 球管模型約束優(yōu)化反演

球管模型[14]考慮了球孔隙和管孔隙在不同組合形式下具有的核磁共振特征,其基本思路是根據(jù)預(yù)先的布點(diǎn)方式確定各孔隙組分的等效球半徑Re,通過(guò)預(yù)先設(shè)定的不同球管配比關(guān)系(引入“路徑”概念),加上表面積約束條件,可以確定新的布點(diǎn)方式,利用新的布點(diǎn)方式優(yōu)化核矩陣進(jìn)行反演,對(duì)比回波擬合曲線,找出擬合方差最小的一組路徑,即得到巖石最可能的理想孔隙結(jié)構(gòu)。

球管模型(見(jiàn)圖2)的描述方程

(14)

式中,S、V和R分別為面積、體積和半徑,下標(biāo)s、c為球孔隙和管孔隙。

圖2 球管模型平面示意圖

表面積約束條件

(15)

求解上述方程得到優(yōu)化的布點(diǎn)方式

(16)

理論上,預(yù)先設(shè)定的路徑有無(wú)數(shù)條,每一個(gè)布點(diǎn)處都可能對(duì)應(yīng)著不同大小的球管模型,即不同的Cd(管半徑與球半徑比值)值(見(jiàn)圖3)。為了方便計(jì)算,簡(jiǎn)化了路徑數(shù)量,給出了8種比較典型的路徑形式,用優(yōu)化的路徑處理實(shí)際數(shù)據(jù),速度大大提高。球管模型的思想可以概括為,將各巖石孔隙網(wǎng)絡(luò)等效為一系列不同大小的球孔隙和管孔隙,對(duì)應(yīng)不同Cd值的球管模型,掃描所有路徑并擇優(yōu),選擇相對(duì)最佳路徑和最優(yōu)解。

圖3 球管模型路徑掃描

圖4 正演雙峰模型與回波衰減曲線

2 正演數(shù)值模擬

構(gòu)造一個(gè)典型的雙峰譜模型,布點(diǎn)時(shí)間0.3～3 000 ms,布點(diǎn)個(gè)數(shù)為30,回波時(shí)間為0.9 ms,回波個(gè)數(shù)500[見(jiàn)圖4(a)]。對(duì)正演的回波[見(jiàn)圖4(b)]加入不同比例的噪音,信噪比定義為第1個(gè)回波幅度與噪音(模型回波衰減曲線與回波擬合曲線之差)方差之比,即

(17)

正演結(jié)果表明,不同程度的噪音影響,CG和TTSVD的2種方法都能反演出模型譜的雙峰結(jié)構(gòu)(見(jiàn)圖5)。噪音水平為5%,2種方法反演譜與模型譜差別不顯著;噪聲水平低于1%,2種方法反演的譜與模型譜幾乎完全重合?？傊?2種方法都能提供較好的反演結(jié)果。TTSVD的反演結(jié)果略優(yōu)于CG法,但反演所需的時(shí)間比CG多。

3 實(shí)際巖心檢驗(yàn)與核磁共振測(cè)井資料處理

圖5 CG與TTSVD正演結(jié)果圖

用這種優(yōu)化核矩陣的方法反演了不同地區(qū)共45塊巖心樣品。其中,除個(gè)別巖樣屬中孔隙度(介于15%～25%)外,其余都屬于低孔隙度(介于10%～15%),大部分屬于特低孔隙度(小于10%)。

圖6 巖心反演結(jié)果對(duì)比

從圖6可知,實(shí)驗(yàn)室給的反演結(jié)果大部分都在1個(gè)孔隙度誤差范圍內(nèi),改進(jìn)的奇異值分解方法(TTSVD)反演結(jié)果都在1個(gè)孔隙度誤差范圍內(nèi),共軛梯度法(CG)反演結(jié)果除個(gè)別點(diǎn)結(jié)果都在1個(gè)孔隙度誤差范圍內(nèi)。2種方法反演結(jié)果都優(yōu)于實(shí)驗(yàn)室給的結(jié)果。圖7是相對(duì)誤差情況,實(shí)驗(yàn)室結(jié)果相對(duì)誤差大部分都在5%以上,平均為7.43%,2種方法的相對(duì)誤差大部分都在5%以內(nèi),平均TTSVD為3.93%、CG為3.54%。特別如黑框所示,巖樣屬特低孔隙度,相對(duì)誤差都在10%以上,最高可達(dá)18%;該方法尤其是CG法基本都在10%以內(nèi)。從T2譜形態(tài)上看,對(duì)于單峰形態(tài),2種方法和實(shí)驗(yàn)室結(jié)果幾乎一致;對(duì)于實(shí)驗(yàn)室給出的不明顯的雙峰形態(tài),2種方法都能給出明顯的雙峰形態(tài)(見(jiàn)圖8)。

圖7 巖心反演誤差分析

圖8 2塊巖樣的反演T2譜對(duì)比

將2種優(yōu)化反演方法用C語(yǔ)言編制成相應(yīng)的計(jì)算機(jī)程序,依托CIFlog測(cè)井解釋平臺(tái),實(shí)現(xiàn)對(duì)實(shí)際核磁共振測(cè)井?dāng)?shù)據(jù)的反演(見(jiàn)圖9),2種方法反演的結(jié)果基本類似,無(wú)明顯差別,僅在T2譜形態(tài)和幅度上略微有不同。

圖9 實(shí)際核磁共振測(cè)井?dāng)?shù)據(jù)反演實(shí)例

4 結(jié) 論

(1) 考慮Tikhonov正則化,提供初始模型;使用截?cái)嗥娈愔捣纸膺M(jìn)行非負(fù)約束迭代,得到最終解,正演模型效果很好。

(2) 共軛梯度算法在T2譜反演中有很好的應(yīng)用。

(3) 以球管模型為基礎(chǔ),優(yōu)化T2弛豫時(shí)間布點(diǎn)方式,改進(jìn)核矩陣,進(jìn)而優(yōu)化T2譜反演,對(duì)比本文方法反演譜和實(shí)驗(yàn)室反演譜有很大程度上的相似性,并且對(duì)T2譜的多峰形態(tài)有更好的反映?？紫抖扔?jì)算結(jié)果,該方法結(jié)果與氣測(cè)孔隙度更加接近,說(shuō)明該算法具有很好的實(shí)用性,可以用于巖心核磁共振T2譜反演,尤其是對(duì)于低孔隙度甚至特低孔隙度巖心有很好的應(yīng)用效果。

(4) 依托CIFlog測(cè)井解釋平臺(tái),實(shí)現(xiàn)了實(shí)際核磁共振測(cè)井?dāng)?shù)據(jù)的反演。

(5) 考慮孔隙結(jié)構(gòu)的核磁共振T2譜反演方法可以作為計(jì)算孔隙結(jié)構(gòu)參數(shù)如分選系數(shù)、歪度、峰度、中值半徑等的基礎(chǔ)。

參考文獻(xiàn):

[1] Kenyon W E. Petrophysical Principles of Applications of NMR Logging [J]. Log Analyst, 1997, 38: 21-43.

[2] Timur A. Producible Porosity and Permeability of Sandstones Investigated Through Nuclear Magnetic Resonance Principles [J]. The Log Analyst, 1969, 10(1): 3-11.

[3] Butler J P, Reeds J A, Dawson S V. Estimating Solutions of First Kind Integral Equations with Nonnegative Constraints and Optimal Smoothing [J]. SIAM Journal on Numerical Analysis, 1981, 18(3): 381-397.

[4] Dunn K J, LaTorraca G A, Warner J L, et al. On the Calculation and Interpretation of NMR Relaxation Time Distributions [J]. SPE Paper, 1994, 28369.

[5] 王為民, 李培. 核磁共振弛豫信號(hào)的多指數(shù)以演 [J]. 中國(guó)科學(xué): A輯, 2001, 31(8): 730-736.

[6] 李鵬舉, 葛成, 孫國(guó)平, 等. 基于奇異值分解的核磁共振測(cè)井T2譜反演方法的改進(jìn) [J]. 測(cè)井技術(shù), 2010, 34(3): 215-218.

[7] 翁愛(ài)華, 李舟波, 莫修文, 等. 核磁共振測(cè)井?dāng)?shù)據(jù)高分辨率反演方法研究 [J]. 測(cè)井技術(shù), 2002, 26(6): 455-459.

[8] 王忠東, 肖立志, 劉堂宴. 核磁共振弛豫信號(hào)多指數(shù)反演新方法及其應(yīng)用 [J]. 中國(guó)科學(xué): G輯, 2003, 33(4): 323-332.

[9] 陳華, 潘克家, 譚永基. 核磁共振弛豫信號(hào)多指數(shù)反演新方法 [J]. 測(cè)井技術(shù), 2009, 33(1): 37-41.

[10] Hansen P C, O’Leary D P. The Use of the L-curve in the Regularization of Discrete Ill-posed Problems [J]. SIAM Journal on Scientific Computing, 1993, 14(6): 1487-1503.

[11] Hansen P C. Analysis of Discrete Ill-posed Problems by Means of the L-curve [J]. SIAM review, 1992, 34(4): 561-580.

[12] 林峰, 王祝文, 劉菁華, 等. 核磁共振T2譜奇異值反演改進(jìn)算法 [J]. 吉林大學(xué)學(xué)報(bào): 地球科學(xué)版, 2009, 39(6): 1150-1155.

[13] 姜瑞忠, 姚彥平, 苗盛, 等. 核磁共振T2譜奇異值分解反演改進(jìn)算法 [J]. 石油學(xué)報(bào), 2006, 26(6): 57-59.

[14] 劉堂晏, 周燦燦, 馬在田, 等. 球管孔隙模型的約束尋優(yōu)及應(yīng)用 [J]. 同濟(jì)大學(xué)學(xué)報(bào): 自然科學(xué)版, 2006, 34(11): 1464-1469.