土壤氧化铁的特征波长选择和高光谱反演*

赵海龙, 甘 淑, 2**, 王俊杰, 胡 琳

(1. 昆明理工大学国土资源工程学院 昆明 650093; 2. 云南省高校高原山区空间信息测绘技术应用工程研究中心 昆明650093)

铁是土壤矿物中主要的元素之一, 是最重要的过渡元素, 而土壤中铁大多是以化合态的形式存在。由于氧化铁的聚集和迁移活动受到化学风化、生物循环等因素影响, 反映了土壤的淋溶过程、风化发育程度以及土壤的地带性分布特征, 常被作为描述土壤发育和土壤分类的指标之一。因此, 掌握土壤中氧化铁含量信息对于揭示土壤环境状况, 指导农作物生产具有重要意义。

常规的氧化铁含量测定技术成本高、周期长, 利用高光谱技术可以快速而高效测定土壤氧化铁含量。高光谱以其光谱分辨率高、波长连续性强的特点, 被广泛用来进行土壤理化属性的反演。目前已有许多学者对氧化铁和反射光谱之间的关系做出了许多成果。何挺等指出土壤中氧化铁含量与反射率呈负相关, 且在400 nm、450 nm、490 nm、700 nm、870 nm等处存在吸收峰。彭杰等指出氧化铁含量与斜率呈线性负相关, 与截距呈线性正相关, 由此说明了利用土壤线的参数预测氧化铁含量具有可行性。熊俊峰等利用偏最小二乘回归法将全反射波段分别与土壤中的全铁、游离铁、无定形铁的含量进行回归建模, 结果表明游离铁和无定形铁精度变化不大, 全铁的估算精度变化较大, 因为全铁中的硅酸铁易受到有机质和土壤深度的干扰。谭洁等指出不同森林土壤氧化铁含量和不同土壤类型光谱曲线在可见光-近红外波段内变化趋势基本一致, 均形似陡坎。基于特征波长筛选方面的研究, 丁海宁等利用相关系数峰值-逐步分析方法筛选出氧化铁的特征波长, 构建的估算模型均较稳定。李双权等利用相关性达<0.05显著性水平的波长构建氧化铁的偏最小二乘回归估算模型, 模型的决定系数达0.791。阳洋等提取的敏感波长是把通过<0.01显著性检验的波长进一步进行多元逐步线性手段剔除共线性波长获取的, 最终构建的模型决定系数达0.790。朱亚星等利用无信息变量消除(UVE)、竞争性自适应重采样(CARS)、无信息变量消除结合竞争性自适应重采样(UVE-CARS)、竞争性自适应重采样结合无信息变量消除(CARS-UVE)这4种变量筛选方法, 得到有机质最佳的变量筛选方法为UVE-CARS, 基于变量筛选后的模型性能优于全波段模型。

综上, 当前学者们对氧化铁反演研究已经取得了较好的成果, 土壤氧化铁的吸收特征都已明确, 但是土壤的理化性质比较复杂, 不同地区其特征波长有所差异, 导致在定量研究氧化铁时, 只根据吸收特征波段来反演氧化铁含量的精度不够理想。同时, 由于光谱数据的高度冗余性, 仅利用相关性分析不能达到去除冗余性的效果。云南禄丰恐龙谷南缘地区分布有多种土壤类型, 其土地利用类型多样, 同时有着特殊的地形构造。因此, 本文以云南禄丰恐龙谷南缘地表土壤氧化铁为研究对象, 在基于相关性分析的基础上, 利用IRIV、CARS和SPA等波长筛选算法对光谱数据进一步的降维, 从而确定该地区土壤氧化铁的特征波长, 并利用随机森林和偏最小二乘回归方法构建土壤氧化铁的反演模型, 用以实现该地区氧化铁含量快速、准确地预测。

1 研究区概况

