祁连山北坡分布着冰川、雪山、森林、草原、沼泽、湿地、戈壁等多种地貌类型，生态地位极其特殊^[1]，是河西走廊多条内陆河的发源地和下游绿洲社会经济永续发展的基础性保障，承担着水源涵养、生物多样性保护等多种生态功能，是西北乃至全国重要的生态屏障区^[2]。受地形、地貌、气候等多种因素的影响，境内土壤侵蚀因子复杂多样，是水力、风力、冻融、重力等多种侵蚀要素的复合体^[3-4]。

长期以来，祁连山区径流泥沙、水文生态、侵蚀机理、侵蚀速率等方面备受广大学者关注，并在不同时空尺度上取得大量研究成果。在研究内容上，揭示了祁连山区径流变化对气候变化的响应^[5-6]、不同植被条件下的水源涵养、蓄水保土等性能^[7-8]。在研究范围上，对河西三大流域均有涉及，例如，汤萃文等^[9]以石羊河上游山区为研究对象，定量分析了各种环境因子与土壤侵蚀之间的关系；郭斌等^[10]对黑河莺落峡控制站自20世纪60年代至21世纪初共50 a时段黑河干流水电站建设前后泥沙变化进行了系统分析，认为水土流失最严重的区域为浅山区；严宇红等^[11]以长序列资料为基础，系统分析疏勒河流域上下游泥沙分布特征，并建立了水沙关系模型。研究方法与手段上，包括水文统计法、地面监测法、原子示踪法、遥感监测法等，在时间尺度上，地质尺度、千年尺度、现代尺度的研究成果均有报道^[12-13]，其中大多以各流域长序列水文实测数据统计分析为主^[14-17]；但长尺度(地质尺度、千年尺度等)的地质构造、气候变化等人类难以有效控制的环境影响因素，对当下环境治理与保护决策难以起到支撑作用。近20 a来，祁连山生态环境经历由资源无序开发向环境有序整改、生态逐渐修复的过程^[18-19]，尤其自2017年开始，甘肃省针对中央就甘肃祁连山国家级自然保护区生态环境问题通报内容开展全方位整改，而目前针对祁连山生态环境问题整改前后流域空间尺度上水沙变化及土壤侵蚀强度的相关研究少见报道。因此，研究近20 a祁连山生态恢复过程中区域水沙动态变化及土壤侵蚀时空特征，可为客观研判生态修复效果、科学防治水土流失、有效提高生态承载能力提供一定的技术依据和参考。

祁连山北坡由东向西依次属于石羊河流域、黑河流域、疏勒河流域3大内陆河流域出山口上游，地势由西南向东北倾斜，南高北低、西高东低。境内高寒山地、高原湖泊、河谷盆地、戈壁平原多种地貌并存，以土石山地为主要地貌特征。高寒山地气候条件属于干旱与半湿润气候，而山前洪积扇为大陆性荒漠气候。寒冷、干旱是境内主要气候特征，降水量自东向西递减，由北向南递增，多年平均降水量在108~335 mm之间。气候与植被的垂直地带性差异明显，由高到低依次为高寒草甸、高寒灌丛、寒温性针叶林、温性灌丛、温性草原和荒漠草原植被。在高寒草甸植被带以上，基本属于高山裸岩或被冰雪覆盖，人迹罕至；在海拔3 000 m左右，部分区域植被覆盖度60%~90%^[20]；在浅山戈壁地带，植被稀疏，覆盖度10%~30%。土壤分布由高到低有高山漠土、高山寒漠土、高山草原土、亚高山草原土、亚高山草甸土、草甸土、盐土、灰漠土、风沙土等^{[12, 21]}。

在全国水土保持区划中，祁连山北坡属于青藏高原区—柴达木盆地及昆仑山北麓高原区—祁连山山地水源涵养保土区，容许土壤侵蚀量为500 t/(km²·a)^[22]。境内气候、地形条件复杂，土壤侵蚀的外营力复杂多元化，从高到低、由南向北，基本随地形和气候的变化呈以冻融侵蚀、重力侵蚀、水力侵蚀、风力侵蚀为主的规律变化。水力侵蚀形成的大小冲沟在山前倾斜平原上呈树枝状展布，部分沟道经过山前洪积扇时，有泥石流发育。在中高山地有水力侵蚀和冻融侵蚀，但范围极小，且大多处于无人居住区，基本属于不可治理的范围。在浅山区为水力、风力混合侵蚀区，分布有小片绿洲，是人口聚集、农牧业生产及生产建设活动比较集中的区域，也是水土流失治理和预防监督的重点区域。

