基于水文模型模拟河流径流对气候变化的响应（Ⅰ）

Posted 2023-04-06 晔烨枼葉

水文模型包括集总式水文模型和分布式水文模型

集总式水文模型是一种基于水文学原理和统计学方法的水文模型，通过对流域内各个子区域水文过程进行集成，预测流域径流量的变化。该模型通常包括降雨输入、蒸散发、土壤水分、地下水和径流等组成部分，并考虑了流域内的时间和空间变异性。常见的集总式水文模型有SAC（Soil and Water Assessment Tool）、SVMM（Storm Water Management Model）和HBV（Hydrologiska Byråns Vattenbalansavdeling）等。传统的集总式水文模型很难直接利用气候模型的输出结果，因此需要用到分布式水文模型。

分布式水文模型是是一种基于地理信息系统（GIS）和数值计算方法的水文模型，将流域划分为多个子流域，并考虑每个子流域内水文过程的时空变化。与集总式水文模型不同，分布式水文模型考虑了流域内地形、土地利用、降雨分布等因素对径流产生的影响，能够更准确地预测流域水循环过程及其变化趋势。常见的分布式水文模型包括TOPMODEL、VIC（Variable Infiltration Capacity）、SWAT（Soil and Water Assessment Tool）、MIKE SHE等。这些模型通过将流域划分为网格状或子流域状单元，在每个单元中模拟降雨入渗、蒸散发、地表径流、地下径流等水文过程，最终计算流域的总径流量。

SWAT水文模型

SWAT水文模型是一种广泛应用于流域水文模拟、土地利用变化和管理措施分析的分布式水文模型。SWAT模型基于物理过程和经验关系，通过在每个子流域内模拟水文循环、作物生长、土壤侵蚀和养分迁移等过程，预测流域的水资源和土地资源变化。

SWAT模型主要包括三个部分：输入数据处理、模型运行和结果输出。输入数据包括气象数据、DEM（数字高程模型）、土地利用数据、土地管理措施、和网数据、水库数据等，其中气象数据可以通过天气站点、遥感卫星等多种途径获取。模型运行部分主要是在每个子流域内计算水文过程、农业活动和土壤侵蚀等过程，并考虑它们之间的相互影响。结果输出包括各个时段内的降雨量、蒸散发量、径流量、土壤侵蚀量、养分输出等多个指标。

SWAT模型具有分布式、综合性、可操作性等优点，被广泛应用于全球范围内的流域和水文问题研究，特别是对于流域的水资源管理和环境保护等方面具有重要意义。

水文响应单元HRU（Hydrological Response Unit）是分布式水文模型中的基本单位之一，通常指一个具有相对均匀的土地利用、土壤类型和地形条件等特征的区域，它可以是一个小流域、一块土地或一个网格单元。在模型中，将流域划分为若干个水文响应单元，可有效地反映流域内不同区域的水文过程和特征。通过对每个水文响应单元进行降雨入渗、蒸散发、土壤水分、地下水和径流等过程的模拟，可以预测流域的总径流量以及各个部分的变化情况。

水文响应单元的选择需要考虑以下因素：

1.地理特征：包括土地利用、土壤类型、地形高程等因素，通常选择具有相似特征的区域作为一个水文响应单元。

2.模型精度：水文响应单元的大小和数量对模型精度有影响，应根据实际情况进行选择和确定。

3.模型性质：不同的水文模型对水文响应单元的要求不同，需要根据具体模型进行选择和设计。

水文响应单元是分布式水文模型的基础，其合理的选择和设计可以提高模型的精度和可靠性。

径流模拟方法选择“日降水数据/径流曲线方法”，以日为单位进行径流演算。

潜在蒸发量的模拟方法选择Penman-Monteith方法。

河道演算采用Muskin-gum方法。

参数敏感性分析

复杂的水文模型需要通过引入大量参数来描述。参数的获取通常存在空间变异性、量测误差等问题致使模拟值与观测值之间吻合程度不高。模型参数灵敏度分析的目的在于分析判断哪些输入参数值的变动性对输出结果的变异程度影响更重要，从而提高模型的可用性。参数敏感度分析采用ArcSWAT提供的参数敏感度分析模块进行分析，此模块采用的方法是LH-OAT灵敏度分析方法，该方法对每一抽样点（LH抽样法）进行OAT灵敏度分析，灵敏度最终值是各局部灵敏度之和的平均值，是综合LH抽样法和OAT灵敏度分析的一种新方法，同时兼备两种方法的优点。

对影响径流的27个参数（SOL_K-土壤饱和导水率、SOL_AWC-土壤田间有效持水量、ESCO-土壤蒸发补偿系数、Alpha_Bf-基流α系数、GWQMN-浅层地下水径流系数、CH_K²-河道水力传导率、CN²-SCS径流曲线系数、BLAI-最大潜在叶面积指数、CANMX-最大覆盖度）进行了敏感性分析，其中SCS径流曲线系数（CN²）、土壤蒸发补偿系数（ESCO）、最大覆盖度（CANMX）、基流α系数（Alpha_Bf）、土壤田间有效持水量（SOL_AWC）等参数对径流影响比较大。

率定

率定是指对水文模型进行参数估计的过程，目的是调整模型参数以使模拟结果与实际观测数据尽可能地吻合。在率定过程中，通过对比模型预测值和实测数据，不断调整模型参数，直至达到最优模拟效果。

