生态系统服务是指人类从生态系统中获得的益处^[1]。土地作为陆地生态系统的载体, 其类型、格局和利用强度都会影响生态系统服务^[2]。首先, 土地类型直接决定了生态系统服务的主导类型, 如农田生态系统粮食供给服务和自然生态系统各种调节服务^[3]；其次, 土地利用的变化引起生态过程发生相应的改变, 进而导致生态系统服务发生量变或质变^[4]；最后, 不同的土地利用强度也会造成所提供生态系统服务的差异, 这种对应关系往往需要较高精度的土地利用数据来反映^[5]。景观格局是由不同的土地利用类型组成并在空间上呈现不同的分布与组合形式^[6-7]。从生态系统服务供给与维持的角度出发, 景观格局的变化会影响生态系统物质流、能量流和信息流, 引起生态过程的改变^[8-9], 最终引起生态系统服务的改变。由于土地利用类型变化往往不能完全反映其布局形态的演化规律, 景观格局指标所展现的空间信息值得被特别关注。

目前, 生态系统服务与景观格局集成是区域环境政策和生态系统管理等跨学科领域的热点问题^[10-11], 但此类研究较多地直接以土地利用类型为基础评估生态系统服务, 通过相关和回归分析探究生态系统服务与景观格局的关系。例如：Su等分析了杭嘉湖4个生态区的景观格局和生态系统服务价值变化, 发现除了水源涵养和土壤保持外其他生态系统服务都与景观指标显著相关^[12]；Ayinuer等发现新疆艾比湖流域的13个景观指标中有8个显示出与生态系统服务有显著的相关性^[13]；Zhang等分析广西壮族自治区喀斯特地区生态系统服务价值变化对景观格局的响应中发现, 总面积、最大斑块指数等指标与生态系统服务呈正相关, 分离度和斑块丰富度则与服务呈负相关^[14]；Zhang等认为斑块大小会造成景观指标的明显差异, 从而影响景观指标与生态系统服务间的关系^[15]。然而, 此类研究多采用土地利用类型赋值的价值当量法评估生态系统服务, 进而解析生态系统服务与景观格局的关系, 忽视了在评估生态系统服务中使用土地利用类型, 而其与景观格局指标高度相关, 导致相关或回归分析在输入数据上存在天然的因果关系。

考虑到区域的差异性和方法的局限性, 仍有必要继续探究重点区域的生态系统服务如何响应于景观格局的变化。其中, 通过生态遥感指标实现生态系统服务制图, 基于独立于土地利用数据的生态系统服务评估结果, 解释景观指标与生态系统服务之间的关联, 可以更有效刻画出景观格局与生态系统服务间的复杂关系。汶川县是我国重要自然保护区——卧龙自然保护区的所在地, 研究区域内保护区与非保护区并存, 自然干扰(地震、滑坡和泥石流等自然灾害)和人为干扰(汶川县作为阿坝州经济发展较快的区域)明显, 研究该区域内生态系统服务对于景观格局及其变化的响应有助于了解保护政策的实施效果, 自然和人为干扰对区域生态环境的影响, 为区域生态保护、景观规划和生态系统服务提升提供理论依据。

1 研究区概况与研究方法 1.1 研究区概况

汶川县位于四川省阿坝藏族羌族自治州的东南部(30°45′—31°43′N, 102°51′—103°44′E), 四川盆地的西北角, 面积约4084 km²。境内地势由西北向东南逐渐降低, 西部海拔较高, 其中著名的四姑娘山海拔约6250 m, 东南部最低海拔约764 m, 河流发育, 属于典型的高山峡谷地形。汶川县属于温带季风气候, 由于地势起伏较大, 研究区域气候的垂直地带性明显, 形成常绿阔叶林、针阔混交等不同的植被带。该区域植被覆盖度较高, 生态系统碳固定作用明显；由于地处降雨丰富的高山峡谷地形, 高植被覆盖可有效增强区域的水源涵养和土壤保持服务。县辖8个镇, 4个乡, 为了保护珍稀动植物种, 国家建立卧龙和草坡自然保护区, 占全县总面积的一半以上(如图 1)。