在涉及河西走廊武威、张掖、酒泉3个市级行政区的祁连山浅山区，选择有完整长系列监测数据的水文测站作为分析站，从东到西依次选择石羊河流域的杂木寺和九条岭水文站、黑河流域的莺落峡水文站、疏勒河流域的昌马堡和党城湾水文站(表 1和图 1)，以各水文站2002—2021年实测径流、泥沙、降水资料为基本数据源；推移质数据以各站多年平均输沙量为基数，按14%的推悬比系数计算^[15]。在分析各测站近20 a水文泥沙变化的基础上，加入推移质输沙量，分析侵蚀强度时空变化。采用水文统计法、距平分析法、相关分析法、方差分析法等研究方法，利用SPSS25对数据进行单因素方差分析和双变量相关分析，利用Origin 2018进行数据绘图，利用AutoCAD 2007软件绘制研究区水文站分布图。为削减水文泥沙年际间的突变、有助于判别分析水文泥沙变化过程，按距平分析方法，将各时段距平值与各站20 a平均值相比较，求得各时段的距平分数(%)，分析径流泥沙变化特征，计算式如下：

表 1(Tab. 1)

表 1 各代表站点基本信息表 Tab. 1 Basic information of each representative station 站名 Station name 位置 Location 记录起始年份 Record start year 控制面积 Controlled area/km2 所属流域 Watershed 所属水系 Water system 多年平均径流量 Average runoff over the years/108 m3 多年平均输沙量 Average sediment load over the years/(104 t·a-1) 推移质量 Amount of push mass/(103 t·a-1) 杂木寺站 Zamusi station E 102°34′, N 37°42′ 1956 851 杂木河 Zamu river 石羊河 Shiyang river 2.365 13.5 19 九条岭站 Jiutiaoling station E 102°03′, N 37°52′ 1957 1 077 西营河 Xiying river 石羊河 Shiyang river 3.195 13.1 18 莺落峡站 Yingluoxia station E 100°11′, N 38°14′ 1956 10 009 黑河 Black river 黑河 Black river 16.414 202.2 283 昌马堡站 Changmabao station E 96°51′, N 39°39′ 1956 10 961 昌马河 Changma river 疏勒河 Shule river 10.271 340.0 476 党城湾站 Dangchengwan station E 94°53′, N 39°39′ 1972 14 325 党河 Dang river 疏勒河 Shule river 3.381 73.8 103

表 1 各代表站点基本信息表 Tab. 1 Basic information of each representative station

图 1(Fig. 1)

图 1 研究区水文站分布图 Fig. 1 Distribution of hydrological stations in the study area

式中，P_i为i时段距平分数；A_i为i时段年实测均值；a为某水文站近20 a径流或泥沙平均值。

根据各代表站近20 a来径流泥沙过程线(图 2)看，各站年际间变化起伏较大，尤其输沙量变化过程极不平稳，且逐年实测数据变化规律不明显，根据20 a内各站径流泥沙值变异系数分析结果，输沙量变异系数远高于径流量变异系数。径流量变异系数排序为：莺落峡站＜党城湾站＜九条岭站＜杂木寺站＜昌马堡站；输沙量变异系数排序为：昌马堡站＜九条岭站＜党城湾站＜杂木寺站＜莺落峡站。莺落峡站在径流量最平稳的状态下，输沙量变率最大(图 2和3)。

图 2(Fig. 2)

图 2 各代表站径流泥沙过程线 Fig. 2 Runoff sediment process line of each representative station

图 3(Fig. 3)

图 3 径流量与输沙量距平图 Fig. 3 Anomaly map of runoff and sediment load

杂木寺站各时段径流量变化幅度均不大，但总体表现为后10 a较前10 a略有减少。最大增减幅度分别出现2002—2006年和2012—2016年，此2个时段间，距平分数相差11.34%，增减突变分别发生在2003年和2013年，大小极值分别为2.958亿和1.488亿m³，相差0.99倍，2014年后开始小幅度增加。泥沙变化总体趋势基本与径流量变化趋势一致，即前多后少，但变化振荡幅度大，其中2002—2006年输沙量大幅度增加，在径流量距平分数仅增加5.23%的情况下，输沙量距平分数剧增42.63%；输沙量减幅最大时段出现在2017—2021年，在径流量略有减少的情况下，输沙量距平分数减幅达到75.98%。总体分析，后10 a与前10 a相比，在径流量减少9.82%的情况下，输沙量减幅达到38.81%。

