无人机观测时间对玉米冠层叶绿素密度估算的影响

周丽丽,冯海宽,聂臣巍,许晓斌,刘 媛,孟 麟,薛贝贝,明 博, 梁齐云,苏 涛,*,金秀良,*

(1.安徽理工大学 空间信息与测绘工程学院,安徽 淮南 232001; 2.中国农业科学院 作物科学研究所,北京 100081; 3.中国农业科学院国家南繁研究院,海南 三亚 572024; 4.国家农业信息化工程技术研究中心,北京 100097;5. 武汉大学 测绘与遥感信息工程国家重点实验室,湖北 武汉 430079)

玉米是重要的粮食作物和工业原料,准确监测玉米的生长状态对农业经济发展具有重要意义[1]。叶绿素是玉米维持光合作用,从而保证产量和籽粒品质的重要色素。及时、准确地监测冠层叶绿素密度(canopy chlorophyll density, CCD)有助于及时掌握作物的生长胁迫状态和营养状况,为农艺学家进行变量施肥提供科学指导[2]。近年来,遥感技术发展迅猛,无人机(unmanned aerial vehicle, UAV)遥感更是以其低成本、可控性高、灵活快速的优势广泛应用于精准农业中[3],可以实现对作物冠层结构参数或生理生化参数的精准估计[4]。

目前,基于无人机遥感技术估算玉米CCD的研究已经在国内外取得了丰富的成果[5-7]。郭松等[8]基于低空无人机对比分析了单因素回归模型、多元线性回归模型和狮群算法优化的核极限学习机模型等各类模型对抽雄期玉米叶绿素的预测能力,结果表明,机器学习算法的效果优于传统回归算法。Qiao等[9]利用无人机RGB数据,建立并讨论了基于人工神经网络的玉米叶绿素含量估算模型。但纵观国内外研究现状,大多数基于无人机遥感平台的研究都是在中午经验性地选择时间进行观测试验的[10-11]。近年来,随着精准农业的不断发展,精细化的试验设计对无人机的作业时间提出了更高的要求。尽管目前很多研究将无人机观测试验的时间严格设定在10：00—14：00[12],但在该时段不同时间获得的冠层反射率与叶绿素密度之间的关系依然可能受光源与传感器观测几何特性的影响[13],无人机飞行时间的随机性仍然可能会影响叶绿素密度的估算精度;然而,关于该时段内不同时间获取的冠层光谱特性如何变化,以及是否会对叶绿素含量的估算产生影响等,尚鲜有研究。

在不考虑多角度观测的情况下,无人机田间试验一般是垂直下视观测,而玉米在同一天几个小时内的CCD变化不会很大[14];因此,在同一天不同观测时间变化最大的通常是太阳天顶角和方位角。PROSAIL模型可以在其他参数条件都固定的情况下通过改变太阳天顶角和方位角来很好地模拟不同天顶角和方位角环境下的植被冠层反射率[15-16],这为本研究提供了模型参考和理论依据。舒方等[17]基于PROSAIL模型分析了归一化植被指数(normalized difference vegetation index, NDVI)对叶绿素含量、叶面积指数(leaf area index, LAI)和太阳天顶角等参数的敏感性,并成功估算了冬小麦的叶绿素含量。何宇航[18]基于PROSAIL模型探究了多角度效应对估算小麦和茶树叶绿素含量的影响。通过计算同一天不同时间的太阳天顶角和方位角,可以模拟理想情况不同观测时间光谱反射率的变化;然而,仅通过模型模拟数据研究不同时间对CCD估算的影响存在一定的局限性,还需要进一步设计玉米不同处理下生长发育差异的田间试验。

综上,本研究围绕无人机不同观测时间对玉米CCD估算的影响展开研究,期望以理论模型和田间试验相结合的方式实现以下主要目标：(1)探究玉米冠层光谱随观测时间的变化规律;(2)挖掘冠层光谱与CCD的相关性在不同观测时间下的差异;(3)探索不同观测时间对玉米CCD估算精度的影响。

1 材料与方法

1.1 研究区与试验设计

