生态系统服务是指生态系统与生态过程所形成及所维持的人类赖以生存的自然环境条件与效用, 是人类社会发展的物质基础。正因此, 生态系统服务已成为国际生态学、经济学和资源环境学等多学科研究的前沿和热点^[1—4]。水源涵养服务是重要的供给服务, 其内涵已由供给水源, 拓展为拦蓄降水、水源供给、调节径流、洪水调蓄、水质净化等多种服务组成的综合体, 其变化不仅直接影响区域水文、植被生长与发育, 还制约着区域社会经济的发展, 是区域生态与环境状况的重要指示器^[5—8]。正因此, 国内外对生态系统水源涵养进行了大量研究^{[1, 9—11]}。

目前为止, 中国生态系统水源涵养研究主要集中在水源涵养内涵、评估方法、重要性评估等方面。关于生态系统水源涵养概念, 龚诗涵等认为早期水源涵养指生态系统对流域径流量的影响^[6]。后来生态系统拦蓄降雨、土壤持水逐渐受到重视^[11], 之后水源涵养概念扩展到生态系统结构、质量与格局在各水文过程中的作用与贡献^{[8, 12]}。生态系统水源涵养评估或计量方法因水源涵养内涵的不同而不同, 常用评估方法有土壤蓄水法、综合蓄水法、水量平衡法、综合指数法等, 如姜志林详细介绍了土壤蓄水法的评估方法与影响因素^[13]；莫菲等采用综合蓄水法研究了东灵山不同森林水源涵养功能, 并用多因子综合评价法评估了不同森林水源涵养功能高低^[14]。刘业轩等采用水量平衡法和综合蓄水法对比研究了福建森林水源涵功能, 发现水量平衡法侧重气象、地形和土地覆被因素, 而综合蓄水法侧重林分条件^[15]。程唱等利用综合指数法对比研究了干旱半干旱不同人工林水源涵养能力, 发现白桦×青杨混交林水源涵养能力最强, 青杨林最弱^[16]。总之, 不同评估方法对应着不同的水源涵养内涵与研究目的, 致使不同方法有不同的数据需求与适用范围。土壤蓄水法、综合蓄水法、综合指数法等因评估精度受研究区实测数据数量的制约, 致使它们主要适用于样地等小尺度研究, 而水量平衡法因是水源涵养机理模型且操作简单、数据可达性强而成为目前最常用的方法。正因此, 本研究采用水量平衡法动态拟合1990—2018年中国生态系统水源涵养量。关于水源涵养演变趋势研究, 目前较多研究几期生态系统水源涵养量来研究区域整体或不同区域水源涵养服务演变规律, 如薛健等通过研究1980—1990年、1990—2000年、2000—2010年和2010—2017年四时段祁连山平均水源涵养量, 揭示该区域不同子流域、不同植被水源涵养服务变化特征^[17]。张福平通过四期数据研究了黑河流域上游水源涵养演变特征^[18]。只有极少数长时间序列动态研究区域生态系统水源涵养演变特征的报道, 如吕乐婷等1980—2016年长时间序列动态研究了三江源国家公园整体水源涵养时空演变特征^[19]；尹云鹤等动态研究了气候变化长时间序列对黄河源区水源涵养的影响^[20]。中国依然缺少长时间序列栅格尺度水源涵养演变规律研究。

中国水资源贫乏, 时空分布不均, 水土资源组合不平衡, 且水资源年际变化大^[21], 探明中国生态系统水源涵养功能时空格局与演变规律, 不仅有助于科学认识和合理保护中国重要生态系统水源涵养功能区, 还可为中国生态系统水源涵养调控与生态系统管理提供数据支撑。正因此, 本研究利用水量平衡模型, 探明中国生态系统水源涵养空间格局, 从全国、分区和栅格尺度揭示中国长时间序列生态系统水源涵养演变特征, 以期为中国生态系统水源涵养保护与建设、生态补偿等提供数据。

1 研究方法 1.1 水源涵养核算

单位面积水源涵养量(即水源涵养能力)采用水量平衡方程核算。水量平衡原理是指在一定的时空内, 水分的运动保持着物质守恒定理, 或输入的水量和输出的水量之间的差额等于系统内蓄水的变化量。

