一种频域高分辨率遥感图像线状特征检测方法

周立国,冯学智,肖鹏峰,马蔚纯

1.复旦大学环境科学与工程系,上海200433;2.南京大学地理与海洋科学学院,江苏南京200091

一种频域高分辨率遥感图像线状特征检测方法

周立国1,冯学智2,肖鹏峰2,马蔚纯1

1.复旦大学环境科学与工程系,上海200433;2.南京大学地理与海洋科学学院,江苏南京200091

高分辨率遥感图像中线状信息过分细节化,相应的对特征的提取带来很大干扰。引入一种基于方向和频率特征的遥感图像频域线状特征检测方法,并通过傅氏变换将图像变换到频率域。在详细分析线状特征和谱线的关系、线状特征和图像频率之间关系的基础上,由分析得到的方向和频率的参数构造Gabor滤波器进行图像线状特征的提取。高分辨率遥感图像试验结果表明,该方法能较好地提取图像的线状特征。

线状特征;频域;Gabor滤波器;特征检测;高分辨率遥感图像

1 引 言

遥感影像中线状特征的自动提取是一个经典研究项目。目前大比例尺遥感制图成为高分辨率遥感图像应用的一大特点[1-2],其核心和关键就在于体现地物形状与结构信息的线状特征提取。另外,线状特征又是道路、机场、河流、山脊线等地物目标的载体,也是各具体地物目标提取的前提和基础。因此,线状特征提取方法的探索,一方面将充实遥感图像特征提取的理论、方法体系;另一方面,将提高遥感图像利用的效率,满足大比例尺遥感制图、道路网更新、军事目标探测、环境监测、影像匹配等实际应用需求。

图像线状特征的检测方法大致可以归为空域和频域两类。空域方法如基于梯度来标识阶跃不连续性的Sober算子、Canny算子等[3-4]。而对于某些在空间域中难于处理或处理起来比较复杂的问题,可以利用傅里叶变换把用空间域表示的图像映射到频率域,利用频域频谱分析、频域滤波等方法进行特征的识别和提取[5]。将图像通过傅氏变换从空域变换到频域可使空域离散的像元信息转变为频域中对应的特征信息,使直接面向频域中的特征信息进行操作成为可能。在频域的处理方法除了常用的小波分析以外,直接在频域进行图像特征提取的研究尚不多见[6]。已有的方法包括利用图像傅里叶变换后在频谱上形成的特征峰值进行识别、提取[7-8],和源于线状信息的相位叠合处的能量最大原理,基于相位一致性进行梯度特征的提取[9-11],或利用完全覆盖频域平面的Gabor滤波器或对数 Gabor滤波器进行线状特征的提取[12-13]等。目前存在的问题主要有:①对空域特征变换到频域的频谱分析还不够清楚,对地物特征到频谱的变换机制研究还有待深入;②往往采用共性的滤波器进行信息提取,还很少有针对于频谱分析结果而设计相应滤波器进行特征提取。

针对以上问题,本文以详细的频谱分析入手,在充分研究遥感影像谱线方向与影像线状特征关系、频率和线状特征局部变化的关系等一些关键理论问题的基础上,基于频谱分析构造具有方向和频率特定参数的 Gabor滤波器进行线状特征提取。

2 方法描述

2.1 基本方法

本文以高空间分辨率卫星图像为试验数据,经过图像融合和直方图均衡化等预处理后,首先对图像进行快速傅里叶变换,通过图像变换分析和数学推导进行线状特征的频谱分析和推论证明,得到线状特征和变换到频域的谱线之间的方向关系;然后通过线状特征的短时傅里叶谱估计,得到线状信息所在的频率段范围;结合频谱的方向特征和频率分布特征设计滤波器,通过匹配Gabor滤波器实现线状特征的滤波提取,如图1所示。

图1 基于频谱分析的线状特征检测过程Fig.1 Linear feature detection based on frequency spectrum analysis

2.2 试验数据

所用的高分辨率遥感数据为 QuickBird图像,成像时间2004-11-21。QuickBird数据设置蓝(0.45μm～0.52μm)、绿(0.52μm～0.60μm)、红(0.63μm～0.69μm)和近红外(0.76μm～0.90μm)四个多光谱波段,空间分辨率为2.5 m;以及全色波段(0.45μm～0.90μm),空间分辨率为0.6 m。本文选用的是图像信息相对丰富的全色(Pan)波段,并在全色波段上截取512×512像素大小的建筑厂房图像进行试验。

