ArcGIS水文分析实战教程（10）河流平均比降计算

本章导读：河流比降指的是河流任意两端点间的高程差与两点间的水平距离之比值，简单的概括为单位河长内的落差，其单位一般都是‰ 。 河流（纵）比降对地质灾害、水运、水电、矿砂和新构造运动研究有实用价值。比降是水文学上重要的一项指标。在绝大部分的GIS通用平台上，都没有实现平均比降的计算，原因在于计算比降的条件比较苛刻，不符合通用GIS的要求。本章同样是一个实操章节，笔者将会利用ArcGIS软件，结合其ArcPy编写一个平均比降计算的工具，在通用GIS软件上实现河流平均比降计算。 BY 李远祥

河流比降是很多水文研究者研究的内容，其水文意义非常的大。这里研究的河流比降指的是河流的纵比降，是河流任意两端点间的高程差与两点间的水平距离之比值，即单位河长的落差。

正是因为直接表达的是单位河长的落差，所以在水文学上其最直接的影响就是流速，带来更多的水动力引发的能量转移，如势能与动能的转化（涉及到水电站的建设）；河水对河床和河岸的冲刷，影响泥沙量；大面积持续降雨时，洪峰到来的时间和强度等。所以，比降在水文学上从来就不缺少话题。

河流的比降分为床面比降和水面比降。床面比降，用以表示河床纵断面地形的变化；水面比降即河流中任意两端点间的瞬时水面高程差与其相应水平距离之比，用以表明河流全程或分段的水面坡度，故又称水力坡度，通常说的河流比降就是河流水面比降。

在恒定均匀流情况下，水面坡度恰好等于河底坡度。所以，在没有大量断面数据支撑的情况下，根据地形和河流线型去研究的河流比降都是床面比降。

笔者在这里研究的都是床面比降，因为只考虑地形与河流本身的几何形状，没有任何的站点数据支撑。但在非洪水情况下，这个床面比降与水面比降基本上是吻合的。

在河流相距较短的两个点上计算河流比降，意义不是十分大。一般需要计算的是整个河流或者某一河段的平均比降，这对于空间跨度较大的河流来说才有研究价值。

本章研究的内容也就是河流（段）的平均比降。研究河流平均比降的方法有几种，比较流行的是【约翰斯通-克罗斯】法。其原理是根据比降大致均匀的原则，将河段划分为若干个子河段，然后去计算整体的平均比降。其公式如下：

公式中个参数如下
S—- 河段平均河底比降
Li—-第i个子河段长度
Si—-第i个子河段河底比降
N—- 子河段数量

另一种方法是【劳伦斯】法，其需要经验系数和流速，相对来说比较复杂，需要结合测站的一些数据。由于笔者没有这些数据，就不作论述。

本章考虑的是在ArcGIS中实现【约翰斯通-克罗斯】法的计算。

利用GIS手段来研究水文，其核心方式是使用DEM数据去模拟水文过程。GIS的空间叠加能力非常的管用，非常简单的就能获取到河流线对应的DEM像元值，也就是河流对应的高程值，这样用于计算比降是非常有用的。最起码解决了河流高程点获取的问题。

河流沿线高程的获取，不能直接使用原始的DEM数据。原因有两点：

1. 如果河流数据与DEM数据不是同一份数据，那么像元点与河流线的重合程度有待商榷。

2. 如果河流数据与DEM数据吻合（或者河流数据直接从DEM数据提取出来的），其高程的获取所依赖的DEM必须是填洼后的无凹陷点DEM。

河流线属于一维数据，DEM在平面上属于二维数据，其相交之后会产生一个像元点的连续像元。如果了解过D8单流向算法的话，就可以很容易理解，如果直接从原始DEM中获取像元值，其值有可能是一个汇点（凹陷点）。然而使用填洼后的DEM，可以确保该点与周边一点的高程持平。

ArcGIS 除了前面章节提及到的水文分析工具及其相关工具之外，还有非常庞大的工具集对空间计算进行支撑。由于没有断面数据的支撑，这里只考虑通过地形和河流线型进行计算，也就是说实现【约翰斯通-克罗斯】法。

根据【约翰斯通-克罗斯】法的计算公式，可以推导出其分析流程

1. 通过河流线提取河流折点，利用这些折点与DEM数据叠加，获取高程值

2. 利用河流折点，获取其坐标值，根据两点距离公式，计算每一段子河段的长度

3. 根据每个点的高程值，计算出相邻两点之间的高差

4. 根据高差与子河段长度，计算出各子河段的河底比降

5. 汇总这些数据计算出整个河段的平均河底比降

在这5各步骤中，只有第1步与ArcGIS工具有关。河流的折点提取，可以使用【数据管理】–【要素】–【要素折点转点】工具来实现，如下图所示：

至于利用折点与DEM数据叠加，获取其高程值，可以使用【数据管理】–【栅格】–【栅格属性】–【获取像元值】工具来获取。这个工具每次只能实现一个坐标点与DEM数据叠加查询。如下图