同属祁连山东段石羊河水系的九条岭站水沙变化则与杂木寺站不完全相同，后10 a较前10 a相比较，在径流量仅减少2.24%的情况下，输沙量增加79.35%，泥沙与径流出现逆向变化过程。径流量在2007—2016年的10 a内为减少时段，2002—2006年和2017—2021年为增加时段；输沙量在2002—2016年的3个时段均为减少时段，最小值出现在2013年，但到2017—2021年，输沙量剧增，突变发生在2018年和2019年，分别达到40.7万和51.4万t，分别高于平均值166.75%和237.05%，大小极值相差14倍。

黑河流域代表站莺落峡站各时段径流量变化表现为持续增加态势，后10 a较前10 a增加11.88%，同期输沙量增加63.96%。输沙量最小值为3.3万t，发生在2021年，最大值为403.7万t，发生年份为2016年，输沙量动态变化振荡幅度大，大小极值相差121倍。

昌马堡站径流量总体呈增加趋势，输沙量呈减少趋势，后10 a较前10 a相比，径流量增加7.45%，输沙量减少14.45%。其中，2002—2006年为输沙量最大时段，2012—2016年为输沙量最小时段。输沙量最大年份出现在2002年，达到1 362万t，最小值为2020年的121万t，大小极值相差10.31倍。该站多年平均输沙率为107.8 kg/s，所代表的昌马河流域为研究区域输沙量最高的流域。

党城湾站径流量前3个时段距平值均为负值，到2017—2021年大幅增高。对应的输沙量表现为先增后减的变化过程，其中2007—2011年内，在径流量无明显增加的情况下，输沙量距平值高达39.64%，突变值在2010年，达到328万t，是20 a年平均值的4.2倍。总体后10 a与前10 a相比，在年均径流量增加11.47%的情况下，输沙量减少40.16%。

从时空分布上分析前后10 a径流动态变化过程发现，祁连山东段石羊河水系的杂木寺站、九条岭站径流量不同程度减少，中段和西段黑河流域、疏勒河流域代表站的径流量均有不同程度增加。前后10 a输沙量对比发现，九条岭站和莺落峡站输沙量不同程度增加，其中九条岭站增幅最高；杂木寺站、昌马堡站、党城湾站输沙量不同程度减少，其中党城湾站减幅最大，昌马堡站减幅最小。

用近20 a输沙量平均值与多年平均值对比分析发现，党城湾站有小幅增加，昌马堡站有大幅增加，莺落峡站有大幅减少(图 4)。

图 4(Fig. 4)

图 4 各代表站输沙量变化特征 Fig. 4 Variation characteristics of sediment load of each representative station

一般情况下，土壤侵蚀强度的大小可直接表现在河流输沙量上，流域产沙与地表径流量的大小及径流形态特征有直接关联，但从各站径流、泥沙动态变化分析结果发现，径流变化与泥沙变化的差异性较大。为进一步探究径流与泥沙之间的关系，对5个代表站径流量和输沙量进行双变量相关分析(表 2)，可以看出，各站输沙量与径流量均呈正相关，其中，九条岭站和昌马堡站输沙量与径流量均呈极显著相关性(P < 0.01)，杂木寺站和党城湾站输沙量与径流量均呈显著相关性(P < 0.05)，莺落峡站输沙量与径流量正相关性不显著(P>0.05)。这说明输沙量随径流量呈正向变化，径流变化是输沙量变化的主要驱动因素之一，但两者的相关性在各代表站间表现出一定的分异性。

表 2(Tab. 2)

表 2 径流泥沙量相关性分析 Tab. 2 Correlation analysis of runoff and sediment load 站名Station name R P 杂木寺站Zamusi station 0.527* 0.017 九条岭站Jiutiaoling station 0.583** 0.007 莺落峡站Yingluoxia station 0.200 0.400 昌马堡站Changmabao station 0.620** 0.004 党城湾站Dangchengwan station 0.452* 0.045 注：*表示P < 0.05，**表示P < 0.01。Notes: * indicates P < 0.05, ** indicates P < 0.01.