3 图像特征的分析与检测

3.1 线状特征的谱线方向

频谱作为物质的能量特征之一,与波谱、能谱、重力、磁力等特征一样,都是用来区别物体属性的重要依据[14]。从傅里叶光学的观点来看,每一种图像结构都有自己特定的光学空间频谱[15],因此要进行线状特征的提取,关键是找到线状特征在频谱中对应的频谱信息,即找到线状特征转换到频域的频谱变换机制[16]。对于线状特征变换到频谱中的方向特性,因线状特征在空域表现为某一方向延伸的特性,故以竖直、水平走向的线状特征为例,分析其变换到频域的谱线方向特征[17]。对于图像处理来说,涉及的图像函数都是离散实函数,设 f(x,y)表示一幅大小为 M×N的图像,图像的离散傅里叶变换公式为

应用欧拉公式展开得

式中

式中,Re(u,v)为 F(u,v)的实部;Im(u,v)为F(u,v)的虚部。下面对含单线的 N×N图像子块进行频谱分析,通过公式推导线状特征傅里叶变换后的谱线方向,以实部Re(u,v)为例推导线状特征的图像变换频谱特征

其中第一项

设AN的第k个行向量用ξk表示,则有

根据正交性质,AN中任意两向量内积有

当图像子块含有竖直走向边缘时,设 f(x,y)表示为

a,b分别为沿边缘走向两侧的像素值,代入式(6)进行乘积运算,得到模拟线状信息的变换结果。

同理,因式(5)都是三角函数内积,故结果也都为式(10)结果的形式。为方便后续计算,一般都将傅里叶变换后的频谱做移频处理,将其移为中心对称。

由公式推导得出,傅里叶谱分布在与其空间域线状走向垂直的方向上。这里再以阶跃型和脊型两种线状特征合成图像的傅里叶变换验证公式推导结果,图2(a)、(c)为阶跃型和脊型线状特征图像,图2(b)、(d)为合成图像的傅里叶变换后将直流分量移到中心的频谱。由变换后的频谱图可以看到,当图像中沿θ方向的线状信息大量存在时,则在频率域内沿θ+π/2,即与θ角方向成直角方向上形成高亮度的谱线。

图2 边缘特征及其频谱方向Fig.2 Linear feature and its frequency spectrum distribution

3.2 图像线状特征的频率分析

线状特征在频谱中的方向特性确定之后,另一个关键是得到特征在频谱的这一方向上的频谱位置。遥感图像目标繁多,其中包含周期性成分、非周期性成分、噪声和背景,利用傅里叶变换可以将图像按频率进行分解,使得不同的周期成分在频谱图中很好地反映出来。傅里叶变换得到的频域,也可狭义地叫做“正弦波域”,是把实际信号分解为一组无限正弦波之叠加,而且这组正弦波是存在且唯一的。其中频谱上的每一点对应这一方向的一个频率值,就是一个由这一方向原信号分解出的某一正弦波的频率值,这一点频谱亮度的大小取决于这个方向上相同频率正弦波幅度能量的叠加。频谱中心为图像的零频分量,代表图像的零次谐波,依次从频谱中心沿各个方向向外,离中心点越近,表示频率越低;距中心点越远,则频率越高。图像的频谱四周为图像的高频所在,反映图像的细节,虽然在能量上较图像中心的低频能量相对较少,但确是图像的线状特征等关键信息所在。

图3(a)为实际遥感图像,图3(b)为图3(a)经傅里叶变换后的图像频谱,从频谱中看出频谱存在两个方向的明显亮线,为和原图像线状信息的角度差为90°的两个方向的频率能量叠加,这也验证了前面频谱方向的结论。除了对角线方向两条明显谱线以外,在水平和垂直方向也有互相垂直的两条谱线,为原图像线状信息在水平和垂直方向的频谱分量。为得出频谱中能量峰值的角度方向即亮线所在方向,采用对频谱进行角向能量采样的方法,图3(c)为图像角向能量曲线,从角向能量曲线得出,图像的能量峰值角度为41°、130°。

图3 图像频谱和图像的角向采样曲线Fig.3 Image magnitude spectrum and angle distribution curves

假设某一线状特征信息可分解出一系列高低频混合的正弦波分量,如果得到这些正弦波的频率值信息,也就得到此线状特征的频率分布。为获得线状特征的频率值分布范围,这里引入信号处理思想,因为线状特征处的灰度分布为脊状或阶跃状局部突变信号,为获得局部短时信号的频率信息,采用对特征信号进行局部加高斯窗的方法,来进行短时信号的傅里叶频谱分析。

