基于高光谱数据的盐渍化土壤含盐量预测模型研究

关键词：高光谱；土壤含盐量；盐渍化土壤；预测模型；多元逐步线性回归

0 引言

土壤盐渍化严重制约农业生产和发展，已成为全球普遍存在的环境问题，是一种常见的土地退化过程[1]。土壤盐渍化是一个复杂的过程，分析盐渍化土壤在时间与空间上的动态变化，能够更好地掌握当地盐渍化的程度与状态，为农业生产、生态环境保护提供理论依据[2-3]。近年来研究人员研究了许多精确和高效的土壤盐分定量分析方法，这些方法对土壤盐渍化的改善和治理、土地资源的保护、农业生产效率的提高作用显著，但对于大范围区域很难做到及时、准确地监测[4-5]。因此，广泛开展盐化土壤动态监测、土壤含盐量精确估算等研究备受关注[6-7]。

在遥感技术日新月异的今天，地物的细微之处在光谱数据中都能得到准确地反映[8]。王飞等[9]对新疆维吾尔自治区天山典型绿洲进行深入研究，利用线性与非线性模型对土壤盐分敏感的原始波段、盐分指数、植被指数等指标进行对比分析，发现SI-T、TGDVI、EEVI、BAND2和BAND5等指标对盐分特别敏感。WHITNEYK等[10]通过对土壤含盐量时空变化的遥感光谱指数和冠层温度进行动态监测，对土壤大尺度盐分时空序列进行数据分析。王丽娜等[11]利用主要成分回归分析方法，基于对盐渍土壤光谱曲线变化规律的研究进行分析，确定了黄河三角洲盐分敏感波段1490～1608和1900～1950nm两个波长范围，并通过估算得到了区域土壤含盐量。在土壤修复过程中，朱赟等[12]通过对相关系数极值和不同相关系数范围分析，从8种光谱数据集中筛选最优波段，运用偏最小二乘法，建立了光谱反演模型。利用野外敏感波段，张同瑞等[13]建立了光谱指数，进一步构建18个模型并验证选定，在黄河三角洲地区盐分条件下，以土壤调节植被指数的线性模型为基础，确定最优的反演分析模型。蔡东全[14]对土壤盐化程度不同的光谱特征规律进行了分析，利用相关分析的方法获得特征波段，以多元线性回归分析方法建立高光谱预测模型，通过波谱角分类法得到了土壤盐化信息，建立了土壤波谱库。黄帅等[15]通过高光谱指数与土壤含盐量的相关性分析，以最小二乘回归和多元渐线性回归的方法，筛选出土壤含盐量的光谱特征波段，建立土壤盐分动态监测模型。蒲智等[16]确定了最优光谱指数和定量预测模型，分析了土壤含盐量与5种高光谱指数的定量关系。亚森江·喀哈尔[17]利用室内高光谱数据建立了渭干河−库车河三角绿洲的PLSR和SVM土壤含盐量反演模型，R2分别为0.762和0.637。为提升模型的反演精度，张晓光等[18]在原始光谱反射率的基础上进行整数阶微分变换，结果表明，一阶和二阶微分变换比原始光谱反射率对土壤含盐量更加敏感，建立的模型精度更高。本研究以内蒙古自治区（简称内蒙古）科尔沁右翼中旗巴彦淖尔苏木为研究区，根据室内分析结果，对盐渍化土壤进行分析，依据高光谱数据，将光谱数据进行多种数学变换，并与土壤含盐量进行相关分析，从而对响应波段进行优选提取。采用地统计分析方法构建盐渍化土壤遥感监测模型，模拟测定含盐数量与光谱数据之间的数学关系，为大规模定量监测盐渍化土壤提供科学依据。

1 材料与方法

1.1 研究区概况

巴彦淖尔苏木隶属于内蒙古兴安盟科尔沁右翼中旗，地处东经121°46′27.5″～122°12′54.8″、北纬44°30′22.3″～44°49′2.7″，海拔高度184m，面积527.1km2。巴彦淖尔苏木属温带大陆性季风气候，春季少雨多风、夏季多雨高温、秋季风大干燥。巴彦淖尔苏木地处科尔沁腹地，地形属沙坨平原地区，风沙土和草甸为主。巴彦淖尔苏木耕地面积1万hm2，由于生态破坏、土地退化、原生与次生盐碱化并存，根据全国第四次沙化监测，巴彦淖尔苏木盐碱地面积1.429万hm2，其中林地面积0.016万hm2、草原面积1.413万hm2，盐碱地面积大、程度重，严重制约当地经济、生态和农牧业的发展[19]。

1.2 数据采集与处理

1.2.1 土壤样品采集与处理

