激光散斑的亚像素位移法计算及比较

王希军

(中国科学院长春光学精密机械与物理研究所，吉林长春130033)

1 引言

数字散斑相关方法(Digital Speckle Correlation Method，DSCM)是一种有效的光学测量手段。与传统的干涉测量方法不同，它直接从物体表面随机分布的人工或自然散斑场中提取位移和变形信息，不需要复杂的干涉条纹处理过程，具有全场测量、非接触等优点，且其光路简单，对测量环境要求不高，因而被广泛应用于实验力学及其他科学研究领域［1-4］。

在数字散斑相关测量中，获取的数字散斑图是离散的灰度信息，经相关运算后，得到整像素为单位的位移值，而物体的真实位移值不一定恰好是整像素，因此，亚像素相关测量方法已成为数字散斑相关测量中的关键技术和研究重点。

目前，国内外常用的亚像素相关测量方法有插值法［5-6］、曲面拟合法［7］、牛顿-拉夫森法［8］、梯度算法［9-10］、频域相关算法［10-11］以及遗传算法［12］等。虽然对这些计算算法的精度分析研究开展得较早，但很少有利用算法分析和计算激光散斑图并比较各种算法精度的工作，因此，根据不同的测量要求选择合适的亚像素算法对于提高激光散斑相关测量技术的精度和计算速度等工作显得尤为重要。

本文逐一比较了插值算法、梯度算法和牛顿-拉夫森迭代算法的精度和效率。分析了处理激光散斑图时这些算法的特点及适用性，为满足激光散斑相关测量数据处理和计算精度的特定要求，提供选择亚像素计算方法的依据。

2 亚像素计算散斑位移图

散斑场中各点周围的邻近区域通常称为子集或窗口。选取散斑场中某一点为中心的窗口作为该点位移和形变的信息，分析和搜索该窗口的移动和变化，可获得该点的位移和形变信息。在进行数字散斑相关测量过程中，首先在物体变形前后拍摄两幅散斑图，如图1所示。在变形前的散斑图(参考图)中以待测点P为中心取m×m的子区A，当被测物体发生形变后，即在物体变形后所拍的散斑图中(目标图)，子区A移至子区B的位置，相应P点移动到P'。由统计学可知，A与B两个样本空间的相关系数最大。因此可以根据相关系数的最大值来确定子区B的位置，从而得到P 点的位移［(x'－x)，(y'－y)］。

图1 数字散斑相关方法测量原理示意图Fig.1 Schematic diagram of digital speckle correlation method

但通过上述方法只能得到整像素位移，而物体实际的位移值往往不是整像素，并且整像素的位移定位精度是远远不够的，需要进行亚像素相关测量。常用的亚像素相关方法有插值法、牛顿-拉夫森法及梯度法。

2.1 相关系数插值法

该方法先用十字搜索法或其它整像素搜索方法算出子区的整像素位移，然后对整像素位移点相邻9×9点的相关系数矩阵进行二维插值。插值函数可以为双线性插值、双三次样条插值函数等，一般双三次样条插值和双线性插值具有较高的精度［13］。相关系数插值法的算法比较容易实现，但计算量较大。本文主要讨论采用双三次样条插值计算激光散斑图的结果。

2.2 梯度算法

基于梯度的亚像素位移算法最初由Davis和Freeman［9］提出，该方法假设物体只发生刚体运动。当物体发生微小变形且选择的计算子区足够小时，则该子区即可近似刚体运动。令f(x，y)，g(x*，y*)分别表示物体变形前后子区的灰度，则:

式中:u，v分别为子区在x，y方向上的整像素位移，Δx，Δy分别为子区在x，y方向上的亚像素位移。选取最小平方距离相关函数作为子区匹配的算子，则:

图2 邻域相关系数矩阵差值示意图Fig.2 Schematic diagram of neighborhood correlation coefficient matrix interpolation

真实的亚像素位移Δx，Δy应使式(2)取得最小值，因此由极值条件:

式中，gx，gy为灰度的一阶导数。

2.3 牛顿-拉夫森迭代算法

牛顿-拉夫森迭代法最初由 Bruck［8］提出，由于在算法中引入了位移的一次导数，因此其可以对具有形变的图像子区进行计算。

令f(x，y)为样品变形前拍摄的散斑图的灰度分布，g(x'，y')为变形后对拍摄的散斑图的灰度分布。对于变形前散斑图的某一子区灰度分布fs(x，y)，设其在变形后对应于灰度分布 gs(x'，y')，则应有 gs(x'，y')=fs(x，y)，且应有以下对应关系:

式中:u为变形后子区中心(x0，y0)在x方向上的位移，v为变形后子区中心(x0，y0)在y方向上的位移，Δx，Δy为点(x，y)到子区中心(x0，y0)的距离，∂u/∂x，∂u/∂y，∂v/∂x，∂v/∂y 为位移的一阶偏导数。选取最小平方距离相关函数:

由式(6)可见，子区 f(x，y)与 g(x'，y')的相关系数C是关于子区位移和形变参数的函数，令P=(u，v，ux，uy，vx，vy)T，即 C=C(P)。显而易见，当P取得子区真实的位移和形变参量时，相关系数应为最小值。因此可以通过式(6)的极值条件:

求出子区的位移和形变参数。式(7)可以用迭代法求解。令P0为迭代初值，P'为式(7)的解，对相关系数C(P')在P0点做泰勒级数展开，得:

