土壤样本采样区域为宁夏银川润丰种业公司试验田(宁夏, 青铜峡市, 38° 06'57″N、 106° 08'16″E), 地处西北内陆, 东部季风区与西部干旱区域的交汇地带, 属中温带大陆性气候, 年降水量177.8 mm, 年平均气温9.8 ℃, 土壤类型为黄土, 种植作物为玉米, 采样区域的中间地势相较于上下两侧偏低。 图1为试验田航拍图和采样点分布图。 采样时间为2023年11月, 采用棋盘式采样法采集土壤, 采样深度为距表层约20 cm左右, 单次采集土壤质量为1 kg。 在玉米不同生长期, 3次采样, 将采集的288份土壤样本在室内晾干、 去杂等步骤后均匀分成2份。 一份用于LIBS光谱数据采集; 一份经过4分法分成2份, 分别用于采集可见-近红外光谱数据和重铬酸钾容量法测定SOM含量。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 土壤采集试验田及采样点分布 (a): 试验田; (b): 采样点Fig.1 Soil collection experimental field and distribution map of sampling points (a): Experimental field; (b): Soil sampling points

在采集LIBS光谱数据前, 将土壤样本烘干、 过筛(100目)、 过滤、 研磨, 最后得到土壤基体粉末。 为进一步去除粉末中的水分, 将土壤基体粉末置于烘箱中80 ℃烘干6 h得到粉末土壤样本。 将粉末土壤样本重新研磨, 并在压片机30 MPa 压力下保持5 mins制得重量约3 g、 直径约25 mm、 厚度约3 mm的圆饼状土壤压片。 在采集可见-近红外光谱数据前, 将土壤样本过1 mm孔筛后, 置于直径75 mm、 深度为10 mm的培养皿中装满, 并用钢尺将其表面刮平。

1.2.1 LIBS试验装置及数据采集

LIBS系统装置如图2(a)所示, 包括中阶梯光栅光谱仪ARYELLE 200(Lasertechnik Berlin, 德国)、 增强电荷耦合器件(intensified charge coupled device, ICCD)探测器、 激光器TINY-200L(卓镭激光有限公司, 北京)、 激光器NASOR-900(卓镭激光有限公司, 北京)、 X-Y-Z 3D样品移动平台(天津三英精密仪器股份有限公司, 天津)、 收集镜、 聚焦透镜、 计算机等。 开始测量前, 为确保光谱仪和ICCD探测器的正常工作, 需调整激光器的能量、 聚焦镜头的位置以及设定延时脉冲装置的时间延迟。 数据采集具体流程: 将土壤压片样本放置在LIBS系统的XYZ样品移动平台上, 将激光器产生的脉冲激光束聚焦到土壤压片样本上, 诱导样本表面产生等离子体, 等离子体发射出的光, 通过光谱仪分光后被ICCD探测器记录, 利用LTB Sophi软件处理采集到的光谱数据, 包括谱线的识别、 背景的校正、 谱线强度的测量, 以识别和量化土壤样本中的元素。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 光谱数据采集装置示意图 (a): LIBS光谱系统; (b): 可见-近红外光谱采集系统Fig.2 Schematic diagram of spectral data acquisition devices (a): LIBS device; (b): Visible-near-infrared spectral acquisition device

1.2.2 可见-近红外光谱试验装置及数据采集

可见-近红外光谱试验装置示意图如图2(b)所示。 为保证光谱数据的准确性, 开始测量前, 需对光谱采集系统预热15 mins, 并对仪器进行暗电流校正和白板校正, 在400~1 200 nm采样间隔为1.4 nm, 重采样间隔为1 nm。 删去两端噪声较大的光谱区间, 保留475~1 100 nm波段, 每个样本采集5次光谱数据后取平均值。

1.2.3 特征波长选择和预测模型构建

将采集的全谱用于建模会导致训练速度下降, 并且光谱数据可能出现信息冗余和光谱重叠, 从而影响模型的准确性^[12]。 连续投影算法(SPA)是一种前向循环选择算法, 可以消除变量间共线性; 竞争自适应重加权采样法(CARS)是一种集成蒙特卡洛抽样和偏最小二乘回归模型(partial least square regression, PLSR)回归系数的降维算法, 可以对无信息变量进行有效去除^{[13, 14]}。 CARS和SPA被广泛应用于光谱信息的数据降维。

