基于机器学习的青藏高原天然草地盖度时空变化特征研究

孟新月， 葛 静， 侯蒙京， 冯琦胜， 金哲人， 高金龙， 梁天刚

(草地农业生态系统国家重点实验室, 兰州大学草地农业科技学院， 甘肃 兰州 730020)

植被盖度是指地上植被占地表总面积的百分比，能够很好地反映植被的生长情况，也是生态系统变化的敏感指标[1]。植被盖度的大小能够反映植被的茂密程度以及植被群落生长动态，同时也可以揭示区域生态环境的变化状况，明确草原植被盖度对评价草地资源、监测草地退化状况具有重要作用[2]。

地表实测法和遥感反演法是获取植被盖度的两种广泛使用的方法[3]。地表实测法以样方尺度来测量，按其测量原理可分为目视估测法、采样法及仪器法3类[4]。但是，传统的地表实测法耗时长、所需精力多、精度低，且只适用于较小的空间尺度[4-6]。近年来，“3S”(Remote Sensing,RS;Global Position System,GPS;Geographic Information System,GIS)技术迅猛发展，具有平台广阔、层次丰富、时相较长等特点，在草地资源监测中发挥着日益重要的作用[7]。吴见等[8]依据像元二分模型法收集了黄山市在各个年度的植物盖度数据，并利用Landsat TM遥感技术影像，分析了黄山市植物盖度空间变化特征，结果表明1988—2018年黄山市植被盖度整体上呈现出增加趋势；赵国忱等[9]利用Landsat-8数据和像元二分模型，并利用归一化植被指数(Normalized difference vegetation index,NDVI)估算了辽宁北票市植被盖度，结果表明2013—2015年间北票市东部为主要增长区域，2015—2017年间，南部和北部有较大增长。宋清洁等[10]通过地面实测数据结合中分辨率成像光谱仪(Moderate-resolution imaging spectroradiometer,MODIS)植被指数数据，分析了增强型植被指数(Enhanced vegetation index,EVI)和NDVI与草地盖度的线性关系，筛选出甘南州草地盖度最优反演模型为基于EVI构建的对数模型，其R2为0.707，RMSE为10.69%。但是，已有研究也存在草地盖度反演误差大、反演模型构建方法有待改进等方面的问题。虽然利用MODIS等遥感数据和像元二分模型方法可以快速模拟植被盖度随时间序列的动态变化，但在估算盖度时尚存在一定缺陷，如不能揭示植被覆盖所体现的生态过程，并容易受纯植被像元的参照光谱值的影响等[11]。由于单一的植被指数模型和多因素线性回归模型精度较低、稳定性差、存在局限性[12]，机器学习模型已逐步取代传统算法成为构建草地盖度模型的主流方法。机器学习模型是基于数据驱动的，其自动检索和解释数据的方法比较灵活，可用于任何训练任务，准确性更高。机器学习模型能够估计多个变量间的复杂关系，与传统模型相比更具稳健性，能有效提高模型的预测精度[13-14]。为了改善传统模型在实际应用中存在的不足，一些学者使用机器学习模型在草地盖度遥感反演方面进行了初步探索，如陈黔等[15]使用Landsat数据与其他数据进行北方沙地灌木覆盖度的估算，研究表明利用机器学习算法可以很好的实现这一估算(R2=0.72，RMSE=13.73%)。Ge等[16]基于MODIS卫星数据和野外实测盖度数据，构建了黄河源区4种遥感反演模型，结果表明支持向量机模型(Support vector machine,SVM)为草地盖度反演的最优模型，比单因素模型的R2提高0.08～0.14。然而，机器学习算法需要大量草地盖度观测数据和多种生态环境变量，存在变量筛选、模型调参等多个环节的大量分析，研究结果受诸多因素的影响，因此在天然草地盖度最优反演算法研究方面仍需要进行深入探索。青藏高原海拔高，地貌丰富，面积广大，草地资源异常丰富，由于自然和人为的双重影响，青藏高原一些地区草地资源呈退化趋势。我国已有多位学者针对青藏高原草地盖度变化开展了大量研究。如邵伟等[17]利用统计资料分析表明西藏高原草地退化严重。高清竹等[18]基于遥感数据对藏北地区草地退化进行遥感监测，结果表明藏北地区1982—2004年草地退化较为严重，重度退化草地面积占草地总面积的8%。陆晴等[19]利用NDVI数据和地面实测数据，对青藏高原1982—2013年高寒草地覆盖时空变化进行了研究，结果表明青藏高原高寒草地生长季NDVI表现为从东南到西北逐渐减少的趋势。但前人研究范围大多限于青藏高原部分区域，从整体上系统进行青藏高原盖度研究较少。

