神东矿区生态系统服务功能评价

刘 英，魏嘉莉，毕银丽，岳 辉，何 雪

(1.西安科技大学 测绘科学与技术学院，陕西 西安 710054; 2.西安科技大学 地质与环境学院，陕西 西安 710054; 3.中国矿业大学(北京) 地球科学与测绘工程学院，北京 100083)

神东矿区是中国已探明储量最大的煤田，然而该矿区生态环境极其脆弱，抗扰动能力差，煤炭开采引起的一系列生态问题如地表沉陷、植被退化以及水土流失等，进一步影响矿区生态系统服务功能，威胁矿区的可持续发展。面对大规模开采和生态保护之间的矛盾，采取科学的方法评估神东矿区的生态系统服务功能(Ecosystem Service,ES)，掌握其时空格局的演变规律及驱动因素，探索煤炭开采对其造成的影响，对兼顾神东矿区煤炭开采和生态环境的协调发展至关重要，可为矿区走资源保护型开采与生态环境治理相协调的绿色矿业之路提供科学支撑。

近年来，随着人类生态保护意识的觉醒，生态系统服务日益被重视，一些生态系统服务评价模型应运而生，如ARIES(Artificial Intelligence for Ecosystem Services)、SoLVES(Social Values for Ecosystem Services)及InVEST(Integrate Valuation of Ecosystem Services and tradeoffs Tool)模型等，其中InVEST模型以其强大的空间表达功能和运行结果可视化功能，迅速成为生态系统服务功能定量评估和动态分析的有效手段，且该模型具有多个模块，对应多种生态系统服务功能。国外学者利用InVEST模型对夏威夷、委内瑞拉、哥伦比亚以及亚洲印度尼西亚等区域[1-3]的ES进行了评估。国内学者集中于InVEST模型的水源供给、土壤保持、碳固定、水质的净化及生物的多样性等模块，实现了多元化区域的生态系统服务功能评估，如横断山区[4]、韶关市[5]、白龙江流域[6]及黄河流域[7-8]等。目前国内对ES的研究日益成熟，大多将单项生态系统服务作为评价标准，缺乏对多种生态系统服务进行综合评价的研究。

因此，笔者基于InVEST模型计算神东矿区1990—2018年的水源供给、土壤保持、碳储存生态系统服务，对3者进行耦合建立综合生态系统服务功能(Comprehensive Ecosystem Service,CES)评价模型，首先基于矿区尺度对神东矿区CES的空间分布特征、变化特征及聚类特征进行分析，并利用地理探测器探究CES变化的驱动因素;然后基于矿井尺度分析神东矿区不同程度开采区、矿井采区和非采区及复垦区CES的差异，探求煤炭开采对矿区CES的影响。本研究旨在从不同尺度综合测度神东矿区CES的变化规律，为促进神东矿区绿色开采和生态可持续的协调发展提供科学依据。

1 研究区概况

神东矿区隶属于神东煤炭基地，地处内蒙古鄂尔多斯和陕西北部交错地带，覆盖黄土高原北缘与毛乌素沙漠过渡地带东部(图1)。矿区属温带干旱、半干旱大陆性季风气候，春季多风，夏季多暴雨，冬季长且寒冷;年降雨量为400 mm左右，年季差异较大，年均蒸散发达1 319 mm;平均风速为2.5～3.0 m/s，地势西北高东南低。矿区以风沙和黄绵土为主，北部、东南部为黄土丘陵沟壑区，梁峁起伏，沟壑纵横，地表支离破碎。地表原生植被种类单调，以耐旱、耐寒的沙生植物、旱生植物为主，呈稀疏灌丛景观。矿区煤层赋存稳定，构造简单，适宜机械化开采，主要以井工矿为主，多采用走向长壁全部垮落法开采，是全国最大的煤炭生产基地。

2 研究方法及数据

2.1 研究方法

2.1.1综合生态系统服务功能评价模型的构建

(1)水源供给功能。InVEST模型中的产水量模块基于水量平衡原理，利用降水量减去实际蒸散量计算每个栅格的水源供给量。计算公式为

(1)

式中，Yxj为第j类土地利用类型栅格x的年产水量，mm;Txj为第j类土地利用类型栅格x的实际年均蒸散发量，mm;Px为栅格x的年均降水量，mm。

(2)土壤保持功能。InVEST模型中的土壤保持模块即泥沙输移比模块以通用土壤流失方程USLE为基础，产生基于栅格单元的土壤保持量。本研究利用InVEST模型中土壤保持模块计算矿区水蚀区的土壤保持量，计算公式为

