CN114076738B - 一种利用遥感构建指数估算秸秆残留农田土壤盐分的方法 - Google Patents

Abstract

本发明公开了一种利用遥感构建指数估算秸秆残留农田土壤盐分的方法。本发明联合卫星数据和地面调查数据，提供了一种利用Sentinel‑2卫星数据和Sentinel‑1雷达数据构建特征指数用于估算秸秆残留农田土壤盐分的方法。通过本发明的方法对干旱区秸秆残留农田土壤盐分含量进行估算后，可得到大空间尺度的干旱区秸秆残留农田土壤盐分含量分布图，弥补了在秸秆残留农田土壤盐分估算中缺乏特征指数的短板，为干旱区秸秆残留农田土壤盐分含量的估算提供了新方法，有利于大区域尺度农田土壤盐分的治理改良政策的制定，具有一定的理论、实践意义和推广应用价值。

Description

一种利用遥感构建指数估算秸秆残留农田土壤盐分的方法

技术领域

本发明属于遥感反演领域，具体涉及一种利用遥感构建特征指数用于估算秸秆残留农田土壤盐分的方法。

背景技术

土壤盐渍化包括原生盐渍化和次生盐渍化，尤其是农田土壤盐渍化严重危害了土壤健康和作物生产。干旱和半干旱地区土壤蒸散量远高于降水量，土壤表面盐分积累结成盐壳，然而，不合理的农田管理(如漫灌洗盐)加剧了二次盐渍化进程。土壤盐渍化的定量估算对于耕地土壤资源的保护和利用具有重要意义。遥感技术因其空间范围覆盖广、观测空间分辨率高、时间重返周期性短的特点，被广泛用于土壤盐分的检测与评估。

在利用遥感手段定量反演土壤盐分中，针对不同的土壤盐渍化类型，以往的研究构建并使用了多个盐分指数、植被指数、遥感波段等用于建模回归。然而，因为地理环境的不同，气候类型、地貌特征、盐渍化类型、土壤类型、植被类型、社会经济要素等本地化特征明显，在应用到不同地区时，不同指数的通用性表现出一定的差异。特别是对于存在秸秆残留的农田土壤而言，此类农田土壤由于表面存在秸秆覆盖等不利因素，因此其在遥感反演盐分含量时缺乏特征指数的问题。因此，针对干旱区秸秆残留农田土壤盐渍化的监测，需要构建并优选一系列指数用于指示农田土壤表层的盐分，从而实现对土壤盐分的高精度定量反演与可视化，为干旱区秸秆残留农田土壤盐分含量的估算提供了新方法

发明内容

本发明的目的在于解决现有技术中存在的问题，联合地面调查数据，提供一种利用Sentinel-2卫星数据和Sentinel-1雷达数据构建特征指数用于估算秸秆残留农田土壤盐分的方法。

为了实现上述发明目的，本发明采用的具体技术方案如下：

第一方面，本发明提供了一种利用遥感构建指数估算秸秆残留农田土壤盐分的方法，其步骤如下：

S1、获取待测区域在待估算时间下对应的Sentinel-2卫星数据、Sentinel-1雷达数据、土壤采样数据集；所述待测区域的土地类型为干旱区秸秆残留农田；所述土壤采样数据集中包含待测区域内不同采样点对应的土壤表层盐分含量；

S2、对Sentinel-2卫星数据中的多光谱影像进行预处理，得到不同波段的表面反射率数据，再通过多种光谱变换得到光谱数据集；所述光谱变换形式包括对数变换、倒数变换和微分变换；所述光谱数据集中包含影像中每一个栅格的光谱集合，所述光谱集合包括原始光谱和经过所述光谱变换后得到的多种变换光谱；对Sentinel-1雷达数据进行预处理，得到包含垂直发射-垂直接收(VV)和垂直发射-水平接收(VH)的雷达后向散射系数的雷达数据集；

S3、遍历土壤表层盐分的敏感波段组合中不同的波段，并将其组合成不同的二维、三维和四维指数，得到多维指数集合；再针对所述土壤采样数据集中的每个采样点，从所述光谱数据集中获取每个采样点对应的光谱集合，基于光谱集合中的每一条光谱分别计算所述多维指数集合中每一种指数的值；然后以每一条光谱下的每一种指数作为待筛选指数，将所述土壤采样数据集中各土壤采样点的土壤表层盐分含量分别与不同的待筛选指数进行相关性分析，筛选出相关性最高的至少一个待筛选指数作为高敏感指数；

