目标跟踪是完成视觉监控、人机交互、车辆导航等诸多视频场景分析和理解任务的基础, 已有大量的跟踪方法被报道，可将这些方法大致分成单源跟踪^[1-5]和多源跟踪^[6-11]两大类。与单源目标跟踪系统相比，多源目标跟踪系统在生存能力、时空覆盖范围、可信度等方面都具有明显的优势，因而被广泛研究，其中最具代表的是红外与可见光的融合跟踪。红外传感器是通过检测目标辐射的热能差异来形成影像，因此在恶劣的光照环境下要比可见光成像更好，但它无法捕获目标对象的颜色和纹理特征。可见光传感器虽然无法感知温度，但在处理多个热目标交汇时，它通常要优于红外传感器，特别当目标对象间有显著的颜色和纹理差异时。因此通过联合它们的数据，能获得比使用单个传感器更好的跟踪性能。文献[6]提出了基于多个空间直方图连乘的红外与可见光目标表示与决策级融合跟踪方法。文献[7]先利用模糊区域动态融合的方式，将红外与可见光图像进行融合，然后在融合图像中提取多个不同特征，并将它们联合来完成红外—可见光目标跟踪。文献[8]提出采用联合稀疏特征表示的方法对红外和可见光目标进行特征级融合跟踪。文献[9]提出了一种基于局部鉴别分析的红外与可见光目标协同跟踪方法。文献[10]提出了先单独采用粒子滤波法跟踪可见光目标，模板匹配法跟踪红外目标，然后将两者的跟踪结果进行联合决策的先跟踪后融合策略。文献[11]提出了基于Kalman滤波和压缩感知的红外与可见光目标时空融合跟踪算法。

以上方法大多是在粒子滤波框架下实现的，其时间复杂度普遍较高。本文拟采用核跟踪策略来实现快速高效的融合跟踪。核跟踪 (KBT)^[12]算法由于使用梯度下降的迭代搜索技术，简单快速，已被广泛应用于实时性要求较高的目标跟踪场合。虽然，经过十几年的发展，核跟踪方法已相对成熟，但是，作为一种轻量级、易实现的算法，它仍广受人们的持续关注和研究^[13-14]。针对原始KBT方法的空间信息丢失和目标模型更新问题，文献[15]提出了空间直方图 (spatiogram) 概念，文献[16]提出了Kalman滤波的目标模型更新机制。本文利用KBT简单快速、空间直方图能记录目标空间结构信息的优点，提出了一种基于联合空间直方图表示的红外与可见光目标自适应融合跟踪算法。论文的主要贡献有：1) 将核跟踪方法拓展到了多源目标跟踪中；2) 导出了多源目标的联动位移公式；3) 建立了融合权值的自适应调节机制；4) 实现了多源目标模型的协同更新。

空间直方图是一种带高阶矩信息的广义直方图。零阶矩对应传统意义上的直方图，而二阶矩直方图则是附加了像素点空间分布的均值和方差信息的直方图，因此能较好地保持目标的空间结构信息^[15]。为了行文方便，以下简称二阶矩空间直方图为二阶直方图。

记$h\left( z \right){\left\{ {\left\langle {{p_u}\left( z \right), {\boldsymbol{\mu} _u}\left( z \right), {\Sigma _u}\left( z \right)} \right\rangle } \right\}^m}_{u = 1}$为中心在$z$点的候选目标的二阶直方图。${p_u}(z)、{\boldsymbol{\mu} _u}(z)$和${\Sigma _u}(z)$分别为第$u$个特征区像素点的概率密度函数以及这些像素点空间分布的均值和协方差矩阵，它们的计算公式为

$ {p_u}\left( z \right) = C\sum\limits_{i = 1}^n {k\left( {{{\left\| {\frac{{{x_i}-z}}{h}} \right\|}^2}} \right)} {\delta _{iu}} $

(1)

$ {\boldsymbol{\mu} _u}\left( z \right) = \frac{1}{{\sum\limits_{j = 1}^n {{\delta _{ju}}} }}\sum\limits_{i = 1}^n {\frac{{{x_i}-z}}{h}} {\delta _{iu}} $

(2)

$ {\Sigma _u}\left( z \right) = \frac{{\sum\limits_{i = 1}^n {\left( {\frac{{{x_i}-z}}{h}-{\boldsymbol{\mu} _u}\left( z \right)} \right)} {{\left( {\frac{{{x_i}-z}}{h} - {\boldsymbol{\mu} _u}\left( z \right)} \right)}^{\rm{T}}}{\delta _{iu}}}}{{\sum\limits_{j = 1}^n {{\delta _{ju}} - 1} }} $

(3)

