基于地面光谱联合SAR多源数据的农田表土氮磷监测

孙宇乐，屈忠义，刘全明

（内蒙古农业大学水利与土木建筑工程学院，呼和浩特 010018）

0 引言

【研究意义】土壤养分是植物生存生长所必需的物质，尤其是在耕种频繁的河套灌区，监测土壤养分状况对农业生产管理具有重要意义。【研究进展】通常国内外土壤氮磷量测定采用的多为有损的化学方法，采样及预处理的复杂烦琐既增加了成本，又延误了时间[1]。高光谱数据具有准确提取和监测土壤信息的技术优势，能够较为完整连续地提供土壤光谱信息，从而可以清晰反映地物光谱的细微特征。大量研究发现，350～2500nm 波段的高光谱反射率及其变换形式可以刻画土壤各参数的微小差别[2-4]。He 等[5]

研究发现种植区土壤中的氮素等指标可以通过红外光谱技术进行有效预测。Dalal 等[6]利用近红外光谱估测了种植区的多种土地利用类型土壤中全氮，发现1 744、1 870 和2 052 nm 波段的估测效果最佳。徐丽华等[7]通过对土壤风干样本光谱曲线进行去包络线处理，构建了土壤中总氮、总磷的偏最小二乘法预测模型。张娟娟等[8]构建了基于特征光谱指数、偏最小二乘法和BP 神经网络的土壤全氮量估算模型，对我国中、东部地区主要土壤类型的土壤全氮量进行了估测。韩磊等[9]基于BP 神经网络建立了土壤养分综合评价模型，对安塞区土壤养分进行了综合评价。【切入点】通过国内外众多专家学者的研究，发现利用光谱分析技术预测土壤氮磷量是可行的，但众多学者对于土壤氮磷量的遥感估测研究多集中在近红外光谱方面，而随着遥感技术的迅猛发展以及地物光谱仪分辨率的提高，现在可以更加精准地对土壤氮磷量进行测定；同时，河套灌区存在着大量的盐碱化土壤，其呈斑状分布，使得表层土壤变化状况较为复杂且土壤结构差异性较大，现有的遥感预测模型不能够准确、便捷地反映不同区域的实际土壤氮磷状况及其变化特征。RADARSAT-2 是一款具有全天候、全天时特点的微波雷达，由于雷达不受雾、云和雨的阻挡，大量研究人员利用其对表层土壤水分进行了研究[10-12]。【拟解决的关键问题】为此，使用RADARSAT-2雷达影像与同时期野外高光谱实测数据结合的方式，探索2种遥感数据源与土壤氮磷之间潜在的关系，以期提高土壤氮磷监测的精度。

1 研究区概况

研究区位于内蒙古河套灌区解放闸灌域境内（40°12′—41°20′N，106°10′—109°30′E）。灌区地形平坦，西南高，东北低，海拔1007～1050 m，坡度0.125‰～0.2‰。属典型的温带大陆性气候，夏季高温干旱、冬季严寒少雪，年降水量100～250mm，蒸发量2400 mm 左右，无霜期120～150 d。河套灌区东西长270 km，南北宽40～75 km，总面积105.33 万余hm2。解放闸灌域位于内蒙古河套灌区西部，南临黄河，北靠阴山，东与永济灌域毗邻，西与乌兰布和沙漠接壤，灌域控制面积22 万hm2，现有灌溉面积14万hm2，排域面积17 万hm2，涉及杭锦后旗、临河区、磴口县、乌拉特中旗、乌拉特后旗共15个镇和1个国营农场的农牧业生产区域。河套灌区降水量小、蒸发量大，土壤养分量地域分布极不均匀，土壤次生盐渍化严重，导致土壤养分淋失、侵蚀流失及气态化损失较大，加之灌区作物对养分吸收与利用的高依赖性，精确迅速地获知灌区春灌（4—6月）前土壤氮磷分布状况对灌区作物选种育种及确定种植和改良土壤方式有着重要意义。

