基于高光谱技术的退耕还林地年限判别

邓永鹏，朱洪芬，丁皓希，孙瑞鹏，毕如田

（山西农业大学资源环境学院，太原 030000）

0 引 言

黄河中游的大部分支流均处于黄土高原区，面积约为34.4×10km，是中国分布面积较大的生态脆弱区。为改善黄河中游的生态环境问题，中国于1999年试行，并于2000年正式开始实施退耕还林（草）工程。随着退耕还林工程长久实施，研究不同年限退耕还林土壤的综合性质对评估退耕还林工程的效果具有重要意义。由于部分地区无法提供有效的退耕资料，因此快速、准确地识别退耕还林地的年限对评估退耕还林工程的生态效益具有重要价值。

近年来，光谱分析技术被广泛应用到中药材种植年限识别、储藏年限识别，野生食用菌储藏年限识别、产地识别，农作物种类区分、茶树种类识别、土壤类型分类等不同领域内。陈泽炎等对目前中药材生长年限鉴定方法进行了总结，发现光谱分析技术可以高效、快速地识别不同种植年限的中药材。卜海博等采用近红外光谱法对不同生长年限的林下山参和园参进行种类和年限识别，结果表明近红外光谱法准确可靠、快速无损。杨天伟等采用傅里叶变换红外光谱技术，准确区分了不同产地、不同年份美味牛肝菌样品，为野生食用菌的鉴别分类提供了技术参考。虞佳维等以黄淮海地区玉米、小麦和杨树三种主要植被为研究对象，采用高光谱技术进行区分，结果表明光谱技术可有效用于不同植被的区分。翁海勇等采用近红外光谱分析技术实现了茶树叶片表没食子儿茶素没食子酸酯（Epigallocatechin Gallate ，EGCG）含量的快速、准确识别，从而为高EGCG含量茶树的育种提供技术支持。史舟等基于中国西藏、新疆、黑龙江、海南等地的16种土类的1 581个土壤样本，以光谱反射率一阶微分主成分数据为自变量进行土壤光谱分类，并根据分类结果建立了中国不同土壤类型的有机质预测模型。赵小敏等以江西红壤地区的7种红壤亚类土壤为研究对象，对不同亚类及土属的土壤高光谱特征进行分析，选取特征变量进行Fastclus聚类分析，结果表明，土壤亚类分类的准确度为86.23%，土属的准确率仅为66.37%。目前，尚未发现开展退耕还林地不同年限识别的研究，应用光谱分析技术可以为不同年限退耕还林地的分类提供技术参考。

本研究以黄河中游大宁县不同年限退耕还林地土壤为研究对象，以不同退耕年限刺槐林地土壤光谱曲线和多种土壤理化数据为基础，分析退耕年限对土壤光谱及其综合属性的影响。利用倒数的对数（Reciprocal of Logarithm，RL）、一阶微分（First Order Differential，FD）等光谱处理方法，结合主成分分析（Principal Component Analysis，PCA）、K 均值聚类（K-means Clustering Algorithm，K-means）、支持向量机（Support Vector Machine，SVM）和线性判别分析（Liner Discriminant Analysis，LDA）等方法探索不同退耕年限土壤的分类结果，实现基于土壤光谱的退耕年限快速区分，以期为快速准确识别退耕年限提供技术参考。

1 材料与方法

1.1 研究区概况

