基于LandTrendr和CCDC算法的神东煤炭基地植被损毁识别对比分析

李 军 张艺藂 张彩月 谢慧真 张成业 杜梦豪 王雅颖

(1.中国矿业大学(北京)地球科学与测绘工程学院,北京 100083;2.煤炭资源与安全开采国家重点实验室,北京 100083)

煤炭基地长时间的采矿活动对土地造成的挖损、压占、用地类型转变等,给生态环境带来了不容忽视的负面影响[1-3],其中最直接的影响之一就是对地表原有植被造成损毁。遥感具有大面积快速观测的优势,利用长时序遥感技术准确识别煤炭基地植被损毁区域和损毁时间,对当地生态环境质量监管具有重要的现实意义[4-5]。

为实现植被损毁遥感识别,国内外学者提出了多种算法。主要包括Sen+Mann-Kendall算法、趋势分析法、LandTrendr (LT)算法、Continuous Change Detection and Classification (CCDC)算法等。钟琪等[6]、王科雯等[7]分别利用Sen +Mann-Kendall算法拟合了大宁矿区和神东矿区植被的损毁趋势,并对变化速率进行了分级。徐佳等[8]、刘英等[9]利用趋势分析法,通过最小二乘原理进行逐像元的一元线性拟合,计算像元回归直线的斜率,实现了对神东矿区植被损毁的趋势识别,并分别预测了矿区植被未来变化趋势,分析了影响植被的驱动因子。但Sen+Mann-Kendall算法和趋势分析法只体现了长时序植被损毁的变化趋势,无法准确判别植被损毁的位置与时间。与之相比,LT算法和CCDC算法能够获取像元尺度的植被变化轨迹,进而确定植被损毁发生的位置、时间等信息,因此被广泛应用于植被损毁识别中。在大尺度森林背景下,殷崎栋等[10]、刘珊珊等[11]利用LT算法分别实现了陕西省和浙江省森林区域的植被损毁识别。肖真等[12]、李贺等[13]利用CCDC算法,分别实现了江西省森林区域和缅甸南部橡胶林区域的植被损毁识别。在小尺度矿区背景下,YANG等[14]、杨张瑜等[15]利用LT算法分别对澳大利亚的库拉格煤矿以及中国的紫金山金铜矿进行了植被损毁的动态识别;ZHANG等[16]利用CCDC算法实现了对江西德兴铜矿的植被损毁动态识别。但需要注意的是,即使在同一区域,不同算法进行植被损毁识别的效果不尽相同[17-18]。因此,需要对LT算法和CCDC算法在同一区域的植被损毁识别结果进行进一步的定量分析与对比,确定其适用性。

近年来,Google Earth Engine (GEE)凭借其海量的数据资源和强大的超算能力成为了快速遥感云计算的利器[19]。利用GEE云计算的优势,LT和CCDC算法计算速度得到了极大的提升,为区域植被损毁识别提供了大规模运用的机会。但是,以往矿业领域的研究主要是针对单个矿区,缺少对LT和CCDC算法以煤炭基地尺度为场景进行植被损毁的研究与应用,未曾定量对比分析这两种算法的适用性,导致相关工作者在选择煤炭基地尺度的植被损毁识别方法时主观性较强,缺少科学的数据支撑。

为了解决上述问题,本研究基于GEE云计算平台,以Landsat长时序卫星影像数据为基础,分别采用LT和CCDC算法对神东煤炭基地进行植被损毁识别,获取了1990—2020年植被损毁识别结果,从植被损毁识别、不同地表区域适用性、植被损毁时间3个角度,对比分析了两种算法识别神东煤炭基地植被损毁的适用性,为两种算法在煤炭基地尺度的进一步应用提供了科学数据参考。

1 研究区与数据

1.1 研究区概况

神东煤炭基地位于晋陕蒙三省交界处,地理范围为东经109°41′～111°36′,北纬38°42′～40°06′,如图1所示。自1985年来,神东煤炭开采行业逐渐发展,“十二五”期间,随着国家的重点建设,形成了千万吨矿井群的生产格局,是目前我国最大的煤炭生产基地和重要的优质动力煤出口基地,承担了国家能源保供的重大责任[20]。然而,该地区位于干旱半干旱地带,属于大陆性季风气候,年均降雨量仅有400 mm,水资源稀少[21]。同时,研究区内晋陕蒙各地城镇建设发展迅速,例如研究区西北部的鄂尔多斯市东胜区,自1990年以来,建成区面积呈明显扩张趋势,1990年城市建成区面积为15.11 km2,2008年达到57.79 km2,2020年达到78 km2,其中2000—2008年城市扩张速率最高[22]。在自然本底条件和人类城镇建设、煤炭工业活动的共同作用下,神东煤炭基地生态环境复杂多变、较为脆弱[23],履行国家能源保供任务的责任与脆弱的生态环境保护之间的矛盾突出。

