羅瑜，王寅，王容，袁雯

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

對于裂縫和孔隙均發(fā)育的碳酸鹽巖滲流特征的研究，通常是構(gòu)建裂縫和具有一般孔隙介質(zhì)的巖塊系統(tǒng)，利用有限元法耦合流體在孔隙和裂縫中的流動特征[1-3]。此種方法運(yùn)算過程較為復(fù)雜，且對于非均質(zhì)性強(qiáng)的樣品所構(gòu)建的模型體量較大，對計算設(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)建的模型，分別計算了3 塊孔隙-裂縫型碳酸鹽巖樣品的孔隙度和滲透率，與常規(guī)物性實(shí)驗所得數(shù)據(jù)進(jìn)行對比，對模型的可靠性進(jì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ù)孔隙和裂縫不同的幾何形狀，計算費(fèi)雷特長寬比（單個孔隙或裂縫的最長直徑與最短直徑之比）、比表面積等參數(shù)，選取合適的閾值，有效分離出孔隙和裂縫（圖1c，圖1d）。

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

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

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

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

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

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

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

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

形狀因子可補(bǔ)充具有不規(guī)則形狀的真實(shí)孔隙空間在提取孔隙網(wǎng)絡(luò)模型時被簡化成節(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)行裂縫中心線的提?。菏紫龋嬎惴蛛x出的裂縫空間數(shù)字巖心模型中的每一個體素到最近巖石骨架的距離；然后，根據(jù)距離圖法，不斷細(xì)化裂縫空間的數(shù)字巖心模型，直至形成一條連續(xù)的體素鏈，該體素鏈即為提取出的裂縫中心線。

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

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

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

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

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

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

3 模型驗證

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

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

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

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

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í)驗結(jié)果存在一定誤差。此外在提取裂縫網(wǎng)絡(luò)的過程中，需要驗證真實(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ò)模型計算孔隙度和絕對滲透率與常規(guī)物性實(shí)驗結(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；

φ——孔隙度。