优化CARS结合PSO-SVM算法农田土壤重金属砷含量高光谱反演分析

袁自然， 魏立飞， 张杨熙， 余 铭， 闫芯茹

湖北大学资源与环境学院， 湖北 武汉 430002

引 言

重金属污染已经严重危害到人类的生活， 其难降解、 易积累、 毒性大等特征对农作物生长、 产量和品质都有重大影响。 农田土壤重金属污染会抑制土壤环境的优势细菌群落， 抑制微生物的生长， 尤其是它可以通过农作物作为媒介通过食物链， 最终威胁人类身体健康[1]。 传统的土壤重金属检测方法有光度法、 化学分析法、 原子荧光光谱法、 电感耦合等离子体发射光谱法、 表面增强拉曼光谱法等[2]， 这类方法虽精度高， 但是仪器设备成本过大， 耗费大量的人力、 物力， 且不具有推广性。 高光谱分析利用连续的高分辨率光谱波段预测土壤中的重金属含量[3]。 它可以实现大面积的快速测定， 避免了复杂的采样步骤， 而且结合采样对比， 高光谱对重金属的测定结果也有很高的可信度。

国内外学者对土壤重金属高光谱反演已经有了很多研究， Gholizadeh等[4]提出了一种基于支持向量机回归(SVMR)交叉验证的多变量校准方法， 建立了可见-近红外(Vis-NIR)区域反射光谱与土壤中锰(Mn)、 铜(Cu)、 镉(Cd)、 锌(Zn)和铅(Pb)之间的关系。 Angelopoulou等[5]通过偏最小二乘法建模发现不同化学监测方法对于土壤重金属Cd， Cr， Cu和Pb反演结果影响不同， 表明土壤反射率完全可以预测土壤中重金属含量。

谭琨等[6]选择中国北方金属矿区和煤矿区的两个研究区域， 利用竞争自适应重加权算法结合偏最小二乘法(CARS-PLSR)选取特征波段， 结果表明与其他线性模型相比， 非线性模型CARS-PLS-SVM通过适当的光谱特征提取， 在实验中表现出较高的精度。 但是目前反演精度受到原始光谱相关性低的限制， 导致反演精度普遍不高， 且现有特征选取方式和反演方法都有待提高。

本文以洪湖市燕窝镇土壤为研究对象， 首先利用CARS智能选取特征波段， 再利用一阶导数、 高斯滤波、 归一化进行特征提高， 获取皮尔逊相关系数大于0.6的波段作为特征波段， 采用3种回归方式PSO-SVM， SVMR和RFR作为对比， 提高建模精度， 找出该区域反演土壤重金属As的最优方法。

1 研究区概况

洪湖市(113°07′—114°05′E， 29°39′—30°12′N)隶属于湖北省荆州市， 该区域土壤具有弱富铝化、 粘化、 酸性的特点， 土壤性质兼有黄壤和棕壤的某些特征。 燕窝镇为洪湖市辖镇， 为传统粮棉油种植大镇， 位于洪湖东北部， 优良的气候条件和充足的水资源， 使得该区域的农业和养殖业得到了快速发展， 但是大量的农田改造和过量施用农药化肥等致使该地区污染逐渐严重， 洪湖水质条件遭到严重破坏， 很多耕地受到As， Cd， Pb， Cd等重金属污染的威胁， 生态环境受到严峻挑战。

图1 研究区区位图

2 研究方法

2.1 竞争性自适应加权算法

竞争性自适应重加权算法(competitive adaptive reweighted sampling)[7]是种用于光谱特征筛选的新型算法， 将每个光谱波段的集合作为单独个体， 采用自适应重加权采样(adaptive reweighted sampling， ARS)， 利用线性模型偏最小二乘法(PLS)作为适应度函数， 和交叉验证不断优化计算， 最终选择出最优子集， 即回归模型精度最高子集， 淘汰误差较大的变量， 经过多次循环采样， 选择出特征波段。

2.2 基于PSO的SVM参数优化

支持向量机(support vector machine， SVM)是种常用非线性负荷预测模型[8]， 粒子群优化(panicle swarm optimization， PSO)算法是一种基于迭代寻优的群计算技术[9]。 使用PSO算法对SVM的参数进行参数寻优， 利用粒子群算法特征， 不断更新粒子适应度， 直到找到全局最优解， 本文采用PSO优化SVM的参数[10]， 其步骤流程如图2所示。

(1)初始化PSO和SVM参数： 设定群体规模n=30、 最大迭代次数N=300、 学习因子C1=1.8和C2=1.2， 粒子的初始位置和速度等。

(2)分别利用每个粒子对所对应的实验样本的光谱和重金属含量进行预测， 得到各粒子当前均方根误差(RMSE)， 将目前最优适应度与粒子本身最优适应度对比， 如果目前最优适应度优于粒子本身最优适应度， 则更新粒子速度和位置， 从而更新个体极值。