本研究以山西省大宁县境内退耕还林地为研究区（110°27′55″～111°0′40″E， 36°16′40″～36°36′25″N）。大宁县地处黄河中游，地势南北高、中间低，属典型的黄土残垣沟壑区，是黄河泥沙的主要输入区之一。2000年大宁县被定为黄河中上游退耕还林（草）生态工程试点县之一，配置有刺槐、柠条等树种。查询国家土壤信息服务平台(http://www.soilinfo.cn/map/)中国1:400万土壤类型图可知，研究区内主要土壤类型为黄绵土和褐土。

1.2 土壤样品采集与处理

2020年11月在研究区内分别选择退耕3 a（RF）、8 a（RF）、14 a（RF）和18 a（RF）的刺槐林地各32处，共128个样点（图1a）。由于黄绵土土质疏松，易受侵蚀，水土流失严重，是实施退耕还林工程的主要土壤类型，因此样点土壤类型均为黄绵土。土壤样品采集选用“S”形布点法，采集0～20 cm表层土壤，每个采样点取5个点的土壤混合后作为该采样点土样。将土样分为两份，一份土样在室内自然风干、研磨并过2 mm筛，用于光谱测试，土壤有机碳SOC含量测试；一份原状土样过2 mm筛用于土壤质地测定。同时用环刀进行取样，测定饱和导水率、土壤含水率。其中SOC含量测定采用重铬酸钾氧化—外加热法测定，土壤含水率采用经典烘干法测定，饱和导水率利用恒定水头法测定。土壤质地数据测定采用激光粒度仪Mastersizer 3000进行测定，测定粒径范围为0.02～2000m，测定所得土壤粒径数据根据美国农部制粒级分级划分为：黏粒（＜0.002 mm），粉粒（0.002～0.05 mm），砂粒（0.05～2 mm）。

1.3 光谱测试及数据处理

采用美国ASD（Analytica Spectra devices，Inc）公司的FiledSpec 4便携式光谱仪测定土壤反射光谱曲线，其波谱范围为350～2 500 nm。土壤光谱测量在暗室中进行，每次测量开始前将光谱仪预热2 h。将土壤样本放于直径为8 cm，深度为2 cm的培养皿中，用直尺轻轻刮平表面，减小土壤表面粗糙度。光源由一个50 W的卤素灯提供，光纤探头视场角为10°，距离土样表面15 cm。每次测试前采用标准白板校正，每个土样测定10条光谱曲线，取平均值作为该土样的实际光谱曲线。

图1 研究区概况及土壤光谱特征参数Fig.1 General situation of the study area and soil spectral characteristic parameters

去除噪声较为强烈的350～399、2 451～2 500 nm波段。同时，为避免数据冗余，对原始光谱反射率数据（Reflectance，）做10 nm重采样，并对其进行Savitzky-Golay平滑（Savitzky-Golay smooth，SG），对平滑后的光谱数据分别做一阶微分、倒数的对数、去包络线（Continuum Removal，CR）和光谱特征参数（Spectral Characteristic Parameter，SCP）处理，最后对、FD、RL等数据进行主成分分析。其中，CR又称为连续统去除法，可以将光谱反射率归一化到0～1之间，有效突出光谱特征。SCP包括：吸收谷（）、面积（）、深度（DP）、吸收波段波长位置（）、斜率（）、宽度（）、对称度（）等（图1b），对称度=/。相关数据处理在The Unscrambler X 10.4、ENVI 5.3、MATLAB R2019b以及SPSS 26.0中完成。

1.4 分类指标选取

分析不同年限退耕还林土壤光谱特征并进行试验，最终以原始反射率主成分（R-PCA）、倒数的对数主成分（RL-PCA）、一阶微分主成分（FD-PCA）、去包络线主成分（CR-PCA）、SCP作为分类指标，其中SCP包含：第1、2、3个吸收谷的深度（DP、DP、DP），第1、2、3个吸收谷的对称度（、、），第一个吸收谷的面积及前两个吸收谷面积之和（、+），495～595、895～995、1 450～1 600 nm这3个波段去包络线的斜率（、、）。

1.5 模型建立与精度评价

本研究采用K-means、SVM、LDA三种方法构建不同年限退耕还林地土壤光谱分类模型，对比分析三种模型的分类精度。K-means优点是处理大样本时精度较高，本研究先对指标标准化以消除指标量级对结果的影响。SVM对解决中小样本、非线性和高维数据具有独特的优势，本研究选用径向基核函数，惩罚参数（）和Gamma参数通过交叉验证方法进行选择，为3，Gamma参数为0.33；LDA常用于降维和分类，核心思想是投影后类内方差最小，类间方差最大，本文LDA变换的时候主要是对自变量数据即光谱数据进行了变换，重新投影到一个超平面上，使得样本有更大的类间距离，更小的类内距离。分类模型构建分别在SPSS 26.0、MATLAB R2019b和The Unscrambler X 10.4中完成。

将不同退耕年限土壤样本分别按3:1划分训练集和验证集，各退耕年限得到24个建模样本和8个验证样本，故本研究建模集为96个样本，验证集为32个样本。使用混淆矩阵验证模型的分类精度，并且使用总精度和Kappa系数对模型精度进行评价，计算公式如下：

式中为混淆矩阵中列的数量；P为混淆矩阵中第行第列的像元数，表示正确分类的个数；P和P分别为第行和第列总样本个数；代表验证样本的总个数。

2 结果与分析

2.1 不同年限退耕土壤的基本理化性质

从表1可知，研究区总样本土壤有机碳介于2.12～30.11 g/kg之间，变异系数为60.52%，属于中等变异程度，不同退耕年限土壤有机碳含量依次为RF＞RF＞RF＞RF。含水率最大值为25.39%，最小值为3.77%。饱和导水率介于0.01～10.98 mm/min之间，变异系数为219.61%，属于强变异程度。黏粒质量百分数介于4.03%～17.81%之间，变异系数为26.46%，属于中等变异程度。砂粒质量百分数最大值为55.95%，最小值为4.70%，不同退耕年限砂粒质量百分数依次为RF＞RF＞RF＞RF。

表1 土壤样品基本理化性质Table 1 Basic physical and chemical properties of soil samples

2.2 不同年限退耕土壤的光谱特征

为了分析不同退耕年限土壤光谱曲线的特征，对不同退耕年限土壤光谱曲线求平均值，得到不同退耕年限土壤光谱反射率平均值曲线（图2a）。从图2a可以看出，不同退耕年限的土壤光谱曲线形状基本相似，总体上呈现递增的趋势。在可见光波段范围内，随着波长的增加，光谱反射率迅速增大，在近红外波段范围内，随着波长的增加，光谱反射率增加的速度逐渐减缓。总体上反射率依次为RF＞RF＞RF＞RF，与SOC含量呈负相关（≤-0.65）。在可见光（400～580 nm）波段则表现为RF＞RF＞RF＞RF，砂粒含量越高，光谱反射率越大。在 610～900 nm 波段范围内，反射率依次为RF＞RF＞RF＞RF，这可能与土壤中氧化铁的含量有关，有研究指出，在可见光波段内氧化铁含量增加会导致反射率明显下降。本研究对不同退耕年限土壤光谱反射率平均值曲线做去包络线处理，从而深入分析不同退耕年限土壤的光谱特征（图2b）。

图2 不同年限退耕土壤反射光谱曲线及去包络线Fig.2 Reflectance spectrum curves and continuum removal of soil with different years of conversion

从图2b中可以看出，在480、900、1 100、1 400、1 900、2 200、2 350 nm处去包络线都出现了明显的吸收特征。在400～600 nm波段范围内，去包络线值依次为RF＞RF＞RF＞RF，与SOC含量呈正相关，并且表现为吸收谷面积随着去包络线值的增大而减小。在800～1 000 nm和1 000～1 200 nm波段范围内，去包络线值表现出的规律相同，均为RF＞RF≈RF＞RF，造成RF和RF去包络线近乎相同的原因可能是RF和RF土壤中氧化铁含量相近，常认为400～1 100 nm处吸收谷的出现与氧化铁的存在有着一定的关系。在1 400 nm处的吸收谷是受到土壤表面吸附水、黏土矿物O-H结构晶体等吸收引起的，1 900 nm处的吸收谷主要是由土壤中水分子的O-H官能基在1 900 nm处发生一级倍频伸缩震动和转角震动造成的，土壤含水率增加，光谱反射率降低，同时土壤含水率的增加也与黏粒的增加有着一定的关系。对比不同年限退耕土壤饱和导水率可以发现，饱和导水率依次为RF＞RF＞RF＞RF，这与1 900 nm处吸收谷最低点去包络线值的规律相同。在2 200～2 300 nm附近存在着Al-OH黏土矿物（高岭石）的吸收带，从而导致了光谱曲线明显的吸收谷。已有研究表明，不同粒径土壤的光谱差异随波长的增加而增大，因此2 000～2 500 nm处的光谱差异很可能由土壤黏粒、粉粒的差异引起。在2 350 nm处吸收谷的最低点去包络线值依次为RF＞RF＞RF＞RF，与土壤黏粒含量呈负相关。对土壤理化性质和光谱曲线进行综合分析，发现不同年限退耕土壤光谱曲线的特征有较大差异，这些差异不仅仅是单一土壤属性造成的，因此不同年限退耕还林土壤的光谱分类不可以仅依据单一土壤属性进行。

2.3 不同年限退耕土壤的分类

2.3.1 基于K-means的分类

从聚类结果（表2）可以看出，不同输入因子构建的聚类模型，聚类精度有很大差距。以CR-PCA作为输入因子时，模型精度最高，总精度达到81.25%，Kappa系数为0.75。其中RF有1个样本错分为RF，RF有1个样本错分为RF，RF有3个样本错分为RF，仅有1个样本错分为RF，并且除RF的聚类精度仅有50.00%外，其余不同年限退耕土壤聚类精度都达到85%以上，说明以CR-PCA为输入因子，K-means为聚类模型时对不同年限退耕土壤的整体聚类效果较好，对RF的聚类效果则较差。以R-PCA作为输入因子时，模型精度最差，总精度为46.88%，Kappa系数仅为0.29，仅有RF的聚类精度达到60%以上，其中RF有2个样本错分为RF，RF有3个样本错分为RF，RF有4个样本错分为RF，RF有2个样本错分为RF。以RL-PCA、FD-PCA、SCP作为输入因子时，模型聚类精度较为均衡，从整体上来看，FD-PCA和CR-PCA是K-means模型中相对较好的分类指标，RF、RF以及RF均以CR-PCA的聚类精度为最高，而RF的聚类精度则以FD-PCA为最高。

2.3.2 基于SVM的分类

从表2中可以看出，以R-PCA、RL-PCA为输入因子时，模型精度较低，总精度和Kappa系数分别为50.00%、46.88%和0.33、0.29，其中RF的分类精度均只有25.00%，为本研究最低的分类精度，并且有3个样本错分为RF，2个样本错分为RF，1个样本错分为RF，说明以SVM作为分类模型，R-PCA和RL-PCA作为输入因子时，RF的分类效果较差。以FD-PCA作为输入因子时，模型精度最高，总精度为84.38%，Kappa系数为0.79，不同年限退耕土壤的分类精度都达到了75.00%及以上，RF的分类精度更是达到100%，其中RF有1个样本错分为RF，有1个样本错分为RF，RF有1个样本错分为RF，有1个样本错分为RF，RF有1个样本错分为RF，说明FD-PCA可以显著地提高不同年限退耕土壤的分类精度。同K-means分类模型相似，以CR-PCA、SCP为输入因子时，模型分类精度较为均衡。总体上来看，以SVM作为分类模型时，FD-PCA会显著提高不同年限退耕土壤的分类精度，CR-PCA、SCP次之。

2.3.3 基于LDA的分类

由表 2可知，分类结果总精度依次为CR-PCA＞FD-PCA=SCP＞R-PCA＞RL-PCA。对比以R-PCA和RL-PCA为输入因子时不同退耕年限的分类精度可以发现，除RF的分类精度达到80%以上，其余退耕年限土壤的分类精度均小于65%，RF的分类精度只有37.5%，并且RF错分的样本大多为RF和RF。以CR-PCA作为输入因子时，除RF的分类精度较低外，其余退耕年限土壤的分类精度均达到85%以上，具有较好的分类效果，其中RF有1个样本错分为RF，RF有1个样本错分为RF，1个样本错分为RF，RF有1个样本错分为RF。以SCP作为输入因子时，RF和RF分类精度差异较大，其中RF分类精度为100.00%，RF分类精度为62.50%，RF中有1个样本错分为RF，2个样本错分为RF。以FD-PCA作为输入因子时，模型分类效果较为均衡。从整体来看，LDA作为分类模型，对RF的分类效果较差；以CR-PCA和FD-PCA为输入因子时，分类精度有较大提升。

表2 不同模型分类结果Table 2 Classification results of different models

2.4 分类精度比较

以R-PCA、RL-PCA为输入因子时，精度最高的模型分别是LDA模型和K-means模型，总精度和Kappa系数分别为59.38%、0.46和62.50%、0.50。以FD-PCA、CR-PCA、SCP为输入因子时，除FD-PCA以SVM模型精度最高，其余均为LDA模型精度最高，总精度分别为84.38%、87.50%、81.25%，Kappa系数分别为0.79、0.83、0.75，一定程度上提高了不同退耕年限土壤光谱分类的精度。在K-means、SVM、LDA模型中，以FD-PCA、CR-PCA为输入量的模型要略优于以其他三种输入量的分类模型。从整体上来看，以CR-PCA为输入因子，基于LDA模型构建的分类模型总精度最高，达到87.50%；以FD-PCA为输入因子，基于SVM构建的分类模型精度次之，也达到了84.38%；以R-PCA为输入因子构建的K-means分类模型和以RL-PCA为输入因子构建的SVM分类模型总精度最差，仅仅只有46.88%；通过对比五种输入因子在三种模型中的分类总精度可以发现，以CR-PCA和FD-PCA为输入因子的分类总精度均达到了75%以上，其中CR-PCA总精度最大达到87.50%，FD-PCA总精度最大为84.38%，并且对比不同年限聚类精度，CR-PCA除K-means模型中的RF聚类精度较低外，其余聚类精度均要优于FD-PCA，因此在本研究中CR-PCA是区分不同年限退耕土壤的最优输入因子。

3 讨 论

土壤光谱曲线是土壤属性的综合体现。对不同年限退耕土壤的光谱曲线进行分析，发现反射光谱曲线形状大致相似，反射率依次为RF＞RF＞RF＞RF，与SOC含量呈负相关，与南锋等研究结果相符，南锋等对黄土高原煤矿复垦区农田SOC进行预测，发现SOC含量与反射率呈负相关。在400～580 nm波段范围内，反射率大小与土壤砂粒含量呈正相关，这与张雅梅等的研究结果不符，张雅梅等对土壤质地不同粒级颗粒含量进行统一预测，发现在380～1 100 nm波段范围内，砂粒含量最高的砂土反射率最低。这可能是SOC含量和砂粒含量共同作用产生的结果。对光谱曲线做去包络线处理，发现在480和900 nm处，随着SOC含量的增加，吸收谷变浅，面积减小，这与赵明松等的研究结果相一致，赵明松等对江苏中部水稻土和潮土的去包络线进行分析，发现480、900 nm处SOC含量增加导致吸收谷减小。在800～1 200 nm波段范围内，土壤光谱曲线也表现出一定的差异，这可能由于氧化铁含量的不同所造成的。去包络线在1 400、1 900 nm处的吸收谷也表现出显著差异，造成这一差异的主要原因可能是由于土壤水分的不同。

目前研究多采用高光谱非成像技术、成像技术及高分影像等对土壤类型进行区分，未能考虑诸如退耕还林地不同年限的区分。本文以地面高光谱数据为基础，实现了退耕还林地不同年限分类研究，可以为高分影像分类时的光谱变换形式和波段选择提供一定参考，并且后续研究中可以将图像与光谱结合，探究分类精度更高的模型。本研究以R-PCA、RL-PCA、FD-PCA、CR-PCA和SCP为输入指标，基于K-means、SVM和LDA构建退耕还林地不同年限分类模型。通过对比不同分类模型的精度，发现以CR-PCA为输入指标构建的LDA分类模型精度最高。造成这一结果的原因可能是LDA属于监督学习的线性分类方法，对线性问题有较好的处理效果。比较不同分类模型中输入指标的精度发现，CR-PCA在K-means、SVM和LDA模型中的分类精度依次为81.25%、78.13%和87.50%，这充分说明对光谱数据进行CR处理可以显著增强土壤属性对光谱曲线产生的综合影响，进而提高模型的精度，其中主要波段范围为405～595、805～995、1 355～1 495、1 835～2 135、2 145～2 285 nm，后期利用高分影像数据时可以着重考虑选择这些波段范围。孟祥添等等采用高分5号影像数据，对原始光谱曲线进行去包络线、主成分分析等处理，进而区分东北典型黑土区土壤类型，发现CR-PCA分类精度较原始反射率提高了9.15%，说明去包络线可以增强光谱差异性，与本研究结果一致。

本研究将高光谱技术应用到退耕还林土壤研究中，并实现了退耕还林地不同年限识别，为研究退耕还林土壤属性变化及影响提供了一种思路和方法。然而本研究构建的模型较为简单，未构建目前广泛应用的机器学习模型，会对分类模型的精度有一定的影响，后续将会对分类模型和输入指标做进一步的研究，为退耕还林对土壤属性影响的深入研究及退耕还林工作奠定基础。

4 结 论

为了明确退耕工程对土壤属性的影响，以及不同年限退耕土壤的光谱特征，实现快速获取退耕年限的目的，本研究以大宁县不同年限退耕还林地土壤为研究对象，获取土壤理化性质及光谱曲线数据，探讨不同年限退耕土壤理化性质的变化及对光谱曲线造成的影响，同时构建分类模型，选取最优分类模型及最优输入因子。主要结论：

1）随着退耕年限的增长，土壤有机碳含量逐渐增加，砂粒含量先增加后减少，土壤有机碳含量与原始反射率大小呈负相关。

2）不同年限退耕土壤原始光谱曲线形态相似，整体呈递增的趋势；进行去包络线处理后显著提升光谱曲线的吸收特征，在480 nm处出现了由土壤有机碳引起的吸收谷，900和1 100 nm的吸收谷可能是由于氧化铁含量不同导致的，1 400和1 900 nm主要是由于土壤水分引起的吸收特征，2 200和2 350 nm处的吸收特征可能是由于不同粒级颗粒含量的不同所导致的。

3）以去包络线主成分为输入量，基于线性判别分析构建的分类模型精度最高，达到87.50%，为本研究最优分类模型；以去包络线主成分为输入因子的分类模型，总精度均达到了75%以上，最大为87.50%，说明在本研究中去包络线主成分是区分不同年限退耕土壤的最优输入因子。