基于铁氧化物特征光谱和改进遗传算法反演土壤Pb 含量

张 霞，王一博,2，孙伟超，黄长平，张 茂,2

（1. 中国科学院空天信息创新研究院，北京 100101；2. 中国科学院大学，北京 100049）

0 引 言

土壤为自然生态系统的运转和人类活动提供了基本保障。随着工业化的不断发展，通过污水排放、大气降尘等途径进入土壤的重金属元素在土壤中持续累积，造成土壤重金属污染[1]。Pb 是土壤中常见的一种重金属污染物，渗透进土壤后移动性较差，残留时间长，且不易被微生物分解，其治理和修复难度大，不仅影响农作物的生长，还会通过食物链威胁人类健康[2]。因此，广泛开展土壤Pb 含量调查和监测，对保障农业生产、保护人类健康和生态系统安全具有重要意义。近年来高光谱遥感技术[3-4]不断发展，其具有光谱分辨率高、波段多且连续性强等特点[5]，能够快速高效地反演土壤重金属含量[6-8]。另外，随着高光谱图像性能的不断提高，土壤重金属含量的大范围制图以及周期性监测成为可能。

对土壤光谱进行分析时，国内外学者常常使用全谱段来对土壤重金属的含量进行预测[9-13]。一般情况下，自然环境中的土壤重金属含量较低，其本身微弱的光谱响应信号又受到多种成像因素的干扰，难以直接通过分析土壤重金属元素的光谱特征来估算其含量[14]。但土壤光谱活性物质（主要指黏土矿物、铁氧化物和有机质等）对重金属元素具有吸附或赋存作用[15]，因而土壤重金属元素对土壤光谱的影响间接反映在这些光谱活性物质的特征光谱波段上。Sun 等[16-17]使用有机质和黏土矿物的特征波段组合建立Zn 和Ni 的反演模型，与全谱段建模相比精度有显著提升。贺军亮等[18]基于土壤有机质敏感波段对应的多种光谱变换指标，建立了Cd 的高光谱间接反演模型。虽然现有研究在间接反演土壤重金属元素含量方面取得了一定进展，但关于提取土壤 Pb 相关光谱活性物质的特征光谱进行Pb 含量反演的可行性还有待验证。

在建模算法的选择方面，遗传算法-偏最小二乘回归（Genetic Algorithm-Partial Least Squares Regression，GA-PLSR）[16-17,19-20]具有特征优选的能力，且能对具有强相关性的自变量进行建模，常被用于高光谱数据模型构建。由于高光谱数据波段众多、相关性强且冗余性高，使用随机寻优的遗传算法对参与建模的波段进行筛选能够减少自变量的个数，降低数据间相关性，提取出光谱的有效信息，最终提高模型估算精度[21-22]。但是遗传算法仍然存在“过早收敛”的问题[23]，算法在迭代时只接受适应度更高的解，所以在开始迭代后便立即向初始种群的局部解收敛，从而容易陷入局部解区域，难以做到全局寻优，对最终模型的精度产生了一定的影响。

以往研究多以矿区作为研究对象，而矿区以外的一般农作区也应是监测的重要方面。本文以河北雄安一般农作区农田采集的土壤样本数据为研究对象，提取土壤中对Pb 起主导吸附作用的土壤光谱活性物质的特征谱段进行Pb 含量反演，以避免全谱段建模的数据冗余问题。为解决GA 算法的“过早收敛”问题，提出改进遗传算法（Improved Genetic Algorithm，IGA），使算法在运行前期能够跳出局部解区域，寻找更优解，并结合偏最小二乘法[24]进行回归建模反演土壤Pb 含量。

1 研究数据与预处理

1.1 研究数据

雄安新区是中国新规划的数字经济创新发展试验区，为关注该地区的环境建设，将其作为研究区。2018 年9 月20 日至9 月23 日，于雄县和安新县的农田开展土壤样本采集与野外原位光谱测量试验，采集农田土壤表层（0～20 cm）作为土壤样本。野外光谱采集使用SVC HR1024i 地物光谱仪，该光谱仪波长范围为 350～2 500 nm，测量时探头离地面距离大约5 cm。采集土壤光谱时每个样点测3～5 条光谱取平均作为样点光谱，共采集了土壤样本及光谱73 组。

土壤样本送实验室进行 Pb 含量测试。测量方法为HNO3-HF-HClO4 消煮法，将样本在微波消解仪进行消解，然后用电感耦合等离子体质谱ICP-MS 测定Pb 元素含量，每个土壤样本测量3 次后取平均值。

1.2 数据预处理

光谱数据处理主要包括光谱维去噪、异常光谱剔除、光谱特征增强及铁氧化物特征谱段提取等。

