0 引言

目标检测和跟踪是当今计算机视觉领域的复杂、热点研究课题之一。随着经济发展，交通工具不断增加，公路承载能力相对不足，交通事故频发^[1-2]。智能交通系统中运动目标检测跟踪技术结合图像处理、计算机和自动化控制技术，实时监控交通目标，提取目标位置、形态结构、运动参数等信息，对改善交通状况提供决策依据。近年来光流技术^[3]在模式识别、计算机视觉以及图像处理领域广泛受到关注。光流法基于图像连续序列，抽取光流场，由于目标和背景存在相对运动，运动速度矢量不同于相邻区域，根据像素点速度矢量特征，进行运动估计，进而检测目标及其位置区域。

Horn等^[4]最早引入光流约束方程，提出二维速度场同灰度相联系，发展至今包括基于梯度的方法、基于匹配的方法、基于能量的方法及基于相位的方法^[5]。基于匹配的方法对相邻帧字块匹配进行运动估计，特征提取和精确匹配较难；基于能量和相位的方法跟踪效果较好, 但计算负荷大, 实时性较差；基于梯度的方法研究最多，以图像灰度不变为条件导出光流约束方程，如Horn-Schunck、Lucas-Kanade(LK)^[6]、Nagel^[7]等算法。Tao等^[8]提出通过样本的稀疏子集估计光流场，以减少计算量。Parvathy等^[9] 定义光流信息的运动轨迹模型，由此检测运动目标的行为。文献[10]提出一种快速计算的基于稀疏光流的物体跟踪方法, 通过计算图像特定像素光流矢量实现物体检测、跟踪。文献[11]提出基于平面假设，通过光流跟踪算法实现帧间小位移定位的卡尔曼方法和特征点匹配融合的算法。张志禹等^[12]提出基于SURF与光流法相结合的增强现实局部跟踪注册方法, 该方法首先通过光流法跟踪目标，然后提取SURF特征，最后建立三维注册矩阵，提高了三维注册的速度和精度，然而该方法仅对光流跟踪更新，没有将其与SURF联合制定更新策略。

基于梯度的光流方法虽实现简单，但存在以下问题：1) 基于光流灰度守恒的制约，该方法不适用于光照变化的情况；2) 光流法仅关注像素点，忽略运动目标和像素点的关联；3) 对环境噪声敏感，跟踪快速移动物体易失败；4) 受摄像机抖动等因素影响，形成的异常信号会导致产生不连续形变场。

针对以上问题，本文以城市道路交通复杂环境为背景，提出基于多分辨率LK光流算法联合快速鲁棒性特征(Speed Up Robust Feature, SURF)的跟踪方法，首先利用高斯背景建模提取前景目标，接着利用SURF局部尺度不变特征变换算法检测前景目标特征点，然后通过构建图像多分辨率小波金字塔改进LK稀疏光流法对特征点进行检测和跟踪，并制定自适应模板实时更新策略，根据模板更新方法实时更新跟踪窗口。新方法对车辆形变、高速、雾霾天气、光照不均、部分遮挡等复杂环境有较高鲁棒性，增强抗噪能力的同时提高运算效率和精度，实现目标稳定跟踪。

1 LK光流法

LK算法属于稀疏光流算法，相比稠密光流法跟踪过程中需要计算图像所有像素点光流矢量，LK算法仅需计算代表图像特征的像素点的光流矢量，在跟踪过程中具有良好的实时性，匹配准确，复杂度较低。算法假设像素邻域空间运动矢量相同, 计算相邻图像中相应像素点的流向信息，通过获取特征点周围小窗口局部信息匹配光流矢量。LK算法如式(1) 所示：

$\sum\limits_{\lambda \in \Omega } {{W^2}\left( x \right){{[\nabla I\left( {x,{\rm{ }}t} \right) \cdot v + {I_t}\left( {x,{\rm{ }}t} \right)]}^2} = 0} $

(1)

其中：邻域内各像素点的权重用W(x)窗口加权函数表示；Ω是光流一致的含有n个点的邻域范围。通过最小二乘法求解式(2) ，最小化能量函数，从而求解光流得式(3) 。

${\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{W}}^2}\mathit{\boldsymbol{A}}V = {A^{\rm{T}}}{\mathit{\boldsymbol{W}}^2}\mathit{\boldsymbol{B}}$

(2)

$V = {({A^{\rm{T}}}{\mathit{\boldsymbol{W}}^2}\mathit{\boldsymbol{A}})^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{W}}^2}\mathit{\boldsymbol{B}}$

(3)