Sx-Retention=Sx-Potential-Sx+Dx

(2)

Sx-Potential=RxKxLx

(3)

Sx=RxKxLxPxCx

(4)

(5)

式中，Sx-Retention,Dx为栅格x的土壤保持量、泥沙持留量，t；Sx-Potential为仅考虑地貌和气候条件下栅格x的潜在土壤侵蚀量，t；Sx和Sy为考虑了管理、工程措施后栅格x及其上坡栅格y的实际侵蚀量，t；Rx，Kx，Lx，Cx，Px分别为栅格x的降雨侵蚀力因子、土壤可蚀性因子、坡长坡度因子、植被覆盖和管理因子和土壤保持措施因子，这些因子的计算参考InVEST用户指导手册[9]；rx为栅格x的泥沙持留效率；ry为上坡栅格y的泥沙持留效率。

神东矿区遭受风蚀的面积约占全区的70%，因而利用土壤通用风蚀方程(Wind Erosion Equation，WEQ)监测神东矿区风力侵蚀区的土壤保持量。该模型通过土壤风蚀速率与各因子之间的经验关系发展而来，被广泛的应用于风蚀区的研究中。计算公式为

Ex=WxKxLxBx

(6)

Ex-Potential=WxKxLx

(7)

Ex-Retention=Ex-Potential-Ex

(8)

式中，E为栅格x的风蚀侵力模数;Ex为栅格x的实际风力侵蚀模数;Ex-Potential为栅格x的潜在风力侵蚀模数;Ex-Retention为栅格x的风力侵蚀土壤保持量;Bx为植被覆盖因子，即研究区的植被覆盖度;Wx为风蚀气候因子，计算参考文献[10-11]。

(3)碳储存功能。碳储量模块基于土地利用类型和其对应的地上生物量碳库(陆地表层上所有存活的植物材料的碳)、地下生物量碳库(植物活的根系统中的碳)、土壤碳库(土壤中的有机碳)和死亡有机碳库(凋落物、倒立或站立着的已死亡的树木的碳)，将这4种碳库的碳储存量相加来评价每个地图单元和整个景观的总碳储存，公式为

Ci=Ci,above+Ci,below+Ci,soil+Ci,dead

(9)

(10)

其中，Ci为土地利用类型i的土壤及生物量总碳密度,t/hm2;Ci,above,Ci,below,Ci,soil,Ci,dead分别为土地利用类型i的地上生物量碳密度、地下生物量碳密度、土壤有机质碳密度和枯落物有机质碳密度;Ctotal为总碳储量，t;Si为土地利用类型i的面积，hm2;n为土地利用类型数。考虑到四大基本碳库已满足研究目的及数据的难获取性，模型中的采伐木材产品或相关木材产品斑块的第五大碳库不予考虑[10]。查阅相关资料和文献[12-17]并采用生物量碳密度和土壤碳密度与降水量、气温的关系校正模型[18-20]得到研究区各土地利用类型的碳密度(表1)。

表1 不同土地利用类型的碳密度Table 1 Carbon density of different land use types

(4)构建综合生态系统服务功能评价模型。考虑到神东矿区的水源供给、土壤保持、碳储存生态系统服务功能具有相同的重要性，赋予其相同的权重以构建综合生态系统服务功能评估模型SCE。计算公式为

(11)

式中，SEw，SEs，SEc分别为标准化后研究区水源供给、土壤保持、碳储存生态系统服务量，取值均为[0，1];SCE为综合生态系统服务量，取值为[0，1]，值越大表示综合生态系统服务功能越高。

2.1.2梯度结构相似度

基于梯度的结构相似度模型(Gradient-based Structural Similarity,GSSIM)[21]从图像组成的3个要素来解释图像的结构信息，公式为

GSSIM(x,y)=[l(x,y)]a[c(s,y)]β[g(x,y)]γ

(12)

其中，l(x，y)，c(x，y)，g(x，y)分别为2期图像x，y的亮度、对比度、梯度结构度比较函数,l(x，y)，c(x，y)分别由x，y的均值、标准差计算得到，g(x，y)通过梯度幅值来计算，梯度幅值则采用Sobel算子计算的梯度得到;常数α，β和γ用于调整各分量的重要程度，取α=β=γ=1。GSSIM值越高，表示图像x，y越相似，2者的变化越小。将GSSIM值分为3个区间:[0,0.25]为突变区，表示2期图像x，y发生了显著变化;(0.25,0.65]为中变区，说明变化处于中等水平;大于0.65的区域定义为低变区，说明两期图像变化并不明显[22]。

