利用機載LiDAR 點云提取復雜城市建筑物面域

李 峰，崔希民，袁德寶，王 強，吳亞軍

(1.防災科技學院，三河 065201;2.中國礦業大學(北京)地球科學與測繪工程學院，北京 100083)

0 引言

利用高分辨率遙感影像自動提取建筑物是攝影測量和計算機視覺領域的熱點之一。3D 建筑物在城市規劃與管理、緊急事態處置以及災害救援等方面有著重要的應用。由于缺少高程信息，以及受到建筑物陰影和樹木遮閉的影響，利用航空或衛星影像數據提取建筑物信息的精度往往不高。最近新興的機載Li-DAR(light detection and ranging)技術為城市建筑物的識別與量測提供了新方法。加載在飛機平臺上的機載LiDAR 系統可以直接獲取地表、建筑物、樹木和汽車等地物對象的大量3D 點云。與航空或衛星影像相比，機載LiDAR 點云不受陰影和投影差的影響，因此在城市建筑物提取方面有著較大優勢。

從LiDAR 點云中提取建筑物的通常做法是:首先，從非地面點云中分離出地面點云并生成數字高程模型(digital elevation model，DEM);然后，用全部點云生成的數字表面模型(digital surface model，DSM)減去DEM 得到包含非地面點云的歸一化的數字表面模型(normalized DEM，nDSM);最后，從nDSM 中區分建筑物點和植被點。該過程用到的點云和影像信息包括高程、面積、形狀、表面的平滑程度、色彩和紋理［1－2］。Zhang 等［3］聯合區域生長法和平面擬合的方法借助高差從非地面點云中分離出建筑物和樹木;Rottensteiner 等［4］將高差、歸一化差值植被指數(normalized difference vegetation index，NDVI)和表面粗糙度參數應用到證據推理理論中，從而提取建筑物;Khoshelham 等［5］將首次和末次回波DSM 的高差以及從彩紅外影像中提取的NDVI應用到貝葉斯決策函數中進行建筑物分類;Zhou等［6］通過點云分布的規則性、水平性、平整度和法線分布等參數，使用支持向量機(support vector machine，SVM)法探測植被點類，根據區域生長法確定建筑物面片;Meng 等［7］基于形態學的方法識別并恢復建筑物區域，以面積和緊湊性剔除非建筑物區;Salah 等［8］引入灰度共生矩陣的同質性和熵來區別建筑物和植被;于海洋等［9］結合面向對象的方法和SVM 技術提取了地震后復雜環境下倒塌建筑物的信息;李婷等［10］采用點云的高程、法向量和平面投影密度信息從車載激光點云中分離建筑物點。

目前多數算法需要結合高分辨率的多光譜影像才能取得較高的精度，計算方法復雜、效率較低;另外，建筑物的投影差、陰影、影像的空間分辨率、時間差異和配準等造成的融合誤差很難消除，建筑物提取精度受到影響。因此，本文首先直接基于原始機載LiDAR 點云，以漸進式不規則三角網(TIN)加密法過濾出非地面點云，構建非地面點的nDSM;然后通過自定義的形態學分割算子斷開樹木與建筑物之間可能的連接，利用區域生成法獲得建筑物和植被的面域;最后根據點云的大坡度密度性質來快速區分建筑物的面域。

1 數據源與研究區概況

采用的機載LiDAR 點云來源于德國攝影測量、遙感與地理信息協會(DGPF)提供的專門用于3D建筑物重建和城市地物分類的ISPRS 測試工程數據。試驗數據由Leica ALS50 掃描儀獲取，相對航高為500 m，點云平均密度為9.2個/m2，點云包括多回波和反射強度信息。試驗區域范圍是德國Vaihingen 市，包含建筑物、低植被(草地)、樹木、汽車、道路及人工地面等地類，有輕微的地形起伏，為了提高算法的執行效率，本文選擇了3個區域進行測試。區域1、區域2 和區域3 的面積分別為23 552.19 m2，30 145.28 m2和30 285.94 m2;區域1 的建筑物以尖頂房為主，大小各異，形狀復雜，植被較少但與建筑物連接緊密;區域2 的地形起伏較大，有4 座大型平頂樓房，樓房形狀復雜，高度呈階梯狀分布，茂密高大的植被包圍了全部建筑物;區域3 以尖頂房為主，大小不一，植被繁多且以中低高度為主。3個區域的數據都很復雜，適合檢測建筑物面域提取算法的效果。試驗區域分布如圖1 所示。

