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

图1

Fig.1

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

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

表1

Table 1

表1(Table 1)

表1 土壤样本养分元素含量基本信息 Table 1 Statistics of nutrient contents in soil samples 养分元素 样本数 最小值/ (g· kg-1) 最大值/ (g· kg-1) 均值/ (g· kg-1) 标准差/ (g· kg-1) 有机质 80 2.138 5.517 3.364 0.535 磷 80 0.452 0.988 0.682 0.087 钾 80 1.875 2.240 2.036 0.062

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

在暗室中进行土壤光谱测量, 测量仪选用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。 在去噪和重采样处理的基础上, 将原始光谱反射率进行一阶微分、 倒数对数等特征变换, 不同的变换形式可以帮助准确、 快速地找到峰谷, 并通过峰谷确定相应的波长, 从而确定敏感波段。

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

ft=fmin+(fmax-fmin)βvit+1=vit+(zit-z* )fizit+1=zit+vit+1(1)

其中, vit和 vit+1分别表示蝙蝠个体i在t和t+1时刻的飞行速度; zit和 zit+1分别表示蝙蝠个体i在t和t+1时刻的位置; z^* 代表全局最优位置。 频率f_i为蝙蝠个体i在搜索时的脉冲频率, 其中, β 为介于[0, 1]间的随机数, (f_min, f_max)为脉冲频率范围。 每次迭代过后, 音强和频度会根据脉冲响度衰减系数和脉冲频度增加系数进行更新计算。

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

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

R2=1-∑i=1n(yi-Yi)2∑i=1n(yi-y¯)2, RMSE=1n∑i=1n(yi-Yi)2(2)

其中n为样本量, y_i为实测值, Y_i为预测值, y¯为实测值的平均值。 预测模型计算的RMSE越低, R²越接近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

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

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

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图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

表2(Table 2)

表2 最大相关系数和敏感波段 Table 2 Maximum correlation coefficients and sensitive bands 元素 光谱变化形式 R R' lg1/R (lg1/R)' 最大值 -0.678 0.796 0.656 -0.778 有机质 对应波段 600 1 370 600 1 140 敏感波段数量 95 85 82 101 最大值 -0.452 0.637 0.458 -0.688 磷元素 对应波段 440 1 390 440 690 敏感波段数量 11 38 13 82 最大值 0.428 -0.436 -0.429 0.553 钾元素 对应波段 450 1 650 700~710 980 敏感波段数量 23 5 26 51

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

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

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

图4

Fig.4

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

为了分析BA-Adaboost模型的预测精度, 分别计算了BA优化前后模型训练集和验证集的R²和RMSE(表3)。 从表中可以看出无论是训练集还是验证集有机质的高光谱预测精度均高于磷元素和钾元素, 经过BA寻优后, Adaboost模型的精度有了较大的提升, 从验证集结果来看, 磷元素的 Rp2变化最为显著, 由0.695提高至0.808; 有机质的RMSEP变化最大, 由0.204 g· kg^-1下降至0.152 g· kg^-1。

表3

Table 3

表3(Table 3)

表3 模型预测精度对比分析 Table 3 Comparison of estimation accuracies 元素 光谱变化形式 预测模型 训练集 验证集 Rc2 RMSEC Rp2 RMSEP 有机质 一阶微分 Adaboost 0.925 0.140 0.756 0.204 BA-Adaboost 0.961 0.101 0.864 0.152 磷元素 倒数一阶微分 Adaboost 0.909 0.025 0.695 0.039 BA-Adaboost 0.944 0.019 0.808 0.031 钾元素 倒数一阶微分 Adaboost 0.877 0.019 0.651 0.029 BA-Adaboost 0.917 0.016 0.761 0.024

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

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

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图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