CN111539965A - 协同多源信息的城市地表水国产高分遥感自适应提取方法 - Google Patents

Abstract

本发明公开了一种协同多源信息的城市地表水国产高分遥感自适应提取方法，包括以下步骤：1)计算高分遥感影像的形态学大/小面积水体指数；2)计算中分遥感影像的改进水体指数；3)对水体指数计算结果，基于GlobeLand30产品的自适应分割，提取城市地表水结果W_MSWI，W_MLWI和W_MNDWI；4)对W_MLWI和W_MNDWI采用决策级融合后处理获得城市地表水后处理结果W_PLN以及对W_MNDWI采用纹理特征后处理获得城市地表水后处理结果W_PS；5)对步城市地表水后处理结果W_PLN和W_PS进行联合处理，获取最终的城市地表水结果。本发明实现城市地表水自动化信息提取，克服了国产高分遥感光谱信息不足的条件下，仍能在复杂的城市背景下进行自动化信息提取。

Description

协同多源信息的城市地表水国产高分遥感自适应提取方法

技术领域

本发明属于高分遥感图像处理及信息提取技术领域，尤其涉及一种协同多源信息的城市地表水国产高分遥感自适应提取方法。

背景技术

城市地表水与人类居住环境和生态系统服务密切相关，在城市规划，环境保护和热岛效应中发挥着重要的作用。然而，伴随着城市化进程的加剧及城市安全应急的需要，传统的监测方法(实地测绘和站点信息等)无法满足及时准确的获取城市大区域地表水信息的要求。随着国产高分遥感数据获取的时效性和便捷性均有大幅度的提升，使其成为城市地表水监测的重要数据源。

目前基于光学遥感的城市地表水提取方法基本为5种类型，具体包括：a)单波段阈值分割；b)多波段光谱水指数；c)有监督或无监督分类方法；d)线性解混方法以及e)基于上述方法的组合方法。在这些方法中，多波段光谱水体指数和阈值分割组合的方法因其具备操作便捷且提取结果可靠的特性被广泛研究与应用。当前应用效果较好的方法是一种基于8波段WorldView-2高分遥感影像的新型水体指数(Xie,C.,Huang,X.,Zeng,W.X.,andFang,X.2016.“A novel water index for urban high-resolution eight-bandWorldView-2imagery”.International Journal of Digital Earth,9,925-941)。这种方法在光谱信息的基础上引入了形态信息，有效的抑制了噪声，在提升了水体特征的显著性上有一定的优势。然而，由于国产高分辨率遥感影像的光谱信息匮乏，该水体指数方法未能应用于国产高分辨率遥感影像中。同时，人工的阈值选择难于保证城市地表水分割结果的稳定性与高效性。

发明内容

本发明解决的技术问题是：克服现有技术的不足，提供了一种协同多源信息的城市地表水国产高分遥感自适应提取方法，在发挥高分辨率遥感影像的空间信息优势基础上，协同中分辨率遥感影像的光谱信息优势，建立多尺度联合特征，在此基础上，引入GlobeLand30信息产品，提出一种协同多源信息的自适应学习方法，进行城市地表水的自动提取，实现城市地表水的快速监测。

本发明目的通过以下技术方案予以实现：一种协同多源信息的城市地表水国产高分遥感自适应提取方法，所述方法包括如下步骤：1)计算高分遥感影像的形态学大面积水体指数MLWI和高分遥感影像的形态学小面积水体指数MSWI；2)计算中分遥感影像的改进水体指数MNDWI；3)基于GlobeLand30产品的水体类型，对步骤1)的高分遥感影像的形态学大面积水体指数MLWI和高分遥感影像的形态学小面积水体指数MSWI和步骤2)的中分遥感影像的改进水体指数MNDWI，分别进行自适应分割，相应的获取城市地表水初始结果W_MSWI、城市地表水初始结果W_MLWI和城市地表水初始结果W_MNDWI；4)对城市地表水初始结果W_MLWI和城市地表水初始结果W_MNDWI采用决策级融合后处理获得城市地表水后处理结果W_PLN；对W_MSWI采用纹理特征后处理获得城市地表水后处理结果W_PS；5)对步骤4)的城市地表水后处理结果W_PLN和城市地表水后处理结果W_PS进行联合处理，获取最终的城市地表水结果W_FR。