圖1 3個測試區域的位置分布Fig.1 Positions of three test areas

2 地面點云的濾波

為了更好地識別建筑物，需要從非地面點云中分離出地面點云。首先，由于激光的多路徑效應、飛鳥、空中漂浮物和激光掃描儀的系統誤差可能造成不必要的高位或低位的粗差點云，所以通過生成LiDAR點云的高程直方圖的方法來確定閾值，剔除低位和高位粗差點云的噪聲點云。然后，采取漸進式TIN 加密法對地面點云進行分類［11］。具體算法:將量測的最大建筑物尺寸作為最大格網的邊長，從每個格網中選擇最低的高程點作為種子點生成初始地表TIN，以待判點到三角面片的迭代距離d 和與結點的迭代角度α 作為判斷依據，對于滿足一定的最大距離閾值dmax和最大角度閾值αmax的點添加到TIN 中，組成新的TIN;對于超過dmax和αmax的地形突變點(如懸崖和建筑物邊緣)，將最近的三角形結點作為對稱中心，生成一個虛擬的鏡像點，計算該點到最近三角形面片的距離，如果該距離大于dmax，則剔除該點，否則保留為地面點類并添加到TIN 中。程序迭代運行直到沒有更多的點被添加到TIN 中為止，最后TIN 中的點即為初步分類的地面點類。每次迭代用到的dmax和αmax分別設定為地形高差和地形坡度的中位數。

漸進式TIN 加密法同時受迭代距離和迭代角度的共同約束，因此會拒絕一些迭代角度較大而迭代距離較小的地面點，并認為它們是非地面點。為了分類出所有的地面點，用初步分類的地面點構建DEM，如果點云與DEM 的高程差小于LiDAR 點云的最大高程誤差0.3 m 時，就將該點云歸類為地面點類，最后將新分類的地面點重新生成DEM。

3 建筑物面域的提取

3.1 濾除無關點云

LiDAR 點云的地面點濾波后，剩余的非地面點類包含草地、灌木、樹木、房屋和汽車等。考慮到城市房屋的高度約3 m，則將與DEM 高度差小于3 m的點簡單地歸類為低植被點，剩余的LiDAR點云主要包括樹木和房屋點類。此時，這些點云中存在著少量的孤立點云，它們組成樹木或屋頂面的幾率很小。為了剔除這些孤立點云，遍歷每個點云周圍其他點云的個數，如果點數小于4，則認為是孤立點，予以刪除。

3.2 建筑物與植被的分割算子

此時剩余的未分類點云用來生成二值化格網，只要空格網最鄰近的點和落入格網內的點為地面點、低植被點或者孤立點中的任意一種，則格網被賦值為0;否則賦值為1。同時，對應于同樣大小的格網生成二維高程數組，如果格網為非空格網，取落入格網內的最高點作為格網高程;如果格網為空格網，則取距離格網最近點的高程作為格網高程。

建筑物和樹木生成的二值化格網仍然存在大量建筑物和植被粘連在一起的現象，為進一步分離二者，需要判斷非0 中心格網周圍8 鄰域格網的分布情況。圖2(a)展示的是水平、垂直和對角線值分布均為0－1－0 的情形，當它們兩側的空白格網分別不同時為0 時，則中心格網賦值0。圖2(b)(c)展示的是非0 中心格網周邊分布值為0－1－0－1 的情形，如果其余5個空白格網中有任意一個格網為0，則中心格網賦值0。

圖2 建筑物與植被的分割算子Fig.2 Segmentation operators of building and vegetation

3.3 高差分割植被與建筑物

