CN102122395B - 一种保持地形特征的自适应尺度dem建模方法 - Google Patents

Abstract

本发明提出了一种保持地形特征的自适应尺度DEM建模方法,即构建一个适用于从广域到局部范围通用的精细DEM数据库，提出自适应尺度数据结构和分析方法，从精细尺度DEM数据库中根据应用需要提取多尺度的地形参数和特征，依据尺度信息保持不同级别地形特征，使得不同尺度下地形骨架基本一致。本发明将自适应尺度技术应用于数字地形分析领域，提高了多尺度DEM生成的效率，保证了多尺度下地形特征的一致性，满足了实际应用中的尺度需求。

Description

一种保持地形特征的自适应尺度DEM建模方法

技术领域

本发明属于数字地形分析领域，特别涉及一种保持地形特征的自适应尺度DEM建模方法。

背景技术

数字高程模型(Digital Elevation Model，DEM)用以表达自然环境中的地形、地貌特征，是空间数据基础设施中的重要组成部分。随着全球范围SRTM(航天飞机雷达地形测绘任务) 90米分辨率DEM数据库和ASTER GDEM(先进星载热发射和反射辐射仪全球数字高程模型)30米分辨率DEM数据库的公开，全国范围1:100万、1:25万、1:5万DEM数据库的建立及区域级1：1万DEM建库工作的相继完成，DEM数据库呈现出多尺度、多分辨率、多数据源的特征。

空间数据处理与分析的尺度效应是地球科学的一项基本问题，被认为是“所有研究的基础”(Quattrochi and Goodchild, 1997，参见背景文献1)。在数字地形模型与分析中，地形参数及特征具有多尺度特性，基于DEM的应用（如地质灾害、环境保护、水土保持等）也具有多尺度需求。当DEM的空间尺度变化时，提取的地形参数如坡度﹑坡向、曲率和地形特征等会出现相应的变化。例如在实际应用中，当研究大流域时，一个平缓的小山丘可以被忽略掉，而若考虑到区域的局部水系特征时，这个小山丘就会显得非常重要。

在纸质地图中，一幅地图数据是单一尺度的。为了解决多尺度问题，出现了1：100万、1：25万、1：5万等不同尺度的纸质地图。但这仅仅是多级固定尺度数据的组合，每一级只考虑了自己的尺度。纸质地图的局限性限制了实际应用。随着国家空间数据基础设施建设的不断完善，国家基础地理信息中各级DEM数据库相继建立。然而，这些DEM数据库仍是针对不同尺度的DEM构建不同的数据库。所谓的多尺度数据库是由多个固定尺度DEM数据库组合而成，从而导致数据冗余、数据更新维护不一致等一系列问题。可以说这只是纸质地图在数字化时代的延续，并没有找到多尺度数字地形表达之间的精髓和内在关系，没有从根本上生成一种适应于多尺度的DEM数据库。

现代DEM的管理和应用需要不同尺度间的协调并行，从固定尺度的应用逐渐演变为任意尺度的应用，进而转化为满足用户需求的自适应尺度应用。例如，对大区域滑坡泥石流灾害的监测往往建立在粗尺度DEM基础上，而真正希望做到滑坡泥石流灾害的预警，则必须能控制在非常精细的尺度上才能达到有的放矢。舟曲的泥石流灾害不仅希望在地方一级做到监测，同时要在国家一级就能进行预警。这就要求根据不同应用，动态构建自适应尺度的实时模型。该模型既涵盖非常大的区域，又可以进行局部精细研究，对于每一个细节都有所控制。这是现有的多级DEM数据库无法解决的问题，受限于现有的DEM分析尺度、建模方法和数据质量等因素。其中，最根本的问题还是现有的多级DEM数据模型所产生的弊病。因此，在实际中迫切需要适应于现代多尺度DEM应用的理论模型。

