雅鲁藏布江下游产水量时空演变及对气候和土地利用变化的响应

赖 明，陈仁升，刘玖芬，刘淑亮，吴浩然，柳晓丹

（1.中国地质大学(武汉)资源学院, 湖北 武汉 430074；2.中国科学院西北生态环境资源研究院, 甘肃 兰州 730000；3.中国地质调查局自然资源综合调查指挥中心, 北京 100055；4.自然资源要素耦合过程与效应重点实验室, 北京 100055；5.中国地质调查局烟台海岸带地质调查中心, 山东 烟台 264000；6.中国地质大学(武汉)海洋学院, 湖北 武汉 430074）

科学界将人类从生态系统获取的各种直接或间接的惠益统称为生态系统服务，包括供给、调节、文化和支持四大类服务[1]，其服务价值无法估量，更是关系着人类福祉[2]，自2005 年以来，有关生态系统服务评估的科学研究数量呈指数级增长[3-4]，其中的相当一部分是关于产水服务的研究。产水服务作为生态系统最具价值、最重要的服务功能之一[5]，是产生与水有关的生态系统服务的属性和过程的关键组成部分[6]，对流域水文状态平衡、区域水循环调节、生物多样性维持等关键生态系统功能具有重要意义[7-9]，影响着地区生物量、碳循环、农业灌溉、人口、生产生活用水、水力发电、旅游等自然和社会经济格局[10-13]。开展产水服务研究，可为区域生态保护修复、可持续发展决策及水资源合理开发、有效利用和综合管理等提供科学指导[14]。

产水过程涉及降水、蒸发、土地利用/覆被、地形、土壤性质、植被蒸腾等众多因素，这使得产水量计算十分复杂且具有很大的不确定性[15-16]。遥感和地理信息技术的发展，允许在生物地球化学过程与水文过程耦合的基础上，开发分布式物理和概念水文模型，对生态系统产水功能进行直观、定量、精细的分析和评价[15,17]，如：MIKE SHE 模型[18]、SWAT模型[19]、TOPMODEL 模型[20]、ARIES 模型[21]等，在数据充分的条件下，这些模型在各种地理和气候环境的水文研究中都能取得良好表现，但在数据匮乏或数据难以获取的情况下，简单的InVEST 模型更加适用[22]。InVEST 模型以其数据输入简单、空间表达能力强、参数设置灵活、操作便利等优点，成为目前国内外生态系统服务研究应用最多、最广泛的模型之一[23-24]。近几十年来，全球气候和土地利用/覆被变化导致产水量呈现出时空差异，InVEST模型能为探究产水量时空动态和影响贡献提供技术支持[25]。

青藏高原是全球气候变化最敏感的地区之一，许多地区的降水量波动很大[26-27]，研究表明，气候变化深刻影响了青藏高原的水文地质过程[28]，而雅鲁藏布江下游位于青藏高原东南部，是整个青藏高原降水最多的地区[29]。该地区具有水量多和地形落差大的条件，水力资源十分丰富，开发潜力巨大[30]。目前，气候和土地利用变化正在逐渐改变产水量的时空格局[31]，这将对地区水能利用产生影响。国内外众多学者已在塔里木河流域[32]、南四湖流域[33]、黄河流域[34]、伊朗[35]、印度尼西亚[36]等地区开展气候和土地利用对产水量的影响研究，但对雅鲁藏布江流域缺乏相关研究。因此，基于InVEST 模型，本文选取2000、2010 和2020 年3 期降水、蒸散发、土地利用、土壤类型等数据，分析雅鲁藏布江下游流域时空演变，探究产水量对气候和土地利用变化的响应，以期为当地水资源开发、利用和管理提供科技支撑。

1 研究区概况

雅鲁藏布江下游流域位于西藏自治区林芝市，包括米林县派镇到出境口的河段，下游河长约496 km，流域面积约为5.08 × 106hm2，地理范围为27°59′～30°59′ N，92°52′～97°05′ E。雅鲁藏布江下游有两条主要支流，位于研究区北部，分别是易贡藏布和帕隆藏布(图1)。流域内高山峡谷林立，平均海拔在3 000 m 以上，最高峰南迦巴瓦峰7 782 m，雅鲁藏布江下游围绕南迦巴瓦峰构成“U”字形的大拐弯，形成世界上最深的峡谷- 雅鲁藏布大峡谷，核心河段的平均切割深度达5 000 m。雅鲁藏布江下游段是青藏高原上最大的水汽通道，使印度洋暖湿气流在该地区形成降水带，下游巴昔卡地区的年均降水量超过4 000 mm，故整个下游地区异常湿润。雅鲁藏布江下游一泻千里，加上丰富的降水，其蕴藏的水能资源仅次于长江。

