基于CARS算法的不同类型土壤有机质高光谱预测

唐海涛，孟祥添，苏循新，马 涛，刘焕军,4，鲍依临，张美薇，张新乐※，霍海志

（1. 东北农业大学公共管理与法学院，哈尔滨 150030；2. 黑龙江省地质资料档案馆，哈尔滨 150030；3. 黑龙江省第五地质勘察院，哈尔滨 150030；4. 中国科学院东北地理与农业生态研究所，长春 130012）

0 引言

土壤有机质（Soil Organic Matter，SOM）可以通过生物合成和分解，改善土壤的物理、化学和生物特性[1]，在控制土壤功能和质量、抵消温室气体排放、完善全球碳循环系统信息等方面发挥着重要作用[2]。高光谱预测模型为实现SOM等土壤属性速测与遥感反演以及表层碳库估算等提供数据信息[3]，并为SOM速测仪器研制、土壤制图与退化监测、精准农业实施等提供数据与技术支持[4]。高光谱技术具有精细的光谱分辨率，可获取地物纳米级的连续光谱信息，SOM具有多种官能团（如羟基、羧基等），分别在红外光谱区域有特征性吸收，且不同波段的吸收强度与该物质的分子结构及浓度存在对应关系，因此，红外光谱可以反映SOM含量，为其定量估算提供了一种有效的手段，为预测SOM提供了可能[5]。黑龙江省海伦市位于世界三大黑土地分布区之一的松嫩平原东北端，土壤类型多样，其中黑土面积达到全市面积1/2以上，且是中国重要的商品粮基地，了解其SOM的分布情况、空间变化规律，有利于科学评价土壤的质量情况并对农场合理施肥提供指导，对耕地资源的可持续利用具有十分重要的实际意义，可为海伦市耕地的可持续利用和土壤质量保护监测提供技术支持，为将来海伦市土地管理建立完整的空间土壤信息系统提供框架。

以往室内高光谱对于SOM的输入变量研究多停留在以全波段反射率或对应的数学变换上，选取相关系数较大的波段进行建模，该方法仅考虑了SOM与光谱间的关系，并没有考虑光谱间的重叠吸收或相互影响[6]。利用光谱指数技术预测SOM的研究成为当前热点，光谱指数是由几个窄波段或宽波段组合而成，可通过分析特定波段间的相互作用，提高对待测属性的敏感程度[7]，有助于挖掘波段间的隐晦信号[8]。SOM空间分布特征受到高程、坡度、坡向等地形因子不同程度的影响，地形条件影响其物质循环过程和强度[9]，通过数字高程模型（Digital Elevation Model，DEM）提取高程作为模型辅助变量参与建模。同时特征波段选择是进行SOM含量预测的一个重要方面，已经引起了越来越多学者的关注。土壤光谱反射数据通过竞争自适应重加权采样（Competitive Adaptive Reweighted Sampling，CARS）筛选出的特征波段不仅将输入波段压缩至全波段数目的一半以下，同时提升了模型估测精度，降低了变量维度和模型复杂度[10]，Vohland等[11]发现，在60个农业样品的土壤属性预测中，CARS算法减少了建模时间，且能够合理、精确、有效的确定特征波段在全波段中的位置。以往的学者多以一种类型的土壤为对象，进行SOM高光谱响应特性研究，但是由于土壤的光谱反射率是土壤内在理化特性光谱行为的综合反应，不同类型土壤的光谱特征不同[12]，因此模型普适性较弱。卢艳丽等[13]利用不同土壤类型分组试验分析了东北平原土壤光谱反射率曲线形状变化，确定了8种不同类型土壤与原始光谱反射率的相关敏感波段并建立了同质性SOM预测线性模型，从而达到简化SOM预测模型的目的。Bao等[14]对比了多种土壤分组策略下SOM的预测精度，同时引入竞争自适应重加权采样方法进行模型输入量的筛选，证实了土壤分类的优势与多输入量降维的必要性。因此，不同类型土壤分别提取输入变量进行高光谱SOM预测将有利于分析各类土壤的理化性质，从而提高预测精度。

已有SOM高光谱预测研究常基于一种土壤类型建立模型或者多种土壤类型进行全局回归建模，且输入变量的类型较为单一，存在SOM预测精度不高的情况[15]。为了充分考虑土壤光谱信息及影响因素，本研究以海伦市为研究区域，根据全国第二次土壤普查结果及对采样点的地理位置对土样进行分类。在土壤分类的前提下，以土壤光谱反射率数据、DEM数据以及光谱指数作为输入变量，建立基于随机森林算法（Random Forest，RF）的分类高光谱SOM预测模型。为了降低输入量之间的共线性，引入CARS算法筛选特征波段，提高不同类型SOM预测的精度，以期实现动态快速预测SOM含量。