(3)检查是否满足极值条件， 即最大迭代次数和精度要求， 如果满足则将最优核函数σ和正则化参数γ， 代入SVM计算， 否则执行步骤2， 继续搜索。

图2 PSO优化SVM参数流程图

2.3 模型验证

选取了决定系数(determination coefficients，R2)、 均方根误差(root mean squared error， RMSE)、 平均绝对误差(mean absolute error， MAE)等3个参数衡量评估模型精度[11]。R2越接近1， 模型的拟合程度越高， RMSE和MAE越小， 模型的预测能力越高， 模型鲁棒性越高。

3 实验部分

3.1 土壤采集制备和理化分析

在洪湖市燕窝镇内进行土壤采样， 该区域土壤为黄棕壤， 定点采集耕层表层土壤土样29份， 深度在0～20 cm， 将土样用聚乙烯袋密封干燥保存。 然后将土样利用搪瓷盘摊开盛放， 放于通风处自然风干， 将风干土壤剔除石块、 植物残体等异物， 用木棍研磨土壤样品， 首先使用2 mm孔径筛筛选土壤， 取出四分之一继续磨细使通过0.15 mm孔径筛作为最终土壤样品， 研磨后的土壤样本分为两份， 分别用于光谱采集和理化分析测试。 理化分析的土样经硝酸， 盐酸， 高氯酸消煮后， 采用硼氢化钾-硝酸银分光光度法进行实验测定， 每份土样测量3次， 从测量结果的算术平均值得出土壤中的重金属As含量。

3.2 土壤光谱测定

采用SVC HR-1024型地物光谱仪对土壤样品进行暗室光谱的测定， 波长范围为350～2 500 nm， 光谱分辨率为： 350～1 000 nm为1.5 nm, 1 000～1 900 nm为3.8 nm, 1 900～2 500 nm为2.5 nm, 总波段数为990。 室内光谱测量具有易控制、 数据结果质量较高等特点。 将土壤放于通体漆黑的培养皿中， 并利用钢尺刮平表层土壤。 在暗室中， 用1 000 W的卤素灯作为光源， 照射方向与垂直方向呈45°夹角， 光源距离土样表面约为30 cm， 光谱仪探头垂直于土壤表层上方， 探头距离土样约为10 cm。 每次测量光谱之前， 需对光谱仪进行白板校正。 为消除测量数据误差和其不稳定性， 每个土样旋转三次， 每次转动90°， 以获得四个方向的反射率， 有4条光谱曲线， 其平均值作为各土壤样本的原始反射率光谱值， 如图3所示。

3.3 数据处理和分析

3.3.1 光谱预处理

由于仪器自身的原因和光谱采集过程中不可避免的受到测试环境、 样品背景、 观测角度、 样品粗糙度、 杂散光等因素的影响， 光谱曲线边缘波段噪声比较大。 为降低外界噪声对实际建模的干扰， 每份土样都去除噪声较大的边缘波段350～399和2 400～2 500 nm， 保留400～2 399 nm波段用于建模分析， 如图4所示。

3.3.2 土壤As含量分析

对29个样点土壤重金属砷含量进行了统计分析， 结果如表1所示。 参考《土壤环境质量标准GB15618—1995》， 洪湖燕窝镇采集的土壤重金属元素As含量的平均值超过三级质量标准(为保障农业生产和植物生长的土壤污染临界值)， 表明该地区土壤的砷污染达到中度和重度污染。 砷对于农作物的生长危害十分严重， 可以通过植物根茎部分， 利用植物生长吸收作用转移到植物各部位， 初期表现为叶片枯萎， 中期为根系停止发育， 最后整株植物枯萎死亡。 砷对于人体而言同样危害巨大， 可以通过食物链累积， 迅速达到临界值， 它可以使红血球溶解， 破坏正常细胞， 并且具有遗传性、 致癌性和致畸性等， 所以对于重金属砷的监测具有重要意义[12]。

图3 原始光谱

图4 去除边缘噪声

表1 各数据集As值的描述性统计

在进行建模之前， 需要将样本分组， 一组用于模型的构建， 即校准集； 另一组用于验证模型的预测能力， 即验证集； 如表1所示。 采用浓度梯度法， 最终选择21个校准样本， 8个验证样本， 校准集和验证集的均值， 标准差和变异系数均较为接近。 因此， 这样划分比较合理， 可以用于后续建模。 对于建立一个满意的预测模型， 校准采样算法和校准集的大小对整组数据集的代表性有着深刻的影响[13]。

4 结果与讨论

4.1 特征提高型CARS特征波段选取

