窟野河流域基于物理成因概念的次降雨径流产沙统计模型研究

苏婷婷，吴志勇，何 海，孙昭敏

（河海大学水文水资源学院，南京 210098）

据研究，林草过牧滥伐、城镇无序扩张等不合理的人类活动是引起土壤侵蚀的主要原因，适度控制人类活动、提高植被覆盖率是土壤侵蚀治理及水土流失防治的重要手段［1，2］。20 世纪 80 年代中期以来，由于黄河流域降雨等自然因素和水利、水土保持工程等人类活动的影响，黄河天然径流量和来沙量急剧减少，水土流失得到改善，同时水沙关系发生重大变化，产生一系列新问题，如下游河槽萎缩，河道输沙能力显著降低，威胁流域的防洪安全，影响区域经济的可持续发展。治理水土流失事关经济社会可持续发展和中华民族长远福祉，因此，研究黄河流域侵蚀产沙是新时期黄河治理亟需解决的基础问题。

现有研究主要关注的是黄河沙量锐减的原因及未来某一降雨情景下的多年平均来沙量预测［3-6］，对场次暴雨可能来沙量预测却少有研究。侵蚀产沙模型是研究流域侵蚀产沙的重要工具。目前，流域产沙模型主要分为经验统计模型和物理成因模型。物理成因模型从大量观测和定量描述入手，基于物理基础建立数学模型模拟土壤侵蚀过程，但物理成因复杂，加上水沙动力过程的复杂性，模拟结果依然存在诸多不确定因素。与物理成因模型相比，经验统计模型结构简单、使用方便，仍然是目前流域侵蚀产沙研究的主要工具［5］。黄土丘陵沟壑区地貌复杂，真实下垫面参数获取难度大［6］，众多学者根据各研究区域的产输沙特点，提取影响土壤侵蚀和流域产沙的主要影响因素，建立了不同的统计模型［4，7-17］。如，江忠善等［7］基于陕北、晋西、陇东南黄土丘陵区沟道小流域实测资料，将洪水径流总量、流域平均坡度、黄土中沙粒和粉粒含量、植被作用系数作为产沙量预报的因子；牟金泽等［8］在岔巴沟流域主要考虑洪水径流模数、洪峰模数、主沟道平均比降及流域长度来构造洪水产沙量综合预报公式。本研究以窟野河流域为对象，利用流域1961—2016年的水文泥沙数据分析水沙关系变化，划分不同下垫面代表时期，确定次含沙量的影响因子，分别构建不同下垫面代表时期的流域次降雨径流产沙经验模型，从不同产沙水平及不同降雨等级2 种角度分析模型模拟效果，为客观认识黄河流域侵蚀产沙特点提供科学支撑。

1 研究区及资料

窟野河是黄河中游水土流失严重的一条多沙粗沙支流，发源于内蒙古南部沙漠地区，于陕西省神木县汇入黄河。干流全长242 km，流域面积为8 706 km2。流域多年平均降水量为410 mm，多年平均径流量为7.47 亿 m3，年均输沙量为1.36 亿t，占黄河总输沙量的6.9%。

该流域气候属寒温带干旱、半干旱大陆性季风气候，降水多为暴雨形式，7—9 月降水占全年的68%。地质地貌为沙质（流域西部及偏北）、砾质（流域北部及东北部）及黄土丘陵（南部，典型黄土丘陵沟壑地貌景观）的组合。

流域内共设有38 个雨量站和6 个水文站，出口控制站为温家川站。水文资料采用窟野河流域1961—2016 年降水摘录表、洪水要素摘录表，包含时段降水量、流量、含沙量等数据，由黄河水利委员会水文局提供。降水量摘录表只记录汛期数据，主要集中在6—9 月。水沙资料主要源于流域出口控制站温家川站实测径流泥沙资料。流域概况见图1。

以窟野河流域温家川站1961—2016 年夏汛期（6—9 月）的场次洪水为研究对象，洪峰流量大于100 m3/s 的洪水全部入选。共计收集到97 场洪水资料，其中参与建模的1961—1979 年共23 场洪水，1980—1996 年共 20 场洪水，1997—2016 年共 30 场洪水，未参与建模的共24 场洪水，用于验证模型模拟效果。

图1 窟野河流域概况

2 次降雨产沙经验模型构建

2.1 不同下垫面代表时期识别

采用降雨径流泥沙双累积曲线法分析流域水沙序列突变点，识别不同下垫面的代表时期。双累积曲线法是利用同期累积降雨量与累积径流量（累积输沙量）的关系曲线斜率变化分析水沙序列的变化趋势，斜率发生变化表明由于人类活动或环境因素导致下垫面产水产沙关系发生改变，转折点为突变点，以突变点为准划分为不同的时期。