1.2 数据源

本研究采用了1990—2020年的Landsat系列卫星地表反射率数据集 (SurfaceReflectance,SR),分别为Landsat 5的Thematic Mapper (TM)、Landsat 7的Enhanced Thematic Mapper Plus (ETM+)、Landsat 8的Operational Land Imager (OLI),其中TM影像2 023景、ETM+影像2 047景、OLI影像943景。所有数据均经过辐射定标、大气校正、去云、拼接、裁剪等处理。

2 技术方法

2.1 技术路线

本研究通过LT和CCDC算法对研究区植被损毁进行长时间序列的监测,提取植被损毁信息。利用Google Earth的高分辨率遥感影像,通过目视解译对神东煤炭基地的植被损毁进行精度验证。从植被损毁识别评价、不同地表区域适用性评价、植被损毁时间提取精度评价3个方面,对两种算法在神东煤炭基地进行植被损毁识别的适用性进行对比分析。技术路线如图2所示。

2.2 LandTrendr算法

由KENNEDY等[24]提出的LT算法是一种时间分割算法,该算法通过对以年为间隔的时间序列数据进行分割、拟合和平滑,生成像元尺度上数据随时间变化的轨迹,并通过轨迹获取植被损毁发生的时间以及植被损毁程度等信息[25]。本研究轨迹构建步骤为:① 去除云或云阴影所引起的异常顶点;② 定义轨迹构造的最大段数;③ 通过去除最小角度的段对初始的时序轨迹进行简化;④ 评估简化的时序轨迹与原始轨迹之间的差异。因此,LT算法需要进行相关参数设置以确保识别质量,本研究LT算法的参数取值见表1。

表1 LT算法参数设置Table 1 Parameters setting of LT algorithm

2.3 CCDC算法

由ZHU等[26]提出的CCDC算法是一种基于Landsat卫星观测数据的像元尺度时间序列分割算法。有别于其他算法,CCDC算法会利用时间序列内所有可利用的观测数据,即通过Fmask算法对云和云阴影进行掩膜处理,当获得新的Landsat卫星观测数据时,CCDC算法通过结合每个像元所有可用的Landsat卫星观测结果来拟合时间序列模型。如果新的连续观测值超过预期范围(3倍均方根误差),则进行断点标记并估计新的时间序列模型,即在断点前后生成两个时间序列段。CCDC算法可以识别出规定时间序列内的所有断点,并获取植被损毁时间等信息。本研究CCDC算法参数设置见表2。

表2 CCDC算法参数设置Table 2 Parameters setting of CCDC algorithm

2.4 基于GEE平台的植被损毁识别

GEE云计算平台集成了海量地理空间数据,具有强大的分析计算和空间可视化能力,是目前最先进的地理大数据分析平台之一[27]。在GEE云计算平台中,利用JavaScript编程语言,能够快速利用LT和CCDC算法实现区域植被损毁识别[28-30]。归一化差值植被指数(Normalized Difference Vegetation Index,NDVI)可描述植被的生长状态,是评定区域植被综合状况的常用指标之一[31]。因此,本研究以Landsat卫星影像为基础,通过NDVI指数表征植被生长状态,基于GEE云计算平台分别使用LT算法和CCDC算法对神东煤炭基地植被进行损毁识别。其中NDVI通过对卫星影像进行如下波段运算得到:

式中,NIR为近红外波段地表反射率;R为红光波段地表反射率。

另外,需要注意的是:① 部分像元不止一次发生植被损毁事件,为便于后续统一进行分析处理,本研究采用像元第一次发生植被损毁的时间;② LT算法会将研究的起始年份(本研究为1990年)识别为植被损毁的开始点,造成错分现象。因此在实际应用中需要将起始年份(1990年)识别到的植被损毁结果进行剔除。

