试验样品为河北省沧州市吴桥县中国农业大学吴桥试验站种植的19个玉米品种: 中单909、 蠡玉31、 沃玉964、 纪元1号、 隆平208、 冀丰223、 登海618、 农大108、 浚单20、 洰丰339、 纪元101、 登海605、 掖单13、 先玉335、 京单28、 隆平206、 蠡玉37、 郑单958、 圣瑞999。 试验品种在同一地块内分不同的小区种植, 种植密度为5 000株· 亩^-1, 每个小区种植9行, 行长5 m, 行距0.66 m, 株距22.2 cm。 玉米生长过程中的田间管理措施按照当地的大田管理方式实行。

在玉米生长至灌浆期和蜡熟期, 分别到田间采集玉米的茎秆带回实验室测量。 在每个小区内随机挑选玉米茎秆长势均匀的6株植株, 每个茎秆样品截取基部的一节茎节, 先进行高光谱图像拍摄, 再用茎秆强度测定仪测试茎秆的穿侧强度和折断力。 高光谱成像采用美国SOC(Surface Optics Corporation)公司的SOC710VP型号高光谱成像光谱仪, 在光学暗箱内进行成像拍摄。 仪器的光谱波长范围374~1 038 nm, 光谱分辨率为4.68 nm, 共128个波段, 成像的图像分辨率为696× 520。 拍摄时使用光源为4个100 W的卤素钕灯, 光谱相机布置在茎秆的正上方, 镜头到拍摄样本的距离为75 cm, 茎秆样本和灰度标准板一同拍摄, 便于后期进行光谱图像反射率校正。 在完成光谱图像拍摄后, 对茎秆目标区域进行图像分割并提取平均光谱曲线用于后续分析。 使用YYD-1型茎秆强度测定仪测量茎秆的穿刺强度和折断力, 测量时设定茎秆的悬空间距为8 cm, 力量的加载位置为茎秆中间位置。

首先通过观察剔除掉因试验误差造成的明显偏离数据中心的样本曲线以及茎秆强度数据缺失的样本, 然后对每一个光谱数据集使用偏最小二乘法和蒙特卡洛交叉验证法进行循环验证, 计算各样本预测误差的均值和方差, 将预测误差均值或方差偏离数据中心较远的样本作为异常值剔除。 然后对每条光谱样本曲线使用多元散射校正(multiplicative scatter correction, MSC)方法进行预处理, 目的是在保留光谱中与化学成分有关信息的同时消除不同样本间因散射影响所导致的基线平移和偏移现象^[6]。 最后对数据集的各波段变量进行标准化处理, 以突出光谱间的特征差异。

选择使用主成分回归(principle component regression, PCR)和偏最小二乘回归(partial least squares regression, PLSR)的方法, 建立不同时期茎秆光谱数据与穿刺强度、 折断力的回归模型。 主成分分析是通过线性变换的方式, 将样本矩阵沿着协方差最大的方向由高维空间向低维空间投影, 得到的各主成分之间相互正交, 可以有效解决光谱数据的多重共线性问题^{[7, 8]}。 主成分回归首先对光谱数据进行主成分分析, 所得到的前f个主成分向量代替原光谱变量, 再进行多元线性回归即得到主成分回归模型。 主成分回归的特点是适合处理共线性特征, 只用几个主成分就能包含原始数据的主要信息^{[9, 10]}。 偏最小二乘回归是光谱数据分析中常用的一种多元统计分析方法, 它在建模时需要同时考虑光谱特征和响应变量的信息, 探究两者之间的线性关系。 PLSR模型通过将光谱特征和响应变量分别投影到相应的低维空间得到相互正交的特征向量, 再在各自的特征向量之间进行线性回归, 并通过包裹式特征选择方法选择合适的特征投影个数, 使模型效果最优, 因此由偏最小二乘回归选择出的特征强调了光谱特征对响应变量的解释和预测作用^{[11, 12, 13]}。

为了保证模型有良好的泛化性能, 使用10折交叉验证法对模型效果进行验证。 使用决定系数(coefficient of determination, R²), 均方根误差(root mean square error, RMSE)和剩余预测残差(residual predictive deviation, RPD)作为模型评价标准, 决定系数和剩余预测残差越高、 均方根误差越小, 说明模型效果越好^{[14, 15, 16]}。

灌浆期和蜡熟期测得的茎秆基部穿刺强度和折断力数据如表1所示。 茎秆在灌浆期和蜡熟期的穿刺强度平均值、 中位数和标准差都比较接近, 说明不同时期的茎秆穿刺强度没有明显变化, 数据内的变异程度也接近。 而蜡熟期的茎秆折断力数据的各项指标均明显高于灌浆期数据, 这主要是后期茎秆脱水、 内部加速纤维化的结果。 数据结果说明茎秆纤维化对折断力指标的影响比较显著, 对穿刺强度影响较小。

