1 引言

中红外（Mid-Infrared，MIR）光谱技术的应用有助于土壤的定性及定量分析，可以节约传统土壤分析需要耗费的大量时间，同时也可以减少分析成本。MIR技术不仅能用于土壤肥力分析，如有机质、pH等^［1］，还能用于检测污染物，如重金属、农药等^［2］，可为土壤改良及施肥提供重要的技术支撑。联合国粮农组织（FAO）于2017年成立了全球土壤实验室网络（GLOSOLAN），并在此框架下于2020年发起了关于土壤光谱学的倡议，旨在促进各国及地区的实验室使用光谱分析土壤^［3］。

近年来，国家^［4］、大陆^［5］和全球尺度^［6］土壤光谱库的建立，为土壤分析提供了数据基础。库包含了土壤样本的理化属性数据、光谱数据，以及坐标、土壤类型、土地利用/覆盖等信息。对于土壤的定量分析（如有机碳），MIR通常比可见近红外（visible–Near Infrared，vis–NIR）具有更佳的能力^［7-10］。但当前土壤光谱库的建设以vis–NIR居多，这是由于MIR传感器对温控的要求较高、仪器昂贵，且便携式仪器选择有限。目前极少数的土壤MIR光谱库已经建成，如全球库^［11］、瑞士库^［12］，土壤有机碳（SOC）、pH等关键土壤属性也得到了较好的定量效果，但MIR光谱库在各国及地区的建立仍待推进。

即便有无限容量的库，全局数据建立的通用关系在应用到局部尺度时面临挑战。这是由于全局模型（Global）包含了多样的土壤变异信息，从而难以准确表征甚至会掩盖局部样本的变异，尤其当待测样本来自多个局部时^［13］。开发新的算法和模型对于提高土壤光谱库的性能是一种有效途径。此外，基于相似的光谱具有相似的土壤成分这一假设，开发“局部化”光谱建模的方法也同样重要。先将少量且具代表性的待测样本入库，令其尽可能囊括待测样本集的主要特征，再从库中匹配到与待测集相似的样本集合用于建模的数据驱动策略，是当前主流。相关方法如Spiking^［14］、RS-Local^［15］、C-KS^［16］等已对土壤vis–NIR光谱库进行了测试并取得积极效果，但尚不清楚适用于vis–NIR的方法是否同样也适用于MIR。而且当库已包含一定的局部样本特征信息时，如何选择恰当的样本数用于建立局部模型以获得较佳的精度尚无定论。

本文探索并提出了一种从土壤MIR光谱库中快速确定局部建模集容量的方法，以期提升MIR光谱分析技术的效果，同时响应FAO发起的GLOSOLAN-Spec倡议。

1　材料与方法

1.1　数据来源

研究数据来自国际土壤参比与信息中心（International Soil Reference and Information Center）发布的全球土壤光谱库（Global Soil Spectral Library，GSSL）^［17］。本研究从库中筛选出同时含有MIR光谱及SOC含量的有效数据共计677个土柱样点（不同土层样本共计3755个），采样点广泛分布在除南极洲之外的六个洲（图1）。经风干、研磨、过筛（＜100 μm）后，土样的MIR光谱由Bruker VERTEX 70红外光谱仪在实验室测得，吸收光谱的波段范围为7498~600 cm^-1，光谱分辨率4 cm^-1。该仪器由液氮制冷，配备金板做谱线校准。

图 1. 全球土壤光谱库中677个包含中红外光谱和土壤有机碳的土样分布

Fig. 1. 677 site locations of with MIR spectra and SOC in the Global Soil Spectral Library

下载图片 查看所有图片

1.2　光谱预处理与集合划分

土壤MIR光谱采用常用的预处理流程：先重采样至16 cm^-1并保留4000~650 cm^-1波段范围，再经Savitzky–Golay平滑法（2阶3窗口）和一阶微分对光谱进行噪声消除和信号增强处理。

光谱数据集的划分使用Kennard–Stone（KS）法，以避免因验证集特征与建模集相近而令预测结果过优，当前也已成为深入研究光谱库时构建代表性样本集合的方法^［18-19］。本研究基于光谱前三个主成分，使用KS法以土柱为单位从总样本中按约5：1选出3167个样本组成土壤光谱库（记作SSL）；剩余588个样本做验证集，用于模型精度的检验（记作Test）。

在光谱分析时，使用主成分分析法（Principal Components Analysis，PCA）对预处理后的光谱数据进行降维，以实现少量不相关的变量来表征样本的主要特征。SSL和Test前三个主成分贡献率分别达55.2%和59.9%。

SSL和Test的SOC含量分布特征如表1所示，前者的均值同后者接近，但SSL的SOC含量范围跨度以及标准差值都显著高于Test。SSL包含的大量土壤变异与土壤样品的成土因素如母质、气候、地形和土地利用类型等密切相关。