所以，为了方便计算，笔者采用ArcPy去写一个脚本工具，实现参数的灵活输入。其中heightList就是用来保存高程值的数组。

    # 获取河流的节点，以用作计算高差和河长
    for row in arcpy.da.SearchCursor(outLine1, ["SHAPE@XY"]):
        # 从点图层数据获取所有的点的XY
        x, y = row[0]
        #从地形数据中获取其高程的象元值
        pointStr = str(x) +" "+str(y)
        result = arcpy.GetCellValue_management(dem,pointStr)
        riverPoints.append(row[0])
        heightList.append(str(result))

第2-5步基本上就是一些数学的运算，利用python就很好实现了。为了让广大水文研究者用起来，笔者将这段代码公开，并且制作成为工具。在代码中笔者已经预留了ArcGIS的工具参数，方便制作成为脚本工具。

河流平均比降【约翰斯通-克罗斯】法的实现代码如下

#############################################################
#                                                           #
#   该工具利用【约翰斯通-克罗斯】方法计算河底平均比降           #
#   在使用的时候必须先选择一条河流，工具会判断输入的河流数量     #
#   大于一条河流时不予以计算                                  #
#   用于计算的DEM为无凹陷点DEM                                #
#   BY 李远祥                                                #
#                                                            #
##############################################################


# -*- coding: gbk -*-
import arcpy
import math
inputLine =arcpy.GetParameterAsText(0) #河流图层

dem = arcpy.GetParameterAsText(1) #无凹陷点DEM
riverLen = 7873.232038

arcpy.AddMessage(u'..................')
arcpy.AddMessage(u'......开始执行...........')

# 记录高程值的数据
heightList = []

iputCout=0
mmCursor = arcpy.da.SearchCursor(inputLine, ["SHAPE@LENGTH"])
for row1 in mmCursor:
    riverLen = row1[0]      #顺便获取河流长度
    iputCout=iputCout+1

arcpy.AddMessage(u'......获取了河流长度......')
#如果选中的是一个线要素，则执行计算，否者不计算多条河流
if iputCout==1:
    ##先将线转为折点
    outLine1 = arcpy.FeatureVerticesToPoints_management(inputLine, "in_memory/outline")
    arcpy.AddMessage(u'......河流线已经转换为点集.......')

    #河流节点
    riverPoints = []

    # 获取河流的节点，以用作计算高差和河长
    for row in arcpy.da.SearchCursor(outLine1, ["SHAPE@XY"]):
        # 从点图层数据获取所有的点的XY
        x, y = row[0]
        #从地形数据中获取其高程的象元值
        pointStr = str(x) +" "+str(y)
        result = arcpy.GetCellValue_management(dem,pointStr)
        riverPoints.append(row[0])
        heightList.append(str(result))

    arcpy.AddMessage(u'......已经获取了折点的高程值......')
    #子河段长度

    subRiverList = []
    hDiffList = []
    # 计算每一个子河段的长度
    for i in range(0,len(riverPoints)-1):
        x,y = riverPoints[i]
        x1,y1 = riverPoints[i+1]
        #计算河段长度
        dis =math.sqrt(pow((x1-x),2)+pow((y1-y),2))
        subRiverList.append(dis)

    arcpy.AddMessage(u'......河流子河段长度计算完成完成.......')

    # 计算从高程从低到高部分的差距
    listLen = len(heightList)
    for j in range(0,listLen-1):
        h1 = float(heightList[j])
        h2 = float(heightList[j+1])
        xheight = 0
        xheight = h1-h2
        hDiffList.append(xheight)

    arcpy.AddMessage(u'......各折点的高程差计算完成......')

    # 计算每个河段的比降
    subbijiang = []
    for i in range(len(hDiffList)):
        Si= hDiffList[i]/subRiverList[i]
        subbijiang.append(Si)

    arcpy.AddMessage(u'......子河段比降计算完成......')
    # 计算最终的比降
    # 计算分子
    sumUp = 0
    for i in range(len(hDiffList)-1):
        sumUp = sumUp + subRiverList[i]*pow(subbijiang[i],0.5)

    S = pow((sumUp/riverLen),2)

    arcpy.AddMessage(u'......该河流的总长为： ' + str(riverLen)+u'   平均比降为: ' +str(S))
    arcpy.AddMessage(u'......计算完毕.......')

else:
    arcpy.AddMessage(u'......您选中了 '+str(iputCout) +  u'  条河流......')
    arcpy.AddMessage(u'......工具已经停止计算,没有算出您想要的结果......')
    arcpy.AddMessage(u'......因为您目前选中了多条河流。只能选中一条河流进行计算!......')

从代码中可以看到，这个工具针对的是一条河流进行计算的。如果想指定河流的某一段进行计算，那么就需要将河流转为折点，选择对应的折点，生成新的河段才能参与计算。这个过程需要用到【点集转线】工具进行重构要分析的河段。

水文分析是一个复杂的过程，需要有深厚的水文学基础作为背景。而GIS是一个通用工具，提供非常庞大的分析功能。二者要结合起来，不但需要学习两个专业的知识，还需要具备知识整合的能力。