黄河源区水文要素变化及驱动因素分析

张 鸣，朱 奎，鲁 帆，戴雁宇，宋昕熠，于嵩彬

（1.中国矿业大学 资源与地球科学学院，江苏 徐州 221008；2.中国水利水电科学研究院 水资源研究所，北京 100038；3.长沙理工大学 水利与环境工程学院，湖南 长沙 410114）

0 引言

黄河是我国仅次于长江的第二大河，多年平均水资源总量为647 亿m3，不到长江的7%，但水资源开发利用率高达80%，远超40%的生态警戒线，水资源短缺是黄河流域最大的矛盾［1－3］。黄河源区唐乃亥水文站的多年平均径流量占黄河水资源总量的34.5%［4］，黄河源区的水资源状况不仅关系当地的经济社会发展和生态环境保护，而且对整个黄河流域水资源开发利用有着重要影响。近年来，已有不少学者对黄河源区径流变化进行了研究，苏中海等［5］、刘希胜等［6］认为黄河源区上游径流量呈现不显著增大趋势，而刘晶等［7］、Yan 等［8］认为黄河源区径流量整体呈减小趋势。各学者研究选取的时间段和数据不同，结论也有所不同。对于径流驱动因素，商滢等［9］采用双累积曲线得出非降水因素是径流变化的主要驱动因素的结论，张士锋等［10］通过建立驱动模型得出降水、潜在蒸散发对径流变化分别起正向、负向驱动作用的结论。上述学者研究径流时多将研究期看作一个整体，较少考虑研究期内径流变化过程，驱动因素分析大都只考虑部分气象要素和人类活动的影响，较少综合考虑气象要素和下垫面的影响并定量分析。黄河源区位于素有“亚洲水塔”和“第三极”之称的青藏高原，对气候变化较敏感，在全球变暖的环境下以往的研究结论已不完全适用，需补充研究论证。

目前采用Budyko 法、水文模型、双累积曲线法均能定量揭示气候变化和人类活动对径流的影响。采用Budyko 法便于从气候变化中分离出不同影响因子。近年来，国内外学者采用Budyko 法建立多种经验公式并在长江和黄河等流域广泛应用［11－15］，但传统Budyko法下垫面参数在不同时段会出现跳跃式的变化，这与实际情况存在一定偏差。本文基于1959—2019 年黄河源区水文气象数据，分析黄河源区水文要素演变规律，运用累积距平曲线和Mann－Kendall 突变检验法确定径流突变年份，采用变参数的Budyko 公式定量分离降水、潜在蒸散发和下垫面变化对径流量变化的贡献率，以期为黄河流域合理开发利用水资源、建设水利工程以及防汛抗旱提供参考。

1 研究区概况和数据来源

1.1 研究区概况

黄河流域位于东经96°—119°、北纬32°—42°，发源于青藏高原巴颜喀拉山北麓的约古宗列盆地，流域地势西高东低，西部河源地区平均高程在4 000 m 以上，源区涉及青海、四川、甘肃3 个省的6 个州18 个县。黄河源区唐乃亥水文站控制集水面积为12.2 万km2，黄河源区是世界上高原水体与湿地集中分布的地区之一。研究区内水系、湖泊、10 个气象站和唐乃亥水文站分布见图1，为增加站网密度，选取的部分气象站位于集水面积之外。

图1 黄河源区范围

1.2 数据来源

使用的唐乃亥水文站数据源自《黄河水文年鉴》《黄河水资源公报》《黄河泥沙公报》，10 个气象站数据源自中国气象数据网，1982—2019 年归一化植被指数（精度为5 km）源自国家地球系统科学数据中心，1980—2020 年土地利用类型遥感监测数据（精度为1 km）源自中国科学院资源环境科学与数据中心。

2 研究方法

2.1 趋势检验和突变检验

趋势检验采用线性回归分析法和Mann－Kendall检验法（简称M－K 检验）［16］。线性回归分析法可以反映变量随时间的变化趋势。M－K 检验是一种非参数检验法，其优点是不需要样本遵从一定的分布规律，也不受少数异常值的干扰。

突变检验采用M－K 检验法和累积距平法，累积距平法是根据时间序列累积值与均值的差值绘制累积距平曲线，能比较直观地判断样本的变化趋势和突变年份。