1.2 数据来源与预处理

本文所采用的土地利用数据(2000、2005、2010、2015年4期)来源于中国科学院资源环境科学数据中心(http://www.resdc.cn/), 分辨率为30m, 其土地利用类型分为林地、草地、耕地、水域、建设用地和未利用地六大类；DEM和EVI数据源于地理空间数据云(http://www.gscloud.cn/search/)提供的GDEMDEM和MODEV1M产品, 分辨率分别为30m和250m；ET和NPP数据是源自美国地质勘察局(USGS)网站提供的MODIS影像MOD16A3产品和MOD17A3H产品, 分辨率500 m, 由于ET与NPP产品在网站中分别缺少2000年和2015年的数据, 在分析中分别2001年作为ET的起始年, 2014年作为NPP变化的终止年；月降水数据来源于国家气象信息中心(http://data.cma.cn/), 采用全国站点在Anusplin 4.2中完成基于地形的薄板样条插值, 并裁切得到汶川县2000—2015年月均降水量；土壤数据源自FAO的HWSD 1.2全球土壤同化数据库(http://westdc.westgis.ac.cn/)。为了保证数据分辨率的统一性, 本文将所有数据的空间分辨率通过三次卷积内插重采样到30 m分辨率。

2 研究方法 2.1 景观格局变化

为了明晰2000—2015年汶川县土地利用变化过程, 本研究采用土地类型变化趋势指数(P_s)为描述指标。该指数的取值在[-1, 1]内, 若P_s＞0, 则表示此类土地类型面积在增加, 越接近1则表示其他地类主要转为该地类；P_s＜0, 则相反, 详见以下公式^[16]：

(1)

式中, P_s为某一类土地类型的变化趋势指数；ΔU_in为其他类型转变为该类型的总面积；ΔU_out为研究T期间内某一类土地类型转变为其他类型的总面积；T为研究时间。

景观格局指数是反映其结构构成和空间特征的简单量化指标, 分为斑块水平、类型水平、景观水平^[17]。本文选择类型和景观2个水平进行分析, 从景观聚集程度、破碎程度、形状复杂程度和多样性等方面描述汶川县的景观格局特征, 在类型水平选取了最大斑块指数(Largest Patch Index, LPI)、斑块形状指数(Landscape Shape Index, LSI)、斑块密度(Patch Density, PD)、平均斑块面积(Mean Patch Area, AREA_MN)、边缘密度(Edge Density, ED)、聚集度指数(Aggregation Index, AI)6个指标；在景观水平选取了最大斑块指数(Largest Patch Index, LPI)、斑块形状指数(Landscape Shape Index, LSI)、斑块密度(Patch Density, PD)、平均斑块面积(Mean Patch Area, AREA_MN)、边缘密度(Edge Density, ED)、聚集度指数(Aggregation Index, AI)、香浓多样性指数(Shannon′s Diversity Index, SHDI)和蔓延度指数(Contagion Index, CONTAG)形成共8个指标。考虑到部分土地类型的覆盖面积较大, 较小的统计尺度会形成大量均质景观单元, 因而将景观格局指数统计到常见的千米尺度进行分析, 包含像元一千余个。利用ArcGIS中Create Fishnet工具得到3 km边长格网, 输入各年份土地利用类型数据, 通过Fragstats 4.2软件计算得到每个格网内的各项景观格局指数。计算景观格局指数使用了林地、草地、耕地、水域、建设用地和未利用地6类土地类型, 但由于研究区水域变化不明显以及建设用地、耕地和未利用地的总面积仅占全域的3.9%, 因此本文在叙述类型水平景观指数时只展示了林地和草地的景观格局指数。

2.2 生态系统服务评估

考虑到生态系统调节服务是自然生态系统最为重要的生态系统服务类型, 对生态系统调节服务的空间分析一直是学界研究的热点议题^[18]。尤其是汶川县大面积处于自然保护区中, 生态系统供给服务和文化服务相对而言不突出, 因而生态系统调节服务对景观格局及其变化的响应成为本研究的关注重点。

基于生态遥感指标, 本研究评估了最为常见的三种生态系统调节服务：碳固定、水源涵养、土壤保持^[19]。碳固定采用MOD17A3H产品中的各年NPP取值予以指代, 单位为g C/m²；水源涵养主要基于地形和土壤指标对产水量进行修正^[20-21], 土壤保持主要依循通用水土流失方程(The Revised Universal Soil Loss Equation, RUSLE)进行计算^[22], 具体过程如下。

2.2.1 水源涵养

在水源涵养量计算中, 首先通过水量平衡原理得到汶川县的产水量, 然后根据土壤的厚度、渗透性、地形和流速系数等因素计算得到涵养量, 具体的计算公式如下^[20]：

(2)

式中, Y、P和ET分别为产水量、降水量和蒸发量。

(3)

式中, WR为多年平均水源涵养量(mm)；V为流速系数(采用模型参数表数据^[21])；D为地形指数, 无量纲；K_soil为土壤饱和导水率(cm/d)。

(4)

式中, Watershed pixel count(Wpc)为集水区栅格数, 无量纲；Soil depth(Sd)为土壤深度(mm)；Percent slope(Ps)为百分比坡度。

2.2.2 土壤保持

采用RUSLE方程计算土壤侵蚀量, 考虑到大尺度人类工程措施无法表征, 默认土壤保持措施因子为1, 方程的公式如下^[22]：

(5)

式中, A为单位时间单位面积的平均土壤流失量(t hm^-2 a^-1)；R为降雨侵蚀力因子(MJ mm hm^-2 h^-1 a^-1)；K 为土壤可蚀性因子(t h MJ^-1 mm^-1))；L为坡长因子；S为坡度因子；C 为作物覆盖与管理因子；P为土壤保持措施因子。

