利用多重-大津閾值算法和掃描電鏡分割CT圖像

高衍武，吳偉，張虔，趙燕紅，邵廣輝，李國(guó)利，毛超杰

(中國(guó)石油集團(tuán)測(cè)井有限公司測(cè)井應(yīng)用研究院，陜西 西安 710000)(中石油玉門油田分公司勘探開(kāi)發(fā)研究院，甘肅 酒泉 735000)(中國(guó)石油集團(tuán)測(cè)井有限公司測(cè)井應(yīng)用研究院，陜西 西安 710000)(中石油玉門油田分公司勘探開(kāi)發(fā)研究院，甘肅 酒泉 735000)

頁(yè)巖復(fù)雜的孔隙分布導(dǎo)致儲(chǔ)層非均質(zhì)性強(qiáng)，孔隙度、滲透率、孔隙大小分布是頁(yè)巖孔隙結(jié)構(gòu)評(píng)價(jià)中重要的參數(shù)[1～4]。作為典型的致密儲(chǔ)層，有必要評(píng)價(jià)頁(yè)巖孔隙結(jié)構(gòu)特征，以便更好地理解頁(yè)巖氣體的儲(chǔ)集、運(yùn)移機(jī)理[5，6]。評(píng)價(jià)頁(yè)巖孔隙結(jié)構(gòu)方法主要有3種：高壓壓汞法、氣體吸附法、場(chǎng)發(fā)射環(huán)境掃描電鏡法(SEM)[7，8]。上述方法可以提供一些局部信息，如孔隙類型、孔隙大小和比表面積等，但都不能從三維尺度評(píng)價(jià)頁(yè)巖孔隙結(jié)構(gòu)，且均對(duì)樣品有破壞作用。作為一種不破壞樣品的測(cè)試方法，XCT(X-ray computed tomography)提供了一種高效的、定量評(píng)價(jià)頁(yè)巖三維孔隙結(jié)構(gòu)的方法，其主要由圖像獲取、圖像處理、結(jié)構(gòu)探測(cè)3個(gè)階段組成。在圖像獲取階段，分析圖像特征來(lái)獲取孔隙結(jié)構(gòu)特征，其主要由圖像分割和孔隙-裂縫提取2部分組成。圖像分割是圖像處理中最重要也最具有挑戰(zhàn)意義的一步[9]，雙掃描算法[10，11]和閾值算法[12～19]都可以用于孔隙結(jié)構(gòu)分割。閾值通常作為圖像分割體的標(biāo)準(zhǔn)來(lái)提取背景圖像中的孔隙結(jié)構(gòu)[20]。閾值技術(shù)又分為全閾值方法和局部閾值方法[21]，局部閾值方法對(duì)灰色信息局部化來(lái)選取多重閾值，該方法處理過(guò)程慢、復(fù)雜；全閾值方法對(duì)整個(gè)圖像使用單一閾值，因此廣泛用于CT(computed tomography)圖像分割中。

多種閾值方法用來(lái)提取CT圖像中的孔隙-裂縫[16,22,19]，如大津閾值算法[14，23]，但對(duì)多礦物組分CT圖像方面的研究很少。利用數(shù)字巖心對(duì)多礦物組分CT圖像進(jìn)行研究，對(duì)于理解不同礦物組分對(duì)流體滲流的影響有重要作用。多礦物組分的CT圖像和掃描電鏡(SEM)圖像也可以通過(guò)人工觀測(cè)的方法分割[24]，但是人工觀察的方法具有耗時(shí)、易受到人為主觀因素影響等缺點(diǎn)。大津閾值算法[25]對(duì)圖像直方圖上的最大類間方差或最小類間方差選擇閾值，是一種對(duì)圖像選擇閾值最常用的方法[26]。分割圖像上的多重元素，則利用多重-大津閾值算法[27]?；诙嘀?大津閾值算法，筆者提出了一種多重閾值算法聯(lián)合場(chǎng)發(fā)射環(huán)境SEM的方法，用以分割CT圖像中的多重元素。

1 樣品和方法

1.1 樣品

