基于无人机高光谱特征参数和株高估算马铃薯地上生物量

刘 杨，冯海宽，黄 珏，杨福芹，吴智超，孙 乾，杨贵军

1. 农业部农业遥感机理与定量遥感重点实验室，北京农业信息技术研究中心，北京 100097 2. 山东科技大学测绘科学与工程学院，山东 青岛 266590 3. 国家农业信息化工程技术研究中心，北京 100097 4. 北京市农业物联网工程技术研究中心，北京 100097 5. 河南工程学院土木工程学院，河南 郑州 451191

引 言

地上生物量(AGB)是作物长势监测和产量预测的重要指标，其周期性变化能够很好地反映作物的生长情况和营养状况[1-3]。因此快速、准确、高效率地获取AGB信息对农业管理实施决策方案起到重要的作用。传统的测量手段虽能高精度获取作物参数信息，但田间调查往往对作物造成破坏，不仅消耗人力和物力，且应用范围有限，难以大面积的推广应用[4]。

遥感技术凭借其大尺度、实时动态和高通量等优势在农业监测方面得到广泛的关注[5]。无人机遥感平台相比航空航天遥感，具有受限条件少且机动性强，获取的影像分辨率更高，运营成本更低等优势，故而在农田管理方面得到重用[6-7]。目前，无人机搭载的多为数码和多光谱相机，其传感器波段少，而高光谱传感器所含波段数量较多，光谱分辨率较高，能够获取作物冠层大量的光谱信息，从而可以有效地估算作物长势，为精准农业管理提供了可靠手段。Yue等[8]利用无人机高光谱影像，基于植被指数以及植被指数融入株高数据采用PLSR对冬小麦每个生育期的AGB进行反演，结果表明融入株高数据的反演模型精度更高。尼加提·卡斯木等[9]利用无人机高光谱影像数据，基于16个光谱特征参数与春小麦叶绿素含量进行相关性分析，挑选出以820～940 nm的最大反射率及反射率总和利用PLSR方法进行建模的精度最高，R2达到了0.80。秦占飞等[10]利用无人机高光谱影像，基于高光谱影像指数和偏最小二乘方法估测引黄灌区水稻叶片全氮含量，结果表明以波段738和522 nm处一阶微分光谱构建的比值光谱指数(ratio spectral index，RSI)进行构建的模型精度最高。Tao等[11]利用无人机高光谱影像数据，基于植被指数、红边参数以及二者结合采用PLSR和SWR估算冬小麦AGB，结果表明植被指数结合红边参数能够提高模型反演精度。刘畅等[12]利用无人机高光谱影像，基于植被指数和纹理特征相结合形成“图-谱”融合指标采用逐步回归法估算冬小麦AGB的模型精度高于单一植被指数或者单一纹理特征所构建的AGB模型。Li等[13]利用无人机高光谱影像，基于植被指数和植被指数中融入株高采用PLSR和RF对马铃薯AGB进行反演，得出2种方法都以融入株高的模型进行反演精度较高。以上研究对于光谱细节变换及以多个特征参数来进行建模相对较少，这可能造成高光谱数据无法充分利用，进而模型精度受到一定的制约。本文以马铃薯为研究对象，利用无人机高光谱冠层数据，提取马铃薯株高，并结合地上实测生物量分析马铃薯不同生育期的光谱特征，利用高光谱特征参数(HCPs)、最优绿边参数(OGEPs)和基于DSM提取的马铃薯株高(Hdsm)3种参数以PLSR和RF进行建模，以期深度挖掘光谱信息并进一步提高马铃薯地上生物量估算精度，为精准农业定量化遥感提供技术手段。

1 实验部分

1.1 研究区域与试验设计