2.2 潜在蒸散发量与植被覆盖度计算

采用FAO Penman－Monteith 方法计算流域潜在蒸散发量：

式中：E0为潜在蒸散发量，Tmean为日均气温，Rn为净辐射，G为土壤热通量，γ为干湿表常数，u2为2 m 高风速，es、ea分别为饱和、实际水气压，Δ为饱和水气压曲线斜率。

李苗苗等［17］基于像元二分模型建立植被覆盖度计算公式：

式中：fc为植被覆盖度，N为归一化植被指数，Nsoil为裸土或者无植被覆盖区域的归一化植被指数，Nveg为完全被植被覆盖区域的归一化植被指数。

由于存在噪声等其他因素的影响，因此取5%的N值为Nsoil，取95%的N值为Nveg。为更好地反映研究区内植被变化情况，归一化植被指数可选取年内最大值。

2.3 径流归因分析

基于Budyko 假设的长时间序列中，水分输入可看作降水量P，能量输入可看作潜在蒸散发量E0，考虑水量平衡和能量平衡的水热耦合方程为

式中：E为流域蒸散发量。

对于闭合流域来说，长期水量平衡方程为

式中：Q为径流量，ΔS为流域水分蓄变量。

考虑到区域气候、地理环境等的差异，国内外学者基于Budyko 理论陆续建立了各种考虑下垫面参数的经验公式［18－21］，本文选用Zhang 等［18］建立的经验公式：

式中：ε为干旱指数，ε＝E0/P；n为下垫面参数。

采用式（5）分析径流变化时，下垫面参数和水文要素会在不同时段呈现跳跃性变化，这与实际情况存在偏差。Jiang 等［15］提出了变参数的Budyko 公式，选取气象和人类活动等要素为协变量构建下垫面参数的相关函数，假定选取m个自变量，滑动平均序列为y1，y2，…，yt，第i（i＝1，2，…，t）年的m个自变量为，，…，，可选用一次函数、二次函数、指数函数等拟合下垫面参数和m个自变量之间的关系，第i年的径流量Qi为

式中：Pi为第i年降水量，εi为第i年干旱指数，βi为第i年回归系数，ei为第i年残差。

假定ei符合正态分布N（0，σ2），则Qi符合正态分布，根据参数β和σ的似然函数求出极大似然估计量和，对应有ni＝Yi，ni为第i年的下垫面参数。

为定量描述气象因子和下垫面参数对径流量变化的影响，引入弹性系数φ来表示径流量对某因子x变化的敏感程度：

假设变化期和基准期内降水量、潜在蒸散发量、下垫面参数变化分别为ΔP、ΔE0、Δn，根据基准期信息推求变化期内径流变化量：

3 水文要素变化

基于唐乃亥水文站实测径流数据，分析1959—2019 年黄河源区年径流量变化情况，见图2（a）。M－K趋势检验结果显示，整体上年径流量呈现减小趋势，线性减小速率约为2.36 亿m3/10 a，多年实测年均径流量为203.9 亿m3，1959 年年径流量（为328.08 亿m3）和2002 年年径流量（为105.74 亿m3）分别为61 a样本系列中径流量的最大值和最小值。利用M－K 检验方法对径流序列进行突变检验，见图2（b），在0.05显著性水平下，UBK曲线与UFK曲线在1987—1990 年置信区间内多次交汇。为明确径流突变年份，引入年径流量累积距平曲线，见图2（c），年径流量在1959—1989 年整体呈现上升趋势，在1989—2019 年整体呈现下降趋势，在1989 年发生一次明显突变，综合判断1989 年为径流突变年份。为进一步反映近年来全球变暖环境下黄河源区的径流变化，绘制1990—2019 年年径流量累积距平曲线，见图2（d），2004 年为第2 个突变年份。总体上，从1959—1989 年、1990—2004 年和2005—2019 年3 个时段来看，年径流量呈现先下降后上升趋势。

图2 1959—2019 年黄河源区径流量变化

