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

趙好好，馮學(xué)智，肖鵬峰

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

1 引 言

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

遙感圖像城市綠地信息識別和提取研究中，光譜特性是主要的研究方法。植被指數(shù)是早期綠地特征識別的主要手段，其中NDVI應(yīng)用最為廣泛。也有學(xué)者則利用光譜混合分析進(jìn)行城市綠地覆蓋類型的識別[2-3]。隨著圖像分辨率的提高，地物紋理、輪廓、形狀等特征得到重視，面向?qū)ο蟮姆椒ㄔ诔鞘芯G地信息提取中開始廣泛應(yīng)用，基于像元的分類方法逐漸向基于特征的提取技術(shù)過渡[4]。

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

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

2 研究數(shù)據(jù)

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

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

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

圖1 研究區(qū)遙感圖像及頻譜特征

3 研究方法

3.1 頻譜分析

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

(1)

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

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

(2)

(3)

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

3.2 特征識別

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

圖2 研究區(qū)頻譜分析

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

圖3 不同中心頻率的濾波結(jié)果

圖4 圖2(d)的一階導(dǎo)數(shù)曲線

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

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

圖5 研究區(qū)截止頻率建立

圖6 研究區(qū)遙感圖像頻譜分析

3.3 濾波器設(shè)計

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

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

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

(4)

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

(5)

二維Gabor函數(shù)的傅里葉變換為：

(6)

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

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

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

(8)

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

根據(jù)頻譜分析結(jié)果，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。

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

假設(shè)理想的階躍邊緣模型為：

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

(9)

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

圖7 兩個方向的空間常數(shù)取值

3.3.2 道路特征濾波器

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

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

4 研究結(jié)果與評價

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

圖8 Gabor濾波結(jié)果

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

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

圖10 研究區(qū)道路綠地覆蓋輪廓

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

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

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

圖11 3種結(jié)果的比較

5 結(jié)束語

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

本文的研究區(qū)位于城市的主、次干道，道路的類型為直線型。其他類型的道路如彎曲形，在進(jìn)行濾波時需要調(diào)整濾波參數(shù)，因此關(guān)于其他類型的城市道路綠地覆蓋還有待進(jìn)一步研究。

參考文獻(xiàn)：

[1] 宮鵬，黎夏，徐冰.高分辨率影像解譯理論與應(yīng)用方法中的一些研究問題[J].遙感學(xué)報，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的城鎮(zhèn)綠地信息動態(tài)監(jiān)測研究[J].上海交通大學(xué)學(xué)報，2009，27(1)：1-6.

[5] 童慶禧，薛永祺.空間遙感技術(shù)發(fā)展與應(yīng)用的深化[C].中國空間科學(xué)學(xué)會第七次學(xué)術(shù)年會，2009.

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

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

[8] 吳桂平，肖鵬峰，馮學(xué)智，等.基于光譜空間變換的遙感圖像目標(biāo)探測方法研究.光譜學(xué)與光譜分析，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] 傅一平，李志能，袁丁.基于優(yōu)化設(shè)計Gabor濾波器的邊緣提取方法[J].計算機(jī)輔助設(shè)計與圖形學(xué)學(xué)報，2004，16(4)：481-486.