基于机载LiDAR数据的建筑物轮廓提取

朱依民，田林亚，毕继鑫，林 松

(1. 河海大学地球科学与工程学院，江苏 南京 211100； 2. 浙江华东测绘与安全技术有限公司，浙江 杭州 310014)

建筑物轮廓线精准提取作为建筑物三维模型重建的关键环节，对建立数字城市、智慧城市具有重要现实意义[1]。近年来，机载激光雷达测量技术得到了较快发展，与无人机倾斜影像立体匹配技术相比，机载LiDAR测量能快速采集包含大量城市建筑物三维坐标信息的点云数据，特别是独特的扫描式测量方法使其在获取丰富建筑物屋顶信息方面有着得天独厚的优势[2]，为快速、精准、自动提取建筑物轮廓提供了可能。

利用机载LiDAR点云数据自动提取建筑物轮廓虽受到了一定关注，但相关研究尚处于起步阶段。文献[3]通过规则化数字表面模型分类地面点和非地面点，基于8邻域搜索非地面点云中的建筑物信息，并对获取的建筑物表面点云采用梯度图边界跟踪法提取建筑物轮廓信息，但该方法建筑物轮廓提取精度受nDSM分类结果的影响较大。文献[4]集成边缘与局部信息的活动轮廓模型，将多波段图像边界提取算法用于LiDAR点云数据中的建筑物轮廓提取，该方法提取的建筑物屋顶轮廓精度受限于点云分类精度，且分类误差需要人工判断与修改。文献[5]基于高差偏度平衡滤波算法和多级最小外接矩形算法提取建筑物轮廓，但滤波算法的窗口大小和种子点个数等阀值需要通过反复试验方可确定。文献[6]结合建筑物点云与配准后的影像提取建筑物轮廓，先利用α-shapes算法和投票机制提取建筑物粗边界，然后使用从影像中提取的建筑物边缘信息对提取的粗轮廓进行修正，该方法涉及多种阈值的设定，且自适应程度不高。文献[7]基于RANSAC算法计算建筑物平面最佳模型参数，进而构建建筑物屋顶面并从中提取其轮廓，但该方法需先进行预处理得到高精度的DSM点云数据，且建筑物的屋顶需是平面，不具备普遍适用性。本文针对以上研究存在的诸多问题，对机载LiDAR扫描的点云数据特征进行分析，采用改进的区域生长算法对原始点云数据中的地面点、植被点和建筑物点进行分类，进一步基于三维Hough变换算法提取机载LiDAR建筑物平面点云，最终采用α-shape算法获取建筑物的轮廓信息，形成一整套从机载LiDAR扫描数据提取建筑物轮廓信息的数据处理体系，为机载LiDAR在数字城市、智慧城市建立中的应用提供技术支撑。

1 基于机载LiDAR扫描数据的建筑物轮廓提取方法

1.1 改进的区域生长算法滤波地面点

飞行器在实际扫描时会受到诸如多路径效应、飞鸟等影响而产生一些噪声点[8]，为了避免噪声点对后期点云数据处理产生影响，需先对原始点云数据进行去噪处理。常用的去噪方法有局部平面拟合法[9-10]、频率域法[11]、高程分布直方图法[12]等。由于噪声点一般较少且在高程空间分布上较孤立，因此本文建议实际应用时选择高程分布直方图法对点云数据进行去噪，从而得到较好的去噪效果。

对原始LiDAR点云数据去噪后，根据由低到高的处理原则对点云数据进行地面点逐步滤波，地面点常用滤波方法有数学形态学滤波算法、基于不规则三角网逐渐加密滤波算法及区域生长的聚类分割算法。数学形态学滤波算法需要通过测区内建筑物尺寸确定形态学窗口尺寸，基于不规则三角网逐渐加密滤波算法需要设置最大建筑物尺寸、最大地形变化角度阈值、角度和距离判别阈值，且需对非地面点进行多次滤波。顾及以上两种滤波算法需要设置诸多参数的问题，本文对区域生长聚类分割算法进行改进，在生长地面点过程只需设定生长的高度和坡度，基于改进的区域生长聚类分割算法对地面点进行滤波的具体步骤如下：

(1) 如图1所示，遍历点云数据，找出高程最小的3个点P1、P2和P3作为种子点。

(2) 计算种子点的平均高程P0，并通过3个种子点生成种子平面。

(3) 利用KD树搜索P0点周围一定数量的点，并从中筛选以圆心P0、半径为r的球形邻域内的点。

