[论文阅读笔记13]Observation-Centric SORT(OCSORT)论文中的公式推导

这篇文章来自CVPR2022, 是我很喜欢的一篇文章, 尝试用更本质的方法优化Kalman滤波.

论文地址: 论文

1. 概述

这篇文章解决的主要问题是，现有的方法对运动预测都是基于线性运动假设，对非线性运动、遮挡、低帧率视频就没有好的处理效果。

作者说SORT有三个缺点：

1. 帧率很高的情况，目标位移的噪声可能就和位移大小本身差不多（因为位移会很小），这样Kalman的方差会很大。
2. 由于遮挡等原因，如果没有新的轨迹和检测匹配，那目标的噪声可能会积累。证明了误差积累关于时间是平方关系。
3. Kalman主要依赖于状态估计，把观测只作为辅助信息。作者认为现在的检测器足够强，应该已经可以把观测作为主要信息了（Observation-Centirc的来源）

针对以上三个问题，作者提出了三个创新：

1. 提出了observation为中心的平滑策略，来减少误差累积。
   具体地就是一个inactive的轨迹要Re-ID的时候，先为这个轨迹建立一个虚拟轨迹。虚拟轨迹的起始为它上一次被发现的位置，终止为它再次被发现的位置。在这两个位置之间作一个平滑。
2. 在cost matrix中加入了轨迹方向的一致性，证明了对相隔比较大的两点之间的方向进行估计可以减少噪声。
3. 为了减少在短时内目标因遮挡成为inactive的问题，作者提出了根据它们新就旧观测进行恢复的策略。

下面对论文里公式推导的细节做一些补充, 并分别介绍以上的三个部分.

1. Kalman滤波

关于Kalman滤波的背景和基础推导可参见Kalman. 下面用比较粗略的方式叙述.

对于 _线性, 高斯_的问题, 我们要对目标的状态进行估计. 我们已经有一个自己的(线性)模型, 以及有一个额外的力量可以对目标的某些性质进行观测(Observation). Kalman滤波要解决的问题是, 如何结合我们自己的模型和观测值, 来更好地对目标的状态进行估计.

简单进行一下推导. 从一维的情形出发, 假设状态变量为x ∈ R x\in\mathbb{R}x ∈R,我们有一个线性模型:

x t ^ = A x t − 1 ^ + B u t y t ^ = C x t ^ \hat{x_t}=A\hat{x_{t-1}}+Bu_t\ \hat{y_t}=C\hat{x_t}x t ​^​=A x t −1 ​^​+B u t ​y t ​^​=C x t ​^​
其中u t u_t u t ​为输入，A A A为状态转移系数（在高维就是矩阵），上标hat表示是模型的估计值，y t y_t y t ​表示输出。我们还有一个测量值，可以把测量值理解成真实模型的一种反映，假设也是现行的：
x t = A x t − 1 + B u t + ω t y t = C x t + v t x_t=Ax_{t-1}+Bu_t+\omega_t\ y_t=C x_t + v_t x t ​=A x t −1 ​+B u t ​+ωt ​y t ​=C x t ​+v t ​
其中ω t , v t \omega_t,v_t ωt ​,v t ​为均值为0的高斯噪声。因为噪声，因此我们假设线性模型估计的输出和测量的输出都是服从高斯分布的(暂时忽略时间下标t t t)：
y ^ ∼ N ( μ 1 , σ 1 2 ) y ∼ N ( μ 2 , σ 2 2 ) \hat{y}\sim N(\mu_1,\sigma_1^2)\ y\sim N(\mu_2,\sigma_2^2)y ^​∼N (μ1 ​,σ1 2 ​)y ∼N (μ2 ​,σ2 2 ​)
Kalman滤波通过将模型输出y , y ^ y,\hat{y}y ,y ^​ 相乘，以达到缩小方差的目的。设相乘后值为y ′ y’y ′, 根据高斯分布的性质:
μ y ′ = μ 1 + σ 1 2 ( μ 2 − μ 1 ) σ 1 2 + σ 2 2 σ y ′ 2 = σ 1 2 − σ 1 4 σ 1 2 + σ 2 2 \mu_{y’}=\mu_1+\frac{\sigma_1^2(\mu_2-\mu_1)}{\sigma_1^2+\sigma_2^2}\ \sigma_{y’}^2=\sigma_1^2-\frac{\sigma_1^4}{\sigma_1^2+\sigma_2^2}μy ′​=μ1 ​+σ1 2 ​+σ2 2 ​σ1 2 ​(μ2 ​−μ1 ​)​σy ′2 ​=σ1 2 ​−σ1 2 ​+σ2 2 ​σ1 4 ​​

