西藏“一江两河”流域生态系统服务变化及权衡与协同关系研究

2022-03-14李景吉王茂生卞梨交黎文婷许文来

水土保持研究 2022年2期

陈安, 李景吉,3, 王茂生, 张晗, 卞梨交, 黎文婷, 许文来,3

(1.成都理工大学生态环境学院, 成都 610059; 2.成都理工大学地球科学学院,成都 610059; 3.地质灾害防治与地质环境保护国家重点实验室, 成都理工大学, 成都 610059)

生态系统服务是指人类从自然或人工生态系统中获得的全部惠益，包括供给服务、调节服务、支持服务和文化服务等[1]。随着联合国千年生态系统评估(Millennium Ecosystem Assessment, MA, 2005)项目实施，生态系统服务和区域发展之间的关系研究已成为生态学研究的热点问题之一[2]。生态系统服务是维持人类社会发展的重要基础，其空间格局存在多样性、非均衡性，加之人类对生态系统服务管理和使用也具有选择性，使得各生态系统服务间有着复杂的相互关系，呈现出相互增益的协同或此消彼长的权衡[3]。研究人类活动影响下的生态系统服务间权衡协同关系，对推进生态文明建设和实现区域可持续发展具有重要意义。

近年来，国内外学者分别从不同研究尺度、方法等方面对生态系统服务间的权衡协同关系开展研究。在研究权衡协同关系的方法上，主要采用了统计学、空间分析、生态系统流动和情景模拟等方法[4]。包蕊等[5]采用多目标线性规划模型，优化了甲积峪流域生态系统服务间的权衡协同关系研究；陈登帅等[6]采用情景模拟法，预测了渭河流域2050年的土地利用格局，并且进行最优生态系统服务间的权衡协同关系研究；尹礼唱等[7]采用相关系数法，量化研究了“两屏三带”区域生态系统服务间权衡协同关系的时空分布。另外，诸多学者也从不同时空尺度下研究了生态系统服务之间的相互作用关系。Schrobback[8]和Dymond[9]等指出，在某些区域实施人工造林会提高土壤保持服务和固碳服务的水平，但也会降低区域的水源涵养服务；Zhou等[10]利用生产性边界 (Production Possibility Frontier, PPF)，深入研究了固碳、农业生产和水土保持3种服务之间的相互作用机制；Howe等[11]通过研究大量生态系统服务关系案例，发现权衡协同关系中的权衡关系是协同关系3倍左右。可见，基于不同研究尺度下，生态系统服务间权衡协同关系不但存在时空差异，而且各服务间的相互作用机制也较为复杂。

一江两河流域地处西藏高原腹地，其中一江是指雅鲁藏布江，两河分别指雅鲁藏布江的一级支流年楚河和拉萨河。该流域面积仅占西藏的5.48%，但人口约占总人口的1/3，经济总量也远超西藏其他区域[12]。目前，国内学者对青藏高原重大流域的研究多聚焦于生态系统服务功能分类、评估与制图，权衡与协同关系等。例如，郑德凤等[13]利用价值模型研究了三江源地区生态系统服务的价值变化及其权衡与协同关系；连喜红等[14]采用InVEST模型和全局莫兰指数研究了青海湖流域生态系统服务的空间变化特征和各服务间的相关关系；卢慧婷等[15]采用InVEST模型揭示了拉萨河流域的生态系统质量变化及其对生态系统服务的影响。然而，国内外对一江两河流域生态系统服务评估的相关研究鲜有报道，各生态系统服务之间的权衡与协同关系不清楚。因此，本文基于InVEST模型，并结合一江两河流域的地域特征，定量研究1990—2018年流域生态系统的土壤保持服务、水源涵养服务和固碳服务的时空变化特征及其权衡与协同关系，以期为流域各县推进生态文明建设、实行更具针对性的生态恢复治理措施提供科学依据。

1 研究区概况

