基于Landsat时间序列的香格里拉市土地覆盖变化检测

陈培高 张加龙 许冬凡 熊登亮

（1. 西南林业大学林学院，云南 昆明 650233；2. 云南省基础测绘技术中心，云南 昆明 650034）

土地利用和土地覆盖作为地球表层系统最突出的景观标志，呈现出与地球表面景观相关的时间和空间动态过程[1]。同样，土地利用和土地覆盖对全球气候、全球循环、土地的质量以及生态系统等有着重要的影响[2]。随着现代化和城市化进程加快，尤其是建设用地的扩张以及自然用地的衰退，对生态环境造成极大的影响[3]。而检测一个区域长期的土地覆盖变化情况对于区域的生态环境保护、发展方针的制定具有重要的科学依据。

目前已有多种有关时间变化检测的算法，在众多方法中选择合适的算法是进行变化检测的关键步骤。BFAST算法通过整合迭代来降低噪声对于变化检测的干扰，DeVries 等[4]提出此算法并用于检测森林扰动，此方法的结果虽精度较高，但需要调整的参数过多且不能检测出森林的恢复情况和反复的扰动事件；连续变化检测和分类（CCDC）算法常用于检测林地变化扰动，能准确识别发生变化的林地[5]，但此算法对于影像的质量要求过高，且不适用于年际变化较大的区域；CMFDA算法也常用于森林扰动情况的检测，Zhu[6]利用此算法检测森林扰动情况，得到了精度较高的结果，此方法虽然能准确地确定扰动发生的时间，但缺点在于预处理数据的精度要求过高。上述算法在数据处理中均需要相当长的计算时间。而在Google Earth Engine （GEE）平台上使用LandTrendr算法则具有较高的效率[7]。LandTrendr相比于其他算法有着独特的优势，可以更加全面地检测出渐变和突变的事件，不仅能检测出植被地物的变化，还能检测出非植被地物的变化[6-8]。

香格里拉市地处三江并流区域，是世界自然遗产地，具有丰富的森林资源，然而近30多年来，该地区已成为滇西北土地利用/覆盖变化最大的区域之一，是我国生态环境较为脆弱的地区之一[3,9-10]。目前针对香格里拉市的长时间序列土地覆被变化研究鲜有报道，本研究基于Landsat时间序列数据并结合LandTrendr算法以检测长期以来香格里拉市的土地覆被变化情况，旨在较为准确地得到香格里拉市土地覆盖变化情况，为香格里拉市的自然资源管理提供科学数据支持。

1 研究区概况

香格里拉市是云南省迪庆藏族自治州的下辖市之一，位于云南省的西北部、迪庆藏族自治州的东北部、青藏高原横断山脉腹地，地处北纬26°52′～28°52′，东经99°22′～100°19′，全市总面积为11 613 km2，是云南省129个县（市）中面积最大的县级市。香格里拉市境内有众多高山亚高山草甸、沼泽、高原湖泊、河流、高山、水域、森林等自然生态系统，生境类型丰富，是我国生态环境较为脆弱的地区之一[10]。

2 材料与方法

2.1 数据来源

筛选了1986—2016年条带号分别为131/041、132/040、132/041的Landsat Level 2级别的影像为研究数据，传感器类型包括TM、ETM和OLI，共计1 939景。影像从美国地质调查局（United States Geological Survey, USGS）地球资源观测与资源中心（EROS）下载（https://glovis.usgs.gov/）。

2.2 土地覆盖分类及变化动态度模型

本研究采取支持向量机的分类方法进行土地覆盖分类，将香格里拉市全境分为林地、草地、建设用地、水体、永久性冰川雪地以及未利用地6种用地类型[11]。对于分类结果，使用基于混淆矩阵的方法进行分类精度评估。

土地覆盖变化动态度指数考虑到了研究区不同时间段内土地覆盖类型间的转移，能准确反应整个区域土地覆盖变化的剧烈程度，便于在不同尺度上找出土地覆盖变化的热点区域[12]。

式中：LC 为土地覆盖变化动态度指数， ΔLUi-j为i时刻与j时 刻的土地覆盖面积之差，n为总地类数，T为变化总年数。

2.3 归一化燃烧指数

Kennedy等[7]提出，对于干扰事件敏感性较大的为归一化燃烧指数（NBR），而非归一化植被指数（NDVI）或是缨帽变化的第3个分量。NBR可增强较大范围的火灾区域，能够突出植被变化前后的情况，可以提高因植被数量变化而造成的裸露地表的提取精度。

式中： NIR 为近红外波段的反射率值， SWIR2为短波红外波段的反射率值。

2.4 LandTrendr变化检测方法与精度验证

