稀疏激光雷达与可见光/红外成像系统的标定方法

刘 宇，蔡 毅，戎宁涛，周云扬，王岭雪

〈系统与设计〉

稀疏激光雷达与可见光/红外成像系统的标定方法

刘 宇，蔡 毅，戎宁涛，周云扬，王岭雪

（北京理工大学 光电学院 光电成像技术与系统教育部重点实验室，北京 100081）

激光雷达与成像系统之间的位姿标定是激光点云与图像像素进行融合的前提。目前主流的离线标定方法中，普通棋盘格标定板用于64线及以上的激光雷达时效果较好，而用于16线激光雷达时由于其数据稀疏而导致误差较大。而且，涉及红外成像系统的标定时，需要特制的棋盘格来获得发射率差异。本文针对稀疏激光雷达点云数据较少的问题，研究了可以同时标定激光雷达与可见光、红外成像系统的方法，设计了菱形九孔标定板，并提出几何约束损失函数来优化特征点的坐标。最后，分别使用红外和可见光成像系统与16线激光雷达进行标定，实验结果表明，平均重投影误差均在3个像素之内，取得了较好的效果。本文方法还能用于稀疏激光雷达与可见光-红外多波段成像系统的标定。

稀疏激光雷达；位姿估计；标定；标定板；红外成像

0 引言

近年来，自动驾驶技术发展迅速，该技术将惯性导航系统、毫米波雷达、超声波雷达、激光雷达、成像系统、全球定位系统等部署在汽车上，配合车载高性能计算平台进行数据融合，在车辆行驶过程中自动完成场景理解、定位、避障，以保障车辆的安全行驶。

激光雷达根据激光测距原理，输出激光点云数据，可以得到目标的距离数据。常用成像系统根据工作波段以及成像原理的不同，分为可见光成像系统和红外成像系统。红外成像系统因可以全天时工作，也逐渐应用于自动驾驶，是对可见光图像的有益补充。激光雷达与可见光和红外成像系统融合，可获得周围环境的三维信息。这些信息在融合之前，需要先对成像系统进行内标定，在此基础上，再完成激光雷达与成像系统之间的外标定。内标定获得成像系统的内参矩阵，包括像素主点、像素焦距、畸变参数等。内参矩阵反映了图像像素与对应空间点的关系；外标定获得成像系统坐标系和激光雷达坐标系之间的刚体变换矩阵，包括旋转矩阵和平移向量等。刚体变换矩阵反映了成像系统与激光雷达之间的位姿关系。

激光雷达与成像系统之间的标定方法主要有两类：离线标定和在线标定。

离线标定又称有目标标定，是在非任务场景下，使用特定目标（通常是标定板），在两个系统中分别找到特定目标的对应特征点。早期，由于激光雷达线数较少，如4线、8线、16线，标定板设计通常比较简单。Rodriguez等[1]设计了挖孔同心圆标定板，特征点为同心圆的圆心；陈东等[2]设计了挖孔圆与棋盘格结合的标定板，特征点为棋盘格的角点。近年来，随着激光雷达技术的发展，激光雷达的线数也越来越多，机械扫描式激光雷达已经达到了128线，固态激光雷达也可以达到等效100线以上的水平。因此，标定任务开始转向使用非挖孔的棋盘格标定板来完成。使用棋盘格标定板的主要好处在于，在稠密的激光点云中可以方便地通过算法自动找到棋盘格平面的角点以及平面方程建立约束求解；棋盘格同时又可用于成像系统内参矩阵的标定。Zhou[3-4]和Geiger[5]等基于从激光点云中找到的棋盘格设计平面约束与直线约束获得刚体变换矩阵。

在线标定又称无目标标定，不需要使用特定目标，可以在执行任务的运动过程中完成。在线标定通常是基于手眼标定（hand-eye calibration）[6]完成的，即Ishikawa[7]等通过分别计算运动过程中激光雷达与成像系统的运动估计，求解方程。近年来，使用卷积神经网络（convolutional neural networks, CNN）进行在线标定成为研究热点，只需输入激光点云与对应的图像，网络即可求解得到刚体变换矩阵。Schneider[8]等提出了RegNet，该网络使用大量标定好的激光点云和图像数据进行训练，执行任务时直接回归计算得到刚体变换矩阵；Ganesh等[9]提出了自监督网络CalibNet，通过最大化点云与图像在几何和光度上的一致性得到刚体变换矩阵。

