黄 茜 杨伟才 魏夏永 毛晓敏
(1.中国农业大学水利与土木工程学院, 北京 100083; 2.农业农村部作物高效用水武威科学观测实验站, 武威 733000)
玉米是黑龙江省种植的主要粮食作物之一,2018年玉米种植面积6.318×106hm2,占全省粮食作物种植面积的44.45%[1]。覆膜种植可以增加土壤温度,从而促进作物生长、提高作物产量[2-4]。叶面积指数(Leaf area index,LAI)是反映作物长势、胁迫程度的重要指标之一[5-7],因此精准监测不同覆膜处理下春玉米LAI,对准确掌握玉米生长状况具有重要意义。传统测量LAI的方法不仅费时、费力,还对作物有一定的破坏性。高光谱遥感技术能从光谱信息中获得作物生长的相关参数,可实现对作物生长状况的实时、无损、快速监测[8-9]。
基于高光谱遥感数据对作物LAI进行估算已有大量的研究成果,形成了较为成熟的LAI反演方法,主要分为统计经验模型反演、机器学习算法模型反演和物理模型反演[10-11]。统计经验反演通过建立作物LAI与光谱变量之间的统计关系求得LAI。袁媛[12]研究确定了夏玉米各个生育期LAI高光谱估算模型的最佳植被指数。LIANG等[13]根据光谱数据收集计算了43种不同的植被指数,经相关性分析,确定了用于估计LAI的最佳植被指数。谭昌伟等[14]研究发现,R810/R560为夏玉米LAI最佳估算植被指数。NGUYEN等[15]分别测定了两个不同品种水稻的冠层反射率,确定了水稻LAI的特征植被指数。ALI等[16]分析比较了不同方法获取的红边位置,并以此为因变量建立了柑桔LAI的最小二乘法估算模型。何小安等[17]基于相关性分析优选47个植被指数,构建了冬小麦生育前期LAI的反演模型。近年来,随着计算机技术的发展,机器学习方法逐渐用于作物生长指标的反演。VERRELST等[18]基于高斯过程回归(Gaussian processes regression,GPR)筛选敏感波段,并用于对LAI的估算。ROCHA等[19]采用支持向量机(Support vector regression,SVR)算法对LAI进行预测。YUAN等[20]采用随机森林(Random forest,RF)、人工神经网络(Artificial neural network,ANN)和SVR对大豆LAI进行回归分析。基于物理学光学模型的方法,例如PROSAIL模型在作物LAI反演研究中得到广泛的应用[21-23]。
已有研究表明,根据高光谱数据建立相应作物LAI的估算模型往往具有较好的效果,但植被指数和光谱变量都具有地域性,在不同地区会表现出较大的差异,建立适用于东北高寒黑土区春玉米LAI估算模型还需进一步探究[24-26]。目前,对玉米冠层高光谱特征研究已有许多成果,但对不同覆膜处理下春玉米冠层光谱曲线形状差异及覆膜对LAI的估算影响的研究鲜见报道[27]。本文利用高光谱数据估算不同覆膜处理下春玉米LAI,研究覆膜因子对LAI估算的影响,探讨利用敏感光谱变量估算混合各覆膜处理LAI的可行性,以期为覆膜春玉米种植长势监测、优化农业生产管理策略提供技术指导。
试验于2019年5—9月在黑龙江省齐齐哈尔市克山县北联镇北部战区空军保障部克山农副基地农场(48°15′N, 125°36′36″E)进行。该地区年平均气温2.4℃,有效积温2 400℃,年降水量为500 mm,无霜期122 d,年平均风速4 m/s。土壤质地为粉壤土,1 m土层平均田间持水量为0.41 cm3/cm3,饱和含水率为0.54 cm3/cm3,干容重为1.23 g/cm3。
春玉米种植品种为富单2号,种植日期2019年4月16日,收获日期为2019年9月16日。种植模式为1.1 m大垄双行南北向种植,播种方式为平播,耕作方式为免耕,播种深度约5 cm,垄行距40 cm,株距20 cm,单株玉米占地面积为0.11 m2;除草、喷洒农药等农艺措施以当地经验为准。底肥由二铵225 kg/hm2、硫酸钾120 kg/hm2和尿素42 kg/hm2混合拌匀施入,其中二铵含氮量不小于18%,含有效磷(P2O5)不小于46%;硫酸钾含有效钾(K2O)不低于50%;尿素含氮量不小于46.4%。相当于施入纯氮量60 kg/hm2,有效磷(P2O5)104 kg/hm2,有效钾(K2O)60 kg/hm2。本试验考虑覆膜、追肥2种因子,设置3种覆膜方式:不覆膜(M0)、覆可降解膜(西安新秦塑料有限责任公司,膜宽2.5 m,膜厚0.012 mm)(M1)、覆普通透明膜(哈尔滨市山水东华硕农膜厂,膜宽1.2 m,膜厚0.008 mm)(M2),均在播种次日人工覆膜。设置4种追氮水平:不追肥;在大喇叭口期追氮30 kg/hm2、吐丝期追氮30 kg/hm2;在大喇叭口期追氮70 kg/hm2、吐丝期追氮70 kg/hm2;在大喇叭口期追氮100 kg/hm2、吐丝期追氮100 kg/hm2。共12种处理,每种处理设置3个重复,共36个小区,小区随机布置。每个小区均为矩形,长度7.5 m,宽度8.8 m。
(1)生育期:将试验区中超过70%的春玉米植株具有一定的某生育期特征,作为进入该生长期的标准,记录每个小区进入各个生长期的时间。
(2)冠层反射率:使用美国Analytical Spectral Devices (ASD)公司生产的Field-Spec HandHeld2型手持式地物光谱仪测量春玉米冠层光谱反射率,光谱仪视场角25°,波段为325~1 075 nm,光谱采样间距为1 nm,光谱分辨率3 nm。采集时间为10:30—14:00,在拔节期、抽雄期、吐丝期、灌浆期、成熟期选择晴天无风无云天气测定玉米的冠层反射率,测定时传感器探头始终保持垂直向下,距玉米冠层1 m[28-29]。每个试验小区选择具有代表性的3个观测点,每个观测点连续采集6个冠层数据,剔除异常数据后,取平均值作为该小区的冠层光谱最终测量值。每个样品测量完毕后,及时对准标准白板进行校正,从而准确测定下一个样品。
(3)LAI:在收集春玉米冠层光谱数据后,在每个采样点取2株玉米,每个试验小区共6株,带回室内测量所有完全展开叶的长和宽,用叶片的实际面积乘以修正系数0.75计算LAI。
植被指数是由高光谱数据线性和非线性组合而成的各个数值,但一般植被指数只利用了特定的波段,没有携带全波段光谱信息,本文选择红光和近红外波段范围反射率的平均值计算植被指数,具体见表1。
表1 植被指数及其计算公式
利用Matlab编程计算原始光谱和一阶微分光谱任意两波段比值指数(RI=Ri/Rj)、归一化指数(NI=(Ri-Rj)/(Ri+Rj))及再比值指数(AI=(Ri/Rj)/[(Ri-Rj)/(Ri+Rj)]),计算3个指数与LAI的决定系数(Coefficient of determination,R2),并绘制3个指数与LAI的R2等势图。
为了探究玉米冠层反射率的特征,以及量化光谱数据与LAI之间的关系,采用Matlab编写计算光谱变量的程序,运用软件IBM SPSS Statistics 20、Excel 2010和OriginPro 9.0分析处理试验数据。
1.5.1光谱数据预处理
使用高光谱仪的光谱处理软件ViewSpecPro将田间测得原始光谱导入,每个小区测得光谱数据进行平均处理后分别将一阶导数[12]和原始光谱转成ASCII文本文档,之后导入Excel 2010,对原始光谱数据进行平滑处理(五点移动平均法)[37]。
1.5.2春玉米LAI估算模型
采用统计回归方法进行建模。
单变量线性模型
y=ax+b
单变量指数模型
y=aebx
单变量幂函数
y=axb
一元二次多项式
y=a+bx+cx2
式中y——预测值x——自变量
a、b、c——系数
1.5.3模型评价指标
采用R2评价春玉米生长指标与模型输出值之间的拟合效果,采用实测值与预测值之间的均方根误差(Root mean square error, RMSE)对模型的预测精度进行评价;采用剩余预测偏差(Residual prediction deviation, RPD)评价模型的预测能力[38],根据MOUAZEN等[39]提出的阈值对RPD范围划分:RPD<1.5表示模型预测能力差,模型不可靠;1.5≤RPD<2表示模型的预测能力可接受,可进行粗略的估算;RPD≥2表示模型稳定预测能力好,可以采用该模型对变量进行预测。
利用R2、RMSE和RPD评价各光谱变量对LAI的估算能力,但这3个指标只是评价了样本总体上与实测值的相关和偏离程度,为了评价光谱变量对所有处理中每一种处理的拟合程度,采用绝对积分差值(Bias of square integral,BSI)进行评价,BSI越小表示各处理的LAI值估算模型一致,反之则表示各处理LAI不适合混合估算[24],计算式为
(1)
f1、f2——两处理LAI的拟合函数
[x1,x2]——积分区间,两处理的LAI估算模型在光谱变量下形成的共同区间的下限与上限
图1为高肥水平下不同覆膜处理春玉米各生育期LAI的变化特征,图中相同小写字母表示在P=0.05水平上差异不显著。从图1可以看出,从苗期到抽雄期春玉米LAI持续增长,到拔节期各覆膜处理春玉米LAI增长速度最大,达0.12 cm2/(cm2·d)左右。到抽雄期各覆膜处理LAI值达最大,分别为4.58、5.16、4.52 cm2/cm2,进入灌浆期后玉米底部叶片变黄脱落,LAI逐渐减小。覆膜处理LAI在抽雄期前显著高于无覆膜小区,但可降解膜和普通透明膜之间差异不显著,在吐丝期到灌浆期期间3个处理玉米LAI间差异不显著,到灌浆期后M0和M2处理之间无显著性差异,但与M1处理之间差异显著。
总体上覆膜与无覆膜处理之间春玉米LAI差异随生育期呈先减小后增大的趋势,苗期最大,M2、M1比M0处理LAI分别增加90.54%和78.24%,到成熟期,M0处理相较于M1、M2处理LAI分别增加9.8%和7.2%。由此可见,覆膜对春玉米LAI的影响主要在抽雄期之前,覆膜增加了生育初期土壤温度,玉米提前出苗,使其峰值到达时间明显早于不覆膜,提前了约14 d,随着气温升高,无覆膜处理迅速生长,导致覆膜与无覆膜处理之间差异减小,之后覆膜处理由于生育期的提前,叶片开始衰败,无覆膜处理LAI高于覆膜处理,各处理之间LAI差异开始增大。整体上看,可降解膜在玉米生育初期的增温效果更好。
图2为各覆膜处理在高肥水平下不同生育期春玉米的冠层光谱曲线。由图2可见,3个覆膜处理的冠层反射率在不同生育期的变化趋势一致,在可见光(波段380~760 nm)范围内,由于叶绿素的强烈吸收,反射率较低,各生育期差异不大;在近红外反射平台(波段760~1 000 nm),反射率均是随生育期的推进呈先增大后下降的趋势,在抽雄期达到最大。与LAI在整个生育期内变化趋势一致,在抽雄期春玉米的LAI也达到最大值。
为了确定各覆膜之间冠层光谱曲线差异,计算了生育初、中、末期同一时间同一追氮水平下不同膜处理(M2、M1、M0)冠层反射率的比值(用百分数表示[40-41]),结果如图3所示。各覆膜处理间冠层反射率的差异由大到小为生育中期、生育末期、生育初期,覆膜增加作物生长速率使其生育期提前,各覆膜间玉米冠层反射率的差异主要是由生育期不同所致。由图2可知,不同覆膜处理冠层反射率均在抽雄期达到最大值,在7月11日,覆膜处理处于抽雄期而无覆膜处理处于大喇叭口期,因此反射率相差最大;在8月12日,覆膜处理处于抽雄灌浆期,叶片已呈减小的趋势,叶片结构等也发生较大的变化,冠层反射率与无覆膜处理差异明显;而6月18日各处理均处于苗期,虽然覆膜处理LAI显著高于无覆膜处理,但苗期各处理的冠层反射率均较低,所以各覆膜玉米冠层反射率间的差异反而小于生育中后期。另外可以看出,各覆膜处理冠层反射率曲线的形状差异主要在波长530、670 nm和波段680~760 nm处,表明覆膜主要影响植物对绿光和红光的吸收和利用。
2.2.1冠层反射率及一阶微分光谱与LAI
为了进一步探究春玉米LAI与原始光谱和一阶微分光谱的定量关系,利用全生育期监测的各覆膜处理的LAI和对应的光谱数据,进行Pearson相关性分析,相关性分析结果见图4。各覆膜处理对应的采用样本数分别为53、33、60。选择显著性水平0.05、0.01进行显著性检验。
由图4可见,M0处理在波段707~1 000 nm范围内光谱反射率与LAI呈显著正相关(P≤0.01),在波段760~800 nm范围内反射率与LAI间相关系数的平均值为0.80,为LAI的敏感波段。M1处理在波段690~720 nm范围内反射率与LAI间的相关系数绝对值均大于0.50,其相关系数的均值为-0.52,在波段400~500 nm范围内相关系数平均值为-0.52,这两个波段范围反射率与LAI均显著相关,但波段350~500 nm范围内的光谱反射率易受土壤背景的影响[11],因此选择波段690~720 nm作为M1处理LAI的敏感波段。M2处理在波段700~725 nm范围内LAI与光谱反射率呈显著相关,其相关系数绝对值均大于0.50。所有处理LAI与波段690~715 nm范围内反射率呈显著相关,其相关系数平均值为-0.35。一阶微分光谱与LAI间的相关系数在一些波段处高于原始光谱与LAI的相关系数,这是由于求导过程去除了背景等噪声对原始光谱的影响,与刘露等[29]、高雨茜[42]得到的结论相同。其中M0、M1、M2处理和最大相关一阶微分光谱波长分别为742、773、546 nm的相关系数绝对值均在0.68以上。
综上各覆膜处理LAI的敏感原始光谱波段均位于近红外波段范围内,但覆膜处理LAI的敏感波段往短波方向移动。
2.2.2植被指数与LAI
为了确定覆膜对植被指数的影响,进行了方差分析,并基于参数η2进行覆膜因子影响植被指数程度的评估,该参数反映了因变量(植被指数)的变异被自变量(覆膜类型)解释的比例[11]。分析结果见表2。从表2可以看出,各覆膜处理下各植被指数之间没有显著性差异(ARVI除外),另外各植被指数对应的η2都小于0.05,说明覆膜因子对各植被指数的影响较小。
表2 覆膜类型因子对春玉米全生育期植被指数值的影响
在方差分析基础上,进行了植被指数与LAI的相关性分析,结果见表3。从整体上看,M0处理LAI与表中所选植被指数均达到了0.01极显著水平;M1处理LAI与RVI、NDVI、RVI/NDVI、ARVI达0.01极显著水平,与MSAVI达0.05显著水平;M2处理相对于M0、M1处理与各个植被指数相关性弱,仅与RVI、NDVI、RVI/NDVI达0.05显著水平;所有处理的LAI与各个植被指数均达0.01极显著水平。从单个植被指数看,RVI、NDVI、RVI/NDVI均达到0.05显著水平,普适性较好。根据相关系数绝对值最大为原则确定M0、M1、M2、所有处理LAI的敏感植被指数分别为MSAVI、NDVI、RVI、NDVI。
表3 植被指数与各覆膜处理LAI之间的相关系数
2.2.3春玉米LAI的特征指数
经对比原始光谱计算的3个指数最大R2相较于一阶微分光谱小,因此本节只展示一阶微分光谱波段400~1 000 nm范围内R2的等高线图。如图5所示,各覆膜处理全生育期春玉米LAI与3个指数的R2较大区域分布一致,集中在波段700~800 nm范围内,其中M0处理计算的3个指数与LAI间R2最大值对应的指数分别为:NI(691,741)、RI(691,749)、AI(715,719),其R2分别为0.90、0.91、0.90。M1处理计算的3个指数与LAI间R2最大值对应的指数分别为:NI(698,755)、RI(728,754)、AI(618,754),其R2分别为0.85、0.87、0.84。M2处理计算的3个指数与LAI间R2最大值对应的指数分别为:NI(723,731)、RI(716,741)、AI(715,719),R2分别为0.69、0.73、0.78。所有处理计算的3个指数与LAI间R2最大值对应的指数分别为:NI(722,731)、RI(722,731)、AI(641,751),其R2分别为0.76、0.77、0.67。
从整体上看与各覆膜处理春玉米LAI间R2最大值对应的特征指数均由波段600~760 nm的一阶微分光谱组成,即大多集中在“红光”范围,这也证明了“红光”可广泛用于对作物长势、生理生态指标的监测[43-45],相较于敏感波段和表1选取的植被指数,特征指数可更好地用于LAI的估算。
以敏感光谱为自变量,LAI为因变量建立各覆膜处理下基于统计方法的LAI单变量估算模型,2/3的数据用于模型建立,1/3的数据进行模型验证,并以R2、RMSE、RPD为模型评价标准,确定各覆膜处理的最佳估算模型,模型拟合和验证结果见表4。从表4可知,M0处理的LAI与各敏感光谱变量建立的估算模型估算精度较高,其拟合和验证的R2在0.76~0.94之间,RPD在2.10~5.01之间,均大于2,说明模型都较为稳定,具有较好的预测能力;M1处理的LAI与特征指数NI(698,755)建立的估算模型能较为精确地估算LAI,其RPD为2.15,其余模型RPD均小于1.5;M2处理LAI与特征指数AI(715,719)建立的模型估算精度理想,RPD为2.11,其余模型均不稳定,预测能力差。
表4 各覆膜处理LAI基于敏感光谱变量的统计估算模型结果
从各敏感光谱变量分析,以特征指数建立的LAI估算模型相较于一阶微分光谱单波段、植被指数模型估算精度有所提高,郑煜[46]也得出相同的结论;各覆膜之间,M0处理与各敏感光谱变量的拟合效果最好,地膜覆盖的土壤表面与无覆膜有着不同的光谱特性。采集的光谱反射信息表现为:地膜覆盖的反演目标有效光谱信息较无覆膜处理少,即存在杂散射等情况,噪声较强,降低了信噪比,影响了光谱相关性,因此覆膜处理LAI反演精度较无覆膜低。姚志华等[47]研究发现覆膜会降低无人机对土壤盐分反演的精度,武江星[48]也发现覆膜处理高粱土壤呼吸速率高光谱反演的精度较无覆膜低。
所有处理全生育期LAI的估算模型拟合和验证的结果见表5和图6,以特征指数NI(722,731)建立的估算模型较为稳定,估算精度较高。进一步计算了采用特征指数NI(722,731)建立的拟合方程y=-50.15x2-23.54x+1.79单独估算各覆膜处理的LAI,其拟合的RMSE分别为0.68、0.68、0.66 cm2/cm2,验证的RMSE分别为0.36、0.71、0.55 cm2/cm2,可以看出,3个处理拟合和验证的RMSE相差不大,说明各处理的LAI值均可用特征指数NI(722,731)建立的估算模型反演,并且反演效果相差不大,精度较高。
表5 所有处理LAI基于敏感光谱变量的统计估算模型结果
为了探究特征指数NI(722,731)对所有处理中每一种处理的拟合程度,采用绝对积分差(BSI)进行评价。经计算M0和M1处理组合间BSI值为0.13,M0和M2处理组合间BSI值为0.09,M1和M2处理组合间BSI值为0.18。各处理之间的组合BSI均小于0.2,说明以特征指数NI(722,731)建立的LAI估算模型对所有覆膜处理均有效,可降低覆膜对LAI的估算影响[24]。
(1)研究了不同覆膜处理下LAI和冠层光谱形状的变化。结果表明,覆膜对春玉米LAI的影响主要在抽雄期之前,相同施肥水平下覆膜与无覆膜处理之间LAI的差异随生育期的推进呈先减小、后增大的趋势,苗期差异最大。各覆膜处理间冠层反射率的差异由大到小为生育中期、生育末期、生育初期,覆膜对春玉米冠层光谱曲线产生一定的影响,且反映在光谱的某些波段上,在波长530、670 nm,波段680~760 nm处,各覆膜处理下冠层反射率曲线形状有明显的变化。
(2)研究了不同覆膜处理下LAI的敏感光谱变量。结果表明,M0、M1、M2处理LAI的原始光谱敏感波段分别为760~800 nm、690~720 nm、700~725 nm,一阶微分光谱敏感波长分别为742、773、546 nm,敏感植被指数分别为MSAVI、NDVI、RVI;一阶微分光谱任意两波段组合的3个指数与各覆膜LAI决定系数R2较大的区域集中在波段700~800 nm范围内,其中AI(715,719)、NI(698,755)、AI(715,719)与M0、M1、M2处理的LAI相关性强,R2分别为0.90、0.85、0.78。
(3)通过对不同LAI估算模型的比较得到最佳估算模型。以特征指数为自变量建立的估算模型稳定性和精度较好,能够对各覆膜处理下玉米LAI进行准确估算,其中无覆膜处理LAI的反演精度高于覆膜处理。
(4)对于3种覆膜处理下的LAI混合反演,以特征指数NI(722,731)为自变量建立的估算模型表现了反演的优越性,模型较为稳定,估算精度较高,可降低覆膜处理对LAI估算的影响。