基于SGA-RF算法的农业土壤镉浓度反演研究

王轩慧 陈建毅 郑西来 朱 成 王轩力 单春芝

(1.中国海洋大学海洋环境与生态教育部重点实验室， 青岛 266100; 2.青岛农业大学理学与信息科学学院， 青岛 266109; 3.中国联合网络通信有限公司济南软件研究院项目管理部， 济南 250000; 4.山西工程技术学院信息工程与自动化系， 阳泉 045000; 5.国家海洋局北海环境监测中心， 青岛 266033)

0 引言

农业土壤重金属污染已成为我国主要的环境问题之一[1-2]。土壤中的镉会在农作物中富集并进入食物链，从而对人类健康造成严重威胁[3-4]。因此，密切监测农业土壤中的重金属镉超标情况对于防止农产品污染，保障人体健康有着重要意义。

目前，常规化学方法对土壤重金属浓度的测定不仅周期长，成本高，而且会产生大量危害环境的废弃物[5]。相比之下，可见/近红外光谱技术具有成本低、效率高、环境污染少等优点[6-7]，逐渐成为有效监测土壤重金属污染情况的重要手段[8-9]。CHEN等[10]采用正交信号校正预处理方法和反向传播神经网络建立污水灌溉区的重金属镉反演模型。GHOLIZADEH等[11]结合各种预处理方法与支持向量机回归建立土壤重金属含量与反射光谱之间的对应关系。TAYEBI等[12]使用可见/近红外反射光谱和偏最小二乘(PLSR)模型来预测农业土壤中重金属铁的浓度，指出二阶微分预处理与PLSR结合建立的模型反演精度较高。

尽管上述方法成功地将可见/近红外光谱技术应用于土壤重金属高光谱反演领域，但由于原始光谱数据中存在大量冗余和不相关特征，从而严重影响了反演模型的性能。特别是农业土壤中重金属镉所对应的可见/近红外光谱波段非常少，因而直接使用原始光谱构建回归反演模型不仅准确率和稳定性偏低，而且运算效率较差。因此，只有通过使用特征选择方法才能选出具有最小冗余度且最具代表性的敏感波段，进一步有效地提高农业土壤镉回归反演模型的预测能力[13]。JIANG等[14]利用遗传算法结合PLSR建立土壤镉含量反演模型。虽然上述方法取得了一定成果，但是该方法对相邻光谱波段的特征值进行了融合，这会导致部分潜在重要特征的丢失。XIA等[15]采用相关系数法结合PLSR成功反演长江流域沉积物中的镉含量。但基于相关系数法的特征选择方法只适用于线性统计方法建立的模型，对于非线性相关及样本分布不均匀的情况，反演精度较低。

本文在重金属镉高光谱反演领域，针对上述特征选择技术的不足，提出一种基于斯皮尔曼等级相关分析的遗传随机森林(Spearman’s rank correlation coefficient-based genetic algorithm using random forest，SGA-RF)特征选择算法。基于大沽河流域表层农业土壤可见/近红外高光谱数据，通过与传统遗传算法和基于K近邻的遗传算法进行比较，探讨SGA-RF算法在农业土壤重金属高光谱反演领域应用的可能性，以期为利用高光谱技术检测农业土壤中低含量重金属污染状况提供一定的理论依据。

1 SGA-RF特征选择算法

图1 SGA-RF特征选择流程图Fig.1 Flow chart of SGA-RF feature selection

SGA-RF特征选择算法流程如图1所示。该算法由2个子算法组成，是一种基于Wrapper[16]型的特征选择方法。首先，为了消除冗余和无关信息，采用斯皮尔曼相关分析特征选择方法(算法1)对2 051个预处理后的初始特征集合进行预选，得到了彼此无关联的特征波段子集作为算法2的输入。其次，在遗传随机森林算法(算法2)中采用基于随机森林袋外误差的适应度函数来评价个体的优劣，最终选出适应度函数最小值所对应的特征子集为最优特征波段子集。最后，利用最优特征波段子集构建随机森林回归(Random forest regression, RFR)反演模型。

