基于数字图像灰度相关性的类岩石材料损伤分形特征研究

邹 飞，李海波，周青春，朱小明，莫振泽，何恩光，赵 羽,

(1. 中国科学院武汉岩土力学研究所 岩土力学与工程国家重点实验室，武汉 430071；2. 贵州省质安交通工程监控检测中心有限责任公司，贵阳 550000；3. 北方重工盾构机公司，沈阳 110041；4. 浙江大学 现代制造工程研究所，杭州 310027)

1 引 言

岩石的损伤和破坏是一个渐进的过程，往往表现为微裂纹的萌生、扩展和贯通，通过微观测量来定量刻画岩石损伤演化所导致的力学性质的劣化是十分必要的。数字散斑相关方法作为一种全场非接触式测量手段，其关键在于所选择基准图像与变形后图像的相关性分析。该方法也被众多研究者[1－7]用于岩石类材料损伤、破坏机制的研究中。马少鹏等[1]采用数字图像灰度相关性分布的标准差为统计指标表征岩石变形和破裂的发展，并通过试验证明灰度性的分布演化与结构的变形局部化演化在空间和时间上都存在对应关系，该方法介于定性描述与定量分析之间。赵永红等[4]采用数字图像相关技术研究了含微裂纹岩石的变形，用相关技术处理灰度分布图获得位移场，对加载和卸载试件表面的灰度分布图进行相关计算，获得位移分布与裂纹分布的关系。朱珍德等[5]通过对红砂岩试样应变响应和表面裂纹图像的同步观测试验以及试样表面细观裂纹萌生、扩展图像信息的采集和量化，提出了用数字图像相关技术处理不同含水状态下红砂岩试样动态变形过程中的灰度分布图，探讨了红砂岩动态损伤状态与其宏观力学之间的关联。在对细观结构劣化的描述中，很多学者采用简单分形维和多重分形谱来进行刻画[8－13]。彭瑞东[8]通过灰岩拉伸过程中表面细观形貌的SEM数字图像，采用分形维数描述损伤的发展过程。曹茂森[9,13]采用多重分形谱理论研究了混凝土结构表面裂纹的分形特征，通过分形特征量来衡量结构损伤的程度。

本文以石膏试件的单轴压缩试验为基础，通过试件表面数字图像灰度相关系数的变化来表征试件表面损伤状态的演化过程。对不同时刻的灰度相关系数图进行图像处理，采用分形维数和多重分形谱定量描述材料表面细观损伤状态在时间和空间上的演化过程。

2 数字图像原理

2.1 数字图像灰度的相关性

图像可以定义为一个二维函数f(x, y)，其中x和y是平面坐标，f是图像在该点的振幅，也称为该点处的亮度。将图像转化为数字形式的结果是获得一个实数矩阵。对一幅图像f(x, y)进行取样后，得到m行n列的数组，这幅图像大小即为m×n，数组的每一个元素称为像素。灰度则是指黑白图像中点的颜色深度，范围为0～255，白色为255，黑色为0。

数字图像的灰度分布 f(x, y)和g(x′, y′)，分别为初始点F(x, y)和移动点G(x′, y′)处的灰度值，在两幅图上的两点间的水平、垂直位移分量u(x, y)和v(x, y)可以由下式[4]给出：

在荷载作用下，变形前后相关点的灰度强度具有如下映射关系：

在基准图像中取一个包含点F(x, y)的具有m×n像素的子集f(m, n)。在变形图像中同样划分m×n的小区域g(m, n)。这两个图像子集的灰度相关性按下式[6]进行计算：

2.2 灰度相关性和损伤的关系

由式（4）可知，相关系数 C是变形前后图像小区域的灰度分布函数 f(m, n)和g(m, n)的函数。在荷载作用下，由于材料的不均一性、原始微裂纹的扩展，局部应力集中产生的颗粒剥离等导致试件拍摄面的图像灰度在时间和空间上发生很大的变化，导致g改变增大，C值减小。文献[1]通过相关性与变形的关系模拟试验证明，大变形与低相关度是对应的，因此，灰度相关性在一定程度上反映了试件表面的变形程度。0

图1 散斑图及其相关系数等值线图Fig.1 Speckle images and contours of correlation coefficient

