基于MSER的無人機圖像建筑區域提取

丁文銳，康傳波，李紅光，劉碩

(1.北京航空航天大學 無人駕駛飛行器設計研究所，北京100191;2.北京航空航天大學 電子信息工程學院，北京100191)

近年來，隨著技術的不斷發展，無人機廣泛應用在軍事和民用領域.建筑區域對于無人機而言，是一類重要的感興趣目標，一方面對其快速檢測，是無人機完成導航、偵察、監測等任務的基礎和重要內容;另一方面，無人機在出現故障等緊急情況時，通過對建筑區域進行準確檢測并及時規避，從而極大地減少或避免人員傷亡和財產損失.

與現有文獻針對的SAR圖像、高光譜圖像及航空數碼圖像不同，無人機機載CCD攝像機拍攝的中低航空高度圖像，主要具有以下特征:

1)在空間分布方面，無人機可拍攝的地域廣闊;

2)從獲取時間上，不同時間可獲取的無人機圖像會受到不同氣象條件的影響;

3)從建筑區域在圖像中的特性而言，建筑區內部及周圍紋理復雜，部分地區建筑分布不規則，區域整體沒有明顯的方向性，不利于紋理特征的檢測.

目前，建筑區域提取的方法主要分為基于監督學習的方法和基于特征統計的方法.基于監督學習的方法主要是融合多種特征訓練分類器，從而完成圖像區域分類和建筑區域提取.文獻［1］提出一種知識規則與支持向量機(SVM，Support Vector Machine)相結合的面向對象建筑物提取方法;文獻［2］利用多核SVM學習，綜合多種特征，并設計基于圖割的算法提高邊界分割精細程度.這類方法需要人工確定合適的訓練特征和分類類別，并且訓練樣本的選取對分類結果具有很大影響.無人機執行任務的地域廣闊，不同的地域圖像特征差異較大及工況條件復雜，以現有無人機圖像作為樣本訓練得到的分類器，往往不能適應實際應用需求.

基于特征統計的方法，利用多種特征進行統計分析，從而判斷圖像中的建筑區域.文獻［3］利用SIFT特征、多模版匹配和圖論的方法進行檢測，取得一定效果，但是時間復雜度高，并且當房屋屋頂與背景對比度不高以及房屋密度過于集中的情況下，檢測結果不理想;文獻［4］利用地面高程模型(DSM，Digital Surface Model)和光譜信息排除非建筑物，該方法依賴DSM高度等信息，應用范圍狹窄;文獻［5］利用知識規則，依靠幾何、灰度及拓撲特征對建筑物進行描述，并不能適應多種多樣的建筑;文獻［6］利用直線檢測形成屋頂矩形幾何形狀結合亮度信息，綜合紋理特征的SVM噪點去除，但該方法在直線特征復雜及對比度低的場景具有局限性.在基于特征統計的方法中，利用紋理進行建筑區域檢測的方法是國內外研究學者關注的重點.由于Gabor濾波器兼具方向性和紋理分析的優點，研究者提出了多種基于Gabor紋理描述的建筑區域檢測方法.文獻［7］利用Gabor濾波器提取圖像特征，再使用最優決策分類方法，取得了較文獻［3］相似的性能，大大減少了計算時間;文獻［8］利用核密度估計生成居民區置信圖像并計算自適應閾值分割置信圖，對候選區幾何形狀進行選擇，較文獻［7］方法計算效率更高，但由于光照視角等因素影響，某些地質區域會呈現居民區域相似紋理，導致準確度下降.紋理特征提取尺度參數仍需依據經驗值給定，缺乏自適應性，并且無人機圖像中建筑區內部紋理復雜，紋理特征往往不具備明顯主方向，影響檢測效果;同時，在無人機中低空偵察圖像中，道路、草地、水域等紋理對建筑區紋理的提取形成很大干擾.

基于上述分析，利用紋理、監督訓練等方法對無人機中低空偵察圖像進行建筑區域檢測并不能完全滿足要求.注意到，無人機圖像中建筑區域表現出高亮性和團聚性，主要是由于自然背景的光反射能力較弱，使得人造區域的亮度值較大，形成局部高亮，同時會團聚成具有一定面積大小、分布集中的區域.本文結合無人機實際應用需求，從無人機圖像中建筑區域的高亮性和團聚性出發，利用最大穩定極值區域(MSER，Maximum Stable Extremal Regions)檢測建筑物，并根據局部房屋密度去除噪點，采用自適應K均值聚類算法得到聚集的若干建筑區域，最后采用Graham算法生成建筑區邊界，從而實現建筑區域的自動快速提取.

