华南地区土壤全钾含量高光谱反演模型研究

彭一平，刘振华，王 璐*，赵 理，胡月明,5

(1. 华南农业大学，广东 广州 510642；2. 国土资源部建设用地再开发重点实验室，广东 广州 510642；3. 广东省土地信息工程技术研究中心，广东 广州 510642；4. 广东省土地利用与整治重点实验室，广东 广州 510642；5. 青海大学农牧学院，青海 西宁 810016)

【研究意义】土壤是农业生产的重要载体，是陆地生态的关键组成部分。土壤钾是土壤中含量最高的营养元素,也是植物生长必需的营养元素之一。在农业生产的过程中，常因土壤钾供应不足而影响作物产量[1-2]。因此，借助现代的技术对土壤全钾含量进行快速的监测尤为重要。【前人研究进展】土壤全钾传统监测方法是采用野外采样调查，结合实验室分析(火焰分光光度法，火焰原子吸收法、流动分析注射仪法和极谱法等方法)对土样采用进行全钾测定[3],最后进行空间插值从而得到区域土壤全钾含量的分布。这些方法虽然测量土壤全钾含量精度较高，但是过程繁琐，需要耗费大量的人力、物力和财力，并且往往具有滞后性。而遥感技术具有快速、大范围动态监测的优势，已有研究表明，土壤钾元素主要以矿物态形式存在，而原生矿物对土壤光谱具有一定的影响，因此利用光谱分析技术预测土壤钾含量，理论上具有一定的可行性[4]。近年来，国内外已有学者尝试利用高光谱遥感技术反演土壤钾含量。Viscarra Rosse等针对不同波长，结合偏最小二乘回归法对多种土壤属性的含量进行反演，结果表明交换性钾在近红外光谱反演效果最佳，其调整R2为0.47，均方根误差(RMSE)为1.84，但其反演效果远不如有机碳、pH等土壤属性[5]。胡芳对光谱数据的4种光谱指标(反射率、反射率倒数的对数、反射率一阶微分和波段深度)采用偏最小二乘回归(PLSR)方法建立反演模型，并对其进行了精度验证。结果表明，波段深度为最优的反演指标，其决定系数R2>0.85，RMSE<0.1[6]。Debaene等采用PLS模型对多种土壤属性(土壤速效钾、速效磷、SOC等)进行估算，结果表明土壤速效钾和速效磷的反演效果不理想，其决定系数R2小于0.5[7]。胡国田等对美国密苏里州8种类型土壤进行正交信号校正(DOSC)处理，以未进行DOSC处理数据为参照，运用偏最小二乘回归(PLSR)分析方法建立土壤磷和钾的反演模型。结果表明，DOSC处理后得到的反演模型，其精度显著提高[8]。上述研究表明利用光谱分析技术估算土壤钾含量的测定具有一定的可行性和准确性，但是由于土壤钾含量在不同的地区，其组成成分和含量差异很大，导致其预测效果随土壤类型、测试环境、颗粒大小、土壤性质等变化，很不稳定，有必要针对不同地区对其进行进一步的研究[4]。【本研究切入点】目前在利用光谱反演土壤全钾含量的研究中，其建立的估算模型通常是不具备普适性。而华南地区作为主要我国农业产区之一，尚未有研究利用高光谱技术构建适用于该地区的土壤全钾预测模型。【拟解决的关键问题】本文的研究目标是分析实测土壤光谱数据与土壤全钾含量之间的关系，构建亚热带典型区域广东省的土壤全钾光谱反演模型，探讨利用高光谱技术进行区域土壤全钾含量快速监测的可行性，为农业有效管理提供新的思路。

1 材料与方法

1.1 研究区概况

广东地处中国大陆最南部，陆地范围位于北纬20°09′～25°31′，东经109°45′～117°20′。自东至西依次与福建省、江西省、湖南省、广西壮族自治区接壤；毗邻香港、澳门特别行政区；西南端隔琼州海峡与海南省相望。根据2016年度土地变更调查，全省陆地面积17.97万km2，约占全国陆地面积的1.87 %。广东地表形态主要分为山地、丘陵、平原、台地4种类型，地形总体呈北高南低之势。山地丘陵居多，草地分布面积较小，平原分为三角洲平原和河谷冲积平原2种类型，约占土地总面积的23.4 %。广东省属于东亚季风区，全省平均日照时数为1745.8 h，年平均气温22.3 ℃，年平均降水量为1777 mm，是中国光、热和水资源最丰富的地区之一。

1.2 数据来源及预处理

