光镊系统中一种亚像素位移测量方法研究

张晓林++赵志丹++王晓玲

摘 要：光镊系统中光阱刚度系数的快速、准确、高精度测量是光镊进行粒子定量操控的关键因素。本文给出了一种利用数字图像相关匹配的方法实现亚像素的位移测量。为提高测量精度，此方法在捕获过程中采集多幅图像并以相邻前一幅图像作为模板去匹配下一幅图像，以此类推，最后将每个增量位移场叠加作为最终的位移，同时在亚像素计算时采用较高精度的梯度法；为了提高计算速度，在整像素点匹配时，采用了动态阈值序贯相似法。最后，通过模拟仿真实验，对此算法的有效性进行了验证，结果表明本算法比普通亚像素算法快1～2倍。

关键词：亚像素位移测量 数字图像相关法 序贯相似法 动态阈值 梯度法

中图分类号：TP2 文献标识码：A 文章编号：1672-3791（2014）01（b）-0001-04

光镊（Optical tweezers）又称为单光束梯度力光阱，是一种利用高度聚焦的激光束形成的三维梯度势阱来捕获、操控微小粒子的技术[1]。光镊自1986年由Arthur Ashkin[2]发明以来，以其非接触、地损伤等优点，已被广泛应用于物理学中的激光冷却、胶体化学、生物医学尤其是分子生物学等领域[3]，成为一项重要的研究工具。

光镊的一个重要功能为微小力的测量，对光阱的刚度进行标定是光镊测力的重要环节。标定光阱刚度有许多方法，常用的有流体力学法、热运动分析法、功率谱法和外加周期驱动力法等[4～6]，文献[7]详细分析比较这四种方法的优缺点，其中流体力学法和热运动分析法均需要CCD跟踪微粒运行轨迹，在微粒运行一段期间内拍摄大量的图像，然后再对些图像进行后期处理。由此可知在硬件条件一定的情况下图像亚像素分析对测量精度有着重大的意义。

亚像素位移测量的算法主要有如下几种：亚像素灰度插值法[8]、曲面拟合法[9]、相关系数插值法，牛顿-拉普森[10]（Newton-Rapshon，简称N-R）、基于梯度的方法；频率相关法，后验概率算法，神经网络方法和基于迭代的最小二乘法[11]。这些算法测量精度所称精度能到0.005～0.1pixel。常用的算法就三种，下面就简单的总结这三种算法的优缺性：亚像素灰度插值法，计算量大，精度较低，一般较少直接使用；牛顿-拉普森是基于最优化的思想，建立合理的位移和变形模型后然后进多次迭代后然后求出其中的参数，由于迭代过程中要用到灰度插值，以及灰度的梯度插值，因此目前这个算法是精度最高的，但耗时也最长。相关系数曲面拟合法，有着很强的抗噪性能，但计算精度比相关系数插值法略低。梯度法（又称微区统计特性梯度法），其基本思想是微小物体的近似刚性位移后微小变形前后点对应点的灰度值保持不变。潘兵等人，指出梯度法与曲面拟合法具有相同的效率，且精度优于相关系数曲面拟合[12]。

本文用图像相关法[12]并利用位移场的连续性，设计了一套分步计算整像素位移、压像素位移，最后将两者叠加作为最终的位移，如图1。在具体的计算过程中将集标准化协方差相关法，动态阈值序贯相似法，梯度法众多优点集于一体。最后通过仿真实验验证算法的有效性和实验。

1 本文算法的计算流程

1.1 整像素位移测量

为了提高计算的效率和速度这里采用了动态阈值序贯相似法（SSDA）[13]。下面给出动态阈值SSDA算法大致流程：

（1）定义绝对误差：

其中，，。

（2）在相邻图像中前一幅中取选定一块大小和位置合适的图像作为模板中心坐标为；

（3）确定在后面一幅图像中的搜索范围（即子图的遍历范围）；

（4）在后面一幅图像中，计算模板图像与初始位置子图中所有像素点的的累加值，并将其作为阈值的初始值；

（5）计算模板和一个位置子图中对应点的并累加记作；

（6）在计算并累加过程中比较与的大小，若在计算完每一行或一列后就立刻比较与，若，则停止计算，并将图像子移动到下一个位置，重复（4）计算，加快匹配速度；

