CN102521884A

CN102521884A - 一种基于LiDAR数据与正射影像的3维屋顶重建方法

Info

Publication number: CN102521884A
Application number: CN2011104235604A
Authority: CN
Inventors: 程亮; 李满春; 刘永学; 童礼华; 张雯; 陈焱明; 蔡文婷; 李�真; 杨康; 潘航; 邹伟
Original assignee: Nanjing University
Current assignee: Nanjing University
Priority date: 2011-12-16
Filing date: 2011-12-16
Publication date: 2012-06-27
Anticipated expiration: 2031-12-16
Also published as: CN102521884B

Abstract

本发明属于影像处理方法领域，公开了一种基于LiDAR数据与正射影像的3维屋顶重建方法。本方法包括以下步骤：（1）基于三角形簇的LiDAR点屋顶面片分割；（2）基于LiDAR数据与正射影像的屋脊线提取；（3）3维屋顶模型重建。本方法面向3维屋顶模型精细重建的需求，集成LiDAR数据与高分辨率影像，综合利用LiDAR数据高程特性与影像高分辨率特性的互补优势，以“屋顶面片分割-屋脊线提取-3维屋顶模型重建”为主线，实现了基于三角形簇的LiDAR点屋顶面片分割算法、基于LiDAR数据与正射影像的屋脊线提取算法，形成了一种3维屋顶模型重建的新方法，实验证明本方法建模的自动化程度高，正确性和完整性较高、定位精度较高，符合实际应用需求。

Description

一种基于LiDAR数据与正射影像的3维屋顶重建方法

技术领域

本方法属于影像处理方法技术领域，特别是涉及一种基于LiDAR数据与正射影像相结合的高精度、自动化的3维屋顶重建方法。

背景技术

3维建筑模型是数字城市的重要组成部分，广泛应用于测绘遥感、虚拟现实、城市规划、科学研究等领域。通过人工数字化、3维手工建模构建3维建筑模型，效率低、成本高，很难满足大范围、实时建筑模型的需求。伴随着摄影测量技术、图像处理技术、计算机视觉技术、地理信息技术等学科的发展，各类建筑模型构建方法不断涌现，这些方法主要可以分为基于影像的自动化重建方法、基于LiDAR数据的自动化重建方法以及LiDAR与影像相结合的自动化重建方法。基于影像的自动化重建是目前最经济快捷的城市建筑物的建模方法，但这种方法存在一定问题：(1)影像数据包含大量的阴影以及噪声，对噪声的相关处理还不是很成熟；(2)像对匹配特征的选择比较困难；(3)利用摄影测量的方法从影像重建3维模型是一个综合性的复杂问题，完全的自动化提取不大现实。

近年来，一种新型的测量技术机载激光雷达发展迅速，并被广泛地应用于遥感测绘以及3维模型构建中，逐渐成为了城市3维模型构建的重要方法。已有一系列学者针对从LiDAR数据出发重建建筑物模型展开了研究并提出了一些方法。这些方法多数是先对LiDAR点进行面片分割，得到初始面片范围，而后进行精化。但由于LiDAR点间距约一米，精度有限，初始的面片分割得到的轮廓不可能是准确轮廓，对于这一问题的解决目前多采用规则化的方法，但规则化方法存在了一定问题：(1)初始面片的构建可能并不完全准确，局部凹凸多边形的出现可能使规则化后的结果出现多余多边形或者局部特征的丢失；(2)单个面片的规则化，对各个面片的拓扑关系的考虑有所欠缺，这样的规则化容易导致轮廓线的偏移。

单纯以影像数据或者是LiDAR数据进行建筑物的自动重建都存在着一定的问题，为了获取更好的建模效果可以考虑二者的结合。目前，已经有一些LiDAR与高分辨率影像结合的建筑物重建方法。2003年，Sohn 等人在International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences of the ISPRS会议中撰文《Building extraction using Lidar DEMs and Ikonos images》，将LiDAR点与IKONOS影像相结合，通过滤波将地面点与非地面点分离，再利用影像计算NDVI指数，剔除植被数据。对建筑物点，使用模型驱动和数据驱动相结合的方法获取轮廓线，最后完成模型的重建。这一方法在建筑物轮廓线的提取上比较独特，但是由于其使用了固定的阈值，对于大区域的数据来说不够灵活。2004年，Chen等人在International archives of Photogrammetry and Remote Sensing会议中撰文《Fusion of LIDAR data and optical imagery for building modeling》，提出从LiDAR点云中检测边缘，再结合高分辨率航空影像，提高边缘的几何精度，组合成闭合的多边形，最后进行模型的重建。这种方法能够很好地修正建筑物的边缘，但建模的精度很大程度依赖于影像分割的精度，未能充分发挥LiDAR点的数据优势。

