基于Google Earth Engine与多源遥感数据的 土地覆盖变化监测

熊元康，张鸿辉*，梁宇哲，罗伟玲，洪 良

（1.广东国地规划科技股份有限公司，广东 广州 510650；2.广东省土地调查规划院，广东 广州 510075）

土地利用/覆被变化（LUCC）记录了人类在地球表面的空间格局活动，是导致生物多样性减少、气候变化、生态环境演变、生物化学循环乃至全球变化的主要因素[1-2]。快速、准确、全方位地获取一个地区的土地覆盖变化信息，可为该地区的社会经济发展、生态环境建设、国土空间规划等提供重要的基础数据。近年来，遥感技术因其具有监测范围广、实时性强以及使用成本低等优点已被广泛应用于土地覆盖变化监测中[3]。

目前，根据变化监测对象的不同，基于遥感影像的土地覆盖变化监测方法可分为3大类：①基于像素级的土地覆盖变化监测方法，直接利用像元的光谱特征值，通过多时相的遥感数据进行影像间的差值、比值等代数运算，进而获取差异影像，再利用经验或通过样本集获得的阈值进行土地覆盖变化监测，如于冰[4]等通过对GlobeLand30数据进行重编码等操作，提出了一种基于像元转换的土地覆盖变化监测方法；②基 于特征级的土地覆盖变化监测方法，通过不同地类所表现出来的光谱特征、纹理特征等进行土地覆盖变化监测，与基于像素级的土地覆盖变化监测方法相比，该方法不易受遥感影像时相变化的影响，差异信息更突出、监测精度更高，是目前土地覆盖变化监测的主流方法之一，如Jose D A[5]等结合Landsat数据以及L波段SAR数据计算得到各类指数（NDVI、EVI等），再利用随机森林分类器成功提取了缅甸南部Tanintharyi地区的土地覆盖动态变化信息；③基于对象级的土地覆盖变化监测方法，主要针对高分辨率的遥感影像，利用多尺度分割生成不同的基元对象，再结合基元对象的特征进行土地覆盖变化监测，与前 两类方法相比，该方法能获取更加丰富的特征信息，便于提高分类精度，如吴田军[6]等基于SPOT 4的多时相数据，利用对象级的土地覆盖变化监测方法成功提取了广东省东莞市东北区域2005－2008年的土地覆盖变化信息。

目前，利用中高分辨率遥感影像在我国南方地区进行大范围土地覆盖变化动态监测的研究较少（四季多雨，无云数据更加难以获取），且多集中于利用统计数据与土地利用数据相结合的方法进行土地覆盖变化动态监测[7]。为了解决大范围土地覆盖变化监测所面临的低分辨率遥感影像混合像元严重、中高分辨率遥感影像监测范围小和重返周期长等问题，本文以广东省为研究区，基于Google Earth Engine （GEE）云平台，结合中等分辨率的Sentinel-1/2和 Landsat 7/8等多源遥感数据，通过构建多元时间序列影像的方式进行土地覆盖变化动态监测。GEE云平台是一个能在大尺度范围下进行遥感数据处理和分析的云平台，提供了一个完备的集成环境。目前，国内外多个研究组已基于该平台开展了各种对地观测研究，如水稻遥感制图[8]、耕地遥感提取[9]、农作物种植结构 提取[10]等。

1 研究区概况与数据来源

1.1 研究区概况

广东省地处我国大陆最南部，东邻福建，北接江西、湖南，南临南海，西连广西；属于热带或亚热带季风气候区。广东省以广州市为中心，东起潮州市、西至雷州市，横跨21个地级市（图1）；陆地面积为17.97万km2；地势表现为北高南低；地貌类型复杂，北部、东北部和西部均有较高山脉，中部和南部沿海地区多为低丘、台地或平原[11]。随着广东省各地区城市化和工业化的快速推进，耕地和林地不断被蚕食，土地可垦率降低，土地后备资源严重不足，导致人地矛盾日益突出[12]。因此，快速、准确、全方位地获取广东省的土地覆盖变化信息，可为其社会经济发展、生态环境建设、国土空间规划等提供重要的基础数据。

