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

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

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

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

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

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

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

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

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

步驟3 計(jì)算每個(gè)像素點(diǎn)到所有邊緣線段的空間距離，并用集合存儲(chǔ)；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 實(shí)驗(yàn)結(jié)果與分析

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

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

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

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

3.2 實(shí)驗(yàn)參數(shù)分析

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

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

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

4 結(jié)語(yǔ)

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

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

[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.