清水河流域场次洪水输沙特性及关键影响因素

付思佳, 张会兰,2, 杨 军

(1.北京林业大学 水土保持学院, 北京 100083; 2.重庆缙云山三峡库区森林生态系统国家定位观测研究站, 重庆 400711; 3.广东省水利水电科学研究院, 广州 510610)

由于气候变化和人类活动的综合影响,流域内降水和产流产沙发生显著变化,导致流域水沙关系发生变化[3],尤其是因暴雨产生的洪水事件。探索洪水过程中的降雨过程、泥沙输移、径流特征及其相互关系,可以为研究流域的洪水和泥沙产生规律以及洪水与泥沙的关系提供重要依据。对于雨洪式河流而言,一场或者几场暴雨事件可能贡献流域绝大部分输沙量[4],所以基于次洪尺度的泥沙输移过程及其主导要素也需引起足够重视。因此,近年来,国内外学者针对场次洪水的水沙输移变化规律及其影响因素的研究颇多[5-7]。比如Bussi等[6]应用水沙关系曲线分析了泰晤士河的悬移质泥沙输移的时空变化,研究结果发现夏季后的第一场洪水在相同的洪量条件下其输沙量是其他次洪输沙量的1.5～2.0倍,产洪输沙量呈现下降的趋势。在暴雨引发的洪水过程中,降水和地表径流为流域泥沙输移过程中主导驱动要素[7]。

清水河流域是黄河在宁夏回族自治区的最大的一级支流,也是宁夏段重点产沙区和侵蚀严重区之一。为保障黄河安澜,加快水土流失综合治理和削减入黄泥沙,20世纪70年代以来,该区域开展了一系列生态工程以及水土保持措施[8]。2000年以后,由于人类活动的增加,清水河流域的径流泥沙产量和水沙关系发生明显变化。李颖曼等[9]对清水河流域近60 a降雨量和入黄泥沙量进行研究,发现流域年降水量分布不均,1990—2000年泥沙量最大,2010年后显著减少。马进国等[10]发现清水河的年内水沙多集中在汛期,且年际水沙量相差悬殊。鉴于清水河流域降水年际年内分布极为不均的气候条件以及丘陵沟壑区沟壑纵横的独特侵蚀地貌,其年土壤侵蚀多由少数几场暴雨洪水主导,年际尺度或者月尺度的研究会弱化甚至掩盖了次洪过程产输沙规律以及动态特性。目前对清水河流域内水沙过程的研究仅聚焦于年际或者年内尺度水沙演变过程,涉及深入次洪尺度的研究相对缺乏,所以明晰清水河流域次洪尺度的径流侵蚀输沙动态以及解析其驱动机制,深入理解流域水沙关系及成因,可以进一步细化对水沙机理的认知,并有利于丰富流域水沙调控基本理论、完善水土保持综合治理方案。因此本文根据清水河流域6个水文站(原州、固原、贺堡、韩府湾、王团以及泉眼山)2006—2016年实测洪水水文数据和表征次洪尺度洪沙输移特征的合理参量,采用随机森林辨析不同区域洪水事件输沙量的影响因素;引入反映单次洪水事件过程属性特征的径流侵蚀功率概念对次洪尺度水沙关系予以表征,于次洪尺度揭示流域径流及输沙特征。

1 研究区概况

清水河流域是宁夏回族自治区黄河上游的一级支流,其发源地为宁夏六盘山北端黑刺沟脑,从中宁县的泉眼山区域流入黄河,该流域是宁夏回族自治区境内最大、侵蚀最为严重的黄河支流,全程长320 km,流域面积为14 480 km2。地跨东经105°00′—107°07′,北纬35°36′—37°37′。清水河流域地势南高北低,流域海拔介于1 182～2 920 m,相对高差达1 738 m,河道平均比降为1.49‰。清水河流域地理位置、水系分布、6个水文站(原州、固原、贺堡、韩府湾、王团以及泉眼山)分布如图1所示。以6个水文站点为出口站,将清水河流域划分为6个呈嵌套关系的子流域,在本文中子流域的面积用水文站的控制面积表示,分别为:原州105 km2,固原105 km2,贺堡200 km2,韩府湾4 935 km2,王团2 430 km2,泉眼山7 115 km2。