式中, Q_{water conservation}为单位面积水源涵养量, mm m^-2 a^-1；P为单位面积降雨量, mm m^-2 a^-1；R为单位面积地表径流量, mm m^-2 a^-1；ET为单位面积年蒸散发量, mm m^-2 a^-1。在没有人为干扰条件下, 降水量减去蒸散量等于产水量, 故在此先利用INVEST软件产水模块核算出产水量, 再利用各类生态系统平均地表径流系数、土地覆被和年降水量数据获得单位面积地表径流量, 进而获得单位面积水源涵养量。具体核算方法详见文献^{[6, 22—23]}。本研究利用CERN台站水文监测数据对核算结果进行验证, 结果表明, 模型拟合结果与台站野外监测结果极显著线性相关, 说明本研究模型拟合结果较科学合理。图 1是2010年CERN台站实测值和拟合值的数量关系。

1.2 数据来源

本研究用到的主要数据包括1990、1995、2000、2010和2015年土地覆被、历年气象、高程、土壤、叶面积指数及文献数据等, 其中五期土地覆被数据和叶面积指数数据来自中国科学院遥感与数字地球研究所；气象数据来自中国气象数据网；潜在蒸散数据利用彭曼公式计算后克里格插值获得；蒸散系数利用LAI指数计算获得；土壤可利用水量利用土壤数据中的质地数据计算获得；根系深度和平均径流系数则利用文献数据获得。本研究主要数据来源与处理方法详见表 1.

1.3 数据分析方法 1.3.1 全国和分区尺度趋势分析

全国和分区尺度水源涵养趋势分析均利用SPSS软件进行回归分析再在Excel中制图。首先在ArcMap中利用Zonal Statistics as Table命令获得历年全国和分区单位面积水源涵养量及水源涵养总量, 再在SPSS统计软件中进行线性回归分析, 最后在Excel中制图^[29]。

1.3.2 栅格尺度趋势分析

1990—2018年中国栅格尺度水源涵养服务趋势分析同样也是进行线性回归分析, 首先在ArcMap中将栅格数据转化成txt数据, 再将txt数据导入MATLAB中进行线性回归分析, 获得各栅格的Slope、F值及其显著性特征, 最后将分析结果转化成栅格数据在ArcMap中制图。

2 结果与分析 2.1 水源涵养功能空间格局与演变 2.1.1 整体空间格局

1990—2018年中国生态系统水源涵养总量及其能力平均分别为15396.2亿m³/a和164.0 mm m^-2 a^-1。中国生态系统水源涵养功能时空异质性显著, 空间格局整体呈现东高西低、南高北低的分布态势, 其中大于1000 mm m^-2 a^-1的髙值区主要分布于热带-亚热带季风区的长江中下游地区, 600—1000 mm m^-2 a^-1的次髙值区主要分布于长江中下游和华南地区, 小于50 mm m^-2 a^-1的低值区广泛分布于西北、青藏高原、内蒙古、华北平原和东北平原地区(图 2)。受年际波动的影响, 中国生态系统水源涵养髙值区主要在长江中下游和华南地区移动, 而东北和青藏高原西南部地区是中国生态系统水源涵养波动较大的区域。1990—2018年6期水源涵养空间分布详见图 2。

深入分析发现, 不同植被多年平均水源涵养能力以常绿阔叶林最高, 约607.7 mm m^-2 a^-1, 园地次之, 约490.1 mm m^-2 a^-1, 随后依次为常绿针叶林、针阔混交林、灌丛、农田、落叶阔叶林、落叶针叶林和草地(图 3)。水源涵养服务在植被间的变化特征是植被分布、年降水量和人类活动等综合作用的结果。

2.1.2 随海拔变化规律

中国地势存在明显的阶梯分布：第一阶梯是青藏高原, 平均海拔≥4000 m；第二阶梯由青藏高原外缘至大兴安岭、太行山、巫山和雪峰山之间, 平均海拔≥1000 m, 主要由广阔的高原和大盆地组成；第三阶梯为中国东部广阔的平原和山地丘陵, 海拔范围主要为＜1000 m地区^[30]。根据上述地貌海拔特征, 将中国高程数据分为3级, 提取不同高程梯度内多年平均单位面积水源涵养量(表 2)。低海拔地区主要为东北平原、黄淮海平原、长江中下游平原以及江南广大地区的丘陵区, 降雨量丰沛, 植被类型丰富且覆盖度高, 该区生态系统水源涵养能力服务均值最高。高海拔地区主要为青藏高寒区, 降雨少, 植被类型为高山灌丛、高寒草甸、高寒草原, 植被覆盖度低, 该区生态系统水源涵养能力均值最低。中国生态系统水源涵养能力随着高程的升高总体呈现降低的趋势(图 4), 但在不同高程阶梯内的变化特征不尽相同, 低海拔地区水源涵养能力随着海拔的升高呈较强的下降趋势, 且降幅较大；中海拔地区生态系统水源涵养能力随海拔升高呈阶梯式波动, 海拔1000—2800 m范围为水源涵养的第一阶梯区, 区内水源涵养能力随着海拔升高而增加, 海拔2800—4000 m范围为水源涵养的第二阶梯, 区内的水源涵养能力也随着海拔升高, 但总体上第二阶梯内的水源涵养能力低于第一阶梯区, 这与山区降水量随地势抬升而增加有关, 第一阶梯主要分布于热带、亚热带区域, 而第二阶梯主要分布于云贵高原与青藏高原的过渡带上, 这里主要受西南季风的影响, 由于地处内陆, 年降水量随海拔升高而增大的幅度小于东部热带-亚热带地区。高海拔区生态系统水源涵养能力随着海拔升高急剧降低。

