机载激光雷达数据的建筑物三维重建方法
【专利摘要】本发明专利——机载激光雷达数据的建筑物三维重建方法的目的是提出一种建筑屋顶边界与屋顶拓扑图相结合的三维重建方法，从而实现建筑物的自动检测和三维重建。首先，对机载LiDAR数据滤波处理得到地面点和非地面点，再结合点云特征信息从非地面点云中提取建筑物点云，继而分割屋顶面、提取边界轮廓线，最后结合建筑物边界与屋顶拓扑图，得到屋顶关键线段，构造各屋顶面的封闭多边形及其组合，从而得到建筑屋顶模型。墙面又可由DTM或地面点的高程信息得到，从而实现建筑物3D模型重建。本专利一定程度上突破了现阶段机载LiDAR数据的应用瓶颈，降低重建过程的复杂性和提高三维重建的灵活性，为城市三维重建提供了一个突破口。
【专利说明】机载激光雷达数据的建筑物三维重建方法 所属技术领域
[0001] 本发明专利属于一种数字摄影测量技术，涉及一种从机载激光雷达数据提取建筑 物点云，并进行建筑物三维模型重建，此方法是一项具有应用价值的方法，建筑物数字模型 是数字城市建设中的重要组成部分。
【背景技术】
[0002] 随着城市的快速发展，许多城市相继推出适合自身发展的"数字城市建设"发展战 略，其中建筑物三维数字重建是其中的重要内容。3D建筑模型在数字城市中是非常重要的 组成部分，在其他应用方面也是如此，如地理信息系统、城市规划、灾害管理、应急响应以及 虚拟/增强现实等。机载LiDAR(Airborne Light Detection And Ranging)技术可以直接、 快速获取建筑顶部高精度、高密度的三维空间信息，是当前城市建筑三维数字重建重要的 数据来源。
[0003] 但是LiDAR数据处理以及3D重建相关算法的发展却滞后于LiDAR数据获取设备的 开发。事实上，对于大范围的建筑自动建模一直都是一项困难而且耗时的工作，尤其结构复 杂的建筑物3D重建是一个具有挑战性的问题，3D重建已经成为高质量点云数据广泛应用的 瓶颈。在过去二十年里，虽然许多学者已经在这一方向做了很多工作，而且也取得了骄人的 成果，但是研究对象一直局限于相对简单的建筑。许多复杂建筑的重建只能由人工或者人 机交互实现。一些常见的简单建筑能够利用现有的方法实现其自动化建模。然而对于结构 复杂或者LiDAR数据不完整的建筑仍难以重建，这也一直是目前的研究热点。由于大部分建 筑比较复杂，城市建筑建模仍然昂贵和费时，能自动处理高质量数据的算法在应用方面一 直很有需求。
【发明内容】
[0004] 本发明专利即利用机载激光雷达数据进行建筑物三维建模，首先利用机载LiDAR 数据特点并结合地物点云特征（回波、几何空间分布等），基于一种层进式方法分离出建筑 物点云;然后基于法向量聚类分割方法或者RANSAC分割算法对建筑物点云进行屋顶面片分 害J;最后基于建筑物边界提供的墙面信息，利用屋顶拓扑图（即屋顶面邻接关系）提取建筑 屋顶关键线段(屋脊线和主要边界线），按照连接规则构造出屋顶面封闭多边形，再结合DTM 或者地面提供的高程信息构造建筑物墙面，通过多边形的组合即可得到完整的建筑物三维 模型。
【附图说明】
[0005] 下面结合附图和实施例对本发明专利进一步说明。
[0006] 图1是基于机载LiDAR数据的建筑物模型重建技术方案
[0007] 图2是建筑物提取流程
[0008] 图3是角度和距离约束
[0009] 图4是α-shape算法提取不意图
[0010] 图5是边界连通性分析 [0011]图6是规则化原则 [0012]图7是屋顶拓扑图 [00 13]图8是最小闭合环搜索
[0014]图9是水平屋脊线(蓝色)拓展到边界(红色）
[0015]图10是水平屋脊线(蓝色)拓展至边界延长线 [0016]图11是屋脊线拓展 [0017] 图12是建筑物模型
【具体实施方式】
[0018]本发明专利总体思路如图1。首先，对机载LiDAR数据滤波处理得到地面点和非地 面点，再结合点云特征信息从非地面点云中提取建筑物点云;其次，针对提取出的建筑屋顶 点云进行屋顶面分割以及边界轮廓线的提取;最后，根据建筑物边界与屋顶拓扑图的结合 得到的屋顶关键线段构造出每个屋顶面的封闭多边形，通过多边形的组合即可得到建筑屋 顶模型。墙面又可由DTM或者地面点的高程信息构造得到，从而实现建筑物3D模型重建。
[0019] (1)建筑物点云提取
[0020] 基于机载LiDAR数据特征建立了一种层进式的建筑点云提取方法（图2)，首先对原 始LiDAR数据采用渐进式形态学滤波算法得到地面点和非地面点，然后结合点云回波次数、 点云法向量等特征逐步从非地面点中提取出初始建筑物区域，最后结合建筑物高程和面积 等几何信息进行精提取得到建筑物点云。
[0021] 首先，对原始点云进行滤波。渐进式形态学滤波过程如下：首先对LiDAR点云进行 格网化，根据点云坐标以及格网尺寸将其分配到对应网格，并对空网格(无 LiDAR点云)进行 内插处理，最终生成格网数据;然后以渐变结构窗口 ^对格网数据进行开运算，当开运算前 后的网格高程之差大于高差阈值时，判定该网格为非地面类网格，否则为地面类网格;最后 通过改变滤波窗口大小进行多次迭代计算，在迭代过程中逐渐增大窗口尺寸(按照线性或 指数形式），最终分离地面点与非地面点。具体算法步骤如下：
L〇〇23J 其中高程差阈值Δ hi计算公式如式（1): dh0 w, < 3
[0024] < s (wi -- ) c + dh0 ^. > 3 (1) 滿醒 吨> dK似'
[0025]式中，dho是初始高差阈值;dhmax为最大高差阈值，一般为最矮建筑物的高度;c为 网格大小；S为区域的平均地形坡度;Wi为第i次（i = 1，2，3，…，M)的窗口大小，Μ为建筑物最 大尺寸，Wl可表示为式(2):
[0026] wi = 21+l (2)
[0027] 然后，基于点云回波次数进行植被点检测。由于激光具有穿透和折射现象，当机载 激光雷达系统扫描地物时会记录点云的回波次数信息，该特征有助于地物的检测。一般植 被以多次回波为主，而建筑物仅有一次回波。据此，利用LiDAR点云的回波次数信息可以检 测出大部分植被点云。但是其中会包含部分建筑物边界点，这会破坏建筑物点云的完整性， 最后导致生成的建筑物模型与真实模型的差异性变大。基于此问题，从地物的同质性出发， 根据回波次数信息计算每个激光点的多回波特性F?:
(3:)
[0029] 式中，？1表示激光点邻域内回波次数为i的点;N为邻域内点的个数；ΣνΡιμ表示该 点邻域内回波次数大于1的点的个数；Σν Ρι表示该点邻域内点的总数。所以上式也可以描述 为激光点为植被点的概率，Fmr越接近于1，则该点为植被点的可能性越大，反之，为非植被点 的可能性越大。
[0030] 其次，利用点云法向量进行建筑物点云粗提取。对于点P，其邻域点集 A = {iV，内，…，Λ ]可以定义为：
[0031] pf-P^<D (4)
[0032] 式中，D表示邻域的最大距离，k表示该点邻域内点的数目。这里选用的邻域点搜索 算法是最邻近点搜索算法(ANN)，该算法采用KD-tree实现，能够快速的进行邻近点搜索。它 有两种途径，一种是k-搜索算法，即搜索与查询点最邻近的k个点；另一种是r-搜索算法，即 搜索与查询点距离在半径r之内的所有k个点。K-搜索算法旨在固定搜索一定数量的邻近 点，有可能将不相关的局外点加入到邻域点集，导致该点集不一定能够真实体现局部区域 特征，据此类点集计算出的点云法向量误差会较大。而r-搜索算法搜索的是固定半径范围 内的邻近点，更能够代表局部特征，故采用该算法进行邻近点搜索，确保点云法向量计算结 果更优。
[0033] 估算曲面上某点法向量的问题可以转化为求取该点的切平面法线的问题。目前常 用的求解方法有最小二乘法、特征值法以及主成分分析法。主成分分析法相较于其它两种 方法，核心思想是将变量从高维空间降到低维空间处理，选取少数最具代表性的不相关变 量来揭示原来变量所包含的信息。本专利对某点邻域内所有点组成的协方差矩阵进行特征 值分解来估算该点的法向量的过程，其本质上是对点云集合进行主成分分析。
[0034] 对于某一点p，其邻域内所有点构造的协方差矩阵Cov为： (5)
[0036] 式中，k表示该点邻域内点的总个数;pj表示邻域内的第j个点；表示该点邻域内 所有点的重心。利用主成分分析法对点P邻域点云进行分析，便可得到其特征值λο,λ^λ^λο 彡入丨彡λ 2)和对应的特征向量eo，ei，e2。三个特征向量反映了点云分布的三个主方向，且eo与 edPe2正交，代表了拟合平面的法线方向，因此把最小的特征值λ〇对应的特征向量 eQ作为该 点的法向量。
[0037] 在对LiDAR点云进行区域增长时，种子点的选取极为重要，关系到增长结果的好 坏。由于曲率能够表示表面的变化程度，因此本专利根据目标点的表面曲率大小来进行种 子点的选取。主成分分析除了可以用于表面法线的估计，也可用于目标点的表面曲率的推 算，利用其解算协方差矩阵Cov得到的最小特征值λ〇可以近似为目标点邻域的曲面变化度。 目标点在该点邻域内沿曲面法线eo的变化度〇定义为：
(6)
[0039] 式中，人〇山，\2(\(^\1彡\ 2)为协方差矩阵〇^的的特征值。目标点曲率值〇越小表 示邻域内的点位于曲面切平面上的概率越大，因此选取曲率比较小的的点作为生长种子 点。具体区域增长算法步骤如下：
[0040] 1)选取未处理点云中曲率小于设定阈值的点作为种子点，并添加种子点集中。
[0041] 2)搜索种子点一定范围内的邻近点云，计算每个邻近点与种子点的法向量夹角， 如果两者的夹角小于设定的阈值，则将该邻近点和种子点归为同质区域。
[0042] 3)比较邻近点的曲率和阈值大小，若小于则将其添加到种子点集中，同时将检测 过的种子点删除。
[0043] 4)对种子点集中每个点重复上述过程2和3。如果点集中的点全部处理完毕，则该 区域增长结束，然后统计同质区域点云数目并判断是否符合最小数目阈值要求。
[0044]对未处理的点云重复上述操作，即可完成建筑物点云的粗提取。
[0045]最后，基于粗提取的结果进行精提取。粗提取的LiDAR点云主要包括建筑物点和部 分植被点。为提取出建筑物点，本专利采用连通成分分析对初始建筑物点云进行欧式聚类。 然后采用反距离加权法对地面点进行空间插值生成DTM。最后结合几何形状(例如面积、高 差)等特征进一步区分建筑物点和植被点，从而提取出建筑物点云。
[0046] (2)建筑屋顶分割
[0047]主要采用两种分割算法:基于点云空间分布特性的聚类增长分割算法和RANSAC算 法。针对利用层进式方法提取出的建筑物点云，两种算法的思路分别如下：
[0048] a、基于点云空间分布特性的聚类增长分割算法:首先利用建筑屋顶平面具有点云 法向量相似性的特点进行聚类得到初始分割结果，然后利用点到面的距离和点云密度将未 分割点正确分配到对应的屋顶面点集，最后结合面片优化策略完成屋顶的准确分割。
[0049]基于法向量的聚类增长算法步骤如下：
[0050]
[0051] b、RANSAC算法:首先随机选择三个点并计算由他们决定的平面方程参数，然后根 据点到该平面的距离，统计满足距离阈值的点云数目，如此重复执行N次，每次的结果和上 一次比较，保存最佳结果(点云数目最大），最得到的就是最佳平面。
[0052] RANSAC算法具体步骤如下：
[0055] 针对上述两种算法的分割结果，建立了如下优策略：①两个面片之间的法向量夹 角小于一定阈值;②两个面片之间的距离小于一定阈值。如图3所示，pl 2、pl2为两个待优化 面片，Pl、p2分别为面片的重心，分别为过点P1、P^法线，g为向量，?平行 于.磧。
[0056] 策略①中，两面片之间的夹角可以用它们的平面法向量夹角Θ表示：
[0057] #= cos 1 ） (7)
[0058] 对于策略②，常用的方法是计算两平面间的距离，即Ad= IcU-cbl，然后通过Ad与 阈值的比较来进行判断。但是计算得到的面片夹角Θ存在一定的误差，并且夹角误差会随着 原点到平面距离的增大而被放大，导致Ad值波动较大，因此很难确定合适的阈值。为此将 平面间距离定义为： / r ~-r· r · V
[0059] = maxl rv?nx , rv?n2 I (8)
[0060] 设置好合适的角度阈值Θ和距离阈值d后，当满足上述两个策略时，将两个面片合 并，即可得到完整的屋顶分割结果。
[0061 ] (3)建筑物模型构建
[0062]结合建筑物边界建提出了一种改进的基于拓扑图的建筑屋顶模型重建方法。首先 根据屋顶面之间的拓扑关系（即邻接关系)对相关屋顶面求交得到屋脊线，再在拓扑图中寻 找最小闭合环来统一关联屋脊线的端点。又建筑物边界代表了建筑物墙面的位置，则屋脊 线的其它端点可通过屋脊线与相邻墙面相交进行调整，同时和每条屋脊线相关联的两个屋 顶面与对应墙面相交又可以得到相关屋顶面的其它边界线。对于没有相交关系的单个屋顶 面，则可以求其外接多边形。在利用拓扑图求得屋顶面的主要边界线后，根据连接规则就可 以得到每个屋顶面的封闭多边形，最后进行多边形组合即可得到整个建筑屋顶模型。上述 过程可以分为两个主要步骤:建筑物边界提取和基于拓扑图的建筑物模型构建。
[0063] (A)建筑物边界提取
[0064]首先对建筑物点云利用α-shape算法（图4)二维平面内进行边界点提取;然后对边 界点集利用RANSAC算法进行边界线段提取，同时进行连通性分析(如图5);最后依据边界规 则化法则(如图6)完成建筑物边界的提取。
[0065] α-shape算法具体用于提取边界点的步骤如下：
[0066] 1)在点集S中选取一点Pi，并搜索其半径2α范围内的空间邻域点集Q，然后选取Q内 任一点Ρ2，并根据圆半径α求过这两点圆的圆心Ρ〇。
[0067] 2)计算点集Q中每个点（除Ρ^ΡΜ)与Ρο的距离L，如果所有的L都大于α，则判定Pl p2都是边界点;反之，有一个L小于〇则停止判断并跳到步骤(3)。
[0068] 3)对Q中下一个点重复过程1)和2)，直至Q中所有点都完成判断。
[0069] 4)选取S中下一个点重复过程1)~3)，当S中所有点完成判断则结束。
[0070] 当目标点集经过上述流程处理后，建筑边界点就被提取出来。为更好表现建筑物 边界范围，提取出的边界点已被投影到水平面。实验表明，当α半径设置为为点云间距1~2 倍时，α-shape算法能完整的提取出建筑物边界。
[0071 ] 在该算法中，圆心Po(xo，y())可以根据Pi(xi，yi)、P2(X2，y2)和半径α计算得到：
m (10)
[0074] 式中，Dpip2为Pi到P2的欧式距尚。
[0075] (B)基于拓扑图的建筑物模型重建
[0076]首先基于拓扑图（图7)提供的屋顶面邻接关系计算得到屋脊线;然后在屋顶拓扑 图中进行最小闭合环搜索（图8)，将闭合环中对应屋脊线的相关端点统一到某一交点;其次 基于屋脊线所提供的相关屋顶面，并结合邻近的建筑物边界线所构造的的墙面，利用线面 相交将屋脊线拓展到边界，利用面面相交完成屋顶面其它边界线段的提取，实现屋脊线的 拓展（图9-图11);最后将屋顶关键线段(屋脊线和主要屋顶面边界线)按照一定规则连接成 封闭多边形，再结合DTM或者地面点提供的高程信息完成建筑物墙面构建，多边形与墙面的 组合即可完成建筑物三维模型重建（图12)。
[0077]利用两个面片LiDAR点云的位置来判断交线的存在(也即拓扑关系的存在）以及交 线端点，判断准则如下：
[0078] 1)如果两面片在交线缓冲区内都存在LiDAR点云则可确定交线存在。缓冲区距离 是根据各自交线来确定，而不是根据整个建筑点云。该距离是一个和两面片中值点间距有 关系的函数:取两面片中点密度较小的面片的中值点间距的两倍。
[0079] 2)选择落在根据1)确定的缓冲区的LiDAR点云来确定交线的长度。交线的端点则 根据相交面片落在缓冲区的LiDAR点云到交线的投影来确定。对于每个面片保留最外层的 投影点位置，这样就可以得到四个投影点。如果两个面片对应的投影部分有重叠，则处在中 间的两个投影点作为交线的端点。交线的最小长度可根据LiDAR点云密度计算得到。因为由 高密度数据得到的交线比由低密度数据得到的交线能更好的代表短的屋脊线，所以再次取 中值点间距的两倍进行计算。在自动化处理过程中通过分析数据来设置参数值极为重要， 因为自适应过程能够减少非自适应(一倍、两倍或者三倍)带来的影响。最后把较长的线段 作为最终的屋脊线。
[0080] 基于屋顶分割面片的边界点云来确定缓冲区点云，进而确定屋脊线的存在及其长 度。
[0081 ]最小闭合环搜索算法如下：
【主权项】
1. 基于机载激光雷达数据的建筑物三维重建方法，其特征主要包括以下几个步骤： (1) 基于原始机载LiDAR数据提取建筑物点云； (2) 建筑屋顶分割； (3) 建筑物模型构建。2. 根据权利要求1所述的方法，其特征在于:利用一种层进式提取方法实现建筑屋顶点 云的快速提取。3. 根据权利要求1所述的方法，其特征在于:利用一种建筑物边界与屋顶拓扑图结合的 方法完成建筑物三维模型重建。
【文档编号】G06T7/00GK106097311SQ201610378366
【公开日】2016年11月9日
【申请日】2016年5月31日
【发明人】王成, 李亮, 习晓环
【申请人】中国科学院遥感与数字地球研究所