本研究于内蒙古自治区通辽市科尔沁区钱家店镇全国玉米绿色高质高效行动示范基地(43°71′N,122°42′E,平均海拔166 m)开展玉米田间试验。研究区属温带半干旱大陆性季风气候,四季分明,雨热同期,气温适中,年平均气温6.1 ℃,年日照时数3 113 h,年均降水量385.1 mm,年均无霜期150 d,年均风速3.6 m·s-1,是种植玉米的黄金地带。

为了体现玉米的长势差异,于2022年设计了两个试验区,分别为氮调控区和氮梯度区,如图1所示,玉米品种均为迪卡159。氮调控区的种植密度为90 000株·hm-2,依据施肥方式的差异设置3类处理,分别为不施肥(ND0)、常规农户施肥(NDT)和施肥调控处理(ND1,ND2,…,ND27)。氮梯度区设置了2个种植密度水平[60 000株·hm-2(4K)和90 000株·hm-2(6K)]和10个施氮水平(N0、N3、N6、N9、N12、N15、N18、N21、N24、N36),二者相互组合,共形成20个处理。上述各处理分别对应于1个小区,每个小区的面积均为240 m2。

图1 研究区的地理位置与试验设计图Fig.1 Diagram of study area and experimental design

在氮调控区：各施肥调控处理统一按照450 kg·hm-2的量施用基肥,从玉米播种后50 d开始,按照试验设计方案(表1)施肥;常规农户施肥按900 kg·hm-2的量施基肥。基肥均使用掺混肥料(N质量分数15%,P2O5质量分数23%,K2O质量分数12%,云南美盛嘉吉化肥有限公司生产),追肥使用尿素(N质量分数46%,内蒙古天润化肥股份有限公司生产)。

表1 施肥方案

在氮梯度区：除N0处理不施肥,N3、N6处理不施基肥外,其他处理统一施基肥76.5 kg·hm-2;除N0外,各处理统一从玉米播种后50 d开始追肥,之后每隔10 d追肥一次,共追肥6次,每次等量追肥,N3～N36处理的追肥量分别为98.75、197.08、128.33、226.67、325.42、423.75、522.5、621.25、1 015.42 kg·hm-2。基肥和追肥所用肥料与氮调控区相同。

1.2 数据获取

1.2.1 无人机数据获取

分别在抽雄吐丝期(VT/R1,7月31日)、籽粒建成期(R2,8月11日)、乳熟期(R3,8月21日)和蜡熟期(R5,9月2日)采集无人机影像。因为上午和下午的太阳入射光亮度有限,所以选择10：00—14：00点进行无人机试验。分别在上述生育时期的10：00—10：59(10AM)、11：00—11：59(11AM)、13：00—13：59(1PM)和14：00—14：59(2PM)4个时间段的晴空无云条件下,采用经纬M300 RTK(深圳市大疆创新科技有限公司)无人机搭载RedEdge-MX(美国MicaSense)多光谱相机(波长400～900 nm)对玉米试验田进行飞行作业。在飞行前后使用自动校准的反射率面板对多光谱传感器进行辐射定标,飞行高度设置为70 m,速度为3 m·s-1,航向重叠和旁向重叠率均为70%。利用Agisoft PhotoScan 1.4.5软件(俄罗斯Agisoft)对无人机多光谱影像进行拼接、辐射校正等预处理,以获得正射反射率影像。

1.2.2 地面农学参数获取

在获取无人机数据的同时,同步进行地面农学参数获取,不同的是,每个生育时期仅获取一次数据。首先,对氮梯度区进行取样,每个小区选取3个位点,每个位点取3株玉米,共计9株;对氮调控区进行取样,每个小区取3个位点,每个位点取2株,共计6株。然后,带回实验室进行茎、叶分离,测量并计算每片叶子的叶面积(利用叶长、叶宽计算展开叶和未展开叶的叶面积时,比例系数分别取0.75和0.50)[19]。用计算得到的叶片面积与种植密度作比,计算LAI。同时,使用打孔器在每片叶片取5～6个完整的圆片(直径6 mm),称量质量后放入试管中,加入10～12 mL无水乙醇浸泡至叶片变白。使用分光光度计测量每个样本在663、645 nm处的吸光值(D663、D645),计算叶绿素含量(8.04D663+20.29D645,单位为mg·L-1),折算叶片叶绿素密度(LCD,单位为μg·cm-2)[20]。最后,对每个小区所有样本的叶片叶绿素密度和叶面积指数求平均值,并计算相应的冠层叶绿素密度(LCD与LAI之积)。

