CN113642663B - 一种卫星遥感影像水体提取方法 - Google Patents

Abstract

本发明提出一种卫星遥感影像水体提取方法，包括：确定样本位置，获取样本点在一定时间序列上的多源遥感影像集，分别完成去云处理后进行合成，获得去云合成影像；使用合成影像的波段数据，使用自定义的归一化差异水体指数(NDWI)进行计算得到指数值；将影像转为RGB图像，输入已经训练完成的分类网络模型，得到该图像的类激活图CAM；融合CAM和NDWI，通过阈值筛选得到marker，将RGB图像转为灰度图并使用随机游走算法进行分割，通过面积阈值筛选得到大面积水体分割结果。本方法通过对仅有图像级标注的遥感影像训练数据集进行学习与训练，融合NDWI和CAM进行分割，在消除对真实水体标注的依赖同时，也加强了水体提取的精度。

Description

一种卫星遥感影像水体提取方法

技术领域

本发明属于水体遥感影像技术领域，主要涉及一种结合水体指数和类激活图的卫星遥感影像的水体提取方法。

背景技术

遥感技术的发展使遥感影像在地理资源研究与调查中的重要性日益增长。在水体提取技术中，使用遥感数据获取水体信息已经成为目前主流的水体提取方式。在日益增长的精度要求下，对高性能的水体提取算法技术的需求也不断增大。目前已经有一些方法取得了不错的效果，然而由于地理、气候环境等因素的影响，目前已有的水体提取方法难以应对复杂的现实状况。因此，研究能基于遥感影像且足以应对多种现实状况的水体提取技术具有重要的意义与价值。

目前单波段阈值法、多波段谱间关系法、水体指数法等方法作为水体提取技术的基础方法，相继得到了应用，产生了一些衍生变种方案。其中，单波段阈值法主要利用红外波段光谱特征的差异，选取合适的阈值进行提取。多波段谱间关系法则利用不同波段的光谱差异，构造一种模型来提取水体。该两种方法在地物复杂多样或者波段相似的情形下不能取得较好的效果。而水体指数法是对多个光谱的特征分析，选取需要的多个波段，构造水体指数计算公式来实现水体提取，其准确度依赖于阈值的选定和计算方案。该方法能较大限度地抑制植被信息带来的干扰，却仍然难以应对城市空间或者特殊地貌下带来的如阴影、土壤成分等复杂要素的干扰。

总的来说，现有的基础方法难以应对复杂的情形，亟需新的方案来对现状作出改进。而近年来，计算机视觉领域正在快速发展，其中一些基于普通图像的深度学习技术值得思考。这给了分析遥感影像进行水体提取一些新的思路与方法。例如可以将遥感图像转换为灰度图像后利用灰度共生矩阵提取纹理特征，从而关注到更多的水体类别；可以使用神经网络的方法对颜色等特征进行进一步的提取映射，如利用语义分割的解码器结构实现像素级别的分割预测，从而对水体颜色信息更加敏感。但大多基于深度学习的方法需要像素级细粒度水体标注对模型进行有监督训练，而获得像素级细粒度水体标注需要耗费大量人力成本。

发明内容

本发明所需解决的技术问题是，如何在缺少像素级细粒度水体标注的情况下对遥感影像中大面积水体区域进行精准水体提取。

为了解决上述的技术问题，本发明所采用的技术方案是提供一种基于归一化差异水体指数和类激活图的遥感影像水体提取的方法，包括以下步骤：

步骤S1：确定样本位置，获取样本点在一定时间序列上的多源遥感影像集，分别完成去云处理后进行合成，获得去云合成影像；

步骤S2：使用合成影像的波段数据，使用自定义的水体指数进行计算，得到指数值；

步骤S3：将影像转为RGB图像，输入已经训练完成的分类网络模型，得到该图像的类激活图CAM；

步骤S4：将NDWI和CAM进行融合得到水体特征图WFM，并设定阈值进行二值筛选，得到该图像的marker；将RGB图像转为灰度图，然后使用随机游走算法对灰度图进行分割；对分割结果进行面积阈值筛选，得到大面积水域的分割结果。

进一步的，步骤S1中对薄云，可以使用同态滤波、光谱差异性或暗通道等方案进行过滤；对厚云，则需要采用不同时间段同位置的影像进行合成互补，重建去云影响，从而保证整个影像完整。

进一步的，步骤S2中所述水体指数为标准的归一化差分水体指数NDWI或其它变种水体指数，可根据具体情况进行选择；所述标准的归一化差分水体指数NDWI使用绿色和近红外(NIR)波段，其表达式为

