閆國(guó)亮，孫建孟，劉學(xué)鋒，姜黎明

（1.中國(guó)石油大學(xué)地球科學(xué)與技術(shù)學(xué)院，山東 青島266580；2.中國(guó)石油大學(xué)理學(xué)院，山東 青島266580；3.中國(guó)石油集團(tuán)測(cè)井有限公司技術(shù)中心，陜西 西安710077）

0 引 言

儲(chǔ)層巖石的滲透率與其微觀孔隙結(jié)構(gòu)密切相關(guān)［1］。數(shù)字巖心作為巖石微觀結(jié)構(gòu)的一種數(shù)學(xué)抽象，可以更加真實(shí)地反映巖石的孔喉大小、形狀和拓?fù)浣Y(jié)構(gòu)，提供了研究巖石微觀孔隙結(jié)構(gòu)和滲透率之間關(guān)系的橋梁。為研究巖石微觀孔隙結(jié)構(gòu)對(duì)滲透率的影響規(guī)律，首先需要對(duì)其進(jìn)行定量表征?；跀?shù)字巖心主要通過(guò)提取微觀孔隙空間的孔隙網(wǎng)絡(luò)模型實(shí)現(xiàn)孔隙結(jié)構(gòu)的定量表征。Zhao等［2］采用多向掃描方法對(duì)數(shù)字巖心孔隙空間進(jìn)行多方向切片掃描搜索孔隙和喉道，但該方法很難準(zhǔn)確探測(cè)孔隙的位置。Shin等［3］采用孔隙居中軸線方法研究了孔隙的微觀結(jié)構(gòu)，應(yīng)用該方法可以合理分割孔隙和喉道，但是算法復(fù)雜，人機(jī)交互較多。Oren等［4－5］采用 Voronoi多面體方法提取了過(guò)程模擬法重建數(shù)字巖心的孔隙網(wǎng)絡(luò)模型并進(jìn)行了孔隙結(jié)構(gòu)表征，但只適用于過(guò)程模擬法建立的數(shù)字巖心而不適用于一般數(shù)字巖心的孔隙網(wǎng)絡(luò)建模和表征。Dong［6］采用最大球算法提取了數(shù)字巖心的孔隙網(wǎng)絡(luò)模型，該方法建模速度較快，可用于過(guò)程模擬法和現(xiàn)有其他方法構(gòu)建的數(shù)字巖心的孔隙網(wǎng)絡(luò)模型建模，可以較合理區(qū)分孔隙和喉道，但確定的孔隙長(zhǎng)度偏大而喉道長(zhǎng)度偏小。

本文以過(guò)程模擬法重建數(shù)字巖心作為輸入數(shù)據(jù)，采用改進(jìn)的最大球算法提取了重建數(shù)字巖心的孔隙網(wǎng)絡(luò)模型并對(duì)其孔隙結(jié)構(gòu)進(jìn)行了定量表征；研究了孔隙半徑和喉道半徑等微觀孔隙結(jié)構(gòu)參數(shù)對(duì)滲透率的影響規(guī)律。

1 數(shù)字巖心重建

