专利名称：一种顾及地表水文上下文的多尺度dem建模方法
技术领域：
本发明涉及一种多尺度DEM建模方法，尤其是涉及一种顾及地表水文上下文的多尺度DEM建模方法。
背景技术：
数字高程模型(Digital Elevation Model，DEM)作为地球表面地形的数字描述和模拟方式，是空间数据基础设施的重要组成部分，在水利、测绘、国防建设等涉及地球三维信息的科学研究领域和生产部门具有广泛的应用。为满足不同层次的需要，各国均采用多种比例尺共存的方式建立多级尺度(或多分辨率)DEM数据库。例如美国联合太空总署(NASA)、国防部国家测绘局(NIMA)以及德国和意大利航天局，建立了 90米分辨率的SRTM3DEM、30米分辨率的SRTMlDEi^P 15米分辨率的ASTER DEM，英国建立了 25米分辨率的Landmap DTM和5米分辨率的Bluesky DTM，我国也相继建立了全国范围的1:100万、1:25万、1:5万DEM以及省市级的1:1万DEM (NFGIS-DEm)0随着合成孔径雷达、激光雷达、激光扫描系统等空间数据获取手段的更新和提高，地形数据的获取变得高效和快捷，一方面极大地满足了人们在各领域进行地球科学研究的潜在需求，另一方面也带来了地形数据在数据量、时效性和复杂性等方面的急剧增长问题。如对地观测系 统E0S(Earth Observing System)通过TERRA、AQUE和QULA等七颗卫星可实现对地球全天候、全天时并且不受云雾干扰的观测，仅TERRA每日获取的遥感数据就达TB级，而TERRA与AQUA卫星的组合还实现了全球每天最少4次的数据更新。由于各级比例尺DEM数据量呈几何级数快速递增，海量数据的存储、更新和一致性维护需求对分比例尺DEM建库这一做法的可持续性提出了严峻的挑战。采用DEM综合技术实现DEM由大比例尺(高精度)向小比例尺(低精度)转换，进而提供尺度灵活的地形信息服务，是现代DEM生产和管理的必然趋势。根据DEM数据组织形式的不同，目前国内外研究DEM综合算法大体上可分为四类:(I)基于离散点的DEM综合。这类方法的核心思想是，认为地表上的高程点在空间布局上与周围一定范围内的其它高程点具有相关性。代表性的研究主要有点面距法、地形描述误差法、空间平面法矢量夹角法。这些方法利用高程点的空间相关性来实现DEM综合，但是相邻点子的相互决定使得点子的处理顺序对综合结果有较大的影响，也不可避免的缺乏对地表形态的整体把握。(2)基于等高线的DEM综合。等高线综合法主要采用曲线化简算法及其改进算法对等高线进行综合，去除细部弯曲，保留主要弯曲，如二维Douglas-Peucker法、L1-Openshaw法、QTM法、小波法、分形综合法、遗传算法。由于相邻等高线通常具有不同程度的相似性，并且等高线以成组的方式来表达地貌，也有学者从等高线间的相互关系对等高线进行结构化综合，如高程带综合法、结构线跟踪法、地表水系结构化法、渐进式综合法。等高线通过等值线上的连续性和等值线间的突变对地表进行“塑形”，等高线的综合过程虽然利用地形结构线予以辅助，但实质是对投影到平面空间的二维信息进行综合，并非真正意义的三维综合，因此，需要对综合后等高线间的拓扑冲突进行进一步的处理。(3)基于规则格网的DEM综合。这类综合核心思想是将规则格网看作是数字灰度图像，位于每一个栅格的灰度值与位于该处或内插于该处的地貌点高程值对应，通过对灰度值的处理实现DEM综合，如滤波法、小波分析方法、信息论法、四叉树局部熵法、地图代数法。由于格网DEM自身数据结构的限制，地形的细节总是被抑制于网格中，因此，在综合时必须加入更多的特征信息。目前大多数基于格网DEM的综合方法都需加入地形特征点或特征线来增强局部栅格与全局地表形态的协调关系。(4)基于剖面线的DEM综合。Mahes (1998)提出以等间距竖直剖面切割地表，对得到的剖面线进行线性化简，保留地 表的主要凹凸特性。这种方法对地形信息的化简提供了新的思路，然而目前相关研究较少，相关技术和应用前景还有待进一步研究验证。在上述研究中，DEM综合通常表现为在数据处理层面上的数据压缩。若要满足人们对地形信息的处理和分析需要，就要求不同比例尺DEM数据能够对地形信息具有尺度适应的“保真”，在操作决策层面上则需要以DEM所蕴含的地貌特征作为地理学依据，并将其作为综合行为的空间上下文与数据模型关联起来进行协同化简。就地貌要素而言，地貌区域地形起伏所反映出来的地貌单元在地质构造、水文特性等方面的地学现象和规律，属于DEM综合的空间上下文范畴。特别是对于发育成熟的地貌，地形的主要谷地和山脊结构线具有明显的集水性和分水性，这既是地理学研究中地表水文分析和水文要素提取的理论基础，也是DEM综合过程中需要考虑的上下文约束。国内外学者的研究虽不同程度的关注了空间上下文对制定综合方案的作用，但仍未形成有效的算法将地理层次的上下文信息与几何层次的数据处理整合到一起。这些算法要么在综合前通过地形结构的层次分析，局限于对被舍弃的次要谷脊局部进行几何平滑，要么在综合后补充未化简的结构线信息，对综合结果进行结构模式的补充。在DEM综合过程中如何关联空间上下文进行一致化的化简，目前尚缺乏适用的算法模型。
发明内容
本发明的上述技术问题主要是通过下述技术方案得以解决的:一种顾及地表水文上下文的多尺度DEM建模方法，其特征是，包括以下步骤:步骤I，对原始高精度DEM进行汇水分析，提取DEM所描述地表形态的集水线和分水线，建立原始高精度DEM的水文网络；步骤2，在步骤I建立的水文网络的基础上解析原始高精度DEM所描述地貌对象的水文上下文:利用地表汇水响应区域的汇入关系建立集水线网络的层次结构，然后基于地表山谷和山脊的耦合关系进行分水线与集水线的匹配，根据集水线网络的层次结构确定分水线网络的层次结构；步骤3，对原始高精度DEM的地表高程点进行水文语义增强:建立地形高程点的水文特征泛化模型，根据步骤2获得的水文上下文计算地表高程点的水文贡献指数；步骤4，根据步骤I所得水文网络对原始高精度DEM按流域进行分解，以流域为单位提取地形特征点集:采用流域边界凸壳和地表高程点构建不规则三角网的方式计算地表高程点存在与否所引起的DEM地形描述误差，比较地表高程点的地形描述误差大小，然后从地表高程点集中提取地貌塑形特征点并建立地形特征点队列，在此基础上根据用户指定尺度确定的综合程度控制参数，获取相应尺度的地形特征点集；步骤5，从地形特征点集中分离出集水线和分水线上的高程点，对原始高精度集水线网络和分水线网络进行相应尺度的化简；步骤6，根据步骤4和步骤5所得到的相应尺度下的地形特征点集和化简的水文系统，建立受约束的不规则三角网TIN，在此基础上内插规则格网作为综合后的DEM ；步骤7，计算综合后DEM的精度，若精度达不到用户设定尺度的综合程度，采用弹性寻位策略自动调整综合程度控制参数，获取相应尺度的地形特征点子集，重复步骤5 步骤7。在上述的一种顾及地表水文上下文的多尺度DEM建模方法，所述步骤I中，建立原始高精度DHM的水文网络的具体步骤如下:步骤1.1，根据DEM规则格网的高程差异建立流向矩阵，以未标记的任意位置栅格作为种子点，按八邻域逆向迭代搜索汇入栅格，以只有流出栅格没有流入栅格的点作为径流源头点，然后按照径流方向追踪，采用改进D8算法建立径流累积矩阵，具体是:首先遍历规则格网，计算每个格网向八领域格网的高程差，把高程差最大的方向标识为水流方向，得到流向矩阵；然后随机选取未标记的格网作为种子点，按八邻域逆向搜索汇入栅格，以只有流出栅格没有流入栅格的点作为径流源头点。再按径流方向追踪，水流路径上某一栅格的径流累积量为上游栅格的径流累积量加I ;对于多个方向汇入的栅格，径流累积量在各汇入方向累积量逻辑和的基础上加I，得到径流累积矩阵；步骤1.2，根据步骤1.1所得径流累积矩阵提取每一个格网点的径流累积量，提取大于设定阈值的稀疏栅格作为集水线栅格，结合步骤1.1中建立的流向矩阵从集水线栅格中提取集水线的源头，然后再根据步骤1.1中建立的径流累积矩阵以集水线的源头作为种子点，逆向延伸集水线直至径流源头；
步骤1.3，对步骤1.2建立的集水线稀疏栅格进行矢量化，然后利用步骤1.2建立的流向矩阵构建集水线网络；步骤1.4，提取以集水线的源头与叉头或者叉头与出口之间的线段作为流域的集水核，然后根据步骤1.1得到的流向矩阵和径流累积矩阵将所有流入集水核的区域标识为对应的汇水响应区域；步骤1.5，对步骤1.4得到的汇水区域边界进行锐化和骨架线提取，建立分水线网。在上述的一种顾及地表水文上下文的多尺度DEM建模方法，所述步骤2中，解析原始高精度DEM所描述地貌对象的水文上下文的具体步骤如下:步骤2.1，以步骤1.4获得的集水核作为集水系统的基本数据组织单元，利用地表汇水响应区域的汇入关系建立集水网络的层次结构树，即有:父子关系的集水核，父级集水核的汇水响应区域包含所有子级集水核的汇水响应区域；以及兄弟关系的集水核，互为兄弟的集水核其汇水响应区域汇入同一父级汇水响应区域；步骤2.2，将分水线的起点与叉点或叉点与出口之间的线段作为流域的分水核，基于分水线与流域的耦合关系，利用汇水响应区域的汇入关系建立分水核网络的层次结构树，即有:父子关系的分水核，父级分水核与父级汇水响应区域的边界重合，子集分水核与子集汇水响应区域的边界重合且与父级汇水响应区域的边界相交；兄弟关系的分水核，互为兄弟关系的分水核其对应的汇水响应区域同为兄弟关系，且均不与父级汇水响应区域的边界重合；在上述的一种顾及地表水文上下文的多尺度DEM建模方法，步骤3对原始高精度DEM的地表高程点进行水文语义增强，具体实现方式是根据步骤2得到的集水线网络与分水线网络的层次结构建立流域指数函数，然后通过计算流域指数和微地貌指数两种水文贡献指数对地表高程点的水文特征进行泛化。在上述的一种顾及地表水文上下文的多尺度DEM建模方法，所述步骤4中，获取地形特征点集的具体实现步骤如下:步骤4.1，将原始高精度DEM按流域划分为若干个特征提取单元，提取流域边界凸壳构建初始不规则三角网TIN，获取除边界外的剩余高程点作为候选特征点集；步骤4.2，随机抽取候选特征点集中的一个候选特征点，结合步骤3得到的地表高程点水文贡献指数计算出其到不规则TIN表面的加权点面距，选取加权点面距最大的高程点作为当前第I特征点，加入地形特征点序列中；步骤4.3，将该 地形特征点插入不规则TIN中进行重构，从候选特征点集中删除该特征点，重复步骤4.2 步骤4.3，直至加权点面距小于给定的阈值或者所有的高程点都进入特征点序列。在上述的一种顾及地表水文上下文的多尺度DEM建模方法，所述步骤5中，化简水文系统的具体实现方式是，根据用户指定尺度确定地形特征点的选取数量，从步骤4得到的特征点序列中获取相应子集，然后从地形特征点子集中分离出位于集水线和分水线上的地形特征点，分别沿着初始集水线和分水线路径，根据子集中的集水线和分水线地形特征点重构新的集水线和分水线，重构过程中保持集水线和分水线与一般地形特征点的拓扑关系不变。在上述的一种顾及地表水文上下文的多尺度DEM建模方法，所述步骤7中，采用弹性寻位策略自动调整综合程度控制参数是通过计算步骤6所得的综合后DEM的精度来进行，具体实现步骤如下:步骤7.1，根据多尺度DEM资源库建立尺度与DEM高程描述差的关系函数F_s，根据用户设定的尺度获取相应的F_s(x)，比较综合前后的DEM高程描述误差，计算出高程描述误差的差值F(X)；步骤7.2，采用弹性寻位策略基于以下约束条件控制DEM综合程度:设定容差阈值为ε ,若I F_s (X)-F(X) |> ε且F_s (x) >F(x)，则增加从步骤4特征点序列中提取的特征点集，反之，若I F_S (χ) -F(X) I > ε且F_S (χ)〈F (χ),则减少从步骤4特征点序列中提取的特征点个数，重复步骤4 步骤7,直至I F_s (χ)-F (χ) I〈 ε结束。因此，本发明具有如下优点:(I)可有效保持原始高精度DEM的地形和水文特征，并能保证多尺度下DEM的地形与水文特征的拓扑和逻辑一致性；(2)提高了多尺度DEM建模的自动化和生成效率，并能满足实际应用中对尺度自适应DEM的需要。
图1是本发明顾及水文上下文的DEM自动综合流程图。图2是本发明集水线逆流延伸流程图。图3a是本发明集水系统示意图。图3b是本发明集水系统Shreve编码图。图3c是本发明集水系统层次结构图。图4a是本发明分水系统及分水线汇水区域剖分示意图。图4b是本发明分水系统层次结构图。图5是本发明特征点选取流程图。图6是本发明基于弹性寻位确定最佳综合阈值的示意图。图7是本发明自动控制DEM综合程度的流程图。图8a是本发明实验结果DEM的等高线图。图8b是水文增强法实验结果DEM的等高线图。图8c是VIP法实验结果DEM的等高线图。图8d是Decimation法实验结果DEM的等高线9是四种方法不同综合程度的地形湿度指数曲线。
具体实施例方式下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。实施例:一、首先，介绍一下本发明涉及的理论知识。1、理论基础。进行DEM数据综合，从较大比例尺(高分辨率)DEM派生得到较小比例尺(低分辨率)DEM,是对地形信息的概括和简化，而不是简单的细节特征弃除。因此，DEM综合尽管在操作层面上通常表现为对数据的压缩，但在操作决策层面上应以DEM所蕴含的地形特征作为地理学依据，将其作为综合行为的空间上下文与数据模型关联起来进行协同化简。就地貌要素而言，地形起伏反映了发育成熟的谷地和山脊具有明显的集水性和分水性，这既是地理学研究中地表水文分析和水文要素提取的理论基础，也应是进行DEM综合操作的空间上下文。DEM综合的本质是地形信息从较大比例尺到较小比例尺的信息流转换。地表形态是一个空间集合体的概念，由不同的几何面所围成，具有一定的体积、大小、高低、形状等外物特征。任何复杂的地形都可认为是由平地(包括高平地、阶地、低洼平地)、坡地和洼地(包括孤立洼地和组合洼地)组成，最终这些基本地形单元又都是通过区域中高程点相对于周围高程点的高低来“定位”与“定性”描述。因此，地形高程点具有集合意义，点集整体构成了地表的基本形态，也构成了地貌空间几何体的支架，包括地貌区域沟谷网络、分水系统。因此，DEM综合的关键在于如何选取对地表形态的描述以及地形上下文特征保持具有较大贡献的高程点。
地形空间上下文对DEM综合操作的约束作用主要体现在对地形特征点的选取过程中。换言之，保持地表水文特征的DEM综合，需关联地形特征点对地形水文特征的贡献程度从原始DEM中选取地形特征点。而当地形特征候选点对地表塑形的贡献程度相近时，需优先选取对地形水文特征表达贡献程度较大的点。基于此，本发明提出并实现了一种新的DEM综合方法。通过对DEM进行径流分析提取地形表面的水文响应区域，建立汇水区域的层次化组织模型，在此基础上根据地表水文系统的层次结构建立地表高程点的水文特征泛化模型，并通过量化地表高程点的水文贡献程度和进行DEM的水文语义增强，选取地形特征点集。可以看出，通过这种方式，DEM的综合过程始终考虑地表水文上下文的约束作用，DEM从大比例尺(高分辨率)向小比例(低分辨率)的数据转换朝着地形信息和水文信息协同化简的方向进行。2、高程点水文特征的泛化.
集水网络和分水网络上的高程点与其他区域的一般高程点具有不同的语义内涵，对DEM数据中的地表水文信息起显式表达的作用。因此，在DEM综合过程中考虑地表水文上下的空间约束，可通过对水文网络系统上的高程点进行语义增强来实现。这需要建立对地表水文上下文进行信息泛化的基础上。我们将水文信息分为两类流域水文信息和水文微地貌信息两部分， 通过定义流域指数和微地貌指数来泛化地表高程点的水文特征值。(I)流域指数由水文学基本理论可知，地表流域具有明显的层次嵌套结构，每一级汇水区域由若干次一级汇水子区域组成。以集水线的源头与叉头或者叉头与出口之间的线段作为流域的集水核，分水线的起点与叉点或叉点与出口之间的线段作为流域的分水核。对任意集水核或分水核单元I，设其对应的汇水区域面积为Area(I)，整个流域的汇水面积为Area(S)，则集水核或分水核单元上任一高程点Pt_l的流域指数a (pt_l)的计算公式为:a (pt_l) =Area(I)/Area(S) (I)其中，集水核和分水核汇入点的流域指数取与之连结的最高级别集水核或分水核高程点贡献值。(2)微地貌指数集水线与分水线的水文微地貌特征同时受局部区域地质构造和地貌营力的作用而形成，其形态可通过地形参数来描述。对任意集水核或分水核片段1，定义I上任一高程
点Pt_l的微地形贡献指数β (pt_l)，为沿水文线高程下降方向的高程变化率pi_/ )
与垂直水文线方向的闻程变化率的矢量和，即:
Pipt _ I) = Slope{pt _ In) + Slope(pt—/丄)(2)3、地形特征点的提取。地表是三维连续的体对象，对于地形信息的化简应在三维立体空间中进行。因此，DEM综合过程中判断是否为地形特征点，主要考察该高程点对地形起伏的贡献量。在DEM质量分析领域，通过DEM模拟地面与实际地形的差异评价DEM的地形描述精度。本发明借鉴该思想，将高程采样点存在与否引起的DEM地形表达差异作为地形特征点的选取依据。考虑到不规则三角网可更灵活的模拟地表起伏，地形特征选取在TIN DEM基础上进行。基本思想是:首先提取区域边界，构建初始TIN，以除边界之外的所有高程点作为候选特征点；然后计算所有候选点到TIN表面的三维点面距，
选取点面距最大的高程点作为第I地形特征点；将第I地形特征点插入TIN中进行重构，然后递归选取和内插第2、3……η地形特征点，直至选出的地形特征点的点面距小于给定的阈值或为O。上述特征点选取过程将地表高程点相对于TIN的点面距作为DEM地形表达差异的量化指标。为了保留地貌形态的集水与分水特征，本发明结合高程点的水文特征指数来计算加权点面距Cli:di=delta_d(i；tin) X (1+ a ^ β ^ (3)其中delta_d(i，tin)为候选特征点i相对于当前TIN的三维点面距，a i和β i为该点的水文贡献指数，分别为流域指数和微地貌指数。集水线和分水线上高程点的^和I通过公式⑴、⑵计算得到，其他区域高程点的a i和β i取值为O。地形特征点选取过程完成后，得到的是对地形由粗到精描述的特征点序列。可根据用户指定的尺度与综合前后DEM高程差描述误差的对应关系，自前向后从序列中选取地形特征点集。二、下面，具体介绍一下结合上述理论知识的具体实施步骤。如附图1所示:(I)根据原始高精度DEM规则格网的高程差异建立流向矩阵M_dir和径流累积矩阵皿_8(^。具体建立过程如下:首先遍历规则格网，计算每个格网向八领域格网的高程差，把高程差最大的方向标识为水流方向，得到M_dir。然后随机选取未标记的格网作为种子点，按八邻域逆向搜索汇入栅格，以只有流出栅格没有流入栅格的点作为径流源头点。再按径流方向追踪，水流路径上某一栅格的径流累积量为上游栅格的径流累积量加I ;对于多个方向汇入的栅格，径流累积量在各汇入方向累积量逻辑和的基础上加1，得到M_acc。该过程是对D8算法的改进，具体实现程序不详细介绍。(2)从径流累积矩中提取汇水量大于给定阈值的栅格，沿水流方向追踪得到初始集水线，然后进行集水线的逆流延伸。具体实现方法如下:①首先提取当前集水线的源头Pt_i作为起点，然后根据流向矩阵M_dir和汇流累积矩阵1_&(^，选取八邻域内流向该起点且汇水量Acc最大者作为第一候选点Pt_a，SPAcc (Pt_a) =max (Acc (流入栅格));②再根据DEM和汇流累计矩阵M_acc，选取八邻域内高程Z不小于Pt_1、汇水量小于Pt」且汇水量最大者作为第二候选点Pt_b，即Z(Pt_b)>Z(Pt_i)，Acc (Pt_b) <Acc (Pt_i),且 Acc ( Pt_b) =max (Acc (八邻域))；③比较Pt_a和Pt_b的汇水量，若两者相等，选Pt_a作为Pt_i的上游延伸点Pt_1-1，否则，选汇水量大者作为上游延伸点即有Acc(Pt_1-l)彡Acc (Pt_a)。④以延伸点Pt_i_l作为新的起点，重复① ④,直至流向延伸点的上游点汇流量为 O,即有 Acc(PtJ)) =0。集水线逆流延伸的实现过程见附图2。该过程可以有效解决径流溯源的分叉问题:当出现多个支流汇入或存在并行径流与谷地延伸方向并肩时，始终取地表水文响应最强的径流路径作为谷地延伸方向。(3)分水网络与集水网络成对偶关系，根据流向矩阵和径流累积矩阵将所有流入集水线的栅格区域标志为对应的水文响应区域，对水文响应区域边界进行锐化和骨架线提取，即可得到对应的分水线。汇水区域边界锐化和骨架线提取属我们设计的基于栅格的地图代数处理方法，具体实现过程不详细介绍。(4)采用Shreve法进行水文系统的组织。将集水线的源头与叉头或者叉头与出口之间的线段作为流域的集水核(附图3a)，分水线的起点与叉点或叉点与出口之间的线段作为流域的分水核(附图4a)，对应于Shreve编码体系中的编码单元(附图3b)。集水核之间的结构关系由汇水区域的汇入关系确定，数据组织模型如附图3c所示。(5)根据汇水区域之间以及汇水区域与其子区域之间可共享分水岭的地貌特点，分水系统的结构关系通过分解分水核涉及的多重汇水区域实现。设分水核AB为集水线a(附图3a)所在子汇水区域的分水岭，若其同时也属于集水核g (附图3a)所在父汇水区域的分水岭，则以分水核AB作为分水系统的最小单元。定义分水网络层次关系与汇水区域层次关系一致，分水系统的组织模型如附图4b所示。(6)在建立集水网络和分水网络的层次结构树基础上，计算各集水核和分水核高程点的水文特征泛化指数，即流域指数α和微地貌指数β。流域指数α的计算公式为:a (pt_l) =Area(I)/Area(S)。其中，I表示某一集水核或分水核，pt_l为对应集水核或分水核上的高程点；Area(l)为集水核或分水核的水文响应区域面积，Area(S)为整个流域的面积。当I在水文网络的层次树结构中为叶节点时，则I的水文响应区域取以I为集水核的流域；当I在水文网络的层次树结构中有子叶节点时，则I的水文响应区域取集水核I的流域与各子集水核流域之和。微地貌指数β 的计算公式为:β( )! 二 Slope[pt + Slope(pt J 0 实中，SlopeiptJ )表示沿水文线高程下降方向的高程变化率，Slopeipf丄)表示垂直水文线方
向的高程变化率。(7)把高程点的水文特征泛化指数作为特征选取的权系数，计算高程点对地表形态的塑形贡献量，建立地形特征点序列。具体过程如下:①首先提取原始高精度DEM的三维边界构建受限初始TINtl,把除构建受限TINtl之外的所有高程点作为候选特征点Cand(η)。②在此基础上，结合流域指数和微地貌指数计算候选地形特征点对地表形态的贡献量，计算公式为Mi=(Ielta-Clatin)X (1+α j+β ^ 公式中的(Ii为候选地形特征点i相对于当前TINi的加权点面距，delta_d(i，tin)为候选地形特征点i相对于当前TINi的三维点面距，^和1为该点的水文贡献指数，分为别流域指数和微地貌指数。集水线和分水线上候选地形特征点的a jP β i在过程(4)中计算得到，其他区域候选地形特征点的a jP β i取值为O。选取加权点面距最大的候选点作为第I地形特征点Pt1，放入地形特征点序列Q中。③然后将选取的地形特征点Pt1插入TINtl中重新建立三角网，得到TIN1,并从候选特征点集中删除该点。④计算剩余候选点对地表形态的贡献量，递归选取、内插第2、3......η地形特征点，
直至选出的地形特征点的点面距为O或小于给定的阈值，或者候选点已全进入地形特征点序列Q中，过程结束。地形特征点选取流程图见附图5。
(8)设地形特征点集的初始数量为ε，把ε作为DEM综合程度控制参数，从地形特征点队列Q中选取特征点子集P，再从P中提取来源于集水网络和分水网络的地形特征点P ’，然后对集水线和分水线进行化简。化简后的水文网络要求保持与选取的特征点集具有拓扑一致性。具体实现方法如下:对于某一集水核I，从源头向末端顺序取3个高程点Ptp1、Pti和Pti+1，Pti^1为源头点、分叉点或集水核上需保留的高程点，记为Pt_A、Pt_B、Pt_C。若Pt_B e P’，或以Pt_A、Pt_B和Pt_C建立三角形八八8(:，有3 1:—6￡ 且？丨_8与三角形Λ ABC有交点，则在集水核I的径流路径上保留Pti,设Pt_A=Pti; PLB=Ptp1 ;否则，在集水核I的径流路径上删除Pti,设Pt_B=Pti+1。取集水核I的下一个高程点Pti+2，设Pt_C=Pti+2，重复本过程，直到Pt_B为集水核I的末端。分水线的化简过程与之相同。(9)用化简后的集水线和分水线、以及选取的地形特征点集P生成受限TIN，再采用原始高精度DEM分辨率内插规则格网得到综合后DEM，计算综合前后DEM高程差的均方根误差RMSE。(10)计算国家地理信息中心中相似地貌类型区域的多尺度DEM数据的高程差RMSE，建立精度评估标准。若多尺度DEM数据库中存在用户指定比例尺的DEM数据，则从精度评估标准中直接取出对应高程差RMSE ;若不存在用户指定比例尺的DEM数据库，则根据多尺度DEM建立尺度与DEM高程差RMSE之间的关系曲线函数F_s，然后根据用户指定的尺度获取相应的高程差RMSE。建立精度评估标准的过程可在建立多尺度DEM数据库的过程中完成，也可在综合过程前期进行。设根据精度评估标准得到的目标尺度高程差RMSE为F_s (χ)，计算综合前后DEM高程差RMSE，设为F (χ)，比较F_s (χ)和F (χ)的大小。(11)根据F_s(x)与F(X)的比较结果，采用弹性寻位的策略自动调整高程点的选取数量。对于初始选取高程点数量ε，设压缩率的变化步长delta_ ε = ε，综合状态记录器kl=l,k2=l。当 F_s (χ)〈F(χ)时，取 k2=l,若 k2==kl,则有 delta_ ε =delta_ ε ;若1^2 古 kl，则 delta_ ε =delta_ ε *0.5 且 kl=k2。当 F_s (x) >F(x)时，取 k2=_l,若 k2==kl,则 delta_ε =delta_ ε ;若 k2 幸 kl,贝丨J delta_ ε =delta_ ε *0.5 且 kl=k2。利用弹性寻位策略自动调整地形特征点集P的数据量ε，当F_s(X)〈F(X)时，ε = ε +delta_ ε ;当F_s (x) >F (x)时，ε = ε -delta_ ε。利用弹性寻位策略调整特征点集数据量ε的示意图如附图6所示。(12)根据特征点集数据量ε从地形特征点队列中选取新的地形特征点集，重复第(8)步 第(12)步,直到F_s (χ) ^ F (χ)或达到跳出循环的条件,过程结束。自动控制DEM综合程度的流程如图7所示。(13)选取地貌发育成熟的西南某高山样区DEM数据进行实验，栅格分辨率为15m*15m，样区范围约为3.6km*4.2km。分别采用本发明提出的顾及地表水文特性DEM综合方法、水文增强综合法、ArcGIS软件提供的VIP综合法以及Decimation综合法对原始精细DEM进行处理。选取5%的地形特征点内插成TIN，然后引绘等高线考察综合结果的地形描述精度。综合前后等高线的叠 加如附图8a_8d所示。从等高线图的叠加效果可以看出，在特征点压缩为原始数据5%的情况下，本发明提出的方法能最大限度的保留原始地形特征，传统方法综合结果的地形谷脊区域等高线结构特征生硬，而VIP法和Decimation法都较大程度的损失了原始数据的地形信息。采用四种算法分别选取50%、25%、10%和5%的地形特征点，然后内插生成格网DEM计算不同综合程度DEM的地形湿度指数，考察四种方法的水文特征保持效果。各算法综合结果地形湿度指数标准差统计曲线如附图9所示。对比不同特征点选取率条件下的地形湿度指数标准差曲线可知，由于在综合过程中，本发明的方法对谷地和山脊的重要性判断以及地形特征点的识别都顾及了地表水文上下文信息，综合结果可较好的保持地形湿度特性。本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的 精神或者超越所附权利要求书所定义的范围。
权利要求
1.一种顾及地表水文上下文的多尺度DEM建模方法，其特征在于，包括以下步骤: 步骤1，对原始高精度DEM进行汇水分析，提取DEM所描述地表形态的集水线和分水线，建立原始高精度DEM的水文网络； 步骤2，在步骤I建立的水文网络的基础上解析原始高精度DEM所描述地貌对象的水文上下文:利用地表汇水响应区域的汇入关系建立集水线网络的层次结构，然后基于地表山谷和山脊的耦合关系进行分水线与集水线的匹配，根据集水线网络的层次结构确定分水线网络的层次结构； 步骤3，对原始高精度DEM的地表高程点进行水文语义增强:建立地形高程点的水文特征泛化模型，根据步骤2获得的水文上下文计算地表高程点的水文贡献指数； 步骤4，根据步骤I所得水文网络对原始高精度DEM按流域进行分解，以流域为单位提取地形特征点集:采用流域边界凸壳和地表高程点构建不规则三角网的方式计算地表高程点存在与否所引起的DEM地形描述误差，比较地表高程点的地形描述误差大小，然后从地表高程点集中提取地貌塑形特征点并建立地形特征点队列，在此基础上根据用户指定尺度确定的综合程度控制参数，获取相应尺度的地形特征点集； 步骤5，从地形特征点集中分离出集水线和分水线上的高程点，对原始高精度集水线网络和分水线网络进行相应尺度的化简； 步骤6，根据步骤4和步骤5所得到的相应尺度下的地形特征点集和化简的水文系统，建立受约束的不规则三角网TIN，在此基础上内插规则格网作为综合后的DEM ； 步骤7，计算综合后DEM的精度，若精度达不到用户设定尺度的综合程度，采用弹性寻位策略自动调整综合程度控制参数，获取相应尺度的地形特征点子集，重复步骤步骤7。
2.根据权利要求1所述的一种顾及地表水文上下文的多尺度DEM建模方法，其特征在于，所述步骤I中，建立原始高精度DHM的水文网络的具体步骤如下: 步骤1.1，根据DEM规则格网的高程差异建立流向矩阵，以未标记的任意位置栅格作为种子点，按八邻域逆向迭代搜索汇入栅格，以只有流出栅格没有流入栅格的点作为径流源头点，然后按照径流方向追踪，采用改进D8算法建立径流累积矩阵，具体是:首先遍历规则格网，计算每个格网向八领域格网的高程差，把高程差最大的方向标识为水流方向，得到流向矩阵；然后随机选取未标记的格网作为种子点，按八邻域逆向搜索汇入栅格，以只有流出栅格没有流入栅格的点作为径流源头点，再按径流方向追踪，水流路径上某一栅格的径流累积量为上游栅格的径流累积量加I ;对于多个方向汇入的栅格，径流累积量在各汇入方向累积量逻辑和的基础上加I，得到径流累积矩阵； 步骤1.2，根据步骤1.1所得径流累积矩阵提取每一个格网点的径流累积量，提取大于设定阈值的稀疏栅格作为集水线栅格，结合步骤1.1中建立的流向矩阵从集水线栅格中提取集水线的源头，然后再根据步骤1.1中建立的径流累积矩阵以集水线的源头作为种子点，逆向延伸集水线直至径流源头； 步骤1.3，对步骤1.2建立的集水线稀疏栅格进行矢量化，然后利用步骤1.2建立的流向矩阵构建集水线网络； 步骤1.4，提取以集水线的源头与叉头或者叉头与出口之间的线段作为流域的集水核，然后根据步骤1.1得到的流向矩阵和径流累积矩阵将所有流入集水核的区域标识为对应的汇水响应区域；步骤1.5，对步骤1.4得到的汇水区域边界进行锐化和骨架线提取，建立分水线网。
3.根据权利要求1所述的一种顾及地表水文上下文的多尺度DEM建模方法，其特征在于，所述步骤2中，解析原始高精度DEM所描述地貌对象的水文上下文的具体步骤如下: 步骤2.1，以步骤1.4获得的集水核作为集水系统的基本数据组织单元，利用地表汇水响应区域的汇入关系建立集水网络的层次结构树，即有: 父子关系的集水核，父级集水核的汇水响应区域包含所有子级集水核的汇水响应区域； 以及兄弟关系的集水核，互为兄弟的集水核其汇水响应区域汇入同一父级汇水响应区域； 步骤2.2，将分水线的起点与叉点或叉点与出口之间的线段作为流域的分水核，基于分水线与流域的耦合关系，利用汇水响应区域的汇入关系建立分水核网络的层次结构树，即有: 父子关系的分水核，父级分水核与父级汇水响应区域的边界重合，子集分水核与子集汇水响应区域的边界重合且与父级汇水响应区域的边界相交； 兄弟关系的分水核，互为兄弟关系的分水核其对应的汇水响应区域同为兄弟关系，且均不与父级汇水响应区域的边界重合。
4.根据权利要求1所述的一种顾及地表水文上下文的多尺度DEM建模方法，其特征在于，步骤3对原始高精度DEM的地表高程点进行水文语义增强，具体实现方式是根据步骤2得到的集水线网络与分水线网络的层次结构建立流域指数函数，然后通过计算流域指数和微地貌指数两种水文贡献指数对地表高程点的水文特征进行泛化。
5.根据权利要求1所述的一种顾及地表水文上下文的多尺度DEM建模方法，其特征在于，所述步骤4中，获取地形特征点集的具体实现步骤如下: 步骤4.1，将原始高精度DEM按流域划分为若干个特征提取单元，提取流域边界凸壳构建初始不规则三角网TIN，获取除边界外的剩余高程点作为候选特征点集； 步骤4.2，随机抽取候选特征点集中的一个候选特征点，结合步骤3得到的地表高程点水文贡献指数计算出其到不规则TIN表面的加权点面距，选取加权点面距最大的高程点作为当前第I特征点，加入地形特征点序列中； 步骤4.3，将该地形特征点插入不规则TIN中进行重构，从候选特征点集中删除该特征点，重复步骤4.2 步骤4.3，直至加权点面距小于给定的阈值或者所有的高程点都进入特征点序列。
6.根据权利要求1所述的一种顾及地表水文上下文的多尺度DEM建模方法，其特征在于，所述步骤5中，化简水文系统的具体实现方式是，根据用户指定尺度确定地形特征点的选取数量，从步骤4得到的特征点序列中获取相应子集，然后从地形特征点子集中分离出位于集水线和分水线上的地形特征点，分别沿着初始集水线和分水线路径，根据子集中的集水线和分水线地形特征点重构新的集水线和分水线，重构过程中保持集水线和分水线与一般地形特征点的拓扑关系不变。
7.根据权利要求1所述的一种顾及地表水文上下文的多尺度DEM建模方法，其特征在于，所述步骤7中，采用弹性寻位策略自动调整综合程度控制参数是通过计算步骤6所得的综合后DEM的精度来进行，具体实现步骤如下:步骤7.1，根据多尺度DEM资源库建立尺度与DEM高程描述差的关系函数F_s，根据用户设定的尺度获取相应的F_s(x)，比较综合前后的DEM高程描述误差，计算出高程描述误差的差值F (X)； 步骤7.2，采用弹性寻位策略基于以下约束条件控制DEM综合程度:设定容差阈值为ε ,若I F_s (X)-F (X) |> ε且F_s (x) >F (χ),则增加从步骤4特征点序列中提取的特征点集，反之，若I F_s (χ) -F(x) I > ε且F_s (χ)〈F(χ),则减少从步骤4特征点序列中提取的特征点个数，重复步骤Γ步骤7， 直至|F_s(X)-F(X)|〈e结束。
全文摘要
本发明属于数字高程模型建模领域，特别涉及一种顾及地表水文上下文的多尺度DEM建模方法。通过对DEM进行径流分析提取地形表面的水文响应区域，建立汇水区域的层次结构树；通过泛化地表高程点的水文特征，计算高程点对于地貌形态的水文贡献指数；通过考察高程点的加权地形描述误差选取地形特征点集；并采用弹性寻位策略自动调整综合程度控制参数，从而实现顾及地表水文上下文的多尺度DEM自动综合。因此，本发明具有如下优点可有效保持原始高精度DEM的地形和水文特征，并能保证多尺度下DEM的地形与水文特征的拓扑和逻辑一致性；提高了多尺度DEM建模的自动化和生成效率，并能满足实际应用中对尺度自适应DEM的需要。
文档编号G06T17/05GK103236086SQ20131014494
公开日2013年8月7日 申请日期2013年4月24日 优先权日2013年4月24日
发明者黄丽娜, 郑斌, 费立凡, 李跃 申请人:武汉大学