上述协同多源信息的城市地表水国产高分遥感自适应提取方法中，在步骤1)中，计算高分遥感影像的形态学大面积水体指数MLWI和高分遥感影像的形态学小面积水体指数MSWI包括如下步骤：

MLWI＝EMPs_MCRmean-EMPs_WTHmean

MSWI＝EMPs_MCRmean+EMPs_WTHmean

其中，EMPs_MCRmean和EMPs_WTHmean的计算方法为：

其中，EMPs_MCR(i,r_Δs)和EMPs_WTH(i,r_Δs)的计算方法为：

EMPs_MCR(i,r_Δs)＝MCR⁽ⁱ⁾(r_Δs)

EMPs_WTH(i,r_Δs)＝WTH⁽ⁱ⁾(r_Δs)

其中，MCR⁽ⁱ⁾(r_Δs)和WTH⁽ⁱ⁾(r_Δs)的计算方法为：

WTH(r)＝NDWI-γ_NDWI(r)

其中，归一化水体指数NDWI的计算方法为：

式中，EMPs_MCRmean为扩展形态学闭重建轮廓序列的平均值；EMPs_WTHmean为扩展形态学形态学顶帽变换轮廓序列的平均值；EMPs_MCR(i,r_Δs)为扩展形态学闭重建轮廓序列，EMPs_WTH(i,r_Δs)为扩展形态学形态学顶帽变换轮廓序列；MCR为形态学闭重建，WTH为形态学顶帽变换；φ和γ分别表示NDWI图像的重建闭合和打开，i表示圆形结构体元素的数量，r表示圆形结构体元素的半径，Δs为圆形结构体元素数量之间的步长距离；ρ_Green和ρ_NIR分别为绿色波段和近红外波段的大气顶部反射率。

上述协同多源信息的城市地表水国产高分遥感自适应提取方法中，在步骤2)中，中分遥感影像的改进水体指数MNDWI为：

式中，ρ_Green和ρ_SWIR分别是绿色波段和SWIR波段的大气顶部反射率。

上述协同多源信息的城市地表水国产高分遥感自适应提取方法中，在步骤3)中，基于GlobeLand30产品的水体类型，对步骤1)的高分遥感影像的形态学大面积水体指数MLWI和高分遥感影像的形态学小面积水体指数MSWI和步骤2)的中分遥感影像的改进水体指数MNDWI，分别进行自适应分割，相应的获取城市地表水初始结果W_MSWI、城市地表水初始结果W_MLWI和城市地表水初始结果W_MNDWI包括如下步骤：

3.1)根据步骤1)的高分遥感影像的形态学大面积水体指数MLWI和高分遥感影像的形态学小面积水体指数MSWI和步骤2)的中分遥感影像的改进水体指数MNDWI，对GlobeLand30的水体类型产品进行K-Means聚类分析，从而获取水体样本信息S；

3.2)基于步骤3.1)的水体样本信息S，分别对高分遥感影像的形态学大面积水体指数MLWI、高分遥感影像的形态学小面积水体指数MSWI和中分遥感影像的改进水体指数MNDWI进行阈值分割，获取城市地表水初始结果W_MSWI、城市地表水初始结果W_MLWI和城市地表水初始结果W_MNDWI。

上述协同多源信息的城市地表水国产高分遥感自适应提取方法中，在步骤3.1)中，水体样本信息S通过以下公式得到：

式中，k是簇类型的数量；Ci是GlobeLand30的水体类型产品；t为GlobeLand30的水体类型产品集合的取值；F为下标，F＝MSWI、MLWI、MNDWI；DN为基于水体指数计算结果F的像素值；

是Ci水体类型产品的平均像素值向量。

上述协同多源信息的城市地表水国产高分遥感自适应提取方法中，在步骤3.2)中，城市地表水初始结果W_MSWI、城市地表水初始结果W_MLWI和城市地表水初始结果W_MNDWI通过以下公式得到：

W_F＝{DN_FI|DN_FI∈[Mean_FS-Stdev_FS,Mean_FS+Stdev_FS]}

