CN104700368B - 基于核函数的数字图像相关方法的位移场自适应平滑方法 - Google Patents

Abstract

本发明公开了一种基于核函数的数字图像相关方法的位移场自适应平滑方法，采用基于核函数的相关函数ρ^*，并利用自适应平滑方法对基于核函数的数字图像相关方法得到的位移场进行处理。本发明的基于核函数的数字图像相关方法的位移场自适应平滑方法可以有效改进位移场的测量精度，定义了一种新的基于核函数的相关函数ρ^*，考虑图像噪声的影响，使得本发明较传统方法而言，具有更好的测量精度。适应平滑方法基于惩罚最小二乘回归技术，对噪声位移场表面的粗糙性进行惩罚，实现有效地消除位移场中的噪声的目的。

Description

基于核函数的数字图像相关方法的位移场自适应平滑方法

技术领域

本发明属于数字图像处理领域，具体涉及一种基于核函数的数字图像相关方法的位移场自适应平滑方法。

背景技术

数字图像相关，最早由Sutton et al.提出，是一种将物体表面随机分布的自然纹理或人工散斑场作为变形信息的载体，通过机器视觉技术获得结构表面在外载荷作用下的全场位移和应变信息的光学测量方法。由于DIC具有诸多优点，例如测量过程简单、测量结果准确性高、非接触、可获得全场数据等，近年来在实验力学领域获得了广泛应用。

作为一种有效的光学测量方法，DIC技术在近二十年中获得了迅速的发展，并在试验力学，材料及其相关领域获得了大量成功应用。这些都归因于DIC方法本身的种种优势，例如非接触式、试验系统搭建便利、可获得全场数据、以及相对较高的测量准确性，等。为此，大量方法被提出用于改进DIC的测量精度和处理速度，推动该方法在更在更大范围内获得成功应用。不幸的是，尽管基于DIC的测量理论近年来获得了充分的研究，传统方法仍然存在众多问题，在面对实际应用时。散斑模式、图像子区尺寸、相关函数、迭代条件等诸多因素会对测量结果噪声影响。对于此类问题，已有一些学者对DIC性能进行过较为细致的分析，并提出了一些指导实际应用的总体性原则。

必须指出，现实应用中图像噪声是不可避免的。而传统DIC方法的测量过程又严重依赖于图像的亮度信息。当噪声存在时，变形前后用于匹配的图像子区中的亮度分布也会随着噪声的分布而发生改变，从而严重影响着测量结果的准确性。为降低图像噪声对测量结果的影响，一些基于前处理或后处理技术来消除测量偏差的方法被提出。尽管获得了一定的效果，但当用于实际测量时，这些方法仍面临缺少平滑参数如何选取的问题。

DIC的基本原理非常简单，记录结构变形的序列图像，并在变形前后的图像子区中搜索最大化某种相关性准则，例如：零均值归一化交叉相关准则，来获取测量点的位移数据。为了提高位移场的计算精度，研究主要关注于通过改善DIC算法来输出高精度的亚像素位移。Lu et al.将高阶梯度引入形函数以实现对复杂变形的描述，从而改善DIC对复杂变形场的测量效果。Cofaru et al.提出用非规则散斑模式来修正DIC，并结合正则化方法来增加位移场的输出精度。Pan的研究表明：采用大小为5*5像素的高斯低通滤波器对散斑图像进行预处理，可有效地降低DIC输出的位移场误差。

在实验力学领域，相比单纯的位移场数据而言，应变场分布信息显得更有价值。尽管DIC技术经历了多年的发展，由于现实条件限制，DIC输出的位移场仍会存在各种偏差，有关DIC系统偏差的分析详见文献。由于计算应变场的所有信息都包含在位移数据当中，若直接利用这些包含噪声的位移数据来计算应变，误差将会被放大，以至于有效的应变场分布规律难以获取。正因如此，基于数据平滑或曲面拟合的后处理技术被用于消除位移场的噪声。然而，这些方法的缺陷是需要人工调整算法的参数。在实际测量过程中，由于没有足够的先验知识来指导参数调整，因而阻碍了其实际应用。

因此，需要一种新的适用于数字图像相关的位移场迭代平滑方法以解决上述问题。

发明内容

本发明的目的是针对现有技术的数字图像相关方法及其位移场平滑方法的不足，提供一种基于核函数的数字图像相关方法的位移场自适应平滑方法。

为实现上述发明目的，本发明基于核函数的数字图像相关方法的位移场自适应平滑方法可采用如下技术方案：

一种基于核函数的数字图像相关方法的位移场自适应平滑方法，基于核函数的数字图像相关方法采用基于核函数的相关函数ρ^*，

其中，c为归一化常数，k(.)为核函数，h为核函数的带宽控制参数，s_p为变形前图像区域S₀中的像素点，f(s_p)为变形前图像中像素点s_p处的图像亮度，g(s_p,P)为变形后图像中与像素点s_p相对应的像素点处的图像亮度；

其中，自适应平滑方法包括以下步骤：