图1 清水河流域地理位置Fig. 1 Location of Qingshui River Basin

清水河流域平均降水量为349 mm,时空分布不均。年降水量最高可达611 mm,最低仅有195.8 mm,极值比达3.0[11];年内70%～75%的降雨集中于汛期[12]。空间上,呈现自上游向下游递减的特征,年际和年内降雨变异大。清水河多年均径流量为2.16亿m3,年输沙量4 940万t,年均含沙量229 kg/m3,流域内水沙特征表现为产流少、产沙多、土壤侵蚀情况严峻[13]。

2 数据与方法

2.1 数据来源

清水河流域6个子流域控制水文站的水沙数据来自《黄河流域水文资料》,其中原州站采用2009—2013年、2013—2016年的次洪尺度径流、泥沙数据;固原站采用2006—2008年的次洪尺度径流、泥沙数据;贺堡站采用2006—2013年、2013—2016年的次洪尺度径流、泥沙数据;韩府湾采用2006—2013年、2013—2016年的次洪尺度径流、泥沙数据;王团采用2015—2016年的次洪尺度径流、泥沙数据;泉眼山采用2006—2013年、2013—2016年的次洪尺度径流、泥沙数据。水文监测、取样以及实验室测定的过程均完全遵循国际方法和标准,并且在发布前经过严格验证。

2.2 研究方法

2.2.1 有效洪水场次筛选方法 依据暴雨—洪水—输沙的逐时过程线选取具有明显涨落段、峰顶、峰谷的完整洪水过程。为剔除低输沙率洪水场次事件,本研究根据水文监测标准,将一次洪水事件定义为径流深超过0.05 mm,洪峰流量大于0.1 m3/s,洪水历时大于150 min的产流产沙过程。本研究选取了6个水文站的155场暴雨洪水事件,对次洪尺度的产流输沙过程特征及其影响因素进行分析。

2.2.2 次洪过程径流及输沙量指标计算方法 本研究选用8个径流泥沙指标表征流域内单次洪水事件内的径流和输沙特征:洪水径流深(H,mm)、输沙模数(SSY,t/km2),洪水历时(T,min)、洪峰流量(Qp,m3/s)、径流变率(FV)、平均径流量(Qm,m3/s)、平均含沙量(SSCm,kg/m3)、最大含沙量(MSSC,kg/m3)。

对于某次特定洪水事件,若水文站所观测的时间间隔为Δt,瞬时流量和含沙量分别为Qt和SSCt,水文站的控制面积为A。单次洪水事件输沙量SY(t)计算如下所示:

由图2(a)可知，随着粒子半径的增加，1.064 μm和3.8 μm的消光效率因子上下振荡逐渐趋近于2，10.6 μm的消光效率因子先增加后减小逐渐趋近于2.由图2(b)可知，整体上看三波段的消光截面，均随着半径的增加而增加且数值十分接近，但由图2(c)可知，1.064 μm和3.8 μm的消光截面曲线波动式上升，且三波段消光截面数值存在一定差别.

(1)

输沙模数SSY可用下式进行计算:

(2)

洪水径流深H(mm)的计算公式如下:

(3)

洪水事件径流变率FV可通过下式计算:

(4)

式中:Δt为水文站观测时间间隔(min);Qt为单次洪水事件的瞬时流量(m3/s);SSCt为单次洪水事件瞬时含沙量(kg/m3);A为水文站的控制面积(km2);Qp为洪峰流量(m3/s);Qm为单次洪水事件的平均流量(m3/s)。

2.2.3 径流侵蚀功率计算方法 本研究采用鲁克新等[5]提出的径流侵蚀功率概念作为单次洪水事件下的侵蚀产沙动力指标,计算公式如下所示:

E=H×Q′p

(5)

式中:E为单次洪水事件侵蚀产沙动力〔mm·m3/(s·km2)〕;H为单次洪水事件径流深(mm);Q′p为单次事件洪峰流量模数〔m3/(s·km2)〕,其大小等于单次洪水事件内的洪峰流量Qp(m3/s)与水文站控制面积A(km2)的比值。