利用精细尺度的DEM综合得到粗尺度DEM是一种解决途径。目前通常采用的DEM综合方法是基于栅格DEM的重采样方法。然而，重采样方法具有“平滑”作用，地形参数随着尺度的变化会发生较大变化，无法保持地形的基本特征。Zhang and Montgomery (1994，参见背景文献2)采用多种分辨率(2–90m)检测了栅格尺寸对地形表达和水文模拟的影响，他们发现DEM分辨率会极大地影响地形参数和水位曲线的计算。Florinsky and Kuryakova (2000，参见背景文献3)利用微观尺度的地形属性（比如坡度、坡向、水平曲率、垂直曲率）相关系数的变化以及地表性质（比如土壤湿度）进行了DEM栅格大小的灵敏度分析，并提出了满足地形研究最小需求的DEM分辨率。Zhou and Liu (2004，参见背景文献4) 研究了DEM分辨率对多种坡度和坡向算法的影响，并推断许多DEM提取因子的不确定性与栅格数据结构紧密相关。研究证明由重采样方法得到粗分辨率DEM时，提取的地形因子被严重低估了。例如，当DEM栅格大小从30'' 增加到 32'，研究区域的平均坡度百分比就从12.4%减少到0.9%，很难直接从粗分辨率DEM计算坡度因子(Zhang et al., 1999，参见背景文献5)。在实际应用中，一个区域的平均坡度不应随着所选用的DEM数据尺度的变化而发生较大变化。尺度变化会损失细节，但会保留重要的地形特征，也就是“骨架”,需要根据尺度而保持不同级别的地形特征。

由于栅格DEM在多尺度下进行重采样后，给提取出的地形参数带来了很大的不确定性，另一种方法是采用不规则三角网（TIN)来表达地形面。Danovaro et al. (2006，参见背景文献6)提出了一种多分辨率表面网络来表达多分辨率地形模型，通过将地形面分解成以特征点（如极小值、极大值和鞍点）和积分线为边界的区域这种几何化简方法，用TIN描述了所有分辨率下的地形。但是寻找特征点和积分线并不容易，特别是对于复杂的地形面。Wu and Amaratunga (2003，参见背景文献7) 提出了小波不规则三角网（WTIN）的方法，该方法基于第二代小波理论，可对数字高程数据进行多分辨率表达。Li (2008，参见背景文献8)描述了一种基于特征点的办法来表达多尺度地形面，它形成了一种建立在多级地形特征点（比如山峰或土坎）上的多尺度数字地面模型结构，这样即使只有这些特征点保留下来，地形面也能得到大部分的保留。许多算法被用于探测沿地形剖面的特征点，比如Douglas-Peucker算法、VIP算法、最大Z-tolerance算法、启发性方法和分层算法。Zhou and Chen (2010，参见背景文献9) 提出了一种由特征点和流域约束线构建TIN的方法，这种TIN受到流域的约束，同时分析了从TIN提取出的地形参数，其结果表明在计算粗尺度下的地形参数方面，TIN结构能够比重采样得到更好的结果。

现有的发明中，自适应尺度技术大多应用于图像处理，与之相关的有：

国家海洋局第二海洋研究所的陈建裕(2007，参见背景文献10)提出了一种尺度自适应的图像分割方法，其发明以图斑自身属性界定图像对象的最优分割尺度，使得不同尺度的对象具有各自适宜的分割尺度。包括如下步骤：

１）根据图像及其变换结果的一个或多个图层确定待分割图像集，并以连续变化的尺度系数分割图像；

２）定义图斑显著性，以连续变化尺度分割下的图斑组成图斑演进曲线，并形成分割演进过程中图斑的显著性曲线、分割初尺度曲线；

３）从显著性曲线中计算极值形成尺度倒序，再根据尺度倒序中对应的极值形成极值尺度图像，最后以极值尺度图像中的优势尺度确定分割图斑。 

杭州电子科技大学的李轶等(2007，参见背景文献11)提出了一种多尺度自适应对比度变换的医学图像增强方法，通过将医学图像分解成金字塔形状排列的分辨率逐步降低的图像集合后，对分解得出的分层系数做调整，最后将调整系数后的各个细节层次图像再合成原图像增强后的图像，使得任意大小的低对比度区域可见度都可以得到较大的提高，且不会产生谷粒效应。

Hughes Aircraft Company的Frankot and Robert(1996，参见背景文献12)提出了一种干涉SAR数据的自适应尺度方法，当干涉SAR数据中高信噪比区域的分辨率需要保留时，就利用自适应尺度滤波对低信噪比区域进行平滑处理。自适应尺度滤波将多种应用于复杂干涉图的线性滤波进行组合，不仅计算效率高，还适合并行实现，同时采用了一种多级的金字塔结构，减少了计算负载和内存需求。