野外光谱测量时，大气中的水汽在1 400 和1 900 nm处存在强吸收，且1 900～2 500 nm 区间的光谱存在较大噪声[25]。为去除大气水汽吸收波段和受噪声影响的波段，并保留尽可能多的光谱区间，将1 800 nm 之后的光谱波段从土壤野外光谱中剔除。为进一步去除野外光谱在 350～1 800 nm 的随机噪声，使用包络线去除及分段Savitzky-Golay（SG）滤波对光谱曲线进行平滑处理。SG滤波中890～1 020 和1 330～1 520 nm 采用窗口大小为15的二次多项式，其余区间采用窗口大小为7 的二次多项式。为排除光谱异常和重金属含量异常对试验结果的影响，通过光谱曲线筛选和z-score 异常点检测方法，去除3 个异常样本，最终保留70 个有效样本用于反演建模。

高光谱数据相邻两波段具有较强的相关性，使得整个波段的反射率具有强多重共线性，给反演带来极大困难。因此采用一阶差分的方式，减弱各波段之间的共线性，增强光谱的特征[4]。差分公式见式（1）

式中Bi′为第i段光谱的一阶差分，Bi为第i段光谱反射率，n为光谱总数。

2 研究方法

2.1 铁氧化物光谱特征提取

Covelo 等[26]的重金属吸附和解析试验表明：土壤对重金属Pb 的吸附作用最强，且土壤有机质、铁氧化物、锰氧化物均对重金属Pb 具有很强的吸附作用，其中铁氧化物在重金属 Pb 的吸附作用中占主导地位，而且随着土壤中不同重金属之间对吸附位置的竞争加剧，铁氧化物在重金属Pb 的吸附作用中的重要性也随之上升。Wu等[27]的研究表明，铁氧化物与重金属含量相关性在所有土壤光谱活性物质中占主导地位，并试验得出铁氧化物的吸收特征在以500 和950 nm 为中心的吸收峰。因此，本文提取500 和950 nm 处的吸收特征谱段预测土壤Pb含量，提取波段区间为[Bm-W/ 4 ,Bm+W/4]，其中Bm是最大的吸收波段，W是吸收区域的宽度。

综上，从70 条有效土壤光谱中提取500 和950 nm为中心的铁氧化物吸收峰，提取出的光谱波段范围为450.7～523.2 以及 9 14.2～1 027.9 nm。

2.2 改进的遗传波段选择算法（IGA）

遗传算法[28-30]仿照物种演化机制，模仿了遗传进化中发生的基因突变和自然选择现象，通过随机选择、交叉和变异等遗传操作，产生适应度高的个体作为最终解。传统遗传算法中存在最严重的问题是“过早收敛”：运行初期会直接向局部解收敛而难以跳出。为解决该问题引入模拟退火算法中的Metropolis 准则[31-32]：若新解优于当前解，则接受新解，否则以某一概率接受新解，且此概率随适应度增大而降低。Metropolis 准则使得算法在运行初期有一定的概率接受较差解而跳出当前局部解区域，从而做到全局寻优。IGA 算法具体流程如图1。

图1 改进的遗传算法流程图Fig.1 Process of improved genetic algorithm

算法的改进核心为生成新种群，其中具体的Metropolis准则为：当新个体目标函数值小于父个体目标函数值时，用新个体取代父个体；否则以 exp （ -k⋅ Δf/E）的概率接受新个体为父个体。其中k为大于0 的实数，E为种群目标函数值的方差。令n为算法总迭代次数，当E大于0.01/n时，k=E+ 0 .01/n；否则k= 0 .01/n。k保证了E大于0.01/N时算法有足够概率能够跳出局部解区域，E小于0.01/N时算法能够保持收敛状态。E定义如式（2）

式中N为种群个体数，fi为个体i的目标函数，f为种群目标函数的平均值。

2.3 偏最小二乘建模与模型评价指标

偏最小二乘法（Partial Least Squares Regression，PLSR）能较好地解决高光谱反演中自变量之间的强多重共线性问题，并且能够实现多对多线性回归建模，特别是在自变量存在多重相关性而且个数较多，而样本数目又远小于自变量的情况下PLSR 仍适用[33-35]。因此，本研究采用PLSR 建模。

使用拟合优度R2，均方根误差（Root Mean Square Error，RMSE）和相对偏差（Relative Percent Difference，RPD）评价模型的优劣。较好的模型通常拥有较高的RPD 和R2及较低的RMSE。模型的评价参考现有的土壤属性含量高光谱估算的评价标准[19,35]：出色模型，R2> 0.9，RPD>3.0；良好模型，0.9 >R2> 0.82，3.0 > RPD > 2.5；近似模型，0.82 >R2>0.65，2.5 > RPD > 2.0；具有一定估算能力，0.65 >R2> 0.50，2.0 > RPD > 1.5；不具备估算能力，0.50 >R2，RPD < 1.5。