LiDAR数据与影像之间是有很强的互补性特征，两者的结合会对重建带来显著的改善。

发明内容

1.发明要解决的技术问题

本发明的目的是提供一种自动化构建3维屋顶模型的方法，该方法根据LiDAR数据进行屋顶面片的分割以及各面片噪声的剔除，再进行屋顶面片的拟合，提高了建模的正确性与完整性，同时结合LiDAR数据与正射影像提取精确的屋脊线、轮廓线，提高了建模的精确性，最后进行3维屋顶模型的自动化构建，实现了建模的自动化。

2. 本发明的技术方案如下：

原理：对LiDAR数据构建Delaunay三角网，利用LiDAR点提供的高程信息使用三角形簇的方法（三角形簇的方法即根据高程信息把邻接的三角形集合组合成不同的组，构成三角形簇，并采用TNI结构形式，绘制时以三角形簇为单位进行调度和构网）提取出屋顶面片；利用高分影像的高分辨率特性提取精确边缘线并构建左右多边形，根据两侧多边形内LiDAR点所属面片情况进行屋脊线、轮廓线的提取；根据屋顶面片内的LiDAR点进行屋顶面片拟合，结合屋脊线、轮廓线进行3维屋顶模型的构建。

一种基于LiDAR数据与正射影像的3维屋顶重建方法，主要包括3个关键步骤：（参照附图1）

步骤1：基于三角形簇的LiDAR点屋顶面片分割，即对LiDAR点构建Delaunay三角网，根据LiDAR点的高程使用三角形簇的方法进行屋顶面片的分割。具体步骤为：（参照附图2）

（1）对LiDAR数据进行预处理，包括重采样和滤波。重采样时把邻阈内最高值赋予采样点的高程，以便有效剔除墙面噪声，保持数据的真实性，保留建筑物的轮廓信息。滤波时采用中值滤波，在剔除噪声的同时有效保留边缘信息。

（2）对预处理后的LiDAR数据构建Delaunay三角网。使用基于格网索引的生长算法进行三角网的实现，并与常规生长算法和基于拓扑结构和凸壳的生成算法进行了比较，比较发现基于格网索引的生长算法效率最高。

（3）对构建好的Delaunay三角网进行边缘线的筛选。传统的筛选是使用三角形法向量法进行筛选，即计算相邻三角形法向量之间的夹角，这种方法以相邻三角形为处理单元，处理单元过小，难以反映整体边缘情况，会保留大量噪声。本发明在三角形法向量基础上，进一步使用三角形簇的点法向量方法筛选，即把三角形簇内的三角形法向量算术平均作为该三角形簇的法向量，计算三角形簇法向量之间的夹角。经过以上处理后，能够剔除大量的非边缘线，但是会保留一部分位于墙面上的法向量变化较大的三角形边，对于这些三角形边，根据边两个端点的高差以及边与水平面之间的夹角进行筛选。

（4）对筛选出的面片边缘线进行面片的生成。Delaunay三角网经过边缘线筛选得到线是间断的且带有分叉的，利用这些线进行面片生成需要进行一定的简化和连接。根据Delaunay三角网边的相邻关系，把这些独立的三角网边组合成边缘线，剔除长度过短的线，同时对于带有分叉的线，保留两个较长的分叉，剔除细小分叉。最后进行边缘线的连接，找到边缘线的端点，再根据端点高程值，进行大致面片划分，高程差在10内厘米可认为属于同一个面片，这些间断点才可以连接，寻找两两最近端点相连接。

步骤2：基于LiDAR数据与正射影像的屋脊线提取，即对正射影像进行边缘检测，再对边缘线构建左右多边形，根据点面空间关系确定落在多边形内的LiDAR点，根据左右多边形内LiDAR点所属面片情况进行屋脊线的筛选，最后进行完整屋脊线的恢复，具体步骤为：（参照附图3）

（1）根据已有的建筑物区域的LiDAR数据确定影像工作区，套合LiDAR数据与正射影像，对LiDAR点内插并外扩，对外扩后的区域做最小外接矩形，以此切割过滤后的影像，使用外扩后的面过滤切割的影像，并得到线段提取的工作影像。

