田志

中國石油遼河油田分公司勘探開發(fā)研究院， 遼寧 盤錦 124010

0 引言

巖石的核磁共振(NMR)弛豫性質(zhì)分析在孔隙度、孔隙結(jié)構(gòu)、滲透率、潤濕性、流體飽和度及黏度等巖石物理參數(shù)評價方面發(fā)揮著重要作用(Coates et al., 1999；鄧克俊和謝然紅，2010；Wang et al., 2018；Liang et al., 2019；王迪等，2019；黨海龍等，2020；王香增等，2020).核磁共振弛豫信號是由流體的分子動力學(xué)和所處的物理化學(xué)環(huán)境共同決定.在連通的孔隙中，流體分子的布朗運(yùn)動會導(dǎo)致其所處的環(huán)境發(fā)生變化，這種變化會反映在核磁共振弛豫信號上.因此，只要通過一定的脈沖序列和量子相干，基于核磁共振技術(shù)就可以得到孔隙的連通性信息.目前，利用核磁共振評價孔隙連通性有兩種方法：一種是Song等(2000)提出通過探測孔隙內(nèi)部梯度磁場的分布(DDIF方法)得到孔隙尺寸分布，結(jié)合壓汞實(shí)驗(yàn)的孔喉分布，評價孔隙連通性.Zhang等(2018)利用DDIF方法結(jié)合CPMG測量弛豫信息，獲取多孔介質(zhì)的孔隙尺寸，評價樣品的非均質(zhì)性.但是，DDIF方法的應(yīng)用存在一定門檻.首先要求外部磁場較高并且非常均勻；其次，該方法是基于曲線形態(tài)的差異定性分析，無法定量評價孔隙連通性.第二種方法是通過探測孔隙間的擴(kuò)散耦合現(xiàn)象，得到孔隙中流體分子的弛豫交換速率，以此表征孔隙的連通性.主要思路分兩種：一種是分析擴(kuò)散耦合對一維T2譜形態(tài)的影響.一些學(xué)者通過建立含大孔和小孔(Ramakrishnan et al., 1999)或含裂縫和基質(zhì)孔(Chi and Heidari, 2015)的數(shù)字巖心模型，運(yùn)用數(shù)值模擬的方法得到巖石的核磁共振響應(yīng)特征，證實(shí)了不同尺度的連通孔隙間存在擴(kuò)散耦合現(xiàn)象.Anand和Hirasaki(2007)、Anand等(2008)通過實(shí)驗(yàn)分析發(fā)現(xiàn)溫度升高會加劇孔隙的擴(kuò)散耦合程度，巖石內(nèi)部黏土的賦存狀態(tài)不同會產(chǎn)生不同程度的擴(kuò)散耦合現(xiàn)象.這類方法的缺點(diǎn)是無法直接得到孔隙流體的弛豫交換速率，很難建立起與連通性的定量關(guān)系；第二種是基于二維核磁共振方法觀測擴(kuò)散耦合現(xiàn)象，定量計(jì)算弛豫交換速率.Lee等(1993)最早提出利用拉普拉斯變換方法反演得到二維的橫向弛豫時間譜，觀測到玻璃狀聚合物的弛豫交換現(xiàn)象.McDonald等(2005)首次設(shè)計(jì)出T2-T2脈沖序列(REXSY)，用該方法成功觀測到水泥的弛豫交換現(xiàn)象，實(shí)驗(yàn)證實(shí)只要流體分子在測量期時的孔隙尺寸環(huán)境發(fā)生變化，在二維譜上就會出現(xiàn)非對角峰.之后，一些學(xué)者陸續(xù)將其應(yīng)用到砂巖(Washburn and Callaghan, 2006)、碳酸鹽巖(Fleury and Soualem, 2009)、生物膜(Codd et al., 2011)等存在多尺度孔隙特征的天然樣品的擴(kuò)散耦合現(xiàn)象的測量分析中.Schwartz等(2013)、Johnson和Schwartz (2014)基于Ramakrishnan等(1999)建立的雙重孔隙數(shù)字巖心模型，運(yùn)用數(shù)值模擬的方法得到了巖石的T2-T2二維譜，結(jié)合雙孔弛豫交換模型的解析公式，分析弛豫交換速率與混合時間的相關(guān)性.Song等(2016)提出一種基于時域的數(shù)據(jù)處理方法分析孔隙擴(kuò)散耦合現(xiàn)象，避免了常規(guī)反演的不確定性.Yu等(2019)基于T2-T2脈沖序列提出一種直接表征多孔介質(zhì)平均孔徑的方法，能夠有效地提取不滿足快擴(kuò)散條件下的孔隙尺寸，有望提升碳酸鹽巖孔徑表征的準(zhǔn)確性.

