一种改进的基于Hankel总体最小二乘法的精子多目标跟踪方法

雷茜茜，黄文明，邓珍荣，冷金强

(桂林电子科技大学 计算机科学与工程学院，广西 桂林541004)

检测跟踪运动目标是图像处理和计算机视觉的重要研究内容，它在军事、视频监控和医疗诊断方面都有广泛的应用。多目标跟踪主要有特征建模、贝叶斯估计、光流场与运动场、概率关联匹配等方法。余东等［1］提出了改进的最邻近搜索法，通过计算帧间精子的运动距离和方向预测其在下一帧图像序列中可能出现的位置，并利用跟踪丢失补偿和预测范围最小原则，与传统最近邻算法相比，提高了跟踪的精度，但对于密度大、运动无规律的精子无法进行有效跟踪。汪忠传等［2］提出了一种基于小窗口的应用于超长时间精子序列图像跟踪的最大距离预测法，通过结合最大距离预测法和最近邻搜索算法对精子进行实时跟踪，但该算法和卡尔曼滤波、粒子滤波算法一样，需要假设先验模型预测下一帧目标所在位置。Ding等［3］提出了一种Hankel矩阵低秩逼近求运动相似度的多目标跟踪算法，在未假设运动先验模型和目标外观信息的条件下实现了多目标跟踪，但该算法采用凸松弛和一般内点(generic interior point)法求解Hankel矩阵的秩，计算复杂度和空间复杂度高，不易于实现。Collins［4］和Andriyenko等［5］提出了高阶运动模型算法，通过利用高阶运动模型实现跟踪外观相似的目标，且未假设先验模型。但这些方法仍有限制，文献［4］的时间复杂度高达O(d2.5)，d为检测区的个数;文献［5］的算法需要调试大量的参数。Dicle等［6］提出了Hankel总体最小二乘法(Hankel total least square，简称HTLS)、广义线性分配法(generalized linear assignment，简称GLA)进行关联的多目标跟踪算法，此算法跟踪准确率较高，不需要目标的外观信息和运动的先验模型，可处理任意高阶运动，且只需2个参数，但该算法漏跟踪率较高，处理时间较长，特别是目标密集的图像序列。鉴于此，提出一种改进的Hankel总体最小二乘法。

1 最近邻算法改进的运动相似性模型

在一组精子图像序列中，精子外观都很相似，且轨迹长度不一，故在进行精子轨迹片段关联的过程中面临没有外观信息、目标遮挡、运动轨迹交叉、运动无规律等，造成了跟踪的困难。为解决此问题，结合最近邻算法和Hankel总体最小二乘法对运动精子进行跟踪。

1.1 最近邻算法

最近邻算法是提出最早、最简单的目标跟踪方法，在某些特定情况下是最有效的。该方法首先设置一个关联门用来限制潜在的决策数目，经关联门将所有回波进行初步筛选得到候选回波。关联门是整个跟踪区域中的一个子区域，其中心位于被跟踪目标的预测位置，范围根据在一定的概率上能收到正确回波而设计。最近邻算法是对落在关联门中的点进行选择，根据距离进行判断，选择距离被跟踪目标预测点最近的点。

在图像处理中应用最近邻算法，计算量小、简单、易于实现。通过计算下一帧图像中所有目标与上一帧图像中被跟踪目标的距离，取距离最近者为被跟踪目标在下一帧图像中的位置来实现跟踪［1］。

1.2 Hankel总体最小二乘法

设一个轨迹片段α由一个有序测量值yk组成，

其中，s≤k≤s+l－1，s为开始帧，l为轨迹片段α的长度。对于足够大的n，线性回归器是通用的近似器［7］，可将轨迹片段α的基本运动模型用线性回归器表示为:

回归器的顺序n为运动的复杂度。在没有噪声的情况下，n=rank()，其中为m(m≥n)列的Hankel矩阵:

于是，2个轨迹片段αi、αj之间的运动相似度Pij可定义为［3］:

对于轨迹片段α的长度l，有运动复杂性n＜l，设分别为无噪声坐标序列和噪声。为了便于表示，使Hankel矩阵［A|b］=［］，［E|f］=［］，其中:

E、f为A与b的误差矩阵。由式(1)、(2)可得，(A+E)x=(b+f)。根据TLS通过最小化误差矩阵E、f可以求解x的思想，可求解总体最小二乘(TLS)问题［8］求出α:

其中:○为Hadamard乘，(A+E)x=(b+f)，［A|b］，［E|f］∈SH。引入Ω恢复丢失的数据。

由于［E|f］和［A|b］是对角线元素相等的Hankel矩阵，可将式(4)改写为:

其中:r(η，x)=b+f－(A+E)x=0;D为对角矩阵，其组成元素包含Hankel矩阵中ηi出现的次数。

结合约束条件和最小化问题［8］，得

其中:π为处罚常量;r(η，x)为

p1=［Om×nIm×m］，m=l－n。线性化r(η，x)可得:

于是，存在矩阵X∈Rm×(m+n－1)满足:

其中，p0=［I(m+n－1)×(m+n－1)O(m+n－1)×1］。因此，最终可将式(6)写为:

此式可通过Hankel最小二乘法求解［9］。

1.3 广义线性分配(GLA)关联

