近30年苦水河水沙变化特征及成因分析

傅彦超, 王友胜, 杨 志, 刘 畅, 张晓明,魏天兴, 辛 艳, 马文涛, 解 刚, 任正龑

(1.北京林业大学 水土保持学院, 北京100083; 2.中国水利水电科学研究院流域水循环模拟与调控国家重点实验室, 北京 100048; 3.宁夏回族自治区水土保持监测总站， 银川750002)

黄河流域是水土保持和流域泥沙研究的重点区域，存在水土流失、泥沙淤积和洪水威胁等问题。特别是近几十年来，黄河流域气温升高、极端事件及自然灾害频发和径流减少等，是中国生态安全保障和经济社会发展的重点和难点地区[1-2]。苦水河是宁夏境内重要的黄河一级支流，流域内蒸发强烈，风大沙多，河流含沙量高，是重要的生态敏感区[3]。流域下游为冲积平原，黏土砾石含量较高，上游为干旱草原区，植被稀疏，且地势较高，发生土壤侵蚀之后极易沿河而下汇入黄河干流。不少研究指出[4-6]，位于黄土高原丘陵沟壑区的黄河干支流输沙量均有大幅减少，在黄河中游产沙区人类活动是输沙量减少的主控因素[7]。在黄河宁夏段，上游来水来沙条件的变化是造成宁夏河段河道淤积的原因[8]。但目前气候变化及人类活动对宁夏干旱区域水沙变化的影响尚不明晰。因此，探究苦水河流域的水沙变化对明晰流域内水沙输移特征和探究宁夏干旱区水土保持理论具有重要意义。

水沙序列是通过水文站实际测量得到的径流及输沙非平稳时间序列[9]，针对流域径流、输沙数据的时间序列分析，可以探究径流泥沙在时间尺度上的变化特征和变化趋势。水沙序列的特征可以采用Mann-Kendall非参数检验，来检测流域径流输沙在特定时间段内的变化趋势以及对变化过程中突变点的检测[9-10]。也有学者开展了径流、输沙集中度研究[11]，该概念最早被用来研究年内降水的分配情况[12]，后也有应用在径流、泥沙相关研究当中[13-14]。水沙关系曲线(sediment rating curve，SRC)是研究流域输沙变化的常用工具，用以表达径流量和输沙量的内在关系，区域内水文循环和泥沙输移可以通过水沙关系曲线的变化来进一步表达[15]。水沙关系曲线有多种形式，包括线性函数、幂函数、分段函数等，以幂函数最为常见[16-17]。幂函数水沙关系满足表达式SSC=aQb(SSC为泥沙浓度，Q为流量)，参数a，b能够反映河流的输沙能力强弱和外界环境对其影响[18]。许多研究发现[19-20]，在历时短、降雨强的洪水事件下，河流输沙和径流会发生迟滞现象，洪水场次下的水沙关系曲线会呈现出不同环路曲线(SSC-Q)类型，这可以用来分析河流的输沙状况和水沙协调关系。

本研究以苦水河流域为研究对象，基于1989—2019年郭家桥水文站长时间序列降水及水沙数据，采用Mann-Kendall非参数检验分析近30年苦水河流域水沙变化情势，根据实测资料，对长时间水沙关系曲线及其参数做趋势分析，并在洪水条件下分析水沙关系曲线的滞回特征。基于双累积曲线法、集中度分析和土地利用转移矩阵等统计方法，识别苦水河流域水沙变化的主要驱动因素，在气候和人类活动的影响下，探讨其对苦水河流域水沙变化的贡献率，研究旨为提高对苦水河流域水沙变化规律及驱动要素的理解，并为未来宁夏干旱区优化水土保持措施提供参考。

1 研究区概况

苦水河为宁夏境内黄河第二大一级支流，位于宁夏回族自治区中东部，源于甘肃省环县甜水堡镇，向北经宁夏吴忠市盐池县、同心县、利通区、红寺堡区、灵武市5个县(市、区)，于灵武市新华桥镇华一村汇入黄河，主河道全长223.8 km。流域位于北纬37°03′—38°05′，东经106°02°—107°15′，宁夏境内流域面积4 942 km2，占苦水河流域总面积的94.7%[21]。流域水系主要有苦水河、甜水河、沙沟、小河4条水系。流域地处腾格里沙漠与毛乌素沙地的交汇处，地貌类型多为风沙丘陵。流域内降雨稀少，多发历时短，强度大的降雨，多年平均降水量248 mm[3]。流域年平均气温处于上升趋势，蒸发强烈。土壤类型以淡灰钙土、黄棉土和风沙土为主，有机质含量低，土质疏松。流域内植被类型属于荒漠草原植被，植被种类少，种群结构单一，植被覆盖度低[22]。自2000年，宁夏开始实施退耕还林还草工程，对南部风沙区的沙漠化土地进行重点治理，并且采取封山禁牧的一系列措施使水土流失和土地沙化问题得到缓解。