表1

Table 1

表1(Table 1)

表1 不同时期茎秆强度数据变化 Table 1 Changes in stalk strength data in different periods 项目 穿刺强度/N 折断力/N 灌浆期 蜡熟期 灌浆期 蜡熟期 平均值 49.20 49.65 350.68 427.92 中位数 49.40 48.65 338.70 409.65 最大值 79.70 86.50 695.60 710.80 最小值 29.90 28.90 168.00 219.70 标准差 9.05 9.26 89.35 102.29

表1 不同时期茎秆强度数据变化 Table 1 Changes in stalk strength data in different periods

首先对光谱图像进行反射率校正。 从图1(a)光谱图像中提取出标准板、 茎秆、 背景区域的光谱反射强度曲线, 如图1(b)所示。 可以看出各区域的光谱反射强度有明显的差别, 选择752.98 nm处的波段图像, 设置反射强度阈值范围[300~1 000]分割出标准板区域, 而反射强度高于2 000的区域为茎秆, 各目标区域具体分割效果如图1(c)所示。 然后在各波段上提取出标准板坐标区域的平均反射强度I_std, 用I_std除以标准板的反射率R_std计算出总入射光的强度I, 再用光谱图像各波段的反射强度除以I完成图像的反射率校正, 其中标准板R_std的数值是已知的。 最后提取出茎秆区域的反射率光谱曲线, 用于对茎秆强度进行建模预测。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 目标区域提取方法 (a): 原始图像; (b): 各区域光谱曲线; (c): 阈值分割结果Fig.1 Target area extraction method (a): The original image; (b): Spectral curve of each region; (c): Threshold segmentation result

经过图像分割提取, 试验在灌浆期和蜡熟期分别采集了114个样本。 人工剔除异常样本后再对每一个光谱数据集使用偏最小二乘法和蒙特卡洛交叉验证法将预测误差均值或方差偏离数据中心较远的样本作为异常值剔除。 如图2(a— d)所示, 大部分样本预测误差的均值和方差都比较集中, 有少数样本偏离的数据中心, 则该样本预测误差的波动比较大。 将灌浆期穿刺强度(CcGjq)的第1、 10、 16、 33号样本, 灌浆期折断力(ZdGjq)第1、 30号样本, 蜡熟期穿刺强度(CcLsq)的第2、 3、 7号样本, 和蜡熟期折断力(ZdLsq)的第3、 10、 11、 28号样本剔除。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 异常样本剔除Fig.2 Abnormal sample removal

经过异常样本剔除后, 绘制两个时期茎秆的平均光谱曲线, 如图3(a)所示。 可以看出不同时期的茎秆光谱曲线变化趋势基本一致, 在679.5 nm附近都存在一个明显的波谷, 在560.9~632.8和731.9~918.6 nm范围内是波峰。 对每条样本曲线进行标准正态变换处理, 然后分析不同时期各波段反射率与茎秆穿刺强度和折断力之间的相关性, 绘制如图3(b)所示的相关系数图。 由图可知在374~479和763~994 nm范围内, 波段反射率与茎秆强度之间是正相关关系, 在499~721 nm范围内则主要是负相关关系。 图3(b)中实线是穿刺强度与波段反射率的相关系数, 在灌浆期560.9 nm处的相关性系数最高为0.48, 蜡熟期时在962 nm处的相关系数最高为0.37。 虚线显示折断力与波段反射率之间的相关关系, 在灌浆期806 nm处的相关性系数最高为0.41, 蜡熟期时在854 nm处的相关系数最高为0.55。 可以看出茎秆穿刺强度的相关系数绝对值在大多数波段上是灌浆期高于蜡熟期, 而折断力的相关系数则是蜡熟期高于灌浆期, 这应该是与玉米茎秆在蜡熟期时开始脱水和茎秆的纤维化相关。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 光谱曲线特征分析 (a): 两个时期的平均光谱曲线; (b): 两个时期各波段的相关系数曲线Fig.3 Characteristic analysis of spectral curve (a): Average spectral curve of two periods; (b): Correlation coefficient curve of each band in the two periods

2.3.1 主成分分析法

主成分分析法(principal component analysis, PCA)只需要使用方差贡献率最大的前几个主成分, 就可以包含原始数据中的主要信息。 以累计方差贡献率达到99.9%为限, 决定待选主成分数量, 选择灌浆期主成分数为16, 蜡熟期主成分数为15。 各生长期的主成分贡献率如表2所示, 由表2看出各生长期均是前2个主成分的贡献率最大, 前3个主成分的累计贡献率都已经超过90%, 说明主成分分析法对光谱特征的降维效果显著。