其中，Mean_FS的计算方法为：

其中，Stdev_FS的计算方法为：

式中，Mean_FS和Stdev_FS分别为水体样本S的均值和标准差；n为水体样本S的数量；F为下标，F＝MSWI、MLWI、MNDWI；W_F为F所对应的城市地表水初始结果，即为W_MSWI，W_MLWI和W_MNDWI。

上述协同多源信息的城市地表水国产高分遥感自适应提取方法中，在步骤4)中，城市地表水后处理结果W_PLN的公式如下：

其中，O为城市地表水初始结果的对象；A是对象O的面积；T1为阈值参数。

上述协同多源信息的城市地表水国产高分遥感自适应提取方法中，在步骤4)中，城市地表水后处理结果W_PS的公式如下：

其中，GLCM的计算方法为：

式中，GLCM为灰度共生矩阵；N是灰度级的数量；g(i,j)是归一化对称灰度级共生矩阵；NDWI为归一化水体指数；O为城市地表水初始结果的对象；(i,j)为像素坐标；T2为阈值参数。

本发明与现有技术相比具有如下有益效果：

本发明的方法通过计算高分遥感影像的形态学水体指数，提升水体在高分遥感影像中的几何显著度，同时协同中分遥感影像的光谱信息，在此基础上补充水体的光谱信息，提升城市地表水的特征显著性；基于GlobeLand30产品，采用聚类的方法，对多水体特征进行自适应阈值分割，较大程度上克服了多维特征阈值难于确定的问题；采用决策级融合和纹理特征分别对不同类型的城市地表水进行类型后处理，去除噪声干扰，提升方法的鲁棒性，最终实现城市地表水的自动化信息提取。

附图说明

通过阅读下文优选实施方式的详细描述，各种其他的优点和益处对于本领域普通技术人员将变得清楚明了。附图仅用于示出优选实施方式的目的，而并不认为是对本发明的限制。而且在整个附图中，用相同的参考符号表示相同的部件。在附图中：

图1是本发明实施例提供的同多源信息的城市地表水国产高分遥感自适应提取方法的流程图；

图2(a)是本发明实施例提供的原始高遥感影像和GlobeLand30的水体类型产品；

图2(b)是本发明实施例提供的原始中遥感影像和GlobeLand30的水体类型产品；

图2(c)是本发明实施例提供的原始分遥感影像和GlobeLand30的水体类型产品；

图3(a)是本发明实施例提供的形态学大面积水体指数(MLWI)特征图；

图3(b)是本发明实施例提供的形态学小面积水体指数(MSWI)特征图；

图3(c)是本发明实施例提供的改进的水体指数(MNDWI)特征图；

图4是本发明实施例提供的基于K-Means聚类分析的GlobeLand30水体样本类型优化结果图；

图5(a)是本发明实施例提供的基于优化样本对MLWI进行自适应分割获取的城市地表水初始结果图；

图5(b)是本发明实施例提供的基于优化样本对MSWI进行自适应分割获取的城市地表水初始结果图；

图5(c)是本发明实施例提供的基于优化样本对MNDWI进行自适应分割获取的城市地表水初始结果图；

图6是本发明实施例提供的对W_MLWI和W_MNDWI采用决策级融合后处理获得城市地表水后处理结果图；

图7是本发明实施例提供的对W_MNDWI采用纹理特征后处理获得城市地表水后处理结果图；

图8是本发明实施例提供的城市地表水最终提取结果图。

具体实施方式

下面将参照附图更详细地描述本公开的示例性实施例。虽然附图中显示了本公开的示例性实施例，然而应当理解，可以以各种形式实现本公开而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本公开，并且能够将本公开的范围完整的传达给本领域的技术人员。需要说明的是，在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本发明。

