HJ-1A HSI与Sentinel-2A遥感数据土壤全氮含量反演精度的对比研究

马 驰

(辽宁省交通高等专科学校 测绘工程系，辽宁 沈阳 110122)

0 引言

氮元素是土壤养分的重要组成部分，是植物生长的必要元素，土壤中氮元素的含量直接影响植物的长势及农作物产量。因此，精确、快速量测土壤氮元素含量对于改善土壤质量、提高农作物产量乃至精细农业的实施均具有重要意义[1-2]。传统的利用化学试剂测试土壤氮元素含量的方法耗费大量的人力与物力，难以实现土壤中氮元素含量的实时监测。遥感技术具有数据获取时间短、遥感影像覆盖范围广、所含信息量丰富等优势。近年来逐渐成熟的高光谱技术由于具有更多的光谱通道和更高的光谱分辨率，能更加精细地反映地物的光谱特征和光谱差异，为土壤成分快速监测提供了更有效的数据源[3-4]。

近年来，国内外诸多学者相继展开了利用遥感技术探测土壤全氮含量等方面的研究。Brunet D等[5]研究表明，将土壤光谱数据进行一阶微分处理后可以有效地提升土壤氮元素含量的反演精度；卢艳丽等[6]利用ASD2500 高光谱仪在实验室测绘土壤样品的光谱曲线，并将光谱反射率进行对数的一阶微分变换后与土壤全氮相关系数在556，1642， 2 491 nm处出现峰值；王世东等[7]利用FieldSpec 3型光谱仪测定永城市矿区土壤光谱曲线，对土壤光谱进行一阶微分、二阶微分、连续统去除数学变换，并利用偏最小二乘分析的方法建立研究区土壤全氮含量的反演模型，模型的决定系数R2达到了0.92；方向等[8]使用OFS-1700型便捷式地物非成像光谱仪对黄山市黄山区和池州市石台县土壤样品进行光谱分析，对土壤光谱进行一阶导数、二阶导数、对数、去趋势校正等数学变换，并利用偏最小二乘分析方法建立研究区土壤速效氮含量的估测模型，模型的决定系数R2达到了0.94。王文俊等[9]以山西典型褐土为研究对象，利用Starter Kit光谱仪在实验室绘制土壤样本的近红外高光谱曲线，采用平均光谱曲线和平均光谱曲线的一阶微分等方法对光谱数据进行预处理，并利用偏最小二乘的方法建立研究区建立土壤全氮含量预测模型，研究结果表明，使用一阶微分进行建模能够获得更好的预测效果。

综上，当前利用遥感技术研究土壤中全氮含量主要集中于氮元素的光谱特征方面，所使用的光谱数据多为利用便携式光谱仪在实验室采集的土壤光谱数据，将其应用于土壤元素的遥感制图仍需辅以其他遥感影像数据，限制了其实用性。2008年9月发射升空的我国自行研制的环境灾害小卫星(HJ-1A)所搭载的超光谱成像仪(Hyperspectral Imaging Radiometer,HSI)，与多光谱遥感影像相比，拥有多达115个探测波段和精细的光谱分辨率(仅为4.3 nm)，在土壤成分探测等方面已获得了广泛应用[10-13]。然而，由于HJ-1A高光谱影像波段过多而造成各波段间的数据冗余度较高、对其进行数据处理的过程也较复杂；影像的空间分辨率较低(100 m)，难以反映地表细部信息。2016年6月升空的Sentinel-2A遥感卫星搭载的多光谱成像仪，与Landsat TM、ETM等多光谱成像仪相比，具有更高的空间分辨率、光谱分辨率和更短的重访周期，必将在土壤成分探测等方面呈现出巨大优势。