針對存在多尺度孔隙的巖石，特別是碳酸鹽巖，受到沉積環(huán)境和成巖作用的影響，儲層具有強(qiáng)非均質(zhì)性，孔隙結(jié)構(gòu)復(fù)雜多樣，常常表現(xiàn)出微裂縫、粒間孔、粒內(nèi)孔和溶蝕孔洞同時發(fā)育的特點(diǎn)，孔隙的特征尺寸從幾十納米到幾厘米不等.在這種復(fù)雜的多尺度孔隙結(jié)構(gòu)背景下，流體分子在孔隙間的核磁共振弛豫交換會產(chǎn)生顯著的擴(kuò)散耦合現(xiàn)象.因此，充分挖掘巖石的核磁共振擴(kuò)散耦合現(xiàn)象與孔隙連通性之間的關(guān)系，對于深入認(rèn)識復(fù)雜孔隙結(jié)構(gòu)對滲透率、有效孔隙度等巖石物理參數(shù)的影響具有重要意義.目前，利用核磁共振評價孔隙連通性還處于早期探索階段，一些學(xué)者嘗試將非對角峰強(qiáng)度與混合時間ts建立關(guān)系，給出了耦合參數(shù)或者弛豫交換速率的定義.但是，關(guān)于非對角峰能量與混合時間ts的定量關(guān)系仍然存在爭議.本文的主要目的有兩點(diǎn)：(1)基于雙孔弛豫交換模型和多尺度孔隙數(shù)字巖石模型，分析T2-T2脈沖序列作用下非對角峰產(chǎn)生機(jī)制及影響因素；(2)計(jì)算表征孔隙間擴(kuò)散耦合強(qiáng)度的弛豫交換速率，分析擴(kuò)散耦合與孔隙連通性的相關(guān)性，為后續(xù)建立擴(kuò)散耦合參數(shù)與孔隙連通性及滲透率之間的定量關(guān)系作鋪墊工作.

1 核磁共振擴(kuò)散耦合探測原理

1.1 脈沖序列及其探測原理

利用T2-T2脈沖序列測量擴(kuò)散耦合現(xiàn)象.如圖1所示，序列由三部分組成，第一部分及第三部分為CPMG脈沖序列，第二部分為一個90°射頻脈沖.基于Torrey-Bloch方程(Bloch,1946；Torrey,1956)推導(dǎo)T2-T2脈沖序列的磁化矢量演化過程，說明測量原理.

圖1 T2-T2脈沖序列Fig.1 T2-T2 pulse sequence

(1)

其中,μ為體弛豫速率，φn為特征函數(shù)，τn為特征值.根據(jù)本征模的正交性，本征模系數(shù)an可表示為：

(2)

對整個孔隙空間進(jìn)行積分，得到總的磁化強(qiáng)度為：

(3)

其中，An為不同弛豫模態(tài)的相對強(qiáng)度.當(dāng)T2-T2脈沖序列測量時，磁化矢量先后經(jīng)歷編輯期t1、混合期ts及測量期t2等三個過程.假設(shè)編輯前流體分子已經(jīng)完全極化，在第一階段編輯期t1內(nèi)，孔隙流體經(jīng)歷橫向弛豫過程，磁化強(qiáng)度為：

(4)

第二階段為混合期ts, 該時期施加另一個90°射頻脈沖將磁化矢量搬轉(zhuǎn)到縱向平面，在混合時間ts內(nèi)，磁化矢量經(jīng)歷縱向弛豫T1過程.第三階段為測量期t2，該部分的CPMG序列將采集記錄回波信號，磁化強(qiáng)度的演化為：

(5)

其中：

(6)

使用量子力學(xué)中的bra-ket符號〈|〉將標(biāo)量積定義為：

(7)

(8)

圖2 T2-T2譜示意圖Fig.2 Sketch map of T2-T2 spectra

1.2 核磁共振響應(yīng)模擬方法

孔隙介質(zhì)中流體的橫向弛豫T2包括表面弛豫、擴(kuò)散弛豫和體弛豫三部分.總的磁化強(qiáng)度M(t)隨時間的變化可表示為(Toumelin et al., 2007; 鄒友龍等，2015；郭江峰等，2016)：

