Initial estimation of digital image correlated deformation based on genetic algorithms

数字图像相关法(digital image correlation, DIC)是由YAMAGUCHI等人^[1-2]提出的一种非接触图像测量方法, 该测量方法光路简单、抗干扰、适合全场测量、使用范围广^[3]。随着数字图像处理技术的发展和图像采集分辨率的提高，DIC方法在变形测量中的应用愈加广泛^[4-8]。DIC方法按其测量精度可分为整像素DIC测量方法和亚像素DIC测量方法，其中整像素DIC方法有粗细搜索法、十字搜索法、人工鱼群算法、粒子群算法等相关搜索方法，亚像素DIC方法有曲面拟合、灰度梯度等。牛顿-拉夫森迭代法(Newton-Raphson, N-R)是一种常用的亚像素搜索方法，其收敛范围一般在几个像素^[4]，但其初值估计的准确度对其收敛速度和计算精度影响较大。测量小变形物体表面时，目标点变形前后只有几个像素的位移量，即使有测量误差，N-R法通过一定的迭代次数仍能找到最优解^[9]，但测量大变形物体表面时，通过数字图像相关法匹配的目标点变形前后像素位移量会骤增，测量误差会被放大，使得N-R迭代次数增加且增大计算量，甚至可能陷入局部峰值^[10-13]。

遗传算法(genetic algorithm, GA)是一种全局搜索方法，该方法搜索速度快，测量结果准确，适用于大变形对象的整像素位移测量^[14-17]。该方法的计算时间受搜索区域增大的影响小，是一种有效的大范围搜索方法，但在末端定位时存在局部震荡现象。本文中提出一种基于遗传算法的数字图像相关变形初值估计法。首先，选取待测数据点邻域内若干估测点，通过一种遗传算法和粗细搜索法的混合方法进行整像素搜索，匹配出变形前后估测点的坐标值变化。等概率选取3组或以上不共线的估测点对进行仿射变换并拟合出具有6个参量的仿射变换模型，依据仿射变换结果得到待测数据点的变形初值估计。最后，将变形初值代入N-R迭代算法中，结合双三次样条插值方法获取最终更精确的亚像素位移值。仿真和橡胶圆柱压缩实验结果表明，本文中的方法有效减少匹配时间，在大变形情况下测量精度仍然保持稳定，与传统方法相比有明显优势。

物体表面常存在自然地斑点纹理，以此作为特征可以估算两幅图片的相似度。DIC方法是以人工制作的细小颗粒状散斑图像为特征，匹配变形前后两幅图像上对应子区相似程度的方法。其基本原理如图 1所示。

Figure 1.  Subset region of image before and after deformation

采集变形前后两幅图像进行对比，变形前图像记为参考图像，变形后图像记为变形图像。在参考图像中选取的一定大小的矩形参考子区，其中心为P(x₀, y₀)，依据变形前后像素点灰度值不变的原则，变形图像必然存在一与参考子区相似度最高的子区图像，记为变形子区，其中心点坐标为P′(x₀′, y₀′)，参考子区和变形子区的相似度通常采用相关函数来度量，依据TONG等人^[18-19]对相关系数的研究可知，零均值归一化互相关函数对图像灰度值的线性变化不敏感，适合在光照环境有微弱变化的环境下使用，一定程度降低外界光照干扰对匹配结果的影响，本文中采用该方法，表达式为:

式中，R为图像子区半径，f(x, y)为参考图像中点(x, y)处的灰度值；g(x′, y′)为变形图像中点(x′, y′)处的灰度值；〈f〉和〈g〉为对应子区的灰度平均值。变形子区和参考子区的匹配度与相关系数C值正相关，C取值在[0, 1]之间，C=1代表二者完全匹配。

为了解决传统DIC方法迭代运算受初值影响较大且容易陷入局部最优的问题，提出基于遗传算法的数字图像相关法。首先，根据待测点坐标选取邻域内3组不共线的估值点，通过遗传算法在大区域匹配各个估值点在变形后的空间域信息，得到待测数据点周围变形前后的若干估值点对，以此为依据根据仿射变换原理求取待测数据点位移变化估计初值和待测数据点变形后坐标值。最后，将数据点位移变化估计初值代入N-R迭代算法求取最终亚像素位移值。算法流程如图 2所示。

Figure 2.  Flow chart of DIC based on GA