2020年6月中旬在研究区内开展野外调查，按照梅花桩取样法（重复采集5次，然后对土壤进行混匀）对研究区表层（0～10cm）土壤样本进行采集，共采集土壤样本100个。样品在实验室内自然风干，充分研磨后过2mm筛。每个样品分2份，用于土壤光谱检测1份，用于水溶性盐总量指标测定1份。土壤水溶性盐用水浴蒸干法测定总量：吸取一定量的土壤浸出液放在瓷蒸发皿上，在水浴上蒸干，用过氧化氢氧化有机质，在105～110°C烘箱中烘干，称质量。

1.2.2 土壤样本光谱测定

光谱检测采用FieldSpec4光谱仪（analyticalspectraldevices），频段范围350～2500nm，采样间隔1nm，输出波段数2151个。在暗室内进行光谱测定，先校正后进行光谱测定。在暗室内对土壤样本进行光谱检测，检测之前进行校正。将处理好的土样装满盛样皿，表面刮平后开始检测，每个样品检测10次，以算术平均值作为实际光谱反射率。

1.3 数据处理

1.3.1 光谱数据平滑去噪

运用Origin8.6软件中的Savitizky-Golay平滑法进行9点加权平均，对测得的光谱反射率进行平滑处理。光谱曲线中噪声较大的波段范围经过平滑处理后被删除，400～2500nm波段的光谱信息得以保留。

1.3.2 实测光谱数据变换

由于原始光谱对噪声的敏感性不强，原始光谱经过数学变换处理后，再与土壤含盐量构建模型，对土壤含盐量进行估算[20]。为了更好地反映土壤含盐量与光谱反射率的相关关系，对原始光谱反射率R分别进行一阶导数R′、二阶导数R″、平方根、平方根一阶导数（）′、平方根二阶导数（）″、对数logR′、对数一阶导数（logR）′、对数二阶导数（logR）″、对数倒数[log（1/R）]、对数倒数一阶导数[log（1/R）]′、对数倒数二阶导数[log（1/R）]″、倒数（1/R）、倒数一阶导数（1/R）′和倒数二阶导数（1/R）″14种常见光谱变换，以选择最佳的光谱处理方法。光谱微分处理计算公式为

1.3.3 相关性分析

对样本光谱数据与土壤含盐量的相关系数进行逐波段测算，对特征波段进行提取。计算相关系数公式为

1.4 建模与检验

1.4.1 模型建立方法

通过原始光谱反射率及其数学变换形式与土壤含盐量进行相关性分析，提取特征波段。供试土壤样本数为100个，以供试土壤样本的70%用于模型建立，30%用于模型的检验，采用多元逐步线性回归方法建立土壤含盐量高光谱预测模型。

1.4.2 模型检验

利用回归决定系数（R2）和均方根误差ERMSE对建模集和验证集进行检验。R2越大、ERMSE越小，表明模型稳定性和预测效果越好。计算公式为

2 结果与分析

2.1 土壤盐渍化等级划分

测得的土壤样本指标值统计特征如表1所示。样本土壤含盐量最大值55.34g/kg、最小值0.52g/kg、平均值13.10g/kg及变异系数97.25%，属强变异，表明研究区土壤含盐量分布很不均匀。根据研究区实测点的情况，参考《中国盐渍土》中土壤分级方法，确定研究区土壤盐渍化分级标准，如表2所示[21]。

2.2 土壤光谱特征分析

将100个土壤样本的反射率依照表2进行分类，在不同盐渍化土壤样本中，取其平均值作为光谱反射率曲线，如图1所示。5条曲线的变化趋势基本一致，供试土壤含盐量越高，光谱的反射率也越大。光谱反射率在350～1300nm范围内呈上升趋势，增加幅度较大，在1400～1850nm范围内反射率不断升高，但增加幅度较缓，1900～2150nm范围内光谱反射率增加且增加速度很快，从2150～2500nm范围内，光谱反射率先逐渐减小，后又快速增加再急剧下降，在接近1400、1900和2200nm处出现3个明显的水分吸收谷。

2.3 最佳波段选取

由于供试土壤含盐量变化较大，对最佳波段的选取会造成一定影响，根据土壤盐渍化分级标准将参试土壤分为两类：非盐土（土壤含盐量≤7g/kg）和盐土（土壤含盐量>7g/kg），计算土壤含盐量与15种数学变换的光谱数据的相关系数，如图2所示。非盐土光谱反射率与土壤含盐量的相关系数最大0.3646，盐土的原始光谱反射率与土壤含盐量的相关系数最大0.4545，经过数学变换后，与土壤含盐量的正负相关性明显增强。非盐土的对数一阶微分处理反射率与土壤含盐量在1490nm处相关性最高，相关系数最−0.5898，盐土的对数一阶微分反射率与土壤含盐量在727nm处相关性最高，相关系数−0.5591。数据中的细微信息经数学变换后得以放大，可以据此筛选出提取pRpR土壤含盐量的特征波段。根据图2中相关分析，R′、R″、（）′、（）″、（logR）′、（logR）″、[log（1/R）]′、[log（1/R）]″、（1/R）′和（1/R）″10组变换下显著相关的波段数较多，而其他变换没有或极少存在与土壤含盐量显著相关的特征波段。因此从以上10种变换中选择特征波段，如表3和表4所示。

