CN105976356A - 一种基于相关熵准则的鲁棒数字图像相关方法 - Google Patents

Abstract

本发明公开了一种基于相关熵准则的鲁棒数字图像相关法，采用基于相关熵准则的相关函数C^*评价变形前后图像子区的相似程度，相关函数表达为本发明的基于相关熵准则的鲁棒数字图像相关法可以很好地处理非连续变形，并且较传统数字图像相关法而言计算精度更高，抗噪性能更强。

Description

一种基于相关熵准则的鲁棒数字图像相关方法

技术领域

本发明属于非接触测量领域，涉及一种基于相关熵准则的鲁棒数字图像相关方法。

背景技术

数字图像相关技术(Digital Image Correlation,DIC)将材料表面变形测量问题转化为材料变形前后图像特征的区域匹配问题，从而获得变形后的全场位移数据，由于DIC具有非接触、测量精度高、全场变形测量便捷、实验台架搭建容易等特点，在力学相关领域的应用越发广泛。影响数字图像相关技术测量精度的因素很多，如散斑特征、图像子区大小、形函数阶次、收敛条件、相机噪声、光源波动、镜头畸变，或是变形状态等，已有一些学者对这些影响因素进行了研究并给出一些指导应用的结论和原则。

对于实际应用过程中的问题，如试件发生断裂变形分析、试件边界变形分析、含孔试件孔洞边缘大变形分析等，传统DIC是一种基于均方误差准则的“全局”算法，该“全局”的含义是整个参考图像子区中所有的计算点以同样的权重和地位参加DIC的搜索匹配运算，对于这些非连续变形问题，显然参考图像子区中存在不符合连续变形模式约束的点，这些“坏点”的参与无疑会干扰位移向量匹配运动的方向和步长。噪声在实际应用中是不可避免的，而由于均方误差准则“全局”的特性，被噪声污染的点也以同样的地位参与了匹配运算，使得传统DIC方法的抗噪性能不强。

因此，需要一种新的数字图像相关方法以解决上述问题。

发明内容

本发明的目的是针对在现有技术中的数字图像相关方法的不足，提供一种基于相关熵准则的鲁棒数字图像相关方法。

为实现上述发明目的，本发明的基于相关熵准则的鲁棒数字图像相关方法可采用如下技术方案：

一种基于相关熵准则的鲁棒数字图像相关方法，采用基于相关熵准则的相关函数评价变形前后图像子区的相似程度，所述基于相关熵准则的相关函数通过下式表达：

$C^{*}\left(P\right)=\frac{1}{c}\underset{x=-M}{\overset{M}{\Σ}}\underset{y=-M}{\overset{M}{\Σ}}\exp \[-{\left(f\right(x,y)-g(x,y;P\left)\right)}^2/2{\mathrm{\σ}}^2\]$

其中，c为归一化常数，M为参考图像子区的半径，exp()为高斯核函数，σ为高斯核函数的带宽，(x,y)为参考图像子区中像素点的坐标，f(x,y)为参考图像子区中像素点(x,y)的灰度值，g(x,y)为目标图像子区中与像素点(x,y)相对应的像素点的灰度值，P为变形参数向量。

更进一步的，包括以下步骤：

(1)、采集被测件表面变形前和变形后的散斑图像，变形前的散斑图像称为参考图像，变形后的散斑图像称为目标图像；

(2)、在参考图像上划分感兴趣的网格点；

(3)、在参考图像上，建立以步骤(2)的网格点为中心的第一矩形区域，称为参考图像子区，第一矩形区域的长度为L，宽度为W；

(4)、在目标图像上，建立以步骤(2)的网格点为中心的第二矩形区域，第二矩形区域的长度为k*L，宽度为k*W，其中，k>1，称为目标图像子区，目标图像子区包含变形后的参考图像子区；

(5)、获得参考图像子区中整像素位置的灰度值和目标图像子区中整像素位置以及亚像素位置的灰度值；

(6)、定义形函数：设定参考图像子区中心点为(x₀,y₀)，参考图像子区中的任意点为(x,y)，变形后目标图像子区中与点(x,y)对应的点为(x',y')，存在一组映射关系ξ使得下式成立：

ξ(x,y)→(x',y')

f(x,y)＝g(x',y')

其中，f(x,y)表示变形前像素点(x,y)处的图像灰度值，g(x',y')表示变形后对应像素点(x',y')处的图像灰度值；

(7)、将步骤(6)的形函数参数化，并用变形参数向量P表示，定义基于相关熵准则的相关函数：