目前主流的棋盘格离线标定方法用于64线及以上的激光雷达效果较好，而用于16线激光雷达时由于其数据稀疏而导致误差较大，且由于棋盘格标定板是一个整体，不存在温度差，无法用于红外成像系统。

因此，本文设计一种标定方法以及标定板，可以满足16线稀疏激光雷达、红外成像系统与可见光成像系统的标定要求。16线激光雷达原理同样是激光测距，但点云信息较少，因此设计了挖孔标定板，且在挖孔特征点之间建立了几何约束；可见光成像系统成像原理是光的反射，为方便在可见光图像中分辨孔和标定板，将标定板整体设计为黑色；红外成像系统的成像原理是热辐射，挖孔标定板的孔和标定板存在温度差，可以在红外图像中分辨。

1 原理

1.1 标定原理

激光雷达与成像系统的标定实质上是坐标系之间的转换。本文使用(l,l,l)表示以激光雷达为中心的雷达坐标系3D点坐标；使用(c,c,c)表示成像系统坐标系下3D点坐标；使用(w,w,w)表示世界坐标系下3D点坐标；使用(,)表示图像的像素坐标系下2D点的坐标。相机的投影过程可以使用式(1)表示，激光雷达与成像系统的坐标转换可以使用式(2)表示：

式中：为成像系统的内参矩阵，可以使用张正友标定法[10]获得，wc和wc是世界坐标系到相机坐标系的旋转矩阵和平移向量，lc和lc是激光雷达坐标系到成像系统坐标系的旋转矩阵和平移向量，组合起来就是刚体变换矩阵。

通过找到对应特征点在激光雷达坐标系的坐标(l,l,l)以及像素坐标系的坐标(,)，就可以通过两种不同方法求解得到lc和lc。

第一种方法不必求解成像系统坐标系下坐标，直接使用2D和3D的对应特征点，来求解PnP（perspective N points）问题。将式(1)和式(2)组合，并使用齐次坐标，将内参和外参组合到一起，得到式(3)。这是一个有12个未知数的方程组，每一对点可以提供两个方程，因此最少需要6对特征点即可解出12个未知数，进而解出旋转矩阵和平移向量的初值。然后将二者作为初始解，使用L-M迭代法，对式(4)的代价函数迭代优化，得到lc和lc的迭代解，其中li与p是一对雷达和像素坐标对应点，≥6。

第二种方法需要求出特征点在成像系统坐标系下的坐标，使用3D和3D对应点来进行求解。根据张正友标定法[10]，在知道标定板的物理尺寸时，假定标定板所在平面的世界坐标为w＝0，就可以将式(1)中的投影过程转换为求解单应性矩阵（homography matrix）的过程，进而求解成像系统的内参矩阵。同理，也可以使用来分解单应性矩阵得到wc和wc，再根据式(1)，计算得到特征点在成像系统坐标系下坐标(c,c,c)。对于多对3D特征点，建立几何空间约束，即可求解方程组得到lc和lc的初值，然后使用ICP（iterative closest point）算法进行迭代求解。空间约束以及迭代方法根据标定板以及实验设备不同都有所不同，将在下一部分进行具体介绍。

1.2 本文3D约束方法

对于稀疏激光雷达，使用普通棋盘格标定板矩形边界点作为特征点，效果不够好，因此需要使用挖孔标定板的内部信息作为特征点。如图1所示，本文标定板由9个半径相等的圆孔组成，特征点对为激光雷达坐标系和图像像素坐标系中9个圆心的坐标。9个圆心之间有严格的几何约束关系，如平行约束、垂直约束、中点约束等；使用这些约束完成特征点坐标的优化计算以及lc和lc的初值计算。

图1 本文设计的标定板

对于图像中的特征点，通过Hough圆变换对图像中的圆进行检测，在每张图像中可以找到9个圆心的像素坐标(,)。

在激光雷达点云中，每条激光扫描线经过圆孔，会产生距离突变，如图2所示。假设有两条雷达扫描线经过该圆，这两条扫描各自会产生一段距离突变，突变的两个端点在圆上，即图中的、、、四点；又因圆上3点可确定一个圆，可将问题转换为求三角形外接圆的问题，使用、、三点即可求解圆心坐标。假设、、三点在激光雷达坐标系下坐标分别为(l1,l1,l1)、(l2,l2,l2)、(l3,l3,l3)，由、、三点确认的平面方程为＋＋＋＝0，圆心坐标为(l,l,l)，半径为（＞0）。可以得到由l、l、l、四个未知数组成的线性方程组，如式(5)所示。解该方程组，即可得到圆心坐标以及半径。为了使结果更准确，本文从、、、四点中任选三点，计算4次，结果求平均值作为圆心坐标。