于2019年3月到7月在北京市昌平区小汤山镇国家精准农业研究示范基地的马铃薯试验田开展试验。该处位于北纬40°10′34″，东经116°26′39″，平均海拔高度为36 m，气候类型为暖温带半湿润大陆性季风气候，年均降水量约为640 mm，年均温度约为10 ℃，年均无霜期180 d。如图1所示，试验区共设密度试验(N区)、氮素试验(S区)、钾肥试验(K区)3个试验区，采用完全随机试验设计，N区和S区均设2个试验品种，分别为中薯5(P1)和中薯3(P2)，均为早熟品种，K区设1个试验品种为中薯3。密度试验设3个水平： 60 000株·hm-2(T1)、72 000株·hm-2(T2)、84 000株·hm-2(T3)，处理6个，每个处理重复3次，共18个试验小区； 氮素试验设4个水平： 0 kg尿素·hm-2(N0)、244.65 kg尿素·hm-2(N1)、489.15 kg尿素·hm-2(N2,正常处理，15 kg纯氮)、733.5 kg尿素·hm-2(N3)，处理8个，每个处理重复3次，共24个试验小区； 钾肥试验设2个水平： 0 kg钾肥·hm-2(K0)、970.5 kg钾肥·hm-2(K1，N区和S区均为K1处理)、1 941 kg钾肥·hm-2(K2)，处理2个，每个处理重复3次，共6个试验小区； 小区总计48个，小区面积为6.5 m×5 m。为了更好的获取试验田的DSM，本试验在试验小区周围均匀布控11个GCP(k1—k11由0.3 m×0.3 m的木板和埋于地下的木桩组成，其上有黑白标志的塑料板，目的是准确确定木板的中心位置)并用差分GPS测定其三维空间位置。

图1 马铃薯试验位置和试验设计Fig.1 The experimental area and design of potatoes

1.2 地面数据获取

分别获取马铃薯现蕾期(2019年5月13日)、块茎形成期(2019年5月28日)、块茎增长期(2019年6月10日)、淀粉积累期(2019年6月20日)、成熟期(2019年7月03日)5个关键时期的株高和生物量数据。其中马铃薯株高的观测方法为： 在每个小区中随机选取具有代表性的4颗植株，分别用直尺测量从茎基到叶顶端的距离，取其均值作为该小区的实测株高。马铃薯地上生物量通过收获法获取，在每个小区中选取代表性的3颗植株，将其茎叶分离，随后用清水洗净，105 ℃杀青，80 ℃烘干48 h以上，直到质量恒定再进行称重。将植株茎、叶的干质量总和进行相加，最后根据种植密度和样本干重换算成单位面积的马铃薯地上生物量[14]。

1.3 无人机高光谱数据获取及处理

利用无人机搭载的德国Cubert公司生产的UHD185机载成像光谱仪分别获取裸土期(2019年4月20日)至成熟期的高光谱影像，该传感器大小为195 nm×67 nm×60 nm，质量470 g，波谱范围为450～950 nm，共有125个光谱通道，采样间隔4 nm，光谱分辨率8 nm，数字分辨率12位。选择晴朗无云的天气进行拍摄，此时太阳光照强度稳定。飞行时间10:30—14:00，无人机的飞行高度为50 m。

无人机获取的高光谱数据处理主要包括几何校正、影像拼接、影像融合和光谱提取4部分组成。首先将航带内的灰度影像和地面11个控制点导入到俄罗斯Agisoft LLC公司的PhotoScan软件中，利用GCP的三维坐标信息进行几何校正，进而生成密集点云完成影像拼接； 其次通过Cubert-Pilot软件，进行高光谱与灰度影像融合； 最后在Arcmap10.2软件中,绘制马铃薯每个小区矢量数据，基于IDL语言统计每个小区平均光谱作为马铃薯冠层光谱，得到高光谱反射率数据。

1.4 植被指数和绿边参数选取

