融合LiDAR点云与正射影像的建筑物图割优化提取方法

杜守基，邹峥嵘，张云生，何 雪，王竞雪

1. 中南大学地球科学与信息物理学院，湖南 长沙 410083； 2. 辽宁工程技术大学测绘与地理科学学院，辽宁 阜新 123000

建筑物是城市空间重要的地理信息要素，从遥感数据中自动提取建筑物对于城市三维建模、灾害评估、数字地图更新等应用具有重要的作用[1-2]。建筑物的自动提取一直是摄影测量与遥感领域一个重要的研究方向，但是由于城市环境复杂，从遥感数据中自动提取建筑物仍然面临许多困难[3]。

近年来，国内外学者对建筑物提取进行了广泛的研究。根据所用数据的不同，可以分为基于影像数据[4-6]、基于LiDAR数据[7-10]，以及联合影像与LiDAR数据[11-15]3类方法。第1类方法主要利用影像的光谱信息。由于影像广泛存在的“同物异谱”、“同谱异物”现象，以及阴影、遮挡等现象，单纯利用影像提取的建筑物有较多的错误，难以实用化[16]。利用LiDAR数据提取建筑物主要利用点云的强度、回波以及几何信息，虽然自动化程度较高，但相较于影像数据，LiDAR数据的空间分辨率较低，提取的建筑物轮廓不够精确。

融合LiDAR点云与影像数据提取建筑物，可以充分利用LiDAR数据精确的三维信息和影像的光谱信息以及分辨率高的优势，二者的融合更加有利于自动提取高精度的建筑物信息。文献[11]提出一种面向对象的航空影像与LiDAR数据融合分类方法，但是该方法依赖于高精度的航空影像分割。文献[12]首先利用光谱信息分离出植被，然后利用高程信息从地面道路和建筑物中分离出建筑物。该方法单独利用光谱信息分离植被，利用高程信息分离建筑物，没有将两种数据进行有效融合。文献[13]采用支持向量机融合光谱与几何特征进行建筑物和植被提取，但是监督分类方法需要大量的标记数据。以上方法主要利用了影像的光谱信息来提高建筑物提取精度，而忽略了影像的高分辨率特征来优化建筑物边缘信息。文献[14]利用超高分辨率影像在LiDAR辅助下提取建筑物轮廓；文献[15]融合LiDAR与正射影像进行建筑物提取和边缘规则化。两种方法都是基于影像上的线特征提取进行边缘信息规则化，但是线特征提取通常不完整，易受阴影区域的影响。针对现有方法存在的问题，基于LiDAR点云与影像数据的特点，本文引入图割算法[17]，提出一种融合几何信息与颜色信息的建筑物提取方法，并利用前后景分割的思想来优化建筑物边缘。在本文方法中，一方面通过图割算法来融合几何与颜色信息，考虑邻域关系以保证提取结果邻域内的一致性；另一方面，通过前后景分割的方法基于影像的颜色信息来优化建筑物边缘，以利用影像分辨率高的优势。本文方法主要流程如图1 所示。

图1 本文方法流程Fig.1 Flowchart of the proposed method

1 融合LiDAR点云与正射影像的建筑物提取

1.1 数据预处理

1.1.1 生成DSM、DTM、nDSM

仪器误差或者航空激光点云获取过程中存在的飞鸟，使得获取的点云中存在粗差点，影响后续的DSM内插以及几何特征计算，首先利用基于统计分析的方法剔除粗差点[18]。在此基础上，采用文献[19]的方法将点云分为地面点和非地面点。然后利用全部点云内插DSM，利用地面点内插DTM，再利用DSM减去DTM获取nDSM。由于本文采用的图割优化过程在DSM格网上实现，因此提出一种保边缘特征的DSM内插方法，主要包括以下几步：

(1) 由点云平均密度估计DSM格网分辨率R。

(2) 含有LiDAR点的格网高程采用其包含的LiDAR点的最大高程值。

