黄山市黄山区仙源弦歌油茶专业合作社种植地(30° 19'N, 118° 9'E, 见图1)研究区地处皖南山区中部, 境内地形地貌多样, 以中、 低山和丘陵为主。 其气候带处于北亚热带和中亚热带过渡区, 属亚热带季风湿润气候。 夏季持续时间长、 温度高, 冬季较冷。 全年雨水充沛。 试验于2021年10月28日至11月2日期间开展, 正值油茶果实成熟收获期, 除了进行锄草, 调查时间的前7个月内无其他经营措施, 林中病虫侵害情况较轻。 研究区内油茶样本具有明显的大小年交替现象, 采集时正处在大年丰产阶段。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 研究区位置图Fig.1 Location of the study area

1.2.1 油茶叶片碳氮数据获取

供试油茶为2012年栽植的长林系列, 样地内随机布点取样120株, 包含三个长林品种(长林27号、 长林40号和长林53号)。 每株油茶采摘冠层中上部的东南西北四个方位, 每个方位2片, 共16片叶子, 带回实验室于105 ℃下杀青30 min, 后在60 ℃下烘干至恒重。 烘干后研磨并过60目筛, 再用锡囊包取2~3 mg样, 最后使用碳氢氮元素分析仪(EURO EA3000, EuroVector Inc., Italy)测定得全碳(total carbon, TC)和全氮(total nitrogen, TN), 单位统一换算成g· kg^-1。

1.2.2 冠层高光谱数据采集与处理

选择晴朗无风或微风的天气, 使用便携式野外地物光谱仪(ASD FieldSpec4 Wide-Res, Analytical Spectral Devices, Boulder, Co., USA), 在北京时间10:00— 14:00间(太阳高度角大于45° )获取油茶冠层高光谱。 仪器先预热20 min, 每株单次测定10条光谱, 除去异常光谱曲线求均值得到最终光谱曲线。 由于油茶植株高度普遍较高, 受限于光纤25° 的前视场角, 需将探头连接光纤跳线再使用伸缩杆抬至冠层正上方2 m位置处。 该设备光谱测定的全波段范围覆盖350~2 500 nm, 其中, 间隔在350~1 000 nm波段范围(visibleand near-infrared, Vis-NIR)的光谱分辨率是3 nm、 1 000~2 500 nm波段范围(short-wave infrared, SWIR)是30 nm, 全波段重采样间隔为1 nm。 除去少量信噪比低的波段(350~400 nm), 采用在400~1 000 nm范围的Vis-NIR光谱特征。 首先, 按照5 nm的重采样间隔设置原始光谱, 形成121个波段总数; 其次, 使用2阶9窗口点数的Savitzky-Golay滤波对每条高光谱曲线因设备暗电流引发的随机高频噪声进行平滑, 同时考虑云层、 地形和冠层结构等众多因素对光谱的影响, 采用平均归一化法[见式(1)]以增加光谱数据内部的可比性; 最后, 将经过上述处理的光谱重新定义为原始光谱(R)。

${{R}_{Ni}}={{R}_{i}}/\left( \frac{1}{n}\overset{n}{\mathop{\underset{i=1}{\mathop \sum }\, }}\, {{R}_{i}} \right)$(1)

式(1)中, R_N是归一化后的光谱反射率, R_i是归一化前不同波长处的光谱反射率, n是波段总数。

1.3.1 原始光谱特征变换

对原始光谱特征进行数学变换将有效消除噪声和突出响应细节, 以便更好地解译光谱信息。 选择用于消除散射水平差异的多元散射校正法(multiplicative scatter correction, MSC)和削弱土壤背景干扰、 基线漂移及分辨重叠峰的一阶导数法(first derivative, FD)进行处理。 相对于单波段参数, 通过两波段线性组合成光谱指数可以放大光谱波段内部相互作用信息、 削弱外部因素的影响, 进而有效提高光谱特征估算理化参数的度量精度。 常见两波段光谱指数组合有差值指数(difference index, DI)、 比值指数(ratio index, RI)和归一化指数(normalized difference index, NDI)三种形式[式(2)-式(4)]。

DI=Ri-Rj(2)

RI=Ri/Rj(3)

NDI=(Ri-Rj)/(Ri-Rj)(4)

式中, R_i、 R_j分别代表在400~1 000 nm范围内第i、 j波长处的光谱反射率。

1.3.2 敏感特征提取和筛选

