BA-Adaboost模型的黑土区土壤养分含量高光谱估测

林 楠，刘海琪，杨佳佳，吴梦红，刘翰霖

1. 吉林建筑大学测绘与勘查工程学院，吉林 长春 130118 2. 吉林大学地球科学学院，吉林 长春 130026 3. 东北大学资源与土木工程学院，辽宁 沈阳 110819 4. 中国地质调查局沈阳地质调查中心，辽宁 沈阳 110034

引 言

土壤养分含量快速估测对作物生长监测和黑土科学管理具有重要的意义。 高光谱数据可以通过数千个紧密排列的波长记录土壤信息。 众多研究表明，土壤高光谱数据与土壤含水率[1]、养分含量[2-3]、重金属含量[4]等土壤指标存在密切关系，利用高光谱技术进行土壤元素含量估测已取得了丰硕的研究成果。

目前，利用机器学习算法建立土壤高光谱特征和元素含量之间的反演模型是土壤元素含量估测的主要方法，其中支持向量机[5]、极限学习机[6]、随机森林[7]等模型都得到了广泛的应用。 机器学习算法通常需要预先设定一些模型参数值，由于预定义的参数值很可能不包含全局最优的参数值，导致机器学习模型达不到最佳效果。 为了克服机器学习模型在寻找最佳模型参数时所存在的问题，遗传算法、粒子群算法等经典优化算法被用来优化机器学习模型的内部参数。 然而粒子群等经典优化算法对初始参数设置比较敏感，寻优过程中容易陷入局部最优解，导致算法后期的收敛速度变慢。 蝙蝠算法是一种新兴的群体智能方法，在寻优过程中，蝙蝠算法模仿蝙蝠声波脉冲响度、频率的自适应调节过程，实现全局寻优和局部寻优过程的自由切换，从而使算法的全局搜索能力和局部搜索能力得到平衡，在模型参数寻优计算方面表现优异[8]。 本工作对黑土土壤有机质、磷元素和钾元素的光谱特征进行详细分析，选择元素最佳光谱变换形式及相关性较高的敏感波段，将蝙蝠算法和Adaboost机器学习模型组合构建土壤养分含量估测模型，利用蝙蝠算法求解Adaboost模型建模中的关键参数，并比较模型参数优化前后的估测精度，为黑土养分含量高光谱估测提供一种高效的新方法。

1 实验部分

1.1 土壤样品采集

土壤样品采集区位于黑龙江省讷河市，该区地处松辽平原的北端，大小兴安岭南缘，是黑土分布的典型区域。 区内耕地面积广阔，农产品种类丰富，是甜菜和马铃薯的盛产地。 2019年4月在该区内挑选典型黑土地块，设置5公里大小的网格作为采样单元，采样点设计结合第二次全国土地调查成果合理分布，所布设的采样点要能够代表采样单元的土壤性质，试验区内共采集土壤样本80个(图1)，采集样品过程中确保采样距离公路至少150 m，以每个采样点位置为中心，在其周围5 m×5 m范围内进行样品采集，每个采样点位置共采集5份黑土样本，取样深度为表层土的15 cm以内，将样本充分混合后装入采样袋。

图1 研究区和样品采集点地理位置Fig.1 Location of study area and sampling collection points

1.2 元素含量测量及样本划分

将样品中的秸秆、砂砾等杂物剔除，风干后研磨过筛，使得土壤的粒径小于0.25 mm; 将样品分成两份，一份用于测定土壤元素含量，另一份用于室内高光谱测量。 结合多目标样品分析质量要求等技术规范，参照《土壤农业化学分析方法》，选用重铬酸钾容量法测定土壤有机质含量，选用X射线荧光光谱法进行土壤磷元素和钾元素含量测量[9]。 对元素含量的测量结果进行统计分析(表1)，将80份土壤样品按元素含量由低到高分为20组，每组随机抽取1个样品放入验证集，共20个样品作为验证集，其余60个样品为训练集。

表1 土壤样本养分元素含量基本信息Table 1 Statistics of nutrient contents in soil samples

1.3 光谱采集及预处理

在暗室中进行土壤光谱测量，测量仪选用ASD FiledSpec4地物波谱仪，为提高光谱测量数据的精度，取5次光谱测量的均值作为土壤样品的反射光谱数据。 由于光谱数据在350～399 nm以及2 451～2 500 nm范围的噪声较大，信噪比低，对土壤元素与反射率关系的分析会造成干扰，所以将400 nm前及2 450 nm后的数据剔除。 光谱仪采样间隔为1 nm，即在400～2 450 nm范围内得到2 051个波段，由于光谱分辨率较高，波段数较多，相邻波段之间可能存在信息重叠，更容易受到噪声的影响，因此，将光谱数据进行重采样处理，设置采样间隔为10 nm。 在去噪和重采样处理的基础上，将原始光谱反射率进行一阶微分、倒数对数等特征变换，不同的变换形式可以帮助准确、快速地找到峰谷，并通过峰谷确定相应的波长，从而确定敏感波段。

