CN113011372B - 一种盐碱化土地自动监测和识别方法 - Google Patents

Abstract

本发明提出一种盐碱化土地自动监测和识别方法，属于环境监测及环境生态保护环境监测及环境生态保护领域。该方法获取研究区域多个采样点的土壤含盐量组成样本数据集；根据各采样点的经纬度提取距离采样月份最近的遥感影像中该经纬度对应栅格的若干波段值，通过构建不同的波段组合形式得到波段组合集合；然后构建机器学习土壤含盐量反演模型，利用样本数据集对模型进行训练，从波段组合集合中选出适合反演研究区域土壤含盐量的波段组合形式构成敏感因子集及最优反演模型；利用该最优反演模型，得到该研究区域内在选定月份的土壤含盐量反演结果。本发明可实际应用于干旱半干旱地区地表土层盐碱化程度的监测和不同程度盐碱化土壤的识别，为当地灌溉规划和农业管理等提供有效的参考依据。

Description

一种盐碱化土地自动监测和识别方法

技术领域

本发明属于环境监测及环境生态保护领域，具体涉及一种盐碱化土地自动监测和识别方法。

背景技术

目前，土壤盐碱化已成为影响农业生产和人类生活的几大自然现象之一，尤其在干旱半干旱地区，由于气候干旱，没有足够的降水淋洗土壤中的盐分，在蒸发和植被根系吸力等影响下，土壤盐碱化现象较为严重；此外，频繁的农业灌溉活动会引发次生盐碱化问题，加剧了对农业生产的威胁。利用遥感影像进行流域范围的盐碱地监测和识别，可以指导地区水资源规划与管理，尤其是可有效指导当地灌溉和农业用水和土壤改良，对生态环境改善和地区可持续发展具有重要意义。

传统的盐碱化监测方法多采用实地测量方式，实地点尺度测量电导率、土壤含盐量等变量，但成本高、周期长、易受地形和地势限制，精度不高且无法实现流域范围的监测。遥感技术的发展为该问题提供了一种低成本、简单易行的方法。多光谱影像、高光谱影像、雷达数据等已逐渐应用于土壤盐碱化的识别和监测领域。但目前已有成果一种土壤盐渍化程度的高光谱遥感判断方法(CN109738380A)研究采用偏最小二乘回归模型，对于高维度、非线性问题无法很好地描述；一种土壤盐碱化的雷达遥感监测方法(200910089204.6)中的土壤介电模型参数较多，调参过程复杂且多靠人为经验确定参数；主要依靠经验公式进行计算，计算过程复杂且可靠性较差。综上，对于土壤盐碱化识别和监测工作来说，一种简单、成本低的基于遥感影像的机器学习方法非常有必要。

发明内容

本发明的目的是为克服目前已有预报方法不能简单、准确地自动监测和识别盐碱化土壤的问题，提出一种盐碱化土地自动监测和识别方法。本发明结合少量的实测土壤含盐量数据，基于遥感数据和机器学习方法，建立自动识别和监测盐碱地演变的模型；可实际应用于干旱半干旱地区作为土地盐碱度的参考值，为当地水资源分配与管理尤其干旱区灌溉规划、农业用水管理、土壤改良和生态修复等提供有效的参考依据。

本发明提出一种盐碱化土地自动监测和识别方法，该方法首先获取研究区域若干个采样点的土壤含盐量数据组成样本数据集；根据每个采样点的经纬度提取该研究区域距离采样月份最近的遥感影像中该经纬度对应栅格的若干个波段值，并通过构建不同的波段组合形式得到波段组合集合；然后构建机器学习土壤含盐量反演模型，利用样本数据集对模型进行训练，从波段组合集合中选出适合反演研究区域土壤含盐量的波段组合形式构成敏感因子集，并得到最优反演模型；获取研究区域在待反演月份对应的遥感影像，逐幅计算该影像中各栅格的敏感因子集；将所有栅格的敏感因子集输入最优反演模型，最终得到该研究区域内的土壤含盐量的反演结果。该方法包括以下步骤：

1)确定研究区域，获取该研究区域n个样本的实测土壤含盐量数据组成样本数据集SC，SC＝[sc₁,sc₂,…,sc_n]，其中sc_j为第j个样本的实测土壤含盐量，j＝1,2,…,n；

获取该研究区域距离采样月份最近的遥感影像；

2)根据每个样本采样点的经纬度提取步骤1)获取的遥感影像中该经纬度对应的栅格的若干个波段值；

