CN111007039B - 中低分辨率遥感图像亚像元级水体自动提取方法和系统 - Google Patents

Abstract

本发明公开了一种中低分辨率遥感图像亚像元级水体自动提取方法和系统，该方法包括：对于给定的中低分辨率遥感图像进行预处理，生成地表反射率图像；对地表反射率图像提取归一化差异水体指数并进行拉伸处理，构建拉伸水体指数图像直方图；使用高斯混合模型描述拉伸水体指数图像直方图，并估算得到水体与陆地参数；根据水体与陆地参数，确定水陆分离阈值；根据水体与陆地参数、水陆分离阈值，提取得到水陆混合像元；对水陆混合像元进行水体丰度解算，得到亚像元级水体提取结果。本发明能够快速、精确地提取水陆分离阈值，且不受人为因素影响；同时对于细小水体在中低分辨率遥感图像上漏提取或过提取问题有了较好的改进。

Description

中低分辨率遥感图像亚像元级水体自动提取方法和系统

技术领域

本发明属于遥感图像技术领域，尤其涉及一种中低分辨率遥感图像亚像元级水体自动提取方法和系统。

背景技术

水体提取对于水资源调查、洪涝灾害监测预警，以及水体环境生态评估等有着重要作用。利用卫星遥感图像数据快速、准确地提取水体信息，成为水资源调查和宏观监测的重要手段。

目前，用于水体提取的遥感图像数据源和水体提取方法较多，如基于不同空间分辨率遥感图像的监督/非监督分类法、波段比值法等。由于高分辨率图像不适合大面积水体提取，因此中低分辨率图像成为水体提取的主要数据源，然而中低分辨率图像则存在混合像元问题。

针对中低分辨率遥感图像，采用水体指数法的水体提取主要表现为以下四个方面的特点：(1)通常采用阈值分割方法，而阈值的选取受图像类型、获取时间、地理位置等因素影响较大，且人工参与度较高，不利于自动化水体提取工作开展；(2)基于像元级水体提取居多，而涉及亚像元级的较少，不利于细小水体的精确提取；(3)对水陆边界简单的区域提取效果较好，但对边界复杂区域，如水产养殖区等水体提取鲜有涉及；(4)水体指数直接使用DN值或表观反射率计算，而这将直接影响水体提取精度。

在这样的形势下，传统的基于水体指数方法经过人工参与确定阈值提取水体的方法出现了新的挑战，而且对于中低分辨率图像中细小水体目标极易出现误提取现象，同时难以适应水陆边界复杂情况下的水体提取。

发明内容

本发明的技术解决问题：克服现有技术的不足，提供一种中低分辨率遥感图像亚像元级水体自动提取方法和系统，能够快速、精确地提取水陆分离阈值，且不受人为因素影响，有效解决水体提取中水陆分离阈值选择难题，同时对于细小水体在中低分辨率遥感图像上漏提取或过提取问题有了较好的改进。

为了解决上述技术问题，本发明公开了一种中低分辨率遥感图像亚像元级水体自动提取方法，包括：

对于给定的中低分辨率遥感图像进行辐射定标和大气校正，生成地表反射率图像；

对地表反射率图像进行波段运算，提取水体指数并进行拉伸处理，并基于拉伸图像构建拉伸水体指数图像直方图；

使用高斯混合模型描述构建的拉伸水体指数图像直方图，并运用EM算法估算得到水体与陆地参数；

根据水体与陆地参数，求取水体与陆地的高斯密度函数交集，并将求取的高斯密度函数交集确定为水陆分离阈值；

根据水体与陆地参数、以及水陆分离阈值，提取得到水陆混合像元；

对提取得到的水陆混合像元进行水体丰度解算，得到亚像元级水体提取结果。

在上述中低分辨率遥感图像亚像元级水体自动提取方法中，对于给定的中低分辨率遥感图像进行辐射定标和大气校正，生成地表反射率图像，包括：

根据传感器的增益和偏置信息，将给定的中低分辨率遥感图像各波段观测的DN值转换为传感器入瞳处的光谱辐亮度；

