苗春丽,伏 帅,刘 洁,高金龙,高宏元,包旭莹,冯琦胜,梁天刚,贺金生,2,钱大文
(1. 兰州大学草地农业科技学院 / 草地农业生态系统国家重点实验室 / 兰州大学农业农村部牧草创新重点实验室 /兰州大学草地农业教育工程研究中心, 甘肃 兰州 730020;2. 北京大学城市与环境学院, 北京 100871;3. 中国科学院西北高原生物研究所, 青海 西宁 810008)
天然草地牧草作为重要的生产资料,其营养价值及品质的好坏直接影响家畜的生长发育以及畜产品的质量和经济价值,对草牧业发展具有非常重要的意义[1]。然而,传统的多光谱遥感因较少的光谱通道及不连续的宽波段在天然草地地上生物量(aboveground biomass, AGB)估测方面有一定的局限性,通常存在空间分布率低、稳定性差、精度低等诸多问题。目前,国内外更多的研究主要使用新一代高分辨率Sentinel-2、Landsat8 和MODIS 等[2-3]多光谱卫星资料在估算AGB,这些卫星数据具有适合监测植被生长状况的红光和近红外等波段。但是,由于受到有限的宽波段的影响,依然存在反演精度不足、误差偏大等问题[4]。
高光谱遥感的光谱分辨率高、波段连续性强、光谱信息量大,具有传统遥感技术无法比拟的优势。高光谱成像是一种获取大量空间位置光谱信息的新技术,使任何一个光谱波长都可以用来制作可识别的图像。高光谱传感器的窄波段通道能够探测地物更加细微的光谱特征,而多光谱传感器的宽波段通道通常难以捕捉到这些特征。目前,高光谱遥感在地质勘探、环境资源监测、生态系统调查等领域发挥着愈加重要的作用。就农业领域而言,国内外已将高光谱遥感技术广泛应用于作物生物生长参量反演、农业生产及农情信息监测等方面[5-9]。草地高光谱遥感在近几年的研究中也得到了广泛的应用,包括草地营养状况、地上生物量和覆盖度的监测、物种识别、健康评价[10-17]。乌兰吐雅等[18]对内蒙古锡林郭勒草原生物量地物高光谱遥感的研究的结果表明:在选用的8 种常用的植被指数中,大气阻抗植被指数ARVI (R945,R620,R430)与生物之间的相关性最强,相关系数为0.874;张亦然等[4]在内蒙古科尔沁通过极限梯度提升(XGboost)算法分别用原始光谱、一阶微分光谱、二阶微分光谱植被指数构建草甸地上生物量模拟估算模型,结果表明原始光谱计算的植被指数模型估算效果最佳,均方根误差(root mean square error, RMSE)为1 402.6 kg·hm-2;张凯等[19]对甘南草地地上生物量地物高光谱遥感估算研究结果显示:在17 个特征变量中,以特征参数光谱反射率723 nm 处一阶微分值(D723)为变量的对数模型是最佳模型,相关系数达到0.896;高宏元等[20]在海北试验站利用地物光谱数据和随机森林算法研究草地AGB 估测模型的结果表明,在4 种特征变量分别构建的草地AGB 估测模型中,基于NDBleaf、OSAVI、TCARI、MTCI、PRI、SIPI 和VARIg7 种植被指数构建的RF 模型精度最高,测试集R2= 0.70,RMSE = 557.87 kg·hm-2。上述光谱变量、VIs 和吸收特征均对草地地上生物量具有一定的敏感性,对监测天然草地地上生物量的空间变异和分布特征有借鉴意义。但是,以上研究都是基于地物高光谱数据构建模型,监测范围有限,而本研究借助UAV 成像高光谱系统在光谱和空间分辨率方面的优势,建立的生物量模型,可对农牧业生产和精细化管理起到很好的支撑作用。
高光谱遥感技术监测草地生物量最常用的方法包括MLR 和PLSR 方法[21]。相较于传统回归方法,随机森林(random forest, RF)机器学习具有避免过度拟合和多重共线性、调优简单、可并行处理、计算速度快等优点,并且能够处理离散数据,也能够处理连续数据,可以增强模型的泛化和抗噪能力。高光谱成像数据具有较高的维数,且包含大量冗余数据,需要去除无关数据,提取特征波长。热力图通过高光谱波段两两组合构建所有植被指数,打破了传统构建植被指数方法的不连续性和不完整性。最小绝对压缩变量筛选方法(least absolute shrinkage and selection operator,LASSO)能够最大程度上对高光谱波段及植被指数进行降维。以往研究,往往是对敏感波段范围,进行对应的卫星光谱响应函数计算或者均值化处理,造成敏感信息的缺失,而本研究通过热力图分析和LASSO 筛选,可精准定位到敏感波段。
基于以上因素的考虑,本研究结合地面实测数据和无人机高光谱遥感技术,利用UAV 搭载的成像高光谱仪对海北试验站天然草地冠层进行光谱成像和分析,筛选光谱敏感波段,结合机器学习方法构建天然草地地上生物量估测模型,为高寒天然草地地上生物量的高精度监测提供理论基础。
本研究的试验区位于青海省海北高寒草地生态系统国家野外科学观测研究站,地理位置为37°37′N,101°19′ E,土壤类型有高山草甸土、高山灌丛草甸土和沼泽土,发育差,砾质性,土层深度60 cm 左右,母质为黄土,其下为洪冲积物。土壤呈现有机质及全量养分丰富而速效养分贫乏的特点。植被类型为青藏高原典型的地带性植被,以金露梅(Potentilla fruticosa)为建群种的高寒灌丛草甸分布于山地阴坡、山麓及河谷低地;以嵩草属植物为建群种的高寒嵩草草甸分布在山地阳坡和滩地;以藏嵩草(Kobresia tibetica)为建群种的沼泽草甸分布于河滩地。群落结构简单、种类组成较少,植物生长期短、生产力较低。家畜主要是藏系绵羊和牦牛。
如图1 所示,试验区包括A、B 两块放牧地,每块放牧地有15 个小区(每个小区面积0.2 hm2),放牧家畜为牦牛,放牧方式为轮牧,放牧时间为每年7 月-9 月。
图1 研究区样地空间分布Figure 1 Spatial distribution of sampling sites of the study area
外业调查工作选在高寒草甸的生长季进行,于2019 年5 月下旬开始,包括9 次监测数据,其中2019 年分别在5 月-9 月的每月下旬开展了5 次观测。调查样本135 个,其中试验区A 有75 个,试验区B 有60 个。2020 年分别在7 月、9 月、10 月的下旬开展了3 次调查,观测样本30 个,均在试验区A。考虑到10 月份研究区高寒草甸已进入枯草期,所以在后续的研究中没有使用该月的观测数据。本项研究采用的7 次野外调查观测样本总计165 个。
在晴天、无云的状态下使用大疆经纬M600(http://www.dji.com)搭载双利合谱 Gaiasky-mini2-VN成像光谱仪(光谱范围为388~1 000 nm,光谱分辨率为 3.5 nm)在航高30 m (高光谱影像地面分辨率为 6.84 cm)开展外业高光谱遥感图像采集。A、B两块放牧地设置大小为40 m × 50 m 的小区各15个,每个小区分别选取0.5 m × 0.5 m 的3 个样方[22]进行观测,主要记录高程、经纬度、草地盖度、草群高度、优势种、物种数、鲜重等指标。为了减少鲜重数据的随机误差,避免局部采样位置空间异质性造成的影响,每个样地内采集的所有样方的草样均在实验室60 ℃下烘干48 h,将每个样地的3 个样方的草样干重取平均值,代表其地上部分的生物量。
使用仪器设备自带的Spec VIEW 高光谱图像采集及数据预处理软件,对采集的成像高光谱数据进行镜头校正、反射率校正和降噪处理以排除环境因素的干扰。将预处理后的影像数据导入ENVI 5.3 软件中转成TIF 格式文件,在ArcGIS 10.2 中对高光谱图像定义地图投影。
1.4.1 植被指数的热力图筛选
利用MATLAB 2021 软件,采用热力图方法[23],分别对常用的与草地AGB 有一定相关性的NDVI、SR、GWDRVI、TSAVI 4 种植被指数算法(表1)进行热力图形式的筛选分析。由于热力图表示的是所有波段进行两两组合,即每种植被指数的波段组合有30 625 种组合形式。为了从众多组合形式中筛选出对反演草地AGB 贡献较大的植被指数,每一种植被指数选择其R2排序前300 的进行筛选分析,并在MATLAB 2021 软件制图。
表1 植被指数列表Table 1 List of vegetation index
1.4.2 草地地上生物量反演模型的变量筛选
采用最小绝对压缩变量筛选方法(least absolute shrinkage and selection operator,LASSO),在R Studio中对165 个样本的176 个原始波段和1 200 个热力图植被指数进行筛选,为了获得稳定的结果,设置重复次数为200 次。LASSO 是最初由Tibshirani[28]提出的一种用于处理多重共线数据的有偏估计方法。该方法通过构造惩罚函数得到一个精细化模型,可以在不受显著影响的情况下将系数压缩到零,进行变量选择和变量估计。在回归系数绝对值之和小于一个常数值的约束条件下,最小化残差平方和,从而生成严格等于零的回归系数,从而得到一个可以解释的模型。本研究使用R Studio 软件中的LASSO 程序进行天然草地AGB 重要波段的筛选。
1.5.1 模型构建
采用RF 方法[29]构建研究区高寒草甸草地AGB反演模型。RF 算法包括ntree (回归树数量)、mtry(每个节点上测试的预测器的数量)和nodesize (树的终端节点的最大最小)3 个主要变量。ntree 是基于样本数所建立的决策树数量,ntree 越大,模拟结果就越稳定,但会导致计算量增加,本研究设置ntree 为1 000;mtry 是随机特征的数量,其默认值是输入变量的平方根,从 2~20 进行实验,间隔为 1,本研究采用默认值设置。RF算法是通过R 软件中自带的 “randomForrst”数据包实现。
1.5.2 精度分析
本研究采用十折交叉验证[30]方法评价模型性能。在建立模型之前,将所有自变量及其对应的因变量以等样本数的方式分为 10 组,每次选取其中一组作为测试数据集用来验证模型的估测能力,其余 9 组作为训练集用来构建模型。每次进行交叉验证,都会计算出预测值与实测值之间的决定系数(R2)和RMSE,重复该过程 10 次,最终用R2、RMSE来评价各个模型的精度。R2表示自变量与因变量的相关性程度,RMSE 用于衡量预测值与真实值之间的误差。R2越大,对应的RMSE 越小,则可以视为该模型的模拟结果精度越高。R2和 RMSE 的计算公式如下:
式中:xi和yi分别为模拟和实测的地上生物量数值,y¯为实测地上生物量的平均值,n为样本数。R2越大,说明模型估测的可信度越高,RMSE 越小,说明模型的拟合度越好。
研究区外业调查样点草地生物量的统计结果如表2 所列。5 月份2 个样区的高寒草甸地上生物量平均值为717.41 kg·hm-2,标准差为240.80 kg·hm-2,6 月份地上生物量平均值为1 609.00 kg·hm-2,7 月份地上生物量较高,平均值为2 602.50 kg·hm-2,8 月草地地上生物量达到最大值,平均值为3 025.70 kg·hm-2。总体来看,5 月-8 月高寒草甸地上生物量逐渐增大,其中8 月份达到最高值,9 月份开始下降。
表2 2019-2020 年观测样点草地生物量统计分析结果(n = 165)Table 2 Statistical analysis table of grassland biomass description (n = 165)kg·hm-2
图2 是成像高光谱数据的原始光谱在5 月-9 月的平均反射率曲线,光谱波段的波长范围介于394~996 nm,反射率在0.01~0.43,反射率最大范围在780~996 nm 的近红外光谱,在394~760 nm的可见光波段内天然草地的反射率较低,且存在两个吸收谷和一个反射峰,即450~490 nm 的蓝光吸收谷、490~580 nm 的绿光反射峰和640~680 nm 的红光吸收谷。这主要是由于天然草地叶片内叶绿素a 和叶绿素 b 在蓝光和红光附近的强吸收和绿光附近的强反射特性所致。在植被生长初期的5 月份,草地盖度较低,裸露土壤占比大,光谱曲线受地表裸土干扰,在可见光区域的波峰波谷不明显,在近红外部分的反射率较高。9 月份草地出现部分枯黄现象,所以反射率最低。6 月、7 月和8 月草地植被生长旺盛,叶绿素含量高,在近红外平台反射率较高,其中7 月份的反射率最高。
图2 高寒草甸成像高光谱反射率曲线Figure 2 Alpine meadow imaging hyperspectral reflectance curve
由于绿光、红光和近红外波段对天然草地地上生物量比较敏感,其光谱反射率差异大,利用这些波段构建的植被指数可以反映草地地上生物量的变化状况。广泛应用的NDVI、SR 等植被指数正是采用绿色植被在红光波段或绿光波段有强吸收和近红外波段具有强反射的特点构建的一些植被指数,此外也有采用绿波段与近红外波段构建的一些植被指数(如TSAVI 等),反演草地地上生物量的时空分布格局。
热力图(图3)可以直观反映由某些波段构建的某种形式的植被指数与草地地上生物量之间的线性相关关系。波段范围介于394~996 nm 的热力图反映出与草地AGB 线性回归决定系数R2位于前300的TSAVI、NDVI、SR 和GWDRVI 4 种植被指数的筛选结果。从图3 可以看出,TSAVI 的高光谱敏感波段位于热力图横轴821~996 nm 和纵轴503~530 nm;NDVI 与草地AGB 线性回归决定系数R2位于前300 的高光谱敏感波段分布范围介于热力图横轴的536~606 nm 和835~919 nm 两个区域,以及纵 轴 的474 ~533 nm 和513 ~523 nm 两 个 区 域;SR 的高光谱敏感波段位于热力图横轴的536~610 nm、779 nm 和810~996 nm 3 个区域,以及纵轴的471~533 nm、743 nm 和712~733 nm 3 个区域;相应的GWDRVI 的高光谱敏感波分布范围在热力图横轴的536~610 nm 和810~996 nm 两个区域,以及纵轴的471~533 nm 和712~736 nm 两个区域。由此可见,热力图中绿光波段、红光波段和近红外波段构成的植被指数与草地地上生物量具有显著的线性关系。
图3 植被指数热力图Figure 3 Thermal map of VIs
表3 和表4 分别表示基于176 个单波段和1 200个植被指数的1 376 个变量LASSO 回归分析结果和利用LASSO 筛选的变量和RF 方法构建的草地地上生物量估测模型十折交叉验证结果。由表4可以看出,在42 种变量组合构建的模型中,模型的训练集R2和验证集R2较接近,其中42 个模型的验证集R2大于 0.77。表明这42 种天然草地地上生物量估测模型均具有较高的精度。
表3 基于176 个单波段和1 200 个植被指数的1 376 个变量LASSO 回归分析结果Table 3 Results of LASSO regression analysis of 1 376 variables including 176 single bands and 1 200 vegetation indices
续表3Table 3 (Continued)
根据精度及稳定性,研究区草地AGB 最优模型为L42,在42 个RF 模型中验证集R2大于0.80 且RMSE 低 于490 kg·hm-2的 模 型 有L42、L31 和L28(表4),其地上生物量实测值和预测值的线性回归系数决定较高,分别达0.79、0.79 和0.78,RMSE 分别为506、502 和510 kg·hm-2(图4)。涉及变量数量分别为6、21 和24 (表3)。综合考虑模型的精度和简约性,最优模型可确定为L42。该模型涉及原始光谱波段536 nm、比值植被指数SR46 (513 nm, 566 nm)、SR47 (510 nm, 573 nm)和SR79 (733 nm, 850 nm),转换的土壤调节植被指数TSAVI287 (803 nm,879 nm)和TSAVI298 (503 nm, 875 nm)这6 个变量。模 型验证集的R2达0.81,RMSE 为489.36 kg·hm-2,训练集R2为0.95,RMSE 为248.70 kg·hm-2。
图4 不同RF 估测模型的实测AGB 与估测AGB 的相关性Figure 4 Correlation between measured AGB and predicted AGB in different RF estimations
表4 基于机器学习模型的草地地上生物量估测精度评价Table 4 Accuracy evaluation of AGB estimation based on machine learning models
图5 是2020 年7 月7 日利用无人机成像高光谱遥感影像和L42 估测模型对青海省海北试验站样区A 高寒草甸草地AGB 反演结果,草地地上生物量大部分介于1 500~2 000 kg·hm-2。从该图可以看出,利用高精度的反演模型和空间分辨率6.84 cm的无人机成像高光谱遥感影像,可以精准反演高寒草甸草地的空间分布格局,也可以清晰反映试验区的围栏、试验小区之间的道路、放牧家畜和裸地等其他地物的细节。通过生物量反演图,可以直观地看出试验区A 生物量的空间分布的异质性,有无人机拍摄的影像基本吻合,说明本项研究的变量选择和模型构建效果良好。
图5 草地地上生物量的空间分布图Figure 5 Spatial distribution of grassland aboveground biomass
利用LASSO 方法,对176 个成像高光谱的单波段和1 200 个植被指数进行筛选的结果表明,与高寒草甸AGB 密切相关的变量包括1 个单波段和5个植被指数,分别是536 nm 和比值植被指数SR46(513 nm, 566 nm)、SR47 (510 nm, 573 nm)、SR79 (733 nm, 850 nm),转换的土壤调节植被指数TSAVI287(803 nm, 879 nm)和TSAVI298 (503 nm, 875 nm),敏感波段如图6 所示。其中,单波段为绿光,组合植被指数涉及的波段分别为绿光、红光和近红外波段,说明这些通道与地上生物量存在重要关系。本研究提取的敏感光谱波段和植被指数对生物量的估测能力较强,高光谱的窄波段和窄带光谱指数相比其他多光谱卫星更能捕捉生物量的空间变化。
图6 AGB 估测变量的敏感光谱波段构成Figure 6 Composition of sensitive spectral bands of AGB estimated variables
传统意义上的NDVI 通常被认为在红光和近红外的波段范围的反射率与草地地上生物量的相关性较好[28]。然而,本研究结果中发现除了近红外波段和红光波段之外,红光和绿光构建的植被指数与高寒草甸AGB 相关性更高。由于传统多光谱遥感因较少的光谱通道及不连续的宽波段在天然草地地上生物量估测方面具有局限性,而UAV 成像高光谱的多通道、窄波段和连续性在估测草地地上生物量方面具有明显优势,其精度更高,效果更好。另外,采用热力图分析敏感波段具有更多优点:一是可以将所有波段进行两两组合,直观地反映决定系数R2的高值区域;二是每一种形式的植被指数都有其一一对应的热力图R2高值区,而传统意义上的植被指数波段往往以一个固定的波段范围来计算其相应的植被指数,这不利于最适波段范围的选取。
目前,对于天然草地地上生物量监测的研究多基于多光谱卫星影像数据,其空间代表性受到限制。而多光谱卫星遥感技术受时空分辨率等因素的影响,难以实现天然草地地上生物量的高精度监测[7]。目前基于传统卫星遥感的遥感反演生物量模型,主要适用于大尺度宏观估测,而借助UAV成像高光谱系统在光谱和空间分辨率方面的优势,建立的生物量模型,可进一步反映出田间尺度家畜采食轨迹,优势物种空间分布等信息,可对农牧业生产和精细化管理起到很好的支撑作用,是开展高精度天然草地地上生物量的动态监测的重要手段[2]。
研究区的高寒草甸生物量空间分布图表明,机载成像高光谱成像技术在样地或者较小空间尺度上对生物量以及不同放牧强度引起的生物量空间异质性具有较强的捕捉能力。但是,由于本研究观测的地域范围小,因此,研究结果的普适性仍然需要进行验证,后续研究应考虑增加更多的样点,并结合卫星遥感数据构建更大区域尺度、更高精度的估测模型。
本研究以青海海北试验站天然草地为研究对象,利用无人机成像高光谱遥感数据和机器学习算法,结合海北天然草地地上生物量实测数据,对比分析了不同时期的高光谱遥感特征,构建了天然草地生物量的估测模型,对天然草地地上生物量进行了反演和评估。主要得出以下结论:
1) 8 月份海北试验站的高寒草甸地上生物量较高,其平均值为3 025.7 kg·hm-2,影响地上生物量的重要变量为窄通道的绿光波段(536 nm)和由绿光波段和近红外波段构建的比值植被指数SR46 (513 nm, 566 nm)、SR47 (510 nm, 573 nm)、SR79 (733 nm,850 nm)以及转换的土壤调节植被指数TSAVI287(803 nm, 879 nm)和TSAVI298 (503 nm, 875 nm)。
2) RF 模型在高寒草甸地上生物量(AGB)估测上有独特的优势,其验证集的R2在 0.77~0.81,RMSE 介于508.825~554.48 kg·hm-2,训练集的R2在0.95~0.96, RMSE 介于237.39~274.44 kg·hm-2。
3)无人机高光谱遥感技术具有高空间分辨率、高光谱分辨率等方面的优势,能够为高寒草甸地上生物量的快速、精准监测提供重要理论依据。研究结果表明基于机器学习方法构建的估测模型的估测效果良好,可以为天然草地地上生物量监测技术提供参考。