表 2 径流泥沙量相关性分析 Tab. 2 Correlation analysis of runoff and sediment load

以各流域实测泥沙数据为基础，加入推移质后计算汇总各站不同时段侵蚀模数，分析其侵蚀强度(表 3)，可知，同一代表站在不同时段的侵蚀模数差异均不显著(P>0.05)，各站点侵蚀模数变化过程与输沙量变化过程一致。同一时段不同站点间侵蚀模数差异均显著(P < 0.05)，昌马堡站侵蚀模数显著高于其他4个代表站，且其他4个代表站间侵蚀模数差异均不显著(P>0.05)。各站近20 a平均水力侵蚀模数由大到小排序为：昌马堡站>九条岭站>杂木寺站>莺落峡站>党城湾站，均小于500 t/(km²·a)，侵蚀速率在0.34~0.05 mm/a之间，属于微度侵蚀，侵蚀强度在容许侵蚀量以内，但昌马堡站侵蚀量接近容许侵蚀量临界值，且年际变化不稳定，在20 a的统计资料中，有6 a超过容许侵蚀量，有4 a侵蚀模数在400~500 t/(km²·a)之间，接近容许侵蚀量。

表 3(Tab. 3)

表 3 不同时段水文站侵蚀模数 Tab. 3 Erosion moduli of hydrological stations in different period  t/(km2·a) 时段 Period 杂木寺站 Zamusi station 九条岭站 Jiutiaoling station 莺落峡站 Yingluoxia station 昌马堡站 Changmabao station 党城湾站 Dangchengwan station 2002—2006 196.60±66.30Ba 109.40±26.75Ba 109.40±42.34Ba 573.20±195.58Aa 67.20±16.62Ba 2004—2011 151.40±28.67Ba 126.80±20.30Ba 70.60±22.92Ba 436.20±43.36Aa 82.80±38.82Ba 2012—2016 178.60±99.76ABa 116.80±21.45Ba 191.00±64.51ABa 403.20±132.01Aa 54.60±14.38Ba 2017—2021 51.60±6.18Ca 280.00±70.39Ba 68.60±10.44Ca 473.00±124.91Aa 41.40±8.79Ca 平均Average 144.55 158.25 109.90 471.40 61.50 注：同一行不同大写字母表示同一时段下不同站点间的差异性(P < 0.05)；同一列不同小写字母表示同一站点不同时段下的差异性(P < 0.05)。Notes: Different capital letters in the same row indicate the difference between different sites in the same time period (P < 0.05); different lowercase letters in the same column indicate the difference between the same site and different time periods (P < 0.05).

表 3 不同时段水文站侵蚀模数 Tab. 3 Erosion moduli of hydrological stations in different period

径流变化是输沙量变化的主要驱动因素之一，但两者的相关性在各代表站间表现出一定的分异性。分异性产生的原因主要包括：1)内陆河流域径流补给来源有降水、冰雪融化、沿途地下水补给等^[23]，除强降雨引发地表侵蚀外，冰雪融水和地下水补给只对沿途河道产生不同程度的冲刷侵蚀，并不产生因雨滴动能击溅地表引发大面积地表侵蚀。这种径流补给方式与外河流域相比，具有明显不同的特点，使径流泥沙相关系数普遍较低；2)水电站工程的施工与建设改变了径流泥沙形态，对二者相关性的分异产生影响。例如，在分析的5个站点中，相关系数最低、分异值最大的莺落峡所在河道上游自2001年后，陆续开发运行引水式、混合式梯级水电站7座；九条岭站上游从2013年开始陆续建设西营河1~3级引水式水电站，但一直未能投运，至2018年在祁连山环保问题整改时被拆除退出。这些电站建设施工及拆除过程中会有一定量泥沙进入河道，同时电站在运营过程中改变河道水文情势，对河道径流及泥沙产生一定影响。

用最近20 a径流资料分析，祁连山北坡径流量变化趋势与以往使用长序列数据分析的结果基本一致，是以往变化趋势的延续，即东段径流量减少、中段和西段径流量增加，主要是受区域内气候变化影响所致^[24]。

