基于机器学习算法的草地地上生物量估测——以祁连山草地为例

张子慧,吴世新,赵子飞,李向义,曾凡江,谢聪慧,侯冠宇,罗格平

1 中国科学院新疆生态与地理研究所荒漠与绿洲生态国家重点实验室,乌鲁木齐 830011

2 中国科学院太空应用重点实验室,北京 100094

3 中国科学院新疆生态与地理研究所新疆荒漠植物根系生态与植被修复重点实验室,乌鲁木齐 830011

4 中国科学院空间应用工程与技术中心,北京 100094

5 中国科学院大学,北京 100049

草地资源是我国陆地上面积最大的生态系统类型,具有防风、固沙、保土、调节气候、净化空气、涵养水源等生态功能,是草地畜牧业发展最重要的物质基础,对于全球生态环境、碳循环和维系生态平衡起重要的作用[1]。其时空分布格局可以反映草地植被的碳汇潜力[2]。草地地上生物量(aboveground biomass,AGB)通常用草地地上部分植被的干重表示,是全球碳循环的重要组成部分,是指导畜牧业生产管理的重要指标[3]。估算草地AGB是草地资源空间格局动态研究的重要内容,也是草畜平衡综合分析的基础,祁连山作为我国西部地区的重要生态保护屏障、黄河流域的重要水源产流地,也是我国生物多样性保护的优先区域,区域土地覆盖以草地为主,其生态保护事关我国西部地区的生态安全及经济社会可持续发展。及时准确地获取祁连山草地产量数据、草地植被生长状况、AGB大小及其空间分布,对区域复合生态系统功能的科学评估具有参考价值,对草地退化机理研究与治理、指导草地畜牧业生产,维护生态服务功能和区域可持续发展具有重要意义[4—5]。

目前,草地AGB反演面临着诸多挑战,数据源方面,遥感影像具有探测范围大、成本低、容易获取等优点,为通过间接测量的、非破坏性的手段反演草地AGB现状,遥感影像已成为生物量反演的主要数据源,目前使用较多的为MODIS、Landsat- 8和Sentinel- 2影像数据[6—9],其中Sentinel- 2数据具有更高的空间分辨率(10m),可以更好地与实测AGB数据相匹配,是目前常用的性能最好的星载遥感数据[7—11]。虽然卫星影像可以进行大范围的监测,考虑到卫星影像与实地采集样方间的尺度差异,有些细节特征还是不能突出,因此许多学者引入UAV数据用于减少多源数据间的尺度效应[12—15]。特征提取方面,国内外学者构建了一系列与AGB密切相关的植被指数用于提高草地AGB遥感反演的精度[10—11,16—18],比如增强型植被指数(Enhanced Vegetation Index,EVI)能够改善归一化植被指数(Normalized Difference Vegetation Index,NDVI)在高植被覆盖区的饱和现象[16]、转换型土壤调节植被指数(Transformed Soil-adjusted Vegetation Index, TSAVI)对土壤背景值敏感[17]、比值植被指数(Ratio Vegetation Index, RVI)在植被覆盖度高的情况下对植被十分敏感[18]等,不同植被指数能够反映不同的植被特征,根据研究区植被特征选取合适的植被指数用于AGB反演建模是极其重要的。反演算法方面,机器学习算法成为当前常用的构建模型的方法,许多学者利用数理方法构建草地AGB参数反演模型[19—24],且已被证实适用于草地AGB的反演中。目前常用的算法包括线性回归算法和非线性回归算法两大类,其中线性回归算法多以一元线性回归和多元线性回归为主[19—20],非线性回归算法多以支持向量机回归(Support Vector Regression, SVR)[21]和DTR[22]、RFR[23]、GBRT[22]和XGBoost[24]等基于决策树形成的回归算法为主,反演AGB过程中可能会遇到变量之间的关系并不是线性的,利用线性回归算法构建模型时受到限制,SVR算法可以用于解决非线性回归问题,但容易出现过拟合现象且计算复杂度高、可解释性不好,基于决策树的算法不仅可以用于线性或非线性的回归问题中,还能够很好的避免过拟合现象,是目前应用最为广泛的反演算法。DTR模型解释性好,可以使用图形化描述模型[22]；RFR是用训练数据随机计算出多颗决策树,实现简单,训练速度快,模型泛化能力强[23]；GBRT对异常值的鲁棒性非常强,在相对较少的调参时间下,预测的准确度较高[22]；XGBoost在代价函数中加入了正则项,能够控制模型的复杂度,可以防止过拟合现象还能够减少计算量[24]。因此,基于4种算法的独特优势,本文探讨了这4种算法反演祁连山草地AGB的适用性。以往的AGB反演研究取得了长足的进展,提出了许多新的方法和理论,综合以上内容依旧存在着三个重要问题,一是如何减少多源数据间的尺度差异；二是如何选取适合研究区内植被特点的植被指数；三是如何选择合适的反演算法用于构建草地AGB反演模型。

