机载激光雷达是一种新型的测量技术，是获取高精度数字地面模型DEM(Digital Elevation Model)的重要手段(Liu，2008；马洪超 等，2008；周淑芳 等，2007)，也是建筑物提取和城市3维建模的有效技术途径(Forlani 等，2006)。建筑物点云提取，是进行建筑物3维识别与重建的前提工作，利用LiDAR点云提取建筑物，难点在于判别、分离与建筑物高度、面积相近的植被(Awrangjeb 等，2012)。

目前，从数据源角度，建筑物点云提取分为仅基于机载LiDAR点云的方法、机载LiDAR点云与遥感影像相结合的方法。基于LiDAR点云，已有文献(Cheng 等，2013；李亮 等，2016；张爱武 等，2016)一般是在点云滤波的基础上，利用高程特性剔除车辆、行人点，根据建筑物点与植被点在高程纹理(高程差、高程变化、坡度等)、法向量、回波次数和强度等特征表现的差异性，利用这些特征、特征组合或衍生特征，基于反向迭代数学形态学、阈值分割、随机森林等方法进行点云的分类或建筑物点云提取。还有文献(陈光洲 等，2016；赖旭东和万幼川，2005)，通过点云内插，生成DSM深度影像，利用边缘检测的方法提取建筑物区域，从而可提取建筑物点云。随着机载雷达技术的发展，在获取3维激光脚点的同时，可获取航空影像，众多学者结合航空影像或卫星遥感影像提取LiDAR建筑物点云。融合遥感影像，主要是利用光谱特征差异解决建筑物与植被难以分离的问题。Lodha等人(2007)利用点云的高程(Height)及其衍射特征、回波强度(Return Intensity)和影像灰度(Image Intensity)等特征参数，采用最大期望EM(Expectation-Maximization)算法进行点云分类，提取建筑物点云。文献(Awrangjeb 等，2012；Hartfield 等，2011；Lee 等，2008)也均利用了光谱特征信息，首先根据高度阈值滤除车辆和低矮植被点，然后融合某一波段或是植被指数等光谱信息，将建筑物点云和高大植被点云进行分离。

以上研究对建筑物点云提取都具有一定的意义和效果，但主要是针对规则、完好建筑物的情况。基于机载LiDAR数据的震后建筑物点云提取或震害识别研究相对较少(Elberink 等，2011)。由于往往缺少震前机载LiDAR数据，难以用变化检测的方法(Vetrivel 等，2015)，多数方法仅限于利用面向对象的方法(Ural 等，2011；于海洋 等，2011)。对于场景复杂的震后灾区，窦爱霞等人(2013)分析了基于机载LiDAR识别建筑物震害的难点。受到地震影响，建筑物单体往往可能呈现倒塌、局部倒塌、未倒塌等多种状态(王晓青 等，2015)，建筑物高度、屋顶发生变化。倒塌或局部倒塌建筑物与植被在纹理特征、坡度变化、法向量以及回波次数等方面表现相似的特性，难以继续利用这些点云特征参数进行二者的分离。此外，房屋建筑种类繁多、形状各式各样、材质复杂(水泥、彩钢棚、青瓦等)，这给利用激光回波强度信息带来了困难。

基于以上困难或不足，本文尝试提出了一种融合航空影像光谱特征信息提取震后机载LiDAR建筑物点云的方法。该方法主要特点是借助航空影像光谱特征信息有效剔除难于与建筑物点区分的植被点，探讨性地利用DBSCAN聚类算法对空间低密度、离散的噪声点(车辆、行人点以及未滤除的少量植被点)进行识别，从而最终实现建筑物点云的提取。

融合航空影像的震后机载LiDAR建筑物点云提取方法，首先对点云滤波处理，分离地面点云和非地面点云(或称“地物点云”，包括建筑物、植被、车辆和行人等)，将航空影像光谱特征信息赋予非地面点云；然后，基于阈值分割的方法，剔除植被点；最后，根据车辆、行人点等噪声点分散、密度小的特点，采用基于密度的DBSCAN算法进行聚类，去除密度过小的噪声点，实现建筑物点云的提取。建筑物点云提取主要步骤有：点云滤波、基于阈值的分割和基于DBSCAN聚类的噪声点去除，整个过程如图1所示。

