径向基神经网络的苏打盐碱地重金属定量反演

付艳华， 刘 晶， 毛亚纯， 曹 旺， 黄家其， 赵占国

1. 东北大学江河建筑学院， 辽宁 沈阳 110819

2. 东北大学资源与土木工程学院， 辽宁 沈阳 110819

3. 中国黄金集团， 北京 100000

引 言

近年来， 随着我国城市化进程加快， 城市周边耕地被大量占用， 人均耕地日益减少， 对盐碱地的改造已成为迫在眉睫的问题。 松嫩平原是世界三大苏打盐碱地分布区之一， 该区域盐碱地总面积约为3.42×106hm2， 70%以上为苏打盐碱土 (主要成分为NaHCO3和Na2CO3)[1]。 由于盐碱地改良及农作物选取方法与盐碱地中重金属含量密切相关， 因此如何准确、 快速确定大区域内盐碱地重金属含量已成为亟待解决的关键问题。 遥感高光谱技术具有快速、 动态获取地物连续光谱信号的优势， 已被广泛应用于河漫滩平原[2]、 矿区[3]、 农田[4]的土壤重金属定量反演等领域。

有效的光谱数据预处理方法和反演方法是实现盐碱地土壤重金属高精度反演的重要保障。 在光谱数据预处理方面， 涂宇龙等[5]在保留土壤样品弱光谱信号的基础上， 发现经标准正态变换的光谱全要素主成分分析后基于逐步回归建模方法能有效提升土壤Cu含量估算精度； 徐丽华等[6]利用Savitzky-Golay平滑(SG)、 多元散射校正(MSC)、 标准正态化(SNV)和克里克滤波(KF)等8种方法对33个水稻土壤样本进行了光谱预处理， 然后基于分区极值法选取的特征波段进行建模， 结果表明建模精度明显提高。 在建模方法方面， 程先锋等[7]开展了兰坪矿区土壤样本的Zn， Pb， As和Cd共4种重金属含量与光谱响应测试实验， 利用逐步回归方法预测土壤重金属含量的反演模型具有较高的精度，R2均高于0.816。 Ma等[8]通过对土壤重金属含量与预处理波段的相关性分析， 建立了土壤重金属含量预测模型， 结果表明， 极限学习机和随机森林的建模精度优于支持向量机。 目前， 国内外学者利用遥感技术对农田、 河漫滩平原土壤的重金属含量研究较多， 但针对盐碱地重金属含量的反演研究较少， 尤其针对该应用中的光谱数据预处理和建模训练样本选择等方面研究不足。

以吉林省白城市镇赉县苏打盐碱地为研究区域， 首先利用Savitzky-Golay (SG)平滑法、 多元散射校正(MSC)和连续统去除法(CR) 3种方法对原始光谱数据进行了预处理， 然后利用差值、 比值及归一化对预处理后的数据提取了光谱指数， 并将光谱指数与重金属含量做相关性分析以选取高相关性的光谱指数， 再以选择的光谱指数作为建模输入的自变量， 盐碱地重金属含量作为因变量， 运用径向基神经网络算法建立重金属元素锰(Mn)、 钴(Co)和铁(Fe)含量的反演模型， 最后通过相关系数等梯度循环建模精度分析确定了光谱指数的最优选取原则和重金属含量的最优反演模型。 研究结果表明： 利用上述方法对Mn， Co和Fe含量进行预测， 相对精度分别达到88.64%， 90.36%和91.78%， 为盐碱地土壤重金属元素含量的反演提供有效技术手段和理论支撑。

1 实验部分

1.1 研究区域概况

镇赉县位于中国吉林省白城市的西北部， 总面积约4 737 km2。 镇赉县百万亩水田具有盐碱地面积大、 可溶性盐含量高等特点， 且锰、 钴、 铁等重金属含量较高。 由于重金属含量决定盐碱地土壤改良方法及种植植被种类， 因此以镇赉县为研究区域， 采样区面积约为15 km2， 每个采样点之间有400 m， 如图1所示。 以均匀抽样的方式共采集65个土样作为实验样本， 开展盐碱地重金属含量反演研究。

图1 研究区概况及采样点布置图

1.2 实验数据获取

首先通过实验室土壤干燥箱对采集的样本进行烘干处理， 然后去除石块和植物根系后研磨， 放在土壤干燥箱中烘干， 过100目筛， 共制成粉末状样本65件， 每个样品分成两份， 每份约150 g， 分别用于化学分析和光谱实验。 采用美国SVC HR-1024便携式地物光谱仪对65件粉末状样品进行光谱测试。 该仪器的基本参数如表1所示。

表1 SVC HR-1024便携式地物光谱仪基本参数