(1)、测量出变形后的位移场U，其中位移场数据U由下式表示：

其中，表示平滑后的位移场，ζ表示测量过程引入的随机误差；

(2)、构造带有惩罚项的二次函数消除随机误差：

其中，||·||为欧式范数，代表数据逼近程度，C为高阶微分算子，α∈[0,1]代表惩罚因子；

对二次函数求导，并令其导数为零，得到其中，I_n为单位对角矩阵，β＝α/(α+1)，C＝VΛV^-1，其中，V为酉矩阵，满足V^T＝V^-1，Λ为C的特征值构成的对角矩阵，

Λ＝diag(λ₁,λ₂，…λ_i…,λ_n)，其中λ_i＝-2+cos[(i-1)π/n]，

其中，V^T和V分别表示离散余弦变换矩阵和逆余弦变换矩阵；

(3)、利用广义交叉验证方法计算惩罚因子β；其中，广义交叉验证方法通过最小化下式得到惩罚因子β

其中，Tr(.)表示矩阵的迹；其中，

(4)、根据步骤(3)得到的惩罚因子β，根据下式计算得到平滑后的位移场

其中，DCT和IDCT分别表示离散余弦变换和逆余弦变换。

更进一步的，步骤4)中C为二阶拉布拉斯算子。

更进一步的，所述基于核函数的数字图像相关方法包括以下步骤：

一)、定义形函数：设定参考图像中的任意点(x₀,y₀)及其周围的邻域S，(x,y)为参考图像中邻域S中的任意像素点的坐标，为目标图像中与像素点(x,y)相对应的像素点的坐标，存在一组映射关系χ使得下式成立：

其中，f(x,y)表示像素点(x,y)处的图像亮度，表示像素点处的图像亮度，映射关系χ即为形函数；其中，形函数χ由下式表示：

其中，u和v分别为变形引起的x和y方向的位移，(x₀,y₀)为邻域S的中心位置坐标，为变形在x和y方向上产生的一阶位移梯度；

二)、将步骤一)的形函数参数化，并用向量P表示，

定义相关函数ρ^*，

其中，c为归一化常数，k(.)为核函数，h为核函数的带宽控制参数，s_p为邻域S中的一个像素点，f(s_p)为变形前图像中像素点s_p处的图像亮度，g(s_p,P)为变形后图像中与像素点s_p相对应的像素点处的图像亮度；

三)、计算得到令相关函数ρ^*最小时的形函数参数P。

更进一步的，步骤三)中通过下式计算相关函数ρ^*的最小值，构造迭代等式：

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

利用迭代等式计算得到令相关函数ρ^*的最小时的解。

更进一步的，所述基于核函数的数字图像相关方法包括以下步骤：

1)、采集物体变形前和变形后的图像，变形前的图像为参考图像，变形后的图像为目标图像；

2)、在参考图像上选取兴趣区域S₁，兴趣区域S₁为参考图像子区；

3)、在目标图像上建立搜索区域S₂，搜索区域S₂为目标图像子区，目标图像子区包含变形后的参考图像子区；

4)、获取目标图像子区和参考图像子区中的点的灰度值；

5)、构建形函数，确定参考图像子区与目标图像子区中对应点的位置关系，其中，(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₂)处的图像亮度，映射关系χ即为形函数；其中，形函数χ由下式表示：

其中，u和v分别为变形引起的x和y方向的位移，(x',y')为参考图像子区的中心位置坐标，为变形在x₁和y₁方向上产生的一阶位移梯度；

6)、将步骤5)的形函数参数化，并用向量P'表示，

定义相关函数ρ^*，

其中，c为归一化常数，k(.)为核函数，h为核函数的带宽控制参数，s_p'为参考图像子区中的一个像素点，f(s_p')为变形前图像中像素点s_p'处的图像亮度，g(s_p',P')为变形后图像中与像素点s_p'相对应的像素点处的图像亮度；

7)、设定向量的初始值P₀'，并设定最大迭代次数n，

8)、将P₀'带入迭代等式：

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

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

根据迭代等式计算得到P₁'；

9)、根据下式判断迭代等式是否收敛，如果收敛，则保存兴趣区域S₁的位移值u'和v'；如果没有收敛，令P₀'＝P₁'，并重复步骤8)和9)；如果迭代次数达到最大迭代次数n，则结束迭代；

10)、输出兴趣区域的位移值。

更进一步的，步骤4)中获取目标图像子区中的点的灰度值通过对目标图像子区进行双三次样条插值，获得目标图像子区中亚像素位置的灰度值。

更进一步的，步骤4)中获取参考图像子区中的点的灰度值通过对参考图像子区进行双三次样条插值，获得参考图像子区中亚像素位置的灰度值。

更进一步的，步骤2)中在参考图像上选取兴趣区域S₁，兴趣区域S₁为参考图像子区通过以下步骤得到：在参考图像上分别划分第一网格，第一网格的个数为M行×N列，在参考图像的长度和宽度方向上第一网格之间的间距分别为l和w，在参考图像上建立以第a个第一网格为中心的兴趣区域S₁，兴趣区域S₁为参考图像子区，参考图像子区的长度和宽度分别为L和W，其中，a＝1,2,3…,M×N。利用在参考图像上通过划分网格的方式选取兴趣区域，简单方便，容易实现。

