青藏高原大通河流域径流变化归因分析

王淑芝， 温得平

（1.青海省水文水资源测报中心，青海 西宁 810001；2.青海省水旱灾害防御服务中心，青海 西宁 810001）

大通河地处青藏高原东北边缘，生态环境脆弱，是我国西北地区重要的生态过渡带，随着水能水资源开发等人类活动的加剧，河流水文过程、水生态环境受到严重影响。另外，受全球气候暖湿化影响，流域降水、径流等关键水循环要素也发生了显著变化，人类活动加剧的水循环不仅会影响到水资源生态环境，还将使下游地区的水资源安全面临诸多风险，因此开展归因分析，研究其水文变化机制对流域水资源科学管理和保护利用具有重要意义。当前，统计分析法、水文模型法和弹性系数法等是定量识别径流变化归因的主要方法［1-2］。吴恒卿等［3］采用分布式水文模型SWAT对大通河流域径流序列进行了模拟。白雁翎等［4］基于Budyko 假设的气候弹性系数法分解了降水、潜在蒸发和其他因素对径流变化的影响分量。王大超［5］采用双累积曲线得出人类活动在1995—2007年、2008—2013年对大通河下游径流的影响率分别为92.5%和135.4%。水文模型参数率定和验证过程复杂，不确定性强［6］，基于Budyko 假设的弹性系数法主要在年尺度上定量区分影响因素对径流变化的影响，部分学者在归因识别时因方法本身、参数选取、计算时段不同往往得出不同的结论。以往研究倾向于采用一种方法来估计径流变化的贡献率，结果缺少不同方法的对比验证，且对径流影响因素分析多为定性描述。王随继等［7］提出的累积量斜率变化率法所需数据简单，累积量在一定程度上消除了人为干扰以及实测数据年际波动的影响，被诸多学者应用于径流归因研究中。笔者基于最新的水文气象资料，首次采用累积量斜率变化率法对大通河流域径流变化的气候变化、人类活动等影响因素进行定量评估［8］，同时利用双累积曲线法对结果进行分析验证，以揭示流域水文过程对环境变化的响应关系，研究方法和成果可为流域综合开发治理提供科学依据。

1 研究区概况

大通河为湟水最大的一级支流，发源于青海省天峻县托勒南山，地理位置介于98°30′～103°15′E、36°30′～38°25′N 之间，河源海拔4520 m，干流全长572 km，河口海拔1727 m，河道平均比降4.19%。河源至尕大滩为上游段，河长308.4 km，水系发育，呈树枝状，主要支流有莫日曲、江仓曲、永安河等7条，主要分布在右岸；尕大滩至连城为中游段，河长223.4 km，水系呈羽状，主要支流有白水河、老虎沟、讨拉沟等；连城以下为下游段，河长40.2 km（图1）。大通河流域属内陆高寒气候区，气候特点为冬长夏短、冰冻期长、气候寒冷、年降水较多、雨热同季，流域内年平均气温-0.3～8.0 ℃，极端最高气温35.8 ℃，极端最低气温-34.1 ℃，多年平均降水量500.6 mm。

图1 大通河流域水系及站点分布Fig.1 Distribution of water system and stations in the Datong River Basin

2 数据与方法

2.1 数据来源

大通河出口处布设有享堂水文站，控制流域面积15126 km2，距入湟口1.9 km，该站设立于1939 年10 月，1940 年开始观测，经黄委整编后的完整连续资料自1950年至今。流域和周围10处气象站的逐日降水量、气温、相对湿度、日照时数、风速等资料来源于中国气象数据网（http://data.cma.cn），面降水量、潜在蒸发量、平均气温采用泰森多边形法计算，考虑各站点资料的时间序列和完整性，水文气象要素分析统一采用1961—2020 年系列。冰川分布［9］和地表植被指数产品［10］来源于国家青藏高原科学数据中心（http://data.tpdc.ac.cn），冰川分布产品原始基础数据为1978—2019年祁连山全境的Landsat 系列影像，时间分辨率为5 a，空间分辨率为30 m，用于冰川变化对流域径流影响定量估计。地表植被指数产品空间分辨率为30 m，采用最大值合成方法实现对地表月度归一化植被指数（Normalized difference vegetation index，NDVI）产品的合成，进而计算植被覆盖度（Fractional vegetation cover，FVC）。

