基于二值圖像處理技術(shù)的砂土顆粒統(tǒng)計(jì)方法

鄒 峰，王炳輝，姜朋明

(江蘇科技大學(xué)土木工程與建筑學(xué)院，江蘇鎮(zhèn)江212003)

土體細(xì)觀組成及結(jié)構(gòu)很大程度上控制了其宏觀的力學(xué)響應(yīng)和破壞機(jī)理及過程［1］，砂土顆粒之間的緊密程度在宏觀上表現(xiàn)出不同的力學(xué)性質(zhì)，顆粒間緊密的砂土(即密砂)在宏觀上表現(xiàn)出剪脹性，相反(如松砂)則表現(xiàn)出壓硬性.砂土的密實(shí)度雖可通過土工常規(guī)試驗(yàn)獲得，但在砂土密實(shí)度的變化過程中，如砂土在剪切過程中剪切帶處砂土的密實(shí)度，則無法采用該方法進(jìn)行測(cè)試.隨著測(cè)試技術(shù)的發(fā)展，通過土體的微細(xì)觀顆粒形狀、分布、運(yùn)移等特征分析獲得土體力學(xué)性質(zhì)的數(shù)字圖像識(shí)別技術(shù)逐漸被學(xué)者所采用.國(guó)內(nèi)李元海等［2－3］較早利用數(shù)字圖像相關(guān)分析開展砂土模型試驗(yàn)的變形場(chǎng)研究，并自行研制了數(shù)字散斑相關(guān)變形量測(cè)實(shí)用軟件系統(tǒng)Photoinfor，為砂土模型局部化變形、剪切帶和漸進(jìn)破壞過程的量化分析提供了一條新途徑.文獻(xiàn)［4］中提出一種補(bǔ)零相關(guān)分析算法，并在此基礎(chǔ)上實(shí)現(xiàn)了土與結(jié)構(gòu)接觸面試驗(yàn)中的土顆粒運(yùn)動(dòng)自動(dòng)化測(cè)量.文獻(xiàn)［5－6］中則在建立砂顆粒與孔隙判別標(biāo)準(zhǔn)的基礎(chǔ)上，采用形態(tài)學(xué)分水嶺分割方法等數(shù)字圖像處理技術(shù)進(jìn)行砂土顆粒及其孔隙識(shí)別的研究，并與離心振動(dòng)臺(tái)試驗(yàn)結(jié)合，從顆粒細(xì)觀組構(gòu)參量的動(dòng)態(tài)變化研究了飽和砂土地基液化的細(xì)觀機(jī)制.土體的內(nèi)部微結(jié)構(gòu)的變化必然引起外部表觀的變化，如砂土顆粒數(shù)量的變化.因此，為了研究砂土內(nèi)部微結(jié)構(gòu)發(fā)生的變化，可以從其表觀的顆粒數(shù)量入手，進(jìn)行顆粒統(tǒng)計(jì)，以顆粒數(shù)量變化來間接探究砂土結(jié)構(gòu)性質(zhì)的演化過程.

1 數(shù)字圖像顆粒統(tǒng)計(jì)方法

1.1 圖像的數(shù)字表示

圖像的數(shù)字表示是數(shù)字圖像在計(jì)算機(jī)中由一系列矩形排列的像素點(diǎn)構(gòu)成，在灰度圖像中每個(gè)像素點(diǎn)對(duì)應(yīng)一個(gè)整數(shù)值，用以表示點(diǎn)的亮度，即灰度.常見的256色，其灰度值為0～255.整個(gè)圖像便由具有不同灰度值的像素點(diǎn)陣構(gòu)成，該像素點(diǎn)陣的灰度值構(gòu)成了一個(gè)離散函數(shù)f(i，j)(i，j分別代表像素點(diǎn)整個(gè)圖像對(duì)應(yīng)像素點(diǎn)陣中的行、列號(hào)).

