基于实测高光谱和Sentinel-2B 影像的银川平原土壤盐分反演

毛鸿欣，贾科利，张 旭

(宁夏大学 地理科学与规划学院，宁夏 银川 750021)

土壤盐渍化已成为中国生态文明建设和全球变化议题中的重要研究内容[1]. 盐渍土肥力低下，影响作物生长，阻碍农业经济可持续发展. 尤其在干旱、半干旱灌溉农业区，长期灌溉引起土壤次生盐渍化，对区域的农业生产和生态环境造成巨大影响，具有严重的生态环境风险，因而盐渍化监测和治理受到各界普遍关注.

近年来，国内外学者利用遥感技术对土壤盐渍化监测和预测开展了较多研究. 如陈俊英等[2]将无人机航拍影像和高分1 号卫星影像用于盐分反演模型构建；再屯古丽·亚库普等[3]引入土壤水分和pH 因子，利用雷达遥感技术实现了土壤盐分预测；厉彦玲等[4]建立了高光谱和多光谱影像的机器学习和统计回归模型；陈实等[5]和安德玉等[6]通过构建土壤指数定量探究了盐分空间格局；吾木提·艾山江等[7]将实测高光谱与WorldView-2 影像进行光谱匹配，提升了土壤盐分的预测精度和制图精度；马驰[8-9]基于Sentinel-1 和Sentinel-2B 影像，探究了哨兵影像反演土壤盐分的可行性，发现效果较好. 为提高模型精度，张贤龙等[10]基于15 种光谱变换构建高光谱指数，能显著提升建模精度；Rocha等[11]利用PCA、一阶微分、和原始反射光谱作为极限学习机（Extreme Learning Machine，ELM）的输入参数，有助于提升模型的估算性能；张智韬等[12]和王海峰等[13]分别利用岭回归和灰度关联-岭回归，可有效消除光谱变量间的共线性，提高反演精度. 此外，众多学者对于融合多光谱影像和实测高光谱数据反演土壤盐分做了较多研究，其一是寻找高-多光谱数据间的敏感参量，进而实现不同尺度数据间的光谱匹配. 如王爽等[1]借助盐分指数SI3实现了实测高光谱到多光谱影像的光谱匹配；姚远等[14]利用SI2指数实现了TM 影像和高光谱的光谱匹配. 另一种方法是利用重采样技术，通过传感器光谱响应函数，将实测光谱的连续波谱重采样至影像的对应波段上，从而进行分析和建模. 如梁静等[15]将美国ASD（Analytical Spectral Devices）公司生产的ASD FieldSpec 3 光谱波段重采样至OLI（Operational Land Imager）影像的7 个波段上；阿尔达克·克里木等[16]利用重采样技术将ASD 与ASTER高光谱影像进行匹配. 但这种方法也存在一定缺点，其本质为利用加权平均方法使实测高光谱数据与多光谱影像的中心波长一一对应，结果难免会使原始高光谱数据失真. 且这些成果多以中空间分辨率影像作为数据源来反演土壤盐分. 以Sentinel-2B影像，并融合实测高光谱数据反演土壤盐分的研究较少，且利用岭回归方法构建土壤盐分反演模型的成果并不多见.

因此，本文以银川平原为研究区，利用数理统计回归实现实测高光谱与Sentinel-2B 影像间的地-星光谱匹配，构建土壤盐分最优预测模型，实现银川平原土壤盐分的定量反演，其结果可为光谱匹配和不同尺度遥感数据联合反演土壤盐分提供理论借鉴和参考.

1 材料与方法