更进一步的，步骤3)中在目标图像上建立搜索区域S₂，搜索区域S₂为目标图像子区，目标图像子区包含变形后的参考图像子区通过以下步骤得到：在目标图像上划分第二网格，第二网格的个数为M行×N列，在目标图像的长度和宽度方向上第二网格之间的间距分别为l和w；在目标图像上建立以第b个第二网格为中心的搜索区域S₂，搜索区域S₂为目标图像子区，目标图像子区的长度和宽度分别为k*L和k*W，其中，k>1，b＝a。利用在目标图像上通过划分网格的方式选取搜索区域，简单方便，容易实现。

有益效果：本发明的基于核函数的数字图像相关方法的位移场自适应平滑方法可以有效改进位移场的测量精度。定义了一种新的基于核函数的相似性系数，考虑图像噪声的影响，使得本发明较传统方法而言，具有更好的测量精度。自适应平滑方法基于惩罚最小二乘回归技术，对噪声位移场表面的粗糙性进行惩罚，实现有效地消除位移场中的噪声的目的。此外，GCV方法和离散余弦变换分别用于从噪声位移场中估计惩罚因子和真实位移场。该方法具有实现简单、计算量小、全自动的优点。

附图说明

图1、是本发明的基于核函数的数字图像相关方法的流程图；

图2、变形前参考图像；

图3、变形后图像(添加4％噪声)；

图4、NR方法计算得到的四组位移误差对比图；

图5、NR方法计算得到的四组位移标准差对比图；

图6、基于核函数的数字图像相关方法的位移场自适应平滑方法计算得到的四组位移误差对比图；

图7、基于核函数的数字图像相关方法的位移场自适应平滑方法计算得到的四组位移标准差对比图；

图8、为试件图；

图9、为基于核函数的数字图像相关方法的位移场自适应平滑方法的处理效果图；

图10、为试件处理原理图；

图11是本发明的位移场自适应平滑方法的流程图；

图12是均匀变形情况下利用常规方法处理DIC输出数据的处理结果；

图13是均匀变形情况下利用本发明的位移场自适应平滑方法处理DIC输出数据得到的结果；

图14是非均匀变形情况下利用常规方法处理DIC输出数据的处理结果；

图15是非均匀变形情况下利用本发明的位移场自适应平滑方法处理DIC输出数据得到的结果；

图16为模拟散斑图像；

图17为散斑图；

图18为散斑原理图；

图19为多孔试件利用常规DIC方法处理得到的拉伸应变场变形测量结果；

图20为多孔试件利用常规DIC方法处理得到的剪切应变场变形测量结果；

图21为多孔试件利用本发明的方法处理得到的拉伸应变场变形测量结果；

图22为多孔试件利用本发明的方法处理得到的剪切应变场变形测量结果。

具体实施方式

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

传统DIC方法

DIC方法的基本原理非常简单，即将结构变形的测量问题转化为试件变形前(参考图像)、后(目标图像)图像的相关性匹配问题进行求解。因此，为描述结构表面的变形情况，首先需要定义形函数。假定参考图像中的任意点(x,y)及其周围的一个小的邻域S，存在一组映射关系χ满足

其中，f(x,y)表示点(x,y)处的图像亮度，表示变形后点处的图像亮度。

映射χ被称为所谓的形函数。若邻域S和变形量足够小，形函数χ可由式(1)来描述，

其中，u和v分别为变形引起的x和y方向的面内位移，(x₀,y₀)为区域S的中心位置坐标。

将形函数写成向量形式，

并定义相关系数，

然后借助于非线性优化方法获得令式(2)最小化的最优解，问题便获得解决。

从式(2)可以看出，当相关函数取极小值时，变形前后图像子区的相似性达到最大值。此时，参数向量P包含的位移参数u和v代表对变形后位移的最佳估计，用同样的方式对所有的测量点进行计算，即可得到全场位移。

最小化ρ有多种求解方法，据作者所知，Newton-Raphson方法因计算精度较高而被众多文献所采用，即构造如下迭代等式

式中，P₀为变形参数初值，▽ρ和▽▽ρ是相关函数ρ的一阶梯度和Hessian矩阵，满足

和

仔细分析式(3-5)可以看出，传统DIC在最小化ρ的过程中，f(s_p)-g(s_p,P)和两项对迭代结果性影响很大，甚至会影响迭代是否收敛。理想情况下(没有噪声)，每迭代一次，形函数的参数都会沿梯度方向向真实变形参数趋近一步。当迭代满足终止条件时，应有趋向于零。然而，当图像噪声存在时，不但f(s_p)-g(s_p,P)项不等于零，项也会因为噪声的影响而放大失真。因此，传统DIC在实际应用时会因噪声的存在而降低测量效果。

