样本采集地位于中国海南省儋州市, 108°56—109°46E, 19°11—19°52N。 该区域属于亚热带气候, 平均气温26.5°。 本研究的对象是中国热带农业科学院试验场广泛种植的“RY-7-33-97”的叶片。 为确保样品丰富度, 在试验场不同区域随机采集成熟、 完整、 无病害的中间小叶119片作为试验样品, 在实验室进行标注和光谱测量, 然后使用化学测量方法得到对应样本的氮含量。 使用python 3.80软件处理光谱数据和建立回归模型。

1.2.1 光谱测量与数据预处理

使用Spectrometer GaiaField-F-N17E光谱仪测量叶片样本的反射率图像, 波长采集范围为840~1 680 nm, 光谱采样间隔平均为3.3 nm。 待仪器效正完成后, 每个样本依次放入移动平台测量叶片正反面的光谱反射率图像, 平均后输出。 由于前30个采集波段信噪比极低, 因此只保留224个波段, 对应的波长范围为942~1 680 nm。 最终获得了119个三维光谱数据。 图1为五幅叶面光谱反射率的灰度图像。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 No.1—5样品在波长1 400 nm时的光谱反射率灰度图像Fig.1 Grayscale images of spectral reflectance of No.1—5 samples image at wavelength 1 400 nm

对每一个波长的反射率图像进一步求平均, 以获取样本的原始平均光谱曲线。 但在原始光谱的采集过程中, 会受到散射和环境中的随机噪声的干扰^[7], 因此, 分别使用了三种广泛使用的预处理方法消除散射和噪声对原始数据中的干扰, 包括多元散射校正(multiplicative scatter correction, MSC), 对数倒数谱变换(Absorbance)和标准正态变换(standard normal variate, SNV)的三种不同方法, 同时分析和对比它们对模型预测精度的提升效果。 经过预处理后的光谱曲线如图2所示。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 不同预处理方法转换后的平均光谱曲线 (a): MSC; (b): Absorbance; (c): SNV; (d): Raw dataFig.2 Transformed and averaged spectra by different pretreatment methods (a): MSC; (b): Absorbance; (c): SNV; (d): Raw data

为了排除异常样本的干扰, 使用了主成分分析(PCA)对每个样本进行降维处理, 将它们投影到三维空间, 再利用隔离森林算法(isolation forest algorithm, IFA)根据每个样本在三维空间的分布情况进行区分。 区分结果如图3所示, 110个正常样本, 9个为异常样本。 IFA是一种常用的异常值检测算法^[8], 该算法通过随机选择特征, 然后随机选择特征的分割值, 递归地生成数据集的分区。 由于和数据集中正常的点相比, 异常值所需的随机分区更少。 因此区分准则为: 异常值是多颗树中路径较短的点, 路径长度计算方法如式(1)和式(2)。

c(n)=2H(n-1)-2n-1n(1)

H(*)=ln(*)+0.5772156649(2)

式中, H(i)为调和数。 c(n)为给定n个样本时的路径长度的平均值, 用于标准化的路径长度h(x)。 因此, 样本的异常得分定义为

s(x, n)=2-E[h(x)]c(n)(3)

式(3)中, E[h(x)]为样本x在一批孤立树中的路径长度的期望。 当E[h(x)]→0时, s→1, 该样本x被判为异常点。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 IFA区分的异常值和正常点的散点图Fig.3 Outliers and normal samples distinguished by IFA

1.2.2 叶片氮含量的化学分析

目前, 凯氏定氮法已广泛用于检测食品、 农作物的氮含量^{[9, 10]}。 利用凯氏定氮法测量剩余110个样品中的叶面氮浓度(LNC), 将得到的数值作为每个叶片样本的标签。 统计直方图为图4所示。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 样品氮浓度统计直方图 横坐标是前后刻度构成叶片氮含量区间, 纵坐标表示频率, 即该浓度下的样本个数Fig.4 Statistical histogram of sample nitrogen concentration The ordinate represents frequency, consists of the times of occurrences for x-axis; Abscissa is the small interval of formed by proceeding and following nitrogen content value

1.3.1 波段变量选择

