CN113763553B - 一种基于数字高程模型的小流域主沟道提取方法 - Google Patents

Abstract

本发明公开了一种基于数字高程模型的小流域主沟道提取方法，主要基于数字高程模型，通过自动捕捉或指定沟源后，采用逐点依次向下自动追索的方法确定主沟道，同时计算并保存主沟道的位置和高程信息，用于后续的数据处理；并在出现局部低洼区时进行异常处理，以最低点为中心扩大搜索范围，直至寻找到低于低洼区最低点，跳出低洼区。本发明采用自动追索的方法，具有计算量小，计算简单，提取的主沟道位置精确，便于后期进一步处理和分析等优点。

Description

一种基于数字高程模型的小流域主沟道提取方法

技术领域

本发明涉及沟道提取技术领域，具体为一种基于数字高程模型的小流域主沟道提取方法。

背景技术

沟道是水文信息的基础要素，准确地提取流域中沟道信息，是水文模型的构建、流域的划分等的必要条件。目前提取沟道通常采用河网识别法。

小流域主沟道提取的方法主要有两种，一种是人工目视描绘提取，这种方法通常基于地形图人工判断描绘沟道，费时费力。另一种是通过专业的水文分析提取，水文分析专业性强，步骤繁琐，计算量大。更重要的是，这两种方法提取的主沟道都只有平面位置信息，没有高程信息，不利于后期的数据自动处理，如生成主沟道纵剖面图，计算沟道加权平均坡降等。

发明内容

针对上述问题，本发明的目的在于提供一种基于数字高程模型的小流域主沟道提取方法，主要基于数字高程模型(下文统一简称“DEM”)，计算简单，提取的主沟道位置精确，便于后期进一步处理和分析。技术方案如下：

一种基于数字高程模型的小流域主沟道提取方法，包括以下步骤：

步骤1：数据预处理

假设流域的DEM已经被正确提取，且DEM共有m行、n列个栅格，每个栅格有一个数值，写成m×n的矩阵G：

其中，流域界线以内的每个栅格都有一个高程值，流域界线以外的栅格在提取时被指定为无数据栅格。

步骤2：计算流域最高点和最低点

计算流域最高点时，把矩阵中无数据栅格的值统一换成很小值，令

A＝{f(i,j)|f(i,j)＝g_ij,1≤i≤m,1≤j≤n}，如果

对于/>

都有

f(r,c)≥f(i,j) (2)

则最高点在DEM中的位置为(r,c)，高程值为g_rc；r、c分别是最高点所在栅格的行号和列号，1≤r≤m,1≤c≤n；

计算流域最低点时，把矩阵中的很小值，统一换成很大值，令

B＝{f(i,j)|f(i,j)＝g_ij,1≤i≤m,1≤j≤n}，如果

对/>

都有

f(p,q)≤f(i,j) (3)

则最低点在DEM中的位置为(p,q)，高程为g_pq；p、q分别是最低点所在栅格的行号和列号，1≤p≤m,1≤q≤n；

步骤3：计算沟源位置

如果沟源是流域最高点，直接以流域内高程最大值所在位置作为主沟道的沟源；如果不是，则指定沟源所在区域，以区域内最高点作为沟源；

步骤4：确定主沟道

主沟道采用逐点追索的方法，以沟源为起始点，记录起始点的位置和高程信息；以起始点为中心寻找其周围8个栅格点中高程最低者作为下一个中心点，同时记录该点的位置和高程信息；然后再以新的中心点为中心寻找其周围8个栅格点中高程最低者作为再下一个中心点，并记录其位置和高程信息；依次向下追索，直至主沟道终点，即流域的最低点；

步骤5：计算主沟道平面位置

基于单个栅格所代表的经纬度值和DEM左上角的坐标，首先计算主沟道上栅格的经纬度坐标，并确定各栅格在平面上的位置，然后依次连接各栅格，得到主沟道在平面上的位置。

进一步的，所述步骤3计算沟源位置具体为：

定义关于沟道上栅格的位置(i,j)和顺序值s的函数H(i,j,s)＝g_ij，用来表示沟道上的栅格位置及该栅格的顺序值与高程的映射关系；s为正整数，并令沟源位置的s＝1，以后各点依次递增；