本发明的DIC方法

为改善传统DIC对噪声的适应性，本文将核函数引入到相似性函数的定义。核函数通常定义为某种沿径向对称的标量函数。例如，高斯核函数，即所谓的径向基函数(RadialBasis Function,RBF)，定义为空间中任一点x到某一中心x_c之间欧氏距离的单调函数，可记作k(||x-x_c||)，形式为

其中，x_c为核函数中心，σ为函数的宽度参数，控制函数的径向作用范围，其作用是局部的，当x远离x_c时函数取值很小。

为利用核函数在径向范围内的控制能力，在给定的邻域S内，通过将变形前后的图像亮度做差的方式，将变形前后图像子区的相似性度量变换到以零为中心的特征空间，并以平方和误差的形式表示，构造出基于核函数的相关系数

其中，c为归一化常数，k(.)为核函数，h为核函数的带宽控制参数。

为简化表达，令

r(s_p,p)＝f(s_p)-g(s_p,p)

式(4)简化为，

ρ的一阶梯度和Hessian矩阵为，

同理，可将式(8)和式(9)带入式(3)，并借助迭代过程获得最小化式(7)的数值解。

两种方法的对比

比较两种方法相似性函数的一阶梯度和Hessian矩阵可以看出，形式上KDIC为TDIC的加权版本，而权函数等价于核函数的一阶导数。因此，不同核函数的类型将会对应不同的加权版本。为直观分析，以最常用的两种核函数为例进行讨论。一种是Epanechnikov核，另一种是高斯核。

Case1：Epanechnikov核

即，

若将Epanechnikov核带入式(8)和(9)，从表达形式上看，KDIC与TDIC是完全等价的。也就是说，采用Epanechnikov核时，TDIC是KDIC的一个特例。

Case2：高斯核

以及

由于高斯函数的导数与原函数具有相同的轮廓，我们可以利用高斯函数的数学特性及其与实际的物理含义之间的关系来实现加权，从而获得对图像噪声的免疫力。假设图像中的噪声为随机分布的白噪声，且方差为σ。当据此，我们可以将核函数中的带宽参数(式(6)中的h)设定为σ的倍数。当偏差值超过该倍数时表明改点处的像素被噪声污染非常严重，分配一个几乎为零的权值可直接消除较大噪声点影响。相反，当偏差值落在该倍数范围内时，通过高斯加权的方式根据噪声污染程度对不同点的偏差值进行加权，作为计算区域中的各个像素点分配贡献强度，从而获得更加准确的测量结果。

实施例1:

基于核函数的数字图像相关方法的位移场自适应平滑方法，包括以下步骤：

1)、采集物体变形前和变形后的图像，变形前的图像为参考图像，变形后的图像为目标图像；

2)、在参考图像上选取兴趣区域S₁，兴趣区域S₁为参考图像子区；通过以下步骤得到：在参考图像上分别划分第一网格，第一网格的个数为M行×N列，在参考图像的长度和宽度方向上第一网格之间的间距分别为l和w，在参考图像上建立以第i个第一网格为中心的兴趣区域S₁，兴趣区域S₁为参考图像子区，参考图像子区的长度和宽度分别为L和W，其中，a＝1,2,3…,M×N。利用在参考图像上通过划分网格的方式选取兴趣区域，简单方便，容易实现。

3)、在目标图像上建立搜索区域S₂，搜索区域S₂为目标图像子区，目标图像子区包含变形后的参考图像子区；通过以下步骤得到：在目标图像上划分第二网格，第二网格的个数为M行×N列，在目标图像的长度和宽度方向上第二网格之间的间距分别为l和w；在目标图像上建立以第b个第二网格为中心的搜索区域S₂，搜索区域S₂为目标图像子区，目标图像子区的长度和宽度分别为k*L和k*W，其中，k>1，b＝a。利用在目标图像上通过划分网格的方式选取搜索区域，简单方便，容易实现。其中，l和w均为3-5像素点。

4)、获取目标图像子区和参考图像子区中的点的灰度值。获取目标图像子区中的点的灰度值通过对目标图像子区进行双三次样条插值，获得目标图像子区中亚像素位置的灰度值。获取参考图像子区中的点的灰度值也通过对参考图像子区进行双三次样条插值，获得参考图像子区中亚像素位置的灰度值。

5)、构建形函数，确定参考图像子区与目标图像子区中对应点的位置关系，其中，(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₂)处的图像亮度，映射关系χ即为形函数；其中，形函数χ由下式表示：

其中，u和v分别为变形引起的x和y方向的位移，(x',y')为参考图像子区的中心位置坐标，为变形在x₁和y₁方向上产生的一阶位移梯度；

6)、将步骤5)的形函数参数化，并用向量P'表示，

定义相关函数ρ^*，

其中，c为归一化常数，k(.)为核函数，h为核函数的带宽控制参数，s_p'为参考图像子区中的一个像素点，f(s_p')为变形前图像中像素点s_p'处的图像亮度，g(s_p',P')为变形后图像中与像素点s_p'相对应的像素点处的图像亮度；

