CN112989639B - 一种基于平均化处理的dem栅格局部排水方向确定方法 - Google Patents

Abstract

本发明提供了一种基于平均化处理的DEM栅格局部排水方向确定方法，其技术方案为：通过将切面曲率作为一项判别依据，在以需要确定排水方向的DEM栅格为中心的3×3窗口内，使用2种划分方法构建2个三角形平面，使用2个平面坡向的平均值作为最终的局部排水方向。本发明的有益效果为：本发明为基于DEM的地表径流过程精确模拟提供了准确模拟排水方向的方法。

Description

一种基于平均化处理的DEM栅格局部排水方向确定方法

技术领域

本发明涉及数字地形分析技术领域，尤其涉及一种基于平均化处理的DEM栅格局部排水方向确定方法。

背景技术

排水方向是指水受到重力影响而在地表流动的方向。目前的水文模型大多基于水的运动轨迹进行模拟计算，因此准确地获取地表各处的排水方向对建模研究非常重要。由于水文模型大多使用数字高程模型(Digital Elevation Model，DEM)代表地形表面的高程分布，排水方向算法也大多以栅格型DEM作为数据基础。

由于在栅格DEM中，每个栅格单元周边有8个单元相邻，因此目前普遍使用中心单元与周边8单元间的高程差异确定中心单元的排水方向。O’Callaghan和Mark(1984)提出的D8方法是最经典的排水方向计算方法，该方法通过比较中心单元与相邻单元间的坡度值，选取坡度最大的相邻单元的中心点作为排水方向。因为D8方法只考虑了8个方向，与现实中水流可以流向0°-360°间任意方向的情况并不相合。Tarboton(1997)提出了一种通过划分8个三角形平面进行排水方向判断的Dinf方法，不过该方法相对固定的三角形平面划分方式并不能得到最优结果。

如何解决上述技术问题为本发明面临的课题。

发明内容

本发明的目的在于，针对现有排水方向确定方法精度不足的问题，而提供一种基于平均化处理的DEM栅格局部排水方向确定方法，通过使用不同的选点方式构建2个三角形平面，并使用2个平面坡向的平均值作为最终的排水方向，确定DEM中任意栅格对应排水方向的方法。

本发明是通过如下措施实现的：一种基于平均化处理的DEM栅格局部排水方向确定方法，其中，所述方法具体步骤如下：

(1)、对于需要确定排水方向的栅格C₀，令其周边相邻的8个栅格按顺序编号为C₁-C₈，并按顺序连接上述9个栅格的中心点，构建出8个三角形平面；

(2)、比较每个三角形平面的另外两个顶点，剔除另外两个顶点高程均高于中心单元C₀的平面，计算剩下的三角形平面的坡度和坡向，若坡向不在三角形平面对应的范围内，则将坡向调整为最接近原方向的允许范围边缘方向，坡度为沿新坡向的坡度，然后选取坡度最大的三角形平面所对应的45°范围作为最终排水方向的可能范围，将此单元的坡向定义为α₁；

(3)、除C₀外，将步骤(2)中选取的坡度最大三角形平面上的另外2个顶点定义为C₀的2个下游单元，计算2个下游单元的切面曲率，并比较其切面曲率大小；

(4)、除步骤(3)选取的2个下游单元外，从C₀的另外6个相邻单元中选取1个新的单元，该单元最接近具有较高切面曲率的下游单元，使用新选取的单元与C₀、具有较低切面曲率的下游单元构建三角形平面，计算新平面的坡向α₂；

(5)、将单元C₀的排水方向α确定为α₁和α₂所夹锐角的角平分线方向，若该角平分线方向不在步骤(2)确定的允许范围内，则将排水方向α定为最接近该角平分线方向的允许范围边缘方向。

作为本发明提供的一种基于平均化处理的DEM栅格局部排水方向确定方法的进一步优化方案，所述步骤(2)具体按以下步骤进行：

(2-1)、设组成三角形平面的三个DEM栅格单元中心点坐标分别为(x₁，y₁，z₁)，(x₂，y₂，z₂)，(x₃，y₃，z₃)，坡度s和坡向αs作为平面固有的几何属性，其数值计算可以使用以下方程(2)和方程(3)计算：

式中：w_x和w_y分别为横轴x和纵轴y方向的坡度，且

上式中，所得到的坡向αs以弧度为单位，正北方向(即正上方)为0rad，顺时针增大；

