Ground Estimation and Point Cloud Segmentation using SpatioTemporal Conditional Random Field （精读）

欢迎访问我的个人博客：zengzeyu.com

 

前言

---

译者按：时空条件随机场（Spatio-Temporal Conditional Random Field，STCRF），所谓“时空”在此论文中：“时”表示，模型与前一时刻 t-1 的地面估计相关；“空”表示与当前时刻 t 的当前空间分布有关。（具体参照正文中的时空模型）

 

摘要（Abstract）

---

无论是生成训练集还是障碍物识别跟踪系统，抑或是合理生成栅格地图，3D点云中地面点的提取分割分类操作是至关重要的一步，整个感知网络都可建立在其有效性上。在点云数据中，地平面估计或模式识别的不合理有可能导致系统性错误和运算量的激增。本文描述一种基于局部点云高度自适应的3D地面点云标记方法。此系统将地面建模为时空条件随机场（Spatio-Temporal Conditional Random Field， STCRF）。分割过程中的时空依赖性由基于 CRF 的概率框架来统一。采用了交叉互连（interconnected）期望最大化（Expection Maximization）算法变量的方法，由并行计算的每个节点一起生成地平面高度参数。该方法为搭载了稠密点云（Velodyne-64）和稀疏点云（IbeoLux）传感器的高速行驶车辆而设计，并且已经在嵌入式系统上进行算了试验（Nvidia Jetson TX1, TK1）。设计的试验场景包括实际各种道路场景（城市，农村，山区道路等），结果都达到了预期效果。

 

Ⅰ. 介绍和相关工作（INTRODUCTION AND RELATED WORK）

---

对于地面点云分割任务，本文总结了4个方法领域：

• 高度地图方法（Elevation map methods）： 将3D点云处理为2D栅格地图，在栅格内使用高度差等方法进行分类。此方法虽然效率很高，但是欠分割情况严重且对悬挂物体点无法处理。
• 地面模型方法（Ground modeling methods）： 基于2D线提取方法估计局部地面，但此方法不适用于弧形地面。之后有学者在极坐标栅格地图下，使用高斯过程回归方法（Gaussian process regression）解决了这个问题。但是因为空间地面被分割为互相独立的栅格，整张地图的地面无法有效进行估计。
• 基于临近点的地面估计方法（Methods based on the relationship between adjacent points）： 使用局部特征的局部凸标准？（Local Convexity criteria）来进行分割。区别于其他方法的是，此方法保留了全局的点云信息。
• 马尔可夫随机场方法（Markov Random Field，MRF）： 有好几篇研究这个方法的论文，这些论文的方法都利用地面的几何梯度信息来构建 MRF 模型，采用置信传播算法（belief propagation，BP）来对不同区域的环境进行分类。这种方法在其中一片参考文献中展示出对任何地面点云（包括山形地图）下的精确地面分割。但是此方法对计算资源消耗较大，不利于实时处理。

本篇论文受论文 “A hybrid conditional random field for estimating the underlying ground surface from airborne lidar data” 启发，并结合 “Ground segmentation based on loopy belief propagation for sparse 3d point clouds” 、“Graph-based 2d road representation of 3d point clouds for intelligent vehicles” 和 “Drivable Road Detection with 3D Point Clouds Based on the MRF for Intelligent Vehicle” 为重要参考文献研究得出。

 

Ⅱ. 模型（PROPOSED MODEL）

---

A. 概述（Overview）

---

地面被建模为CRF并且在XY平面用2D点阵（lattice）来表示。点阵中每一个节点（Node）$N_i = (n_{x_i},n_{y_i})$与一个连续隐随机变量（continuous hidden random variable） $G_i$相关：

$$G = (G_i, ..., G_n)$$

$$G_i = （h_i, s_{x_i}, s_{y_i}）^T \epsilon \ \mathcal R^3$$

  $h_i$ — 估计高度（estimated elevation）
  $s_{x_i}, s_{y_i}$ — $x$，$y$ 方向的估计坡度（slope）

每个点的测量数据为 $(x_j, y_j, z_j)\ \epsilon \ \mathcal R^3$，节点 $N_i$ 的相关的观察索引点集合为 $\mathcal M_i$。对于点 $j$ 定义一个二进制随机变量（random variable） $C_j \ \epsilon \ \{ 0, 1 \}$，表示该点属于地面点（0）或非地面点（1）。 $Z$ 代表所有点的测量值集合， $C$ 代表所有的与测量值相关的二值标记。 $N_i$ 节点邻域内相关点集合为 $\mathcal N_i$。随机变量 $G_i^{\ t-1}$ 与节点 $N_i$ 有关的上一个时间相关。
目标：根据已经观察到的变量 $Z$ 和上一个时间点的随机分布变量 $G_i^{\ t-1}$，推导出 $G$ 和 $C$。

B. 地面模型（Ground Model）

---