2 数据来源及研究方法

2.1 数据来源

1989—2019年逐日径流量、输沙量取自苦水河流域水文控制站——郭家桥水文站，降雨数据来自吴忠气象站。郭家桥水文站控制面积5 216 km2，年平均径流量0.933 8亿m3。土地利用遥感数据从PIE-engine遥感与地理信息云服务平台获取，产品为中国30 m年度覆盖数据(annual China Land Cover Dataset, CLCD)，该数据基于5463个独立参考样本，产品总体精度为79.31%。

2.2 研究方法

2.2.1 Mann-Kendall趋势分析 利用Mann-Kendall(下文简称M-K)趋势检验的方法对郭家桥水文站1989—2019年的径流、输沙数据做趋势性分析。M-K检验法对长时间序列的数据变化趋势具有良好的识别能力，在检验气候和水文要素的变化趋势上应用广泛。该方法不要求数据遵循特定的分布，同时个别异常值的干扰对结果影响不大，因此，可以相对客观地反映出长时间水文时间序列变化情况[23]。也可以进一步用此方法检验径流泥沙量时间序列的突变情况。具体计算方法参考文献[24]。

2.2.2 集中度分析 将集中度指数[25]应用于流域内径流和输沙时间序列，来分析年内径流、输沙的集中情况，集中度指数可以用来表达日降雨对年内总降水的贡献程度。计算步骤：首先将年径流量标准化，之后按照等级划分(0～0.000 9，0.001～0.009，0.01～0.09，0.1～0.9)，然后在逐个径流等级下计算出发生天数和径流量之和，并计算各个径流等级内的累积径流量百分比(Y，%)和累积径流天数的百分比(X，%)，用同样的方法计算输沙量的年内集中度CI。其中Y，X符合以下的洛伦兹曲线分布:

Y=aXebX

可以通过最小二乘法进行计算a和b的值，

两个常数被确定后，洛伦兹曲线在0到100之间的定积分就是曲线下的面积S可以通过以下公式计算得到：

最后计算面积S占总面积的百分比得到计算年份的集中度CI:

通过上述步骤得到的集中度介于0～1，靠近1则说明代表年径流、输沙量集中分布于一段时间内，反之，靠近0则代表年径流和输沙量分布较为分散[9]。

2.2.3 水沙关系曲线及其参数 水沙关系曲线常用来描述径流和泥沙之间的关系，是基于水文站实测径流和泥沙数据的一种经验回归关系[26]，在各种回归关系中，幂函数关系最为常见[27]，流量(Q)满足在一定参数下对悬移质泥沙输沙率(SSC)的函数。表达形式为:

SSC=aQb

式中:a为系数,b为指数。a和b的值表示沙源供应情况和水流影响下的泥沙输移的变化情况[28]。有学者对SRC函数关系中的参数a,b赋予了一定的实际意义并被广泛应用于水沙关系的研究中。系数a代表了流域内的侵蚀严重程度,反映了流域泥沙的可侵蚀性,系数a越大表明越易被侵蚀,同时,系数a的变化容易受到外界环境和人为活动干扰的影响;指数b与河流本身的性质有关,受河床形态、水流速度、水沙配比的影响,代表了河流自身的侵蚀能力,其值越大表明随流量的增加含沙量或输沙率增加的越快[29-30]。

2.2.4 贡献度分析 流域水沙情势变化主要受到气候变化和人类活动的影响[31]，为分析这两种因素对苦水河流域的影响贡献度，采用双累计曲线法定量计算。双累计曲线法常用于水文气象要素的趋势变化分析，用于判断人为扰动是否造成了流域要素的趋势性变化[32]。将逐年累计降雨量和径流量分别与逐年累计降水在径流、输沙突变前后建立线性回归直线，如果二者斜率差别较大，则说明人类活动对径流输沙变化的贡献度较高。并通过计算土地利用转移矩阵分析近30年间的土地利用变化情况，进一步判断人类活动的主要作用形式。

