3种亚像素位移测量算法的比较研究

资新运， 耿 帅， 赵姝帆， 夏军剑， 舒 放

（1.军事交通学院工程实验中心，天津300161； 2.天津财经大学理工学院，天津300222）

3种亚像素位移测量算法的比较研究

资新运1， 耿 帅1， 赵姝帆1， 夏军剑1， 舒 放2

（1.军事交通学院工程实验中心，天津300161； 2.天津财经大学理工学院，天津300222）

梯度法、曲面拟合法和N-R法是数字散斑相关方法中提高位移测量精度的3种主要亚像素位移算法。研究了3种算法的原理并进行了数值模拟仿真。首先使用计算机模拟生成一系列位移为0.01 pixel的散斑图，并在生成的图像中添加方差不同的高斯噪声，对3种算法在噪声条件下的精度进行比较，结果显示在精度要求为0.01 pixel时，梯度法的噪声方差上限为0.006，N-R法的噪声方差上限为0.000 5，而曲面拟合法的计算精度具有较大不确定性，无法保证精度，梯度法更适用于工程应用。

计量学；数字散斑相关；曲面拟合；梯度法；N-R法；亚像素

1 引 言

数字散斑相关方法（Digital Speckle Correlation Method，DSCM）又称数字图像相关方法（Digital Image Correlation Method，DICM），是上世纪80年代提出的一种可以测量物体变形和应变的光学非接触测量方法［1］，它通过分析变形前后物体表面散斑图像的灰度来获得被测物体的力学性能。该方法具有非接触、全场测量、实验设备简单，对环境要求低等优点［2］，经过多年的研究，其在理论和应用方面取得了较大的发展。近年来，基于数字图像相关方法的测量软件已经面世，应用领域包括材料力学、振动力学、冲击力学、破坏力学、微纳米力学、疲劳及可靠性等［3～5］，并可实时监测零部件应力分布情况［6，7］，预防故障发生。

数字散斑相关技术通过记录物体变形前后的图像并进行相关计算得到物体的位移和变形。由于该方法处理的是数字化的图像，因此获得的位移只能是像素的整数倍，但实际的位移值通常不恰好为整像素，而且整像素的位移精度在固体材料或构件表面的变形测量中是远远不够的。为获得更高测量精度，在不提高光学系统放大倍数和CCD分辨率的前提下，需要对数字散斑进行亚像素处理。目前，相关研究人员提出了亚像素灰度插值法、曲面拟合法、梯度法、N-R法（Newton-rapshon method）［8～10］等亚像素搜索算法，其中应用最为广泛的是曲面拟合法、梯度法和N-R法。在无噪声条件下3种算法的测量精度已有定论［11］，然而在工程测量中，图像采集和传输过程中，电器系统和外界影响都可产生噪声。由系统外部环境干扰引起的噪声可以通过改善外部环境来消除或减小，由光、电、元器件本身及系统内部设备电路引起的噪声，有些内部噪声是可以减小的，但有些是不可避免的。噪声会使散斑图像表面产生不同的灰度特征，从而对数字图像相关计算结果产生一定影响，因此有必要对噪声情况下的3种亚像素位移算法进行比较研究。

2 数字散斑相关方法的基本原理

在外力作用下，物体表面发生应变位移，散斑的灰度值保持不变，因此可以计算散斑位置的变化，从而求得物体表面的位移。利用数字图像相关方法时，首先采集试件表面变形前后的两幅散斑图，在变形前的图像中，取以待求点P（x1，y1）为中心的一块矩形图像子区f（x，y），这块区域的灰度特征保持不变，因此，根据灰度信息在变形后的目标图像中寻找相同的矩形区域g（x′，y′）和它的中心点P（x2，y2），计算得到两中心点的位移u，v，如图1所示。在寻找过程中是通过一定的搜索方法，并按某一相关函数来进行相关计算，求得相关系数C（u，v）取最大值的矩形区域的［4］。

式中：x′＝x＋u，y′＝y＋v，u和v分别为中心点在x方向和y方向上的位移；Π是刻画f（x，y）和g（x′，y′）在某种程度上相似的函数。目前相关函数有近10种，可以根据要求选择合适的相关函数。

求相关函数最值的过程就是求子区中心点位移和变形的过程，选择不同的相关函数，两区域完全相关时相关系数取最大或最小。根据峰值相关系数计算该子区域在变形后的位置，由此得到该区块的变形量。对变形前散斑中的所有子区域进行类似运算，就可以得到整个位移场。

图1 数字图像相关原理示意图

