CN110457772A - 一种结合平面曲率和最陡下坡方向的dem流向估计方法 - Google Patents

Abstract

本发明提供了一种结合平面曲率和最陡下坡方向的DEM流向估计方法，包括以下步骤：步骤S1：加载数字高程模型；步骤S2：计算最陡下坡方向所在三角形平面的3个顶点单元的平面曲率；步骤S3：使用步骤S2计算得到的3个单元的平面曲率的比值对最陡下坡方向进行调整。本发明的有益效果为：本发明结合了平面曲率和最陡下坡方向，发明了一种估算数字高程模型水流方向的方法，相较于传统的直接使用最陡下坡方向的方法，该方法具有精度更高，对收敛的主要河道还原效果更好，考虑地形起伏变化的优势，同时该方法相对传统方法增加的步骤不多，不会引起太多的计算机运行时间的差异。

Description

一种结合平面曲率和最陡下坡方向的DEM流向估计方法

技术领域

本发明涉及数字地形分析技术领域，尤其涉及一种结合平面曲率和最陡下坡方向的DEM流向估计方法。

背景技术

水流方向，指的是水、土壤、溶质在重力作用下在地球表面的移动方向，是许多水文学、地貌学领域的数值模拟中必须确定的参数。由于真实的复杂地形无法完全存储进计算机，因此目前的水文地貌模型主要以数字高程模型作为地形数据源，水流方向的计算也普遍基于数字高程模型进行。

O’Callaghan和Mark于1984年提出了最早的水流方向计算方法，通过比较中心单元与8个相邻单元间的坡度，选取与中心单元间坡度最大的相邻单元作为下游单元，将中心单元的流向指向该单元中点。由于这种方法只允许从8个方向中选取流向，因此被称为八方向(D8)方法。由于真实的水流方向可能不是8个允许方向之一，因此D8方法存在较大误差。为了降低D8方法的误差，Tarboton(1997)提出以中心单元及8个相邻单元构造8个三角形平面，确定最陡下坡方向作为水流方向的方法，使用该方法计算得到的水流方向可以为0°-360°之间的任何值，精度相比D8方法大大提高。该方法被命名为无穷流向(Dinf)方法。然而由于Dinf方法在构造三角形平面的过程中忽略了地形存在的弯曲起伏，因此在地形并不平坦是同样存在误差，尤其是当地形弯曲程度较大时误差尤为严重。Hooshyar等人于2016年提出了使用切面曲率调整最陡下坡方向的方法，但是他们的方法只适用于优化收敛型地形的流向。由于完全通过数学方程还原地形表面并使用方向导数获取流向的理论方法在现实中难以实现，因此只能使用一定方法将最陡下坡方向调整得更接近真实流向。

目前常用的地形曲率主要有平面曲率、剖面曲率、切面曲率3种，其中平面曲率能够很好地反映地形发散、收敛的程度，包含着大量Dinf方法确定最陡下坡方向过程中无法获得的信息。借助平面曲率和最陡下坡方向，确定的水流方向就同时考虑到了局部坡度和地形起伏信息。

发明内容

本发明的目的在于，为了解决现有无穷流向方法在弯曲地形上确定的流向存在一定误差的问题，而提供一种结合平面曲率和最陡下坡方向的DEM流向估计方法，采用了平面曲率对最陡下坡方向进行调整，进而确定水流方向，是估计数字高程模型中栅格单元水流方向的方法。

本发明是通过如下措施实现的：一种结合平面曲率和最陡下坡方向的DEM流向估计方法，其中，包括以下步骤：

步骤S1：加载数字高程模型，将需要确定水流方向的栅格单元作为中心单元，使用中心单元及与其相邻的8个单元的中心点划分出8个三角形平面，比较每个平面的最陡坡度，选取最陡坡度值最大的平面，记录该平面，并将该平面的最陡坡度对应的坡度方向作为中心单元的最陡下坡方向；

步骤S2：计算最陡下坡方向所在三角形平面的3个顶点单元的平面曲率；