1.3 试验模型

1.3.1 随机森林模型

随机森林(random forest, RF)模型是一种基于决策树的非参数回归方法,通过建立多棵决策树对样本进行训练和预测,可以处理复杂的非线性问题,且高效稳定,在参数反演方面比其他方法更有优势[21-23]。本研究利用RF构建不同观测时间的CCD估算模型,以探索不同观测时间对CCD估算的影响。

RF包含3个重要的超参数：每个决策树变量的个数(mtry)、决策树的数量(ntree)和决策树节点最小个数(nodesize)。通过多次试验,本研究将上述超参数分别设置为5、800和1。由于本研究侧重分析不同观测时间对模型估算精度的影响,因此选择同样的超参数组合构建玉米CCD估算模型即可满足研究目标的需要。

1.3.2 PROSAIL模型

PROSAIL模型是由PROSPECT模型和SAIL模型耦合而成的冠层辐射传输模型[24]。该模型可以根据叶片生理生化参数、冠层结构参数和太阳角度等模拟植被的冠层反射率[24]。相比其他冠层辐射传输模型,该模型的参量更少,且容易获得[25];因此,本研究选择PROSAIL模型模拟玉米冠层反射率。

通过RedEdge-MX红边相机提供的波谱响应函数将模拟得到的高光谱反射率数据重采样到多光谱对应的宽波段,以便与实测的多光谱数据进行比较。依据试验日期、观测时间和经纬度计算当地不同观测时间下的太阳天顶角和方位角(表2)。固定PROSAIL模型的其他参数,基于太阳角度参数即可模拟同一天不同时间的光谱反射率,进而可以计算理想状况下的植被指数。模型中的输入参数(结构参数、含水量和干物质等)依据Li等[16]的研究设定,简述如下：叶片结构参数(N),取值为1.5;LCD,取值范围为10～80;叶片类胡萝卜素含量(Lcar,单位为μg·cm-2),取值范围为5～20;叶片花青素含量(Lant,单位为μg·cm-2),取值范围为0～10;叶片褐色素含量(Lbrown,单位为μg·cm-2),取值范围为0～1;等效水厚度(Ewt,单位为cm),取值范围为0～0.06;干物质含量(Dm,单位为g·cm-2),取值范围为0.002 8～0.005 4;LAI,取值范围为1～7;热点效应因子(hspot),取值范围为0～1;土壤因子(Psoil),取值范围为0～1;观测天顶角(VZA),取值为0。

表2 PROSAIL模型中太阳天顶角和方位角参数的具体设置

1.4 植被指数

根据前人的研究成果选取了13个与叶绿素相关的植被指数,用于分析不同观测时间光谱特征的变化及对玉米CCD估算的影响,具体包括NDVI[26]、绿色归一化植被指数(GNDVI)[27]、红边叶绿素指数(CIrededge)[28]、结构不敏感色素指数(SIPI)[29]、红边归一化植被指数(NDVIrededge)[30]、MERIS陆地叶绿素指数(MTCI)[31]、改良叶绿素吸收指数(MCARI)[32]、转化叶绿素吸收反射指数(TCARI)[33]、差异植被指数(DVI)[34]、增强植被指数(EVI)[35]、优化土壤调节植被指数(OSAVI)[36]、植被色素比率(PPR)[37],以及由TCARI和OSAVI之比构建的组合指数(TCARI/OSAVI)[38]。

1.5 精度评价

本研究以玉米每个生育时期每个时间段获取的植被指数对CCD进行建模和验证,样本总计49个,随机采取35个用于模型的构建,剩余样本用于模型的精度验证。使用决定系数(R2)、均方根误差(RMSE)和归一化均方根误差(NRMSE)对模型进行综合评价。

2 结果与分析

2.1 不同时间地物光谱反射率对比

对比抽雄吐丝期、籽粒建成期、乳熟期和蜡熟期10AM、11AM、1PM和2PM试验小区内所有像元的平均反射率可以发现,不同观测时间玉米的绿、红边和近红外的波段反射率差异较大(图2),其中近红外的差异最大,以乳熟期为例,10AM近红外反射率最低的为0.47,11AM近红外反射率最高的为0.84,二者之差达到0.37。不同时间,蓝、红波段的反射率差异相对较小。