以上发明只适用于图像处理与分析领域，在本发明涉及的数字地形分析领域中，由于多尺度DEM数据库带来的数据冗余、不一致性以及无法满足实际应用的尺度需求，需要解决这些问题，同时要保持不同尺度下地形的基本特征。

背景文献：

1.      Quattrochi, D.A., Goodchild, M.F. (eds.), 1997. Scale in Remote Sensing and GIS, CRC Press。

2.      Zhang, W., Montgomery, D.R., 1994. Digital elevation model grid size, landscape representation, and hydrologic simulations. Water Resources Research 30(4), 1019-1028；

3.      Florinsky, I.V., Kuryakova, G.A., 2000. Determination of grid size for digital terrain modelling in landscape investigation–exemplified by soil moisture distribution at a micro-scale. International Journal of Geographical Information Science 14(8), 815-832；

4.      Zhou, Q., Liu, X., 2004. Error analysis on grid-based slope and aspect algorithms. Photogrammetric Engineering and Remote Sensing 70(8), 957-962；

5.      Zhang, X., Drake, N.A., Wainwright, J., Mulligan, M., 1999. Comparison of slope estimates from low resolution DEMs: scaling issues and a fractal method for their solution. Earth Surface Processes and Landforms 24(9), 763-779。

6.      Danovaro, E., Floriani, L.D., Papaleo, L., Vitali, M., 2006. A multi-resolution representation for terrain morphology. Lecture Notes in Computer Science 4197, 33-46；

7.      Wu, J., Amaratunga, K., 2003. Wavelet triangulated irregular networks. International Journal of Geographical Information Science 17(3), 273-289；

8.      Li, Z., 2008. Multi-scale digital terrain modelling and analysis.Advances in Digital Terrain Analysis. 59-83；

9.      Zhou, Q., Chen, Y., 2011. Generalization of DEM for terrain analysis using a compound method. ISPRS Journal of Photogrammetry and Remote Sensing.Vol 66,1. 38-45。

10.  国家海洋局第二海洋研究所，陈建裕，一种尺度自适应的图像分割方法，中国专利, 200710068214 , 2007-09-19；

11.  杭州电子科技大学，李轶，范影乐，庞全，多尺度自适应对比度变换的医学图像增强方法，中国专利,200710067693，2007-09-05；

12.  Hughes Aircraft Company，Frankot,et al.Multi-scale adaptive filter for interferometric SAR data，United States Patent，5488374，1996-1-30。

发明内容

    本发明所要解决的是当前多尺度DEM存在的数据冗余、不一致性、无法满足实际尺度需求以及在多尺度下无法保持基本地形特征等问题，提供一种能保持地形特征的自适应尺度DEM建模方法(Self-adaptive Scale DEM，SS-DEM)。

    本发明的技术方案所提供保持地形特征的自适应尺度DEM建模方法，包括如下步骤：

步骤1，根据原始精细尺度的数字高程模型构建不规则三角网，构建时采取由特征点和流域线构建不规则三角网的方式，包括根据最大高程误差z的不同取值从原始精细尺度的数字高程模型中提取地表的特征点，从原始精细尺度的数字高程模型提取流域线并简化，综合地表的特征点和流域线构建不规则三角网；

步骤 2，比较原始精细尺度的数字高程模型和步骤1所得不规则三角网，计算出最大高程误差z在不同取值下的均方根误差RMSE，对最大高程误差z和均方根误差RMSE进行曲线拟合，得到最大高程误差z和均方根误差RMSE之间的最佳函数解析表达式；

步骤3，建立不同尺度的特征点集并进行优化，具体步骤如下，

步骤3.1，根据制图规范中比例尺和等高距的关系，获取与尺度对应的最大高程误差z取值范围和均方根误差RMSE取值范围，并根据步骤2所得最大高程误差z和均方根误差RMSE之间的最佳函数表达式，确定最大高程误差z的准确取值范围，采用由特征点和流域线构建不规则三角网的方式提取出不同尺度下的特征点集， 

对任一尺度提取特征点集的实现过程为，根据与尺度对应的最大高程误差z在准确取值范围的不同取值，从原始精细尺度的数字高程模型中提取地表的特征点，然后从原始精细尺度的数字高程模型提取流域线并按照尺度进行简化，综合地表的特征点和流域线得到特征点集；

步骤3.2，如果在某尺度下，存在有两个特征点的间距小于该尺度下栅格数字高程模型的格网单元间距，就只保留其中最大高程误差z取值较大的特征点；