如果沟源是流域的最高点，即g_gy＝g_rc，则以此栅格为沟源位置，沟源记为H(g,y,1)＝g_gy，位置为(g,y)，是沟道上第1个点，高程值为g_gy；

如果沟源不在流域最高点，则以行、列号的方式指定沟源所在矩形区域，计算区域内的最大值，即区域内的最高点，并以此栅格点为沟源位置；

令C＝{f(i,j)|f(i,j)＝g_ij,ru≤i≤rl,cl≤j≤cr}，如果

对/>

都有

f(g,y)≥f(i,j) (4)

则沟源栅格的位置为(g,y)，高程为g_gy，记H(g,y,1)＝g_gy；g、y分别是沟源所在栅格的行号和列号，且ru≤g≤rl,cl≤y≤cr；ru为指定的矩形区域顶行行号，rl为指定的矩形区域底行行号；cl为指定的矩形区域左边界列号，cr为指定的矩形区域右边界列号。

更进一步的，所述步骤4确定主沟道的具体方法为：

定义点集L，令已提取的沟道上的栅格点(i,j)∈L，以沟源点g_ij为中心，寻找其周围8个栅格的高程最小值g_kl；

令D＝{f(u,v)|f(u,v)＝g_uv,(i-1)≤u≤(i+1),(j-1)≤v≤(j+1),且f(u,v)的定义域不包括(i,j)}，如果

对/>

都有

f(k,l)≤f(u,v) (5)

其中，k、l分别是高程最小值所在栅格点的行号和列号，(i-1)≤k≤(i+1),(j-1)≤l≤(j+1)，k≠i，l≠j；

对于得到的点(k,l)判断其是否在点集L中：

当(k,l)∈L时，表明该点是沟道上的点，且被重复确定为中心点；为避免该点再次成为中心点，令

g′_kl＝g_kl+Δ (6)

其中，g′_kl为该点修正后的高程值，Δ为修正量；

当

时，记(k,l)∈L，H(k,l,s)＝g_kl；

重复以上步骤，寻找主沟道上的下一个点，直至主沟道终点。

更进一步的，所述步骤4中，在追索主沟道终点的过程中，若存在无法继续追索下去的低洼区的最低点，则通过异常处理算法，以该最低点为中心扩大搜索范围，直至寻找到低于低洼区最低点的点，再以此点为中心，继续向下游追索。

更进一步的，所述异常处理算法具体为：

令E_b＝{f(h,f)|f(h,f)＝g_hf,(i-b)≤h≤(i+b),(j-b)≤f≤(j+b),b∈N^*}，

当b＝1时，E₁为不存在异常的集合，此时E＝E₁；

当b>1时，E＝E_b-E_b-1，且

对/>

都有

f(k,l)≤f(h,f) (7)

如果g_kl<g_ij，记(k,l)∈L，H(k,l,s)＝g_kl；否则，令b加1并重复以上步骤直至g_kl<g_ij。

更进一步的，所述步骤5具体为：

所述单个栅格所代表的经纬度值通过整个DEM左右边界经度差，除以栅格的总列数得到

L_raster＝(L_r-L_l)/n (8)

其中，L_raster为单个栅格所代表的经度宽度或纬度高度值；L_r为DEM右边界的经度值；

L_l为DEM左边界的经度值；n为DEM的栅格总列数；

假设DEM左上角坐标为(lon,lat)，即g₁₁栅格的经度为lon，纬度为lat，则栅格g_ij的地理坐标为

L_i＝lat-L_raster×(i-1)

L_j＝lon+L_raster×(j-1) (9)

其中，L_i为栅格g_ij的纬度值；L_j为栅格g_ij的经度值；

通过上式求得主沟道上每个栅格点H(k,l,s)的经纬度，根据s值依次连接各点，所得曲线即为主沟道平面位置图。

本发明的有益效果是：本发明提出的小流域主沟道自动提取算法，主要基于数字高程模型，通过自动捕捉或指定沟源后自动追索主沟道，同时计算并保存主沟道的位置和高程信息，用于后续的数据处理。该算法采用自动追索的方法，具有计算量小，计算简单，提取的主沟道位置精确，便于后期进一步处理和分析等优点。

附图说明

图1为小流域主沟道提取流程图。