1.2.1 土壤样本采集及处理 土壤样本的采集时间为2017年5月,在广东省境内共采集了75个土壤样本，每个样本重约300 g，所有采样都是在土壤表层(0～20 cm)完成，采样分布点如图1所示。土壤样本经过剔除杂质，风干后，进行碾磨并过0.2 mm土壤筛，最后将每个样品均分为2份，分别用于测定土壤光谱反射率与土壤全钾含量。

图1 研究区土壤样点分布Fig.1 Distribution of soil samples in the study area

表1 土壤全钾含量描述性统计

为了测定土壤全钾含量，每个样本分别称取土壤0.2 g，用H2SO4-HNO3-KMnO4法消解，测定方法为三酸消化-ICP测定法。测定的土壤全钾含量描述性统计如表1所示。其变异系数为72.80 %，表明研究区的土壤全钾含量为中等变异性。在随后的数据处理中将75个样本数据按照4∶1的比例随机分为建模样本与验证样本。

1.2.2 土壤光谱放射率测定及处理 本研究使用荷兰Avantes公司生产的AvaField便携式光谱仪对土壤光谱数据进行测定，波段范围为340.316～2511.179 nm，采样间隔为0.6 nm，测定光源为50 W的卤素灯，视场角为10°，将探头垂直接触样本进行光谱数据采集，每个土样的光谱测定重复5次，每次自动采集10条数据，利用AvaReader软件剔除异常曲线后，取光谱反射率的算数平均值作为土样的光谱反射率。考虑到波段测试范围两端为光谱不稳定区域，因此在本文中不予考虑。

1.3 研究方法

1.3.1 光谱数据分析 由于光谱仪波段之间对能量响应上的差别，会使光谱曲线上存在许多突兀的“毛刺”噪声，为了降低噪声的影响，选用Savitzky-Golay平滑方法对光谱数据进行平滑处理[6]。

由图2可知，研究区域的土壤反射光谱曲线波形大致相同，呈上凸的抛物线型，在400～1100 nm光谱区间，随着波长的增加，反射光谱曲线单调上

图2 土壤光谱反射率Fig.2 Spectrum of soil samples

升，在1100 nm附近达到峰值后，光谱曲线变化较为平缓。在1400、1900和2200 nm附近存在3个明显的吸收峰，深度略有差别，符合土壤光谱的曲线特征。

为了突出光谱曲线的吸收和反射特征，更好的筛选特征波段，本研究对光谱数据进行了连续统去除处理。由图3可以看出，光谱曲线吸收带主要集中4个波段附近(500、1000、1400、1900和2200 nm)。根据前人研究可知，500 nm附近是主要是有机质的光谱吸收区，1000 nm附近的吸收主要是矿物引起的，1400、1900和2200 nm附近的吸收主要是水分子振动引起的[9]。

1.3.2 构建光谱指标 通过对光谱进行适当变换，可消除或减弱背景噪声的影响，突出光谱曲线的吸收和反射特征，减小土壤表面散射和吸收引起的信号强度的改变，提高光谱敏感度，从而提高模型的估算能力和稳定性[10]。因此，本研究以平滑处理后的光谱曲线作为土壤样品的光谱变换的原始光谱反射率，对其进行3种光谱变换处理(一阶微分、二阶微分以及倒数对数)，增强反射率与土壤全钾含量的相关性，更有效的筛选波段。

1.3.3 模型建立与验证 利用皮尔逊相关分析进行特征指标的筛选，以经过筛选后的光谱数据作为函数的自变量，实验室测得的土壤全钾含量作为因变量，对土壤全钾含量进行曲线拟合。利用决定系数(R2)和均方根误差(RMSE)验证模型效果，其中决定系数(R2)用于衡量模型的稳定性，决定系数越大，代表模型越稳定；模型精度是由均方根误差(RMSE)来反映，均方根误差越小，预测能力越好[12]，R2与RMSE的计算公式分别如下式(1)、(2)所示：

图3 土壤连续统去除光谱反射率Fig.3 Continuum removal of soil spectrum

表2 光谱指标描述

a) 光滑处理， b) 一阶微分，c)二阶微分 和d)倒数对数a) REF; b) FDR; c)SDR and d)RL图4 土壤光谱指标与土壤全钾含量的相关系数Fig.4 Correlation coefficient between soil total potassium content and spectral indices in soil

(1)

(2)

表3 光谱指标与土壤全钾含量的特征波段及相关系数

注：*表示相关性在0.05层上显著(双尾)，**表示相关性在0.01层上显著(双尾)。

Note：* Correlation is significant atP<0.05 level; **Correlation is significant atP<0.01 level.

2 结果与分析