步骤3.3，如果一个特征点和简化后流域线的间距与简化后流域线长度的比值小于设定的阈值，去除该特征点；

步骤4，构建自适应尺度DEM数据结构，在应用中根据用户指定尺度获取特征点和流域线，生成数字地形模型并提取出地形参数，具体步骤如下，

步骤4.1，根据步骤3所得不同尺度的特征点集构建自适应尺度DEM数据结构，包括按照特征点出现的最粗尺度赋予各特征点尺度属性，得到包含所有特征点的尺度属性图；

步骤4.2，根据用户指定的尺度获取特征点和流域线，获取方式为，

如果用户指定尺度存在于尺度属性图中，则该尺度下的特征点由尺度属性图中所有具有该尺度以及更粗尺度属性的特征点组成；如果用户指定尺度不存在于尺度属性图中，则重复步骤3提取出该尺度下的特征点，并按照步骤4.1对尺度属性图进行更新，该尺度下的特征点由更新后的尺度属性图中所有具有该尺度以及更粗尺度属性的特征点组成；

根据用户指定尺度对流域线进行简化，得到简化后的流域线；

步骤4.3，根据步骤4.2所得用户指定尺度下的特征点和流域线生成受流域约束的不规则三角网作为数字地形模型，并提取出该用户指定尺度下的地形参数。

而且，在步骤4后执行步骤5进行精度评估，具体步骤如下，

步骤5.1，评估生成的数字地形模型的精度，包括将步骤4生成的数字地形模型与国家基础地理信息中心现有的多尺度数字地形模型进行比较，计算出均方根误差；

步骤5.2，评估提取出的地形参数和主要地形特征，所述地形参数包括平均坡度、平均地表粗糙度和曲率，所述主要地形特征包括流线匹配度和流线匹配误差；

步骤5.3，评估流域分析不同算法的精度，包括选取一种数学模拟曲面计算其理论单位汇水面积，再对该数学模拟曲面采取其他算法计算出单位汇水面积，将单位汇水面积与理论单位汇水面积进行比较，计算出均方根误差。

而且，采用Douglas-Peucker算法实现对流域线进行简化。

本发明创造性地将自适应尺度技术应用于数字地形分析领域，并且自适应尺度技术只在原始精细尺度DEM的基础上增加了一个索引结构，而没有更改原始的数据源，从而解决了多尺度DEM数据库带来的数据冗余、不一致性以及无法满足实际应用的尺度需求等问题，并能最大程度保持不同尺度下的基本地形特征。因此本发明特别适合应用于数字地形分析领域，实现自适应尺度DEM的地形建模与分析。

附图说明

图1为本发明实施例的流程图。

图2为本发明实施例步骤1的子流程图。

图3为本发明实施例步骤2的子流程图。

图4为本发明实施例步骤3的子流程图。

图5为本发明实施例步骤4的子流程图。

图6为本发明实施例步骤5的子流程图。

具体实施方法

本发明要解决的核心问题是：构建一个适用于从广域到局部范围通用的精细DEM数据库，根据应用需要从该数据库中提取多尺度的地形参数和特征，依据尺度信息保持不同级别地形特征，使得不同尺度下地形骨架基本一致。从理论上说，为进行多尺度地形分析，地表高程数据集需要满足三个基本条件：

1.      采用一个精细尺度的DEM数据库，表现所需的任何尺度的数据。

2.      地形参数 (平均坡度、粗糙度) 随着尺度的变化不能发生巨变，需要保持地形基本特征。

3.      所提取的流域网络、地形要素等不能发生巨变。即地貌特征, 流水线，几何形状和整体结构不能发生巨大变化。

实施例的流程参见附图1，该方法通过基于尺度的由特征点和流域线构建TIN的方法(scale-deterministic Compound Point Extraction，S-CPE)，即从精细尺度DEM中提取出多种更粗尺度下的特征点，建立一个自适应尺度DEM数据结构，以满足范围内任何尺度的应用需求。根据用户指定的尺度提取出的特征点，生成数字地形模型，并对模型本身以及提取出的地形参数进行验证。以下分步骤详细说明实施例的具体实施过程：

步骤1，对原始精细尺度的DEM采取现有由特征点和流域线构建TIN的方法（Compound Point Extraction，CPE），构建过程包括根据最大高程误差z的不同取值从原始精细尺度的数字高程模型中提取地表的特征点，从原始精细尺度的数字高程模型提取流域线并简化，综合地表的特征点和流域线构建不规则三角网。

