Event-based Stereo Visual Odometry（双目事件相机里程计）论文学习

本文详细介绍一篇双目事件相机里程计的论文，Event-based Stereo Visual Odometry，港科大沈邵劼团队Yi Zhou和TU Berlin的Guillermo Gallego共同完成，并公布了代码。我准备在接下来一段儿时间认真研究一下这篇论文与代码，欢迎想深入研究这个方法的朋友与我交流。转载请注明出处

论文下载：arXiv, https://arxiv.org/abs/2007.15548
github代码：https://github.com/HKUST-Aerial-Robotics/ESVO
柯翰同学在公众号的一篇简单介绍：基于双目事件相机的视觉里程计

1. 算法流程

算法总体流程如上图所示，主要包括三个重要模块（虚线矩形框），分别是：数据处理、Tracking和Mapping。暂时不去管初始化部分，我们关注数据流的方向，可以发现，在Tracking时，传入的是Time-Surface（本文往后缩写为TS）和Local Map，得到运动变换位姿pose；之后，pose和TS等共同完成完整的Mapping过程。

2. 基础知识

2.1 Time-Surface 时间表面

是一种事件相机数据表达方式，请看之前我博客中对这个进行的总结：【事件相机整理】角点检测与跟踪总结 - Time surface / SAE

2.2 Student t 分布

wiki介绍：https://en.wikipedia.org/wiki/Student%27s_t-distribution
作者通过对数据集测试，发现左右目相事件机中，匹配像素的残差$r$统计上符合“学生t分布”（中文表述不清，原文为：$r$ denotes the temporal difference between two corresponding pixels $x_{i,1},x_{i,2}$ inside neighborhoods）由此实现了多帧图像对深度进行融合。

2.3 ZNCC 零均值归一化相似性

是双目像素匹配时利用的一种度量方法，在极线约束上寻找匹配的像素点。在算法初始化部分使用。

2.4 IRLS 迭代重权重最小二乘

wiki介绍: https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares
一种非线性迭代优化算法，在Mapping过程中利用IRLS得到最优逆深度

2.5 Compositional LK

LK光流跟踪方法的一种变形，增量式而不是绝对式，避免了重复计算完整的雅克比矩阵。在Tracking过程中不断更新计算Warpping Parameter，从而得到运动轨迹。注意每次只是计算更新的部分$W(x,\rho ;\theta )\leftarrow W(x,\rho ;\theta )\cdot W(x,\rho ;\Delta \theta )$，即式中的$\Delta \theta$

3. Tracking部分

3.1 TS negative

论文并没有使用原始的TS进行Tracking，而是使用了 TS negative（后面简写为TSn)，$\overline{\tau }(x)=1-\tau (x)$。作者指出，由于TS的一侧往往是斜坡状，另一侧是陡峭的，所以沿着斜坡可以搜寻到当前Edge在图像中的位置。由此，当使用$(1-\tau )$时，值越小意味着距离真实Edge越近。

3.2 目标函数