在实际中，可能模型的输出与测量值的量纲不同。假设有一转换系数H H H, 使得H μ 1 H\mu_1 H μ1 ​和μ 2 \mu_2 μ2 ​的量纲相同，这样才有意义 （别忘了均值μ \mu μ 是我们关心的结果） 。因此对上式要加以修改：

H μ y ′ = H μ 1 + H 2 σ 1 2 ( μ 2 − H μ 1 ) H 2 σ 1 2 + σ 2 2 H 2 σ y ′ 2 = H 2 σ 1 2 − H 4 σ 1 4 H 2 σ 1 2 + σ 2 2 H\mu_{y’}=H\mu_1+\frac{H^2\sigma_1^2(\mu_2-H\mu_1)}{H^2\sigma_1^2+\sigma_2^2}\ H^2\sigma_{y’}^2=H^2\sigma_1^2-\frac{H^4\sigma_1^4}{H^2\sigma_1^2+\sigma_2^2}H μy ′​=H μ1 ​+H 2 σ1 2 ​+σ2 2 ​H 2 σ1 2 ​(μ2 ​−H μ1 ​)​H 2 σy ′2 ​=H 2 σ1 2 ​−H 2 σ1 2 ​+σ2 2 ​H 4 σ1 4 ​​
也即
μ y ′ = μ 1 + H σ 1 2 ( μ 2 − H μ 1 ) H 2 σ 1 2 + σ 2 2 σ y ′ 2 = σ 1 2 − H 2 σ 1 4 H 2 σ 1 2 + σ 2 2 \mu_{y’}=\mu_1+\frac{H\sigma_1^2(\mu_2-H\mu_1)}{H^2\sigma_1^2+\sigma_2^2}\ \sigma_{y’}^2=\sigma_1^2-\frac{H^2\sigma_1^4}{H^2\sigma_1^2+\sigma_2^2}μy ′​=μ1 ​+H 2 σ1 2 ​+σ2 2 ​H σ1 2 ​(μ2 ​−H μ1 ​)​σy ′2 ​=σ1 2 ​−H 2 σ1 2 ​+σ2 2 ​H 2 σ1 4 ​​

令
K = σ 1 2 H H 2 σ 1 2 + σ 2 2 K=\frac{\sigma_1^2H}{H^2\sigma_1^2+\sigma_2^2}K =H 2 σ1 2 ​+σ2 2 ​σ1 2 ​H ​
则
μ y ′ = μ 1 + K ( μ 2 − H μ 1 ) σ y ′ 2 = σ 1 2 − K H σ 1 2 \mu_{y’}=\mu_1+K(\mu_2-H\mu_1)\ \sigma_{y’}^2=\sigma_1^2-KH\sigma_1^2 μy ′​=μ1 ​+K (μ2 ​−H μ1 ​)σy ′2 ​=σ1 2 ​−K H σ1 2 ​

μ y ′ \mu_{y’}μy ′​就是Kalman对最终结果的预测值

将上式推广到高维，时变情形：,在第t − 1 t-1 t −1到t t t的状态转移中，线性模型的输出值为
x ^ t ∣ t − 1 = F t x ^ t − 1 ∣ t − 1 + B t u t \hat{x}{t|t-1}=F_t\hat{x}{t-1|t-1}+B_tu_t x ^t ∣t −1 ​=F t ​x ^t −1 ∣t −1 ​+B t ​u t ​

其中x ^ t − 1 ∣ t − 1 \hat{x}_{t-1|t-1}x ^t −1 ∣t −1 ​为上一时刻的Kalman输出值, u t u_t u t ​为输入.

假设x ^ t ∣ t − 1 \hat{x}{t|t-1}x ^t ∣t −1 ​的协方差矩阵为P t ∣ t − 1 P{t|t-1}P t ∣t −1 ​.

