CN110188759A - 一种在数字图像相关法中应变场子区动态选择方法 - Google Patents

Abstract

一种在数字图像相关法中应变场子区动态选择方法，提出基于位移场局部梯度强度的评价参数，用于描述局部非均匀变形的强度。通过权衡位移场噪声强度和局部非均匀变形强度，在非均匀变形梯度强度较小的地方，选择较大的应变计算子区，从而有效抑制位移场噪声；而在非均匀变形梯度强度较大的地方，选择较小的应变计算子区，从而保留局部变形特征。实现应变场计算子区大小的动态选择，使得有效抑制噪声的同时最大程度保留局部变形特征。本发明中基于位移场自身特征的子区大小动态选择过程与传统基于散斑灰度梯度信息的位移场计算子区大小动态选择过程不同，可以有效提高数字图像相关法在非均匀变形中应变场的测量精度。

Description

一种在数字图像相关法中应变场子区动态选择方法

技术领域

本发明属于数字图像处理领域，具体涉及一种在数字图像相关法中应变场子区动态选择方法。

背景技术

数字图像相关法(Digital Image Correlation，DIC)以数字图像处理技术为基础，是先进的一种现代光测力学方法，属于实验力学的前沿范畴。随着近年来工业相机制造技术和计算机处理性能的不断提升，DIC技术以其非接触式、自然光源、全场测量等优势特点，逐步在众多测试领域得到应用与实践。它除了应用在传统结构与材料力学性能测试外，在微纳电子系统、生物力学等领域具备了前所未有的解决问题能力。

影响测量精度的诸多因素被分为两类，一类是环境噪声，如相机噪声、环境光源波动、环境振动等；二类是算法误差，主要有位移场计算误差、应变场计算误差等。本发明提出一种应变场计算子区动态选择方法，用于优化应变计算过程，尤其在大位移梯度情况下效果更为明显。

传统DIC应变子区设置为固定值，即全图所有节点的应变计算采用相同的子区大小。这一设置在非均匀变形区域的计算是欠妥的。如裂纹近场应变计算过程中，裂纹近场的位移梯度极大，而随着与裂尖距离增加，位移梯度迅速下降。如果采用统一大小应变子区计算应变场，是十分不合理的。传统做法是将全局的子区设置为接近计算下限值，从而有效计算裂尖近场数据，然而这样又损失了距离裂纹尖端较远区域的应变计算精度。据此本发明提出一种基于位移梯度的数字图像相关应变场子区大小动态选择方法，用于改善DIC应变计算精度。

发明内容

本发明针对现有技术中的不足，提供一种在数字图像相关法中应变场子区动态选择方法，着重改善在非均匀变形测量中DIC应变计算精度。

为实现上述目的，本发明采用以下技术方案：

一种在数字图像相关法中应变场子区动态选择方法，其特征在于，包括：

在加载测试前首先采集若干静态图像，用于评价系统静态噪声与位移场跟踪精度，将系统静态的位移场跟踪精度等效为应变计算子区大小选择的上限值；

根据相关函数求解得到的亚像素精度位移场，提出局部变形梯度强度评价参数，用于描述局部非均匀变形的强度；

根据局部变形梯度强度不同程度地调低应变计算子区大小，实现动态选择过程，同时保留更为真实的局部变形特征。

为优化上述技术方案，采取的具体措施还包括：

进一步地，为权衡位移场噪声抑制效果与局部变形特征的描述能力，在非均匀变形梯度强度较小的地方，选择较大的应变计算子区，从而有效抑制位移场噪声；而在非均匀变形梯度强度较大的地方，选择较小的应变计算子区，从而保留局部变形特征。

进一步地，具体步骤如下：

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

2)在试验加载前，采集20张以上静态图像，用于评价系统静态噪声与位移场跟踪精度，从而确定应变计算子区大小选择的上限值；

3)在参考图像上划分感兴趣计算分析区域，称为ROI；

4)在ROI内划分用于计算位移场和应变场的网格节点，网格节点间距为5～10像素；

5)在参考图像中，以步骤3)的ROI内按照步骤4)所划分的任意网格节点为中心，分别建立边长为Subset size的相关子区，称为参考图像子区；

6)在目标图像中，以步骤5)的参考图像子区为相关函数匹配范围，求解对应节点的目标图像子区；

7)定义形函数，用于描述加载过程中从参考图像子区到目标图像子区的变形情况；

8)定义灰度插值函数，用于构造对应参考图像子区灰度矩阵和形函数参数的目标图像子区灰度矩阵；

