基于Landsat8 与高分数据的矿山植被动态监测研究

2022-09-01吴凤敏郑稚棚张灵犀杨光谱赵珍妮

地理空间信息 2022年8期

吴凤敏，郑稚棚，余静，梁星，张灵犀，杨光谱，余洋，赵珍妮

（1. 重庆市地理信息和遥感应用中心，重庆 401147）

目前，对于植被遥感动态监测有大量研究[1-9]，其中Landsat 数据应用最为广泛，该类数据可以免费下载，且适用于大尺度范围的长时间序列监测[10-22]。基于此，研究以重庆市中梁山片区为实验区，利用Landsat8 数据对矿山植被长时间序列情况进行分析，同时结合高分数据对矿山植被短时间序列进行动态监测，并对于不同分辨率数据的植被监测结果及关联性进行分析，希望能够为后续矿山修复治理过程中植被状况动态监管提供技术支撑。

1 研究区概况及研究方法

1.1 研究区概况

中梁山片区地理坐标为106°23′20″～106°24′25″E，29°26′10″～29°31′50″N，总面积为28.63 km2，最高点位于中部山顶，最高高程为695.8 m，最低点位于山脚，最低高程为270 m。研究区内矿山有70个，总面积为191.36 hm2，位于中梁山背斜之南倾伏端，以低山地貌为主。研究选择中梁山片区为示范区域，主要由于该区域露天矿山数量较多，存在已完成修复治理、正在修复治理和关闭未治理等多种类型，在矿山植被状况动态监测研究中具有一定代表性（图1）。

图1 中梁山片区研究区域

研究利用landsat8 数据对2012—2021 年中梁山片区露天矿山植被指数进行变化分析，利用高分数据进行2017—2021 年短时间序列变化分析。Landsat8 数据来源于地理空间数据云，高分数据包括高分1号和高分6 号数据，来源于重庆市地理信息和遥感应用中心。为尽量保证植被监测年度生长一致，Landsat8数据主要采用夏天时期影像，数据获取时间分别为2013-08-19、2014-07-21、2015-07-08、2016-07-11、2017-08-14、2018-05-13、2019-03-29、2020-09-07、2021-06-06，根据1∶5 000DEM数据计算该区域平均高程为0.37 km。高分数据数据获取从2017—2021年，部分数据因为不完全覆盖实验区需要进行拼接，每景数据平均高程同样根据1∶5 000DEM进行计算（表1）。

表1 高分卫星数据使用总体情况表

1.2 主要技术流程

研究采用ENVI5.3 软件首先对Landsat8 数据和高分数据进行预处理，分别进行植被指数计算，并对矿山植被状况进行现状及变化分析，最后对两者的关联性进行详细对比分析（图2）。研究选择NDVI 作为矿山植被监测指数，由于矿山生态修复过程中植被生长变化较大，NDVI 能够较好地反映植被指数从低值到高值区间的动态变化，与植被分布密度呈线性相关。数据预处理主要包括辐射定标、大气校正、正射校正等流程。辐射定标处理主要指将影像的灰度（D）值转化成辐射亮度值或表观反射率；大气校正主要是将辐射定标后的表观反射率或辐射亮度值转化为地表反射率的过程；正射校正需基于地形高程模型（DEM）来对影像中所有的像元执行地形变形的校正，让图像满足正射投影的要求，可根据影像自带的RPC参数和ENVI中的DEM信息，对影像数据进行正射校正。

图2 Landsat8数据和高分数据处理及关联性分析技术流程

采用ENVI软件进行辐射定标时，Landsat8数据可以直接采用自带模型Radiometric Calibration进行辐射定标，高分数据需要网上下载并安装China Satellites Support插件，同时根据中国资源卫星应用中心发布的辐射定标参数（http://www.cresda.com/CN/index.shtml）。由于NDVI计算取值范围为[-1.0，1.0]，研究对该范围外的作为异常值进行剔除。根据已有学者研究成果，裸土NDVI 理论上接近于0，但由于受众多因素影响，其取值范围在-0.1～0.2之间。

2 Landsat植被动态监测

2.1 现状情况分析

