一种海岸带水体遥感信息全自动提取方法

摘要

本发明公开了一种海岸带水体遥感信息全自动提取方法，属于遥感信息全自动提取方法领域。其步骤包括遥感影像分割、水体信息粗提取和水体信息精提取三个阶段，过程中进行了两次尺度转换，其一为从像元到对象的转换，属于自下向上的尺度转换，其二为从全域到局部的转换，属于自上向下的尺度转换。相较现有技术，本发明实现了一种在尺度转换框架中融入地学知识和数据挖掘相结合的方法，整个方法具有零样本、零参数的特性，完全自动化运行。本发明能够适应各海域多类型海岸带环境，具有较好的稳定性，其提取结果精度较高，对细节信息的提取完整性和连续性均优于经典方法，能够直接应用于国家各级基础地理信息数据库遥感专题信息的维护与更新。

权利要求

1.一种海岸带水体遥感信息全自动提取方法，包括以下步骤：

步骤1：通过卫星遥感技术手段，选择美国陆地探测卫星Landsat 7所装备的Enhanced Thematic Mapper Plus设备获取描述地表物体对太阳光照电磁波的反射和其自身发射的电磁波的遥感影像，结合水体在遥感影像中具有特定的光谱显示特征，采用基于非参数核密度估计的Mean Shift聚类算法进行空间滤波，减轻噪声对目标和背景信息的干扰：具体对特征空间中的所有数据点，依次计算其在核窗口中的均值偏移向量，并沿均值偏移向量的方向移动该点到新的位置，然后以此为新的起点，继续计算均值偏移向量并移动，直到满足迭代条件结束运算过程；

步骤2：在Mean Shift算法滤波的基础上，采用种子生长算法进行区域标号，将具有相似性质的像素进行合并构成一个完整区域，完成遥感影像分割，得到由一系列区域对象构成的分割图像，从而实现从遥感影像像元到区域对象的转换；

步骤3：对Landsat ETM+遥感影像分割所形成的各区域对象计算水体判别依据指数，具体应用影像所拥有的TM1～TM6共6个波段的波段值，按照公式①~⑦分别计算归一化水体指数NDWI、改进归一化水体指数MNDWI、周氏归一化水体指数ZNDWI、缨帽变换湿度分量TCW四个水体指数，归一化植被指数NDVI、土壤调整植被指数SAVI两个植被指数，以及归一化建筑指数NDBI；

NDWI＝(TM2-TM4)/(TM2+TM4) 公式①MNDWI＝(TM2-TM5)/(TM2+TM5) 公式②ZNDWI＝(TM2+TM3)-(TM4+TM5) 公式③TCW＝0.1446×TM1+0.1761×TM2+0.10022×TM3

+0.10096×TM4-0.6210×TM5-0.4186×TM6 公式④

NDVI＝(TM4-TM3)/(TM4+TM3) 公式⑤SAVI＝(TM4-TM3)/(TM4+TM3+0.5)×1.5 公式⑥NDBI＝(TM5-TM4)/(TM5+TM4) 公式⑦步骤4：通过运用公式①~⑦计算所得的7个水体判别依据指数，构建一组决策规则实现遥感影像中水体区域的判别，其中，约束性条件共有五条，见不等式方程组⑧，反映了水体遥感信息的三个特点，即水体指数大于零、水体指数大于植被指数、水体指数大于建筑指数；非约束性条件共有三条，均为水体指数的阈值判别规则，具体方法为：先对各水体指数进行直方图统计，再用OTSU算法计算各水体指数的阈值，进行阈值分割，计算各水体指数的类间方差，选取类间方差最大，即区分度最好的三个指数建立阈值判别规则；本步骤中水体判别流程为：

1）对于分割图像中的任意对象，如果其特征空间满足决策规则中的所有约束性条件，则将其判别为不确定水体，如果有任意一条约束性条件不能满足，则将其判别为否定水体；

2）对于判别出的任意不确定水体对象采用了投票策略，如果其特征空间满足多数非约束性条件，则将其判别为确定水体，否则仍判别为不确定水体，最终并将水体提取结果划分为确定区、不确定区和否定区；

方程组⑧