一江两河流域位于雅鲁藏布江中游，行政区划上包括拉萨市等18个市县(图1)，地理范围为87°40′—92°37′E，28°60′—30°30′N，面积为6.67万km2。流域位于高原温带季风半干旱气候区，年均温度5～10℃，年累计降雨量200～500 mm，年均相对湿度40%～50%[16]。区内地势差异大、降雨季节性强，主要集中在5—9月，是西藏重点的生态脆弱区和治理区[17]。全区植被类型以高山草甸为主，土壤类型以沙壤土和轻壤土为主，土壤有机质含量较低。自20世纪80年代，一江两河流域开始实施大规模治沙造林工程，河谷两岸的阶地、河滩地和山地缓坡上分布大量人工林，树种以北京杨(Populus×beijingensisW. Y. Hsu)和旱柳(Salixmatsudana)为主。

图1 一江两河流域位置

2 数据与方法

2.1 数据来源

气象数据来自中国气象数据共享网 (http:∥data.cma.cn)，选取一江两河流域内及周围15个气象站点的气温、降水、风速、湿度和日照时数等数据；陆地生态系统类型数据来自中国科学院资源环境数据中心(http:∥www.resdc.cn)，该数据基于1∶10万土地利用数据，将各土地覆被类型分类合并成农田、森林、草地、水体与湿地、聚落、荒漠生态系统类型，空间分辨率为30 m；土壤数据来源于世界土壤数据库 (Harmonized World Soil Database)的中国土壤数据集(v1.2)(westdc.westgis.ac.cn)，分辨率为1 km；地形数据来自地理空间数据云平台(http:∥www.gscloud.cn)，分辨率为30 m。

2.2 研究方法

2.2.1 土壤保持本研究采用InVEST(Integrated Valuation of Environmental Services and Tradeoffs)模型3.8.2版本进行生态系统服务分析与评价。土壤保持服务利用InVEST模型中Sediment Delivery Ratio模块进行计算，该模块是在广为应用的土壤流失方程(USLE)的基础上，考虑了地块本身对上游沉积物的拦截能力，使计算结果更为准确。计算公式如下：

RKLS=R×K×LS

(1)

USLE=R×K×LS×C×P

(2)

SD=RKLS-USLE

(3)

式中:RKLS为特定地貌及气候条件下研究区的潜在土壤侵蚀量(t);USLE是考虑了植被拦截作用、实施管理工程措施后的实际土壤侵蚀量(t);SD为土壤保持量(t);R为降雨侵蚀力[MJ·mm/(hm2·h·a)],基于Wischmeier等[18]月尺度计算公式,利用克里金法插值得到;K为土壤可蚀性[t·hm2·h/(hm2·MJ·mm)],利用EPIC模型[19]计算得到;LS为坡长坡度因子;C,P为植被覆盖因子和土壤保持措施因子,参考前人[20-21]研究结果得到。

2.2.2 水源涵养水源涵养服务通过InVEST模型中Water Yield模块进行计算。该模块是基于水量平衡原理，考虑研究区不同植被、地貌和气候条件等因素，利用降水量减去实际蒸发量估算出地区产水量。相关公式如下：

(4)

(5)

(6)

(7)

式中:Yx为生态系统类型栅格x的年产水量(mm);AETx为栅格x的年实际蒸散发量(mm);Px为栅格x的年降水量(mm);PETx为栅格x的潜在蒸散量(mm);Kx为栅格x上特定植被的蒸散发系数,由前人研究[22-23]及联合国粮农组织(FAO)的作物蒸发系数指南综合确定;ET0为参考蒸散量,由Modified-Hargreaves[24]公式计算得到;AWCx为植被可利用含水量(mm),依据周文佐等[25]确定的拟合模型计算;wx为非物理参数;Z为经验参数。

2.2.3 固碳固碳服务采用InVEST模型中Carbon模块进行计算。该模块以各生态系统类型或植被类型为评价单元，通过输入对应的碳密度参数计算地区生态系统碳储量。本研究考虑4个基本碳库(表1)，即地上生物碳、地下生物碳、土壤有机碳和枯落物有机碳。计算公式如下：

C=Cabove+Cbelow+Csoil+Cdead

(8)