式中:n，m分別為圖像中所包含的像素點(diǎn)陣的行數(shù)和列數(shù).

對(duì)于一個(gè)二值化圖像(只有黑白兩色)其灰度值僅有0(代表黑色)和255(代表白色)兩種.因此，經(jīng)過一定的圖像處理生成的數(shù)字圖像的各個(gè)灰度值代表了圖像所包含的不同信息，也即反映了材料的細(xì)觀內(nèi)部結(jié)構(gòu)特征.

1.2 圖像的類型轉(zhuǎn)化與分割

圖像類型轉(zhuǎn)換最常用的是圖像灰度化和圖像二值化.圖像二值化在圖像分割中說明，這里主要說明圖像灰度化的理論基礎(chǔ)及其Matlab實(shí)現(xiàn)方法.將彩色RGB圖像轉(zhuǎn)換成灰度圖像的過程叫圖像的灰度化.把灰度圖像中不同灰度值像素個(gè)數(shù)所占比例的柱狀圖，稱為灰度直方圖.由于R，G，B的取值范圍是0～255，所以灰度級(jí)別只有256級(jí)、灰度圖像僅能表現(xiàn)256種灰度.常用灰度化處理方法有如下4種:分量法、最大值法、平均值法、加權(quán)平均法.文中采用平均值法，將彩色圖像三分量亮度的平均值作為灰度值，即:

圖像的分割主要運(yùn)用到了閾值分割法，閾值分割法是傳統(tǒng)的圖像分割方法，其優(yōu)點(diǎn)是實(shí)現(xiàn)簡(jiǎn)單、計(jì)算量小、性能比較穩(wěn)定，在圖像分割中得到廣泛應(yīng)用.閾值分割時(shí)，圖像灰度被分成不同等級(jí)，通過設(shè)定不同的閾值把像素點(diǎn)分為若干類.當(dāng)小于閾值的灰度像素構(gòu)成目標(biāo)時(shí)，則大于閾值的灰度像素就構(gòu)成背景.于是，邊界上的點(diǎn)都至少有一個(gè)相鄰點(diǎn)不屬于目標(biāo)對(duì)象.所選閾值要使目標(biāo)和背景的錯(cuò)誤分割像素個(gè)數(shù)最少.

閾值分割法中，二值化方法最為常用.所謂圖像二值化，就是將圖像中任一位置的灰度值設(shè)置為0或1，整個(gè)圖像呈現(xiàn)出明顯的黑白效果.假設(shè)灰度直方圖對(duì)應(yīng)圖像f(x，y)，由較暗的背景和較亮的目標(biāo)所組成.f(x，y)≥T的點(diǎn)(x，y)稱為對(duì)象點(diǎn)，其他點(diǎn)則稱為背景點(diǎn).經(jīng)閾值處理后的圖像g(x，y)可定義為:

式中:f(x，y)為二維灰度函數(shù)，T為閾值.

閾值T確定方法是得到二值圖像的關(guān)鍵，常用方法有目視檢查法、反復(fù)試驗(yàn)法、自動(dòng)迭代法、Otus自適應(yīng)法等4種.文中使用了目視檢查法，該方法適用于背景、目標(biāo)灰度級(jí)差別較為明顯的情況.

1.3 Matlab軟件處理圖像

運(yùn)用編制的Matlab圖像處理程序?qū)D像進(jìn)行數(shù)字化分析.整個(gè)程序流程分3個(gè)步驟:圖像分割、顆粒搜索、圖像還原.首先試驗(yàn)利用編制的程序來處理圓粒圖像，統(tǒng)計(jì)該圖像中黑色的圓顆粒個(gè)數(shù).如圖1，此圖像中顆粒形狀都接近圓形，且可以直觀的得到黑色顆粒個(gè)數(shù)，便于驗(yàn)證程序統(tǒng)計(jì)的準(zhǔn)確性.接著應(yīng)用程序于實(shí)際砂土圖像中.圖像是拍攝混合有相同粒徑的染色顆粒的砂土而得［7－8］.試驗(yàn)將標(biāo)準(zhǔn)砂進(jìn)行篩分，把篩分好的粒徑為0.5～1 mm的顆粒進(jìn)行染色，染色完全并烘干后把染色顆粒倒入其他未染色的砂土顆粒中混合均勻，如圖2，拍攝砂土圖像，最后把圖像導(dǎo)入計(jì)算機(jī)進(jìn)行程序處理.