图1 研究区采样空间分布图Fig.1 Sampling diagram of study area

2 数据处理与分析

2.1 数据获取与处理

于2017年4月6—10日完成土壤采样与调研工作，该期为非种植期，地面植被覆盖较少，采样、光谱收集以及卫星影像所测均为裸露土壤。样点布设综合考虑了研究区种植结构、不同土地利用方式、灌溉方式、地下水位等情况，取样土壤深度为表层0～10 cm，99个采样点分布见图1。

野外采样时通过手持GPS 接收机测量采样点WGS84 坐标系统下的经纬度坐标，获取采样点地表粗糙度等数据，把每个土壤样本编号装入铝盒，在距采样点5 m的两侧位置，以同样的方式进行2次重复采样，用于精度校验。使用半微量凯氏法测定土壤全氮，高氯酸-硫酸法、钼锑抗比色法测定全磷量。野外采集数据时，试验区均未进行灌溉或降雨。

本研究所使用的RADARSAT-2 影像数据，为C波段合成孔径雷达系统。购置同时期RADARSAT-2精细四极化SLC（single look complex）格式雷达影像一景（地面分辨率8 m，相幅25 km×25 km）。使用ENVI 5.3 软件的SAR Scape 模块对原始SLC 影像进行多视、滤波、地理编码及辐射定标等图像预处理措施，在Google 影像上选取GCP（ground control point）对影像进行几何配准后，利用Output ROIs to ASCII功能将所需的采样点后向散射系数输出成文本文档以供使用。

同时，使用美国Analytical Spectral Device（ASD）公司生产的Field Spec4 便携式光谱仪进行采样点土壤反射率光谱数据采集，该光谱仪可测波段范围为350～2500nm，其中，350～1000nm 波段的光谱分辨率为3nm，光谱采样间隔为1.4nm；1000～2500nm 波段的光谱分辨率为10nm，光谱采样间隔为2nm。测量时要保证天气晴朗，光照稳定，且避开阴影和邻近的运动物体；在测试前调试和检验光谱仪与参考板，以保证测试结果的准确性；测试前先用白板定标，白板应小于10°水平放置；光谱仪探头距离待测土壤的高度统一固定在1m 处，每个采样点测量25次，算术平均后得到该采样点土壤的实际反射光谱数据；由于导数光谱有助于限制低频噪声的影响，可在一定程度上提高土壤氮磷量和光谱反射率的相关性。通过光谱仪自带的ViewSpec Pro Version 6.0将原始光谱反射率R进行了3种不同形式的变换：对数（lgR）、一阶导数（R'）、二阶导数（R''），并将上述3种变换光谱及土壤反射率作为本研究的输入光谱，选择和氮磷相关性最大的波段作为特征波段，取其值和SAR 后向散射系数值一同与氮磷回归。

2.2 相关分析法

通过对土壤氮磷量与土壤光谱反射率及其3种变换形式各数据组进行相关分析，得到各数据组与土壤氮磷量之间的相关系数（R）。选取相关系数较大的数据变换形式，进行多元逐步线性回归分析，分别建立土壤光谱反射率及其一阶导数、二阶导数及对数与土壤氮磷量之间的回归方程。计算与分析使用Matlab 软件。

2.3 小波分析法

小波变换可以通过伸缩平移运算对信号逐步进行多尺度细化，最终达到高频处对时间细分，低频处对频率细分的目的，在局部可做到更精确地描述与分离[13]。小波分解可以生成多层包含影响土壤高光谱反射率因素的特征图谱，其中低频系数有利于提取可以决定整个光谱形状的原始光谱的吸收特征，高频系数则更有利于提取原始光谱的噪声及微小的吸收特征[14-15]，重构后的低频分量能够一定程度上剔除由光谱仪精度、测试条件等不确定因素影响的高频噪音，而重构后的高频分量则可以突显光谱信息特征变化明显、波动较为剧烈的细节信息。计算部分使用Matlab 软件自带小波分析工具箱完成。