为此，本文试验以中国的环境灾害小卫星高光谱影像(HJ-1A HSI)和Sentinel-2A遥感影像为数据源，辅以研究区土壤采样的全氮含量实验室化验数据，利用SPSS软件统计光谱反射率与全氮含量的相关系数，筛选HJ-1A HSI影像全氮含量的敏感波段，建立两种遥感影像土壤全氮含量的反演模型；通过比较两种遥感影像反演研究区土壤全氮含量的模型精度，分析两种遥感影像在土壤全氮含量方面的探测能力，为Sentinel-2A遥感影像在土壤元素探测方面的应用研究提供参考。

1 材料与方法

1.1 研究区概况

农安县隶属吉林省长春市，东经124°31′～125°45′，北纬43°55′～44°55′，属典型大陆性季风气候，年平均温度为4.9 ℃、年平均降雨量为500 mm。农安县地势平坦，平均海拔约250 m，地貌为冲积、湖积平原区，土壤类型主要包括黑土、黑钙土、草甸土、风沙土、盐碱土等。

1.2 土壤采样及化验

2017年4月29日—4月30日在农安县进行土壤采样，共采集土壤样品82个，其采样点分布如图1所示。土壤采样过程中选择裸土区域，利用GNSS接收机测量采样位置的经纬度坐标，在采样位置20 m×20 m范围内采集表层土壤的多个土样并混合，收集约0.5 kg土壤样品装入土壤采集袋；将土壤样品置于干燥阴凉处风干，剔除土壤样品中的植物茎叶与根须、昆虫躯壳、小石子等杂物并过筛。土壤全氮含量测定采用凯氏蒸馏法。

图1 土壤采样点分布Fig.1 Distribution of soil sampling points

1.3 遥感数据的选取与预处理

选取HJ-1A HSI遥感影像6景(成像时间为2017年4月30日)、Sentinel-2A遥感影像2景(成像时间为2017年4月29日)，遥感影像的云覆盖率均小于2%。将HJ-1A HSI和Sentinel-2A两种遥感影像进行FLAASH辐射校正和大气校正，用以消除两种遥感卫星的传感器在成像过程中大气的辐射误差；分别对HJ-1A HSI和Sentinel-2A遥感影像进行几何精校正、遥感影像的拼接和裁剪等工作，获得覆盖农安县的遥感影像；根据土壤采样时利用GNSS接收机测量获得的采样点坐标，在HJ-1A HSI和Sentinel-2A遥感影像中提取采样点反射率数据。

1.4 相关性分析

将HJ-1A HSI和Sentinel-2A遥感影像的反射率与研究区土壤全氮含量逐波段进行相关性分析如式(1)所示，用以提取研究区土壤全氮含量敏感波段，并获得建立全氮含量反演模型的最佳组合波段。

(1)

为了降低遥感成像中的噪声对土壤光谱的影响，提高HJ-1A HSI和Sentinel-2A遥感影像反射率与土壤元素含量的相关性[5-9，12-14]，本文将两种遥感影像的反射率进行适当的数学变换，包括对数lnR、倒数1/R、一阶微分R′等，并与研究区土壤全氮含量在SPSS软件中进行相关性分析，再筛选出研究区土壤全氮含量的敏感波段。本试验中遥感影像反射率的一阶微分采用光谱差分的方法估算如下：

(2)

式中，R′i为第i波段反射率的一阶微分；Ri为第i波段的反射率；Δλ为第i波段至第i+1波段的光谱间隔。

1.5 土壤全氮含量反演模型的建立与模型检验

将82个土壤样本中的72个作为建模样本、其余的10个作为模型精度的检验样本。将相关性分析所提取的敏感波段反射率与农安县土壤全氮含量进行逐步回归分析，建立农安县裸土全氮含量的反演模型。为了对反演模型进行检验，试验引入模型的决定系数R2以衡量模型的稳定性，引入模型的均方根误差(Root Mean Square Erro，RMSE)以衡量模型的精度。选择具有较大模型决定系数和较小的模型均方根误差的土壤全氮含量反演模型作为研究区最优反演模型，用以反演研究区裸土的全氮含量。

2 结果与分析

2.1 相关性分析与敏感波段的选取

