基于Landsat 5影像的青海省大通县土壤表层全氮空间格局反演

雷浩川 ， 刘尊方， 于晓晶， 曹金莲， 吴笑天

(1.青海省基础地理信息中心，青海 西宁 810016； 2.青海大学地质工程系，青海 西宁 810016)

土壤是人类生产生活中重要的自然资源，是人类生存环境的重要载体。土壤作为土地资源的主要物质基础，其固有的不可再生性决定了承载容量的有限性[1,2]。在数字化农业快速发展的时代，准确、快速、动态按需获取土壤信息是现代精准农业的保证。土壤全氮含量是土壤氮素供应状况的重要指标，也是衡量土壤肥力的重要因素之一。氮素是植被生长过程中所需较多的营养元素[3]。随着现代农业的发展，如何在有限的时间里获取所需的土壤信息，进而及时评价土地肥力，指导农业生产的测土配方科学施肥，已是研究人员面临的重大课题。如今，实时获取土壤全氮含量，利用遥感影像进行大范围土壤全氮反演是现代精准农业的迫切要求。从 20 世纪 70年代开始，土壤光谱与氮素的关系已经是国内外许多学者的研究热点[4]。徐丽华等[5]构建PLSR模型反演土壤全氮含量，结果显示土壤全氮含量实测值与预测值的相关系数r为0.672，预测效果较好。李燕丽等[6]研究遥感因子与土壤全氮含量变化的关系，研究结果显示，土壤全氮含量预测的RMSE为0.35，表明土壤全氮含量反演模型的预测精度较高；王世东等[7]建立了PLSR、BPNN、RF等模型，并将PLSR与其他两个模型相结合共同反演土壤全氮含量，结果表明，PLSR-BPNN模型的回归决定系数R2可以达到0.92，说明模型具有较好的反演能力；祁亚琴等[8]构建了土壤全氮含量的指数函数估算模型，该模型的预测效果最好，其回归决定系数R2最高，为0.7982；杨越超等[9]构建PLSR模型和BPNN模型反演了黑土土壤全氮含量，BPNN模型的土壤全氮含量预测精度比PLSR模型的预测精度高6.5%；张强等[10]基于NDI估算土壤全氮含量，其中预测效果最好的模型是指数函数(回归决定系数R2为0.7982)；袁石林等[11]采用的最小二乘支持向量机(LS-SVM)的耦合模型对土壤氮素含量的反演精度较高；张娟娟等[12]研究了我国中、东部地区5种土壤，结果表明，经过Norris滤波平滑之后，根据反射率的一阶导数光谱构建的模型预测效果最好；邱壑[13]研究反演了SOM、速效氮等养分含量，通过相关性分析提取相关波段，并根据实测值与预测值的偏差选出最优的反演模型；肖云飞[14]的研究表明，随机森林模型预测土壤全氮含量的精度高于PLSR模型，PLSR模型建模具有一定的局限性。

综上，基于卫星遥感影像的土壤全氮含量预测反演仍有很大的研究空间。本文以Landsat 5 TM多光谱遥感影像为数据源，分析了青海省西宁市大通回族土族自治县的土壤表层全氮空间分布格局。利用PLSR模型和BP神经网络模型研究该地区的土壤全氮含量，为土壤碳氮循环研究、土壤质量评价、粮食估产等工作提供数据支撑和理论参考。

1 材料与方法

1.1研究区概况青海省西宁市大通回族土族自治县位于青海省东部河湟谷地(100°51′～101°56′E，36°43′～37°23′N)，祁连山南麓，湟水河上游北川河流域，是青藏高原和黄土高原的过渡地带。海拔2 280～4 622 m，地势西北高、东南低。研究区属于高原大陆性气候，年平均气温4.9 ℃，最值气温分别达到35.6 ℃和-26.1 ℃，8月份降水量最多。研究区行政区划图及采样点分布图如图1所示。

图1 研究区行政区划图及采样点分布图Fig.1 Administrative division and sampling point distribution of the study area

1.2样本采集在研究区的农用地块以0～20 cm 深度采集土壤样品73个，各约2 kg。在实验室将样品自然风干、磨细。每个样品分为2份，分别用于光谱测试和土壤成分化学分析。土壤全氮质量分数的测定采用半微量凯氏定氮法[15]。53个样品用于建立模型，剩余的20个样品用于验证模型。53个建模样品的全氮质量分数最大值为3.375 g/kg，最小值为0.639 g/kg，平均值为 1.857 g/kg，标准差为0.578 g/kg，变异系数为18%。

