基于Landsat时序数据的高潜水位煤矿区植被扰动分析

方良成，陈永春，安士凯，徐燕飞，殷梦杰，李志辉，赵 萍，陈业禹

(1.淮南矿业(集团)有限责任公司，安徽 淮南 232001；2.平安煤炭开采工程技术研究院有限责任公司，安徽 淮南 232001；3.合肥工业大学 资源与环境工程学院 空间信息智能分析与应用研究所，安徽 合肥 230009；4.安徽大学 资源与环境工程学院 安徽省矿山生态修复工程实验室，安徽 合肥 230601)

我国是全球最大的煤炭开采和消费国[1]，煤炭在我国的一次性能源结构中长期保持主体能源地位。然而，煤炭资源的开发利用也不可避免地造成了地表塌陷、大气污染、水污染、土地侵蚀、农作物减产等一系列生态环境问题，阻碍矿区社会经济的可持续发展[2-4]。

植被作为生态环境的重要组成部分，其长势是衡量矿区生态环境的重要指标之一[5-8]。近年来，Landsat系列遥感数据的免费开放为植被长时序变化分析提供了数据源。归一化差值植被指数NDVI(Normalized Difference Vegetation Index)是近红外波段与红光波段反射率的差与和之比，可以灵敏地对植被进行检测，反映地表植被生长状况、植被覆盖度、生物量及环境要素影响[6-9]。众多学者基于NDVI对煤矿区生态环境监测做了大量研究，蔡宗磊等[10]基于GF-1与SPOT6数据，以NDVI等多个植被指数为自变量，对北方露天煤矿区植被覆盖度进行估测，发现基于NDVI的模型精度最高；黄海[11]利用MODIS和Landsat遥感数据反演矿区NDVI等多个植被指数，分析矿区地表植被和土壤湿度的空间时序变化特征和演化机理，有效地揭示西部黄土矿区煤炭开采对地表植被与土壤湿度的影响；黄翌等[12]基于像元二分法，以NDVI值为参数计算半干旱煤矿区植被覆盖度，探讨矿区植被覆盖度的空间相关性和空间异质性的变化；李晶等[13]结合NDVI和综合森林(IFZ)重建阿巴拉契亚地区1987年以来27 a间LUCC动态变化特征，探究采矿扰动随开采规模变化的关系。上述研究对象多为露天煤矿或干旱半干旱矿区，针对我国东部水热条件较好的高潜水位煤矿区研究尚不多见[14]。笔者以安徽淮南顾桥矿采煤沉陷区为研究对象，利用2007—2018年Landsat影像构建NDVI时间序列数据，通过NDVI年际变化分析、热点分析、聚类与异常值分析、剖面线分析，研究煤炭开采沉陷对高潜水位平原矿区植被的扰动效应，以期为矿区生态环境评价提供参考。

1 研究区概况

淮南矿区位于安徽省中北部的淮南市，地理坐标为东经116o20'～117o14'，北纬32o22'～33o22'，处于暖温带和亚热带的过渡地区[15]。矿区地理位置优越，气候适宜，煤炭储量丰富，植被以农作物为主，成片林地与草地较少，是我国华东地区重要的煤粮生产基地。由于地下水位埋藏浅，可采煤层厚度大、层数多、埋深大，重复开采扰动多，淮南矿区煤矿地表沉陷具有规模大、积水深、影响时间长等特点，导致了突出的生态环境问题。顾桥矿位于淮南矿区中西部，正式投产于2007年4月28日，是亚洲地下开采规模最大的矿井之一，煤炭开采形成的采煤沉陷区具有典型性，因此，选取顾桥矿采煤沉陷区为研究对象，探讨煤炭开采对周围植被造成的扰动效应，研究区域地理位置及范围如图1所示。

2 数据获取与处理

2.1 Landsat时间序列NDVI数据的构建

图1 研究区地理位置和范围Fig.1 Geographical location and scope of the research area