式中，$n$是候选目标的像素点个数，${\boldsymbol{x}_i}$是第$i$个像素的二维坐标向量。${\delta _{iu}}$是Kronecker德塔函数，如果第$i$个像素落在第$u$个区间，则${\delta _{iu}}$=1，否则为零。归一化常数$C = 1/\sum\limits_{i = 1}^n {k\left( {{{\left\| {\frac{{{x_i}-z}}{h}} \right\|}^2}} \right)} $。$m$为特征区个数。$h$为带宽，代表候选目标的大小。$k$($x$) 为Epanechni-kov核函数的轮廓函数，计算公式为^[12]

$ k\left( x \right) = \left\{ {\begin{array}{*{20}{l}} {{c^{-1}}_d\left( {d + 2} \right)\left( {1-x} \right)/2}&{x \le 1}\\ 0&{其他} \end{array}} \right. $

(4)

式中，$d$=2是核函数维度，$c_d$为$d$维单位球的体积。

设$\mathit{\boldsymbol{\tilde h}} = {\{ \left\langle {{q_u}, {{\mathit{\boldsymbol{\tilde \mu }}}_u}, {{\tilde \Sigma }_u}} \right\rangle \} ^\mathit{\boldsymbol{m}}}_{\mathit{\boldsymbol{u}} = 1}$为目标模型的二阶直方图，则候选目标与目标模型之间的相似度为

$ \mathcal{B}\left[{\mathit{\boldsymbol{h}}\left( z \right), \mathit{\boldsymbol{\tilde h}}} \right] = \sum\limits_{u = 1}^m {{\psi _u}} \left( z \right)\sqrt {{p_u}\left( z \right){q_u}} $

(5)

式中，$\sqrt {{\boldsymbol{p_u}}\left( z \right){\boldsymbol{q_u}}} $可理解为计算候选目标与目标模型在特征空间中的相似度，而$\boldsymbol{{\psi _u}}$($z$) 则用于计算两者在空间分布上的相似度，其计算公式为^[17]

$ \begin{array}{l} {\psi _u}\left( z \right) = 8{{\rm{ \mathsf{ π} }}^4}\sqrt {|{\Sigma _u}(z){{\tilde \Sigma }_u}|} \\ \mathcal{N}\left( {{{\mathit{\boldsymbol{\tilde \mu }}}_u}\left( z \right);{{\mathit{\boldsymbol{\tilde \mu }}}_u}, {{\hat \Sigma }_u}(z)} \right) \end{array} $

(6)

式中，$\mathcal{N}\left( \cdot \right)$表示正态分布，${{\hat \Sigma }_u}(z) = 2\left( {{\Sigma _u}(z) + {{\tilde \Sigma }_u}} \right)$。

为了表述方便，文中用$\boldsymbol{I}$代表红外，用$\boldsymbol{J}$代表可见光。记$\boldsymbol{I}\left( z \right) = {\{ \left\langle {{p^I}_u\left( z \right), {\boldsymbol{\mu} ^I}_u\left( z \right), {\Sigma ^I}_u\left( z \right)} \right\rangle \} ^m}_{u = 1}$为中心在$z$点的红外候选目标二阶直方图，$\mathit{\boldsymbol{\tilde I}} = {\{ {q^I}_u, {{\mathit{\boldsymbol{\tilde \mu }}}^I}_u, {{\mathit{\tilde \Sigma }}^I}_u\} ^m}_{u = 1}$为红外目标模型。相应地，记$\mathit{\boldsymbol{\tilde J}} = {\{ \left\langle {{q^J}_u, {{\mathit{\boldsymbol{\tilde \mu }}}^J}_u, {{\mathit{\tilde \Sigma }}^J}_u} \right\rangle \} ^m}_{u = 1}$为可见光目标模型，$\mathit{\boldsymbol{J}}\left( z \right) = {\{ \left\langle {{p^J}_u\left( z \right), {\mathit{\boldsymbol{\mu }}^J}_u\left( z \right), {\mathit{\Sigma }^J}_u\left( z \right)} \right\rangle \} ^m}_{u = 1}$为中心在$z$点的可见光候选目标的二阶直方图；仿照式 (5)，可得红外与可见光的候选目标与其目标模型的相似度分别为

$ \mathcal{B}\left[{\mathit{\boldsymbol{I}}\left( z \right), \mathit{\boldsymbol{\tilde I}}} \right] = \sum\limits_{u = 1}^m {{\psi ^I}_u\left( z \right)} \sqrt {{p^I}_u\left( z \right){q^I}_u} $

(7)

$ \mathcal{B}\left[{\mathit{\boldsymbol{J}}\left( z \right), \mathit{\boldsymbol{\tilde J}}} \right] = \sum\limits_{u = 1}^m {{\psi ^J}_u\left( z \right)} \sqrt {{p^J}_u\left( z \right){q^J}_u} $

(8)

要特别说明的是，在经配准后的可见光图像和红外图像中，同一运动对象的尺寸应该是一致的，因此记它在两图像中的像素点个数均为$n$、带宽均为$h$。

同时考虑可见光相似度与红外相似度对目标评判的影响，则可构建联合目标函数为

$ \rho \left( z \right) = \alpha \mathcal{B}\left[{\mathit{\boldsymbol{I}}\left( z \right), \mathit{\boldsymbol{\tilde I}}} \right] + \beta \mathcal{B}[\mathit{\boldsymbol{J}}\left( z \right), \mathit{\boldsymbol{\tilde J}}] $

(9)

式中，0≤$\alpha $, $\beta $≤1为权重系数，用于调节红外相似度与可见光相似度在目标函数中所占的比重，且有$\alpha $+$\beta $=1，0≤$\rho $($z$)≤1。

设目标在前一帧中的位置为$z$₀。将式 (7) (8) 代入式 (9) 中，并在$p_u^I$($z$₀)、$\boldsymbol{\mu} _u^I$($z$₀)、$p_u^I$($z$₀) 和$\boldsymbol{\mu} _u^I$($z$₀)，处对$\rho $($z$) 进行泰勒展开，得到其线性逼近形式

$ \begin{array}{l} \rho \left( z \right) \approx \sum\limits_{u = 1}^m {{w^I}_u\left( {{z_0}} \right)} \left[{{p^I}_u\left( z \right)-{p^I}_u\left( {{z_0}} \right)} \right] + \ldots + \\ \sum\limits_{u = 1}^m {{w^I}_u\left( {{z_0}} \right)} \left[{{\mathit{\boldsymbol{\mu }}^I}_u\left( z \right)-{\mathit{\boldsymbol{\mu }}^I}_u\left( {{z_0}} \right)} \right] + \ldots + \\ \sum\limits_{u = 1}^m {{w^I}_u\left( {{z_0}} \right)} \left[{{p^I}_u\left( z \right)-{p^I}_u\left( {{z_0}} \right)} \right] + \ldots + \\ \sum\limits_{u = 1}^m {{w^I}_u\left( {{z_0}} \right)} \left[{{\mathit{\boldsymbol{\mu }}^J}_u\left( z \right)-{\mathit{\boldsymbol{\mu }}^J}_u\left( {{z_0}} \right)} \right] + \rho ({z_0}) \end{array} $

(10)

式中

$ \begin{array}{l} {w^I}_u\left( {{z_0}} \right) = \frac{\alpha }{2}{\psi ^I}_u\left( {{z_0}} \right)\sqrt {\frac{{{q^I}_u}}{{{p^I}_u\left( {{z_0}} \right)}}}, {w^J}_u\left( {{z_0}} \right) = \\ \frac{\beta }{2}{\psi ^J}_u\left( {{z_0}} \right)\sqrt {\frac{{{q^J}_u}}{{{p^J}_u\left( {{z_0}} \right)}}}, \\ {w^I}_u\left( {{z_0}} \right) = \alpha {\psi ^I}_u\left( {{z_0}} \right)\sqrt {{p^I}_u\left( {{z_0}} \right){q^I}_u} \times \\ ({({{\mathit{\hat \Sigma }}^I}u(z0))^{-1}}(\mathit{\boldsymbol{\tilde \mu }}_u^I-\mathit{\boldsymbol{\mu }}_u^I(z0))\\ {w^J}_u\left( {{z_0}} \right) = \beta {\psi ^J}_u\left( {{z_0}} \right)\sqrt {{p^J}_u\left( {{z_0}} \right){q^J}_u} \times \\ ({({{\mathit{\hat \Sigma }}^J}u(z0))^{-1}}(\mathit{\boldsymbol{\tilde \mu }}_u^I - \mathit{\boldsymbol{\mu }}_u^I(z0)) \end{array} $

求式 (10) 中$\rho $($z$) 关于$z$的导数有

$ \begin{array}{l} \frac{{\partial \rho \left( z \right)}}{{\partial z}} = \sum\limits_{i = 1}^n {{v^I}_i{k^\prime }} \left( {{{\left\| {\frac{{z-{x_i}}}{h}} \right\|}^2}} \right)\left( {z-{x_i}} \right)-\sum\limits_{u = 1}^m {{v^I}_u} + \\ \sum\limits_{i = 1}^n {{v^J}_i{k^\prime }} \left( {{{\left\| {\frac{{z - {y_i}}}{h}} \right\|}^2}} \right)\left( {z - {y_i}} \right) - \sum\limits_{u = 1}^m {{v^I}_u} \end{array} $

(11)

式中

$ \left\{ \begin{array}{l} {v^I}_i = \frac{C}{{{h^2}}}\sum\limits_{u = 1}^m {{w^I}_u} \left( {{z_0}} \right){\delta _{iu}}, {v^I}_u = \frac{{{w^I}_u\left( {{z_0}} \right)}}{h}\\ {v^J}_i = \frac{C}{{{h^2}}}\sum\limits_{u = 1}^m {{w^I}_u} \left( {{z_0}} \right){\delta _{iu}}, {v^J}_u = \frac{{{w^J}_u\left( {{z_0}} \right)}}{h} \end{array} \right. $

(12)

令$\partial \rho $($z$)=0，则可得到从当前位置$z$₀到新位置$z$₁的关系式

$ {z_1} = \frac{{\left[\begin{array}{l} \sum\limits_{i = 1}^n {\left[{{v^I}_ig\left( {{{\left\| {\frac{{{z_0}-{x_i}}}{h}} \right\|}^2}} \right){x_i} + } \right.} \\ \left. {{v^J}_ig\left( {{{\left\| {\frac{{{z_0}-{y_i}}}{h}} \right\|}^2}} \right){y_i}} \right] -\\ \sum\limits_{u = 1}^m {\left( {{v^I}_u + {v^J}_u} \right)} \end{array} \right]}}{{\left[\begin{array}{l} \sum\limits_{i = 1}^n {\left[{{v^I}_ig\left( {{{\left\| {\frac{{{z_0}-{x_i}}}{h}} \right\|}^2}} \right) + } \right.} \\ \left. {{v^J}_ig\left( {{{\left\| {\frac{{{z_0}-{y_i}}}{h}} \right\|}^2}} \right)} \right] \end{array} \right]}} $

(13)

式中，$g$($x$)=-$k$′($x$)。式 (13) 表明目标的位置是由红外图像与可见光图像共同决定。由于Epanechnikov核函数的导数为常量，因此式 (13) 可以进一步简化为

$ {z_1} = \frac{{\sum\limits_{i = 1}^n {\left( {{v^I}_i{x_i} + {v^J}_i{y_i}} \right)}-\sum\limits_{u = 1}^m {\left( {{v^I}_u + {v^J}_u} \right)} }}{{\sum\limits_{i = 1}^n {({v^I}_i + {v^J}_i)} }} $

(14)

如前所述，权重系数$\alpha $、$\beta $被用于调节可见光相似度与红外相似度在目标函数中所占的比重，很显然，它们的值应该是动态变化的。一般地，相邻帧间目标相似度不会发生较大变化，基于这一事实，本文通过前一帧的相似度值来决定当前帧的权重系数。假设前一帧中红外与可见光的最佳目标与其目标模型间的相似度分别为$\mathcal{B}[\mathit{\boldsymbol{I}}\left( z \right), \mathit{\boldsymbol{\tilde I}}]$和$\mathcal{B}[\mathit{\boldsymbol{J}}\left( z \right), \mathit{\boldsymbol{\tilde J}}]$，则当前帧的$\alpha $、$\beta $分别为

$ \left\{ \begin{array}{l} \alpha = B[\mathit{\boldsymbol{I}}\left( z \right), \mathit{\boldsymbol{\tilde I}}]/\left\{ {B[\mathit{\boldsymbol{I}}\left( z \right), \mathit{\boldsymbol{\tilde I}}] + B[\mathit{\boldsymbol{J}}\left( z \right), \mathit{\boldsymbol{\tilde J}}]} \right\}\\ \beta = B[\mathit{\boldsymbol{J}}\left( z \right), \mathit{\boldsymbol{\tilde J}}]/\left\{ {B[\mathit{\boldsymbol{I}}\left( z \right), \mathit{\boldsymbol{\tilde I}}] + B[\mathit{\boldsymbol{J}}\left( z \right), \mathit{\boldsymbol{\tilde J}}]} \right\} \end{array} \right. $

(15)

受光照、遮挡等因素的影响，目标表观可能会发生变化，因此必须及时更新其目标模型。在为数不多的研究核跟踪框架下模型更新的文献中，Peng等人^[16]提出的基于Kalman滤波的模型更新方法最有效，但Kalman滤波的高斯线性假设并不能很好地表述运动目标表观变化的高度随机性。因此，在借鉴其滤波器更新思想的基础上，本文拟采用粒子滤波方式对空间直方图模型进行在线更新。

称在当前帧$t$进行匹配时所采用的目标模型为“当前模型”，记为${{\mathit{\boldsymbol{\tilde h}}}^{t-1}} = {\left\{ {\tilde h_u^{t-1} = \left\langle {q_u^{t-1}, \mathit{\boldsymbol{\tilde \mu }}_u^{t - 1}, \mathit{\tilde \Sigma }_u^{t - 1}} \right\rangle } \right\}^m}_{u = 1}$，而由跟踪算法在最佳收敛点处得到目标图像的空间直方图称为“观测模型”，记为${\mathit{\boldsymbol{h}}^t} = {\left\{ {h_u^t = \left\langle {q_u^t, \mathit{\boldsymbol{\mu }}_u^t, \mathit{\Sigma }_u^t} \right\rangle } \right\}^m}_{u = 1}$。粒子滤波器在这两个模型中进行最优估计输出的模型称为“候选模型”，记为${{\mathit{\boldsymbol{\hat h}}}^t} = {\left\{ {\hat h_u^t = \left\langle {\hat q_u^t, \mathit{\boldsymbol{\hat \mu }}_u^t, \mathit{\hat \Sigma }_u^t} \right\rangle } \right\}^m}_{u = 1}$。由于直方图各分量是相互独立的，因此可分别用独立的滤波器对每个分量进行滤波更新，由此，定义空间直方图第$u$个分量的状态方程和观测方程为

$ \begin{array}{l} 状态方程:\tilde h_u^t = \tilde h_u^{t-1} + \omega _u^{t-1}\\ 观测方程:h_u^t = \tilde h_u^t + v_u^t \end{array} $

(16)

式中，${\omega ^{t-1}}_u$为状态噪声，${\nu ^t}_u$为观测噪声。由于协方差矩阵$\sum $是2×2的对角矩阵，因此只对其对角元素更新，从而使式 (16) 中状态变量的维数降成5维 (1$q$+2$\mu $+2$\sum $)。

在滤波的过程中，粒子的重要性是通过该粒子对应的观测与当前模型之间的欧氏距离来度量，即${{\tilde \pi }^k}_u = {\rm{exp}}\left\{ {-{{\left\| {\tilde h_u^k-\tilde h_u^{t-1}} \right\|}^2}/{\sigma ^2}} \right\}/\sqrt {2{\rm{\pi }}} \sigma $，式中${\tilde h_u^k}$为第$k$个粒子。${\omega ^{t-1}}_u$和${\nu ^t}_u$均以高斯白噪声模拟，虽然也是高斯线性假设，但与文献[16]不同，本文方法并不局限于高斯线性的假设。

为了避免模型过更新，通过计算第$t$帧时“观测模型”与“当前模型”之间的相似度来确定是否采纳目标“候选模型”${{\boldsymbol{\hat h}}^t}$。相似度计算公式为

$ \begin{array}{c} \mathcal{B}\left[{{\mathit{\boldsymbol{h}}^t}, {{\mathit{\boldsymbol{\tilde h}}}^{t-1}}} \right] = \\ \sum\limits_{u = 1}^m {{\psi ^t}_u} \left( {{\mathit{\boldsymbol{\mu }}^t}_u, {\mathit{\Sigma }^t}_u, \mathit{\boldsymbol{\tilde \mu }}_u^{t -1}, \mathit{\tilde \Sigma }_u^{t -1}} \right)\sqrt {{p^t}_u\tilde q_u^{t -1}} \end{array} $

(17)

如果$\mathcal{B}$很小说明本次观测模型受到了剧烈的干扰，不应当采用滤波结果更新模型，而仍然采用当前的模型，于是得到更新准则函数

$ {{\mathit{\boldsymbol{\tilde h}}}^t} = \left\{ {\begin{array}{*{20}{l}} {{{\mathit{\boldsymbol{\tilde h}}}^{t - 1}}} & {\mathcal{B}\left[ {{\mathit{\boldsymbol{h}}^t},{{\mathit{\boldsymbol{\tilde h}}}^{t - 1}}} \right] < \gamma }\\ {{{\mathit{\boldsymbol{\hat h}}}^t}} & {\mathcal{B}\left[ {{\mathit{\boldsymbol{h}}^t},{{\mathit{\boldsymbol{\tilde h}}}^{t - 1}}} \right] \ge \gamma } \end{array}} \right. $

(18)

式中，$\gamma $为阈值。实验表明取0.5就能很好满足需求。

算法1 基于粒子滤波的目标模型更新算法

输入: ${{\hat h}^{t-1}}$

1)  for $u$=1, ..., $m$ do

2) 粒子采样:根据状态转移方程生成采样粒子$\left\{ {\tilde h_u^k} \right\}_{k = 1}^K$;

3) 权值计算:根据观测方程得到每个粒子的观测，由此计算重要性权值$\left\{ {\tilde \pi _u^k} \right\}_{k = 1}^K$，并进行归一化;

4) 重采样:计算${N_{{\rm{eff}}}} = 1/\sum\limits_k {{{\left( {\tilde \pi _u^k} \right)}^2}} $，如果$N$_eff < $N$_th，则根据系统重采样方法得到新的粒子集$\left\{ {\bar h_u^k, \bar \pi _u^k} \right\}_{k = 1}^K$; 否则，$\left\{ {\bar h_u^k, \bar \pi _u^k} \right\}_{k = 1}^K = \left\{ {\tilde h_u^k, \tilde \pi _u^k} \right\}_{k = 1}^K$;

5) 状态估计: $\hat h_u^t = \sum\limits_k {\bar \pi _u^k\bar h_u^k} $;

6)  end for