式中:C为总碳储量[t/(hm2·a)];Cabove为地上生物碳储量[t/(hm2·a)];Cbelow为地下生物碳储量[t/(hm2·a)];Csoil为土壤有机碳储量[t/(hm2·a)];Cdead为枯落物碳储量[t/(hm2·a)]。碳密度参数依据文献[26-29]得到。

表1 一江两河流域各生态系统的碳密度 t/hm2

2.2.4 生态系统权衡与协同关系本研究采用相关分析法，分别从流域尺度和县域尺度下探讨生态系统服务之间的权衡与协同关系，如果两类生态系统服务间相关系数为正，说明二者之间呈协同关系，如果相关系数为负，说明两者之间呈权衡关系。具体做法为：流域尺度下以1990—2018年各生态系统服务栅格数据为基础，利用ArcGIS 10.2软件设置分区采样点间隔距离为10 km，然后将采样点导出到SPSS 22.0软件中进行相关分析；县域尺度下以2000年和2010年各生态系统服务栅格数据为基础，设置分区采样点间隔距离为3 km，同理进行后续的分析。

通过对生态系统服务热点地区进行识别，可以更好地分析不同地域生态系统服务供给能力的强弱，从而对流域进行划分，对生态系统服务和人类活动提出更具针对性的对策及建议。例如，草地既可以提供碳储存能力，又可以提供土壤保持能力，但是其提供的这两种服务能力大小不同，也就导致其在单位面积上提供服务能力的不同。本研究将土壤保持服务、水源涵养服务和固碳服务超过各自平均值的地区进行叠加分析后得到一江两河流域多重生态系统服务的热点地区(表2)。

表2 热点地区分类和定义

3 结果与分析

3.1 生态系统结构时空变化特征

从1990—2018年一江两河流域土地覆被类型分布(附图1)来看，流域西北部有大片草地变为裸土地或裸岩地，主要分布在谢通门和南木林二县，但在南木林县香曲流域、谢通门县那东曲流域、雅鲁藏布江拉孜段和日喀则段的林地面积增加。流域东部墨竹工卡县有部分高植被覆盖草地退化为低植被覆盖草地，达孜区、贡嘎县、扎囊县、乃东区和桑日县疏林地和灌木林面积明显增加。近28 a流域整体的生态系统结构已经出现重大变化：相比1990年、2018年流域的农田面积减少2.09%，森林面积增加476.67%，草地面积减少31.64%，水体与湿地面积增加28.29%，聚落面积增加108.75%，荒漠面积增加87.46%。

3.2 生态系统服务时空变化特征

3.2.1 土壤保持服务 1990—2018年一江两河流域生态系统土壤保持量空间分布格局变化较小，仅在南木林县香曲流域和谢通门县那东曲流域的土壤保持量出现明显增加现象(图2)。土壤保持量较高区域主要分布在拉萨河及支流的河谷两侧，该区域林地和草地覆盖率相对较大，并呈相间分布，土壤保持能力强；较低区域分布在雅鲁藏布江上游和年楚河两侧部分地区，该地区地势相对低平，农田、建设用地等单一地类分布较多，土壤保持能力弱。统计结果显示(表3)，流域土壤保持量总体呈上升趋势，年平均值由1990年的106.65 t/(hm2·a)提高到2018年的155.03 t/(hm2·a)，增幅为45.36%。

图2 1990-2018年一江两河流域土壤保持量空间分布

表3 1990-2018年一江两河流域各项生态系统单位面积值

3.2.2 水源涵养服务一江两河流域的产水量总体呈下降趋势，年平均值由1990年的91.32 mm/(m2·a)下降到2018年的67.99 mm/(m2·a)，降幅为25.55%(表3)。其中1990—2000年和2010—2018年的产水量为增长趋势，2000—2010年为减少趋势。产水量的高值区域主要分布在流域东部的拉萨河流域，低值区域分布在年楚河流域及雅鲁藏布江拉孜段(图3)。由于一江两河流域气候干燥、光照强烈，因此根据模型计算出的产水量与各生态系统类型的蒸发能力密切相关。流域东部林地的产水量普遍较低，裸地的产水量较高，这是因为林地蒸发作用比裸地强，且林地主要分布在降水相对偏少的拉萨河下游地区，而裸地部分位于降水量偏多的拉萨河上游地区。整体上研究区的各子流域产水量基本受降水控制，降水量高的地方，产水量也高；其次受各生态系统类型的植被覆盖度和蒸发量的影响。