（2）对影响工作区进行边缘检测和Hough变换，获取矢量的边缘线数据。对影像工作区使用Edison边缘检测方法得到栅格化的边缘信息，经过Hough变换得到矢量化的线性边缘数据。

（3）对从影像上获取的边缘线进行屋脊线的提取。对所有边缘线构建左右多边形，根据点面空间关系确定落在左右多边形里面的LiDAR点，根据LiDAR点所属面片情况进行边缘线的筛选，得到属于屋脊线或者位于屋脊线附近的边缘线。

（4）对筛选出来的屋脊线进行编组及融合。在某些屋顶面片交界处存在多条线段。将相互平行且距离小于阈值的线段编组。当组内存在多条线段时，根据这些线段的端点拟合新的线段，完成屋脊线的融合处理。

（5）利用融合后的屋脊线进行屋脊线的恢复。对屋脊线进行延长，与轮廓线相交，把屋顶分割成多个细小多边形，确定落在这些多边形内LiDAR点。两两比较多边形的LiDAR点，如果所属面片基本相同则合并为同一个多边形。

步骤3：3维屋顶模型重建，即根据已得到的屋顶面片，确定落在各屋顶面片内的LiDAR点，拟合各个屋顶面片，再根据屋脊线和轮廓线确定的每个屋顶面片的精确范围，构建3维屋顶模型，具体步骤为：

（1）对于第一步骤得到的屋顶面片，根据点面空间关系，确定各面片内的LiDAR点，使用随机采样一致性方法拟合面片。

（2）把屋脊线与拟合的面片进行套合，确定面片与屋脊线的对应关系，根据拟合的面片计算屋脊线各角点的3维坐标。

（3）在ArcEngine中使用MultiPatch类进行3维模型的构建。

3. 有益效果

本发明集成LiDAR数据与高分辨率影像，综合利用LiDAR数据高程特性与影像高分辨率特性的互补优势，以“屋顶面片分割-屋脊线提取-3维屋顶模型重建”为主线，形成了基于三角形簇的LiDAR点屋顶面片分割方法、基于LiDAR数据与正射影像的屋脊线提取方法，并最终构成了一种3维屋顶模型重建的新方法。

相比于现有技术，该方法主要优势在于：

（1）本发明构建的3维屋顶模型正确性和完整性较高。基于三角形簇的LiDAR点屋顶面片分割过程中，使用三角形簇方法进行边缘线提取相比于三角形法向量方法，能够剔除更多的噪声，同时有效保留边缘信息，提取的结果经过高差、倾斜度筛选、线的简化和线的连接操作后，能够有效地形成建筑物各个面片，正确性和完整性较高；

（2）本发明构建的3维屋顶模型定位精度较高。基于LiDAR与正射影像的屋脊线提取过程中，从高分辨率影像中提取边缘线再进行屋脊线的筛选，充分发挥了高分影像的高分辨率特性以及LiDAR数据高程特性的优势，相比于规则化方法，精度较高；

（3）本发明的自动化程度高、建模速度快。在Delaunay三角网生成过程中，研究了常规生长算法、基于格网的生长算法和基于拓扑结构和凸壳的三角网生成算法，发现基于格网的生长算法效率最高，时间消耗最小。在屋脊线提取的过程中根据LiDAR进行影像工作区的建立，在影像工作区内进行后续处理，大量减少了工作量；

综上，本发明能够快速自动化地实现3维屋顶模型的构建，构建模型的正确性、完整性和定位精度较高，能够广泛应用于测绘遥感、虚拟现实、城市规划、科学研究等领域。

附图说明

附图1研究技术路线图；

附图2基于三角形簇的LiDAR点屋顶面片分割流程图；

附图3基于LiDAR数据与正射影像的屋脊线提取流程图；

附图4实施例的LiDAR数据与正射影像图；

附图5三角形簇方法提取边缘线结果图；

附图6高差和倾斜度筛选结果图；

附图7分叉处理结果图；

附图8面片分割结果图；

附图9影响工作区的建立示意图；

附图10影像线段提取结果图；

附图11屋脊线筛选结果；

附图12屋脊线融合结果；

附图13完整屋脊线的恢复；

附图14分割后屋顶面片与屋脊线数据套合；

附图15建模结果图；

附图16为3维屋顶模型平面误差向量图。

