一種頻域高分辨率遙感圖像線狀特征檢測方法

周立國,馮學智,肖鵬峰,馬蔚純

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