基于頻域特征的遙感圖像城市道路綠地覆蓋輪廓提取

趙好好，馮學智，肖鵬峰

(1.南京大學 地理信息科學系，南京 210023；2.南京信息工程大學 遙感學院，南京 210044)

1 引 言

道路綠地是城市綠地的重要組成部分，具有改善城市生態環境，調節小氣候，提高城市居民生活舒適度等功能，城市道路綠地覆蓋信息的提取對城市建設和管理具有重要意義。隨著衛星遙感圖像空間分辨率的提高，基于像元的遙感圖像分類技術開始向基于對象的特征識別和提取技術過渡，這使得城市地物特征信息的提取成為可能[1]。

遙感圖像城市綠地信息識別和提取研究中，光譜特性是主要的研究方法。植被指數是早期綠地特征識別的主要手段，其中NDVI應用最為廣泛。也有學者則利用光譜混合分析進行城市綠地覆蓋類型的識別[2-3]。隨著圖像分辨率的提高，地物紋理、輪廓、形狀等特征得到重視，面向對象的方法在城市綠地信息提取中開始廣泛應用，基于像元的分類方法逐漸向基于特征的提取技術過渡[4]。

頻譜特征是進行圖像識別的重要依據，通過建立空間紋理結構特征的頻譜分析模型來計算不同大小地物類型的紋理頻譜，可有效提高圖像地物識別能力[5]。通過對不同地物進行頻譜特征分析，建立地物頻譜識別標志，可以有效提高地物信息提取的精度和自動化水平[6]。

在頻率域遙感圖像分析中，王珂等對城市河道線性信息進行了識別和提取[7]，吳桂平等對城市人工建筑進行了識別[8]。而城市道路綠地覆蓋，與城市道路之間具有特殊的幾何關系，本文研究在幾何約束條件下的城市地物(道路綠地覆蓋)頻域信息提取方法。通過對高分辨率遙感圖像進行傅里葉變換，建立對應道路綠地覆蓋輪廓的截止頻率，設計合適的濾波器，并借助道路中心線幾何輔助信息，實現道路綠地覆蓋的信息提取。

2 研究數據

本文所使用的高分辨率遙感數據為南京市主城區2007年7月獲取的QuickBird圖像，包含4個多光譜波段和1個全色波段，融合后的空間分辨率為0.61m。圖像成像質量良好，無云層覆蓋。

圖1(a)為研究區近紅外波段圖像，位于中山北路與云南路交叉地帶，圖像大小為558×558像素。道路綠地與道路方向一致，在圖像中沿著兩個方向分布，且基本垂直。在幾何關系上，道路綠地位于道路紅線以內，如果道路中心線能夠確定，做道路中心線的緩沖區，可提取道路綠地覆蓋的位置信息。

圖1(a)中，近紅外波段道路綠地與道路的光譜差異顯著，道路綠地的亮度值平均為650，而道路亮度均值為310。道路信息包含道路色調信息和道路邊緣信息。本實驗區中，由于道路被行道樹覆蓋比較嚴重，所以道路邊緣信息不顯著。而道路的色調信息在圖像上亮度最低，與建筑、行道樹差別顯著，較易區別。

圖1 研究區遙感圖像及頻譜特征

3 研究方法

3.1 頻譜分析

假設一幅大小為M×N的遙感圖像，其二維離散傅里葉變換如下式表示：

(1)

其中，x，y為空域坐標，u，v為頻率。圖像的頻譜值可表示為|F(u,v)|。

頻譜分析通常在極坐標系中進行，假設極坐標系頻譜函數為S(r,θ)，r和θ是坐標變量。通過輻射掃描法可分析徑向分布特征和角向分布特征，如下式所示：

(2)

(3)