R通过采用Wischmeier提出的相关公式^[23]计算得出, 如下：

(6)

式中, R为降水可持续性因子, P为年降水量(mm), P_i为第i月的降水量(mm)。

K采用EPIC模型^[24]得到, 计算公式如下：

(7)

(8)

式中, K为土壤可侵蚀性；SAN、SIL、CLA为砂粒、粉粒、粘粒含量；SNI为非砂粒所占百分比；C为土壤有机含量(%)。

坡度因子S根据坡度θ^[25]可计算得到：

(9)

坡长L根据以下公式计算得到^[26-28]：

(10)

(11)

(12)

式中, λ、m、β分别为水平坡长、坡长指数、细沟侵蚀。

考虑到EVI比常用的NDVI不易饱和, 更适用于植被密集地区, 因而根据EVI换算植被覆盖度c, 从而得到植被因子C^[29]：

(13)

其中EVI为12月的平均值。

土壤保持量应为潜在土壤侵蚀量与现实侵蚀量的差值, 其中潜在土壤侵蚀量应排除地表覆被类型, 即C=1：

(14)

A_r为土壤保持量, A_m为潜在土壤侵蚀量, A为现实土壤侵蚀量。

2.3 Mann-Kendall (M-K)检验

利用M-K显著性检验分析生态系统服务在连续时间内的变化趋势, 该方法常用于连续时间序列的趋势分析^[30-31], 同时采用Theil-Sen (TS)过程辅助, 其计算过程如下^[32-33]：

(15)

式中, Slpoe > 0表示增加, Slpoe < 0表示减小；n是时间序列长度, x_i和x_j分别是时间i和j的生态系统服务量；σ是标准偏差, | Z | > 1.64则认为变化趋势显著。

3 结果与分析 3.1 景观格局