7)、设定向量的初始值P₀'，并设定最大迭代次数n，

8)、将P₀'带入迭代等式：

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

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

根据迭代等式计算得到P₁'；

9)、根据下式判断迭代等式是否收敛，如果收敛，则保存兴趣区域S₁的位移值u'和v'；如果没有收敛，令P₀'＝P₁'，并重复步骤8)和9)；如果迭代次数达到最大迭代次数n，则结束迭代；

10)、输出兴趣区域的位移值。

本算法输入的是采集的试件变形前后的散斑图，利用图像的灰度特征，最终输出的是试件变形的位移场；计算时，在变形前试件的感兴趣区域划分网格点，目的就是求取这些网格点在变形后图像的位置，从而可以作差求出这些网格点处的位移。

以划分的网格点为中心划分参考图像子区，以同样的中心点在变形后图像中划分目标图像子区，大小为参考图像子区的k倍，k>1。为了找到变形之后网格点所在的位置，将参考图像子区在目标图像子区中进行逐点搜索，找到最相似的子区，那么该子区的中心点位置即为网格点在变形后图像中的位置。

利用图像灰度特征构建相似性函数，如最小平方距离函数等。为提高该相似性函数的抗噪性，为该函数添加权值项，式中，k可选择Gaussian函数。

式中，f为参考图像子区中各点的灰度值，g为变形后图像子区中的各点的灰度值。变形后图像子区中的点的位置(x',y')由形函数描述，并且(x',y')为亚像素位置值，因此需对目标子区进行插值，如双三次样条插值，来获得亚像素位置的灰度值。P＝[U,V,Ux,Vx,Uy,Vy]^T为形函数中的未知量，实际上，该相似函数是P的函数，最小化该函数，可解得P的值。

相似函数是一个多变量的函数，可由非线性优化方法，如Newton-Raphson方法，进行迭代求解。

迭代等式的迭代收敛条件：可取ε＜0.1，在保证收敛的前提下，ε越小，输出位移场越准确。或设置收敛次数，如50次，在保证收敛的前提下，该值越大，输出位移场越准确。

模拟分析：为评估上述方法的有效性，数字散斑图像模拟测量试件变形的方法被采用。由于变形参数已知，可客观地评价处理效果的优劣。参考图像为随机生成的数字散斑图像，分辨率为256×256pixels，散斑颗粒数为2000，散斑颗粒半径为4，如图2所示。通过预先设定的变形参数，沿y方向，每隔0.5pixel生成一幅变形散斑图，共生成20张。向这组位移为0-1pixel的图像分别添加1％、2％和4％的噪声，加上无噪声的一组散斑图，共获得4组散斑图。其中，添加4％噪声水平的散斑图如图3所示。

为考察不同级别噪声的影响，首先利用传统DIC方法处理四组散斑图，计算得到位移值将采用两种方法进行评价，即位移平均误差和标准差，表达式如下：

其中，表示同一幅图中所有计算点位移的平均值，即v_true代表理论变形参数；N表示所有计算点数。

图4和5分别表示了由NR方法计算位移场后分析得到的误差和标准差。从图4可以看到，随机噪声对位移场计算影响显著。随着噪声等级的增加，位移误差整体呈上升趋势，从最大误差值来看，从0.0077pixel增加到0.0994pixel。

同样，从图5可以看到，所得位移场的标准差也是随着噪声等级的增加而增加，最大标准差从0.0029pixel增加到0.0133pixel。

为验证本文改进DIC方法的有效性，KDIC方法被用来计算四组仿真散斑图。同样，位移平均误差和标准差两种评价指标被用来评估方法的有效性。计算结果如图6和图7所示。

从图7中可以看到，对比传统DIC计算得到的含4％噪声水平的位移场误差，由KDIC方法计算得到的所有噪声水平的误差均较小，其中以4％噪声水平计算误差为例分析，最大误差从0.0094pixel降到0.0271，新方法效果显著。

为便于比较，4组位移计算结果的最大误差和平均标准差被列在表1中。结果证明，基于核函数的DIC方法能有效的降低位移场的计算误差，并且，噪声水平越大效果越显著。但是，从标准差的计算结果来看，处理效果不是很明显。对于1％噪声水平而言，即图像信噪比为40dB，位移场的计算误差在6‰，而实际用于采集图像的工业相机基本可以达到40dB的信噪比甚至更高，因此该方法用于处理实际采集的图像时，千分位的精度完全可以满足要求。

Table 1TDIC和KDIC方法计算结果最大误差和平均标准差对比

以上仿真分析证明了本文改进DIC方法降低位移场计算误差的有效性。为证明该方法可有效处理实际采集的图像，下面将以KDIC方法进行计算分析。实际图片如图8所示，试件为铝制材料，进行y方向拉伸。处理效果如图9所示。

从处理结果可以看到，经KDIC方法计算得到的位移场较为平滑，可有效减少后期由位移场差分得到的应变场计算的误差，该结果证明了本文改进DIC方法用于实际实验分析的可行性和有效性。