3)对每一个栅格，分别从提取的波段值中进行任意组合运算，得到各栅格的M种波段组合形式构成的波段组合集合BC，BC＝(bc₁,bc₂,…,bc_M)，其中bc_m为任一样本采样点对应栅格的第m种波段组合；每个栅格M种波段组合形式均相同；将波段组合集合BC作为当前波段组合集合；

4)从步骤3)得到的波段组合集合BC中选出适合反演研究区域土壤含盐量的波段组合形式构成敏感因子集，该因子集包含z个因子即z种波段组合形式；具体步骤如下：

4-1)记因子序号i＝0；波段组合遍历序号m＝1；采用k-折交叉验证方法，将集合SC中的所有样本随机平均分成k份样本集，记为：SC₁,SC₂,…,SC_k；其中，k大于等于3，且n能被k整除；

设定机器学习土壤含盐量反演模型参数的优化范围和步长；

4-2)初始化一个机器学习土壤含盐量反演模型并作为当前模型；

4-3)令验证集序号q＝1，则当前模型的验证集SC_q为SC₁，将剩余(k-1)份样本集共同构成当前模型的训练集；

4-4)将训练集的每个样本对应的i个因子各自的波段组合数据向量和第m种波段组合数据向量bc_mt共同作为当前模型的输入，bc_mt表示训练集中所有样本对应的当前波段组合集合中第m种波段组合值组成的向量，下标t表示训练集；

当前模型输出为训练集各样本对应的土壤含盐量的反演值向量SC_t’；

4-5)将训练集中各样本的实测土壤含盐量组成的向量记为SC_t，利用SC_t与SC_t’进行对比，计算本次训练相关系数R₁和训练相对偏差rBias₁评价当前模型：

式中，sc_j表示训练集中第j个样本的实测土壤含盐量数据，sc_ave表示训练集中所有样本实测土壤含盐量的平均值；sc'_j表示训练集中第j个样本对应的土壤含盐量的反演值；sc'_ave表示训练集中所有样本对应的土壤含盐量的反演值的平均值；n₁表示训练集样本总数，n₁＝((k-1)/k)n；

4-6)将验证集的每个样本对应的i个因子各自的波段组合数据向量和第m种波段组合数据向量BC(m)_v作为当前模型的输入，BC(m)_v表示验证集中所有样本对应的当前波段组合集合中第m种波段组合组成的向量，下标v表示验证集；当前模型输出为验证集各样本的土壤含盐量的反演值向量SC_v’；

将验证集每个样本的实测土壤含盐量组成的向量记为SC_v，利用SC_v与SC_v’进行对比，计算本次验证相关系数R₂和验证相对偏差rBias₂评价当前模型：

式中，sc_l表示SC_v中第l个样本的实测土壤含盐量数据，l＝1,2,…,n₂；sc_ave2表示SC_v中所有样本实测土壤含盐量数据的平均值；sc'_l表示SC_v’中第l个样本对应的反演值，sc'_ave2表示SC_v’中所有样本对应的反演值的平均值；n₂表示验证集样本总数，n₂＝(1/k)n；

4-7)对q进行判定：

若q＜k，则令验证集序号q＝q+1，更新当前模型的验证集SC_q，同时将剩余(k-1)份样本集共同构成更新后的训练集，重复步骤4-4)-4-6)，计算训练集和验证集更新后的当前模型的R₁、R₂、rBias₁和rBias₂；

若q≥k，则验证集遍历完毕，得到第m种波段组合在不同验证集和训练集组合下分别对应的当前模型的R₁、R₂、rBias₁和rBias₂；然后将当前模型运行k次得到的R₁、R₂、rBias₁和rBias₂各自的平均值分别作为第m种波段组合在当前模型的最终评价指标R₁、R₂、rBias₁和rBias₂的值，然后进入步骤4-8)；

4-8)根据模型参数的优化范围和步长，调整模型参数，更新当前模型，然后重新返回4-3)，直至得到第m种波段组合在每一组模型参数对应的模型下的最终评价指标R₁、R₂、rBias₁和rBias₂的值；

4-9)令m＝m+1，然后重复步骤4-2)至4-9)，直至得到当前波段组合集合中每种波段组合在每一组模型参数对应的模型下的最终评价指标R₁、R₂、rBias₁和rBias₂的值，然后进入步骤4-10)；