圖1 圓粒圖像Fig.1 Rounded grain image

圖2 經(jīng)過染色的砂土Fig.2 Dyed sand

將圖像讀入程序，把圖像轉(zhuǎn)化為灰度圖像.顏色較深顆粒所包括的像素的灰度值都比較小，在0到某一閾值范圍內(nèi)，選取閾值把灰度圖像轉(zhuǎn)化為二值圖像，即小于該閾值的顯示為黑色，灰度值為0，大于該閾值的顯示為白色，灰度值為1，從而將要找的顆粒劃分出來.

將圖1轉(zhuǎn)化為灰度圖像，如圖3.由于黑色顆粒所包括的像素灰度值都較小在0～45范圍內(nèi)，與其他顆粒像素相差大，所以根據(jù)上述閾值分割法中的目視檢查法來選取閾值，輸入閾值45，獲得二值圖像，如圖4.

圖3 灰度處理過的圖像Fig.3 Gradation processed image

圖4 閾值分割后的圖像Fig.4 Threshold split-image

1.4 連通域搜索

將二值圖像中黑色顆粒找出，并統(tǒng)計(jì).由于圖像中的顆粒形狀比較相近，搜索到一塊0像素組成的獨(dú)立區(qū)域即為一個(gè)顆粒的像素集合，因此只要統(tǒng)計(jì)出0像素區(qū)域個(gè)數(shù)就為黑色顆粒個(gè)數(shù).搜索區(qū)域采用的是八鄰域深度搜索的方法即DFSea深度搜索算法，該算法圍繞0像素組成的區(qū)域邊界進(jìn)行搜索并可以統(tǒng)計(jì)出區(qū)域內(nèi)所有0像素個(gè)數(shù).八鄰域算法的運(yùn)算規(guī)則為將某一個(gè)方位的鄰域作為起始搜索點(diǎn)，環(huán)繞該點(diǎn)八個(gè)方向，按自行編制的順時(shí)針或是逆時(shí)針方向搜索.如點(diǎn)P(x，y)為0像素區(qū)域的一個(gè)邊界點(diǎn)，則P(x，y)的下一邊界點(diǎn)必在其八鄰域內(nèi)，因此可以根據(jù)程序首先找到該點(diǎn)左下角的邊界點(diǎn)作為搜索起點(diǎn)，按逆時(shí)針方向，自下而上、從左至右，搜索其八鄰域，找到下一邊界點(diǎn)，然后以此邊界點(diǎn)為當(dāng)前點(diǎn)繼續(xù)搜索，重復(fù)此過程.搜索時(shí)遇到一個(gè)0先統(tǒng)計(jì)值加1，然后把像素置為1，這樣就能將區(qū)域內(nèi)所有0像素統(tǒng)計(jì)出來且不會(huì)造成死循環(huán)一直重復(fù)統(tǒng)計(jì).

文中程序結(jié)合八鄰域搜索規(guī)則［9］按照順時(shí)針的方向，在全局范圍內(nèi)來尋找0像素區(qū)域.各個(gè)方向的坐標(biāo)表示和方向碼如圖5，6所示.

圖5 八鄰域的坐標(biāo)表示Fig.5 Eight neighborhood coordinates expressions

圖6 八鄰域方向碼Fig.6 Eight neighborhood direction code