X射線CT法是建立數(shù)字巖心最直接、最準(zhǔn)確的方法，但由于實(shí)驗(yàn)成本過(guò)高，限制了該項(xiàng)技術(shù)的廣泛應(yīng)用。一般采用數(shù)值重建算法得到數(shù)字巖心。數(shù)值重建算法是通過(guò)巖心二維薄片的統(tǒng)計(jì)分析，得到孔隙度、兩點(diǎn)相關(guān)函數(shù)、粒度分布等統(tǒng)計(jì)信息和黏土含量、礦物組成等巖石特性，利用這些信息重建數(shù)字巖心。數(shù)值重建方法主要有隨機(jī)法和過(guò)程模擬法2種。隨機(jī)法包括高斯場(chǎng)法［7］，模擬退火法［8－9］，順序指示模擬法［10－12］，多點(diǎn)地質(zhì)統(tǒng)計(jì)法［13－15］和馬爾可夫鏈法［16］。當(dāng)孔隙度較小時(shí)，除多點(diǎn)地質(zhì)統(tǒng)計(jì)法和馬爾可夫鏈法之外，隨機(jī)法重建數(shù)字巖心的孔隙連通性較差。過(guò)程模擬法通過(guò)模擬巖石的形成過(guò)程重建數(shù)字巖心，包括沉積過(guò)程、壓實(shí)過(guò)程和成巖過(guò)程的模擬。Oren和 Bakke［5，17］應(yīng)用該方法建立了Fontainebleau砂巖和Berea砂巖的數(shù)字巖心，與對(duì)應(yīng)的X－CT掃描得到的數(shù)字巖心的滲流屬性計(jì)算結(jié)果比較表明，過(guò)程模擬法建立的數(shù)字巖心可以重現(xiàn)真實(shí)巖心的孔隙結(jié)構(gòu)和連通特性。筆者采用Oren等提出的過(guò)程模擬方法，以某地區(qū)砂巖（稱為S1砂巖）巖心的二維薄片為基礎(chǔ)，應(yīng)用自主開(kāi)發(fā)的過(guò)程模擬法程序重建了該砂巖的數(shù)字巖心。

S1砂巖的二維薄片如圖1所示。圖1中白色部分為巖石顆粒，黑色部分為孔隙。采用過(guò)程模擬法重建S1砂巖的數(shù)字巖心結(jié)果如圖2（a）所示，圖2中藍(lán)色部分為巖石骨架，紅色部分為孔隙空間。

圖1 S1砂巖巖心的二維薄片

圖2 S1砂巖的數(shù)字巖心及其孔隙網(wǎng)絡(luò)模型示意圖

2 孔隙網(wǎng)絡(luò)模型提取及孔隙結(jié)構(gòu)定量表征

2.1 孔隙網(wǎng)絡(luò)模型提取

基于過(guò)程模擬法重建的S1砂巖的數(shù)字巖心應(yīng)用改進(jìn)的最大球算法提取數(shù)字巖心的孔隙網(wǎng)絡(luò)模型。首先對(duì)涉及到的基本概念進(jìn)行解釋。

體素：在三維空間中，用以進(jìn)行空間信息的數(shù)據(jù)記錄、處理、表示等所采用的具有一定大小的最小體積單元。

孔隙體素：數(shù)字巖心中表示孔隙部分的最小體積單元。

骨架體素：數(shù)字巖心中表示骨架部分的最小體積單元。

內(nèi)切球：以孔隙體素為球心向四周等速延伸，直到碰到最近的骨架體素，延伸區(qū)域中所有體素的集合。

冗余球：設(shè)B表示內(nèi)切球，如果存在A為任一內(nèi)切球，使得B?A，則稱B為冗余球。

最大球：設(shè)所有內(nèi)切球的集合用I表示，冗余球的集合用M表示，則I－M為最大球集合，其中每個(gè)元素稱為最大球。任意一個(gè)最大球至少包含1個(gè)其他最大球沒(méi)有的體素。

Dong等［6］提出的最大球算法存在的主要問(wèn)題是提取的孔隙網(wǎng)絡(luò)模型的孔隙長(zhǎng)度偏大而喉道長(zhǎng)度偏小，間接影響其他孔喉參數(shù)的確定且計(jì)算的滲透率偏大。筆者針對(duì)這一問(wèn)題，采用計(jì)算機(jī)圖形處理學(xué)中的幾何變換技術(shù)和應(yīng)用統(tǒng)計(jì)學(xué)中的判別分析方法，確定孔隙網(wǎng)絡(luò)模型的孔隙長(zhǎng)度和喉道長(zhǎng)度，從而確定孔隙和喉道的其他參數(shù)。

改進(jìn)后的最大球算法提取數(shù)字巖心孔隙網(wǎng)絡(luò)模型的主要步驟為：

（1）搜尋內(nèi)切球?；跀?shù)字巖心，首先采用擴(kuò)張算法，從第1個(gè)孔隙體素開(kāi)始，以其為中心，從26個(gè)方向（見(jiàn)圖3，紅色方塊表示孔隙體素，灰色方塊表示與其相鄰的26個(gè)體素）尋找距離最近的骨架體素，然后以找到的骨架體素與孔隙體素的距離為最大范圍，采用收縮算法逐一檢查該范圍內(nèi)的體素，確定該孔隙體素對(duì)應(yīng)的內(nèi)切球。

