基于河道断面坡度值的Mann-Kendall突变检验的河道宽度提取方法

摘要

本发明公开了一种基于河道断面坡度值的Mann‑Kendall突变检验的河道宽度提取方法，包括以下主要步骤：基于DEM数据计算得到坡度栅格；在坡度栅格的基础上构建坡度值的点阵；数字化河道中心线；沿河道中心线布设河道断面基线；结合坡度点矩阵，在河道断面基线上进行坡度值突变检验，判定河道栅格；结合河道栅格以及断面基线与水平方向的夹角计算断面宽度。本发明具有数据来源稳定可靠、计算效率高、结果客观合理等优点，有利于流域河道宽度的快速估算和提取。

说明书

技术领域

本发明属于水文技术领域，具体涉及一种基于河道断面坡度值的Mann-Kendall突变检验的河道宽度提取方法。

背景技术

随着遥感、地理信息以及数字流域等技术的发展，基于栅格数字高程模型(DEM，Digital Elevation Model)的分布式水文模型以其充分考虑降雨和下垫面条件空间变化的特点，现已成为流域水文模型的发展趋势，尤其是在地形复杂的山区性中小流域中，分布式水文模型比传统的不能考虑流域内地形变化的集总式模型更有优势。在构建分布式的水文模型时，通常将流域划分为若干个正交的网格，以每一个网格作为计算单元。在每一个计算单元中通过产流模型计算产生的径流，然后按照一定的汇流顺序逐单元进行汇流计算，在汇流计算中十分重要的一环就是河道演算。当采用水力学方法进行河道演算时需要知道每一个栅格单元的河道宽度，这也是分布式水文模型建模的重点和难点之一。

为了进一步促进流域水文模型的发展，需要更深入研究河道宽度的计算方法。

在水流的长期作用下，河床有可能达到与所在河段的水力条件相应的平衡状态。在这种平衡状态下水深、河宽、比降等各种河相物理量之间可能存在某种稳定的函数关系，称之为河相关系。国内外的许多学者对河相关系开展过广泛的研究，这些研究中河宽的提取方法大致可以分为两类：一类在常规的河道水文数据的基础上，将河道宽度与断面常年的流量、输沙量等水文要素之间建立起经验性质的统计关系，即统计拟合法，这类方法主要基于数据之间的统计分析，进而拟合出河道宽度与所选水文要素之间的函数关系。该方法较为简单，因而开展研究的时间也较早。近十年来，随着地理信息系统的完善以及遥感技术的发展，出现了基于遥感影像来提取河道宽度的新方法。这一类方法不再基于河道宽度与某个或多个因子之间的函数关系，而在遥感影像上对河道所占的栅格进行识别与提取，直接对河道的宽度进行测量。结合这一类的方法的特点，可以将其称之为影像测量法。

统计拟合法原理简单，易于操作，但是经验性较强，难以取得一个普遍的通用公式，因此局限性较大。与河宽建立函数关系的最为主要而普遍的因子为上游集水面积，一般认为该因子与河道宽度成正相关系，集水面积越大，则河道越宽，即河道从上游至下游会因为集水面积的不断增大而逐渐加宽，然而这与实际并不相符，实际情况下，下游的局部河段的河道宽度也会小于上游。影像测量法具有形象直接、不受地域限制的特点。但是影像测量法主要基于河道水面像素与河道两侧坡地像素的差异来提取河道宽度，由于卫星影像易受到云层的影响，为得到全图清晰的影像，需要将多幅不同时间点的卫星图片和航拍图片拼接，因此无法避免不同时期水动力条件的变化对河道宽度提取结果的影响。因而对河道宽度估算的不准确限制了河道宽度这个参数在分布式建模中的应用，不利于国内分布式水文模型的发展。

针对以上不足，如何提出一种综合两种方法的优点，符合客观实际，能够结合有限的资料简便可靠地估算河道断面宽度的新方法，正是发明人需要解决的问题。

发明内容

为了解决现有技术中存在的不足，本发明提供了一种基于河道断面坡度值的Mann-Kendall突变检验的河道宽度提取方法，本发明具有数据来源稳定可靠、计算效率高、结果客观合理等优点，有利于流域河道宽度的快速估算和提取，值得推广。

为解决上述问题，本发明具体采用以下技术方案：

一种基于河道断面坡度值的Mann-Kendall突变检验的河道宽度提取方法，其特征在于，包括以下步骤：