本发明实施例中，所采用的高空间分辨率遥感影像和中空间分辨率遥感影像分别是2016年4月9日和2016年5月4日的高分一号高空间分辨率遥感影像和Landsat-8中空间分辨率遥感影像(见图2(a)和图2(b))。其中，高分一号卫星遥感数据的多光谱和全色分辨率分别为8米和2米，含有四个波段，分别是蓝波段(0.45-0.52μm),绿波段(0.52-0.59μm),红波段(0.63-0.69μm),近红外波段(0.77-0.89μm)，辐射量化等级为16bit。Landsat-8卫星遥感数据的多光谱和全色分辨率分别为30米和15米，主要应用六个波段，分别是蓝波段(0.45-0.52μm),绿波段(0.52-0.59μm),红波段(0.63-0.69μm),近红外波段(0.77-0.89μm)，中红外波段(1.56-1.66μm)以及全色波段(0.500–0.680)。此外，本实施例中所采用的城市地表水参考产品数据为2010年GlobaLand30土地覆盖产品中的水体产品结果(见图2(c))。

如图1所示，本发明提供了一种基于多源数据协同的城市地表水高分遥感自适应提取方法的具体步骤如下：

1)计算高分遥感影像的形态学大/小面积水体指数(MLWI/MSWI)的具体方法是，

MLWI＝EMPs_MCRmean-EMPs_WTHmean

MSWI＝EMPs_MCRmean+EMPs_WTHmean

其中，EMPs_MCRmean和EMPs_WTHmean的计算方法为：

其中，EMPs_MCR(i,r_Δs)和EMPs_WTH(i,r_Δs)的计算方法为：

EMPs_MCR(i,r_Δs)＝MCR⁽ⁱ⁾(r_Δs)

EMPs_WTH(i,r_Δs)＝WTH⁽ⁱ⁾(r_Δs)

其中，MCR⁽ⁱ⁾(r_Δs)和WTH⁽ⁱ⁾(r_Δs)的计算方法为：

WTH(r)＝NDWI-γ_NDWI(r)

其中，NDWI的计算方法为：

式中，EMPs_MCRmean和EMPs_WTHmean被定义为扩展形态学轮廓(EMP)的平均值；i表示MCR的数量，其参数(i_min，i_max和Δi)分别设置为(1，21和1)；WTH是圆形结构元素(SE)；φ和γ分别表示高分一号影像NDWI的重建闭合和打开，r表示圆形SE的半径，其参数与i的相同。本实施例大水体指数(MLWI)和形态学小水体指数(MSWI)的结果见图3(a)和图3(b)。

2)计算中分遥感影像的改进水体指数(MNDWI)的具体方法是，

式中，ρ_Green和ρ_SWIR分别是绿色波段和SWIR波段的大气顶部(ToA)反射率。本实施例改进的水体指数(MNDWI)特征的结果见图3(c)。

3)根据步骤1)和步骤2)的水体指数计算结果MLWI、MSWI和MNDWI，对GlobeLand30的水体类型产品进行K-Means聚类分析，从而获取纯度较高的水体样本信息S，计算方法为：

式中，k是簇类型的数量；Ci是GlobeLand30的水体类型产品；t为GlobeLand30的水体类型产品集合的取值；F为水体指数计算结果，其取值为{MSWI，MLWI和MNDWI}；DN为基于水体指数计算结果F的像素值；

是Ci水体类型产品的平均像素值向量。k簇聚集成水体和非水体两类，其中，水体样本信息即为S；

基于步骤3.1)的水体样本结果S，分别对水体指数计算结果MLWI、MSWI和MNDWI特征进行阈值分割，获取城市地表水初始结果W_MSWI，W_MLWI和W_MNDWI，计算方法为：

W_F＝{DN_FI|DN_FI∈[Mean_FS-Stdev_FS,Mean_FS+Stdev_FS]}

其中，Mean_FS的计算方法为：

其中，Stdev_FS的计算方法为：

式中，Mean_FS和Stdev_FS分别为水体样本S的均值和标准差；n为水体样本S的数量；F为水体指数计算结果，其取值为{MSWI，MLWI和MNDWI}；W_F为F取值的城市地表水初始结果，即为W_MSWI，W_MLWI和W_MNDWI。本实施例城市地表水初始W_MSWI，W_MLWI和W_MNDWI的结果见图5(a)，图5(b)和图5(c)。

4)采用面积指标对W_MLWI和W_MNDWI进行决策级融合后处理，并获得城市地表水后处理结果W_RLN，计算方法是：

