基于恢复力的半干旱矿区生态系统退化风险空间评估

2021-06-18吴秦豫张绍良杨永均侯湖平陈东兴

煤炭学报 2021年5期

吴秦豫，张绍良，杨永均，侯湖平,陈东兴

(1.中国矿业大学公共管理学院，江苏徐州 221116; 2.矿山生态修复教育部工程研究中心，江苏徐州 221116; 3.山东省采煤塌陷地与采空区治理工程研究中心，山东济宁 272100)

2020年，中共中央政治局召开会议审议《黄河流域生态保护与高质量发展规划纲要》，指出黄河流域生态保护与高质量发展是事关中华民族伟大复兴的千秋大计，其中黄河流域中段的晋陕蒙煤炭基地众多，面积达数千平方千米。煤炭开采导致大面积地表沉陷和地裂缝、地下水位下降、泉水干涸、河川基流量衰减和表生生态环境恶化等一系列环境问题[1]，生态保护和退化问题受到广泛重视，现有的技术和方法还不能支撑黄河流域煤矿开采与生态环境的协调发展[2]，亟待实施空间精准的生态退化控制与保护。

针对生态脆弱区大型煤炭基地的生态退化控制与保护问题，已经有较多学者取得了研究成果。在煤炭开采的生态影响方面:有学者分析了不同地域矿区煤炭开采引起的生态效应及地域分异现象[1];LIU等[3]研究了自然条件和采矿活动对干旱半干旱矿区植被变化的影响，发现地下水埋深是矿区植被生长的关键控制因素;侯湖平等[4]对煤矿区的植被碳储量进行测算，研究发现矿区碳储量的平均值变化量与采矿生产能力呈现高度负相关;在生态风险控制方面:主要针对水这一关键生态要素，提出了控制方法。如，范立民[5]研究了榆神矿区工程地质条件分区，为保水采煤方法选择提供了地质基础依据;李涛等[6]提出了基于采动分水岭的地表水保水煤柱留设方法和采煤对地表径流的影响方式，认为在重点水体处应合理留设煤柱保护地表水。在生态恢复方面:胡振琪[7]在分析采前地面地形特征的基础上，结合煤矿开采计划，对地表动态沉陷过程进行了模拟，分析煤炭开采影响下地面土地利用变化的动态过程，以确定复垦边界，实现复垦空间精准化;毕银丽等[8]提出了通过接种丛枝根菌真菌提高煤矸石山复垦土地上的植被成活率;郭洋楠等[9]提出基于黄土沟壑区地形地貌建立沟道防护林，在沟道自上而下分为侵蚀区、流过区、沉积区3个不同地段，配合沟道内的工程防护措施，分别布置不同林木保证矿区治理工程措施的正常运行。

然而，已有研究成果主要是针对水、土壤等单一生态要素展开的。地表生态系统的状态实际上是由环境因素和非环境因素共同决定的。这些关键因素作用结果体现在生态系统的恢复力，即生态系统面临扰动时保存其状态的能力[10]。目前，矿区生态恢复力研究已有一些探索。在恢复力理论方面，有学者指出研究复垦区域的生态系统问题，需要考虑生态系统稳定性(包括阻力、恢复力和持续力)，恢复力可用系统受扰动后恢复到平衡态所需时间表示[11]。也有学者提出恢复力在恢复半干旱矿山生态系统和创建可持续生态系统具有重要作用[12]。实现矿业社会-生态系统的可持续发展，需要从恢复力视角管理扰动和脆弱性[13]。在恢复力评估方面，有一些学者针对矿区生态系统恢复力开展了实证研究，如以植被覆盖度作为弹性分值，并结合土地利用数据对矿区生态系统弹性进行评估[14]，也有学者在区域尺度上利用多因素加权方法解释了矿区对自然与人为扰动的恢复力的空间差异[15]。综合来看，构建一个自维持、具有恢复力的生态系统已经成为矿山生态修复实践的普遍需求。然而，现有的矿区生态保护修复规划还较少考虑恢复力。笔者通过定量评估恢复力，进而识别矿山生态退化风险，有助于在矿区尺度上制定精准的生态保护与退化控制措施。