遗传算法是模拟自然生物进化过程的算法，本文中选取100组初始解S =[u, v]作为初始种群，其中u, v分别代表x, y方向上的位移值，其取值为[-25, 25]间随机数。通过一定方式对其编码，对个体进行选择、交叉、变异等操作，依据适应度函数计算数值筛选出性能优良的个体并保留其基因至下一代，代代之间不断优化，直至满足最优条件。

轮盘赌法是遗传算法中最常用的选择方法之一，但该方法选择过程中有一定概率跳过最优解。本文中将轮盘赌法和最佳保留选择法结合使用，经过轮盘赌法选择后如果当代中最优个体没有满足要求，则进行最佳保留策略，选出n个个体中最差的5个，用当代最优的5个个体将其替换。

遗传算法中通过两个体之间染色体交叉产生新个体。交叉算子有简单交叉、启发式交叉、算术交叉等，此处采用算数交叉，交叉产生的两个新的个体表示为:

式中, 问题的规模i=1, 2，…, n；S_i为种群中编号为i的个体，β为[0, 1]之间随机数。此处交叉发生概率设为P₁=0.8。

变异操作发生概率较小，此处设其概率为P₂=0.01。变异算子常有均匀变异、非均匀变异、多维正态变异、边界变异等，此处采用非均匀变异操作。取第l次迭代时种群中的个体S^l进行变异, 变异生成的新个体表达式为:

式中, L为最大迭代次数；Ν_M×M为对角矩阵且其对角线元素为[0, 1]间随机数；M为个体染色体向量的维数；D_u为染色体分量最大值，D_l为染色体分量最小值；b代表对迭代数依赖程度。

经过上述选择、交叉、变异产生的下一代群体，若最优个体的适应度值不满足要求，则重新进行个体适应度评估并进入下一轮循环，直到最优个体满足要求。遗传过程的停止条件为两代种群中的最优个体和当代种群中的最优个体相同或者进行到最大迭代次数。

传统遗传算法进行整像素搜索时，实现了全局寻优搜索，但其收敛位置受到交叉和变异操作的影响，有一定概率偏离真实位置1~2个像素，造成N-R迭代收敛效率的降低。本文中将遗传算法和粗细搜索法结合，如图 3所示。首先通过遗传算法在大范围内进行全局搜索，借助其快速收敛的优势，在短时间内找到真实位置的初步估算位置A。遗传算法存在末端收敛局部震荡问题，直接使用遗传算法的搜索结果会存在一定像素偏差，故引入粗细搜索思想，在最终末端定位时以A为中心，选取周围5×5像素的子区逐点搜索使得相关函数达到最大值的点作为最终结果，有效消除末端收敛的局部震荡现象，提升整像素搜索的稳定性。

Figure 3.  Schematic of hybrid algorithm

针对遗传算法有一定概率陷入局部最优解的现象，本文中对比前后两个待测点检测位移值之差，若超出设定阈值，则重新计算，从而有效避免陷入局部最优。

在参考图像中选取感兴趣区域，以51×51^[11]像素大小对该区域进行网格划分，得到若干窗口，每个窗口的中心即为一个待测点。以其中一个待测点为例，取其邻域内的若干待测点作为估值点，通过遗传算法计算得到第i个估值点变形前后的位置坐标分别为h_i=(x_i, y_i)^T和h_i′=(x_i′, y_i′)^T。随机选取3组估值点对代入仿射变换，其表达式为:

表达式中通过3组估值点对的坐标信息既可求解包含6个参量a₁, a₂, a₃, a₄, a₅, a₆的仿射变换。通过仿射变换即可求得变形后待测点的精确坐标位置，从而得到相对变形前待测点的像素位移估值，并以此作为N-R迭代的初始值。该方法相较于传统方法，计算出的位移初值更接近待测点的实际位移，有效降低了后续N-R迭代运算次数，缩短收敛时间。

为提高测量精度，常需要采用曲面拟合法、灰度梯度法、N-R迭代等方法计算亚像素位移值。其中N-R迭代法是一种常用的获取亚像素位移值的方法，其具有精度高、计算结果稳定可靠等优点，因此，本文中选用该算法计算亚像素位移值。进行N-R迭代之前需要通过灰度差值方法获取亚像素位置的灰度值和灰度梯度等信息，本文中采用精度较高的双三次样条插值法，其表达式为:

式中, f(x, y)为插值区域中位于(x, y)处的灰度值，a_ij为插值系数，对f(x, y)求偏导即得到各方向灰度梯度值。通过N-R迭代运算后求取的亚像素位移值即为最终结果。