1 材料与方法

1.1 研究区概况

海伦市位于松嫩平原的中心地带，地理位置在46°58"N～47°52"N，126°14"E～127°45"E之间，属温带大陆性季风气候，地势平坦，土质肥沃，耕地面积广阔，是国家重要的商品粮基地。其土壤类型主要为黑土、草甸土和沼泽土，在该研究区内还有少量的水稻土、暗棕壤及白浆土。黑土土层深厚，结构良好，富含SOM和腐殖质，自然肥力高。沼泽土所处的地势大都比较低洼，SOM累积明显。由于该区地形高程差较大，加上耕地的长期粗放利用导致土壤侵蚀严重，降水将地势较高的土壤冲积到地势较低的草甸土表面，导致表层草甸土性质较为复杂多样[16]。海伦市主要土壤类型（全国第二次土壤普查结果）及采样点分布图和海伦市30 m空间分辨率的DEM数据见图1。

1.2 土样样品

2019年5月15—20日，于作物出苗前，沿主要乡级以上道路，在黑龙江省海伦市全市进行样本采集。选择土壤裸露的地区作为样区，考虑土地利用类型和土壤类型采集0～20 cm耕层土壤。为保证采样点的有机质含量能够代表采样点附近一定空间内的SOM水平，采用四分法收集样品，同时利用GPS记录采样点经纬度，总共采集土壤样本548个。采集的样品经过风干，研磨，过2 mm筛。每个样品分2份，一份用于光谱测量；一份用于SOM含量分析。SOM含量用高温外热重铬酸钾氧化容量法测定[17]。

1.3 光谱测量及数据预处理

采用ASD FieldSpec○R3便携式光谱仪在暗室内对风干土进行光谱测试。光谱测试流程详见文献[18]。由于反射率波谱在400～430和2400～2500 nm范围内噪声较为强烈，为减少高频噪声的干扰，本文选取光谱反射率波谱范围为430～2400 nm，并对其进行9点平滑、10 nm重采样处理，此过程分别在EXCEL和ENVI 5.3中实现。

不考虑土壤空间差异性，将整个土壤样本作为全局回归预测数据集。同时，土壤样本根据全国第二次土壤普查图，利用ArcGIS 10.1中的工具箱提取每个土壤样本的土壤类型，将土壤样本划分为不同土壤类型，同一种土壤具有相同光谱表现特征的土壤样本集。根据中国土壤分类，土壤类型可分为黑土、草甸土、沼泽土，然后针对不同分类样本进行局部回归预测建模。

1.4 输入量提取

国内外学者进行SOM高光谱估测时，输入量多选择为高光谱反射率或光谱吸收特征建立模型，输入变量类型结构单一，容易忽略土壤高光谱反射率之间的高度共线性[19]。本研究通过CARS算法挑选的特征变量、光谱指数结合DEM数据作为模型输入变量。

1.4.1 光谱指数

在高光谱数据预测SOM的研究中，为了确定敏感的波段，必须从SOM含量信息中获取深度信号，因此光谱指数常作为一个重要指标[20]。本文探讨归一化指数（Normalized Difference Index，NDI）、再归一化指数（Renormalized Difference Vegetation Index，RDVI）、比值指数（Ratio Index，RI）与SOM含量之间的关系。

表1 光谱指数及公式 Table 1 Spectral indices and formula

1.4.2 地形因素

地表微气候、土壤中的水分运动以及物质的重新分配进程，都受到地形的影响[25]。在美国地质勘探局网站（http://www.usgs.gov/）下载DEM数据，其空间分辨率为30 m。在ArcGIS 10.1中，利用Spatial Analyst Tools中的Extract Multi Values to Points工具，提取出每个采样点的DEM，将DEM作为模型的输入变量。

1.4.3 竞争性自适应加权算法

土壤高光谱数据量大、存在光谱信息冗余和重叠现象，通过CARS算法挑选特征变量可以降低光谱波段之间的高度共线性问题，从而提高预测模型的精度及速度。CARS算法将各波段变量作为单一个体，在进行个体选择的过程中，保留具有较强适应能力的个体。其具体步骤为：首先，随机抽取固定比率的样本作为校正集建立PLS模型，计算回归系数的绝对值和每个波段点对应的权重，然后利用指数衰减函数（Exponentially Decreasing Function，EDP）和自适应重加权采样法（Adaptive Reweighted Sampling，ARS）对变量进行选择，通过交叉验证的方法计算交叉验证均方根误差（Root Mean Square Error of Cross-Validation，RMSECV），N次蒙特卡罗采样后选择N个子集，得到N个RMSECV，选择RMSECV最小的波段子集，该子集所包含的变量即为最优变量组合[14,26]。本次试验在MATLAB 2014a软件环境中运行CARS算法。由蒙特卡罗交叉验证法选择最优潜在波段变量，其中将蒙特卡罗采样次数设定为100，对采样次数进行反复迭代，通过对比各次采样的RMSECV值，当其值最小时，相应采样次数的变量被筛选为最优变量子集。

