CN108596139A - 一种基于Gabor特征显著性的遥感影像城市区域提取方法 - Google Patents

Abstract

本发明公开了一种基于Gabor特征显著性的遥感影像城市区域提取方法，包括：S1确定Gabor滤波器参数构造二维Gabor滤波器核函数；S2利用构造好的Gabor滤波器对图像进行多方向多次卷积运算；S3多方向卷积结果综合；S4方向综合后的特征图进行规则化；S5特征图进行平滑运算；S6确定城市区域分割阈值；S7形态学综合消除噪声点和孤立建筑物并将城区连通为整体。本发明方法需人工干预，可以对不同类型传感器、不同空间分辨率、不同场景的影像，均能实现准确的城市区域自动化提取，并且提取精度高，满足海量遥感影像自动化生产的需要。

Description

一种基于Gabor特征显著性的遥感影像城市区域提取方法

技术领域

本发明属于遥感图像处理技术领域，特别涉及一种基于Gabor特征显著性的遥感影像城市区域提取方法。

背景技术

随着国内外对地观测技术的蓬勃发展，对地观测系统日趋完善，光学遥感影像的数据量迅速增加，市场化程度逐年提高，已广泛应用于城镇化研究、土地调查、防灾减灾等领域。在这种应用背景下，研究针对海量遥感影像的自动化城市区域提取方法，有着重要的意义。

遥感影像场景复杂，包含信息量巨大，其中纹理丰富、边缘明显的人工目标更容易引起人类视觉系统的注意。基于认知心理学原理的视觉显著性分析方法，常被用于遥感影像建筑物、船只、飞机等目标的识别和提取。不同学者立足于不同的应用，给出了各不相同的显著性计算方式。一般思路是：找到能够反映人工目标的底层特征作为显著性的主要测度，对特征进行提取和量化表达，然后通过特定的分割或投票策略，对视觉显著区域进行定性标记，从而判定某一个像素或超像素是否为目标。然而，由于人类活动的多样性，人工目标的形态复杂多变，且在不同尺度空间中的特征表现各不相同，与此同时，自然地表和云层在某些情况下，也会表现出丰富的纹理和明显的边缘，给目标的识别造成干扰。

发明内容

本发明的目的是针对对地观测光学遥感卫星影像，提供一种高效而精确的基于Gabor特征显著性的遥感影像城市区域提取方法。

为达到上述目的，本发明采用如下技术方案：一种基于Gabor特征显著性的遥感影像城市区域提取方法，包括如下步骤：

步骤1，确定Gabor滤波器参数构造二维Gabor滤波器核函数；

步骤2，利用构造好的Gabor滤波器对遥感影像进行多方向多次卷积运算，得到对应的单方向Gabor特征图；

步骤3，多方向Gabor特征图综合；

步骤4，方向综合后的特征图进行规则化；

步骤5，使用二维高斯滤波器对方向综合后的特征图进行平滑运算；

步骤6，确定城市区域分割阈值，根据此阈值对Gabor平滑后的特征图进行分割获得城市区域；

步骤7，形态学综合消除噪声点和孤立建筑物并将城市区域连通为整体，从而获得遥感影像城市区域。

进一步的，步骤1中二维Gabor滤波器核函数的表达式为，

其中i为虚数单位，(x,y)为像素位置坐标；λ为正弦函数的波长；σ为高斯函数标准差；γ为高斯函数纵横比；θ为Gabor核函数的方向；φ为正弦波的初始相位，x'＝xcosθ+ysinθ,y'＝-xsinθ+ycosθ；波长、标准差、纵横比、方向和初始相位的确定方式如下，

步骤1.1，根据遥感影像的空间分辨率，确定Gabor滤波器的尺度参数，包括波长λ、标准差σ、纵横比γ和初始相位φ，其中，Gsd为遥感影像空间分辨率，γ固定取值1，φ固定为0°；

步骤1.2，以22.5°为采样间隔，将[0°,180°)区间划分出8个不同的方向作为方向参数θ的取值；