进一步的，变种水体指数的一种为NWI:

其中Band₁为蓝光波段；Band₄近红外短波；Band₅为短波红外1；Band₇为短波红外2。

进一步的，步骤S3中的具体步骤为：

步骤3.1：选取红绿蓝波段，通过真彩色合成将影像转为RGB图像。

步骤3.2：将RGB图像输入已经训练完成的分类网络模型，取出分类网络的最后一个卷积层的输出，使用grad-CAM算法得到该图像的类激活图CAM。

更进一步的，所述分类网络模型的训练方法为：

训练集中的图像X带有图像级标签y，即表示图像中是否存在水体；将X逆时针旋转90°，得到X’；将X和X’分别输入结构相同的分类网络C和C’中，输出二分类结果o和o’，并与真实标签y比较，得出loss1进行反向传播优化参数；分别取出C和C’最后一个卷积层输出的特征图M和M’，分别使用grad-CAM算法得到cam和cam’，将cam’顺时针旋转90°，并与cam比较，得出loss2进行反向传播优化参数。

再进一步的，loss1和loss2的计算方式如下：

采用cross entropy作为loss1，其计算公式为loss1＝ylog(o)+(1-y)log(1-o)；采用mse作为loss2，其计算公式为

其中

表示归一化，ij表示第i行第j列的像素点。

grad-CAM算法是计算类激活图CAM的一种算法，其输入为分类网络最后一层卷积层的输出特征图M和类激活图CAM所对应的类的输出概率o，输出为类激活图CAM，即cam＝grad_CAM(M，o)。

进一步的，步骤S4中的具体步骤为：

步骤4.1：选取NDWI和CAM的权重α_D和α_C，将NDWI和CAM进行线性加权融合得到水体特征图，NDWI和CAM融合公式为WFM_ij＝α_DNDWI_ij+α_CCAM_ij，其中ij表示第i行第j列的像素点，也可以使用其他融合方式。

步骤4.2：设定阈值threshold_WFM进行二值筛选，得到该图像的marker，规定WFM小于threshold_WFM的部分为0，其余部分为1。

步骤4.3：将RGB图像转为灰度图Gray，转换公式为

其中R表示红色通道，G表示绿色通道，B表示蓝色通道，ij表示第i行第j列的像素点，

表示对x向下取整。然后使用随机游走算法random walk算法对灰度图Gray进行分割。

更进一步的，随机游走算法random walk算法是一种基于图论的分割算法；算法首先需要给出种子点，即标记点；输入为灰度图Gray和种子点marker，输出为分割结果Seg，即Seg＝randomWalk(Gray，marker)；其分割思想如下：

以图像的像素为图的顶点，相邻像素之间的四邻域关系为图的边，并根据像素属性及相邻像素之间特征的相似性定义图中各边的权值，以此构建网络图。然后根据指定的前景和背景标记，即前景物体和背景物体的种子像素，以边上的权重为转移概率，未标记像素节点为初始点，计算每个未标记节点首次到达各种子像素的概率；根据概率大小，划分未标记节点，得到最终分割结果。

步骤4.4：对分割结果中的每个连通多边形进行面积阈值threshold_S筛选，得到大面积水域的分割结果；面积阈值threshold_S可根据具体需要分割水体区域的大小选取。

本发明通过旋转的数据增强方式对没有真实水体标注的遥感影像训练数据集进行学习，利用类激活图CAM对水体区域进行初步的定位，使用阈值筛选出种子点并用随机游走算法进行水体的分割，得到水体提取结果。

本发明与现有技术相比有如下提升：

(1)本发明结合了传统的NDWI和类激活图CAM的方法，将两种水体特征相融合，提升了水体特征提取的效果，解决了只根据NDWI作为单一水体特征的特征信息不足的问题，使得水体提取的结果更加准确。

(2)本发明通过旋转的数据增强和CAM的增强约束的方法来对模型的训练进行监督，在不需要真实的水体标注的情况下能对水体进行准确定位。

(3)本发明使用深度学习可视化中类激活图CAM来生成水体特征图，使深度学习模型具有一定的可解释性。

附图说明

图1是本发明的方法流程图；

图2是本发明分类网络类激活图约束训练流程图；

图3是本发明的输入图像示意图；

图4是本发明的提取结果示意图。

具体实施方式

下面结合附图对本发明的具体实施方式以及工作原理作进一步详细说明。此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

实施例1

本实施例对鄱阳湖水域进行水体提取，使用的数据为Landsat遥感图像数据，如图3所示。

