基于微粒子群优化算法的数字散斑图像相关方法

梁智锦，王开福，顾国庆，张成斌

（南京航空航天大学航空宇航学院，南京210016）

基于微粒子群优化算法的数字散斑图像相关方法

梁智锦，王开福*，顾国庆，张成斌

（南京航空航天大学航空宇航学院，南京210016）

为了能够通过一步搜索同时得到数字散斑图像中所测点的整像素和亚像素位移信息，采用灰度插值的方法构造了亚像素子区，改进了基于微粒子群算法的数字图像散斑相关方法。对含有平移信息的模拟散斑图和具有应变的模拟散斑图进行相关计算，验证了该方法的适用性；在对具有微小面内位移转动的试件进行测量时，比较了整像素的微粒子群算法和不同量级的灰度插值下的亚像素微粒子群算法。结果表明，基于微粒子群算法的亚像素数字散斑图像相关方法在测量小位移方面具有一定的优越性。

图像处理；数字散斑相关；微粒子群优化算法；灰度插值

引 言

微粒子群算法是由KENNEDY和EBERHART于1995年提出来的一种模仿鸟类捕食的智能搜索算法［1］，其算法结构简单、易于实现。微粒子群算法的主要思想是：在被测目标的周围散落着一定数量的粒子，这群粒子像鸟一样的在被测区域内移动，每个粒子分别记录了自己所经历的最佳位置并与其它粒子的最佳位置进行比较以记录所有粒子经历的全局最佳位置，之后粒子将根据自己的历史最佳位置和全局最佳位置决定自身的下一步移动目标，通过这样不断的迭代使得粒子向全局最佳位置靠近并最终达到最佳值［2］。

数字散斑相关是数字散斑计量技术的一种，在20世纪80年代由YAMAGUCHI［3］和PETERS等人［4］提出来的位移和变形测量技术。与数字散斑条纹技术不同，其实质是根据物体变形前后散斑场的互相关性来获取物体的位移和变形信息，根据寻求相关系数的极值位置，在变形后的散斑场中识别出对应于变形前的散斑场［5-6］。数字散斑相关技术的关键在于图像相关算法。相关算法从最开始的粗细搜索法发展到现在，已经有很多种，其中的典型算法主要有：双参量法、牛顿迭代法、梯度搜索法、遗传算法和差分进化算法等［7-9］。DU和WANG曾经在2012年发表过基于粒子群算法的整像素数字图像相关方法［10］，但是在很多微小位移的测量中，整像素的测量范围并不能满足测量的需要。2009年，天津大学CHANG在其硕士毕业论文中提到使用微粒子群算法进行图像相关的亚像素搜索，但他所采用的方法是在得到整像素位移的情况下通过拟合整像素的相关系数曲线后利用微粒子群算法进行曲线极值的搜索［11］。作者是在微粒子群算法应用于数字散斑图像相关运算基础上，采用灰度插值的方式使得微粒子群算法可以直接进行精确到亚像素的数字图像相关运算。

1 基于微粒子群算法的数字图像相关算法的基本原理

1.1 微粒子群算法的基本原理

在n维搜索空间中，设粒子群总数为k，则第i个粒子的位置可以表示为Xi=［xi1，xi2，xi3，…，xin］T（T表示转置），速度vi=［vi1，vi2，vi3，…，vin］T，其所经历过程的最佳位置坐标为pi，整个粒子群所经历过程的全局最佳位置坐标为p。每个粒子所处的位置坐标信息都表示一个候选解，将候选解的值代入目标函数计算其适度值以确定候选解的优劣度，并据实时更新pi和p。

微粒子群算法可表示为［2］：

式中，vid，t+1表示迭代至（t+1）次时第i个粒子第d维的速度，其中t=1，2，3…N，N为最大迭代次数，d=1，2，3…n；w为惯性常数，其具体值根据代数t的增加由公式得出；wmax一般取值1.4，wmin取值0；c1和c2为加速常数，一般取值均为2；r1和r2为随机数，其范围为［0，1］；pd，t为截止到第t代时全局最佳位置点的第d维坐标值；pid，t为截止到第t代时第i个粒子所处的最佳位置的第d维坐标值；xid，t表示第t代时，第i个粒子第d维坐标值。

