基于MODIS数据与机器学习的青藏高原草地地上生物量研究

2022-10-24金哲人冯琦胜王瑞泾梁天刚

草业学报 2022年10期

金哲人，冯琦胜，王瑞泾，梁天刚

（兰州大学草地农业生态系统国家重点实验室，兰州大学农业农村部草牧业创新重点实验室，兰州大学草地农业教育部工程研究中心，兰州大学草地农业科技学院，甘肃兰州 730020）

我国幅员辽阔，拥有约960万km2的陆地面积，其中有400万km2的草原，占陆地面积的40%；同时，我国包含有被称为“世界屋脊”的青藏高原，占陆地面积的1/4，其中约70%是草原［1-3］。草地资源不仅是陆地上绿植中所占面积最大的可再生资源，而且对畜牧业的发展也起到重要的基础作用［4］。草地地上生物量（aboveground biomass，AGB）作为评估草地生态环境的重要载体可监测草地退化程度。草地地上生物量在评估区域碳循环中作用重大，是用来判断草地能否可持续健康发展的重要指标［5］。研究青藏高原草地地上生物量对当地草地资源的监测与良性发展具有积极影响。

运用遥感技术分析草地资源拥有比野外勘测更加高效和便捷的优势，可以对大范围的草地资源进行监测。国内外利用卫星影像对草地资源的研究已有很多，很多研究采用实测数据与植被指数结合进行草地生物量反演，对草地生物量变化进行评估。赵慧芳等［6］研究发现基于MOD13Q1数据的青海省草地地上生物量与植被指数NDVI具有从东南部向西北部逐级递减的空间分布特征。乔蕻强等［7］利用MODIS-NDVI和MODIS-EVI通过植被指数对肃南县草原生产能力进行估算，建立研究区草地估产模型。杜玉娥等［8］结合实测数据并利用MODIS植被指数NDVI和EVI资料，估算了2002-2009年青藏高原草地生物量变化。马青青等［9］结合玛曲2016年地面实测数据和MODIS影像资料，得出NDVI和地上生物量在6-9月相关性较高，其中9月的拟合水平（R2=0.5102）最高。陆荫等［10］基于实测数据和MOD13Q1植被指数产品，反演了2000-2019年甘南州草地地上生物量时空分布特征，结果表明MODIS EVI植被指数适合于当地草地生物量变化研究，模型决定系数R2=0.52149，RMSE=527.9 kg·hm-2。姚兴成等［11］利用MOD13Q1数据，建立对云南省草地AGB遥感估算模型，之后引入高度和盖度采用优化的估算模型分析其空间分布格局，模型的估算精度由原来的35.0%提升到43.7%。Gao等［12］利用MOD13Q1数据和随机森林模型与多因子（NDVI、纬度、经度、海拔、草地类型）结合，估算了2000-2017年青藏高原高寒草地地上生物量和盖度，并提高了估算精度。Zeng等［13］利用随机森林算法结合野外观测数据、MOD13A2遥感植被指数（NDVI、EVI）、气象数据和地形数据，对2000-2014年青藏高原草地地上生物量进行估算，所建立的随机森林模型在草地地上生物量估算中表现良好，能解释观测数据86%的变化。Zhou等［14］采用支持向量机（support vector machine，SVM）、随机森林（random forest，RF）和高精度地表建模（high accuracy surface modeling，HASM）3种机器学习模型对2001-2019年中国三江源地区草地AGB进行了模拟，HASM模型比RF和SVM模型表现更为优秀（R2=0.8459＞0.72＞0.5858；RMSE=29＜41＜56）。Chen等［15］使用多层感知器神经网络的机器学习方法，结合现场观测和气候数据，从Sentinel-2影像中估算奶牛场牧草地上生物量，大约60%的生物量变化得到了解释。通过以往研究发现，目前所选用的遥感数据多为月数据，实测数据与遥感影像时间匹配跨度大，会存在实测数据与遥感数据不完全匹配的问题；并且构建草地地上生物量的模型精度较低。