将农安县土壤全氮含量与HJ-1A HSI影像反射率(反射率变换)逐波段进行相关性分析如图2和图3所示，获得农安县裸土全氮含量的敏感波段如表1所示。相关性分析结果显示,研究区土壤的全氮含量与HJ-1A HSI影像反射率呈负相关，并在第28波段(中心波长为524 nm)达到峰值，为r=-0.499。将反射率进行适当的数学变换后，与研究区裸土全氮含量的相关性有所提升。其中，反射率的对数变换与研究区土壤全氮含量的相关系数在第92波段(中心波长为782 nm)达到峰值，为r=0.605。反射率(反射率变换)经过一阶微分后与农安县裸土全氮含量的相关性进一步得到改善，HJ-1A HSI影像反射率的一阶微分变换与土壤全氮含量的相关系数在第92波段(中心波长为782 nm)达到峰值，相关系数为r=0.628，遥感影像的反射率指数一阶微分变换与研究区土壤全氮含量的相关系数在第92波段达到峰值，相关系数达到0.695。

(a) 反射率

(a) 反射率的一阶微分

表1 建模波段反射率(反射率变换)与裸土全氮含量的相关系数

2.2 反演模型的建立

根据相关性分析结果，选择建模波段建立研究区土壤全氮含量的反演模型。建模波段的选取原则为：波段反射率(反射率变换)与研究区土壤全氮含量相关性较好、波段所含信息量丰富；波段间相距较远、波段间的相关性较小、数据的冗余度较低。本试验依据建模波段的选取原则，参照HJ-1A HSI影像反射率(反射率变换)与农安县裸土全氮含量相关性分析结果，选取全氮含量反演模型的建模波段(见表1)。将随机选取的72个建模样本所对应遥感影像的建模波段反射率(反射率变换)作为自变量，以建模样本全氮含量的化验值为因变量，利用SPSS软件进行逐步回归分析，建立农安县裸土区全氮含量的反演模型如表2所示。

表2 HJ-1A HSI土壤全氮含量的反演模型

对比研究区土壤全氮含量反演模型发现，利用HJ-1A HSI影像反射率建立的反演模型的决定系数R2为0.568，RMSE为0.296 g/kg ；利用反射率的数学变换建立的全氮含量反演模型较反射率模型，模型的决定系数普遍有所提升。其中，利用反射率的幂变换建立的土壤全氮含量反演模型，模型决定系数R2达到0.754，模型的RMSE为0.218 g/kg。反射率(反射率变换)经过一阶微分后建立的反演模型的决定系数普遍高于非微分形式反演模型的决定系数。其中，反射率幂变换的一阶微分反演模型的决定系数R2达到了0.806，模型的RMSE减小为0.185 g/kg。

2.3 模型精度检验

利用反射率幂的一阶微分反演模型Y=2.45+0.517X15+4.062X85+7.281X92+2.575X105计算检验样本全氮含量的反演值，与检验样本全氮含量的化验值建立散点图。散点图4显示，研究区土壤检验样本全氮含量的反演值和实测值较均匀的分布于1∶1直线两侧，其拟合模型为Y=1.062X-0.097，拟合模型的决定系数R2=0.912，模型的RMSE=0.162 g/kg。

图4 反演值与实测值散点图Fig.4 Scatter diagram of inversion value and measured value

相对误差散点图5显示，10个检验样本中有6个检验样本的相对误差δ在(-0.1,0.1)之间，占检验样本总数的60%，有2个检验样本的相对误差δ在(-0.2,-0.1)或(0.1,0.2)之间，表明HJ-1A HSI反射率幂的一阶微分反演模型Y=2.45+0.517X15+4.062X85+7.281X92+2.575X105具有较好的预测效果，可以应用于研究区土壤全氮含量的反演。

图5 相对误差散点图Fig.5 Scatter diagram of relative error

2.4 多光谱建模