步骤S3：以最陡下坡方向所在三角形平面限制的4/π角度为允许变化范围，使用步骤S2计算得到的两个下游单元的平面曲率与中心单元平面曲率差值的比值对最陡下坡方向进行调整。

进一步地，作为本发明提供的一种结合平面曲率和最陡下坡方向的数字高程模型水流方向估计方法的进一步优化方案，所述步骤S2计算平面曲率的方程为：

其中z₁,…,z₉的分布如图4所示，其中z1,…,z9表示当前计算单元及其周边8个相邻单元的高程，K为中心单元z₅的平面曲率。

进一步地，作为本发明提供的一种结合平面曲率和最陡下坡方向的数字高程模型水流方向估计方法的进一步优化方案，所述步骤S3对步骤S2中计算得到的两个下游单元的平面曲率与中心单元平面曲率差值的比值对最陡下坡方向进行调整依照下面的方程进行：

当平面曲率K₁和K₂的值为一正一负时，修改后的水流方向α^*指向正平面曲率对应的单元中心，即水流放弃发散单元指向收敛单元，即：若K₁>0并且K₂≤0：

若K₂>0并且K₁≤0：

式中i是步骤S2选取三角形平面的编号，当平面曲率K₁和K₂的值同为正或同为负时，使用下游单元与中心单元的平面曲率差值K₁-K₀和K₂-K₀的比值控制转动幅度，确定旋转后水流方向α^*的方程为：

该方向以正上方为0rad，并顺时针增大；针对α^*的值可能超出允许范围的情况，需对其进行修正，当时，当时，经过修正后最终流向即确定。

进一步地，作为本发明提供的一种结合平面曲率和最陡下坡方向的数字高程模型水流方向估计方法的进一步优化方案，所述步骤S1的三角形平面的最陡坡度s及对应的坡度方向α可以由以下方程得到：

式中：

其中(x₁，y₁，z₁)，(x₂，y₂，z₂)，(x₃，y₃，z₃)为所计算的三角形平面的三个顶点的三维坐标，由于从中心点出发，最陡坡度方向可能落于三角面范围外，对坡度方向进行修正，将坡度方向调整到与已计算方向最近的三角面边缘方向，并同步将坡度修改为该方向坡度，修改后的方向α_m为：

其中，式中i是步骤S2选取三角形平面的编号；针对α^*的值可能超出允许范围的情况，需对其进行修正，当时，当时，经过修正后确定最终流向。

所述DEM流向即为数字高程模型水流方向。

本发明的有益效果为：本发明结合了平面曲率和最陡下坡方向，利用不同下坡位置平面曲率的不同对最陡下坡方向进行调整，提供了一种获取数字高程模型中水流方向的方法；相较于传统的完全依赖最陡下坡方向的方法，该方法具有精度更高，对收敛的主要河道还原效果更好，考虑地形起伏变化的优势；本发明为进行水、土壤、溶解质输移模拟的水文地貌模型提供了获取更高精度水流方向的方法，同样依靠数字高程模型就能得到比以往更高精度的水流方向估计结果；同时该方法相对传统方法增加的步骤不多，不会引起太多的计算机运行时间的差异。

附图说明

图1为本发明的整体流程图；

图2为本发明实施例中最陡下坡方向计算时的3×3窗口划分示意图，图中的中心点P₀二维坐标为(0，0)；

图3为本发明实施例中切面曲率计算时的3×3窗口高程示意图；

图4为本发明实施例中最陡下坡方向确定后需计算的下坡单元平面曲率K1和K2对应的位置图；

图5为本发明实施例中的数字高程模型；

图6为本发明实施例中需要计算的中心单元及周边相邻单元组成的八个三角面各自的最大坡度示意图；

图7为本发明实施例中使用的屯溪流域DEM地形图；

图8(a)为本发明实施例中使用Tarboton的方法得到的主要河道图；

图8(b)为本发明实施例中使用本发明方法得到的主要河道图。

具体实施方式

为了能清楚说明本方案的技术特点，下面通过具体实施方式，对本方案进行阐述。

实施例1