步骤1，基于DEM数据计算得到坡度栅格；

步骤2，在坡度栅格的基础上构建坡度值的点阵；

步骤3，数字化河道中心线；

步骤4，沿河道中心线布设河道断面基线；

步骤5，结合坡度值的点阵，在河道断面基线上对坡度值Mann-Kendall突变检验，判定河道断面点；

步骤6，结合河道断面点以及断面基线与水平方向的夹角计算断面宽度。

前述的一种基于河道断面坡度值的Mann-Kendall突变检验的河道宽度提取方法，其特征在于，所述步骤1中基于DEM数据计算得到坡度栅格，具体包括以下步骤：

步骤1.1，以栅格单元Cell为中心，通过周围栅格的高程值与该栅格单元的高程值的对比，找出与其相比最低的栅格单元Cell_D；

步骤1.2，计算Cell和Cell_D之间的高程差DH_max和距离Dis；

步骤1.3，基于DH_max和Dis计算栅格单元Cell的坡度S；

S＝DH_max/Dis

步骤1.4，逐栅格循环，计算得到每一个栅格单元中的坡度S。

前述的一种基于河道断面坡度值的Mann-Kendall突变检验的河道宽度提取方法，其特征在于，所述步骤2中在坡度栅格的基础上构建坡度值的点阵，具体包括以下步骤：

步骤2.1，以栅格单元的几何中心为该栅格单元对应的点的位置；

步骤2.2，以栅格单元的坡度值为该栅格单元对应的点的坡度值；

步骤2.3，按照步骤2.1和步骤2.2中的方法将每一个栅格单元转换成点，构建坡度值的点阵。

前述的一种基于河道断面坡度值的Mann-Kendall突变检验的河道宽度提取方法，其特征在于，所述步骤3中数字化河道中心线，具体包括以下步骤：

步骤3.1，在地理信息系统中加载Landsat卫星影像，并将Landsat卫星影像中的7号波段、4号波段以及1号波段分别与光谱中的红、绿、蓝相对应，突出显示河道水面的轮廓；

步骤3.2，在地理信息系统中扫描河道水面的影像，沿着河道中心进行数字化，生成河道中心线。

前述的一种基于河道断面坡度值的Mann-Kendall突变检验的河道宽度提取方法，其特征在于，所述步骤4中沿河道中心线布设河道断面基线，具体包括以下步骤：

步骤4.1，对Landsat卫星影像进行目视检译，大致判断河道的最宽处，并结合比例尺测量其宽度得到最大河道宽度MaxB；

步骤4.2，结合4.1测量的最大河道宽度MaxB计算得到河道断面基线的长度B：

B＝MaxB+10D

式中：B为河道断面基线的长度；MaxB为河道的最宽处的宽度；D为DEM数据的分辨率，即以米为单位的栅格单元的边长；

步骤4.3，沿步骤3中提取出的河道中心线，在需要提取河道宽度处，布设若干条垂直于河道中心线的河道断面基线，其长度均为步骤4.2中计算得到的B。

前述的一种基于河道断面坡度值的Mann-Kendall突变检验的河道宽度提取方法，其特征在于，所述步骤5中结合坡度值的点阵，在河道断面基线上对坡度值Mann-Kendall突变检验，判定河道断面点的步骤包括：

步骤5.1，在步骤2构建的坡度值点阵的基础上，提取出与河道断面基线垂向距离在以内的点；

步骤5.2，对于一条河道断面基线，以河道断面基线与河道中心线的交点为原点，从原点出发，按照空间距离由近及远的方式，将河道中心线两侧的点的坡度值依次排列成数据序列，分别构建出向左岸延伸的数据序列Point_L和向右岸延伸的数据序列Point_R，其中左岸为面向河道下游方向时的左侧河岸，右岸为面向河道下游方向时的右侧河岸；

步骤5.3，分别对数据序列Point_L和数据序列Point_R进行Mann-Kendall突变检验，确定数据序列Point_L中发生坡度值突变的点的位置Site_L和数据序列Point_R中发生坡度值突变的点的位置Site_R，

在Mann-Kendall突变检验中首先构造一个秩序列d_k：

其中：d_k为秩序列，n为数据序列中数据的数目，k为取值在2到n之间的数，r_i为一个计顺序的数，若序列中的第i个数比它之前的数都大，则r_i＝1，否则r_i＝0；