信号的短时傅里叶变换可定义为

信号 f(t)在时间t的短时傅里叶变换就是信号乘上一个以t为中心的分析窗g(t-u)所作的傅里叶变换[18]。可以称短时傅里叶变换为信号f(t)在时间t附近的“局部谱”。

为模拟线状信号的频率分布,在图3(a)的能量峰值方向41°取一条灰度剖面线,以剖面线为对象,用剖面线的频率表征图像的频率变化,以实现面向对象的特征识别思想。图4(a)为模拟线状信号,图4(b)为实际剖面线信号,将其变换后取频谱中虚部的正频率部分表示,图4(c)、(d)分别为模拟信号和实际边缘特征的频谱,整体如图4所示。

图4(b)剖面线为变换前的线状特征信号,信号长度600点,将每个像素作为一个采样点,采样周期 T=1/600 s,这样横轴总时间长度为1 s,为图像的横轴坐标范围,得到的频率也为实际频率,这样就将灰度剖面线的频率和图像频率关联起来。从模拟和实际频谱图上看出,图像的阶跃边缘处在频率分布平面中表现为竖直频谱,空域的阶跃线状处和竖直的频谱分布的集中位置对应准确。在灰度变化率相对较小的线状处,频率分布为相对高频;随着灰度变化率增加,在灰度变化率最高的线状处,对应的频率分布为绝对高频,线状特征的高频几乎可达到整个图像的最高频率。从试验分析的极限类推中可得出一个结论,如果图像很平滑,没有线状特征存在,那么频谱就只有直流分量,随着线状尖锐程度的增加,线状特征对应的频谱分布向高频延伸,直到线状变为垂直线状时,灰度变化率最大,这时对应线状所具有的高频成分也达到最高,因此,经过线状特征的频率分布谱分析,基本可以得到某一线状特征处的频率值分布范围,并可以从频率纵坐标轴上取一个截止频率值,将感兴趣的线状特征包括进来。这样也就使得线状特征在频域中提取具有便捷性和自适应性。

图4 线状特征的频率分布谱Fig.4 Image’s profile curve and frequency spectrum expression of linear feature

3.3 G abor滤波器设计

Gabor滤波器是1964年由 Gabor引入的窗口傅里叶变换。Gabor滤波器因能够模拟人眼视觉系统来分析图像的特征且具有检测方向性特征的功能而被广泛应用于图像的增强、提取等方面[19]。基于上文关于线状特征频谱方向和频率分析获得的方向和频率参数,构造特定的 Gabor滤波器进行线状特征的提取,Gabor滤波器的数学表达式为

式中,u和v分别为沿x和y坐标轴的频率

上式是空间域坐标旋转的结果。应用φ为高斯函数g(x,y)的方向角,高斯函数的长轴方向是沿 x轴正向顺时针旋转,当φ=0时,高斯函数的长轴平行于x轴。g(x,y)为下式所示的高斯函数

式中,σx和σy为高斯函数的方差,决定滤波器的带宽;λ是高斯环的长径比(如果λ=1,g(x,y)为圆对称);σ是尺度参数。若令σx=σ,σy=σ λ,那么上式可以改写为如下形式

h(x,y)的傅里叶变换 H(u,v)为

4 线状特征检测及结果分析

Gabor滤波器构建的核心是滤波参数的选择,根据前面给出的 Gabor滤波器表达式,需要确定的参数为σ、θ、λ、F,F= u20+v20,比较关键的参数是方向参数和频率参数。利用前面频谱分析确定出的最佳参数,代入 Gabor滤波器的频域表达式,用频域滤波器和图像频谱做乘积运算,再将运算结果反变换到空域即得到滤波结果图像。

由于Gabor滤波器的实部对纹理敏感,而虚部对线状特征敏感[20],因此选择 Gabor滤波器的虚部进行特征的滤波提取。但在有效的特征区域,Gabor滤波器虚部的一些滤波结果会出现模值相对较大的负数,因此对Gabor滤波器的滤波结果取模值。另外,由于 Gabor滤波器的方向和空域图像的特征方向相同时,基本没有滤波结果,只有和滤波器垂直时才会取得最佳效果,因此需要确定适当的提取角度,这里取θ=41°、130°,根据频率分析结果,取 F=0.4,据求解[21],其中为倍频程,求得σ= 1.2,λ为滤波器长短轴比,根据多次迭代尝试,取λ=0.6。用这些参数设计出来的 Gabor滤波器相应的滤波提取结果如图5所示。

