基于图像处理的混凝土骨料边界提取及有限元分析

基于图像处理的混凝土骨料边界提取及有限元分析

1混凝土细观非均值化的应用混凝土是由砂浆、骨料和空隙组成的多相化合物。在详细观察方面，我们可以观察混凝土的材料组成和力学性质的不均匀性。混凝土细观数值模型大致可以分为两种:(1)根据真实骨料形状及其分布规律应用统计学理论随机生成数值骨料模型;(2)通过对数字图像处理重构的混凝土细观结构。第一种模型建立在混凝土细观数值模拟研究的基础上,是由简单的几何体组合而成,且骨料位置随机分布,因此生成的模型与真实骨料结构之间或多或少存在差别,如假定骨料为圆形或球形,即使是任意多边形(多面体),也与真实骨料形状有一定差别,从而导致数值模拟的力学性能与真实结果有一定差距。第二种模型可以得到真实的骨料形状和分布情况,能够很好表征混凝土细观的非均值性。目前很多学者已就重构混凝土模型进行了研究,Buyukozturk通过对比CT方法、热红外线法、微波法和声发射法等,得出X射线的CT方法是研究混凝土内部结构的有效方法;Mora等研究了数字图像处理技术在骨料生成上的应用;John等发现数字图像关联术(DIC)适合于观察混凝土表面的小裂纹,X射线CT描述混凝土内部大裂纹更有效,并根据混凝土破裂后的CT图像讨论了骨料形状、裂纹形状对混凝土强度和韧性的影响;Yang等对由扫描电镜得到的照片进行处理,得到了从照片中分离出骨料元素的通用方法;Yue等发展了一种基于数字图像的有限元方法(DIP-FEM),并研究了沥青混凝土中集料形状和空间分布对混合料力学性质的影响;田威等利用CT图像信息研究了混凝土细观破坏过程;姜袁等、戚永乐等对二维CT图像处理,重构混凝土的二维、三维模型,并进行有限元分析。于庆磊等依据处理后的数字图片建立混凝土损伤数值模型,模拟了混凝土单轴载荷作用下的破坏过程。以上学者多以CT图像为研究对象。CT图像反映了材料的密度分布,混凝土材料中骨料、砂浆及空隙的密度差别较大,通过CT图像可以有效地分辨出各组成的分布,但CT图像的采样成本较高,对试件尺寸、形状和环境有限制,广泛应用成本较高。而由数码相机获得混凝土截面数字照片反映了材料的颜色分布,从照片中提取出混凝土各相信息成本低,操作简便,可以广泛使用。但混凝土材料中骨料相灰度值跨度大,砂浆相灰度介于其间,且骨料粒径分布不均,直接依据混凝土数字照片识别真实的混凝土细观结构有一定的难度,因此针对混凝土截面数字照片处理的研究具有非常重要的意义。本文以混凝土截面数字图(照)片为研究对象,应用数字图像技术和MATLAB图像处理工具针对二维混凝土骨料的边缘提取进行研究,然后依据图片处理结果建立混凝土细观层次力学模型,并对比随机骨料模型进行四点弯拉和轴拉断裂破坏数值分析。2混凝土在数字化处理时工作原理图1为二维混凝土截面的数字照片。首先从混凝土试件中选取一块100mm×100mm的混凝土截面做为本文的研究对象,图片的像素为1000×1000。以其中一块作为样本给出图像处理过程,样本的像素为460×400。从样本图片中可以看出骨料颜色深浅不一,同一块骨料亦由不同颜色的材料组成。砂浆颜色较单一,但是骨料灰度值分布于砂浆灰度值两端。在数字化处理时作如下假定:(1)将混凝土看作是由骨料和砂浆组成的二相材料,认为骨料和砂浆为单质材料,不考虑空隙和裂缝等;(2)处理过程中,主要考虑大骨料(骨料区域在图中的总像素大于100),忽略掉部分小骨料;(3)不考虑照片采样、量化过程引起的误差。依据以上假定在空间域上对混凝土截面照片进行边界拾取。由于混凝土骨料和砂浆内部纹理较多,且某些部位骨料和砂浆的边界不明显,因此首先对图片进行二值化、滤波等处理,以得到骨料和砂浆之间对比度高的图像,减少噪声等对边界提取的影响。2.1骨料和砂浆的二维荧光定量相关算法由于只考虑骨料和砂浆两相单质材料,所以首先对图像进行二值化处理,将骨料和砂浆直接以黑白区分。这样可以消除骨料和砂浆各自内部的纹理,同时可以减小图像的数据量,减少边界提取的复杂程度。首先绘制样本的直方图(如图2所示)。依据直方图可知,骨料的灰度值主要分布在55~100之间,砂浆灰度值主要为140~170。MATLAB提供了直接进行二值化的函数其中:f(i,j)为图像中第i行第j例像素的灰度值;ξ为骨料和砂浆灰度值的分界阈值。即按式(1)以两个波峰之间的波谷为分界阈值,将阈值两边的灰度值分别置0和1,处理结果如图3所示。由于骨料中含有灰度值较高的材料,因此直接二值化后骨料内部被空蚀较大,边界信息也造成很大的缺失,这使得提取的边界具有很大的误差,因此将图像按下式对图片进行分段变换其中:ξ1、ξ2为骨料和砂浆灰度值的分界阈值。确定骨料的灰度值区域后将大于和小于此灰度值区域的部分提取并设定为砂浆。由于骨料和砂浆的灰度值有互相侵入的现象,如果分段变换着重保留骨料信息,避免空蚀,则砂浆部分将引入更多的噪声。因此根据原图像和灰度直方图最终确定ξ1、ξ2分别为120和210,处理结果如图4所示。2.2中值滤波med由于机电噪声、拍摄环境以及上述分段变换等因素,生成的图片不可避免存在噪声的影响,因此需对分段变换结果进行滤波除噪,去除骨料和砂浆区域的杂质信息。经过上述处理后,骨料区域和砂浆区域为连续的单值块状区域,因此滤波过程要求还原两区域内部信息且平滑两者边界。本文主要考虑空间域内的滤波方法,按上述要求选取中值滤波法进行除噪。中值滤波的原理是将图像中某点灰度值用该点的领域(即单位处理区域)中各点值的中值代替。若设定单位处理区域为A,则中值滤波器输出为:其中:f(i,j)为滤波后第i行第j例像素的灰度值;Med{}为取中值函数;xij为数字图像中原来各点的灰度值;r、s为单位处理区域A的长和宽,单位为像素。中值滤波可解释为:若图像中以某点为中心的单位处理区域内骨料元素居多(像素个数超过半数)则令该点为骨料元素,砂浆元素亦然。本文采用方形的单位处理区域,并取边长L分别为4、10、16、24、30(单位:像素)对分段变换后的图像进行中值滤波(如图5所示)。从图5可以看出,L=4和L=10时中值滤波效果较差,砂浆区域仍存在大量噪声,且骨料部分仍有空蚀现象;L=30时中值滤波去除了砂浆区域多余信息,但是粒径较小的骨料在滤波后消失,而且在边界区域和骨料之间,骨料边界被扩大,识别效果较差;L=16和L=24时图像噪声得到了较好的去除,处理后骨料的大致范围已经确定,砂浆的边界得到平滑,且保留了较小粒径的骨料信息,但L=24时骨料与边界的连接范围增大,边界拾取时将产生误差,根据骨料的尺寸和内部空蚀情况等最终选取16×16像素的方形模块作为中值滤波的单位处理区域,处理结果如图5(c)所示。2.3骨料表面杂质目前研究者就CT图像的处理多进行了以上两步的处理,但从图4、图5可以看出,对混凝土截面照片进行以上处理后骨料内部仍存在杂质,且骨料边界区域较粗糙。这是由于骨料内部杂质较大且颜色与砂浆相近,分段变换后边界信息缺失较大,又由于中值滤波单位处理区域大小受到骨料粒径的限制。因此本文应用形态学理论继续对照片进行处理。2.3.1结构元素b编码基于膨胀腐蚀运算进一步对图片进行形态学处理,以消除图片边缘部分的多余信息。膨胀为对图像中的对象进行“加长”和“加粗”的操作;腐蚀则是“收缩”或“细化”图像中的对象。若用结构元素B对图像f进行膨胀操作,记作f⊕B,其集合运算定义为:f⊕B={z|||(B)̂z∩f≠∅}(4)其中:z为平移量;为结构元素B映像平移z后的集合。若用结构元素B对图像f进行腐蚀操作,记作f⊙B,其集合运算定义为:其中:(B)z为结构元素B平移z后的集合;fc为f的补集。直接应用膨胀或者腐蚀无法解决文中问题,且会在边界产生很大的误差,因此本文应用膨胀和腐蚀的组合操作进行处理。开运算将图像经结构元素腐蚀后再用结构元素对图像进行膨胀操作,记作fue049B,可以删除图像中不包含结构元素的对象区域,平滑对象的轮廓定义为:闭运算则是先进行膨胀操作再进行腐蚀操作,记作f·B,可以将狭窄的缺口连接起来,并填充比结构元素小的洞,定义为:2.3.2开运算和闭运算本文应用MATLAB中的IPT函数strel构造结构元素B,文中首先以“disk”构造半径为R的圆形结构元素,记作(“disk”,R),分别取半径R等于2、4、6、8(单位:像素)对滤波后图像先作开运算(如图6所示),再执行闭运算(如图7所示)。从图6、图7可以看出:当R=2时,开运算和闭运算对图像作用较小,图中仍存在多余信息;当R=8时,开运算和闭运算后图像中骨料之间的区域被连通,且骨料的形状被膨胀,与实际情况不符;当R=4、6时,开运算和闭运算结果较好,能够很好的提取出骨料的形状和范围,且去除了图像边界的多余信息。依据原图像中骨料的形状和大小,文中选取半径为6的“disk”构造最终结构元素。2.3.3骨料处理后的表面应力为说明结构元素形状对图像处理的影响,又分别选取(“line”,12,0)、(“square”,6)、(“diamond”,6)构造结构元素对滤波后图像进行开运算和闭运算,结果如图8所示。结合图7(c)、图8可以看出:(“line”,12,0)构造的结构元素对边界竖向的多余信息有很好的处理效果,但没有能够消除上下边界的多余信息,且小骨料处理后被拉伸;(“square”,6)、(“diamond”,6)同样能够消除图中的多余信息,提取骨料的形状,但处理后的小骨料分别趋近于正方形和菱形,在之后的力学分析中这些区域的边界容易形成应力集中,产生较大的误差;而(“disk”,6)则能得到很好的骨料边界形状,且简化了下一步的边界拾取。2.4图像边界提取经过上述处理后骨料区域和砂浆区域的多余信息已经基本去除,骨料的形状从原始照片中分离了出来,因此接下来对图像边界进行提取。选取处理边缘连续性方面较好的Canny边缘检测方法进行边界提取,处理结果如图9所示。边界清晰,数据连续且闭合,能够较好地表征混凝土细观结构。按照上述处理方法,本文对混凝土截面整体图像进行处理得到了混凝土截面的边界结果(图10)。3混凝土薄观察层机械分析3.1数值模型建立对混凝土截面数字照片进行处理后得到了骨料和砂浆的边界,其结果与真实结构较吻合,适合做进一步的力学分析。将前述边界结果导入商用软件ABAQUS,首先依据骨料边界数据生成骨料区域数值模型,然后设置一定厚度的界面层,最后将剩余部分生成砂浆模型,对整个界面区域网格划分如图11所示。为分析模型的力学特性,本文按图12、图13所示建立混凝土数值试件,并进行弯拉和轴拉数值分析。为减少计算量,只考虑混凝土试件的中间100mm长的区段为非均质区,其余部位看作均质各向同性的线弹性混凝土材料。3.2应变参数选取为研究混凝土模型受拉断裂过程,取应力-裂缝模型对混凝土的力学性能进行描述和模拟,并采用弧长法进行断裂过程分析。在混凝土试件的数值计算中若材料本构模型采用应力-应变关系,计算结果在一定程度上受计算网格大小的影响。Hillerborg提出了断裂能Gf的概念,即单位面积裂缝扩展需要的能量,并将其作为材料参数。由断裂能与裂缝宽度的关系,混凝土的软化性质就可以用拉应力-裂缝宽度关系曲线来代替应力-应变曲线从而缓解仿真结果对有限元网格尺寸的依赖。本文假定混凝土各相材料在达到峰值应力前为线弹性材料,峰值后以拉应力-裂缝宽度曲线来代替应力-应变全曲线的下降段(如图14所示)表示混凝土的应变软化行为。计算时选取表1中的数据做为混凝土各组分的力学参数,由于目前计算条件的限制,文中选取界面层厚度为0.2mm。图14中,ft为材料的抗拉强度,Wc为材料的极限张开位移。图15、图16是弯拉和轴拉力学分析后非均值区的破坏情况。从图中可以看出,混凝土试件受拉破坏时主要沿砂浆和骨料间的界面层裂开,最终形成一条贯穿试件的裂缝,破坏结果与实际试验分析结果相似。3.3骨料模型计算本文提出骨料面积级配来衡量骨料的的含量和分布情况,以生成与真实骨料截面相近的圆形骨料和任意凸多边形骨料试件。设骨料的截面面积为S,则定义面积级配P(S)为所有面积小于等于S的骨料的面积之和与骨料总面积的比例,表达式如下式所示其中:Sj为试件截面内骨料的面积;Si为面积小于等于S的骨料面积;m为面积小于等于Si的骨料的个数;n为骨料的总个数。以真实骨料的面积级配为参照,应用文献中随机骨料生成方法,分别生成了3组面积级配相近的圆形骨料和任意凸多边形骨料试件,同种类型骨料试件的面积级配相同而骨料位置不同。各类型骨料模型计算参数和面积级配曲线如表2和图17所示,图18和图19分别为试件圆形1和多边形1的非均值区。上述计算结果存在差异的主要原因有:(1)3种骨料模型的面积级配具有一定的差异;(2)随机骨料模型相对真实骨料模型骨料分布较均匀,承载时混凝土截面的应力分布则更加均匀,所以随机骨料的承载能力较大;(3)多边形骨料为各向异性,所以骨料位置分布不同时离散性大;(4)真实骨料的形状曲直不一,与多边形骨料具有较好的相似性,因此加载曲线多边形骨料更像真实骨料;(5)计算模型的网格的疏密和网格畸变对模型的加载结果也有一定的影响。4骨料分析结果本文应用数字图像处理技术和MATLAB图形工具对混凝土截面数字照片进行处理,通过图形分段变换、中值滤波、膨胀腐蚀等操作得到了骨料和砂浆之间的边界,对骨料元素和砂浆元素进行了提取。从提取结果看,本文提出的处理方法能够较好地从混凝土截面照片中提取骨料和砂浆元素,边界信息与原图像