一种结合切面曲率和Dinf方法的DEM水流方向模拟方法与流程

本发明涉及数字地形分析技术领域，尤其涉及一种结合切面曲率和dinf方法的dem水流方向模拟方法。

背景技术：

基于dem结构的水流方向被广泛应用于水文预报、土壤厚度演化、污染物运移等研究领域，然而目前还没有公认的最佳水流方向模拟方法，已有的方法或多或少存在一定问题，这严重限制了以水流方向为基础的相关研究的成果精度。

早期的经典d8方法(o’callaghan和mark，1984)由于只允许水流向8个可行方向之一，与现实世界水流在0°至360°中任意方向流动严重不符，进而导致了较大的误差存在。因此，tarboton于1997年提出了允许水流方向指向任意方向的dinf方法，初步解决了d8方法所受的限制。然而由于dinf方法是通过把所计算单元周边区域划分出8个三角形平面来实现水流方向计算的，这种将起伏的复杂地形简化为平面的做法会在弯曲起伏的地形造成一定误差。这种情况最好的解决方法就是将地形起伏纳入水流方向计算的考虑范围。

目前反映地形变化主要使用地形曲率来完成，包括平面曲率、剖面曲率、切面曲率，这其中又以切面曲率与水流过程关系最为密切，能更为准确地反映地形对水流的收敛、发散影响。hooshyar等(2016)就曾提出利用切面曲率来优化dinf流向的方法，然而他们的方法只能在切面曲率数值为正的收敛地形下适用。

技术实现要素：

本发明的目的在于为了解决现有无穷流向方法由于地形起伏导致得到的水流方向精度不够，不能做到精确汇流，最终获取的河网不够准确的问题，采用切面曲率对dinf方向进行修正，进而模拟水流方向，而提出种结合切面曲率和dinf方法的dem水流方向模拟方法。

本发明是通过如下措施实现的：一种结合切面曲率和dinf方法的dem水流方向模拟方法，包括以下步骤：

步骤s1：加载dem，检索需要计算水流方向的栅格单元及其周围的8个相邻单元，连接它们的中心点划分出dinf方法需要的8个三角形平面，并计算dinf流向，设置dinf流向以正北方向为0°顺时针增大；

步骤s2：如果dinf流向在两个三角形平面的交线方向，则直接选择dinf方向作为最终流向，不再进行后续步骤，否则确定dinf流向所在的三角形平面，确定该平面上除了中心点p0外的另外2个顶点对应栅格单元的切面曲率；

步骤s3：使用切面曲率对dinf流向进行校正，将dinf流向向高切面曲率，即是更为收敛的下坡方向旋转一定角度。

作为本发明提供的一种结合切面曲率和dinf方法的dem水流方向模拟方法进一步优化方案，所述步骤s3中更为收敛的下坡方向旋转一定角的具体步骤如下：

(1)、如果2个切面曲率中只有1个属于收敛型，即是只有1个的切面曲率值大于0，最终确定的方向指向大于0的切面曲率值对应的栅格单元中心；

(2)、如果2个切面曲率同为正值或同为负值，则最终水流方向比dinf流向更接近其中收敛度大或发散度低的地形，将最终流向确定为切面曲率值更大的单元中心所在方向和dinf流向的角平分线，该角平分线方向α计算方程为：

其中αk为高切面曲率单元中心所在方向，αd为dinf流向；设置方向指向正北时都为0°，当αk＝0且αd＞315°时，使用αk＝360°代入上述角平分线方向α计算方程中。

作为本发明提供的一种结合切面曲率和dinf方法的dem水流方向模拟方法进一步优化方案，所述步骤2中切面曲率的具体计算方法如下：

假如p0的dinf流向如图4所示的箭头方向，即其所在的三角形平面由p0、p1、p2构成，需要计算切面曲率的点就是p1和p2；以需要计算的点的高程为z5，切面曲率计算需要的3×3窗口各栅格高程如图5所示，切面曲率的方程为：

其中：

其中△x是dem栅格的边长。

本发明的有益效果为：

(1)本发明使用下游方向切面曲率的区别对dinf流向进行调整，发明了一种新型的模拟dem水流方向的算法，该方法实现步骤简单，并未在dinf方法的基础上增加过多步骤。相比于dinf方法，本发明具有在部分地区还原的水系网络更为真实的优势。

(2)本发明通过将dinf流向向更收敛的下坡方向转动，使得得到的水流方向更贴近真实情况下水流向山谷汇聚的特征；经验证，本发明提供的水流方向比dinf方法在部分地形更加精确，本发明为水文模型获取水流方向以及更进一步的水系河网提供了新的方法。

