倾斜影像中提取高精度DEM的方法研究

李 莹，苏国中，林宗坚

(1. 辽宁工程技术大学测绘与地理科学学院，辽宁 阜新 123000； 2. 中国测绘科学研究院地理信息工程国家测绘地理信息局重点实验室，北京 100039)

倾斜影像中提取高精度DEM的方法研究

李 莹1，2，苏国中2，林宗坚2

(1. 辽宁工程技术大学测绘与地理科学学院，辽宁 阜新 123000； 2. 中国测绘科学研究院地理信息工程国家测绘地理信息局重点实验室，北京 100039)

针对传统航空影像获取的DSM在立面及局部地面、建筑物屋顶空间信息的不足，获取高精度DEM较为困难的问题，提出了基于倾斜影像提取高精度DEM的方法。首先对倾斜影像获取的点云DSM结构进行分析，得出了DSM具有几何约束特点，能够在城区很好地区分地面点和地物点；然后指出对DSM滤波处理是获取高精度数字高程模型(DEM)的关键技术，提出了基于法向量差值区域生长分割TIN的滤波方法；最后选取吉林省敦化市的倾斜影像数据进行了滤波试验和算法验证。试验结果表明，该方法能够快速、有效地滤除不同尺寸的建筑物、植被和其他地物，获取高精度DEM。

倾斜影像；DSM；TIN；区域增长；滤波；高精度DEM

在当今卫星遥感和普通航空遥感蓬勃发展的形势下，轻小型航空遥感，即无人飞行器低空航测系统(UAV-MAP)是不可缺少的补充技术手段[1-3]。传统的UAV上装载单相机进行航测技术获取的数字地表模型(DSM)，通常为2.5维，尤其在城市地区，很难区分地物点和地面点，获取数字高程模型(DEM)较为困难。但是UAV装载组合宽角相机拍摄的倾斜影像，利用4个不同视图方向DSM来弥补垂直DSM在立面及局部地面、建筑物屋顶面空间信息的不足，来获取全方位的DSM，很容易检测和剔除建筑物等高出地面的地物，获得真实地形DEM。

倾斜影像获取的DSM具有几何约束特点，影像密集匹配生成全方位点云DSM。利用已知良好匹配点建立不规则三角网格(TIN)，以实现纹理自适应的匹配传播，并利用新匹配点对TIN进行动态细化，有利于减少因纹理贫乏因素而导致的错误匹配，提高DSM的精度；这种点、三角网及纹理几何约束的关系，决定了DSM具有一定几何约束特点。

滤波是DSM获取高精度DEM的关键技术[4-13]，必须从DSM检测和剔除建筑物等高出地面的地物，获得真实地形DEM。现有滤波的方法很多，由于数据源的不同，滤波方法也大不相同。根据DSM的三角网格具有几何约束特点，需要在基于TIN模型上进行滤波处理。在基于TIN算法用于点云数据滤波原理方面，Krzystek[14]建议,先获取局部区域最低点构建一个粗糙的TIN凸壳,以高差临界值判别地面点,利用有限元分析法调整TIN模型,该算法在具有不同树木结构的森林地区得到了较好的应用。Haugerud[15]直接用原始点云数据建立了TIN模型,然后通过TIN模型计算地形曲率参数,认为曲率急剧变化的点为非地面点，该方法不适合地形起伏地形。以上方法侧重于解决特定地形条件下的点云滤波问题，每一个算法都有局限性，不能准确有效地过滤建筑物、植被及其他地物等。因此，以倾斜影像获取的点云DSM为主要研究对象，本文在基于TIN的分割方法的基础之上，提出基于法向量差值区域生长分割TIN的滤波算法，该算法首先以规则格网将TIN进行分块，进而在每个区域内将若干个地面点和地物点作为种子，用区域增长的方法进行生长，计算三角形的法向量，认为法向量变化到一定阈值时进行分割；其次利用分割段高程、法向量作为约束条件，利用分割后的点云实现了建筑物、植被和其他地物的识别分类，并在获得地面点后将其插值生成数字高程模型。

1 基于区域生长分割TIN的滤波算法

1.1 等格网DSM分块

DSM是点云构成的不规则三角网格,数据量十分庞大。通常在滤波之前，考虑将复杂的整体DSM进行分块，采用“先分块后融合”的思路构建子网格，最后融合子网格恢复原始实体的表面形状。这样就能把复杂的问题简单化，提高了提取DEM的效率。