1.1 研究区概况银川平原位于宁夏回族自治区北部，地理坐标为105°45′～106°56′E，37°46′～39°23′N，南起青铜峡冲积扇口，北达石嘴山麻黄沟峡谷，西至贺兰山脊宁蒙分界，东达鄂尔多斯台地西南缘（图1（a））[17]. 银川平原由贺兰山洪积扇和黄河冲积平原两大地貌单元构成，东西宽42～60 km，南北长165 km[18-19]，面积7 615 km2. 地势自西向东、由南向北微倾，坡降比为0.6‰～1‰，形成一个巨大的“凹”字型地貌（图1（b））. 该区属温带半干旱大陆性气候，干旱少雨，蒸发强烈，年均降水量180～200 mm. 但湖泊密布，灌溉渠系纵横，成为水稻、小麦等作物的主产区. 长期的引黄灌溉和地下水埋深，加之强烈的地表蒸发，使盐渍化土壤分布广泛.

1.2 数据获取与处理

1.2.1 光谱数据采集与处理 基于ArcGIS 软件，大致以5 km 等间隔网格布设样点（图1（c）），利用美国ASD 公司生产的FieldSpec 4 野外便携式地物光谱仪进行光谱数据测定，其光谱范围为350～2 500 nm，重采样间隔为1 nm. 土壤光谱反射率采集于2019 年5 月12—22 日10:00—14:30 进行，期间晴朗无风，光线充足. 采集时先用白板进行校正，选取地表面平坦光滑，远离植被阴影的区域，将光纤探头垂直对准地表进行光谱数据测定. 为使数据准确且具有代表性，每个样点采集5 条光谱曲线取其均值作为该样点的光谱数据. 异常光谱曲线去除和光谱均值输出均在ViewSpecPro 软件中进行，并利用Origin 2018 和ENVI 5.5 软件进行光谱数据平滑去噪和连续统去除（Continuum Removal，CR）预处理.

图1 银川平原位置、地形和采样点分布Fig. 1 Location，terrain，and distribution of sampling points of Yinchuan Plain

1.2.2 土壤样本采集与处理 光谱采集后同点采集0～20 cm 表层土装入自封袋中密封带回，利用手持GPS 记录每个样点的经纬度坐标、采样日期、土样与光谱的对应编号及样点周围环境信息.再将采集的土样剔除砖块、碎石、杂草等杂质，经自然风干、研磨、过2 mm（10 目）孔筛[20]，按照土水比1∶5 配置浸提液[21]，利用电导率法[22-23]计算其盐分含量. 经对测定结果及光谱数据进行分析，剔除异常样点后得到167 个样本供试.

参考文献[17]，根据盐分含量进行土壤样本盐渍化程度分级. 银川平原盐渍化土壤分级和描述性统计分析如表1. 样本中盐渍化土壤样本共100 个，约占样本的60%. 轻度盐渍化土壤样本占盐渍化样本比例为34%，盐土样本个数为44 个，中度和重度盐渍化样本相对较少，约占总样本的13%. 样本土壤含盐量总体变异系数超过160%，变异性较强，表明样本分布离散，对土壤盐渍化反演模型建立具有普适性.

表1 银川平原土壤盐渍化程度及样本统计Tab. 1 Soil salinization degree and sample statistics in Yinchuan Plain

1.2.3 影像数据获取与处理 研究使用的Sentinel-2B 影像从欧洲航天局数据网站（https://scihub.copernicus.eu/）获取，影像有13 个波段，光谱范围为400～2 400 nm，基准面为WGS84，地图投影为UTM 投影，空间分辨率为10 m.

选用覆盖研究区的5 景Sentinel-2 L1C 级数据产品，轨道号为61，数据时间为2019 年5 月22 日，与野外采样时间对应. 数据已经过像元级几何精校正和正射校正，为实现大气表观反射率到地表辐射率间的转换，利用ESA 自带插件Sen2cor 进行Sentinel-2B 数据生产和大气校正，于SNAP 软件中进行重采样（重采样后剔除了卷云波段）和文件格式转换，最后在ENVI 5.5 中进行波段合成和无缝镶嵌，利用研究区矢量边界裁剪拼接好的影像，得到研究区影像数据.

1.3 研究方法