恐龙谷南缘地区位于云南省中部滇中高原的楚雄彝族自治州禄丰市境内(24°55′25″～25°22′05″N, 102°00′00″～102°9′00″E)。研究区东西宽约6 km、南北长约8 km, 面积约为32 km, 整体呈现四周高, 中间低的凹陷坑, 最高海拔为2200 m, 最低海拔为1302 m。研究区所在的禄丰市属亚热带低纬度高原季风气候。该地区是一个小型的中生代红色沉积盆地, 属于奥陶系下统红石页岩组, 简要的岩性为紫红色、灰绿色粉砂岩。据禄丰市1982年及1985年两次土壤调查资料显示, 市内有棕壤、黄棕壤、红壤、紫色土和水稻土5个土类、10个亚类、20个土属、40个土种。其中紫色土占土地面积的56.9%, 为辖区内最主要的土壤类型; 其次为红壤, 占土地面积的22.8%; 黄棕壤占7.8%; 水稻土占6.3%。

2 研究方法

2.1 样本采集与数据获取

土壤样品于2021年7月底采于云南省楚雄彝族自治州禄丰市恐龙谷南缘山地, 根据地势地貌差异进行采样点布设。每个样本点采样范围为5 m×5 m。在采样范围内, 根据5点采样法, 采集5～20 cm的表层土壤, 取约1 kg土壤装袋保存, 并用手持GPS记录该点的经纬度坐标, 共135份样品, 采样点分布如图1所示。采集的土壤类型有紫壤、红壤和黄棕壤。将采集回来的土壤首先去除杂草、石块等杂质, 然后进行自然风干处理, 最后用玛瑙球碾磨机研磨并过100目筛。将处理后的每份样品分成两份, 一份用于测定土壤氧化铁的含量, 另一份进行土壤的高光谱数据测量。结合样品分析质量要求等技术规范, 在满足样品的检出限、准确度、精密度的前提下, 参照《土壤农业化学分析方法》, 选用X射线荧光光谱法进行土壤氧化铁含量的测定。

图1 土壤采样点分布图Fig.1 Distribution map of soil sampling points

土壤光谱测量在暗室中进行, 光谱仪为ASD Field Spec 3 地物光谱分析仪, 光源为探头内置卤素光源, 探头内径为21 mm, 前视场角为25°, 波段范围为350～2500 nm。在300～1000 nm间隔内, 光谱采样间隔和分辨率分别为1.4 nm和3 nm; 在1000～2500 nm内, 光谱采样间隔和分辨率分别为2 nm和10 nm。将光谱采样间隔重采样为1 nm后, 得到的波长数为2151个。光谱测量时, 把土壤样品放在宽10 cm、高2 cm的容器内, 并把土壤样品刮平, 以减少土壤样品粗糙度对光谱测量的影响。探头与土壤样品保持2 cm的高度, 并垂直对准样品。每个样品采集5条光谱曲线, 对5条光谱曲线取平均后作为样本的实际光谱反射率曲线。

2.2 数据预处理

受仪器自身的影响, 300～399 nm、2451～2500 nm处信噪比较低, 因此把该两段光谱数据作去除处理, 共得到2051个波长。为了消除由于周围环境引起的干扰误差, 去除随机因素的干扰, 采用窗口数为9, 多项式阶数为2的Savizky-Golay平滑之后的曲线作为原始光谱曲线(origin spectral reflectance, OR)。为了进一步降低环境对光谱数据的干扰, 加强土壤氧化铁含量与光谱反射率之间的相关性, 对OR进行一阶微分(first-order differential reflectance, FD)、倒数的对数(reciprocal logarithm reflectance, RL)变换处理。

2.3 特征波长选取

由于高光谱数据冗余度较大, 在进行回归分析时, 并不是所有的波长都对建模精度的提高有益, 如果把所有的波长都进行建模分析, 不仅运算量巨大, 而且会降低建模精度。因此在建模之前进行特征波长选取是必要的。

利用相关系数法(correlation coefficient method, CC)对氧化铁含量与原始光谱及其各变换光谱进行皮尔森相关性分析, 计算每个波长与土壤样本氧化铁含量的相关系数, 把通过=0.01显著性检验的波长作为特征波长。相关系数的公式如下:

