CN104616272A - 一种适用于数字图像相关的位移场迭代平滑方法 - Google Patents

Abstract

本发明公开了一种适用于数字图像相关的位移场迭代平滑方法，在原最小二乘法的基础上，添加粗糙惩罚项，使解的准确性和解的粗糙性之间达到一个平衡，其中，平滑参数α可以自适应求解，而不用人为地选择该参数，方便、实用，计算输出的平滑位移场，结果较准确。

Description

一种适用于数字图像相关的位移场迭代平滑方法

技术领域

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

背景技术

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

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

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

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

发明内容

本发明的目的是针对现有技术中适用于数字图像相关的位移场平滑方法的不足，提供一种适用于数字图像相关的位移场迭代平滑方法。

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

一种适用于数字图像相关的位移场迭代平滑方法，包括以下步骤：

1)、测量出变形后的位移场数据dⁿ，其中位移场数据dⁿ由下式表示：dⁿ＝d+n；

其中，d为平滑后的位移场数据，n为位移场噪声；

2)、构建二次函数 $D(d,\mathrm{\α})=\mathrm{\α}{\left|\right|d-d^n\left|\right|}_2^2+(1-\mathrm{\α}){\left|\right|\mathrm{Cd}\left|\right|}_2^2$

其中，dⁿ为变形后的位移场数据，d为平滑后的位移场数据，n为位移场噪声，α为平滑参数，C为高阶算子；

求解二次函数的一阶偏导数，得到

$\frac{\∂D(d,\mathrm{\α})}{\∂d}=2\mathrm{\α}(d-d^n)+2(1-\mathrm{\α})*4\left(\mathrm{Cd}\right),$

$\frac{\∂D(d,\mathrm{\α})}{\∂\mathrm{\α}}={\left|\right|d-d^n\left|\right|}_2^2-{\left|\right|\mathrm{Cd}\left|\right|}_2^2$

其中，C为高阶算子，

将d和α变为列向量并令Q＝[d1,d2,d3...dm,α]，其中，m为点数；则可得

$\▿D\left(Q\right)=\left[\begin{array}{c}2\mathrm{\α}(d_1-{d_1}^n)+2(1-\mathrm{\α})*4\left(C{d}_1\right)\\ 2\mathrm{\α}(d_2-{d_2}^n)+2(1-\mathrm{\α})*4\left(C{d}_2\right)\\ \·\·\·\·\·\·\\ 2\mathrm{\α}(d_m-{d_m}^n)+2(1-\mathrm{\α})*4\left({\mathrm{Cd}}_m\right)\\ {\left|\right|d-d^n\left|\right|}_2^2-{\left|\right|\mathrm{Cd}\left|\right|}_2^2\end{array}\right];$

3)、根据步骤2)的▽D(Q)求解二次函数的二阶偏导数▽▽D(Q)，

$\frac{\∂D^2\left(Q\right)}{\∂d_i\∂d_j}=\left[\begin{array}{cccc}32-32\mathrm{\α}& & & \\ & 32-30\mathrm{\α}& & \\ & & \·\·\·\·\·\·& \\ & & & 32-30\mathrm{\α}\end{array}\right]$

$\frac{\∂D^2\left(Q\right)}{\∂d_i\∂\mathrm{\α}}=2(d-d^n)-8\left(\mathrm{Cd}\right)$

$\frac{\∂D^2\left(Q\right)}{\∂\mathrm{\α}\∂d_j}=2(d-d^n)-8\left(\mathrm{Cd}\right)$

$\frac{\∂D^2\left(Q\right)}{\∂\mathrm{\α}\∂\mathrm{\α}}=0$

将上述矩阵组合，得到▽▽D(Q)；

4)、利用步骤2)和3)得到的▽D(Q)和▽▽D(Q)，对下式进行迭代计算：

$P_{k+1}=P_k-\frac{\▿D\left(Q\right)}{\▿\▿D\left(Q\right)},$

