基于数字图像的粒状材料细观组构特征分析技术

第30卷 第10期 岩 土 工 程 学 报 Vol.30 No.10 2008年 10月 Chinese Journal of Geotechnical Engineering Oct., 2008

基于数字图像的粒状材料细观组构特征分析技术

张连卫

1，3

，张建民

1，2

，张 嘎

1，2

（1. 清华大学岩土工程研究所，北京 100084；2. 清华大学水沙科学与水利水电工程国家重点试验室，北京 100084；

3. 北京科技大学土木与环境工程学院，北京 100083）

摘 要：针对由特殊截面形状的金属棒所组成的理想二维粒状材料，基于数字图像技术，提出了一种细观组构特征分析方法，并编制了程序IPFA，实现了颗粒识别、接触搜索与组构分析等功能。该方法首先对粒状材料试样的原始数字图像进行增强，以消除噪声干扰。在此基础上采用模板匹配技术进行颗粒识别，并将识别结果用于颗粒间的接触搜索，最后给出试样内颗粒长轴方向与接触法线方向等组构特征的统计分析结果。该方法特别适用于由规则颗粒组成的二维粒状材料，识别精度与效率均较高，可作为粒状材料的细观组构特征及其演化规律分析的有效工具。 关键词：粒状材料；细观组构；数字图像；颗粒识别

中图分类号：TU43 文献标识码：A 文章编号：1000–4548(2008)10–1555–05

作者简介：张连卫(1978– )，男，山东无棣人，讲师，主要从事粒状材料各向异性方面的研究。E-mail: zhanglianwei01@。

Microfabric analysis technique for granular materials based on digital images

ZHANG Lian-wei13，ZHANG Jian-min12，ZHANG Ga1

，

，

，2

(1. Institute of Geotechnical Engineering, Tsinghua University, Beijing 100084, China; 2. State Key Lab of Hydroscience and Engineering, Tsinghua University, Beijing 100084, China; 3. Civil and Environmental Engineering School, University of Science and Technology

Beijing, Beijing 100083, China)

Abstract: Based on the digital image processing technique, a method was presented for microfabric analysis for ideal granular materials composed of metal bars with specified sections. A code named IPFA was then developed to implement the functions such as particle detection, contact searching and microfabric analysis. By use of this method, the digital images of dense samples were firstly enhanced to clear the noises. The particles were then detected from the enhanced images using template matching algorithm presented here, and contact normals between particles were searched based on the particle detection results. Finally, the statistical analysis could be processed on two important elements of microscopic fabric, including directions of the long axis of particles and the contact normals between particles. The proposed method was especially suitable for ideal two-dimensional granular materials, and could be used as a powerful tool for the analysis of microfabric of granular materials and its evaluation because of its high precision and efficiency of particle detection. Key words: granular material; microfabric; digital image; particle detection

0 引 言

粒状材料的细观组构主要指颗粒与相应孔隙在空

其对粒状材料的宏观力学行间中的几何排列方式[1-3]，

为具有非常重要的影响。目前，在土力学与粒状介质力学研究领域，不少学者都试图揭示细观组构对宏观变形的影响规律。这有赖于对细观组构特征及其演化规律的定量分析。

组构通常不涉及颗粒之间的胶结等相互作用，因此可通过分析粒状材料试样的图像得到。随着测量技术的迅速发展，这种细部图像可通过多种方法获得[4-14]

。早期多采用光学照相放大技术[4]，针对黏土等颗

粒较细的材料则主要采用光学显微与电子显微技术。

近几十年数字照相技术发展迅速，分辨率提高，成本降低，逐渐被广泛用于获得粒状材料试样的细部图像。

数字照相技术的进步也促进了数字图像分析技术的发展，主要有两类不同的技术开始应用于粒状材料细观测量：基于相关分析的技术与基于图像分割的技术。前者主要通过比较试样在不同时刻的图像，分析颗粒位移，可得到试样内部的连续位移场[7-13]。该技

───────

基金项目：国家自然科学基金项目（50478016）；北京市自然科学 基金重点项目（8061003） 收稿日期：2007–11–02

1556 岩 土 工 程 学 报 2008年