偏最小二乘回归(PLSR)能够很好地从土壤光谱信息估测土壤的理化信息^[15]。 反向传播人工神经网络(back propagation neural network, BPNN)是一种通过模仿人的神经元激活和传输过程对误差进行反向传播训练, 而输出结果按照前向传播的多层前馈网络^[16]。 本研究将R² (R-squared)、 RMSE (root mean square error) 和RPD (ratio of performance to deviation)作为评价指标。 其中R²越接近1, 模型的预测能力越好; RMSE值越低, 模型的预测准确性越高; 在土壤科学中, RPD值小于1.4的模型被认为不适合预测, RPD值在1.4到2之间的模型被认为具有良好的预测能力, 而RPD值大于2的模型被认为具有优秀的预测能力^[12]。

用重铬酸钾容量法进行化学分析的288份样本的测量结果, 采用SPSS 27进行Shapiro-Wilk测试, 结果如图3所示。 土壤样本的SOM含量的范围为6.79~25.86 g· kg^-1, 均值为15.70 g· kg^-1, 标准差为2.82 g· kg^-1, K-S检验值是0.212(p> 0.05), 表明SOM含量数据符合正态分布。 宁夏青铜峡市的SOM含量平均值为14.9 g· kg^-1, 与之相比, 本研究的土壤样本的SOM含量平均值略高。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 土壤有机质含量的统计结果Fig.3 Statistical description of SOM content in soil

2.2.1 土壤样本的LIBS光谱预处理

为了抑制背景干扰和噪声的影响, air-PLS多次迭代优化过程, 逐步拟合出光谱数据中的信号成分^[17], 在迭代过程中, 算法会根据当前基线估计对数据进行加权, 以增强信号与基线的区分度。 为了避免过拟合, air-PLS引入了惩罚函数, 使得基线估计更加平滑。

图4为土壤样本在300~800 nm的预处理前后的平均激光诱导等离子光谱曲线。 土壤样本的光谱在经过air-PLS基线矫正后的谱线变得平滑, 光谱数据中的背景干扰和噪声被很好地去除。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 LIBS平均光谱预处理前后的曲线Fig.4 LIBS average spectra before and after spectral pretreatment

2.2.2 土壤样本的可见-近红外光谱预处理

为降低仪器和环境噪声、 土壤样本中散射效应的影响, 采用SG卷积平滑(Savitzky Golay convolution smoothing, SGCS)、 一阶导数变换、 多元散射矫正(multivariate scattering correction, MSC)对光谱数据进行预处理。 图5(a)— (d)分别为475~1 100 nm的原始可见-近红外光谱反射曲线、 SG卷积平滑后的光谱曲线、 一阶导数变换后的光谱曲线、 MSC后的光谱曲线。 表1为不同光谱预处理方法的PLSR模型性能评价。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 原始光谱曲线和预处理后的光谱对比 (a): 原始可见-近红外光谱曲线; (b): SG卷积平滑; (c): 一阶导数变换; (d): 多元散射校正Fig.5 The original spectral curve and the spectral curve after different pretreatment methods (a): Raw visible-near-infrared spectral curves; (b): SG convolution smoothing; (c): The first derivative transformation; (d): MSC

表1

Table 1

表1(Table 1)

表1 不同光谱预处理方法的PLSR模型性能评价 Table 1 Performance evaluation of PLSR models with different spectral preprocessing methods 光谱预处 理方法 Rc2 RMSEC/ (g· kg-1) Rp2 RMSEP/ (g· kg-1) RPDv 未预处理 0.675 1.722 0.525 1.449 1.453 SG 0.844 1.015 0.707 0.870 1.768 D 0.972 0.425 0.685 1.031 1.775 MSC 0.897 0.823 0.759 0.910 1.859

表1 不同光谱预处理方法的PLSR模型性能评价 Table 1 Performance evaluation of PLSR models with different spectral preprocessing methods

由表1可知, 3种光谱预处理方法均可以提高模型性能。 经过MSC预处理的光谱数据的PLSR模型效果最好。 一阶导数法可以增强光谱的细节特征, 但不可避免地会加大噪声; MSC通过调整光谱的基线和斜率, 使其更接近于无散射或最小散射的条件, 可以减少由样品的物理属性(如粒度、 包装密度等)引起的散射效应, 本研究将MSC预处理后的光谱数据作为后续可见-近红外光谱的特征波长选择的输入。

2.2.3 土壤样本SOM含量的异常值剔除