(3)本发明为各类水文模型获取高精度水流方向，从而借此提取高精度水系河网提供了更好的解决方案。

附图说明

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

图2为本发明的dinf方法计算水流方向时划分周围区域的示意图；

图3为本发明的方向与角度对应图；

图4为本发明中dinf方法模拟效果；

图5为本发明的切面曲率计算时的3×3窗口；

图6(a)为本发明的实施例一使用的dem；

图6(b)为本发明的实施例一得到的最终水流方向；

图7为本发明测试的的鱼梁流域dem地形图；

图8(a)为本发明的用于对比的使用dinf方法还原的鱼梁流域水系；

图8(b)为本发明的使用本发明方法还原的鱼梁流域水系。

具体实施方式

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

参见图1至图8，本发明是：一种结合切面曲率和dinf方法的dem水流方向模拟方法，包括以下步骤：

步骤s1、加载dem，检索需要计算水流方向的栅格单元及其周围的8个相邻单元，连接它们的中心点划分出dinf方法需要的8个三角形平面(如图2，其中的p0为需要计算水流方向的单元)，并计算dinf流向，为便于后续步骤的方程表达，此处如图3设置dinf流向以正北方向为0°顺时针增大；

步骤s2、如果dinf流向在两个三角形平面的交线方向，也即dinf流向指向某个相邻单元的中心点，则直接选择dinf方向作为最终流向，不再进行后续步骤，否则确定dinf流向所在的三角形平面，确定该平面上除了中心点p0外的另外2个顶点对应栅格单元的切面曲率，具体为：

如图4，假如p0的dinf流向如图4中箭头方向，即其所在的三角形平面由p0、p1、p2构成，需要计算切面曲率的点就是p1和p2；以需要计算的点的高程为z5，切面曲率计算需要的3×3窗口各栅格高程如图5所示，切面曲率的方程为：

其中：

其中△x是dem栅格的边长；

步骤s3、使用切面曲率对dinf流向进行校正，将dinf流向向高切面曲率，也就是更为收敛的下坡方向旋转一定角度，具体步骤如下：

(1)如果2个切面曲率中只有1个属于收敛型，也就是只有1个的切面曲率值大于0，此时最终确定的方向指向大于0的切面曲率值对应的栅格单元中心，如图4中的情形，若k1>0，k2≤0，则最终水流方向从p0指向对应栅格单元中心p1点；

(2)如果2个切面曲率同为正值或同为负值，即同为收敛地形或发散地形，则将认为最终水流方向比dinf流向更接近其中收敛度大或发散度低的地形，将最终流向确定为切面曲率值更大的单元中心所在方向和dinf流向的角平分线；如图4所示，若k1>k2>0，则最终流向为点p1相对p0所在的0°方向和黑色箭头表示的dinf流向的角平分线方向；该角平分线方向α计算方程为：

其中αk为高切面曲率单元中心所在方向，αd为dinf流向；由于0°和360°的重合可能导致方程(7)的计算结果存在问题，因此设置方向指向正北时都统一为0°，当αk＝0且αd＞315°时使用αk＝360°代入上面方程；

此步骤确定的水流方向就是本发明得到的最终水流方向。

为了更好地验证本发明模拟水流方向的过程和优越性，以下两个实际例子进行详细说明。

(一)

以图6(a)展示的dem作为例子，计算其中灰色单元的水流方向。

根据步骤s1使用dinf方法得到的水流方向αd＝162.7°，具体方向如图6(b)中黑色箭头所示；

根据步骤s2中dinf方向所在三角面需要计算切面曲率的两个单元分别为高程为43.860m和43.103m的单元，二者切面曲率为-0.00927和-0.00906；

根据步骤s3，由于位于灰色单元东南方向的单元切面曲率值更大，故αk＝135°(图6(b)中虚线箭头方向)，所以最终确定的水流方向(图6(b)中虚白色箭头方向)。

(二)

使用图7所示的鱼梁流域dem验证本发明的应用效果，使用tarbonton(1997)提出的方法对dinf流向和本发明的水流方向进行汇流计算，然后提取出汇流面积大于4.5km²的dem单元作为主河道。dinf方法得到的河道如图8(a)，本发明得到的河道如图8(b)，可见二者间的主要区别在框出的位置，dinf方法提供的河道在该处过于直，而本发明很好地还原了河道的弯曲，并且与实地勘测结果相近。

综上所述，本发明使用下游方向切面曲率的区别对dinf流向进行调整，发明了一种新型的模拟dem水流方向的算法，该方法实现步骤简单，并未在dinf方法的基础上增加过多步骤，相比于dinf方法，本方法还原的水系网络更为真实。

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