根据表1中HJ-1A HSI影像的土壤全氮含量敏感波段的分析结果，参考Sentinel-2A影像各波段的中心波长及光谱范围，选择建模波段并采用回归分析的方法，建立农安县裸土全氮含量的Sentinel-2A多波段反演模型，如表3所示。建模结果显示，非微分模型中，Sentinel-2A遥感影像的对数模型精度较高，反演模型的决定系数R2=0.634，RMSE=0.235 g/kg；而一阶微分反演模型按反射率倒数模型、反射率模型、反射率指数模型、反射率对数模型、反射率幂模型的精度顺序递增，且模型间的反演精度差异较大，即：反射率幂的一阶微分反演模型的预测能力最强，模型决定系数R2最大，达到0.768，模型的RMSE最小，为0.195 g/kg；反射率倒数的一阶微分反演模型的预测能力最差，模型的决定系数R2仅为0.541。

表3 Sentinel-2A土壤全氮含量的反演模型

利用Sentinel-2A遥感影像反射率幂的一阶微分反演模型Y=-0.254+2.6687X2+0.508X6-5.164X7-0.150X8A反演农安县土壤检验样本的全氮含量，与实验室土壤全氮含量的实测值建立散点图,如图6所示。由图6可以看出，检验样本全氮含量的反演值与实测值较均匀的分布于1∶1直线两侧，拟合模型为Y=0.918X+0.264，拟合模型的决定系数R2=0.892，RMSE=0.196 g/kg。相对误差散点图(图6)显示，10个检验样本中有5个样本的相对误差δ在(-0.1,0.1)之间，占检验样本总数的50%，有4个检验样本的相对误差δ在(-0.2,-0.1)或(0.1,0.2)之间，表明利用Sentinel-2A多光谱遥感影像反演研究区土壤全氮含量是可行的。

图6 反演值与实测值散点图Fig.6 Scatter diagram of inversion value and measured value

图7 相对误差散点图Fig.7 Scatter diagram of relative error

2.5 土壤全氮含量制图

根据表2和表3建模结果，选择HJ-1A HSI及Sentinel-2A土壤全氮含量的最优反演模型，反演研究区土壤全氮含量，如图7和图8所示。研究区全氮含量的反演结果图显示，利用HJ-1A HSI和Sentinel-2A遥感影像反演的研究区土壤全氮含量具有相似的空间分布，即：土壤全氮含量大于3 g/kg的区域主要集中于农安县东部和南部，土壤全氮含量小于2 g/kg的区域主要集中于农安县西部，农安县的中部和北部地区土壤全氮含量多介于2～3 g/kg。对研究区进行调研，土壤采样时发现，农安县东部、南部与我国东北黑土带毗邻，土质肥沃、土壤全氮含量较高，农安县西临松原市，为黑土、黑钙土与盐碱土过渡地带，土壤含盐量升高而全氮含量降低，在遥感影像中表现出较高的反射率。

图8 HJ-1A HSI 全氮含量反演结果Fig.8 Total nitrogen inversion results chart of HJ-1A HSI

图9 Sentinel-2A 全氮含量反演结果Fig.9 Total nitrogen inversion results chart of Sentinel-2A

3 讨论

研究结果显示，HJ-1A HSI遥感影像在可见光波谱范围与研究区土壤全氮含量具有良好的相关性，且在红光、绿光、蓝光区域均存在峰值，与文献[4，6，15]的研究结论相同或相近。将HJ-1A HSI遥感影像反射率进行倒数、对数、指数、幂、一阶微分等数学变换后可以有效改善与土壤全氮含量的相关性，与文献[5-9，16-18]的研究结论相同。其中，HJ-1A HSI遥感影像反射率的对数变换与研究区土壤全氮含量的相关性在第92波段(中心波长为782 nm)达到最大值，为r=0.605；经过一阶微分变换后的遥感影像反射率后可以有效提升与农安县裸土全氮含量的相关性。其中，反射率指数一阶微分与土壤全氮含量在第92波段达相关性最好，为r=0.695。利用回归分析方法建立的HJ-1A HSI影像反射率幂的一阶微分反演模型(模型的判定系数为R2=0.806，RMSE=0.185 g/kg)，其模型精度略优于Sentinel-2A影像反射率幂的一阶微分反演模型(模型判定系数R2=0.768，RMSE=0.195 g/kg)，究其原因，可能与HJ-1A HSI遥感影像的光谱分辨率远高于Sentinel-2A遥感影像有关，例如，HJ-1A HSI影像的全氮含量敏感波段第28波段，波谱范围为5 nm，而与之相对应的Sentinel-2A遥感影像绿光波段的波谱范围为35 nm，HJ-1A HSI影像的精细光谱大大提高了其对土壤全氮的识别能力。

