点云数据配准作为地面激光雷达(light detection and ranging，LiDAR)点云数据处理中的重要环节，在现实工作中仍然存在人工干预多、效率低的问题^[1]，实现点云数据快速自动配准是目前地面LiDAR领域研究的热点^[2-4]。在实际工程中，人们较多地采用标靶(如：球面标靶^{[5, 6]}、平面标靶^{[7, 8]}等)辅助点云数据配准，其中球面标靶以各向同性、配准精度高^{[9, 10]}、配置灵活等特点获得了广泛应用，它以同名球面标靶的球心点作为同名点，弥补了因扫描仪分辨率、遮挡等因素较难形成同名点的缺陷^[10]。

目前，球面标靶探测方法主要有基于曲率的方法、解析法、分割法和特殊材质法。蓝秋萍等^[6]根据曲率和法向确定球状目标可能存在的位置，该方法理论上可识别多种半径的球面标靶，但是曲率、法向计算易受噪声影响^[11]；Franaszek等^[10]在点云数据集中三点滤波的基础上确定候选球面中心点，根据球体与测站几何关系估算球心位置、确定候选球面点，拟合球体参数，该方法理论简单明了，易于实施，但是该方法以整个数据集作为探测范围，存在效率低下，易受扫描分辨率、噪声等因素影响；石宏斌等^[12]提出了一种多几何基元的点云数据分割方法，因而也可用于球面标靶探测。但是，为获得有限个球面标靶，要先将整个场景分割，存在时间开销大，效率低的问题；张立朔等^[13]根据球面标靶的特殊材质所形成的特殊反射特性，根据经验值确定候选标靶区域，利用圆度和周长确定球面标靶区域。在地面激光三维扫描中，由于受距离、入射角等因素的影响，激光点的强度值存在“曝光”不一致的问题^[14]，同一材质不同的入射角、扫描距离，其强度值差异较大。因此，直接根据经验值提取球面标靶，缺少理论支持，算法通用性不强。针对球面标靶检测过程中自动化程度低、效率低和易受噪声影响等问题，本文在分析地面LiDAR数据采集特点的基础上，提出了一种基于遮挡边界检测的球面标靶快速探测方法。

1 地面LiDAR扫描中的球面标靶

地面LiDAR以预设的经差、纬差做穹形阵列式扫描，因而单站数据具有近似的栅格结构^[15]。利用球面标靶辅助点云数据配准，要求扫描仪与球面标靶之间是通视的，球面标靶的放置客观上造成场景的局部不连续，产生遮挡与被遮挡边界(图 1(a)和图 1(b))。由于球面标靶形体较小和扫描经差、纬差近似相等，遮挡边界在单站栅格结构中形状为近圆形(图 1(c))。在经过球体的扫描线中，根据扫描线与球面标靶空间位置关系，可将它们的关系分为相交和相切两种(图 2(a)和图 2(b))。假设球体的起始和终止扫描线均与球体相切(图 2(a))，则可建立图 2(c)的几何关系。

图 2中，S为扫描仪中心(为坐标原点)，O为球体中心(坐标为x_o, y_o, z_o)，R为球体半径，SA和SB分别与球体相切于A、B两点。因此，有sinγ=R/l(其中$l=|S O|=\sqrt{x_{o}^{2}+y_{o}^{2}+z_{o}^{2}}$，为扫描仪中心到球心的距离)，由于扫描范围限制及球面标靶半径较小，因此$R \ll l$，所以有γ=sinγ。根据这一几何关系，完成球体扫描需要2γ/δ条扫描线(δ为相邻扫描线间夹角)。在二维栅格结构中，球体所经过的扫描线数目即为二维栅格结构中近圆的直径估值，因此圆半径估值r_est=γ/δ，对于图 2(b)中的相交情况，其扫描线数目比相切少1；另外，球面标靶在二维栅格结构中的近圆区域，是由球面点填充而成，球面点数目也就是它的面积，因此，可用近圆直径的估值来估算近圆区域的面积，也就是球面点数目。从图 2可得出以下结论：

1)$r \in\left[r_{\min }, r_{\max }\right]$，其中$r_{\min }=\left\lfloor r_{\mathrm{est}}-1\left(r_{\min }>0\right)\right.$，${{r}_{\max }}=\left\ulcorner {{r}_{\text{est}}} \right\rceil +1$；

