Dynamic monitoring of loess terraces based on Google Earth Engine and machine learning

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

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

Figure 1.  Flowchart of dynamic monitoring of loess terraces

使用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景影像。

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

地类仅分为梯田和其他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年厘米级高清遥感影像人工勾绘以及实地验证。

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

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

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

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

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

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

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

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

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

Figure 2.  Classification probability of the original results and the results of using LandTrendr algorithm at 3 different positions

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

Figure 3.  Variations of annual terraces area and annual mean fractional vegetation cover in the research area from 1988−2019

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

Figure 4.  Distribution of time to use terraces in the research area from 1988−2019

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

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

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