中小流域时空变源混合产流模型及参数区域化方法研究

刘昌军，周 剑，文 磊，3，马 强，郭 良，丁留谦，孙东亚

（1.中国水利水电科学研究院 减灾中心，北京 100038；2.中国科学院 西北生态环境资源研究院，甘肃 兰州 730000；3.河海大学 水文水资源学院，江苏 南京 210098）

1 研究背景

小流域产汇流理论是水文模拟的重要基础，其中产流计算是关键环节，产流量的大小和过程会对汇流计算造成直接影响。早期的产流计算主要采用降雨径流相关图，1960年代以后，产流理论研究开始迅速发展，赵人俊等［1］通过大量实验研究得出了蓄满产流和超渗产流机制，指出蓄满产流在较为湿润的地方应用较好，超渗产流在较为干旱或半干旱的地方应用较好。随着产流理论的不断应用与深入研究，人们发现对于一个实际流域而言并不存在绝对的超渗或蓄满的单一产流模式，往往是多种产流模式并存。1985年，于维忠［2］针对界面产流规律提出了多种产流模式，认为地表面某一点的产流机制是随实际情况实时变化的；1996年，芮孝芳［3］通过研究各种径流是如何形成的，发现径流成分是在不同透水性的界面上产生的，用界面产流规律来统一不同径流成分的产流机制。

中小流域受流域几何特征及下垫面条件影响显著，产汇流机制复杂多变，多以超渗/蓄满混合产流模式为主［4］。相比传统单一产流模式模拟，目前的流域混合产流模拟方法主要有两种：一种是面积比例法，通过把流域按面积比例划分为超渗和蓄满两种产流区，利用相应产流公式分别计算不同区域产流量然后相加求得流域总产流量。这种方法计算简单、概念直观，但由于该方法多以常数面积比进行计算，无法准确模拟中小流域产流过程中蓄满和超渗产流区面积的时空变化。另一种是垂向混合法，通过将蓄满和超渗两种产流模式进行垂向组合计算产流量［5］。该方法利用通过空间分布的下渗曲线将流域内的落地净雨量分解为地面产流量和下渗水量。地面产流量直接产流，下渗水量在补给土壤缺水量之后再产流。相比面积比例法，该方法中流域蓄满和超渗产流区的面积比例是随前期土壤含水量和实际下渗量实时变化的，在对流域混合产流机理的模拟上更加贴近实际情况［6］。然而由于包气带中土壤水分下渗是受多种因素影响制约的复杂物理过程，极难对其进行精准详细的数值模拟，因此目前的垂向混合产流理论仍存在不足，尤其是对流域蓄满和超渗产流模式的时空转化机制仍缺乏清晰和定量地描述。

随着水文模型和大数据方法研究的不断进步，基于一定有资料流域的数据，通过机器学习等人工智能技术，对缺资料流域进行模型参数区域化已成为当今水信息学及防洪减灾技术研究中的热点之一。Oudin等［7］分别使用集总式水文模型（GR4J）和半分布式水文模型（TOPMODEL）对法国913个流域进行参数移植分析，发现基于属性相似法所得参数区域移植模拟结果明显优于回归分析法所得结果；在进行流域水文模型参数区域化时所选模型对流域产流过程物理要素描述精确性及所选参证流域在研究区内的分布数量也是决定参数移植成功率的重要条件。Garambois等［8］在利用分布式水文模型（MARINE）对法国地中海区域16个流域进行山洪模拟参数区域化研究后，发现参证流域数据质量、详细程度以及目标流域水文规律变化的波动性是影响模型参数移植效果的重要因素。Ragettli等［9］利用基于PRMS模型的机器学习CART方法对我国35个流域进行参数区域化研究后发现，CART方法比传统基于物理-气候相似法进行的参数区域化模拟效果提升20%。

本文选择我国不同水文分区的15个小流域开展模型验证，通过与国内外8个知名水文模型进行对比，旨在验证时空变源混合产流模型在我国山丘区小流域洪水模拟中的适用性。同时采用机器学习CART方法对河南省19个小流域进行参数区域化研究，通过缺一法检验CART方法对无资料小流域模型参数移植的有效性，为该方法向全国其他无资料小流域推广提供一定的理论基础。

