碳酸鹽巖孔隙-裂縫雙重網(wǎng)絡(luò)模型構(gòu)建及分析

羅瑜，王寅，王容，袁雯

（1.中國石油 西南油氣田分公司 勘探開發(fā)研究院，成都 610093；2.科吉思石油技術(shù)咨詢有限公司，北京 100176）

對于裂縫和孔隙均發(fā)育的碳酸鹽巖滲流特征的研究，通常是構(gòu)建裂縫和具有一般孔隙介質(zhì)的巖塊系統(tǒng)，利用有限元法耦合流體在孔隙和裂縫中的流動(dòng)特征[1-3]。此種方法運(yùn)算過程較為復(fù)雜，且對于非均質(zhì)性強(qiáng)的樣品所構(gòu)建的模型體量較大，對計(jì)算設(shè)備要求較高，每塊樣品的研究周期較長，無法被廣泛應(yīng)用于實(shí)際生產(chǎn)中。

孔隙網(wǎng)絡(luò)模型自引入巖石分析領(lǐng)域以來，已成為表征巖石微觀孔隙結(jié)構(gòu)以及預(yù)測微觀滲流特征的重要工具。以往的孔隙網(wǎng)絡(luò)模型主要用于多孔介質(zhì)的研究，少有針對裂縫介質(zhì)的研究。本文將孔隙-裂縫型碳酸鹽巖的孔隙空間和裂縫空間分離，利用最大球填充法[4]提取孔隙網(wǎng)絡(luò)模型，中心線提取法與最大球填充法相結(jié)合提取裂縫網(wǎng)絡(luò)模型，再將孔隙網(wǎng)絡(luò)模型與裂縫網(wǎng)絡(luò)模型融合[5]，構(gòu)建了孔隙-裂縫雙重網(wǎng)絡(luò)模型。最后，基于構(gòu)建的模型，分別計(jì)算了3 塊孔隙-裂縫型碳酸鹽巖樣品的孔隙度和滲透率，與常規(guī)物性實(shí)驗(yàn)所得數(shù)據(jù)進(jìn)行對比，對模型的可靠性進(jìn)行了驗(yàn)證。

1 模型構(gòu)建

1.1 數(shù)字巖心模型構(gòu)建

通常利用CT 掃描獲得孔隙-裂縫型碳酸鹽巖樣品的灰度圖像，該方法可以在不損壞樣品的情況下，觀察樣品的內(nèi)部微觀結(jié)構(gòu)[6-7]。利用CT掃描獲得的灰度圖像建立數(shù)字巖心模型的步驟如下：①利用CT掃描設(shè)備對經(jīng)過鉆取、切割等初步處理后的碳酸鹽巖樣品進(jìn)行掃描，獲得其投影數(shù)據(jù)；②利用圖像重建方法，將獲取的投影數(shù)據(jù)重建成灰度圖像（圖1a）；③對灰度圖像進(jìn)行預(yù)處理后，利用圖像分割技術(shù)，將圖像中的孔隙和裂縫與巖石骨架分離，建立數(shù)字巖心模型（圖1b）。

圖1 基于CT掃描灰度圖像分割得到的孔隙和裂縫Fig.1.Pores and fractures from segmented CT images

根據(jù)孔隙和裂縫不同的幾何形狀，計(jì)算費(fèi)雷特長寬比（單個(gè)孔隙或裂縫的最長直徑與最短直徑之比）、比表面積等參數(shù)，選取合適的閾值，有效分離出孔隙和裂縫（圖1c，圖1d）。

1.2 孔隙網(wǎng)絡(luò)模型構(gòu)建

在分離出孔隙-裂縫型碳酸鹽巖樣品的孔隙空間后，利用最大球填充法，對孔隙空間的網(wǎng)絡(luò)模型進(jìn)行提取[4]。該方法以分割出的孔隙空間中任一體素為基本點(diǎn)，采用擴(kuò)張算法確定可能的搜索空間，確定出距離每個(gè)體素最近的骨架體素范圍，采用收縮算法識(shí)別出中心體素對應(yīng)的內(nèi)切球，計(jì)算出該內(nèi)切球半徑的上限和下限。