研究樣品主要取自2口井，S1樣品取自一口井，S2、S3樣品來(lái)自另一口井。研究區(qū)位于渤海灣盆地東南部的東營(yíng)凹陷，頁(yè)巖樣品深度在2500～3000m，位于古近系沙河街組，是東營(yíng)凹陷頁(yè)巖油勘探開(kāi)發(fā)的主要層位。3塊樣品都切削為直徑25mm的圓柱狀，利用氦氣擴(kuò)散法測(cè)量干巖樣在實(shí)驗(yàn)室條件下的有效孔隙度，同時(shí)對(duì)樣品碎屑進(jìn)行了XRD(X-ray diffraction)礦物組分分析和總有機(jī)碳質(zhì)量分?jǐn)?shù)(w(TOC))測(cè)試分析(見(jiàn)表1)。該次研究中S1、S2樣品利用多重-大津閾值算法和SEM進(jìn)行CT圖像分割，S3樣品用于驗(yàn)證該方法。

表1 3塊樣品礦物組分質(zhì)量分?jǐn)?shù)、w(TOC)及有效孔隙度統(tǒng)計(jì)

1.2 試驗(yàn)

使用美國(guó)Xradia公司制造的Micro XCT-200儀器進(jìn)行XCT測(cè)量，其對(duì)長(zhǎng)度1mm物體的掃描精度可達(dá)0.7μm。圓柱狀頁(yè)巖樣品垂直于樣品臺(tái)放置在掃描儀的中央位置。掃描樣品的X射線源電壓為150kV、功率10W。S1號(hào)樣品長(zhǎng)度5.2cm，直徑2.512cm；S2號(hào)樣品長(zhǎng)度6.1cm，直徑2.512cm；S3號(hào)樣品長(zhǎng)度4.8cm，直徑2.512cm。樣品的重構(gòu)切片圖像分辨率為2048×2048，切片長(zhǎng)度為13.7135μm，產(chǎn)生的立體像素的體積為2578.96μm3?？傮w積為16.8275cm3的3個(gè)頁(yè)巖樣品在該試驗(yàn)條件下掃描產(chǎn)生1571個(gè)二維切片圖像，并對(duì)二維切片圖像疊加產(chǎn)生三維數(shù)據(jù)。

為了獲取頁(yè)巖樣品礦物組分的灰度成分分布，對(duì)S1、S2樣品進(jìn)行SEM測(cè)試分析。在SEM觀察之前，2種頁(yè)巖樣品巖心碎屑首先進(jìn)行氬離子拋光，根據(jù)GB/T 5162—1997標(biāo)準(zhǔn)對(duì)巖心碎屑進(jìn)行涂碳，目的是提高巖心表面的導(dǎo)電性進(jìn)而增加圖像的質(zhì)量和分辨率。利用FEI Quanta 200F場(chǎng)發(fā)射環(huán)境掃描電鏡進(jìn)行觀測(cè)，圖像被放大8000倍，觀測(cè)分辨率可達(dá)10nm，通過(guò)搭載能譜儀的X射線光譜儀定量分析頁(yè)巖礦物組分。

1.3 方法

1.3.1 CT和背散射電子(BSE)成像技術(shù)原理

X射線照射到不同礦物時(shí)，由于光電效應(yīng)、康普頓效應(yīng)、電子對(duì)效應(yīng)等復(fù)雜物理過(guò)程的影響，造成其信號(hào)衰減幅度不同[18,28]，可以用Lamber-Beer定律來(lái)表示：

I=I0e-μΔx

(1)

(2)

式中：I為傳輸強(qiáng)度，cd；I0為離子束發(fā)射強(qiáng)度,cd；μ為線性衰減幅度系數(shù),1；Δx為X射線照射的樣品厚度,cm；ρb為物體的骨架密度,g/cm3；a為Kling-Nishian系數(shù),1；b取值9.8×10-24；Z是原子數(shù)目,1；E為X射線能量,kV。

在低X射線能量(<100kV)下，μ是一個(gè)與Z有主要關(guān)系的函數(shù)；在高X射線能量(>100kV)下，μ主要是與E有關(guān)的函數(shù)[29]。該次研究中X射線電壓設(shè)置在150kV，因此X射線衰減幅度主要取決于物體的密度。不同的物質(zhì)由于其密度不同產(chǎn)生不同衰減幅度的信號(hào)。CT值是用來(lái)定量評(píng)價(jià)X射線穿透物質(zhì)時(shí)的衰減幅度，在高密度礦物中，X射線衰減量大，CT值也相應(yīng)變大。CT值計(jì)算公式為：

