陳 洪， 陶 超， 鄒崢嶸， 邵 磊

中南大學(xué)地球科學(xué)與信息物理學(xué)院，長沙410083

隨著遙感影像分辨率的逐步提高，居民區(qū)提取已成為一個熱門課題.準(zhǔn)確而快速地提取居民區(qū)可以為土地管理、城市規(guī)劃等市政部門在做土地利用現(xiàn)狀調(diào)查和宏觀規(guī)劃等方面的工作時提供重要的決策支持.在遙感影像中，居民區(qū)不僅具有覆蓋面廣、地物信息豐富的特點，而且它是一個隨著時間動態(tài)變化的區(qū)域.人工提取居民區(qū)雖然可以保持較高的準(zhǔn)確度，但是僅僅依靠人工的監(jiān)測和分割提取居民區(qū)不但耗時而且成本高.因此，自動居民區(qū)提取方法越來越受到關(guān)注.

目前，許多研究學(xué)者提出了關(guān)于居民區(qū)提取的方法.文獻[1]首先通過數(shù)學(xué)形態(tài)學(xué)方法從衛(wèi)星影像上提取結(jié)構(gòu)信息，然后利用神經(jīng)網(wǎng)絡(luò)分類居民區(qū)，該方法需要人工訓(xùn)練樣本集的支持.文獻[2]融合傅里葉變換和神經(jīng)網(wǎng)絡(luò)分類算法提出了一種居民區(qū)提取方法，該方法同樣需要構(gòu)建訓(xùn)練樣本.文獻[3]采用灰度共現(xiàn)矩陣方法從高分辨率SAR影像上提取居民區(qū).文獻[4]利用多級分類器在Ikonos灰度影像上提取居民區(qū).文獻[5]結(jié)合結(jié)構(gòu)特征與圖形學(xué)理論從衛(wèi)星影像提取居民區(qū)，該方法使用了統(tǒng)計分類器.文獻[6]利用建筑物密度分類居民區(qū)，該方法是基于影像紋理信息而提出的.文獻[7]以紋理單元為模型檢測高分辨率遙感影像上的居民區(qū)，該方法利用支持向量機進行分類.文獻[8]融合了多種不同的分類器從全色衛(wèi)星影像上提取居民區(qū).文獻[9]利用局部邊緣分布信息檢測高分辨率遙感影像上的居民區(qū).分析現(xiàn)有居民區(qū)提取方法，發(fā)現(xiàn)存在以下兩方面的局限性：1)大多數(shù)方法基于監(jiān)督分類機制，需要大量訓(xùn)練樣本保證分類精度，干擾因素較多，自動化程度有限；2)現(xiàn)有的居民區(qū)提取方法強調(diào)影像整體的紋理、光譜等特征，卻忽略了居民區(qū)的局部特征，如邊緣特征.在影像中，居民區(qū)內(nèi)部包含大量的人工地物(如建筑物、道路等)，相對于其他非居民區(qū)(如耕地、林地等)具有更為明顯和豐富的邊緣特征，因此可以根據(jù)邊緣密度差異提取居民區(qū).

本文提出一種無監(jiān)督的高分辨率遙感影像居民區(qū)自動提取方法，首先將居民區(qū)內(nèi)的邊緣特征直線化表達，然后以高斯函數(shù)為模型將邊緣密度特征量化，最后通過閾值分割提取居民區(qū).

1 居民區(qū)邊緣密度特征提取

居民區(qū)與非居民區(qū)不同，其內(nèi)部包含了密集的建筑物和道路等地物，這些地物都具有很明顯的邊緣特征；非居民區(qū)(如耕地、林地等)內(nèi)的紋理較單一，沒有明顯的邊緣特征.因此，居民區(qū)具有更高的邊緣密度，根據(jù)邊緣密度特征差異就可以區(qū)分識別居民區(qū)與非居民區(qū).本文按照以下4個步驟提取影像上的邊緣密度特征：

步驟1 采用Mean Shift算法平滑原始影像，并檢測平滑影像上的邊緣；