图2 激光雷达数据计算圆心坐标

由于测量过程存在误差，在求得2D和3D的9个圆心特征点的坐标之后，还需利用标定板的几何约束对两个坐标系下特征点坐标进行优化。本文设计了损失函数（loss function），优化loss，使坐标的误差尽可能小。如图1，本文标定板存在中点约束，比如为中点等；垂直约束，比如^等。本文使用lossmid表示中点约束误差，如式(6)所示，即点与中点的距离；本文使用lossvertical表示垂直约束误差，如式(7)所示，即两个向量之间点积的绝对值；整体loss如式(8)所示，每个点均使用了4次，所占权重相等以免引入新的误差。本文使用随机搜索（random search）的方法，不断优化loss接近于0，降低9个点的误差。

分别对2D特征点和3D特征点坐标进行优化后，使用本文1.1节介绍的第二种方法来求lc和lc的初值，先将2D特征点转换为成像系统坐标系下3D坐标，然后与激光雷达坐标系下的3D特征点建立约束。与第一种求解PnP问题，直接得到外参矩阵的方法相比，第二种方法对旋转矩阵和平移向量是分别求解的，可以更好地减少二者的误差，因此效果更好。在第2章的仿真实验中，本文也将对比这两种方法的效果。首先建立约束来求解lc，相比于点坐标，空间向量只受到lc的约束，因此，使用直线的方向向量与平面的法向量来建立约束，如式(9)(10)所示：

对于两个坐标系下对应的9个特征点，求6条线的方向向量：、、CD、、、建立约束，分别使用、三点和、三点求解两次平面法向量建立约束；这样可以保证9个点都用了两次，所占权重一致。一个位置下的9组对应点可以建立8个约束，通过奇异值分解（singular value decom- position，SVD）的方法即可解出lc的初值，使用多个位置下的多组对应点，建立更多约束，可得到更精确的初值结果。

使用上面得到的lc的初值作为lc，计算lc的初值，使用的约束有平面约束、直线约束以及垂直约束，如式(11)(12)(13)所示：

对于一个位置下的9组对应点，求出了2个平面法向量和6个直线方向向量。使用9个点的质心和两个平面法向量代入式(11)得到2个方程；使用6条线段的中点和方向向量代入式(12)得到18个方程；使用雷达坐标系下、、、四个点在成像系统坐标系的投影点，分别代入式(13)，建立^、^、^、^四个垂直约束，得到4个方程，一共得到24个方程。因为平移向量lc有3个未知数，24个方程构成的是超定线性方程组，可以求得lc的线性最小二乘解作为初值。

在求得lc和lc的初值之后，就可以使用3D-2D点迭代或3D-3D点迭代的方法来求解lc和lc的优化解，本文激光雷达与成像系统标定方法的整体流程如下：

输入：多组图像与对应的激光雷达点云数据；

输出：激光雷达到成像系统的旋转矩阵lc和平移向量lc；

1） 使用张正友标定法[10]获得相机内参；

2） 每组数据在图像中和点云中分别找到对应的9个圆心点坐标；

3） 使用设计的loss优化点的坐标（2D和3D）；

4） 根据像素坐标得到相机坐标系下坐标；

5） 建立相机和雷达坐标系下的3D约束，求得lc和lc初值；

6） 使用3D-2D点迭代或3D-3D点迭代，得到lc和lc的优化解。

2 仿真

实际情况下，lc和lc的真实值是未知的，通常使用特征点的重投影误差来评价lc和lc的准确性。但这个指标默认了特征点没有误差，为了更好地验证本文方法的有效性，本文进行了仿真实验。

2.1 初值结果与PnP初值结果对比

