一种表现力强的地貌晕渲图洼谷地形填补方法与流程

本发明涉及一种地貌晕渲图洼谷地形填补方法，特别涉及一种表现力强的地貌晕渲图洼谷地形填补方法，属于地貌晕渲图制作技术领域。

背景技术：

地貌是地球表面各形态的总称，有多种空间实物形体，分为自然形体和人工形体两大类。地球表面自形成以来一直在不断的变化，地球表面各种地貌形态是内力和外力共同作用的结果，其中内力因素决定地貌主要形态，外力因素的侵蚀作用形成地貌的次级形态。地球表面面积巨大，地貌形态千姿百态，在不同比例尺的地图中合理正确的表示不同的地貌形态，一直是地图制图的重点难点。地图中地貌表现方法主要有晕滃法、地貌写景法、等高线法、分层设色法、地貌晕渲法，地貌晕渲法根据光照原理，采用明暗和冷暖对比表现地形起伏，该方法表达直观生动，能够容易的建立三维立体感。相对于其它地貌表示方法优势明显，是平面地图中表示地貌形态的最主要方法，得到了广泛的研究和应用。

在手工制图时代，制作晕渲图是一项复杂的制图过程，需要有丰富的专业制图知识，熟悉各种地貌类型，有很好的三维立体感，同时需要制图人员有相当的绘图和色彩运用能力。手工晕渲制图不仅非常依赖制图人员的制图水平和经验，而且是一个非常复杂艰苦的过程，一张手工晕渲图需要耗费大量的精力，这是手工制作晕渲图的巨大瓶颈。

当前计算机技术飞速发展，其在地图制图领域里也得到广泛的应用，数字地图制图技术彻底变换了传统的地图制图的方式。作为数字地貌晕渲图的基本信息源，数字高程模型获取方式不断丰富，数据精度不断提高，这也促进了数字地貌晕渲技术的发展。数字地貌晕渲技术有许多优势：一是成图速度快，大大缩短了数字制图周期；二是制图成本低，数字制图成图周期降低，人工成本减少，制图成本直线下降；三是受制图人员主观因素影响小。数字晕渲已是制作晕渲图的主要方法，很多制图软件如arcgis等都能制作数字地貌晕渲图。

但数字地貌晕渲制图也有明显不足之处，地貌晕渲没有数量概念，只表现地貌高低起伏，数字地貌晕渲是以数字高程模型为基本信息源，但高精度的数字高程模型数据包含太多的地形细节信息，这些零碎的地形细节信息干扰了数字晕渲图中主体地貌特征的表达，淡化了晕渲图的一览性，降低了地貌晕渲图的表现力。

为制作直观和表现力强的数字地貌晕渲图，需要对地貌进行科学合理的综合，在保持或夸张主体地貌特征基础上，综合不必要的地形细节信息，突出地貌的主体特征。因此，地貌综合是数字地貌晕渲图制作必不可少的环节和重要的手段，地貌综合的好坏直接关系到地貌晕渲图的质量高低，特别是中小比例尺的地貌晕渲图更为明显，所以地貌综合是地貌晕渲技术中关键问题。

现有技术中，数字地貌晕渲中一张晕渲图的光照方向一般都是固定不变，为了改善数字地貌晕渲的效果，采用模拟手工地貌晕渲中局部调整光照的方法，同时为了在地貌晕渲图中突出一些关键地形特征，计算地形单元受光量是考量平面曲率等地形因素，但这些方法效果都不够理想。以数字高程模型为基础信息源制作数字地貌晕渲图时，晕渲图中的主体地貌特征往往受到数字高程模型过度精细的地形细节干扰。

数字地貌晕渲技术是以数字高程模型为基础信息源，在制作中小比例尺地貌晕渲图时，零碎的地形细节阻碍了关键地形特征的表达，需要对数字高程模型进行综合，以删减不必要的地形细节信息，而且随着数字高程模型数据获取的精度的提高，数字地貌晕渲中进行地貌综合显得尤为必要，现有技术出现了空间滤波算法对数字高程模型进行综合，并且尽量避免关键地形特征被过渡平滑。对于制作小比例尺晕渲图可以采用重采样方法综合数字高程模型，同时对数字高程模型进行局部的地形夸张，通过分辨率互补的方法融合实验区域高低分辨率的数字高程模型数据，但对于不同地貌形态综合效果不一,特别是陡坡和大面积平地效果不好。

综合来看，针对现有技术存在的一些缺陷，本发明拟解决以下问题：一是在制作数字晕渲图特别是中小比例尺的晕渲图中，数字高程模型中包含的过度精细的地形细节会降低晕渲图的一览性与三维立体感，所以需要从数字高程模型中删减这些不必要的地形细节信息，数字地貌晕渲中往往采用地貌综合的方法删减数字高程模型中过多的地形细节，但是这也不可避免地导致地形的主要特征被过渡平滑与削弱，针对综合掉数字高程模型中零碎的地形细节的同时保持地形的主要特征的难题；二是删减次级地形细节与保持主体地形特征的矛盾问题，数字地貌晕渲是以数字高程模型为基础数据信息源,晕渲图的主体地貌特征的描述往往被数字高程模型所包含的过度零碎的地形细节所干扰,这些零碎的地形细节信息干扰了数字晕渲图中主体地貌特征的描述，特别是对于中小比例尺的晕渲图的影响尤为明显，往往导致关键地形特征被综合漏掉；三是现有技术洼谷地形填补的数字高程模型综合方法属于结构线综合方法，该方法填补待删减的次级洼谷地形，达到删减地形细节综合数字高程模型的目的，该方法先确定洼谷地形区域与边界，计算洼谷地形的大致走向，最后根据洼谷地形线的法线方向延伸到洼谷地形边界，在洼谷地形边界两点之间进行线性填补。这样填补的洼谷地形等高线沿着洼谷地形线的法线方向延伸，比较平直，与洼谷地形周边的等高线不协调，无法使等高线平滑过渡，地貌晕渲图的整体效果不够好。