植被指数(vegetation indexs,VIs)就是将高光谱不同波段反射率进行组合，来消除或者降低背景环境对冠层光谱信息的影响。根据文献构建了SAVI*SR，所得效果较好，由于指数SRI和SR性质相似，且对于作物冠层光谱信息较为敏感，本文新构建了植被指数SAVI*SRI。可见-近红外形成的蓝边、绿边、黄边和红边等特征参数是区别植被与背景的特有参数，因此可以用来估测作物的长势参数[15]。本文以502～554 nm作为绿边区域，提取绿边最大反射率、绿边反射率总和、绿边位置、绿边面积、绿边振幅、最小振幅和绿边振幅/最小振幅，将这7个绿边参数与挑选的12个植被指数一起作为模型因子估算马铃薯的AGB，选取的模型参数如表1所示。

表1 选取的模型参数Table 1 Modeling parameters selected

1.5 分析方法

马铃薯株高的提取，即利用试验田的灰度影像，结合GCP，借助PhotoScan软件生成试验田的DSM，通过将马铃薯不同生育期的DSMi(i=1,2,3,4,5)与播种后至出苗前的裸土时期之间的DSM0作差运算，可以得到相应生育期马铃薯的株高Hi。计算公式见式(1)

Hi=DSMi-DSM0(i=1,2,3,4,5)

(1)

式(1)中： DSM0，DSM1，DSM2，DSM3，DSM4和DSM5分别代表裸土期、现蕾期、块茎形成期、块茎增长期、淀粉积累期、成熟期的高程；H1，H2，H3，H4和H5分别代表各生育期的提取株高。

本文估算马铃薯不同生育期AGB时，所采用的方法为偏最小二乘(partial least square regression,PLSR)和随机森林(random forest,RF)。PLSR是利用了多元线性回归、典型相关分析和主成分分析结合为一体，可以提供一种多对多的线性回归建模方法，能够消除自变量之间的相关性，使得用较少数据来估测因变量。RF是一种机器学习方式，通过多个决策树建立分类器，以bootstrap取样方式对训练样本集进行抽样调查，最后以多个决策树对训练样本进行预测。本工作的输入数据为3种参数，输出结果为马铃薯的地上生物量。

1.6 精度评价

每个生育期分别挑选2/3样本数据(32个，重复1和重复2)作为建模集，1/3样本数据(16个，重复3)作为验证集以此来构建马铃薯生物量估算模型。采用决定系数(coefficient of determination,R2)、均方根误差(root mean square error,RMSE)、标准均方根误差(normalized root mean square error,NRMSE)来评价模型的精度。R2越接近于1，RMSE和NRMSE越低，估测模型的精度就越高。

2 结果与讨论

2.1 马铃薯株高的提取

根据DSM结合GCP信息，分别得到马铃薯现蕾期、块茎形成期、块茎增长期、淀粉积累期和成熟期的估测株高(Hdsm)，其结果分布如图2所示。5个生育期共得到240个马铃薯株高数据，将提取的株高与实测的马铃薯株高进行对比分析可知，基于DSM提取的马铃薯株高Hdsm和实测株高H，R2达到了0.84，NRMSE为15.67%，其值处于10%～20%之间，说明提取的株高Hdsm精度较高，对马铃薯生物量能够较好的估测。因此，利用基于DSM提取的马铃薯株高进行相关研究。

图2 马铃薯不同生育期的株高空间分布Fig.2 Spatial distribution of height at different growth stages of potato

2.2 相关性分析

2.2.1 一阶微分光谱