為了保證建筑物與樹木之間最大程度的分離，參照LiDAR 點云之間的高差，采用圖像處理中的區域生長算法繼續分割房屋與樹木。首先，搜索值為1 的格網作為種子點，將種子點的8 鄰域范圍內點的高程與種子點的高程進行對比，如果小于一定的高差閾值，則將該點包含到種子點區域，新搜索的點重新作為種子點重復以上搜索及合并的過程，直到沒有更多的點可以合并為止。高差閾值dh可由房頂最大坡度smax和點云平均間距c 計算得到，有時點云平均間距小于實際間距，所以擴大2 倍以保證得到更準確的高差;計算的高差再與LiDAR 點云的高程誤差dh0比較，取其中較大者作為高差閾值，即

式中dh0一般在0.15～0.3 m 之間。

3.4 提取建筑物面域

在建筑物和植被面域被分割出來后，需要判斷各個面域的歸屬類別。統計各個面域的面積，如果小于預先設置的最小面積閾值，就刪除此塊面域。另外，引入LiDAR 點云大坡度密度數來區分建筑物面域和植被面域。因為規則建筑物的大坡度值多數分布在建筑物的邊緣處，而植被點云的高程差異相對較大，它的大坡度值比建筑物來更容易向中心擴展;因此，LiDAR 點云的大坡度密度數可以用來分辨二者。首先，依據格網高程計算LiDAR 點云的坡度值g，即

式中:d z/d x 表示坡度在水平方向上的變化率;d z/d y表示坡度在垂直方向上的變化率。

然后，統計各個面域內坡度值大于61°的坡度數，用坡度數除以此面域的面積可得大坡度密度數ρs。對于ρs小于密度閾值的面域，就識別為建筑物面域;否則認為是植被面域。經試驗獲知，ρs處于0.5～2.5個/m2的范圍內。

3.5 填充與平滑建筑物面域

識別出的建筑物面域中仍然存在由采集過程產生的少量數據空白，為了填充這些小空白區域，需要預先量測出最大數據空白的面積，以這個填充面積除以單位格網的面積，取整后可得數學形態學閉運算所需要的窗口(結構元素)B 的大小。數學形態學閉運算具有填充比窗口小的缺口或孔洞，連通小間隔的間斷以及平滑邊界的功能，閉運算還可以部分彌補建筑物面域的缺失。閉運算是用結構元素B對初始提取的建筑物面域A 進行膨脹，緊接著再進行一次腐蝕的結果，其公式為

式中:A 為提取建筑物面域的集合;B 為結構元素;符號·表示閉運算;⊕表示膨脹運算;Θ 表示腐蝕運算。

4 試驗結果與分析

測試算法時所用的3個區域參數和測試結果分別見表1 和圖3。

表1 建筑物面域提取參數Tab.1 Parameters for extracting building regions

圖3 建筑物面域提取結果Fig.3 Results of extracting building regions

為了評價建筑物提取的精度，使用以下評價因子來估計提取結果精度［12］，即

式中:T 表示提取數據與參考數據相匹配的部分;F 表示未與參考數據匹配的提取數據;N 表示未與提取數據匹配的參考數據;C 表示正確提取的結果占參考數據的比率;Q 表示正確提取的結果占全部數據的比率;B 表示未匹配的提取結果占正確提取結果的比率;M 表示未匹配的參考數據占正確提取結果的比率。

結合點云的高程和反射強度信息，在高密度的點云上人工繪制參考建筑物區域的線性輪廓，再用本文算法自動提取的建筑物面域共同計算以上質量評價因子，具體結果見表2。

表2 3個區域建筑物面域提取結果的精度Tab.2 Precision of extracting building regions from three test areas (%)

從表2 看出，3個區域的建筑物提取精度結果中，區域1 的Q 最好，B 最高，而M 最小;區域2 的Q 最差，B 一般，但M 最大;區域3 的Q 一般，B 最小，而M 居中。綜合來看，區域1 的提取效果最好，區域2 的提取效果最差，區域3 的提取效果一般。區域1 均存在多余提取和丟失建筑物的情況，區域2 和區域3 的多余提取較少而丟失率較大。結合圖3(b)(c)可以看出，主要是因為區域2 和區域3 中均存在未探測出的建筑物。經查看原始點云和影像數據后發現，有2 棟未探測出的建筑物的屋頂為高吸收率的材料(如瀝青)，所以導致屋頂點云數量異常稀少，算法無法進行識別。假設不存在這2 棟高吸收率的房屋，則區域2 和區域3 的M 將會很小，質量將有大的提升。區域1 的多余探測率較高，是因為其中大部分房屋被茂密的植被緊密包圍，造成算法無法準確區分植被與建筑物。

