基于高光谱成像的烟田土壤pH估测

张 恒,梁太波,宋效东,江 鸿,郭文孟,戴华鑫,翟 振,冯长春,张艳玲

(1. 中国烟草总公司郑州烟草研究院,郑州 450001; 2. 中国科学院南京土壤研究所,南京 210008; 3. 四川省烟草科学研究所,成都 610041)

【研究意义】土壤pH是评估土壤质量的重要指标[1]。土壤的物理、化学和生物特性均会受到土壤pH的影响[2],土壤pH不仅影响土壤矿物质和其他养分的生物有效性,也影响微生物活动,改变土壤有机质分解速度,进而影响土壤养分释放。因此,土壤pH的准确测定对烟田土壤质量监测和土壤改良至关重要。传统的土壤pH实验室测定方法费时费力,且易受测定条件的影响,难以满足智慧农业对田间管理的实时性需求[3]。高光谱遥感技术的发展为土壤成分快速检测提供了新方法,高光谱技术因其快速、无损且可获得被检测物体较高精度的信息,目前已被广泛用于定量估测土壤中的水分含量[4]、养分状况[5-6]、有机质(碳)含量[7-8]以及土壤重金属[9-10]等土壤属性信息。因此可通过构建土壤pH高光谱定量估测模型实现土壤pH的准确快速测定,为烟田土壤精准施肥提供技术支撑。【前人研究进展】近年来,国内外学者在构建土壤pH和高光谱反射率定量估测模型方面进行了一些研究。如蔡海辉等[11]以棉田土壤为研究区,分析9种单一预处方法处理后的光谱数据与pH的相关性,发现光谱全波段经二阶微分处理后结合随机森林算法建模是估测棉田土壤pH的最佳模型;沈从旺等[12]以重庆涪陵区土壤为研究对象,利用相关分析法筛选特征波段构建pH估测模型,发现支持向量机是估测水稻土pH的最佳建模方法,偏最小二乘回归是估测紫色土pH的最佳建模方法;Ji等[13]以浙江省水稻田土壤为研究对象,结果发现原始光谱全波段经Savitzky-Golay平滑处理后,使用偏最小二乘回归建立的土壤pH估测模型效果较好;Yang等[1]以长江中下游平原水稻土为研究对象,在4种建模方法下,对比了使用光谱全波段和遗传算法筛选特征波段对构建土壤pH估测模型的影响。前人选择的研究区域多集中于单个县(区),研究对象多局限于单一土壤类型,所取土壤样品类型和空间分布差异较小,在更大区域尺度上运用高光谱技术对土壤pH进行建模估测的研究较少;且较多研究在估测不同地区土壤pH时,多使用光谱全波段构建土壤pH估测模型,关于主成分分析对光谱数据降维以及竞争自适应重加权采样(Competitive adaptive reweighted sampling, CARS)算法在筛选土壤pH特征波段方面的研究鲜有报道。此外,由于高光谱成像技术可获得被扫描样品区域内所有像素点的光谱数据[14],且能同时扫描分析多个样品,所以高光谱成像技术比使用地物光谱仪基于点获取的高光谱数据更快更准确,然而到目前为止,利用高光谱成像技术估测烟田土壤pH的研究鲜见报道。【本研究切入点】以我国四川省烟田土壤为研究对象,利用高光谱成像技术获取土壤样品的高光谱数据,采用12种光谱预处理方法对原始光谱进行处理,结合原始光谱全波段、主成分分析降维得到的特征光谱以及CARS筛选的特征波段,运用4种建模方法分别构建土壤pH估测模型。【拟解决的关键问题】旨在探索区域尺度下运用高光谱成像技术估测土壤pH的可行性,为四川省烟田土壤pH快速测定提供理论依据和方法参考。

1 材料与方法

1.1 土壤样品采集与处理

2021年在四川省4个地区(凉山州、攀枝花、宜宾、泸州)18个县(区、县级市)共采集土壤样品296份,采样点分布见图1。具体采集方法:采用5点取样法采集烟田0～20 cm耕层土壤样品,四分法保留2 kg,在实验室土壤样品风干箱内风干,剔除杂物,研磨后过2 mm筛,分成两部分,一部分用于土壤高光谱数据的采集,另一部分用于测定土壤pH。用于采集高光谱数据的土壤样品,在采集前平摊于通风条件良好的实验室内72 h,以降低土壤水分对高光谱数据采集过程中的不利影响。土壤pH测定参照《土壤农业化学分析方法》[15],采用水浸提电位法测定pH,水土比为2.5∶1.0。

图1 采样点分布Fig.1 Distribution of sampling points