其中：j为第1个数到第i个数之间数据的编号，x_i为编号为i的数，x_j为编号为j的数；

其中：E(d_k)为秩序列的数学期望；Var(d_k)为秩序列的方差；

定义统计量UF_k：

再将数据序列逆序，构造逆序统计量UB_k′:

UB_k^′＝-UF_k

式中：k′＝n+1-k

由UF_k和UB_k′分别绘制曲线。按从左到右的方向，两条曲线首次出现交点，则第1个交点对应的数据点的位置就是突变开始的位置；

对于数据序列Point_L和数据序列Point_R分别进行以上操作，分别判断出发生突变的位置：Site_L和Site_R；

步骤5.4，提取位于Site_L和Site_R之间的点，这些点均位于河道断面上，即为河道断面点。

前述的一种基于河道断面坡度值的Mann-Kendall突变检验的河道宽度提取方法，其特征在于，所述步骤6中结合河道断面点以及断面基线与水平方向的夹角计算断面宽度的步骤包括：

步骤6.1，过河道断面基线与河流右岸的交点作一条南北方向的直线，称其为河道断面投影基线；

步骤6.2，将河道断面点投影到河道断面投影基线上，并统计投影后的点的数目μ；当多个河道断面点投影后在河道断面投影基线上的位置重合，则在计数时将重合的多个河道断面点记为一个；

步骤6.3，过河道断面基线与河流左岸的交点作一条东西方向的水平直线，度量该水平直线与河道断面基线之间的锐角α；

步骤6.4，结合μ和α计算河道断面的宽度

本发明的有益效果：本发明提供的一种基于河道断面坡度值的Mann-Kendall突变检验的河道宽度提取方法，以坡度为基础，沿着河道中心线布设垂向的河道断面基线，提取了断面基线上的坡度分布，通过突变性检验判断出河岸的位置，并经过投影变换构建了栅格单元的边长与河道断面之间的函数关系，进而提取了河道宽度，这样既保证了计算结果的精度与可靠性，同时解决了缺少观测资料的河道宽度的计算问题。且本方法主要应用流域数字高程模型，数据来源稳定可靠，方法中变量之间的函数关系明确，有利于流域河道宽度的自动化生成，通过数字流域技术以简化提取步骤，同时，保证了结果的客观合理性，有利于分布式水文模型的直接调用，可以进一步促进数字水文学以及分布式模型的深入发展。

附图说明

图1是本发明的计算流程示意图；

图2为本发明计算出的坡度栅格示意图；

图3为本发明中布设的河道中心线和河道断面基线示意图；

图4为本发明中一条示例河道断面基线上提取出的坡度值分布示意图；

图5为本发明中一条示例河道断面基线上河道中心线左侧坡度值突变位置判断示意图；

图6为本发明中一条示例河道断面基线上河道中心线右侧坡度值突变位置判断示意图；

图7为本发明中一条示例河道断面基线上计算河道断面宽度时的几何关系示意图。

具体施方式

下面结合附图和具体实施例对本发明作进一步描述。

如图1所示，本发明提供的一种基于河道断面坡度值的Mann-Kendall突变检验的河道宽度提取方法，包括以下步骤：

步骤1，基于DEM数据计算得到坡度栅格，如图2所示，具体包括以下步骤：

步骤1.1，以栅格单元Cell为中心，通过周围栅格的高程值与该栅格单元的高程值的对比，找出与其相比最低的栅格单元Cell_D；

步骤1.2，计算Cell和Cell_D之间的高程差DH_max和距离Dis；

步骤1.3，基于DH_max和Dis计算栅格单元Cell的坡度S；

S＝DH_max/Dis

步骤1.4，逐栅格循环，计算得到每一个栅格单元中的坡度S。

步骤2，在坡度栅格的基础上构建坡度值的点阵，具体包括以下步骤：

步骤2.1，以栅格单元的几何中心为该栅格单元对应的点的位置；

步骤2.2，以栅格单元的坡度值为该栅格单元对应的点的坡度值；

步骤2.3，按照步骤2.1和步骤2.2中的方法将每一个栅格单元转换成点，构建坡度值的点阵。

步骤3，数字化河道中心线，具体包括以下步骤：

步骤3.1，在地理信息系统中加载Landsat卫星影像，并将Landsat卫星影像中的7号波段、4号波段以及1号波段分别与光谱中的红、绿、蓝相对应，突出显示河道水面的轮廓；