在提取出一系列以體素為球心的最大球后，由于部分球是其他球的子集，需要對提取的最大球進(jìn)行篩選。篩選方法如下：①根據(jù)每個(gè)體素對應(yīng)的最大球半徑的平方大小，以鏈表的形式表示體素；②當(dāng)循環(huán)計(jì)算至相關(guān)系數(shù)為0 時(shí)，停止運(yùn)算；③對一個(gè)球心為O1、半徑為r1的最大球，在生成的鏈表中搜索另一個(gè)球心為O2的最大球，這2個(gè)最大球球心之間的距離為d，是否滿足條件，以此判斷是否和指定的球存在包含關(guān)系；④將存在包含關(guān)系的最大球從鏈表中刪除，在鏈表中留下的就是沒有被其他球所包含的最大球集合。

采用成簇算法對最大球集合中的元素進(jìn)行孔隙和喉道的劃分[4]（圖2）。經(jīng)過成簇算法的處理，代表孔隙的最大球集合被轉(zhuǎn)換成相互連接的簇，每一個(gè)簇中的主球代表孔隙空間中的一個(gè)孔隙體，不同簇之間相連的節(jié)點(diǎn)代表該孔隙空間中的喉道。當(dāng)確定一個(gè)喉道之后，即可得到該喉道分別通往兩端主球的鏈路，而鏈路則是由喉道中心與通往主球之間所有最大球球心的鏈接組成，這些鏈路構(gòu)成了整個(gè)孔隙網(wǎng)絡(luò)模型的骨架。

獲取孔隙網(wǎng)絡(luò)模型骨架后，需要對每個(gè)孔隙和喉道的幾何特征進(jìn)行表征，包括喉道長度、孔隙半徑（最大球半徑）、喉道半徑、形狀因子等。

喉道半徑計(jì)算方法如下：首先，在每個(gè)喉道中軸線體素的局部范圍內(nèi)，過該體素作垂直于中軸線的垂面，垂面范圍不超過喉道范圍，任意一個(gè)垂面都與喉道巖石壁面相交得到圖中所示的喉道橫截面。在橫截面內(nèi)以喉道中軸線上的體素為中心，采用最大球填充法中的球體膨脹方法得到該橫截面內(nèi)的喉道最大內(nèi)切圓半徑（圖3）。依次計(jì)算出所有橫截面的最大內(nèi)切圓半徑，最終選取所有最大內(nèi)切圓半徑的最小值，即為喉道半徑。

圖2 成簇算法劃分孔隙及喉道Fig.2.Pores and throats from clustering maximum balls

圖3 喉道半徑計(jì)算示意Fig.3.Schematic diagram of throat radius

形狀因子可補(bǔ)充具有不規(guī)則形狀的真實(shí)孔隙空間在提取孔隙網(wǎng)絡(luò)模型時(shí)被簡化成節(jié)點(diǎn)-鏈接形式后的形狀信息[8]，用公式表示為

1.3 裂縫網(wǎng)絡(luò)模型構(gòu)建

由于裂縫空間可以近似看作由大量連通孔隙組成的空間，因此需要盡量避免孔隙網(wǎng)絡(luò)提取對裂縫空間簡化造成的信息丟失，本文采用中心線提取法與最大球填充法相結(jié)合提取裂縫網(wǎng)絡(luò)模型。

依據(jù)距離圖法[9]，結(jié)合細(xì)化算法[10]，進(jìn)行裂縫中心線的提?。菏紫龋?jì)算分離出的裂縫空間數(shù)字巖心模型中的每一個(gè)體素到最近巖石骨架的距離；然后，根據(jù)距離圖法，不斷細(xì)化裂縫空間的數(shù)字巖心模型，直至形成一條連續(xù)的體素鏈，該體素鏈即為提取出的裂縫中心線。

基于提取的裂縫的中心線體素坐標(biāo)，采用最大球填充法，以每一個(gè)中心線體素為基本點(diǎn)，填充裂縫空間，從而得到碳酸鹽巖樣品的裂縫網(wǎng)絡(luò)模型。該模型將裂縫空間近似等效為密集的連通孔隙空間，用一系列節(jié)點(diǎn)鏈接來表征裂縫空間的結(jié)構(gòu)，并計(jì)算裂縫開度、裂縫長度、裂縫形狀等參數(shù)，進(jìn)而分析其滲流特征。

1.4 孔隙-裂縫雙重網(wǎng)絡(luò)模型構(gòu)建