本文1.1中介绍过，使用3D约束求解初值比PnP求解初值效果更好，本节将进行仿真实验验证结果。首先使用随机的尺寸生成标定板，并随机生成特征点的激光雷达坐标、lc和lc；使用固定的内参矩阵，得到成像系统坐标系和像素坐标系下坐标，分别使用3D约束和PnP方法求解lc和lc初值，与实际值进行对比，多次仿真结果如表1所示。lc使用旋转轴的单位向量和旋转角度表示，lc使用向量表示。表中表示使用3D约束的结果，表示使用PnP方法的初值结果，旋转轴误差和平移误差均为旋转轴向量与平移向量相对于实际的1误差，角误差为旋转角与实际旋转角的误差，单位是弧度。可以看到，使用3D约束方法求初值的误差相对于PnP方法的误差小了几个数量级，因此使用3D约束求初值可以大大减少后续迭代求解所需要的迭代次数，而PnP方法的优点在于简单。

2.2 本文3D约束方法的仿真

2.1节仿真了理想情况下求解初值的结果，验证了本文方法的有效性。在实际测量中，得到的特征点坐标存在一定的误差，需要对2.1节生成的激光雷达坐标系数据和成像系统坐标系数据随机增加误差，来模拟测量值，然后使用本文方法得到lc和lc的结果，仿真结果如表2所示。

表1 两种方法初值误差对比

表2 本文方法的仿真结果

从表2中可以看出，本文方法可以将特征点的loss大幅降低，以此方法解得的lc和lc误差与实际的误差非常小。仿真实验验证了本文方法的有效性，下一章将进行实际实验。

3 实验

3.1 实验设备

本次实验使用激光雷达分别与可见光成像系统、红外成像系统进行标定，以验证本文方法的有效性。

可见光成像系统使用LT-USB1080P型号的高清工业摄像头，成像分辨率为1920×1080像素，焦距为6mm；红外成像系统使用艾睿光电的IR-Pilot 640型号热像仪，成像分辨率为640×512像素，焦距6.9mm，工作波段为8～14mm；激光雷达使用Velodyne VLP-16型号的16线稀疏激光雷达，垂直视场角为－15°～15°，分辨率为2°；水平视场角为360°，分辨率最高为0.1°，两种实验系统如图3(a)(b)所示。根据前面的原理介绍，标定板的每个圆孔需要有至少两条激光扫描线穿过，结合激光雷达与成像系统的参数以及作用距离，设定标定板的尺寸为1.2m×1.35m，每个圆孔的半径为9cm。

3.2 实验过程

将标定板置于激光雷达和成像系统视场的重合区域，距离系统大概2m左右，使用实验系统采集多组数据，根据1.2节的本文的算法流程进行实验，图4(a)(b)展示了一组可见光图像和激光点云的对应数据。图4(c)(d)分别展示了红外实验系统的实验场景以及红外成像系统得到的红外图像。实验中使用Matlab的相机标定工具箱，获取成像系统的内参矩阵，使用Opencv中对应的函数来完成图像中的圆检测以及畸变校正，方法的其他部分使用Python进行编写。根据标定板的距离以及激光雷达的水平分辨率，可以计算得到，相邻激光扫描点之间的距离大约为0.5cm；也就是说理想情况下判断圆上点时的误差最大为0.5cm，考虑到一些激光点可能无法接收到返回值，以及激光雷达的误差和算法的误差，认为雷达坐标系中圆心点误差最大为5cm，像素坐标系中圆心点最大误差为10像素。根据这两个参数来优化loss，取多组数据建立3D约束求得初值，然后分别使用3D-3D的ICP算法和3D-2D的LM迭代法得到优化解。

3.3 实验结果

在可见光和红外两个波段系统的实验中，本文使用相同数量的特征点，迭代相同次数，得到最后的结果。

可见光系统以及红外系统的多次实验结果如表3、表4所示，其中表示使用ICP算法进行3D-3D迭代的结果，表示使用L-M迭代法进行3D-2D迭代的结果，均使用3D到2D特征点在和方向的平均重投影误差（单位为像素）来评价方法的准确性。

图3 标定实验装置

图4 实验场景及数据

表3 可见光系统实验结果

对于可见光系统，可以看到，本文方法在使用3D-3D迭代情况下，重投影误差平均在4个像素左右，而使用3D-2D迭代情况下，重投影误差平均小于3个像素，都取得了很好的效果。

对于红外系统，本文方法在使用3D-3D迭代情况下，重投影误差平均在4个像素以内，使用3D-2D迭代，重投影误差平均在3个像素以内，效果很好。

图5(a)(b)展示了两种迭代方法的特征点重投影误差散点图。可以看到，除了部分特征点之外，重投影误差都比较小，误差较大的点可能是因为激光雷达、图像算法等引入了测量误差。

