归一化植被指数（NDVI）在草原露天煤矿区的适用性分析

2023-01-05惠嘉伟白中科刘凯杰王子昊

工程科学学报 2023年1期

惠嘉伟，白中科，刘凯杰，王子昊

1) 中国地质大学（北京）土地科学技术学院，北京 100083 2) 自然资源部土地整治重点实验室，北京 100035 3) 自然资源部矿区生态修复工程技术创新中心，北京 100083

煤炭在我国能源消费结构中占据主体地位，煤炭生产是我国重要的能源产业.2020 年，全国原煤产量达39 亿t，较上年增加1.4%.大量的煤炭开采不可避免地会对矿区的生态环境造成严重的破坏[1-2]，对采煤矿区进行植被恢复已经成为各主要国家的共识[3].对矿区植被进行监测是矿区植被恢复工作的重要组成部分，对矿区植被工作的设计实施与管理维护都具有重要作用[4].

归一化植被指数（NDVI）是目前在生态监测工作中最常用的植被指数[5-6].NDVI 虽然可以反映地表植被覆盖情况，但其作为一个指数并不具有明确的生态学含义.目前的研究中，大多采用植被指数转换模型将NDVI 转化为植被覆盖度后用于各类生态监测工作[7].但植被转化模型的本质是对NDVI 进行线性运算，其反映植被分布的功能与NDVI 完全相同.NDVI 植被覆盖度模型被广泛运用于矿区植被监测工作[7-12].此外，NDVI 也常被用于各类其他生态环境监测工作[13-18].

NDVI 的特点是对于低植被覆盖度非常敏感，但在植被覆盖度很高时，其检测能力逐渐下降[6,17].随着遥感技术的发展，各波段反射率被以不同形式进行组合来消除外在的影响因素，如遥感器定标、大气、观测和照明几何条件等，发展形成了多种植被指数[6]，其中有很多植被指数是基于NDVI的原理或形式进行改进而来的.邱洁等[12]采用植被恢复光谱指数（NBR）对幕府山矿区植被恢复工作进行了监测和评价.徐霞[19]等构建了温度植被干旱指数（TVDI）与NDVI 相配合对布尔台矿区进行了生态监测.Wang 等[20]基于改进型方根植被指数（GRNDVI）对巴音孟克纳源矿区的植被采矿及恢复过程进行了监测与评价.

在生态监测过程中，针对研究区域选择合适的评价指数是生态环境监测工作的基础[21-23].鉴于NDVI 在生态监测工作中的广泛使用[21]，对其在不同自然区域下的适用程度和误差进行研究显得尤为必要，这对矿区植被监测方法的完善具有借鉴意义，对矿区植被监测工作具有参考价值.

1 研究区概况

胜利矿区位于内蒙古自治区锡林浩特市西侧，煤炭年产量2000 万t，是全国储量最大的褐煤煤田.本文研究区域主要为胜利煤矿西二号露天矿及周边矩形自然植被，地理坐标为115°52′43″E～115°58′28″E，43°55′58″N～43°59′20″N，总面积约42.36 km2.研究区周边地貌以高平原为主，属于中温带干旱半干旱大陆性气候，年均降水量309 mm，主要植被为温带干旱草原.

同时采用平朔矿区安家岭露天矿、安太堡露天矿采坑、排土场及周边自然植被形成的矩形区域作对比，地理坐标为112°17′47″E～112°26′23″E，39°26′36″N～39°31′41″N，总面积约108.34 km2.平朔矿区位于山西省朔州市西北侧，属于温带半干旱季风气候区，年均降雨量450mm.研究区内地貌多为黄土低山丘陵，地带性植被为温带半干旱草原和半湿润森林草原.研究区位置见图1.

图1 研究区位置Fig.1 Location of the study area

2 数据来源与研究方法

2.1 数据来源

本文采取遥感波段反演的方式计算各植被指数，选取的遥感卫星为Sentinel-2 系列卫星.Sentinel-2 卫星分为Sentinel-2A 和Sentinel-2B 两颗卫星，重访周期为5 d.Sentinel-2 卫星携带一枚多光谱成像仪（MSI），具有13 个光谱波段，空间分辨率分别为10、20 和60 m，具有三个红边波段，适合植被健康监测.本文选取胜利矿区2019 年和2020 年7～9 月共6 景遥感数据进行植被情况监测，并选取平朔矿区2 景遥感数据进行对比与讨论，遥感数据的详细信息见表1.遥感数据来源为欧洲航天局（https://scihub.copernicus.eu/dhus/#/home，访问时间：2021 年7 月20 日）.