输出:${{\hat h}^t} = \left\{ {\hat h_u^t} \right\}_{u = 1}^m$

算法1给出了第$t$帧时的目标模型更新流程。显然，滤波器的个数等于直方图分量数$m$，所以更新算法的时间复杂度是恒定的。而且，在实际应用中只处理非零分量，因此整个更新算法的计算量很少。

算法2 总结了基于联合空间直方图的跟踪 (JST) 流程。容易看出，当将空间相似度全部设为1时，JST就退化为使用普通直方图的联合跟踪算法，叫其基于联合直方图的跟踪 (JHT)。

算法2 基于联合空间直方图的跟踪算法

输入:初始的目标位置$z$₁、目标模型${\tilde I}$₁和${\tilde J}$₁；

1)  for $t$=2, ..., $T$ do

2) 根据式 (1)—式 (3) 计算当前位置候选目标的$\boldsymbol{I}$($z$_$t$－1) 和$\boldsymbol{J}$($z$_$t$－1)，并根据式 (9) 计算$\rho $($z$_$t$－1)；

3) 根据式 (12) 计算权值${\left\{ {{v^I}_i, {v^J}_i} \right\}^n}_{i = 1}, {\left\{ {{v^I}_u, {v^J}_u} \right\}^m}_{u = 1}$；