以斯皮尔曼等级相关分析作为SGA-RF算法的特征预选方法，主要目的是有效减少初始特征集维数并剔除原始光谱特征波段之间的冗余和无关信息。该算法包括：①数据准备阶段(步骤1～步骤2)，通过计算两两特征波段之间的斯皮尔曼等级相关系数将强相关的特征波段存入二维矩阵X中。②算法的核心组成部分(步骤3～步骤8)，通过遍历矩阵X构造所有相关特征的分组Y。③算法输出部分(步骤9～步骤11)，通过对特征波段集合Q与Y中每组代表波段求并集得到最优光谱特征波段子集。算法的具体描述如算法1所示。

(1)算法1——基于斯皮尔曼等级相关分析的特征选择算法

输入：①预处理后的光谱特征集合M。②循环变量i1。③循环变量j1。④循环变量b1。⑤循环变量c1。⑥最大波长数L2 051。

输出：所有相关光谱特征分组集合Y(所有相关的分组)；最优光谱特征子集F。

步骤1: 计算Mi与Mj两个特征波段之间的斯皮尔曼等级相关系数，即

(1)

式中N——样本数量

d——两列成对变量的等级差分集合

在式(1)中,对Mi与Mj进行排序(同时为升序或降序)，得到两个元素等级集合m和n，其中mi为mi在Mi中的等级，ni为ni在Mj中的等级。将集合m、n中的元素对应相减得到一个排行差分集合d，其中di=mi-ni，1≤i≤N。如果ρ的绝对值(|ρ|)大于或等于0.8且标准偏差p小于10%，则特征值Mi与Mj被认为是显著相关的，将Mi与Mj存入矩阵X中，X为二维矩阵。

步骤2:重复执行步骤1，直至j达到最大波长数L。

步骤3:把i赋值给j，重复执行步骤1～步骤2，直至i达到最大波长数L。

步骤4：把X中的第一行两个元素赋值给Yb。设置ismatch为0，a为1。

步骤5:取出矩阵Xa中第a行的第1个元素Xa1。如果矩阵Yb中不含有Xa1，转入步骤6；否则设置ismatch为1，再判断矩阵Y中是否含有Xa2，如果矩阵Y中不含Xa2，将Xa2存入Yb中,转入步骤7。

步骤6:取出矩阵Xa中第a行的第2个元素Xa2。如果矩阵Yb中含有Xa2，设置ismatch为1，并判断矩阵Yb中是否含Xa1，如果矩阵Yb中不含有Xa1，将Xa1存入Yb中。

步骤7:如果ismatch为1，在X中剔除Xa，得到新的矩阵Xa与Yb。

步骤8:重复执行步骤5～步骤7，直到遍历完Xa的所有元素。Yb中存储一组显著相关的特征波段，Xa为去掉与Yb相关特征波长之后的剩余特征波长变量。若Yb中含有不相邻波段，对该组进行标注。

步骤9:重复执行步骤4～步骤8，直到遍历完X的所有元素，输出所有相关特征的分组Y，通过计算MY(Y在M中的余集)得到特征波段集合Q。

步骤10:计算Yc这一组内部所有斯皮尔曼相关系数的平均值和相应的标准偏差，保留与待测组分相关性最大的特征波段作为这一组的代表性特征波段。

步骤11:重复执行步骤10，直至处理完Y的所有分组，每组代表波段组成集合P。

步骤12:通过计算集合Q∪P，输出最优光谱特征子集F。

采用算法1的特征选择方法，不仅能够剔除原始全光谱的所有相关特征值，而且大幅度地缩减了下一步特征精选(算法2)的筛选范围，能够达到提高算法反演精度和执行效率的双重目的。对于波长特征选择，遗传算法是一种非常有效的方法。适应度函数[17]的设计在遗传特征选择算法中起着至关重要的作用，但是常规的适应度函数与具体的解决问题领域联系不够紧密，导致最终筛选出的特征子集建立的模型预测精度较低。为此提出了一种新型的适应度函数，综合考虑了特征子集的回归性能和算法的运行效率，该算法的具体描述如算法2所示。

(2)算法2——遗传随机森林的特征选择算法

