基于事件对比度最大化的双目视觉里程计方法

1.本发明涉及机器人视觉技术领域，具体涉及基于事件对比度最大化的双目视觉里程计方法。


背景技术：

2.vo(visual odometry，视觉里程计)是vslam(visual slam，视觉slam)系统中最核心的部分，它是通过相机在运动过程中拍摄的一系列连续的图像来解算出相机的位姿和运动轨迹。然而，只依赖于标准相机的vo仍存在着许多问题。例如：在纹理少的区域，特征点数量不足，vo无法工作；因为对于具有高动态范围或高速运动的物体，基于帧的相机无法获得清晰的图像，视觉前端提取特征点困难；在动态场景中，移动的物体和行人会导致特征点匹配错误；因此，在这种情况下，基于帧的相机很难实现特征跟踪和相机位姿估计。
3.此外，基于帧的相机在静态场景中捕获冗余信息，这不仅会造成存储资源的浪费，而且在处理过程中也会消耗大量额外的计算资源。仿生事件相机，例如动态视觉传感器(dvs)，克服了基于帧的相机的上述限制。作为一种新型的生物启发式视觉传感器，事件相机具有与传统相机完全不同的模式，事件相机仅输出依赖像素亮度的变化。当每个像素的强度变化达到一定阈值时，将触发事件，事件携带像素坐标、时间戳和极性等信息。而imu(inertial measurement unit，惯性测量单元)则可以很好地捕捉高速运动时相机在三个坐标轴上的加速度和旋转信息，辅助相机的定位；同时，相机的位姿信息也有助于减少imu漂移的影响；此外，imu具有集成度高和轻巧、耐用、低价的优势。因此，将相机信息和imu信息融合起来，构成vio(visual-inertial odometry，视觉-惯性里程计)，可以使二者的优缺点互补，提高了系统的精准度和可靠性。同时，融合相机与imu信息的vio也极大地促进了vislam(visual-inertial slam)在救援机器人、ar/vr快速3d重建、自动驾驶等方面的应用。
4.目前，主流的vio框架都是基于单目相机和imu的融合，但受限于单目相机的尺度不确定性，该算法在恒定加速度运动条件下会失效，导致视觉-惯性里程计的精度和鲁棒性较差。此外，虽然单目相机位姿的绝对尺度可通过与imu信息对齐来获取，但imu通常会受到噪声的干扰，且地面机器人采用的往往是廉价的imu，对齐的结果较差且不够稳定。


技术实现要素：

