張吉群, 胡長(zhǎng)軍, 和冬梅, 常軍華, 李心浩, 李華

(1.北京科技大學(xué), 北京 100083; 2.中國(guó)石油規(guī)劃總院, 北京 100083; 3.中國(guó)石油勘探開(kāi)發(fā)研究院, 北京 100083)

0 引 言

國(guó)內(nèi)外研究孔隙結(jié)構(gòu)的方法主要分為3類(lèi)。①室內(nèi)實(shí)驗(yàn)法。該方法主要包括毛細(xì)管力曲線法、鑄體薄片法、掃描電鏡法、CT掃描法,這4種方法如果不借助第3類(lèi)方法,一般都是采用人工鑒定,分析速度慢,準(zhǔn)確度取決于個(gè)人鑒定水平。②測(cè)井資料現(xiàn)場(chǎng)評(píng)價(jià)法。該方法主要是用電阻率測(cè)井資料或核磁共振測(cè)井研究巖石孔隙結(jié)構(gòu),但該類(lèi)方法主要是反映大面積、縱向上巖石宏觀孔隙性質(zhì)。③圖像分析法。該類(lèi)方法主要針對(duì)第1類(lèi)方法中的鑄體薄片法、掃描電鏡法、CT掃描法所生成的圖像,應(yīng)用數(shù)學(xué)形態(tài)學(xué)等學(xué)科分析孔隙的微觀結(jié)構(gòu)。

圖像分析法對(duì)算法的準(zhǔn)確性要求較高,特別是對(duì)識(shí)別喉道算法要求極高。因?yàn)楹淼雷R(shí)別的準(zhǔn)確性直接影響每個(gè)孔隙的屬性數(shù)據(jù)(包括面積、周長(zhǎng)、配位數(shù)等),以至影響巖石孔隙結(jié)構(gòu)參數(shù)。識(shí)別喉道的數(shù)學(xué)形態(tài)學(xué)的方法主要有流域分割方法[1-2]、中軸線算法[3]和采用數(shù)學(xué)形態(tài)學(xué)的灰度膨脹和腐蝕方法[4]等3種。這3種識(shí)別喉道算法的共同特點(diǎn)是先由極限腐蝕算法[1]確定孔隙位置,進(jìn)而尋找連接孔隙間喉道的位置。每種方法都有其優(yōu)勢(shì)和局限性。前2種方法對(duì)圖像的噪聲較為敏感、抗干擾性差,因而過(guò)分割嚴(yán)重,第3種方法的誤分割較多。

本文提出了一種改進(jìn)的識(shí)別喉道算法,并結(jié)合計(jì)算機(jī)圖形學(xué)算法分析計(jì)算孔隙的各項(xiàng)屬性數(shù)據(jù),提高識(shí)別喉道和分析孔隙各項(xiàng)屬性的準(zhǔn)確性和速度,并以巖石背散射掃描電鏡圖像為例,說(shuō)明其他圖像如鑄體薄片圖像等也可做類(lèi)似處理。

1 算法流程與數(shù)學(xué)描述

對(duì)巖石孔隙結(jié)構(gòu)分析的算法流程如圖1所示。

圖1 巖石孔隙結(jié)構(gòu)分析算法流程圖

1.1 去除巖石孔隙圖像的噪聲

巖石背散射掃描電鏡圖像受離散脈沖和椒鹽噪聲的影響較為嚴(yán)重。去除脈沖干擾及椒鹽噪聲最常用的算法是圖像處理學(xué)中的中值濾波[5]。中值濾波是一種去除噪聲的非線性處理方法,在某些條件下既可去除噪聲又可保護(hù)圖像細(xì)節(jié)和邊緣,能獲得較好的圖像效果。

1.2 提取巖石中的孔隙-圖像分割