先利用Pearson相关分析得到光谱原始及变换特征与油茶叶片C/N的响应关系, 再将其中达到极显著水平的波长(p-value≤ 0.01)初筛为敏感波段。 为进一步细化模型输入变量, 使用变量组合集群分析(variable combination population analysis, VCPA)获取光谱特征变量子集。 VCPA^[11]是基于指数递减函数(exponentially decreasing function, EDF)、 二进制矩阵采样(binary matrix sampling, BMS)和模型集群分析(model population analysis, MPA)进行不断迭代筛选的最佳变量子集。 其中, EDF基于达尔文进化理论中“ 适合生存” 的简单有效原则使可变空间连续缩小至小而优化的空间; BMS通过生成二进制矩阵确保以采样几率相等的方式对所有变量随机组合; MPA用于统计BMS产生的大量随机可变组合模型。 为最大程度减少变量组合内部的多重共线性, 对VCPA结果还需要计算两两变量间的相关系数, 设定阈值0.9, 所有超过阈值的各对变量仅保留其中与叶片C/N相关程度最高的, 以得到最优变量组合子集。

1.4.1 数据集划分

首先, 利用留出法按照7∶ 3比例随机将样本数据划分为训练集(Training set)和测试集(Test set); 其次, 利用5折交叉验证对训练集按照平均交叉验证均方误差最小进行调参; 然后, 选取最优超参数组合重新在训练集上训练模型; 最后, 在测试集上评估模型泛化能力。

1.4.2 机器学习模型

选择随机森林、 支持向量机和BP神经网络三种机器学习算法进行油茶叶片C/N估算模型构建。

随机森林(random forest, RF)^[12]是集成算法中Bagging架构的具体实现, 并行构造多个CART决策树, 用多个决策树学习结果的算术平均作为回归结果进行最终的模型输出。 由于每棵树选取特征的随机性, 因此, RF整体模型具有较好的泛化能力和抗过拟合性能。 通常来说, RF训练前需要配置的超参数较多, 其中的决策树个数和深度对模型拟合能力具有较大影响, 所以, 主要对决策树数目M和最大决策分割数N两项参数进行调优。

支持向量机(support vector machine, SVM)^[13]是一种基于Vapnik-Chervonenkis维理论和结构风险最小化原则的监督式学习算法, 能够利用有限的训练样本信息寻求模型复杂度与学习能力之间的最佳平衡。 SVM理论的关键在于核技巧, 即通过构造核函数避免特征从低维空间映射到高维空间的复杂运算。 本工作选取鲁棒性较好的径向基函数(radial basis function, RBF)作为核函数, 并用网格搜索法对惩罚系数C和核参数γ 两个超参数加以调优。

BP神经网络(back propagation neural network, BPNN)^[14]是一种基于Widrow-Hoff学习规则广泛应用的多层前馈神经网络, 能够模拟多种非线性关系。 BPNN的拓扑结构包括输入层、 隐含层和输出层, 通过设定网络运行参数(迭代数、 学习率等)计算各层权重和偏置, 使用梯度下降法通过误差逆反馈对权重和偏置进行不断修正, 最终使网络的误差平方和达到最小。 BPNN通常包括隐藏层个数和各隐藏层节点数这两个重要的超参数, 配置网络时需要指定这两个参数的值, 层数和节点数的最优配置仍然是目前具有挑战的未知问题。 为了获取泛化能力更好、 预测精度更高的神经网络模型, 将隐藏层数设定为3层, 依据式(5)并采用网格搜索法配置各层的最佳节点数。

$n=\sqrt{I+O}+\alpha $(5)

式(5)中, n是隐藏层节点数, I和O分别是输入层和输出层的变量维数, α 是1~10之间的调节常数。

1.4.3 模型评价与比较

选取决定系数R²(coefficient of determination)、 均方根误差RMSE(root mean square error)和相对分析误差RPD(relative percent deviation)作为模型评估指标, 计算公式见式(6)— 式(8)。

${{R}^{2}}=\overset{n}{\mathop{\underset{i=1}{\mathop \sum }\, }}\, {{({{{\hat{y}}}_{i}}-{\bar{y}})}^{2}}/\overset{n}{\mathop{\underset{i=1}{\mathop \sum }\, }}\, {{({{y}_{i}}-{\bar{y}})}^{2}}$ (6)

$\text{RMSE}=\sqrt{\overset{n}{\mathop{\underset{i=1}{\mathop \sum }\, }}\, {{({{{{\hat{y}}}}_{i}}-{{y}_{i}})}^{2}}/n}$ (7)

RPD=SDy/RMSE(8)

式中, y_i、${{{\hat{y}}}_{i}}$、${\bar{y}}$分别表示油茶叶片C/N的观测值、 预测值和平均值, SDy表示观测值的标准差, n为样本容量。 通常来说, R²越大且RMSE越小, 模型效果越好; RPD越大, 精度越好, RPD≥ 1.4时表示模型具有较好的估算能力^[15]。 所有相关操作均在Matlab(Version R2020b, MathWorks.Inc)编译环境中完成。

