CN110033499B - 一种基于Landsat OLI影像水体制图法 - Google Patents

Abstract

本发明公开了一种基于Landsat OLI遥感影像的区域性水体制图方法，包括以下内容：1）普适性水体指数模型构建；2）最佳水体指数阈值自动选取；3）区域性水体提取；4）水体提取结果验证；5）区域性水体制图。本发明在现有水体指数模型基础上，发展了一种适合于Landsat OLI遥感影像的区域水体提取指数模型，采用该模型波段运算后，可以将水体和非水体类型的值域进行有效分离，其中水体表现较为集中的值域范围，并采用了图像熵算法克服时空差异带来的节地物反射率异质性，该算法自动选取最佳水体提取阈值，减少人工确定阈值的不确定性和主观性，提高了水体制图精度和效率。

Description

一种基于Landsat OLI影像水体制图法

技术领域

本发明涉及的是一种区域性水体遥感制图方法，特别是一种基于Landsat OLI遥感影像的区域性水体制图方法，属于遥感制图技术领域。

背景技术

随着现代水体制图方法研究的深入和进步，水体提取指数模型的建构方法也在不断提升和改进：1996年Mcfeeters第一次利用绿光(Green，0.52-0.60um)波段和近红外(NIR，0.76-0.96um)波段构建了归一化水体指数(NDWI)，2004年Roger和Kearney利用红光(Red，0.62-0.69um)波段和短波红外(SWIR1， 1.55～1.75μm)波段提出了另一种水体提取的指数模型，该模型可以有效鉴别水体和裸土；2005年徐涵秋教授利用绿光(Green，0.52～0.60μm)和短波红外(SWIR1， 1.55～1.75μm)构建了归一化差异水体指数模型(MNDWI)，该模型可以有效鉴别水体和低反射率的人工建筑(如老旧房屋和柏油马路等)；2014年Feyisa等利用绿光(Green，0.525-0.600μm)波段、短波红外(SWIR1：1.560-1.660μm，SWIR2：2.100-2.300μm)波段、近红外(0.845-0.885μm)波段和蓝光(Blue，0.450-0.515μm) 波段构建了自动水体提取指数模型，该模型可以有效鉴别水体和阴影，适合在城市地区和山地地区的水体提取。现有主要水体提取指数模型如表1所示：

表1现有主要水体提取指数模型

然而，上述水体提取指数模型的建构存在如下缺陷：

(1)现有水体指数模型较少考虑到不同地理空间的遥感影像对水体反射值的影响，大部分水体指数模型阈值选择采用人工方式，导致水体提取结果具有较强主观性，且阈值试探过程费时费工；

(2)现有水体指数模型大多没有考虑到季节变化对水体反射值的影响，由于不同季节相同地物的反射率存在一定的差异，采用不变的水体指数阈值提取水体时，将导致提取的水体精度不一致，波动性大；

(3)现有指数模型阈值选取时主观性强，水体提取精度随选取阈值的变化而变化。

Landsat OLI遥感影像是Landsat系列卫星中Landsat 8搭载的传感器--运营性陆地成像仪(Operational Land Imager：OLI)采集的影像，是目前Landsat系列卫星中最新的在轨运行多光谱遥感成像仪。利用Landsat遥感影像提取水体方法较多，其中运算复杂度较低而效率较高的方法发展多种水体指数然后阈值化操作提取水体。因此，将Landsat OLI遥感影像作为水体模型建构和水体图像绘制的基础，能够获得更加精确的水体模型，特别是在区域性水体的建模过程中能够体现其显著的优势。

发明内容

本发明的目的在于克服现有水体指数模型的基础上，在其基础上提出一种基于Landsat OLI遥感影像的普适性水体指数模型。该模型充分考虑到地理空间位置变化和季节变化对地物反射率的影响，采用了图像熵方法自动调节水体指数阈值，减少人为干预，提高水体提取的精度、效率和稳定性。

本发明的技术解决方案：一种基于Landsat OLI影像水体制图法，其具体步骤如下：

1)普适性水体指数模型构建；定义普适性水体指数模型UWI的指数方程为：