圖像分割是指根據(jù)灰度、顏色、紋理、幾何形狀等信息把圖像分成各具特點(diǎn)、互不重疊的區(qū)域并提取出感興趣目標(biāo)的技術(shù)和過(guò)程。即在一幅畫(huà)圖像中把目標(biāo)從背景中分離出來(lái),以便于對(duì)感興趣目標(biāo)進(jìn)行更進(jìn)一步分析。圖2是所要處理的原始巖石孔隙結(jié)構(gòu)圖像,目標(biāo)是提取出圖像中的孔隙。

圖2 原始圖像(巖石背散射掃描電鏡圖像)

圖3 原始圖像的灰度直方圖

圖3是圖2的灰度直方圖,從直方圖中可以看到明顯的雙峰,但雙峰的灰度值相差極大,有寬且平的谷底。

根據(jù)直方圖的特征,采用基于迭代的圖像閾值分割算法[5]對(duì)原始孔隙圖像進(jìn)行分割,得到圖像中的孔隙。迭代的圖像閾值分割算法步驟如下。

第1步:求出圖像最小灰度值Rmin和最大灰度值Rmax,計(jì)算初始閾值

T0=(Rmin-Rmax)/2

(1)

第2步:根據(jù)閾值將圖像分割成目標(biāo)和背景2部分,示出2部分的平均灰度值。

(2)

(3)

式中,R(i,j)為圖像上(i,j)點(diǎn)的灰度值;N(i,j)為(i,j)點(diǎn)的權(quán)重系數(shù),一般N(i,j)為R(i,j)的個(gè)數(shù);TK為閾值。

第3步:重新選擇閾值TK+1。

TK+1=(R0-RG)/2

(4)

第4步:循環(huán)第2步到第4步,當(dāng)TK=TK+1則結(jié)束,即可獲得最佳閾值來(lái)對(duì)進(jìn)行分割。

圖4是圖2提取孔隙后的圖像,其中綠色部分表示孔隙,非綠色部分(灰色部分)表示巖石。

圖4 提取孔隙

1.3 去除孔隙二值圖像的噪聲

通常,有噪聲的圖像用閾值二值化所得到的邊界往往是很不平滑的,物體區(qū)域具有一些錯(cuò)判孔洞,背景區(qū)域上則散布著一些小的噪聲物體。這些噪聲被稱(chēng)為前景噪聲(如沙眼噪聲)和背景噪聲(如胡椒狀噪聲)。在識(shí)別喉道前,先要對(duì)孔隙的二值圖像去噪。采用數(shù)學(xué)形態(tài)學(xué)的開(kāi)運(yùn)算[1][先腐蝕,見(jiàn)式(5),后膨脹,見(jiàn)式(6)、式(7)]去除胡椒狀噪聲

A?B={x:B+x?A}

(5)

A⊕B=[AC?(-B)]C

(6)

A°B=[A?B]⊕B

(7)

采用閉運(yùn)算[1](先膨脹,后腐蝕)[見(jiàn)式(8)],去除沙眼噪聲

A·B=[A⊕(-B)]?(-B)

(8)

如果連續(xù)采用開(kāi)和閉運(yùn)算可以顯著改善二值圖像中的前、背景噪聲的情況。

1.4 識(shí)別喉道-分割孔隙

巖石圖像中,多個(gè)孔隙被喉道連接成連通區(qū)域,應(yīng)用改進(jìn)的識(shí)別喉道算法,把連通的孔隙分割開(kāi),以便能更準(zhǔn)確計(jì)算每個(gè)孔隙的各項(xiàng)屬性數(shù)據(jù)。改進(jìn)識(shí)別喉道中軸線算法流程圖見(jiàn)圖5。

圖5 改進(jìn)的識(shí)別喉道中軸線算法流程圖

A*B=(A?E)∩(AC?F),E∩F=?

(9)

B∶A=∪{(B+h)∩A∶h∈H}

(10)

采用上述算法識(shí)別圖4中孔隙的喉道,識(shí)別結(jié)果如圖6所示。