研究區遙感圖像的頻譜圖如圖1(b)所示。從頻譜圖上可以看到，圖像的中高頻豐富，在與道路綠地垂直的方向上，存在兩條顯著的譜線。另外在水平與豎直方向上，也存在少量譜線，這與圖像中的建筑等地物有關。圖1(c)徑向分布曲線表明，隨著頻率的增加，頻譜逐漸降低，沒有出現顯著的諧波成分。在頻率20周/像元，如圖中箭頭所示，能量下降的幅度出現了顯著轉折。圖像的色調信息主要集中在直流中心，20周/像元反映了低頻位置。圖1(d)角向分布曲線中，出現了兩個峰值能量方向，分別為45°和135°，兩個方向均垂直道路綠地的方向。圖像的徑向分布和角向分布主要反映了在頻率和方向上的能量累積值，而與道路綠地覆蓋輪廓特征相關的頻率細節信息則無法體現。

3.2 特征識別

為了詳細分析對應道路綠地覆蓋輪廓的截止頻率，取圖1(a)中小區域256×256尺寸的圖像T1，如圖2(a)所示。圖2(b)為其頻譜圖，從角向分布曲線圖2(c)可以發現，其峰值角向能量為135°，與道路綠地方向垂直。現做頻譜圖135°方向的剖面線，因頻譜圖的對稱性質，這里僅取其一半顯示。如圖2(d)圓圈所示，在頻率20、37、47周/圖像出現了顯著的波谷位置。頻譜曲線中的波谷位置，反映了此頻率處的頻譜能量值較低，同時也是能量變化顯著的區域。

圖2 研究區頻譜分析

進一步分析這3個頻率所反映的圖像特征，分別以20、37、47周/圖像為中心頻率，對原圖像進行半徑為5的頻域理想濾波，并將濾波后的結果與原圖像進行對比。取原圖像T1的第89行以及濾波后圖像的第89行比較，如圖3所示，發現中心頻率為37周/圖像的濾波結果與道路綠地的輪廓位置符合最好。圖2(d)中頻譜剖面曲線的一階導數曲線，如圖4所示，頻率37周/圖像是曲線的最高峰值頻率。

圖3 不同中心頻率的濾波結果

圖4 圖2(d)的一階導數曲線

用同樣的方法，截取原圖像1(a)的另外3幅小區域256×256大小圖像T2、T3、T4進行分析，如圖5所示。發現3幅圖像對應道路綠地覆蓋輪廓的最佳截止頻率分別為64、44、32周/圖像。圖5(b)、圖5(d)、圖5(f)分別為圖像T2、T3、T4頻譜最大能量方向一階導數曲線，這幾個頻率均處于曲線的較高峰值位置。據此判斷，做頻譜圖中最大能量方向一階導數曲線，曲線的峰值位置可以作為判斷圖像輪廓截止頻率的依據。37、64、44、32周/圖像，位于頻譜半徑的中頻0.32～0.64階段，也即對應道路綠地覆蓋輪廓的圖像截止頻率應主要分布在此頻段范圍。

對研究區原圖像(圖1(a))進行頻譜特征識別，在45°、135°兩個方向的頻譜曲線如圖6(a)、圖6(b)所示，其一階導數曲線如圖6(c)、圖6(d)所示。取中高頻區域，也即頻率半徑的0.32～0.64頻段的頻譜能量峰值，那么45°方向的在頻率122周/圖像出現最高峰值；135°方向在頻率為156周/圖像出現最高峰值。因此，45°方向對應道路綠地覆蓋輪廓的截止頻率為122周/圖像；135°方向截止頻率為156周/圖像。

圖5 研究區截止頻率建立

圖6 研究區遙感圖像頻譜分析

3.3 濾波器設計

3.3.1 道路綠地覆蓋輪廓特征濾波器

為了提取特定方向、特定頻率的圖像信息，需要設計適合圖像特征的濾波器。Gabor濾波器與哺乳動物視網膜神經細胞的接收場模型相吻合，它能夠同時在時域和頻域獲得最佳的局部化[9-10]，因此利用Gabor濾波器提取道路綠地的輪廓信息。一般情況下，二維Gabor函數h(x,y)可表示為：

h(x,y)=g(x′,y′)exp[2πj(Ux,Vy)]

(4)

其中，(U，V)表示特定的空間中心頻率，g(x，y)為高斯函數，(x′，y′)是(x，y)旋轉θ度角，即：

(5)

二維Gabor函數的傅里葉變換為：

(6)

其中，(U′，V′)和(u′，v′)表示類似于式(5)的角度旋轉。