1.2 室内高光谱图像采集与数据提取

高光谱图像采集装置主要包括消色差镜头(HSIA-OLE23)、可见—近红外高光谱相机(GaiaField-V10E-AZ4)、电动载物台、计算机和Spec View图像采集软件等。光谱采集方法:将处理好的土壤样品放置于直径10 cm、深2.0 cm被黑布包裹的培养皿中,刮平土壤表面后放置在电动载物台上,通过操作电脑中的Spec View图像采集软件,使电动载物台以1.35 cm/s的速度带动土壤样本进入图像采集箱内,在暗箱环境下完成对土壤样品高光谱图像的采集。光谱数据采集范围为390～2561 nm。在390～1030 nm范围内,光谱分辨率为2.6 nm,采集250个波段;在967～2561 nm范围内,光谱分辨率为5.4 nm,采集288个波段。高光谱图像采集后,需要对采集的高光谱图像进行黑白版校正,以消除仪器本身在采集过程中带来的光谱噪声。由公式(1)对原始高光谱图像进行校正。

(1)

式中,I0为采集的原始高光谱图像,W为采集反射率为99%的标准校正白板获取白板标定图像,B为采集反射率为0%的内置黑板获取黑板标定图像,I为校正后原始高光谱图像。

对于每个土壤样本校正后的高光谱图像,选择培养皿内土壤样本区域作为感兴趣区域(Region of interest,ROI)。提取土壤样本ROI内所有像素点的光谱数据,然后对ROI内所有像素点的光谱反射率值进行平均,以此生成平均光谱作为此样品的原始高光谱反射率数据。

1.3 高光谱数据的预处理

由于每段光谱数据的头部和尾部均存在噪声,因此每个土壤样本只使用408～1007和1012～2500 nm的数据作为进一步分析和建模的高光谱数据。采用最大最小归一化(Max-min scaling, MMS)、多元散射校正(Multiplicative scatter correction, MSC)和标准正态分布(Standard normal variate, SNV)3种单一预处理方法对原始光谱(Reflectance,R)进行预处理[16]。同时为比较不同组合光谱预处理方法对建模精度的影响,在一阶导数(First derivative, D1)、Savitzky-Golay平滑(Savitzky-Golay smooth, SG)和趋势校正(DT)基础上,对光谱数据再分别进行MMS、SNV和MSC预处理。

1.4 高光谱数据特征波段提取与降维

为减少高光谱数据的冗余,降低模型的复杂性和计算量。研究中使用CARS算法筛选特征波段和主成分分析(Principal component analysis,PCA)对光谱数据进行降维处理。PCA主要通过某种线性投影将高维数据映射到低维的空间中,可将众多具有相关性的指标重新组合一组新的无相关的综合指标来代替,以达到对原始特征进行降维的目的[17]。CARS是通过适应重加权采样和指数衰减函数筛选出偏最小二乘模型中回归系数绝对值大的波段,去掉权重较小的波段,利用交互验证选出交互验证均方根误差(Root mean square error of cross validation,RMSECV)最小时的波段子集,从而有效寻找出最优变量组合。研究中CARS算法设置采样次数为50次[18]。

1.5 建模方法

采用偏最小二乘回归、岭回归、核岭回归和支持向量机4种方法构建土壤pH估测模型。偏最小二乘回归(Partial least square regression,PLSR)是一种广泛用于土壤高光谱定量分析的线性回归模型[19],可有效降低自变量间多重共线性问题,避免模型过度拟合以及优化模型估测的精度。岭回归(Ridge regression,RR)是常见包含惩罚项的线性学习算法,可减少线性回归模型的方差[20]。核岭回归(Kernel ridge regression,KRR)是岭回归(线性最小二乘与L2范数正则化)与核函数结合起来的算法,不同的核函数会导致不同的空间函数,研究中核岭回归使用高斯径向基(RBF)为核函数进行建模。支持向量机(Support vector machine,SVM)是基于统计学习理论的机器学习算法,通过核函数将数据转换到高维特征空间中去,已被用于解决线性和非线性问题,研究中SVM使用高斯径向基(RBF)核函数[21]。

建模过程中采用系统网格搜索方法对模型参数进行优化选择,通过10折交叉验证,以建模集均方根误差(Root mean square error of training set,RMSET)最小时确定最佳估测模型。

1.6 模型评估

