采矿对植被变化的影响提取与生态累积效应量化分析

2022-08-01李全生许亚玲张成业郭俊廷佘长超宋子恒

煤炭学报 2022年6期

李全生，许亚玲，李军，张成业，郭俊廷，佘长超，宋子恒，王菲

(1.煤炭开采水资源保护与利用国家重点实验室，北京 102209；2.中国矿业大学(北京)地球科学与测绘工程学院，北京 100083；3.国家能源投资集团有限责任公司，北京 100011；4.北京低碳清洁能源研究院，北京 102211；5.神华北电胜利能源有限公司，内蒙古锡林浩特 026015)

矿区生态系统是一个受气候环境、地形地貌、农牧活动、城镇扩张以及采矿等工业活动综合影响的复杂耦合系统，其生态系统具有时空累积性和多因素综合驱动的典型特征。传统观念认为，矿区生态的退化主要是采矿活动引起，然而近几十年的实践表明，开采活动确实会对矿区生态产生影响，但影响范围和程度有限。为科学分析采矿对生态系统的影响，量化采矿对生态系统的时间累积效应和空间累积效应，并对矿区生态系统进行有效修复和防治，迫切需要一种在多因素综合作用下分离出采矿活动造成的生态影响的方法。

在众多生态要素中，植被是最重要的要素之一，综合表征了一个地区生态环境状况，因此煤炭开采活动对植被的影响成为学者们关注的重点。国内外学者对此展开了较多的研究，主要分为两大类方法。一类是基于地面调查采样数据的研究，例如许传阳等、王锐等通过实地调查矿区地表土壤物化性质、植被光合速率、植被产量等研究了沉陷裂缝对作物的影响。杨金芳等、马守臣等、顾鑫等通过实地采样研究了矿区特有的矿井水、煤炭腐殖酸等对植被的影响。另一类方法是利用遥感技术监测随采矿过程推进时的植被状态变化，例如，王双明等基于Landsat影像通过拟合TVDI变化趋势线反解煤炭开采活动对矿区生态的影响边界。YANG等利用多时相Landsat影像定义开采扰动模板，并对煤炭开采的草地扰动进行识别。此外，还有许多学者将实地采样与遥感影像结合起来，评价采矿过程在不同开采尺度和阶段对不同植被类型的影响。

然而，现有研究均假定直接监测到的植被状态变化代表了采矿的影响，但事实上，植被状态变化是自然条件、采矿及其他人类活动等综合作用的结果，因此当前研究鲜有将采矿活动和其他因素产生的生态扰动分离开，导致无法揭示采矿活动单独对植被的影响规律和特性。

针对上述问题，笔者基于长时序的植被遥感监测数据、气候气象数据、人类活动数据等多源信息，提出了一种提取和量化煤炭开采对植被影响的方法，即提取采矿干扰两阶段法TS-EMD(Two Stage-Extracting Mining Disturbance)。首先，利用采矿前的长时序多源数据构建植被变化的驱动模型；然后，利用上述模型预测无采矿活动条件下的植被演变过程，进而与遥感监测的采矿活动背景下实际的植被演变进行对比，提取出采矿活动对植被变化的影响。以内蒙古自治区胜利矿区为研究区，利用TS-EMD方法在胜利矿区提取和量化了煤炭开采活动对植被的影响，剖析了采矿影响的时空分布差异和生态累积效应，并揭示了不同采矿阶段的演变规律，为矿区的生态环境保护和生态修复决策提供理论数据支持。

1 研究区与数据

1.1 研究区概况

内蒙古自治区锡林郭勒盟锡林浩特市，位于43.03°N—44.87°N，115.3°E—117.1°E，如图1所示。

图1 研究区地理位置

锡林浩特市的地势南高北低，海拔970～1 212 m，属于中温带半干旱大陆性气候。辖区内草原、矿产等自然资源富集，是国家重要的绿色畜产品生产加工输出基地与煤电基地。

笔者以2020年的胜利一号和西二露天矿区边界以外5 km为研究区，探索区域内采矿活动对周边植被的影响。

根据煤炭开采数据，将1995—2003年作为无采矿活动阶段，将2004—2020年作为有采矿活动阶段。

1.2 数据来源与处理