迭代保留信息变量(iteratively retaining informative variables, IRIV)是一种通过随机组合来考虑变量之间可能存在的交互作用的策略, 基于模型集群分析方法, 逐个波长变量计算包含或者不包含该变量时的均方根误差平均值, 得到两者之差(DMEAN)和非参数检验方法曼-惠特尼检验的值, 从而确定该变量的重要性。每次迭代都保留强信息变量和弱信息变量, 直至消除无信息变量和干扰变量, 最后逆向消除获得最佳的特征变量。通过相关性分析获得的波长中, 并不是所有波长都对建模效果有益, 有些波长信息量可能对建模效果有干扰。因此, 进一步利用IRIV处理可以获得对建模效果有益的强信息变量和弱信息变量。

竞争性自适应重加权算法(competitive adaptive reweighted sampling, CARS)是一种基于蒙特卡洛采样和PLS回归系数的特征波长选择方法, 把每个变量当成一个个体, 选择适应性能力较强的个体。具体步骤为: 随机选取固定比率的样本作为校正集并建立PLS模型, 然后计算模型的回归系数的绝对值和每个波长对应的权重, 利用指数衰减函数和自适应重加权采样法对变量进行选择, 同时计算交叉验证均方根误差,次采样后, 选择均方根误差最小的子集作为最优变量子集。本研究中, 对通过相关性分析获得的波长中, 进一步利用CARS可以消除无信息变量, 最终选择对模型起关键作用的变量。

连续投影法(successive projections algorithm, SPA)是一种前向变量选择算法, 运用变量投影操作寻找冗余度最低、共线性最小的光谱特征变量, 从而提高模型的稳定性和准确性。本研究中, 通过相关性分析获得的波长之间具有共线性, 利用连续投影法, 可以消除波长之间的共线性。

2.4 模型建立与评价

随机森林回归(random forest regression, RF)是一种基于决策树的集成算法, 对于数据是非线性以及较大时应用较多。它本身不是一个单独的机器学习算法, 而是通过在数据上构建多个模型, 集成所有模型的建模结果, 以此来获取比单个模型更好的结果。随机森林模型构建在python第三方库scikit-learn中实现。

偏最小二乘回归(partial least squares regression, PLSR)集合了多元线性回归分析、相关性分析和主成分分析的优点, 可以很好地解决变量之间的多重共线性问题, 同时能够揭示光谱数据中最大氧化铁含量变化的主成分。在模型构建时采用交叉检验可以确定最佳的主成分因数。偏最小二乘回归利用Minitab分析软件实现。

采用Kennard-Stone (K-S)算法进行建模集和验证集的划分, 选用70%的样本划分为建模集, 剩下的30%为验证集。根据CC和基于CC的CARS、IRIV、SPA算法筛选出的波长分别作为自变量, 土壤氧化铁含量作为因变量, 建立反演模型。由于土壤样本呈现非正态分布, 采用RPIQ可以给出模型更真实的评价。因此反演模型的精度指标采用决定系数(R)、均方根误差(RMSE)、性能与四分位数间距比(RPIQ)和1∶1线4个参数来衡量。

3 结果与分析

3.1 氧化铁含量统计分析

利用K-S法选出来的建模集共有95个样本, 验证集共有40个样本。利用Origin 2021经过统计分析可以看出研究区土壤氧化铁最小值为18.293 g∙kg, 最大值为66.978 g∙kg, 均值为41.201 g∙kg, 变异系数为28.393%, 属于中等变异。氧化铁统计特征如表1所示。氧化铁含量与诸多因素有关, 地类、成土母岩以及地质地貌都会影响土壤氧化铁含量, 而其中影响最大的是成土母岩。该地区主要的岩性为紫红色、灰绿色粉砂岩。

表1 研究区氧化铁含量统计特征Table1 Statistical characteristics of iron oxide content of the study area

3.2 土壤光谱特征

图2为研究区全部土壤样本的反射率光谱曲线, 其反射率整体较低, 范围为0.02～0.49, 但其光谱曲线变化趋势大致相同。

图2 研究区全部土壤样品的光谱反射率曲线Fig.2 Soil spectral reflectance curves of all soil samples of the study area