图 1 建筑物点云提取流程图 Figure 1 Flow chart of building point cloud extraction

点云滤波实现地面点云和非地面点云的分离，目前常用的滤波算法主要有基于不规则三角网(TIN)滤波法(Axelsson，2000)、基于坡度/斜率滤波法(Susaki，2012)、数学形态学滤波(Zhang 等，2003)、移动曲面拟合法(苏伟 等，2009；张小红，2007)和一些其他的滤波算法等。本文采用Zhang等人(2016)提出的布料模拟滤波(CSF)方法。该滤波方法具有应用于各种地形效果较好、设置参数个数少等优点。算法基本原理来源于一个简单的物理模拟过程：假设一块足够柔软的布料放置在地形之上，由于重力作用，布料下沉，并逐渐吸附在地行表面，最后布料的形状即为数字表面模型(DSM)；若将地形倒置，具有一定硬度的布料放置其上，布料的形状即为数字地形模型(DTM)，如图2所示。同理，将原始点云倒转，布料覆盖在该倒置点云表面，通过分析布料质点与激光脚点的相互作用，最终确定布料的形状。该物理过程采用“质点—弹簧(Mass-Spring)模型”的“布料模拟(Cloth Simulation)”技术实现。

图 2 CSF基本原理(Zhang等人(2016)) Figure 2 Basic principle of CSF (Zhang et al, 2016)

阈值分割用于剔除非地面点云中的植被点。根据绿色植被吸收红、蓝光，对绿光具有反射作用，在近红外波段范围有强反射作用的光谱特性(林文鹏 等，2010)，对波段间进行数值运算，构建植被指数，用于植被识别或提取。常用且特征明显的归一化植被指数(NDVI)需近红外波段参与运算，而本文航空影像仅有红(R)、绿(G)、蓝(B)共3个波段，无法构建NDVI。因此本文考虑基于单波段灰度值的阈值分割。采集建筑物、植被点云训练样本，提取红、绿、蓝3个波段灰度值，确定训练样本的灰度分布范围，基于核密度估计绘制两者训练样本在各单波段中的密度分布图(图3)。

从图3可看出，建筑物、植被点云训练样本的密度曲线在红、绿、蓝3个波段均有相交，都会出现混淆的情况。统计相交部分混淆点云，按照式(1)计算总错分概率：

图 3 建筑物、植被训练样本灰度值分布图 Figure 3 Grey value distribution of training samples of buildings and vegetation

式中，P(t)为总错分概率，t为分割阈值，θ为建筑物点云占总训练样本的比例，P_B(t)表示分割阈值为t时植被点错分为建筑物点云的概率，P_V(t)表示分割阈值为t时建筑物点错分为植被点云的概率。

最佳波段选择时，取定各波段的t值，以建筑物点和植被点的错分概率之和最小为原则。当选定最佳波段后，分割阈值的确定，以该波段中建筑物点和植被点的错分概率之和最小为原则，即当P(t)最小时对应的t值，其解算表达式为

当且仅当

时，P(t)取得最小值。

DBSCAN聚类算法用于去除车辆、行人和极少的植被点，提取最终的建筑物点。DBSCAN算法(Ester 等，1996) 进行建筑物点云空间密度聚类，其基本思想：点云空间2维平面内组成簇(或称为“类”)的每个激光脚点, 在给定其一定邻域ε范围内，包含的激光脚点个数应大于给定的阈值T_min。该算法具有无需指定簇的数目、能有效处理噪声点和发现任意形状空间聚类的簇等优点。下面介绍DBSCAN聚类，首先给出算法中所用几个参数的定义：

定义1　 ε邻域：给定一个激光脚点，其半径为ε的区域称为该脚点的ε邻域。

定义2　核心脚点：若一个脚点的ε邻域内，包含的脚点个数大于或等于T_min，则该脚点为核心脚点。

定义3　直接密度可达：在点云集合D中，若核心脚点p的ε邻域内有一脚点q，则称脚点q是从脚点p出发直接密度可达的(图4(a))。