(3) 缺少LiDAR点的格网高程通过邻近格网计算。首先将当前格网周围8邻域内的有效格网(含有LiDAR点的格网)的高程按照1.5 m间隔统计直方图。当某一个区间内含有的格网数大于有效格网数的一半时，直接将这一区间格网高程的平均值赋予当前格网，否则利用8邻域格网中所有的LiDAR点，按照与格网中心距离反距离加权计算当前格网高程。

(4) 如果8邻域内没有有效格网，则继续扩大邻域范围，直到内插出有效高程。

为节省后续处理时间，将非地面点云投影到DSM格网上，生成一张非地面掩膜Mng以指示非地面区域，也就是建筑物提取候选区域。在Mng中，“0”指示地面区域，“1”指示非地面区域。

1.1.2 正射影像阴影增强

NDVI常用于区分植被和非植被区域，本文也将NDVI作为区分植被和建筑物的一个特征。但是阴影区域较暗的颜色会对NDVI的计算产生较大影响[20]，对此，首先采用文献[21]中的方法检测阴影区域，然后对阴影区域进行增强。具体做法是利用连通分析将阴影区域聚类，然后对每个聚类的对象分别进行阴影增强，即将影像的RGB转换到HSI颜色空间，然后按照式(1)对S和I通道进行拉伸，再转换到RGB空间，以实现对阴影区域的增强。

(1)

1.2 特征计算

本文主要采用了基于点云的3个几何特征：点云平整度、法向量分布、基于nDSM影像的纹理一致性，以及基于正射影像的颜色特征NDVI。下面对采用的特征进行介绍：

1.2.1 点云平整度

通常情况下，建筑物由多个平面组成，而植被区域表面不规则，因此建筑物上的LiDAR点云具有局部平整性。本文采用协方差矩阵分析来计算点云这一特征[22]。

令NP=pii=1,2,…,n表示非地面点，pi=(xi,yi,zi)为其中一个采样点，Np=pj|pj∈NP是pi的k邻近点集(本文k取经验值为15)。Np的3×3协方差矩阵C3×3定义如下

(2)

通过特征分解计算协方差矩阵C3×3的特征值λ0、λ1、λ2(0≤λ0≤λ1≤λ2)，点pi的平整度fp定义如下

(3)

fp越小，表明pi更可能位于平面上，也就是更可能属于建筑物点。计算完所有非地面点的特征fp后，每一个DSM格网的平整度取格网内所有LiDAR点的fp均值，对于没有对应LiDAR点的格网采用1.1节所提的内插方法进行内插。为保证效率，只对Mng中值为“1”的格网进行内插。

1.2.2 点云法向量分布

式中

(4)

1.2.3 基于nDSM影像的纹理特征

1.2.4 颜色相关特征

影像具有丰富的颜色信息，对于建筑物与植被具有较好的区分度。本文采用NDVI来区分植被与建筑物，定义如下

(5)

式中，IR、R分别为对应颜色波段的值，采用阴影增强后的影像计算。计算完正射影像的NDVI后，将其赋予每个DSM格网。此外，生成一张基于NDVI的植被区域掩膜Mv，将fNDVI>0.3的格网值设为“1”，表示植被区域，其他区域设为“0”表示非植被区域。

以上介绍的4种特征fp、fN、fTH、fNDVI取值范围不同，因此利用Logistic方程进行归一化，定义如下

(6)

式中，x0表示归一化的中心；k控制Logistic方程的深度；k取值对结果影响较小，根据方程曲线的形态将4种特征分别设置为-35、2.0、0.2、-10；x0通过试验进行确定，即单独利用某一个特征进行建筑物提取，通过调整x0值，提取结果最优的为该特征的x0取值，对于4个特征，分别设置为0.06、0.8、18、0.12。

1.3 基于图割算法的建筑物提取

上节介绍的4个特征仅仅描述一个点或者一个格网的特征，并没有考虑上下文关系，因此本文引入图割算法来考虑邻域关系以保证邻域内的一致性。图割算法基本思想是构建一个权重图，每一条边的权表示相应的能量值，然后利用最大流/最小割算法寻找图的最小割[17]。图割算法的目的是通过最小化如式(7)定义的能量函数为每一个结点赋予一个标记lp。