步驟2 將所有邊緣擬合成直線段；

步驟3 計算每個像素點到所有邊緣線段的空間距離，并用集合存儲；

步驟4 根據(jù)距離集合，利用高斯函數(shù)構(gòu)建空間投票矩陣，最終獲取影像上的邊緣密度特征.

在復(fù)雜的高分辨率遙感影像中具有大量的無關(guān)噪聲，因此在提取邊緣之前需要做平滑去噪處理.Mean Shift濾波與常用的去噪算法(如中值濾波、高斯濾波等)相比具有良好的保邊緣性[10-11]，于是本文采用Mean Shift濾波平滑原始影像，平滑結(jié)果如圖1中的(b)所示.由圖1可知，Mean Shift具有良好的保邊緣性，因此在濾波后的影像中，大量紋理噪聲被有效地抑制，同時居民區(qū)內(nèi)的邊緣細節(jié)也保留得較好.在有效地濾除掉影像中的噪聲后，采用Canny算子提取影像中的邊緣.圖1中的(c)和(d)分別顯示了在原始影像和Mean Shift濾波影像中的邊緣提取結(jié)果，從中可以看出：直接對原始影像做邊緣提取，無論是居民區(qū)還是非居民區(qū)，其邊緣都非常密集，邊緣密度沒有區(qū)分度；而在Mean Shift濾波后影像中，非居民區(qū)內(nèi)的邊緣比較稀疏，但居民區(qū)仍然擁有較高的邊緣密度，因此可以利用影像中邊緣密度分布推斷影像中居民區(qū)的位置.

圖1 平滑及邊緣影像Figur e 1 Smoothing and edge images

假設(shè)經(jīng)過Canny邊緣提取得到一組邊緣特征集合E0={e1,e2,···,en}，影像中邊緣密度特征提取分為以下3個步驟：

步驟1 如圖2所示，將檢測到的邊緣特征擬合成一條或多條直線段，依次按照先左至右，后從上至下的順序，記錄下每條邊緣上的所有像素點，并令P0,P1,···,Pn為某一邊緣上的像素點.根據(jù)以下步驟進行直線擬合：

步驟1-1 令邊緣的起始像素點P0為固定點，P2為浮動點，兩點構(gòu)成直線段P0P2，計算頂點P1到P0P2的垂直距離a；

步驟1-2 設(shè)定最大容許距離值為ε，若a≤ε，則P3定義為新的浮動點，計算P1、P2到P0P3垂直距離的最大值，假設(shè)為b，并將其與ε再次進行比較；若a＞ε，則用直線段P0P2替代邊緣折線P0P1P2，然后將P2定義為新的固定點；

步驟1-3 重復(fù)執(zhí)行上一步直至遍歷所有邊緣像素點，最終獲取邊緣的直線表示形式并記錄下每條直線段兩個端點的像素坐標(biāo).直線擬合后的邊緣集合表示為，其中表示第i條邊緣上的第j條直線段，分別將它的兩個端點坐標(biāo)記為和

圖2 直線擬合算法Figure 2 Line straight algorithm

如圖3所示，將所有邊緣擬合成直線段后，每條邊緣均用隨機的顏色顯示.邊緣特征經(jīng)過直線擬合后，邊緣像素點總數(shù)減少了，邊緣直線表達也得到了簡化，因此出現(xiàn)了局部邊緣信息缺失的現(xiàn)象，但總體而言，居民區(qū)內(nèi)的邊緣密度仍然保持較高的水平，滿足居民區(qū)提取的要求.

圖3 邊緣直線擬合影像Figure 3 Edge straight line image

步驟2 計算圖像上每個像素點(假設(shè)共有R個像素點)到所有邊緣直線段(假設(shè)共有K條邊緣線段)的距離，用集合D={d11,d12,···,d1K,d21,d22,···,d2K,···,dR1,dR2,···,dRK}表示，其中dij表示第i個像素點到第j條邊緣的距離，計算公式如下：