以决定系数(Coefficient of determination,R2)、建模集均方根误差(RMSET)、验证集均方根误差(Root mean square error of validation,RMSEV)、相对分析误差(Residual prediction deviation,RPD)等参数对建立的pH估测模型性能进行评估[22]。模型的均方根误差越小,R2越接近1,表明模型的准确性越高。当RPD≥2时,表明模型可以较准确的对土壤pH进行估测;当1.4≤RPD<2时,表明模型可以粗略对土壤pH进行估测。模型评估参数具体计算公式参考文献[22]。

1.7 数据分析

高光谱数据的采集和校正均使用 Spec View软件完成;感兴趣区选取、光谱数据提取、光谱数据预处理、特征波段筛选以及模型的建立均使用Python3.9软件;使用Microsoft Excel 2016软件制图。

2 结果与分析

2.1 土壤pH统计特征

为保证建模集和验证集样本中pH分布的均匀性,采用梯度质量法[16]对土壤样本集进行划分,首先将所有样本按照土壤pH进行升序排序,然后从小到大每隔3个样品取1个作为验证集样本,把所有样本以3∶1的比例划分为建模集和验证集,如表1所示。土壤总样本中pH变化范围为4.360～8.715,均值为6.299。建模集和验证集的描述统计特征相近,与总样本各统计特征基本在同一水平,且验证集样本pH分布范围在建模集样本pH分布范围之内,说明验证集样本在建模集中分布均匀。

表1 土壤样本pH描述统计

2.2 不同土壤pH高光谱反射率特征

将总样本按照土壤pH大小升序排列,均分为6组后求每组对应pH和光谱数据的平均值(图2)。不同pH的土壤光谱反射率曲线虽各有不同,但变化趋势较为相似。光谱反射率在400～800 nm可见光范围内迅速上升且光谱反射率值相对较小,在800～2500 nm范围内光谱反射率曲线波动上升,在1400、1900和2200 nm处均存在明显的水分吸收谷。此外,在600～2500 nm范围内,不同pH的光谱反射率曲线差异明显,土壤pH越高,土壤光谱反射率越低,这种差异的存在是后续运用高光谱建立土壤pH估测模型的基础。

图2 不同pH土壤原始光谱反射率曲线Fig.2 Original spectral reflectance curves of soil with different pH

2.3 基于光谱全波段的土壤pH估测模型构建

由表2～3可知,与使用原始光谱全波段建模相比,使用单一预处理或组合预处理方法处理后的光谱全波段,结合不同建模方法建立的土壤pH估测模型性能均有不同程度提升。其中使用PLSR方法建模时,基于DT-MMS预处理方法使用全波段建立的土壤pH估测模型精度最高;分别使用SVM、RR以及KRR 3种建模方法时,基于D1-SNV预处理方法使用光谱全波段建立的pH估测模型精度均最高。整体看,不同建模方法得到的土壤pH最优估测模型,均使用光谱组合预处理方法处理后的光谱全波段,说明在相同建模方法下,选择合适的光谱组合预处理方法对原始光谱进行处理,有助于提高土壤pH估测模型的准确度。原始光谱经D1-SNV预处理方法处理后,结合4种建模方法建立的pH估测模型精度均较高,因此在后续分析中选择D1-SNV对原始光谱进行预处理,以降低数据噪声,提高光谱数据有效性。

表2 基于光谱全波段的PLSR和SVM模型构建及验证

表3 基于光谱全波段的RR和KRR模型构建及验证

2.4 基于CARS筛选特征波段和PCA降维的pH估测模型构建

2.4.1 CARS算法筛选特征波段 利用CARS算法筛选D1-SNV预处理后的特征波段。由图3-a可知,随着运算次数增加CARS算法筛选出的特征波段变量个数逐渐减少。经过前10次的运算,筛选的特征波段变量个数由504个降低至182,呈明显降低趋势,随后特征变量个数呈平稳降低趋势。由图3-b可知,在0～50次运算过程中RMSECV波动变化明显,在0～16运算次数过程中RMSECV整体呈降低趋势,在16～50运算次数过程中RMSECV呈明显升高趋势,说明在16～50运算次数过程中与土壤pH无关的噪声信息增多。当运算次数为第16次时,RMSECV最小,筛选的特征光谱波段组合最优。结合图3-a可知当运算次数为16时,对应的特征波段变量数目为93,其中有35个波段分布在可见光400～760 nm范围内,有58个波段分布在近红外800～2500 nm。

2.4.2 确定PCA最佳维度数 图4显示基于D1-SNV预处理方法,不同建模方法下模型评估参数RPD随维度数(1～222)增长的变化规律。可见针对不同建模方法,运用PCA降维过程中得到的最佳维度数差异明显。当维度数为60时,估测土壤pH的PLSR模型验证集RPD达最大值;当维度数是180时,土壤pH的SVM估测模型验证集PRD达最大值;当维度数为最大值176时,估测土壤pH的RR模型验证集PRD达最大值;当维度数为最大值222时,估测土壤pH的KRR模型验证集PRD达最大值。因此选择60、180、176和222分别作为PLSR、SVM、RR和KRR 4种建模方法在后续构建土壤pH估测模型过程中PCA降维的最佳维度数。

