EN中文

Journal of Geo-information Science >

2023 , Vol. 25 >Issue 8: 1531 - 1545

DOI: https://doi.org/10.12082/dqxxkx.2023.220927

Complex Roof Structure Reconstruction by 3D Primitive Fitting from Point Clouds

  • ZHANG Wenyuan ,
  • CHEN Jiangyuan , * ,
  • TAN Guoxin
Expand
  • National Research Center of Cultural Industries, Central China Normal University, Wuhan 430079, China
*CHEN Jiangyuan, E-mail:

Received date: 2022-11-28

  Revised date: 2023-04-15

  Online published: 2023-07-14

Supported by

National Natural Science Foundation of China(41801295)

National Culture and Tourism Science and Technology Innovation Project(2019-008)

Fold

Abstract

Geometric and semantic integration of 3D building models are important infrastructure data for smart city, they are conducive for promoting the refined management and intelligent application of building facilities. However, most of the existing point cloud-based modeling methods focus on the reconstruction of geometric models with simple roof structure, and semantic and topological relations are ignored. Moreover, these methods are sensitive to noise, which are difficult to assure topological consistency and geometric accuracy. To solve these problems, this paper proposes a 3D primitive fitting algorithm for automatically reconstructing building models with complex roof structure from point clouds. Firstly, a 3D building primitive library is designed, including various 3D building primitives with simple and complex roof types. Secondly, an individual building point cloud input is segmented into multiple planes using RANSAC algorithm. The Roof Topology Graph (RTG) is then generated according to the relationship of roof planes, and the roof type of point cloud is subsequently recognized by comparison of RTG between point cloud and building primitives. Thirdly, the reconstruction is formulated as an optimization problem that minimizes the Point-to-Mesh Distance (PMD) between the point cloud and the candidate meshed building primitive. The sequential quadratic programming optimization algorithm with necessary constraints is adopted to perform holistically primitive fitting, so as to estimate the shape and position parameters of a 3D primitive. Finally, the parameterized model is automatically converted into City Geography Markup Language (CityGML) building model based on the prior 3D building primitive. The generated CityGML LoD2 (second level of detail) models are different from mesh models created by conventional building modeling methods, which are represented with geometric, semantic, and topological information. To evaluate the quality and performance of the proposed approach, airborne lidar and photogrammetric building point clouds with different roof structures are collected from public datasets for test. Several building models with complex roof types are successfully reconstructed by using this approach, and the average PMD of five models is 0.17 m. The proposed algorithm is also compared with three other methods. Experimental results indicate that the proposed method achieves the best geometric accuracy, because the average PMD of each model is less than that of other methods. Moreover, this automatic primitive fitting method is efficient, and it is also robust to noise and local data missing. This study demonstrates that the resulting building models can well fit the input point cloud with topologic integrity and rich semantic. This method provides great potential for accurate and rapid reconstruction of geometric-semantic coherent building models with complex roof condition.

Key words： building; complex roof; point cloud; 3D reconstruction; 3D primitive; CityGML; semantic; model-driven

Cite this article

ZHANG Wenyuan , CHEN Jiangyuan , TAN Guoxin . Complex Roof Structure Reconstruction by 3D Primitive Fitting from Point Clouds[J]. Journal of Geo-information Science, 2023 , 25(8) : 1531 -1545 . DOI: 10.12082/dqxxkx.2023.220927

1 引言

