基于亚像元定位的遥感影像水体提取方法与流程

本发明涉及遥感影像技术领域，具体涉及一种基于亚像元定位的遥感影像水体提取方法。

背景技术：

水是地球上包括人类在内的所有生物生存必不可少的重要物质，也是工农业生产、经济发展及环境改善所必须的自然资源。随着经济的快速发展，城市化、人口的急剧增加，对水资源的需求量也急剧增加，水资源安全问题日益突出，世界面临着严峻的水资源威胁。-。如何设计一套有效的方法，对水体进行科学、准确、快速的监测与评价，实时掌握水体变化信息，继而采用对应的策略，是水资源的有效利用、合理规划以及保护的关键因素，对社会经济和生态环境的可持续发展具有重要意义。

传统获取水体信息的方法是通过实地调查的方式。由于遥感数据具有观测范围大、实时性和动态性的优势，利用遥感技术方法获取水体信息，进而对水资源进行调查和监测，现已有广泛的研究及丰富的研究成果。高精度的水体遥感监测是对地表水体进行空间识别和定位、水体面积计算，或者水体时空变化动态模拟的基础。现阶段，利用遥感技术提取水体信息主要利用计算机自动解译，许多学者提出了相应的遥感影像水体提取方法，如密度分割法、lbv变换等，其中比较常用的方法包括单波段阈值法、谱间关系法、监督分类法和水体指数法等。然而，常用的基于光谱分析的水体提取方法对小河道水体提取整体精度不高，而且这些方法与具体应用和研究区域有关，可移植性较差，依赖于专家的经验和已有的参考数据，对地表复杂多样的区域水体提取效果不太理想；需要人工确定阈值，阈值需要多次实验来确定，且容易受到人为主观影响。目前，亚像元定位的思想已经广泛应用于遥感影像分类，但针对水体特殊的地理对象，如何利用同时考虑像元间空间相关性以及像元内部亚像元间相关性的亚像元定位思想，自动地提取遥感影像亚像元尺度的水体信息的研究与应用相对较少。

技术实现要素：

针对现有技术的不足，本发明的目的在于提供一种基于亚像元定位的遥感影像水体提取方法，通过引入亚像元定位的思想，实现遥感影像亚像元尺度水体提取的目的。

为了实现上述目的，本发明采取的技术方案是：

一种基于亚像元定位的遥感影像水体提取方法，其特征在于：

步骤1：获取遥感影像素材文件并进行数据预处理，并通过最小噪声分离变换和纯净像元指数计算的过程得到建设用地、植被、裸土以及水体端元，结合线性光谱解混模型并采用最小二乘法求解每个像元中各端元所占的比例，所采用的数学模型为:式中，i＝1,2…,m，m为光谱波段数；n为端元数目；ri是波段i的反射率；fk是端元k在像元中所占的面积比例；rik是波段i上端元k的光谱反射率；eri是波段i的残差；

步骤2：利用陆地成像仪的绿波段和短波红外波段计算遥感影像的归一化水体指数，即式中，ρgreen和ρswir1分别为绿波段和短波红外波段反射率；

步骤3：在归一化水体指数图中利用最大类间方差法自动提取区分水体和非水体的阈值，即

式中，μ1和μ0分别为影像中属于水体和非水体像元的平均值，μ为整幅影像像元平均值，w1和w0分别为影像中水体和非水体的比例；

步骤4：结合步骤2的归一化水体指数和步骤3的阈值提取纯水体，归一化水体指数大于阈值的像元视为纯水体像元，得到纯水体图；

步骤5：采用步骤4得到的纯水体图像对水体比例进行优化处理，如果该像元是纯水体，那么将线性光谱解混结果的水体比例设置为1；

步骤6：通过考虑像元间的空间相关性对优化的线性光谱解混结果中的水体进行亚像元初步定位，计算每一个亚像元到周围粗像元的距离的倒数与粗像元水体比例的乘积的累加和，即式中，n为周围粗像元数目，f(pk)为粗像元pk对应的水体比例，dk为亚像元的几何中心和粗像元pk的几何中心的欧式距离，最后对每一个粗像元中的亚像元对应的累加和从大到小进行排序，根据当前粗像元的水体比例确定亚像元的数目及其空间分布；

步骤7：通过步骤6获取初步亚像元定位结果后，通过像元交换技术对像元内亚像元间水体和非水体交换，实现亚像元尺度上空间相关性达到最大从而得到水体最正确的位置分布，对每个亚像元pij，其受到邻域亚像元pk中属于水体的引力之和为，其中，nb为邻域亚像元个数，λk为亚像元间的空间相关性权重，x1k为亚像元的像元值，表达式为，式中，a0为参数，dk为亚像元pij的几何中心和pk几何中心间的欧式距离，在计算出水体所有亚像元的a1,ij后，像元交换技术以粗像元为单位进行操作，具体过程包括以下三步：

a、选择粗像元，找出所有属于水体的亚像元中最小引力值的空间位置和属于非水体的亚像元中最大引力值的空间位置，如果最小引力值小于最大引力值，则将两处的亚像元类别进行交换，否则不作任何改变；

