CN103729552B - 一种自动识别流域数字化过程中出口点的方法 - Google Patents

Abstract

本发明公开了一种自动识别流域数字化过程中出口点的方法，包括如下步骤：读取研究区域的数字高程模型的原始数据并进行预处理；对水文站点在河道栅格矩阵内的位置进行定位，并判断水文站点是否位于河道栅格上；根据栅格控制面积矩阵，提取流域潜在出口点的控制面积，并计算该面积与水文站集水面积之间的相对误差；根据设定的面积许可误差，判断的绝对值是否小于该面积许可误差；根据流向矩阵，对流域的潜在出口点进行修正；识别过程结束。本发明可以对流域潜在出口点进行自动识别和自动修正，实现对流域出口点进行识别的自动化处理，保证了出口点的定位精度，同时避免了人为干预，提高了流域数字化的效率。

Description

一种自动识别流域数字化过程中出口点的方法

技术领域

本发明涉及流域数字化，尤其是识别流域数字化过程中出口点的方法。

背景技术

为了更好地描述水文现象、探索和揭示水文规律，水文学研究正逐步迈向数字化，近年来，随着大范围、高精度数字高程模型（DEM）的出现，数字水文学的研究得到了进一步发展。利用DEM对流域进行数字化也就成为了构建流域水文模型、进行水文模拟与研究的前提。

DEM原始数据是进行流域数字化的基础数据，在利用该基础数据进行流域数字化时，其关键是如何精确、便捷地识别出流域出口点。流域出口点的准确定位对数字化结果有着重要影响，在流域数字化过程中，即使是位置相近的几个点，也会生成完全不同的研究流域边界和集水区域。因此，如果流域出口点定位不准确，不仅限制了水文模拟的精度，也不利于真实掌握研究流域的降雨-径流响应。

目前，流域数字化过程中对于流域出口点的识别方法主要是依据水文站点的经纬度坐标，但是由于水文站点坐标定位时本身存在误差，导致了站点所对应位置往往并不是流域的真实出口点，现有方法则需要多次人为干预对站点位置进行修正，才能确定出流域出口，生成流域集水区域。

但是，人为干预修正流域出口点非常繁琐复杂，尤其当研究流域较多时，数字化过程也将更加的费时费力，不利于数字水文学的发展。

发明内容

发明目的：提供一种自动识别流域数字化过程中出口点的方法，以解决现有技术需要人为修正出口点，效率低下的问题。

技术方案：一种自动识别流域数字化过程中出口点的方法，包括如下步骤：

S1、读取研究区域的数字高程模型的原始数据并进行预处理，依次获得流向矩阵、累积汇水面积矩阵、栅格控制面积矩阵与河道栅格矩阵；

S2、基于水文站点的经纬度坐标，对水文站点在河道栅格矩阵内的位置进行定位，并判断水文站点是否位于河道栅格上，若不位于河道栅格上，则执行S3；若位于河道栅格上，该河道栅格即为流域潜在出口点，则执行S4；

S3、根据流向矩阵和水文站点，识别流域潜在出口点，并进行标记；

S4、根据栅格控制面积矩阵，提取流域潜在出口点的控制面积，并计算该面积与水文站集水面积之间的相对误差δ；

S5、根据设定的面积许可误差，判断δ的绝对值是否小于该面积许可误差，若判断结果为否，则执行S6；若判断结果为是，则执行S7；

S6、根据流向矩阵，对流域的潜在出口点进行修正，然后重复S4和S5；

S7、自动识别过程结束，当前的流域潜在出口点即为流域出口点。

所述数字高程模型的原始数据采用ArcInfo ASCII格式，地理坐标系统采用的是WGS84经纬度坐标。

栅格控制面积矩阵的提取步骤包括：

1）由数字高程模型分辨率计算得到单个栅格单元的面积；

2）以累积汇水面积矩阵为基础，将每个栅格单元的上游累积栅格数乘以单个栅格单元的面积，得到栅格控制面积矩阵。