(2-2)、若坡向不在三角形平面对应的范围内，则将坡向调整为最接近原方向的允许范围边缘方向，坡度为沿新坡向的坡度，设三角形平面允许的角度范围的最大值和最小值分别为α_max和α_min，若方程(4)得到的αs方向不在步骤(2)确定的允许范围内，将排水方向αs调整为最接近原方向的允许范围边缘方向，设原方向为

则最终αs方向的方程为：

此时新的坡度为中心单元C₀跟新坡向指向单元间的高差△h与两个单元中心间距离l的比值，若新坡向指向单元在C₀的上、下、左、右4个方向之一，则l为栅格单元的尺寸△x，若新坡向指向单元在C₀的左上、左下、右上、右下4个方向之一，则l为栅格单元的尺寸△x的

倍；

(2-3)、若中心单元C₀与C₂、C₃组成的三角形平面坡度最大，则将该三角形平面的角度范围(即为图3的灰色区域)作为最终排水方向的允许范围，并记录该平面的坡向为α₁。

作为本发明提供的一种基于平均化处理的DEM栅格局部排水方向确定方法的进一步优化方案，所述步骤(3)具体内容如下：

对于栅格尺寸为△x的DEM，假设z₅为需要求取切面曲率的中心单元，如图4以其为中心的3×3窗口内9个单元的高程编号分别为左上单元z₁、正上单元z₂、右上单元z₃、正左单元z₄、正中单元z₅、正右单元z₆、左下单元z₇、正下单元z₈、右下单元z₉，步骤(3)中切面曲率的计算具体可以采用如下方程：

其中，f_x和f_xx分别为高程函数f(x,y)的x方向的一阶和二阶导数，f_y和f_yy分别为高程函数f(x,y)的y方向的一阶和二阶导数，f_xy为高程函数f(x,y)的x、y交互方向的二阶导数；根据3×3窗口内的高程分布(z₁-z₉)，f_x、f_y、f_xx、f_yy、f_xy的数值可以用如下方程确定：

作为本发明提供的一种基于平均化处理的DEM栅格局部排水方向确定方法的进一步优化方案，所述步骤(4)的具体内容如下：

假设步骤(2)选取的坡度最大三角形平面如图3中灰色区域所示，即步骤(3)选取的下游单元为C₂和C₃(图5中(a)，图中6(a))，C₂和C₃的切面曲率分别为K₂和K₃，若K₂>K₃，则从除C₂和C₃外的C₀的6个邻近单元中新选取最接近C₂的单元C₁，连接C₁，C₃，C₀三点构建三角形平面，由方程(3)确定新平面的坡向α₂；若K₂<K₃，则从除C₂和C₃外的C₀的6个邻近单元中新选取最接近C₃的单元C₅，连接C₂，C₅，C₀三点构建三角形平面，由方程(3)确定新平面的坡向α₂。

作为本发明提供的一种基于平均化处理的DEM栅格局部排水方向确定方法的进一步优化方案，所述步骤(5)具体按以下步骤进行：

(5-1)、排水方向α作为α₁和α₂所夹锐角的角平分线方向(图5中(c)，图6中(c))，其数值可以使用如下方程计算：

(5-2)、设步骤(2)确定的允许的角度范围的最大值和最小值分别为α_max和α_min，若方程(11)得到的α方向不在步骤(2)确定的允许范围内，将排水方向α调整为最接近原方向的允许范围边缘方向，设原方向为

则最终α方向的方程为：

本发明的有益效果为：本发明以需要确定排水方向的DEM栅格单元为中心，结合与其相邻的8个栅格单元，使用2种三角面连接方式构建2个三角形平面，使用2个三角形平面的平均坡向，即2个面各自坡向的角平分线方向作为DEM栅格单元排水方向；本发明借助切面曲率辅助三角形平面的构建，使用两个三角形平面的坡向的平均值作为模拟方向，相较于传统方法，使用本发明得到的排水方向与理论方向间的偏差更小；本发明为基于DEM的水文建模提供了更准确地模拟地表水流过程的方法。

附图说明

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

图2为本发明实施例中8个三角形平面的划分方法示意图。

图3为本发明实施例中允许范围和最陡坡向确定示意图。

图4为本发明实施例中切面曲率计算使用的3×3窗口内的DEM单元标识示意图。

图5为本发明中单元C2的切面曲率较高时排水方向的确定示意图。

图6为本发明中单元C3的切面曲率较高时排水方向的确定示意图。

图7为本发明实施例中使用的DEM栅格局部排水方向确定方法的示意图。

图8为本发明中单元C2的切面曲率较高时排水方向的确定示意图。

具体实施方式

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