2 时空变源混合产流模型

在山丘区小流域山洪产流过程模拟中，表层土壤由非饱和到饱和的动态变化过程是模拟的难点。传统的水文模型多采用基于经验的产流系数或Green-Ampt方程来进行模拟，由于这些方法忽略了土壤基质吸力对下渗的影响，对流域产流机制的物理变化过程描述简化，模拟精度不高。时空变源混合产流模型基于山洪灾害调查评价成果中小流域划分数据及基础属性数据划分地貌水文响应单元，根据响应单元的产流机制确定不同区域的产流模型算法，建立小流域平面混合产流模型；对在空间上识别为混合产流机理的地貌水文响应单元构建垂向混合产流模型，垂向混合产流模型采用组合概念水库方法来模拟包气带与饱和带土壤水量交换（分水源）。该模型从产流模式的平面混合、垂向混合和时段混合三个方面进行小流域时空变源混合产流模拟。

2.1 平面混合产流模型结构 平面混合产流模型构建的基本方法是利用小流域及下垫面数据进行地貌水文响应单元划分。根据响应单元的产流机制，确定不同的产流模型算法，建立小流域平面混合产流模型。通过查找文献和数值模拟，定性和定量确定各种类型土壤的下渗能力和对产流的影响，初步建立了小流域地貌水文响应单元的划分标准和框架。为快速识别小流域主导产流机制，将小流域地貌水文响应单元初分为快速响应单元、慢速响应单元、滞后响应单元和贡献较小单元，对于不同类型的响应单元匹配不同的主导产流机制。为了给小流域不同产流模式单元匹配相对适合的产流模型，以便能够较为准确地描述不同下垫面产流过程，合理反映产流机制，提高产流模拟精度，采用非饱和数值模拟方法，计算完成了9种典型土壤最大影响深度［10］、非饱和下渗过程和下渗速率等，根据这些计算结果将快速响应单元和滞后响应单元进一步细分为快、中、慢三种类型，确定的小流域水文响应单元划分框架及其主要对应产流模式见表1，后续可根据流域内不同区域的气象、水文、地形特征、土壤类型、土地利用和植被类型等资料进行调整修改。

表1 小流域水文响应单元划分框架及其对应的主要产流模式

2.2 垂向混合产流模型结构 垂向混合产流模型采用不同的概念水库来模拟土壤包气带与饱和带之间的水量交换［11］，从上至下依次是毛细水库、重力水库和变动面积饱和产流水库以及地下水库。为了更好模拟降雨入渗过程中的表层土壤水动态变化过程，将毛细水库所模拟的土壤区又细分为浅层土壤区和深层土壤区，在浅层土壤区根据瞬时下渗率、下渗能力、土壤含水量和储水能力判定超渗/蓄满过程，动态计算超渗产流量和因浅层土壤暂态饱和形成的蓄满产流量，深层土壤区则以蓄满产流机制计算。时空变源混合产流模型采用模块化的建模方式构建［12］，模型根据每个计算时段内模拟得到的水文响应单元土壤含水量和累计下渗量，计算超渗与蓄满产流的面积变化，同时通过对比瞬时下渗率与下垫面入渗能力，实现超渗/蓄满产流在每个地貌水文响应单元的时空转化。时空变源混合产流模型垂向结构见图1，参数及参数取值见表2。

表2 时空变源混合产流模型参数

图1 时空变源混合产流模型垂向混合产流模型结构

垂向混合产流模型中表层土壤下渗能力和土壤含水量变化兼具线性和非线性［13］。在以蓄满产流模式为主区域，包气带土壤下渗为稳定下渗；而在混合产流模式为主区域，时变的包气带土壤下渗率f(t)采用GARTO模型计算。在小流域实际产流过程中往往同时存在蓄满和超渗两种模式，因此，在时空变源混合产流模型模拟中，针对判定为混合产流为主的响应单元定义Rh为超渗部分总产流量，Rd为蓄满部分总产流量。蓄满部分总产流量由优先流层蓄满产流量Rd1和下层土壤蓄满后产流量Rd2组成：