S4、以所述土壤采样数据集中的各采样点为回归样本，以所述高敏感指数作为自变量，以土壤表层盐分含量作为因变量，建立回归模型；并从所述光谱数据集和所述雷达数据集中获取待测区域内每一个栅格对应的所述高敏感指数的值，根据所述回归模型计算得到每一个栅格处的土壤表层盐分含量，最终形成待测区域内秸秆残留农田土壤盐分含量的空间分布图。

作为优选，所述Sentinel-2卫星数据中的多光谱影像的预处理包括辐射定标和大气校正，且其中大气校正所选用的模型为FLAASH。

作为优选，所述Sentinel-1雷达数据的预处理包括轨道校正、热噪声去除、辐射校正、Lee滤波和多普勒地形校正，并选用水云模型去除植被含水量对土壤后向散射系数的影响后，将垂直发射-垂直接收(VV)和垂直发射-水平接收(VH)的雷达后向散射系数转换并得到分贝(dB)格式的雷达数据集。

作为优选，所述土壤采样数据集中，每一个采样点对应的土壤表层盐分含量由土壤表面0～0.2m的土层通过电导率法测定得到。

作为优选，所述Sentinel-2卫星数据和Sentinel-1雷达数据的空间分辨率为10米。

作为优选，所述敏感波段组合由光学遥感敏感波段和雷达敏感波段组成，其中所述光学遥感敏感波段为Sentinel-2卫星数据中的蓝(Blue)波段、绿(Green)波段、红(Red)波段、近红外(NIR)波段、短波红外1波段(SWIR 1)、短波红外2波段(SWIR 2)和3个植被红边指数(Red Edge 1，Red Edge 2，Red Edge3)波段，所述雷达敏感波段为Sentinel-1雷达数据中垂直发射-垂直接收(VV)和垂直发射-水平接收(VH)的雷达后向散射系数，。

作为优选，所述微分变换包括一阶微分变换、二阶微分变换和三阶微分变换，影像中每一个栅格的光谱集合均包含6条光谱，分别为原始光谱、对数变换光谱、倒数变换光谱、一阶微分变换光谱、二阶微分变换光谱和三阶微分变换光谱。

作为优选，所述多维指数集合中，二维、三维和四维指数均由相应维数的波段通过组合计算得到，组合计算的形式为差值计算、比值计算和幂计算中的一种或多种结合。

作为优选，所述回归模型为多项式回归模型，其形式为：

EC=a×ERSSI²+b×ERSSI+C

式中：EC为土壤表层盐分含量，a、b、c分别为回归系数。

需说明的是，上述第一方面中所提供的各技术方案，可以应用于任意的待测区域和时间，只要多维指数集合中所涵盖的指数可选范围足够广泛，一般都可以筛选得到高敏感指数并最终构建得到回归模型。

在本发明中，通过后续实施例，筛选到了一种高敏感指数，该指数是能有效指示干旱区秸秆残留农田土壤盐分的新指数增强型秸秆残留盐分指数(Enhanced Residues SoilSalinity Index，ERSSI)，该高敏感指数为基于二阶微分变换光谱计算的三维指数，三维指数形式为

其中ERSSI为高敏感指数，Green为Sentinel-2卫星数据的绿(Green)波段，Blue为Sentinel-2卫星数据的蓝(Blue)波段，SWIR₁为Sentinel-2卫星数据的短波红外1波段(SWIR 1)。

基于这种高敏感指数，本发明进一步对第一方面的方法进行了简化，提出了下述第二方面的技术方案。

第二方面，本发明提供了一种利用遥感构建指数估算秸秆残留农田土壤盐分的方法，其步骤如下：

1)、获取待测区域在待估算时间下对应的Sentinel-2卫星数据；所述待测区域的土地类型为干旱区秸秆残留农田；

2)、对Sentinel-2卫星数据中的多光谱影像进行预处理，得到不同波段的表面反射率数据，再对每个栅格的光谱进行二阶微分变换得到二阶微分变换光谱；

3)、基于待测区域内每个栅格的二阶微分变换光谱计算对应的三维指数，三维指数形式为

