罗为东,甘 淑,袁希平,高 莎,毕 瑞,袁新悦

(1.昆明理工大学国土资源工程学院,云南 昆明 650093；2.云南省高校高原山地空间信息测绘技术应用工程研究中心,云南 昆明 650093；3.滇西应用技术大学地球科学与工程技术学院,云南 大理671006)

1 引 言

近年来随着科学技术的发展,复杂山区地形地貌提取分析朝着精细化方向发展,如何高效、快速获取高精度地形地貌空间数据成为了分析地形地貌特征的关键因素。无人机航测系统在对生产地面表面模型(Digital Surface Model,DSM)、数字正射影像图(Digital Orthophoto Map,DOM)、构建数字高程模型(Digital Elevation Model,DEM)等具有快速、高效的特点。地形数学产品均是对地形信息的详细表达。其中数字高程模型的快速获取,并满足相应的精度要求,特别是在困难山区对数字高程模型的处理,是山区数字化发展的重要趋势。目前,国内外学者利用无人机遥感技术采集构建DEM对地形地貌等空间特征分析展开了研究。不难发现,构建DEM对于无人机航测技术,主要是要对密集影像匹配(Dense image matching,DIM)点云进行滤波处理,对此GordanaJakovljevic[1]等通过LiDAR点云和密集匹配点云进行滤波处理,构建DEM用于降低洪水危险,并计算DEM的RMSE,结果表明平原地区适宜应用LiDAR采集点云,而丘陵和山区密集匹配点云构建DEM的RMSE低于LiDAR点云。Mustafa Zeybek[2]等应用四种滤波算法结果提取地形特征,结果表明粗糙度和复杂地貌提取四种滤波算法表现出相似性,CSF滤波算法在平坦地形区分地面点和非地面点正确率达到93 %。叶立志[3]针对DIM点云本身的复杂性,提出了一种基于交叉线元分割的密集匹配点云滤波算法,结果表明滤波结果构建DEM中误差接近1m,相对于商业软件Terrasolid的滤波结果和滤波效率更优。JiXian[4]等应用ISPRS第三个工作组的七个数据集,首先构建完初试TIN后,再SUSC扩展得到更多地面点,最后进行对TIN进行迭代加密,结果表明与PTD相比,所提出的方法能够预景观的不连续性,并将遗漏误差和总误差分别减少了大约10 %和6 %,降低了后期校正结果所需的人工操作成本。ShengNie[5]等对PTD滤波算法中建立TIN进行改进并改变原有的迭代判断标准,在ISPRS提供LiDAR点云数据集进行处理,结果表明改进的PTD方法能够把I类误差、Ⅱ误差和总误差分别减少10.26 %、0.79 %和8.07 %。为过滤机载LiDAR点云离散数据集滤波供了决方案。高广[6]等针对山区LiDAR点云的特点,对PTD滤波首先应用随机格网搜索算法获取更多精准初试种子点,其次利用地形预测角的判断准则提高地线附近滤波精度,最后在有大量地形断裂线山区进行验证,结果表明改进PTD滤波算法有效地保留了山区的地形断裂线特征。

对于禄丰恐龙谷南缘山区特殊地形地貌,采用低空无人机测量技术获取的影像数据,并构建DIM点云,首先考虑山顶、山体两侧和地面存在高程差,选择渐进加密三角网滤波算法进行处理。其次针对实验区特殊地理位置和地貌形态,改进传统渐进加密三角网滤波算法,最后对两种滤波算法处理结果进行直观分析。较为完整得到实验区地面点数据可为后续解译地貌特征信息和提取地形特征信息奠定基础。

2 PTD改进方法研究

2.1 DIM点云PTD滤波方法