术的应用研究发展较快，已经被广泛应用于粒状材料的细观测量。但这类技术无法获得粒状材料的组构特征。基于图像分割的技术通过边缘检测等算法分割图像中特征不同的成分，用于分析粒状材料的细观组构特征。Hryciw和Raschke[8]以及Guler等[9]采用这种技术实现了对砂土颗粒的识别。岳中琦等[14]将该技术应用于岩石试样表面不同矿物颗粒的区分，并在此基础上建立岩石试样的有限元模型进行数字模拟，显示了这类技术的广阔应用前景。但这种技术应用仍然较少，其原因一方面在于粒状材料本身组构特征的复杂性，另一方面则在于图像分割技术尚不够完善。目前的图像分割多通过边缘搜索实现[14]，对图像质量要求比较高，搜索速度较慢。模板匹配是另一类图像分割技术，特别适合于颗粒形状规则的二维粒状材料。本文采用模板匹配技术实现了一种新的颗粒识别方法，并编制了适用于理想二维粒状材料的细观组构特征分析程序（IPFA）。该程序具有图像增强、颗粒识别、接触搜索与组构分析等功能。

Gr(x,y)=Gt(x,y)=

255×h(x,y) h(x+1,y)MAXh(x0,y0) h(x0+1,y0)255×h(x,y+1) h(x,y)MAXh(x0,y0+1) h(x0,y0)

， (5) ， (6) 。 (7)

Gb(x,y)=

255×h(x,y) h(x,y+1)MAXh(x0,y0) h(x0,y0+1)

其中MAX( )表示针对整个数字图像取最大值。式（4）～（7）等式右侧的分母分别为整个数字图像在四个方向上的最大灰度梯度。为Heavyside函数。以G(x,y)表示Gl(x,y)，Gr(x,y)，Gt(x,y)与Gb(x,y)的最大值，按式（8）变换该像素点的灰度值。

， (min9(x,y)>200)， 223+h(x,y)/8

F(x,y)= h(x,y) ， (else)，(8)

h(x,y)/2+G(x,y) ， (max9(x,y)<200)。 式中 F(x,y)为变换后该像素点灰度值；max9(x,y)

与min9(x,y)则表示该像素点9邻域内各像素点灰度值的最大值和最小值。其中223和200是根据本文数字照片的特点给定的阈值。当数字图像特征不同时，该阈值需要重新设定。

图1给出了按上述方法进行图像增强的效果，其

、（b）分别为原始灰度图像与图像增强后中图1（a）

的效果。图1显示，上述图像增强方法能够有效消除

图像中的噪声以及颗粒表面的斑点。

1 图像增强

在试验所获得的原始数字图像中不可避免的包含一定噪声，影响颗粒识别的精度，因此在颗粒识别前需要增强颗粒本身的特征，消除噪声的影响。这一步骤通常被称作图像增强。本文首先采用式（1）所示的9邻域卷积方法对图像进行模糊处理。

g9(x,y)=∑∑ωijg(x+i,y+j) ， (1)

i= 1j= 11

1

其中，g(x,y)与g9(x,y)分别表示像素点(x,y)的原灰度值与模糊处理后的灰度值；ωij为加权系数，定义如下：

121 1

ωij= 242 。 (2) 16

121 位于颗粒边缘的像素点灰度值通常集中在灰度较小的一侧，因此有必要将这一区域的灰度差放大。本文采用式（3）所示的五次函数实现这一目的。

h(x,y)=255 255×[1 g9(x,y)/255] 。 (3)

5

颗粒表面本身还可能有颜色不同的斑点。这并非数字图像采集过程中产生的噪声，但同样会对颗粒识别造成一定影响，并且这些斑点无法通过上述方法消除。考虑到颗粒内部像素点的灰度梯度较颗粒边缘处像素点的灰度梯度小，本文基于灰度梯度阈值对这些斑点作了处理。按式（4）～（7）分别计算像素坐标为(x,y)的像素点在左、右、上、下4个方向上的归一化灰度梯度Gl(x,y)，Gr(x,y)，Gt(x,y)与Gb(x,y)。

255×h(x+1,y) h(x,y)

Gl(x,y)= ， (4)

MAXh(x0+1,y0) h(x0,y0)

图1 试样数字照片局部

Fig. 1 Local zone of the digital image of the specimen

第10期 张连卫，等. 基于数字图像的粒状材料细观组构特征分析技术 1557

2 颗粒识别

2.1 基本原理

