CN115527108A - 一种基于多时相Sentinel-2快速识别水土流失人为扰动地块的方法 - Google Patents

Abstract

本发明提供了一种基于多时相Sentinel‑2快速识别水土流失人为扰动地块的方法，本发明基于中高空间分辨率影像和时间序列Sentinel‑2数据，大大提高了人为扰动区域时空分布监测的效率和提取结果的精度，经过特征优选和最优参数选择后的模型逼近不仅能确保数据不会过度拟合，而且可以得到性能更佳的预测模型；同时结合鲁棒性好，泛化能力强的随机森林算法，解决了人为扰动区域监测所面临的高时空数据不足、提取精度较低以及高效、准确提取人为扰动区域的有效方法较少等瓶颈问题。以上方法在保证精度的情况下，实现人为扰动区域时空分布精细化提取识别。

Description

一种基于多时相Sentinel-2快速识别水土流失人为扰动地块 的方法

技术领域

本发明涉及快速、高效、高精度遥感人为扰动区域监测领域，具体涉及一种基于多时相Sentinel-2快速识别水土流失人为扰动地块的方法。

背景技术

水土流失是当前中国最亟待解决的生态环境问题之一，据最新发布的中国水土保持公报，我国水土流失面积达269.27万平方公里，占全国国土面积的 28.05％。除自然因素外，农林开发建设、矿产资源开发、城镇建设等各类地表覆盖发生改变的人为扰动项目增多也是加剧水土流失的重要原因。然而，目前针对人为扰动区域的自动识别与提取研究较少，其监查方法仍主要依靠人工目视解译并结合实地勘查，该方法周期长且成本高，无法高效率地实现大范围人为扰动区域时空分布提取，存在诸多局限性。遥感技术具有大面积、同步性、周期性对地观测的优势，是实现长时序大尺度人为扰动区域提取的重要手段。如何充分挖掘遥感数据用于人为扰动区域的高效识别，加强对生产建设项目的监管，从而降低人为因素引起的水土流失，是目前急需解决的问题。

遥感技术为人为扰动区域识别提取的问题提供了一定的解决思路，中高分辨率遥感影像因其具有空间分辨率相对较高，纹理信息丰富的特点，被广泛应用于特定人为扰动区域制图。然而，截至目前，该方面的应用局限于小范围研究区及特定扰动类型，在面对大尺度区域和多种扰动类型时，中高分辨率影像数据量巨大，结合面向对象方法，也直接增大了影像分割、样本训练、类型判定的耗时，对于快速提取大尺度区域的人为扰动并不友好，难以满足实际工作中对生产建设项目的高效实时监管的需要。

相比之下，中低分辨率影像在大区域尺度上的快速提取更有优势，但影像时间分辨率和空间分辨率等条件难以精准完成人为扰动提取工作。例如， MODIS的空间分辨率有限，难以满足地块破碎化地区的精细化制图要求； Landsat影像尽管有了更高的空间分辨率，但其时间分辨率较低，难以通过完整的时序数据获取准确的扰动信息。欧空局哨兵系列卫星的出现为大范围、长时序、高精度的人为扰动监测提供了契机，Sentinel-2数据拥有10m的多光谱波段分辨率和双星2-5d的重返周期，提供陆地卫星相似的光谱信息和空间信息，能在保持较高提取精度的同时，大幅提高运算效率，有效解决中高分辨率卫星、 MODIS与Landsat在人为扰动区域快速提取上遇到的困境。目前，基于Sentinel-2 时序影像进行人为扰动提取区域提取的应用研究较少，仍缺乏覆盖范围广、精度高、类型全面的人为扰动高效提取方法，因此，有必要探究Sentinel-2在该问题上的制图方法与制图潜力。

因此，针对人为扰动区域识别提取目前面临的：依赖人工解译、大区域多类型扰动提取方法研究较少、高分影像结合面向对象方法运算效率低等问题，基于Sentinel-2时序数据，通过选择最优分类特征组合及最优分类参数，利用随机森林的方法实现高精度、低成本、大范围的人为扰动区域进行逐年识别提取，能够为后续人为扰动区域提取及水土流失相关问题提供理论支撑和方法参考。

发明内容

本发明的目的在于提供一种能够高效、快速、高精度进行人为扰动监测的基于随机森林和特征优选的人为扰动区域时空分布信息提取方法。

本发明提供的基于多时相Sentinel-2快速识别水土流失人为扰动地块的方法，包括如下步骤：