为此，笔者基于改进的突变级数法，考虑生态系统多个关键参数及其作用关系，开展恢复力空间评估，并基于不同生态系统类型的恢复力指数划分退化风险等级，为黄河流域煤炭基地生态系统保护提供决策依据。

1 研究区域与研究方法

1.1 研究区概况

陕西省神木市位于黄河流域中段，晋陕蒙3省交界处，东临黄河，地理位置为东经109.62°～101.93°，北纬38.20°～39.45°。总面积7 470 km2。

神木市地理区位表现出明显的过渡性。以明长城为界，明长城以南主要为丘陵沟壑区，明长城以北主要为风沙草滩区。市内拥有丰富的地形地貌类型，包括沙地、丘陵、高原、草滩、湖泊等(图1)。市内年均气温8.6 ℃，年均降水量425.5 mm，年际差异很大。70%的降水集中在7—9月，属于半干旱气候。

图1 神木市的地理区位Fig.1 Geographical location of Shenmu City

神木市境的河流有黄河、窟野河和由流入红碱淖河流组成的内陆水系。窟野河与秃尾河以黄河峡谷为其侵蚀基准，在新构造上升的配合下，河流下切剧烈，有些河段已切入基岩。黄河在地质因素控制下，河床切入三叠系。在黄河及水土流失等因素影响下，河流多泥沙。

神木市土壤类型以风沙土和黄绵土为主，风沙土主要分布在县域北部及西北部的风沙草滩区，土壤结构松散，受风蚀影响成土过程不稳定，容易形成沙化。该地区的植被以沙生植物和隐喻性植被为主，主要有沙蒿、沙柳、杨柴等。黄绵土主要分布在县域东南部黄土丘陵地区，土质绵软松散，有机质含量低，容易受水蚀和风蚀影响，土壤侵蚀情况严重。该区域多分布草本植物，优势种主要有沙葱、早熟禾、尔泰狗娃花等，生态脆弱。

神木市煤炭主要分布在县境西北部，储煤面积达4 500 km2，占全县总面积的59%。其中富煤区的煤炭储量超过1 000万t/km2，探明储量450亿t。2019年底，全市共有各类煤矿115处，设计生产能力2.95亿t/a，形成产能2.55亿t/a。

1.2 研究方法

1.2.1生态系统恢复力理论

恢复力(Resilience)的起源是拉丁文resilio，表示再次跳回的动作。自20世纪70年代以来，恢复力的概念被延伸为系统在压力下恢复到最初原始状态的能力。HOLLING[16]最先将恢复力引入生态学领域，他对恢复力做出的定义为:系统吸纳状态变量、驱动变量和参数变化，并持续存在的能力。

恢复力是矿山生态系统的基本属性，使之在有限强度的扰动后，系统具有的恢复到平衡点、保持原状态和定性结构的能力。恢复力的产生原因是系统具有复杂结构，且是自组织的，在面临扰动时表现出自维持能力。参数变量决定了恢复力大小，其阈值空间是测度恢复力的直接指标[10]。

当生态系统遭受的扰动超过了一定阈值，任何生态系统都可能受到破坏而进入另一种不理想或者不可逆的状态。这种生态系统结构和功能发生的与原有平衡状态或进化方向相反的跃迁，就是生态系统的退化。因此，恢复力为评估生态系统的退化风险奠定了基础。

1.2.2生态系统恢复力评价指标

生态系统包括地形、土壤、植被和水文等环境要素，这些要素决定了生态系统的功能、结构与恢复力，因此本文从中选择有效指标来刻画恢复力。根据生态系统恢复力理论，结合神木市的实际情况，综合考虑生态系统多变量之间的相互控制关系，笔者选取的指标见表1。

表1 生态系统恢复力评价指标含义、表达因素与数据来源Table 1 Meaning,expression factors and data sources of evaluation indicators for ecosystem resilience

坡向、坡度是地形条件的表达指标。研究区沟壑纵横，植被盖度低，气候干燥地表蒸发量大，煤矿开采引起的地表沉陷会导致地表垂直变形和水平变形，从而改变地表的坡向和坡度。坡向通过影响光照、温度及土壤养分进而影响群落植物种的组成和分布。坡度则决定了水土流失的程度，从而影响土壤的质量。进一步影响生态系统在扰动后恢复的能力。