参见图1至图8，本发明是：一种基于平均化处理的DEM栅格局部排水方向确定方法，其中，所述方法具体步骤如下：

(1)、对于需要确定排水方向的栅格C₀，令其周边相邻的8个栅格按顺序编号为C₁-C₈，并按顺序连接上述9个栅格的中心点，构建出8个三角形平面；

(2)、比较每个三角形平面的另外两个顶点，剔除另外两个顶点高程均高于中心单元C₀的平面，计算剩下的三角形平面的坡度和坡向，若坡向不在三角形平面对应的范围内，则将坡向调整为最接近原方向的允许范围边缘方向，坡度为沿新坡向的坡度，然后选取坡度最大的三角形平面所对应的45°范围作为最终排水方向的可能范围，将此单元的坡向定义为α₁；

(3)、除C₀外，将步骤(2)中选取的坡度最大三角形平面上的另外2个顶点定义为C₀的2个下游单元，计算2个下游单元的切面曲率，并比较其切面曲率大小；

(4)、除步骤(3)选取的2个下游单元外，从C₀的另外6个相邻单元中选取1个新的单元，该单元最接近具有较高切面曲率的下游单元，使用新选取的单元与C₀、具有较低切面曲率的下游单元构建三角形平面，计算新平面的坡向α₂；

(5)、将单元C₀的排水方向α确定为α₁和α₂所夹锐角的角平分线方向，若该角平分线方向不在步骤(2)确定的允许范围内，则将排水方向α定为最接近该角平分线方向的允许范围边缘方向。

具体地，所述步骤(2)具体按以下步骤进行：

(2-1)、设组成三角形平面的三个DEM栅格单元中心点坐标分别为(x₁，y₁，z₁)，(x₂，y₂，z₂)，(x₃，y₃，z₃)，坡度s和坡向αs作为平面固有的几何属性，其数值计算可以使用以下方程(2)和方程(3)计算：

式中：w_x和w_y分别为横轴x和纵轴y方向的坡度，且

上式中，所得到的坡向αs以弧度为单位，正北方向(即正上方)为0rad，顺时针增大；

(2-2)、若坡向不在三角形平面对应的范围内，则将坡向调整为最接近原方向的允许范围边缘方向，坡度为沿新坡向的坡度，设三角形平面允许的角度范围的最大值和最小值分别为α_max和α_min，若方程(4)得到的αs方向不在步骤(2)确定的允许范围内，将排水方向αs调整为最接近原方向的允许范围边缘方向，设原方向为

则最终αs方向的方程为：

此时新的坡度为中心单元C₀跟新坡向指向单元间的高差△h与两个单元中心间距离l的比值，若新坡向指向单元在C₀的上、下、左、右4个方向之一，则l为栅格单元的尺寸△x，若新坡向指向单元在C₀的左上、左下、右上、右下4个方向之一，则l为栅格单元的尺寸△x的

倍；

(2-3)、若中心单元C₀与C₂、C₃组成的三角形平面坡度最大，则将该三角形平面的角度范围(即为图3的灰色区域)作为最终排水方向的允许范围，并记录该平面的坡向为α₁。

具体地，所述步骤(3)具体内容如下：

对于栅格尺寸为△x的DEM，假设z₅为需要求取切面曲率的中心单元，如图4以其为中心的3×3窗口内9个单元的高程编号分别为左上单元z₁、正上单元z₂、右上单元z₃、正左单元z₄、正中单元z₅、正右单元z₆、左下单元z₇、正下单元z₈、右下单元z₉，步骤(3)中切面曲率的计算具体可以采用如下方程：

其中，f_x和f_xx分别为高程函数f(x,y)的x方向的一阶和二阶导数，f_y和f_yy分别为高程函数f(x,y)的y方向的一阶和二阶导数，f_xy为高程函数f(x,y)的x、y交互方向的二阶导数；根据3×3窗口内的高程分布(z₁-z₉)，f_x、f_y、f_xx、f_yy、f_xy的数值可以用如下方程确定：

具体地，所述步骤(4)的具体内容如下：

假设步骤(2)选取的坡度最大三角形平面如图3中灰色区域所示，即步骤(3)选取的下游单元为C₂和C₃(图5中(a)，图6中(a))，C₂和C₃的切面曲率分别为K₂和K₃，若K₂>K₃，则从除C₂和C₃外的C₀的6个邻近单元中新选取最接近C₂的单元C₁，连接C₁，C₃，C₀三点构建三角形平面，由方程(3)确定新平面的坡向α₂；若K₂<K₃，则从除C₂和C₃外的C₀的6个邻近单元中新选取最接近C₃的单元C₅，连接C₂，C₅，C₀三点构建三角形平面，由方程(3)确定新平面的坡向α₂。