(3)

式中：μw為水的衰減幅度系數(shù),1。

該次研究的初始CT切片圖像被轉(zhuǎn)換成具有256位灰色范圍值的灰度圖像，在該范圍中0代表最黑色，255代表最亮色。CT圖像像素的灰度值主要與巖石樣品密度有關(guān)：對(duì)于高密度物體，其圖像為亮色；對(duì)于低密度物體，圖像為暗色。

場(chǎng)發(fā)射環(huán)境SEM的原理是利用高能量電子轟擊樣品表面積產(chǎn)生次生電子和背散射電子，次生電子對(duì)樣品表面敏感，可以有效揭示其微孔圖像特征，因此被廣泛用于孔隙研究。背散射電子是指從樣品表面反射角度小于90°的原子束，其反射強(qiáng)度依賴于樣品密度，物體的密度越高，產(chǎn)生的背散射電子反射強(qiáng)度越高，其生成的背散射電子圖像就越亮。背散射電子圖像像素的灰度值與巖石樣品密度有關(guān)，因此可以用來(lái)印證CT圖像不同礦物組分灰度值的分布。

1.3.2 多重-大津閾值算法

大津閾值算法是對(duì)圖像選擇閾值的一種最常用方法[25，26]，其核心思想是發(fā)現(xiàn)最大類間方差的閾值。一個(gè)圖像可以描述為I(x，y)，其灰度值為0～L-1，L為灰度區(qū)間最大值，灰度值為I的像素值為ni，灰色圖像中像素總數(shù)為n。該灰度值出現(xiàn)的概率Pi為：

(4)

整個(gè)圖像的平均灰度范圍μT為：

(5)

當(dāng)一個(gè)單一的閾值t使用時(shí)，一個(gè)圖像被分成2個(gè)部分——D0和D1，D0像素的灰度區(qū)間為[0，t]，D1像素的灰度區(qū)間為[t+1，L-1]。P0(t)、P1(t)表示D0和D1的累加概率，μo(t)和μ1(t)表示D0和D1的灰度平均值:

(6)

(7)

(8)

(9)

對(duì)于閾值t，其類間方差σB2(t)定義為：

σB2(t) =P0(t)(μ0(t)-μT)2+P1(t)(μ1(t)-μT)2

(10)

最優(yōu)化閾值T為：

T= argmaxσB2(t)

(11)

當(dāng)一個(gè)圖像被分割成k+1部分時(shí)，k個(gè)不同的閾值(t1，t2，t3，…，tk)被使用，變量k的不同閾值的類間方差為：

(12)

多重-大津閾值算法中最優(yōu)化閾值為：

(13)

1.3.3 CT圖像和SEM圖像處理

使用ImageJ數(shù)字圖像處理分析軟件對(duì)單一的BSE圖像和CT圖像進(jìn)行分析，利用BSE圖像中灰度值模式識(shí)別不同的元素相(包括孔隙-裂縫、有機(jī)質(zhì)、礦物)，識(shí)別結(jié)果與XRD、EDX(搭載能譜儀的X射線光譜儀)測(cè)試結(jié)果進(jìn)行相互驗(yàn)證。根據(jù)直方圖法求取不同元素相的平均灰度值，利用多重-大津閾值算法求取CT切片圖像中不同元素相的閾值；CT圖像被分割之后，在ImageJ軟件中分析所有元素相，分割的圖像在ImageJ軟件中導(dǎo)出后疊加用于重構(gòu)三維數(shù)字巖心。

1.4 礦物和有機(jī)質(zhì)體積分?jǐn)?shù)

通過(guò)XRD分析測(cè)試得到的礦物組分含量是在沒(méi)有有機(jī)質(zhì)和流體的存在下礦物的質(zhì)量分?jǐn)?shù)，但是w(TOC)作為有機(jī)質(zhì)在巖石中占據(jù)一定的比例，通過(guò)CT圖像疊加分析計(jì)算的礦物組分和有機(jī)質(zhì)含量是其占巖石總體積的體積分?jǐn)?shù)，因此礦物組分和有機(jī)質(zhì)含量應(yīng)轉(zhuǎn)換為體積分?jǐn)?shù)，然后歸一化處理。

