生态系统弹性力描述了生态系统受到环境冲击之后保持其原有结构和功能的能力, 它对于区域生态系统安全评估具有重要意义^[1-2]。生态系统弹性力被定义为在维持系统结构、功能反馈等不变的前提下, 通过调整系统状态变量和驱动变量等参数, 系统能吸收的扰动量, 用以表示生态系统在面对外部干扰时保持生态系统功能的能力^[3-5]。Holling等^[3]将这种能力量化为一个生态系统能够承受的变化幅度。通常在受到较小的环境干扰时, 生态系统更容易恢复, 当生态系统遭受到严重的扰动时会趋向远离其均衡稳态, 进而转变为另一种状态, 最终导致系统受损而退化^[5]。在没有人为干扰的自然生态系统中, 生态弹性力也可用来反映自然生态系统承载力, 是区域生态承载力的研究基础和描述脆弱生态区状态的综合变量^[6-8]。

量化生态系统弹性力是预测生态系统应对全球变化的关键^[9]。国内外学者采用不同的方法对生态系统弹性力展开了大量的研究, Ashutosh等^[10]采用水分效率模型来评估印度的生态系统弹性力, 发现生态系统弹性力可能受土地覆盖、气候、水文要素的影响；Stephen等^[11]利用全局敏感性和不确定性分析法评估了生态系统弹性力, 为生态系统弹性力定量评价提供了有效方案；Luciano等^[12]基于气候生态位理论提出了测量生态系统弹性力的方法, 发现持续的气候变化可通过促进森林生物多样性的侵蚀并导致植被覆盖度降低而加速生态系统弹性力的丧失；Wu等^[9]用叶面积指数来评估陆地生态系统的弹性力, 并将其动态变化归因于气候与环境因素；生态弹性力模型^[13-14]结合主成分分析、模糊数学的综合指标评估^[15-16]方案目前也被广泛应用。对于生态系统弹性力的影响因素, 已有研究表明生态系统弹性力与植被生产力、生物多样性呈正相关关系^[17]。此外, 气候变化和人类活动的累计效应影响着生态系统的结构和动态, 进而显著影响生态系统的恢复力^[18]。虽然大量学者对生态系统弹性力的影响因素进行了探究, 但目前仍缺少一致的观点和表征生态系统弹性力的指标体系^[19]。总的来说, 现有研究方法多应用于区域尺度上, 且考虑的因素还不够全面, 针对工矿活动下的流域生态系统弹性力的研究还鲜有报道。

陕北榆神矿区作为我国重要的煤炭生产基地之一^[20], 煤炭资源丰富, 人类活动剧烈, 持续的工矿活动引发了矿区局部水土流失^[21]、植被退化^[22]以及地下水位下降^[23-24]等一系列生态环境问题, 导致该地区生态环境更加脆弱^[25]。自1999年起, 该地区退耕还林还草、天然林保护和封山禁牧等生态修复工程的实施虽然整体上改善了区域的生态环境^[26-27], 但工矿活动的直接影响和间接影响区内的生态环境变化特征尚不清晰, 同时对近年来工矿影响区的生态系统稳定性状况缺乏系统的了解。煤炭开发区生态环境复杂, 影响其生态系统稳定性因素较多, 如何找准生态影响因子, 客观合理地评价矿区生态现状及其恢复力显得尤为重要。本文以榆神矿区所在流域为研究对象, 在了解研究区整体土地利用和煤炭矿区分布格局的基础上, 依据研究区生态环境背景, 选取水文、土壤、植被3个方面5个指标, 运用综合指标法和主成分分析法对研究区生态弹性力进行评估, 探讨研究区2009—2018年生态系统弹性力时空演变特征, 分析煤炭开发对研究区生态系统弹性力的影响。本研究通过对矿区生态系统弹性力的评估来量化矿区生态系统稳定性和可调节能力, 以期为优化矿区生态格局和进一步完善矿区可持续发展策略提供科学依据。

1 研究区概况

