时空变源混合产流模型及其在小流域洪水模拟中的应用<br/>——以四川省、甘肃省的4个小流域为例

时空变源混合产流模型及其在小流域洪水模拟中的应用
——以四川省、甘肃省的4个小流域为例

2021-08-25张相芝

水资源与水工程学报 2021年3期

张相芝，周剑，文磊

(1.中国科学院西北生态环境资源研究院，甘肃兰州 730000； 2.中国科学院大学，北京 100049； 3.河海大学水文水资源学院，江苏南京 210098)

1 研究背景

山地和丘陵覆盖了中国69%的陆地面积，山丘区降雨集中，山洪灾害频发[1]。仅2018年夏季，全国洪水造成的经济损失就高达1 615.47×108元，四川、甘肃等6省市的山洪经济损失占全国山洪损失的70.2%[2-3]，山洪灾害已经成为制约当地经济发展的突出问题。山洪灾害预警是减少人员伤亡和财产损失的重要手段，在暴雨洪涝灾害预报中，流域水文模型和基于实测的降雨径流方法仍然是进行洪水实时预报和模拟的主要方式方法[4-5]。

我国山丘区的山洪预报预警系统尚处于发展阶段，相关学者在中小流域水文模型应用方面进行了诸多探索：利用多种气象水文模型耦合流域气候水文过程[6]；利用概念式水文模型如PRMS、MIKE等模拟雨洪过程[7-8]；利用水动力学模型FloodArea获取不同风险等级的临界雨量阈值[9]；以及利用数据驱动方法与传统水文模型结合来提高模拟精度[10]。国内山丘区小流域洪水预报的主要挑战是建立合适的模型来模拟短历时强降雨下的山洪过程。对于这一问题，刘昌军等[11]、王璐等[12]、黄艳等[13]进行了多种水文模型的示范应用和比较研究，Samadi等[14]在国外进行了多个流域的洪水模拟验证。上述研究表明，现有的大部分降雨径流模型通常需要根据特定流域和多年实测数据进行大幅度的校正，因此迫切需要提出适用于山区洪水过程的产流理论，并将之推广应用于暴雨洪水模拟中[15-16]。

由于强降雨条件下的混合产流计算较为复杂，目前没有针对于山丘区短历时条件下的洪水演算水文模型。已有相关研究基于垂向混合模型理论[17]和山洪评价中应对小流域高效划分的基本原则[18]，建立了基于山坡地貌水文响应单元的时空变源混合产流模型[19-20]，但对于时空变源混合产流模型的示范性研究很少，缺乏典型地貌区的混合产流特点研究。本研究选取4个小流域，利用时空变源混合产流模型模拟洪水过程，评估模型精度和下垫面信息对小流域产流成分的影响，为山丘区洪水灾害防治提供决策支持。

2 时空变源混合产流模型

时空变源混合产流模型(spatio-temporal variable source mixed model, SVSM)是在垂向混合产流模型的基础上，通过GIS和遥感技术识别小流域基础地貌信息，针对不同的地貌水文响应单元选择不同产流计算方法，利用离散湿润锋含水量代替Richard方程计算水在非饱和土壤中的运动，结合Brooks-Corey(BC)模型和Van-Genuchten(VG)模型计算土壤水分特征参数，构建地貌水文响应单元下的非饱和土壤下渗过程[19]。SVSM模型改进了暴雨条件下的河道汇流滞后算法，建立了一种新的垂向降水入渗和水分再分配的计算方法。

2.1 垂向混合产流模型结构

时空变源混合产流模型的架构采用垂向混合产流模型的基本结构[21]。垂向混合产流模型是Bao等[17]为半干旱流域产流模式提出的，以蓄水容量曲线为参考将超渗产流和蓄满产流结合起来[21-23]，以不同的概念水库概化土壤在不同程度蓄水量时的水文运动过程和产流形式。垂向混合产流模型结构如图1所示，时空变源混合产流模型参数和主要变量如表1所示。

表1 时空变源混合产流模型参数和主要变量