为了降低光谱数据采集和SOM含量值测定过程中的误差对建模的影响, 采用蒙特卡洛交叉验证法(Monte carlo cross validation, MCCV)剔除异常样本。 蒙特卡洛迭代次数设置为2000次, 训练集与验证集的比例设置为3∶ 1, 判定异常的阈值设置为1.0 g· kg^-1, 采用SPXY法划分样本集, 并采用PLSR模型对训练集进行训练, 将平均预测误差作为模型性能评估指标, 将平均预测误差大于2 g· kg^-1的样本视为异常样本, 共剔除25个异常样本, 如图6所示。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 288份土壤样本的平均预测误差值Fig.6 The average prediction error of 288 soil samples

将NIST库的元素谱线信息和实际采集的4种金属元素(Na、 K、 Mg、 Mn)的光谱信息相结合, 根据谱线发射强度、 自发跃迁概率、 激发能、 高能级和低能级的简并度, 选取Na Ⅰ 589.59 nm、 K Ⅰ 769.89 nm、 Mg Ⅰ 518.36 nm、 Mn Ⅰ 445.30 nm为元素特征谱线, 图7为4种金属元素的特征谱线图。 分别将Na特征谱线附近589.26~589.92 nm、 K特征谱线附近769.70~769.99 nm、 Mg特征谱线附近518.17~518.62 nm、 Mn特征谱线附近445.10~445.65 nm范围内的谱线信号用于后续建模。

图7

Fig.7

	Figure Option ViewDownloadNew Window
	图7 Na、 K、 Mg、 Mn4种金属元素的特征谱线 Ⅰ 表示元素的原子谱线Fig.7 Characteristic spectral lines of Na, K, Mg, and Mn Ⅰ : atomic spectral lines

为了提高建模效率, 降低模型复杂度, 采用CARS和SPA算法对475~1 100 nm的光谱数据进行特征波长选择。 将经过MCCV、 MSC后的光谱数据和SOM含量作为CARS算法的输入, 在第19次迭代时, 即RMSECV取最小值0.639时, 选择出77个特征波长用于后续预测模型建立。 图8(a)为CARS选择特征波长的分布。

图8

Fig.8

	Figure Option ViewDownloadNew Window
	图8 特征波长的分布 (a): CARS; (b): SPAFig.8 Distribution of characteristic wavelengths (a): CARS; (b): SPA

将经过MCCV、 MSC处理后的光谱数据和SOM含量作为SPA算法的输入, 选择RMSE值最小的变量子集作为特征波长, 最终, 选择了479、 484、 488、 640、 811、 874、 996、 1 005、 1 013和1 045 nm共10个特征波长。 图8(b)为SPA选择特征波长的分布。

由表2可知, 不同特征波长选择算法建立的模型效果: CARS> SPA> 未选择。 该结果与筛选出的特征波长数量相关, SPA选择的波长较少, 一定程度上剔除了原始光谱中部分相关性较多的信息, 使得SPA选择后建立的模型的均方根误差高于原始光谱建立的模型。 与SPA相比, CARS选择的波长较多且分布较为均匀, 保留了原始数据中绝大部分重要信息, 模型达到了更好的效果; 在500~780nm的可见光波段, 光谱的吸收受有机质特性、 Fe²⁺、 Fe³⁺和Mn³⁺的影响, 这些因素同时也对土壤的反射光强度有影响。 CARS算法在可见光波段也提取了一些波长。 将CARS算法选择的特征波长用于后续SOM含量的建模。

表2

Table 2

表2(Table 2)

表2 不同特征波长提取算法的PLSR模型评价 Table 2 Evaluation of PLSR models based on different feature wavelength extraction algorithms 波长选 择算法 Rc2 RMSEC/ (g· kg-1) Rp2 RMSEP/ (g· kg-1) RPDv 未选择 0.897 0.823 0.741 0.910 1.859 CARS 0.920 0.729 0.779 0.776 2.011 SPA 0.830 1.057 0.678 1.973 1.727

表2 不同特征波长提取算法的PLSR模型评价 Table 2 Evaluation of PLSR models based on different feature wavelength extraction algorithms

2.5.1 土壤样本金属元素预测模型的建立和分析

选取已知元素真实含量的11个国家标准土壤样本, 将标准土样和种田土样经过统一的光谱预处理后, 在同一实验条件下采集LIBS光谱, 每一个标准土样测量20次LIBS光谱数据, 根据220组标准土样的4种金属元素特征谱线附近的数据和已测定的元素含量真实值建立PLSR预测模型。 将288个种田土壤样本的四种金属元素光谱数据分别作为利用国家标准土壤样本建立的PLSR模型的输入, 得到288个种田土壤样本的不同金属元素的预测值。 图9为标准土样的不同金属元素的PLSR模型结果。

