基于光流法和卡尔曼滤波的多目标跟踪方法与流程

本发明涉及一种目标跟踪算法，尤其涉及一种基于光流法和卡尔曼滤波的多目标跟踪方法，属于计算机视觉领域。
背景技术：
：视觉是人类从自然环境中获取信息的主要手段之一，在可见光的照射下，周围环境在人眼的视网膜上形成图像，经感光细胞转换为神经脉冲信号并通过神经纤维传递至大脑皮层进行处理和分析。而计算机视觉则是计算机模拟人眼的视觉功能，从图像中提取有效信息，对环境中的背景和目标进行形态学和运动学方面的分析。基于图像处理技术的目标跟踪是当前计算机视觉领域的前沿课题之一，其研究内容涉及图像分析、计算机视觉、模式识别、机器学习、自动化控制等多个领域。同时，目标跟踪有着非常广阔的应用前景，主要涉及智能视频监控、智能车辆辅助驾驶、视觉导航、虚拟现实、武器制导等方面。因此，对目标跟踪的研究具有很重要的理论和实际意义。随着计算机和各种多媒体设备的普及为研究者们使用各种图像采集设备获取外界环境信息并转换成便于存储的数字信号提供了便利，从而使得电脑可以像人脑一样处理各种图像信息，进而代替人类完成一些与视觉信息相关的工作。目标跟踪是指通过确定目标在视频图像中的位置信息，对目标的运动状态进行估计，同时预测其在下一帧中的运动特征，进而根据所获取的特征对视频中的目标进行匹配，从而达到连续准确跟踪目标的目的。就当前目标跟踪的方法而言，根据是否依赖先验知识可以分为两大类，即依赖目标的先验知识和不依赖目标的先验知识。前者需要为目标建模，然后在后续图像序列中进行模板匹配，进而找到感兴趣目标。而后者则可以直接从图像序列中检测到运动目标，结合相应的搜索算法进而实现对运动目标的有效跟踪。技术实现要素：鉴于以上缺陷，本发明研究了光流法和卡尔曼滤波。其中，光流法是一种对真实运动场的近似估计，在比较理想的情况下，它不需要预先知道场景的任何信息就能够检测独立运动的目标。此外，由于连续帧的时间间隔很短，各个运动目标的运动可近似为匀速运动，而且卡尔曼滤波是一种使用状态方程和观测方程来对系统的状态序列进行线性最小方差估计的算法。它基于系统以前的状态序列对下一个状态作最优估计，预测时具有无偏、稳定、最优、计算量小等特点，并且可以准确地预测目标的位置和速度。因此，本发明结合基于图像金字塔的Lucas-Kanade光流法和卡尔曼滤波两种方法的优势将它们加以融合用于目标跟踪，过程中使用改进的中值滤波算法对获取的图像进行滤波处理以便更准确的获取所跟踪目标。为了实现上述目的本发明采用了如下技术方案：基于光流法和卡尔曼滤波的多目标跟踪方法，包括以下步骤：(1)读取视频帧，采用基于图像金字塔的L-K光流法计算每帧图像的光流。(2)对光流进行聚类处理得到若干光流类，获得标有光流信息的光流检测效果图。(3)对光流检测效果图采用改进的中值滤波方法进行去噪。(4)光流类的个数为运动目标的个数，将这些运动目标的信息传递给卡尔曼滤波进行目标跟踪处理。在上述方案中，所述基于图像金字塔的L-K光流法的具体过程包括：在最高一层的图像上计算得出光流；将计算的结果作为下一层图像的初始值，在这个初始值的基础上计算本层的光流；重复这一过程，直到传递给最后一层，即原始图像层。在本发明的具体实施例中，所述出光流的计算方法为：在一个以a点为中心的局部邻域上定义以下函数，并使该函数值最小：其中，Ω表示点a的局部邻域，W(x,y)表示权函数，表示图像在点a处的梯度，Va表示点a的光流，It表示点a＝(x,y)的灰度值为I＝(x,y,t)图像在t时刻的时域导数。对上式进行最优化求解可得：表示包含n个点的梯度的列向量；W＝diag[W(x1),W(x2),...,W(xn)]表示包含n个点的权值的对角矩阵；b＝-[It(x1),It(x2),...,It(xn)]Τ表示包含t时刻n个点的时域导数的向量；V＝[AΤW2A]-1AΤW2b表示所求的光流信息；表示图像在第n个点处的梯度，W(xn)表示图像在第n个点处的窗口权重，It(xn)图像第n个点的灰度值在t时刻的时域导数，其中xi∈Ω且i＝1,2,3...n。所述聚类处理包括：设定一个阈值Dth，比较两个光流矢量的度量函数D，若此度量值小于阈值则将这两个光流合并，作为一个光流类，计算该光流类的平均光流；再将其余光流矢量与平均光流的度量函数与阈值比较，小于阈值则将该光流并入该光流类中。在同一光流类中，若某点(xi,yi)距离本光流类中心的距离di≥(μ+2σ)，则删除该点。di＝(xi-cx)2+(yi-cy)2n表示同一光流类中像素点的总数。此外，μ是均值也叫期望，同时是正态分布的位置参数，描述正态分布的集中趋势位置，正态分布以x＝μ对称。σ表示随机变量的标准差，也表示正态分布数据的分散程度。σ越大表示数据越分散，σ越小表示数据越集中。同时，σ也是正态分布曲线的形状参数，其值越大，分布曲线也就越扁平，其值越小，曲线也就越高。上述度量函数D为(u1,v1)和(u2,v2)分别表示两个光流矢量值。具体地，所述步骤(3)中还包括以下步骤：将光流检测效果图转化为灰度图；通过改进的中值滤波对灰度图像进行去噪；将去噪后的图像转化为二值化图像；通过膨胀处理和形态学闭的操作获取最终的运动目标。进一步，所述改进的中值滤波采用快速的并行中值滤波方法：对滑动窗口内的每一列或每一行分别计算最大值、中值和最小值后得到滑动窗口行数或者列数组数据；滑动窗口的最大值是最大值组中的最大值，滑动窗口的最小值则是最小值组中的最小值，滑动窗口的中值是最大值组中的最小值、中值组中的中值和最小值组中的最大值这三个数的中间值。本发明的优点：其一，采用了非依赖目标先验知识的方法，所以无需训练分类器和构建目标模板。其二，光流法可以很好地捕获运动目标，通过聚类处理后能够更好的用光流信息标记运动目标。其三，通过改进中值滤波的处理，能够很好的滤出部分噪声，而且比传统中值滤波速度快了近2倍。最后就是结合卡尔曼滤波的预测和修正过程的不断迭代从而达到动态平衡且迭代速度较快，而且可达到很好的跟踪效果。附图说明图1为本发明具体流程图；图2为运动目标经处理后的光流图；图3为中值滤波效果；图4为部分跟踪效果图。具体实施方式下面结合附图对本发明的实施方案进行详细说明。本发明方法算法流程如图1所示。具体实施流程如下：步骤1：使用基于图像金字塔的Lucas-Kanade光流法处理读入的视频帧。光源照射在物体的表面，会使其表面的灰度呈现一定的空间分布，人们在观察运动物体时会在视网膜上形成连续变化的图像即为光流。光流表达图像的变化，包含了运动目标的信息，可用来确定目标的运动。L-K光流首先假设在一个小的空间邻域上运动矢量保持恒定，然后使用加权最小二乘法来估计光流场。假设在t时刻点a＝(x,y)的灰度值为I＝(x,y,t)，可以导出光流约束方程：其中，u和v是该点光流的水平和垂直分量，代表着目标的运动信息，对式(1)进行简化：Ixu+Iyv+It＝0(2)Ix，Iy分别表示点a＝(x,y)的灰度值为I＝(x,y,t)在x轴和y轴方向的空域梯度场，It表示点a＝(x,y)的灰度值为I＝(x,y,t)图像在t时刻的时域导数。将式(2)写成式(3)的形式：用点a处的梯度向量和点a的光流向量Va的点乘再加上It的和为零来表示：其中，表示图像在点a处的梯度；Va＝(u,v)是点a的光流。假定在一个以a点为中心的局部邻域中，各点的光流相同，搜索位移量应使得与之对应的点的相关邻域的匹配误差最小，即在此邻域上定义函数为式(4)，并使其函数值最小：其中，Ω表示点a的局部邻域，W(x,y)表示窗口权函数，它使邻域中心部分对光流约束产生的影响比外围部分更大。对式(4)进行最优化求解可得：表示包含n个点的梯度的列向量；W＝diag[W(x1),W(x2),...,W(xn)]表示包含n个点的权值的对角矩阵；b＝-[It(x1),It(x2),...,It(xn)]Τ表示包含t时刻n个点的时域导数的向量；V＝[AΤW2A]-1AΤW2b表示所求的光流信息；表示图像在第n个点处的梯度，W(xn)表示图像在第n个点处的窗口权重，It(xn)图像第n个点的灰度值在t时刻的时域导数，其中xi∈Ω且i＝1,2,3...n。若单独使用L-K光流法，虽然能够求得光流且对噪声不敏感，但当目标运动尺度过大或者出现遮挡等情况的时候，会因为运动尺度过大而产生较大的计算误差，不但会影响算法的准确度，同时也降低了整体的运算速度。对此，本发明采用了基于金字塔的L-K光流法，其原理描述如下：首先，在最高一层的图像上计算得出光流和仿射变换矩阵。下一层图像的初始值是上一层的计算结果，下一层的图像在这个初始值的基础上，计算本层的光流和仿射变换矩阵。然后重复这一过程，直到传递给最后一层，即原始图像层。最终结果为在原始图像层上计算出来的光流和仿射变换矩阵，最终实现由粗到精的筛选过程。步骤2：对光流进行聚类处理的依据是同一个运动目标上的光流近似且其分布呈现一定规律性，而不同目标上的光流不同，进而可以得到若干类。假设存在两个光流矢量值为(u1,v1)和(u2,v2)，设置衡量它们相似性的度量函数如下：光流聚类可分为粗聚类和精聚类。粗聚类时需要设定一个小的阈值Dth，比较两个光流矢量的度量函数D，若此度量值小于设定阈值就将这两个光流合并，并计算该类的平均光流。再判断其他光流矢量是否并入某类中时，需衡量它们与该类平均光流的度量函数值的大小。另外，精聚类时需要处理已获取的若干光流类。然而，各种噪声的存在使得每类中可能存在一定的离群点，所以就有必要去掉离群点。假设某一光流类中的像素点如：(x1,y1),(x2,y2),...,(xn,yn)由于同一类中的点到类中心的距离服从参数为(μ,σ)的正态分布，若某点(xi,yi)距离本类中心的距离di≥(μ+2σ)，则删除该点，其中各参数如下所示：di＝(xi-cx)2+(yi-cy)2n表示同一类中像素点的总数。此外，μ是均值也叫期望，同时是正态分布的位置参数，描述正态分布的集中趋势位置，正态分布以x＝μ对称。σ表示随机变量的标准差，也表示正态分布数据的分散程度。σ越大表示数据越分散，σ越小表示数据越集中。同时，σ也是正态分布曲线的形状参数，其值越大，分布曲线也就越扁平，其值越小，曲线也就越高。运动目标经处理后能被光流信息所明确标记，效果图如图2所示。步骤3：图像处理技术具有很强的针对性，由于应用的不同、需求的不同，处理方法也就不同。此外，影响图像变化的因素很多，如：天气、光照、摄像机抖动、各种噪声等，所以就有必要对图像进行去噪。本发明选择了中值滤波方法进行去噪。首先，在此过程中并非是对光流进行处理，而是对光流检测效果图(如图2所示)的图像转化为灰度图。其次，通过中值滤波对灰度图像进行去噪。然后将去噪后的图像转化为二值化图像。最后，通过膨胀处理和形态学闭的操作来获取最终的目标。在一定条件下，中值滤波可以处理线性滤波器产生的图像细节模糊问题，特别是对扫描噪声或滤波脉冲的处理有非常好的效果。传统的中值滤波采用含有奇数点的滑动窗口，用窗口中各点灰度值的中值来替代指定点的灰度值。对于奇数元素，中值为大小排序后中间的数值，对于偶数元素，中值为排序后中间两个元素灰度值的平均值，而排序算法又是中值滤波的关键。传统算法使用冒泡排序法，假设窗口中像素点的个数为m，则每个窗口需要比较m(m-2)/2次，时间复杂度为O(m2)。此外，每移动一次窗口都进行一次排序，当图像的大小为N*N时，时间复杂度则为O(m2N2)。当图像较大时，计算量很大，较为费时。对此，本发明使用一种快速的并行中值滤波方法。为便于说明，假设一个3*3的窗口，其内部各个像素如表1所示：表13*3窗口内各像素像素第1列第2列第3列第1行F0F1F2第2行F3F4F5第3行F6F7F8对窗口内的每一列分别计算最大值、中值和最小值后得到了3组数据。其中，max表示取最大值，med表示取中值，min表示取最小值。最大值组：Max1＝max[F0,F3,F6]，Max2＝max[F1,F4,F7]，Max3＝max[F2,F5,F8]中值组：Med1＝med[F0,F3,F6]，Med2＝med[F1,F4,F7]，Med3＝med[F2,F5,F8]最小值组：Min1＝min[F0,F3,F6]，Min2＝min[F1,F4,F7]，Min3＝min[F2,F5,F8]此滑动窗口中的最大值是最大值组中的最大值，最小值则是最小值组中的最小值。而中值组中的最小值至少小于5个像素值，中值组中的最大值至少大于5个像素值。此时，剩下3个要比较的像素：最大值组中的最小值、中值组中的中值、最小值组中的最大值，而这三个数的中间值即为该窗口的中值，其表达式如下：MininMax＝min[Max1,Max2,Max3]MedinMed＝med[Med1,Med2,Med3]MaxinMin＝max[Min1,Min2,Min3]Med＝med[MininMax,MedinMed,MaxinMin]该中值滤波方法比传统算法的运算速度提高了近2倍，经处理后的效果图如图3所示，图中白色连通区域为运动目标，黑色区域为背景区域。步骤4：经上述过程处理后，通过光流聚类，得到类的个数为运动目标的个数。得到类中光流矢量的均值为运动目标的速度。将检测到的运动目标信息传递给卡尔曼滤波进行后续的跟踪处理。卡尔曼滤波以状态方程和观测方程为基础，运用递归方法来预测线性系统变化的方法。其状态方程和观测方程如下：xk＝Ak,k-1xk-1+ξk-1(6)zk＝Hkzk+ηk(7)其中，xk-1和xk分别表示k-1时刻和k时刻的状态向量；zk表示k时刻的观测向量；Ak,k-1表示从k-1时刻到k时刻的状态转移矩阵；Hk表示观测矩阵；ξk或ξk-1表示不同时刻系统噪声，且ξk∈N(0,Qk)；ηk表示观测噪声，且ηk∈N(0,Rk)；Qk和Rk分别是系统噪声ξk和观测噪声ηk的方差。卡尔曼滤波可概述为状态预测过程和状态修正过程，涉及如下部分：状态向量预测方程：误差协方差预测方程：卡尔曼滤波增益：修正状态向量：修正误差协方差矩阵：Pk＝Pk,k-1-KkHkPk,k-1(12)其中式(8)和式(9)为状态预测过程，式(9)、式(10)及式(11)为状态修正过程。预测以状态方程(6)为基础，求取状态预测向量和误差协方差预测向量Pk,k-1。而修正则以观测方程(7)为基础，修正状态预测向量，得到向量并通过计算求出最小误差协方差矩阵。视频中相邻帧的目标之间运动可近似为匀速运动，运动学公式如下，其中△t为相邻帧的时间间隔：Sk＝Sk-1+△tVk-1(13)Vk＝Vk-1(14)Sk表示目标经过时间间隔△t后的位移，Vk表示目标当前的运动速度，由于是近似匀速运动，所以每一时刻速度基本一致。卡尔曼滤波的状态向量如下：其中，x(k)和y(k)是运动目标的中心点坐标，Vx(k)和Vy(k)是目标中心在两个坐标轴方向上的速度分量。由此，其状态转移矩阵如下：根据目标中心点的坐标信息可以得到观测向量和观测矩阵如下，xz(k)和yz(k)表示第k帧中标记目标的矩形框的中心点坐标：将式(15)至(18)带入式(6)和(7)整理后便可重新描述状态方程和观测方程。其中ξk-1和ηk分别为4*1维和2*1维的系统噪声和观测噪声，它们也都是均值为0的高斯白噪声。系统噪声和观测噪声的协方差矩阵Q和R分别设定如下：设定式(20)中观测噪声两个分量的方差则观测噪声协方差矩阵R即为单位矩阵。此外，误差协方差矩阵初始值可设定如下：因相邻帧中目标的运动可近似为匀速运动，为便于计算，可令帧间间隔△t＝1，则状态向量初始值可表示如下：x0为卡尔曼滤波的初始状态，x(0)和y(0)表示为初始状态时目标的中心点坐标，x(k)-x(k-1)和y(k)-y(k-1)分别表示第K时刻与第K-1时刻x轴方向和y轴方向的位移。卡尔曼滤波跟踪是状态预测和状态修正的不断迭代过程中寻求一种动态平衡，当迭代次数达到设定值或者取得局部最优时，输出最终的结果，部分效果图如图4所示。步骤5：判断视频是否结束。假设待处理视频有N帧，当前帧号为i，通过判断当前帧号i与N的大小。如果i小于N，即视频没有结束，则转到步骤1，然后继续循环重复上述过程。如果i大于或者等于N，则表示视频已经被处理完毕，并结束程序。当前第1页1 2 3