1.2 数字图像相关的基本原理

数字图像相关计算是指位于散斑场上的两个子区进行匹配运算的相关方法。当两个子区位于同一个散斑场时，图像相关计算被称为自相关。自相关算法主要运用于粒子图像测速中。当子区位于不同散斑场时，图像相关算法被称为互相关算法，也就是一般所指的数字图像相关。

通过CCD记录被测物体变形或者运动前后的散斑图，在散斑图上每点的位置信息用该点周围区域的散斑分布来表示，其散斑分布的实质就是点周围区域（称为子区）组成的灰度值矩阵。在变形前的散斑图上选取被测点子区I，在变形后的散斑图上选取待测点子区Ii′，两者进行相关运算得到相关系数C［5］：

式中，〈…〉表示系综平均。相关系数表明两个子区的相关程度；当相关系数为1时，说明两个子区完全相关，当相关系数为0时，说明两个子区完全不相关。

1.3 基于微粒子群算法的数字散斑图像相关算法流程

微粒子群算法的相关流程是：（1）在搜索域内随机的分布k个粒子，并将表示代数的t设置为1；（2）通过适应度函数判断这所有粒子的适应度，比较获得每个粒子的适应度历史最大值并记录其位置坐标为pi；（3）通过比较pi的适应度值大小，得到其中的最大适应度值的位置坐标，即为p；（4）通过（1）式计算每个粒子的速度，并确定下一代粒子的新坐标；（5）将p坐标下的适应度值和适应度函数的阈值进行比较，如果此时的p满足阈值要求，则其坐标值即为所求目标。若不满足，则比较此时代数t和最大迭代数N，如t=N，则此时的p坐标也为此次所求目标；如t＜N，则令t=t+1，将所得到的新粒子作为第t+1代粒子并转入步骤（2）中进行循环迭代。

将微粒子群算法应用于数字散斑图像相关运算中的主要改造过程见下。

（1）微粒子群的每一个粒子在图像相关中代表一个子区Ii′（xi，yi），该子区用其中心坐标（xi，yi）来表示。而图像相关算法中的搜索域在变形或运动后的图像中已所测点坐标为中心而设置的一定大小的区域。

（2）相关运算中的微粒子群算法适应度函数是在物体的变形或运动之前的图像中选取以被测点为中心坐标的子区I（x，y）与粒子子区Ii′（xi，yi）之间利用（2）式求取相关系数C。在图像相关运算中C值越大两个子区相关程度越高，当C=1的时候表示两个子区完全相关，C=0则表示两个子区完全不相关；因此，粒子的适应度可以直接用C的大小表示，同时适应度的阈值也被设为1。

（3）通过（1）式计算粒子运动速度的时候，由于公式中的参量非整且存在除法，粒子的速度也不可能全是整数。这导致了后面的粒子子区不会全为整像素子区。如果使用微粒子群整像素图像相关运算的话，就需要将得到的粒子运动速度进行取整运算以确保每个粒子子区均为整像素子区。对于亚像素子区的获取，主要采取灰度插值的方式。本文中采用的是双线性插值法，其插值过程可表示为［5］：

式中，⌊x⌋表示向下取整；I0，I1，I2和I3表示中心坐标离坐标（x，y）最近的4个子区，因为对子区I（x，y）中任意一个灰度点来说，它都包含在这4个相邻子区的对应点中且距离左上角对应点的距离均为Δx和Δy，所以可以直接用子区来进行灰度插值运算而不需要对每个点单独插值计算。同时为了防止粒子的跳跃性太大，超出范围。粒子的速率被限制在±5个像素之间，当粒子的速率超出这个范围时，粒子的速率将被定为边界值。

2 对微粒子群图像相关算法的实验

2.1 微粒子群算法对模拟散斑图的测量