4) 根据式 (14) 找到新的候选目标位置$z_t$；

5) 根据式 (1)—式 (3) 计算$\boldsymbol{I}\left( {{z_t}} \right), \boldsymbol{J}({z_t})$和$\rho ({z_t})$；

6)  while $\rho \left( {{z_t}} \right) < \rho ({z_{t-1}})$ do

7) 将$\left( {{z_{t-1}} + {z_t}} \right)/2 \Rightarrow {z_t}$，并计算$\boldsymbol{I}\left( {{z_t}} \right), \boldsymbol{J}({z_t})$和$\rho ({z_t})$；

8)  end while

9) 如果$\left\| {{x_t}-{z_{t-1}}} \right\| < \in $则转向步骤10)，否则${z_t} \Rightarrow {z_{t-1}}$，转向步骤3)；

10) 计算$\mathcal{B}[\mathit{\boldsymbol{I}}\left( {{z_t}} \right), {{\mathit{\boldsymbol{\tilde I}}}_{t-1}}], \mathcal{B}[\mathit{\boldsymbol{J}}\left( {{z_t}} \right), {{\mathit{\boldsymbol{\tilde J}}}_{t-1}}]$，并更新$\alpha $和$\beta $；

11) 如果$\mathcal{B}[\mathit{\boldsymbol{I}}\left( {{z_t}} \right), {{\mathit{\boldsymbol{\tilde I}}}_{t-1}}] < \gamma $，则将${{\mathit{\boldsymbol{\tilde I}}}_{t-1}} \Rightarrow {{\mathit{\boldsymbol{\tilde I}}}_t}$，否则${{\mathit{\boldsymbol{\tilde I}}}_{t-1}} \Rightarrow {{\mathit{\boldsymbol{\tilde h}}}^{t-1}}$，启动算法1得到${{\mathit{\boldsymbol{\hat h}}}^t}$，并将${{\mathit{\boldsymbol{\hat h}}}^t} \Rightarrow {{\mathit{\boldsymbol{\tilde I}}}_t}$；