3.2.3 固碳服务一江两河流域碳储量在逐年上升，年平均值由1990年的129.81 t/hm2上升到2018年的143.95 t/hm2，增幅为10.9%(表3)。碳储量高值区主要分布在拉萨河流域的河谷地带和雅鲁藏布江两侧沿岸(图4)，该地区的生态系统类型以森林和草地为主，生物多样性丰富，固碳服务水平高。碳储量低值区多出现在流域西北部南木林县和谢通门县的觉母曲流域和布曲藏布流域，该地区裸地面积较大，且在1990年后出现大量草地退化成裸地现象，导致草地生态系统的固碳能力下降，区域生物多样性降低，有机物积累减少。

3.3 生态系统服务间权衡与协同关系

3.3.1 相关性分析对整个一江两河流域土壤保持服务、水源涵养服务与固碳服务3个生态系统服务进行两两之间的相关分析，结果见表4。土壤保持服务与水源涵养服务、土壤保持服务与固碳服务的相关系数显著为正(p<0.05)，并且均呈波动减小趋势，说明在流域尺度下这两对生态系统服务之间存在显著的协同关系，但协同关系不断减弱。此外，水源涵养服务与固碳服务之间的相关系数在2010年表现出显著负相关关系(p<0.01)，其余年份均为显著正相关关系(p<0.01)，这表明一江两河流域水源涵养服务与固碳服务的权衡协同关系存在时空差异性。

图3 1990-2018年一江两河流域产水量空间分布

图4 1990-2018年一江两河流域碳储量空间分布

表4 1990-2018年一江两河流域生态系统间的相关系数

在县域尺度下对3种生态系统服务进行相关分析，结果见表5—6。2000年一江两河流域各区县3种生态系统服务间相关系数基本为正，表明各服务间的主导关系为协同关系，其中土壤保持服务与水源涵养服务协同作用程度最高(18个区县中15个通过0.05水平的显著性检验)，以墨竹工卡县和拉孜县最为显著。但2010年各区县土壤保持服务和水源涵养服务的相关系数普遍减小，协同关系减弱，这与上文得出的在流域尺度下协同关系减弱的结果相一致，并在谢通门县和林周县两种服务表现出较强的协同关系。另外，2000年和2010年各区县水源涵养服务与固碳服务之间的权衡与协同关系基本相反，二者关系在南木林县和尼木县的变化最大。

表5 2000年一江两河流域县域尺度生态系统服务间的相关系数

表6 2010年一江两河流域县域尺度生态系统服务间的相关系数

3.3.2 生态系统服务热点分析对附图2中各类生态系统服务热点地区进行统计可知，1990—2018年，0类服务区主要分布在雅鲁藏布江和拉萨河河谷两侧的沙地和盐碱地附近。由此可见，沙地和盐碱地的总体生态系统服务水平较低，并且0类服务区在雅鲁藏布江上游地区有逐年增加趋势，所占比例由1990年的1.45%增加到2018年的3.44%。1类服务区面积占比由52.13%增加到52.87%，期间呈先减后增的波动特征，该类型的生态系统主要为低植被覆盖的草地，产水量高而土壤保持量和碳储量低。2类服务区主要分布在各子流域部分林地和高植被覆盖草地周边，两种生态系统服务多为土壤保持服务和固碳服务。2类服务区面积由37.63%增加到39.48%，主要增加区域为南木林县的香曲流域。3类服务区的面积占比由8.79%逐年减少至4.22%，主要分布在流域的高植被覆盖草地和灌木林区。3类服务区的时空变化特征显著，1990—2000年3类服务区主要集中于流域东部的高植被覆盖草地区；但在2000年后，流域东部大部分高植被覆盖草地退化为中植被覆盖草地，而且随着雅鲁藏布江上、下游两侧及拉萨河河谷地带的灌木林增多，3类服务区开始出现向该区域转移现象，以墨竹工卡县、谢通门县、南木林县、扎囊县和乃东区几个重点治沙造林县最明显。