时空变源混合产流模型包括时空变源和混合产流两个方面，时空变源是指土壤含水量在降雨入渗、蒸发和重力及基质吸力的共同作用下的时空变化过程；混合产流是指在土壤含水量时空变化下的流域产流呈现超渗/蓄满时空动态组合的过程。模型根据每个计算时段内计算的水文响应单元含水量和累计下渗量计算超渗与蓄满产流的面积变化，并判别雨强与下垫面入渗能力的关系，实现超渗/蓄满产流在每个地貌水文响应单元时空转化。

模型采用垂向混合法进行产流计算，包括超渗产流机制和蓄满产流机制。超渗产流定义为降水超过土壤瞬时渗透速率时的产流；蓄满产流定义为降水超过土壤饱和稳定下渗速率时形成的产流。设定模型在蓄满产流的面积上土壤下渗能力为稳定下渗率；在超渗产流的面积上采用包气带非饱和下渗模型计算产流。根据下渗速率将地表分为透水和不透水两部分，透水部分又划分为优先流层和非优先流层，超渗产流在两种土壤表面有不同的计算方式。总超渗产流计算公式如下：

Rh=PE·perip+max(PE·(1.0-perip)-

f·Δt, 0)

(1)

式中：第1项为不透水面积下的超渗量，mm；第2项为透水面积的超渗产流量，mm；Rh为总产流量，mm；PE为净雨量，mm；perip为不透水面积比；f为瞬时土壤下渗率，m/s； Δt为计算时间步长，h。

上层土壤优先流层的出流由出流系数和优先流层含水量决定：

(2)

式中：Wp为优先流含水层的含水量, m3；kp1为优先流出流的线性系数；kp2为优先流出流非线性系数。

深层土壤和变动面积饱和产流水库以蓄满产流机制计算，当重力含水层出流使得优先流含水层饱和时生成蓄满产流，蓄满产流分为优先流层蓄满径流和深层蓄满后向上补给产生的蓄满径流。根据优先流区域的下渗量计算得到优先流层蓄满径流为：

Rd1=max(0,Wpt-1+perp·PE·

(1.0-perip)-Rh2-Wpmax)

(3)

当重力含水层出流流入优先流含水层使得优先流含水层饱和时，产生深层蓄满后向上补给产生的蓄满径流量Rd2为：

Rd2=max(0,Wp+Wgrvf-Wpmax)

(4)

式中：Wpt-1为优先流层上一时刻土壤含水量，m3；perp为优先流面积百分比，%；Wpmax为优先流层土壤最大含水量，m3；Wp为优先流层土壤含水量，m3；Wgrvf为深层土壤达到田间持水量后向优先流的排水量,m3。

下层土壤产生的壤中流由壤中流出流系数决定：

(5)

式中：Sl为下层土壤的含水量，m3；ki1为壤中流出流的线性系数；ki2为壤中流出流的非线性系数。

地下水水库的出流(基流)由出流系数和水库含水量决定：

Rg=Sg·kg

(6)

式中：kg为地下水水库出流系数;Sg为地下水水库含水量，m3。

2.2 非饱和土壤下渗

在垂向混合产流模型基础上，时空变源混合产流模型采用了一种新的垂向降水入渗和水分再分配的计算方法，即非饱和土壤下渗。将含水量范围离散为具有恒定含水量的区间，假定区间之间土壤含水量是连续的且有水力联系。湿润锋在区间内推进造成的水分缺失可以由毛细吸力为驱动力的侧向土壤水分再分布来补充。

计算主要思路是用单个含水量离散区间的湿润锋模拟水分渗透到土壤过程中整体湿润锋的区间流动过程。具体为在一个特定的含水量离散区间里，湿润锋首先在重力和毛细力以及表面边界条件的相互作用下(比如渗透或蒸发)垂直移动，然后计算地下水位的上升或下降，计算饱和含水量离散区间与地下水位接触部分的压力导致的毛细上升高度，最后在不同含水量土柱中通过与其他含水量离散区间吸力的相互作用计算水分再分配过程，使水在不同的含水量离散区间内横向转移，并导致湿润锋变化。

使用BC模型和VG模型计算含水量离散区间的非饱和渗透系数K(θd)和含水量离散区间的毛细压力Ψ(θd)，计算公式为：

(7)

式中：θL、θR分别为含水量离散区间内左、右边界处的含水率，Δθ=θR-θL。