汶川县土地类型主要是在林地、草地和耕地之间转换其中在2000—2015年间, 林地增加了112.21 km²、草地减少了139.58 km²、耕地增加了20.34 km², 而水域、建设用地和未利用地的土地变化相对较小, 分别为0、3.23 km²和0。为了强调研究区景观格局的变化效应, 考虑到研究区的水域、建设用地和未利用地的占地面积以及变化面积较小, 因此在后续分析中主要关注耕地、林地和草地变化对生态系统服务的影响。具体而言, 汶川县2000—2015年内耕地和林地面积处于缓慢增长的趋势, 草地面积变化趋势则相反。其中, 不同时间序列变化趋势指数P_s(表 1)表明, 2000—2005年间耕地和林地的变化趋势指数P_s均小于0, 即其面积处于减少趋势, 而草地的变化趋势指数P_s大于0, 面积处于增长趋势；2005—2010年耕地、林地的变化趋势指数均大于0, 面积处于增长趋势, 尤其是耕地增长量较大；草地相反, P_s低至-0.42；2010—2015年林地、草地和耕地的变化趋势和动态度指数P_s波动较小(-0.01、0.02、0.11), 变化均不明显。

2000—2015年景观指数的变化趋势见图 2, 分析结果表明研究时段内汶川县整体景观更加破碎化, 斑块形状复杂化, 斑块异质性增加。在景观水平, 整体景观的AREA_MN从452.57减少到424.37、AI从97.91减少到97.79, 而LSI从23.82增加到25.13、LPI从50.84增加到52.76、CONTAG从71.35增加到73.67、SHDI从0.81增加到0.83。草地的景观变化也显示出与整体景观较为一致趋势, 其AREA_MN、AI和LPI均表现出下降趋势, PD、ED和LSI有所增加, 表明景观破碎度增加, 聚集程度减弱, 最大斑块面积在逐渐下降, 形状趋于复杂。但是, 林地景观的变化规律则与上述两者不同, 林地的AREA_MN从2000年的1116.61增加到2015年的1870.84, LPI从50.81增加到52.76、AI从98.21增加到98.33, 但PD从0.05下降到0.03, 说明林地破碎化减弱, 团聚度增强。

3.2 生态系统服务变化

由图 3中可知, 研究区2000—2015年绝大部分地区水源涵养和土壤保持量呈增长趋势, 而NPP呈下降趋势。其中, 2000—2015年内汶川县水源涵养量波动上升, 除卧龙自然保护区的西部外, 大部分地区水源涵养呈增长趋势；另外, 草坡和卧龙保护区的水源涵养量主要在300mm左右波动上升, 而非保护区的水源涵养在250mm处波动上升, 即保护区的水源涵养要明显高于非保护区, 表明保护区的生物多样性保护起到一定的积极作用(图 3)。汶川县土壤保持量在15年间波动较大, 出现多处低谷(2000、2004、2007、2015年)和峰值(2001、2006、2008、2013年)；大部分地区土壤保持处于增长趋势, 但在城镇比较聚集的地区土壤保持量明显减少(图 3)。碳固定服务在2000—2015年呈减少趋势, 研究发现三个区域中NPP变化趋势一致, 呈现同升同降的规律, 可推测碳固定服务的变化受更宏观的气候条件影响较大(图 3)。三个区域的生态系统服务比较分析表明, 卧龙和草坡保护区内的水源涵养和土壤保持有明显的提高, 但在卧龙保护区与县域的交界处, 水源涵养和碳固定服务呈明显的下降趋势。

3.3 生态系统服务对景观格局的响应

将2015年生态系统服务评估结果以分位数法等分为低、较低、中等、较高、高5个等级, 根据每个等级上景观指数的波动范围, 得到不同生态系统服务等级所对应的景观指数(图 4和图 5)。

3.3.1 景观水平