本研究使用的主要数据有遥感数据、气候数据、放牧数据、煤炭开采数据和基础地理信息数据。其详细信息见表1。

表1 研究数据集

Landsat数据主要用于土地利用分类和植被覆盖度反演。根据国家现行土地利用分类标准(GB/T 21010—2017)及胜利矿区的区域特点，将土地利用类型分为工矿用地、城镇建设用地、草地、林地、耕地、裸地、水体7类。总体实地采样数量为260，采用神经网络方法对所选取的训练样本进行训练，然后对胜利矿区的遥感影像进行信息提取，并对分类结果目视判读，研究区工矿用地和城镇建设用地较为集中，城镇建设用地集中在研究区东南角，通过遥感影像对其进行人工修正，最后得到胜利矿区26 a的土地利用分类，利用混淆矩阵对分类结果进行精度评价，分类总体精度为93.47%，部分年份分类结果如图2所示。

图2 土地利用分类结果

反演植被覆盖度需要对Landsat数据做时间和云量筛选、去云处理、指数计算、镶嵌、裁剪等处理。首先选择研究区云量小于5%的影像，时间周期按照每年7月1日至9月1日进行筛选，利用GEE(Google Earth Engine)遥感云计算平台对Landsat数据进行去云处理，在线计算研究区1995—2020年逐年的归一化差值植被指数(Normalized Difference Vegetation Index，NDVI)，利用最大值合成算法获取NDVI最大值空间分布结果。植被覆盖度(Fractional Vegetation Cover，FVC)指一个区域内植被在地面的垂直投影面积占区域总面积的比例，笔者利用像元二分模型计算FVC，如式(1)所示。笔者得到了1995—2020年共26 a逐年的FVC遥感监测数据，部分年份结果如图3所示。

图3 FVC遥感监测结果

(1)

其中，FVC为像元的植被覆盖度；NDVI为该像元的归一化植被指数；NDVI，NDVI分别为影像中纯净植被像元和纯净裸土像元的NDVI。本研究以植被完全覆盖区域生长旺季观测值的95%代表NDVI，以完全处于裸土状态、无植被覆盖区域全年NDVI观测值的95%代表NDVI。根据气象站点数据，研究时段无极端天气变化情况。地形数据需要按研究区边界镶嵌裁剪。放牧数据为锡林浩特市历年以嘎查(村)和分场为统计单元的牛、羊和马等牲畜头数。根据锡林浩特市林草局草畜平衡政策，将牛和马的数量换算成羊的数量，并以单位面积内羊的数量将放牧数据空间化，部分年份的结果如图4所示。

图4 放牧数据空间化结果

1.3 研究区样点选取

根据土地利用分类数据，提取研究区1995—2020年草地、林地部分(耕地由于人工灌溉不考虑在内)，城镇建设用地等土地利用类型不参与样点选取，以避免干扰。在研究区内每隔400 m选取一个样点，样点分布如图5所示。

图5 研究区样点

2 研究方法

胜利矿区植被覆盖度受到自然条件、人类活动等多种驱动因素的耦合影响，在矿区扩张、采煤强度提高的前提下，会导致植被和社会经济结构发生改变，进而导致牲畜数量变化而对植被影响产生变化，因此为避免其耦合影响，以单位面积羊的数量作为放牧强度因子参与时空地理加权回归，可以剔除由煤炭开采—植被变化—社会经济结构变化—牲畜数量变化—植被变化链式反应的影响。在数据处理阶段，通过逐步回归方法对气温、降水、放牧、人口、GDP、煤炭开采等众多因子进行筛选，得到对FVC影响最大的4个因子：气温、降水、放牧、煤炭开采，因此重点对这4个因子进行分析。这些因子在不同地理位置、不同时间的影响程度不同，存在地理空间异质性与时间异质性。在复杂的矿区生态场景中，植被受多种因素的耦合影响，为了分析煤炭开采对植被变化的单独影响，本文提出一种提取和量化煤炭开采对植被单独影响的方法，即提取采矿干扰两阶段法TS-EMD。该方法具体实现步骤如图6所示。

图6 提取采矿干扰两阶段法