4-10)从步骤4-9)得到的所有R₂中选出R₂的绝对值最高值对应的模型作为第i+1个因子选取过程中得到的最优模型，并将该模型对应的波段组合作为第i+1个因子；同时，从当前波段组合集合中删除该波段组合，更新当前波段组合集合；

4-11)对i进行判定：

若i≤1，则令i＝i+1，然后重新返回步骤4-2)；

若i＞1，则计算ΔR₂，ΔR₂＝R_2,i-R_2,i-1，判定ΔR₂是否已经小于设定阈值：若否，则令i＝i+1，然后重新返回4-2)；若是，则z＝i+1，当前选出的i+1个因子即为敏感因子集的z个因子，将第i+1个因子选取过程中得到的最优模型作为最终的土壤含盐量反演模型，并记录该反演模型的参数；

5)获取研究区域在待反演月份对应的遥感影像，对遥感影像逐幅进行所覆盖范围土壤含盐量的反演，得到研究区域内的土壤含盐量的反演结果；具体步骤如下：

5-1)获取研究区域在待反演月份对应的遥感影像，该遥感影像与步骤1)的遥感影像类型相同；逐幅读取遥感影像中每个栅格的各个波段值，对每个栅格按照步骤4-11)选出的z个因子对应的波段组合形式进行波段之间的组合计算，得到各栅格对应的敏感因子集；

5-2)将步骤5-1)得到的各栅格对应的敏感因子集输入步骤4-11)得到的反演模型中，逐幅得到遥感影像所覆盖范围土壤含盐量的反演结果；

5-3)将步骤5-2)得到的反演结果进行影像拼接和裁剪提取，得到研究区域内的土壤含盐量的反演结果。

本发明的特点及有益效果在于：

本发明提出的一种盐碱化土地自动监测和识别方法，利用易获取的Landsat 8OLI影像数据，借助机器学习方法，巧妙地找出地表含盐量最敏感的波段组合形式，并建立可反演含盐量的机器学习模型，从而探索出一种低成本、易操作的盐碱地自动监测和识别方法。基于此，可实际应用于干旱半干旱区流域尺度的土壤盐碱化程度监测，合理制定农业灌溉计划，减少次生盐碱化问题，为当地水资源分配与管理、农业灌溉规划和生态修复等提供有效的参考依据。

具体实施方式

本发明提出一种盐碱化土地自动监测和识别方法，下面结合具体实施例进一步详细说明如下。

本发明提出一种盐碱化土地自动监测和识别方法，该方法首先获取研究区域若干个采样点的土壤含盐量数据组成样本数据集；根据每个采样点的经纬度提取该研究区域距离采样月份最近的遥感影像中该经纬度对应栅格的若干个波段值，并通过构建不同的波段组合形式得到波段组合集合；然后构建机器学习土壤含盐量反演模型，利用样本数据集对模型进行训练，从波段组合集合中选出适合反演研究区域土壤含盐量的波段组合形式构成敏感因子集，并得到最优反演模型；获取研究区域在待反演月份对应的遥感影像，逐幅计算该影像中各栅格的敏感因子集；将所有栅格的敏感因子集输入最优反演模型，最终得到该研究区域内的土壤含盐量的反演结果。该方法包括以下步骤：

1)确定研究区域，获取该研究区域n个样本的土壤含盐量数据组成样本数据集SC(每个样本的采样点在研究区域均匀分布)，SC＝[sc₁,sc₂,…,sc_n]，其中sc_j为第j个样本的实测土壤含盐量，j＝1,2,…,n。一般地，n≥50；获取该研究区域距离采样月份最近的遥感影像，本实施例采用Landsat8 OLI level2影像；

2)根据每个样本采样点的经纬度提取步骤1)获取的遥感影像中该经纬度对应的栅格的若干个波段值(本实施在每个栅格提取前七个波段，即B1～B7)；

3)对每一个栅格，分别从提取的波段值中进行任意两个、三个或四个波段的加、减、乘、除运算，得到各栅格的M种波段组合形式构成的波段组合集合BC，BC＝(bc₁,bc₂,…,bc_M)，其中bc_m为任一样本采样点对应栅格的第m种波段组合；每个栅格M种波段组合形式均相同；将波段组合集合BC作为当前波段组合集合；

4)从步骤3)得到的波段组合集合BC中选出适合反演研究区域土壤含盐量的波段组合形式构成敏感因子集，该因子集包含z(

下同，z值需通过模型确定)个因子即z种波段组合形式，也就是机器学习模型有z个输入变量；具体步骤如下：