(4) 计算球形邻域内所有点到种子平面的垂直距离ρi，以及每个点和3个种子点连线与种子平面构成的最大夹角βi。

(5) 将计算的球形邻域内每个点的ρi、βi与设定的距离阈值ρ0、角度阈值β0作比较，将同时小于两个阈值的点放入候选种子点集合中。

(6) 为了继续生长，将初始种子点中任意点归为地面点(此处设置为P1)，储存到地面点集合中，从候选种子点集合中选出最小值ρi所对应的点作为新的种子点P4。

(7) 以P2、P3和P4作为种子点重复步骤(2)—(6)，直至所有点均有所判断，最终得到的地面点集合即为通过区域生长算法生长得到的地面点。

采用区域生长算法滤波地面点后，剩余点云数据中包含的低矮植被与车辆，以及较高的植被和建筑物等非地面点，为了在后期三维Hough变换提取建筑物平面过程减少噪声点对提取结果的影响，以滤波的地面点为基准，将非地面点的高程进行归一化处理，再通过一定的高度阈值将低矮植被点和车辆等非地面点滤除，经过归一化高度阈值分类后的点云数据将只包含高植被点和建筑物点。

1.2 三维Hough变换提取建筑物平面

当将一条直线用参数方程r=xcosθ+ysinθ表示时，可以建立一个二维Hough空间，同理将平面用参数方程r=xsinθcosφ+ysinθsinφ+zcosθ表示，便可建立一个三维Hough空间(r，θ，φ)。其中，r为原点到平面的距离，θ∈(0，π)为平面法向量方向的天顶角，φ∈(-π，π)为平面法向量的坐标方位角。

本文将三维Hough变换算法用于机载LiDAR建筑物平面提取，将所有点云数据视为待检测的一系列空间点，以(r，θ，φ)参数的取值间隔将三维Hough空间分割为紧密相连的三维格网，以一定的θ和φ的取值间隔遍历每个点云计算其对应的r值，当计算的r值落入某个格子内便使该格子的计数加1。对所有待检测点云进行三维Hough变换后，取计数最大的格子所对应的(r，θ，φ)并计算出相应的平面，此平面便为提取的建筑物屋顶平面。

1.3 基于α-shape算法提取建筑物轮廓

α-shape算法是指通过将一个半径固定的圆绕着点集滚动进而构建点集轮廓的边缘点检测算法[13-14]，目前该算法已被用于网格生成、医学图像分析等领域。如图2所示，本文将α-shape算法用于建筑物屋顶面轮廓的边缘提取，并将圆半径α设置为三维Hough变换点集中各点平均间距的2～3倍，滚动圆的圆心坐标为

(1)

(2)

式中，(Ox，Oy)为滚动点的圆心坐标；(x1，y1)与(x2，y2)分别为待判断点P1与P2的平面坐标。

基于α-shape算法提取建筑物轮廓的具体步骤如下：

(1) 将点云数据集合S内的所有点投影到XY平面，从中任取一点P1作为圆心，以2α为半径，搜索圆内所有点并构成新的点集S1，从S1中任取一点P2，由P1、P2构成平面圆，计算平面圆的圆心P0。

(2) 计算点集S1内其他点到点P0的距离，当所有点的距离都大于α时，将P1、P2定义为轮廓点，并转至步骤(4)；当出现点的距离小于α的情况时，则转至步骤(3)。

(3) 对点集S1内下一点重复步骤(1)—(2)，直至S1内所有点均有所判断方可结束。

(4) 对点集S内未被判断的点重复步骤(1)—(3)，直至S内所有点均有所判断方可结束。

2 实例计算

为了验证本文从机载LiDAR扫描数据中提取建筑物轮廓方法的可行性，选取某测区的机载雷达数据进行试验，测区使用RIEGLVQ-1560i机载激光雷达测量系统进行测量，整个测区面积约为549 388 m2，共扫描点云3 823 127个。限于篇幅，本文对整个测区进行分块处理，选取其中一块面积为16 941.9 m2、点云数量为190 338个、点云密度为11.23个/m2的分块测区进行试验。测区高低起伏复杂，含有地面、道路、低植被、中等植被、高大植被、建筑物等丰富的地物，完全满足试验需要。经去噪后的测区点云三维显示如图3所示。

利用改进的区域生长算法对去噪后的点云数据进行地面点生长，根据测区实际情况，设置KD树搜索点的数量为500个，生长的距离阈值ρ0为0.5 m，角度阈值β0为15°。经过区域生长后得到的地面点和非地面点如图4所示。