其中ERSSI为高敏感指数，Green为Sentinel-2卫星数据的绿(Green)波段，Blue为Sentinel-2卫星数据的蓝(Blue)波段，SWIR₁为Sentinel-2卫星数据的短波红外1波段(SWIR 1)；

4)、将待测区域内每个栅格的三维指数作为自变量，基于预先针对待测区域构建的回归模型估计每个栅格对应的土壤表层盐分含量，最终形成待测区域内秸秆残留农田土壤盐分含量的空间分布图。

本发明相对于现有技术而言，具有以下有益效果：

本发明联合卫星数据和地面调查数据，提供了一种利用Sentinel-2卫星数据和Sentinel-1雷达数据构建特征指数用于估算秸秆残留农田土壤盐分的方法。通过本发明的方法对干旱区秸秆残留农田土壤盐分含量进行估算后，弥补了在秸秆残留农田土壤盐分估算中缺乏特征指数的短板，为干旱区秸秆残留农田土壤盐分含量的估算提供了新方法，有利于大区域尺度农田土壤盐分的治理改良政策的制定，具有一定的理论、实践意义和推广应用价值。

另外，本发明筛选得到了一种有效指示干旱区秸秆残留农田土壤盐分的新指数增强型秸秆残留盐分指数(Enhanced Residues Soil Salinity Index，ERSSI)，以Sentinel-2卫星遥感数据和新指数ERSSI构建了多项式线性回归反演干旱区秸秆残留农田的土壤盐分含量的方法，最终可以得到空间分辨率为10米的高空间分辨率、高质量的土壤盐分含量空间变异结果。

附图说明

图1为使用增强型秸秆残留盐分指数(Enhanced Residues Soil SalinityIndex，ERSSI)和多项式回归的土壤盐分含量估算值在Q1(建模集)和Q2(独立验证集)中表现的散点图；

图2为本实施方式估算得到的2020年南疆干旱区秸秆残留农田土壤的盐分含量分布图。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步阐述和说明。本发明中各个实施方式的技术特征在没有相互冲突的前提下，均可进行相应组合。

在本发明的一个较佳实施例中，提供了一种利用遥感构建指数估算秸秆残留农田土壤盐分的方法，该方法用于对土地类型为干旱区秸秆残留农田的待测区域进行土壤盐分的估算。该方法的具体步骤如下：

S1、获取待测区域在待估算时间下对应的Sentinel-2卫星数据、Sentinel-1雷达数据、土壤采样数据集。

其中Sentinel-2卫星数据是一种光学遥感影像数据，Sentinel-1雷达数据是一种雷达数据，土壤采样数据集是一种地面调查数据，三者来源不同，属于多源数据组合。为了保证估算的准确性，Sentinel-2卫星数据、Sentinel-1雷达数据、土壤采样数据集这三类数据均需要是在待估算时间下同步采集的。例如，待估算时间是某一个时期，那么Sentinel-2卫星数据、Sentinel-1雷达数据、土壤采样数据集均需要是这一个时期同步采集的数据。但由于土壤盐分含量以及冠层信息在短时间内的相对稳定性，因此本发明中不同数据源的待估算时间允许存在一定的时间偏差。

Sentinel-2卫星数据、Sentinel-1雷达数据均可通过相应的卫星数据发布方处下载。Sentinel-2卫星数据来源于哨兵2号卫星，该卫星携带一枚多光谱成像仪(MSI)，高度为786km，可覆盖13个光谱波段。Sentinel-1雷达数据来源于哨兵1号卫星，由两颗卫星组成，载有C波段合成孔径雷达，可提供连续图像(白天、夜晚和各种天气)。本实施例中采用的Sentinel-2卫星数据和Sentinel-1雷达数据的空间分辨率为10米。而土壤采样数据集是通过对待测区域进行地面调查获得的，地面调查应当在待测区域内进行采样点的布设，并在每个采样点处采集相应的土壤表层样品(0～0.2m的土层)，然后对土壤表层样品测定其土壤盐分含量，由此土壤采样数据集中最终将包含待测区域内不同采样点对应的土壤表层盐分含量。

在本实施例中，土壤采样数据集中每一个采样点对应的土壤表层盐分含量由土壤表层样品通过电导率法测定得到。电导率法的测定过程如下：将土壤样品风干、研磨并过2毫米筛，以土水比1：5得到土壤浸出液，过滤取上清液测定其电导率EC_1:5，得到土壤样品的电导率数据，作为土壤表层盐分含量的表征值。