请参阅图11所示，本发明的适用于DIC的位移场自适应平滑方法。根据物体表面的光滑性事实，建立惩罚最小二乘回归目标函数，并借助于GCV方法从噪声数据中估计惩罚因子来惩罚解的粗糙性，从而实现位移场的有效平滑处理。该方法具有实现简单、计算量小和全自动的优点。

利用DIC进行结构表面的变形场测量：

将DIC用于结构表面的变形场测量问题，主要包括三个重要环节。首先，构造合适的形函数描述结构在外载荷作用下的变形；然后，建立某种相关性指标，用来定量评价被测结构变形前后的图像亮度分布的相似性程度；最后，通过多变量优化方法求解出令相似性准则最大化的形函数参数，从而间接测量出变形后的位移场数据。

给定未变形图像中的任意点(x,y)及其周围的一个小的邻域S，存在一组映射关系χ满足且f(x,y)表示点(x,y)处的图像亮度，表示变形后相应坐标处的图像亮度。若邻域S足够小，映射关系χ可由式(1)来描述，

其中，u和v分别为x和y方向的面内位移，(x₀,y₀)为区域S的中心位置。

令参数向量并定义相关系数

从上式可以看出，当相关函数取极小值时，变形前后图像子区的相似性达到最大值，此时，参数向量P包含的位移参数u和v代表对变形后位移的最佳估计，用同样的方式对所有的测量点进行计算，即可得到全场位移。

为最小化ρ，可通过求解式(11)一阶梯度为零的解，

有很多方法可用来求解式(12)，本文采用Newton-Raphson方法迭代求解，有

式中，P₀为变形参数初值，▽▽C(P₀)是相关函数ρ的Hessian矩阵，满足

从以上过程可以看出，多种因素会对位移场测量结果的准确性造成影响，例如子区大小、评价准则、迭代收敛条件等。虽然DIC方法得到了大量的研究，但指导消除各种系统测量偏差的有效方法还很缺少。如希望获得有价值的应变场分布数据，需要对含有随机误差的位移场数据进行精心处理。

利用GCV方法对位移场进行平滑处理：

假设U代表DIC输出的位移场数据，从数学意义上说，U可以由如下模型来进行描述，

其中，表示真实的位移场数据，ζ表示测量过程引入的随机误差，一般假设其为满足零均值的高斯分布。

为消除随机误差，可通过最小化如下带有惩罚项的二次函数

其中，||.||为欧式范数，代表数据逼近程度，C一般为高阶微分算子，用来描述解的粗糙型，在本文中取为二阶拉布拉斯算子，α∈[0,1]代表惩罚因子，取值大小代表对粗糙性惩罚的程度。

对上式求导，整理得

其中，I_n为单位对角矩阵，β＝α/(α+1)。

上式虽然可由数值计算方法求解，但非常耗时。如位移场数据为等间隔的，可采用一种有效的方法实现快速求解。实际应用中，图像数据通常以像素为单位，因此等间隔假设是成立的。若假设位移测量点之间的间隔为1，且针对一维数据，C可表示为如下方阵，

对C进行特征值分解，有

C＝VΛV^-1 (18)

其中，V为酉矩阵，满足V^T＝V-¹，Λ表示C的特征值构成的对角矩阵，有

Λ＝diag(λ₁,…,λ_n)withλ_i＝-2+cos[(i-1)π/n]

V^T和V分别表示离散余弦变换矩阵和逆余弦变换矩阵，据此，将式(18)代入式(17)，并结合V^T和V描述的矩阵变换关系，有

其中，DCT和IDCT分别表示离散余弦变换和逆余弦变换。

从式(10)可以看出，惩罚因子β取值(式(16)中的α)对于准确估计至关重要，过大或者过小都会面临过度或欠平滑的风险。本发明选择广义交叉验证(GCV)方法来实现对β的自动估计。

GCV方法通过最小化式(11)来获得理想的惩罚因子β。

其中，Tr(.)表示矩阵的迹。

对于Tr(I_n+βC^TC)项而言，有

而对于项而言，有

综上，我们可以借助式(11)从含噪声的位移场U中先估计出惩罚因子β，再代入式(19)计算出平滑后的位移场估计整个计算过程都非常简洁、高效。实验证明：

结合仿真分析和试验验证来评估上述方法的有效性。首先，在计算机仿真分析中，利用数字散斑图像来模拟测量对象的变形情况。由于变形参数已知，可客观地评价处理效果的优劣。

选择两组变形参数，分别模拟均匀和非均匀变形。数字散斑图像大小为500×500pixels，散斑颗粒数为4000，散斑颗粒半径为4，散斑信噪比为40db。模拟散斑图像参见图16。

实施例1：均匀变形

变形参数取P＝(0,0.3,0.001,0,0,0)^T，这代表在x方向有1000με的均匀拉伸，而其他方向上的应变为零，拉伸后的散斑图像如图17所示。采用式(13)给出的Newton-Raphson方法(传统DIC方法)先计算出位移场，再通过差分运算获得的应变场数据如图12所示，而经过本文平滑处理后的结果在图13给出。