河道栅格矩阵的提取步骤包括：

1）设定一个生成河道栅格的阈值；

2）利用栅格控制面积矩阵，按照给定的阈值对每个栅格单元进行标记，若栅格单元的控制面积小于给定的阈值，将该单元标记为无数据区；若栅格单元的控制面积大于给定的阈值，将该单元标记为河道栅格，得到河道栅格矩阵。

对水文站点在河道栅格矩阵内的位置进行定位包括如下步骤：

1）由水文站点的纬度坐标计算该站点在河道栅格矩阵内所处的行数：

SNx＝nrows-INT[(SLat-yllcorner)×60×(60/(3600×cellsize))]；

式中，SNx为水文站点在河道栅格矩阵内所处的行数；SLat为水文站点的纬度坐标；nrows为河道栅格矩阵的总行数；yllcorner为河道栅格矩阵西南角栅格单元对应的纬度坐标；cellsize为DEM分辨率；INT[]为取整函数；

2）由水文站点的经度坐标计算该站点在河道栅格矩阵内所处的列数：

SNy＝INT[(SLon-xllcorner)×60×(60/(3600×cellsize))]+1；

式中，SNy为水文站点在河道栅格矩阵内所处的列数；SLon为水文站点的经度坐标；xllcorner为河道栅格矩阵西南角栅格单元对应的经度坐标；

3）根据计算的SNx与SNy，在河道栅格矩阵内搜索与该行列数对应的栅格单元，将其标记为水文站点，实现对站点位置的定位；若该单元同时被标记为河道栅格，表明水文站点位于河道上，该单元即为流域的潜在出口点；若该单元同时被标记为无数据区，表明水文站点不位于河道上。

对流域的潜在出口点进行识别包括如下步骤：

1）根据流向矩阵，从水文站点所在栅格单元开始，沿着水流方向往下游搜索，搜索到第一个河道栅格后将其标记为C₁；

2）以水文站点所在栅格单元为中心，建立3×3的窗口，并搜索窗口内的河道栅格，当窗口内只有一个河道栅格时，则将其标记为C₂；当窗口内有多个河道栅格时，C₂则为距离水文站点最近的河道栅格；当窗口内没有河道栅格时，则扩大窗口重新搜索，直至C₂标记完成；

3）比较C₁与C₂的高程，并由流向矩阵判断较高的栅格单元是否流经较低栅格单元，如果判断结果为C₁大于等于C₂，则将C₁作为流域的潜在出口点；反之，则将C₂作为流域的潜在出口点。

S4中所述的相对误差δ的计算公式为：

$\mathrm{\δ}=\frac{A_i-A_g}{A_g}\×100\%;$

式中，A_i为流域潜在出口点的控制面积；A_g为水文站集水面积。

对流域的潜在出口点进行修正包括如下步骤：

1）当A_g＞A_i时，从修正前的潜在出口点开始，按照流向矩阵往下游搜索，搜索到的第一个河道栅格作为修正后的潜在出口点；

2）当A_g＜A_i时，从修正前的潜在出口点开始，以距离水文站点最近为原则，按照流向矩阵往上游搜索，搜索到的第一个河道栅格作为修正后的潜在出口点。

有益效果：与现有技术相比，本发明不仅利用了水文站点经纬度坐标对流域的潜在出口点进行自动识别，还引入了水文站集水面积对流域的潜在出口点进行自动修正，最终确定出数字化流域的出口点，实现了对流域出口点进行识别的自动化处理，既保证了出口点的定位精度，同时避免了人为干预，提高了流域数字化的效率，可以进一步促进数字水文学的深入发展。

附图说明

图1是本发明的示意图。

图2为本发明DEM原始数据头文件信息示意图。

图3为本发明预处理后的DEM示意图。

图4为本发明提取的流向示意图。

图5为本发明提取的栅格控制面积示意图。

图6为本发明提取的河道栅格示意图。

图7为本发明水文站点位置以及识别出的流域潜在出口点对应的集水区域示意图。