(7)

如式(7)所示，能量函数由数据项和平滑项构成，本文建筑物提取的结果是将DSM格网分为两类，一类是建筑物(building)，其他为非建筑物(non-building)，即lp∈{building,non-building}，因此数据项定义如下：

(8)

式中，λ1、λ2分别为几何和颜色特征的权重，λ11、λ12、λ13分别为几何特征中fp、fN、fTH的权重，且满足λ1+λ2=1、λ11+λ12+λ13=1。对于Mng中值为“0”的格网，直接赋予标记为“non-building”，值为“1”的格网，根据式(8)设置数据项。

由于DSM格网高程值在局部区域内具有较高的一致性，而在地物边缘处则有较明显的差异，尤其是建筑物边缘，因此本文采用格网高程值来构建平滑项。此外，在建筑物附近有高程相近的植被的情况下，可能会出现过度传播现象，易出现植被误检为建筑物或者建筑物被漏检的问题。因此本文在平滑项中引入基于NDVI指数的惩罚因子δ，平滑项设置如下

(9)

式中，hp和hq分别为相邻格网的高程值；常数ε保证分母不为0，可由LiDAR数据的高程精度确定，本文设为0.2 m；δ由植被区域掩膜Mv确定是否给予平滑项，当相邻格网在Mv中的值不同时，给予惩罚，避免分为一类，本文δ取值0.1；平滑项系数β与建筑物提取数据(如建筑物高度、建筑物屋顶复杂度等)有关，本文将其设置为1。

能量函数定义后，采用标准最大流/最小割算法[24]求解能量函数的最小割，然后确定每一个DSM格网的标记lp。

1.4 初始结果优化

上一节图割算法优化后的结果中可能存在一些低矮地物，如灌木、篱笆等，因此进一步采用高程阈值Th和面积阈值Ta去除这些地物，即高程小于Th的格网直接去除，然后采用连通分析将提取结果聚为不同的对象，当对象面积小于Ta时直接去除。本文Th和Ta分别设置为1.5 m和2.5 m2。

2 建筑物边缘优化

经过图割优化后提取的建筑物结果可以视为一张二值图。在影像信息的辅助下，本文引入前后景分割的思想，利用图割算法对建筑物的边缘提取结果进一步优化，能量函数依然采用式(7)，但数据项和平滑项的设置与初始建筑物提取阶段有所区别。首先，将初始检测的建筑物区域设为前景，非建筑物区域设为后景，然后检测初始建筑物区域的边缘像素，并在边缘像素附近划定缓冲区。通过图割优化，将前景和后景的标记传播到缓冲区，从而实现建筑物边缘的精确提取。3个区域的数据项分别设置如下

(10)

平滑项由对应正射影像颜色信息构建，定义如下

Vp,q(lp,lq)=

(11)

式中，R、G、IR分别为对应颜色波段的值；εcol设置为0.1。定义能量函数后，采用标准最大流/最小割算法[24]解算能量函数对缓冲区域的格网重新标记。从能量函数的数据项可以看出，本文将初始提取的建筑物区域以及非建筑物区域设置为固定的值，防止在能量最小化过程中标记发生变化，而主要依赖于平滑项将建筑物区域和非建筑区域的标记传播到缓冲区，从而实现缓冲区的精确标记。在此基础上，进一步利用形态学滤波来消除建筑物边缘的锯齿现象。

3 试验与分析

3.1 试验数据

为了验证本文方法，采用如图2所示的ISPRS Vaihingen的数据进行测试，数据采集自德国的Vaihingen地区，具有正射影像和激光点云数据，正射影像如图2(a)所示，分辨率为0.09 m，点云平均密度为4 pts/m2，内插的DSM如图2(b)所示，覆盖面积约2.0 km2，拥有超过1800栋建筑物。首先采用如图2(c)、(d)、(e)所示的3个区域基准数据测试，正射影像黄色线内为测试区域，参考数据由ISPRS官方提供，其中区域1为形状复杂的密集历史建筑物区域(图2(c))，区域2为树木环绕的高耸居民楼区域(图2(d))，区域3为带有小附属建筑物的居民区(图2(e))。整个数据区域没有官方的参考数据，笔者采用手工标记的形式在正射影像上标记了所有的房屋区域，以用于验证本文方法的有效性。