其中

分别对应遥感影像选取的波段的地表反射值，Mean(b_i+b_j)表示对于选取波段地表反射值的均值和；b_i选取Landsat OLI影像的波段七：band7， SWIR2：2.100-2.300μm或波段六：band6，SWIR1：1.560-1.660μm，b_j选取Landsat OLI影像的波段四：band4，Red：0.630-0.680μm或波段三：band3，Green：0.525- 0.600μm或波段二：band2，Blue：0.450-0.515μm；

公式(1)可以组合为六种指数，分别为：

2)最佳水体指数阈值自动选取；针对遥感影像存在的时空差异性，采用图像熵的图像分割方法，对步骤1)中水体指数模型运算后不同时空的影像进行阈值的自动选取，通过迭代选取适合该影像的最佳阈值，实现较高精度水体的提取，具体操作为将图像分割为前景为水体，逻辑值1，背景为非水体，逻辑值为0，其具体步骤如下：

(1)将图像划分为k个灰度级别，给定一个初始阈值T₀，0≤T₀≤k-1；

(2)用阈值T₀将图像分成两个区域，即前景区域R₁和背景区域R₀，小于等于T₀为背景，大于T₀为前景，其每个区域中各灰度级别的概率密度可表示为：

其中p(0)，p(1)……p(k-1)，p(T₀)表示各灰度等级的概率，而p₀(T₀)和p₁(T₀) 分别表示T₀阈值对应像元的灰度级所在背景色和前景色出现概率；

(3)计算背景中各灰度级别像素的累计概率：

则前景各灰度级别的像素的累计为：

P₁(T₀)＝1-P(T₀)； (6)

(4)计算背景和前景各对应的图像熵：

(5)当阈值为T₀时，计算图像总熵：

H(T₀)＝H₀(T₀)+H₁(T₀)； (8)

(6)随着阈值变化，分别计算各阈值时对应的图像总熵；

(7)比较各阈值时图像总熵，求最大图像总熵为H(T_j)(0≤j≤k-1)，其中k-1为图像被分割的灰度级别数，T_j为第j个灰度级别对应的阈值，将最大总熵对应的阈值T_j作为最佳阈值；图像中灰度大于T_j的像元分割为前景，小于等于T_j的像元为背景。

3)区域性水体提取；根据不同地理空间区域和不同季节的地物反射率存在的差异，通过步骤2)中得到的图像熵自动选取最佳阈值T_j，对不同区域不同季节的影像自动调节阈值，提取水体的条件为：

4)水体提取结果验证；

(1)验证设计

采用2.5m高分辨率Google Earth遥感影像对提取的水体进行验证：在 LandsatOLI原始影像区域中随机选取300个点，将这些随机生成点的坐标导入对应年份的GoogleEarth影像中，通过目视解译判别对于该点的真实地物类型为水体或非水体，并进行计数统计。

(2)精度分析

对不同时空的遥感影像，随机生成300个验证点，提取各点坐标导入到对应年份的Google Earth影像中，当验证点落在云上或阴影上，验证时将这些点作为非有效点排除统计之外；从Google Earth影像上识别所有效验证点的类别为水体或非水体，分别统计落在水体和非水体中有效验证点总数，将提取后的影像与 Google Earth影像配准后逐点验证，统计出水体提取影像中包含落在水体上的正确点数和落在非水体上的正确点数，分别计算并比较各区域水体提取精度，总体精度和Kappa系数。

5)区域性水体制图：将对象区域的维度变化和季节变化对水体反射率的影响作为重要考虑因素，进行区域水体提取及制图：提取及制图过程采用Python编程语言，逐步实现原始Landsat OLI遥感图像的读取，水体指数模型波段运算，图像熵算法实现最佳水体提取阈值选取，水体提取和图像输出；其中，水体提取的效果验证采用人工解译计数方式进行，利用ArcGIS软件对原始遥感影像自动生成300个验证点，将各点的坐标导入到GoogleEarth中，将300个点地物类别在Google Earth中通过人工解译并计数各地物类别的点数，将点的坐标导入到提取的水体图像中，识别提取的水体中点的个数，并建立混淆矩阵，计算提取精度和Kappa系数。