为了消除冗余特征和共线特征对预测模型参数的有害影响并对模型精度进行有针对性的二次优化, 采用了基于IMF优化框架的特征选择方法。 改进后的MF框架是高效的一体式框架, 由全局筛选和局部优化组成: CARS作为全局变量筛选算法, 用于消除光谱特征间的冗余性和共线性并选出在所有波段中表现好的波段作为初选子集。 贪心算法近邻搜索(NNS)作为局部优化算法, 利用初选子集的邻近波段对已选出的波段进行更新, 进行二次优化, 提升模型精度。 该框架的核心思想是对光谱波段从全局到局部, 整体到个体的比较与筛选, 达到优化精度的目的。 具体的算法框架如图5所示。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 算法流程图Fig.5 Method flowchart

1.3.2 全局筛选

从整体角度出发, 全局筛选的目的是从全波段中选出能提供更高预测精度的变量, 作为初选子集。 不同于原始的模因算法(memtic algorithm, MA), 使用竞争性自适应重加权采样(CARS)取代基于基因传播和演化的变量筛选过程, 作为新的全局优化算法。 其原因为, 相比于传统的基因传播过程, CARS-PLSR能够更快的收敛并且受初始化的影响更小。 竞争性自适应重加权采样法(CARS)是一种在光谱波段筛选中广泛使用的启发式方法, 其与蒙特卡罗采样(monte carlo, MC)和局部最小二乘回归(partial least squares regression, PLSR)相结合。 每次采样时, 具有更大的回归系数绝对权值的光谱波段将被保留为新子集。 大量的研究已证明将CARS作为变量选择算法对模型精度提升十分明显且不需要花费大量计算时间^[11]。

1.3.3 局部优化

局部优化的过程可以看作是对选出的波段个体的进一步优化。 由于邻近的光谱波段之间存在着高度的共线性和相关性, 通过引入了皮尔逊系数(person coefficient, PC)以描述邻近波段之间的相关性强度, 其数值越接近于1说明变量之间相关性越强。 图6证实了光谱波长越接近, 两者相关程度会更高。 分析可知两个邻近的光谱波段变量拥有最强的相关性和相似性。 所以使用相邻波段对初选波段进行逐个替换, 有一定的机率能产生新的表现更好的模型。 因此, 使用了一种近邻搜索算法(NNS), 有方向性和针对性地对波段子集进行二次优化, 而不需要消耗大规模的计算资源就能明显地提高模型精度, 是一种代价小收益高的贪心算法。 近邻搜索的具体步骤如图7所示。 在获得全局算法选出的子集后, 再对原始子集中的逐个光谱波段个体进行更新。 如果更新后的模型性能提高, 保留本次更新, 否则不进行更新, 直接进入下一个循环。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 光谱子波段之间的相关性程度热力图Fig.6 A heat map of the degree of correlation between spectral bands

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 局部优化算法流程图Fig.7 Flowchart of local optimization algorithm

如何定义邻近域范围是使用NNS算法的核心问题。 选择较大的范围对初选子集进行优化可找到更优秀的特征子集解, 但会导致大规模计算和时间消耗, 甚至出现不收敛的现象。 因此, 引入了一种新的参数K, 用于构造邻近区域[-K, K]。 分析图6可知, 随着波长间隔的增加, 两个波长特征之间的系数减小。 当两个波段的索引间隔超过3时, 两者的相关系数小于等于0.5, 相关程度明显下降。 因此令K等于3, 即, 定义邻近区域由左右三个波段(共六个)组合而成。 算法的详细流程如图7所示。

叶片样本的光谱反射率数据集被分为校准集(calibration set)和预测集(prediction set), 其中预测集占比0.3, 校准集占比0.7。 校准数据集一般用于调整模型的参数使得在校准集上表现最好, 并使用预测数据集进行测试模型的性能。 校准集与预测集分布尽可能的接近, 可以更真实地反映模型预测能力。 因此, 使用了以下的数据集拆分方式, 如表1所示。

表1

Table 1

表1(Table 1)

表1 数据集划分与氮含量分布 Table 1 Data split and LNC distribution Data Set (数据集划分) Size (样本数量) LNC Range /% LNC Mean /% CV /% Total 110 2.02~4.23 3.11 15.28 Calibration Set 77 2.02~4.23 3.03 15.77 Prediction Set 33 2.21~4.08 3.10 13.97

表1 数据集划分与氮含量分布 Table 1 Data split and LNC distribution

模型精度应该根据校准集和预测集来评估。 均方根误差(RMSE)和决定系数(R²)用于评价预测模型的性能, 两个评价指标的方程列在式(4)和式(5)中。