在可见光波段, 光谱反射率较低, 但增加较快;在近红外波段, 光谱反射率较高, 光谱曲线变化较为平缓; 在短波红外波段, 光谱反射率略有下降。同时, 在波长500 nm、900 nm处有微弱氧化铁吸收峰; 在波长1400 nm、1900 nm和2210 nm处分别是OH和HO两者的合成峰、HO吸收峰和黏土矿物中金属-OH振动合频产生的吸收峰。为了探明不同氧化铁含量的光谱特征, 将样本集的氧化铁含量先从高到低排序, 并分成6份, 计算每一部分的平均氧化铁含量和平均光谱反射率曲线, 得到如图3所示的不同氧化铁含量的光谱反射率曲线图。可以发现土壤的氧化铁含量越高，相应的光谱反射率越低。

图3 不同氧化铁含量的土壤样本的光谱反射率曲线Fig.3 Spectral reflectance curves of soil samples with different iron oxide contents

3.3 特征波长选取

如果将全波段直接利用变量优选算法去筛选波长, 不仅效率太低, 且可能会丢失对建模有用的信息。因此, 本文先利用相关性分析进行变量的粗选, 然后利用IRIV、CARS和SPA算法进行变量的精选。具体步骤为: 首先将训练集样本的原始光谱(OR)、OR一阶微分(FD)以及OR的倒数对数(RL)与氧化铁含量进行相关性分析, 把通过<0.01显著性检验的波长作为粗选的特征波长, 然后利用IRIV、CARS和SPA对光谱数据进一步的降维。

OR、FD以及RL与氧化铁含量相关性分析后, 通过<0.01显著性检验的波长分别达2051个、1252个和2051个, 即原始光谱和倒数的对数光谱全波长都通过了<0.01显著性检验(图4)。

图4 原始光谱及其变换光谱与土壤氧化铁含量的相关系数Fig.4 Correlation coefficient of original spectrum and its transformed spectrum with iron oxide content

经过相关性筛选后的波长仍然较多, 波长之间的共线性强, 因此需要利用变量优选算法进一步减少波长的数量。特征波长选择的结果如图5所示。

图5 不同算法选择的特征波长Fig.5 Feature bands selected by different algorithms

以FD为例, 本文设置IRIV波长筛选算法为5折交叉验证, 最大主成分为10。随着迭代次数的增加, 保留的变量逐渐减少, 之后进行反向消除从而保留52个变量, 迭代次数为8。选择的波长在可见光波段有431 nm、436 nm和461 nm, 其他波长集中在近红外和短波红外波段。

设置CARS算法为5折交叉验证, 最大主成分为10, 重采样次数为50。在指数衰减函数的作用下, 选择的变量随着运行次数的增加呈指数减少, 其交叉验证均方根误差值呈现先减小后增加的趋势, 在第22次采样有最小的交叉验证均方根误差值, 即选择的变量子集为最佳。最终选择的变量数为79个, 选择的波长在447 nm、1100 nm、1440 nm、1900 nm、2000 nm和2300 nm附近。

SPA对FD进行特征波长筛选时, 随着筛选变量的增加, RMSE先是迅速下降, 当变量数为92时, RMSE值达到稳定。SPA筛选出来的波长分布较广, 大致分布在400～476 nm、804 nm、838 nm、842 nm、1080 nm、1354 nm等波长。总体来说, 选择的特征波长分布在近红外和短波红外波段的居多, 其余区间分布较少。

3.4 反演模型构建

3.4.1 随机森林反演模型构建

随机森林回归将相关性分析和在相关性分析基础上进一步进行IRIV、CARS和SPA算法获得的特征波长作为输入自变量, 氧化铁含量作为因变量进行回归建模。为了模型更稳定和结果更可靠, 模型训练过程中使用5折交叉验证, 同时对随机森林中的最重要的两个参数: n_estimators (决策树的数量)和max_depth (树的最大深度)调优利用网格搜索法实现。随机森林反演结果如表2所示。