参见图1至图8，本发明是：结合平面曲率和最陡下坡方向的DEM流向估计方法，其中，包括以下步骤：

步骤S1：加载数字高程模型，将需要确定水流方向的栅格单元作为中心单元，使用中心单元及与其相邻的8个单元的中心点划分出8个三角形平面，比较每个平面的最陡坡度，选取最陡坡度值最大的平面，记录该平面，并将该平面的最陡坡度对应的坡度方向作为中心单元的最陡下坡方向；

步骤S2：计算最陡下坡方向所在三角形平面的3个顶点单元的平面曲率；

步骤S3：以最陡下坡方向所在三角形平面限制的4/π角度为允许变化范围，使用步骤S2计算得到的两个下游单元的平面曲率与中心单元平面曲率差值的比值对最陡下坡方向进行调整。

其中，所述步骤S2计算平面曲率的方程为：

其中z₁,…,z₉的分布如图4所示，其中z1,…,z9表示当前计算单元及其周边8个相邻单元的高程，K为中心单元z₅的平面曲率。

其中，所述步骤S3对步骤S2中计算得到的两个下游单元的平面曲率与中心单元平面曲率差值的比值对最陡下坡方向进行调整依照下面的方程进行：

当平面曲率K₁和K₂的值为一正一负时，修改后的水流方向α^*指向正平面曲率对应的单元中心，即水流放弃发散单元指向收敛单元，即：若K₁>0并且K₂≤0：

若K₂>0并且K₁≤0：

式中i是步骤S2选取三角形平面的编号，当平面曲率K₁和K₂的值同为正或同为负时，使用下游单元与中心单元的平面曲率差值K₁-K₀和K₂-K₀的比值控制转动幅度，确定旋转后水流方向α^*的方程为：

该方向以正上方为0rad，并顺时针增大；针对α^*的值可能超出允许范围的情况，需对其进行修正，当时，当时，经过修正后最终流向即确定。

其中，所述步骤S1的三角形平面的最陡坡度s及对应的坡度方向α可以由以下方程得到：

式中：

其中(x₁，y₁，z₁)，(x₂，y₂，z₂)，(x₃，y₃，z₃)为所计算的三角形平面的三个顶点的三维坐标，由于从中心点出发，最陡坡度方向可能落于三角面范围外，对坡度方向进行修正，将坡度方向调整到与已计算方向最近的三角面边缘方向，并同步将坡度修改为该方向坡度，修改后的方向α_m为：

其中，式中i是步骤S2选取三角形平面的编号；针对α^*的值可能超出允许范围的情况，需对其进行修正，当时，当时，经过修正后确定最终流向。

实施例2

参见图1至图8，本发明还提供结合平面曲率和最陡下坡方向的DEM流向估计方法，具体包括以下步骤：

步骤S1：加载数字高程模型，将需要确定水流方向的栅格单元作为中心单元，使用中心单元及与其相邻的8个单元的中心点划分出如图2所示的8个三角形平面，比较每个平面的最陡坡度，选取最陡坡度值最大的平面，记录该平面，并将该平面的最陡坡度对应的坡度方向作为中心单元的最陡下坡方向；

步骤S1的具体内容为：将当前计算的栅格单元作为中心单元，如图2所示，连接中心单元的中心点P₀与其周围8个相邻单元的中心点(P₁、P₂、…、P₈)，划分出如图2中粗实线框出的8个三角形平面，并如图为这8个平面从1至8编号；

以P₀点的二维平面XY横纵坐标为(0，0)，另外8个点的二维平面坐标设置如图2所示；此时将每个栅格单元具备的高程值赋予其中心点，即知晓P₀至P₈这9个点的高程Z，则全部9个点的XYZ三维坐标已知，也即知晓每个三角形平面的三个顶点的三维坐标；此时可用下面的方程计算每个平面的最陡坡度s及对应的坡向α：

式中：

其中(x₁，y₁，z₁)，(x₂，y₂，z₂)，(x₃，y₃，z₃)为所计算的三角形平面的三个顶点的三维坐标，此处使用的3个三维坐标的顺序不会影响结果，因此无需考虑三个顶点坐标的顺序；方程(2)得到的坡度方向为弧度制，以正上方向为0rad，顺时针增加；