RMSE=∑i=1n(yi-ŷi)2n(4)

R2=∑i=1n(yi-y¯i)(ŷi-ŷ¯i)∑i=1n(yi-y¯i)2∑i=1n(ŷi-ŷ¯i)22(5)

式中, y_i是第i个样本的预测值;ŷi是第i个样本的真实值;y¯i是平均预测值;ŷ¯i为平均测量值; n是样本的总数。 一般来说, 理想模型应该表现出较高的Rp2值和较低的RMSE_p值。

数据预处理方法和特征提取算法都会对预测模型的精度有影响, 需要分别评估它们对预测结果的影响。 因此, 使用RMSE_c和Rc2用于评价模型在校准集上的表现, 使用RMSE_p和Rp2考察模型在预测集上的表现。

2.3.1 光谱预处理方法评估

当进一步可视化异常值消除的效果时, 经IFA算法去除异常点后的模型结果比没有异常点消除的模型拥有更符合预期的精度, 如表2所示, RMSE_p提高了42.2%,Rp2明显提高了14.9%。 分析表2可知, 对比其他的预处理方法, 使用MSC对原始光谱进行转换后的结果最佳, 预测集的Rp2平均值等于0.841, RMSE_p为0.188。Rp2比未经过预处理的光谱建立的预测模型的精度提高了0.104。 因此, 选择使用MSC算法对原始数据进行转换, 然后对转换后的数据建立预测模型。

表2

Table 2

表2(Table 2)

表2 异常点消除对预测结果的影响 Table 2 The effect of outliers elimination on prediction results 异常点 消除方法 预处理方法 校准集 预测集 RMSEc Rc2 RMSEp Rp2 Non Raw data 0.512 0.742 0.388 0.607 SNV 0.500 0.779 0.368 0.611 Log(1/r) 0.478 0.794 0.296 0.707 MSC 0.415 0.821 0.270 0.724 IFA Raw data 0.421 0.818 0.256 0.736 SNV 0.311 0.872 0.242 0.765 Log(1/r) 0.262 0.913 0.200 0.824 MSC 0.236 0.908 0.188 0.841

表2 异常点消除对预测结果的影响 Table 2 The effect of outliers elimination on prediction results

2.3.2 光谱波段选择算法评估

由图8可以分析出, 对比使用MSC转换后的光谱和不同波段选择算法建立的预测模型, 基于改进后的MF的CARS-NNS的特征选择方法的预测效果最好, 拥有最低的RMSE_p和最高的Rp2, 分别为0.116和0.951。 同时, 箱线图9更直观地反映了不同特征选择算法的预测决定系数值(Rp2), 相比于分别使用CARS和MA算法建立模型, 使用由两者组成的IMF框架得到的模型评价指标显著提升。 分析表3的最后两列, 在使用相同的预处理方法的情况下, 基于MF的CARS-NNS被证明是最佳的波段选择方法。 综上所述, 由MSC预处理方法和IMF框架训练后生成的预测模型是橡胶叶片LNC定量估算的最佳选择。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 MSC预处理对不同算法的预测结果评估的散点图 (a): Full bands; (b): SPA; (c): CARS; (d): IMF; (e): Basic MAFig.8 The scatter diagrams of the prediction evaluation for different selection approaches by using the same MSC method (a): Full bands; (b): SPA; (c): CARS; (d): MF; (e): Basic MA

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 不同波段选择算法的预测结果箱线图Fig.9 The box diagram of prediction results using different band selection algorithms

表3

Table 3

表3(Table 3)

表3 所有模型的预测结果 Table 3 Result evaluation of all prediction models in our work 预处理方法 波段选择算法 校准集 预测集 RMSEc Rc2 RMSEp Rp2 FB-PLSR 0.442 0.771 0.301 0.662 SPA-PLSR 0.253 0.834 0.265 0.738 Raw data CARS-PLSR 0.237 0.872 0.207 0.821 IMF(ours) 0.211 0.911 0.204 0.825 MA 0.516 0.704 0.304 0.636 FB-PLSR 0.226 0.825 0.281 0.677 SPA-PLSR 0.15 0.909 0.228 0.802 SNV CARS-PLSR 0.127 0.934 0.202 0.829 IMF(ours) 0.117 0.971 0.161 0.895 MA 0.489 0.722 0.340 0.626 FB-PLSR 0.155 0.901 0.241 0.768 SPA-PLSR 0.163 0.917 0.187 0.856 Log(1/r) CARS-PLSR 0.121 0.966 0.167 0.886 IMF(ours) 0.041 0.991 0.131 0.934 MA 0.384 0.790 0.278 0.680 FB-PLSR 0.222 0.858 0.224 0.789 SPA-PLSR 0.127 0.934 0.190 0.850 MSC CARS-PLSR 0.129 0.955 0.143 0.929 IMF(ours) 0.043 0.987 0.116 0.951 MA 0.213 0.808 0.267 0.681