$C^{*}\left(P\right)=\frac{1}{c}\underset{x=-M}{\overset{M}{\Σ}}\underset{y=-M}{\overset{M}{\Σ}}\exp \[-{\left(f\right(x,y)-g(x,y;P\left)\right)}^2/2{\mathrm{\σ}}^2\]$

其中，c为归一化常数，M为参考图像子区的半径，exp()为高斯核函数，σ为高斯核函数的带宽，(x,y)为参考图像子区中像素点的坐标，f(x,y)为参考图像子区中像素点(x,y)的灰度值，g(x,y)为目标图像子区中与像素点(x,y)相对应的像素点的灰度值，P为变形参数向量；

(8)、利用灰度直方图计算得到高斯核函数的带宽σ；

(9)、将步骤(8)得到的高斯核函数的带宽σ代入所述基于相关熵准则的相关函数，计算得到令基于相关熵准则的相关函数C^*(P)最大时的变形参数向量P。

更进一步的，步骤(1)中被测件表面的散斑为自然纹理特征的散斑或人工制备的散斑。

更进一步的，步骤(5)中通过插值方法获得目标图像子区亚像素位置的灰度值。

更进一步的，所述插值方法为最近邻插值法、双线性插值法或双三次样条插值法。

更进一步的，步骤(6)中的映射关系ξ为一阶形函数或二阶形函数，

其中，一阶形函数通过下式表达：