光谱测试实验于10:00—14:00进行， 该时段可有效降低气溶胶及太阳辐射传播路径对光谱测试的影响。 将被测样品放置于直径为6 cm的圆形黑色小盒中， 保证样品表面平整， 可有效避免测量背景对光谱实验的影响。 光谱仪镜头垂直于样品观测面， 取3次测试的反射率均值作为样品最终反射率。

光谱测试结束后， 对样品进行化学测试确定各个实验样本的锰、 钴和铁的含量， 实验土壤样本中重金属含量描述统计分析如表2所示。

表2 重金属含量描述统计分析(mg·kg-1)

1.3 光谱特征分析

由于边缘波段测试效果不佳， 因此在分析中去除了350～399和2 401～2 500 nm两段数据， 最终获得890个波段。 65件土壤样品的可见光-近红外波段光谱反射率曲线如图2所示， 由图可知， 光谱反射率介于5%～60%之间， 在400～2 400 nm之间呈现一定的波动性； 在1 400， 1 900和2 200 nm波段附近水分吸收特征较为明显。 在400～1 000 nm之间， 光谱反射率随着波长的增加呈现迅速增大趋势， 最高值达到47%； 在1000～1 400 nm之间， 光谱反射率仍保持随着波长的增加而增大， 但增加曲线的斜率逐渐变小； 在1 400和1 900 nm这两个水分吸收带之间， 光谱曲线仍然呈现上升的趋势， 趋势较为平缓， 累计增幅约为5%； 在2 100 nm附近， 反射率总体呈下降趋势。

图2 实验样品的可见-近红外光谱曲线

对于土壤光谱而言， 其吸收特征与某些特定的土壤属性有关。 在400～600 nm区域形状较陡， 斜率较大， 这与土壤中所含的铁有关。 在可见光(400～700 nm)和短近红外(700～1 000 nm)波段范围， 土壤光谱的吸收特征主要由于金属离子(如Fe2+， Fe3+和Mn3+等)的电子跃迁形成的， 且Mn的主要吸收波段为588.5～591.4， 595.7～601.4， 604.3和615.7 nm[9]。

2 结果与讨论

2.1 光谱数据预处理

由于光谱数据在获取过程中易受到环境变化、 仪器稳定性的影响， 因此获得的光谱数据中包括地物光谱数据和一定的噪声， 而通过光谱数据预处理可以有效降低噪声的影响， 并增强土壤原始光谱的细节特征， 更好地识别地物和土壤重金属信息， 增强反演模型的鲁棒性[5]。 因此， 对原始光谱数据进行Savitzky-Golay平滑、 多元散射校正、 连续统去除3种预处理变换。

(1)Savitzky-Golay平滑法(SG)

Savitzky-Golay算法是最早在1964年引入的基本平滑方法之一， 它利用多项式来对移动窗口内的数据进行最小二乘拟合， 其实质是一种加权平均法。 该算法不受样本数据限制， 适用于各种信号的平滑去噪, 能够保留频谱的峰谷特征。 与传统算法相比， 该算法具有更稳定、 误差更小的平滑去噪效果[10]。 计算公式如式(1)所示

(1)

(2)多元散射校正(MSC)

多元散射校正算法[11]常用于消除因样本颗粒分布不均匀造成的散射影响, 在农业土壤研究领域的应用较为广泛。 该方法可有效提高原始光谱的信噪比， 消除光谱数据的线性散射干扰。 具体处理过程如下：

计算待校正光谱的平均光谱

(2)

一元线性回归

(3)

得出多元散射校正

(4)

(3)连续统去除法(CR)

连续统去除法又称去包络线法， 可有效去除光谱数据中的不相关信息， 是一种可有效增强光谱特征的预处理方法[13]。 连续统去除法可以将光谱数据归一化， 不仅能保留和增强光谱的吸收特征， 还能减轻混合物引起的非线性影响， 以便于进行光谱吸收特征分析和光谱特征波段选择。 其公式如式(5)所示

Rc(λ)=Rcr(λ)-R(λ)

(5)

式(5)中，Rc，Rcr和R分别是去包络(特征吸收)、 包络线和光谱反射率值，λ是波长。 经过对光谱反射率曲线去包络后， 能够清晰地看到特征吸收峰。

2.2 光谱指数构建

为提高重金属含量反演精度， 分别构建了比值(RI)、 差值(DI)和归一化差值(NDI)三种光谱指数[9]。 对预处理后的光谱数据所建立的指数与重金属含量作相关性分析， 各指数计算公式如式(6)—式(8)

(6)

DI=Rλm-Rλn

(7)

(8)

式(6)—式(8)中：λm为m点对应的波长；λn为n点对应的波长；Rλm为m点波长对应的反射率值；Rλn为n点波长处对应的反射率值。

利用3种光谱预处理方法处理后的光谱数据构建RI， DI和NDI三种光谱指数， 并分析重金属含量与光谱指数的Spearman秩相关系数。 对锰原始光谱数据的DI、 连续统去除后的DI、 Savitzky-Golay平滑后的DI和多元散射校正后的DI与锰重金属含量相关性分布如图3所示。

