杨莎莎,郝 龙,李 瑛,傅少君,贺昌海
(1. 延安大学 石油工程与环境工程学院,陕西 延安 716000; 2.西京学院 陕西省混凝土结构安全与耐久性重点实验室, 陕西 西安 710123; 3. 中铁二十局集团第六工程公司,陕西 西安 710032; 4. 武汉大学 土木建筑工程学院,湖北 武汉 430072; 5. 武汉大学 水利水电学院,湖北 武汉 430072)
混凝土结构裂缝的产生及延展严重危害建筑结构的健康与安全。大体积混凝土结构内部大量的水泥水化热使得结构内部温度升高,较大的内外温度差形成复杂的温度应力,最终导致裂缝产生;加上其他各种因素的作用,大体积混凝土的裂缝呈现分布范围大、形态不规则、数量繁多且交错复杂的特征[1-3],因此,识别难度大大增加。这些复杂裂缝如果没有被及时的发现和处理,极易发展成贯穿裂缝而对结构安全产生严重影响[4]。因此,快速、高效地检测大体积混凝土裂缝,并采取措施及时修补,对确保混凝土建筑物结构安全和正常运行有着十分重要的作用[5-9]。
人工检测是混凝土裂缝最常用的检测方法,主要依靠裂缝刻度尺和放大镜检测裂缝表面特征,该方法耗时费力,检测结果受人为因素影响大,精度明显不足[10],在震后鉴定中还具有一定的危险性[11]。随着计算机性能的不断提高,借助数字图像处理技术进行混凝土裂缝识别成为可能。亓大鹏[12]建立了分布式集群,实现了隧道裂缝图像分布式存储与分布式处理,解决了海量隧道裂缝图像存储与处理的问题;TONG Xuhang等[13]和XU Bing等[14]对比了基于自适应阈值和基于区域增长的图像分割算法,得到了图像中的裂缝宽度值;方志等[15]研究了采用正拍和斜拍图像检测混凝土裂缝的方法;刘娜等[16]运用形态学和最大熵图像分割方法来检测路面裂缝,较好地实现了裂缝与背景间的分离;QU Zhong等[17]提出的改进基于渗流模型的图像检测算法能精确检测到不清晰的裂缝区域;卫军等[18]采用遍历边缘点的最短距离法求裂缝宽度,但在求取最短距离时无效计算偏多;许薛军等[19]提出了基于斜率法的裂缝宽度检测技术,但该技术不适用于裂缝拐角过多的情况。
笔者基于数字图像处理技术,在MATLAB平台上开发了一个混凝土裂缝识别程序,可快速、准确检测裂缝的位置、宽度和长度,该程序可作为混凝土结构工程健康监测的一个新方法。
1.1.1 图像灰度化
对于由三基色红(R)、绿(G)、蓝(B)相互叠加的RGB彩色照片,R、G、B分量取值范围为0~255,每个像素点有256 × 256 × 256种颜色。RGB图像的每个像素点需用8个字节存放信息,而灰度图像的每个像素点只需1个字节存放信息。由于混凝土表面裂缝无须彩色信息,因此,为了减小处理计算工作量,提高运算效率,可将RGB彩色图像转换为灰度图像。
当RGB彩色图像的红、绿、蓝3个分量等值时即为灰度图像。目前,图像灰度化方法主要有最大值法、平均值法及加权平均值法。人眼通常对绿色最为敏感,红色次之,蓝色最弱。相对而言,加权平均值法算出的各像素灰度值较为合理。因此,笔者采用加权平均值法,即利用YUV空间与RGB空间的对应关系来对彩色图像进行灰度化。YUV色彩空间中的Y表示图像亮度,灰度图的相关信息包含在Y分量中,转换关系如式(1):
(1)
从而,可得灰度计算式(2):
Y=0.299R+ 0.587G+ 0.114B
(2)
1.1.2 灰度校正
由于自然光照不充足,或者相机缺陷,或者混凝土表面杂物等影响,灰度化图像存在失真问题。因此,可用数学方法来改变图像灰度值的分布,使其分布的更加均匀,从而改善图像的对比度。
一般认为灰度化的图像由噪声、背景及裂缝3部分组成。其中,裂缝的灰度值相对较小,背景的灰度值相对较大,噪声的灰度值介于二者之间。如果不对接近裂缝灰度值的背景及噪声区域进行处理,在后续操作中就会把背景或噪声处理成裂缝,使参数的计算结果出现较大误差。笔者采用线性变换方法进行校正:用一个线性函数来改变图像灰度值,当线性函数的斜率k> 1时,输出的图像对比度增大;当k≤1时,输出的图像对比度减小。对整个灰度范围进行线性变换的方法不利于将裂缝和背景的对比度很好地区分,而分段线性变换能较好地解决这一问题。分段线性变换的步骤如下:
Step 1计算图像的大小,建立矩阵来储存变换图像。
Step 2遍历原图像的所有像素点,统计其灰度值。
Step 3对于每一个像素点,用分段线性变换式(3)来校正灰度值:
(3)
含裂缝的混凝土结构表面,图像会受到多个因素影响而产生噪声。如:在成像过程中系统本身产生的噪声,日照不充足或不均匀产生的噪声,结构表面平整度差产生的噪声,结构表面附着杂物引起的噪声等。这些噪声在灰度图上表现为一系列孤立的像素点或像素块,导致干扰信息增加、图像质量降低,极有可能误将这些噪声识别为裂缝。
图像滤波处理就是消除或抑制存在于图像中的噪声信息,使裂缝特征信息更加明显,同时尽可能不改变图像的原始细节特征。经过滤波后的图像,裂缝和背景的界面更易区分,更有利于图像分割处理以及裂缝各个参数的准确计算。
通常用到的滤波方法有均值滤波法和中值滤波法。由于中值滤波法是一种基于像素点灰度值大小排序的非线性滤波方法,它能有效地消除图像中存在的孤立噪声点,所产生的图像细节模糊程度比均值滤波要少很多。因此,笔者采用中值滤波法进行滤波处理。中值滤波法的基本思路是:①选定合适的含有若干个像素点的模板窗口,将模板窗口中心位置放于待处理像素点;②遍历模板窗口包含的像素点,首先计算出对应的亮度值,然后按从小到大的顺序排列;③待处理像素点新的灰度值为排序最中间像素点的灰度值即窗口像素点个数为奇数时,待处理像素点的灰度值取序列中间的灰度值,而窗口像素点个数为偶数时,待处理像素点灰度值取序列中间两个像素点灰度值的平均。
如图1(a),设待处理像素点的原灰度值为22,其与周围相邻几个像素点的灰度值存在显著差异,可看成噪声点。采用3 × 3窗口选出包含待处理的9个像素点,并按灰度值从小到大排列:22、39、56、59、63、65、71、88、94,则待处理像素点灰度值取排列中间值63,如图1(b)。
图1 原灰度值及中值滤波后灰度值Fig. 1 Original gray value and median-filtered gray value
图2为将同一图像经过3×3、5×5、7×7不同模板窗口进行滤波处理的结果。可以看出:中值滤波的效果受模板窗口尺度影响明显,窗口过小不能很好地消除噪声,而窗口过大图像又变得模糊。与7×7窗口相比,3×3窗口的滤波效果较差,但模糊度较小。因此,在图像识别中应考虑选择合适的模板窗口大小。笔者采取先小尺度模板窗口处理,再逐渐增大模板窗口处理,可获得很好的效果。
图2 原灰度图像及3×3、5×5、7×7中值滤波后图像Fig. 2 Original grayscale image and 3×3, 5×5 and 7×7 image after median filtering
采用阈值分割方法进行图像分割处理,即从图像中提取裂缝。首先,确定图像分割的合理阈值t;然后,将阈值t与图像所有像素点进行比较;最后,区分裂缝区域和背景区域。经过图像分割处理后,灰度为0~255的原灰度图像变成了0~1的二值图像,从而大大减少了图像占用的内存,并且方便了裂缝参数的提取。阈值分割法的基本表达式如式(4):
(4)
阈值t的确定对图像分割效果至关重要。若t过大,则会将本不被纳入的背景像素点当成裂缝区域点;若t过小,则会将本不被纳入的裂缝像素点当成背景区域点。为解决这一问题,笔者采用了迭代法、自定义阈值法及Otsu法3种算法来确定阈值t。
1.3.1 迭代法
迭代法的基本思想是预先设置一个初始阈值t0,按照某种迭代规则不断循环直至最终阈值的误差小于预设的允许误差值。迭代法在计算过程无需人为设定参数,是一种自动求取阈值的算法,简单易操作、准确性较高。迭代法的基本流程如图3。
1.3.2 自定义阈值法
自定义阈值法,即提前预设1个阈值t,然后统计所有像素点的灰度值大小,将其逐一与设定阈值作比较,小于阈值的像素点归于裂缝区域,大于阈值的像素点归于背景区域。
图3 迭代法的基本流程Fig. 3 Basic flow of iterative algorithm
1.3.3 Otsu法
假定灰度值范围为[1,t-1]的像素组成C0区域,灰度范围为[t,m]的像素组成C1区域,原图像灰度值平均值为μ,则
C0出现概率σ0及均值μ0为
(5)
C1出现概率σ1及均值μ1为
(6)
C0与C1区域的总方差为
σ2=σ0(μ0-μ)2+σ1(μ1-μ)2
(7)
当σ2达到最大时,灰度值t就是最佳阈值。
笔者对某灰度化图像〔图4(a)〕,分别采用迭代法、自定义阈值法(阈值t= 128)、Otsu法进行分割处理,处理后的效果如图4(b)~(d)。可以看出,当t=128时,灰度图像分割效果最佳。
图4 图像分割效果Fig. 4 Image segmentation effect
对RGB彩色图像进行灰度化、滤波及分割处理后得到了二值灰度图像,但该二值灰度图像仍然不可避免会存在孤立点或孤立块等残余噪声,从而影响裂缝形态及参数的识别精度。针对混凝土结构表面裂缝特点,笔者采取以下方法来消除影响裂缝形态和参数识别的残余噪声:
1)用标记连通域的方法对裂缝对象进行逐个标记,得到3个参数:疑似裂缝的面积S、最小外接椭圆长短轴比β、裂缝区域与最小外接矩形区域的灰度平均值之差δ。由于裂缝为细长状,裂缝区域与非裂缝区域的3个参数相差显著,因而可通过设定这3个参数的限值将孤立点或孤立块等噪声消除。
2)通过对混凝土结构表面裂缝的统计分析,裂缝连通区域的最小外接椭圆的长短轴比值β基本上均大于6.5;裂缝区域的灰度平均值与最小外接矩形区域的灰度平均值之差δ> 25的情况占90%。因此,若疑似裂缝连通区域的最小外接椭圆的长短轴比值β< 6.5,或者疑似裂缝区域的灰度平均值与最小外接矩形区域的灰度平均值之差δ< 25,则可将该裂缝区域划分为背景区。
3)当完成上述裂缝判断后,有时会存在裂缝区域片段或不连续,则需对裂缝进行片段连接从而得到连续裂缝。笔者采用张国旗[21]提出的模型进行裂缝片段连接,模拟效果如图5。
图5 裂缝判断结果及断点连接效果Fig. 5 Crack judgment and connection effect of breakpoints
2.2.1 像素标定
由于图像要素的几何特征(如裂缝宽度、长度)是通过像素个数来度量的,因此须确定单个像素代表的实际长度大小,也就是单个像素尺寸的标定。
拍照前,先在拟拍照的混凝土结构表面刻画已知尺寸的参照标记,通过参照标记尺寸反推单个像素尺寸大小。设参照标记的实际宽度为s,对参照标记图像进行图像处理,可得到二值灰度图像,再求出参照标记宽度的像素个数a,即可得到单个像素代表的实际尺寸s/a。
2.2.2 裂缝面积S
对于二值灰度图,白像素点代表裂缝区域,灰度值为1;黑像素点代表背景区域,灰度值为0。设白像素点总数为N,则裂缝的实际面积S=N(s/a)2。
2.2.3 裂缝长度L
假设裂缝长度为裂缝两条边缘线长度的平均值,当裂缝两条边缘线的像素个数共计为M,则裂缝长度L=Ms/(2a)。
2.2.4 裂缝宽度D
裂缝宽度D可近似地看成裂缝面积与长度的比值,即D=2Ns/(Ma)。裂缝宽度计算示意如图6,计算步骤如下:
1)设裂缝图像由m行n列个像素组成,首先统计每列包含像素点个数B(j),同时记录对应的位置坐标〔C(k,j),j〕,j=1~n,k=1~B(j)。
3)第j列的裂缝实际宽度为D(j)=B(j)cosθjs/a。
4)裂缝的最大宽度即D(j)中的最大值。
3.1.1 程序总体框架
笔者利用MATLAB开发了基于数字图像的混凝土结构表面裂缝识别程序identification program of crack width(IPCW)。IPCW总框架如图7。
图7 IPCW总框架Fig. 7 Overall framework of IPCW
3.1.2 程序功能
IPCW的主要功能是:读入数码相机采集的混凝土结构表面裂缝图像;对RGB图像进行灰度化处理以增强对比度;对灰度化图像进行滤波降噪处理;分割处理为二值灰度图像;裂缝判断和断点连接处理;计算混凝土结构表面裂缝的面积、长度及最大宽度等参数。上述功能均为模块化,既可根据需要独立运行也可依序运行。
为了检验IPCW的可靠性及裂缝参数识别精度,笔者构建了如下算例:
1)用AutoCAD画出10条不同宽度的实线A,在实线旁绘制1个正方形块,并将其打印出来。
2)用显微镜读出实线的最大宽度,每条实线读取10次,取平均值作为实线的实测最大宽度,结果见表1。
3)将打印出的纸张固定在墙壁上,用固定在三脚架上、分辨率为2 448 × 3 264的数码相机进行拍摄获取RGB彩色照片。
4)运行IPCW,读入照片,依次进行图像处理,如图8。实线A最大宽度的识别结果见表1,而正方形块被IPCW当成噪声去除。
图8 RGB照片及图像处理Fig. 8 RGB photo and image processing
表1 实线A最大宽度实测值与IPCW计算值Table 1 The measured value of the maximum width of solid line A and the calculated value of IPCW
由表1可以看出,实线A(即裂缝)最大宽度IPCW识别值与实测值的相对误差基本随实测值的减小而增大,变化范围为1.342%~8.333%。原因是随着裂缝宽度的减小,一方面显微镜读数误差会变大,另一方面裂缝宽度方向的像素点偏少,由此引起的误差也会增大。
在我国西南地区,某混凝土坝工程竣工时,现场检查发现,在坝体下游表面存在一些较长的水平裂缝和不规则裂缝,经凿槽细查,确定均为浅表性裂缝。坝体下游表面裂缝宽度D的实测值与IPCW计算值见表2。
表2 裂缝宽度实测值与IPCW计算值Table 2 The measured value of crack width and the calculated value of IPCW
由表2可以看出,裂缝宽度IPCW计算值与实测值的相对误差总体在11%以内,精度满足工程要求,证明IPCW在实际工程大范围搜索裂缝中具有明显优势,适用于工程的裂缝初查。
采用数字图像处理方法,构建了混凝土结构裂缝识别算法;在MATLAB平台上开发了裂缝识别程序IPCW并进行了检验;通过工程实例验证了IPCW的可靠性及精度。研究结果表明:IPCW对大范围搜索裂缝具有较大优势,可应用于土木工程结构的裂缝初查。