1 最大穩定極值區域提取算法

Matas等［9］所提出的最大穩定極值區域主要基于魯棒性的寬基線立體重建，使用地形中分水嶺的概念來求解穩定局部區域.文獻［10］證明，MSER檢測算子在眾多仿射不變特征區域檢測中，多數情況下具有最佳性能.實驗證明，該方法具有仿射變化不變性、良好的穩定性、計算簡單高效，并且可以檢測不同精細程度的區域.其基本思想如下:選取一系列閾值 t={0，1，…，255}對圖像進行二值分割，低于閾值的像素置為黑色，高于或等于閾值的像素置為白色.在閾值t從0增大到255的過程中，圖像中會形成閉合區域，其中在一定閾值范圍內面積變化最小的區域即為最大極值穩定區域.

最大穩定極值穩定區域提取算法步驟包括以下幾個步驟［11］.

1.1 MSER 檢測

1)對給定的圖像根據灰度值排序，如果是彩色圖像，需要將彩色圖像轉換成灰度圖像

2)按照降序或升序將排序后的像素放入圖像中，鏈接區域的列表和面積使用高效的合并-查找(UnionFind)［12］算法來維護.

3)Qi表示二值化閾值對應的二值圖像中的任意連通區域，當二值化閾值在［i－Δ，i+Δ］變化時，連通區域也相應地變為 Qi+Δ和 Qi－Δ.在這個變化范圍內具有極小變化率q(i)區域就被認為是MSER.

為了能同時得到最小灰度(即黑色)和最大灰度(即白色)的MSER，在求得最小灰度MSER+后，需要將原始圖像進行反轉:Ir=Imax－I，再利用反轉后的圖像求取最大灰度的MSER－.

1.2 MSER區域擬合及歸一化

完成圖像中MSER的區域檢測后，不規則區域不便于進行歸一化和提取特征描述，必須對它們進行橢圓擬合.為了能夠構造具有仿射不變性的特征描述子，首先將MSER對應的橢圓擬合區域擴大成抽取特征用的橢圓特征測量區，測量區與擬合區同中心，一般測量區尺度是擬合區尺度的3倍.然后將測量區歸一化為指定大小的區域(成為歸一化區)，再在歸一化區的圖像上提取特征描述子.區域歸一化的方法主要包括3步［13］:

1)將MSER擬合區對應的測量區進行仿射歸一化，得到一個半徑指定大小的歸一化區，歸一化區的半徑尺寸為20～30像素.

2)在歸一化區內進行圖像梯度直方圖統計，其目的是找出該直方圖的最大值，并將該最大值對應的方向作為歸一化區圖像梯度的主方向.

3)根據主方向對歸一化區圖像再進行旋轉歸一化，旋轉后圖像的梯度主方向為零度.

2 建筑區域提取算法

2.1 方法框架和流程

如圖1所示，本方法的總體思路為:首先對無人機航拍圖像進行預處理，去除模糊，然后利用MSER算法提取建筑房屋，進而計算每個房屋中心點鄰域內的房屋密度，房屋密度較低的區域被認為是干擾點予以剔除，完成點的選擇與濾波后，再利用自適應K均值聚類算法，將離散點聚成若干建筑區，最后采用Graham算法，生成每個建筑區的包絡，完成建筑區的自動提取與標記，并在最后輸出區域的數量K，若為K=0，則表示不存在建筑區.

圖1 算法流程圖Fig.1 Algorithm flow chart

2.2 算法原理

2.2.1 無人機圖像預處理

無人機圖像容易受不同氣象條件的影響，尤其在霧天等不良天氣狀況下拍攝的圖像模糊不清，圖像質量差，導致提取的特征數量較少，給建筑區域檢測帶來影響.為了使算法既簡潔又高效，采用直方圖均衡和GAMMA調節對圖像進行預處理，擴展圖像的動態范圍，同時增強圖像的亮度，實驗表明經過預處理之后，檢測的特征數量明顯增加，進一步提高建筑區的可分辨性.預處理前后的結果對比如圖2和圖3所示.

圖2 原始圖像Fig.2 Initial image

圖3 預處理后的圖像Fig.3 Preprocessed image

2.2.2 最大穩定極值區域檢測

屋頂等人造物體，光反射能力比自然背景更加強烈，呈現團塊高亮特點，如圖4所示，每個圖對應的下方圖像為其二值化圖像，對比度高，并且在一定閾值范圍內面積變化最小，符合MSER特征.