S1、获取单时相中高空间分辨率Sentinel-2遥感影像并进行正射校正、辐射校正以及几何配准数据预处理，得到包含研究区域的单时相中高空间分辨率 Sentinel-2遥感影像；

S2、多张所述包含研究区域的单时相中高空间分辨率Sentinel-2遥感影像形成L2A级别的多时相时间序列Sentinel-2影像，将L2A级别的多时相时间序列 Sentinel-2影像通过欧空局SNAP平台进行格式转换输出与批量裁剪处理，得到研究区域的中高空间分辨率时序数据；

S3、对于步骤S2中得到的研究区域的中高空间分辨率时序数据，计算归一化植被指数、归一化水体指数、归一化建筑指数、比值居民地指数、植被稀疏与旺盛期的NDVI差值、对比度、方差和熵，从而生成中高空间分辨率指数时间序列数据集，并结合时序数据中红边波段R、绿边波段G、蓝边波段B、近红外波段Nir，生成中高空间分辨率待选特征数据集；

S4、根据步骤S3得到的中高空间分辨率待选特征数据集，结合野外采样各地物的样本点，构建分类样本库并进行分层抽样，从而形成训练集和验证集；

S5、基于所述训练集形成训练数据和测试数据，采用所述训练数据先利用随机森林算法进行模型训练，然后通过测试数据的表现确定最优模型参数和最优特征组合；

S6、将所述最优模型参数和最优特征组合应用至研究区域，实现研究区域域植被、水体、耕地、不透水层和人为扰动地块的识别，完成研究区域土地利用类型的提取工作；

S7、利用所述验证集对步骤S6的提取结果进行精度评价，若是精度评价结果不满足预设的期望精度，则返回步骤S5，精度评价结果满足预设的期望精度，完成土地利用类型识别，提取其中的人为扰动区域，最终完成人为扰动区域时空分布提取工作。

进一步地，步骤S1中，在进行几何配准时，依托于开源数据，手动选择若干控制点进行几何纠正，从而实现几何配准。

进一步地，在将L2A级别的多时相时间序列Sentinel-2影像通过欧空局 SNAP平台进行格式转换输出与批量裁剪处理之前，还基于云量对影像进行筛选：通过整景云量小于5％或者研究区云量小于2％来进行数据的筛选；所利用到的多时相时间序列Sentinel-2影像是有蓝、绿、红以及近红外四个波段组成，空间分辨率为10米。

进一步地，归一化植被指数、归一化水体指数、归一化建筑指数、比值居民地指数、植被稀疏与旺盛期的NDVI差值、对比度、方差和熵的计算公式如下：

归一化植被指数NDVI：

NDVI＝(Nir-Red)/(Nir+Red)

归一化水体指数NDWI：

NDWI＝(Green-Nir)/(Green+Nir)

归一化建筑指数NDBI：

NDBI＝(Mir-Nir)/(Mir+Nir)

比值居民地指数RRI：

RRI＝Blue/Nir

植被稀疏与旺盛期的NDVI差值dNDVI：

dNDVI＝NDVI_旺盛期-NDVI_稀疏期

对比度Contrast：

方差Variance：

熵Entropy：

式中，Nir、Red、Green、Mir、Blue分别为哨兵影像的近红外波段、红光波段、绿光波段、中红外波段和蓝色波段；i、j是像元在图像中的行列坐标，P (i，j)是灰度联合概率矩阵，μ为像元均值。

进一步地，所述步骤S5中随机森林具体包括：

对中高空间分辨率待选特征数据集中的12个特征进行优选，采用所述训练数据利用基于CART的随机森林方法计算出待选特征的相对重要性得分，并对12个待选特征依据相对重要性进行排序；依据重要性从高到低选取排名位于前 1/3的特征作为基础特征，并进行组合，得到待验证的特征基础组合；分别依据重要性从大到小，以此在前一组合的基础上增加一个重要性较高的特征，最终得到9组待验证特征组合；再利用随机森林算法产生的袋外数据，对每组待验证特征组合进行OOB误差评估，并进行内部的交叉验证；得到各特征组合的 OOB误差，其中误差最小的组合即为对应模型参数下的初步最优特征组合；

利用测试数据对不同模型参数下训练后的随机森林模型进行准确率的验证，得到最优模型参数，最优模型参数对应的的初步最优特征组合即为最终的优特征组合。

进一步地，所述模型参数包括：随机树的数量，叶子节点的数量以及最大叶深度。

进一步地，步骤S7进行精度评价具体采用指标选择总体准确率accuracy、用户精度Precision、生产者精度Recall以及Kappa系数四个指标进行评价，具体计算公式如下：