通过计算每株油茶叶片TC与TN的比值可以得到一个无量纲参数叶片C/N, 三者的统计描述如表1所示。 由表1可知, 叶片TC含量为弱变异, 变异系数仅有1.71%, 说明TC总体稳定离散分布在511.92 g· kg^-1附近, 而叶片TN和叶片C/N均为中等变异且十分接近, 体现出收获期油茶氮元素含量在不同单株间存在较大差异。 叶片TC偏度值为负值, 呈现左偏态分布; 而TN和C/N偏度值为正值, 均为右偏态分布。

表1

Table 1

表1(Table 1)

表1 TC、 TN和叶片C/N的描述统计 Table 1 Descriptive statistics of TC, TN, and leaf C/N 变量 TC/(g· kg-1) TN/(g· kg-1) C/N 平均值 511.92 11.30 45.96 标准差 8.77 1.36 5.49 变异系数/% 1.71 12.08 11.95 偏度 -0.29 0.28 0.33 峰度 0.26 0.18 0.16

表1 TC、 TN和叶片C/N的描述统计 Table 1 Descriptive statistics of TC, TN, and leaf C/N

图2为油茶叶片C/N分别与叶片TC、 TN的响应关系, 由图2可知, 油茶叶片C/N幅度变化与TN的关联程度明显大于TC, 且呈极显著的负线性相关(r=-0.98, p=0.00< 0.01), 而TC累积量的振荡现象揭示, 主要是氮元素在碳氮代谢协调中充当限制因子。 油茶收获期的叶片C/N差异在一定程度上也反映出树体营养器官与生殖器官之间对同化产物的争夺和转化情况。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 叶片C/N与TC、 TN响应关系的双Y轴图Fig.2 Double Y-axis diagram of the response relationship between leaf C/N and TC, TN

两种特征变换引起原始光谱曲线(R)的变化情况(图3)表明, MSC光谱[图3(b)]总体波峰走势基本与R光谱[图3(a)]保持一致, 仅对绿光波段500 nm处至近红外波段900 nm处之间的反射率标准差(standard deviation, S.D.)有明显减弱, 而对其余波段特征改变不明显。 和R与MSC光谱相比较, FD光谱[图3(c)]出现了较多的凹凸区域, 在可见光范围内由峰谷特征构成的光谱曲线起伏均有不同程度的放大, 尤其体现在近红外波段过渡区域的红边增强效应最直观。 三种光谱曲线反射率在近红外波段900nm处均有较大的波动, 这与植物叶片的内部构造差异有关。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 不同特征变换的光谱反射率曲线 (a): R光谱; (b): MSC光谱; (c): FD光谱Fig.3 Spectral reflectance curves of different feature transformations (a): R spectra; (b): MSC spectra; (c): FD spectra

进一步用油茶叶片C/N与R、 MSC和FD光谱反射率分别构造相关系数曲线, 由图4可知, 油茶叶片C/N与R光谱反射率在400~745 nm之间呈正相关, 而在其他波段呈负相关, 最高相关性的波长位置在725 nm(r=0.61), 其中敏感波段(p≤ 0.01)范围主要位于515~660、 690~740和760~940 nm。 MSC光谱反射率明显增加了在400~520 nm范围内与叶片C/N的关联程度, 在该区域维持较高的负线性相关; 同时, 相对于R相关曲线, MSC凸显了在650~680 nm的负相关, 但是, 在750 nm之后光谱特征响应的敏感程度则明显下降, 其最高相关性的位置与R保持一致, 均在725 nm处(r=0.61)。 FD光谱反射率与叶片C/N的敏感特征响应区间主要在可见光范围(400~780 nm), 其中, 在485~555和680~710 nm呈现极显著的正相关, 而在570~675和720~775 nm为极显著的负相关, 且仅在近红外区域的840~870和945~950 nm有极显著的线性相关。 使用VCPA方法对R、 MSC和FD光谱反射率的敏感特征进一步遴选并去除共线性, 得到如图4所示最优特征变量子集位置分布, 其中, R光谱的变量组合集中在红边范围, 而MSC和FD光谱的变量组合子集在多个敏感特征响应区域均有分布。 这些结果表明, MSC和FD特征变换处理可以有效提高VCPA挖掘潜在信息的能力。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 R、 MSC和FD与叶片C/N的相关系数Fig.4 Correlation coefficients of R, MSC, and FD with leaf C/N