3 结果与分析

3.1 样本集划分

将 70 个样本分为训练集和测试集。依据重金属 Pb的含量从小到大进行排序，对样本集按照2∶1 的比例进行分层抽样，获得47 个样本的训练集和23 个样本的测试集，样本集的数据统计见表1。

表1 样本集中Pb 含量数据统计Table 1 Data statistics of Pb content in the sample set

3.2 模型构建

为验证 IGA 波段选择算法改进的有效性，本文对GA-PLSR 和IGA-PLSR 这2 种算法进行建模。2 种算法的参数设置相同，个体数为 20，个体编码方式为二进制编码，代际差为90%。将RMSE 作为波段选择算法的目标函数，即算法目标是最小化 RMSE。将选择算法运行结果中RMSE 最小的个体作为当前最优解输出，并将其作为PLSR 的自变量建立最终的回归预测模型。最后使用测试集数据计算模型评价指标。

3.3 土壤Pb 含量估算模型评价及分析

根据2.1 节的理论与方法，提取出铁氧化物特征谱段（450.7～523.2 nm 和914.2～1 027.9 nm）作为输入，使用GA-PLSR 及其改进算法IGA-PLSR 对土壤Pb 含量进行反演建模。为对比两算法的运行效果，选取迭代次数为100、250、500 和1 000，每种迭代次数下运行5 次算法，评价参数及运行时间取5 次运算的平均值。

试验结果如表 2 所示，可以看到随着迭代次数的增加，2 种算法的平均耗时增加，同时平均精度也在不断提高，IGA-PLSR 因加入了Metropolis 准则的筛选过程，相同迭代次数下耗时略大于GA-PLSR 算法。由模型精度随迭代次数的变化可得，2 种算法的最高精度对应的迭代次数均为1 000 次，这与Wang 等[36]应用GA 反演土壤重金属的迭代次数设置一致。在迭代次数为 1000 次时，GA-PLSR 算法的反演精度R2和 RPD 平均值分别为0.782、2.117，RMSE 的平均值为2.487 mg/kg，满足了近似模型的条件，能够近似估计Pb 含量；而IGA-PLSR 的3 种评价参数均要优于改进前的算法，模型的R2和RPD为 0.822、2.377，RMSE 为 2.221 mg/kg，其中R2已经达到了良好模型的标准，表明模型对土壤Pb 含量具有良好的估算性能。因此，IGA 算法相对于传统的 GA 算法具有更优的波段选择能力，显著提升了反演模型的精度。

表2 基于铁氧化物特征谱段建模的反演精度Table 2 Inversion accuracy based on iron oxide characteristic spectrum modeling

图2 分别展示了2 种算法在1 000 次迭代时的5 次运算中的最优结果。可以看到IGA-PLSR的R2、RPD为0.837、2.403，RMSE 是 2.063 mg/kg，显著优于传统算法 GA-PLSR（R2、RPD 分别为 0.788、2.140，RMSE 为 2.447 mg/kg），且通过散点的分布可以看到IGA-PLSR 模型中的样本点更加集中于1∶1 线附近。因此，IGA-PLSR 算法的反演能力显著优于改进前的GA-PLSR 算法。

图2 最优Pb 含量估算结果Fig.2 Optimal results of Pb content pridiction

3.4 铁氧化物特征谱段的有效性分析

为对比分析铁氧化物特征谱段的有效性，基于350～1 800 nm 的全部光谱波段对土壤重金Pb 含量进行建模反演，根据3.3 节模型分析结果，设置迭代次数为1 000 次，表3 展示了试验结果。

表3 基于全谱段建模的反演精度(迭代次数设为1 000 次)Table 3 Inversion accuracy based on full spectrum modeling(1 000 iterations)

对比表2（铁氧化物特征谱段）和表3（全谱段）的结果可以看出，2 种建模方法在采用全谱段建模时模型精度较低。GA-PLSR 模型和IGA-PLSR 模型的平均R2分别为 0.292、0.440，平均 RPD 分别为 1.216、1.366。IGA-PLSR模型有较高的精度和反演能力，表明IGA 波段选择算法在光谱波段较多、冗余度较高的情况下依然能选择出相对有效的光谱波段。但此时 2 种模型均未达到估算的精度要求，而表 2 中使用铁氧化物特征谱段建模的结果要显著优于表 3 全谱段的结果，说明铁氧化物特征波段的提取能切实有效地提高Pb 含量反演的精度。