accuracy＝(TP+TN)/(TP+FN+FP+TN)

precision＝TP/(TP+FP)

recall＝TP/(TP+FN)

Kappa＝(p-p_e)/(1-p_e)

其中，

p_e＝(a₁*b₁+a₂*b₂+…+a_i*b_i)/(n*n)

式中，TP为实际是正样本预测为正样本的样本数，FP为实际是负样本预测为正样本的样本数，TN为实际是负样本预测为负样本的样本数，FN为实际是正样本预测为负样本的样本数；p_o代表总体分类精度，a_i表示i类地物的真实样本个数，b_i代表样本的预测值，n代表样本总数，p表示所有样本中分类正确的概率。

本发明提供的技术方案带来的有益效果是：本发明基于中高空间分辨率影像和时间序列Sentinel-2数据，大大提高了人为扰动区域时空分布监测的效率和提取结果的精度，经过特征优选和最优参数选择后的模型逼近不仅能确保数据不会过度拟合，而且可以得到性能更佳的预测模型；同时结合鲁棒性好，泛化能力强的随机森林算法，解决了人为扰动区域监测所面临的高时空数据不足、提取精度较低以及高效、准确提取人为扰动区域的有效方法较少等瓶颈问题。以上方法在保证精度的情况下，实现人为扰动区域时空分布精细化提取识别。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将对本发明实施方式作进一步地描述。

本发明的基于多时相Sentinel-2快速识别水土流失人为扰动地块的方法，包括如下步骤：

S1、获取单时相中高空间分辨率Sentinel-2遥感影像并进行正射校正、辐射校正以及几何配准数据预处理，得到包含研究区域的单时相中高空间分辨率 Sentinel-2遥感影像。辐射校正通过不同传感器的定标系数来完成；正射校正利用全球30米数字高程模型数据辅助完成；几何配准需要基于谷歌地球、天地图等中高分辨率无偏移辅助开源数据，手动选择地面控制点完成。

S2、多张所述包含研究区域的单时相中高空间分辨率Sentinel-2遥感影像形成L2A级别的多时相时间序列Sentinel-2影像，将L2A级别的多时相时间序列 Sentinel-2影像通过欧空局SNAP平台进行格式转换输出与批量裁剪处理，得到研究区域的中高空间分辨率时序数据。

其中格式转换输出可以输出tif文件。多时相时间序列Sentinel-2影像是指在年周期内获取不同月份的遥感影像，时相范围应考虑研究区不同作物类型的生长物候时间。在将L2A级别的多时相时间序列Sentinel-2影像通过欧空局 SNAP平台进行格式转换输出与批量裁剪处理之前，还基于云量对影像进行筛选：通过整景云量小于5％或者研究区云量小于2％来进行数据的筛选；所利用到的多时相时间序列Sentinel-2影像是有蓝、绿、红以及近红外四个波段组成，空间分辨率为10米。

S3、对于步骤S2中得到的研究区域的中高空间分辨率时序数据，计算归一化植被指数、归一化水体指数、归一化建筑指数、比值居民地指数、植被稀疏与旺盛期(3月与8月)的NDVI差值、对比度、方差和熵，从而生成中高空间分辨率指数时间序列数据集，并结合时序数据中红边波段R、绿边波段G、蓝边波段B、近红外波段Nir，一共12个特征，生成中高空间分辨率待选特征数据集。

归一化植被指数、归一化水体指数、归一化建筑指数、比值居民地指数、植被稀疏与旺盛期的NDVI差值、对比度、方差和熵的计算公式如下：

归一化植被指数NDVI：

NDVI＝(Nir-Red)/(Nir+Red)

归一化水体指数NDWI：

NDWI＝(Green-Nir)/(Green+Nir)

归一化建筑指数NDBI：

NDBI＝(Mir-Nir)/(Mir+Nir)

比值居民地指数RRI：

RRI＝Blue/Nir

植被稀疏与旺盛期(3月与8月)的NDVI差值dNDVI：

dNDVI＝NDVI_旺盛期-NDVI_稀疏期

对比度Contrast：

方差Variance：

熵Entropy：

式中，Nir、Red、Green、Mir、Blue分别为哨兵影像的近红外波段、红光波段、绿光波段、中红外波段和蓝色波段；i、j是像元在图像中的行列坐标，P (i，j)是灰度联合概率矩阵，μ为像元均值。

S4、根据步骤S3得到的中高空间分辨率待选特征数据集，结合野外采样各地物的样本点，构建分类样本库并按照一定比例进行分层抽样，从而形成训练集和验证集。