基于上述原因，本研究基于MODIS逐日反射率数据构建的多种植被指数，比较分析了不同机器学习方法反演青藏高原草地地上生物量的适宜性，评估了2000-2020年青藏高原草地地上生物量的变化，研究结果有助于对草地资源变化状况进行了解，为青藏高原地区草地生物量估算提供科学依据，服务于对青藏高原草地地上生物量的相关研究。

1 材料与方法

1.1 研究区域

青藏高原（25°59′26″-40°1′6″N，67°40′37″-104°40′43″E）位于我国西部，东至祁连山东端、横断山等山脉东缘，南抵喜马拉雅山等山脉南麓，西起兴都库什山脉和帕米尔高原西缘，北至帕米尔高原北缘、西昆仑山和祁连山山脉北麓，平均海拔约4320 m，总面积为3.08×106km2，其中我国境内的青藏高原面积约2.58×106km2［16］。以青藏高原地区为研究区域（图1），整个青藏高原地区天然草地面积约为1.5×108hm2，在高原气候的特殊影响与作用下，整个地区的草原类型以高寒草甸和高寒草原为主［17］。青藏高原不仅影响我国的气候与生态环境，甚至对世界的环境变化也具有重要作用，其中的草地资源尤其突出表现在维持该地区生态系统的稳定性方面，发挥了高原生态屏障的保护功能［18］。青藏高原的高寒草甸、高寒草原等大面积的草地资源作为高原本身的一道生态屏障，对防控当地土地沙化和水土流失具有积极的影响［19］。

1.2 数据采集与预处理生物量数据集

本研究中所用的草地外业调查数据来自2003-2020年5-9月青藏高原的实地观测调查。样地设置在地势平坦、草地植物群落均一性较好的地段，每块样地随机选择3～5个大小为0.5 m×0.5 m的样方。通过GPS测定样方的经纬度和海拔，并记录样地所隶属的行政区、样方内的物种名称、盖度、高度等植物群落信息。地上生物量采用齐地面刈割的方法，现场称其鲜重，之后带回实验室在70℃烘干至恒重称其干重。在本研究中，共采用4823条采样点的实测数据来评估该地区的草地生物量（图1），其中有1338条采样点的实测数据是使用Xia等［20］估算青藏高原草地生物量和周转时间的研究中所使用的数据。

气象数据来源于国家科技基础条件平台-国家地球系统科学数据中心（http：//www.geodata.cn）提供的2000-2020年逐月气象差值数据集。该气象数据集包括逐月平均气温（℃）、逐月最低气温（℃）、逐月最高气温（℃）和逐月降水量（mm）数据，空间分辨率为0.0083333°（约1 km）。数据是据CRU（气候研究中心Climatic Research Unit）发布的全球0.5°气候数据以及WorldClim发布的全球高分辨率气候数据，通过Delta空间降尺度方案在中国地区降尺度生成。并用496个独立气象观测点数据进行验证，验证结果可信［21］。

表1 本研究中使用的植被指数Table 1 The vegetation index used in this study

1.3 机器学习方法

本研究主要运用R语言的caret包（classification and regression training）进行机器学习模型的构建，该包提供的工具有数据分割、预处理、特征选择、使用重采样的模型调整、变量重要性评估等。本研究使用了caret包提供的机器学习方法包括：基于样条函数的广义加性模型（generalized additive model using splines，bam）、贝叶斯规整化神经网络（bayesian regularized neural networks，brnn）、提升树（boosted tree，bstTree）、弹性网络（elasticnet，enet）、基于径向基函数核的高斯过程（gaussian process with radial basis function kernel，gaussprRadial）、广义线性模型网（generalized linear model net，glmnet）、k-临近算法（k-nearest neighbors，kknn）、基于正向选择的线性回归（linear regression with forward selection，leapForward）、单调多层感知机神经网络（monotone multi-layer perceptron neural network，monmlp）、非负最小二乘（non-negative least squares，nnls）、岭回归（ridge regression，ridge）、基于径向基函数核的支持向量机（support vector machines with radial basis function kernel，svmRadialSigma）、稀疏偏最小二乘（sparse partial least squares，spls）、条件推断随机森林（conditional inference random forest，cforest）、平行随机森林（parallel random forest，parRF）、分位数随机森林（quantile random forest，qrf）、基于ranger的随机森林（random forest，ranger）、基于Rborist的随机森林（random forest，Rborist）、基于randomForest的随机森林（random forest，rf）、基于randomForest和RRF的正则化随机森林（regularized random forest，RRF）、基于RRF的正则化随机森林（regularized random forest，RRFglobal）、极端的梯度提升1（extreme gradient boosting 1，xgbDART）和极端的梯度提升2（extreme gradient boosting 2，xgbTree）共23种，探索适合特殊地区青藏高原的机器学习方法。各个方法具体描述见caret包的说明文档（https：//topepo.github.io/caret/index.html）。

