论文链接：https://arxiv.org/pdf/2103.01627.pdf
代码链接：https://github.com/hku-mars/livox_camera_calib
主要内容
本文提出了一种不依赖于棋盘格等辅助标定物体，实现像素级相机和激光雷达自动标定的方法。
方法
直接从点云中提取3D边特征，一避免遮挡问题，并且使用了精确度更高的深度连续边。
文中首先指出：以下四种变换方式区分性能不好，提供的约束效果不好：
进一步指出将点云投影到图像平面再提取特征的问题主要在于相机和激光雷达事业遮挡形成的零值映射和多值映射：
因此本文直接从点云中提取线特征。
目前的大部分工作是针对深度不连续的线特征展开的，因为它们根据深度信息的不同可以被轻松的区别出来，但是这种方法的准确性不如深度连续的线特征，由于激光雷达发射的是光束(beam)而不是理想的点，会存在foreground inflation和bleeding points现象：
边缘提取
因此本文采用深度连续的边特征，首先将点云分成特定尺寸的小体素。针对每一个体素，我们重复地使用RANSAC来适配和提取包含在体素中的点进行平面拟合，我们保留夹角在30-150度之间的，连接的平面对，并且求出其交线，如图4中的深度连续边，作为其边特征。我们的方法可以提取出或平行或垂直的边缘。更进一步，选择适当的体素大小，我们的方法甚至可以提取曲线边缘。
对于图像中的边缘提取，我们采用了Canny算法，边缘的像素被保存到了k-D tree中，进行匹配。
边缘匹配
激光雷达提取的边需要和图像中的边做相应的匹配，激光雷达的点到相机平面的变换可以用如下的式子进行表达：
C
P
i
=
L
C
P
i
(
L
P
i
)
∈
R
3
^{C}P_i=^{C}_{L}P_i(^{L}P_i)\in \mathbb{R}^3
CPi=LCPi(LPi)∈R3
进一步将
C
P
i
^{C}P_i
CPi投影到图像坐标系下：
C
p
i
=
π
(
C
P
i
)
^{C}p_i=\pi(^{C}P_i)
Cpi=π(CPi)
π
(
P
)
\pi(P)
π(P)是针孔相机的投影模型，经过畸变处理后的投影点真实位置
p
i
=
(
u
i
,
v
i
)
p_i=(u_i,v_i)
pi=(ui,vi):
p
i
=
f
(
C
p
i
)
p_i=f(^{C}p_i)
pi=f(Cpi)
这里
f
(
p
)
f(p)
f(p)是相机的畸变模型。
我们在图像边构建好的k-d树中寻找
p
i
p_i
pi的前k个最近邻，记为
Q
i
=
{
q
i
j
;
j
=
1
,
.
.
.
,
k
}
Q_i=\{q_i^j;j=1,...,k\}
Qi={qij;j=1,...,k},并且记：
q
i
=
1
k
∑
j
=
1
k
q
i
j
;
S
i
=
1
k
∑
j
=
1
k
(
q
i
j
−
q
i
)
(
q
i
j
−
q
i
)
T
q_i=\frac{1}{k}\sum_{j=1}^kq_i^j;\quad S_i=\frac{1}{k}\sum_{j=1}^k(q_i^j-q_i)(q_i^j-q_i)^T
qi=k1j=1∑kqij;Si=k1j=1∑k(qij−qi)(qij−qi)T
由
Q
i
Q_i
Qi组成的线，被
q
i
q_i
qi和
n
i
n_i
ni两个参数所确定，
n
i
n_i
ni是
S
i
S_i
Si最小特征值对应的特征向量。
同时，将雷达点云边的方向信息投影到图像平面上，验证其和
n
i
n_i
ni的正交性，可以有效减少两条接近的非平行直线的误匹配。匹配结果如下图所示：
外参标定
测量噪声
提取出的雷达点
L
P
i
^LP_i
LPi和和对应的边特征
(
n
i
,
q
i
)
(n_i,q_i)
(ni,qi)在图像平面中，会受到测量噪声的影响。记
I
w
i
∈
N
(
0
,
I
∑
i
)
^Iw_i\in \mathcal{N}(0,^I\sum_i)
Iwi∈N(0,I∑i)为图像边缘提取中关于
q
i
q_i
qi的噪声，它的协方差记为
I
∑
i
=
δ
2
I
2
×
2
^I\sum_i=\delta^2 I_{2×2}
I∑i=δ2I2×2，
δ
I
=
1.5
\delta_I=1.5
δI=1.5表示像素误差。
对于雷达点
L
P
i
^LP_i
LPi，记
L
w
i
^Lw_i
Lwi为它的测量噪声。事实上，激光雷达通过扫描电机的编码器来计算方向，通过光的飞行时间计算深度。记
w
i
∈
S
w_i\in \mathbb{S}
wi∈S为方向，
δ
w
i
∼
N
(
0
2
×
1
,
∑
w
i
)
\delta{w_i}\sim\mathcal{N}(0_{2×1},\sum_{w_i})
δwi∼N(02×1,∑wi)是
w
i
w_i
wi的正切平面上的测量噪声。使用在
S
\mathbb{S}
S中封闭的运算符
田
田
田，得到了
w
i
w_i
wi和真值
w
i
g
t
w_i^{gt}
wigt之间的关系：
w
i
g
t
=
w
i
田
S
2
δ
w
i
≜
e
⌊
N
(
w
i
δ
w
i
)
×
⌋
w
i
w_i^{gt}=w_i 田_{\mathbb{S}^2} \;\delta_{w_i}\triangleq e^{{\left \lfloor N(w_i\delta_{w_i})× \right \rfloor}_{w_i}}
wigt=wi田S2δwi≜e⌊N(wiδwi)×⌋wi
这里，
N
(
w
i
)
=
[
N
1
N
2
]
∈
R
3
×
2
N(w_i)=[N_1\quad N_2]\in \mathcal{R}^{3×2}
N(wi)=[N1N2]∈R3×2是
w
i
w_i
wi处切平面的正交基，
⌊
×
⌋
\left \lfloor × \right \rfloor
⌊×⌋表示向量叉积对应的偏对称矩阵，
田
田
田表示在
w
i
w_i
wi的切平面关于
δ
w
i
\delta_{w_i}
δwi轴旋转
w
i
w_i
wi的操作。
类似的，记
d
i
d_i
di为关于深度的测量，
δ
d
i
∼
N
(
0
,
∑
d
i
)
\delta_{d_i}\sim\mathcal{N}(0,\sum_{d_i})
δdi∼N(0,∑di)为其误差值，然后深度的真值可以表示为：
d
i
g
t
=
d
i
+
δ
d
i
d_i^{gt}=d_i+\delta_{d_i}
digt=di+δdi
结合上面两个式子，我们可以得到点位置的真值
L
P
i
g
t
^LP_i^{gt}
LPigt和测量值
L
P
i
^LP_i
LPi之间的关系：
所以：
这个误差模型可以用于产生一致性的外参标定。
标定的形式化和优化
从雷达点云中提取出来的边缘点
L
P
i
^LP_i
LPi可以用投影到图像中的对应边上点
q
i
∈
R
2
q_i\in \mathbb{R}^2
qi∈R2以及法线向量
n
i
∈
S
1
n_i\in \mathbb{S}^1
ni∈S1.进行表示，对
L
P
i
^LP_i
LPi的噪声进行补并且使用真实外参进行补偿，
L
P
i
^LP_i
LPi的投影应该刚好落在
(
n
i
,
q
i
)
(n_i,q_i)
(ni,qi)
上式说明一个雷达边缘点对外参提供一个约束，和上述的一个边特征对外参提供两个约束一致，因为一条边特征包含两个独立的点。同时，上式为外参
L
C
T
^C_LT
LCT，
L
P
i
,
n
i
,
q
i
^LP_i,n_i,q_i
LPi,ni,qi和未知的噪声
L
w
i
,
I
w
i
^Lw_i,^Iw_i
Lwi,Iwi提供了一个非线性不等式。这个非线性不等式可以通过迭代求解：记
L
C
T
‾
^C_L\overline{T}
LCT为当前的外参估计，使用运算符
田
田
田在
L
C
T
‾
^C_L\overline{T}
LCT的切平面参数化
L
C
T
^C_LT
LCT：
L
C
T
=
L
C
T
‾
田
S
E
(
3
)
δ
T
≜
E
x
p
(
δ
T
)
⋅
L
C
T
‾
^C_LT=^C_L\overline{T} 田 _{SE(3)}\delta T\triangleq Exp(\delta T) \cdot ^C_L\overline{T}
LCT=LCT田SE(3)δT≜Exp(δT)⋅LCT
这里：
将上一个等式代入上上个等式可得：
这里：
定义了一条边对应的约束，使用全部的
N
N
N条边可以得到下面的等式：
这里：
上式说明：
根据上式，我们提出了极大似然的外参估计（也是最小化差值v）如下：
优化的方向是：
对当前的
L
C
T
‾
^C_L\overline{T}
LCT进行更新：
L
C
T
‾
←
L
C
T
‾
田
S
E
(
3
)
δ
T
∗
^C_L\overline{T}\gets ^C_L\overline{T}田_{SE(3)}\; \delta T^*
LCT←LCT田SE(3)δT∗
不断进行迭代，直至
∣
∣
δ
T
∗
∣
∣
<
ϵ
||\delta T^*||<\epsilon
∣∣δT∗∣∣<ϵ，
L
C
T
‾
^C_L\overline{T}
LCT就是所求的标定外参。
标定不确定性
除了对外参进行估计，计算真实值与估计值之间误差的协方差也很有必要。为了达到这个目的，我们在全部的
N
N
N条边的等式约束两边同乘,
J
T
T
(
J
w
∑
J
w
T
)
J_T^T(J_w\sum J_w^T)
JTT(Jw∑JwT)然后解得
δ
T
\delta T
δT为：
进一步求得协方差矩阵：
边缘分布对标定结果的影响的分析
实验结果：
参考链接：泡泡机器人