采用LandTrendr时间分割算法对影像进行处理，通过算法能得出干扰发生时的年份、干扰量、干扰持续时间、干扰发生前的光谱值以及恢复相关的信息，核心算法主要包括3个部分：轨迹分割、简化迭代模型以及去除噪声[6-7,13]。本研究基于GEE平台，利用NBR作为LandTrendr时间分割算法的光谱数值，对香格里拉市进行土地覆盖变化检测。主要的控制参数设置标准见表1。

表 1 LandTrendr 时间分割算法主要控制参数Table 1 Main control parameters of LandTrendr algorithm

采用对节点验证的方式对LandTrendr时间分割算法的结果进行验证。随机选取431个训练样本，用于验证LandTrendr算法的精度，并借助Google Earth历史影像，在土地覆盖变化发生的突变时间前后查看Google Earth历史影像，目视判别是否土地发生变化。若土地发生干扰现象，在突变时间前后，通过Google Earth历史影像验证，得出土地由林地变成裸地或者建筑用地等非林地类型，则算法精度准确。

以北纬27.784 485°，东经100.070 757°这块训练样本为例，原影像中植被在1986—2007年处于自然生长状态，2008年受到干扰，从2008—2017年植被又缓慢增长；根据图1观察其NBR变化，蓝色折线为样点的NBR值变化轨迹，红色折线为LandTrendr拟合后的NBR的变化轨迹。

图 1 样点1986—2016年NBR数值Fig. 1 Sample point NBR value from 1986 to 2016

图 2 1986—2016年香格里拉市林地干扰和恢复时空格局Fig. 2 The temporal and spatial pattern of disturbance and recovery in Shangri-La during 1986-2016

2.5 干扰与恢复等级划分

整理香格里拉市的相关资料，得出研究区域内的干扰因素主要有人为砍伐、森林火灾、森林病虫害、地质灾害等[14-16]。干扰、恢复、稳定代表着土地变化的状态，只有干扰和恢复这2种类型能根本的改变土地覆盖类型，即发生干扰和恢复的区域是土地覆盖发生变化的地方。因此，将香格里拉市的土地覆盖变化按干扰的严重程度分为严重干扰、中度干扰和轻微干扰3个级别；与之相反，按恢复程度可以分为重度恢复、中度恢复以及轻微恢复3个等级。各级别的分类与所对应的dNBR阈值见表2。

表 2 干扰与恢复等级划分Table 2 Definition of different levels of disturbances and recoveries

3 结果与分析

3.1 土地覆盖变化分析

使用基于混淆矩阵精度验证的方法对分类结果进行精度验证，得到1986年的总体分类精度为86.55%，Kappa系数为0.836；2016年的总体分类精度为86.81%，Kappa系数为0.836。分类精度满足后续研究的要求。各类型的土地覆盖面积及变化动态度见表3。

表 3 土地覆盖类型面积汇总表Table 3 Summary of the area of each land cover type

1986年，林地面积最大，占比达到87.23%，永久性冰川雪地、未利用地与草地的面积相近，3类总占比为12%，而水体与建设用地面积较少。2016年内，林地面积仍然最大，占比为84.56%，而永久性冰川雪地、未利用地的面积减少较多，草地、建设用地与水体面积得以增长。

1986—2016年，面积减少的地类为林地、未利用地与永久性冰川雪地，面积增加的地类为草地、建设用地与水体；未利用地多数转变为草地，而永久性冰川雪地不断消融导致水体面积持续增加，建设用地以较稳定速度持续增长。

1986—2016年，林地为香格里拉市最主要的土地覆盖类型，虽面积逐年减少，但仍然是香格里拉最主要的土地覆盖类型，变化动态度在各用地类型中最低；而未利用地与永久性冰川雪地是变化动态度最大的土地覆盖类型，此两类土地覆盖类型的减少相对较为剧烈；草地、水体与建设用地以较为稳定的速率逐年增长。

3.2 变化检测结果精度验证

基于TimeSync算法，利用Google Earth历史影像对1986—2016年LandTrendr 时间分割算法的结果进行验证，对431个随机的训练样本进行检验，验证精度见表4。

表 4 基于随机样点目视判读的验证精度Table 4 Accuracy assessment result based on random sample points

训练样本中有166块土地发生干扰，144块土地发生恢复以及121块土地持续稳定状况。其中，经过目视验证得到发生干扰的数量是145个，发生恢复的数量是121个，保持稳定的数量是103个。变化检测总体精度为85.62%，Kappa系数为0.78。在生产者精度方面，干扰生产者精度为87.35%，恢复生产者精度为84.03%，稳定生产者精度为85.12%；在用户精度方面，干扰用户精度为87.88%，恢复用户精度为87.68%，稳定用户精度为80.47%。

3.3 干扰与恢复情况分析