特征波段选取流程如图5所示， 首先在特征粗选阶段利用CARS对暗室实测光谱进行粗选， 设定蒙特卡罗采样为100次， 从图6可以看出经过交叉检验(CV)所得到的RMSECV的变化趋势图, RMSECV的变化曲线先由大到小至最小后再逐渐变大的变化， 当采样循环次数为60时， RMSECV值达到最小值， 表明此时模型精度最高， 而随着采样次数变大， RMSECV上升， 模型效果变差。 因此， 当采样次数为60， RMSECV值达到最小值时， 所选择波段子集为最优。 但在CARS算法过程中， 采用偏最小二乘法为适应度函数, 这种线性模型面对高维数据， 会出现精度不高等问题，且CARS算法未改变原始数据， 面对相关性较低数据， 往往无法使用， 所以有必要提高特征。 在特征提高阶段， 从图7可以看出CARS筛选出来的原始数据与重金属As之间的相关性较低， 而经过高斯滤波、 一阶导数， 再次高斯滤波后相关性大大提高。

图5 CARS特征波段选择流程图

图6 CARS关键变量选择

图7 不同预处理情况下相关系数变化

在特征精选阶段， 利用皮尔逊相关系数(pearson correlation coefficient, PCC)求取预处理后的光谱指标与土壤重金属As之间的相关系数， 获取绝对值相关性大于0.6的波段作为特征波段。

4.2 回归预测

特征提高型CARS特征波段选取选择了7个特征波段用于建模， SVM模型的预测效果与核函数σ和正则化参数γ有密切的关系， 极大程度影响着预测精度， 如何通过优化获取最优参数是非常重要的步骤。 因此选取了基于PSO优化SVM算法来选取最优参数获得最佳模型， 在PSO优化模型中， 粒子群的种群数N=30， 最大的迭代次数Gmax=300， 学习因子C1和C2为1.8和1.2， 核函数形状参数σ取值范围为[0.1～1 000]， 正则化参数γ取值范围为[0.01～100]。 分别利用每个粒子对所对应的实验样本的光谱和重金属含量进行预测， 得到各粒子当前均方根误差(RMSE)， 将目前最优适应度与粒子本身最优适应度进行对比， 如果目前最优适应度优于粒子本身最优适应度， 则更新速度和位置， 最终获得核函数为径向基核函数(RBF)， 核函数形状参数σ为100， 正则化参数γ为0.357 3。 将最优核函数和正则化参数代入到SVM模型中， 用于土壤重金属As含量高光谱预测。 从图8中可以看出， 适应度值逐渐向最优值接近， 一旦达到最优值后就趋于稳定， 利用PSO-SVM模型可以有效预测土壤重金属As含量， 其中R2， RMSE， MAE达到0.982 3， 0.521 6和0.416 4， 再利用SVMR和RFR回归作为对比， 从图9可以看出验证集基本上呈现线性均匀分布， 表现出较好的预测水平和稳定性， 满足了实际预测要求。 其中SVMR表现最差， 偏离较大， 而PSO-SVM离1∶1线偏离最小， 拟合程度最高。

图8 适应度曲线变化

对比传统化学分析法， 利用土壤反射率高光谱反演土壤重金属As含量的方法， 有效率高、 易操作、 环境友好型等优点， 对于前期生态环境监测和后续生态修复提供数据支持具有重要意义。

图9 实测值与预测值关系

研究中分别利用三种方法建立了高光谱重金属As含量预测模型， 从表2可以看出SVMR预测效果较差， 而PSO-SVM验证集的反演精度表现最优。 PSO优化SVM的回归模型精度上有了极大地提高， 利用优化得到的核函数σ可以保证模型反演精度， 正则化参数γ可以使模型具有较好的泛化能力， 为预测其他重金属元素提供了参考依据。

表2 算法精度对比

5 结 论

以江汉平原典型区域洪湖市燕窝镇污染土壤为研究对象， 采集总计29个土壤样本， 21为校准集， 8个为验证集， 基于三种不同模型构建方法， 重点讨论了利用土壤反射率光谱可以较好的反演土壤重金属As含量， 主要结论如下：

(1) 特征提高型CARS算法， 可以有效剔除无关冗余信息， 极大地提高了光谱与重金属As含量之间相关性， 且选取波段少， 模型简单， 大大提高了效率。

(2) 对比支持向量机回归(SVMR)、 随机森林回归(RFR)等建模方式， 采用PSO-SVM通过迭代更新个体极值和群体极值， 获得最佳参数， 具有较高的精度。 PSO-SVM的验证集的R2为0.982 3、 RMSE为0.521 6、 MAE为0.416 4。 说明模型具有较高的稳定性和预测精度， 可以满足实际预测需求， 为今后大规模反演洪湖地区土壤重金属含量提供依据。 此外， 区域构建高光谱模型也可为无人机、 航空、 航天遥感提供参考。