定义4　密度可达：在点云集合D中，若存在一个脚点链p₁，p₂， $ \cdots $ ，p_n，p_i∈D(1≤i≤n)，给定ε邻域和T_min，脚点p_i+1是脚点p_i的直接密度可达，则脚点p_i是从脚点p₁的密度可达(图4(b))。

定义5　密度相连：给定邻域ε和半径阈值T_min，若存在脚点o∈D，脚点o到脚点p和脚点q均是密度可达的，则脚点p和脚点q是关于邻域ε和半径阈值T_min密度相连的(图4(c))。

图 4 直接密度可达、密度可达、密度相连 Figure 4 Directly density-reachable, density-reachable, density-connected

DBSCAN算法就是找到从核心脚点出发的最大密度相连脚点，这些最大集合脚点对象形成“簇”，不在簇中的脚点则为噪声点。算法的具体聚类过程如下：

(1) 在点云集合D中，随机选择未被经处理的激光脚点p，脚点p的ε邻域包含的脚点个数大于或等于T_min，即选择的脚点对象为核心脚点。将脚点p的ε邻域内所有脚点形成一个簇，标记为“C”。

(2) 遍历C中所有未处理的激光脚点q，若脚点q也是给定邻域半径ε和阈值T_min条件下的核心脚点，则将未被规到任何簇的脚点q加入C。若脚点q不是核心对象，则将脚点q暂时标记为“噪声”。

(3) 重复步骤(2)，遍历C中未被处理的脚点，直到没有新的激光脚点加入C。

(4) 重复步骤(1)到(3)，直到点云集合D中所有激光脚点均被“处理”、没有脚点加入任何簇时，整个过程结束。最终未被加入簇的脚点，即为噪声脚点。

DBSCAN聚类算法，设置的参数较少，仅有邻域半径ε和最小脚点个数阈值T_min。合适ε和T_min的选取十分重要，ε取值过大，导致很多非建筑物点归到同一个簇中，过小则导致建筑物点被错误地标记为噪声点；T_min取值越大，簇的密度越高，建筑物边缘脚点将视作噪声点被去除，过小则导致大量核心脚点，从而将噪声归入簇。

本文研究区位于海地共和国首都太子港，2010年1月12日发生里氏7.0级地震，震中位于首都太子港以西约25 km处，震源深度约10 km(由USGS提供)，地震造成大量房屋建筑损毁和人员伤亡。实验数据受世界银行全球减灾与灾后恢复机构资助GFDDR (Global Facility for Disaster Deduction and Recovery)，罗切斯特理工学院和ImageCat公司采集完成，获取时间为2010年1月21日，主要包括机载LiDAR点云，平均点云密度约为2 points/m²；高分辨航空影像(仅3个波段)，空间分辨率为0.13 m。本文利用总统府及附近数据进行实验测试，面积为260×307 m²，点云个数为163215个，该区域场景复杂，包含建筑物、高大植被、道路、草坪、车辆以及行人等地物(图5)。

图 5 实验数据 Figure 5 The data of study site

数据预处理主要包括噪声脚点去除和影像匹配两部分。点云去噪是点云滤波、分类的前提工作，是指剔除由空中飞鸟、激光多路径效应、飞行姿态以及扫描系统误差等引起的高程异常点云，包括极高点和极低点，可通过查看点云高程分布直方图来剔除。点云通过最邻近算法内插生成数字表面模型(DSM)，利用DSM进行航空影像的正射校正，从而完成航空影像向点云数据的配准。点云内插，根据一个格网内至少包含一个点云的原则，采样间隔取 $1/\sqrt n $ ，n为单位面积内激光脚点数目。

基于上述方法和流程对建筑物点云进行提取。在数据预处理的基础上，利用CSF算法进行点云滤波，很好地实现了地面点云与非地面点云的分离(图6(a)、图6(b))。

图 6 地震灾区建筑物点云提取 Figure 6 Extraction of building point cloud in earthquake disaster area

在确定最佳波段时，按照式(1)最小为原则，不妨把t的值设为各波段中两密度曲线相交点对应的灰度值，即