2) 球面点数估值$S_{e s t} \in\left[S_{\min }, S_{\max }\right]$，其中${S_{\min }} = \left\lfloor {{\rm{ \mathsf{ π} }}{r_{{{\min }^2}}}, {S_{\max }} = \ulcorner {\rm{ \mathsf{ π} }}{r_{{{\max }^2}}}} \right\rceil $。

因此，本文球面标靶检测的方法可描述为以下过程：

(1) 遮挡边界检测。在单站栅格数据中，检测各有效点的八邻域是否存在距离突变，存在距离突变则视为遮挡边界点；

(2) 遮挡边界空间聚类。根据遮挡边界点到扫描仪中心距离近似相等原理和遮挡边界点的空间位置分布，以较小空间阈值对遮挡边界进行空间聚类；

(3) 球体参数探测。在各候选聚类结点中利用随机采样一致性(random sample consensus, RANSAC)方法探测边界中的圆形边界，标示球体可能的存在区域，在各圆形区域所对应的三维点集中探测球体模型，利用几何关系约束获得最终球体参数模型。

2 球面标靶自动检测 2.1 遮挡边界检测

在地面LiDAR数据采集中，遮挡与被遮挡是最为常见的现象，它一般存在于不同对象或对象内部不同部件间，具有深度不连续的特点，可分为遮挡边界和被遮挡边界，其中位于遮挡产生处的前置对象边界称为遮挡边界，位于后置对象边界称为被遮挡边界(如图 1(b)所示)。利用数据集的栅格结构，计算当前非空点p_{m, n}=(x, y, z)到其八邻域N={p_j=(x_j, y_j, z_j)|0≤j < 8}的深度差d_j=|p_{m, n}|－|p_j|，当d_j大于所设定的阈值d_threshold，则p_{m, n}为被遮挡点，当d_j < 0且|d_j|>d_threshold，p_{m, n}为遮挡点。本文认为遮挡产生的地方即是球面标靶可能存在的地方。因此设d_threshold=R_default(R_default为预设球面标靶的半径)。

2.2 遮挡边界空间聚类

由于球面标靶各向同性，位于球面标靶的遮挡边界点到扫描仪中心的距离近似相等。因此，可以以较小的距离阈值对遮挡边界进行空间聚类，以加快后续的圆检测。由于单站数据为具有二维栅格特性的三维点集，所提取的边界点仍然具有该特性。因此，遮挡边界空间聚类方法可表述为：

(1) 标识所有边界点为待分类点；

(2) 在栅格结构中选取一待分类点e_{m, n}，由该点创建聚类结点c_i，改变该点标识为分类点；

(3) 计算该点八邻域中待分类点p_i到扫描仪的距离与结点内第一个边界点到扫描仪的距离的差值，小于所设定阈值d_th1，则将p_i加入c_i，并改变其分类标识；

(4) 对c_i中所有成员执行步骤(3)，对c_i扩展，直到c_i成员不再增加。则聚类集合C=C∪c_i；

(5) 对当前所有待分类点执行步骤(2)至步骤(4)，直到待分类点数目为0。由于球面标靶各向同性和对称原理，位于球面的遮挡边界点到扫描仪中心距离近似相等，其差值小于$2\delta \sqrt{x^{2}+y^{2}+z^{2}}$，因此设$d_{{th1}}=2 \delta \sqrt{x^{2}+y^{2}+z^{2}}$。

2.3 球面标靶探测

在遮挡边界各聚类结点中，首先在二维空间中进行圆模型检测，根据§1中结论1)滤除不合理模型，获得候选圆模型集合。在各候选圆区域对应的三维点集中探测球体模型。两个模型检测均为在含有噪声的数据集中探测预设数学模型，鉴于RANSAC方法的低内存消耗、易扩展等特性^[12]，本文的模型检测采用RANSAC方法进行。

1) 圆模型检测。在当前结点中随机抽取三个非共线点p₁，p₂，p₃构建三角形Δp₁p₂p₃及其外接圆M_{(o, r)}=(x－x₀)²+(y－y₀)²=r²作为假设模型(x₀, y₀为M_{(o, r)}在二维栅格结构中的圆心o的坐标，r为圆的半径)，计算各二维边界点p_i到圆M_{(o, r)}的距离d_M，若满足式(1)，则加入M_{(o, r)}的一致集P_M，以P_M的数目作为M_{(o, r)}的打分S_M，以迭代的方式获得当前结点最优模型，半径满足§1中结论1)的，进行球体模型探测。式(1)中d_th2为边界点作为圆模型一致集的距离允许误差。