3 结果与分析

3.1 水沙序列的时间变化

3.1.1 年际尺度特征 苦水河流域1989—2019年均降雨198.10 mm，年均径流量14.16×107m3，年输沙量5.83×106t，降雨少，输沙量大。为探究其在年际尺度上变化趋势。对苦水河郭家桥水文站1989—2019年径流、泥沙数据进行时间序列分析，见图1，径流和泥沙的变化趋势相近，径流量大的年份往往伴随着高输沙量，但在2012—2019年径流量上升的过程中，输沙量并未显著提高，说明流域“大水大沙”的情况有所改善。进一步，对年径流量和年输沙量数据进行线性拟合分析，结果显示在近30年间，苦水河流域年径流量和输沙量均有减少的趋势。

进一步对苦水河近30年水沙数据做 Mann-Kendall趋势检验。研究时段内，苦水河径流量(p<0.1)和输沙量(p<0.05)锐减趋势达到显著水平。为判断年径流和输沙数据变化拐点，对径流输沙数据进行Mann-Kendall突变检验。结果显示，年径流时间序列Z值为-1.67，年输沙时间序列Z值为-3.37，均通过了显著性检验(p<0.05)，年径流量在2002年发生减少突变，而年输沙量减少突变在2004年发生。以突变点为时间划分节点，节点之前的时间段为基准时期，节点之后的时间段为变化时期，进一步结合径流输沙量的突变点将郭家桥水文站水沙序列做时间段划分，见表1。郭家桥水文站变化期较基准期年径流及年输沙量分别减少10.16%和87.26%。此结果表明，在1989—2019年，年径流量并没有大幅减少，与M-K检验结果减少趋势不显著相一致，但年输沙量却有明显的减少趋势。该结果为进一步判断气候变化和人类活动对苦水河流域水沙情势变化的影响提供了研究基础。

图1 1989-2019年径流输沙年际变化

表1 苦水河流域径流泥沙时间序列分析

3.1.2 月际尺度特征 对1989—2019年每个月的径流量和输沙量数据用做统计，见图2，利用箱线图分析径流和输沙的年内分配情况。结果显示，近30 a总径流量为4.45×109m3，总输沙量为1.93×108t，径流量和输沙量的年内分配比较集中，每年的5—8月为丰水期，丰水期径流量3.04×109m3，占总量的68.36%，同时这一期间也伴随着大量的泥沙输移，丰水期4个月内的输沙量为1.75×108t，占到总量的90.49%。每年的12月至3月为枯水期，输沙量为0.18×108t，占比为9.51%，输沙量较低。

3.1.3 日际尺度特征 为分析年径流量和年输沙量在日尺度上的分配情况，引入集中度指数(concentration index)。集中度指数最早应用于经济学领域，用来表征贫富差距，近些年来，有研究人员将这一指标应用到年内降雨在日尺度上的集中度上。本文拟采用这一概念来分析年径流和泥沙的分布情况。首先，用洛伦兹曲线拟合1989年郭家桥水文站的累计径流百分比和累计天数百分比数据，并验其拟合优度，发现决定系数R2为0.927 8，说明拟合效果较好。统计后续年份，发现决定系数均在0.9以上。

图2 径流、输沙的年内分配情况

1989—2019年径流和输沙集中度见图3，年径流集中度均小于年输沙集中度，年径流集中度均值为0.53，年输沙集中度均值为0.911。径流集中度指数的变化趋势采用M-K检验法进行研究，发现径流和输沙集中度指数的Z值分别是-4.22和-2.59，p值均小于0.01，均有显著下降趋势(p<0.01)，说明径流和输沙在年内的分布均有趋向于均匀的趋势，但输沙集中度仍然维持在较高的水平。

图3 年径流、输沙集中度分布

观察数据发现，近30年输沙量CI值的最大值发生于1997年，CI值达到0.99，并且在1999年、2001年CI值也达到了0.98，表明在这些年份苦水河流域在某一时间段输沙量比较集中。在1997年(图4)，7月30-7月31日仅两天的输沙量就达到年内总输沙量的74.08%。尽管苦水河流域近30 a输沙量大幅减少并且集中度指数有显著的(p<0.01)下降趋势，但其值仍旧较高，这表明苦水河流域内极端产沙事件频发，需要对流域洪水事件提高重视，防范未来洪水下的产沙事件。