其中：$\mathit{\boldsymbol{A}} = {(\nabla I({x_1}),\nabla I({x_2}),\nabla I({x_3}), \cdots \nabla I({x_n}))^{\rm{T}}}$, $\mathit{\boldsymbol{W}} = {\rm{diag}}(W({x_1}),{\rm{ }}W({x_2}),{\rm{ }} \cdots ,W({x_n}))$，$\mathit{\boldsymbol{B}} = - {({I_t}({x_1}),{I_t}({x_2}), \cdots ,{I_t}({x_n}))^{\rm{T}}}$。A^TW²A表达如式(4) ：

${\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{W}}^2}\mathit{\boldsymbol{A}} = \left[ \begin{array}{l} \begin{array}{*{20}{c}} {\sum {{W^2}\left( x \right)I_x^2\left( x \right)} } & {\sum {{W^2}\left( x \right){I_x}\left( x \right){I_y}\left( x \right)} } \end{array}\\ \begin{array}{*{20}{c}} {\sum {{W^2}\left( x \right){I_y}\left( x \right){I_x}\left( x \right)} } & {\sum {{W^2}\left( x \right)I_y^2\left( x \right)} } \end{array} \end{array} \right]$

(4)

利用线性系统求解得到光流场。令δI为相邻图像I_t+1(x, y)、I_t(x, y)在(x, y)处的差，I_t+1(x, y)、I_t(x, y)分别为t+1和t时刻像素(x, y)的灰度值。δI=I_t+1(x-d_x, y-d_y)-I_t(x, y)，δ表示在t时刻，邻域像素在t+1产生位移量为[dx dy]^T的移动，当δ取最小值时即得位移量。光流残差函数δ表示如下：

$\[\delta =\iint_{W}{{{[{{I}_{t+1}}(x-{{d}_{x}},\text{y-}{{\text{d}}_{\text{y}}})-{{I}_{t}}\left( x,\text{y} \right)]}^{2}}\mathsf{dxdy}}\]$

(5)

虽然LK算法有良好的实时性, 但通用性不强，如在跟踪过程中摄像机运动造成抖动, 使物体在相邻两帧之间位移较大，或跟踪快速移动物体时会出现目标误匹配的现象, 导致丢失目标，跟踪精度下降，其抗噪能力也有待增强。

2 联合跟踪方法 2.1 多分辨率位移改进

由于LK算法中存在同一像素点在相邻两帧间存在较大位移量易导致跟踪失败的问题，基于多分辨率思想, 本文构建图像多分辨率小波金字塔，在不同尺度求取位移量，实施LK算法，利用小波金子塔层与层间的位移传递计算，使得不同尺度上的位移量相比原始位移量缩小。

基于卷积的Mallat 离散小波变换算法利用图像水平、垂直方向上交替使用低通及高通滤波实现，其计算复杂度高，实时性差。本文采用提升小波变换算法，信号的高频分量通过多项式插补求得，信号的低频分量通过构建多尺度函数求得，不依赖傅里叶变换，将其引入多分辨率金字塔中，可以将二阶滤波变换或离散小波变换分解为系列提升步骤，通过分解现有的小波滤波器完成小波变换，逆小波变换实现也快速、简单，将正向变换步骤求反即可。图像n层提升小波分解变换后可得不同尺度的1个低频分量及3n个高频分量。

通过提升小波算法在多尺度分解图像, 构造提升小波金字塔，令I^n-1(x, y)为第n-1层的金字塔图像，Iⁿ(x, y)为上层图像，n为金字塔总层数。令g_m为金字塔中m层初始估计的光流矢量，Δf_m为光流计算矢量值，最上层光流初始值为g_n-1=0。金字塔的各层之间光流映射关系如式(6) 所示，其中k为层间映射关系系数。根据上层图像对下层金字塔图像进行光流计算，得到原图像的运动参数如式(7) 所示：

${g_{m - 1}} = k({g_m} + \Delta {f_m})$

(6)

$f = {g_0} + \Delta {f_0} = \sum\limits_{m = 0}^{n - 1} {{k_m}{f_m}} $

(7)

本文选取3层金字塔，顶层光流量g₂=0，根据V=(A^TW²A)^-1A^TW²B可以求取Δf₂的光流矢量值，代入g_m-1=k(g_m+Δf_m)中，可以得到2层光流量g₁=kΔf₂，依次迭代，进一步可以求得最底层的光流量g₀以及Δf₀的值，f=g₀+Δf₀，因此初始图像的光流值为：

$v\left( {x,{\rm{ }}y} \right) = \sum\limits_{m = 0}^2 {{k_m}\Delta {f_m}} $

(8)

