数字图像相关方法(Digital Image Correlation, DIC)由于具有非接触、系统组成简单、全场测量等优点在众多领域中的应用越来越广泛^[1-5]。DIC是基于数字图像处理的一种全场变形测量技术，其基本原理是选取变形前感兴趣的像素点的灰度子区，通过相关计算，跟踪变形后感兴趣点的位移变化。全场计算时，需要对每个感兴趣区域的种子点进行整像素或者亚像素匹配以达到全场测量的目的^[6]。为了确保测量点变形前后的唯一性以及测量结果的准确性，一般测量前需要对被测物或被测区域随机喷涂两种或多种颜色差异较大的涂料形成散斑^[7]。如文献[1]中提到在高温测量下不同温度段最优成像对应的散斑图样。文献[4]人车碰撞实验中会对撞击物与车前盖喷涂黑白散斑，利用DIC方法计算碰撞后的头部伤害指数(HIC)。

由于DIC的计算效率受到测量点个数的影响，测量点个数越多，精度越高，但所花费的时间也越长。如果不事先选择散斑区域约束计算点的范围，会产生很多无用结果，浪费了大量的计算时间，而且在实际拍摄的图像中散斑区域一般只占图像的一部分。长期以来，国内外学者在实验过程中往往会人工选择散斑区域并设定一个种子点进行计算^[8]，目的是限制相关搜索范围以减少计算点的数量，但面对大量不同的测量图片以及复杂的轮廓，手动选取散斑区域无疑增加了时间成本和人工成本。近几年散斑质量评价方法得到广泛研究^[9-10]，但目前提出的散斑评价方法前提都是整幅图像全是散斑，而在实际应用中，由于背景的影响，直接对图片进行散斑质量评价显然不再适用。在其他研究领域仍以抑制散斑为主要研究课题^[11-13]。在DIC测量中实现图像散斑区域的自动提取显得尤为重要。这有助于提升DIC的应用效果及实现DIC自动化。

为了实现DIC散斑区域的自动提取，本文通过分析散斑区域特征，对比常用的边缘检测算法，提出了一种二阶梯度熵函数用于实现对散斑区域的判别，并根据设定的连通区域获取真正的散斑区域。最后通过实验拍摄的散斑图像来验证该方法的有效性。

DIC散斑区域的自动提取方法相当于一种特征提取技术，散斑区域即所需要提取的区域，散斑区域具有以下几个明显特征：(1)散斑区域内的像素点周围子区灰度梯度大；(2)散斑区域子区内最大灰度值与最小灰度值差值较大；在图像上该特征的表现方式是散斑对比度明显，并具有随机的灰度分布特征；(3)散斑区域内的点构成的区域具有连通性，即符合(1)(2)特征的像素点并不是孤立的点，在周围必定有相似的像素点分布。利用上述3个特征，可以将非散斑区域与散斑区域进行分割。

在图像特征提取中，常用的边缘检测算子一般分为3类：一阶梯度算子、二阶梯度算子与多级检测算子^[14]。代表类型有Sobel算子、Laplacian算子与Canny算子。通过分析散斑区域的特征不难发现，散斑区域内的灰度变化曲线非常陡峭，对所有边缘检测算子都很敏感，所以不能通过边缘检测算子对散斑区域直接进行分割。图 1(a)为一张复杂背景下的散斑原始图像。灰度均匀变化的背景，经过一阶梯度算子计算的结果如图 1(b)所示。可见，其不能很好地将散斑区域与背景区域分离，结果不能让人满意；多级检测的Canny算子计算的结果如图 1(c)所示，其通过调节高阈值与低阈值可以很好地消除图片背景的伪边缘，但是由于对边缘过于敏感，且最终输出图片的形式为二值化，不仅没有消除背景带来的影响，还将散斑的特征减弱了; 而经过二阶梯度算子计算的结果如图 1(d)所示，消除了部分背景的影响，将图像灰度值转为梯度值，适合用作散斑提取前的预处理。

图  1  几种边缘检测算子处理结果对比

Figure 1.  Comparison of processing results from several edge detection operators

本文用拉普拉斯算子处理采集到的散斑图像，目的是保留散斑区域的特征，剔除背景的影响。拉普拉斯算子是二阶梯度算子，对于连续可微的图像灰度函数或曲面f(x, y)而言，其Laplacian图像g(x, y)定义为：

在数字图像中，图像像素点以离散形式表达为：