1.5 模型构建与验证

RF是基于决策树分类集成算法，其中每一棵树都依赖于一个随机向量，通过对数据集的列变量和行变量观测进行随机化，生成多个分类树，最终将分类树结果进行汇总。RF对于非线性问题有很好的解释能力，降低了运算量的同时也提高了预测精度[27]。本试验在R语言中，利用‘Random Forest’工具包进行预测，在进行拟合前，分别对需要生成树的数量（ntree）参数设定为500，每个节点用于分割节点的预测变量树（mtry）参数设定为1/3总变量数[28]。

模型构建按照建模集与验证集2∶1的比例选取样本。以CARS筛选后土壤高光谱反射率数据、DEM以及光谱指数为自变量，SOM含量作为因变量，运用RF，构建SOM预测模型。使用调整后决定系数（R2adj）、均方根误差（RMSE）以及性能与四分位间隔距离的比率（Ratio of Performance to Interquartile distance，RPIQ）为精度评价指标。R2adj越大、表明模型越稳定；RMSE越小、表明模型精度越高；RPIQ同时考虑了预测误差和观测值的变化，提供了一个更客观、更容易在模型验证研究中进行比较的模型有效性度量。RPIQ越大，模型的预测能力越强。与残差预测偏差不同，RPIQ对观测值的分布没有任何假设[29]，其公式如下：

式中IQ是第三和第一个四分位数之间的差值。

2 结果与分析

2.1 SOM描述统计

土壤样本SOM含量统计特征见表2，质量分数最大值为11.38%，最小值为0.98%，土壤样品SOM差异较大，这为全面解析SOM反射光谱特性研究提供了较完整的样本数据。根据土壤样本SOM描述统计表的偏度和峰度值可以判断SOM含量数据呈现非正态分布。在SOM相关的研究中可知SOM质量分数达到2%以上，对土壤光谱特征起主导作用[30]，SOM质量分数小于2%的土壤，其光谱曲线特征易受其他母质等成分的影响，而本次研究中SOM平均含量（质量分数）4.5%以上，能够充分说明SOM的含量决定了土壤光谱的特征。

表2 土壤样本有机质含量统计结果 Table 2 Statistical results of organic matter content in soil samples

2.2 CARS算法筛选特征波段

3种土壤类型以及未分类整体在指数衰减函数的作用下，优选变量的数量均随迭代次数的增加呈指数减少，其RMSECV值整体均呈现先减后升的趋势。以黑土为例（图 2），从图2a可以看出，随着运行次数增加，被优选出的波段变量数逐渐减少，前5次采样过程有明显递减，此后逐渐平稳。图2b 整体上在1～47次采样中，RMSECV值不断降低，表明筛选过程中剔除的变量与SOM去除量无关，而47次采样迭代以后，RMSECV值呈回升趋势，表明反射率光谱中与SOM无关的大量信息或噪声被添加，从而导致RMSECV值上升。图2c为所有变量在每次采样过程中的回归系数路径变化图，图中各线表示随着运行次数的增加各波段变量回归系数的变化趋势。结合图2b分析发现当采样次数为第47次时，RMSECV值最小即所选择的光谱变量子集最优。草甸土、沼泽土以及未分类整体的RMSECV最小值、相应运行次数及特征波段见表3。

表3 CARS下基于不同土壤类型的特征波段，运行次数和最小交叉验证均方根误差 Table 3 Characteristic wavebands, number of sampling runs and minimal RMSECV of different soil types under CARS

