在经济高速发展的影响下, 土地资源的交易和占用愈发活跃, 促使土地使用类型和利用方式快速转变^{[1, 2]}, 对生态系统服务价值(ESV, Ecosystem Services Value)产生了显著影响^[3]。自1997年Costanza在《Nature》提出ESV估算原理及方法后^[4], 我国学者开始寻找更适合于中国ESV评估的方法。比较有代表性研究有：欧阳志云等^[5]综合运用生态学及经济学方法, 探讨了大区域生态服务功能的内涵与评价方法。陈仲新等^[6]、毕晓丽等^[7]采用Costanza等人提出的生态系统类型单位价值对中国生态系统功能与效益进行了价值估算。何浩等^[8]结合遥感技术计算了2000年中国陆地ESV。谢高地等^[9]参考Costanza等人部分成果, 同时综合中国专业人士的生态问卷调查结果, 建立了中国陆地生态系统单位面积服务价值表, 在此基础上, 于2006年对中国500位具有生态学背景的专业人员进行问卷调查, 得出新的生态系统服务评估单价体系^[10], 并于2015年进行修订^[11], 构建了年内逐月和空间上分省的时空动态当量因子表。虽然后续也有不少学者进行ESV方法的改进与应用研究^[12-13], 其实质仍以谢高地等于2003年提出的基于单位面积价值当量表为基础进行。

由于计算ESV的价格基础不同, 出现基于单位面积ESV当量的计算结果远远高于基于单位面积价值量的计算结果^{[6, 11, 14]}, 直接弱化ESV在土地管理实践中的应用^[15]。张舟等^[15]通过分析Zang等^[16]两种ESV当量下的ESV变化情况, 发现ESV对当量参数敏感性小的特点, 并认为边际变化方法可用于ESV变化研究。因此, 现有文献常对两个时期ESV做差值, 进而讨论ESV的时序演变情况^{[8, 12, 17]}。事实上, 在排除各ESV当量系数影响敏感性小情况下, 单位面积价值成为引起ESV变化的主要因子。现有单位面积价值用研究区主要粮食单位面积产值的1/7进行替代, 该部分计算并未考虑物价上涨因素, 造成后续进行边际变化分析的ESV时序差值仍不是同一定基, 从而导致即使生态系统实际价值相同也会表现出以评估年份价格为表征的价值量不同现象^{[11, 14]}, 直接影响边际变化分析结果准确性。另外, 已有研究基于土地面积进行ESV变化研究, 更多关注于土地资源面积变化问题, 而忽略土地利用面积变化所引发的结构与格局的变化^{[8, 13]}, 特别是土地利用过程变化。即使一些学者研究农地非农化的代价性损失^[18], 仍未对土地利用时空变化过程与ESV变化的影响关系进行分析与探讨。

淮河流域地处我国南北气候过渡带和中原腹地。至2010年, 流域以不足全国13%的耕地面积和19%的播种面积, 承载了全国20.4%的乡村人口, 生产了全国23%的粮食, 提供了全国14%的农业产值^[19]。在此自然资源敏感性和经济需求迫切性并存的背景下, 本研究选取淮河上游区域, 进行土地利用及景观格局演变对ESV影响机理研究, 提出基于单位面积生态系统服务价值CPI(居民消费价格指数)系数修订方法, 并结合ESV损益方法, 进行ESV过程变化研究。旨在从宏观和微观空间尺度上, 定量评估土地利用变化引发的生态效应, 以期为淮河上游未来土地资源可持续利用和社会经济可持续发展提供理论和方法支持。

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

淮河上游处于我国南北气候的过渡带。在气候区划中, 以淮河和苏北灌溉总渠为界, 北部属暖温带半湿润地区, 南部属亚热带湿润区, 年平均气温为11—16℃, 并呈现出自南向北逐渐递减趋势。流域内降水时空分布特征明显, 具有南多北少, 山区高于平原, 年内降水多发生于汛期, 年际降水量变化大的特点。淮河上游南部由西向东依次为桐柏山和大别山, 山地丘陵占优势, 海拨1000m以上的中山面积约2000km², 海拔1000m以下的低山丘陵面积近30000km^2[20]。淮河上游北部主要为平原, 是我国重要的粮食、棉花、油料生产基地, 农业生产水平和土地利用率较高。