除了常见的波段组合形成的植被指数外，光谱微分等变换形式也可以形成光谱特征参数。一阶微分代表反射率的变化，即波长的斜率，可以去除部分线性背景噪声对作物冠层光谱的影响，高光谱含有波段数量较多，适合对其进行一阶微分处理。马铃薯5个生育期的一阶微分光谱与AGB的相关性如图3所示。现蕾期，在554～582，862～950和702～774 nm波段区间内AGB分别与一阶微分光谱呈极显著负、正相关，相关性最好的波段为554，906和742 nm，其相关系数值分别为-0.59，-0.65和0.65。由于与AGB相关的光谱反射率波段为可见光波段，因此，在554～582和702～774 nm区间内，选取波谷554 nm和波峰742 nm作为敏感特征参数估算马铃薯的AGB。块茎形成期，在550～578，782～810，838～862，866～950 nm和658～674，698～774 nm波段区间内AGB分别与一阶微分光谱呈极显著负、正相关，相关性最好的波段为558，802，850，946 nm和670，718 nm，其相关系数值分别为-0.66，-0.59，-0.60，-0.69和0.49，0.72。在其他波段处，还存在一些离散点呈极显著负相关。因此，在可见光550～578，658～674和698～774 nm区间内，选取波谷558 nm，波峰670nm和波峰718 nm作为敏感特征参数估算马铃薯的AGB。块茎增长期，在550～646，786～830，846～950 nm和518～530，686～778 nm波段区间内AGB分别与一阶微分光谱呈极显著负、正相关，相关性最好的波段为570，806，950 nm和526，742 nm，其相关系数值分别为-0.69，-0.61，-0.73和0.52，0.75。在其他波段处，还存在一些离散点呈极显著负相关。因此，在可见光518～530，550～646和686～778 nm区间内，选取波峰526，波谷570 nm和波峰742 nm作为敏感特征参数估算马铃薯的AGB。淀粉积累期，在558～626，634～666，814～826，858～950 nm和514～538，682～786 nm波段区间内AGB分别与一阶微分光谱呈极显著负、正相关，相关性最好的波段为570，642，818，922和518 nm和722 nm，其相关系数值分别为-0.68，-0.67，-0.57，-0.76和0.49和0.76。因此，在可见光514～538，682～786，558～626和634～666 nm区间内，选取波峰518 nm，波峰722 nm，波谷570 nm和波谷642 nm作为敏感特征参数估算马铃薯的AGB。成熟期，在522～542 nm范围内AGB与一阶微分光谱呈极显著负相关，相关性最好的波段为526 nm，相关系数为-0.44。在570～682 nm范围内，除了个别波段未达到极显著相关外，其余波段的一阶微分光谱与AGB呈极显著正相关，相关性最好的波段为622 nm，相关系数为0.59。

图3 马铃薯AGB与不同光谱参数的相关性Fig.3 Correlation analysis of different spectral parameters with potato AGB

2.2.2 植被指数和绿边参数

将挑选的12种植被指数和7个绿边参数，分别与马铃薯不同生育期的AGB进行相关性分析，结果见表2所示。根据表2可知，19种光谱参数与马铃薯5个生育期的AGB绝大多数达到极显著水平(p<0.01)，各生育期光谱参数与AGB的相关系数有所波动，差异性明显，整体上5个生育期的绿边参数与AGB的相关性要远弱于植被指数的相关性。现蕾期，所有植被指数均呈现显著相关(p<0.05)，相关性最高的植被指数为MSR，相关系数是0.675，本工作构建的植被指数SAVI*SRI的相关系数是0.669，相关性位于第3位； 绿边参数除了GEP，SDr和Dr未达到显著相关外，其余绿边参数均达到极显著相关，相关性最好的绿边参数为Rsum，相关系数是-0.586。块茎形成期，所有植被指数均达到极显著相关(p<0.01)，相关性最好的植被指数是PSRIII，相关系数为-0.727，新构建的植被指数的相关系数为0.672，仅次于极显著最高相关的PSRIII； 绿边参数除了GEP，SDr和Dr/Drmin未达到显著相关外，其余参数均达到显著相关，相关性最好的绿边参数为Drmin，相关系数是-0.545。块茎增长期，所有植被指数与AGB均达到极显著相关，相关性最高的植被指数为MSAVI，相关系数是0.730，新构建的植被指数相关系数为0.710，相关性位于第9位； 绿边参数除了GEP和Dr未达到显著相关外，其余参数均达到显著相关，相关性最高的参数为Rsum，相关系数是-0.668。淀粉积累期，所有植被指数与AGB均达到极显著