本文基于2007—2018年Landsat遥感影像进行长时间序列NDVI指数的构建。为尽可能保证时相一致，选用的12期影像数据采集时间均为每年的4、5月份，影像分辨率为30 m，质量较好，基本无云覆盖，具体情况见表1。

表1 Landsat 系列数据基本信息Table 1 Basic information of Landsat series data

在对遥感影像进行辐射定标、大气校正、交叉定标[16]、裁剪等图像预处理后，根据如下NDVI计算公式，计算得到每一期影像的NDVI值。

式中：NIR、R分别代表近红外、红光波段的像元反射率值，分别对应TM或ETM+影像4、3波段，或OLI影像5、4波段。NDVI值的范围为[–1,1]，反映了植被的生长状态；负值与云、雪、水等相关，0一般代表裸土或岩石，在有植被生长的区域，正值越接近1表示植被覆盖度越高[17]。

2.2 沉陷积水区提取

水体对太阳光具有强吸收性，在近红外、短波红外的波长范围几乎可吸收全部的入射能量。可以利用水体在可见光波段和近红外或短波红外波段的反差来构建水体指数[18]。MNDWI(Modified Normalized Difference Water Index)除能抑制植被的背景噪声外，还可有效抑制土壤及建筑物对水体提取的干扰[19]。MNDWI自提出至今已在国内外得到广泛应用与认可，因此，本次研究采用MNDWI水体指数来进行沉陷积水区的提取，MNDWI的计算公式为：

式中：GREEN、MIR分别代表绿光、短波红外波段的像元反射率值，分别对应TM或ETM+影像2、5波段，或OLI影像3、6波段。

3 NDVI时序变化特征分析

将各年份NDVI影像按如下标准(表2)进行分级，结果如图2所示。

表2 NDVI等级划分标准Table 2 NDVI classification standard

对各年份不同等级NDVI进行统计，如图3所示，可知研究区NDVI值主要分布在0.4～1，所占百分比均在70%以上。随着煤炭的开采，地表发生沉陷积水，NDVI值小于0.1的区域所占百分比逐年增加。研究区内各年份沉陷积水区的扩张(图2)与开采工作面一致，表现为南北向延伸和东向扩展。西北部虽也有少量开采工作面，但由于地表为生产和城镇生活区，开采过程中根据相关规范留设大量保护煤柱，因此，沉陷区未进一步向西扩展[20]。

利用MNDWI提取的沉陷积水区范围制作掩膜，统计掩膜后研究区年际NDVI均值与变异系数，如图4所示。年际NDVI均值在不同阶段有所波动，变异系数总体呈上升趋势。变异系数最小的时间为2007年与2008年，表明煤炭开采初期，研究区内植被的总体生长状况未发生明显改变。随着煤炭的开采，整体的NDVI离散程度逐渐增大。年均值最小的也是2007年和2008年，由文献[5]可知，矿区存在大量耕地，植被NDVI值变化主要受研究区内农作物耕种和收割活动影响，5—6月为冬小麦成熟及收获季节，NDVI呈下降趋势。结合当地历史气温数据及降水数据发现，2007、2008年由于水热条件好，农作物成熟期提前，农作物大多被收割，导致其NDVI与其他年份相比较低。尤其是2007年，年平均气温为17.4℃，高出常年1.8℃，全年除7月份平均气温低于常年0.5℃外，其余11个月份平均气温均比往年偏高且降雨量丰富，雨水偏多60%以上[21]。

图2 2007—2018年NDVI分级Fig.2 NDVI classification from 2007 to 2018

图3 2007—2018年NDVI占比Fig.3 The proportion of NDVI from 2007 to 2018

图4 研究区2007—2018年NDVI均值与变异系数变化Fig.4 The variety of mean value and variation coefficient from 2007 to 2018

4 NDVI空间分异特征

