基于多源遥感协同反演的区域性土壤盐渍化监测

冯雪力 刘全明

(1.内蒙古建筑职业技术学院，呼和浩特 010070； 2.内蒙古农业大学土木建筑工程学院，呼和浩特 010018)

0 引言

土壤盐渍化是干旱、半干旱农业灌区面临的主要土地退化问题。内蒙古河套灌区是西北地区典型的寒旱灌区，传统灌溉多为含盐黄河水漫灌，且受降水量少、蒸发量大、排水排盐不畅、地势平坦而地下水位居高等自然客观条件限制，造成灌区盐渍化和生态环境恶化，从而影响了农作物生长和粮食生产安全。因此，精准快速地获取盐渍土盐分时空分布信息，对灌区的盐渍土防治、中低产田改造与农业可持续发展具有重大的指导作用。传统野外土壤定点调查盐渍化方式不仅耗费时力、破坏性强，而且测点少，难以反映区域尺度的盐渍化分布，而多传感器、多波段、多时相的遥感技术为大面积实时动态监测盐渍土状况提供了可能。

国外从20 世纪70 年代[1-4]就开始研究盐碱土遥感监测，国内从20世纪80 年代开始集中在利用可见光波段进行遥感盐分反演。盐渍化土壤光谱反射率和土壤参数间的转换函数是复杂非线性的关系，故人工神经网络逐渐被更多地用于获取土壤盐渍化参数[5-9]。吐尔逊·艾山[10]以新疆渭干河-库车河三角洲绿洲作为研究区，利用实测光谱数据和地下水埋深、矿化度等因子构建了基于BP人工神经网络(Back propagation artificial neural networks，BPANN)的土壤盐分反演模型，模型精度为80.77%。王静等[11]以松嫩平原西部长岭县为例，利用光谱导数变换，获取了表征盐渍土盐分信息的最佳波段，所建立的3 层神经网络盐渍土盐分模型预测结果优于传统回归模型。丁建丽等[12]研究发现，实测高光谱土壤、植被一阶微分光谱变换对土壤盐渍化响应度最敏感，基于实测综合光谱指数的盐渍化监测高光谱模型可以准确提取土壤盐渍化信息，明显优于传统遥感方法中单纯利用植被指数或者土壤盐分指数的模型。姚远等[13]以新疆渭干河-库车河三角洲绿洲不同盐渍化程度的土壤为研究对象，以电磁感应仪(EM38)测得的盐渍土电导率数据和ASD光谱仪测得的盐渍土高光谱数据为基础数据源。王爽等[14]以新疆渭干河-库车河三角洲绿洲区域范围内不同盐渍化程度的土壤高光谱反射率及其土壤含盐量为基础数据源，构建了地表光谱建模的最佳土壤盐渍化监测模型。刘全明等[15]以受盐渍化影响较严重的内蒙古河套灌区解放闸灌域为试验区，构建了BPANN土壤含盐量的定量反演模型。彭杰等[16]通过分析新疆、浙江、吉林3 个不同地区盐渍化土壤的高光谱特征，研究了盐渍化土壤高光谱特征的区域异质性，并构建高精度的跨区域土壤盐分高光谱定量反演模型。彭杰等[17]以南疆地区水稻土为研究对象，通过分析土样的高光谱数据和室内测定的盐分与电导率数据，研究了耕作土壤含盐量与电导率的关系，并比较了含盐量和电导率与不同光谱指标的相关性以及二者高光谱反演的精度。研究者们使用高光谱和电磁感应数据进行土壤盐分模拟，获得了不少成果[18-22]。

在前期研究基础上，本文优选高光谱土壤盐分特征波段，融合微波后向散射特性，构建基于多源遥感数据的土壤盐分预测模型，为利用多源遥感数据协同反演盐渍土时空分布提供科学基础。为此，在高光谱盐分特征波段优选方法对比的基础上，创建基于多源遥感数据的土壤盐分经验模型和人工智能非线性预测模型，从而将雷达后向散射系数、光谱反射率及其变换值等参数转换为土壤盐分，以提高盐渍化监测的精度与广适性。

1 研究区概况