测量值（观测）值为z t z_t z t ​，转换矩阵为H t H_t H t ​, 则：

1. 一维的μ y ′ \mu_{y’}μy ′​就是Kalman最终的预测值，对应：
   x ^ t ∣ t = x ^ t ∣ t − 1 + K t ( z t − H t x ^ t ∣ t − 1 ) \hat{x}{t|t}=\hat{x}{t|t-1}+K_t(z_t-H_t\hat{x}_{t|t-1})x ^t ∣t ​=x ^t ∣t −1 ​+K t ​(z t ​−H t ​x ^t ∣t −1 ​)
2. 其中的Kalman矩阵：
   K t = P t ∣ t − 1 H t T H t T P t ∣ t − 1 H t + R t K_t=\frac{P_{t|t-1}H_t^T}{H_t^TP_{t|t-1}H_t+R_t}K t ​=H t T ​P t ∣t −1 ​H t ​+R t ​P t ∣t −1 ​H t T ​​
3. 一维的σ y ′ 2 \sigma_{y’}^2 σy ′2 ​对应Kalman预测值的协方差：
   P t ∣ t = P t ∣ t − 1 − K t H t P t ∣ t − 1 P_{t|t}=P_{t|t-1}-K_tH_tP_{t|t-1}P t ∣t ​=P t ∣t −1 ​−K t ​H t ​P t ∣t −1 ​

现在只差一个P t ∣ t − 1 P_{t|t-1}P t ∣t −1 ​没有求了。怎么求呢？

通过式x ^ t ∣ t − 1 = F t x ^ t − 1 ∣ t − 1 + B t u t \hat{x}{t|t-1}=F_t\hat{x}{t-1|t-1}+B_tu_t x ^t ∣t −1 ​=F t ​x ^t −1 ∣t −1 ​+B t ​u t ​入手。假如真实值为
x t ∣ t − 1 = F t x t − 1 ∣ t − 1 + B t u t + ω t x_{t|t-1}=F_tx_{t-1|t-1}+B_tu_t+\omega_t x t ∣t −1 ​=F t ​x t −1 ∣t −1 ​+B t ​u t ​+ωt ​，其中ω t ∼ N ( 0 , Q ) \omega_t\sim N(0,Q)ωt ​∼N (0 ,Q )为多维高斯噪声, Q为其自相关矩阵。

有：

P t ∣ t − 1 = E [ ( x t ∣ t − 1 − x ^ t ∣ t − 1 ) ( x t ∣ t − 1 − x ^ t ∣ t − 1 ) T ] = E [ ( F t ( x ^ t − 1 ∣ t − 1 − x t − 1 ∣ t − 1 ) + ω t ) ( ( x ^ t − 1 ∣ t − 1 − x t − 1 ∣ t − 1 ) T F t T + ω t T ) ] = F t E [ ( x ^ t − 1 ∣ t − 1 − x t − 1 ∣ t − 1 ) ( x ^ t − 1 ∣ t − 1 − x t − 1 ∣ t − 1 ) T ] F t T + E [ F t ( x ^ t − 1 ∣ t − 1 − x t − 1 ∣ t − 1 ) ω t T ] + E [ ω t ( x ^ t − 1 ∣ t − 1 − x t − 1 ∣ t − 1 ) T F t T ] + E [ ω t ω t T ] = F t P t − 1 ∣ t − 1 F t + Q P_{t|t-1}=E[(x_{t|t-1}-\hat{x}{t|t-1})(x{t|t-1}-\hat{x}{t|t-1})^T]\ =E[(F_t(\hat{x}{t-1|t-1}-x_{t-1|t-1})+\omega_t)((\hat{x}{t-1|t-1}-x{t-1|t-1})^TF_t^T+\omega_t^T)]\ =F_tE[(\hat{x}{t-1|t-1}-x{t-1|t-1})(\hat{x}{t-1|t-1}-x{t-1|t-1})^T]F_t^T\ +E[F_t(\hat{x}{t-1|t-1}-x{t-1|t-1})\omega_t^T]+E[\omega_t(\hat{x}{t-1|t-1}-x{t-1|t-1})^TF_t^T]\ +E[\omega_t\omega_t^T]\ =F_tP_{t-1|t-1}F_t+Q P t ∣t −1 ​=E [(x t ∣t −1 ​−x ^t ∣t −1 ​)(x t ∣t −1 ​−x ^t ∣t −1 ​)T ]=E [(F t ​(x ^t −1 ∣t −1 ​−x t −1 ∣t −1 ​)+ωt ​)((x ^t −1 ∣t −1 ​−x t −1 ∣t −1 ​)T F t T ​+ωt T ​)]=F t ​E [(x ^t −1 ∣t −1 ​−x t −1 ∣t −1 ​)(x ^t −1 ∣t −1 ​−x t −1 ∣t −1 ​)T ]F t T ​+E [F t ​(x ^t −1 ∣t −1 ​−x t −1 ∣t −1 ​)ωt T ​]+E [ωt ​(x ^t −1 ∣t −1 ​−x t −1 ∣t −1 ​)T F t T ​]+E [ωt ​ωt T ​]=F t ​P t −1 ∣t −1 ​F t ​+Q