设光流在计算时每次处理的移动距离D_MAX最大为1像素，利用提升小波分辨率金字塔计算后可以将像素最大移动距离提升为(2^m+1-1) D_MAX，当m取3时，构造三层金字塔，代入D_MAX=1，则最大位移量是15像素，因此，通过引入提升小波构造多分辨率变换金字塔，可以解决LK算法不能跟踪快速运动目标的问题。

2.2 联合SURF

SURF^[13]算法保持了尺度不变特征变换(Scale-Invariant Feature Transform, SIFT)算法尺度、仿射不变性的特点，同时匹配效率得到很大提高。本文高斯背景建模后，提取运动目标，将SURF算法提取的目标特征点用多分辨率小波金字塔改进的LK稀疏光流法计算运动估计。由于仅对局部运动矢量特征点处理，可以减少特征点的提取、匹配区域，相对全局光流算法，极大降低冗余数据，减少了光流计算量，增强了系统实时性。此外，通过制定自适应模板更新策略，可以有效解决因长时间跟踪致使累积误差导致的跟踪失败问题。

常用的Harris、Susan、Moravec角点检测算子的描述信息比较单一，特征点集中在局部区域，对尺度、光照、旋转、噪声、缩放等较为敏感。SURF算法在借鉴SIFT简化近似的思想的同时计算速度提高近3倍，该算法简化了高斯二阶微分模板，使用积分图像生成特征矢量，从而避免特征矢量在生成的时候的重复运算，包括以下步骤：

1) 构建尺度空间检测特征点。

使用Hessian矩阵行列式在尺度空间寻找局部极大值。图像某点的Hessian矩阵定义如下：

$\mathit{\boldsymbol{H}}\left( {x,{\rm{ }}\sigma } \right) = \left[ \begin{array}{l} \begin{array}{*{20}{c}} {{L_{xx}}\left( {x,{\rm{ }}\sigma } \right)} & {{L_{xy}}\left( {x,{\rm{ }}\sigma } \right)} \end{array}\\ \begin{array}{*{20}{c}} {{L_{xy}}\left( {x,{\rm{ }}\sigma } \right)} & {{L_{yy}}\left( {x,{\rm{ }}\sigma } \right)} \end{array} \end{array} \right]$

(9)

其中:σ为尺度，L_xx(x, σ)为x点高斯函数二阶导数$\frac{{{\partial ^2}}}{{\partial {x^2}}}g\left( \sigma \right)$和图像的卷积。L_xy(x, σ)，L_yy(x, σ)类似。如图 1所示，图 1(a)为SIFT尺度空间构建的金字塔，图像尺寸变化，滤波器大小不变；图 1(b)为SURF尺度空间中滤波器大小改变，图像保持不变。

令行列式为det(H)=L_xx×L_yy-L_xy²，简化高斯二阶微分模板，使用积分图像卷积。取σ=1.2，建立9×9尺度空间大小的盒式滤波模板对图像滤波，检测特征点。图 2为y方向以及x和y方向的盒式滤波器的近似表示。使用滤波器后的Hessian矩阵行列式近似为:

${\rm{det}}\left( \mathit{\boldsymbol{H}} \right) = {D_{xx}} \times {D_{yy}} - {(0.9 \times {D_{xy}})^2}$

(10)

D_xx、D_xy、D_yy为滤波模板和图像卷积后的结果。选择不同尺度Hessian检测算子的盒式滤波器构建多尺度空间金字塔图像，寻找不同尺度的特征点。

2) 主方向与特征描述算子。

特征点作中心，设s是特征点所在的尺度值，计算半径6s(s=1.2) 邻域内的点在x和y方向的Haar小波响应，同时设置高斯权重系数。为了求取主方向，以步长0.2弧度左右旋转π/3范围扇形滑动窗口，对范围内的Haar小波响应值累加，得到新的矢量，整个圆形区域遍历后主方向为最长矢量对应的方向，即最大Haar小波响应累加值对应的方向，选取过程如图 3所示。

将坐标轴旋转至主方向，围绕特征点，取边长为20s的方形区域，将其划分成4×4的子窗口，如图 4所示，每个子窗口有5 s×5 s个像元，利用尺度是2 s的小波对子窗口图像计算响应值，采样后得到垂直主方向的dx和沿主方向的dy，以特征点作为中心进行高斯加权计算，计算各个子窗口的响应值，其矢量表示如下：

${V_{子窗户}} = \left[ {\sum dx,{\rm{ }}\sum \left| {dx} \right|,{\rm{ }}\sum dy,{\rm{ }}\sum |dy|} \right]$

(11)