2.4 BP 神经网络

BP 神经网络由大量神经元构成，即输入层、输出层和隐藏层，其结构的关键并非神经元，而是神经元之间的连接线，每条连接线通过训练被赋予不同的权值[16]。它的学习规则是使用最速下降法，利用梯度搜索技术，经过反复学习训练，以期使网络的实际输出值和期望输出值的误差均方差为最小。线性分析及神经网络计算使用Matlab 软件及其工具箱编程实现。

3 模型构建

3.1 土壤氮磷量与光谱反射率的相关性分析

图2（a）显示了反射率及其3种变化形式和全氮量的相关关系在不同波段上的变化。全氮量与土壤原始反射率在整个光谱波段负相关；与对数变换光谱在350～1850nm 基本正相关，在1850～2450nm 负相关；与一阶、二阶导数变换光谱的相关系数正负交叉，波动剧烈，介于0.5 与-0.5之间，其峰值和相关性大幅提升。4种不同形式的光谱与全氮量的最大相关系数值分别达到0.468、-0.426、0.412，都出自一阶、二阶导数变换光谱。这表明，一阶、二阶导数变换光谱相较原始反射率与对数变换光谱更适合土壤全氮量的敏感波段筛选，同时，最大相关系数绝对值均小于0.50，相关性均不明显，这也说明利用光谱反射率与其3种不同形式的变换光谱通过相关分析难以提取有关土壤全氮量的敏感波段。图2（b）所示全磷量的相关关系与图2（a）相似的规律，不同之处在于光谱反射率在660～1000nm 呈稳定的负相关，相关系数约为-0.18；一阶导数变换光谱在1660～1690nm 与2070～2100nm 波段都有相关系数绝对值大于0.5的敏感反应。

图2 土壤反射率及其不同变换形式与土壤氮磷量的相关性分析Fig.2 Correlation analysis of soil reflectance and its different transformation forms with Nitrogen and Phosphorus in soil

3.2 土壤氮磷模型拟合

分别对土壤氮磷量与光谱的4种形式之间进行相关性分析，将各形式中与对应氮磷量相关系数较大的敏感波段和土壤氮磷量分别作为预测土壤氮磷的自变量和因变量，多元线性回归模型和BP 神经网络模型拟合所得的土壤氮磷实测值与预测值的对比分析结果如表1所示。从表1可以看出，反射率一阶导数变换的预测值与实测值之间的拟合系数综合高于其他变换，在氮磷两种模型中分别为0.346 与0.349，其次为对数变换的0.332 与0.228，光谱反射率与二阶导数变换的拟合度最低。通过多元线性回归的拟合系数均小于0.5，结果都不理想，加入神经网络后结果有所提升，但依然欠佳。

3.3 小波变换后的模型拟合

为提取具有与土壤氮磷量更高相关性的光谱信息，通过Matlab 中小波分析模块对反射率特异值进行了适当剔除，采用Biorthogonal 小波系的Bior1.3函数分别对上文中综合拟合最优的反射率一阶导数光谱进行1～8层的小波变换。图3（a）为与土壤全氮相关系数最高的1层低频重构光谱、5层和6层高频重构光谱的相关系数。1层低频重构光谱（CA1）在剔除了部分噪音的基础上，和原输入光谱（反射率一阶导数光谱）的相关性较高，尤其是在1480～1489、1514～1523、2050～2059 和2314～2323nm 4个波段，5层高频重构光谱（CD5）在419～434、1123～1154、1987～2010nm 3个波段，6层高频重构光谱（CD6）在1143～1473nm 1个波段均表现出与土壤全氮量大于0.4的相关，这些波段与未进行小波分析的原输入光谱的特征波段近似但却不一致，说明小波变换在压缩波长自变量的同时，增加了各组分之间的区别，可以显著提高模型的反演精度。与土壤全氮不同，土壤全磷与2层和3层低频重构光谱、6层高频重构光谱的相关系数最高，如图3（b）所示。2层与3层低频重构光谱相关曲线几近重合，体现出较高的相关性，其中，2层低频重构光谱在1679～1688、2133～2142、2319～2328nm 3个波段，3层低频重构光谱在2083～2092、2131～2140、2321～2330nm3个波段以及 6层高频重构光谱在 1443～1452、2275～2284、2371～2380nm 3个波段均表现出与土壤全磷量大于0.42的相关，尤其是2层低频重构光谱的第3个特征波段相关系数超过了0.5。