图5 不同滤波参数的滤波结果Fig.5 Linear feature detect result based on different parameter

首先考察滤波器提取结果的参数影响,图5为原图像和两个不同方向滤波器的线状特征检测结果,固定方向的滤波器将和其垂直的线状信息完全提取,可见应用确定方向和频率的滤波器不仅减少盲目计算,而且减少其他方向非特征信息的干扰。图6(c)为两个不同方向滤波结果叠加得到的本文设计方法提取结果。

图6 不同滤波器检测结果比较Fig.6 Linear feature detect result comparison with Sobel and Canny operator

图6应用本文确定的最佳参数提取结果和Sobel算子、Canny算子提取结果对比,从图6(a)中发现,Sobel算子的检测结果受图像局部对比度的影响最为明显,对比度大的线状信息产生很强的梯度响应,对比度小的线状信息则产生弱的梯度响应。而图6(b)中,Canny算子因引入高斯平滑函数,容易造成图像模糊,检测得到的线状特征较宽。本文算法提取的线状特征较细致、尖锐,并且对阶跃和屋脊状线状特征都有明显地响应,边界连续光滑,在边缘交汇的角点处,因来自于两个方向的能量重复叠加,故呈现出局部亮点。源于频率选择策略,非主要信息被适当弱化,突出了感兴趣线状特征。

5 结 论

应用QuickBird Pan图像进行线状特征提取试验的提取效果表明本方法的主要优势有:①通过频谱分析得到谱线的方向特征和频率分布特征,使得之后的频域滤波变得有理有据,匹配的滤波器简化了运算,也使本方法具有空域不具备的速度优势,算法的进一步扩展可实现海量数据信息的自动、快速提取;②根据角向采样和频谱分布的谱估计,可实现滤波参数的自动化获取,方法的自适应性较强;③根据频谱分析得到的方向和频率参数,针对不同需要可通过滤波器参数的改变提取不同方向、尺度的线状特征,具有多方向、多尺度特性。

[1] WIL KINSON G G.Recent Development in Remote Sensing Technology and the Importance of Computer Vision Analysis Techniques[C]∥Proceedings of Concerted Action MAVIRIC.Kingston:[s.n.],1999.

[2] CHENG Chengqi,MA Ting.Automatic Recognition of Landscape Linear Features from High-resolution Satellite Images[J].Journal of Remote Sensing,2003,7(1):26-30.(程承旗,马廷.高分辨率卫星影像上地物线性特征的自动识别[J].遥感学报,2003,7(1):26-30.)

[3] SOBEL I.Neighbourhood Coding of Binary Images for Fast Contour Following and General Array Binary Processing [J].Computer Graphics and Image Processing,1978,8: 127-135.

[4] CANNYJ F.A Computational Approach to Edge Detection [J].IEEE Transactions on Pattern Analysis and Machine Intelligence,1986,8(6):679-698.

[5] HOSSEINI R.Fourier Analysis of Plain Weave Fabric Appearance[J].Textile Research Journal,1995,65(11): 676-683.

[6] KOVESI P.Invariant Measures of Image Features from Phase Information[D].Perth:University ofWestern Australia,1996.

[7] DELENNE C,RABATEL G,AGURTO V,et al.Vine Plot Detection in Aerial Images Using Fourier Analysis[C]∥Proceedings of1stInternationalConference OBIA. Salzburg:ISPRS,2006.

[8] WASSENAAR T,ROBBEZ-MASSON J M,ANDRIEUX P,et al.Vineyard Identification and Description of Spatial Crop Structure by Per-field Frequency Analysis[J].International JournalofRemote Sensing,2002,23(17): 3311-3325.

[9] LUO Y,SAUDI A,MARHOON M.Generalized Hilbert Transform and Its Applications in Geophysics[J].The Leading Eedge,2003,22(3):198-202.

[10] MORRONE M C,BURR D C,ROSS J,et al.Mach Bands Are Phase Dependent[J].Nature,1986,324: 250-253.

[11] XIAO Pengfeng,FENG Xuezhi,ZHAO Shuhe,et al. Segmentation of High-resolution Remotely Sensed Imagery Based on Phase Congruency[J].Acta Geodaetica et Cartographica Sinica,2007,36(2):146-151.(肖鹏峰,冯学智,赵书河,等.基于相位一致的高分辨率遥感图像分割方法[J].测绘学报,2007,36(2):146-151.)

