基于Google Earth Engine与机器学习的黄土梯田动态监测

李万源，田 佳，马 琴，金学娟，杨泽康，杨鹏辉

（宁夏大学 农学院，宁夏 银川 750021）

宁夏回族自治区固原市位于黄土高塬沟壑区[1]。长期的过度放牧、不合理耕作，导致该地区植被稀疏、水土流失加剧[2]，严重影响了当地社会经济发展和生态安全。梯田有效缓解了农业生产带来的水土流失问题[3]，从20世纪80年代开始，固原市实施了大面积的坡改梯工程[4]。加之2000年开始实施的国家退耕还林还草工程[5]，该地区的水土流失问题有所缓解，生态环境持续向好[6]。随着遥感技术的快速发展，如何从遥感影像中高效、准确、大尺度地获取梯田时空分布信息，对于指导农业生产、水土保持监测和防治水土流失具有重要的意义。传统的梯田遥感识别主要采用目视解译[7]，该方法精度较高，但存在耗时耗力、成本高、方法复用性差等问题，目前更多用来采集机器学习(machine learning)的样本[8]。近年来，大部分学者采用面向对象或基于像元的监督识别技术，利用决策树(CART)、随机森林(RF)、支持向量机(SVM)、深度学习(DL)等[9−11]机器学习算法，先学习采集的样本，然后利用学习好的模型识别新的样本。面向对象技术较基于像元识别技术，不仅依靠地物的光谱特征，还利用像元和像元之间的关系提高识别精度，识别过程更加复杂，影像分辨率要求更高[7]。但是，无论采用哪种方法进行梯田遥感识别，基本上都是基于单机处理，普遍存在遥感数据获取困难、预处理复杂、性能限制等问题[9]，难以开展大尺度的遥感识别研究。为了解决这些问题，Google公司借助其强大的计算资源与海量数据存储，推出了遥感云平台Google Earth Engine(GEE)[12]。借助该平台，研究人员可以极大扩展自身原有研究的覆盖范围，提供国家乃至全球尺度的研究成果[13]。目前，GEE在大尺度森林变化监测、土地利用类型分类、人类居住地动态监测等[14−16]方面应用广泛，但大尺度梯田遥感识别未见相关报道。为此，本研究在GEE平台支持下，利用Landsat时间序列数据和SRTM数字高程模型(digital elevation model，DEM)，建立每年时间序列影像的百分位数特征。对比3种机器学习算法的分类精度大小，选择分类精度最高的识别结果，应用LandTrendr时序算法逐像元拟合修正时间序列，实现固原市1988−2019年度梯田动态监测的目的。研究结果可为黄土丘陵地区梯田的高效、准确识别和水土保持监测、评价提供参考。

1 研究区概况

固原市 (35°14′～36°31′N，105°19′～106°57′E)位于宁夏回族自治区南部的六盘山地区，辖原州区、西吉县、隆德县、泾源县、彭阳县，国土面积1.05 万km2。属大陆暖温带半干旱气候，年均气温6.3 ℃，年均降水量493.5 mm，降水量由东南向西北递减，年均蒸发量1 472.9 mm，年均无霜期152.0 d。域内地形南高北低，沟壑纵横，黄土丘陵面积达67.9%。地带性土壤以黑垆土为主，但严重的土壤侵蚀导致土壤母质层出露，黄绵土广布。植被总体上由东南的半湿润森林草原区向西北的干旱半干旱草原区过渡[4]。

2 数据源与研究方法

黄土梯田动态监测的流程可分为4个主要功能模块：遥感数据加载、数据预处理、分类算法优选、序列优化。各模块从上到下，层层递进，最终实现黄土梯田动态监测(图1)。

图1 黄土梯田动态监测流程图Figure 1 Flowchart of dynamic monitoring of loess terraces

2.1 数据源

2.1.1 Landsat影像 使用 T1级别 (质量最高)的 Landsat地表反射率数据 (surface reflectance, SR)。该数据产品已经过几何校正、辐射校正和大气校正，空间分辨率30 m，时间分辨率16 d。由于Landsat 5/7/8卫星的服务年限不同，1988−2011年使用Landsat 5影像，2012年使用 Landsat 7影像，2013−2019年使用Landsat 8影像，共使用1 690景影像。

2.1.2 高程数据 采用 30 m 空间分辨率的数字高程模型，具体编号为SRTMGL1_003。