根据图像获取时间和太阳天顶角，将光谱辐亮度转换为大气顶层表观反射率，完成图像的辐射定标；

使用大气校正模型，将光谱辐亮度或大气顶层表观反射率转换为地表反射率，完成图像的大气校正，并最终生成地表反射率图像。

在上述中低分辨率遥感图像亚像元级水体自动提取方法中，对地表反射率图像进行波段运算，提取水体指数并进行拉伸处理，并基于拉伸图像构建拉伸水体指数图像直方图，包括：

对地表反射率图像进行波段运算，使用绿光谱段和短波红外谱段生成归一化差异水体指数MNDWI；

将归一化差异水体指数MNDWI线性拉伸到[0,255]范围，并重新定义为线性差异水体指数SMNDWI；

对线性差异水体指数SMNDWI构建直方图，得到伸水体指数图像直方图。

在上述中低分辨率遥感图像亚像元级水体自动提取方法中，使用高斯混合模型描述构建的拉伸水体指数图像直方图，并运用EM算法估算得到水体与陆地参数，包括：

使用高斯混合模型描述构建的拉伸水体指数图像直方图，通过EM迭代算法求取最优化水体与陆地的均值和标准差参数。

在上述中低分辨率遥感图像亚像元级水体自动提取方法中，水陆分离阈值介于水体均值与陆地均值之间。

在上述中低分辨率遥感图像亚像元级水体自动提取方法中，根据水体与陆地参数、以及水陆分离阈值，提取得到水陆混合像元，包括：

根据水体与陆地参数、以及水陆分离阈值，采用统计样本分布理论与正态分布区间估计规律，确定混合像元在拉伸水体指数图像直方图的分布范围；

根据确定的混合像元在拉伸水体指数图像直方图的分布范围，分离出纯净水体、纯净陆地像元和水陆混合像元。

在上述中低分辨率遥感图像亚像元级水体自动提取方法中，对提取得到的水陆混合像元进行水体丰度解算，得到亚像元级水体提取结果，包括：

根据水陆混合像元在拉伸水体指数图像直方图中距离纯净水体和纯净陆地均值的位置关系估算水体丰度，得到混合像元水体丰度；

根据提取的纯净水体、纯净陆地和丰度解算后的水陆混合像元，得到亚像元级水体提取结果。

本发明还公开了一种中低分辨率遥感图像亚像元级水体自动提取系统，包括：

地表反射率图像生成模块，用于对于给定的中低分辨率遥感图像进行辐射定标和大气校正，生成地表反射率图像；

直方图生成模块，用于对地表反射率图像进行波段运算，提取水体指数并进行拉伸处理，并基于拉伸图像构建拉伸水体指数图像直方图；

估算模块，用于使用高斯混合模型描述构建的拉伸水体指数图像直方图，并运用EM算法估算得到水体与陆地参数；

阈值确定模块，用于根据水体与陆地参数，求取水体与陆地的高斯密度函数交集，并将求取的高斯密度函数交集确定为水陆分离阈值；

提取模块，用于根据水体与陆地参数、以及水陆分离阈值，提取得到水陆混合像元；

解算模块，用于对提取得到的水陆混合像元进行水体丰度解算，得到亚像元级水体提取结果。

本发明具有以下优点：

本发明能够快速、精确地提取阈值，且不受人为因素影响，有效解决水体提取中阈值选择难题，同时对于细小水体在中低分辨率遥感图像上漏提取或过提取问题有了较好的改进。

附图说明

图1是本发明实施例中一种中低分辨率遥感图像亚像元级水体自动提取方法的步骤流程图；

图2是本发明实施例中一种中分辨率水陆边界的遥感图像亚像元级水体提取实验结果示意图；

图3是本发明实施例中一种低分辨率水陆边界的遥感图像亚像元级水体提取实验结果示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明公开的实施方式作进一步详细描述。

实施例1