2.1.3局部莫兰指数

探究生态系统服务的空间聚类特征可有效把控其强弱分布的走向。空间聚类的高-高聚类要求要素具有高值且被其他具有高值的要素所包围，低-低聚类则要求要素具有低值且被其他具有低值的要素所包围。局部莫兰指数(local Moran’s I)可将空间关联结构模式可视化以衡量局域空间自相关性，对ES的空间聚集性进行测度[23]，其中空间联系局域指标ILSA(Local Indicators of Spatial Assoction)的公式为

(13)

2.1.4地理探测器

地理探测器[24]中的因子探测器可探测自变量对因变量空间分异性的解释程度，用q值度量，q值越大表示自变量对因变量的解释力越强，反之则越弱;交互作用探测器可以识别不同影响因子的交互作用，评估2个自变量协同作用对因变量的解释力。本研究利用地理探测器探究神东矿区生态系统服务空间分异的驱动因子。基于自然和人为因素2个方面，选取降水、坡度、高程、植被覆盖度及土地利用类型5个影响因子。地理探测器要求输入变量为类别数据，参考王劲峰提出的数据离散化方法并根据相关文献的先验知识[6,24-25]，将坡度按照<5°，5°～10°，10°～15°，15°～20°，20°～25°，25°～30°，30°～35°，>35°分为8级，植被覆盖度按照<0.3，0.3～0.4，0.4～0.5，0.5～0.6，0.6～0.7，0.7～0.8，0.8～0.9，0.9～1.0分为8级，将高程、多年平均降雨量按等距分类法分为9类，土地利用类型依据类别进行分类。王劲峰指出提取栅格格网点密度越大计算精度越高，考虑到软件可容纳的最大行为32 767，因此基于ArcGIS随机提取神东矿区30 000个点，作为地理探测器的运行数据。

2.2 数据来源

InVEST模型与WEQ模型的运行以多源数据为依托，本研究的数据来源与说明见表2。纵观神东矿区的开采历程，经历了4个阶段：① 1985—1990年，处于开发起步阶段;② 1991—2000年，矿区确立了“高起点、高技术、高质量、高效率、高效益”的矿井建造方针，处于探索提高阶段;③ 2001—2008年，矿区对小煤矿进行整合，进入集成创新开发阶段;④ 2009—2018年，矿区进入快速发展跨越阶段，高产、高效、高清洁的开采方式成为主流。因此，选取1990,2000,2005,2010,2015,2018年探究神东矿区各开采周期的综合生态系统服务功能的演变规律。另外，本研究用到的遥感影像有1990年的Landsat 5和2018年的Landsat 8影像，对其进行辐射定标、大气校正等预处理。

表2 数据来源与说明Table 2 Data source and description

3 结果与分析

3.1 基于矿区尺度的生态系统服务功能评价

3.1.1神东矿区CES时空分布特征

图2为神东矿区各年份的SCE，神东矿区SCE由0.448 6(1990年)下降到0.382 5(2000年)后持续增加至0.471 6(2015年)，2018年又下降到0.453 2，整体以每年0.009的速度呈增加趋势。1990年为神东矿区未开采期，自然生态系统受人为干扰较少，生态系统服务功能为较高水平;开采时期2000年生态系统服务功能减弱，但又在2005—2015年不断增强，且2010年超出1990年水平，这与神东矿区实行煤炭开采与矿区生态环境治理工作并重的方针政策密切相关，1990年虽为未开采时期，但是其植被覆盖极其稀少，神东矿区开采时期实行植被修复，植被覆盖率由开发初的3%～11%提高到60%以上。岳辉等[34]的研究也表明神东矿区生态环境质量在2005—2015年总体趋于稳定并向生态环境质量好的方向发展。

图2 神东矿区1990—2018年SCEFig.2 SCE of Shendong mining area from 1990 to 2018