具体地，所述步骤(5)具体按以下步骤进行：

(5-1)、排水方向α作为α₁和α₂所夹锐角的角平分线方向(图5中(c)，图6中(c))，其数值可以使用如下方程计算：

(5-2)、设步骤(2)确定的允许的角度范围的最大值和最小值分别为α_max和α_min，若方程(11)得到的α方向不在步骤(2)确定的允许范围内，将排水方向α调整为最接近原方向的允许范围边缘方向，设原方向为

则最终α方向的方程为：

为了更好地实现本发明目的，本发明以下具体实例进行进一步说明：

以图7所示的栅格DEM为例，计算高程为17.027的栅格单元的排水方向，该DEM每个栅格单元的高程z使用方程

确定，整数x和y分别为横向和纵向坐标，且6≤x≤10，6≤y≤10；当坐标系以正上方为0rad，角度顺时针增大时，根据高程的定义方程，该DEM上所有栅格的理想流向都指向π/2+arctan(1/5)rad(约2.034rad)；

步骤(1)：以DEM中需要计算排水方向的高程为17.027的单元为中心单元C₀，连接中心单元与周边8个单元的中心点，划分8个三角形平面，如图8中(a)；

步骤(2)：比较每个三角形平面的另外两个顶点，剔除另外两个顶点高程均高于中心单元C₀的平面，随后计算剩下的三角形平面的坡度和坡向，若坡向不在三角形平面对应的范围内，则将坡向调整为最接近原方向的允许范围边缘方向，坡度为沿新坡向的坡度，除被剔除的三角形平面外剩下的三角形平面的坡度值如图8中(b)所示，选取坡度值最大的三角形平面，然后选取坡度最大的三角形平面所对应的范围作为最终排水方向的可能范围，即图8中(b)中坡度为0.899的灰色三角形平面，将此单元的坡向定义为α₁，即α₁＝1.761rad；

步骤(3)：根据方程(5)计算图8中(b)中灰色三角形平面的2个下游单元的切面曲率，得到高程为16.144的单元切面曲率为-0.0134，高程为15.974的单元切面曲率为-0.0133；

步骤(4)：由于步骤(3)中高程为15.974的切面曲率较大，因此从中心单元的另外6个相邻单元中选取最接近高程15.974的单元，即高程为16.846的单元，使用该单元与高程为17.027的中心单元、高程为16.144的单元构建三角形平面(图8中(c))，由方程(3)求得新平面的坡向α₂＝2.558rad；

步骤(5)：最终确定的高程为17.027的栅格单元的排水方向α为α₁和α₂所夹锐角的角平分线方向，根据方程(11)得到α＝2.160rad。

对于图7中的实例，传统的D8方法将中心单元的排水方向指向高程为16.144的单元的中心点，得到排水方向数值为1.571rad；Dinf方法只接受图2中的固定三角面划分方法，使用最陡三角面的坡向作为排水方向(即等同于本发明中的α₁)，得到排水方向数值为1.761rad；相比D8和Dinf这两种传统方法，使用本发明得到的2.160rad排水方向更接近栅格处2.034rad的理论方向。

综上所述，本发明借助切面曲率辅助三角形平面的构建，使用两个三角形平面的坡向的平均值作为模拟方向，发明了一种确定DEM栅格单元局部排水方向的方法，相较于传统方法，使用本发明得到的排水方向与理论方向间的偏差更小。

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

Claims (5)

1.一种基于平均化处理的DEM栅格局部排水方向确定方法，其特征在于，所述方法具体步骤如下：

(1)、对于需要确定排水方向的栅格C₀，令其周边相邻的8个栅格按顺序编号为C₁-C₈，并按顺序连接上述9个栅格的中心点，构建出8个三角形平面；

(2)、比较每个三角形平面的另外两个顶点，剔除另外两个顶点高程均高于中心单元C₀的平面，计算剩下的三角形平面的坡度和坡向，若坡向不在三角形平面对应的范围内，则将坡向调整为最接近原方向的允许范围边缘方向，坡度为沿新坡向的坡度，然后选取坡度最大的三角形平面所对应的45°范围作为最终排水方向的可能范围，将此单元的坡向定义为α₁；

(3)、除C₀外，将步骤(2)中选取的坡度最大三角形平面上的另外2个顶点定义为C₀的2个下游单元，计算2个下游单元的切面曲率，并比较其切面曲率大小；