12) 如果$\mathcal{B}[\mathit{\boldsymbol{J}}\left( {{z_t}} \right), {{\mathit{\boldsymbol{\tilde J}}}_{t-1}}] < \gamma $，则将${{\mathit{\boldsymbol{\tilde J}}}_{t-1}} \Rightarrow {{\mathit{\boldsymbol{\tilde J}}}_t}$，否则${{\mathit{\boldsymbol{\tilde J}}}_{t-1}} \Rightarrow {{\mathit{\boldsymbol{\tilde h}}}^{t-1}}$，启动算法1得到${{\mathit{\boldsymbol{\hat h}}}^t}$，并将${{\mathit{\boldsymbol{\hat h}}}^t} \Rightarrow {{\mathit{\boldsymbol{\tilde J}}}_t}$；

13)  end for

算法2主要涉及式 (1)—式 (3) 的空间直方图计算、式 (5)(6) 的相似度计算以及式 (14) 的中心位置计算，而权值调整和模型更新部分的计算量相对较小。设式 (1)(2)(3) 的平均计算量分别为${C_p}、{C_\mu }和{C_\varepsilon }$，则空间直方图计算量为$m({C_p} + {C_\mu } + {C_\varepsilon })$。设${\psi _u}$的计算量为${C_\varphi }$，则式 (5) 的相似度计算量为$m({C_\varphi } + 2{C_m} + {C_a} + {C_s})$，式中${C_m}、{C_a}、{C_s}$分别为乘法、加法和开方运算量。${\boldsymbol{w}^I}_u$的计算量为$(2{C_m} + {C_s}), {\boldsymbol{w}^I}_u$的计算量为$(7{C_m} + {C_s} + 4{C_a} + {C_d})$，式中${C_d}$为矩阵求逆运算，${\boldsymbol{w}^J}_u$和${\boldsymbol{w}^J}_u$的计算与${\boldsymbol{w}^J}_u$和${\boldsymbol{w}^J}_u$类似，则式 (14) 计算量约为$[\left( {4mn + 4n + 14m} \right){C_m} + 4m{C_s} + \left( {2n + 8m} \right){C_a} + 2m{C_d}]$。算法2第2) 步要计算两个空间直方图和两个相似度，步骤3)—9) 的一次迭代要计算一个中心位置、4个空间直方图和4个相似度。假设算法在每一帧中的平均迭代次数为$M$，其复杂度可近似地估计为${C_{ms}} \approx (2 + 4M)m[({C_p} + {C_\mu } + {C_\varepsilon }) + $$({C_\varphi } + 2{C_m} + {C_a} + {C_s})] + M \times [(4mn + 4n + 14m){C_m} + $$4m{C_s} + (2n + 8m){C_a} + 2m{C_d}]$。