参考已有研究，按照无植被覆盖区[-1.0，0.2），低植被覆盖区[0.2，0.3），中等植被覆盖区[0.3，0.4），高植被覆盖区[0.4，1.0）进行分级。结果显示，2021年中梁山片区矿山高植被覆盖区面积为81.09 hm2，占中梁山片区矿山总面积的42.38%；中等植被覆盖区面积为29.35 hm2，占比15.34%；低植被覆盖区面积为30.72%，占比16.05%，表明中梁山片区矿山植被整体情况较好，但存在部分矿山还在修复治理过程中，植被指数较小或者为推土情况无植被覆盖（表2）。高植被覆盖区中，NDVI 以0.4～0.5 区间为主，该区间矿山面积占高植被覆盖区面积的50%以上，＞0.6的面积占比为46.88%。以单个矿山来看，中梁山片区70 个矿山中，平均NDVI 超过0.4 的矿山数量有36 个，平均NDVI 为0.3～0.4 的矿山为19 个，平均NDVI 为中等的矿山共有55 个，说明随着矿山修复治理工作的不断开展，矿山的植被指数状况总体较好，大部分矿山植被恢复较好，部分矿山由于正处于修复治理过程中，可能存在部分区域推填土或者刚种植植被，植被指数较低。

表2 2013—2021年中梁山片区矿山植被指数NDVI分级表

2.2 年际变化分析

研究以每隔1 a NDVI值对矿山范围内数据进行对比分析，2013—2021年中梁山片区植被指数NDVI呈现先减小后增加趋势（图3、4），主要波谷在2017年左右，由于矿山开采时期植被破坏较大，植被指数较低，在修复治理过程中，也可能存在部分区域推填土破坏植被的情况；随着矿山不断修复，前期的植被指数从较低的覆盖指数逐渐增加，最终修复完成后植被不断恢复，植被指数不断增加，表明矿山生态修复治理取得了较好的成绩。其中，矿山高植被覆盖区（NVDI＞0.4）从2013 年的104.81 hm2下降至2017 年的42.17 hm2，减少面积为62.64 hm2；从2017—2021 年面积又逐渐增加，2021年的高植被覆盖区面积达到81.09 hm2。低植被覆盖区（0.2＜NDVI≤0.3）具有同样的规律，2013—2017 年面积减少12.64 hm2，2017 年到2021 年面积增加17.54 hm2。2017—2021年间高植被覆盖区增加的趋势较低植被覆盖区更大，说明这个期间为矿山植被恢复的主要时期，最开始以低植被覆盖的治理初期随着时间的推移，植被逐渐变好，并且变为高植被覆盖区较多。

图3 2013年和2021年中梁山片区植被指数NDVI分布图

研究发现，存在个别矿山植被指数到2021年总体情况较低，比如中梁山片区萌特矿山，经实地调查，该矿山于2019年6月开始工程治理，目前正处于修复治理过程中，项目周期为66 个月，到2025 年为止，存在部分区域因修复治理过程中对周边土地的踩踏以及修复治理最开始种植小树苗在影像上不明显，因此总体上植被指数较低。

图4 中梁山片区矿山2013—2021年植被指数变化情况

2.3 变化斜率分析

研究采用一元线性回归方法模拟每个栅格的NDVI变化趋势，重点分析矿山范围内不同修复治理区域植被变化情况。对于NDVI 时间序列数据，每个像元对应有若干年的时间序列数值，对这些数值进行线性拟合，获取直线的斜率，像元斜率表明了在该时间序列中该像素所代表的植被指数的演化趋势。主要计算公式如下：

式中，k为斜率；xi为具体年份；yi为每年的NDVI值。

2017—2021 年中梁山片区矿山范围内斜率＞0 的像元总面积为85.15 hm2，占矿山总面积的44.49%，表明这些区域内植被逐渐变好；斜率＜0 的像元面积为106.21 hm2，表明这些区域植被状况有变差的趋势。参考已有研究[13]，将NDVI 变化程度分为7 级（表3），图5 为以30 m×30 m 像元大小，统计中梁山片区矿山范围内所有像元的斜率值分布情况。中梁山片区矿山范围内像元明显改善面积有55.22 hm2，占像元总面积的28.85%；轻微改善及中度改善的面积合计为26.14 hm2，占比13.66%，存在一般以上斜率为负值，植被情况较差的情况，主要原因是中梁山片区面积最大的矿山重庆萌特矿山处于修复治理过程中。虽然有一部分区域在种植树木，但由于数据还在恢复中，植被指数变化较小，且治理过程中存在部分推填土情况，因此植被指数总体变小。

表3 2013—2021年中梁山片区矿山像元NDVI斜率分级（Landsat8数据）

图5 中梁山片区矿山范围内单个像元斜率分布情况

研究通过ArcGIS 工具将计算斜率值与中梁山片区矿山像元进行空间关联，可以看出中梁山片区中部及东部矿山修复情况较好，植被斜率大部分为正数，右下角萌特矿山和中部部分矿山目前还处于修复治理过程中或者关闭未治理情况，植被覆盖较差，如图6所示。

图6 2013—2021年中梁山片区矿山范围内像元植被指数斜率分布图

3 高分数据植被动态监测

3.1 现状情况分析

研究按照与Landsat8 数据同样分析方法对高分数据计算的植被指数进行分析，由于分辨率较高（8 m），矿山植被指数监测结果较landsat8 数据更精细。中梁山片区整体植被情况较好，大部分区域均为高植被覆盖区域（NDVI ＞0.4），70个矿山总体上位于中高植被覆盖区域范围内，2021 年有植被覆盖区域面积为150.49 hm2，其中高植被覆盖区域94.84 hm2，占矿山总面积的49.56%。与landsat8 计算结果相比，高植被覆盖区面积较大，说明高分数据对于植被指数的响应更敏感。中等植被覆盖区面积为29.48 hm2，占比为15.41%（图7）。

图7 2021年中梁山片区高分数据NDVI监测分布图

3.2 变化斜率分析

由于高分数据只有2017—2021 年5 a 的数据，因此进行短时间序列斜率变化分析（表4，图8），图8为以8 m×8 m像元大小，统计中梁山片区矿山范围内所有像元的斜率值分布情况。结果表明中梁山片区矿山范围内NDVI斜率＞0的像元面积为78.94 hm2，占矿山总面积的41.25%；斜率＜0的像元面积为112.42 hm2，总体上超过矿山总面积的一半，NDVI 斜率的分布情况与Landsat8 数据分析结论基本一致，仅部分分级数值略有差异。

图8 中梁山片区矿山范围内单个像元斜率分布情况（高分数据）

表4 2017—2021年中梁山片区矿山像元NDVI斜率分级（GF数据）

从分级情况看，矿山范围内严重退化区域面积最大，为108.51 hm2，主要是由于中梁山片区萌特矿山正在修复治理过程中造成。明显改善区域面积居第二，占矿山总面积的39.45%，表明其他大部分矿山区域的植被都是变好的，说明矿山修复治理取得了良好的效果。

4 植被指数关联性分析

研究选择2019 年数据进行关联性分析，利用Landsat8 （2019-03-29）和高分数据（2019-07-28）分别计算NDVI，得到两者空间分布图。结果显示中梁山片区范围内高分数据在高植被覆盖区域分级粒度更细，NDVI 大于0.8 在图上分布较为明显，小于0.2区域面积总体较少，主要由于高分数据分辨率为8 m，而Landsat8分辨率为30 m，在像元尺度上部分高值和低值区域进行了一定综合，因此在中间部分值的像元数量较多（图9，表5）。

表5 中梁山片区矿山2019年Landsat8和高分数据计算NDVI分级情况对比

图9 2019年Landsat8和高分数据计算NDVI空间分布图对比

通过对中梁山片区矿山范围内的两种影像计算的DNVI 进行计算得出，高分数据有植被覆盖区域（NDVI ＞0.2）为82.68 hm2，明显小于Landsat8计算出的有植被覆盖区域面积（132.93 hm2），分析造成趋势情况不一致的原因有以下几种：一是分辨率原因导致部分Landsat8 像元综合了周围的NDVI 值较高的植被区域导致像元内植被指数整体较高；二是由于两类NDVI数据在进行正射纠正时均采用ENVI自带的DEM数据，可能存在部分像元位置有一定偏移；三是两类NDVI 数据的时相不一致，对于植被生长期情况的计算结果本身也有一定差异。

研究利用中梁山片区范围对2 种数据进行植被指数关联性分析，通过ArcGIS 中Spatial Analyst Tools 工具进行计算。计算结果显示，通过Landsat8 和高分数据计算NDVI 相关系数为0.739 3，表明2 种遥感影像在计算结果上呈现较高的正相关性。同时，研究以Landsat8 像元为网格（30 m × 30 m），计算高分数据NDVI 平均值，分析空间上同一像元NDVI 值差异情况，结果表明，高分数据与Landsat8 数据计算NDVI值呈现花瓣状相关特征，两者在高值区和低值区呈聚拢状态，在中间值区呈离散状态，总体上表现为中间宽两头窄的花瓣形。以线性关系式表达两者特征，决定系数R2=0.692，表明两者虽然对于同一像元的值存在差异，但总体上还是具有较好的线性相关性（图10）。