表2 现蕾期、块茎形成期、块茎增长期、淀粉积累期、成熟期的相关性Table 2 Correlation of bud period,tuber form period,tuber grow period,starch store and maturity stages

相关，相关性最高的是SPVI，相关系数为0.756，本文构建的植被指数相关系数为0.735，相关性位于第6位； 绿边参数除了GEP未达到显著相关外，其余参数均达到极显著相关，相关性最高的绿边参数为Rsum，相关系数是-0.644。成熟期，植被指数与AGB的相关性在整个生育期最低，除了MSAVI，SRI，SAVI和MSR未达到显著相关外，其余植被指数均达到显著相关，相关性最高的是本文所构建的SAVI*SRI，相关系数为0.529； 绿边参数只有Rmax，SDr和Dr与AGB达到显著相关，相关性最高的是SDr，相关系数为-0.447。

2.2.3 敏感波段及光谱参数确定

为了尽可能避免高光谱特征参数之间的共线性，将每个生育期的一阶微分光谱、VIs和GEPs与马铃薯AGB相关系数绝对值按从大到小依次排列，筛选出每个生育期与AGB达到极显著水平的前7个高光谱特征参数和最优绿边参数，结果见表3所示。由表可知，构建的SAVI*SRI除了在块茎增长期未能作为模型因子外，在其他生育期均参与模型构建，表明本文构建的植被指数与AGB相关性较好。相比不同生育期，成熟期的前7个高光谱特征参数和绿边参数与AGB的相关性明显弱于前4生育期； 且每个生育期通过一阶微分光谱分别得到可见光范围内的较优微分光谱，说明微分光谱技术适合高光谱数据的处理。

表3 敏感光谱参数的相关性Table 3 Correlation of sensitive spectral parameters

2.3 马铃薯生物量模型建立与验证

每个生育期以马铃薯AGB为因变量，通过PLSR和RF两种方法利用HCPs(x1)、HCPs加OGEPs(x2)以及HCPs、OGEPs和Hdsm(x3)3种模型自变量分别构建马铃薯现蕾期、块茎形成期、块茎增长期、淀粉积累期和成熟期的AGB估算模型，并计算每种模型的精度评价指标，结果见表4所示。由表可知，每个生育期利用PLSR和RF两种方法以x1，x2和x3为自变量估算AGB，其效果逐渐变好，2种方法都以x3为因子的模型精度最高。从AGB模型的精度评价指标看，每个生育期利用两种方法构建的模型NRMSE值绝大多数在20%以内，说明构建的模型稳定。现蕾期到块茎增长期利用PLSR和RF以x3为变量构建的模型R2范围分别是0.61～0.75，0.56～0.71，RMSE范围是140.77～110.34和187.31～140.35 kg·hm-2，NRMSE范围是15.49%～11.02%，18.39%～16.84%，两种方法的验证R2也是不断增大，RMSE和NRMSE都是逐渐减小； 淀粉积累期到成熟期，R2范围分别是0.65～0.63和0.60～0.58，RMSE范围是112.32～233.96和193.40～206.79 kg·hm-2，NRMSE范围是13.93%～14.25%和17.98%～19.43%，验证的R2逐渐降低，RMSE和NRMSE都是逐渐增加。综合建模和验证集可知，每个生育期以x3为自变量利用PLSR估算AGB效果最好，其中块茎增长期构建的模型精度最高，建模最佳R2为0.75(RMSE=110.34 kg·hm-2，NRMSE=11.02%)，验证最佳R2为0.95(RMSE=86.50kg·hm-2，NRMSE=9.21%)，结果如图4所示。

