黄河源区水沙变化特征分析

董晓宁，田世民，张 丽

（1.华北水利水电大学，河南 郑州450046；2.黄河水利科学研究院 河南省黄河水生态环境工程技术研究中心，河南 郑州450003；3.黄河水利科学研究院 河南省黄河流域生态环境保护与修复重点实验室，河南 郑州450003）

受全球气候变暖的影响，水循环过程发生了剧烈变化。黄河源区是黄河流域的重要产流区和水源涵养区，多年（1956—2017年）平均年径流量约198.2亿m3，占黄河年径流量的34.1%，多年平均年输沙量为0.1亿t，水资源状况对整个黄河流域有极大的影响，且其输沙量对上游梯级水库的寿命有极大影响。黄河源区地处青藏高原腹地，地形复杂，海拔落差较大，且流域内还分布有冰川、冻土、湖泊、湿地等，生态系统多样且脆弱，人类活动痕迹较少，几乎处于自然状态，水文过程对气候变化敏感程度较高，气候变暖导致的水循环剧烈变化在这里得到放大，多种影响因素叠加导致水沙过程较为复杂。

由于黄河源区对西部地区及黄河流域的重要性，因此许多专家学者深入分析了其径流演变机制。流域内水沙过程在时空变化上主要受降水影响，且降水转化率偏低［1］，此观点得到较多学者的认同［2-4］。除此之外，气温变化引起的蒸发变化［5-6］、冰川消退［7］和冻土的冻融过程［8］也对驱动水沙变化有一定的贡献。地表改造和引水等活动也是水沙特征改变的一大驱动因素［9］。湖泊、湿地的调蓄作用对流域内部分地区的水沙过程也存在影响［10］。降水与水沙变化成正比，气温与水沙变化成反比［11］。有学者对黄河源区近几十年径流减少各影响因素作用占比进行分析，发现20世纪90年代到2005年左右径流减少主要受人类活动影响，而2008年以后径流减少则主要受气候因素影响［12］。也有学者预测随着冰川积雪的减少，融雪径流将会减少，部分河流径流量会出现由增转减的“拐点”［13］，流域内多年冻土也会持续退化为季节冻土［14］。总之，黄河源区气候、下垫面、人类活动等与水文过程之间存在着相当复杂的关系。

要探究黄河源区水沙变化规律，仅从气候变化和人类活动等外在因素分析黄河源区径流、泥沙变化特征还不够全面，空间上流域内各区域的水沙变化特征及其对整个流域的水沙变化贡献率也值得研究，蓝永超等［15］认为在对黄河源区径流变化进行研究时，不应忽略流域内不同区域径流变化的差异。2019年9月习近平总书记提出要推进黄河流域生态保护和高质量发展，并对三江源、祁连山、秦岭等重点区域生态保护建设提出要求，以三江源、祁连山、甘南黄河上游水源涵养区等为重点，推进实施一批重大生态保护修复和建设工程，提升水源涵养能力。笔者运用统计学分析方法分析黄河源区把口站唐乃亥站径流泥沙演变趋势及特征，并将黄河源区分为黄河沿以上、黄河沿—吉迈、吉迈—玛曲、玛曲—唐乃亥4个区间，从空间上分析径流泥沙在各区间不同时间段水沙变化规律及区间产流产沙对唐乃亥站水沙变化的贡献率，以期为黄河源区生态系统恢复和水源涵养能力提高的区域性规划提供参考。

1 资料来源和研究方法

1.1 资料来源

选取黄河源区黄河沿、吉迈、玛曲、唐乃亥4个水文站的径流和泥沙数据，这4个水文站具有较长的数据系列，且资料少有丢失，可以代表黄河源区流域内的水沙变化。由于历史原因，黄河沿站1968—1975年水文资料缺失，下游吉迈站与黄河沿站的水文资料具有较好的相关关系，因此使用吉迈站资料对黄河沿站缺失资料进行插补。径流、泥沙资料来自水利部黄河水利委员会刊印的黄河流域水文年鉴。

