矿区煤炭开采活动通过改变生态系统结构和功能, 影响生态系统服务价值, 严重威胁矿区生态环境和以煤为生城市社会经济的可持续发展。生态系统服务是人类能够从生态环境中得到的各种恩惠, 包括有形物质产品和无形服务两方面, 可分为供给服务、调节服务、支持服务和文化服务^[1-3]。生态系统服务的概念早在1970年联合国大会上发表的《人类对全球环境的影响》报告中被首次提出^[4]。之后, 西方学者^[5]进展缓慢的开始对全球生态系统服务价值进行概略的评估, 没有形成相应的评估理论和方法体系。直到1997年Constanza等^[6]发表了《The Value of the World′ s Ecosystem Services and Nature Capital》一文, 在国内外学者中引起了很大的反响, 生态系统服务价值研究开始成为生态领域研究的热点。国内学者20世纪80年代开始对生态系统服务的概念和计算方法进行关注, 至21世纪初谢高地等人^[7]在Constanza等人研究成果的基础上, 采用问卷调查法制定出了适合中国陆地生态系统服务价值评估的单位面积价值当量因子, 有利的推动了我国不同尺度(全国、流域、省域或市域、县域)、不同生态系统类型(森林、农田、湿地、水域、草地、城市等)^[8-19]下生态系统服务价值的研究。目前, 对中小尺度矿区生态服务价值研究相对较少。对矿区生态系统服务价值的研究主要集中在采矿前后、复垦前后生态系统服务价值的变化情况, 以价值大小来反映采矿活动或土地复垦对生态环境的影响^[20-24], 对采矿过程中生态系统服务价值的变化研究少见。

生态系统服务价值评估方法大致可分为基于单位服务功能价格的方法和基于单位面积价值当量因子的方法。基于单位面积价值当量因子的方法要求数据少, 简单易算, 比较适用于区域和大尺度各类生态系统服务价值的评估^[25-26]。矿区是以矿业作业区为核心的一个独特的人工、半人工的区域生态系统^[27], 为人类提供多种产品和服务。与其他区域生态系统相比, 矿区是被人类扰动较大的土地受损区, 各区域土地受损程度不同^[28-29], 各生态系统具有很强的空间异质性。基于单位面积价值当量因子法评估矿区生态系统服务价值, 考虑不到矿区各生态系统内部服务的空间异质性问题, 评估出价值的结论可能存在着一定的偏差。一些学者采用不同方法, 针对自己研究区特点对单位面积价值当量因子进行了适当的修正。谢高地等^[30-31]利用某区某月/某年某生态系统的NPP、平均单位面积降水量和土壤保持量与全国范围内年均NPP、年均单位面积降水量和单位面积平均土壤保持量的比值对不同类型生态系统单位面积上各类服务功能年均价值当量进行了修正。韩增林等^[32-33]引入植被覆盖度系数, 在栅格尺度上根据植被覆盖度和NDVI的对应关系对生态系统服务进行修订。那么, 对于人工扰动强烈的、土地受损程度严重、空间异质性较强的矿区, 如何采用单位面积价值量法评估出更符合矿区实际的生态系统服务价值？开采时序、开采强度等开采活动对矿区生态系统服务价值又有何影响？

综上, 煤矿区是生态环境受损的集中区域^[34], 有必要对单位面积价值当量因子进行修正, 评估出矿区受损生态系统的生态系统服务价值, 从生态系统服务的视角研究采煤活动对矿区生态环境的影响。鹤岗矿区为我国典型的东北多煤层老矿区, 煤层厚、开采历史长。2013年沉陷面积达7932.12hm², 严重制约鹤岗城市的发展。因而, 该研究对单位面积生态系统服务价值当量因子进行修正, 对鹤岗矿区不同开采时期的生态系统服务价值进行估算, 分析采煤过程中鹤岗矿区生态系统服务价值的基本特征, 以期为生态系统服务价值的计算和煤矿区生态系统服务的恢复提供参考。

1 研究区概况

