生态系统服务是指生物体从维持地球生命支持的生态系统中获得的利益^[1], 其可持续供给是经济社会可持续发展的基础^[2]。生态系统服务价值(Ecosystem Services Values，ESV)作为评价区域生态系统服务强弱的一种方法^[3], 开展ESV的科学评估并对其时空演化特征进行定量描述, 对区域制定合理的生态保护政策具有重要的意义^[4]。土地利用/覆被变化是人类活动与自然环境相互作用最直接的表现形式^[5], 是引起全球环境变化的重要因素^[6]。土地利用变化影响生态系统格局与过程, 改变着生态系统提供产品与服务的能力^[7], 从而影响着ESV的变化。近年来, 随着3S技术的发展及运用, 土地利用变化与生态系统服务价值的影响研究成为热点。

自Costanza等^[8]实现了全球ESV的量化后, 掀起了ESV估算热潮^[9]。国内, 谢高地等^[10-11]基于Costanza等^[8]研究建立了中国生态系统单位面积服务价值当量表, 并被国内学者广泛使用。由于生态系统存在异质性和复杂性, 直接引用当量因子法或者只考虑ESV的静态估算价值, 将难以完全解释区域ESV的动态变化^[12]。因此, 修正当量表与研究区当量表之间的差异成为了研究热点。目前, 有学者在对谢高地等制定的当量因子进行修正的基础上开展了区域ESV的研究。如姜栋栋等^[13]利用植被净初级生产力(Net Primary Productivity, NPP)进行区域修正；王永琪等^[14]利用谢高地等制作的生物量因子表进行修正；刘倩等^[15]采用区域差异系数、社会发展系数进行调整；雷军成等^[16]利用植被覆盖度指数进行修正；谢高地等^[17]、潘洪义等^[12]利用净初级生产力、降水量、土壤保持量数据计算其调整系数。但在修正过程中, 同时考虑自然和社会经济因素对ESV的影响仍然关注较少。随着ESV研究的不断深入, 针对不同的空间尺度(行政区域、流域等)、生态系统(森林、河流、草地等)以及功能都有了大量研究^[18], 但鲜有专门针对山地等典型地区的研究, 尤其是我国西南生态环境脆弱山区ESV的时空分布特征及演变趋势仍难以完全被表征。

岷江上游地处我国地形阶梯中第一级阶梯与第二级阶梯过渡地带^[19], 是重要的生态屏障区、生物多样性区和典型生态环境脆弱区。除此之外, 岷江上游也是集资源供给、生态服务、环境调节等多种功能与经济社会发展于一体的复杂系统, 对整个岷江流域乃至长江上游地区的生态平衡具有重要意义^[20]。鉴于上述情况, 利用空间异质系数、社会经济调整系数、资源稀缺系数对当量进行修正, 构建了更能反映研究区实际情况的ESV模型。且基于遥感解译的2000—2015年土地利用数据, 利用CA-Markov模型预测了2035年岷江上游土地利用格局, 进一步分析岷江上游土地利用与生态系统服务价值的时空演变趋势, 为区域可持续发展提供科学依据。

1 研究区及数据来源 1.1 研究区概况

岷江上游介于102°32′—104°15′E, 30°45′—33°09′N之间, 位于四川西北部, 地处我国青藏高原东麓、横断山脉北端与川西北高山峡谷东侧的结合部(图 1)^[20]。地形结构复杂, 以高山峡谷地貌为主, 河谷深切, 山高坡陡。区域内多地震、滑坡、泥石流等自然灾害, 生态环境极其脆弱。该区域是长江上游生态保护屏障的重要组成部分, 也是成都平原生态保护的生命线^[20]。岷江上游流经的行政区域大部分在阿坝藏族羌族自治州, 主要包括茂县、松潘、汶川、理县和黑水五县, 面积24741.05 km²。根据阿坝州2016年年鉴, 2015年岷江上游总人口约为39.01×10⁴人, 其中农业人口约28.77×10⁴人, 非农业人口约10.24×10⁴人。区域内地广人稀, 羌、藏、回等少数民族聚居, 是民族迁移的“走廊”和交汇地带^[21]。地区经济发展整体落后, 是中国传统农耕文化区向游牧文化区的交汇面^[21]。

