范宜仁, 徐擁軍, 范卓穎, 葛新民

(1.中國石油大學地球科學與技術(shù)學院, 山東 青島 266555; 2.中國石油大學CNPC測井重點實驗室, 山東 青島 266555)

0 引 言

鑄體薄片是將染色樹脂或液態(tài)膠在真空狀態(tài)下灌注到巖石的孔隙空間,在一定的溫度和壓力下使樹脂或液態(tài)膠固結(jié),然后磨制成巖石薄片,再在顯微鏡下觀察孔隙、喉道及其相互連通、配合的二維空間結(jié)構(gòu)的一種孔隙結(jié)構(gòu)評價方法。鑄體薄片對樣品形狀、尺寸要求低,且處理簡單[1-5]。本文介紹了閾值分割、數(shù)學形態(tài)學濾波等算法在鑄體薄片孔隙信息提取、圖像增強等方面的應用,引入兩點相關(guān)函數(shù)法對骨架-孔隙兩相介質(zhì)進行分析,得到了面孔隙度、比表面、平均孔隙半徑、平均顆粒半徑等孔隙參數(shù),并將計算結(jié)果與其他方法進行比較。在此基礎(chǔ)上,重點研究了面孔隙度、比表面、平均顆粒-孔隙半徑比對飽和度指數(shù)的影響,得到與微觀孔隙結(jié)構(gòu)有關(guān)的飽和度指數(shù)模型。

1 鑄體薄片處理及圖像增強方法

根據(jù)鑄體與礦物、黏土的顏色差異及其對光的吸收和敏感性,選擇合適的放大倍數(shù)對鑄體薄片進行成像分析即可得到樣品的形貌照片。鑄體薄片通常以RGB的格式保存,由紅、綠、藍三基色構(gòu)成。由圖像學的理論可知,對一副RGB格式的鑄體薄片其灰度值可以表示為[6]

Ic=fc(x,y) 0≤fc(x,y)≤255

(1)

式中,c表示為R、G、B。

孔隙、骨架和黏土在鑄體薄片中具有不同的灰度特征。國際上通用藍色液態(tài)膠進行薄片的灌注和壓制。如圖1所示,孔隙部分被注入藍色膠體后呈藍色特征。采用閾值分割算法,可將特征明顯的巖石組分從二維圖像中提取出來,閾值分割的操作算子T可表示[6]為

T(tlow,thigh)fc(x,y)=

(2)

式中,tlow、thigh分別表示某巖石組分灰度下、上限值。

圖1 孔隙、骨架及黏土在鑄體薄片中的像素特征

對于孔隙參數(shù)的定量計算,只需將孔隙信息從鑄體薄片中提取出來即可。通過閾值分割后,鑄體薄片變成0或1的像素矩陣,其中0代表骨架,1代表孔隙。

閾值分割算法的缺陷在于,當目標被分割在不同的子塊時會存在明顯的塊狀效應,將孔隙與骨架信息模糊化。數(shù)字形態(tài)學濾波是一種很好的去除塊效應的方法,其基本思想是將圖像看成是點的集合,用結(jié)構(gòu)元素對其進行移位、交、并等集合運算。它具有膨脹、腐蝕、開啟和閉合等4種最基本的運算[7-8]。

設(shè)A為圖像集合,B為結(jié)構(gòu)元素,則結(jié)構(gòu)元素B對圖像集合A的膨脹算法為

A⊕B{B+x∶?A}

(3)

式(3)是將結(jié)構(gòu)元素B對圖像集合A的所有像素作平移,然后對平移得到的結(jié)果作并運算。膨脹運算保留了原圖的基本形狀,但填平了原圖邊界上不平滑的凹陷部分。當結(jié)構(gòu)元素和像素值為1時,即結(jié)構(gòu)元素擊中圖像,輸出1,否則輸出0。

腐蝕算法為

AΘB{x∶B+x?A}

(4)