对1959—2019 年黄河源区年内径流量进行分析，1981 年9 月月径流量最大（为92.02 亿m3），2003 年1月月径流量最小（为2.33 亿m3），每年7 月月均径流量最大（为35.04 亿m3），每年1 月月均径流量最小（为4.58 亿m3），7—9 月径流量占全年总径流量的46.97%，而1—3 月径流量仅占全年总径流量的7.32%，由此可见，黄河源区径流量季节分配不均。径流集中度指月径流量以矢量形式累积，其各分矢量之和占年径流量的百分比，反映年内径流量的集中程度。集中期反映全年径流量集中的月份［22］。唐乃亥水文站年内径流集中度和集中期计算结果见图3，集中度为35%～55%，整体上呈现下降趋势，下降率为0.716%/10 a。集中期为7—9 月，整体上呈现略微提前的趋势，对于唐乃亥水文站和下游水库来说，洪峰有提前的趋势。

图3 唐乃亥水文站年内径流集中度和集中期计算结果

基于10 个气象站降水数据分析黄河源区1959—2019 年降水量变化趋势，见图4。1959—2019 年多年平均降水量为504.83 mm，1962 年降水量（为398.24 mm）和2018 年降水量（为674.21 mm）分别为61 a 样本系列中的最小值和最大值。M－K 趋势检验结果显示，统计量Z为3.72＞2.56，降水量呈现显著增加趋势（显著性水平α＝0.01），平均线性增加速率为9.1 mm/10 a。

图4 1959—2019 年降水量变化

径流系数代表流域产流能力，综合反映流域内自然地理要素对降水—径流关系的影响。1959—2019年径流系数变化情况见图5，1959—1989 年、1990—2004 年和2005—2019 年对应的径流系数均值分别为0.36、0.28 和0.32，径流系数呈现先减小后增大的趋势，与径流量变化趋势相同。

图5 1959—2019 年径流系数变化情况

4 驱动因素

4.1 降水、气温和潜在蒸散发分析

基于黄河源区10 个气象站数据对降水量进行Kriging 插值分析，得出多年平均降水量的空间分布，见图6（a），黄河源区1959—2019 年多年平均降水量在空间上呈现由西北向东南递增的趋势。M－K 突变检验结果［见图6（b）］和降水量累积距平曲线［见图6（c）］显示，2002 年和2016 年为降水序列的突变点。1959—2002 年、2003—2016 年和2017—2019 年多年平均降水量分别为490.67、525.97、613.84 mm，尤其是近年来降水量增加趋势愈发明显。

图6 降水量变化

基于10 个气象站实测数据分析黄河源区1959—2019 年年均气温变化趋势，见图7（a）。多年平均气温为－0.66 ℃，1965 年年均气温最低（为－2.02 ℃），2016 年年均气温最高（为0.73 ℃），M－K 趋势检验统计量Z为11.40＞2.56，气温呈现显著上升趋势（α＝0.01），平均线性上升速率为0.3 ℃/10 a。由气温累积距平曲线［见图7（b）］可以看出，气温累积距平值在1997 年发生转折，结合趋势检验结果，气温在1959—1997 年呈现不显著上升趋势，在1997 年之后开始显著上升。

图7 1959—2019 年气温变化趋势与气温累积距平曲线

基于10 个气象站实测气温、风速、日照和相对湿度数据，利用FAO Penman－Monteith 方法计算各气象站潜在蒸散发量，采用Thiessen 多边形方法计算出黄河源区年平均潜在蒸散发量，见图8（a）。1959—2019年多年平均潜在蒸散发量为724.8 mm，1965 年潜在蒸散发量最小（为672.4 mm），2010 年潜在蒸散发量最大（为784.5 mm）。M－K 趋势检验统计量Z为7.89＞2.56，潜在蒸散发量呈现显著增加趋势（α＝0.01），平均线性增加率为8.2 mm/10 a。潜在蒸散发量主要受日照时数、风速、气温和相对湿度的影响，因此利用多元线性回归法对每个气象站的潜在蒸散发量进行归因分析。结果表明，气温是90%气象站潜在蒸散发量的主要影响因素，风速是10%气象站潜在蒸散发量的主要影响因素。M－K 突变检验结果［见图8（b）］显示，潜在蒸散发量在1997 年发生一次突变，突变时间点与气温的相同，说明气温与潜在蒸散发量之间具有较强的相关性。