K(θd)值也可以用公式(5)来计算，即用K(θk)和K(θ)分别代替Ψ(θk)和Ψ(θ)计算K(θd)。

根据Green-Ampt方程，当下渗初始含水率θ0增加到θd时，得到湿润锋在k含水量离散区间内Z方向的速度：

(8)

式中：θd为含水量离散区间中最大孔隙的土壤含水率，m3/m3；θ0为土壤初始含水率，m3/m3；K(θd)为θd含水量离散区间的非饱和渗透系数，m/s；Ψ(θd)为θd含水量离散区间的毛细压力水头，m；Zk为k含水量离散区间的湿润锋到地表面的饱和部分的长度，m。

水分首先从含水率最低的含水量离散区间内计算，含水量从地表到地下水阶段中离散区间非饱和时，得到任意含水区湿润锋垂向计算公式为：

(9)

随着θ的增加，非饱和渗透系数会快速增加，由于含水量离散区间施加的毛细管的作用，水分再分配由右向左开始，水的再分配按照N个含水量离散区间中Ψ(θj)值进行，得到湿润锋位移横向再分配公式为：

(10)

2.3 河道汇流滞后算法

时空变源混合产流模型中河道汇流主要通过将河网概化为河道、水库和节点3部分来实现，河网的每个部分接受来自其他部分的3个上游入流。河道还可以接受来自坡面的侧向入流。将上游汇流单元的入流和单位河道长度侧向入流作为河道汇流的输入，而在河道汇流中往往出现降雨过程线和出流过程线存在一定的时间差，因此，SVSM采用汇流滞后算法进行校正。河道出流过程不发生坦化现象，形状与入流过程相同，只是往后平移了Lag(Lag为滞后时间)时间，滞后算法的出流和入流的关系如下：

(11)

式中：Qt为t时刻河道的出流量，m3/s；It为t时刻河道的入流量，m3/s；Lag为滞后时间，h。

3 研究区域与数据来源

3.1 研究区概况

本研究数据基于中国水利水电科学研究院“全国山洪灾害防治项目”数据集，选取剑阁、汉王场、临洮、尧甸4个小流域作为研究区，各研究区域处于不同的气候分区。

剑阁流域和汉王场流域处于四川盆地，属于亚热带季风气候，雨量充沛，年降雨量超过800 mm，受到地理位置和多变地貌影响，垂直气候明显，其中四川剑阁流域面积245 km2，属于岷江水系的支流, 汉王场流域面积401.28 km2，属于长江干流；甘肃临洮流域和尧甸流域处于黄土高原和青藏高原交汇地带，气候兼具温带大陆性气候和地形气候的特征，冬、春季少雨干旱多风沙，夏、秋季高温多暴雨，年降雨量小于400 mm，山地气候垂直变化明显，其中临洮流域面积585 km2，属于黄河干流；尧甸流域面积463 km2，属于黄河干流。4个小流域的地理位置及水系分布如图2所示。

图2 研究区4个小流域的地理位置及水系分布

3.2 土壤数据和土地利用数据

从全国山洪数据集中找出各个流域的土壤数据绘制土壤质地类型分布图，如图3所示。由图3可看出，四川省两个流域土壤主要以壤土和黏土为主，甘肃省两个流域土壤主要以沙壤土和黏壤土为主。时空变源混合产流模型采用标准化国际制土壤类别。为适用模型需要将各个流域的土壤分为12个大类，按照饱和导水率和土壤田间持水量将12类土壤类型分为易下渗型、较易下渗型、中等下渗型和不易下渗型，然后划分出地貌水文响应单元类型。

图3 研究区各小流域土壤质地类型分布

采用Landsat/TM影像，通过目视解译和非监督分类等方法生成基本的土地利用类型图，包括阴影、林地、草地、农田、房屋建筑等5类。在GIS中通过图层叠加和裁剪的方式得到研究流域的土地利用类型图，见图4。由图4可见，四川省两个流域土地利用类型主要为林地和草地，甘肃省两个流域土地利用类型主要为耕地。土地利用数据和土壤质地数据作为提取参数的输入文件导入山洪分析软件。

图4 研究区各小流域土地利用类型

3.3 降水流量数据