第1步：在1995—2003年无煤炭开采活动时期，利用长时序逐年的FVC、降水、气温、放牧等数据，选取地理时空加权回归模型(Geographically and Temporally Weighted Regression，GTWR)，建立降水、气温、放牧3种驱动因素与FVC的定量关联关系。研究区小范围内气温、降水在空间尺度上可近似看作是均一的，但在时间尺度上是变化的。GTWR在权重矩阵中集成了时间和空间信息，以捕获空间和时间的异质性，即在某一样本点位置处搜索其时空范围内的邻近样本点共同参与回归以构建该点的模型，这样可以充分反映气温、降水的变化(尤其是时间上的变化)导致的植被覆盖度的变化。GTWR建模方法如式(2)所示。

FVC=(,,)+(,,)+

(,,)+(=1,2,…,)

(2)

其中，FVC为第个样本点处的植被覆盖度;为降水量，mm;为气温，℃；为放牧量，头/km；(,,)为时刻(,)位置处各自变量的系数；为样本总数；为误差项。

第2步：利用训练得到式(2)所示的定量关系，输入2004—2020年的降水、气温、放牧数据，可以预测出没有煤矿开采情况下研究区逐年的FVC。事实上，2004—2020年研究区存在不同强度的煤炭开采活动，所以预测的FVC与遥感监测的实际FVC差值可认为是剔除气温、降水、放牧等影响后，由煤炭开采导致的FVC变化量，将其定义为采矿植被覆盖度扰动量FVC-MD(Fractional Vegetation Cover-Mining Disturbance)。利用FVC-MD可以剖析煤炭开采对植被影响的时空分布差异和演变量化特征，进而发现煤炭开采对植被影响的客观规律。

3 结果与讨论

3.1 植被覆盖度的时空变化特征

图7为1995—2020年期间研究区的平均FVC，其取值在0.38～1.00波动。对这26 a的FVC值进行线性拟合，斜率为-0.011 95，斜率为负说明这26 a间的FVC总体上呈现下降趋势，植被发生退化。计算每年的FVC标准差并进行线性拟合，斜率为0.001 61，斜率为正说明矿区的FVC值标准差总体上呈现上升趋势，区域内植被覆盖度的空间异质性增强。

图7 研究区1995—2020年植被覆盖度及标准差变化趋势

为进一步分析胜利矿区FVC的时空变化特征，采用绿度变化率(RCG)量化FVC在1995—2020年期间的年际变化趋势。针对植被覆盖度FVC，计算RCG的方法如式(3)所示，RCG计算结果如图8(a)所示。

图8 1995—2020年植被覆盖度变化趋势

(3)

其中，为监测年份，取值为1～26；为总年数，设为26 a；FVC为第年的FVC值。如果RCG为正值，说明植被覆盖度存在上升趋势，相应区域的植被呈现改善状态；如果RCG为负值，说明植被覆盖度存在下降趋势，相应区域的植被呈现退化状态。

根据RCG的计算结果，将胜利矿区的植被变化趋势划分为5种类型：极显著退化(RCG<0,<0.01)、显著退化(RCG<0,0.01<≤0.05)、变化不显著(>0.05)、显著改善(RCG>0,0.01<≤0.05)、极显著改善(RCG>0,≤0.01)，如图8(b)所示。

在空间分布上，露天煤矿的采场和锡林浩特市的城区扩张，对土地有明显的直接压占，造成植被破坏，表现为植被的“极显著退化”；“显著退化”区域主要分布在矿区和城市周边一定范围的平原地带；“显著改善”和“极显著改善”的地区很少，大都零散地分布在城镇内部；其他大片区域变化为“变化不显著”。

露天采矿具有特殊的地貌形态，除了露天采场外，还具有大型外排土场，采场对植被覆盖的影响是通过压占土地，直接改变土地利用类型，导致采场处的植被退化为零，影响机理简单清晰。而外排土场处的植被受到煤炭开采、人工修复等复杂影响，根据胜利矿区生态修复6步法“一排、二覆、三沙障、四种、五灌、六养护”可知，外排土场影响机理较为复杂，因此将其纳入研究，提取了胜利一号矿和西二矿排土场的FVC年际变化曲线，如图9所示。由于2004年以前不存在排土场，因此1995—2003年统计的是2004年排土场范围内的FVC均值。

图9 排土场FVC年际变化曲线