2.4 土壤含盐量估测模型构建

根据相关分析结果得到的特征光谱波段，建立土壤含盐量的多元逐步线性回归模型。运用DPS数据处理软件将非盐土47个样本、盐土53个样本中的70%（非盐土32个样本，盐土36个样本）用于建模，30%（非盐土15个样本，盐土17个样本）用于模型验证。以实测土壤含盐量作为因变量，以特征波段的光谱反射率作为自变量，在保证R2尽可能大的前提下，建立土壤含盐量的多元逐步线性回归模型，如表3和表4所示。

利用多元逐步线性回归分析方法建立的盐渍化土壤含盐量预测模型精度较高，非盐土用于建模的R2平均值0.7306、最大值0.7685、最小值0.6514，ERMSE平均值0.97、最大值1.10、最小值0.90；盐土用于建模的R2平均值0.5557、最大值0.7682、最小值0.3466，ERMSE平均值7.76、最大值9.50、最小值5.66。根据R2越大、ERMSE越小，模型越稳定的原则，非盐土的原始二阶方程建模的R2=0.7685、ERMSE=0.90，模型稳定性最好；盐土的对数二阶方程建模的R2=0.7682、ERMSE=5.66，模型稳定性最好。

2.5 模型精度检验

模型的验证主要考虑稳定性和预判能力两个方面。测试模型的预测能力，以未参加建模的土壤样本中土壤含盐量为准真值。由表3和表4可知，非盐土检验模型检测的R2最大值0.7292、最小值0.4059，ERMSE最大值1.91、最小值0.95；非盐土的原始二阶方程检测的R2=0.7292、ERMSE=0.95，模型稳定性最好。盐土检验模型检测的R2最大值0.8718、最小值0.3944，ERMSE最大值9.77、最小值5.80，盐土的对数二阶方程检测的R2=0.7682、ERMSE=5.66，模型稳定性最好。

由图3可知，建模样本基本上聚集在对角线附近，非盐土原始光谱反射率二阶微分预测模型实测值与预测值间的样本拟合系数R2=0.8246，盐土原始光谱反射率对数二阶微分预测模型实测值与预测值之间的样本拟合系数R2=0.7666，预测值和实测值具有较好的一致性。本研究建立的土壤盐渍化遥感监测模型是有效的，经检验合格，可在一定程度上用于预测土壤含盐量。

3 讨论

在构建模型时发现，当对全部供试土壤反射率及其数学变换形式进行整体构建模型时，回归决定系数（R2）整体不高，模型稳定性较差，预测效果不好，当将供试土壤按照土壤含盐量进行分类，再对其反射率及其数学变换形式进行建模时，R2值整体升高，模型的稳定性和预测效果增强。这是由于供试土壤实测土壤含盐量最大值55.34g/kg、最小值0.52g/kg、平均值13.10g/kg及变异系数97.25%，研究区域内土壤含盐量分布极不均匀，导致在特征波段选取上存在差别。

本研究是对采自内蒙古自治区科尔沁右翼中旗巴彦淖尔苏木典型盐渍化土壤进行光谱反射率测定，并对其进行均方根、对数、对数倒数、倒数及其一阶、二阶微分等数学变换，利用多元逐步线性回归分析方法，建立盐渍化土壤含盐量预测模型，并对预测模型进行对比。研究表明，非盐土的原始光谱对数一阶微分预测模型最优，盐土的对数一阶微分预测模型最优，这也说明选取不同的数学变换形式在模型预测上是不同的。因此本研究筛选的光谱变换形态所建立的预测模型，对于其他地区盐土的适用性尚待考证。

4 结束语

（1）从形态上看，不同盐渍化程度的土壤光谱特征曲线基本一致；土壤含盐量越高的土壤样本，其光谱反射率就越高。

（2）光谱反射率经数学变换后与土壤含盐量的正负相关性增强，特别是经过一阶和二阶微分变换后，相关性显著增强：非盐土的原始光谱对数一阶微分反射率与土壤含盐量相关性最高，在1490nm处相关系数最大为−0.5898，在727nm处盐土原始光谱对数一阶微分反射率最大为−0.5591。

（3）利用多元逐步线性回归建立的预测模型，非盐土以原始光谱的二阶微分模型为最优，检测R2=0.7292，盐土以对数二阶微分模型为最优，检测R2=0.8718，可快速测定盐碱土土壤含盐量。