S5、基于所述训练集形成训练数据和测试数据，采用所述训练数据先利用随机森林算法进行模型训练，然后通过测试数据的表现确定最优模型参数和最优特征组合。

其中，基于所述训练集，有放回地抽取70％作为训练数据，30％作为测试数据。

对中高空间分辨率待选特征数据集中的12个特征进行优选，采用所述训练数据利用基于CART的随机森林方法计算出待选特征的相对重要性得分，并对 12个待选特征依据相对重要性进行排序；依据重要性从高到低选取排名位于前 1/3的特征作为基础特征，并进行组合，得到待验证的特征基础组合；分别依据重要性从大到小，以此在前一组合的基础上增加一个重要性较高的特征，最终得到9组待验证特征组合；再利用随机森林算法产生的袋外数据(Out of Bag, OOB)，即利用Bootstrap抽样时剩余未被抽取的训练数据，对每组待验证特征组合进行OOB误差评估，并进行内部的交叉验证；得到各特征组合的OOB误差，其中误差最小的组合即为对应模型参数下的初步最优特征组合；

利用测试数据对不同模型参数下训练后的随机森林模型进行准确率的验证，得到最优模型参数，最优模型参数对应的的初步最优特征组合即为最终的优特征组合，模型参数包括：随机树的数量，叶子节点的数量以及最大叶深度。

S6、将所述最优模型参数和最优特征组合应用至研究区域，实现研究区域域植被、水体、耕地、不透水层和人为扰动地块的识别，完成研究区域土地利用类型(植被、水体、耕地、不透水层和人为扰动地块)的提取工作。

S7、利用所述验证集对步骤S6的提取结果进行精度评价，若是精度评价结果不满足预设的期望精度，则返回步骤S5，精度评价结果满足预设的期望精度，完成土地利用类型识别，提取其中的人为扰动区域，最终完成人为扰动区域时空分布提取工作。

进行精度评价具体采用指标选择总体准确率accuracy、用户精度Precision、生产者精度Recall以及Kappa系数四个指标进行评价，在四项指标均满足预设的期望精度时，从而完成土地利用类型识别。具体计算公式如下：

accuracy＝(TP+TN)/(TP+FN+FP+TN)

precision＝TP/(TP+FP)

recall＝TP/(TP+FN)

Kappa＝(p-p_e)/(1-p_e)

其中，

p_e＝(a₁*b₁+a₂*b₂+…+a_i*b_i)/(n*n)

式中，TP为实际是正样本预测为正样本的样本数，FP为实际是负样本预测为正样本的样本数，TN为实际是负样本预测为负样本的样本数，FN为实际是正样本预测为负样本的样本数；p_o代表总体分类精度，a_i表示i类地物的真实样本个数，b_i代表样本的预测值，n代表样本总数，p表示所有样本中分类正确的概率。

在本文中，应当理解，所述方位词的使用不应限制本申请请求保护的范围。

在不冲突的情况下，本文中上述实施例及实施例中的特征可以相互结合。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于多时相Sentinel-2快速识别水土流失人为扰动地块的方法，其特征在于，包括如下步骤：

S1、获取单时相中高空间分辨率Sentinel-2遥感影像并进行正射校正、辐射校正以及几何配准数据预处理，得到包含研究区域的单时相中高空间分辨率Sentinel-2遥感影像；

S2、多张所述包含研究区域的单时相中高空间分辨率Sentinel-2遥感影像形成L2A级别的多时相时间序列Sentinel-2影像，将L2A级别的多时相时间序列Sentinel-2影像通过欧空局SNAP平台进行格式转换输出与批量裁剪处理，得到研究区域的中高空间分辨率时序数据；

S3、对于步骤S2中得到的研究区域的中高空间分辨率时序数据，计算归一化植被指数、归一化水体指数、归一化建筑指数、比值居民地指数、植被稀疏与旺盛期的NDVI差值、对比度、方差和熵，从而生成中高空间分辨率指数时间序列数据集，并结合时序数据中红边波段R、绿边波段G、蓝边波段B、近红外波段Nir，生成中高空间分辨率待选特征数据集；

S4、根据步骤S3得到的中高空间分辨率待选特征数据集，结合野外采样各地物的样本点，构建分类样本库并进行分层抽样，从而形成训练集和验证集；

S5、基于所述训练集形成训练数据和测试数据，采用所述训练数据先利用随机森林算法进行模型训练，然后通过测试数据的表现确定最优模型参数和最优特征组合；

