青藏高原草地地上生物量估算

姚雨微,任鸿瑞

太原理工大学测绘科学与技术系,太原 030024

草地是陆地生态系统的重要组成部分,在生物多样性保护和畜牧业经济中发挥着关键作用[1—2]。草地地上生物量(Aboveground Biomass,AGB)是表征植被生长状况、固碳潜力的重要指标,可以有效衡量草地生产功能,是草地利用决策和资源管理的基础[3—4]。明晰草地AGB现状对草地资源管理和生态保护具有重要意义。

目前,草地AGB监测方法主要包括传统的地面测量和遥感监测。地面测量即直接收割法,通过割草、干燥、称重获得AGB,但需要耗费大量时间与人力,部分地区容易受地形和环境制约,不适合用于长时间序列、大区域尺度的生物量估算[5]。遥感监测以其宏观性、综合性等特点,克服了地面测量的不足,被视为大范围监测的有效工具[6]。早期的草地AGB遥感监测主要是构建植被指数与实测草地AGB的经验关系,并根据该关系反演区域尺度的草地AGB[7]。植被指数是通过对植被在不同波段部分的反射率数据进行组合和计算获得[8—9],可以有效反映植被生长状况和变化信息。常用的植被指数有归一化差值植被指数(normalized difference vegetation index,NDVI)[10]、增强型植被指数(enhanced vegetation index,EVI)[11]和土壤调整植被指数(soil adjusted vegetation index,SAVI)[12]等。但植被指数自身局限性可能会影响AGB估算准确性,例如在覆盖度高的草地上,NDVI敏感性降低,易出现饱和现象[13]。此外,草地AGB的变化受外部环境因素影响较大[14]。因此基于单因素的回归模型在估算草地生物量时存在一定的误差和不确定性。对于具有多种变量的模型,机器学习模型比传统的回归模型更适用[15—16]。机器学习模型可以整合草地AGB的多种影响因素,学习高度复杂的非线性映射,从而获得更好的拟合结果[17]。与传统的回归模型相比,梯度增强回归树是一种更强大的经验建模算法。梯度增强回归树是一种迭代决策树,利用加性模型和前向分布算法将若干个弱学习器组合成强学习器[18]。该模型可灵活处理各种数据,预测性能好,模型复杂度较低,泛化能力较强[19]。Yang等[20]探索了三种机器学习方法,分别是人工神经网络、多元自适应回归样条和梯度增强回归树来建立估算模型,结果表明梯度增强回归树是生成1 km空间分辨率的全球森林AGB分布图的最优模型,模型验证R2达到0.90。Yu等[18]对比了梯度增强回归树、随机森林和极端随机树三种模型在生成中国三江源地区1982—2018年草地AGB数据集方面的有效性,结果表明梯度增强回归树模型精度最高。以上研究表明梯度增强回归树模型在估算草地AGB方面具有很大潜力。

青藏高原是世界上最高的独特自然地理单元,有“地球第三极”之称,是我国重要的生态安全屏障,也是全球气候变化的敏感区[21—22]。草地是青藏高原的主要生态系统,由于青藏高原特殊的自然环境以及多年来气候变化和人类活动影响,草地生态系统十分脆弱,局部退化严重[1,17,23]。及时准确的估算青藏高原草地AGB能够为草地资源的可持续管理提供参考意义。然而青藏高原地形地貌复杂多样,气候差异显著[24],基于单一植被指数的回归模型无法准确反映青藏高原生态环境的影响。环境因素与草地AGB密切相关[25],地形影响土壤湿度、土壤有机质含量、接收的太阳辐射强度和温度,间接影响到物种生长和草地AGB分布格局[26]。降水气温等关键因素的变化将会影响物种的分解、蒸散和光合作用等生理生态过程从而影响草地AGB[27]。因此构建基于气象、地形和遥感因子的多因素模型是准确估算青藏高原草地AGB的可能途径。

为此,本研究将青藏高原草地作为研究区,借助Google Earth Engine(GEE)平台,探索遥感植被指数表征草地AGB信息的有效性,评估气象、地形信息对草地AGB估算的影响,构建基于气象、地形和遥感因子的青藏高原草地AGB估算模型—梯度增强回归树,实现青藏高原草地AGB空间分布制图,为准确估算青藏高原草地产草量提供重要依据。