图1 研究区位置与验证样本点的空间分布 （审图号：GS(2016)2929）

1.2 数据来源

1.2.1 遥感数据

由于广东省四季多雨，无云数据难以获取，且中高分辨率遥感影像监测范围小、重返周期长，采用单一数据源的中高分辨率影像不能满足大范围土地覆盖变化监测的需求，因此本文结合中高分辨率的 Sentinel-1/2和Landsat7/8等多源遥感数据进行研究区内的土地覆盖变化监测。遥感影像数量如表1所示。

表1 本文采用的遥感影像数量

1.2.2 地面参考数据

本文结合野外实地考察、航拍影像（0.2 m）、Google高清影像、Planet高清影像（5 m）以及土地利用数据，将研究区内的土地覆盖类型分为6大类，如表2所示，并根据上述数据集获得训练数据集。

表2 土地覆盖类型分类表

1.3 数据预处理

本文利用Sentinel-2、Landsat 8、Landsat 7自身的质量评价波段（QA）进行相应的去云处理。同时，由于不同传感器在设计上的差异，导致不同传感器具有不同的光谱响应函数（SRF），如Sentinel-2和Landsat 7数据在红光波段的均方根误差（RMSE）超过了8%[13]，因此为了减小不同光谱响应函数带来的误差，本文以Landsat 8数据为基础，对Sentinel-2和Landsat 7数据相应波段进行相关的线性转换[14-15]，并将所有数据的分辨率重采样至30 m。

2 研究方法

2.1 多元时序特征选择

土地覆盖变化监测与分析主要的依据是遥感影像时间序列的变化特征信息，因此在构建时间序列时，时序特征的选择显得尤为重要。通过对研究区内土地覆盖变化类型进行分析发现，研究区内主要的土地覆盖变化过程为：耕地或林地转化为建设用地、耕地转化为鱼塘或鱼塘转化为耕地、耕地转化为园林地等。鉴于此，本文选择的时序特征为：

1）归一化植被指数（NDVI）。NDVI常被作为特征参数来评估地表植被的生长状况[16]，能反映红光波段（植物吸收强烈）与近红外波段（植被反射强烈）之间的关系，因此能良好地区分植被与建筑用地等。其计算公式为[17]：

式中，ρNIR为近红外波段的反射率；ρRed为红光波段的反射率。

2）归一化建筑指数（NDBI）。NDBI常被用于城镇建设用地提取，能反映城镇建设用地在近红外和中红外波段具有高反射率值，而其他地物具有低反射率值的趋势。其计算公式为[18]：

式中，ρMIR为中红外波段的反射率；ρNIR为近红外波段的反射率。

3）改进的归一化差异水体指数（MNDWI）。MNDWI是徐涵秋[19]为了改进NDWI[20]提取水体信息的局限性（影像中有建筑阴影的水体提取效果不佳）而提出的，能明显增强水体与建筑物的反差，减少背景误差，因此能良好地区分水体与建筑等。其计算公式为：

式中，ρGreen为绿光波段的反射率；ρMIR为中红外波段的反射率。

4）SAR影像时序特征。针对研究区内长时间序列无云数据难以获取等问题，结合合成孔径雷达（SAR）不受云雨雾等自然条件影响以及全天候等特性，利用不同地类的不同散射机理[21]来区分水体、植被以及建筑等土地覆盖类型。本文将Sentinel-1的VV极化数据作为研究区土地覆盖变化监测的雷达数据。为了将VV极化数据与NDVI、NDBI以及MNDWI数据更好地结合起来，本文利用式（4）对其进行归一化处理。

式中，maxρVV为研究区内SAR长时序数据上的最大VV极化散射系数；minρVV为研究区内SAR长时序数据上的最小VV极化散射系数；ρVV为研究区内长时序上某一时期的VV极化散射系数。

2.2 分类特征光谱响应分析