1.3.1 光谱变换及光谱指数选取 为减弱和消除背景噪声，突出拐点和最值特征[24-27]，增强吸收和反射特征[28-32]，提高光谱敏感性与相关性，本文基于实测高光谱反射率做了对数、倒数、一阶微分、二阶微分、CR、SG（Savitzky-Golay，SG）等多种光谱变换，并与土壤盐分进行相关性分析.

借助土壤光谱指数探究盐渍土对其上覆植被特征的响应和进行盐渍化信息提取已成为一种重要技术手段[33-35]. 植被长势是土壤盐分很好的指示因子，植被指数高低能在一定程度上反映土壤含盐量的多少. 参考已有研究成果，选取常用光谱指数见表2.

表2 植被指数和土壤盐分指数计算Tab. 2 Calculation of vegetation indices and soil salinity indices

1.3.2 模型构建及精度评价 岭回归模型（Ridge Regression，RR）是以牺牲无偏估计和部分建模精度，获取更为稳定可靠模型的一种改良最小二乘法，因其具有变量自动筛选和共线性消除功能，被广泛应用于光谱模型构建[12-13]. 偏最小二乘回归模型（Partial Least Squares Regression，PLSR）兼具多元逐步线性回归模型（Multivariable linear Stepwise Regression，MLSR）、主成分模型（Principal Component Analysis，PCA）以及简单线性回归的优势，既能依据相关性进行变量自动筛选，又可保证模型稳定性，近年来被广泛应用于高光谱数据建模[36-38]. 本文RR 模型的构建在SPSS 22 中编写代码实现，PLSR模型的构建在SIMCA-P 12.0 软件中进行. 利用决定系数R2、RMSE 评价模型的拟合度，R2越大、RMSE越小则模型预测效果越佳，将样本按含盐量从小到大进行排列，按留二选一的方法划分全部样本，得到110 个建模样本和57 个验证样本.

2 结果分析

2.1 不同程度盐渍化土壤光谱特征分析由图2（a）不同盐渍化程度土壤高光谱反射光谱曲线可知，400～1 300 nm 光谱反射率随着波长增加而逐渐增加，其中在400～750 nm 曲线的斜率最大，形成了反射的“陡坡”. 在1 450 nm 和1 950 nm 附近有2 个明显的吸收谷，为水体的吸收带. 在2 200～2 500 nm，反射率逐渐减小，从不同程度盐渍化土壤光谱反射率看，盐土光谱在整个波段范围都具有最大反射率. 非、轻、中、重盐渍化土壤在450～1 800 nm，盐分含量越高其反射率越大. 在 1 800～2 500 nm，不同程度盐渍化土壤光谱反射率没有明显规律，但随着盐分的增加，反射率有下降的趋势.

利用ArcGIS“Extract Multi values to points”工具获取各样点在Sentinel-2B 影像上逐波段光谱反射率，并按盐分等级进行均值处理，得到多光谱反射率曲线（图2（b））. 由图2（b）可知，在可见光范围内，光谱反射率随波长的增加持续上升，于780 nm附近达到反射顶峰后开始逐渐下降，并于2 200 nm处到达吸收谷，曲线走势与高光谱反射光谱特征一致. 但由于多光谱影像光谱分辨率较小，难以识别地物的细微光谱特征变化，在1 640 nm 后重度盐渍化和盐土曲线下降速度较快，没有呈现出明显规律性变化.

图2 银川平原不同程度盐渍土光谱曲线Fig. 2 Spectral curves of soil with different degrees of salinity in Yinchuan Plain

2.2 土壤盐分与实测高光谱反射率相关性分析与建模对土壤实测光谱反射率做对数、倒数、一阶微分、二阶微分、CR、SG 等变换，与土壤含盐量进行相关性分析. 由图3 可知，R、[lg（R）]′、[lg（1/R）]′、SG 和（SG）0.5与土壤含盐量的最大相关系数绝对值依次为0.632、0.703、0.696、0.629、0.614，相关性较好，所含特征波段较多. 其中[lg（R）]′为最优光谱变换，其最大相关系数达到了0.703，相关性较大的波段为435、500、499、498 nm 和501 nm. 故本文以这些波段为特征波段构建模型.