按照两波段差值(DI)、 比值(RI)和归一化(NDI)组合方式, 进一步使用R、 MSC和FD光谱反射率计算光谱指数与油茶叶片C/N的相关系数矩阵(图5)。 结果显示, R光谱在三种光谱指数中的敏感区域基本一致, 其中, 波长在550和720 nm附近的反射率分别与可见光、 近红外波段组合呈现一正一负两种响应关系, 而在920 nm附近由波长反射率自身的组合表现出显著的正线性相关, 同时其RI和NDI的处理效果非常相似。 MSC光谱在三种光谱指数中的敏感区域相较于R光谱均有所拓展, 图5(e)和(f)表明, MSC处理可提高可见光与近红外波段组合的响应程度。 由图5(g— i)中可知, FD的光谱指数从DI、 RI到NDI逐步细化敏感特征的响应区域, 其中FD-NDI光谱形成大量的敏感细碎斑块。 总体来看, 两波段组合相对于单波段变换能够提高光谱特征与油茶叶片C/N最大的响应程度, 但效果不明显。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 不同单波段变换的DI、 RI和NDI与叶片C/N的相关系数矩阵Fig.5 The correlation coefficients matrix of DI, RI and NDI with leaf C/N in different single band transformations

使用VCPA方法面向三种单波段变换的DI、 RI和NDI光谱指数, 从对应处理的敏感区域反射率组合值中遴选出变量子集并去除共线性得到如表2所示的最优组合分布情况。 由表2可知, FD-RI与FD-NDI处理能够保留较多可见光与近红外波段组合的响应信息。

表2

Table 2

表2(Table 2)

表2 两波段最优光谱特征变量组合子集 Table 2 Optimal spectral characteristic variable combination subset of two bands R [i, j] MSC [i, j] FD [i, j] DI [725, 695] [755, 735] [920, 440] DI [575, 570] [755, 735] [875, 830] [985, 930] DI [540, 405] [975, 745] RI [580, 455] [615, 520] [630, 530] [755, 740] RI [465, 460] [565, 550] [755, 740] RI [635, 465] [640, 540] [785, 745] [800, 705] [820, 750] [840, 710] [880, 495] [890, 410] NDI [550, 480] [750, 735] NDI [470, 445] [590, 565] [730, 700] [750, 740] NDI [605, 415] [665, 475] [755, 545] [755, 565] [890, 490] [930, 535] [945, 845] [955, 480] [960, 920]

表2 两波段最优光谱特征变量组合子集 Table 2 Optimal spectral characteristic variable combination subset of two bands

用R、 MSC和FD及对应两波段组合指数分别构建三种机器学习模型(RF、 SVM与BPNN)的效果见表3— 表5, 各表中均相应列出经过交叉验证调参的参数配置。 由表3可知, RF模型在训练集中的R²范围为0.33~0.69, RMSE为2.85~4.20; 在测试集中的R²为0.36~0.60, RMSE为3.85~4.91。 结果显示, 经过MSC处理的RF模型表现最佳, 其测试集R²为0.60且RPD达到1.61, 可以说明该RF模型具有较好的叶片C/N估算能力。

表3

Table 3

表3(Table 3)

表3 不同处理下RF模型的评价指标汇总 Table 3 Summary of evaluation indexes of RF models under different treatments 处理 超参数 训练集(n=84) 测试集(n=36) M N R2 RMSE R2 RMSE RPD R 14 1 0.36 4.11 0.41 4.68 1.32 R-DI 18 2 0.40 4.00 0.49 4.38 1.42 R-RI 62 8 0.57 3.37 0.51 4.27 1.45 R-NDI 164 8 0.53 3.55 0.36 4.91 1.26 MSC 6 32 0.62 3.16 0.60 3.85 1.61 MSC-DI 8 1 0.33 4.20 0.40 4.72 1.31 MSC-RI 5 4 0.46 3.77 0.46 4.47 1.39 MSC-NDI 82 1 0.34 4.17 0.39 4.77 1.30 FD 21 2 0.40 3.97 0.49 4.35 1.42 FD-DI 4 8 0.49 3.69 0.55 4.09 1.51 FD-RI 4 8 0.63 3.12 0.47 4.46 1.39 FD-NDI 49 16 0.69 2.85 0.50 4.34 1.43

表3 不同处理下RF模型的评价指标汇总 Table 3 Summary of evaluation indexes of RF models under different treatments

表4

Table 4

表4(Table 4)