1 研究区概况及数据处理

1.1 研究区概况

研究对象为祁连山的草地生态系统,位于青藏高原东北部(图1,94°33 ′10.8″ —100°39′ 7.2″ E, 36°51 ′25.2 ″—39°54′ 37.2″ N),是河西走廊地区重要的生态保护屏障和主要的水资源涵养地[26]。研究区内平均海拔约为3800m,占地面积约为128300km2,年平均温度为0.6℃,年平均降水量约为400—700mm[27]。光照充足,冬季长而寒冷,夏季短而温和,干湿分明,常年受西风环流的影响,属于典型的高原大陆性气候[28]。水热条件随海拔的升高而发生变化,气候类型和植被类型表现为明显的垂直地带性差异[29]。主要植被类型有木本猪毛菜(SalsolaarbusculaPall.)、高山嵩草(KobresiaPygmaeaC. B. Clarke)、针茅(StipacapillataL.)、鬼箭锦鸡儿(CaraganaJubata(Pall.) Poir.)、二裂委陵菜(PotentillabifurcaLinn.)等[30]。主要草地类型有荒漠草原、高寒荒漠草原、高寒苔草草原、高寒嵩草草甸、高寒针茅草原等[31—33]。

图1 研究区示意图

1.2 数据来源与处理

1.2.1地面实测AGB数据

祁连山草地实测AGB数据获取时间为2021年7月20日—8月10日,正值植被生长期。共布设23个能够代表祁连山草原植被、地形及土壤等特征的采样区(如图1所示),按一定方向在样地中间位置设100m×100m样线(图2),记录每个采样区的基本信息,包括经度、纬度、海拔高度、盖度、株高等。为方便与Sentinel- 2 影像数据的分辨率(10m)相匹配,保证样方在采样区内覆盖主要草地类型且均匀分布,沿样线布置5个1m×1m样方,每个样方间隔距离为20m(图2),共计115个样方。收集样方内所有植物的地上部分,称其鲜重后在实验室内置于105℃下杀青,采用85℃烘干至恒重,获得每个样方的草干重,分别计算每个样方的生物量,获取115个样方的AGB实测数据。用于构建祁连山草地AGB反演的训练集和测试集。经实地考察,研究区群落由西北至东南方向具有明显差异,西北区域植被较稀疏,以荒漠草原为主,东南区域植被分布均匀且覆盖度均达85%以上,以嵩草草原为主。

图2 采样区示意图

1.2.2图像数据与处理