淮河上游流域包括信阳市、平顶山市、南阳市和驻马店市部分地区, 为便于社会经济资料的获取、分析及成果的应用, 本文把政区面积60%以上在淮河上游流域边界内的县(市、区)纳入研究范围。

1.2 数据来源

研究所用流域行政区划数据和6个年份的土地利用数据, 均来自于国家科技基础条件平台——国家地球系统科学数据共享平台——黄河下游科学数据中心, 30m DEM数据来源于中国科学院数据云平台。计算各生态系统单位面积服务价值所用耕地数据、粮食产值数据来源于1990年、1995年、2000年、2005年、2010年和2013年河南统计年鉴。消除物价上涨所用的CPI指数数据, 来源于国家统计局网站。土地利用数据源为landsatTM影像, 采用人机交互解译方法获得, 数据准确率达到90%, 满足研究需要。

1.3 数据处理方法

土地利用数据主要借助于ENVI 5.1和ArcGIS 10.2平台进行影像处理和空间分析计算, 景观指数采用Fragstats 4.2软件进行获取, 参与计算的栅格数据单元为30m×30m。相关系数计算基于SPSS统计软件。土地利用变化分析所用到的变化率、转移矩阵计算公式参照文献17, 景观格局指数计算方法详细说明参照文献22。下面着重介绍生态服务价值修订及计算方法。

(1) 生态服务价值

(1)

式中, ESV为研究区生态系统服务总价值(元);n为土地利用类型的数量;A_k为研究区内土地利用类型k的分布面积(km²);VC_k为第k类土地利用类型的ESV系数, 即单位面积上土地利用类型k的ESV(元 km^-2 a^-1)。

在式(1) 的基础上, 得出ESV损益计算公式：

(2)

式中, P_ij为第i类土地利用类型转化为第j类土地利用类型后的ESV损益;VC_i, VC_j分别为第i、j类土地利用类型的ESV系数;A_ij为第i类土地利用类型转化为第j类土地利用类型面积。

(2) 单位面积ESV修订

(3)

(4)

(5)

(6)

式中, CPI为居民消费指数;a_i为去除单位CPI, 为第一期环比CPI, 为现期同比定基CPI;当量数据参考文献表 1^[14];单位面积价值为研究区当年主要粮食作物单位面积产值的1/7。

2 土地利用结构和空间格局时空变化分析 2.1 土地利用结构的时间变化特征分析

土地利用变化主要从土地利用类型面积变化和变化速率两个方面进行。由表 1可以看出, 23年间研究区建筑用地面积比例持续上升, 耕地和草地面积呈下降趋势。至2013年底, 耕地占研究区67.2%的面积, 仍据主导地位, 草地所占比例由1990年的3.02%下降至2013年2.16%。以每5a为间隔计算各时间段的土地利用综合变化动态度, 分别得到0.48%(1990—1995年)、0.36%(1995—2000年)、0.09%(2000—2005年)、0.38%(2005—2010年)和0.11%(2010—2013年), 从中发现, 1995年以前是一个地类快速转化的时期, 1995年之后地类变化速度减慢, 2005年之后逐渐加速转化。

表 2地类转移矩阵说明了1990—2013年各地类的互转情况, 其中, B%为初期转出比重, C%为末期转入比重。23年间耕地转出面积最大, 为1801.11km², 按转出面积多少依次为：草地379.82km², 林地340.61km², 水域63.51km², 建筑用地31.46km²;转出比重最大的为草地, 转出率为34.52%, 其次为耕地7.12%。转入面积最大的为建筑用地, 达1034.93km², 转入率为16.66%, 其他地类依次为：林地791.10km², 耕地496.34km², 水域229.73km²和草地64.35km²。因为未利用地本身基数小, 微小的数量变化都会导致较大的比重差异, 23年间只有耕地有少量转为未利用地外, 未利用地并未发生其他变化, 因此, 忽略未利用地的转移分析。