内蒙古河套灌区地处内蒙古自治区的巴彦淖尔市，位于40°12′～41°20′N，106°10′～109°30′E。灌区地形平坦，西南高，东北低，海拔1 007～1 050 m，土壤以盐渍化浅色草甸土和盐土为主。夏季高温干旱、冬季严寒少雪，年降水量130～215 mm，蒸发量高达2 100～2 300 mm，无霜期120～150 d，封冻期长，是典型的温带大陆性气候。河套灌区东西长250 km，南北宽40～60 km，土地总面积18 890 km2，其中灌溉面积5 583 km2。河套灌区土壤养分含量地域分布极不均匀，降水量小、蒸发量大、土壤次生盐渍化严重，快速了解灌区夏灌(4—6月)前土壤盐分分布状况对灌区灌溉制度和土壤改良有着深刻影响和重要意义。

2 材料与方法

2.1 数据获取与处理

试验区选在河套灌区内的解放闸灌域。为保证野外采样时间与雷达成像时间对应一致，提前购置2016年4月6日RADARSAT-2精细四极化SLC(Single look complex)格式雷达影像一景，影像尺寸为25 km×25 km，其地面分辨率为8 m。考虑到土壤盐渍化空间分布的不均匀性，在野外采样前使用Google Earth对灌域的地形地貌进行分析，结合前期土壤盐分统计资料并制定野外采样路线，89个采样点分布如图1红色五角星所示。野外通过手持GPS接收机测量采样点在WGS84坐标系下的经纬度坐标，获取采样点地表粗糙度等数据，并同时在每个样点采集3份土壤表层土，取平均值以消除取样代表性误差，把每个土壤样本编号装入铝盒，采用干燥法测定土壤水溶性盐分含量。同时使用带有厘米格网的剖面板方法来测量地表粗糙度。

图1 试验区雷达影像图Fig.1 Radar image of experimental area

试验区所使用的RADARSAT-2影像数据是由C波段合成孔径雷达系统获取，其影像处理主要包括辐射定标、几何校正、斜距转地距和滤波等。本研究使用ENVI软件的SAR Scape模块对原始SLC影像进行多视、滤波、地理编码和辐射定标等处理，在Google影像上选取地面控制点(Ground control point,GCP)对影像进行几何配准后，得到标准四极化后向散射系数影像数据，表1所示为试验区部分采样点四极化后向散射系数与土样全盐量、地表组合粗糙度。表中SHH、SHV、SVH、SVV分别代表HH(水平极化方式发射，水平极化方式接收)、HV(水平极化方式发射，垂直极化方式接收)、VH(垂直极化方式发射，水平极化方式接收)、VV(垂直极化方式发射，垂直极化方式接收)4种极化方式的后向散射系数。

表1 雷达影像后向散射系数与土壤全盐量原始数据Tab.1 Data of radar back-scattering coefficients and soil saltcontents

野外采样点土壤反射率光谱数据的采集使用美国Analytical Spectral Device公司生产的Field Spec4便携式光谱仪，可测量350～2 500 nm波段，其中，350～1 000 nm波段的光谱采样间隔为1.4 nm，光谱分辨率为3 nm；1 000～2 500 nm波段的光谱采样间隔为2 nm，光谱分辨率为10 nm。利用光谱仪自带的ViewSpec Pro Version 6.0软件将原始光谱反射率R进行了3种不同形式的变换：对数(lgR)、一阶导数(R′)、二阶导数(R″)，并将光谱反射率及上述3种变换光谱作为本研究的输入光谱数据，后期分别将上述4种输入值与对应采样点盐分值进行相关性分析。

2.2 主成分回归

主成分回归是考察多个变量间相关性的一种多元统计方法，研究如何通过少数几个主成分来揭示多个变量间的内部结构，即从原始变量中导出少数几个主成分，使它们尽可能多地保留原始变量的信息，且彼此间互不相关。通常数学上的处理就是将原来多个指标作线性组合，作为新的综合指标。

2.3 逐步回归分析

逐步回归分析是将变量逐个引入模型，每引入一个解释变量后都要进行F检验，并对已经选入的解释变量逐个进行t检验，当原来引入的解释变量由于后面解释变量的引入变得不再显著时，则将其删除。以确保每次引入新的变量之前回归方程中只包含显著性变量。这是一个反复的过程，直到既没有显著的解释变量选入回归方程，也没有不显著的解释变量从回归方程中剔除为止，以保证最后所得到的解释变量集是最优的。