表1 遥感数据信息Table 1 Remote sensing data information

遥感数据预处理包括辐射定标、大气校正等，预处理步骤主要借助哨兵数据应用平台软件（SNAP）和Sen2Cor 软件（欧洲航天局，2020）进行，然后根据研究区范围将其裁剪为合适大小的区域.

2.2 研究方法

2.2.1 NDVI 计算

本文使用遥感波段反演的方法计算研究区的NDVI 值.计算公式为[5]：

式中，NDVI 为归一化植被指数；NIR、R分别为地物在近红外、红波段的反射率，在本文使用的Sentinel-2 数据中，二者分别为波段8 和波段4.

2.2.2 矿区生态监测

采用植被指数转换模型将NDVI 转化为可比性较强的植被覆盖度.公式为：

式中，f为植被覆盖度；NDVI 为归一化植被指数；NDVImin为裸地的NDVI 值；NDVImax为完全植被覆盖地区的NDVI 值.

每个研究区的两期植被覆盖度进行差值处理，即可得到在这两个年份之间研究区内植被覆盖的变化的空间分布.

2.2.3 NDVI 特性研究

采取实证对比的方式验证NDVI 在矿区植被监测工作中表现出的特征.收集平朔矿区的同来源近似时间的遥感影像，并采用相同的方法计算平朔矿区植被覆盖变化，以验证NDVI 在其他地区的矿区植被覆盖监测工作中是否会表现出相同的特征.

3 结果

3.1 胜利矿区植被覆盖监测

图2 是胜利矿区2019 年、2020 年的植被覆盖度分布.由图2 可知，2019 年胜利矿区植被覆盖度情况总体较差，7～9 月NDVI 均值分别为0.136、0.182 和0.163，2020 年7～9 月NDVI 均值分别为0.186、0.261 和0.240；2019 年7～9 月植被覆盖度均值分别为0.421、0.438 和0.421，2020 年7～9 月植被覆盖度均值分别为0.442、0.510 和0.490.从NDVI 的植被覆盖度可以反映出胜利矿区的月际植被变化及年际植被变化：就月份而言，8～9 月是研究区植被覆盖最高的时间段，7 月相对较低；就年际而言，相较于2019 年，2020 年研究区内矿山、排土场及周边自然植被覆盖度均出现了明显的上升，表明生态环境出现了明显的改善.

图2 胜利矿区2019、2020 年夏季植被覆盖度分布.（a）2019 年7 月；（b）2019 年8 月；（c）2019 年9 月；（d）2020 年7 月；（e）2020 年8 月；（f）2020 年9 月Fig.2 Distribution of vegetation coverage in the summer of 2019 and 2020 in Shengli Coal Mine: (a) July 2019;(b) August 2019;(c) September 2019;(d) July 2020;(e) August 2020;(f) September 2020

植被密度的空间分布方面，植被覆盖度较低的区域主要分布在研究区中部、北部排土场及采坑边缘区域.植被覆盖度较高的区域主要分布在矿区南部及西南部的排土场区域，成规则的条带状分布，具有明显的人工干预的特征.矿区周边的自然植被的植被覆盖度明显低于矿区南部条带状复垦区域，这表明人工干预可以有效恢复矿区生态质量.矿区南部具有较高的植被覆盖度的人工复垦区域在2020 年出现了明显的扩张，其扩张范围在8 月最为明显.由矿区南部植被优化区域的位置可以判断研究区的土地复垦与植被修复工作是自南向北进行的.

3.2 植被覆盖误差区域

经过植被覆盖分布与遥感影像的对比，可以得到植被覆盖度在空间分布中的一个异常现象.在矿区及排土场内部出现了两个较为明显的植被覆盖度较高的区域.此现象在2019 年7～9 月及2020 年9 月的研究均较为明显，在2020 年7、8 月的对应区域也有类似现象，但并不明显.采用目视解译的方式对这两个区域的范围进行大致圈定.这个植被覆盖度计算结果远高于周边自然植被，与矿区南部人工干预下的植被恢复区域水平相当.异常现象与遥感影像对比见图3.