本发明公开一种中低分辨率遥感图像亚像元级水体自动提取方法，利用高斯混合模型表达水体指数图像直方图，进而精确表达水体和非水体的分布特征，同时根据水体与非水体参数自动确定水陆分离，其与以往方法相比，最大的不同是能够自动提取水陆分离，且水陆分离具有高可靠性；在提取的水陆分离基础上，通过分析水陆分离附近混合像元的分布规律，智能判断并确定混合像元在直方图中的分布范围，进而根据混合像元分布区间对水体丰度进行解算，该方法尤其适用于水陆边界复杂区域的水体提取；另外，所使用的水体指数是基于地表反射率图像上获得的，基于水体指数法提取水体信息，其基本思想是利用绿光波段和近红外(或短波红外)波段的反差，增强水体信息的同时抑制非水体信息，由于遥感图像所记录的DN值与目标本身辐射的能量不一致，尤其是用于计算水体指数的绿光波段和红外波段受气溶胶和水汽等大气条件影响较大，为提取精确的水体指数，有必要通过辐射定标和大气校正来消除各种失真的影响，而以往的方法多数基于原始DN值或表观反射率图像。

在本实施例中，该中低分辨率遥感图像亚像元级水体自动提取方法，包括：

步骤101，对于给定的中低分辨率遥感图像进行辐射定标和大气校正，生成地表反射率图像。

在本实施例中，可以根据传感器的增益和偏置信息，将给定的中低分辨率遥感图像各波段观测的DN值转换为传感器入瞳处的光谱辐亮度；根据图像获取时间和太阳天顶角，将光谱辐亮度转换为大气顶层表观反射率，完成图像的辐射定标；使用大气校正模型，将光谱辐亮度或大气顶层表观反射率转换为地表反射率，完成图像的大气校正，并最终生成地表反射率图像。

步骤102，对地表反射率图像进行波段运算，提取水体指数并进行拉伸处理，并基于拉伸图像构建拉伸水体指数图像直方图。

在本实施例中，步可以对地表反射率图像进行波段运算，使用绿光谱段和短波红外谱段生成归一化差异水体指数MNDWI；将归一化差异水体指数MNDWI线性拉伸到[0,255]范围，并重新定义为线性差异水体指数SMNDWI；对线性差异水体指数SMNDWI构建直方图，得到伸水体指数图像直方图。

步骤103，使用高斯混合模型描述构建的拉伸水体指数图像直方图，并运用EM算法估算得到水体与陆地参数。

在本实施例中，可以采用高斯混合模型描述构建的拉伸水体指数图像直方图，通过EM迭代算法求取最优化水体与陆地的均值和标准差参数。

步骤104，根据水体与陆地参数，求取水体与陆地的高斯密度函数交集，并将求取的高斯密度函数交集确定为水陆分离阈值。

在本实施例中，水陆分离阈值介于水体均值与陆地均值之间。

步骤105，根据水体与陆地参数、以及水陆分离阈值，提取得到水陆混合像元。

在本实施例中，可以根据水体与陆地参数、以及水陆分离阈值，采用统计样本分布理论与正态分布区间估计规律，确定混合像元在拉伸水体指数图像直方图的分布范围；根据确定的混合像元在拉伸水体指数图像直方图的分布范围，分离出纯净水体、纯净陆地像元和水陆混合像元。

步骤106，对提取得到的水陆混合像元进行水体丰度解算，得到亚像元级水体提取结果。

在本实施例中，可以根据水陆混合像元在拉伸水体指数图像直方图中距离纯净水体和纯净陆地均值的位置关系估算水体丰度，得到混合像元水体丰度；根据提取的纯净水体、纯净陆地和丰度解算后的水陆混合像元，得到亚像元级水体提取结果。

实施例2

在上述实施例的基础上，下面结合一个具体实例进行说明。基于高斯混合模型进行中低分辨率遥感图像亚像元级水体自动提取的流程如下：

1)对于给定的中低分辨率遥感图像进行辐射定标和大气校正，生成地表反射率图像。

(1.1)根据公式(1)将图像各波段的DN值转换为传感器入瞳处的光谱辐亮度：

L_λ＝DN_λ·g_λ+Lo_λ···(1)