2.2 土地利用变化的空间变化特征分析

借助于ENVI 5.1对两个时期的土地利用类型图进行土地利用变化监测制图, 找出变化区域, 并在ArcGIS中按5种坡度类别统计, 得出表 3。表 3表明随坡度增加, 各地类转移量均有明显下降, 并呈现地类自然属性效应。其中, 耕地转林地, 草地转林地两类转移更多发生在6—15°坡度区间。而林地转耕地、耕地转建筑用地则更多发生在0—2°区域。

根据表 1和表 3, 平缓区域的耕地不断转为建设用地, 而草地、林地又补给耕地, 造成耕地不断向山区扩张现象。该现象在2005年之前尤为明显, 2010年左右出现耕地转林地、耕地转水域现象, 林地、水域面积有所增加, 说明为了维护生态效应开始退耕还林、还草这样一个人类生态修复重建过程。

2.3 景观格局时空演变分析

景观指数高度浓缩景观格局信息, 能够反映景观结构组成和空间配置等方面特征^[21], 是定量描述景观格局演变及生态过程影响的重要方法^[22]。由于对斑块的分析常常是以地类为单位, 分析某类型景观的斑块特征, 因此, 本文选取类别层次的地类形状指标(FARC)和斑块聚合度指标(AI), 景观层次的蔓延度(CONTAG)、Shannon景观多样性(SHDI)、Shannon均匀度(SHEI)和散布与并列指标(IJI)进行研究区23年景观格局变化分析。

从图 2各地类FARC指数来看, 历年耕地分维数最高, 表明耕地斑块形状不规则程度最高, 原因是：耕地分布最广, 且转出量最多, 导致破碎程度高, 但其面积大(67%以上), 所以仍显示出高的聚合度。从时间演变来看, 耕地、草地、建筑用地、林地斑块形状随时间推移愈加不规则, 水域和未利用地形状则呈相反发展趋势;林地在2010年后聚合度下降、而水域和草地呈高的聚集状态。

从表 4景观尺度指标来看, 研究区CONTAG呈逐年下降趋势, 由1990年的68.10降为2013年的67.59。SHDI、SHEI和IJI均逐年上升, 且IJI上升幅度高于SHDI, SHDI高于SHEI, 这是水域、草地等地类形状指标趋于规则、聚集程度上升共同作用的结果。SHDI和SHEI值上升, 表明23年间地类破碎现象一直在延续, 优势地类不断减弱, 地类分布复杂化。CONTAG逐年下降, 即为优势度下降所导致的景观间连接性和通达性的下降。

3 ESV变化分析

生态系统的异质性、生态系统服务的复杂性, 造成ESV测量困难性^[15]。Gret-Regamey^[23]等认为ESV研究不是给区域生态环境加上一个准确的、绝对的价格标签, 而是探讨单位自然资产增加或减少时其服务价值的边际变化, 进而探讨一种合理的土地利用格局。已有研究表明, 边际变化分析更适合于研究ESV变化^[24-26]。

3.1 ESV时间变化特征分析

参照河南省1990—2013年统计年鉴, 进行主要粮食作物播种面积及产值的计算。在此基础上, 根据各年份CPI和式(3) —(6) , 得到剔除通货膨胀部分的单位面积价值, 结合各生态系统面积数据和式(1) , 计算研究区各年份ESV(表 5) 。

由表 5可知, 1995—2005年ESV波动下降, 2005年后呈明显上升趋势。11种生态系统服务功能亦呈现上述年际变化特征。其中调节服务功能增值最多, 其次为支持服务功能, 分别增加154.24亿元和29.79亿元。结合表 6, 6类ESV在1995—2005年间均有不同程度下降, 2005年之后明显上升。林地、水域和农田是ESV的主要构成部分, 23年间, 水域增值最大, 增值97.55亿元, 其次为林地、农田, 分别增加72.70亿元和39.12亿元, 湿地为唯一负增长项, 减少4.64亿元。