5 結論

在有茂密植被覆蓋區域的復雜城市環境下，本文嘗試探測并提取了高度在3 m 以上的建筑物面域，選取的3個典型區域的提取精度都達到了91%以上，可以滿足自動識別建筑物的要求，但不能識別出覆蓋有高吸收率材料的屋頂。此外，本文未考慮高度在3 m 以下的點云少量小建筑物存在的情形。因此，在未來的工作中，可以結合高分辨率的多光譜影像來彌補點云缺失所造成的建筑物丟失的問題，測試算法對大面積城市區域的適應性，繼續研究3 m以下建筑物的識別算法，以求達到更準確地提取城市建筑物面域的目標。

志謝:感謝德國DGPF 為本文提供了試驗數據支持。

［1］Awrangjeb M，Zhang C，Fraser C S.Improved building detection using texture information［C］//International Archives of Photogrammetry，Remote Sensing and Spatial Information Sciences，2011，38(3/W22):143－148.

［2］成曉倩，樊良新，趙紅強.基于圖像分割技術的城區機載Li-DAR 數據濾波方法［J］.國土資源遙感，2012，24(3):29－32.Cheng X Q，Fan L X，Zhao H Q.Filtering of airborne LiDAR data for cityscapes based on segmentation［J］.Remote Sensing for Land and Resources，2012，24(3):29－32.

［3］Zhang K，Yan J，Chen S C.Automatic construction of building footprints from airborne LiDAR data［J］.IEEE Transactions on Geoscience and Remote Sensing，2006，44(9):2523－2533.

［4］Rottensteiner F，Trinder J，Clode S，et al.Building detection by fusion of airborne laser scanner data and multi－spectral images:Performance evaluation and sensitivity analysis［J］.ISPRS Journal of Photogrammetry and Remote Sensing，2007，62(2):135－149.

［5］Khoshelham K，Nedkov S，Nardinocchi C.A comparison of bayesian and evidence－ based fusion methods for automated building detection in aerial data［C］//International Archives of the Photogrammetry，Remote Sensing and Spatial Information Sciences，2008，37(B7):1183－1188.

［6］Zhou Q Y，Neumann U.Fast and extensible building modeling from airborne LiDAR data［C］//Proceedings of the 16th ACM Sigspatial International Conference on Advances in Geographic Information Systems.New York，2008－11－05(07).

［7］Meng X，Wang L，Currit N.Morphology－based building detection from airborne LiDAR data［J］.Photogrammetric Engineering and Remote Sensing，2009，75(4):437－442.

［8］Salah M，Trinder J，Shaker A.Evaluation of the self－ organizingmap classifier for building detection from LiDAR data and multispectral aerial images［J］.Journal of Spatial Science，2009，54(2):1－20.

［9］于海洋，程 鋼，張育民，等.基于LiDAR 和航空影像的地震災害倒塌建筑物信息提取［J］.國土資源遙感，2011，23(3):77－81.Yu H Y，Cheng G，Zhang Y M，et al.The detection of earthquake－ caused collapsed building information from LiDAR data and aerophotograph［J］.Remote Sensing for Land and Resources，2011，23(3):77－81.

［10］李 婷，詹慶明，喻 亮.基于地物特征提取的車載激光點云數據分類方法［J］.國土資源遙感，2012，24(1):17－21.Li T，Zhan Q M，Yu L.A classification method for mobile laser scanning data based on object feature extraction［J］.Remote Sensing for Land and Resources，2012，24(1):17－21.

［11］Axelsson P E.DEM generation from laser scanner data using adaptive TIN models［C］//International Archives of Photogrammetry，Remote Sensing and Spatial Information Sciences，2000，33(B4):110－117.

［12］Hermosilla T，Ruiz L A，Recio J A，et al.Evaluation of automatic building detection approaches combining high resolution images and LiDAR data［J］.Remote Sensing，2011，3(6):1188－1210.