图3 银川平原高光谱反射率数学变换与土壤盐分的相关性Fig. 3 Correlation between mathematical transformation of hyperspectral reflectance and soil salinity in Yinchuan Plain

基于敏感波段衍生计算高光谱指数，并与含盐量进行相关分析，依据Pearson 相关系数选取特征光谱指数，进行土壤盐分与实测光谱反演模型的构建. 表3 所列出的6 个特征高光谱指数与盐分相关性较好，且均在0.01 置信水平上显著，部分指数相关性较低，文中未详细列出. 为与多光谱指数对应，选用SI1、SI2、SI3、S3、Int1、Int2为建模所用的特征光谱指数.

表3 高光谱指数与含盐量的相关性Tab. 3 The correlation between hyperspectral indices and soil salt content

利用RR 和PLSR 方法，分别基于高光谱特征波段和特征光谱指数进行盐分反演模型构建. 表4表明，基于实测光谱数据构建的RR 和PLSR 盐分反演模型决定系数均通过0.01 置信水平检验，显著性相差不大. 其中，基于特征波段构建的PLSR模型拟合度最高.

表4 基于特征波段和衍生高光谱指数的盐分反演模型构建Tab. 4 Construction of salt inversion model based on sensitive band and derived hyperspectral index

2.3 土壤盐分与Sentinel-2B 影像反射率相关性分析与建模在ENVI 软件中利用“Band Math”计算表3 所示植被指数和盐分指数，并在ArcGIS 软件中运用“Extract Multi values to points”工具提取各样点对应多光谱指数值. 将提取的多光谱反射率和多光谱指数进行数学变换，并与含盐量进行相关性分析，筛选多光谱最优变换、特征波段和特征光谱指数.

表5 表明，土壤盐分与影像各波段反射率的相关性在0.01 置信水平上通过检验. 其中，光谱倒数变换与含盐量呈正相关关系，其余3 种变换与盐分呈现负相关关系. 多光谱反射率经1/R变换后与盐分相关性较好，B6波段与盐分的相关系数达到最大值0.502. 经R、R0.5、lg（R）变换后与盐分的最大相关系数分别为0.473、0.475、0.488，对应波段分别为B11、B11、B8. lg（R）变换特征波段最多，有6个波段相关系数绝对值大于0.45.

表5 含盐量与影像各波段反射率的相关性Tab. 5 Correlation between soil salinity and reflectance of multispectral bands

表6 表明，光谱指数SI1、S3、Int1、Int2和冠层响 应 植 被 指 数（Canopy Response Salinity Index，CRSI）与土壤盐分相关性较好. 除1/R变换外，其余变换的SI1、S3和Int2均在0.01 置信水平显著. 其中S3指数对盐分的响应程度最好，除1/R变换外相关系数绝对值均大于0.5，经lg（R）变换后与含盐量相关性最好（0.535）. 其余光谱指数与盐分相关性较低，相关系数未通过置信度检验.

表6 含盐量与多光谱指数的相关性Tab. 6 Correlation between soil salinity and multispectral indices

基于上述分析，选取1/R变换后的B6、B7、B8、B8a波段作为特征波段，lg（R）变换后的S3、Int1、Int2指数作为特征光谱指数，构建多光谱RR 和PLSR 盐分反演模型. 如表7 所示，多光谱RR 和PLSR 模型显著性相差不大，但特征光谱指数模型整体略高于特征波段模型，其中基于特征光谱指数的PLSR 模型RMSE 最小，拟合度最高. 由于采样时间为2019 年5 月中旬至下旬，此时植被长势良好，反映在多光谱影像上多为混合像元，采样点提取的土壤光谱反射率不可避免地掺杂了植被光谱特征. 为减小植被等的影响，本文基于特征光谱指数进行高-多光谱数据融合，以精度较高的高光谱反演模型修正多光谱模型，以提升模型反演精度.

