基于Landsat影像的采煤对地表植被覆盖度时空变化的影响研究

李瑞军， 杨志文， 吴士文， 郝仕龙,3， 邓荣鑫， 尹卿芳

(1.阳泉市大阳泉煤炭有限责任公司，山西 阳泉 045000； 2.华北水利水电大学 测绘与地理信息学院，河南 郑州 450046； 3.南阳师范学院 南水北调中线水源区水安全河南省协同创新中心，河南 南阳 473061； 4.北京京水建设集团有限公司，北京 100193)

植被是陆地生态系统的重要组成部分，是连接土壤圈、大气圈、水圈、生物圈的自然纽带[1]。植被覆盖度作为表征地表生态环境和植被分布特征的主要参数，在水土保持、调节区域小气候、涵养水源等方面发挥着重要作用，具有明显的季候和年际变化特征[2-3]。传统植被覆盖度测量，多数是在研究区选择有限的几个或几十个典型样点进行实地测量，虽能够保证数据观测的精度，但费时费力、成本高、难度大，并且难以满足区域地表植被覆盖度的需求[4]；而遥感技术具有快速实时、监测范围广、人力物力耗费较少、非破坏性等特点，为估算区域尺度的陆表植被覆盖度提供了强有力的手段[5]。植被指数是遥感领域中用来表征植被生长、长势的有效指标，在众多植被指数中，归一化植被指数(Normalized Difference Vegetation Index,NDVI)与植被覆盖度相关性最高[6-8]，已有不少研究将NDVI的像元二分模型应用于植被覆盖动态变化监测分析[3-4,9]。

煤炭作为我国的主要能源，是我国经济和社会发展的重要基础和基本保障[10]。及时、高效地监测矿区生态环境状况，为矿区资源的合理开发提供决策参考，对于实现经济发展和环境效益的有机统一尤为重要。国内外研究发现，煤炭资源的大规模、高强度开采虽然促进了当地经济发展，也不可避免地造成矿区内地表塌陷、水土流失、植被破坏、地下水系紊乱，进而导致生态环境结构、功能被破坏等问题，对地区生态环境和经济社会可持续发展造成了严重影响[11-15]，在生态环境较为脆弱的地区，该问题表现得尤为突出。另有研究结果表明，煤炭资源的合理开发能使矿区的土地沙漠化逆转、植被覆盖率提高，整个生态环境向良性方向发展[16-18]。但以上研究均是针对其它地区，对阳泉煤矿的研究大多集中在煤炭开采、景观生态效应和生态修复等领域[19-21]，在煤矿开采对地表植被覆盖度的时空动态变化及其影响因素研究方面还未见报道。本文以阳泉市大阳泉煤矿不同开采年份工作面为研究对象，基于Landsat TM/OLI数据并结合像元二分模型，对研究区2009年和2019年地表植被覆盖度的分布格局及时空动态变化进行分析；将植被覆盖图与开采年份、各地形因子叠加统计，探究影响地表植被覆盖度变化的主要因子，以期为阳泉市煤炭资源的开发和地表植被生态管理作出科学、准确的评价。

1 材料与方法

1.1 研究区概况

大阳泉煤炭有限责任公司为阳泉市地方国有企业，位于阳泉市西南郊(东经113°30′54″～113°33′6″，北纬37°48′54″～37°50′12″)。1984年开始开采，研究区地理位置及不同开采年份工作面分布如图1所示。阳泉市位于沁水煤田的东北角，其煤炭资源储量约1.7×1010t，埋藏深度150～500 m，易于开发，品质优良；它是我国最大的无烟煤生产基地，也是山西中东部的政治、经济、文化中心。境内以山地为主，地形复杂，大部分地区海拔700 m以上，地势西高东低，属暖温带半湿润大陆性季风气候区，年均气温8～12 ℃，无霜期130～180 d，年降水量450～550 mm。

图1 研究区位置示意图及不同开采年份工作面分布图

1.2 数据来源与处理