步骤5：对遥感影像划分格网，搜索与各不确定水体对象相交的网格单元，如果其相交的网格单元内既包含确定水体对象，又包含否定水体对象，则认为当前不确定水体对象局部特征显著，分别自动选取其邻域内的确定水体对象和否定水体对象作为样本，使用朴素贝叶斯分类器对不确定水体对象进行局部地物类别判定，明确该不确定水体对象最终为确定水体对象或否定水体对象；如果其相交的网格单元内全为确定水体对象或否定水体对象，则认为其局部特征不显著，不进行局部地物类别判定；此后，遍历遥感分割影像中所有不确定水体对象，以当前对象作为应用数据集，使用朴素贝叶斯分类器进行全域地物类别判定，确定当前不确定水体对象究竟为确定水体对象或否定水体对象；

步骤6：输出遥感影像中水体识别的提取结果图。

说明书

一种海岸带水体遥感信息全自动提取方法

技术领域

[0001] 本发明涉及一种遥感影像数据的水体信息提取方法，特别是涉及一种多层次的海岸带水体遥感信息全自动提取方法。

背景技术

[0002] 水体无论作为一种独立的环境因子，还是作为一种资源，都受到格外的重视，对各种水体进行空间识别、定位及定量计算面积、体积或模拟水体动态变化，是环境与资源遥感监测与分析的重要内容。随着遥感技术的改进和遥感应用的深入，水体提取方法不断改进，已在应用中取得了较为广泛的应用。目前，常用的水体信息提取方法主要有监督分类法和阈值分割法两大类。

[0003] 利用监督分类法提取水体信息，主要是针对多光谱图像数据，采用数据挖掘的最大似然法、决策树法和人工神经网络法等技术方法，这些方法取得了很好的效果。例如：秦其明等于2001年第20期的《地理研究》中撰文“卫星图像中不同水体类型识别研究”采用最大似然法对卫星影像进行处理，实现水体类型的识别；都金康等于2001年第5期《遥感学报》中发表论文“SPOT卫星影像的水体提取方法及分类研究”通过决策树法实现水体的提取与分类；王晶晶等在2005年《海洋技术》24期的“盐田水体遥感分类方法研究”一文中以人工神经网络法方法实现水体的监督分类。在现有的监督分类方法中，虽然样本的采集是基于专家的先验知识，但在分类过程中仅利用了这些样本自身的统计特征，而没有充分利用水体目标光谱特征和空间特征，难以达到更高的精度，更重要的是，监督分类方法需要人工采集较多的样本，不利于开展大量水体信息的自动提取。

[0004] 利用阈值分割法提取水体信息，主要是指选取水陆差异比较明显的某个波段或指数，通过阈值分割的方法实现水体信息的提取。目前，已见诸文献的水体信息提取波段或指数有1984年C.Eric等在期刊《Photogram metric Engineering and Remote Sensing》第3卷的“Application of the Tasseled Cap concept to simulated Thematic Mapper data”一文中提出的缨帽变换湿度分量（TCW）；1985年，D.B.Jupp等在《Landsat Based Interpretation of the Cairns Section of the Great BarrierReef Marine Park》一书中提出的中红外波段和近红外波段；1996年，S.K.McFeeters等在期刊《International Journal of Remote Sensing》的 17 卷“The Use of Normalized Difference Water Index(NDWI)in the Delineation of Open Water Features”论文中提出的归一化水体指数（NDWI）；2005年，徐涵秋在期刊《遥感学报》9卷的“利用改进的归一化差异水体指数（MNDWI）提取水体信息的研究”论文中提出的改进归一化水体指数（MNDWI）。由于下垫面环境、大气环境和遥感影像成像情况的差异，没有一种指数具有普适性，并且采用单一的阈值进行分割，很难同时避免地物误提与细小水体漏提的现象，而采用多种指数结合的方法进行提取，尽管能够有效提高精度，但也需要人工比较分析和反复试验，同样不利于开展大量水体信息的自动提取。