本发明的优点：

(1)提出一种普适性区域水体指数模型，该模型可以根据选择的波段不同，组合六种指数模型，且每种指数模型均可以有效进行区域性水体提取；

(2)有效解决了地理空间变化和季节变化引起水体反射率的时空差异，解决了采用同一阈值对不同影像提取水体时存在精度差异的问题；

(3)采用了图像熵算法自动识别指数阈值，克服了人工选择的带来的不确定性和主观性；

(4)提出了区域水体制图的思想，将地理空间和季节变化对水体反射率的影响因素纳入水体提取及制图策略中，采用自动指数阈值调节方案，克服不同地理空间和不同季节的时空差异导致水体提取及制图精度差异。

附图说明

附图1是本发明基于Landsat OLI影像水体制图法流程图。

附图2是实施例中五个地区的水体的Google Earth遥感影像图。

附图3是采用本发明技术方案获得的哈尔滨地区水体制图。

附图4是采用本发明技术方案获得的云南大理地区水体制图。

具体实施方式

下面结合附图进一步说明本发明的技术方案。

如图1所示的一种基于Landsat OLI影像水体制图法，其具体步骤如下：

1)普适性水体指数模型(UWI)构建：

其中

分别对应遥感影像选取的波段的地表反射值，Mean(b_i+b_j) 表示对于选取波段地表反射值的均值和，本发明中b_i的选取为Landsat OLI影像的波段为波段七(band7，SWIR2：2.100-2.300μm)或波段六(band6， SWIR1：1.560-1.660μm)，b_j可以选取为波段四(band4，Red：0.630- 0.680μm)、波段三(band3，Green：0.525-0.600μm)、波段二(band2，Blue： 0.450-0.515μm)，该模型可以组合为六种指数，分别为：

2)最佳水体指数阈值自动选取

本发明针对遥感影像存在的时空差异性导致采用同一水体指数提取精度不一致问题，采用了图像熵的图像分割方法，对上述水体模型运算后不同时空的影像进行阈值的自动选取，通过迭代选取适合该影像的最佳阈值，实现较高精度水体的提取。本发明采用图像熵进行图像分割的目标是：将图像分割为前景为水体，逻辑值1，背景为非水体，逻辑值为0。其具体步骤如下：

(1)将图像划分为k个灰度级别，给定一个初始阈值T₀，0≤T₀≤k- 1；

(2)用阈值T₀将图像分成两个区域，即前景区域R₁和背景区域R₀，小于等于T₀为背景，大于T₀为前景，其每个区域中各灰度级别的概率密度可表示为：

其中p(0)，p(1)……p(k-1)，p(T₀)表示各灰度等级的概率，而p₀(T₀)和p₁(T₀) 分别表示T₀阈值对应像元的灰度级所在背景色和前景色出现概率；

(3)计算背景中各灰度级别像素的累计概率：

则前景各灰度级别的像素的累计为：

P₁(T₀)＝1-P(T₀)，

(4)计算背景和前景各对应的图像熵：

(5)当阈值为T₀时，计算图像总熵：H(T₀)＝H₀(T₀)+H₁(T₀)；

(6)随着阈值变化，分别计算各阈值时对应的图像总熵；

(7)比较各阈值时图像总熵，求最大图像总熵为H(T_j)(0≤j≤k- 1)，其中k-1为图像被分割的灰度级别数，T_j为第j个灰度级别对应的阈值，将最大总熵对应的阈值T_j作为最佳阈值；图像中灰度大于T_j的像元分割为前景，小于等于T_j的像元为背景。

3)区域性水体提取

不同地理空间区域和不同季节的地物反射率存在一定的差异，本发明通过图像熵自动选取最佳阈值T_j，对不同区域不同季节的影像自动调节阈值，提取水体的条件为：

4)水体提取结果验证

(1)验证设计：