针对经过增强后的粒状材料试样的数字图像，本文采用模板匹配技术识别其中具有规则形状的颗粒，即用模板图像T与从待识别图像S0中切取的与模板图像大小相同的图像S进行匹配，主要包括3个步骤：①建立用于匹配的模板图像T；②从待识别图像S0中切取待匹配图像S；③判断S与T是否匹配。图2给出了本文模板匹配过程的示意图。

mST=

∑F(x,y)×255 。 (10) F(x,y)

STT

显然，当mST=255时，图像S与模板T完全匹配；

本文在0到255之间设当mST=0时两者完全不匹配。

置域值mv，仅当mST>mv时认为两者匹配。其中域值mv需根据实际匹配效果调整。

综上，模板匹配过程可用图3所示的流程图表示。其中为减少计算量，本文在匹配成功后将该颗粒所占据的像素灰度值设为0，以便从图像S0中抹除已识别的颗粒。需要指出的是，图3所示的模板匹配流程仅针对一种倾角的颗粒，因此待该过程完成后需要变换颗粒的倾角，重新生成模板图像，重复图3所示的模板匹配过程，直到遍历所有可能的颗粒倾角。为提高颗粒识别效率，模板图像T中颗粒的初始倾角可预设

为根据图像S0估计得到的颗粒长轴优势方向的倾角。

图2 模板匹配过程示意图

Fig. 2 Diagrammatic sketch of the template matching process

模板图像T须反映颗粒的典型特征，故应根据待识别图像选定模板类型，确定模板参数。本文主要采用椭圆形模板以匹配椭圆颗粒。用于描述这种椭圆形颗粒几何信息的参数共5个，包括颗粒长轴与短轴尺寸、中心位置以及长轴方向的倾角。模板参数则包括颗粒长轴与短轴尺寸以及长轴方向，可根据图2（a）所示的待识别图像S0估计。图2（b）即为按照估计的颗粒长短轴尺寸，设颗粒倾角为0°建立的模板图像T，其中颗粒与背景所占据的像素的灰度值分别设为255和0。IPFA提供了手动匹配功能，可用于估计椭圆颗粒的长轴与短轴大小。当材料中包含多种粒径或形状的颗粒时，IPFA允许分别设定相应的模板参数。

获得模板图像T后，可将S0中的像素点(x,y)作为左上角点，按照模板图像T的大小从S0中切取用于匹配的图像S，如图2（c）所示。图2（a）中虚线矩形框显示了图2（c）在S0中的位置。待识别图像S0中一般包含有多个目标。本文通过变化S在S0中的位置，逐像素点匹配，以搜索整个原始图像S0中与设定模板参数相同的椭圆形颗粒。 2.2 模板匹配判据

为定量描述图像S与模板T之间的匹配程度，本文按式（9）生成图像ST，如图2（d）所示，其效果相当于透过模板图像T观察图像S。

FST(x,y)=FS(x,y) FT(x,y)/255 ， (9)

其中FS(x,y)，FT(x,y)与FST(x,y)分别为图像S，T以及ST中(x,y)处的像素点灰度值。模板T与图像S之间的匹配度则由式（10）计算。

图3 模板匹配算法流程图

Fig. 3 Flow chart of the template matching

2.3 颗粒识别效果

图4给出IPFA针对图1（a）进行颗粒识别的效果，依据模板参数绘制了识别出的椭圆形颗粒轮廓，表明IPFA可正确识别其中多数颗粒，剩余大孔隙为

小颗粒占据，受图像质量限制，目前还无法完全识别。

图4 椭圆形颗粒识别与接触搜索结果

Fig. 4 Results of particle identification and contact searching

1558 岩 土 工 程 学 报 2008年

本文还针对以上算法在理想条件下的颗粒识别效果作了统计分析。按照随机生成的椭圆颗粒参数（包括位置参数与形状参数）可建立仅包含一个理想椭圆颗粒的灰度图像S0，其中颗粒与背景所占据的像素点灰度分别设为255与0。本文针对大量这样建立的图像进行了颗粒识别，并将识别结果与S0的参数比较，以计算识别误差。

图5所示为IPFA对识别误差的统计分析结果，其中A为目标颗粒的实际长半轴尺寸，a为模板颗粒的长半轴大小，σ为识别误差的标准方差。图5中颗粒长轴倾角识别误差的标准方差σ以角度为单位，其纵坐标轴位于图5右侧；其余参数误差的标准方差以像素为单位，纵坐标轴位于左侧。