榆神矿区位于陕晋蒙接壤地带, 是陕北煤炭基地的主要矿区之一, 矿区面积约5265 km², 2013年, 榆神矿区煤炭产量已达到2.78亿t^[24]。矿区总体被划分为四个规划区, 根据不同的资源勘探程度, 开启了前三期规划区的开发, 第四期规划区正在准备阶段^[23]。本文以榆神矿区所在的3个流域为研究区(图 1), 包含秃尾河流域(占流域总面积44.7%)、榆溪河流域(占45%)和佳芦河流域(占10.3%), 研究区面积约为1.07万km², 处于东经109°—110.52°, 北纬37.98°—39.34°之间。研究区东南部为黄土丘陵沟壑区, 西北部为风沙草滩区, 地势由东南向西北递增, 主要植被有沙柳、柠条、沙蒿、锦鸡儿、苜蓿等天然植被和柳树、杨树、枣树等人工植被, 主要土壤类型为风沙土、黄绵土和冲积土^[23]。该区地处西北内陆干旱区, 降水年际变化大, 多年年均降水约400 mm, 水资源量短缺, 多年年均潜在蒸发量约1700 mm；区域内年均气温在6—14℃之间, 多年均温约10℃^[28]。由于研究区气候干燥、植被稀疏、地形坡度大等自然地理特征使其生态环境脆弱, 煤炭资源的开采, 导致研究区地表植被破坏、水土流失、土地荒漠化程度加剧、地下水位下降、地面塌陷等生态环境问题的产生和发展^[4], 严重影响该区域生态系统的稳定性^[29]。

2 数据与方法 2.1 数据来源与处理

本研究选取2009、2012、2015和2018年4期Landsat TM 30 m分辨率的遥感影像数据(行列号127/33、34), 数据都集中于6月份, 且云量小于5%。对遥感影像进行辐射定标、大气校正、波段融合、影像镶嵌和裁剪等预处理工作, 依据全国生态系统分类体系^[30]和研究区土地利用类型的实际情况, 将研究区的土地利用类型划分为耕地、草地、林地、建设用地、水域、未利用地和工矿用地7种。通过人机交互式解译获得研究区土地利用类型空间分布图, 最后利用目视解译结果与整体分类结果建立混淆矩阵, 计算Kappa系数来验证分类精度, 四期土地利用分类精度达80%。

