矿产资源开采为国民经济的发展提供了保障, 但破坏了生态环境, 引发了土地损毁、植被破坏等一系列生态问题和生态危机^[1], 威胁区域的生态安全和可持续发展。矿产资源的开采会导致自然生态系统区域的减少或改变^[2], 进而影响生态系统服务的供给^[3]。生态系统服务与人类福祉密切相关, 矿区生态环境的恢复治理最终要服务于矿区的生态系统服务^[4]。厘清生态系统服务在生态环境变化下的响应机理和变化特征对指导矿区的生态重建工作, 协调区域经济发展和生态保护之间的矛盾具有重要意义。

作为矿区生产建设活动的空间载体, 土地对维系矿区生态系统功能起着决定性作用^[5], 分析土地利用变化过程, 并探索相应的生态系统服务变化规律成为研究矿区生态系统服务应对外部生态环境变化响应机制的重要途径。Townsend等^[6]基于土地覆被研究了美国东部阿巴拉契亚山脉露天采矿和复垦对生物多样性和美学景观等生态系统服务的影响。Xiao等^[7]将土地利用变化与生态系统服务价值评估方法相结合, 监测了兖州矿区采矿前后的生态系统服务变化。胡振琪等^[8]利用土地利用类型的时空演变特征对榆林矿区造成的土地动态变化进行了相关探索。因此, 当前对矿区生态系统服务的研究主要聚焦于土地利用变化下价值的测算及时空变化特征的分析, 缺少不同生态系统服务间耦合关系的探索。生态系统服务间关系密切且复杂, 分为此消彼长的权衡关系及同时增强或减弱的协同关系。明晰不同服务间的关系, 有利于区域的利益最大化^[9]。生态系统服务权衡与协同的研究最早出现于20世纪90年代, 其目的在于解决生物多样性保护的问题^[10]。进入21世纪后, 联合国发布的《千年生态系统服务评估报告》指出生态系统供给服务上升可能会导致其他类型服务的下降^[11], 此后权衡与协同的研究得到广泛关注, 其研究的领域也被进一步拓宽, 主要围绕理论框架构建^[12—13]、时空演变^[14—15]、驱动力^[16—17]和尺度效应^[18]等方面展开讨论。Delphin等^[15]以美国佛罗里达州为研究对象, 揭示了城市化是如何影响生态系统服务的时空演变以及权衡关系。Wuletawu等^[19]分析了埃塞俄比亚西南部地区农业生产与其他类型生态系统服务的权衡与协同关系, 进而为农业生产最佳投入、环境保护和管理工作提供依据。Li等^[20]的研究则表明, 在规划和政策的影响下, 关中盆地各项生态系统服务的关系容易在权衡与协同之间相互转换。由此可以看出, 近年来与土地利用选择决策相关的生态系统服务权衡与协同研究已成为土地利用管理的重点, 这种生态系统服务间的交互作用可作为今后国土空间布局优化的重要参考依据。

矿农城复合区是指在工业化、城镇化及城乡一体化的背景下, 配合资源城镇开发与转型的以地下矿产资源开发、加工与利用、地上农业生产及农业生物资源利用为主要产业发展形式的资源、社会、政治共同体^[21—22]。与自然生态系统相比, 矿农城复合区土地利用类型变化更加频繁^[23—25], 对生态系统服务的影响更为明显, 因此在生态系统演替过程中更加需要人的调控和参与。那么, 在矿产资源开采的背景下, 矿农城复合区的生态系统服务价值如何变化, 不同的生态系统服务到底呈现怎样的权衡与协同关系, 以及针对这种关系采取怎样的管理措施来保证区域生态系统服务的持续供给, 成为亟需回答与解决的问题。本文以平朔矿农城复合区为研究对象, 结合价值当量法和FLUS模型模拟研究区生态系统服务价值时空变化, 进而分析研究区目标年生态系统服务权衡与协同关系, 提出土地利用结构调整的建议和对策, 以期为矿农城复合区后续的生态系统服务管理工作提供参考。

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