等格网DSM分块的基本思路是：设定固定格网大小，对DSM进行规则格网划分，根据格网单元的范围确定每个区域中的三维点，即构建数据索引，从而实现对DSM数据的分块。如图1所示，实现了基于邻域关系的快速查找，在各自网格内进行区域生长，从而提高滤波算法的效率，便于区域生长分割TIN。

图1

1.2 基于区域生长分割TIN方法

本文采用基于区域生长分割TIN方法。区域生长的基本思想是：首先将三角网内的三角形进行拓扑构建。采用自动选取种子的方式，在三角网平坦的地方，依据三角面片的几何原则，选择种子点；以一个种子点开始，根据生长原则选出可以进行生长的三角面片进行区域生长，把未在任何类中的三角面片元素放入种子类中，反复运行直到种子队列为空，停止分割，实现了三角网分割成段。

图2为区域生长过程示意图，图中颜色部分与非颜色部分分别代表区域生长过程中的已生长区域与未生长区域。该算法的流程如图3所示。

图2 区域生长过程

图3 算法流程

1.2.1 算法原理

(1) 种子选取：种子点的选择是进行区域生长的关键前提，对最后的分割结果影响很大。如果种子点选择不当，将导致整个分割过程的判断依据失真，最终得出错误的分割结果。因此，选择种子的几何原则如下：

三角形Tt和其所有相关三角形中的一个三角形Ti法向量的差值为

(1)

三角形Tt和其所有相关三角形法向量的期望为E(Tt)，方差为D(Tt),如下

(2)

(3)

保存一个三角形和其所有相关三角形法向量的最大值与最小值的差值记为

(4)

(2) 区域生长原理：DSM一般是由地面、建筑物、植被及其他地物组成。因此，TIN模型中只要存在两个面法向量差值较大的明显不是一类，根据DSM特点，确定相邻三角面片法向量差值作为区域生长的驱动信号。

1.2.2 算法流程

(1) 构建点云数据三角网格的拓扑关系，使得点、线及三角形有一定的拓扑关系。

(2) 种子的选择采用候选种子的方法，考虑面平坦的地方，根据一定的几何原则，将符合要求的三角形压入一个候选种子队列中；然后在候选种子队列中从队列头开始查找，如果已经在类中则弹出，直到找到一个没有在任何类中的三角形，把它选为种子，把种子初始化为一个类。

(3) 从类的边界三角形队列中，根据生长原则，开始生长。

(4) 同样根据生长原则选出类中边界三角形的相关三角形进行生长，在生长的同时，记录下生长中类的新边界，并将它们压入到一个暂时队列中。当类边界三角形队列为空时，把类中边界三角形队列更新。

(5) 一个类边界三角形的队列为空时，如果类中三角形总体的数目同生长前没有发生变化，类就停止生长。

(6) 如果没有达到停止条件，转到步骤(2)。如果达到了停止的条件，生成了一个类，并把它压入存放类的队列中，并转到步骤1。

(7) 反复顺序执行步骤(1)—步骤(5)直到候选种子队列为空。把没有在任何类中的三角形按照期望从小到大的顺序压入种子候选队列二中，再执行步骤(1)—步骤(5)，直到候选种子队列二的队列为空，结束分割。

1.3 基于聚类分析的DSM分类

由于基于法向量差值分割TIN具有局限性，还是可能会出现同一特征被划分为不同区域的过分割现象。在区域生长分割后，会由于过分割现象、模型误差而产生小区域，因此，根据相邻分割段法向量的朝向小于一定的阈值，进行合并处理。通过聚类分析的方法，将合并后的分割段，根据法向量、高度的约束条件，将TIN分类成地面、建筑物、植被及其他地物，方便以后的过滤。

1.3.1 建筑物判别

由于地面部分相对平坦，几乎没有突变；绝大部分建筑物的顶部是一个平面。根据分割段与段之间的几何拓扑，首先按照相邻分割段的平均法向量的夹角φ区分建筑物的墙面。类的平均法向量为P1(A1,B1,C1),相邻类的平均法向量为P2(A2,B2,C2)，则两个法向量之间的夹角为

(5)即φ>φ1或φ<φ2为墙面值。(其中，φ1=0.392 5，φ2=1.777 5为经验值)。