步骤3.2，在地理信息系统中扫描河道水面的影像，沿着河道中心进行数字化，生成河道中心线。

步骤4，沿河道中心线布设河道断面基线，具体包括以下步骤：

步骤4.1，对Landsat卫星影像进行目视检译，大致判断河道的最宽处，并结合比例尺测量其宽度得到最大河道宽度MaxB；

步骤4.2，结合步骤4.1测量的最大河道宽度MaxB计算得到河道断面基线的长度B：

B＝MaxB+10D

式中：B为河道断面基线的长度；MaxB为河道的最宽处的宽度；D为DEM数据的分辨率，即以米为单位的栅格单元的边长；

步骤4.3，沿步骤3中提取出的河道中心线，在需要提取河道宽度处，布设若干条垂直于河道中心线的河道断面基线，如图3所示，其长度均为步骤4.2中计算得到的B。

步骤5，结合坡度值的点阵，在河道断面基线上对坡度值Mann-Kendall突变检验，判定河道断面点，具体包括以下步骤：

步骤5.1，在步骤2构建的坡度值点阵的基础上，提取出与河道断面基线垂向距离在以内的点，如图4所示；

步骤5.2，对于一条河道断面基线，以河道断面基线与河道中心线的交点为原点，从原点出发，按照空间距离由近及远的方式，将河道中心线两侧的点的坡度值依次排列成数据序列，分别构建出向左岸延伸的数据序列Point_L和向右岸延伸的数据序列Point_R，其中左岸为面向河道下游方向时的左侧河岸，右岸为面向河道下游方向时的右侧河岸；

步骤5.3，分别对数据序列Point_L和数据序列Point_R进行Mann-Kendall突变检验，确定数据序列Point_L中发生坡度值突变的点的位置Site_L和数据序列Point_R中发生坡度值突变的点的位置Site_R，分别如图5、图6所示，

在Mann-Kendall突变检验中首先构造一个秩序列d_k：

其中：d_k为秩序列，n为数据序列中数据的数目，k为取值在2到n之间的数，r_i为一个计顺序的数，若序列中的第i个数比它之前的数都大，则r_i＝1，否则r_i＝0；

其中：j为第1个数到第i个数之间数据的编号，x_i为编号为i的数，x_j为编号为j的数；

其中：E(d_k)为秩序列的数学期望；Var(d_k)为秩序列的方差；

定义统计量UF_k：

再将数据序列逆序，构造逆序统计量UB_k′:

UB_k′＝-UF_k

式中：k′＝n+1-k

由UF_k和UB_k′分别绘制曲线。按从左到右的方向，两条曲线首次出现交点，则第1个交点对应的数据点的位置就是突变开始的位置。

对于数据序列Point_L和数据序列Point_R分别进行以上操作，分别判断出发生突变的位置：Site_L和Site_R；

步骤5.4，提取位于Site_L和Site_R之间的点，这些点均位于河道断面上，即为河道断面点。

步骤6，结合河道断面点以及断面基线与水平方向的夹角计算断面宽度，如图7所示，具体包括以下步骤：

步骤6.1，过河道断面基线与河流右岸的交点作一条南北方向的直线，称其为河道断面投影基线；

步骤6.2，将河道断面点投影到河道断面投影基线上，并统计投影后的点的数目μ；当多个河道断面点投影后在河道断面投影基线上的位置重合，则在计数时将重合的多个河道断面点记为一个；

步骤6.3，过河道断面基线与河流左岸的交点作一条东西方向的水平直线，度量该水平直线与河道断面基线之间的锐角α；

步骤6.4，结合μ和α计算河道断面的宽度

以宜昌市至宜都市的长江河段为例，本例的研究区中的数字高程数据(DEM)采用美国太空总署(NASA)与国防部国家测绘局(NIMA)联合提供的90m分辨率的STRIM(ShuttleRadar Topography Mission)数据。

步骤一、基于DEM数据计算得到坡度栅格，包括以下步骤：

1)以栅格单元Cell为中心，通过周围栅格的高程值与该栅格单元的高程值的对比，找出与其相比最低的栅格单元Cell_D；

2)计算Cell和Cell_D之间的高程差DH_max和距离Dis；

3)基于DH_max和Dis计算栅格单元Cell的坡度S；

S＝DH_max/Dis