需说明的是，由于土壤样品的电导率与土壤表层盐分含量存在正相关关系，因此本实施例中可以直接用电导率来表征待测区域的土壤盐渍化水平。但是在其他实施例中，也可以根据两者的相关关系进行换算后得到土壤表层盐分含量的绝对值，土壤表层可溶性盐含量c的换公式为c(g/kg)＝0.0275*EC_1:5-0.0573，其中EC_1:5是指土壤的电导率，单位为mS/m。因此，土壤采样数据集中采样点对应的土壤表层盐分含量既可以用电导率EC_1:5代表相对含量，也可以用土壤表层可溶性盐含量c来代表绝对含量，具体可根据实际需要进行选择。

S2：对多源数据的融合和处理，将S1获取的Sentinel-2卫星数据、Sentinel-1雷达数据分别进行预处理和计算，具体依次进行如下步骤：

S21：对Sentinel-2卫星数据中的多光谱影像(MSI图像)进行预处理计算，得到不同波段的表面反射率数据，其中预处理包括辐射定标和大气校正，且其中大气校正所选用的模型为FLAASH。另外，也对Sentinel-1雷达数据进行预处理，预处理包括轨道校正、热噪声去除、辐射校正、Lee滤波和多普勒地形校正，应用水云模型去除植被含水量对土壤后向散射系数的影响，将垂直发射-垂直接收(VV)和垂直发射-水平接收(VH)的雷达后向散射系数转换并得到分贝(dB)格式的雷达数据集，得到地表反射特征数据。

S22：对经过预处理的Sentinel-2卫星数据的各个波段进行光谱变换，得到光谱变换后的光谱数据集。本步骤中的光谱变换采用逐波段计算方法，而光谱变换形式包括对数变换、倒数变换和微分变换三类，微分变换又包含一阶微分变换、二阶微分变换和三阶微分变换。由此，对Sentinel-2影像逐个波段地进行对数变换、倒数变换、一阶微分变换、二阶微分变换和三阶微分变换，得到原始光谱曲线以及五种数学计算变换后的光谱曲线。因此，最终影像中每一个栅格的光谱集合均包含6条光谱，分别为原始光谱、对数变换光谱、倒数变换光谱、一阶微分变换光谱、二阶微分变换光谱和三阶微分变换光谱。

S3、遍历土壤表层盐分的敏感波段组合中不同的波段，并将其组合成不同的二维、三维和四维指数，得到多维指数集合。

其中土壤表层盐分的敏感波段组合是指与土壤盐分相关的所有波段组成的集合，其中哪些波段与土壤盐分相关可以根据相关性分析确定，亦可通过相关的参考文献和历史研究进行确定。由于本发明中遥感数据包含了Sentinel-2卫星数据和Sentinel-1雷达数据，因此土壤表层盐分的敏感波段组合也需要由光学遥感敏感波段和雷达敏感波段组成。本发明在历史研究报道的基础上，采用归纳统计方法进行敏感波段筛选，根据影响土壤盐分监测精度的土壤、作物、水分因子，以及盐渍土的光谱反射特性，总结归纳现有的盐分指数、植被指数和土壤相关指数等使用的波段，在光学遥感数据中筛选得到蓝、绿、红、近红外、短波红外、植被红边指数为敏感波段，一共得到9个光学敏感波段；在雷达数据中筛选得到雷达后向散射系数VV、VH为敏感波段，一共得到2个雷达敏感波段。本实施例中，共得到11个对土壤盐分敏感的波段构成敏感波段组合，其中光学遥感敏感波段为Sentinel-2卫星数据中的蓝波段(Blue，波长458-523nm)、绿波段(Green，波长543-578nm)、红波段(Red，波长650-680nm)、近红外波段(NIR，波长785-900nm)、短波红外1波段(SWIR 1，波长1565-1655nm)、短波红外2波段(SWIR 2，波长2100-2280nm)和3个植被红边指数波段(Red Edge，波长分别为698-713nm、733-748nm、773-793nm)，而雷达敏感波段为Sentinel-1雷达数据中垂直发射-垂直接收(VV)和垂直发射-水平接收(VH)的雷达后向散射系数。

