多尺度孔隙巖石的核磁共振擴(kuò)散耦合現(xiàn)象及其探測(cè)方法

田志

中國(guó)石油遼河油田分公司勘探開(kāi)發(fā)研究院， 遼寧 盤(pán)錦 124010

0 引言

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

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

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

1.1 脈沖序列及其探測(cè)原理

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

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

(1)

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

(2)

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

(3)

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

(4)

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

(5)

其中：

(6)

使用量子力學(xué)中的bra-ket符號(hào)〈|〉將標(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)隨時(shí)間的變化可表示為(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時(shí)刻與體弛豫、表面弛豫和擴(kuò)散弛豫相關(guān)的磁化矢量信號(hào)強(qiáng)度.

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

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

(10)

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

(11)

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

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

(12)

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

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

(13)

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

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

μGC模型如圖3所示，是具有雙峰孔徑分布特征的三維模型，模型大小設(shè)置為300 μm×300 μm×300 μm.采用實(shí)體球的空間堆積產(chǎn)生，實(shí)體的球形顆粒半徑為7.5 μm，圖中白色小圓點(diǎn)所示.通過(guò)均勻排列堆積實(shí)體的球形顆粒，形成一個(gè)立方體，這樣球形顆粒之間產(chǎn)生小尺寸孔隙，圖中黑色小圓點(diǎn)所示.為形成另外一種大尺寸孔隙，設(shè)置8組球心在立方體頂點(diǎn)，半徑為175 μm的大球，通過(guò)空間集合運(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ì)弛豫交換模型是理論模型，只存在兩個(gè)弛豫模態(tài)，公式推導(dǎo)相對(duì)簡(jiǎn)單，可以用解析解方式描述.模型假設(shè)滿足快擴(kuò)散條件，即認(rèn)為孔隙內(nèi)的磁化矢量總是均勻分布.小孔A的磁化矢量大小為SA，大孔B的磁化矢量大小為SB，初始時(shí)，總磁化矢量可表示為：

(14)

(15)

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

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

(16)

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

(17)

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

(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 不同混合時(shí)間下的模擬結(jié)果

圖4 基于μGC模型的六組混合時(shí)間不同的二維譜結(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 基于雙孔弛豫交換模型的六組混合時(shí)間不同的二維譜結(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ù)值解與解析解對(duì)比

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

2.4 非對(duì)角峰影響因素分析

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

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

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

圖9 不同混合時(shí)間的回波串信號(hào)對(duì)比Fig.9 Comparison of echo train signals with different mixing time

3 碳酸鹽巖實(shí)例

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

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

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

圖12 碳酸鹽巖模型二維非對(duì)角峰值強(qiáng)度與混合時(shí)間之間的關(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脈沖序列可以定量評(píng)價(jià)多尺度孔隙介質(zhì)的擴(kuò)散耦合現(xiàn)象.通過(guò)理論推導(dǎo)、數(shù)值模擬驗(yàn)證得到如下結(jié)論：

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

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

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

致謝感謝兩位匿名評(píng)審專家的寶貴意見(jiàn)！