基于机载LiDAR数据的建筑物屋顶边界提取

朱军桃,王 雷,赵苗兴,龚朝飞,罗 乐

(桂林理工大学 a.测绘地理信息学院;b.广西空间信息与测绘重点实验室,广西 桂林 541006)

随着机载激光雷达测量技术的不断完善, 其广阔应用前景受到人们的极大关注, 机载LiDAR技术已渗入到生产、 生活的各个领域。 数字地球及数字城市建设产生了构建建筑物三维模型的极大需求, 而对于构建建筑物三维模型来说, 需要对建筑物边界进行精确提取,才能提高三维建模的精度。现今,对建筑物边界提取的方法主要分为两类: 1)结合其他数据源与点云数据对建筑物边界进行提取[1-3](此类方法对数据类型要求高)； 2)仅使用LiDAR点云数据进行建筑物边界提取。 李江雄[4]通过将散乱点云网格化, 获得边界网孔并求得每个边界网孔中的最小凸包来提取点云边界, 该方法对点云密度要求高, 对机载LiDAR数据限制性强;曾齐红等[5]通过构建TIN模型来判断点云边界, 该方法由于受到凸包性限制, 可能使边缘点提取不完整; 沈蔚等[6]采用Alpha Shape算法提取建筑物边界, 该方法可以避免因凸包性造成的影响, 但判别时要遍历所有激光点, 同时要求解圆心, 导致程序运行效率低; 符小俐等[7]采用带距离控制的卷包裹算法提取建筑物边缘, 此方法对距离控制要求很高, 且提取边精度不高。本文提出了一种基于机载点云数据提取建筑物屋顶边界的方法，与基于不规则三角网模型的边缘提取方法相比，获得了更多边缘点，且避免了由于凸包性导致对建筑物边缘提取造成关键点遗漏问题，对多种形状的建筑物适应性强，运行效率高。

1 原理与方法

现实中,绝大多数房屋都属于精心设计过的、转角为直角的规则或不规则多边形。针对此特点,本文提出了一种直接采用点云数据提取房屋屋顶边线的方法。 首先,以细小网格划分点云,从包围网的4个方向对指定房屋进行边缘点初步提取,其结果基本反映了建筑物的具体形状;然后,采用改进的带距离控制的卷包裹算法提取建筑物角点，并依据相邻点之间的关系二次提取边缘点,提高提取精度;最后,对边缘点进行最小二乘直线拟合、规则化及扩展，重现了建筑物边界线,具体流程见图1。

图1 基于机载LiDAR点云数据的房屋轮廓线提取流程Fig.1 Extraction process of building outlines based on airborne LiDAR data

1.1 边界初步提取

将屋顶点云数据投影至XOY平面, 对屋顶面投影的二维点集构建囊括点云簇的四边形(大于最小外包四边形); 以边长接近点云平均间距的正方形网格将四边形分割, 分割网格能分割所有点云，并使每个网格内具有1～2个激光点。 至此, 每个网格近似代表了一个激光点, 并且相邻网格间没有间隙, 使该分割结果确保对初始边界搜索时遍历所有的点云。 四边形内部的网格分为两类: 一类是内部有激光点，称为“实孔”; 另一类是内部没有激光点，称为“空孔”。 以四边形某条边的最外侧网孔开始, 对垂直于四边形边界的某一列网孔顺序判断: 如果第一个网孔为空孔，则将空孔计数器赋值为1, 如果下一个网孔仍为空孔，计数器增加1, 按序逐一判断； 如果标记空孔的计数器超过某一阈值, 且遇到下一个网孔为实孔则认为其内部的点为边缘点, 并将计数器赋值为0; 若连续出现空孔数量并未超过阈值,不记录下次出现实孔内的点,并将计数器重新赋值为0。按照此方法,从4个方向逐次判断所有网孔。如图2所示, “●”为按箭头方向搜索到的边缘点。由于点云分布并不完全均匀,若内部出现较大间距,可能将属于建筑物内部的非边缘点标记为建筑物的边缘点,所以在所有方向的边缘点提取完成之后,删除距离其他边缘点较远的孤立点和数量极少的点簇。

图2 边缘点初步提取Fig.2 Preliminary extraction of edge points

1.2 角点提取

在初步提取边缘点云后,由于点云分布的离散性会导致初步提取的相邻边缘点间出现严重锯齿现象,该现象导致传统依赖相邻激光点之间夹角来判断边缘角点的方法失效。带距离控制的卷包裹算法是从某一边界上的固定点开始，并以一定长度旋转,触碰到的第一个点必是外包多边形的一个顶点,并以此点为固定点继续旋转,直到回到起始点。方法步骤具体为[8]: 1)在点集中依次找出X和Y坐标最大点与最小点，取其中一点作为出发点，记为P1； 2)从P1出发作一条射线，顺时针或逆时针方向旋转，直至触碰到点集中某点，记为P2； 3)重复步骤2)，若触碰点回到P1，则结束。