当得到上述包含11个敏感波段的敏感波段组合后，即可将这些敏感波段组合成多维的指数。N维指数是指用N个敏感波段构建的指数，指数的形式不限，可以采用不同类型的计算式。在本实施例的多维指数集合中，指数的维数包括二维、三维和四维三种，分别由N＝3个、N＝4个、N＝5个敏感波段通过计算式组合计算得到。组合计算的计算式形式可以为差值计算、比值计算和幂计算中的一种，差值计算形式为(A-B)，比值计算形式为(A/B)，幂计算形式为(A^k)，A和B各自为一个敏感波段，k为常数。当然，组合计算的计算式形式也可以是差值计算、比值计算和幂计算中两种或三种的结合，即将差值计算、比值计算和幂计算进行结合构建为更复杂的计算式，如(A-B)/C，[(A-B)/(A+B)]^k，[(A-B)/(A+C)]^k等等，C也为一个敏感波段。对于任意一种形式的N维指数，所用的N个敏感波段可从11个敏感波段中进行遍历采样，因此本发明的多维指数集合中指数存在指数形式和敏感波段两个区分维度，相同形式的指数因为所选择的敏感波段不同还存在众多变种，每一个N维指数通过选择不同的N个敏感波段可以进一步衍生出一系列的指数，这些指数形式不同或者敏感波段不同的指数都被纳入多维指数集合中并分配唯一的ID，作为用于后续筛选的指数。

当得到多维指数集合后，再进行指数的筛选，从中找到与秸秆残留农田土壤盐分含量最相关的指数。由于多维指数集合中的每一个指数都是基于一个或多个波段进行计算的，而波段的波段值需要从光谱曲线上获取。在前述S22步骤中得到了6条光谱组成的光谱数据集，因此，在进行最敏感指数筛选时，上述多维指数集合中的每一个指数还需要分别从6条光谱上进行采样。具体进行筛选时，可针对前述S1的土壤采样数据集中的每个采样点，从前述S22的光谱数据集中获取每个采样点对应的光谱集合(一共包含6条光谱曲线)，基于光谱集合中的每一条光谱分别计算前述多维指数集合中每一种指数的值；然后以每一条光谱下的每一种指数作为待筛选指数，将上述土壤采样数据集中各土壤采样点的土壤表层盐分含量分别与不同的待筛选指数进行相关性分析。不同的待筛选指数之间可能是采样的光谱不同，也可能是指数的ID不同。每一次相关性分析时，以土壤采样数据集中各土壤采样点的土壤表层盐分含量作为第一条数据序列，以基于土壤采样数据集中各土壤采样点的某一条光谱计算得到的某一种指数作为第二条数据序列，计算两条数据序列的相关性指数。最终，筛选出相关性最高的一个待筛选指数作为高敏感指数。

S4、当确定高敏感指数后，即可以土壤采样数据集中的各采样点为回归样本，以采样点处的高敏感指数作为自变量，以采样点处的土壤表层盐分含量作为因变量，通过回归分析建立回归模型。本实施例中，所采用的回归模型为多项式回归模型，其形式为：

EC=a×ERSSI²+b×ERSSI+C

式中：EC为土壤表层盐分含量，a、b、c分别为回归系数。

当构建得到回归模型后，从前述S2得到的光谱数据集和雷达数据集中获取待测区域内每一个栅格对应的高敏感指数的值，将其输入构建得到的回归模型中进而计算得到每一个栅格处的土壤表层盐分含量，所有栅格计算完毕后，最终形成待测区域内秸秆残留农田土壤盐分含量的空间分布图。

在本发明的后续实施例中，筛选到了一个高敏感指数，即为基于二阶微分变换光谱计算的三维指数，三维指数形式为

其中ERSSI为高敏感指数，Green为Sentinel-2卫星数据的绿(Green)波段，Blue为Sentinel-2卫星数据的蓝(Blue)波段，SWIR₁为Sentinel-2卫星数据的短波红外1波段(SWIR1)，Green、Blue和SWIR₁波段值都需要从二阶微分变换光谱上进行采样。当然，如果构建的多维指数集合不同，亦可能可以筛选到其他的高敏感指数，对此不做限定。

需说明的是，前述的三维指数ERSSI可以用于直接预测秸秆残留农田土壤盐分，对于干旱区秸秆残留农而言，仅依靠该指数即可实现农田土壤盐分的反演，无需依赖于其他的波段。因此，本发明中进一步提供了一种利用遥感构建指数估算秸秆残留农田土壤盐分的方法，其步骤如下：