为直观显示铁氧化物特征波段对土壤Pb 含量预测的贡献度，图3 展示出全谱段IGA-PLSR 模型中IGA 算法的波段筛选结果。波段选择算法的作用是筛选出能有效反演土壤重金属含量的波段，可以看到筛选的波段在铁氧化物的特征谱段（500 和950 nm 吸收峰）分布最为密集，表明对研究区土壤Pb 含量进行建模反演时，铁氧化物特征谱段相比其他波段更加有效。

以上试验均说明了提取相关土壤活性物质的光谱特征能够在一定程度上减少光谱信息的冗余，从而提高Pb含量反演的精度。本试验为基于铁氧化物特征波段反演土壤Pb 含量提供了试验依据，同时为反演土壤中其他重金属含量提供方法参考。

图3 全谱段IGA-PLSR 模型波段筛选结果Fig.3 Band selecting results of the IGA-PLSR model with full spectrum

3.5 算法收敛性及不确定性分析

为进一步分析 IGA-PLSR 算法对“过早收敛”的改进，选择GA-PLSR 和IGA -PLSR 迭代1 000 次时的模型，绘制代价函数 RMSE 随迭代次数变化的曲线，如图 4。该曲线图将算法寻优过程可视化，以便对比 2 种算法运算过程中的异同。向局部最优解收敛的特点为缓慢逼近，由图4 可看到GA-PLSR 模型由于仅接受更优解，整体呈单调下降趋势，在250 次迭代后RMSE 基本稳定，随后400、600、720 及900 次迭代邻域中都仅仅在某一次迭代后RMSE 值有极小的降低，表明已陷入局部最优，后续仅是逼近局部最优解；IGA-PLSR 算法在迭代初期有较高的概率接受较差解，RMSE 波动较大，在 300 次后开始收敛，每次迭代都有一定的概率接受较差解，从而跳出当前的局部解区域，如在380 和640 次迭代附近接受较差解后RMSE 都会有大幅降低出现，说明跳出了当前局部解区域从而进一步寻找更优解。RMSE 曲线图的对比表明了 IGA-PLSR 算法提高反演精度的原因，凸显了其在全局寻优方面的优势。

改进前后的GA 和IGA 都是随机寻优算法，故存在一定的不确定性。主要体现在随机产生初始解、随机选择、交叉和变异等一系列遗传操作中，同时引入的Metropolis 准则中依概率接受较差解也增加了模型的不确定性。因此在实际应用时，应多次运行算法后选择最优的模型进行重金属含量的反演。

图4 不同算法RMSE 随迭代次数的变化Fig.4 Changes of RMSE with iterations of different algorithms

4 结论与展望

本文基于河北雄安一般农作区采集的70 个土壤样本的野外光谱数据，研究具有机理性的特征波段的选取方法对土壤Pb 含量进行反演。首先，为规避全谱段建模的数据冗余和机理性不足问题，提取出土壤中对Pb 具有主导吸附性作用的铁氧化物的特征波段用于Pb 含量反演，在此基础上，通过引入Metropolis 准则对GA 波段选择算法改进并提出改进遗传算法（Improved Genetic Algorithm，IGA），解决遗传算法的“过早收敛”的问题，利用IGA 算法筛选铁氧化物特征谱段中更加有效的光谱波段组合，采用偏最小二乘法（Partial Least Squares Regression，PLSR）构建反演模型。结论如下：

1）IGA-PLSR 模型相对GA-PLSR 拥有更高的反演精度，其最优模型的R2、RPD 分别为0.837、2.403，RMSE为2.063 mg/kg，其中R2达到了良好模型的标准，可用于雄安新区农田土壤Pb 含量的定量估算。

2）铁氧化物特征谱段的反演精度要显著优于全谱段的精度，说明了依据土壤Pb 的吸附机理而提取相关土壤光谱活性物质的特征谱段能够有效减少全谱段建模的冗余信息，从而提高Pb 含量反演精度。

3）通过对2 种算法运行过程中RMSE 的变化曲线分析，表明IGA-PLSR 能够有效地解决GA 算法的“过早收敛”问题，能够跳出局部解区域寻找更加有效的光谱波段组合，提高PLSR 的建模精度。

综上所述，在应用于实际区域土壤Pb 含量快速检测时，可以从两方面提高土壤元素含量反演模型的精度和适用性：提取出与土壤元素相关的光谱活性物质的特征光谱代替全谱段进行建模反演；改进现有 GA 算法使用IGA-PLSR 算法进行反演建模。但是，在野外光谱及航空/卫星高光谱图像获取过程中，不可避免地会受到土壤粒径、含水率等环境因素影响，增加土壤重金属高光谱反演的不确定性，因此，下一步将研究野外光谱乃至像元光谱的环境因素去除方法，以期提高土壤重金属高光谱遥感反演方法的稳定性。