如果采用粒子滤波的方法实现目标搜索，则虽不需要计算位移公式，但需要进行2$N$个空间直方图和相似度计算，其中$N$为粒子个数，即其计算量约为${C_{pf}} \approx 2Nm[({C_p} + {C_\mu } + {C_\varepsilon }) + ({C_\varphi } + 2{C_m} + {C_a} + {C_s})]$。通常，${C_a}:{C_m}:{C_d}:{C_s}:{C_p}:{C_\mu }:{C_\varepsilon }:{C_\varphi }$=1:10:70:400:60:60:600: 500，$M$=5，$n$=30×30，$m$=16×16×16。当粒子数$N$=300时，粒子滤波搜索和均值漂移搜索能获得同样的跟踪精度，此时前者的时间约是后者的${C_{pf}}:{C_{ms}} \approx 4.4$倍。采用Matlab和MEX对算法2进行混合编程，在仿真环境为主频2×3.2 GHz、内存4 GB的PC机上对JST进行了多次测试，结果表明对帧频为30帧/s的视频数据，算法可同时在线跟踪到4个目标。

采用4段视频来测试JST的性能，并将其与₁跟踪器 (L1T)^[5]、基于区域模糊动态融合 (FRD) 的跟踪器^[7]，以及JHT这3个算法作了比较和分析。所有跟踪器的初始目标相同，均由手工在第1帧中标出。初始化$\alpha = \beta = 0.5$。

