基于时序多光谱影像的干旱草原区开采扰动信息提取方法
李晶, 邓晓娟, 杨震, 刘乾龙, 王媛, 崔绿园
中国矿业大学(北京)地球科学与测绘工程学院, 北京 100083

作者简介: 李 晶, 1975年生, 中国矿业大学(北京)地球科学与测绘工程学院教授 e-mail: lijing@cumtb.edu.cn

摘要

露天开采会彻底改变原有土地利用景观格局, 直接破坏当地生态环境, 甚至会影响附近居民的生产和生活, 因此越来越多的学者开始关注开采扰动。 先前有关利用时序多光谱影像提取开采扰动的研究区集中于扰动形式单一的森林区。 而我国露天煤矿大多集中于草原区, 且我国东北部的草原矿区因其脆弱的生态环境以及其他多种扰动形式的存在, 使得开采扰动识别更加困难。 为明确我国东北部生态脆弱区草原露采场的开采扰动, 以胜利矿区为例, 利用1986年—2017年27期Landsat多光谱遥感影像, 基于归一化植被指数NDVI(normalized difference vegetation index)的长时间序列轨迹变化特征(为了去除物候、 云和阴影等对时序多光谱影像的影响, 利用BISE-WT滤波器对原始NDVI时间序列进行滤波处理, 有效地去除时序NDVI数据中的噪声并同时保留有效信息), 经过样本点训练, 获得CV阈值(变异系数coefficient of variation)和Max阈值(植被阈值), 构建CV-Max扰动识别模型, 提取研究区的扰动分布。 并利用植被阈值, 分析NDVI时序轨迹, 获得扰动年际信息, 重构扰动历史地图; 进而通过分析研究区典型地物的光谱特征, 构建裸煤提取规则, 以此来提取研究区的裸煤分布; 最后通过构建裸煤及扰动区两者间的拓扑关系, 进行空间拓扑叠置分析, 从而获得开采扰动信息。 经过精度验证, 开采扰动的提取精度达到93.17%(Kappa系数=0.85), 扰动年际信息提取精度达到83.35%(Kappa系数=0.81)。 结果表明: 在研究期间, 空间上, 开采扰动面积占研究区总面积的8.90%; 时间上, 开采扰动的发生集中于2000年—2009年, 期间开采扰动像元占开采扰动总像元的76.70%; 1988年—1998年矿区属于土地损毁初始期, 2000年—2005年矿区属于土地损毁加速期, 2006年—2009年矿区属于土地损毁高峰期, 2010年—2017年开采扰动像元占比趋势比较平缓且持续处于较低水平, 矿区土地损毁范围基本稳定。 所提出的针对我国东北部生态脆弱性草原矿区, 基于时序多光谱影像, 利用植被指数NDVI和裸煤光谱特征提取开采扰动信息的方法是可行的, 该研究结果可为干旱、 半干旱草原露天矿区的可持续发展提供数据和理论方法支撑。

关键词: Landsat时序影像; NDVI; 光谱特征; 开采扰动; CV-Max模型
中图分类号:TP753 文献标志码:A
A Method of Extracting Mining Disturbance in Arid Grassland Based on Time Series Multispectral Images
LI Jing, DENG Xiao-juan, YANG Zhen, LIU Qian-long, WANG Yuan, CUI Lü-yuan
College of Geoscience and Surveying Engineering, China University of Mining and Technology(Beijing), Beijing 100083, China
Abstract