2.2 研究方法

2.2.1 潜在蒸发量采用Penman-Monteith 公式[11]计算流域潜在蒸发量：

式中：ET 为潜在蒸发量（mm·d-1）；Δ 为饱和水汽压曲线斜率；Rn为地面净辐射（MJ·m-2·d-1）；G为土壤热通量（MJ·m-2·d-1）；γ为湿度计常数（kPa·℃-1）；Tmean为日平均气温（℃）；u2为2 m高处风速（m·s-1）；es为饱合水汽压（kPa）；ea为实际水汽压（kPa）。

2.2.2 趋势分析及突变检验采用线性倾向估计［12］分析气象水文要素的变化趋势。根据变异诊断分析方法，采用过程线、累积距平进行年径流量的初步诊断，确定序列是否发生变异，通过Mann-Kendall、有序聚类、Pettitt 检验［13-14］详细诊断识别径流的突变年份。采用集中度、集中期、不均匀系数［15-16］等指标反映河川径流年内分配均匀程度。

（1）线性倾向估计

用yi表示样本为n的变量，xi为yi所对应的时间，建立yi与xi的线性回归方程：

式中：a为回归系数；b为回归常数；a、b值可采用最小二乘法估计。

式中：yˉ为径流量序列的平均值；xˉ为时间序列序号的平均值。可求出时间xi与变量yi之间的相关系数r：

（2）有序聚类检验

其实质在于求算最优分割点，即寻找序列突变前后系列离差平方和综和最小的点。

式中：Vτ、Vn-τ分别为分割点前后离差平方和；τ为可能分割点；xt为序列值；为突变点τ前的序列均值；为突变点τ后的序列均值。

式中：Sn(τ)为总离差平方和，其取最小值时对应的τ即为最优分割点。

（3）Pettitt检验

Pettitt 检验为非参数检验法，物理意义清晰，能较好地识别水文气象时间序列中的突变点。对于具有样本量为n的时间序列x，可以构造一秩序列为：

其中：

式中：sk,P为统计量；k为样本总数；P为显著性水平；xi、xj为相同时间序列的2个样本值；ri为符号函数，即第i数值大于j时刻数值个数的累计数。

Pettitt 法是直接利用秩序列来检验变异点的，若t0时刻满足：

式中：kt0为最大统计量，对应的t0点为变异点，计算统计量当P≤0.5 时，则认为检验出的变异点在统计意义上是显著的。

2.2.3 小波周期分析小波分析[17]是在傅里叶变化基础上发展起来的一个数学分支，能清晰揭示出隐含在时间序列中的多种变化周期，被广泛应用于水文气象序列周期识别。小波方差可识别各种尺度的震荡强弱和周期变化特征，方差值越大，表明对应尺度的周期越明显。对于小波函数ψ(t)，径流时序x(t)∈L2(R)的连续小波变换为：

式中：wf(a,b)为小波变换系数；a为尺度因子；b为平移因子；R为实数域；x(t)为径流序列；ψˉ(x)为ψ(x)的共轭。

2.2.4 径流变化归因识别

（1）累积量斜率变化率法

累积量斜率变化率法是在基准期和突变期分别建立影响因子与年份之间的线性方程，从而通过斜率变化率的比值计算出影响因子的贡献率。累积降水量、潜在蒸发量和径流量的斜率变化率分别为：