2.5 适用性分析

本研究从植被损毁识别精度、不同地表区域的适用性、植被损毁时间提取精度3个方面,以总体精度、生产者精度、用户精度、错分误差、遗漏误差等作为评价指标,对比分析了LT和CCDC算法在神东煤炭基地植被损毁识别中的适用性。其中,总体精度是被正确分类的像元总和除以总像元数,即混淆矩阵中对角线上的各类别像元数总和除以影像总像元数目。生产者精度是某一类别被正确分类的像元数(混淆矩阵对角线上该类别的数量)除以该类别的地面真实像元总数。用户精度是某一类别被正确分类的像元数(混淆矩阵对角线上该类别的数量)除以被算法识别为该类别的像元总数。错分误差是在未发生植被损毁地区的样本区域中,被错误识别为损毁区域的面积与样本区域总面积的比值。遗漏误差是在发生植被损毁地区的样本区域中,未被识别为损毁区域的面积与样本区域总面积的比值。

(1)植被损毁识别精度评价。本研究针对两种算法识别结果,分别在识别到植被损毁的区域和未识别到植被损毁的区域随机选择了250个像元作为样本。对于每个样本像元,通过Google Earth提供的历史高分辨率影像进行目视判读,判断该像元是否发生了植被损毁,作为精度评价的真值,生成混淆矩阵,计算总体精度、生产者精度和用户精度。

(2)不同地表区域的适用性评价。本研究以Google Earth中神东煤炭基地原始高分辨率影像为基础,通过目视解译勾画了水体、林地、城市扩张用地、露天采场4类典型样本区域(图1)。水体和林地样本区域未发生植被损毁事件,在此区域识别的植被损毁为算法造成的错分所致。城市和露天采场区域地表植被完全破坏,在此区域未识别到的植被损毁为算法造成的遗漏所致。通过样本区域内两种算法识别的植被损毁面积来判断算法的错分误差和遗漏误差。

(3)植被损毁时间精度评价。本研究面向算法识别出的植被损毁事件,每年随机选取50个样本点(共1 500个),通过Google Earth中的历史高分辨率影像进行目视解译,确定植被损毁的实际时间,作为精度评价的真值,生成混淆矩阵,计算总体精度、生产者精度和用户精度。

3 结果与分析

3.1 植被损毁识别结果

基于GEE云计算平台,使用LT算法和CCDC算法对神东煤炭基地进行了植被损毁识别,结果分别如图3、图4所示。图中A、B区域分别为圈划的鄂尔多斯市东胜区的城市建成区和研究区露天采场集中分布区域。

图3 LT算法植被损毁识别结果Fig.3 Identification results of vegetation damage of LT algorithm

图4 CCDC算法植被损毁识别结果Fig.4 Identification results of vegetation damage of CCDC algorithm

对比图3和图4可知:① 总体而言,LT和CCDC算法在研究区西部识别到的植被损毁多于东部地区。LT算法能够识别到的植被损毁明显少于CCDC算法。② 在鄂尔多斯市东胜区的城市建成区(A区域)中,LT算法基本无法识别植被损毁,而CCDC算法的植被损毁识别效果较好,且算法所识别的植被损毁多发生于2000—2008年,这一现象与东胜区城市扩张历史相吻合。③ 在B区域中,LT算法能够识别到的植被损毁比A区域显著增多,但仍有部分露天采场区域的植被损毁无法识别(图3(e)中圈定区域)。相较之下,CCDC算法相同区域(图4(e)中圈定区域)识别到的植被损毁能够更好地吻合露天采场边界。

3.2 植被损毁识别精度评价

利用选择的样本像元,对两种算法的植被损毁识别结果进行了验证,精度分别如表3和表4所示。

表3 LT算法识别结果的混淆矩阵Table 3 Confusion matrix of the identifying results with LT algorithm

由表3、表4可知:LT算法的总体精度为73.6%,CCDC算法的总体精度为84.4%,LT算法比CCDC算法低10.8个百分点。LT算法“损毁像元”的生产者精度和“非损毁像元”识别的用户精度均比较低,分别为66.9%和54%。CCDC算法“非损毁像元”的生产者精度以及“损毁像元”的生产者精度、用户精度都达到了80%以上,“非损毁像元”的用户精度为78.8%,也显著高于LT算法对应的54.0%。通过与已有研究结果进行对比发现,本研究所得总体精度略低于以往的研究结果(LT算法74.4%～93%,CCDC算法85.1%～99.4%)[10-16]。