共4×4个子窗口，特征描述算子形成4×(4×4) =64维特征矢量后再对向量进行归一化处理，相比SIFT算法，维度减少一半，大幅度提高了匹配效率。

3) 特征匹配。

本文在匹配SURF特征矢量过程中采用改进的K-维树(K-dimensional tree, K-ｄ树)最近邻优先(Best Bin First, BBF)查询机制提高查询效率。BBF将查询的节点按分割超平面和查询点的距离进行排序，能确保优先检索到含有最近节点可能性高的空间，从最高优先级的节点回溯检查，当优先级队列中节点都检查完毕或超时的时候，返回找到的最优结果作为最近邻。前后两帧中计算A图中某特征点和B图所有特征点的欧氏距离，当最近和次近距离比小于设定阈值0.6时，则特征点匹配。SURF提取的特征点集通过随机抽样一致性(RANdom SAmple Consensus, RANSAC)算法利用内点建立参数模型，可以有效去除不准确的特征点。

针对求解单应性变换矩阵H的非线性优化问题，为了避免照度变换影响，对RANSAC计算的单应变换矩阵H构造梯度误差函数，设(x′, y′)和(x, y)是匹配的特征点对，变换模型和梯度误差函数如式(12) 、(13) 所示：

$\left[ {\begin{array}{*{20}{c}} {{{x'}_a}}\\ {{{y'}_b}}\\ 1 \end{array}} \right] = \mathit{\boldsymbol{H}}\left[ {\begin{array}{*{20}{c}} {{x_a}}\\ {{y_b}}\\ 1 \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{h_{11}}} & {{h_{12}}} & {{h_{13}}}\\ {{h_{21}}} & {{h_{22}}} & {{h_{23}}}\\ {{h_{31}}} & {{h_{32}}} & 1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{x_a}}\\ {{y_b}}\\ 1 \end{array}} \right]$

(12)