4-1)记因子序号i＝0；波段组合遍历序号m＝1；采用k-折交叉验证方法，将集合SC中的所有样本随机平均分成k份样本集，记为：SC₁,SC₂,…,SC_k；其中，k值自己设定，大于等于3，且n能被k整除；将每份样本集SC_q轮流作为当前模型的验证集对应的输出(即一组模型参数下，运行k次)；同样地，剩余(k-1)份样本集(SC₁,SC₂,…,SC_q-1,SC_q+1,…,SC_k)共同构成模型的训练集，q＝1,2,…,k；

设定机器学习土壤含盐量反演模型参数的优化范围和步长；对于神经网络模型，则需要确定模型层数(一个输入层、一个或多个隐含层、一个输出层)、传递函数和激活函数；若为支持向量机，需要确定核函数及模型的惩罚因子、松弛变量和核函数中的参数；若为其他模型，则根据所用模型确定需要率定或优化的参数；

4-2)初始化一个机器学习土壤含盐量反演模型(简称机器学习模型)并作为当前模型(机器学习模型为“黑箱”模型，确定输入因子后，通过调参进行模型优化，模型会输出对目标变量SC的模拟值向量SC’)

4-3)令验证集序号q＝1，则当前模型的验证集SC_q为SC₁，将剩余(k-1)份样本集共同构成当前模型的训练集；

4-4)将训练集的每个样本对应的i个因子各自的波段组合数据向量和第m种波段组合数据向量bc_mt共同作为当前模型的输入，bc_mt表示训练集中所有样本对应的当前波段组合集合中第m种波段组合值组成的向量，下标t表示训练集；

当前模型输出为训练集各样本对应的土壤含盐量的反演值向量SC_t’；

4-5)将训练集中各样本的土壤含盐量组成的向量记为SC_t，利用SC_t与SC_t’进行对比，计算本次训练相关系数R₁和训练相对偏差rBias₁评价当前模型：

式中，sc_j表示训练集中第j个样本的土壤含盐量数据，sc_ave表示训练集中所有样本实测土壤含盐量的平均值；sc'_j表示训练集中第j个样本对应的土壤含盐量的反演值；sc'_ave表示训练集中所有样本对应的土壤含盐量的反演值的平均值；n₁表示训练集样本总数，n₁＝((k-1)/k)n；

4-6)将验证集的每个样本对应的i个因子各自的波段组合数据向量和第m种波段组合数据向量BC(m)_v作为当前模型的输入，BC(m)_v表示验证集中所有样本对应的当前波段组合集合中第m种波段组合组成的向量，下标v表示验证集；当前模型输出为验证集各样本的土壤含盐量的反演值向量SC_v’；

将验证集每个样本的实测土壤含盐量组成的向量记为SC_v，利用SC_v与SC_v’进行对比，计算本次验证相关系数R₂和验证相对偏差rBias₂评价当前模型：

式中，sc_l表示SC_v中第l个样本的实测土壤含盐量数据，l＝1,2,…,n₂；sc_ave2表示SC_v中所有样本实测土壤含盐量数据的平均值；sc'_l表示SC_v’中第l个样本对应的土壤含盐量的反演值，sc'_ave2表示SC_v’中所有样本对应的土壤含盐量的反演值的平均值；n₂表示验证集样本总数，n₂＝(1/k)n；

4-7)对q进行判定：

若q＜k，则令验证集序号q＝q+1，更新当前模型的验证集SC_q，同时将剩余(k-1)份样本集共同构成模型更新后的训练集，重复步骤4-4)-4-6)，计算训练集和验证集更新后的当前模型的R₁、R₂、rBias₁和rBias₂；

若q≥k，则验证集遍历完毕，得到第m种波段组合在不同验证集和训练集组合下分别对应的当前模型的R₁、R₂、rBias₁和rBias₂；然后将当前模型运行k次得到的评价指标R₁、R₂、rBias₁和rBias₂各自的平均值分别作为第m种波段组合在当前模型的最终评价指标R₁、R₂、rBias₁和rBias₂的值，然后进入步骤4-8)；

4-8)根据模型参数的优化范围和步长，调整模型参数，更新当前模型，然后重新返回4-3)，直至得到第m种波段组合在每一组模型参数对应的模型下的最终评价指标R₁、R₂、rBias₁和rBias₂的值；

4-9)令m＝m+1，然后重复步骤4-2)至4-9)，直至得到当前波段组合集合中每种波段组合在每一组模型参数对应的模型下的最终评价指标R₁、R₂、rBias₁和rBias₂的值，然后进入步骤4-10)；