3 数字散斑相关中的亚像素搜索算法

3.1 曲面拟合法

整像素运算后得到整像素位移结果（u，v），同时也得到了参考窗口与变形后图像中目标窗口间的相关系数值，C（u，v）为全部相关系数值中的最大值。曲面拟合法是在整像素目标点及其邻近点的相关系数阵的基础上，对其进行曲面拟合，拟合曲面的极值点位置即为所求图像区域的中心位置。常用的曲面拟合法有高斯曲面拟合、二次曲面拟合、二维拉格朗日曲面拟合等方法。本文选取二次曲面拟合法来实现相关系数场的模拟，再通过求导的方法获得模拟函数的极大值，得到亚像素位移。

以整像素搜索算法得到的相关系数峰值点为中心，取3×3 pixels拟合窗口，通过这9点相关系数值及坐标位置进行曲面拟合。二次曲面的拟合公式一般采用二元二次多项式：

式中有6个待定系数，A＝［a0，a1，a2，a3，a4，a5］T。以相关系数峰值点为中心，取3×3 pixels拟合窗口进行拟合，可得到包含9个式（2）的方程组，因此可利用最小二乘法来求解。

对于拟合函数C，拟合曲面极值点满足对x和y的偏导数为0，即：

由式（3）可知曲面极值点位置为：

3.2 梯度法

基于梯度的亚像素位移算法是假设图像子区做微小位移时，若子区足够小，则图像子区运动可近似为面内刚体运动［11］，此时有：

式中：f（x，y）和g（x′，y′）分别表示位移前、位移后的子区图像；u和v分别为与x、y方向上的整像素位移；Δu，Δv为相应的亚像素位移。

式（5）可通过一阶的泰勒级数展开，同时舍弃高阶小量，可得：

gx和gy为子区灰度一阶梯度，这里选用抗噪声能力强的标准化协方差互相关函数的平方：

3.3 N-R法

N-R法考虑子区变形的情况，假设变形前子区灰度分布为f（x，y），变形后相应分布为g（x′，y′），则存在以下对应关系：

式中：u和v为图像子区中心在x和y方向上的位移；Δx和Δy为点（x，y）到子区中心（x0，y0）的距离；∂u／∂x，∂u／∂y，∂v／∂x，∂v／∂y为位移一阶偏导数。

为了提高精度，选择零均值归一化的最小平方距离相关函数：

两子区最相似即为相关系数C取最小值，令P＝（u，v，∂u／∂x，∂u／∂y，∂v／∂x，∂v／∂y）T，即公式（21）可以表示为C＝C（P），参数可以通过对该式求极值获得：

6个未知参数构成一个非线性方程，取迭代的估计初值为：

其中u0，v0为整像素搜索得到的x和y方向上的位移。

通过N-R法求解：

其中▽▽C（P）为相关系数的二次偏导数，也被称为Hessian矩阵。Hessian矩阵可近似为［12］：

因此有：

经推导，可得迭代关系为：

其分量形式为：

4 亚像素位移测量算法的数值模拟验证

4.1 模拟散斑生成方法

为了避免在实际图像采集过程中引入误差，更好地验证算法，可采用计算机模拟生成数字散斑，进行仿真实验验证。文献［13］对模拟散斑的颗粒大小、颗粒数目等做了研究，并得到以下结论：最佳散斑颗粒尺寸在3～6 pixels，最佳计算窗口为31×31～51×51 pixels。借助模拟数字散斑，通过精确的控制位移量、自由添加噪声，可以仿真不同位移、噪声环境下的变形实验。

根据文献［13］提供的算法，生成一对位移为（u，v）的模拟散斑，位移分量的表达式为：

式中：u0，v0为刚体位移量；ux，uy，vx，vy为位移梯度。

假设物体变形前后的数字散斑呈现高斯分布，以高斯函数模拟散斑颗粒光强分布，则变形前后图像灰度分布函数分别为I1（x，y）和I2（x，y）。

式中：s为散斑颗粒数量；d为散斑颗粒尺寸；I0为散斑图像背景光强，计算中取值为1；（xk，yk）为散斑颗粒的位置，在程序中可由随机分布函数生成。

根据上述方法，可生成模拟数字散斑图如图2所示，其中散斑图像大小为256×256 pixels，散斑颗粒尺寸为3 pixels，颗粒数为1 000个。

图2 数字散斑图像及其灰度分布图

图3 理想情况的数字散斑图平移0～0.1 pixel时3种亚像素算法测量误差

4.2 理想条件下亚像素位移模拟验证