强降雨是诱发输沙量剧增的主要动力。研究发现，月最大输沙量总是伴随月最大降水量出现，如莺落峡站2016年8月记录降水量57 mm，是研究时段内该站的最大值，对应的本月输沙率为994 kg/s，也是研究时段内的最大值；九条岭站2018年8月记录降水量147.4 mm，是研究时段内该站月降雨量最大值，对应的当月输沙率62.5 kg/s为同期的最大值。但年最大输沙量与年最大降水量并不相伴出现，如九条岭站年最大降水量发生出现在2014年，为432 mm，但年输沙率最大值16.3 kg/s却出现在2019年。由此分析，九条岭站和莺落峡站后10 a输沙量增加，主要是强降雨天气增多的原因。

减弱或消除人为干扰因素，可有效减少流域输沙量。自2017年以来，当地对祁连山生态环境问题进行全方位整改，包括治理超载过牧、保障河道生态流量、清退旅游设施、关停非法矿山、恢复矿山植被等措施，仅在祁连山西段党城湾站上游的甘肃盐池湾国家级自然保护区内关停整改探矿点40个、生产矿山5个，冻结或注销探矿权，清理矿渣等废弃物，全面实施生态恢复工程。祁连山3个代表站输沙量减少，应与生态环境问题整改效果的发挥有关，说明流域内植被增加、生态环境好转。

不同站点之间侵蚀模数的差异性主要与所在流域下垫面物质组成有关，尤其与植被覆盖度关系密切。侵蚀模数最大的昌马河流域表现为裸岩山地、戈壁平原、天然草地嵌套的复合地貌景观。境内植被稀疏，多为风蚀劣地。河槽两岸为坡积、洪积碎石土、风沙土等松散堆积物构成的昌马洪积扇，为研究区域最大的洪积扇，并且在洪积扇上多有泥石流沟道发育。在强降雨条件下，洪水容易携带堆积物进入河道，增加侵蚀源。侵蚀模数最小的党河流域上游为盐池湾自然保护区，有灌丛草甸、草原化草甸，草甸化草场、盐生草甸草场、沼泽化草甸草场等优质牧草地，尤其在平草湖一带水草丰美，植被覆盖度超过80%，对土壤侵蚀有极强的抑制作用。在党城湾镇周围的山间盆地倾斜平原上有绿洲农田分布，为当地的主要饲草基地，植被覆盖度高，侵蚀强度小。有研究显示，土壤侵蚀过程变化与所研究对象的时空尺度及研究方法有密切关系，本研究以近20 a短序列水文泥沙数据统计分析为主，获得的侵蚀强度结果与以往研究结果不尽一致。这可能与研究的时空尺度及研究方法不同有关^[25-28]。祁连山区侵蚀因子复杂多样，以河道输沙量为基础数据分析，各控制站的侵蚀量均在容许侵蚀范围内，但水文分析法没有考虑风力侵蚀因素，所计算的侵蚀模数中没有包含风蚀量，以后研究中尚需加强。针对研究区水力侵蚀的危害，根据具体的保护对象，开展如洪积扇泥石流沟道治理、倾斜平原侵蚀沟道治理、中小河流治理等工程及生态清洁型小流域建设。针对研究区风力侵蚀危害，应结合人居环境改善与美丽乡村建设目标，以主要城镇为中心，围绕山间盆地及倾斜平原的绿洲饲草基地，开展农田防护林建设、牧区水利建设等内容^[29-30]，巩固退牧还草成果，大范围保护与小范围治理相结合，达到建设小绿洲、保护大草原的目的。

1) 与前10 a相比，最近10 a祁连山东段石羊河水系杂木寺站和九条岭站径流量均出现不同程度的减少，减少幅度分别为9.82%和2.24%；而中西段径流量不同程度增加，莺落峡、昌马堡和党城湾站径流量的增幅分别为11.88%、7.45%和11.47%。

2) 与前10 a相比，近10 a杂木寺、昌马堡、党城湾站输沙量均有不同程度减少，减少幅度分别为38.81%、14.45%和40.16%，表明流域内植被覆盖度增加，生态环境呈现好转趋势。

3) 与前10 a相比，近10 a九条岭和莺落峡2站的输沙量有不同程度增加，增幅分别为79.35%和63.96%，主要与强降雨有关，也与前后时段水电站的盲目建设与拆除整改有关。

4) 各流域侵蚀强度均在容许侵蚀范围之内，但昌马河流域洪积扇发育，植被覆盖条件差，侵蚀量年际变化不稳定，在强降雨条件下易发生超过容许侵蚀量的水土流失，需重点关注，并有效治理。