依據(jù)上述方法分別獲得孔隙網(wǎng)絡(luò)模型和裂縫網(wǎng)絡(luò)模型，在表征孔隙和裂縫的最大球之間添加連通喉道，實(shí)現(xiàn)孔隙-裂縫網(wǎng)絡(luò)模型構(gòu)建，方法如下：①對孔隙與裂縫的接觸面進(jìn)行標(biāo)記，并找到該接觸面對應(yīng)的孔隙和裂縫的最大球編號；②根據(jù)接觸面的數(shù)量、面積以及對應(yīng)的孔隙最大球編號，建立孔隙與裂縫之間連通的喉道半徑及配位數(shù)分布關(guān)系；③根據(jù)喉道半徑及配位數(shù)分布關(guān)系，對接觸面對應(yīng)的孔隙和裂縫的最大球，選取合適半徑和數(shù)量的喉道將兩者連接，從而構(gòu)建出孔隙-裂縫雙重網(wǎng)絡(luò)模型（圖4）。

圖4 孔隙-裂縫型碳酸鹽巖的孔隙-裂縫雙重網(wǎng)絡(luò)模型Fig.4.Pore?fracture network model of dual?porosity carbonate rocks

2 模型巖石物性參數(shù)計(jì)算

孔隙-裂縫雙重網(wǎng)絡(luò)模型的孔隙度可以通過計(jì)算孔隙空間體積和裂縫空間體積之和與樣品總體積的百分比得到[11]：

孔隙-裂縫雙重網(wǎng)絡(luò)模型為擬靜態(tài)模型，即所有流體流動(dòng)受毛細(xì)管力控制[12]，由黏滯力造成的壓降忽略不計(jì)。計(jì)算絕對滲透率時(shí)，假設(shè)喉道中飽和一種流體，在相連的孔隙i與孔隙j之間的喉道兩端分別施加驅(qū)動(dòng)壓力pi和pj，流量與壓降的關(guān)系為[13-15]

根據(jù)泊肅葉定律，傳導(dǎo)率gij為

對于不可壓縮流體的穩(wěn)態(tài)流，模型中進(jìn)入每個(gè)孔隙中的流量應(yīng)與流出的流量相等，因此每個(gè)孔隙的質(zhì)量守恒可表示為

將（3）式—（5）式聯(lián)立求解，可得到每一個(gè)孔隙的壓力，以此應(yīng)用于所有孔隙，進(jìn)一步可以求出模型兩端壓差下的總流量Q。根據(jù)達(dá)西公式，可計(jì)算出孔隙-裂縫雙重網(wǎng)絡(luò)模型的絕對滲透率：

3 模型驗(yàn)證

為檢驗(yàn)孔隙-裂縫雙重網(wǎng)絡(luò)模型構(gòu)建方法的可靠性，選取四川盆地雙魚石區(qū)塊的孔隙-裂縫型碳酸鹽巖樣品3 塊，將基于孔隙-裂縫雙重網(wǎng)絡(luò)模型計(jì)算得出的孔隙度和絕對滲透率，與常規(guī)物性實(shí)驗(yàn)得出的孔隙度和絕對滲透率進(jìn)行對比。

對3塊樣品分別進(jìn)行CT掃描，掃描精度為7.832 μm，視域范圍為13 mm×13 mm×10 mm，通過CT 掃描灰度圖像，可以看到樣品內(nèi)部孔隙和裂縫均發(fā)育（圖5）。

對CT 掃描圖像進(jìn)行初步分割，將孔隙空間和裂縫空間與巖石骨架進(jìn)行初步劃分。再根據(jù)幾何形態(tài)差異，進(jìn)一步劃分孔隙和裂縫，分別得到每塊樣品的孔隙和裂縫數(shù)字巖心模型，并利用本文方法分別提取孔隙網(wǎng)絡(luò)模型和裂縫網(wǎng)絡(luò)模型（表1）。

圖5 四川盆地雙魚石區(qū)塊孔隙-裂縫型碳酸鹽巖樣品CT掃描圖像Fig.5.Micro?CT images of dual?porosity carbonate rocks

分別計(jì)算3 塊樣品的裂縫數(shù)字巖心模型和裂縫網(wǎng)絡(luò)模型的裂縫體積，樣品A、樣品B 和樣品C 的裂縫數(shù)字巖心模型裂縫體積分別為1.08×109μm3、1.13×109μm3和4.88×108μm3，與裂縫數(shù)字巖心模型計(jì)算所得均一致，表明利用中心線提取法結(jié)合最大球填充法建立的裂縫網(wǎng)絡(luò)模型提取的裂縫空間較為完整，證明了該方法的可靠性。