图3中的横坐标和纵坐标是样品的光谱波长， 颜色表示光谱指数与该点样品的Mn含量之间的相关系数的绝对值。 颜色由蓝色逐渐增加到红色， 说明二者之间相关系数的绝对值由0增加到1， 红色越深， 说明光谱指数与重金属含量的相关性越强。 图3中(a)， (b)， (c)和(d)分别为原始光谱指数DI、 CR处理后的光谱指数DI、 SG平滑处理后的光谱指数DI、 MSC处理后的光谱指数DI与Mn含量的相关性。 比较上述三种预处理方法， CR处理后光谱指数DI与Mn的相关系数基本低于0.5， 勉强达到弱相关， 光谱预处理效果较差； SG平滑处理后的光谱指数DI与Mn的相关系数可达到0.8， 能达到强相关， 预处理效果较为理想； MSC处理后的光谱指数DI与Mn的相关系数和原始光谱数据的光谱指数DI与Mn的相关系数基本一致。

图3 三种预处理方法结合差值指数与锰含量的相关分布图

2.3 模型建立与验证

各种数据处理方法增强了部分光谱特征， 但也残留了部分冗余信息， 数据处理后的光谱指数选取过多或过少， 均会造成建模精度的降低， 因此需要确定一个合理的数据集来保证建模的反演精度， 具体流程如图4所示。

图4 土壤重金属含量反演建模流程图

将上述共计65件样本， 按照重金属含量从低到高排列， 抽样选取50件样本作为训练集， 15件样品作为测试集。 训练集的光谱指数按照秩相关系数r>τ进行等梯度选取输入数据集， 利用径向基神经网络算法对重金属含量进行建模预测， 输出预测的重金属含量。τ的初始取值为0.6， 然后按Δτ=0.1的梯度递增， 选取特定输入数据集进行多次循环建模预测， 在分析中为验证所建模型的精度， 以决定系数R2、 均方根误差RMSE和平均相对精度作为模型评价指标。R2越接近1， 均方根误差RMSE值越小， 平均相对精度越接近1， 表明模型精度越高[14]。 基于反演精度分析最终确定重金属锰、 钴和铁含量的最优选取原则和最优反演模型。

根据建模后的决定系数R2、 均方根误差RMSE和平均相对精度这三个评价指标， 确定了最优选取原则， 如表3所示。 Mn选取相关系数r>0.70， Co选取相关系数r>0.80， Fe选取相关系数r>0.80， 并分别选取了108组、 690组和31组。 根据最优选取原则， 对比三种光谱指数(DI， RI和NDI)， 其中RI和NDI被选取的组数多， 表明这两种光谱指数算法可有效增强光谱特征， 提升建模精度。

表3 重金属元素光谱指数最优选取原则

根据最优选取原则选取后光谱指数组合作为输入数据， 利用径向基神经网络建立了最优重金属含量反演模型， 其预测值与实际重金属含量对比如图5， 图6和图7所示， 其中Mn， Co和Fe的预测值与实测值的回归曲线R2分别为0.703 4， 0.897 6和0.848 4， 均方根误差RMSE分别为53.007 3， 1.059 2和0.363 4， 平均相对精度达到88.64%， 90.36%和91.78%。

图5 径向基神经网络预测的Mn含量与实测Mn含量对比图

图6 径向基神经网络预测的Co含量与实测Co含量对比图

图7 径向基神经网络预测的Fe含量与实测Fe含量对比图

3 结 论

以中国吉林省白城市镇赉县盐碱地65件样本的土壤可见光-近红外光谱数据和Co、 Fe和Mn元素的含量为数据源， 采用了Savitzky-Golay平滑、 多元散射校正、 连续统去除共3种数据预处理方法对原始光谱数据进行了处理， 并构建了差值指数、 比值指数和归一化指数， 同时提出了适用于Mn， Co和Fe含量反演的光谱指数最优选取原则。 然后， 利用径向基神经网络算法建立了盐碱地重金属Mn， Co和Fe含量最优反演模型， 结论如下：

(1)对比三种预处理方法， 利用MSC算法对盐碱地土壤光谱数据进行预处理、 对光谱数据中干扰信息的去除效果最为明显。 对比三种光谱指数算法， RI和NDI两种光谱指数算法可有效增强光谱特征， 提升建模精度。

(2)根据提出的光谱指数最优选取原则， 基于径向基神经网络建立的盐碱地重金属Mn， Co和Fe反演模型预测精度较高， 其预测值与实测值的平均相对精度分别为88.64%， 90.36%和91.78%。

研究对盐碱地重金属含量精确快速反演提供了一种有效方法。 对含重金属的盐渍化土壤治理具有重要的现实意义。