基于X射线荧光光谱与多特征串联策略的土壤重金属含量预测

任 顺，张 雄，任 东，杨信廷,3，张 力

(1.三峡大学 计算机与信息学院，湖北 宜昌 443002；2.三峡大学 湖北省农田环境监测工程技术研究中心，湖北 宜昌 443002；3.农产品质量安全追溯技术及应用国家工程实验室，北京 100097)

近年来，有毒重金属对我国土壤环境质量危害日益严重，耕地、农产品等污染问题亟待解决,因此，重金属监测已成为农产品保护和安全生产的重要任务[1]。目前，土壤重金属污染检测方法有传统实验室检测法与速测法，其中传统实验室检测方法主要有原子吸收与原子荧光光谱法、电感耦合等离子质谱法等，该类方法虽然准确度和精确度高，但仪器运行条件要求高且不适用于现场快速判断；速测法有X射线荧光光谱法(XRF)、生物传感器技术和激光诱导击穿光谱法[2-3]等。其中，XRF具有检测速度快、成本低且可同时检测多种元素等优点，已被广泛用于土壤重金属污染和农产品检测。相比实验室测定方法，XRF法前处理简单，对主要的重金属可有效监测和快速筛查，可为农田和农作物重金属污染的监测和防控、制定合理的农业发展规划提供科学依据[4]。XRF法快速检测重金属时，针对不同类型和性质的土壤，需建立重金属浓度预测模型，再根据实际检测结果对预测模型优化。由于土壤成分复杂，采用XRF仪获取的光谱数据具有很高的空间复杂程度。此外，所得到的光谱数据往往由数以千计的波长点组成，而样本数量常受制于实验条件和成本，在进行化学建模算法时，常会出现统计学称为维数灾祸的多项式复杂程度的非确定性(NP)难题。模型集群分析策略强调，要最大限度分析已有样本集信息，通过随机采样利用大量子集模型信息，获得数据集的内在结构，而不仅仅依靠单一模型信息，可更好地避免模型对样本的依赖[5-7]。

本研究以土壤重金属的X射线荧光光谱为研究对象，将0～26 keV范围内的光谱数据与浓度梯度法配制的铜(Cu)、锌(Zn)、砷(As)、铅(Pb)、铬(Cr)5种重金属污染土壤样本含量值进行关联，将多特征串联策略区间组合优化-竞争适应性重加权采样-连续投影(ICO-CARS-SPA)算法所提取出的特征变量作为输入，建立偏最小二乘(PLS)模型预测土壤样本中5种重金属含量，并与一些常见变量选择方法及其组合对比，结合光谱特征变量的优选和模型的预测性能，实现了对农田土壤重金属含量高效、精准的定量检测，可为土壤重金属治理和污染防控以及制定科学的农业发展规划提供依据。

图1 ICO算法流程示意图Fig.1 Flowchart of ICO algorithm

1 原理与算法

1.1 试验区间组合优化算法

区间组合优化(ICO)是基于模型集群分析策略(MPA)框架下的化学建模算法[8]，MPA强调从数据集中统计分析时，应通过随机采样，分析大量随机产生的子集模型信息，从不同角度考察数据集的内在性质，获得数据集的内在结构。ICO算法流程如图1所示。

1.2 竞争性自适应重加权采样算法

竞争性自适应重加权采样(CARS)算法选择波数的方法是基于PLS模型回归系数，衡量波长选择的重要依据是PLS回归系数绝对值的大小，通过挑选出交叉验证均方根误差(RMSECV)值最小的子集来对应获得与目标值最优的波长组合。计算则采用指数衰减函数，逐次保留绝对值大的回归系数对应的波数点[9-10]。

1.3 连续投影算法

连续投影算法(SPA)是一种前向变量循环特征波长选取方法，该算法利用向量的投影分析，能找到含有最低限度冗余信息的变量组。SPA一般能够较为有效地消除波长之间共线性，同时提高模型精确度[11-12]。

1.4 多特征串联算法

X射线荧光光谱的高维度无法规避地含有较多噪声和冗余。ICO采用波长区间进行变量优选，能排除大量无用信息和噪声波长点，且在每次优化后的联合区间上进行局部搜索策略，更好地体现了柔性收缩策略优势，减少偶然误差的发生，最终可得到一组位置、组合、宽度等都经过优化的有效波长区间。但ICO算法选择的波长数量一般较多，区间内部不可避免地存在共线性和冗余问题，最终选中的变量区间可采用单一变量选择算法进行精简[13]。