整体上看，胜利一号矿和西二矿排土场FVC变化趋势一致，均表现为先降低后波动上升。2004—2006年排土场区域由于土地压占导致植被退化，FVC呈下降趋势；2007年后FVC逐渐上升，至2013年完成第1轮生态修复，胜利一号矿排土场FVC稳定在0.4左右，西二矿排土场FVC在0.2～0.4波动，FVC均大于开采前期，说明排土场区域生态修复效果明显。

3.2 无采矿时的植被变化驱动过程

将FVC与降水、气温分别进行了相关性分析和显著性检验，见表2。在降水方面，FVC与当年7月份的降水量相关系数最高，达到0.863。在气温方面，FVC与7—9月积温的相关性最高，达到-0.942。因此，笔者选择7月份的降水量与7—9月积温作为驱动因子。式(2)中的放牧因子采用换算后的羊数量。利用GTWR模型方法建立研究区1995—2003年FVC与降水、气温、放牧等驱动因子的定量关系，结果如图10所示。

表2 FVC与降水、气温的相关性

由于时空异质性的存在，每一个位置在不同年份都有各自对应的驱动模型，利用交叉验证法计算得到总体的建模精度为84.7%。

本模型是用1995—2003年的数据构建，为了验证模型的预测精度，利用训练得到的GTWR模型预测同一区域1991年和1992年的FVC，并与实际FVC做比较，计算得到总体预测精度为78.3%。

图10(a)，(b)分别为1996，1997，1999年驱动模型中降水和气温的系数分布。由图10(a)可知，研究区FVC与降水均呈正相关(系数均为正值)，中部地区降水的系数较高，植被受降水影响较大；降水的系数从1996—1999年逐渐增大，说明降水影响程度逐渐加大。由图10(b)可知，研究区FVC与气温均呈负相关(系数均为负值)，北部地区气温的系数绝对值较高，说明北部地区FVC受气温影响较大，而东、西部受气温影响较小；且与降水影响相似，从1996—1999年气温影响程度逐渐加大。

图10(c)为部分样点的植被FVC变化的驱动方程，同一个样点在不同时期具有不同的回归方程，该方程可以用于预测在无采矿背景下的FVC值。例如：1997年1号样点回归方程为：FVC=1.983-0.000 624+0.001 84-0.024 0，表示1997年在该位置处，降水量每升高10 mm，FVC会上升0.018 4，即植被覆盖度增加1.84%；气温每上升5 ℃，FVC会下降0.12，即植被覆盖度减少12%；每平方千米增加10只羊，FVC会减少0.006 24，即植被覆盖度减少0.624%。2号样点位于锡林浩特城区边界，没有放牧活动，FVC仅受气温和降水因子影响。上述结果表明，GTWR方法能够建立植被FVC与降水量、气温、放牧量之间的定量关系，从而为采矿活动背景下剔除降水、气温和放牧活动对FVC的影响提供定量模型基础。

图10 GTWR建模获取的FVC与降水量、气温、放牧量之间的定量关系

3.3 采矿植被覆盖度扰动量(FVC-MD)提取及其时空分布规律

3.3.1 FVC-MD提取结果

使用TS-EMD第2步方法，输入2004—2009年的实际降水、气温、放牧数据，预测在没有煤炭开采情况下研究区逐年的FVC，并与遥感监测到的有煤炭开采情况下的实际FVC比较，即可提取得到植被变化采矿影响——采矿植被覆盖度扰动量FVC-MD。

以2015年为例展示FVC-MD提取过程：图11(a)是由模型预测得到的2015年研究区每个样点的FVC，图11(b)为遥感监测到的实际FVC，2者做差运算即可得到2015年的FVC-MD，如图11(c)所示。

图11 FVC-MD提取结果示意(以2015年为例)

3.3.2 FVC-MD时空变化特征