从表2可以看出, 不同的光谱变换以及不同的波长筛选方法, 精度有一定的差异。从光谱变换对精度的影响来看, FD的验证集精度总体比其他光谱高, OR和RL均较低, 不能很好地预测土壤氧化铁含量, 这与FD能够消除背景噪声的干扰同时对重叠混合光谱进行分解有关。从波长筛选对精度的影响来看, OR和FD中, CC-CARS比CC、CC-IRIV和CC-SPA的验证集精度高, 且变量数量也比较少, 具有较好的预测性和稳定性。OR中, CC-CARS的验证集最高, 达0.5, 但只能粗略估算氧化铁含量。而CC、CCIRIV和CC-SPA的均未达0.5, 不能进行该地区的氧化铁含量估算; FD中, 4种波长筛选方法均能预测土壤氧化铁的含量, 其中CC-CARS精度最高, 建模集达0.929, RMSE为2.832 g∙kg, 验证集为0.709, RMSE为7.252 g∙kg, RPIQ值为2.794, 具有很强的预测能力; RL中, 4种波长筛选方法的验证集均未达0.5, 不能进行该地区的氧化铁含量估算。

表2 土壤氧化铁含量随机森林反演结果Table2 Random forest inversion results for soil iron oxide content

3.4.2 偏最小二乘反演模型构建

偏最小二乘回归建模输入的自变量和因变量和上述随机森林回归建模相同, 建模过程中的交叉验证采用逐一剔除方法。偏最小二乘回归建模结果如表3所示。

表3 土壤氧化铁含量偏最小二乘反演结果Table3 Partial least squares inversion results for soil iron oxide content

从光谱变换对精度的影响来看, RL的验证集中除CC-IRIV之外, 其他波长筛选方法的精度均比OR和FD要高。这与RL能够放大吸收特征, 降低背景噪声的影响有关。从特征波长筛选对精度的影响来看, OR中, CC-IRIV的验证集最高, 为0.570, 其次为CC-SPA、CC和CC-CARS, 验证集分别为0.569、0.540和0.498; FD中, CC-IRIV的验证集最高, 达0.609, 其次为CC、CC-CARS和CC-SPA, 验证集分别为0.496、0.485和0.447; RL中, CC-CARS处理后精度最高, 其建模集和验证集分别为0.833和0.826, RMSE分别为4.361 g·kg和5.600 g·kg, RPIQ为3.618, 具有优秀的模型稳定性和预测能力, 其次为CC-SPA、CC和CC-IRIV, 其验证集分别为0.745、0.632和0.420。所有模型中, OR-CC-CARS-PLSR、FDCC-PLSR、FD-CC-CARS-PLSR、FD-CC-SPA-PLSR和RL-CC-IRIV-PLSR验证集均低于0.5, 且其RPIQ值分别为2.128、2.125、2.102、2.029和1.980, 不能预测土壤中氧化铁含量。

3.4.3 模型对比

为了对比土壤氧化铁在不同模型的反演效果, 选择两种建模方法中预测效果最好的模型, 绘制了验证集的氧化铁实测值与预测值的1∶1散点图。如图6所示。

图6 验证集土壤氧化铁含量实测值与预测值散点图Fig.6 Scatter plot of the observed and predicted values of soil iron oxide content in the validation set

从图中可以看出, 两种模型的预测值和实测值基本分布在1∶1线周围, 其中RL-CC-CARS-PLSR更接近1∶1线, 模型的预测精度最高。通过综合比较建模集和验证集的精度, 认为RL-CC-CARS-PLSR模型的预测精度更高, 模型更加稳定, 因此可作为预测该地区土壤氧化铁含量的最佳模型。

4 讨论

本研究以云南禄丰恐龙谷南缘地表土壤的室内高光谱数据结合氧化铁含量的化验数据, 反演该地区的氧化铁含量。研究表明, 原始光谱反射率与氧化铁含量呈负相关性, 这与何挺等的研究结果一致。

RF建模中, FD的整体精度要比OR和RL高。而PLSR建模中, RL整体精度比其他光谱高。这说明, 在面对土壤的复杂性、氧化铁含量变异系数较高时, 组合不同的光谱变换建模方法, 可通过比较最后的精度选择最佳的光谱变换和建模方法。总体来说, PLSR比RF具有更高的预测精度, 因为PLSR可以运用主成分分析提取到对土壤氧化铁含量相关性最大的波长, 从而减少了光谱维度。这与众多研究成果相似。