平朔矿农城复合区位于山西省朔州市平鲁区, 地处东经111°52′—112°41′, 北纬39°21′—39°58′, 东西宽67.91 km, 南北长69.50 km, 总面积为517.48 km²。如图 1所示, 平朔矿农城复合区位于下面高乡、向阳堡乡、白堂乡、榆岭乡、陶村乡和井坪镇6个乡镇内, 包括安太堡、安家岭、和东露天三个大型露天矿和3个现代化井工矿, 煤炭储量约为127.15亿吨, 属于晋北煤炭基地。研究区所在地属于典型的温带半干旱大陆性季风气候, 冬春干旱少雨, 夏秋降水集中；年平均气温5.5℃, 年均降雨量410.66 mm。

作为典型矿农城复合区和人工生态系统, 研究区反映了大多数矿区亟待解决的生态问题, 因此研究该区域土地利用变化对生态系统服务的影响对于其他复合区生态系统服务权衡与协同的研究具有典型的示范性。

1.2 数据来源

本文以研究区1990年10月25日、1996年6月25日、2000年5月22日、2009年9月20日的Landsat 5遥感影像与2014年11月5日、2018年5月8日的Landsat 8遥感影像(空间分辨率为30m)作为主要数据源, 在系统辐射校正和几何校正以及影像裁剪的基础上, 参考《土地利用现状分类标准(GB/T 21010—2017)》和考虑复合区实际情况, 将遥感影像解译为耕地、林地、草地、水域、建设用地和损毁地共6类土地利用类型。其中, 损毁地包括露采场、排土场、剥离区以及与采矿活动相配套的工业场地。影像分类结果的Kappa系数在0.85以上, 基本满足本研究对数据的精度要求；DEM数据来自于地理空间数据云(http://www.gscloud.cn/search), 空间分辨率为30m。社会经济统计数据主要来源于《朔州市统计年鉴》(2018年)、《山西统计年鉴》(2018年)以及中国粮食年鉴(2017年)。研究区污染物排放情况和收费标准主要参考《朔州市主要污染物控制十三五规划》(2016—2020)和国家发展计划委员会、财政部等颁布的31号令(自2003年7月1日起施行)。

1.3 研究方法 1.3.1 土地利用类型情景模拟

(1) 模型简介

FLUS模型是由华东师范大学黎夏教授及其团队研发的一种土地利用模拟模型。与传统的元胞自动机(Cellular Automata, CA)模型相比, FLUS模型考虑到区域宏观需求、政治规划和气候变化, 将自上而下的系统动力学(System Dynamics, SD)需求预测模型与自下而上的CA本地分配模型进行整合, 进而提高了模型模拟未来土地利用格局的准确度。该模型主要包括土地利用需求预测和土地利用变化模拟两个模块。其中, SD模型将人类活动因素和自然生态因素相结合, 用于预测不同情景下多种土地利用需求；CA模型则被用于模拟未来的土地利用格局, 具体可分为两步：一是基于人工神经网络(Artificial Neural Network, ANN)计算各种土地利用类型出现在具体栅格上的可能性。二是利用自适应惯性竞争机制(Self—Adaptive Inertia and Competition Mechanism)解决不同土地利用类型间的竞争和交互。通过这两步, 具体栅格单元的所有土地利用类型的组合概率得以计算。在最终分配过程中, 具体的土地栅格是保留当前的土地利用类型还是转变成其他的土地类型取决于它们的组合类型和轮盘选择(Roulette Selection, RS)。本文对这部分的原理不再详细阐述, 具体参考Liu等^[26]的研究。

(2) 目标年选取与情景设定

在使用CA模块前, 需设定各土地利用类型变化的数量目标, 本文基于SD模型中的马尔科夫链(Markov Chain, MC)预测未来情景下土地利用类型的数量需求。由于马尔科夫过程是一类离散随机过程, 因此预测的土地利用需求时间点并不连续, 其时间间隔取决于需要输入的初始数据的年份。受到相关数据获取的限制, 本文输入初始数据为2014年和2018年土地利用类型数据, 所以可以预测2018年后以4年为间隔的时间点土地利用需求量。考虑到本文的研究目标是针对研究区近中期的生态系统服务变化情况提出合理的生态修复意见, 因此选定2026年作为目标年。

在研究区, 受到采矿活动和城镇扩张等因素的扰动，土地利用空间格局发生了变化，进而影响了生态系统服务间的权衡与协同关系。在不同的政策背景和保护目标下, 土地利用空间格局的演替方向也会产生分化。因此, 本文模拟了在自然、生态和经济发展情景下的未来3种土地利用空间格局。

① 情景一：自然发展情景。在此情景下, 土地利用格局的演替更多遵循历史时期土地利用本身的变化规律, 较少考虑到宏观政策的调控, 因此除水域外各土地利用类型间可以相互转换。

② 情景二：经济发展情景。经济发展情景以经济效益的最大化作为发展目标。在此情景下, 资本和劳动难以流向前期投入较大的土地复垦工程, 因此限制损毁地向其他用地类型的转换。另外, 考虑到农产品是当地农民的主要收入来源, 基本农田也将作为限制区域, 其用地类型在研究期间内不发生改变。

③ 情景三：生态发展情景。在此情景下, 采取较为严格的生态保护措施, 研究区的生态保护红线区域将被作为限制区, 其用地类型保持不变。另外, 在一定程度上限制林地向其他用地类型的转换。

成本矩阵代表各类用地相互转换的规则。如果某类用地可以转换为另一类用地, 则赋值为1；如果两类用地之间无法转换, 则赋值为0。基于3种不同的情景, 相应设计3种转换成本矩阵, 具体参数设置见表 1。

(3) 邻域权重设置

邻域权重参数用于设定各类土地利用类型的邻域因子参数, 参数范围为0—1, 越接近1, 代表该土地类型的扩张能力越强。借鉴人类扰动较为剧烈的快速城市化地区的研究成果^[27], 并结合2014—2018年研究区土地利用类型转换情况, 得到邻域权重参数表(表 2)。

(4) 土地利用变化驱动因子

土地利用变化是自然、社会、经济、政治等诸多影响因子复合作用的结果^[28—30]。在研究区, 受到上述各因素的共同影响, 土地利用类型在短时间内变化较大, 可能会对生态环境产生较强的负面效应^[31]。参考已有的土地利用变化驱动因素研究^[32—33]并结合矿区实际情况, 选取高程、坡度、坡向3个自然因素和到矿场、到中心城镇、到主要道路、到水域距离4个社会经济因素作为土地利用变化的主要驱动因子(图 2)。在FLUS模型中, 这7个驱动因素用于估算ANN中各种土地利用类型出现的概率。

(5) 模拟预测

在使用FLUS模型进行未来发展情景模拟前需进行有效性检验。首先, 利用2014年土地利用数据模拟2018年土地利用状况, 将得到的模拟结果与2018年真实土地利用数据通过FLUS模型精度验证(Precision Validation, PV)模块进行精度验证。结果表明Kappa系数为0.77, 总体精确度为0.808, 基本上满足精度需求。其次, 利用2018年数据并结合FLUS模型MC模块预测未来土地利用需求量。最后, 依据不同发展情景的限制条件, 模拟2026年各发展情景的土地利用状况。基于修正后的价值当量表, 对研究区的ESV变化进行分析。

1.3.2 ESV评估

生态系统服务评估方法主要包括物质量评估和价值量评估^[34]。物质量评估主要是指通过定点观测得到的监测数据, 再结合相关生态学模型对生态系统服务进行评估^[35]。物质量评估的优势在于通过动态物质量反映生态系统的生态过程, 进而反映生态系统服务的可持续性^[36], 但该方法不便于计算ESV的空间变化^[37], 且难以把不同类型的生态系统服务放在在同一维度比较。因此, 本文采用价值量的方法对各项生态系统服务进行评估。目前, 谢高地等人^[38]将Robert Costanza提出的17项生态服务功能合并成9项生态服务功能, 并制定了适合我国实际情况的陆地生态系统服务价值当量表, 其研究成果被我国学者广泛使用。因此, 本文参考价值当量表对耕地、林地、草地和水域这四种地类的生态系统服务价值进行评估。另外, Robert Costanza等人认为城镇生态系统服务价值为“0”, 而谢高地等人的研究成果则未提及城镇生态系统服务价值^[37—40]。本文的研究对象是矿农城复合区, 其中涉及到大量的建设用地和损毁地。考虑到研究区用地类型的特殊性, 因此参考张建军^[41]的相关研究, 采用成本替代法对这两种地类的生态系统服务价值单独进行计算。

(1) 耕地、林地、草地和水域ESV评估方法

考虑社会经济条件的差异性, 本文首先依据朔州市2018年主要粮食作物(玉米和大豆)的产量、粮食价格及种植面积对单位面积农田粮食作物的经济价值进行修正。单位面积耕地食物生产生态服务价值(1个标准当量)相当于研究区域当年平均粮食单产市场价值的1/7^[42], 修订方法如下：

(1)

式中, E_a为研究区单位耕地生态系统提供食物生产服务功能的经济价值(元/hm²)；i代表作物种类；M_i为主要粮食作物的经济价值(元)；m_i为主要粮食作物的种植面积(hm²)。另外, 考虑生态系统条件的时空差异性, 利用全国各省域农田生态系统的生物量因子表^[43]对研究区的生态系统服务价值系数进行修正(山西省农田生态系统生物量因子取0.46), 得出研究区耕地、林地、草地和水域的单位面积生态系统服务价值系数(表 4)。

(2) 建设用地和损毁地ESV评估方法

本文主要采用成本替代法对建设用地和损毁地的气体调节、水源涵养和废物处理这三项生态系统服务的价值进行计算, 计算方式如表 3所示。其中, 气体调节主要考虑研究区排放的空气污染物, 水源涵养主要考虑用水和排水对生态系统造成的负面影响, 而损毁地废物处理价值的核算主要考虑煤矸石和粉煤灰这两种研究区排弃最多的固体污染物。另外, 建设用地和损毁地仍为完整的系统, 其内部和外部的能流、物流和信息流始终也处在不断的交换中, 对生物多样性产生一定影响。因此, 生态系统的生物多样性保护服务功能价值参照荒漠生态系统对应项的数值。其他服务功能在建设用地和损毁地中体现较弱^[40], 故参考Robert Costanza等人的研究成果中的城镇生态系统对应服务功能的数值, 本文将忽略不计, 取“0”值。最后得到的建设用地和损毁地的生态系统服务价值系数如表 4所示。

为了更加直观的分析土地利用变化对ESV的影响差异, 本文选取格网作为评价单元进行权衡与协同的研究。目前, 以矿区或矿业城市作为研究对象所构建的格网单元主要包括500 m×500 m^[44]和1 km×1 km^[45]两种。其中, Bian等^[44]通过构建500 m×500 m的格网分析了徐州市贾汪区ESV变化趋势, 其研究对象和空间尺度与本研究有一定的相似性。另外, 景观生态学理论认为, 景观样本选取面积为斑块平均面积的2—5倍时能综合反映样区内景观的格局信息^[46]。在此基础上, 杨庚等^[47]研究表明平朔矿农城复合区在500 m×500 m尺度下各地类的面积变异系数较小, 格网同质性较一致。基于上述分析, 本文最后选择500 m×500 m作为基本的格网单元。

1.3.3 土地复垦收益的评估方法

国务院2011年颁布的《土地复垦条例》将土地复垦定义为：土地复垦是指对生产建设活动和自然灾害损毁的土地, 采取整治措施, 使其达到可供利用状态的活动。在矿区, 被复垦的对象主要是指因采矿而造成的损毁地。为简化分析过程, 本文通过ESV增量来衡量土地复垦带来的生态收益。在本文中, 土地复垦划分为两个阶段：损毁地转为复垦地(耕地、林地、草地和建设用地等)被划分为第一阶段；损毁地转为复垦地后, 低ESV的复垦地向高ESV的复垦地的转型(如草地转为林地)可以被认为是损毁地转为复垦地的后续生态收益, 被划分为第二阶段。

考虑到平朔矿区土地复垦工作自20世纪80年代后期展开, 因此将1990年作为研究的初始年, 将这一年的土地利用类型图与1996年的土地利用类型图相交, 即可得到1990—1996年土地利用类型的转换情况。其中, 损毁地转为其他地类作为1990—1996土地复垦第一阶段带来的ESV增量, 而转换后的新地类则作为1990—1996年的新增复垦地用来研究后续ESV增量。将1996和2000年的土地利用类型图相交, 可计算得到1996—2000年第一阶段带来的ESV增量。另外, 利用1990—1996年的新增复垦地图层裁剪1996—2000年土地利用类型转换图层, 通过计算可得到1996—2000年第二阶段带来的ESV增量。以此类推, 依次得到2000—2004年, 2004—2009年, 2009—2014年以及2014—2018年的土地复垦第一阶段和二阶段带来的ESV增量。

1.3.4 不同情景下生态系统服务权衡与协同分析

以划分的500 m×500 m格网为基础, 对供给、调节、支持和文化4类生态系统服务进行分区统计, 并利用SPSS 23和GeoDa 1.18.0 10软件对4项生态系统服务进行成对的相关性分析。SPSS 23的皮尔逊(Pearson)指数和Geoda 1.18.0 10的莫兰(Moran′s I)指数用来衡量不同类型的生态系统服务在空间上表现出的整体相关性, 主要是通过数值反映权衡或协同的关联性大小。另外, 为了进一步反映生态系统服务权衡与协同在局部空间表现出的差异性, 本文通过Geoda1.18.0 10软件进行双变量局部自相关分析。相较于前两种指数形式的非空间化表达, 双变量局部自相关分析可以借助局部空间相关性指标(Local Indicators of Spatial Association, LISA)图来更加直观地呈现研究区域的局部关联性情况。其中高高集聚和低低集聚表示协同关系, 高低集聚和低高集聚表示权衡关系。

2 结果和分析 2.1 1990—2026年研究区ESV变化特征分析

结合图 3可以看出, 研究区的整体ESV在1990—2018年呈现出先降低后升高的趋势。其中, 最高值和最低值出现在1990年和2014年, 分别为3.73亿元和0.96亿元。在单项ESV变化中, 调节服务的变化最明显, 其原因是调节服务价值的变化主要取决于损毁地数量的增减。1990—2014年间, 随着三大露天矿区的相继开发, 损毁地数量的快速增加, 导致调节服务的ESV快速下降, 并在2014年前后达到最低点。2014—2018年间, 虽然东露天矿区的规模继续扩大, 但安太堡和安家岭矿区有较多的损毁地被复垦成农用地, 因此研究区调节服务的ESV首次呈现出上升趋势。从图 4可以看出, 1990—2018年间废物处理和水源涵养的ESV下降趋势明显, 分别从0.31亿元和0.36亿元下降至-1.13亿元和-0.23亿元, 这说明矿区规模的扩大导致矿区排弃物处理成本大幅度增加, 废物处理和水源涵养这两项生态系统服务不仅不能提供正面的生态效益, 反而对周边生态环境产生危害, 严重威胁区域生态安全。相较于调节服务, 供给、支持、文化生态系统服务的价值变化并不明显, 原因是提供这两项服务的地类主要为林地、耕地等农用地。这些农用地本身数量较多, 因此小范围内用地类型转换不会引起这三项服务ESV的大幅度波动。

2018—2026年研究区土地利用格局变化显著, 主要表现在农用地内部的结构调整以及损毁地与其他地类的相互转换上。结合表 4可计算出, 2026年自然、经济和生态发展情景下的整体ESV分别为1.87、1.69和2.01亿元。经济发展情景在3个情景中价值最低, 原因是损毁地和耕地被限制转换为其他高ESV用地类型, 因此不能通过土地复垦和农用地内部结构调整提高整体ESV。生态发展情景下, 生态保护红线的划定一定程度上限制了该情景的林地和草地转为损毁地和其他农用地, 因此2026年生态发展情景表现出较高的ESV。

结合图 5可以看出, 1990—2018年间, 研究区ESV变化最为明显的区域是中部和东北部的矿区以及西部的城区, 这表明采矿活动和城镇扩张是引起研究区ESV变化的主要原因。另外, 在自然和生态发展情景下, 2018—2026年研究区中部的部分损毁地被复垦成农用地, 这将成为改善区域生态系统服务状况的关键因素。在研究区, ESV较为稳定的高值区主要分布在西南部林地, 这与黄土丘陵区本身的地貌特征密切相关。西南部是整个研究区海拔最高, 坡度最陡的区域, 由于地表起伏较大, 无论是露天采矿、城镇发展或是粮食耕种都会尽量避免向这个方向延伸, 受到人类的扰动较小, 因此该区域的ESV一直保持较为稳定的状态。

2.2 1990—2018年复垦地ESV变化特征分析

结合表 5可以看出, 1990—2018年间土地复垦共带来ESV增量约为1.96亿元。其中, 第一阶段和第二阶段的增量占比分比为95.55%和4.45%, 说明当前研究区ESV增量来源主要是损毁地转为复垦地。通过对比不同年份的数据可以看出, 第一阶段带来的ESV增量持续增长, 由2002.21万元增加至5924.69万元, 说明研究区土地复垦规模不断扩大对改善区域生态质量状况起到关键作用。另外, 第二阶段带来的ESV增量虽相对较少, 但其增长的速度相对较快, 且在总量中的占比不断提高, 因此有较大的发展潜力。

为进一步分析土地复垦第二阶段的ESV增量来源, 本文统计了1990—2018年间第二阶段的各类用地转换面积占比和ESV增量占比。由图 6可知, 土地复垦第二阶段的地类转换主要是草地转林地(397.64 hm²)、耕地转林地(278.13 hm²)以及耕地转草地(222.63 hm²), 分别占第二阶段地类转换的42.41%、29.66%和23.74%。这说明损毁地被复垦后, 耕地、草地和林地这三种农用地间的相互转换十分频繁。另外, 这期间ESV增量的主要来源是草地转林地和耕地转林地, 分别带来430.03万元和404.61万元的ESV增量, 共占第二阶段ESV增量的95.86%。

2.3 生态系统服务权衡与协同关系特征

为探究研究区在不同情景下生态系统服务权衡与协同关系, 本文以500 m×500 m格网为基础评价单元, 获取3个情景下每个格网单元单一服务的单位栅格面积的ESV, 然后进行SPSS相关性、双变量全局自相关分析。由表 6可知, 在不同的发展情景下, 研究区的4项生态系统服务均表现为协同关系。按照Pearson指数大小进行排序, 显著程度依次为供给-支持>支持-文化>供给-文化>供给-调节>调节-支持>调节-文化；若按照Moran′s I指数进行排序, 则为供给-支持>支持-文化>供给-文化>供给-调节>调节-支持>调节-文化。

为深入了解研究区不同生态系统服务在内部空间上权衡与协同关系的分异特征, 本文以自然发展情景为例, 使用GeoDa 1.18.0 10软件对成对的生态系统服务进行双变量局部自相关分析, 生成的局部空间相关性指标图可反映权衡与协同的空间聚类情况。从图 7可以看出, 生态系统服务之间协同与权衡关系的空间聚集特征总体较为显著。供给、调节、支持和文化四种服务在空间上主要表现为协同关系, 主要分布在研究区中部、北部、东北部等地区。其中, 供给—调节、供给—支持、供给—文化、调节—支持、调节—文化和支持—文化表现为协同关系的格网数分别为719、787、1048、1078、645和1078；表现为权衡关系的格网数目分别为47、32、79、145、482和49。北部和西南部地区协同关系表现为高高集聚特征, 是因为该区域林地数量较多且人类扰动较少, 生态本底质量较高；中部、东南部和西部的协同关系表现为低低集聚, 主要由于矿区和城区向外扩展导致；权衡关系大多聚集在矿区和城区的周边区域以及东南部的耕地, 由于人类扰动破坏稳定的生态系统结构, 导致表现出权衡特征。

3 讨论与结论 3.1 讨论

通过上述分析可以看出, 相较于供给、支持和文化服务, 调节服务受到采矿活动的影响尤为明显。一般而言, 调节服务的价值是间接的, 难以直接用货币衡量^[41]。因此, 本文采用替代成本法对调节服务的ESV进行核算, 这也是调节服务比较常用的核算方法之一^[48]。矿区属于人工/半人工生态系统, 人类的活动对其演替的方向影响较大。在“可持续发展”的战略背景下, 如何解决经济效益和生态效益之间的矛盾成为矿农城复合区未来转型发展道路上不可避免的问题。在露天矿区, 生态修复的主要方式是土地复垦, 因此本文通过划分不同的复垦阶段进而挖掘研究区ESV变化的驱动因素。从结果来看, 当前研究区生态系统服务状况的改善仍主要依赖于土地复垦第一阶段, 即损毁地向其他土地类型的转换。尤其在2014—2018年间, 随着东露天矿规模不断扩大, 平朔矿区采矿的重心开始由安太堡矿和安家岭矿向东露天矿转移, 安太堡和安家岭的损毁地被大面积复垦, 研究区生态系统服务状况得到显著改善。相较于第一阶段, 第二阶段虽带来的ESV增量少, 但增长速度较快, 表现出较大的发展潜力。另外, 土地复垦第二阶段的收益是发生在第一阶段地类转换后的后续收益, 因此带来生态收益的时间点明显滞后于第一阶段, 这导致在对比同一年不同复垦阶段带来的ESV增量时, 难以直接变现新增复垦地第二阶段带来的生态收益。随着复垦地规模的不断扩大, 土地复垦第二阶段带来的收益将对改善区域生态质量状况起到重要作用。结合图 6可以看出, 第二阶段损毁地被复垦成耕地后转为林地和草地的现象比较明显, 说明当前研究区对复垦地用途的确定仍存在着一些问题, 未来仍需综合考虑区域发展需求、复垦地土壤质量状况以及地形地貌等因素, 复垦地应宜林则林、宜草则草、宜耕则耕, 合理确定土地用途。

生态系统服务以多种复杂方式相互作用, 此消彼长的权衡、相互增益的协同是生态系统服务之间的典型关系^[49—50]。生态系统服务权衡与协同的研究旨在厘清各项生态系统服务变化的耦合关系, 基于未来的发展情景改进当前的管理与决策, 促进首要目标和其他多种效益间的平衡^[51], 实现经济发展与生态保护的并行不悖, 这既是IUCN提出的《基于自然解决方案全球标准》原则之一, 同时也是适应性管理的核心内容。本文研究结果表明调节、支持与文化服务表现为协同关系, 这与前人的研究成果较为一致^[52], 但供给服务与其他三项服务表现为协同关系的结论却有别于部分研究^[53]。一方面, 本研究区含有大量的损毁地, 该地类为各项生态系统服务的低值区, 这种负向的协同关系是使研究结果不同于其他研究的重要原因。另一方面, 生态系统服务的权衡与协同与空间异质性、时空尺度和研究方法等因素息息相关^[54—55]。当前许多研究得出的结论很大程度上只是反映特定区域和时间节点的现象, 并不具有普适性。另外, 调节—文化服务在局部区域出现较强的权衡关系, 说明土地利用格局存在着竞争关系, 需充分考虑各方利益进一步化解生态风险, 保证区域生态安全。结合图 7可知, 安太堡、安家岭及东露天矿ESV低—低协同特征明显, 因此未来生态服务功能的提升仍需围绕露天矿区的土地复垦工作展开。研究区西部和西南部地区表现为ESV高—高协同的特征, 因此需要加强保护, 尽可能减少人类扰动对其产生负面影响, 避免单项生态服务功能的退化引起其他生态服务功能的衰退。

3.2 结论

(1) 土地利用变化模拟预测结果显示2018—2026年间研究区整体ESV呈上升趋势, 在自然发展情景、经济发展情景和生态发展情景下的ESV分别为1.87亿元、1.69亿元和2.01亿元。这说明生态保护红线, 基本农田的划定以及矿区土地复垦等因素会对研究区未来ESV变化产生重要影响。

(2) 分析研究区供给、调节、支持、文化4种单项ESV变化可知, 1990—2026年调节服务的变化最明显, 且通过土地复垦可以极大地提升研究区的生态系统服务状况。1990—2018年间, 土地复垦第一阶段和第二阶段总共带来ESV增量1.95亿元。相较于第一阶段, 第二阶段目前来说占比小, 但其增长速度较快, 有较大的发展潜力。

(3) 在不同发展情景下, 矿农城复合区各项生态系统服务的权衡与协通关系表现较为一致, 均存在明显的协同关系。空间上, 损毁地的“低低协同”是使各项生态系统服务表现为协同特征的重要原因。此外, 研究区东南部在调节—文化服务上表现出较为明显的“高低权衡”关系。研究结论可为矿农城复合区的生态系统服务管理提供科学依据。