式中：SP、SET、SR分别为累积降水量、潜在蒸发量、径流量变化期相对于天然时期的变化率（%）；KPa、KPb，KETa、KETb，KRa、KRb分别为累积降水量、潜在蒸发量、径流量突变前后线性关系式的斜率（mm·a-1），正值（负值）表示斜率增大（减小）。

通常情况下径流量与降水量呈正相关，潜在蒸发量与径流量呈负相关：

式中：CP、CET分别为降水量、潜在蒸发量对径流变化的贡献率。

人类活动对径流变化的贡献率CH：

（2）双累积曲线法

采用双累积曲线法对大通河流域径流变化归因结果进行验证。双累积曲线［18-19］是水文学中常用来反映降水径流系统变化的特征曲线，天然情况下累积降水量和累积径流量应为一条斜率固定的直线，当斜率发生显著变化时说明人类活动对径流的影响持续增加，该年份可判定为突变点。采用基准期降水径流关系和变化期面降水量计算仅受降水影响的理论径流量，与基准期的径流量差值为降水变化的影响量，扣除降水因素导致的径流影响量即为人类活动对径流变化的影响量。

3 结果与分析

3.1 年际变化趋势分析

近60 a大通河流域气候呈暖湿化方向发展（表1、图2a～b），年平均气温表现出非常明显的上升趋势，增幅0.42 ℃·(10a)-1明显高于全球［0.26 ℃·(10a)-1］［20］，通过0.001的显著性检验；年降水量呈增加趋势，倾向率为8.9 mm·(10a)-1。潜在蒸发量呈增加趋势，增幅5.6 mm·(10a)-1。年径流量呈减少趋势，变化速率0.67×108m3·(10a)-1，最大值50.19×108m3（1989 年），最小值19.95×108m3（2015 年），极值比2.52，变差系数为0.20，年降水量、潜在蒸发量、年径流量均通过0.05的显著性检验。

表1 大通河流域各要素趋势检验Tab.1 Trend test of elements in the Datong River Basin

图2 各要素序列变化趋势Fig.2 Sequence change trend of each element

3.2 径流年内变化特征

大通河流域径流年内分配的集中度、不均匀系数（图3a～b）呈微弱下降趋势，径流年内分配趋于均匀，汛期（5—10月）径流量占全年比重呈减小趋势，而枯季（11月至翌年4月）径流量显著增加，反映出流域生态环境好转，蓄水和水源涵养能力逐步增强。随时间推移，集中期有推迟趋势，延迟速率为3.0 d·(10a)-1，2011—2020年最大径流量出现时间较20世纪60年代平均延迟18 d。从季节分析，夏季径流量（占比48.1%）远大于其他季节，呈现夏季（6—8月）＞秋季（9—11月）＞春季（3—5月）＞冬季（12至翌年2月）。

图3 享堂站径流量指标年际变化Fig.3 Interannual variation of runoff indicators at Xiangtang Station

3.3 径流周期性分析

享堂站小波系数实部等值线和小波方法曲线如图4a～b 所示，大通河年径流主要存在12 a、24 a、44 a的周期震荡。3个时间尺度中44 a峰值最明显，说明周期震荡最强，为第一主周期，其次为24 a 和12 a。在44 a 左右的特征时间尺度下，年径流经历了2 次枯-丰转化，当前的丰水期直至2021 年等值线仍未闭合，表明偏丰状态仍将持续。

图4 径流小波周期分析Fig.4 Period analysis of runoff wavelet

3.4 径流突变特征分析

从大通河流域降水量、径流量累积距平曲线（图5）初步判断，年径流量在1990年左右发生突变，由丰水期转为枯水期。Mann-Kendall检验表明突变发生在1995 年；有序聚类检验曲线离差平方和在1990、1994 年达到最小值，出现明显的突变；Pettitt检验中1990 年径流突变点统计量P=0.063＜0.5，为有效突变点。综合各检测方法结果，确定大通河年径流在1990年发生突变。20世纪90年代末大通河流域水电开发持续加快，1994年引大入秦调水工程建成通水是引起流域径流突变的主要原因，后外调水量逐年增加，流域径流量明显减少。相较基准期，变化期径流量减少3.52×108m3，减幅12.1%；降水量增加14.5 mm，增幅3.0%；潜在蒸发量增加11.1 mm，增幅1.2%。