基于以上考虑，本研究利用2003—2018年草地植被生长季野外实测数据，首先用最小绝对压缩变量筛选方法(Least absolute shrinkage and selection operator，LASSO)筛选合适的变量，再分别建立单因素回归模型和多因素机器学习模型，基于最优模型进一步分析青藏高原地区2001—2019年生长季草地盖度的时空变化特征。本研究旨在通过分析草地植被盖度的时空变化特征，为青藏高原地区草地资源的长期利用和动态监测提供科学依据[20]。

1 数据及方法

1.1 研究区概况

青藏高原(26°50′～39°19′N，78°25′～103°04′E)包括西藏自治区、青海省和四川省、甘肃省、新疆维吾尔自治区及云南省的部分地区，地域辽阔，平均海拔在4 000 m以上，有“世界屋脊”之称[21]。气候具有辐射强、昼夜温差大等特点，年平均气温范围为-6℃～3℃；降水多集中在6—9月，占全年降水量的70%左右，且东西部降水量差距较大，东南部降水最高可达1 000 mm，而西北部降水多为40～100 mm[22]。青藏高原草地总面积约为121 742 km2，占高原总面积的60%[23]。因独特的气候，高寒草甸、高寒草原、高寒荒漠为青藏高原主要的三种草地类型(图1)[24]。

图1 研究区草地类型及外业调查样点空间分布Fig.1 Grassland types and spatial distribution of field investigation samples in the study area

1.2 样地数据获取

地面实测数据的获取时间为2003—2018年的7—9月，主要获取了生长季的草地植被盖度数据，总计获得4 376个样点(图1)。在研究区内设定面积为100 m×100 m的样地，选择植被生长状况均一、地势平坦、具有代表性的典型群落地段[25]。依据样地内草地植被的变化状况，每个样地选择3～5个样方采用目测法估计草地植被的盖度。在采样过程中，记录通过GPS定位的每一样方中心点的经纬度、高程和优势种、草地植被盖度、草层高度等指标[26]。一个样地的草地盖度值由5个样方观测值的平均值代表。剔除明显异常的数据，使数据采样点信息一致，最终得到有效样本记录4 355条用于模型构建(表1)。

表1 2003-2018年青藏高原地区天然草地盖度统计结果Table 1 Statistical results of natural grassland coverage on the Tibetan Plateau from 2003 to 2018

1.3 MODIS数据获取及处理

本研究下载了覆盖整个研究区的2001—2019年5—9月的MCD43A4反射率产品，时间分辨率为每天，空间分辨率为500 m。本研究共利用7个波段，每个波段所对应的波谱范围依次是620～670 nm(Band1)，841～876 nm(Band2)，459～479 nm(Band3)，545～565 nm(Band4)，1 230～1 250 nm(Band5)，1 628～1 652 nm(Band6)和2 105～2 155 nm(Band7)。

利用MODIS数据重投影工具进行数据格式转换和投影定义，将数据处理为TIFF格式和WGS84投影。在ArcMap 10.2软件中利用extract by mask工具裁剪出研究区范围，用于后续分析，同时，采用最大值合成法(Maximum value composite,MVC)，得到了月最大反射率数据集[27]。

1.4 生态环境指标及预处理