1.4 变量的选择

本研究在选取变量的时候考虑了两个方面，一方面是植被指数，一方面是温度降水，本研究认为温度降水对植被敏感，所以本研究用了温度和降水数据，通过查阅相关文献使用了大家常用的反映植被生长状况的植被指数。

将2000-2020年的青藏高原各个相关数据集整合，变量共包含MCD43A4产品的7个波段（Band1、Band2、Band3、Band4、Band5、Band6、Band7）、10种植被指数（NDVI、EVI、EVI2、DVI、RVI、SAVI、MSAVI、OSAVI、SATVI、NDPI）、1-12月的每月降水量（prec_01、prec_02、prec_03、prec_04、prec_05、prec_06、prec_07、prec_08、prec_09、prec_10、prec_11、prec_12）、年降水量（prec_all）、1-12月的每月最高温度（tmax_01、tmax_02、tmax_03、tmax_04、tmax_05、tmax_06、tmax_07、tmax_08、tmax_09、tmax_10、tmax_11、tmax_12）、1-12月的每月最低温度（tmin_01、tmin_02、tmin_03、tmin_04、tmin_05、tmin_06、tmin_07、tmin_08、tmin_09、tmin_10、tmin_11、tmin_12）、

1-12月的每月平均温度（tp_01、tp_02、tp_03、tp_04、tp_05、tp_06、tp_07、tp_08、tp_09、tp_10、tp_11、tp_12）和年平均温度（tp_year），一共67个变量。

本研究采用R语言的CAST包中的前向特征选择方法（the forward feature selection，ffs）进行变量筛选，此筛选方法先基于全变量中任意2个变量训练模型，选出精度最高的模型所用变量，之后逐个添加特征变量，并选取模型精度提高最大的特征变量，当模型精度不再提高时停止变量筛选，最终得到最优特征变量组合。

1.5 草地地上生物量变化趋势分析

采用Theil-Sen median趋势分析、Mann-Kendall以及Hurst指数方法，研究青藏高原覆盖区域草地生物量的空间分布特征、时间变化趋势、变化趋势特征和可持续性特征［32-33］。Theil-Sen median是一种非参数统计的趋势计算方法，效果稳健并可以减少数据异常值的影响［32］。Mann-Kendall作为一种非参数统计检验方法，优点是不需要样本服从特定的分布，可以直接检验时间序列的单调变化趋势［34-36］。Theil-Sen median趋势分析和Mann-Kendall趋势检验结合可判断长时间序列数据趋势，能够有效分析青藏高原草地生物量时间序列的变化趋势［37］。Hurst指数是定量描述时间序列信息长期依赖性的有效方法［37］。Theil-Sen median趋势分析、Mann-Kendall以及Hurst指数（即H值）方法的计算公式详见文献［32］。将结果划分为轻微降低、显著降低、稳定不变、轻微增加和显著增加5种类别。H值取值越接近1，持续性越强；H=0.5，说明不存在长期相关性；H越接近0，反持续性越强。

2 结果与分析

2.1 采样点数据分析