屋顶和地面均相对平坦，很难区分。因此，需要根据墙面特征来判断屋顶，判断的条件为：参照图6为例，以一个墙面开始，找到该墙面的相邻分割段，从中找出高出该墙面平均高度的分割段，则该分割段为屋顶；低于墙面平均高度的分割段为地面点。

1.3.2 植被判别

植被通常表现为相邻点之间空间分布离散，且有较大的高差，因此植被通常会被分为破碎的分割段。这些分割段面积小、数量多、高差大，检测到这种分割段则认为是植被。

对于有一定坡度地面上的植被来说，由于通常滤波方法都是基于地物点明显高于地面点实现的，因此这种假设对于地形起伏较大时，地形点可能和植被点具有相同的高度从而造成滤波失败。当然，低矮的植被是低矮地物的一种，与低矮地物一样，仅仅依靠几何信息时，往往难以滤除。对于这种难以分辨植被的情况，采取人工选择地面种子点，进行区域生长，获取地面点。

1.3.3 其他地物判别

按照尺寸和几何特征检测出建筑物和植被，剩下的基本是地面点，但是也包含一些附在地面的低矮地物，如小汽车、人等。其他地物分割段之间的所有可能总结如下：①待定分割段比它周围的分割段都要高；②待定分割段面积小。

经过上述分类，将TIN分为地面点和地物点。将获取的地面点直接进行三角剖分，得到高精度DEM。

2 试验及分析

2.1 滤波试验

本文采用C++来实现基于分割TIN的滤波算法，获取高精度DEM，并用若干个复杂城市地形的倾斜影像测试算法。下面是其中一个试验区的试验结果及分析。

该试验区地形复杂，右下角是一个明显高低起伏的山地，整个区域覆盖了不同高度的植被，一个明显的居民区位于左部分，居民区内的建筑物大小不一。试验数据为吉林省敦化市的倾斜影像匹配的DSM数据，大小为1000 m×600 m。点云数据如图4所示，TIN数据如图5所示。按照本文的算法，在全区范围内选取若干个种子点，取相邻三角形法向量差值进行区域生长分割三角网，根据分割段的平均法向量和平均高度将三角网分类成地面、建筑物和其他地物，图6为分类的晕渲图；根据分类的结果，过滤地物点，获取地面点，图7为全区滤波后的晕渲图，该晕渲图的颜色与图6一致，即两幅图中相同的高程具有相同的颜色；图7将获取的地面点生成DEM，图8为内插DEM结果晕渲图。

图4 点云数据

图5 TIN数据

图6 分类的晕渲图

图7 全区滤波晕渲图

图8 内插DEM结果晕渲图

2.2 精度评定

DEM作为基础地理信息数据，其精度至关重要。DEM提取之后，需要对其质量进行精度评定。本试验采用检查点方法对DEM进行精度评定。

精度评定的目的是检验DEM的高程精度，常规检查点的方法是一个检查点对应图上一个点，这种一一对应的方法适用于数字测图的评定；本文是检验DEM的高程精度，因此需要一个检查点对应DEM局部位置，以进行评定。在已有DSM上选取5个GPS控制点作为检查点，高度为h；找到新生成的DEM所对应点的位置，在此附近选点(10～15个点)，高程值为hi，逐一和检查点进行比较，计算高程的平方差值(用RMSE表示)。检查点的实际高程与提取结果的均方误差值见表1，该方法能够很好地对DEM精度进行评定。

(6)

经过精度评定，最大的RMSE为0.45 m，能够有效地滤除地物点，获取DEM精度较高。

造成误差的原因有：①倾斜影像精度及匹配生成DSM方法的精度影响DEM的精度；②区域本身的众多地形因素如地形复杂度、地形坡度等影响DEM的精度；③滤波算法影响DEM的精度。

3 结束语

(1) 本文直接基于倾斜影像匹配生成的DSM进行滤波，采用分块的DSM数据组织形式和基于分割的滤波算法，提高滤波的精度和效率的同时，能有效滤除不同尺寸的建筑物、低矮的植被和其他地物，能较好地保持地形特征，快速提取高精度DEM。

(2) 本文对于大面积、密度低矮的植被，在提取地面点的同时，算法会将非常靠近地面的低矮植被误认为地面点，目前解决的方法是通过人工选择种子点进行区域增长获取地面点。如何动态调整植被法向量差值的阈值，以提高算法的自适应能力，需要进一步研究。