1. SORT

SORT中的状态变量有7维, 即: x = [ u , v , s , r , u ˙ , v ˙ , s ˙ ] T x=[u, v, s, r, \dot{u}, \dot{v}, \dot{s}]^T x =[u ,v ,s ,r ,u ˙,v ˙,s ˙]T. 各部分含义如下:

u中心点x坐标v中心点y坐标s面积r长宽比(通常为常数)u一导x方向速度v一导y方向速度s一导面积变化率

SORT的线性模型是不含输入的，只有从过去状态到现在状态的估计。因此根据第2部分推导的Kalman滤波， 预测步骤为：

x ^ t ∣ t − 1 = F t x ^ t − 1 ∣ t − 1 P t ∣ t − 1 = F t P t − 1 ∣ t − 1 F t T + Q \hat{x}{t|t-1}=F_t\hat{x}{t-1|t-1}\ P_{t|t-1}=F_tP_{t-1|t-1}F_t^T+Q x ^t ∣t −1 ​=F t ​x ^t −1 ∣t −1 ​P t ∣t −1 ​=F t ​P t −1 ∣t −1 ​F t T ​+Q

更新步骤为：

K t = P t ∣ t − 1 H t T H t T P t ∣ t − 1 H t + R t x ^ t ∣ t = x ^ t ∣ t − 1 + K t ( z t − H t x ^ t ∣ t − 1 ) P t ∣ t = P t ∣ t − 1 − K t H t P t ∣ t − 1 K_t=\frac{P_{t|t-1}H_t^T}{H_t^TP_{t|t-1}H_t+R_t}\ \hat{x}{t|t}=\hat{x}{t|t-1}+K_t(z_t-H_t\hat{x}{t|t-1})\ P{t|t}=P_{t|t-1}-K_tH_tP_{t|t-1}K t ​=H t T ​P t ∣t −1 ​H t ​+R t ​P t ∣t −1 ​H t T ​​x ^t ∣t ​=x ^t ∣t −1 ​+K t ​(z t ​−H t ​x ^t ∣t −1 ​)P t ∣t ​=P t ∣t −1 ​−K t ​H t ​P t ∣t −1 ​

预测步骤与更新步骤循环进行，以上即为论文的(1)(2)式。

SORT采用线性模型, 即
u t + Δ t = u t + u ˙ Δ t u_{t+\Delta t}=u_t+\dot{u}\Delta t u t +Δt ​=u t ​+u ˙Δt
, v , s v,s v ,s同理, 注意一般认为宽高比a a a和各种速率u ˙ \dot{u}u ˙等不变, 因此线性模型对应的更新矩阵为

F t = [ 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 ] 7 × 7 F_t =\begin{bmatrix} 1 & 0 & 0 & 0& 1& 0& 0\ 0 & 1 & 0 & 0& 0& 1& 0\ 0 & 0 & 1 & 0& 0& 0& 1\ 0 & 0 & 0 & 1& 0& 0& 0\ 0 & 0 & 0 & 0& 1& 0& 0\ 0 & 0 & 0 & 0& 0& 1& 0\ 0 & 0 & 0 & 0& 0& 0& 1 \end{bmatrix}_{7\times 7}F t ​=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​1 0 0 0 0 0 0 ​0 1 0 0 0 0 0 ​0 0 1 0 0 0 0 ​0 0 0 1 0 0 0 ​1 0 0 0 1 0 0 ​0 1 0 0 0 1 0 ​0 0 1 0 0 0 1 ​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​7 ×7 ​