表3 所有模型的预测结果 Table 3 Result evaluation of all prediction models in our work

不同方法提取的敏感波段的数量和分布存在明显差异, 选取的波段数量越少越有利于减少冗余信息, 越有利于降低模型的复杂度和运算时间。 正如图10所示, 基于IMF的CARS-NNS有效地提取了28个波段用于后续的回归建模, 是所有波段提取算法中数量最少的, SPA算法提取了34个波段。 而MA算法选取的特征波段数量为78个, 压缩效果不太理想, 如图10(d)。 分析和比较所选光谱波段的分布位置可以看出: 由不同算法训练出的模型所选取的敏感带不同。 使用MA算法所得到的波段子集数量比其他三种特征选取算法更多, 如图10(d)所示。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 不同算法对敏感光谱波段提取的结果 (a): SPA; (b): CARS; (c): IMF-CARS-NNS; (d): MAFig.10 Distribution of sensitive spectral wavebands extracted by different algorithms (a): SPA; (b): CARS; (c): IMF-CARS-NNS; (d): MA

使用基于IMF框架的CARS-NNS的特征搜索方式, 其优势可以被归纳为以下几点: (1), IMF框架的局部优化是对全局算法CARS选出的波段子集进行了更精细, 更个体化地再次优化, 由此可以保证模型精度果不逊于传统的CARS。 (2) 局部优化的NNS是一种代价小收获高的贪心算法, 能够很好地利用光谱波段之间高相似性建立邻域范围, 针对初选子集进行高效地二次优化, 提升模型精度。 (3) IMF框架由两个子模块组成, 因此, 它能一次性完成消除冗余特征和有针对性的二次优化的两个任务, 是一种成熟的一体式框架。 最后, 本工作所提出的基于IMF框架的波段选择算法需要进行调整的超参数比较少, 只有全局优化中的采样次数N和局部优化中的邻近范围K。 超参数的数量也是衡量算法优劣的一个重要部分, 越多的超参数意味着想要模型到达最佳预测效果, 需要花费更多的时间去对超参数一一调整^[12], 这些时间显然是研究者不想要花费的。

研究发现对PLSR的依赖度非常高。 主要原因是PLSR是一种非常适合在样本少, 特征多的情况下的回归建模的算法。 在耗时和算法鲁棒性的方面, 本方法比基于遗传算法(genetic algorithm, GA)的传统MA算法要好得多。 可以从表3中得出此结论。 在使用传统的MA算法建立预测模型的时候, MA算法需要花费更多的时间找到优秀特征集合, 更多的没有方向的遍历和搜索。 与此相比, 使用CARS作为全局筛选算法, 能够更快地找到全局优秀特征集合, 还能保证预测模型的精度下限。 此外, 由于传统的MA算法拥有大量的超参数, 例如, 初始化种群数量, 交叉率和变异系数等等参数, 如果初始化的参数不理想, 会导致结果不可信, 具体体现在使用MA进行预测时决定系数R²会小于0。

尽管本研究有了一些收获, 但仍然存在一些局限性。 在采集的时候无法判断样本的氮含量状况, 直接导致了叶片的氮含量的方差很大, 因此造成了对数据集划分产生不利影响。 如果不经过人为的调整, 校准集和预测集之间的分布存在明显的差距, 必须人为地进行调整才能使得两个集合分布情况尽可能的相同。 此外, 全局优化对结果的影响大于局部优化。 具体地, 基于IMF框架的混合算法的速度和精度直接受CARS结果精度的影响, 因为该算法具有随机性。 今后的研究将着重于直接通过高光谱反射率图像(hyper spectral image, HSI)提高预测精度。 虽然基于高维数据的回归和分类是一项极具挑战性的任务, 但高维数据中包含了丰富的空-谱信息, 有利于训练出更好的模型^[13]。 毫无疑问, 今后将在深度学习领域涉及更强的模型, 以处理大规模的高维数据。