输入：①光谱特征集合F。②每一个体的染色体长度(基因)m集合F中元素的数量。③染色体的基因位数p1。④初始群体包含的个体数f。⑤最大迭代次数Generations200。⑥停滞代数StallGenLimit100。⑦交叉概率pn0.8。⑧变异算子pm0.5。

输出：最优光谱特征子集E。

步骤1(参数编码): 用一串含有mp个0/1字符(基因)的字符串(染色体串)来表示每种区间组合。若基因为1表示该特征被选中，基因为0表示该特征未被选中。

步骤2(种群初始化):随机产生fm个字符作为初始群体。

步骤3(新的适应度函数)：本算法中的适应度函数被定义为随机森林袋外误差和特征维度的加权和，即

(2)

式中α——RFR袋外误差，是衡量RFR模型性能的重要指标

β——特征子集含有的特征值数量，直接决定算法的运行效率

ω——回归精度在适应度函数中所占的比重

m——染色体长度(原始特征值数量)

步骤4(复制)：采用精英保留策略。设置精英个数c1为2，按照式(2)计算所有个体的适应度，选出适应度最低的前两名个体不进行配对交叉自动推到下一代。

步骤5(选择): 采用锦标赛选择策略。当精英子代从当前种群移除后，对剩余m-2个个体执行锦标赛选择算法。设置参加锦标赛的个体个数为2，依据式(2)评估所有个体的适应度，从种群中选出2个最优个体作为父代个体。

步骤6(交叉): 对锦标赛策略生成的父代个体进行交叉操作。选用的交叉函数是算术交叉类型的异或运算，对两个父类染色体进行异或操作产生交叉子代。产生交叉子代的数量为c2=round((m-c1)pn)。

步骤7(变异): 变异操作采用均匀变异的方式，最终产生变异子代的数量c3=m-c1-c2。

步骤8(生成子代个体): 精英子代，交叉子代和变异子代的所有个体组成了新一代种群。利用式(2)的适应度函数评价新一代种群的适应性。

步骤9:迭代步骤4～步骤8，直到满足停止条件。停止条件包含如下两类：①最大迭代次数达到设定值Generations。 ②停滞代数达到设定值StallGenLimit。

步骤10:输出最优光谱特征子集E。

传统的适应度函数一般只考虑回归模型中的交叉验证均方根误差和相关系数，提出的适应度函数综合考虑特征子集的回归性能和算法的运行效率。由式(2)可以得出，当随机森林的袋外误差取得最小值，并且选出的特征值数量最少时，本文算法设计的适应度函数取得最小值。

在变异操作算法中，首先生成符合某一范围内均匀分布的随机数集合(RD)，每个随机数的值与染色体中每个基因座上的原有基因相关联。然后从左到右扫描染色体，将每个RD的值与变异算子pm进行比较。如果位置i处的RD小于pm，则位置i处的基因被翻转；否则，该基因不被翻转。

2 实验

2.1 数据来源

在青岛市大沽河流域采集具有代表性的农田表层土壤样品124个，摊开并置于阴凉处，待自风干后，去除土样中的石子和动植物残体等明显杂物，用研钵研磨，再过100目筛子保存于自封袋中待测。准确称取0.250 0 g样品于150℃下以“HNO3-HF-H2O2”消解体系和密封高压釜消解罐法进行消解，利用电感耦合等离子体发射光谱仪(ICP-OES，Optima 8000，PerkinElmer，USA)来检测124个样品的镉含量，并保留每份土样剩余部分用于高光谱数据采集。

2.2 光谱采集及预处理

本次实验选取ASD FieldSpec 3型便携式光谱仪测量土壤样品的光谱反射率。将土壤样品放入直径10 cm、厚度约为2 cm的培养皿中。在50 W卤素灯作为光源的暗室中测量土壤样品的光谱，光源距样品35 cm，天顶角为30°，光谱仪探头安装在样品垂直上方15 cm处。每次实验开始之前和每测量10组样品之后使用具有100%反射率的标准化白板进行校正。每个土样测量时转动培养皿3次，每次转动90°，每个方向上取样10次，4个方向上共得到40条光谱测量值，算术平均后作为该样品的最终光谱反射率。

