(上海理工大學(xué) 能源與動(dòng)力工程學(xué)院; 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室，上海 200093)

巖心是礦物開采、石油勘探等地質(zhì)工程研究中極其重要的實(shí)物地質(zhì)資料。巖心分析是指針對(duì)采用各種儀器設(shè)備(如電鏡掃描儀、偏光與電子顯微鏡等)獲取的巖心圖像進(jìn)行觀測(cè)與模擬試驗(yàn)分析，在礦產(chǎn)開發(fā)與油氣地質(zhì)研究中具有重要作用。巖心顆粒圖像分割在巖心分析過程中占據(jù)著重要地位。巖心顆粒與背景區(qū)域的分割準(zhǔn)確性直接決定了顆粒目標(biāo)識(shí)別和粒度分布信息的提取的準(zhǔn)確性。巖心分析結(jié)果可作為油氣盆地地層和沉積環(huán)境的探索與開發(fā)的直接依據(jù)[1-2]。

國(guó)內(nèi)外學(xué)者提出了許多種巖心圖像分割算法和技術(shù)來(lái)獲取巖層物性特征、孔隙結(jié)構(gòu)和顆粒粒度分布等信息。Zhou等[3]結(jié)合形態(tài)學(xué)運(yùn)算的優(yōu)點(diǎn)，提出一種改進(jìn)的分水嶺算法，對(duì)巖心顆粒分割結(jié)果區(qū)域合并避免了過分割現(xiàn)象；Samet等[4]針對(duì)巖石圖像數(shù)據(jù)定義了一種新的模糊規(guī)則圖像分割技術(shù)，實(shí)現(xiàn)了顆粒邊緣檢測(cè)；陳本廷等[5]應(yīng)用超像素對(duì)巖心偏振圖像進(jìn)行初始化分割，再利用超像素區(qū)域的顏色信息進(jìn)行K均值聚類分割；Fazekas等[6]利用完全卷積神經(jīng)網(wǎng)絡(luò)(FCN)進(jìn)行巖心圖像分割，并基于深度學(xué)習(xí)的方法去除巖心圖像二次色彩的偽影。

雖然針對(duì)巖心圖像分割方法的研究已取得一定的成果，但大都在分割過程很少考慮到圖像的紋理特征。因此，本文中充分利用顆粒圖像中豐富的紋理、顏色和空間位置特征信息，針對(duì)巖心顆粒的顯微偏振圖像和電鏡掃描圖像進(jìn)行巖心顆粒的彩色圖像分割算法的研究，從而提高巖心顆粒的粒度分布信息提取的準(zhǔn)確性。

1 彩色圖像分割的原理與流程

在巖心分析過程中，需要從巖心顆粒彩色圖像中提取相關(guān)特征信息，確定圖像的像素點(diǎn)是否屬于同種特征或特征相似，進(jìn)而把圖像像素點(diǎn)分成不同的子集，判定像素集屬于巖心顆?；虮尘啊５珒H僅利用色彩信息對(duì)巖心顆粒彩色圖像進(jìn)行分割，往往會(huì)導(dǎo)致欠分割現(xiàn)象，巖心顆粒無(wú)法完全識(shí)別出來(lái)。因此，必須綜合利用局部紋理特性和空間位置關(guān)系，確保屬于同一顆粒的像素點(diǎn)的顏色、紋理特征均相似、空間位置上相鄰。

1.1 顏色特征空間的選取

顏色特征的描述依賴于所選取的色彩空間,針對(duì)圖像處理常用的色彩空間RGB、HSV、HSI和CIELAB等，研究學(xué)者做過大量對(duì)比性圖像分割試驗(yàn)，分析了在各類顏色空間上圖像分割的適用性[7-9]。

CIELAB色彩模型是一種基于生理特征的均勻色彩空間。CIELAB色彩空間上不同像素點(diǎn)的亮度信息分量為L(zhǎng)*，顏色特征信息分量分別為a*和b*，直接使用這3個(gè)分量之間的幾何距離來(lái)實(shí)現(xiàn)不同顏色之間的對(duì)比, 通過數(shù)字化形式直觀地描述了人眼的視覺感知，可以有效地應(yīng)用于測(cè)量微小色差。因此選取CIELAB色彩空間對(duì)巖心顆粒彩色圖像進(jìn)行分割研究。

在CIELAB色彩空間中，亮度信息L*對(duì)于顆粒識(shí)別具有一定的重要性，以特征向量[L*,a*,b*]的形式表示顏色信息可以通過亮度信息優(yōu)化那些僅利用顏色信息無(wú)法完全分割的顆粒圖像。它把顆粒圖像中每個(gè)像素點(diǎn)近似為三維空間上的一點(diǎn)，并將兩像素點(diǎn)之間的歐式距離ΔE作為顏色特征相似性測(cè)量判據(jù)：