2.1.3 随坡度变化规律

根据《森林资源规划设计调查主要技术规定》(2003年4月), 将中国坡度数据分为6级, 提取不同坡度梯度内的多年平均水源涵养量(表 3)。结果表明, 中国水源涵养能力以平坡地区最小, 以险坡地区最大。平坡地区主要分布于年降水量较小的西北和年降水量较大但主要种植农作物的华北平原、东北平原、成都平原和长江中下游平原等地, 而险坡主要分布于降水量较大的南方山地丘陵区, 如福建武夷山地区坡陡降水量又大, 此外, 云岭、南岭、四川盆周山地、秦岭、长白山等降水量较大的山区, 也是中国陡坡险坡主要分布地区。同时将坡度与多年平均单位面积水源涵养量进行叠加分析, 制作散点图, 进而分析水源涵养能力与坡度的关系(图 5)。水源涵养能力与坡度关系因不同坡度梯度而不相同, 其中0°—5°坡度范围内生态系统水源涵养量随着坡度的升高而增大, 且变化幅度较大；5°—35°坡度范围内生态系统水源涵养能力变化较小, 相对稳定, 生态系统水源涵养能力平均为200 mm m^-2 a^-1；＞35°坡度范围水源涵养能力呈散点状分布, 且值域波动范围较大, 生态系统水源涵养能力最小值和最大值均位于该坡度范围内, 这是因为该坡度梯度内既有降水量大的南方山地丘陵区, 又有降水量较少的西北干旱区。

2.2 中国水源涵养服务时间尺度演变趋势 2.2.1 全国尺度

1990—2018年, 中国生态系统水源涵养总量在1.35×10¹²—1.92×10¹² m³/a波动, 平均为1.54×10¹² m³/a, 水源涵养能力变化范围为116.3—198.8 mm m^-2 a^-1, 平均为164.0 mm m^-2 a^-1。趋势分析结果表明, 中国生态系统水源涵养总量和水源涵养能力均呈显著线性降低态势(图 6)。

2.2.2 分区尺度

1990—2018年, 中国生态系统水源涵养服务平均为1.20×10¹² m³/a, 以热带-亚热带季风区占主导, 约占中国水源涵养总量的77.9%, 青藏高寒区次之, 约占总量的10.5%, 随后是温带季风区, 温带大陆气候区最低, 仅约占水源涵养总量的2.5%。与水源涵养总量不同, 水源涵养能力从大到小排序依次为热带-亚热带季风区＞温带季风区＞青藏高寒区＞温带大陆气候区, 平均分别为483.7 mm m^-2 a^-1、82.5 mm m^-2 a^-1、63.9 mm m^-2 a^-1和14.2 mm m^-2 a^-1。回归分析表明, 除温带大陆气候区1990—2018年生态系统水源涵养服务呈显著降低外, 其余分区水源涵养服务均呈无显著规律的波动变化。各区生态系统水源涵养服务演变特征详见图 7。

2.2.3 栅格尺度

1990—2018年, 中国绝大部分栅格生态系统水源涵养服务变化不显著, 平均约占中国陆地总面积的84.0%, 仅约16.0%的栅格水源涵养发生显著变化, 其中又以极显著降低和显著降低占主导, 分别约占陆地总面积的7.6%和5.3%, 二者之和约占显著变化栅格总数的80.7%, 主要分布于藏东南、川西滇北、滇中与滇西、冀东北与浑善达克沙地周边、天山东麓和台湾省。极显著升高或显著升高区域约占陆地总面积的1.4%和1.7%, 它们主要分布青藏高寒区中东部地区、秦岭、三江平原和大兴安岭北部寒温带地区。1990—2018年中国栅格尺度生态系统水源涵养演变特征详见图 8和表 4。