受以煤炭开发为主的多种因素共同作用，矿区植被覆盖度表现出明显的局部依赖性和异质性[22]。热点分析(Getis-Ord Gi*)能够运用Getis-Ord Gi*统计对数据集中的每个要素进行计算，并在空间上获取发生聚类的位置，寻找具有显著统计学意义的热点，即要素在具有高值特征的同时，邻近要素也同样具有高值特征[21]，可以有效揭示植被扰动中的局部效应。聚类与异常值分析(Anselin Loacal Moran I)根据要素位置和要素值对空间的自相关性进行判断，表达的是空间邻近位置上数值间的相互趋同或背离的倾向。因此，空间关联指数Getis-Ord Gi*、Anselin Loacal Moran I可根据空间位置的高低值簇分布情况，有效揭示研究区空间聚簇特征，对矿区植被扰动研究具有重要意义。

4.1 热点分析

热点分析(Getis-Ord Gi*)用来识别具有空间统计性的高值聚集与低值聚集，即热点(Hot Spots)与冷点(Cold Spots)的分布。热点的分布表示NDVI高值集中出现的位置，冷点的分布表示NDVI低值集中出现的位置，从而，可根据NDVI热点和冷点区域的空间演变探测人类活动对矿区植被生长的影响。表达式和原理如下：

式中：Gi*为Getis-Ord Gi*统计；wi,j为以距离规则定义的空间权重，同样空间范围相邻为1，不相邻为0；xj为要素j的属性值；n为要素总和。Gi*统计出的结果即为Z得分。对于具有显著统计学意义的正值Z得分，Z得分表示标准差的倍数，Z得分越高，热点的聚类就越紧密，即NDVI高值聚集的情况越明显。反之，负值的Z得分越高，冷点的聚类就越紧密，即NDVI低值聚集的情况越明显。按照显著性水平分别为“0.01水平(双侧)上显著相关”、“0.05水平(双侧)上显著相关”、“0.1水平(双侧)上显著相关”与“不相关”，将结果显示为冷点区(Cold Spot～99% Confidence)、较冷点区(Cold Spot～95%Confidence)、次冷点区(Cold Spot～90% Confidence)、无显著性(Not Significant)、热点区(Hot Spot～99%Confidence)、较热点区(Hot Spot～95% Confidence)、次热点区(Hot Spot～90% Confidence)。

从图5可知，2007—2018年间热点聚集和冷点聚集区域总体上均在增大。具体来看，2007—2008年，农作物为成熟季节，大部分农作物被收割，热点聚集和冷点聚集的区域较小；2009—2012年热点聚集和冷点聚集的区域明显增大；2013—2014年间热点区和较热点区的数量明显减少，次热点区和冷点区的数量有所增加。2015年以后，热点和冷点聚集区域又开始增加，且逐渐趋于稳定。

从空间分布上看，随着煤炭的开采，积水沉陷面积逐年变大，冷点围绕沉陷区逐年增加。同时，在永幸河与德上高速沿线冷点区增加的现象也较为明显。德上高速与永幸河是矿区煤炭运输的主要道路和居民生活的轴线，植被覆盖度变化十分剧烈，人类活动对地表植被扰动程度显著，大面积的建筑修建使该处的植被覆盖度急剧降低。2007—2008年间热点区分布相对较少，主要分布于研究区北部，2009年后热点分布增多，几乎覆盖整个研究区；2013—2014年间热点分布区减少，大多转化为较热点和次热点；2015年后热点区又开始在北部聚集。

4.2 聚类与异常值分析

聚类与异常值分析(Anselin Loacal Moran I)能有效识别具有统计显著性的空间异常值，发现研究区内NDVI高值与低值聚集的边界以及高低值分布的异常模式，揭示植被的干扰效应。其工作原理如下：

式中：Ii为Local Moran’s I指数；xi为要素i的属性；为对应属性的平均值。

Anselin Loacal Moran I可区分具有显著统计性0.05水平(双侧)上显著相关的“高–高”聚类(High-High Cluster)、“低值被高值包围”类(Low-High Outlier)、无显著性类(Not Significant)、“高值被低值包围”类(High-Low Outlier)及“低–低”聚类(Low-Low Cluster)。

图5 2007—2018年NDVI热点分析Fig.5 Hot spots analysis of NDVI from 2007 to 2018