上圖中方向碼中心即為起始點(diǎn)，0，1，2，3，4，5，6，7為起始點(diǎn)的八個(gè)鄰域的方向碼，分別代表了起始點(diǎn)的上、下、左、右、左上、右上、左下、右下8個(gè)方向.圖5運(yùn)用數(shù)學(xué)公式表示各個(gè)方向，坐標(biāo)起始點(diǎn)為(x，y)，i為定義的方向碼.編制的程序是從左下開始，逆時(shí)針運(yùn)行搜索找到區(qū)域邊界并由外向里，直至全部覆蓋.當(dāng)一塊區(qū)域內(nèi)所有0像素都統(tǒng)計(jì)出來并依次置為1后，程序會(huì)繼續(xù)全局搜索，直至找到下個(gè)區(qū)域.調(diào)用深度搜索的算法流程及具體DFSea算法流程如圖7，8.

運(yùn)行程序?qū)D1處理得到圖中的黑色顆粒有78個(gè)，與原圖像中的黑色顆粒個(gè)數(shù)相符.

圖7 調(diào)用DFSea算法流程Fig.7 Call DFSea algorithm flow

圖8 DFSea算法流程Fig.8 DFSea algorithm flow

1.5 基于聚類算法的顆粒統(tǒng)計(jì)

運(yùn)用聚類算法［10］來還原出原圖像并進(jìn)行對(duì)比，從而可以實(shí)現(xiàn)驗(yàn)證統(tǒng)計(jì)的顆粒數(shù)量是否準(zhǔn)確以及判斷程序的有效性.砂土雖然形狀各異，但是仍然可以將它們看成是一粒粒圓粒狀，并且在染色前通過篩分器已經(jīng)將大小相近、形狀相似的同級(jí)配顆粒篩分出來，更加適合聚類算法.K-means算法是基于距離的聚類算法，采用距離作為相似性的評(píng)價(jià)指標(biāo)，即認(rèn)為兩個(gè)對(duì)象的距離越近，其相似度就越大.該算法認(rèn)為簇是由距離靠近的對(duì)象組成的，因此把得到緊湊且獨(dú)立的簇作為最終目標(biāo).k個(gè)初始類聚類中心點(diǎn)的選取對(duì)聚類結(jié)果具有較大的公式影響，因?yàn)樵谠撍惴ǖ谝徊街惺请S機(jī)的選取任意k個(gè)對(duì)象作為初始聚類的中心，初始地代表一個(gè)簇.該算法在每次迭代中對(duì)數(shù)據(jù)集中剩余的每個(gè)對(duì)象，根據(jù)其與各個(gè)簇中心的距離將每個(gè)對(duì)象重新賦給最近的簇.當(dāng)考察完所有數(shù)據(jù)對(duì)象后，一次迭代運(yùn)算完成，新的聚類中心被計(jì)算出來.如果在一次迭代前后，原聚類中心的值沒有發(fā)生變化，說明算法已經(jīng)收斂.算法過程如下:

1)從搜索到的連通域中隨機(jī)選取K個(gè)點(diǎn)作為質(zhì)心，K值為連通域內(nèi)搜索顆粒時(shí)統(tǒng)計(jì)出的顆粒個(gè)數(shù).

2)對(duì)剩余的每個(gè)像素測(cè)量其到每個(gè)質(zhì)心的距離，并把它歸到最近的質(zhì)心的類.

3)用坐標(biāo)均值的算法重新計(jì)算已經(jīng)得到的各個(gè)類的質(zhì)心.

4)迭代2～3步直至新的質(zhì)心與原質(zhì)心相等或小于指定閾值，算法結(jié)束.

求得質(zhì)心后，利用每個(gè)顆粒包括的像素?cái)?shù)除以圓周率求的顆粒半徑，在全圖范圍內(nèi)通過質(zhì)心與半徑還原顆粒，呈現(xiàn)染色顆粒的分布情況.圖像準(zhǔn)確的還原了黑色顆粒的分布，與原圖像對(duì)比如圖9，10.

圖9 原圖像Fig.9 Original image