为客观地评价，分别将传统DIC和平滑DIC的处理结果与变形参数的理论计算值进行比较，采用均方根误差公式，

其中，σ为理论值，为估计值，n表示数据总数。

计算结果列于表1，可以看出，无论对于x、y方向的应变，还是剪切应变，本文方法在RMS评价指标上，分别比传统DIC方法提高了39％，48％和52％，表明本文方法的对于均匀变形情况是有效的。

表1 均匀变形的均方根误差对比

实施例2：非均匀变形

实际变形往往是非均匀的，为评估非均匀变形工况下的效果，选择正弦函数u(x,y)＝Asin(2πx/p)来描述结构的变形情况，其中幅值A取0.1，周期参数p取200，其余变形参数同实施例1。计算结果如图14和图15所示，左边是传统DIC输出的位移场数据差分计算出的应变场，右边是平滑DIC计算出的应变场数据。可以看到，传统DIC方法计算得到的应变场较为粗糙，受位移场噪声影响较为严重。而利用本文方法对位移场进行平滑处理后，得到的应变场数据平滑，更加接近于理论值。非均匀变形工况的RMS指标对比数据详见表2，可以看出本文方法在非均匀变形情况下也是有效的。

表2 非均匀变形的均方根误差对比

实施例3：孔洞实验:

请参阅图16、17和18所示，为了进一步验证本文方法的实用性，此处将对实际的试件变形情况进行处理。实验试件为一带孔洞试件，其示意图如图18所示，试件编号为BD1W2L3，其中B代表钢，D＝1为孔洞直径，W＝2为横向孔洞间隔，L＝3为纵向孔洞间隔，试件总长200mm，宽为20mm。测试时，采用工业AVT相机，型号F-125B/C，采集变形前后图像。

传统DIC和平滑DIC输出的应变云图分别在图19、20、21和22中给出。经对比可知，两种方法计算结果总体规律一致，即发生应力集中的位置相吻合。但是，传统DIC方法计算的结果由于受噪声干扰，数据表达的应变分布规律存在较大的失真情况，表现为等高线呈现出锯齿状或不连续，这与试件的实际受力情况不符。而本文方法计算得到的应变场分布更符合实际变形规律，等高线连续且平滑，应变梯度和对称性分布规律与理论值更加接近。

结论：

综上所述,针对传统DIC直接输出的位移场数据，提出了一种自适应的位移场平滑方法，该方法基于GCV技术，能自动地从噪声场数据中计算出惩罚因子，并估计出平滑的位移场数据，从而为应变场的计算提供更加可信的位移场数据。仿真分析和试验分析结果表明，该方法能实现全自动的平滑过程，较传统方法能给出更加合理的应变云图。

Claims (9)

1.一种基于核函数的数字图像相关方法的位移场自适应平滑方法，其特征在于：基于核函数的数字图像相关方法采用基于核函数的相关函数ρ^*，

其中，c为归一化常数，k(.)为核函数，h为核函数的带宽控制参数，s_p为变形前图像邻域S中的一个像素点，f(s_p)为像素点s_p处的图像亮度，g(s_p,P)为变形后图像中与像素点s_p相对应的像素点处的图像亮度；P为形函数参数；

其中，自适应平滑方法包括以下步骤：

(1)、测量出变形后的位移场U，其中，位移场U由下式表示：

其中，表示平滑后的位移场，ξ表示测量过程引入的随机误差；

(2)、构造带有惩罚项的二次函数消除随机误差：

其中，||.||为欧式范数，代表数据逼近程度，C为高阶微分算子，α∈[0,1]代表惩罚因子；

对二次函数求导，并令其导数为零，得到其中，I_n为单位对角矩阵，β＝α/(α+1)，C＝VΛV^-1，其中，V为酉矩阵，满足V^T＝V^-1，Λ为C的特征值构成的对角矩阵，

Λ＝diag(λ₁,λ₂，…λ_i…,λ_n)，其中λ_i＝-2+cos[(i-1)π/n]，

其中，V^T和V分别表示离散余弦变换矩阵和逆余弦变换矩阵；

(3)、利用广义交叉验证方法计算惩罚因子β；其中，广义交叉验证方法通过最小化下式得到惩罚因子β

其中，Tr(.)表示矩阵的迹；其中，

(4)、根据步骤(3)得到的惩罚因子β，根据下式计算得到平滑后的位移场

其中，DCT和IDCT分别表示离散余弦变换和逆余弦变换。

2.如权利要求1所述的基于核函数的数字图像相关方法的位移场自适应平滑方法，其特征在于：步骤(4)中C为二阶拉布拉斯算子。

3.如权利要求1所述的基于核函数的数字图像相关方法的位移场自适应平滑方法，其特征在于：所述基于核函数的数字图像相关方法包括以下步骤：