為了消除有機(jī)質(zhì)含量的影響，首先對(duì)礦物組分的質(zhì)量分?jǐn)?shù)進(jìn)行歸一化處理：

wnm,a=wm,a×(100-wom)

(14)

式中：wnm,a為礦物歸一化后的質(zhì)量分?jǐn)?shù),%；wm,a為XRD分析測(cè)試的礦物質(zhì)量分?jǐn)?shù),%；wom為有機(jī)質(zhì)質(zhì)量分?jǐn)?shù),%;下標(biāo)a為不同種類的礦物。

礦物組分和有機(jī)質(zhì)含量根據(jù)其密度歸一化后的體積分?jǐn)?shù)為：

(15)

(16)

式中：Vnm,a、Vnom分別為歸一化后礦物組分和有機(jī)質(zhì)的體積分?jǐn)?shù),%；ρa(bǔ)為礦物組分的密度,g/cm3；ρom為有機(jī)質(zhì)的密度,g/cm3；φ為孔隙度,%。

礦物組分、有機(jī)質(zhì)、孔隙-裂縫歸一化后的體積分?jǐn)?shù)見(jiàn)表2。

表2 3塊樣品礦物組分、有機(jī)質(zhì)及孔隙-裂縫歸一化體積分?jǐn)?shù)統(tǒng)計(jì)

2 結(jié)果和討論

2.1 礦物成分鑒定及其灰度值

在BSE圖像中，黃鐵礦(最亮的多邊形顆粒)、有機(jī)質(zhì)(最暗的無(wú)定形顆粒)、孔隙-裂縫由于其自身特別高或低的密度被識(shí)別出來(lái)，其他的礦物組分,如石英、碳酸鹽礦物、長(zhǎng)石、菱鐵礦、黏土礦物等由于其具有相同的元素組分，不能僅僅根據(jù)其自身的灰度值而被區(qū)分出來(lái)，圖像中特別暗或亮的區(qū)域可能是由于部分礦物顆粒導(dǎo)電性較弱或表面不規(guī)則性而導(dǎo)致局部電子電荷釋放造成的。

根據(jù)BSE圖像、EDX圖像和XRD分析測(cè)試結(jié)果確定其他礦物，每一種礦物和有機(jī)質(zhì)根據(jù)其自身形狀和灰度值來(lái)識(shí)別(見(jiàn)圖1，其中(a)～(c)為S1樣品；(d)～(f)為S2樣品)。最亮的多邊形顆粒為黃鐵礦，在BSE圖像中呈草莓狀(見(jiàn)圖1(c)～(e))；有機(jī)質(zhì)和孔隙-裂縫在BSE圖像中以無(wú)定形態(tài)呈黑色區(qū)域(見(jiàn)圖1(b)～(f))；菱鐵礦在BSE圖像中呈最亮的平行四邊形，由于其存在Fe元素，利用EDX可以將其區(qū)分出來(lái)(見(jiàn)圖1(c)、(e)、(f))；碳酸鹽礦物大多為方解石，在BSE圖像中呈較亮的灰色不規(guī)則礦物顆粒(見(jiàn)圖1(a)、(b)、(d)～(f))；石英呈黑色圓形顆粒(見(jiàn)圖1(b)、(d)、(e))；正長(zhǎng)石呈亮的灰色圓狀顆粒(見(jiàn)圖1(a)、(e))；斜長(zhǎng)石在EDX圖像中呈暗的灰色圓狀顆粒(見(jiàn)圖1(a)、(f))；黏土礦物通常為暗的灰色不規(guī)則長(zhǎng)條狀的集合體，其內(nèi)部發(fā)育微構(gòu)造(見(jiàn)圖1(a)、(d))。

圖1 S1、S2樣品在BSE模式下的礦物組分、有機(jī)質(zhì)和孔隙-裂縫及其灰度值