CARS算法作为一种高效的波长选择算法，通过选择少量波长变量子集，可以给出更佳或相当的预测效果。由于该算法在迭代过程中同时引入了蒙特卡罗采样(MCS)和自适应加权采样(ARS)两个随机因素，其单独使用的稳定性难以令人满意[14-15]。

SPA算法生成的每一个波长组合中，新入选的波长点与上一个波长点的相关性均为最低，这也体现了SPA算法能够有效消除波长之间的共线性[16]。但对于X射线荧光光谱而言，有效变量间的投影距离并不一定最大，通过SPA算法挑选出变量子集很可能存在部分无关信息变量甚至是干扰变量。

将ICO选择的优化组合区间结合与竞争适应性重加权采样法(CARS)和连续投影算法(SPA)串联(ICO-CARS-SPA)进行特征波长提取，既能有效地锁定特征波段区间，同时优选出的最优波长子集又能够充分减少最优区间内的共线性和冗余，形成优势互补。

2 实验部分

2.1 材 料

在周边1 km范围内无污染源的农田采集制样土壤，为防止实验器具影响和干扰样本重金属的浓度，样本制作过程采用陶瓷用具。将土壤烘干研磨过0.45 mm孔筛后，参考《土壤环境质量标准》[17]中对Ⅰ、Ⅱ、Ⅲ类土壤中各元素含量的规定并结合农田土壤实际情况制作实验样本。每个样本制作30 g土壤，共91个样本。经初步测量，所采集农田土壤中As、Cu、Cr、Pb、Zn的含量依次为3.0、15、30、18、39 mg/kg。综合国家土壤等级划分标准和所采集土壤的具体情况，设I为各个元素初始浓度梯度，其中1～40号样本的浓度梯度为I；41～50号样本的浓度梯度为2*I；51～60号样本浓度梯度为3*I；61～70号样本浓度梯度为4*I；71～91号样本浓度梯度为J*I；其中J=1,2,3,4,……,n，具体实施方案参照表1。按土壤需求浓度，将重金属标准溶液用丙酮稀释后混入土壤中，混合均匀后置于通风橱中自然挥发，待土壤完全风干后，再研磨混匀，按指定标号放入样品盒中，密封保存。

表1 土壤样本配制方案Table 1 Soil sample allocation scheme

2.2 光谱数据采集

每次取一个样本置于便携式X射线土壤重金属检测仪(由三峡大学协同创新中心与北京农业质量标准与检测技术研究中心联合研制)上，获取样本在0～26 keV范围内共4 096个通道内的光谱信息，每个样本经过转动方向测量3次并求取其平均光谱。由于36和37号样本的光谱数据遗失，同时根据配制的样品浓度绘制91个样本的实际浓度曲线，发现60号样品中的Zn浓度和79号样品中Pb浓度分布与配制浓度差异较大，分析原因可能是溶液配制时出现了问题，故剔除这4个样本。剩余87个样本的平均光谱图见图2，其中配置Ag靶微型X光管、电子冷却Si-PIN探测器，设置电压为30 kV，电流为40 μA，积分时间为300 ms。

图2 87个样本的平均光谱图Fig.2 Average spectrogram of 87 samples

2.3 建模方法与评价

2.3.1 光谱数据集的划分该方法以样本集被测指标的理化参照值作为划分标准，将剔除4个异常样本后的87个不同浓度重金属的土壤样本通过浓度梯度法进行校正集和验证集的划分，将样本按照2∶1分成2组，得校正集58个，验证集29个，样本分布如图3。

图3 校正集(A)和验证集(B)样本分布图Fig.3 Distribution figure of calibration set(A) and prediction set(B)

3 结果与讨论

3.1 ICO算法选择特征波长子区间

采用ICO算法进行特征变量选择，并建立PLS模型预测重金属含量。不同区间划分下试验结果如表2所示，当子区间划分数为20时，5种重金属均获得最小的均方根误差值。

表2 不同子区间数PLS建模对比Table 2 Comparison of PLS modeling for different subinterval numbers

(续表2)

