基于高光谱遥感技术对土壤氮磷钾的估算

文章编号： 1006-8023（2018）06-0025-08

Estimates of Soil N， P， K Concentration by Using Hyperspectral Remote Sensing Technology

QIAO Lu1， CHEN Lixin2*， DONG Chengming1

（1. College of Pharmacy， Henan University of Traditional Chinese Medicine， Zhengzhou 450000; 2. College of Forestry， Northeast Forestry University， Harbin 150040）

Abstract： The purpose of this paper is to study the feasibility of real-time monitoring soil total nitrogen （TN）， total phosphorus （TP） and total potassium （TK） in large area by using Hyperspectral Remote Sensing technology. Based on the combined methods of field investigation and indoor data from hyperspectral measurement（350-2 500 nm）， the soils in four counties（zhaoyuan county， zhaozhou county， dumeng county， and lindian county）， and daqing urban area， Heilongjiang province were studied and analyzed for the characteristics of soil spectrum and soil nutrient index. The variation in soil spectral reflectance was calculated in the form of first-order differential （R）′， reciprocal （1/R）， reciprocal first-order differential （1/R）′， logarithm （log R） and so on. Then the estimating models of TN， TP and TK were established by unary linear regression（ULR）， multivariate linear regression （MLR） and partial least squares （PLSR）， and calibration model is verified by independent samples. The results showed that the optimum absorption bands of TN， TP and TK were TN （802 nm， 1 108 nm）， TP （784 nm， 1 085 nm）， TK （1 079 nm， 1 578 nm）， respectively， the precision of partial least squares （PLSR） was higher than that of Multivariate linear regression（MLR） model and unary linear regression model（ULR）， which were TN （R2 =  0.831， RMSE = 2.506g/kg）， TP （ R2 = 0.687， RMSE = 0.844 g/kg ）， TK （R2 = 0.832， RMSE = 0.097g/kg）. The prediction accuracy of TN and TK was higher， while the prediction result of TP was relatively lower， but it can also be used for rough estimation. Meanwhile， MODIS images were used to map the concentrations of TN， TP and TK in soil. This study confirmed that hyperspectral technology combined with specific algorithms can better predict the concentrations of TN， TP and TK in soils with larger differences， and the visualization of predictive information can be realized. It is of great practical significance for real-time and rapidly monitoring large-scale changes in soil environment， predicting the trend of soil information changes， monitoring ecological environment， establishing soil nutrient database in China and reducing the monitoring cost of soil composition.

Keywords： Hyperspectra; TN; TP; TK; multivariate linear regression （MLR）; partial least squares regression （PLSR）; MODIS image

0 引言

氮磷钾作为土壤中的大量元素，是植物生长发育的必要营养成分。及时有效地了解土壤中氮磷钾的潜在储备，对植物的生长发育、土壤生态环境监测和精准农业发展等方面至关重要[1]。近年来高光谱遥感技术的发展，对及时有效的获取土壤信息起到极大的帮助作用。