Gabor函數實際上是復數形式，將其分解可表示為：

h(x,y)=R(x,y)+jI(x,y)R(x,y)=g(x,y)cos[?(xcosθ+ysinθ)]I(x,y)=g(x,y)sin[?(xcosθ+ysinθ)]

(7)

如式(7)所示，Gabor函數可以分成奇函數和偶函數兩個部分，實部是個偶函數，用實部對圖像濾波僅能起到圖像平滑，為了減少運算量，實現邊緣的快速提取，僅使用Gabor函數的奇部進行道路綠地覆蓋輪廓的圖像濾波。一個θ方向的二維Gabor奇函數可具體表示為：

(8)

其中，ω為角頻率矢量，ω=2πf，f為頻率，σ為空間常數。

根據頻譜分析結果，45°方向邊緣截止頻率為122周/圖像，為了計算方便，將其除以頻率半徑，取其歸一化頻率，即122/279=0.44周/像元，那么奇Gabor濾波器角頻率取值為ω=2πf=2π×0.44=2.76。135°方向的邊緣截止頻率為156/279=0.56周/像元，則奇Gabor濾波器角頻率取值為ω=2πf=2π×0.56=3.52。

除了方向和角頻率，空間常數也是奇Gabor濾波器設置的一個重要參數。

假設理想的階躍邊緣模型為：

U(x,y)=1(y≤0)U(x,y)=0(y>0)

(9)

做奇Gabor函數與理想邊緣的卷積，并固定角頻率ω和方向θ，則卷積函數隨著空間常數σ變化曲線如圖7所示。45°方向卷積函數隨著空間常數的變化曲線見圖7(a)，當σ=0.70，函數取最大值，因此45°方向濾波器參數設置為σ=0.70，ω=2.76，θ=π/4；135°方向卷積函數隨著空間常數的變化曲線見圖7(b)，當σ=0.55時函數取得最大值，因此135°方向濾波器參數為：σ=0.55，ω=3.52，θ=3π/4。

圖7 兩個方向的空間常數取值

3.3.2 道路特征濾波器

根據研究區近紅外波段圖像道路的基本特征，道路信息提取將主要提取色調信息，然后提取其脊線作為道路的中心線。據參考文獻[7]，色調信息集中在頻譜圖的低頻直流中心，因此道路色調信息通過對低頻信息進行濾波得出。結合圖1(a)，道路信息的直線方向特征顯著，但色調最暗。因此，本文對原圖像進行反變換，突出道路色調信息，然后進行低頻濾波提取道路色調信息。

用式(6)Gabor濾波器進行頻域濾波，濾波器中心設置為(0，0)，方向設置分別為45°和135°。在圖1(c)中看到，隨著頻率增加，頻譜能量迅速下降，當頻率為20周/圖像時，下降速度顯著減慢，因此沿著長軸方向的σu取值為20周/圖像。為了突出線性信息，令垂直長軸方向的σv為6周/圖像。

4 研究結果與評價

對圖像進行基于3.3.1參數設置的頻域濾波，結果如圖8所示。從濾波結果可以看出，道路綠地覆蓋方向上的輪廓信息提取較為完整，閉合性好。同時部分非道路綠地、人工建筑等其他地物頻率與確立的截止頻率接近，也有部分提取出來。

圖8 Gabor濾波結果

為了限制道路綠地的位置信息，利用3.3.2關于道路特征的濾波器，進行道路中心線輔助信息的提取。濾波結果如圖9所示，圖9(c)中道路的線性特征顯著，通過二值分割和脊線提取，道路中心線如圖9(d)所示。文中方法提取的道路中心線連續，沒有出現斷點，是直線型道路中心線提取較為實用的方法。

根據實地調查，研究區道路45°方向寬20m，約32像素，135°方向寬27m，約45像素。分別在道路中心線兩側設置緩沖區，對邊緣提取圖像進行緩沖區限制之后的結果如圖10所示。圖10(a)和圖10(b)為兩個方向的道路綠地覆蓋輪廓，圖10(c)為兩個方向之和，圖10(d)為二值化后的結果。從圖中看出，建筑、道路等其他地物完全去除，道路綠地很好地檢測出來，45°方向輪廓完整清晰，整體閉合性較好；135°方向輪廓清晰，較為完整，受圖像道路綠地覆蓋特征影響，右下角有部分邊緣沒有閉合，如圖中圓圈所示。