9)定义相关函数，结合步骤6)-8)所定义的形函数与灰度插值函数，对步骤3)的ROI内按照步骤4)所划分的网格节点，迭代求解亚像素精度位移；

10)按照所提相关函数，逐一求解步骤4)所划分的全部网格节点的位移，构成步骤3)的ROI内全部网格节点位移场；

11)基于局部最小二乘原理的应变场求解方法，构建从步骤10)所得位移场到应变场的计算过程；

12)在参考图像中，以任意网格节点为中心，建立边长为Filter size的应变计算子区，用于步骤11)的计算过程；

13)定义局部变形梯度强度评价参数，用于描述局部非均匀变形的强度，并求解在不同Filter size下对应函数值，Filter size的取值范围为5～21；

14)步骤12)中Filter size的上限由步骤2)确定，逐步向下调整，在步骤13)所得的函数表中查找最适值；

15)将步骤14)中所得的Filter size值代入步骤11)，逐一求解网格节点应变向量，构建集合得到应变场。

进一步地，步骤7)中，形函数定义如下：

一阶形函数：

二阶形函数：

其中，u(x，y)、v(x，y)分别表示平面内坐标(x，y)处的水平和竖直方向位移，a₀₀、a₁₀、a₀₁、a₂₀、a₁₁、a₀₂、b₀₀、b₁₀、b₀₁、b₂₀、b₁₁、b₀₂为形函数各项系数；形函数是描述子区变形特征的函数，由形函数系数确定其具体变形情况，若已知形函数系数，则输入坐标(x，y)，即可得对应点的位移u(x，y)、v(x，y)；反之，若已知一系列坐标与对应点位移关系，则可通过最小二乘方法拟合求解形函数系数；该形函数定义不仅用于位移场计算过程，而且用于应变场计算过程。

进一步地，步骤8)中，插值方法为双线性插值法或双三次样条插值法。

进一步地，步骤9)中，相关函数为归一化的标准互相关函数

或者归一化的最小平方距离相关函数

其中表示相关函数，f(x，y)是参考图像子区中坐标为(x，y)点处灰度值，g(x′，y′)是目标图像子区中对应点坐标为(x′，y′)点处灰度值；M表示位移场计算子区内的以像素为单位的局部坐标，其大小由位移场计算子区半径确定；

和分别表示当前计算参考图像子区和目标图像子区的灰度均值。

进一步地，步骤11)中，DIC应变场计算采用基于子区局部最小二乘方法，通过调整全局子区大小来进行平滑效果的控制；

一阶形函数时，计算公式如下：

XA＝U

上式对应权利要求4中形函数公式的简化表达，式中，X表示应变计算子区内节点的坐标矩阵，A表示形函数系数向量，U表示应变计算子区内节点的位移向量；

其展开形式为：

其中，m表示应变计算子区内的以节点为单位的局部坐标，其大小由应变计算子区半径确定；X矩阵为(2m+1)²×3的矩阵，表示子区内有(2m+1)²个节点，第二列表示对应节点的水平方向坐标，第三列表示对应节点的竖直方向坐标；

通过伪逆求解得A＝(X^TX)^-1X^TU，同理得到B＝(X^TX)^-1X^TV，由此得到形函数系数A、B，V表示应变计算子区内节点的位移向量；在小变形情况下的柯西应变分量按如下计算：

ε_x、ε_y、γ_xy这三个变量共同组成二维平面应变张量，用来表示平面内坐标(x，y)处的应变状态，其中，ε_x表示x方向的正应变，ε_y表示y方向的正应变，γ_xy表示切应变；

有理论推导证明，计算应变场时，采用一阶形函数与二阶形函数描述下的应变计算输出是完全相同的。

进一步地，步骤13)中，局部变形梯度强度评价参数是对位移梯度的量化以及用于应变计算子区大小选择的重要参考，计算过程如下：

式中，r表示应变计算子区的半径，即Filter size＝(2r+1)；dis_x、dis_y分别表示x方向位移和y方向位移，Ω表示应变计算子区内全体网格节点组成的集合，Ω_r表示半径r的应变计算子区内全体网格节点组成的集合，λ表示局部变形梯度强度评价参数，下标x、y分别表示x方向和y方向。

本发明的有益效果是：主要针对现有DIC技术在应变场计算算法中的不足，考虑到非均匀变形计算不宜采用全局统一应变子区大小的实际情况，提出一种基于位移梯度的应变子区大小动态选择方法，较传统DIC计算结果在非均匀变形应变场计算中有更好的精度和准确性。

附图说明

图1为应变场子区大小的动态选择示意图。

图2为局部位移场误差与非均匀性关系示意图。

