姜黎明, 趙強先, 李兵, 蔡成定, 冀昆, 閆國亮

(1.中國石油集團測井有限公司, 陜西 西安 710077; 2.中國石油天然氣集團公司測井重點實驗室,陜西 西安 710077; 3.中國石油勘探開發(fā)研究院西北分院, 甘肅 蘭州 730020)

0 引 言

孔隙結構對巖石的電性、滲流、核磁共振等物理特性具有重要的影響[1-5],實驗室內對巖石微觀孔隙結構的描述和評價方法(毛細管壓力曲線、鑄體薄片、掃描電鏡等)已經日趨成熟,但存在一些缺陷。例如毛細管壓力曲線法不夠形象直觀,特別是壓汞對巖心樣品造成的損耗嚴重;鑄體薄片和掃描電鏡得到的二維圖像空間延續(xù)性差。隨著CT掃描技術在石油行業(yè)的廣泛應用,形成了基于CT成像的三維高精度儲層表征技術,其對于儲層微觀孔隙結構的研究具有重大意義。

對巖石三維復雜孔隙結構進行定量表征最新的方法是基于數字巖心,通過提取微觀孔隙空間的孔隙網絡模型進行定量表征。Zhao等[6-7]采用多向掃描方法對數字巖心孔隙空間進行多方向切片掃描搜索孔隙和喉道,該方法很難準確探測孔隙。Lindquist等[8-9]采用孔隙居中軸線方法研究了孔隙的微觀結構,應用該方法可以合理分割孔隙和喉道,但算法復雜,人機交互繁瑣。?ren等[10-12]采用Voronoi多面體方法提取了過程模擬法重建數字巖心的孔隙網絡模型并進行了孔隙結構表征,但只適用于過程模擬法建立的數字巖心而不適用于一般數字巖心的孔隙網絡建模和表征。Dong等[13-14]采用最大球算法提取了數字巖心的孔隙網絡模型,該方法建模速度較快,能劃分出孔隙和喉道,但最終確定的孔隙長度偏大而喉道長度偏小,間接影響其他孔隙和喉道參數的確定。本文在最大球算法的基礎上通過引入圖形幾何變換和判別分析2項關鍵技術,實現(xiàn)了孔隙、喉道的合理劃分,準確提取巖石的孔徑分布。

圖1 X射線CT掃描樣品展示(從左起砂巖、碳酸鹽巖、頁巖)

1 利用CT掃描構建三維數字巖心

研究中采用的X射線CT設備為美國GE公司v|tome|xs 180型微米級CT掃描儀系統(tǒng),最高分辨率可達1 μm。對巖心進行CT掃描,具體流程:巖樣準備—儀器抽真空及預熱—建立工程文件—調整巖樣檢測位置—設置X射線參數—進行探測器校準—設置掃描參數—巖樣掃描及數據保存—掃描數據重建—導出掃描圖片—編寫測試報告。

(1) 巖心準備。將待掃描的巖樣進行洗油、洗鹽處理,去除結晶在巖心喉道里面的鹽。該儀器提供了一系列尺寸的樣品夾持器,樣品的大小與X射線CT掃描分辨率有關。掃描分辨率越高,掃描樣品的尺寸越小,可以根據研究需求選擇合適的分辨率及樣品尺寸。

(2) 儀器預熱、載物臺位置調節(jié)以及X射線參數設置。

(3) 分辨率優(yōu)選。掃描分辨率的選擇是利用CT掃描重構巖心三維圖像關鍵的一步,直接關系在三維圖像中能否準確直觀反映巖心的孔隙結構。若分辨率過高,則掃描樣品的尺寸太小,導致三維圖像中無法反映真實巖心中的大尺寸孔隙。相反,若分辨率過低,則無法識別巖心中的小尺寸孔隙。

分辨率優(yōu)選的方法一般有2種,對于相對比較均質的巖心,可以統(tǒng)計不同分辨率下三維巖心圖像的孔隙度與實驗孔隙進行對比,從而確定最佳掃描分辨率。對于非均質性強的巖心,不僅需要通過孔隙度選擇掃描分辨率,還應對比三維數字巖心的滲透率模擬結果與實驗結果。