上述方法获取的边缘点密度较低且无法提取凸包边缘。 本文改进了带距离控制的卷包裹算法, 并采用该算法进行角点提取, 过程如下: ① 顺序连接1.1节中提取的边缘点, 将边缘点所围成图形分为内外两侧; ② 在初步提取的边缘点中寻找距离最大的两个点, 认为这两个点一定属于边界上的点; ③ 以其中一点P开始，沿点云外侧定义一条方向边,并以此边开始以一定长度L旋转, 在长度L范围内碰到第一个点P1时停止旋转并记录P1; ④ 以P与P1所连直线为方向边,再次以P1为中心，长度L为半径旋转,碰到第一个点P2时停止旋转,并记录夹角α。当Pn-1、Pn、Pn+1所构成的夹角αn∈180°±θ时(θ为允许角度变化的范围),Pn为边线上的点; 当αn∈ 270°±θ时,Pn为角点(图3); 当αn∉(180°±θ)∩αn∉(270°±θ)时, 将长度L增加D(D为所增加的长度), 再次旋转, 若触碰点改变, 则此点为边线上的点,以此处为起始继续步骤③(Pn+2点处); ⑤ 重复步骤④, 直到回到起始点P处; ⑥ 以点P开始从点云内侧定义一条方向边,继续步骤③～⑤。

卷包裹算法只能判断具有“凸”性质的角点,而对内“凹”性质的角点无法判断。因初步提取的建筑物边缘点已形成了闭合环状点云,从内侧进行改进算法是将凹角点转化为凸角点的处理过程。同时,步骤④中对长度L增加D，避免了从内侧旋转过程中角点的错误判断。

1.3 边缘点二次提取

初步提取中得到的边缘点呈现锯齿状主要原因是:(1)点云数据本身的特性,由于激光点云是对物体表面的一种不规则分布点位置的描述,所以造成锯齿现象与多种因素有关,处于最外侧网孔的激光点未必都为边缘点;(2)由于网孔大小和点云密度分布不均的影响,导致每个最外侧网孔中的激光点未必只有一个，所以需要对边缘点进行二次提取。依据建筑物的相邻边缘点间所连线段某侧不存在激光点的性质,提取相邻角点间的边缘点,此边缘点集合同属于一条房屋边线上。二次提取方法为:将两个角点存入集合list中,获得两角点连线的斜率K,向边缘点外侧任意距离构建一条斜率为K的直线S,将此相邻角点及角点间的边缘点投影至S上,设角点A投影点为A′,距其最近的投影点B′,次近的投影点为C′，分别计算与投影点所对应的角点A,边缘点B、C到直线S的距离DIS:

图3 角点提取示意图Fig.3 Sketch map of corner points extraction

若DISA≥DISB且DISB

若DISA≥DISB且DISB>DISC,将B存入集合list中,并以C点作为下次判断起始;

若DISADISC,剔除B点,并以C点作为下次判断起始点;

若DISA

重复上述方法判断两角点间所有边缘点,并根据边缘点的粗糙程度及数量将集合list中的边缘点设置合适的迭代次数进行迭代。此方法是将两角点连线模拟为房屋边界线,逐步剔除不属于边缘点的激光点,采用上述方法可以极大地弱化两个角点所连直线与真正房屋边界线斜率的差异所造成同一边线上边缘点错误提取的问题。

图4为某建筑屋顶的边缘点提取结果图,其中图4a为投影到XOY平面的建筑物点云, 图4b为初步提取的边缘点集, 图4c为二次提取的边缘点集。

可以看出， 经二次提取后其边界粗糙程度明显降低。

1.4 房屋边线拟合、规则化及扩展

针对经过二次提取后的边缘点集, 采用最小二乘方法将角点间的边缘点拟合为直线。 所拟合的边界线虽然已经能够很好地反映出建筑物的具体形状, 但相邻的边界线仍存在不垂直的关系, 为了更好地反映屋顶的真实形状, 将各相邻边界线规则化为直角, 并对规则化后的边线进行扩展处理。 本文采用文献[9]所述的扩展方法, 计算某规则化后的边线与参与此边线拟合的外侧边缘点的距离, 并将此边平移至距离最大处,最大程度地还原建筑物的原本大小及形状。

图4 边缘点提取结果Fig.4 Extraction results of edge points

2 实验分析

由于实验区域地势平坦且建筑物较高, 因此屋顶点云的获取直接根据点云的高程分布[7],高程分布密度最大区域为地面点云所在集合, 去除此部分及以下点云数据得到建筑物点云, 使用CloudCompare软件剔除墙面点, 最终得到建筑物屋顶点云。

为保证获得外侧边缘的“实孔”,所构外包四边形边长为Lt=Lmin t+Gli×Gle+1,其中:Lt为外包四边形边长,Lmin t为最小外包四边形边长,Gli为空孔阈值,Gle为网孔边长。

实验参数:编译工具为Matlab R2014b，点云平均密度为3.4个/m2，网孔边长0.5 m，空孔阈值为6。