表4 不同处理下SVM模型的评价指标汇总 Table 4 Summary of evaluation indexes of SVM model under different treatments 处理 超参数 训练集(n=84) 测试集(n=36) C γ R2 RMSE R2 RMSE RPD R 1.48 0.53 0.32 4.23 0.47 4.47 1.39 R-DI 2.04 0.49 0.34 4.18 0.50 4.31 1.44 R-RI 1.51 0.55 0.32 4.24 0.44 4.55 1.36 R-NDI 5.75 0.92 0.35 4.16 0.50 4.31 1.44 MSC 2.69 0.21 0.34 4.20 0.48 4.40 1.41 MSC-DI 8.71 0.25 0.48 3.72 0.54 4.16 1.49 MSC-RI 5.75 0.78 0.41 3.95 0.45 4.52 1.37 MSC-NDI 5.01 0.56 0.38 4.06 0.44 4.57 1.36 FD 7.59 0.27 0.49 3.67 0.54 4.15 1.49 FD-DI 5.50 1.08 0.35 4.15 0.50 4.33 1.43 FD-RI 1.29 0.18 0.32 4.23 0.29 5.15 1.20 FD-NDI 3.02 0.17 0.50 3.65 0.29 5.16 1.20

表4 不同处理下SVM模型的评价指标汇总 Table 4 Summary of evaluation indexes of SVM model under different treatments

表5

Table 5

表5(Table 5)

表5 不同处理下BPNN模型的评价指标汇总 Table 5 Summary of evaluation indexes of BPNN model with different treatments 处理 节点数 训练集(n=84) 测试集(n=36) R2 RMSE R2 RMSE RPD R 5, 8, 4 0.68 2.90 0.57 4.01 1.54 R-DI 3, 5, 9 0.62 3.17 0.51 4.28 1.45 R-RI 3, 8, 12 0.58 3.32 0.56 4.04 1.53 R-NDI 4, 5, 7 0.51 3.60 0.48 4.40 1.41 MSC 3, 11, 7 0.65 3.05 0.62 3.76 1.65 MSC-DI 3, 8, 11 0.58 3.32 0.47 4.47 1.39 MSC-RI 3, 5, 3 0.51 3.62 0.51 4.27 1.45 MSC-NDI 4, 6, 5 0.52 3.56 0.51 4.30 1.44 FD 3, 10, 11 0.67 2.95 0.52 4.25 1.46 FD-DI 3, 6, 5 0.52 3.56 0.51 4.29 1.44 FD-RI 4, 4, 11 0.60 3.26 0.53 4.20 1.48 FD-NDI 4, 11, 5 0.71 2.76 0.66 3.56 1.74

表5 不同处理下BPNN模型的评价指标汇总 Table 5 Summary of evaluation indexes of BPNN model with different treatments

表4为经过不同处理构建的SVM模型评价效果, SVM模型在训练集中的R²范围为0.32~0.50, RMSE为3.65~4.24; 在测试集中的R²为0.29~0.54, RMSE为4.15~5.16。 结果显示, MSC-DI和FD处理的SVM模型具有较好的估算能力, 其测试集的R²均为0.54, RPD均达到1.49。

表5为光谱经过不同处理构建的BPNN模型评价效果, BPNN模型在训练集中的R²范围为0.51~0.71, RMSE为2.76~3.62; 在测试集中的R²为0.47~0.66, RMSE为3.56~4.47。 其中, 经过FD-NDI处理的BPNN模型在测试集的R²为0.66, RPD达到1.74, 说明该BPNN模型相对来说具有更高的准确性和泛化能力, 可以用此模型进行收获期油茶叶片C/N的估算。

上述结果分析可知, 由不同处理产生的变量组合子集内部差异会对模型估算精度有较大影响, 而三种机器学习模型整体效果最好的是BPNN, 其次是RF和SVM。 从计算机对算法运行的效率来看, 网络结构更加复杂的BPNN需要的运行时间更长。 为了更好地解释模型的估算特征, 选取RPD评价效果最好的三种模型进行展示, 即FD-NDI-BPNN[图6(a)]、 MSC-BPNN[图6(b)]和MSC-RF[图6(c)]。 由图6可知, 这三个模型在测试集中的回归线斜率均低于1∶ 1线, 因此, 模型倾向于过高估计低C/N和过低估计高C/N。 同时, 数据回归分布的95%置信区间(confidence interval, C.I.)特点表明, 即使在同一生长时期C/N水平相近的油茶个体, 其模型估算精度也会有较为明显的差异, 这可能与选取范围内光谱携带信息的上限有关, 需要考虑其他影响油茶生长状态的生理生化指标来综合提升模型估算精度。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 选取模型的叶片C/N观测值与预测值比较Fig.6 Comparison of the observed and predicted values of leaf C/N with different models