率定通常需要考虑以下几个方面：

1.参数选择：根据模型类型和流域特征等因素确定要进行率定的参数，同时要注意避免过度调整或调整无关参数。

2.观测数据处理：包括数据的质量检查、插补和平滑等预处理方法，以确保模型输入数据的可靠性和一致性。

3.评价指标的选择：常用的评价指标包括误差平方和（SSE）、均方根误差（RMSE）、相关系数（R²）等。

4.算法选择：常用的参数估计算法包括试错法、梯度下降法、遗传算法等，需根据模型类型和数据特征等选择适当的算法。SWAT模型常用的参数估计算法有试错法、单因素敏感性分析、自适应随机扰动优化法等。

评价指标值中Nash-Sutcliffe模型效率系数ENS计算公式如下：

式（1）中，为观测值；为模拟值；为平均观测值；n为观测的次数。当=时，ENS=1；如果ENS为负值，说明平均观测值比模拟值更具有代表性。

当ENS>0.75时，可以认为模拟效果好；0.36≤ENS≤0.75时模拟效果为令人满意；ENS<0.36时，模拟效果不好。

SWAT对年径流量的模拟能够较为真实的反映径流量的变化情况。 

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

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

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

河流比降

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

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

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

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

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

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

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

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

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

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

河流比降与GIS的关系

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

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

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

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

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

下图是一个测试数据，黄色部分为河流线，点为河流的节点。包含两份DEM，一份是原始的DEM，一份是填洼过的DEM(Fill_DEM) 。河流是使用了前面章节介绍到的河流提取方法提取出来的线状河流。读者可以参考《ArcGIS水文分析实战教程（6）河流提取与河网分级》的做法去提取。
测试的主要内容是分两种情况打印河流各节点的高程值。一是获取原始DEM的高程值，而是获取填洼后的DEM的高程值。通过高程值进行对比选用不同数据引起的差异。


第一张图是原始DEM的沿河高程值，第二张图是填洼后的DEM的沿河高程值。可以看到使用原始DEM的时候，有可能出现下游高程比上游高程还要高的情况，也就是说这个上游地势低的点刚好就是一个汇（凹陷点）。而填洼后的DEM则会将这个点给填平，确保不会断流，从而使得下游高程最少不会高于上游高程。对比其他高程值，是没有变化的，变化的部分就只是出现【汇】的点。

如果原始DEM的精度足够高的话，这种情况一般很少碰到。但为了稳妥起见，笔者还是建议对原始DEM做填洼，生成无凹陷点的DEM。读者可以参考之前章节 《ArcGIS水文分析实战教程（4）地形预处理》

如果所使用的河流是从DEM中提取出来的，那么在矢量河流提取的时候使用【栅格河流矢量化】时是否需要使用【简化折线】参数，这个参数是默认选上的，如下图所示：

不选择简化折线，河流线会基本上就是按照流量线走向，河流折点密度会比较高；选择简化的话，河流线的折点会被抽希一部分。这个怎么取舍？如果在地形变化比较大的情况下，不建议简化；如果地势趋平，那就直接简化线。这种方式会直接影响到最终的比降数值。

下图黄线部分是简化节点后的河流线，红线是没有简化节点的河流线，其形状基本上是重合的，而且在关键的河流节点上，也是完全重合。

也就是说，简化与不简化，其参与计算的河段数量是不一样的。上图数据中简化后的子河段数为32，不简化子河段的数量是92 。如果人为计算的话，那么可定是后者的计算量比较大了。

ArcGIS中实现河流平均比降计算

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'......因为您目前选中了多条河流。只能选中一条河流进行计算!......')

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

河流平均比降工具使用

可以将上面的python代码制作成为脚本工具，方便在ArcGIS中进行可视化的调用。工具的制作可以参考笔者另一篇博文《ArcMap自定义脚本工具制作》，里面有详细的论述。

该工具用到的ArcGIS工具涉及到ArcInfo级别 和 Spatial 扩展模块，请先确保ArcGIS软件具备相应的授权。

制作好工具之后可以双击调用该工具了，【约翰斯通-克罗斯】法需要用到无凹陷点DEM和河流线型数据。

从上面代码可以看到，工具会判断选中的河流数量，如果数量大于1，就不作计算，确保计算的是某一河段。ArcGIS的交互性非常的优秀，只要选中某一条河流，就会将选中的部分作为参数输入。

对应工具的流程，这里第一步需要选中要分析的一条河流。如下图

接着双击调用工具，两个参数分别是河流图层和无凹陷点DEM。一定要注意，这个DEM是填洼过的DEM，不是原始DEM数据。如下图

点击确定就可以执行计算结果，结果会在最终的提示框中找到，如下图所示：

至此，复杂的河流平均比降就计算完毕了！

总结

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

代码和制作ArcGIS工具的方法已经在文中贴出来了。读者可以自行参考制作。当然想省事的话，笔者也已经上传制作好的工具，可以点击跳转到 下载页面 。为了加强交流，防止伸手，多动手操作，笔者已经将下载分值设置为最高分10分。 积分不够的读者请多动手制作 ^_

如果读者对水文分析感兴趣的话，可以持续关注CSDN的GIS制图乐园，以及微信公众号【GIS制图乐园】。BY 李远祥