具体实施方式

以下通过实施例对本发明作进一步描述。

实施例

以附图4所示的LiDAR数据与正射影像的两个实验区为例说明。正射影像数据分辨率为5cm，LiDAR数据以不同的灰度值来表示不同的高程，平均点间距为1m，高程精度为15cm，平面精度约为30cm。实验区内包含了不同屋顶结构、不同走向、不同复杂程度的多个建筑物。

参照技术方案描述的三个步骤来对LiDAR数据和正射影像进行处理，具体流程步骤如下：

步骤一：基于三角形簇的LiDAR点屋顶面片分割

（1）对LiDAR数据进行预处理，包括重采样和滤波。

①重采样

常用重采样方法包括最邻近法、双线性插值法和三次卷积方法，不过这三种方法都不太适合LiDAR数据。最邻近法重采样会保留大量墙面点，带来大量的噪声；使用双线性插值法和三次卷积的方法重采样的结果较为平滑，淡化建筑物边缘信息。这里在计算采样点的高程时，把邻域内最高点的Z值赋给采样点，有效剔除墙面点信息，保持数据真实性，保留轮廓信息，重采样间隔为原始数据的平均点间距1m。

②滤波

由于LiDAR数据存在大量噪声，为了有效剔除噪声，同时保留建筑物的边缘信息，采用中值滤波方法，这里选用3×3的滤波器进行滤波。

（2）对预处理后的LiDAR数据构建Delaunay三角网。

从LiDAR点提取建筑物面片主要的运算量在于TIN的生成，其生成效率直接关系到整个自动化提取的效率，所以有必要选择一种效率最高的生成算法。发明中采用了基于格网索引的生长算法进行Delaunay三角网的实现，并把该方法与常规生长算法以及基于拓扑和凸壳的生成算法进行了对比，比较结果如表1。

表1 TIN生成算法比较

Figure 2011104235604100002DEST_PATH_IMAGE001

（3）对构建好的Delaunay三角网进行边缘线的筛选。

从Delaunay三角网中筛选边缘线通常做法是使用三角形法向量的方法，即计算相邻三角形法向量的夹角，计算方法如公式（1），式中Ni、Nj为相邻两三角形的法向量。

Figure 2011104235604100002DEST_PATH_IMAGE002

公式（1）

使用三角形法向量方法搜索变化边缘，搜索单元过小，难以反映整体的边缘情况，同时容易带来大量的多余边，给后面处理带来麻烦。

本发明在三角形法向量搜索的基础上扩大搜索范围，计算相邻顶点周围两个三角形簇的法向量的夹角。（参照附图5）具体做法是，首先计算出单个三角形的法向量，对于相邻顶点分别搜索到与其相邻的三角形，将这些三角形法向量平均作为三角形簇的法向量，最后计算相邻三角形簇的夹角。

核心代码如下：

public double CalAngleESOD(Edge edTemp, List<Pdata> pd, List<Edge> ed, List<Tri> tr,int Ctri)

{

#region //计算左边三角形簇的法向量

int[] ids = new int[3];

Get3Vertex(tr[edTemp.LeftTri], ref ids, ed, tr); //得到每个三角形的顶点

int N = 0;

for (int i = 0; i < 3; i++)

if ((ids[i] - edTemp.LeftTri) * (ids[i] - edTemp.RightTri) != 0)

N = ids[i]; //左边三角形对应的顶点

List<int> ringIDs = new List<int>();

for(int i = 0; i < Ctri;i++)

if ((ed[tr[i].Et[0]].Start - N) * (ed[tr[i].Et[0]].End - N) * (ed[tr[i].Et[1]].Start - N) *(ed[tr[i].Et[1]].End - N) * (ed[tr[i].Et[2]].Start - N) * (ed[tr[i].Et[2]].End - N) == 0)

ringIDs.Add(i); //得到邻接三角形

Vector leftTri = new Vector(); //左边图形的法向量

for (int i = 0; i < ringIDs.Count; i++)

{

ids = new int[3];

Get3Vertex(tr[ringIDs[i]], ref ids, ed, tr); //得到每个三角形的顶点

Vector a = new Vector();

a.x = pd[ids[0]].X - pd[ids[1]].X;

a.y = pd[ids[0]].Y - pd[ids[1]].Y;

a.z = pd[ids[0]].Z - pd[ids[1]].Z;

Vector b = new Vector();

b.x = pd[ids[1]].X - pd[ids[2]].X;

b.y = pd[ids[1]].Y - pd[ids[2]].Y;

b.z = pd[ids[1]].Z - pd[ids[2]].Z;

Vector ab = new Vector(); //三角形法向量

ab.x = a.y * b.z - a.z * b.y;

ab.y = a.x * b.z - a.z * b.x;

ab.z = a.x * b.y - a.y * b.x;

leftTri.x = leftTri.x + ab.x;

leftTri.y = leftTri.y + ab.y;

leftTri.z = leftTri.z + ab.z;

}

if (ringIDs.Count != 0)

{

leftTri.x = leftTri.x / ringIDs.Count;

leftTri.y = leftTri.y / ringIDs.Count;

leftTri.z = leftTri.z / ringIDs.Count;

}

#endregion

Vector rightTri = new Vector(); //左边图形的法向量

//类似方法计算右三角形簇的法向量，省略

return CalVectorAngle(leftTri, rightTri); //计算两法向量之间的夹角

}