由于从中心点P₀出发，每个三角面仅有π/4的允许范围，而最陡坡度方向可能落于允许范围外，因此需对坡度方向进行修正，将坡度方向调整到与已计算方向最近的三角面边缘方向，并同步将坡度修改为该方向坡度，修改后的方向α_m为：

式中：n为对应三角形平面的编号；此时对应的坡度s_m方程为

式中：h₀为中心单元高程；h_a为α_m指向的目标单元的高程；r为中心单元中点与α_m目标单元中点的距离，对于位于中心单元正方向的目标单元，r值等于DEM的栅格边长，对于位于斜方向的目标单元，r值等于的栅格边长。

分别计算8个三角形平面的s_m值并进行比较，选取其中s_m值最大的三角形平面，以该平面的α_m值作为中心单元的最陡下坡方向α_m ^*，并记录该三角形平面编号；如果存在多个三角形平面具备最大的s_m值，则选择编号最小的三角形平面；

步骤S2：计算最陡下坡方向所在三角形平面的3个顶点单元的平面曲率，具体计算方法为：

同样使用以所计算单元及其周围单元组成的3×3窗口计算平面曲率，对于如图3所示的栅格高程分布，中心点Z₅的平面曲率K为：

式中：其中z1,…,z9表示当前计算单元及其周边8个相邻单元的高程，Z₁-Z₉是图3中3×3窗口内中心单元及其周边8个相邻单元的高程，Δx是栅格的边长；

步骤S3、以最陡下坡方向所在三角形平面限制的4/π角度为允许变化范围，使用步骤S2计算得到的3个单元的平面曲率对最陡坡度方向进行调整，得到最终的水流方向，具体为：

中心单元的平面曲率为K₀，除中心单元外，在步骤S2中所选取的三角形平面的2个顶点，必有1个处于中心单元的上下左右4个正方向之一，另一个处于左上、左下、右上、右下4个斜方向之一；以正方向的顶点单元的平面曲率记为K₁，斜方向单元平面曲率为K₂；如图4，若最陡下坡方向α_m ^*处于P₁、P₂、P₀所组成的三角形平面，K₁则为P₁的平面曲率，K₂为P₂的平面曲率，K₀为P₀的平面曲率；由于平面曲率表现了地形的发散与收敛，平面曲率数值越大地形越收敛，因此需要将最陡下坡方向向高平面曲率的下游方向转动。

当平面曲率K₁和K₂的值为一正一负时，修改后的水流方向α^*指向正平面曲率对应的单元中心，即水流放弃发散单元指向收敛单元，即是：若K₁>0并且K₂≤0：

若K₂>0并且K₁≤0：

式中i是步骤S2选取三角形平面的编号。当平面曲率K₁和K₂的值同为正或同为负时，使用下游单元与中心单元的平面曲率差值K₁-K₀和K₂-K₀的比值控制转动幅度，具体确定旋转后水流方向α^*的方程为：

式中i是步骤S2选取三角形平面的编号；针对α^*的值可能超出允许范围的情况，需对其进行修正，当时，当时，经过修正后确定最终流向。

为了更清楚说明本发明的技术特点，下面通过具体实例，对本方案进行详细阐述本发明确定水流方向的流程和优势如下：

以图5所示的数字高程模型为例，根据步骤1，将灰色单元和其周围8个单元连接划分三角形平面，计算得到各平面的最陡坡度如图6所示。可见8个平面的最大坡度中，编号为4的三角形平面最大坡度的数值最大，因此最陡下坡方向为该平面的方向，即α_m ^*＝2.840rad。

选取编号为4的三角形平面除中心点外的两个顶点计算平面曲率，位于中心单元正向，高程43.860m的单元的平面曲率对应K₁，位于中心单元斜向，高程43.103m的单元的平面曲率对应K₂，根据方程(5)计算得到，K₀＝-0.319，K₁＝-0.271，K₂＝-0.257。使用步骤S3中的方程(6)得到最终的水流方向α^*为2.752rad。