2.2.4 随机森林 洪水事件输沙量受不同径流因子指标的耦合作用。本研究采用随机森林算法来量化影响流域内输沙量的驱动因子,其计算过程基于R语言环境调用Random Forest函数包实现。随机森林是一种基于分类树(classification tree)的统计学习理论[14],主要是从原始样本中采用Bootstrap重抽样方法抽取若干样本然后分别构建分类树,接着组合分类树的预测并通过投票方式得出最终结果。随机森林对于变量之间的多重共线性不敏感,不容易出现过拟合的现象,大量的理论和实证都证实了随机森林算法对异常值、噪声具有很好地容忍度甚至能够解析和预测高达上千条解释变量的作用,所以被应用于诸多领域[15]。次洪尺度流域输沙量的驱动因子涉及多指标变量,随机森林算法对多变量数据处理具有优越性能,可以用于各指标重要性评价中。

3 结果与分析

3.1 洪水产流输沙特征

图2为清水河流域事件输沙模数(SSY)量级分布及其所占比例。根据输沙模数将单次洪水事件划分为5个等级[16],即特大量级洪水事件(SSY>10 000 t/km2)、大量级洪水事件(SSY介于5 000～10 000 t/km2)、中量级洪水事件(SSY介于1 000～5 000 t/km2)、小量级洪水事件(SSY介于100～1 000 t/km2)以及特小量级洪水事件(SSY<100 t/km2)。从图2可以看出,清水河流域内并无特大量级以及大量级洪水事件发生,原州、固原以及贺堡特小量级洪水事件发生最为频繁,比例分别为86%,100%,67%,而贺堡中量级洪水事件发生次数也较多,占比为19%。而韩府湾、王团以及泉眼山特小量级、小量级以及中量级洪水事件均有分布,分布比例均在30%附近。由此可以看出,流域基于事件的产洪输沙呈现显著的空间分异特性,即从上游到下游径流泥沙呈现递增的趋势,位于上游的原州以及固原子流域洪水事件主要为特小量级洪水,其径流输沙量较小,显著低于中下游区域;而位于中下游的韩府湾—王团区间以及王团—泉眼山区间中等量级洪水发生相对频繁,为流域洪水事件产洪输沙最为强烈的区域。

图2 2006-2016年清水流域洪水事件输沙模数量级分布及其比例Fig. 2 The frequency distribution of flood specific sediment yield from 2000-2006 in Qingshui River Basin

不同子流域基于事件的洪水径流泥沙特征对比分析如图3所示(不同字母表示为存在显著差异,p<0.05)。从图3可知,韩府湾、王团以及泉眼山的洪水历时(T)显著高于原州等3个子流域;韩府湾的洪水径流深(H)显著低于贺堡以及王团,与原州、泉眼山没有显著差异,从全流域尺度上看,流域基于事件的洪水径流深差异相对较小;王团基于事件的洪峰流量(Qp)为最大,原州与固原基于事件的洪峰流量(Qp)显著低于贺堡、韩府湾以及泉眼山;洪水平均流量(Qm)呈现出与洪峰流量相似的空间分异规律;基于事件的径流变率(FV)的空间分异规律与洪水历时的规律相似,即韩府湾、王团以及泉眼山的径流变率显著高于原州等3个子流域,已有研究表明,对于输沙量级较大的洪水事件,其输沙量主要根据径流总量决定,对其而言水沙关系已经相对趋于稳定[17],所以径流变率相对较小,但是对于小型洪水事件频繁的原州、固原以及贺堡,其水沙关系相对多变不稳定,所以导致其径流变率显著低于韩府湾等子流域;基于事件的输沙模数(SSY)的空间分异特征与径流深的相似;基于事件的平均含沙量(SSCm)呈现从上游向下游递增的趋势,王团的基于事件的含沙量要显著高于其他5个小流域,而最大含沙量(MSSC)与平均含沙量空间分异规律相似。

图3 清水河流域洪水事件的径流输沙空间特征Fig. 3 Spatial characteristics of runoff and sediment transport during flood events in Qingshui River Basin

3.2 基于事件的输沙量驱动要素辨析