图9

Fig.9

	Figure Option ViewDownloadNew Window
	图9 4种金属元素的PLSR预测模型结果 (a): Na; (b): K; (c): Mg; (d): MnFig.9 PLSR prediction results for different metal elements (a): Na; (b): K; (c): Mg; (d): Mn

2.5.2 土壤样本SOM预测模型的建立和分析

根据SOM含量真实值和CARS算法选择的77个特征波长分别建立PLSR、 BPNN模型。 将2份种田土壤样本按照3∶ 1划分为训练集和验证集。 对于PLSR模型, 将主成分因子数设置为5。 对于BPNN模型, 在MATLAB 2022b软件建立了3层人工神经网络, 经试验确定输入层、 隐含层、 输出层神经元分别为77、 18、 1, 将迭代次数和训练精度设置为1 000、 0.01, 建立BPNN模型。 由图10可知, PLSR预测模型效果较好。

图10

Fig.10

	Figure Option ViewDownloadNew Window
	图10 模型的预测结果 (a): PLSR模型; (b): BPNN模型Fig.10 Prediction results of different models (a): PLSR model; (b): BPNN model

对于PLSR模型, 其训练集样本点在高值区域内的样本存在不同程度的低估, 而其中值区域内的样本则能够大致地分布在直线y=x的两侧; 对于BPNN模型, 其训练集样本点在低值区域的样本存在不同程度的高估, 而其中高值区域内的样本则能够大致地分布在直线y=x的两侧。 这种问题的产生主要是因为SOM含量真实值的高值较少, 大部分样点SOM含量都位于中值区域, 因此模型在中值区域拟合的效果较低、 高值区域更好。

PLSR效果比BPNN好的原因主要是因为PLSR能够有效利用光谱数据中的关键信息, 并且对数据分布适应性较强, 其通过自动提取与目标变量最相关的潜在变量来平衡不同区间的样本, 从而提升模型的预测精度; BPNN对输入特征较为依赖, 并且对数据分布较为敏感, SOM含量样本分布不均, BPNN可能会过度拟合中值样本, 从而对低、 高值样本的预测效果较差。

将栅格地图图像转化为灰度图, 并通过阈值分割和轮廓检测来识别二维图像中标定的取样点并记录其坐标, 将288个取样点的预测值一一对应于二维地图上的288个取样点坐标, 随后在 MATLAB 2022b利用已知样本点的值和其与未采样点的空间关系, 构建克里金系统并求解, 用得到的权重对数据值进行插值(克里金插值), 绘制出不同金属元素预测值的分布图。 四种金属元素含量预测值的分布图如图11(a)— (d)所示。 同绘制金属元素含量预测值的分布图的方法, 绘制出SOM含量真实值、 PLSR预测值、 BPNN预测值的分布图, 如图11(e)— (g)所示。

图11

Fig.11

Figure Option ViewDownloadNew Window

图11 4种金属元素的PLSR模型预测值、 SOM含量真实值和PLSR、 BPNN模型预测值的空间可视化分布 (a): Na; (b): K; (c): Mg; (d): Mn; (e): SOM实测值; (f): SOM的PLSR模型预测值; (g): SOM的BPNN模型预测值; (h): 采样点分布Fig.11 Spatial visualization distribution of the predicted values of PLSR model, true value of SOM content, and predicted values of PLSR and BPNN models of four metal elements (a): Na; (b): K; (c): Mg; (d): Mn; (e): SOM true value; (f): Predicted SOM value by PLSR model; (g): Predicted SOM value by BPNN model; (h): Distribution of sample points

由图11(e)— (g)可知, 采样区域的SOM含量集中在14~18 g· kg^-1, 采样区域上下两侧SOM含量总体偏低, 采样区域的左上方土壤颜色较深, SOM含量最高。 本研究的采样区域中间的地势相较于上下两侧较低, 而斜坡上的水流的侵蚀作用会导致SOM流失, SOM在斜坡上的积累通常少于平坦地区; 水分是影响SOM分解和植物生长的关键因素, 水分充足的地区通常有较高的植被生产力, 一定程度上可以增加SOM含量, 所以采样区域的中间位置可能因为水分较多而有较高的SOM含量。