建筑物三维建模近年来一直是地理信息系统（Geographic Information System, GIS）、测绘遥感和计算机视觉等领域的研究热点。基于点云的建筑物三维重建不仅可实现自动化建模，而且在建模精度、语义信息提取、几何语义一体化表达等方面都有优势，因而成为了建筑物高精度重建的一种主流技术[1-2]。近年来国内外大量学者对其进行了研究，杨必胜等[3⇓-5]对点云数据处理和几何模型重建进展进行了总结。
当前基于点云的建筑物三维重建主要分为数据驱动(data-driven)和模型驱动(model-driven) 2类建模方法[6]。基于数据驱动的三维重建是目前建筑物点云三维建模的一种主流方法，主要有区域增长、随机抽样一致性(Random Sample Consensus, RANSAC)[7]、以PointNet[8]为代表的深度学习网络模型等几类典型算法。赵传等[9]提出一种区域增长与RANSAC相结合的机载LiDAR点云分割方法，用以分割不同复杂程度的建筑物屋顶面，能有效识别出面积较小的屋顶面。Li等[10]利用区域增长算法先将点云分割为不规则三角网区域组成的屋顶面片集合，再将其投影到二维格网地图，并引入能量优化方法对格网区域进行标记，最终利用屋顶轮廓线向量优化引导方法实现屋顶面片精细分割及屋顶三维重建。然而，由于点云数据量大、分布不均匀、结构复杂等特性，从大规模三维点云中精确分割出属于同一个平面或曲面的点仍具有挑战性[11]。除了需要从杂乱无序的点云中高质量提取平面对象之外，建筑物三维模型的边界正则化处理也是一大挑战。此外，这类算法分割得到的往往是一系列平面或多边形集合，只能构建几何格网（Mesh）模型，缺乏建筑物各个表面的语义信息，限制了模型的应用范围。
模型驱动方法采用自上而下的重建策略，首先建立基元（Primitive）库，将建筑物几何要素和拓扑关系隐含在基元的参数化描述中，再采用最小二乘优化等方法将点云与基元库中预定义的某类基元进行模型匹配，最终重建出建筑物三维模型。鉴于建筑物屋顶风格的多样性和复杂性，屋顶结构的三维重建是整个建筑物三维建模的核心和难点，目前大量研究集中在建筑物屋顶三维重建。Zheng和Weng[12]在建筑物二维轮廓数据辅助下，将不规则建筑物轮廓划分为互不重叠的系列四边形，再通过决策树匹配确定每个四边形对应的简单屋顶基元类型，然后使用静力矩方程从LiDAR点云中求解每个基元的长、宽和主方向等参数，最后将重建的相邻简单屋顶基元组合形成复杂的三维屋顶模型。Xiong等[13]通过对建筑屋顶面片无向图的分解和认知，以节点、边和最小循环数为基础构建出屋顶拓扑图(Roof Topology Graphs, RTG)，利用RTG来识别屋顶并选取相应的基元，并采用最小二乘法来完成点云与屋顶平面基元的拟合，最后对优化后的屋顶平面进行拓扑运算来获得完整的建筑物模型。除了采用数据驱动、RTG等方法识别屋顶结构外，部分学者尝试将深度学习模型引入屋顶点云类型识别。Huang 等[14]设计了11种屋顶基元类型，采用参数化表达的二维拓扑图方式进行定义，然后对点云进行栅格化，在影像空间通过规则进行基元合并与融合，最终实现复杂屋顶三维重建。Li等[15]和Zhang等[16]定义了几种常见的屋顶结构三维基元，使用PointNet++[17]网络模型对屋顶点云数据进行高精度分类，在此基础上提出了一种模型驱动的参数化建模方法，应用于欧美地区不同风格居民建筑的大规模三维重建。模型驱动方法融入了先验假设，能够保证重建模型拓扑关系的正确性，也有利于语义信息提取，但是还存在点云数据与基元类型匹配错误、模型拟合精度不高等问题。
当前基于点云的建筑三维重建算法大多是从几何维度关注模型精度和抽象程度，而对模型缺少必要的语义描述，其建模结果通常用Mesh模型或参数化模型等几何体表达，仅限于可视化展示，不利于建筑模型及局部构件的查询、编辑与分析，限制了其深层次应用。从大规模非结构化点云中自动重建几何和拓扑均正确的建筑物三维模型仍然是点云重建 的重点和挑战[18-19]。随着城市地理标记语言(City Geography Markup Language, CityGML)国际开放标准在GIS领域的推广和应用[20]，建筑物几何语义一体化建模成为了近年来一个新的研究方向。CityGML定义的建筑物三维模型强调几何、拓扑和语义表达的一致性，包含屋顶、墙面、地面、门窗、室内家具等丰富的语义信息，可弥补传统几何模型在数据共享和互操作等方面的不足，因而在城市规划、三维地籍、建筑设施管理等领域有广泛应用[21]。
建筑物几何语义一体化三维建模是对纯几何建模的一种有益补充，可大大增强三维模型数据的价值，促进城市建筑的精细化管理和智能化应用。针对当前点云三维重建方法难以适应复杂屋顶结构、几何精度不高，且模型语义信息严重缺乏等问题，本文提出了一种数据驱动与模型驱动相结合的复杂屋顶点云几何语义一体化三维重建算法。该算法以设计一组具有分段光滑表面的建筑3D实体基元为基础，先采用数据驱动方法进行复杂屋顶基元类型识别，然后利用点到3D基元格网表面的最短平均距离为目标函数，通过支持约束条件的最小二乘优化算法动态估计候选基元的正确参数，最后将参数化几何模型转换为几何语义一体化表达的CityGML模型。本算法能够让点云中的所有点对象在统一框架下同时参与3D基元拟合，从而能够获取全局最优参数，使得重建的三维模型具有较高的几何精度。同时，结合模型驱动方法和CityGML的优点，本算法能够降低对点云噪声和密度的敏感度，保证重建模型语义和拓扑关系的正确性。

2 建筑物3D基元库设计

基元是模型驱动方法的基础，为了支持各类几何形状的屋顶三维重建，本文设计了一套类型丰富、覆盖面广、代表性强的建筑物基元库，并对每个基元进行合理的参数化表达。

2.1 3D基元库设计

传统的基元库往往采用平面、长方体、圆柱体和圆锥体等相对简单的几何模型，难以代表风格多样的屋顶结构。本文根据建筑屋顶造型以及建筑形状的不同，设计了一套具有一定代表性的3D基元库，包含简单屋顶和复杂屋顶2类基元。其中，简单屋顶包括平屋顶、双坡顶、四坡顶、四角攒尖顶等类型；复杂屋顶主要以双坡顶和四坡顶的组合为主，有L型和T型2类组合。每类基元都是对建筑物屋顶、立面和地面的整体表达，构成封闭的3D实体对象，如图1所示。
图1 建筑物3D基元库

Fig. 1 3D primitive library of buildings

2.2 3D基元参数化表达

针对基元库中的每个3D基元，统一用形状(shape)和位置(position) 2类参数来表示。形状参数 θ s通常包括长度(l)、宽度(w)、高度(h)等。位置参数 θ p则描述基元在局部坐标系下的位置和方向，通常由一个绕Z轴旋转的方位角( κ)和3个平移参数( t x , t y , t z)来表示。基元库中所有基元的位置参数数量相同，但是其形状参数的种类与数量则因几何形状的差异而各不相同。以图2所示的双坡L型3D基元为例，可以用6个形状参数来表示： X轴方向建筑长度( l x)和宽度( w x)、 Y轴方向建筑长度( l y)和宽度( w y)、Z轴方向屋顶高度( h r)和墙体高度( h w)，即 θ s = { l x , l y , w x , w y , h r , h w }。假设局部坐标系原点O为基元2条屋脊线的交点在底面投影，根据基元的形状参数可确定其所有顶点在局部坐标系下的初始坐标，如 V 0 = l x - w x 2 , w y 2 , 0。
图2 双坡L型建筑3D基元参数

Fig. 2 Parameterization of a L-type building primitive with two gables