1.2 研究方法

1.2.1 趋势检验

趋势检验使用Mann-Kendall非参数趋势检验法，根据统计量Z值的大小来确定序列变化趋势以及趋势是否显著。Z值大于0，序列呈上升趋势；Z值小于0，序列呈下降趋势。当Z的绝对值大于1.28、1.64、2.32时，分别表示趋势通过了显著性水平为90%、95%、99%的检验。Sen’s斜率估计是Sen于1968年提出并发展的一种非参数检验法，估计n个样本中N对数据的趋势斜率，其值可以反映数据趋势的陡峭程度。

1.2.2 突变检验

综合累积距平法和Mann-Kendall突变检验确定突变点，再用t检验确定突变点是否合理。距平是离散数据与其均值的偏差，累积距平是距平值的叠加。累积距平值的起伏变化可以用来判断序列的长期变化趋势以及持续时间，起伏拐点处则可能为序列突变点。Mann-Kendall突变检验是根据时间序列数据构造秩序列，根据秩序列计算统计量UFk和UBk，给定显著性水平，统计量UFk和UBk在临界值内的交点即为序列突变点，若UFk或UBk的值大于0，则表明序列呈上升趋势，小于0则表明呈下降趋势。当它们超过临界直线时，表明上升或下降趋势显著。

1.2.3 周期检验

周期性检验使用Morlet小波分析方法，绘制小波系数实部等值线图和小波方差图，峰值对应的点即为周期，最大峰值所对应的周期即为主周期，以此确定水沙序列的变化周期。

1.2.4 趋势持续性

用重标极差法（R/S分析法）计算Hurst指数、预测径流泥沙在未来一段时间内的变化趋势。当H=0.5时，表示时间序列不存在记忆性；当H＞0.5时，时间序列具有状态持续性，未来趋势与过去相似；当H＜0.5时，时间序列具有反持久性，未来趋势与过去相反。Hurst指数分级见表1。

表1 Hurst指数分级

2 研究结果分析

2.1 水沙变化趋势

唐乃亥站水沙变化曲线及累积距平曲线见图1。1956—2017年黄河源区唐乃亥站多年平均流量为628.58 m3/s，多年平均输沙率为369.94 kg/s。流量呈波动下降趋势，下降速率为1.21 m3/（s·a），年均流量最大值和最小值出现在1989年和2002年，其值分别为1 035.57 m3/s和334.51 m3/s；输沙率变化同样呈波动下降趋势，下降速率为1.62 kg/（s·a），年均输沙率最大值与年均流量最大值出现年份相同，为1989年，其值为1 298.49 kg/s，而最小值则出现在2008年，其值为87.29 kg/s。累积距平曲线表明径流和泥沙变化过程均经历了1989年之前的上升阶段和之后的下降阶段。

图1 唐乃亥站水沙变化曲线及累积距平曲线

计算得黄河源区水沙变化相关系数为0.85，说明径流在很大程度上影响泥沙过程。用变差系数CV来表示径流、泥沙的变异程度，泥沙过程的变差系数CV和极值比均大于径流过程（见表2）。使用Mann-Kendall趋势检验法对黄河源区径流、泥沙过程进行分析，得到统计量Z分别为-1.02和-0.98，均未通过置信度为90%的显著性检验，说明黄河源区径流泥沙虽一直呈下降趋势，但并未达到显著性水平。

表2 径流泥沙变异程度

2.2 水沙突变和周期性检验

2.2.1 突变检验

使用Mann-Kendall突变检验法（以下简称MK突变检验）对黄河源区1956—2017年径流、泥沙过程进行分析，如图2所示。结合前文所述累积距平曲线所展现出的特征，确定几个较为可能的突变点，再对可能突变点进行t检验，根据统计量t值的大小决定突变的年份。此处给定显著性水平α=0.05，MK突变检验上下临界值为±1.96。