圖6 識(shí)別喉道后的圖像(紅色像素表示喉道)

1.5 計(jì)算孔隙各項(xiàng)屬性-圖像分析與理解

1.5.1 計(jì)算孔隙總數(shù)及孔隙面積

計(jì)算孔隙總數(shù)及孔隙面積的算法是用不同的顏色表示不同的孔隙,這為計(jì)算每個(gè)孔隙的周長(zhǎng)及喉道提供基礎(chǔ)。

(1) 查找未計(jì)算的孔隙像素。從上一個(gè)種子點(diǎn)(若剛開(kāi)始掃描圖像,則從圖像的(0,0)像素點(diǎn))開(kāi)始對(duì)圖像進(jìn)行由左向右、由上向下掃描,如果掃描到未計(jì)算的孔隙像素(綠色像素)時(shí),設(shè)置該像素為種子點(diǎn),且孔隙總數(shù)加1,并執(zhí)行(2),否則表示所有的孔隙已經(jīng)計(jì)算完成。

(2) 計(jì)算孔隙顏色。計(jì)算得到一個(gè)圖像中未使用過(guò)的顏色,且該顏色不能為綠色、紅色和灰度色。

(3) 計(jì)算孔隙面積。根據(jù)(1)中查找到的種子點(diǎn),使用(2)中計(jì)算出的孔隙顏色,采用計(jì)算機(jī)圖形學(xué)的種子填充算法,記錄下孔隙被填充的像素個(gè)數(shù)即為孔隙包含的像素個(gè)數(shù)。根據(jù)像素個(gè)數(shù)及比例值計(jì)算得到孔隙面積值。

1.5.2 計(jì)算每個(gè)孔隙的配位數(shù)及喉道長(zhǎng)度

一個(gè)孔隙的配位數(shù)就是與一個(gè)孔隙連接的喉道個(gè)數(shù)。

第1步:查找喉道。從上一個(gè)紅色像素(若剛開(kāi)始掃描圖像,則從圖像(0,0)像素)開(kāi)始對(duì)圖像進(jìn)行由左向右、由上向下掃描,如果掃描到喉道像素(紅色像素)時(shí),執(zhí)行第2步,否則繼續(xù)掃描圖像,直至掃描完整個(gè)圖像。

第2步:得到該喉道分割的孔隙。得到第1步找到的紅色像素的4個(gè)近鄰的顏色值,從這4個(gè)顏色值中找到該喉道分割開(kāi)的孔隙顏色值,對(duì)這2個(gè)孔隙的喉道各加1,并根據(jù)喉道的2個(gè)端點(diǎn)值計(jì)算出喉道的長(zhǎng)度。

第3步:刪除喉道。把圖像中該喉道所經(jīng)過(guò)的像素點(diǎn)的顏色設(shè)置為白色(或設(shè)置為原始巖石孔隙結(jié)構(gòu)圖像中的顏色)。返回執(zhí)行第1步。

1.5.3 計(jì)算孔隙周長(zhǎng)

運(yùn)用數(shù)學(xué)形態(tài)學(xué)的邊緣跟蹤(又稱(chēng)輪廓跟蹤)算法計(jì)算每個(gè)孔隙的周長(zhǎng)。邊緣跟蹤算法如下。

第1步:按從上到下、從左到右的順序掃描圖像,尋找沒(méi)有標(biāo)記跟蹤結(jié)束記號(hào)的第1個(gè)邊界起始點(diǎn)A0,A0是具有最小行和列值的邊界點(diǎn)。定義一個(gè)掃描方向變量dir,該變量用于記錄上一步中沿著前一個(gè)邊界點(diǎn)到當(dāng)前邊界點(diǎn)的移動(dòng)方向,其初始化取值:對(duì)4連通區(qū)域取dir=3;對(duì)8連通區(qū)域取dir=7。