2.4 偏最小二乘回归

偏最小二乘回归是一种多因变量对多自变量的回归建模方法，可以较好地解决许多用普通多元回归无法解决的问题。和PCR一样，两种方法都采用得分因子作为原始预测变量线性组合的依据，用于建立预测模型的得分因子之间必须线性无关。而PLSR与PCR的不同之处在于得分因子的提取方法不同，PCR产生的权重矩阵反映的是预测变量之间的协方差，PLSR产生的权重矩阵反映的是预测变量与响应变量之间的协方差。在建模当中，PLSR产生了权重矩阵，矩阵的列向量用于计算变量列向量的得分矩阵。通过不断地迭代计算这些权重，使得响应与其相应的得分因子之间的协方差达到最大。PLSR是使用波段的权重来表示波段对于土壤实测值的敏感性，进而选取对盐分敏感的特征波段，并进行土壤盐分的回归拟合。

2.5 BP人工神经网络

BP人工神经网络能够辨识复杂系统输入和输出数据组间的非线性关系，不需任何假定或构建数学方程，可在一定程度上消除特异值的影响，具备强大的适应性、并行处理及抗干扰能力和突出的非线性拟合的特点，现已被各学科领域广泛应用。本文经综合优化试验，确定选用3 层结构的BP模型，该模型由输入层、隐含层和输出层组成。考虑土壤含盐量与高光谱光学特征波段及SLC影像的四极化后向散射系数、地表粗糙度有着显著的响应关系，最终确定输入层由11个神经元组成，包括4个光学特征波段反射率，6个后向散射系数及组合值SHH、SVV、SHV、SVH、SHH/SVV、SHV/SVH和地表组合粗糙度；隐含层神经元数的确定由试算法完成，本文经试算法确定隐含层神经元为10 个；输出层神经元为1 个，对应为采样点土样的全盐量。计算采用Matlab自带神经网络工具箱(nprtool)，采用原始数据的70%参与学习与训练网络模型，其余30%作为验证。

3 特征波段选取建模与土壤盐分预测

3.1 特征波段选取建模

图2反映了土壤全盐量与反射率的不同变换形式的相关关系。不难发现，相较于原始光谱和对数变换，光谱的一、二阶导数可以获得更好的相关性。在土壤全盐量的相关性中，在一阶导数变换的1 088～1 092 nm、1 806～1 810 nm、2 138～2 142 nm、2 289～2 295 nm 这4个波段表现最佳，相关系数分别为0.29、0.36、0.3、0.27，在二阶导数变换的618～622 nm、1 802～1 806 nm、2 169～2 173 nm、2 344～2 348 nm这4个波段相关性最高，相关系数分别为0.37、0.28、0.39、0.27。

图2 反射率的不同变换形式与土壤全盐量的相关性分析Fig.2 Analysis on correlation of reflectance transforms and soil salinity contents

利用PCR对光谱的一、二阶导数中的特征波段进行主成分分析，结果如图3所示。可以发现，无论是一阶还是二阶导数变换的特征波段，都可以从中选取3特征波段来涵括，累积贡献率分别为99.31%和99.87%；其中，第一主成分一阶导数变换特征波段1 088～1 092 nm与二阶导数变换特征波段618～622 nm的特征值和贡献因子分别为4.07×10-7、1.51×10-7和72.19、74.08，皆具有较高代表性。

图3 特征波段的主成分分析Fig.3 Principal component analysis of characteristic bands

利用MSR对光谱的一、二阶导数中的特征波段进行逐步回归分析，结果如图4所示。其中图4a中X1～X4分别表示1 088～1 092 nm、1 806～1 810 nm、2 138～2 142 nm和2 289～2 295 nm波段，图4b中X1～X4分别表示618～622 nm、1 802～1 806 nm、2 169～2 173 nm、2 344～2 348 nm，研究表明，光谱一阶导数的1 806～1 810 nm、2 138～2 142 nm 2个特征波段和二阶导数的618～622 nm、1 802～1 806 nm、2 344～2 348 nm 3个特征波段更具有代表性。