（7）若再遍历模板图像与该位置子图的所有像素点后，有，则用T更新，并记录此子图中心点的坐标。

下面给出算法的具体流程如图2所示。

1.2 亚像素位移的测量

为了获取更高的精度，需要在正像素结果的基础上进一步进行亚像素位移的求解，设变形前的图像为，变形后的图像为，分别为对应于原图像中所求位移点在变形的图像中对应点的整像素位移，为对应于整像素位移结果的亚像素位移。

当选物体作微小位移时，且物体表面上任意一点在周围的领域内的元面积足够小，则小面元可以看成近似刚体运动，亦元面内所有的点的均匀相同的位移量。根据数字图像基本假设，在微区内，和有下面的关系：

同时，考虑到相关搜索对应变的不敏感性，设真实位移为：

定义是微区内的刻画与的相似程度函数：

在微区内，和相似程度最大，应该满足（4）式，此时应该取驻值。

将（4）式待入中并进一步表示的函数如下：

对应真实的微小变形应有：

本文中选取了公式（9）作为相似程度函数：

式中，为模板在点点处的灰度值，是变形子区在点处灰度值；分别是模板区域与变形子区的中像素灰度值的平均值。详细的推导过程，见文献[14]。

1.3 位移场的叠加

位移场叠加时考虑相邻图像之间的位置传递。例如，根据图像1和图像2，可以算出图像2相对图像1的位移增量为 ，其中是以图像1为参考系的坐标。同理由图像2、图像3，得出，其中是以图像2为坐标系，简单推导后得到图像3相对图像1的位移场为：

据此类推，可以得到任意一个图像相对于第一副图像的位移。

2 数值模拟

采用斯坦福大学Peng Zhou等人提出的算法[15]，生成标准散斑图如图3，每一幅图大小为512×512，散斑的光强程高斯分布，散斑尺寸大小为4个像素，散斑数为1200，在水平方向移动0～0.1pixel像素内以0.01步长生成9幅图像，在0.1～1像素位移范围内以0.1为步长，生成9幅图像。同时在每幅图像上取上5个不同位置的采用41×41模板，然后经行统计分析，如图4所示。endprint

由图4可以看出，在理想的条件下，本文的亚像素有0.005pixel精度，单从精度要求上，满足光镊光镊中对亚像素速位移精度的要求。同时发现在0～0.5pixel时，计算位移值大于设定位移值，而在0.5～1pixel时计算位移小于设定的位移值。设定位移在0.5pixel时出现反转现象。

第二组实验，生成10幅散斑图像，每幅图像沿y轴上移动移动3.25 pixel，在x轴上外加一个随机0.1*rand的小位移量作为步进电机的运作时的扰动，同时兼顾电磁噪声和其他噪声，在后续的图像中加均值为0，方差为0.01+0.01*rand高斯噪声和噪声密度0.01+rand*0.01的椒盐噪声如图5所示。表1，采用相邻的模板匹配得到的不同帧数上的匹配点的坐标。图5，微粒的不同时刻的相对位置显示。

由第一幅到最后一幅在y轴理论位移量为256+9×3.25=285.250；本文提出的算法最后结果为：285.195，误差值为0.055pixel；而直接用第一幅和最后一幅相关算出的位移为：285.097；与理论值相差为0.1530pixel。显然本文提出的算法有较强的抗干扰性。

第三组实验，是比较在整像素点搜索采用动态阈值的SSDA的和普通的搜索的时间，而亚像素点的计算采用相同算法。对比，实验结果如下表2。表格中的运行时间是在CPU为Intel（R） Core（TM）2 T5870，主频为2.00 GHz处理器，内存大小为2 G，计算5次匹配所画的时间。实际上程序运行的时间取决于程序效率，编程语言、以及计算机硬件设备。本文算法基本比普通亚像素算法快1～2倍。倍数相差不明显的原因是：亚像素计算的时间占整个计算时间很大的部分。

3 结论