从表3可知，通过CARS算法，黑土、草甸土、沼泽土以及整体未分类分别筛选出23、30、21和9个特征波段，输入波段压缩至全波段数目的16%以下。黑土特征波段的分布主要在1280～2230 nm近红外光谱区域，这是由于受到NH，CH和CO等基团的分子振动的倍频与合频吸收影响[31]，草甸土在可见光-近红外光谱区域均有波段选中，其中1700～1790 nm处SOM响应可能是由氧化铝影响的光谱变化引起的。沼泽土筛选的特征波段在1300～2000 nm比较均匀分布，这主要是由于沼泽土中的大量三氧化物被还原。值得注意的是，波段1450、1470、2150 nm在3种土壤类型中均被选择，这是由于SOM在1400 nm附近受到土壤黏土矿物质中所含羟基的影响，2220 nm附近存在一个与SOM相关的烷烃特征峰和存在的氢氧化铝黏土矿物吸收带影响[32]。沼泽土、草甸土筛选的430、440、530、550、670 nm少量特征波段位于可见光波段，这是由于受到了土壤发色团和SOM本身黑色的影响，可见光波段存在较宽的吸收波段。

2.3 光谱指数的选取依据

光谱指数通过迭代运算，充分考虑波段之间的协同作用，同时最小化无关波段的影响[33]。研究选取的光谱指数是通过文献查阅，选择可用来估测SOM的一系列物理和化学参数的相关光谱指数，并结合本次实际采样点数据进行相关性计算得出。3种土壤类型原始反射率数据与SOM之间的NDI、RDVI、RI指数的相关性均较高，且均通过了P=0.01水平上的极显著性检验（表4）。黑土RI指数与SOM的相关性最高，相关系数为0.757，草甸土RDVI指数与SOM的相关性最高，相关系数为-0.784，沼泽土RDVI指数与SOM的相关性最高，相关系数为0.922。图3是不同土壤类型的3种光谱指数与SOM含量的二维相关系数矩阵图。3种土壤类型的SOM敏感波段区域主要集中于短波红外部分，主要集中在1000、1900和2200 nm附近。

表4 土壤有机质含量与最佳光谱指数的关系 Table 4 Relationship between soil organic matter content and optimal spectral index

2.4 SOM光谱预测模型

由表5可知，黑土、草甸土、沼泽土的验证集调整后决定系数依次为0.678、0.674、0.768，其中沼泽土精度最高，草甸土精度最低，这是由于沼泽土在积水条件下，空气隔绝，微生物活动受到强烈抑制，植物残体不能充分分解，而以粗SOM和半腐烂SOM的形式积累于地表。全局回归模型R2adj达到0.742，局部回归模型R2adj达到0.777。通过局部回归，在一定程度上提高了SOM的预测精度。无论是单一土壤类型，还是整体SOM预测，其R2adj均达到0.67以上，RPIQ均大于1.8，表明该模型能较好实现SOM预测。

表5 不同土壤类型随机森林预测模型精度 Table 5 Prediction model accuracy of random forest for different soil types

在高光谱SOM预测研究中，波段筛选是一个关键方面。本研究通过CARS算法筛选波段与已往学者利用相关分析取相关系数大于0.65筛选出的波段[34]进行建模比较，研究发现CARS算法不仅极大地降低土壤高光谱变量维度和计算复杂程度，验证集R2adj提高了0.167，精度有一定程度的提升。

表6 不同波段筛选方式随机森林预测模型精度 Table 6 Accuracy of random forest prediction model with different band screening methods

3 讨论

本研究将不同土壤类型（黑土、草甸土、沼泽土）分别进行SOM的预测，取得了较高精度。通过土壤分类进行SOM预测，消除了不同土壤类型由于“向邻性”导致的反射光谱曲线相似的影响，从而有利于提高预测精度。由于不同类型土壤中矿物成分与SOM含量的差异，造成反射光谱间存在显著的区别，通过土壤分类，将有利于提取不同类型土壤光谱参数进行SOM预测。陆龙妹等[35]通过全局回归与局部回归进行SOM预测比较，依照传统土壤类型建立各自的有机质光谱预测模型精度并不好，这是由于砂姜黑土和黄褐土2种土壤类型的黏土矿物都存在蒙脱石且含量较高，SOM含量接近，所以2种土壤类型之间光谱曲线特征相似，造成SOM全局回归精度低。而黑土、草甸土、沼泽土之间黏土矿物存在着较大的差异，因此通过全局回归与局部回归比较，全局回归能够提高有效信息的获取程度提高模型精度。其沼泽土的预测精度高于草甸土，这是由于沼泽土土壤湿、土层紧且富有弹性，有机质含量丰富、土体酸碱度从微酸到碱性、土壤颜色深，而草甸土土壤表层砂砾化、有浮沙覆盖、有机质含量较低、土体呈碱性、质地较粗、细粒物质少、土壤色泽浅有一定的关系。