Surface mining will completely change the original landscape pattern of land use, directly destroy the local ecological environment, and even affect the production and life of the nearby residents; therefore, more and more scholars have begun to pay attention to mining disturbance. Previous studies on extracting mining disturbances from temporal multispectral images focused on forest areas with single disturbance form. However, most surface mines in China are concentrated in grassland areas, and in the grassland mining areas in northeastern China, mining disturbances are more difficult to be identified because of their fragile ecological environment and the existence of various other forms of disturbance. In order to clarify the mining disturbance of grassland open stope in ecologically fragile areas in northeastern China, the authors taking Shengli mining area as an example, firstly employs 27 Landsat multi-spectral remote sensing images from 1986 to 2017, and bases the study on the long time series trajectory change characteristics of NDVI (normalized difference vegetation index). (In order to remove the effects of phenology, cloud and shadow on time series multispectral images, BISE-WT filter is used to filter the original NDVI time series to effectively remove the noise in the time series NDVI data and retain the effective information at the same time). After sample point training, CV threshold (Coefficient of Variation) and Max vegetation threshold are obtained. The Max vegetation threshold (vegetation threshold) is then used to construct the CV-Max disturbance recognition model and extract the disturbance distribution in the study area. Furthermore, using vegetation threshold, NDVI time series trajectory is analyzed to obtain disturbance interannual information and reconstruct disturbance history map. Then, by analyzing the spectral characteristics of typical terrain in the study area, bare coal extraction rules are constructed to extract the distribution of bare coal in the study area. Finally, the topological relationship between bare coal and disturbance area is constructed and a spatial topological overlay analysis is conducted to obtain mining disturbance information. The accuracy verification reveals the extraction accuracy of mining disturbance is 93.17% (Kappa coefficient=0.85) and the extraction accuracy of disturbance interannual information is 83.35% (Kappa coefficient=0.81) respectively. The results show that during the study period, the mining disturbance area accounts for 8.90% of the total area of the study area in space; in terms of time, the occurrence of mining disturbance concentrated in 2000—2009, during which the mining disturbance pixels accounted for 76.70% of the total mining disturbance pixels; the years from 1988 to 1998 witness the initial period of land destruction, and in 2000—2005, land destruction increased in the mining area, and in 2006—2009, the land destruction the mining area reached the peak. The proportion trend of mining disturbance pixels in 2010—2017 is relatively flat and continues to be at a low level, and the scope of land damage in mining area is basically stable. In view of the ecologically fragile grassland mining area in northeastern China, the method of extracting mining disturbance information by using NDVI and bare coal spectral features based on time series multispectral images is feasible. The research results can provide data and theoretical method support for the sustainable development of arid and semi-arid grassland surface mining area.

Keyword: Landsat multi-temporal; NDVI; Spectral features; Mining disturbance; CV-Max model
引 言

中国露天煤矿多处于干旱、 半干旱草原生态脆弱区, 煤炭开采会彻底改变原有土地利用景观格局, 对矿区生态环境造成破坏狭义的。 开采扰动区是指露天开采过程中导致的直接表土层剥离所裸露的露天采场、 排土场、 井工矿区塌陷地或煤矸石堆放地, 矿区工业场地等一系列采矿活动直接影响的空间区域。 因此, 提取开采扰动的空间位置及时间信息, 可为煤炭开采及矿区土地复垦制定正确决策提供科学依据和理论支撑。

传统获取地表覆盖信息的方法主要依据野外实际测定和调查分析, 工作量大、 经费多、 周期长以及效率低等问题凸显。 近年来, 随着遥感数据的不断积累及相关理论与技术的不断革新, 使得基于遥感影像的中长期变化检测技术兴起。 在时序遥感获取信息的技术上, 许多国内外学者在森林扰动及土地利用变化方面进行了大量研究[1, 2], 提出了基于LandTrendr算法[3]、 BFAST算法[4, 5, 6]、 DI指数法[7]、 VCT算法[8]、 CMFDA算法[9]等多种方法进行森林干扰和恢复趋势检测。 时序遥感在开采扰动信息提取方面也获得了新的进展, Sen[10]等利用1984年— 2008年的Landsat TM/ETM+遥感影像, 对美国App alachia矿区开采扰动进行识别, 研究结果表明: NDVI和TC G/B等遥感指数能较好地用于识别开采扰动, 面向对象的分类方法可以提高开采扰动识别精度。 Li[11]等通过选取特征阈值, 基于NDVI长时间序列轨迹分析了像元尺度的土地损毁和复垦过程特征, 将土地类型分为持续森林覆盖、 开采扰动和其他扰动三类, 并提取扰动区与扰动年际信息。