鹤岗市是一座缘煤而建的城市。鹤岗矿区(47°08′—47°23′, 130°13′—130°22′)行政区隶属于黑龙江省东北部的鹤岗市管辖, 位于鹤岗市辖区南部, 兴山区、向阳区、工农区、南山区的东部一带, 兴安区中部，地处三江(黑龙江、松花江和小兴安岭)平原交汇的金三角地带(图 1)。地貌单元属砂砾石丘陵台地, 海拔约100—140m。鹤岗矿区(主要包括兴山、益新、振兴、鸟山、南山、新陆、富力、兴安和峻德煤矿)煤炭开采始于1917年, 至今已有百年开采历史。现矿区开采范围北起兴山矿, 南至峻德矿, 南北长27km, 东西宽4km, 面积108km²。核定生产能力为1351万t/a。鹤岗矿区是我国典型的多煤层老矿区, 煤层最多的益新煤矿层组数达36个, 截至2012年末已开采或局部开采28个煤层；煤层最少的振兴煤矿层组数共计13个, 已开采或局部开采5个煤层。煤矿煤层多属于中厚及厚煤层, 总厚在44—90m, 富集在中生界白垩系下统石头河子含煤组的中部含煤段。经开采沉陷预计, 到闭矿, 总沉陷面积11794.02hm², 最大下沉值达39.4m, 严重制约鹤岗城市的发展。

2 研究方法 2.1 数据来源

土地利用矢量数据通过对1993年9月、2000年9月和2013年9月Landsat TM遥感卫星影像(30m分辨率)进行非监督分类获得。结合Goggle earth对识别不清的地类进行校正。Kappa系数依次为80.74%、93.06%和84.95%, 满足精度要求。基于Landsat TM数据, 参照中科院资源环境科学数据库分类法将1993、2000和2013年土地利用类型分为耕地、建设用地、林地、草地、水域和未利用地六大类。各期土地利用数据利用ArcGIS软件进行统计分析。其他社会经济数据均从黑龙江省统计年鉴和鹤岗历年统计年鉴中获得。矿区植被净初级生产力(NPP)通过CASA(Carnegie-Ames-Stanford Approach)模型模拟估算获得。降雨数据由鹤岗气象站提供。鹤岗矿区栅格尺度上的土壤保持量通过InVEST(Integrated Valuation of Ecosystem Services and Trade-offs)模型估算获得。

2.2 生态系统服务价值评估

采用基于单位面积价值当量因子法^[35], 经农田为基准的地区修正法对各生态系统服务进行修正^[36], 将计算全国生态系统服务价值的单位面积价值当量因子修正为黑龙江省的单位面积价值当量因子。再经矿区NPP、降雨量和沉陷因子将黑龙江单位面积价值当量因子修正为矿区单位面积价值当量因子。经公式(5)和(6)计算出鹤岗矿区不同类型生态系统和各服务的价值量。具体做法为：经公式(1)、(2)、(3)和(4)分别计算矿区1993、2000和2013年生态系统服务当量的修订系数。之后, 依据中国陆地生态系统单位面积生态服务价值当量表, 以2013年黑龙江省粮食平均收购价格(2.41元/kg)和1993—2013年农田单位面积平均产量(3779.38kg/hm²)为标准, 按照在没有人力投入的自然生态系统提供的经济价值是现有单位面积农田提供食物生产服务经济价值的1/7^[7], 计算矿区不同生态系统的单项生态功能服务价值, 其中建设用地生态服务价值参考袁兴中^[37]等的研究结果予以赋值。总服务价值经公式(7)进行敏感性指数检验。

2.2.1 地区修正法

采用以农田为基准的地区修正法, 修正公式如下：

(1)

式中, E_i为第i类土地利用类型经地区修订后的生态系统服务价值当量；Q和Q₀分别为研究区和全国的农田单位面积粮食产量；E_0i为第i类土地利用类型全国平均的生态系统服务价值当量, 其中i=1, 2, …, 6分别代表林、草、耕、水域、建设和其他土地。

2.2.2 NPP因子修正

NPP因子修正食物生产、原材料、气体调节、气候调节、废物处理、生活多样性保护和娱乐文化服务功能^[38-39]。计算公式如(2)所示。

侯湖平等^[38]以矿区耕地为例, 将CASA模型模拟四期矿区耕地的NPP值与其对应检验值(产量转为碳储量的算法作为检验值^[40])进行对比、验证。结果发现, 四期模拟值的平均值与其检验值的平均值残差为8.04gC (m²月)^-1, 并进一步对四期NPP模拟值与其检验值进行了相关性分析, 相关系数为0.75, 并对矿区NPP进行了模拟, 证明采用CASA模型模拟矿区NPP值可行。肖武等^[24]采用CASA模型对矿区NPP进行了模拟, 并将其计算结果与朱文泉的^[41]进行比对, 也在合理范围之内。本文NPP采用CASA(Carnegie-Ames-Stanford Approach)模型进行估算。