图8为本发明水文站点位置以及最终识别出的流域出口点对应的集水区域示意图。

具体实施方式

如图1至图8所示，本发明一种流域数字化过程中出口点自动识别的方法，包括以下步骤：

S1、读取研究区WGS84经纬度坐标系下ArcInfo ASCII格式的DEM原始数据，并对其进行预处理；再按照D8算法提取研究区流向矩阵；通过该水流方向数据提取研究区累积汇水面积矩阵；以此为基础，将每个栅格单元的上游累积栅格数乘以单个栅格单元的面积，提取栅格控制面积矩阵；最后按照给定的阈值对栅格控制面积矩阵进行重分类处理，筛选出河道栅格，提取河道栅格矩阵；

S2、基于水文站点的经纬度坐标，对水文站点在河道栅格矩阵内的位置进行定位，即利用水文站点的经纬度坐标计算该站点在河道栅格矩阵内所处的行列数：

SNx＝nrows-INT[(SLat-yllcorner)×60×(60/(3600×cellsize))]；

SNy＝INT[(SLon-xllcorner)×60×(60/(3600×cellsize))]+1；

式中，SNx为水文站点在河道栅格矩阵内所处的行数；SLat为水文站点的纬度坐标；nrows为河道栅格矩阵的总行数；yllcorner为河道栅格矩阵西南角栅格单元对应的纬度坐标；cellsize为DEM分辨率；INT[]为取整函数；SNy为水文站点在河道栅格矩阵内所处的列数；SLon为水文站点的经度坐标；xllcorner为河道栅格矩阵西南角栅格单元对应的经度坐标；根据计算的SNx与SNy，在河道栅格矩阵内搜索与该行列数对应的栅格单元，将其标记为水文站点，实现对站点位置的定位；若该单元同时被标记为无数据区，表明水文站点不位于河道上，在此情况下进入S3；

S21、若该单元同时被标记为河道栅格，表明水文站点位于河道上，该单元即为流域的潜在出口点，在此情况下进入S4；

S3、根据流向矩阵，对流域的潜在出口点进行识别，包括如下步骤；

1）根据流向矩阵，从水文站点所在栅格单元开始，沿着水流方向往下游搜索，搜索到第一个河道栅格后将其标记为C₁；

2）以水文站点所在栅格单元为中心，建立3×3的窗口，并搜索窗口内的河道栅格，当窗口内只有一个河道栅格时，则将其标记为C₂；当窗口内有多个河道栅格时，C₂则为距离水文站点最近的河道栅格；当窗口内没有河道栅格时，则扩大窗口重新搜索，直至C₂标记完成；

3）比较C₁与C₂的高程，并由流向矩阵判断较高的栅格单元是否流经较低栅格单元，如果判断结果为是，则将C₁作为流域的潜在出口点；如果判断结果为否，则将C₂作为流域的潜在出口点；

S4、根据栅格控制面积矩阵，提取流域潜在出口点的控制面积，并计算该面积与水文站集水面积之间的相对误差δ，即：

$\mathrm{\δ}=\frac{A_i-A_g}{A_g}\×100\%;$

式中，A_i为流域潜在出口点的控制面积；A_g为水文站的集水面积；

S5、根据设定的面积许可误差，判断δ的绝对值是否小于该面积许可误差，若判断结果为否，则进入步骤六；若判断结果为是，则进入步骤七；

S6、根据流向矩阵，对流域的潜在出口点进行修正，修正完成后重复S4和S5；其中，对流域的潜在出口点进行修正的具体方法如下：

1）当A_g＞A_i时，从修正前的潜在出口点开始，按照流向矩阵往下游搜索，搜索到的第一个河道栅格作为修正后的潜在出口点；

2）当A_g＜A_i时，从修正前的潜在出口点开始，以距离水文站点最近为原则，按照流向矩阵往上游搜索，搜索到的第一个河道栅格作为修正后的潜在出口点；