具体实施时，由特征点和流域线构建TIN的方法可参见Zhou, Q., Chen, Y., 2011. Generalization of DEM for terrain analysis using a compound method. ISPRS Journal of Photogrammetry and Remote Sensing.Vol 66,1. 38-45。

根据这种CPE方法，实施例首先利用现有最大z-tolerance算法，z-tolerance指定了由该等级下的特征点集生成的TIN容忍范围内的最大高程误差，随着z值的变化从原始精细尺度的DEM中检索出不同级别的地表特征点。为了突出关键流域特征，再将补充的流域线从原始精细尺度的DEM中识别出来，并加入到特征点集，最后生成受流域约束的TIN。

为了便于实施，在此介绍最大z-tolerance算法的实现过程：最大z-tolerance算法采用一种迭代过程。首先选取DEM的四个角点，使得所构建的第一个TIN覆盖整个DEM的范围，然后连接一条对角线，将TIN划分为两个三角形，接着对TIN的每个三角形，计算格网中的每一点与所落入三角形面的高差，并确定差值最大的点。如果差值大于指定的z-tolerance，便标记该点并将其添加到TIN中，现存TIN中的每个三角形都被检测后，选中添加的点重新计算三角网，这个过程持续进行直到格网里的所有点都在指定的最大z-tolerance之内。

但是，最大z-tolerance算法所提取出的特征点集不能确定在地形分析中十分重要的流域类型，所以CPE方法增加了额外的流域特征线。为了获取这些流域线，可利用D8算法从DEM中进行提取。

参见图2，实施例中，为了减少数据的冗余，首先比较由最大z-tolerance算法提取出的特征点和D8算法提取出的流域线，若存在有重叠点，就由最大z-tolerance算法提取出的特征点中去除重叠点，不存在则直接保留所有最大z-tolerance算法提取出的特征点，结果记为点集A。然后根据现有技术中的Douglas-Peucker算法对流域线进行简化以提取出更有意义的流域特征线，如果原始DEM的分辨率为5m，可将Douglas-Peucker算法的阈值设为5m。经过Douglas-Peucker算法去除流域线上冗余的特征点后，简化后的流域线记为点集B，最终得到的特征点集包括点集A加上点集B，基于该特征点集（A+B）即可生成TIN。

步骤2，比较原始精细尺度的数字高程模型和步骤1所得不规则三角网，计算出最大高程误差z在不同取值下的均方根误差RMSE，对最大高程误差z和均方根误差RMSE进行曲线拟合，得到最大高程误差z和均方根误差RMSE之间的最佳函数解析表达式。

有两种属性与特征点集相关，也就是z和RMSE。z值指定了由该等级下的特征点集生成的TIN容忍范围内的最大高程误差。z值越大，特征点就越显著。RMSE描述了生成的TIN表面和原始DEM间的可容许的平均高程误差。描述这两种属性关系的最佳函数式可通过回归分析获取，其相关系数R ²应大于0.99。根据最佳函数式，在CPE中获取的z值大小能够根据指定的RMSE计算出来。

参见附图3，实施例比较原始DEM和步骤1得到的TIN，通过精度分析计算出不同z值下的均方根误差（RMSE），对z和RMSE进行曲线拟合，得到它们之间的最佳函数解析表达式。具体操作为，根据各尺度下生成的TIN生成DEM，将其与原始DEM进行精度分析比较，计算出RMSE。利用统计回归分析对z值和RMSE进行曲线拟合，得到一个最佳函数解析表达式，比如                                                

(a = 2.166, b = 0.173, c = 1.096)，

代表RMSE，

代表z。此时得到的相关系数R ²大于0.999，满足要求。根据得到的最佳函数式可将与每个特征点相关的z值与相应的RMSE进行转换。

步骤3，建立不同尺度的特征点集。

本步骤也利用了现有由特征点和流域线构建TIN的方法，并改进以建立不同尺度的特征点集，可称为基于尺度的由特征点和流域线构建TIN的方法(scale-deterministic Compound Point Extraction，S-CPE)。

参见附图4，实施例具体步骤如下：

    步骤3.1，根据据制图规范（例如2001年国家质量技术监督局发布的《数字地形图系列和基本要求》）中比例尺和等高距的关系，获取与尺度对应的z和RMSE取值范围，并根据步骤2得到的z和RMSE之间的最佳函数表达式，确定z的准确取值范围，再次利用CPE提取出不同尺度下的特征点。