土壤有机质质量分数对植被的生长起着最直接的作用，有机质质量分数可作为土壤质量直接有效的表达指标。神木市南部受到地形限制，人口耕作集中于河谷川道，长期培肥使土壤有机质质量分数提升，而被不稳定风沙土所覆盖的丘陵区土壤有机质的质量分数则较低。北部大规模煤炭资源开发、运输等活动也导致周围土壤受到不同程度的重金属累积性污染，可改变土壤有机质的矿化速率，影响土壤有机质的积累与分配。进而影响生态系统恢复力。

水文条件是半干旱地区生态系统最大的制约因素。神木市风沙区地形平坦不易形成地表径流，河谷川地的汇水能力则较强，土壤水分含量较高。煤矿开采导致的地裂缝和挖损引起地下水位下降，土壤湿度减小。土壤水分条件对植被养分吸收具有显著影响，因此选择土壤湿度指数作为水文条件的表达指标。

植被盖度和多样性可以作为植被条件的表达指标。采煤塌陷直接导致植被景观被破坏，导致地表变形而产生地下水位改变、土壤侵蚀、水土流失，间接阻碍植被对水分和养分的吸收，破坏植被生长。生态系统应对环境变化的恢复力是由其生物和自然资源决定的，随着植被盖度的下降，区域维持水分存储和养分循环的能力降低，恢复力水平就越低。而植被多样性则意味着可持续的生态系统包含着多种功能群，为生态系统提供一定程度的功能冗余。当生态系统面临物种丧失或环境变化时，功能冗余通过维持生产力的方式提供恢复力。多样性可用Simpson指数表示。

1.2.3数据标准化

为了消除量纲与量纲单位的影响，在决策与排序之前，应首先将评价指标进行无量纲化处理[17]，即数据标准化。植被种类越多、覆盖度越大、土壤有机质和水分含量越高的生态系统抵御外界胁迫的能力越大，地形起伏度越小的生态系统抵抗外界干扰的能力也越大，因此，可将恢复力替代指标分为“效益型”与“成本型”，其标准化方法分别为

式中,Cjmax,Cjmin分别为第j个指标的最大值和最小值;Ci-j为第j个指标的第i个属性值。

1.2.4均方差决策法

突变级数法避免了其他评价方法需要主观确定权重的问题，但仍需对指标给予重要性排序，笔者采用均方差决策法确定恢复力指标的权重对突变级数法进行改进。均方差决策法是一种客观赋权方法，不依赖人的主观判断，具有概念清晰、计算简便的特点。设多指标综合评价问题中方案集合为A={A1,A2,…,An}，指标集为G={G1,G2,…,Gn}。若指标Gn对所有决策方案而言均无差别，则其对方案决策和排序不起作用，可令其权系数为0，若指标Gj使所有决策方案的属性值存在较大差异，则其对方案决策与排序将起重要作用，应给予较大的权数。即在多指标决策与排序的情况下，各指标相对权重系数的大小取决于在该指标下各方案属性值的相对离散程度，离散程度越大，权系数越大。均方差决策法确定评价指标离散程度最常用的指标是均方差。首先求出各方案在各指标下无量纲化属性值的均方差，公式为

(3)

(4)

其中，Pi-j为网格i中指标j的比例；k=1/lnn满足ej≥0，x′i-j为网格i中指标j的标准化值；ej为指标j的熵值。得到所有网格对第j个指标的总贡献值，将其归一化，结果即为各指标的权重系数。

1.2.5突变级数法

采用突变级数法评价生态系统恢复力。突变理论研究系统在不同稳态之间转化的现象与规律，是将系统的质变用数学模型描述出来，帮助人们认识并理解系统发生变化或者突然中断的现象[18]。突变级数法是在突变理论的基础上，与模糊数学相结合产生模糊隶属函数，从而进行多目标评价的一种研究方法。