使用Kennard-Stone[18]算法从原始光谱数据集(包含124个样本)中选出100个样本作为校验集，剩余的24个样本作为预测集用于模型验证。剔除具有低信噪比的350～399 nm、2 451～2 500 nm两个边缘波段，在400～2 450 nm波段上对原始光谱数据采用多种预处理方法。基于随机森林回归模型，不同预处理方法所建模型的预测性能如表1所示。

由表1可得，一阶微分处理方法所建立的RFR模型结果最优，此时预测集相关系数、相对分析误差最大，分别为0.897 3和1.26，均方根误差达到最小值0.079 4。可见，一阶微分预处理最适于获取含量较少组分的光谱信息，与文献[19]研究结果一致。因此，后续特征选择与建模均在一阶微分分析基础上进行。

大沽河流域124个土壤样本的原始可见/近红外光谱曲线如图2a所示。从图2a可以看出，所有样本的反射率曲线特征基本相似。与文献[20]报道的一致，大多数土壤光谱在波长1 400 nm、1 800～2 000 nm和2 200～2 400 nm附近都会出现明显的水分吸收峰。图2b为最优预处理方式(一阶微分)后的光谱曲线。从图2b可以看出，一阶微分预处理不仅能够加强原始光谱的3个强吸收峰，而且能够明显增强原始光谱中1 000 nm和高于2 200 nm的弱吸收峰。这表明一阶微分预处理方法能够增强样品之间的光谱特征差异，较适应于农业土壤镉元素的高光谱响应。

表1 不同光谱预处理方法的随机森林建模结果Tab.1 Prediction results of RFR modeling by using different pre-processing methods

注：D1表示对光谱矩阵求一阶导，D2表示对光谱矩阵求二阶导，SNV(Standard normal variate transformation)表示对光谱矩阵进行标准正态变量变换，S-G(Savitzky-Golay)表示对光谱矩阵进行Savitzky-Golay卷积平滑。

3 结果与讨论

3.1 基于斯皮尔曼等级相关分析的特征选择算法结果分析

一般来说，特征值之间的相关性是通过相关系数来衡量的。由于斯皮尔曼等级相关系数对数据条件的要求没有皮尔逊相关系数严格，并且可以使用单调函数来描述变量之间的相关性，因此SGA-RF算法选用斯皮尔曼等级相关分析作为特征预选方法剔除冗余特征。对于初始的2 051个特征波段应用斯皮尔曼等级相关系数分析，在校验集上筛选的特征波长变量分布如图3所示。

图2 124个土壤反射光谱Fig.2 Soil reflectance spectra of 124 soil samples

图3 算法1筛选的波长变量分布图(不相邻波段之间的相关性用红色竖线表示)Fig.3 Distribution diagrams of spectral features selected by algorithm 1 (there was also a correlation between nonadjacent bands, which were shown on red long string)

算法1将原始全光谱变量从2 051个减少到108组特征波段,每组内部选出与镉浓度相关性最强的特征值为该组代表波段，所有代表波段在400～2 450 nm的分布如图3a所示。为了更加清晰地观察不相邻波段之间是否存在相关性，将图3a中蓝色框内的局部波长变量放大于图3b中，将图3a中紫色框内的局部波长变量放大于图3c中。其中，对特征值的标注已在算法1中步骤8完成。在图3c中以2 446 nm为代表波段的这一组被标识为黑色，表明该组内部所有成员均为相邻波段。在图3a中以999 nm为代表波段的这一组被标识为红色，表明该组内部包含不相邻波段。以特征波段999 nm与1 829 nm的相关性为例，算法1判定二者具有相关性；但是采用传统遗传算法的区间分割法，二者不具有相关性，原因是999 nm和1 829 nm肯定不会被划分到同一个子区间内部。因此，传统遗传算法按照固定区间大小均分原始光谱，会导致部分潜在重要特征值的丢失。