第2步:按逆時(shí)針?lè)较蛩阉鳟?dāng)前象素的3×3鄰域,其起始搜索方向設(shè)定:對(duì)4連通區(qū)域取(dir+3) mod 4;對(duì)8連通區(qū)域,若dir為奇數(shù)取(dir+7) mod 8;若dir為偶數(shù)去(dir+6) mod 8;

在3×3鄰域中搜索到的第1個(gè)與當(dāng)前像素值相同的像素便為新的邊界點(diǎn)An,同時(shí)更新變量dir為新的方向值。

第3步:如果An等于第2個(gè)邊界點(diǎn)A1且前1個(gè)邊界點(diǎn)An-1等于第1個(gè)邊界點(diǎn)A0,則停止搜索,結(jié)束跟蹤,否則重復(fù)步驟2繼續(xù)搜索。

第4步:由邊界點(diǎn)A0、A1、A2、…、An-2構(gòu)成的邊界便為要跟蹤的邊界。

圖7是以圖6為基礎(chǔ)計(jì)算孔隙各項(xiàng)屬性后的圖像。

圖7 計(jì)算孔隙各項(xiàng)屬性后的圖像

2 巖石孔隙結(jié)構(gòu)分析結(jié)果

根據(jù)每個(gè)孔隙的各項(xiàng)屬性數(shù)據(jù),可計(jì)算出如表1所示的18項(xiàng)孔隙結(jié)構(gòu)參數(shù)表、表2所示的孔徑分布測(cè)定數(shù)據(jù)。

表1 孔隙結(jié)構(gòu)參數(shù)

表2 孔徑分布測(cè)定數(shù)據(jù)表

3 在巖石粒度分析中的應(yīng)用

該算法除了可以分析巖石孔隙圖像中的孔隙外,還可分析巖石顆粒。

圖8是圖2提取巖石顆粒的結(jié)果。圖9為分割巖石顆粒的結(jié)果。計(jì)算巖石顆粒的各項(xiàng)屬性包括顆粒面積、周長(zhǎng)、粒徑,圖10中一種顏色表示1個(gè)顆粒。

圖8 提取巖石顆粒(綠色表示顆粒)

根據(jù)每個(gè)巖石顆粒的屬性值可計(jì)算出如表3所示的9項(xiàng)粒度結(jié)構(gòu)參數(shù)以及圖11所示的粒度面積分布頻率圖。

表3 粒度結(jié)構(gòu)參數(shù)

圖9 巖石分割后的圖像(紅色像素為分割線)

圖10 計(jì)算巖石顆粒各項(xiàng)屬性后的圖像

圖11 粒度面積分布頻率直方圖

4 結(jié)束語(yǔ)

提出了一套完整的分析巖石孔隙結(jié)構(gòu)圖像的算法流程,特別是改進(jìn)了識(shí)別喉道的算法,大大加快了計(jì)算速度,如在聯(lián)想W500(處理器為Intel(R) Core(TM)2 Duo CPU P8600 @ 2.40 GHz,內(nèi)存2 GB)的筆記本上,對(duì)圖4(1 024像素×768像素)識(shí)別喉道,該算法用時(shí)0.853 s,流域分割方法用時(shí)2.156 s,中軸線算法用時(shí)2.785 s,灰度膨脹和腐蝕方法用時(shí)4.267 s。

應(yīng)用本文算法開(kāi)發(fā)的孔隙結(jié)構(gòu)圖像分析軟件已在中國(guó)石油勘探開(kāi)發(fā)研究院使用多年,分析各類(lèi)樣品數(shù)千個(gè)。實(shí)例應(yīng)用表明,這些算法抗干擾能力強(qiáng),可快速處理各種復(fù)雜類(lèi)型的巖石圖像,并準(zhǔn)確獲取孔隙結(jié)構(gòu)、粒度參數(shù)等信息,有效支撐了相關(guān)科學(xué)研究與生產(chǎn)任務(wù)。