国内外学者先后利用高光谱遥感技术实现并发展了对土壤氮磷钾的预测，Rossel等建立土壤光谱PLSR [2]，发现近红外波段NIR对磷（P）的估算较好，而中红外波段MIR能更精确地估计钾（K）含量。A. M. Mouazen等用PLSR联合完整的交叉验证技术对法国北部和比利时土壤中的有效磷建模效果显著[3]。M. R. Maleki等研究显示[4]，室外土壤光谱测定所提取磷的吸收带高于室内风干土测定的结果。A. M. Mouazen等比较了3种模型PLSR、主成分分析（PCA）和反向传输神经网络（BPNN-LVs）[5]，指出BPNN-LVs能更好的反演土壤氮磷钾含量，并证实BPNN-LVs模型反演磷精度要高于钾模型精度。Shao Yong-ni研究发现最小二乘法-支持向量机（LS-SVM）结合近红外光谱能更有效的获取土壤氮磷钾信息[6]。徐永明利用土壤光谱吸收带特征与总氮含量进行逐步回归运算，确定了与氮元素关系比较密切的几个吸收带，包括一阶导数、倒数 （1/R）、倒数之对数（log（1/R））、波段深度（Depth）[7]。袁石林利用均值法、卷积滤波确定了光谱数据[8]，采用PLSR和（LS-SVM）分别建立了浙江省潮化盐土土壤总氮和总磷的近红外光谱模型。陈鹏飞利用傅里叶变换光谱技术和PLSR建立了北京大兴地区土壤总氮和总磷含量的近红外光谱校正模型 [9]。宋海燕采用PLSR建立了一阶微分光谱的光谱吸光度与速效磷、速效钾含量和pH值之间的定量分析模型，结果表明pH值的预测效果比磷、钾好[10]。高洪智利用连续投影算法能够提取土壤总氮的近红外特征波长 [11]。本文拟通过利用土壤反射光谱样本建立土壤全氮、全磷和全钾含量的一元线性回归、多元线性回归和偏最小二乘法模型，估算大庆土壤全氮、全磷和全钾养分含量，对比3种模型的精度差异。为高光谱在土壤养分的监测和估算提供理论依据。

1 实验部分

1.1 样品采集及分析

大庆市位于黑龙江省松嫩平原中部。地理位置在北緯45°46′ ～ 46°55′，东经124°19′ ～ 125°12′之间。土壤类型包括黑钙土、草甸土、盐碱土和沼泽土。在大庆五区和四县（图1）（肇源县、大庆市区、肇州县、杜蒙和林甸县）按照土地利用方式划分为耕地、草地、荒地和林地，分别采集土壤耕作层（0 ～ 20 cm）90个样本。采集后样本剔除树叶，根枝，自然风干后研磨、过2 mm筛，一部分用于土壤全氮、全磷和全钾的测定，一部分用于室内土壤光谱的采集。土壤全氮用凯式定氮法测定;土壤全磷用钼锑抗比色法测定;土壤全钾用氢氟酸-高氯酸消煮法[12-13]。

1.2 光谱测定

采用美国ASD FieldSpec-FRTM背挂式野外光谱辐射仪。测量时，分别设置350 ～ 1 000 nm波段，光谱分辨率为3 nm;1 000 ～ 2 500 nm波段，光谱分辨率为10 nm。光谱测试在暗室进行，每次测量前后都进行标准版的测量以校准光谱仪。为减少误差，每个样本重复测量10次，最后取平均值。

1.3光谱数据预处理

利用SAVITZKY-GOLAY（SG）滤波器进行对原始的土壤样本光谱进行平滑和消除包络线等预处理（图2），去除包含在光谱数据中的噪声。土壤全氮（TN）、全磷（TP）、全钾（TK）在土壤反射光谱中没有明显的特征，难以提取特征波段。对光谱反射率做微分、倒数、对数等数学形式变换可降低部分噪声、背景、地形和光照等因素对光谱数据的影响。提高光谱数据与土壤氮磷钾含量的相关性。本论文的光谱反射率的处理方式主要包括：光谱反射率对数logR，反射率倒数1/R，光谱反射率对数的倒数（1/lgR）、光谱反射率倒数之一阶微分（1/R）′、光谱反射率对数的倒数之一阶微分（1/lgR）′5种形式。所有这些预处理都是在Matlab2007， Origin7.5软件中进行。

1.4 土壤养分的估算模型

将土壤样品按照3∶1的比例分成建模样本和检验样本，在分析光谱数据与TN、TP、TK含量相关性的基础上，应用一元线性回归、多元线性回归、偏最小二乘法模型进行拟合。决定系数R2 用于判定模型的稳定性，均方根误差RMSE用于评估模型的反演能力。较高的判定系数和较小的均方根误差表明该模型更加稳定和精确[14]。

2 结果与分析