圖4 單個房屋特征Fig.4 Features of single house

利用第1節所述最大穩定極值區域提取算法，計算并標記圖像穩定區域，數量為1 358個區域，圖5中橢圓標注的區域為穩定區域MSER+，圖6中橢圓標注的區域為代表圖像反轉后得到的穩定區域MSER－.

圖5 MSER+區域Fig.5 MSER+area

圖6 MSER-區域Fig.6 MSER-area

為了確定圖像中建筑房屋的大小，還需知道圖像的像素地面分辨率，在垂直下視的情況下，無人機高度為H(單位:m)，拍攝的圖像大小為w×h(單位:像素)，攝像機的視角為水平 α×垂直β(單位:°)，計算方法如式(1)～式(4).選取實飛圖像數據，飛行高度為5 km，攝像機視場角為14.38°× 10.59°，圖像大小為1392像素×1040像素，像素地面分辨率計算如圖7所示.水平距離 =2×5 000×tan(7.19)=1 261.15 m=1392像素，因為1個像素≈0.906m，同樣，垂直方向:1個像素≈0.891 m.

圖7 像素地面分辨率計算Fig.7 Calculation of ground pixel resolution

式中，Lhor，Lver分別表示水平和垂直距離;Phor，Pver分別表示水平和垂直方向上單位像素的實際地理距離.

郊區或農村房屋的面積S一般在30～200 m2，根據圖7計算的像素地面分辨率結果可知，在圖像中一間房屋長度和寬度所占像素個數的范圍分別在［8，20］和［5，12］之間，以此限定MSER檢測的穩定區域大小，檢測結果如圖8所示.

圖8 建筑區域檢測Fig.8 Building area detection

2.2.3 基于區域密度的建筑區篩選方法

在實際圖像中，MSER所檢測得到的最大穩定極值區域并不都是建筑區，如圖7所示，田野中也存在部分最大穩定極值區域，形成干擾.鑒于通過平滑濾波或無重疊框遍歷圖像的方法來排除干擾不夠準確，本方法利用建筑區特征豐富密集的特點，通過計算最大穩定極值區域密度篩選出建筑區.

確定MSER檢測的最大穩定極值區域的中心，并用來代替最大穩定極值區域計算區域密度.如圖9所示，黑點代表特征區域，白色為背景.采用m×m大小的窗口，依次計算每個特征位置(x，y)的區域密度 f(x，y).

圖9 建筑區檢測點Fig.9 Detected points of building area

圖10 濾波去除噪點Fig.10 Filter to remove noise

2.2.4 自適應K均值聚類

因為聚集性和分群而居的特點，同一幅圖像中可能出現多個聚集的建筑區域，為了準確標記建筑區，需要對上述特征區域中心點進行聚類.K均值聚類算法［14］處理這種聚類問題簡單高效，但是無人機圖像中建筑區域的個數并不確定，因此無法自動確定算法的參數.由于無人機視場較小，在單幅無人機圖像中，出現的建筑區域個數并不會太多，因此本文在K均值算法的基礎上提出了一種自適應策略，給定K一個較大的初始值K=6，從而達到自動確定建筑區數量的目的.具體流程如圖11所示.

圖11 自適應K均值聚類計算流程圖Fig.11 Flow chart of adaptive K-means clustering calculation

限于圖像實際地理范圍，給定初始K值為6，自適應聚類步驟如下:

1)利用 K均值聚類算法，對離散點進行聚類;

2)分別計算K個區域的中心點;

3)K個區域兩兩組合，如圖12所示，計算區域間最短距離CD;

4)若距離CD＞V或者K=1，則認為全部都是獨立區域，否則K值減1，繼續循環計算，閾值V的選取根據像素所代表的實際距離計算，V=，按照圖7條件，V值取150.

基于此方法，對圖10處理結果為K=1，存在一個建筑區.

圖12 區域間距離計算方法Fig.12 Calculation methods of distance between areas

2.2.5 生成建筑區包絡

根據2.2.4節得到的K個離散區域，采用經典的Graham算法［15］生成建筑區凸殼包絡，該算法為最優凸殼算法，時間復雜度約為O(nlogn).