图2 为场次洪水累积雨量-场次洪水累积沙量双累积曲线，图3 为年径流量和年输沙量的累积距平曲线，2 条曲线均在1979 年和1997 年出现转折点，以此将径流和输沙序列划分为3 个阶段。将1961—1979 年作为天然时期，1997 年之后流域产沙环境大体稳定，故1997—2016 年作为现状下垫面代表时期。据研究，黄河流域产沙环境在20 世纪70 年代以前持续恶化，70 年代以后逐渐好转，1997 年以后显著改善，因此1997 年之后产沙量和产沙强度大幅降低［15］，此结论与本研究相符。

图2 窟野河流域场次洪水累积雨量-场次洪水累积沙量双累积曲线

图3 1960—2016 年窟野河流域年径流量和年输沙量累积距平曲线

2.2 水沙关系分析

降雨产生径流，径流是泥沙的载体。将流域内场次降雨面雨量与次洪沙量一一对应，生成面雨量和次洪沙量关系散点图（图4）。通过对比1961—2016 年3 个不同时期面雨量和次洪沙量关系散点位置，可直观反映流域的面雨量和次洪沙量关系的变化。若面雨量和次洪沙量关系散点位置偏上，在面雨量相同的前提下，则表明该时期的输沙量较多，反之则较少。分析不同时期面雨量与次洪沙量关系（图4）可知，随时间推移，点群向右下偏移，呈现3 个带状分布，说明在现状下垫面条件下，同样量级的降雨，与天然时期相比产沙量明显减少。以2012 年7月和2016 年8 月发生的2 场洪水为例，次洪沙量点据明显向右下偏离历史洪水形成的趋势线。

图4 不同时期面雨量与次洪沙量的关系

分析1961—2016 年3 个不同时期次洪水量与次洪沙量关系（图5）可知，不同时期次洪水量与次洪沙量关系发生明显变化，点群逐渐向X 轴方向偏移，近年来相同量级场次洪水的含沙量与天然时期相比显著减小。

图5 不同时期次洪水量与次洪沙量的关系

2.3 次洪含沙量影响因素分析

由次洪输沙量的定义可知，径流事件中一次暴雨径流的水沙有如下关系：

式中，M为次洪产沙量（kg），C为次洪含沙量（kg/m3），Q为次洪径流量（m3），等式两边同除以流域面积得：

式中，Ms为输沙模数（t/km2）；Cs为次洪平均含沙量（kg/m3）；H为径流深（mm）。

流域侵蚀产沙是个复杂的物理过程，是流域下垫面包括地貌、土壤、植被、水保措施等以及降水因子综合作用的结果，而流域侵蚀产沙会输出径流量和次洪含沙量，因此次洪平均含沙量和径流深可以完整地表达流域产沙［18］，所以式（2）利用径流深和影响次洪含沙量的因子的线性关系来模拟产沙是一种可行的研究方法。

次洪含沙量与次洪水流挟沙能力和流域可被携带土壤分散量有关。由于黄河流域沙源丰富，故水流挟沙能力决定了次洪含沙量的大小。水动力学研究表明，水流挟沙能力与流量关系密切［19-22］。参与相关分析的因子有反映洪水径流特征的平均流量、洪峰流量、径流总历时及反映径流输沙特征的最大含沙量。采用1997—2016 年30 场洪水参与相关分析，采用SPSS 软件进行相关分析，结果（表1）表明，次洪含沙量与洪峰流量、最大含沙量、径流总历时在0.01 水平上显著相关，故选用这3 个要素指标和径流深模拟产沙。

表1 窟野河流域次洪平均含沙量与其影响因素的相关系数

2.4 模型构建

基于公式（2），采用径流深和与次洪含沙量相关的最大含沙量、洪峰流量、径流总历时的非线性函数形式模拟产沙，利用SPSS 软件中的非线性回归过程模块确定输沙模数与各因子之间的关系，建立流域输沙模型。

由输沙模数与影响因素相关关系（图6）可知，输沙模数与各影响因素间并非简单的线性相关。输沙模数与最大含沙量的对数形式拟合较好，与洪峰流量和径流总历时的多项式形式拟合较好。基于输沙模数影响因素相关关系及流域输沙关系表达式假设输沙模数与各因子间的模型表达式为：

式中，Ms为输沙模数（t/km2）；H为径流深（mm）；Cm为最大含沙量（kg/m3）；Qm为洪峰流量（m3/s）；T为径流总历时（min），a1、b1、b2、b3为模型参数。

图6 窟野河流域输沙模数与其影响因素的关系

建模具体步骤如下：①为所有未知参数指定一个初始值，将原方程按泰勒级数展开，并只取一阶各项作为线性函数的逼近，其余各项均归入误差；②采用最小二乘法对模型中的参数进行估计，用参数估计值替代初始值，得到1 个新的曲线函数；③将新得到的曲线函数展开，线性化，求出参数估计值；④如此反复，直至参数估计值收敛，得到最优解。

基于现状下垫面代表时期（1997—2016年）30场洪水的模拟结果见表2。由表2可知，Ms与H、Qm的拟合关系式的R2=0.960，均高于其余2个模型，可见该模型能解释96.0%的变异，模型的拟合效果较好。最终确定模型结构为Ms=109.174H+0.113HQm-10.225。

表2 基于现状下垫面（1997—2016 年）代表时期的产沙模型拟合关系式