(1)

1.2 二維Gabor濾波器原理

目前，紋理特征提取方法應(yīng)用廣泛的主要有灰度共生矩陣(GLCM)、局部二進(jìn)制模式(LBP)、高斯馬爾可夫隨機(jī)場(chǎng)(GMRF)和Gabor濾波變換等。其中，Gabor濾波器具有時(shí)域(或空域)和頻域的聯(lián)合最佳分辨率,能夠很好地模擬人類視覺感知系統(tǒng)的特點(diǎn)[10]，其最重要的優(yōu)點(diǎn)是對(duì)于旋轉(zhuǎn)、縮放和平移的不變性，而且對(duì)圖像噪聲干擾和光照條件改變具有很強(qiáng)的魯棒性[11-12],因而選用Gabor濾波器組在頻域上從不同尺度和方向提取巖心顆?；叶葓D像中的相關(guān)紋理特征。

二維Gabor濾波器函數(shù)實(shí)質(zhì)是二維高斯核函數(shù)和正弦平面波的乘積，其數(shù)學(xué)表達(dá)式[13]為

g(x,y,λ,θ,ψ,σ,γ)=

(2)

x′=xcosθ+ysinθ,

y′=-xsinθ+ycosθ,

(3)

式中：x和y分別為高斯核的坐標(biāo)信息；λ為正弦函數(shù)波長(zhǎng)；θ為Gabor核函數(shù)的方向；ψ為相位偏移；σ為高斯函數(shù)標(biāo)準(zhǔn)差；γ為空間寬高比。

將二維Gabor濾波器函數(shù)對(duì)應(yīng)的實(shí)數(shù)部分gR與虛數(shù)部分gI分別與巖心灰度圖像f(x，y)進(jìn)行卷積運(yùn)算，可得到濾波后Gabor幅值(能量)特征圖像S(x，y)的數(shù)學(xué)表達(dá)式為

(4)

1.3 Gabor紋理特征提取流程

在提取巖心顆粒圖像紋理特征時(shí)，需選取一個(gè)合適的Gabor濾波器組，在頻域的不同尺度、 不同方向上提取相關(guān)特征。 由于各種不同方向和尺度的Gabor濾波器是自相似的，因此可以從母小波函數(shù)g(x，y)經(jīng)過尺度變換和旋轉(zhuǎn)變換來(lái)構(gòu)建Gabor濾波器組。

選取母小波濾波函數(shù)的波長(zhǎng)初始值為4個(gè)像素值，Gabor核函數(shù)方向?yàn)?；然后母小波函數(shù)經(jīng)過尺度變換和旋轉(zhuǎn)變換構(gòu)建5個(gè)尺度、 8個(gè)方向的Gabor濾波器組，其中波長(zhǎng)范圍為4～16像素，方向?yàn)?～2π。 隨著卷積模板的增大，Gabor濾波器變換提取的特征逐漸由全局特征向局部特征變化，對(duì)比濾波結(jié)果選擇卷積模板大小為35×35像素。將巖心顆粒的彩色圖像從不同尺度和方向上利用Gabor濾波器組提取紋理特征的流程如圖1所示。

圖1 Gabor濾波器組提取紋理特征流程圖Fig.1 Flow chart of Gabor filter banks used to extract texture features

首先，將巖心顆粒的彩色圖像轉(zhuǎn)換為灰度圖像f(x，y), 彩色圖像像素大小為576×768，如圖2所示。

然后利用5個(gè)尺度、8個(gè)方向的Gabor濾波器組對(duì)巖心灰度圖像進(jìn)行卷積運(yùn)算，Gabor濾波器組輸出的40張幅值特征圖像如圖3所示。

最后，由于采用40個(gè)Gabor濾波器提取巖心顆粒圖像的紋理信息，對(duì)于圖2中大小為576×768像素的巖心灰度圖像，獲取的紋理特征向量維度為576×768×40，維數(shù)太大且存在許多不需要的特征信息，因此對(duì)圖3所示的多尺度、多方向Gabor濾波器組輸出濾波結(jié)果求均值，將得到的特征均值向量歸一化為零均值和單位方差，作為巖心顆粒圖像紋理特征T，圖4為Gabor濾波器組最終輸出的幅值特征均值圖像。

1.4 多維特征KFCM聚類巖心圖像分割流程

核模糊C均值(KFCM)[14]聚類算法中基于核函數(shù)思想，把巖心顆粒彩色圖像像素點(diǎn)特征信息映射至高維特征空間后，原線性不可分?jǐn)?shù)據(jù)信息變得可分或者接近于線性可分，改善了FCM聚類分割思想對(duì)于巖心圖像中非球形分布數(shù)據(jù)的性能，進(jìn)而實(shí)現(xiàn)巖心顆粒圖像數(shù)據(jù)更準(zhǔn)確的聚類分析。當(dāng)選用高斯核函數(shù)時(shí)，在高維特征空間中KFCM聚類的目標(biāo)函數(shù)[15]為