3.3 不同地表区域适用性评价

两种算法对于研究区植被损毁识别的遗漏误差、错分误差分别见表5、表6。

表5 植被损毁遗漏误差Table 5 Omission errors of vegetation damage

表6 植被损毁错分误差Table 6 Commission errors of vegetation damage

由表5、表6可知:在城市扩张区域,CCDC算法的遗漏误差为21.6%,LT算法的遗漏误差为91.0%,比CCDC算法高69.4个百分点。在露天采场区域,LT算法的遗漏误差为63.0%,CCDC算法的遗漏误差为8.6%,比LT算法低54.4个百分点。由此可见:LT算法仅能识别部分露天采场区域的植被损毁,对于城市扩张所造成的植被损毁基本无法识别。在林地区域,LT和CCDC算法的错分误差分别为0.7%和11.6%,两种算法都能有效地避免错分现象。但在水体区域,LT和CCDC算法的错分误差分别为3.7%和76.0%,CCDC算法难以避免水体区域的错分现象。

为分析LT和CCDC算法在不同地表类型区域造成遗漏的原因,本研究通过GEE云计算平台提取了两种算法拟合的像元NDVI时序曲线,典型示例分别如图5和图6所示。图5所示的像元在2007年成为露天采场,NDVI值出现明显降低,但LT算法拟合的曲线并未在此处进行分段,没有成功识别出该损毁事件。图6所示的像元于2000年左右经历了城市快速扩张,2000年后未再出现图中虚线框中每年7、8月的周期性高值现象,NDVI数值的时序分布明显改变,但该变化未被CCDC算法识别出。总体而言,在干旱半干旱地区发生植被损毁事件后,尽管NDVI值的分布明显改变,但是变化幅度仍然较小,与植被茂盛地区(如以往研究中的森林区域)NDVI的“断崖式”变化存在较大差异。这导致算法在植被稀疏的干旱半干旱地区仅能识别出NDVI整体呈降低的变化趋势,并未实现对断点变化的识别。受此影响,LT算法仅能识别部分露天采场区域的植被损毁,对于城市扩张所造成的植被损毁基本无法识别。CCDC算法通过均方根误差对时间序列模型进行分段,因此该算法在露天采场和城市扩张用地区域的识别效果较好,相较于LT算法更适用于神东煤炭基地的植被损毁识别。

图5 LT算法损毁识别的遗漏误差Fig.5 Omission errors of the damage identification of LT algorithm

图6 CCDC算法损毁识别的遗漏误差Fig.6 Omission errors of the damage identification of CCDC algorithm

总体上,两种算法在神东煤炭基地的识别精度略低于以往研究的原因在于,以往研究均针对森林区域或小区域的单个矿区,而神东煤炭基地覆盖面积大,且位于干旱半干旱地区,植被覆盖度较低、生态环境的地理空间异质性强。一方面,植被覆盖度低容易导致植被损毁时NDVI时序曲线的下降断点不明显,算法识别不出植被损毁现象。另一方面,在煤炭基地尺度的大面积植被损毁自动识别中,LT和CCDC算法自适应性弱,一个算法在整个区域采用一套相同的参数设置,无法针对不同位置异质性强的特点自适应优化调整。总之,两种算法在神东煤炭基地进行大面积植被损毁识别中效果仍然较好,只是精度稍有降低,根据相关原因分析,两种算法仍存在改进的空间。

3.4 植被损毁时间精度评价

对于能够识别到植被损毁的像元,本研究进一步评估了LT和CCDC算法在识别损毁时间上的精度,结果分别如图7和图8所示。

图7 LT算法植被损毁时间的精度验证结果Fig.7 Precision verification results of vegetation damage time based on LT algorithm

图8 CCDC算法植被损毁时间的精度验证结果Fig.8 Precision verification results of vegetation damage time based on CCDC algorithm

两种算法识别植被损毁时间的总体精度分别为89.2%和79%,LT算法总体精度比CCDC算法提高了10.2个百分点。以“十二五”规划的起始时间(2011年)对植被损毁时间序列进行划分,2011—2020年LT算法总体精度由2011之前的87.1%提高至92%,而CCDC算法总体精度由2011之前的74.1%提高至88.8%,总体精度提升幅度较大。