1.4 模型概述

蝙蝠算法(Bat algorithm，BA)是一种启发式搜索算法，它模拟蝙蝠使用声纳探测猎物和躲避障碍物，通过模拟蝙蝠飞行寻找猎物的过程来模拟优化搜索过程，在计算过程中利用求解问题的适应度值来选取蝙蝠的位置，利用优胜劣汰的进化过程来模拟较优可行解代替较差可行解的迭代搜索过程[10]。 基于BA的基本原理，在算法各项参数被初始化后，从d维搜索空间中的一个随机位置zl开始启发式搜索。 以固定的频率、不同的波长和音强搜索猎物，搜索过程中，蝙蝠根据接近猎物的距离自动调整波长的大小。 经全局搜索后更新每只蝙蝠的飞行速度和空间位置，并计算目标函数的适应度值，速度和空间位置更新公式如式(1)[11]

ft=fmin+(fmax-fmin)β

(1)

Adaboost算法基本思想是针对同一训练集训练能力一般的弱回归器，通过叠加方法构建一个训练能力很强的强回归器。 其算法本身是通过改变样本分布权重来实现的，根据每次训练集中每个样本的预测精度以及上次的总体预测精度来计算每个弱回归器的权重，同时更新每个样本的分布权重，最后将每次训练得到的回归器结果加权求和，作为强回归器最后的输出结果[12]。 建模过程中，弱学习器的最大迭代次数n和弱学习器的权重缩减系数v是两个重要的参数，迭代次数n设置过小会导致模型拟合不充分，n设置太大则会导致模型拟合过度，而较小的权重缩减系数v意味着需要更多的弱学习器的选代次数，所以通常这两个参数要一起进行优化调整[13]。

1.5 模型精度检验

为了评价模型预测能力及模型的稳定性，选择决定系数(R2)和均方根误差(RMSE)对模型的建模效果进行评价

(2)

2 结果与讨论

2.1 光谱特征变换

经处理后的室内土壤样品光谱曲线如图2(a)所示，从图中可以看出，采集的黑土样品原始反射率在0～0.9之间，各样品光谱曲线波动形状相似，在可见光波段，随着波长的增加，反射率逐渐增加，直到1 200 nm处反射率趋于稳定。 近红外区土壤光谱反射率总体高于可见光区，两个明显的波谷分布在1 400和1 900 nm附近，主要是土壤中残留的水分和空气中水蒸气的吸收造成的，而在2 200 nm处有轻微的凹陷，是受土壤中存在的粘土矿物的影响。 图2(b)—(d)分别为原始反射率经一阶微分(R′)、倒数对数(lg1/R)、倒数对数一阶微分[(lg1/R)′]变换后的光谱曲线，从变化结果可以看出，一阶微分变换可以对原始光谱变化起到放大作用，经变换后反射率在1 400，1 900和2 200 nm处的波动更大。

图2 反射率及其变换的土壤光谱曲线Fig.2 Soil spectral reflectance curves and its transformations

2.2 相关性分析及特征波段选择

分别计算了土壤有机质、磷、钾含量与土壤反射率的相关系数，并绘制相关系数曲线图(图3)。 从图3(a)原始光谱的相关系数可以看出，土壤有机质和磷元素含量与光谱反射率呈负相关，而钾元素含量则相反。 与原始光谱反射率相比，变换后的光谱数据与土壤有机质、磷、钾含量的相关性更高，其中，一阶微分变换形式与土壤有机质、磷、钾含量的相关系数表现为正、负交叉，波峰和波谷较多，而且经一阶微分变换后，各元素的最高相关系数显著提高。

图3 土壤元素含量与光谱相关系数曲线Fig.3 Correlation coefficient curves of soil element content and spectral reflectance

选取相关系数大于0.4的敏感波段作为预测模型的样本输入数据(表2)。 由统计结果可以看出，经过不同的光谱特征变换，土壤元素含量与光谱反射率的相关系数有所提高，其中，土壤有机质对应的最佳变换形式为一阶微分，与光谱反射率数据的相关系数最高为0.796，波长范围在1 370 mm附近； 土壤磷元素对应的最佳变换形式为倒数对数的一阶微分，相关系数最高为-0.688，波长范围在690 mm附近； 土壤钾元素对应的最佳变换形式为倒数对数一阶微分，相关系数最高为0.553，波长范围在980 mm附近。

