岩土工程材料的数字图像有限元分析

1、第23卷 第6期岩石力学与工程学报 23(6:8898972004年3月 Chinese Journal of Rock Mechanics and Engineering March ,20042002年8月28日收到初稿,2002年10月19日收到修改稿。* 特别行政区政府研究基金委员会(CA99100,EG01,HKU7029102E和香港赛马会慈善基金资助项目。作者 岳中琦 简介:男,40岁,博士,1983年毕业于北京大学,现任副教授、博士生导师,主要从事岩土工程、道路工程、滑坡调查、斜坡加固和应用力学等方面的研究工作。E-mail :yueqzqhkucc.hku.hk 。岩土工程材

2、料的数字图像有限元分析*岳中琦 陈 沙 郑 宏 谭国涣(香港大学土木工程系 香港摘要 综合数字图像处理理论、几何矢量转换技术与有限元网格自动生成原理,提出了岩土工程材料的数字图像有限元分析方法。以岩土工程材料的图像为研究对象,先采用数字图像算法获得材料的真实细观结构;然后,通过几何矢量转换技术将二元图像的细观结构转化为矢量化的细观结构;在矢量化细观结构的基础上,再利用网格自动生成技术生成材料细观结构的有限元网格;最后,采用传统的有限元计算理论,对细观结构中的不同材料赋予相应的材料参数,从而实现了岩土工程材料的非均质力学分析,真实地反映了岩土工程材料的力学性能。以花岗岩为例,详细阐明了实现数字图

3、像数值分析方法的整个过程。Brazilian 试验的数值模拟表明,材料的细观结构对Brazilian 圆盘中轴线上拉应力的分布有着明显的影响,数字图像数值分析方法可真实地实现材料的非均质研究。 关键词 数值分析,岩土工程材料,数字图像,细观结构,几何矢量化分类号 O 241 文献标识码 A 文章编号 1000-6915(200406-0889-09DIGITAL IMAGE PROCEEDING BASED ON FINITE ELEMENT METHODFOR GEOMATERIALSYue Z Q ,Chen S ,Zheng H ,Tham L G(Department of Civil

4、 Engineering ,The University of Hong Kong , Hong Kong ,China Abstract This paper presents an innovative digital image processing method based on a method to examine the mechanical responses of geomaterials by taking into account the actual material inhomogeneity. The proposed method incorporates the

5、 digital image processing procedures ,geometry vectorization algorithms and automatic finite element mesh generation techniques into the conventional finite element methods. Using a granitic rock as an example, the paper illustrates the innovative numerical method for mechanical analyses of rock mat

6、erials. Numerical results presented in the paper show that this method can be used to calculate the mechanical responses of inhomogeneous geomaterials with high efficiency and accuracy.Key words numerical analysis ,geomaterial ,digital image ,microstructure ,geometry vectorization1 概 述土、岩石、沥青混凝土和水泥混

7、凝土等岩土工程材料是天然材料或天然材料所组成的非均质材料,它们常常由孔洞、砂砾和石头等不同物质组成。这些物质的物理力学性质各不相同,当岩土工程材料受到外力作用时,它们对外力的响应有着很大的差异,同时,它们之间的相互作用也非常复杂。显然,岩土工程材料的力学性能,如应力分布、裂隙扩展、破坏模式等与材料的非均质性及细观结构有着密切的联系。 890 岩石力学与工程学报 2004年文献调查表明,目前大多数岩土工程材料的研究都是在传统表象理论基础上开展的宏观分析,众多力学理论和数值模型忽略了材料的细观结构,在假定材料为均质或分片均质的基础上,岩土工程材料的传统力学性能受到了较多的重视,而其细观力学性质却被

8、大大忽略。最近,文18试图开展岩土工程材料的非均质性数值分析,这些研究主要是基于在统计学层次上生成的虚拟细观结构,使生成的细观结构中各种材料的形状和分布与材料的真实情形尽可能一致。虽然以上研究得到了岩土工程材料的一些细观力学性质,然而,似乎很少有文献来说明这类虚拟细观结构能在多大程度上真实地体现岩土工程材料的非均质性。笔者认为,建立在虚拟细观结构上的力学分析是很难得到材料的真实力学性质的。要研究非均质的岩土工程材料,真实地反映材料的细观力学性能,需要开发新的数值计算方法。以往难以对材料实现真实的非均质性研究是由于缺少有效的设备与工具来观测材料的内部结构,以期得到其内部不同性质材料的真实几何分布

9、。由于岩土工程材料在转换为图像时,其内部不同物质可以通过灰色度的变化在图像中得到再现,图像很好地反映了材料的细观结构,保留了材料的非均质性。所以,数字图像处理方法为研究岩土工程材料的力学性能提供了一种新途径。数字图像处理是一门系统地研究数字图像理论、技术和应用的新学科,其主要过程是将研究的物体转化为储存在计算机中的数字图像,并运用计算机对图像信息进行分析和处理从而得到所需要的研究结果。文912表明,数字图像处理已经在岩土工程领域,如道路状况评估、混凝土细观结构定量观察中得到运用;文13,14应用数字图像处理方法成功地对沥青混凝土中骨料的方位、分布和形状进行了定量的分析;文15,16应用数字图像

10、处理技术分析了水泥混凝土中粗骨料形状分布,得到了骨料的厚薄率与延伸率;文5运用二次复型技术,采用光学显微镜及扫描电镜,研究了香港地区的花岗岩材料的细观破裂模式。目前,在岩土工程领域,数字图像处理主要应用于定量的、几何上的统计分析与观测,对材料的非均质性力学研究尚未得到实现。本文将综合数字图像处理方法、几何矢量转换技术与有限元网格自动生成原理,提出基于数字图像的有限元力学分析方法。本方法可用于研究岩土工程材料内部细观结构对材料力学响应的影响。本文的数字图像数值分析方法是基于二维空间的,但根据立体逻辑变换原理17,可将其应用到三维空间的有限元分析。2 图像数字化先用圆锯将现场得到的试样或试验室制成

11、的试块锯成垂直或水平的光滑横截面,然后通过数码相机、普通照相机或扫描仪将所研究的试件转换为储存在计算机中的物理图像文件。可将一测量尺放置在试块横截面下,并将试块的实际尺寸拍摄纪录在图像之中。就其自然格式而言,计算机中的物理图像文件是不能直接用来加工与处理的(用已有的图像软件则不属于此类,它必须先由图形格式转化为数字形式。数字图像在计算机中是由矩形排列的一个个图像元素,亦称像素点组成的。每个像素点是横向和纵向的扫描线组成的交叉区域,这些扫描线均有相同的宽度h。每个像素点有一对应的整数值来代表该像素点的亮度,亦称为该像素点的灰色度。对于常见的256色或二值图像,灰色度值分别为0255和01。实际上

12、,整个图像由有着不同灰色度的像素点阵组成,这个像素点阵的灰色度构成了一个离散函数f (i,j(i,j代表着笛卡尔系统中的位置,i = 1N,j = 1M:=(2(1(2(22(12(1(21(11(MNfNfNfMfffMfffjif,LMMMLL(1图1(a为岩石试块横截面的图像。岩石断面的实际直径为52 mm。它由816条横向扫描线和816条纵向扫描线组成。每个扫描线宽度为0.063 7 mm,每个像素点的实际面积是0.004 mm2。图1(b为部分数字图像的离散函数f (i,j示例。图像可以储存为很多种物理格式,如BMP, JPG,GIF等。根据不同的储存格式,编写了程序从物理图像中提取

13、出离散数据函数f (i,j,这些数据即是进行下一步数字图像处理的基础。3 材料的细观结构图1(a为香港地区典型花岗岩的试样图像,其主要由石英、长石和黑云母组成。3种矿物有不同的力学性能,它们的弹性模量分别为90,70和20第23卷 第6期 岳中琦等. 岩土工程材料的数字图像有限元分析 891 (a 岩石试块横截面图像 (b 局部坐标下部分离散函数f (i ,j 示例图1 香港地区典型花岗岩图像Fig.1 Typical digital image for cylindrical granite of Hong Kong GPa 。由于石英和长石的弹性模量相近并比黑云母的弹性模量大很多,因此,可

14、以将黑云母看成一种物质,而将其他的矿物看成另一种物质来进行力学分析。得到花岗岩内部材料的真实分布是实现数字图像有限元计算的第1步,在本节中,将采用2种不同的算法来提取材料的真实细观结构,得到的细观结构实际上是内部材料分界面二值图(其中,灰色度为0的像素点为分界面点。 3.1 区域分割算法从图1(b的离散函数f (i ,j 中可以看到,黑云母的灰色度值比较低,而其他矿物的灰色度值比较高。图2为图1(a的灰色度直方图,其横坐标为0255,代表着图像的256个灰色度。对于每一灰色度,相应纵坐标上的数字表示在整个图像中有多少个像素点拥有该灰色度值。图2的灰色度直方图曲线上有一明显的波谷,从相应的灰色度

15、位置(灰色度值为100左右,可以将整个图像分成2个部分。根据下面的算式,得到一个二值图像:<=H j i f H j i f j i g ( ( 01(, (2式中:H 为一整数阈值(对应波谷处的灰色度值。H 值的选定影响着花岗岩的细观结构,如果阈图2 图1(a灰色度直方图Fig.2 Histogram of image gray level in fig.1 (a值H 设置太低,一部分黑云母将被认定为其他矿物;相反,如果阈值H 设定太高,黑云母成分可能会被扩大。因此,可以逐步尝试直到得到一理想的结果。图3为设定H 为100后得到的二值图,中间部分灰色度为0的黑色物体代表着黑云母矿物。与

16、图1(a比较,可以看出,图3正确地检测出了大部分的黑云母矿物。在现阶段,可以只考虑较大面积的黑云母,而将图3中一些小黑色物体转换成其他矿物。如前面所说,图像中每一个像素点的图像面积为1,它相当于h 2的实际面积。每个黑色物体的面积可以通过计算它的整个像素点数而得到:io69696764647189113 131 138 141 14772757065667594119 138 144 144 151677269707788107132 152 159 157 16273778088100115135160 180 177 174 1788990103109119143161172 187 18

17、2 187 194120120136146156173182185 194 195 197 199146149164176184195199199 203 202 202 200153164175183188195201206 210 205 204 200161178181186188190197205 207 207 207 202178193192195195192196200 196 201 203 199198208203206205199201202 193 195 199 195211218210212212206209211 200 199 203 196ojii 892 岩石

18、力学与工程学报 2004年 图3 区域分割算法结果(中间黑色区域代表黑云母矿物 Fig.3 Result of region growing method (black areas represent biotite2nhA=(3 式中:n为每一黑色物体内的总的像素点数。设定一临界值来保留面积较大的黑云母,而去掉较小块部分。比如设定临界值n = 37,由于本例中每一像素点代表着0.004 mm2的实际面积,相应的实际临界面积为0.148 mm2。利用式(3将图3中小块的黑云母转换成其他矿物后,由于图3为二值图,可以很容易地找出黑云母与其他材料的分界面,从而得到花岗岩真实的细观结构(图4。与图1

19、(a比较,图4真实地反映了花岗岩的内部材料分布。 图4 花岗岩真实的细观结构图Fig.4 Real microstructure of granite3.2边界检测算法不同的材料在图像中往往有不同的灰色度,其分界面处的灰色度值变化也很大。一般而言,材料的分界面处灰色度的1阶导数往往是一较大数值。根据这一原理,可以通过计算灰色度的导数而检测材料的分界面,这种算法在数字图像理论中称为边界检测算法。点(i,j处灰色度的1阶导数可以通过下面数学公式计算:1(1(21jifjiffifi,+=(4a1(1(21+=jifjiffjfj, (4b 式(4a和(4b分别代表i和j方向的1阶导数,也可用2个方

20、向的1阶导数矢量和来综合考虑:22+=jfiff(5 式(4仅考虑了相邻的2个像素点,对于数字图像,实际上可以采用4个相邻的像素点来计算,相应的数学公式为+=11(11(41jifjiffifi,11(11(+jifjif,(6a+=11(11(41jifjiffjfj,11(11(+jifjif,(6b亦可综合式(4和(6来计算i和j方向的1阶导数:+=1(1(2121jifjiffifi,+11(41jif,+11(jif,+11(11(jifjif,(7a+=1(1(2121jifjiffjfj,+11(11(41jifjif,+11(11(jifjif,(7b 当灰色度的1阶导数计算出

21、之后,可以设定一阈值来找出不同材料的分界面。如果像素点的1阶导数较大,认为这个像素点为材料的分界面,计算后产生的边界图像为第23卷 第6期 岳中琦等. 岩土工程材料的数字图像有限元分析 893 |(| |(| 10(K f K f j i g <=, (8式中:K 为一非负的阈值。图5为采用边界检测算法后得到的内部材料分界面图。必须指出的是,由于图像噪音的大量存在以及在图像数字化过程中的一些不可避免的因素,如光线的不均匀照射等,要实现完美无暇的边界检测是非常困难的。在自动检测的基础上可用已有图像处理软件(photoshop,对分界面图像中不完善的地方进行加工,去掉多余的线段,闭合分界面线

22、,最终形成的分界图像如图6所示。通过对图1和6的比较,可以看出,经过很小的加工,花岗岩的细观结构得到了真实的反映。由于采用了不同的算法,图6的细观结构与图4略有不同。为了简便起见, 图5 边界检测算法结果 Fig.5 Result of edge detection method 图6 边界检测修正结果Fig.6 Revised image of edge detection result将采用图4的花岗岩细观结构作为进一步分析的基础。以上介绍了采用数字图像技术得到岩土工程材料细观结构的2种方法,实际上,在数字图像理论中还存在其他算法来获取材料的细观结构。以边界检测方法为例,除了应用1阶导数外

23、,还可应用2阶导数检测边界,很多具体的算法,如文18,19算法等也是检测边界的有效方法。由于图像的多样性,对于不同的图像,常常需要用不同的数字图像算法来获得其细观结构。而在数字图像自动检测的基础上,最终的人工修正亦是不可避免的过程。4 边界图形几何矢量化4.1 几何矢量转换算法检测后的分界面图像,是一个二元图。闭合分界面就是由灰色度为0的像素点互相连接而成。为实现下一步的力学分析,必须生成有限元网格。从分界面图像是无法来直接生成有限元网格的,原因在于:(1 目前,没有任何有限元软件可接受图像作为输入的数据;(2 若将闭合分界面的数据按一个一个像素点输入,则根据有限元网格的生成原理,每一个像素将

24、成为一个有限元单元的节点,这将导致网格生成非常困难和进一步导致有限元力学分析的失败。因此,必须把离散的闭合分界面像素点数据转换为几何矢量信息,才能实现有限元网格的生成。仔细观察分界面的形状,可以看到大多数呈多边形状,只有少数分界面是圆形的,但如果多边形的边数足够多,那么用其来代替圆形的分界面的误差是可以忽略的。因此,可用矢量多边形数据来代替图像中离散的闭合分界面结构,构造几何矢量转换算法如下:(1 设定一个阈值t ;(2 找出边界中相距最远的2个像素点; (3 连接这2个像素点,则该直线将闭合分界面分成2个部分,先考虑任一部分;(4 求出该部分(区间所有像素点到分割线的距离,记录下最大距离及相

25、应的像素点;(5 若最大距离小于阈值t ,则认为该段直线可以代表这一部分(区间的边界,可考虑另一部分;(6 若最大距离大于阈值t ,则将记录下的像素点与前分割线的2个端点连线,将这一部分(区间 894 岩石力学与工程学报 2004 年 再分为 2 个小区间，重复(4，(5的过程； (7 不断循环，直到计算出每个区间内像素点 到分割该区间的直线距离均小于阈值 t 为止，那么 此时形成的闭合的分割线便是代表分界面图像的多 边形。 设定不同阈值 t，将得到不同的多边形，阈值 t 越小，则多边形与分界面图像越相近。图 7 和 8 为 设置不同阈值的二值图像细观结构的矢量转换结 果。可以看到，当 t =

26、 3 时，矢量多边形已经非常接 近于实际的细观结构图像。 边长除以该边长在图像中的长度。将几何矢量化后 的多边形参数(纵横坐标乘以系数 h，则可以将图像 尺寸转换为实际试块尺寸，这个转变对下一步的有 限元计算是非常重要的。 得到的多边形数据可以直接用来生成有限元 网格或者根据不同格式输入到已有的有限元软件 中。图形交换文件 DXF 是 AutoCAD 的标准 ASCII 文件，由于其已被大多数的有限元软件采纳为输入 及输出数据的文件格式，故又编写了 Autolisp 程序， 根据多边形的数据可以自动地将矢量细观结构读入 并显示在 AutoCAD 中，将其导出，存为 DXF 文件 格式，这样就建

27、立了边界矢量数据与有限元软件的 数据输入接口。 5 有限元网格自动生成 图7 Fig.7 of 3 阈值为 3 的矢量转换结果 Geometry vector transformation result with a threshold 有限元网格的自动生成是一个比较成熟的研 究课题，主要有 Octree 算法20、Delaunay 算法21、 Advancing Front 算法22等。本文采用 Delaunay 算法 来自动生成有限元网格，其适用于生成二维与三维 的三角形与四边形网格，生成过程是： (1 先产生边界上的网格顶点； (2 再根据 Delaunay 判断准则(任何网格顶点 的外

28、接圆只能包含 3 个网格顶点来生成区域内有 限元网格。 编写了计算机程序读入包含了矢量数据的 DXF 文件，再运用 Delaunay 算法来自动生成有限 元网格。图 9 为自动生成的有限元网格，其总共包 括了 16 361 个三角形单元和 8 327 个网格顶点。可 以看出，网格生成得非常均匀，图 9 的有限元网格 可作为下一步的力学计算的基础。 图8 Fig.8 of 1 阈值为 1 的矢量转换结果 Geometry vector transformation result with a threshold 4.2 比例变换及有限元软件接口的实现 当花岗岩试块切面转换为图像时，试块的尺寸 也

29、转变为图像的尺寸，转换比例 h 即为实际试块的 Fig.9 图9 自动生成的有限元网格 Automatically generated finite element meshes 第 23 卷 第6期 岳中琦等. 岩土工程材料的数字图像有限元分析 895 6 应 用 式中： 为垂直于 AB 轴应力， 为平行于 AB 轴 应力，F 为施加的荷载，r0 为圆盘的半径，t 为圆盘 厚度，r 为离圆盘中心的距离，2 为力 F 在圆周上 均匀分布的角度。 在式(9a中，取 r 为 0，则可得 F sin 2 F = 0 1 r0t r0t F 。 r0t 本文采用 Brazilian 圆盘试验来进行数字

30、有限元 数值模拟计算。Brazilian 圆盘试验是测量岩石抗拉 强度的室内试验，如图 10 所示，对称荷载 F 均匀 地施加在圆盘上(AB，假定圆盘材料为均质线弹性， 则可得到轴线 AB 上的应力理论解23(假定压应力为 正为 = (10 式中： 0 = 定义无量纲化的应力为 0 = 0 = o (11a (11b 岩石材料的抗拉强度即可采用式(10来计算。 以图 9 的有限元网格为基础，采用数字图像数 值分析方法模拟均质材料和非均质材料的 Brazilian 试验。假定花岗岩由均质材料组成，材料参数为： E = 6 GPa 和 = 0.2，1 kPa 均匀力作用在圆周上， 其在圆周上的分布

31、角度为 3° 。 圆盘中心 x 和 y 方向无量纲的应力 和 的 解析解分别为0.998 和 2.996，数字图像有限元的 和 计算结果分别为0.099 4 和 2.999；计算结 果与理论解吻合得非常好，最大误差仅 0.4%。 考虑到花岗岩材料非均质性，取黑云母的材料 参数为： E = 2 GPa 和 = 0.25，其他矿物的材料参 数为： E = 8 GPa 和 = 0.25 进行计算。在分辨不同 材料区域时，对于每一种不同的材料，已给其一个 代号。在有限元分析时，根据不同的代号便可以指 定相应的材料参数，再利用图 9 的有限元网格，从 而实现了岩土工程材料的非均质性分析。图 1

32、114 为均质与非均质材料的 Brazilian 圆盘中最大和最小 的主应力(无量纲式(11等值曲线图。可以看到，材 料的非均匀性明显影响了在均匀材料中是对称分布 的应力曲线，并且在不同材料的界面处，出现了大 量应力集中现象。 图 10 Brazilian 试验示意图 Fig.10 Sketch diagram of Brazilian tension test r 2 1 sin 2 r0 F = 2 4 r0t r r 1 2 r cos 2 + r 0 0 r 2 1 + r 0 tan 1 tan 2 1 r r0 2 r 1 sin 2 r0 F + = 2 4 r0t r r 1

33、2 r cos 2 + r 0 0 r 2 1 + r 0 tan 1 tan 2 1 r r0 (9a 7 结束语 综合数字图像处理理论、几何矢量转换技术与 有限元网格自动生成原理，本文提出了岩土工程材 (9b 料数字图像的数值分析方法。应用数字图像算法， 首先，通过数字图像数值分析方法得到非均质材料 的真实细观结构；然后，采用几何矢量转换技术和 896 岩石力学与工程学报 2004 年 图 11 Fig.11 均质材料花岗岩最小主应力等值曲线 图 13 均质材料花岗岩最大主应力等值曲线 Contours of the minor principal stress in homogeneou

34、s granite cross-section under Brazilian indirect tension test Fig.13 Contours of the major principal stress in homogeneous granite cross-section under Brazilian indirect tension test 图 12 非均质材料花岗岩最小主应力等值曲线 图 14 非均质材料花岗岩最大主应力等值曲线 Fig.12 Contours of the minor principal stress in inhomogeneous granite

35、cross-section under Brazilian indirect tension test Fig.14 Contours of the major principal stress in inhomogeneous granite cross-section under Brazilian indirect tension test 参 考 文 献 网格自动生成方法，自动地形成了材料的矢量细观 结构及细观结构的有限元网格；最后，与传统有限 元分析相结合，实现了岩土工程材料的非均质分析， 真实地反映了非均质材料的力学性能。Brazilian 计 算实例表明，本文采用的数字图像数值分

36、析方法计 算正确，具有较强的优越性，是一种值得推广的新 方法。 在本文二维数字数值分析方法的基础上，亦可 通过对多层二维图像的分析，利用立体逻辑变换原 理， 实现三维空间有限元的计算。 作者将在以后的文 章中详细介绍三维数字图像数值分析的实现过程。 3 2 1 Sepher K，Svec O J，Yue Z Q，et al. Finite element modelling of asphalt concrete microstructureA. In ： Localized Damge ： Computer- Aided Assessment and ControlC. Southampto

37、n，Boston： Computational Mechanics Publications，225232 Wang Z M，Kwan A K H，Chan H C. Mesoscopic study of concrete ： generation of random aggregate structure and finite element meshJ. Computers & Structure，1999，70：533544 Tang C A，Liu H，Lee P K K，et al. Numerical tests on micro-macro relationship o

38、f rock failure under uniaxial compressionpart I： effect of heterogeneityJ. International Journal of Rock Mechanics and Mining Science，2000，37(4：555569 第 23 卷 4 第6期 岳中琦等. 岩土工程材料的数字图像有限元分析 5360 14 897 Tang C A，Tham L G，Lee P K K，et al. Numerical tests on micro-macro relationship of rock failure under

39、uniaxial compression part II： constraint， slenderness and size effectsJ. International Journal of Rock Mechanics and Mining Science，2000，37(4：571583 Yue Z Q，Morin I. Digital image processing for aggregate orientation in asphalt concrete mixturesJ. Canadian Journal of Civil Engineering，1996，23：479489

40、 5 李 炼，李启光，徐 钺等. 一种用于岩石材料微裂纹观察的复型 技术J. 岩石力学与工程学报，2002，21(6：797802 15 Kwan A K H，Mora C F，Chan H C. Particle shape analysis of coarse aggregate using digital image processingJ. Cement & Concrete Research，1999，29：1 4031 401 6 李 炼，徐 钺，李启光等. 花岗岩渐进破坏的微观研究J. 岩石 16 力学与工程学报，2002，21(7：940947 7 Li L. Micr

41、oscopic study and numerical simulation of the failure process of granitePhD DissertationD. Hong Kong：Department of Civil Engineering，The University of Hong Kong，2001 8 Li L， Tham L G， Tsui Y， al. Numerical investigation on the failure et process of Hong Kong graniteA. In：Proceedings of the 2001 ISRM

42、 International SymposiumSecond Asian Rock Mechanics SymposiumC. Beijing：s. n.，2001，303305 9 Frost J D ， Wright J R. Techniques and applications in civil engineeringA. In：Proceedings of EF/NSF Conference on Digital Image ProcessingC. Hawaii：A.S.C.E.，1993，139146 10 Lee H，Chou E. Survey of image proces

43、sing applications in civil engineeringA. In：Proceedings of EF/NSF Conference on Digital Image Processing：Techniques and Applications in Civil EngineeringC. Hawaii：A.S.C.E.，1993，203210 11 Chermant J L. Why automatic image analysis? an introduction to this issueJ. Cement & Concrete Composites，2001

44、，23：127131 12 Chermant J L，Chermant L， Coster M，et al. Some fields of applications of automatic image analysis in civil engineeringJ. Cement & Concrete Composites，2001，23：157169 13 Yue Z Q，Bekking W，Morin I. Application of digital image processing to quantitative study of asphalt concrete microstructureR. Mora