圖10 還原圖像Fig.10 Restore image

上述利用程序?qū)Ρ容^規(guī)整的圓顆粒進(jìn)行了搜索還原，得出的黑色顆粒個(gè)數(shù)與已知的黑色顆粒個(gè)數(shù)相同，都為78個(gè)，并且還原的黑色顆粒分布的圖像與原圖一致，驗(yàn)證了程序的準(zhǔn)確性與有效性，為程序處理更復(fù)雜的砂土顆粒奠定了基礎(chǔ).

2 應(yīng)用實(shí)例

通過實(shí)驗(yàn)獲取混有染色顆粒的砂土圖像，如圖11;將圖像灰度處理，如圖12.

圖11 原圖像Fig.11 Original image

圖12 灰度圖像Fig.12 Grayscale image

由圖像可以知道，黑色砂土顆粒與圓粒圖像中黑色顆粒有顯著不同，圓粒圖像中顆粒形狀基本都是圓形，形狀都比較相似，而染色砂土顆粒的形狀各異，大小不一.另外，染色砂土顆粒之間還會(huì)相互疊加，導(dǎo)致一塊區(qū)域中不只有一個(gè)顆粒，與之前處理圓粒圖像情況不同.因此針對(duì)砂土顆粒的形狀不一，相互疊加的復(fù)雜情況，應(yīng)對(duì)程序進(jìn)行優(yōu)化.

根據(jù)砂土形狀的無規(guī)則性，另外還有相互疊加的顆粒，統(tǒng)計(jì)顆粒時(shí)就不能僅僅將像素相加，故引入連通域的概念.考慮到顆粒既有一個(gè)顆粒，又有幾個(gè)顆粒疊在一起，故只需要關(guān)注找到由顆粒組成的連通域，再由連通域來確定顆粒個(gè)數(shù).運(yùn)行程序?qū)⒒叶葓D像(圖12)閾值分割，根據(jù)目視檢查法選取閾值為30，結(jié)合原圖像，設(shè)定每個(gè)染色顆粒包含有255個(gè)像素.得到二值圖像，如圖13.最后通過聚類算法得到模擬圖像，如圖14.

圖14基本還原了原圖像中染色顆粒分布.經(jīng)過程序處理可以清楚地觀察到該圖像的染色顆粒分布情況，并且得出在圖像中染色顆粒個(gè)數(shù)為857個(gè).

圖13 二值圖像Fig.13 Binary image

圖14 模擬圖像Fig.14 Simulation image

3 結(jié)論

1)結(jié)合實(shí)驗(yàn)方法和數(shù)字圖像技術(shù)，獲得染色顆粒在圖像中的分布，基本還原了真實(shí)砂土顆粒分布情況，通過程序得到了具體顆粒個(gè)數(shù)，為進(jìn)一步研究砂土內(nèi)部結(jié)構(gòu)性質(zhì)變化與顆粒變化的內(nèi)在聯(lián)系作鋪墊.

2)結(jié)合深度搜索算法和聚類算法而編制的Matlab程序能夠有效處理砂土顆粒圖像，搜索并還原出所要得到的對(duì)象，更加支持了利用數(shù)字圖像技術(shù)處理研究砂土內(nèi)部結(jié)構(gòu)性質(zhì)的研究思想.

3)根據(jù)砂土特性對(duì)編制的Matlab軟件程序進(jìn)一步的優(yōu)化可以在處理圖像時(shí)更加有效與準(zhǔn)確.

References)

［1］ 徐文杰，岳中琦，胡瑞林.基于數(shù)字圖像的土、巖和混凝土內(nèi)部結(jié)構(gòu)定量分析和力學(xué)數(shù)值計(jì)算的研究進(jìn)展［J］.工程地質(zhì)學(xué)報(bào)，2007，15(3):289－313.Xu Wenjie，Yue Zhongqi，Hu Ruilin.Research progress in the numerical calculation of digital images of the internal structure of rock and soil，the quantitative analysis of concrete and based on Mechanics［J］.Journal of Engineering Geology，2007，15(3):289 －313.(in Chinese)