b、对每个像元重复过程a；

c、重复步骤a和b至一定的迭代次数或者前后两次迭代的结果相差非常小为止。

相比现有技术，本发明实现了半自动遥感影像提取亚像元尺度水体的流程，引入了亚像元定位的思想，提高了水体提取精度。具体有益效果如下：

第一，本发明提高了线性光谱解混水体的提取结果。现有的线性光谱解混模型提取水体比例时容易受到端元选取与数值求解的影响，导致精度不高，特别是在河流的水体提取上，容易受到水质的影响以及数值求解影响导致水体比例不为1；而本方法通过otsu方法自动获取纯水体图像，然后再对线性光谱解混的水体比例进行处理，大大提高水体比例的提取精度。

第二，本发明通过引入亚像元定位的思想，提出一种充分考虑像元间空间相关性以及像元内部亚像元间相关性的定位方法，确定了水体在混合像元中的空间分布情况，将水体信息从像元级深入到亚像元级，获取连续的亚像元尺度水体信息。

第三，本发明的结果为亚像元尺度的水体-背景二值影像，可以便捷地实现栅格向矢量的转换，生成连续的区域水体，为更新基础地理数据库中的陆地水体包括湖泊，河流等提供更精细尺度的数据。

综上所述，本发明提高了现有陆地水体提取方法的提取效率和正确率，降低了水体提取后续处理的复杂度，提取结果可以十分方便地进行数据格式转换。实际应用表明，通过本发明能够快速、准确地提取高分辨率的水体信息，降低人工数字化采集水体信息的成本，降低基础地理信息数据库的更新成本，提高数据更新的效率。

附图说明

图1：待提取河流的遥感影像图；

图2：线性光谱解混结果图；

图3：归一化水体指数图；

图4：通过otsu阈值提取的纯水体像元结果图；

图5：亚像元尺度的水体提取结果图。

具体实施方式

下面结合附图和实施例对本发明作详细说明。

采用附图1所示遥感影像作为待提取河流影像，影像大小为401行，615列，成像时间为2015年10月18日，成像传感器为landsat8oli，空间分辨率为30米，总共9个波段。采用envi扩展工具——全约束最小二乘线性光谱解混工具实现河流比例提取的步骤；在这个基础上，结合arcgis软件，采用c#编程语言在visualstudio2012平台下实现本方法的其它6个步骤。

步骤1：利用envi软件对遥感影像进行大气校正、和几何校正等预处理；通过对预处理后的影像进行最小噪声分离变换、纯净像元指数计算以及端元选取等流程获取光谱曲线图，采用envi的全约束最小二乘线性光谱解混工具解混遥感影像，得到的结果中包括建设用地、植被、裸土和水体这四类地物比例，如附图2所示。

步骤2：利用landsat8影像波段3(绿波段)和波段6(短波红外波段)计算遥感影像的归一化水体指数(附图3)，具体核心代码如下：

步骤3：在附图3中利用otsu算法自动寻找区分纯水体与陆地的最佳阈值，具体核心代码如下：

步骤4：结合附图3的归一化水体指数，利用otsu计算的阈值提取纯水体，归一化水体指数大于阈值的像元认为是纯水体像元，如附图4所示，具体核心代码如下：

步骤5：尽管全约束最小二乘线性光谱解混方法可以较好地提取水体，但是由于端元选择及数值求解问题使得水体的提取仍存在较大误差。为了进一步改善水体的提取精度，采用归一化水体指数对水体比例进行优化处理。结合附图4中的纯水体像元对线性光谱解混结果进行修改。如果该像元是纯水体，那么将解混结果的水体比例设置为1。具体核心代码如下：

步骤6：设定亚像元空间分辨率为6米，设置搜索窗口为5×5，建立一个数组存储每一个亚像元到窗口内landsat粗像元(空间分辨率为30米)的距离的倒数与粗像元水利比例的乘积的累加和，并对这个数组存储的结果按从大到小进行排序，最后根据当前粗像元的水体比例确定亚像元的数目及其空间分布。核心代码如下：

步骤7：通过步骤6获取初步亚像元定位结果后，通过像元交换(psa)技术进行像元内亚像元间水体和非水体的交换，使得按照亚像元尺度上空间相关性达到最大的趋势得到水体最正确的位置分布。首先选择粗像元pab，计算像元pab所有属于水体的亚像元引力值以及属于非水体亚像元的引力值；然后找出属于水体的亚像元中最小引力值的空间位置和属于非水体的亚像元中最大引力值的空间位置；如果最小引力值小于最大引力值，则将这两个水体亚像元和非水体亚像元进行交换，否则，不作任何修改；对每个粗像元重复上述步骤，重复上述过程至一定的迭代次数或者前后两次迭代结果相差非常小为止，通常差值应小于0.001，这里设置迭代次数为50次。核心代码如下：

通过上述7个步骤最终得到二值化的亚像元尺度的水体提取结果(附图5)。

上列详细说明是针对本发明可行实施例的具体说明，该实施例并非用以限制本发明的专利范围，凡未脱离本发明所为的等效实施或变更，均应包含于本案的专利范围中。