为进一步分析降雨对输沙量的影响，在集中度不同的年份下选取降雨量大于10 mm、大于20 mm的降雨频次和日输沙量选取大于前5%(21 772.8 t)的输沙频次，旨在分析强降雨下的输沙情况。由表3看出，高集中度的年份在降雨发生频率上均高于低集中度年份。在强降雨(>10 mm、>20 mm)事件上则更加显著，例如在输沙集中度较高的2010年发生6次大于10 mm降雨事件，而在输沙集中度较低2011年则没有发生大于10 mm降雨事件；在输沙集中度较高的2002年发生4次大于10 mm降雨事件，而在输沙集中度较低2005年则仅发生1次大于10 mm降雨事件。在强降水场次多的年份其输沙量较大的天数也增多，表明强降雨事件对极端产沙事件具有显著的影响。

图4 1997年逐日输沙量分布

表2 极端产沙年份的强降雨频数

3.2 水沙关系的时间变化

3.2.1 水沙关系曲线及参数 拟合1989—2019年径流量和输沙量(图5)，符合幂函数分布且其拟合效果较好，决定系数R2为0.67。在水沙关系曲线的参数中，参数a值为5.48×10-6，参数b值为4.24。

图5 年径流量-年输沙量的水沙关系

为进一步探究参数a和b在近30年的变化趋势和规律，按照公式(6)将其逐年拟合计算参数a和b。计算结果见图6，a总体均值为9.07×10-5，在1989—2003年、2004—2019时间段内的均值分别为1.17×10-4，6.63×10-5；b总体均值为2.65，在1989—2003年、2004—2019时间段内的均值分别为2.61，2.70。然后对a和b的时间序列做M-K检验，序列a对应的Z值为-3.355 9，p值为0.000 8(<0.01)，说明近30年来，a值有显著下降的趋势，表明说明流域内侵蚀强度降低，泥沙可蚀性降低。由于系数a受到人类活动的影响较大，表明近30年间大规模的生态修复工程和退耕还林工程(2000年)发挥了较大的作用，人类活动对流域产水产沙的影响加大;序列b对应的Z值为1.02，p值为0.307 8(>0.1)，b值呈上升趋势但并不显著，说明河流本身的能量增加并不显著、河流的输沙能力和输沙特性变化不大。

图6 1989-2019年参数a和b的变化趋势

3.2.2 SSC-Q环路曲线 悬浮泥沙运移是一个与降雨特征、径流流量和泥沙可用性相关的复杂过程，受到水文和地貌过程时空异质性的显著影响[33]。悬浮泥沙浓度和流量(SSC-Q)滞后类形可用于解释悬浮泥沙动力学。滞后环的形状包括线性、顺时针形、逆时针形、正8字形、逆8字形、线性环路和复杂形，同时有助于确定排水系统内单个事件的沉积物可用性和沉积物源的近似空间分布[34-35]。

分析苦水河流域23场典型洪水过程，发现以下5种SSC-Q滞后曲线类型：正“8”字形、逆“8”字形、顺时针形、逆时针形和线形环路(图7)。逆时针滞回环路表示涨水期含沙量大于落水期含沙量，含沙量早于径流量达到峰值，此时泥沙物源比较丰富且主要来源于流域上游；相反的，顺时针滞回环路通常表示涨水期含沙量小于落水期含沙量，径流量早于含沙量达到峰值，说明泥沙补给不充分，泥沙主要来自于水文站附近[36]；逆“8”字形表示洪水过程中，滞回环路在径流量较低时表现为逆时针，同时在径流量较大时表现为顺时针，代表了短而急促的泥沙；正“8”字形与逆“8”字相反，代表了平缓且长久的泥沙；线形环路代表含沙量和径流量的输移时间和变化比例同步。

图7 5种 SSC-Q环路类型及其水沙趋势

统计不同时段26场典型洪水事件中5种环路类型的发生频率，结果见表3。近30年，顺时针滞回类型大于逆时针滞回类型。分阶段来看，2002年之后顺时针滞回类型的发生频率进一步增加。这表明苦水河流域的河流泥沙补给程度降低。同时，在2002年之前的时间段出现一次逆“8”形环路曲线，而在2002年之后的时间段出现正“8”形环路曲线，说明河流由短而急促的输沙模式渐渐转向平缓长久的输沙模式，进一步验证了近30年来苦水河输沙能力减弱，泥沙传播速率降低的特点。