对各个子流域洪水事件驱动的输沙量和上述基于事件的径流泥沙指标进行随机森林分析,用以量化各指标对输沙量的影响(图4),横坐标为重要性,反映响应变量受各解释变量的影响大小,其数值越大表征对响应变量的影响程度越大,同时计算出各随机森林模型的拟合精度(R2)以及对全模型进行显著性检验。本研究中所有子流域随机森林模型的拟合效果较好,其拟合精度均超过0.59,并且所有随机森林模型全部通过显著性检验(p<0.001)。由图4可知,各子流域的产洪输沙特征主要影响因素不同,主要体现为在小量级洪水频发的原州、固原以及贺堡子流域解释变量重要性前三名分别为径流深(H)、洪水历时(T)和洪峰流量(Qp),而在中量级洪水频发的韩府湾、王团以及泉眼山子流域解释变量重要性前三名分别为径流深(H)、洪峰流量(Qp)和平均流量(Qm),受洪水历时(T)的影响较小。

图4 清水河流域不同流域基于事件的输沙量驱动要素重要性排序Fig. 4 Importance ranking of event-based drivers of sediment in the different sub-basins of Qingshui River Basin

地表径流是引发水力侵蚀以及协同土壤颗粒分离和泥沙搬运的基本驱动力,径流总量(深)与地表径流直接相关,是洪水侵蚀产沙的综合潜力的指标。可以看出,在流域下垫面等基本条件差异较大的情况下,流域次洪过程的输沙量主要受控于洪水径流总量、洪峰流量(Qp)及变率的影响(FV)。对于小量级洪水事件而言,洪水历时(体现为径流侵蚀能量对流域侵蚀产沙的时效性)对其影响较大,而中量级洪水事件几乎不受制于洪水历时的影响。小型洪水事件中,其水沙关系并不稳定,形式多变,而对于中型洪水水沙关系相对趋于稳定,洪水历时的影响相对甚微。这表明,对于中型洪水或者高含沙水流而言,因为其水沙关系相对趋于稳定,输沙量的增大主要原因是径流量的增加。

3.3 基于事件的水沙关系

径流总量(深)可以作为水沙分离土壤颗粒以及搬运泥沙的综合潜力的有效指标,尤其对于量级较大的洪水事件[18]。洪峰流量(Qp)和平均流量(Qm)是基于次洪过程的径流指标,可以有效表征水流的潜在侵蚀能力和挟沙能力的大小,且对基于事件的平均含沙量(SSCm)和最大含沙量(MSSC)具有十分重要的影响。因此,次洪过程中径流所驱动的侵蚀速率是洪水径流总量和水沙强度综合作用后的结果。但是如果仅仅使用径流总量(深)单一指标作为揭示流域洪沙关系以及对侵蚀产沙量进行预测,其结果可能会对洪水事件的产沙量预估产生较大偏差。所以,除径流深之外,还需要引入反映单次洪水事件过程属性特征的多维指标用以表征次洪过程水沙关系以及进行预测泥沙。在解释变量中引入洪峰流量(Qp)作为次洪过程的典型指标,能够在一定程度上消解由径流变率(FV)导致的不确定性,表现出水流造成的土壤侵蚀与泥沙输移的特有过程。径流侵蚀功率(E)是基于径流深和基于过程的变量洪峰流量这两个指标来概化事件尺度的洪水径流过程,其值是由洪峰流量和单次洪水事件径流深的乘积得到,该概念代表了径流深和水流强度的累积效应和瞬时作用强度对次洪水过程中侵蚀产沙的耦合效应,也综合反映了侵蚀过程中洪水侵蚀能量的连续性和不连续性的统一。图5分别采用洪水径流深(H)以及洪水径流侵蚀功率(E)进行水沙关系分析,结果表明径流参量指标与输沙量均呈幂函数关系,同时可以看出采用径流侵蚀功率与次洪尺度的输沙量其拟合优度(R2)要明显高于径流深与输沙量的拟合精度,尤其对于小型洪水,其拟合优度均在0.9之上。相比于广泛使用的单一指标径流深,采用径流侵蚀功率所拟合的决定系数(R2)比传统的径流量(深)-输沙关系较高,可以更好地表征次洪过程中的泥沙输移特征。

注:H为洪水径流深;SY为与输沙量;E为洪水径流侵蚀功率;SY为输沙量。图5 清水河流域不同子流域基于事件的水沙关系Fig. 5 Event-based discharge-sediment relationships of different sub-basins in Qingshui River Basin

3.4 径流侵蚀功率与水沙关系