定量比较中使用6个评价指标，即位置误差$\varepsilon = \sqrt {{{({x_{\rm{G}}}-{x_{\rm{T}}})}^2} + {{({y_{\rm{G}}}-{y_{\rm{T}}})}^2}} $、平均位置误差$\bar \varepsilon = \sum\limits_i {{\varepsilon _i}/n} $、重叠率$\in = area({\boldsymbol{R}_{\rm{G}}} \cap {\boldsymbol{R}_{\rm{T}}})/area({\boldsymbol{R}_{\rm{G}}} \cup {\boldsymbol{R}_{\rm{T}}})$、平均重叠率$\bar \in = \sum\limits_i {{ \in _i}/n} $、成功率$sr = \sum\limits_i {\delta \left( {{ \in _i} > 0.5} \right)/n} $以及平均计算时间${\bar t}$，式中，$n$为总帧数，$({x_{\rm{G}}}, {y_{\rm{G}}}, {\boldsymbol{R}_{\rm{G}}})$分别为真实目标的中心和边界盒，$({x_{\rm{T}}}, {y_{\rm{T}}}, {\boldsymbol{R}_{\rm{T}}})$分别为跟踪器给出的目标中心和边界盒。一个理想的跟踪器应该有较小的位置误差、较快的跟踪速度、较高的重叠率和成功率。JST、JHT和FRD均为融合跟踪方法，都是通过对红外与可见光数据源进行联合来完成跟踪的，因此得到的跟踪结果都只有一个，相应的每种性能指标也都只有一个。而L1T为非融合跟踪方法，它是对红外和可见光视频单独跟踪，因此每种性能指标都有两个，为了区别，用L1TIR表示跟踪红外目标，而L1TVS表示跟踪可见光目标。表 1列举了定量比较结果。不难发现，相比而言，本文所提出的方法在位置误差、重叠率、成功率以及计算时间指标上获得了较好的综合性能。从图 1和图 2也能看出，JST方法的逐帧位置误差和重叠率整体较好。