图3为带状非均匀变形模拟散斑图。

图4为位移幅值α＝1/π动态子区选择应变场计算仿真图。

图5为位移幅值α＝1动态子区选择应变场计算仿真图。

具体实施方式

现在结合附图对本发明作进一步详细的说明。

图1为应变场子区大小的动态选择示意图，由图1可知，在较大变形梯度区域应选择较小的应变计算子区r_L，在较小变形梯度区域应选择较大的应变计算子区r_S。

图2为局部位移场误差与非均匀性关系示意图，水平面坐标分别表示网格节点尺度，垂向坐标表示对应节点的位移量，图中可以看出：(1)在子区中心的位移相比子区上下边缘的位移要小；(2)子区边缘位移非平滑，而是有锯齿状的波动，表示位移场测量误差；表明随着子区尺寸的增大，子区内非均匀变形成分超过位移场误差量级，此时再用大子区尺寸进行平滑降噪求解应变场的方法不再适宜，需要选择较小的子区从而协调子区内非均匀变形的影响。其量化指标为局部变形梯度强度评价参数，是子区尺寸的函数。随着子区尺寸的增加，在均匀变形区域，函数值增长十分缓慢；在非均匀变形区域，函数值增长急剧。结合由环境本底噪声等效的最大子区及子区上限值，可以逐步确定任意节点出最适子区尺寸，从而实现应变计算子区动态选择的过程。

图3为带状非均匀变形模拟散斑图，基于反向映射法的模拟散斑生成方法生成变形后仿真散斑图，在水平方向位移场模拟加入按正弦规律变化的非均匀变形，带状模拟散斑图可以呈现其多个周期性规律。在水平方向模拟加入的位移场由下式给出：

U(X)＝αsin(2πτX/q)

其中α表示水平位移峰值，q表示正弦规律变化位移场的周期；当周期q固定的情况下，幅值α反映非均匀变形的强度。设定周期q＝200像素，设定幅值α＝1/π和α＝1，分两组进行仿真验证分析。对应的，ε_xmax＝2πα/q＝1％表示较小程度的非均匀变形，而ε_xmax＝2πα/q≈3.14％表示较大程度的非均匀变形。如图4所示，位移幅值α＝1/π情况下x方向位移场误差的仿真结果中，虚线表示未经算法修正的原始误差，实线表示经过算法修正后的误差量。如图5所示表示位移幅值α＝1情况下的仿真结果。

由上可知，本发明基于位移场自身特征进行子区大小动态选择，可以有效提高数字图像相关法在非均匀变形中应变场的测量精度。在工程实践中，非均匀变形测量问题更加广泛，例如裂纹近场变形测量等，因而该方法具有普遍意义。

需要注意的是，发明中所引用的如“上”、“下”、“左”、“右”、“前”、“后”等的用语，亦仅为便于叙述的明了，而非用以限定本发明可实施的范围，其相对关系的改变或调整，在无实质变更技术内容下，当亦视为本发明可实施的范畴。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims (8)

1.一种在数字图像相关法中应变场子区动态选择方法，其特征在于，包括：

在加载测试前首先采集若干静态图像，用于评价系统静态噪声与位移场跟踪精度，将系统静态的位移场跟踪精度等效为应变计算子区大小选择的上限值；

根据相关函数求解得到的亚像素精度位移场，提出局部变形梯度强度评价参数，用于描述局部非均匀变形的强度；

根据局部变形梯度强度不同程度地调低应变计算子区大小，实现动态选择过程。

2.如权利要求1所述的一种在数字图像相关法中应变场子区动态选择方法，其特征在于：在非均匀变形梯度强度较小的地方，选择较大的应变计算子区；而在非均匀变形梯度强度较大的地方，选择较小的应变计算子区。

3.如权利要求1所述的一种在数字图像相关法中应变场子区动态选择方法，其特征在于：具体步骤如下：

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

2)在试验加载前，采集20张以上静态图像，用于评价系统静态噪声与位移场跟踪精度，从而确定应变计算子区大小选择的上限值；

3)在参考图像上划分感兴趣计算分析区域，称为ROI；

4)在ROI内划分用于计算位移场和应变场的网格节点，网格节点间距为5～10像素；

5)在参考图像中，以步骤3)的ROI内按照步骤4)所划分的任意网格节点为中心，分别建立边长为Subset size的相关子区，称为参考图像子区；

6)在目标图像中，以步骤5)的参考图像子区为相关函数匹配范围，求解对应节点的目标图像子区；

7)定义形函数，用于描述加载过程中从参考图像子区到目标图像子区的变形情况；