$ d_{M}=| | p_{i}-o|-r| \leqslant d_{t h 2} $

(1)

2) 球体模型检测。通过半径约束的M_{(o, r)}，在其圆形区域对应的三维点集中探测球体模型，以文献[12]中的方法构建球体模型，即首先对三维点集计算法向量，然后随机抽取两个采样点p₁，p₂(其法向为n₁，n₂)，以(p₁，n₁)、(p₂，n₂)分别确定两个空间直线l₁，l₂，以l₁，l₂间最短垂线段中点作为球心o，以(|op₁+op₂|)/2为半径R，确定球模型M_{(o, R)}，计算三维点p到球面距离D_sph及其法向n与该点在球面投影处法向n′间的夹角θ，设定距离允许误差d_th、法向夹角允许偏差θ_th，满足式(2)为模型M_{(o, R)}的一致集，以一致集数目作为M_{(o, R)}的打分S_M，以迭代的方式确定当前三维点集中的最优球体模型M_{(o, R)}及其一致集。然后对一致集利用最小二乘法重新拟合球体模型参数M_{(o, R)}，若M_{(o, R)}中半径R满足式(3)且S_M∈[S_min, S_max]，则将M_{(o, R)}作为球面标靶所对应的球体模型。

$ D_{\mathrm{sph}} <d_{t h}, \theta=a \cos \left(n \times n^{\prime}\right) \leqslant \theta_{t h} $

(2)

$ \left|R-R_{\text { default }}\right| \leqslant R_{t h} $

(3)

式中，R_th为所探测到的球体半径与预设球面标靶半径之间的允许误差。

3 球面标靶检测实验

为验证本文算法，选择扫描仪R EGL VZ-400作为数据采集设备，对某建筑物作全景扫描(5 056×2 151)，扫描距离约为6 m，测距误差约为±2 mm，扫描线夹角近似为0.071 2°，场景中放置4个标称半径为0.072 m的球面标靶，标靶到扫描中心的平均距离为3 m，获得数据的强度反射图像如图 3所示。算法采用C++语言，以VS.NET2013为开发平台，运行平台为Windows 7 64bit专业版，内存为16 GB，CPU i7-4790，主频3.6 GHz。

3.1 球面标靶探测

首先对实验数据进行遮挡边界检测，距离阈值d_threshold=R_default=0.072 m，边界检测效果如图 4(a)所示，从结果看，本文方法能够有效提取扫描场景中普通存在的遮挡边界。在这些遮挡边界当中，由球面标靶遮挡产生的边界比例较小，应用§2.2中的方法对遮挡边界进行空间聚类，聚类结果如图 4(b)所示，由图 4可以看出，该聚类能够有效将球面标靶产生的遮挡边界与其他对象产生的边界区分开，从而将RANSAC方法的采样样本限制在局部空间，从而大大减少圆模型、球体模型检测过程中的时间开销。在各聚类集合中，利用RANSAC方法分别进行了圆、球体模型的检测，这一过程中，RANSAC方法的参数设置如表 1所示，利用§1结论1)和结论2)分别对圆、球模型检测结果进行约束，球体半径满足式(3)的，作为最终的球体模型参数，圆模型和球体模型探测结果如图 4(c)和图 4(d)所示，球体模型对应的探测过程参数如表 2所示(其中Pt(o)为模型M_{(o, r)}的圆心所对应的三维点)。

3.2 球面标靶探测效率

在3.1球面标靶检测过程中，共检测到464 206个遮挡边界点、24 603个聚类结点、126个圆形区域和4套球面标靶参数，共耗时16 411 ms，其中遮挡边界点检测796 ms，边界点空间聚类219 ms，圆检测11 119 ms，球体检测4 277 ms。实际工作中，在千万级点云中由人工识别这样的4个球面标靶，约需10 min时间。因此本文方法能够实现球面标靶的快速自动识别。

4 结束语

本文针对球面标靶检测中自动化程度低、效率低和易受噪声影响等问题，提出了基于遮挡边界检测的球面标靶自动检测方法，实验结果表明，本文所提方法具有较高的效率，适合于具有阵列式扫描的站载激光点云数据，在当前工程应用当中，具有很大的实用价值。