1) 夜间跟踪。实验1通过跟踪video1序列对来测试JST在夜晚跟踪目标的性能，此序列中，可见光目标被淹没在夜幕中。该序列对有红外、可见光图像各270帧，帧大小为320×240像素，其跟踪结果如图 3所示，其中，第1行为可见光序列，第2行为红外序列，下同。图 1(a)和图 2(a)分别展示了各算法在该序列对上的逐帧位置误差和重叠率。可以看到JST、FRD、JHT这些融合跟踪方法自始至终都能较好地跟踪到运动目标，而L1T跟踪器对红外目标跟踪比较稳定，但是对可见光目标跟踪却非常糟糕，因为在夜间，可见光目标基本掩埋在黑暗中，而热红外传感不受光照影响。这也充分说明，通过红外与可见光的联合，可以弥补可见光在夜间成像的不足，从而拓展跟踪器的时间生存空间。

2) 背景阴影。实验2通过跟踪video2序列对来测试当目标所处背景出现阴影时JST的性能，此序列中，可见光目标的整体颜色与背景阴影非常相近。该序列对有红外、可见光图像各216帧，帧大小为320×240像素。图 4、图 1(b)和图 2(b)分别展示了各算法在该序列对上的部分跟踪结果、逐帧位置误差和重叠率。可以看到目标颜色与背景阴影相近的情况和目标处于夜间的情况比较类似，本质上都是目标与背景间的对比度不高，但两种情况对热红外成像来说没有影响。因此，与video1相似，JST、FRD、JHT这些融合跟踪方法自始至终都能较好地跟踪到运动目标，而L1T跟踪器对红外目标跟踪比较稳定，但是对可见光目标跟踪却非常糟糕。

3) 目标交汇与拥簇。实验3通过跟踪video3序列对来测试当背景拥簇时JST的性能，此序列中出现了多个热目标的交汇拥簇。该序列对有红外、可见光图像各300帧，帧大小为640×480像素。图 5、图 1(c)和图 2(c)分别展示了各算法在该序列对上的部分跟踪结果、逐帧位置误差和重叠率。可以看出JST自始至终都成功地跟踪到了目标对象，这主要是因为该方法将多模融合、目标空间结构保持和表观更新进行了有机整合。与JST相比，JHT的直方图模型中没有保留目标的空间结构信息，因此当出现热目标交互时 (大约在第220帧) 就会被相似的图像区域吸引过去了。FRD方法采用先融合后检测的跟踪策略，由于本视频中出现了多个位置相近的目标，且这些目标在红外成像上外观极为相似，因此FRD检测到其他相似目标的概率就会大大增加。值得注意的是，在该实验中，L1T能很好地跟踪到可见光目标，却跟丢了红外目标，这主要是因为在可见光模式下，真实目标与其左边的同伴在颜色纹理上反差较大，而在红外模式下，两人区分度极低。当然，从第270帧开始，可见光目标也跟丢了，因为此时目标进入了一个拥簇的背景环境中。

4) 目标遮挡。实验4通过跟踪video4序列对来测试JST在处理目标遮挡时的性能，此序列中，目标对象被多次部分或完全遮挡。该序列对有红外、可见光图像各150帧，帧大小为288×240像素。图 6、图 1(d)和图 2(d)分别展示了各算法在该序列对上的部分跟踪结果、逐帧位置误差和重叠率。不难发现，在第80帧以前，所有跟踪方法都比较稳定的跟踪到了男子头部，但之后男子头部被树叶全部遮挡，导致所有跟踪器都失效了，这主要是因为这些跟踪器都缺乏遮挡处理机制。从红外图像中可以清晰的看到，再之后 (大约是130140帧之间) 男子头部又短暂地重现视场中，FRD方法又跟踪到了目标，这是因为该方法是通过全图搜索来锁定目标的，因此只要目标再次出现，它就能检测到。由于缺乏遮挡处理机制，因此本文方法目前还是无法处理目标遮挡情况。

论文以空间直方图为目标表示模型，并通过红外目标相似度和可见光目标相似度的加权联合，实现了核框架下的多源目标特征级融合跟踪。算法在目标空间结构保持、模型更新、权值自适应调节和实时性方面实现了有机统一，是将传统单核跟踪向多核融合跟踪的有效拓展。通过对多组红外与可见光序列对的测试结果表明，本文所设计的跟踪器除了能跟踪良好环境下的运动目标之外，在处理夜间环境、背景阴影以及背景拥簇方面也展现出了较好的性能。目前，本文跟踪器处理目标遮挡的能力较弱，因此下一步我们将深入研究这个问题。