图3 基于NDVI 的植被覆盖度误差区域.(a)植被覆盖度分布；（b）遥感影像(2020 年9 月19 日，真彩色合成)Fig.3 Vegetation coverage error area based on NDVI: (a) vegetation coverage distribution;(b) remote sensing image (September 19,2020,true color synthesis)

结合目视解译可知，基于NDVI 的植被覆盖度识别工作在自然植被、未修复排土场和已修复排土场等具有植被覆盖的区域识别比较准确，但在部分无植被覆盖的区域识别结果是完全错误的.这些区域具有以下特征：（1）集中于矿坑内部被煤炭完全覆盖的地区；（2）由于被煤炭完全覆盖，其光谱特征在可见光波段往往反射率极低，在遥感影像中呈现为黑色；（3）根据矿坑形状特点，这些区域的形状一般比较规则，或具有一定的纹理特征；（4）数值不稳定，植被覆盖度计算结果容易出现大幅度异常的上升或下降.

3.3 NDVI 的矿区归一化误差

借助ENVI5.3 抽样测定误差区域、人工恢复区域和自然植被的光谱曲线，以确定出现误差区域的原因.不同区域的部分光谱曲线见图4.因为本文涉及的Sentinel-2 数据的红、近红外波段的波长分别为665 和842 nm，因而只截取了500～1000 nm的光谱曲线.

图4 胜利矿区各区域抽样光谱曲线.（a）误差区域；（b）人工恢复区域；（c）自然植被区域Fig.4 Sampling spectrum curve of each area in Shengli mining area: (a) error area;(b) artificial restoration area;(c) natural vegetation area

由图4（a）可知，误差区域在可见光、近红外波段的反射率仅为其他区域反射率的1/10，且随着波长增大，其反射率也在逐渐增高，表现出除拐点外与自然植被近似的光谱特征.而人工恢复区域具有植被所具有的典型的“绿峰”和“红谷”，在近红外波段反射率迅速升高[24]，如图4（b）所示.此外，植被覆盖度较人工恢复区域较低，受到土壤背景值的影响更大，因而草原自然植被并没有“红谷”现象，而是反射率随着波长增大不断增大，与误差区域的光谱曲线类似，见图4（c）.

绿色植物光谱曲线在可见光波段有一个小的反射峰，两侧有两个吸收带，即在蓝和红波段为低谷，这是因为叶绿素对蓝光和红光吸收作用强，而对绿光反射作用强造成的.但是在近红外波段有一反射的“陡坡”[25].NDVI 的设计就是基于绿色植物在红波段和近红外波段的“陡坡”，利用地物在红波段和近红外波段的反射率之差（NIR-R）反映地表覆盖情况，并利用地物在红波段和近红外波段的反射率之和（NIR+R）进行归一化，归一化过程的设计目的是为了消除季节性太阳角度的差异，并最小化大气衰减的影响[5].

选取“NDVI 的矿区归一化误差”最明显的时间阶段，2019 年8 月和2020 年9 月，分别在3 个区域抽样计算NIR-R和NIR+R的均值，结果见表2.由表2 可知，NIR-R从大到小依次为：人工恢复区域、自然植被和误差区域；NIR+R从大到小依次为：自然植被、人工恢复区域和误差区域.NIRR的结果与实际植被覆盖基本相符，但归一化过程改变了不同区域间NIR-R的相对大小，导致误差区域的NDVI 值出现明显异常.同时，误差区域的NIR-R和NIR+R相对其他区域都非常小，导致该区域基于NDVI 的植被覆盖度的年际变化非常不稳定，容易出现研究区内植被覆盖度变化的最值.

表2 不同区域NDVI 计算参数均值Table 2 Calculation parameters of NDVI in different areas

这可能是因为煤矿地表覆盖特性，NDVI 在草原煤矿植被监测工作中有可能在裸煤覆盖区域呈现误差的计算结果，主要表现为矿坑内部NDVI 值异常偏高，被识别为中高植被覆盖区域.本文将这种现象称为“NDVI 的矿区归一化误差”.