3.2 试验结果

本文中的λ1、λ2以及λ11、λ12、λ13通过试验确定，分别取λ1=λ2=0.5、λ11=0.25、λ12=0.5、λ13=0.25。对于ISPRS Vaihingen基准数据和全部数据，采用相同的参数设置。

3.2.1 ISPRS Vaihingen基准数据结果

图3列举了区域1边缘优化前后建筑物区域在正射影像上套合的结果。图中右上角的区域为图中心黄色方框内区域放大显示结果。从图中可以看出，经过优化后，边缘更准确，消除了锯齿现象，并在一定程度上恢复了拐角信息。

为了进一步评价本文方法，采用ISPRS官方提供的参考数据进行定量评价。采用的指标包括基于面积(per-area)、目标(per-object)、面积大于50 m2的目标(per-object>50 m2)的完整率(%)，正确率(%)以及质量(%)(见文献[25])。结果如表1所示，与ISPRS测试网站上其他方法结果的对比如表2所示。

表1 基准数据提取结果

表2 与其他方法的结果对比

注：指标取3个区域的平均值，加粗的数值为该列最佳值。

从表1和表2可以看出，本文方法具有一定的优势。对于基于面积的指标，本文方法取得了最优的完整率，质量指标最高。基于目标的指标中，正确率为100%，说明本文方法没有误检目标，且完整率和质量指标也较优异。基于目标面积大于50 m2的指标，全部达到了100%。此外，通过表1可以看出，本文方法对于3个环境差异较大的区域都取得了较好的提取结果，也说明了本文方法的稳定性。基于面积的评价结果如图4所示，黄色为正确提取，红色为误检，蓝色为漏检。从视觉上看，本文方法取得了较好的结果。误检主要来自于建筑物附近的植被区域，主要是因为这些区域处于阴影之中，即使进行了阴影增强，也难以有效区分。漏检建筑物主要为较低矮的建筑物，漏检主要原因为：有些建筑物较矮小，只含有少量激光点而无法被提取，有些建筑物部分被植物遮挡而无法提取。此外，一些屋顶结构复杂的建筑物以及在地面点滤波时错分为地面的建筑物也被漏检。本文优化建筑物边缘信息时，采用的是正射影像，因此，正射影像建筑物边缘纠正的效果也对本文方法的结果产生一定影响。

商圈利用共享的WIFI网络形成与个体的直接联系，进而对商圈APP的下载使用进行推广，通过电子地图、商圈导购和优惠停车等多项便民服务加强APP的用户黏度，为商圈内商铺的广告和活动推广提供优质载体。同时,通过与第三方APP(包括大众点评等团购、外卖平台)的信息进行对接，提升用户体验。

3.2.2 ISPRS Vaihingen全部数据结果

由于ISPRS Vaihingen全区域没有提供参考数据，因此笔者通过人机交互方式进行了标记。对于Vaihingen全区域测试了基于面积(per-area)、基于目标(per-object)两个指标，结果如表3 所示，基于面积的评价结果如图5所示。

表3ISPRSVaihingen全部数据提取结果

Tab.3BuildingextractionresultsforwholeareaofISPRSVaihingen(%)

从表3和图5可以看出，本文方法在全部大范围数据上也表现出了较好的性能，说明本文方法具有一定的适应性。从图5中还可以看出，主要的漏检建筑物面积较小，另外部分建筑物屋顶覆盖了较多的植被也被漏检。除了部分误检是由植被区域导致以外，还有部分误检是由于地面点在点云滤波时被错分为非地面，这部分区域被误检为建筑物。该区域右下角为工厂，有大量货物堆积在空旷地面上，这些货物也被误检为建筑物。