2.1 土壤TN、TP、TK含量分析

表1是大庆市90个土壤样品中全氮、全磷和全钾的常规分析结果及标准差，从表1中看出各营养成分含量的变化幅度较大，土壤类型较多，这对于借助高光谱遥感技术大面积土壤养分模型的建立起到一定的作用。

2.2 特征波段选择

对经过预处理后得到的光谱反射率R、 logR、 1/R、（1/lgR）、（1/R）′和 （1/lgR）′分别与土壤的TN、 TP、TK含量进行相关分析，得到结果如图3所示。从图3（a）中可以看出，土壤TN与光谱反射率相关系数变化范围分别-0.45 ～ 0.63，在波段在600 ～ 700 nm，1 100 ～ 1 600 nm具有显著的相关性。图3（b）中，TP 在整个波段反射光谱均表现出了较低相关性。只在波段600 ～ 800 nm内相关系数达到0.5。经过一阶微分后的反射率对数在波段2 100 ～ 2 300 nm的正相关显著的提升。图3（c）中，TK与光谱反射率对数在波段1 300 ～ 1 500 nm，1 900 ～ 2 100 nm呈显著正负相关。据分析结果，TN入选波段为802、1 108、1 558 nm，P为784、1 046、1 085 nm，TK为1 070、1 578 nm。

2.3模型建立与预测

2.3.1 一元线性的逐步回归模型

将采集90个样本随机提取22个样本作为检验模型精度，68个样本建立氮磷钾的一元线性模型。从表2中可以看出，土壤TN的预测模型中，采用R1108nm 变量建立的模型精度最高 R2 = 0.821、RMSE = 1.821 g/kg，采用logR802nm 变量的模型精度达到R2 = 0.771、 RMSE = 1.771 g/kg; TP的预测模型中，以R784nm 、1/logR1085变量建立的模型，精度较高分别是 R2 = 0.622、RMSE = 1.222 g/kg和R2 =  0.595、RMSE = 1.295 g/kg;TK的预测模型中，以R1079nm 、1/logR1578nm 为变量的模型的精度分别为 R2 = 0.642、 RMSE = 0.156 g/kg和 R2 =  0.668、 RMSE = 0.147 g/kg。这与Rossel[2]发现近红外波段NIR对磷（P）的估算较好、而中红外波段MIR能更精确地估计钾（K）结论相一致。由于选取变量较多，计算方差膨胀因子（Variance Inflation Factor，VIF）均小于10，说明变量间不存在多重共线性。

2.3.2多元线性的逐步回归模型

利用相关分析的结果，提取较高变量，建立多元线性回归模型（表3），预测模型的结果为TN（R2 = 0.824、 RMSE = 2.01 g/kg），入选波段802、1 558 nm;TP （R2 = 0.655、RMSE = 0.985 g/kg），在784 ～ 1 050 nm 内的波段入选回归方程数量最多，说明该吸收带范围的波段与P元素的关系最为密切。TK （R2 =  0.722、RMSE = 0.134 g/kg），TK的入选波段在1 070、1 578 nm波段。研究中可以看出，提取多波段，多种光谱反射率变化形式共同作为自变量，所建立的MLR精度高于ULR，且模型更加稳定。

2.3.3偏最小二乘法模型与验证

偏最小二乘法（PLSR）是集多元线性回归法、主成分分析法、相关性分析法等功能为一体的综合分析法，它可以很好地解决自变量之间多重共线性、样本容量小的问题。利用PLSR建立的预测模型结果显示：TN（R2 = 0.831，RMSE = 2.509 g/kg），TP（R2 = 0.687，RMSE = 0.844 g/kg）TK（R2 = 0.832，RMSE = 0.097 g/kg）（表4）。预测精度R2和RMSE均高于MLR和ULR。将22个未参与建模样本代入模型。绘制预测值与实测值的散点，如图3所示。 检验样本集中在趋势线y = x附近，表明模型精度高，稳定性强。图3中预测值和实测值的相关系数为：TN（R2 = 0.893，RMSE = 2.556 g/kg），TP（R2 = 0.858 8，RMSE = 0.951 g/kg），TK（R2 = 0.907 1，RMSE = 0.083 g/kg），这也说明验证样本的预测的效果较好。