以生物毒性显著的重金属元素Cr为特征提取对象的结果如下：采用ICO算法，进行最优联合区间的选取，利用5折交叉验证法建立PLS模型选择特征变量。选择最佳波长区间划分数量20，PLS模型中最大潜变量数为10，所选子模型的比例为0.05，加权自举采样次数为1 000。在ICO算法迭代过程中，随着迭代次数的增加，每个波长区间的采样权重变化情况见图4A，采样权重值越接近1时，颜色越接近深红色；采样权重值越接近0时，颜色越接近深蓝色；颜色介于深蓝色和深红色之间表明采样权重值处于0和1之间。图中还可见，第6个波长区间的采样权重值在第2次迭代过程中已经约为1，其采样权重值在第3次迭代过程中依然有机会变得小于1，最终该区间由于权重过小被剔除；第9个波长区间在第1～5次迭代过程中颜色偏绿蓝，权重系数约为0.3～0.5，但在第6次迭代过程中权重系数仍有机会变大，最终被选中。由ICO每次迭代中提取的子模型的RMSECV图可见，经过11次迭代后均方根误差趋于稳定(图4B)，此时其值为22.624 5，最终挑选的波长区间为[4，9，13，20](图4C)。在ICO最终选中的联合区间内引入了局部搜索策略，进行宽度的自动优化，最终挑选了805个特征波长。

图4 重金属Cr进行ICO变量选择(区间数20)Fig.4 ICO variable selection for heavy metal Cr(interval number 20)A.sampling weight change graph of each interval during the iteration process(迭代过程中各区间的采样权重变化图),B.RMSECV value of the sub-model extracted in each iteration of the ICO(ICO 每次迭代中提取的子模型的RMSECV值),C.ICO algorithm for wavelength range selection(ICO算法进行波长区间选择)

图5 CARS运算提取变量原理图Fig.5 Schematic diagram of CARS operation extraction variablesA.variation trend of the wavelength variable(波长变量的变化趋势),B.variation trend of the RMSEcv(RMSEcv的变化趋势),C.trend of wavelength regression coefficient(波长回归系数的变化趋势)

3.2 ICO-CARS算法选择特征波长

将ICO算法挑选出的805个特征波长进一步使用CARS剔除区间波长中的无关变量。由于CARS方法在迭代过程引入蒙特卡罗采样(MCS)和自适应加权采样(ARS)两个随机因素，造成每次运行挑选的波长结果不尽相同，为增加波长选择的稳定性和可靠性，对CARS进行100次重复计算，得到其RMSECV平均值为8.642 0，标准偏差(STD)为0.686 3。图5为基于CARS算法的波长变量筛选过程图，其MCS抽样运行次数N为50，PLS主成分数为10，采用5折交叉验证。

由CARS算法筛选有效波长变量的变化趋势可看到变量数呈指数函数下降(图5A)，选择变量的个数也从急剧减少到缓慢递减，最终趋于稳定。波长选择过程采用5折交叉验证得到的RMSECV变化趋势，若RMSECV值减小，表明剔除了无关信息变量；若RMSECV值增大则表明剔除了有效信息变量。由图5B可见，在1～20次的采样过程中，RMSECV呈递减趋势，第20次时的RMSECV值最低，此后开始增加，表明采样20次时有效剔除了ICO提取光谱区间中的无关信息。图5C表示各波长回归系数随着采样次数增加而变化，图中“*”为RMSECV值最低点。由图可见，当运行20次时，RMSECV值最低，此时保留的变量为79个。

3.3 ICO-CARS-SPA算法选择特征波长

经过ICO-CARS变量筛选后，再通过SPA算法进行波长精简(图6)。通过最小误差均方根值(RMSE)来确定最终所选特征波长个数，RMSE值越小，表明模型稳定性越好、精度越高。

图6 ICO-CARS-SPA变量选择Fig.6 ICO-CARS-SPA variable selectionA.RMSE of SPA model varies the number of variables(SPA模型RMSE随变量个数的变化),B.optimal characteristic wavelength selected by SPA(SPA模型选择的最优特征波长)

结果显示，RMSE随着波长个数的增加呈逐渐减小的趋势，当波长个数大于33时，RMSE值变化不再显著，此时RMSE值为15.439 9(图6A)。说明最终优选的33个特征波长保留了更多有效信息，作为最优的波长点个数，选取的波长点在实验光谱中的索引如图6B所示，所选波段占原始光谱信息的0.81%。