为了验证以上3种亚像素测量算法的性能，用计算机模拟生成一幅没有噪声的数字散斑图像作为基准，然后依次平移0.01 pixel，连续生成10幅数字散斑序列图作为目标图像。计算点阵为11×11共121个点，计算子区大小均采用41×41。使用本文的3种亚像素处理算法计算位移，得到位移的均值误差和标准差如图3所示。

由图3可见，在理想情况下，N-R法计算精度最高，梯度法计算精度稍劣于N-R方法，曲面拟合法最差。这与文献［11］中得到的结论相符。

4.3 噪声影响下亚像素位移模拟验证

在实际工程应用中，测量过程难免会受到噪声的影响，根据文献［14］，噪声主要的来源有：CCD暗电流影响、CCD复位噪声、图像卡噪声、外部环境振动等。为模拟更加真实的情况，讨论噪声对3种亚像素算法的影响，分别为基准图像添加均值为0，方差σ为0.000 1、0.001、0.01、0.1的高斯噪声，添加噪声后的数字散斑图像如图4所示。

以相同的方法，分别以每幅含噪声的数字散斑图为新的基准图像，依次平移0.01 pixel，连续生成10幅数字散斑序列图作为目标图像。计算其均值误差与标准差结果如图5所示。图5（a）、（b）分别为σ＝0.000 1时3种算法均值误差和标准差，图5（c）、（d）分别为σ＝0.001时3种算法均值误差和标准差，图5（e）、（f）分别为σ＝0.01时3种算法均值误差和标准差。

由图5（a）中可以看到，当σ＝0.000 1时梯度法和N-R法的测量精度均能达到0.01 pixel，曲面拟合法已达不到精度要求；由图5（c）中可以看到，当σ＝0.001时梯度法的测量精度能达到0.01 pixel，N-R法也达不到精度要求；由图5（e）中可以看到，当σ＝0.01时梯度法也达不到精度要求。3种亚像素算法的均值误差和标准差均呈随着噪声增大而增加的趋势，其中N-R对噪声最敏感，而梯度法受噪声影响最小，曲面拟合法介于两者之间。当噪声方差在0.01时，3种算法计算误差均大于了0.01 pixel。对于精度为0.01 pixel的测量而言，3种算法都无法实现，继续增加噪声方差已无意义，本文不再计算在噪声方差为0.1时3种算法的测量精度。

从图5（a）和（c）中可以大致判断，当测量精度要求为0.01时，N-R法对应的噪声方差上限σ应在0.000 1～0.001之间。为了确定N-R法的噪声方差上限值，现将一张模拟散斑图依次加入均值为0，方差大小为0.000 1～0.001的高斯噪声，每次增量为0.000 1，得到10加入噪声后的图像作为目标图像，原始图像为参考图像。按照同样的方法计算均值误差和标准差。

由图6（a）中可以看到，N-R法在x方向上要满足0.01 pixel的测量精度，噪声方差应控制在0.000 6以下，由图6（c）中可以看到，在y方向上要满足0.01 pixel的测量精度，噪声方差应控制在0.000 5以下。因此x、y方向要同时满足测量精度，噪声方差应控制在0.000 5以下。

从图5（c）和（e）中可以大致判断，梯度法对应的噪声方差上限σ应在0.001～0.01之间，为了确定梯度法的噪声方差上限值，运用同样的方法在一张模拟散斑图依次加入均值为0，方差大小为0.001～0.01的高斯噪声，每次增量为0.001。

由图7（a）中可以看到，梯度法在x方向上要满足0.01 pixel的测量精度，噪声方差应控制在0.007以下，由图7（c）中可以看到，在y方向上要满足0.01 pixel的测量精度，噪声方差应控制在0.006以下。因此x、y方向要同时满足测量精度，噪声方差应控制在0.006以下。

图5 噪声影响的数字散斑图平移0～0.1 pixel时3种亚像素算法测量误差

图6 噪声对N-R法计算精度的影响

图7 噪声对梯度法计算精度的影响

5 结 论

在无噪声条件下，N-R法计算精度最高，梯度法次之，曲面拟合法最差。在精度要求为0.01 pixel时，梯度法的噪声方差上限为0.006，N-R法的噪声方差上限为0.000 5，而曲面拟合法的计算精度具有较大不确定性，无法保证精度。3种方法的均值偏差都随噪声增大而增大，其中梯度法的均值偏差最小，曲面拟合法次之，N-R的均值误差最大；3种方法的标准差随着噪声的增大而增加，N-R法的标准差增加的幅度最大。由此可见，梯度法的抗噪能力最强，曲面拟合法次之，N-R法抗噪能力最弱。由此可见，N-R法在无噪声情况下精度很高，但是计算速度慢并且对噪声非常敏感，所以N-R法适用于实验条件下采集的图像质量较高的情况；曲面拟合法的计算速度快但精度低，因此适用于计算速度要求高、精度要求一般的情况；梯度法的抗噪声能力较好，精度较高，计算速度在两者之间，可广泛运用于工程实践。