对“NDVI 的矿区归一化误差”出现的原因进行分析，主要有：（1）煤矿矿坑内部裸地区域的光谱曲线与低覆盖植被的光谱曲线具有相似的特征，近红外波段的反射率均远大于红波段；（2）NDVI的归一化过程改变了煤矿矿坑内部裸地区域和低覆盖植被区域的NIR-R值的相对大小，导致出现了矿坑内部裸地的NDVI 值大于草地NDVI 值的现象.

4 讨论

4.1 胜利矿区与平朔矿区NDVI 特征对比

为了验证“NDVI 的矿区归一化误差”的科学性，有必要在其他矿区进行相同的实验验证.选取与胜利矿区的遥感影像成像时间相近的Sentinel-2平朔露天矿遥感影像，采用相同的方式计算NDVI，结果见图5.

图5 平朔矿区NDVI 分布.（a）2019 年；（b）2020 年Fig.5 NDVI distribution in Pingshuo mining area: (a) 2019;(b) 2020

从遥感影像看到，平朔露天矿的NDVI 的高值主要分布在已经实施植被恢复工程多年的矿区西侧排土场处，低值主要分布在未实施植被恢复工作的排土场内部和建设用地处，新复垦排土场的NDVI 值介于二者之间.

与胜利矿区出现误差区域不同，平朔矿区并没有明显的大范围误差区域.但通过与遥感影像的对比发现，在未修复排土场内部同样具有少量具有明显线性纹理特征的异常区域，这些异常区域大多为采煤运输通道或采煤矿坑作业面，地表没有植被生长，完全被煤炭覆盖呈现黑色.值得注意的是这些区域的NDVI 值明显高于周边具有少量植被的未修复排土场，与已修复草地接近.部分异常区域与人工修复草地区域对比见图6.

图6（c）、（d）和（e）中的具有明显纹理的黄色条带状区域即为误差区域，其NDVI 值明显高于周边.其纹理特征与遥感影像中完全被煤炭覆盖的道路、矿坑完全一致.但不同于胜利矿区，平朔矿区的误差区域呈线状纹理分布，宽度较小无法划定边界，所以采用抽样的方式计算误差区域即其他区域的NDVI 均值.在误差区域、植被人工恢复完成区域、近期植被恢复区域及自然植被区域分别各抽取10 个样点测定其NDVI 均值，结果为：误差区域0.283、近期植被恢复区域0.294、植被恢复完成区域0.802、自然植被区域0.825.误差区域和以低密度草地为主的近期植被恢复区域的NDVI值相近，远低于以林地为主的植被恢复完成区域和自然植被区域.

图6 2020 年平朔矿区NDVI 误差区域与人工修复区域对比.（a）2020 年平朔矿区NDVI 分布；（b）近期植被恢复区域；（c）误差区域示例1；（d）误差区域示例2；（e）误差区域示例3Fig.6 Comparison of NDVI error area and manual restoration area in Pingshuo mining area in 2020: (a) NDVI distribution in Pingshuo mining area in 2020;(b) recent vegetation restoration area;(c) error area example 1;(d) error area example 2;(e) error Area example 3

选取平朔矿区4 个区域抽样分析其光谱曲线，结果见图7.对比分析图4 与图7 各区域的光谱曲线的特征.两个地区的误差区域的光谱曲线类似，在可见光和近红外波段反射率总体较低.随着波长增加，其反射率逐渐增大，在近红外波段达到最高随后出现波动，近红外波段的反射率远大于红波段.平朔矿区近期植被恢复区域和胜利矿区的自然植被区域在可见光至近红外波段光谱曲线特征相近，二者自然植被均以低密度草地为主，受到土壤背景的影响，导致其“绿峰”和“红谷”现象均不明显.总体上呈现随着波长增加反射率逐渐增大，近红外波段反射率远大于红波段的特征，这个特征与误差区域十分相似.因此只依据NDVI 的植被覆盖度测算容易导致误差区域被判定为低密度草地.

图7 平朔矿区各区域抽样光谱曲线.（a）误差区域；（b）植被恢复完成区域；（c）近期植被恢复区域；（d）自然植被区域Fig.7 Sampling spectrum curve of each area in Pingshuo mining area: (a) error area;(b) vegetation restoration completed area;(c) recent vegetation restoration area;(d) natural vegetation area