(4)、除步骤(3)选取的2个下游单元外，从C₀的另外6个相邻单元中选取1个新的单元，该单元最接近具有较高切面曲率的下游单元，使用新选取的单元与C₀、具有较低切面曲率的下游单元构建三角形平面，计算新平面的坡向α₂；

(5)、将单元C₀的排水方向α确定为α₁和α₂所夹锐角的角平分线方向，若该角平分线方向不在步骤(2)确定的允许范围内，则将排水方向α定为最接近该角平分线方向的允许范围边缘方向。

2.根据权利要求1所述的基于平均化处理的DEM栅格局部排水方向确定方法，其特征在于，所述步骤(2)具体按以下步骤进行：

(2-1)、设组成三角形平面的三个DEM栅格单元中心点坐标分别为(x₁，y₁，z₁)，(x₂，y₂，z₂)，(x₃，y₃，z₃)，坡度s和坡向αs作为平面固有的几何属性，其数值计算使用如下方程(2)和方程(3)计算：

式中：w_x和w_y分别为横轴x和纵轴y方向的坡度，且

上式中，所得到的坡向αs以弧度为单位，正北方向为0rad，顺时针增大；

(2-2)、若坡向不在三角形平面对应的范围内，则将坡向调整为最接近原方向的允许范围边缘方向，坡度为沿新坡向的坡度，设三角形平面允许的角度范围的最大值和最小值分别为α_max和α_min，若方程(4)得到的αs方向不在步骤(2)确定的允许范围内，将排水方向αs调整为最接近原方向的允许范围边缘方向，设原方向为

则最终αs方向的方程为：

此时新的坡度为中心单元C₀跟新坡向指向单元间的高差△h与两个单元中心间距离l的比值，若新坡向指向单元在C₀的上、下、左、右4个方向之一，则l为栅格单元的尺寸△x，若新坡向指向单元在C₀的左上、左下、右上、右下4个方向之一，则l为栅格单元的尺寸△x的

倍；

(2-3)、若中心单元C₀与C₂、C₃组成的三角形平面坡度最大，则将该三角形平面的角度范围作为最终排水方向的允许范围，并记录该平面的坡向为α₁。

3.根据权利要求2所述的基于平均化处理的DEM栅格局部排水方向确定方法，其特征在于，所述步骤(3)具体内容如下：

对于栅格尺寸为△x的DEM，假设z₅为需要求取切面曲率的中心单元，以其为中心的3×3窗口内9个单元的高程编号分别为左上单元z₁、正上单元z₂、右上单元z₃、正左单元z₄、正中单元z₅、正右单元z₆、左下单元z₇、正下单元z₈、右下单元z₉，步骤(3)中切面曲率的计算具体采用如下方程：

其中，f_x和f_xx分别为高程函数f(x,y)的x方向的一阶和二阶导数，f_y和f_yy分别为高程函数f(x,y)的y方向的一阶和二阶导数，f_xy为高程函数f(x,y)的x、y交互方向的二阶导数；根据3×3窗口内的高程分布(z₁-z₉)，f_x、f_y、f_xx、f_yy、f_xy的数值用如下方程确定：

4.根据权利要求3所述的基于平均化处理的DEM栅格局部排水方向确定方法，其特征在于，所述步骤(4)的具体内容如下：

假设步骤(2)选取的坡度最大三角形平面，即步骤(3)选取的下游单元为C₂和C₃，C₂和C₃的切面曲率分别为K₂和K₃，若K₂>K₃，则从除C₂和C₃外的C₀的6个邻近单元中新选取最接近C₂的单元C₁，连接C₁，C₃，C₀三点构建三角形平面，由方程(3)确定新平面的坡向α₂；若K₂<K₃，则从除C₂和C₃外的C₀的6个邻近单元中新选取最接近C₃的单元C₅，连接C₂，C₅，C₀三点构建三角形平面，由方程(3)确定新平面的坡向α₂。

5.根据权利要求4所述的基于平均化处理的DEM栅格局部排水方向确定方法，其特征在于，所述步骤(5)具体按以下步骤进行：

(5-1)、排水方向α作为α₁和α₂所夹锐角的角平分线方向，其数值使用如下方程计算：

(5-2)、设步骤(2)确定的允许的角度范围的最大值和最小值分别为α_max和α_min，若方程(11)得到的α方向不在步骤(2)确定的允许范围内，将排水方向α调整为最接近原方向的允许范围边缘方向，设原方向为

则最终α方向的方程为：

Images