图5为提取某地区高程230～260 m、260～308 m、308～325 m具有内外边界的不规则多边形(图5a)、 内转角为直角的不规则多边形(图5b)及规则四边形(图5c)3种最具代表性建筑物。左侧房屋内边界中曲线部分仍按本文方法进行二次提取,但依然完好地保留边界特征。在拟合相邻角点间的直线时,对于距离非常近的角点不进行拟合,图5b中圆孔非房屋内边界,因此未进行边界拟合。

对各二次提取设置不同迭代次数后所得到结果作比较:

(1)

式中: (xi,yi)为属于同一边线上边缘点集合中的点;f(xi)为xi处经最小二乘拟合的值;Δ值反映边缘点在所拟合直线两侧分布的离散程度。

表1各边号对应图5中所在边,可以看出绝大多数边在增加迭代次数后Δ值会降低,但少量数据存在异常：4、9、13、19、28、32、35号边Δ值在降低到一定程度就不再降低,造成这种情况的原因是边缘点集形成了一种光滑的、不再具有锯齿状的形状; 2、 14、 18、 21、 25、 27、 39号边Δ值在降低到一定程度后增大了, 造成这种情况的原因是设置太高的迭代次数过滤掉过多的边缘点。 上述两种情况均出现在一些长度较短的边上, 此类边缘增加迭代次数后最终因太少的边缘点影响了拟合效果, 同时因扫描线方向、 点云分布不均等原因造成初始迭代的Δ值各不相同。 经试验, 本文增加迭代次数找到最低Δ值, 并增加限制条件:边缘点最少不能小于0.5个/m, 该条件在保证足够数量的边缘点基础上达到最佳提取效果。 最终, 规则化及扩展后的屋顶边线如图6所示。

将使用不规则三角网(triangulated irregular network，TIN)模型所提取的边缘点(图7)与本文方法提取的边缘点(图8)作比较:由于受凸包性限制,TIN所构成的三角形可能会掩盖其他边缘点,同时对于属于凹角的边线转角处所构成的三角形会掩盖转角点,如图7所示,其中各凹角处均有不同程度的边缘点损失,而采用本文方法可以有效避免了此缺陷。

表2比较了3种不同轮廓形状屋顶用不同方法提取屋顶边缘点的效果,可以看出，采用本文方法对各类房屋的提取效果均优于基于构网的提取效果。

图5 最小二乘拟合屋顶边线Fig.5 Building outlines by Least Square Fitting method

表1 不同迭代次数比较Table 1 Comparison of different iterations

图6 规则化及扩展后的屋顶边线Fig.6 Building outlines after regularization and expansion

图7 构网提取的边缘点Fig.7 Edge points by delaunay triangles

图8 本文方法提取的边缘点Fig.8 Edge points through the algorithm in this paper

表2 两种方法提取效果对比
Table 2 Comparison of extraction results of two methods

房屋形状点云数量算法边缘点数量相同点数量重复率/% 转角为直角规则四边形3 900本文构网144968588.54 转角为直角不规则多边形9 429本文构网31517514683.43 具有内外边界的不规则多边形6 125本文构网36620119295.52

表3对比两种方法程序运行时间,采取本文方法对各形状建筑物边界屋顶提取所花费时间均少于构网， 所以采取本文提取方法可以高效且充分的提取建筑物边缘点。

表3 两种方法运行时间对比Table 3 Operation time comparison of two methods

为充分验证本文提取结果的有效性,采用文献[10]所述方法对本文提取结果予以验证,以质量因子Q、形状相似度因子Ra(面积差)、Rp(周长差)及位置精度因子Dtrc为验证指标,定量评价本文提取效果。

(2)

(3)

(4)

其中:TP为正确分割激光点数量;FN为漏提激光点数量;FP为错误分割激光点数量;Ae为本文提取方法所获面积;Ar为标准数据所获面积;Pe为本文提取方法所获周长;Pr为标准数据所获周长;Xe、Ye为本文方法提取的建筑物中心坐标;Xr、Yr为标准数据所获建筑物中心坐标。最终精度结果如表4。

表4 精度评定Table 4 Precision assessment

3 结 论

本文提出了一种基于机载LiDAR点云数据建筑物屋顶边界提取方法,相较于基于不规则三角网算法提取边缘点的繁琐构网过程,该方法运行效率更高,并且规避了凸包性缺陷对边缘点及角点提取的影响。经试验分析,该方法对各类型建筑物轮廓提取效果均优于构建不规则三角网的提取算法,所提取建筑物的边缘较为完整,有较强适应性。

采用的边缘点初步提取环节根据简单的提取规则,有效地获得了边缘的原始点集,为二次提取的运行压缩了数据量;针对初步提取的边缘点进行的角点提取并非遍历所有边缘点,提高算法运行效率;二次提取提高边缘点提取的精确度;最终采用最小二乘法拟合出房屋屋顶边线,并对拟合出的直线进行规则化及扩展处理,重现了房屋边界线,为建筑物边界提取提出一种新的思路。本文方法仍有不足之处:若初步提取的边缘点过于粗糙,则可能造成提取的角点异常;同时,本文算法在处理弧形边界的建筑物时适用性较差。笔者将在后续的研究中逐步完善此类问题。