1.3影像获取及预处理本研究所用的技术流程如图2所示，Landsat 5 TM遥感影像数据从谷歌地球引擎(GEE)遥感云平台快速获取。观测参数及波段特征如表1所示。影像成像时间分别是2008年9月11日和2008年9月20日，与土壤样本的采集时间基本一致。9月左右农作物已经收割完毕，土壤表面没有大量植被和作物遮挡，获取土壤光谱信息较为方便。通过实地采样养分化验数据和影像各个波段的相关性分析，得到土壤全氮含量的最佳反演波段。

图2 技术路线图Fig.2 Technology roadmap

表1 Landsat 5 TM观测参数及波段特征Tab.1 Observation parameters of Landsat 5 TM and band characteristics

几何校正、辐射校正(辐射定标、大气校正等)、影像去云、影像拼接和影像裁剪等是遥感影像预处理的主要过程。经过影像预处理，原始影像的数字量化值(Digital Number，DN)转化为地表真实的反射率(Surface Reflectance，SR)。谷歌地球引擎(GEE)遥感云平台可以提供大气表观反射率TOA影像集和地表反射率SR影像集，为处理数据节省大量时间。影像预处理阶段从GEE平台使用函数ee.ImageCollection调用Landsat 5地表反射率影像(LANDSAT/L5/C01/T1_SR)，对于调用的影像只进行影像去云操作、影像拼接和影像裁剪的预处理。依据前人的研究结论可知，仅仅通过原始波段的反射率提取特征波段有一定的困难，但是通过处理原始反射率可以有效提取全氮的特征波段(光谱预处理包括反射率的对数(logR)、倒数(1/R)、倒数的对数log(1/R)和对数的倒数1/logR等处理。)

土壤全氮含量与波段反射率及变换形式间的关系如表2所示，可以看出土壤全氮含量与波段原始反射率的相关性并不显著，最大的负相关系数-0.534出现在B2波段，并且可见光波段的反射率与土壤全氮含量的相关性高于其他波段。不过经过数学变换之后，原始反射率与土壤全氮含量的相关性明显有所提升，尤其是经过倒数处理的反射率与土壤全氮含量相关性最高，最大相关系数r达到了0.584。通过以上分析可知，与土壤全氮含量最敏感的是可见光波段(B1、B2和B3)。

表2 土壤全氮含量与波段反射率及变换形式间的关系Tab.2 Relationship between soil total nitrogen and band reflectance and transformation forms

1.4 模型构建与验证

1.4.1 偏最小二乘回归模型 Wold等[16]首次提出了偏最小二乘回归(Partial Least Squares Regression，PLSR)算法。此算法涵盖了普通的多元回归分析、主成分分析和相关性分析，同时保留了这三种回归分析的优点，是以往线性回归的优化算法。先以TM影像多光谱波段的原始反射率或经过数学变换后的反射率为基础，通过KL主成分变换得到各个主成分，再以研究区土壤全氮含量为因变量，53个采样点的地表真实反射率值或经过数学变换之后的反射率为自变量，借助SPSS软件建立PLSR模型，偏最小二乘回归方程如式(1)所示：

Y=a1X1+a2X2+…+a2Xm

(1)

式中:Y代表土壤全氮含量，X代表回归方程的自变量，系数ai对应本文所研究的遥感影像多光谱波段的原始反射率和经过数学变换之后的反射率，取自每个主成分对应的系数值。对每个主成分进行回归分析，得到回归决定系数R2和均方根误差RMSE，取R2最大及RMSE最小时对应的主成分建立PLSR回归模型。

1.4.2 神经网络模型 神经网络(Neural Network)是以工程技术手段模拟人脑神经元网络结构与功能的系统,由大量简单的非线性处理单元组成。人工神经网络通过网络中神经元群体的相互作用体现自身的处理功能,可以处理模糊的、非线性的、含有噪声的资料，特别适用于处理非线性问题，因而在模式识别、图像处理和自动控制等方面获得了广泛应用。神经网络模型由输入层、隐藏层和输出层组成，其中输入层和输出层各1个[17]，隐藏层可以为多个。本文采用的BP神经网络模型是比较常用的神经网络模型之一，其结构示意图如图3所示。

图3 BP神经网络模型Fig.3 BP neural network model

本文构建的BP神经网络模型包含3层，其中只包含1个隐藏层。输入层的神经元个数为通过相关性分析得到特征波段数目，以土壤全氮含量分析得到的特征波段确定神经网络的输入层的神经元的节点数，本研究中输出层的神经元数为1个，即土壤全氮含量，输入层节点数为光谱波段个数，即6个神经元数，隐藏层的节点数为10个。设置模型学习率为0.01，最大8 000次的迭代次数，BP神经网络模型的传递函数采用sigmoid函数计算隐藏层的输出和最终输出层的结果。BP神经网络模型的各个参数需要经过反复试验后比较运行结果而确定。