4-10)从步骤4-9)得到的所有R₂中选出R₂的绝对值最高值对应的模型作为第i+1个因子选取过程中得到的最优模型，并将该模型对应的波段组合作为选出的第i+1个输入变量(因子)，即选出了第i个指标；同时，从当前波段组合集合中删除该波段组合形式，更新当前波段组合集合；

4-11)对i进行判定：

若i≤1，则令i＝i+1，然后重新返回步骤4-2)；

若i＞1，则计算ΔR₂，ΔR₂＝R_2,i-R_2,i-1，判定ΔR₂是否已经小于设定阈值(本实施例为0.1)：若否，则令i＝i+1，然后重新返回4-2)；若是，则z＝i+1，当前选出的i+1个因子即为敏感因子集的z个因子，也就是最终的z个输入变量；将第i+1个因子选取过程中得到的最优模型作为最终的土壤含盐量反演模型，并记录该反演模型的参数设置；

5)获取研究区域在待反演月份对应的遥感影像，对遥感影像逐幅进行所覆盖范围土壤含盐量的反演，得到研究区域内的土壤含盐量的反演结果；具体步骤如下：

5-1)获取研究区域在待反演月份对应的遥感影像，该遥感影像与步骤1)的遥感影像类型相同(本实施例为Landsat 8OLI影像数据)，逐幅读取影像中每个栅格的各个波段值，对每个栅格按照步骤4-11)选出的z个因子对应的波段组合形式进行波段之间的组合计算，得到各栅格对应的敏感因子集；

5-2)将步骤5-1)得到的各栅格对应的敏感因子集输入到步骤4-11)得到的反演模型中(即第z个输入变量对应的最优模型)，逐幅得到遥感影像所覆盖范围土壤含盐量的反演结果；

5-3)将步骤5-2)逐幅的反演结果进行影像拼接和裁剪提取，得到研究区域内的土壤含盐量的反演结果，从而实现对该研究区域土壤盐碱化的监测和识别。

下面结合一个具体实施例对本发明进一步详细说明如下。

本实施例提出的一种盐碱化土地自动监测和识别方法，包括以下步骤：

1)确定研究区域；本实施例研究区域可为一个灌区或一个流域A；

获取研究区域n个样本的实测土壤含盐量数据组成样本数据集SC，在本实施例中，在A流域共选取60个样本采样点进行土块采集和实验室分析，得到每个样本的实测土壤含盐量；获取该研究区域距离采样月份最近的遥感影像，本实施例采用Landsat8 OLI level2影像；

2)根据每个样本采样点的经纬度提取步骤1)获取的遥感影像中该经纬度对应的栅格的若干个波段值(本实施在每个栅格提取前七个波段，即B1～B7)；

本实施例Landsat 8 OLI影像前七个波段信息如表1所示：

表1 Landsat 8 OLI影像的光谱信息

获取方法具体如下：2-1)在美国地质调查局(United States GeologicalSurvey，USGS)网站(https://earthexplorer.usgs.gov/)获取研究区域的Landsat 8 OLIlevel2影像；

2-2)将步骤2-1)得到每幅影像B1～B7的波段值读取成R×L×7的三维矩阵形式；R和L分别表示每幅影像所含栅格的行数和列数；

2-3)根据每个样本采样点的经纬度信息提取影像上相应栅格的B1～B7值；

3)对每一个栅格，分别从提取的波段值中进行任意两个、三个或四个波段的加、减、乘、除运算，得到各栅格的M种波段组合形式构成的波段组合集合BC，BC＝(bc₁,bc₂,…,bc_M)，其中bc_m为任一样本采样点对应栅格的第m种波段组合；每个栅格M种波段组合形式均相同；将波段组合集合BC作为当前波段组合集合；

本实施例中，M种波段组合形式如表2所示；

表2本实施例的波段组合形式

其中，下标a,b,c,d代表波段的序号；

4)从步骤3)得到的波段组合集合BC中选出适合反演研究区域土壤含盐量的波段组合形式构成敏感因子集，该因子集包含z(

下同，z值需通过模型确定)个因子即z种波段组合形式，也就是机器学习模型有z个输入变量；具体步骤如下：4-1)记因子序号i＝0；波段组合遍历序号m＝1；采用k-折交叉验证方法，将集合SC中的所有样本随机平均分成k份样本集，记为：SC₁,SC₂,…,SC_k；其中，k值自己设定，大于等于3，且n能被k整除；将每份样本集SC_q轮流作为当前模型的验证集对应的输出(即一组模型参数下，运行k次)；同样地，剩余(k-1)份样本集(SC₁,SC₂,…,SC_q-1,SC_q+1,…,SC_k)共同构成模型的训练集，q＝1,2,…,k；