(4) 掃描過程。上述參數設定完畢后,開始掃描。

(5) 三維圖像重建。首先通過X射線CT掃描得到與樣品X射線吸收系數有關的投影數據,然后通過重建算法,將其轉換為一系列的巖心橫截面二維圖像,將巖心的二維橫截面圖像組合起來便得到巖心的三維灰度圖像,最后通過對灰度圖像進行圖像分割,得到研究所需的三維二值圖像即三維數組。圖1給出了砂巖、碳酸鹽巖、壓裂后頁巖利用X射線CT掃描圖片,樣品尺寸分別為5、25、25 mm,掃描分辨率分別為3.6、14.5、14.8 μm/像素點。

2 孔隙網絡模型提取算法改進

巖石孔隙結構參數的計算通過基于X射線CT掃描三維圖像提取孔隙網絡模型實現(xiàn),采用方法為最大球算法。具體實現(xiàn)過程:在三維數字巖心的孔隙像素中任選1點,以其為球心不斷增大球半徑向四周延伸,直至球表面觸碰到最近的骨架像素為止,形成的區(qū)域中所包含的所有像素的集合稱為最大球。1個最大球可以重疊1個半徑小于它的相鄰最大球,這個相鄰的最大球又可以重疊半徑小于它的相鄰最大球,這樣形成1個最大球多簇。將簇中所有被包含的重復的球體刪除,保證任意1個最大球至少包含1個其他最大球沒有的體素。此時,數字巖心中所有的孔隙就被不同的最大球多簇填滿。

圖2 最大球多簇定義的孔隙、喉道示意圖

圖3 任意平面繞直線P1P2旋轉θ角示意圖

將最大球多簇中半徑最大的球體定義為該簇的祖先,其限定該簇所占據的孔隙。如果1個多簇中某一最大球具有2個祖先,這個公共的最大球便被定義為喉道。喉道一旦確定,便形成了2個樹狀結構(亦稱作多簇)(見圖2)。在孔隙喉道鏈中,一般以祖先半徑的0.7倍作為孔隙部分和喉道部分的界限,所有鏈中小于0.7倍祖先半徑的部分設定為喉道部分,大于等于0.7倍祖先半徑部分被設定為孔隙部分。

最大球算法存在的問題:①0.7這個數的選擇比較隨意,缺乏依據;②最終確定的孔隙長度偏大而喉道長度偏小,間接影響其他孔隙和喉道參數確定。

(1) 圖形幾何變換。圖形的幾何變換是指對圖形的幾何信息經過平移、縮放、旋轉等變換后產生新的圖形[15]。本文在計算孔隙長度時需要利用平面對孔隙進行剖切處理,剖切過程中平面要連續(xù)旋轉,由于平面在三維空間中展布的任意性,這給實際問題的處理帶來很大的不便,利用三維復合變換可以解決這一問題。

設直線P1P2上2點P1、P2的坐標分別為(x1,y1,z1)、(x2,y2,z2),平面π1上的任意1點P(x,y,z)繞直線P1P2旋轉θ角后的點P′的坐標為(x′,y′,z′)(見圖3),則可以用式(1)描述變換前后的坐標關系

(1)

式中,Rp為待求的三維復合變換矩陣。

第1步:把直線P1P2平移至原點,相應的平移變換矩陣為

(2)

cosφ=c/d1sinφ=b/d1

(3)

則繞x軸的旋轉變換矩陣為

Rx(φ)=

(4)

第3步:繞y軸旋轉,使直線與z軸重合。

cosα=d1/d2sinα=-a/d2

(5)

則繞y軸的旋轉變換矩陣為

Ry(α)=

(6)

經該變化后,直線P1P2已經于z軸重合。

第4步:繞z軸旋轉θ角。