式中Kgs为地下水库的渗漏系数。

2.3 非线性土壤下渗计算方法 高效准确地计算包气带土壤下渗过程是模拟小流域混合产流的基础与关键。2015年，Lai等［14］综合T-O算法和Smith等［15］基于Green-Ampt方程提出的土壤水再分配模型（GAR模型）［16］，建立了通过离散化土壤含水量来计算土壤内湿润锋下移的GARTO数值模拟模型［17］。该模型不但解决了非线性Richard方程求解过程中数值计算量大、不易收敛等难题，而且克服了传统Green-Ampt方程在对透水性较好土壤下渗模拟及双峰降雨过程中土壤下渗模拟精度不高的局限性。

GARTO模型将随深度分布的土壤含水量曲线对含水率θ进行离散［18］，分成n个含水率为Δθ的区间，则第k个区间内的湿润锋下移速度可表示为：

式中：θk-1和θk为第k-1个离散区间和第k个离散区间的含水率；G( )θi，θn为湿润锋整体基质吸力；θn为离散化后最后一个湿润锋的含水率。

其中土壤水的再分配基于GAR模型进行计算：

式中：r为间歇期降雨强度（r＜Ks）； β为形状系数； p为经验参数，该模型假设达西流作用于土壤表面至土壤深度1/p处之下。

将通过Brooks-Corey模型［19］求得的非饱和导水率曲线 K(θ)以及通过Van-Genuchten模型［20］计算得到的土壤水分特征曲线Ψ(θ)带入基于Philip方程［21］定义的湿润锋整体基质吸力计算公式，则可求解表层土壤的时变下渗能力：

式中：S2为土壤吸附性；θr为土壤残余含水率。

2.4 模型模块及相关参数 时空变源混合产流模型基于水文过程模块化搭建，是基于物理机制且支持日模式和暴雨洪水模式联合计算的水文模型［22］。模型的输入包括降水量，最小和最大空气温度，短波太阳辐射时间序列值等。输入数据的时间序列可以是日步长或小时步长，也可以是两种数据的混合格式。模型会自动判断数据的时间间隔。如果输入数据是日间隔，启动日模型相应的计算模块，如果输入的数据是小时间隔则启动洪水模型相应的计算模块。日模式的主要模块包括降水插值、截流及蒸发计算、蓄满产流计算、土壤水分计算、壤中流计算、地下水库计算和河道汇流计算等；洪水模型的主要计算模块包括降水插值、蒸发计算、时空变源混合产流计算、坡面径流运动波计算以及河道汇流运动波计算等。

2.5 模型对比验证 为了对比验证时空变源混合产流模型层（SKBY）在我国山丘区小流域洪水模拟中的应用效果，收集了我国不同水文分区（半干旱、半湿润、湿润）的15个中小流域暴雨洪水实测资料（表3），开展了时空变源混合产流模型（SKBY）与国内外知名的8个水文模型的模拟结果对比验证。所选8个模型包括：API模型（API）、大伙房模型（DHF）、初损后损法（HEC1）、SCS模型（HEC2），TOPMODEL模型（TOPMODEL）、新安江模型（XAJ）、中国山洪模型（CNFF）、PRMS模型（PRMS）。计算得到的不同模型场次洪水模拟结果如图2所示，模拟结果洪峰误差百分比和纳什系数见图3。

图2 多模型模拟结果对比

图3 9个模型计算15个中小流域的模拟结果

表3 所选小流域基本信息