1)、获取待测区域在待估算时间下对应的Sentinel-2卫星数据；其中待测区域的土地类型为干旱区秸秆残留农田，即位于干旱区且地表残留有秸秆的农田。

2)、对Sentinel-2卫星数据中的多光谱影像进行预处理，得到不同波段的表面反射率数据，再对每个栅格的光谱进行二阶微分变换得到二阶微分变换光谱。多光谱影像的预处理可参见前述方法。

3)、基于待测区域内每个栅格的二阶微分变换光谱计算对应的三维指数ERSSI，计算公式如前所述。

4)、将待测区域内每个栅格的三维指数作为自变量，基于预先针对待测区域构建的回归模型估计每个栅格对应的土壤表层盐分含量，最终形成待测区域内秸秆残留农田土壤盐分含量的空间分布图。

为了进一步便于理解本发明的优点，下面将上述实施例中S1～S4步骤的利用遥感构建指数估算秸秆残留农田土壤盐分的方法应用于一个具体的案例中，以便于展示具体的技术效果。

实施例

选取新疆维吾尔自治区的棉花种植田块(81°17'-81°22'E，40°28'-40°31'N)作为待测区域，利用2020年11月1日地面调查获取的土壤样点数据得到的0.2米土壤盐分含量作为因变量，以2020年11月1日的单景Sentinel-2卫星遥感影像数据及2020年11月2日的Sentinel-1数据作为自变量，基于光谱变换和二维、三维、四维指数的构建，得到与干旱区秸秆残留农田土壤盐分的高敏感型指数增强型秸秆残留盐分指数(Enhanced ResiduesSoil Salinity Index，ERSSI)。并以ERSSI为自变量，以土壤盐分含量为因变量，利用线性多项式回归模型估算，得到待测区域的秸秆残留农田土壤盐分含量分布图。该指数及估算方法具体如下：

步骤1)数据获取：以干旱区秸秆残留农田为估算对象，根据待估算的时间，获取待测区域同一时期的Sentinel-2卫星数据、Sentinel-1雷达数据、土壤表层盐分含量数据集；Sentinel-2卫星数据和Sentinel-1雷达数据的空间分辨率为10米。

其中，Sentinel-2数据是欧洲航天局(ESA)发布的L1级产品。Sentinel-2由两颗极轨卫星(Sentinel-2A和Sentinel-2B)组成，两颗卫星数据互补，重访周期为5天，轨道周期为100分钟，轨道高度为786公里。扫描宽度为290公里，轨道倾角为98.62°。多光谱扫描成像数据(MSI)在可见光至红色边缘区域提供一个空间分辨率为60m的波段、三个10m波段和两个30m波段，该数据可在美国地质调查局(USGS)下载免费下载，该步骤使用了2020年11月1日的Sentinel-2B数据；Sentinel-1卫星是欧洲航天局哥白尼计划(GMES)发射的一颗地球观测卫星，由两颗卫星(Sentinel-1A和Sentinel-1B)组成，携带C波段合成孔径雷达，允许以6天的时间分辨率提供全天图像。该步骤使用了2020年11月2日的Sentinel-1A数据，因为该时期的影像与野外实验日期接近，土壤和冠层条件几乎相同。Sentinel-1A数据上升方向和干涉宽(IW)模式下获得，具有两个双偏振，空间分辨率为10m。

土壤表层盐分含量数据集是通过对不同点位的0～0.2m土层进行土壤样品采样，并运回实验室进行分析得到的额。获取的土壤样品采用电导率法测定其电导率：将土壤样品风干、研磨并过2毫米筛，以土水比1：5得到土壤浸出液，过滤取上清液测定其电导率，得到土壤样品的电导率数据，用于表征待测区域的土壤盐渍化水平。采用地理分区法将土壤样点数据分为两个区域，根据样点分布的地理位置不同分为两个子区域Q1、Q2，其中，Q1的电导率数据集作为建模集，用于新指数的构建、相关性评估和多项式回归方程的构建，Q2的电导率数据集作为独立验证集，用于独立验证指数和多项式回归方程的有效性。

步骤2)数据预处理：基于多源数据的融合和处理，将步骤1)获取的光学遥感影像数据、雷达数据和土壤样品数据进行预处理和计算。