本文采用的遥感数据为2009年5月31日的Landsat-5 TM和2019年5月20日的Landsat-8 OLI影像数据，下载自中国科学院遥感与数字地球研究所对地观测数据共享计划网站(http://ids.ceode.ac.cn)，影像质量好，空间分辨率为30 m。运用ENVI 5.1软件对Landsat影像依次进行辐射定标、大气校正和几何校正，将影像的像元亮度值转换为地表真实反射率值。Landsat-5 TM影像为四级产品，基于该影像完成Landsat-8 OLI影像的几何校正。最后，利用研究区的边界矢量数据进行裁切，得到研究区的遥感影像。

DEM数据由阳泉市大阳泉煤矿有限责任公司提供，运用ArcGIS 10.2软件的3D分析模块，生成空间分辨率为30 m的高程图、坡度图及坡向图，通过重新调整高程、坡度和坡向分级并进行重分类，再与植被覆盖度分布图叠加，探讨地形因子及开采年份对研究区植被覆盖度分布的影响；利用ArcGIS 10.2软件，在研究区随机生成359个采样点，提取2009—2019年植被覆盖度的变化、开采时间、高程等信息，利用地理探测器探测各影响因子的重要性。结合研究区实际情况，对各地形因子进行如下等级划分。

1)高程：研究区最大高程1 125 m，最低高程766 m，为方便研究，将研究区高程分为8类，即[766 m,800 m]、(800 m,850 m]、(850 m,900 m]、(900 m,950 m]、(950 m,1 000 m]、(1 000 m,1 050 m]、(1 050 m,1 100 m]、(1 100 m,1 130 m]。

2)坡度：依据《水土保持综合治理规划通则》(GB/T 15772—1995)，结合研究区地形及其对植被长势的影响，将坡度分级改为(0°,5°]、(5°,8°]、(8°,15°]、(15°,25°]、(25°,35°]、(35°,70°)[22]。

3)坡向：把研究区分为8个自然坡向，分别为(0°,22.5°]和(337.5°,360°)(北坡)、(22.5°,67.5°](东北坡)、(67.5°,112.5°](东坡)、(112.5°,157.5°](东南坡)、(157.5°,202.5°](南坡)、(202.5°,247.5°](西南坡)、(247.5°,292.5°](西坡)、(292.5°,337.5°](西北坡)。

1.3 研究方法

1.3.1 植被覆盖度反演

像元二分模型[23]假定单个像元地表S由植被覆盖地表和裸土地表构成。若假定像元地表完全由植被覆盖的纯像元遥感信息为Sveg，完全由裸土覆盖的纯像元遥感信息为Ssoil，则像元植被覆盖度f为：

f=(S-Ssoil)/(Sveg-Ssoil) 。

(1)

要估算像元植被覆盖度，只需要计算土壤和植被的纯像元遥感信息即可。该模型简单实用，通过引入参数减弱了土壤背景、大气以及植被类型等的影响。

归一化植被指数(NDVI)是遥感图像中近红外波段反射值与红光波段反射值之差与二者之和的比值，是反映植被生长状态和植物空间分布密度的重要参数。其比值形式也能够消除大部分与地形、云阴影、太阳角和仪器定标等引起的辐照度变化，以加强对植被的响应。把NDVI代入像元二分模型，其公式为：

fc=(NDVI-NDVIsoil)/ (NDVIveg-NDVIsoil)。

(2)

式中：fc为植被覆盖度；NDVI为归一化植被指数；NDVIsoil为裸土或无植被覆盖区域像元的NDVI值；NDVIveg为完全植被覆盖区域像元的NDVI值，即纯植被像元的NDVI值。

理论上，NDVIsoil的值应接近于0，NDVIveg代表全植被覆盖像元的最大值。但是，由于受光照条件、影像质量、植被类型等众多因素影响，影像不可避免地存在噪声，两者的值均会发生变化。文中，NDVIsoil取影像中NDVI累计概率在5%附近的值；NDVIveg取影像中NDVI累计概率在95%附近的值[3]，以此可反演得到研究区两个时期的植被覆盖度。

1.3.2 地理探测器

地理探测器是王劲峰等提出的一种基于空间分异理论，揭示其背后驱动因子的新的统计学方法[24-25]，主要包括因子探测器、交互作用探测器、风险区探测器等。传统植被变化驱动因子分析方法对模型假设条件和数据要求较多，当因子过多时，计算过程较为繁琐[26]；而地理探测器模型假设制约条件较少，已有学者将其用于荒漠化、NDVI、植被覆盖度等时空演变的驱动因子研究[27-30]。因子探测器的具体公式为：

