青藏高原保护区植被与生态系统服务功能变化的一致性分析

2018-02-13程琳琳刘焱序

水土保持通报 2018年6期

程琳琳, 刘华, 刘焱序

〔1.中国矿业大学(北京) 地球科学与测绘工程学院, 北京 100083;2.北京师范大学地理科学学部地表过程与资源生态国家重点实验室, 北京100875〕

植被是陆地生态系统的主体，是调节全球和区域物质和能量平衡的关键要素，也是生态系统提供服务的重要载体[1-2]。在近年研究中，遥感植被指数已被广泛用于陆地植被时空变化的刻画，其中归一化差值植被指数(NDVI， the normalized difference vegetation index)是衡量植被变化的最常用遥感监测参数[3-6]。不同数据源的NDVI产品交互验证表明，近30 a来中国植被变化普遍呈现增加趋势，但不同数据产品和趋势提取方法所呈现的区域植被变化状态略有差异[7-9]。然而，NDVI在高覆盖植被区的取值往往存在饱和现象，导致不利于监测高密度植被的时空演化。为解决这一问题，增强型植被指数(EVI， enhanced vegetation index)被逐步应用于植被变化研究中，并已成为监测大尺度植被变化特征的重要工具[10-12]。

青藏高原被誉为亚洲水塔，是世界多条大江大河的发源地。青藏高原生态系统的碳固定、土壤保持、水源涵养、栖息地保护等生态系统服务功能对中国和亚洲人类生活具有重要的影响。作为评估生态系统服务功能的关键指标，青藏高原植被变化一直备受关注[13-15]。1982—2011年，青藏高原生态系统碳汇的增加占全国增加碳汇的10%左右，但也有部分地区出现了植被退化[16]。对应了生态十年评估，是前人研究往往采用的时段，并且2000年以来，中分辨率成像光谱仪(MODIS， moderate-resolution imaging spectroradiometer)系列产品的发布为更高精度的刻画植被变化、评估生态系统服务功能提供了新的数据支持。在此背景下，本研究聚焦青藏高原各类保护区，试图回答以下两项科学问题： ①基于MODIS EVI产品的植被变化制图是否反映了青藏高原植被整体好转局部退化的特征； ②青藏高原保护区内2000—2010年的生态系统服务变化是否与同期的植被变化一致。研究结果将有助于进一步明晰生态系统服务功能和遥感植被指数的关系，为区域生态系统服务功能评估与生态保护提供科学参考。

1 数据与方法

1.1 研究区与数据来源

本研究所采用的青藏高原边界源自张镱锂等的划定结果，源自全球变化科学研究数据出版系统[17-18]。青藏高原的保护区边界源自世界保护区数据库(WDPA, world database on protected areas)(https:∥www.protectedplanet.net/)。考虑到数据库中保护区矢量图层相互交叠，以取最大面积为规则，对图层拓扑关系进行调整，共提取保护区矢量区块105个。2000—2010年生态系统服务功能制图数据源自Zhang等的评估结果，包括碳固定、土壤保持、水源涵养、栖息地保护以及总生态系统服务功能5项[19]。所有生态系统服务功能图层标准化至0～1，并统计各保护区的平均值。

MODIS EVI数据源自地理空间数据云中国月合成产品(http:∥www.gscloud.cn)。在月合成影像向年值转化中，最大值合成法较为常见[20-23]。同时，也有研究采用生长季的平均NDVI作为各个季节的植被变化[24-26]。本研究拟同时采用以上两种方法，进行相互对比。由于青藏高原东南与西北的生长季时段差异较大，本研究将EVI小于0.1的部分作为非植被区域予以去除，而将其余部分作为有植被覆盖的生长季EVI，从而得到基于EVI年最大值和年平均值。

1.2 变化趋势识别

为避免最小二乘法方法对异常点拟合而产生的参数误差，Theil-Sen(TS)中位数斜率与Mann-Kendall(M-K)单调性检验作为一种非参数趋势识别方法逐渐被引入植被变化分析中[8,27-28]。TS-MK算法对植被变化的具体形式可表征为[29-30]：

(2)

(3)

(4)

(5)

式中：slope>0表示植被增加; slope<0表示植被减少；ti，j——年份; EVIi，j——第i和j年的观测像元EVI值； var(S)表示方差；n——时间长度;Z的大小表示趋势的显著性，基于S指数判定，若|Z|>1.64则趋势是显著的。

1.3 生态系统服务功能评估

当前，生态系统服务功能的评估方式尚未统一，一些评估模型在运算中采用土地利用参数赋值的方式，掩盖了同种土地利用类型中不同植被覆盖条件的生态系统服务功能差异。Zhang等改进了一种针对生态系统服务供给的生态系统容量评估方法，其以CASA(carnegie-ames-stanford approach)模型评估逐月净初级生产力(net primary productivity, NPP)，并以年NPP为基础指标，将不同生态系统服务功能表示为NPP及其对应的环境因子加权(表1)。由于NPP与植被覆盖程度直接关联，这种生态系统评估方式可以更有效的体现出年际变化，更有助于观测生态系统服务功能对植被变化的响应情况，并服务于生态系统动态评估。由于荒漠、水体难以估算其NPP，本研究生态系统服务估算主要以森林、草地、耕地等有植被覆盖的生态系统为对象。

表1 不同生态系统服务功能的生态系统容量评估方法

2 研究结果

2.1 青藏高原EVI变化