由图7和图8可知:在生产者精度方面,LT算法的生产者精度主要分布于75%～98%,均值为88.7%。CCDC算法的生产者精度变化幅度较大,均值为81.6%。在某些年份中,两种算法的生产者精度较低,例如1999年,LT算法和CCDC算法的生产者精度分别为68.1%和41.2%,均为历史最低值。在用户精度方面,LT算法的用户精度主要分布于80%～96%,均值为89.4%,用户精度低值出现在2010年,为78%。CCDC算法的用户精度均值为79%,用户精度变化幅度较大,低值出现于1992年和2000年,不足60%。

为分析总体精度提高原因,本研究统计了两种算法在通过目视解译不同地表区域(水体、林地、城市扩张用地、露天采场)典型样本区中识别的植被损毁面积与时间信息,如图9所示。

图9 1991—2020年不同地表区域损毁面积统计Fig.9 Statistics of damaged area of different surface areas from 1991 to 2020

由图9可知:2011年以后,露天采场区域植被损毁面积显著提升,这一现象与“十二五”规划实施后神东煤炭基地的蓬勃发展有着密切联系。LT和CCDC算法在露天采场典型样本区域的植被损毁识别遗漏误差优于城市扩张典型样本区域,因此露天采场面积的增加带动了总体精度得到不同程度的提高。

本研究进一步统计了CCDC算法和LT算法识别损毁时间的误差分布,以“误差=实际损毁时间-算法识别植被损毁时间”方法计算了误差并取绝对值,统计结果如图10所示。根据统计信息,在LT算法中,准确识别损毁时间的比例为88.7%,误差绝对值为1 a的比例为7%;对于CCDC算法,准确识别损毁时间的比例为79.1%,误差绝对值为1 a的比例为10.9%。换言之,LT和CCDC算法识别植被损毁时间的误差≤1 a的比例分别为95.7%和90%。上述结果说明,对于两种算法已经识别出的植被损毁区域,其损毁时间的识别精度均较高。

图10 LT和CCDC算法误差分布Fig.10 Error distribution of LT and CCDC algorithms

4 结论与讨论

本研究基于GEE云计算平台,从植被损毁识别精度评价、不同地表区域的适用性评价、植被损毁时间精度评价等3个方面定量对比分析了LT和CCDC算法在神东煤炭基地进行植被损毁识别的适用性。主要得到以下结论:

(1)在植被损毁识别上,LT和CCDC算法总体精度分别为73.6%和84.4%。从总体精度角度上看,识别效果仍然较好,但是在神东煤炭基地尺度的精度略低于以往研究在其他区域的精度。

(2)在不同地表区域适用性分析中,LT和CCDC算法都能较好地避免林地区域的错分误差。LT算法仅能识别露天采场区域的部分植被损毁,遗漏误差较大,且基本无法识别到城市扩张所造成的植被损毁,但是CCDC算法在这两类区域的识别效果较好。在水体区域,LT算法显著优于CCDC算法。

(3)对于植被损毁时间,LT和CCDC算法识别的总体精度分别为89.2%和79%,误差在1 a内的结果分别占95.7%和90%,对于损毁时间识别效果很好。

总体而言,相较于LT算法,CCDC算法更适用于城市扩张明显、水体面积很少的神东煤炭基地植被损毁识别。本研究成果为神东煤炭基地生态环境质量监管提供了数据参考,更为两种算法在其他煤炭基地尺度的进一步应用提供了方法优选借鉴。

LT和CCDC算法不依赖于影像分类、植被物候模型等条件,仅使用长时序Landsat影像即可快速、自动化地实现煤炭基地尺度上大范围长时序的植被损毁识别,获取植被损毁发生的时间、位置等信息。但两种算法都依赖于植被时序轨迹的变化幅度,本身都存在一些局限性,具体来说:① 在干旱半干旱地区,原始植被覆盖度低,部分像元发生植被损毁时两种算法拟合的时序轨迹变化幅度不明显,导致损毁植被识别的遗漏;② 两种算法只能识别植被损毁事件,无法对自然因素引起的波动和人类活动导致的突变进行区分。针对上述问题,在未来研究中,需要面向采矿活动对植被损毁的时序场景特点,研究一种能够准确自适应识别煤炭基地内大范围露天矿群植被损毁事件的新方法。