[0005] 海岸带地区是自然界的特殊区域，不但受到海洋、陆地、大气等自然环境的综合影响，更受到人类活动剧烈、频繁、持久的作用。与内陆地区相比，海岸带地区的水体具有差异性大、噪声强的特点，一方面，除海水外，海岸带地区往往还有河网、湿地、潮沟、坑塘水面、养殖水面、盐田等不同类型水体分布，其本身的光谱特性比较复杂，另一方面，海岸带地区土地类型多样性较强，景观格局破碎度较高，各种背景元素交错分布，对水体遥感信息具有较大的干扰。简言之，海岸带的复杂性对该区域的水体遥感信息自动提取方法提出了更高的要求。

[0006] 针对以上问题，已有学者开展了进一步的研究。2008年，吴文渊等在学术期刊《科技通报》第24卷的“基于Landsat ETM+影像的水体信息综合提取方法”一文在先利用谱间关系和缨帽变换湿度分量建立提取水体的最优空间掩膜，再使用最大似然法去除最优空间掩膜中其他的地物信息，有效提高了水体提取的相对精度。2009年，骆剑承等在期刊《遥感学报》第13卷的“分步迭代的多光谱遥感水体信息高精度自动提取”一文中提出了“全域-局部”的分步迭代空间尺度转换机制，将全域分割、全域分类、局部分割与分类等计算过程有机地结合起来，分阶段地融合了水体信息提取所需的不同层次知识，并建立迭代算法实现了水体最佳边缘的逐步逼近，获得了高精度的水体信息提取。上述研究分别从地学知识和数据挖掘相结合及尺度转换的角度取得了突破性的进展。基于监督分类法和阈值分割法的水体遥感信息提取方法均较为经典，已在科学研究和生产实践中得到了广泛地应用，但仍存在一些共性问题。第一，两种方法均需要人工干预，难以实现水体遥感信息的全自动提取；第二，两种方法均忽略了地理学研究中的一个核心问题，即尺度问题，在不同的尺度上研究地理现象，其得出的结果往往相去甚远，这两种经典方法都是在较大的尺度上进行分析，即将整个遥感影像作为一个整体，建立统一的水体信息提取模型，在这种情况下，遥感影像统计数据反映的是研究区内不同类型不同水域光谱特征的差异性，而细节的水体信息遭到掩盖，提取的结果难免趋于粗放。

发明内容

[0007] 1.发明要解决的技术问题

[0008] 本发明的目的是提供一种海岸带水体遥感信息全自动提取方法，结合海岸带水体的特点，以尺度转换方法为框架，通过研究海岸带水体信息在不同层次水平的变化，开展多层次的海岸带水体信息提取，即先按总体结构特征建立分类标识，再对每一层次设计相应的分类决策规则，对影像中的单元进行逐层的判别分类，在尺度转换框架中融入地学知识和数据挖掘相结合的方法，通过地学知识分析在全域范围内进行粗提取，再通过数据挖掘分析在局部范围内进行精提取，以实现海岸带水体遥感信息全自动提取。

[0009] 2.本发明的技术方案如下：

[0010] 原理：本发明在对经典的水体遥感信息提取方法分析的基础上，针对海岸带水体的特点，通过“像元-对象”及“全域-局部”的双重尺度转换，综合利用了海岸带水体遥感信息的光谱特征、空间特征和统计特征，有机结合了地学知识与数据挖掘方法，实现了海岸带水体遥感信息的全自动提取。

[0011] 一种海岸带水体遥感信息全自动提取方法，包括遥感影像分割、水体信息粗提取和水体信息精提取三个阶段。过程中进行了两次尺度转换，其一为从像元到对象的转换，属于自下向上的尺度转换；其二为从全域到局部的转换，属于自上向下的尺度转换。具体技术方案为：

[0012] 一种海岸带水体遥感信息全自动提取方法，其步骤为：

[0013] 步骤1：通过卫星遥感技术手段，选择美国陆地探测卫星Landsat 7所装备的Enhanced Thematic Mapper Plus(ETM+)设备获取描述地表物体对太阳光照等电磁波的反射和其自身发射的电磁波的遥感影像，结合水体在遥感影像中具有特定的光谱显示特征，采用基于非参数核密度估计的Mean Shift聚类算法进行空间滤波，减轻噪声对目标和背景信息的干扰：具体对特征空间中的所有数据点，依次计算其在核窗口中的均值偏移向量，并沿均值偏移向量的方向移动该点到新的位置，然后以此为新的起点，继续计算均值偏移向量并移动，直到满足迭代条件结束运算过程，示意图见附图2；