图6(a)(b)分别展示了可见光系统和红外系统使用实验结果lc和lc将点云数据投影到图像上的融合结果，可以看到，除了在一些距离突变很大的边缘点有一点偏移之外，16条激光扫描线都很连贯，说明标定的效果很好，得到的融合数据可以用于后续的其他任务中。

表4 红外系统实验结果

图5 重投影误差

4 结论

本文针对稀疏激光雷达与可见光及红外成像系统之间的标定问题，设计了一种菱形九孔标定板，并根据标定板的约束，设计了几何约束loss来优化特征点的测量坐标，减小误差；根据3D几何约束得到lc和lc的初值，使用迭代法得到lc和lc的迭代解。

本文分别使用可见光成像系统、红外成像系统与16线激光雷达进行标定实验。在可见光系统中平均重投影误差在3个像素之内，在红外系统中平均重投影误差在3个像素之内。

图6 激光雷达点云投影结果

本文方法还适用于可见光-红外多波段成像系统，只需采集一组激光点云数据，以此为基准分别与多波段成像系统多个相机进行标定，方便了多波段的成像信息与激光雷达点云数据的融合过程。

[1] Fremont V, Bonnifait P. Extrinsic calibration between a multi-layer lidar and a camera[C]//, 2008: 214-219.

[2] 陈东. 三维激光和单目视觉间的联合标定与数据融合[D]. 大连: 大连理工大学, 2009.

CHEN Dong. Calibration and Data Fusion between 3D-Laser Scanner and Monocular Vision[D]. Dalian: Dalian University of Technology, 2009.

[3] ZHOU L, DENG Z. Extrinsic calibration of a camera and a lidar based on decoupling the rotation from the translation[C]//, 2012: 642-648..

[4] ZHOU L, LI Z, Kaess M. Automatic extrinsic calibration of a camera and a 3d lidar using line and plane correspondences[C]//(IROS), 2018: 5562-5569.

[5] Geiger A, Moosmann F, Car O , et al. Automatic camera and range sensor calibration using a single shot[C]//, 2012: 3936-3943.

[6] Horaud R, Dornaika F. Hand-eye calibration[J]., 1995, 14(3): 195-210.

[7] Ishikawa R, Oishi T, Ikeuchi K. Lidar and camera calibration using motions estimated by sensor fusion odometry[C]//(IROS), 2018: 7342-7349.

[8] Schneider N, Piewak F, Stiller C, et al. RegNet: multimodal sensor registration using deep neural networks[C]/(IV), 2017: 1803-1810.

[9] Iyer G, Ram R K, Murthy J K, et al. CalibNet: geometrically supervised extrinsic calibration using 3D spatial transformer networks[C]//(IROS), 2018: 1110-1117.

[10] ZHANG Z. A flexible new technique for camera calibration[J]., 2000, 22(11): 1330-1334.

Calibration Between Sparse LIDAR and Visible/Infrared Imaging Systems

LIU Yu，CAI Yi，RONG Ningtao，ZHOU Yunyang，WANG Lingxue

(,,,100081,)

Pose estimation between LIDAR and imaging system is the prerequisite for the data fusion. Among current mainstream off-line calibration methods,common checkerboard is generally effective for 64-line and above LIDAR, but not for 16-line LIDAR due to its sparse data and will lead to large error. Furthermore, when involving calibration of infrared imaging system, specially-made checkerboard is needed to produce difference of emissivity. Aiming at the problem of less information provided by sparse LIDARs, we propose a new calibration method that can jointly calibrate LIDAR and visible/infrared imaging systems. A novel diamond-shaped nine-hole calibration board is designed, and a geometric constraint loss function is proposed to optimize the coordinates of feature points. Finally, the infrared and visible light imaging systems are used respectively, to calibrate with 16-line LIDAR. Good results are achieved and show that, all the average reprojection error is within 3 pixels.The proposed method can also be used in calibration of multi-band imaging systems that include sparse LIDAR, visible imaging system and infrared imaging system.

sparse LIDAR, pose estimation, calibration, calibration board, infrared imaging

TP249

A

1001-8891(2022)12-1264-09

2022-03-12；

2022-04-19.

刘宇（1997-），男，辽宁沈阳人，硕士研究生，主要从事激光雷达与图像处理的研究。E-mail：andyliu0881@163.com。

王岭雪（1973-），女，云南石屏人，副教授，博士，主要从事红外成像、图像处理和红外光谱等方面的研究。E-mail：neobull@bit.edu.cn。