按草地类型统计青藏高原逐年实地采样数据的结果表明（表2）：在采样数据中，风干重由高到低排序，前3位依次是山地草甸类、沼泽类和高寒草甸类，分别达到2001.48、1812.50和1533.51 kg·hm-2。风干重较低的3位是高寒荒漠草原、温性荒漠和高寒荒漠，分别为406.76、393.85和304.06 kg·hm-2，同时高寒荒漠的标准偏差也最低，只有89.62 kg·hm-2。山地草甸的盖度达到87.13%，高寒草甸的盖度达到83.15%，沼泽的盖度达到75.50%，是所有草地类型中盖度较高的前3位。盖度最低的3位是温性草甸草原、高寒荒漠草原和温性荒漠，分别为35.62%、29.16%和25.24%。各草地类型的高度大多集中在4～15 cm，低地草甸平均高度最高，达到了26.12 cm，高寒草甸草原平均高度最低，只有4.74 cm。

(3) 电网输配电效益下降问题。分布式电源使电网输配电效益下降几乎是不可避免的,主要表现为:分布式发电影响电网的售电量,从而影响电网的效益;增大公共电网网损率,降低公共电网的设备利用率。

表2 不同草地类型的调查样本数据基本情况Table 2 Survey sample data and basic statistics of different grassland type

2.2 变量相关分析

利用ffs算法筛选变量的结果如图2所示。经过筛选得出变量“prec_05”、“prec_06”、“tp_12”、“NDPI”、“prec_04”、“tmax_01”、“prec_08”、“prec_12”与生物量关系较大，其中筛选出的8个变量中有5个与月降水量有关，它们各自之间的相关系数关系如图3所示。结果表明：与生物量相关性最大的变量是“NDPI”，相关系数值为0.589，其次是“prec_05”和“prec_04”，即4、5月的降水量，与生物量的相关性达到了0.400以上。

2.3 基于机器学习的最优模型

运用R语言中caret包，将67个变量全部带入所选用机器学习模型中运算与生物量的关系，结果如图4所示。全变量与生物量决定系数最高的模型是rf，R2达到了0.6360。决定系数较低的后3位模型是glmnet、nnls和leapForward，R2都在0.46以下。按照R2由高到低排序，前8位模型依次为rf、RRF、RRFglobal、Rborist、parRF、ranger、qrf、cforest，它们都属于随机森林类方法，R2都在0.61以上。其次R2较高的模型是xgbTree、xgbDART、kknn、svmRadialSigma和bstTree，这些模型的R2也都达到了0.57以上。说明对于研究全变量和生物量的关系，随机森林rf模型为最优模型。

基于筛选出来的8个变量构建机器学习模型，根据图4可知，在所用的机器学习的方法中，R2由高到低前5位依次是Rborist、RRF、ranger、RRFglobal和parRF，R2分别为0.6484、0.6477、0.6474、0.6473和0.6460。R2在0.40以下的模型有spls、glmnet、enet、ridge、leapForward和nnls模型，说明这些方法运行效果欠佳，并不适合研究筛选后的变量与生物量的关系。结果表明Rborist、RRF、ranger、RRFglobal、parRF、rf、qrf和kknn这8种方法基于筛选变量的模型比全变量模型的R2值高，说明减少变量可以消除数据冗余，从而提高模型精度。由此本研究选取出最优的机器学习模型为Rborist。筛选前后各机器学习方法参数取值如表3所示。

2.4 青藏高原2000-2020年生物量变化分析

基于选择出的最优机器学习模型，本研究模拟了青藏高原2000-2020年草地地上生物量（图5）。结果表明，青藏高原的草地地上生物量由西北向东南方向逐步增加，东南部的草地地上生物量普遍在900 kg·hm-2以上，要明显好于西北部。

采用Theil-Sen median趋势分析和Mann-Kendall方法，可有效反映出青藏高原地区2000-2020年的草地生物量的时空变化趋势，如图6所示。研究结果表明，青藏高原草地生物量轻微降低区域占草地总面积的9.83%、显著降低区域占比1.33%、稳定不变区域占比64.33%、轻微增加区域占比18.24%，显著增加区域占比6.26%。草地显著降低的面积最少，稳定不变的面积最大，并且从图中可看出草地生物量降低的区域主要集中在青藏高原东部和中部，增加的区域分布在青藏高原东北部。