1 材料与方法

1.1 研究区域

青藏高原地处26°00′12″—39°46′50″N,73°18′52″—104°46′59″E之间,面积约为257.24万km2,占中国陆地面积的26.8%[28],主要分布在西藏自治区、青海省、甘肃省、四川省、云南省和新疆维吾尔自治区。其地势西北高,东南低,平均海拔在4000 m以上,气候敏感独特,东南部温暖湿润,西北部寒冷干旱。由于气候变化和地形特征的共同作用,青藏高原植被类型复杂,其中草地为主要植被类型,草地空间分布数据(图1)来源于全球30 m地表覆盖数据产品GlobeLand30[29—30]。

图1 青藏高原草地AGB样点分布图Fig.1 Distribution of grassland AGB sample sites on the Qinghai-Tibet PlateauAGB:地上生物量 Aboveground Biomass

1.2 数据来源与处理

1.2.1草地AGB数据

本研究中青藏高原草地AGB野外样点数据共210个(图1),主要来源包括:(1)野外实测,共13个。野外调查时间在2019年8月,选择物种组成均匀、地势平整、面积足够的典型样地,记录经纬度、海拔和主要植被类型等信息,根据代表性原则在每个样地中均匀布设3个样方,样方大小为1 m×1 m,将样方中的地上部分沿地面刈割,除去附着的土壤等杂物后,烘干至恒重,该样地草地AGB即为3个样方的干重取平均值。(2)文献查阅。基于中国知网(https://www.cnki.net/)和Web of science(https://www.webofscience.com/),以“草地生物量”“生产力”“Aboveground biomass”为关键词检索,筛选样本点必须包含准确的采样时间、经纬度坐标和植被类型,共收集到符合标准且公开发表的学术论文34篇,实测样点197个。文献中样点采样时间为2004—2020年7—8月,采样通常选择地势平坦且具有空间代表性的草地样地,其中布设2—9个重复样方,样方大小为0.25 m×0.25 m—1 m×1 m,将样方内的草地地上部分齐地面刈割,放入烘箱烘干,称其干重。

1.2.2数据资料

遥感资料采用空间分辨率为30 m,时间分辨率为16 d的Landsat遥感卫星数据,包括Landsat 5、Landsat 7和Landsat 8。1984年3月发射的陆地卫星项目第五颗卫星,即Landsat 5卫星,携带TM传感器。1999年4月发射的Landsat 7卫星,与Landsat 5在空间分辨率和光谱特性方面基本一致。该卫星携带ETM+传感器,较TM传感器多了一个15 m空间分辨率的全色波段。2013年2月发射的Landsat 8 OLI传感器较ETM+传感器新增了深蓝和卷云波段,波段光谱范围划分更加精细。Landsat 5 TM、Landsat 7 ETM+和Landsat 8 OLI各个波段的具体信息见表1。实测草地AGB样点采样时间为2004—2020年7—8月,对应的影像选取时间分别为2004—2011(Landsat 5 TM)、2012(Landsat 7 ETM+)和2013—2020(Landsat 8 OLI),通过GEE平台获取Landsat影像,经去云、裁剪等处理后,根据样点采样时间提取相应的单波段反射率值和植被指数数据,由于去云后导致部分影像数据不完整,选取采样时间临近的影像代替。

表1 Landsat 5 TM、Landsat 7 ETM+和Landsat 8 OLI波段信息Table 1 Landsat 5 TM, Landsat 7 ETM+ and Landsat 8 OLI band information

辅助数据包括地形(高程、坡度、坡向)和气象(年平均气温、年降水量)资料。高程数据来自Shuttle Radar Topography Mission(SRTM),空间分辨率为30 m,由云计算平台GEE提供。坡度和坡向数据通过高程数据计算生成。气象数据来自国家科技基础条件平台—国家地球系统科学数据中心(http://www.geodata.cn)的逐月气温、降水数据,空间分辨率为1000 m[31]。下载2004—2020年逐月气象数据,对每年12个月份的气温数据取平均值得到年平均气温,12个月份的降水数据求和得到年降水量,将结果重采样至30 m空间分辨率,提取与各样点采样年份相对应的年平均气温和年降水量。

1.3 研究方法

1.3.1建模变量集构建

研究构建的建模变量集共包含20个变量,其中遥感数据包括单波段反射率Blue、Green、Red、NIR、SWIR1和SWIR2(表1),常用植被指数NDVI、EVI、GNDVI、NIRV、RVI、SAVI、DVI、OSAVI和MSAVI(表2)。地形数据包括高程(dem)、坡度(slope)和坡向(aspect)。气象数据包括年平均气温(temperature)和年降水量(precipitation)。

表2 本研究采用的植被指数Table 2 Vegetation indices used in this study

1.3.2变量筛选

特征变量间的多重共线性会影响模型预测精度,而且过多的特征变量会使模型更加复杂,计算效率低。为提高模型准确性,本研究选用递归特征消除法(RFE)优选特征变量,将所有变量输入到模型中进行训练,根据模型给出的特征变量重要性,删除最差的特征,用剩下的特征变量继续输入模型训练,再删除最差的特征,重复上述步骤,直至所有特征都被删除,根据删除特征的顺序,得到最优特征子集。

1.3.3模型构建

梯度增强回归树(Gradient Boosting Regression Tree,GBRT)[38]是基于集成学习中Boosting算法的一种改进。基本思想是以回归树作为弱学习器,每个学习器的训练都是基于上一个学习器的预测偏差结果,沿着减少残差的梯度方向构建新的学习器,从而使预测偏差不断减小,模型精度提高,经过多次迭代后,最终结果由所有学习器互相联系得到[19,39]。GBRT能够适应复杂的非线性关系,产生具有竞争力、高度鲁棒、可解释的程序,尤其适用于挖掘冗杂的数据[38]。

本研究利用GEE平台提供的GBRT算法估算草地AGB。为提高草地AGB估测准确性,分别建立遥感、遥感+气象、遥感+地形、遥感+气象+地形四种GBRT模型。基于实测AGB样点,提取对应采样时间的遥感、地形和气象信息,用于模型训练。并结合递归特征消除法筛选最优变量,比较每种GBRT模型变量筛选前后的精度。经反复训练选取拟合效果最好的模型反演青藏高原草地AGB空间分布。

1.3.4精度评价

选择样点数据的70%作为训练集,30%作为测试集,选用五折交叉验证进行准确性评估。同时,选取决定系数R2(Coefficient of Determination)和均方根误差RMSE(Root Mean Square Error)对模型拟合效果进行精度评价。

(1)

(2)

2 研究结果

2.1 实测草地AGB统计分析

青藏高原实测草地AGB统计分析结果表明(表3):210个实测样点的AGB值在草地生长季存在一定差异,最小值为2.95 g/m2、最大值为440.46 g/m2、平均值为135.09 g/m2,变异系数为0.70。其中,文献收集样点共197个,AGB平均值为134.89 g/m2,变异系数为0.71。野外实测样点共13个,AGB最大值为336.95 g/m2、最小值为31.13 g/m2、平均值为138.12 g/m2。

表3 草地AGB观测值的统计分析Table 3 Statistical analysis of grassland AGB observations

2.2 单变量回归分析

将单变量与实测草地AGB进行线性回归分析,对比分析结果如下(图2)。除坡度、坡向和年平均气温外,其余变量均通过显著性检验(P<0.01)。遥感因子与草地AGB间的决定系数R2在0.08—0.40范围内,其中,GNDVI与AGB之间相关性最高,R2为0.40。NDVI(R2=0.39)估算AGB的能力略低于GNDVI,NIR与AGB相关性最低,R2为0.08。气象和地形因子中,年降水量、高程与AGB之间相关性较高,R2分别为0.28和0.12。经比较,GNDVI与AGB的拟合效果较好,但仅能解释青藏高原40%的草地AGB变化。

图2 单变量与实测草地AGB的关系Fig.2 Relationship between single variables and measured grassland AGBNIR:近红外波段;SWIR1:短波红外1波段;SWIR2:短波红外2波段;GNDVI:绿色归一化植被指数;EVI:增强型植被指数;NDVI:归一化差值植被指数;DVI:差值植被指数;RVI:比值植被指数;NIRV:植被近红外反射率指数;SAVI:土壤调整植被指数;OSAVI:优化土壤调整植被指数;MSAVI:修改型土壤调节植被指数

2.3 模型性能评估

基于GEE平台构建的遥感、遥感+气象、遥感+地形、遥感+气象+地形四种GBRT模型性能如表4、表5所示。未经过RFE特征优选之前(表4),基于全部遥感因子构建的GBRT模型测试R2为0.55,RMSE为61.75 g/m2。较基于单一植被指数的回归模型精度有明显提高,可解释55%草地AGB变化情况。在遥感因子的基础上分别添加气象、地形因子,两种模型的测试R2分别为0.59和0.60。综合遥感、气象和地形数据后,模型测试精度达到0.62,RMSE为56.49 g/m2,模型性能得到了有效提升。

表4 未经过特征筛选的GBRT模型评估Table 4 Evaluation of GBRT models without feature filtering

表5 经过特征筛选的GBRT模型评估Table 5 Evaluation of GBRT models with feature filtering

基于RFE优选的特征变量构建的GBRT模型精度优于输入全部变量的模型。由表5可知,以GNDVI、Green、SWIR2、EVI、NDVI、DVI、Red、SWIR1、NIR、Blue、SAVI和MSAVI为输入变量的遥感模型测试R2为0.57,RMSE为60.27 g/m2。基于筛选后特征变量构建的遥感+气象、遥感+地形模型性能提高,测试R2分别为0.62和0.63,RMSE分别降低3.93 g/m2和4.49 g/m2。以GNDVI、precipitation、slope、temperature、dem、aspect、SWIR2、EVI、Blue、Green、NIR、SAVI和SWIR1为基础构建的遥感+气象+地形GBRT模型,测试集R2为0.66,RMSE为53.64 g/m2,比经过特征筛选后的遥感GBRT模型估测精度提高了0.09,RMSE降低了6.63 g/m2,模型训练集R2也从0.70提高至0.79。

经综合考虑气象、地形及遥感数据和通过RFE特征优选后,模型拟合精度逐步提高。以GNDVI、precipitation、slope、temperature、dem、aspect、SWIR2、EVI、Blue、Green、NIR、SAVI和SWIR1为基础构建的GBRT模型具有最佳性能,可解释青藏高原草地AGB 66%的变化(测试数据集R2=0.66,RMSE=53.64 g/m2,P<0.01)。由该模型估测的草地AGB与实测草地AGB的线性关系以及数据分布情况(图3),可得两者线性关系较好,但实测AGB中有些较高的值被估计的偏低,存在低估情况。从数据分布而言,实测AGB分布范围大于估测草地AGB,实测AGB数据四分位数范围在50—180 g/m2,估测AGB数据四分位数范围在80—160 g/m2。两者在中位数、平均数上具有相似的分布特征。

图3 实测AGB与最优GBRT模型估测AGB之间的线性关系和小提琴图Fig.3 Linear and violin plots between the measured AGB and the estimated AGB from the optimal GBRT model

2.4 青藏高原草地AGB估算

基于最优GBRT模型反演了青藏高原2010、2015和2020年草地AGB空间分布图(图4),草地AGB平均值范围为93.63—100.78 g/m2,分别是2010年94.58 g/m2,2015年93.63 g/m2,2020年100.78 g/m2。青藏高原草地AGB整体呈现东南高,西北低的空间格局。其中,西藏自治区阿里地区草地AGB较低,2010年基本小于60 g/m2,2015、2020年AGB有明显上升;青海省果洛藏族自治州南部、黄南藏族自治州和四川省西北部草地AGB相对较高,与实际空间分布基本一致。

图4 青藏高原2010、2015和2020年草地AGB空间分布图Fig.4 Spatial distribution of grassland AGB on the Qinghai-Tibet Plateau in 2010, 2015 and 2020

3 讨论

3.1 最优GBRT模型分析

研究[40—41]表明草地AGB受到多种环境因素影响,其中,气象、地形因素与草地生长状况密不可分。青藏高原气候敏感、地形险峻、生态环境复杂,为提高草地AGB估算准确度,基于遥感、气象和地形信息构建了GBRT模型,测试集R2达到0.66,表明构建多因素模型能够更好的解释草地AGB变化情况。这一研究结果与Liang等[42]、Tang等[16]一致,Liang等[42]基于地理位置(经度、纬度)和草地覆盖度构建了三江源地区多因素非线性回归模型,与基于NDVI的最佳单因素模型相比,模型RMSE降低了20%。Tang等[16]基于经度、植被指数、光照和土壤因子等环境变量构建随机森林模型估算黄河流域草地AGB,结果表明,该模型性能高于NDVI构建的单变量回归模型。

图5为最优GBRT模型中特征变量重要性排序,重要性得分越高,说明该特征对模型估算结果贡献越大。遥感因子特征重要性占比62.06%,其中,GNDVI最高,该指数使用绿波段代替了NDVI中的红波段,对叶绿素含量更加敏感,赵翊含等[43]研究表明由绿波段构建的植被指数对高寒草地AGB有重要作用。气象与地形因子特征重要性占比37.94%,年降水量特征重要性最高,说明年降水量与AGB密切相关。大量研究也表明[17,44]在青藏高原地区,降水比气温对草地AGB影响更大。气象、地形信息是影响草地AGB的重要环境因素,降水影响土壤含水量、空气湿度,温度影响植被的光合作用、蒸腾作用、呼吸作用[27],地形主要影响气温、湿度和光照度[26],以及影响土壤中物质、能量的再分配影响植被生长。因此,综合考虑了气象、地形和遥感多个因素的最优GBRT模型计算效率高,预测准确性强,在估算青藏高原草地AGB方面具有良好应用前景。

图5 最优GBRT模型中特征变量重要性Fig.5 Importance of feature variables in the optimal GBRT model

3.2 青藏高原草地AGB比较

本研究使用最优GBRT模型反演了三年青藏高原草地AGB,空间变化呈现自西北向东南方向增加的趋势,研究结果与Zhang等[41]、金哲人等[45]一致。AGB平均值范围为93.63—100.78 g/m2,略高于Zeng等[17](77.12 g/m2),与Zhang等[41](102.40 g/m2)最为接近。研究[17,41]表明过去十几年青藏高原草地AGB呈上升趋势,本研究中AGB平均值虽有较小波动,总体仍呈上升趋势。不同研究中存在AGB估算差异的原因可能是:数据来源不同,包括实测AGB数据、遥感影像数据等。Zeng等[17]研究所用遥感数据是空间分辨率为250 m的MODIS数据,与本研究中数据空间分辨率不同。机器学习模型效果易受样点数量和分布情况影响[46]。而且,不同研究所用模型不同,也会导致估算差异[47]。

本研究中最优GBRT模型具有较高的精度,但机器学习模型性能依赖于训练数据。由于交通不便和海拔限制,青藏高原西北部样点较少,这可能会影响该区域草地AGB估算的准确性;此外,AGB样点大多收集于不同文献,尽管操作人员在测量时严格遵守了相关规范,但不同研究的测量方法、样本选择可能存在不同,实测数据中的不确定性可能会导致估算误差。未来可继续扩大样点数量,尝试不同的深度学习算法及特征优选方法,提升草地AGB遥感监测和建模的准确性。

4 结论

本研究利用GEE平台和机器学习算法,结合详实的2004—2020年期间青藏高原实测草地AGB数据,评估了基于不同特征变量构建的GBRT模型的性能,估算了青藏高原草地AGB,得到主要结论如下:

(1)基于单一遥感因子的线性回归模型性能较差,仅能解释草地AGB变化的8%—40%。基于全部遥感因子构建的GBRT模型测试数据集R2为0.55。

(2)基于遥感+气象、遥感+地形的GBRT模型拟合效果较好,经RFE特征筛选后,测试数据集R2分别为0.62和0.63,对草地AGB的解释能力有所增强。

(3)在青藏高原地区,基于气象、地形和遥感信息构建GBRT模型能够有效提高草地AGB估算精度,模型验证表明,R2=0.66,RMSE=53.64 g/m2,P<0.01,估测值与实测值吻合性较好,适用于青藏高原草地AGB模拟。

(4)利用最优GBRT模型反演青藏高原草地AGB分布图,呈现东南高,西北低的空间格局。