式中，t为灰度值，其中f_B(t)表示灰度值为t的建筑物点频率值，f_V(t)表示灰度值为t的植被点频率值。计算各波段建筑物点和植被点的总错分概率分别为：0.92%、2.98%、1.96%，可见基于红波段灰度阈值分割的总错分概率最小、不易混淆，因此选择红波段灰度值作为植被点分离的光谱特征信息参数。确定红波段为最佳波段后，通过统计训练样本的红波段灰度值，根据式(3)计算，当分割阈值为209时，总错分概率达到最小。利用该阈值，基于红波段灰度阈值分割的方法，得到图6(c)所示结果，从图6(c)中可见非地面点云中的植被点被有效地剔除，剩下建筑物点、离散分布在道路上的车辆行人点以及极少的植被点(图6(d))。最后，采用基于密度的DBSCAN算法进行聚类，其主要是对半径ε邻域参数和最小点云个数阈值T_min的设置。根据实验区建筑物面积和点云密度情况，对参数ε在3—8 m、T_min在20—50进行试验。当ε和T_min分别取值为4和35时，漏检概率与虚警概率之和最小(图7)。经过聚类后，高密度的建筑物点几乎不变，离散的低密度噪声点(如道路上的车辆、行人以及未分割去除的植被点)被有效去除(图6(e))。

图 7 参数取值 Figure 7 Parameter setting

为定量评价提取效果，验证该方法对震后建筑物点云提取的正确率，通过与手工提取的建筑物点云结果进行比对。建筑物点云提取误差来源包括正确建筑物点未被准确提取和其他类别点云划分为建筑物点云，二者分别用漏检概率和虚警概率衡量。建筑物点云提取误差统计见表1。

表 1(Table 1)

表 1 建筑物点云提取误差 Table 1 Extraction error of building points 参考建筑物点 提取建筑物点 实际建筑物点 漏检点 漏检概率/% 虚警点/% 虚警概率/% 25140 23106 21216 3924 15.61 15.61 7.52

表 1 建筑物点云提取误差 Table 1 Extraction error of building points

表1表明，漏检概率和虚警概率分别为15.61%、7.52%。漏检主要原因有二：一是建筑物毁坏严重，在进行点云滤波时，把倒塌建筑物点当做地面点滤除；二是阴影区的建筑物点，在基于红波段阈值分割时错分为植被点。虚警点产生主要是局部地区车辆等地物点过于密集，基于DBSCAN算法聚类难以去除。漏检概率和虚警概率均较小，建筑物点云提取总精度达到84.39%，满足一定精度范围的要求。因此，该方法提取震后灾区建筑物点云是比较准确有效的，具有一定的实用性，为震后建筑物点云提取建立了一种可靠的方法。

地震灾区建筑物震害识别是建筑物损失评估的基础，利用机载LiDAR点云数据识别建筑物震害的前提是提取震后建筑物点云。建筑物点云提取难点在于其与植被点云分离，而地震灾区场景更加复杂，建筑物点云，特别是遭受毁坏的建筑物点云与植被点云在法向量、回波次数等特征表现更为相似。许多关于建筑物点云提取的研究多针对规则、完好的建筑物开展；此外，基于机载LiDAR点云识别建筑物震害的研究多仅用面向对象的方法，缺乏基于建筑物点云自身特征识别震害的研究。

因此，本文提出了融合航空遥感数据光谱信息提取震后机载LiDAR建筑物点云的方法。该方法实现简单、流程清晰，为地震灾区建筑物点云提取提供了一种思路。其充分利用了遥感数据资源，通过融合航空影像红波段信息去除难于与建筑物点分离的植被点，基于密度空间聚类的思想剔除低密度噪声点，最终实现灾区建筑物点云的提取，总体精度可达到84.39%。这表明在一定的精度要求范围内，该方法提取地震灾区建筑物点云是可行的。

本文航空影像仅有3个波段，提取建筑物点云的误差来源之一是基于红波段阈值分割时，建筑物点与植被点易混淆。随着激光雷达技术的发展，机载激光雷达遥感可获取多光谱甚至高光谱数据，从而探讨利用NDVI或其他波段组合进行阈值分割，提高分辨植被点云与建筑物点云的能力。此外，DBSCAN算法聚类参数设置有赖于人工试验确定，工作量较大，有待于改进DBSCAN算法使参数自适应确定。建筑物点云提取是建筑物震害识别的前提工作，下一步的研究是根据提取的建筑物点云自身特征进行建筑物震害的识别。