2.1.3 样本数据 地类仅分为梯田和其他 2类。通过Google Earth Pro提供的高清历史影像，利用目视解译法采集样本数据。样本数据包括样点数据和斑块数据。样点数据按时间分为2010−2014年地类属性相同和2000年的样点，以满足Landsat 5/7/8不同卫星分别进行机器学习样本训练的需求。样点采集遵循以下原则：①在研究区生成5 km方形格网，以使样点分布均匀；②保持样点100 m以内属性相同。样点数据共2 673个，梯田样点1 040个，其他样点1 633个。斑块数据为6个随机分布的5 km×5 km 正方形区域，参考 Google Earth Pro中 2019年厘米级高清遥感影像人工勾绘以及实地验证。

2.2 研究方法

2.2.1 合成影像 选择 Landsat对应卫星影像的红波段 (Br)、绿波段 (Bg)、蓝波段(Bb)、近红外 (Bnir)、短波红外1(Bswir1)、短波红外2(Bswir2)6个光谱波段；再经裁边(坏像元)、光谱指数计算(计算方法如表1)、去云后，针对黄土梯田全年季相变化特点[17]，统计每年度内时间序列影像百分位数特征融合影像[18]，即逐像元对某一波段1 a内所有观测值取其10%、25%、50%、75%、90%百分位数，获得该像元位置该波段对应的5个指标波段；再与6个地形特征波段组合，即由数字高程计算得到的海拔、坡度、坡向，以及 3个 3×3、7×7、11×11像元窗口内地形起伏度波段。共计61个特征波段。

表1 光谱指数计算方法Table 1 Calculation methods of spectral index

2.2.2 机器学习 3种机器学习算法为随机森林、决策树、支持向量机，GEE均有内建，可直接调用。另外，针对不同卫星分别进行机器学习，把样点数据分年度映射到对应合成影像并汇总(如Landsat 5包括2000、2010和2011年的样本)，再按9∶1划分样本，90%的样本用于分类器训练，10%的样本用于精度验证。

2.2.3 LandTrendr算法 LandTrendr算法将以年时间序列的值进行分割、逐段拟合、平滑[19]，获取单个像元在整个研究时间段内的整体变化特征。具体介绍参考文献[19]。

2.2.4 识别结果优化 应用前文分类精度最高的机器学习算法，对研究区1988−2019年逐年进行梯田遥感识别。为减少极端天气和人类活动导致识别错误，利用地类在时间序列上连续、稳定的特征，使用LandTrendr算法[19]对识别结果的时间序列(概率为0～1的浮点)拟合平滑处理。参考中国水土保持措施分类[20]，提取坡度＞2°和坡度＜25°区域的梯田，以减少沟壑地及塬地的误分。

2.2.5 精度验证方法 采用混淆矩阵的方法，以总体精度、Kappa系数、生产者精度和用户精度等指标作为识别精度评价依据。具体计算方法参考文献[18]。

2.3 植被覆盖度计算

植被覆盖度(fractional vegetation cover, FVC)采用归一化植被指数和像元二分模型计算。具体计算方法参考文献[21]。

3 结果与分析

3.1 不同机器学习算法的精度

表2为随机抽取的1 051个样点的验证结果。4种精度指标均为随机森林算法最高，决策树算法次之，支持向量机算法最小。随机森林算法基于样点检验的精度分别为：梯田的生产者精度94.46%、梯田的用户精度89.03%、总体精度94.10%、Kappa系数为0.87，都远大于另外2种算法。因此，后文采用随机森林机器学习算法进行梯田遥感识别。

表2 不同机器学习算法识别结果的样点验证精度Table 2 Sample points verification accuracy of the results of different machine learning algorithms

3.2 识别结果与人工勾绘斑块的验证

表3显示：去除交界100 m缓冲区后的验证精度高于未去除时(0 m)的验证精度。另外，经LandTrendr处理后梯田的生产者精度、梯田的用户精度、总体精度和Kappa系数分别为：81.75%、85.97%、93.33%、0.80，均大于LandTrendr处理前的验证精度。

3.3 LandTrendr拟合效果示例

选择3个不同位置来展示LandTrendr算法拟合效果(图2)，位置A原始识别结果在1994、2002、2004年被错误识别为其他类型，位置B原始识别结果在1997年被错误识别为其他类型，在2015年被错误识别为梯田类型。经LandTrendr算法处理后，这些错误类型均被校正。位置C原始识别结果与经LandTrendr算法处理后的结果均为其他类型，识别类型没有变化。