步骤1.3，确定滤波器窗口大小W，进一步包括，

步骤1.3.1，根据如下公式设置初始滤波窗口W₀：

步骤1.3.2，计算当前二维Gabor滤波器最外围一圈核函数值的绝对值之和S；

步骤1.3.3，若S<10^-8，则以此时窗口的大小W为最终的采用值；否则，令W＝W+2，重复1.3.2-1.3.3，直至满足条件。

进一步的，步骤3中多方向卷积结果综合的计算公式如下，

式中，g_θ(x,y)表示像点(x,y)在θ方向上的Gabor特征值；表示对应像点综合了8个方向之后的特征值，它是一个非负值。

进一步的，步骤4中规则化的计算公式如下，

式中，表示像点(x,y)处规则化后的Gabor特征值；表示方向综合后的非负特征值；G(x,y)表示影像亮度值；f(i,j)表示二维Gabor滤波器在(i,j)处的核函数的值；W为滤波器窗口大小；即为滤波器权重之和。

进一步的，步骤5中二维高斯滤波器标准差的计算公式为，

σ_gauss＝2σ_gabor

式中，σ_gabor为步骤1.1中计算出的Gabor滤波器标准差σ。

进一步的，步骤6中将经验阈值与Otsu自动阈值分割法相结合确定城市区域分割阈值，

步骤6.1，观察分析得到经验双阈值，低阈值T_L＝0.01，高阈值T_H＝0.03；

步骤6.2，像素亮度值低于T_L的为非城区，高于T_H的为城区；在T_L与T_H之间的像素进行Otsu阈值分割，得到分割阈值T，高于阈值T的像素认为是视觉显著区域，即城区，其余为背景。

进一步的，步骤7进一步包括如下子步骤，

步骤7.1，消除面积小于S₁的目标，其中S₁＝(20σ_gabor)²，σ_gabor为步骤1.1中计算出的Gabor滤波器标准差σ；

步骤7.2，在步骤7.1的基础上执行尺度为S₂形态学膨胀运算，其中S₂＝15σ_gabor；

步骤7.3，在步骤7.2的基础上执行尺度为S₃形态学腐蚀运算，其中S₃＝10σ_gabor。

本发明无需人工干预，可以对不同类型传感器、不同空间分辨率、不同场景的影像，均能实现准确的城市区域自动化提取。

本发明具有如下优点和有益效果：

(1)不受传感器的限制，适用于不同部分辨率不同场景下影像；

(2)提取出的城区精确度高，漏检率和误检率低；

(3)可实现全自动化处理，效率高，满足海量遥感影像自动化生产的需要。

附图说明

图1为本发明实施例遥感影像城区提取流程图。

图2为本发明实施例不同波长及初相的Gabor滤波器(上)和特征图(下)。

图3为本发明实施例不同标准差及纵横比的Gabor滤波器(上)和特征图(下)。

图4为本发明实施例不同方向的Gabor滤波器(上)和特征图(下)。

具体实施方式

为了更好地理解本发明技术方案，下面将结合附图和实施例对本发明做进一步详细说明。

参见图1，本发明实施例所提供方法包含以下步骤。

步骤1，确定Gabor滤波器参数构造二维Gabor滤波器核函数。

在空间域中，一个二维Gabor函数是一个由正弦平面波调制的高斯核函数，可以看作是一个高斯函数和一个复指数函数相乘，通常的形式为：

其中i为虚数单位，(x,y)为像素位置坐标；λ为正弦函数的波长；σ为高斯函数标准差；γ为高斯函数纵横比；θ为Gabor核函数的方向；φ为正弦波的初始相位，x'＝xcosθ+ysinθ,y'＝-xsinθ+ycosθ。波长、标准差、纵横比、方向和初始相位需要根据影像信息进行确定。

本步骤的具体实施方式如下：