Sentinel- 2是高分辨率多光谱成像卫星,覆盖可见光、近红外等13个波段,其中包含3个红边波段(表1),能够反映植被的反射率光谱特征,重采样后分辨率可达到10m[6—9]。本研究中所需的Sentinel- 2遥感产品通过欧洲航天局(European Space Agency, ESA)官网(https://scihub.copernicus.eu/dhus/#/home)下载云覆盖率低于10%、与实地采样时间相近的21景影像,以及2016年植被生长期的影像。

表1 Sentinel- 2 MSI数据的参数信息

UAV数据通过大疆Phantom 4 Pro获取,无人机按照DJI GO 4设定的飞行路线自动飞行至覆盖整个采样区,设置航向重叠率为80%、旁向重叠率为70%,飞行高度为15m,并按一定的时间间隔(1.5s)垂直拍摄地表,共拍摄23个样地的UAV数据,空间分辨率为0.4cm。基于Sentinel- 2影像数据对UAV数据进行图像方向配准、尺度配准,提取两种影像数据特征向量,从而将UAV数据配准至Sentinel- 2数据。UAV数据作为联系Sentinel- 2遥感影像和1m×1m实地样方的媒介,可以在一定程度上减少Sentinel- 2 影像数据和地面样方之间存在的尺度差异。研究区DEM数据来源于美国地质勘探局(United States Geological Survey,USGS)提供的30m分辨率SRTM数据(https://lpdaac.usgs.gov/products/srtmgl1v003/)。

教学结束后，采用自制调查问卷对两组护生职业核心能力和人文素养各项指标进行测评。问卷包括主动学习力、学习兴趣、团队协作、语言表达能力、沟通能力、人文素养6个方面，满分100分，分值越高表明能力越强。问卷以不记名方式填写。共发放问卷127份，回收127份，有效回收率100%。

利用SNAP、ArcGIS软件对Sentinel- 2数据进行预处理,包括重采样、镶嵌、裁剪等(图3)。使用Pix4D Mapper对UAV数据进行拼接、地理位置叠加、空三加密等操作,得到数字正射影像,最后导入ENVI图像处理软件中进行影像裁剪等处理(图3)。

图3 数据处理流程图

1.2.3植被指数选取

为了证实Sentinel- 2数据在反演AGB方面的潜力,比较分析相关研究,在UAV数据辅助下,对处理后Sentinel- 2 影像数据提取了EVI、RVI、红边归一化植被指数(Red Edge Normalized Difference Vegetation Index, NDVI705)、TSAVI、红边简单比值(Red Edge Simple Ratio,SRre)和改进型简单比值(Modified Simple Ratio, MSR)6种植被指数(表2),用于反演祁连山AGB的空间分布。

表2 选取的植被指数

2 研究方法

建立祁连山草地AGB反演模型,进行模型的适用性分析,技术路线图见图4,图中y、y-y1、y-y1-y2-…-yn-1为单颗决策树的相关AGB输入值；y1、y2、yn、y1+Ω(f1)、y2+Ω(f2)、y1+Ω(fn)分别为单颗决策树的AGB输出值；F(x)为算法的输出值。

图4 技术路线图

本文利用了DTR、RFR、CBRT和XGBoost 4种算法构建祁连山草地AGB反演模型。这4种算法均基于分类与回归树(Classification and Regression Trees, CART)实现预测,CART采用二分递归分割技术,通过构建一个二叉树将样本进行递归划分,使用样本的最小方差作为分裂节点的依据。

2.1 回归模型算法

DTR在训练数据集所在的输入空间中,递归地将每个区域划分为两个子区域并决定每个子区域上的输出值,构建二叉决策树,根据样本中不同特征的信息纯度形成树的各个结点,输入样本数据流经整颗决策树,最终在叶子结点输出最终的结果[22]。RFR作为一种集成学习方法,采用了Bagging的思想,在原始训练样本集中有放回地重复随机抽取新的训练样本集合,基于多颗决策树引入了随机属性选择对数据进行预测,输出所有决策树输出的平均值,可以避免单颗决策树的局限性,通过降低模型方差来提高性能[21,23]。GBRT基于CART决策树及Boosting技术的集成学习算法,以串行的方式建立多颗决策回归树,是一种迭代的决策树算法,能够抑制决策树的复杂性,降低单颗决策树的拟合能力,再通过梯度提升的方法集成多个决策树,使用一阶泰勒公式优化模型,通过降低残差来提高性能[22,34]。XGBoost基于GBRT进行了改进,在传统的Boosting的基础上,使用了二阶泰勒公式优化模型,并在代价函数里加入了正则化项(公式(1)),能够降低模型的方差,用于控制模型的复杂度,使学习出来的模型更加简单,提高模型泛化能力[24,35]。

(1)

式中:Ω(ft)为第t颗树的正则化项；T为第t棵树的叶子结点数；wj为叶子结点j的权重向量l2范数；Υ和λ为XGBoost算法的参数。

2.2 精度评定

交叉验证的基本思想是将原始数据进行分组,其中80%的数据作为训练集,20%的数据作为测试集,首先用训练集对回归器进行训练,再利用验证集来测试训练得到的模型,作为评价回归器的性能指标,使用交叉验证是为了得到可靠稳定的模型,常见的交叉验证形式有K折交叉验证、留一交叉验证等。

为了分析4种算法反演祁连山草地AGB的适用性,本文利用5折交叉验证对模型进行评估,利用平均绝对误差(Mean Absolute Deviation,MAE,公式(2))、均方根误差(Root Mean Squared Error,RMSE,公式(3))和决定系数(R-squared,R2,公式(4)) 3个指标[21—24]用于模型精度评定。

(2)

(3)

(4)

3 结果与分析

3.1 变量相关系数矩阵分析

计算了植被指数与实测AGB间的Pearson相关系数(图5),发现AGB与各植被指数均表现为显著正相关(P≤0.01)。其中,EVI、TSAVI和AGB的相关性相对较低；TSAVI与其他植被指数的相关性也较低。造成该结果的原因是祁连山草地植被覆盖度较高,受土壤背景影响低,而EVI和TSAVI更适用于低覆盖区的植被研究中,能够降低土壤背景的影响,提高对植被的敏感性,主要与EVI和TSAVI本身所反映的植被特征有关。

图5 植被指数与AGB的相关系数矩阵

3.2 模型的构建与因子重要性评估

本文评估了4种算法反演祁连山草地AGB的因子重要性(图6),每种算法中各因子重要性表现也有其异同点。在RFR模型中,各因子间重要性差异较小,利用该模型反演AGB时考虑的因子较多,模型复杂度较高；而DTR模型中,各因子间重要性差异较大,利用该模型反演AGB时考虑的因子较少,可能会出现影响因子考虑不全的结果。RVI在4种算法中均占较高权重,对于反演祁连山AGB状况有着不可代替的作用；并且RVI对土壤背景变化的敏感性较弱,在植被覆盖度较高的情况下对植被信息反映敏感,表明祁连山草地的覆盖度较高； EVI和TSAVI 2种植被指数的重要性均较低也能够表明研究区植被覆盖度较高的现状,因为这两个植被指数更适合对低覆盖度的植被监测,所表现出的结果与它们和AGB间的相关性也有关。在XGBoost算法中,降低了EVI的影响,说明该植被指数不适用于祁连山草地AGB的反演,不能有效地反映出研究区草地的植被特征；明显突出了SRre的重要性,该指数与近红外与红边波段信息有关,这两个波段对植被信息反映比较敏感。

图6 AGB与各植被指数的因子重要性评估结果

3.3 模型适用性评价

通过比较4种算法中反演祁连山草地AGB最优模型的MAE、RMSE、R2和模型精度进行模型的适用性评价(表3)。其中,XGBoost模型具有最低的MAE(71.60 kg/hm2)和RMSE(99.74 kg/hm2),最高的R2(0.78)和精度(74.75%),反演效果最好,其次分别为GBRT、RFR、DTR。因为RFR、GBRT和XGBoost 3种算法均是在DTR的基础上进行了改进,其中,RFR采用的是bagging思想,而GBRT采用的是boosting思想,以串行的方式建立决策树；GBRT在优化的时候只用到一阶导数信息,XGBoost在代价函数中加入了正则项,与其他算法有其独特的先进性。XGBoost算法中突出的SRre植被指数,更适用于祁连山草地AGB的反演,对植被信息反映敏感,因此本文将选用XGBoost算法对祁连山草地AGB进行反演,用于进一步分析其空间分布状况。

表3 4种回归模型算法的适用性比较

3.4 地上生物量的空间分布

综合上述模型适用性的分析,本文利用XGBoost反演祁连山草地AGB,总体上呈现出由西北向东南增加的趋势(图7),平均AGB为925.43 kg/hm2,AGB分布主要受到水热条件的影响,靠近水源的区域AGB更大,水分条件差的区域AGB较低。东南区域靠近青海湖,属于高原大陆性气候,光照充足,具有良好的适宜植被生长的水热条件,同时在植被生长峰值期间受到西南季风的影响,带来的部分水汽适宜植被生长,因而AGB较高；西北区域水分条件较差,受季风影响弱,主要受到山间径流的影响,使得AGB沿水资源分布状况而发生变化,与东南区域相比,AGB较低。

图7 祁连山草地AGB空间分布图

4 讨论

上述研究结果表明,在UAV数据的基础上,利用Sentienl- 2 影像数据提取出的植被指数,选用XGBoost模型反演祁连山草地AGB空间分布的方法具有很强的应用潜力,在交叉验证的基础上,模型精度也在可接受范围,与其他研究相比[23],本文针对祁连山草地AGB的空间制图结果是可信的,该研究能够提供一定的可信度与理论支撑,但依旧存在四类重要问题值得讨论分析。

4.1 降水量的年度波动对AGB反演的影响

在空间分布上受到降水量差异影响,草地AGB表现出空间分异性,由于研究区内降水量存在年际差异,大范围下草地AGB状况会存在一定的滞后性,以往研究结论也提出相关论证[36—37]。本研究及将在后续工作中继续探讨降水量的年波动是否会在很大程度上影响XGBoost模型反演祁连山草地AGB空间分布的效果,以及祁连山草地AGB呈现出怎样的动态变化模式。

4.2 草地地上生物量反演模型的不确定性与误差分析

本文在反演祁连山草地AGB过程中主要存在以下3各个方面的不确定性:首先是多源数据间的尺度效应[11,17],Sentinel- 2数据的空间分辨率为10m,UAV数据分辨率为0.4cm,实地采集样方大小为1m×1m,使用Sentinel- 2数据和UAV数据协同反演祁连山草地AGB能够在一定程度上提高模型精度,但依旧会存在一些处理上以及算法上的误差。其次是在研究区内采样点布设较少,既对研究区AGB的反演会存在一定误差[20,38]；也可能会造成过拟合现象,由于XGBoost算法中添加了正则化项来控制模型的复杂度,与其它3种算法相比,能够提高模型反演的精度,但不能完全消除反演时存在的高估和低估问题[19]。第三是本文中所使用的Sentinel- 2数据由于受到云覆盖的限制,选择了与采样时间相近的影像用于研究,虽然时间不能完全匹配,但造成的这种不确定性对于祁连山草地AGB变化的影响是可以忽略不计的。

4.3 空天地一体化的耦合效应

利用UAV数据作为实现空天地一体化监测技术的中间媒介,将Sentinel- 2影像数据和地面采样数据进行像元尺度的匹配,尤其是对于西北区域,建群种为木本猪毛菜,散布在该区域且自身生物量较重,导致在实地采样过程中会受到较大影响,而Sentinel- 2数据的像元尺度是10m,计算出的植被指数为10m×10m范围内的均值,本研究中通过提取UAV数据的特征波段,分析其与Sentinel- 2影像数据间的差异,发现UAV数据能够观测到Sentinel- 2数据所不能表现出的细节特征。因此,UAV数据的引入能够降低匹配Sentinel- 2数据和实地样方间的误差,该结论与孙世泽等[20]、刘沁茹等[39]的研究结果一致。UAV数据可以减少尺度匹配过程中的误差,但这种误差是不能完全消除的,UAV数据获取过程中会受到GPS误差的影响,以及对UAV数据处理时也会受到系统误差的影响。

4.4 XGBoost模型反演草地AGB的适用性

本文所选用的XGBoost模型被证实适用于反演研究区草地AGB状况。在空间分布适用性上,XGBoost模型精度为74.75%,MAE为71.60 kg/hm2,RMSE为99.74 kg/hm2,R2为0.78。除多等[19]基于MODIS数据,利用GIS空间数据处理技术和多元统计分析方法反演草地AGB,R2介于0.1295—0.6929之间。黄家兴等[40]基于Sentinel- 2植被指数数据,利用二次曲线回归模型反演天祝县草地AGB,R2介于0.4019—0.6983之间,郭超凡等[41]基于Sentinel- 2影像数据,利用随机森林回归模型反演草地AGB,R2为0.74。与前人的研究结果相比,本研究所利用的XGBoost模型在反演祁连山草地AGB中具有一定的优势,反演效果较好,模型适用性较高。

5 结论

本文以Sentinel- 2 遥感影像计算的6种植被指数为特征,结合实拍UAV数据以及地面的115个样本点实测AGB数据,采用DTR、RFR 、GBRT和XGBoost 4种算法进行了祁连山草地地上生物量的遥感反演,以期为绿色推动祁连山草地畜牧业发展、维护生态平衡与生物多样性提供科学依据和理论支持。主要结论如下:

(1)6种植被指数均与实测AGB表现为显著正相关,其中RVI、NDVI705、SRre和MSR这4个植被指数与AGB的相关性更强,主要因为它们是由红色、红边和近红外波段对AGB较敏感的波段信息构成的；通过评估6种植被指数的因子重要性,RVI和SRre在祁连山草地AGB反演中占据重要地位,RVI本身对高覆盖区植被监测敏感,红边波段(SRre)的加入能够在一定程度上提高AGB反演模型的精度,减少植被在红色波段的饱和效应。

(2)通过对4种反演算法进行适用性评价,得出XGBoost算法反演祁连山草地AGB的效果最好,具有最强的适用性,MAE为71.60 kg/hm2、RMSE为99.74 kg/hm2,R2为0.78、模型精度为74.75%,且适用于祁连山草地AGB的长时间序列反演。

(3)UAV数据的应用能够减少Sentinel- 2数据(10m)与实地采样(1m×1m)间尺度效应的影响,UAV数据可以提供更详细的空间细节特征,便于Sentinel- 2影像中提取出的植被指数与地面实测AGB数据相匹配。

(4)利用XGBoost算法对祁连山草地AGB进行反演制图,祁连山草地平均AGB为925.43 kg/hm2,其空间分布上表现为由西北向东南递增的趋势,主要原因是受到水资源分布的影响。