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

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

(中國石油集團測井有限公司測井應用研究院，陜西 西安 710000)(中石油玉門油田分公司勘探開發研究院，甘肅 酒泉 735000)(中國石油集團測井有限公司測井應用研究院，陜西 西安 710000)(中石油玉門油田分公司勘探開發研究院，甘肅 酒泉 735000)

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

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

1 樣品和方法

1.1 樣品

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

表1 3塊樣品礦物組分質量分數、w(TOC)及有效孔隙度統計

1.2 試驗

使用美國Xradia公司制造的Micro XCT-200儀器進行XCT測量，其對長度1mm物體的掃描精度可達0.7μm。圓柱狀頁巖樣品垂直于樣品臺放置在掃描儀的中央位置。掃描樣品的X射線源電壓為150kV、功率10W。S1號樣品長度5.2cm，直徑2.512cm；S2號樣品長度6.1cm，直徑2.512cm；S3號樣品長度4.8cm，直徑2.512cm。樣品的重構切片圖像分辨率為2048×2048，切片長度為13.7135μm，產生的立體像素的體積為2578.96μm3。總體積為16.8275cm3的3個頁巖樣品在該試驗條件下掃描產生1571個二維切片圖像，并對二維切片圖像疊加產生三維數據。

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

1.3 方法

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

X射線照射到不同礦物時，由于光電效應、康普頓效應、電子對效應等復雜物理過程的影響，造成其信號衰減幅度不同[18,28]，可以用Lamber-Beer定律來表示：

I=I0e-μΔx

(1)

(2)

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

在低X射線能量(<100kV)下，μ是一個與Z有主要關系的函數；在高X射線能量(>100kV)下，μ主要是與E有關的函數[29]。該次研究中X射線電壓設置在150kV，因此X射線衰減幅度主要取決于物體的密度。不同的物質由于其密度不同產生不同衰減幅度的信號。CT值是用來定量評價X射線穿透物質時的衰減幅度，在高密度礦物中，X射線衰減量大，CT值也相應變大。CT值計算公式為：

(3)

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

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

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

1.3.2 多重-大津閾值算法

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

(4)

整個圖像的平均灰度范圍μT為：

(5)

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

(6)

(7)

(8)

(9)

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

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

(10)

最優化閾值T為：

T= argmaxσB2(t)

(11)

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

(12)

多重-大津閾值算法中最優化閾值為：

(13)

1.3.3 CT圖像和SEM圖像處理

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

1.4 礦物和有機質體積分數

通過XRD分析測試得到的礦物組分含量是在沒有有機質和流體的存在下礦物的質量分數，但是w(TOC)作為有機質在巖石中占據一定的比例，通過CT圖像疊加分析計算的礦物組分和有機質含量是其占巖石總體積的體積分數，因此礦物組分和有機質含量應轉換為體積分數，然后歸一化處理。

為了消除有機質含量的影響，首先對礦物組分的質量分數進行歸一化處理：

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

(14)

式中：wnm,a為礦物歸一化后的質量分數,%；wm,a為XRD分析測試的礦物質量分數,%；wom為有機質質量分數,%;下標a為不同種類的礦物。

礦物組分和有機質含量根據其密度歸一化后的體積分數為：

(15)

(16)

式中：Vnm,a、Vnom分別為歸一化后礦物組分和有機質的體積分數,%；ρa為礦物組分的密度,g/cm3；ρom為有機質的密度,g/cm3；φ為孔隙度,%。

礦物組分、有機質、孔隙-裂縫歸一化后的體積分數見表2。

表2 3塊樣品礦物組分、有機質及孔隙-裂縫歸一化體積分數統計

2 結果和討論

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

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

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

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

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

2.2 CT圖像分割

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

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

(17)

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

圖2 CT切片圖像切割過程

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

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

(18)

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

利用ImageJ離子分析工具分別計算出S1、S2樣品的6個不同元素部分的體積分數(見圖3(a)、(b))，與XRD、氣測孔隙度等實際測量結果相比可以看出,CT計算值與CT實際測量值相關性較好，說明多重-大津閾值算法求取的閾值可以準確進行CT圖像分割。

圖3 頁巖樣品CT計算值與實際測量值關系圖

3 S3樣品驗證

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

4 三維模型的重構

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

圖4 6種不同元素部分的三維重構模型

CT圖像的多元素分割方法可以準確地確定6種不同元素部分的灰度值分布和閾值，該方法比人工觀測更高效，能夠消除人為因素的影響，但是也存在部分缺陷：不能有效區分CT圖像中每一種礦物組分，特別是孔隙-裂縫與有機質，因為兩者的灰度值分布相重疊[13]；每一個CT圖像中，像素的灰度值代表了表面物體體積的衰減特征，當物體體積由多種不同物質組成時，其灰度值代表了平均衰減特征[19]，因此，部分孔隙體積效應對頁巖礦物組分的鑒別帶來困難。

5 結語

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