［2］ 李元海，朱合華，上野勝，等.基于圖像相關(guān)分析的砂土模型試驗(yàn)變形場(chǎng)量測(cè)［J］.巖土工程學(xué)報(bào)，2004，26(1):36－41.Li Yuanhai，Zhu Hehua，Shang Yesheng，et al.Sand model test based on image correlation analysis of deformation field measurement［J］.Chinese Journal of Geotechnical Engineering，2004，26(1):36 －41.(in Chinese)

［3］ 李元海，靖洪文.基于數(shù)字散斑相關(guān)法的變形量測(cè)軟件研制及應(yīng)用［J］.中國(guó)礦業(yè)大學(xué)學(xué)報(bào)，2008，37(5):635－640.Li Yuanhai，Jing Hongwen.Based on digital speckle correlation method distortion gauging software development and application［J］.Journal of China University of Mining and Technology，2008，37(5):635 －640.(in Chinese)

［4］ 張嘎，張建民，梁東方.土與結(jié)構(gòu)接觸面試驗(yàn)中的土顆粒細(xì)觀運(yùn)動(dòng)測(cè)量［J］.巖土工程學(xué)報(bào)，2005，27(8):903－907.Zhang Ga，Zhang Jianming，Liang Dongfang.Soil particle surface in the experiment of meso movement measuring the interface between soil and structure［J］.Chinese Journal of Geotechnical Engineering，2005，27(8):903 －907.(in Chinese)

［5］ 周健，余榮傳，賈敏才.基于數(shù)字圖像技術(shù)的砂土模型試驗(yàn)細(xì)觀結(jié)構(gòu)參數(shù)測(cè)量［J］.巖土工程學(xué)報(bào)，2006，28(12):2047－2052.Zhou Jian，Yu Rongquan，Jia Mincai.Sand model test based on digital image measurement of microstructure parameters［J］.Chinese Journal of Geotechnical Engineering，2006，28(12):2047 －2052.(in Chinese)

［6］ 周健，陳小亮，賈敏才，等.有地下結(jié)構(gòu)的飽和砂土液化宏細(xì)觀離心機(jī)試驗(yàn)［J］.巖土工程學(xué)報(bào)，2012，34(3):392－399.Zhou Jian，Chen Xiaoliang，Jia Mincai，et al.Underground structures in saturated sand liquefaction macro and micro centrifuge［J］.Chinese Journal of Geotechnical Engineering，2012，34(3):392 － 399.(in Chinese)

［7］ 邱琳，吳華山，陳效民，等.利用染色示蹤和圖像處理技術(shù)對(duì)土壤大孔隙進(jìn)行定量研究［J］.土壤，2007(4):621－626.Qiu Lin，Wu Huashan，Chen Xiaoming，et al.Quantitative study on soil macropores by dye tracing and image processing technology［J］.Soil，2007(4):621 － 626.(in Chinese)

［8］ 劉永祿，邵龍?zhí)叮耸?三軸實(shí)驗(yàn)土樣變形數(shù)字圖像測(cè)量(DDIM)的改進(jìn)算法［C］∥第十一屆全國(guó)實(shí)驗(yàn)力學(xué)學(xué)術(shù)會(huì)議，大連，2005:1107－1112.

［9］ Rahhn Xavier J，Camanho P P.High strain rate characterization of unidirectional carbon-epoxy IM7-8552 in transverse compression and in-plane shear using digital image correlation ［J］.Mechanics of Materials，2010，42:1004－1019.

［10］ Hosseini A，Mostofinejad M.Effect of groove characteristics on CFRP-to-concrete bond behavior of EBROG joints:experimental study using particle image velocimetry(PIV)［J］.Construction and Building Materials，2013，49:364 －373.