10AM、11AM、1PM、2PM对应的时间分别为10：00—10：59、11：00—11：59、13：00—13：59、14：00—14：59。下同。图a～d分别对应于抽雄吐丝期(VT/R1)、籽粒建成期(R2)、乳熟期(R3)、蜡熟期(R5)。10AM, 11AM, 1PM, 2PM represent 10：00—10：59, 11：00—11：59, 13：00—13：59, 14：00—14：59, respectively. The same as below. Fig. a-d represent the tasseling silking stage (VT/R1), blister stage (R2), milk stage (R3), and dough stage (R5), respectively.图2 不同观测时间下的波段反射率变化Fig.2 Changes of waveband reflectance at different observation time

不同生育时期,各观测时间的影响差异有所不同,以籽粒建成期的差异最小。原因可能是,籽粒建成期的风速相较于其他3个时期最小,平均风速为2.37 m·s-1[数据来自欧洲中期天气预报中心(ECMWF)]。冠层光谱数据不仅受观测几何条件的影响,也会受到冠层结构的影响。风会引起冠层波动,风速较大的时候会导致冠层结构产生较大的波动;而风速较小的时候,冠层结构相对稳定,冠层光谱差异较小。PROSAIL模拟的结果也证实了这一结论。综上所述,玉米的光谱反射率在不同的观测时间存在差异,且不同波段的响应不同。

2.2 不同时期CCD相关特征筛选

虽然波段反射率在不同观测时间存在差异,但并没有显示出一定的规律性,且受背景、天气等因素影响较大;为此,本研究选用13个植被指数,分析其在不同生育时期与CCD的相关性(表3),每个生育时期各196个样本。在不同生育时期,同一植被指数与CCD的相关系数有所不同,但相关性一致。4个生育时期,与CCD相关性最高的前5个植被指数均是NDVIrededge、MTCI、CIrededge、GNDVI和NDVI,证明在玉米全生育时期红边波段均对CCD有较好的响应。后续,重点对这5个植被指数随观测时间的变化及其对CCD估算的影响进行分析。

表3 植被指数与冠层叶绿素密度(CCD)的相关性

2.3 植被指数的动态变化

2.3.1 模拟植被指数的动态变化

利用PROSAIL模型对不同观测时间选定的5种植被指数进行模拟,发现模型模拟的植被指数(以下简称为模拟植被指数)在同一生育时期的不同观测时间下存在微小的差异(图3)。一方面,这说明由不同观测时间所引起的玉米冠层光谱的微小变化可被PROSAIL模型模拟到;另一方面,这说明不同观测时间对PROSAIL模型模拟的植被指数的影响不明显。尽管不同时间模拟植被指数的差异较小,但其整体变化趋势是越接近中午,值越低。进一步比较不同生育时期同一观测时间模拟植被指数的变化,发现包含红边波段的NDVIrededge,CIrededge和MTCI等指数在不同生育时期间差异较大,不含红边波段的指数差异略小,说明红边波段在体现不同生育时期玉米长势差异方面更有优势,这点在前人的研究中也有所体现[39]。

NDVI,归一化植被指数;GNDVI,绿色归一化植被指数;CIrededge,红边叶绿素指数;NDVIrededge,红边归一化植被指数;MTCI,MERIS陆地叶绿素指数。下同。NDVI, Normalized difference vegetation index; GNDVI, Green normalized difference vegetation index; CIrededge, Red edge chlorophyll index; NDVIrededge, Red edge normalized difference vegetation index; MTCI, MERIS terrestrial chlorophyll index. The same as below.图3 基于PROSAIL模型模拟的植被指数在不同观测时间下的变化Fig.3 Changes of vegetation indexes simulated by PROSAIL model under different observation time

2.3.2 实测植被指数的动态变化