Graham算法的基本思想是:以y坐標最小點為初始點P1，計算該點與其他點連線的水平夾角，按夾角大小和到初始點P1的距離進行排序，得到一系列點集合 P={P1，P2，…，Pn}，其中Pn+1=P1，對點集合{P1，P2，…，Pn，Pn+1}依據凸多邊形的各頂點必在該多邊形的任意一條邊的同一側的方法原理，刪除不屬于凸殼上的點，得到新的點集合P′={P1，P2，…，Pm}，最后根據點集合 P′標記邊界，實現建筑區的自動提取，效果如圖13所示.

圖13 建筑區包絡圖Fig.13 Envelop of building area

3 實驗結果與分析

對部分圖像人工標注出建筑區域作為基準數據，并將提取結果與人工標注結果相比較.圖14展示了不同測試數據的提取結果.

圖14 測試結果Fig.14 Experiment results

從圖14直觀上看，人工標記的結果與自動提取的結果大體上比較接近，由于Graham算法是凸包算法，只能生成區域凸殼，因此人工標注結果與自動提取結果有部分位置的偏差.為了更好地描述算法效果，采用無人機圖像數據集，對自動提取的像素數目做統計，并以人工提取的像素數目作為標準，統計準確率、誤檢率和漏檢率及時間性能.本文進行兩個方面的測試:一方面是MSER性能分析;另一方面是分析圖像預處理和自適應K均值聚類方法對結果的影響.

本方法以無人機中低空偵察圖像集為對象，無人機高度為5 km，攝像機視場角為14.38×10.59°，圖像大小為 1 392像素 ×1 040像素，實驗計算機配置為Intel Core 2 Duo雙核處理器，主頻2.20 GHz，內存 2.00 GB，平臺為 Microsoft Visual C++6.0.

3.1 MSER性能分析

基于特征統計的建筑區域提取方法中，利用DSM、光譜等信息有自身應用的局限性，而根據SIFT、紋理特征有更為廣泛的適用性.本部分討論MSER算法相比基于Gabor變換的紋理特征、SIFT特征點的提取算法的時間性能分析，Gabor濾波器選擇4方向，設定濾波窗口大小為尺度為5×5.

時間性能方面，如表1所示，針對大小為1392像素 ×1 040像素的圖像集做特征檢測，MSER算法平均檢測1 358個特征區域耗時1.092 s，SIFT算法平均檢測1 175個特征點耗時3.073 s，紋理特征檢測算法平均耗時 2.509 s.因此，在相同特征點數量的情況下，MSER算法時間性能優于SIFT，在圖像大小相同的情況下，MSER算法優于紋理檢測方法，計算簡單高效.

表1 特征檢測方法分析Table1 Analysis of feature detection methods

3.2 實驗結果分析

表2為本算法提取結果，采取4組圖像集測試，共3722個建筑數量，平均正確率達到92.25%，誤檢率也較低.同時，分析對比有無圖像預處理和自適應K均值聚類兩部分的算法結果，如表3所示，實驗1為本文算法結果，實驗2為沒有圖像預處理的提取結果，實驗3代表了沒有自適應K均值聚類的提取結果.實驗1和2對比，在沒有預處理的情況下，雖然時間性能得到提高僅有0.06s，但正確率下降了3.5%同時造成誤檢率上升;實驗1和3對比，自適應K均值聚類對正確率和漏檢率基本無影響，但有效降低了誤檢率，非自適應K均值聚類的情況下，雖然運算時間減少了 0.02s(約 1.69%)，但誤檢率高.

表2 算法提取結果分析Table2 Analysis of algorithm extraction results

表3 提取結果對比分析Table3 Comparison analysis of extraction results

因此，本方法在保證時間性能的前提下取得了較高的正確檢測率，大多數圖像中的建筑區域能夠被檢出.另外，也會造成區域漏檢和誤檢，例如只檢測出部分區域的情況，造成區域的漏檢，以及檢測出不屬于建筑區的區域，造成誤檢，但漏檢率和誤檢率比較低.這種情況的發生，主要是因為若干建筑區域面積較小不易被檢出，或者居民區距離較近，此時居民區之間的空地被連入建筑區范圍內，造成誤檢，誤檢率在1.1%左右.文獻［8-9］中使用紋理檢測原理檢測建筑區，使用中低分辨率圖像，圖像大小為235像素×265像素，算法耗時分別為1.99 s和0.42 s，本文算法使用1392像素×1 040像素大小圖像耗時1.18 s.綜上所述，相比SIFT和紋理檢測方法，本算法耗時較短，處理效率較高.

4 結論

通過分析計算，現得出結論如下:

1)最大穩定極值區域提取算法具有良好的穩定性，計算簡單高效;