式中，(xi,yi)為第i個像素點的坐標(biāo)，和分別為第j條邊緣直線段的兩個端點的像素坐標(biāo).

步驟3 利用步驟2計算得到的距離集合衡量邊緣密度特征.假設(shè)原始影像大小為M×N，定義投票矩陣V(M,N).如果像素點p(xi,yi)是居民區(qū)內(nèi)的像素點，由于居民區(qū)內(nèi)的邊緣更為密集，因此從統(tǒng)計學(xué)的角度來說，該像素點到所有邊緣直線段的空間距離的總和會相對較小，然后通過式(2)將空間距離轉(zhuǎn)化為投票數(shù)V(xi,yi)，則V(xi,yi)的值就更大；而非居民區(qū)內(nèi)的像素點到所有邊緣線段的空間距離的總和相對較大，故投票數(shù)就更少

式中，K為邊緣總條數(shù)，dij為根據(jù)式(1)計算得到的像素點到邊緣直線段的距離，σ為鄰近投票參數(shù)，3.2節(jié)對該參數(shù)進行了具體分析.

圖4為根據(jù)上述方法計算得到的邊緣密度特征圖，圖中紅色區(qū)域為擁有投票數(shù)較多的區(qū)域，即為邊緣密度較高的地方；藍色區(qū)域獲得的投票數(shù)最少，即為邊緣密度較低的地方.由圖4可知，居民區(qū)內(nèi)具有更高的投票數(shù)，離居民區(qū)中心越遠，投票數(shù)就越少.

圖4 邊緣密度特征圖Figure 4 Edge density features image

2 基于邊緣密度特征居民區(qū)提取

如圖4所示，邊緣密度特征圖中投票值較大的區(qū)域就可能為居民區(qū)，而非居民內(nèi)投票值普遍偏低，所以影像上所有像素點的投票值服從一個雙峰分布，則在兩個峰值之間必然存在一個最佳的分割閾值，本文參考Ostu方法[12]自動獲取該閾值：

1)對所有投票值進行排序，用集合Vs={v1,v2,···,vL}表示，令投票值為vi的像素點的數(shù)目為ni，則總的像素數(shù)位N=n1+n2+···+nL；

2)將集合Vs拆分兩個類V0={v1,v2,···,vk}和V1={vk+1,vk+2,···,vL}，然后分別計算兩個類的類內(nèi)方差、類間方差和總方差；

3)進行循環(huán)運算，直至式(3)取得最大值時判定取得最佳閾值k?

通過上述方法獲取邊緣密度特征圖的最佳分割閾值，然后進行圖像分割，如圖5中的(a)為閾值分割后的二值圖，白色區(qū)域即為居民區(qū)，5中的(b)用紅色曲線標(biāo)記居民區(qū)外圍邊界，可見提取結(jié)果與影像上的居民區(qū)具有很高的吻合度，說明本文方法具有可行性.另外，對于不同的實驗影像，利用Ostu方法獲取的分割閾值各不相同，因此本文方法是一個自適應(yīng)的居民區(qū)分割提取方法.

圖5 居民區(qū)二值及邊界線圖Figur e 5 Built-up area binary and boundary line images

3 實驗結(jié)果與分析

3.1 居民區(qū)提取結(jié)果分析及評價

整個算法流程已由前面的例子完整演示.為驗證本文算法的準(zhǔn)確性及魯棒性，下面采用多幅高分辨率遙感影像進行實驗.另外通過人工方式準(zhǔn)確地標(biāo)記出居民區(qū).與此同時，采用Paul Rosin的基于邊緣特征的顯著區(qū)域檢測方法[13]進行對比試驗，實驗結(jié)果如圖6所示，黃色區(qū)域人工標(biāo)記居民區(qū)結(jié)果，紅色區(qū)域為本文方法實驗結(jié)果，藍色區(qū)域為Paul Rosin方法實驗結(jié)果.用于實驗的10幅影像分辨率在0.5～4.0 m之間，場景均較復(fù)雜，包含各種具有密集邊緣特征的地形地貌情況，如規(guī)則田地、山地起伏地區(qū)，以及林地、耕地、河流交織區(qū)域，居民區(qū)提取結(jié)果良好，說明本文方法具有良好的實際性能.