其中，Q的初值为[dⁿ,α]^T，其中，dⁿ为变形后的位移场数据，当上式收敛后得到的d即为平滑后的位移值。

更进一步的，所述α的范围为(0,1)。

更进一步的，步骤4)的收敛条件为其中，ε≤0.1。

更进一步的，步骤4)中最大迭代次数N≤50。在保证收敛的前提下，迭代次数N值越大，输出位移场越准确。

有益效果：本发明的适用于数字图像相关的位移场迭代平滑方法在在原最小二乘法的基础上，添加粗糙惩罚项，使解的准确性和解的粗糙性之间达到一个平衡，其中，平滑参数α可以自适应求解，而不用人为地选择该参数，方便、实用，计算输出的平滑位移场，结果较准确。

附图说明

图1是本发明的适用于数字图像相关的位移场迭代平滑方法的流程图；

图2是不同α初值条件下，α迭代变化情况；

图3是N-R迭代计算位移值、平滑后位移值和真实位移值的对比图；

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

图5是非均匀变形情况下利用本发明的适用于数字图像相关的位移场迭代平滑方法处理DIC输出数据得到的结果；

图6为模拟散斑图像；

图7为散斑图；

图8为散斑原理图；

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

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

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

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

具体实施方式

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

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

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

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

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

$\begin{array}{c}\tilde{x}=x+u(x_0,y_0)+\frac{\∂u}{\∂x}{|}_{x_0,y_0}(x-x_0)+\frac{\∂u}{\∂y}{|}_{x_0,y_0}(y-y_0)\\ \tilde{y}=x+v(x_0,y_0)+\frac{\∂v}{\∂x}{|}_{x_0,y_0}(x-x_0)+\frac{\∂v}{\∂y}{|}_{x_0,y_0}(y-y_0)\end{array}---\left(1\right)$

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

令参数向量 $P=(u,v,\frac{\∂u}{\∂x},\frac{\∂v}{\∂x},\frac{\∂u}{\∂y},\frac{\∂v}{\∂y}),$ 并定义相关系数

$\mathrm{\ρ}=\frac{\underset{(x,y)\∈S}{\mathrm{\Σ}}{[f(x,y)-g(x,y,P)]}^2}{\underset{(x,y)\∈S}{\mathrm{\Σ}}{f}^2(x,y)}---\left(2\right)$

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

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

$\begin{array}{c}\▿\mathrm{\ρ}={\left(\frac{\∂\mathrm{\ρ}}{\∂P_i}\right)}_{i=1,\·\·\·,6}\\ =\frac{2}{\underset{(x,y)\∈S}{\mathrm{\Σ}}{f}^2(x,y)}{\left\{\underset{(x,y)\∈S}{\mathrm{\Σ}}\right[f(x,y)-g(x,y,P)\left]\frac{\∂g(x,y,P)}{\∂P}\right\}}_{i=1,\·\·\·,6}=0\end{array}---\left(3\right)$

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

$P=P_0-\frac{\▿\mathrm{\ρ}\left(P_0\right)}{\▿\▿\mathrm{\ρ}\left(P_0\right)}---\left(4\right)$

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

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

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

利用N-R迭代方法对位移场进行平滑处理：

求解出令相关系数ρ最小的形函数参数，根据得到的形函数参数测量出变形后的位移场数据dⁿ，其中位移场数据dⁿ由下式表示：dⁿ＝d+n；

其中，d为理想真实位移值，n为位移场噪声；

消除变形后的位移场数据U的随机误差，通过最小化下式消除随机误差

$D(d,\mathrm{\α})=\mathrm{\α}{\left|\right|d-d^n\left|\right|}_2^2+(1-\mathrm{\α}){\left|\right|\mathrm{Cd}\left|\right|}_2^2$

其中，dⁿ为变形后的位移场数据，d为理想真实位移值，n为位移场噪声，α为平滑参数；对上式求导，得到

$\frac{\∂D(d,\mathrm{\α})}{\∂d}=2\mathrm{\α}(d-d^n)+2(1-\mathrm{\α})*4\left(\mathrm{Cd}\right),$