按相关性公式（4），在基准图像（见图 1(a)）中间区域选择一个20×20像素点的方格区域，在图1(a)、图 1(b)中进行匹配。在对两幅散斑图进行灰度相关性计算后，由图 1(c)可知，在图像中心区域有很好的单峰性，即在该区域有最佳匹配点，其灰度相关系数峰值为 1，说明采用的相关系数计算公式有很高的精度。图1(b)是在图1(a)基础上添加随机噪声获得的。采用该方法来模拟荷载作用下材料的劣化引起的灰度变化，由图1(d)可知，峰值明显降低，但单峰性依然明显，其峰值为0.7。可见式（4）对荷载作用下由于材料的劣化所引起灰度变化在空间上仍然有很好的匹配定位性。从而使得基于灰度相关系数峰值改变量和描述岩石空间损伤的损伤状态之间建立起对应关系成为可能。

3 试验结果及分析

3.1 试件、试验装置

为了获得材料力学性质相对均一的试样，采用高强度石膏在钢模中浇注试件（保证试件表面较高的平整度），试件尺寸为5 cm×5 cm×10 cm的立方柱。为了使整个拍摄面的灰度反差较大有利于空间匹配，在试件表面用黑色玻璃微珠漆进行喷涂，形成随机人工散斑场。采用RMT-150型岩石力学试验系统进行单轴压缩加载，加载速率为0.005 mm/s。在整个试件加载过程中用 CCD相机实时记录试件表面散斑场的变化，采集分析软件为DSCM观测系统[1－3]，该系统的核心是图像灰度相关性分析，对整个图像先进行像素搜索，后进行亚像素搜索，获得最佳匹配点，并获得拍摄面的位移场和应变场。图像采集系统的采集速率可控，最高为17帧/s，整个试验过程共记录约450幅图像。

3.2 相关性分布演化与荷载对应关系

单轴压缩应力-应变曲线如图2所示，整个曲线阶段性明显，将曲线分为以下5个阶段：

图2 加载过程中应力-应变曲线及St演化曲线Fig.2 Evolution curve of Stand stress-strain curve during loading process

（1）OA：弹性阶段，试件变形符合弹性变形，此段近似直线，符合线弹性本构关系；

（2）AB：硬化阶段，A点为屈服极限，从该点开始试件产生塑性变形；

（3）BC：快速软化阶段，应变迅速增加，而应力不断减小；

（4）CD：近似理想塑性阶段，其特点是应力无明显变化，应变不断增加，该阶段应力-应变呈理想塑性本构关系；

（5）DE：破坏阶段，应力迅速跌落，试件压碎破坏，仅保持一定的残余强度。

在采集图像中截取感兴趣的代表性区域（ROI）进行计算，截取过程中尽量靠近试件的中间部位，又必须包含尽可能多的物面信息，减小加载过程中端部效应的影响。

采用文献[1]的方法，通过不同时刻图像与基准图像对比计算，以获得的灰度相关性系数的标?准差St为统计指标，其定义为

式中：Ckl为t时刻变形图像中第k个小区域的相关系数值为t时刻变形图像中m×n个小区域的相关系数值的均值。

从图2中可知，在加载初期由于试件原始裂隙的压实导致试件表面灰度调整，使St值有轻微扰动。在 OA、AB段标准差值变化不是很明显；BC段中期由于微裂纹的出现，St值有明显的陡增；CD段由于微裂纹的大量产生和扩展，St值迅速增加；DE段，试件变形急剧增大，表面宏观裂纹产生和贯通，导致C值大幅度降低，St值以更大的速率增加。可见该统计指标能有效反映试件损伤变形的各个阶段的特征。

3.3 损伤阈值的选取

假定OA阶段为完全弹性阶段，以AB段起点时刻A点（弹性应力上限）的试件表面灰度图相关系数图为阈值选取的基准图，如图3(a)所示（以O点时刻的灰度图为基准图像计算获得）。为了获得该时刻灰度相关系数的概率分布密度，采用核密度估计法。

一元连续的总体样本在任意点x处的总体密度函数f(x)的核函数密度估计定义为

式中：K()称为核函数；h称为窗宽。

具体方法为：将通过DSCM系统计算获得的灰度相关系数图由二维矩阵转换为一维矩阵。将该一维矩阵作为总体样本，采用 Gaussian核函数，h=0.0011，求得核密度估计曲线如图3(b)所示。核密度估计曲线与灰度相关系数的频率直方图符合的很好，为了进行对比，另外采用正态分布估计进行拟合，发现灰度相关系数的频率直方图与均值为0.9793，标准差为0.0101的正态分布曲线也非常的接近。