4)逐栅格循环，计算得到每一个栅格单元中的坡度S。

步骤二、在坡度栅格的基础上构建坡度值的点阵，包括以下步骤：

1)以栅格单元的几何中心为该栅格单元对应的点的位置；

2)以栅格单元的坡度值为该栅格单元对应的点的坡度值；

3)按照(1)和(2)中的方法将每一个栅格单元转换成点，构建坡度值的点阵。

步骤三、数字化河道中心线，包括以下步骤：

1)在地理信息系统中加载Landsat卫星影像，并将Landsat卫星影像中的7号波段、4号波段以及1号波段分别与光谱中的红、绿、蓝相对应，突出显示河道水面的轮廓；

2)在地理信息系统中扫描河道水面的影像，沿着河道中心进行数字化，生成河道中心线。

步骤四、沿河道中心线布设河道断面基线，包括以下步骤：

1)对Landsat卫星影像进行目视检译，大致判断河道的最宽处，并结合比例尺测量其宽度为1523.81米；

2)结合(1)中测量得到的最大河道宽度计算得到河道断面基线的长度B，B＝1523.81+10×90＝2423.81米；

3)沿(1)中提取出的河道中心线，在需要提取河道宽度处，布设若干条垂直于河道中心线的河道断面基线，其长度均为(2)中计算得到的B；

步骤五、结合坡度值的点阵，在河道断面基线上对坡度值Mann-Kendall突变检验，判定河道断面点的步骤包括：

1)在步骤3构建的坡度值点阵的基础上，提取出与河道断面基线垂向距离在45米以内的点；

2)对于从坐标点(经度：111.509，纬度：30.289)到坐标点(经度：111.525，纬度：30.306)的河道断面基线，以河道断面基线与河道中心线的交点为原点，其坐标为(经度：111.516，纬度：30.297)，从原点出发，按照空间距离由近及远的方式，将河道中心线两侧的点的坡度值依次排列成数据序列，分别构建出向左岸延伸的数据序列Point_L【00035411010】和向右岸延伸的数据序列Point_R【000006600】，其中左岸为面向河道下游方向时的左侧河岸，右岸为面向河道下游方向时的右侧河岸；

3)分别对数据序列Point_L和数据序列Point_R进行Mann-Kendall突变检验，确定数据序列Point_L中发生坡度值突变的点的位置Site_L，和数据序列Point_R中发生坡度值突变的点的位置Site_R，

在Mann-Kendall突变检验中首先构造一个秩序列d_k：

其中：d_k秩序列，n为数据序列中数据的数目，k为取值在2到n之间的数，r_i为一个计顺序的数，若序列中的第i个数比它之前的数都大，则r_i＝1，否则r_i＝0；

其中：j为第1个数到第i个数之间数据的编号，x_i为编号为i的数，x_j为编号为j的数；

其中：E(d_k)为秩序列的数学期望；Var(d_k)为秩序列的方差；

定义统计量UF_k：

再将数据序列逆序，构造逆序统计量UB_k′:

UB_k′＝-UF_k

式中：k′＝n+1-k

由UF_k和UB_k′分别绘制曲线。按从左到右的方向，两条曲线首次出现交点，则第1个交点对应的数据点的位置就是突变开始的位置。对于数据序列Point_L和数据序列Point_R分别进行以上操作：Site_L中发生突变的位置在第4个数字0到第5个数字6之间，而在Site_R中发生突变的位置在第5个数字0到第6个数字6之间。

4)提取位于Site_L和Site_R之间的点，这些点均位于河道断面上，即为河道断面点。

步骤六、结合河道断面点以及断面基线与水平方向的夹角计算断面宽度的步骤包括：

1)过河道断面基线与河流右岸的交点作一条南北方向的直线，称其为河道断面投影基线；

2)将河道断面点投影到河道断面投影基线上，并统计投影后的点的数目μ，其值为10。若多个河道断面点投影后在投影基线上的位置重合，则在计数时将重合的多个点记为一个；

3)过河道断面基线与河流左岸的交点作一条东西方向的水平直线，度量该水平直线与河道断面基线之间的锐角α，其值为40.59度；

4)结合μ和α计算河道断面的宽度

在Google Earth上经过测量得到该河道断面的宽度约为1353.64米，相对误差为2.1％，可以认为在该断面上本方法对于河宽的计算具有很好的精度。

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

以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。