为了评价研究区的生态系统弹性力, 本文收集了研究区2009—2018年生态环境因子时空数据。其中地表径流和地下径流数据基于SWAT(Soil and Water Assessment Tool)模型在水文响应单元上的模拟结果, 利用地理信息系统进行数据格式转换, 得到2009—2018年分辨率1 km的地表、地下径流空间分布图。植被覆盖度是SPOT5-VGT2卫星采集的旬值数据, 空间分辨率为1 km；植被生产力(NPP)数据来源于美国航空航天局(NASA)(https://earthdata.nasa.gov/)提供的MOD17A3数据产品, 空间分辨率为1 km。土壤保持量基于通用土壤流失方程(revised universal soil loss equation, RUSLE)计算, 利用地形、土壤数据并结合月降雨量、植被覆盖度计算得到, 其中地形数据采用地理空间数据云SRTM 90 m数字高程模型(DEM)产品(http://www.gscloud.cn/), 土壤数据来源于寒区旱区科学数据中心(http://westdc.westgis.ac.cn)提供的世界土壤数据库, 空间分辨率为1 km, 所需数据包括土壤机械组成和土壤有机碳数据。

2.2 评价方法 2.2.1 生态弹性力指标体系构建

煤炭矿区生态弹性力反映了生态系统对矿区社会经济活动压力的恢复能力^[31], 根据矿区实际生态环境状况, 选取水文、土壤和植被三个方面的指标(表 1), 通过构建生态弹性力评价指标体系, 结合空间主成分分析法量化研究区的生态弹性力。根据研究目标, 按照目标分层法建立指标体系, 其中生态弹性力作为目标层, 水文、土壤和植被三个要素作为准则层, 基于准则层选取相关指标作为指标层。具体步骤如下：

① 为解决各个指标量纲不同的问题, 首先需要对原始数据进行归一化处理；

(1)

式中, Y_i为指标标准值；X_i为第i项指标的值；max(X_i)与min(X_i)为第i项指标的最大值和最小值, i的范围为1—n, n表示栅格个数。

② 根据归一化数据判定指标之间的相关性, 建立相关系数矩阵, 指标之间的相关系数越大表明有必要对数据进行主成分分析；

③ 计算相关系数矩阵的特征值、主成分贡献率和累计贡献率, 确定各主成分个数, 并计算主成分权重；

(2)

式中, G_j为主成分j的权重, λ为特征值的方差贡献率, m为主成分个数(一般特征值λ_m的累计贡献率达80%—95%)；

④ 计算生态弹性力指数；

(3)

式中, Z为生态弹性力指数, F_j为主成分j的综合得分。

2.2.2 水文要素计算

SWAT模型是对流域尺度水文过程进行连续时段模拟的物理模型, 能够模拟多种土地利用、土壤和管理条件下的水文过程^[32]。本文采用该模型对研究区的降雨径流进行模拟, 使用逐月径流对模型进行率定和验证。选用决定系数R²、百分数偏差PBIAS评价SWAT模型在率定期和验证期的表现, 具体的计算公式参照已有文献[33]。基于Moriasi等^[34]人2007年提出的标准, 两个水文站点的模拟结果在率定期和验证期凭借R²>0.50, |PBIAS|<25%被判定为良好(表 2)。

2.2.3 土壤保持量计算

土壤保持是指生态系统防止水土流失的侵蚀调控及对泥沙的储积保持能力, 可以用一个地区的土壤保持量来衡量土壤保持能力^[35]。目前往往采用修正后的通用土壤流失方程评估土壤保持能力^[36], 认为潜在土壤侵蚀量与实际侵蚀量的差值是在植被和人为管理作用下的土壤保持量, 其计算公式为：

(4)

(5)

(6)

式中, SR为土壤保持量(t hm^-2a^-1)；AP_n为潜在侵蚀量；AR_n为在植被覆盖和水土保持措施下发生的实际侵蚀量；R为降雨侵蚀力因子(MJ mm h^-1 hm^-2 a^-1)；K为土壤可蚀性因子(t hm² h hm^-2 MJ^-1 mm^-1)；LS为坡长坡度因子；C为植被覆盖与管理因子(0—1)；P为水土保持措施因子(0—1)。

降雨侵蚀力因子(R)可以反映降雨引起土壤侵蚀的潜在能力^[37]。本研究采用年降雨侵蚀力模型计算降雨侵蚀力因子^[38], 其计算公式为：

(7)

式中, Pre_i为月降雨量(mm), Pre为年总降雨量(mm)。

土壤可侵蚀性因子(K)反映了土壤的抗侵蚀能力, 其大小与土壤结构和有机碳含量密切相关, 具体计算公式^[37]如下：

(8)

式中, sa为沙粒含量(%), si为粉粒含量(%), cl为粘粒含量(%), c为有机碳含量(%), sn=1-sa/100。

植被覆盖因子(C)反映的是植被覆盖及相关管理措施对土壤的综合作用, 主要受植被覆盖度和土地利用类型等因素的影响, C值越小, 表明地面植被覆盖度越高。其计算公式^[39]为：

(9)

式中, f_c为植被覆盖度。

坡长坡度因子(LS)是影响坡面土壤侵蚀的重要因素, 其计算公式^[39]如下：

(10)

(11)

(12)

式中, λ为坡长(m), θ为坡度, n为坡长指数, S为坡度指数。

水土保持措施因子(P)是特定的水土保持措施下土壤流失量与未实施水土保持措施时的土壤流失量比值, 研究发现水土保持措施因子与坡度存在线性关系^[40], 其计算公式为：

(13)

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

研究区土地利用类型以林地、草地、未利用地为主(图 2、表 3), 耕地、林地和草地三者面积之和约占整个研究区的70%, 未利用地面积约占研究区的30%。研究期间, 未利用地面积不断减少, 减少了780.2 km², 主要表现在研究区的中部和西北部；其他土地利用类型都表现出了不同程度的增加, 其中草地、林地增加幅度最大, 分别增加了351.2、217.7 km²；建设用地增加了59.1%, 主要体现在榆阳区和神木县市区。对比2009年和2018年土地利用变化, 土地利用转移主要发生耕地、林地、草地和未利用地之间, 同时, 89 km²的未利用地和112 km²的草地转换成建设用地, 3.4 km²的未利用地和4.4 km²的草地转换成工矿用地。总体来看, 随着生态恢复工程的实施和社会经济的快速发展, 研究区土地利用空间格局发生了很大的变化, 生态环境明显改善, 城市发生了进一步的扩张^[41]。

研究期间, 工矿用地面积明显增加(图 3), 从8.4 km²增加到14.2 km²。为了解煤炭开发对矿区生态环境的影响, 以工矿用地为中心分别设置半径1 km和3 km的缓冲区, 1 km缓冲区作为煤炭开发的直接影响区, 1—3 km缓冲区作为间接影响区。2009年、2012年、2015年、2018年直接影响区分别占整个区域面积的14.8%、14.7%、19.6%、24%, 间接影响区面积也逐渐增加, 约占研究区面积的45%。

3.2 生态弹性力评价 3.2.1 生态弹性力评价指标分析

经过分析, 得到了研究区2009—2018年生态弹性力评价指标时间趋势(图 4)和空间分布(图 5)。

2009—2018年, 研究区地表、地下径流量均呈现显著增加的趋势, 增加速率分别为6.27、8.36 mm/a。对于长期来说, 随着植被的增加研究区径流量呈下降趋势, 但是近年来, 由于降雨量的显著增加导致径流量增加, 这与实测资料分析结果一致^[42]；此外, 在植被增加的同时研究区建设用地面积也增加, 通过统计不同土地利用对应的地表径流量(表 4), 发现林、草地仅是建设用地地表径流量的50%左右。从径流量在水文响应单元尺度上的空间分布特征来看, 研究区地表径流量主要处于6—30 mm之间, 在研究区北部红碱淖湖周围地表径流量最小(< 6.5 mm), 总体上呈现北部低西南部高的空间分布特征。研究区地下径流量较地表径流量大, 主要由于研究区大面积为风沙区, 地表下渗量大；地下径流量在空间上表现出东北部高、西部偏低的特征, 东南部佳芦河流域地下径流量最小, 主要处于29—70 mm之间。

2009—2018年土壤保持量总体上呈现增加趋势, 局部波动幅度较大, 主要归因于降雨侵蚀力因子, 说明在该地区降雨是土壤侵蚀的主导因素。空间上, 50%以上的研究区域土壤保持量小于10 t/hm², 呈现东南部明显高于西北部的特征, 由于西北部属于流域上游地势高、植被覆盖低且土壤类型多为砂土, 土壤易发生侵蚀, 致使土壤保持量较低。

研究期间, 植被覆盖情况显著改善, 增加速率为0.63%/a；植被生产力变化趋势基本与植被覆盖度一致, 其增加速率为6.85 gC/m²。空间上, 研究区植被覆盖度处在0—0.29之间, 呈现出东南部明显高于西北部的特征。研究区植被生产力多年均值处在8—355 gC/m²之间, 东南部和北部林草覆盖较多的区域植被生产力较高, 且沿着河道的耕地区植被生产力最高(> 200 gC/m²), 由于研究区的西北部主要为未利用地, 植被生产力也处于较低的水平, 主要集中在28—140 gC/m²之间。

3.2.2 生态系统弹性力分析

根据生态弹性力评价指标间存在相关系数大于0.5的结果得知, 指标之间存在一定的相关性, 适合采用主成分分析法确定权重。依据选取原则确定了两个主成分, 其权重如表 5所示。根据主成分荷载, 主成分一主要支配指标为地下径流量、植被覆盖度, 主成分二主要支配指标为地表径流量、植被生产力。

研究区2009—2018年生态弹性力的时间变化趋势如图 6所示, 生态弹性力指数总体呈现增加趋势, 从2009年的0.53增加到2018年的0.60, 但中间过程波动幅度较大, 主要受径流量的影响, 这是由于近年来陕北地区降雨量显著增加且年际变化大。2009—2017年, 生态弹性力呈现“W”型变化特征, 其中2009—2011年, 弹性力由0.53增加至0.56, 基本呈现平稳趋势, 2012年后弹性力呈现先增加后减小的变化趋势, 2015年降低至研究期间的最低值, 2015—2017年弹性力指数显著增加, 2017年生态弹性力达到最大值0.77, 2017—2018年弹性力又表现出明显的减小。总体上, 研究区生态系统抗干扰和自我调节能力不断提高但不稳定, 与郑欣等^[16]在鄂尔多斯的研究结果相似。

从各个主成分来看, 第一主成分与植被覆盖度和地下径流量都存在较高的正相关关系。由于第一主成分的贡献率较大, 其变化趋势与生态弹性力基本一致。因此, 第一主成分对生态弹性力起决定性作用, 而植被覆盖度和地下径流量主要决定了研究区生态弹性力的变化趋势。从图 6中可以看出, 主成分二减缓了生态弹性力的波动幅度, 而主成分二的变化主要取决于地表径流量和植被生产力。因此, 地表径流量和植被生产力对减缓研究区生态系统稳定性的变化起到正向作用。分析表明, 本研究区生态弹性力对水文、植被条件变化比较敏感, 主要受植被覆盖度和地下径流量的影响。

由于考虑了多个空间指标, 使得生态弹性力空间分布差异较大(图 6), 总体呈现东南部偏高西北部低的特征, 表明研究区东南部的生态系统较西北部稳定。由于研究区东南部植被覆盖度大, 又属于流域的下游, 水资源相对丰富, 生态弹性力较高, 而西北部受气候和土壤条件影响, 裸地分布较广, 植被覆盖度和植被生产力都较低, 因此生态弹性力较低。局部来看, 在榆阳区和神木市县区周围及沿线生态弹性力较低, 主要归因于城市建成区人类活动频繁, 生态环境较差。

3.2.3 榆神矿区生态弹性力变化

统计不同土地利用对应的生态弹性力(图 7), 得到不同土地利用生态弹性力排序：林地>草地>耕地>建设用地>工矿用地, 由于林地的植被生产力、植被覆盖度高且水土保持能力强, 因此生态系统稳定性较高。2009—2015年, 生态弹性力呈现下降趋势, 林草地的弹性力下降最为明显。虽然近年来随着生态修复措施的实施, 研究区植被覆盖度和植被生产力增加, 但生态弹性力受水文要素的影响很大, 由于降雨量的降低使得径流量减少, 导致生态弹性力降低。

煤炭矿区直接影响区是煤炭开发的内围区域, 生态系统弹性力较间接影响区低, 平均低9.23%。其中2009年两个影响区生态系统弹性力差距最大, 间接影响区比直接影响区弹性力高0.07, 2015年差距最小, 整个区域的弹性力介于直接影响区和间接影响区之间(图 8)。总体而言, 由于耕地的植被覆盖度和植被生产力较高, 生态系统弹性力也较高, 而直接影响区中耕地面积比例较大, 导致直接影响区的生态弹性力偏高。由于间接影响区中有很大比例是裸地, 裸地对应的生态弹性力低。因此, 间接影响区相比直接影响区的弹性力虽有增加, 但差距不明显。

4 结论与讨论 4.1 讨论 4.1.1 煤矿开采对矿区生态系统的影响

随着矿区煤炭基地的建成, 工矿用地侵占草地和灌丛面积, 加上生态环保政策落实不到位, 加剧了矿区植被退化。由于工矿用地人类活动频繁, 煤矿开采直接影响区生态系统稳定性较差。大规模开采活动可能造成矿区地面塌陷、地裂缝, 导致植物根系直接拉断, 并且土壤物理结构的破坏和水分养分的流失也影响植物生长, 造成开采区植被退化和生产力降低。随着采煤面积的增大可能会引起地裂隙加大、加深, 地表水通过裂缝和破碎区连通的位置渗入地下, 造成地表水流失；同时裂隙可能穿透上部隔水层, 从而造成地下水渗漏, 地下水位下降。在煤矿开采过程涉及和影响的土地范围内, 原生森林植被生态系统、土壤结构特征以及其天然的水土保持功能被严重削弱, 表层土壤和植被的水土保持能力大幅下降。因此, 在煤炭开采区应及时采取生态修复治理措施, 降低采矿活动引起的地质环境问题, 从而提高矿区生态系统弹性力。

4.1.2 生态弹性力指标选取合理性分析

本文从水文、土壤、植被3个方面选取指标综合评估榆神矿区生态系统弹性力, 而气候条件也是决定生态弹性力的主要因素。通常在区域尺度较大的情况下评估弹性力需要考虑气候要素, 而本研究区气象要素空间差异较小, 因此没有直接选取气象指标。由于矿区本身属于水资源缺乏地区, 水资源的供给能力是生态弹性力提高的主要限制因素, 选用了地表、地下径流两个水文指标。考虑到主成分分析时各指标间需要存在一定的相关性, 选用了土壤保持量、植被覆盖度和生产力等相关指标来反映研究区生态环境变化。植被覆盖度和生产力均可表征生态系统的自我恢复能力。植被覆盖度对区域生态环境变化有重要指示作用, 能够反映区域的水土保持、气候调节以及生态系统的稳定性强弱^[43], 多项研究直接从植被覆盖度、叶面积指数等遥感数据中检测植被异常以测量生态系统的恢复力^{[9, 44]}。而植被生产力可表征一个区域地形地貌、水分热力状况以及植被的生产适宜性, 可通过植被生产力来客观反应生态系统各要素的健康程度^{[43, 45]}, Frazier等利用遥感反演的NPP作为生态系统健康指标监测生态系统恢复情况^[46]。植被灌层通过拦截降雨、减少雨滴击溅、阻挡泥沙运输来减少水土流失, 不同的植被和类型都会影响一个区域的水土保持能力^[47], 而对于陕北黄土沟壑地区, 有研究表明降雨是影响水土流失的主要因素, 植被覆盖与土壤侵蚀的关系不明显^[48]。同时土壤组成也是土壤保持量的决定性因素, 研究发现不同性质的土壤生态弹性力存在很大差异^[49]。

4.1.3 优势与不足

利用空间二维栅格数据进行空间主成分分析来评估生态弹性力, 高度保留了各指标的空间信息, 从而可以反映区域生态弹性力的空间差异, 便于识别生态问题严重的区域, 可为生态环境修复提供参考依据。以往的研究大多考虑了生态弹性力的时间趋势, 反映的是区域整体生态弹性力相对变化情况。本研究方法更能反映生态系统弹性力区域空间格局情况。

由于数据分辨率以及指标获取等限制, 本研究仍存在一定局限性：数据的分辨率较低导致综合评估得到的生态弹性力空间细节描述较弱；水文模型校准时采用的观测数据年份较早, 且在研究区的模拟效果只达到了满意, 得到的径流量存在不确定性；用于表征生态环境变化的因素还有很多, 本研究中由于数据可获取性问题, 所选取的指标不够全面, 我们认为尽可能全面的评估对矿区开发规划更具有指导意义, 但本研究在区域生态环境健康评价方面提供了一套有效的思路和方法。

4.2 结论

本研究通过遥感解译获得了榆神矿区土地利用分布格局, 在充分明确矿区生态环境情况的基础上评估了榆神矿区生态系统稳定性, 分析了煤炭开发对矿区生态环境质量的影响程度。主要结论：

近10年来, 随着生态修复措施的实施和社会经济的快速发展, 研究区生态环境明显改善, 林地草地面积明显增加, 但是煤炭开发用地也大幅增加, 从8.4 km²增加至14.2 km²。

研究期间, 榆神矿区生态系统抗干扰能力和可调节能力总体不断提高但伴有波动, 生态系统弹性力从2009的0.53增加至2018年的0.60, 空间上呈现东南部高西北部低的特征, 该区生态系统弹性力主要受植被覆盖度、地下径流量的影响。

煤炭开采活动引发的地面塌陷、地裂缝等地质环境问题致使煤矿开采区植被退化、地下径流、地表径流量降低以及水土保持能力下降, 导致煤炭开采直接影响区的生态弹性力低于间接影响区, 平均低9.23%。因此, 及时采取修复治理措施降低煤炭开发引起的地质环境问题是提高矿区生态弹性力有效途径。