设定机器学习土壤含盐量反演模型参数的优化范围和步长；本案例以支持向量机模型为例，核函数采用径向基核函数，模型需要率定的参数是惩罚因子C和核函数中的参数g；取C和g的优化范围为[0.5,20]，步长均取0.5；进行网格化搜索；

4-2)构建一个机器学习土壤含盐量反演模型(简称机器学习模型)并作为当前模型(机器学习模型为“黑箱”模型，确定输入因子后，通过调参进行模型优化，模型会输出对目标变量SC的模拟值向量SC’)

4-3)令验证集序号q＝1，则当前模型的验证集SC_q为SC₁，将剩余(k-1)份样本集共同构成当前模型的训练集；

4-4)将训练集的每个样本对应的i个因子各自的波段组合数据向量和第m种波段组合数据向量bc_mt共同作为当前模型的输入，bc_mt表示训练集中所有样本对应的当前波段组合集合中第m种波段组合值组成的向量，下标t表示训练集；

当前模型输出为训练集各样本对应的土壤含盐量的反演值向量SC_t’；

4-5)将训练集中各样本的实测土壤含盐量组成的向量记为SC_t，利用SC_t与SC_t’进行对比，计算本次训练相关系数R₁和训练相对偏差rBias₁评价当前模型：

式中，sc_j表示训练集中第j个样本的实测土壤含盐量数据，sc_ave表示训练集中所有样本实测土壤含盐量的平均值；sc'_j表示训练集中第j个样本对应的土壤含盐量的反演值；sc'_ave表示训练集中所有样本对应的土壤含盐量的反演值的平均值；n₁表示训练集样本总数，n₁＝((k-1)/k)n；

4-6)将验证集的每个样本对应的i个因子各自的波段组合数据向量和第m种波段组合数据向量BC(m)_v作为当前模型的输入，BC(m)_v表示验证集中所有样本对应的当前波段组合集合中第m种波段组合组成的向量，下标v表示验证集；当前模型输出为验证集各样本土壤含盐量的反演值向量SC_v’；

将验证集每个样本的实测土壤含盐量组成的向量记为SC_v，利用SC_v与SC_v’进行对比，计算本次验证相关系数R₂和验证相对偏差rBias₂评价当前模型：

式中，sc_l表示SC_v中第l个样本的实测土壤含盐量数据，l＝1,2,…,n₂；sc_ave2表示SC_v中所有样本实测土壤含盐量数据的平均值；sc'_l表示SC_v’中第l个样本对应的土壤含盐量的反演值，sc'_ave2表示SC_v’中所有样本对应的土壤含盐量的反演值的平均值；n₂表示验证集样本总数，n₂＝(1/k)n；

4-7)对q进行判定：

若q＜k，则令验证集序号q＝q+1，更新当前模型的验证集SC_q，同时将剩余(k-1)份样本集共同构成模型更新后的训练集，重复步骤4-4)-4-6)，计算训练集和验证集更新后的当前模型的R₁、R₂、rBias₁和rBias₂；

若q≥k，则验证集遍历完毕，得到第m种波段组合在不同验证集和训练集组合下分别对应的当前模型的R₁、R₂、rBias₁和rBias₂；然后将当前模型运行k次得到的评价指标R₁、R₂、rBias₁和rBias₂各自的平均值分别作为第m种波段组合在当前模型的最终评价指标R₁、R₂、rBias₁和rBias₂的值，然后进入步骤4-8)；

4-8)根据模型参数的优化范围和步长，调整模型参数，更新当前模型，然后重新返回4-3)，直至得到第m种波段组合在每一组模型参数对应的模型下的最终评价指标R₁、R₂、rBias₁和rBias₂的值；

4-9)令m＝m+1，然后重复步骤4-2)至4-9)，直至得到当前波段组合集合中每种波段组合在每一组模型参数对应的模型下的最终评价指标R₁、R₂、rBias₁和rBias₂的值，然后进入步骤4-10)；

4-10)从步骤4-9)得到的所有R₂中选出R₂的绝对值最高值对应的模型作为第i+1个因子选取过程中得到的最优模型，并将该模型对应的波段组合作为选出的第i+1个输入变量(因子)，即选出了第i个指标；同时，从当前波段组合集合中删除该波段组合形式，更新当前波段组合集合；