表1 土壤氮磷多元回归及神经网络分析结果Table1 Test results of MLR model and neural network for Nitrogen and Phosphorus contents of soil

图3 小波分解重构分量与土壤氮磷量的相关性分析Fig.3 Analysis on correlation of between wavelet decomposition and reconstruction and Nitrogen and Phosphorus contents of soil

图4 土壤全氮量反演值与实测值拟合分析Fig.4 Fitting analysis of inversion value and measured value of soil total nitrogen content

通过选取筛选后1层低频重构光谱、5层和6层 高频重构光谱中的特征波段，结合微波雷达RADARSAT-2 四极化后向散射系数VV、VH、HH、HV及其2种组合HH/VV、HV/VH以及地表粗糙度ZS建立线性回归模型，结果如式（1）所示。其模型反演与实测值决定系数R2为 0.46，均方差RMSE=0.1729，如图4（a）所示。将线性回归模型中的参量带入神经网络中，其模型反演与实测值决定系数R2为0.7496，均方差RMSE=0.1102，如图4（b）所示。

式中：X1480、X1514、X2314分别为1层低频重构光谱1480～1489、1514～1523、2314～2323nm的波段的算术平均数；X1123、X1987为 5层高频重构光谱1123～1132、1987～1996nm的波段的算术平均数；X1443为六层高频重构光谱1443～1452nm的波段的算术平均数，ZS为地表组合粗糙度。

通过选取筛选后2层和3层低频重构光谱、6层高频重构光谱中的特征波段，结合微波雷达及地表参数建立线性回归模型，结果如式2所示。其模型反演与实测值决定系数R2为0.3954，均方差RMSE=0.1921，如图5（a）所示。将线性回归模型中的参量带入神经网络中，其模型反演与实测值决定系数R2为0.7592，均方差RMSE=0.1102，如图5（b）所示。

式中：X1679、X2133、X2319分别为2层低频重构光谱1679～1688、2133～2142、2319～2328nm的波段的算术平均数；X1443、X2275、X2371分别为6层高频重构光谱1443～1452、2275～2284、2371～2380nm的波段的算术平均数，ZS为地表组合粗糙度。

图5 土壤全磷量反演值与实测值拟合分析Fig.5 Fitting analysis of inversion value and measured value of total phosphorus content in soil

图6 土壤全氮与土壤全磷量反演预测结果Fig.6 Inversion results of Nitrogen (a) and Phosphorus (b) content of soil

3.4 土壤氮磷的预测

图6所示为小波重构特征波段和雷达参数建立的经验回归模型模拟的氮磷分布，先将特征光谱波段中心反射率一阶导数和地表粗糙度用ArcGIS 软件地统计模块克立格插值生成栅格文件，同时提取生成SAR 后向散射系数栅格文件，最后基于上述式（1）和式（2），采用ENVI 遥感图像软件进行栅格文件波段运算获取氮磷空间分布。从结果看出研究区西部养分贫瘠，其主要原因在于研究区西部主要为土壤盐碱化严重地区和多年撂荒地，土地未得到合理的治理和使用，这正是今后规划土地利用和土壤研究治理关注的重点。

4 讨论