基于每个生育时期不同观测时间实际获取的无人机影像测算的植被指数(以下简称为实测植被指数)与基于PROSAIL模型模拟得到的植被指数差异明显,但整体上实测植被指数的值同样表现为越靠近12：00值越低的趋势(图4)。在抽雄吐丝期,实测植被指数值最低的是1PM时段;在籽粒建成期,NDVI、GNDVI、NDVIrededge和CIrededge值最低的是1PM时段,MTCI值最低的是2PM时段;在乳熟期,实测植被指数值最低的是11AM时段;在蜡熟期,NDVI、NDVIrededge,CIrededge和MTCI值最低的是1PM时段,而GNDVI值最低的是11AM时段。比较不同生育时期植被指数的变化,我们发现,NDVI、GNDVI、NDVIrededge、CIrededge和MTCI均在抽雄吐丝期达到最高值,籽粒建成期与乳熟期的数值相近,蜡熟期的值最低。这可能是因为,在抽雄吐丝期玉米的生长最为旺盛,此后随着玉米生长转入生殖生长阶段,与长势相关的植被指数值开始降低。

图4 基于无人机影像实测的植被指数在不同观测时间下的变化Fig.4 Changes of measured vegetation indexes based on unmanned aerial vehicle images under different observation time

2.4 不同时间植被指数与CCD的相关性变化

2.4.1 模拟植被指数与CCD的相关性变化

测算模拟植被指数与玉米CCD的相关性(图5),结果显示,整体上各植被指数与CCD的相关性均较高(均大于0.7),其中,包含红边波段的CIrededge,NDVIrededge和MTCI等指数与CCD的相关性相对更高,表明玉米反射光谱的红边位置与CCD更为敏感,这与Miller等[40]的研究一致。综合玉米所有生育时期的5个植被指数发现,模拟植被指数与CCD的相关性在同一生育时期不同观测时间下的差异不明显。主要原因是,理想情况下模拟植被指数随不同观测时间的变化差异较小。因此,仅通过模拟数据来研究不同观测时间玉米光谱反射特性与CCD相关性的差异存在一定的局限性。

2.4.2 实测植被指数与CCD的相关性变化

与模拟植被指数不同,不同观测时间下实测植被指数与CCD的相关系数存在较大的变化,且不同生育时期和不同指数间的差异不一致(图6)。MTCI、CIrededge和NDVIrededge与CCD的相关性随时间的变化规律基本一致,主要是由于这3种植被指数均是近红外和红边波段的比值或归一化形式,而波段反射率随不同观测时间的变化主要表现在近红外波段,红边波段变化较小。NDVI、GNDVI与CCD的相关性随时间变化的规律同上述3种指数并不一致,主要是由于这两种植被指数所用的光谱波段随不同观测时间的变化不一致。NDVI使用近红外和红波段的组合方式,GNDVI将NDVI中的红波段替换成绿波段。在不同观测时间,近红外波段反射率变化最明显,绿波段也有一定的变化,而红波段几乎不变;因此,这2种植被指数的变化是不一致的。在抽雄吐丝期,这两种指数的值在不同观测时间存在一定的波动。这可能是因为该时期玉米还没完全结束营养生长,不同施肥和密度处理的玉米长势存在差异,长势慢的小区冠层光谱受非绿色组分的影响较大,导致了植被指数均值与CCD相关性的波动。在抽雄吐丝期、籽粒建成期和乳熟期,包含红边波段的CIrededge,NDVIrededge和MTCI与CCD的相关性随时间波动不明显,表明在不同观测时间玉米CCD的光谱响应方面,包含红边波段的植被指数具有更好的稳定性。在蜡熟期,10AM、11AM和1PM时段的NDVI、GNDVI、CIrededge、NDVIrededge、MTCI与CCD的相关性差异相对较小,在2PM时段上述相关性均降低。这可能是由于2PM时段的太阳天顶角达到44.35°,高天顶角导致大气散射效应增强,太阳光照强度变弱,限制了太阳光线在玉米冠层内部的反射和吸收,导致传感器获取的冠层光谱携带的叶绿素信息变少[18],进而导致此时测得的植被指数与CCD的相关性降低。不同时间相关性变化最小的是籽粒建成期和乳熟期。这可能是因为,这2个时期玉米均进入到生殖生长阶段,叶绿素含量不再增加,不同处理的玉米长势差异缩小,非绿色组分及其他因素的反射对冠层光谱造成的影响变小。总的来看,选定的5种植被指数的实测值与CCD的相关性在所有生育时期均较高,表明这些植被指数在玉米不同生育时期的变异性较小,适用于玉米全生育期的叶绿素含量监测研究[41]。