综上，“NDVI 的矿区归一化误差”导致的误差区域在胜利矿区和平朔矿区均存在.误差区域的光谱曲线与受到土壤背景值影响的草地光谱曲线特征类似，计算得到二者的NDVI 值也相近.这导致在基于NDVI 的矿区植被监测工作中，误差区域（被煤炭完全覆盖的裸地）与草地容易产生混淆.对于不同研究区域，该现象的影响程度不同.胜利矿区位于干旱草原区，研究区内主要植被为牧草，植被覆盖度不高.误差区域容易被误判为草地，被归类为研究区内植被覆盖度较高的区域，对矿区植被监测工作影响较大.而平朔矿区位于季风区，研究区内优势植被为乔木、灌木，植被覆盖度较高.虽然误差区域的植被覆盖度被误判为比实际情况稍高，但误判后的结果也仅接近于低覆盖草地，在研究区内属于低植被覆盖区域，在实践中往往被视为计算误差，并不会造成较大影响.因此，“NDVI 的矿区归一化误差”对以草本植物为主的干旱区草原矿区影响较大，且影响区域较为集中，这些区域植被覆盖度几乎为0，进行掩膜处理并不会对矿区植被监测工作产生较大影响.因此建议在矿区植被监测的实践中根据上述特点对可能存在的误差区域进行掩膜处理.

4.2 其他研究成果中对比

NDVI 被广泛运用于矿区生态监测，如果“NDVI的矿区归一化误差”是存在的，那么误差区域也可能会出现在其他研究成果中.例如：图8 是帕提古丽·如则等[2]利用NDVI 对准北煤田和什托洛盖矿区的植被覆盖进行反演的结果与遥感影像的对比.图8（a）和（d）为文章中的原图，其使用NDVI对矿区地表覆盖情况进行反演.图8（d）中红框区域被识别为4 级高覆盖度植被，对比图8（a）中红框区域，二者位置及形状相似，可能是一个 “NDVI的矿区归一化误差”区域.

图8 其他研究成果中的“NDVI 的矿区归一化误差”现象[2].（a）原文研究地区；（b）Landsat 8OLI 影像（2019 年6 月1 日，真彩色合成）；（c）Google Earth 高清影像（2019 年8 月27 日）；（d）原文植被覆盖度分布（2019 年）；（e）NDVIFig.8 Phenomenon of "NDVI mining area normalization error" in other research results[2]: (a) research area in the citation;(b) landsat 8OLI image (June 1,2019,true color composite);(c) Google earth image (August 27,2019);(d) distribution of vegetation coverage in citations (2019);(e) NDVI

为了验证该结论是否成立，本文选取与原文相同来源的2019 年6 月1 日的准北煤田和什托洛盖矿区Landsat 8OLI 数据并计算NDVI，如图8（e）.为了确定该区域的地表覆盖性质，本文选取了该地区2019 年8 月27 日的Google Earth 高清遥感影像，如图8（c）.由图8（c）可知，该区域的地表覆盖以裸煤为主，形状位置与图8（a）和（d）高度相似.由图8 可知，在帕提古丽·如则等人的准北煤田和什托洛盖矿区的矿区植被覆盖反演过程中，部分影像中裸煤覆盖区域是研究区内的NDVI 值的高点，出现了典型的“NDVI 的矿区归一化误差”现象，可以证明该现象在矿区生态监测工作中是存在的并且具有一定的影响.

4.3 NDVI 适用性分析

综合对比本文研究区域与其他文献，可以发现“NDVI 的矿区归一化误差”现象随着纬度的降低逐渐变弱，在本文涉及的3 个研究区域中以平朔矿区影响最低.随着纬度降低，区域主要植被类型逐渐向乔木转变，低密度草地逐渐减少，“NDVI的矿区归一化误差”的影响逐渐减弱，所以NDVI的植被覆盖监测比较适合中低纬度地区露天煤矿.

此外，还有其他因素会影响NDVI 对于露天矿区的反演精度，例如在混合像元现象较为严重的城镇区域、大范围绿色地物覆盖区域，NDVI 的反演精度会受到严重影响.在相关研究区的遥感监测工作建议更换遥感指数进行.