基于二值圖像處理技術的砂土顆粒統計方法

鄒 峰，王炳輝，姜朋明

(江蘇科技大學土木工程與建筑學院，江蘇鎮江212003)

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

1 數字圖像顆粒統計方法

1.1 圖像的數字表示

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

式中:n，m分別為圖像中所包含的像素點陣的行數和列數.

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

1.2 圖像的類型轉化與分割

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

圖像的分割主要運用到了閾值分割法，閾值分割法是傳統的圖像分割方法，其優點是實現簡單、計算量小、性能比較穩定，在圖像分割中得到廣泛應用.閾值分割時，圖像灰度被分成不同等級，通過設定不同的閾值把像素點分為若干類.當小于閾值的灰度像素構成目標時，則大于閾值的灰度像素就構成背景.于是，邊界上的點都至少有一個相鄰點不屬于目標對象.所選閾值要使目標和背景的錯誤分割像素個數最少.

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

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

閾值T確定方法是得到二值圖像的關鍵，常用方法有目視檢查法、反復試驗法、自動迭代法、Otus自適應法等4種.文中使用了目視檢查法，該方法適用于背景、目標灰度級差別較為明顯的情況.

1.3 Matlab軟件處理圖像

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

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

圖2 經過染色的砂土Fig.2 Dyed sand

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

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

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

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

1.4 連通域搜索

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

文中程序結合八鄰域搜索規則［9］按照順時針的方向，在全局范圍內來尋找0像素區域.各個方向的坐標表示和方向碼如圖5，6所示.

圖5 八鄰域的坐標表示Fig.5 Eight neighborhood coordinates expressions

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

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

運行程序對圖1處理得到圖中的黑色顆粒有78個，與原圖像中的黑色顆粒個數相符.

圖7 調用DFSea算法流程Fig.7 Call DFSea algorithm flow

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

1.5 基于聚類算法的顆粒統計

運用聚類算法［10］來還原出原圖像并進行對比，從而可以實現驗證統計的顆粒數量是否準確以及判斷程序的有效性.砂土雖然形狀各異，但是仍然可以將它們看成是一粒粒圓粒狀，并且在染色前通過篩分器已經將大小相近、形狀相似的同級配顆粒篩分出來，更加適合聚類算法.K-means算法是基于距離的聚類算法，采用距離作為相似性的評價指標，即認為兩個對象的距離越近，其相似度就越大.該算法認為簇是由距離靠近的對象組成的，因此把得到緊湊且獨立的簇作為最終目標.k個初始類聚類中心點的選取對聚類結果具有較大的公式影響，因為在該算法第一步中是隨機的選取任意k個對象作為初始聚類的中心，初始地代表一個簇.該算法在每次迭代中對數據集中剩余的每個對象，根據其與各個簇中心的距離將每個對象重新賦給最近的簇.當考察完所有數據對象后，一次迭代運算完成，新的聚類中心被計算出來.如果在一次迭代前后，原聚類中心的值沒有發生變化，說明算法已經收斂.算法過程如下:

1)從搜索到的連通域中隨機選取K個點作為質心，K值為連通域內搜索顆粒時統計出的顆粒個數.

2)對剩余的每個像素測量其到每個質心的距離，并把它歸到最近的質心的類.

3)用坐標均值的算法重新計算已經得到的各個類的質心.

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

求得質心后，利用每個顆粒包括的像素數除以圓周率求的顆粒半徑，在全圖范圍內通過質心與半徑還原顆粒，呈現染色顆粒的分布情況.圖像準確的還原了黑色顆粒的分布，與原圖像對比如圖9，10.

圖9 原圖像Fig.9 Original image

圖10 還原圖像Fig.10 Restore image

上述利用程序對比較規整的圓顆粒進行了搜索還原，得出的黑色顆粒個數與已知的黑色顆粒個數相同，都為78個，并且還原的黑色顆粒分布的圖像與原圖一致，驗證了程序的準確性與有效性，為程序處理更復雜的砂土顆粒奠定了基礎.

2 應用實例

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

圖11 原圖像Fig.11 Original image

圖12 灰度圖像Fig.12 Grayscale image

由圖像可以知道，黑色砂土顆粒與圓粒圖像中黑色顆粒有顯著不同，圓粒圖像中顆粒形狀基本都是圓形，形狀都比較相似，而染色砂土顆粒的形狀各異，大小不一.另外，染色砂土顆粒之間還會相互疊加，導致一塊區域中不只有一個顆粒，與之前處理圓粒圖像情況不同.因此針對砂土顆粒的形狀不一，相互疊加的復雜情況，應對程序進行優化.

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

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

圖13 二值圖像Fig.13 Binary image

圖14 模擬圖像Fig.14 Simulation image

3 結論

1)結合實驗方法和數字圖像技術，獲得染色顆粒在圖像中的分布，基本還原了真實砂土顆粒分布情況，通過程序得到了具體顆粒個數，為進一步研究砂土內部結構性質變化與顆粒變化的內在聯系作鋪墊.

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

3)根據砂土特性對編制的Matlab軟件程序進一步的優化可以在處理圖像時更加有效與準確.

References)

［1］ 徐文杰，岳中琦，胡瑞林.基于數字圖像的土、巖和混凝土內部結構定量分析和力學數值計算的研究進展［J］.工程地質學報，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］ 李元海，朱合華，上野勝，等.基于圖像相關分析的砂土模型試驗變形場量測［J］.巖土工程學報，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］ 李元海，靖洪文.基于數字散斑相關法的變形量測軟件研制及應用［J］.中國礦業大學學報，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］ 張嘎，張建民，梁東方.土與結構接觸面試驗中的土顆粒細觀運動測量［J］.巖土工程學報，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］ 周健，余榮傳，賈敏才.基于數字圖像技術的砂土模型試驗細觀結構參數測量［J］.巖土工程學報，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］ 周健，陳小亮，賈敏才，等.有地下結構的飽和砂土液化宏細觀離心機試驗［J］.巖土工程學報，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］ 邱琳，吳華山，陳效民，等.利用染色示蹤和圖像處理技術對土壤大孔隙進行定量研究［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］ 劉永祿，邵龍潭，潘石.三軸實驗土樣變形數字圖像測量(DDIM)的改進算法［C］∥第十一屆全國實驗力學學術會議，大連，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.