分析结果如图6所示，2007—2010年间，顾桥矿附近的“高–高”聚类逐渐消失，“低–低”聚类现象有明显增加。2011—2012年“高–高”聚类又再次集中出现；2013—2016年“高–高”聚类效应明显减弱，基本只有“低–低”聚类现象出现且集中出现在沉陷积水区、岗河、永幸河和德上高速附近；2017年以后，“高–高”聚类的区域面积又明显增大。在整个研究时间段内，研究区出现的空间聚集情况主要为“高–高”聚类及“低–低”聚类，无异常出现，表明矿区植被局部空间差异程度小，且聚集现象的增加和减少都是整体性的。

4.3 剖面线分析

顾桥矿正式投产于2007年4月28日，开采方式为地下开采，2008年地表开始出现沉陷。为排除其他人类活动干扰，分析煤矿开采沉陷对周边植被长势的影响，以2008年沉陷积水区几何中心为起点,选择沉陷积水区正北方向剖面线(图1中AB线)对NDVI进行分析。该方向上地表覆盖类型以耕地为主，且变化较小。根据NDVI变化趋势，分别分析煤炭开采不同阶段对沉陷积水区周边植被的影响，将NDVI增减趋势与上一年相比无明显变化的点定义为NDVI稳定值(图7和表3)，图中的虚线表示沉陷起点到影响边界(NDVI值稳定的位置)的距离。

图6 2007—2018年NDVI聚类与异常值分析Fig.6 Clustering and outliner analysis of NDVI from 2007 to 2018

图7 2007—2018年NDVI变化剖面线分析Fig.7 NDVI variation profile analysis from 2007 to 2018

由图7和表3可知，2007—2009年，采煤造成的沉陷积水区显著增大，NDVI值急剧下降的范围由沉陷积水边缘外的360 m扩大到420 m。2010—2013年，沉陷积水区缓慢增大，NDVI受影响的范围由150 m缓慢增大至210 m，且NDVI逐年下降。2013年，因近岸水生植物生长，810 m范围内NDVI呈上升趋势，之后较2012年仍呈下降趋势。2014—2018年，沉陷积水区几乎没有扩张，NDVI影响范围保持为510 m，但在720～1 230 m距离区间内，NDVI仍为逐年下降趋势。其中，2014年，1 650～1 800 m出现了较大的NDVI谷值，经比对高分辨率卫星影像，此处为蔬菜大棚。2017年，650～800 m的NDVI值较2016年也因近岸水生植物的生长出现增长；2018年，800～900 m 范围内NDVI发生急剧下降，经比对该处为鱼塘。

表3 2007—2018年沉陷积水影响分析数据Table 3 The impact of subsidence water from 2007 to 2018

综合上述分析可知，煤炭开采当年对植被扰动影响较小，开采一年后，随着沉陷积水的出现和面积的增长，扰动影响范围不断增加；当沉陷积水区稳定后，其影响范围仍随时间推移持续增加，最后趋于稳定，说明采煤沉陷对植被的扰动存在时序滞后性和时空累积性。

5 结论

a.通过分析不同等级NDVI所占比例及矿区各年份NDVI的变异系数发现，NDVI值大于0.4的植被在研究区占主导地位，所占百分比均在70%以上，但NDVI值的变异系数总体呈上升趋势，表明植被生长状况整体良好，但植被覆盖的离散程度在增大。

b.NDVI热点分析、聚类和异常值分析结果显示，2007—2018年间，研究区内NDVI具有明显的聚集特征，局部空间差异程度小，均为“高–高”聚类或“低–低”聚类，且聚集现象的增加和减少都是整体性的。总体表现为热点在减少，冷点在增加，二者间的转化主要发生在沉陷积水区、德上高速和永幸河附近，主要由采煤沉陷、工矿和居民区建设引起。

c.受煤炭开采影响，沉陷积水区周围一定范围内存在明显的植被扰动效应，随着沉陷积水区范围的增长和时间的推移，影响范围逐渐增大，最后趋于稳定，具有时序滞后性和时空累积性特征。