基于二值图像处理技术的砂土颗粒统计方法

邹 峰，王炳辉，姜朋明

(江苏科技大学土木工程与建筑学院，江苏镇江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.