以往许多学者们采用相关分析法研究SOM与土壤光谱反射率（或其不同数学变换形式）的关系，将相关系数高的波段作为SOM敏感波段。而后，越来越多的学者采用CARS变量优选方法，从全波段中滤除无效变量或冗余变量，优选出敏感波段。本研究基于CARS算法，黑土、草甸土、沼泽土分别选择23、30、21个特征变量，占全波段数目的11.6%、15.2%、10.6%，极大地缩减了波段信息，解决了SOM预测研究中波段数目多，计算任务繁重的问题。结果表明，CARS筛选的最优子集存在一定的规律性，波段主要集中在1100～2400 nm之间，这主要由于受到羰基、酰胺和羟基等基团的分子振动的倍频与合频吸收影响。其中，黑土筛选的特征波段少位于1000 nm以下，这是由于CARS是通过利用线性模型偏最小二乘法作为适应度函数，及交叉验证不断优化计算，最终选择出最优子集而不是常用的相关性分析确定特征波段。已有的相关研究表明：SOM在整个NIR-SWIR范围比较敏感，李稳冠等[26]将栗钙土、黑钙土、灰钙土、山地草甸土等土壤光谱曲线通过CARS挑选的特征波段，变量主要分布在1900～2400 nm的近红外光谱区域，在可见-近红外光谱区域均有分布。CARS对原始光谱进行特征变量筛选，在保证模型精度的同时显著减少构建模型的变量数。Bao等[14]对黑土、黑钙土、风沙土、草甸土4种土壤类型通过CARS算法筛选最优变量子集，其波段大多位于1350～2400 nm范围内，少量位于400～1200 nm。因此，通过CARS算法筛选的特征波段，与已有研究SOM的反射光谱响应波段相吻合。不同土壤类型通过CARS筛选的最优子集也存在差异，其选择的特征变量具有不稳定性。

通过耦合敏感波段的反射率数值进行数学变换所计算得到的光谱指数，避免了由于原始反射率作为输入量所造成的数据冗余，以及产生的共线性问题。黑土筛选出的波段主要为1030、1910、1940、1950 nm，草甸土在1420、1340、2150、2230 nm，沼泽土集中在1920和1930 nm。3种土壤类型的筛选的波段都位于NIR-SWIR范围，这是由于羰基基团的基频振动和其在NIR-SWIR范围所对应的酰胺、羟基等基团倍频和合频吸收影响，也与以往的研究一致[36]。因此通过将不同类型土壤分别，以CARS筛选的特征波段、DEM数据和光谱指数作为数据源，建立的RF模型能够有效实现SOM预测，使精度有着显著的提升。然而，本次研究仍存在不足之处：土壤的光谱反射率还会受到土壤的成土母质、矿物成分、土壤表面粗糙度、粒径、水分等因素的影响，因此，后续研究在原土室外光谱的基础上，将考虑更多的影响因素，加强原土室外光谱SOM的估测模型研究，以提升SOM的预测精度。

4 结论

为了解决不同类型土壤预测有机质（Soil Organic Matter，SOM）输入量类型单一造成精度偏低的问题，本文以海伦市3种土壤类型（黑土、草甸土、沼泽土）的室内光谱反射率为研究对象，结合数字高程模型（Digital Elevation Model，DEM）以及光谱指数作为输入量，运用随机森林算法（Random Forest，RF）进行SOM预测，得出以下结论：

1）通过竞争自适应重加权采样（Competitive Adaptive Reweighted Sampling，CARS）算法，筛选出的特征波段不仅将输入波段压缩至全波段数目的16%以下，而且能够在很大程度上降低土壤高光谱变量维度和计算复杂程度，从而提高了模型的预测能力。光谱变量经CARS算法筛选后模型调整后决定系数提高0.167，估测效果更好。说明CARS算法在提取特征关键波段变量、优化模型结构方面起到关键作用。

2）通过土壤分类进行SOM预测，不同土壤类型的SOM调整后决定系数存在差异，沼泽土的调整后决定系数最高为0.768，黑土次之，草甸土的预测精度最低，只有0.674，运用RF对3类土壤的SOM预测性能与四分位间隔距离的比率均大于1.8，说明无论是黑土、草甸土还是沼泽土，该模型都有一定的可信度，具有较好的预测能力。

3）通过将CARS筛选的特征波段、DEM以及光谱指数作为输入量，运用RF模型，SOM的局部回归模型验证集精度最优，调整后决定系数为0.777，且RPIQ达到2.689，与全局回归模型相比，模型的验证精度提高了0.035。研究表明，3种类型的输入量，进行单一土壤类型分别建模和全局回归建模，其均具有较好的预测能力，在一定程度上可为以后不同土壤类型SOM预测时输入量的选择提供帮助，从而促进区域不同类型土壤进行SOM预测研究的进展，为农业和环境领域SOM的动态监测和建模提供理论支撑。