S6、将所述最优模型参数和最优特征组合应用至研究区域，实现研究区域域植被、水体、耕地、不透水层和人为扰动地块的识别，完成研究区域土地利用类型的提取工作；

S7、利用所述验证集对步骤S6的提取结果进行精度评价，若是精度评价结果不满足预设的期望精度，则返回步骤S5，精度评价结果满足预设的期望精度，完成土地利用类型识别，提取其中的人为扰动区域，最终完成人为扰动区域时空分布提取工作。

2.根据权利要求1所述的基于多时相Sentinel-2快速识别水土流失人为扰动地块的方法，其特征在于，步骤S1中，在进行几何配准时，依托于开源数据，手动选择若干控制点进行几何纠正，从而实现几何配准。

3.根据权利要求1所述的基于多时相Sentinel-2快速识别水土流失人为扰动地块的方法，其特征在于，在将L2A级别的多时相时间序列Sentinel-2影像通过欧空局SNAP平台进行格式转换输出与批量裁剪处理之前，还基于云量对影像进行筛选：通过整景云量小于5％或者研究区云量小于2％来进行数据的筛选；所利用到的多时相时间序列Sentinel-2影像是有蓝、绿、红以及近红外四个波段组成，空间分辨率为10米。

4.根据权利要求1所述的基于多时相Sentinel-2快速识别水土流失人为扰动地块的方法，其特征在于，归一化植被指数、归一化水体指数、归一化建筑指数、比值居民地指数、植被稀疏与旺盛期的NDVI差值、对比度、方差和熵的计算公式如下：

归一化植被指数NDVI：

NDVI＝(Nir-Red)/(Nir+Red)

归一化水体指数NDWI：

NDWI＝(Green-Nir)/(Green+Nir)

归一化建筑指数NDBI：

NDBI＝(Mir-Nir)/(Mir+Nir)

比值居民地指数RRI：

RRI＝Blue/Nir

植被稀疏与旺盛期的NDVI差值dNDVI：

dNDVI＝NDVI_旺盛期-NDVI_稀疏期

对比度Contrast：

方差Variance：

熵Entropy：

式中，Nir、Red、Green、Mir、Blue分别为哨兵影像的近红外波段、红光波段、绿光波段、中红外波段和蓝色波段；i、j是像元在图像中的行列坐标，P(i，j)是灰度联合概率矩阵，μ为像元均值。

5.根据权利要求1所述的基于多时相Sentinel-2快速识别水土流失人为扰动地块的方法，其特征在于，所述步骤S5中随机森林具体包括：

对中高空间分辨率待选特征数据集中的12个特征进行优选，采用所述训练数据利用基于CART的随机森林方法计算出待选特征的相对重要性得分，并对12个待选特征依据相对重要性进行排序；依据重要性从高到低选取排名位于前1/3的特征作为基础特征，并进行组合，得到待验证的特征基础组合；分别依据重要性从大到小，以此在前一组合的基础上增加一个重要性较高的特征，最终得到9组待验证特征组合；再利用随机森林算法产生的袋外数据，对每组待验证特征组合进行OOB误差评估，并进行内部的交叉验证；得到各特征组合的OOB误差，其中误差最小的组合即为对应模型参数下的初步最优特征组合；

利用测试数据对不同模型参数下训练后的随机森林模型进行准确率的验证，得到最优模型参数，最优模型参数对应的的初步最优特征组合即为最终的优特征组合。

6.根据权利要求5所述的基于多时相Sentinel-2快速识别水土流失人为扰动地块的方法，其特征在于，所述模型参数包括：随机树的数量，叶子节点的数量以及最大叶深度。

7.根据权利要求1所述的基于多时相Sentinel-2快速识别水土流失人为扰动地块的方法，其特征在于，步骤S7进行精度评价具体采用指标选择总体准确率accuracy、用户精度Precision、生产者精度Recall以及Kappa系数四个指标进行评价，具体计算公式如下：

accuracy＝(TP+TN)/(TP+FN+FP+TN)

precision＝TP/(TP+FP)

recall＝TP/(TP+FN)

Kappa＝(p-p_e)/(1-p_e)

其中，

p_e＝(a1*b1+a2*b2+…+ai*bi)/(n*n)

式中，TP为实际是正样本预测为正样本的样本数，FP为实际是负样本预测为正样本的样本数，TN为实际是负样本预测为负样本的样本数，FN为实际是正样本预测为负样本的样本数；p_o代表总体分类精度，a_i表示i类地物的真实样本个数，b_i代表样本的预测值，n代表样本总数，p表示所有样本中分类正确的概率。