由图3可知，神东矿区CES在空间上具有一致性，均以西北部最弱，整体呈由西北向东南逐渐增强的格局。1990年矿区CES在西北部高低值混合分布，东南部CES则以高值为主;2000年的CES较1990年有所减弱，尤以西北部减弱明显，CES低值区连片性有所增强，东南部CES也呈现退化趋势;2005，2010，2015年矿区CES持续好转，2005,2010年西北部矿井密集处的CES轻微增强，2010年矿区西南部CES有所增强，2015年矿区西北部的CES变化并不明显，集中在东南部呈现增强态势;2018年CES全区较2015年有所减弱，与1990年相比，矿区西北部CES呈减弱趋势，南部CES得到增强。整体来看，神东矿区CES往好的方向发展。神东矿区CES呈现西北部最弱，且由西北向东南逐渐增强的格局，可能是由于神东矿区西北部自1990年后工矿用地不断扩张，加之一部分未利用地稳定分布于此，两者的各项CES均较弱，导致西北部CES处于低水平，而东南部有部分草地和耕地覆盖，可在一定程度上抑制土壤流失、增强碳汇;另外，矿区降雨自西北向东南逐渐增多，水源供给受降雨分布的影响较为明显，且东南区域土壤侵蚀以水蚀为主，高降雨量区域的降雨侵蚀力较高会导致土壤潜在侵蚀大，土壤保持量因此表现为高值，因此水源供给和土壤保持均表现为西北弱东南强的特征。

图3 神东矿区1990—2018年生态系统服务空间分布Fig.3 Spatial distribution of ecosystem services in Shendong mining area from 1990 to 2018

3.1.2神东矿区CES空间变化特征

为定量分析神东矿区CES的空间变化特征，基于MATLAB平台利用GSSIM模型计算神东矿区1990—2000，2000—2005，2005—2010，2010—2015，2015—2018，1990—2018年6个时期的GSSIM并生成空间分布图(图4)。统计各变化区所占比例并结合图4可知，神东矿区CES空间变化在1990—2000年以低变区为主(占比为45.69%)，突变区占比16.12%，在矿区西北部分布连片，该区域在1990—2000年之间建成投产了寸草塔、大柳塔等矿井，这可能是导致此区域CES发生突变的原因;神东矿区CES空间变化在2000—2005，2005—2010年以低变区为主(占比为64.14%，51.37%)，在全区均有分布，中变区与突变区呈细碎状分布在矿区西部以及矿井密集的西北部，这是由于2000—2010年，矿区大部分的植被有所改善，植被退化的区域主要集中在矿区的西部和西北部地区;2010—2015年矿区东南部突变区明显增多，矿区CES发生明显变化的区域由西北部转移到东南部，这是由于2010年之后，矿区植被覆盖逐渐呈自西向东逐渐增加的规律，加之2013年矿区内降雨量骤然增加(降雨量达661 mm，高出平均51.2%)，东南部的植被状况长势越发良好，水源供给与碳储量均有所增强，与2010年相比形成突变;2015—2018年全区以低变区为主，突变有所减少仍分布在矿区东南部;1990—2018年仍以低变区为主(40.04%)，全区均有分布，突变区次之(30.72%)，呈细条斑块状分布，西北部较其余区域相对密集。

图4 1990—2018年各时段GSSIM空间分布Fig.4 Spatial distribution of GSSIM in various periods from 1990 to 2018

为进一步分析神东矿区突变区和中变区CES的变化状态，以1990—2018年GSSIM影像为例在图上任意选取突变区3个样点区(A，B，C)与中变区3个样点区(E，F，G);并利用1990，2018年Landsat影像和1990—2018年的生态系统服务影像差值图进行验证(图5)。由图5可知，1990—2018年，突变样点区A的植被覆盖得到了明显改善，结合1990和2018年的CES差值图可发现，此样点区的CES有所增强;突变样点区B由低植被覆盖区变成了工矿用地，差值图显示此样点区的CES减弱;突变样点区C由低值被覆盖区变为露天开采区的排土场，差值图显示其CES有所减弱;中变样点区D，E的地表覆盖发生了明显改变，均由低植被覆盖区变为露天开采区或建筑区，差值图显示中变样点区D，E的CES均呈减弱趋势;中变样点区F的植被覆盖有所增加，CES为增加状态。由此可知，神东矿区CES减弱的区域往往伴随着工矿用地扩张对植被覆盖区的侵占，而神东矿区CES增强的区域植被覆盖度呈增加趋势，因此神东矿区CES的空间变化与植被覆盖和用地类型的变化密切相关。

图5 突变样点区(A，B，C)和中变样点区(D，E，F)Fig.5 Mutation sample point area (A,B,C) and intermediate variation sample point area (D,E,F)