腐蝕運算將結(jié)構(gòu)元素B平移x,但仍然包含在圖像集合A內(nèi)的所有點x組成。結(jié)構(gòu)元素在整個圖像區(qū)域平移的過程中,結(jié)構(gòu)元素的原點到所有可能的圖像像素點檢查圖像,當結(jié)構(gòu)元素與圖像的前景部分完全匹配時,即結(jié)構(gòu)元素適合圖像,輸出1,否則輸出0。

開啟算法可以表示為

AoB=(AΘB)A⊕B

(5)

A被B開啟運算就是圖像集合A被結(jié)構(gòu)元素B腐蝕后的結(jié)果再被B膨脹。開啟運算可去掉圖像中的一些孤立子域、毛刺和凸出部分,使圖像更加平滑。

閉合運算可以表示為

A.B=A⊕BΘB

(6)

A被B閉合運算就是圖像集合A被結(jié)構(gòu)元素B膨脹后的結(jié)果再被B腐蝕。閉合運算沿著圖像的外邊緣填充,去除了凸向圖像內(nèi)部的尖角,填平縫隙,彌合孔洞和裂縫,能去除圖像上的小孔和凸部,使斷線連接,而總的位置和形態(tài)不變。

一個完整的形態(tài)學濾波過程包括開啟運算和閉合運算,鑄體薄片從孔隙信息提取到形態(tài)學濾波的過程見圖2。

圖2 鑄體薄片數(shù)字圖像處理及形態(tài)學濾波流程

2 相關(guān)函數(shù)理論及孔隙參數(shù)提取

巖石可近似成由骨架和孔隙組成的兩相介質(zhì),可以采用指示函數(shù)法進行骨架和孔隙表征。d維歐氏空間(d)的孔隙-骨架兩相介質(zhì)的指示函數(shù)可定義為[9-14]

(7)

式中,Vi∈d表示歐氏空間內(nèi)介質(zhì)為孔隙;Vi?d表示歐氏空間內(nèi)介質(zhì)為骨架;x為d空間的坐標。

對于用指示函數(shù)表示的骨架-孔隙兩相系統(tǒng),采用相關(guān)函數(shù)表示某相同時存在于空間內(nèi)不同點的概率,對于孔隙,其n點相關(guān)函數(shù)表達式為

Sn(x1,x2,…,xn)=

(8)

式中,< >表示為總體平均。

對于各向同性介質(zhì),n點相關(guān)函數(shù)與各點的絕對坐標無關(guān),僅與各點之間的相對位移有關(guān),則式(8)可轉(zhuǎn)換為[9]

Sn(x1,x2,…,xn)=Sn(x12,…,x1n)

(9)

式中,xij=xj-xi。根據(jù)多點相關(guān)函數(shù)的定義可知,對于孔隙,若n=1,則有

S1==φ

(10)

式中,φ為孔隙度。若n=2,則有

S2(x1,x2)=S2(r)==

(11)

式中,r=x2-x1。

Berryman等[9,11]研究表明,當r=0時,有

S2(0)=S1=φ

(12)

當r→∞時,有

φ2

(13)

從兩點相關(guān)函數(shù)中還可得到比表面積s、平均孔隙半徑rc、平均顆粒半徑rg等。其計算公式分別為

(14)

(15)

Stephen C Blair等[12]認為,相關(guān)函數(shù)值最小時所對應的相關(guān)長度為1.3rg。研究采用該方法確定平均顆粒半徑。

3 實驗分析

經(jīng)閾值分割及數(shù)字形態(tài)學濾波后,鑄體薄片的灰度值為0或1。兩點相關(guān)函數(shù)的離散形式可寫成

(16)

式中,M、N為灰度圖像在x和y方向的點數(shù)。結(jié)合式(12)至式(15)即可得到樣品的相關(guān)函數(shù)圖及面孔隙度、比表面s、平均孔隙半徑rc、平均顆粒半徑rg等。