其中，λ为图像的波段号，L_λ为波段λ的光谱辐亮度(单位为w·m^-2·sr^-1·μm^-1)，DN_λ为波段λ的DN值，g_λ和Lo_λ分别为波段λ的增益(Gains)与偏置(Biases)。

(1.2)根据公式(2)，将辐亮度转换为大气顶层表观反射率：

其中，ρ_λ为波段λ的大气顶层表观反射率，π为常量，D为日地天文单位距离，ESUN_λ为波段λ在大气顶层的太阳平均光谱辐射值，θ为太阳天顶角。

(1.3)使用大气校正模型(如FLAASH、6S等)将辐亮度或表观反射率转换为地表反射率。其中模型所需的大气参数可从MODIS Terra大气标准数据产品(包含气溶胶、水汽含量和臭氧含量等)获取，为保证大气校正后的图像质量，图像获取时间与MODIS Terra的过境时间尽可能接近。

2)对上述得到的地表反射率数据进行波段运算，提取水体指数并进行拉伸处理，并基于拉伸图像构建直方图。

(2.1)根据公式(3)，计算改进的归一化差异水体指数MNDWI(ModifiedNormalized Difference Water Index)，定义为：

其中，Green为绿光谱段，MIR为短波红外谱段。

(2.2)根据公式(4)，将MNDWI线性拉伸到[0,255]范围，定义为SMNDWI(StretchedMNDWI)。计算方法如下：

其中，max和min分别为MNDWI的最大值和最小值，

表示向下取整。

3)对上述得到的直方图使用高斯混合模型表达，并使用EM算法估算模型参数，得到水体与非水体的均值和标准差。

(3.1)高斯混合模型由两个高斯概率密度分布混合而成，分别表示陆地和水体的像元分布，其数学表达式为：

其中，Θ＝{a_k,μ_k,σ_k}是一组混合模型的参数，a_k表示第k个高斯成分权重，μ_k和σ_k表示第k个高斯成分的均值和标准差，f(x|μ_k,σ_k)表示第k个高斯成分的正态分布概率密度函数。

a_k满足：

0≤a_k≤1，(k＝1，2)···(6)

f(x|μ_k,σ_k)表达式为：

(3.2)使用期望最大化(Expectation Maximization，EM)算法，通过迭代方式求解参数。算法包括两部分：E-Step和M-Step，前一步计算参数的期望值和最大化高斯混合模型似然函数，后一步重新估计并计算参数，并返回前一步迭代计算，直到满足条件停止。

(3.3)经过EM算法迭代后提取的〈μ₁,σ₁>和〈μ₂,σ₂>值即分别为陆地和水体的均值与标准差。

4)对上述得到的水体与非水体参数，计算水体与非水体概率密度函数的交集，自动确定水陆分离阈值t。

(4.1)得到高斯混合模型参数后，通过计算水体与非水体正态分布概率密度函数的交点，获取水陆分离阈值t，即：

f(t|μ₁,σ₁)＝f(t|μ₂,σ₂)···(9)

其中，求解的t值介于μ₁和μ₂之间。

5)根据上述水陆参数和获取的水陆分离阈值t，确定纯净水体、纯净陆地和水陆混合像元。

(5.1)水陆混合像元通常分布在水陆分离阈值t附近。根据统计分析中样本分布理论与正态分布区间估计规律，样本落入不同置信区间内的概率不同，例如置信区间μ±σ、μ±1.65×σ、μ±1.96×σ、μ±2.58×σ内样本落入概率分别为68％、90％、95％、99％。因此，对于水陆边界较为明显、过渡区域较为简单的遥感图像，为确保纯净水体和纯净陆地样本像元落入各自区间内的概率在95％以上，一般可选取(μ₁+2.58×σ₁)～(μ₂-2.58×σ₂)或(μ₁+1.96×σ₁)～(μ₂-1.96×σ₂)区间为水陆混合像元。对于水陆边界模糊、过渡区域较复杂的图像，水陆混合像元的选取则需要判断t值与〈μ₁,σ₁>和<μ₂,σ₂>之间的关系。

(5.2)水陆混合像元区间的判别方法如下：