一)、定义形函数：设定参考图像中的任意点(x₀,y₀)及其周围的邻域S，(x,y)为参考图像中邻域S中的任意像素点的坐标，为目标图像中与像素点(x,y)相对应的像素点的坐标，存在一组映射关系χ使得下式成立：

其中，f(x,y)表示像素点(x,y)处的图像亮度，表示像素点处的图像亮度，映射关系χ即为形函数；其中，形函数χ由下式表示：

其中，u和v分别为变形引起的x和y方向的位移，(x₀,y₀)为邻域S的中心位置坐标，为变形在x和y方向上产生的一阶位移梯度；

二)、将步骤一)的形函数参数化，并用向量P表示，

定义相关函数ρ^*，

其中，c为归一化常数，k(.)为核函数，h为核函数的带宽控制参数，s_p为邻域S中的一个像素点，f(s_p)为像素点s_p处的图像亮度，g(s_p,P)为变形后图像中与像素点s_p相对应的像素点处的图像亮度；

三)、计算得到令相关函数ρ^*最小时的形函数参数P。

4.如权利要求3所述的基于核函数的数字图像相关方法的位移场自适应平滑方法，其特征在于：步骤三)中通过下式计算相关函数ρ^*的最小值，构造迭代等式：

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

其中，P_i和P_j分别为形函数参数P中第i和j个参数；

利用迭代等式计算得到令相关函数ρ^*最小时的解。

5.如权利要求1所述的基于核函数的数字图像相关方法的位移场自适应平滑方法，其特征在于：所述基于核函数的数字图像相关方法包括以下步骤：

1)、采集物体变形前和变形后的图像，变形前的图像为参考图像，变形后的图像为目标图像；

2)、在参考图像上选取兴趣区域S₁，兴趣区域S₁为参考图像子区；

3)、在目标图像上建立搜索区域S₂，搜索区域S₂为目标图像子区，目标图像子区包含变形后的参考图像子区；

4)、获取目标图像子区和参考图像子区中的点的灰度值；

5)、构建形函数，确定参考图像子区与目标图像子区中对应点的位置关系，其中，(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₂)处的图像亮度，映射关系χ即为形函数；其中，形函数χ由下式表示：

其中，u'和v'分别为变形引起的x'和y'方向的位移，(x',y')为参考图像子区的中心位置坐标，为变形在x₁和y₁方向上产生的一阶位移梯度；

6)、将步骤5)的形函数参数化，并用向量P'表示，

定义相关函数ρ^*'，

其中，c为归一化常数，k(.)为核函数，h为核函数的带宽控制参数，s_p'为参考图像子区中的一个像素点，f(s_p')为像素点s_p'处的图像亮度，g(s_p',P')为变形后图像中与像素点s_p'相对应的像素点处的图像亮度；

7)、设定向量的初始值P₀'，并设定最大迭代次数n'，

8)、将P₀'带入迭代等式：

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

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

根据迭代等式计算得到P₁'；

9)、根据下式判断迭代等式是否收敛，如果收敛，则保存兴趣区域S₁的位移值u'和v'；如果没有收敛，令P₀'＝P₁'，并重复步骤8)和9)；如果迭代次数达到最大迭代次数n'，则结束迭代；

10)、输出兴趣区域的位移值。

6.如权利要求5所述的基于核函数的数字图像相关方法的位移场自适应平滑方法，其特征在于：步骤4)中获取目标图像子区中的点的灰度值通过对目标图像子区进行双三次样条插值，获得目标图像子区中亚像素位置的灰度值。

7.如权利要求5所述的基于核函数的数字图像相关方法的位移场自适应平滑方法，其特征在于：步骤4)中获取参考图像子区中的点的灰度值通过对参考图像子区进行双三次样条插值，获得参考图像子区中亚像素位置的灰度值。

8.如权利要求5所述的基于核函数的数字图像相关方法的位移场自适应平滑方法，其特征在于：步骤2)中在参考图像上选取兴趣区域S₁，兴趣区域S₁为参考图像子区通过以下步骤得到：在参考图像上分别划分第一网格，第一网格的个数为M行×N列，在参考图像的长度和宽度方向上第一网格之间的间距分别为l和w，在参考图像上建立以第a个第一网格为中心的兴趣区域S₁，兴趣区域S₁为参考图像子区，参考图像子区的长度和宽度分别为L和W，其中，a＝1,2,3…,M×N。

9.如权利要求5所述的基于核函数的数字图像相关方法的位移场自适应平滑方法，其特征在于：步骤3)中在目标图像上建立搜索区域S₂，搜索区域S₂为目标图像子区，目标图像子区包含变形后的参考图像子区通过以下步骤得到：在目标图像上划分第二网格，第二网格的个数为M行×N列，在目标图像的长度和宽度方向上第二网格之间的间距分别为l和w；在目标图像上建立以第b个第二网格为中心的搜索区域S₂，搜索区域S₂为目标图像子区，目标图像子区的长度和宽度分别为k*L和k*W，其中，k>1，b＝a，L和W分别为参考图像子区的长度和宽度，其中，a＝1,2,3…,M×N。