图3 相关系数等值线图及频率直方图、密度分布图Fig.3 Contour, frequency histogram, density distribution of correlation coefficient

为了建立损伤和灰度相关系数改变量之间的关系，依据灰度相关系数的密度分布图（见图3(b)），并认为占总数0.05的最小值部分为初始损伤和采集过程中噪声污染之和，该部分对应的最大 C值为0.9605。以0.9605为二值化转换阈值对整个灰度相关系数图进行图像处理，C＜0.9605则认为该点发生大的变形，处于损伤或破坏状态，在图像矩阵中将其赋值为1；C＞0.9605则认为未损伤，将其赋值为 0。建立以时间为序列的损伤演化二值图（二值图中1代表白色像素，0代表黑色像素）。从AB、BC、CD、DE段中间时刻各选取一幅具有代表性的图，如图4所示。

图4 损伤演化二值图Fig.4 Binary images of damage evolution

这里将发生损伤或破坏的微元数Nd（白色像素总数）与微元总数N（像素总数）的比值定义为损伤变量 Dd，其范围为 0～1。则试件表面的损伤演化如图5所示，在加载初期，处于弹性阶段，试件无明显损伤，随着应变的增加，损伤也随之明显增加，在破坏时有明显的急速上升。

图5 加载过程中应力-应变曲线及Dd演化曲线Fig.5 Evolution curve of Ddand stress-strain curve during loading process

4 损伤分形特征

分形几何是研究自然界不规则现象及其内在规律的学科。岩石材料损伤和破坏过程是其内部缺陷不断萌生发育、扩展、聚集和贯通的最终结果，这个从细观损伤发展到宏观破碎的过程具有分形性质[8]。但整个发展过程并不严格满足自相似性的分形，往往只具有统计自相似性。对于上节损伤演化二值图中损伤微元的分布，本文从统计自相似性出发，用分维数定量描述试件表面的损伤程度。

4.1 分形测定方法及特征

分维数的测定采用基于图像处理的盒维数法[9]，具体方法如下：将相关系数分布图进行阈值处理以后，获得二值图（本文中图像像素为177×399），用边长为l的方格（在图像中为l个像素点），计算图像被边长为l的方格覆盖占有的方格数N(l)，不断改变方格边长l的大小，并计算其相应的N(l)，可得到一组N(l)和l数据。根据盒维数的基本定义，损伤分维数D可由下式计算：

lgN(l)-lg(1/l)的曲线如图 6所示，两者呈线性关系，相关系数 R=0.9865，表明损伤微元分布具有自相似性，可以用分形进行表述。

为获得图2中各个阶段的损伤与应变的对应关系，从塑性变形开始往后各阶段的中间部分（避免转折点处的突变）以5 s为时间间隔，获得各个阶段损伤演化二值图。AB段选择7幅图，BC阶段9幅图，CD阶段9幅图，DE阶段6幅图。分别计算各阶段在不同应变条件下的分形维数。

在对不同阶段下相应损伤演化二值图进行分形分析后，提取各个所选二值图的分形维数。其维数的变化如图7所示，随着应变增加，各阶段分形维数的范围为 1.011～1.055、1.092～1.208、1.210～1.362、1.365～1.469，最小值为 1.011，最大值为1.469。可见随着应变的增大，材料表面细观结构的损伤程度增加，其对应的损伤演化二值图上的损伤微元越来越多，分形维数也随之增大；损伤微元由整体随机分布，向局部集中发展。宏观上表现为微裂纹的萌生、扩展，宏观裂纹的出现。盒维数的大小反应了损伤的集中程度。不同阶段分形维数和应变的一元线性回归拟合系数分别为 0.0305、0.0567、0.0738、0.0930，由此可知，在不同阶段损伤的发展是一个逐渐加快的过程，也从另外一个方面反应了试件加载作用下，力学性质加速劣化直到破坏的过程。

图7 损伤分形维数与应变关系图Fig.7 Relationships between fractal dimension and strain

4.2 多重分形测定方法

简单的分形维数对试件表面损伤的演化只能作整体性、平均性的描述与表征。多重分形是定义在分形结构上的由无穷多个标度指数所组成的集合，是通过一个谱函数来描述分形结构上不同的局域条件或分形结构在演化过程中不同层次所导致的特殊结构行为与特征，是从系统的局部出发来研究其整体的特征，并借助统计物理学的方法来讨论特征参量概率测度的分布规律[10]。