目前无人机航测制作的主要数据之一是点云数据,以点集的形式表达被测物体的属性、位置信息和纹理信息。影像密集匹配获取点云,即DIM点云,不可透过地表的地物获取地面点,获取的地表点云与建筑物、信号塔和植被等地物构成地面表面模型(Digital Surface Model,DSM),存在噪声点及高程异常点。密集匹配点云滤波,在点云分类、三维建模、单体化提取等方面有重要作用。国内外研究者[7-14]在点云滤波方面开展了大量工作,依据主要滤波原理形成了基于坡度、基于形态学、基于曲面拟合、基于分割、基于机器学习、布料模拟滤波和渐进加密不规则三角网算法等主要滤波方法。基于不规则三角网的滤波算法,在丘陵地区、林区和山区等对点云数据有较好的适用性。经典的渐进加密不规则三角网(Progressive TIN Densification,PTD)滤波算法获取实验区地面点。PTD算法在商用软件Terrasolid的TerraScan模块中采用,对实验区滤波主要进行四步:①把实验区数据划分规则格网,认定格网最低点为地面点,在实验区选取初试种子点；②利用种子点构建不规则三角网(Triangulated Irregular Network,TIN)；③判断地面点,如图1所示,P点为非种子点,底面的三角平面顶点分别是V1、V2、V3,d表示非种子点到三角面的距离,α1、α2、α3是非种子点与底面三角形顶点连线的夹角,通过计算点P到三角形的距离d及α1、α2、α3的最大值分别与距离和角度的阈值比较,若小于阈值则定为地面点,归入三角网中。④迭代循环步骤②、③直至分类完所有实验区点云。

图1 判断地面点

2.2 PTD算法的改进技术

PTD滤波算法每次迭代过程中,都对其余各点到所在三角形的反复角和反复距离进行阈值判断,将满足条件的点加到TIN,初始构建的TIN对后续滤波判断影响较大。PTD滤波算法把尺寸大于研究区最大建筑物尺寸的格网中的最低点视为地面点,从该角度出发又考虑实验区山体两侧沟壑丛生,且对于山地滤波,山体与地面存在高程差,山体的地面点更易被误识别为非地面点。因此在PTD滤波算法选取种子点时,加入山脊线和山谷线交汇的地形特征点,通过地形特征点使山体两侧构建密集且符合实际地形的三角网,使山体两侧的点易通过角度和距离的阈值。针对实验区特殊地形分布情况,改进种子点选择方法,使点位分类在构建TIN中符合实际情况。如2图所示,C1～C5是初始种子点,A和B是在脊-谷的位置选择特征点,AB连线是构建三角网底边,点G1和G2是三角网下部的点,与B点在同一三角网但高程值低于A点。G2点为非地面点,G1为地面点,所以仍不能将三角网的下部点均视为地面点。

为此需要不断剔除最远点和增加点构建移动曲面,确保点云规则传递并覆盖整个实验区,如3图示,主要进行:①把实验区固定格网按建筑物最大尺寸划分；②确定实验区最外围格网高程最低点,以此为地面种子点；③点A、B、C、D是研究区四个角格网中心,以此构网,域内全部点为中心建立移动格网；④如图4示,进行多尺度迭代把移动格网内最低高程点作为地面种子点；⑤剔除重复点,并和格网外地面种子点建立实验区地面种子点。

图2 表示待定点识别

图3 移动格网构建

图4 格网多尺度迭代

为避免格网缺失产生的错误分类,格网ABCD外部不再构建移动格网。每个移动格网均依照实验区建筑物最大尺寸划分,既保证了所有格网均有地面点,又确保获取地面种子点无需进一步优化,直接进入三角网迭代判断的步骤。

3 试验区及试验数据采集

3.1 试验区

研究区域位于云南禄丰恐龙国家地质公园南缘山区,隶属云南楚雄彝族自治州,地理坐标为N 24°51′33″～25°30′45″,E 101°38′06″～102°24′34″。测区范围内地貌类型复杂多样,以构造侵蚀地貌、方山地貌为主,地势东北高,西南低,最高海拔2200 m,最低海拔1302 m,干湿分明。由于测区内存在小型中生代红色沉积盆地,成土母岩由中生代紫色砂页岩和元古代的碳酸盐岩交错分布,受母岩影响,土壤类型呈带分布,紫色和红色土壤分布广泛,测区位置如图5所示。