景观水平上生态系统服务主要与LSI、SHDI、AI、LPI等指数密切相关(图 4)。以景观格局指数的中位数为对象, 当LSI分别在小于1.95、小于1.88时, 所对应的土壤保持、水源涵养量处于较高水平, 当LSI范围小于1.87时, 碳固定服务处于较低水平；当SHDI小于0.55时, 水源涵养和土壤保持均处于较高水平, 而碳固定则反之；AI在98.36以上时, 土壤保持和水源涵养量相对较高, 碳固定值较低；LPI大于85时, 土壤保持和水源涵养量较高, 而碳固定则反之。但三种生态系统服务对CONTAG并不敏感, CONTAG值一直在60左右波动。由此可知, 当LSI和SHDI较低、AI和LPI较高时, 土壤保持和水源涵养量较高, 而碳固定相对较低。此外, 箱图中显示的奇异值被界定为大于上四分位数1.5倍四分位数差的值, 或者小于下四分位数1.5倍四分位数差的值, 结果中奇异值较多存在与LSI、AI、LPI、SHDI等景观指数, 体现了格网中的这些指数变化幅度较大。说明汶川县景观格局具有较强的异质性, 而这种异质性特征很可能导致了生态系统服务空间分布的异质性。

3.3.2 类型水平

以景观指数的中位数作为参考, 生态系统服务对类型水平景观格局的响应与景观水平不同。由图 5可知, 对于土壤保持, 林地和草地斑块形状较为简单时, 会使得景观的土壤保持能力提高, 主要表现为林地和草地的LSI对土壤保持服务的响应与景观水平较为一致, LSI分别小于1.90和1.79时, 林地和草地的土壤保持量处于较高水平(景观水平为1.88)。增大草地的最大斑块面积是提高景观土壤保持能力的有效途径, 主要体现在LPI分别在17.83和81.85时, 林地和草地土壤保持达到较高水平；说明草地LPI对土壤保持服务的响应与景观水平(82.40)较为一致, 而林地则不一致。草地斑块破碎度较低时, 土壤保持量较高, 而林地土壤保持对PD的响应并不明显, 表现为PD在0.2上下波动。

水源涵养对景观指数的响应与土壤保持大多一致, 其明显区别在于, 增大林地的最大斑块面积可有效提高景观水源涵养能力的有效因素。其主要表现在林地与草地LPI在81.13和12.14时, 水源涵养处于较高水平, 而林地LPI与景观水平(82.40)范围大致相当。对于碳固定服务, 林地和草地斑块形状较为复杂时的碳固定相对较高。主要表现在林地与草地LSI分别低于1.89和1.42时, 碳固定服务处于低水平, 与景观水平LSI(1.47)对碳固定的响应大致相当。增大草地的最大斑块面积是提高景观碳固定能力的有效方式, 主要表现在林地和草地LPI分别在18.19、90.04时, 碳固定处于低水平, 其对草地LPI的响应与景观水平(89.97)一致。

3.4 生态系统服务变化对景观格局变化的响应

2000—2015年生态系统服务与景观格局变化量的对应关系表明, 只有碳固定变化明显响应于景观格局变化, 景观变化会引起研究区大面积范围碳固定值的降低, 而水源涵养和土壤保持在研究区中并没有明显响应于景观格局的变化。不同地类的景观指数变化对碳固定服务变化的敏感性存在差异(图 6)。以碳固定服务显著变化像元所占格网面积百分比的中位数为对象, 在景观水平上, 当LSI、LPI和AI发生明显变化时, 会造成约50%以上的碳固定服务显著减少, 而在LPI、LSI和PD稳定区间内该面积百分比中位数仅为29.98%、33.25%和32.36%。

在类型水平上, 生态系统服务响应景观格局现象更具有针对性。例如, 林地最大斑块面积显著减少、斑块形状趋向简化、聚集度减弱时, 引起碳固定在大面积范围内显著减少。主要体现在LPI、PD在减少区间内, 对应的碳固定显著减少所占面积百分比分均在约55%以上, 高于其他区间。草地的最大斑块面积增加、形状趋于复杂、聚集程度增加时会造成碳固定在较大范围内显著减少, 表现在LPI、LSI、PD在增加区间内, 碳固定显著减少所占面积百分比中位数均在60%以上, 而其稳定区间的百分比中位数低于40%。