表7 多光谱RR 和PLSR 反演模型Tab. 7 Multispectral RR and PLSR inversion model

2.4 高-多光谱数据联合及模型验证Sentinel-2B影像的光谱混合像元难以准确探测地物细微的光谱特征，所构建模型难以实现空间范围内的准确反演. 实测高光谱反射率光谱特征连续，可视为遥感影像的纯净像元，能够实现地表参量的精细探测.为提升土壤盐分反演精度，本文对3 种特征光谱指数进行曲线估计，选取最优拟合度模型进行高光谱实测端元和多光谱像元间的光谱匹配，实现地面点到空间面域尺度的土壤盐渍化反演.

由图4 可知，光谱指数S3和Int2以三次多项式曲线模型为最优拟合度模型，Int1指数的最优拟合度模型为幂模型. 3 种高-多光谱指数间均为正相关关系，其中Int2指数（R2=0.740，RMSE=0.023 g·kg-1）拟合度最高，光谱数据融合效果最好.

图4 银川平原高-多光谱敏感光谱指数融合Fig. 4 Fusion of hyperspectral-multispectral sensitive indices in Yinchuan Plain

利用光谱匹配后的多光谱指数构建RR 和PLSR 盐分反演回归模型（表8），拟合度比表9 中单独影像建模（R2=0.412，RMSE=6.941）有了较大提升，R2最大值提升了0.309，RMSE 减小了2.085 g·kg-1.其中，PLSR 模型较RR 模型拟合度更好，其R2更大，RMSE 更小，表明基于光谱匹配后建立的PLSR模型为最优盐分反演模型.

表8 光谱匹配后RR 和PLSR 盐分反演模型Tab. 8 Salt inversion model of RR and PLSR after spectral matching

利用验证集样本对拟合度最高的高光谱反演模型、多光谱反演模型和光谱匹配后反演模型进行验证，结果如图5 所示. 各反演模型均能通过验证，其中高光谱反演模型验证精度最高，多光谱最低，经光谱匹配后的反演模型验证精度较多光谱模型也有较大提升. 可见，基于特征光谱指数进行地表实测高光谱与影像间光谱匹配，进而提升单独利用影像构建模型的精度是可行的. 因此，本文基于光谱匹配后反演模型预测银川平原空间尺度上的盐渍化分布格局.

图5 银川平原土壤盐分的不同光谱数据模型验证（单位：g·kg-1）Fig. 5 Validation of different spectral data models of soil salinity in Yinchuan Plain（unit：g·kg-1）

2.5 土壤盐分反演及结果统计利用光谱匹配后

PLSR 反演模型对银川平原土壤盐分进行预测，并按盐分含量进行分级和统计，结果如图6 和表9.银川平原盐分空间格局总体呈现南轻北重分布趋势. 原因可能是平原北部地势低洼、地下埋水较深且蒸发强烈、引黄河水和沟渠灌溉且排水不畅，盐分极易于土壤表层聚积形成盐结皮、盐壳、和龟裂碱土等现象，如平原北部平罗县、大武口区，以及黄河沿岸主要分布盐土，且面积较大，占盐渍化总面积的59.38%. 而平原南部地势高，且作为黄河入水口排水顺畅，携带盐基离子外流，土壤积盐较少，如平原南部永宁县、青铜峡市和利通区主要分布非盐渍化和轻度盐渍化土，占区域总面积的29.87%. 反演结果与实际采样情况较为一致，表明光谱匹配后的反演模型可用于该区土壤盐分预测.

表9 银川平原不同等级盐分反演像元数及占比Tab. 9 Number and proportion of inversion pixels of different grades of salinity in Yinchuan Plain