根据表 1和表 6, 水域和林地所占区域面积比例分别为3.26%、17.62%, ESV贡献量却占有33.67%和36.00%, 表明水域和林地是引起ESV提升的主要土地利用类型。耕地、草地面积的减少, 其ESV却在增加, 这一现象的产生在于研究区主要粮食单位产值的增加, 致使各生态系统单位面积价值增加, 抵消了面积减少带来的负面效应。

3.2 生态系统服务损益分析

基于式2和1990—2013年土地转移矩阵进行ESV损益计算, 得到表 7。由表 7可知23年间, ESV增加205.93亿元, 主要由农田转为林地、农田转为水域带来的生态系统服务增值构成;所有地类转水域均为正流向, 增值178.52亿元。ESV损失主要由林地、草地、水域转农田造成, 林地转农田损失最多, 达138.05亿元。所有地类转建筑用地均呈负流向, 损失80.34亿元, 该部分损失主要由农田转建设用地构成。整体来看, 1990—2013年ESV呈收益大于损失情况, 生态增值主要来源于林地和水域生态价值贡献。

3.3 ESV空间变化特征分析

从研究区地形分布来看, 沿西北-东南向, 由南至北, 依次呈山地、丘陵、平原三个梯段, 图 3各县市ESV亦呈现梯段变化规律：山区面积居多的泌阳县、桐柏县、信阳市、罗山县、新县和商城县6个时期的ESV普遍最高, 丘陵面积居多的确山县、光山县、潢川县和固始县6个时期的ESV较高, 平原区域的县市ESV较低。其中, 汝南县因其水资源丰富, ESV偏高于其他平原区域县市。各县市1995年和2005年ESV均呈明显低谷区, 2005年之后增值迅速, 2010年之后, 山区县市仍呈高速增长趋势, 其他区域县市2013年ESV基本与2010年相平。驻马店市是唯一的负增长区域, 究其原因, 城市建设迅速, 地类多样性下降严重。

结合表 7和图 3, 23年间生态系统服务总值增加, 增值主要来源于林地与水域生态价值贡献, 而林地与水域的分布多依赖于地形因素, 因此出现生态系统服务总值增加, 但又空间分布不平衡现象, 且耕地不断被转移为建筑用地所造成的负面效应已经呈现, 因此未来应更加重视自然资源保护和生态服务价值的空间均衡化, 使得人人得以享受生态环境改善所带来的增值服务。

4 ESV与土地类型变化、土地景观变化的响应关系分析 4.1 土地利用变化与ESV变化的相关性分析

采用研究区1990年、1995年、2000年、2005年、2010年和2013年6个年份的一级地类面积与相应ESV进行相关性分析, 相关系数为0.35, 在0.10水平上达到显著。表明经过CPI单位价值系数修正后的ESV与6大地类面积存在一定的相关性。那么6大地类的空间分布情况与其ESV存在怎样的相关性？以各县市为统计单元(图 3) , 以1990年和2013年的土地利用面积变化量与ESV变化量进行相关性分析, 其中, 变化量为前、后两个时期的差值, 两组变化量的相关系数为0.38, 在0.01水平上达到显著, 相关程度明显高于两个总量的相关性, 表明ESV在随面积变化影响的基础上, 主要受地类面积的空间分布影响。那么地类的时空演变过程是否与ESV存在一定联系？仍以县市为统计单元, 以1990年和2013年的土地利用变化率与ESV变化率进行相关性分析, 变化率参照文献^[17](式(3-3) )获得, 相关系数为0.28, 在0.05水平上达到显著。该值略低于ESV与空间分布变化量相关性, 表明地类变化的时间效应对ESV的影响程度要低于地类的空间位置变化。

但研究区地形差异大, 处于南、北气候交接带, 自古以来是我国重要的粮食、棉花、油料生产基地, 农业生产水平和土地利用率较高, 各地类转移频繁, 不同地类的景观指数与ESV又表现为怎样的相关性？基于此, 以下从土地景观格局角度, 分析不同地类景观指标与ESV的响应关系。

4.2 景观格局变化与ESV变化的相关性分析