選取10塊樣品開展了孔隙度、巖電、鑄體薄片及壓汞實驗。實驗環(huán)境為常溫常壓(25 ℃,1 MPa),飽和溶液為NaCl,為降低附加導電對巖電參數(shù)的影響,其濃度配置為50 000 mg/L,孔隙度的測量采用氦氣法,電阻率的測量采用二極法,采用高速冷凍離心機進行巖石失水。圖3為1號樣品的閾值分割、形態(tài)學濾波、兩點相關(guān)函數(shù)計算的全過程。應用同樣的流程對其他樣品進行處理并進行孔隙參數(shù)提取。表1為實驗樣品的常規(guī)物性、巖電及孔隙參數(shù)處理結(jié)果。

圖3 1號樣品鑄體薄片處理及兩點相關(guān)函數(shù)計算結(jié)果

樣號?/%S/%nrp/ms/m-1rc/mrg/mDH/m117.1325.021.467.230.004182.15580.77242.93218.8427.761.281.240.01080.54178.46111.50319.0821.021.431.570.00887.84150.77111.22425.4138.101.355.370.009104.35376.15168.57521.7941.691.221.510.01187.13483.85149.44616.2015.881.576.250.004140.64575.38167.20718.8032.501.725.140.007119.87467.69177.57818.7037.081.491.850.01093.32483.85148.32920.3735.131.486.140.008118.69519.23182.96107.252.242.210.340.00328.76278.4629.42

圖4為10塊樣品的鑄體薄片處理孔隙度與氣測孔隙度的對比。從圖4可知,兩者呈良好的線性關(guān)系,復相關(guān)系數(shù)為0.756,擬合公式為

φs=2.286φ-14.326

(17)

式中,φ為氣測孔隙度,%;φs為薄片處理面孔率,%。

圖4 鑄體薄片處理孔隙度與氣測孔隙度對比

圖5 鑄體薄片處理平均孔隙半徑與壓汞平均孔喉半徑對比

圖5是10塊樣品的鑄體薄片處理的平均孔隙半徑與壓汞法得到的平均孔喉半徑對比。從圖5可知,兩者的線性關(guān)系較好,其復相關(guān)系數(shù)為0.75。由于研究尺度的差別,鑄體薄片處理的平均孔隙半徑要遠遠大于壓汞法平均孔喉半徑。

rc=13.865×rmp+53.529

(18)

式中,rc為鑄體薄片處理平均孔隙半徑,μm;rmp為壓汞法平均孔喉半徑,μm。

以上分析表明,通過數(shù)字圖像處理和兩點相關(guān)函數(shù)法得到的孔隙參數(shù)與其他方法得到的孔隙參數(shù)具有一致性,可用表1的處理結(jié)果進行巖石的電學性質(zhì)分析。圖6是飽和度指數(shù)與比表面的相關(guān)關(guān)系圖。從圖6可知,飽和度指數(shù)隨著比表面的增大而減小,兩者呈良好的冪函數(shù)關(guān)系,其復相關(guān)系數(shù)為0.56,擬合公式為

n=0.38/s0.275

(19)

圖6 飽和度指數(shù)與比表面相關(guān)關(guān)系

巖石的導電及滲流性質(zhì)還受到孔隙半徑、顆粒半徑等因素的共同影響。為了綜合考慮這些參數(shù)對飽和度指數(shù)的貢獻,構(gòu)造平均顆粒-孔隙半徑比表征孔喉的連通程度。圖7是飽和度指數(shù)與平均顆粒-孔隙半徑比的關(guān)系圖。從圖7可知,飽和度指數(shù)隨著平均顆粒-孔隙半徑比的增大而減小,兩者呈良好的線性關(guān)系,復相關(guān)系數(shù)為0.53。經(jīng)擬合,飽和度指數(shù)與平均顆粒-孔隙半徑比的關(guān)系為

(20)

圖7 飽和度指數(shù)與平均顆粒-孔隙半徑比關(guān)系

圖8 飽和度指數(shù)與面孔隙度相關(guān)關(guān)系

圖8是飽和度指數(shù)與面孔隙度的相關(guān)關(guān)系圖。從圖8可知,飽和度指數(shù)隨著孔隙度的增大而減小。經(jīng)擬合,得到飽和度指數(shù)與面孔隙度的關(guān)系為

(21)

結(jié)合式(19)至式(21),將飽和度指數(shù)寫為

