基于GEE平台的国家公园土地覆盖变化遥感检测方法构建

毛丽君，李海涛，薛晓明，李建伟，李明诗

(1.南京林业大学林学院，江苏 南京 210037；2.南京森林警察学院刑事科学技术学院，野生动植物物证技术国家林业和草原局重点实验室，江苏 南京 210023；3.临沂市城市管理综合服务中心，山东 临沂 276000；4.云南省林业调查规划院大理分院，云南 大理 671000；5.南京林业大学，南方现代林业协同创新中心，江苏 南京 210037)

2017年9月我国《建立国家公园体制总体方案》公布(http://www.gov.cn/zhengce/2017-09/26/content_5227713.htm)，至今已陆续建立11个国家公园体制试点区。其中，钱江源国家公园体制试点区(简称钱江源国家公园)是长三角经济发达地区唯一的试点区，集体林占比高，存在自然保护与社区发展利益相冲突的问题[1-2]。国家公园的设立和管理能否发挥保护成效，直接影响我国自然保护地体系建设的成败[3]，而土地覆盖变化能够直接反映国家公园建立前后的保护与干扰状况，对于评估管理有效性、确定脆弱区域并制定保护战略至关重要[4-5]。目前我国大部分国家公园体制试点区都由原先的保护区及周边连接地带整合而成[1]，重新整合后国家公园的自然生态系统势必受到内部或相邻地区人类活动的影响，因此，高效、可靠地监测土地覆盖动态变化成为国家公园保护与管理的紧迫需求。

Landsat长时间序列遥感影像的开放获取为研究土地覆盖变化新方法提供了基础[6-7]，过去30年发展了许多基于遥感影像的变化检测技术[4,8-9]，其中传统的分类后比较法简单直观，最常用于自然保护地监测[7]，但主要限制因素是缺乏可靠的历史训练样本数据。在森林占主导地位的国家公园中，受制于严格的管理制度及其复杂地形，收集和识别足够的地类样本往往存在挑战[10-11]。目视解译方法利用高分辨率影像收集样本点或验证分类精度，但是通常难以获取可用的高分辨率历史影像。Collect Earth软件提供免费的增强型图像解译和样本收集环境，将公共调查的易用要素与专业结构化目视图像解译任务相结合，实现多种不同空间和时间分辨率图像的地理同步显示，包括DigitalGlobe、Sentinel-2、Landsat和MODIS影像、百度地图和必应地图等[12]。因此，运用基于增强型目视解译的稳定地类样本选取方法能够提升国家公园历史遥感影像分类训练样本的识别能力。而分类后比较法提取变化信息的可靠性取决于每幅图像的分类精度。近年来机器学习分类算法发展迅速，其中随机森林(random forest，RF)算法计算精度高，模型训练时间少，对训练样本数量和质量的敏感度较低[13-14]，已被成功应用于自然保护地分类研究中[7,15]。但是仅基于光谱信息的分类通常无法准确描述土地覆盖变化，使用RF算法结合光谱指数、纹理、高程等特征能够提高分类准确性。此外，季节性合成图像或多时相的密集图层堆栈能够去除单幅图像中的云以及其他缺失数据，并且能有效利用物候信息，从而提升土地覆盖制图精度[11,15-16]，但这在单机应用环境中会给计算资源带来巨大压力。GEE(Google Earth Engine)是一个提供全球尺度地球观测数据存储和用户友好界面的开放数据访问平台[17]，已被广泛应用于遥感研究领域[7]。笔者整合上述系统的优势，通过Landsat数据和GEE平台绘制、监测和量化钱江源国家公园实施生态保护政策以来的土地覆盖变化，以期为国家公园高效、有针对性及低成本的管理提供方法参考。

1 材料与方法

1.1 研究区概况及数据来源

钱江源国家公园位于浙江省衢州市开化县境内(118°03′～118°21′E，29°10′～29°26′N)，总面积约252 km2。空间范围包括古田山国家级自然保护区、钱江源省级风景名胜区及两者之间的连接地带[18]。钱江源国家公园属于亚热带季风性湿润气候，适合多种农作物生长，但平原地貌稀缺，典型植被为低海拔亚热带常绿阔叶林。所在开化县年均降水量为1 963 mm，年均气温16.2 ℃，无霜期252 d，年均日照1 334.1 h[2]。2004年浙江省全面启动森林生态效益补偿制度后，当地生产方式从传统农耕转变为发展森林旅游、林下经济及农家乐等[19]。大量年轻人外出寻求就业机会，早期耕作的农田逐渐被废弃。研究区域如图1所示，底图来自SRTM(shuttle radar topography mission)数字高程数据[20]。

图1 钱江源国家公园空间范围示意Fig.1 Locations of the study area

研究使用的数据包括Landsat Collection 1 Tier 1 Raw数据、SRTM数据、OSM[21](open street map)开源矢量数据、实地调查数据(2019年11月)等。其中Landsat数据来源于GEE(https://developers.google.com)平台；OSM数据下载地址为https://www.openstreetmap.org(下载日期为2019年12月15日)；其他资料来自钱江源国家公园官方网站(http://www.qjynp.gov.cn)。结合遥感数据分辨率及研究区土地覆盖特点，将土地覆盖类型划分为耕地、森林、草地、水体、人造地表和裸地6类。

1.2 研究方法

1.2.1 历史样本点解译

Collect Earth与GEE代码编辑器的集成使用户能够通过自动处理连续多年的卫星图像提供派生数据的概览，作为高分辨率影像解译的补充，从而超越了传统的目视解译过程。设计Collect Earth数据收集表单，联合当地专家使用Google Earth、Planet地图和GEE Apps执行样本点增强型目视解译，直观识别土地覆盖类型的变化。综合使用实地调查数据、Google Earth高分辨率影像、OSM矢量数据、钱江源国家公园官方网站资料选取训练和验证样本。均匀选择样本点以消除样本间的共线性，共选择了830个样本点作为独立样本集。将样本点随机分割为70%和30%两部分，分别用于分类器训练和精度验证。

1.2.2 分类特征输入

首先，在GEE平台筛选研究区2001年、2009年、2017年春季(3月1日—5月31日)和秋季(9月1日—11月30日)所有Landsat可用影像。尽量选择拥有全色波段的卫星传感器产品(即Landsat-7和Landsat-8)，但由于2009年Landsat-7 ETM+传感器在研究区范围内存在条带噪声，故采用Landsat-5数据替代。所有可用Landsat影像基本信息如表1所示。然后使用简单复合算法分别创建春季和秋季的无云Landsat复合数据。高通滤波法具有理想的视觉增强效果，并能较好地保持原始图像的定量信息[22]，运用高通滤波法对Landsat季节性复合数据进行全色锐化处理，提高分类输入数据的空间分辨率。由于2009年Landsat-5影像缺乏全色波段，故用最邻近点法将其重采样为15 m 空间分辨率。

表1 Landsat数据筛选结果Table 1 Data screening results of Landsat

利用GEE平台中预处理后的季节性复合图像派生分类特征，包括所有可用光谱波段、8种光谱指数(归一化植被指数、增强型植被指数、绿色归一化指数、绿叶指数、可见光大气阻抗指数、绿色大气阻抗指数、改进的土壤调节植被指数和改进的归一化水指数，其简称依次为NDVI、EVI、GNDVI、GLI、VARI、GARI、MSAVI和MNDWI，具体公式见表2)、基于灰度共生矩阵的纹理均值，并且纳入SRTM数据派生的地形特征(包括高程、坡度和坡向)作为辅助分类数据。

表2 研究中使用的光谱指数Table 2 The spectral indexes used in the study

1.2.3 图像分类及变化检测

使用RF分类器生成3个时期的土地覆盖分类图，首先使用分类样本点对选取的实验特征组合数据进行采样，用于训练和验证。其次训练样本，设置RF算法中的两个参数，ntree为500，mtry默认为输入特征总数的平方根；然后，根据训练好的数据对选定的特征组合图像进行分类。生成分类图像后，利用众数移动窗口滤波器进行图像平滑，最小化混合像素及边界效应造成的误分类。最后，利用QGIS(Quantum GIS)软件绘制2001—2009年和2009—2017年期间的土地覆盖类型转移矩阵和土地覆盖类型转化图。在分析土地覆盖类型转移矩阵的基础上将转化图像进行分类整合，描述基于Landsat影像得到的主要土地覆盖变化类型，包括造林活动、耕地废弃和开发建设。根据国家公园土地覆盖变化情况制定相应的像元转换规则：早期绘制为耕地、草地或裸地的像元，后期转换为森林像元，则表示造林活动；在早期标记为耕地而后期为草地的像元表示耕地废弃；开发建设则是早期标记为非人造地表类型(耕地、森林、草地和裸地)，后期图像中变为人造地表类型的像元。

2 结果与分析

2.1 样本点解译结果分析

解译过程主要依据2014—2017年Google Earth影像、2009年Planet高分辨率影像，以及GEE中不同时期的Landsat影像及其NDVI指数变化曲线，植被的季节性特征在多光谱假彩色合成图像中最为显著(表3)。

表3 各类土地覆盖样本点数量Table 3 The number of sample points of different land cover types

如表3所示，土地覆盖样本点数量变化较大，其中耕地减少和草地增加趋势明显，2001年裸地样本数量较少。综合多种数据源执行增强型目视解译时，可以直接观察国家公园土地覆盖经历的主要变化过程，包括人造地表增加、水体干涸等。

2.2 土地覆盖分类制图结果分析

通过RF算法分类得到2001年、2009年和2017年钱江源国家公园的土地覆盖分类结果，其精度评估结果见表4。

表4 土地覆盖分类的精度评估结果Table 4 Accuracy assessment results for the classified maps

为了提高土地覆盖变化地图的可靠性，通过人工编辑消除明显的分类误差，尤其是山体阴影导致人造地表和水体的错分类。利用30%的验证样本估计总体分类精度见表4，研究区土地覆盖分类结果见图2。

图2 基于RF算法的土地覆盖分类图Fig.2 Land cover classification maps generated from RF algorithm

如表4所示，研究期间土地覆盖分类结果的总体精度均较高，其中2009年分类图的总体精度最低，可能是由于采用的Landsat-5数据未进行全色锐化处理，导致采样数据中存在混合像元数量多所致。而使用Landsat-8作为输入数据的2017年分类结果总体精度最高。在6种土地覆盖类型中，对森林的提取始终比较准确(表4的用户精度和生产者精度)；分类结果混淆主要发生在耕地、草地、人造地表和裸地之间。从总体分类结果看，耕地类型的用户精度低而生产者精度高，代表被过度分类。相反，人造地表和裸地类型很可能分类不足。由于裸地的训练样本数量过少，因此在2001年未能有效评估。

2.3 土地覆盖变化结果分析

研究期间造林活动、耕地废弃和开发建设变化类型的空间分布分别见图3。由图3A可见，造林活动广泛分布于核心保护区之外区域，证明保护区的设立、扩区、升级以及国家公园体制试点的开展对维持森林生态系统的原真性和完整性起到显著推动作用。

另外，2001—2009年期间造林面积为1 940.55 hm2，是2009—2017年(1 062.49 hm2)的近2倍(表5)。研究期间森林覆盖率持续增加，表明当地在实施森林生态效益补偿制度后开展了大规模造林及森林保护活动。其中，52.79%的造林活动分布于连接地带，该区域包括大横、真子坑、高田坑、西山等多个古村落，这些地域原本长期处于传统农耕活动。钱江源省级风景名胜区范围内造林比例最低，该区域造林活动集中于枫楼坑、仁宗坑等古村落。

图3 2001—2017年研究区造林活动、废弃耕地和开发建设分布图Fig.3 Maps of afforestation activities,cropland abandonment,and development and construction generated from the spatial overlay analysis during 2001-2017

表5 2001—2017年各分区造林活动与耕地废弃面积及占比Table 5 The area and proportion of afforestation activities,cropland abandonment in different zones during 2001-2017

如图3B所示，研究区内耕地普遍发生退化现象，未变化的耕地大多位于居民区附近。表5显示，2001—2009年期间耕地废弃面积达到1 800.36 hm2，2009—2017年这种现象持续发展。耕地废弃现象反映出在实施保护生态环境政策的同时，当地产业结构转型以及劳动力流失现象明显。其中，60%左右的耕地废弃现象发生在连接地带，大片传统农耕区被弃置。钱江源省级风景名胜区范围内耕地弃置比例最低，分别为15.20%(2001—2009年)和17.35%(2009—2017年)。该区域依托齐溪水库和莲花溪沿线景观大力发展生态旅游、农家乐和民宿等产业，因此景区周边耕地得以保留。古田山自然保护区也呈现类似的耕地废弃模式，耕地废弃现象主要发生在靠近保护区核心区域的一侧。

图3C显示了研究期间的开发建设进展。许多建设项目沿河道以线性方式分布，主要分布在钱江源风景名胜区以及连接地带，如西莲线公路建设工程及其沿线景点打造项目、枫岭头隧道建设工程，以及连接地带各个古村落与外界通行道路的建设。齐溪水库东侧出现一些新建建筑，用途主要为民宿和餐饮等旅游配套设施。相对而言，古田山自然保护区范围内的开发建设面积最少。由于受到地形限制，山区的人造地表覆盖范围较小，且大部分道路宽度不足一个Landsat像元大小，因此新建道路呈不连续的破碎状分布。

3 讨 论

本研究采用GEE、Collect Earth、QGIS等开源软件平台，基于开放获取的Landsat、SRTM和OSM数据分析研究区的土地覆盖变化特征。研究证明了利用Landsat历史档案数据对空间范围广阔的国家公园进行长期动态监测的可行性，这种做法在国内经济遭受新冠疫情冲击的背景下显得尤其重要，因为当前各级财政资源有限，科研与行政事务需要充分利用免费开放的数据资源、产品和技术平台。基于RF算法融合不同季节的多时相Landsat图像分类能够提高土地覆盖分类精度，但由于该地区景观异质性较高，地形导致的山体阴影仍然是使人造地表、水体、裸地等类型产生识别混淆的主要原因。与Tsai等[15]对梵净山国家级自然保护区的土地覆盖制图结果相比，本研究中采用Landsat系列数据衍生的各类土地覆盖类型制图精度普遍更高，说明土地覆盖分类图的总体质量可靠，适用于土地覆盖变化检测。当前Landsat-8数据具有广阔的应用前景，Landsat-5数据有助于在使用一致性数据的基础上进行历史回顾分析，然而扫描线校正器故障造成的图像条带限制了Landsat-7数据的使用[31]。虽然Sentinel-2影像时空分辨率、光谱分辨率比Landsat系列更具优势，但数据可用的时间跨度相对较短，目前最适合在景观尺度上进行长时间序列变化检测分析的仍然是Landsat系列数据。

由于国家公园地形复杂以及对核心保护区管理严格，导致无法进入其内部收集足够的样本数据。对于年代久远的遥感影像，根据近期的实地观测数据难以可靠地重建过去的地类样本[5]。Collect Earth提供的增强型目视解译方法能够结合各类高分辨率影像、长期档案数据与专家经验，高效、可靠地选取历史样本点[32]。Collect Earth在许多地点呈现多幅不同采集日期的图像，提供的历史图像足以评估近期发生的土地覆盖变化，同时Collect Earth也在不断扩充其档案，并扩大卫星图像的覆盖范围。使用甚高分辨率影像收集的信息对基于半自动方法的土地覆盖分类研究极为重要，因为半自动分类是通过直接观测到的少量样点信息训练算法，然后将算法应用于相对大量的未直接观测区域的分类过程。Collect Earth是目前唯一可以同时访问所有免费可用影像档案的工具[12]，降低了用于图像采集的成本和时间。免费甚高分辨率影像不可避免地存在时空缺口，而Collect Earth与GEE的集成方便地利用中低分辨率Sentinel-2、Landsat和MODIS图像及其光谱指数的时间序列变化趋势填补了这些空白信息[12]。在半自动土地覆盖分类过程中，集成多个甚高分辨率图像档案可以促进训练样本数据的生成和精度评估，节省购买高分辨率影像所需的高昂潜在支出。

我国已经从不同角度开展国家公园试点区生物多样性监测研究[2]，但是对于重新整合的国家公园，尚未开展针对性的土地覆盖变化监测和评估。为了填补这一空白并据此检验国家公园的管护效果，及时了解土地覆盖变化状态及趋势非常必要。国外已有研究利用分类后比较法识别毁林[33]或农业扩张[10]状况，据此判断国家公园的设立与管理是否有效。近年来我国生态文明建设水平不断提升，为世界森林覆盖率增长作出重要贡献，国家公园中森林和耕地的变化趋势往往与国外研究相反，因此应关注生态系统变化中的造林过程，通过遥感技术识别国家公园建设以及生态修复工程(如退耕还林)实施的有效性[15]。本研究结果表明，随着生态保护力度的加强、产业结构转型和劳动力外流，研究区大量耕地被荒废，退化为茅根草地。事实上，废弃耕地也是一种可以合理利用的资源[34]，在土地稀缺和受到严格保护的国家公园能够有效解决保护和利用的矛盾，因此如何进行废弃耕地的修复和利用是当前的重要课题，比如开展退耕还林或发展旅游相关产业。钱江源国家公园现已开展集体林地地役权改革，针对公园范围内的山林权属归村集体所有的森林、林木和林地，在不改变权属的基础上，将上述自然资源交由钱江源国家公园管委会统一管理，管委会每年支付一定的补偿金[35]，这一举措可为森林生态系统及废弃耕地的生态修复工作奠定良好的基础。本研究结果也能够为国家公园范围内废弃耕地修复项目的规划和实施提供空间明确的分布地图。