利用機(jī)載LiDAR 點(diǎn)云提取復(fù)雜城市建筑物面域

李 峰，崔希民，袁德寶，王 強(qiáng)，吳亞軍

(1.防災(zāi)科技學(xué)院，三河 065201;2.中國(guó)礦業(yè)大學(xué)(北京)地球科學(xué)與測(cè)繪工程學(xué)院，北京 100083)

0 引言

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

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

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

1 數(shù)據(jù)源與研究區(qū)概況

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

圖1 3個(gè)測(cè)試區(qū)域的位置分布Fig.1 Positions of three test areas

2 地面點(diǎn)云的濾波

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

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

3 建筑物面域的提取

3.1 濾除無(wú)關(guān)點(diǎn)云

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

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

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

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

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

3.3 高差分割植被與建筑物

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

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

3.4 提取建筑物面域

在建筑物和植被面域被分割出來(lái)后，需要判斷各個(gè)面域的歸屬類別。統(tǒng)計(jì)各個(gè)面域的面積，如果小于預(yù)先設(shè)置的最小面積閾值，就刪除此塊面域。另外，引入LiDAR 點(diǎn)云大坡度密度數(shù)來(lái)區(qū)分建筑物面域和植被面域。因?yàn)橐?guī)則建筑物的大坡度值多數(shù)分布在建筑物的邊緣處，而植被點(diǎn)云的高程差異相對(duì)較大，它的大坡度值比建筑物來(lái)更容易向中心擴(kuò)展;因此，LiDAR 點(diǎn)云的大坡度密度數(shù)可以用來(lái)分辨二者。首先，依據(jù)格網(wǎng)高程計(jì)算LiDAR 點(diǎn)云的坡度值g，即

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

然后，統(tǒng)計(jì)各個(gè)面域內(nèi)坡度值大于61°的坡度數(shù)，用坡度數(shù)除以此面域的面積可得大坡度密度數(shù)ρs。對(duì)于ρs小于密度閾值的面域，就識(shí)別為建筑物面域;否則認(rèn)為是植被面域。經(jīng)試驗(yàn)獲知，ρs處于0.5～2.5個(gè)/m2的范圍內(nèi)。

3.5 填充與平滑建筑物面域

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

式中:A 為提取建筑物面域的集合;B 為結(jié)構(gòu)元素;符號(hào)·表示閉運(yùn)算;⊕表示膨脹運(yùn)算;Θ 表示腐蝕運(yùn)算。

4 試驗(yàn)結(jié)果與分析

測(cè)試算法時(shí)所用的3個(gè)區(qū)域參數(shù)和測(cè)試結(jié)果分別見(jiàn)表1 和圖3。

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

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

為了評(píng)價(jià)建筑物提取的精度，使用以下評(píng)價(jià)因子來(lái)估計(jì)提取結(jié)果精度［12］，即

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

結(jié)合點(diǎn)云的高程和反射強(qiáng)度信息，在高密度的點(diǎn)云上人工繪制參考建筑物區(qū)域的線性輪廓，再用本文算法自動(dòng)提取的建筑物面域共同計(jì)算以上質(zhì)量評(píng)價(jià)因子，具體結(jié)果見(jiàn)表2。

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

從表2 看出，3個(gè)區(qū)域的建筑物提取精度結(jié)果中，區(qū)域1 的Q 最好，B 最高，而M 最小;區(qū)域2 的Q 最差，B 一般，但M 最大;區(qū)域3 的Q 一般，B 最小，而M 居中。綜合來(lái)看，區(qū)域1 的提取效果最好，區(qū)域2 的提取效果最差，區(qū)域3 的提取效果一般。區(qū)域1 均存在多余提取和丟失建筑物的情況，區(qū)域2 和區(qū)域3 的多余提取較少而丟失率較大。結(jié)合圖3(b)(c)可以看出，主要是因?yàn)閰^(qū)域2 和區(qū)域3 中均存在未探測(cè)出的建筑物。經(jīng)查看原始點(diǎn)云和影像數(shù)據(jù)后發(fā)現(xiàn)，有2 棟未探測(cè)出的建筑物的屋頂為高吸收率的材料(如瀝青)，所以導(dǎo)致屋頂點(diǎn)云數(shù)量異常稀少，算法無(wú)法進(jìn)行識(shí)別。假設(shè)不存在這2 棟高吸收率的房屋，則區(qū)域2 和區(qū)域3 的M 將會(huì)很小，質(zhì)量將有大的提升。區(qū)域1 的多余探測(cè)率較高，是因?yàn)槠渲写蟛糠址课荼幻艿闹脖痪o密包圍，造成算法無(wú)法準(zhǔn)確區(qū)分植被與建筑物。

5 結(jié)論

在有茂密植被覆蓋區(qū)域的復(fù)雜城市環(huán)境下，本文嘗試探測(cè)并提取了高度在3 m 以上的建筑物面域，選取的3個(gè)典型區(qū)域的提取精度都達(dá)到了91%以上，可以滿足自動(dòng)識(shí)別建筑物的要求，但不能識(shí)別出覆蓋有高吸收率材料的屋頂。此外，本文未考慮高度在3 m 以下的點(diǎn)云少量小建筑物存在的情形。因此，在未來(lái)的工作中，可以結(jié)合高分辨率的多光譜影像來(lái)彌補(bǔ)點(diǎn)云缺失所造成的建筑物丟失的問(wèn)題，測(cè)試算法對(duì)大面積城市區(qū)域的適應(yīng)性，繼續(xù)研究3 m以下建筑物的識(shí)別算法，以求達(dá)到更準(zhǔn)確地提取城市建筑物面域的目標(biāo)。

志謝:感謝德國(guó)DGPF 為本文提供了試驗(yàn)數(shù)據(jù)支持。

［1］Awrangjeb M，Zhang C，F(xiàn)raser 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］成曉倩，樊良新，趙紅強(qiáng).基于圖像分割技術(shù)的城區(qū)機(jī)載Li-DAR 數(shù)據(jù)濾波方法［J］.國(guó)土資源遙感，2012，24(3):29－32.Cheng X Q，F(xiàn)an 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 和航空影像的地震災(zāi)害倒塌建筑物信息提?。跩］.國(guó)土資源遙感，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］李 婷，詹慶明，喻 亮.基于地物特征提取的車載激光點(diǎn)云數(shù)據(jù)分類方法［J］.國(guó)土資源遙感，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.