式中，f(x, y)代表图片离散坐标(x, y)处的灰度值，g(x, y)代表利用拉普拉斯算子预处理过后对该坐标重新赋予的灰度值。从上述的表达式可以看出，如果图像中像素点的灰度值与相邻像素点的灰度均值间的差值越大，则该点灰度值会越大，即通过拉普拉斯算子处理后，会基本保留图像边缘与灰度值跳跃明显的区域。如图 2所示，在散斑区域和背景区域中各取一部分，统计拉式变换前后的灰度直方图。在灰度直方图中可以直观看出，利用该算子进行图像预处理，处理前后散斑区域的灰度分布几乎保持不变，而背景区域灰度分布从0~175缩小至0~40。

图  2  拉式变换后散斑区域与背景区域灰度直方图的对比

Figure 2.  Comparison of gray histograms between speckle area and background area after laplacian transformation

由于散斑区域较其他区域的特征对比度更高，灰度分布更广，所包含的信息更多。以点(x, y)为中心的一定邻域内，散斑判定值R_s(x, y)与散斑信息量I_s(x, y)有如式(3)所示关系，可以认为其R_s(x, y)值越大，则该点越符合人们所认知的散斑内的点的要求。

为了统计散斑所含信息量的大小，本文引入熵的概念。熵一般是指统计混乱的程度，在不同的学科中也引申出了不同的定义，熵的概念最先由Rudolf Clausius提出，并应用在热力学中。后来Shannon第一次将熵的概念引入到信息论中。在信息论中，设x_n为信源中的随机变量，其出现的概率设为p(x_n)，则定义信源的信息熵为^[15]：

将熵的概念推广到图像中，即可统计图像信息量的大小, 一幅8位BMP图像，设灰度值i出现的概率为ρ_i, 则图像信息熵可定义为：

式(5)反映了图像平均信息量的大小，对于纯黑或纯白的图像，显然图像信息熵值为零，则可以认为图像所带的信息量很小。但是当灰度直方图平稳，即图片中灰度值出现的概率都相等，那么图像信息熵的值相对比较大，此时可以认为图像所含信息量较大。

一个适合用于计算图像信息量的函数应该具有较好的单峰性，较高的灵敏度。基于图像灰度分布的图像信息熵函数虽然能很好地统计出图像信息量的大小，但受实验条件的影响，背景区域的信息量是未知的，因此无法将散斑区域与背景进行有效分离。考虑到图像信息熵这一缺陷，本文从研究散斑特征出发，提出一个适用于提取散斑特征的函数：二阶梯度熵函数Q_M。消除背景影响，利用Q_M对图像中的散斑区域进行判定。

其中，M代表检验子区的边长大小，ρ_{∇ (i, j)}是当前像素点(i, j)的二阶灰度梯度值在M×M大小子区中出现的概率。以M×M的子区遍历图片，可以获得图像的二阶梯度熵值L_Q序列。

从L_Q序列中可以找到最大值Q_max, 设定散斑所在的熵值范围处于[Q_max-k, Q_max], 当图像中以点(x, y)为中心的邻域对应的二阶梯度熵值Q_{(x, y)}处于[Q_max-k, Q_max]时，将该点置为1，否则将该点置为零，这样可获得分割后的二值化图像L_b。以图 1(a)为例，计算其二阶梯度熵，取k为1.25，得到的二值化图像如图 3所示。

图  3  二阶梯度熵值分割结果

Figure 3.  Second order gradient entropy segmentation results

在二值化图像L_b中，根据数值1所在的区域可获得连通区域的面积数组[S₁, S₂, …, S_n]，将连通区域的面积从大到小排列，若设定散斑区域个数为N，最终散斑所在的区域为[S₁, S₂, …, S_N]。以图 3为例，设定散斑区域数量为1，散斑区域提取结果如图 4所示。

图  4  散斑区域提取结果

Figure 4.  Extraction results of speckle area

综上所述，基于二阶梯度熵函数的散斑区域提取算法实现流程如图 5所示。

图  5  散斑区域提取算法的实现流程图

Figure 5.  Implementation flow chart for speckle area extraction algorithm

对于式(6)提出的二阶梯度熵函数，显然，子区范围越大时，Q值的判定结果越准确，但子区范围过大会导致计算量增加，因此需要确定一个合适的子区范围，既能保证计算速度也能确保结果的准确性。本文利用计算机生成的模拟散斑来对子区大小进行对比分析，模拟散斑生成函数如下^[16]：

其中，s是散斑的个数，a表示为模拟散斑的像素大小，C表示模拟的背景光强，一般设为255，(x_n, y_n)表示散斑颗粒的位置。

前期有学者认为散斑颗粒尺寸为5个像素是最优的散斑尺寸^[17]，因此选取a为5，散斑个数为2 000个，模拟散斑图片的大小为512×512，利用图像信息熵函数(式5)与二阶梯度熵函