为了验证本文中提出算法在计算精度和计算性能上面的提升，在计算机上使用本文中算法和其它传统算法对模拟散斑图像进行特征点匹配，计算亚像素位移值。参考ZHOU^[20]提出的模拟散斑图模拟物体表面的变形过程，变形前后的模拟散斑图灰度表示为:

式中，s为散斑颗粒数量，r为散斑颗粒半径，(x₀, y₀)为散斑图中心位置，I₀为光强分布; u, u_x, u_y, v, v_x, v_y为散斑的变形参量，决定图像的位移量和应变量。为验证不同的纹理对测量的影响，将散斑颗粒大小r、颗粒数目s分别设置为3pixels、2500和2pixels、4500两组，如图 4所示。每组v_x应变值按照0.005的步距由0.001过渡到0.05。

Figure 4.  The simulated speckle image generated by MATLAB

分别利用传统N-R迭代法、曲面拟合法与本文中算法对目标图像中随机数据点进行匹配。比较不同散斑颗粒尺寸、不同散斑颗粒数目和不同应变值对测量平均误差和标准差的影响，结果如图 5、图 6所示。最后对迭代次数和匹配时间进行对比，如表 1所示。

Figure 5.  Average measurement error

Figure 6.  Standard deviation of measurement error

模拟散斑实验结果表明，随着应变的增加，3种算法的误差平均值均呈现上升趋势，且曲面拟合方法的误差较大，本文中算法生成较为准确的初值估计，一定程度上降低了平均误差值，但在应变值小于0.01时误差值稍大于传统N-R迭代法，这是仿射变换过程引入了计算误差的缘故。本文中算法和传统N-R迭代法的标准差较为稳定，曲面拟合法标准差数值较大且有震荡现象。应变进一步增大后，本文中算法的标准差值相对传统N-R迭代法的标准差值降低，说明在大变形情况下本文中算法能够保持稳定性。同时，对比结果可看出散斑颗粒大小和颗粒数量的变化对变形测量的影响不大。从表 1可看出，本文中算法一定程度上也减少了N-R迭代运算次数，对比两种方法的平均计算时间可发现，本文中算法的匹配时间相对传统方法平均降低37.54%, 在计算效率上有一定的提高。

为验证本文中算法在实际应用中的可靠性，选取直径60mm、高度150mm的橡胶圆柱棒材进行压缩变形测量实验。在试件表面制作散斑图案，将试件放置在位移分辨率为0.001mm的精密伺服压力机上，以0.1mm为步距对试件压缩，试件压缩量分别为0.5mm, 0.6mm, …, 3.5mm。采用Baumer的LXG120M型号CCD相机采集图像，帧频为10frame/s。具体实验装置如图 7a所示。

Figure 7.  Experimental device and sample images

为了获得材料在压缩后的变形场信息，在一系列压缩图像中选取两幅图像作为图像处理的样本，如图 7b、图 7c所示。采用本文中算法和传统N-R迭代法对橡胶压缩变形表面进行区域位移场分布计算，计算过程中数据点之间的间距设置为51个像素，计算结果如图 8所示。结果显示位移场中矢量箭头朝向试件的右下方，这是由于所选测量区域位于橡胶中部偏右，轴向的压缩位移和径向的膨胀位移同时作用。将位移场沿x和y方向分解，结果显示y方向为主要变形方向，x方向由左向右位移值逐渐增大，符合橡胶圆柱压缩过程中轴向变形为主，膨胀变形为辅的规律，等值线分布图如图 9所示。通过位移场可以计算出应变分布，本文中算法感兴趣区域内测量应变均值为0.0061，与采用ANSYS 17.0仿真分析中该区域0.0069应变值接近，说明测量算法在实际应用中可靠。

Figure 8.  Displacement field distribution

Figure 9.  Displacement field distribution of x and y direction

提出一种基于遗传算法的数字图像相关法，结合遗传算法与粗细搜索混合算法的全局搜索性能准确匹配出整像素位移, 后续通过仿射变换原理根据待测点周围若干变形前后的估值点对推算出待测点更加准确的位移估值，以此代入N-R迭代算法，加快其收敛速度和精度。模拟实验表明，本文中采取的混合算法使得各算法优势优势互补，对不同的散斑图像均有适应性，且相较于传统DIC方法，在大应变测量中能够有效提高测量精度且保证测量结果稳定性。橡胶圆柱压缩实验结果验证了本文中算法在实际大变形测量中的可行性。后期将使用该方法进行金属材料的大变形测量和振动抗干扰分析。