牛顿-拉夫森迭代法需要用其他算法获得整像素位移值作为迭代初值，否则将严重影响收敛速度。

3 计算结果及讨论

采用数值模拟的仿真散斑图对上述3种算法的精度进行了比较。激光散斑图来自于磁流体干燥过程实验，如图 3所示［14］。

图3 磁流体干燥实验中的激光散斑图Fig.3 Typical laser speckle diagram of the magnet fluid drying

在相关搜索算法中需要首先选择模板窗口，选取大模板窗口可以降低噪声的影响，提高计算精度，但大的窗口意味着要计算更多的点，因此计算效率会降低。

Zhang在文献［15］中分析了窗口大小对计算精度和计算效率的影响，他指出窗口选取在31 pixel×31 pixel到51 pixel×51 pixel之间可以在精度和效率之间获得较好的平衡。本文选取窗口为 41 pixel×41 pixel。

类似文献［9］方法，截图以512×512为参考，不失一般性地模拟出平移0.01～0.09 pixel(行方向上)，间隔0.01 pixel的 9幅散斑图和平移0.1～0.9 pixel(列方向)，间隔为0.1 pixel的9幅位移散斑图，共模拟出1幅基准散斑图和18幅位移散斑图。

然后，分别用3种算法计算18幅相对于基准散斑图的位移散斑图，得到了3种亚像素算法偏差均值图如图4所示和标准差图如图5所示。

从图4中可以看出，在这3种算法中梯度算法的偏差均值在整个预设位移范围0.01～0.9 pixels内绝对值都是最小的，在±0.001 pixels之内;牛顿-拉夫森迭代算法在位移为0.01～0.1 pixels时较插值算法的均值误差小，在位移为0.1～0.9 pixels时，牛顿-拉夫森迭代算法和插值算法的平均误差相当，达到±0.003 pixels。但就处理激光散斑图的算法而言，梯度算法处理精度最好，插值法的处理精度最差，牛顿-拉夫森迭代算法的处理精度居中。标准差表示数据的误差分布，标准值越小，误差分布越小，故可用它来评价算法的稳定性。从图5标准差曲线可以看出，梯度算法在3种算法中标准差值最小，但在0.5 pixel位移图计算时是个例处，插值算法的标准差最大。因此，在3种算法中，梯度算法最稳定，而插值算法最不稳定。至于梯度算法在位移0.5 pixels处产生了较大的标准差，潘兵在文献［16］给出的解释:式(3)在做泰勒展开时的截断误差为o(h)2，这是位移为0.5 pixel时达到最大值所致。

图4 3种算法的均值误差Fig.4 Means bias error of three methods

图5 3种算法的标准差Fig.5 Standard deviation of three methods

表1 3种算法的计算时间(以插值算法的用时为归一单位)Tab.1 Computing consuption of three methods

表1给出了3种算法的相对时间比较，3种算法是用Matlab编程实现的，计算机的CPU为Intel Pentium D，主频2.8 GHz，内存2 G。通过比较中可以看出，梯度算法和插值算法的用时相当，略高于插值算法;而牛顿-拉夫森迭代算法用时最长，约是插值算法的9倍。这是由于牛顿-拉夫森算法中每次迭代都要对相关系数矩阵进行插值运算，而在插值算法和基于梯度的算法中只进行了一次插值运算。综合比较3种算法的精度、稳定性和用时，可以看出，在处理激光散斑图时，梯度算法优于其它两种算法。

4 结论