$E\left( h \right) = \sum\limits_{i = 1}^n {{{[\nabla I'({x_i}) - \nabla I({x_i})]}^2}} $

(13)

其中:I(x_i)、I′(x_i)分别为变换前后图像，每一对(x_i, y_i)和对应点(x′_i, y′_i)之间的误差为▽I′(x_i)-▽I(x_i)，迭代更新单应性矩阵H，将误差函数E(h)收敛于指定阈值，如0.1，从而得到优化后的变换矩阵。

图 5为不同环境中，两帧通过改进SURF算法匹配的结果：图 5(a)中存在明显仿射变换，且右图增加了噪声，分辨率低；图 5(b)光照变化明显，且左右两图拍摄视角不同；图 5(c)中车辆分别位于阴影处和光照处，且存在旋转变化；图 5(d)中左图车辆前半部分区域被树遮挡，但由于SURF算法无需提供全局信息，因此未遮挡区域仍可较好匹配，体现了SURF较好抗复杂环境的能力。

表 1是SIFT和SURF两种算法对以上不同环境两幅图匹配后的特征点个数和计算时间的比较结果，经过RANSAC优化，两种算法特征点匹配准确率及数量差异不大，但SURF算法计算效率明显有较大优势，时间节省近一半。SIFT算法实现时，需要重设图像尺寸，且每层图像依赖于前一层结果，运算量很大；而SURF算法在多尺度空间同时被处理，无需二次抽样，运算性能得到增强。

2.3 模板更新

在跟踪目标过程中，环境背景时常变化，如果初始模型不进行更新，遇到照度变化、目标形变、旋转、遮挡等情况时因模型偏差极大，影响跟踪效果，因此本文通过制定模板更新方法以提高跟踪算法适用性。

在模板更新过程中，无需每帧都更新，否则会增加系统运行负担。利用高斯背景建模后提取前景目标确定跟踪窗口，对窗口内的SURF特征点分配权值W_j，接近目标的权值较大，接近背景的权值较小，当特征点估计位置超出跟踪窗口一定比例时，更新模板。令窗口内外特征点数分别为N_in和N_out，若N_out/(N_in+N_out)＞K时，则更新模板，文中K取0.25。更新模板时，令跟踪窗口包含所有光流估计特征点，将检测到的SURF特征点和前一帧的跟踪模板进行特征匹配，匹配后按照下式更新所有特征点权值：

$W_{j}^{'}\left( t \right)=C\left( \alpha ({{W}_{j}}\left( t \right)+1/M)+\left( 1-\alpha \right){{W}_{j}}\left( t-1 \right) \right)$

(14)

其中：M为特征点总个数，C为归一化系数，α是为了使特征模板在更新时，降低对背景环境的变化的敏感程度引入的折中系数，本文取0.7。更新结束后，删除掉权值低的N个特征点，提高特征点纯度，实现模板更新，对抗复杂环境有较好鲁棒性。

2.4 算法总体流程

算法步骤：

步骤1 高斯背景建模得到前景目标，建立初始模板，获得跟踪窗口。

步骤2 检测SURF特征点，利用RANSAC算法对内点建立参数模型，去除不准确的特征点，并对RANSAC计算的单应变换矩阵H构造梯度误差函数, 迭代更新单应性矩阵。

步骤3 对SURF算法提取的特征点利用提升小波构造多分辨率变换金字塔改进的LK稀疏光流法计算运动估计，实现LK算法对快速运动目标的跟踪。

步骤4 判断是否最后一帧，若是，结束跟踪；若不是，进行模板更新判断。若不需模板更新，转步骤3继续跟踪；若需要模板更新，根据上文模板更新方法更新跟踪窗口，检测SURF特征点，更新特征点权值，删除低权值特征点，转步骤3继续跟踪。

3 本文算法验证

在交通环境中，利用监控视频在Windows XP SP4、2.4 GHz CPU、2 GB内存、Matlab R2012a环境的计算机中进行实验，视频图像大小为400×320。利用高斯背景建模确定初始跟踪目标，通过优化后的SURF算法求取目标局部不变特征点，对生成的矢量特征点进行LK光流运动估计，同时运用提升小波多分辨率金字塔改进光流运动矢量，并自适应更新模板，完成目标跟踪。

图 6为不同算法在不同环境中的跟踪结果比较。图 6(a)~(e)分别为本文算法测试的高速路段车辆光流场矢量图、运动车辆二值图、本文高速环境跟踪图以及全局SURF算法和传统LK算法高速环境跟踪图，从左到右依次为视频中35帧、76帧、97帧和121帧。图 6(b)中，高斯背景建模进行背景差分提取前景车辆，并用滤波、形态学腐蚀开闭去掉噪声和道路信息。图 6(c)中，本文算法基于多分辨率思想，在不同尺度求取位移量，实施LK算法，利用小波金子塔层与层间的位移传递计算，使得不同尺度上的位移量相比原始位移量缩小，在相邻两帧间存在较大位移量的环境中能够稳定跟踪。图 6(d)、(e)图中全局SURF算法和传统LK算法在高速环境下易丢失目标。图 6(f)~(k)分别为不同算法在雾霾环境和车辆遇到树木遮挡情况下的跟踪图，由于本文算法利用局部SURF算法提取的特征点进行LK稀疏光流法计算运动估计，且自适应更新特征模板，对雾霾环境抗噪能力强，当遮挡情况发生时，未被遮挡的车身仍可检测出来，汽车重现时继续跟踪，有较高的抗遮挡性能。由于全局SURF算法耗时大，虽具有一定抗复杂环境能力，但实时跟踪难以保证，跟踪过程不稳定。传统LK算法抗复杂环境能力差，且在车辆帧间位移较大的情况下，易导致跟踪精度下降，致使目标最终丢失。

表 2为对高噪、雾霾、部分遮挡、照度不均等复杂环境的6组交通视频，每组分成100段，通过全局SURF逐帧跟踪、传统LK算法及本文LK联合SURF局部跟踪方法后进行跟踪平均耗时及跟踪准确率比较。SURF逐帧跟踪对每一帧提取SURF特征，并通过匹配帧间特征点偏移向量跟踪目标，耗时大。LK稀疏光流算法有较高的实时性，但对快速运动、低分辨率、照度变化等情况鲁棒性低。本文算法实时性和LK相差不多，但准确率得到很大程度的提高，抗复杂环境能力较强。

4 结语

传统基于稠密光流的跟踪方法计算量大，难以满足实时性要求；LK算法复杂度较低、实时性好，但存在的对噪声敏感、大尺度运动跟踪易失败的问题。为了将其应用于复杂交通环境车辆的检测跟踪，本文采用提升小波多分辨率算法构建图像多分辨率小波金字塔, 用改进LK算法解决跟踪快速移动物体时会出现目标误匹配的现象, 同时联合SURF特征点匹配算法，制定模板更新规则，在检测特征点的基础上，降低外点概率，通过多分辨率小波金字塔光流算法跟踪这些局部特征点，解决户外交通中由于车辆尺度改变、光照不均、雾霾天气、噪声大等复杂环境引起的跟踪失败问题，提高了运算效率，增强抗干扰能力，鲁棒性高，能稳定跟踪快速运动车辆，在智能交通系统目标跟踪中有良好应用前景。