为定量分析突变区、中变区及低变区CES的变化规律，将各时间段CES变化差值图与对应时段突变区、中变区及低变区的矢量叠加，统计不同变化区生态系统服务量增加和减少区域所占面积百分比(表3)。由表3可知，1990—2000年CES的下降主要来源于低变区(占比为32.71%)，中变区次之(28.64%);2000—2005年、2005—2010年，尽管突变区和中变区增加区域面积在持续增加，但是CES的增加均以低变区为主(占比分别为53.14%，42.43%);2010—2015年，突变区、中变区和低变区CES增加区域比例依次为19.19%，21.15%，20.85%，此阶段突变和中变增加的区域对SCE增加的贡献增大;2015—2018年，低变区再次成为CES的下降主要贡献者。1990—2018年，突变区和中变区增加和减少的区域并无明显差异，低变区增加的区域大于减少的区域，且低变区占据全区的40.04%，说明1990—2018年CES的变化仍以低变区变化为主。总体上，神东矿区CES的变化主要是由于矿区内低变区的变化导致，中变区和突变区对神东矿区CES变化的贡献力次之。

表3 生态系统服务量增加和减少区域所占面积百分比Table 3 Increase and decrease percentage of CES in study area

3.1.3神东矿区CES空间聚类分析

本研究基于GeoDa空间分析软件，利用局部莫兰指数对神东矿区CES的空间聚类进行测度，以有效把控矿区CES强弱分布走向。由图6可知，神东矿区CES的空间聚类特征以不显著为主，低-低聚集次之，高-低聚集最少。不显著区域分布在矿区中部及东南部区域，低-低聚集区由1990年的西北部、南部及河流两岸破碎分布转变为2000—2018年的西北部连片分布，在1990年矿区开采初期，西北部CES持续处于最弱水平;高-高聚集区的面积较小，零落分布在河流两岸和矿区南部的部分区域，由土地利用类型图可知河流沿岸分布着CES较强的耕地，所以神东矿区河流沿岸的CES水平较高;1990年低-高聚集区的分布处于矿区西北部，2000年以后主要分布在矿区东南部，该区域1990年被草地和耕地覆盖，2000年以后被局部的工矿用地和建设用地破碎分割，因此呈低值区被高值区包围的格局。整体上看，神东矿区的空间聚类特征明显，呈“西北部最弱，河流沿岸最强，东南部不显著”的格局。此格局可能是气候、植被、人类活动等综合因素的影响所致，矿区降雨自西北向东南逐渐增多，西北部因采矿造成的植被覆盖破坏均为西北部低-低聚类区提供了动力。因此，应将矿区西北部作为实行“煤炭开采与矿区生态环境治理工作并重”政策的倾斜区，同时对河流沿岸的高-高聚类区实施保护政策，稳定此格局，并采用提高植被覆盖、及时复垦废弃矿区及改进绿色开采技术等策略，提高神东矿区的生态系统服务功能。

3.1.4神东矿区CES变化影响因素分析

由地理探测器结果可知，对神东矿区CES解释力最大的因子为土地利用(q=0.173)，其次为坡度(q=0.144)，然后依次为高程(q=0.106)，降雨量(q=0.102)、植被覆盖度(0.057)，且各影响因子均通过了置信水平99%的显著性检验。土地利用的解释力最高是由于不同地类的固碳能力存在差异，其蒸散能力与泥沙持留效果的不同导致水源供给和土壤保持功能也存在差异;坡度和高程对生态系统服务功能的影响主要体现在坡度小高程低的区域人为活动较为频繁，对生态系统的扰动较大;降雨量对水源供给的作用较为明显，土壤保持和碳储量受到间接影响，所以解释能力次于坡度与高程;植被覆盖度高的区域碳汇较强，可抵御降雨侵蚀及其他物理性侵蚀，产水能力高而蒸散系数低，对生态系统功能均有正向促进作用，然而尽管神东矿区的植被覆盖在逐渐改善，但沙地植被占绝对优势，弱化了对生态系统服务的促进作用。因此，神东矿区应着重进行以提高植被覆盖为导向的土地利用结构优化方式，人工调控植被演替方向与速度，优化草本为主、草灌结合的林分结构，采用微生物复垦等科学方法复垦矿区沉陷地，对未利用地加大绿化力度等。

交互探测器的统计结果(表4)表明，神东矿区生态系统服务量的影响因子交互作用均高于单因子。土地利用与其他因子交互作用居于首位，且与高程的交互作用最为明显(q=0.276);植被覆盖度与其余因子交互作用提升明显，与高程、坡度和土地利用的交互均为植被覆盖度单因子的3倍以上，与土地利用的交互作用提升最多;高程与坡度的交互作用为单因子的1.5～2倍;降雨量与其他因子的交互作用提升并不明显，与土地利用的交互作用相对较高。分析可知，神东矿区的生态治理应以优化土地利用格局为先导，根据矿区高程的分布特征实施草地修复、沉陷地复垦等政策;神东矿区在坡度30°以下的区域随坡度的递加CES不断增强，30°以上发生减弱，有研究指出[6]随着坡度的递增，土壤保持功能整体呈先增强后减弱的趋势，因此神东矿区应对坡度与植被覆盖、人类活动进行动态监测，掌握CES强度分布的大致走向，将工作重点放在植被覆盖度较低的陡坡(30°以上)区域，也可在缓坡区域(15°以下)实施生态管理措施以平衡采矿及其他人类活动的地表扰动，以期提高神东矿区生态系统服务能力。