圖3 孔隙體素的26個(gè)搜索方向

采用相同方法尋找下一孔隙體素的內(nèi)切球，直到所有孔隙體素都被遍歷，即可得到所有孔隙體素對(duì)應(yīng)的內(nèi)切球集合。

（2）刪除冗余球。設(shè)A和B分別為內(nèi)切球，它們的球心和半徑分別為CA、CB、RA和RB，且RA＞RB，如果滿足條件

則球B為冗余球，從內(nèi)切球集合中刪除。式（1）中dist（CA，CB）表示內(nèi)切球A和球B兩者球心的距離。

采用上述判別法則從內(nèi)切球集合中刪除所有冗余球。

（3）孔隙喉道識(shí)別。去掉冗余球的內(nèi)切球集合稱為最大球集合，應(yīng)用排序算法和成簇算法［6］識(shí)別孔隙和喉道。采用排序算法將最大球集合中的所有元素按照半徑從大到小排序，根據(jù)半徑將其劃分為一系列的子集，每一個(gè)子集中最大球的半徑相同。對(duì)每一子集中的最大球采用成簇算法確定其屬于孔隙或喉道并得到相應(yīng)的孔隙和喉道的半徑。

（4）孔隙網(wǎng)絡(luò)模型參數(shù)計(jì)算。對(duì)于每一個(gè)孔隙，以該孔隙對(duì)應(yīng)的最大球球心體素為原點(diǎn)，選擇過(guò)該孔隙中心體素的任一直線作為旋轉(zhuǎn)軸。設(shè)平面β繞該旋轉(zhuǎn)軸每隔一定角度旋轉(zhuǎn)，直到轉(zhuǎn)過(guò)180°。應(yīng)用幾何變換技術(shù)可以得到平面β對(duì)該孔隙局部空間剖切的剖面。在剖面內(nèi)每隔一定角度發(fā)出1條射線，用于測(cè)量該方向上孔隙中心體素到巖石骨架體素的距離。對(duì)記錄的所有距離應(yīng)用判別分析方法確定出孔隙的長(zhǎng)度。確定孔隙長(zhǎng)度后即可確定孔隙的體積、形狀因子（形狀因子是指孔隙或喉道截面面積與其周長(zhǎng)平方之比）和喉道的長(zhǎng)度等參數(shù)。喉道體積及形狀因子的確定方法與孔隙體積及形狀因子的確定方法類似。

采用改進(jìn)的最大球算法提取的S1砂巖數(shù)字巖心的孔隙網(wǎng)絡(luò)模型［見(jiàn)圖2（b）］，用球表示孔隙，管表示喉道。根據(jù)形狀因子的不同，實(shí)際計(jì)算中孔隙網(wǎng)絡(luò)模型的孔隙和喉道均為不同截面形狀的柱體。

2.2 孔隙結(jié)構(gòu)定量表征

基于孔隙網(wǎng)絡(luò)模型，統(tǒng)計(jì)孔隙半徑和喉道半徑的頻率分布，得到S1砂巖數(shù)字巖心的孔隙結(jié)構(gòu)參數(shù)分布圖（見(jiàn)圖4和圖5）。由于缺少S1砂巖的恒速壓汞實(shí)驗(yàn)數(shù)據(jù)，因此得到的孔隙結(jié)構(gòu)參數(shù)分布沒(méi)有與實(shí)驗(yàn)數(shù)據(jù)對(duì)比，但經(jīng)過(guò)數(shù)據(jù)擬合分析，S1砂巖數(shù)字巖心的孔隙半徑和喉道半徑等微觀孔隙結(jié)構(gòu)參數(shù)滿足對(duì)數(shù)正態(tài)分布，這種分布特征與Morrow［18］得到的結(jié)論一致。

圖4 S1砂巖數(shù)字巖心的孔隙半徑頻率分布

圖5 S1砂巖數(shù)字巖心的喉道半徑頻率分布

3 滲透率數(shù)值模擬