将Theil-Sen median趋势分析和Mann-Kendall检验结果与Hurst指数相结合，得到青藏高原草地地上生物量变化趋势与持续性的耦合信息，如图7所示。研究结果表明，青藏高原草地地上生物量变化趋势不稳定的区域占比61.38%，持续性轻微降低的区域占比4.67%，持续性明显降低的区域占比1.19%，持续性稳定不变的区域占比19.71%，持续性轻微增加的区域占比7.87%，持续性明显增加的区域占比5.18%。总的来说，青藏高原2000-2020年草地地上生物量增加的区域要大于降低的区域，总体向好发展。

将草地类型按青藏高原行政区划统计（表4），结果表明各省份不同草地类型变化趋势不确定的区域面积占比较大，其中甘肃省、青海省和四川省持续性增加的草地区域面积都达到了20%以上，四川省持续性降低的草地区域面积超过了10%。山地草甸类在各省份均持续性轻微降低，其中西藏持续性轻微降低的区域面积占比最大，达到15.57%。甘肃省持续性明显增加的草地类型前3位依次是温性荒漠草原类、高寒草原类和温性荒漠类。西藏各草地类型生物量均有所降低，其中温性草甸草原类持续性轻微降低的区域面积达到20%以上，持续性增加总面积比率与退化的比例相持平。整体上各个省份大部分区域的各草地类型在持续性增加。

表4 2000-2020年青藏高原各类草地变化特征分析Table 4 The change characteristic analysis of different grassland from 2000-2020 on the Tibetan Plateau（%）

续表Continued Table

3 讨论

3.1 机器学习模型的优劣

本研究基于地面实测数据和MCD43A4产品构建的植被指数，使用ffs筛选得到变量“prec_05”、“prec_06”、“tp_12”、“NDPI”、“prec_04”、“tmax_01”、“prec_08”、“prec_12”和生物量的关系较大，运用最优模型Rborist，对2000-2020年青藏高原地区的草地生物量进行反演，并对青藏高原2000-2020年的草地地上生物量时空变化趋势进行分析。Meng等［38］基于遥感数据和机器学习算法对2011-2016年甘南地区草地地上生物量建模，用随机森林算法构建了最优估计模型，该模型可解释牧草生长季AGB的89.41%的变化。Zeng等［13］利用随机森林算法对2000-2014年青藏高原草地生物量反演，该研究在RF建模过程中选择MOD13A2遥感植被指数（NDVI和EVI）、气象变量（年平均温度和年平均降水量）和地形变量（海拔、坡度和坡向），所建立的随机森林模型在草地地上生物量估计中表现良好，能解释观测数据86%的变化。Lin等［39］利用优化算法和多维特征改进高寒草地植被覆盖度估算，其中随机森林算法（R2：0.861，RMSE：9.5%）的植被覆盖度反演效果最好。Tang等［40］利用机器学习方法对2001-2020年黄河源地区反演草地生物量，研究结果表明随机森林模型优于其他3种模型。本研究主要是利用多种机器学习方法对草地生物量进行反演。同样发现随机森林方法效果最佳，R2由高到低排序前10位中有8位是随机森林的方法，前5位都是随机森林的方法。随机森林优点是可以准确描述生物物理参数与复杂环境因素之间的复杂非线性关系，并且易于实现，因为它的参数简单，预测变量不需要预处理，且具有较好的泛化性和鲁棒性［41］。本研究所选用的最优模型是随机森林方法中基于Rborist包的随机森林（random forest，Rborist）。Rborist包能够实现随机森林算法，并且拥有较高性能。这个包是Arborist项目的衍生品，Arborist项目是一个针对各种决策树方法的多语言项目。Arborist在实现随机森林的基础上具有可扩展和高性能特点，支持回归和分类决策树模型，能够在训练和测试中利用多核并行运算，并且具有集群友好性。