4 讨论

本研究结果表明, 并不是所有景观格局变化均会引起生态系统服务变化。以显著变化的面积占到格网的百分比的中位数为对象, 土壤保持服务的变化量对景观格局的变化呈现出震荡波动的状态(图 7)。例如, 景观水平LSI、LPI、PD的不同变化区间内, 土壤保持显著增加所占面积百分比均在15%上下波动。在类型水平也是如此, 如草地也是在LSI、LPI、PD的不同变化区间内, 面积百分比也在15%上下波动。造成这种现象的原因可能是：虽然景观格局的改变影响土壤流失方程中的C值, 但在30m分辨率植被覆盖度遥感监测中这种变化不易被识别, 尤其是对高植被覆盖的保护区而言林地和草地均呈现较高的植被覆盖度, 导致遥感评估得到的土壤保持服务变化往往不依赖于景观格局变化。

不同尺度也会影响生态系统服务变化对景观格局变化的响应程度。在全县范围内, 水源涵养变化与景观格局变化并没有明显的响应, 但是, 在局部区域景观剧烈变化时, 水源涵养变化随景观格局变化的响应依然有明显响应。如在卧龙保护区的东南部选择34个格网单元的景观指数如LPI、LSI、PD等景观指数变化相对明显：中位数分别为18.29、-0.61、-0.35, 其水源涵养量显著增加面积中位数达56.34%；而卧龙保护区的东北部25个格网单元的景观格局变化较小：LPI、LSI、PD中位数分别为-0.89、0.07、0.11, 对应的水源涵养显著增加面积中位数仅为1.41%。

就研究方法而言, 本文揭示了生态系统服务与景观格局关系的测度直接依赖于生态系统服务评估手段, 高精度的生态过程地面观测数据补充将有助于进一步揭示生态系统服务与景观格局之间的复杂联系。此外, 卧龙、草坡保护区与其他地区的交界处往往是生态系统服务易于变化的区域, 因此, 保护区与非保护区生态系统服务供应存在明显差异。草地对研究区生态系统服务的影响较大, 优化区域区内合理的林地、草地规模配比与空间布局将更有效的提升研究区的生态系统服务水平。

5 结论

针对重点区域的生态系统服务如何响应于景观格局变化的科学问题, 本文识别了汶川县土壤保持、水源涵养和碳固定对景观格局及其变化的响应, 得到以下三点发现：

2000—2015年汶川县林地、草地和耕地之间存在转变, 整体景观更加破碎化、斑块形状复杂化, 斑块异质性增加, 整体景观格局的变化受草地影响较大。研究区碳固定服务有所减少, 而大部分区域水源涵养和土壤保持呈增长趋势, 在卧龙和草坡保护区内的水源涵养和土壤保持有明显提高, 但在卧龙保护区与县域交界处的生态系统服务呈下降趋势。

不同生态系统服务对景观格局的响应存在差异, 其不仅表现在指标上, 也体现程度上。在景观水平, 土壤保持、水源涵养和碳固定服务主要对LSI、SHDI、LPI等指数较为敏感, 而对CONTAG响应并不明显；在类型水平, 生态系统服务对类型水平景观格局的响应与景观水平有所不同, 林地类型水平与景观水平差异较大, 而草地类型水平与景观水平大致相当。例如土壤保持高值区间对应林地和草地LPI中位数分别为17.83和81.85, 而景观水平LPI中位数82.40。

汶川县生态系统服务与景观格局变化量的对应关系发现, 只有碳固定变化明显响应于景观格局变化, 景观的显著变化会引起研究区大面积范围碳固定值的显著降低, 而水源涵养和土壤保持在研究区中并没有观测到这一响应方式。在景观水平, 景观指数变化直接影响碳固定服务变化；在类型水平上, 林地LPI和PD的减少和草地LPI、LSI和PD的增加, 会导致碳固定服务的显著减少。