为检验此亚像素插值算法对整像素位移的测量正确性，用计算机制作模拟散斑图并加入整像素位移值，如图1所示。其中图1a为散斑参考图，图1b为加入了整像素位移值的模拟散斑图，图中的U场和V场位移量为8pixel。采用微粒子群整像素图像相关运算算法对其中的U场位移进行测量可以得到图2所示结果。图2中，U场位移的浅灰色平面中的深黑色属于算法中所测得的坏点，这些坏点的表现形式为其搜索出来的值远远超过周围区域所搜索出的正常值，通过重复性的使用算法对这些坏点进行搜索可以得到正确值，说明其出现的主要原因为在搜索区域中对随机粒子的选取没有做到足够均匀且散斑场的部分区域比较相似造成的。但需要说明的是这些坏点在所测区域100pixel×100pixel的10000个点中不超过500个，此算法的正确率达到了95%以上，所以该算法在图像相关中是适用的。图3所显示的U场位移是采用微粒子群亚像素图像相关算法对图1a和图1b进行相关运算所得到结果，在所得的结果中同样的存在不超过500个的坏点，算法的正确率也达到了95%以上。

为检测此算法对物体变形信息测量的正确性，用计算机制作模拟散斑图（见图4），其中图4b为图4a加入V场位移关于x方向的应变0.05（单位为1）后的模拟散斑图。亚像素微粒子群算法程序对图4a和图4b中的50pixel×20pixel区域内进行数字图像相关运算，所得V场位移结果如图5所示，由于相关运算算法中没有考虑到子区本身也发生了一定变形，图中所得的数据呈现出一定得起伏，但其数值的总体计算结果仍然能够得到x方向上的V场应变0.05（单位为1）。

2.2 微粒子群算法对面内旋转的数字散斑图的测量

采用如图6所示实验系统，当试件发生面内转动（即物面围绕试件表面的法线转动）时，利用CCD相机记录试件的面内位移。

实验中所获得的散斑图如图7所示，其中图7b为图7a中的试件绕垂直于试件表面的z轴进行顺时针面内旋转后的试件散斑图（其中旋转中心在散斑图的右上方）。

根据面内旋转位移的特点，V场位移和x坐标具有与旋转角度有关的线性关系［11］。因此旋转的角度可以根据下列公式求得［12］：

对这两幅散斑图分别进行整像素的微粒子群算法测量其所得到的结果如图8所示。由图8可以看出，由于旋转的角度非常小，V场的位移也非常小，这时整像素的图像相关算法只能在200pixel× 200pixel的计算区域中获得两个阶梯值，同时由于V场的变化缓慢，两个阶梯值的边界也模糊不清。

在采用0.1pixel插值的微粒子群亚像素图像相关算法对图7a和图7b进行测量其所得到的结果如图9所示。虽然在200pixel×200pixel的计算区域内所获得的也是阶梯型的结果，但由于所获得的阶梯较多，已经能够计算出V场位移与x的线性关系，将数据代入（4）式中经过计算可得转动角度θ= 0.003846rad。

图10为采用0.01pixel插值的亚像素数字图像相关运算所得到的结果，在这个量级下就能很明显的得到转动角度θ=0.003855rad，其与上面0.1pixel所得结果误差不超过0.3%。

图11显示是在0.001pixel灰度插值下的V场位移图。可以发现图10和图11相比较来说在整体趋势上和单体数据上都是基本一致的。

对于任意散斑图来说，其插值精度取决于测量的精度要求。对本文中测量面内旋转角度来说，在200pixel×200pixel测量范围内散斑图在处理后所得到的位移信息需要含有多个阶梯值，而整像素测量无法达到这个要求，就需要插值到0.1pixel。如果所获得的测量数据更少，以至于0.1pixel插值不能得到两个以上阶梯值，那就要通过进一步的插值来获取位移信息。

3 结 论

（1）为了扩展微粒子群算法在数字散斑图像相关中的应用，采用灰度插值的方法构造了亚像素的子区，并同时编写单独测量整像素和在测量整像素的同时测量亚像素的MATLAB程序，应用带有面内位移信息模拟散斑图验证了两个程序的正确性。

（2）在应用带有应变信息的模拟散斑图对亚像素微粒子群算法程序进行验证的时候，证明了亚像素微粒子群算法在对具有微小应变的情况下依然适用。