1.4.3 精度评价指标 本文选用标准回归评价中的回归决定系数R2和误差指数评价中的均方根误差RMSE评价模型的稳定性和预测精度。R2取值为[0，1]，模型的拟合程度越高，即实测值与预测值越接近，R2越接近于1；RMSE的值越小，说明模型预测值与实测值偏差较小，证明预测精度越高，土壤全氮含量反演误差越小[18,19]。具体公式如(式2和式3)。

(2)

(3)

2 结果与分析

2.1 偏最小二乘回归

最佳的主成分数、模型建模样本及验证样本的回归决定系数R2和均方根误差RMSE如表3所示。本文构建的两个模型R2都在0.5～0.6，根据反射率倒数构建的模型R2最大为0.604，RMSE为0.285，通过分析得出的全氮反演模型可以对土壤全氮含量进行预测。

表3 全氮PLSR模型评价结果比较Tab.3 Comparison of evaluation results of total nitrogen PLSR model

通过以上分析绘制出土壤全氮含量实测值与根据反射率倒数构建的模型反演预测值的散点图，如图4所示。根据各个样本点的离散程度可知，以1/R为自变量构建的模型预测效果最佳，对于土壤全氮含量的预测有较高的精度。

图4 PLSR 1/R模型土壤全氮含量预测值与实测值比较Fig.4 Comparison between predicted and measured total nitrogen content of PLSR 1/R model

2.2 BP神经网络

全氮BP神经网络模型反演评价结果如表4所示。BP神经网络模型的预测能力较好，其均方根误差RMSE最小值为0.246。根据4种反射率构建的模型回归决定系数R2均大于0.6，可以对土壤全氮含量进行定量反演。验证样本的回归决定系数R2均略低于建模样本的回归决定系数R2，原因可能是采样点的土壤全氮含量分布不均匀。总体而言1/R模型的预测能力高，运用该模型绘制了土壤全氮含量预测值与实测值的散点图(图5)。相较而言，BP神经网络模型具有很大的优势，其预测能力和稳定性都高于PLSR，这与肖云飞[14]的研究结果一致。

图5 BP神经网络1/R模型土壤全氮含量预测值与实测值比较Fig.5 Comparison between predicted and measured total nitrogen content of BP neural network 1/R model

表4 全氮BP神经网络模型反演评价结果Tab.4 Inversion evaluation results of total nitrogen BP neural network model

2.3 研究区土壤全氮含量反演

利用所建立的偏最小二乘回归模型和BP神经网络模型比较预测精度，结合遥感影像并选择最优模型反演研究区土壤全氮含量，以精度最高的BP神经网络模型对土壤全氮含量反演，如图6所示。从图中可以看出研究区西北部的土壤全氮含量较高，南部城区的土壤全氮含量相对较少。

图6 研究区土壤全氮含量分布图Fig.6 Distribution of soil total nitrogen content

3 讨论与结论

土壤全氮含量是衡量土壤肥力的重要因素之一。实时获取土壤全氮含量，利用遥感影像进行大范围土壤全氮反演，对于满足现代精准农业的要求、为土壤质量评价和粮食估产等工作提供基础数据具有重要意义。近年来，国内外许多学者关注土壤全氮研究以及土壤光谱与氮素的关系。高灯州等[20]研究了闽江口湿地的土壤全氮含量，发现在波长0.5 μm附近土壤全氮含量的负相关系数最高；吴明珠等[21]研究亚热带地区土壤全氮含量与光谱反射率之间的相关性，结果表明在波长0.35～2.5 μm光谱反射率与土壤全氮含量均呈现负相关关系。彭杰等[22]分析研究了4种土壤，所构建的多元逐步回归模型预测精度高于一元线性回归模型，其中根据对数倒数的导数所构建的模型土壤全氮含量预测值与实测值的回归决定系数R2达到0.837 2。本文以大通县Landsat 5遥感影像和采样地数据为基础，基于多元线性回归分析方法和主成分分析法，分别构建了偏最小二乘回归模型和BP神经网络模型，反演出整个研究区的土壤全氮含量空间分布状况；利用统计软件对采样点表层土壤全氮含量与对应点的TM影像各波段SR值进行了相关分析，并得出了以下主要结论：

共有7个参与相关性分析的Landsat 5光谱响应波段，其中B2波段最佳，对原始光谱反射率进行数学变换预处理之后，B2波段的反射率1/R的相关系数r达到了0.584，P<0.01；对比偏最小二乘回归模型和BP神经网络模型可知，非线性模型的预测精度优于线性模型，根据反射率倒数构建的BP神经网络模型回归决定系数R2为0.792，均方根误差RMSE为0.246，明显优于偏最小二乘回归模型，预测能力较好。