M(t)=M0[MB(t)Ms(t)MG(t)]，

(9)

式中M0為初始磁化強(qiáng)度，MB(t)、MS(t)、MG(t)分別代表t時刻與體弛豫、表面弛豫和擴(kuò)散弛豫相關(guān)的磁化矢量信號強(qiáng)度.

采用隨機(jī)游走(Random Walk)方法模擬大量自旋粒子在受限空間內(nèi)的布朗運(yùn)動求解擴(kuò)散方程(1).模型假設(shè)孔隙中飽和水，背景磁場為均勻場，回波間隔較小，可忽略擴(kuò)散弛豫項(xiàng)MG(t)的影響.由公式(9)可知，總的磁化矢量隨時間t的變化可通過模擬表面弛豫和體弛豫的信號衰減進(jìn)行表征.

初始時，隨機(jī)在孔隙空間內(nèi)分布大量游走粒子,假設(shè)所有粒子均完全極化.設(shè)定粒子的游走步長為ε，其對應(yīng)的時間步長為：

(10)

其中，D0為流體的自由擴(kuò)散系數(shù)，反映流體的布朗運(yùn)動.設(shè)粒子的位置坐標(biāo)為R=(xold,yold,zold)經(jīng)過時間步長Δt后，粒子的新位置為：

(11)

其中，φ的取值范圍為[0，π]，θ的取值范圍為[0，2π]，在范圍內(nèi)隨機(jī)取值，確保粒子運(yùn)動的隨機(jī)性.

粒子在擴(kuò)散過程中，會不斷地和孔隙與骨架的界面發(fā)生碰撞，磁化矢量會被界面部分吸收，該過程造成的衰減代表公式(9)中的表面弛豫MS(t)的變化，也指邊界條件.磁化矢量被吸收的強(qiáng)度與表面弛豫率ρ的大小有關(guān).通常有兩種模擬算法，一種是粒子以一定的概率被“殺死”，若沒被殺死，則粒徑反彈回原位置或鏡像對應(yīng)位置，反彈的概率可表示為：

(12)

在數(shù)據(jù)采集過程中，通過計(jì)算時刻t沒被“殺死”的粒子數(shù)與初始狀態(tài)時的粒子總數(shù)的比值來表征表面弛豫過程產(chǎn)生的磁化矢量信號的衰減.

另一種則不“殺死”，而是認(rèn)為只要粒子碰到邊界，磁化矢量就會衰減，衰減速率由與表面弛豫率ρ有關(guān)的指數(shù)函數(shù)給出，表面弛豫時間T2s可以由經(jīng)驗(yàn)公式給出：

(13)

體弛豫產(chǎn)生的信號衰減與孔隙結(jié)構(gòu)無關(guān)，只受流體本身性質(zhì)的影響，溫度一定時，水的體弛豫時間T2B是固定值，本文設(shè)定為3.0 s.模擬過程中，通過乘以一個與T2B相關(guān)的指數(shù)衰減項(xiàng)就可表征體弛豫MB(t)造成的信號衰減.

1.3 擴(kuò)散耦合模型建立

μGC模型如圖3所示，是具有雙峰孔徑分布特征的三維模型，模型大小設(shè)置為300 μm×300 μm×300 μm.采用實(shí)體球的空間堆積產(chǎn)生，實(shí)體的球形顆粒半徑為7.5 μm，圖中白色小圓點(diǎn)所示.通過均勻排列堆積實(shí)體的球形顆粒，形成一個立方體，這樣球形顆粒之間產(chǎn)生小尺寸孔隙，圖中黑色小圓點(diǎn)所示.為形成另外一種大尺寸孔隙，設(shè)置8組球心在立方體頂點(diǎn)，半徑為175 μm的大球，通過空間集合運(yùn)算，保留大球范圍內(nèi)的小球形顆粒，舍棄范圍外的小球形顆粒，得到立方體中心的大尺寸孔隙，圖中間黑色斑塊部分.建模主要參數(shù)設(shè)置如表1所示.

圖3 μGC模型及其切片(黑色為孔隙，白色為骨架)Fig.3 μGC model and its slice(black is pore. white is matrix)

表1 μGC模型參數(shù)設(shè)置Table 1 Parameter settings of μGC model

雙孔介質(zhì)弛豫交換模型是理論模型，只存在兩個弛豫模態(tài)，公式推導(dǎo)相對簡單，可以用解析解方式描述.模型假設(shè)滿足快擴(kuò)散條件，即認(rèn)為孔隙內(nèi)的磁化矢量總是均勻分布.小孔A的磁化矢量大小為SA，大孔B的磁化矢量大小為SB，初始時，總磁化矢量可表示為：