(5)

(6)

式中：c為聚類數(shù)目；n為像素點(diǎn)個(gè)數(shù)；m為加權(quán)指數(shù)；vj為第j個(gè)聚類中心；uij為像素點(diǎn)xi被分配在第j類的隸屬度；K為高斯核函數(shù)。

a 巖心顆粒原圖b 灰度圖像圖2 巖心顆粒的彩色圖像轉(zhuǎn)換為灰度圖像Fig.2 Color image of core particles converted into grayscale image

圖3 Gabor濾波器組輸出的40張幅值特征圖像Fig.3 40 value feature images output of Gabor filter bank

圖4 Gabor濾波器組輸出的幅值特征均值圖像Fig.4 Gabor filter bank output amplitude characteristic mean image

KFCM聚類的目標(biāo)函數(shù)經(jīng)拉格朗日乘子法最小化，推導(dǎo)得出算法迭代過程中的聚類中心vj和隸屬度uij為

(7)

(8)

綜合從巖心顆粒彩色圖像中提取的顏色信息(L*、a*、b*)、 Gabor紋理特征T和空間坐標(biāo)(Sx，Sy)3種特征信息，巖心顆粒彩色圖像多維特征KFCM聚類分割算法具體過程如下：

1)將巖心顆粒彩色圖像像素點(diǎn)包含的特征信息量化為特征向量，每一像素點(diǎn)表示為特征向量[L*,a*,b*,T,Sx,Sy]，作為樣本數(shù)據(jù)輸入；

2)參數(shù)初始化設(shè)置，設(shè)置算法的最大迭代次數(shù)和結(jié)束閾值ε，確定聚類數(shù)目c，加權(quán)指數(shù)m和高斯核尺度參數(shù)σ；

3)使用0～1之間的隨機(jī)數(shù)初始化隸屬度，且滿足式(6)的約束條件；

4)根據(jù)式(7)計(jì)算聚類中心vj；

5)由公式(5)獲得KFCM聚類分割的目標(biāo)函數(shù)值。如果此時(shí)超過最大迭代次數(shù)，或者相對(duì)上次迭代目標(biāo)函數(shù)值的變化量低于閾值ε，那么算法迭代過程結(jié)束；

6)否則，使用公式(8)計(jì)算新的隸屬度函數(shù)并返回步驟4繼續(xù)執(zhí)行迭代。

巖心顆粒圖像的多維特征KFCM聚類分割算法流程如圖5所示。

按照?qǐng)D5所示的圖像分割流程，在巖心顆粒的顯微偏振圖像和掃描儀采集圖像上進(jìn)行分割實(shí)驗(yàn)。分割過程中遵循同一聚類內(nèi)相似度盡可能大、不同聚類之間相似度盡可能小的原則，比較在不同聚類數(shù)目下目標(biāo)損失函數(shù)，以此衡量聚類分割后的結(jié)果，從而確定適用于當(dāng)前類型的巖心顆粒彩色圖像的最佳聚類數(shù)目。

圖5 巖心顆粒圖像的多維特征KFCM聚類分割算法流程Fig.5 Multi-dimensional feature KFCM clustering segmentation algorithm for core particle image

2 仿真試驗(yàn)與結(jié)果

2.1 試驗(yàn)?zāi)康暮蜅l件

為驗(yàn)證算法可行性，首先對(duì)比分析多維特征KFCM聚類分割算法、超像素優(yōu)化的基于顏色信息的K均值聚類分割方法[5]和基于HSV空間V分量分割方法[1]這3種巖心圖像分割算法的分割效果；然后，應(yīng)用多維特征KFCM聚類分割算法在不同類型的巖心顆粒彩色圖像上進(jìn)行分割試驗(yàn)。不同類型的巖心顆粒彩色圖像均利用偏光顯微鏡、巖心掃描儀采集而得。

試驗(yàn)的硬件環(huán)境條件為：Intel(R) Core(TM) i5-3230M型CPU，內(nèi)存4G；軟件環(huán)境為：Window10操作系統(tǒng)，Matlab R2017a。

2.2 不同算法分割效果的對(duì)比

分別應(yīng)用多維特征KFCM聚類分割算法、超像素優(yōu)化的基于顏色信息的K均值聚類分割方法和基于HSV空間V分量分割方法，在圖2所示的巖心顆粒彩色圖像上進(jìn)行分割試驗(yàn)，不同算法類型的分割結(jié)果如圖6所示。

a KFCM聚類分割算法b K均值聚類分割算法c V分量分割算法圖6 不同算法的彩色圖像分割效果Fig.6 Results of different algorithm types of color image segmentation