利用HJ-1A HSI遥感影像反射率幂的一阶微分建立的研究区土壤全氮含量的最优反演模型Y=2.45+0.517X15+4.062X85+7.281X92+2.575X105与参考HJ-1A HSI全氮含量敏感波段而建立的Sentinel-2A全氮含量的最优反演模型Y=-0.254+2.6687X2+0.508X6-5.164X7-0.150X8A均有较高的模型精度，绘制的全氮含量反演结果图具有相似的全氮含量空间分布。究其原因，可能与以下因素有关：① HJ-1A HSI遥感影像具有较高的光谱分辨率，大大提高了对土壤成分的识别能力；② 土壤采样时间为4月29日—30日，与遥感影像的获取时间几乎同步，此时刻研究区地表已无冰雪及绿色植被，遥感影像能够真实的反映土壤采样时刻的裸土信息；③ 通过对HJ-1A HSI、Sentinel-2A遥感影像进行了大气校正，消除了传感器成像时噪声对反射率的影响；④ Sentinel-2A遥感影像光谱分辨率虽然低于HJ-1A HSI遥感影像，但Sentinel-2A遥感影像具有更高的空间分辨率，更容易显示地表的细部特征。

然而，本试验中的检验样本全氮含量反演值与实测值之间仍然存在一定误差，可能与以下因素有关：① 在利用FLAASH模型对两种遥感影像进行大气校正时，由于所输入的大气能见度、大气气溶胶浓度、研究区地表高程等参数与实际存在一定误差，使校正后影像的反射率与地表真实反射率存在误差；② 本试验未考虑土壤中有机质、盐分等对土壤光谱的影响，使土壤全氮含量反演值存在一定误差，将在以后研究中进一步探究。

4 结束语

本文试验利用Sentinel-2A遥感影像代替HJ-1A HSI遥感影像反演农安县土壤全氮含量，获得了以下结论：

(1) HJ-1A HSI遥感影像在可见光波段的反射率与农安县土壤全氮含量具有较强的负相关性，并在第28波段(中心波长为524 nm)相关性最好。

(2) 经过适当的数学变换后，HJ-1A HSI遥感影像反射率与农安县裸土全氮含量的相关性，较变换前有所改善；反射率(反射率变换)经过一阶微分变换后可以进一步提升与裸土全氮含量的相关性。

(3) HJ-1A HSI遥感影像反射率经过幂的一阶微分变换后所建立的农安县土壤全氮含量反演模型，其模型精度最高，模型判定系数R2达到0.806；参考HJ-1A HSI遥感影像的土壤全氮含量敏感波段，建立Sentinel-2A遥感影像的土壤全氮含量反演模型，最优模型的判定系数R2达到0.768，2个模型具有相近的模型精度，表明用于反演研究区土壤全氮含量的Sentinel-2A遥感影像可以代替HJ-1A HSI遥感影像。

(4) 利用HJ-1A HSI与Sentinel-2A遥感影像的土壤全氮含量最优反演模型，反演研究区土壤全氮含量并制图，均获得了较好的反演效果。两幅反演结果图显示，研究区土壤全氮含量具有相同的空间分布，即：农安县东部、南部土壤全氮含量较高，中部、西部土壤全氮含量较低。