1、跟踪基础知识简介
首先要说明一点，现在多目标跟踪算法的效果，与目标检测的结果息息相关，因为主流的多目标跟踪算法都是TBD（Tracking-by-Detecton）策略，SORT同样使用的是TBD，也就是说先检测，再跟踪。这也是跟踪领域的主流方法。所以，检测器的好坏将决定跟踪的效果。
本文抛开目标检测（YOLO V3）不谈，主要看SORT的跟踪思路。SORT采用的是在线跟踪的方式，不使用未来帧的信息。在保持100fps以上的帧率的同时，也获得了较高的MOTA（在当时16年的结果中）。
1.1、SORT的3个主要贡献
- 利用强大的CNN检测器的检测结果来进行多目标跟踪；
- 使用基于卡尔曼滤波（Kalman filter）与匈牙利算法（Hungarian algorithm）的方法来进行跟踪；
- 开源了代码，为MOT领域提供一个新的baseline。
1.2、SORT算法的简单理解
- 跟踪之前，对所有目标已经完成检测；
- 第一帧进来时，以检测到的目标初始化并创建新的跟踪器，标注id；
- 后面帧进来时，先到卡尔曼滤波器（Kalman Filter）中得到由前面帧box产生的状态预测和协方差预测。求跟踪器所有目标状态与本帧检测的Box的IOU，通过匈牙利算法（Hungarian Algorithm），得到IOU最大的唯一匹配（数据关联部分），在去掉匹配值小于IOU_threshold的匹配对；
- 用本帧中匹配到的目标检测Box去更新卡尔曼跟踪器，计算卡尔曼增益，状态更新和协方差更新。并将状态更新值输出，作为本帧的跟踪Box。对于本帧中没有匹配到的目标重新初始化跟踪器，卡尔曼跟踪器联合了历史跟踪记录，调节历史Box与本帧Box的残差，更好地匹配跟踪id。
- 其中卡尔曼滤波（Kalman filter）与匈牙利算法（Hungarian algorithm）对于大家来说可能是两个新名词。先简单解释一下，匈牙利算法是一种寻找二分图的最大匹配的算法，在多目标跟踪问题中可以简单理解为寻找前后两帧的若干目标的匹配最优解的一种算法。而卡尔曼滤波可以看作是一种运动模型，用来对目标的轨迹进行预测，并且使用确信度较高的跟踪结果进行预测结果的修正，是控制领域常用的一种算法。
1.3、卡尔曼滤波
1.3.1、卡尔曼滤波算法有两个基本假设
(1)、信息过程的足够精确的模型，是由白噪声所激发的线性(也可以是时变的) 动态系统;
(2)、每次的测量信号都包含着附加的白噪声分量。
当满足以上假设时，可以应用卡尔曼滤波算法。
1.3.2、问题描述与定义
定义一个随机离散时间过程的状态向量 ，该过程用一个离散随机差分方程描述：
其中n维向量 为k时刻的系统状态变量，n维向量是k-1时刻的系统状态变量。A是状态转移矩阵或者过程增益矩阵，是n×n阶方阵，它将k-1时刻状态和当前的k时刻状态联系起来。B是可选的控制输入 u∈Rl的增益，在大多数实际情况下并没有控制增益，所以这一项很愉快的变成零了。是n维向量，代表过程激励噪声，它对应了中每个分量的噪声，是期望为0，协方差为Q的高斯白噪声，。
再定义一个观测变量 ，得到观测方程：
其中观测值是m阶向量，状态变量是n阶向量。H是m×n阶矩阵，代表状态变量对测量变量的增益。观测噪声是期望为0，协方差为R的高斯白噪声，~N(0,R)。
1.3.3、卡尔曼滤波算法步骤
卡尔曼滤波器可以分为时间更新方程和测量更新方程。时间更新方程（即预测阶段）根据前一时刻的状态估计值推算当前时刻的状态变量先验估计值和误差协方差先验估计值；测量更新方程（即更新阶段）负责将先验估计和新的测量变量结合起来构造改进的后验估计。时间更新方程和测量更新方程也被称为预测方程和校正方程。因此卡尔曼算法是一个递归的预测——校正方法。
(1)、预测(Prediction)：根据上一时刻(k-1时刻)的后验估计值来估计当前时刻(k时刻)的状态，得到k时刻的先验估计值，卡尔曼滤波器时间更新方程：
(2)、 更新(Update)：使用当前时刻的测量值来更正预测阶段估计值， 得到当前时刻的后验估计值，卡尔曼滤波 器状态更新方程：
下面来一个个详细剖析每个参数：
1、和 ：分别表示k-1时刻和k时刻的后验状态估计值，是滤波的结果之一，即更新后的结果，也叫最优估计(估计的状态，根据理论，我们不可能知道每时刻状态的确切结果所以叫估计)。
2、：k时刻的先验状态估计值，是滤波的中间计算结果，即根据上一时刻(k-1时刻)的最优估计预测的k时 刻的结果，是预测方程的结果。
3、和：分别表示k-1时刻和k时刻的后验估计协方差(即和的协方差，表示状态的不确定度)，是滤波的结果之一。
4、:k时刻的先验估计协方差(的协方差)，是滤波的中间计算结果。
5、H：是状态变量到测量(观测)的转换矩阵，表示将状态和观测连接起来的关系，卡尔曼滤波里为线性关系， 它负责将m维的测量值转换到n维，使之符合状态变量的数学形式，是滤波的前提条件之一。
6、：测量值(观测值)，是滤波的输入。
7、：滤波增益矩阵，是滤波的中间计算结果，卡尔曼增益，或卡尔曼系数。
8、A：状态转移矩阵，实际上是对目标状态转换的一种猜想模型。例如在机动目标跟踪中，状态转移矩阵常常用来对目标的运动建模，其模型可能为匀速直线运动或者匀加速运动。当状态转移矩阵不符合目标的状态转换模型时，滤波会很快发散。
9、Q：过程激励噪声协方差(系统过程的协方差)。该参数被用来表示状态转换矩阵与实际过程之间的误差。因为我们无法直接观测到过程信号，所以Q的取值是很难确定的。一般有两种思路：一是在某些稳定的过程可以假定它是固定的矩阵，通过寻找最优的Q值使滤波器获得更好的性能，这是调整滤波器参数的主要手段，Q一般是对角 阵，且对角线上的值很小，便于快速收敛；二是在自适应卡尔曼滤波（AKF）中Q矩阵是随时间变化的。
10、R：表示测量噪声协方差，它是一个数值，这是和仪器相关的一个特性，作为已知条件输入滤波器。需要注意的是这个值过大过小都会使滤波效果变差，且R取值越小收敛越快，所以可以通过实验手段寻找合适的R值再利用它进行真实的滤波。
11、B：是将输入转换为状态的矩阵。
12、：实际观测和预测观测的残差，和卡尔曼增益一起修正先验(预测)，得到后验。
1.4、卡尔曼滤波的推导
- 更详细文末见文档，文档笔记中已经进行推导
1.5、卡尔曼滤波波代码注释
import numpy as np import scipy.linalg chi2inv95 = { 1: 3.8415, 2: 5.9915, 3: 7.8147, 4: 9.4877, 5: 11.070, 6: 12.592, 7: 14.067, 8: 15.507, 9: 16.919 } class KalmanFilter(object): def __init__(self): ndim, dt = 4, 1. # 构建Kalman Filter的矩阵model,并初始化初始状态 self._motion_mat = np.eye(2 * ndim, 2 * ndim) for i in range(ndim): self._motion_mat[i, ndim + i] = dt self._update_mat = np.eye(ndim, 2 * ndim) self._std_weight_position = 1. / 20 self._std_weight_velocity = 1. / 160 def initiate(self, measurement): mean_pos = measurement mean_vel = np.zeros_like(mean_pos) mean = np.r_[mean_pos, mean_vel] std = [ 2 * self._std_weight_position * measurement[3], 2 * self._std_weight_position * measurement[3], 1e-2, 2 * self._std_weight_position * measurement[3], 10 * self._std_weight_velocity * measurement[3], 10 * self._std_weight_velocity * measurement[3], 1e-5, 10 * self._std_weight_velocity * measurement[3] ] covariance = np.diag(np.square(std)) return mean, covariance def predict(self, mean, covariance): # 相当于得到t时刻估计值 # Q 预测过程中噪声协方差 std_pos = [ self._std_weight_position * mean[3], self._std_weight_position * mean[3], 1e-2, self._std_weight_position * mean[3] ] std_vel = [ self._std_weight_velocity * mean[3], self._std_weight_velocity * mean[3], 1e-5, self._std_weight_velocity * mean[3] ] # np.r_ 按列连接两个矩阵 # 初始化噪声矩阵Q motion_cov = np.diag(np.square(np.r_[std_pos, std_vel])) # x' = Fx mean = np.dot(self._motion_mat, mean) # P' = FPF^T+Q covariance = np.linalg.multi_dot((self._motion_mat, covariance, self._motion_mat.T)) + motion_cov return mean, covariance def project(self, mean, covariance): # 将状态分布投影到测量空间 std = [ self._std_weight_position * mean[3], self._std_weight_position * mean[3], 1e-1, self._std_weight_position * mean[3] ] # 初始化噪声矩阵R innovation_cov = np.diag(np.square(std)) # 将均值向量映射到检测空间，即Hx' mean = np.dot(self._update_mat, mean) # 将协方差矩阵映射到检测空间，即HP'H^T covariance = np.linalg.multi_dot((self._update_mat, covariance, self._update_mat.T)) return mean, covariance + innovation_cov def update(self, mean, covariance, measurement): # 计算 Kalman Filter的Correction Step==>Update Step # 将mean和covariance映射到检测空间，得到Hx'和S projected_mean, projected_cov = self.project(mean, covariance) # 矩阵分解 chol_factor, lower = scipy.linalg.cho_factor(projected_cov, lower=True, check_finite=False) # 计算卡尔曼增益K kalman_gain = scipy.linalg.cho_solve((chol_factor, lower), np.dot(covariance, self._update_mat.T).T, check_finite=False).T # z - Hx' innovation = measurement - projected_mean # x = x' + Ky new_mean = mean + np.dot(innovation, kalman_gain.T) # P = (I - KH)P' new_covariance = covariance - np.linalg.multi_dot((kalman_gain, projected_cov, kalman_gain.T)) return new_mean, new_covariance def gating_distance(self, mean, covariance, measurements, only_position=False): # 计算系统的预测优化值 mean, covariance = self.project(mean, covariance) if only_position: mean, covariance = mean[:2], covariance[:2, :2] measurements = measurements[:, :2] cholesky_factor = np.linalg.cholesky(covariance) d = measurements - mean z = scipy.linalg.solve_triangular(cholesky_factor, d.T, lower=True, check_finite=False, overwrite_b=True) squared_maha = np.sum(z * z, axis=0) return squared_maha
1.6、匈牙利算法
匈牙利算法是一种在多项式时间内求解任务分配问题的组合优化算法，并推动了后来的原始对偶方法。美国数学家哈罗德·库恩于 1955 年提出该算法。此算法之所以被称作匈牙利算法，是因为算法很大一部分是基于以前匈牙利数学家 Dénes Kőnig 和 Jenő Egerváry 的工作之上创建起来的。
说明一下匈牙利算法的流程：
以上图为例，假设左边的四张图是我们在第N帧检测到的目标(U)，右边四张图是我们在第N+1帧检测到的目标(V)。红线连起来的图，是我们的算法认为是同一行人可能性较大的目标。由于算法并不是绝对理想的，因此并不一定会保证每张图都有一对一的匹配，一对二甚至一对多，再甚至多对多的情况都时有发生。这时我们怎么获得最终的一对一跟踪结果呢？我们来看匈牙利算法是怎么做的。
第一步：
首先给左1进行匹配，发现第一个与其相连的右1还未匹配，将其配对，连上一条蓝线。
第二步：
接着匹配左2，发现与其相连的第一个目标右2还未匹配，将其配对。
第三步：
接下来是左3，发现最优先的目标右1已经匹配完成了，怎么办呢？我们给之前右1的匹配对象左1分配另一 个对象(黄色表示这条边被临时拆掉)
可以与左1匹配的第二个目标是右2，但右2也已经有了匹配对象，怎么办呢？我们再给之前右2的匹配对
象左2分配另一个对象（注意这个步骤和上面是一样的，这是一个递归的过程）。
此时发现左2还能匹配右3，那么之前的问题迎刃而解了，回溯回去。左2对右3，左1对右2，左3对右1。
所以第三步最后的结果就是：
第四步：最后是左4，很遗憾，按照第三步的节奏我们没法给左4腾出来一个匹配对象，只能放弃对左4的匹配，匈牙利算法流程至此结束。蓝线就是我们最后的匹配结果。至此我们找到了这个二分图的一个最大匹配。
最终的结果是我们匹配出了三对目标，由于候选的匹配目标中包含了许多错误的匹配红线（边），所以匹配准确率并不高。可见匈牙利算法对红线连接的准确率要求很高，也就是要求我们运动模型、外观模型等部件必须进行较为精准的预测，或者预设较高的阈值，只将置信度较高的边才送入匈牙利算法进行匹配，这样才能得到较好的结果。
匈牙利算法的流程大家看到了，有一个很明显的问题相信大家也发现了，按这个思路找到的最大匹配往往不是我们心中的最优。匈牙利算法将每个匹配对象的地位视为相同，在这个前提下求解最大匹配。这个和我们研究的多目标跟踪问题有些不合，因为每个匹配对象不可能是同等地位的，总有一个真实目标是我们要找的最佳匹配，而这个真实目标应该拥有更高的权重，在此基础上匹配的结果才能更贴近真实情况。