联合使用地面实测光谱和SAR 后向散射特性进行土壤氮磷反演试验。全氮量与土壤原始反射率在整个光谱波段负相关；与对数变换光谱在350～1 850 nm基本正相关，在1 850～2 450 nm 负相关；与一阶、二阶导数变换光谱的相关系数正负交叉，波动剧烈，介于0.5 与-0.5之间，这与万余庆等[17]的研究一致。全磷量与光谱反射率在660～1 000 nm 呈稳定的负相关，相关系数约为-0.18；一阶导数变换光谱在1 660～1 690 nm 与2 070～2 100 nm 波段都有相关系数绝对值大于0.5的敏感反应，这与徐丽华等[18]研究一致。分析中对光谱值进行了不同形式的数学变化，发现光谱反射率一阶导数与土壤氮磷的相关性优于原始数据及其他变换，这与王静等[19]研究一致。在进行土壤氮磷模型拟合时分别使用了一元线性回归、多元线性回归和BP 神经网络等不同的拟合分析方法，通过BP 神经网络所取得的拟合结果较其他方法更佳，BP 神经网络表现出在复杂非线性函数逼近方面的优越性，这与吴善玉等[20]研究一致。通过对影像原始数据进行小波变换，从而获得的影像数据的有效融合，可以实现多源遥感数据的特征互补，这与杨晓蕊等[21]研究成果一致。本文更进一步的发现高光谱数据低频分量的1～3层、高频分量的4～6层经小波变换后的相关性较其他各层效果更佳。土壤全氮量与1层低频重构光谱（CA1）在1 480～1 489、1 514～1 523、2 050～2 059 nm和2 314～2 323 nm 4个波段，5层高频重构光谱（CD5）在419～434、1 123～1 154 nm 和1 987～2 010 nm 3个波段，6层高频重构光谱（CD6）在1 143～1 473 nm 1个波段均表现出大于0.4的相关。土壤全磷与2层和3层低频重构光谱、6层高频重构光谱的相关系数最高，2层与3层低频重构光谱相关性较高，其中，2层低频重构光谱在1 679～1 688、2 133～2 142、2 319～ 2 328 nm 3个波段，3层低频重构光谱在2 083～2 092、2 131～2 140、2 321～2 330 nm 3个波段以及6层高频重构光谱在1 443～1 452、2 275～2 284、2 371～2 380 nm 3个波段均表现出与土壤全磷量大于0.42的相关，尤其是2层低频重构光谱的第3个特征波段相关系数超过了0.5。这些波段与未进行小波分析的原输入光谱的特征波段近似但却不一致，说明小波变换在压缩波长自变量的同时，增加了各组分之间的区别，可以显著提高模型的反演精度，这与莫才健等[22]研究一致。

本文拟合的经验模型依赖于地表实验数据且具有一定的地域限制性。如何从试验方案的制定，依据高光谱成像理论出发，进一步寻找特征波段与土壤氮磷之间的联系，以完善土壤氮磷遥感监测模型的准确性和普适性是今后研究的重点。

5 结论

1）不同数学形式变换可以提高光谱值与土壤氮磷量数据的相关性，其中光谱反射率一阶导数变换模型的拟合相关系数高于其他变换。

2）对高光谱数据的1～8层进行小波分解与重构，可有效提高光谱反射率及其3种变换值与土壤氮磷的相关性，降低外界环境对光谱值引起的干扰，尤其是对低频分量的1～3层、高频分量的4～6层的效果更好。

3）筛选小波变换后低频重构光谱和高频重构光谱中的氮磷反演特征波段，结合微波雷达四极化后向散射系数及其组合值，考虑地表粗糙度因子，所建立的氮磷多源遥感神经网络模型，其决定系数均超过0.74，总均方差RMSE均为0.11。

4）BP 神经网络可以有效逼近光谱数据和微波雷达数据组合与土壤氮磷之间复杂的非线性关系，从而克服线性回归拟合度不高的局限。这既证实了采用非传统、非暗室条件下野外测定土壤高光谱的可行性，同时也为河套灌区土壤的监测提出了新方案，可以为日后灌区的土壤改良起到一定的指导意义。若进一步减少野外光谱测定时环境的干扰因素，将有利于该模型的完善与实践的过程简化。