為了度量居民區(qū)提取精度，參考文獻[14]的評價方法，首先分別統(tǒng)計以上10幅影像經(jīng)本文方法及Paul Rosin方法所提取的居民區(qū)像素總數(shù)Sauto，人工提取的理想居民區(qū)像素總數(shù)Smanual，兩者相同的像素總數(shù)Scommon，然后分別以Scommon與Smanual之比為準(zhǔn)確率Pd，以兩種方法錯檢像素總數(shù)(Sauto-Scommon)與Sauto之比為虛警率因子Pf，評價結(jié)果如圖7所示：

圖7中的(a)為本文方法和Paul Rosin方法提取以上10幅影像的居民區(qū)的精度對比圖，本文方法提取精度穩(wěn)定保持在一個較高的水平，而Paul Rosin方法則波動性較大；圖7中的(b)為對應(yīng)的虛警率因子圖，Paul Rosin方法虛警率普遍高于本文方法，說明本文方法具有較高的準(zhǔn)確度性.

3.2 實驗參數(shù)分析

圖7 提取結(jié)果精度及虛警率圖Figur e 7 Images of accuracy and false alarm rate

本文涉及的可調(diào)參數(shù)可用于邊緣特征量化的鄰近投票參數(shù)σ、Mean Shift平滑參數(shù)(包括空間域帶寬Bs和頻率域帶寬Br)以及Canny邊緣檢測參數(shù).其中Mean Shift平滑主要應(yīng)用在影像預(yù)處理階段，因此相關(guān)的平滑參數(shù)取值選擇在合理的取值范圍內(nèi)，提取結(jié)果精度能穩(wěn)定保持在較高的水平.另外，本文方法不要求獲取準(zhǔn)確的邊緣特征，且要求居民區(qū)內(nèi)邊緣密度明顯比非居民區(qū)高，因此Canny邊緣檢測的參數(shù)選擇對本文提取結(jié)果影響不大.綜上所述，對本文算法提取結(jié)果影響最大的參數(shù)為用于量化邊緣特征密度的鄰近投票參數(shù)σ.

以實驗的10幅影像為例，分辨率為0.5～4.0 m，逐一分析上述各參數(shù)對居民區(qū)檢測成功率的影響：1)圖8(a)給出了鄰近投票參數(shù)σ取值為10～50時，10幅實驗影像的平均居民區(qū)檢測成功率的變化情況，σ分別取值34時，平均居民區(qū)檢測成功率取得最大值，σ取值為26～34時，居民區(qū)檢測成功率變化不明顯，穩(wěn)定保持在88%以上.σ取值不能過大或者過小，因為不同分辨率的影像均有其相適應(yīng)的最佳σ取值范圍；2)圖8(b)為空間域帶寬Bs和頻率域帶寬Br與10幅實驗影像平均檢測民區(qū)成功率的關(guān)系，其中紅色折線表示Br=20,Bs取值10～30之間時，平均檢測成功率的變化情況，藍色折線表示Bs=18,Br取值10～30之間時平均檢測成功率的變化情況，分析可知，兩者取值在16～22之間時，居民區(qū)檢測成功率變化不大；3)圖8(c)為Canny邊緣檢測的相關(guān)參數(shù)與10幅影像的平均居民區(qū)檢測成功率之間的關(guān)系，邊緣檢測的閾值分割參數(shù)在MATLAB編程環(huán)境下自動獲取，另外一個可變參數(shù)為高斯濾波標(biāo)準(zhǔn)方差α，默認取值為，圖8(c)中α取值在0.5～3.0之間時，平均檢測率的變化情況，整體波動較小，均能保持在88%以上，說明參數(shù)α對實驗結(jié)果影響不大.綜合上述分析，當(dāng)以上參數(shù)選取在合理的取值區(qū)間內(nèi)時，本文方法提取居民區(qū)的精度能夠穩(wěn)定保持在較高的水平，分析結(jié)果表明本文算法對于參數(shù)變化具有良好的魯棒性.

