野外试验在以色列内盖夫沙漠北部的Migda实验农场和加利利地区东部的Karei-Deshe实验农场开展, 土壤类型分别为沙质黄壤土和玄武岩基棕壤土。 不同的土地利用方式, 如小麦种植、 放牧、 休耕等使土壤的AN含量有较大差异。 分别在两个农场随机选取了48和28个采样点, 共76个。 野外光谱试验时间选在干旱季节的末期, 能够最大程度减小环境因子如水分、 植被和云层对试验的影响^[13]。

采用美国ASD公司的Field Spec^® 便携式光谱仪进行野外原位土壤光谱的测量。 为了减小自然光的影响, 利用自带光源的高密度接触式反射探头直接垂直紧贴土壤表面进行光谱测量。 每个采样点采集4次用以得到平均光谱, 同时每次测量均利用标准白板校正, 得到采样点的反射率。 该光谱仪波长范围为350~2 500 nm, 光谱分辨率在350~1 000 nm为3 nm, 在1 001~2 500 nm波段为8~10 nm, 经重采样后, 得到分辨率为1 nm的反射率光谱。

光谱采集后, 利用土钻在每个采样点0~15 cm深处进行土壤采样, 并封存送到实验室进行风干、 研磨和过筛。 利用称重法测量土壤水分, 所有的样本含水率均不超过0.2%, 因此可以忽略水分对光谱的影响。 最后使用氯化钾溶液提取法对土壤AN的含量进行测量, 测量结果统计数据如表1所示。

表1

Table 1

表1(Table 1)

表1 土壤速效氮的描述统计量 Table 1 Descriptive statistics of soil AN (mg· 100 gr-1) Category Min Max Mean SD No. Migda 10.92 35.12 20.33 6.79 48 Karei-Deshe 23.25 101.90 43.27 16.99 28 Total 10.92 101.90 28.78 16.02 76

表1 土壤速效氮的描述统计量 Table 1 Descriptive statistics of soil AN (mg· 100 gr-1)

在两组数据集中各随机选取70%的样本(共54个)作为训练集结合建模, 以提高预测模型的适用性; 剩余的22个样本作为测试集, 验证模型的泛化效果。

光谱反射率曲线在首尾区域有较低的信噪比, 影响建模预测精度, 因此仅使用400~2 400 nm的光谱作为原始光谱 (raw spectra, RS)进行进一步分析。 使用以下组合方法对光谱进行预处理与变换: (1)首先使用2阶15点多项式Savitzky-Golay卷积平滑算法(SG)对光谱曲线进行平滑去噪; (2)再对光谱进行一阶导数(FD)变换去除基线和背景的影响; (3)然后使用标准正态变换(SNV)来消除固体颗粒大小、 表面散射以及光程变化对光谱的影响; (4)建立与应用滤波模型或回归模型之前, 分别对光谱矩阵和AN浓度向量做均值中心化处理, 以保持相同标度。 预处理变换后的光谱称为PPT光谱。

浓度相似或相同的样品, 在外界因子的影响下, 光谱会出现一定的差异。 GLSW以差异矩阵为参考, 建立加权滤波模型, 用以消除外界因子对光谱的影响^[10]。 但当样品的浓度差异较大时, 差异光谱矩阵会携带与被测对象(即AN含量)相关的特征信息, 而不应该被滤掉。 为此, Zorzetti等^[14]提出了基于Y-梯度变化的广义最小二乘加权算法(Y-gradient general least square weighting, Y-GLSW)。 算法主要思路如下:

(1)按照被测对象浓度y梯度递增的顺序同时对y和光谱矩阵X进行排序, 数据每1行代表1个样本;

(2)利用1阶5点多项式的SG滤波1阶求导, 对X矩阵的每1列(即每个光谱特征)进行导数变换, 得到光谱差异矩阵X_d, 该矩阵每1行(即每个样本)代表了它与邻近4个样本的平均差异, 同时y也执行类似操作, 得到被测对象差异向量y_d;

(3)将y_d转化为对角再加权矩阵W, 其第i个对角元素计算如下: log₂(w_i)=-y_{d, i}s_yd, 其中y_{d, i}是y_d中的第i个元素, s_yd是y_d的标准差;

(4)利用W对X_d再加权, 并计算协方差矩阵: C= XdTW²X_d;

(5)协方差矩阵奇异值分解: C=VS²V^T, 其中V是左特征向量, S是奇异值的对角矩阵;

(6)计算加权和变换脊后的对角矩阵: S_w= S2α+I; 其中I是单位对角矩阵, α 为正则化参数;

(7)利用特征值和加权奇异值的逆来计算滤波矩阵: G=V Sw-1V^T;

(8)利用滤波矩阵对需要校正的光谱即输入光谱X_input进行投影滤波: X_output=X_inputG。

正则化参数α 为标量, 可以用来控制奇异值的增益大小, 进而调节滤波程度。 α 越小(如0.000 1等), S_w趋近于0, 滤波程度越大。 反之, α 越大, S_w趋近于1, 滤波程度越小。

偏最小二乘回归(PLS-R)算法利用对自变量和因变量系统中的数据信息进行分解和筛选的方式, 提取对因变量的解释性最强的综合变量, 辨识系统中的有用信息和无关信息, 克服变量多重相关性的影响, 从而建立适当的预测模型。 建模过程中, 潜在变量(又称主成分)个数是模型优化的关键, 本工作采用10-折交叉验证的方式, 以训练集最小平均误差为标准确定最佳潜在变量个数。