1.3　局部建模集优化

光谱经预处理后，从SSL中优化局部建模集（Local）的过程分两步：光谱相似度计算和建模子集构建（图2）。

图 2. 局部建模集优化流程

Fig. 2. The scheme of optimizing local calibrations

下载图片 查看所有图片

本文考察了欧式距离（Euclidean Distances，ED）、马氏距离（Mahalanobis Distances，MD）和余弦距离（Cosine Distance，CD）三种经典方法来衡量Test样本与SSL样本之间的相似度。其中，ED和MD得到的距离值越小表示越相似，反之越相异；CD计算的是夹角余弦，值越接近0表示越相异，越接近1表示越相似^［20］。计算过程如图2左下所示，Test（n个样本）中每一个待测样本都与SSL（m个样本）的所有样本进行相似度计算，得到的n×m个距离值经排序之后组成距离矩阵（ED和MD按升序排，CD按降序）。随着距离矩阵行数的逐渐增加（从1到m），累积从SSL匹配到无重复的相似样本数（记作“库容”）也随之增加。

连续统去除也称去包络线，常用于提取连续变量的局部变化特征^［6］。如图2右下所示，先从距离矩阵提取库容曲线，再对该曲线进行连续统去除，然后通过设定可选阈值p和必选阈值q分别约束矩阵行数∆x和库容∆y，提取库容曲线上陡变之处的“特征点”。p值代表距离矩阵的行间隔，相邻行（p = 1）包含的样本信息存在显著冗余；q值过小则有更多的变化特征被选入，反之则更少。最后，用每个特征点的位置（第xi行）对应其从库中累积（前xi行）匹配到无重复样本，组成一系列的Local集。

矩阵第一行是从库中匹配到的最相似样本，样本数位于库容曲线的起点，所以该点位不会出现在连续统去除的曲线上。本文以距离矩阵第一行作为考察的起点，将其库容对应的样本所组成的建模集记作“Local₀”；以全库作为考察终点，记作“Global”。

1.4　建模方法及评价指标

本研究用经典的偏最小二乘回归（Partial Least-Squares Regression，PLSR）建立土壤MIR光谱与SOC含量间的定量关系模型，并用留一法交叉检验以避免模型校正过拟合^［21］。评估和比较模型的性能使用决定系数（Coefficient of Determination，R²），评估模型的误差使用均方根误差（Root Mean Square Errors，RMSE）。较高的R²和较低的RMSE表明模型较好。当预测精度差异不显著时，使用更少的库容则更佳。在R软件中，光谱数据处理使用prospectr包，PLSR使用pls包。

2　结果与讨论

2.1　光谱分析

图3（a）和3（b）分别是SSL与Test的平均光谱及95%置信区间。土壤MIR吸收光谱在4 000~2 000 cm^-1范围内表现剧烈，而在2 000~650 cm^-1呈现较多峰谷但幅度较小。在2 920 cm^-1和2 850 cm^-1与-CH₂伸缩振动有关，而在1 400 cm^-1附近的吸收由有机化合物引起^［22］。SSL和Test的平均光谱曲线无显著差异，但前者的置信区间更宽，如3 500 cm^-1和1 500~650 cm^-1附近。图3（c）和3（d）分别是SSL和Test光谱的前三个主成分特征向量曲线，共同解释了超55%的总体变异。PC1的贡献最大，且两者的PC1曲线波形相似，前述2 920 cm^-1、2 850 cm^-1和1 400 cm^-1附近与C有关的吸收也体现在该曲线上；PC2和PC3的曲线波形仅在部分波段表现出相似的波形。这表明，SSL与Test相似的特征信息主要保存在PC1中。

图 3. SSL（a，c）和Test（b，d）的平均光谱曲线及95%置信区间

Fig. 3. Average absorbance with 95% confidence intervals of SSL （a，c） and Test （b，d）

下载图片 查看所有图片

2.2　阈值设定与预测精度对比

随着行数的增加（图4（a）），ED、MD和CD的库容曲线均呈快速增长直至在约200行开始趋缓。放大观察可发现三条曲线并不完全重合，并随着各行样本量的加入，曲线斜率也发生相应变化，细节如图4（b）所示。在约200行之前，ED、MD和CD存在多处高于5的峰，在600~800行范围也存在少量的高峰，峰值在850行之后均低于5。三者连续统去除曲线在2 300或2 400行附近陆续开始平直，表明局部建模集的库容已达全库水平。

图 4. 距离矩阵所有行的库容（a）、连续统去除（b）及模型精度（c，d）的变化

Fig. 4. The capacity （a），continuum-removal （b） and model accuracy （c，d） for all rows of the distance matrix

下载图片 查看所有图片