上面提到，当TSn越小，这意味着图像这个位置对应了真实世界的Edge。在Tracking时，我们需要根据已有的Mapping结果（Local Map），将已知逆深度的Edge点投影到Frame上，计算对应的TSn的数值，进行相加，越小这意味着投影的参数约正确。例如，当位姿完全正确时，真实世界的Edge投影到Frame时，对应的TSn应该全接近0。从而我们可以写出目标函数为：$${\theta }^{\ast }=\underset{\theta }{\mathrm{argmin}}\sum _{x\in S}(\overline{\tau }(W(x,\rho ;\theta ){)}^2,$$ 即寻找一个最佳的Warpping Parameter $\theta$，使Frame中所有有逆深度的点($x\in S$)经过Warp后对应的TSn尽可能小。其中Warp就是定义为反投影后经过6DoF变换再投影，$$W(x,\rho ;\theta )={\pi }_{left}(T({\pi }_{ref}^{-1}(x,\rho ),G(\theta ))),$$ 即根据参考位置反投影，经过运动参数$G(\theta )$变换再投影到左目。论文指出，只利用左目进行计算就够了，加上右目没有明显的提升。

3.3 优化过程

论文指出，利用Compositional LK算法，不断迭代更新，得到最佳的T。

4. Mapping

Mapping部分主要也有两部分，首先对event的逆深度进行估计，之后再进行半稠密重建。

4.1 Event的深度估计

4.1.1 目标函数

深度估计时，需要知道Tracking的轨迹，进行运动补偿（所以需要一个初始化部分，才能够保证Tracking部分能够首先进行），估计一个event的逆深度$\rho$。

当我们估计在$t$时刻一个event的深度$\rho$时，我们可以得到这个event在左右目的坐标$x_1^t,x_2^t$，如果深度估计的正确，这两者对应的TS应该是相同，且之前一段儿时间$\delta t$(100Hz)内对应的TS应该也是相同的。所以目标函数定义为： $${\rho }^{\ast }=\underset{\rho }{\mathrm{argmin}}C(x,\rho ,{\mathrm{T}}_{left}(,t),{\mathrm{T}}_{right}(,t),T_{t-\delta t,t})$$ 其中${\mathrm{T}}_{left/right}$为该像素点在这段儿时间内的TS对应的值，而$C$具体为：$$C(...)=\sum _{x_{1,i}\in W_1,x_{2,i}\in W_2}|{\tau }_{left}^t(x_{1,i})-{\tau }_{right}^t(x_{2,i}){|}_2^2$$ 即对应的TS值的差。而$x_1,x_2$为：$$x_1=\pi (T_{c_{t-ϵ}}^{c_t}\cdot {\pi }^{-1}(x,{\rho }_k),$$ 即经过当前假设逆深度投影后，经过对应时间戳投影到指定时刻的Frame。

4.1.2 优化

优化采用任何一种非线性优化方法即可，例如GN方法

4.1.3 初始化

初始化时，由于不知道 $T$ 的轨迹，所以直接在极线上搜索匹配，只用ZNCC即可。不知道$T$所以无法得到 $T_{c_{t-ϵ}}$，所以不是很精确，但作为初始化已经足够了。

4.2 半稠密重建

由4.1，我们得到了一个稀疏的深度图，下一步通过融合的方式，得到半稠密的深度图。在这个过程中，并不是利用一次观测，而是利用了多次观测的结果，所以同时收敛了深度的取值。这部分内容类似深度滤波，但提出了新的融合方法。

如图，每次观测得到一个稀疏深度图，经过一段儿时间后（M次观测）融合得到更为精确、半稠密的深度图。

作者通过真值，统计两个数据集估计逆深度的残差分布，发现是个Student t分布（下文缩写为St）。当某个像素有多次观测时，多个逆深度进行融合，得到新的St分布。

接下来用到了IRLS进行优化，但我还没能理解这部分是在干啥，就先跳过了。

在融合时，将每次观测计算得到的有深度的像素点投射到下一次观测的像平面，由于不是整数像素，所以在新的像平面中考虑周围4像素是否是同一个数据点。1. 如果是同一个数据点，则利用St分布叠加方式将上一次的逆深度分布作为先验，叠加最新的观测逆深度分布获得后验逆深度分布；2. 如果不是同一个数据点（由二者逆深度分布的均值和方差判断），则取方差较小的作为当前像素点的逆深度；3. 如果最新的测量中没有得到这个像素点的逆深度，则直接定为上次的分布。由此，实现了多次测量的逆深度的融合，也得到了更稠密的深度图。

5. 小结

本文我认为最有创新之处就是多次测量的融合方法。

其他很多地方作者考虑的很细致，也有很多技巧。参考文献有五六十篇非常之多，只能说港科和Guillermo Gallego在EB和Stereo VO上有很深的积累。

(Running results of source codes)