2.4.3 土壤pH估测模型的构建 由表4可知,使用PLSR和RR方法建模时,运用CARS筛选特征波段建立的pH估测模型精度均较高,PLSR模型验证集R2、RMSEV和PRD分别为0.758、0.555和2.034;RR模型验证集R2、RMSEV和PRD分别为0.751、0.563和2.002;使用KRR和SVM方法建模时,运用全波段和PCA降维建立的土壤pH估测模型精度均较高,模型验证集的R2均大于0.750,RMSEV均小于0.570,RPD均大于2.0。

整体看,使用全波段建模情况下,土壤pH估测模型精度排序为KRR>SVM>RR>PLSR,采用KRR模型估测效果最佳(RPD=2.105)。经过CARS筛选特征波段后,pH估测模型精度排序为PLSR>RR>KRR>SVM,采用CARS-PLSR模型估测效果最佳(RPD=2.034),RPD较全波段下的PLSR模型提升了0.29。经过PCA对光谱全波段降维后,pH估测模型精度排序为SVM>KRR>RR>PLSR,采用PCA-SVM模型估测效果最佳(RPD=2.040)。可见在全波段、PCA降维和CARS筛选特征波段建模情况下,得到的pH最优估测模型精度相近。但由于CARS-PLSR模型入选变量数目较少、建模方法解释性较强和模型精度较高,故在D1-SNV预处理方法下,结合CARS筛选特征波段使用PLSR方法建立的土壤pH估测模型(D1-SNV-CARS-PLSR)效果最佳。

3 讨 论

3.1 不同光谱预处理方法对土壤pH建模的影响

构建估测模型前对原始光谱进行适当的预处理,可明显降低外界环境(温度、湿度、光线等)和土壤自身(质地、颗粒大小以及土壤类型等)对采集光谱数据产生的干扰,有助于提高模型的估测精度[3]。本研究发现,原始光谱经单一预处理或组合预处理后,建立的土壤pH估测模型精度均有所提升,且不同建模方法得到的估测土壤pH最佳模型,均使用的是光谱组合预处理方法。这与Lin等[23]研究结果一致,其研究表明基于组合预处理方法建立的土壤pH估测模型精度较高。同时研究还发现基于DT-MMS预处理后光谱全波段建立的土壤pH PLSR估测模型均方根误差较小且精度较高,基于D1-SNV预处理后光谱全波段建立的土壤pH SVM、RR和KRR估测模型精度均较高。原因可能是:趋势校正(DT)有效校正光谱数据中的倾斜和曲线漂移[24];最大最小归一化(MMS)有效消除了尺度差异带来的不良影响[25];一阶导数(D1)对重叠光谱进行拆分,增大了光谱曲线中波谷和波峰的特征[22];标准正态分布(SNV)有效降低了因颗粒分布不均或颗粒大小不同产生的散射影响[26];本研究将上述预处理方法组合起来使用,有效消除了原始光谱曲线中的曲线漂移和散射效应,增大了光谱数据的信噪比,从而显著提高了模型估测土壤pH的准确度和稳定性。

图3 CARS算法筛选特征波段Fig.3 Feature bands screened by CARS algorithm

n为PCA降维过程中的维度数。n is the number of dimensions in PCA dimension reduction.

表4 不同变量筛选方法的土壤pH估测模型构建及验证

3.2 CARS筛选的土壤pH特征波段

CARS算法是一种有效的特征波段筛选方法,主要通过逐步去除不重要变量和冗余,从全波段中筛选出用于建模的最优特征波段变量集合。本研究中,与全光谱相比,CARS算法筛选出了93个土壤pH特征波段,占光谱全波段数目的18%左右,其中筛选的特征波段在可见光和近红外区域内均有分布,主要集中在近红外波段800～2500 nm。与Wang等[9]研究结果相似,其研究发现CARS算法筛选土壤pH的特征波段主要在1903 nm以下的可见光和近红外区域;与Vohland等[27]发现CARS算法筛选土壤pH特征波段只限于1915 nm以外的结果相悖。究其原因可能是由于光谱活性土壤成分(如有机质和粘土)的协变性在不同区域土壤中存在显著差异,导致在不同地区使用同一特征波段筛选方法筛选出的pH特征波段存在差异。