为了对节点间的相互作用进行建模，此处采用了条件随机场（Conditional Random Field， CRF）方法。变量 $G_i$ 与其三个邻域变量相关：空间邻域（spatial neighborhood） $\{ G_k | k \ \epsilon \ \mathcal N_i \}$ 、时间邻域（temporal neighborhood） $G_i^{t-1}$ 以及观测值（observations）和隐藏变量（hidden variables） $\{ Z_j, C_j | j \ \epsilon \ \mathcal M_i \}$。

STCRF示意图如上图所示，图中的各子势函数（clique potentials）的吉布斯分布形式（Gibbs distribution form）定义如下：

• 观测势函数（measurement potentials）：

$$\psi = \frac{1}{W_1} exp(-\sum_{i} \alpha \ \sum_{j\ \epsilon\ \mathcal M_i} c_j||z_j - H_{ij}G_j||^2)$$

$$W_1 = \sum_G exp(-\sum_i \alpha \sum_{j\ \epsilon \ \mathcal M_i} c_j|| z_j - H_{ij}G_j ||^2)$$

$$H_{ij} = [\ \ 1\ \ \ x_i-nx_i\ \ \ y_i-ny_i \ \ ]$$

   $W_1$ —— 配分函数（partition function）
   $\alpha$ —— （标量）观测值的置信度（confidence in the measurement）
  $H_{ij}$ —— 节点和观测值相对水平位置有关的观测值矩阵（measurement matrix）

• 节点平滑势函数（Node smoothness potential）

$$\phi = \frac{1}{W_2} exp(-\sum_{i} \beta \ \sum_{j\ \epsilon\ \mathcal N_i} ||G_i - F_{ij}G_j||^2)$$

$$W_2 = \sum_G exp(-\sum_i \beta \sum_{j\ \epsilon \ \mathcal N_i} || G_i - F_{ij}G_j ||^2)$$

$$F_{ij} = \begin{bmatrix} 1 & nx_j-nx_i & ny_j-ny_i \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix}$$

   $W_2$ —— 配分函数（partition function）
   $\beta$ —— （标量）邻近单元的高度估计值置信度（confidence in the elevation estimates）
  $F_{ij}$ —— 相对位置有关的矩阵（measurement matrix）

• 时间交叉势函数（Temporal interaction potential）

  $$\xi = \frac{1}{W_3} exp(-\sum_{i} \gamma \ ||G_i - Q_iG_j^{t-1}||^2)$$

$$W_3 = \sum_G exp(-\sum_i \gamma || G_i - Q_iG_j^{t-1} ||^2)$$

$$Q_i = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix}$$

   $W_3$ —— 配分函数（partition function）
   $\gamma$ —— （标量）上一个时间节点的高度估计值置信度（confidence in the elevation estimates）
  $Q_i$ —— 上一个时间节点到当前时间节点的转换节点矩阵矩阵（transform matrix）

C. 车辆运动模型（Vehicle motion model）

---

为了得到上一个时间段的节点的 $G_i^{t-1}$，并且将前一帧的栅格地图被转换到当前的相对坐标系下， 本文的转换由车辆的位置和方向变化（displacement and orientation change）计算得出，这些值由创传感器 IMU，GPS 和里程计数据得出。对于每一个节点 $G_i$ 上个时间节点的值，由插值计算得出。在与新探测到的区域相关的部分，上一个时间点的值被设为未定义，取出的区域部分可以简单的给去掉或者储存在一个地图中。（此句翻译不太准，In the areas corresponding to newly-discovered regions, the values at previous time step are set undefined, while taken out areas can simply be forgotten, or stored in a long-term map.）

 

Ⅲ. 预测（INFERENCE）

---

A. 高斯模型和迭代期望最大化方法（Gaussian Model and iterative Expectation Maximization method）

---

求解后验概率（posterior inference）的最大值来对单元变量 $\{ G_i \}_{i=1}^n$ 进行估计。地面高度估计值的最优解值 $\{ \widetilde g_i \}$ 由求解后验概率的最大值 $\{ \widetilde g_i \} \ =\ argmax_{\{ g_i \}}(P(G = \{ g_i \}|C,Z,G^{t-1},\theta))$。
由于求解上述式子的精确解是非常繁琐的，因为意味着，计算过程中对所有可能值（在连续空间中，in continuous space）和途中所有可能的点集（cliques）进行求和计算。
为了简化预测过程，引入高斯模型（Gaussian model）来计算地面的状态参数。每个节点 $G_i$ 的分布由一个均值向量（mean vector） $\widetilde G_i$ 和协方差矩阵(convariance matrix) $M_i$ 来代替进行计算。接下来，为了对方程进行简化，使用了信息向量（information matrix） $X_i$ 和信息矩阵（information matrix） $P_i$ ，它们的定义如下：

$$P_i=M_i^{-1}$$

$$X_i=M_i^{-1} \widetilde G_i$$