本文介绍了数字散斑图像相关法的位移测量原理以及较常用的插值算法、梯度算法和牛顿-拉夫森迭代算法3种亚像素位移算法及特点。同时利用算法得到了激光散斑图计算的偏差均值曲线和标准差曲线。利用偏差均值的物理意义，比较出梯度算法在亚像素计算散斑图时的精度较高。通过标准差曲线得出梯度算法在亚像素计算散斑图时的稳定性较好。尽管梯度算法用时略高于插值法，但对于激光散斑测量小位移实验或连续刚性位移实验，梯度算法比较有效，应作为优先考虑的算法。

当像素位移量高达0.5时，梯度算法产生了较大误差。此时牛顿-拉夫森迭代法和插值法在0.5～0.9 pixels区间的偏差均值较小，对于激光散斑测量大位移或振动实验，它们的标准差在0.003 pixels之内。

［1］ PETERS W H，RANSON W F.Digital imaging techniques in experimental stress analysis［J］.Opt Eng.，1981，21:427-431.

［2］ 邢冬梅，李鸿琦，李林安，等.梯度功能材料稳态温度场下的热应力分析［J］.功能材料，2002，33(5):505-507.XING D M，LI H Q，LI L A，et al..Study on thermal strain field of functional gradient material in steady temperature［J］.Functional Material，2002，33(5):505-507.(in Chinese)

［3］ 云礼宁，马少鹏，李霞镇，等.用数字散斑相关方法研究竹材在拉伸载荷下的断裂行为［J］.北京理工大学学报，2011，31(3):258-261.YUN L N，MA SH P，LI X ZH，et al..Experimental study of the fraxturc behavior of bamboo under tension based on digital speckle correlation method［J］.Transactions of Beijing Institute of Technology，2011，31(3):258-261.(in Chinese)

［4］ CHU T C，RANSON W F，SUTTON M A，et al..Applications of digital image-correlation techniques to experimental mechanics［J］.Exp.Mech.，1985，25:232-244.

［5］ PEARSON J J，HINES D C，GOLOSMAN S.Video-rate image correlation processor［J］.SPIE，1977，119:197-205.

［6］ DVORCHENKO V N.Bounds on(deterministic)correlation functions with applications to registration［J］.IEEE Trans Pattern Anal Mach Intell，1983，5(2):206-213.

［7］ 潘兵，续伯钦.数字图像相关中亚像素位移测量的曲面拟合法［J］.计量学报，2005，26(2):128-134.PAN B，XU B Q.Subpixel registration using qradratic surface fitting in digital image correlation［J］.Acta Metrologica Sinica，2005，26(2):128-134.(in Chinese)

［8］ BRUCK H A，MCNEILL S R，SUTTON M A，et al..Digital image correlation using Newton Raphson method of partial di6erential correction［J］.Exp Mech，1989，29:261-267.

［9］ AVIS C Q，FREEMAN D M.Statistics of subpixel registration algorithms based on spatiotemporal gradients or block matching［J］.Opt Eng.，1998，37(4):1290-1298.

［10］ ZHOU P，GOODSON K E.Subpixel displacement and deformation gradient measurement using digital image/speckle correlation(DISC)［J］.Opt.Eng.，2001，40(8):1613-1620.

［11］ CHEN D J，CHIANG F P，TAN Y S，et al..Digital speckle-displacement measurement using a complex spectrum method［J］.Appl.Opt.，1993，32:1839-1849.

［12］ ORIAT L，LANTZ E.Subpixel detection of the center of an object using a spectral phase algorithm on the image［J］.Pattern Recognition，1998，31:761-771.

［13］ 李庆扬，王能超.易大义编著，数值分析［M］.北京:清华大学出版社，2008.LI Q Y，WANG N CH，YI D Y.Numerical Analysis［M］.Beijing:Tsinghua University Press，2008.(in Chinese)

［14］ 李新忠，岱钦，王希军，等.用动态散斑法研究磁流体薄膜的干燥过程［J］.光电子·激光，2007，18(2):959-962.LI X ZH，DAI Q，WANG X J，et al..Dynamic speckle method for studying the drying process of magnet fluid film［J］.J.Optoelectronics·Laser，2007，18(2):959-962.(in Chinese)

［15］ ZHANG J，JIN G CH，MA SH P，et al..Application of an improved subpixel registration algorithm on digital speckle correlation measurement［J］.Optics Laser Technology，2003，35:533-542.

［16］ 潘兵，谢惠民，戴福隆.数字图像相关中亚像素位移测量算法的研究［J］.力学学报，2007，39(2):245-252.PAN B，XIE H M，DAI F L.An investigation of sub-pixel displacements registration algorithms in digital image correlation［J］.Chinese J.Theoretical and Appl.Mechanics，2007，39(2):245-252.(in Chinese)