3 偏最小二乘法模型在MODIS影像上的推广

土壤养分含量高光谱反演地面模型无法在大面积上进行估算，需要结合遥感影像，然而直接将地面模型推广到遥感影像上，往往又十分困难，因为光谱分辨率存在不一致的问题。研究采用多光谱MODIS影像数据进行土壤TN、TP、TK含量的專项制图。中分辨率成像光谱仪（Moderate-resolution Imaging Spectroradiometer，MODIS）影像的传感器主要是搭载Terra和Aqua两颗卫星上，通过两颗卫星相互配合平均每24 ～ 48 h对整个地球表面重复观测，MODIS影像一共有36个波段的多光谱数据，可覆盖研究区。采用MODIS影像的响应函数对地面高光谱数据进行重采样，将土壤养分含量的地面模型推广到MODIS影像上。

3.1数据获取及处理

根据大庆地区的地理坐标，从NASA网站免费下载MODIS影像（表5）：标准陆地2级数据产品，其内容为表面反射，空间分辨率为250 m。

下载时间为2008年第281 d，研究中使用的MODIS数据已经过几何校正，辐射校正。用ENVI软件的FLAASH模块来完成MODIS影像的大气校正（图5）。目的是消除大气和光照等因素对地物反射的影响，获得地物反射率和辐射率、地表温度等真实物理模型参数。利用ENVI提取MODIS影像（01-07）的各波段的响应函数（图6），将其与地面测定的高光谱反射率数据相乘，对地面光谱数据进行重采样，获得与MODIS影像（01-07）波段相对应的地面光谱数据。将重采样后的高光谱数据的7个波段反射率进行数学运算，包括：光谱倒数1/R，对数logR，对数倒数1/logR，NDVI运算、平方根等土壤光谱指数。

3.2 PLSR模型与检验

选取地面模型中精度最高的PLSR推演到MODIS影像中。随机以50个样本作为建模样本，40个检验样本，建立PLSR模型（表6）。由图7可以看出，大多土壤指标的实测值与预测值在趋势线y = x处聚集，个别点较分散，精度较高，RMSE较低，模型稳定。TP的建模和检验精度相对较低，精度只有0.566，只能粗略估算。

3.3 TN、 TP、TK的专项制图

根据PLSR模型的结果：

YTN = 1.702-10.145×RNDVI3+2.29×（1/logb7）+16.649×R6+0.11×（1/R1）。

YTP = 0.899+0.728×（1/logb7）+3.798×R6。

YTK = 44.789+85.689×（1/logb6）+23.631×5- 136.106RNDVI6-62.190×（1/logb7）+8.517RNDVI4。

用以上这些模型对土壤TN、TP、TK专项制图，如图8所示。从图8中可以看出，土壤TN含量较分散，但具有一定的规律性。为达到细分专项制图的效果，根据所测土壤TN含量最大值、最小值及平均值分段统计，共分为8个区间，对应不同颜色使数据可视化。大庆地区的土壤TN主要含量在1.62 ～ 1.9、1.9 ～ 2.28、2.28 ～ 2.66、02.66 ～ 3.049 g/kg含量。TN含量的地区差异：肇源县>肇州县>大庆地区>林甸县>杜蒙。根据TK所测含量分段统计，共分为6各区间，TK主要含量在18 ～ 20、20 ～ 22、22 ～ 24 g/kg 3个范围。TK含量的地区差异：杜蒙>肇源县>肇州县>大庆市区>林甸县。根据TP所测含量分段统计，共分为6个区间，土壤TP主要含量在0.65 ～ 0.78、0.78 ～ 0.91 g/kg 之间， TP含量的地区差异：肇源县>大庆市区>肇州县>杜蒙>林甸县。