4 結(jié)語

本文利用居民區(qū)邊緣密度特征并結(jié)合空間投票機制，提出一種新的高分辨率遙感影像居民區(qū)自動提取方法.實驗表明，該算法具有較好的準(zhǔn)確性和魯棒性，可以從場景復(fù)雜的高分辨率遙感影像上提取居民區(qū)，主要創(chuàng)新點是以邊緣特征為基礎(chǔ)，發(fā)展了一種新的直線擬合方法，并以高斯函數(shù)為模型構(gòu)建了一種新的空間投票機制，從而準(zhǔn)確有效地度量了直線擬合后的邊緣密度特征.

對于地物信息過于復(fù)雜的影像，例如當(dāng)居民區(qū)周圍存在過多零散分布的灌木叢或獨立樹木時，或者當(dāng)實驗影像上居民區(qū)與非居民區(qū)比例相差較大時，本文算法的提取結(jié)果仍然存在一些錯判問題.因此，下一步工作考慮在居民區(qū)提取過程中加入除邊緣特征外更多的特征信息以提高檢測的準(zhǔn)確度，同時進一步提高算法的運行效率，以實現(xiàn)既快速又準(zhǔn)確的高分辨率遙感影像居民區(qū)自動提取方法.

[1]BENEDIKTSSON J A,PESARESI M,ARNASON K.Classif ication and feature extraction for remote sensing images from built-up areas based on morphological transformations[J].IEEE Transactions on Geoscience and Remote Sensing,2003,41(9):1940-1949.

[2]WEIGL K,GIRAUDON G,BERTHOD M.Application of projection learning to the detection of builtup areas in spot satellite images[R].Rapports de recherche-INRIA,1993.

[3]FANG W,CHAO W,HONG Z.Residential area information extraction by combining China airborne SAR and optical images[C]//IEEE International Geoscience and Remote Sensing Symposium,2004.

[4]ZHONG P,WANG R.A multiple conditional random f ields ensemble model for built-up area detection in remote sensing optical images[J].IEEE Transactions on Geoscience and Remote Sensing,2007,45(12):3978-3988.

[5]UNSALAN C,BOYER K L.A theoretical and experimental investigation of graph theoretical measures for land development in satellite imagery[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2005,27(4):575-589.

[6]KARATHANASSI V,IOSSIFIDIS C,ROKOS D.A texture-based classif ication method for classifying built areas according to their density[J].International Journal of Remote Sensing,2000,21(9):1807-1823.

[7]BRUZZONE L,CARLINL.A multilevel context-based system for classif ication of very high spatial resolution images[J].IEEE Transactions on Geoscience and Remote Sensing,2006,44(9):2587-2600.

[8]FAUVEL M,CHANUSSOT J,BENEDIKTSSONJ A.Decision fusion for the classif ication of urban remote sensing images[J].IEEE Transactions on Geoscience and Remote Sensing,2006,44(10):2828-2838.

[9]HU Xiangyun,SHENJiajie,SHAN Jie,PAN Li.Local edge distributions for detection of salient structure textures and objects[J].IEEE Geoscience and Remote Sensing Letters.

[10]COMANICIU D.Mean Shift:a robust approach toward feature space analysis[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2002,24(5):1-18.

[11]COMANICIU D,RAMESH V,MEER P.Real-time tracking of non-rigid objects using Mean Shift[C]//IEEE Conference on Computer Vision and Pattern Recognition,2000.

[12]OTSU N.A threshold selection method from graylevel histograms[J].IEEE Transactions on Systems,Man and Cybernetics,1979,SMC-9(1):62-66.

[13]ROSIN P L.A simple method for detecting salient regions[J].Pattern recognition,2009,42(11):263-2371.

[14]SIRMACEK B,UNSALAN C.Built-up area detection using local feature pointsand spatial voting[J].IEEE Geoscience and Remote Sensing Letters,2010,7(1):146-150.