（3）同时应用整像素微粒子群算法程序和亚像素微粒子群算法程序对面内微小转动的试件进行了测量，在整像素微粒子群算法无法计算出转动的微小角度时成功地运用亚像素微粒子群算法测量出了面内转动的角度；从而比较说明了面对微小位移时整像素微粒子群算法的局限性和亚像素微粒子群算法的有效性。同时对不同量级的灰度插值的亚像素微粒子群算法进行了比较，表明了在测量微小位移信息时亚像素微粒子群的灵活性和多样性。

［1］ KENNEDY J，EBERHART R C.Particle swarm optimization［C］//The University of Western Australia.IEEE International Conference on Neural Networks.Perth，Australia：The University of Western Australia，1995：1942-1948.

［2］ XIE X F，ZHANG W J，YANG ZH L.Overview of particle swarm optimization［J］.Control and Decision，2003，18（2）：129-134（in Chinese）.

［3］ YAMAGUCHI I.A laser speckle strain gauge［J］.Journal of Physics，1981，E14（11）：1270-1273.

［4］ PETERS W H，RANSON W F.Digital imaging technique in experimental stress analysis［J］.Optical Engineering，1982，21（3）：427-431.

［5］ WANG K F，GAO M H，ZHOU K Y.Modern photo mechanics［M］.Harbin：Harbin Institute of Technology Press，2009：58-61，154-156（in Chinese）.

［6］ XU X，WANG K，GU G Q，et al.Measurement of internal material flaws based on out-of-plane displacement digital speckle pattern interferometry［J］.Laser Technology，2012，36（4）：548-552（in China）.

［7］ JIN G C.Computer-aided optical measurement［M］.2nd ed.Beijing：Tsinghua University Press，2007：141-147（in Chinese）.

［8］ LIU M.The improved genetic algorithms for digital image correlation technique［D］.Tianjin：Tianjin University，2005：4-6（in Chinese）.

［9］ PAN B，XIE H M.Digital image correlation method with differential evolution［J］.Journal of Optoelectronics·Laser，2007，18（1）：100-103（in Chinese）.

［10］ DU Y Zh，WANG X B.Digital image correlation method based on particle swarm optimization algorithm without sub-pixel interpolation［J］.Computer Engineering and Applications，2012，48（6）：200-204（in Chinese）.

［11］ CHANG Y M.Research and application of sub-pixel search algorithm for digital image/speckle correlation method［D］.Tianjin：Tianjin University，2009：12-14（in Chinese）.

［12］ GU G Q，WANG K F，XU X.Measurement of in-plane microrotations and rotation-center location of a rigid body by using digital image correlation［J］.Chinese Journal of Lasers，2012，39（1）：0108004（in Chinese）.

Digital speckle image correlation method base on particle swarm optimization algorithm

LIANG Zhijin，WANG Kaifu，GU Guoqing，ZHANG Chengbin
（College of Aerospace Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China）

In order to get the integer pixel displacement information and the sub-pixel displacement information simultaneously at one time，the sub-pixel subarea was constructed by using the method of gray interpolation and the digital image correlation method based on particle swarm optimization algorithm was proposed.The applicability of the method was verified by measuring both the simulated speckle pattern with the translational information and with the strain information.And then，the integer pixel particle swarm algorithm and the sub-pixel of different magnitude interpolation gray particle swarm algorithm were compared by measuring the specimen with tiny rotational displacement.The results show that the particle swarm algorithm based on sub-pixel digital scattered spot image correlation method has advantages for small displacement measurement.

image processing；digital speckle correlation；particle swarm optimization algorithm；gray interpolation

TN911.73

A

10.7510/jgjs.issn.1001-3806.2014.05.006

1001-3806（2014）05-0603-05

国防技术基础科研资助项目；江苏省高校优势学科建设工程资助项目

梁智锦（1989-），男，硕士研究生，研究方向为光测力学及图像处理。

*通讯联系人。E-mail：kfwang@nuaa.edu.cn

2013-09-16；

2013-12-02