图12展示了2004—2020年的FVC-MD分布。从时间上看，2004—2009年，西二、三矿还没有开采，高扰动(FVC-MD的高值)主要集中在胜利一号矿附近；2009—2015年，西二、三矿周边FVC-MD升高，说明该区域植被受煤炭开采影响逐渐增大；2016和2017年，研究区中部部分地区FVC-MD降低，说明矿区生态修复显见成效。2019年后，由于采矿活动的继续加大以及采矿对植被影响的时间累积效应，FVC-MD逐渐增大，到2020年大面积区域的FVC-MD达到0.4。从空间上看，FVC-MD高值多集中在采场周边有限范围内，低值主要分布在离矿较远的区域。西二、西三矿南侧FVC-MD值总体上高于北侧，胜利一号矿东侧FVC-MD值总体上高于西侧，说明煤炭开采对植被的扰动在不同方向上可能存在异质性。图12中的绿色区域表示FVC的实际值高于无采矿条件下的预测值，主要位于矿区西北方向的山地和东北方向锡林河的河谷湿地，其主要原因如下：本文的放牧数据是从锡林浩特市纸质统计年鉴资料中获取，空间化后的最小统计单元为嘎查(村级)和分场。在山地和河谷湿地区域，放牧强度采用了其所在嘎查和分场的均值来量化表达，但是本研究经过实地核查发现这些区域的放牧活动很少；换言之，在山地和河谷湿地区域，放牧强度量化值偏高，会导致FVC的预测值偏低，此外，其他小型人类活动也会对预测精度造成影响。因此，从数据本身以及建模方法考虑，尚有进一步改善之处：一是采用了放牧数据参与建模，放牧数据以嘎查(村级)为统计单位，空间分辨率低，会对建模精度产生一定影响，但就建模精度84.7%来看仍能满足建模要求，在离矿区较远的局部区域(山地和河谷)，通过改善放牧数据的空间分辨率将有助于进一步提升植被覆盖度预测精度。二是两阶段法TS-EMD模型采用线性回归，假定气温和降水与植被覆盖度之间呈线性相关，且作用效果不随时间发生变化，但各影响因子相互耦合，影响模式复杂，在后续研究中还可采用非线性建模进一步研究。此外，用历史数据建模预测未来从而分离出采矿因素的影响需要一些理想化的假设，假定在研究时段内，除了参与建模的影响因子外，其他条件均不变，由于现实与理想假设的偏差可能会引入误差。本文重点关注采矿对植被的单独影响，各因素之间的差异性研究也十分必要，需要在后续的研究中更加深入的探讨。

图12 FVC-MD时空分布

图13为研究区FVC-MD均值、开采原煤量随时间的变化关系，根据开采原煤量的年际变化将整个采矿周期分为3个阶段，2004—2011年为第1阶段(开采量迅速攀升)，FVC-MD均值随着开采量的增大上升，于2011年达到0.27。2012—2016年为第2阶段(开采量呈下降趋势)，在此期间，FVC-MD均值在0.17～0.25波动，原因可能在于：当前开采量下降对扰动量的负影响与历史采矿生态累积效应的正影响共同作用。2017—2020年为第3阶段(开采量回升)，随着煤炭行业的回温，开采量继续增大，FCV-MD也继续增加，且速率增快。

图13 FVC-MD与原煤量年际变化曲线

3.3.3 FVC-MD距离衰减特性

以2020年胜利一号矿采场边界为起点，沿西北方向划定一条样区，如图14所示，统计区域内不同时间阶段(2004—2011年、2012—2016年、2017—2020年)的FVC-MD均值，提取统计结果。图15为胜利一号露天矿采场西北方向各个开采阶段里FVC-MD与离矿距离关系，可以看出FVC-MD与离矿的距离呈现明显的负相关，最优拟合曲线为对数函数曲线，3个阶段拟合均大于0.8。由图15可知，随着离矿距离的增长，采矿对植被覆盖度的扰动量降低，且降低速率变缓。从第1阶段到第3阶段，依次变大(0.83<0.86<0.94)，说明随着采煤时间的推移，FVC-MD与离矿距离之间的对数函数关系越来越显著。

图14 样区示意

图15 FVC-MD与离矿距离的关系

第1阶段FVC-MD最大值为0.36，然后随距离逐渐降低，在距离矿边界4.5 km处趋于稳定，稳定值在0.08附近；第2阶段FVC-MD最大值为0.47，在距离矿边界5.2 km处趋于稳定，稳定值为0.1；第3阶段，FVC-MD最大值为0.57，在距离矿边界6.7 km处还未达到稳定状态。上述结果说明，随着时间的推移，越靠后的阶段FVC-MD越高、衰减距离越远，采矿对植被影响表现出一定的时间累积效应。