由于生态系统的多稳态机制，任何外部干扰都可能导致系统状态的突变，进入到管理者所不希望的状态[19]，生态系统状态是由多种因素相互控制关系决定的，恢复力描述了系统在面临扰动时保存这种综合状态的能力。突变级数法可以定量刻画多因素对系统综合状态的作用，因而，突变级数法适用于评价生态系统恢复力。根据恢复力替代指标的变化特点，选择3种突变模型，分别为尖点突变系统，燕尾突变系统与蝴蝶突变系统，模型结构见表2，每个指标的突变模型见表3。

表3 生态系统恢复力评价的突变模型Table 3 Catastrophe model of ecological resilience evaluation

表2中，f(x)是系统中的一个状态变量势函数;a,b,c,d为势函数的系数，表示该状态变量的控制变量。2者是相互矛盾的，但是2者相互作用相互牵制。

表2 突变模型结构Table 2 Structure and icon of mutation model

由于x,a,b,c和d取值不统一，也为了计算方便，需要将状态变量与控制变量的取值限制在0～1，通过对分叉集进行变化和推导得到归一公式。归一公式最终归统一参数同一状态。

尖点突变归一化公式为

xa=a1/2,xb=b1/3

燕尾突变归一化公式为

xa=a1/2,xb=b1/3,xc=c1/4

蝴蝶突变归一化公式为

xa=a1/2,xb=b1/3,xc=c1/4,xd=d1/5

恢复力指数计算规则:若变量之间存在相关性，则遵循计算指标平均值的“互补原则”，否则遵循取最小指标值的“非互补原则”。随后还需要自下而上进行递归运算，以隶属函数的综合值进行全面的评估与分析，取其平均值作为恢复力指数。

2 结果与分析

2.1评价指标的空间格局

本研究基础数据空间分辨率不一，为了相互匹配，借助ArcGIS的空间重采样方法获得1km分辨率的指标集，虽然最小单元尺寸的增加伴随着平均作用，但是从总体空间格局来看，这种平均作用所带来的误差是可以接受的。

神木市的土壤有机质质量分数高值分布在北部沙漠草滩区，土壤有机质质量分数普遍高于0.8%，其次是北部土山石区，土壤有机质质量分数在0.5%左右，低值分布在中部的黄土丘陵区。地统计分析发现，神木市土壤有机质分布属正偏态分布，偏度为1.45，有机质质量分数低的区域分布广。

神木市植被盖度高的区域位于中部的丘陵沟壑区及南部的沙漠草滩区，北部的土山石区较低。地统计分析发现，研究区地表植被盖度的均值为60.07%，分布为正偏态，偏度为0.44，低覆盖区域面积大。覆盖度大于50%的区域仅为14.8%，大部分位于最南部的草滩区。该区地势起伏不大，滩地地下水与地表水补给充足，是神木市农牧业的集中区域。

神木市植物多样性高值区分布在中部黄土丘陵沟壑区和南部河谷川道区，低值处位于北部沙漠草滩区。地统计分析发现，研究区陆地生态系统每平方米Simpson指数的均值为0.44，属于负偏态，偏度为-0.32，说明多样性较高的区域所占面积较大。

神木市的水文条件通过土壤湿度指数表达，高值区主要分布在河流、湿地以及汇水能力较强的沟谷地带。通过地统计分析发现，研究区土壤湿度指数分布属正偏态分布，偏度为1.36，土壤湿度指数低的区域分布较广。

神木市东南部土山石区山大沟深，地面倾斜度大，沿河谷两岸地形狭窄、山体陡峭，山体朝向多为北、西北、西。中部丘陵沟壑区梁多峁少。梁面呈鱼脊形，以10°～20°向两侧沟谷倾斜。沟边缘线以下谷坡陡峭。北部沙漠草滩区地势较为平坦。

各评价指标的空间格局如图2所示。

图2 研究区生态指标的分布格局Fig.2 Distribution pattern of ecological indicators in the study area

2.2生态系统恢复力空间分布

利用均方差决策法分别计算表4中指标的权重，准则层中，植被、水文、地形、土壤均方差分别为0.47，0.33，0.18，0.02。地形指标中，坡向和坡度的权重分别为0.72，0.28，植被指标中，多样性和植被盖度分别为0.63，0.37。据此，利用突变级数法的归一化公式和“平均值原则”计算得到研究区生态系统恢复力指数，结果如图3(a)所示。