步骤1.1，确定Gabor滤波器的尺度参数，包括波长λ、标准差σ、纵横比γ和初相φ。初始相位φ影响的是滤波器的中心是亮条纹还是暗条纹，如图2所示。对目标提取不会造成影响，固定取值为0°即可。γ影响的是滤波器的形状是否为圆形，如图3所示，固定取值为1。λ和σ对滤波器条纹的粗细和间隔大小有极大影响，如图2和图3所示不同波长和标准差的取值对滤波器的影响不同，二者与影像空间分辨率有关，计算公式如下：

其中，Gsd为地面空间分辨率，一般Gsd取[0.5m,10m]，这是主流光学遥感影像常用的分辨率范围。

步骤1.2，以22.5°为采样间隔，将[0°,180°)区间划分出8个不同的方向作为方向参数θ的取值。θ取不同值时Gabor滤波器及其滤波后特征图如图4所示。

步骤1.3，确定滤波器窗口大小。本具体实施中策略如下：

步骤1.3.1设置初始滤波窗口W₀：

步骤1.3.2计算当前二维Gabor滤波器最外围一圈核函数值的绝对值之和S；

步骤1.3.3若S<10^-8，则以此时窗口的大小W为最终的采用值；否则，令W＝W+2，重复步骤1.3.2至1.3.3，直至满足条件。

步骤2，利用构造好的Gabor滤波器对图像进行多方向多次卷积运算。

基于步骤1中确定了参数的8个不同方向的Gabor滤波器，分别对图像进行卷积运算得到对应的单方向Gabor特征图。

步骤3，多方向卷积结果综合。将8个方向的Gabor特征图进行综合，计算公式如下：

式中，g_θ(x,y)表示像点(x,y)在θ方向上的Gabor特征值；表示对应像点综合了8个方向之后的特征值，它是一个非负值。

步骤4，特征图规则化处理。按照如下公式对方向综合后的特征图进行规则化。

式中，表示像点(x,y)处规则化后的Gabor特征值；表示方向综合后的非负特征值；G(x,y)表示影像亮度值；f(i,j)表示二维Gabor滤波器在(i,j)处的核函数的值；W为滤波器窗口大小；即为滤波器权重之和。

步骤5，特征图平滑处理。使用二维高斯滤波器对方向综合后的特征图进行平滑运算，其中二维高斯滤波器标准差的计算公式为：

σ_gauss＝2σ_gabor

式中，σ_gabor为Gabor滤波器标准差σ。

步骤6，确定城市区域分割阈值。

将经验阈值与Otsu自动阈值分割法相结合进行城市区域的提取，主要步骤包括：

步骤6.1，通过观察并分析一定数目不同传感器、不同场景的特征图，得到固定经验双阈值，低阈值T_L＝0.01，高阈值T_H＝0.03。

步骤6.2，像素亮度值低于T_L的为非城区，高于T_H的为城区。在T_L与T_H之间的像素进行Otsu阈值分割，得到分割阈值T，高于阈值T的像素认为是视觉显著区域，即城区，其余为背景。

步骤7，对上一步得到的影像进行形态学综合运算。

为了消除小面积噪声点和孤立建筑，同时使城市区域联通趋于片状整体，本步骤使用形态学运算的方式对上一步的分割结果进行综合。具体步骤包括：

步骤7.1，消除面积小于S₁的目标，其中S₁＝(20σ_gabor)²，σ_gabor为步骤1.1中计算出的Gabor滤波器标准差σ。

步骤7.2，在步骤7.1的基础上执行尺度为S₂形态学膨胀运算，其中S₂＝15σ_gabor；

步骤7.3，在步骤7.2的基础上执行尺度为S₃形态学腐蚀运算，其中S₃＝10σ_gabor。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (7)

1.一种基于Gabor特征显著性的遥感影像城市区域提取方法，其特征在于，包括如下步骤：

步骤1，确定Gabor滤波器参数构造二维Gabor滤波器核函数；

步骤2，利用构造好的Gabor滤波器对遥感影像进行多方向多次卷积运算，得到对应的单方向Gabor特征图；

步骤3，多方向Gabor特征图综合；

步骤4，方向综合后的特征图进行规则化；

步骤5，使用二维高斯滤波器对方向综合后的特征图进行平滑运算；