图1 研究区概况图Figure 1 Overview of the study area

2 方法及数据

2.1 InVEST 产水量模型

InVEST 模型全称“生态系统服务功能与权衡综合 评 价 模 型”(integrated valuation of ecosystem services and tradeoffs)，是由美国斯坦福大学、大自然保护协会和世界自然基金会于2007 年联合开发的开源模型，用于绘制、评估维持和满足人类生活的自然产品和服务。产水量模型是InVEST 模型的子模块之一，基于一个简单的水量平衡，即假设降水量中除开蒸散损失的部分，其他全都汇集到流域出口，则流域中每个栅格的产水量为降水量减去实际蒸发量，公式如下：

式中：Y(x)代表年产水量；AET(x)代表栅格单元x 的年实际蒸散量；P(x)代表栅格单元x的年降水量。该公式中的AET(x)/P(x)基于Budyko 水热耦合平衡假设，后来由傅抱璞[37]和Zhang 等[38]改编，公式如下：

式中：PET(x)代表潜在蒸散量；ω(x)代表自然气候-土壤特性的非物理参数。其中：

式中：Kc(lx)代表栅格单元x中特定土地利用/覆被类型的植被蒸散系数；ET0(x)代表栅格单元x的植被参考蒸散量。

ω(x)是一个经验参数，在InVEST 模型中，ω(x)采用Donohue 等[39]提出的公式表达：

式中：Z代表季节性因子，是一个经验常数，表征区域降水分布及其他水文地质特征，通常取值在1 到30 之间。AWC(x)代表土壤有效含水量(mm)，由土壤质地和有效生根深度决定，用于确定土壤为植物生长储存和提供的总水量，计算公式如下：

式中：Ds是最大根系埋藏深度，代表植物根系能在土壤中延伸到达的最大深度，Dr是植被根系深度，代表特定植被95%的根系所存在的深度。PAWC代表植被可利用含水量，是土壤中供植被生长的水量所占的比例。

2.2 情景分析法

气候变化和土地利用变化是区域产水量变化的主要影响因素，情景分析法即假定一种要素在一段时间内没有变化而其他要素变化，本文设立两个情景，情景1：保持2000 年土地利用不变，气候要素变化；情景2：保持2000 年气候要素不变，土地利用变化。通过变化贡献率可以判断不同要素对产水量变化的影响程度，公式如下：

式中：RP和RL分别代表气候变化和土地利用变化对产水量变化的贡献率；ΔP和ΔL分别代表气候变化情景和土地利用变化情景的产水量变化量。

2.3 数据来源及处理

InVEST 产水量模型使用当前或可预见的未来的土地利用和气候数据，对每个栅格单元的产水量进行空间制图，需要的基础数据具体包括：多边形流域图、年降水量(mm)、年平均参考蒸散量(mm)、最大根系埋藏深度(mm)、植被可利用含水率、土地利用数据、植被蒸散系数和植被根系深度。