(3)

式中：q为度量空间分异性的指标，表示某因子解释了q×100%的植被覆盖度；h为影响因子的分层数，h=1、2、…、L；Nh和N分别为影响因子的层h和全区的样本数；δh和δ分别为层h和研究区的植被覆盖度的方差。

q的值域为[0,1]，q值越大，表明该因子对研究区植被覆盖度变化的解释力越强，反之则越弱。

交互探测器用于评估因子共同作用是否会增加或减弱对植被覆盖度的解释力。首先，计算两种因子X1和X2对因变量的q值；其次，计算因子交互的q值，并对q(X1)、q(X2)与q(X1∩X2)进行比较。

2 结果与分析

2.1 植被覆盖度空间分布特征

为更好地研究开采时间对地表植被生长的影响情况，选取远离研究区且没开采过的区域作为对照。因植被覆盖度为[0,1]，为便于进行植被覆盖动态变化分析，结合植被覆盖实际情况，将其分为5级：低植被覆盖度[0.0,0.2]、中低植被覆盖度(0.2,0.4]、中植被覆盖度(0.4,0.6]、中高植被覆盖度(0.6,0.8]、高植被覆盖度(0.8,1.0]，通过密度分割处理，生成植被覆盖度分布图。

研究区2009年和2019年植被覆盖度的空间分布如图2所示，研究区内的空值区域属于非植被覆盖区，已掩膜去除。

图2 研究区2009年和2019年植被覆盖度的空间分布

由图2可知：2009年研究区整体覆盖度较好，多数区域达到了中高植被覆盖度等级，但东北方向区域(2000年以前开采)的植被覆盖度明显低于其它区域，个别区域植被覆盖度甚至低于0.2；到2019年，研究区整体植被覆盖度略有改善，东北方向(2000年以前开采)植被覆盖度也明显提高，植被生长状况明显改善。东北方向地表植被覆盖度较差，可能是由于该区距离城区较近，受人类活动影响较大；近年来，随着生态保护和地表植被恢复工程的开展，到2019年，该区域植被覆盖度有一定提高，但仍略低于其它区域。

表1给出了研究区2009年和2019年植被覆盖度的统计结果。由表1可知：①近10年来，研究区中植被覆盖度及其以上的面积占比均超过了88.86%，植被覆盖总体情况良好；中高植被覆盖度和高植被覆盖度的面积均呈增加趋势。②2009年中植被覆盖度及其以下的面积占比为41.37%，所占比例远高于2019年的20.60%，多出了4.70×105m2；而中高植被覆盖度的面积在2009年占55.05%，到2019年，占比增加到了64.76%；高植被覆盖度面积增加比例最为显著，2019年的面积占比比2009年增加了308.89%，达到了14.64%。

表1 2009年和2019年各等级植被覆盖度的面积及占比

2.2 地形因子对植被覆盖度的影响

山区地形特征是一个多维变量，地表植被除了受降水、气温等气象因子影响外，在小范围区域，地形因子对植被覆盖度空间分布的影响尤为重要。地形因子主要影响不同区域的光、热、水、土等资源，进而影响植物生长。为探究矿区植被覆盖度随地形因子变化的空间分布规律，本文提取了高程、坡度、坡向等地形因子，并划分等级，分析地形因子对植被覆盖度空间分布的影响。

2.2.1 植被覆盖度随高程的变化特征

在一定高程范围内，山地存在气温垂直递减和相对湿度随高程升高而增加的规律；同时，随着高程的变化，不同高度上接受到的太阳辐射、热量，气温以及土壤性质等也会发生变化，进而影响地表植被覆盖度的分布[31]。表2为研究区高程分级统计表，在[766 m,1 130 m]范围内，不同高程所占面积比例随高程的增加而呈先上升后降低的趋势；研究区高程主要分布在(950 m,1 100 m]，占总面积的67.99%。

表2 研究区高程分级统计表