图5 大通河流域年径流突变检验Fig.5 Abrupt change of annual runoff in the Datong River Basin

3.5 径流变化定量评估

绘制大通河流域年径流量、降水量和潜在蒸发量的累积曲线，对累积曲线的基准期、变化期分别建立累积量与年份的线性关系，拟合方程（图6），各序列在1990年前后斜率变化率有明显的变化特征，累积降水量、潜在蒸发量的斜率变化率为正，累积径流深的斜率变化率为负，说明降水量、潜在蒸发量序列呈增加趋势，径流深序列为下降趋势。由式（11）～（13）计算各累积序列影响期的斜率变化率，由式（14）～（16）计算气候变化和人类活动对径流变化的贡献率，计算结果见表2。结果表明，气候变化中降水量增加对径流量减少呈负贡献，贡献率-31.4%，潜在蒸发量增加对径流量减少呈正贡献，贡献率为13.7%；人类活动导致径流量减小的贡献率达117.7%。因此人类活动是导致流域径流变化的主导因素，结论与杨大文等［21］的研究结果一致，结果可信。

表2 气候变化和人类活动对大通河流域径流影响Tab.2 Impacts of climate change and human activities on runoff in the Datong River Basin

图6 大通河流域累积年径流深、降水量和潜在蒸发量与年份的拟合关系Fig.6 Fitting relationship between cumulative annual runoff depth,precipitation,evapotranspiration and year in the Datong River Basin

3.6 双累积曲线结果验证

采用双累积曲线验证径流变化归因成果的合理性（图7a～b），突变后的降水量-径流量关系曲线斜率明显小于天然时期，产流量减小。经计算，大通河流域1991—2020 年径流量相比基准期减小3.52×108m3，气候变化导致径流量增加1.40×108m3，影响率为-39.8%，而人类活动导致径流量减小4.92×108m3，影响率达139.8%，即人类活动是导致流域径流变化的主要因素，与累积量斜率变化率法计算的贡献率接近，计算结果可信。

图7 大通河流域降水量-径流量双累积曲线及其关系Fig.7 Double-mass curve and relationship of precipitation-runoff in the Datong River Basin

3.7 径流变化影响因素分析

3.7.1 气象因子对径流的影响分析影响径流变化的主要气象因子包括降水、气温和蒸散发等。大通河流域径流与降水的年内分配规律一致，汛期5—10 月径流量占年径流量的80.4%，降水对径流表现为正贡献，是流域来水的主要补给来源，蒸散发则表现为负贡献。气温升高一方面通过影响蒸散发减少河川径流量，另一方面则引起流域内积雪和冰川消融，间接增加河川径流量。解译流域内冰川面积（表3）发现，1986—2019 年冰川面积整体呈减小趋势，1986年流域冰川面积为37.0 km2，2005年减小至21.3 km2，2015 年以后面积稳定在16.5 km2左右。绘制年径流与年平均气温相关关系为负相关，分析流域平均气温对年径流表现为负贡献。因此，气候条件中降水对径流变化起主导作用，气温、蒸散发、冰川面积影响径流变化但不显著。

表3 大通河流域冰川面积及植被覆盖度统计Tab.3 Statistics of glacier area and fractional vegetation cover(FVC)in the Datong River Basin