表2

Table 2

表2(Table 2)

表2 各时期光谱数据主成分的贡献率 Table 2 Contribution rate of principal components of spectral data in each period 主成分 灌浆期 蜡熟期 贡献率 /% 累计贡献率 /% 贡献率 /% 累计贡献率 /% PC1 67.89 67.89 78.80 78.80 PC2 17.27 85.16 10.69 89.49 PC3 5.46 90.62 4.9 94.39 PC4 4.78 95.4 2.39 96.78 PC5 1.68 97.08 1.97 98.75 PC6 1.13 98.21 0.48 99.23 PC7 0.45 98.66 0.28 99.51 PC8 0.39 99.05 0.11 99.62 PC9 0.24 99.29 0.08 99.7 PC10 0.19 99.48 0.06 99.76 PC11 0.14 99.63 0.05 99.81 PC12 0.11 99.73 0.03 99.84 PC13 0.08 99.81 0.03 99.87 PC14 0.05 99.86 0.02 99.89 PC15 0.03 99.89 0.02 99.91 PC16 0.02 99.91 - -

表2 各时期光谱数据主成分的贡献率 Table 2 Contribution rate of principal components of spectral data in each period

2.3.2 包裹式特征选择方法

使用包裹式特征选择方法选择PLSR模型的最佳建模特征。 对各生长期的光谱数据, 以5为增量分别设置特征投影个数为5~60个, 训练偏最小二乘回归模型, 并使用10折交叉验证法检验模型, 用RMSE和R²评价模型预测性能, 各数据集的训练效果如图4(a, b)所示。 从图4可以看出随着特征投影数量的增加, 各训练模型的性能逐渐变强。 灌浆期穿刺强度和折断力数据训练模型的均方根误差RMSE普遍低于蜡熟期数据模型, 决定系数R²也高于蜡熟期, 说明模型对灌浆期的数据更敏感。 看到当投影个数为30时, 各训练模型的R²均在0.9以上, RMSE也较低, 因此选择30个投影特征为建模的最少特征个数。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 投影特征个数对PLSR模型性能的影响 (a): 各穿刺强度训练模型的性能; (b): 各折断力训练模型的性能Fig.4 The influence of the number of projection features on the performance of the PLSR model (a): The performance of each puncture intensity training model; (b): The performance of each breaking force training model

将灌浆期和蜡熟期时PCR模型和PLSR模型的最佳参数分别带入相应的光谱数据集建立预测模型, 各模型的预测效果如图5(a— d)所示, 各模型评价指标见表3。 偏最小二乘模型在不同生长期对穿刺强度和折断力的模型预测指标均优于主成分回归模型。 PCR模型对各茎秆强度指标预测模型的决定系数R²在0.36~0.64之间, 模型的预测值与真实值相关度一般, 剩余预测残差RPD在1.25~1.67之间, 说明模型的预测效果较差。 PLSR模型对各茎秆强度指标预测模型的决定系数R²在0.88~0.94之间, 模型的预测值与真实值密切相关, 剩余预测残差RPD在2.90~3.93之间, 说明模型有很好的预测效果, 可以用于定量分析预测。 通过对比发现, PCR模型和PLSR模型在不同生长期对茎秆折断力的预测效果均优于对茎秆穿刺强度的预测效果, 说明光谱数据对茎秆折断力的特征反映比对穿刺强度的反映更加敏感; 且两类模型对蜡熟期穿刺强度的预测效果均为各自模型中最低, 说明玉米在不同生长期光谱数据的茎秆穿刺强度特征在减弱。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 训练集模型预测效果 (a): 灌浆期-穿刺强度; (b): 灌浆期-折断力; (c): 蜡熟期-穿刺强度; (d): 蜡熟期-折断力Fig.5 The prediction effect of the training set model (a): Puncture strengthduring grouting period; (b): Breaking force during grouting period; (c): Puncture strength during wax maturity period; (d): Breaking force during wax maturity period

表3

Table 3

表3(Table 3)

表3 各模型性能评价 Table 3 Performance evaluation of each model PCR RMSE R2 RPD 灌浆期 穿刺强度 6.10 0.48 1.38 折断力 49.91 0.64 1.67 蜡熟期 穿刺强度 6.83 0.36 1.25 折断力 71.36 0.49 1.40 PLSR RMSE R2 RPD 灌浆期 穿刺强度 2.43 0.92 3.47 折断力 21.73 0.93 3.84 蜡熟期 穿刺强度 2.95 0.88 2.90 折断力 25.45 0.94 3.93

表3 各模型性能评价 Table 3 Performance evaluation of each model