以同样的方法在1961—1979 年（共挑选23 场洪水参与建模）、1980—1996 年（挑选20 场洪水参与建模）分别构建模型，模型表达式见表3。

表3 不同时期产沙模型拟合关系式

3 模型验证

3.1 不同产沙水平模拟效果分析

为了验证不同产沙水平的模型模拟效果，从窟野河1961—2016 年的97 场长系列输沙资料中选择未参与模型构建的24 场洪水，其中1961—1979 年11 场，1980—1996 年 7 场，1997—2016 年 6 场。以常规分类方法将产沙水平分为3类，分别为<500 t/km2、500～1 000 t/km2、>1 000 t/km2，验证模型对不同产沙水平次降雨产沙量的模拟效果，模拟结果如图7 所示。由图7 可知，模型在不同的产沙水平的模拟误差上表现出一定的变异特性，总体而言，平均相对误差绝对值为30.6%，相对误差最小达9.7%，随着产沙水平的增大，模型模拟相对误差逐渐减小。3个时期不同产沙水平的平均相对误差绝对值见表4，模拟规律与长系列输沙资料一致，具体表现为随着产沙水平的增大，模型模拟相对误差逐渐减小。上述规律说明小产沙事件与大产沙事件的产沙机制有所差异，由于黄土高原丘陵沟壑区高含沙水流的存在，使得不同产沙水平下流量含沙量之间的作用发生变化。综上可知，模型对大产沙事件的模拟效果优于小产沙事件。

图7 1961—2016 年窟野河流域24 场不同水平产沙事件与产沙模型的相对误差

表4 不同时期不同产沙水平模型模拟平均相对误差绝对值

3.2 不同降雨等级模型模拟效果分析

为验证所建立的次降雨产沙经验模型模拟不同等级降雨侵蚀产沙的效果，从窟野河1961—2016 年的97 场长系列输沙资料中选择未参与模型构建的24 场洪水，其中 1961—1979 年 13 场，1980—1996 年6 场，1997—2016 年5 场。按不同降雨等级将洪水进行分类分析，分类标准：中小雨，1 d（24 h）降雨量<25 mm；大雨，1 d（24 h）降雨量在25～50 mm；暴雨，1 d（24 h）降雨量在50～100 mm。用3 个不同时期的模型分别模拟对应时期的场次洪水产沙量，将实测值与模型模拟值作对比，模拟相对误差见图8，平均相对误差见表5。总体而言，模型模拟平均相对误差为32.8%，相对误差最小达0.03%，模型对大雨及暴雨的模拟效果略优于中小雨。以19910720 次（1991年7 月20 日）洪水为例，本时期模型模拟输沙模数为273 t/km2，1960—1979 年模型模拟输沙模数为725.3 t/km2，同样等级的降雨不同时期模型模拟值的差异除了模型本身误差，可认为是不同时期下垫面差异导致的。

图8 1961—2016 年24 场不同降雨等级典型产沙事件模拟结果相对误差

表5 不同时期不同降雨等级模型模拟平均相对误差绝对值

3.3 不同时期模型参数差异分析

不同时期次降雨径流产沙模型参数差异的原因：①分析输沙模数与模型自变量（洪峰流量、径流深）在1961—1996 年和1997—2016 年的变化趋势可知（图9），1961—1996 年洪峰流量、输沙模数及径流深总体变化不大，1997 年之后显著减小，变化明显；②分析不同时期的模型参数，发现1961—1979 年与1980—1996 年的模型参数相差无几，1997—2016 年模型项的系数明显增大，导致输沙模数急剧减小。据研究［23］，窟野河流域 1979 年以前，年径流量和输沙量处于上升趋势，该阶段受人类活动影响较小；1980—1996 年，年径流量和输沙量呈微弱下降趋势，主要是20 世纪80 年代开展了水土保持治理及80 年代末开始的煤矿开采等造成的影响；1997 年以后，年径流量和输沙量大幅下降，究其原因主要是20 世纪末大规模的煤矿开采对流域水文地质造成一定的影响，进而影响径流和输沙过程。

图9 输沙模数、洪峰流量及径流深在不同时期的变化趋势

4 小结

1）双累积曲线、累积距平曲线均在1979 年和1997 年出现转折点，以此将径流和输沙序列划分为3 个阶段。将 1961—1979 年作为天然时期，1997 年之后流域产沙环境大体稳定，1997—2016 年为现状下垫面代表时期。

2）研究基于Ms=CsH，采用非线性回归法，构建以径流深与洪峰流量为变量的流域场次洪水产沙统计模型。模型R2均大于0.9，平均相对误差绝对值在30%左右；通过模拟不同时期不同产沙水平的场次洪水可知，随着产沙水平的增大，模型模拟相对误差逐渐减小，模型对大产沙事件的模拟效果优于小产沙事件；通过模拟不同等级降雨侵蚀产沙发现，模型对大雨及暴雨的模拟效果略优于中小雨，为客观认识黄河流域侵蚀产沙特点提供科学支撑。