采用ICO-CARS-SPA最终优选出的33个敏感波长为输入变量，利用PLS对重金属元素Cr含量建模，并按照相同特征提取步骤对Cu、Zn、As、Pb进行特征优选和PLS建模。结果表明，采用多特征串联策略ICO-CARS-SPA-PLS建模后5种重金属的R2、RMSE和MRB均达到了满意的效果(表3)。

3.4 不同特征波长选择方法比较

本文建立的ICO-CARS-SPA方法与ICO、最小绝对值收敛和选择算法(LASSO)、CARS、ICO-CARS和ICO-SPA等[18-19]其他波长选择算法的对比结果见表3。其中LASSO是一种系数压缩变量选择方法，即通过引入惩罚项使影响较小或者无影响的自变量系数趋近于零，进而可实现只保留与响应变量最相关的解释变量，将大量无关变量剔除，但高维数据中常会面临统计学称为维数灾祸的NP难题，因此不可避免地造成LASSO特征波长选择个数将小于样本个数，无法保证数据集中所有重要信息被完全保留。由于土壤中重金属的X射线荧光光谱数据样本数据集的光谱波长数远大于样本个数，LASSO算法不适合直接应用于波长选择问题。

表3 不同波长选择算法PLS建模性能比较Table 3 Comparison of PLS modeling performance of different wavelength selection

由表3可知，尽管ICO算法可将局部搜索策略应用于最终挑选出的最优联合区间进行宽度的自动优化，但区间挑选算法选择的波长数量一般较多，不利于快速预测，且忽视了联合区间内部光谱间存在的相关性和冗余问题，因此可对优选出的联合区间进行进一步波长筛选。而ICO-CARS建模效果优于ICO-SPA是因为对于X射线荧光光谱而言，有效变量间的投影距离并不一定最大，在初选算法波长数量较多的情况下，SPA筛选出的变量子集中可能包含一些无关信息甚至是干扰变量。

通过对ICO、LASSO、CARS、ICO-CARS、ICO-SPA和ICO-CARS-SPA波长优选对比分析可知，采用波长优选算法建模后5种重金属的R2、RMSE和MRB相对于单一PLS建模大多有不同程度地提升，其中以ICO-CARS-SPA-PLS提升最明显。表明ICO-CARS-SPA-PLS的多特征串联提升策略从整体上比其他模型优化方法建模效果更好，更能兼顾到其他元素，能够增强上一轮变量选择的过程与下一轮变量选择的关联性，根据每轮变量选择的好坏进行动态调整，在一定程度上解决了变量区间选择的“筑巢效应”。

3.5 模型广适性比较研究

为验证ICO-CARS-SPA变量筛选模型具有较好的广适性，加入100次随机分组进行实验，5种重金属元素的评价参数及标准偏差结果见表4。

表4 5种重金属100次随机分组ICO-CARS-SPA-PLS定量模型预测结果Table 4 Prediction results of ICO-CARS-SPA-PLS quantitative model of 100 random groups of five heavy metals

由表3～4可见，由于Zn元素在土壤本底中的浓度较大，导致在制备Zn样本时，其配制浓度跨度较大，建模结果显示验证集上Zn的决定系数提高，但RMSEp效果较差，这可能是由于出现了过拟合现象。因此在后续实验中通过扩大样本数量，对样本中的重金属浓度划分更细致，对Zn元素做进一步改进，可提高其预测效果。由表4可知，ICO-CARS-SPA-PLS模型具有较好的广适性，能够实现对样本中5种重金属的定量检测。

4 结 论

通过多特征串联策略的ICO-CARS-SPA算法，建立了X射线荧光光谱重金属含量的定量检测模型。实验制备了91个土壤样本，剔除4个问题样本后，选取了87个样本，并对比分析了几种常见的光谱特征变量选择方法及其组合。结果表明，以ICO-CARS-SPA算法提取的特征波长建立的PLS重金属含量预测模型精度高、误差小、模型推广性能好。基于X射线荧光光谱结合多特征串联策略算法ICO-CARS-SPA是有效的特征提取方法，降低了模型复杂度，有利于推动X射线荧光光谱技术在农业信息领域上的进一步应用。