降水和流量基础数据来源于中国水利水电科学研究院。将四川省和甘肃省的4个小流域的降水和流量资料进行整理核对，剔除不符合要求的资料，通过降水和流量插值小程序，将整理好的时间序列数据按照插值程序要求格式订正进行插值，转换为逐小时站点数据序列，然后将数据导入洪水分析软件，利用反距离平方加权法将各个雨量站数据转换为流域面雨量数据,分布于各个子流域[24]。

4 结果与分析

4.1 时空变源混合产流模型构建

为了构建暴雨模式下的变源混合产流模型，必须提供模型输入数据。时空变源混合产流模型的输入包括降水量、最低和最高气温、短波太阳辐射的时间序列值，输入数据的时间序列为小时步长。时空变源混合产流模型以地貌水文响应单元为计算单元，在每一个子流域单元考虑的产流过程包括截流填洼、流域蒸发、土壤蒸发、下渗、超渗、蓄满、优先流、壤中流和地下水径流。

时空变源混合产流模型的模块包括雨量流量插值模块、截流模块、包气带土壤下渗代理计算模块、混合产流模块、坡面径流和河道汇流模块等。在考虑DEM且土壤和植被覆盖空间变异性前提下，将小流域离散成地貌水文响应单元计算，并进行参数自动提取和分配。模型评价指标选取洪峰误差百分比、洪量误差百分比和纳什效率系数3种。

采用SCE(shuffled complex evolution)全局搜索算法对参数进行率定。该方法将确定性的复合搜索技术和生物竞争进化原理结合，引入种群的思想，很好地解决了非线性约束最优化问题[25]，是针对水文模型率定的高效算法。算法基本思路为给定维数，将样本点构成一个子代，子代向多个方向进行搜索，达到给定的维数后，群体向最优点收敛，得到最优的参数范围，然后进行多次迭代，最终产生新的若干数据取值。

4.2 时空变源混合产流模型参数率定结果

采用SCE全局搜索算法进行参数自动率定，考虑目标函数即纳什效率系数，在参数取值范围内对关键参数进行多次迭代寻找最佳参数，然后手动修改参数范围确保参数在取值范围区间内，得到的4个流域模拟较好的一套参数值见表2。

通过模拟得到的汉王场流域1994-2004年、剑阁流域1966和1988年、临洮流域1976年、尧甸流域1976年的逐小时流量数据与通过插值得到的逐小时实际数据绘制模拟流量过程线和实测流量过程线以及相应的降雨过程如图5所示。

由图5可看出，4个流域的水文站流量模拟值与实测值吻合较好，与降雨量情况基本拟合。然而对4个流域进行细致比较发现，汉王场流域的洪水最大峰值模拟结果不如其余3个流域，临洮流域洪水产生和退水过程与实际观测出现时间上的偏差，尧甸流域1976-08-04场次洪水在退水过程中出现二次洪峰，与实际观测不符。之所以导致这样的误差，可能的原因是汉王场流域年降雨量大于800 mm，土壤湿度大，体现在参数中则为表层土壤厚度小；临洮流域由于降雨集中且历时短，汇流滞后现象弱，因而出现模拟洪峰历时比实际缩短，尧甸流域由于降雨持续，且无优先流面积，因此出现二次洪峰。

图5 时空变源模型率定期模拟和实测流量过程线及相应降雨过程

基于参数的物理特性和功能属性，将参数进行如下分组：不透水面积比、土壤深度和表层土壤深度为超渗产流响应参数；优先流面积比、优先流出流线性和非线性系数为蓄满产流响应参数；饱和导水度和洼地平均蓄水量为下渗参数；其余参数为土壤排水参数。根据参数分组对照表2可看出，四川省两个小流域以壤中流占主导；甘肃省两个小流域以超渗产流占主导。

表2 各研究流域时空变源混合产流模型率定关键参数汇总表

对率定期的各类产流量和评价指标进行统计分析，结果见表3。由表3可见，参数率定期内的各个流域洪水场次洪峰误差百分比小于15%，合格率超过90%，纳什效率系数大于0.8，由此说明，本研究的数据预处理是合理的，水文模型估算的模拟参数准确有效，模型在4个小流域内的模拟结果良好。

表3 各研究流域各类产流量及特征值模拟结果汇总

4.3 参数敏感性分析