[0014] 步骤2：在Mean Shift算法滤波的基础上，采用种子生长算法进行区域标号，将具有相似性质的像素进行合并构成一个完整区域，完成遥感影像分割，得到由一系列区域对象构成的分割图像，从而实现从遥感影像像元到区域对象的转换；

[0015] 步骤3：对Landsat ETM+遥感影像分割所形成的各区域对象计算水体判别依据指数，具体应用影像所拥有的TM1~TM6共6个波段的波段值，按照公式①~⑦分别计算归一化水体指数（NDWI）、改进归一化水体指数（MNDWI）、周氏归一化水体指数（ZNDWI）、缨帽变换湿度分量（TCW）四个水体指数，归一化植被指数（NDVI）、土壤调整植被指数（SAVI）两个植被指数，以及归一化建筑指数（NDBI）；

[0016] NDWI＝(TM2-TM4)/(TM2+TM4) 公式①

[0017] MNDWI＝(TM2-TM5)/(TM2+TM5) 公式②

[0018] ZNDWI＝(TM2+TM3)-(TM4+TM5) 公式③

[0019] TCW＝0.1446×TM1+0.1761×TM2+0.10022×TM3

[0020] +0.10096×TM4-0.6210×TM5-0.4186×TM6 公式④

[0021] NDVI＝(TM4-TM3)/(TM4+TM3) 公式⑤

[0022] SAVI＝(TM4-TM3)/(TM4+TM3+0.5)×1.5 公式⑥

[0023] NDBI＝(TM5-TM4)/(TM5+TM4) 公式⑦

[0024] 步骤4：通过运用公式①~⑦计算所得的7个水体判别依据指数，构建一组决策规则实现遥感影像中水体区域的判别，其中，约束性条件共有五条（见不等式方程组⑧），反映了水体遥感信息的三个特点，即水体指数大于零、水体指数大于植被指数、水体指数大于建筑指数；非约束性条件共有三条，均为水体指数的阈值判别规则，具体方法为：先对各水体指数进行直方图统计，再用OTSU算法计算各水体指数的阈值，进行阈值分割，计算各水体指数的类间方差，选取类间方差最大，即区分度最好的三个指数建立阈值判别规则；本步骤中水体判别流程为：

[0025] 1）对于分割图像中的任意对象，如果其特征空间满足决策规则中的所有约束性条件，则将其判别为不确定水体，如果有任意一条约束性条件不能满足，则将其判别为否定水体；

[0026] 2）对于判别出的任意不确定水体对象采用了投票策略，如果其特征空间满足多数非约束性条件，则将其判别为确定水体，否则仍判别为不确定水体，最终并将水体提取结果划分为确定区、不确定区和否定区；

[0027] （方程组⑧）

[0028] 步骤5：对遥感影像划分格网，搜索与各不确定水体对象相交的网格单元，如果其相交的网格单元内既包含确定水体对象，又包含否定水体对象，则认为当前不确定水体对象局部特征显著，分别自动选取其邻域内的确定水体对象和否定水体对象作为样本，使用朴素贝叶斯分类器（Naive Bayes Classifier）对不确定水体对象进行局部地物类别判定，明确该不确定水体对象最终为确定水体对象或否定水体对象；如果其相交的网格单元内全为确定水体对象或否定水体对象，则认为其局部特征不显著，不进行局部地物类别判定；此后，遍历遥感分割影像中所有不确定水体对象，以当前对象作为应用数据集，使用朴素贝叶斯分类器进行全域地物类别判定，确定当前不确定水体对象究竟为确定水体对象或否定水体对象；

[0029] 步骤6：输出遥感影像中水体识别的提取结果图。

[0030] 具体的海岸带水体遥感信息全自动提取方法框架见附图1。

[0031] 3.有益效果

[0032] 相比现有技术，本发明实现了一种在尺度转换框架中融入地学知识和数据挖掘相结合的海岸带水体遥感信息全自动方法。该方法通过面向像元的全域遥感影像分割、面向对象的全域水体信息提取和面向对象的局部水体信息提取三个过程实现了海岸带水体遥感信息的精确提取，整个方法具有零样本、零参数的特性，完全自动化运行。具体有益效果如下：