完成上述一系列平移和旋轉變換后,平面π1上的任意一點P(x,y,z)繞直線P1P2旋轉θ角就變成繞z軸旋轉θ角。繞z軸旋轉θ角的變換矩陣為

(7)

第5步:執(zhí)行步驟(3)、步驟(2)、步驟(1)的逆變換。

圖4 孔隙剖面和孔隙長度示意圖

完成第4步的旋轉變換后,還需要將旋轉軸變回到原來P1P2的位置。此時只需要求Ry(α)、Rx(φ)、T3(-x1,-y1,-z1)的逆矩陣Ry(-α)、Rx(-φ)、T3(x1,y1,z1)即可。

(8)

(9)

(10)

第6步:求得三維復合變換矩陣Rp。

綜上,繞直線P1P2旋轉θ角的變換矩陣Rp為

Rp=T3(-x1,-y1,-z1)·Rx(φ)·Ry(α)·

Rz(θ)·Ry(-α)·Rx(-φ)·T3(x1,y1,z1)

(11)

平面π1上的任意一點P(x,y,z)繞直線P1P2旋轉θ角后的點P′的坐標為(x′,y′,z′)均可由式(1)、式(10)和式(11)求得。

(2) 判別分析法。實驗或計算得到的一系列數據樣品劃分成2組,需要確定合理的閾值,則需要用到判別分析方法。判別分析方法是將1組數據按照從小到大的順序排列,選定一個初始閾值,將數據分成2組,計算2組數據的組內方差和組間方差,使得組內方差和組間方差之比為最大以確定合適的閾值。

設給定的閾值為k,將數據分成2組。計算第1組數據的樣品數N1(k)、平均值μ1(k)和方差σ1(k);計算第2組數據的樣品數N2(k)、平均值μ2(k)和方差σ2(k)。同時計算所有樣品的平均值μK,則組內方差σW和組間方差σB分別為

(12)

(13)

如果閾值k使得σW/σB達到最大,則k即為所求的合理閾值。在下文孔隙長度分割的過程中將采用判別分析方法確定孔隙長度的閾值。

在樹狀結構的孔隙空間中選取一平面對數字巖心孔隙空間進行剖切處理,應用幾何變換技術可以把剖切得到的局部孔隙空間記錄下來。圖4(a)就是剖切過程中得到的1個切面,圖4(a)中紅星表示孔隙中心。對于剖切得到的每1個平面,從孔隙中心每隔一定角度發(fā)射1條射線,射線不斷延伸直到碰到骨架體素為止[見圖4(b)]。計算每1條射線段的長度,這樣就可以形成1個度量孔隙局部空間尺度的數據集合?？梢钥闯?總有一些射線可以穿過與當前孔隙相連的喉道進入相鄰孔隙,這類射線段的長度較大,不能準確表征當前孔隙的長度。因此,在確定孔隙長度的過程中,不能對所有射線段取平均值,而應該采用判別分析法確定孔隙喉道的長度。

3 效果驗證

利用改進后的最大球方法得到了某砂巖孔隙和喉道長度分布(見圖5)。通過與原算法計算結果對比發(fā)現(xiàn),改進后方法得到的孔隙長度分布峰值左移,喉道長度分布峰右移,說明改進后方法得到的孔隙長度變短,喉道長度變長。因此,改進后方法克服了原最大球方法在孔隙空間分割上存在的問題。進一步對算法處理效果進行對比,利用改進前后算法分別對具有雙重孔隙空間的巖心進行了孔徑提取,計算結果見圖6。從圖6上可以看出改進后算法計算孔徑精度明顯優(yōu)于原算法。