S7、自动识别过程结束，当前的流域潜在出口点即为流域出口点，流经该出口点的所有上游栅格组成该流域对应的集水区域。

以皖南山区临溪流域为例，研究区DEM原始数据采用美国太空总署（NASA）与国防部国家测绘局（NIMA）联合提供的SRTM（Shuttle Radar Topography Mission）数据；软件平台采用Visual Basic6.0；临溪水文站的站点经纬度坐标及其集水面积来源于中华人民共和国水利部水文局2012年12月刊印的《浙闽台河流水文资料》（中华人民共和国水文年鉴，2011，第7卷，第1册）。

步骤一：读取研究区DEM原始数据，并对其进行预处理，依次提取相应的流向矩阵、累积汇水面积矩阵、栅格控制面积矩阵与河道栅格矩阵，具体包括：

1）读取研究区WGS84经纬度坐标系下ArcInfo ASCII格式的DEM原始数据，实施例中该数据的头文件信息如图2所示；

2）对研究区DEM原始数据中的洼地进行填充，得到预处理后的DEM，如图3所示；

3）利用D8算法，即最陡坡度算法，提取出研究区的流向矩阵，如图4所示；

4）根据流向矩阵，对每个栅格单元的上游累积栅格数进行统计，得到研究区的累积汇水面积矩阵，再由DEM分辨率计算得到单个栅格单元的面积，将每个栅格单元的上游累积栅格数乘以单个栅格单元的面积，得到研究区的栅格控制面积矩阵，如图5所示；

5）设定实施例中生成河道栅格的阈值10km²，并对栅格控制面积矩阵进行重分类处理，即当栅格单元的控制面积小于10km²时，该单元则被赋值-9999，代表无数据区；而当栅格单元的控制面积大于10km²时，该单元赋值为1，代表河道栅格；由此即可提取出河道栅格矩阵，如图6所示。

步骤二：基于水文站点的经纬度坐标，对水文站点在河道栅格矩阵内的位置进行定位，即利用水文站点的经纬度坐标计算该站点在河道栅格矩阵内所处的行列数：

SNx＝nrows-INT[(SLat-yllcorner)×60×(60/(3600×cellsize))]；

SNy＝INT[(SLon-xllcorner)×60×(60/(3600×cellsize))]+1；

式中，SNx为水文站点在河道栅格矩阵内所处的行数；SLat为水文站点的纬度坐标；nrows为河道栅格矩阵的总行数；yllcorner为河道栅格矩阵西南角栅格单元对应的纬度坐标；cellsize为DEM分辨率；INT[]为取整函数；SNy为水文站点在河道栅格矩阵内所处的列数；SLon为水文站点的经度坐标；xllcorner为河道栅格矩阵西南角栅格单元对应的经度坐标；根据计算的SNx与SNy，在河道栅格矩阵内搜索与该行列数对应的栅格单元，将其标记为水文站点，实现对站点位置的定位；若该单元同时被标记为无数据区，表明水文站点不位于河道上，在此情况下进入步骤三；若该单元同时被标记为河道栅格，表明水文站点位于河道上，该单元即为流域的潜在出口点，在此情况下进入步骤四；本该实施例中，临溪水文站的经度坐标为118.533，纬度坐标为29.983，结合图2中相关信息的取值，可以计算出临溪水文站对应的SNx为381，SNy为400，其位置如图7或图8中箭头所示，由河道栅格矩阵可知，该对应单元并非河道栅格，在此情况下需要进入步骤三。