[0033] 第一，通过遥感影像分割实现从像元到对象的尺度转换。该方法对水体遥感信息提取具有三重意义：首先，使处理的对象已从像元过渡到了基元（对象）,更接近人们观测数据的思维逻辑；其次，使同一个对象中的不同像元具备相同的光谱特征，在一定程度上克服了同物异谱的问题，同时减轻了信息提取时的椒盐噪声效应；再次，面向对象分割实现了像元计算的集群化，减少了进一步分析时处理的数据量，有利于深入开展高性能遥感地学计算。

[0034] 第二，水体信息精提取中通过从全域到局部的尺度转换对水体遥感信息提取。该方法具有两重意义：首先，在局部工作区域内，由于水体与陆地的范围大致均匀，同时该单元的水体也不受周围其他水体的干扰，数据统计特性更为单纯，从而降低了局部提取模型的不确定性；其次，进一步分析的重点被聚焦在水体不确定区，使得局部水体分类模型更具针对性，同时水体确定区得以保留，从而避免了陷入局部最优，兼顾了宏观和微观的水体特征。

[0035] 综上，本发明结合海岸带水体的特点，提出一种在尺度转换的框架下融合地学知识与数据挖掘的海岸带水体遥感信息全自动提取方法。该方法通过面向像元的全域遥感影像分割、面向对象的全域水体信息提取和面向对象的局部水体信息提取三个过程实现了海岸带水体遥感信息的精确提取，整个方法具有零样本、零参数的特性，完全自动化运行。试验证明，该方法能够适应各海域多类型海岸带环境，具有较好的稳定性，其提取结果精度较高，对细节信息的提取完整性和连续性均优于经典方法。本发明的方法框架具有一定的普适价值，在实际应用中，具体的分割算法、特征空间构建方法、分类器算法等均可进行调整。

附图说明

[0036] 图1为海岸带水体遥感信息全自动提取方法框架图；

[0037] 图2为Mean Shift算法迭代过程示意图；

[0038] 图3为各实施例区域遥感影像原始图像（RGB:5,4,3）；

[0039] 图4为各实施例区域影像粗提取结果（黑色区域：否定水体；灰色区域：不确定水体；白色区域：确定水体）；

[0040] 图5为不确定水体对象邻域示意图；

[0041] 图6为实施例区域影像精提取结果（黑色区域：非水体；白色区域：水体）。

具体实施方式

[0042] 采用附图3所示的一组遥感影像作为海岸带水体遥感信息全自动提取方法的实施例遥感影像数据，并以此数据来验证海岸带水体遥感信息全自动提取方法的有效性。该组影像数据分别是我国四大海域中不同类型的海岸带区域。其中，渤海实施例区域位于河北省唐山市曹妃甸开发区附近海岸带（图3(a)），区域内盐场及养殖水面规则分布；黄海实施例区域位于江苏省盐城市大丰丹顶鹤国家级自然保护区附近海岸带（图3(b)），区域内滩涂湿地分布较多；东海实施例区域位于浙江省宁波市北仑港附近海岸带（图3(c)），区域内城市和细小水体交错复杂；南海实施例区域位于海南省海口市东寨港红树林国家级自然保护区附近海岸带（图3(d)），区域内植被茂密，大小湖泊零星分布。

[0043] 各实施例区域所使用的影像数据均为Landsat ETM+遥感影像，成像时间分别为2001年9月1日、2002年5月26日、2000年6月24日及1999年12月24日，影像大小均为500×500像元。采用标准C++编程语言在VC 6.0平台下实现了上述方法，应用该方法对四个实施例区域进行了海岸带水体遥感信息提取，并将其结果与最大似然法监督分类提取结果及NDWI阈值分割提取结果进行比较。其中，遥感影像数据的读写操作通过开源地理数据格式转换类库GDAL 1.60实现。具体实施步骤如下：

[0044] 步骤1：以GDAL为影像数据读写工具，利用GDAL.Open方法读取遥感影像，选用Landsat ETM+遥感影像所拥有的TM1~TM6共6个波段的波段值作为计算对象，构建一个长度为1500000的一维数组（500行×500列×6个波段），依次读取各波段遥感影像数据并存入该数组。