利用Fragstats 4.2软件进行斑块个数(NP)、斑块平均面积(AREA_MN)、最大斑块所占景观面积的比例(LPI)、形状(SHAPE_MN)、形状维数(FRAC_MN)、聚集度(AI)6个景观指数的计算, 分别得到出6期农田、林地、草地、水域、湿地和荒地6个地类的景观指数。在此基础上, 对6个地类的每个指数与各期ESV进行相关性分析。得出下表：

表 8表明ESV与景观丰度和地类聚集度存在明显相关性, 与地类斑块大小相关性偏弱, 与各斑块形状复杂程度呈较强负相关性, 说明了景观破碎化对ESV的影响程度。

4.3 ESV时空变化的影响因素分析

已有研究表明, 土地利用变化影响生态系统服务, 但对影响的过程和机理认识仍然不足^[27]。相关分析是定量研究空间关系问题, 分析空间格局的重要方法和有效手段。运用相关系数量化土地利用指标与ESV的关系, 有助于凸现ESV动态变化规律, 进而探讨构建更为合理的生态系统空间格局。研究表明：(1) ESV与地类空间分布相关性显著, 表 3和图 3的空间变化一致性亦验证了该结论, 也与胡和兵等^[28]和欧阳志云等^[5]研究结果相一致。均表现为受地形地貌影响的水热条件好、利于植被生长的区域, 其生态系统的单位面积价值量较高, 而往往这些区域景观多样性和生态完整性也较好, 因此ESV与NP和AI亦呈明显相关性。(2) ESV与地类年际变化相关性偏高。从表 1可知, 1995—2005年地类转移最为频繁, 大面积耕地转为建筑用地, 而林地和草地转出补给耕地, 2010年后, 部分耕地转出为林地和草地。而ESV在相应年份分别出现波动下降, 和明显上升趋势(表 6) , 两者时间节点变化一致。虽然林地与水域面积比例小(表 1) , 受ESV当量影响, 两者对区域ESV总值贡献比重大, 成为ESV的关键“热点”区^[28-29]。也验证了, 对于不同类型的生态系统, 面积变化与其ESV相关性偏弱这一结论。(3) ESV与AREA_MN和LPI相关性偏高, 与FRAC_MN存在负相关, 与郝芳华等^[30]研究结果相一致。这是因为, 23年间景观破碎化在逐年加重, 景观优势度逐年降低, 这一现象在1995年左右表现最为明显(图 3) , 生态系统“热点”区域不断被破坏, 也说明了人为活动是影响土地景观格局的主要因素, 而城镇扩张引发生态用地的减少往往是导致ESV总值损失的主要原因^[31]。(4) 研究地类结构和格局时空演变特征, 细化地类结构、格局指标, 并按各指标与ESV的相关程度进行梳理, 分析自然和人为因素在ESV过程变化中的作用机理, 较好的克服了已往评价尺度过大, 评价内容过泛^[28], 结论针对性不强的问题。成为探究土地利用变化下生态系统服务响应特征研究的一种有益尝试。今后仍需加强土地利用、景观空间格局、CPI对最终的ESV影响定量判定研究, 以深入揭示土地利用变化对生态系统服务的影响过程。

5 结论

ESV与地类的的空间分布属性相关性最为显著, 与地类的变化速度相关偏弱;ESV与土地景观异质性和破碎性相关性显著, 对景观的形状特征相关性偏弱。

采用生态系统价值损益分析方法, 使价值变化更多关注在自然资源的过渡性损失方面, 弥补了现有生态系统服务在过程变化评价方面的欠缺, 更利于评价结果在土地利用实际管理中的运用。

采用CPI指数修订ESV系数, 能够消除时间序列内各年份物价上涨带来的价值变化, 避免了价格影响下的ESV研究结果不尽一致问题, 使得计算结果更贴近实际生态系统价值, 提高了多时间序列生态系统服务评价的准确性和可比性。

致谢: 感谢国家地球系统科学数据共享平台——黄河下游科学数据中心(http://henu.geodata.cn)提供数据支撑。