从图3a可以看出，预选出的108个特征波段主要位于1 550～1 850 nm和2 300～2 450 nm的范围内。分析图3b 可得，1 550～1 720 nm范围内共有9组特征波段内部含有不相邻波段。由图3c可得，1 800～2 450 nm范围内共有7组特征波段内部含有不相邻波段。108组相关特征中共有17组内部存在不相邻波段之间具有相关性的现象。可见，相邻波段之间存在高相关性的概率较大，但是不相邻波段之间也有可能存在相关性。因此，基于斯皮尔曼等级相关分析的特征选择方法克服了传统遗传算法只对相邻波段特征值进行融合的局限性，采用一种更加合理的方法缩减遗传算法的自变量个数，使得遗传算法能够在冗余度最小且最具代表性的敏感波段中进行全局搜索，最终达到优化反演模型，提高预测能力的目的。

3.2 系数对适应度函数的影响

适应度函数的表达式(式(2))中包含一个权重参数ω，取值范围为[0,1]，该参数表示袋外误差在整个适应度函数中所占比重。为了找到适应度函数取最小值时所对应的权重参数ω，需要考察不同权重对反演结果的影响。表2给出了不同权重参数ω所对应的适应度以及相对应的相关系数、均方根误差和预测相对分析误差。

从表2可以看出，当权重参数ω为0.3，适应度函数取得最小值。因此，最优权重参数ω为0.3。

3.3 特征选择结果分析

在算法2中，遗传算法种群规模被定义为100，运行总代数为200，基因位数设为1。

表2 不同权重参数ω所对应的适应度以及相对应的相关系数、均方根误差和预测相对分析误差Tab.2 Values of fitness function (fi) and corresponding R, RMSE and RPD according to several ω values

图4 SGA-RF算法特征选择过程适应度函数进化曲线Fig.4 Evolution curves of fitness function during feature selection process of SGA-RF

由图4可得，从第109代到第200代，适应度在0.014～0.016 5之间，变动幅度很小。最佳个体的适应度为0.015 397 7，平均值约为0.019 616 8。因此，可以判断SGA-RF算法没有出现过早收敛现象，原因是SGA-RF算法设置了正确有效的适应度，同时采用精英保留策略和锦标赛选择算法确保了该算法的必然收敛。

采用SGA-RF特征波段选择方法最终选出37个敏感波段，它们在400～2 450 nm波段的分布图如图5所示。图中黑色曲线为某个土壤样本的原始光谱曲线，红色空心圆圈表示特征波段所在位置。图5b和图5c分别放大了图5a中矩形框和正方形框内的敏感波段密集区域。

图5 基于SGA-RF算法的镉可见/近红外光谱特征波段分布图Fig.5 Visible/NIR spectra distribution diagrams for soil Cr based on variables selected by SGA-RF (shown in red hollow circle markers)

由图5可知，重金属镉的敏感波段主要位于1 600～1 800 nm和2 350～2 450 nm的范围内。已有研究表明，位于这2个范围的吸收峰主要与有机化合物中的O—H、 N—H和C—H基团有关[21]。这可能与重金属镉的有机化合物吸附有关，吸附机理有待进一步深入研究。

3.4 与其他特征选择算法的比较与分析

为了验证算法1，分别对不采用算法1的遗传随机森林算法(GA-RF)和采用算法1做预选的SGA-RF算法进行建模性能比较，结果如表3所示。两种算法均采用RFR建立回归模型。

表3 评估算法1对本文算法预测效果的影响Tab.3 Impact of algorithm 1 on predictive effect by using SGA-RF

由表3可以看出，与GA-RF算法相比，SGA-RF算法不仅建模时间短，特征值数量少，而且模型性能明显较高。因为算法1在有效剔除初始特征集冗余信息的同时能够保留与镉元素相关性最强的重要特征，大幅度减少特征精选阶段的输入特征值维数，不仅减少了计算复杂度，而且提高了模型反演能力。

为了进一步验证SGA-RF特征选择算法的有效性，本文算法与其他3类特征选择方法进行对比，结果如表4所示。这3类特征选择方法分别为：全光谱不做特征波段选择方法、单一特征选择方法和以K最近邻法构建适应度函数的遗传算法 (SGA-KNN)。单一特征选择方法包含两种方法：原始遗传算法与偏最小二乘算法结合的特征波段选择方法(GA-PLSR)和斯皮尔曼等级相关分析(Spearman Rank)方法。这3类特征选择算法均采用RFR作为回归模型，并以选出特征波段的数量、建模时间、相关系数、均方根误差和预测相对分析误差作为模型的评价指标。