经过LandTrendr时间分割算法后的变化像元数量及比例见表5，可见研究区在31年内的林地变化以中度干扰和轻微恢复为主，整体变化进程较为平缓。

表 5 干扰和恢复数量及比例Table 5 The number and proportion of disturbances and recoveries

在干扰层面上，严重干扰占据的比例很小，占总干扰数的6.93%，而大部分的干扰类型为中度干扰，占总干扰数的81.22%，轻微干扰类型占总干扰数的11.85%。这表明林地的减少主要体现为自然灾害或者人为择伐使得茂密林分变为疏林分，林分的自扰衰退也占据一定比例。

在林地恢复层面，轻微恢复类型占比95.95%，中度恢复类型占4.05%，无重度恢复类型。这表明林地的恢复主要体现为林分的自然生长，而人为的大面积植树造林现象没有发生。

按持续时间的不同，统计香格里拉市恢复、干扰和无变化这3种类型的面积及占比，见表6。

表 6 按持续时间划分香格里拉市干扰和恢复类型统计Table 6 Statistics of interference and recovery type by duration in Shangri-La

在1986—2016年，香格里拉市的土地覆盖未发生变化的占比为81.62%，发生持续恢复的占比为9.28%，发生持续干扰的占比为9.10%，这表明大部分林地未发生变化，表现为恢复或者干扰的面积相近。

按持续时间来看恢复类型，恢复持续21～31 a为主要的恢复类型，这表明大部分的林地恢复体现为林分的自然增长。按持续时间来看干扰类型，持续1～5 a的干扰类型面积最大，持续21～31 a的也较多，这表明香格里拉的林地减少主要体现为自然灾害、人为砍伐，而林分的自扰衰退也占据一定比例。

1986—2016年香格里拉市林地干扰和恢复时空格局见图2。由图2可知，1～5 a的干扰类型和21～31 a的恢复类型是香格里拉市土地覆盖的主要变化类型。1～5 a的干扰类型在香格里拉市周边区域分布较多，而这些区域多为河谷，农业活动更为频繁，故而发生较多的干扰；21～31 a的干扰类型主要分布于香格里拉市西南部，这一区域牧场分布较多，故而产生了长期的干扰；21～31 a的恢复类型在香格里拉市全境均有广泛的分布，这表明香格里拉市多数林地资源都得到了较好的保护，整体呈现自然增长的趋势。

4 结论与讨论

本研究基于Landsat时间序列数据，提取了NBR为参数，使用LandTrendr算法进行变化检测，使用TimeSync方法进行精度验证，所得结论如下：

1）1986—2016年，林地、未利用地、永久性冰川雪地面积持续减少，草地、水体和建设用地面积持续增加；林地为香格里拉的主要土地覆盖类型，多年来面积虽不断减少，但面积占比始终达到80%以上。

2）在6类土地覆盖类型中，未利用地与永久性冰川雪地的变化动态度最高，分别为-7.19%和-4.76%，林地的变化动态度最低，为-0.11%，其余地类变化动态度较为相近；香格里拉土地覆盖的变化动态度总体较低，土地覆盖变化过程较为平缓。

3）香格里拉市的林地变化进程较为平缓，大部分林地未发生变化，在31 a内的林地变化以中度干扰和轻微恢复为主，主要的干扰类型为中度干扰，表现为自然灾害或者人为择伐使得茂密林分变为疏林分；主要的恢复类型为轻微恢复，表现为林分的自然增长。

4）香格里拉市林地总面积多年来持续减少，按干扰和恢复的持续时间来看，1986—2016年间林地发生变化的主要类型为持续1～5 a的干扰和持续21～31 a的恢复；自然灾害与人为采伐的消耗速度远大于林分的自然增长速度，如果不健全林地保护的措施，林地资源将日益稀少。

在方法上，基于Google Earth Engine平台使用LandTrendr算法的确省去了大量的数据处理时间，相较其他模型，能有效地基于长时间序列进行林地扰动分析[13,17-20]。实验的总体精度为85.615%，Kappa系数为0.783，保持了良好的精度，研究结果较为可信。

采用LandTrendr时间分割算法，并利用NBR作为算法参数，检测土地覆盖的变化，对于发生扰动后的恢复情况探测更为敏感，尤其是植被覆盖区域的变化检测效果较好，为森林资源的年度数据更新提供了林地的变化量以及地理坐标位置，以及为恢复区域的植被提供数据的支持，同样给外业开展林地调查提供一定的参考依据。

虽然NBR对于植被的检测比较友好，但其对于建筑物的检测存在不足，可以进一步对比其他遥感因子或者其他有针对性的遥感因子来进行土地覆盖变化检测，比如可以用归一化建筑指数（NDBI）来单独检测香格里拉城市部分，细调LandTrendr算法控制参数的阈值，以提高中纬高寒地区的变化检测精度。