4-11)对i进行判定：

若i≤1，则令i＝i+1，然后重新返回步骤4-2)；

若i＞1，则计算ΔR₂，ΔR₂＝R_2,i-R_2,i-1，判定ΔR₂是否已经小于设定阈值(本实施例为0.1)：若否，则令i＝i+1，然后重新返回4-2)；若是，则z＝i+1，当前选出的i+1个因子即为最终的z个输入变量，将第i+1个因子选取过程中得到的最优模型作为最终的土壤含盐量反演模型，并记录该反演模型的参数设置；

5)获取研究区域在待反演月份对应的Landsat 8 OLI影像数据，对遥感影像逐幅进行所覆盖范围土壤含盐量的反演，得到研究区域内的土壤含盐量的反演结果；具体步骤如下：

5-1)获取研究区域在待反演月份对应的Landsat 8 OLI影像数据，逐幅读取影像中每个栅格的各个波段值，对每个栅格按照步骤4-11)选出的z个波段组合形式进行波段之间的组合计算，得到各栅格对应的敏感因子集；

5-2)将步骤5-1)得到的各栅格对应的敏感因子集输入到步骤4-11)得到的反演模型中(即第z个输入变量对应的最优模型)，逐幅遥感影像得到所覆盖范围土壤含盐量的反演结果；

5-3)将步骤5-2)逐幅的反演结果进行影像拼接和裁剪提取，得到研究区域内的土壤含盐量的反演结果，从而实现对该研究区域土壤盐碱化的监测和识别。

Claims (1)

1.一种盐碱化土地自动监测和识别方法，其特征在于，该方法首先获取研究区域若干个采样点的土壤含盐量数据组成样本数据集；根据每个采样点的经纬度提取该研究区域距离采样月份最近的遥感影像中该经纬度对应栅格的若干个波段值，并通过构建不同的波段组合形式得到波段组合集合；然后构建机器学习土壤含盐量反演模型，利用样本数据集对模型进行训练，从波段组合集合中选出适合反演研究区域土壤含盐量的波段组合形式构成敏感因子集，并得到最优反演模型；获取研究区域在待反演月份对应的遥感影像，逐幅计算该影像中各栅格的敏感因子集；将所有栅格的敏感因子集输入最优反演模型，最终得到该研究区域内的土壤含盐量的反演结果；该方法包括以下步骤：

1)确定研究区域，获取该研究区域n个样本的实测土壤含盐量数据组成样本数据集SC，SC＝[sc₁,sc₂,…,sc_n]，其中sc_j为第j个样本的实测土壤含盐量，j＝1,2,…,n；

获取该研究区域距离采样月份最近的遥感影像；

2)根据每个样本采样点的经纬度提取步骤1)获取的遥感影像中该经纬度对应的栅格的若干个波段值；

3)对每一个栅格，分别从提取的波段值中进行任意组合运算，得到各栅格的M种波段组合形式构成的波段组合集合BC，BC＝(bc₁,bc₂,…,bc_M)，其中bc_m为任一样本采样点对应栅格的第m种波段组合；每个栅格M种波段组合形式均相同；将波段组合集合BC作为当前波段组合集合；

4)从步骤3)得到的波段组合集合BC中选出适合反演研究区域土壤含盐量的波段组合形式构成敏感因子集，该因子集包含z个因子即z种波段组合形式；具体步骤如下：

4-1)记因子序号i＝0；波段组合遍历序号m＝1；采用k-折交叉验证方法，将集合SC中的所有样本随机平均分成k份样本集，记为：SC₁,SC₂,…,SC_k；其中，k大于等于3，且n能被k整除；

设定机器学习土壤含盐量反演模型参数的优化范围和步长；

4-2)初始化一个机器学习土壤含盐量反演模型并作为当前模型；

4-3)令验证集序号q＝1，则当前模型的验证集SC_q为SC₁，将剩余(k-1)份样本集共同构成当前模型的训练集；

4-4)将训练集的每个样本对应的i个因子各自的波段组合数据向量和第m种波段组合数据向量bc_mt共同作为当前模型的输入，bc_mt表示训练集中所有样本对应的当前波段组合集合中第m种波段组合值组成的向量，下标t表示训练集；

当前模型输出为训练集各样本对应的土壤含盐量的反演值向量SC_t’；