图2为栅格g_ij周围8个栅格分布位置示意图。

图3为扩大搜索范围时参与计算最小值的栅格示意图；(a)扩大搜索范围时参与计算的栅格(灰色部分)(b)b＝2时参与计算的栅格(灰色部分)。

图4为瓦格沟影像图(a)和瓦格沟数字高程模型(DEM)图(b)。

图5为根据最高点(a)和指定沟源区域(b)追索的主沟道图，为展示方便,图中行列号已经换算成地理坐标。

图6为利用主沟道数据集自动绘制的瓦格沟流域主沟道纵剖面图。

具体实施方式

下面结合附图和具体实施例对本发明做进一步详细说明。

1、计算思路

小流域主沟道提取主要涉及沟源位置的确定、沟道提取、沟道终点的确定和异常处理几部分算法。

(1)沟源位置的确定：流域主沟道沟源大多是流域最高点，但也可能不是。如果是流域最高点，直接以流域内高程最大值所在位置作为主沟道的沟源；如果不是，需要指定沟源所在区域(无需指定具体位置，给出一个包含沟源在内的区域即可)，以区域内最高点作为沟源。

(2)主沟道的确定：主沟道采用逐点追索的方法确定。以沟源为起始点，记录起始点的位置和高程信息；以起始点为中心寻找其周围8个栅格点中高程最低者作为下一个中心点，同时记录该点的位置和高程信息；然后再以新的中心点为中心寻找其周围8个栅格点中高程最低者作为再下一个中心点，并记录其位置和高程信息；依次向下追索，直至主沟道终点。

(3)主沟道终点的确定：主沟道终点通常定在主沟道与主河交汇处，是流域的最低点。当主沟道追索到流域最低点时即可结束，并以此点作为主沟道的终点。

(4)异常处理：在主沟道上有时会出现局部低洼区，更多的是由于DEM数据质量问题造成的，一旦追索到低洼区的最低点，很难“跃”出低洼处，导致无法继续追索下去，此时就需要进行异常处理。处理方法是以最低点为中心扩大搜索范围，直至寻找到低于低洼区最低点的点，再以此点为中心，继续向下游追索。

本发明的算法是基于DEM，文中的“点”代表DEM的“栅格”，一个点代表一个栅格，点的位置即为栅格中心的位置。

具体的提取思路和流程详见图1。

2、数据预处理

假设流域的DEM已经被正确提取，且DEM共有m行、n列个栅格，每个栅格有一个数值，可写成m×n的矩阵G(m,n)

矩阵G(m,n)中数值可能是高程数据，也可能不是有效的高程数据。流域界线以内的每个栅格都有一个高程值；流域界线以外的栅格在提取时被指定为无数据(NoData)栅格，通常用一个很大或很小的值表示，如32767或-32767。

流域内最高点通常是沟源，最低点是主沟与主河交汇处，需要首先计算出来。

2.1最高点的算法

为了计算流域最高点的信息，需要把矩阵中很大值(无数据栅格的值，也许没有很大值)，统一换成很小值(如-32767)，即

如果g_ij>9000，则令g_ij＝-32767。

采用以下算法寻找流域最高点

令A＝{f(i,j)|f(i,j)＝g_ij,1≤i≤m,1≤j≤n}，如果

对于/>

都有

f(r,c)≥f(i,j) (2)

则最高点在DEM中的位置为(r,c)，高程为g_rc。

这里，r、c分别是最高点所在栅格的行号和列号。

DEM的每个栅格都有一个值(也许是有效高程，也许不是)，当栅格的值大于珠峰高度(8848m)或低于死海高程(-400m)时，就不是有效的高程。为了避免错误通常会用很大或很小的值来表示无效高程，这样以来当DEM栅格的值很大(如32767)或很小(-如32767)时，就被视为“无数据(NoData)栅格”。因此，为了寻找流域内最高点，就需要把流域界线以外的很大值换成很小值(负值)，这样矩阵G中的最大值就是流域最高点。流域最低点采用相反的处理方法来确定。

2.2最低点的算法

为了计算流域最低点的信息，需要把矩阵中的小值，统一换成很大值(如32767)，即

如果g_ij≤-32767，则令g_ij＝32767。

采用以下算法寻找流域最低点