根据选取的时序特征，通过Google的高清影像，本文选取研究区内典型的土地覆盖变化类型进行时序光谱响应分析，结果如图2所示，可以看出，在耕地转化为建设用地的过程中，NDVI和VV后向散射系数在不断下降，而MNDWI和NDBI在地类表现为耕地时，随时间的变化出现了较大的浮动，在地类转化为建设用地时，随时间的变化表现较为平稳；在林地转化为建设用地的过程中，NDVI和MNDWI随时间的变化呈下降趋势，而NDBI和VV后向散射系数随时间变化呈逐步上升趋势；在林地转化为裸地的过程中，NDVI、MNDWI和VV后向散射系数随时间变化呈逐步下降趋势，而NDBI呈逐步上升趋势；在水体转化为裸地的过程中，NDVI、NDBI和VV后向散射系数随时间变化呈上升趋势，而MNDWI呈下降 趋势。

图 2 各种土地覆盖类型的时间序列

2.3 土地覆盖类型提取方法

本文分别利用随机森林分类（RF）、支持向量机（SVM）以及回归分类决策树（CART）进行土地覆盖类型提取，其中RF和SVM已被广泛应用于土地覆盖类型分类中，并取得了良好的分类结果，而CART易于表达且能很好地解释某些特定的规则[22]。

RF是一种集成学习方法，采用Bootstrap抽样技术从原始数据集中抽取训练集，并通过训练集构建CART决策树。影响RF分类器性能和效率的主要参数为决策树数量、候选特征子集以及叶节点最小样本数。本文将这3个参数设置为：RF分类器的决策树数量=100；候选特征子集=4；叶节点最小样本数=1。SVM是一种非参数机器学习算法，核心是找到一个最优的超平面作为高维空间中的决策函数，进而将输入向量分成为不同的类别。影响SVM性能和效率的主要参数包括核函数的选择和cost参数。本文选择线性函数作为SVM的核函数，并将cost参数设置为100。CART是一个在二进制递归分区过程中成长的决策树，通过子集内变量的最大方差和最小方差将训练数据集分成不同的类别。树的最大深度参数决定了CART模型的复杂性，大的深度可能具有更高的精度，但也会增加过度拟合的风险。综合以往的研究和实际需求，本文将树的深度参数设置 为10。

2.4 精度评价

对遥感影像分类结果进行精度评价是一个十分重要的工作[10]。为了保证精度验证的准确性，本文选择总体精度（OA）和F-Score作为土地覆盖变化监测结果的精度评价指标。这两个指标均来源于混淆矩阵，其中OA用以评价整体算法的有效性，而F-Score用以评价每一类的分类精度。

3 研究结果与分析

3.1 精度分析

为了对土地覆盖类型变化监测结果进行精度评价，本文通过分层随机采样在研究区内生成了1 000个验证样本点，如图1所示，并通过人工目视解译判读其是否为土地覆盖类型发生变化的区域，进而获得587个非变化样本点和413个变化样本点。由精度评价结果可知（表3），RF方法变化监测的OA为83.5%，F-Score为0.82；SVM方法变化监测的OA为73.5%，F-Score为0.58；CART方法变化监测的OA为70.9%，F-Score为0.51。因此，RF表现最优，SVM次之，而CART结果最不理想。

表3 分类精度统计

3.2 土地覆盖类型提取结果

根据多元分类特征，在GEE云平台上分别融合生成2017年、2018年研究区内的多元分类特征影像（49个特征波段、12个NDVI时间序列波段、12个MNDWI时间序列波段、12个NDBI时间序列波段、12个VV后向散射系数时间序列波段，1个Slope波段），并同时利用收集到的训练数据集进行RF、SVM以及CART训练，以建立研究区内不同土地覆盖类型的监测模型。在此基础上，利用训练好的RF、SVM以及CART分类模型进行2017年、2018年研究区内的土地覆盖类型监测，结果如图3所示。

图3 广东省2017－2018年土地覆盖类型监测结果

3.3 土地覆盖变化提取结果