表3 不同时间段水沙滞回类型频率

3.2.3 水沙变化的驱动因素 为进一步探究苦水河流域水沙减少的驱动因素，结合年径流量与年输沙量时间序列变化突变点，采用双累积曲线法分别分析气候条件和人类活动对苦水河流域水沙变化的影响(表4)。

表4 降雨与人类活动对径流和泥沙减少的贡献率

由图8看出，在2002年前后，拟合累计泥沙曲线的直线斜率明显变小，由0.061变化至0.006。而在2004年前后，拟合累计径流曲线的直线斜率并未发生明显的变化，由0.085减小至0.059。

图8 人类活动对径流和输沙变化的影响

结果表明：人类活动是苦水河的水沙减少的主要因素。径流变化百分比为23.49%，其中气候变化贡献率为26.35%，人类活动贡献率达73.65%；输沙变化百分比为87.26%，其中气候变化贡献率为4.33%，人类活动贡献率达95.67%。输沙变化相较于径流变化受人类活动的影响较大。

为进一步分析人类活动对水沙变化的作用机制，基于landsat遥感数据分析苦水河流域近30年的土地利用变化情况。通过分析1990年、2019年苦水河流域土地利用类型可知，草地是流域内优势土地利用类型，占流域内面积80%以上，其次为农田，占比为0.14%。采用叠加分析的方法分析近30年苦水河流域土地利用类型的变化情况，从空间上看，发生土地利用类型变化主要聚集在西部、北部、东南部，表明苦水河上、下游均发生了变化，且从变化面积和强度上看，下游(西北部)大于上游(南部)(图9)。

图9 1990-2019年土地利用类型的变化情况

流域土地利用类型转移矩阵表5表明：建设用地、林地、农田和水域面积增加，相应地，草地和荒地面积减少。1990—2019年，在总共407.80 km2农田中，有229.97 km2转化为草地，0.29 km2转化为林地，转化区域集中在流域上游；同时有519.45 km2由草地转化为农田，转化区域集中在流域下游；上游有146.38 km2荒地转化为草地；建设用地和水域的增加集中在流域中部，位于暖泉湖周边。

综合上述分析，在发生土地利用类型变化的区域在空间上具有明显的集聚效应，土地利用类型的变化主要集中在南部山区和西北部沿黄河平原。用地类型的变化主要是草地、农田和荒地之间的转化。南部山区土层浅薄，土壤贫瘠且水资源匮乏，不适合大面积耕作，而退耕还草可以有效地防风固沙，减少水土流失。而在水土条件较好的流域下游，增加耕种面积可以维持耕地的平衡，促进经济社会的发展。这种土地利用情况的转变促使苦水河流域水沙格局发生变化，特别是南部山区大面积的耕地转变为草地使得苦水河流域上游来沙进一步减少，这与前述苦水河流域输沙量减少相吻合。

表5 1990-2019年苦水河流域土地利用类型转移矩阵 单位：km2

4 结 论

(1) 苦水河年输沙量显著减少，年径流量有减少趋势但并不显著。输沙量和径流量分别在2002年和2004年发生减少突变，突变点前后变化量分别为10.16%和87.26%。年径流和泥沙集中在5—8月，且从日尺度上看，输沙集中度很高，需要防范汛期洪水引发的极端产沙事件。

(2) 通过分析水沙关系曲线，参数a有显著下降的趋势，表明流域内侵蚀程度降低，参数b有不显著的上升趋势，说明河流自身的输沙特性变化不大。洪水场次下的水沙关系呈现出以下5种环路曲线类型：顺时针、逆时针、正“8”字、逆 “8”字、线性。顺时针和正“8”字滞回类型增加，表明河流泥沙补给程度降低，输沙能力减弱，并由短而急促的输沙模式渐渐转向平缓长久的输沙模式。

(3) 近30 a苦水河流域内径流量和输沙量的减少，人类活动都是首要驱动因素，贡献度分别达到73.65%和95.67%。南部山区大范围的农田转变为草地有效降低了苦水河流域上游来沙量，退耕还草措施对促进水土保持、协调水沙关系具有积极意义。