技术实现要素：

针对现有技术的不足，本发明提供的一种表现力强的地貌晕渲图洼谷地形填补方法，在删减地形细节信息时能较好的保持地形特征，该方法构建数字高程模型金字塔，通过填补不同级别的洼谷地形生成不同综合程度的数字高程模型图层，调整各个数字高程模型图层的权重系数来突出或则削弱某一图层，在删减地貌碎部的同时较好的保持地形特征，而且还能通过权重系数的设置突出地形的关键地形特征，制作出一览性好重点突出且三维立体感强的数字地貌晕渲图。填补洼谷地形时考量周围等高线在洼谷地形内走势，并分段插值洼谷地形内网格单元高程值，使得填补后的洼谷地形内地形与周围地形平滑过渡，数字地貌晕渲图的洼谷地形表达更加准确自然，数字地貌晕渲图更加清楚美观。

为达到以上技术效果，本发明所采用的技术方案如下：

一种表现力强的地貌晕渲图洼谷地形填补方法，采用洼谷地形填补方法综合数字高程模型，填补次级洼谷地形删减地形细节信息，依次填补不同级别的洼谷地形生成不同综合程度的数字高程模型图层，即数字高程模型金字塔，并为各层数字高程模型图层设置相应的权重系数，通过调整权重系数突出重建地形中的关键地形特征；在删减过度零碎的地形细节信息时避免关键地形特征被淡化，通过数字高程模型图层权重系数的调整适当夸张关键地形特征；填补洼谷地形时考量周围等高线的趋势，采用三次多项式插值法填补待删减的洼谷地形；

一种表现力强的地貌晕渲图洼谷地形填补方法的具体步骤为：

第一步，提取洼谷地形线：采用物理流水模拟法提取集水线；

第二步，洼谷地形分级编码：基于霍顿分级的编码原则，编码唯一且包含洼谷地形重要信息；

第三步，洼谷地形识别：基于洼谷地形分级编码和网格单元的流向确定洼谷地形的范围和洼谷地形边界；

第四步，考量填补洼谷地形周围等高线趋势：根据等高线与洼谷地形边界的交点和在交点处等高线切线模拟等高线的走势；

第五步，三次多项式插值：根据等高线与洼谷地形边界交点和在交点处的等高线的切线方向进行三次多项式插值填补洼谷地形；

第六步，插值等高线间插值：采用线性插值法变换插值新生成的等高线间的洼谷地形单元；

第七步，构建数字高程模型金字塔：填补不同级别洼谷地形获得不同详细程度的数字高程模型图层，通过各图层的权重系数突出或削弱某一级别的地形特征。

一种表现力强的地貌晕渲图洼谷地形填补方法，进一步的，第一步，提取洼谷地形线时，物理流水模拟法提取集水线的过程分为五个步骤：

第1步，网格单元的流向确定；

第2步，网格数字高程模型数据中的洼地和平地处理；

第3步，网格单元汇水量计算；

第4步，汇水流量阈值确定；

第5步，集水点连接。

一种表现力强的地貌晕渲图洼谷地形填补方法，进一步的，本发明提出基于霍顿分级的洼谷地形编码方法,洼谷地形根据霍顿定律分级,每条洼谷地形都对应一个唯一的洼谷地形编码，且包含洼谷地形间的连接关系、洼谷地形级别和洼谷地形在洼谷地形系统中编号信息。

一种表现力强的地貌晕渲图洼谷地形填补方法，进一步的，洼谷地形编码后，各条洼谷地形都有一个唯一且包含了洼谷地形信息的洼谷地形编码，根据洼谷地形编码进行洼谷地形集水区域识别和洼谷地形边界的识别；

洼谷地形集水区域识别中，流向洼谷地形所在单元的网格单元属于当前洼谷地形的集水区域；依次遍历洼谷地形数据的存储结构，获得洼谷地形的编码和点集，查询洼谷地形的点集，找到集水量最大网格单元，洼谷地形集水量最大的单元即洼谷地形的出水口；如果洼谷地形集水量最大的点不是连接点，则把该集水量最大的网格单元压入一个先进先出的堆栈，弹出堆栈的堆首单元，并以该单元为中心建立3*3的滑动窗口，把其中流向中心单元的网格单元设置当前洼谷地形编码并压入堆栈，表明该单元为当前洼谷地形的集水区域，循环直至先进先出的堆栈为空结束，数字高程模型的网格单元中洼谷地形编码为当前洼谷地形编码的区域就是当前洼谷地形的集水区域；当洼谷地形集水量最大的点是一个连接点即与父级洼谷地形相交时，从洼谷地形的点集中找出集水量第二高的点，以这个点做为洼谷地形的出水口，进行递归迭代计算出当前洼谷地形的集水区域；

洼谷地形边界识别中，建立一个3*3的滑动窗口，使滑动窗口滑动遍历数字高程模型所有的网格单元，如果中心单元的洼谷地形编码不为空且为洼谷地形集水区域的边界，则以该单元为起始单元顺时针沿洼谷地形集水区域的边界遍历，沿着固定方向把当前洼谷地形集水区域的边界单元依次加入到边界数组中；接下来用一个3*3滑动窗口从头开始遍历洼谷地形集水区域边界数组，如果窗口中存在二个以上的相邻边界单元，则比较这几个边界相邻的单元，如果其中单元有重复，只保留其中的一个边界单元，其余的边界单元类型从边界改为洼谷地形集水区域内部单元。