(2)

式中, N_i为第i类土地利用类型经NPP修正后的生态系统服务价值当量；P_{i_y}为研究区第y年第i类土地利用类型的单位面积NPP；y指1993年、2000年和2013年；E_0i和i同上。

2.2.3 降水因子修正

降水因子修正水源涵养服务功能。计算公式如(3)所示。

韩鹏等利用当量因子方法计算干旱半干旱区生态系统服务价值时认为降雨因素会对水源涵养价值产生影响^{[30, 42]}。

(3)

式中, P_i为第i类土地利用类型经降水因子修正后的生态系统服务价值当量；W和W₀分别为研究区和全国年平均降雨量；E_0i和i同上。

2.2.4 沉陷因子修正

采煤沉陷影响土壤保持量, 栅格尺度上土壤保持量修正土壤形成与保护服务功能^[30]。计算公式如下：

(4)

式中, s_i为第i类土地利用类型经沉陷因子修正后的土壤保持服务价值当量；S_{i_y}为研究区第y年第i类土地利用类型土壤保持量；E_0i和i同上。

2.2.5 生态系统服务价值计算

(5)

(6)

式中, ESV_f为第f项生态服务价值；ESV为总的服务价值；A_k为研究区第k类型的土地利用面积(hm²)；VC_k为第k类型土地利用单位面积的服务价值系数(元hm^-2 a^-1)；VC_fk为第k类型土地利用对应的第f项生态功能的单位面积服务价值系数(元hm^-2 a^-1)。

2.2.6 敏感性指数检验

敏感性指数的计算公式如下：

(7)

式中, ESV为生态系统服务价值, E为生态系统服务价值当量, i, j分别表示初始的生态系统服务价值和生态系统服务价值当量调整后的值。

敏感性指数检验, 意指将各类土地利用类型的价值指数分别调整50%, 衡量总生态系统服务价值的变化情况。若CS>1, 表明生态系统服务价值(ESV)相当于当量系数(E)是富有弹性的；若CS < 1, ESV则被认为是缺乏弹性的。比值越大, 表明当量系数的准确性越关键。

3 结果与分析 3.1 土地利用变化特征

1993—2013年, 耕地和建设用地为鹤岗矿区的主要用地类型, 平均比重依次为47.18%和44.85%。耕地和建设用地变化最大, 耕地依次减少, 减少量共计910.16hm², 变化率为-14.99%；建设用地依次增加, 增加量共计776.69hm², 变化率为15.98%(表 1)。经土地利用转移矩阵分析, 2013年土地利用类型的减少或增加均是1993年土地类型的得失转换。建设用地的扩张主要是由耕地的减少转化而来, 2013年1078.91hm²耕地转化为建设用地(图 2)。

3.2 生态系统服务价值分析 3.2.1 生态系统服务总价值

由1993—2013年鹤岗矿区不同生态系统服务类型的价值变化来看(表 2), 1993、2000和2013年鹤岗矿区生态系统服务的总价值量分别为2219.21、1025.15和3531.95万元。耕地的平均总服务价值最大, 约占平均总价值的189.35%；其次是林地(46.50%)、水域(26.83%)；草地、未利用地和建设用地所占比重偏低, 建设用地平均总服务价值为-3925.51万元, 约占平均总价值的-173.79%。

1993—2000年, 总价值变化量为-1194.06万元, 其中耕地变化量最大, 为-1059.36万元, 变化率为-22.44%；其次是水域, 变化价值量为-54.53万元, 变化率为-9.45%。2000—2013年, 总价值增加2506.80万元, 其中林地、草地、耕地、水域和建设用地价值均在增加, 耕地和建设用地价值增加最大, 分别为789.97和1037.68万元, 变化率最大的是林地和草地, 分别为43.94%和39.97%。1993—2013年间, 总价值增加1312.74万元, 其中建设用地增加量最大, 为1010.60万元, 变化率为-23.76%(表 2)。1993—2013年价值先减少后增加主要是随着采矿活动, 耕地面积在不断下降, 建设用地、林地和水域面积在不断增加, 而采煤沉陷地土地生态状况呈现出先下降后回升的趋势。建设用地面积增加到一定程度, 会使总服务价值提高。

3.2.2 生态系统服务类型