表2 最大相关系数和敏感波段Table 2 Maximum correlation coefficients and sensitive bands

2.3 基于BA-AdaBoost模型的预测

将选取60个训练样本最佳光谱变换形式的敏感波段和对应元素含量值作为Adaboost模型的建模数据，3种元素的建模波段数分别为85个、82个和51个，选择CART决策树为Adaboost模型的弱回归学习器，采用蝙蝠搜索算法对Adaboost模型的最大迭代次数n和弱学习器权重缩减系数v进行寻优。

基于BA的基本原理，首先需对模型各项参数进行初始化，BA初始化参数较多，但除了迭代次数外其他参数敏感性不强，可选择默认缺省参数进行初始化。 BA算法的搜索空间是由最大迭代次数n和弱学习器权重缩减系数v为坐标轴组成的二维空间，迭代搜索过程从搜索空间内的L个随机位置开始，在每次迭代过程中，利用每只蝙蝠在空间中位置的二维坐标(n，v)作为Adaboost模型的初始化参数，然后基于Adaboost模型对样本数据训练建模，并计算预测结果。 根据计算结果选取R2值最大时对应的位置作为蝙蝠当前最优位置，利用式(1)更新每只蝙蝠的空间位置。 设置种群大小L=20，脉冲频率范围fmin=0、fmax=1、脉冲音强范围Amin=0，Amax=1，脉冲响度衰减系数α=0.9，脉冲频度增加系数γ=0.9，作为缺省参数，并计算了不同迭代次数对应的R2值(图4)。 从图中可以看出，随着迭代次数T的增加，R2值逐渐增大，3种元素中，钾元素收敛速度最快，当迭代数达到10次的时候，R2值达到了最大，有机质和磷元素分别是当迭代数达到12次和16次达到最大。 建模估测时，选取R2最大时对应的n和v的值作为Adaboost模型的建模参数。

图4 BA-Adaboost模型R2变化曲线Fig.4 R2 variation curves of BA-Adaboost

2.4 模型精度分析

表3 模型预测精度对比分析Table 3 Comparison of estimation accuracies

利用模型的预测值和实测值绘制预测结果拟合图，进一步对比分析3种元素模型拟合效果(图5)。 从图中可以看出，与Adaboost相比，优化后的BA-Adaboost模型的预测值更紧密地分布在1∶1左右，数据拟合能力和稳定性优于Adaboost模型。 同时可以看出训练集的拟合效果比验证集更优异，验证集样本点分布相对分散，说明机器学习模型对训练样本数据的可靠性要求较高。

3 结 论

以黑龙江省讷河市80个黑土样品为研究对象，对黑土土壤有机质、磷元素和钾元素的光谱特征进行分析，构建BA-Adaboost模型对元素含量进行估测，主要结论如下：

(1)光谱数据的特征变换处理可以增强光谱特征，一阶微分变化有效地突出了光谱曲线的波峰和波谷，提高了光谱反射率与元素含量的相关系数，有机质、磷元素、钾元素分别在1 370，690和980 mm处相关系数达到最大值。

图5 元素含量实测值和预测值(a): OM含量Adaboost模型； (b): OM含量BA-Adaboost模型； (c): P含量Adaboost模型； (d): P含量BA-Adaboost模型； (e): K含量Adaboost模型； (f): K含量BA-Adaboost模型Fig.5 Measured and predicted values of element content(a): Adaboost of OM content; (b): BA-Adaboost of OM content; (c): Adaboost of P content; (d): BA-Adaboost of P content; (e): Adaboost of K content; (f): BA-Adaboost of K content

(2)蝙蝠算法在参数寻优过程中能够动态控制全局寻优和局部寻优过程之间的自由切换，利用蝙蝠算法对Adaboost模型参数进行优化计算，避免了估测模型陷入局部极小值，与传统参数寻优方法相比，该方法可更快更准确地找到全局最优参数。

(3)将BA与Adaboost模型相结合，构建BA-Adaboost土壤含量估测模型，该组合模型只需设置搜索空间，然后自动搜索模型最优参数值； 对比BA算法优化前后的估测精度可以看出，优化后BA-Adaboost模型R2增大、RMSE变小，估测精度明显提高，说明BA-Adaboost模型在土壤元素含量高光谱估测中具有一定的适用性，扩展了机器学习模型在土壤成分估测研究中的应用。