以上论文(3)式.

对于检测器的观测来说, 通常包含五种信息: 中心点坐标, 长宽和置信度, 记为z = [ u , v , w , h , c ] T z=[u,v,w,h,c]^T z =[u ,v ,w ,h ,c ]T.

随后文章证明了SORT的两个问题:1. 状态对噪声敏感 2. 没有观测的情况下, 噪声会随时间二次方累积. 下面简单对此证明.

3.1 状态对噪声敏感

假定状态x x x的u , v u,v u ,v分量也符合高斯分布（合理，因为x x x服从高斯分布），设u ∼ N ( μ u , σ u 2 ) u\sim N(\mu_u, \sigma_u^2)u ∼N (μu ​,σu 2 ​), 同理v ∼ N ( μ v , σ v 2 ) v\sim N(\mu_v, \sigma_v^2)v ∼N (μv ​,σv 2 ​). 因为是线性运动模型，我们可以这样估计速度：
u ˙ = u t + Δ t − u t Δ t v ˙ = v t + Δ t − v t Δ t \dot{u}=\frac{u_{t+\Delta t}-u_t}{\Delta t}\ \dot{v}=\frac{v_{t+\Delta t}-v_t}{\Delta t}u ˙=Δt u t +Δt ​−u t ​​v ˙=Δt v t +Δt ​−v t ​​

以上论文(5)式.

根据高斯分布相加的性质，有u ˙ ∼ N ( 0 , 2 σ u 2 / ( Δ t ) 2 ) , v ˙ ∼ N ( 0 , 2 σ v 2 / ( Δ t ) 2 ) \dot{u}\sim N(0, 2\sigma_u^2/(\Delta t)^2), \dot{v}\sim N(0, 2\sigma_v^2/(\Delta t)^2)u ˙∼N (0 ,2 σu 2 ​/(Δt )2 ),v ˙∼N (0 ,2 σv 2 ​/(Δt )2 ). 其中Δ t \Delta t Δt是帧数的差，我们考虑Δ t = 1 \Delta t=1 Δt =1, 也就是u ˙ ∼ N ( 0 , 2 σ u 2 ) , v ˙ ∼ N ( 0 , 2 σ v 2 ) \dot{u}\sim N(0, 2\sigma_u^2), \dot{v}\sim N(0, 2\sigma_v^2)u ˙∼N (0 ,2 σu 2 ​),v ˙∼N (0 ,2 σv 2 ​).

所以说，对速度的估计实际上把原有的方差放大了。但是通常没事，是因为即使线性模型预测的位置有问题导致方差变大，其实也可以通过detector来校正。

3.2 没有观测的情况下，误差累积速度为时间的平方

假设连续T T T帧没有观测，而且只依赖Kalman进行纯运动线索的更新。我们有：
u t + T = u t + u t ˙ T v t + T = u t + v t ˙ T u_{t+T}=u_t+\dot{u_t}T \ v_{t+T}=u_t+\dot{v_t}T u t +T ​=u t ​+u t ​˙​T v t +T ​=u t ​+v t ​˙​T

以上论文(6)式

注意，因为没有观测，因此我们只能利用其消失前的最后一帧的状态。

根据高斯分布相加的性质，有
u t + T ∼ N ( μ u , σ u 2 + ( T 2 ⋅ 2 σ u 2 ) ) v t + T ∼ N ( μ v , σ v 2 + ( T 2 ⋅ 2 σ v 2 ) ) u_{t+T} \sim N(\mu_u, \sigma_u^2+(T^2·2\sigma_u^2))\ v_{t+T} \sim N(\mu_v, \sigma_v^2+(T^2·2\sigma_v^2))u t +T ​∼N (μu ​,σu 2 ​+(T 2 ⋅2 σu 2 ​))v t +T ​∼N (μv ​,σv 2 ​+(T 2 ⋅2 σv 2 ​))

论文里这里应该是推错了

因此这时对位置估计的方差是时间差T T T的二次方