从图5可见，随a/A增大，颗粒识别精度逐渐提高，表现为标准方差σ逐渐减小。其中长半轴和短半轴尺寸的识别误差受a/A影响不大，σ均在0.3像素以内；颗粒中心位置则随a/A变化明显，并且当a/A大于0.8时σ均在1像素以内；颗粒长轴倾角的识别

并且a/A大于0.8时标准方差σ精度也受a/A的影响，

在3°以内。这样的识别精度已能够满足组构分析的基本要求。增加a/A虽然能够提高颗粒识别精度，但当a/A超过0.95时，部分颗粒无法识别，导致识别率降低。因此保持a/A在0.9～0.95之间较为合适。

为检查本文模板匹配的效率，本文还针对一幅包含几千个颗粒的实际图像进行了颗粒识别。该图像的像素尺寸为1536×2048，包含3种直径的椭圆形颗粒。本文采用IPFA对其中一种直径尺寸较大的颗粒进行了识别，其短半轴的像素尺寸分别为8.4和4.1，长短轴的尺寸之比为2∶1。所采用的处理器为AMD Opteron Processor 246（1.99 GHz）。IPFA程序可在20 min左右的时间内识别出约2600个颗粒。

此上述颗粒识别算法仅能识别粒径1 mm以上的颗粒。随数字照相技术分辨率水平提高，该算法识别能力将逐渐增强。此外，考虑到实际材料中的颗粒并非理想椭圆形，通过本文模板匹配算法进行颗粒识别可能出现颗粒嵌入的现象，但这种现象对长轴方向的分析影响不大。

3 接触搜索

接触法线的分布特征是造成诱发各向异性的重要因素。本文开发了接触搜索算法，用于在上述颗粒识别结果基础上判断相邻颗粒之间是否发生接触并搜索接触点位置与接触法线方向。

如图6所示，P1和P2为两个相邻的椭圆形颗粒，C1和C2分别为这两个颗粒边缘上的点，其像素坐标分别为(xC1,yC1)、(xC2,yC2)，且相应的外法线单位向量nC1与nC2满足nC1=-nC2。若将满足这种关系的点C1和C2称为一个点对，以dc表示该点对的距离，则由几何关系可知当C1和C2恰好是两个颗粒的接触点时，dc最小。因此，可通过求dc极值的方法搜索颗粒之间的接触点并给出相应的接触点法线方向。

图6 接触搜索原理示意图

Fig. 6 Schematic diagram of the contact searching method

考虑到颗粒识别的结果中可能存在颗粒相互嵌入的现象，为避免错误的将颗粒轮廓的交叉点作为接触点对，本文通过式（11）计算dc。

dc= ， (11)

其中，k为待定系数。当C1与C2同时位于另外一个颗粒的内部时，k= 1；否则k=1。这样可保证搜索到正确的接触点位置与接触法线方向。为判断两个颗粒是否发生接触，本文定义了域值dv，以dc<dv作为相邻颗粒发生接触的判据。

图4显示了针对图1（a）进行接触搜索的结果，

可看出IPFA可正确给出接触点位置与接触法线方向。

图5 不同初始模板长轴时的模板匹配精度

Fig. 5 Precision of template matching varied with initial long axis

of the elliptical template

4 组构分析

图7为根据对试样的数字照片颗粒识别与接触搜索结果，采用IPFA颗粒长轴与接触法线的定向性分别进行统计分析的结果。其中，试样由椭圆形颗粒组成，水平制样，与竖直方向沉积的效果相当。图7（a）显示了颗粒长轴倾角的分布情况，表明颗粒长轴基本沿水平方向；图7（b）给出了接触法线倾角的分布情

当颗粒长轴像素尺寸小于4时，用于计算匹配度mST的像素点数目过少，使得该算法无法给出正确的匹配结果。限于当前数字照相设备的技术水平，粒状材料数字相片的分辨率可达到10像素/mm左右，因

第10期 张连卫，等. 基于数字图像的粒状材料细观组构特征分析技术 1559

况，表明接触法线基本沿竖直方向。该结论符合当前对粒状材料细观组构特征的一般认识。

[5] ROSCOE K, ARTHUR J R F, JAMES R G. Determination of