将步骤1)获取的Sentinel-2遥感数据进行预处理：利用Sentinel ApplicationPlatform(SNAP)包中的Sen2Cor模块对Sentinel-2数据进行辐射校准和大气校正，将MSI图像转换为表面反射率格式输出，得到分辨率为10米的Sentinel-2影像。另外，将步骤1)获取的Sentinel-1雷达数据进行预处理：依次进行轨道校正、热噪声去除、辐射校正、Lee滤波和多普勒地形校正，应用水云模型去除植被含水量对土壤后向散射系数的影响，将垂直发射-垂直接收(VV)和垂直发射-水平接收(VH)的雷达后向散射系数转换并得到分贝(dB)格式的雷达数据集，得到表面反射率数据。

将经过预处理的Sentinel-2卫星数据的各个波段采用逐波段计算方法进行光谱变换，逐个波段进行对数变换、倒数变换、一阶微分变换、二阶微分变换和三阶微分变换一共6中数学计算光谱变换，得到原始光谱曲线以及五种数学计算变换后的光谱曲线，每一个栅格的6条光谱作为该栅格对应的光谱集合，所有栅格的光谱集合作为光谱数据集。

步骤3)新指数的构建：遍历土壤表层盐分的敏感波段组合中不同的波段，并将其组合成不同的二维、三维和四维指数，得到多维指数集合。

在本实施例中，的对土壤盐分敏感的敏感波段组合是采用归纳统计方法筛选的，根据影响土壤盐分监测精度的土壤、作物、水分因子，以及盐渍土的光谱反射特性，总结归纳现有的盐分指数、植被指数和土壤相关指数等使用的波段，在光学遥感数据中筛选得到蓝、绿、红、近红外、短波红外、植被红边指数为敏感波段，在Sentinel-2原始光谱及五种数学计算变换后的光谱中提取以上波段，得到9个光学敏感波段；在雷达数据中筛选得到土壤的雷达后向散射系数VV、VH为敏感波段，在Sentinel-1数据中计算得到2个雷达敏感波段。共得到11个对土壤盐分敏感的波段。采用归纳、增减维度方法进行新指数的构建。首先，对现有盐分指数的计算形式进行总结，在指数的构建形式上归纳得到差值、比值、幂计算三种计算形式，在此基础上，对三种计算形式中的变量改变维度、减少和增加维度，进行维度的改变，用于构建新指数，以11个敏感波段为变量，构建27个指数形式，包括9个二维指数、8个三维指数、10个四维指数，指数及其计算式如表1所示；

表1二维、三维、四维指数及其计算式

其中，B_i，B_j，B_k，B_h分别为蓝、绿、红、近红外、短波红外、植被红边指数波段以及雷达后向散射系数VV、VH这11个波段中的任意四种。每一个指数形式均需要遍历整个敏感波段集合，选择不同的波段构成不同的指数加入多维指数集合中，用于后续的筛选。

最后，再针对Q1中的每个采样点，从光谱数据集中获取每个采样点对应的光谱集合，基于光谱集合中的每一条光谱(分别为原始光谱、对数变换光谱、倒数变换光谱、一阶微分变换光谱、二阶微分变换光谱和三阶微分变换光谱)分别计算多维指数集合中每一种指数的值。然后以每一条光谱下的每一种指数作为待筛选指数，将土壤采样数据集中各土壤采样点的土壤表层盐分含量分别与不同的待筛选指数进行相关性分析，筛选出相关性最高的至少一个待筛选指数作为高敏感指数。本实施例中，采用相关性分析评估不同光谱下不同指数(共六组，每组27种指数形式)与Q1中土壤电导率数据集的相关性(如表2所示)，同时在Q2中土壤电导率数据集进行独立验证(如表3所示)，得到相关性最高的指数。

表2由六种形式的光谱构建的二维、三维、四维指数在Q1中的最高相关性(r)

表3由六种形式的光谱构建的二维、三维、四维指数在Q2中的最高相关性(r)

在Q1和Q2的相关性分析中，得到对干旱区秸秆覆盖农田土壤盐分相关性最高的指数为二阶导数变换后构建的三维指数，命名为增强型秸秆残留盐分指数(EnhancedResidues Soil Salinity Index，ERSSI)，该新指数的计算公式为：

其中ERSSI为高敏感指数，Green为Sentinel-2卫星数据的绿(Green)波段，Blue为Sentinel-2卫星数据的蓝(Blue)波段，SWIR₁为Sentinel-2卫星数据的短波红外1波段(SWIR-1)，Green、Blue和SWIR₁波段值都需要从二阶微分变换光谱上进行采样。