图8 1959—2019 年潜在蒸散发量变化趋势与M－K 突变检验

为分析1959—2019 年黄河源区降水量、气温、潜在蒸散发量与径流量的相关性，对其进行单因素相关性检验，结果显示，径流量与降水量之间的皮尔逊系数为0.778，具有极强的相关性，远远大于气温的和潜在蒸散发量的。

4.2 下垫面变化分析

下垫面主要指与大气直接接触的地球表面，本文主要分析黄河源区植被、冰雪和冻土的变化情况。植被覆盖度是衡量地表植被变化的重要指标，黄河源区植被覆盖度变化趋势与累积距平曲线见图9。M－K 趋势检验统计量Z为2.01＜2.56，植被覆盖度呈现不显著上升趋势（α＝0.01），这与管晓祥等［23］、刘启兴等［24］、韩思淇等［25］的研究结果一致。在空间上植被覆盖度呈现由西北向东南递增的趋势，与降水量的空间分布一致，这与高思琦等［26］认为黄河源区植被的空间分布主要受降水影响的观点一致。

为分析黄河源区内不同植被类型的变化情况，对1980—2020 年土地利用类型数据进行分析，见图10。源区内土地利用类型主要为草地，城镇用地和湖泊面积变化不大，而林地、草地分别在1990—1995 年、1996—2000 年呈现略微退化趋势。2000—2013 年未利用土地面积呈现显著减小趋势，而林地面积和草地面积呈现显著增大趋势。植被变化是人类活动和气候因子共同作用的结果，在生态环境较为敏感的青藏高原地区植被变化更为显著。气候变化对植被的影响主要通过降水量和气温的变化改变土壤湿度和温度，调节植物生长过程中酶的活性，影响植物光合作用和呼吸作用［27－29］。降水量增加会促使高寒草原和荒漠植被生长，增加植被面积，气温升高则促使植被覆盖度提高［30］，植被的改善会增强区域蒸散发，进而增加降水量。

图10 1980—2020 年土地利用类型变化

归一化植被指数是表征区域植被变化的重要指数，1982—2019 年黄河源区降水量与归一化植被指数的相关性分析及显著性检验结果见图11。78%的区域降水量对归一化植被指数起正向促进作用，两者成显著正相关的区域主要集中在降水相对贫乏的西北部地区和东北部地区；21%的区域降水量与归一化植被指数负相关，两者成显著负相关的区域很少，植被类型以林地和耕地为主，降水较充沛，这些地区的植被变化可能受人类活动影响较大，这与史丹丹等［31］的研究结果一致。管晓祥等［23］研究发现黄河源区气温与归一化植被指数显著正相关，且春季的影响最大。黄河源区多年平均降水量为504.83 mm，降水较充沛，因此水量不是限制植被生长的主要条件。黄河源区高程普遍在4 000 m以上，热量不足限制植被生长，因此气温是影响黄河源区植被变化的主要因子［23，27－28，32］。此外，黄河源区的人类活动也对植被变化产生一定影响，如过度砍伐和放牧、退耕还林、减畜工程和成立三江源国家自然保护区等。

图11 1982—2019 年降水量与归一化植被指数相关性分析

黄河源区冰雪面积变化情况见图12，累积距平曲线显示，1982—1997 年冰雪面积呈现增加趋势，1997—2015 年冰雪呈现消融趋势，冰雪面积变化的突变时间为1997 年，与气温的一致，表明气温是控制冰雪变化的主要因素［33］。

图12 1982—2015 年冰雪面积占比累积距平曲线

关于黄河源区冻土的变化情况，曹慧宇［34］的研究表明2003—2019 年多年冻土面积由3.59 万km2减小到3.42 万km2，4.82%的多年冻土转化为季节性冻土，季节性冻土最大冻结深度减小16 cm，活动层厚度增加9.28 cm。冯雨晴［35］的研究表明1961—2015 年黄河源区冻土融化夹层厚度和活动层厚度呈现增加趋势，1997 年之后呈现显著增大趋势。多年冻土的退化以及活动层厚度的增加导致土壤渗透性增强，地表径流更多入渗成为地下水，土壤含水量整体增加，蒸散发量增加，流域地下水储量增加，从而减少夏季径流量，增加冬季径流量和基流量，调节径流量的季节性分布［36－40］。此外，活动层底部的冻结层上水可能侧向补给江河湖泊，黄河源区地下冰融水对地表径流的补给比例为13.2%～16.7%，但其融水对径流的改变较小［41－44］。