步骤三：对流域的潜在出口点进行识别：首先，根据流向矩阵，从水文站点所在栅格单元开始，沿着水流方向往下游搜索，搜索到第一个河道栅格后将其标记为C₁；其次，以水文站点所在栅格单元为中心，建立3×3的窗口，并搜索窗口内的河道栅格，当窗口内只有一个河道栅格时，则将其标记为C₂；当窗口内有多个河道栅格时，C₂则为距离水文站点最近的河道栅格；当窗口内没有河道栅格时，则扩大窗口重新搜索，直至C₂标记完成；最后，比较C₁与C₂的高程，并由流向矩阵判断较高的栅格单元是否流经较低栅格单元，如果判断结果为是，则将C₁作为流域的潜在出口点；如果判断结果为否，则将C₂作为流域的潜在出口点；按照上述步骤，即可识别出实施例中临溪流域的潜在出口点，位于河道栅格矩阵的第376行、第400列，流经该潜在出口点的所有上游栅格组成与该点对应的集水区域，如图7中灰色区域所示。

步骤四：根据栅格控制面积矩阵，提取流域潜在出口点的控制面积，并计算该面积与水文站集水面积之间的相对误差δ，即：

$\mathrm{\δ}=\frac{A_i-A_g}{A_g}\×100\%;$

式中，A_i为流域潜在出口点的控制面积；A_g为水文站的集水面积；实施例中A_i为229km²，临溪水文站的集水面积A_g为585km²，计算的δ为-60.9%。

步骤五：根据设定的面积许可误差，判断δ的绝对值是否小于该面积许可误差，若判断结果为否，则进入步骤六；若判断结果为是，则进入步骤七；实施例中设定的面积许可误差为2.0%，然而δ的绝对值为60.9%，明显大于面积许可误差，在此情况下需要进入步骤六。

步骤六：对流域的潜在出口点进行修正，即当A_g＞A_i时，从修正前的潜在出口点开始，按照流向矩阵往下游搜索，搜索到的第一个河道栅格作为修正后的潜在出口点；当A_g＜A_i时，从修正前的潜在出口点开始，以距离水文站点最近为原则，按照流向矩阵往上游搜索，搜索到的第一个河道栅格作为修正后的潜在出口点；修正完成后，重复步骤四和步骤五，直至步骤五中δ的绝对值小于设定的面积许可误差，进入步骤七。

步骤七：自动识别过程结束，当前的流域潜在出口点即为流域出口点，流经该出口点的所有上游栅格组成该流域对应的集水区域；实施例中，最终识别出的临溪流域出口点位于河道栅格矩阵的第376行、第399列，对应的栅格控制面积为587km²，临溪流域的集水区域如图8中灰色区域所示。

以上详细描述了本发明的优选实施方式，但是，本发明并不限于上述实施方式中的具体细节，在本发明的技术构思范围内，可以对本发明的技术方案进行多种等同变换，这些等同变换均属于本发明的保护范围。

另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明对各种可能的组合方式不再另行说明。

此外，本发明的各种不同的实施方式之间也可以进行任意组合，只要其不违背本发明的思想，其同样应当视为本发明所公开的内容。

Claims (7)

1.一种自动识别流域数字化过程中出口点的方法，其特征在于，包括如下步骤：

S1、读取研究区域的数字高程模型的原始数据并进行预处理，依次获得流向矩阵、累积汇水面积矩阵、栅格控制面积矩阵与河道栅格矩阵；

S2、基于水文站点的经纬度坐标，对水文站点在河道栅格矩阵内的位置进行定位，并判断水文站点是否位于河道栅格上，若不位于河道栅格上，则执行S3；若位于河道栅格上，该河道栅格即为流域潜在出口点，则执行S4；

S3、根据流向矩阵和水文站点，识别流域潜在出口点，并进行标记；

S4、根据栅格控制面积矩阵，提取流域潜在出口点的控制面积，并计算该面积与水文站集水面积之间的相对误差δ；

S5、根据设定的面积许可误差，判断δ的绝对值是否小于该面积许可误差，若判断结果为否，则执行S6；若判断结果为是，则执行S7；

S6、根据流向矩阵，对流域的潜在出口点进行修正，然后重复S4和S5；

S7、自动识别过程结束，当前的流域潜在出口点即为流域出口点；

其中，对水文站点在河道栅格矩阵内的位置进行定位包括如下步骤：

1)由水文站点的纬度坐标计算该站点在河道栅格矩阵内所处的行数：