1.2 数据来源与处理

岷江上游2000、2005、2010和2015年的四期遥感影像来源于地理空间数据云网站(http://www.gscloud.cn/)。选取了云量最少的5—9月份的Landsat5/7和Landsat 8遥感影像数据为数据源, 轨道号为130/37、130/38和130/39。参考中国科学院的土地利用和覆被变化分类方法, 结合岷江上游实际情况, 分为耕地、林地、草地、水域、建设用地和未利用地。利用ENVI软件进行预处理, 具体步骤包括：几何纠正、图像拼接和裁剪等。再通过目视解译, 得到四期岷江上游土地利用分类图, 经过随机抽样检查Kappa检验精度达到96%以上。

生态系统服务价值估算所用到的耕地面积、粮食作物产量、CPI和修正系数计算所需要的社会经济数据来源于全国和地方统计年鉴以及国民经济和社会发展公报(http://www.stats.gov.cn/tjsj/)。2015年粮食平均收购价来源于《中国粮食年鉴》, NPP计算所需要的全国和地方的降水和气温资料来源于中国国家气象信息中心和中国环境状况公报以及水资源公报。驱动因子数据来源于资源环境科学与数据中心(http://www.resdc.cn/)、地理空间数据云(http://www.gscloud.cn/)和OSM(https://www.openstreetmap.org)。

2 研究方法 2.1 土地利用变化特征分析

土地利用变化幅度指各土地利用类型的面积在数量上的变化, 有助于了解研究区土地利用面积的变化趋势^[22]。计算公式如下：

(1)

式中, J为研究期某种土地利用类型的年变化量(hm²)；U_a和U_b分别为研究期初和研究期末某种土地利用类型的面积(hm²)；T为年数。

土地利用转移图谱, 利用ArcGIS空间分析模块中地图代数叠加运算, 可以反映不同时期土地利用类型转换情况。计算公式如下：

(2)

式中, Y为土地利用图谱单元类型；P为前期的土地利用类型；Q为后期的土地利用类型。从而生成两个年份间的土地利用转移图谱。

2.2 生态系统服务价值评估

依据Costanza等^[8]的研究结果, 生态系统服务可分为9种^[10]。本文参考谢高地等2007年修订的“中国陆地生态系统单位面积生态系统服务价值当量表”^[10]。由于谢高地等人的研究是基于全国生态系统的平均水平, 未能有效反映生态系统的空间异质性、生态系统资源和服务功能的区域稀缺性与需求性^[23]、区域社会经济发展水平的不均衡性。因此, 本文参考已有研究^[24-25], 并结合研究区域的实际情况, 通过空间异质系数、资源稀缺系数和社会经济调整系数进行修正。

2.2.1 单位面积生态系统服务的经济价值确定

参考谢高地等^[26]的研究, 单位面积农田食物生产生态系统服务价值(1个标准当量因子)相当于全国平均粮食单产市场价值的1/7。计算公式如下:

(3)

式中, E为单位面积农田生态系统提供食物生产服务功能的经济价值(元/hm²)；M为研究区每年的粮食总价值(元)；H为研究区耕地面积(hm²)。

为了减少不同年份物价差异带来的影响, 以2015年粮食平均收购价为基准, 运用CPI指数对不同年份价值进行修正。

(4)

式中, VC为生态系统服务价值系数(元/hm²)；X_ij为i类生态系统第j类生态系统服务的当量因子。耕地、林地、草地、水域、未利用地的价值系数参考谢高地的当量因子表计算, 耕地对应农田、林地对应森林、草地对应草地、水域对应河流/湖泊、未利用地对应荒漠。建设用地参考文献取值^[27]。

2.2.2 生态系统服务价值系数修正

(1) 空间异质性修正

生态系统服务的强弱直接影响生态系统服务价值, 而生态系统服务与生物量有着密切联系^[26]。一般情况, 生物量越大, 生态系统服务也就越强^[28]。因此, 本文用NPP反映生物量并作为空间异质系数进行修正^[15]。生产力的计算要考虑到数据获取难易程度且又能说明变化, 本文选用Thornthwaite Memorial模型^[29]。计算公式如下：

(5)

(6)

(7)

(8)

式中, N表示空间异质系数；NPP表示植被净初级生产力(kg hm^-2 a^-1)；NPP′表示研究区植被净初级生产力；NPP″表示全国植被净初级生产力；V表示实际蒸散量(mm)；r为年降水量(mm); L为平均蒸发量(mm)；T为年平均气温(℃)。

(2) 社会经济因素修正

各地区的自然和社会经济条件存在差异, 导致社会经济发展水平和生态系统服务以及人们对生态系统服务价值的认知具有差异性^[30]。因此, 本文从支付意愿和支付能力两个方面对生态系统服务价值系数进行修正。计算公式如下：

(9)

式中, D为社会经济调整系数；R为ESV的支付意愿；W为居民对ESV的支付能力。

居民的支付意愿会随着社会经济发展水平发生变化。利用皮尔生长曲线模型可以表示社会发展水平^[31]。计算公式如下：

(10)

(11)

(12)

式中, Z表示社会发展水平；Z′_m、Z′_g分别表示研究区、全国社会发展水平；E′_n为恩格尔系数；E_a、E_b分别为城镇、农村恩格尔系数；P_a、P_b分别为城镇、农村人口在总人口的比重；h、H为研究区、全国城市化水平。

支付能力是个人或者国家经济能力的表现^[32]。利用研究区和全国人均国内生产总值的比值进行ESV支付能力修正。计算公式如下：

(13)

式中, GDP'表示研究区人均国内生产总值；GDP″表示全国人均国内生产总值。

(3) 资源稀缺性修正

区域的人口越多, 人口密度越大, 人们对于资源的需求也就会越大。因此, 利用人口密度反映系数修正^[33]。计算公式如下：

(14)

式中, F表示岷江上游资源稀缺系数；g为岷江上游的人口密度；G为全国人口密度。

综上, 2000—2015年的生态系统服务价值调整系数见表 1。

修正后的生态系统服务价值计算公式如下：

(15)

式中, ESV生态系统服务价值(元)；A为各类土地利用类型的面积(hm²)；VC年份生态系统服务价值系数(元/hm²)；N空间异质系数；D社会经济调整系数；F资源稀缺系数。

2.3 生态系统服务价值的敏感性分析

对敏感性进行分析的目的是确定时间变化后价值系数对生态系统服务价值的影响程度。计算公式如下：

(16)

式中, CS是敏感性指数；VC_nk和VC_mk分别为调整前后k地类的生态系统服务价值系数(元/hm²)；ESV_nk和ESV_mk 分别为调整前后k地类的生态系统服务价值(元)。

若CS < 1, 说明ESV对VC不敏感, 缺乏弹性, 即结合区域特征修正后的ESV是合理^[34]。

2.4 土地利用变化模拟

元胞自动机(Cellular Automata, CA)是一种基于不连续的时空动力学模型, 其特点是时间、空间和状态都是离散的, 具有强大的空间运算能力, 可以有效地模拟系统的空间变化^[35-36]。其模型可表示为：

(17)

式中, S为元胞有限、离散的状态集合；U为元胞的邻域；t和t+1表示不同的时刻；f为局部空间元胞状态的转化规则。

马尔科夫(Markov)模型, 是根据马尔科夫随机过程理论形成的方法, 研究和分析随机事件的变化规律, 并预测未来变化^[37]。其模型可表示为：

(18)

式中, S_t和S_t+1为该时刻土地利用系统的状态；Z_ij为状态转移矩阵。

CA-Markov模型综合了CA模型所具有的空间动态演化优势及Markov模型的长期预测优势, 可以有效地表现土地利用时间和空间的变化情况^[38]。

3 结果分析 3.1 土地利用变化 3.1.1 模拟验证

本文利用CA-Markov模型, 对耕地、建设用地进行限制, 选用地貌、气象、距离因素作为影响土地利用变化的驱动因素, 包括海拔、坡度、降水、气温、GDP、人口、河流距离、道路距离和中心城镇距离, 制作适宜性图集并对土地利用格局进行预测。为了验证CA-Markov的准确性, 以2010年为基年, 预测了岷江上游2015年土地利用格局, 并与实际解译的2015年土地利用状态进行比较, 结果见图 2。检验得出, 总体Kappa指数为0.9462, 随机指数为0.9578、位置指数为0.9543、分层区位指数为0.9543, 说明预测结果可信度较高。通过检验后, 以2015年为基年, 迭代次数为20, 选择5×5连续性滤波器, 模拟2035年岷江上游土地利用格局(图 3)。

3.1.2 土地利用类型面积变化

由表 2和图 3可知, 岷江上游土地利用类型以林地和草地为主, 建设用地和未利用地面积最少。2000—2015年林地面积持续减少, 年均下降438.738 hm², 2015—2035年林地面积年均增加8 593.691 hm²。建设用地和耕地面积持续增加, 尤其是耕地在2015—2035年年均增加了2456.901 hm²。草地、水域和未利用地面积呈现波动变化, 2000—2015年草地面积年均增加33.906 hm², 但在2015—2035年大幅度下降, 年均下降11066.008 hm²。水域和未利用地在2010—2015年均分别下降59.706 hm²和40.554 hm², 在2015—2035年下降81.729 hm²和0.810hm²。

3.1.3 土地利用转移图谱分析

由图 4可知, 林地和草地在与其它土地类型转换中较为剧烈。2000—2015年林地转出到草地、耕地面积较大。建设用地的增加主要来源于林地和耕地的转入, 耕地的增加主要来源于林地和草地的转入。因此, 这一阶段主要是林地面积在减少。未利用地的增加主要来源于林地和草地的转入, 说明区域的林地和草地在一定程度上出现退化的现象。从空间上看, 林地转耕地和耕地转建设用地主要分布在松潘县和汶川县东部, 林地转草地主要分布在松潘和黑水县, 林地转建设用地和草地转耕地主要分布在松潘和茂县。草地和林地转未利用地主要分布在松潘县。

2015—2035年, 建设用地的增加主要来源于草地和耕地的转入, 草地转建设用地主要分布在松潘县和汶川县, 耕地转建设用地主要分布在松潘和理县。耕地的增加主要来源于草地和水域的转入, 大致沿着岷江上游河谷低地分布, 主要分布在黑水、茂县、松潘县以及理县和汶川县东部。林地增加, 从转移来看得益于草地和水域的转入, 主要分布在松潘县境内。水域主要转入林地、耕地和建设用地。

总体来说, 无论是2000—2015年还是预测的2035年, 建设用地和耕地的持续增加与林地、草地、水域的转出关系最为密切, 呈现出正相关。岷江上游地形复杂, 低地和人口多集中在河谷地区, 这些区域成为了耕地、建设用地主要的侵占地, 也说明随着岷江上游经济的发展, 城镇化和工业化的快速推进, 人类活动占用生态空间, 导致植被破坏, 生态退化。因此, 随着人口的增加, 协调人类需求与生态保护会变得愈加重要。

3.2 生态系统服务价值的变化 3.2.1 生态系统服务价值时间变化

由表 3可知, 所有土地利用类型产生的ESV中, 林地和草地贡献最大, 对岷江上游ESV有着重要影响。2000—2015年岷江上游的ESV总体增加了254756.926×10⁴元, 一方面虽然这一阶段林地面积减少, 但转出主要是向价值较高的生态系统, 如草地。另一个方面社会经济调整系数和空间异质系数的变化对其造成影响。2010—2015年ESV下降了166743.189×10⁴元, 这与林地、水域减少, 建设用地、耕地增加带来的生态负效应以及单位面积农田食物生产的经济价值差异有关。2015—2035年, ESV增加了56494.557×10⁴元, 这与林地面积增加带来的生态正效应密切相关。

由表 4可知, 各单项ESV由高到低依次为：维持生物多样性>保持土壤>气体调节>水文调节>气候调节>原材料生产>废物处理>提供美学景观>食物生产。说明岷江上游生态系统在维持生物多样性、保持土壤、气体调节、水文调节和气候调节这五个方面发挥着重要作用, 在未来仍然是区域生态系统提供的重要功能以及保护的对象。总体上, 2000—2015年和2015—2035年各单项ESV都呈现增长的趋势。但是, 2015—2035年各单项ESV增长的幅度明显低于2000—2015年, 说明未来研究区土地生态系统提供各种服务的能力在降低。

3.2.2 生态系统服务价值空间变化

由图 5可知, 以村界为基础, 根据自然断点法从低到高对研究区2000—2035年ESV划分了5个等级。ESV较高的地区集中在岷江上游西部, 尤其是松潘县的西北部, 黑水县、理县和汶川县西部。这些区域海拔高, 林地和草地覆盖面积较大, 人口及经济活动少。岷江上游的ESV低值区主要沿着河谷地带延伸, 集中在中东部地区。由于岷江上游地形地质条件复杂, 河谷地带地势低平, 光热条件较好, 人口集中, 人为扰动较大, 土地利用变更也较西部更快。因此, 这些区域ESV较低。

由图 6可知, 对2000—2035年岷江上游ESV变化做冷热点分析。2000—2015年岷江上游ESV增加热点区域主要分布在西部地区, 损失冷点区域集中分布在中东部河谷地带。热点地区主要分布在人类活动较少的山区, 冷点地区主要分布在人类活动较为频繁的河谷地带, 说明人类活动对ESV的变化影响较大。2015—2035年ESV增加热点区域集中在西北部, 损失冷点区域仍然分布在东部, 相比2000—2015年冷点区域向理县东部和汶川东部集中, 黑水东部的冷点区域有一定规模的缩小。

综上, 2035年ESV的高低值分布大致与2000—2015年范围相同, 说明2035年的预测数据能在一定程度上反映岷江上游ESV的情况。并且岷江上游中东部地区和河谷地带是人口集中, 人类活动最为频繁的地区, 也是低值ESV和损失冷点的集中分布区, 尤其是黑水、汶川和茂县的东部地区。因此, 这些区域的生态修复是岷江上游地区生态建设的重点。

3.3 敏感性分析

本文对各类土地利用类型的生态系统服务价值系数分别上下调整50%, 根据敏感性指数公式, 分别求出2000年、2005年、2010年、2015年和2035年的敏感性指数。由表 5可知, 其中林地和草地在各时期的敏感性指数均较高, 说明两者的变化对ESV影响很大。主要原因是林地和草地面积在所有地类中占比最大。不同土地利用类型的敏感性指数, 在同一时期具有较明显的差异。但是每一种土地利用类型在不同时期, 敏感性指数变化都比较小, 且都小于1, 说明价值系数的调整对生态系统服务价值影响不大, 估算出的生态系统服务价值结果适用于岷江上游区域。

4 结论与讨论 4.1 讨论

岷江上游作为高原和盆地的过渡地带, 地质结构复杂, 自然资源丰富, 生态环境脆弱, 其地理环境在川西高原极具代表性。在气候变化和人类活动的双重作用下, 山地水土要素趋向失调、失衡, 对区域生态服务产生影响^[39]。因此, 加强过渡带研究不仅具有科学内涵^[40], 对于长江上游生态屏障建设具有重要战略意义。目前, 国内少数研究者对岷江上游的ESV进行了研究^[41-45]且与本文部分结果呈现出一致性。比如林地和草地对ESV的贡献率最大, 林地ESV变化最大, ESV的空间分布呈现西高东低的趋势等。本文在当量因子法的基础上进行了系数修正并且从时空尺度对土地利用格局和ESV动态趋势进行研究, 发现土地利用的变化会导致ESV的协同变化, 但修正系数的差异也会对ESV产生影响。比如, 2000—2005年和2000—2015年林地面积减少幅度较大, 但ESV增加幅度也较大, 造成这种现象的主要原因是修正系数的变化。因此, 合理确定ESV的修正系数, 对真实反应研究区ESV的变化具有重要的作用^[46]。

林地、草地是岷江上游高ESV的代表^[41-43], 且单项功能价值的变化与草地、林地面积变化的相关性大, 为了确保岷江上游生态系统发挥生态平衡的作用, 建议对这两种地类进行保护或修复, 引导形成合理且利于提高ESV的土地利用格局。以耕作和城市化等人类活动为主导的土地利用方式往往导致区域ESV呈现下降或者减少趋势^[47]。岷江上游河谷低地是耕地和建设用地增加的主要区域也是ESV低值分布的主要区域。因此, 建议通过制定合理的国土空间规划^[46], 加强生态保护红线的管控, 开展生态系统修复工作, 对ESV损失较为严重区域的生态系统进行修复并加以保护, 从而提高整个区域生态系统提供服务的能力。

就研究方法而言, 本文在使用当量法的过程中, 为了减少当量表与研究区现状的脱节, 结合区域实际情况进行了系数修正。虽然综合考虑了影响ESV的自然和社会经济因素, 但岷江上游是青藏高原东缘高山峡谷区, 自然地理环境复杂, 影响其ESV的因素有很多, 如, 地貌、自然灾害、政策因素等^[48]。因此, 还需要进一步的探索, 以便构建更为精准的ESV估算模型。基于预测的土利用数据模拟ESV变化的过程中, 2035年的修正系数沿用了2015年。虽然修正系数会对ESV产生影响, 比如2000—2015年ESV的变化。但由于经过敏感性分析, 生态系统服务价值对系数缺乏弹性。因此, 其影响在可接受的范围内, 仍然能对2035年ESV的变化做一定程度的反映。

4.2 结论

本文利用自然和社会经济修正因子, 构建了符合岷江上游实际情况的ESV估算模型。并且利用CA-Markov模型, 模拟了2035年岷江上游土地利用空间格局。从时空尺度分析了西南生态脆弱山区—岷江上游土地利用及ESV的动态变化。具体结论如下：

(1) 2000—2035年, 岷江上游土地利用类型以林地和草地为主。建设用地和耕地面积的持续增加与林地、草地和水域的转出呈现正相关, 草地、水域和未利用地面积呈现波动变化。从空间上看, 未来, 岷江上游人类活动的主要区域仍然在河谷地区。

(2) 岷江上游生态系统服务中, 林地和草地产生的ESV最高。2000—2015年岷江上游的ESV增加了254756.926×10⁴元, 2015—2035年ESV增加了56494.557×10⁴元, 其变化与地类的变化密切相关, 尤其是林地、草地和耕地的变化。各单项生态系统服务中, 维持生物多样性、保持土壤、气体调节、水文调节和气候调节这五项服务在岷江上游生态系统中起主要的作用, 是岷江上游生态系统保护应该关注的服务类型。

(3) 从空间上看, 2000—2035年ESV较高区域和变化热点区域集中在岷江上游西部, 岷江上游的ESV低值区和变化冷点区域主要沿着河谷地带延伸, 集中在东部地区。未来, 中东部及河谷地区的生态保护应该加以重视。

(4) 通过敏感性分析, 各土地利用类型的敏感性指数均小于1, 且每一类土地利用类型的敏感性指数变幅较小, 说明生态系统服务价值对价值系数不敏感、缺乏弹性, 研究结果可信, 并且林地和草地是生态系统服务价值变化的敏感因子。