根据指定的尺度确定相应的RMSE。按照国家质量技术监督局的标准(2001)，二分之一等高线间距定义为丘陵地区的的RMSE，而且z不能大于RMSE的两倍，所以要得到RMSE和z的有效取值范围，就要先确定指定尺度下的等高距。由于国家制图标准只规定了一些固定尺度的等高距，参见附图4，需要进行内插得到尺度和等高距间的连续分段函数解析式，这样对于任意一种尺度就能计算出相应的等高距，从而得到相应的RMSE和z的取值范围。然后根据步骤2得到的z和RMSE之间的最佳函数表达式，确定z的准确取值范围，接着采用CPE提取出特征点，从而得到该尺度下的特征点集。

在实施例中，如果比例尺为1:100000，得到丘陵地区的等高距为20m，则RMSE大于10m，z大于20m的特征点集应该出现在1:100000的尺度下，再通过步骤2的z和RMSE之间的最佳函数式得到更准确的z值范围，从而利用CPE提取出该尺度下的特征点集。而在1:200000比例尺下，由于制图标准中不存在该尺度，需要根据内插得到的尺度和等高距间的连续分段函数解析式进行计算，最终得到等高距为40m，则RMSE大于20m，z大于40m的特征点应出现在1:200000的尺度下，同样根据最佳函数式得到准确的z值范围，利用CPE提取出1:200000尺度下的特征点集。

因为同样采用CPE方法，本步骤实现与步骤1类似：对任一尺度提取特征点集的实现过程为，根据与尺度对应的最大高程误差z在准确取值范围的不同取值，从原始精细尺度的数字高程模型中提取地表的特征点，然后从原始精细尺度的数字高程模型提取流域线并按照尺度进行简化，综合地表的特征点和流域线得到特征点集。也就是说，在提取某一尺度的特征点集，实施例首先由最大z-tolerance算法提取出的特征点构成点集，然后比较由D8算法提取出的流域线，若有重叠点就从点集中去除；然后用Douglas-Peucker算法对D8算法提取出的流域线进行简化，将简化结果加入点集。提取不同尺度的特征点集时，区别只在于：（1）对应于不同尺度，z的取值范围不同。（2）对应于不同尺度，用Douglas-Peucker算法对D8算法提取出的流域线进行简化时，Douglas-Peucker算法的阈值根据与尺度相应的分辨率设定，因此取值不同。所以，只需从原始精细尺度的数字高程模型采用D8算法提取一遍流域线，然后在不同尺度下根据不同的阈值采取Douglas-Peucker算法进行简化即可。在后续步骤4.2中，也只需根据用户指定尺度对流域线进行简化，得到简化后的流域线。

    步骤3.2，如果两个特征点的间距小于该尺度下栅格DEM的格网单元间距，就只选择其中z值较大的点。

    在实例中，1:100,000尺度应用下两特征点的间距应该大于50m，而在1:250,000尺度下距离会增加到125m，如果两个特征点的间距小于该尺度下栅格DEM的格网单元间距，就只保留拥有较大z值的特征点。

    步骤3.3，如果一个特征点和流域线的间距与流域线长度的比值小于设定的阈值，去除该特征点。

本步骤所述流域线是指简化后的流域线。在实施例中，根据公式

，t表示点与流域线的距离指数，d表示点与流域线的间距，L表示流域线段的长度， Dt表示阈值，该阈值与尺度无关，一般设为0.5%。如果一点的距离指数小于阈值，就去除该特征点。

步骤3.1所得不同尺度的特征点集经步骤3.2和步骤3.3优化后，成为下一步骤的处理对象。

步骤4，构建自适应尺度DEM数据结构，在应用中根据用户指定的尺度获取特征点和流域线，生成数字地形模型并提取出地形参数。参见附图5，实施例具体步骤如下：

    步骤4.1，根据步骤3所得不同尺度的特征点集构建自适应尺度DEM数据结构。实现过程为，首先对各特征点赋予尺度属性，某个特征点的尺度属性为所有出现该特征点的尺度中的最粗尺度，表示从最精细尺度到该最粗尺度都包含该特征点；从而得到包含所有特征点的尺度属性图，构建出自适应尺度DEM数据结构。