3.3.4 FVC-MD的方向异质性

为探究采矿影响在不同方向上的差异，以每个阶段的最大采场边界为中心、1.3 km为间隔向外依次建立缓冲区，按照不同矿、不同方向将缓冲区划分为16个子样区(图16)，方向为北方向(E)、西北方向(NW)、西方向(W)、西南方向(SW)、南方向(S)、东南方向(SE)、东方向(E)、东北方向(NE)。其中胜利一号矿的东南方向和东方向第2，3阶段由于被城镇扩张而出现土地压占，因此不做分析。统计胜利一号矿和西二矿不同方向(不同子样区)、不同阶段的FVC-MD均值，并以雷达图的形式展示，如图17所示。

图16 缓冲区与子样区示意

由图17可以看出，在距离矿边界同一距离下，采矿对植被的影响程度在不同方向上表现为不同的特征。胜利一号矿表现为：以南北方向为分界，西方向FVC-MD小于东方向；西二矿表现为：以东北—西南方向分界，西北方向FVC-MD明显小于东南方向，其中南方向植被受采矿影响最严重；随着时间的推移，2个矿区越靠后的阶段FVC-MD越大。

图17 不同方向FVC-MD的雷达

同时，FVC-MD在任意方向都表现出随离矿距离衰减的特性。其中在第2，3阶段，胜利一号矿的西南方向在距矿边界2.6～3.9 km内FVC-MD大于1.3～2.6 km的区域，这是由于2.6～3.9 km的区域靠近西二矿，同时受到胜利矿和西二矿的叠加影响，FVC-MD较高。同理，西二矿的东北方向也是如此。

为了探究采矿在不同方向上对植被扰动差异产生的原因，以40 m为间隔将研究区地形划分为如图18(a)所示的7个高程带，统计每个高程带内的FVC-MD均值，如图18(b)所示。FVC-MD与海拔高程整体上呈明显的负相关关系，高程越高，FVC-MD越小。7个高程带中，950～990 m平均FVC-MD最高，值为0.165，高程1 190～1 320 m区域平均FVC-MD最低，值为0.089，其间1 150～1 190 m高程带FVC-MD有所上升，值为0.115。结合不同方向FVC-MD特征以及FVC-MD随高程变化特征可知，海拔高程可能是导致采矿在不同方向上对植被扰动差异产生的主要原因，地势低的区域植被受煤炭开采影响较严重。此外，其他土壤湿度和其他小型人类活动也会导致采矿在不同方向上对植被扰动产生差异。

图18 研究区FVC-MD与海拔高程的关系

3.3.5 FVC-MD的空间累积效应

为了量化分析多个矿对植被影响在空间上的累积效应，按图19所示的样区测线计算每个阶段内的FVC-MD均值，结果如图20所示，沿西二矿-胜利一号矿的样区测线，FVC-MD空间变化特性表现为：靠

图19 样区示意

图20 西二矿-胜利一号矿样线FVC-MD变化

近矿边界的位置FVC-MD高，随距离递减，在两矿中间位置FVC-MD增大，存在高扰动区。

在第1阶段，叠加影响区的FVC-MD在0～0.2，在距西二矿2.29～4.25 km，距胜利一号矿采场边界3.27～5.22 km存在植被高扰动区；第2阶段FVC-MD在0.04～0.40，在距西二矿1.63～3.59 km，距胜利一号矿采场边界2.61～4.57 km存在植被高扰动区；第3阶段叠加影响区FVC-MD在0.22～0.49，两矿之间植被退化量均较高；由此可以看出，随着时间的推移，FVC-MD在增大，高扰动区的范围也在扩大。由于单矿对植被覆盖度的扰动量与离矿距离呈对数函数关系，拟合大于0.8。而西二矿和胜利一号矿2个采场之间的FVC-MD变化曲线在数学上恰巧吻合2个对数函数的叠加，因此这种现象很可能是由于采矿活动在空间上的累积效应引起的，但要想完全确认这种因果关系，需在后续研究中做更加深入的空白对比研究。高扰动区虽然离各个矿都有一定距离，但由于2个矿在空间上对该区域的叠加影响，植被受采矿影响较为严重，需要重点关注。