地形数据来源于地理空间数据云(http://www.gscloud.cn)，下载航天飞机雷达地形任务(SRTM) V4.1 TIFF数据，将SRTM数字高程模型(DEM)数据重采样至500 m空间分辨率，以和MCD43A4产品匹配。利用ArcMap 10.2软件计算得到研究区经度、纬度、海拔、坡向和坡度的空间分布数据[28]。

土壤数据来源于中国土壤特征数据集(http://globalchange.bnu.edu.cn/research/soil)，其中包括土壤深度0～30 cm的碱解氮(Alkali-hydrolysable N,AN)、有效钾(Available K,AK)、有效磷(Available P,AP)、容重(Bulk density,BD)、砂土含量(Clay fraction,CL)、砾石含量(Gravel,GRAV)、PH值(PH)、孔隙度(Porosity,POR)、含泥量(Silt fraction,SI)、含砂量(Sand fraction,SA)、土壤有机质(Soil organic matter,SOM)、全钾(Total potassium,TK)、全氮(Total nitrogen,TN)及全磷(Total phosphorus,TP)。同时使用本研究区边界对土壤数据进行裁剪，以便进一步分析。

温度和降水数据来源于国家地球系统科学数据中心共享服务平台(http://www.geodata.cn)。为了与遥感植被指数产品MCD43A4的时间分辨率相匹配，对应该产品周期将日数据处理为月累积降水和月平均温度。利用该数据使用插值软件ANUSPLIN Version 4.3 进行空间插值。空间插值数据选择500 m空间分辨率和Albers投影。拟合降水表面时，将经度、纬度和海拔作为三个独立的样条变量；拟合温度表面时，高程作为独立协变量。插值结果通过裁剪得到本研究区的气候空间数据集，进而提取采样地对应的降水和温度值。

1.5 变量及筛选

本研究所涉及的变量包括4类，共55个。其中，遥感植被指数有31个(表2)，气象指标包括年降水量(prec)、年平均温度(tmp)、大于0度年积温(GDD)、1月至采样月累积降水(per_samp)和1月至采样月累积温度(tem_samp)，地形指标包括海拔(DEM)、坡位(TPI)、坡向(Aspect)、曲率(Curvature)和坡度(Slope)，土壤指标包括AK，AN，AP，BD，CL，GRAV，pH，POR，SA，SI，SOM，TK，TN及TP。为了精简变量，使用最小绝对压缩变量筛选(Least absolute shrinkage and selection operator，LASSO)方法来选择敏感变量。这种方法可以显著减少变量的数量，降低共线性，提高模型精度，还可以消除模型的冗余性，精简模型。LASSO回归分析是基于最小二乘法的原理，利用模型中各变量的回归系数构造一个函数作为惩罚函数，并对各变量的系数进行压缩，直到模型的残差平方和最小，去掉系数降为0的变量。因此，LASSO回归分析方法被广泛应用于高维数据的变量筛选[29-30]。

表2 基于MODIS数据的草地植被指数变量Table 2 Grassland vegetation index variables based on MODIS data

1.6 模型构建与分析

本研究采用单因素回归、人工神经网络(Artificial neural network,ANN)、SVM和随机森林(Random forest,RF)4种模型对草地盖度进行建模。根据帕累托原则[31-32]，首先将4 355个样本数据以8∶2的比例分为两部分，其中3 484个样本用于建模，871个样本用于空间模拟结果的精度验证，再将3 484个样本数据以8∶2的比例分为两部分，其中2 787个样本作为训练集数据，697个样本作为测试集数据(图1)。

1.6.1传统回归模型 本项研究首先使用回归分析法构建草地盖度模型，利用LASSO筛选出的变量作为自变量，分别与草地盖度建立对数、线性、乘幂和指数4类回归模型[27]。

1.6.2机器学习模型 本文共采用3种机器学习模型，分别为RF，SVM，ANN模型。RF模型是基于分类树算法的一种算法。RF模型利用bootstrap进行抽样，抽出的样本用于回归树的构建。对训练样本进行连续筛选，得到最小残差平方和，最后形成一棵完整的树[33]。本研究在MATLAB 2019a中，通过Tree Bagger工具包构建RF模型，主要测试决策树的数目，在ntree中设置以100为间隔，从100开始增加到2 000，共迭代20次，当十折交叉验证误差最小时，确定为最优ntree，用于最终模型的构建。SVM是一种具有相关学习算法的监督学习模型，由高维或无限维空间中的一组超平面构成，可以用于分类、回归和其他任务。本研究在MATLAB 2019a中构建SVM回归，选择径向基函数(RBF)作为核函数，采用遗传算法(GA)智能寻找最佳cost (bestc)和gamma (bestg)参数。ANN是由输入层、输出层和一个或多个隐含层组成的多层网络结构[34]。本研究通过MATLAB软件的神经网络工具箱完成BP神经网络建模及其验证，设置隐含层神经元个数以10为间隔，从10增加到100，共迭代10次，取交叉验证误差最小时的neuron个数为最优neuron，构建最终模型。

1.6.3模型评价指标 模型的适应度由训练数据集的决定系数(R2)来衡量，预测精度的评价指标为相关系数(r)、均方根误差(RMSE)和均方根误差变异系数(CVRMSE)[35-36]。具体计算公式如下：

式中，E(yi)和yi分别为y的实测值和预测值，n等于样本量(这里n=697)。RMSE越低，模型精度越高[37]。

1.7 草地盖度空间分布与动态变化

2001—2019年最大草地盖度的时间变化采用坡度线性趋势模型分析，计算各像元的植被盖度变化[40-41]。当slope>0时，草地盖度呈增加趋势；当slope<0时，草地盖度呈降低趋势。计算公式为：

式中，n等于总年数(19)；i是2001年至2019年的年份序号(即1～19)；其中，coveri为最优模型模拟的第i年年最大草地盖度。

利用F检验分析草地变化的显著性[42]。如果F>F0.05(1,n-2)，在95%置信水平下变化显著，此处F0.05(1,14)=4.60。研究区草地盖度变化分为4类：显著增加(slope>0,F>4.60)，增加(slope>0,F<4.60)，减少(slope<0,F<4.60)和显著减少(slope<0,F>4.60)。F检验公式为:

式中，n等于总年数19，R2为每个像元的草地盖度与时间序列的相关系数r的平方；是1～19的平均值，此处等于10。coveri为最优模型模拟的第i年年最大草地盖度，其中，为2001—2019年平均最大草地盖度。

利用RF最优模型，反演出每年生长季草地盖度最大值。结合slope线性趋势分析和F检验分析2001—2019年草地盖度的变化。在ArcMap 10.2中，利用Reclassify工具对19年草地盖度增减情况进行分级，分为增加、显著增加、减少、显著减少四级[25]。

2 结果与分析

2.1 变量筛选结果

通过LASSO筛选，共选出10个敏感变量用于天然草地植被盖度模型的构建与分析。其中，包括2个植被指数，分别为LSWI，RVI；2个气象指标，分别为prec和tem_samp；1个地形参数，为Slope；5个土壤属性指标，分别为AN，AP，BD，POR，TP。由此可见，土壤、地形、气象和植被指数这四类变量与草地盖度均具有密切关系。

2.2 单因子模型分析

表3是利用LASSO筛选出的变量建立的回归模型及精度分析结果。由表3可以看出，在所有草地盖度单因子参数模型中，基于LSWI和RVI的单因子模型精度优于其他模型，其R2均为0.52，RMSE介于15.56%～18.21%；其次为基于气象指标的单因子模型，R2介于0.02～0.39，RMSE介于17.77%～23.06%；而基于土壤、地形等因素的草地盖度单因子参数模型精度均较低，R2介于0.1～0.26，RMSE介于19.35%～23.12%。

表3 草地盖度与植被指数的单因子模型精度评价结果Table 3 Evaluation results of single-factor model accuracy of grassland coverage and vegetation index

LSWI和RVI与草地植被盖度之间有较好的线性相关关系。基于RVI的线性模型的拟合决定系数(R2=0.52)高于指数模型(R2=0.50)、对数模型(R2=0.47)及乘幂模型(R2=0.42)，基于LSWI的线性模型的拟合决定系数(R2=0.52)高于指数模型(R2=0.47)，而基于RVI的均方根误差低于LSWI。由此可见，基于RVI的线性模型为青藏高原草地植被盖度的最优单因子模型。

2.3 机器学习模型分析

从分析结果(表4)可以看出，3种机器学习模型精度均高于单因子最优模型，R2提高了0.09～0.16，RMSE降低了1.52%～2.81%。RF模型的决定系数(R2=0.68)高于SVM模型(R2=0.66)和ANN模型(R2=0.61)，而RF模型的均方根误差最低(RMSE=12.75%)，同时RF模型的CVRMSE为16.81%，代表其有较好的预测能力。研究结果说明，本研究建立的RF模型预测值与实地测量的草地盖度值非常接近，可见RF模型优于其他模型，为青藏高原草地植被盖度的最优估测模型。

2.4 草地盖度空间动态变化分析

图2为利用RF模型反演的2001—2019年草地年最大盖度的平均结果。总体上看，青藏高原草地植被覆盖状况较好，东部地区草地覆盖度高，西部地区覆盖度低。低盖度(<40%)草地主要分布在西藏自治区的中部高海拔地区，其面积占比为16%；较低盖度(40%～60%)草地多分布在青海省西部和西藏自治区的周边，其面积占比为33%；较高盖度(60%～80%)草地则主要分布在青藏高原中部地区，与中低盖度草地交错分布，其面积占比为36%；高盖度(>80%)草地多分布在四川省、甘肃省和青海省的东部地区，面积占比为15%。高和较高植被覆盖区域面积占青藏高原草地总面积的15%和36%，共占青藏高原草地总面积的51%，说明青藏高原植被生长季大部分区域的草地植被覆盖率较高，且植被覆盖状况良好。

表4 青藏高原草地植被盖度机器学习模型精度评价结果Table 4 Accuracy evaluation results of machine learning model for grassland vegetation coverage in Qinghai-Tibet Plateau

利用最大值合成法，通过对7—9月的最大草地盖度趋势分析的结果表明(图3)，2001—2019年青藏高原大部分地区草地盖度呈现增加趋势，呈增加趋势区域占比为55.4%，减少区域占比为44.6%。草地盖度减少大部分发生于研究区中心地带，少数分布于四川省部分地区；草地盖度显著减少的地区主要在西藏自治区的东南部地区和西藏自治区北部的少部分地区，面积占比为5.2%；青藏高原的东、西部地区多为草地盖度增加区域；青海省的东北地区为盖度显著增加区域，另有少数分布在西藏自治区的西南部地区，面积占比为8.1%。

图3 2001—2019年青藏高原草地盖度变化趋势Fig.3 Trends of grassland coverage in the Tibetan Plateau from 2001 to 2019

3 讨论

本研究对比分析了单变量参数模型与非参数模型性能，确定了青藏高原地区最优草地盖度反演模型。单因素模型统计分析的结果表明，土壤、地形、气象和植被指数这4类变量均与草地盖度具有密切的关系。土壤、地形、气象因素对草地盖度有一定程度的影响，土壤属性指标与草地盖度的R2介于0.01～0.26，地形参数R2为0.02，气象因素R2达0.38。由此可见，单因素模型不能反映多种因素的综合影响，因此在大范围的草地盖度分析时基于单因素构建的回归模型的误差偏大。本研究采用了3种常见的机器学习模型。RF模型对参数的调整变化并不敏感，具有极佳的性能。经研究发现，噪声数据对RF模型影响较小，具有鲁棒性[43-45]。与ANN，SVM模型相比，RF模型的训练集和测试集精度存在一定差异。SVM模型更适合处理小样本数据，但当样本量过大时，SVM模型需耗费更多运算时间用于参数调整。参数调整对ANN模型的影响较大。最优机器学习的RF模型自变量包括LSWI和RVI 2个植被指数，prec和tem_samp 2个气象指标，1个坡度地形参数，以及AN，AP，BD，POR，TP5个土壤属性指标。这些变量涵盖了可代表草地生态环境的多种属性，各指数之间具有优势互补的特点，因此利用LASSO变量筛选方法选出的最优变量组合方式构建的RF模型可以更好的反演整个青藏高原天然草地盖度的变化。然而，基于RF的最优模型仍然存在一定的局限性和不确定性。首先，青藏高原地区草地面积较大、地形复杂，大部分样地分布在高原东部地区，而中西部地区由于受道路和海拔的限制，样地数量有限，而地面采样点和遥感数据的时间匹配也存在很大不确定性，使得构建的模型存在一定的误差。同时，青藏高原气象站分布不均匀，气象资料的空间插值存在一定误差[46]。其次，与单变量参数模型相比，RF模型算法更灵活，由基于高维数据训练的大样本决策树组成，具有较强的数据误差容忍度，但通常需要大量的标记和地面测量数据。因此，模型仍有一定的局限性和不确定性。

本研究分析了2001—2019年青藏高原草地盖度空间分布格局与动态变化特征，针对青藏高原草地盖度变化，学者已展开了大量研究。如唐志光等[47]研究表明三江源地区草地盖度空间分布呈现出东高西低的特点；于伯华等[48]的研究结果表明，雅鲁藏布江流域草地盖度总体上呈现出上升趋势。丁明军等[49]基于NDVI反演了1982—2009年青藏高原草地植被盖度的年际变化，高原大部分地区草地盖度呈现增加趋势。本研究结果与前人在该地区所做研究结果基本一致，均得出该地区2001—2010年内草地盖度呈现整体增加的趋势。马琳雅等[50]利用MODIS植被指数数据，分析了甘南州2001—2011年草地盖度空间变化特征，基于MODIS-EVI的对数模型为最优模型，R2为0.47。本研究与之相比精度较高，这得益于我们不仅考虑了遥感植被指数，还充分考虑了土壤、地形、气象等变量，能更好的反演草地盖度的时空变化。草地盖度不仅受气候变化的影响，同时人为因素也会对草地植被盖度产生很大的影响。2001—2019年青藏高原地区草地年平均最大盖度空间分布整体上呈现自西向东、自北向南递增的趋势，并且增加区域大于减少区域。这与近年来国家加强对草地退化的重视有关。自1999年以来，国家先后实施了一系列以草原综合治理和恢复为重点的工程措施，如退牧还草、草原生态保护补助奖励等措施[51]。因此，合理的人为干预对草地植被的恢复具有显著的作用。

在本项研究基础上，为提高模型精度，可通过科学手段改进现有模型，如:合理扩大观测范围、增加更多采样点以提高卫星影像的空间匹配性等。另外，传统的地表实测法费时、费力、精度低，且不能反映较大空间尺度的植被结构和空间变化等信息，近年来无人机(Unmanned aerial vehicle,UAV)技术普遍应用于地表观测，具有体积小、质量轻、灵活性高、可在特殊地区探测等优点[52-53]，因此利用基于UAV的大范围多样点的草地盖度快速监测方法，增加建模样本的数量和时空代表性，是未来改进研究区天然草地盖度反演模型的重要研究方向。此外，青藏高原由于地形复杂，东西部气温、降水变化强烈，受气候等诸多因素的影响，东西部的人口也存在巨大差距，气象因子和土壤、地形因子也存在一定的空间差异性。因此，深入研究和分析气象因素和人为因素对草地盖度的影响也是今后研究的重要内容。

4 结论

本研究基于2003—2018年青藏高原地区地形、土壤、气候和植被指数数据，对比分析了青藏高原天然草地植被盖度单因素回归模型和非参数模型的精度，研究了近19年(2001—2019)青藏高原草地盖度时空动态变化。结果表明，在单因子草地盖度遥感模型中，MODIS RVI与草地植被盖度的相关性最好，基于RVI构建的线性模型的R2为0.52，RMSE为15.56；在3种机器学习方法中，RF算法是青藏高原草地植被盖度反演的最优模型，R2达0.68，RMSE为12.75；RF模型的自变量包括2个植被指数、2个气象指标、1个地形参数及5个土壤属性指标，不仅有基于红光和近红外波段构建的传统植被指数RVI，也包括可以反映干旱、半干旱区域植被生长状况的LSWI，这些变量涵盖了优势互补的植被指数，因此RF模型在反演类型多样、地形、气候和土壤异质性强的青藏高原天然草地盖度变化方面具有独特的优势；2001—2019年间青藏高原地区草地植被盖度整体呈现出由西向东、由北到南的增加趋势，55.4%区域呈增加趋势，44.6%区域呈减少趋势。