應(yīng)用逾滲網(wǎng)絡(luò)理論［19］可求得孔隙網(wǎng)絡(luò)模型i方向（可取x，y，z）的流量，根據(jù)Darcy定律計(jì)算滲透率

式中，Ki為孔隙網(wǎng)絡(luò)模型i方向的滲透率，mD＊非法定計(jì)量單位，1mD＝9.87×10－4μm2，下同；Li為i方向長(zhǎng)度，cm；Ai為i方向截面積，cm2；pI－pO為i方向壓差，×10－1MPa；μ為流體黏度，mPa·s；Qi為i方向流量，cm3／s。

通過(guò)對(duì)比S1砂巖實(shí)驗(yàn)測(cè)量的滲透率和改進(jìn)后最大球算法提取的孔隙網(wǎng)絡(luò)模型計(jì)算的滲透率評(píng)價(jià)上述方法的準(zhǔn)確性。S1砂巖實(shí)驗(yàn)測(cè)量的孔隙度為19.70%，平均滲透率為1 286mD。S1砂巖孔隙網(wǎng)絡(luò)模型的孔隙度為19.78%，計(jì)算的平均滲透率為1 362mD，與S1砂巖實(shí)驗(yàn)測(cè)量結(jié)果接近，說(shuō)明了該方法的正確性。同時(shí)，原有最大球算法提取的S1砂巖的孔隙網(wǎng)絡(luò)模型的孔隙度為19.78%，計(jì)算的平均滲透率為1 531mD，與改進(jìn)后的結(jié)果相比，改進(jìn)后最大球算法提取的孔隙網(wǎng)絡(luò)模型的滲透率更接近實(shí)驗(yàn)測(cè)量結(jié)果。

根據(jù)孔隙結(jié)構(gòu)表征部分孔隙半徑和喉道半徑頻率分布，采用式（3）計(jì)算得到滲透率貢獻(xiàn)值曲線［20］

式中，ΔKj為第j個(gè)區(qū)間的滲透率貢獻(xiàn)值；rj為第j個(gè)區(qū)間的孔隙或喉道半徑值；αj為第j個(gè)區(qū)間的孔隙半徑或喉道半徑頻率。

選取滲透率貢獻(xiàn)最大值對(duì)應(yīng)的孔隙半徑和喉道半徑作為S1砂巖孔隙和喉道半徑的特征值。分別改變孔隙網(wǎng)絡(luò)模型的孔隙半徑分布和喉道半徑分布，得到圖6和圖7所示的S1砂巖的孔隙半徑特征值和喉道半徑特征值與滲透率的關(guān)系曲線?？梢?jiàn)喉道半徑特征值對(duì)滲透率的影響大于孔隙半徑特征值對(duì)滲透率的影響。滲透率隨孔隙半徑特征值和喉道半徑特征值的變化規(guī)律均可以用Logistic函數(shù)關(guān)系描述

式中，K為滲透率，mD；r為孔隙半徑或喉道半徑特征值，μm；c為曲線趨于穩(wěn)定段的滲透率，mD；a為滲透率快速增加段中點(diǎn)對(duì)應(yīng)的孔隙半徑或喉道半徑特征值，μm；p表征滲透率隨孔隙半徑或喉道半徑特征值變化的快慢，取值一般在2～3之間，無(wú)量綱。

圖6 S1砂巖孔隙半徑特征值與滲透率關(guān)系

圖7 S1砂巖喉道半徑特征值與滲透率關(guān)系

4 結(jié) 論

（1）在巖心二維薄片的基礎(chǔ)上，采用過(guò)程模擬法重建的數(shù)字巖心不僅能夠反映巖石的微觀孔隙結(jié)構(gòu)特征，而且具有與真實(shí)巖心相近的滲流特性。

（2）通過(guò)孔隙半徑和喉道半徑等孔隙結(jié)構(gòu)參數(shù)的定量表征，表明S1砂巖微觀孔隙結(jié)構(gòu)參數(shù)的分布基本滿足對(duì)數(shù)正態(tài)函數(shù)分布。