3.2 研究结果的一致性

植被指数对地表植被状况敏感，张雅等［42］基于Landsat 8数据提取常用的植被指数，利用统计分析方法建立紫泥泉牧场阴坡和阳坡的草原生物量遥感估算模型，并进行了生物量空间反演与验证，其中NDVI与生物量的相关性最高。葛静等［43］基于ADC（农业数字相机agricultural digital camera）和MODIS遥感数据对黄河源区的高寒草地地上生物量进行监测研究，结合MODIS NDVI构建草地地上生物量反演模型。Wang等［31］提出的NDPI植被指数表现优异。Gan等［44］利用MODIS MOD09A1产品提取的植被指数结合多种方法监测2009-2013年黄淮冬小麦（Triticum aestivum）返青期（green-up dates，GUDs），结果表明NDPI与地面数据的一致性优于其他植被指数。Su等［45］利用Sentinel-2卫星遥感影像，通过对植被光谱反射率与降尘吸收率（the amount of dust absorption，ADA）的相关性分析，建立降尘反演模型，获取常绿植被的降尘分布，结果表明NDPI更适合建立反演模型。An等［46］基于MODIS MOD09A1和不同算法分别监测青藏高原植被生长季的开始和结束时间，研究结果指出基于EVI和NDPI的物候指标时空稳定性更好。本研究选用的10种植被指数与草地生物量之间的相关系数如图8所示，可以看出NDPI与草地生物量关系更密切。原因在于NDPI使用加权的red-SWIR组合来代替NDVI中的红波段，克服了土壤背景异质性的不利影响；同时NDPI通过添加SWIR波段来响应叶片水分含量，由于水分在草本植被总重量中所占比例较高，因此NDPI能够捕捉叶片水分含量的变化，考虑了草地叶片水分变化对生物量的影响［47］。本研究结果与Xu等［47］认为NDPI在草地地上生物量估算中具有优越性相一致。

本研究经过筛选得出8个变量中有5个与月降水量有关，说明降水量对青藏高原草地生物量影响很大。卓嘎等［48］对2000-2016年青藏高原植被覆盖时空变化及其对气候因子的响应的研究结果表明降水量对植被覆盖的影响大于平均气温对植被覆盖的影响。Zhang等［49］基于MODIS MOD09A1产品和气象数据，监测了青藏高原的绿度变化趋势及其驱动因素，结果表明青藏高原东北部植被对降水相对敏感，而高原南部植被对温度相对敏感。Zeng等［13］利用随机森林算法估算青藏高原草地地上生物量，研究指出在整个研究区草地AGB与气温、降水呈显著正相关。草地AGB与年平均降水量的相关性为0.54（P＜0.05），远高于年平均温度（R=0.38，P＜0.05）。Shen等［50］利用NDVI和气象数据研究生长季前温度和降水对青藏高原中东部草地春季物候的影响，研究结果表明温度比降水对青藏高原草地返青期空间格局的影响更大。韩炳宏等［51］基于2000-2018年MODIS NDVI数据及气象数据，探讨了青藏高原地区NDVI的时空变化趋势及其与气温、降水的关系，研究发现青藏高原地区NDVI与同期气温和降水均显著相关，但与气温的关系更密切。综上，目前对影响青藏高原地区草地生物量的气象因素，有学者认为主要是气温，有学者认为是降水。出现差异的原因很可能是青藏高原不同地区对气温和降水敏感度不一样，而且研究的年份、使用的数据和研究的方法不同都可能导致结果不同［49］。

3.3 研究展望

本研究证明了NDPI在青藏高原地区反演草地生物量有出色的表现，并且证实了Rborist方法对青藏高原草地生物量建模的有效性。今后将继续深入研究NDPI的适用性，并进一步分析温度和降水与草地生物量的关系。在未来的工作中，考虑把土壤性质等其他环境因素纳入建模过程，以提高预测性能，并且最优模型的选取是否适合运用在其他地方还有待进一步验证。