试验结果证明：本文的算法适合倾斜影像，执行效率高、算法稳定,能有效地滤除不同尺寸的建筑物、植被和其他地物,地面点丢失较少，可获取高精度的DEM，具有一定的可靠性与实用性。

[1] 桂德竹.基于组合宽角相机低空影像的城市建筑物三维模型构建研究[D].北京：中国矿业大学(北京),2010.

[2] 解斐斐.基于无人飞艇低空航测系统建筑物纹理获取与处理技术[D].武汉: 武汉大学,2014.

[3] 张平,刘怡，蒋红兵.基于倾斜摄影测量技术的“数字资阳”三维建模及精度评定[J].测绘,2014,37(3):115-118.

[4] 袁枫.机载 LiDAR 数据处理与土地利用分类研究[D].北京：中国矿业大学(北京)，2010.

[5] 张小红.机载激光扫描测高数据滤波及地物提取[D].武汉：武汉大学，2004.

[6] 胡永杰，程朋根，陈晓勇，等.机载激光雷达点云滤波算法分析与比较[J].测绘科学技术学报，2015，32(1):72-77.

[7] 隋立春，张熠斌，枊艳，等.基于改进的数学形态学算法的LiDAR点云数据滤波[J].测绘学报，2010,39(4)：390-396.

[8] 隋立春，张熠斌，张硕，等.基于渐进三角网的机载LiDAR点云数据滤波[J].武汉大学学报(信息科学版)，2011,36(10):1159-1162.

[9] 秦川. 无人机影像匹配点云滤波处理与三维重建[D].成都：西南交通大学,2015.

[10] 马树发．基于改进虚拟格网的记载数据的形态学滤波[D].西安：西安电子科技大学，2014.

[11] 杨必胜,董震,魏征,等. 从车载激光扫描数据中提取复杂建筑物立面的方法[J]. 测绘学报,2013,42(3):411-417.

[12] SITHOLE G. Segmentation and Classification of Airborne Laser Scanner Data[D].TU Delft: Aerospace Engineering，2005:35-68．

[13] VOSSELMAN G. Slope Based Filtering of Laser Altimetry Data[C]∥International Archives of the Photogrammetry,Remote Sensing and SpatialInformation Sciences.Amsterdam:[s.n.]，2000.

[14] KRZYSTEK P. Filtering of Laser Scanning Data in Forest Areas Using Finite Elements[C]∥ISPRS Working Group Ⅲ/3 Workshop/3-D Reconstruction from Airborne Lasers Cannerand InSAR Data.[S.l.]：Dredsen University of Technology,2003.

[15] HAUGERUD R，HARDING A. Some Algorithms for Virtual Deforestation(VDF) of LIDAR Topographic Survey Data [J]. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences,2001,34(W4):211-215.

Study on the Method of Extracting High Precision DEM from Oblique Image

LI Ying1，2，SU Guozhong2，LIN Zongjian2

(1. School of Geomatics,Liaoning Technical University, Fuxin 123000, China; 2. Key Laboratory of National Administration of Surveying,Mapping and Geoinformation, Chinese Academy of Surveying and Mapping, Beijing 100039, China)

For the DSM of the traditional aerial images acquisition, the lack information of the facade,the local floor and the roof of the building, the problem which difficult to obtain high precision DEM, the method of extracting high precision DEM from oblique image is proposed in this paper. Firstly, the structure of cloud point of oblique image acquisition is analyzed, the DSM has a geometric constraint characteristics which can well distinguish between ground points and ground feature points in the city. Then it points out that the filtering is the key technology to obtain the high precision digital elevation model (DEM), and proposes a filtering method based on region growing segmentation TIN.Finally, the filtering experiment and algorithm validation are carried out by the oblique image data of Dunhua City, Jilin Province.The experimental results show that the method can quickly and effectively remove the different sizes of buildings,vegetation and other objects, and obtain high precision DEM.

oblique image; DSM；TIN；regional growth；filtering；high precision DEM

李莹，苏国中，林宗坚.倾斜影像中提取高精度DEM的方法研究[J].测绘通报，2017(2)：30-34.

10.13474/j.cnki.11-2246.2017.0043.

2016-05-28；

2016-09-07

国家自然科学基金(41371425；41271451)

李 莹(1991—)，女，硕士生，主要从事倾斜摄影城市三维建模方面的研究。E-mail：ly89832623@sina.com

P23

A

0494-0911(2017)02-0030-05