图6 银川平原土壤盐分反演Fig. 6 Inversion result of soil salinity in Yinchuan Plain

3 讨论与结论

3.1 讨论融合多源卫星影像数据进行优势互补，从而反演土壤盐分已成为盐分预测和空间可视化表达的重要手段. 如利用WorldView-2[39]、IKONOS[40]等高分辨率的多光谱影像与Hyperion 高光谱影像数据融合，与OLI[4]中分辨率多光谱影像和HIS 高光谱影像的数据融合能有效地提升盐分预测精度.ASD 地表实测高光谱数据的光谱连续，像元相对纯净，能够实现地表参量的精细探测，已被广泛用于卫星影像融合和模型精度提升. 如贾萍萍[41-42]和尚天浩等[43]利用重采样法实现了ASD 高光谱和OLI 影像的数据融合，构建了宁夏银北地区土壤盐分的高精度预测模型；王爽[1]和姚远等[14]借助统计回归法进行了ASD 数据与TM 影像数据的尺度转换. 同时，数理统计回归法也被认为是衔接和转换地-星光谱关系的一种有效途径，因此本文利用数理统计回归法进行地-星光谱匹配具有可行性和普适性.

利用统计回归方法进行光谱数据匹配，对于实现由点及面的方法实践，解决单纯利用影像建模精度较差等问题具有重要意义. 研究表明，利用敏感波段构建的特征光谱参量[44]和进行光谱数学变换[10]能有效增强光谱响应特征，提升模型精度. 本文以3 种光谱响应指数S3、Int1、Int2为匹配参量，通过回归方法实现地-星遥感数据光谱匹配，发现数据匹配后的模型能有效地提升模型精度，相比单独利用影像建模，其R2提升了0.309，RMSE 减小了2.085 g·kg-1，且反演结果与野外实地采样具有一致性，均呈现出盐分北高南低的空间格局特征，这与孙亚楠等[34]的研究结果一致.

由于土壤盐分受到地下水、地势、降水量、灌溉渠系分布等空间因子的影响，本文利用光谱反射率和土壤盐分较好的光谱响应关系，构建统计回归模型对盐分空间分布格局进行了定量估测，未考虑这些因素对盐渍化的影响；且Sentinel-2B 多光谱影像为中等空间分辨率数据，所构建模型的普适性有待进一步验证. 下一步将机器学习和深度学习等方法用于土壤盐渍化反演，期望进一步提高模型的反演精度和普适性.

3.2 结论本文基于实测高光谱和Sentinel-2B 影像数据，通过光谱变换，选取特征波段和特征光谱指数，构建RR 和PLSR 盐分反演模型，并以特征光谱指数为光谱参量进行地-星数据光谱匹配，利用数据匹配后的模型对银川平原土壤盐分进行了定量反演. 主要结论如下：

（1）光谱变换能显著地提高盐分与光谱反射率敏感程度. 高光谱最优变换为[lg（R）]′，特征波段为ρ435、ρ500、ρ499、ρ498、和ρ501. 1/R变换后B6、B7、B8和B8a为多光谱特征波段，lg（R）变换后S3、Int1、Int2指数为多光谱特征光谱指数.

（2）以特征光谱指数为敏感参量，对于实现高光谱实测端元到多光谱影像像元尺度间的过度，提升高-多光谱融合反演模型精度是行之有效的.

（3）经光谱数据匹配后构建的PLSR 模型为最优盐分反演模型，其精度比单独影像建模有了较大提升，R2提升了0.309，RMSE 减小了2.085 g·kg-1.

（4）利用光谱数据匹配后的模型对银川平原空间范围内土壤盐分进行预测，反演结果与野外实地调研具有较好的一致性. 银川平原盐分空间格局表现为南轻北重，盐渍化土壤类型多为重度盐渍化土和盐土，南部主要分布非盐渍化和轻度盐渍化土.