一种表现力强的地貌晕渲图洼谷地形填补方法，进一步的，本发明在填补洼谷地形时考量周围登高线的趋势，根据填补洼谷地形的级别提取相应等高距的等高线，计算等高线与填补洼谷地形边界的交点；采用等高线与洼谷地形边界交点处的斜率表示等高线的趋势，根据同一条等高线与洼谷地形边界的二个交点和二个交点处的斜率进行三次多项式插值，将交点间的插值单元赋值等高线的高程；填补洼谷地形内高程没有修改的单元沿着洼谷地形的大致走势进行线性内插，直到所有填补洼谷地形内的单元高程都被修改为止。

一种表现力强的地貌晕渲图洼谷地形填补方法，进一步的，本发明提出根据等高线与洼谷地形边界交点和在交点处的等高线的切线方向进行三次多项式插值填补洼谷地形的方法；

三次多项式插值假设平面上有n个已知点(xi，yi)(i＝1,2…n)，要曲线拟合出一条光滑的曲线y＝f(x)，使得yi＝f(xi)(i＝1,2…n)。根据已知点三次多项式插值出的曲线需要满足二个要求：一是yi＝f(xi)(i＝1,2…n)，即曲线通过所有的已知点；二是曲线光滑，直观上线条不能有明显的折角，数学上是曲线方程任意一点都存在一阶导数yi’＝f’(xi)。三次多项式插值是用三次多项式无限逼近相邻两点，设有二个相邻的点di和di+1，在二个点之间进行三次多项式插值，其多项式曲线方程：

x＝d0+d1z+d2z²+d3z³

y＝e0+e1z+e2z²+e3z³

式中di、ei(i＝0、1、2、3)为常数，z为可变参数，当点di沿着三次多项式曲线变化到点di+1时，参数z是从0变化到1；

通过确定等高线与待填补洼谷地形边界的交点，计算在交点处等高线的近似切线与x轴的夹角；计算di和ei(i＝0、1、2、3)，使参数z从0变化到1进行点di和点di+1间的三次多项式插值；

本发明提出在计算插值之前先计算二个洼谷地形边界交点之间的距离的网格单元数量，记为f，参数z从0变化到1的插值间隔dz＝1/2f，插值结束后获得二个洼谷地形边界交点之间的插值数组，数组存储二个洼谷地形边界交点之间连接的洼谷地形单元，把洼谷地形内的插值单元高程值都赋值为等高线的高程值。

一种表现力强的地貌晕渲图洼谷地形填补方法，进一步的，填补洼谷地形内没有变换过的网格单元线性内插前，确定洼谷地形的粗略走势，查看洼谷地形线的上游点是否到达洼谷地形源头，如果洼谷地形线上游点没有到达洼谷地形源头，即上游点的集水量不等于0，需要进行洼谷地形线的延长；洼谷地形线的延长方法为：

步骤一，所有流向洼谷地形上游点的相邻网格单元为候选上游点；

步骤二，从候选的网格单元中选出集水流量最大的点为新的洼谷地形上游点，如果集水流量最大的网格单元不止一个，则选择最邻近的网格单元为新的洼谷地形上游点；

重复上述步骤一和步骤二，直到上游点的集水流量等于0，则洼谷地形线的延长工作结束；

本发明采用洼谷地形线各单元流向分段确定洼谷地形粗略走势的方法，分别统计新插值生成等高线间的洼谷地形线单元的流向个数，数量最多的流向为插值生成的等高线间洼谷地形的粗略走势。

一种表现力强的地貌晕渲图洼谷地形填补方法，进一步的，通过三次多项式插值法生成待填补洼谷地形新的等高线和确定等高线间洼谷地形的粗略走势后，进行待填补洼谷地形内没有变换的洼谷地形单元的线性填补；具体方法为：依次遍历所有的网格单元，如果该单元为待填补洼谷地形内，且高程值没有变换，则向插值新生成的等高线间洼谷地形的粗略走势的方向两端扩充，直至遇到当前洼谷地形的边界或插值生成等高线的单元结束，获得一个首尾洼谷地形边界单元或新生成等高线单元的数组，计算数组两端的高程差，并对数组内其它单元进行线性内插，使得首尾两端的单元平滑过渡。

一种表现力强的地貌晕渲图洼谷地形填补方法，进一步的，采用基于洼谷地形填补的方法进行地貌综合，综合过滤地貌晕渲中过度零碎的地形细节信息，采用不同综合程度生成不同详细程度的数字高程模型图层，不同详细程度的图层构成一个地形详细程度递减的数字高程模型金字塔，通过设置各图层的权重系数，调整相应精细程度的图层在重建地形中的比重，突出关键地形特征；

本发明提出填补某一级别所有的洼谷地形生成相应详细程度的数字高程模型图层，假设提取的洼谷地形有m个级别，填补第i级洼谷地形生成数字高程模型图层ai(i＝1、2…m-1)，当填补的洼谷地形级别i大于1时，高级别洼谷地形内经常包含有与其相连的次级洼谷地形，在洼谷地形识别和确定洼谷地形边界时要做特别处理；高级别洼谷地形内的与其相连的低级别洼谷地形属于当前高级别洼谷地形，在洼谷地形识别时根据洼谷地形编码确定高级洼谷地形和低级洼谷地形的关系，本发明提出的洼谷地形编码原则是与其相连的父级洼谷地形编码加上本身洼谷地形级别和洼谷地形编号，只要低级别洼谷地形编码的头几位与高级别洼谷地形编码相同则确定低级别洼谷地形属于当前高级别洼谷地形；依次填补不同级别洼谷地形生成不同详细程度的数字高程模型图层ai(i＝0,1,2,3…m)，其中a0表示没有任何综合的原始地形图层，各个数字高程模型图层设置一个权重系数vi(i＝0，1,2,3…n)，重建的地形g公式如下所示：