本文方法采用C++实现，运行环境为台式机(Intel i5-3470CPU，主频：3.20 GHz，内存8 GB)，对于3个测试区域的提取时间分别为12、15和14 s，全部区域的提取时间为660 s，其中耗时最多的步骤为点云几何特征计算。总体上看，本文方法具有较好的效率。

图2 ISPRS Vaihingen试验数据Fig.2 Experiment data of ISPRS Vaihingen

图3 建筑物边缘优化前后对比Fig.3 The effect of building boundary optimization

图4 基准数据基于面积的评价结果Fig.4 Per-area evaluation results of benchmark data

图5 ISPRS Vaihingen全部数据基于面积的评价结果Fig.5 Per-area evaluati on result for whole area of ISPRS Vaihingen

4 结 论

本文基于图割算法，提出融合LiDAR点云与正射影像的建筑物自动提取方法。本文方法充分利用LiDAR点云与影像提供的几何以及颜色特征，通过图割算法进行优化，获得了邻域平滑一致的结果。综合DSM和NDVI构建能量函数平滑项，有效防止了建筑物与植被区域的标记过度传播。利用影像分辨率高的特点，基于前后景分割的方法进一步优化建筑物边缘，获得了更加准确的边缘信息。通过ISPRS Vaihingen地区的数据进行测试，并与其他方法进行比较，证实了本文方法可以有效地提取建筑物，精度较高，具有一定的实用推广价值。进一步顾及建筑物边缘的规则特性提高建筑物的检测精度是下一步的研究内容。

参考文献：