strains in soils by X-ray method[J]. Civil Engineering, 1963, 58(685): 1009–1012.

[6] ODA M, KAZAMA H. Microstructure of shear bands and its

relation to the mechanisms of dilatancy and failure of dense granular soils[J]. Geotechnique, 1998, 48(4): 465–481. [7] DONOHOE G W, BOCABELLA M F. Segmentation and

object tracking for the microstructure analysis of soil[J]. IEEE, 1991, 1/2(2): 1215–1219.

[8] HRYCIW R D, RASCHKE S A. Development of computer

vision technique for in situ soil characterization[J]. Transportation Research Record 1526, Transportation Research Board, Washington D C, 1996: 86–97.

[9] GULER M, EDIL T B, BOSSCHER P J. Measurement of

particle movement in granular soils using image analysis[J]. Journal of Computing in Civil Engineering, ASCE, 1999, 13(2): 116–122.

[10] 孔宪宾, 何卫忠, 佘跃心. 测量土体内部位移场的一种新

方法[J]. 力学与实践, 2000, 22(3): 41–44. (KONG

图7 试样组构特征分析结果

Fig. 7 Results of microfabric analysis of the specimen

Xian-bin, HE Weizhong, SHE Yuexin. A new method in measuring the soil displacement field[J]. Mechanics and Engineering, 2000, 22(3): 41–44. (in Chinese))

[11] 张 嘎, 张建民, 梁东方. 土与结构接触面试验中的土颗

粒细观运动测量[J]. 岩土工程学报, 2005, 27(8): 903–907. (ZHANG Ga, ZHANG Jian-min, LIANG Dong-fang. Measurement of soil particle movement in soil-structure interface test[J]. Chinese Jounal of Geotechnical Engineering, 2005, 27(8): 907–907. (in Chinese))

[12] ZHANG G, LIANG D, ZHANG J M. Image analysis

measurement of soil particle movement during a soil-structure interface test[J]. Computers and Geotechnics, 2005, 33(4/5): 248–259.

[13] 李元海, 靖洪文, 朱合华, 马险峰. 数字照相测量在砂土

地基离心试验中的应用[J]. 岩土工程学报, 2006, 28(3): 306–311. (LI Yuan-hai, JING Hong-wen, ZHU He-hua, MA Xian-feng. Experimental investigation on progressive deformation patterns of sand foundation in centrifuge test using digital photogrammetry[J]. Chinese Journal of Geotec- hnical Engineering, 2006, 28(3): 306–311. (in Chinese)) [14] 岳中琦, 陈 沙, 郑 宏, 谭国涣. 岩土工程材料的数字

图像有限元分析[J]. 岩石力学与工程学报, 2004, 23(6): 889–897. (YUE Zhong-qi, CHEN Sha, ZHENG Hong, THAM L G. Digital image proceeding based on finite element method for geo-materials[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(6): 889–897. (in Chinese))

5 结 语

本文基于模板匹配技术提出了一种颗粒识别方法，并编制了相应的程序（IPFA），实现了对粒状材料细观组构特征的自动分析，可作为组构特征分析的有效工具。其特点在于：①适用于颗粒形状规则的理想二维粒状材料；②具有图像增强、颗粒识别、接触搜索和组构分析等多种功能；③颗粒识别具有较高的识别率，精度和效率均较高；④不需要破坏试样，能够连续跟踪粒状材料试样的细观组构演化。

参考文献：

[1] BREWER R. Fabric and mineral analysis of soils[M]. New

York: J Wiley, 1965.

[2] ODA M. Initial fabrics and their relations to mechanical

properties of granular materials[J]. Soils and Foundations, 1972, 12(1): 17–36.

[3] 张建民, 谢定义. 饱和砂土动本构理论研究进展[J]. 力学

进展, 1994, 24(2): 187–201. (ZHANG Jian-min, XIE Ding-yi. Advances on the research of dynamic constitutive theory of saturated sand[J]. Advances in Mechanics, 1994, 24(2): 187–201. (in Chinese))

[4] PARKIN A K, GERRARD C M, WILLOUSGHBY D R.

Discussion of deformation of sand in shear[J]. Journal of the Soil Mechanics and Foundations Division, ASCE, 1968, 94(S1): 336–340.

本文来源：https://www.bwwdw.com/article/mvb1.html