摘 要：
长江上游镇江关流域作为生态脆弱区，其径流驱动机制研究对于流域的发展十分关键。当前对该流域径流分析的研究多集中于水沙关系和径流时空变化分析方面，对于径流归因分析的研究不足。针对镇江关流域径流变化成因复杂的问题，为了明确河川径流变化的主导因素，基于Budyko假设建立了适用于镇江关流域的水热耦合模型。根据下垫面变化强弱划分径流序列的不同阶段，计算反映镇江关流域下垫面特征的参数ω,进而量化气候变化和下垫面变化对径流产生的影响。结果表明，1985—2015年间镇江关流域径流量呈现总体上升趋势，其中2003—2008年间气候变化、下垫面变化的贡献率分别为65.6%和34.4%;2009—2015年间气候变化、下垫面变化的贡献率分别为57.6%和42.4%,气候变化仍是影响径流变化的主导因素。研究结果有利于深入认识变化环境对流域径流变化的影响，对未来水资源系统规划设计与流域水土保持工作具有重要意义。
关键词：
Budyko假设;归因分析;径流变化;气候变化;下垫面变化;人类活动;长江流域;数值模拟;
作者简介：
张成(1982—),男，副教授，博士，主要从事水力学及河流动力学研究。
*张尚弘(1977—),男，教授，博士，主要从事流域水文与水动力学模拟、数字流域仿真、决策支持系统等方面的研究。
基金：
国家自然科学基金项目(52061135104);
中央高校基本科研业务费专项资金(2019MS030);
引用：
张成, 吴传森, 蒯升阳, 等. 基于水热耦合模型的镇江关流域径流驱动机制研究[ J] . 水利水电技术(中英文), 2022, 53(8): 78- 87.
ZHANG Cheng, WU Chuansen, KUAI Shengyang, et al. Water-heat coupling model-based study on runoff driving mechanism of Zhenjiangguan Watershed[J]. Water Resources and Hydropower Engineering, 2022, 53(8): 78- 87.
0 引 言
水资源作为人类必不可少的资源之一，是支撑人类社会发展的基础。随着气候变化和人类活动影响的加剧，水资源问题变得更加突出。在气候变化方面，全球平均气温的持续升高会导致冰川融化，加速水资源的区域性迁移，使得水资源在时空分布上发生较大改变；人类活动对流域水文过程的影响主要是通过土地利用、植被覆盖变化、水利工程建设等方式改变下垫面特征，使其产流机制发生改变，因此通过对下垫面变化的研究能够间接反映人类活动对水文过程的影响。气候变化和下垫面变化对水文水资源的影响过程十分复杂,开展变化环境下径流响应特征与机制的研究能够厘清气候变化和下垫面变化对流域径流变化的贡献关系，为水资源优化配置给予科学的理论依据，进而对水资源综合管理和发展方向提供决策支持。
镇江关流域位于长江上游，是长江上游生态屏障的重要组成部分，生态环境地位显著。近几十年来气候条件的持续变化和人类活动的不断加强，使流域内的河川径流规律发生了很大的变化。其中气候变化主要体现在由于气温和降水的变化使得流域内径流量随之减少，年径流量1975年为22.82亿m3,2002年减少到10.17亿m3。而人类活动体现在因人为导致的土地利用和植被覆盖变化，自1958年后，人类过度开垦、城市建设等行为使得流域内的植被覆盖率急剧减少以及土地利用类型大面积变化；至80年代中期，开始实施水土保持措施，植被覆盖情况有所好转，水土流失现象也得以改善。因此，研究以及量化气候变化和人类活动对镇江关流域径流的影响，能够为流域水资源分配提供科学理论支持，对促进长江上游区域可持续发展具有现实意义。
径流变化是水文循环过程演变的综合体现，也是气候变化和下垫面变化综合影响的结果。为了揭示径流演变规律，学者们采用了许多方法量化气候变化和下垫面变化对径流的影响，其方法大致可分为三类：基于数理统计的经验对比法、基于弹性系数的分析法和基于物理过程的水文模型法。其中基于Budyko假设的弹性系数分析法使用降水、潜在蒸发量、土地利用要素、植被覆盖要素等参数并结合水量平衡方程的方式，其物理意义明显，参数容易获取，能够有效量化气候变化和下垫面变化对径流的影响，被许多学者广泛使用。赵娜娜等基于Budyko假设，应用敏感性分析方法，对若尔盖流域径流变化进行归因分析发现流域下垫面特征变化对若尔盖流域变化的贡献达90%以上；商放泽等分析了三江源地区的径流变化，发现降雨量是三江源地区气候变化影响径流变化的主要控制因子，植被覆盖是下垫面变化影响径流变化的主要控制因子；WANG等基于Budyko假设使用55 a数据定量分析了美国地区下垫面变化和气候变化对径流的贡献率；LI等将水量平衡原理和Budyko假设相结合，将其运用于不同丰枯情境下气候变化和下垫面变化对径流变化的影响分离。结果表明，气候变化和下垫面变化对河川径流变化的影响因地而异，且区域差异性明显。在镇江关流域，国内学者对于径流变化的研究大多停留在径流的变化趋势和空间特征分析上，关于径流变化成因分析方面的研究，仅有张会兰等通过建立BPCC分布式水文模型研究了径流对气候波动和覆被变化的响应。因缺少足够相关的案例进行比较分析，对于是何种因素导致的径流变化尚不明确。同时相比较基于Budyko的水热耦合模型，分布式模型的计算需要更多参数的支持，因此本文对镇江关流域构建Budyko模型尝试在少量参数的条件下进行径流变化归因分析，比较不同研究方法下的结果，有利于更深入地认识该流域的径流变化驱动因素。
本文基于Budyko假设建立镇江关流域水热耦合平衡模型，通过傅式参数ω研究气候和下垫面变化对径流的影响，厘清二者对径流变化的贡献率，为后续的流域综合管理和发展方向提供参考。
1 研究区域和使用数据1.1 研究区域
镇江关流域位于四川省松潘县境内，是岷江的发源地，属于长江支流的一部分(见图1)。流域面积4 486 km2,海拔2 445～5 360 m, 主河道长103.7 km。受东南季风的影响，全年降雨量时空分布不均，5—10月的降雨量占全年的80%,多年平均降雨量570 mm。流域出口的镇江关水文站多年平均流量为55.0 m3/s, 最大流量为410 m3/s, 最小流量为9.2 m3/s。流域内植被类型主要有草地和林地，各占流域面积的60%和30%以上，其余为少量的建筑用地、水域、耕地和裸岩。截至2019年，镇江关镇户籍人口为3 564人，产业主要以农业、畜牧业为主。流域上游河道比降大，土壤结构松散，对自然环境变化的反应十分敏感，是我国重要的生态系统脆弱区域。
图1 镇江关流域
1.2 使用数据
本文采用30 m精度的DEM数据、0.05°精度的全国NDVI数据、250 m精度的全球土地利用类型图作为空间数据。水文数据节选自《藏南滇西河流水文年鉴》中1985—2015年的径流数据。气象数据为1985—2015年的日尺度数据，具体信息如表1所列。
2.1 径流序列突变点分析
本文通过分析流域的下垫面变化寻找下垫面特征突变的年份，同时采用Mann-Kendall分析法分析镇江关流域径流的趋势变化特征。综合考虑下垫面特征前后变化较大的时间点和径流序列的突变点来确定时间序列中的可能突变点，并将该突变点作为径流变化归因分析的时间突变点。通过时间突变点将研究时间序列划分成两个时期，将突变点前的序列称为自然期(下垫面变化较少，河流保持自然条件的状态),突变点后的序列称为影响期(下垫面变化影响明显，径流量发生改变),进而分析河道径流从自然期过渡到影响期发生的变化。
2.1.1 MK趋势检验
Mann-Kendall趋势检验法是一种应用广泛的非参数统计检验方法。其基本原理如下，定义统计量Sk表示第i个样本xi>xj(1≤j≤x)的累计数，计算公式如下
统计量S的方差Var[S]为
定义趋势分析统计量Z通过下式计算
在双边趋势检验中，对于给定的置信水平α,若统计量|Z|≥Z1−α/2,则在置信水平α上时间序列Xn具有明显的上升或下降趋势。
2.1.2 MK突变检验
径流变化是一个不稳定且不连续的变化过程，而检验其变化的常用方法之一就是Mann-kendall突变检验方法,检验方法如下。
定义统计量Sk为第i个样本xi大于xj(1<j<i)的个数累积量
Sk的均值和方差分别为
定义Uk为Sk标准化后的统计量为
将此方法应用到时间序列Xn的逆序列中计算UBk。给定显著性水平α,将UFk和UBk两个统计量曲线和显著性水平线绘在同一张图上，若UFk和UBk两条曲线出现交点，则交点对应的时刻即是突变可能出现的时间。
2.2 建立Budyko水热耦合平衡模型
建立Budyko水热耦合平衡模型首先计算镇江关流域的潜在蒸发量，国内外许多学者研究表明，联合国粮食及农业组织(FAO)在1998年推荐的Penman-Monteith 公式是目前计算潜在蒸发适用性较好的方法,因此本文选用Penman-Monteith 公式进行计算。利用Budyko假设的解析式——傅抱璞公式试算得到镇江关流域逐年的下垫面参数ω,结合M-K分析得到的时间序列突变点分别计算自然期和影响期的参数ω,使得水热耦合模型能够模拟不同时期的水文变化过程，进而分析由气候变化和下垫面变化引起的径流变化。
2.2.1 Budyko假设
Budyko认为在一定的气候及下垫面条件下，流域长期(年及年尺度以上)的水分和能量之间存在耦合平衡关系，即水热耦合平衡。在年或多年尺度上，用降水量代表地表水分供应条件，用潜在蒸散发代表地表能量供应条件，并以此为理论基础建立了水热耦合平衡模型，其表示方式如下
式中，E为年总实际蒸散发量(mm);P为年总降水量(mm);E0为年总潜在蒸散发量(mm)。
在Budyko研究的基础上，傅抱璞[26]根据流域水文循环的机制提出了一个基于Budyko假设的解析表达式
式中，ω为一个无量纲参数，代表流域特征，能够反映地形、土壤性质、植被之间复杂的相互作用。其余参数含义同上。
采用Penman-Monteith公式计算每日的潜在蒸发量Ei,进而得到年总潜在蒸发量E0
式中，Ei为计算得到的每日潜在蒸发量(mm);Δ为饱和水汽压温度曲线上的斜率(kPa/℃);Rn为作物表面的净辐射[MJ/(m2·d)];G为土壤热量通量密度[MJ/(m2·d)],当计算尺度为天或以上时，G的取值为0,T为2 m高的日均空气温度(℃);U2为2 m高的风速(m/s);es为饱和蒸气压(kPa/℃);ea为实际蒸气压(kPa/℃);γ为湿度计常数[MJ/(m2·d)]。
利用径流资料计算流域逐年的年径流深，在年尺度上忽略流域储水量的年际变化，由水量平衡方程计算出流域的年实际蒸发量，公式如下
式中，R为年总径流深(mm)。其余参数含义同上。
使用多年年降水量、年径流深、由水量平衡方法的得到的年实际蒸发量和由P-M公式计算得到的年潜在蒸发量，根据式(11)拟合流域的参数ω并分析其变化规律。
2.3 气候与下垫面变化下的径流响应特征分析
基于Budyko假设的傅抱璞公式，可以研究气候变化和下垫面变化对径流变化的贡献。由于傅式参数ω是表征流域下垫面状态的参数，基于此假设使用自然期计算的参数ω来代表自然期的下垫面条件，代入式(11)来模拟自然期的蒸散发量ES,并利用水量平衡公式式(13)计算自然期的径流量RS,在保证参数ω不变的条件下通过改变影响期和自然期气候条件计算相应时期的径流量，使其相减得到因气候变化导致的径流量变化量。公式如下
同时因为模拟径流量RS使用同时期的气候数据计算，其气候条件与同时期实测径流量R相同，在保证气候条件不变的前提下使用自然期参数ω计算的模拟值减去径流量实测值得到因下垫面变化导致的径流变化量。公式如下
式中，C气候变化为气候变化因素对径流变化的贡献率；R气候变化为因气候因素导致的径流改变量(mm);C下垫面变化为下垫面变化因素对径流变化的贡献率；R下垫面变化为因下垫面因素导致的径流改变量(mm);R总变化为自然期至影响期的径流改变量；RS-自然为使用自然期参数ω计算的自然期模拟径流量(mm);RS-影响为使用自然期参数ω计算的影响期模拟径流量(mm);R实测为实测径流量(mm);R气候变化为因气候变化导致的径流减少量(mm);R下垫面变化为因下垫面变化导致的径流减少量(mm)。
基于以上原理即可分离气候变化和下垫面变化对径流量的影响，进而量化气候与下垫面对径流变化的贡献。
3 结 果3.1 径流序列突变点分析
3.1.1 流域气候及下垫面变化分析
图2为镇江关流域1985—2015年历史降雨及平均气温曲线，流域内降水及气温呈总体上升趋势。其中降水序列的最大值为794 mm(2012年),最小值为446 mm(1985年),多年平均降雨量为577 mm; 平均气温最大值为8.87 ℃(2011年),最小值为6.91 ℃(1992年),多年平均气温为7.85 ℃。
图2 镇江关流域历史降雨及气温曲线
通过研究镇江关流域1995—2015年的土地利用变化发现，林地是镇江关流域最大的土地利用组成成分，约占总面积的60%,主要分布于镇江关流域的下游部分(见图3)。其次是草地，约占总面积的34%,耕地和灌木占总面积3%和2%,城市面积占比不足1%。根据土地利用类型遥感解译成果可知，2010年首次出现了较为明显的城市利用类型，2005年以前该区域城市化发展较慢。
图3 镇江关流域土地利用变化
通过与1995年土地利用类型比较(见图4),2000—2015年间耕地的土地利用类型增加，最大涨幅达到了流域总面积的1.1%,同时随时间呈减少的趋势变化；草地面积略微减少，林地面积持续增加，意味着从1995年后有部分的草地和耕地转化为了林地。
图4 镇江关流域土地利用变化比例
为了更加深入镇江关流域下垫面变化的研究，本文在考虑了土地利用变化的同时又选用归一化植被指数(NDVI)来分析流域植被变化(见图5)。
图5 镇江关流域NDVI变化
图中为镇江关流域的年平均NDVI变化图，以五年为一个跨度计算NDVI的年际变化(1985代表1985—1989年的多年平均NDVI值),从1985—2015年，镇江关流域的NDVI值增加了0.061。空间分布上，NDVI高值主要集中在流域下游区域，这与林草地土地利用类型分布相符；时间分布上，流域的NDVI值呈现总体增加的趋势，在1995年区间有减少的倾向，并且在1995年后大幅度增加。原因在于对比2000年前后的多年平均降雨量，2000年后的降雨量增加了12%,降雨量的增多促进了林草地的增长使得植被覆盖度增加；同时在2000年后该地区相继实施了“天然保护林”和“退耕还林”等一系列重大环保措施，促进了区域森林和草地资源的保护与重建；现有研究表明20世纪90年代末期，为了促进工业发展，流域内开始建设水利工程设施，流域总装机容量和水电站数量增长，水利工程的修建改变了流域内的生态格局;全球变暖使得区域内气温上升也加剧了高海拔地区的冰川消融，为该地区植被生长提供了丰富水源，一定程度上促进了植被生长。
3.1.2 径流变化趋势及突变分析
图6为镇江关流域1985—2015年历史径流曲线，图6中表明1985—2015年镇江关流域的年径流深呈总体上升变化趋势。与1985年相比径流深最低年份(2002年)总计下降了98 mm, 1985年与径流深最高年份(2012年)相比上升了257 mm。对径流数据进行M-K趋势分析得到径流数据的趋势分析特征值为1.52,表明径流呈现不显著的上升趋势。
图6 镇江关流域历史径流曲线
使用M-K突变分析的检测方法(见图7),从图上看出UF与UB曲线在2002年出现交点，因此2002年被检测为镇江关流域径流变化的可能突变点。
图7 镇江关流域径流突变分析
对比2000—2005年区间，流域内的土地利用空间分布和类型比例都有所改变，同时植被指数NDVI在2000年后大幅度增加，综合考虑下垫面因素以及径流序列可能突变点，选择2002年作为划分自然期和影响期序列的突变点。将研究年份1985—2015年分割为两个时期：1985—2002年和2003—2015年。1985—2002年称为自然期，该时期下垫面变化和气候变化影响较小，水文循环长时间保持在自然状态；2003—2015年称为影响期，该时期较为频繁的下垫面变化(“天然保护林”“退耕还林”政策，修建水利工程，城市化等)对水文过程产生较大影响。同时考虑到该地区2008年发生特大地质灾害，短期内生态环境变化较大，因此为了科学分析气候变化和土地利用变化对水文过程的具体影响方式，将研究时期分为三个阶段：1985—2002年(自然期)、2003—2008年(影响期Ⅰ)、2009—2015年(影响期Ⅱ)用于分析气候变化和下垫面变化对镇江关流域水文平衡的影响。
3.2 傅式公式参数ω的区域化
本文利用P-M公式计算流域潜在蒸散发量，得到镇江关流域1985—2015年的潜在蒸散发估算量。计算表明，1985—2015年镇江关流域潜在蒸散发量的平均值为1 229.9 mm, 最大值为1992年的1 379.5 mm, 最小值为2000年的1 079.4 mm, 估算的潜在蒸发量随时间呈逐年递减的趋势。
将估算的潜在蒸散发计算结果及镇江关流域1985—2015年的降水径流数据代入水量平衡公式和傅抱璞公式试算得到镇江关流域傅式参数ω。计算结果如表2所列。
结果表明，1985—2015年镇江关流域的参数ω的平均值为1.286,最大值为2007年的1.467,最小值为1988年的1.132。其中自然期、影响期Ⅰ、影响期Ⅱ的参数ω平均值分别为1.25、1.33、1.336三个时段的参数都有明显不同且随时间递增。
3.3 量化气候和下垫面变化对径流的影响
2003—2008年与1985—2002年相比径流增加了32.0 mm(见表3)。使用1985—2002年计算的参数ω模拟的径流量和蒸发量分别为377.0 mm和207.0 mm。
因此由于气候变化和下垫面变化导致的径流变化分别增加了21.0 mm和11.0 mm。与1985—2002年相比，气候变化和下垫面变化分别占镇江关流域流量增加的65.6%和34.4%。
2009—2015年与1985—2002年相比径流增加了68.7 mm。使用1985—2002年计算的参数ω模拟的径流量和蒸发量分别为395.6 mm和241.6 mm。因此由于气候变化和下垫面变化导致的径流变化分别增加了39.6 mm和29.1 mm。与1985—2002年相比，气候变化和下垫面变化分别占镇江关流域流量增加的57.6%和42.4%。因此得出结论，与1985—2002年相比，2003—2008年和2009—2015年时间段中气候变化都是径流增加的主要原因。
表3中，R为实测径流深(mm);ETS为使用自然期参数ω估模拟的实际蒸发量(mm);RS为使用自然期参数ω估模拟的径流深(mm);R下垫面为下垫面变化导致的径流变化(mm)。
4 结果讨论4.1 土地利用-植被覆盖变化的关系
土地利用和植被覆盖变化主要由人类活动引起，包括退耕还林、禁牧等水土保持措施。1999年，四川省率先实施了退耕还林政策，导致耕地类型逐步向林地和草地转移，耕地面积逐渐减少。从NDVI的变化上也能体现土地利用的变化，NDVI从1985年的0.383增加到2015年的0.479,植被更加茂盛导致被植被吸收的水分更多，最终被蒸发蒸腾消耗。近年来全球气候变暖加剧了高山地区的冰川融化，为植被生长提供了充足的水源，同时冰川融水流入河道导致了河道径流的增加。
4.2 气候变化和下垫面变化对径流变化的综合影响
径流变化的归因对径流序列突变点选择十分敏感，根据检验径流序列突变点和分析流域下垫面变化，选择2002年为研究时间序列突变点。同时考虑到2008年该地区发生了特大地震，地质灾害一定程度上会改变下垫面条件，为了对其进行更详细的分析，增加2008年为另一突变点。对镇江关流域进行参数ω区域化，1985—2015年镇江关流域的多年平均参数ω为1.286。同时2002年以前的参数ω变化比较平缓，在2003—2008和2009—2015年期间参数ω逐渐增加，这与土地利用和植被覆盖变化相一致，证明参数ω的变化与下垫面变化有明显因果关系。利用参数ω进行径流变化的归因分析，得出影响期Ⅰ与影响期Ⅱ气候都是影响径流变化的主要因素。通过分析气候及下垫面变化发现，流域内降雨量的增加是径流增加的直接原因，同时“退耕还林”、“禁牧”等政策的实施使得流域内的植被覆盖率提高，水土流失现象有所改善，但因流域的人口因素以及所处地理环境导致该区域存在土地利用变化不明显、城市化进展缓慢等现象，一定程度上降低了下垫面对径流的影响，使得气候因素占据主动。而2008年特大地震的发生又使下垫面的影响变得强烈，导致影响期Ⅱ下垫面变化占比有所增加。张会兰等通过建立BPCC分布式水文模型在镇江关流域研究径流对气候及覆被变化的响应，其结果为气候因素和下垫面因素对径流量的贡献率分别为84.5%和15.9%,气候是影响流域径流变化的主要原因，这与本文的结果一致，但气候与下垫面的影响占比不同，这是由于研究方法和研究尺度的差异所导致结论的不同。
5 结 论镇江关流域位于我国西南部的生态环境脆弱区，对自然环境的变化十分敏感。本文基于径流、气候、土地利用和植被覆盖资料，利用M-K分析法研究了1985—2015年的径流数据并分析相应时期的下垫面变化。以Budyko假设为理论基础，在参数ω区域化的基础上，分析了镇江关流域气候变化和下垫面变化对水文平衡的响应以及径流变化的归因，结论如下：
(1)在下垫面变化方面，由土地利用遥感影像分析可知，2000年前该地区的城市化发展比较缓慢，植被覆盖指数NDVI在2002年前后出现了大幅度上升的现象，结合径流序列的突变分析，认为2002年为镇江关流域下垫面干预流域径流变化明显的突变点。
(2)通过建立镇江关流域水热耦合模型，对傅式公式参数ω与下垫面变化的关系分析，参数ω的变化与下垫面因素的变化保持一致。
(3)通过径流变化归因分析，与自然期相比，影响期Ⅰ阶段气候变化和下垫面变化分别占镇江关流域径流增加的65.6%和34.4%,而影响期Ⅱ阶段气候变化和下垫面变化分别占流域径流增加的57.6%和42.4%,表明气候变化是镇江关流域径流增加的主要原因。另外若下垫面条件在短时间内发生突变(如地震、滑坡等地质灾害),也会剧烈影响流域的径流过程，导致生态环境的恶化。
本文使用基于Budyko假设的物理模型，考虑了气候因素、下垫面因素包括下垫面突变情况下(如地震等自然灾害)的影响，分析了气候变化和下垫面变化影响下的径流响应过程。研究中将气候条件和下垫面条件看作是两个相互独立的因素，但实际上气候与下垫面存在相互促进相互制约的关系，下一步的研究应尝试考虑森林植被与气候条件的协同机制，建立植被时空变化与气候因子序列的对应耦合关系，实现气候与下垫面条件协同变化下的径流响应分析。
水利水电技术（中英文）
水利部《水利水电技术（中英文）》杂志是中国水利水电行业的综合性技术期刊（月刊），为全国中文核心期刊，面向国内外公开发行。本刊以介绍我国水资源的开发、利用、治理、配置、节约和保护，以及水利水电工程的勘测、设计、施工、运行管理和科学研究等方面的技术经验为主，同时也报道国外的先进技术。期刊主要栏目有:水文水资源、水工建筑、工程施工、工程基础、水力学、机电技术、泥沙研究、水环境与水生态、运行管理、试验研究、工程地质、金属结构、水利经济、水利规划、防汛抗旱、建设管理、新能源、城市水利、农村水利、水土保持、水库移民、水利现代化、国际水利等。