本文采用Ghhabra和Jensen[11]于1989年提出的计算多重分形谱的方法，首先用边长为l的方格对图像进行分网，损伤分布的总像素点为N，第i个方格所包括的损伤区域（白色像素点）为Ni，从而得到每个尺度为l的方格中分形图像的概率测度：

Pi(l)的归一化概率测度为

式中：q为权重因子，不同的q表示不同大小的概率测度 Pi(l)在配分函数pi(l)q中所具有的比重。

a(q)为描述概率测度局部奇异性强度的参数，表示多重分形集整体奇异性的平均值，其表达式为

f(a)为用奇异性标度指数标识的分形子集的Hausdorff维数，其表达式为

奇异谱 f[ a (q)]为 q∈(-∞,+∞)的单参数曲线,由于 f[ a (q)]与a(q)有惟一对应性，不同的a(q)及其所对应的 f[ a (q)]便构成多重分形的维数谱。

多重分形谱绘制的具体步骤为：首先程序读入损伤二值图，统计损伤微元总数N（图像中白色像素点总数），再对二值图像进行网格划分（2n，n=1～8），程序对图像中划分的网格进行逐个识别统计，获得每一方格损伤微元数目Ni，建立概率测度Pi=Ni/N，读入q值，通过式（9）建立 ui(q, l)，计算获得 ∑ui(q, l) lg[ui(q, l)]和 ∑ui(q, l) lg[pi(l)]的数值，并分别与lgl建立数据对，对上述两组数据对采用最小二乘法进行一元线性回归拟合，拟合出a(q)和 f[ a (q)]值，即可绘制多重分形谱图。基于以上思路，用Matlab编制了基于数字图像处理的多重分形谱计算程序。

4.3 损伤演化多重分形特征

图8 ε=0.35%时多重分形谱Fig.8 Mutifractal spectrum with ε=0.35%

表1中仅列出应变为0.35%、0.70%、0.88%时的多重分形谱值。实际计算中 q取值范围为-10～10，步长为 0.1。多重分形谱中大的a(q)反映的是小概率测度区域的性质，小的a(q)反映的是大概率测度区域的性质。图9为图4中4个不同应变条件下损伤演化二值图的多重分形谱图，其曲线由开始的对称单峰形式，逐渐发展为左侧相对集中，右侧稀疏的钩状。 a(q)-f[ a(q)]曲线上的数据点在左侧逐渐相对集中，表明大概率分布逐渐占主导地位，损伤演化从整体随机均匀分布，到局部集中的过程。

图9 多重分形谱Fig.9 Multifractal spectrum

多重分形谱比简单分形包含了更丰富的结构信息，可以从分形谱中提取低维特征量来表征损伤的内在分形机制。当q=0、1、2时，D0为容量维，D1为信息维，D2为关联维，其数值见表 2，三者的关系为D0>D1>D2。

表2 不同应变下D0、D1和D2的值Table 2 Values of D0, D1and D2at different strains

多重分形谱的宽度Δa=a (q)max-a(q)min的大小则反映了整个分形结构上概率测度分布的不均匀性的程度。Δa越大，表示损伤分布越不均匀，Δa越小，则表示损伤分布相对集中。ε=0.35%时，谱宽为0.381， a(q)-f(q)曲线呈对称分布；ε=1.00%时，谱宽为 0.345，a(q)-f[ a (q)]曲线相对集中，其多重分形性相对减弱，损伤演化的奇异性也减小（见表3）。

表3 不同应变下Δα的值Table 3 Values of Δα at different strains

5 结 论

（1）荷载作用下，试件表面的数字图像灰度相关性和试件表面的损伤演化在时间和空间上存在对应关系，并可以用灰度相关系数C值的变化来对损伤状态进行表述。

（2）通过简单分形维数对由灰度相关系数建立的损伤演化二值图进行刻画，发现在各个不同阶段，损伤分维数与应变呈线性关系。不同阶段的增长速率不同，越接近破坏，增长速率越大。反映了材料在荷载作用下的加速劣化、破坏过程。

（3）通过多重分形谱反映了试件表面损伤演化在分形结构上不均匀分布的性质。其曲线形状从开始的对称光滑单峰形式，逐渐向左侧集中发展，其谱宽逐渐减小，大概率分布逐渐占主导地位，反映了损伤演化从整体随机分布，到局部集中的过程。