本发明设计了一种基于深度学习的航拍图像水体提取方法与系统，方法流程图如图1所示。

步骤S1：

选取一个时间序列上的原始多源遥感影像集，进行去云处理后进行合成。

本发明需要将遥感图像输入到去云合成模块中。

本实施例选用Landsat遥感图像数据，利用其质量评估波段(QA)波段进行去云，其具体步骤为：

首先提取给定影像的质量评估波段，在Landsat影像质量评估波段数据中，第三位表示云阴影，第五位表示云层，第六到七位表示其云层置信度。移除影像中第三、五位被界定为云，且置信度中等的像元，保存去云影像序列。

对整个去云影像序列，取序列上最后最新的像元，合并得到一张去云遥感图像。

步骤S2：

使用合成影像的波段数据，使用自定义的归一化差异水体指数进行计算。

本实施例采用新型水体指数NWI，选定阈值为-0.6。NWI指数结合了水体在近红外和中红外波段具有强吸收的特点，根据波段计算得到的NWI值大于等于-0.6的像元即认定为水体，值越大，置信度越高。

在Landsat TM中，结合NWI值的计算公式为

其中Band₁为BLUE，即蓝光波段，可用于水体穿透、分辨土壤植被等；Band₄为NIR，即近红外短波，可用于估算生物量、分辨潮湿土壤等；Band₅为SWIR1，即短波红外1，可用于分辨道路、裸露土壤、水，还能在不同植被之间有好的对比度，并且有较好的大气、云雾分辨能力；Band₇为SWIR2，即短波红外2，可用于岩石、矿物的分辨，也可用于辨识植被覆盖和湿润土壤。

步骤S3：

对遥感图像进行真彩色合成。真彩色合成是指多光谱遥感图像彩色合成处理时，选取与红、绿、蓝三种原色的三个波段，合成一个颜色近似于地面景物真实颜色的RGB图像。

本发明需要将RGB图像输入到训练好的分类模型中，所以要先对分类模型进行训练。

本实施例选用Alexnet模型。应当指出，本实施例选用Alexnet仅用于说明本发明，但本发明的分类网络不限于Alexnet。

对Alexnet分类网络进行训练，其具体步骤为：

训练集中的图像X带有图像级标签y，即表示图像中是否存在水体，0表示无水体，1表示有水体；将X逆时针旋转90°，得到X’；将X和X’分别输入结构相同的Alexnet分类网络C和C’中，输出二分类结果o和o’，并与真实标签y比较，得出loss1进行反向传播优化参数；分别取出C和C’最后一个卷积层输出的特征图M和M’，分别使用grad-CAM算法得到cam和cam’，将cam’顺时针旋转90°，并与cam比较，得出loss2进行反向传播优化参数，CAM约束部分的流程如图2所示。

loss1和loss2的计算方式如下：

采用cross entropy作为loss1，其计算公式为loss1＝ylog(o)+(1-y)log(1-o)；采用mse作为loss2，其计算公式为

cam′_ij)²，其中

表示归一化，ij表示第i行第j列的像素点。

Alexnet训练完成后，将RGB图像输入已经训练完成的分类网络模型，取出Alexnet的最后一个卷积层输出的特征图M，使用grad-CAM算法得到该图像的类激活图CAM。

grad-CAM算法是计算类激活图CAM的一种算法，其输入为分类网络最后一层卷积层的输出特征图M和类激活图CAM所对应的类的输出概率o，输出为类激活图CAM，即cam＝grad_CAM(M，o)。

步骤S4：

选取NWDI和CAM的权重α_D和α_C，将NDWI和CAM进行线性加权融合得到水体特征图。NDWI和CAM融合公式为WFM_ij＝α_DNDWI_ij+α_CCAM_ij，其中ij表示第i行第j列的像素点。本实施例选取α_D＝0.6，α_C＝0.4。

设定阈值threshold_WFM进行二值筛选，得到该图像的marker，规定WFM小于threshold_WFM的部分为0，其余部分为1。本实施例选取threshold_WFM的值为0.6。

将RGB图像转为灰度图Gray，转换公式为

其中R表示红色通道，G表示绿色通道，B表示蓝色通道，ij表示第i行第j列的像素点，

表示对x向下取整。然后使用随机游走算法random walk算法对灰度图Gray进行分割，输入为灰度图Gray和种子点marker，输出为分割结果Seg，即Seg＝randomWalk(Gray，marker)。

随机游走算法random walk算法是一种基于图论的分割算法，其分割思想如下：