假设已知N条轨迹片段{α1，α2，…，αN}，GLA问题可表述为最优化问题［10-11］:

1.4 改进的轨迹片段运动相似度模型

针对Hankel算法计算复杂问题，结合最近邻算法，提出最近邻算法改进的Hankel总体最小二乘运动相似性计算方法。首先计算相邻帧间所有轨迹片段αi和αj(i，j=1，2，…，N)的欧式距离Di，Di为αi与下一帧所有轨迹片段的距离向量。对Di按从小至大顺序进行排序，并记录对应索引Ii:

然后采用式(3)计算距离最近的轨迹片段αi与αIi(1)的相似度，若相似度为零，则计算αi与距离次近的轨迹片段αIi(2)的相似度，以此类推，直至相似度P为1停止计算。最后利用式(10)对相似度为1的轨迹片段进行关联，实现对动物精子的多目标跟踪。

2 实验结果与分析

算法实验在配有Pentium IV 2.93 GHz和2 GB内存的计算机上使用Matlab仿真实现，以采样频率为25帧每秒、经分割识别后的动物精子图像序列为实验对象，其初始图像如图1(精子已编号)［1］。

图1 首帧精子图像Fig.1 The first picture of sperm

为进行有效对比，实验分别运行最近邻算法、文献［6］算法和本算法对25帧动物精子图像序列进行多目标跟踪，获得的精子运动轨迹如图2所示，其中图2(a)、(b)和(c)分别为最近邻算法、文献［6］算法及本算法的跟踪轨迹图，图2(d)为动物精子的实际运动轨迹图。对比分析图2的跟踪轨迹可看出，最近邻算法与精子实际运动轨迹的差别最大，改进算法与精子实际运动轨迹最接近。

各算法的处理时间和准确率结果如表1所示。从表1可见，3种算法中，最近邻算法仅需计算精子质心间的距离，快速简单，跟踪速度最快，仅为3.59 s，但对于精子密集度较高、目标遮挡、运动发生突变的情况无法进行有效处理，使得准确率为90.60%最低。文献［6］算法的跟踪线索为运动状态，不依赖目标外观，也不受外界因素的干扰，准确率较高，达95.58%，但该计算相似度的过程很复杂，算法运行时间较长。本算法结合两者的优点，取距离最近的待跟踪目标与跟踪目标进行相似性匹配，跟踪时间缩短至3.94 s，降低了11.06%，且精子漏跟踪率下降了一定幅度，更快速准确地实现了精子多目标跟踪。

表1 算法运行时间与准确率Tab.1 Time and accuracy of three algorithms

3 结束语

Hankel总体最小二乘法为相似外观的多目标跟踪提供了可靠的跟踪线索。鉴于动物精子外观的相似性，利用文献［6］算法，并在进行运动相似性匹配的同时结合最近邻算法缩小匹配范围，提高跟踪速度。实验结果表明，本算法能稳定、准确且快速地跟踪动物精子的运动，为后续精子质量参数分析的准确性提供了可靠保证。但本算法在进行轨迹跟踪的过程中未解决精子弹性碰撞问题。为满足实际应用的要求，进一步提高准确率仍需今后更深入研究。

图2 3种算法跟踪及实际轨迹Fig.2 Tracked trajectories of three algorithms and real trajectories

［1］余东.猪精子质量自动分析系统应用研究［D］.桂林:桂林电子科技大学，2012:27-39.

［2］汪传忠，任明响，武海燕.最大距离预测法在超长时间精子序列图像跟踪中的应用［J］.测试技术学报，2014，28(2):132-136.

［3］Ding T，Sznaier M，Camps O.Fast track matching and event detection［C］//Anchorage:the IEEE Conference on Computer Vision and Pattern Recognition，2008:1-8.

［4］Collins R.Multi-target data association with higher-order motion models［C］//Providence:the IEEE Conference on Computer Vision and Pattern Recognition，2012:1744-1751.

［5］Andriyenko A，Schindler K，Roth S.Discrete-continuous optimization for multi-target tracking［C］//Providence:the IEEE Conference on Computer Vision and Pattern Recognition，2012:1926-1933.

［6］Caglayan D，Mario S，Octavia C.The way they move:tracking multiple targets with similar appearance［C］//IEEE International Conference on Computer Vision，2013:2304-2311.

［7］Ross G，Soland R.A branch and bound algorithm for the generalized assignment problem［J］.Mathematical Programming，1975，8(1):91-103.

［8］Shmoys D，Tardos E.An approximation algorithm for the generalized assignment problem［J］.Mathematical Programming，1993，62(1):461-474.

［9］Gold S，Rangarajan A.Softmax to softassign:neural network algorithms for combinatorial optimization［J］.Journal of Artificial Neural Networks，1995，2(4):381-399.

［10］Breiman L.Hinging hyperplanes for regression，classification and function approximation［C］//IEEE Transaction on Information Theory，1993:999-1011.

［11］Park H，Zhang L，Rosen J B.Low rank approximation of a Hankel matrix by structured total least norm［J］.BIT Numerical Mathematics，1999，39(4):757-779.

［12］Markovsky I，Huffel S V.Overview of total least squares methods［J］.Signal Processing，2007，87(10):2283-2302.