各图层的权重系数根据综合尺度和不同地形特征调整，若综合的尺度较小，则低级别图层的权重系数相应大一些；若要突出地形的主要特征，则设置高级别图层的权重系数大一些。

与现有技术相比，本发明的优点和创新点在于：

一是本发明提供的一种表现力强的地貌晕渲图洼谷地形填补方法，从结构线综合数字高程模型的方法出发，采用填补洼谷地形的方法综合数字高程模型，生成不同精细程度的数字高程模型图层，通过设置不同详细程度层次的数字高程模型的重要系数来突出地形的关键地形特征。同时，为了改善填补洼谷地形综合数字高程模型方法的地貌综合效果，本发明提出模拟手工地貌综合的方法填补待删减的洼谷地形，保持地形的主体特征同时又删减不必要地形细节，地貌晕渲图的整体效果好。

二是本发明提供的一种表现力强的地貌晕渲图洼谷地形填补方法，基于洼谷地形填补数字高程模型综合方法属于结构化综合方法，在删减地形细节信息时能较好的保持地形特征。该方法构建数字高程模型金字塔，通过填补不同级别的洼谷地形生成不同综合程度的数字高程模型图层，调整各个数字高程模型图层的权重系数来突出或则削弱某一图层，在删减地貌碎部的同时较好的保持地形特征，在综合数字高程模型的同时较好地保持地形的主体特征不被过渡淡化，而且还能通过权重系数的设置突出地形的关键地形特征，制作出一览性好重点突出且三维立体感强的数字地貌晕渲图。

三是本发明提供的一种表现力强的地貌晕渲图洼谷地形填补方法，线性填补的洼谷地形填补方法根据洼谷地形走势的法线方向进行洼谷地形填补，该填补方法简单易行，本发明填补洼谷地形时考量周围等高线在洼谷地形内走势，并分段插值洼谷地形内网格单元高程值，用三次二项式插值代替线性插值，计算出平滑的插值结果，从而使洼谷地形边界之间平滑过渡，数字地貌晕渲图的洼谷地形表达更加准确自然，数字地貌晕渲图更加清楚美观。

附图说明

图1是本发明的流向图式与流向编码示意图。

图2是本发明网格单元汇水量计算过程示意图。

图3是本发明洼谷地形的分级编码示意图。

图4是本发明三次多项式插值后生成的阶梯示意图。

图5是本发明分段统计流向确定洼谷地形粗略走势方法示意图。

图6是本发明等高线间的线性填补示意图。

具体实施方式

下面结合附图，对本发明提供的一种表现力强的地貌晕渲图洼谷地形填补方法的技术方案进行进一步的描述，使本领域的技术人员可以更好的理解本发明并能予以实施。

本发明提供的一种表现力强的地貌晕渲图洼谷地形填补方法，采用洼谷地形填补方法综合数字高程模型，填补次级洼谷地形删减地形细节信息，依次填补不同级别的洼谷地形生成不同综合程度的数字高程模型图层，即数字高程模型金字塔，并为各层数字高程模型图层设置相应的权重系数，通过调整权重系数突出重建地形中的关键地形特征。本发明提出的方法在删减过度零碎的地形细节信息时避免了关键地形特征被淡化，而且通过数字高程模型图层权重系数的调整可以适当夸张关键地形特征。同时，为了使填补后洼谷地形内等高线平滑过渡，即填补后谷地内地形与周围地形协调过渡，本发明填补洼谷地形时考量周围等高线的趋势，并采用三次多项式插值法填补待删减的洼谷地形；

本发明提出的一种表现力强的地貌晕渲图洼谷地形填补方法的具体步骤为：

第一步，提取洼谷地形线：采用物理流水模拟法提取集水线；

第二步，洼谷地形分级编码：基于霍顿分级的编码原则，编码唯一且包含洼谷地形重要信息；

第三步，洼谷地形识别：基于洼谷地形分级编码和网格单元的流向确定洼谷地形的范围和洼谷地形边界；

第四步，考量填补洼谷地形周围等高线趋势：根据等高线与洼谷地形边界的交点和在交点处等高线切线模拟等高线的走势；

第五步，三次多项式插值：为了使填补的洼谷地形内等高线平滑过渡，根据等高线与洼谷地形边界交点和在交点处的等高线的切线方向进行三次多项式插值填补洼谷地形；

第六步，插值等高线间插值：采用线性插值法变换插值新生成的等高线间的洼谷地形单元；

第七步，构建数字高程模型金字塔：填补不同级别洼谷地形获得不同详细程度的数字高程模型图层，通过各图层的权重系数突出或削弱某一级别的地形特征。

一、地形特征线提取

本发明提出的地貌晕渲图洼谷地形填补方法，首先采用物理流水模拟法提取集水线，并对集水线进行分级编码，本发明提出的编码方法有利于后续洼谷地形的自动识别与填补。

(一)物理流水模拟法

物理流水模拟法提取集水线的过程分为五个步骤：

第1步，网格单元的流向确定；

第2步，网格数字高程模型数据中的洼地和平地处理；

第3步，网格单元汇水量计算；

第4步，汇水流量阈值确定；

第5步，集水点连接。

(1)网格单元的流向确定

物理流水模拟法采用d8算法确定每个网格单元的流向，d8算法中的流向图式和流向编码如图1所示。d8算法中数字高程模型的每个网格单元的降水量即初始流量相同，并设置为1，建立一个3*3滑动窗口，依次遍历所有的数字高程模型网格单元，在与中心单元zd相邻的八个单元中找出高程最小的单元xd，如果xd的高程小于zd的高程则根据xd的位置给zd赋以流向编码，如果xd的高程大于等于zd的高程，则给zd的流向编码为0，表示该单元的水不流向其它单元。