图2 黄河源区水沙MK突变检验

由图2可以看出，径流的拐点在1989年、2008年，泥沙的拐点在1989年。MK突变检验中径流过程UFk与UBk曲线的交点在1993年、2008年、2012年，泥沙过程UFk与UBk曲线的交点在2006年、2013年。将两种方法所得的可能突变年份设为基准年，对其前后两个序列做t检验，得到统计量t并判断其显著性，检验结果见表3。

表3 突变点显著性t检验结果

由于本文研究的时间尺度基于年代际变化，因此上述可能突变点距离时间系列终止年份不足10 a的不再考虑，相邻两个突变年份相差不超过10 a的，取统计量t值高的突变点为突变年份。总结表3检验结果，1989年、1993年均检测到径流突变，但两个突变年份仅相差5 a，取突变更显著的1989年为径流突变年份；1989年、2013年均检测到泥沙突变，但2013年距资料系列终止年份仅有4 a，取1989年为泥沙突变年份。径流、泥沙过程的突变年份均在1989年，区别在于径流过程相对泥沙过程突变结果更为显著［16］。

2.2.2 周期检验

对1956—2017年唐乃亥站径流、泥沙过程进行小波分析，分别绘制小波系数实部等值线和小波方差，如图3所示。径流、泥沙过程存在周期性，从图3可以看出，径流和泥沙在周期性变化上存在相似性，径流过程存在7、12、25 a的周期和45 a的准周期，泥沙过程存在5、12、25 a的周期和34、44 a的准周期。径流过程的7 a周期主要发生在1966—1980年，12 a周期发生在1956—2000年，25 a周期发生在1980—2010年。泥沙过程的5 a周期主要发生在1980—1995年，12 a周期发生在1956—2000年，25 a周期主要发生在1980年以后，从图3可以看出泥沙过程在1965—1980年也存在7 a的周期。结合小波系数实部等值线图的周期密集性和小波方差图的峰值大小来看，在5～7 a和12 a周期，泥沙过程震荡强度大于径流过程的，尤其在12 a周期的震荡强度，而在25 a周期尺度上泥沙过程的震荡强度则小于径流过程的。

图3 水沙小波分析

2.3 水沙变化的区间趋势

黄河源区水文站点稀少，有50 a以上资料的站点仅有黄河沿、吉迈、玛曲、唐乃亥4个，因此本文选用这4个站点1959—2017年的水沙资料来研究黄河源区的水沙变化特征，并将黄河源区分为黄河沿以上、黄河沿—吉迈、吉迈—玛曲、玛曲—唐乃亥4个区间，按区间上下游顺序记为区间一、二、三、四，将相邻上、下游两站实测径流泥沙数据之差作为区间产流产沙量。由表4可以看出，黄河源区的最大产流区为区间三，多年平均产流量占总产流量的50.80%；最大的产沙区为区间四，多年平均产沙量占总产沙量的62.53%。使用MK趋势检验法对各区间水沙趋势进行分析，由表5得出，4个区间除区间二径流过程表现出不显著上升趋势外，其余3个区间都表现出下降趋势，且主产流区区间三表现出置信度为99%的显著下降趋势；4个区间的泥沙量都表现出下降趋势，区间一和区间四均通过置信度为90%的显著水平检验，而区间二则通过置信度为99%的显著水平检验。

表4 黄河源区水沙区间占比

表5 水沙MK趋势检验

根据上述突变点选取标准，对表6可能突变点显著性检验结果总结后得出，区间一径流过程突变年份为1990年、2006年，泥沙过程突变年份为1994年；区间二径流过程突变年份为1986年、2004年，泥沙过程突变年份为1985年；区间三径流过程突变年份为1989年，泥沙过程突变年份为1994年；区间四径流过程突变年份为1989年、2007年，泥沙过程突变年份为1989年。

2.4 不同时期的水沙区间变化