3.3 PCA光谱数据降维对土壤pH建模的影响

主成分分析是最常用的数据降维方法,可以保证在信息量尽量不丢失的情况下对高维高光谱数据进行降维处理,以减少建模过程中的计算量和建模时间。研究中发现与使用光谱全波段建模相比,使用PCA对光谱降维后结合PLSR和RR建模方法建立的土壤pH估测模型性能有所提升,但提升程度并不显著。应该注意到运用PCA对光谱降维建模虽然模型性能提升不明显,但可以有效降低建模过程中的运算量,提高模型的拟合效果[12]。此外,使用SVM和KRR 2种建模方法时,基于光谱全波段建立的土壤pH估测模型较PCA 降维建模取得了更好的估测性能,造成这一结果的因素可能是由于PCA降维过程中虽然减少了用于建模的数据量,但也损失了部分重要光谱信息,导致所建土壤pH估测模型效果不理想[28],这一点在本研究中也得到了证实,研究中基于PCA筛选特征波段运用KRR方法建模时,模型评估参数RPD随维度数增长而增加,但由于降维过程中最大维度数限制导致模型精度低于使用光谱全波段建立的KRR估测模型。

3.4 不同建模方法对土壤pH建模的影响

在4种建模方法中,无论使用光谱全波段还是PCA降维后的特征光谱,KRR和SVM 2种建模方法都取得了比PLSR和RR更好的估测准确度。这可能是由于土壤pH与有机质、粘土等光谱活性成分之间关系较密切,且主要存在内在非线性关系[2],而偏最小二乘回归和岭回归属于线性模型,在解决非线性关系时存在一些不足,故导致模型精度相对偏低。但使用CARS算法筛选的特征波段建模时,PLSR和RR 2种建模方法较KRR和SVM取得了更好的估测精度,可能是由于CARS算法主要通过偏最小二乘回归模型从光谱全波段数据中筛选出光谱波段变量组合,故CARS筛选的特征波段与土壤pH之间主要存在线性关系。另外,通过对比分析发现基于D1-SNV预处理方法,使用CARS筛选特征波段建立的土壤pH PLSR估测模型精度较高且计算量较小,模型验证集R2、RMSEV和RPD分别为0.758、0.555和2.034,说明在区域尺度上运用高光谱成像技术估测土壤pH具有可行性。但与蔡海辉等[11]在阿拉尔市对棉田土壤pH建立的最优模型相比,本研究建立的最佳估测模型精度相对偏低,原因可能是由于本研究所选研究区域尺度较大,所取土壤样品理化性质和土壤类型复杂多样,故导致建立的pH最佳估测模型精度偏低,但也应该注意到本研究所建立的土壤pH估测模型适用区域较大。

因此,在未来运用高光谱成像技术对土壤化学成分进行估测研究时,需要充分考虑土壤样本自身的光谱特性和土壤类型,深入分析土壤化学成分与光谱数据间的内在关系,加强对光谱数据预处理方法、特征波段筛选方法以及机器学习建模算法方面的研究[23],以提高高光谱估测土壤pH的准确度。此外,虽然研究中建立的土壤pH估测模型精度相对不太高,但所建模型的适用区域较广,模型的泛化能力较强,可以有效地对四川省植烟土壤pH进行估测,是否对于其它省份或区域的植烟土壤具有普适性,还需要采集更多具有代表性的植烟土壤样品对模型进行完善。

4 结 论

(1)研究区域耕层土壤pH变化范围为4.360～8.715,均值为6.299,变异系数为17.998%,土壤pH呈中等变异性。

(2)使用单一预处理或组合预处理方法处理后的光谱波段建模,有利于提高土壤pH估测模型的精度;其中原始光谱经D1-SNV预处理方法处理后,结合4种建模方法建立的pH估测模型精度均较高。

(3)CARS算法筛选出了93个土壤pH特征波段,筛选的特征波段较均匀的分布在可见光和近红外区域内,且多数集中在近红外波段800～2500 nm。原始光谱经D1-SNV预处理后,基于偏最小二乘回归和岭回归两种线性建模方法建模时,使用PCA降维或光谱全波段建立的土壤pH估测模型精度均低于CARS筛选特征波段建立的估测模型;基于核岭回归和支持向量机两种非线性建模方法建模时,使用PCA降维或光谱全波段建立的土壤pH估测模型精度均高于CARS筛选特征波段建立的估测模型。综合比较下D1-SNV-CARS-PLSR模型估测土壤pH的效果最好,模型建模集与验证集的R2分别为0.802和0.758,RMSET和RMSEV分别为0.504和0.555,验证集RPD为2.034。