由1993—2013年鹤岗矿区不同服务功能的价值变化可知(表 3), 1993—2013年, 土壤形成与保护价值最大, 分别占该年总价值量的145.18%、192.97%和154.86%；其次是气候调节服务和生物多样性服务, 平均价值比分别为40.54%和37.79%；食物生产的平均价值比为22.36%；其他服务价值占比较小。由于鹤岗矿区建设用地面积较大, 水源涵养服务价值为负值, 分别占该年总服务价值量的-155.99%、-308.65%和-127.68%。

1993—2000年, 土壤形成与保护价值减少最大, 为-1243.66万元, 年变化率为-5.51%。气候调节和水源涵养服务价值均处于增加状态。2000—2013年间, 土壤形成与保护价值增加最大, 为3491.17万元, 年变化率为13.58%；其次为气体调节和废物处理服务, 增加值分别为129.34和277.79万元。1993—2013年间, 气体调节价值年变化率最大, 分别为5.70%、44.04%和42.04%(表 3)。1993—2013年, 土壤形成与保护价值先减少后增加, 主要是随煤炭开采, 耕地和建设用地的土壤保持量呈先下降后上升、建设用地面积大大增加的结果。

3.2.3 生态系统服务价值的空间差异及其影响因素

鹤岗矿区煤炭开采方向自西向东行进, 且北部煤矿开发较早于南部。1993—2013年, 生态系统服务价值西部小、东部大, 北部小、南部大, 自西向东递增。西部一带为过去各煤矿煤炭开采初始地段, 土地利用类型多由耕地转变为建设用地, 生态系统服务价值较低；东部一带为各煤矿开发较晚或还未开发地段, 土地利用类型多为耕地, 生态系统服务价值较高。鸟山煤矿2017年投产使用, 受采煤扰动较少, 生态系统服务价值最高, 占平均总价值的60.19%。北部煤矿开发程度较深, 兴山、益新、振兴和南山煤矿生态系统服务价值均值之和占平均总价值的-11.65%；南部峻德煤矿, 服务价值均值占平均总价值的48.08%, 服务单位价值最高的地区主要分布在峻德矿沉陷区范围, 地型为林地(图 3)。这种分布趋势主要受由采煤活动引起的土地利用类型、面积及质量变化的影响。

3.2.4 鹤岗矿区生态系统服务价值修正前后比较

1993—2013年, 鹤岗矿区总生态系统服务价值分别为2219.21、1025.15和3531.95万元, 与无考虑空间异质性因素情况相比, 1993和2000年服务总价值分别减少了152.96和580.12万元, 2013年增加了781.63万元。随着煤炭开采, 1993—2000年, 鹤岗矿区生态环境受采煤影响较大, 到2013年, 生态环境呈好转的趋势。

4 结论与讨论

生态系统服务与人类福祉关系极其密切。充分评估矿区沉陷地生态服务价值是生态系统服务恢复的科学依据。本文主要采用修正的单位面积生态系统价值当量因子的方法, 对鹤岗矿区生态系统提供的9种生态系统服务类型进行评估, 得出以下结论：

(1) 1993、2000和2013年鹤岗矿区总服务价值量分别为2219.21、1025.15和3531.95万元。随着采矿活动的进行, 这与耕地面积不断减少、建设用地不断增加有关, 也与采煤沉陷地土地生态指标呈现先下降后回升有关。建设用地面积增加到一定程度, 土壤形成与保护价值的大大增加会使总服务价值提高。

(2) 按照生态系统类型划分, 1993—2013年, 鹤岗矿区耕地的平均服务价值最高, 占总平均价值的189.35%；其次是林地和水域, 分别占总平均价值的46.50%和26.83%；按照生态系统服务类别划分, 支持功能服务价值最高, 平均价值占总平均价值的186.67%, 其次是调节服务, 占45.34%, 文化服务占5.73%, 供给服务占-137.74%。

(3) 1993、2000和2013年, 鹤岗矿区生态系统服务价值自西向东递增, 北部小、南部大, 与采煤方向由西向东行进、开采程度差异有关。生态系统服务价值空间差异主要受由采矿活动引起的土地利用类型、面积、质量变化的影响。

(4) 与无考虑空间异质性因素情况相比, 鹤岗矿区1993和2000年生态系统服务总价值分别减少了152.96和580.12万元, 2013年增加了781.63万元。随着煤炭开采, 1993—2000年, 鹤岗矿区生态环境受采煤影响较大, 到2013年, 自然生态环境状况好转。