应用时空变源混合产流模型对15个小流域共计202场洪水进行模拟，计算得到的15个流域场次平均纳什系数在0.34～0.91，均值为0.78，明显高于其它8个模型模拟平均纳什系数值（0.32～0.60），除甘肃马街流域纳什系数较低外（0.34），本模型对其它14个小流域洪水场次模拟纳什系数均高于0.65。而就甘肃马街而言，本模型模拟精度仍高于其他模型（均低于0.3）。进一步分析时空变源混合产流模型计算得到的平均洪峰误差，本模型平均洪峰误差在0.3%～16%，均值仅为7%，精度远高于其他8个模型相对洪峰误差均值（16%～148%）；在甘肃马街流域，本模型模拟相对洪峰误差仅为1%。从上述结果可看出，时空变源混合产流模型在中小流域暴雨洪水模拟中明显优于其他水文模型。

以流域洪峰相对误差和纳什系数组合作为评价指标将模拟结果分为三类：模拟结果较好（纳什系数大于0.5，洪峰相对误差小于20%），模拟结果可接受（纳什系数大于0.5，洪峰相对误差大于20%）和模拟结果较差（纳什系数小于0.5）。如图4所示，本模型在全部所选小流域场次洪水模拟中模拟较好场次约占场次总数的60%，相比其他模型模拟提高约20%。在湿润区流域场次洪水模拟中无较差模拟场次，模拟较好场次约占总场次数的80%。在半湿润半干旱区模拟较好场次约占总场次50%，仍明显高于其他模型。时空变源混合产流模型在我国不同水文分区中小流域洪水模拟中具有明显优势，尤其对半湿润半干旱中小流域场次洪水的模拟精度明显优于其他水文模型。

图4 不同模型在不同水文分区模拟结果占比

为了进一步验证时空变源混合产流模型模拟适用性和模拟产流机制的时空变化机理，以四川剑阁小流域为例，选取时空变源混合产流模型与蓄满产流模式的PRMS模型计算结果进行对比分析。四川剑阁小流域面积245.6 km2，流域内有雨量站和水文站各1个，本研究搜集整理了1982—2012年5场洪水及1998—2012年长系列降雨流量资料利用全国山洪灾害调查评价成果中小流域土地利用类型、土壤质地和小流域划分及属性等成果数据，基于地貌水文响应单元划分标准对剑阁小流域进行了地貌水文响应单元划分（图5），并建立了针对剑阁小流域的时空变源混合产流模拟模型。

图5 剑阁流域土地利用、土壤质地及地貌响应单元

利用长系列实测降雨径流资料验证本模型与PRMS模型的模拟结果。时空变源混合产流模型采用日模型和场次洪水模型同时计算模式。模型计算选取1998—2012年系列降雨径流资料。PRMS模型与时空变源混合产流模型模拟结果如图6。PRMS模拟剑阁流域多年洪水过程纳什系数仅为0.54，最大洪峰误差为-27.07%；时空变源混合产流模型模拟纳什系数为0.73，最大洪峰误差6.42%；时空变源混合产流模型计算长系列洪水模拟效果较好。

图6 混合产流模型与蓄满产流模型多场次洪水过程模拟结果

对比本模型与PRMS模型对场次洪水模拟结果。选取19820726、20120720等2场洪水模拟结果进行分析，计算结果见表4。从表4可以看出，两个模型均能较好地模拟剑阁流域的2场场次洪水，本模型计算得到2场场次洪水模拟结果较好，纳什确定系数分别为0.88和0.83，洪峰相对误差均小于20%；PRMS模型模拟结果可接受，纳什确定系数分别为0.75和0.63，但洪峰相对误差均大于20%；与蓄满产流模式的PRMS模型相比，本模型模拟的洪峰和洪量更加接近实测值，模拟精度明显提高。

表4 时空变源混合产流模型（SKBY）与PRMS模型模拟结果对比

另外，为了研究剑阁流域产流机制时空变化过程，以剑阁流域20120720场次洪水计算结果为例，分析两种模型对流域不同产流组分的模拟效果（见图7）。结果表明，PRMS模拟的洪峰、洪量主要由壤中流和地表径流组成，未能准确反映洪水的产流组分。这主要是因为该模型产流计算模式单一，对小流域产流过程中的物理机制变化描述不清，仅是通过参数优化方案从数学上寻找最优解，因此该模型对小流域产流过程中复杂物理机制的变化模拟效果不佳。时空变源混合产流模型模拟的产流组分既包括由于降雨强度超过表层土壤下渗能力而产生的超渗产流，也包括由浅层土壤暂态饱和后地表产生的蓄满产流，所模拟的物理机制变化更加清晰，因此能够准确反映山丘区中小流域暴雨洪水超渗蓄满时空混合产流过程。