本文采用概率统计的方法，建立图像灰度相关系数的改变和试件表面损伤状态之间的关系，并对材料表面的细观结构破坏进行了表述。文中仅仅对损伤状态进行表述，更深层次的损伤机制及程度则需要更进一步的研究。

[1] 马少鹏, 刘善军, 赵永红. 数字图像灰度相关性用以描述岩石试件损伤演化的研究[J]. 岩石力学与工程学报,2006, 25(3): 590－595.MA Shao-peng, LIU Shan-jun, ZHAO Yong-hong. Gray correlation of digital images from loaded rock specimen surface to evaluate its damage evolution[J]. Chinese Journal of Rock Mechanics and Engineering, 2006,25(3): 590－595.

[2] 马少鹏, 周辉. 岩石破坏过程中试件表面应变场演化特征研究[J]. 岩石力学与工程学报, 2008, 27(8): 1667－1673.MA Shao-peng, ZHOU Hui. Surface strain field evolution of rock specimen during failure process[J]. Chinese Journal of Rock Mechanics and Engineering, 2008,27(8): 1667－1673.

[3] MA S P, WANG L G, JIN G C. Damage evolution inspection of rock using digital speckle correlation method[J]. Key Engineering Materials, 2006, 326－328:1117－1120.

[4] 赵永红, 梁海华, 熊春阳, 等. 用数字图像相关技术进行岩石损伤的变形分析[J]. 岩石力学与工程学报, 2002,21(1): 73－76.ZHAO Yong-hong, LIANG Hai-hua, XIONG Chun-yang,et al. Deformation measurement of rock damage by digital image correlation method[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(1):73－76.

[5] 朱珍德, 张勇, 李术才, 等. 用数字图像相关技术进行红砂岩细观裂纹损伤特性研究[J]. 岩石力学与工程学报, 2005, 24(7): 1123－1128.ZHU Zhen-de, ZHANG Yong, LI Shu-cai, et al. Studies of microcosmic crack damage properties of a red sandstone with digital image technique[J]. Chinese Journal of Rock Mechanics and Engineering, 2005,24(7): 1123－1128.

[6] CHU T C, RANSON W F, SUTTON M A. Applications of digital-image-correlation techniques to experimental mechanics[J]. Experimental Mechanics, 1985, 25(3):232－244.

[7] MA S P, JIN G C. Digital speckle correlation method improved by genetic algorithm[J]. Acta Mechanica Solida Sinica, 2003, 16(4): 366－373.

[8] 彭瑞东, 鞠杨, 谢和平. 灰岩拉伸过程中细观结构演化的分形特征[J]. 岩土力学, 2007, 28(12): 2579－2582.PENG Rui-dong, JU Yang, XIE He-ping. Fractal characterization of meso-structural evolution during tension of limestone[J]. Rock and Soil Mechanics, 2007,28(12): 2579－2582.

[9] 曹茂森, 任青文, 翟爱良, 等. 混凝土结构损伤的分形特征实验分析[J]. 岩土力学. 2005, 26(10): 1570－1574.CAO Mao-sen, REN Qing-wen, ZHAI Ai-liang, et al.Experimental study of fractal characterization in damages of concrete structures[J]. Rock and Soil Mechanics,2005, 26(10): 1570－1574.

[10] 李平, 胡可乐, 汪秉宏. 多重分形谱在材料分析中的应用研究[J]. 南京航空航天大学学报, 2004, 36(1): 77－81.LI Ping, HU Ke-le, WANG Bing-hong. Design and application about computing program of material multifractal spectrum[J]. Journal of Nanjing University of Aeronautics & Astronautics, 2004, 36(1): 77－81.

[11] CHBABRA A, JENSEN R V. Direct determination on the f(α)singularity spectrum[J]. Physical Review Letters,1989, 62(12): 1327－1330.

[12] 王金安, 谢和平. 岩石断裂面的各向异性分形和多重分形研究[J]. 岩土工程学报, 1998, 20(6): 16－21.WANG Jin-an, XIE He-ping. On anistropic fractal and multifractal properties of rock fracture surface[J].Chinese Journal of Geotechnical Engineering, 1998,20(6): 16－21.

[13] 曹茂森, 任青文. 钢筋混凝土结构损伤检测的分行特征因子[J]. 土木工程学报, 2005, 38(12): 59－64.CAO Mao-sen, REN Qing-wen. Damage detection of reinforced concrete structures based on fractal characteristic factor[J]. China Civil Engineering Journal, 2005, 38(12): 59－64.