4 讨论与结论

4.1 讨论

受气候变化和人类活动影响，1990—2018年一江两河流域生态系统结构已经出现重大变化，这直接导致其服务水平发生改变。本研究基于3种生态系统服务计算模型，得到一江两河流域生态系统土壤保持服务和固碳服务水平整体上升、水源涵养服务整体下降的结果。卢慧婷等[15]的研究表明，1990—2015年拉萨河流域由于灌木、草原等NDVI的上升，土壤保持量和碳储量增加，而产水量降低；Wang等[30]的研究指出，在干旱、半干旱地区实施造林不但没有增加降雨量，而且还导致地表植被蒸发作用加强，大量水资源被浪费，区域产水量下降。以上成果与研究结果基本吻合。一江两河流域属高原温带季风半干旱气候，近28 a来拉萨河等子流域河谷两岸的阶地、河滩地和山地缓坡上增加了大量人工灌木林，全区土壤保持量和碳储量增加区域均与林地增加区域保持高度一致，而作为提供水源涵养服务较高的草地发生不同程度的退化，造成了流域产水量整体下降。但本文与黄麟等[31]得到的西藏水源涵养服务1990—2010年呈上升趋势的结果有所差异，这可能是研究区范围、分析尺度及生态系统服务测算方法不同造成的。

研究一江两河流域多重生态系统服务的权衡协同关系发现，各生态系统服务间的权衡协同关系存在空间异质性，这是由生态系统服务间复杂相互作用和共同的驱动因素(地貌、气候条件或人为因素造成的土地覆被变化)引起的。1990—2018年水源涵养服务与土壤保持服务、水源涵养服务与固碳服务的相关系数波动减小，协同作用程度下降，一方面这是由于流域气候变化而导致生态系统服务水平发生改变造成的，另一方面与人类活动的影响密切相关。近28 a来一江两河流域城镇化进程加快、基础设施扩建以及不合理的旅游开发造成草地等植被面积呈减少趋势[32]，各生态系统服务水平逐年下降，生态系统服务间权衡或协同关系受到影响。但随着一江两河流域治沙造林工程实施，造林区的土壤保持服务和固碳服务协同关系得到增强。另外，水源涵养服务与固碳服务之间的相关系数在2010年表现出显著负相关，其余年份均呈显著正相关，这是由于流域内森林面积的增加提高了生态系统的固碳服务水平，同时也使部分林区水汽凝结成雨，提高生态系统的水源涵养服务，两种服务之间表现为协同关系。但是，2010年一江两河流域的降雨量少、地表蒸发量增大[33-34]，流域内植被蒸腾作用加强导致其产水量普遍降低，而这些植被又对流域生态系统固碳服务贡献较大，所以该年水源涵养服务与固碳服务表现出权衡关系。

由于一江两河流域各区县的生态系统结构和自然资源禀赋存在较大差异，因此本文计算出县域尺度下各生态系统服务数值上的相关系数普遍偏低。县域尺度下各生态系统服务间权衡协同关系与流域尺度下的不完全一致，表明流域尺度的权衡协同关系不代表次一级尺度下也存在相同关系，生态系统服务间权衡协同关系具有尺度效应。此外，本文得出森林和草地相间分布地区土壤保持量的数值相对较高，而单一生态系统类型聚集的农田或草地数值偏低，表明不同生态系统对各生态系统服务水平可能存在组合效应。

生态系统服务分为多种类型，文中仅分析了一江两河流域的3种调节服务，未来可对其他生态系统服务进行评估，为区域生态环境治理和管理提供更全面的决策依据。科学问题方面，下一步将探讨一江两河流域典型生态工程区的多种生态系统服务权衡协同关系，深入研究不同生态修复措施带来的多重生态效益。