以图像的像素为图的顶点，相邻像素之间的四邻域关系为图的边，并根据像素属性及相邻像素之间特征的相似性定义图中各边的权值，以此构建网络图。然后根据指定的前景和背景标记，即前景物体和背景物体的种子像素，以边上的权重为转移概率，未标记像素节点为初始点，计算每个未标记节点首次到达各种子像素的概率。根据概率大小，划分未标记节点，得到最终分割结果。

对分割结果中的每个连通多边形进行面积阈值threshold_s筛选，得到大面积水域的分割结果。面积阈值threshold_S可根据具体需要分割水体区域的大小选取。本实施例选取threshold_S的值为10²，得到的水体提取结果如附图4所示。从结果可以看出，无论是大面积水体还是细小的水体，本发明都能够精确地进行提取。

本发明结合传统水体指数和深度学习的方法，分类网络在训练过程中能够学到水体特征，与传统的基于水体指数的水体提取方法相比，具有更强的泛化能力和更高的提取精度。同时，分类网络的训练不依赖真实的水体标注，与普通的深度学习方法相比，节省了获取真实水体标注的昂贵成本。此外，CAM的使用也使模型具有了一定的可解释性。

以上实施例仅用于说明本发明，而并非对本发明的限制。应当指出，本领域技术人员在不脱离本发明的精神和范围的情况下还可以做出各种变化和变型，因此所有不脱离本发明范畴所做出的改进和修改也在本发明的保护范围之内。

Claims (8)

1.一种卫星遥感影像水体提取方法，其特征在于，包括以下步骤：

步骤S1：确定样本位置，获取样本点在一定时间序列上的多源遥感影像集，分别完成去云处理后进行合成，获得去云合成影像；

步骤S2：使用合成影像的波段数据，使用水体指数进行计算，得到指数值；

步骤S3：将影像转为RGB图像，输入已经训练完成的分类网络模型，得到该图像的类激活图CAM；

步骤S4：将NDWI和CAM进行融合得到水体特征图WFM，并设定阈值进行二值筛选，得到该图像的marker；将RGB图像转为灰度图，使用随机游走算法对灰度图进行分割；对分割结果进行面积阈值筛选，得到大面积水域的分割结果。

2.根据权利要求1所述的卫星遥感影像水体提取方法，其特征在于，步骤S2中所述水体指数为标准的归一化差分水体指数NDWI或其它变种水体指数。

3.根据权利要求2所述的卫星遥感影像水体提取方法，其特征在于，变种水体指数的一种为NWI:

其中Band₁为蓝光波段；Band₄近红外短波；Band₅为短波红外1；Band₇为短波红外2。

4.根据权利要求1所述的卫星遥感影像水体提取方法，其特征在于，步骤S3所述分类网络模型的训练方法为：训练集中的图像X带有图像级标签y，表示图像中是否存在水体；将X逆时针旋转90°，得到X’；将X和X’分别输入结构相同的分类网络C和C’中，输出二分类结果o和o’，并与真实标签y比较，得出loss1进行反向传播优化参数；分别取出C和C’最后一个卷积层输出的特征图M和M’，分别使用grad-CAM算法得到cam和cam’，将cam’顺时针旋转90°，并与cam比较，得出loss2进行反向传播优化参数。

5.根据权利要求4所述的卫星遥感影像水体提取方法，其特征在于，所述loss1和loss2的计算方法为：采用cross entropy作为loss1，其计算公式为loss1＝ylog(o)+(1-y)log(1-o)；采用mse作为loss2，其计算公式为

其中

表示归一化，ij表示第i行第j列的像素点。

6.根据权利要求1所述的卫星遥感影像水体提取方法，其特征在于，其特征在于，步骤S3得到类激活图的方法为：取出分类网络的最后一个卷积层的输出，使用grad-CAM算法得到该图像的类激活图CAM；grad-CAM算法其输入为分类网络最后一层卷积层的输出特征图M和类激活图CAM所对应的类的输出概率o，输出为类激活图CAM，即cam＝grad_CAM(M，o)。

7.根据权利要求1所述的卫星遥感影像水体提取方法，其特征在于，步骤S4中：选取NWDI和CAM的权重，将NDWI和CAM进行融合得到水体特征图。

8.根据权利要求1所述的卫星遥感影像水体提取方法，其特征在于，步骤S4中使用随机游走算法对灰度图进行分割：算法首先给出种子点，即标记点；输入为灰度图Gray和种子点marker，输出为分割结果Seg，即Seg＝randomWalk(Gray，marker)。

Images