以(μ₁+2.58×σ₁)～(μ₂-2.58×σ₂)为水陆混合像元区间；

若(μ₁+2.58×σ₁)≥(μ₂-2.58×σ₂)，则选取(μ₁+1.96×σ₁)～(μ₂-1.96×σ₂)为水陆混合像元区间；

若(μ₁+1.96×σ₁)≥(μ₂-1.96×σ₂)，则选取(μ₁+1.65×σ₁)～(μ₂-1.65×σ₂)为水陆混合像元区间；

若(μ₁+1.65×σ₁)≥(μ₂-1.65×σ₂)，则选取(μ₁+σ₁)～(μ₂-σ₂)为水陆混合像元区间，以适应更大范围的水陆混合过渡区域。

6)对上述生成的水陆混合像元进行水体丰度解算，得到亚像元级水体提取结果。

(6.1)在线性混合像元分解模型中，混合像元光谱被认为是由各端元光谱根据端元在混合像元中所占比例线性加权组成。在本发明中涉及的亚像元水体提取中，同样将水陆混合像元光谱当作由水体端元与非水体端元线性组合而成。需满足的两个约束条件为：a)端元丰度之和为1；b)端元丰度取值范围为0-1。

以混合像元处在(μ₁+2.58×σ₁)～(μ₂-2.58×σ₂)为例，混合像元丰度解算公式为：

0.01[(μ₂-2.58×σ₂)-(μ₁+2.58×σ₁)]×(p_mix-μ₁-2.58×σ₁)···(10)

其中，p_mix为混合像元的SMNDWI值。

(6.2)利用步骤(5)中得到的纯净水体和纯净陆地像元，结合对混合像元水体丰度解算结果得到最终的水体提取结果。

如图2，为中分辨率(Landsat5 TM数据，空间分辨率为30米)水陆边界较简单的遥感图像亚像元级水体提取实验，其中，(2a)为原始图像(波段4)，(2b)为提取的水陆混合像元，(2c)为亚像元级水体提取结果；如图3，为低分辨率(天宫二号宽波段成像仪数据，空间分辨率为100米)水陆边界较复杂的遥感图像亚像元级水体提取实验，其中，(3a)为原始图像(波段2)，(3b)为提取的水陆混合像元，(3c)为亚像元级水体提取结果。

实施例3

在上述实施例的基础上，本发明还公开了一种中低分辨率遥感图像亚像元级水体自动提取系统，包括：地表反射率图像生成模块，用于对于给定的中低分辨率遥感图像进行辐射定标和大气校正，生成地表反射率图像；直方图生成模块，用于对地表反射率图像进行波段运算，提取水体指数并进行拉伸处理，并基于拉伸图像构建拉伸水体指数图像直方图；估算模块，用于使用高斯混合模型描述构建的拉伸水体指数图像直方图，并运用EM算法估算得到水体与陆地参数；阈值确定模块，用于根据水体与陆地参数，求取水体与陆地的高斯密度函数交集，并将求取的高斯密度函数交集确定为水陆分离阈值；提取模块，用于根据水体与陆地参数、以及水陆分离阈值，提取得到水陆混合像元；解算模块，用于对提取得到的水陆混合像元进行水体丰度解算，得到亚像元级水体提取结果。

对于系统实施例而言，由于其与方法实施例相对应，所以描述的比较简单，相关之处参见方法实施例部分的说明即可。

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

本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。

Claims (5)

1.一种中低分辨率遥感图像亚像元级水体自动提取方法，其特征在于，包括：

对于给定的中低分辨率遥感图像进行辐射定标和大气校正，生成地表反射率图像；

对地表反射率图像进行波段运算，提取水体指数并进行拉伸处理，并基于拉伸图像构建拉伸水体指数图像直方图；

使用高斯混合模型描述构建的拉伸水体指数图像直方图，并运用EM算法估算得到水体与陆地参数；

根据水体与陆地参数，求取水体与陆地的高斯密度函数交集，并将求取的高斯密度函数交集确定为水陆分离阈值；

根据水体与陆地参数、以及水陆分离阈值，提取得到水陆混合像元；