2)將MSER檢測結果與區域大小、密度計分布特征相結合，采用自適應K均值聚類的方法可穩定地提取建筑區域.

實驗數據表明，該方法提取精度為92.25%，具有較好的準確率，有效可行.但是，該算法采用Graham算法生成的凸殼邊界較為簡單，這也是造成誤檢率很重要的一個原因;另外，在森林地區由于樹木之間的遮擋造成局部存在穩定區域，對建筑區域檢測提取造成誤擾.如何改進邊界生成的方法和增強本文算法魯棒性是日后工作的重點和方向.

References)

［1］ 劉海飛，常慶瑞，李粉玲.高分辨率影像城區建筑物提取研究［J］.西北農林科技大學學報:自然科學版，2013，41(10):221-227.Liu H F，Chang Q R，Li F L.Urban building extraction of highresolution images［J］.Northwest Agriculture and Forestry University of Science and Technology:Natural Science Edition，2013，41(10):221-227(in Chinese).

［2］ Tao C，Tan Y H，Yu J G，et al.Urban area detection using multiple Kernel Learning and graph cut［C］//2012 32nd IEEE International Geoscience and Remote Sensing Symposium.Piscataway，NJ:IEEE，2012:83-86.

［3］ Sirmacek B，Unsalan C.Urban-area and building detection using SIFT keypoints and graph theory［J］.Geoscience and Remote Sensing，IEEE Transactions on，2009，47(4):1156-1167.

［4］ 張立民，張建廷，徐濤.基于對象的最優尺度建筑物信息提取方法［J］.計算機應用研究，2012，29(12):4789-4792.Zhang L M，Zhang J T，Xu T.Object-based of optimal scale building information extraction method［J］.Application Research of Computers，2012，29(12):4789-4792(in Chinese).

［5］ 黃金庫，馮險峰，徐秀莉，等.基于知識規則構建和形態學修復的建筑物提取研究［J］.地理與地理信息科學，2011，27(4):28-31.Huang J K，Feng X F，Xu X L，et al.Building extraction based on knowledge rule and morphological restoration［J］.Geography and Geo-Information Science，2011，27(4):28-31(in Chinese).

［6］ 楊萍，姜志國，劉濱濤.一種遙感圖像建筑物檢測新方法［J］.航天返回與遙感，2013，34(5):70-77.Yang P，Jiang Z G，Liu B T.A new building detection method in remote sensing image［J］.Aerospace ＆ Remote Sensing，2013，34(5):70-77(in Chinese).

［7］ Sirmacek B，Unsalan C.Urban area detection using local feature points and spatial voting［J］.IEEE Geoscience and Remote Sensing Letters，2010，7(1):146-150.

［8］ 谷多玉，郭江，李書曉，等.基于Gabor濾波器的航空圖像居民區域提?。跩］.北京航空航天大學學報，2012，38(1):106-110.Gu D Y，Guo J，Li S X，et al.Resident region extraction using Gabor filter in aerial imagery［J］.Journal of Beijing University of Aeronautics and Astronautics，2012，38(1):106-110(in Chinese).

［9］ Matas J，Chum O，Urban M，et al.Robust wide-baseline stereo from maximally stable extremal regions［J］.Image and Vision Computing，2004，22(10):761-767.

［10］ Mikolajczyk K，Tuytelaars T，Schmid C，et al.A comparison of affine region detectors［J］.International Journal of Computer Vision，2005，65(1-2):43-72.

［11］ 王永明，王貴錦.圖像局部不變性特征與描述［M］.北京:國防工業出版社，2010:102-127.Wang Y M，Wang G J.Image local invariant features and description［M］.Beijing:Defense Industry Press，2010:102-127(in Chinese).

［12］ Murphy-Chutorian E，Trivedi M.N-tree disjoint-set forests for maximally stable extremal regions［C］//2006 17th British Machine Vision Conference.Edinburgh，United Kingdom:British Machine Vision Association，2006:739-748.

［13］ Mikolajczyk K，Schmid C.A performance evaluation of local descriptors［J］.Pattern Analysis and Machine Intelligence，IEEE Transactions on，2005，27(10):1615-1630.

［14］ MacQueen J.Some methods for classification and analysis of multivariate observations［C］//Proceedings of the fifth Berkeley Symposium on Mathematical Statistics and Probability.Berkeley，California:University of California Press，1967:281-297.

［15］ Graham R L.An efficient algorithm for determining the convex hull of a finite planar set［J］.Information Processing Letters，1972，1(4):132-133.