符號(hào)說(shuō)明。T0:初始閾值;TK:迭代第K步閾值;Rmin:圖像最小灰度值;Rmax:圖像最大灰度值;R0:第K步圖像目標(biāo)的平均灰度值;RG:第K步圖像背景的平均灰度值;R(i,j):圖像上(i,j)點(diǎn)的灰度值;N(i,j):圖像上(i,j)點(diǎn)的權(quán)重系數(shù);A:輸入圖像;AC:為A的補(bǔ)集;B:結(jié)構(gòu)元素或結(jié)構(gòu)元素對(duì)(E,F);-B:將B相對(duì)原點(diǎn)旋轉(zhuǎn)180°;E:探測(cè)圖像內(nèi)部的結(jié)構(gòu)元素;F:探測(cè)圖像外部的結(jié)構(gòu)元素;H:為A的子集,第1次膨脹的種子點(diǎn)。

參考文獻(xiàn):

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

[2] 印勇, 李阿瓊. 一種粘連血細(xì)胞圖像分割新方法 [J]. 計(jì)算機(jī)工程與應(yīng)用, 2009, 45(35): 173-175.

[3] HU Dong, Mustafa Touati, Martin J Blunt. Pore Network Modeling: Analysis of Pore Size Distribution of Arabian Core Samples [J]. SPE 105156.

[4] 劉莉莉, 王錚. 一種適合血細(xì)胞圖像分割的改進(jìn)流域分割算法 [J]. 微電子學(xué)與計(jì)算機(jī), 2010, 27(11).

[5] 王志明. 數(shù)字圖像處理與分析 [M]. 北京: 清華大學(xué)出版社, 2012.

[6] 王雪莉, 張吉群. 高壓地層微觀模型系統(tǒng)中的背景校正 [J]. 科學(xué)技術(shù)與工程, 2012, 12(25): 6498-6502.

[7] Dmitry B Silin, Guodong Jin, Tad W Patzek. Robust Determination of the Pore Space Morphology in Sedimentary Rocks [J]. SPE 84296.

[8] AI Ibrahim M A, Hurley N F, Zhao W, et al. An Automated Petrographic Image Analysis System: Capillary Pressure Curves Using Confocal Microscopy [J]. SPE 159180.

[9] 桑卡, 赫拉瓦卡, 博伊爾. 圖像處理、分析與機(jī)器視覺(jué) [M]. 艾海舟, 蘇延超, 譯. 3版. 北京: 清華大學(xué)出版社, 2011.

[10] Suicmez V S, Touati M. Pore Network Modeling: A New Technology for SCAL Predictions and Interpretations [J]. SPE 110961.

[11] Liu Jianjun, Lin Lijun, Ji Youjun. Using Rock SEM Image to Create Pore-scale Finite Element Calculation Mesh [C]∥2011 International Conference on Physics Science and Technology (ICPST 2011).

[12] 熊承仁, 唐輝明, 劉寶琛, 等. 利用SEM照片獲取土的孔隙結(jié)構(gòu)參數(shù) [J]. 中國(guó)地質(zhì)大學(xué)學(xué)報(bào): 地球科學(xué), 2007, 32(3): 415-419.

[13] 尼克松, 阿瓜多. 特征提取與圖像處理 [M]. 李實(shí)英, 楊高波, 譯. 2版. 北京: 電子工業(yè)出版社, 2010.

[14] 帕科爾. 景麗譯. 圖像處理與計(jì)算機(jī)視覺(jué)算法及應(yīng)用 [M]. 2版. 北京: 清華大學(xué)出版社, 2011.

[15] 赫恩, 巴克. 計(jì)算機(jī)圖形學(xué) [M]. 宋繼強(qiáng), 蔡敏, 譯. 3版. 北京: 電子工業(yè)出版社, 2011.

[16] 汪燦, 劉艷敏, 祝艷波. SEM照片孔隙參數(shù)提取技術(shù)研究 [J]. 安全與環(huán)境工程, 2011, 18(3): 117-120.