对提取得到的水陆混合像元进行水体丰度解算，得到亚像元级水体提取结果；

其中：

使用高斯混合模型描述构建的拉伸水体指数图像直方图，并运用EM算法估算得到水体与陆地参数，包括：使用高斯混合模型描述构建的拉伸水体指数图像直方图，通过EM迭代算法求取最优化水体与陆地的均值和标准差参数，具体包括：

(3.1)高斯混合模型由两个高斯概率密度分布混合而成，分别表示陆地和水体的像元分布，数学表达式为：

其中，Θ＝{a_k,μ_k,σ_k}是一组混合模型的参数，a_k表示第k个高斯成分权重，μ_k和σ_k表示第k个高斯成分的均值和标准差，f(x|μ_k,σ_k)表示第k个高斯成分的正态分布概率密度函数；

a_k满足：

0≤a_k≤1…(6)

k＝1，2；

f(x|μ_k,σ_k)表达式为：

(3.2)使用EM算法，通过迭代方式求解参数；EM算法包括两部分：E-Step和M-Step，E-Step计算参数的期望值和最大化高斯混合模型似然函数，M-Step重新估计并计算参数，并返回前一步迭代计算，直到满足条件停止；

(3.3)经过EM算法迭代后提取的<μ₁,σ₁>和<μ₂,σ₂>值即分别为陆地和水体的均值与标准差；

根据水体与陆地参数，求取水体与陆地的高斯密度函数交集，并将求取的高斯密度函数交集确定为水陆分离阈值，包括：

得到高斯混合模型参数后，通过计算水体与非水体正态分布概率密度函数的交点，获取水陆分离阈值t，即：

f(t|μ₁,σ₁)＝f(t|μ₂,σ₂)…(9)

其中，求解的t值介于μ₁和μ₂之间。

2.根据权利要求1所述的中低分辨率遥感图像亚像元级水体自动提取方法，其特征在于，对于给定的中低分辨率遥感图像进行辐射定标和大气校正，生成地表反射率图像，包括：

根据传感器的增益和偏置信息，将给定的中低分辨率遥感图像各波段观测的DN值转换为传感器入瞳处的光谱辐亮度；

根据图像获取时间和太阳天顶角，将光谱辐亮度转换为大气顶层表观反射率，完成图像的辐射定标；

使用大气校正模型，将光谱辐亮度或大气顶层表观反射率转换为地表反射率，完成图像的大气校正，并最终生成地表反射率图像。

3.根据权利要求1所述的中低分辨率遥感图像亚像元级水体自动提取方法，其特征在于，对地表反射率图像进行波段运算，提取水体指数并进行拉伸处理，并基于拉伸图像构建拉伸水体指数图像直方图，包括：

对地表反射率图像进行波段运算，使用绿光谱段和短波红外谱段生成归一化差异水体指数MNDWI；

将归一化差异水体指数MNDWI线性拉伸到[0,255]范围，并重新定义为线性差异水体指数SMNDWI；

对线性差异水体指数SMNDWI构建直方图，得到伸水体指数图像直方图。

4.根据权利要求1所述的中低分辨率遥感图像亚像元级水体自动提取方法，其特征在于，根据水体与陆地参数、以及水陆分离阈值，提取得到水陆混合像元，包括：

根据水体与陆地参数、以及水陆分离阈值，采用统计样本分布理论与正态分布区间估计规律，确定混合像元在拉伸水体指数图像直方图的分布范围；

根据确定的混合像元在拉伸水体指数图像直方图的分布范围，分离出纯净水体、纯净陆地像元和水陆混合像元。

5.根据权利要求4所述的中低分辨率遥感图像亚像元级水体自动提取方法，其特征在于，对提取得到的水陆混合像元进行水体丰度解算，得到亚像元级水体提取结果，包括：

根据水陆混合像元在拉伸水体指数图像直方图中距离纯净水体和纯净陆地均值的位置关系估算水体丰度，得到混合像元水体丰度；

根据提取的纯净水体、纯净陆地和丰度解算后的水陆混合像元，得到亚像元级水体提取结果。

Images