[12] KOVESI P.Image Features from Phase Congruency[J]. Journal of Computer Vision Research,1999,1(3):1-26.

[13] BEZALEL E,EFRON U.Efficient Face Recognition Method Using a Combined Phase Congruency/Gabor Wavelet Technique[C]∥Proceedings of the SPIE.San Diego:SPIE,2005:437-444.

[14] CHEN Shupeng.Investigation Studies on Geo-informatic Tupu[M].Beijing:Commercial Press,2001.(陈述彭.地学信息图谱探索研究[M].北京:商务印书馆,2001.)

[15] TAN T N.Texture Edge Detection by Modeling Visual Cortical Channels[J].Pattern Recognition,1995,28(9): 1283-1298.

[16] KIM D S,LEE S U.Image Vector Quantizer Based on a Classification in the DCT Domain[J].IEEE Transactions on Communications,1991,39(4):549-556.

[17] TAN Yanying,DONG Zhixin.Spectrum Analyzing of Image Edges and Application in Transform Coding[J].Acta Electronica Sinica,1995,23(9):61-65.(谭雁英,董志信.图像边缘信息的谱分析及其在变换编码中的应用[J].电子学报,1995,23(9):61-65.)

[18] CHEN Xuehua,HE Zhenhua,HUANG Deji.Seismic Data Edge Detection Based on Higher-order Pseudo Hilbert Transform[J].Progress in Geophysics,2008,23(4): 1106-1110.(陈学华,贺振华,黄德济.地震资料的高阶伪希尔伯特变换边缘检测[J].地球物理学进展,2008,23 (4):1106-1110.)

[19] DEL ENN E C,RABA TEL G,DESHA YES M.An AutomatizedFrequency Analysis for Vine Plot Detection and Delineation in Remote Sensing[J].IEEE Geoscience and Remote Sensing Letters,2008,5(3):341-345.

[20] FU Yipping,LI Zhineng,YUAN Ding.Edge Detection with Optimized Gabor Filter[J].Journal of Computeraided Design and Computer Graphics,2004,16(4):481-486.(傅一平,李志能,袁丁.基于优化设计 Gabor滤波器的边缘提取方法[J].计算机辅助设计与图形学学报, 2004,16(4):481-486.)

[21] WU Gaohong ZHAN G Yujin L IN Xinggang.Optimal Gabor Filter Design for Bi-textured Image Segmentation [J].Acta Electronica Sinica,2001,29(1):48-50.(吴高洪,章毓晋,林行刚.分割双纹理图像的最佳Gabor滤波器设计方法[J].电子学报,2001,29(1):48-50.)

Linear Feature Detection for High-resolution Remotely Sensed Imagery in Frequency Domain

ZHOU Liguo1,FENG Xuezhi2,XIAO Pengfeng2,MA Weichun1
1.Department of Environmental Science and Engineering,Fudan University,Shanghai 200433,China;2.School of Geographic and Oceanographic Sciences,Nanjing University,Nanjing 200091,China

Information of linear features in high-resolution remotely sensed imagery is very prolific,consequently the immoderate specializations and influence of noise disturbances are serious,this made it difficult to detect and recognize the linear features of landscape targets.A method is discussed based on direction and frequency characters through designing frequency domain filter for the linear feature automatic recognition and detection, owed to Fourier transformation method transform the image into frequency domain.Analyzed the relations between the linear feature and its spectrum line,also the relations between the linear feature and its frequency character, according of the direction and frequency selected above the Gabor wave filter the image linearity feature extracting was constructed.The relevant extract experiments were tested with QuickBird image,the result indicates that the linearity characteristic extracted in this method is fairly good,it provided a reference for high-resolution remotely sensed imagery linear feature extracting.

linear feature;frequency domain;Gabor filter;characteristic detect;high-resolution satellite imagery

ZHOU Liguo(1980—),male,PhD,lecturer, majors in digital image processing of remote sensing and environmental remote sensing.

1001-1595(2011)03-0312-06

TP751

A

国家863计划(2008AA12Z106);国家自然科学基金(40801166;41001234)

(责任编辑:宋启凡)

2009-12-17

2010-10-28

周立国(1980—),男,博士,讲师,主要研究方向为遥感数字图像处理及环境遥感。

E-mail:Zhouli-guo@tom.com