图3给出了植被覆盖度与高程的叠加分析结果。由图3可以看出：2019年所有高程范围内的植被覆盖度均高于2009年的；随着高程的增加，植被覆盖度逐渐升高，达到900 m时逐渐平稳，1 050 m之后又逐渐降低，(900 m,1 050 m]高程范围内的区域是植被生长状况良好的区域。这是因为低高程区，虽然水分、热量和气温都适合植被的生长需要，但该区域受人类活动干扰大；随着高程升高，人类活动影响减少，热量和水分仍适合植被生长，植被覆盖度开始上升，在(900 m,1 050 m]高程范围内，植被生长最好；当高程大于1 050 m后，随着高程升高，气温明显下降，热量成为限制植被生长的主要因素，则植被覆盖度下降[22]。

图3 2009年和2019年植被覆盖度随高程的变化

2.2.2 植被覆盖度随坡度的变化特征

坡度表示局部地表的倾斜程度，直接影响物质能量的转换方式与程度，与热量、水分条件、土壤厚度等因素关系密切，可对土壤水分和养分的分布进行调控，改变土壤的基本属性，在很大程度上影响地表植被的分布[32]。研究区坡度分级统计情况见表3。由表3可以看出：研究区地表坡度主要在(8°,35°]范围内，该坡度范围的面积占研究区总面积的比例为83.14%；而坡度为(35°,70°)的陡坡区域面积只占研究区总面积的7.12%，说明整个研究区的地势较为平坦，起伏不大。

表3 研究区坡度分级统计表

图4给出了研究区植被覆盖度与坡度的叠加分析结果。由图4可以看出：2019年研究区的植被生长状况在所有坡度区均优于2009年的；坡度较低区域的植被覆盖度变化趋势较为平缓，坡度大于8°后逐渐升高，在坡度为(25°,35°]时植被覆盖度达到最大值，2019年的植被覆盖度最高达到了0.73；坡度大于35°之后，植被覆盖度开始降低。

图4 2009年和2019年植被覆盖度随坡度的变化

总体来讲，除坡度为(35°,70°)的陡坡区之外，坡度越大，植被覆盖度越高。这是由于地势平坦区域是人类生产活动影响较大的区域，植被覆盖度较小；随着坡度增加，坡度较大区域受人类活动影响较小，该区域土壤持水能力和营养成分都适合植被生长，尤其有利于草地、灌木等覆盖度较高的植被生长；而坡度为(35°,70°)的地区，地势过于陡峭，水土流失严重，土壤养分也不易保存，植物生长受限。

2.2.3 植被覆盖度随坡向的变化特征

坡向表示每一个栅格高程值变化量的改变方向，主要影响太阳辐射和土壤含水量，从而改变植被的生长和分布特征[22]。表4为研究区坡向分级统计表。由表4可以看出：除南坡的面积占比仅6.32%以外，其余坡向的面积占比为10.02%～16.71%，各个坡向上的分布面积比较均匀。图5给出了研究区植被覆盖度与坡向的叠加分析结果。

表4 研究区坡向分级统计表

图5 2009年和2019年植被覆盖度随坡向变化

由图5可以看出：2019年植被生长状况在各个坡向上均优于2009年的；2009年和2019年各坡向植被覆盖度最大值与最小值差值均为0.06，各坡向间差异不大，总体表现出东坡、南坡和东南坡植被生长状况最差，西北坡生长状况最好的特点。这是由于研究区属半湿润大陆性季风气候，年降雨量较少，仅450～550 mm，东坡、南坡和东南坡方向虽光热条件充足，但蒸发量过大，易使土壤水分含量不足；西北坡接受太阳辐射较少，土壤水分易于保持，光热条件适宜，更有利于植被生长。因此，西北坡植被生长状况优于东坡、南坡和东南坡的[29]。

2.3 植被覆盖度的时空演变及其驱动因素

2.3.1 植被覆盖度的时空变化趋势