[0045] 步骤2：对于特征空间中的所有数据点（即遥感影像的所有像元）进行Mean Shift方法的空间滤波：依次计算其在核窗口中的均值偏移向量，并沿均值偏移向量的方向移动该点到新的位置，然后以此为新的起点，继续计算均值偏移向量并移动，直到满足迭代收敛条件后结束迭代运算过程。在空间滤波的基础上采用种子生长算法将具有相似性质的像素集合起来构成区域，完成遥感影像分割，得到由一系列对象构成的遥感影像分割结果。

[0046] 步骤3：对影像分割结果TM1~TM6共6个波段的数值进行7个水体判别依据指数计算，具体计算公式见本发明的技术方案中公式①~⑦，分别计算出归一化水体指数（NDWI）、改进归一化水体指数（MNDWI）、周氏归一化水体指数（ZNDWI）、缨帽变换湿度分量（TCW）、归一化植被指数（NDVI）、土壤调整植被指数（SAVI）两个植被指数以及归一化建筑指数（NDBI）。

[0047] 步骤4：根据所计算出的7个水体判别依据指数，构建一组决策规则实现遥感影像中水体区域的判别，其中，有约束性条件共有五条，为水体指数大于零、水体指数大于植被指数、水体指数大于建筑指数，具体见本发明的技术方案中公式⑧；非约束性条件共有三条，均为水体指数的阈值判别规则，具体方法为：先对各水体指数进行直方图统计，再用OTSU算法计算各水体指数的阈值，进行阈值分割，计算各水体指数的类间方差，选取类间方差最大，即区分度最好的三个指数建立阈值判别规则。具体操作，对于分割图像中的任意对象，如果其特征空间满足决策规则中的所有约束性条件，则将其判别为不确定水体，如果有任意一条约束性条件不能满足，则将其判别为否定水体；对于判别出的任意不确定水体对象采用了投票策略，如果其特征空间满足多数非约束性条件，则将其判别为确定水体，否则仍判别为不确定水体，最终并将水体提取结果划分为确定区、不确定区和否定区。此时完成面向对象的全域水体信息的粗提取结果，包括水体的确定区、不确定区和否定区，结果见图4。

[0048] OTSU算法核心代码如下所示：

[0049]

[0050]

[0051] 步骤5：构建网格单元的映射关系。将整个遥感影像划分为1km×1km的格网，与对象相交的网格单元的集合即为对象的邻域，如图5所示，不确定水体对象O的邻域为网格单元A、B、C和D，此种情况为不确定水体对象O的邻域内，既包含确定水体对象（C、B、D），又包含否定水体对象（A、B、C），则认为当前不确定水体对象局部特征显著，实现从全域到局部的转换，分别自动选取其邻域内的确定水体对象和否定水体对象作为样本，使用朴素贝叶斯分类器（Naive Bayes Classifier）对不确定水体对象进行局部分类，明确该不确定水体对象最终为确定水体对象或否定水体对象。如果其相交的网格单元内全为确定水体对象或否定水体对象，则认为其局部特征不显著，不进行局部分类。完成所有不确定水体对象的局部分类后，遍历遥感分割影像中所有不确定水体对象，以当前对象作为应用数据集，使用朴素贝叶斯分类器进行全域分类，确定当前不确定水体对象究竟为确定水体对象或否定水体对象，完成影像分割结果下的水体信息提取。

[0052] 该步骤核心代码如下：

[0053]

[0054]

[0055]

[0056]

[0057]

[0058] 步骤6：完成海岸带水体遥感信息全自动提取，并对分类结果进行图像输出，结果如图6所示，图中白色为水体区域，黑色为非水体区域。

[0059] 将其结果与最大似然法监督分类提取结果及NDWI阈值分割提取结果进行比较，提取精度结果如表1所示。

[0060] 表1精度比较结果

[0061]

[0062] 测试方法为，从各试验区分别随机抽取1000个样本点，对提取结果进行精度评价，结果表明，运用本发明进行海岸带水体遥感信息自动提取，提取结果的精度明显高于最大似然法监督分类及NDWI阈值分割。经分析发现，本发明未能正确提取的部分大多是光谱混合现象较为严重的像元。