本发明采用2.5m高分辨率Google Earth遥感影像对提取的水体进行验证。在Landsat OLI原始影像区域中随机选取300个点，将这些随机生成点的坐标导入对应年份的Google Earth影像中，通过目视解译判别对于该点的真实地物类型(水体或非水体)，并进行计数统计。

在地理空间位置上选择了中国五个不同地理空间维度的具有代表性区域，并对这五个区域的不同季影像采用本发明水体指数模型和自动选取指数阈值提取水体进行测试验证。这五个区域分别为：哈尔滨地区(东经：126.428°～ 129.453°，北纬：44.952°～47.090°)，太湖地区(东经：119.160°～ 121.583°，北纬：30.680°～32.790°)，鄱阳湖地区(东经：115.355°～ 117.703°，北纬：27.809°～29.916°)，青海湖地区(东经：99.055°～101.664°，北纬：36.408°～38.525°)，大理地区(东经：99.222°～ 101.494°，北纬：24.935°～27.038°)，测试遥感影像由美国USGS网站提供(地面反射率数据)。测试后提取的水体如图2所示。

(2)精度分析：

对不同时空的遥感影像，随机生成300个验证点，提取各点坐标导入到对应年份的Google Earth影像中，但Landsat OLI影像中有云和阴影的影响，当验证点落在云上或阴影上，验证时将这些点作为非有效点排除统计之外，如表2 所示。从Google Earth影像上识别所有效验证点的类别(水体或非水体)，分别统计落在水体和非水体中有效验证点总数，将提取后的影像与Google Earth影像配准后逐点验证，统计出水体提取影像中包含落在水体上点的个数(正确点数)和落在非水体上的正确点数，分别计算各区域水体提取精度，总体精度和 Kappa系数，如表3所示。从提取精度上，五个地区的都达到了96.6％以上，总体精度达到了97％以上，平均精度为97.79％。上述六种指数在相同影像中提取水体效果也存在一定的差异性。经验证，UWI_s1b有利于鉴别城市高反射人工建筑和裸土，UWI_s2r有利于鉴别植被、裸土，UWI_s1r、UWI_s1g、UWI_s1b比 UWI_s2r、UWI_s2g、UWI_s2b指数更有效鉴别城市中高反射人工建筑，在水体提取精度上，UWI_s1r、UWI_s2r指数略高于其他指数提取效果，提取精度高达98％。

表2验证点统计情况表

表3提取精度计算结果表

5)区域性水体制图

区域水体提取及制图中，将对象区域的维度变化和季节变化对水体反射率的影响作为重要考虑因素。提取及制图过程采用Python编程语言，逐步实现原始Landsat OLI遥感图像的读取，水体指数模型波段运算，图像熵算法实现最佳水体提取阈值选取，水体提取和图像输出；水体提取的效果验证采用人工解译计数方式进行，利用ArcGIS软件对原始遥感影像自动生成300个验证点，将各点的坐标导入到Google Earth中，将300个点地物类别在Google Earth中通过人工解译并计数各地物类别的点数，将点的坐标导入到提取的水体图像中，识别提取的水体中点的个数，并建立混淆矩阵，计算提取精度和Kappa系数。以哈尔滨区域和云南大理地区为例进行水体制图，其结果如图3、图4所示。

本发明设计了一种基于Landsat OLI遥感影像的普适性区域水体指数模型 UWI，通过设计该模型的算式和选取的提取波段组合，并采用图像熵分割方法自动提取最佳水体指数阈值，根据水体模型提取后的影像建立灰度多级级别分割，通过图像总熵的最大值确定最佳指数阈值，克服不同地区各季节的水体指数模型提取的时空差异。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明 构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (1)

1.一种基于Landsat OLI影像水体制图法，其特征是具体包括如下步骤：

1)普适性水体指数模型构建；

2)最佳水体指数阈值自动选取；

3)区域性水体提取；

4)水体提取结果验证；

5)区域性水体制图；

所述步骤1)具体包括如下步骤：

定义普适性水体指数模型UWI的指数方程为：

其中