4 讨论

在3种模型中，偏最小二乘法得到TN、TP、TK的反演模型精度，均高于一元线性模型和多元线性回归模型，其决定系数R2 最高，均方根误差RMSE最小，TN（R2 = 0.831、RMSE = 2.509 g/kg），TP（R2 = 0.687、RMSE = 0.844 g/kg）TK（R2 = 0.832， RMSE = 0.097 g/kg）。TN、TK估算效果较好，是由于土壤中全氮绝大多数贮藏在土壤有机态含氮化合物中，而土壤中钾多来自于正长石、微斜长石和白云母等原生矿物及伊利石等次生矿物，占全钾的90%以上，利用光谱较容易定量估算。磷含量模型精度相对较低，可能是因为土壤中全磷含量较低，且难溶性磷化合物占土壤无机磷的绝大部分，属植物难以利用的迟效磷，不易吸收。但与已有研究相比，本研究得到的全磷光谱模型精度较高，是因为已有研究大多以单波段反射率或多波段反射率组合的方式作为模型自变量[15]，而本文采用的是多波段、多种反射率变换形式相组合作为模型自变量，所得参数是基于一元线性回归模型获得，所建立的模型稳定性更强。

已有的研究主要是针对特定的土壤利用方式，或是单一的土壤类型所建立的预测模型，模型精度较高，但只适用于局部地区土壤的研究，不适合大面积土壤养分的估算。本论文主要研究区土壤检测面积广，土壤类型复杂，利用方式多样化，养分含量差异较大，与已有的研究相比较，模型精度相对较低，但对于模型的通用性的建立提供了理论依据。

5 结论

本文通过对土壤光谱和TN、 TP、TK含量的相关性分析，确定了与土壤氮、磷、钾含量相关性较高的吸收带，利用一元线性回归、多元线性回归、偏最小二乘法模型对土壤全氮、全磷、全鉀含量建模，得到结果如下：

（1）土壤氮磷钾含量与土壤光谱中700 ～  1 200 nm，1 500 ～ 1 600 nm存在显著的相关性，经对数、一阶微分转换后，光谱反射率与氮磷钾含量相关性得到显著增强。表明对数、微分光谱有助于提取土壤元素信息。全氮的主要吸收带在波段（802、1 108、1 558 nm），全磷（784、1 046 ～ 1 100 nm），全钾（1 070 ～ 1 100、1 550 ～  1 600 nm）。

（2）以偏最小二乘法建立的TN、 TP、TK模型精度优于一元线性模型和多元线性回归模型，其结果显示：全氮N（R2 = 0.831， RMSE = 2.509 g/kg），全磷P（R2 = 0.687， RMSE = 0.844 g/kg），全钾K（R2 = 0.832，RMSE = 0.097 g/kg）。

（3）土壤TN、TK的估算效果较好，TP的模型精度相对较低，高光谱遥感可以有效的估算土壤TN、 TP、TK含量。

（4）利用高光谱波段模拟MODIS多波段，实现“窄波段”变“宽波段”。建立PLSR土壤估算模型，大多数土壤指标模型估算精度能够达到0.7以上，可实现精确的估算;TP模拟R2相对较低，可用于粗略估算。同时提供了大庆地区土壤TN、 TP、TK养分的MODIS专项制图，实现了预测信息可视化。其对于实时快速监测大面积土壤环境变化、预测土壤信息变化趋势、预警生态灾难、建立我国土壤养分数据库和降低土壤成分监测成本具有重要的现实意义。

【参 考 文 献】

[1]童庆禧， 张宾，郑兰芬.高光谱遥感的多学科运用[M]. 北京：电子工业出版社， 2006.

TONG Q X， ZHANG B， ZHENG L F. Hyperspectral remote sensing