圖10 研究區道路綠地覆蓋輪廓

將二值化道路綠地輪廓圖像經過數學形態學開運算，去除斑點噪聲，結合圖像NDVI，最后結果如圖11(a)所示。道路綠地覆蓋輪廓曲線閉合，完整，內部基本無噪聲影響。非道路綠地部分已被較好地去除。

為了進行比較分析，做圖像的Canny邊緣提取(圖11(b))。由于Canny算法不是針對圖像特征提取的邊緣檢測算法，沒有顧及道路綠地覆蓋輪廓的線性特征，因此檢測結果中，部分受建筑物影響的道路綠地覆蓋輪廓不連續，沒有較好地提取出來，如圖11(b)中圓圈所示部分。而且Canny算子檢測的是圖像的邊緣點，這不僅包含道路綠地覆蓋輪廓，還包含許多具有邊緣特征的內部噪聲。

研究區NDVI結果的二值化圖像如圖11(c)所示。NDVI主要依據綠地的光譜特征，能將大多數與綠地光譜差異大的地物去除。但是沒有道路中心線輔助信息，利用NDVI提取的綠地結果中，存在大片的非道路綠地。盡管NDVI在中低分辨率綠地覆蓋信息提取中具有優勢，但對于高分辨率遙感圖像說來，NDVI結果中的綠地輪廓過于光滑，不能充分顯示與圖像結構相關的細節信息。

圖11 3種結果的比較

5 結束語

本文通過對高分辨率遙感圖像城市道路綠地覆蓋進行頻譜分析，找到對應道路綠地覆蓋輪廓的截止頻率，設計奇Gabor濾波器提取兩個方向上的線性輪廓信息。奇Gabor濾波器提取出的道路綠地輪廓線性特征顯著，但在提取過程中受到了噪聲影響，攜帶了部分其他地物方向特征信息。通過道路中心線進行道路綠地覆蓋范圍的限制，最終得到了道路綠地覆蓋輪廓圖。本文的研究方法還可以用到其他受線性條件約束的地物特征提取中。

本文的研究區位于城市的主、次干道，道路的類型為直線型。其他類型的道路如彎曲形，在進行濾波時需要調整濾波參數，因此關于其他類型的城市道路綠地覆蓋還有待進一步研究。

參考文獻：

[1] 宮鵬，黎夏，徐冰.高分辨率影像解譯理論與應用方法中的一些研究問題[J].遙感學報，2006，10(1)：1-5.

[2] SMALL C.High spatial resolution spectral mixture analysis of urban reflectance [J].Remote Sensing of Environment，2003，88(1)：170-186.

[3] TOOKE T R，COOPS N C，GOODWIN N R.Extracting urban vegetation characteristics using spectral mixture analysis and decision tree classifications[J].Remote Sensing of Environment，2009，113(2)：398-407.

[4] 申廣榮，錢振華，徐敬敬，等.基于eCognition的城鎮綠地信息動態監測研究[J].上海交通大學學報，2009，27(1)：1-6.

[5] 童慶禧，薛永祺.空間遙感技術發展與應用的深化[C].中國空間科學學會第七次學術年會，2009.

[6] 葉澤田.頻譜段圖像及其應用的探討[J].環境遙感，1993，8(2)：139-146.

[7] 王珂，肖鵬峰，馮學智，等.基于頻域濾波的高分辨率遙感圖像城市河道信息提取.遙感學報，2013，17(2)：269-285.

[8] 吳桂平，肖鵬峰，馮學智，等.基于光譜空間變換的遙感圖像目標探測方法研究.光譜學與光譜分析，2013，33(3)：741-745.

[9] MEHROTRA R，NAMUDURI K R，RANGANATHAN N.Gabor filter-based edge detection [J].Pattern Recognition，1992，25(12)：1479-1494.

[10] 傅一平，李志能，袁丁.基于優化設計Gabor濾波器的邊緣提取方法[J].計算機輔助設計與圖形學學報，2004，16(4)：481-486.