程序中主要构建了PData类、Edge类、Tri类。PData类是存储处理LiDAR点数据的类，主要包括了点的编号、坐标以及点所属屋顶面片信息等。Edge类是处理线的类，主要包括了线的编号、起始LiDAR点的编号、终止LiDAR点的编号、Delaunay三角网中左右三角形的编号等。Tri类是存储处理Delaunay三角形的类，主要包括了三角形的编号、构成三角形的边的编号、三角形顶点等。

使用以上方法计算得到三角形簇的法向量，根据处理数据的建筑物屋顶特点确定阈值，这里阈值取10°，保留夹角大于阈值的两三角形簇的公共边，处理结果见附图5。

对Delaunay三角网使用三角形簇方法筛选后，除了保留了建筑物屋顶面片的边缘线外，还保留了墙面上一些角度变化较大的线，计算边缘线两端点之间的高差及其倾斜度，根据建筑物屋顶特征确定阈值，这里高差阈值设为0.3m，剔除高差大于阈值的线，倾斜度阈值取10°，剔除倾斜度大于阈值的线。处理结果见附图6。

（4）对筛选出的面片边缘线进行面片的生成。

使用三角形簇方法从Delaunay三角网中筛选边缘线，得到的结果是离散的、缺失拓扑关系的线，从这些线段中进行面片的生成还需要进行以下处理。

首先根据线段之间顶点的邻接关系，把相邻接的线段合并成一条折线段。具体做法是：构建折线段链表，对于链表内的每一条Delaunay边缘线，遍历提取出来的边缘线，判断链表中的边缘线是否与遍历所用的边缘线具有相同端点，如果有则把该边缘线添加到折线段链表中。

将相邻边缘线合并成同一折线段后，需要根据折线段所包含的边缘线的数量剔除数量较少的折线段。经过以上处理得到的折线段可能是带有分叉的。对于分叉边的处理，首先找到各个分叉点，即周围有超过两条边缘线的点，然后计算从分叉点沿着每条边缘线到折线段端点或者其他分叉点的边缘线的数量，保留数量最多的两条边，其余的删除。处理结果见附图7。

经过以上处理后，已经能够提取出边缘线了，但是这些线段是间断的，并不构成封闭的面，为了得到各个面片，还需要连接这些间断线形成各个面片。具体的做法分为三步，首先找到各个间断点，即只与一条边相连的点，再根据这些点的高程值，进行大致面片划分，认为高程差在10cm范围内，才属于同一个面片，这些端点才可以连接，最后对于被划分到各个面片的端点，寻找两两最近的端点并进行连接。面片分割结果见附图8。

步骤二：基于LiDAR数据与正射影像的屋脊线提取

（1）根据已有的建筑物区域的LiDAR数据确定影像工作区。

影像工作区的建立是指从整个原始影像中获取建筑物区域的影像，进行影像工作区的建立，使得后续的处理都在工作区影像中进行，能够有效地减少工作量，提高系统效率。影像工作区建立具体过程如下：套合LiDAR数据与影像数据，即根据两套数据的坐标位置进行叠合；对LiDAR点进行向内插值向外扩张，即对LiDAR点进行插值，对插值的结果进行一定的扩张，以保证所有轮廓信息都在影像范围内，阈值取LiDAR数据平均点距的2倍；根据外扩的多边形生成外接矩形，使用该矩形切割影像；最后使用外扩的多边形过滤切割后的影像，即可得到工作区影像（见附图9）。