(2)网格数字高程模型数据中的洼地和平地处理

地表上的平地或洼地在降雨时会形成集水区，采用d8算法确定数字高程模型网格单元的流向时，如果高程最小的单元xd的高程大于等于中心单元zd的高程，则说明xd是平地或洼地单元，需要对数字高程模型中的洼地和平地单元进行处理以确保各个网格单元都有流出的流向。采用d8算法计算网格单元后会形成一个流向矩阵，其中流向编码为0且相邻的单元为同一个平地或洼地，确定平地或洼地的边界，并找出进水口和出水口单元，根据进水口和出水口单元确定平地或洼地内单元的流向，填平洼地确定洼地内流向。

(3)网格单元汇水量计算

经过上述两步后，可得到网格单元的流向编码矩阵，且其中的所有网格单元都有流出的流向编码，接下来根据流向编码矩阵计算每个网格单元的汇水量。单元汇水量等于自身汇水量和流入单元汇水量之和，如果流入单元汇水量还没有计算，则先计算流入单元汇水量，直到所有的单元汇水量都确定，汇水量计算过程如图2所示。

(4)汇水流量阈值确定

汇水流量计算完后会得到格网单元的汇水流量矩阵，需要根据汇水量确定集水点，但如果汇水流量阈值设置过大，则提取的集水线不完整；如果汇水流量阈值过小，则提取的集水线过于详细。因为河流水系的自相似性和平均分叉比趋近于常数，根据汇水流量阈值与提取的河流水系分叉比的关系确定的阈值方法，可以避免阈值过大过小，能较好的保持洼谷地形结构特征。

(5)集水点连接

根据汇水流量阈值确定洼谷地形集水点后，把这些离散的洼谷地形集水点连接成洼谷地形线。从没有连接过得洼谷地形集水点中选择汇水量最小的开始连接，根据洼谷地形流向确定下一个洼谷地形集水点，直到与已经连接的洼谷地形线相交或者没有下一个连接点，循环连接直到所有的洼谷地形集水点都被连接。

二、洼谷地形分级编码

采用物理流水模拟法提取完洼谷地形线之后进行洼谷地形的分级编码,洼谷地形根据霍顿定律分级,本发明提出基于霍顿分级的洼谷地形编码方法,该洼谷地形编码唯一,且隐含主支流之间的连接关系和支流的数量信息,是洼谷地形自动识别和填平的基础。

(一)洼谷地形的霍顿分级

水系结构具有明显的自相似性，数量和结构关系可用霍顿定律进行描述。霍顿编码特点有：一是对河系洼谷地形的地理实体编码，二是反应了河系洼谷地形的等级关系，三是反应河系洼谷地形的深度。

实施例对洼谷地形进行霍顿分级编码，使用arcgis中的“水文分析”工具提取集水线，编码方法只有strahler和shreve编码方法，进行strahler编码到霍顿编码的转换。编码转换先从strahler编码中级别最高且没有转换过的洼谷地形开始搜索，搜索过程中遇到连接点时首先判断连接洼谷地形的级别，如果连接洼谷地形的级别不相等，则把级别高的洼谷地形并入当前搜索洼谷地形。

(二)洼谷地形的分级编码

本发明基于霍顿分级编码提出了新的洼谷地形分级编码方法，该洼谷地形分级编码方法在洼谷地形系统中唯一，且包含了洼谷地形间的连接关系、洼谷地形级别和洼谷地形在洼谷地形系统中编号等重要信息。每条洼谷地形都对应一个唯一的洼谷地形编码，通过洼谷地形编码解译出很多关于洼谷地形的重要信息，为后续的计算机自动识别洼谷地形汇水区域、确定洼谷地形边界和填补洼谷地形相关操作提供了极大的方便。

在进行洼谷地形编码前先要遍历所有洼谷地形，统计出不同级别洼谷地形a1、a2…an的数量{b1、b2…bn}，根据相应级别的洼谷地形数量确定表示洼谷地形编号的长度，如10<b1<100,则第一级洼谷地形的编号用两位数字表示，依次类推。本发明把与低级别洼谷地形相连的高级别洼谷地形称为低级别洼谷地形的父级洼谷地形，洼谷地形之间是父子关系。洼谷地形编码是由父级洼谷地形编码和当前子洼谷地形编码组成，如果没有父级洼谷地形就是只由当前洼谷地形编码组成，图3为洼谷地形的分级编码示意图。

三、洼谷地形识别

洼谷地形编码后，各条洼谷地形都有一个唯一且包含了洼谷地形信息的洼谷地形编码，根据洼谷地形编码进行洼谷地形集水区域和洼谷地形边界的识别。

(一)洼谷地形集水区域识别

洼谷地形集水区域是一个综合性的指标，一定程度上是洼谷分布密度、地形等级、长度和宽度多个因素的集成。本发明根据上述步骤确定的洼谷地形编码和物理流水模拟法计算的网格单元的流向确定洼谷地形的集水区域，自然界中洼谷地形有集水作用，即洼谷地形两侧山坡的水流向谷底汇聚在一起流出，所以本发明中流向洼谷地形所在单元的网格单元属于当前洼谷地形的集水区域。

依次遍历洼谷地形数据的存储结构，获得洼谷地形的编码和点集，查询洼谷地形的点集，找到集水量最大网格单元，从洼谷地形地貌形态可知，洼谷地形出口的集水量等于整个集水区域的集水量，所以洼谷地形集水量最大的单元即洼谷地形的出水口。如果该洼谷地形集水量最大的点不是连接点，则把该集水量最大的网格单元压入一个先进先出的堆栈，弹出堆栈的堆首单元，并以该单元为中心建立3*3的滑动窗口，把其中流向中心单元的网格单元设置当前洼谷地形编码并压入堆栈，表明该单元为当前洼谷地形的集水区域，循环直至先进先出的堆栈为空结束，这时在数字高程模型的网格单元中洼谷地形编码为当前洼谷地形编码的区域就是当前洼谷地形的集水区域。