式中，O为城市地表水初始结果的对象；A是对象O的面积；T1为阈值参数，本研究通过实验将其设置为0.7。本实施例的采用面积指标的大水体后处理结果图6；

采用纹理特征对W_MSWI进行后处理的计算方法是：

其中，GLCM的计算方法为：

式中，GLCM为灰度共生矩阵；N是灰度级的数量；g(i,j)是归一化对称灰度级共生矩阵；NDWI为归一化水体指数；O为城市地表水初始结果的对象；(I，j)为像素坐标；T2为阈值参数，本研究通过实验将其设置为0.4。本实施例的基于纹理特征的小水体后处理结果图7。

5)对4)的城市地表水后处理结果W_PLN和W_PS进行联合处理，获取最终的城市地表水结果W_FR，应用ArcGIS软件中的联合工具。本实施例的城市地表水最终提取结果图8。

本发明的方法利用多尺度形态学方法表达高分遥感的两类城市地表水空间形态特征，提出了形态学大/小面积水体指数，增强了不同形态的城市地表水表达能力；引入GlobeLand30参考产品数据，提出基于聚类的阈值区间自动选择方法，克服了多波段阈值难于确定的问题；由于国产高分遥感影像的光谱信息匮乏，导致城市地表水易于出现与裸地、建筑物阴影等混分问题，为此，协同具备光谱信息丰富的中分Landsat-8数据，采用决策级融合和纹理特征的组合方法，抑制不同类型噪声的干扰，最终现实城市场景下的城市地表水自动化提取。

本发明虽然已以较佳实施例公开如上，但其并不是用来限定本发明，任何本领域技术人员在不脱离本发明的精神和范围内，都可以利用上述揭示的方法和技术内容对本发明技术方案做出可能的变动和修改，因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化及修饰，均属于本发明技术方案的保护范围。

Claims (8)

1.一种协同多源信息的城市地表水国产高分遥感自适应提取方法，其特征在于，所述方法包括如下步骤：

1)计算高分遥感影像的形态学大面积水体指数MLWI和高分遥感影像的形态学小面积水体指数MSWI；

2)计算中分遥感影像的改进水体指数MNDWI；

3)基于GlobeLand30产品的水体类型，对步骤1)的高分遥感影像的形态学大面积水体指数MLWI和高分遥感影像的形态学小面积水体指数MSWI和步骤2)的中分遥感影像的改进水体指数MNDWI，分别进行自适应分割，相应的获取城市地表水初始结果W_MSWI、城市地表水初始结果W_MLWI和城市地表水初始结果W_MNDWI；

4)对城市地表水初始结果W_MLWI和城市地表水初始结果W_MNDWI采用决策级融合后处理获得城市地表水后处理结果W_PLN；对W_MSWI采用纹理特征后处理获得城市地表水后处理结果W_PS；

5)对步骤4)的城市地表水后处理结果W_PLN和城市地表水后处理结果W_PS进行联合处理，获取最终的城市地表水结果W_FR。

2.根据权利要求1所述的协同多源信息的城市地表水国产高分遥感自适应提取方法，其特征在于：在步骤1)中，计算高分遥感影像的形态学大面积水体指数MLWI和高分遥感影像的形态学小面积水体指数MSWI包括如下步骤：

MLWI＝EMPs_MCRmean-EMPs_WTHmean

MSWI＝EMPs_MCRmean+EMPs_WTHmean

其中，EMPs_MCRmean和EMPs_WTHmean的计算方法为：

其中，EMPs_MCR(i,r_Δs)和EMPs_WTH(i,r_Δs)的计算方法为：

EMPs_MCR(i,r_Δs)＝MCR⁽ⁱ⁾(r_Δs)

EMPs_WTH(i,r_Δs)＝WTH⁽ⁱ⁾(r_Δs)

其中，MCR⁽ⁱ⁾(r_Δs)和WTH⁽ⁱ⁾(r_Δs)的计算方法为：

WTH(r)＝NDWI-γ_NDWI(r)

其中，归一化水体指数NDWI的计算方法为：