5.本发明实施例的目的在于提供了基于事件对比度最大化的双目视觉里程计方法，解决了单目基于帧的相机与imu结合的情况下不能稳定运行的问题。
6.本发明的技术方案如下：
7.为实现上述目的，本发明提供一种基于事件对比度最大化的双目视觉里程计方法，包括以下步骤：
8.基于事件对比度最大化的双目视觉里程计方法，包括：
9.s1：对来自标准相机的图片进行预处理，建立事件对比度最大化的双目视觉里程
计模型；
10.s2：对来自事件相机的事件流和图片进行时间戳上的同步，根据事件的对应周期选择时空窗口；
11.s3：对每一个时空窗口中的事件进行对比度最大化，并求模板对应事件之间的光流；
12.s4：利用imu数据和计算出的光流去校正和更新模板边的位置；
13.s5：利用beta-gaussian的滤波器对模板边的深度进行估计，并得到3d空间中的模板边的位置；
14.s6：通过icp位姿求解算法去估计事件相机的运动轨迹。
15.进一步地，上述步骤s2具体如下：
16.假设在时间t0处从图像帧中检测到特征点x，则特征点的运动可以描述为：
[0017][0018][0019]
公式(1)中，为表示特征点x在时间s的速度；从事件流中选取一组事件，选取特征点x的时间作为初始时间t0对应的时空窗口；
[0020]
公式(2)中，w为时空窗内的事件集；ei表示时空窗口中的第i个事件；表示事件发生的时间；n表示事件的数量；由于[t0，t1]是第一个子时间间隔，可以通过设置时空窗口的大小来计算t1的值；
[0021]
为了实现特征点跟踪方法的异步性，子时间间隔的大小由该方法在实时运行时确定，具体计算过程如下：
[0022]
得到本次迭代中所有特征点的光流后，通过光流计算下一个子时间间隔的大小；定义xi来表示第i个特征点i＝{1，
…
，m}，特征点的个数为m；是在第n个子时间间隔[t
n-1
，tn]内特征点xi的光流；给定子时间间隔[t
n-1
，tn]内所有特征点的光流可以计算出下一个子时间间隔[tn，t
n 1
]；
[0023][0024]
在公式(3)中，数字3的单位是像素；表示子时间间隔[t
n-1
，tn]内所有特征点的平均光流；通过公式(3)计算出t
n 1
，将上一个时间间隔内的特征点平均移动3个像素所需的时间作为当前间隔的估计值。
[0025]
进一步地，上述步骤s3具体如下：
[0026]
通过对应特征点确定时空窗口后，使用对比度最大化算法将特征点x周围的事件集w与模板点集进行匹配；假设所有模板点与特征点x具有相同的光流，区域内所有像素的光流θ相同，并且特征点的光流在子时间间隔[t0，t1]内是恒定的；将时间间隔[t0，t1]内特
征点x的光流定义为v；对于w中的事件ci；使用扭曲事件的图像(iwe)计算其在时间t0的位置x’k
，公式如下：
[0027][0028][0029]
公式(5)中，x’kj
是第k个事件沿第j组光流扭曲后的位置；ne表示事件的个数；δ表示狄拉克函数；p
kj
表示第k个事件属于第j组光流的概率；ij(x)表示光流对应的iwe；
[0030]
事件通过图像对比度对齐，图像的对比度由锐度比色散度量定义，比如方差
[0031]
var(ij)＝∫
ω
(ij(x)-μj)2dx
ꢀꢀ
(6)
[0032]
公式(6)中，ω是图像平面；μj是扭曲事件图像的平均值；
[0033][0034]
公式(7)中，μ表示步长；n
l
表示簇数；θ为对应时空窗口中的光流。
[0035]
进一步地，上述步骤s4具体如下：
[0036]
得到特征点的光流后，更新特征点和模板边缘；然而，当相机旋转时，模板边缘和模板点以明显不同的速度移动，并且远离旋转中心的点移动得更快；更新模板边缘的过程分为两个步骤；使用光流更新模板边缘的位置，然后使用imu数据校正模板边缘的位置；
[0037]
光流θ用于更新特征点x的位置和模板边缘xj的对应位置；假设在子时间区间[t0，t1]内，特征点x的光流是恒定的，所以利用光流来更新特征点的位置，公式如下：
[0038]
x(t1)＝x(t0) v.(t
1-t0)
ꢀꢀ
(8)
[0039]
公式(8)中，x(t0)是最初提取的特征点的位置；x(t1)是更新后的特征点的位置；
[0040]
为了消除旋转对模板边缘更新的影响，引入imu数据来校正模板边缘位置；校正后的位置更接近模板边缘位置的真实值；在t0时刻计算模板点相对于特征点的相对位置，然后利用imu数据校正相对位置；对于模板点，xi我们定义它的相对位置：
[0041][0042]
定义符号xj和x分别表示相机坐标系中模板点xj和特征点x对应的3d坐标；
[0043][0044][0045]
公式(10)和(11)中，r和t是相机的旋转矩阵和平移向量，k表示相机的内参矩阵；将上述公式代入(9)并对3d坐标进行归一化；
[0046][0047]
xj在t1时刻的位置是通过特征点的位置和当时的相对位置相加得到的：
[0048][0049]
最后，通过公式即可完成所有模板点的更新。
[0050]
进一步地，上述步骤s5具体如下：
[0051]
得到每一帧边缘图的对应位置后，通过不同边缘图之间边缘上像素点的匹配关系，可以恢复出边缘在空间中的3d位置坐标；采用三角测量法，利用不同视角下观察到的3d点的像素位置，恢复3d点在空间中的坐标；同一个3d点可以在多帧中观察到，任意两帧都可以计算出对应的空间3d点坐标；在对同一3d点的最终估计中，该策略基于均匀的高斯混合分布；
[0052]
首先，使用三角测量方法来恢复边缘图中像素的深度；假设边图的某个像素点p0和边图的像素点p1是一对匹配点，两个像素点对应空间中的同一个空间3d点，公式(14)成立如下：
[0053]
y0k-1
p0＝y1rk-1
p1 t
ꢀꢀ
(14)
[0054]
其中，r和t分别是事件相机从边缘图到边缘图的旋转矩阵和平移向量，y0和y1分别表示事件相机坐标系中3d点在时间ti和时间tj的深度；公式(14)进一步简化为：
[0055][0056]
然后可以从两幅图像中的对应点计算出一个深度值yi；y的分布可以用高斯分布和均匀分布联合表示；
[0057][0058]
公式中是以真实值为中心的高斯分布，是其方差，π是内部概率，良好跟踪的特征的π接近1，y
min
和y
max
是均匀分布的下限和上限的特征。
[0059]
进一步地，上述步骤s6具体如下：
[0060]
给定场景的稀疏地图，通过icp算法获得当前相机位姿(r，t)；可以通过步骤s4和步骤s5计算对应的3d点集p二{p1，p2，
…
，pn}和q＝{q1，q2，
…
，qn}，其中对应点对的数量为n；然后通过icp算法计算出最优坐标变换，即旋转矩阵r和平移向量t，问题可以描述如下：
[0061][0062]
公式(17)中wi表示每个点的权重；r和t是我们需要的旋转矩阵和平移向量；为了对异常值具有鲁棒性，使用钟形tukey权重函数，如下所示：
[0063][0064]
公式(18)中且b＝5；此外，将长时间跟踪良好的特征被赋予最高权重，而丢失跟踪的特征通常由于较大的重投影误差而被删除。
[0065]
本发明的实施例提供的技术方案可以包括以下有益效果：
[0066]
1)本发明与传统的逐事件跟踪方法相比，提出了一种对比度最大化算法来解决事件和图像的数据关联，大大提高了事件流的计算速度。
[0067]
2)本发明由于对比度最大化算法高度依赖于场景的深度，因此提出了一种鲁棒的beta-gaussian分布深度滤波器，以获得比仅使用三角测量的深度估计更准确的线段模板深度。
[0068]
3)本发明应用于公共事件相机数据集的评估实验，与orb的视觉里程计算法相比，所提出的算法可以实现更好的性能并获得更低延迟的相机轨迹。
附图说明
[0069]
图1是本发明的流程图；
[0070]
图2(a)是本发明在图像帧上进行角点检测；
[0071]
图2(b)是本发明在图像帧上进行边缘检测；
[0072]
图3是本发明选择的第一个时间窗口的事件累积图；
[0073]
图4是本发明线段模板的更新示意图；
[0074]
图5是本发明方差收敛曲线图；
[0075]
图6是本发明三角化深度估计的示意图；
[0076]
图7是本发明和不同方法进行位姿估计后的轨迹与真实轨迹的对比图；
[0077]
图8是本发明与不同方法的绝对位姿误差对比图。
具体实施方式
[0078]
下面结合附图和具体实施方式对本发明进行详细说明。
[0079]
如图1所示，基于事件对比度最大化的双目视觉里程计方法，具体如下：
[0080]
步骤1，对来自标准相机的图片进行预处理，建立事件对比度最大化的双目视觉里程计模型；在本发明中我们使用了davis事件相机，davis由同一像素阵列中的基于事件的动态视觉传感器(dvs)和基于帧的有源像素传感器(aps)组成，考虑到场景边缘区域的事件比低纹理区域的事件更频繁地触发，我们设计了图像与事件结合的特征跟踪，去实现基于事件对比度最大化的双目视觉里程计，如图一所示；首先，从来自标准相机中的图像通过harris角点检测和canny边缘检测提取特征点和边缘图，然后在每个特征点周围选取指定区域内的边作为特征点的模板边，如图2(a)和图(b)所示。
[0081]
步骤2，对来自事件相机的事件流和图片进行时间戳上的同步，根据事件的对应周期选择时空窗口，如图3所示；假设在时间t0处从图像帧中检测到特征点x，则特征点的运动
可以描述为：
[0082][0083][0084]
为了实现特征点跟踪方法的异步性，子时间间隔的大小由该方法在实时运行时确定，具体计算过程如下：
[0085]
得到本次迭代中所有特征点的光流后，通过光流计算下一个子时间间隔的大小；定义xi来表示第i个特征点i＝{1，
…
，m}，特征点的个数为m；是在第n个子时间间隔[t
n-1
，tn]内特征点xi的光流；给定子时间间隔[t
n-1
，tn]内所有特征点的光流可以计算出下一个子时间间隔[tn，t
n 1
]；
[0086][0087]
步骤3，对每一个时空窗口中的事件进行对比度最大化，并求模板对应事件之间的光流，如图4所示；通过对应特征点确定时空窗口后，使用对比度最大化算法将特征点x周围的事件集w与模板点集进行匹配；对于w中的事件ei；使用扭曲事件的图像(iwe)计算其在时间t0的位置x’k
，公式如下：
[0088][0089][0090]
事件通过图像对比度对齐，图像的对比度由锐度比色散度量定义，比如方差
[0091]
var(ij)＝∫
ω
(ij(x)-μj)2dx
ꢀꢀ
(6)
[0092][0093]
θ为对应时空窗口中的光流。
[0094]
步骤4，利用imu数据和计算出的光流去校正和更新模板边的位置，如图5所示；得到特征点的光流后，更新特征点和模板边缘；更新模板边缘的过程分为两个步骤；使用光流更新模板边缘的位置，然后使用imu数据校正模板边缘的位置；
[0095]
光流θ用于更新特征点x的位置和模板边缘xj的对应位置；假设在子时间区间[t0，t1]内，特征点x的光流是恒定的，所以利用光流来更新特征点的位置，公式如下：
[0096]
x(t1)＝x(t0) v.(t
1-t0)
ꢀꢀ
(8)
[0097]
引入imu数据来校正模板边缘位置，校正后的位置更接近模板边缘位置的真实值；
在t0时刻计算模板点相对于特征点的相对位置，然后利用imu数据校正相对位置；对于模板点，xj我们定义它的相对位置：
[0098][0099]
定义符号xj和x分别表示相机坐标系中模板点xj和特征点x对应的3d坐标；
[0100][0101][0102]
公式(10)和(11)中，r和t是相机的旋转矩阵和平移向量，k表示相机的内参矩阵；将上述公式代入(9)并对3d坐标进行归一化；
[0103][0104]
xj在t1时刻的位置是通过特征点的位置和当时的相对位置相加得到的：
[0105][0106]
最后，通过公式即可完成所有模板点的更新。
[0107]
步骤5，利用beta-gaussian的滤波器对模板边的深度进行估计，并得到3d空间中的模板边的位置，如图6所示；得到每一帧边缘图的对应位置后，通过不同边缘图之间边缘上像素点的匹配关系，可以恢复出边缘在空间中的3d位置坐标；采用三角测量法，利用不同视角下观察到的3d点的像素位置，恢复3d点在空间中的坐标；在对同一3d点的最终估计中，该策略基于均匀的高斯混合分布；
[0108]
首先，使用三角测量方法来恢复边缘图中像素的深度；假设边图的某个像素点p0和边图的像素点p1是一对匹配点，两个像素点对应空间中的同一个空间3d点，公式(14)成立如下：
[0109]
y0k-1
p0＝y1rk-1
p1 t
ꢀꢀ
(14)
[0110]
公式(14)进一步简化为：
[0111][0112]
然后可以从两幅图像中的对应点计算出一个深度值yi；y的分布可以用高斯分布和均匀分布联合表示；
[0113][0114]
步骤6，通过icp位姿求解算法去估计事件相机的运动轨迹并画出他与0rb估计的相机轨迹的绝对位姿误差，如图7和图8所示所示。给定场景的稀疏地图，通过icp算法获得当前相机位姿(r，t)；可以通过步骤s4和步骤s5计算对应的3d点集p＝{p1，p2，
…
，pn}和q＝{q1，q2，
…
，qn}，其中对应点对的数量为n；然后通过icp算法计算出最优坐标变换，即旋转矩
阵r和平移向量t，问题可以描述如下：
[0115][0116]
公式(17)中wi表示每个点的权重；r和t是我们需要的旋转矩阵和平移向量；为了对异常值具有鲁棒性，使用钟形tukey权重函数，如下所示：
[0117][0118]
公式(18)中且b＝5；此外，将长时间跟踪良好的特征被赋予最高权重，而丢失跟踪的特征通常由于较大的重投影误差而被删除。
[0119]
最终事件相机的运动轨迹被icp算法求出。
[0120]
在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。此外，本领域的技术人员可以将本说明书中描述的不同实施例或示例进行接合和组合。
[0121]
本领域技术人员在考虑说明书及实践这里公开的发明后，将容易想到本发明的其它实施方案。本技术旨在涵盖本发明的任何变型、用途或者适应性变化，这些变型、用途或者适应性变化遵循本发明的一般性原理并包括本发明未公开的本技术领域中的公知常识或惯用技术手段。说明书和实施例仅被视为示例性的，本发明的真正范围和精神由所附的权利要求指出。