表4 各项生态系统服务影响因子交互作用q值Table 4 q values of interaction of impact factors

3.2 矿井尺度CES差异分析

3.2.1不同开采强度下的CES差异分析

神东矿区煤层构造简单，主要以井工矿开采为主，是我国最大的井工煤矿开采基地，为了实现分区探究矿区的生态系统服务功能，通过不同区域的开采强度进行进一步分析。范立民等[35]指出煤矿的开采强度可由平面开采强度和立面开采强度决定，即开采强度主要取决于工作面长度和采厚大小2个因素。本研究尝试探究矿区内不同开采强度下CES的差异及变化规律。神东矿区开采强度的划分依据文献[35-36]的研究成果，将神东矿区划分为开采强度极高、强度高、开采强度中等和低强度4级，得到神东矿区开采强度分区图(图7)。极高强度开采区位于乌兰木伦河和窟野河之间，此处分布着大柳塔、活鸡兔等千万吨矿井;高强度采区的特点是工作面较为连续，采煤高度大(1.3～4.5 m)，单位面积采出煤炭量高;中低强度采区的煤炭作业方式并未全部实现机械化开采，少部分为整合煤矿，产能低，工作面较为破碎。

图7 神东矿区不同开采强度分布Fig.7 Distribution map of different mining intensity in Shendong mining area

将不同程度开采区的矢量与矿区生态系统服务功能图叠加，对不同程度开采区的SCE进行统计(图8)。从各程度开采区的SCE差异来看，6个年份SCE大小顺序:中强度开采区>低强度开采区>未开采区>高强度开采区>极高强度开采区，中强度开采区的综合生态系统服务最强，多年SCE平均值达到0.502 8;极高强度开采区的综合生态系统服务最弱，SCE平均值为0.430 8。各个程度开采区的SCE均在1990—2000年和2015—2018年呈下降趋势，1990—2000年未开采区SCE下降量为各区下降总量的22.34%，该值在2015—2018年减少为15.09%，而高强度采区和极高强度采区CES下降量在1990—2000年为各区下降总量的7.52%，18.19%，2015—2018年上升为21.67%，26.94%;2015—2018年与1990—2000年相比，导致矿区SCE减少的未开采区的贡献率变小，而高强度、极高强度的贡献率变大了，说明神东矿区在各个开采区实行的生态保护措施并未完全抵消掉采矿活动带来的负面影响。而且，在2010—2015年，极高强度的SCE发生下降，其他程度采区则处于上升趋势，说明极高强度开采区对SCE的增长具有抑制作用。张文凯等[36]的研究也表明，神东矿区中强度开采区的NPP最高，极强度开采区的NPP最低，与本研究结果一致。

图8 不同开采强度区SCEFig.8 SCE quantities in different mining intensity zones

3.2.2采区和非采区的CES差异分析

为进一步探究采矿活动对生态系统服务的影响，选取神东矿区的7个主要矿井对其采区和非采区的生态系统服务量及变化比进行统计(表5)。从采区与非采区的差异来看，补连塔矿井的CES均值由采区大于非采区，其变化比由2005年的34.48%大幅下降为2010年的0.92%之后保持低值，说明2010年之后补连塔矿井采区和非采区的差距大幅缩小;大柳塔矿井采区均值均小于非采区，变化比介于-12.32%～-5.30%，说明采区与非采区差异较为稳定;活鸡兔矿井采区均值大于非采区，2015年之后小于非采区，且2015年之后采区非采区差异增大;哈拉沟矿井采区均值小于非采区，变化比持续下降至2015年的-10.67%后上升，说明哈拉沟矿井采区、非采区的差异先变小后变大;乌兰木伦矿井采区均值小于非采区，采区非采区差异十分显著(变化比为-35.41%～-26.88%)，且该矿井的CES最弱，可能是由于该矿井位于干旱少雨、植被覆盖低的毛乌素沙漠边缘;榆家梁矿井由最初的采区大于非采区在2010年后逆转为采区小于非采区，变化比介于-3.85%～1.96%;石圪台矿井采区均值小于非采区，变化比由2010年的-14.03%上升为-5.49%，说明非采区和采区的差异在变小。