再令B＝{f(i,j)|f(i,j)＝g_ij,1≤i≤m,1≤j≤n}，如果

对/>

都有

F(p,q)≤f(i,j) (3)

则最低点在DEM中的位置为(p,q)高程为g_pq。

这里，p、q分别是最低点所在栅格的行号和列号。

3、核心算法

3.1沟源位置的算法

现定义一个函数H(i,j,s)＝g_ij，该函数用来表示沟道上的栅格位置及该栅格的顺序值与高程的(映射)关系。这里，H(i,j,s)是关于沟道上栅格的位置(i,j)和顺序值s的函数，s为正整数并令沟源位置的s＝1，以后各点依次递增；g_ij为栅格(i,j)的高程值。

如果沟源(g_gy)是流域的最高点，即g_gy＝g_rc，则以此栅格为沟源位置。沟源记为H(g,y,1)＝g_gy，位置为(g,y)，是沟道上第1个点，高程值为g_gy。

如果沟源不在流域最高点，则需要以行、列号的方式指定沟源所在矩形区域，计算区域内的最大值，即区域内的最高点，并以此栅格点为沟源位置。

令C＝{f(i,j)|f(i,j)＝g_ij,ru≤i≤rl,cl≤j≤cr}，如果

对/>

都有

f(g,y)≥f(i,j) (4)

则沟源栅格的位置为(g,y)，高程为g_gy，记H(g,y,1)＝g_gy。

这里，g、y分别是沟源所在栅格的行号和列号；

ru-指定的矩形区域顶行行号，rl-指定的矩形区域底行行号；

cl-指定的矩形区域左边界列号，cr-指定的矩形区域右边界列号。

其中，沟源的判断主要依赖于人的经验，根据流域和沟道特征综合确定。通常情况下，以最长的沟道，或者从流域最高点发育下来的沟道作为主沟道。实际上，通过指定沟源所在区域这一功能，可以把整个流域内所有的支沟都能提取出来。

3.2主沟道提取算法

现定义一个点集L，令已提取的沟道上的栅格点(i,j)∈L。假设以g_ij为中心，寻找其周围8个栅格(图2)的高程最小值g_kl。

令D＝{f(u,v)|f(u,v)＝g_uv,(i-1)≤u≤(i+1),(j-1)≤v≤(j+1),且f(u,v)的定义域不包括(i,j)}，如果

对/>

都有

f(k,l)≤f(u,v) (5)

这里，k、l分别是高程最小值所在栅格点的行号和列号。

当(k,l)∈L时，表明该点是沟道上的点，且被重复确定为中心点(最小值)。为了避免该点再次成为中心点，令

g′_kl＝g_kl+Δ (6)

g′_kl为该点修正后的高程值，Δ为修正量

当

时，记(k,l)∈L，H(k,l,s)＝g_kl。

重复以上步骤，寻找主沟道上的下一个点，直至主沟道终点。

当寻找到的最低点已经存在于点集L中时，该点与上一个中心点(上一个最小值)会被交替确定为中心点，且一直交替循环下去。为避免这种情况，利用公式(6)把该点加上一个数值(比如100)，使其不再是局部的最小值。同时，该点不再重复记入点集L中。如果寻找到的最低点不在点集L中，则把该点记入点集。

3.3异常处理算法

当中心点落入低洼区，即寻找到的最小值(g_kl)大于中心点的值(g_ij)时，则需要通过扩大搜索范围的办法“跳”出低洼区。

令E_b＝{f(h,f)|f(h,f)＝g_hf,(i-b)≤h≤(i+b),(j-b)≤f≤(j+b),b∈N^*}，

当b＝1时，E₁为不存在异常的集合，此时E＝E₁。

当b>1时，E＝E_b-E_b-1，(图3(a))，且

对/>

都有

f(k,l)≤f(h,f) (7)

如果g_kl<g_ij，记(k,l)∈L，H(k,l,s)＝g_kl。否则，令b加1并重复以上步骤直至g_kl<g_ij。

例如，当b＝2时(图3(b))，E＝E₂-E₁，则

对/>

都有

f(k,l)≤f(h,f) (8)

如果g_kl<g_ij，记(k,l)∈L，H(k,l,s)＝g_kl。

否则，令b＝3，…，重复上一步，直至g_kl<g_ij。