图5 实验区概况

3.2 试验数据采集

试验区影像数据获取是通过DJI Phantom 4 RTK进行影像数据采集。首先根据测区实际周边地理环境,交通状况,结合Google earth平台对测区进行航线规划,其次通过前期航线规划以及后期数据质量要求,选择合适的参数设置。鉴于本次数据获取主要面向微地貌特征分析应用,故本次飞行参数设置为:航向重叠度80 %,旁向重叠度80 %,平均飞行高度150 m。本次航测天气条件良好,共采集392张影像,影像平均分辨率为0.07 m。

4 实验分析

4.1 基于DIM点云可视化分析

通过对实测影像进行提取同名像点、特征匹配、计算旋转矩阵、对选择矩阵进行畸变校正和光束法平差迭代完成影像密集匹配获取DIM点云。点云密度为154.66个/m2,标准差为16.13 m。因无人机航测影像仅包含可见光波段,不具备透过地表已有地物如:植被、建筑物、信号塔等,即DIM点云含有噪声点和高程异常点,初步获得地表DIM点云生成DSM,如图6所示。

图6 数字表面模型(DSM)

由图6可知,恐龙谷南缘山区受亚热带低纬度高原季风气候影响山体沟壑交错,常年干湿分明,也导致山体右侧坡度通过目视解译判断接近90°。DSM中明显高程异常有T1、T2和T3的三个区域,在数字正射影像(Digital Orthophoto Map:DOM)图7中选定同样三个区域,DOM包含真实色彩信息,可判断出T1区域表示实测区域植被覆盖区域,主要是低矮独立树和灌木。T2区域是山顶的信号塔,信号塔在一定程度确保了免像控无人机对山区进行航测时获取POS的精度,但因信号塔自身高程值,后续构建DEM同样需先进行点云滤波处理。T3区域是山脚的建筑物和蔬菜大棚,经实地勘察该区域建筑物不超两层,除上述三个区域的地表地物,实验区西北角和东南角分布的低矮灌木和农作物均需要滤波处理。

图7 正射影像图(DOM)

4.2 基于PTD改进的点云滤波处理成效对比分析

4.2.1 PTD滤波处理及去噪成效分析

实验通过“试错法”调整阈值,选定最大地形坡度角为88°,迭代角度阈值设置为6°,迭代距离设置为1.4 m。70个初试种子点选取如图8(a)示,种子点除了边界点,其余均匀分布于实验区,由70个种子点构建TIN如图8(b)示,构建三角网较为规整。

(a) 滤波种子点

(b) 种子点构建不规则三角网

实验区采集DIM点云有17271784个点,由初试种子点构建TIN进行PTD滤波,剔除非地面点并进行冗余数据抽希,最后保留地面点5723274个点。

4.2.2 改进PTD滤波处理技术应用

对于上述PTD滤波算法对实验区造成山体部分区域地面点被误识别为非地面点,实验区位于滇中高原,山体顶部两侧与山底存在高程差,山顶和山脊区域易误选为非地面点。因此,在PTD滤波算法基础上选取种子点时,依据实验区特性,在ArcGIS中基于地表水流分析和几何分析相结合的方法提取山脊线和山谷线,二者交汇的脊-谷地形特征点为改进PTD滤波算法的种子点。图9为主要山脊线、山谷线分布和脊-谷交汇点示意。

图9 主要脊-谷线示意

对于地形特征点的选取精度评价需构建一个评价模型。评价模型的构建主要包含:①剔除冗余数据,保留主要脊-谷特征线；②统计地形特征点数量,用M表示；③主要脊-谷特征线以1 m、5 m、10 m和15 m为线性单位构建缓冲区,计算落入缓冲区范围内地形特征点数量,压盖边缘的点同样算入,用Ni(i=1,5,10,15)表示其数量；④Qi(i=1,5,10,15)表示在缓冲区内地形特征点数量占总点数的比,具体如式1,以此为地形特征点选取精度标准。

