2. 国土资源部城市土地资源监测与仿真重点实验室, 广东 深圳 518000;
3. 西南交通大学地球科学与环境工程学院, 四川 成都 611756;
4. 香港中文大学太空与地球信息科学研究所, 香港;
5. 首都师范大学成像技术高精尖创新中心, 北京 100037;
6. 国家基础地理信息中心, 北京 100083;
7. 四川省第二测绘地理信息工程院, 四川 成都 610100;
8. 中国测绘科学研究院, 北京 100083;
9. 黑龙江地理信息工程院, 黑龙江 哈尔滨 150086;
10. 首都师范大学资源环境与旅游学院, 北京 100037
2. Key Laboratory of Urban Land Resources Monitoring and Simulation, Ministry of Land and Resources, Shenzhen 518000, China;
3. Faculty of Geosciences and Environmental Engineering, Southwest Jiaotong University, Chengdu 611756, China;
4. Institute of Space and Earth Information Science, the Chinese University of Hongkong, Hongkong, China;
5. Beijing Advanced Innovation Center for Imaging Technology, Capital Normal University, Beijing 100037, China;
6. National Geomatics Center of China, Beijing 100830, China;
7. The Third Geoinformation Mapping Institute of National Administration of Surveying, Mapping and Geoinformation, Chengdu 610100, China;
8. Chinese Academy of Survey and Mapping, Beijing 100830;
9. Heilongjiang Institute of Geomatics Engineering, Haerbin 150086, China;
10. College of Environment and Tourism, Capital Normal University, Beijing 100037, China
“一带一路”倡议和“走出去”国家战略急需全球地理信息资源保障，由于全球地理信息资源的缺乏，我国急需开展“一带一路”重点区域地理信息资源建设[1]。为实现全球重点区域的高精度数字高程模型(DEM)覆盖，如何利用国产高分辨率卫星影像高效生产全球DEM已经成为我国全球地理信息资源建设工程的重大任务。
全球尺度的三维地理信息获取通常采用光学立体影像或雷达干涉测量的方式，例如我国的资源三号(ZY-3)三线阵卫星[2-3]与德国TanDEM-X雷达卫星[4]。欧洲空中客车防务及航天公司(Airbus Defence and Space)，利用TanDEM-X雷达卫星在2010年底至2014年采集的高分雷达生成DSM数据，通过对地形及水域等典型区域的人工交互式编辑，生产了覆盖全球的高精度WorldDEM产品[5](格网分辨率12 m，垂直精度为2 m(相对)/4 m(绝对))。我国国家测绘地理信息局于2015年启动全球地理信息资源建设工程，已利用资源三号卫星(ZY-3)影像数据生产了境内外DOM和DSM产品[6]。针对我国全球DEM生产任务，拟利用ZY-3影像的DSM产品，通过滤波编辑进行高精度DEM生产。然而，我国各生产单位采用的传统立体测图的方式进行DEM生产，一幅标准的1:5万DEM产品，耗时需1人/月，其中立体视图下采集高程点需1~3 d，剩余90%的工作量则主要耗费在植被、水域、城区等特殊地形地貌特征的交互式修复、修补等后处理工作。DEM生产的后处理操作量大，费时费力，且生产效率及产品精度完全依靠作业员的经验和熟练程度，效率和精度难以保证。技术上仍面临以下挑战：
(1) 传统点云滤波方法的瓶颈问题。传统点云滤波算法原理上需假设在局部邻域内的最低点为地面点，通过建立临时的地面结构，利用阈值参数对地面点进行加密，区分地面与非地面点。全球地形起伏变化、地表结构形态复杂(如图 1所示)，在点云密度和数据特性上均有较大差异。现有滤波算法通常对同一区域采用一套参数，难以适应错综复杂的地面-非地面地形结构，且对参数变化较为敏感，导致误分类问题严重[7-8]。当区域内地形地貌种类较多时，无法批量处理整块区域、滤波参数顾此失彼，需大量的交互式后处理操作。例如在处理图 1(c)所示的密集城区时，常出现大量房屋无法滤除或因过度滤波导致的山地地形特征缺失或河岸线消退等错误结果。自动、高效、可靠的点云滤波算法依旧是学术界尚未解决的难题，更是生产实践普遍面临的技术瓶颈[9-10]。
|
| 图 1 不同地形地表的DSM数据 Fig. 1 DSM data for different scenarios |
(2) 利用ZY-3卫星影像的DSM数据进行DEM生产的难题。光学卫星影像的DSM，存在明显噪声及缺漏现象。采用光学立体卫星影像，通过影像密集方式获取密集匹配点云，并进行DSM重建的过程，主要步骤包括：利用影像灰度信息，获取影像像素点的坐标对应关系，并通过空间前方交会方式获取地面/地表点的三维坐标。不同于激光雷达(LiDAR)点云，影像密集匹配点云受限于影像纹理信息，在无纹理特征、纹理特征较弱、纹理特征重复出现的区域，由于像素灰度信息难以可靠恢复影像对应关系，难以获取有效的同名点，导致匹配噪声多甚至匹配失败。在利用DSM数据进行DEM生产的过程中，此类噪声以及缺失的地面信息，均需要后续交互式的平滑、修补处理等后处理操作。例如，无纹理特征的水域(如图 2所示)，传统立体测图方式下，DSM通常需进行人工交互式编辑，在DSM自动滤波或交互式编辑中，因水域高程低于河岸，常会出现河岸扩张、湖中岛缺失现象。因此，构建顾及点云噪声影响及保地形特征的点云滤波方法，是目前亟待解决的问题。
|
| 图 2 水域处理效果 Fig. 2 Processing results for the water bodies |
不同于LiDAR点云，光学影像密集匹配点云无多次回波效应，无法穿透植被，因此在大范围植被覆盖区域，密集匹配点云中不包含地面点信息[11]。传统DEM生产过程，对于大范围林地覆盖的区域，多是采用人工交互方式，将DSM降低固定的植被高度值后再做整体平滑操作，可靠性无法保证[12]。因此对于ZY-3 DSM，在林地等不包含地面信息的区域，需引入全球土地利用信息，如GLC30、正射影像分类信息，利用植被覆盖信息，进行定向精细滤波，降低植被高度。并且，由于林地通常出现于山区，在定向滤波的同时，需同时保留山地的典型地形特征，如山脊、山谷等特征线。
最后，DSM在高楼密集区域会因为遮挡和视差断裂原因导致地面点缺失问题，导致在建筑物边缘出现缺漏，使建筑物的直角特征缺失、变平滑。这类地形区域的自动滤波结果通常呈现为某些或成片建筑无法过滤完整，或建筑物区域存在显著高于周围地表、却低于建筑物高度的离散格网点。因此针对上述现象，需要针对性的能获取平滑地表且可降低建筑物高度的后处理交互式编辑方法。
综上所述，我国拟用ZY-3卫星影像的DSM产品进行全球DEM生产。ZY-3影像的DSM产品具有如下特点：①在低纹理、水面等区域存在明显噪声、缺漏现象；②在成片林地或其他大面积植被覆盖区域，不包含地面信息，因此违背了滤波算法的局部邻域内最低高程为地面高程的基本假设，需要引入全球土地利用信息进行定向滤波；③在建筑物覆盖密集的城市区域，会由于遮挡、视差断裂等因素，导致在建筑物边缘出现缺漏，使建筑物的直角特征缺失、变平滑，该特性违背了点云滤波的地面平滑的基本假设，需要针对性的后处理交互式编辑。针对上述问题，笔者提出了一种面向全球DEM生产的点云智能滤波与泊松编辑方法，①顾及弯曲能量的点云自适应滤波，以弯曲能量系统显式定量刻画复杂的地形地表结构特征，并以此驱动滤波模型参数的自适应优选；②构建多边界约束的地形泊松编辑框架，引入分类、地形地貌特征线等边界信息，支撑并实现城区、林地、山地、水域等复杂地形的交互式保特征泊松编辑。
1 面向全球DEM生产的点云智能滤波与泊松编辑原理 1.1 顾及弯曲能量的点云自适应滤波方法顾及弯曲能量的点云自适应滤波方法(图 3)，主要是针对DSM点云数据，在局部邻域的最低点为地面点与地表平滑这两个前提下，采用金字塔的滤波策略，逐级加密非地面点，从而逼近真实地表。在金字塔滤波的过程中，通过带有抗噪性顾及弯曲能量正则化约束的DEM内插算法与自适应的参数优选算法，提高自动滤波算法对噪声与参数选择的稳健性。
|
| 图 3 顾及弯曲能量的点云自适应滤波算法 Fig. 3 Adaptive surface filter by means of bending energies |
1.1.1 渐进金字塔滤波策略
基于地面最低高程假设，在构建点云金字塔时，通常采用局部邻域内的最低点作为金字塔对应分辨率的格网数值，因此在从粗到精的金字塔滤波策略中，每层的点云属于地面点的可能性将逐渐变小。若采用传统的四叉树的金字塔结构，如图 4(a)所示，每层级之间的点云数量将为1:4，待分类的点过多，易导致地面点误分，而错误点将在逐级滤波过程中被不断放大。因此，本方法设计提出了渐进金字塔的构建方法，即金字塔的每层级之间分辨率渐进加密，增加层级数目，提高分类可靠性[13]，如图 4(b)所示。
|
| 图 4 渐进金字塔构建 Fig. 4 Construction of the progressive pyramid |
1.1.2 顾及弯曲能量正则化约束的局部抗噪性内插方法
传统的DEM内插方法，通常是通过拟合地面控制点，P={pi=(xi, yi, zi)|i=1, 2, 3, …, n}，得到一张参数化的曲面，z=f(x, y)，在拟合过程中，由于地面控制点中不包含粗差，因此仅需要保持对地面控制点即可有较好的拟合度，即数据拟合度εdata(data term)具有较小的值，或完全经过所有控制点，如式(1)所示
(1)
在实际点云滤波过程中，地面点不可避免地会引入非地面粗差点，因此此时即使εdata为0，也并不意味着对地面有较好的拟合程度。现有的研究已经证明，采用地表一致性约束，可以一定程度上，降低噪声敏感性，提高内插算法的抗差性。现有的聚类、平面拟合的一致性约束方法并不通用，且计算耗时。弯曲能量εsmooth见式(2)
(2)
式(2)可以表达一张参数化曲面的二阶连续程度，因此若将εsmooth作为一个正则化约束引入地面内插过程中，则可隐式地蕴含着地表一致性约束[14]，抵抗噪声点的干扰。即通过优化如下能量函数，进行DEM内插，如式(3)所示
(3)
该能量函数可以很方便地采用带正则化约束的薄板样条函数来拟合求解[15]。
1.1.3 顾及弯曲能量的滤波参数自适应优选方法在DEM内插过程中，每个格网点对应于一个参数化的薄板样条函数曲面，因此也可以计算其弯曲能量值。考虑到在山脊、地面断裂处，在由粗到精的点云滤波过程中，内插得到的地表通常会低于真实地表值，如图 5所示。一方面，由于这些区域的地形结构弯曲起伏更为剧烈，通常表现为较大的弯曲能量，因此在对这些区域的点云进行滤波处理时，依据弯曲能量的大小补偿滤波参数阈值，以此显式约束滤波参数阈值自适应调整。另一方面，在金字塔滤波策略中，由于窗口大小不同，因此为了顾及坡度变化的影响，也需要考虑到尺度变化的影响，滤波阈值应随尺度的增大而增加。因此最终的滤波阈值Zt，如式(4)所示
|
| 图 5 顾及弯曲能量的滤波参数自适应优选方法 Fig. 5 Automatic parameter tuning by bending energies |
(4)
式中，t为人为指定的初始阈值；s为尺度信息；b为弯曲能量信息。通过以上信息实现滤波阈值参数随地形变化、尺度的自适应优选。
本文方法是在笔者的已有研究[13]基础上进行了扩展。前期研究主要是针对Lidar点云滤波瓶颈，提出了顾及弯曲能量的Lidar点云自适应滤波方法，该方法利用ISPRS提供的17种标准数据集进行了多组试验分析，并设计了误差评价体系，定量分析了弯曲能量方法的可靠性和自适应能力。
1.2 多边界约束的泊松编辑方法 1.2.1 泊松编辑方法的基本原理DSM数据中可能由于遮挡、大范围植被、水域等因素导致自动滤波算法效率降低，交互式的后处理DEM编辑必不可少。在DEM编辑中需要同时顾及DEM的平滑特性与水域、山脊特征线等特征信息，传统简单的DEM编辑方法并不能满足上述需求，因此本方法引入一种泊松地形编辑(Poisson terrain editing)的策略[16]。泊松编辑的含义是，通过优化域中的一个已知导向信息v，并给定边界条件f|∂Ω，通过下列积分方程，求解未知函数f。定义影像范围
(5)
对于DEM编辑而言，优化域即DEM上的一个可带洞的闭合多边形区域，在区域内，其导向域条件是依不同滤波算法给定的，该优化域必须闭合，即被边界条件完全包含。由于优化域可为带洞多边形，因此在优化域内部也可以包含任意边界条件，如图 6所示。一方面，由于边界条件在优化过程中保持不变，因此可以用于保留DEM中的地形结构特征信息，如山脊线、山谷线、河岸线、河流区域等。另一方面，通过调整优化域Ω中的导向信息的定义方式，实现不同类型的滤波效果，如平滑、恢复梯度信息等，如图 7所示。
|
| 图 6 泊松地形编辑的导向域与边界条件 Fig. 6 The guided field and boundary conditions for Poisson terrain editing |
|
| 图 7 DEM编辑中的边界条件 Fig. 7 Boundary conditions in DEM editing process |
1.2.2 泊松地形编辑方法
本方法中，泊松地形编辑方法包含3个主要步骤：①依据滤波需求创建多边形区域，定义优化域与包围优化域的边界条件；②在优化域中定义合适的导向信息以实现不同的滤波效果，并给定合适的边界条件的初值；③求解方程，获取泊松方程的解，即滤波后的DEM。
1.2.2.1 交互式编辑创建多边形如图 6和图 7所示，在二维影像域S中，待编辑的优化域Ω及其边界区域∂Ω，都可用多边形进行表达，其中∂Ω可为Ω的边界线或Ω内部的空洞，此在本方法中采用通用的带洞多边形来支持优化域与边界条件的创建与编辑。本方法主要支持两种方式的多边形创建方法：①交互式勾绘创建多边形。②采用区域增长的方式创建多边形区域。对于前者，可以采用折线或流形的方式操作，其中折线为每次点击创建一个点，最后闭合多边形，流形通过拖动鼠标绘制平滑曲线并闭合。而区域增长可根据DSM/DEM高程值或分类信息进行区域增长，并通过等值线提取算法获取闭合多边形区域。
1.2.2.2 滤波算法创建在创建了优化域Ω与边界区域∂Ω后，泊松方程的设计矩阵和结构实质上已经可以确定，但为了获取最终优化域Ω中的未知函数f，仍然需要创建优化域中的导向信息v与边界区域的初值f*|∂Ω，此信息决定了不同的滤波效果。目前本方法已集成实现了十几种不同的滤波效果，在此列举几个说明导向信息与初值的创建方式。
(1) 复制DSM滤波：此滤波为恢复原始DSM的形状，用于在某些滤波过程中由于滤波选择不恰当导致编辑DEM失真，此滤波可以保证在DEM边界高程不变的情况下，恢复DSM的地形起伏特征。值得注意的是，此算法并非直接恢复DSM的高程值，因为这样会导致编辑的边缘产生突变。对此算法，导向域是直接从DSM中采用拉普拉斯算子计算所得的，而初值为DEM边界区域中的高程值。
(2) 林地高程降低滤波：在林地，可通过分类信息创建优化域与边界区域，并给定一个待降低的高程值，整体降低林地高度。然而若直接整体降低高程，势必会在边缘区域产生突变效应，因此对于此滤波，导向域为当前DEM的拉普拉斯算子计算所得，即保证当前DEM的起伏特性不变，边界区域的初值，在编辑区域的边缘采用DEM的高程值，而在内部随机创建一系列点状的边界区域，其初值为DEM减去降低高程值。采用此方法，可保证编辑后植被高度正常降低且不会产生显著突变效果。
(3) 均值/中值/最低值平滑：同样，若简单地采用高斯平滑，难免会在编辑的边缘处产生突变效果，因此本软件采用泊松编辑的方式实现平滑，其核心思路是采用0值作为导向域信息，而采用均值/中值/最低值作为内部边界区域的初值。该方法等价于薄膜样条函数内插(membrane interpolation)[17]。
1.2.2.3 泊松方程求解在构建了泊松方程后，在离散的二维影像域中可将式(5)确定的初值积分方程，转化为一个大规模的稀疏方程AX=B，其中A为可逆方阵，而X为最终优化域中每个格网点的值，该方程可采用稀疏矩阵的LU分解[18]或共轭梯度优化方法进行求解(conjugate gradient)[16, 19]。由于传统基于CPU的优化算法时间较长，难以满足交互式的实时需求，因此本文基于通用的并行GPU加速计算框架CUDA与OpenCL实现了共轭梯度法的并行算法，速率较传统CPU算法提高10倍，可满足交互式的实时需求。
2 试验分析基于本文所述两项关键技术，自主研制了点云智能滤波与泊松编辑软件LINK。在试验分析部分，本文将着重介绍利用不同类型不同地区的全球DEM生产性试验情况，并将试验结果与当前我国其他两款主流软件GEOWAY CIPS和PixelGrid进行对比分析。
2.1 数据试验数据包括全球地理信息资源建设工程前期生产的DSM、DOM成果，均采用WGS-84坐标系，以1:5万标准分幅为存储单元。基于DSM生产DEM的难点主要在于大面积植被覆盖区域和建筑物覆盖区域，因此试验数据涵盖了建筑区、森林、水域等典型数据具体见表 1。
| 国家 | 图号 | 地形类别 | 地物要素特征 |
| 缅甸 | NC47E006011 | 平地 | 植被、水域、裸地等 |
| 缅甸 | NE47E019001 | 平地 | 建筑、植被、水域、裸地等 |
| 缅甸 | NE47E017002 | 丘陵地 | 建筑、植被、水域、裸地等, 其中大面积植被覆盖区域 |
| 阿富汗 | NH41E005013 | 丘陵地 | 存在大面积裸露地表 |
| 缅甸 | NE47E015008 | 山地 | 存在大面积植被、水域区域 |
| 阿富汗 | NJ42E020023 | 高山地 | 建筑、水域、裸地等 |
试验数据分布情况如图 8所示。
|
| 图 8 试验数据分布 Fig. 8 The distribution of test data |
2.2 试验结果质量定量评定与分析
本次试验进行了多种检查方法对DEM成果数据进行全面地比较、质量评定与分析，其中包括：①采用人工采集的检查点进行质量评定；②利用原始DSM与DEM成果数据的差值检查滤波可靠性；③在立体环境下检查DEM成果与空三加密成果符合程度；④以传统的立体采集地貌方式生产的DEM数据为参考检查滤波成果。以下重点介绍采用人工采集检查点进行质量评定与分析的结果。
(1) 采用人工采集的检查点进行质量评定。人工采集的检查点是在DOM上选取位于建筑物顶部、植被顶部等区域的点位，其X、Y值从DOM获取，Z值在DSM的高程值基础上，立体环境下读取建筑物、植被的高度，对DSM高程值进行修正，从而得到最终的Z值。在建筑和植被区域，均匀采集了大量检查点；在裸地、水域等区域，主要采集在能反映地形特征位置的检查点，具体如图 9所示。
|
| 图 9 检查点分布 Fig. 9 Check points distributions |
(2) 精度检查。利用检查点数据对3种软件的滤波结果进行精度检查，结果见表 2。
| 试验数据 | 检查点个数 | 中误差/m | |||
| LINK | PixelGrid | GEOWAY CIPS | |||
| 建筑 | NE47E019001 | 364 | 3.430 | 3.539 | 3.509 |
| NE47E017002 | 88 | 2.621 | 2.797 | 2.787 | |
| NJ42E020023 | 43 | 4.007 | 4.182 | 4.310 | |
| 植被 | NE47E006011 | 98 | 1.573 | 3.256 | 4.306 |
| NE47E017002 | 202 | 1.865 | 2.271 | 2.139 | |
| NE47E015008 | 279 | 2.646 | 16.124 | 6.726 | |
| 裸地 | NH41E005013 | 423 | 1.618 | 2.133 | 2.628 |
从表 2可以看出，针对建筑、植被及裸地区域，本方法的滤波精度明显优于其他两款软件，对于植被区域的滤波，泊松滤波方法优势更突出。
2.3 与传统人工立体采集生产的DEM对比为了将DSM滤波生产的DEM与传统方式生产的DEM进行效果比较，本次试验采用立体采集地貌的方式生产了6幅试验数据的DEM成果，并对两种生产方式的DEM效果、精度进行了详细的对比分析，如图 10所示。
|
| 图 10 成果DEM对比 Fig. 10 Comparison between different DEM results |
从表 3可以看出，利用DSM滤波生产DEM的方式相比传统生产方式而言，精度明显提高。证明本文方法在进行DEM生产时，避免了人工采集误差、格网点内插误差等影响因素，可以尽量保持实际地貌细节。本方法相比传统生产方式而言，生产效率大大提高，约提高至少2~3倍。尤其是在山地、高山地等DEM生产困难地区，本方法的优势更加明显。
| 图幅 | 地形 | 检查点个数 | 中误差/m | |
| 传统方式DEM成果 | LINK生产DEM成果 | |||
| NE47E006011 | 平地 | 98 | 1.871 | 1.573 |
| NE47E019001 | 平地 | 364 | 4.574 | 3.430 |
| NE47E017002 | 丘陵地 | 290 | 5.105 | 2.652 |
| NH41E005013 | 丘陵地 | 423 | 4.265 | 1.618 |
| NE47E015008 | 山地 | 258 | 6.984 | 2.646 |
| NJ42E020023 | 高山地 | 43 | 5.141 | 4.007 |
2.4 与其他主流软件的试验结果对比分析
在对比分析LINK与当前我国其他两款主流软件GEOWAY CIPS和PixelGrid的滤波结果时，将重点分析建筑、植被和水域的处理效果：建筑、植被是否有效滤掉，是否将非地面点有效降至地面；在水域等地貌突变区域，DEM数据是否准确保留了地貌细节、地形结构特征，比如湖中小岛，河流边缘纹理，细窄的沟壑，纹理变化复杂的山地区域等。
2.4.1 建筑区滤波效果对比建筑区滤波效果对比，如图 11所示，可见其他两款软件均存在大量建筑物未滤除，LINK对建筑物区域的滤波效果较好。LINK软件设计了建筑物区域的泊松滤波模块，可根据自动分类提取出的建筑区域多边形，实现建筑区域的局部自适应滤波；根据地形特征可将建筑区分为平地区域、丘陵区域及山地高山地区域，通过调整邻域大小及坡度阈值达到不同地形建筑区的最佳滤波效果。
|
| 图 11 建筑物滤波效果 Fig. 11 Filtering results for building areas |
2.4.2 植被区域滤波效果对比
由于DSM成果的格网间距为10 m，地物地貌细节体现不足，常规统一滤波算法难以将大面积植被区域对应的非地面点降至地面。LINK针对性设计了面向植被区域的泊松滤波模块，即根据自动分类提取出的植被区域多边形及其植被高度属性值，结合林地区域所在的地形特征，通过调整平滑大小和降低高程两个参数值，在保持地形特征的同时将植被区域非地面点高程降至地面。在大面积植被覆盖区域，通过滤波处理后，滤波结果的高程值普遍小于原始DSM高程值。因此试验将原始DSM数据减去3种软件的滤波结果，对栅格相减的值进行比较，效果如下(白色为DSM-DEM≤0，黑色为DSM-DEM≥0)，如图 12所示。
|
| 图 12 大面积植被覆盖区域DSM-DEM结果 Fig. 12 Filtering results for vegetation areas |
通过对比3种软件DSM与DEM差值栅格图可以看出，LINK软件可以根据圈定的植被范围，有针对性地将植被区域的非地面点降至地面，而CIPS和PixelGrid采用的整体滤波方法，难以识别大片植被覆盖区域，滤波效果不理想。
2.4.3 水域滤波效果对比对于水域区域，另两款软件均无特定的处理模块，滤波结果均需要进行大量的交互式后处理操作。LINK软件针对性的设计了水域泊松滤波模块。DSM成果数据的湖泊等静止水域范围是整体置平的，河流等流动水域是分段置平的，针对这一特点，LINK提供了水域自动提取功能，基于自动提取的水域范围设置缓冲区域，并将水域范围在后续处理中会作为泊松滤波的控制多边形存在，从而使河岸线等水域特征不受周围数据编辑处理的影响。
(1) 湖心岛滤波效果对比。湖心岛滤波结果如图 13所示。
|
| 图 13 湖心岛滤波结果 Fig. 13 Filtering results for island |
(2) 河岸线滤波效果对比。用原始DSM数据与3种软件的滤波结果进行栅格减运算，通过河岸线区域的变化值来说明3种软件的滤波方法对河岸线的保持效果。由图 14可见，LINK能够完整保持河岸线特征，而其他两款软件会出现河岸线消退的情况。
|
| 图 14 河岸线特征保持情况(DSM-DEM) Fig. 14 Filtering results for river bank |
3 结论
针对全球地形地表结构具有多样性和复杂性特点，现有依靠单一滤波模型或有限滤波规则的点云滤波方法均难以根据复杂多样的地表结构特征进行自适应调整，滤波结果的准确性和可靠性问题十分突出，导致大量人工交互后处理工作量，DEM的质量和精度难以保证。针对上述问题，本文提出了一种智能滤波与泊松编辑方法：以弯曲能量定量显式刻画地表形状特征，以此驱动滤波模型的自适应调整，提高滤波结果的可靠性；设计提出了多边界约束的泊松地形编辑框架，引入分类、地形地貌特征线等边界信息，实现建筑区、森林和水域等区域的定向智能精准编辑。经过大量不同类型不同地区的全球DEM生产性试验，证明了本文方法的可靠性和有效性，为全球DEM大规模生产提供了有力的技术方法支撑。
| [1] |
王东华, 刘建军.
国家基础地理信息数据库动态更新总体技术[J]. 测绘学报, 2015, 44(7): 822–825.
WANG Donghua, LIU Jianjun. Key techniques for dynamic updating of national fundamental geographic information database[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(7): 822–825. DOI:10.11947/j.AGCS.2015.20150089 |
| [2] |
李德仁.
我国第一颗民用三线阵立体测图卫星——资源三号测绘卫星[J]. 测绘学报, 2012, 41(3): 317–322.
LI Deren. China's first civilian three-line-array stereo mapping satellite:ZY-3[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(3): 317–322. |
| [3] |
唐新明, 张过, 祝小勇, 等.
资源三号测绘卫星三线阵成像几何模型构建与精度初步验证[J]. 测绘学报, 2012, 41(2): 191–198.
TANG Xinming, ZHANG Guo, ZHU Xiaoyong, et al. Triple linear-array imaging geometry model of Ziyuan-3 surveying satellite and its validation[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(2): 191–198. |
| [4] | KRIEGER G, MOREIRA A, FIEDLER H, et al. TanDEM-X:a satellite formation for high-resolution SAR interferometry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45(11): 3317–3341. DOI:10.1109/TGRS.2007.900693 |
| [5] | RIEGLER G, HENNIG S D, WEBER M. WorldDEM-A novel global foundation layer[J]. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 2015, 40(3): 183–187. |
| [6] |
王密, 杨博, 李德仁, 等.
资源三号全国无控制整体区域网平差关键技术及应用[J]. 武汉大学学报(信息科学版), 2017, 42(4): 427–433.
WANG Mi, YANG Bo, LI Deren, et al. Technologies and applications of block adjustment without control for ZY-3 images covering China[J]. Geomatics and Information Science of Wuhan University, 2017, 42(4): 427–433. |
| [7] | VOSSELMAN G. Automated planimetric quality control in high accuracy airborne laser scanning surveys[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2012(74): 90–100. |
| [8] | XU S, VOSSELMAN G, ELBERINK S O. Multiple-entity based classification of airborne laser scanning data in urban areas[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2014(88): 1–15. |
| [9] | PINGEL T J, CLARKE K C, MCBRIDE W A. An improved simple morphological filter for the terrain classification of airborne LiDAR data[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2013(77): 21–30. |
| [10] | SITHOLE G, VOSSELMAN G. Experimental comparison of filter algorithms for bare-earth extraction from airborne laser scanning point clouds[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2004, 59(1-2): 85–101. DOI:10.1016/j.isprsjprs.2004.05.004 |
| [11] | WANG Cheng, GLENN N F. Integrating LiDAR intensity and elevation data for terrain characterization in a forested area[J]. IEEE Geoscience and Remote Sensing Letters, 2009, 6(3): 463–466. DOI:10.1109/LGRS.2009.2016986 |
| [12] | ZHU Qing, LI Yuan, HU Han, et al. Robust point cloud classification based on multi-level semantic relationships for urban scenes[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2017(129): 86–102. |
| [13] | HU Han, DING Yulin, ZHU Qing, et al. An adaptive surface filter for airborne laser scanning point clouds by means of regularization and bending energy[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2014(92): 98–111. |
| [14] | SZELISKI R. Computer vision:algorithms and applications[M]. London: Springer, 2011: 812. |
| [15] | BOOKSTEIN F L. Principal warps:Thin-plate splines and the decomposition of deformation[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1989, 11(6): 567–585. DOI:10.1109/34.24792 |
| [16] | MORSE B S, YOO T S, RHEINGANS P, et al. Interpolating implicit surfaces from scattered surface data using compactly supported radial basis functions[C]//Proceedings of the International Conference on Shape Modeling & Applications. Washington, DC: IEEE Computer Society, 2001. https://erie.nlm.nih.gov/~dave/clsrbf.pdf |
| [17] | PÉREZ P, GANGNET M, BLAKE A. Poisson image editing[C]//SIGGRAPH'03 ACM SIGGRAPH 2003 Papers. San Diego, California: ACM, 2003. |
| [18] | FARBMAN Z, FATTAL R, LISCHINSKI D. Convolution pyramids[J]. ACM Transactions on Graphics (TOG), 2011, 30(6): 175. |
| [19] | WAECHTER M, MOEHRLE N, GOESELE M. Let there be color! large-scale texturing of 3D reconstructions[M]//FLEET D, PAJDLA T, SCHIELE B, et al. Computer Vision-ECCV 2014. Cham: Springer, 2014: 836-850. |
| [20] | KAZHDAN M, BOLITHO M, HOPPE H. Poisson surface reconstruction[C]//Symposium on Geometry Processing.[S.l.]: The Eurographics Association, 2006. http://hhoppe.com/poissonrecon.pdf |