以图6所示的韩府湾和泉眼山为例,可以看出当径流量低于某一临界值时,含沙量变化很大,径流量与输沙率呈正比;当径流量超过某一临界值时,输沙量变化很小并趋于稳定。目前众多关于黄土高原丘陵沟壑区洪沙关系的研究中,发现次洪事件下的含沙量在流域洪水超过某一固定值时将趋向于饱和,并不会随着洪水量的增加而进一步增加[17],这表明高含沙水流受自身分离能力限制,在侵蚀过程中存在着一定的极限含沙量(侵蚀阈值),并且因所处流域内部位或者环境等的不同而有所差异[18]。在此基础上,Zheng等[19]在研究大理河流域次洪事件下径流量以及含沙量关系时亦指出,洪峰流量通常被视为影响次洪事件输沙量重要驱动因素,但单次洪水的产输沙能力会随着径流量的不断增加而达到一定某一固定值。

图6 清水河流域两个子流域基于事件的平均流量与含沙量的关系Fig. 6 Event-based SSC-Qm sediment relationships of two sub-basins in Qingshui River Basin

基于事件的含沙量及其径流侵蚀功率进行回归分析,其结果表明含沙量和径流侵蚀功率之间存在对数函数关系,这反映出受到水流分离能力的限制,在沙源丰富的区域存在一定的侵蚀阈值现象。Zhang等[20]的研究亦发现这一现象,洪水事件下的径流-输沙系统一旦突破侵蚀阈值状态后,基于事件的含沙量也将趋于稳定。一般而言,当径流侵蚀功率在较低水平时,其所释放的侵蚀能量也相对较低,在此阶段,次洪过程的含沙量随着径流侵蚀释放的能量变大而急速增大。但是一旦径流侵蚀功率位于较高水平时,基于事件的含沙量则维持相对稳定的状态,即不再随着径流侵蚀功率的进一步增加而变大(图7)。刘尚武等[4]在对嘉陵江场次洪水输沙规律的研究中发现径流侵蚀功率与次洪输沙模数之间存在较好幂函数关系,且当径流侵蚀功率逐渐增加,次洪输沙模数的变化量逐渐减小,发生这一变化的主要原因是人类活动影响了次洪输沙量。自2000年以来清水河流域实施淤地坝建设工程及人工林草等水土保持措施[9],影响了流域的水沙关系。在黄土高原地区,因其重力侵蚀等多种侵蚀营力可为侵蚀过程提供充足的泥沙供给,大部分侵蚀产沙事件在全坡面尺度就已经到了极限含沙量,即受分离能力受限的最大含沙量,随后进入流域系统之后含沙量可保持稳定不变,高含沙水流更为如此。综合而言,与传统的径流深、平均流量或者洪峰流量等单一径流参量指标相比,径流侵蚀功率可以反映单次洪水事件过程属性特征,更适宜表征次洪尺度的水沙关系并对次洪尺度下泥沙输移的基本动力来源进行解释。

4 结 论

(1) 根据洪水事件的输沙量级的差异对不同子流域内155场次洪水进行分类可知,流域并未发生特大量级和大量级洪水事件,原州、固原以及贺堡特小量级洪水事件发生最为频繁;而韩府湾、王团以及泉眼山特小量级、小量级以及中量级洪水事件均有分布。流域产洪输沙量呈现显著的空间分异特性,即从上游到下游径流泥沙呈现递增的趋势。

(2) 影响不同子流域洪水输沙过程的主要径流参量有所区别,在小量级洪水频发的原州、固原以及贺堡子流域解释变量重要性前三名分别为径流深、洪水历时与洪峰流量,而在中量级洪水频发的韩府湾、王团以及泉眼山子流域解释变量重要性前三名分别为径流深、洪峰流量与平均流量。从整体上看,随着流域空间尺度的不断推移,次洪过程输沙量的主控因素逐渐趋于单一,并且愈发受到流量过程变化的影响。

(3) 在侵蚀产沙的解释变量中引入表征径流过程的复合指标更能综合体现天然降水和流域下垫面属性对次洪过程侵蚀产沙的综合影响,其拟合精度相较于传统的径流-输沙关系更高。相较于单一径流参量指标而言,径流侵蚀功率更适宜对次洪尺度水沙关系予以表征,更能合理解释次洪尺度下泥沙输移的基本动力来源。