当洼谷地形集水量最大的点是一个连接点即与父级洼谷地形相交时，不能简单的以该点做起始单元通过流向判断洼谷地形的集水区域，因为这个单元不仅是当前洼谷地形出水口，也是与当前洼谷地形相连的父级洼谷地形中一个洼谷地形单元，流向这个单元的不仅包括当前洼谷地形的集水区域内的网格单元，还包括，父级洼谷地形中这个单元上游所有集水区域内的网格单元，如果以这个单元计算当前洼谷地形的集水区域出现图中的错误现象。所以从洼谷地形的点集中找出集水量第二高的点，以这个点做为洼谷地形的出水口，进行递归迭代计算出当前洼谷地形的集水区域。

(二)洼谷地形边界识别

上述步骤根据洼谷地形编码和网格单元的流向确定了每个洼谷地形的集水区域，并给洼谷地形集水区域内所有的网格单元赋以洼谷地形编码标识。接下来确定洼谷地形集水区域的边界，洼谷地形集水区域的边界标明了流入当前洼谷地形的范围，这个范围的网格单元流入其它洼谷地形。洼谷地形集水区域的边界也称为分水岭，分水岭包含很多重要的特征点，对后续洼谷地形填补十分重要。

建立一个3*3的滑动窗口，使滑动窗口滑动遍历数字高程模型所有的网格单元，如果中心单元的洼谷地形编码不为空，则查看相邻的8个网格单元，若其中有洼谷地形编码为空的单元或与中心单元编码不同的网格单元，则该中心单元为洼谷地形集水区域的边界，因为它与不是本洼谷地形集水区域的单元相邻；如果中心单元的编码为空，则说明该单元不属于任何一个洼谷地形集水区域，继续遍历，直到所有的单元都遍历完。按照上述方法只需要遍历一次数字高程模型的网格单元就可以很快确定洼谷地形集水区域的边界，但是这样确定的洼谷地形集水区域的边界有一些不必要的重复单元，考量到洼谷地形集水区域边界单元可以删减，所以对这些洼谷地形集水区域边界单元进行细化处理。

建立一个3*3滑动窗口，遍历数字高程模型的所有网格单元，如果中心单元的洼谷地形编码不为空且为洼谷地形集水区域的边界，则以该单元为起始单元顺时针沿洼谷地形集水区域的边界遍历，沿着固定方向把当前洼谷地形集水区域的边界单元依次加入到边界数组中。接下来用一个3*3滑动窗口从头开始遍历洼谷地形集水区域边界数组，如果窗口中存在二个以上的相邻边界单元，则比较这几个边界相邻的单元，如果其中单元有重复，只保留其中的一个边界单元，其余的边界单元类型从边界改为洼谷地形集水区域内部单元。

四、次级洼谷地形的填补

经过上述步骤，采用物理流水模拟法提取了集水线，并对集水线进行了霍顿分级和编码，基于本发明提出的编码方法对洼谷地形集水区域范围的识别和边界的确定，接下来进行次级洼谷地形的填补。

现有技术的洼谷地形填补方法是先确定要填补洼谷地形的范围和边界，延长洼谷地形线并统计出洼谷地形的大致走势，在洼谷地形大致走势的法线方向上的洼谷地形边界的二个点之间洼谷地形内的点进行高程值的线性内插。这样的洼谷地形填补方法简单，但是这种沿着洼谷地形大致走势法线方向线性填补的方式，会造成填补的等高线沿着洼谷地形大致走势法线方向平直延伸，与周边等高线不协调。针对填补后等高线过度平直，与周围等高线不协调的问题，本发明做出相应改进。

本发明在填补洼谷地形时考量了周围登高线的趋势，根据填补洼谷地形的级别提取相应等高距的等高线，计算等高线与填补洼谷地形边界的交点，实施例为了计算方便采用等高线与洼谷地形边界交点处的斜率表示等高线的趋势。根据同一条等高线与洼谷地形边界的二个交点和二个交点处的斜率进行三次多项式插值，该插值法可以使二个交点处圆滑过渡，将交点间的插值单元赋值等高线的高程。填补洼谷地形内高程没有修改的单元沿着洼谷地形的大致走势进行线性内插，直到所有填补洼谷地形内的单元高程都被修改为止。

(一)确定洼谷地形边界与等高线交点

在确定填补洼谷地形边界与等高线交点前，需要基于网格数字高程模型提取等高线，提取的等高线的等高距要根据填补洼谷地形来确定，如果填补的洼谷地形级别较低，且洼谷地形多是细小狭窄的洼谷地形，则提取的等高线的等高距要相应的小一些，因为等高距太大不能很好的反应填补洼谷地形周围地形的走势；如果填补的洼谷地形级别较高，且洼谷地形多为较长的宽阔洼谷地形，则提取的等高线的等高距要相应大一些，因为等高距太大会使计算量陡增。实施例采用arcgis中的“表面分析”工具提取等高线，arcgis采用矩形网格线性内插法提取规则网格数字高程模型中的等高线，该方法提取等高线首先确定规则网格数字高程模型中的等值点，之后把这些等值点连接成等值线，规则网格数字高程模型中等高线分闭合曲线和开曲线，其中开曲线与数据边界相交，所以从边界开始搜索等值点。