基于无人机高光谱灰度影像结合GCP生成了试验田的DSM，提取了马铃薯株高Hdsm，并将其与5个生育期的240个实测株高H进行了对比，其R2达到0.86，表明基于DSM提取马铃薯株高Hdsm具有较高的精度，所得结果与他人提取作物株高结果基本一致[18]。提取的Hdsm相对于实测株高H整体偏低，这是因为实测马铃薯株高最高点冠层空间结构小，在进行三维重建时，可能被当作噪声清除，从而导致马铃薯顶层空间信息丢失，与实测株高H相比，基于DSM提取的株高Hdsm较小； 此外，无人机拍摄的马铃薯冠层影像中包括了大量的裸土像元，从而缺失了一部分马铃薯冠层结构信息，这也导致基于DSM提取的Hdsm偏低。因此，确保马铃薯冠层空间信息的重建精度，对提高基于DSM提取的Hdsm精度尤为关键。

分析每个生育期的高光谱特征参数和绿边参数与马铃薯AGB的相关性，并挑选出前7个HCPs和OGEP作为模型变量，经对比可知，成熟期的光谱参数与AGB的相关性绝对值普遍低于前4期，主要原因由于生长后期，地上部有机物不断向地下输送，外加连续多天大雨，造成地上部叶片迅速枯黄脱落，植被覆盖度明显低于前4期，冠层光谱信息的提取更多受到地面土壤的影响，造成冠层光谱信息不能够充分表达和生物量之间的联系，进而使光谱参数与AGB的相关性普遍偏低。每个生育期采用PLSR和RF方法构建马铃薯AGB估算模型，通过PLSR构建的估算模型和验证模型的R2相较于RF较高，RMSE和NRMSE相较于RF也都较低，表明PLSR方法拟合效果较好，得到的模型稳定性较强。

表4 PLSR和RF以不同变量的AGB估算与验证Table 4 AGB estimated and verified by different variables with PLSR and RF

图4 块茎增长期基于PLSR以x3估算AGB的建模与验证Fig.4 Modeling and verification of estimated AGB based on PLSR with x3 at tuber grow stage

而以RF进行建模和验证的精度较低，主要原因一是机器学习适用于较大数据集，而本研究用于建模数据32个，验证数据16个，都属于较小数据集。二是所选的8个参数存在一定的多重共线性，而RF对多重共线性不敏感。本工作得出各个生育期融入Hdsm能够提高模型精度，这与牛庆林等[18]将提取的株高和植被指数融合后估算玉米的LAI，精度明显提高，具有一致的结论。本工作每个生育期通过PLSR方法以x3为自变量估算马铃薯AGB精度最高，模型最为稳定，这与陶惠林等[19]关于估测冬小麦LAI的精度结果一致。本文仅用一年的马铃薯无人机高光谱影像数据构建模型，还需要进一步分析并利用多个年份和不同地点的马铃薯数据进行验证，以此得到一个更具有普适性的AGB估算模型。

3 结 论

(1)基于无人机高光谱灰度影像结合GCP生成试验田的DSM，基于DSM提取的株高Hdsm与实测株高H高度拟合(R2=0.84，RMSE=6.85 cm，NRMSE=15.67%)，说明利用DSM进行马铃薯株高提取方法可行。

(2)现蕾期，一阶微分光谱与AGB极显著相关的波段为554和742 nm； 块茎形成期，与AGB极显著相关的波段为558，670和718 nm； 块茎增长期，与AGB极显著相关的波段为526，570和742 nm； 淀粉积累期，与AGB极显著相关的波段为518，570，642和722 nm； 成熟期，与AGB极显著相关的波段为526和622 nm。

(3)每个生育期得到的最优绿边参数不完全相同，现蕾期、块茎增长期和淀粉积累期OGEPs为Rsum，块茎形成期和成熟期OGEPs分别为Drmin和SDr。

(4)每个生育期，与仅使用HCPs估算AGB相比，使用HCPs加入OGEPs，HCPs加入OGEPs和Hdsm可以提高AGB估算精度，且以后者为自变量提高精度的幅度更大。每个生育期利用PLSR和RF以前7个HCPs以及OGEPs和Hdsm估算AGB的建模和验证R2在块茎增长期表现效果最好，其中PLSR方法得到的模型效果优于RF模型。