鑄體薄片孔隙參數(shù)提取及其對(duì)飽和度指數(shù)的影響

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

(1.中國石油大學(xué)地球科學(xué)與技術(shù)學(xué)院, 山東 青島 266555; 2.中國石油大學(xué)CNPC測(cè)井重點(diǎn)實(shí)驗(yàn)室, 山東 青島 266555)

0 引 言

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

1 鑄體薄片處理及圖像增強(qiáng)方法

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

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

(1)

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

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

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

(2)

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

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

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

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

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

A⊕B{B+x∶?A}

(3)

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

腐蝕算法為

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

(4)

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

開啟算法可以表示為

AoB=(AΘB)A⊕B

(5)

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

閉合運(yùn)算可以表示為

A.B=A⊕BΘB

(6)

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

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

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

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

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

(7)

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

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

Sn(x1,x2,…,xn)=

(8)

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

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

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

(9)

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

S1==φ

(10)

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

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

(11)

式中,r=x2-x1。

Berryman等[9,11]研究表明,當(dāng)r=0時(shí),有

S2(0)=S1=φ

(12)

當(dāng)r→∞時(shí),有

φ2

(13)

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

(14)

(15)

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

3 實(shí)驗(yàn)分析

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

(16)

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

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

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

樣號(hào)?/%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塊樣品的鑄體薄片處理孔隙度與氣測(cè)孔隙度的對(duì)比。從圖4可知,兩者呈良好的線性關(guān)系,復(fù)相關(guān)系數(shù)為0.756,擬合公式為

φs=2.286φ-14.326

(17)

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

圖4 鑄體薄片處理孔隙度與氣測(cè)孔隙度對(duì)比

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

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

rc=13.865×rmp+53.529

(18)

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

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

n=0.38/s0.275

(19)

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

巖石的導(dǎo)電及滲流性質(zhì)還受到孔隙半徑、顆粒半徑等因素的共同影響。為了綜合考慮這些參數(shù)對(duì)飽和度指數(shù)的貢獻(xiàn),構(gòu)造平均顆粒-孔隙半徑比表征孔喉的連通程度。圖7是飽和度指數(shù)與平均顆粒-孔隙半徑比的關(guān)系圖。從圖7可知,飽和度指數(shù)隨著平均顆粒-孔隙半徑比的增大而減小,兩者呈良好的線性關(guān)系,復(fù)相關(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í)驗(yàn)測(cè)量得到的飽和度指數(shù)對(duì)比圖,兩者分布在45 °線左右,經(jīng)統(tǒng)計(jì),其平均絕對(duì)誤差為0.067,平均相對(duì)誤差為4.7%,精度非常高。

圖9 計(jì)算飽和度指數(shù)與實(shí)測(cè)飽和度指數(shù)對(duì)比

4 結(jié) 論

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

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

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

參考文獻(xiàn):

[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ù)學(xué)形態(tài)學(xué)方法應(yīng)用 [M]. 北京: 科學(xué)出版社, 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.