在计算提取等高线与填补洼谷地形边界的交点前需进行若干数据转换，规则网格数字高程模型数据格式是一个规则的行列式矩阵，本发明用一个二维数组存储规则网格数字高程模型数据，行列号表示了网格单元的位置。行列式矩阵中，原点在左上角，行号向右增加，列号向下增加。基于规则网格数字高程模型提取的等高线数据格式是shape格式，提取的矢量格式等高线数据坐标系是世界坐标，为减少等高线与填补洼谷地形边界交点的计算量，把等高线的世界坐标转换成行列式表示，即通过行列号表示点的位置，先获得等高线的世界坐标范围，由于世界坐标系x轴向右增加，y轴向上增加，把世界坐标点转换成行列式矩阵的公式如下：

其中(zx，zy)为等高线的左上角世界坐标，(dx，dy)是当前转换点的坐标，(jx，jy)是规则网格数字高程模型的坐标间隔,(r，c)是转换后的等高线点对应对的行列号。坐标转换完之后，依次遍历等高线所有点集，根据行列号查看网格数字高程模型单元，如果该单元是要填补洼谷地形的边界则标记该网格单元为洼谷地形边界与等高线的交点，同时标识等高线的识别码，方便后续的插值与填补操作。按照上述方法，只需遍历一次等高线就可以计算出所有与带填补洼谷地形边界的交点。

(二)确定洼谷地形边界与等高线交点处的等高线走势

人工地貌综合删减洼谷地形时，需要拉直洼谷地形周围的等高线，为了使洼谷地形内的等高线圆滑过渡，且与周边地形保持一致，需要事先确定待综合洼谷地形周围等高线的趋势，这样才能使综合后的地形与实际相符，判断待删减洼谷地形周围的等高线趋势是主观认知过程，所以自动识别判断待填补洼谷地形周围等高线的趋势是本发明需要解决的问题。

本发明采用等高线与带填补洼谷地形边界交点处的切线反应等高线的趋势，在网格数字高程模型试验区提取的等高线有两种曲线状态：闭合曲线和开曲线，开曲线是与网格数字高程模型边界相交的等高线，闭合曲线是封闭的曲线。求曲线的切线通过对曲线方程求导计算出相应曲线点的切线，提取的等高线的曲线方程是未知的，通过等高线与待填补洼谷地形边界交点处附近的若干点曲线拟合生成曲线方程，求出曲线方程，但计算量大不利于大数据量的处理。

(三)洼谷地形边界交点之间采用分段三次多项式插值

现有技术的洼谷地形填补是沿着洼谷地形法线方向线性填补，这种填补方法简单，但是填补后洼谷地形的等高线太过平直，与填补洼谷地形周围的地形走势不协调。为了改善填补后地形与周围地形协调性，本发明提出根据等高线与洼谷地形边界交点和在交点处的等高线的切线方向来进行三次多项式插值填补洼谷地形的方法。

三次多项式插值假设平面上有n个已知点(xi，yi)(i＝1,2…n)，要曲线拟合出一条光滑的曲线y＝f(x)，使得yi＝f(xi)(i＝1,2…n)。根据已知点三次多项式插值出的曲线需要满足二个要求：一是yi＝f(xi)(i＝1,2…n)，即曲线通过所有的已知点；二是曲线光滑，直观上线条不能有明显的折角，数学上是曲线方程任意一点都存在一阶导数yi’＝f’(xi)。三次多项式插值是用三次多项式无限逼近相邻两点，设有二个相邻的点di和di+1，在二个点之间进行三次多项式插值，其多项式曲线方程：

x＝d0+d1z+d2z²+d3z³

y＝e0+e1z+e2z²+e3z³

式中di、ei(i＝0、1、2、3)为常数，z为可变参数，当点di沿着三次多项式曲线变化到点di+1时，参数z是从0变化到1。

通过确定等高线与待填补洼谷地形边界的交点，计算在交点处等高线的近似切线与x轴的夹角，已经获得在二个同一条等高线在洼谷地形边界二个交点间进行三次多项式插值所需的条件。计算di和ei(i＝0、1、2、3)，使参数z从0变化到1进行点di和点di+1间的三次多项式插值，在插值之前需要进行坐标变换。

因为本发明网格单元的位置是采用行列号表示，行列式矩阵的列号向右增加，行号向下增加，这与三次曲线方程的坐标系不符，需要把行以原点为中心旋转180°，即把行号变为负值。坐标转换后需根据三次多项式的公式计算二个交点之间的插值点坐标，从点di插值到点di+1时，参数z从0变到1，但是不同洼谷地形边界交点之间的距离不同，如果参数z从0变化到1的间隔过大，有可能二个交点之间的插值不完全，插值点之间出现间断不连续的情况；如果参数z从0变到1的间隔过小，则会计算出很多重复的插值，造成计算资源浪费，影响计算速度。

所以本发明提出在计算插值之前先计算二个洼谷地形边界交点之间的距离有多少个网格单元，记为f，参数z从0变化到1的插值间隔dz＝1/2f，这可以满足大多数洼谷地形边界交点之间的三次多项式插值。插值结束后获得二个洼谷地形边界交点之间的插值数组，数组存储二个洼谷地形边界交点之间连接的洼谷地形单元，把洼谷地形内的插值单元高程值都赋值为等高线的高程值。

(四)内插等高线间的洼谷地形内部点

填补洼谷地形时考量周围等高线的趋势，确定等高线与待填补洼谷地形边界的交点和交点处等高线的近似切线夹角，根据洼谷地形边界中的交点和等高线在交点处的切线的夹角进行三次多项式插值，插值生成网格单元的高程值被赋以相应等高线的高程，这在填补洼谷地形内会形成若干新等高线，如图4所示。这些新生成的插值等高线形似一系列不同程度弯曲的阶梯，但是不同插值等高线间的待填补洼谷地形单元的高程值没有变换，为了使洼谷地形内网格单元在新生成的插值等高线间平滑过渡，需要对没有变换的洼谷地形内网格单元进行线性内插。

