刘 建 兴
(江西省核工业地质调查院,330038,南昌)
机载激光雷达(Airborne Light Detection And Ranging, LiDAR)已经成为构建智慧城市获取三维地理空间数据不可或缺的工具之一,并且在建筑物提取[1]、城市变化检测[2]、城市三维建模[3]等取得了广泛的应用。建筑物作为构建智慧城市最主要、最基本、最重要的组成部分,建筑物三维建模俨然已经成为当前研究的一个热点。建筑物是由多个屋顶面片构成的,如何快速、精确地对建筑物屋顶面分割是三维建模必要的步骤。
近几十年来,许多的学者对屋顶面分割进行了很多深入的研究,他们的研究大致可以分为两类:特征聚类法、模型驱动法。特征聚类法主要方法有区域增长算法及其变种[4]、K-means[5]、DBSCAN[6]等。如Sampath等[5]使用K-means从非地面点中分割建筑物,提取建筑物屋顶面片。周钦坤等[7]通过主成分分析计算点云可靠性指标,利用k-means算法利用法向量进行聚类提取建筑物屋顶面。Xu等[8]提出一种改进DBSCAN的平面分割算法,该方法使用自适应的阈值改进DBSCAN算法,然后使用法向量精细化分割平面。这类方法是首先计算点云的特征,如法向量、曲率、密度等特征,然后根据相似性进行分割数据。这类方法在特征计算准确的时候,能够有着较好的分割结果,但是局限性也正是在这里,很多时候需要依赖点搜索邻域值的设定,同时这类方法对噪声非常敏感。
模型驱动法主要是随机采样一致性(Random Sample Consensus,RANSAC)和Hough变换(Hough Transform,HT)。如艾效夷等[9]使用3D Hough变换从机载LiDAR数据中提取平面特征。章大勇等[10]利用以3D Hough变换基础,加入了参数空间的对偶性,提出了一种基于对偶空间分割的三维Hough 变换算法,进行点云平面地标的提取。李娜等[11]引入角度、距离约束使用RANSAC算法对建筑物立面进行提取。何明等[12]使用形态学滤波改进RANSAC算法从点云中提取平面。周嘉俊等[13]提出了使用高差阈值来改进RANSAC算法来进行建筑物平面的提取。这2种方法利用平面数学模型进行屋顶面片的提取,其对噪声和异常值具有稳健性,可以同时得到优化后的屋顶面片模型参数,但是这类方法对分割参数的设置非常敏感,容易得到伪平面,且算法耗时长。
针对上述方法所存在的问题,本文基于5种不同复杂程度的建筑机载点云数据,提出了一种结合改进区域增长和RANSAC由粗到精的建筑物分割方法。该方法首先基于LRSCPK计算点云法向量,然后利用最小曲率区域增长算法进行屋顶面粗分割,最后利用RANSAC进行小平面的分割和屋顶面的优化。
本文方法主要包括3个步骤,点云预处理、屋顶面粗分割、屋顶面精分割,具体的分割流程如图1所示。LiDAR数据是采用有人机搭载Leica ALS60激光雷达扫描系统获取的某城区点云数据,然后通过Leica CloudPro软件进行点云解算与误差检校。
图1 屋顶面分割流程
预处理的目的是将LiDAR数据中建筑物屋顶点云提取出来。目前国内外针对从机载LiDAR数据中提取建筑物的问题已经进行了大量的研究。本文采用文献[1]提出的反向迭代数学形态学(RIMM)算法进行建筑物点云的提取,该算法相较于著名的点云数据处理软件Terrasoild中传统的建筑物点云提取算法更为方便,该方法不需要对LiDAR数据进行滤波处理,且参数设置更为简单,该方法不需要设置窗口大小阈值,只需要设定一个高程差阈值,后续窗口大小阈值与高程差阈值在迭代过程中自动计算。进行建筑物提取之后,将建筑物划分格网投影至二维图像,存在点云数据的格网对应赋值为1,反之赋值为0,最后生成一张二值图。之后再进行连通性分析,实现建筑物点云的单体化分割。
1.2.1 基于LRSCPK的法向量计算 法向量是描述屋顶面中一个重要的特征,在分割屋顶面中起着至关重要的作用。许多的算法中都需要计算点云的法向量,所以某种程度上来说法向量计算的好坏也会影响到数据处理的精度。目前,应用最广泛的法向量计算方法主要是主成分分析法(principle component analysis,PCA),但是该方法计算出来的法向量在面片交界处并不准确(如图2),所以在进行屋顶面分割的时候,在屋脊线附近数据不能准确划分[14]。
(a)PCA法向量计算
(b)LRSCPK法向量计算图2 法向量估算对比示意图
针对这个问题,本文采用一种具有先验知识的低秩子空间聚类框架的法向量估算方法(low-rank subspace clustering framework with prior-knowledge, LRSCPK)[15],该方法在将主成分分析法估算的法向量作为先验知识进行领域聚类,然后利用引导矩阵(式(1))的低秩子空间聚类将各向异性邻域分割成几个各向同性邻域,并为当前点识别一致的子邻域。最后将拟合一致子邻域的平面的法线作为当前点在尖锐特征附近的法线。
min‖Z‖*+β‖PΩ(Z)‖1+γ‖E‖2,1
(1)
式中:‖·‖*表示奇异值之和,‖·‖2,1表示范数之和,Z是向量矩阵,Ω是引导矩阵(guiding matrix),E是融合约束条件,β和γ是2个参数,文献[15]中设置为1。
1.2.2 最小曲率区域增长法 在估算出点云稳健的法向量之后,本文采用PCL库中区域增长算法(Region Growing,RG)进行建筑物屋顶面的分割。该算法主要思路是对点云按照曲率值进行排序,将最小曲率值作为种子点开始增长,直到没有标记的点为止。详细步骤如下。
1)输入单栋建筑物点云数据及其法向量n和曲率c,设置曲率阈值cth、法向量夹角阈值θth、搜索邻域点数目k,根据曲率从小到大进行排序生成队列VP,将最小曲率值作为初始种子点,初始标记为0。
2)对于每个种子点,若该点标记不为-1,则初始化一个种子点队列SC,利用Kd-tree找到它的相邻点{BC},依次遍历{BC},逐一计算相邻点与种子点的曲率阈值cth和法向量夹角阈值θth,将符合条件的点添加到队列SC中且进行标记。
3)遍历种子点队列SC,重复步骤2)直至队列SC为空。
4)再次从队列VP中查找种子点,重复步骤2)、3)直到遍历完队列VP。
1.3.1 RANSAC分割小面片 在进行区域增长之后,大面积的屋顶面片基本会分割出来,但是还存在一部分小面积的屋顶面片不能分割出来,同时一些变化角度小的屋顶面片区域增长算法也不能识别。针对这两点问题,采用RANSAC算法对上一步分割后的结果进行优化。具体步骤如下。
1)上一节方法分割后数据可以分为两类:分割出的数据集{Ri}和分割后留下的数据Q,遍历数据集{Ri},利用RANSAC进行屋顶面片分割,确定分割出的数据为单个平面数据;反之,则将其分割成多个面片,将提取的面片数据添加到面片集{Pi}中。
2)遍历的面片集{Pi},利用点到平面的距离从数据Q中将符合条件的点添加至该集合。
3)利用RANSAC算法从分割后的数据屋顶面片分割,将小面积的屋顶面片分割出来,将其加入面片集{Pi},直到平面分割完成,输出面片集{Pi}。
1.3.2 屋顶面优化 通过前面的步骤中基本能够将数据中屋顶面片全部分割出来,但是可能分割出来的数量会比真实的屋顶面片数量多,还需要进行屋顶面合并工作。由于上一步已经获取到每个平面的参数方程,可以通过平面的法向量和距离来进行合并屋顶面。具体步骤如下。
1)遍历面片集{Pi},依次计算Pi与集合中其他面片的法向量夹角,如果存在Pj与Pi法向量夹角小于阈值,则进入步骤2。
2)计算面片Pj与Pi之间的距离,如果距离小于阈值,则同时标记面片Pj与Pi,将这2个面片合并成一个面片添加至面片集{Pi}out,。反之则标记面片Pi,将添加到面片集{Pi}out。
3)继续遍历面片集{Pi},直至所有面片标记完成,输出{Pi}out将其作为最终分割结果。
为了验证本文方法的有效性与适用性,本次实验数据选取T字型建筑物、人型建筑、复杂建筑等具有代表性的建筑物(数据详细描述见表1)。本文实验的硬件环境:Intel(R)CoreTMi7-11800H处理器、16G内存,软件环境:Visual Studio2015 C++和PCL1.8.1(Point Cloud Libary)。同时将本文算法分割结果与RANSAC[16]、区域增长[17]进行对比分析。本文以舒敏等[14]人的研究成果,不断调整算法参数进行试验,最终确定参数如下:搜索邻近点数目为20,区域增长法向量夹角阈值θth为0.03°,曲率阈值cth为0.1,屋面优化中点到面的距离为0.1 m。
表1 代表性建筑物相关信息
为了定量地评价屋顶分割结果,本文采用文献[18]的方法与人工分割结果对比,计算分割后的屋顶面的完整性(Comp)、正确性(Corr)和质量(Quality)来评估分割点结果:
(2)
其中:TP(True Positive)是指分割结果中发现的真实平面数值,FN(False Negative)是指分割结果中发现的错误平面数值,FP(False Positive)是指分割结果中未发现的平面数值。
本文了利用3种不同方法对5栋具有代表性的建筑物点云数据进行了分割,其结果如图3所示。图3中原始数据按高程进行着色,其中白色实线为建筑物屋脊线(参考分割线)。
图3 分割结果对比
在屋顶面片法向量夹角较小的建筑物上,如B1对应的分割结果,可以发现区域增长算法可能会存在2个屋顶面无法区分的情况;本文算法与 RANSAC算法却能够较好地分割。在常见的T型建筑物中,3种算法都能够比较好地对建筑物进行分割,但是在细节方面,本文的算法分割效果更好,B2中红色方框部分,由于数据法向量之间的角度偏差较小,区域增长算法和RANSAC算法将其分割为同一屋顶面数据。在规则对称的数据上,如B3对应的分割结果,可以发现区域增长算法和RANSAC算法能够很好地将面积较大的屋顶面分割出来,但是在较小的屋顶面却不能较好地完成分割,本文算法却能够较好地将这部分小面积的屋顶面分割出来,但是对于较小的烟囱,这部分数据量较少,也不能进行分割。在人型建筑物较多且屋顶面片法向量夹角较小的建筑物上,如B4对应分割结果,区域增长算法与本文算法分割效果都较好,但是RANSAC算法会存在分割错误的情况。对于屋顶面多且大小不一的复杂建筑物,如B5对应的分割结果,区域增长算法能够较好地分割出面积较大的屋顶面,对于面积较小的屋顶面则会失败,但是本算法却能够很好地将这一部分分割出来。
总的来说,与区域增长算法相比,在夹角较小的屋顶面上,区域增长算法可能会分割失败,将其识别为同一个平面,如B1和B5所示,由于本算法能够精确计算每个点多法向量,能够精确地分割屋顶面,也能够为每个屋顶面分割得到更多的点。针对小范围的屋顶面数据,还利用RANSAC对区域增长后数据进行了优化,也能够精确到将这部分分割出来。与RANSAC算法相比,RANSAC可能会存在屋顶面竞争的问题,如B3分割后效果所示,由于该算法需要随机采样确定起始平面,这样先分割的屋顶面会将其他屋顶面的数据分割到该平面中,造成局部错误分割,而本方法是先进行区域增长,再利用RANSAC算法对区域增长后的数据进行优化,由于分割后数据一般只存在1个或者2个屋顶面,这时候就能减少许多屋顶面竞争的问题。
为了定量评估不同方法的屋面风格性能,以屋顶面为单元统计了分割的完整性(Comp)、正确性(Corr)和质量(Quality),如表2所示。
表2 精度对比
由表2可知,从整体上来看,本文方法分割相较于区域增长算法和RANSAC算法有着较好的分割效果,完整性、正确性、质量的平均值分别为:100%、94.6%、94.6%。反观区域增长和RANSAC算法在完整性、正确性、质量的平均值数值均不高,分别为:84.9%/84%、70.2%/63.6%、65.2%/60.2%。从细节上来看,本文算法在B1、B3、B4数据上能够达到100%的完整性、正确性、质量,对于B2、B5中存在着点云数目较小的人型烟囱,完整性还是100%,但是正确性和质量分别为77.8%和95.2%,这是数据本身缺陷造成的。区域增长和RANSAC算法在单个数据集上并未有很好的表现。
本文提出了一种结合改进区域增长和RANSAC由粗到精的建筑物分割方法。该方法首先基于LRSCPK计算点云法向量,然后利用最小曲率区域增长算法进行屋顶面粗分割,最后利用RANSAC进行小平面的分割和屋顶面的优化。试验表明,所提出的方法能够有效地分割不同程度的复杂建筑物,且在小面积的屋顶面有着较好的分割效果。本文的方法未能考虑到精分割中RANSAC分割小面片中相邻面片可能会存在面片竞争的情况,这也是下一步研究需要优化改进的地方。同时由于机载LiDAR对建筑物侧面信息的获取能力有限,本文方法未应用于建筑侧面点云数据分割,下一步可考虑通过多源点云数据进行融合,获取空地一体化建筑点云进行建筑物分割。