在实施例中，假如选择的尺度范围是250m到30m，以10m递减，那么先对250m的最粗尺度采取步骤3提取出特征点集，并将这些特征点标记上250m的尺度属性。然后同样对240m的尺度采取步骤3进行提取，此时得到的特征点集包含了250m尺度下的特征点，将这些在上一级已经出现的特征点去除后，得到的就是新出现的特征点，将它们标记为240m的尺度属性。依次对所有尺度进行这样的处理，直到在30m的最精细尺度下新出现的特征点也进行了标记，不是特征点的单元赋予0，这样就得到了一幅包含所有特征点的尺度属性图。

    步骤4.2，根据用户指定的尺度获取特征点和流域线。如果用户指定的尺度存在于尺度属性图中，则该尺度下的特征点就由尺度属性图中所有该尺度以及更粗尺度属性的特征点组成；如果用户指定的尺度不存在，则重复步骤3通过S-CPE方法提取出该尺度下的特征点，并按照步骤4.1对尺度属性图进行更新。同时根据用户指定尺度采用Douglas-Peucker算法对流域线进行过滤，得到简化后的更有意义的流域线。

在实施例中，如果用户指定的尺度为100m，而尺度属性图中恰好存在该尺度值，则该尺度下的特征点集就由100m到250m尺度属性下的所有特征点组成，直接从尺度属性图中检索出这些特征点即可；如果用户指定的尺度为95m，而尺度属性图中不存在该尺度值，则重复步骤3提取出95m尺度下的特征点，并按照步骤4.1对尺度属性图进行更新。为提高效率起见，具体实施时也可以不对所有特征点重新赋予尺度属性，只对尺度属性为用户指定尺度的特征点进行调整，即只需过滤掉100m尺度下的特征点，对于剩余的特征点，将它们在尺度属性图中的尺度属性更改为95m即可。根据用户指定的尺度，如100m，采用Douglas-Peucker算法对流域线进行过滤时，其阈值可设为该尺度下的DEM分辨率大小，即100m。

    步骤4.3，根据步骤4.2得到的用户指定尺度下的特征点和流域线生成受流域约束的TIN，即可得到数字地形模型。并提取出该尺度下的地形参数，如坡度、坡向，曲率等。

    为了反映保持地形特征的效果，本发明实施例还设计了精度评估步骤。通过计算均方根误差可以从整体意义上描述地形参数和其真值的离散程度。

步骤5，对于生成的数字地形模型、提取出的地形参数以及流域分析中不同算法的精度，采用统计分析的方法评估精度。参见附图5，具体步骤如下：

步骤5.1，评估生成的数字地形模型的精度。将生成的数字地形模型与国家基础地理信息中心现有的多尺度数字地形模型进行比较，计算出均方根误差，该值越小说明精度越高。

步骤5.2，评估提取出的地形参数和主要地形特征，比如平均坡度、平均地表粗糙度、曲率等地形参数，以及流线匹配度和流线匹配误差等提供主要地形特征的流域保持指数。

对于地形分析的应用，关键是保持地形参数和主要地形特征。因而不仅对地形参数如平均坡度和平均地表粗糙度，也对流域保持指数做了比较，比如自适应尺度DEM方法和重采样方法的结果之间进行比较得到的流线匹配度（SMR）和流线匹配误差（SME），它们评价了在不同尺度下综合后的DEM与原始DEM提取得出的流域线之间的匹配度。

在实例中，计算平均坡度

和平均地表粗糙度K的公式如下：

 

其中，S表示坡度，A表示投影面积，n表示单元总数，下标i用于标识第i个单元，sec表示正割三角函数。

为了计算SMR，首先对从原始DEM提取出的河网生成缓冲区，缓冲半径依尺度而定，然后将综合后的DEM提取出的河网和该缓冲区进行叠置分析，最后得到SMR。计算公式为：

，其中

表示综合后的DEM提取出的河网落在缓冲区内的长度，L表示河网的总长。

   为计算SME，需要生成从原始DEM和综合后的DEM提取出的河网间所交叉的面积。计算公式为：，其中DA表示原始DEM和综合后的DEM提取出的河网间所交叉的面积，L表示河网的总长。

