【2D与3D SLAM中的扫描匹配算法全面解析】
引言
扫描匹配(Scan Matching)是同步定位与地图构建(SLAM)系统中的核心组件,它通过对齐连续的传感器观测数据来估计机器人的运动。本文将深入探讨2D和3D SLAM中的各种扫描匹配算法,包括数学原理、实现细节以及实际应用中的性能对比,特别关注激光雷达SLAM中的典型方法。
一、扫描匹配数学基础与核心原理
1.1 刚体变换的数学表示
扫描匹配的核心是求解刚体变换,在2D和3D空间中有不同的数学表示:
2D刚体变换:
由旋转矩阵R和平移向量t组成,可表示为:
T = [ cos θ − sin θ t x sin θ cos θ t y 0 0 1 ] T = \begin{bmatrix} \cos\theta & -\sin\theta & t_x \\ \sin\theta & \cos\theta & t_y \\ 0 & 0 & 1 \end{bmatrix} T= cosθsinθ0−sinθcosθ0txty1
其中 θ θ θ 为旋转角度, ( t x , t y ) (t_x, t_y) (tx,ty) 为平移量。
3D刚体变换:
使用齐次坐标表示:
T = [ R 3 × 3 t 3 × 1 0 1 × 3 1 ] T = \begin{bmatrix} R_{3×3} & t_{3×1} \\ 0_{1×3} & 1 \end{bmatrix} T=[R3×301×3t3×11]
其中 R ∈ S O ( 3 ) R∈SO(3) R∈SO(3) 为旋转矩阵, t ∈ R 3 t∈ℝ³ t∈R3 为平移向量。
1.2 优化问题的构建
扫描匹配本质上是非线性最小二乘问题:
min T ∑ i ρ ( ∥ T ∘ p i − q i ∥ 2 ) \min_T \sum_i \rho( \| T \circ p_i - q_i \|^2 ) Tmini∑ρ(∥T∘pi−qi∥2)
其中 ρ ρ ρ 为鲁棒核函数(如Huber、Tukey),用于抑制异常值影响。
目标函数线性化:
通过一阶泰勒展开近似:
e ( x + Δ x ) ≈ e ( x ) + J ( x ) Δ x e(x+\Delta x) ≈ e(x) + J(x)\Delta x e(x+Δx)≈e(x)+J(x)Δx
其中 J J J 为雅可比矩阵, Δ x Δx Δx 为参数增量。
二、2D扫描匹配算法原理详解
2.1 点到点ICP的数学推导
误差函数:
对于单个点对 ( p i , q i ) (p_i, q_i) (pi,qi),误差为:
e i = R p i + t − q i e_i = R p_i + t - q_i ei=Rpi+t−qi
雅可比矩阵计算:
对2D变换参数 x = [ t x , t y , θ ] T x=[t_x, t_y, θ]^T x=[tx,ty,θ]T 求导:
J i = ∂ e i ∂ x = [ 1 0 − p i x sin θ − p i y cos θ 0 1 p i x cos θ − p i y sin θ ] J_i = \frac{\partial e_i}{\partial x} = \begin{bmatrix} 1 & 0 & -p_i^x \sin\theta - p_i^y \cos\theta \\ 0 & 1 & p_i^x \cos\theta - p_i^y \sin\theta \end{bmatrix} Ji=∂x∂ei=[1001−pixsinθ−piycosθpixcosθ−piysinθ]
高斯牛顿法更新:
通过求解正规方程获得参数更新:
( J T J ) Δ x = − J T e (J^T J) \Delta x = -J^T e (JTJ)Δx=−JTe
2.2 点到线ICP的几何原理
线段表示:
目标扫描中的线段由两点 q 1 q₁ q1, q 2 q₂ q2 确定,其直线方程为:
( q 2 − q 1 ) × ( p − q 1 ) = 0 (q₂-q₁) \times (p - q₁) = 0 (q2−q1)×(p−q1)=0
距离计算:
点 p p p 到线段的距离:
d = ∣ ( q 2 − q 1 ) × ( p − q 1 ) ∣ ∣ q 2 − q 1 ∣ d = \frac{| (q₂-q₁) \times (p-q₁) |}{| q₂-q₁ |} d=∣q2−q1∣∣(q2−q1)×(p−q1)∣
误差函数:
e i = ( q 2 − q 1 ) × ( R p i + t − q 1 ) ∣ q 2 − q 1 ∣ e_i = \frac{(q₂-q₁) \times (R p_i + t - q₁)}{| q₂-q₁ |} ei=∣q2−q1∣(q2−q1)×(Rpi+t−q1)
雅可比矩阵:
对旋转角度 θ θ θ 的导数为:
∂ e i ∂ θ = ( q 2 − q 1 ) × ( ∂ R ∂ θ p i ) ∣ q 2 − q 1 ∣ \frac{\partial e_i}{\partial \theta} = \frac{(q₂-q₁) \times ( \frac{\partial R}{\partial \theta} p_i )}{| q₂-q₁ |} ∂θ∂ei=∣q2−q1∣(q2−q1)×(∂θ∂Rpi)
2.3 似然场法的概率解释
似然场构建:
- 将目标扫描转换为二值栅格地图
- 对每个栅格计算到最近障碍物的距离 d d d
- 定义概率场:
P ( p ) = exp ( − d 2 2 σ 2 ) P(p) = \exp(-\frac{d^2}{2\sigma^2}) P(p)=exp(−2σ2d2)
优化目标:
最大化源扫描点在似然场中的总概率:
max T ∏ i P ( T ∘ p i ) \max_T \prod_i P(T \circ p_i) Tmaxi∏P(T∘pi)
取负对数后转化为最小化问题:
min T ∑ i d ( T ∘ p i ) 2 2 σ 2 \min_T \sum_i \frac{d(T \circ p_i)^2}{2\sigma^2} Tmini∑2σ2d(T∘pi)2
三、3D扫描匹配算法
3.1 3D点到点ICP
原理延续性:
3D点到点ICP是2D版本的直接扩展,核心思想完全一致:
- 建立点对点对应关系(最近邻搜索)
- 最小化对应点间的欧氏距离
- 通过SVD或非线性优化求解刚体变换
数学形式扩展:
- 旋转矩阵 R ∈ S O ( 3 ) R ∈ SO(3) R∈SO(3)(从2D的 θ θ θ 扩展为3D旋转)
- 平移向量 t ∈ R 3 t ∈ ℝ³ t∈R3( z z z 方向扩展)
- 误差函数: e i = R p i + t − q i e_i = R p_i + t - q_i ei=Rpi+t−qi
关键差异:
-
最近邻搜索复杂度:
- 2D:可用KD-tree或栅格搜索(O(n log n))
- 3D:必须使用KD-tree/Octree(O(n log n)但常数项更大)
-
变换求解方法:
-
2D:可直接解析求解或高斯牛顿法
-
3D:通常采用SVD分解(更稳定):
W = ∑ ( p i − p ˉ ) ( q i − q ˉ ) T W = \sum (p_i - \bar{p})(q_i - \bar{q})^T W=∑(pi−pˉ)(qi−qˉ)T
U Σ V T = SVD ( W ) UΣV^T = \text{SVD}(W) UΣVT=SVD(W)
R = V U T , t = q ˉ − R p ˉ R = V U^T, \quad t = \bar{q} - R \bar{p} R=VUT,t=qˉ−Rpˉ
-
PCL实现示例:
pcl::IterativeClosestPoint<pcl::PointXYZ, pcl::PointXYZ> icp;
icp.setInputSource(source_cloud);
icp.setMaxCorrespondenceDistance(0.05); // 比2D更小的距离阈值
icp.align(output_cloud);
3.2 3D点到线ICP
原理延续性:
与2D点到线ICP思想相同,但:
- 3D中的"线"通常来自边缘特征提取
- 距离计算使用3D空间中的点线距离公式
数学形式:
给定目标扫描中的直线(方向向量 v v v,通过点 q q q ):
d = ∥ ( R p i + t − q ) × v ∥ / ∥ v ∥ d = \| (R p_i + t - q) \times v \| / \| v \| d=∥(Rpi+t−q)×v∥/∥v∥
实现关键点:
-
线特征提取:
pcl::PointCloud<pcl::PointXYZ>::Ptr edge_points(new pcl::PointCloud<pcl::PointXYZ>); pcl::extractEdgePoints(source_cloud, edge_points); // 基于曲率或PCA方法
-
误差雅可比:
J = ∂ d ∂ [ t , ω ] = v × ∥ v ∥ ⋅ [ I 3 × 3 , − [ R p i ] × ] J = \frac{\partial d}{\partial [t, \omega]} = \frac{v \times}{\|v\|} \cdot [I_{3×3}, -[R p_i]_\times] J=∂[t,ω]∂d=∥v∥v×⋅[I3×3,−[Rpi]×]
(其中 ω ω ω 为角速度的李代数表示)
适用场景:
- 室内结构化环境(门框、桌腿等)
- 激光SLAM中的边缘跟踪(如LOAM)
3.3 3D点到面ICP
平面表示:
目标点云中的平面由点 q q q 和法向量 n n n 确定,平面方程为:
n T ( p − q ) = 0 n^T (p - q) = 0 nT(p−q)=0
误差函数:
e i = n j T ( R p i + t − q j ) e_i = n_j^T (R p_i + t - q_j) ei=njT(Rpi+t−qj)
雅可比矩阵:
对旋转部分采用李代数 s o ( 3 ) so(3) so(3) 参数化:
∂ e i ∂ ξ = − n j T [ R p i ] × \frac{\partial e_i}{\partial \xi} = -n_j^T [ R p_i ]_\times ∂ξ∂ei=−njT[Rpi]×
其中 [ ⋅ ] × [·]× [⋅]× 表示向量的叉积矩阵, ξ ∈ s o ( 3 ) ξ∈so(3) ξ∈so(3) 为旋转的李代数表示。
3.4 NDT算法的统计原理
体素划分:
将目标点云划分为体素网格,对每个体素内的点 q k {q_k} qk 计算:
- 均值: μ = ( 1 / n ) ∑ q k μ = (1/n)∑q_k μ=(1/n)∑qk
- 协方差: Σ = ( 1 / n ) ∑ ( q k − μ ) ( q k − μ ) T Σ = (1/n)∑(q_k-μ)(q_k-μ)^T Σ=(1/n)∑(qk−μ)(qk−μ)T
概率密度函数:
p ( p ) = exp ( − ( p − μ ) T Σ − 1 ( p − μ ) 2 ) p(p) = \exp \left( -\frac{(p-μ)^T Σ^{-1} (p-μ)}{2} \right) p(p)=exp(−2(p−μ)TΣ−1(p−μ))
优化目标:
最大化源点云在NDT场中的总概率:
min T ∑ i ( T p i − μ i ) T Σ i − 1 ( T p i − μ i ) \min_T \sum_i (T p_i - μ_i)^T Σ_i^{-1} (T p_i - μ_i) Tmini∑(Tpi−μi)TΣi−1(Tpi−μi)
Hessian矩阵计算:
NDT的高效实现依赖于精确计算Hessian矩阵:
H = J T Σ − 1 J H = J^T Σ^{-1} J H=JTΣ−1J
四、PCL算法实现原理对比
4.1 点到面ICP的实现优化
PCL中的pcl::IterativeClosestPointWithNormals
通过以下优化提升性能:
-
法线估计加速:
- 使用KD-tree近邻搜索
- 基于PCA计算法向量
pcl::NormalEstimation<pcl::PointXYZ, pcl::Normal> ne; ne.setInputCloud(cloud); ne.setSearchMethod(tree); ne.setKSearch(30); // 使用30个最近邻计算法线
-
对应点筛选策略:
- 法线夹角阈值过滤
- 距离阈值过滤
4.2 GICP的协方差传播
广义ICP(GICP)通过考虑局部几何结构提升精度:
源点协方差:
Σ i s r c = R Σ i s r c , l o c a l R T Σ_i^{src} = R Σ_i^{src,local} R^T Σisrc=RΣisrc,localRT
目标点协方差:
Σ i t g t = Σ i t g t , l o c a l Σ_i^{tgt} = Σ_i^{tgt,local} Σitgt=Σitgt,local
马氏距离:
d i = ( T p i − q i ) T ( Σ i t g t + Σ i s r c ) − 1 ( T p i − q i ) d_i = (T p_i - q_i)^T (Σ_i^{tgt} + Σ_i^{src})^{-1} (T p_i - q_i) di=(Tpi−qi)T(Σitgt+Σisrc)−1(Tpi−qi)
4.3 NDT的多分辨率策略
PCL中的NDT实现采用多分辨率加速:
- 初始使用大体素进行粗配准
- 逐步缩小体素尺寸提高精度
ndt.setResolution(2.0); // 初始分辨率
ndt.setStepSize(0.1); // 步长控制
ndt.setOulierRatio(0.55);// 异常值比例
4.4 对比总结
方法类型 | PCL类名 | 优点 | 缺点 | 适用场景 |
---|---|---|---|---|
点到点ICP | pcl::IterativeClosestPoint | 实现简单,通用性强 | 收敛慢,对初始位置敏感 | 通用3D配准 |
点到面ICP | pcl::IterativeClosestPointWithNormals | 收敛快,精度高 | 需要法线计算 | 结构化环境 |
GICP | pcl::GeneralizedIterativeClosestPoint | 结合点和面信息,更鲁棒 | 计算复杂度高 | 复杂环境 |
NDT | pcl::NormalDistributionsTransform | 对初始位置不敏感,效率高 | 依赖网格分辨率设置 | 大场景,自动驾驶 |
Feature-based | pcl::SampleConsensusInitialAlignment | 利用特征,全局配准 | 依赖特征提取质量 | 无初始猜测的情况 |
五、算法选择的数学依据
5.1 收敛性分析
方法 | 收敛半径 | 收敛速度 | 数学特性 |
---|---|---|---|
点到点ICP | 小(~1m) | 线性 | 非凸,局部极值多 |
点到面ICP | 中等(~3m) | 超线性 | 利用几何约束,更凸 |
NDT | 大(~10m) | 二次 | 平滑概率场,全局性更好 |
5.2 计算复杂度比较
算法时间复杂度分析:
- ICP类算法: O ( n m k ) O(nmk) O(nmk),其中 n n n 为源点数, m m m 为目标点数, k k k 为迭代次数
- NDT: O ( n + v ) O(n + v) O(n+v), v v v 为体素数量,预处理阶段 O ( m ) O(m) O(m)
- 特征匹配: O ( n + m + f ) O(n + m + f) O(n+m+f), f f f 为特征匹配耗时
这部分的内容还是很重要的,我之前学得很浅,然后到后面实际应用的时候总是搞不明白扫描匹配这一步,所以还是回过头来重新学了。