当基元根据其位置参数进行整体平移和旋转时，变换后的顶点坐标( V ')可以通过基元的位置参数 θ p = { t x , t y , t z , κ }来计算，其公式如下：
V ' = R k V + [ t x , t y , t z ] T
式中：V表示局部坐标系下基元某顶点三维坐标，对应图2中 V 0 ~ V 14 中 任 意 一 个； t x 、 t y 、 t z分别表示基元在局部坐标系下沿X、Y和Z轴方向的平移量； k表示方位角； R k为沿Z轴逆时针旋转 k角度的旋转矩阵，表示为：
R k = c o s ( κ ) - s i n ( κ ) 0 s i n ( κ ) c o s ( κ ) 0 0 0 1
从上述定义可以看出，利用基元的形状和位置参数就可以确定3D基元所有顶点在局部坐标系下的唯一坐标。后续通过点云与3D基元的拟合来求解基元的正确参数，就能实现各类屋顶的三维重建。建筑物顶点、边、面片以及它们之间的拓扑关系都在3D基元的参数化表达中进行了明确定义，包括建筑物各表面之间的内在几何约束，如对称、平行或垂直等关系也隐含在3D基元中。此外，建筑物各个表面的语义信息也隐含在这种参数化表达的3D基元中。以图2为例， v 6 , v 7 , v 14 , v 12和 v 7 , v 8 , v 13 , v 14构成的2个平面属于屋顶面(语义：RoofSurface)， v 0 , v 1 , v 7 , v 6为其中一个墙面(语义：WallSurface)， v 0 , v 1 , v 2 , v 3 , v 4 , v 5则形成地面(语义：GroundSurface)。因此，设计的这种3D实体基元有利于快速构建几何、语义和拓扑结构均正确的建筑物三维模型。

3 建筑物点云基元类型识别

定义了可参数化表达的3D基元库之后，采用模型驱动方法对点云数据进行三维重建的首要任务是识别点云对应的3D基元类型。本文首先从点云中提取出屋顶面片，然后建立屋顶面片之间的拓扑关系，最后通过对比点云与3D基元的屋顶拓扑图来判别点云所属分类。

3.1 点云基元类型识别总体流程

本文设计的点云与建筑物3D基元类型匹配流程如图3所示。首先，对基元库中的所有基元使用邻接矩阵描述其屋顶拓扑关系。其次，对建筑点云数据，使用RANSAC算法提取屋顶面片，经过面片优化、几何中心求解、屋脊线提取等过程，建立描述点云屋顶拓扑关系的邻接矩阵。最后，将点云的邻接矩阵与各基元的邻接矩阵逐一比较，找到最匹配的候选基元。下面重点介绍屋顶面片提取和屋顶拓扑关系构建。
图3 点云基元类型识别流程

Fig. 3 Procedure of recognizing the roof structure from a point cloud

3.2 基于RANSAC分割的屋顶面片提取

使用RANSAC算法对屋顶点云进行分割。由于RANSAC分割结果依赖阈值和迭代次数等参数设置，容易出现过分割或欠分割等情况，因而还需要对分割得到的面片进行后处理。针对过分割结果，参考RANSAC划分内点的思路，合并应该属于同一屋顶面的面片。具体步骤如下：
（1）设置面片法向量夹角阈值 θ和面片中心点之间的空间欧式距离阈值 d。对于未被合并的面片集 { P 1 , P 2 , ⋯ , P n }中每一个面片 P i，计算其法向量与其他面片 P j ( j ≠ i )法向量的夹角 ϑ i j以及面片中心点之间距离 d i j。如果 ϑ i j < θ，或者 ϑ i j < 2 θ 且 d i j < d，标记 { P i , P j }为合并候选面片对。否则，认为面片 P i不需要合并，继续判断下一面片 P i + 1。
（2）对于提取出的候选面片对 { P i , P j }，分别计算面片 P i(或 P j)所有点到另一面片 P j(或 P i)的垂直距离，并统计距离小于阈值 d v的点在该面片所有点中的占比 r i和 r j，如果占比均超过50%或占比之和大于某一阈值，则认为此面片对属于同一个面片，将其合并，标记面片 P j为被合并的面片，并用合并后的所有点重新拟合平面，将得到的平面参数更新至面片 P i，参与下一面片的合并判断。重复以上两步，直至所有面片合并完毕。
针对没有考虑点云空间位置关系等原因造成的面片内孤立点的情况，则基于该点到各面片的垂直距离和其邻域范围内各面片的点数进行判断。点到面片 P j的垂直距离越小，邻域点中属于面片 P j的点数越多，则该点属于面片 P j的概率越大。对屋顶点云中属于面片 P i的任一点 p k i，计算其到各面片的垂直距离 { d k 1 , d k 2 , ⋯ , d k n }，并统计其邻域内屋顶面片相应的点数 { s k 1 , s k 2 , ⋯ , s k n }，记录其距离最小且邻域点数量最多的面片 P j，若 j ≠ i，则更新该点所属的面片为 P j。对所有点进行优化处理之后，为避免屋檐、天窗和烟囱等屋顶附属结构的影响，可进一步根据平面系数剔除近似水平面片，最终获取能正确代表屋顶结构的面片。由于RANSAC分割仅用于识别主要屋顶面片拓扑结构，点云高程误差和部分噪声点造成的分割不精确并不会影响后续拓扑关系建立。

3.3 屋顶拓扑关系建立

正确提取出屋顶面片之后，还需要利用建筑结构先验知识来寻找点云与3D基元的共性表示，以此判断点云与基元之间的匹配关系。本文参考文献[13]提出的屋顶拓扑图(RTG)来实现点云所属基元类型的识别。
为了准确唯一表示屋顶结构，根据屋顶面片之间的交线（屋脊）属性将屋顶拓扑图进一步划分为带属性的RTG。首先根据屋脊的方向划分为水平和倾斜屋脊，然后根据相邻面片之间形成的夹角，划分为向上的凸角和向下的凹角。因此，该属性可以划分为：仅有一个水平面的水平脊，2个屋顶面形成凸角和凹角的凸水平脊和凹水平脊，以及由具有凸角和凹角的倾斜屋顶面产生的凸倾斜脊和凹倾斜脊，如图4所示，边的颜色代表该属性。依据该规则建立基元库中各3D基元对应的带属性RTG，见图5。这种带属性的RTG可以唯一确定屋顶类型，通过逐一比较点云和所有3D基元的RTG来识别候选基元类型。
图4 使用带属性的边表示2个相邻屋顶面片的拓扑关系

Fig. 4 Adjacent relationships between two roof planes represented by edges with attributes

图5 建筑物3D基元屋顶拓扑图

Fig. 5 Roof topology graphs of the predefined 3D building primitives

4 基于3D基元拟合的复杂屋顶几何语义一体化建模

基于3D基元拟合的复杂屋顶几何语义一体化重建流程如图6所示。其中，点云与候选3D基元的拟合是重建核心，涉及拟合目标函数定义、目标函数约束条件设置和参数初始化等关键步骤。在此基础上通过最小二乘优化算法估计出基元最佳参数，即可形成参数化表达的几何模型。最后，再提取参数化模型各个表面的语义信息，按照CityGML标准转换为几何语义一体化表达的建筑物模型。
图6 复杂屋顶几何语义一体化CityGML模型重建流程

Fig. 6 Workflow of reconstructing geometric-semantic coherent CityGML building model with complex roof

4.1 点云与3D基元拟合目标函数

利用候选3D基元对点云进行建模，本质上是在特定的目标函数下估计3D基元参数的最优值。本文利用点云中所有点对象到3D基元表面的最小距离平方的平均值作为目标函数，寻找点云与3D基元之间的最佳拟合状态。目标函数定义如下：
J θ = m i n 1 N ∑ k = 1 N d i s t p k , M 2
式中： p k代表点云中的第k个点； N为点云总数量； M代表3D基元的Mesh模型； θ为3D基元的参数集。以双坡L型基元为例， θ就是6个形状参数与4个位置参数的集合，即：
θ = l x , l y , w x , w y , h r , h w , t x , t y , t z , κ
从目标函数定义可以看出，求解该目标函数的关键是如何计算点到3D基元表面的最短距离PMD (Point-Mesh-Distance)，其具体计算过程可参考文献[16]。使用PMD作为目标函数，适用于各种3D基元的参数寻优。因为无论基元表面几何形状是否规则，都可以转换为三角网(Triangular Mesh)来表示，这样便只需要通过求解点到基元所有三角形的最短距离，即可求得点到该基元的最短距离。另一方面，利用PMD作为目标函数可以使所有点同时参与优化，与RANSAC这类平面拟合算法每次只有局部点参与优化相比，本文所建立的目标函数适用分段光滑曲面，每次迭代能够让3D基元的所有平面同时参与拟合，从而实现PMD度量下的参数全局寻优。

4.2 目标函数约束条件

为了充分利用3D基元隐含的先验知识，并获得更加稳定和精确的优化结果，目标函数可引入额外的约束条件。
以图2所示的双坡L型基元为例，对于二维边界的约束，首先由Alpha Shapes[22]或凸包算法检测点云的二维投影边界，并计算投影后的凸多边形面积(S)，然后用该面积约束基元的长宽形状参数，即：
l x ∙ w y + l y - w y ∙ w x ≤ S
类似地，针对点云的垂直边界 ( h m i n , h m a x )，可通过式(6)进行约束。
h r + h w ≤ h m a x - h m i n
考虑到建筑物宽度通常小于其长度，且屋顶高度通常低于墙体高度，还可增加额外约束，定义如下：
w x ≤ l x , w y ≤ l x w x ≤ l y , w y ≤ l y h r ≤ h w
综上，为了对双坡L型3D基元进行更加准确的参数估计，可针对目标函数（式（3））定义如下约束条件。
S - l x ∙ w y + l y - w y ∙ w x ≥ 0 h m a x - h m i n - h r - h w ≥ 0 l x - w x ≥ 0 l x - w y ≥ 0 l y - w x ≥ 0 l y - w y ≥ 0 h w - h r ≥ 0
不同3D基元的几何结构不同，其目标函数的约束条件也会有所不同，以图7所示的四坡顶基元为例，其屋脊线长度一般会小于其建筑长度，可为L型四坡顶和T型四坡顶分别增加约束，如式（9）和式（10）所示。
图7 四坡顶3D基元形状参数示例

Fig. 7 The shape parameters of building primitives with hip roofs

r x ≤ l x - 1 2 w x r y ≤ l y - 1 2 w y
r x ≤ 1 2 l x r y ≤ l y - 1 2 w y

4.3 目标函数参数初始化

针对上述带约束条件的目标函数，参数初值的合理性会影响目标函数优化过程，一定程度上影响参数估计的精确性和算法效率，尤其是针对复杂屋顶点云。因此，在使用优化算法求解之前，还需要对3D基元参数进行合理的初始化。
为了估计出相对准确的基元形状参数初始值，需要以点云的分布方向作为局部坐标轴方向，即让局部坐标系的X轴与建筑主方向基本平行，以两屋脊线交点为原点，在此基础上根据点云的分布范围估计基元的形状参数初值。
考虑到屋脊点一般位于建筑物顶部，首先从点云中选择一定数量的Z坐标较大点作为屋脊候选点。随后使用K-Means算法将候选点分为两类，剔除可能是天窗、烟囱的噪声点之后，对候选点使用RANSAC算法进行直线拟合。虽然拟合的直线不能精确代表屋脊线的准确位置，但2条直线的夹角可以作为方位角 κ的初值。
双坡与四坡屋顶基元形状参数初值的设置见表1，二者共有的6个形状参数 θ s = { l x , l y , w x , w y , h r , h w }中， X轴方向建筑长度( l x)和Y轴方向建筑长度( l y)，取该方向坐标的最大值和最小值之差作为其初值 l x 0和 l y 0。对于 X轴方向建筑宽度( w x)和Y轴方向建筑宽度( w y)，首先在其垂直方向上取一定数量的最大值作为边界候选点，然后计算边界候选点在该坐标轴方向坐标的最大值和最小值之差作为其初始值 w x 0和 w y 0。针对屋顶垂直高度( h r)和墙体垂直高度( h w) 2个参数，取 Z轴方向的最大值和最小值之差，按3:7的比例分配为初始值 h r 0和 h w 0。在此基础上，四坡基元还有 X轴方向屋脊长度( r x)和Y轴方向屋脊长度( r y)的形状参数，取该方向屋脊线上候选点的坐标最大值作为其初始值 r x 0和 r y 0。
表1 建筑物3D基元形状参数初值

Tab.1 Initial values of shape parameters for some building primitives

形状参数 数学表达 描述
l x 0 l x 0 = m a x x - m i n ( x ) x坐标的最大值与最小值之差
l y 0 l y 0 = m a x y - m i n ( y ) y坐标的最大值与最小值之差
w x 0 w x 0 = m a x y m a x x - m i n y m a x ( x ) 取一定数量的 y坐标最大的点，求其 x坐标最大值与最小值之差
w y 0 w y 0 = m a x x m a x y - m i n x m a x ( y ) 取一定数量的 x坐标最大的点，求其 y坐标最大值与最小值之差
h r 0 h r 0 = 0.3 × ( m a x z - m i n ( z ) ) z坐标最大值与最小值之差×0.3
h w 0 h w 0 = 0.7 × ( m a x z - m i n ( z ) ) z坐标最大值与最小值之差×0.7
r x 0 r x 0 = m a x r i d g e _ x x X轴方向上的屋脊线候选点的最大 x坐标
r y 0 r y 0 = m a x r i d g e _ y y Y轴方向上的屋脊线候选点的最大 y坐标

4.4 基于最小二乘优化的参数估计

从式（3）定义的目标函数可知，基元参数估计相当于一个非线性最小二乘问题求解。考虑到4.2节定义的相关约束，本文采用最小二乘优化算法中的序列二次规划(Sequential Quadratic Programming, SQP)算法[23]求解约束最优化问题。SQP算法是目前公认的求解约束非线性问题的最有效方法之一，具有全局收敛性好、边界搜索能力强和计算效率高等优点。除设置约束条件外，还可以根据点云分布情况定义相关参数的取值范围。对目标函数的各项参数进行合理的初始化之后，SQP算法会根据这些约束和边界限制对基元的所有参数通过多次迭代不断更新。当目标函数获得最小值时，点云与3D基元达到最佳拟合状态，此时的参数即为待求解的全局最优参数。

4.5 几何语义一体化模型构建

利用SQP算法求解候选3D基元的正确参数后，即可形成参数化表达的几何模型。为了生成几何语义一体化的建筑物模型，还需要按照CityGML标准进行各类语义对象的存储和表达。CityGML定义的建筑物模型包含多细节层次(Levels of Detail, LoD)表达，最新的CityGML 3.0标准将建筑物模型分为4级LoD，其LoD2模型对屋顶样式进行了显式表达[24]。因此，本文主要将参数化表达的几何模型转换为CityGML LoD2模型。CityGML使用边界表示法(Boundary Representation, B-Rep)描述各个对象的面、边缘、顶点及它们之间的关系，进而表达三维对象的几何形状。CityGML LoD2建筑物模型需要显式表达的语义对象主要有：RoofSurface、WallSurface、GroundSurface等，而预定义的3D基元已经隐含了各个平面的语义信息，只需按照逆时针顺序将几何模型每个平面的顶点坐标串 (posList) 存储到相应语义平面的多边形对象(Polygon)中即可。由于语义对象及其坐标可以从参数化3D基元直接提取，因而CityGML LoD2模型能利用程序快速自动化构建。

5 实验与分析

5.1 实验数据

使用2种不同类型的点云公开数据集进行实验，第一个是国际摄影测量与遥感学会（International Society for Photogrammetry and Remote Sensing， ISPRS）提供的来自德国Vaihingen地区的 LiDAR点云标准数据集[25]，整个数据集由10条机载激光雷达条带组成，通过Leica ALS50系统获得，点云平均密度为4 points/m2。第二个是senseFly提供的来自瑞士Merlishachen地区的密集摄影测量点云数据集[26]，使用Pix4D软件从eBee无人机影像生成，点云平均密度为64.67 points/m2。首先使用CloudCompare开源软件从2个数据集中人工选择并分割出几个代表不同屋顶类型的单体建筑物点云，然后对单体建筑物点云进行噪声滤波处理。

5.2 实验结果及分析

5.2.1 复杂屋顶基元类型识别

采用点云和3D基元的RTG邻接矩阵比较法来确定点云对应的屋顶基元类型，其结果如表2所示。考虑到sensefly_2、sensefly_3等少量建筑存在烟囱，先滤除高程值最大的100个点之后，再从中提取屋脊候选点，这样有利于屋脊线提取和RTG构建的准确性。从点云识别出的屋顶类型与参考影像对比可以看出，RANSAC分割算法具有一定的鲁棒性，提取的面片经优化处理之后能准确代表屋顶结构，有利于点云数据与3D基元的正确匹配。
表2 点云屋顶类型识别结果

Tab. 2 Roof type recognition results of different point clouds

ID 点云分割结果 邻接矩阵 建筑类型 参考影像
sensefly_1 0 1 3 0 1 0 0 2 3 0 0 1 0 2 1 0 双坡L型
sensefly_2 0 2 1 0 0 2 2 0 0 2 1 0 1 0 0 0 3 2 0 2 0 0 2 0 0 1 3 2 0 0 2 0 2 0 0 0 四坡L型
sensefly_3 0 1 2 0 2 0 1 0 0 0 2 3 2 0 0 2 0 1 0 0 2 0 0 2 2 2 0 0 0 0 0 3 1 2 0 0 四坡L型
ISPRS_1 0 1 3 0 1 0 0 2 3 0 0 1 0 2 1 0 双坡L型
ISPRS_2 0 1 0 0 2 0 3 1 0 0 0 2 0 3 0 0 0 2 0 2 1 0 0 2 0 0 0 2 2 2 0 0 0 0 0 0 0 2 0 0 0 2 3 3 1 2 0 2 0 四坡T型

5.2.2 3D基元参数优化

以双坡L型建筑3D基元为例，图8展示了其参数初始化流程，输入点云数据来自senseFly数据集。首先从去噪后的点云中选取100个Z值较大的点作为屋脊候选点，用K-means算法将其分为 2类。接下来使用RANSAC直线拟合算法从两类候选点中分别提取2条屋脊线，计算屋脊线交点和夹角，据此对点云做局部坐标变换以获取其形状参数的初始值。图8(a)为选取100个点作为屋脊候选点的提取结果，图8(b)是使用K-means算法得到的聚类结果，图8(c)展示了RANSAC拟合出的 2条屋脊线，图8(d)为根据屋脊线夹角和交点进行局部坐标变换后的点云，以及用其确定的基元形状参数和位置参数初值构建的3D基元。对于其他类型的建筑，同样根据识别的点云建筑类型从基元库中提取候选基元，执行相同的基元初始化过程。
图8 3D基元参数初始化实验流程

Fig. 8 The work flow of parameter initialization for a gable-L type building primitive

从表3可以看出，尽管有部分3D基元的初始状态与点云有较大差异，但是通过多次迭代优化后，基元与点云之间的差异越来越小，当点云到3D基元表面的距离最小时，表明点云与基元达到了最佳拟合状态。利用该状态下获取的参数构建的建筑物三维模型就是点云几何重建结果。图10中虽然有几个建筑物存在烟囱或天窗，这些构成烟囱或天窗的点由于所占比例小，且几何结构面积不大，没有影响建筑物主体结构的拟合，这也说明本文使用的模型驱动方法对噪声和局部点云缺失不敏感。
表3 点云与3D基元拟合结果

Tab. 3 3D primitive fitting result of point clouds with different roof types

ID 点云数据 点云与基元拟合结果 参考影像
sensefly_1
sensefly_2
sensefly_3
ISPRS_1
ISPRS_2
点云对象PMD值/m
采用PMD平均值作为量化指标来评估重建模型的几何精度，上述5个建筑点云的精度统计结果见表4。可以看出所有建筑基本都实现了较为精确的重建，除ISPRS_1外，大部分建筑的PMD值小于0.20 m。仔细观察ISPRS_1点云和参考影像，发现其PMD值较大主要是因为该建筑物天窗较多，存在大量噪声点。定性和定量分析实验结果表明，采用建筑物3D基元对点云进行整体拟合的三维重建方法是有效且鲁棒的，不仅能避免噪声和点云空洞等影响，而且重建模型的几何结构和拓扑关系均正确。
表4 建筑点云三维重建结果的几何精度统计

Tab. 4 Geometric accuracy statistics of building reconstruction from several point clouds

建筑ID 点数/个 屋顶类型 平均PMD/m
sensefly_1 5 463 双坡L型 0.086 0
sensefly_2 5 879 四坡L型 0.102 0
sensefly_3 12 082 四坡L型 0.102 7
ISPRS_1 3 864 双坡L型 0.386 6
ISPRS_2 1 665 四坡T型 0.164 5

5.2.3 几何语义一体化建模

将参数优化后获得的建筑物几何模型按照基元隐含的语义信息转换为CityGML LoD2模型，其结果如表5所示。
表5 CityGML LoD2模型构建结果

Tab. 5 Generated CityGML LoD2 building models from parameterized models

ID 参数化几何模型 CityGML LoD2模型
sensefly_1
sensefly_2
sensefly_3
ISPRS_1
ISPRS_2

注：红色表示屋顶面，灰色表示墙面和地面。

6 讨论

6.1 约束条件对重建精度的影响分析

目标函数约束条件会影响最小二乘优化算法的正确收敛，进而影响模型重建精度。以sensefly_2 四坡L型建筑点云为例，无约束条件的优化结果见图9(b)。可以看出，其中决定四坡屋顶的一个屋脊点不仅没有得到正确估计，甚至超出了建筑的边界范围，从而导致错误的基元拟合参数。若增加四坡顶的屋脊线长度不大于其建筑长度这一约束条件，见式（11）。
图9 有无约束条件拟合结果对比

Fig. 9 Comparison of primitive fitting results without and with constraints

r x ≤ l x - 1 2 w x r y ≤ l y - 1 2 w y
获得的拟合结果如图9(c)所示。可以看到屋顶各个坡面均取得较好拟合，平均PMD值从0.187 3 m减小到0.102 0 m，几何精度有较大提高。

6.2 不同建模方法对比

为了进一步验证本方法的有效性，将其与经典的Delaunay算法、RANSAC平面拟合算法[7]，以及Li[15]提出的基于模型驱动的屋顶重建算法进行对比分析。其中，Delaunay方法计算点云在最佳平面的2.5D三角剖分，并构建Mesh表面模型；RANSAC算法主要进行平面基元提取，每个平面包含的最小点数设置为20；Li[15]针对三种类型的简单屋顶基元，该方法设计的目标函数为点到基元平面的垂直距离(Point-Surface-Distance, PSD)和二维投影边界的交集与并集之比(Intersection over Union, IoU)的线性组合。
点云中所有点的PMD平均值是衡量模型重建几何精度的一种有效指标，采用该指标对表3展示的5个建筑点云的4种不同方法建模结果进行精度对比分析，其结果见表6。
表6 不同建模方法的几何精度比较

Tab. 6 Comparison of geometric accuracy for different methods

建筑ID 平均PMD/m
Delaunay RANSAC Li[15]方法 本文方法
sensefly_1 0.989 5 0.097 2 0.420 7 0.086 0
sensefly_2 0.511 1 0.138 2 1.491 6 0.102 0
sensefly_3 0.491 8 0.106 9 0.297 9 0.102 7
ISPRS_1 0.593 8 0.596 3 1.025 0 0.386 6
ISPRS_2 0.275 0 0.827 3 1.618 6 0.164 5
由表6可看出，本文方法得到的PMD均小于其他3种方法，说明本文方法的几何精度最高。RANSAC算法虽然得到的PMD值也相对较小，但是其只能提取离散平面集合，不能构建封闭的建筑物模型，且存在大量拓扑错误。Delaunay方法也是通过最佳平面拟合方式构建点云三角网，由于墙面和地名点云密度较低，因此也不能生成完整的建筑物模型，其精度相对较低。Li[15]使用PSD和IoU的线性组合作为度量指标，优化时要同时兼顾2个指标，并且还要受指标权重影响。包含地面点、建筑边界噪声和其他地物噪声点的建筑点云数据一定程度上会影响IoU计算。如果目标函数包含IoU指标，通常会通过扩大模型的边界来提高IoU，导致点到模型表面的平均距离增大，从而使得优化后的参数模型难以与点云实现最佳拟合。本文仅采用PMD一个度量指标作为目标函数，不仅比PSD具有更加严格的限制，而且优化时只考虑全局PMD最小，不受多个指标权重参数影响，其灵活性和拟合效果更好。
本文算法使用Python语言编程实现，测试计算机硬件配置为Intel(R) Core(TM) i9-10900K CPU @3.70 GHz，32 GB内存，4种方法的运行效率对比如表7所示。
表7 不同建模方法的效率对比

Tab. 7 Comparison of efficiency for different methods

建筑ID 耗时/s
Delaunay RANSAC Li[15]方法 本文方法
sensefly_1 0.70 1.49 517.79 84.27
sensefly_2 1.02 1.35 1 274.70 162.99
sensefly_3 2.98 1.68 2 158.06 302.87
ISPRS_1 0.96 1.15 219.10 44.76
ISPRS_2 0.83 0.98 99.98 47.47
从表7可以看出，Delaunay和RANSAC方法效率最高，主要是因为这2种方法只考虑平面拟合，未考虑平面之间的约束关系和拓扑错误修正，得到的结果存在大量几何和拓扑错误。本算法和Li[15]的耗时主要体现在基元参数寻优环节，但本文算法效率要明显高于Li[15]，因为本文的优化只考虑PMD一个度量指标计算，而Li[15]需要同时计算PSD和IoU，并且本算法增加了各种约束条件，可以更快收敛。本方法针对sensefly_3点云的耗时最长，因为该数据包含的点最多(12 082)，且基元是四坡L型，每次参与PMD距离计算的表面也最多，故参数优化需要更长时间。
综合几何精度和运行效率对比来看，本算法原理简单清晰，普适性强，准确性和效率均较高，能够满足各类复杂屋顶三维重建需求。但是，本方法也存在不确定性，点云分割精度、点云与基元RTG匹配准确率、目标函数约束条件和参数初值设置合理性等因素一定程度上会影响重建精度。

7 结论

从大规模点云重建几何精确、拓扑正确和语义丰富的建筑物三维模型是当前点云三维建模的难点和重点，尤其是针对复杂屋顶结构。本文提出了一种基于3D基元拟合的复杂屋顶点云三维自动重建方法，采用具有不同屋顶形状的LiDAR点云和摄影测量点云进行三维重建实验，定性分析和定量对比结果证明了该方法的有效性和准确性。本方法创新之处体现在：
（1）定义了一套可参数化表达的建筑物3D封闭实体基元库，包含简单屋顶和复合屋顶等各种类型，可灵活扩展。
（2）设计的点云与3D基元拟合目标函数普适性强，而且无需手动调节参数，模型参数求解精度和算法自动化程度较高，可用于简单屋顶和复杂屋顶等各类几何结构不规则的建筑物三维重建。
（3）基于3D基元拟合的重建方法可以克服数据驱动等基于简单平面或连续光滑表面拟合方法的局限性，不受面片分割精度影响，对噪声、局部点云缺失具有鲁棒性。
（4） 3D基元隐含的对称性、平行和正交关系、面片语义信息等先验知识保证了重建结果具有正确的拓扑关系和丰富的语义特征。
本方法综合了数据驱动和模型驱动方法的优点，且重建结果为几何语义一体化表达的CityGML模型，增强了建筑物三维模型的互操作、空间查询与分析能力，有利于促进建筑设施的精细化管理和智能化应用。尽管本文设计的重建算法能够实现几类组合屋顶的建模，但该方法还存在一定的局限性。首先，基元库中预定义的3D基元有限，难以覆盖所有复杂屋顶类型。其次，如何自动准确地从点云中识别复杂多变的建筑类型也是一项挑战。未来将进一步扩展基元库，引入更多类型的复杂屋顶基元。由于建筑结构复杂化导致基元参数增加的情况，还需要研究如何确定参数初值和约束条件，并提高大规模点云数据参数估计效率。

References

Publishing order | Descend order by publishing year | Descend order by cited within
[1]
单杰, 李志鑫, 张文元. 大规模三维城市建模进展[J]. 测绘学报, 2019, 48(12):1523-1541.

DOI

[Shan J, Li Z X, Zhang W Y. Recent progress in large-scale 3D city modeling[J]. Acta Geodaetica et Cartographica Sinica, 2019, 48(12):1523-1541.] DOI:10.11947/j.AGCS.2019.20190471

DOI

[2]
文学东, 陈为民, 谢洪, 等. 一种融合多源特征的建筑物三维模型重建方法[J]. 武汉大学学报·信息科学版, 2019, 44(5):731-736,764.

[Wen X D, Chen W M, Xie H, et al. A method for building model reconstruction based on multi-source feature fusion[J]. Geomatics and Information Science of Wuhan University, 2019, 44(5):731-736,764.] DOI:10.13203/j.whugis20180320

DOI

[3]
杨必胜, 梁福逊, 黄荣刚. 三维激光扫描点云数据处理研究进展、挑战与趋势[J]. 测绘学报, 2017, 46(10):1509-1516.

DOI

[Yang B S, Liang F X, Huang R G. Progress, challenges and perspectives of 3D LiDAR point cloud processing[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(10):1509-1516.] DOI:10.11947/j.AGCS.2017.20170351

DOI

[4]
杜建丽, 陈动, 张振鑫, 等. 建筑点云几何模型重建方法研究进展[J]. 遥感学报, 2019, 23(3):374-391.

[Du J L, Chen D, Zhang Z X, et al. Research progress of building reconstruction via airborne point clouds[J]. Journal of Remote Sensing, 2019, 23(3):374-391.] DOI:10.11834/jrs.20188199

DOI

[5]
Xu Y S, Stilla U. Toward building and civil infrastructure reconstruction from point clouds: A review on data and key techniques[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2021, 14:2857-2885. DOI:10.1109/JSTARS.2021.3060568

DOI

[6]
Maas H G, Vosselman G. Two algorithms for extracting building models from raw laser altimetry data[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 1999, 54(2/3):153-163. DOI:10.1016/S0924-2716(99)00004-0

DOI

[7]
Schnabel R, Wahl R, Klein R. Efficient RANSAC for point-cloud shape detection[J]. Computer Graphics Forum, 2007, 26(2):214-226. DOI:10.1111/j.1467-8659.2007.01016.x

DOI

[8]
Qi C R, Su H, Mo K C, et al. PointNet: deep learning on point sets for 3D classification and segmentation[C]// 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR). IEEE, 2017:77-85. DOI:10.1109/CVPR.2017.16

DOI

[9]
赵传, 郭海涛, 卢俊, 等. 结合区域增长与RANSAC的机载LiDAR点云屋顶面分割[J]. 测绘学报, 2021, 50(5):621-633.

DOI

[Zhao C, Guo H T, Lu J, et al. Roof segmentation from airborne LiDAR by combining region growing with random sample consensus[J]. Acta Geodaetica et Cartographica Sinica, 2021, 50(5):621-633.] DOI:10.11947/j.AGCS.2021.20200270

DOI

[10]
Li M L, Rottensteiner F, Heipke C. Modelling of buildings from aerial LiDAR point clouds using TINs and label maps[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2019, 154:127-138. DOI:10.1016/j.isprsjprs.2019.06.003

DOI

[11]
张良培, 张云, 陈震中, 等. 基于分裂合并的多模型拟合方法在点云分割中的应用[J]. 测绘学报, 2018, 47(6):833-843.

DOI

[Zhang L P, Zhang Y, Chen Z Z, et al. Splitting and merging based multi-model fitting for point cloud segmentation[J]. Acta Geodaetica et Cartographica Sinica, 2018, 47(6):833-843.] DOI:10.11947/j.AGCS.2018.20180131

DOI

[12]
Zheng Y F, Weng Q H. Model-driven reconstruction of 3-D buildings using LiDAR data[J]. IEEE Geoscience and Remote Sensing Letters, 2015, 12(7):1541-1545. DOI:10.1109/LGRS.2015.2412535

DOI

[13]
Xiong B, Jancosek M, Elberink S O, et al. Flexible building primitives for 3D building modeling[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2015, 101:275-290. DOI:10.1016/j.isprsjprs.2015.01.002

DOI

[14]
Huang H, Brenner C, Sester M. A generative statistical approach to automatic 3D building roof reconstruction from laser scanning data[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2013, 79:29-43. DOI:10.1016/j.isprsjprs.2013.02.004

DOI

[15]
Li Z X, Zhang W Y, Shan J. Holistic parametric reconstruction of building models from point clouds[J]. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2020, XLIII-B2-2020:689-695. DOI:10.5194/isprs-archives-xliii-b2-2020-689-2020

DOI

[16]
Zhang W Y, Li Z X, Shan J. Optimal model fitting for building reconstruction from point clouds[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2021, 14:9636-9650. DOI:10.1109/JSTARS.2021.3110429

DOI

[17]
Qi C R, Yi L, Su H, et al. PointNet++: Deep hierarchical feature learning on point sets in a metric space[C]// Proceedings of the 31st International Conference on Neural Information Processing Systems. New York: ACM, 2017:5105-5114. DOI:10.5555/3295222.3295263

DOI

[18]
Iman Zolanvari S M, Laefer D F. Slicing method for curved façade and window extraction from point clouds[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2016, 119:334-346. DOI:10.1016/j.isprsjprs.2016.06.011

DOI

[19]
Cao R J, Zhang Y J, Liu X Y, et al. 3D building roof reconstruction from airborne LiDAR point clouds: A framework based on a spatial database[J]. International Journal of Geographical Information Science, 2017, 31(7):1359-1380. DOI:10.1080/13658816.2017.1301456

DOI

[20]
Gröger G, Plümer L. CityGML - Interoperable semantic 3D city models[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2012, 71:12-33. DOI:10.1016/j.isprsjprs.2012.04.004

DOI

[21]
Biljecki F, Stoter J, Ledoux H, et al. Applications of 3D city models: State of the art review[J]. ISPRS International Journal of Geo-Information, 2015, 4(4):2842-2889. DOI:10.3390/ijgi4042842

DOI

[22]
Edelsbrunner H, Mücke E P. Three-dimensional alpha shapes[J]. ACM Transactions on Graphics, 1994, 13(1):43-72. DOI:10.1145/174462.156635

DOI

[23]
Boggs P T, Tolle J W. Sequential quadratic programming[J]. Acta Numerica, 1995, 4:1-51. DOI:10.1017/s0962492900002518

DOI

[24]
Kolbe T H, Kutzner T, Smyth C S, et al. OGC City Geography Markup Language (CityGML) part 1: conceptual model standard, Version 3.0.0[S]. Open Geospatial Consortium Inc., 2021.

[25]
ISPRS. ISPRS test project on urban classification and 3D building reconstruction[EB/OL]. [2022-3-16]. https://www2.isprs.org/commissions/comm3/wg4/detection-and-reconstruction.html

[26]
Sensnfly. Explore how senseFly drone solutions are employed around the globe[EB/OL]. [2022-3-8]. https://www.sensefly.com/education/datasets/?dataset=1419

Options
Outlines

/

〈 〉