依据文献[21]，GA-PLSR算法采用2～20 nm不同波段间隔对原始光谱进行重采样，每个子区间取平均数作为新的光谱自变量，以交叉验证均方根误差(RMSECV)最小值作为最优特征谱区筛选的标准,最终得到12 nm为最优光谱区间划分间隔。GA-PLSR算法使用PLS-Toolbox[22]重复运行50次，以PLSR算法的RMSECV作为适应度函数。SGA-KNN算法采用文献[23]中的基于K最近邻分类法(KNN)的适应度函数，其中参数K选取3。上述2种算法的种群规模、运行总代数和SGA-RF算法保持一致。

表4 本文特征选择算法与其他特征选择算法性能比较Tab.4 Comparison of performance between proposed and other feature selection methods

由表4可知，与其他3类特征选择算法相比，本文算法的建模时间最短，相关系数和预测相对分析误差最高，取得最佳预测效果。全光谱和斯皮尔曼等级相关分析建模的预测相对分析误差均在1.4以下，说明这2种方法建立土壤镉含量预测模型的性能较差，原因是全光谱中含有大量的冗余和无关信息；斯皮尔曼等级相关分析只能去除特征波段之间的冗余，无法剔除那些与待测组分无关的特征值，导致模型精度较低。GA-PLSR算法和SGA-KNN算法建模的预测相对分析误差均在1.0以下，说明以GA-PLSR算法和SGA-KNN算法作为特征波段选择方法无法建立土壤镉含量预测模型。原因为：① GA-PLSR算法不仅忽略了不相邻波段之间的相关性，同时也没有考虑特征波段与待测组分之间的相关性，而且该算法用于评估每个染色体的PLSR模型为线性模型。②由于KNN算法对参数K的取值，距离的度量和多重共线性都非常敏感，导致该算法对特征的回归能力较差，并不适用于构造遗传算法的适应度函数。③SGA-RF算法首先采用斯皮尔曼等级相关分析进行特征预选，去掉所有特征波段之间的冗余信息；其次，结合遗传算法的全局搜索能力和RFR的较高反演能力，优选出适应度最低的特征子集。其中，在适应性函数中采用非线性模型RFR评估每个个体的性能，而RFR对于参数设置与多重共线性问题均不敏感。

在SGA-RF算法特征筛选结果的基础上，以RFR作为回归模型，模型校验集和预测集样本的估测值和实测值之间的关系如图6所示。

图6 SGA-RF算法结合RFR模型建模样本、检验样本的实测值与估测值比较Fig.6 Comparison of measured Cd contents and estimated values of modeling and testing samples through SGA-RF algorithm with RFR model

由图6可得，采用SGA-RF特征选择算法在校验集和预测集均取得了较好的反演效果，验证了该算法与随机森林回归模型相结合能够有效解决农业土壤中痕量级重金属含量的反演问题。

4 结论

(1)基于大沽河流域农业土壤样本真实数据集，SGA-RF算法能够将2 051个原始特征波段优选至37个。以选出的37个敏感波段作为自变量建立大沽河流域农业土壤镉高光谱反演模型，该模型具有较低的预测均方根误差(0.060 1)，较高的相关系数(0.950 2)和预测相对分析误差(2.03)。建模结果证实了该算法在降低计算复杂度的基础上提高了土壤镉含量反演模型的精度。

(2)斯皮尔曼等级相关分析将原始2 051波段缩减到108组特征子集，通过分析每组特征子集的组内相关性，发现在不相邻波段之间也存在相关性。说明基于斯皮尔曼等级的特征预选方法能够克服传统遗传算法的局限性，最大限度地缩减遗传算法的自变量个数，是SGA-RF算法的重要组成部分。

(3)以相关系数、均方根误差和预测分析相对误差作为评价标准，本文算法分别与其他4类特征选择方法进行了对比。比较结果说明，SGA-RF算法可以有效剔除冗余波长变量，提取最具代表性的特征波长变量，减少建模时间，提高预测效果。

(4)SGA-RF算法有效解决了土壤重金属镉高光谱反演领域建模时间长，预测精度差的问题，可用于其他农业土壤的重金属镉污染监测与识别。