4.3 Budyko 假设定量分离

根据径流突变结果，以1989 年和2004 年为突变年份将研究期分为3 个时段，3 个时段的多年平均径流深、降水量、潜在蒸散发量、气温见表1。

表1 黄河源区各时段径流深和协变量

为消除气候变化及地下水储量对水量平衡的影响，构建下垫面参数时滑动平均步长取9 a，利用式（5）计算滑动平均后的下垫面参数n。基于Jiang等［15］提出的变参数Budyko 公式，选取降水量、潜在蒸散发量和气温作为协变量，选用线性函数进行拟合，下垫面参数n与降水量、潜在蒸散发量和气温的关系式为

式（9）标准化系数显示降水量是下垫面参数的主要影响因素，其次是气温。根据式（5）和式（4）利用气象数据和径流数据推求气候变化造成的径流变化量dQc，再根据dQh＝dQ－dQc计算分析人类活动对径流的影响（dQh为人类活动对径流的改变量，dQ为径流总改变量）。结果显示，在第1 个突变期（1959—2004年）气候变化对径流量减小的贡献率为86.63%、人类活动对径流量减小的贡献率为13.37%，在第2 个突变期（1990—2019 年）气候变化对径流量增加的贡献率为101.18%、人类活动对径流量增加的贡献率为－1.18%。在整个研究期内气候变化是径流量变化的主导因素，人类活动对径流量的增大一直起负向作用，但对径流量变化的影响程度明显降低。

为定量分离不同气象因子和下垫面变化的贡献率，分别计算降水量、潜在蒸散发量和下垫面参数在不同时期对应的弹性系数，结果见表2。基准期和变化期（分别以1959—1989 年为基准期、1990—2004 年为变化期，1990—2004 年为基准期、2005—2019 年为变化期）降水量的弹性系数始终大于潜在蒸散发量的和下垫面参数的，说明径流对降水的响应最为敏感。

表2 黄河源区下垫面参数和各协变量对应的弹性系数

根据式（8）计算突变期各因子对径流量变化的贡献率，在第1 个突变期内下垫面参数n对径流量变化的贡献率达到60.69%，是径流量减小的主要驱动因素，降水其次（为31.75%），潜在蒸散发的贡献率最小（为7.56%）；第2 个突变期内降水对径流量变化的贡献率达到103.91%，是径流量增加的主要驱动因素，潜在蒸散发其次（为－10.37%），下垫面变化的贡献率最小（为6.46%）。此外，计算得出2 个突变期内各因子引起的模拟径流变化量之和与实测径流变化量的误差分别为3.98%和0.1%，误差较小，表明结果具有准确性。

5 结论

1）黄河源区1959—2019 年径流整体呈现先下降后上升的变化趋势，多年实测年均径流量为203.9 亿m3，年内径流集中度为35%～55%，集中期位于7—9 月，年内径流量季节分配不均，1989 年和2004 年为径流量突变年份。1959—2019 年年降水量和潜在蒸散发量都呈现显著增加趋势，气温在1997 年之后显著上升。1982—2019 年植被覆盖度呈现不显著上升趋势，在空间上植被覆盖度呈现由西北向东南递增的趋势，与降水量的空间分布一致。

2）基于变参数的Budyko 方法定量分离不同因子对径流量变化的影响，气候变化是黄河源区1959—2019 年径流量变化的主要影响因素，径流对降水的响应最为敏感。1959—2004 年下垫面参数n对径流量变化的贡献率达到60.69%，是径流减小的主要驱动因素，降水其次，潜在蒸散发的贡献率最小。1990—2019年降水对径流量变化的贡献率达到103.91%，是径流增大的主要驱动因素，潜在蒸散发次之，下垫面的贡献率最小。

3）径流量与降水量之间的皮尔逊相关系数为0.778，远远大于气温的和潜在蒸散发量的，因此1959—2019 年黄河源区径流量改变最主要的影响因素是降水，其次为气温。