(1)

统计分析如表1所示,本次共提取486个地形特征点,以1 m为距离的缓冲区中涵盖27个点,仅占总数的5 %,5 m为距离的缓冲区包含132个点,占总数的27.16 %,15 m为距离的缓冲区有429个点,占88.27 %。后续试验较少的种子点会造成滤波精度低,较多的种子点在处理过程中导致处理效率不高,故以10 m为缓冲距离覆盖地形特征点为改进PTD滤波算法的种子点。

表1 地形特征点提取精度评价

以图8(a)初始种子点并叠加提取的脊-谷交汇的地形特征点,二者共有300个点为种子点,如图10所示。通过此类改进使选取种子点贴合实验区地形分布,并对选取种子点完成TIN构建如图10(b)所示。通过改进方法构建三角网,位于山体两侧相对于初始构建,三角网更为密集。改进滤波算法并进行冗余点抽希保留6692859个地面点,相对于PTD滤波算法保留地面点更充足。

图10 (a)地形特征点

图10 (b)特征点构建不规则三角网

4.2.3 改进PTD去噪结果对比分析

对两种滤波方法保留的地面点进行目视解译,对PTD滤波算法得到地面点结果和脊-谷交汇特征点叠加初始选取种子点结合构建不规则三角网点云滤波结果均显示真实色彩信息,以图11和图12示。对比分析两种点云滤波结果,结果表示:①两种方法均保留山体整体轮廓特征,上述可视化分析地表地物T1区域的低矮植被PTD滤波算法结果基本保留原有特征而改进PTD滤波算法剔除部分低矮植被,并未完全清除,植被覆盖区域点云出现小面积空洞现象；T2区域两种滤波算法均把信号塔剔除；T3区域建筑物被两种算法过滤,但PTD滤波算法对蔬菜大棚仍有冗余信息,而改进PTD相对来说剔除大部分蔬菜大棚。②实验区山体按道路划分左侧K1区域和右侧K2、K3区域PTD滤波算法均误把山体面点识别为非地面点剔除,从而在山体两侧出现点云空洞的情况。而改进PTD滤波算法对K1、K2和K3区域的山体地面信息都较为完整保留。

图11 PTD滤波算法结果

图12 改进PTD滤波算法结果

5 结 论

无人机DIM点云滤波处理不仅在地物分类、地物单体化提取和地形特征分析中起决定性作用,剔除地物点对后续生成相关产品也是核心步骤。对于恐龙谷南缘山区DIM点云,考虑实验区位于滇中高原和山体本身与地面就存在高程差,选择使用PTD滤波算法对山区进行滤波处理并依据地形特征点对PTD滤波算法进行改进试验研究。结果表明:①明显地物在被剔除的同时,较为完整保留整个实验区,明显地物中山体两侧低矮植被和山脚蔬菜大棚基本未被剔除,山顶信号塔和山脚建筑物全部清除,但山体部分的地面点易被识别为非地面点而在出现山体K1、K2、K3区域的空洞现象。②对此,针对恐龙谷南缘山区复杂地形,提出以脊-谷交汇地形特征点为改进PTD滤波算法的种子点,在山体两侧精细构网,山体低矮植被部分清除,相对于PTD滤波算法蔬菜大棚大面积被清除。且山体两侧未出现明显点云空洞的现象。

DIM点云通过脊-谷交汇地形特征点改进PTD算法对山区点云滤波相对于经典PTD算法更为精准区分地面点和非地面点,为后续地形地貌提取和分析提供了空间数据支持依据。但对于复杂地形和地形起伏度较大区域,对于地形特征点作为种子点的PTD滤波算后续研究中对提高提取地形特征点的精确性有待进一步加强。