（2）对影响工作区进行边缘检测和Hough变换。

利用Edison算子，作用于工作影像，进行边缘检测。在边缘检测的基础上，利用Hough变换，将栅格化的边缘线转化为矢量边缘，提取结果如附图10所示。

（3）对从影像上获取的边缘线，进行屋脊线筛选。

对获取的边缘线进行屋脊线的筛选，首先对每条边缘线垂直于线段方向做左右多边形，通过点面空间关系确定落在多边形内的LiDAR点，根据线段与LiDAR点的关系来筛选屋脊线，方法如下：

①若左右多边形中一个有点而另一个没有点，则说明该线段是建筑物外围轮廓线，使用建筑物外围轮廓线生成方法，该方法已有很多学者进行了详细的介绍和实现，此处不再赘述。

②若左右矩形的LiDAR都来自于同一个屋顶面片，则说明该线段不在屋脊线附近，则剔除。

③若左右多边形来自于不同的点集，则有三种情况。左矩形内点全部来自屋顶1且右矩形中点全部来自屋顶2，说明线段准确的落在屋顶1、2之间的屋脊上，线段保留；一侧矩形内点全部来自某一屋顶，而另一侧矩形内点来自两个屋顶，则说明线段处于屋顶面片交界处，但与精确屋脊有一定的距离偏差，计算该侧矩形内两个面片点的比例，若小于阈值（20%），保留线段，否则剔除；两侧矩形内点均来自两个屋顶，则说明线段处于屋顶面片交界处，但与精确屋脊有一定的方向偏差，依上述计算比例的方法，若小于阈值（20%）,保留线段，否则剔除。

屋脊线筛选结果参照附图11。

（4）对筛选出来的屋脊线进行编组以及屋脊线的融合。

根据屋脊线筛选结果，如附图11，在某些屋顶面片交界处（如，1与2），只有一条线段，在某些屋顶面片交界处（如，3与4、3与5、4与5、7与8），存在多条线段。将相互平行且距离小于阈值（垂直距离小于1m、平行距离小于2m）的线段编组。当组内存在多条线段时，根据这些线段的端点，利用随机采样一致性算法拟合出一条新的线段，完成屋脊线的融合处理，得到如附图12中的线段。

（5）对融合后的屋脊线进行屋脊线的恢复。

融合后的屋脊线是断裂而非连续闭合的，拓扑关系丢失，使用“LiDAR数据支撑下改进的分割－合并－成型算法”，完成其拓扑关系恢复。具体做法是，将所有屋脊线线段自延长，其与轮廓线将建筑物分割成多个多边形；分析每个多边形内LiDAR点的归属，如果相邻多边形内的LiDAR点来源于同一面片，则合并多边形，即完成了屋脊线的拓扑关系重建，也即实现了完整屋脊线的提取（见附图13）。

步骤三：3维屋顶模型重建

（1）对屋顶各面片进行面片拟合。

对于分割好的屋顶面片数据，根据点面空间关系确定落在各个面片内的LiDAR点，根据这些点拟合屋顶各个面片，这里采用随机采样一致性算法进行拟合。随机抽样一致性方法只使用正确的数据进行拟合，忽略少数错误的数据，对于有较多误差的数据可以经过多次迭代得到较好的结果。

RANSAC算法过程如下：

①选取一定的模型，这里采用最小二乘模型；

②在原始数据中随机选取一定数量的初始点作为局内点，并使用这些点进行最小二乘拟合；

③使用新生成的拟合模型对剩余的原始数据点进行误差检验，如果误差小于给定的阈值，则把该数据点加入局内点；

④对所有数据检验完毕后如果局内点的点数大于一定的数量，则认为该模型是合理的，使用这些点重新生成模型，并把该模型记入最优模型，进行下一次迭代；

⑤在迭代的过程中当有新的最优模型产生时，比较其与先前的最优模型的优劣，评估规则是局内点与模型的错误率；

（2）计算屋脊线各角点的3维坐标。

拟合屋顶面片后就可以根据各个面片计算屋脊线各角点的3维坐标。首先将屋脊线数据与屋顶面片数据套合（如附图14所示），计算每个分割的屋顶面片与屋脊线形成的面片相交的面积，如果相交面积占屋顶面片比例达到85%，则说明该屋脊线确定的面片对应于此屋顶面片，计算屋脊线高程时就使用该屋顶面片的拟合方程计算。核心代码如下：