根据2017年、2018年广东省土地覆盖类型分类结果，本文采用基于空间叠置分析和统计原理分析的分类后变化检测方法进行以a为单位的土地覆盖变化监测，结果如图4、5所示。在此基础上，本文对研究区内的土地覆盖变化区域的空间分布状况进行统计发现，研究区内土地覆盖类型发生变化的区域主要集中在珠三角地区（广州市等9个地区），其次为粤东地区（潮州市等4个地区），最后为粤西地区（湛江市等3个地区）和粤北山区（清远市等5个地区）。由于珠三角地区人口的集聚和工业的发展，其土地覆盖变化类型主要为耕地或林地转化为建设用地或工业用地；由于粤东部分地区人口众多，其土地覆盖变化类型主要为耕地或林地转化为建设用地；由于粤西地区和粤北山区是重要的农业生产空间，其土地覆盖变化类型主要为耕地转化为鱼塘或鱼塘转化为耕地，以及耕地转化为园林地。

图4 广东省2017－2018年土地覆盖变化监测结果

图5 珠江口2017－2018年土地覆盖变化监测结果

为了发掘影响土地覆盖变化的潜在因素，结合2004－2018年广东省的各类统计数据，本文对广东省的土地覆盖变化进行分析发现，影响广东省土地覆盖变化空间分布的主要原因为人口增长、区域社会经济发展以及政策变化，如图6所示。在人口增长方面（图6a），珠三角地区人口年均增长99.2万人，粤东地区人口年均增长4.6万人，粤西地区人口年均增长 7.8万人，粤北山区人口年均增长8.8万人；在区域社会经济发展方面（图6b），珠三角地区的年均生产总值增长为5 031亿元，粤东地区的年均生产总值增长为386亿元，粤西地区的年均生产总值增长为 438亿元，粤北山区的年均生产总值增长为281亿元； 在政策变化方面，主要体现为政府通过制定发布相关的政策来干预、调整土地的使用（如加强基础设施建设等），由研究区内各地区的固定资产投资（基本建设投资、更新改造投资、房地产开发投资）情况可知（图6c），珠三角地区的年均固定资产投资增长为 1 833亿元，粤东地区的年均固定资产投资增长为 397亿元，粤西地区的年均固定资产投资增长为290亿元， 粤北山区的年均固定资产投资增长为242亿元，进而导致研究区内建设用地与耕地面积比值的不断增加，即建设用地面积不断增加，而耕地面积不断减 少（图6d）。

图6 影响广东省土地覆盖变化的潜在因素

4 结 语

本文以广东省为研究区，通过分析研究区内的主要土地覆盖变化类型，选取NDVI、MNDWI、NDBI、VV后向散射系数的时序影像以及DSM影像作为研究区土地覆盖类型监测的特征波段影像，同时结合训练数据集、RF、SVM以及CART对研究区内 2017年、2018年的土地覆盖类型进行了监测，并在此基础上进行了研究区内的土地覆盖变化监测。通过对土地覆盖变化监测结果进行分析，本文得出的主要结 论为：

1）构建的49个分类特征波段以及RF分类器能较好地适应于研究区内的土地覆盖类型监测。较SVM和CART而言，RF具有最高的土地覆盖变化监测精度，其土地覆盖变化监测OA为83.5%，F-Score为 0.82。

2）对研究区内土地覆盖变化监测结果进行分析发现，研究区内土地覆盖发生变化的区域主要集中在珠三角地区，其土地覆盖变化类型主要为耕地或林地转化为建设用地或工业用地；其次为粤东地区，其土地覆盖变化类型主要为耕地或林地转化为建设用地；最后为粤西地区和粤北山区，其土地覆盖变化类型主要为耕地转化为鱼塘或鱼塘转化为耕地，以及耕地转化为园林地。

3）结合统计年鉴等相关数据，对研究区内的土地覆盖变化进行分析发现，影响研究区内土地覆盖变化空间分布的主要原因为人口增长、区域社会经济发展和政策变化。