基于改进的区域生长算法对测区地面点进行滤波后，得到的非地面点包含低矮植被、车辆，以及高大的植被和建筑物点云。为了尽量减少非地面点中的非建筑物点对后期建筑物平面提取的影响，本文以生长得到的地面点为基准，将非地面点的高程进行归一化处理，设置3 m的高度阈值将低矮植被和车辆等非地面点去除，此时非地面点由高植被点云和建筑物点云构成。经过高度阈值滤波处理后的非地面点如图5所示。

利用三维Hough变换算法提取建筑物平面，获取计算平面时所需的参数(r，θ，φ)，再计算各点到平面的距离，并将距离小于阈值所对应的点归于同一平面。需要注意的是该步骤会将某些高于建筑物的高植被点归类至所提取的平面点集中，本文通过判断局部邻域内点的法向量方向的一致性去除三维Hough变换提取平面后所含的高植被噪声点。

经过Hough变换提取的平面中会包含建筑物立面，但在后期提取建筑物轮廓的过程要将建筑物点云投影至XY平面，此时建筑物的立面点会被建筑物屋顶面点云所覆盖。因此，为了提高后期提取建筑物轮廓算法的运算效率，本文通过点的法向量方向信息将建筑物立面点去除，提取的建筑物屋顶面点云数据如图6所示。

将提取的建筑物平面点云数据投影至XY平面，再利用α-shape算法提取建筑物轮廓点云。设置滚动圆半径α为0.5 m得到的建筑物轮廓点云和轮廓线如图7所示。

3 精度评定

为了对建筑物轮廓提取精度进行定量评价，参考文献[15]并结合实际情况,本文提出匹配度、形状相似度及位置精度共3个精度评价指标来评定建筑物轮廓的提取精度。

3.1 匹配度

本文采用图像分类精度评估参数中的质量因子Q作为匹配度的评价指标，即

(3)

式中，A为本文算法提取的位于原始建筑物轮廓点集中的点个数；B为属于本文算法提取出的轮廓点但不属于原始的建筑物轮廓点的个数；C为属于原始的建筑物轮廓点但不属于本文算法提取出的轮廓点的个数。

3.2 形状相似度

通过定义面积差Ad和周长差Pd描述提取的建筑物形状相似度,即

(4)

式中，Se为本文算法提取的建筑物面积；Sr为原始的建筑物面积；Ce为本文算法提取的建筑物周长；Cr为原始的建筑物周长。

3.3 位置精度

选用建筑物轮廓拐点平面坐标的平均距离差Mdd对位置精度进行评价，即

(5)

式中，n为每个建筑物拐点的数量；(xei，yei)为提取的建筑物拐点坐标；(xri，yri)为参考建筑物拐点的坐标。

对整个测区实例(如图8(a)所示)所有建筑物采用上述方法进行建筑物轮廓提取，提取结果如图8(b)所示，同时计算上述3个建筑物轮廓提取评价指标，结果见表1。

表1 建筑物轮廓提取精度评价

分析表1可知，匹配度质量因子Q的平均值达到0.848，说明本文研究的建筑物轮廓提取方法具有较高的提取准确度；面积差Ad和周长差Pd的平均值分别为0.209和0.134，表明本文方法可以较为完整地提取建筑物轮廓形状；分析位置精度的计算结果可知，建筑物轮廓拐点平面坐标的平均距离差Mdd优于0.7 m，表明提取的建筑物和建筑物的实际位置在空间上具有很小的差异性，能够为后期建筑物三维重建提供准确的基础依据；此外，匹配度质量因子、面积差和周长差的标准差均优于0.2，而Mdd的标准差大于0.5，究其原因是本文研究实例中多栋建筑物连接在一起，致使提取建筑物拐点具有一定的难度。因此，表1中的标准差计算结果表明本文方法在提取不同建筑物轮廓时具有稳定性和普遍适用性。

4 结 语

建筑物轮廓的准确提取对建筑物三维重建至关重要，本文在深入分析已有算法优缺点基础上，提出了一种综合改进的区域生长算法、三维Hough变换算法和α-shape算法的建筑物轮廓提取方法，研究了对机载LiDAR点云进行去噪、滤波、平面提取、轮廓提取及精度评定等一系列环节，形成了一套完整的从机载LiDAR扫描数据中提取建筑物轮廓信息的数据处理模式与方法。随着建筑物设计概念的多元化，现有规则建筑物轮廓提取方法已无法满足异形建筑物轮廓提取的需要，对于不规则建筑物的轮廓提取笔者将另文研究。