多边形流域图由DEM 图经ArcGIS 软件水文分析获得，DEM 数据来源于地理空间数据云网站(http://www.gscloud.cn/)；年降水量数据由国家地球科学数据中心(http://www.geodata.cn/)下载的中国1 km 分辨率逐月降水量累加合成得到，该数据是采用CRU 发布的全球0.5°气候数据(https://crudata.uea.ac.uk/cru/data/hrg/)以及WorldClim 发布的全球高分辨率气候数据(http://www.worldclim.org/)，通过Delta空间降尺度方案在中国地区降尺度生成，并用496 个独立气象观测点数据进行验证，验证结果可信；年平均参考蒸散量由Modified-Hargreaves 方程导出(公式8)[38]，相关的气象数据来源于中国气象数据网(http://data.cma.cn/)；最大根系埋藏深度用土壤深度代替，植被可利用含水量采用Zhou 等[40]的经验公式(公式9)，利用土壤质地数据计算获得，相关土壤数据来源于世界土壤数据，其中的中国大陆区域数据由南京土壤研究所提供，数据下载自国家青藏高原科学数据中心(http://data.tpdc.ac.cn/)；土地利用数据来源于资源环境科学与数据中心(https://www.resdc.cn/)；植被蒸散系数参考InVEST 模型用户手册和联合国粮食和农业组织蒸散系数，并根据研究区土地利用和植被类型调整；植被根系深度参考联合国粮食和农业组织FAO-56 指南的标准作物系数[41]。

式中：ET0代表参考蒸散量；RA代表天文辐射量[MJ·(m2·d)-1]；Tavg代表平均日最高气温和平均日最低气温的均值(℃)；TD代表平均日最高气温和平均日最低气温的差值(℃)；P代表月平均降水量(mm)。

一是水位变化的不规则性，往往使闸引水流量无法准确控制。过闸流量受上下游水位组合及闸门开高双重因素影响，闸门不能及时跟踪调整，流量易失去控制。

式中：PAWC代表植被可利用含水量；sand代表土壤沙粒含量(%)；silt代表土壤粉粒含量(%)；clay代表土壤粘粒含量(%)；OM代表土壤有机质含量(%)。

由于本研究缺乏雅鲁藏布江下游实际径流量数据，加上该区域有相当一部分的径流量来自于冰川积雪融水，因此不适合用实际径流量对比验证Z值。据相关研究，以夏季降水为主的地区，Z值接近1；以冬季降水为主的地区，Z值接近10[42]。雅鲁藏布江流域的年内降水量分布不均，5 月-9 月的降水量可达全年降水量的80%[43]，经验证，产水量与Z值成反比，本文取最大产水量模拟值，即Z= 1。

3 结果与分析

3.1 产水量空间分异特征

2000、2010 和2020 年雅鲁藏布江下游流域的产水深度范围分别为9.6～2 543 mm、0～2 694 mm 和0～2 395 mm，各年份的产水量空间分布如图2 所示。整体上，2000-2020 年间的产水量空间分布特征相类似，呈现南高北低的格局，2 000 mm 以上的高值区分布于雅鲁藏布江出境口至雅鲁藏布江大峡谷的拐弯处，500 mm 以下的低值区分布于北部帕隆藏布和易贡藏布的流域范围内，500～1 000 mm的次低值区广泛分布于研究区北部和中部；从时间上看，2000 和2010 年的产水量高值区分布基本一致，而2010 年的低值区面积明显增加，主要分布于易贡藏布流域，2020 年的低值区面积进一步增加，分布范围扩大至研究区整个北部区域，并且高值区分布面积减少，向南部集中。

图2 2000、2010 和2020 年雅鲁藏布江下游流域产水量空间分布图Figure 2 Spatial distribution of water yield in the lower reaches of the Yarlung Zangbo River in 2000, 2010 and 2020

垂向上，雅鲁藏布江下游流经雅鲁藏布江大峡谷后，海拔由高到低急剧下降，水面落差超过2 700 m，“天上之水”实至名归；整个流域范围内，海拔从南迦巴瓦峰超过7 700 m 下降到雅鲁藏布江出境处低于200 m，海拔高差超过7 000 m，高山深壑纵横交错。海拔的巨大落差使产水量随海拔高度变化而呈现明显垂向差异，各个年份的平均产水深度随海拔升高的变化趋势完全一致，均是先急剧减少后缓慢增加，在海拔低于500 m 的区域内，各年份的平均产水深度均超过2 000 mm，最高值为2 341 mm，2000 和2010 年在海拔4 001～5 000 m 时，平均产水深度最低，2020 年的最低值出现在海拔3 501～4 000 m 的区域内，最低值为478 mm (图3)；在面积和产水量方面，以2020 年为例，海拔1 500 m 以下和3 000 m 以下的区域面积分别约占研究区总面积的11%和29%，而产水量占研究区总产水量的比例分别达到28%和50%，由此可知，低海拔区的面积集中且对产水量的贡献大，高海拔区的面积广阔，但对产水量的贡献相对不足(图3)。

图3 研究区各海拔区间的产水量情况Figure 3 Water yield at different altitudes in the study area

3.2 产水量时间变化分析

产水量模型模拟结果为2000、2010 和2020 年雅鲁藏布江下游流域的平均产水深度分别为809、772 和709 mm，产水量分别为411 × 108、392 × 108和360 × 108m3(产水量等于产 水深度与面积的乘积)。产水量总体呈下降趋势，且下降幅度呈扩大趋势，20 年间平均产水深度下降100 mm (产水量下降51 × 108m3)，其中，2000-2010 年平均产水深度下降37 mm (产 水 量下降19 × 108m3)，2010-2020 年平均产水深度下降63 mm (产水量下降32 × 108m3)。

将不同年份产水量分布图相减，获取2000-2010、2010-2020 和2000-2020 年的产水量时空变化图(图4)，其中，将2000-2010 年的产水量增减区间划分为3 类，即小于-70 mm (明显减少)、-70～70 mm(少量变动)和大于70 mm (明显增加)，由于2010-2020 年间的产水量下降幅度较大，因此2010-2020和2000-2020 年的时空变化图增加小于-200 mm(大量减少)的变化区间。各年间的产水量时空变化呈现不同的分布特征。2000-2010 年，产水量明显减少区域面积(面积占比约为29%)大于明显增加区域面积(8%)，产水量明显减少区主要分布于研究区西北部(易贡藏布流域)，明显增加区则集中于研究区南部(雅鲁藏布江末端)；2010-2020年，产水量明显减少区域面积进一步扩大(39%)，分布状况与前10 年相反，产水量明显减少区主要分布于研究区东部(帕隆藏布流域)和南部，并且南部区域的产水量大量减少，最大减少量超过1 100 mm，明显增加区则分散分布于研究区北部；2000-2020年总体上，产水量明显减少区域面积占比达到74%，遍布于整个研究区，而少量变动区和明显增加区分布范围较小，呈分散分布状。

图4 2000-2010 年、2010-2020 年和2000-2020 年产水深度时空变化图Figure 4 Temporal and spatial variation of water yield in 2000-2010, 2010-2020, and 2000-2020

3.3 产水量对气候和土地利用变化的响应

表1 不同情景下平均产水深度Table 1 Average water yield depth under different scenarios mm

在以2000 年土地利用为基准而气候要素数据改变的情景下(情景1)，2010 和2020 年的平均产水深度分别较2000 年的真实情景减少了37 和120 mm，相应产水量减少约19 × 108和61 × 108m3，约占2000年产水量的5%和15%。模拟结果表明，在土地利用不变而气候要素改变的情况下，气候要素对研究区产水量具有较大的影响。

在以2000 年气候要素为基准而土地利用数据改变的情景下(情景2)，2010 年的平均产水深度较2000 年的真实情景没有明显变化，在产水量上仅增加约1.1 × 105m3，占比低于0.01%，可以忽略不计。2020 年的平均产水深度较2000 年的真实情景增加了16 mm，产水量增加约8 × 108m3，约占2000 年产水量的2%，如图5 所示，在情景2 下，2020 年的产水量增加部分主要来自于林地和草地，减少部分主要来自于水体和未利用地；其中，较2000 年，2020年林地和草地的面积占比分别增加约7.3%和11.6%，相应的产水量分别增加约22 × 108和41 ×108m3，水体和未利用地的面积占比分别减少约6.9%和11.9%，相应产水量分别减少约16 × 108和37 × 108m3，4 类土地利用类型的面积增减相抵消，而产水量增加了约10 × 108m3。模拟结果表明，相较于情景1，在气候要素不变而土地利用改变的情况下，林地和草地的扩张使得产水量明显增加，水体和未利用地的缩减致使产水量明显减少，但总体上，土地利用变化对研究区产水量影响较小。

图5 土地利用变化情景下(情景2)各地类产水量和面积占比变化Figure 5 Change of water yield and area proportion of different types under the land use change scenario (scenario 2)

在2000-2010 年间，气候要素变化和土地利用变化对研究区产水量变化的贡献率分别为-100%和0，2000-2020 年间，气候要素变化和土地利用变化对研究区产水量变化的贡献率分别为-115%和15%，说明雅鲁藏布江下游流域产水量的主要影响因素是气候要素，而土地利用的影响贡献相对不足。另外，如表1 所列，在不同年份土地利用和气候要素的组合情景中，产水量最高的是2000 年气候要素和2020 年土地利用的组合，产水量最低的是2020 年气候要素分别与2000 和2010 年土地利用的组合，通过情景组合能够初步划定未来产水量的变化范围。

4 讨论

由于缺乏数据和融水干扰，模拟结果未能得到验证，因此本文利用雅鲁藏布江中上游多年产水量模拟结果验证模型适用性和参数准确性。采用相同的数据和方法，及Z取值为1 的条件下，获取雅鲁藏布江中上游流域的逐年产水量(图6)，2000-2010年的多年平均产水量约为470 × 108m3，据相关文献，相近时间段内的雅鲁藏布江中上游流域多年平均水文站实测径流量约为643 × 108m3[44]，近20 年的冰川积雪融化对雅鲁藏布江中上游流域实际径流的平均贡献率约为26%[45]，除去冰川积雪融化的多年平均实际径流量约为476 × 108m3，与多年平均模拟产水量相差6 × 108m3，仅约1%，说明当Z=1 时，模型能够很好模拟雅鲁藏布江流域的产水量。

图6 2000-2010 年雅鲁藏布江中上游流域逐年产水量Figure 6 Annual water yield in the middle and upper reaches of Yarlung Zangbo River from 2000 to 2010

产水量时空变化与降水量和蒸散发的时空差异密切相关，降水量与产水量呈正相关关系，蒸散发与产水量呈负相关关系，土地利用则通过影响实际蒸散发从而影响产水量[46-47]。研究区产水量由南向北迅速递减，并且平均产水深度相差巨大，这是由于研究区南部受印度洋气流影响，该区域降水量十分丰富，但气流在北上途中受到高海拔山脉阻挡，研究区整个北部降水量很低，从而使得研究区产水量北部与南端形成鲜明对比；经验证，蒸散发的年际变化很小，而2000-2020 年研究区降水量呈整体下降趋势(图7)，这与产水量的整体趋势相同，说明降水量是产水量年际变化的主要影响因素，并排除了2000、2010 和2020 年属于极端年份的可能性；土地利用对产水量的影响主要是通过土地利用类型变化来控制实际蒸散发变化，但研究区土地利用在2000-2020 年变化较小，尤其是2000-2010 年，且土地利用变化引起的产水量增减基本相抵消，因此，土地利用变化的影响相对不足。总之，研究区产水量时空演变与降水量时空变化的趋势一致，土地利用变化的影响基本被降水量稀释或掩盖。

图7 2000-2020 年雅鲁藏布江下游流域平均降水量Figure 7 Average precipitation in the lower reaches of Yarlung Zangbo River Basin from 2000 to 2020

InVEST 模型的产水量模块基于简单的水平衡模型，虽然能够对各种地区的产水量进行很好地模拟，但也存在一定的局限性和不确定性[48-49]，地形海拔是水资源变化的一个重要影响因素，尤其是具有复杂地形的山区，但产水量模块对此没有考虑[50]，本文研究区地形极其复杂，海拔高差巨大，对产水量分布产生重要影响，虽然本研究有做描述，但还有待完善。在有融水补给径流的区域，模拟产水量不能代表真实径流量，而是需要去掉融水部分。植被蒸散系数和植被根系深度参考相应文献和InVEST模型使用手册，并根据研究区特征进行调整，存在一定不确定性，另外，数据精度、参数确定和校验、参考蒸散发和植被可利用含水量计算公式的选取，这些都对最终评估结果精度产生影响，但总体上不会改变研究区产水量格局。

5 结论

本文基于InVEST 模型产水量模块计算了2000、2010 和2020 年雅鲁藏布江下游流域的产水量，分析研究区产水量的时空演变特征，并采用情景分析法探究气候变化和土地利用变化对产水量的影响，得到以下结论：

1)雅鲁藏布江下游流域各年产水量分布特征相似，研究区内产水量差别大，平均产水深度最高值超过2 500 mm，2 000 mm 以上高值区主要分布于研究区南端，面积较小，1 000 m 以下的低值区和次低值区主要分布于研究区北部，面积广阔。垂向上，平均产水深度随海拔升高而迅速递减，海拔3 000 m以下的区域以29%的面积占比拥有了全区50%的产水量，低海拔区表现为面积小而产水量高。

2)雅鲁藏布江下游流域2000、2010 和2020 年的产水量分别为411 × 108、392 × 108和360 × 108m3，呈整体下降趋势，且下降幅度扩大。2000-2010 年间，研究区西北部产水量明显减少，南部地区产水量明显增加，2010-2020 年间与前10 年相反，研究区南部地区产水量明显减少，且减少量很大，明显增加的区域分布于北部的少数地区，2000-2020 年产水量在整个研究区呈明显减少趋势。

3) 2000-2020 年，气候要素变化(情景1)和土地利用变化(情景2)对雅鲁藏布江下游流域的产水量变化的贡献率分别为-115%和15%，其中，气候要素以降水为主导。表明研究区产水量变化的主要影响因素是降水量变化，而土地利用变化的影响较小。