］

［1］ PetersW H，Ranson W F.Digital imaging techniques in experimental stress analysis［J］.Optical Engineering，1982，21（3）：427-431.

［2］ 潘兵，谢惠民，续伯钦，等.数字图像相关中的亚像素位移定位算法进展［J］.力学进展，2005，35（3）：345-352.

［3］ 宋义敏，马少鹏，杨小彬，等.岩石变形破坏的数字散斑相关方法研究［J］.岩石力学与工程学报，2011，30（1）：170-175.

［4］ Dong Y，Kakisawa H，Kagawa Y.Optical system for microscopic observation and strain measurement at high temperature［J］.Measurement Science and Technology，2014，25（2）：025002.

［5］ 俞海，郭荣鑫，夏海廷，等.数字图像相关法在WC／Cu复合材料线膨胀系数测量中的应用［J］.光学精密工程，2013，21（10）：2696-2703.

［6］ Cofaru C，PhilipsW，Van Paepegem W.A novel speckle pattern—adaptive digital image correlation approach with robust strain calculation［J］.Optics and Lasers in Engineering，2012，50（2）：187-198.

［7］ Sicard J，Sirohi J.Measurement of the deformation of an extremely flexible rotor blade using digital image correlation［J］.Measurement Science and Technology，2013，24（6）：065203.

［8］ 潘兵，续伯钦，陈丁，等.数字图像相关中亚像素位移测量的曲面拟合法［J］.计量学报，2005，26（2）：128-134.

［9］ 潘兵，吴大方，谢惠民，等.基于梯度的数字体图像相关方法测量物体内部变形［J］光学学报，2011，31（6）：0612005

［10］ Cofaru C，Philips W，Van Paepegem W.Improved Newton-Raphson digital image correlation method for full-field displacement and strain calculation［J］.Applied Optics，2010，49（33）：6472-6484.

［11］ 潘兵，谢惠民，戴福隆.数字图像相关中亚像素位移测量算法的研究［J］.力学学报，2007，39（2）：245-252.

［12］ Vendroux G，KnaussW G.Submicron deformation field measurements：Part 2.Improved digital image correlation［J］.ExperimentalMechanics.1998，38（2）：86-92.

［13］ Zhou P，Goodson K E.Sub-pixel displacement and deformation gradient measurement using digital image／speckle correlation［J］.Optical Engineering，2001，40（8）：1613-1620.

［14］ 杨勇，王琰蕾，李明.高精度数字图像相关测量系统及其技术研究［J］.光学学报，2006，26（2）：197-201.

Com parison of Three Sub-pixel Displacement Algorithm

ZIXin-yun1， GENG Shuai1， ZHAO Shu-fan1， XIA Jun-jian1， SHU Fang2
（1.Engineering Experiment Center，Military Transportation University，Tianjin 300161，China；2.Institute of Technology，Tianjin University of Finance，Tianjin 300222，China）

Gradientmethod，surface fitting and N-R method are 3 main sub-pixel displacement algorithms of digital speckle correlationmethod to improve the precision of themeasuring system.The principles of three kinds of algorithm were studied and numerical simulation was carried out.To compare the different precision from 3 methods under noise，the computermodels generate a series of speckle gram whose displacement is 0.01 pixels and different Gaussian noise was added in generated graphics.The results showed that in the precision requirement in 0.01 pixels，the upper bounds on noise variance of gradientmethod and N-R method are respectively 0.006 and 0.000 5，and the calculation precision of surface fittingmethod has huge uncertainty，gradientmethod ismore suitable for engineering application.

Metrology；Digital speckle correlation；Quadratic surface fitting；Gradientmethod；N-R method；Subpixel

TB93

：A

：1000-1158（2015）03-0260-08

10.3969／j.issn.1000-1158.2015.03.09

2014-08-20；

：2014-11-27

国家863计划项目（2013AA065303）；国家自然科学基金（91120306）；天津市自然科学基金（14JCQNJC01600）

资新运（1971-），男，湖南衡阳人，军事交通学院教授，博士生导师，主要研究方向为传感与检测技术、数字图像处理。lookttlook＠163.com