CN111260735A - 一种单次拍摄的lidar与全景相机的外参数标定方法 - Google Patents

一种单次拍摄的lidar与全景相机的外参数标定方法 Download PDF

Info

Publication number
CN111260735A
CN111260735A CN202010034949.9A CN202010034949A CN111260735A CN 111260735 A CN111260735 A CN 111260735A CN 202010034949 A CN202010034949 A CN 202010034949A CN 111260735 A CN111260735 A CN 111260735A
Authority
CN
China
Prior art keywords
point
chessboard
point cloud
lidar
points
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010034949.9A
Other languages
English (en)
Other versions
CN111260735B (zh
Inventor
何炳蔚
邓清康
胡誉生
林立雄
张立伟
陈彦杰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Fuzhou University
Original Assignee
Fuzhou University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Fuzhou University filed Critical Fuzhou University
Priority to CN202010034949.9A priority Critical patent/CN111260735B/zh
Publication of CN111260735A publication Critical patent/CN111260735A/zh
Application granted granted Critical
Publication of CN111260735B publication Critical patent/CN111260735B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/80—Analysis of captured images to determine intrinsic or extrinsic camera parameters, i.e. camera calibration
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/22—Matching criteria, e.g. proximity measures
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06V—IMAGE OR VIDEO RECOGNITION OR UNDERSTANDING
    • G06V10/00—Arrangements for image or video recognition or understanding
    • G06V10/70—Arrangements for image or video recognition or understanding using pattern recognition or machine learning
    • G06V10/74—Image or video pattern matching; Proximity measures in feature spaces
    • G06V10/75—Organisation of the matching processes, e.g. simultaneous or sequential comparisons of image or video features; Coarse-fine approaches, e.g. multi-scale approaches; using context analysis; Selection of dictionaries
    • G06V10/757—Matching configurations of points or features
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10028—Range image; Depth image; 3D point clouds
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/20—Special algorithmic details
    • G06T2207/20084—Artificial neural networks [ANN]

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Data Mining & Analysis (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Computing Systems (AREA)
  • Databases & Information Systems (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Multimedia (AREA)
  • Image Processing (AREA)

Abstract

本发明涉及一种单次拍摄的LIDAR与全景相机的外参数标定方法，包括如下内容将LIDAR与全景相机固定在Robotnik移动机器人上。然后将多个棋盘放置于LIDAR与全景相机的共同视场下，一次拍摄收集单帧的全景图像与该帧全景图像对应的点云数据；接着，利用生长的棋盘角点检测算法，检测出全景图像的棋盘角点；对点云数据进行预处理，分割去除点云地面，分割点云平面、提取棋盘点云；基于点云的反射强度，估计出棋盘点云的棋盘角点；最后，通过定义从棋盘左下侧开始的角点共同计数顺序，建立全景图像的棋盘角点与点云的棋盘角点的几何约束方程，求解出外部校准参数。只需要一次拍摄，就能实现LIDAR和全景相机的外参数标定。

Description

一种单次拍摄的LIDAR与全景相机的外参数标定方法
技术领域
本发明涉及传感器标定方法领域，特别是一种单次拍摄的LIDAR与全景相机的外参数标定方法。
背景技术
近些年来，随着机器人技术的发展成熟，广泛应用于资源勘探开发、救灾排险、家庭娱乐、定位导航等各类领域。为了使机器人在环境中感测更多可用信息，一般需要配备多类传感器。最常见的是LIDAR与全景相机的组合，全景相机能够获取丰富的颜色、形状、纹理等环境信息，但无法获得环境目标的距离信息；而LIDAR正好相反,它能够获取广泛范围的环境目标的位置信息与距离信息，但是无法得到颜色、形状、纹理等信息。因此，基于LIDAR与全景相机传感器的互补特性，将两类传感器数据融合，可以获得更为精确可用的目标信息。所以说，融合来自LIDAR和全景相机数据信息的关键步骤是准确、快捷的外部校准。
在移动机器人上搭载LIDAR和全景相机，LIDAR与全景相机传感器通过外部校准以便在公共坐标系中表示感测信息。为了确定3D LIDAR和全景相机两者之间的位置关系，通过建立LIDAR和全景相机各自采集目标特征之间的几何约束关系来求解。因此，外部校准方法可以分为基于点云特征线与图像特征线或者特征面、基于点云特征点与基于图像特征线或者特征面以及基于点云特征点与图像特征点的三种几何对应约束关系的方法。一般来说，基于点云特征点与图像特征点对应的几何约束关系的方法，比基于点云特征线与图像特征线或者特征面、基于点云特征点与基于图像特征线或者特征面的方法精度更高，但是在点云中，特征点相比于特征线与特征面更加难以捕捉。在大多数的激光校准工作中，经常需要手动干预校准过程，比如手动的选择点、线或者面，除此之外，实验中需要多次使用相机采集图像，激光采集点云数据，因此，实验过程繁琐。
发明内容
有鉴于此，本发明的目的是提供一种单次拍摄的LIDAR与全景相机的外参数标定方法，克服大多数技术需要手动干预，多次采集数据，以及校准过程繁琐不足的问题。
本发明采用以下方案实现：一种单次拍摄的LIDAR与全景相机的外参数标定方法，提供一Robotnik移动机器人，其特征在于：包括以下步骤：
步骤S1：将LIDAR(Velodyne-HDL-64e)与全景相机(Ladybug5)固定在Robotnik移动机器人上；然后将m个棋盘放置于LIDAR与全景相机的共同视场下，一次拍摄收集单帧的全景图像与该帧全景图像对应的点云数据；
步骤S2：利用生长的棋盘角点检测算法，检测出全景图像的棋盘角点Ic；
步骤S3：对点云数据进行预处理，分割去除点云地面，分割点云平面、提取棋盘点云；
步骤S4：基于点云的反射强度，估计出点云棋盘角点pL；
步骤S5：通过定义从棋盘左下侧开始的角点共同计数顺序，建立全景图像的棋盘角点Ic与点云的棋盘角点pL的几何约束方程，求解出外参数R*,t*。
进一步地，所述步骤S1具体包括以下步骤：
步骤S11：通过螺栓连接将LIDAR(Velodyne-HDL-64e)与全景相机(Ladybug5)固定在Robotnik移动机器人上；
步骤S12：构建一个室外场景，在场景中放置m块标定棋盘，每块棋盘的大小为600mm×450mm，棋盘中每个正方形的大小为75mm×75mm，并且满足多个棋盘在LIDAR与全景相机共同视场下的要求，其中，m取值为3、4、...11、12，m为整数；
步骤S13：利用步骤S11固定在移动机器人上的LIDAR与全景相机，使用全景相机收集一帧步骤S12构建场景的全景图像，LIDAR收集这帧全景图像对应的点云数据。
进一步地，步骤S2中所述检测出全景图像的棋盘角点Ic的具体内容为：
步骤S21：粗定位棋盘格角点的位置：首先定义两种不同类型的角点原型，原型1是一种和坐标轴平行的角点，原型2是一种与坐标轴成45°的角点；每个原型分别由4个卷积核组成，令原型1由四个卷积核K1,K2,K3,K4组成，原型2由四个卷积核K5,K6,K7,K8组成，分别用于与全景图像进行卷积操作；
通过两个角点原型来定义全景图像中每个像素点与角点的相似程度；
Figure BDA0002365136510000041
Figure BDA0002365136510000042
Figure BDA0002365136510000043
Figure BDA0002365136510000044
Figure BDA0002365136510000045
Figure BDA0002365136510000046
Figure BDA0002365136510000047
其中
Figure BDA0002365136510000048
表示与原型1相似程度的两种可能，
Figure BDA0002365136510000049
表示与原型2相似程度的两种可能，原型1与原型2的相似程度的两种可能相同，表示的是左对角线为黑，右对角线为白，或者左对角线为白，右对角线为黑，
Figure BDA00023651365100000410
表示卷积核K1,K2,K3,K4原型1在某个像素的卷积值，
Figure BDA00023651365100000411
表示卷积核K5,K6,K7,K8原型2在某个像素的卷积值，c表示图像棋盘角点的最大相似程度；通过计算角点相似程度，得到大致的角点范围；然后通过非极大值抑止算法来获得候选角点cp；
步骤S22：令c是理想的角点位置，p是c局部邻域的一个像素点，Gp是p点的图像梯度向量，此时满足如下式子：
Figure BDA00023651365100000412
由于实际图像中不止一个局部领域的像素点，因此在候选角点cp的邻域N(cp)内满足下列公式条件的就是所需要的棋盘角点Ic；
Figure BDA00023651365100000413
进一步地，所述步骤S3具体包括以下步骤：
步骤S31：估计点云棋盘角点之前，对点云数据进行预处理；通过PCL中的直通滤波器模块，将点云PcL＝{(x,y,z)}中X,Y方向超过8m远的点剔除；
Figure BDA0002365136510000051
其中，pi＝(x,y,z)是点云PcL中的一点；
步骤S32：根据步骤S31将点云远点剔除后，基于形态学算法将点云地面分离，用以减小在点云平面分割时，地面点的干扰；
Figure BDA0002365136510000052
Figure BDA0002365136510000053
其中，将点云测量点p(x,y,z)，x,y处的高度z定义为膨胀因子dp，腐蚀因子ep与膨胀因子相对应；w是测量点p的邻近窗口大小；
wj＝2jb+1 (13)
对于w窗口大小根据上述公式，线性增加窗口大小；其中j＝1,2,3,...,360，j为整数，b是初始窗口大小；
定义一个因子s来判断切除深度阈值dh；令物体坡度恒定，则最大坡度差为dhmax(t),k，因子s与窗口大小wk存在关系：
Figure BDA0002365136510000054
Figure BDA0002365136510000061
其中，dhT,j表示第j窗口切除深度阈值，将点云中的点依次带入公式(11)-(15)，计算出切除深度阈值dhT,j，若腐蚀因子大于切除深度阈值，则将该点移除，否则该点保存，其中T表示阈值，j表示第j个窗口大小；
步骤S33：根据步骤S32将点云地面移除后，剩余的点云集合为PrL；利用区域增长算法对点云平面进行分割；首先将剩余点云中每个点的曲率值
Figure BDA0002365136510000062
从小到大排序，曲率值越小，表示所在区域越平坦；然后，将曲率值最小点Pmin加入种子点，搜索它的K个最近邻点
Figure BDA0002365136510000063
计算每个点的法向量
Figure BDA0002365136510000064
并与最小法向量Nmin相比，若不大于平滑阈值Tth，则该点加入点云平面RL:
Figure BDA0002365136510000065
若小于曲率阈值cth，则将该点加入种子点，生长每个区域直到它收敛，并从PrL移除，如此循环，直到PrL中没有剩余的点；
最后使用棋盘平面度、形状大小条件，提取出棋盘点云
Figure BDA0002365136510000066
m表示棋盘的数量；
其中，提取棋盘点云
Figure BDA0002365136510000067
的具体内容：
由提取的点云平面RL组成矩阵Mn×3，沿着三个基矢量Mb＝(ux,uy,uz)T的方向分解，每个基矢量上的分量比为λ1,λ2,λ3；当最小比率λ3小于0.05并且单个棋盘的宽度dW和高度和dH满足公式(17)的点云片段被认为是棋盘点云；
Figure BDA0002365136510000071
其中W,H是单个棋盘的宽度和高度。
进一步地，所述步骤S4具体包括以下步骤：
步骤S41：将第k块棋盘点云
Figure BDA0002365136510000072
降低到二维平面并与第k块图像棋盘对应：使用主成分分析法，通过旋转矩阵
Figure BDA0002365136510000073
与平移矩阵
Figure BDA0002365136510000074
如公式(18)所示，将第k块棋盘点云
Figure BDA0002365136510000075
转化为与LIDAR坐标系一致的棋盘平面坐标系，其它棋盘点云跟随旋转平移变化；此过程中，求解出Cov矩阵的三个特征值(γ1,γ2,γ3)，得到对应的三个特征向量(v1,v2,v3)，旋转矩阵
Figure BDA0002365136510000076
定义为(v1,v2,v3)；其中，k表示第k块棋盘，取值范围是1-m；
Figure BDA0002365136510000077
Figure BDA0002365136510000078
Figure BDA0002365136510000079
步骤S42：第k块棋盘点云转化为与LIDAR坐标系一致的棋盘平面坐标系后，利用黑白棋盘格图案反射强度的对应关系，设置阈值[λL,λH]，小于λL表示从黑色图案反射的低强度，而大于λH表示从白色图案反射的高强度；
反射强度值在[λL,λH]为点云棋盘角点所在区域；通过黑白棋盘格反射强度的对应关系来制定成本函数，如公式(21)所示，从而估计出点云棋盘角点pL；
Figure BDA0002365136510000081
Figure BDA0002365136510000082
是
Figure BDA0002365136510000083
降维到XOY平面的第k块棋盘点云，其中，i表示第i点，coi表示
Figure BDA0002365136510000084
落入的图案颜色，令黑色为0，白色为1；ri是点云反射强度的第i点；{V1,V2,V3,V4}表示棋盘的四个顶点；Fg(ri)确定点是否落在阈值[λL,λH]中，
Figure BDA0002365136510000085
表示顶点{V}的棋盘是否包含
Figure BDA0002365136510000086
点，
Figure BDA0002365136510000087
表示
Figure BDA0002365136510000088
点距离棋盘边缘X,Y方向的最小距离之和；
Figure BDA0002365136510000089
Figure BDA00023651365100000810
Figure BDA00023651365100000811
Figure BDA00023651365100000812
步骤S43：其余(m-1)块棋盘点云角点仍然根据步骤S41，步骤S42得到。
进一步地，所述步骤S5具体包括以下步骤：
步骤S51：从棋盘左下侧开始共同计数顺序，检测到的图像棋盘角点Ic与点云棋盘角点pL的角点对应；
步骤S52：利用公式(26)计算图像棋盘角点与点云棋盘角点的迭代次数，，然后选择4个不共面的控制点；通过公式(27)，将成像平面坐标系下的棋盘角点Ic转化到相机坐标系，得到了相机坐标系下的棋盘角点pc；
Figure BDA0002365136510000091
Figure BDA0002365136510000092
其中，(fx,fy)为相机焦距，(u0,v0)为相机主点，s为畸变系数；
步骤S53：如公式(28)所示，计算出pc,pL的中心，通过
Figure BDA0002365136510000093
得到矩阵
Figure BDA0002365136510000094
然后利用奇异值分解法H＝U∑VT，公式(29)求解出R,t；
Figure BDA0002365136510000095
Figure BDA0002365136510000096
步骤S54：将求解出来的R,t带入下式误差函数，计算出角点误差；选择最小的角点误差所对应的外部参数，为最终的外部参数R*,t*。
Figure BDA0002365136510000097
与现有技术相比，本发明具有以下有益效果：
(1)本发明仅需将多个棋盘放置于LIDAR和全景相机的共同视场下，只需要一次拍摄，即只需要一张全景图像与对应的点云数据，与以往激光相机校准方法相比，数据采集更加简单快捷。
(2)本发明通过公式(26)-(30)构建出点云特征点与图像特征点的几何约束方程，比通过点云特征线与图像特征线，或者点云特征面与图像特征面建立的几何约束方程进行外部校准的误差更小，更加准确。
(3)本发明外部校准过程，全程计算机计算，无需手动选择对应的棋盘角点。
附图说明
图1为本发明实施例的流程图。
图2为本发明实施例的安装LIDAR与全景相机的移动机器人图。
图3为本发明实施例的一个实验场景图。
图4为本发明实施例的图像棋盘角点相似程度的两种可能图，其中图4(a)为左对角线为黑，右对角线为白图，图4(b)为左对角线为白，右对角线为黑图。
图5为本发明实施例的云数据处理的实际效果图，其中图5(a)为去除x,y方向远距离点图，图5(b)为去除地面点云图，图5(c)为分割点云平面图，图5(d)为棋盘点云平面提取图。
图6为本发明实施例的棋盘点云平面估计出棋盘角点的实际效果图。
图7为本发明实施例的棋盘点云的降维过程图。
图8为本发明实施例的2D-3D棋盘角点建立过程图。
图9为本发明实施例的外部校准结果的实际投影效果图。
具体实施方式
下面结合附图及实施例对本发明做进一步说明。
应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
如图1所示，本实施例提供一种单次拍摄的LIDAR与全景相机的外参数标定方法，提供一Robotnik移动机器人，包括以下步骤：
步骤S1：将LIDAR(Velodyne-HDL-64e)与全景相机(Ladybug5)固定在Robotnik移动机器人上；然后将m个棋盘放置于LIDAR与全景相机的共同视场下，一次拍摄收集单帧的全景图像与该帧全景图像对应的点云数据；(本实例采用的是五个棋盘)
步骤S2：利用生长的棋盘角点检测算法，检测出全景图像的棋盘角点Ic；
步骤S3：对点云数据进行预处理，分割去除点云地面，分割点云平面、提取棋盘点云；
步骤S4：基于点云的反射强度，估计出点云棋盘角点pL；
步骤S5：通过定义从棋盘左下侧开始的角点共同计数顺序，建立全景图像的棋盘角点Ic与点云的棋盘角点pL的几何约束方程，求解出外参数R*,t*。
上述步骤S1是将LIDAR与全景相机安装在Robotnik移动机器人上，安装完成效果如图2所示。图3所示的是一个LIDAR与全景相机外部校准的实验场景，通过装载LIDAR与全景相机的移动机器人获取单帧的全景图像与此帧图像对应的点云数据。
在本实施例中，所述步骤S1具体包括以下步骤：
步骤S11：通过螺栓连接将LIDAR(Velodyne-HDL-64e)与全景相机(Ladybug5)固定在Robotnik移动机器人上；
步骤S12：构建一个室外场景，在场景中放置m块标定棋盘，每块棋盘的大小为600mm×450mm，棋盘中每个正方形的大小为75mm×75mm，并且满足多个棋盘在LIDAR与全景相机共同视场下的要求，其中，m取值为3、4、...11、12，m为整数；
步骤S13：利用步骤S11固定在移动机器人上的LIDAR与全景相机，使用全景相机收集一帧步骤S12构建场景的全景图像，LIDAR收集这帧全景图像对应的点云数据。
在本实施例中，步骤S2中所述检测出全景图像的棋盘角点Ic的具体内容为：
一般可以将生长的图像棋盘角点检测算法分为两个步骤：粗定位棋盘角点位置与进一步确定棋盘角点位置。
步骤S21：粗定位棋盘格角点的位置：首先定义两种不同类型的角点原型，原型1是一种和坐标轴平行的角点，原型2是一种与坐标轴成45°的角点；每个原型分别由4个卷积核组成，令原型1由四个卷积核K1,K2,K3,K4组成，原型2由四个卷积核K5,K6,K7,K8组成，分别用于与全景图像进行卷积操作；
通过两个角点原型来定义全景图像中每个像素点与角点的相似程度；
Figure BDA0002365136510000131
Figure BDA0002365136510000132
Figure BDA0002365136510000133
Figure BDA0002365136510000134
Figure BDA0002365136510000135
Figure BDA0002365136510000136
Figure BDA0002365136510000137
其中，如图4所示，
Figure BDA0002365136510000138
表示与原型1相似程度的两种可能，
Figure BDA0002365136510000139
表示与原型2相似程度的两种可能，原型1与原型2的相似程度的两种可能相同，表示的是左对角线为黑，右对角线为白，或者左对角线为白，右对角线为黑，
Figure BDA00023651365100001310
表示卷积核K1,K2,K3,K4原型1在某个像素的卷积值，
Figure BDA00023651365100001311
表示卷积核K5,K6,K7,K8原型2在某个像素的卷积值，c表示图像棋盘角点的最大相似程度；通过计算角点相似程度，得到大致的角点范围；然后通过非极大值抑止算法来获得候选角点cp；
步骤S22：由于步骤S21得到的候选角点不是很精确，需要进一步确定角点的位置；令c是理想的角点位置，p是c局部邻域的一个像素点，Gp是p点的图像梯度向量，此时满足如下式子：
Figure BDA0002365136510000141
由于实际图像中不止一个局部领域的像素点，因此在候选角点cp的邻域N(cp)内满足下列公式条件的就是所需要的棋盘角点Ic；
Figure BDA0002365136510000142
在本实施例中，对点云数据进行预处理，包括去除X,Y方向远距离点云，分割去除点云地面，分割点云平面，提取棋盘点云，图5展示的是点云数据处理的具体效果。
所述步骤S3具体包括以下步骤：
步骤S31：估计点云棋盘角点之前，对点云数据进行预处理；通过PCL中的直通滤波器模块，将点云PcL＝{(x,y,z)}中X,Y方向超过8m远的点剔除；
Figure BDA0002365136510000143
其中，pi＝(x,y,z)是点云PcL中的一点；
步骤S32：根据步骤S31将点云远点剔除后，基于形态学算法将点云地面分离，用以减小在点云平面分割时，地面点的干扰；形态学算法是根据设计出的膨胀因子和腐蚀因子，经过系列组合处理分割出地面点云。
Figure BDA0002365136510000144
Figure BDA0002365136510000145
其中，将点云测量点p(x,y,z)，x,y处的高度z定义为膨胀因子dp，腐蚀因子ep与膨胀因子相对应；w是测量点p的邻近窗口大小；
wj＝2jb+1 (13)
对于w窗口大小根据上述公式，线性增加窗口大小；其中j＝1,2,3,...,360，j为整数，b是初始窗口大小；
定义一个因子s来判断切除深度阈值dh；令物体坡度恒定，则最大坡度差为dhmax(t),j，因子s与窗口大小wj存在关系：
Figure BDA0002365136510000151
Figure BDA0002365136510000152
其中，dhT,j表示第j窗口切除深度阈值，将点云中的点依次带入公式(11)-(15)，计算出切除深度阈值dhT,j，若腐蚀因子大于切除深度阈值，则将该点移除，否则该点保存，其中T表示阈值，j表示第j个窗口大小；
步骤S33：根据步骤S32将点云地面移除后，剩余的点云集合为PrL；利用区域增长算法对点云平面进行分割；首先将剩余点云中每个点的曲率值
Figure BDA0002365136510000153
从小到大排序，曲率值越小，表示所在区域越平坦；然后，将曲率值最小点Pmin加入种子点，搜索它的K个最近邻点
Figure BDA0002365136510000154
计算每个点的法向量
Figure BDA0002365136510000155
并与最小法向量Nmin相比，若不大于平滑阈值Tth，令
Figure BDA0002365136510000161
π是圆周率,则该点加入点云平面RL:
Figure BDA0002365136510000162
若小于曲率阈值cth，令cth＝1.0，则将该点加入种子点，生长每个区域直到它收敛，并从PrL移除，如此循环，直到PrL中没有剩余的点；
最后使用棋盘平面度、形状大小条件，提取出棋盘点云
Figure BDA0002365136510000163
m表示棋盘的数量；
其中，提取棋盘点云
Figure BDA0002365136510000164
的具体内容：
由提取的点云平面RL组成矩阵Mn×3，沿着三个基矢量Mb＝(ux,uy,uz)T的方向分解，每个基矢量上的分量比为λ1,λ2,λ3；当最小比率λ3小于0.05并且单个棋盘的宽度dW和高度和dH满足公式(17)的点云片段被认为是棋盘点云；
Figure BDA0002365136510000165
其中W,H是单个棋盘的宽度和高度。
在本实施例中单个棋盘的宽度和高度为(600mm×450mm)。
在本实施例中，基于点云的反射强度，估计出棋盘点云角点pL，
图6展示的是估计出的棋盘角点。所述步骤S4具体包括以下步骤：
步骤S41：图7所示的是棋盘点云的降维过程，其中黑色表示低强度点云，白色表示高强度点云。将第k块棋盘点云
Figure BDA0002365136510000166
降低到二维平面并与第k块图像棋盘对应：使用主成分分析法，通过旋转矩阵
Figure BDA0002365136510000167
与平移矩阵
Figure BDA0002365136510000168
如公式(18)所示，将第k块棋盘点云
Figure BDA0002365136510000169
转化为与LIDAR坐标系一致的棋盘平面坐标系，其它棋盘点云跟随旋转平移变化；此过程中，求解出Cov矩阵的三个特征值(γ1,γ2,γ3)，得到对应的三个特征向量(v1,v2,v3)，旋转矩阵
Figure BDA0002365136510000171
定义为(v1,v2,v3)；其中，k表示第k块棋盘，取值范围是1-m；
Figure BDA0002365136510000172
Figure BDA0002365136510000173
Figure BDA0002365136510000174
步骤S42：第k块棋盘点云转化为与LIDAR坐标系一致的棋盘平面坐标系后，利用黑白棋盘格图案反射强度的对应关系，设置阈值[λL,λH]，设为[2.5,59]，小于λL表示从黑色图案反射的低强度，而大于λH表示从白色图案反射的高强度；
反射强度值在[λL,λH]为点云棋盘角点所在区域；通过黑白棋盘格反射强度的对应关系来制定成本函数，如公式(21)所示，从而估计出点云棋盘角点pL；
Figure BDA0002365136510000175
Figure BDA0002365136510000176
是
Figure BDA0002365136510000177
降维到XOY平面的第k块棋盘点云，其中，i表示第i点，coi表示
Figure BDA0002365136510000178
落入的图案颜色，令黑色为0，白色为1；ri是点云反射强度的第i点；{V1,V2,V3,V4}表示棋盘的四个顶点；Fg(ri)确定点是否落在阈值[λL,λH]中，
Figure BDA0002365136510000179
表示顶点{V}的棋盘是否包含
Figure BDA00023651365100001710
点，
Figure BDA00023651365100001711
表示
Figure BDA0002365136510000181
点距离棋盘边缘X,Y方向的最小距离之和；
Figure BDA0002365136510000182
Figure BDA0002365136510000183
Figure BDA0002365136510000184
Figure BDA0002365136510000185
步骤S43：其余(m-1)块棋盘点云角点仍然根据步骤S41，步骤S42得到。
如图8所示，在本实施例中，建立全景图像的棋盘角点Ic与点云的棋盘角点pL的几何约束方程，求解出外参数R*,t*。所述步骤S5具体包括以下步骤：
Figure BDA0002365136510000186
Figure BDA0002365136510000191
步骤S51：从棋盘左下侧开始共同计数顺序，检测到的图像棋盘角点Ic与点云棋盘角点pL的角点对应；
步骤S52：利用公式(26)计算图像棋盘角点与点云棋盘角点的迭代次数，，然后选择4个不共面的控制点；通过公式(27)，将成像平面坐标系下的棋盘角点Ic转化到相机坐标系，得到了相机坐标系下的棋盘角点pc；
Figure BDA0002365136510000192
Figure BDA0002365136510000193
其中，(fx,fy)为相机焦距，(u0,v0)为相机主点，s为畸变系数；
步骤S53：如公式(28)所示，计算出pc,pL的中心，通过
Figure BDA0002365136510000194
得到矩阵
Figure BDA0002365136510000195
然后利用奇异值分解法H＝U∑VT，公式(29)求解出R,t；
Figure BDA0002365136510000196
Figure BDA0002365136510000197
步骤S54：将求解出来的R,t带入下式误差函数，计算出角点误差；选择最小的角点误差所对应的外部参数，为最终的外部参数R*,t*。
Figure BDA0002365136510000201
在本实施例中，R是一个3*3矩阵，Rx是绕x轴的矩阵，Ry是绕y轴的矩阵，Rz是绕z轴的矩阵，旋转角度θ＝(θx,θy,θz)分别对应x,y,z轴的旋转角度，T＝(tx,ty,tz)分别为x,y,z轴的平移矢量。
R(θ)＝Rz(θz)Ry(θy)Rx(θx)
Figure BDA0002365136510000202
Figure BDA0002365136510000203
Figure BDA0002365136510000204
所以求解出来的R*,t*，是包括六个外参数Rx,Ry,Rz；tx,ty,tz。
较佳的，如图9所示，本实施例是通过构建点云特征点与图像特征点的几何约束关系来进行外部校准，更加准确。本实施例基于生长的图像棋盘角点检测算法与基于点云反射强度估计出棋盘角点，将外参数标定问题转换成2D-3D棋盘角点匹配的几何约束问题，只需要一次拍摄，就能实现LIDAR和全景相机的外参数标定。并且仅需将多个棋盘放置于LIDAR和全景相机的共同视场下，只进行一次拍摄，就能够基于生长的图像棋盘角点检测方法获取图像棋盘角点，根据点云反射强度估计出点云棋盘角点，然后建立2D-3D图像棋盘角点与点云棋盘角点的几何约束方程，自动计算出外部校准参数，实现单次LIDAR与全景相机传感器的外部校准。
以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。

Claims (6)

1.一种单次拍摄的LIDAR与全景相机的外参数标定方法，提供一Robotnik移动机器人，其特征在于：包括以下步骤：
步骤S1：将LIDAR与全景相机固定在所述Robotnik移动机器人上；然后将m个棋盘放置于LIDAR与全景相机的共同视场下，一次拍摄收集单帧的全景图像与该帧全景图像对应的点云数据；
步骤S2：利用生长的棋盘角点检测算法，检测出全景图像的棋盘角点Ic；
步骤S3：对点云数据进行预处理，分割去除点云地面，分割点云平面、提取棋盘点云；
步骤S4：基于点云的反射强度，估计出点云棋盘角点pL；
步骤S5：通过定义从棋盘左下侧开始的角点共同计数顺序，建立全景图像的棋盘角点Ic与点云的棋盘角点pL的几何约束方程，求解出外参数R*,t*。
2.根据权利要求1所述的一种单次拍摄的LIDAR与全景相机的外参数标定方法，其特征在于：所述步骤S1具体包括以下步骤：
步骤S11：通过螺栓连接将LIDAR与全景相机固定在Robotnik移动机器人上；
步骤S12：构建一个室外场景，在场景中放置m块标定棋盘，每块棋盘的大小为600mm×450mm，棋盘中每个正方形的大小为75mm×75mm，并且满足多个棋盘在LIDAR与全景相机共同视场下的要求，其中，m取值为3、4、...11、12，m为整数；
步骤S13：利用步骤S11固定在移动机器人上的LIDAR与全景相机，使用全景相机收集一帧步骤S12构建场景的全景图像，LIDAR收集这帧全景图像对应的点云数据。
3.根据权利要求1所述的一种单次拍摄的LIDAR与全景相机的外参数标定方法，其特征在于：步骤S2中所述检测出全景图像的棋盘角点Ic的具体内容为：
步骤S21：粗定位棋盘格角点的位置：首先定义两种不同类型的角点原型，原型1是一种和坐标轴平行的角点，原型2是一种与坐标轴成45°的角点；每个原型分别由4个卷积核组成，令原型1由四个卷积核K1,K2,K3,K4组成，原型2由四个卷积核K5,K6,K7,K8组成，分别用于与全景图像进行卷积操作；
通过两个角点原型来定义全景图像中每个像素点与角点的相似程度；
Figure FDA0002365136500000021
Figure FDA0002365136500000022
Figure FDA0002365136500000023
Figure FDA0002365136500000024
Figure FDA0002365136500000025
Figure FDA0002365136500000026
Figure FDA0002365136500000027
其中
Figure FDA0002365136500000028
表示与原型1相似程度的两种可能，
Figure FDA0002365136500000029
表示与原型2相似程度的两种可能，原型1与原型2的相似程度的两种可能相同，表示的是左对角线为黑，右对角线为白，或者左对角线为白，右对角线为黑，
Figure FDA00023651365000000210
表示卷积核K1,K2,K3,K4原型1在某个像素的卷积值，
Figure FDA0002365136500000031
表示卷积核K5,K6,K7,K8原型2在某个像素的卷积值，c表示图像棋盘角点的最大相似程度；通过计算角点相似程度，得到大致的角点范围；然后通过非极大值抑止算法来获得候选角点cp；
步骤S22：令c是理想的角点位置，p是c局部邻域的一个像素点，Gp是p点的图像梯度向量，此时满足如下式子：
Figure FDA0002365136500000032
由于实际图像中不止一个局部领域的像素点，因此在候选角点cp的邻域N(cp)内满足下列公式条件的就是所需要的棋盘角点Ic；
Figure FDA0002365136500000033
4.根据权利要求1所述的一种单次拍摄的LIDAR与全景相机的外参数标定方法，其特征在于：所述步骤S3具体包括以下步骤：
步骤S31：估计点云棋盘角点之前，对点云数据进行预处理；通过PCL中的直通滤波器模块，将点云PcL＝{(x,y,z)}中X,Y方向超过8m远的点剔除；
Figure FDA0002365136500000034
其中，pi＝(x,y,z)是点云PcL中的一点；
步骤S32：根据步骤S31将点云远点剔除后，基于形态学算法将点云地面分离，用以减小在点云平面分割时，地面点的干扰；
Figure FDA0002365136500000035
Figure FDA0002365136500000041
其中，将点云测量点p(x,y,z)，x,y处的高度z定义为膨胀因子dp，腐蚀因子ep与膨胀因子相对应；w是测量点p的邻近窗口大小；
wj＝2jb+1 (13)
对于w窗口大小根据上述公式，线性增加窗口大小；其中j＝1,2,3,...,360，j为整数，b是初始窗口大小；
定义一个因子s来判断切除深度阈值dh；令物体坡度恒定，则最大坡度差为dhmax(t),j，因子s与窗口大小wj存在关系：
Figure FDA0002365136500000042
Figure FDA0002365136500000043
其中，dhT,j表示第j窗口切除深度阈值，将点云中的点依次带入公式(11)-(15)，计算出切除深度阈值dhT,j，若腐蚀因子大于切除深度阈值，则将该点移除，否则该点保存，其中T表示阈值，j表示第j个窗口大小；
步骤S33：根据步骤S32将点云地面移除后，剩余的点云集合为PrL；利用区域增长算法对点云平面进行分割；首先将剩余点云中每个点的曲率值
Figure FDA0002365136500000044
从小到大排序，曲率值越小，表示所在区域越平坦；然后，将曲率值最小点Pmin加入种子点，搜索它的K个最近邻点
Figure FDA0002365136500000045
计算每个点的法向量
Figure FDA0002365136500000051
并与最小法向量Nmin相比，若不大于平滑阈值Tth，则该点加入点云平面RL:
Figure FDA0002365136500000052
若小于曲率阈值cth，则将该点加入种子点，生长每个区域直到它收敛，并从PrL移除，如此循环，直到PrL中没有剩余的点；
最后使用棋盘平面度、形状大小条件，提取出棋盘点云
Figure FDA0002365136500000053
m表示棋盘的数量；
其中，提取棋盘点云
Figure FDA0002365136500000054
的具体内容：
由提取的点云平面RL组成矩阵Mn×3，沿着三个基矢量Mb＝(ux,uy,uz)T的方向分解，每个基矢量上的分量比为λ1,λ2,λ3；当最小比率λ3小于0.05并且单个棋盘的宽度dW和高度和dH满足公式(17)的点云片段被认为是棋盘点云；
Figure FDA0002365136500000055
其中W,H是单个棋盘的宽度和高度。
5.根据权利要求1所述的一种单次拍摄的LIDAR与全景相机的外参数标定方法，其特征在于：所述步骤S4具体包括以下步骤：
步骤S41：将第k块棋盘点云
Figure FDA0002365136500000056
降低到二维平面并与第k块图像棋盘对应：使用主成分分析法，通过旋转矩阵
Figure FDA0002365136500000057
与平移矩阵
Figure FDA0002365136500000058
如公式(18)所示，将第k块棋盘点云
Figure FDA0002365136500000059
转化为与LIDAR坐标系一致的棋盘平面坐标系，其它棋盘点云跟随旋转平移变化；此过程中，求解出Cov矩阵的三个特征值(γ1,γ2,γ3)，得到对应的三个特征向量(v1,v2,v3)，旋转矩阵
Figure FDA00023651365000000510
定义为(v1,v2,v3)；其中，k表示第k块棋盘，取值范围是1-m；
Figure FDA0002365136500000061
Figure FDA0002365136500000062
Figure FDA0002365136500000063
步骤S42：第k块棋盘点云转化为与LIDAR坐标系一致的棋盘平面坐标系后，利用黑白棋盘格图案反射强度的对应关系，设置阈值[λL,λH]，小于λL表示从黑色图案反射的低强度，而大于λH表示从白色图案反射的高强度；
反射强度值在[λL,λH]为点云棋盘角点所在区域；通过黑白棋盘格反射强度的对应关系来制定成本函数，如公式(21)所示，从而估计出点云棋盘角点pL；
Figure FDA0002365136500000064
Figure FDA00023651365000000611
是
Figure FDA0002365136500000065
降维到XOY平面的第k块棋盘点云，其中，i表示第i点，coi表示
Figure FDA0002365136500000066
落入的图案颜色，令黑色为0，白色为1；ri是点云反射强度的第i点；{V1,V2,V3,V4}表示棋盘的四个顶点；Fg(ri)确定点是否落在阈值[λL,λH]中，
Figure FDA0002365136500000067
表示顶点{V}的棋盘是否包含
Figure FDA0002365136500000068
点，
Figure FDA0002365136500000069
表示
Figure FDA00023651365000000610
点距离棋盘边缘X,Y方向的最小距离之和；
Figure FDA0002365136500000071
Figure FDA0002365136500000072
Figure FDA0002365136500000073
Figure FDA0002365136500000074
步骤S43：其余(m-1)块棋盘点云角点仍然根据步骤S41，步骤S42得到。
6.根据权利要求1所述的一种单次拍摄的LIDAR与全景相机的外参数标定方法，其特征在于：所述步骤S5具体包括以下步骤：
步骤S51：从棋盘左下侧开始共同计数顺序，检测到的图像棋盘角点Ic与点云棋盘角点pL的角点对应；
步骤S52：利用公式(26)计算图像棋盘角点与点云棋盘角点的迭代次数，，然后选择4个不共面的控制点；通过公式(27)，将成像平面坐标系下的棋盘角点Ic转化到相机坐标系，得到了相机坐标系下的棋盘角点pc；
Figure FDA0002365136500000075
Figure FDA0002365136500000076
其中，(fx,fy)为相机焦距，(u0,v0)为相机主点，s为畸变系数；
步骤S53：如公式(28)所示，计算出pc,pL的中心，通过
Figure FDA0002365136500000081
得到矩阵
Figure FDA0002365136500000082
然后利用奇异值分解法H＝U∑VT，公式(29)求解出R,t；
Figure FDA0002365136500000083
Figure FDA0002365136500000084
步骤S54：将求解出来的R,t带入下式误差函数，计算出角点误差；选择最小的角点误差所对应的外部参数，为最终的外部参数R*,t*。
Figure FDA0002365136500000085
CN202010034949.9A 2020-01-13 2020-01-13 一种单次拍摄的lidar与全景相机的外参数标定方法 Active CN111260735B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010034949.9A CN111260735B (zh) 2020-01-13 2020-01-13 一种单次拍摄的lidar与全景相机的外参数标定方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010034949.9A CN111260735B (zh) 2020-01-13 2020-01-13 一种单次拍摄的lidar与全景相机的外参数标定方法

Publications (2)

Publication Number Publication Date
CN111260735A true CN111260735A (zh) 2020-06-09
CN111260735B CN111260735B (zh) 2022-07-01

Family

ID=70948727

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010034949.9A Active CN111260735B (zh) 2020-01-13 2020-01-13 一种单次拍摄的lidar与全景相机的外参数标定方法

Country Status (1)

Country Link
CN (1) CN111260735B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113281723A (zh) * 2021-05-07 2021-08-20 北京航空航天大学 一种基于AR tag的3D激光雷达与相机间结构参数的标定方法
EP4086846A1 (en) * 2021-05-03 2022-11-09 The Boeing Company Automatic detection of a calibration standard in unstructured lidar point clouds

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103049912A (zh) * 2012-12-21 2013-04-17 浙江大学 一种基于任意三面体的雷达-相机系统外部参数标定方法
CN105096317A (zh) * 2015-07-03 2015-11-25 吴晓军 一种复杂背景中的高性能相机全自动标定方法
US20180096493A1 (en) * 2017-12-04 2018-04-05 GM Global Technology Operations LLC Detection and recalibration for a camera system using lidar data
CN110161485A (zh) * 2019-06-13 2019-08-23 同济大学 一种激光雷达与视觉相机的外参标定装置及标定方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103049912A (zh) * 2012-12-21 2013-04-17 浙江大学 一种基于任意三面体的雷达-相机系统外部参数标定方法
CN105096317A (zh) * 2015-07-03 2015-11-25 吴晓军 一种复杂背景中的高性能相机全自动标定方法
US20180096493A1 (en) * 2017-12-04 2018-04-05 GM Global Technology Operations LLC Detection and recalibration for a camera system using lidar data
CN110161485A (zh) * 2019-06-13 2019-08-23 同济大学 一种激光雷达与视觉相机的外参标定装置及标定方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
SURABHI VERMA ET AL.: "Automatic extrinsic calibration between a camera and a 3D Lidar using 3D point and plane correspondences", 《2019 IEEE INTELLIGENT TRANSPORTATION SYSTEMS CONFERENCE (ITSC)》 *
ZOU,CHENG ET AL.: "Learning motion field of LiDAR point cloud with convolutional networks", 《PATTERN RECOGNITION LETTERS》 *
黎云飞 等: "基于fmincon法的单线激光雷达与单目相机外参数标定法", 《工业控制计算机》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP4086846A1 (en) * 2021-05-03 2022-11-09 The Boeing Company Automatic detection of a calibration standard in unstructured lidar point clouds
CN113281723A (zh) * 2021-05-07 2021-08-20 北京航空航天大学 一种基于AR tag的3D激光雷达与相机间结构参数的标定方法

Also Published As

Publication number Publication date
CN111260735B (zh) 2022-07-01

Similar Documents

Publication Publication Date Title
CN104574347B (zh) 基于多源遥感数据的在轨卫星图像几何定位精度评价方法
CN112949478B (zh) 基于云台相机的目标检测方法
CN109785379A (zh) 一种对称物体尺寸及重量的测量方法及测量系统
CN103218787B (zh) 多源异构遥感影像控制点自动采集方法
US7747106B2 (en) Method and system for filtering, registering, and matching 2.5D normal maps
CN110473221B (zh) 一种目标物体自动扫描系统及方法
CN109118544B (zh) 基于透视变换的合成孔径成像方法
CN107977996B (zh) 基于靶标标定定位模型的空间目标定位方法
CN113012234B (zh) 基于平面变换的高精度相机标定方法
CN111739031B (zh) 一种基于深度信息的作物冠层分割方法
CN111260735B (zh) 一种单次拍摄的lidar与全景相机的外参数标定方法
CN110930382A (zh) 基于标定板特征点提取的点云拼接精度评估方法及系统
CN113642463B (zh) 一种视频监控和遥感图像的天地多视图对齐方法
CN109214254A (zh) 一种确定机器人位移的方法及装置
CN115201883A (zh) 一种运动目标视频定位测速系统及方法
CN115456870A (zh) 基于外参估计的多图像拼接方法
CN113658279B (zh) 相机内参和外参估算方法、装置、计算机设备及存储介质
CN113589263B (zh) 一种多个同源传感器联合标定方法及系统
CN114549780A (zh) 一种基于点云数据的大型复杂构件智能化检测方法
CN114299153A (zh) 一种超大电力设备的相机阵列同步标定方法及系统
CN111735447B (zh) 一种仿星敏式室内相对位姿测量系统及其工作方法
EP2791865A1 (en) System and method for estimating target size
CN116468621A (zh) 一种一键式数字化航空影像数据处理方法
CN117036235A (zh) 一种继电保护柜端子排线序检测方法
Jarron et al. Automatic detection and labelling of photogrammetric control points in a calibration test field

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant