儲層巖石微觀孔隙結構特征及其對滲透率影響

閆國亮，孫建孟，劉學鋒，姜黎明

（1.中國石油大學地球科學與技術學院，山東 青島266580；2.中國石油大學理學院，山東 青島266580；3.中國石油集團測井有限公司技術中心，陜西 西安710077）

0 引 言

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

本文以過程模擬法重建數字巖心作為輸入數據，采用改進的最大球算法提取了重建數字巖心的孔隙網絡模型并對其孔隙結構進行了定量表征；研究了孔隙半徑和喉道半徑等微觀孔隙結構參數對滲透率的影響規律。

1 數字巖心重建

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

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

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

圖2 S1砂巖的數字巖心及其孔隙網絡模型示意圖

2 孔隙網絡模型提取及孔隙結構定量表征

2.1 孔隙網絡模型提取

基于過程模擬法重建的S1砂巖的數字巖心應用改進的最大球算法提取數字巖心的孔隙網絡模型。首先對涉及到的基本概念進行解釋。

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

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

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

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

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

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

Dong等［6］提出的最大球算法存在的主要問題是提取的孔隙網絡模型的孔隙長度偏大而喉道長度偏小，間接影響其他孔喉參數的確定且計算的滲透率偏大。筆者針對這一問題，采用計算機圖形處理學中的幾何變換技術和應用統計學中的判別分析方法，確定孔隙網絡模型的孔隙長度和喉道長度，從而確定孔隙和喉道的其他參數。

改進后的最大球算法提取數字巖心孔隙網絡模型的主要步驟為：

（1）搜尋內切球?；跀底謳r心，首先采用擴張算法，從第1個孔隙體素開始，以其為中心，從26個方向（見圖3，紅色方塊表示孔隙體素，灰色方塊表示與其相鄰的26個體素）尋找距離最近的骨架體素，然后以找到的骨架體素與孔隙體素的距離為最大范圍，采用收縮算法逐一檢查該范圍內的體素，確定該孔隙體素對應的內切球。

圖3 孔隙體素的26個搜索方向

采用相同方法尋找下一孔隙體素的內切球，直到所有孔隙體素都被遍歷，即可得到所有孔隙體素對應的內切球集合。

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

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

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

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

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

采用改進的最大球算法提取的S1砂巖數字巖心的孔隙網絡模型［見圖2（b）］，用球表示孔隙，管表示喉道。根據形狀因子的不同，實際計算中孔隙網絡模型的孔隙和喉道均為不同截面形狀的柱體。

2.2 孔隙結構定量表征

基于孔隙網絡模型，統計孔隙半徑和喉道半徑的頻率分布，得到S1砂巖數字巖心的孔隙結構參數分布圖（見圖4和圖5）。由于缺少S1砂巖的恒速壓汞實驗數據，因此得到的孔隙結構參數分布沒有與實驗數據對比，但經過數據擬合分析，S1砂巖數字巖心的孔隙半徑和喉道半徑等微觀孔隙結構參數滿足對數正態分布，這種分布特征與Morrow［18］得到的結論一致。

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

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

3 滲透率數值模擬

應用逾滲網絡理論［19］可求得孔隙網絡模型i方向（可取x，y，z）的流量，根據Darcy定律計算滲透率

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

通過對比S1砂巖實驗測量的滲透率和改進后最大球算法提取的孔隙網絡模型計算的滲透率評價上述方法的準確性。S1砂巖實驗測量的孔隙度為19.70%，平均滲透率為1 286mD。S1砂巖孔隙網絡模型的孔隙度為19.78%，計算的平均滲透率為1 362mD，與S1砂巖實驗測量結果接近，說明了該方法的正確性。同時，原有最大球算法提取的S1砂巖的孔隙網絡模型的孔隙度為19.78%，計算的平均滲透率為1 531mD，與改進后的結果相比，改進后最大球算法提取的孔隙網絡模型的滲透率更接近實驗測量結果。

根據孔隙結構表征部分孔隙半徑和喉道半徑頻率分布，采用式（3）計算得到滲透率貢獻值曲線［20］

式中，ΔKj為第j個區間的滲透率貢獻值；rj為第j個區間的孔隙或喉道半徑值；αj為第j個區間的孔隙半徑或喉道半徑頻率。

選取滲透率貢獻最大值對應的孔隙半徑和喉道半徑作為S1砂巖孔隙和喉道半徑的特征值。分別改變孔隙網絡模型的孔隙半徑分布和喉道半徑分布，得到圖6和圖7所示的S1砂巖的孔隙半徑特征值和喉道半徑特征值與滲透率的關系曲線?？梢姾淼腊霃教卣髦祵B透率的影響大于孔隙半徑特征值對滲透率的影響。滲透率隨孔隙半徑特征值和喉道半徑特征值的變化規律均可以用Logistic函數關系描述

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

圖6 S1砂巖孔隙半徑特征值與滲透率關系

圖7 S1砂巖喉道半徑特征值與滲透率關系

4 結 論

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

（2）通過孔隙半徑和喉道半徑等孔隙結構參數的定量表征，表明S1砂巖微觀孔隙結構參數的分布基本滿足對數正態函數分布。

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

［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］ 趙秀才，姚軍，陶軍，等.基于模擬退火算法的數字巖心建模方法 ［J］.高校應用數學學報：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］ 朱益華，陶果.順序指示模擬技術及其在3D數字巖心建模中的應用 ［J］.測井技術，2007，31（2）：112－115.

［12］ 劉學鋒，孫建孟，王海濤，等.順序指示模擬重建三維數字巖心的準確性評價 ［J］.石油學報，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］ 張挺，盧德唐，李道倫.基于二維圖像和多點統計方法的多孔介質重構研究 ［J］.中國科技大學學報，2010，40（3）：271－277.

［15］ 張麗，孫建孟，孫志強，等.多點地質統計學在三維巖心孔隙分布建模中的應用 ［J］.中國石油大學學報：自然科學版，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］ 羅蟄潭，王允誠.油氣儲集層的孔隙結構 ［M］.北京：科學出版社，1986.