根據(jù)圖1中每種礦物組分的平均灰度值，將7種礦物組分、有機(jī)質(zhì)和孔隙-裂縫分成6種不同元素部分——有機(jī)質(zhì)和孔隙-裂縫、石英、黏土礦物和斜長(zhǎng)石、碳酸鹽礦物和正長(zhǎng)石、菱鐵礦、黃鐵礦。S1樣品6種不同元素部分根據(jù)其平均灰度值依降序排列為：黃鐵礦(平均灰度值254)>菱鐵礦(平均灰度值237)>碳酸鹽礦物和正長(zhǎng)石(平均灰度值221～233)>石英(平均灰度值197～199)>黏土礦物和斜長(zhǎng)石(平均灰度值191～193)>有機(jī)質(zhì)和孔隙-裂縫(平均灰度值<78)。S2樣品6種不同元素部分根據(jù)其平均灰度值依降序排列為：黃鐵礦(平均灰度值為253～255)>菱鐵礦(平均灰度值230～235)>碳酸鹽礦物和正長(zhǎng)石(平均灰度值190～199)>黏土礦物和斜長(zhǎng)石(平均灰度值165～167)>石英(平均灰度值160～161)>有機(jī)質(zhì)和孔隙-裂縫(平均灰度值<56)。

2.2 CT圖像分割

在原始CT切片圖像中央提取像素為1300×1300的局部圖像以消除巖心外部區(qū)域造成的誤差，為了區(qū)分6種不同元素部分，確定出5種不同的全閾值(T1、T2、T3、T4、T5)。根據(jù)大津閾值算法求取了5種全閾值并獲取了分割圖像(見(jiàn)圖2)，為了獲取5個(gè)全閾值，對(duì)每50個(gè)切片圖像進(jìn)行多閾值分割，每一個(gè)CT圖像被分割為30個(gè)切片圖像，最終通過(guò)30個(gè)切片圖像中的平均閾值求取CT圖像中每個(gè)元素的全閾值(見(jiàn)表3)。每個(gè)元素的閾值可以通過(guò)BSE圖像中灰度值分布圖來(lái)確定。如S1樣品，黃鐵礦、菱鐵礦、碳酸鹽礦物和正長(zhǎng)石、石英、黏土礦物和斜長(zhǎng)石、有機(jī)質(zhì)和孔隙-裂縫的閾值分別為58～255、47～57、45～46、43～44、41～42、0～40。

根據(jù)上述閾值，利用ImageJ軟件計(jì)算出6種不同元素部分所占圖像的比例，對(duì)單一切片進(jìn)行處理產(chǎn)生平面含量，最后對(duì)所有切片的平均平面含量進(jìn)行計(jì)算，產(chǎn)生6種不同元素部分的體積分?jǐn)?shù)。在原始圖像中提取局部圖像后，總共有1571張像素為1300×1300的灰色圖像被用來(lái)作為樣品。在灰色圖像中，每個(gè)像素的灰度值范圍為0～255，每個(gè)元素部分有一定的灰度區(qū)間值，利用Boolean方法計(jì)算每個(gè)元素部分的灰色值[19]：

(17)

式中：g(i,j)為圖像的二進(jìn)制值；f(i,j)為坐標(biāo)(i,j)的像素灰度值；Ti-1、Ti為每一元素部分的閾值區(qū)間范圍。

圖2 CT切片圖像切割過(guò)程

表3 CT圖像中每個(gè)元素部分的全閾值

根據(jù)Boolean 方法，CT圖像被轉(zhuǎn)換為二進(jìn)制的黑白圖像，黑色圖像(g(i,j)=0)代表需要研究的部分(region of interest，ROI)，通過(guò)計(jì)算黑色區(qū)域的像素值和圖像中的像素總和進(jìn)一步求取黑色區(qū)域平面含量：

(18)

式中：Ci為切片中每一元素部分的平面含量，%；NROI為黑色區(qū)域像素值，pt；NT為切片中所有像素總和,pt。

利用ImageJ離子分析工具分別計(jì)算出S1、S2樣品的6個(gè)不同元素部分的體積分?jǐn)?shù)(見(jiàn)圖3(a)、(b))，與XRD、氣測(cè)孔隙度等實(shí)際測(cè)量結(jié)果相比可以看出,CT計(jì)算值與CT實(shí)際測(cè)量值相關(guān)性較好，說(shuō)明多重-大津閾值算法求取的閾值可以準(zhǔn)確進(jìn)行CT圖像分割。