对比2000—2010年青藏高原EVI变化如图1所示。在栅格尺度上，EVI年内平均值合成与EVI年内最大值合成的趋势提取结果具有差异性。具体而言，EVI年内平均值合成结果中，7.44%的像元显著增加，63.49%的像元呈现非显著增加趋势；而在EVI年内最大值合成结果中，显著增加像元仅占3.38%，非显著增加像元也降至54.14%。相比之下，EVI年内平均值合成结果中27.98%的像元呈现非显著减少，而在EVI年内最大值合成结果中该比例提升至40.83%。同时，两组提取结果中呈现显著减少的像元极少，仅有EVI年内平均值合成结果中的青海湖区域呈现集中的EVI显著下降。值得注意的是，EVI非显著增加的像元和非显著减少的像元往往相互交叠，并不呈集中连片分布。因此，并不能有充足的证据将大面积的植被增加或减少直接归因于某一局部区域的生态恢复或破坏过程。将EVI统计至保护区尺度上，EVI年内平均值合成与EVI年内最大值合成的趋势提取结果总体上具有一致性。其中，EVI年内平均值合成结果中有3个区块出现EVI显著减少，此外有非显著减少区块37个，非显著增加区块64个。而在EVI年内最大值合成结果中，不存在EVI显著变化的样本，非显著减少和增加的区块分别为35个和69个。其中，面积最大的四个区块，即羌塘、阿尔金山、三江源、可可西里，全部呈现出EVI的非显著增加。总体而言，EVI变化趋势在像元和保护区2个尺度上均反映出青藏高原植被整体好转、局部退化的特征。

2.2 保护区生态系统服务功能变化

图2为青藏高原2000—2010年生态系统服务功能保护区平均值。如图2所示，青藏高原生态系统服务功能的保护区平均值体现为碳固定最高，土壤保持次之，栖息地保护最低。这是由于青藏高原大部分保护区如羌塘、阿尔金山等地尽管面积较大，却大部分属于荒漠生态系统，不能提供高质量的栖息地。水源涵养相对低也是由于荒漠生态系统广布所致。就时间变化而言，水源涵养和栖息地保护服务略有增加，而碳固定和土壤保持服务功能的变化不明显。

虽然遥感植被指数是估算生态系统服务功能的重要参数，但后者更多的考虑了区域实际的气候状况，二者不能等同。图3所反映的青藏高原保护区生态系统服务功能变化表明，植被增加并不一定意味着生态系统服务功能增加。具体而言，在保护区尺度上，青藏高原生态系统服务功能减少的区块数量大于服务功能增加的区块数量，与植被整体好转、局部退化的总体现象相反。总生态系统服务功能呈现显著减少的区块共9个，非显著减少的区块共66个，而非显著增加的区块仅29个。面积较大的生态系统服务功能显著减少区包括色林错和雅鲁藏布江中游河谷黑颈鹤保护区等。就各类生态系统服务功能而言，碳固定服务功能与土壤保持服务功能的变化规律比较相似，显著减少区块均为7个，非显著减少区块分别为61和60个，而非显著增加区块分别为36和37个。水源涵养服务功能的显著减少区块高达12个，非显著减少区块78个，而非显著增加区块仅为14个；与之对应，栖息地保护服务功能的显著减少区块仅为4个，非显著减少区块72个，而非显著增加区块为28个。同时，羌塘、阿尔金山、三江源、可可西里4个面积最大的区块均呈现出栖息地保护服务功能的增加，该服务功能与EVI的对应关系相比其他服务功能更为紧密，但在东南部其他区块仍显现出与EVI相反的变化趋势特征。该结果证明，2000—2010年青藏高原保护区的各项生态系统服务功能并未随植被增加有明显的增加趋势。

图1 青藏高原2000-2010年EVI变化

图2 青藏高原2000-2010年生态系统服务功能保护区平均值

3 讨论

目前，遥感植被指数已经成为表述植被变化的常用手段，虽然不同的遥感产品、不同的植被指数、不同的年内合成方法势必影响最终结果的表征，但全球和中国近年来的植被总体增加已经在学界达成共识[21,8,28]。本研究也证实，尽管EVI与NDVI在取值上确有差异，不同的年内合成方法也导致了一定幅度的结果改变，但不影响青藏高原植被总体增加的结论。然而，与以往认知有所差异的是，本研究认为植被增加并不一定导致除水生态系统以外的陆地生态系统服务的增加。该现象一方面可能源于生态系统服务自身评估方法的不确定性，而另一方面并更有可能是由于生态系统服务功能评估需要考虑当前的气候状况，未经过模型处理的遥感参数信息并不应直接作为生态过程、功能的表征途径。

近年来，InVEST(integrated valuation of ecosystem services and tradeoffs)等生态系统服务功能评估模型在生态系统服务评估中被广泛采用。然而，此类模型的部分模块如碳固定、生境质量等对土地利用数据的依赖较强，并不能细致刻画植被的疏密特征，以及气候的年际变化。本研究采用以NPP为基础数据的生态系统服务供给的生态系统容量评估方法，虽然在权重系数的刻画上较为简单，但可以有效反映气候变化下几项关键生态系统服务功能的时空变化格局。尤其是NPP本身估算采用CASA模型，该模型已相对成熟并在近年研究中取得了大量应用。NPP所反映的碳固定生态系统服务功能呈现非显著退化，即反映了碳固定能力的变化并不能等同于植被变化，这也是符合植被年内季节生长的碳循环规律的。在未来研究中，可以考虑将生态模型的计算结果与遥感参数观测结果进行进一步对比，从而深入验证遥感植被指数的实际应用能力。