这里，k、l分别是高程最小值所在栅格点的行号和列号。

3.4主沟道平面位置算法

在地理坐标系统下，DEM像素单元为正方形栅格，即单个栅格在经度方向和纬度方向代表的长度是一样的。因此，单个栅格所代表的经纬度值可通过整个DEM左右边界经度差，除以栅格的总列数得到

L_raster＝(L_r-L_l)/n (9)

这里，L_raster-单个栅格所代表的经度宽度或纬度高度值(°)；

L_r-DEM右边界的经度值(°)；

L_l-DEM左边界的经度值(°)；

n-DEM的栅格总列数。

假设DEM左上角坐标为(lon,lat)，即g₁₁栅格的经度为lon，纬度为lat，则栅格g_ij的地理坐标为

L_i＝lat-L_raster×(i-1)

L_j＝lon+L_raster×(j-1) (10)

这里，L_i-栅格g_ij的纬度值(°)；

L_j-栅格g_ij的经度值(°)。

通过式(10)可求得主沟道上每个栅格点H(k,l,s)的经纬度，根据s值依次连接各点，所得曲线即为主沟道平面位置图。

“单个栅格所代表的经纬度值”是一个计量单位，可理解成单个格栅的长和宽。如果要把主沟道绘制出来，首先需要确定每个栅格的经纬度坐标(不是栅格的长和宽)，然后依次连接各点可得到主沟道平面位置图。

4、应用实例

本发明以瓦格沟为实例，进一步说明该发明的原理和算法。瓦格沟为金沙江上游右岸的一条一级支沟(图4)，行政区划上属西藏自治区昌都市江达县汪布顶乡，流域面积约37.37km²，沟口地理坐标为98°17'8.173"E，32°9'51.806"N。

4.1数据预处理

瓦格沟流域的DEM见图4(b)，DEM占据的矩形区域共有285行、335列，各栅格点的高程可写成矩阵G(285,335)

矩阵G中流域界线以内的每个栅格都有一个高程值；流域界线以外的栅格在提取时被指定为无数据(NoData)栅格，即-32767。

根据计算，流域内最高点高程为4964m，栅格位于148行，12列；最低点高程为3140m，栅格位于第6行，333列。最低点为主沟与主河交汇处，即主沟道终点。

4.2、沟源位置的计算

假设以该流域的最高点作为沟源，则沟源为H(148,12,1)＝4964，位置为(148,12)，是沟道上第1个点，高程值为4964m。以此点作为沟源计算得到的主沟道如图5(a)所示。

由于该流域沟源并不在最高点，需要指定沟源所在区域。根据影像和DEM可确定流域沟源大致位于顶行号(ru)270，底行号(rl)285；左界列号(cl)40，右界列号(cr)70所限定的区域，见图5(b)。通过计算，实际沟源为H(282,65,1)＝4759，位置为(282,65)，是沟道上第1个点，高程值为4759m。

4.3主沟道提取

通过指定沟源所在区域，计算并确定了沟源位置为(282,65)(图5)；然后，以此点为起点，根据式(5)依次向下游追索主沟道上的各点，直到追索至H(6,333,349)＝3140处，即流域最低点结束，结果见图5(b)。

4.4异常处理

当中心点落入低洼区，即寻找到的最小值大于中心点的值时，则需要通过扩大搜索范围的办法“跳”出低洼区。如图5(a)，当从最高点追索主沟道时，中游沟道H(140,188,186)＝3982处开始出现“低洼区”，利用式(7)的计算方法“跳”出低洼区至H(129,202,187)＝3981后，继续向下游追索。

4.5主沟道平面位置计算

瓦格沟流域的DEM共有栅格285行，335列。左上角的地理坐标为(98.1936590778°E,32.1655015497°N)，右下角为(98.2867146334°E,32.0863348831°N)。根据式(9)计算得到单个栅格所代表的经度宽度和纬度高度值(L_raster)均为0.0002778°。

以DEM左上角坐标为计算参考点，根据式(10)计算得到主沟道上每个栅格点H(k,l,s)的经纬度；根据s依次连接各点，得到如图5(a)所示的主沟道平面位置图。

4.6主沟道剖面图计算