$\left\{\begin{array}{c}{x}^{\′}=x+P_1+P_3\mathrm{\Δ}x+P_5\mathrm{\Δ}y\\ y^{\′}=y+P_2+P_4\mathrm{\Δ}x+P_6\mathrm{\Δ}y\end{array}\right.$

式中，变形参数向量P＝(P₁,P₂,P₃,P₄,P₅,P₆)^T＝(u,v,u_x,v_x,u_y,v_y)^T，u和v分别为x和y方向的位移，(u_x,v_x,u_y,v_y)为位移梯度；Δx＝x-x₀，Δy＝y-y₀，(x₀,y₀)为参考图像子区的中心位置；

二阶形函数通过下式表达：

$\left\{\begin{array}{c}{x}^{\′}=x+P_1+P_3\mathrm{\Δ}x+P_5\mathrm{\Δ}y+P_7{\mathrm{\Δx}}^2+P_9\mathrm{\Δ}x\mathrm{\Δ}y+P_{11}{\mathrm{\Δy}}^2\\ y^{\′}=y+P_2+P_4\mathrm{\Δ}x+P_6\mathrm{\Δ}y+P_8{\mathrm{\Δx}}^2+P_{10}\mathrm{\Δ}x\mathrm{\Δ}y+P_{12}{\mathrm{\Δy}}^2\end{array}\right.$

式中，变形参数向量

P＝(P₁,P₂,P₃,P₄,P₅,P₆,P₇,P₈,P₉,P₁₀,P₁₁,P₁₂)^T＝(u,v,u_x,v_x,u_y,v_y,u_xx,v_xx,u_xy,v_xy,u_yy,v_yy)^T，u和v分别为x和y方向的位移，(u_x,v_x,u_y,v_y)为位移梯度，(u_xx,v_xx,u_xy,v_xy,u_yy,v_yy)为位移场的二阶偏导数项；Δx＝x-x₀，Δy＝y-y₀，(x₀,y₀)为参考图像子区的中心位置。

更进一步的，步骤(8)中统计得到参考图像子区与目标图像子区中对应区域的灰度差的直方图后，利用下式计算得到带宽：

σ＝sum(P*I)

式中，P和I分别表示子区灰度差值落入直方图不同亮度区间的概率和每个亮度区间的中心亮度值。鲁棒数字图像相关方法中，带宽参数σ的选择很关键，一方面σ控制着被剔除点的数目，另一方面σ关系到算法迭代的收敛性。由于实际实验时“坏点”的数目和位置，以及其引起的灰度差的幅值未知，因此无法根据先验知识确定带宽参数的值。带宽参数σ是用来评估子区中各点灰度差的有效与否，因此可以根据统计规律利用灰度直方图来统计各个灰度差级别出现的概率。

更进一步的，步骤(9)中利用Newton-Raphson方法计算得到令基于相关熵准则的相关函数C^*(P)最大时的变形参数向量P，包括以下步骤：

a、构建迭代等式：

$P=P_0-\frac{\▿C^{*}\left(P_0\right)}{\▿\▿C^{*}\left(P_0\right)}$

式中，P₀为变形参数初值，和是相关函数C^*的一阶梯度和Hessian矩阵，其中和的公式如下：

$\begin{array}{c}\▿C^{*}={\left(\frac{\∂C^{*}}{\∂P_i}\right)}_{i=1,\mathrm{...},6}\\ =\frac{-2}{c}\underset{x=-M}{\overset{M}{\Σ}}\underset{y=-M}{\overset{M}{\Σ}}{\left\{\begin{array}{c}\[f(x,y)-g(x,y;P)\]\frac{\∂g(x,y;P)}{\∂P_i}\\ \exp (-{\left(f\left(x,y\right)-g\left(x,y;P\right)\right)}^2/2{\mathrm{\σ}}^2)\end{array}\right\}}_{i=1,\mathrm{...},6}\end{array}$

$\begin{array}{c}\▿\▿C^{*}={\left(\frac{{\∂}^2{C}^{*}}{\∂P_i\∂P_j}\right)}_{\underset{j=1,\mathrm{...},6}{i=1,\mathrm{...},6}}\\ \≈\frac{2}{c}\underset{x=-M}{\overset{M}{\Σ}}\underset{y=-M}{\overset{M}{\Σ}}{\left\{\begin{array}{c}\frac{\∂g(x,y;P)}{\∂P_i}\frac{\∂g(x,y;P)}{\∂P_j}\\ \exp (-{\left(f\left(x,y\right)-g\left(x,y;P\right)\right)}^2/2{\mathrm{\σ}}^2)\end{array}\right\}}_{\underset{j=1,\mathrm{...},6}{i=1,\mathrm{...},6}}\end{array}$

式中，P_i和P_j分别为变形参数向量P的第i个和第j个元素；

b、利用迭代等式计算得到令相关函数C^*(P)最大时的变形参数向量P。

有益效果：本发明的基于相关熵准则的鲁棒数字图像相关方法相对于传统的数字图像相关方法可以更好地处理非连续变形，并且较传统数字图像相关法而言计算精度更高，抗噪性能更强。

附图说明

图1本发明的基于相关熵准则的鲁棒数字图像相关方法的计算流程图；

图2变形前后图像子区亮度差分布直方图；

图3模拟相机噪声情况下，SDIC和RDIC计算变形的位移误差对比图；

图4模拟相机噪声情况下，SDIC和RDIC计算变形的位移标准差对比图；

图5试件非连续变形示意图；

图6 SDIC计算得到非连续变形位移场；

图7 RDIC计算得到非连续变形位移场；

图8 SDIC计算得到非连续变形应变场；

图9 RDIC计算得到非连续变形应变场；

图10 SDIC计算得到含孔试件孔洞周围位移场；

图11 RDIC计算得到含孔试件孔洞周围位移场；

图12含孔试件孔洞边缘子区有无效点标记示意图；

图13人工加权值后，SDIC计算得到含孔试件孔洞周围位移场；

图14碳纤维实验件断裂图；

图15 SDIC计算裂缝附近位移场(左侧裂缝附近)；

图16 RDIC计算裂缝附近位移场(左侧裂缝附近)；

图17 SDIC计算裂缝附近位移场(右侧裂缝附近)；

图18 RDIC计算裂缝附近位移场(右侧裂缝附近)。

具体实施方式

下面结合附图和具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

传统数字图像相关法(Standard Digital Image Correlation，SDIC)

数字图像相关法以试件表面变形前后的图像作为形变的载体，通过匹配变形前后图像获得试件表面变形的位移场。以变形前图像作为参考图像并在其中划分感兴趣点，以感兴趣点为中心划分参考图像子区，将该参考图像子区在变形后图像中移动并搜索，最终找到最为相似的区域，该区域的中心即为感兴趣点在变形后图像中的位置，从而获得变形位移场。子区匹配搜索最相似区域的过程可以通过优化某相似目标函数来实现。传统DIC是基于均方误差准则构建的二次函数：

$C\left(P\right)=\frac{1}{c}\underset{x=-M}{\overset{M}{\Σ}}\underset{y=-M}{\overset{M}{\Σ}}{\[f(x,y)-g(x,y;P)\]}^2---\left(1\right)$

式中，变形参数向量P＝(u,v,u_x,v_x,u_y,v_y)^T，u,v分别为x和y方向的位移，其余为位移梯度；c为归一化系数常量；M代表参考图像子区半径；(x,y)为参考图像子区中像素点的坐标，f(x,y)为参考图像子区中像素点(x,y)的灰度值，g(x,y,P)为目标图像子区中与像素点(x,y)相对应的像素点的灰度值。变形后图像中对应位置点的坐标(x',y')可由(x,y)和变形参数向量P表示为

$\{\begin{array}{c}{x}^{\′}=x+P_1+P_3(x-x_0)+P_5(y-y_0)\\ y^{\′}=y+P_2+P_4(x-x_0)+P_6(y-y_0)\end{array}---\left(2\right)$

式(2)中，(x₀,y₀)为参考图像子区的中心位置。(x',y')可能为亚像素位置点，因此，可借助插值法如双三次样条插值法，获得亚像素位置对应的灰度值g(x',y')。

由式(1)可知，相似函数实际上是P＝(u,v,u_x,v_x,u_y,v_y)^T的函数，其解可以利用非线性优化方法获得，如Newton-Raphson法。可以获得参数P的迭代算式如下：

$P=P_0-\frac{\▿C\left(P_0\right)}{\▿\▿C\left(P_0\right)}---\left(3\right)$

式中，P₀为变形参数初值估计，P为迭代值，为式(1)的一阶梯度向量，是式(1)的二阶偏导矩阵，又称作Hessian矩阵，分别表示为

$\begin{array}{c}\▿C={\left(\frac{\∂C}{\∂P_i}\right)}_{i=1,\mathrm{...},6}\\ =\frac{-2}{c}{\{\underset{x=-M}{\overset{M}{\Σ}}\underset{y=-M}{\overset{M}{\Σ}}\[f(x,y)-g(x,y;P)\]\frac{\∂g(x,y;P)}{\∂P_i}\}}_{i=1,\mathrm{...},6}\end{array}---\left(4\right)$

$\begin{array}{c}\▿\▿C={\left(\frac{{\∂}^{2}C}{\∂P_i\∂P_j}\right)}_{\underset{j=1,\mathrm{...},6}{i=1,\mathrm{...},6}}\\ \≈\frac{2}{c}{\left\{\underset{x=-M}{\overset{M}{\Σ}}\underset{y=-M}{\overset{M}{\Σ}}\frac{\∂g(x,y;P)}{\∂P_i}\frac{\∂g(x,y;P)}{\∂P_j}\right\}}_{\underset{j=1,\mathrm{...},6}{i=1,\mathrm{...},6}}\end{array}---\left(5\right)$

根据式(4)-(5)可以看出，参考图像子区中所有的点均以同等的权重参与运算，因此，该算法本身不具备区分“坏点”的功能，因此传统DIC的鲁棒性较差，即易受到“噪声”点的干扰从而影响DIC匹配的精度。

本发明的鲁棒数字图像相关法(Robust Digital Image Correlation，RDIC)

传统DIC中基于均方误差准则(Mean Square Error,MSE)构建相似函数，MSE是一种基于二阶统计量的全局相关准则，其全局特性要求待计算的图像子区具有连续性，否则将会对变形参数的计算带来影响。此外，图像采集系统中的随机噪声也会对子区之间的精确匹配造成干扰。

本文基于信息论中的信息熵和核方法构建随机变量X和Y之间的相似函数：

V_σ(X,Y)＝E[k_σ(X,Y)] (6)

式中，k_σ(.)为核函数，σ为核函数的带宽，E表示k_σ(X,Y)的期望，对于有限的N个数据相关熵可表示为

${\hat{V}}_{N,\mathrm{\σ}}=\frac{1}{N}\underset{i=1}{\overset{N}{\Σ}}{k}_{\mathrm{\σ}}(x_i-y_i)---\left(7\right)$

选取常用高斯核，将输入空间通过核函数映射到高维空间，再结合带宽σ的选择控制去除信号中的大偏差值，从而减小野值点的影响。

基于局部相关熵准则建立鲁棒相关函数：

$C^{*}\left(P\right)=\frac{1}{c}\underset{x=-M}{\overset{M}{\Σ}}\underset{y=-M}{\overset{M}{\Σ}}\exp \[-{\left(f\right(x,y)-g(x,y;P\left)\right)}^2/2{\mathrm{\σ}}^2\]---\left(8\right)$

其中，c为归一化常数，M为子区半径大小，exp()为高斯核函数，σ为高斯核函数的带宽，(x,y)为参考图像子区中像素点的坐标，f(x,y)为参考图像子区中像素点(x,y)的灰度值，g(x,y)为目标图像子区中与像素点(x,y)相对应的像素点的灰度值，P为变形参数向量。

最大化该相关熵准则，建立NR迭代式求得最优变形参数P：

$P=P_0-\frac{\▿C^{*}\left(P_0\right)}{\▿\▿C^{*}\left(P_0\right)}---\left(9\right)$

其中，式(8)的一阶偏导和Hessian矩阵为

$\begin{array}{c}\▿C^{*}={\left(\frac{\∂C^{*}}{\∂P_i}\right)}_{i=1,\mathrm{...},6}\\ =\frac{-2}{c}\underset{x=-M}{\overset{M}{\Σ}}\underset{y=-M}{\overset{M}{\Σ}}{\left\{\begin{array}{c}\[f(x,y)-g(x,y;P)\]\frac{\∂g(x,y;P)}{\∂P_i}\\ \exp (-{\left(f\left(x,y\right)-g\left(x,y;P\right)\right)}^2/2{\mathrm{\σ}}^2)\end{array}\right\}}_{i=1,\mathrm{...},6}\end{array}---\left(10\right)$

$\begin{array}{c}\▿\▿C^{*}={\left(\frac{{\∂}^2{C}^{*}}{\∂P_i\∂P_j}\right)}_{\underset{j=1,\mathrm{...},6}{i=1,\mathrm{...},6}}\\ \≈\frac{2}{c}\underset{x=-M}{\overset{M}{\Σ}}\underset{y=-M}{\overset{M}{\Σ}}{\left\{\begin{array}{c}\frac{\∂g(x,y;P)}{\∂P_i}\frac{\∂g(x,y;P)}{\∂P_j}\\ \exp (-{\left(f\left(x,y\right)-g\left(x,y;P\right)\right)}^2/2{\mathrm{\σ}}^2)\end{array}\right\}}_{\underset{j=1,\mathrm{...},6}{i=1,\mathrm{...},6}}\end{array}---\left(11\right)$

若令K＝exp(-(f(x,y)-g(x,y；P))²/2σ²)，式(10)和式(11)表示为

$\begin{array}{c}\▿C^{*}={\left(\frac{\∂C^{*}}{\∂P_i}\right)}_{i=1,\mathrm{...},6}\\ =\frac{-2}{c}\underset{x=-M}{\overset{M}{\Σ}}\underset{y=-M}{\overset{M}{\Σ}}{\{K\[f(x,y)-g(x,y;P)\]\frac{\∂g(x,y;P)}{\∂P_i}\}}_{i=1,\mathrm{...},6}\end{array}---\left(12\right)$

$\begin{array}{c}\▿\▿C^{*}={\left(\frac{{\∂}^2{C}^{*}}{\∂P_i\∂P_j}\right)}_{\underset{j=1,\mathrm{...},6}{i=1,\mathrm{...},6}}\\ \≈\frac{2}{c}\underset{x=-M}{\overset{M}{\Σ}}\underset{y=-M}{\overset{M}{\Σ}}{\left\{K\frac{\∂g(x,y;P)}{\∂P_i}\frac{\∂g(x,y;P)}{\∂P_j}\right\}}_{\underset{j=1,\mathrm{...},6}{i=1,\mathrm{...},6}}\end{array}---\left(13\right)$

与式(4)和式(5)相比，相当于在子区亮度差上添加了权值K。权值的大小与对应位置点上的灰度差以及带宽的大小有关。对于子区中出现的野值点，这些点对应的灰度差较大，超出了带宽允许的通过范围，则这些野值点对应的灰度差被赋予很小的权值，降低其对迭代过程的影响。该做法从原理上保证了算法的抗噪性和处理试件不连续变形的可行性。

鲁棒数字图像相关方法中，带宽参数σ的选择很关键，一方面σ控制着被剔除点的数目，另一方面σ关系到算法迭代的收敛性。由于实际实验时“坏点”的数目和位置，以及其引起的灰度差的幅值未知，因此无法根据先验知识确定带宽参数的值。带宽参数σ是用来评估子区中各点灰度差的有效与否，因此可以根据统计规律利用灰度直方图来统计各个灰度差级别出现的概率。图像的灰度范围在0-255，所以子区中各点的灰度差也在0-255的范围内，灰度差的统计直方图如图2所示。带宽参数σ可由式(13)计算：

σ＝sum(P*I) (13)

式中，P和I分别表示亮度差值落入不同亮度范围的概率和每个亮度区间的中心亮度值。从图2中可以看到，灰度差主要集中在50以内，而对于不含噪声的图像，初始迭代时灰度差均在20以内，图中灰度差的分布也说明即使受噪声干扰，子区中大多数点的灰度差集中在较小的偏差值内，只有少数点受到较大随机噪声的干扰，σ控制将这些受噪声干扰严重的点剔除掉，而剩余的点可保证迭代方向的正确性和迭代过程的收敛性。

仿真分析实例1：模拟相机噪声

为评估鲁棒数字图像相关法的抗噪性能，在模拟散斑图中添加噪声来进行分析。考虑到高斯和椒盐噪声为图像系统中最为常见的两种噪声源，特在此仿真生成的散斑图中添加1％等级的高斯噪声和密度d＝0.001的椒盐噪声。散斑图的大小为500×500pixels，散斑颗粒半径为3，散斑颗粒数目为4000。对于变形后图像，沿y方向生成sin变形图像，满足v＝Asin(2πy/T),A＝1,T＝200。则分别利用传统数字图像相关方法(SDIC)和鲁棒数字图像相关方法(RDIC)来计算分析这组受噪声污染的sin变形散斑图，计算位移场的误差和标准差分别如图3和4所示。

从图3可以看到，RDIC计算得到的位移场的误差较SDIC小，前者最大误差在0.04像素，而后者最大误差在0.08像素，并且SDIC计算的位移场误差较为波动，这直接反应在位移场标准差SDIC的较大。可见，RDIC对于图像噪声的抗干扰能力较强，相比于SDIC，RDIC具有更高的计算精度和计算稳定性。

仿真分析实例2：模拟非连续变形

仿真生成的散斑图大小为800×400像素，变形分为两部分，上半部分沿X向发生2000微应变的均匀变形，下半部分X向发生600微应变的均匀变形，变形示意图如图5所示。分别利用SDIC和RDIC计算这组不连续变形散斑图，计算得到的位移场和应变场分布分别如图6-7以及图8-9所示。

对于散斑图中均匀变形部分，RDIC计算的位移场与SDIC计算效果类似，表明两种方法均可以较好地处理连续变形情况。对于由2000微应变到600微应变不连续变形的过渡带，RDIC的过渡带较窄，即RDIC具有区分不同变形状态的能力，而SDIC由于“全局”计算的限制，过渡带中不同变形的存在使得匹配过程出现偏差，因而过渡带较宽。对于从位移场计算得到的应变场，RDIC的应变场计算误差较小。可见，RDIC计算得到的较高精度的位移场为后续应变场准确的计算提供保障。对于不连续变形问题，RDIC具有自动识别子区中不满足连续变形约束的点的能力，从而降低其对匹配过程的干扰，提高不连续变形边界点的位移场计算精度。

实验分析实例1：含孔试件孔洞边缘变形分析

本组实验中对碳纤维含孔试件进行单向拉伸实验，通过在碳纤维试件的表面喷制人工散斑构造随机特征，实验过程中，以0.3mm/min位移加载方式对试件进行y向拉伸，同时以5帧/s的速度采集试件变形过程中的图片，采集图片结束后，分别以SDIC和RDIC对碳纤维试件空洞周围区域进行计算分析，计算得到的位移场如图10和11所示。

计算出位移场后发现，SDIC计算得到的结果在孔洞左下部分出现计算错误点，而由RDIC计算得到的结果过度平滑，分布规律符合实际孔洞变形的规律。经分析，位于孔洞边缘的点进行匹配计算时，子区中包含了非试件上的孔洞中的背景点，这部分点已不满足连续变形的约束条件，因此影响了计算结果。而RDIC由于其自身具有去除这部分点的能力，所以孔洞中的点不会对计算结果产生影响。

为验证该思路，现通过人为赋权值的方式，将SDIC子区中落在孔洞中的点赋予0权值，而其他试件上的点赋予1权值，人为进行背景点剔除的操作，该操作过程可由图12表示。对同样的数据点利用加权SDIC的方式计算，结果如图13所示，可以明显地看到，原来计算错误的点得到了纠正，证明了确实是由于孔洞中的点导致了计算结果的偏差。加权SDIC的计算结果与RDIC的计算结果分布规律完全吻合，这也证明了RDIC具有自动区分“坏点”并添加权值将坏点自动去除的能力。

实验分析实例2：含孔试件发生拉伸断裂后裂纹周围变形分析

碳纤维含孔试件在y向拉伸过程中发生断裂，实际采集得到的断裂效果如图14所示。碳纤维材料属于脆性材料，因此在拉伸一段时间后，碳纤维材料直接沿着孔洞截面最小受应力最大的方向发生断裂，左侧未完全裂开，裸露出里面的纤维结构，而右侧则完全断开，并且裂纹附近有少许区域表面喷制的散斑脱落。分别以RDIC和SDIC处理这种实际发生的不连续变形裂纹周围的变形状态，结果如图15-18所示。

图中叉叉表示计算迭代不收敛的点，箭头的指向和长短分别表示y向位移场的方向和大小。计算时左右两侧共计算7455点，SDIC有1043个点不收敛，RDIC有744个点不收敛。不收敛的点基本集中在完全断裂区域或由于撕裂而裸露碳纤维试件内部纤维结构上，这些区域由于严重“变形”已不具备对应可匹配的特征，所以迭代不收敛是合理的。但是相比于SDIC，RDIC可以较好的处理试件断裂边缘的点，或是由于撕裂导致表面部分小区域表面散斑脱落的点，以及圆孔边缘的点。只要子区中存在足够的匹配特征，无论子区中是否存在“坏点”，RDIC都可以将这些坏点剔除，然后利用剩下的优质特征点进行匹配搜索，从而实现正确的定位，获得准确的位移场结果分布。

Claims (8)

1.一种基于相关熵准则的鲁棒数字图像相关法，其特征在于：采用基于相关熵准则的相关函数评价变形前后图像子区的相似程度，所述基于相关熵准则的相关函数通过下式表达：

$C^{*}\left(P\right)=\frac{1}{c}\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-M}{\overset{M}{\mathrm{\Σ}}}\exp \[-{\left(f\left(x,y\right)-g\left(x,y;P\right)\right)}^2/2{\mathrm{\σ}}^2\]$

其中，c为归一化常数，M为参考图像子区的半径，exp()为高斯核函数，σ为高斯核函数的带宽，(x,y)为参考图像子区中像素点的坐标，f(x,y)为参考图像子区中像素点(x,y)的灰度值，g(x,y)为目标图像子区中与像素点(x,y)相对应的像素点的灰度值，P为变形参数向量。

2.如权利要求1所述的基于相关熵准则的鲁棒数字图像相关方法，其特征在于：包括以下步骤：

(1)、采集被测件表面变形前和变形后的散斑图像，变形前的散斑图像称为参考图像，变形后的散斑图像称为目标图像；

(2)、在参考图像上划分感兴趣的网格点；

(3)、在参考图像上，建立以步骤(2)的网格点为中心的第一矩形区域，称为参考图像子区，第一矩形区域的长度为L，宽度为W；

(4)、在目标图像上，建立以步骤(2)的网格点为中心的第二矩形区域，第二矩形区域的长度为k*L，宽度为k*W，其中，k>1，称为目标图像子区，目标图像子区包含变形后的参考图像子区；

(5)、获得参考图像子区中整像素位置的灰度值和目标图像子区中整像素位置以及亚像素位置的灰度值；

(6)、定义形函数：设定参考图像子区中心点为(x₀,y₀)，参考图像子区中的任意点为(x,y)，变形后目标图像子区中与点(x,y)对应的点为(x',y')，存在一组映射关系ξ使得下式成立：

ξ(x,y)→(x',y')

f(x,y)＝g(x',y')

其中，f(x,y)表示变形前像素点(x,y)处的图像灰度值，g(x',y')表示变形后对应像素点(x',y')处的图像灰度值；

(7)、将步骤(6)的形函数参数化，并用变形参数向量P表示，定义基于相关熵准则的相关函数：

$C^{*}\left(P\right)=\frac{1}{c}\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-M}{\overset{M}{\mathrm{\Σ}}}\exp \[-{\left(f\left(x,y\right)-g\left(x,y;P\right)\right)}^2/2{\mathrm{\σ}}^2\]$

其中，c为归一化常数，M为参考图像子区的半径，exp()为高斯核函数，σ为高斯核函数的带宽，(x,y)为参考图像子区中像素点的坐标，f(x,y)为参考图像子区中像素点(x,y)的灰度值，g(x,y)为目标图像子区中与像素点(x,y)相对应的像素点的灰度值，P为变形参数向量；

(8)、利用灰度直方图计算得到高斯核函数的带宽σ；

(9)、将步骤(8)得到的高斯核函数的带宽σ代入所述基于相关熵准则的相关函数，计算得到令基于相关熵准则的相关函数C^*(P)最大时的变形参数向量P。

3.如权利要求2所述的基于相关熵准则的鲁棒数字图像相关法，其特征在于：步骤(1)中被测件表面的散斑为自然纹理特征的散斑或人工制备的散斑。

4.如权利要求2所述的基于相关熵准则的鲁棒数字图像相关法，其特征在于：步骤(5)中通过插值方法获得目标图像子区亚像素位置的灰度值。

5.如权利要求4所述的基于相关熵准则的鲁棒数字图像相关法，其特征在于：所述插值方法为最近邻插值法、双线性插值法或双三次样条插值法。

6.如权利要求2所述的基于相关熵准则的鲁棒数字图像相关法，其特征在于：步骤(6)中的映射关系ξ为一阶形函数或二阶形函数，

其中，一阶形函数通过下式表达：

$\left\{\begin{array}{c}{x}^{\′}=x+P_1+P_3\mathrm{\Δ}x+P_5\mathrm{\Δ}y\\ y^{\′}=y+P_2+P_4\mathrm{\Δ}x+P_6\mathrm{\Δ}y\end{array}\right.$

式中，变形参数向量P＝(P₁,P₂,P₃,P₄,P₅,P₆)^T＝(u,v,u_x,v_x,u_y,v_y)^T，u和v分别为x和y方向的位移，(u_x,v_x,u_y,v_y)为位移梯度；Δx＝x-x₀，Δy＝y-y₀，(x₀,y₀)为参考图像子区的中心位置；

二阶形函数通过下式表达：

$\left\{\begin{array}{c}{x}^{\′}=x+P_1+P_3\mathrm{\Δ}x+P_5\mathrm{\Δ}y+P_7\mathrm{\Δ}{x}^2+P_9\mathrm{\Δ}x\mathrm{\Δ}y+P_{11}\mathrm{\Δ}{y}^2\\ y^{\′}=y+P_2+P_4\mathrm{\Δ}x+P_6\mathrm{\Δ}y+P_8\mathrm{\Δ}{x}^2+P_{10}\mathrm{\Δ}x\mathrm{\Δ}y+P_{12}\mathrm{\Δ}{y}^2\end{array}\right.$

式中，变形参数向量

P＝(P₁,P₂,P₃,P₄,P₅,P₆,P₇,P₈,P₉,P₁₀,P₁₁,P₁₂)^T＝(u_,v_,u_x,v_x,u_y,v_y,u_xx,v_xx,u_xy,v_xy,u_yy,v_yy)^T，u和v分别为x和y方向的位移，(u_x,v_x,u_y,v_y)为位移梯度，(u_xx,v_xx,u_xy,v_xy,u_yy,v_yy)为位移场二阶偏导数项；Δx＝x-x₀，Δy＝y-y₀，(x₀,y₀)为参考图像子区中心位置。

7.如权利要求2所述的基于相关熵准则的鲁棒数字图像相关法，其特征在于：步骤(8)中统计得到参考图像子区与目标图像子区中对应区域的灰度差的直方图后，利用下式计算得到带宽：

σ＝sum(P*I)

式中，P和I分别表示子区灰度差值落入直方图不同亮度区间的概率和每个亮度区间的中心亮度值。

8.如权利要求2所述的基于相关熵准则的鲁棒数字图像相关法，其特征在于：步骤(9)中利用Newton-Raphson方法计算得到令基于相关熵准则的相关函数C^*(P)最大时的变形参数向量P，包括以下步骤：

a、构建迭代等式：

$P=P_0-\frac{\▿C^{*}\left(P_0\right)}{\▿\▿C^{*}\left(P_0\right)}$

式中，P₀为变形参数初值，和是相关函数C^*的一阶梯度和Hessian矩阵，其中和的公式如下：

$\begin{array}{c}\▿C^{*}={\left(\frac{\∂C^{*}}{\∂P_i}\right)}_{i=1,...,6}\\ =\frac{-2}{c}\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-M}{\overset{M}{\mathrm{\Σ}}}{\left\{\begin{array}{c}\[f\left(x,y\right)-g\left(x,y;P\right)\]\frac{\∂g\left(x,y;P\right)}{\∂P_i}\\ \exp \left(-{\left(f\left(x,y\right)-g\left(x,y;P\right)\right)}^2/2{\mathrm{\σ}}^2\right)\end{array}\right\}}_{i=1,...,6}\end{array}$

$\begin{array}{c}\▿\▿C^{*}={\left(\frac{{\∂}^2{C}^{*}}{\∂P_i\∂P_j}\right)}_{\underset{j=1,...,6}{i=1,...,6}}\\ \≈\frac{2}{c}\underset{x=-M}{\overset{M}{\mathrm{\Σ}}}\underset{y=-M}{\overset{M}{\mathrm{\Σ}}}{\left\{\begin{array}{c}\frac{\∂g\left(x,y;P\right)}{\∂P_i}\frac{\∂g\left(x,y;P\right)}{\∂P_j}\\ \exp \left(-{\left(f\left(x,y\right)-g\left(x,y;P\right)\right)}^2/2{\mathrm{\σ}}^2\right)\end{array}\right\}}_{\underset{j=1,...,6}{i=1,...,6}}\end{array}$

式中，P_i和P_j分别为变形参数向量P的第i个和第j个元素；

b、利用迭代等式计算得到令函数C^*(P)最大的变形参数向量P。