步骤6，确定城市区域分割阈值，根据此阈值对Gabor平滑后的特征图进行分割获得城市区域；

步骤7，形态学综合消除噪声点和孤立建筑物并将城市区域连通为整体，从而获得遥感影像城市区域。

2.如权利要求1所述的一种基于Gabor特征显著性的遥感影像城市区域提取方法，其特征在于：步骤1中二维Gabor滤波器核函数的表达式为，

其中i为虚数单位，(x,y)为像素位置坐标；λ为正弦函数的波长；σ为高斯函数标准差；γ为高斯函数纵横比；θ为Gabor核函数的方向；φ为正弦波的初始相位，x'＝xcosθ+ysinθ,y'＝-xsinθ+ycosθ；波长、标准差、纵横比、方向和初始相位的确定方式如下，

步骤1.1，根据遥感影像的空间分辨率，确定Gabor滤波器的尺度参数，包括波长λ、标准差σ、纵横比γ和初始相位φ，其中，Gsd为遥感影像空间分辨率，γ固定取值1，φ固定为0°；

步骤1.2，以22.5°为采样间隔，将[0°,180°)区间划分出8个不同的方向作为方向参数θ的取值；

步骤1.3，确定滤波器窗口大小W，进一步包括，

步骤1.3.1，根据如下公式设置初始滤波窗口W₀：

步骤1.3.2，计算当前二维Gabor滤波器最外围一圈核函数值的绝对值之和S；

步骤1.3.3，若S<10^-8，则以此时窗口的大小W为最终的采用值；否则，令W＝W+2，重复1.3.2-1.3.3，直至满足条件。

3.如权利要求2所述的一种基于Gabor特征显著性的遥感影像城市区域提取方法，其特征在于：步骤3中多方向卷积结果综合的计算公式如下，

θ∈{0°,22.5°,45°,67.5°,90°,112.5°,135°,157.5°}

式中，g_θ(x,y)表示像点(x,y)在θ方向上的Gabor特征值；表示对应像点综合了8个方向之后的特征值，它是一个非负值。

4.如权利要求3所述的一种基于Gabor特征显著性的遥感影像城市区域提取方法，其特征在于：步骤4中规则化的计算公式如下，

式中，表示像点(x,y)处规则化后的Gabor特征值；表示方向综合后的非负特征值；G(x,y)表示影像亮度值；f(i,j)表示二维Gabor滤波器在(i,j)处的核函数的值；W为滤波器窗口大小；即为滤波器权重之和。

5.如权利要求4所述的一种基于Gabor特征显著性的遥感影像城市区域提取方法，其特征在于：步骤5中二维高斯滤波器标准差的计算公式为，

σ_gauss＝2σ_gabor

式中，σ_gabor为步骤1.1中计算出的Gabor滤波器标准差σ。

6.如权利要求5所述的一种基于Gabor特征显著性的遥感影像城市区域提取方法，其特征在于：步骤6中将经验阈值与Otsu自动阈值分割法相结合确定城市区域分割阈值，

步骤6.1，观察分析得到经验双阈值，低阈值T_L＝0.01，高阈值T_H＝0.03；

步骤6.2，像素亮度值低于T_L的为非城区，高于T_H的为城区；在T_L与T_H之间的像素进行Otsu阈值分割，得到分割阈值T，高于阈值T的像素认为是视觉显著区域，即城区，其余为背景。

7.如权利要求6所述的一种基于Gabor特征显著性的遥感影像城市区域提取方法，其特征在于：步骤7进一步包括如下子步骤，

步骤7.1，消除面积小于S₁的目标，其中S₁＝(20σ_gabor)²，σ_gabor为步骤1.1中计算出的Gabor滤波器标准差σ；

步骤7.2，在步骤7.1的基础上执行尺度为S₂形态学膨胀运算，其中S₂＝15σ_gabor；

步骤7.3，在步骤7.2的基础上执行尺度为S₃形态学腐蚀运算，其中S₃＝10σ_gabor。