主沟道数据集包含沟道上各点完整的位置和高程信息，除了能够自动生成主沟道平面位置图，还可用于后续的数据处理，如自动绘制沟道纵剖面图等。

首先需要根据该流域所在地理(坐标)位置及其投影参数，分别计算单个栅格经纬度方向的长度(该计算涉及地理坐标“弧长”的算法，是一种成熟的计算方法且与本发明算法无关，具体计算方法不再罗列)。根据瓦格沟流域所在区域的地理坐标和DEM的分辨率，计算得到该流域单个栅格经度方向所代表的长度为26.1580265m，纬度方向所代表的长度为30.8874796m。

利用主沟道函数H(i,j,s)构建的数据集，根据单个栅格经纬度方向所代表的长度，计算得到瓦格沟流域主沟长12.29km。以沟源为起点，利用主沟道数据集自动绘制得到该流域主沟道纵剖面图，如图6所示。

4.7沟床加权平均坡降的计算

利用主沟道数据集还可以计算沟床坡降：平均沟床坡降和加权平均沟床坡降。

平均沟床坡降直接通过主沟道最高点和最低点的高差除以主沟长度得到。瓦格沟主沟沟源高程4759m，主沟与主河交汇处高程为3140m，主沟道长度为12.29km；根据计算，平均沟床坡降为131.73‰。

加权平均沟床坡降的计算，既可以利用主沟道数据集计算沟道上任意两个栅格间的坡度，也可以根据需要计算任意高差(或水平距离)的沟床坡度，然后再进行加权平均计算。本实例以沟源为起点，计算了高差为20m的主沟道上各段的沟床坡度，然后计算了沟床加权(高差)平均坡降。根据计算，瓦格沟主沟床加权平均沟床坡降为143.31‰。

Claims (6)

1.一种基于数字高程模型的小流域主沟道提取方法，其特征在于，包括以下步骤：

步骤1：数据预处理

假设流域的DEM已经被正确提取，且DEM共有m行、n列个栅格，每个栅格有一个数值，写成m×n的矩阵G：

其中，流域界线以内的每个栅格都有一个高程值，流域界线以外的栅格在提取时被指定为无数据栅格；

步骤2：计算流域最高点和最低点

计算流域最高点时，把矩阵中无数据栅格的值统一换成很小值，令A＝{f(i,j)|f(i,j)＝g_ij,1≤i≤m,1≤j≤n}，如果

对于/>

都有

f(r,c)≥f(i,j) (2)

则最高点在DEM中的位置为(r,c)，高程值为g_rc；r、c分别是最高点所在栅格的行号和列号，1≤r≤m,1≤c≤n；

计算流域最低点时，把矩阵中的很小值，统一换成很大值，令B＝{f(i,j)|f(i,j)＝g_ij,1≤i≤m,1≤j≤n}，如果

对/>

都有

f(p,q)≤f(i,j) (3)

则最低点在DEM中的位置为(p,q)，高程为g_pq；p、q分别是最低点所在栅格的行号和列号，1≤p≤m,1≤q≤n；

步骤3：计算沟源位置

如果沟源是流域最高点，直接以流域内高程最大值所在位置作为主沟道的沟源；如果不是，则指定沟源所在区域，以区域内最高点作为沟源；

步骤4：确定主沟道

主沟道采用逐点追索的方法，以沟源为起始点，记录起始点的位置和高程信息；以起始点为中心寻找其周围8个栅格点中高程最低者作为下一个中心点，同时记录该点的位置和高程信息；然后再以新的中心点为中心寻找其周围8个栅格点中高程最低者作为再下一个中心点，并记录其位置和高程信息；依次向下追索，直至主沟道终点，即流域的最低点；

步骤5：计算主沟道平面位置

基于单个栅格所代表的经纬度值和DEM左上角的坐标，首先计算主沟道上栅格的经纬度坐标，并确定各栅格在平面上的位置，然后依次连接各栅格，得到主沟道在平面上的位置。

2.根据权利要求1所述的基于数字高程模型的小流域主沟道提取方法，其特征在于，所述步骤3计算沟源位置具体为：

定义关于沟道上栅格的位置(i,j)和顺序值s的函数H(i,j,s)＝g_ij，用来表示沟道上的栅格位置及该栅格的顺序值与高程的映射关系；s为正整数，并令沟源位置的s＝1，以后各点依次递增；