8)定义灰度插值函数，用于构造对应参考图像子区灰度矩阵和形函数参数的目标图像子区灰度矩阵；

9)定义相关函数，结合步骤6)-8)所定义的形函数与灰度插值函数，对步骤3)的ROI内按照步骤4)所划分的网格节点，迭代求解亚像素精度位移；

10)按照所提相关函数，逐一求解步骤4)所划分的全部网格节点的位移，构成步骤3)的ROI内全部网格节点位移场；

11)基于局部最小二乘原理的应变场求解方法，构建从步骤10)所得位移场到应变场的计算过程；

12)在参考图像中，以任意网格节点为中心，建立边长为Filter size的应变计算子区，用于步骤11)的计算过程；

13)定义局部变形梯度强度评价参数，用于描述局部非均匀变形的强度，并求解在不同FiLter size下对应函数值，Filter size的取值范围为5～21；

14)步骤12)中Filter size的上限由步骤2)确定，逐步向下调整，在步骤13)所得的函数表中查找最适值；

15)将步骤14)中所得的Filter size值代入步骤11)，逐一求解网格节点应变向量，构建集合得到应变场。

4.如权利要求3所述的一种在数字图像相关法中应变场子区动态选择方法，其特征在于：步骤7)中，形函数定义如下：

一阶形函数：

二阶形函数：

其中，u(x，y)、v(x，y)分别表示平面内坐标(x，y)处的水平和竖直方向位移，a₀₀、a₁₀、a₀₁、a₂₀、a₁₁、a₀₂、b₀₀、b₁₀、b₀₁、b₂₀、b₁₁、b₀₂为形函数各项系数；形函数是描述子区变形特征的函数，由形函数系数确定其具体变形情况，若已知形函数系数，则输入坐标(x，y)，即可得对应点的位移u(x，y)、v(x，y)；反之，若已知一系列坐标与对应点位移关系，则可通过最小二乘方法拟合求解形函数系数。

5.如权利要求3所述的一种在数字图像相关法中应变场子区动态选择方法，其特征在于：步骤8)中，插值方法为双线性插值法或双三次样条插值法。

6.如权利要求3所述的一种在数字图像相关法中应变场子区动态选择方法，其特征在于：步骤9)中，相关函数为归一化的标准互相关函数

或者归一化的最小平方距离相关函数

其中表示相关函数，f(x，y)是参考图像子区中坐标为(x，y)点处灰度值，g(x′，y′)是目标图像子区中对应点坐标为(x′，y′)点处灰度值；M表示位移场计算子区内的以像素为单位的局部坐标，其大小由位移场计算子区半径确定；

和分别表示当前计算参考图像子区和目标图像子区的灰度均值。

7.如权利要求4所述的一种在数字图像相关法中应变场子区动态选择方法，其特征在于：步骤11)中，DIC应变场计算采用基于子区局部最小二乘方法，通过调整全局子区大小来进行平滑效果的控制；

一阶形函数时，计算公式如下：

XA＝U

式中，X表示应变计算子区内节点的坐标矩阵，A表示形函数系数向量，U表示应变计算子区内节点的位移向量；

其展开形式为：

其中，m表示应变计算子区内的以节点为单位的局部坐标，其大小由应变计算子区半径确定；X矩阵为(2m+1)²×3的矩阵，表示子区内有(2m+1)²个节点，第二列表示对应节点的水平方向坐标，第三列表示对应节点的竖直方向坐标；

通过伪逆求解得A＝(X^TX)^-1X^TU，同理得到B＝(X^TX)^-1X^Tv，由此得到形函数系数A、B，V表示应变计算子区内节点的位移向量；在小变形情况下的柯西应变分量按如下计算：

ε_x、ε_y、γ_xy这三个变量共同组成二维平面应变张量，用来表示平面内坐标(x，y)处的应变状态，其中，ε_x表示x方向的正应变，ε_y表示y方向的正应变，γ_xy表示切应变。

8.如权利要求3所述的一种在数字图像相关法中应变场子区动态选择方法，其特征在于：步骤13)中，局部变形梯度强度评价参数是对位移梯度的量化以及用于应变计算子区大小选择的重要参考，计算过程如下：

式中，r表示应变计算子区的半径，即Filter size＝(2r+1)；dis_x、dis_y分别表示x方向位移和y方向位移，Ω表示应变计算子区内全体网格节点组成的集合，Ω_r表示半径r的应变计算子区内全体网格节点组成的集合，λ表示局部变形梯度强度评价参数，下标x、y分别表示x方向和y方向。