3.7.2 植被变化对径流的影响分析大通河流域人烟稀少，开发利用程度低，土地利用类型以林地和草地为主，面积达70%，短时期内可认为下垫面地形、地质和土壤条件相对稳定，因此下垫面变化更多的归因于植被功能变化。植被变化对径流的作用，部分学者认为植被覆盖度（FVC）增加会减少区域径流量，也有学者认为FVC增加会提高区域的产流能力。从表3 可知，1986—2019 年流域植被变化趋势不明显，相关关系分析流域FVC 年值、生长季均值（6—8 月）与径流相关性弱。考虑上述因素得出1990—2020 年大通河流域植被变化对径流的影响较小。

3.7.3 水利工程对径流的影响分析

（1）蓄水工程

经调查，大通河流域蓄水工程主要为水电站，无淤地坝工程。流域内已建水电站40座（含水库两座），总装机容量566.97 MW，分布在中下游河段，建成时间集中在1990—2016 年（图8a），其中干流23座，装机548.2 MW，老虎沟等8 条支流17 座，装机18.8 MW。除2014、2016 年相继建成的两座大型水库纳子峡和石头峡水库对流域径流、洪水有调蓄作用外，其余径流引水式电站对来水调蓄影响较小，水库运行对河道径流的影响表现为对天然来水进行调节，蓄丰补枯，水库下游汛期径流比重、集中度减小，枯季径流增加（图8b）。

图8 水电站装机容量累积及不同年代月径流占比Fig.8 Accumulation of installed capacity of hydropower stations and distribution of monthly runoff in different years

（2）调水工程

根据大通河水量分配方案，大通河设计外调水量14.43×108m3，其中引大济西2.50×108m3，引大济湟7.50×108m3，引大入秦4.43×108m3。随着引大济西一期引硫济金、引大入秦工程相继建成发挥效益，流域外调水量逐年增加（图9a）。其中引硫济金工程从大通河支流硫磺沟调水至甘肃省金昌市，设计引水规模0.40×108m3，2003 年5 月建成通水，2003—2021年平均调水量0.25×108m3，2011年达到设计调水量。引大入秦是在天堂水文站下游500 m将大通河水引至秦王川，设计引水量4.43×108m3，1994 年10 月工程建成通水，1994—2021 年平均调水量1.91×108m3，年最大调水量3.70×108m3（2021年）。引大济湟工程目前正处于建设阶段，2019—2021 年向湟水上游宝库河平均调水1.17×108m3。通过对享堂站实测径流量进行还原，经肯德尔秩次检验，享堂站天然年径流呈显著增加趋势（图9b）。2019、2020年流域外调水量占享堂站径流还原量的90%以上，因此跨流域调水是导致流域径流量减少最主要的关键因素。

图9 调水工程引水量及享堂站天然年径流量Fig.9 Diversion of water diversion project and natural annual runoff hydrograph of Xiangtang Station

4 结论

本文利用大通河流域气象水文等数据，分析了近60 a径流演变特征，运用累积量斜率变化率法和双累积曲线法对径流变化进行归因分析，得出主要结论如下：

（1）大通河流域气候暖湿化明显，气温、降水量、潜在蒸发量增幅分别为0.42°C·(10a)-1、8.9 mm·(10a)-1和5.6 mm·(10a)-1，年径流量呈减少趋势，倾向率为0.67×108m3·(10a)-1，在1990 年发生突变，突变后径流量减少3.52×108m3，减幅12.1%。

（2）大通河流域径流集中度和不均匀系数呈微弱下降趋势，枯季径流增加趋势明显，径流年内分配趋于均匀，集中期有推迟趋势，延迟速率为3.0 d·(10a)-1。

（3）大通河年径流具有12 a、24 a及44 a左右的震荡周期，在44 a左右尺度上周期震荡最为显著，为第一主周期。最近丰水期的小波等值线还未闭合，表明未来一段时间年径流将持续偏丰。

（4）气候、人类活动对大通河径流减少的贡献率分别为-17.7%和117.7%，降水是流域来水的主要补给来源，跨流域调水是引起径流减小的主要人为因素。