2.1 土壤全钾含量与光谱指标的相关性分析

利用皮尔逊相关分析分别将反射率光谱平滑(REF)、一阶微分(FDR)、二阶微分(SDR)、倒数对数(RL)4种光谱指标与土壤全钾含量进行相关性分析，结果如图4所示，可知FDR与土壤全钾含量的相关性明显高于其他处理，其特征波段主要在500和1000 nm附近。

从图4a～c可看出，4种光谱指标与土壤全钾含量的相关性大小为：FDR>SDR>RL>REF，其对应的最高相关系数为0.50，0.25，0.22，0.20，由此表明反射率一阶微分(FDR)与土壤全钾含量的相关性最好，为最优光谱指标。由于不同光谱波段的信息可以互相补充，为了提高土壤全钾与光谱指标之间的相关性，本文对光谱数据进行了波段组合，由表3可知，与未进行波段组合处理相比，相关系数有了显著的提高，其中反射率一阶微分的波段组合相关系数达到0.634，其最优的波段组合为R442.459nm+R1739.447nm+R965.768nm。

2.2 模型构建

将筛选所得的特征波段和土壤全钾含量分别作为自变量x与因变量y，采用6种常用的回归模型，除去部分因自变量为负值而无法利用对数函数和幂函数拟合，其余的分别建立反射率光谱平滑(REF)、一阶微分(FDR)、二阶微分(SDR)、倒数对数(RL)最优波段组合的土壤全钾反演模型，如表4所示。

表4 反演模型对比分析

表5 土壤全钾含量预测值与实测值的比较

从表4可以看出，不同光谱指标的预测效果差异较大。在建模阶段，R2介于0.12～0.65，RMSE介于4.850～7.978 g/kg；验证阶段，R2介于0.02～0.63；RMSE介于4.347～9.449 g/kg。通过比较分析，基于反射率一阶微分(FDR)波段组合构建的预测模型反演效果明显优于其他3种。进一步比较反射率一阶微分(FDR)波段组合所建立预测模型，其中，一元二次函数模型y=2E+06x2+11328x+16.372拟合效果最好，R2达到0.64，RMSE为4.850 g/kg，预测结果如图5所示。

通过计算验证数据集(n=15)的精度参数(R2和RMSE)，进一步验证最优模型的预测效果。由表5可知，验证样本的预测平均值为10.30 g/kg，实测平均值11.95 g/kg，其对应的标准差分别为4.81和7.29。验证样本R2为0.61，RMSE为6.87 g/kg。图6为土壤全钾含量实测值与预测值的比较，散点分布接近于1∶1线，可知此模型反演结果可靠。

3 讨 论

本研究所得土壤全钾含量与光谱变量最大相关系数与前人的研究有较大的差异[3-4,13-14]，主要原因是针对不同的地区，土壤类型差异大，其土壤属性与反射率差异较大。但是由于经费、时间等限制，本文仅选取广东省作为华南地区典型区域，进行土壤全钾含量的高光谱反演研究，今后还需要在其他地区增加采样点，以提高模型的可靠性；此外，本研究是利用采样点实测光谱数据与土壤全钾含量数据，进而建立样点尺度上的土壤全钾预测模型，今后可考虑利用高光谱遥感影像，对土壤全钾含量进行区域尺度上的反演。

图5 土壤全钾含量最优反演模型Fig.5 Optimal inversion model of soil total potassium content

4 结 论

本文对土壤光谱反射率数据进行Savitzky-Golay平滑(REF)、一阶微分(FDR)、二阶微分(SDR)、倒数对数(RL)4种光谱变换处理后，利用皮尔逊相关性分析，筛选出各光谱变换处理的特征波段。为了充分利用光谱信息，对特征波段进行了波段组合，结果表明，与单波段相比，波段组合能显著提高光谱变量和土壤全钾含量之间的相关性，其中反射率一阶微分的波段组合R(442.459nm)+R(1739.447nm)+R(965.768nm)相关性最高，其相关系数达到0.634。本研究建立6种回归模型对土壤全钾含量进行高光谱反演，通过对比分析，最佳反演模型是以反射率一阶微分(FDR)的波段组合R(442.459nm)+R(1739.447nm)+R(965.768nm)为光谱指标所建立的一元二次函数模型y=2E+06x2+11328x+16.372，建模时决定系数R2为0.64，RMSE为4.850 g/kg，验证时决定系数R2为0.61，RMSE为6.874 g/kg。研究结果为土壤养分含量的高精度预测提供了可靠的参考依据。

图6 土壤全钾实测值与预测值Fig.6 The predicted values against the measured contents