图6 基于无人机影像实测的植被指数与冠层叶绿素密度(CCD)的绝对相关系数Fig.6 Absolute correlation coefficient within the measured vegetation indexes based on unmanned aerial vehicle images and canopy chlorophyll density (CCD)

2.5 不同观测时间对CCD估算的影响

基于设定的PROSAIL模型的参数范围,随机生成10 000组模型输入参数,并利用生成的LCD和LAI计算CCD,使用PROSAIL模型模拟这些参数对应的5种植被指数,将获得的模拟的植被指数和CCD数据组合用于训练RF模型。基于每个生育时期的实测数据和不同观测时间的太阳天顶角和方位角,利用PROSAIL模型模拟对应的植被指数,对RF模型进行验证(图7)。结果显示,10 000组模拟数据的建模精度达到了极显著水平(P<0.01),R2为0.99。综合所有生育时期和观测时间,RMSE和NRMSE均以11AM和1PM时段最低,且二者相差不大,这表明,理想情况下,基于同一天内不同时间的模拟数据估算CCD的精度存在差异,验证误差随不同时间变化的规律整体表现为越接近12：00越小。

RMSE,均方根误差;NRMSE,归一化均方根误差。下同。RMSE, Root mean square error; NRMSE, Normalized root mean squared error. The same as below.图7 基于PROSAIL模型模拟的植被指数对冠层叶绿素密度(CCD)估算的影响Fig.7 Influence of the vegetation indexes simulated by PROSAIL model on canopy chlorophyll density (CCD) estimation

基于不同时间的观测数据,利用RF方法构建玉米CCD的估算模型(图8)。综合所有生育时期和观测时间,模型精度均达到了极显著水平(P<0.01),R2最低的为0.53,最高的为0.80。用实测数据构建的CCD估算模型的精度在不同观测时间下有所差异,但整体规律较为一致,表现为越接近12：00的时段精度越高,这与PROSAIL模型模拟的结果较为一致。在拟合精度上,抽雄吐丝期,乳熟期和籽粒建成期的趋势较为一致,以1PM时段最好。蜡熟期的整体精度较差,这可能与玉米后期的衰老有关。玉米的衰老一般是从下往上开始的,导致整体的CCD降低,而冠层光谱由于过高的覆盖度仅探测到上层未衰老的信息。抽雄吐丝期、乳熟期和蜡熟期模型的误差趋势在不同时间基本一致,整体表现为11AM时段的误差最低。籽粒建成期的误差趋势表现为10AM时段的误差最低。综上所述,不同观测时间对CCD估算模型的精度存在一定的影响,整体表现为越接近12：00,精度越高。

R2,决定系数。R2,Coefficient of determination.图8 基于无人机影像实测的植被指数对冠层叶绿素密度(CCD)估算的影响Fig.8 Influence of the measured vegetation indexes based on unmanned aerial vehicle images on canopy chlorophyll density (CCD) estimation

3 结论

本研究通过设计无人机多频次的田间观测试验,结合PROSAIL模型重点分析了不同观测时段下植被指数的变化规律及其对CCD估算精度的影响。主要结论如下：(1)不同观测时间的波段反射率和植被指数存在差异,其中,近红外波段反射率的变化最明显,与CCD相关性强的植被指数在越接近12：00的时段值越低。(2)在同一生育时期的不同观测时间,模拟植被指数与CCD之间的相关性差异不明显,但实测植被指数与CCD的相关性存在较大差异,且在不同生育时期的差异不一致。(3)UAV数据获取时间对CCD的估算结果有影响,数据获取时间越接近12：00,估算模型的精度越高,但不同生育时期的最佳获取时间有一定差别,抽雄吐丝期、乳熟期和蜡熟期宜利用接近12：00的时段获取的数据估算CCD,籽粒建成期在10AM时段获取的光谱估算CCD最优。本文结果初步揭示了无人机多光谱数据观测时间对CCD估算的影响,后续将深入探讨如何消除由于观测时间导致的估算误差,为构建一个对太阳天顶角和方位角不敏感的CCD精准估算模型提供基础,从而为作物长势的精准监测提供技术支撑。