由上面的分析可以得出，本文给出一种了利用图像相关分析法通过整数像素位移场计算、亚像素位移计算、和位移场的叠加，来实现大位移场的高精度的亚像素位移测量的方法。在整像素的计算时采的动态的阈值的序贯相似法能使整个匹配过程所花时间节省1～2倍；第二组实验中看出，粒子运动时提高采集图像的频率能的改善测量精度，本文直接用起始位置和终点位置的图像计算出的位移误差是采集多张图像和计算相邻的位移最后按1.3所述叠加后位移出差的近似9倍。此外，虽然本文尽力考虑了各种的干扰影响，步进电机的震动，电磁脉冲的等，但是实际实验中的影响远远不如此，如细胞各自的布朗运动、焦平面的变化、光照的不均匀等。下一步工作应在光镊捕获微粒实验中检验本算法的性能。

参考文献

[1] 余娜，蔡志岗，梁业旺，等.光镊系统的组建及光阱效应的观察[J].大学物理，2010，29（3）：59-62.

[2] Ashkin A，Dziedzic JM， Bjorkholm JE，Chu S.Observation of a single-beam gradient force optical trap for dielectric particles[J].Opt Lett，1986，11（5）：288-290.

[3] 任洪亮，庄礼辉，李银妹.双光镊测量胶体微粒间相互作用势[J].中国激光，2008，35（1）：151-155.

[4] Florin E.-L，Pralle A，Stelzer E.H.K，et al.Photonic force microscope calibration by theraml noise analysis[J].J.Appl.Phys.A.1998，66：S75-S78.

[5] John Bechhoefer，Scott WilsonFaster cheaper，safer optical tweezers for the undergraduate laboratory[J].Am.J.Phys.2002，70（4）：393-400.

[6] 龚契，陈洪涛，李银妹.四种光阱刚度测量法的实验研究与比较[J].中国科大学报，2005，35：601.

[7] Bamea D I Silveman H E A class of algorithms for digital image reg istration[J].IEEE，1972，C-21（2）：179-186.

[8] 于起峰.基于图像的精密测量与运动测量[M].北京：科学出版社，2002.

[9] 王怀文，亢一郎，谢和平.数字散斑相关方法与应用研究进展[J].力学进展，2005，32（2）：195-203.

[10] Bruck HA McNeil SR ，Sutton MA，et，al.Bigital Image correlation using Newton-Rapshon method of partial differential correction[J].Exprimental Mechanics，1989，29（3）：261-267.

[11] 潘兵，吴大方，国宝桥.数字体图像相关方法中基于迭代最小二乘法的物体内部变形测量[J].实验力学，2011，26（6）：666-673.

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

[13] 杜德生，叶建平，等.一种新的自适应阈值SSDA算法[J].现在电子技术，2010，6：135-139.

[14] 张军，金光昌，马少鹏，等.基于微区域统计特性的数字散斑相关测量亚像素位移梯度算法[J].光学技术，2003，29（4）：467-472.

[15] Zhou P，Goodson KE.Subpixel displacement and deformation gradient measurement using digital image speckle correlation（DISC）[J].Opt.Eng.2001，40（8）：1613-1620（August 2001）.endprint

由图4可以看出，在理想的条件下，本文的亚像素有0.005pixel精度，单从精度要求上，满足光镊光镊中对亚像素速位移精度的要求。同时发现在0～0.5pixel时，计算位移值大于设定位移值，而在0.5～1pixel时计算位移小于设定的位移值。设定位移在0.5pixel时出现反转现象。

第二组实验，生成10幅散斑图像，每幅图像沿y轴上移动移动3.25 pixel，在x轴上外加一个随机0.1*rand的小位移量作为步进电机的运作时的扰动，同时兼顾电磁噪声和其他噪声，在后续的图像中加均值为0，方差为0.01+0.01*rand高斯噪声和噪声密度0.01+rand*0.01的椒盐噪声如图5所示。表1，采用相邻的模板匹配得到的不同帧数上的匹配点的坐标。图5，微粒的不同时刻的相对位置显示。

由第一幅到最后一幅在y轴理论位移量为256+9×3.25=285.250；本文提出的算法最后结果为：285.195，误差值为0.055pixel；而直接用第一幅和最后一幅相关算出的位移为：285.097；与理论值相差为0.1530pixel。显然本文提出的算法有较强的抗干扰性。