SNx＝nrows-INT[(SLat-yllcorner)×60×(60/(3600×cellsize))]；

式中，SNx为水文站点在河道栅格矩阵内所处的行数；SLat为水文站点的纬度坐标；nrows为河道栅格矩阵的总行数；yllcorner为河道栅格矩阵西南角栅格单元对应的纬度坐标；cellsize为DEM分辨率；INT[]为取整函数；

2)由水文站点的经度坐标计算该站点在河道栅格矩阵内所处的列数：

SNy＝INT[(SLon-xllcorner)×60×(60/(3600×cellsize))]+1；

式中，SNy为水文站点在河道栅格矩阵内所处的列数；SLon为水文站点的经度坐标；xllcorner为河道栅格矩阵西南角栅格单元对应的经度坐标；

3)根据计算的SNx与SNy，在河道栅格矩阵内搜索与该行列数对应的栅格单元，将其标记为水文站点，实现对站点位置的定位；若该单元同时被标记为河道栅格，表明水文站点位于河道上，该单元即为流域的潜在出口点；若该单元同时被标记为无数据区，表明水文站点不位于河道上。

2.如权利要求1所述的自动识别流域数字化过程中出口点的方法，其特征在于，所述数字高程模型的原始数据采用ArcInfo ASCII格式，地理坐标系统采用的是WGS84经纬度坐标。

3.如权利要求1或2所述的自动识别流域数字化过程中出口点的方法，其特征在于，栅格控制面积矩阵的提取步骤包括：

1)由数字高程模型分辨率计算得到单个栅格单元的面积；

2)以累积汇水面积矩阵为基础，将每个栅格单元的上游累积栅格数乘以单个栅格单元的面积，得到栅格控制面积矩阵。

4.如权利要求1所述的自动识别流域数字化过程中出口点的方法，其特征在于，河道栅格矩阵的提取步骤包括：

1)设定一个生成河道栅格的阈值；

2)利用栅格控制面积矩阵，按照给定的阈值对每个栅格单元进行标记，若栅格单元的控制面积小于给定的阈值，将该单元标记为无数据区；若栅格单元的控制面积大于给定的阈值，将该单元标记为河道栅格，得到河道栅格矩阵。

5.如权利要求1所述的自动识别流域数字化过程中出口点的方法，其特征在于，对流域的潜在出口点进行识别包括如下步骤：

1)根据流向矩阵，从水文站点所在栅格单元开始，沿着水流方向往下游搜索，搜索到第一个河道栅格后将其标记为C₁；

2)以水文站点所在栅格单元为中心，建立3×3的窗口，并搜索窗口内的河道栅格，当窗口内只有一个河道栅格时，则将其标记为C₂；当窗口内有多个河道栅格时，C₂则为距离水文站点最近的河道栅格；当窗口内没有河道栅格时，则扩大窗口重新搜索，直至C₂标记完成；

3)比较C₁与C₂的高程，并由流向矩阵判断较高的栅格单元是否流经较低栅格单元，如果判断结果为C₁大于等于C₂，则将C₁作为流域的潜在出口点；反之，则将C₂作为流域的潜在出口点。

6.如权利要求1或2所述的自动识别流域数字化过程中出口点的方法，其特征在于，S4中所述的相对误差δ的计算公式为：

$\mathrm{\δ}=\frac{A_i-A_g}{A_g}\×100\%;$

式中，A_i为流域潜在出口点的控制面积；A_g为水文站集水面积。

7.如权利要求1所述的自动识别流域数字化过程中出口点的方法，其特征在于，对流域的潜在出口点进行修正包括如下步骤：

1)当A_g>A_i时，从修正前的潜在出口点开始，按照流向矩阵往下游搜索，搜索到的第一个河道栅格作为修正后的潜在出口点；

2)当A_g<A_i时，从修正前的潜在出口点开始，以距离水文站点最近为原则，按照流向矩阵往上游搜索，搜索到的第一个河道栅格作为修正后的潜在出口点。