图2 3个不同位置的原始识别结果及使用LandTrendr算法处理后的概率Figure 2 Classification probability of the original results and the results of using LandTrendr algorithm at 3 different positions

3.4 近30年梯田面积和植被覆盖度

经LandTrendr算法处理后的研究区梯田面积(图 3)变化趋势更稳定，从 1988年 5 816.59 km2减少到 2019年 3 146.72 km2，年均减少 90.85 km2·a−1。1988−2019年，研究区植被覆盖度则呈现不断增加的趋势，与梯田面积变化趋势相反。另外，处理前、处理后的梯田面积与植被覆盖度极显著(P＜0.001)相关，其相关系数分别为−0.50和−0.75。

图3 1988−2019年研究区梯田面积与植被覆盖度变化Figure 3 Variations of annual terraces area and annual mean fractional vegetation cover in the research area from 1988−2019

3.5 梯田使用时间分析

图4显示了研究区1988−2019年梯田使用时间长短的分布。从整体上来看，梯田主要分布在六盘山山脉两侧，且西部的梯田使用时间较东部更长。从局部来看，南部的泾源县区域，梯田零星分布，使用时间相对较短；西部西吉县的沟谷条带、中部的六盘山山脉、北部原州区清水河的河谷冲积平原(红色部分)能明显区分出来。

图4 1988−2019年研究区梯田使用时间分布示意图Figure 4 Distribution of time to use terraces in the research area from 1988−2019

4 讨论

4.1 机器学习识别精度的影响与优化

已有的梯田遥感监测研究[3−4]受限于单机处理性能和准确的历史样本采集，其研究内容往往时间短、区域小，限制了长时间序列、大尺度遥感监测的应用与发展。本研究使用模型迁移法，针对每一个传感器独立训练机器学习分类模型，减少了样本采集的难度，得以实现黄土梯田动态监测。然而，机器学习的识别精度主要受样本量、特征、机器学习算法的影响[7]。本研究利用多年采样法增加样本量，选取最优机器学习算法，得到较高的识别精度。另外，关于特征选取，我们前期使用了最大值、最小值、众数、中位数、平均数等多种特征融合方法，但识别精度均低于本研究的百分位数特征融合。而对于深度学习，我们在本地电脑使用相同样本集，多次构建深度学习模型并训练，然而识别精度也低于本研究的随机森林。最后引入LandTrendr算法逐像元拟合时间序列轨迹，有效校正了如图2中的异常值，提高了识别精度。而且，在斑块验证数据与样点采集时同样保留100 m空间误差时，消除2种利用类型相邻区域地理配准误差带来的系统错误后，基于2019年斑块检验总体精度93.33%，与样点验证总体精度94.10%相当，说明训练好的模型随时间迁移应用，识别性能不会降低。

4.2 梯田面积变化分析

整体来看，研究区1988−2019年梯田面积呈减少趋势，植被覆盖度则呈现逐步增长趋势，梯田面积与植被覆盖度极显著相关(P＜0.01)，说明梯田面积减少有助于生态环境向好发展。局部来看，研究区在1988−1996年梯田面积年均减少69.02 km2·a−1，与该时期宁南山区逐步退耕还林还草时间一致；1997−2000年梯田面积年均增长 91.60 km2·a−1，与该时期耕地面积增长趋势相同[22]；在 2001−2005年梯田面积下降较快，梯田面积年均减少 250.51 km2·a−1，远高于 1988−2019年年均减少速率 90.85 km2·a−1，且2001−2005年植被覆盖度年均增长速率是1988−2019年植被覆盖度年均增长速度的4倍，与宁夏“退耕还林工程”生态政策大力实施的时间节点相符；从2007年开始，为巩固退耕还林工作，持续推进生态文明建设，研究区梯田面积下降减缓[6]。另外，研究区西部的梯田使用时间较东部更长，这可能与东部年降水量达650 mm，而西部年降雨量不到450 mm[4]，在东部进行梯田退耕后有利于提高植被成活率有关。

5 结论

基于GEE云平台，使用随机森林机器学习算法与LandTrendr算法，可以高效、准确地实现长时间序列、大尺度的黄土梯田动态监测。相比1988年，研究区2019年梯田面积减少45.90%，植被覆盖度增长52.44%，说明近30 a梯田农业比例逐渐降低，生态环境持续向好发展。