第三组实验，是比较在整像素点搜索采用动态阈值的SSDA的和普通的搜索的时间，而亚像素点的计算采用相同算法。对比，实验结果如下表2。表格中的运行时间是在CPU为Intel（R） Core（TM）2 T5870，主频为2.00 GHz处理器，内存大小为2 G，计算5次匹配所画的时间。实际上程序运行的时间取决于程序效率，编程语言、以及计算机硬件设备。本文算法基本比普通亚像素算法快1～2倍。倍数相差不明显的原因是：亚像素计算的时间占整个计算时间很大的部分。

3 结论

由上面的分析可以得出，本文给出一种了利用图像相关分析法通过整数像素位移场计算、亚像素位移计算、和位移场的叠加，来实现大位移场的高精度的亚像素位移测量的方法。在整像素的计算时采的动态的阈值的序贯相似法能使整个匹配过程所花时间节省1～2倍；第二组实验中看出，粒子运动时提高采集图像的频率能的改善测量精度，本文直接用起始位置和终点位置的图像计算出的位移误差是采集多张图像和计算相邻的位移最后按1.3所述叠加后位移出差的近似9倍。此外，虽然本文尽力考虑了各种的干扰影响，步进电机的震动，电磁脉冲的等，但是实际实验中的影响远远不如此，如细胞各自的布朗运动、焦平面的变化、光照的不均匀等。下一步工作应在光镊捕获微粒实验中检验本算法的性能。

参考文献

[1] 余娜，蔡志岗，梁业旺，等.光镊系统的组建及光阱效应的观察[J].大学物理，2010，29（3）：59-62.

[2] Ashkin A，Dziedzic JM， Bjorkholm JE，Chu S.Observation of a single-beam gradient force optical trap for dielectric particles[J].Opt Lett，1986，11（5）：288-290.

[3] 任洪亮，庄礼辉，李银妹.双光镊测量胶体微粒间相互作用势[J].中国激光，2008，35（1）：151-155.

[4] Florin E.-L，Pralle A，Stelzer E.H.K，et al.Photonic force microscope calibration by theraml noise analysis[J].J.Appl.Phys.A.1998，66：S75-S78.

[5] John Bechhoefer，Scott WilsonFaster cheaper，safer optical tweezers for the undergraduate laboratory[J].Am.J.Phys.2002，70（4）：393-400.

[6] 龚契，陈洪涛，李银妹.四种光阱刚度测量法的实验研究与比较[J].中国科大学报，2005，35：601.

[7] Bamea D I Silveman H E A class of algorithms for digital image reg istration[J].IEEE，1972，C-21（2）：179-186.

[8] 于起峰.基于图像的精密测量与运动测量[M].北京：科学出版社，2002.

[9] 王怀文，亢一郎，谢和平.数字散斑相关方法与应用研究进展[J].力学进展，2005，32（2）：195-203.

[10] Bruck HA McNeil SR ，Sutton MA，et，al.Bigital Image correlation using Newton-Rapshon method of partial differential correction[J].Exprimental Mechanics，1989，29（3）：261-267.

[11] 潘兵，吴大方，国宝桥.数字体图像相关方法中基于迭代最小二乘法的物体内部变形测量[J].实验力学，2011，26（6）：666-673.

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

[13] 杜德生，叶建平，等.一种新的自适应阈值SSDA算法[J].现在电子技术，2010，6：135-139.

[14] 张军，金光昌，马少鹏，等.基于微区域统计特性的数字散斑相关测量亚像素位移梯度算法[J].光学技术，2003，29（4）：467-472.

[15] Zhou P，Goodson KE.Subpixel displacement and deformation gradient measurement using digital image speckle correlation（DISC）[J].Opt.Eng.2001，40（8）：1613-1620（August 2001）.endprint

由图4可以看出，在理想的条件下，本文的亚像素有0.005pixel精度，单从精度要求上，满足光镊光镊中对亚像素速位移精度的要求。同时发现在0～0.5pixel时，计算位移值大于设定位移值，而在0.5～1pixel时计算位移小于设定的位移值。设定位移在0.5pixel时出现反转现象。

第二组实验，生成10幅散斑图像，每幅图像沿y轴上移动移动3.25 pixel，在x轴上外加一个随机0.1*rand的小位移量作为步进电机的运作时的扰动，同时兼顾电磁噪声和其他噪声，在后续的图像中加均值为0，方差为0.01+0.01*rand高斯噪声和噪声密度0.01+rand*0.01的椒盐噪声如图5所示。表1，采用相邻的模板匹配得到的不同帧数上的匹配点的坐标。图5，微粒的不同时刻的相对位置显示。