式中，EMPs_MCRmean为扩展形态学闭重建轮廓序列的平均值；EMPs_WTHmean为扩展形态学形态学顶帽变换轮廓序列的平均值；EMPs_MCR(i,r_Δs)为扩展形态学闭重建轮廓序列，EMPs_WTH(i,r_Δs)为扩展形态学形态学顶帽变换轮廓序列；MCR为形态学闭重建，WTH为形态学顶帽变换；φ和γ分别表示NDWI图像的重建闭合和打开，i表示圆形结构体元素的数量，r表示圆形结构体元素的半径，Δs为圆形结构体元素数量之间的步长距离；ρ_Green和ρ_NIR分别为绿色波段和近红外波段的大气顶部反射率。

3.根据权利要求1所述的协同多源信息的城市地表水国产高分遥感自适应提取方法，其特征在于：在步骤2)中，中分遥感影像的改进水体指数MNDWI为：

式中，ρ_Green和ρ_SWIR分别是绿色波段和SWIR波段的大气顶部反射率。

4.根据权利要求1所述的协同多源信息的城市地表水国产高分遥感自适应提取方法，其特征在于：在步骤3)中，基于GlobeLand30产品的水体类型，对步骤1)的高分遥感影像的形态学大面积水体指数MLWI和高分遥感影像的形态学小面积水体指数MSWI和步骤2)的中分遥感影像的改进水体指数MNDWI，分别进行自适应分割，相应的获取城市地表水初始结果W_MSWI、城市地表水初始结果W_MLWI和城市地表水初始结果W_MNDWI包括如下步骤：

3.1)根据步骤1)的高分遥感影像的形态学大面积水体指数MLWI和高分遥感影像的形态学小面积水体指数MSWI和步骤2)的中分遥感影像的改进水体指数MNDWI，对GlobeLand30的水体类型产品进行K-Means聚类分析，从而获取水体样本信息S；

3.2)基于步骤3.1)的水体样本信息S，分别对高分遥感影像的形态学大面积水体指数MLWI、高分遥感影像的形态学小面积水体指数MSWI和中分遥感影像的改进水体指数MNDWI进行阈值分割，获取城市地表水初始结果W_MSWI、城市地表水初始结果W_MLWI和城市地表水初始结果W_MNDWI。

5.根据权利要求4所述的协同多源信息的城市地表水国产高分遥感自适应提取方法，其特征在于：在步骤3.1)中，水体样本信息S通过以下公式得到：

式中，k是簇类型的数量；Ci是GlobeLand30的水体类型产品；t为GlobeLand30的水体类型产品集合的取值；F为下标，F＝MSWI、MLWI、MNDWI；DN为基于水体指数计算结果F的像素值；

是Ci水体类型产品的平均像素值向量。

6.根据权利要求4所述的协同多源信息的城市地表水国产高分遥感自适应提取方法，其特征在于：在步骤3.2)中，城市地表水初始结果W_MSWI、城市地表水初始结果W_MLWI和城市地表水初始结果W_MNDWI通过以下公式得到：

W_F＝{DN_FI|DN_FI∈[Mean_FS-Stdev_FS,Mean_FS+Stdev_FS]}

其中，Mean_FS的计算方法为：

其中，Stdev_FS的计算方法为：

式中，Mean_FS和Stdev_FS分别为水体样本S的均值和标准差；n为水体样本S的数量；F为下标，F＝MSWI、MLWI、MNDWI；W_F为F所对应的城市地表水初始结果，即为W_MSWI，W_MLWI和W_MNDWI。

7.根据权利要求1所述的协同多源信息的城市地表水国产高分遥感自适应提取方法，其特征在于：在步骤4)中，城市地表水后处理结果W_PLN的公式如下：

其中，O为城市地表水初始结果的对象；A是对象O的面积；T1为阈值参数。

8.根据权利要求7所述的协同多源信息的城市地表水国产高分遥感自适应提取方法，其特征在于：在步骤4)中，城市地表水后处理结果W_PS的公式如下：

其中，GLCM的计算方法为：

式中，GLCM为灰度共生矩阵；N是灰度级的数量；g(i,j)是归一化对称灰度级共生矩阵；NDWI为归一化水体指数；O为城市地表水初始结果的对象；(i,j)为像素坐标；T2为阈值参数。

Images