4-5)将训练集中各样本的实测土壤含盐量组成的向量记为SC_t，利用SC_t与SC_t’进行对比，计算本次训练相关系数R₁和训练相对偏差rBias₁评价当前模型：

式中，sc_j表示训练集中第j个样本的实测土壤含盐量数据，sc_ave表示训练集中所有样本实测土壤含盐量的平均值；sc′_j表示训练集中第j个样本对应的土壤含盐量的反演值；sc′_ave表示训练集中所有样本对应的土壤含盐量的反演值的平均值；n₁表示训练集样本总数，n₁＝((k-1)/k)n；

4-6)将验证集的每个样本对应的i个因子各自的波段组合数据向量和第m种波段组合数据向量BC(m)_v作为当前模型的输入，BC(m)_v表示验证集中所有样本对应的当前波段组合集合中第m种波段组合组成的向量，下标v表示验证集；当前模型输出为验证集各样本的土壤含盐量的反演值向量SC_v’；

将验证集每个样本的实测土壤含盐量组成的向量记为SC_v，利用SC_v与SC_v’进行对比，计算本次验证相关系数R₂和验证相对偏差rBias₂评价当前模型：

式中，sc_l表示SC_v中第l个样本的实测土壤含盐量数据，l＝1,2,…,n₂；sc_ave2表示SC_v中所有样本实测土壤含盐量数据的平均值；sc′_l表示SC_v’中第l个样本对应的反演值，sc′_ave2表示SC_v’中所有样本对应的反演值的平均值；n₂表示验证集样本总数，n₂＝(1/k)n；

4-7)对q进行判定：

若q＜k，则令验证集序号q＝q+1，更新当前模型的验证集SC_q，同时将剩余(k-1)份样本集共同构成更新后的训练集，重复步骤4-4)-4-6)，计算训练集和验证集更新后的当前模型的R₁、R₂、rBias₁和rBias₂；

若q≥k，则验证集遍历完毕，得到第m种波段组合在不同验证集和训练集组合下分别对应的当前模型的R₁、R₂、rBias₁和rBias₂；然后将当前模型运行k次得到的R₁、R₂、rBias₁和rBias₂各自的平均值分别作为第m种波段组合在当前模型的最终评价指标R₁、R₂、rBias₁和rBias₂的值，然后进入步骤4-8)；

4-8)根据模型参数的优化范围和步长，调整模型参数，更新当前模型，然后重新返回4-3)，直至得到第m种波段组合在每一组模型参数对应的模型下的最终评价指标R₁、R₂、rBias₁和rBias₂的值；

4-9)令m＝m+1，然后重复步骤4-2)至4-9)，直至得到当前波段组合集合中每种波段组合在每一组模型参数对应的模型下的最终评价指标R₁、R₂、rBias₁和rBias₂的值，然后进入步骤4-10)；

4-10)从步骤4-9)得到的所有R₂中选出R₂的绝对值最高值对应的模型作为第i+1个因子选取过程中得到的最优模型，并将该模型对应的波段组合作为第i+1个因子；同时，从当前波段组合集合中删除该波段组合，更新当前波段组合集合；

4-11)对i进行判定：

若i≤1，则令i＝i+1，然后重新返回步骤4-2)；

若i＞1，则计算ΔR₂，ΔR₂＝R_2,i-R_2,i-1，判定ΔR₂是否已经小于设定阈值：若否，则令i＝i+1，然后重新返回4-2)；若是，则z＝i+1，当前选出的i+1个因子即为敏感因子集的z个因子，将第i+1个因子选取过程中得到的最优模型作为最终的土壤含盐量反演模型，并记录该反演模型的参数；

5)获取研究区域在待反演月份对应的遥感影像，对遥感影像逐幅进行所覆盖范围土壤含盐量的反演，得到研究区域内的土壤含盐量的反演结果；具体步骤如下：

5-1)获取研究区域在待反演月份对应的遥感影像，该遥感影像与步骤1)的遥感影像类型相同；逐幅读取遥感影像中每个栅格的各个波段值，对每个栅格按照步骤4-11)选出的z个因子对应的波段组合形式进行波段之间的组合计算，得到各栅格对应的敏感因子集；

5-2)将步骤5-1)得到的各栅格对应的敏感因子集输入步骤4-11)得到的反演模型中，逐幅得到遥感影像所覆盖范围土壤含盐量的反演结果；

5-3)将步骤5-2)得到的反演结果进行影像拼接和裁剪提取，得到研究区域内的土壤含盐量的反演结果。