项方林,李鑫格,马吉锋,刘小军,田永超,朱艳,曹卫星,曹强
(南京农业大学/国家信息农业工程技术中心/农业农村部农作物系统分析与决策重点实验室/智慧农业教育部工程研究中心/江苏省信息农业重点实验室,南京 210095)
【研究意义】实时、准确地预测作物单产是进行作物精确管理的关键环节。时序植被指数能够反映作物全生育时期的长势信息,从而为田间快速、准确地预测作物单产提供技术支撑。【前人研究进展】基于卫星遥感影像,研究者利用作物关键生育时期获取的冠层植被指数构建了不同的估产方法[1]。然而,基于单时相遥感影像信息所构建的单产预测模型易受环境影响,普适性欠佳。相较于单一生育时期的冠层光谱信息,时序植被指数蕴含了作物全生育时期的生长发育动态信息,具有良好的单产预测的潜力[2-5]。目前,研究者多基于卫星遥感获取高时间分辨率的植被指数信息,如基于AVHRR和MODIS卫星获取时序归一化植被指数(normalized difference vegetation index,NDVI)等,其在植被分类[6]、物候期提取[7]、作物种植面积估测[8]等方面得到了广泛应用。前人在估产方面也进行了研究,如陈鹏飞等[9]通过提取基于环境减灾卫星拟合的时序 NDVI曲线特征参数进行小麦估产,结果表明利用多元逐步线性回归建立的估产模型效果较好。然而,由于卫星影像处理的复杂性、时效性等问题限制了该技术为当季生产管理决策提供支持的潜力[10]。随着近地面传感器的发展,其低成本、易操作和快解译等优点成为作物生长监测不可或缺的技术途径[11]。已有研究者利用近地面传感器获取时序植被指数并基于局部加权回归、双 Logistic模型等进行数据重构,从而估测作物生育时期及相关农学参数[7,12]。上述研究多基于植被指数NDVI,然而,有研究表明NDVI在植被覆盖度较高时具有明显的饱和效应[13]。近年研究发现,基于红边波段拟合的植被指数与作物氮素营养状况关系更为密切[14]。如THOMPSON等[15]研究揭示了基于红边波段的归一化红边植被指数(normalized difference red edge,NDRE)能有效解决饱和问题。相关学者还利用红边和近红外波段组合的比值植被指数以及包含有红边波段的3波段指数有效地估测了小麦等氮素营养状况[16-17]。RapidSCAN CS-45(Holland Scientific Inc.,Lincoln,USA)是一种包含红边波段在内的主动冠层传感器,已有学者利用该传感器获取的植被指数NDRE和NDVI在诊断作物氮素营养状况及估测产量时取得较好的效果[18-20]。【本研究切入点】基于时序植被指数进行作物单产估测的研究多是利用卫星遥感数据,但其时间分辨率和空间分辨率相对较低,且易受外界条件影响,而应用主动冠层传感器能有效避免这些不足。【拟解决的关键问题】本研究拟应用主动冠层传感器 RapidSCAN CS-45获取冬小麦冠层时序植被指数,并基于时序曲线特征参数构建冬小麦单产预测模型,揭示不同品种及氮素水平下植被指数的时序变化特征,构建并评价基于时序植被指数的冬小麦单产预测模型。
试验1:于2017—2019年在江苏省兴化市万亩高产高效粮食产业园(119°53′ E, 33°48′ N)进行。试验区土壤有机质 26.5 g·kg-1,全氮 2.07 g·kg-1,有效磷52.54 mg·kg-1,速效钾 93.66 mg·kg-1。供试品种为镇麦12号、扬麦23号、宁麦13号,均为江苏省主栽品种。施氮水平分别为纯氮0(N0)、90 kg·hm-2(N1)、180 kg·hm-2(N2)、270 kg·hm-2(N3)、360 kg·hm-2(N4)。小麦行距为0.25 m,小区面积为63 m2(7 m×9 m),将每个小区一分为二,一侧为测试区,另一侧为测产区(图 1)。试验采用完全随机设计,重复3次,分别于2017年11月8至9日和2018年11月1日于南北行向人工播种,基本苗为2.25×106株/hm2。结合整地一次性施用105 kg·hm-2P2O5,钾肥共施用120 kg·hm-2K2O,氮肥、钾肥基追比5∶5,基肥在播种前施入,追肥在拔节期施入。在作物生长阶段,及时除草防治病虫害,其他栽培措施同一般高产田。
试验2:于2018—2019年在相同园区进行,供试品种为扬麦23号。试验设置2个播期(11月4日和24日),3个密度(1.8、2.7、3.6×106株/hm2)以及4个氮素水平(0、180、240、300 kg·hm-2)。其中0氮水平不设重复,其他处理3次重复,小区面积为30 m2(5 m×6 m),其他试验设计、栽培管理措施与试验1一致。
图1 试验1田块分布图Fig. 1 Map of experiment 1
1.2.1 主动冠层传感器获取光谱数据 本研究使用的RapidSCAN CS-45包括红光(R,670 nm)、红边(RE,730 nm)、近红外(NIR,780 nm)3个波段。可以自动收集每个波段的光谱反射率和地理坐标数据,并以2.5 Hz(每0.4 s读取一次)的频率将其记录在传感器的存储模块中,数据可以通过 PC软件导出为.csv文件。仪器视场角为横向45°,纵向10°,仪器自带光源,且不受天气条件的影响,较方便轻巧的特点有利于测试。在冬小麦越冬—返青期,每隔15 d测试1次;返青期之后,每隔3—5 d测试1次,2个冬小麦生长季共测试63次(2017—2018冬小麦生长季测试30次,2018—2019冬小麦生长季测试33次)。测定时,在每个待测小区选择中间两行手持传感器匀速前行,并使传感器距离冬小麦冠层的高度约100 cm,且与地面保持平行,每行约20—30个测试点,每个小区测定的平均值代表该小区的测试值,本研究基于NIR、R和RE波段计算了各小区的NDRE和NDVI(公式(1—2))。
式中,NIR表示近红外波段反射率,RE表示红边波段反射率,R表示红光波段反射率。
1.2.2 冬小麦单产的测定 冬小麦收获时,在各试验小区选取3个1 m2样方,单独收割、晾晒、脱粒并称重,取平均值作为该小区的产量。
累积生长度日数(accumulated growing degree days,AGDD)是描述时序植被指数动态变化的时间驱动因子,由从播种日期到各测试日期生长度日数值累加而成[21-22],如公式(3)所示。本研究所采用的温度数据来源于试验点所在的气象站。
式中,i为冬小麦播种后天数,Tmax和 Tmin分别为当日最高温度和最低温度,Tbase为冬小麦生长的基点温度,本研究中设为0℃[23-24]。
相对AGDD(relative AGDD,RAGDD)由AGDD归一化计算得到,如公式(4);相对NDRE(Relative NDRE,RNDRE)和相对 NDVI(Relative NDVI,RNDVI)分别由NDRE和NDVI归一化计算得到,如公式(5)和公式(6)。
式中,AGDDi为测试当日的 AGDD值,AGDDharvest为收获时的AGDD值,即全生育时期的累积生长度日数。NDREi和 NDVIi为测试当日该处理的 NDRE和NDVI值,NDREmax和NDVImax为该品种全生育时期NDRE和NDVI最大值。
本研究中,采用双Logistic函数模型拟合时序植被指数的动态变化,建立了 RAGDD和 RNDRE、RNDVI的定量关系。1stOpt15PRO(7D-Soft High Technology Inc.,Beijing,China)主要用于时间序列数据的分析,利用双Logistic曲线拟合模型,并计算出模型中的拟合参数。模型拟合参数由拟合方程在1stOpt15PRO软件下输出,参照Fischer[25]模型,如公式(7)计算:
式中,第1个Logistic函数代表冬小麦的生长过程,第2个Logistic函数代表冬小麦的衰老过程。y为拟合后的植被指数值,A1为归一化后的植被指数最大值,A2为A1与成熟期植被指数的差值;a和c分别表示冬小麦生长和衰老过程中2个Logistic曲线模型的拐点处斜率,即植被指数增长或下降的最大速度,与这 2个拐点对应的RAGDD分别为b、d。
基于拟合后的时序RNDRE和RNDVI曲线,本研究提取的冬小麦生长季特征参数如下:
(1)RNDRE和RNDVI最大值
冬小麦RNDRE和RNDVI最大值一般出现在其孕穗—抽穗期,其值大小反映了最终形成营养器官的状态,体现了后期向籽粒转运物质的能力。
(2)累积RNDRE和RNDVI
冬小麦分蘖越冬期到成熟期的RNDRE和RNDVI曲线积分即累积RNDRE和RNDVI,反映了冬小麦各个生育时期的平均生长状况。
(3)RNDRE和RNDVI增长速率
RNDRE和RNDVI增长速率是指从分蘖越冬期到最大值的相对变化速率,通过公式(8)和(9)计算,反映了冬小麦在这一过程中的生长发育状态。
式中,RNDREincrease和 RNDVIincrease表示生长过程RNDRE和 RNDVI的增长速率,RNDREmax和RNDVImax表示冬小麦生长季RNDRE和RNDVI最大值,RNDREinitial和RNDVIinitial表示冬小麦分蘖越冬期的RNDRE和RNDVI值;RAGDDm表示冬小麦达到RNDRE和RNDVI最大值的RAGDD,RAGDDi表示分蘖越冬期出现的RAGDD。
(4)RNDRE和RNDVI下降速率
RNDRE和RNDVI下降速率是指从抽穗—开花期到成熟期的相对变化速率,通过公式(10)和(11)计算,反映了冬小麦衰老过程的生长发育状态。
式中,RNDREdecrease和 RNDVIdecrease表示衰老过程的RNDRE和RNDVI变化速率,RNDREmature和RNDVImature表示冬小麦成熟期的RNDRE和RNDVI值,RNDREmax和RNDVImax表示冬小麦生长季RNDRE和RNDVI最大值;RAGDDmature表示冬小麦成熟期的RAGDD值,RAGDDm表示冬小麦达到RNDRE和RNDVI最大值的RAGDD。
本研究利用了试验1的数据集进行模型的构建,利用试验2的数据集进行模型的验证。首先,构建基于单生育时期 NDRE和 NDVI的冬小麦单产预测模型,并通过多元逐步线性回归法构建基于多生育时期NDRE和NDVI的多变量冬小麦单产预测模型。其次,基于时序RNDRE和RNDVI曲线提取的各特征参数建立单产预测模型,利用独立试验数据对所构建的各单产预测模型进行验证比较。
通过计算每个模型的决定系数(R2,coefficient of determination)来评价模型的表现,但R2不适用于评估不同形式、不同参数个数的模型之间的预测和拟合性能,因为R2较高但参数较多的模型存在过度拟合的风险。因此,本研究进一步使用赤池信息准则(akaike information criterion, AIC)对模型进行比较分析,AIC的基本观点是所选择的模型是为了准确预测未来的数据,而不是为了推断校准数据的“真实分布”[26-27]。由于拟合模型的真实预测性能在很大程度上取决于模型的自由参数的数量,因此使用AIC选择模型剔除大量参数并阻止过度拟合[28]。AIC在R语言3.5中计算,其计算方法如公式(12)所示:
式中,k是模型中参数的数量,Lˆ是模型中极大似然函数值似然函数。在进行预测同一变量时,通常选取AIC值最小的模型。
通过比较预测值与实测值的R2、相对均方根误差(RRMSE,%)和相对误差(RE,%),来评价模型的验证效果。R2越高,RRMSE和RE越低,表示该模型预测单产的精度和准确性越高。RRMSE和RE的计算公式如公式(13)和(14)所示:
式中,Pi、Oi分别表示单产的预测值和实测值,表示单产实测值的平均值,n表示模型测试样本数。
本研究采用 Microsoft Excel 2016(Microsoft Corporation,Redmond,Washington,USA)进行数据整理,结合SPSS 25.0(SPSS Inc., Chicago, IL, USA)软件进行统计和分析;利用 Origin 2016(OriginLab Corporation,Northampton,MA,USA)作图。
从表1可以看出,植被指数NDRE和NDVI与小麦单产的相关关系在不同品种及生育时期表现并不相同,总体而言,在拔节期相关性较差,在孕穗期—灌浆期相关性较好。综合所有数据,植被指数NDRE与冬小麦单产在孕穗期—灌浆期都表现出稳定的相关关系(R2>0.81);而 NDVI与单产在孕穗期—灌浆期的关系略差于 NDRE,R2为0.76—0.80。可见,植被指数 NDRE和 NDVI具有良好的估测单产的潜力,NDRE略优于NDVI,且在孕穗期—抽穗期最佳。
表1 不同生育时期植被指数与冬小麦单产的相关关系Table 1 Coefficients of determination (R2) for the relationship between vegetation indices (NDRE and NDVI) and yield at different growth stages across varieties and years
为提高预测单产的准确性和稳定性,本研究尝试通过多元逐步线性回归方法,利用多生育时期NDRE和NDVI构建多变量冬小麦单产预测模型(表2)。结果表明,对于同一植被指数而言,无论利用几个时期的数据,其R2、RRMSE和RE没有明显的差异。因此,表明利用较少的时期就可以满足单产预测结果的准确性。对于NDRE而言,当利用拔节期和孕穗期的数据时,其AIC最小,模型表现最优。对于NDVI而言,当利用孕穗期、抽穗期和开花期3个生育时期的数据时,AIC最小,模型表现较优。
冬小麦全生育时期冠层植被指数 NDRE(图 2)和NDVI(图3)在不同品种和氮水平间变化趋势基本相同,随着AGDD的增加,植被指数NDRE和NDVI呈现先增加后降低的趋势。在冬小麦生长前期,由于冬小麦未封行,受土壤背景的影响不同氮水平之间植被指数差异较小,而当AGDD大于650,冬小麦封行之后(即拔节期),植被指数值随着施氮量的增加而升高,处理间差异较为明显。
图2 不同氮素水平下基于AGDD的NDRE变化动态Fig. 2 Dynamic changes of NDRE based on AGDD with different nitrogen rates
结合作物生长发育规律,基于双Logistic模型拟合冬小麦生长发育过程中植被指数RNDRE和RNDVI的动态变化,图4表明,随着RAGDD的增加,时序RNDRE和RNDVI呈现出相似的趋势,即先增加后下降的动态变化规律。由表3可知,冬小麦生长过程中最大速度a和衰老过程中的最大速度c因施氮水平而异,2种植被指数均以N3、N4水平最高,在N0水平下最低。且2种植被指数的生长拐点b和衰老拐点d均在N0水平下先出现,表明在低N水平下的冬小麦生育进程缩短,生长发育差,容易造成冬小麦单产降低。从模型拟合精度来看,施氮水平越高,模型的精度也越高,N3和N4水平下的RNDRE和RNDVI曲线拟合R2达到0.95以上,RMSE小于0.045,且N3、N4的曲线基本重合(图4),可能是因为氮素过量,造成冬小麦群体过大,植被指数趋于饱和。相较于其他氮处理,N0处理的曲线拟合精度略低,R2<0.80,RMSE<0.066。总体而言,时序 RNDRE曲线拟合精度较 RNDVI曲线拟合精度高,可能是因为红边波段的植被指数更为稳定。
图3 不同氮素水平下基于AGDD的NDVI变化动态Fig. 3 Dynamic changes of NDVI based on AGDD with different nitrogen rates
表2 基于多元逐步线性回归基的不同生育时期植被指数单产预测模型Table 2 Stepwise multiple linear regression models based on vegetation indices for estimating winter wheat yield at different growth stages
图4 不同氮素水平下RNDRE(A)和RNDVI(B)时序曲线变化动态Fig. 4 The dynamic changes of time-series curve of RNDRE (A) and RNDVI (B) with different nitrogen rates
表3 不同氮处理下基于RAGDD的时序RNDRE和RNDVI曲线拟合参数Table 3 The fitting parameters of time-series curve of RNDRE and RNDVI based on RAGDD with different nitrogen rates
由表4可知,基于双Logistic函数拟合提取的曲线特征参数值因不同施氮水平而异,随着施氮水平的增加,RNDRE和RNDVI最大值、累积值、增长速率和下降速率值也是上升的。RNDRE和RNDVI曲线拟合提取的4个特征参数均在N0水平下达到最小值,在高氮水平(N3、N4)下可以达到最大值。从表 5可以看出,时序曲线特征参数与冬小麦单产具有稳定的相关性,其中最大值、累积值和增长速率与单产相关性较好。对于RNDRE和RNDVI而言,AIC最小的模型均是基于增长速率建立的单产模型。
表4 不同氮素水平下RNDRE和RNDVI时序曲线特征参数Table 4 The characteristic parameters of time-series curve of RNDRE and RNDVI with different nitrogen levels
为了探明单产预测模型的准确性和效果,用独立试验数据对上述模型进行验证。基于单生育时期NDRE和NDVI构建的单产预测模型,如表6所示,其在拔节期的验证结果较差,R2小于0.6;其在孕穗期—灌浆期的验证结果较好,R2大于 0.68,RRMSE和RE小于20%。基于多元逐步线性回归构建的基于NDRE多生育时期单产预测模型验证结果略优于NDVI,但与单生育时期模型相比无显著提高。从表7可以看出,基于时序RNDRE曲线特征参数构建单产模型,其最大值和累积值所构建的单产模型验证结果较好,R2大于0.8,RRMSE和RE小于10%;而增长速率和下降速率构建的单产模型验证结果较差,R2小于0.5。基于时序RNDVI曲线特征参数构建的单产模型也展示了相似的结果,但略差于时序RNDRE曲线特征参数所构建的模型(图 6)。综合考虑,基于RNDRE最大值和累积值构建的单产模型验证结果较好,其验证效果优于利用单生育时期或多生育时期的预测模型,且优于基于NDVI构建的单产模型。
表5 基于RNDRE和RNDVI时序曲线特征参数的冬小麦单产预测模型Table 5 The characteristic parameters of time-series RNDRE and RNDVI curve based prediction models of winter wheat yield
表6 基于NDRE和NDVI的冬小麦单产预测模型验证结果Table 6 Validation results of yield prediction model based on NDRE and NDVI in winter wheat
表7 基于RNDRE和RNDVI时序曲线特征参数的冬小麦单产模型验证结果Table 7 Validation results of evaluation with winter wheat yield prediction models based on the characteristic parameters of time-series RNDRE and RNDVI curve
图5 基于 RNDRE(A:最大值;B:累积值)和 RNDVI(C:最大值;D:累积值)时序曲线特征的单产预测模型验证结果Fig. 5 Validation results of evaluation with winter wheat yield prediction models based on the characteristic parameters of time-series RNDRE (A: Maximum; B: Accumulative value) and RNDVI (C: Maximum; D: Accumulative value) curve
植被指数被广泛应用在基于卫星遥感的单产预测研究中,研究表明获取于MODIS和HJ-1A/1B影像的植被指数NDVI等能有效预测油菜、小麦等作物的产量[1,29-30]。近地面遥感也有利用植被指数进行单产预测的研究,HASSAN等[31]利用无人机搭载多光谱传感器获取 NDVI,结果表明在小麦灌浆期预测单产效果最好。前人研究大多基于传统植被指数NDVI等进行单产预测,而本研究使用了基于红边波段的植被指数NDRE,有关研究表明基于红边波段的植被指数能有效改善其与单产的相关性[32]。与NDVI和增强型植被指数(EVI)相比,前人研究表明宽动态范围植被指数(WDRVI)能准确预测小麦产量[33]。本研究比较了不同生育时期植被指数NDRE和NDVI与冬小麦单产的相关性(表1),除拔节期关系较差外,其他生育时期均较好,且总体而言,基于NDRE的结果优于NDVI,原因可能是NDRE带有红边波段,其植被指数稳定性较高,也有研究表明带有红边波段的植被指数预测产量时效果更好[17]。通过多元逐步回归分析结果表明,利用多个生育时期的植被指数所构建的模型可进一步提高单产预测的精度。上述研究结果表明,多生育时期的植被指数可以更为全面反映产量的形成过程,利用其预测小麦产量效果更佳。
利用单时相数据预测作物物候期和单产易导致一定误差,而多时相的时序数据能够较准确预测作物物候期,具有提升估产精度的潜力。前人利用环境减灾卫星获取时序植被指数并通过重构植被指数时序数据获取曲线特征参数,利用其构建的估产模型效果较好[9]。已有研究基于获取的时序EVI并利用S-G滤波重构时序EVI曲线,其曲线特征参数中最大值与产量相关性较好,也有研究基于NDRE构建了水稻早中期动态模型,其曲线特征平台值与产量关系较好[34-35]。本研究通过近地面主动冠层传感器获取的植被指数,构建时序曲线并提取特征参数预测冬小麦单产。其中,时序曲线最大值可以反映冬小麦群体生长的潜力,基于其构建的单产预测模型准确性较高;时序曲线增长速率可以反映冬小麦生长季的营养生长期状态;时序曲线下降速率反映了冬小麦生殖生长期的状态;累积植被指数反映出冬小麦全生育时期的生长状态。本研究结果表明,RNDRE时序曲线最大值、累积值与冬小麦单产具有良好的线性关系,其结果略优于基于RNDVI时序曲线最大值和累积值构建的单产模型,原因可能是带有红边波段的植被指数RNDRE稳定性较高,饱和程度较低,且产量的形成过程是一个不断累积的过程,因此,利用全生育时期RNDRE曲线特征参数具有良好估测小麦单产的潜力(图 6)。前人关于时序植被指数的研究多是基于卫星遥感获取的,并且多基于时序 NDVI,而本文是基于手持式传感器每隔3—5 d天测试1次,意在获取更丰富的时序植被指数数据,可以不受环境条件的影响。
作物生长过程具有时序动态性,充分利用不同时期间差异以及多时期的综合信息,可以提高对作物分类、生育期估测和单产预测等评价精度。已有研究基于 MODIS-NDVI时序数据进行作物分类、预测冬小麦生育期均表现较好[36-37]。JOHNSON[38]基于MODIS时序遥感产品对玉米等单产预测R2达到0.7以上。时序植被指数被同化到作物模型中,也能够提高作物估产的精度[39-40]。作物不同时期和氮水平的长势显著影响植被指数,进而会影响时序植被指数模型的参数。本研究采用双Logistic函数拟合时序植被指数曲线,能较好地反映冬小麦冠层植被指数的动态轨迹。前人研究多基于单生育时期或几个生育时期的光谱数据预测小麦单产,其结果可能较好,但并不能反映产量的形成过程。产量形成是一个不断累积的过程,全生育时期的时序光谱可以更加全面反映这一形成过程,因此时序曲线具有更好的预测小麦产量的潜力。作物时序植被指数的动态模型采用播种后天数、每年中的天数(days of year)等作为时间驱动因子的研究较多[41-42],由于作物不同条件下的温光反应、基因型效应与阶段发育之间的动态关系使得作物全生育时期天数差异较大,导致模型缺乏动态性和广适性。本研究采用RAGDD为时间驱动因子可以较好地消除不同年份、品种对冬小麦时序植被指数模型的影响,而且模型参数较少,生物学意义比较明确。在冬小麦实际生产中获取时序光谱,通过构建本地化时序植被指数模型对当地冬小麦单产进行预测。由于影响冬小麦单产的因素很多,诸如品种、栽培条件、土壤背景、气候因素等,在未来研究中,需融合上述因素建立单产预测模型,从而提高单产预测的准确性和精度。本试验研究只在兴化市一个试验点进行,还需增加其他生态点试验来进行模型的验证与调优,从而提高模型的稳定性。
冬小麦生长季植被指数NDRE和NDVI具有良好的估测单产的潜力,NDRE略优于NDVI,且在孕穗期—抽穗期最佳。基于多元逐步线性回归构建的基于NDRE多生育时期单产预测模型验证结果略优于NDVI,与单生育时期模型相比无显著提高。基于RNDRE时序曲线最大值和累积值所构建的单产模型效果较好,其结果略优于基于RNDVI时序曲线最大值和累积值构建的单产模型。综合考虑,表明RNDRE最大值和累积值具有良好估测单产的潜力。研究结果为田间进行实时、准确的单产预测提供了技术参考,但是冬小麦单产受土壤、生态点、气候及生产水平影响较大,本研究的时序植被指数单产预测模型还需在不同生态区域检验其适用性和可靠性。