while (pFeature1 != null)

{

IGeometry pPolygon1 = (IPolygon)pFeature1.Shape;

IFeatureClass pFeatureClass2 = pFeatureLayer2.FeatureClass;

ISpatialFilter pSpatialFilter = new SpatialFilterClass();

pSpatialFilter.Geometry = pFeature1.Shape;

pSpatialFilter.GeometryField = pFeatureClass2.ShapeFieldName;

pSpatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;

IFeatureCursor pFeatureCursor2 = pFeatureClass2.Search(pSpatialFilter, false);

IFeature pFeature2 = pFeatureCursor2.NextFeature();

while (pFeature2 != null)

{

IGeometry pPolygon2 = (IPolygon)pFeature2.Shape;

IGeometry pPolygon = new PolygonClass();

pTopoOper = (ITopologicalOperator)pPolygon1;

pPolygon = pTopoOper.Intersect(pPolygon2,esriGeometryDimension.esriGeometryNoDimension); //计算相交部分

if (pPolygon.IsEmpty == false)

AddPolygonFeature(pFeatureClass, (IPolygon)pPolygon);

CalArea(pFeatureClass, 3); //计算相交部分面积

pFeature2 = pFeatureCursor2.NextFeature();

}

pFeature1 = pFeatureCursor1.NextFeature();

}

确定了屋顶面片与屋脊线面片之间的对应关系，根据屋顶面片的拟合方程以及屋脊线角点的平面坐标，就可以计算得到屋脊线角点的3维坐标。

（3）在ArcEngine中使用MultiPatch类进行3维模型的构建。

ArcEngine的MultiPatch类提供了3维建模的方法，这里主要用到的类库包括：IGeometryCollection、IMultiPatch、IPointCollection、IRing、IZAware等。具体过程如下：对于每个屋顶面片，首先根据屋脊线角点3维坐标构建3维屋顶面片，再根据角点与其在地面上的投影点构建各个墙面，最后把屋顶面与墙面相结合构成房屋3维模型。核心代码如下：

IGeometryCollection multiPatchGeometryCollection = new MultiPatchClass();

IMultiPatch multiPatch = multiPatchGeometryCollection as IMultiPatch;

IPointCollection topPointCollection = new RingClass();

for (int j = 0; j <= planeList[i].PointIDs.Count; j++) //构造屋顶面片

{

ESRI.ArcGIS.Geometry.IPoint point = new ESRI.ArcGIS.Geometry.PointClass();

IZAware pointZA = point as IZAware;

pointZA.ZAware = true;

point.X = pd[planeList[i].PointIDs[j % planeList[i].PointIDs.Count]].XData;

point.Y = pd[planeList[i].PointIDs[j % planeList[i].PointIDs.Count]].YData;

point.Z = pd[planeList[i].PointIDs[j % planeList[i].PointIDs.Count]].Z;

topPointCollection.AddPoint(point, ref missing, ref missing);

}

multiPatchGeometryCollection.AddGeometry(topPointCollection as IGeometry, ref missing, ref missing);

multiPatch.PutRingType(topPointCollection as IRing, esriMultiPatchRingType.esriMultiPatchOuterRing);

for (int j = 0; j < planeList[i].PointIDs.Count; j++) //构建各个墙面

{

IPointCollection wallPointCollection = new RingClass();

ESRI.ArcGIS.Geometry.IPoint point1 = new ESRI.ArcGIS.Geometry.PointClass();

IZAware pointZA1 = point1 as IZAware;

pointZA1.ZAware = true;

point1.X = pd[planeList[i].PointIDs[j % planeList[i].PointIDs.Count]].XData;

point1.Y = pd[planeList[i].PointIDs[j % planeList[i].PointIDs.Count]].YData;

point1.Z = pd[planeList[i].PointIDs[j % planeList[i].PointIDs.Count]].Z;

wallPointCollection.AddPoint(point1, ref missing, ref missing); //加入第一起始墙面点

//类似于加入起始墙面点，加入墙面上其他的墙面点，省略

wallPointCollection.AddPoint(point1, ref missing, ref missing); //再加入起始墙面点

multiPatchGeometryCollection.AddGeometry(wallPointCollection as IGeometry, ref missing, ref missing);

multiPatch.PutRingType(wallPointCollection as IRing, esriMultiPatchRingType.esriMultiPatchOuterRing);

}