（3）S1砂巖的孔隙半徑和喉道半徑特征值與滲透率之間的關(guān)系可以用Logistic函數(shù)關(guān)系描述，且喉道半徑特征值對(duì)滲透率的影響大于孔隙半徑特征值對(duì)滲透率的影響。

［1］ Dullien F A.Porous Media Fluid Transport and Pore Structure［M］.New York：Academic Press，1992.

［2］ Zhao H Q，Macdonald I F，Kwiecien M J.Multi－orientation Scanning：A Necessity in the Identification of Pore Necks in Porous Media by 3－D Computer Reconstruction from Serial Section Data ［J］.Journal of Colloid and Interface Science，1994.162（2）：390－401.

［3］ Shin H，Lindquist W B，Sahagian D L，et al.Analysis of the Vesicular Structure of Basalts［J］.Computers and Geosciences，2005，31：473－487.

［4］ Bakke S，Oren P E.3－D Pore－scale Modeling of Sandstones and Flow Simulations in the Pore Networks［C］∥SPE 35479，1997.

［5］ Oren P E，Bakke S.Reconstruction of Berea Sandstone and Pore－scale Modeling of Wettability Effects［J］.Journal of Petroleum Science and Engineering，2003，39（3－4）：177－199.

［6］ Dong H.Micro－CT Imaging and Pore Network Extraction［D］.London：Imperial College，2007.

［7］ Joshi M.A Class Three－dimensional Modeling Technique for Studying Porous Media［D］.Kansas：University of Kansas，1974.

［8］ Hazlett R D.Statistical Characterization and Stochastic Modeling of Pore Networks in Relation to Fluid Flow ［J］.Mathematical geology，1997，29（6）：801－822.

［9］ 趙秀才，姚軍，陶軍，等.基于模擬退火算法的數(shù)字巖心建模方法 ［J］.高校應(yīng)用數(shù)學(xué)學(xué)報(bào)：A輯，2007，22（2）：127－133.

［10］ Keehm Y.Computational Rock Physics：Transport Properties in Porous Media and Applications［D］.Stanford：Stanford University，2003.

［11］ 朱益華，陶果.順序指示模擬技術(shù)及其在3D數(shù)字巖心建模中的應(yīng)用 ［J］.測(cè)井技術(shù)，2007，31（2）：112－115.

［12］ 劉學(xué)鋒，孫建孟，王海濤，等.順序指示模擬重建三維數(shù)字巖心的準(zhǔn)確性評(píng)價(jià) ［J］.石油學(xué)報(bào)，2009，30（3）：391－395.

［13］ Okabe H，Blunt M J. Pore Space Reconstruction Using Multiple－point Statistics［J］.Journal of Petroleum Science and Engineering，2005，46：121－137.

［14］ 張挺，盧德唐，李道倫.基于二維圖像和多點(diǎn)統(tǒng)計(jì)方法的多孔介質(zhì)重構(gòu)研究 ［J］.中國(guó)科技大學(xué)學(xué)報(bào)，2010，40（3）：271－277.

［15］ 張麗，孫建孟，孫志強(qiáng)，等.多點(diǎn)地質(zhì)統(tǒng)計(jì)學(xué)在三維巖心孔隙分布建模中的應(yīng)用 ［J］.中國(guó)石油大學(xué)學(xué)報(bào)：自然科學(xué)版，2012，36（2）：105－109.

［16］ Wu K，Nunan N，Crawford J W，et al.An Efficient Markov Chain Model for the Simulation of Heterogeneous Soil Structure ［J］.Soil Science Society of America，2004，68：346－351.

［17］ Oren P E，Bakke S.Process Based Reconstruction of Sandstones and Predictions of Transport Properties［J］.Transport in Porous Media，2002，46（2－3）：311－343.

［18］ Morrow N R.Interfacial Phenomena in Petroleum Recovery［M］.New York：Marcel Dekker Inc，1991.

［19］ Valvatne P H.Predictive Pore－scale Modeling of Multiphase Flow［D］.London：Imperial College，2004.

［20］ 羅蟄潭，王允誠(chéng).油氣儲(chǔ)集層的孔隙結(jié)構(gòu) ［M］.北京：科學(xué)出版社，1986.