(14)

(15)

其中，KA和KB表示孔隙間的等效弛豫交換速率.

假設(shè)孔隙A和B的體積分別為VA和VB，當(dāng)自旋密度SA/VA=SB/VB時，孔隙A到B的弛豫交換量與孔隙B到A的弛豫交換量相等，故可定義弛豫交換速率為：

(16)

方程的特征值(Schwartz et al., 2013)可以表示為：

(17)

其中，uA,B=1/VA,B.在弱擴(kuò)散耦合機(jī)制下，即弛豫交換速率K趨近于0時：

(18)

在強(qiáng)擴(kuò)散耦合機(jī)制下：

+O(1/K)，

(19)

(20)

2 模擬結(jié)果分析與討論

表2 雙孔弛豫交換模型參數(shù)設(shè)置Table 2 Parameter settings of two-site model

2.1 不同混合時間下的模擬結(jié)果

圖4 基于μGC模型的六組混合時間不同的二維譜結(jié)果(a) 0 ms； (b) 100 ms； (c) 300 ms； (d) 600 ms； (e) 800 ms； (f) 1000 ms.Fig.4 T2-T2 maps of μGC model with different mixing time

圖5 基于雙孔弛豫交換模型的六組混合時間不同的二維譜結(jié)果(a) 0 ms； (b) 100 ms； (c) 300 ms； (d) 600 ms； (e) 800 ms； (f) 1000 ms.Fig.5 T2-T2 maps of two-site model with different mixing time

圖6 弛豫交換速率K的計(jì)算方法Fig.6 Calculation of relaxation exchange rate K

2.2 弛豫交換速率計(jì)算

2.3 數(shù)值解與解析解對比

圖7 μGC模型二維峰值強(qiáng)度與混合時間之間的關(guān)系Fig.7 The relationship between the 2D peak intensities and mixing time of μGC model

2.4 非對角峰影響因素分析

圖8 擴(kuò)散耦合通道關(guān)閉時的二維譜(a) μGC模型； (b) 雙孔弛豫交換模型.Fig.8 T2-T2 maps when the diffusional coupling channels closed(a) μGC model; (b) Two-site model.

進(jìn)一步分析μGC模型中非對角峰產(chǎn)生原因，保持上述實(shí)驗(yàn)條件不變，考察在T2-T2脈沖序列作用整個過程中磁化矢量的演化.分別模擬混合時間ts為0 ms、50 ms、300 ms的回波串信號，如果在混合期過后，磁化矢量發(fā)生額外衰減，則證實(shí)出現(xiàn)與擴(kuò)散耦合不相關(guān)的擴(kuò)散平均效應(yīng).圖9中長虛線代表混合時間為0的回波信號，實(shí)線代表混合時間不為0的回波信號.橫折線前為編輯期的回波信號，橫折線處為混合期，磁化強(qiáng)度不變.橫折線后為測量期，即T2-T2脈沖序列的數(shù)據(jù)采集期.將測量期的混合時間為0的回波信號平移，使其與混合時間不為0的回波信號在測量期開始處重合(圖中點(diǎn)劃線)，通過對比混合時間不為0的回波信號是否發(fā)生額外衰減，驗(yàn)證孔隙內(nèi)部是否產(chǎn)生擴(kuò)散平均效應(yīng).

圖9a是ts為50 ms與0 ms的對比圖，圖9b是300 ms與0 ms對比圖.在編輯期，孔隙流體發(fā)生正常的橫向弛豫過程，流體分子總的磁化強(qiáng)度不受混合時間長短影響，實(shí)線與長虛線重合，衰減速率一樣；在混合時期內(nèi)，由于孔隙間的通道關(guān)閉，不會發(fā)生擴(kuò)散耦合，流體分子只在孔隙內(nèi)部運(yùn)動，發(fā)生孔隙內(nèi)部的擴(kuò)散平均效應(yīng).由于不考慮縱向弛豫作用，所以該階段總的磁化強(qiáng)度保持不變；進(jìn)入測量期時，孔隙流體重新開始橫向弛豫過程，ts為300 ms的回波串信號的衰減明顯加快，而ts為50 ms的回波串信號衰減不明顯.這是因?yàn)?，相?0 ms,混合時間為300 ms時，自旋粒子充分?jǐn)U散，導(dǎo)致整體上大孔內(nèi)的自旋密度分布更為平均，擴(kuò)散平均效應(yīng)更顯著.因此，最終導(dǎo)致非對角峰的信號有一部分來自于該部分的貢獻(xiàn)，出現(xiàn)對角峰信號比理論值偏小，非對角峰偏大的現(xiàn)象.在孔隙尺度跨度大的巖石中，特別是碳酸鹽巖中，部分大孔隙處于慢擴(kuò)散的范疇，這種情況比較常見，在研究擴(kuò)散耦合現(xiàn)象時需要注意.