$\frac{\∂D(d,\mathrm{\α})}{\∂\mathrm{\α}}={\left|\right|d-d^n\left|\right|}_2^2-{\left|\right|\mathrm{Cd}\left|\right|}_2^2$

其中，C为拉普拉斯算子，

将d和α变为列向量并令Q＝[d1,d2,d3...dm,α]，其中，m为d中元素的个数；则可得

$\▿D\left(Q\right)=\left[\begin{array}{c}2\mathrm{\α}(d_i-d^n)+2(1-\mathrm{\α})*4\left({\mathrm{Cd}}_i\right)\\ ......\\ 2\mathrm{\α}(d_m-d^n)+2(1-\mathrm{\α})*4\left(C{d}_m\right)\\ \mathrm{\α}\end{array}\right];$

根据▽D(Q)求得二阶偏导数▽▽D(Q)，

$\frac{\∂D^2\left(Q\right)}{\∂d_i\∂d_j}=\left[\begin{array}{cccc}32-32\mathrm{\α}& & & \\ & 32-30\mathrm{\α}& & \\ & & \·\·\·\·\·\·& \\ & & & 32-30\mathrm{\α}\end{array}\right]$

$\frac{\∂D^2\left(Q\right)}{\∂d_i\∂\mathrm{\α}}=2(d-d^n)-8\left(\mathrm{Cd}\right)$

$\frac{\∂D^2\left(Q\right)}{\∂\mathrm{\α}\∂d_j}=2(d-d^n)-8\left(\mathrm{Cd}\right)$

$\frac{\∂D^2\left(Q\right)}{\∂\mathrm{\α}\∂\mathrm{\α}}=0$

将上述矩阵组合，得到▽▽D(Q)；

利用得到的▽D(Q)和▽▽D(Q)，对下式进行迭代计算：

$P_{k+1}=P_k-\frac{\▿D\left(Q\right)}{\▿\▿D\left(Q\right)},$

其中，Q的初值为[dⁿ,α]^T，收敛条件为当上式收敛后得到的d即为理想真实位移值。

仿真分析：仿真散斑图参数，500pixel×500pixel，散斑颗粒半径4，散斑颗粒数4000

(1)均匀变形，x方向1000微应变，y方向位移0.3图像中添加噪声。分两种情况：

第一种情况：信噪比40db

表1位移、应变by N-R与本文计算结果误差比较

第一种情况：信噪比35db

表2位移、应变by N-R与本文计算结果误差比较

改变α的初值，进行迭代计算，α变化图2所示，图2为不同α初值条件下，α迭代变化情况

(2)非均匀变形：仿真散斑图参数，500pixel×500pixel，散斑颗粒半径4，散斑颗粒数4000

变形情况：u(x,y)＝0.1sin(2*pi*x/200)分两种情况：第一，信噪比30；第二，信噪比35

第一种情况：信噪比30db

请参阅图3、图4和图5所示，

表1、信噪比30db条件下，两种方法输出结果误差比较

方法	Std_v/pixel	Std_v y /με
			N-R迭代方法	0.02488	1227.8411
区分平滑位移t	0.02299	770.0738

第二种情况，信噪比35db

请参阅图6、图7和图8，

表2、信噪比35db条件下，两种方法输出结果误差比较

实验分析：铆钉板试件应变计算

实验时，铆钉板下端固定，上端拉伸，理论分析应该是铆钉下端受力最大。图9和图10给出的是铆钉周围区域纵向变形情况(可再补充其他方向变形情况)

(2)多孔试件实验

请参阅图11和图12所示，

结论：

经分析可知，要获得准确的应变场，需满足两个条件：

(1)N-R计算输出的位移场较准确。如果N-R输出的位移场偏差太大，那计算的应变结果必定不准确。位移场偏差很大，平滑后的数据峰值不会降很多，因为数据准确性和平滑性之间存在一个平衡；如果位移场数据偏小，那平滑后的值会更小，误差也大。因此，N-R较准确位移场的输出是必须的。