(22)

式中,a、b、c、d均為擬合參數(shù),無量綱,可用最小二乘擬合法得到。圖9是根據(jù)式(22)擬合得到的飽和度指數(shù)與實驗測量得到的飽和度指數(shù)對比圖,兩者分布在45 °線左右,經(jīng)統(tǒng)計,其平均絕對誤差為0.067,平均相對誤差為4.7%,精度非常高。

圖9 計算飽和度指數(shù)與實測飽和度指數(shù)對比

4 結(jié) 論

(1) 鑄體薄片含有豐富的孔隙參數(shù)等信息,通過數(shù)字圖像處理與兩點相關(guān)函數(shù)理論得到的面孔隙度、平均孔隙半徑等與其它方法所得的參數(shù)呈良好的線性正相關(guān)關(guān)系。

(2) 飽和度指數(shù)與面孔隙度、巖石比表面積、平均顆粒-孔隙半徑比等密切相關(guān),隨著面孔隙度、巖石比表面積和平均顆粒-孔隙半徑的增大而呈冪函數(shù)衰減。

(3) 通過引入面孔隙度、巖石比表面積及平均顆粒-孔隙半徑比等微觀參數(shù)建立的飽和度指數(shù)模型精度較高,計算的飽和度指數(shù)與實測飽和度指數(shù)基本相等。

參考文獻:

[1] Shelley D. Optical Mineralogy [M]. 2nd ed. New Zealand: Elsevier Press, 1985.

[2] Barber D J. Demountable Polished Extra-thin Sections and Their Use in Transmission Electron Microscopy [J]. Mineralogical Magazine, 1981, 44: 357-359.

[3] Anderson A N, McBratney A B, FitzPatrick E A. Soil mass, Surface, and Spectral Fractal Dimensions Estimated from Thin Section Photographs [J]. Soil Sci Soc Am J, 1996, 60: 962-969.

[4] Grove C, Jerram D A. An ImageJ Macro to Quantify Total Optical Porosity from Blue-stained Thin Sections [J]. Journal of Computers and Geosciences, 2011, 37, jPOR: 1850-1859.

[5] 張麗霞. 高蒙皂石砂巖鑄體薄片制片新方法 [J]. 新疆石油地質(zhì), 2009, 30(6): 709-780.

[6] 岡薩雷斯. 數(shù)字圖像處理 [M]. 北京: 電子工業(yè)出版社, 2003.

[7] 崔屹. 圖象處理與分析-數(shù)學形態(tài)學方法應用 [M]. 北京: 科學出版社, 2000.

[8] Serra J. Image Analysis and Mathematical Morphology [M]. London: Academic Press, 1982.

[9] Berryman J G. Measurement of Spatial Correlation Functions Using Image Processing Techniques [J]. J Appl Phys, 1985, 57: 2374-2384.

[10] Berryman J G, Blair S C. Use of Digital Image Analysis to Estimate Fluid Permeability of Porous Materials: Application of Two Point Correlation Functions [J]. J Appl Phys, 1986, 60: 1930-1938.

[11] Berryman J G. Relationship Between Specific Surface Area and Spatial Correlation Functions for Anisotropic Porous Media [J]. J Math Phys, 1987, 28: 244-245.

[12] Stephen C Blair, Patricia A Berge, James G Berryman. Using Two-point Correlation Functions to Characterize Microgeometry and Estimate Permeabilities of Sandstones and Porous Glass [J]. Journal of Geophysical Research, 1996, 101(B9): 20359-20375.

[13] Ridgway R, Irfanoglu O, Machiraju R, et al. Image Segmentation with Tensor-based Classification ofn-point Correlation Functions [C]∥ Proceedings of the Microscopic Image Analysis with Applications in Biology (MIAAB) Workshop in MICCAI, Dec 2006.

[14] Jiao Y, Stillinger F H, Torquato S. Modeling Heterogeneous Materials Via Two-point Correlation Functions: Basic Principles [J]. Physical Review E (Statistical, Nonlinear, and Soft Matter Physics), 2007, 76(3): 31110-31146.