上述研究多利用时序光谱遥感影像探索森林扰动或者森林矿区开采扰动信息的提取, 鲜有干旱、 半干旱草原矿区的研究。 而干旱、 半干旱草原矿区多处于生态脆弱区, 且存在包括城镇建设、 工业开发等其他扰动, 导致利用时序光谱遥感影像识别干旱、 半干旱草原矿区的露天煤矿开采扰动更加困难。 针对这一问题, 本文以胜利露天煤矿为例, 利用1986年— 2017年27期Landsat TM/ETM+/OLI多光谱遥感影像, 构造了该区的长时间NDVI序列, 结合扰动像元的轨迹变化特征、 地物光谱特征及对象间的拓扑关系, 准确地提取了研究区的开采扰动区以及开采扰动年际信息。 研究结果可为矿区生态重建提供数据支撑和科学决策。

1 研究区概况和数据来源
1.1 研究区概况

胜利煤田位于内蒙古自治区锡林郭勒盟锡林浩特市西北部胜利苏木境内, 地处内蒙古高原的东部, 大兴安岭西延的北坡。 地理坐标为东经115° 30'— 116° 26', 北纬43° 57'— 44° 14', 见图1所示。 属于中温带半干旱、 干旱大陆性季风气候, 草原资源丰富, 草场种类繁多。 坐落于我国生态安全“ 两屏三带” 的北部放沙地区, 受气候条件制约和地区产业发展影响, 具有酷寒、 半干旱、 土壤贫瘠等生态脆弱性特征。 年平均气温14℃, 年平均降水150400 mm, 由东南向西北递减; 年平均相对湿度在60%以下, 蒸发量为1 5002 700 mm, 由东向西递增; 年日照时数2 8003 200 h, 无霜期为90120 d, 全年平均风速3.55.3 m· s-1; 研究区总面积653 km2, 煤田含煤面积342 km2, 地质储量22 442 Mt, 是全国最大、 煤层最厚的褐煤田。

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

1.2 数据来源