图4 特征波段的多元逐步回归分析Fig.4 MSR analysis of characteristic bands

再使用PLSR选取特征波段并拟合土壤全盐量，PLSR使用波段的权重来表示波段对于土壤实测值的敏感性。如图5所示，PLSR筛选的波相较PCR以及MSR选取的波段延后，且存在极为敏感的波段，如一阶导数变换中1 850～1 854 nm波段权重可达4，二阶导数变换中1 853～1 857 nm与2 483～2 487 nm权重分别可达15和5.7，高权重波段与上文中分析的高相关性波段也均有延后现象。

图5 反射率导数变换与土壤全盐量的相关性及PLSR权重分析Fig.5 Analysis on correlation and PLSR weight between reflectance derivative transformation and soil salinity

PCR和MSR的拟合结果对比如表2所示。不难发现，二阶导数变换的拟合程度高于一阶导数变换，均方根误差也都有所下降，说明二阶导数变换更适合于土壤全盐量的拟合；综合比较使用光谱二阶导数变换特征波段的MSR经验模型模拟全盐量效果最优，其决定系数(R2)和均方根误差(RMSE)分别为0.343和9.414 g/kg。

表2 多模型拟合结果对比Tab.2 Results contrast of multi-model

图6 土壤全盐量实测值与模拟值Fig.6 Measured and simulated values of soil salinity

使用MSR选取的3个特征中心波段反射率、6个SAR后向散射系数及其组合值和地表组合粗糙度，建立土壤全盐量与光谱二阶导数的经验模型

式中X618、X1 802、X2 346——618～622 nm、1 802～1 806 nm、2 344～2 348 nm波段的反射率二阶导数算术平均值

同时建立神经网络模型，由图6可见，神经网络模型R2为0.890 8，远优于MSR模型，体现了神经网络模型的优越性。

3.2 土壤盐分预测

根据表3的河套灌区土壤盐渍化分类标准，将采样点所覆盖的试验区域分为5类盐渍化类型，即非盐渍土、轻度盐渍土、中度盐渍土、强盐渍土和盐碱土。图7所示为MSR选取的光谱二阶导数的3个特征波段和雷达参数建立的回归方程模型模拟的全盐量预测结果，先将MSR提取的3个特征光谱波段中心反射率二阶导数和地表组合粗糙度用ArcGIS软件地质统计模块克立格插值生成5个栅格文件，同时提取生成6个SAR后向散射系数栅格文件，最后基于光谱二阶导数的回归方程，采用ENVI遥感图像软件进行11个栅格文件运算，其模拟结果的统计分析见表4，可见受到盐渍化影响的面积占96.23%，盐碱土占21.47%，应对重点区域采取科学的盐碱土改良措施进行治理，以保证农业生产的可持续发展。

表3 土壤盐渍化分类Tab.3 Classification of soil salinity

注：Ⅰ～Ⅴ级分别代表非盐渍土、弱盐渍土、中盐渍土、强盐渍土和盐碱土。

图7 研究区土壤全盐量预测图Fig.7 Prediction image of soil salinity in study area

表4 土壤盐渍化预测分类占比Tab.4 Proportion of all levels of soil salinity prediction classification

4 结论

(1)通过不同形式的光谱变换处理可以使光谱原反射率与土壤盐分含量的关系得以提升。分别利用PCR、MSR和PLSR选取特征光学波段，发现光谱的一、二阶导数具有更好的相关性，尤其是二阶导数变换的618～622 nm、1 802～1 806 nm、2 169～2 173 nm、2 344～2 348 nm这4个波段相关性较强；PLSR筛选的波段相较相关系数法偏后，其二阶导数变换模型拟合度小于MSR。

(2)通过协同光谱特征波段反射率二阶导数、雷达后向散射系数和地表组合粗糙度诸影响因子，对比土壤盐分的PCR、MSR和PLSR模型，其中MSR模型的决定系数和均方根误差分别为0.343和9.414 g/kg，优于另两种模型，而BP人工神经网络模型为最佳预测模型，其模型预测R2为0.890 8，模型的稳定性和预测精度较好。融合高光谱特征波段与SAR后向散射系数及地表组合粗糙度参数而构建的多源遥感数据源神经网络模型具有明显的优势。