(2)合适的位移降噪方法。因为应变由位移差分得到，因此位移的波动会给应变计算带来较大误差。

Claims (4)

1.一种适用于数字图像相关的位移场迭代平滑方法，其特征在于：包括以下步骤：

1)、测量出变形后的位移场数据dⁿ，其中位移场数据dⁿ由下式表示：dⁿ＝d+n；其中，d为平滑后的位移场，n为位移场噪声；

2)、构建二次函数 $D(d,\mathrm{\α})=\mathrm{\α}{\left|\right|d-d^n\left|\right|}_2^2+(1-\mathrm{\α}){\left|\right|\mathrm{Cd}\left|\right|}_2^2$

其中，dⁿ为变形后的位移场数据，d为平滑后的位移场数据，n为位移场噪声，α为平滑参数，C为高阶算子；

求解二次函数的一阶偏导数，得到

$\frac{\∂D(d,\mathrm{\α})}{\∂d}=2\mathrm{\α}(d-d^n)+2(1-\mathrm{\α})*4\left(\mathrm{Cd}\right)$

$\frac{\∂D(d,\mathrm{\α})}{\∂\mathrm{\α}}={\left|\right|d-d^n\left|\right|}_2^2-{\left|\right|\mathrm{Cd}\left|\right|}_2^2$

其中，C为高阶算子，

将d和α变为列向量并令Q＝[d1,d2,d3...dm,α]，其中，m为点数；则可得

$\▿D\left(Q\right)=\left[\begin{array}{c}2\mathrm{\α}(d_1-{d_1}^n)+2(1-\mathrm{\α})*4\left({\mathrm{Cd}}_1\right)\\ 2\mathrm{\α}(d_2-{d_2}^n)+2(1-\mathrm{\α})*4\left({\mathrm{Cd}}_2\right)\\ ......\\ 2\mathrm{\α}(d_m-{d_m}^n)+2(1-\mathrm{\α})*4\left({\mathrm{Cd}}_m\right)\\ {\left|\right|d-d^n\left|\right|}_2^2-{\left|\right|\mathrm{Cd}\left|\right|}_2^2\end{array}\right];$

3)、根据步骤2)的▽D(Q)求解二次函数的二阶偏导数▽▽D(Q)，

$\frac{{\∂D}^2\left(Q\right)}{{\∂d}_i{\∂d}_j}=\left[\begin{array}{cccc}32-30\mathrm{\α}& & & \\ & 32-30\mathrm{\α}& & \\ & & ......& \\ & & & 32-30\mathrm{\α}\end{array}\right]$

$\frac{{\∂D}^2\left(Q\right)}{{\∂d}_i\∂\mathrm{\α}}=2(d-d^n)-8\left(\mathrm{Cd}\right)$

$\frac{{\∂D}^2}{{\∂\mathrm{\α}\∂d}_i}=2(d-d^n)-8\left(\mathrm{Cd}\right)$

$\frac{{\∂D}^2\left(Q\right)}{\∂\mathrm{\α}\∂\mathrm{\α}}=0$

将上述矩阵组合，得到▽▽D(Q)；

4)、利用步骤2)和3)得到的▽D(Q)和▽▽D(Q)，对下式进行迭代计算：

$P_{k+1}=P_k-\frac{\▿D\left(Q\right)}{\▿\▿D\left(Q\right)},$

其中，Q的初值为[dⁿ,α]^T，其中，dⁿ为变形后的位移场数据，当上式收敛后得到的d即为理想真实位移值。

2.如权利要求1所述的适用于数字图像相关的位移场迭代平滑方法，其特征在于：所述α的范围为(0,1)。

3.如权利要求1所述的适用于数字图像相关的位移场迭代平滑方法，其特征在于：步骤4)的收敛条件为其中，ε≤0.1。

4.如权利要求1所述的适用于数字图像相关的位移场迭代平滑方法，其特征在于：步骤4)中最大迭代次数N≤50。