最后將孔隙網(wǎng)絡(luò)模型和裂縫網(wǎng)絡(luò)模型融合為孔隙-裂縫雙重網(wǎng)絡(luò)模型（圖6）。將基于孔隙-裂縫雙重網(wǎng)絡(luò)模型計(jì)算得出的孔隙度和絕對滲透率與常規(guī)物性實(shí)驗(yàn)所得的數(shù)據(jù)進(jìn)行對比，可以看出，模型計(jì)算結(jié)果與實(shí)驗(yàn)測試結(jié)果較為接近（表2），表明孔隙-裂縫雙重網(wǎng)絡(luò)模型可較為準(zhǔn)確地表征孔隙和裂縫的特征。

此外，由于將裂縫近似等效為密集的節(jié)點(diǎn)-鏈接形式的網(wǎng)絡(luò)系統(tǒng)，可以根據(jù)侵入-逾滲理論，進(jìn)一步分析該網(wǎng)絡(luò)系統(tǒng)的兩相以及多相滲流特征，這對于雙重介質(zhì)碳酸鹽巖儲(chǔ)集層的研究具有重要意義。

4 結(jié)論

（1）通過CT掃描及圖像分割可以將孔隙-裂縫型碳酸鹽巖中的孔隙和裂縫進(jìn)行區(qū)分并加以提取，分別建立了孔隙和裂縫的節(jié)點(diǎn)-鏈接形式網(wǎng)絡(luò)模型，從而實(shí)現(xiàn)了對孔隙和裂縫的精準(zhǔn)表征。

（2）基于裂縫數(shù)字巖心模型，利用中心線提取法結(jié)合最大球填充法，能夠建立精細(xì)表征巖石內(nèi)部真實(shí)裂縫空間的裂縫網(wǎng)絡(luò)模型。

（3）利用節(jié)點(diǎn)-鏈接形式的網(wǎng)絡(luò)模型，可以近似表征裂縫的開度等參數(shù)以及裂縫的三維空間展布特征，該模型可用于滲流特征機(jī)理研究。

（4）裂縫網(wǎng)絡(luò)模型采用節(jié)點(diǎn)-鏈接網(wǎng)絡(luò)真實(shí)表征了裂縫空間結(jié)構(gòu)，由于模型的滲流模擬方法為準(zhǔn)靜態(tài)模擬，得到的結(jié)果與常規(guī)物性實(shí)驗(yàn)結(jié)果存在一定誤差。此外在提取裂縫網(wǎng)絡(luò)的過程中，需要驗(yàn)證真實(shí)裂縫信息的保留是否完整，以保證模型的可靠性。

表1 四川盆地雙魚石區(qū)塊孔隙-裂縫型碳酸鹽巖樣品的孔隙和裂縫數(shù)字巖心模型以及孔隙和裂縫網(wǎng)絡(luò)模型Table 1.Digital pore models,digital fracture models,and pore?fracture network models of three carbonate rock samples from Shuangyushi area,Sichuan basin

圖6 四川盆地雙魚石區(qū)塊孔隙-裂縫型碳酸鹽巖樣品的孔隙-裂縫雙重網(wǎng)絡(luò)模型Fig.6.Pore?fracture network models of three carbonate rock samples from Shuangyushi area,Sichuan basin

表2 網(wǎng)絡(luò)模型計(jì)算孔隙度和絕對滲透率與常規(guī)物性實(shí)驗(yàn)結(jié)果對比Table 2.Porosity and absolute permeability calculated from a network model and experimental measurements

符號注釋

A——孔隙橫截面積，μm2；

gij——連通孔隙i和孔隙j的喉道的傳導(dǎo)率，cm3/（mPa·s)；

G——形狀因子；

K——滲透率，mD；

l——喉道長度，cm；

L——模型長度，cm；

pi——孔隙i端壓力，MPa；

pj——孔隙j端壓力，MPa；

P——孔隙橫截面周長，μm；

qij——孔隙i流入孔隙j的流體流量，cm3/s；

Q——模型內(nèi)流體總流量，cm3/s；

r——喉道半徑，cm；

Vb——樣品總體積，μm3；

Vf——裂縫空間體積，μm3；

Vp——孔隙空間體積，μm3；

μ——流體黏度，mPa·s；

φ——孔隙度。