1986年以来的Landsat TM/ETM+/OLI多光谱遥感影像数据, 来自于美国地质调查局(USGS)资源探测数据集(http://glovis.usgs.gov/), 空间分辨率为30 m, 影像轨道行列号为P124/R29, 每年优选一期云层遮盖小于10%的植被主要生长期(数据采集时间多为7月— 9月)影像, 1990年、 1999年、 2001年、 2012年及2016年植物生长季采集的影像均为大范围云层覆盖, 32年间共选用27期影像, 影像获取类型及日期见表1

表1 影像类型及获取日期 Table 1 List of Landsat image dates and their sensor types
2 实验部分

本研究具体技术流程见图2。 对所选的Landsat多光谱遥感数据进行预处理以后, 计算每期影像的植被指数(NDVI), 通过BISE-WT[12]滤波器处理NDVI时间序列, 获得高质量的NDVI时间序列。 经过样本点训练, 利用CV阈值和最大植被阈值构建CV-Max[13]扰动提取模型, 以此提取扰动像元。 利用植被阈值, 分析NDVI时序变化轨迹获取扰动年际信息[11]; 分析裸煤光谱特征, 构建裸煤提取规则, 实现利用光谱遥感影像自动化提取裸煤[14]。 最后通过构建两者间的拓扑关系, 进行空间拓扑叠置分析, 提取开采扰动信息。

图2 流程图Fig.2 Flow chart

2.1 数据预处理

对Landsat TM/ETM+/OLI多光谱遥感影像逐一进行波段合成后, 进行云及阴影掩膜等, 应用矿区边界矢量文件进行研究区边界裁剪。

植被指数NDVI是植被生长状态及植被覆盖度的最佳指示因子, 常广泛应用于检测矿区土地损毁变化方面[15, 16]。 计算研究区每期影像的NDVI, 其计算公式见式(1)

NDVI=ρNIR-ρREDρNIR+ρRED(1)

式中: ρ NIRρ RED分别代表近红外波段和红波段的反射率, NDVI的值介于-1和1之间。

将获得的1986年— 2017年NDVI时间序列, 利用BISE-WT滤波器进行滤波处理, 该滤波器是由Yang等于2018年提出[12], 该方法结合最佳斜率提取算法(BISE)和小波变换(WT)对年际Landsat NDVI时序数据进行去噪处理。 图3为BISE-WT滤波示例, 可以看出BISE-WT滤波器在降噪的同时保留有效信息, 得到高质量NVDI时间序列。

图3 不同时序NDVI经BISE-WT滤波前后对比图Fig.3 Contrast charts of NDVI with different sequences before and after BISE-WT filtering

2.2 扰动区与扰动年际信息识别

2.1.1 扰动区识别

利用CV阈值(变异系数)和Max阈值(在本文被解译为植被阈值)构建CV-Max扰动提取模型。 Max阈值是非植被与植被的NDVI临界值, 用来确定某个像元属于植被或非植被。 CV阈值可以反映NDVI时间序列的离散程度, 用来确定某个像元在研究期间的任何时间点是否受到过扰动, 计算公式见式(2):

CV(x, y)=std(NDVI(x, y))mean(NDVI(x, y))(2)

式中: (x, y)是样本点位置; std是某个样本点时间序列的标准偏差; mean是某个样本点时间序列的NDVI平均值。

为了获取CV阈值和Max阈值, 在Landsat多光谱遥感影像上随机选择386个样本点, 通过实地调查结合Google Earth目视判读确定每个样本点属于植被或者非植被。 并对于判定为植被的点进一步判断在研究期间是否受到过扰动。 获取样本点的NDVI值, 并作为训练样本点, 得到研究区的植被阈值为0.28, 因此, 所有NDVI值> 0.28的像元都被解读为植被, 反之为非植被。 CV阈值为0.235, 所有CV值> 0.235的像元被认为在研究期间植被受到破坏或之前的某个时间点被破坏而研究期间植被经历恢复。

基于CV阈值和植被阈值, 构建CV-Max分类模型, 将研究区的每一个像元点归类为低植被或非植被、 未受干扰的植被与受干扰的植被三种类型中的一种, 具体分类如图4所示。 图4(a)代表整个研究期间无植被覆盖或低植被覆盖的像元, NDVI时序轨迹始终处于比较低的状态; 图4(b)代表整个研究期间始终没有受到过扰动的像元, NDVI时序轨迹始终处于比较高的状态; 图4(c)和(d)代表在研究期间或之前受到过扰动的像元, NDVI时序轨迹变化幅度较大, 这部分像元就是CV-Max模型所要提取的扰动区。

图4 CV-Max分类模型Fig.4 CV-Max classification model

2.1.2 扰动年际信息识别

将典型开采扰动NDVI时序轨迹分为4个阶段, 见图5所示。 第一阶段: t0-t1是开采前, 地表植被未被破坏, NDVI值始终保持在较高水平并仅在小范围内波动; 第二阶段: t1-t2表示地表植被受到开采干扰, 植被去除地表裸露, NDVI明显下降; 第三阶段: t2-t3表示开采期间以及复垦前, NDVI始终低于植被阈值并处于相对较低的水平, 整体趋势缓慢上升; 第四阶段: t3-t4表示植被恢复阶段, NDVI高于植被阈值并逐渐上升。

图5 典型开采扰动NDVI时序轨迹Fig.5 Typical mining disturbed NDVI trajectory

由于开采扰动事件发生于t1-t2时间段内, 属于采矿活动初始阶段, 由于遥感影像重访周期、 研究区气候及人类活动的限制, 无法精确地确定开采扰动年际信息, 所以将扰动时间定义为NDVI首次由大于植被阈值变为小于0.28的年份。

2.3 开采扰动识别

通过分析研究区典型地物(植被、 建设用地、 裸土、 水体及裸煤)的光谱特征, 如图6所示, 可以得知裸煤(指裸露在地表的煤层或煤堆)区别于其他地物的光谱特征规则, 见式(3)— 式(5)。 基于上述规则自动化提取每景Landsat TM/ETM+/OLI影像的裸煤, 然后将27期影像的裸煤进行叠加即为后续分析所用的裸煤。

图6 各地物光谱特征Fig.6 Spectral characteristics of various objects

ρSWIR1-ρNIR> 0(3)ρSWIR2-ρSWIR1> 0ρSWIR2< 0.15(5)(4)

综合利用Google Earth、 空间拓扑及空间叠加分析构建空间拓扑叠加分析提取开采扰动, 具体步骤为:

(1)利用ArcGIS将扰动区与裸煤转为面矢量文件;

(2)将扰动区和裸煤区面文件进行叠加, 通过目视分析和野外调查, 若一个扰动面包含一个或者多个大范围裸煤多边形(面积> 0.1 km2), 确定该扰动为开采扰动区; 剩余扰动明确为其他扰动。

3 结果与讨论
3.1 精度评价

使用混淆矩阵来验证该方法提取的开采扰动信息精度。 为验证开采扰动的分类识别精度, 应用ArcGIS随机生成600个验证点, 其中有235个开采扰动点以及365个其他扰动点, 所有验证点都通过遥感影像、 Google Earth以及实地调查目测验证。 本文开采扰动识别方法的总体精度达到了93.17%, Kappa系数值为0.85。 验证精度整体较高, 分类结果较理想。 错误分类像元多分布在在扰动多边形的边缘, NDVI时间序列轨迹震荡比较复杂, 不容易分辨。

为验证开采扰动年际信息的识别精度, 随机生成931个验证点, 通过验证得到开采扰动年的总体分类精度为83.35%, Kappa系数值为0.81, 分类效果较理想。 错误分类像元大多在研究期间经历两次或两次以上的扰动, 或是像元仅在某一年没有植被覆盖后续受到扰动植被持续处于较低状态, 发生识别错误。

3.2 开采扰动分析

提取的开采扰动区如图7所示。 研究区总面积为653.359 0 km2, 提取的开采扰动区面积为58.124 2 km2, 占总面积的8.90%。 由图可知: 开采扰动主要集中分布在四个区域, 经实地调查, 四个扰动聚合区是胜利煤矿的四个露天矿坑。 其他扰动包括由于物候条件影响的植被年际变化以及城镇建设等。 除露天矿之外的较大范围裸煤区经实际调查发现是洗煤厂, 所以将其明确为其他扰动。

图7 开采扰动提取Fig.7 Mining disturbance extraction

图8所示是胜利四个露天矿坑的扰动年际信息分类。 2006年和2009年开采扰动像元最多, 分别占开采扰动总像元的22.07%和23.81%。 2000年— 2009年开采扰动像元比重为76.70%, 可知, 胜利矿区开采时间基本集中于2000年— 2009年, 十年间, 采矿作业区的面积占比呈显著上升趋势, 且增幅较大。 1988年— 1998年间矿区处于土地损毁初始期, 开采扰动损毁区域增长趋势较稳定, 2000年— 2005年处于土地损毁加速期, 开采扰动像元加速增长, 2006年— 2009年矿区属于土地损毁高峰期, 2010年— 2017年开采扰动像元占比趋势比较平缓且持续处于较低水平, 矿区土地损毁范围基本稳定。 综合可知, 开采扰动区随着时间的推移逐渐扩大, 露天煤炭开采对干旱、 半干旱草原矿区及周边的土地利用影响较严重。

图8 开采扰动年Fig.8 Mining disturbance year

总体而言, 基于时序多光谱影像的开采扰动信息提取方法在干旱、 半干旱草原地区的露天矿应用效果比较好, 可以将开采扰动与其他扰动有效区别, 并准确提取开采扰动信息。 然而, 对于井工矿开采影响边界识别, 指标的选取以及方法的应用都有待于进一步研究。

4 结 论

利用27期Landsat TM/ETM+/OLI多光谱遥感影像, 通过分析光谱指数NDVI时序变化轨迹、 裸煤遥感光谱特征以及对象间的拓扑关系, 准确提取了干旱草原矿区的开采扰动信息, 为矿区生态重建提供科学支持。 具体结论如下:

(1)经混淆矩阵精度评价, 该方法识别开采扰动的总体精度为93.17%, Kappa系数值为0.85, 扰动年际信息提取总体精度为83.35%, Kappa系数值为0.81。 整体分类精度较高, 研究表明该方法适用于干旱、 半干旱草原矿区的土地动态变化检测。

(2)空间上, 研究区在1986年— 2017年开采扰动面积为58.124 2 km2, 占总面积的8.90%; 时间上, 研究区开采发生基本集中于2000年— 2009年, 占开采扰动总像元的76.70%。 十年间, 采矿扰动区的面积及占比呈显著上升趋势, 且增幅较大。 综合可知, 开采扰动区随着时间的推移逐渐扩大, 可见煤炭开采对干旱、 半干旱草原矿区生态脆弱区的土地利用影响增强。

参考文献
[1] Michael S, Richard L, Peter B, et al. Remote Sensing of Environment, 2015, 158(1): 156. [本文引用:1]
[2] Huang C B, Dian Y Y, Zhou Z X, et al. Journal of Remote Sensing, 2015, 19(4): 657. [本文引用:1]
[3] Kennedy R E, Yang Z Q, Cohen W B. Remote Sensing of Environment, 2010, 114(12): 2897. [本文引用:1]
[4] Verbesselt J, Hyndman R, Newnham G, et al. Remote Sensing of Environment, 2010, 114(1): 106. [本文引用:1]
[5] Verbesselt J, Hyndman R, Zeileis A, et al. Remote Sensing of Environment, 2010, 114(12): 2970. [本文引用:1]
[6] Verbesselt J, Zeileis A, Herold M. Remote Sensing of Environment, 2012, 123: 98. [本文引用:1]
[7] Yang C, Shen R P, Yu D W, et al. Journal of Remote Sensing, 2013, 17(5): 1246. [本文引用:1]
[8] Huang C, Goward S N, Masek J G, et al. Remote Sensing of Environment, 2010, 114(1): 183. [本文引用:1]
[9] Zhu Z, Woodcock C E, Olofsson P. Remote Sensing of Environment, 2012, 122: 75. [本文引用:1]
[10] Sen S, Zipper C E, Wynne R H, et al. Photogrammetric Engineering and Remote Sensing, 2012, 78(3): 223. [本文引用:1]
[11] Li J, Zipper C E, Donovan P F, et al. Environ. Monit. Assess, 2015, 187(9): 557. [本文引用:2]
[12] Yang Z, Li J, Shen Y Y, et al. International Journal of Remote Sensing, 2018, 39(12): 3816. [本文引用:2]
[13] Yang Z, Li J, Zipper Carl E, et al. Science of the Total Environment, 2018, 644(10): 916. [本文引用:1]
[14] Yang Z, Li J, Yin S Q, et al. Remote Sensing Letters, 2018, 9(12): 1224. [本文引用:1]
[15] LI Jing, Zipper C E, LI Song, et al(李晶Zipper Carl E李松,)Transactions of the Chinese Society of Agricultural Engineering(农业工程学报), 2015, 31(16): 251. [本文引用:1]
[16] LI Jing, JIAO Li-peng, SHEN Ying-ying, et al(李晶,焦利鹏,申莹莹,)Journal of China Coal Society(煤炭学报), 2016, 41(11): 2822. [本文引用:1]