从图3(a)可以看出，神木市恢复力指数的高值区主要分布在汇水条件好、土壤有机质质量分数较高的河谷川地。从乡镇上来看，尔林兔镇、大堡当镇、中鸡镇和孙家岔镇的恢复力水平较高，马镇、麻家塔镇的水平较低。地统计分析发现，研究区恢复力指数的平均值为0.18，数据分布属于正偏态，偏度为0.59，恢复力指数较低的区域面积大。恢复力指数大于0.5的区域不到10%，说明研究区生态系统恢复力总体不强，生态脆弱。

图3 研究区生态系统恢复力指数与生态退化风险等级Fig.3 Ecosystem resilience index and ecological degradation risk grade in the study area

2.3生态系统退化风险分区

一般情况下，在生态系统的自我调节能力范围内，生态系统多样性越好、生态冗余越多，对外界干扰和胁迫的抵抗能力越强[20]，恢复力指数就越大，生态退化风险就越低。从整体上来看，恢复力指数越大意味着恢复力越强，但是对于不同类型的生态系统，在面临同种强度的干扰时，其退化风险不同。基于此，根据恢复力指数对不同的生态系统类型划分了生态系统的退化风险等级。

根据中国科学院环境保护部2015年发布的《全国生态功能区划》，神木市的主要生态类型包括草原、荒漠、森林、湿地。

根据研究区生态系统退化实际情况[21-23]，确定了神木市不同类型生态系统各评价因子的分级标准，见表4。根据以上确定的神木市生态系统退化风险评价指标分级标准，利用突变级数方法得到的不同退化风险等级的最小恢复力指数，可得到神木市生态系统退化风险评价等级标准，见表5，计算出的评价值范围在0～1，1表示最佳状态，0表示最差状态。由图3(b)可以看出生态退化风险在空间上相间分布，自西向东，由南向北呈现逐渐增强的趋势。其中3级中等退化风险区范围最大，占32.05%，5级和4级严重退化风险区次之，分别占19.94%和16.71%，1，2级轻微风险区最少，分别占6.95%和11.93%。从区域上看，退化风险等级由高到低排序为风沙草滩区>黄土丘陵区>黄河沿岸土山石区。研究区不同类型的生态系统的退化风险差异很大，荒漠生态系统的退化风险等级整体较高，说明植被盖度与土壤水分是引起生态退化的关键因素。湿地生态系统的退化风险整体较低。占研究区总面积75%的草原生态系统和森林生态系统受地形与土壤的影响，退化风险在空间上的分布特征为沟谷低，山脊及梁峁高。

表4 神木市典型生态系统退化风险评价指标分值Table 4 Degradation risk assessment index scores of typical ecosystems in Shenmu City

表5 基于恢复力指数的生态退化风险等级Table 5 Ecological degradation risk grade resilience value

3 讨论

3.1矿区生态系统退化的评估

矿区生态系统退化的评估已经有了一些研究和实践，例如，选择植被盖度作为神东矿区遥感监测和土地沙化评估分析的定量因子[24]，用MODIS影像监测神东矿区土地湿度变化，以此判别采矿扰动地表程度[25]。也有学者通过原位测试和遥感研究了煤矿开采导致的潜水水位变化及其生态影响[26]。可以看出，已有的研究成果主要是针对植被、土壤、地下水等单一生态要素展开的。在评估方法方面，目前运用较多的包括层次分析法[27]、灰色关联分析法[28]、模糊评价法[29]等。地表生态系统实际上是由环境因素和非环境因素共同决定的，对于生态系统恢复力评价这类多指标集成问题，无论采用何种技术，只有减少权重赋值的主观性才能体现评价结论的科学性。