敏感性分析用于研究模型输入的不确定性对于模型输出不确定性的贡献，模拟中不仅要考虑率定结果，而且要保证率定参数在合理的取值范围内，因此需要对参数分布和敏感性进行研究。本研究选取率定的14个关键参数进行敏感性分析，根据参数取值做出参数敏感性对照图，结果见图6。

由图6可见，尧甸流域的主要敏感参数为ks、wfd和Epan，且wfd>Epan>ks，说明尧甸流域产流机制主要受饱和水力传导度、洼地平均蓄水能力和蒸发能力影响；临洮流域的主要敏感参数为ks和wfd，说明其产流机制受饱和水力传导度和洼地平均蓄水能力两个参数的影响较大；剑阁流域和汉王场流域的主要敏感参数为perip、kp2、kp1和Hl，即其产流机制受不透水面积比、壤中流和优先流出流系数以及表层土壤厚度的影响较大。

图6 模型模拟的各流域参数敏感性占比分析结果图7 各流域产流成分汇总占比分析结果

4.4 主导产流成分

在现有不同地貌特征暴雨水文过程经验的基础上识别小流域主导产流成分。模型将地貌水文响应单元下的产流成分分为超渗产流、蓄满产流、壤中流和基流。鉴于模型将土壤划分为浅层土壤和深层土壤，在不同土壤层上受到重力和毛细力等的作用，则将壤中流细分为快速壤中流和慢速壤中流。

根据模拟结果得到各个流域的产流主导成分，如图7所示。

由图7可得出，四川省的两个流域洪峰主要由蓄满产流和壤中流叠加构成，其中蓄满产流占比最大，洪水地表产流主要为蓄满的形式；甘肃省两个流域洪峰由超渗产流和壤中流叠加形成，且两个流域超渗产流占比达70%～80%，该两流域洪水地表产流主要为超渗的形式。

5 讨论

5.1 模型预报检验

为验证时空变源混合产流模型在山丘区小流域的适用性，需要测试模型对于小流域洪水场次的预测精准性[26-27]。

本文采用率定后得到的参数，通过模型参数文件替换方案和4个小流域水文站获得的实测数据进行模拟结果的验证。图8给出了4个小流域验证期的模拟水文过程线、实测水文过程线以及相应的降雨过程。

图8 时空变源模型验证期各流域模拟和实测流量过程线及相应降雨过程

由图8可看出，在验证期间内，汉王场流域模拟洪峰与实测洪峰不完全相符，剑阁流域洪峰出现时间晚于实测数据，原因可能均为在调节参数过程中将线性出流系数调节过小和汇流参数调节不当所致。相关研究也表明，垂向混合法中对模型影响最敏感参数为汇流参数和土壤排水参数[17,23]。

临洮流域和尧甸流域的模拟峰现时间与实测相吻合，但退水曲线与实测数据不完全匹配，其原因是实测数据为通过插值获得的数据，研究表明降雨数据插值算法会影响降雨量精度。

由《水文情报预报规范》(GB/T 22482—2008)[28]可知，洪水预报中以实测洪峰流量的20%作为许可误差,本文中各个场次洪峰误差百分比均小于20%,且对于验证水文模拟结果的纳什效率系数(NSE)来说,全部场次的NSE值均大于0.8,预报合格率大于90%，因此本文采用的时空变源混合产流模型来计算洪水是可行的。从断面模拟过程线分析出，本模型对于多峰值洪水的模拟效果良好，生成的参数文件能够推广应用于无资料地区洪水的预报中。

5.2 下垫面特征信息对山洪特性的影响

5.2.1 土地利用类型对产流特征的影响山丘区流域下垫面特征复杂，其土地利用类型对于暴雨条件下产流的响应状况不同[29]。根据模型模拟出的4个小流域主要土地利用类型对产流成分占比的影响见表4。由表4可看出，土地利用类型中耕地面积占比最大的流域其超渗产流占比也最大。土地利用类型中有林地时，由于林地对于洪峰的削减作用，产流成分中壤中流占比最大。流域中耕地面积直接影响了地下水储量，而基流主要由地下水补给，因而耕地面积的增大导致产流成分中基流减小。相关文献研究也表示，耕地在涵养水源方面的能力明显弱于林地[30]。