3 结论与讨论 3.1 讨论 3.1.1 水源涵养功能与格局

尽管中国生态系统水源涵养研究较多, 但采用水量平衡法动态研究全国生态系统水源涵养服务的较少。龚诗涵等研究发现, 2010年中国林地、灌丛、草地、园地和湿地水源涵养总量的12224.33亿m³, 其中水源涵养量较高的区域主要集中在武夷山脉、南岭、武陵山区、大巴山区, 四川盆地, 其次为云贵高原量较高的区域, 再次为长白山与大小兴安岭地带, 蒙古高原、青藏高原以及新疆西北部的水源涵养量较低^[6]。Zhang等利用250 m分辨率土地覆被数据研究发现, 2010年中国生态系统水源涵养空间格局为东高西低、南高北低, 水源涵养总量为12932.2亿m³, 平均为136.1 mm m^-2 a^-1[31]。本研究发现, 1990—2018年中国生态系统水源涵养分布格局为南高北低、东高西低, 水源涵养总量平均为15396.2亿m³/a, 水源涵养能力平均为164.0 mm m^-2 a^-1。因此, 水源涵养空间格局与前人研究完全一致, 仅水源涵养功能大小略有差异。造成生态系统水源涵养服务差异主要有以下原因：(1)数据分辨率所致, 本研究土地覆被数据的分辨率为1km, 而龚诗涵等^[6]土地覆被数据分辨率为90 m, Zhang等^[31]用到的土地覆被数据分辨率为250 m, 即本研究分辨率均低于龚诗涵等和Zhang等的数据分辨率, 必将导致本研究植被系统面积小于Zhang等和龚诗涵等植被系统面积, 而植被系统面积的增加必导致植被蒸散量的增加, 进而导致水源涵养量的降低, 这可能是本研究水源涵功能高于龚诗涵等和Zhang等的重要原因；(2)研究时段差异所致。龚诗涵等^[6]和zhang等^[31]的研究时段均为2000—2010年, 而本研究的研究时段为1990—2018年。由于水源涵养数值与年降水量极显著相关, 而中国降水量年际波动较大。(3)研究对象差异, 龚诗涵等仅研究了一些植被系统, 而本研究研究了所有生态系统。因此, 剔除研究时段、研究对象和数据精度差异的影响, 本研究结果与前人研究结果相符。就不同植被水源涵养量而言, 本研究结果与龚诗涵等^[6]研究结果基本一致, 如都是常绿阔叶林＞常绿针叶林＞针阔混交林＞落叶看阔叶林＞落叶针叶林, 且不同森林类型水源涵养能力也与之结果较接近, 园地水源涵养能力(490.1 mm m^-2 a^-1)略大于龚诗涵等研究结果(467.9 mm m^-2 a^-1)。常绿阔叶林、常绿针叶林、落叶阔叶林和混交林水源涵养能力分别为607.7 mm m^-2 a^-1、439.4 mm m^-2 a^-1、147.2 mm m^-2 a^-1和380.3 mm m^-2 a^-1, 森林水源涵养能力小于刘业轩等福建研究结果(分别约为660 mm m^-2 a^-1、610 mm m^-2 a^-1、600 mm m^-2 a^-1和560 mm m^-2 a^-1)^[15], 却大于余新晓等北京森林平均水源涵养量(75 mm m^-2 a^-1)^[32]、薛健等祁连山森林平均水源涵养量(102.1 mm m^-2 a^-1)等研究结果^[17], 这主要由研究区域年降水量差异所致。

3.1.2 水源涵养演变趋势

关于中国水源涵养功能演变的研究, 主要在区域/流域尺度开展几期水源涵养变化研究, 缺乏栅格尺度水源涵养功能长时间系列的动态演变研究。本研究发现, 1990—2018年, 约16.0%中国陆地栅格水源涵养功能发生显著变化, 其中主要为极显著或显著降低, 约占显著变化区域面积的80.7%, 主要分布于藏东南、川西滇北、滇中、滇西、冀东北与浑善达克沙地周边、天山东麓和台湾省；仅约19.3%的显著变化区域为极显著或显著升高, 它们主要零星分布青藏高寒区中东部地区、河套平原、秦岭和东北平原等地区。中国已有不少区域/流域水源涵养下降的报道, 如1960—2020年川西亚高山区域水源涵养功能在以2.83 mm/a的速率下降^[33]；2001—2015年祁连山上游黑河流域水源涵养服务以7.6 mm/a的速度递减^[18]；1980—2005年三江源地区水源涵养功能整体呈下降态势^[34]；2000—2010年陕北黄土高原大部分区域水源涵养服务有不同程度下降, 其中以中部地区下降的尤为显著^[25]。