图4（c）和4（d）展示了三种距离算法各自3167个模型预测同一Test集的精度。以图4（c）为例，三条曲线的R²值从第一行开始增长并在200行内波动且维持较高位，至200行附近出现陡降，随后在600 ~ 800行范围又出现了小幅跃升之后趋于平缓，最终收敛至全库的R²。表2显示，ED、MD和CD在第一行（Local₀）的库容均达300个左右。虽然MD的RMSE值高于全库（Global）之外，但其R²明显高于ED、CD以及Global。总体而言，三者Local₀的各指标均优于Global，因此，从前200行快速定位出精度同时优于Local₀和Global的行位置是关键。

图5更详细地展现了前200行的曲线变化情况。经连续统去除后（图5（b）），ED、MD和CD的第一个峰均高于5，结合图4（b）的全局变化，本研究将图2中的必选阈值q设置为5。连续统去除后所有大于5的“特征点”以相应颜色的圆点标于图5各曲线上。而预测精度曲线（图5（c）、5（d））显示，相邻行特征点的精度几乎一致。因此，本文建议将可选阈值p设置为1，即当p > 1时非相邻的特征点可被保留。本研究在局部建模集的优化过程中对p和q的阈值设定仅供参考，在今后的研究中其普适性还需更多方的实验加以印证。

图 5. 距离矩阵前200行的库容（a）、连续统去除（b）和模型预测精度（c，d）的变化

Fig. 5. The capacity （a），continuum-removal （b） and model accuracy （c，d） for the first 200 rows of the distance matrix

下载图片 查看所有图片

以图5（c）为例，ED、MD、CD的R²曲线总体在前50行和100 ~ 150行区间范围存在较高值。图5（a）显示，50行附近三者的库容已增长至约75%，150行时则已接近全库。因此，本研究继续从前50行寻找合适的局部建模集。图5（c）、5（d）显示，ED、MD和CD由第一特征点组成建模集（Local₁），其精度均优于各自的Local₀，即比Local₀的R²更高且RMSE更低。结合表2，三者Local₁的精度也均优于Global。

MD在Local₁的精度优于ED和MD，且430个样本的库容相较于ED的926个和CD的1038个明显更少。而ED和CD此处的精度差异不显著，两点在图5（c）、5（d）中几乎贴在一起。因此，MD更适用于对土壤MIR光谱进行相似度的衡量，这与前人对土壤vis–NIR光谱进行的相关研究结论不同，他们认为CD更佳^{［20，23-24］}。这说明MIR与vis–NIR不仅数据表征不同，还存在方法上的适用性差别。

值得注意的是，三者第一特征点的位置并非对应前5行R²最高点和RMSE最低点，这可能是由于该数据包含了大量的土壤变异信息，前三个主成分累计贡献也未超过80%所致。即便如此，寻找并确定Local₁的过程简单、快速，且其精度优于首行样本集和全库，已达到本研究的目的。

2.3　主成分空间分布

以MD为例，图6中的红色点是其Local₀、前4个特征点（Local₁ ~ Local₄）和Global的库容在主成分空间的分布情况，蓝色点代表待测样本集Test。从Local₀开始，局部建模集主成分空间的覆盖范围（以阴影表示）随着库容的增大而随之扩大并发生形变。Local₀的库容占8%，其阴影范围接近Test但未完全覆盖。结合表2说明，尽管Local₀是从SSL中匹配出的与Test最相似的样本，但由于Local₀未能完全表征Test，预测精度虽然较高但仍未达到最佳水平，这与彭青青等的研究吻合^［20］。

图 6. 基于MD法组建的五个Local集（红点）投影到Test集（蓝点）主成分空间的分布

Fig. 6. The PCA score plot from five Locals （red dots） projected onto Test （blue dots） by MD

下载图片 查看所有图片

占库14%的Local₁的阴影范围与Local₀接近，但前者因样本的加入而使其范围更大，如PC12的第三象限、PC13的第二象限、PC23的第一象限都或多或少超越了Test，从而使Local₁的精度优于Local₀。随着库容的逐渐增大伴随着阴影范围显著超越Test，Local₁之后的局部建模集相应的精度也开始下降，这种现象表明，越来越多与Test相异的样本被加入到了Local集合中^［25］。因此，Local₁是当前优化框架下得到的较为理想的局部建模集。

3　结论

本研究针对土壤MIR光谱库提出了一种快速构建最优局部建模集的方法，包含两步：光谱相似度计算和建模子集构建。通过比较ED、MD和CD三种距离算法来衡量光谱相似度并得到距离矩阵，使用连续统去除法来提取库容曲线上的特征点，并考察了每种距离算法各自3167个模型对同一个待测集合的预测和检验。结果表明：

（1） ED、MD和CD三种方法的第一特征点均可得到较佳的预测精度；

（2） 对库容曲线进行连续统去除是提取特征点的有效手段；

（3） 待测样本通过MD法从土壤MIR光谱库构建的最优局部建模集预测精度较ED和CD更优，所用库样本数最少。