结果

使用本发明方法分别对如附图4所示的实验区域1和实验区域2进行了处理，建模结果如附图15。

以3维模型的每块屋顶面片为单位，根据公式2，评价所重建3维屋顶模型准确性与完整性。

Figure 2011104235604100002DEST_PATH_IMAGE003

式中TP指屋顶面片重建正确的数目，FN指屋顶面片遗漏重建的数目，FP指屋顶面片重建错误的数目，Completeness指所提取屋顶与实际屋顶相比较的完整程度，Correctness指所提取的屋顶与实际情况相比的正确程度。对于屋顶面片重建正确、遗漏、错误的判定如下：以人工建模结果为参考数据，将所重建屋顶面片与人工建模结果投影到二维平面分别形成多边形A、B，获取两多边形的重叠区域C，计算重叠区域C与多边形A的面积比值，若大于80%，则认为屋顶面片重建正确，否则视为重建错误。遗漏是指本身存在的屋顶却未被重建。区域1、2中所重建3维屋顶模型的正确率与完整率数据见表2。

表2 模型重建完整性与准确性表

实验区	正确	错误	遗漏	正确率	完整率
						区域1	52	7	9	86%	85%
区域2	38	4	5	90%	88%

3维屋顶模型的定位精度包括平面精度与高程精度，这里重点考察其平面精度。在实验区的区域1、2内，通过评价屋顶面片各角点的平面精度，来评价模型总体的精度。精度评价中分别采用平均误差、均方根误差和最大偏移三项指标。从重建正确的屋顶面片入手，挑选出一定数量的屋顶角点（区域1、2分别挑选25、20个角点），以真正射影像图为地面真值，将其与真正射影像图上的对应点进行比较，衡量本文方案所重建3维屋顶模型的定位精度。精度评价结果如表3所示，屋顶角点的平面误差向量图如附图16所示（图中加粗短线表示屋顶角点定位误差，其方向表示误差方向，其长度表示误差大小）。

表3 屋顶模型重建的平面精度

实验区	平均误差(m)	均方根误差(m)	最大偏移(m)
				区域1	0.35	0.39	0.74
区域2	0.28	0.31	0.63

Claims

1.一种基于LiDAR数据与正射影像的3维屋顶重建方法，包括以下步骤：

步骤1：基于三角形簇的LiDAR点屋顶面片分割，首先对LiDAR数据进行预处理，包括重采样和中值滤波；对于预处理后的LiDAR点，使用基于格网索引的生长算法，实现对大量LiDAR点的不规则三角网的快速生成；对于生成的Delaunay三角网，在三角形法向量法进行筛选的基础上，进一步采用基于三角形簇的点法向量估计技术，通过相邻点法向量的差异分析，探测建筑物顶部变化之处，并结合三角网边的高差和倾斜度信息剔除部分噪声；为了从复杂的线数据中有效形成屋顶各个面片，采用基于分叉点的简化方法以及基于间断点的连接方法，最终生成屋顶面片；

步骤2：基于LiDAR数据与正射影像的屋脊线提取，把LiDAR数据与影像套合，根据LiDAR内插外扩的结果进行影像工作区的建立；利用Edison算子和Hough变换提取建筑物范围内的边缘线数据；通过对每条边缘线数据构建左右多边形，根据点面空间关系确定落在左右多边形内的LiDAR点，对比左右多边形的LiDAR点分布情况，选取出位于建筑物屋脊线附近的线；为了进一步简化屋脊线，对于近似平行且距离较近的屋脊线进行融合；通过屋脊线的自延伸与轮廓线相交构成一个个面片，两两比较面片内的LiDAR点，若都来自于同一屋顶面片，则合并两个面，最终可以得到准确的屋脊线数据；

步骤3：3维屋顶模型重建，获取落在每个屋顶面片内的LiDAR点，使用随机采样一致性算法拟合3维屋顶面片；屋顶面片数据与屋脊线面片数据套合，实现两种面片数据的一一对应，利用拟合面片方程以及屋脊线角点的平面坐标计算屋脊线角点的3维坐标；使用ArcEngine的MultiPatch类实现3维屋顶模型的重建。

2.根据权利要求1所述的一种基于LiDAR数据与正射影像的3维屋顶重建方法，其特征在于步骤1中进行的以下处理：在预处理重采样的过程中，采样点的高程被赋予了邻阈内最高值，有效剔除了部分墙面噪声。