分析各区间水沙突变年份可以发现，黄河源区各区间水沙突变年份和变化趋势具有很大的相似性，根据突变点年份再综合水沙贡献率区间占比，选择1989年和2007年作为此次研究的“拐点”，由此将黄河源区水沙变化趋势分为1959—1989年、1990—2007年和2008—2017年3个时段来研究各区间的水沙变化，按时间顺序分别将3个时期定义为T1、T2、T3。各区间年均产流产沙量和不同时期年均水沙变化率分别见表7和表8。

表6 各区间水沙突变点显著性t检验结果

从产流时期来看，所有区间都经历了T2时期的径流减少和T3时期的径流回升，上游两区间T3时期的年均产流量超过了T1时期，而下游两区间T3时期径流量虽有所回升，但并未达到T1时期的产流量。从T1—T3时期的年均产流量和变化率来看，区间三既是产流量减少最多又是变化率最大的区间。

表8 各区间不同时期年均水沙变化率 %

黄河源区产沙区域主要集中在区间三和区间四，其中约60%的泥沙来自区间四、约30%的泥沙来自区间三，而区间一和区间二的产沙量约占总产沙量的10%。与径流变化相似，各区间产沙量也经历了T2时期的减少，但在T3时期只有上游3个区间产沙量有所回升，区间四的产沙量及其占比均在持续下降。从年均产沙量来看，产沙量减少最多的是区间四，但从泥沙变化率来看，变化最大的是区间一。

2.5 水沙变化的区间贡献率

将不同时期各区间对促进黄河源区水沙变化的作用占比定义为水沙变化区间贡献率。黄河源区不同时期的年均水沙产量变化情况和区间贡献率见表9。

根据表9，T1到T2时期黄河源区年均产流量减少48.14亿m3，区间三和区间四对径流减少的贡献率较大，分别是45.60%和35.56%，区间一和区间二对径流减少的贡献率比较接近，均不到10%；T2到T3时期黄河源区年均产流量增加29.09亿m3，区间四对径流增加的贡献率为55.92%，区间二和区间一对径流增加的贡献率分别为25.03%和17.98%，区间三对径流增加的贡献率则仅有1.07%；T1到T3时期黄河源区年均产流量减少19.05亿m3，区间三对径流减少的贡献率达113.60%，区间一和区间二对径流减少反而起到牵制作用。

表7 各区间年均产流产沙量

表9 不同时期水沙变化情况和区间贡献率

T1到T2时期黄河源区年均输沙量减少492.71万t，区间四对输沙量减少的贡献率最大，超过其他3个区间对输沙量贡献率的总和；T2到T3时期黄河源区年均输沙量减少2.28万t，区间四对输沙量减少的贡献率最大，达2 872.81%，而其余3个区间则均对输沙量减少起牵制作用；T1到T3时期黄河源区年均输沙量减少494.99万t，区间四对输沙量减少的贡献率为81.81%，其余3个区间对输沙量减少的贡献率总和不到20%。

2.6 区间输沙模数

通常用输沙模数来表示流域内水土流失程度，从表10可以看出，从上游至下游水土流失的严重程度逐渐加重，区间四水土流失最为严重。但可以看到区间四在整个时期内输沙模数在减小，此结果与表7年均产沙量的变化一致。

表10 各区间不同时期输沙模数 t/（km2·a）

2.7 水沙趋势持续性

用R/S分析法估算Hurst指数以确定未来一段时间水沙变化趋势，估算的H值见表11。

表11 黄河源区各区间H值

由表11可知，各区间水沙过程的H值均小于0.5，表示黄河源区水沙变化存在反持续性，且泥沙的反持续性稍强于径流的，对比上文所得MK趋势检验结果，未来一段时间黄河源区各区间水沙将会呈现出上升的趋势，这一结果与蓝永超等［17］所得未来径流是先扬后抑的结果相一致。再将表11与Hurst指数分级表对照，得到黄河源区各区间径流的反持续性的强度均为较弱，而泥沙则是区间一和区间四的反持续性较弱，其余区间的反持续性较强。