表5 各矿井采区与非采区生态系统服务量及变化比Table 5 Ecosystem service quantity and change ratio between mining area and non-mining area

采区与非采区生态系统服务量的绝对差异仅能表明2者的CES是否存在明显差异，差异的来源可能是采区与非采区生态环境的差异性或矿区生态环境治理的影响。因此为深入分析地下采矿活动的影响，需进一步对比两者的相对差异。神东矿区大规模建设开始于20世纪90年代，以1990年为采矿前的初始状态(此时采矿活动并未进行，采区、非采区未受采矿活动的影响)，以此时矿井内采区与非采区初始差异变化比作为参考基准，采后状态选取2018年，若采矿后的差异变化比小于参考基准，说明矿井采区受到了地下开采活动的负面影响，相反，则可能受到了生态保护措施的正面影响[37-38]。

将神东矿区的7个主要矿井采区与非采区范围与1990年、2018年的生态系统服务量叠加即可得到采区与非采区的初始差异和采后差异变化。由表6可知，① 补连塔矿井各个年份的采后变化比均高于初始变化比(差值为0.51%～7.74%)，说明采矿并未对采区CES产生显著的负面影响，可能受到了生态治理的正面影响;② 大柳塔矿井采后变化比在2010年之前大于采前变化比，2010年之后小于采前变化比，因此该矿井在2010年以后受到了采矿活动的负面影响;③ 活鸡兔矿井采后变化比在2010年之前大于采前变化比，2010年之后小于采前变化比，差值由1.54%(2005年)下降为-15.60%(2010)年，该矿井在2010年以后受到了采矿活动的负面影响;④ 乌兰木伦矿井的SCE最低，采前采区SCE大于非采区，采后采区SCE小于非采区，且初始与采后状态的差异变化最显著(采后变化比较采前下降23.10%～36.28%)，说明采区与非采区差异在采后大幅增加，采矿活动对采区CES的负面影响较为显著;⑤ 哈拉沟矿井采前采区SCE大于非采区，采后采区SCE小于非采区，且采后变化比小于采前变化比(差值为-10.98%～-8.90%)，说明哈拉沟矿井的CES受到采矿活动的负面影响;⑥ 榆家梁矿井采后变化比小于采前变化比(差值为-4.08%～-0.57%)，且该矿井未采取植被改善措施，说明地下采矿活动对采区CES产生了负面影响;⑦ 石圪台矿井采区CES低于非采区，且采后变化比小于采前变化比(差值为-6.51%～-4.58%)，说明石圪台采区受到了采矿活动的负面影响。总体上，补连塔矿井采区的CES受到了生态治理措施的正面影响，大柳塔、活鸡兔矿井采区的CES在2010年以后受到了采矿的负面影响，乌兰木伦矿井采区CES受到采矿的显著负面影响，石圪台、哈拉沟及榆家梁矿井采区的CES自采矿以来持续受到采矿的负面影响，或采矿的负面影响高于生态保护措施的正面影响。

表6 各矿井采区和非采区相对变化比Table 6 Relative change ratio of mining area and non-mining area in each mine

为缓解矿区干旱、贫瘠、地表沉降植被受损等生态问题，神东矿区于2008年开始和中国矿业大学(北京)合作开展微生物复垦关键技术研究与试验，主要是通过从当地土壤中筛选出适宜的丛枝菌根真菌，培养后接种于植物根系，扩大植物根系对土壤水分和养分吸收能力与吸收量，利用菌丝修复了断根，以提高植物成活率和植株生长量。

本研究为探究矿区微生物复垦工作对生态系统服务功能的影响，选取大柳塔矿井采煤沉陷地微生物复垦试验区，探究其CES对微生物复垦的响应。将大柳塔采区分为复垦试验区和其他区域，统计SCE的值和变化比(表7)。由表7可知，2000—2005年复垦试验区SCE与采区其他区域的变化比呈下降趋势(变化比分别为2.64%，2.03%)，2010，2015，2018年变化比呈持续增加态势(变化比分别为5.10%，6.36%，9.45%)，说明复垦实验区的SCE与采区其他区域的差距越来越大，生态系统服务功能往变强的方向发展。矿区微生物修复是一个生态持续演变过程[39]，在2008年之后矿区微生物复垦对CES持续起正向促进作用。毕银丽等[39]指出微生物复垦后，土壤质地发生改变、土壤保水能力增强，生态效应也得到了改善。

表7 大柳塔采区与复垦试验区的变化比Table 7 Change ratio between Daliuta mining area and reclamation experimental area