圖9 不同混合時間的回波串信號對比Fig.9 Comparison of echo train signals with different mixing time

3 碳酸鹽巖實(shí)例

采用基于沉積過程的數(shù)字巖石建模方法(田志等，2019)建立含多尺度孔隙特征的碳酸鹽巖數(shù)字巖石模型，探索不同尺度孔隙間的擴(kuò)散耦合現(xiàn)象，分析其對NMR橫向弛豫分布測量的影響，評價孔隙連通性.該建模方法通過模擬礦物顆粒在重力、內(nèi)摩擦力及浮力等作用下發(fā)生沉積作用，建立微觀尺度數(shù)字巖石模型.建模過程中，采取粒徑分布為雙峰特征的礦物顆粒模擬地層的沉積過程，粒徑范圍為40～110 μm，初始沉積完成后，進(jìn)行壓實(shí)和隨機(jī)溶蝕算法處理，處理后部分區(qū)域發(fā)育較大尺度的溶蝕孔，同時采取自仿射分形插值算法模擬裂縫的形成，在模型中加入部分微裂縫，最后對數(shù)據(jù)體離散化，形成碳酸鹽巖數(shù)字巖石模型(圖10)，圖中淺色代表固體骨架，黑色部分代表孔隙空間，包括微裂縫、粒間小孔、溶蝕大孔等三種類型孔隙，總孔隙度為20.8%.從圖11a的一維T2譜可以看出，孔徑分布具有多峰特征，微裂縫的孔隙度較小，粒間孔和溶蝕孔的孔隙度相當(dāng).

圖10 碳酸鹽巖數(shù)字巖石模型Fig.10 Digital rock model of carbonate rock

圖11 碳酸鹽巖數(shù)字巖心模型的四組混合時間不同的二維譜結(jié)果(a) 0 ms； (b) 100 ms； (c) 200 ms； (d) 300 ms. maps of carbonate rock model with different mixing time

圖12 碳酸鹽巖模型二維非對角峰值強(qiáng)度與混合時間之間的關(guān)系Fig.12 The relationship between the intensities of 2D non-diagonal peaks and mixing time of the carbonate rock model

4 結(jié)論

結(jié)果表明，T2-T2脈沖序列可以定量評價多尺度孔隙介質(zhì)的擴(kuò)散耦合現(xiàn)象.通過理論推導(dǎo)、數(shù)值模擬驗(yàn)證得到如下結(jié)論：

(1)針對雙尺度孔隙巖石模型，孔隙間的擴(kuò)散耦合強(qiáng)度與混合時間呈正相關(guān)性，本文提出的弛豫交換速率計(jì)算圖版評價擴(kuò)散耦合強(qiáng)度準(zhǔn)確可靠，具有一定的推廣應(yīng)用價值.

(3)針對孔隙尺寸跨度大的碳酸鹽巖模型，模擬結(jié)果表明微裂縫、小孔隙和大孔隙之間均存在擴(kuò)散耦合現(xiàn)象.隨混合時間的增加，代表不同類型孔隙的弛豫組分在T2-T2二維譜上信號能量的變化趨勢不同，同時T2譜的分布形態(tài)畸變程度加重，反映了孔隙間的連通性.

本文提出的T2-T2二維譜數(shù)據(jù)解釋方法對孔隙尺寸差異顯著的雙尺度孔隙巖石應(yīng)用效果較好，對孔徑連續(xù)分布或孔隙尺寸差異較小的巖石，定量計(jì)算弛豫交換速率存在一定困難.后續(xù)研究工作，將針對實(shí)際巖心設(shè)計(jì)實(shí)驗(yàn)，改進(jìn)數(shù)據(jù)的解釋方法，有效地提取擴(kuò)散耦合信息，進(jìn)一步探索弛豫交換強(qiáng)度與孔隙連通性及巖心滲透率之間的定量關(guān)系.

致謝感謝兩位匿名評審專家的寶貴意見！