3.3 为什么可以主要依靠观测而不是线性模型

一句话，作者认为检测器的性能已经足够好，以至于可以保证检测器的方差小于线性模型预测结果的方差。

1. OCSORT的三方面改进

4.1 观测为中心的平滑(OOS)

该部分为了解决的问题是: 如果一个目标在untracked状态一段时间后要恢复, 如何减少误差的累积.

具体的做法是建立了一个虚拟轨迹, 轨迹的开始是目标最后一次的位置, 结束是目标再次被发现的位置. 以往的算法是直接根据当前位置预测什么的, 然而作者做的是结合当前的位置和最后一次被发现的位置做一个平滑. 我们假设估计的这个虚拟轨迹为最后一次位置, 再次发现位置和时间的函数:

z ^ t = T r a j v i r t u a l ( z t 1 , z t 2 , t ) , t 1 < t < t 2 \hat{z}t=Traj{virtual}(z_{t_1}, z_{t_2}, t),~~~t_1

论文(7)式

在实现中，作者也是用的恒定速度假设去建立这个虚轨迹的。根据匀速运动模型x ( t ) = x ( t 0 ) + v ( t − t 0 ) x(t)=x(t_0)+v(t-t_0)x (t )=x (t 0 ​)+v (t −t 0 ​)，有：
z ^ t = z t 1 + t − t 1 t 2 − t 1 ( z t 2 − z t 1 ) , t 1 < t < t 2 \hat{z}t=z{t_1}+\frac{t-t_1}{t_2-t_1}(z_{t_2}-z_{t_1}),~~~~~~~t_1

论文(10)式

代入Kalman的更新式子，我们有：
x ^ t ∣ t = x ^ t ∣ t − 1 + K t ( z ^ t − H t x ^ t ∣ t − 1 ) = F t x ^ t − 1 ∣ t − 1 + K t ( z ^ t − H t F t x ^ t − 1 ∣ t − 1 ) \hat{x}{t|t}=\hat{x}{t|t-1}+K_t(\hat{z}t-H_t\hat{x}{t|t-1})\ =F_t\hat{x}{t-1|t-1}+K_t(\hat{z}_t-H_tF_t\hat{x}{t-1|t-1})x ^t ∣t ​=x ^t ∣t −1 ​+K t ​(z ^t ​−H t ​x ^t ∣t −1 ​)=F t ​x ^t −1 ∣t −1 ​+K t ​(z ^t ​−H t ​F t ​x ^t −1 ∣t −1 ​)

这里有一个问题，在untracked目标被再次发现以后，通过建立虚拟轨迹来对再次发现后的状态做一个修正。文中说”Given the supervision of the virtual observation trajectory, error raised in state estimation would not be accumulated any more”, 但是没有给出证明.

4.2 以观测为中心的动量(OCM)

大多数MOT方法在建立cost matrix的时候基本依据两点: 运动特征和外观特征. 作者在相似度矩阵中还融入了一点, 那就是目标的运动趋势, 也就是动量.

因此总的cost matrix就写为:
C ( X ^ , Z ) = C i o u ( ( X ^ , Z ) + λ C v ( X ^ , Z , V ) C(\hat{X},Z)=C_{iou}((\hat{X},Z)+\lambda C_v(\hat{X},Z,V)C (X ^,Z )=C i o u ​((X ^,Z )+λC v ​(X ^,Z ,V )

其中X X X是状态, Z Z Z是观测, V V V是轨迹的运动方向, 实际上应该是按照θ = arctan ⁡ ( Δ y / Δ x ) \theta=\arctan(\Delta y/\Delta x)θ=arctan (Δy /Δx )得到的.

作者在附录证明了, 实际上估计方向时取的间距越大, 方差就会越小.

4.3 以观测为中心的恢复(OCR)

这个部分就是一个策略, 一句话:
一旦一条轨迹在正常关联阶段之后仍然没有被跟踪，我们尝试将 这条轨迹的最后一次观测与 新到来的时间步上的观测进行关联。

在伪代码中, 就是这部分:

Original: https://blog.csdn.net/wjpwjpwjp0831/article/details/124767905
Author: wjpwjpwjp0831
Title: [论文阅读笔记13]Observation-Centric SORT(OCSORT)论文中的公式推导