为揭示近10年来植被生长状态的空间变化特征，利用ArcGIS软件对2009年和2019年的植被覆盖度图进行栅格叠加运算，得到研究区近10年来植被覆盖度的时空变化分布图，如图6所示。从图中可以看出，2009—2019年研究区内植被覆盖度的变化整体上呈增加趋势。表5为植被覆盖度变化等级的面积统计，近10年来，植被覆盖度增加的区域面积占研究区总面积的81.70%，植被覆盖度降低的区域面积仅占总面积的18.30%，植被覆盖度增加的区域面积远高于降低区域的面积，且以[0.0,0.2]的植被覆盖度变化等级所占面积比例最大，达到了71.64%。结果表明，近10年来，研究区植被覆盖度整体呈上升趋势，植被生长状况明显改善。

图6 2009—2019年研究区植被覆盖度的时空变化

表5 2009—2019年植被覆盖度变化的面积统计

为探究开采时间对近10年来植被覆盖度变化的影响，从图6中提取了2009—2019年不同开采年份工作面的植被覆盖度变化量；同时，为降低不同时期影像的误差影响，提取的不同开采年份工作面的植被覆盖度变化量需同时减去对照区的变化量。图7为不同开采年份工作面2009—2019年植被覆盖度变化量。由图7可知：植被覆盖度的变化以2007年为界限分为两段，开采年份早于2007年的煤矿，除1984年、1990年、1998—2000年开采工作面的植被覆盖度变化量为负值外，其余开采年份植被覆盖度变化量均大于0，而2007年后开采的所有煤矿开采工作面的植被覆盖度变化量均小于0；与其它年份相比，1984年开采工作面的植被覆盖度变化明显异常。这可能是由于其工作面距城区较近，受人类活动影响较大，且工作面面积较小，也增大了分析误差。综上可知，煤矿开采后，虽然短期内会对地表植被生长状况造成影响，但自然生态系统的自我调节能力辅以人工保护措施，会使生态系统逐渐向良性循环方向发展，植被生长状况也会不断得到改善[16-18]。

图7 2009—2019年不同开采年份工作面植被覆盖度变化量

2.3.2 植被覆盖度影响因子的地理探测

由影响因子的探测结果发现，各因子q值排序为：开采时间>高程>坡向>坡度。开采时间是影响植被覆盖度时空变化的最主要因素，其q值为0.14；高程为次要影响因素，其q值为0.12；坡度和坡向虽然是影响植物生长的重要因素，但这2个因子的q值均低于0.05，且坡度的q值最低，说明2009—2019年植被覆盖度变化主要受开采时间和高程的影响，坡度和坡向的直接影响较小，这与前人的研究结论一致[30,33]。

表6为不同影响因子的交互作用的探测结果，从表6中可以看出，不同影响因子的两两交互作用均大大增加了单因子对植被覆盖度变化的解释力。交互作用q值排序为：开采时间∩高程>高程∩坡向>开采时间∩坡向>高程∩坡度>开采时间∩坡度>坡度∩坡向。开采时间、高程和坡向的两两交互作用较强，q值均在0.33以上，坡度与坡向的交互作用q值最低，仅0.18。

表6 不同影响因子的交互作用的探测结果

3 结语

1)2009年和2019年，研究区植被生长总体情况良好，中植被覆盖度及其以上面积占比均超过了88.86%，中高植被覆盖度和高植被覆盖度面积均呈增加趋势；空间分布上，东北方向2000年前开采工作面的地表植被覆盖度明显低于其它区域。

2)从地形因子对植被覆盖度的影响分析发现，植被覆盖度总体随高程和坡度增加呈现先增加后减小的趋势，在高程(900 m,1 050 m]和坡度(25°, 35°]的地带达到最大值；在坡向上呈现出西北坡最高，而东坡、南坡和东南坡相对较低的趋势，且高程对植被覆盖度分布的影响大于坡度和坡向。

3)对比2009—2019年植被覆盖度的时空变化发现，2007年以前开采的大多区域植被覆盖度有明显改善，而2007年后开采的区域植被覆盖状况呈退化趋势，表明虽然短期内地下采煤会对植被覆盖度造成影响，但生态系统的自我调节能力再辅以人工保护措施，会使生态系统逐渐向良性循环方向发展，植被生长状况逐步改善。

4)因子探测结果同样表明，开采时间是影响研究区植被覆盖度时空变化的最主要因素，高程次之，坡度的q值最低；交互作用探测中，以高程和开采时间的共同作用最强，说明人为因素与自然因素的协同能够增强对植被覆盖度变化的影响力。