笔者在基于生态系统恢复力，识别生态系统状态的决定因素，统筹考虑生态系统的状态，多变量之间的相关性以及生态系统突变型，确定了生态系统恢复力评价指标。并选择突变级数法进行综合评估，它只需要按照指标间的内在逻辑关系对其重要程度进行排序，避免了直接使用难以确定且主观性较大的权重系数。此外，该方法有效的考虑了生态系统这类内部机理未知的复杂大系统多变量之间的控制关系和综合状态。对比国内外研究，有学者采用RRM模型评价神木县土地整治规划生态风险[30]，也有学者对生态系统的关键要素进行评估，比如通过植被变异系数检测毛乌素沙地东南缘的2000—2018年的植被覆盖变化特征[31]，采用压力-状态-响应模型构建了西部干旱矿区的土壤侵蚀因子等[32]，均反映了东南高、西北低的分布趋势，和本文的总体趋势一致。而笔者采用多因素综合评价，并按照不同的生态系统类型划分退化风险等级，更加精细地反映了不同地形、生态系统类型的退化风险空间格局。

在受到干扰时，生态系统具有自组织能力，通过自组织恢复到原有的状态或进入新的状态。但是如果外来干扰超过一定的限度，任何生态系统都可能受到破坏而进入另一种不同的状态。另外，生态系统原有状态下恢复力的丧失而导致的生态系统状态转化，常转入不理想或不可逆的状态。这种生态系统结构和功能发生的与原有平衡状态或进化方向相反的跃迁，就是生态系统的退化。笔者计算的恢复力指数本质上是多个指标的综合值。从宏观角度来看，恢复力指数越高就代表生境条件越好，受到扰动与胁迫时就越不容易发生退化。但是，不同类型的生态系统发生退化的恢复力阈值不同，一般而言，在生态系统自我调节能力范围内，生态系统各营养级的生物种类越多，营养结构越复杂，对外界干扰和胁迫的抵抗能力就越强;比如在水蚀、风沙等不利自然因素的影响下，林地比草地有更强的抵抗能力。基于此，笔者提出了生态系统退化风险分级标准，按照不同的生态系统类型来确定生态系统发生退化时的恢复力指数阈值，例如相同的恢复力指数，荒漠比林地的退化风险更大。

未来的研究应集中在矿区退化风险的控制与保护方面，基于生态退化风险的评价结果，针对各个区域影响恢复力的关键因素，开发保护修复的具体措施与方法，在空间上组合使用，达到矿区生态修复与保护的最终目的。

3.2 生态退化风险评估的应用

研究区内的神府新民矿区和榆神矿区自下而上共有10余层可采及局部可采煤层，位于含煤岩层上部的2-2煤层为该区域主要开采煤层，该煤层厚度为0.26～12.50 m，平均厚6.50 m(图4(a))，煤层埋深特征的为东浅西深，大部分煤层埋藏深度小于150 m，煤层上层覆盖的基岩厚度小于100 m，属于典型的浅埋藏煤层(图4(b))。

图4 神木市矿区2-2煤层厚度及埋藏深度分布与矿区生态系统保护措施Fig.4 Spatial pattern of coal seam thickness and burial septh in the mining area and Ecosystem protection measures in mining areas of Shenmu City

根据研究区的最主要开采煤层赋存规律，结合生态退化风险评价结果，对矿区生态进行保护与控制。沟谷地带因水源充足生境条件好，退化风险较低，应当持续保护。煤层埋深<100 m或者在100～200 m，生态退化风险达到4,5级的高风险区划分为人工引导修复区。这种区域煤层与地表之间的距离过近，且地下含水层不连续，采煤产生的导水裂隙带发育到含水层底部直达地表，引起地表生态恶化。煤炭资源开发、运输等活动也导致周围土壤受到重金属不同程度的累积性污染。生态系统受到的扰动已经或将要超过其所能承受的阈值，生态系统的恢复力部分或完全丧失，生态系统的退化程度比较严重。可以加强物理手段进行干预，注重土壤和植被这些关键要素的保护，开展土壤改良、植被重建等，通过土地复垦工程技术的手段帮助极度退化的工矿区的生态系统恢复重建。生态退化风险等级为1～3级的低风险区主要分布在新民矿区、榆神矿区东源，这种区域萨拉乌苏组厚度小，分布不连续，地下水贫乏，采后通过简单充填、平整后即可恢复。生态系统的恢复力较强，可以通过生态系统的自我恢复能力实现生态修复，本研究划分的矿区控制与保护分区如图4(c)所示。