图7 时空变源混合产流模型与PRMS模型计算产流组分

图8（a）为本模型模拟剑阁流域20120720场次洪水中不同地貌响应单元对洪水贡献。快速响应（慢）单元产流量对流域总产流量贡献最大，快速响应（快）单元产流量对流域总产流量贡献最小，模拟结果基本符合流域产流的物理机制。图8（b）、图8（c）和图8（d）为本模型模拟的不同响应单元产流过程中的产流组分分布。同一时段不同地貌响应单元的产流组分所占比例不同。地表产流所占比例为快速响应（快）＞快速响应（中）＞快速响应（慢）。不同地貌响应单元超渗与蓄满产流所占比例也不尽相同。就本场洪水而言，超渗产流占主导地位。同一种地貌响应单元产流机制随时间变化。在降雨开始阶段以下渗为主，地表不产流，当雨强超过土壤下渗能力后开始出现超渗产流。随着土壤含水量的增加，部分地区产流机制由超渗转为蓄满，但不同地貌响应单元产流机制的转换时间存在差异。模型根据下垫面物理产流过程和流域土壤含水量分布，通过判别瞬时下渗率与下垫面入渗能力的关系和土壤含水量与储水能力关系，计算超渗与蓄满产流的产流量变化，实现超渗/蓄满产流在每个水文响应单元时空转化，从物理机理方面刻画超渗和蓄满的时空转化过程。

图8 时空变源混合产流模型中不同响应单元产流组分分析

3 参数区域化及其验证

3.1 参数区域化方法 机器学习CART模型由Breiman等［23］在1984年提出，是目前应用最广泛的决策树模型之一。CART模型可以对复杂样本进行处理和分析，但是在实际应用中，其仍存在不能准确识别主要显性控制因素及过度拟合的问题。为了最大程度降低过度拟合对最终结果的影响，在实际应用中应该对其复杂程度加以限制，即只针对样本中的重要信息进行决策分析。本文采用主成分分析法（PCA，Principal Component Analysis）对样本流域属性信息进行提取，提取出相互关联度较低的独立主成分代替原始变量输入机器学习样本库，同时利用时空变源混合产流模型对流域进行模拟，对率定后的流域参数进行流域之间的随机参数移植，计算移植模拟的纳什系数作为预测标签输入机器学习样本库。把构建的样本库数据输入CART模型，将树节点限制在3～9个，最大限度地提高交叉验证测量的精度，以建立流域主成分与预测标签之间的关系，确定相似流域判别标准。最后使用“缺一法”逐一对样本流域进行验证，具体流程见图9。

图9 基于机器学习的小流域参数区域化方法与流程

3.2 参数区域化应用 选择我国河南省水文气象资料比较齐全的19个代表性中小流域开展参数区域化方法研究。所选19个中小流域分属黄河（6个）、淮河（6个）和长江（7个）三大流域，具体流域分布及基本信息见图10和表5。

图10 19个小流域分布及降雨分区

表5 河南19个小流域信息

时空变源混合产流模型在河南省19个小流域多场次洪水模拟率定平均纳什系数高于0.7以上的流域共计15个（图11），且平均洪峰误差均小于16%。相比其它15个流域流域，时空变源混合产流模型在东湾、口子河、米坪、石寺4个流域模拟精度稍低（纳什系数0.66～0.68，洪峰误差5%～14%）。利用时空变源混合产流模型对19个小流域进行相互参数移植模拟，共进行19×19=361次模拟，分别计算纳什系数。

图11 河南19个流域率定、随机移植和机器学习纳什系数值对比

将19个流域的属性特征按对其应物理意义进行分类，提取6类共计33个流域属性进行主成分分析。首先针对不同特征属性单位不同统一的问题，对不同单位特征进行标准化处理，之后通过主成分分析，设定累计方差贡献率之和80%作为阈值来提取该类别的主成分。主成分分析输入因子及提取主成分见表6。