图1为使用ERSSI和多项式回归的土壤盐分含量估算值在Q1(建模集)和Q2(独立验证集)中表现的散点图，表面本发明估算得到的盐分数据相对于地面实测值的拟合程度较高。

ERSSI以比率的形式组合了绿、蓝和短波红外波段。绿波段对植被类型敏感，而对大气影响不太敏感，并使用蓝波段解决了初步大气校正后残留气溶胶引起的植被指数衰减问题。在处理土壤背景、大气和饱和噪声对盐分监测的影响方面，ERSSI中蓝绿波段组合可以有效应对干扰。从评价结果可以看出，待测区域的覆盖类型为秸秆残留，盐侵状态包括盐壳，常规植被敏感波段不适用于干旱区秸秆残留农田。为了解决植物残留的影响，使用SWIR₁代替NIR，可以有效识别缺水环境下植被的细胞结构和土壤含水量，并利用绿和SWIR₁波段的组合可以表征铁含量和土壤水分含量低和土壤状况，从而间接表明盐渍化。待测区域土壤盐渍化主要含有钙离子和硫酸根离子，两者在1443-1745nm范围内反射强烈，SWIR₁波段对这些物质敏感，使用SWIR₁波段为监测新疆秸秆残留农田土壤盐渍化提供了机会。

步骤4)干旱区秸秆残留的农田土壤盐分估算：根据步骤3)得到的相关性评价结果，以Sentinel-2遥感影像的二阶导数变换计算得到的增强型秸秆残留盐分指数(Enhanced Residues Soil Salinity Index，ERSSI)为自变量，以步骤2)得到的Q1中土壤电导率数据集为因变量，在R语言中构建线性多项式回归方程：

EC=0.00072×ERSSI²+0.072×ERSSI+0.92

其中，EC为干旱区秸秆残留农田的土壤表层盐分含量(dS/m)。该多项式回归方程在步骤2)得到的Q1(建模训练集)中精度R²为0.63；同时，将增强型秸秆残留盐分指数(ERSSI)和多项式回归方程在步骤2)得到的Q2(独立验证集)中进行独立验证，验证精度R²为0.64；

根据增强型秸秆残留盐分指数(ERSSI)的计算公式，基于Sentinel-2光学遥感影像，对待测区域的每一个栅格点进行ERSSI指数计算，将构建的多项式回归模型应用到待测区域中，得到待测区域的逐像元的土壤盐分含量。

选取新疆维吾尔自治区的棉花种植田块作为待测区域，以本实施方式估算得到的2020年土壤盐分含量分布图如图2所示，空间分辨率为10米。需说明的是，该土壤盐分含量分布图为以电导率作为表征的相对盐分含量分布图，如需转换为绝对盐分含量的分布图，则将每个像元的值按照电导率与土壤盐分含量的转换公式进行转换即可。

以上所述的实施例只是本发明的一种较佳的方案，然其并非用以限制本发明。有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型。因此凡采取等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。

Claims (1)

1.一种利用遥感构建指数估算秸秆残留农田土壤盐分的方法，其特征在于，步骤如下：

1)、获取待测区域在待估算时间下对应的Sentinel-2卫星数据；所述待测区域的土地类型为干旱区秸秆残留农田；

2)、对Sentinel-2卫星数据中的多光谱影像进行预处理，得到不同波段的表面反射率数据，再对每个栅格的光谱进行二阶微分变换得到二阶微分变换光谱；

3)、基于待测区域内每个栅格的二阶微分变换光谱计算对应的三维指数，三维指数形式为

其中ERSSI为高敏感指数，Green为Sentinel-2卫星数据的绿Green波段，Blue为Sentinel-2卫星数据的蓝Blue波段，SWIR₁为Sentinel-2卫星数据的短波红外1波段SWIR 1；

4)、将待测区域内每个栅格的三维指数作为自变量，基于预先针对待测区域构建的回归模型估计每个栅格对应的土壤表层盐分含量，最终形成待测区域内秸秆残留农田土壤盐分含量的空间分布图；

所述回归模型为多项式回归模型，其形式为：

EC＝a×ERSSI²+b×ERSSI+c

式中：EC为土壤表层盐分含量，a、b、c分别为回归系数。

Images