分别对应遥感影像选取的波段的地表反射值，Mean(b_i+b_j)表示对于选取波段地表反射值的均值和；b_i选取Landsat OLI影像的波段七：band7，SWIR2：2.100–2.300μm或波段六：band6，SWIR1：1.560–1.660μm，b_j选取Landsat OLI影像的波段四：band4，Red：0.630–0.680μm或波段三：band3，Green：0.525–0.600μm或波段二：band2，Blue：0.450–0.515μm；

公式(1)可以组合为六种指数，分别为：

所述步骤2)具体包括如下步骤：

针对遥感影像存在的时空差异性，采用图像熵的图像分割方法，对步骤1)中水体指数模型运算后不同时空的影像进行阈值的自动选取，通过迭代选取适合该影像的最佳阈值，实现较高精度水体的提取，具体操作为将图像分割为前景为水体，逻辑值1，背景为非水体，逻辑值为0，其具体步骤如下：

(1)将图像划分为k个灰度级别，给定一个初始阈值T₀，0≤T₀≤k-1；

(2)用阈值T₀将图像分成两个区域，即前景区域R₁和背景区域R₀，小于等于T₀为背景，大于T₀为前景，其每个区域中各灰度级别的概率密度可表示为：

其中p(0)，p(1)……p(k-1)，p(T₀)表示各灰度等级的概率，而p₀(T₀)和p₁(T₀)分别表示T₀阈值对应像元的灰度级所在背景色和前景色出现概率；

(3)计算背景中各灰度级别像素的累计概率：

则前景各灰度级别的像素的累计为：

P₁(T₀)＝1-P(T₀)； (6)

(4)计算背景和前景各对应的图像熵：

(5)当阈值为T₀时，计算图像总熵：

H(T₀)＝H₀(T₀)+H₁(T₀)； (8)

(6)随着阈值变化，分别计算各阈值时对应的图像总熵；

(7)比较各阈值时图像总熵，求最大图像总熵为H(T_j)，0≤j≤k-1，其中k-1为图像被分割的灰度级别数，T_j为第j个灰度级别对应的阈值，将最大总熵对应的阈值T_j作为最佳阈值；图像中灰度大于T_j的像元分割为前景，小于等于T_j的像元为背景；

所述步骤3)具体包括如下步骤：

根据不同地理空间区域和不同季节的地物反射率存在的差异，通过步骤2)图像熵算法自动选取得到的最佳阈值T_j，对不同区域不同季节的影像自动调节阈值，提取水体的条件为：

所述步骤4)具体包括如下步骤：

(1)验证设计

采用2.5m高分辨率Google Earth遥感影像对提取的水体进行验证：在Landsat OLI原始影像区域中随机选取300个点，将这些随机生成点的坐标导入对应年份的Google Earth影像中，通过目视解译判别对于该点的真实地物类型为水体或非水体，并进行计数统计；

(2)精度分析

对不同时空的遥感影像，随机生成300个验证点，提取各点坐标导入到对应年份的Google Earth影像中，当验证点落在云上或阴影上，验证时将这些点作为非有效点排除统计之外；从Google Earth影像上识别所有效验证点的类别为水体或非水体，分别统计落在水体和非水体中有效验证点总数，将提取后的影像与Google Earth影像配准后逐点验证，统计出水体提取影像中包含落在水体上的正确点数和落在非水体上的正确点数，分别计算并比较各区域水体提取精度，总体精度和Kappa系数；

所述步骤5)具体包括如下步骤：

将对象区域的地理空间变化和季节变化对水体反射率的影响作为重要考虑因素，进行区域水体提取及制图：提取及制图过程采用Python编程语言，逐步实现原始Landsat OLI遥感图像的读取，水体指数模型波段运算，图像熵算法实现最佳水体提取阈值选取，水体提取和图像输出；其中，水体提取的效果验证采用人工解译计数方式进行，利用ArcGIS软件对原始遥感影像自动生成300个验证点，将各点的坐标导入到Google Earth中，将300个点地物类别在Google Earth中通过人工解译并计数各地物类别的点数，将点的坐标导入到提取的水体图像中，识别提取的水体中点的个数，并建立混淆矩阵，计算提取精度和Kappa系数。

Images