表6 河南19个流域主成分分析输入及输出因子

选择一个流域作为目标流域，另一个流域作为参证流域，对目标流域的提取主成分与参证流域提取主成分求差后，对应以该参证流域对该目标流域进行参数移植后模拟纳什系数，作为一组机器学习样本。在剔除对应纳什系数低于0.5的样本后，本次河南19流域机器学习参数区域化研究样本共计342个，建立机器学习样本库。将机器学习样本输入CART模型，所得结果如图12所示，可以看出土壤质地（ΔPC3_soil）被确定为最主要的划分属性，地形属性、土地利用类型和排水属性也均为重要的划分属性，是确定19个流域相似的重要指标。按照平均纳什系数预测值高于0.5选择分类路径，符合条件的路径有NSE=0.92、0.83、0.77、0.72、0.62，共5条路径。

图12 河南19个流域CART分类树结果

采用缺一法对机器学习CART方法在河南19流域依次进行验证，统计机器学习参数区域化模拟的平均纳什系数值并与随机参数移植模拟平均纳什系数值进行对比，结果如图11所示。对河南19个流域进行随机参数移植后模拟结果平均纳什系数为0.25～0.69之间，机器学习参数移植后19个流域模拟平均纳什系数为0.42～0.85，模拟精度明显上升。定义移植后模拟平均纳什系数值大于0.5为移植成功，则随机移植成功率为74%，机器学习成功率为95%。

图13线段上下端为对该流域进行随机参数移植后模拟纳什系数值上下限，灰色柱为机器学习模拟纳什系数范围。如图11所示经机器学习参数区域化后模拟效果明显提高，对比随机参数移植模拟，平均纳什系数由0.52提高至0.70，提高约34%。

图13 河南小流域参数区域化模型模拟结果

同时对比Ragettli等采用PRMS模型和CART决策树与回归树方法进行的全国35个流域参数区域化结果。由于PRMS模型是蓄满产流机制，该模型参数区域化在南方湿润地区效果较好，在半湿润半干旱地区较差，移植后成功率约34%（见文献［9］）。时空变源混合产流模型参数区域化后移植成功率约95%，与PRMS模型参数区域化结果相比，基于物理机制模型和机器学习方法的参数移植效果提升明显。

4 结论

本文提出了适用于我国山丘区中小流域暴雨洪水模拟的时空变源混合产流模型，并针对缺资料小流域洪水模拟难题，采用机器学习CART方法在河南省19个小流域进行了参数区域化研究，结果表明：

（1）提出了考虑超渗/蓄满产流机制时空组合的小流域时空变源混合产流模型，该模型物理机制清晰，能够反映超渗/蓄满时空变化过程，可以较为准确地模拟非线性产流过程和机理；该模型参数物理意义明确，不确定性较小。通过与国内外8个水文模型在我国不同水文分区的15个小流域进行对比验证，该模型优于其他8个模型，且小流域洪水模拟精度提高20%以上，适用于湿润、半湿润和半干旱地区。

（2）研究了基于机器学习CART方法和时空变源混合产流模型的参数区域化方法，并对比分析了随机移植和机器学习移植的效果。通过分析河南19个小流域参数区域化结果发现，本文提出的参数区域化方法移植成功率比随机移植提高约20%。采用超渗/蓄满机制时空转换的时空变源混合产流模型比单一蓄满机制模型参数移植成功率提高约60%，因模型物理机制清晰、参数不确定性小，为缺资料小流域参数区域化和洪水模拟提供了新的技术手段和思路。

（3）提出的时空变源混合产流模型和参数区域化方法在典型流域取得了较好应用效果，但针对我国复杂地貌类型区，在地貌水文响应单元划分标准、缺资料地区参数区域化、不同类型区产流机制时空变化规律和模型机理等方面还需结合观测试验、数值模拟、大数据分析和人工智能等手段进一步研究，进而完善相关理论模型和算法。