圖3 頁(yè)巖樣品CT計(jì)算值與實(shí)際測(cè)量值關(guān)系圖

3 S3樣品驗(yàn)證

S3樣品用于驗(yàn)證聯(lián)合多重-大津閾值算法、多元素分割方法、BSE圖像等方法對(duì)CT圖像進(jìn)行處理的結(jié)果。S2、S3樣品取自同一口井、相同層位、有相同的物源和礦物組分，因此S3樣品與S2樣品的不同元素灰度值分布圖類似。利用上述方法求取了S3樣品的黃鐵礦、菱鐵礦、碳酸鹽礦物和正長(zhǎng)石、黏土礦物和斜長(zhǎng)石、石英、有機(jī)質(zhì)和孔隙-裂縫的閾值和體積分?jǐn)?shù)(見(jiàn)表3，圖3(c))。由圖3(c)可以看出，S3樣品的CT計(jì)算值與CT實(shí)際測(cè)量值有較好的正相關(guān)性(R2=0.90)，說(shuō)明相同物源和礦物組分的樣品在CT圖像不同元素部分中有相同的灰度值分布，同時(shí)說(shuō)明該分割方法可以用來(lái)準(zhǔn)確、定量評(píng)價(jià)CT圖像。

4 三維模型的重構(gòu)

為了獲取3塊樣品中6種不同元素部分在三維狀態(tài)下的分布特征，對(duì)CT圖像進(jìn)行疊加重構(gòu)三維模型?？紤]到計(jì)算機(jī)運(yùn)行處理速度，在疊加的CT圖像中提取像素500×500×500的立體模型，根據(jù)6個(gè)不同元素部分的閾值把CT切片圖像轉(zhuǎn)換為分割圖像，在ImageJ軟件中處理分割圖像產(chǎn)生不同元素部分的三維模型(見(jiàn)圖4)。分割圖像重構(gòu)的三維模型可以直觀地展現(xiàn)不同元素的分布特征和非均質(zhì)性特征。由圖4可以看出，S1樣品有較小的非均質(zhì)性，孔隙、有機(jī)質(zhì)、黏土礦物在空間分布上有明顯的非均質(zhì)性弱的特點(diǎn)(見(jiàn)圖4(a))；S2、S3樣品有較強(qiáng)的非均質(zhì)性，孔隙-裂縫、有機(jī)質(zhì)、黏土礦物分布范圍差異較大(見(jiàn)圖4(b)、(c))。

圖4 6種不同元素部分的三維重構(gòu)模型

CT圖像的多元素分割方法可以準(zhǔn)確地確定6種不同元素部分的灰度值分布和閾值，該方法比人工觀測(cè)更高效，能夠消除人為因素的影響，但是也存在部分缺陷：不能有效區(qū)分CT圖像中每一種礦物組分，特別是孔隙-裂縫與有機(jī)質(zhì)，因?yàn)閮烧叩幕叶戎捣植枷嘀丿B[13]；每一個(gè)CT圖像中，像素的灰度值代表了表面物體體積的衰減特征，當(dāng)物體體積由多種不同物質(zhì)組成時(shí)，其灰度值代表了平均衰減特征[19]，因此，部分孔隙體積效應(yīng)對(duì)頁(yè)巖礦物組分的鑒別帶來(lái)困難。

5 結(jié)語(yǔ)

筆者根據(jù)BSE圖像求取每個(gè)礦物組分灰度值分布范圍，確定出每個(gè)礦物的閾值，然后根據(jù)多重-大津閾值算法分割CT圖像。利用多組分分割方法疊加圖像可有效切割CT圖像，S1、S2樣品計(jì)算的CT值與實(shí)際測(cè)量值相符。利用S3樣品的CT圖像驗(yàn)證該方法，結(jié)果表明：有相同物源和礦物組分的樣品在BSE和CT圖像中有相同的灰度值分布，利用該方法可以準(zhǔn)確、高效地切割CT圖像。