為進一步檢驗改進最大球算法提取巖心孔徑分布的效果,選取了某油田9塊低孔隙度低滲透率、中高孔隙度與滲透率砂巖巖心,孔隙度為7.5%～18.5%,平均孔隙度為13.2%,滲透率為4.69～523 mD*非法定計量單位,1 mD=9.87×10-4μm2,下同,平均為158.6 mD。對該批巖樣進行了孔徑計算,并把計算結果與核磁共振T2結果進行了比較。圖7給出了2塊巖心的計算結果。通過對比發(fā)現(xiàn)數字巖心計算結果與核磁共振實驗結果基本吻合,低孔隙度低滲透率巖心誤差在15%以內,中高孔隙度與滲透率巖心誤差控制在10%以(見圖8)。改進的最大球算法計算的孔徑分布曲線比原最大球算法計算精度有了很大的提升。

圖5 孔喉長度計算結果比較

圖6 孔徑分布計算結果比較

圖7 孔徑分布計算結果比較

圖8 孔徑分布計算誤差分析

4 結 論

(1) 給出了CT掃描均質性和非均質性巖心時選擇最佳分辨率的方法。對于均質的巖心,可以通過圖像孔隙度與實驗孔隙比對的方式確定最佳掃描分辨率。對于非均質性強的巖心,同時考察孔隙度和滲透率的圖像模擬結果與實驗結果的對比,確定最佳掃描分辨率。

(2) 在最大球算法的基礎上,引入圖形幾何變換和判別分析2項關鍵算法,實現(xiàn)了孔隙、喉道的合理劃分,準確提取巖石的孔徑分布。通過把新方法計算結果與原算法計算結果、核磁共振實驗結果進行對比,發(fā)現(xiàn)巖心孔徑計算精度有了明顯提高。這為CT掃描圖像數據在測井解釋評價中的應用奠定了基礎。

參考文獻:

[1] SUN J M, LIU X F, WANG H T. Improved Numerical Simulation of Electrical Properties of Reservoir Rock Using Morphology [C]∥SPWLA 50th Annual Logging Symposium, June 21-24, 2009.

[2] JIANG L M, SUN J M, LIU X F, et al. Study of Different Factors Affecting the Electrical Properties of Natural Gas Reservoir Rocks Based on Digital Cores [J]. Journal of Geophysics and Engineering, 2011, 8(2): 366-371.

[3] 毛志強, 高楚橋. 孔隙結構與含油巖石電阻率性質理論模擬研究 [J]. 石油勘探與開發(fā), 2000(2): 87-90.

[4] BLUNT M J. Flow in Porous Media Pore-network Models and Multiphase Flow [J]. Current Opinion in Colloid & Interface Science, 2001, V6(3): 197-207.

[5] 王克文, 關繼騰, 范業(yè)活, 等. 孔隙網絡模型在滲流力學研究中的應用 [J]. 力學進展, 2005(3): 353-360.

[6] 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.

[7] LIANG Z, IOANNIDIS M A, CHATZIS I. Geometric and Topological Analysis of Three-dimensional Porous Media: Pore Space Partitioning Based on Morphological Skeletonization [J]. Journal of Colloid and Interface Science, 2000, 221: 13-24.

[8] LINDQUIST W B, LEE S M, COKER D, et al. Medial Axis Analysis of Void Structure in Three-dimensional Tomographic Images of Porous Media [J]. Journal of Geophysical Research, 1996, 101B: 8297.

[9] LINDQUIST W B, VENKATARANGAN A. Investigating 3D Geometry of Porous Media from High Resolution Images [J]. Phys. Chem. Earth(A), 1999, 25(7): 593-599.

[10] 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.

[11] 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.

[12] DELERUE J F C, PERRIER E. DXSoil, a Library for 3d Image Analysis in Soil Science [J]. Computers & Geosciences, 2002, 28: 1041-1050.

[13] SILIN D B, JIN G, PATZEK T W. Robust Determination of Pore Space Morphology in Sedimentary Rocks [C]∥Paper SPE 84296, Proceedings of SPE Annual Technical Conference and Exhibition, Denver, Colorado, U. S. A, 2003.

[14] DONG H, BLUNT M J. Pore-network Extraction from Micro Computerized Tomography Images [J]. Physical Review E, 2009, 80(3): 36307.

[15] 朱虹. 數字圖像處理基礎 [M]. 北京: 科學出版社, 2005.