[1] KHOSHELHAM K, NARDINOCCHI C, FRONTONI E, et al. Performance Evaluation of Automated Approaches to Building Detection in Multi-source Aerial Data[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2010, 65(1): 123-133.

[2] TOMLJENOVIC I, HÖFLE B, TIEDE D, et al. Building Extraction from Airborne Laser Scanning Data: An Analysis of the State of the Art[J]. Remote Sensing, 2015, 7(4): 3826-3862.

[4] MANNO-KOVACS A, SZIRANYI T. Orientation-selective Building Detection in Aerial Images[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2015, 108(1): 94-112.

[5] HUANG Xin, ZHANG Liangpei. Morphological Building/Shadow Index for Building Extraction from High-resolution Imagery over Urban Areas[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2012, 5(1): 161-172.

[6] 胡荣明, 黄小兵, 黄远程. 增强形态学建筑物指数应用于高分辨率遥感影像中建筑物提取[J]. 测绘学报, 2014, 43(5): 514-520. DOI: 10.13485/j.cnki.11-2089.2014.0084.

HU Rongming, HUANG Xiaobing, HUANG Yuancheng. An Enhanced Morphological Building Index for Building Extraction from High-resolution Images[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(5): 514-520. DOI: 10.13485/j.cnki.11-2089.2014.0084.

[7] ZHANG Jixian, LIN Xiangguo, NING Xiaogang. SVM-based Classification of Segmented Airborne LiDAR Point Clouds in Urban Areas[J]. Remote Sensing, 2013, 5(8): 3749-3775.

[8] 徐文学, 杨必胜, 魏征, 等. 多标记点过程的LiDAR点云数据建筑物和树冠提取[J]. 测绘学报, 2013, 42(1): 51-58.

XU Wenxue, YANG Bisheng, WEI Zheng, et al. Building and Tree Crown Extraction from LiDAR Point Cloud Data Based on Multi-marked Point Process[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(1): 51-58.

[9] 徐文学, 杨必胜, 董震, 等. 标记点过程用于点云建筑物提取[J]. 武汉大学学报(信息科学版), 2014, 39(5): 520-525.

XU Wenxue, YANG Bisheng, DONG Zhen, et al. Building Extraction from Point Cloud Using Marked Point Process[J]. Geomatics and Information Science of Wuhan University, 2014, 39(5): 520-525.

[11] 管海燕, 邓非, 张剑清, 等. 面向对象的航空影像与LiDAR数据融合分类[J]. 武汉大学学报(信息科学版), 2009, 34(7): 830-833.

GUAN Haiyan, DENG Fei, ZHANG Jianqing, et al. Object-based Fusion and Classification of Airborne Laser Scanning Data and Aerial Images[J]. Geomatics and Information Science of Wuhan University, 2009, 34(7): 830-833.

[12] 程效军, 程小龙, 胡敏捷, 等. 融合航空影像和LIDAR点云的建筑物探测及轮廓提取[J]. 中国激光, 2016, 43(5): 0514002.

CHENG Xiaojun, CHENG Xiaolong, HU Minjie, et al. Buildings Detection and Contour Extraction by Fusion of Aerial Images and LiDAR Point Cloud[J]. Chinese Journal of Lasers, 2016, 43(5): 0514002.

[13] ZAREA A, MOHAMMADZADEH A. A Novel Building and Tree Detection Method from LiDAR Data and Aerial Images[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2016, 9(5): 1864-1875.

[14] 程亮, 龚健雅. LiDAR辅助下利用超高分辨率影像提取建筑物轮廓方法[J]. 测绘学报, 2008, 37(3): 391-393, 399. DOI: 10.3321/j.issn:1001-1595.2008.03.021.

CHENG Liang, GONG Jianya. Building Boundary Extraction Using Very High Resolution Images and LiDAR[J]. Acta Geodaetica et Cartographica Sinica, 2008, 37(3): 391-393, 399. DOI: 10.3321/j.issn:1001-1595.2008.03.021.

[15] GILANI S A N, AWRANGJEB M, LU Guojun. An Automatic Building Extraction and Regularization Technique Using LiDAR Point Cloud Data and Orthoimage[J]. Remote Sensing, 2016, 8(3): 258.

[16] AWRANGJEB M, FRASER C S. Automatic Segmentation of Raw LiDAR Data for Extraction of Building Roofs[J]. Remote Sensing, 2014, 6(5): 3716-3751.

[17] BOYKOV Y Y, JOLLY M P. Interactive Graph Cuts for Optimal Boundary & Region Segmentation of Objects in N-D Images[C]∥Proceedings of the 8th International Conference on Computer Vision. Vancouver: IEEE, 2001: 105-112.

[18] FERRARI S, FERRIGNO G, PIURI V, et al. Reducing and Filtering Point Clouds with Enhanced Vector Quantization[J]. IEEE Transactions on Neural Networks, 2007, 18(1): 161-177.

[19] HU Han, DING Yulin, ZHU Qing, et al. An Adaptive Surface Filter for Airborne Laser Scanning Point Clouds by Means of Regularization and Bending Energy[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, 92: 98-111.

[20] QIN Rongjun, FANG Wei. A Hierarchical Building Detection Method for Very High Resolution Remotely Sensed Images Combined with DSM Using Graph Cut Optimization[J]. Photogrammetric Engineering & Remote Sensing, 2014, 80(9): 873-883.

[21] 帅滔, 张洪艳, 张良培. 面向对象的高分辨遥感影像阴影探测方法[J]. 光子学报, 2015, 44(12): 1228002.

SHUAI Tao, ZHANG Hongyan, ZHANG Liangpei. The Object-based Method of Shadow Detection in High-resolution Remote Sensing Imagery[J]. Acta Photonica Sinica, 2015, 44(12): 1228002.

[22] ZHOU Qianyi, NEUMANN U. Fast and Extensible Building Modeling from Airborne LiDAR Data[C]∥Proceedings of the 16th ACM SIGSPATIAL International Conference on Advances in Geographic Information Systems. Irvine: ACM, 2008: 1-8.

[23] HARALICK R M, SHANMUGAM K, DINSTEIN I. Textural Features for Image Classification[J]. IEEE Transactions on Systems, Man, and Cybernetics, 1973, SMC-3(6): 610-621.

[24] BOYKOV Y, KOLMOGOROV V. An Experimental Comparison of Min-cut/Max-flow Algorithms for Energy Minimization in Vision[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2004, 26(9): 1124-1137.

[25] ROTTENSTEINER F, SOHN G, GERKE M, et al. Results of the ISPRS Benchmark on Urban Object Detection and 3D Building Reconstruction[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2014, 93: 256-271.