另一个具体实例为如下：

选取屯溪流域的DEM，如图7所示，分别使用Tarboton(1997)的最陡下坡方向方法和本发明的方法确定DEM单元流向，然后使用Tarboton(1997)针对此类流向提出的汇流计算方法确定流域汇流值，并提取出图8所示的流域内集水面积大于4.5km²的河道。使用两种方法确定水流方向后得到的河道在图中方框区域存在明显差异，使用Tarboton(1997)的方法得到的河道过于比值，且右侧河道与实地勘查中的向西大幅弯曲后折回的走向不合，而本发明得到的流向还原的河道更接近真实河道。

综上所述，本发明结合了平面曲率和最陡下坡方向，发明了一种估算数字高程模型水流方向的方法，相较于传统的直接使用最陡下坡方向的方法，该方法具有精度更高，对收敛的主要河道还原效果更好，考虑地形起伏变化的优势，同时该方法相对传统方法增加的步骤不多，不会引起太多的计算机运行时间的差异。

本发明未经描述的技术特征可以通过或采用现有技术实现，在此不再赘述，当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的普通技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (4)

1.一种结合平面曲率和最陡下坡方向的DEM流向估计方法，其特征在于，包括以下步骤：

步骤S1：加载数字高程模型，将需要确定水流方向的栅格单元作为中心单元，使用中心单元及与其相邻的8个单元的中心点划分出8个三角形平面，比较每个平面的最陡坡度，选取最陡坡度值最大的平面，记录该平面，并将该平面的最陡坡度对应的坡度方向作为中心单元的最陡下坡方向；

步骤S2：计算最陡下坡方向所在三角形平面的3个顶点单元的平面曲率；

步骤S3：以最陡下坡方向所在三角形平面限制的π/4角度为允许变化范围，使用步骤S2计算得到的两个下游单元的平面曲率的比值对最陡下坡方向进行调整。

2.根据权利要求1所述的结合平面曲率和最陡下坡方向的DEM流向估计方法，其特征在于，所述步骤S2计算平面曲率的方程为：

其中z₁,…,z₉表示当前计算单元及其周边8个相邻单元的高程，K为中心单元z₅的平面曲率。

3.根据权利要求1或2所述的结合平面曲率和最陡下坡方向的DEM流向估计方法，其特征在于，所述步骤S3对步骤S2中计算得到的两个下游单元的平面曲率与中心单元平面曲率差值的比值对最陡下坡方向进行调整依照下面的方程进行：

当平面曲率K₁和K₂的值为一正一负时，修改后的水流方向α^*指向正平面曲率对应的单元中心，即水流放弃发散单元指向收敛单元，即：若K₁>0并且K₂≤0：

若K₂>0并且K₁≤0：

式中i是步骤S2选取三角形平面的编号，当平面曲率K₁和K₂的值同为正或同为负时，使用下游单元与中心单元的平面曲率差值K₁-K₀和K₂-K₀的比值控制转动幅度，确定旋转后水流方向α^*的方程为：

该方向以正上方为0rad，并顺时针增大；针对α^*的值可能超出允许范围的情况，需对其进行修正，当时，当时，经过修正后最终流向即确定。

4.根据权利要求1-3任一项所述的结合平面曲率和最陡下坡方向的DEM流向估计方法，其特征在于，所述步骤S1的三角形平面的最陡坡度s及对应的坡度方向α可以由以下方程得到：

式中：

其中(x₁，y₁，z₁)，(x₂，y₂，z₂)，(x₃，y₃，z₃)为所计算的三角形平面的三个顶点的三维坐标，由于从中心点出发，最陡坡度方向可能落于三角面范围外，对坡度方向进行修正，将坡度方向调整到与已计算方向最近的三角面边缘方向，并同步将坡度修改为该方向坡度，修改后的方向α_m为：

其中，式中i是步骤S2选取三角形平面的编号；针对α^*的值可能超出允许范围的情况，需对其进行修正，当时，当时，经过修正后确定最终流向。