如果沟源是流域的最高点，即g_gy＝g_rc，则以此栅格为沟源位置，沟源记为H(g,y,1)＝g_gy，位置为(g,y)，是沟道上第1个点，高程值为g_gy；

如果沟源不在流域最高点，则以行、列号的方式指定沟源所在矩形区域，计算区域内的最大值，即区域内的最高点，并以此栅格点为沟源位置；

令C＝{f(i,j)|f(i,j)＝g_ij,ru≤i≤rl,cl≤j≤cr}，如果

对/>

都有

f(g,y)≥f(i,j) (4)

则沟源栅格的位置为(g,y)，高程为g_gy，记H(g,y,1)＝g_gy；g、y分别是沟源所在栅格的行号和列号，且ru≤g≤rl,cl≤y≤cr；ru为指定的矩形区域顶行行号，rl为指定的矩形区域底行行号；cl为指定的矩形区域左边界列号，cr为指定的矩形区域右边界列号。

3.根据权利要求1所述的基于数字高程模型的小流域主沟道提取方法，其特征在于，所述步骤4确定主沟道的具体方法为：

定义点集L，令已提取的沟道上的栅格点(i,j)∈L，以沟源点g_ij为中心，寻找其周围8个栅格的高程最小值g_kl；

令D＝{f(u,v)|f(u,v)＝g_uv,(i-1)≤u≤(i+1),(j-1)≤v≤(j+1),且f(u,v)的定义域不包括(i,j)}，如果

对/>

都有

f(k,l)≤f(u,v) (5)

其中，k、l分别是高程最小值所在栅格点的行号和列号，(i-1)≤k≤(i+1),(j-1)≤l≤(j+1)，k≠i，l≠j；

对于得到的点(k,l)判断其是否在点集L中：

当(k,l)∈L时，表明该点是沟道上的点，且被重复确定为中心点；为避免该点再次成为中心点，令

g′_kl＝g_kl+Δ (6)

其中，g′_kl为该点修正后的高程值，Δ为修正量；

当

时，记(k,l)∈L，H(k,l,s)＝g_kl；

重复以上步骤，寻找主沟道上的下一个点，直至主沟道终点。

4.根据权利要求1所述的基于数字高程模型的小流域主沟道提取方法，其特征在于，所述步骤4中，在追索主沟道终点的过程中，若存在无法继续追索下去的低洼区的最低点，则通过异常处理算法，以该最低点为中心扩大搜索范围，直至寻找到低于低洼区最低点的点，再以此点为中心，继续向下游追索。

5.根据权利要求4所述的基于数字高程模型的小流域主沟道提取方法，其特征在于，所述异常处理算法具体为：

令E_b＝{f(h,f)|f(h,f)＝g_hf,(i-b)≤h≤(i+b),(j-b)≤f≤(j+b),b∈N^*}，

当b＝1时，E₁为不存在异常的集合，此时E＝E₁；

当b>1时，E＝E_b-E_b-1，且

对/>

都有

f(k,l)≤f(h,f) (7)

如果g_kl<g_ij，记(k,l)∈L，H(k,l,s)＝g_kl；否则，令b加1并重复以上步骤直至g_kl<g_ij。

6.根据权利要求1所述的基于数字高程模型的小流域主沟道提取方法，其特征在于，所述步骤5具体为：

所述单个栅格所代表的经纬度值通过整个DEM左右边界经度差，除以栅格的总列数得到

L_raster＝(L_r-L_l)/n (8)

其中，L_raster为单个栅格所代表的经度宽度或纬度高度值；L_r为DEM右边界的经度值；

L_l为DEM左边界的经度值；n为DEM的栅格总列数；

假设DEM左上角坐标为(lon,lat)，即g₁₁栅格的经度为lon，纬度为lat，则栅格g_ij的经纬度为

L_i＝lat-L_raster×(i-1)

L_j＝lon+L_raster×(j-1) (9)

其中，L_i为栅格g_ij的纬度值；L_j为栅格g_ij的经度值；

通过上式求得主沟道上每个栅格点H(k,l,s)的经纬度，根据s值依次连接各点，所得曲线即为主沟道平面位置图。

Images