(1)确定填补洼谷地形粗略走势

填补洼谷地形内没有变换过的网格线性内插前，需确定洼谷地形的粗略走势，按照洼谷地形的粗略走势内插可使填补后的地形与实际地形一致，保持洼谷地形从高处向洼谷地形的低处过渡的趋势。洼谷地形线的趋势能很好表示洼谷地形的粗略走势，但很多洼谷地形的上游点没有到达洼谷地形的源头，因为采用物理流水模拟法提取集水线时根据集水流量阈值确定集水点，洼谷地形源头的集水流量等于0，而集水流量阈值不能过小，导致提取的集水线上游点不能到达洼谷地形源头。

所以在确定洼谷地形粗略走势前，需要查看洼谷地形线的上游点是否到达洼谷地形源头，如果洼谷地形线上游点没有到达洼谷地形源头，即上游点的集水量不等于0，需要进行洼谷地形线的延长。洼谷地形线的延长方法为：

步骤一，所有流向洼谷地形上游点的相邻网格单元为候选上游点；

步骤二，从候选的网格单元中选出集水流量最大的点为新的洼谷地形上游点，如果集水流量最大的网格单元不止一个，则选择最邻近的网格单元为新的洼谷地形上游点；

重复上述步骤一和步骤二直到上游点的集水流量等于0，则洼谷地形线的延长结束。

现有技术根据洼谷地形线各单元流向确定洼谷地形粗略走势的方法为：统计洼谷地形线各单元流向的个数，其中数量最多的流向反应洼谷地形的粗略走势。当洼谷地形没有过度弯曲时该方法能反应洼谷地形的走势，但是当洼谷地形弯曲程度很大时，这样确定洼谷地形走势的方法在局部无法反应洼谷地形的走势。所以本发明采用洼谷地形线各单元流向分段确定洼谷地形粗略走势的方法，如图5所示，分别统计新插值生成等高线间的洼谷地形线单元的流向个数，数量最多的流向为插值生成的等高线间洼谷地形的粗略走势。

(2)等高线间的线性填补

通过三次多项式插值法生成待填补洼谷地形新的等高线和确定等高线间洼谷地形的粗略走势后，进行待填补洼谷地形内没有变换的洼谷地形单元的线性填补。具体方法为：依次遍历所有的网格单元，如果该单元为待填补洼谷地形内，且高程值没有变换，则向插值新生成的等高线间洼谷地形的粗略走势的方向两端扩充，如图6所示，直至遇到当前洼谷地形的边界或插值生成等高线的单元结束，获得一个首尾洼谷地形边界单元或新生成等高线单元的数组，计算数组两端的高程差，并对数组内其它单元进行线性内插，使得首尾两端的单元平滑过渡。

五、洼谷地形填补的数字高程模型金字塔

为制作直观和表现力强的数字地貌晕渲图，需对地貌进行合理综合，在保持或夸张主体地貌特征的基础上，综合不必要的地形细节信息，突出地貌的主体特征。

本发明采用基于洼谷地形填补的方法进行地貌综合，综合过滤地貌晕渲中过度零碎的地形细节信息。采用不同综合程度生成不同详细程度的数字高程模型图层，不同详细程度的图层构成一个地形详细程度递减的数字高程模型金字塔，通过设置各图层的权重系数，调整相应精细程度的图层在重建地形中的比重，突出关键地形特征。

为使生成的不同详细程度的数字高程模型图层连续且方法简单可行，本发明提出填补某一级别所有的洼谷地形生成相应详细程度的数字高程模型图层，假设提取的洼谷地形有m个级别，填补第i级洼谷地形生成数字高程模型图层ai(i＝1、2…m-1)，当填补的洼谷地形级别i大于1时，洼谷地形的识别与洼谷地形边界的确定与第一级洼谷地形不同，因为高级别洼谷地形内经常包含有与其相连的次级洼谷地形，所以在洼谷地形识别和确定洼谷地形边界时要做特别处理。高级别洼谷地形内的与其相连的低级别洼谷地形属于当前高级别洼谷地形，所以在洼谷地形识别时根据洼谷地形编码确定高级洼谷地形和低级洼谷地形的关系，因为本发明提出的洼谷地形编码原则是与其相连的父级洼谷地形编码加上本身洼谷地形级别和洼谷地形编号，所以只要低级别洼谷地形编码的头几位与高级别洼谷地形编码相同则确定低级别洼谷地形属于当前高级别洼谷地形。依次填补不同级别洼谷地形生成不同详细程度的数字高程模型图层ai(i＝0,1,2,3…m)，其中a0表示没有任何综合的原始地形图层，各个数字高程模型图层设置一个权重系数vi(i＝0，1,2,3…n)，重建的地形g公式如下所示：

各图层的权重系数根据综合尺度和不同地形特征调整，若综合的尺度较小，则低级别图层的权重系数相应大一些；若要突出地形的主要特征，则设置高级别图层的权重系数大一些。

本发明提供的一种表现力强的地貌晕渲图洼谷地形填补方法，在删减地形细节信息时能较好的保持地形特征，该方法构建数字高程模型金字塔，通过填补不同级别的洼谷地形生成不同综合程度的数字高程模型图层，调整各个数字高程模型图层的权重系数来突出或则削弱某一图层，在删减地貌碎部的同时较好的保持地形特征，而且还能通过权重系数的设置突出地形的关键地形特征，制作出一览性好重点突出且三维立体感强的数字地貌晕渲图。本发明填补洼谷地形时考量周围等高线在洼谷地形内走势，并分段插值洼谷地形内网格单元高程值，使得填补后的洼谷地形内地形与周围地形平滑过渡，数字地貌晕渲图的洼谷地形表达更加准确自然，数字地貌晕渲图更加清楚美观。