本研究发现, 川西滇北地区水源涵养功能极显著降低, 祁连山黑河上游地区大部分区域有不同程度下降黄土高原和三江源地区大部分区域水源涵养有不同程度增加, 这与川西亚高山区域和祁连山上游黑河流域水源涵养呈下降态势的报道空间上完全吻合, 而与三江源区和陕北黄土高原水源涵养降低的报告不符, 其原因主要有：(1)研究时段差异所致, 本研究时段为1990—2018年, 三江源区和陕北黄土高原区研究时段分别为1980—2005年和2000—2010年。研究时段的不同必然导致主要受降雨影响的水源涵养功能变化的差异；(2)数据差异, 本研究用到的是中科院遥感所解译的1km土地覆被栅格数据, 陕北黄土高原区域研究土地利用分辨率为30 m, 三江源区研究用的是中国科学院地理科学与资源研究所刘纪元团队解译的1km土地利用数据, 数据源和数据精度均有差异；(3)研究方法差异, 本研究方法是基于线性拟合做出的趋势判定, 而三江源区和陕北黄土高原区的研究结果是根据不同时期差值判定获得。

中国也有不少区域水源涵养功能增强的报道, 如张彪等发现2000—2018年京津风沙源治理工程区的水源涵养服务有所增强^[35]；薛健等发现1980—2017年祁连山水源涵养呈上升趋势^[17]；宁亚洲等得出2000—2015年秦岭水源涵养服务以0.26 mm/a的速度增加^[36]。本研究发现, 秦岭地区和祁连山南部生态系统水源涵养有不同程度的增加, 且不少栅格增加达极显著或显著水平, 而浑善达克山地周边极显著或显著降低。本研究结果与李亚洲等和薛健等研究结果较一致, 而与张彪等研究结果不同, 研究结果差异主要由以下三方面原因所致：(1)研究时段差异, 本研究和张彪等研究时段分别为1990—2018年和2000—2018年, 相差了10年数据；(2)数据差异, 本研究1 km土地覆被数据来自中国科学院遥感所遥感解译, 而张彪等用的是中国科学院地理科学与资源研究所资源环境科学数据中心30 m土地覆被数据, 数据源及数据精度均存在差异；(3)蒸散发的计量方法。本研究采用彭曼公式计算潜在蒸散量, 再根据植被蒸散系数计算实际蒸散量, 而张彪等根据评估区温度计算最大蒸散量, 再引入土地覆被调整系数获得蒸散发, 实际蒸散发计量方法不同。

3.2 结论

基于水量平衡模型, 从全国、分区和栅格三个尺度对中国1990—2018年生态系统水源涵养格局与演变展开研究, 得出以下主要结论：

(1) 中国生态系统水源涵养整体呈南高北低、东高西低的分布格局。水源涵养总量及水源涵养能力均呈显著下降态势, 平均分别为15396.2亿m³/a和164.0 mm m^-2 a^-1。水源涵养能力在植被类型、海拔梯度和坡度梯度间差异显著, 分别以常绿阔叶林、＜1000m的低海拔和＞45°险坡最高。

(2) 中国生态系统水源涵及其能力均以热带-亚热带季风区最高, 温带大陆气候区最低, 分别约提供了中国水源涵养总量的77.9%和2.5%。除温带大陆气候区水源涵养服务显著降低外, 其余分区水源涵养服务在近30年间变化不显著。

(3) 84.0%的中国陆地栅格水源涵养变化不显著, 仅约16.0%陆地栅格水源涵养变化显著, 主要为极显著降低和显著降低类型, 分别约显著变化栅格总数的47.6%和33.1%, 主要分布于藏东南、川西滇北、滇中西、冀东北、浑善达克沙地、天山东麓和台湾省；极显著或显著增加的栅格约占显著变化总量的19.3%, 主要分布于青藏高寒区东北部、秦岭、三江平原和大兴安岭寒温带地区。本研究将为中国生态系统水源涵养服务保护与调控区域遴选提供支撑。