3 讨 论

黄河源区水沙变化的影响因素主要可以归结为气候变化和人类活动两种，气候变化主要通过降水和气温影响水沙变化，人类活动带来的影响主要是引水用水和对土地覆被的改造等。

径流对降水变化最为敏感，从降水空间分布来看，黄河源区降水最丰沛的地区在东南部［18］，而区间三主要处于黄河源区的中部和东南部，降水空间上可以解释区间三是黄河源区的主产流区间，且区间三同时有来自黑河和白河两大支流的径流补充，黑河多年平均产流量占区间三产流的16.1%，白河多年平均产流量占区间三产流量的34.3%。

黄河源区T1到T2时期径流减少的原因已被多位学者探讨且得出结论，径流减少是降水减少、气温升高和人类活动加剧共同作用的结果［19］，其中降水是径流变化的主导因素［20］，尤其在T2时期降水强度和空间变化对径流减少起到很大作用［21-22］。

观察到区间二产流量的MK趋势检验结果与其余区间不同，呈微弱上升趋势，且T1到T2时期年均产流量的变化率是最小的。对此也有学者给出相应的研究结果，黄河源区的冰川冻土主要集中在吉迈以上区域，阿尼玛卿山冰川1966—2000年面积萎缩21.7 km2，冰川储量亏损2.66 km3，冰川水资源损失23.9亿m3［7，23］。20世纪60年代到90年代约有3.1万km2的多年冻土转化为季节冻土［24］，而多年冻土退化对径流的影响在多年冻土覆盖率大于40%的流域较为显著［8］。因此，区间二产流量在T1到T2时期变化率较小是因为受到冰川冻土融水的补给。

区间四在T1到T2时期和T2到T3时期产沙量均在减少，虽然T1到T2时期气温升高导致的植被覆盖度降低和人类活动对地表的破坏会导致产沙量增加［25-26］，但对输沙量变化影响最大、最显著的因素是径流，因此导致T1到T2时期输沙量减少的主要原因是径流量减少。而在T2到T3时期区间四产流量回升，而输沙量依然减少，是因为2005年国家三江源生态修复计划的实施，2005—2012年三江源区草地覆盖率明显提高，尤其在兴海县北部和玛多县草地覆盖率提升10%以上［27］，且有研究表明随着地表植被覆盖度的提高，流域产沙量会减少［28］。

4 结 论

（1）1956—2017年黄河源区水沙变化均呈下降趋势，但并未达到显著性水平，径流下降速率为1.21 m3/（s·a），泥沙下降速率为1.62 kg/（s·a）。

（2）使用MK、累积距平、滑动t检验结合的方法，检验出水沙过程的突变点均为1989年。使用小波分析法对水沙过程进行周期性检验，水沙周期变化有很强的相似性，均存在5～7 a、12 a和25 a的周期变化，区别在于泥沙过程12 a的周期震荡强度远大于径流过程的。

（3）对水沙过程进行区间分析，60多a来黄河源区各区间产流量均经历了T1到T2时期的减少和T2到T3时期的回升；区间一、区间二和区间三输沙量同样经历了T1到T2时期的减少和T2到T3时期的回升，而区间四在3个时期内产沙量在持续减少。

（4）T1到T2时期对径流量减少贡献率最大的是区间三，对输沙量减少贡献率最大的是区间四；T2到T3时期对径流量增加贡献率最大的是区间四，对输沙量减少贡献率最大的是区间四。

（5）结合R/S法和MK趋势检验，径流量、输沙量在未来一段时间内将会呈上升趋势。

【责任编辑 张 帅】

［20］ ZHENG H，YUAN J，CHEN L.Short-Term Load Forecasting Using EMD-LSTM Neural Networks with a Xgboost Algorithm for Feature Importance Evaluation［J］.Energies，2017，10（8）:1168-1172.