由于高光谱的波长冗余性, 进行特征波长提取是必要的。不同研究者在特征波长提取方法中存在差异, 包括利用不同的光谱变换与氧化铁含量进行相关性分析提取特征波长, 或者在相关性分析基础上利用逐步回归和主成分分析进行特征波长的选取。本文首先利用相关性分析进行特征波长的粗选, 然后在相关性分析的基础上利用IRIV、CARS和SPA算法进一步进行特征波长的精选。基于此方法的RF和PLSR模型均比仅基于相关性筛选变量的模型精度高, 且筛选出来的波长数远比基于相关性筛选的波长少, 不仅降低了模型的复杂度, 而且提高了模型精度。不同的波长筛选方法, 模型精度有差异。最好的波长筛选方法为CC-CARS, 筛选出的波长大部分都分布在500 nm、700 nm和900 nm处, 这与铁的吸收峰一致, 其他波段分布在1400 nm、2200 nm处, 这是由于受到各种功能基团的影响, 与前人的研究一致。因此该方法可以作为一种有效的波长筛选方法。

土壤光谱是土壤属性的综合反映, 本文在估算氧化铁含量时, 未考虑有机质对氧化铁的影响。根据图6的实测值与预测值散点图发现, FD-CC-CARSRF模型在氧化铁含量低于30 g·kg时, 预测值不稳定, 使预测值偏高。这可能是氧化铁含量较低, 光谱反射率受有机质影响导致的。氧化铁含量较低时, 光谱反射率将由有机质主导。根据本次研究区氧化铁含量统计分析, 存在部分点位氧化铁含量较低, 对于有机质是否对该点位土壤的光谱反射率造成影响暂不可知。因此下一步研究要考虑有机质对氧化铁的影响, 以实现更高精度的反演。

对于研究区不同土壤类型, 仅利用原始光谱进行分析达不到精度要求, 往往需要进行光谱变换。结合相关系数法和波长选择算法, 可提升模型的预测精度。对比发现RL-CC-CARS-PLSR模型的预测精度最高, 这与倒数的对数光谱可以放大光谱的吸收特征和降低背景噪声的影响有关, 同时, 由于氧化铁含量与光谱反射率呈负相关, 且CARS算法是根据PLS回归系数来进行波长的选择, 因此线性模型在进行土壤氧化铁含量的反演有着更好的表现。该方法可为不同土壤类型的氧化铁含量反演提供参考。

5 结论

本研究以云南禄丰恐龙谷南缘为研究区, 以研究区的紫壤、红壤和黄棕壤为研究对象, 利用实测高光谱土壤数据和土壤氧化铁含量信息, 选用3种光谱和4种变量筛选方法, 结合随机森林和偏最小二乘回归对氧化铁含量进行定量反演。结果表明, 将OR、FD和RL与氧化铁含量分别进行相关性分析, 发现利用CC进行特征波长的筛选后, 保留的波长数仍然过多, 这对后续模型的精度带来负影响。在相关性分析基础上利用IRIV、CARS和SPA算法可以使特征波长数量进一步减少。

随机森林和偏最小二乘回归建模, 分别以筛选后的特征波长为自变量, 氧化铁含量为因变量构建氧化铁的反演模型。通过结合光谱变换、波长筛选方法以及建模方法, 反演模型结果有所不同。随机森林构建的模型中, FD结合特征波长选择算法构建的模型较优, 其中FD-CC-CARS-RF的建模集为0.929, 验证集为0.709, RMSE为7.252 g·kg, RPIQ为2.794。偏最小二乘回归中, RL结合特征波长选择算法构建的模型较优, 其中RL-CC-CARS-PLSR的建模集为0.833, RMSE为4.361 g·kg, 验证集为0.826, RMSE为5.600 g·kg, RPIQ为3.618。综合1∶1散点图, 发现RL-CC-CARS-PLSR的实测值和预测值均比较靠近1∶1线, 具有较好的预测结果。因此以RL-CC-CARS-PLSR为该地区的氧化铁最佳反演模型。