由第一幅到最后一幅在y轴理论位移量为256+9×3.25=285.250；本文提出的算法最后结果为：285.195，误差值为0.055pixel；而直接用第一幅和最后一幅相关算出的位移为：285.097；与理论值相差为0.1530pixel。显然本文提出的算法有较强的抗干扰性。

第三组实验，是比较在整像素点搜索采用动态阈值的SSDA的和普通的搜索的时间，而亚像素点的计算采用相同算法。对比，实验结果如下表2。表格中的运行时间是在CPU为Intel（R） Core（TM）2 T5870，主频为2.00 GHz处理器，内存大小为2 G，计算5次匹配所画的时间。实际上程序运行的时间取决于程序效率，编程语言、以及计算机硬件设备。本文算法基本比普通亚像素算法快1～2倍。倍数相差不明显的原因是：亚像素计算的时间占整个计算时间很大的部分。

3 结论

由上面的分析可以得出，本文给出一种了利用图像相关分析法通过整数像素位移场计算、亚像素位移计算、和位移场的叠加，来实现大位移场的高精度的亚像素位移测量的方法。在整像素的计算时采的动态的阈值的序贯相似法能使整个匹配过程所花时间节省1～2倍；第二组实验中看出，粒子运动时提高采集图像的频率能的改善测量精度，本文直接用起始位置和终点位置的图像计算出的位移误差是采集多张图像和计算相邻的位移最后按1.3所述叠加后位移出差的近似9倍。此外，虽然本文尽力考虑了各种的干扰影响，步进电机的震动，电磁脉冲的等，但是实际实验中的影响远远不如此，如细胞各自的布朗运动、焦平面的变化、光照的不均匀等。下一步工作应在光镊捕获微粒实验中检验本算法的性能。

参考文献

[1] 余娜，蔡志岗，梁业旺，等.光镊系统的组建及光阱效应的观察[J].大学物理，2010，29（3）：59-62.

[2] Ashkin A，Dziedzic JM， Bjorkholm JE，Chu S.Observation of a single-beam gradient force optical trap for dielectric particles[J].Opt Lett，1986，11（5）：288-290.

[3] 任洪亮，庄礼辉，李银妹.双光镊测量胶体微粒间相互作用势[J].中国激光，2008，35（1）：151-155.

[4] Florin E.-L，Pralle A，Stelzer E.H.K，et al.Photonic force microscope calibration by theraml noise analysis[J].J.Appl.Phys.A.1998，66：S75-S78.

[5] John Bechhoefer，Scott WilsonFaster cheaper，safer optical tweezers for the undergraduate laboratory[J].Am.J.Phys.2002，70（4）：393-400.

[6] 龚契，陈洪涛，李银妹.四种光阱刚度测量法的实验研究与比较[J].中国科大学报，2005，35：601.

[7] Bamea D I Silveman H E A class of algorithms for digital image reg istration[J].IEEE，1972，C-21（2）：179-186.

[8] 于起峰.基于图像的精密测量与运动测量[M].北京：科学出版社，2002.

[9] 王怀文，亢一郎，谢和平.数字散斑相关方法与应用研究进展[J].力学进展，2005，32（2）：195-203.

[10] Bruck HA McNeil SR ，Sutton MA，et，al.Bigital Image correlation using Newton-Rapshon method of partial differential correction[J].Exprimental Mechanics，1989，29（3）：261-267.

[11] 潘兵，吴大方，国宝桥.数字体图像相关方法中基于迭代最小二乘法的物体内部变形测量[J].实验力学，2011，26（6）：666-673.

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

[13] 杜德生，叶建平，等.一种新的自适应阈值SSDA算法[J].现在电子技术，2010，6：135-139.

[14] 张军，金光昌，马少鹏，等.基于微区域统计特性的数字散斑相关测量亚像素位移梯度算法[J].光学技术，2003，29（4）：467-472.

[15] Zhou P，Goodson KE.Subpixel displacement and deformation gradient measurement using digital image speckle correlation（DISC）[J].Opt.Eng.2001，40（8）：1613-1620（August 2001）.endprint