步骤5.3，评估流域分析不同算法的精度。选取一种数学模拟曲面，如椭球面和马鞍面，计算其理论单位汇水面积。再对该数学模拟曲面采取其他较为普遍使用的算法（例如D8，DEMON，Dinf，FMFD等）计算出单位汇水面积，将其与理论单位汇水面积进行比较，计算出均方根误差，该值越小说明精度越高。

以上所述仅为本发明中的一个实施例，并不用于限制本发明。凡在本发明的精神与原则之内，所做的任何修改，改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种保持地形特征的自适应尺度DEM建模方法，其特征在于，包括以下步骤：

步骤1，根据原始精细尺度的数字高程模型构建不规则三角网，构建时采取由特征点和流域线构建不规则三角网的方式，包括根据最大高程误差z的不同取值从原始精细尺度的数字高程模型中提取地表的特征点，从原始精细尺度的数字高程模型提取流域线并简化，综合地表的特征点和流域线构建不规则三角网；

步骤2，比较原始精细尺度的数字高程模型和步骤1所得不规则三角网，计算出最大高程误差z在不同取值下的均方根误差RMSE，对最大高程误差z和均方根误差RMSE进行曲线拟合，得到最大高程误差z和均方根误差RMSE之间的最佳函数解析表达式；

步骤3，建立不同尺度的特征点集并进行优化，具体步骤如下：

步骤3.1，根据制图规范中比例尺和等高距的关系，获取与尺度对应的最大高程误差z取值范围和均方根误差RMSE取值范围，并根据步骤2所得最大高程误差z和均方根误差RMSE之间的最佳函数表达式，确定最大高程误差z的准确取值范围，采用由特征点和流域线构建不规则三角网的方式提取出不同尺度下的特征点集，

对任一尺度提取特征点集的实现过程为，根据与尺度对应的最大高程误差z在准确取值范围的不同取值，从原始精细尺度的数字高程模型中提取地表的特征点，然后从原始精细尺度的数字高程模型提取流域线并按照尺度进行简化，综合地表的特征点和流域线得到特征点集；

步骤3.2，如果在某尺度下，存在有两个特征点的间距小于该尺度下栅格数字高程模型的格网单元间距，就只保留其中最大高程误差z取值较大的特征点；

步骤3.3，如果一个特征点和简化后流域线的间距与简化后流域线长度的比值小于设定的阈值，则去除该特征点；

步骤4，构建自适应尺度DEM数据结构，在应用中根据用户指定尺度获取特征点和流域线，生成数字地形模型并提取出地形参数，具体步骤如下：

步骤4.1，根据步骤3所得不同尺度的特征点集构建自适应尺度DEM数据结构，包括按照特征点出现的最粗尺度赋予各特征点尺度属性，得到包含所有特征点的尺度属性图；

步骤4.2，根据用户指定的尺度获取特征点和流域线，获取方式为：

如果用户指定尺度存在于尺度属性图中，则该尺度下的特征点由尺度属性图中所有具有该尺度以及更粗尺度属性的特征点组成；如果用户指定尺度不存在于尺度属性图中，则重复步骤3提取出该尺度下的特征点，并按照步骤4.1对尺度属性图进行更新，该尺度下的特征点由更新后的尺度属性图中所有具有该尺度以及更粗尺度属性的特征点组成；

根据用户指定尺度对流域线进行简化，得到简化后的流域线；

步骤4.3，根据步骤4.2所得用户指定尺度下的特征点和流域线生成受流域约束的不规则三角网作为数字地形模型，并提取出该用户指定尺度下的地形参数。

2.根据权利要求1所述自适应尺度DEM建模方法，其特征在于：在步骤4后执行步骤5进行精度评估，具体步骤如下：

步骤5.1，评估生成的数字地形模型的精度，包括将步骤4生成的数字地形模型与国家基础地理信息中心现有的多尺度数字地形模型进行比较，计算出均方根误差；

步骤5.2，评估提取出的地形参数和主要地形特征，所述地形参数包括平均坡度、平均地表粗糙度和曲率，所述主要地形特征包括流线匹配度和流线匹配误差；

步骤5.3，评估流域分析不同算法的精度，包括选取一种数学模拟曲面计算其理论单位汇水面积，再对该数学模拟曲面采取其他算法计算出单位汇水面积，将单位汇水面积与理论单位汇水面积进行比较，计算出均方根误差。

3.根据权利要求1或2所述自适应尺度DEM建模方法，其特征在于：采用Douglas-Peucker算法实现对流域线进行简化。

Images