多維特征KFCM聚類分割算法和超像素優(yōu)化的K均值聚類分割算法均是基于聚類分析思想，圖6b所示超像素優(yōu)化的K均值聚類分割方法首先將巖心圖像進(jìn)行超像素預(yù)分割，選定超像素個(gè)數(shù)為10 000個(gè)；然后在分割過程中此兩種算法聚類數(shù)c均為8，最大迭代次數(shù)為100次，聚類結(jié)果根據(jù)聚類數(shù)目分層后，選取相應(yīng)的顆粒層合并獲得最終分割結(jié)果。對(duì)比圖6a、6b可知，K均值聚類分割算法經(jīng)超像素預(yù)分割后僅依賴于顏色特征信息，當(dāng)背景顏色與顆粒像素點(diǎn)的顏色接近時(shí)會(huì)被分配至同一個(gè)聚類中，存在欠分割現(xiàn)象，顆粒間發(fā)生一定程度的粘連，白線框內(nèi)為一些粘連或欠分割顆粒。而融合紋理、顏色與空間信息的KFCM聚類分割算法增強(qiáng)了顆粒識(shí)別的穩(wěn)定性，分割結(jié)果更準(zhǔn)確。

如圖6c所示，V分量分割算法首先將原始巖心圖像轉(zhuǎn)換至HSV顏色空間下，利用亮度信息V分量進(jìn)行閾值分割，V分量迭代最佳閾值為0.630 69；分割后部分顆粒失真，存在一定孔洞現(xiàn)象，紅線框中為一些未能完整識(shí)別的巖心顆粒，其原因主要是由于亮度相同而顏色不同的像素點(diǎn)無(wú)法正確識(shí)別。

綜上，在3種算法中，多維特征KFCM聚類分割算法的分割結(jié)果最優(yōu)。

2.3 不同類型巖心顆粒彩色圖像的分割效果

圖7為2種典型的巖心顆粒顯微偏振圖像，由于顆粒樣本的差異，圖像的質(zhì)量存在差異，特別是顆粒與背景之間的顏色和對(duì)比度差別。

圖7a、7b分割過程中聚類數(shù)c分別為8和6，巖心顯微偏振圖像轉(zhuǎn)換成灰度圖像后，應(yīng)用Gabor濾波器組從巖心灰度圖像提取的紋理特征均值圖像如圖8所示。

a 巖心原圖1b 巖心原圖2圖7 2種巖心顆粒顯微偏振圖像Fig.7 Microscopic polarization image of core particles

a Gabor紋理圖像1b Gabor紋理圖像2圖8 巖心顆粒顯微偏振圖像的Gabor紋理均值圖像Fig.8 Gabor texture mean image of microscopic polarization image of core particles

融合從巖心顯微偏振圖像中提取的顏色、空間坐標(biāo)以及Gabor紋理特征，進(jìn)行KFCM均值聚類，然后根據(jù)聚類數(shù)目分層后選取相應(yīng)的顆粒層合并，最終分割結(jié)果如圖9所示。

同理，遵照上述的多維特征KFCM聚類分割算法的過程，圖10為從巖心掃描儀中采集的巖心顆粒掃描圖像的分割過程，分割過程中聚類數(shù)為5。

針對(duì)巖心顆粒的顯微偏振圖像和掃描圖像進(jìn)行多維特征核模糊C均值聚類分割試驗(yàn)，均獲得較好的分割結(jié)果，能清晰地識(shí)別每個(gè)巖心顆粒的整體輪廓。 圖7a、 7b中大部分巖心顆粒著色均勻，不同顆粒以及顆粒與背景之間的顏色對(duì)比度大，顆粒分割效果最好，小粒徑顆粒均能被識(shí)別。 圖10a巖心掃描圖像雖然存在噪聲影響，巖心顆粒與背景顏色相接近，但也能較好地將顆粒與背景分離。

3 結(jié)論與展望

將巖心顆粒彩色圖像包含的顏色、空間位置和紋理特征信息融合后進(jìn)行多維特征KFCM均值聚類分割，能更準(zhǔn)確地實(shí)現(xiàn)巖心顆粒彩色圖像的分割，獲得可靠的顆粒粒度分布信息。

多維特征KFCM均值聚類數(shù)目和隨機(jī)的初始化聚類中心的確定方法，對(duì)分割結(jié)果的穩(wěn)定性存在一定的影響，這一問題可作為后期進(jìn)一步研究的方向。

a 分割結(jié)果1b 分割結(jié)果2圖9 巖心顆粒顯微偏振圖像最終分割結(jié)果Fig.9 Final segmentation results of microscopic polarization image of core particles

a 巖心掃描圖像b Gabor紋理特征均值圖像c 分割結(jié)果圖10 巖心掃描圖像分割過程Fig.10 Segmentation process of core scan images