4 讨 论

鉴于已有生态系统服务的文献缺乏对生态系统服务的综合评价，本研究通过耦合单项生态系统服务综合测度神东矿区的生态系统服务量，并从矿区和矿井尺度考量了研究区生态系统服务的变化规律，弥补了对神东矿区生态系统服务多尺度研究的空缺。

大量学者对神东矿区及主要矿井的生态环境质量、植被覆盖度及土壤湿度进行了研究。笔者[40]指出神东矿区植被状况2000—2015年有改善趋势且未来会有较弱的退化趋势，本研究中神东矿区生态系统服务在2000—2015年逐渐增强，2018年发生减弱，与前人研究结果一致并相互验证。基于矿井尺度，笔者[37]指出2010年以前，采矿活动并未对大柳塔矿井采区的植被产生影响，这与本研究得到的大柳塔矿井2010年之后逐渐受到采矿负面影响的结论一致;岳辉等[34]指出乌兰木伦采区生态环境质量最差，本研究乌兰木伦CES也处于最低水平;雷少刚[38]指出榆家梁矿井在2005年前受到采矿活动的负面影响，本研究结果也表明地下采矿活动对榆家梁采区CES产生了负面影响。本研究更新了基于神东矿区矿井尺度的研究年限，进一步厘清了主要矿井是否受采矿活动的干扰，客观评价了矿井生态保护措施的实施成果。笔者等[22]指出矿区土壤湿度变化与地表覆盖高度相关且与采矿活动有关，这与本文得到的神东矿区CES与植被覆盖、土地利用类型变化密切相关保持一致。

本文构建了神东矿区综合生态系统服务评价模型，通过前人的研究进行间接验证，证明了该模型的综合评估结果能够定量反映神东矿区生态系统服务的演变规律。建议神东矿区将西北部作为实行“煤炭开采与矿区生态环境治理工作并重”政策的倾斜区，并着重进行以提高植被覆盖为导向的土地利用结构优化，利用微生物复垦等科学方法复垦矿区废弃土地等。研究成果可为综合评价区域生态系统服务抛砖引玉，而区域的多项生态系统服务的耦合机制仍然是一个具有挑战性的难题，应进一步强化决策过程与机理耦合，依托多源数据获取生态系统服务功能供给、需求、目标管理等参数信息，致力于多项生态系统服务多尺度优化集成模型的开发与应用。因而本研究下一步将选取多种生态系统服务功能，寻求多种耦合机制，以求得更加精确的实验结果，促使生态系统服务评估为区域生态系统服务可持续管理提供更精准可靠的决策支持。

CES受自然、社会和经济多种复杂因素的影响，探究生态系统服务的驱动因素时，本研究所选因子多基于自然因素，进一步研究可融入人口密度、经济指标、煤炭产量等，以探究自然社会和经济因素对神东矿区生态系统功能变化的综合影响。另外，基于生态系统服务功能与人类福祉的紧密联系性，也可探寻生态系统服务与社会经济系统的耦合关系，为区域提升生态系统服务功能、实现可持续发展提供生态决策。

5 结 论

(1)矿区尺度上，① 神东矿区SCE由1990年的0.448 6下降为2000年的0.382 5，持续增加至2015年的0.471 6又下降为2018年的0.453 2，整体呈增加趋势，速率为0.009 a-1;空间上神东矿区西北部CES最弱，整体呈由西北向东南部逐渐增强的格局;② 神东矿区CES空间变化特征明显，矿区CES变化主要由矿区内低变区的变化导致，中变区和突变区次之，神东矿区CES空间变化与植被覆盖、用地类型的变化密切相关;③ 神东矿区CES空间聚类表现为“西北部最弱、河流沿岸最强、东南部不显著”的格局;④ 神东矿区CES变化的影响因子解释力大小为:土地利用类型>坡度>高程>降雨>植被覆盖度，且因子两两交互的驱动力较单因子提升明显。

(2)矿井尺度上，① 神东矿区不同程度开采区CES大小依次为中强度开采区、低强度开采区、未开采区、高强度开采区、极高强度开采区;② 补连塔矿井采区的CES受生态治理措施的正面影响，使得采区CES大于非采区;大柳塔、活鸡兔矿井采区CES在2010年以后受到了采矿的负面影响;乌兰木伦、石圪台、哈拉沟及榆家梁矿井采区的CES自采矿以来持续受到采矿的负面影响，或采矿带来的负面影响高于生态保护措施带来的正面影响;③ 大柳塔复垦实验区CES往变强的方向发展。