模型预测性能使用相对分析误差RPD(标准偏差除以均方根误差, SD/RMSE), 作为评价标准; 同时利用解释总方差占实际总方差之比SSR/SST(计算方式为: ∑i=1m(f(xi)-y̅)2/∑i=1m(yi-y̅)2)来评价模型的解释能力^[15]。 较好的模型应有较高的RPD和SSR/SST。 根据文献[1]及前人的研究, RPD大于1.4并且SSR/SST大于0.5, 模型才有值得信服的预测能力。

使用MATLAB R2016a (MathWorks, Natick, Massachusetts, USA)和PLS-Toolbox 8.0 (Eigenvector Research Incorporated)进行数据处理分析。

两个实验农场的原始土壤光谱如图1所示, 可以看出两组光谱均有随着波长的增加反射率逐渐增加的趋势, 但Migda农场的土壤反射率明显高于Karei农场。 Migda农场的沙质黄壤土成分主要包括20%的黏土, 38%的壤土和42%的沙土, 而Karei农场的玄武岩基棕壤土却是65%的黏土, 11%的壤土和24%的沙土。 不同的土壤成分使光谱反射率特征产生了较大差异, 且黏土矿物成分及其含量的不同, 使1 400, 1 900和2 200 nm的水分吸收峰差生了形状和深度上的差异。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 土壤样本的反射率光谱Fig.1 Spectral reflectance of soil samples

使用欧氏距离法, 在训练集中选择20个代表性样品进行梯度排序, 按照1.3节中的方法建立滤波模型。 正则化参数α 直接控制滤波程度, 对建模效果有较大影响。 为了选择最优的滤波参数, 先依次从[10^-4, 10^-3, …, 10³]中取值分别建立滤波模型, 对光谱矩阵进行滤波校正, 再使用PLS-R建立回归模型, 测试集的预测结果如图2所示。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 不同α 取值下的模型性能Fig.2 Model performance with various α

从图2可以看出, 随着α 的增加, 即滤波程度的从强到弱, 回归模型的预测精度呈现先增高再降低的趋势, 模型解释能力总体也有类似的结果。 这是因为, 当滤波程度太强的时候也会将一些与土壤AN含量相关的重要特征信息过滤掉, 从而严重影响模型的性能。 当滤波程度太小时, 一些噪声信号依然存在, 对模型产生一定的干扰。 因此选择α =0.01建立滤波模型, 对光谱进行校正处理。

预处理变换后的光谱和利用Y-GLSW滤波校正后的光谱如图3所示。 可以看出, Y-GLSW可以有效的对光谱数据进行滤波, 原始光谱信号中的数据被一定程度上加权收缩^[10], 一些特征峰被移除, 并生成了新的特征峰。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 预处理变换后的光谱和滤波校正后的光谱Fig.3 Preprocessing and transformed spectra and filtering corrected spectra

基于PLS-R算法分别使用原始光谱RW、 预处理变换后的光谱PPT和滤波校正后的光谱Y-GLSW对土壤AN含量建立回归模型, 训练集和测试集的结果如表2所示。 从训练集结果可以看出, 经预处理和变换后, 模型的RPD值和SSR/SST值, 都得到了显著提高, 再经过Y-GLSW滤波校正处理, 两个数值都进一步增加。 在测试集中, RW模型RPD值为1.64, 甚至超出了训练集的结果, 但SSR/SST仅为0.46, 因此预测模型是不可靠的。 PPT模型的测试集RPD值降低较多, 仅为1.41, SSR/SST值却达到了0.57, 因此PPT模型具有一定的可靠性。 利用Y-GLSW滤波变换后的光谱数据进行建模预测, RPD值达到了2.07, 按照文献[1]中的标准, Y-GLSW模型预测能力等级达到了A类。 且模型解释能力也相对于预处理变换光谱模型也得到了明显提高。 因此, 可以确定图3中Y-GLSW滤波去除的信息数据为环境中或系统中的无关噪声信号数据。

表2

Table 2

表2(Table 2)

表2 不同建模方法的模型性能 Table 2 Model performance of different methods Methods LVs Training set Test set RPD SSR/SST RPD SSR/SST RS 6 1.49 0.54 1.64 0.46 PPT 7 4.86 0.96 1.41 0.57 Y-GLSW 11 11.05 0.99 2.07 0.69

表2 不同建模方法的模型性能 Table 2 Model performance of different methods

图4为22个测试集样本的土壤AN实际理化测量值与三种光谱模型预测值的对比散点图。 从图4可以看出, 三种模型在AN含量较低的时候预测结果差异不大, 但当AN含量增加时, RS模型和PPT模型的预测结果就明显变差。 相对来说, Y-GLSW模型的预测数据较为均匀的分布在y=x这条直线的两侧, 说明Y-GLSW模型具有更稳健的预测能力。 这与表2中的结果也是相符合的。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 不同建模方法实测值与预测值的散点图Fig.4 Scatter plot of measured values versus predicted value of different methods

表2中Y-GLSW模型使用了11个潜在变量, 这可能是由于AN在土壤中以无机态氮为主, 而氮素在VIS-NIR光谱的直接响应波段主要是有机态的N— H键发生能级跃迁的倍频峰和合频峰; 且AN的含量较低, 一般仅占总氮的5%^[4], 因此确实需要更多的间接光谱数据信息来反应AN的光谱特征, 这也就解释了模型解释能力不够强, 即SSR/SST值不是很高的原因。 在PPT模型中, 由于噪声信号的存在对模型的训练产生了影响, 虽然仅使用7个潜在变量, 训练集就达到了较好的效果, 但模型在测试集的泛化能力却相对Y-GLSW模型降低了很多。