然后，为了对全部 $G_i$ 和 $C_i$ 分布进行近似计算，使用迭代期望最大化算法（iterative Expectation-Maximization algorithm，EM）进行求解。在期望步（Expectation step）对点的分类进行概率分布（probability distribution）估计，在最大化步（Maximation step）使用这个分布来对地面状态分布进行估计。这两步可选择进行多次迭代来求解结果。

B. EM 模型初始化（EM initialization）

---

开始时，当汽车开始在栅格地图的边缘开始运动，地面的初始状态是未知的，为了初始化 EM 模型过程，最小的地面状态：高度被设为 0，并且没有任何坡度。协方差矩阵被初始化为大对角系数，使这些初始值对最终的影响降到最小。

C. E-step

---

E-step 计算每个测量值 $z_j$ 为地面点的概率。概率 $c_j=p(C_j=1)$ 的值大小通过对最近节点 $i$ 的状态 $G_i$ 的相关计算得出（由上一步 M-step 迭代得出）。
此概率值由局部地表面与观测值的相对高度差 $dz_j$ 计算得出：

$$dz_j=z_i-H_{ij}X_i$$

$$c_j = \begin{cases} exp(-\frac {dz^2} {2.\delta_{up}^2}), &if \ dz \ \ge 0 \\ exp(-\frac {dz^2} {2.\delta_{down}^2}), &if \ dz \lt 0 \\ \end{cases}$$

E-step 中用来计算观测值属于地面店的概率的似然函数如下图所示。该似然函数只依赖于变量 $dz$ 。相对于 0 点的不对称性使得相对低点和障碍物的底部聚集点更有机会被分类为地面点。

D. M-step

---

M-step 通过测量值、邻域的分类来估计每个节点地面的状态分布。为了对分布进行估计，每个节点的信息向量 $X_i$ 和信息矩阵 $P_i$ 由下列方程计算更新：

$$X_i^{k+1} = \alpha \sum_{j \ \epsilon \ \mathcal M_i} c_j z_j H_{ij} ^ T + \beta \sum_{j \ \epsilon \ \mathcal N_i} F_{ij}^T X_j^k + \gamma Q_i^T X_i^{t-1}$$

$$P_i^{k+1} = \alpha \sum_{j \ \epsilon \ \mathcal M_i} c_j H_{ij}^T H_{ij} + \beta \sum_{j \ \epsilon \ \mathcal N_i} F_{ij}^T P_j^k F_{ij} + \gamma Q_i^T P_i^{t-1} Q_i$$

 

Ⅳ. 实验结果（EXPERIMENTAL RESULTS）

---

A. 实验平台（Experimental Platform）

---

B. 运行时间（Runtime）

---

硬件配置与相关参数如下表所示：

本文算法的计算效率对各种不同 LiDAR 传感器，结果如上表所示。这些数值由中等条件的设备计算得出。在另外一个试验平台上，装备有 Titan X，Tegra X1 和 K1，本文提出的方法能广泛满足实时性的需求。算法的效率被提升对了很大空间，通过使用同一个核心中的 E-steo 与 M-step 间的迭代和 EM step 之后的内部区域的快速同步 GPU 通信。

C. 定性分析（Qualitative Analysis）

---

参数设置：$\ \alpha = 1, \beta = 0.5, \gamma = 0.2, \sigma_{up} = 0.05m, \sigma_{down} = 0.5m$
图中地面点为绿色，非地面点为紫色，红色网格代表估计地面高度，网格内每个点与局部区域内节点的平均高度相关。需要注意的是，对于一段时间内没有数据的区域，平均高度设为0，变化为无线大，系统就能默认这些值是没有价值的。本文图1展示的为最简单的路面情况，即路面为水平面，并且障碍物很高。而上图所示情况就稍微复杂一点：障碍物高度低、边缘平滑，会对地面估计带来不好影响。下图展示了不同传感器探测到，有坡度的山区农村公路。如果Velodyne雷达的点云数据分割有效，那么 IBEO 雷达数据发生有意思的情形： 虽然输入较少数据，由于时空条件的调节，系统仍然能够正确的可行驶区区域和障碍物。

下图中展示了模型的空间分割功能强大的泛化性，特别是在障碍物遮挡干扰情况下对地面的估计。

下图展示了在 KITTI 数据集上本文算法结果。

由上述结果可看出，效果令人信服，并且在这一领域广泛应用的部分打标的数据，将会引领更多的定量评估和比较。

 

Ⅴ. 结论与展望（CONCLUSION AND FUTURE WORK）

---

本篇论文提出了一种基于动态估计的局部地面高度与坡度估计的3D 地面点云分割方法。此方法将地面建模为 STCRF， 将周围环境划分受局部观测和时空条件的，互相相关的高度值单元。同时每个单元的地面高度参数，由一个互联的 EM 算法来进行动态更新。
未来的工作包括一下内容：
将模型中移动的障碍物集成到一起，并使用这些结果自动生成占用网格，重新定义一个智能传感器模型。

以上。

---

参考文献：
1. Ground Estimation and Point Cloud Segmentation using SpatioTemporal Conditional Random Field