朱依民，田林亞，畢繼鑫，林 松

(1. 河海大學(xué)地球科學(xué)與工程學(xué)院，江蘇 南京 211100； 2. 浙江華東測(cè)繪與安全技術(shù)有限公司，浙江 杭州 310014)

建筑物輪廓線精準(zhǔn)提取作為建筑物三維模型重建的關(guān)鍵環(huán)節(jié)，對(duì)建立數(shù)字城市、智慧城市具有重要現(xiàn)實(shí)意義[1]。近年來，機(jī)載激光雷達(dá)測(cè)量技術(shù)得到了較快發(fā)展，與無人機(jī)傾斜影像立體匹配技術(shù)相比，機(jī)載LiDAR測(cè)量能快速采集包含大量城市建筑物三維坐標(biāo)信息的點(diǎn)云數(shù)據(jù)，特別是獨(dú)特的掃描式測(cè)量方法使其在獲取豐富建筑物屋頂信息方面有著得天獨(dú)厚的優(yōu)勢(shì)[2]，為快速、精準(zhǔn)、自動(dòng)提取建筑物輪廓提供了可能。

利用機(jī)載LiDAR點(diǎn)云數(shù)據(jù)自動(dòng)提取建筑物輪廓雖受到了一定關(guān)注，但相關(guān)研究尚處于起步階段。文獻(xiàn)[3]通過規(guī)則化數(shù)字表面模型分類地面點(diǎn)和非地面點(diǎn)，基于8鄰域搜索非地面點(diǎn)云中的建筑物信息，并對(duì)獲取的建筑物表面點(diǎn)云采用梯度圖邊界跟蹤法提取建筑物輪廓信息，但該方法建筑物輪廓提取精度受nDSM分類結(jié)果的影響較大。文獻(xiàn)[4]集成邊緣與局部信息的活動(dòng)輪廓模型，將多波段圖像邊界提取算法用于LiDAR點(diǎn)云數(shù)據(jù)中的建筑物輪廓提取，該方法提取的建筑物屋頂輪廓精度受限于點(diǎn)云分類精度，且分類誤差需要人工判斷與修改。文獻(xiàn)[5]基于高差偏度平衡濾波算法和多級(jí)最小外接矩形算法提取建筑物輪廓，但濾波算法的窗口大小和種子點(diǎn)個(gè)數(shù)等閥值需要通過反復(fù)試驗(yàn)方可確定。文獻(xiàn)[6]結(jié)合建筑物點(diǎn)云與配準(zhǔn)后的影像提取建筑物輪廓，先利用α-shapes算法和投票機(jī)制提取建筑物粗邊界，然后使用從影像中提取的建筑物邊緣信息對(duì)提取的粗輪廓進(jìn)行修正，該方法涉及多種閾值的設(shè)定，且自適應(yīng)程度不高。文獻(xiàn)[7]基于RANSAC算法計(jì)算建筑物平面最佳模型參數(shù)，進(jìn)而構(gòu)建建筑物屋頂面并從中提取其輪廓，但該方法需先進(jìn)行預(yù)處理得到高精度的DSM點(diǎn)云數(shù)據(jù)，且建筑物的屋頂需是平面，不具備普遍適用性。本文針對(duì)以上研究存在的諸多問題，對(duì)機(jī)載LiDAR掃描的點(diǎn)云數(shù)據(jù)特征進(jìn)行分析，采用改進(jìn)的區(qū)域生長(zhǎng)算法對(duì)原始點(diǎn)云數(shù)據(jù)中的地面點(diǎn)、植被點(diǎn)和建筑物點(diǎn)進(jìn)行分類，進(jìn)一步基于三維Hough變換算法提取機(jī)載LiDAR建筑物平面點(diǎn)云，最終采用α-shape算法獲取建筑物的輪廓信息，形成一整套從機(jī)載LiDAR掃描數(shù)據(jù)提取建筑物輪廓信息的數(shù)據(jù)處理體系，為機(jī)載LiDAR在數(shù)字城市、智慧城市建立中的應(yīng)用提供技術(shù)支撐。

1 基于機(jī)載LiDAR掃描數(shù)據(jù)的建筑物輪廓提取方法

1.1 改進(jìn)的區(qū)域生長(zhǎng)算法濾波地面點(diǎn)

飛行器在實(shí)際掃描時(shí)會(huì)受到諸如多路徑效應(yīng)、飛鳥等影響而產(chǎn)生一些噪聲點(diǎn)[8]，為了避免噪聲點(diǎn)對(duì)后期點(diǎn)云數(shù)據(jù)處理產(chǎn)生影響，需先對(duì)原始點(diǎn)云數(shù)據(jù)進(jìn)行去噪處理。常用的去噪方法有局部平面擬合法[9-10]、頻率域法[11]、高程分布直方圖法[12]等。由于噪聲點(diǎn)一般較少且在高程空間分布上較孤立，因此本文建議實(shí)際應(yīng)用時(shí)選擇高程分布直方圖法對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行去噪，從而得到較好的去噪效果。

對(duì)原始LiDAR點(diǎn)云數(shù)據(jù)去噪后，根據(jù)由低到高的處理原則對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行地面點(diǎn)逐步濾波，地面點(diǎn)常用濾波方法有數(shù)學(xué)形態(tài)學(xué)濾波算法、基于不規(guī)則三角網(wǎng)逐漸加密濾波算法及區(qū)域生長(zhǎng)的聚類分割算法。數(shù)學(xué)形態(tài)學(xué)濾波算法需要通過測(cè)區(qū)內(nèi)建筑物尺寸確定形態(tài)學(xué)窗口尺寸，基于不規(guī)則三角網(wǎng)逐漸加密濾波算法需要設(shè)置最大建筑物尺寸、最大地形變化角度閾值、角度和距離判別閾值，且需對(duì)非地面點(diǎn)進(jìn)行多次濾波。顧及以上兩種濾波算法需要設(shè)置諸多參數(shù)的問題，本文對(duì)區(qū)域生長(zhǎng)聚類分割算法進(jìn)行改進(jìn)，在生長(zhǎng)地面點(diǎn)過程只需設(shè)定生長(zhǎng)的高度和坡度，基于改進(jìn)的區(qū)域生長(zhǎng)聚類分割算法對(duì)地面點(diǎn)進(jìn)行濾波的具體步驟如下：

(1) 如圖1所示，遍歷點(diǎn)云數(shù)據(jù)，找出高程最小的3個(gè)點(diǎn)P1、P2和P3作為種子點(diǎn)。

(2) 計(jì)算種子點(diǎn)的平均高程P0，并通過3個(gè)種子點(diǎn)生成種子平面。

(3) 利用KD樹搜索P0點(diǎn)周圍一定數(shù)量的點(diǎn)，并從中篩選以圓心P0、半徑為r的球形鄰域內(nèi)的點(diǎn)。

(4) 計(jì)算球形鄰域內(nèi)所有點(diǎn)到種子平面的垂直距離ρi，以及每個(gè)點(diǎn)和3個(gè)種子點(diǎn)連線與種子平面構(gòu)成的最大夾角βi。

(5) 將計(jì)算的球形鄰域內(nèi)每個(gè)點(diǎn)的ρi、βi與設(shè)定的距離閾值ρ0、角度閾值β0作比較，將同時(shí)小于兩個(gè)閾值的點(diǎn)放入候選種子點(diǎn)集合中。

(6) 為了繼續(xù)生長(zhǎng)，將初始種子點(diǎn)中任意點(diǎn)歸為地面點(diǎn)(此處設(shè)置為P1)，儲(chǔ)存到地面點(diǎn)集合中，從候選種子點(diǎn)集合中選出最小值ρi所對(duì)應(yīng)的點(diǎn)作為新的種子點(diǎn)P4。

(7) 以P2、P3和P4作為種子點(diǎn)重復(fù)步驟(2)—(6)，直至所有點(diǎn)均有所判斷，最終得到的地面點(diǎn)集合即為通過區(qū)域生長(zhǎng)算法生長(zhǎng)得到的地面點(diǎn)。

采用區(qū)域生長(zhǎng)算法濾波地面點(diǎn)后，剩余點(diǎn)云數(shù)據(jù)中包含的低矮植被與車輛，以及較高的植被和建筑物等非地面點(diǎn)，為了在后期三維Hough變換提取建筑物平面過程減少噪聲點(diǎn)對(duì)提取結(jié)果的影響，以濾波的地面點(diǎn)為基準(zhǔn)，將非地面點(diǎn)的高程進(jìn)行歸一化處理，再通過一定的高度閾值將低矮植被點(diǎn)和車輛等非地面點(diǎn)濾除，經(jīng)過歸一化高度閾值分類后的點(diǎn)云數(shù)據(jù)將只包含高植被點(diǎn)和建筑物點(diǎn)。

1.2 三維Hough變換提取建筑物平面

當(dāng)將一條直線用參數(shù)方程r=xcosθ+ysinθ表示時(shí)，可以建立一個(gè)二維Hough空間，同理將平面用參數(shù)方程r=xsinθcosφ+ysinθsinφ+zcosθ表示，便可建立一個(gè)三維Hough空間(r，θ，φ)。其中，r為原點(diǎn)到平面的距離，θ∈(0，π)為平面法向量方向的天頂角，φ∈(-π，π)為平面法向量的坐標(biāo)方位角。

本文將三維Hough變換算法用于機(jī)載LiDAR建筑物平面提取，將所有點(diǎn)云數(shù)據(jù)視為待檢測(cè)的一系列空間點(diǎn)，以(r，θ，φ)參數(shù)的取值間隔將三維Hough空間分割為緊密相連的三維格網(wǎng)，以一定的θ和φ的取值間隔遍歷每個(gè)點(diǎn)云計(jì)算其對(duì)應(yīng)的r值，當(dāng)計(jì)算的r值落入某個(gè)格子內(nèi)便使該格子的計(jì)數(shù)加1。對(duì)所有待檢測(cè)點(diǎn)云進(jìn)行三維Hough變換后，取計(jì)數(shù)最大的格子所對(duì)應(yīng)的(r，θ，φ)并計(jì)算出相應(yīng)的平面，此平面便為提取的建筑物屋頂平面。

1.3 基于α-shape算法提取建筑物輪廓

α-shape算法是指通過將一個(gè)半徑固定的圓繞著點(diǎn)集滾動(dòng)進(jìn)而構(gòu)建點(diǎn)集輪廓的邊緣點(diǎn)檢測(cè)算法[13-14]，目前該算法已被用于網(wǎng)格生成、醫(yī)學(xué)圖像分析等領(lǐng)域。如圖2所示，本文將α-shape算法用于建筑物屋頂面輪廓的邊緣提取，并將圓半徑α設(shè)置為三維Hough變換點(diǎn)集中各點(diǎn)平均間距的2～3倍，滾動(dòng)圓的圓心坐標(biāo)為

(1)

(2)

式中，(Ox，Oy)為滾動(dòng)點(diǎn)的圓心坐標(biāo)；(x1，y1)與(x2，y2)分別為待判斷點(diǎn)P1與P2的平面坐標(biāo)。

基于α-shape算法提取建筑物輪廓的具體步驟如下：

(1) 將點(diǎn)云數(shù)據(jù)集合S內(nèi)的所有點(diǎn)投影到XY平面，從中任取一點(diǎn)P1作為圓心，以2α為半徑，搜索圓內(nèi)所有點(diǎn)并構(gòu)成新的點(diǎn)集S1，從S1中任取一點(diǎn)P2，由P1、P2構(gòu)成平面圓，計(jì)算平面圓的圓心P0。

(2) 計(jì)算點(diǎn)集S1內(nèi)其他點(diǎn)到點(diǎn)P0的距離，當(dāng)所有點(diǎn)的距離都大于α?xí)r，將P1、P2定義為輪廓點(diǎn)，并轉(zhuǎn)至步驟(4)；當(dāng)出現(xiàn)點(diǎn)的距離小于α的情況時(shí)，則轉(zhuǎn)至步驟(3)。

(3) 對(duì)點(diǎn)集S1內(nèi)下一點(diǎn)重復(fù)步驟(1)—(2)，直至S1內(nèi)所有點(diǎn)均有所判斷方可結(jié)束。

(4) 對(duì)點(diǎn)集S內(nèi)未被判斷的點(diǎn)重復(fù)步驟(1)—(3)，直至S內(nèi)所有點(diǎn)均有所判斷方可結(jié)束。

2 實(shí)例計(jì)算

為了驗(yàn)證本文從機(jī)載LiDAR掃描數(shù)據(jù)中提取建筑物輪廓方法的可行性，選取某測(cè)區(qū)的機(jī)載雷達(dá)數(shù)據(jù)進(jìn)行試驗(yàn)，測(cè)區(qū)使用RIEGLVQ-1560i機(jī)載激光雷達(dá)測(cè)量系統(tǒng)進(jìn)行測(cè)量，整個(gè)測(cè)區(qū)面積約為549 388 m2，共掃描點(diǎn)云3 823 127個(gè)。限于篇幅，本文對(duì)整個(gè)測(cè)區(qū)進(jìn)行分塊處理，選取其中一塊面積為16 941.9 m2、點(diǎn)云數(shù)量為190 338個(gè)、點(diǎn)云密度為11.23個(gè)/m2的分塊測(cè)區(qū)進(jìn)行試驗(yàn)。測(cè)區(qū)高低起伏復(fù)雜，含有地面、道路、低植被、中等植被、高大植被、建筑物等豐富的地物，完全滿足試驗(yàn)需要。經(jīng)去噪后的測(cè)區(qū)點(diǎn)云三維顯示如圖3所示。

利用改進(jìn)的區(qū)域生長(zhǎng)算法對(duì)去噪后的點(diǎn)云數(shù)據(jù)進(jìn)行地面點(diǎn)生長(zhǎng)，根據(jù)測(cè)區(qū)實(shí)際情況，設(shè)置KD樹搜索點(diǎn)的數(shù)量為500個(gè)，生長(zhǎng)的距離閾值ρ0為0.5 m，角度閾值β0為15°。經(jīng)過區(qū)域生長(zhǎng)后得到的地面點(diǎn)和非地面點(diǎn)如圖4所示。

基于改進(jìn)的區(qū)域生長(zhǎng)算法對(duì)測(cè)區(qū)地面點(diǎn)進(jìn)行濾波后，得到的非地面點(diǎn)包含低矮植被、車輛，以及高大的植被和建筑物點(diǎn)云。為了盡量減少非地面點(diǎn)中的非建筑物點(diǎn)對(duì)后期建筑物平面提取的影響，本文以生長(zhǎng)得到的地面點(diǎn)為基準(zhǔn)，將非地面點(diǎn)的高程進(jìn)行歸一化處理，設(shè)置3 m的高度閾值將低矮植被和車輛等非地面點(diǎn)去除，此時(shí)非地面點(diǎn)由高植被點(diǎn)云和建筑物點(diǎn)云構(gòu)成。經(jīng)過高度閾值濾波處理后的非地面點(diǎn)如圖5所示。

利用三維Hough變換算法提取建筑物平面，獲取計(jì)算平面時(shí)所需的參數(shù)(r，θ，φ)，再計(jì)算各點(diǎn)到平面的距離，并將距離小于閾值所對(duì)應(yīng)的點(diǎn)歸于同一平面。需要注意的是該步驟會(huì)將某些高于建筑物的高植被點(diǎn)歸類至所提取的平面點(diǎn)集中，本文通過判斷局部鄰域內(nèi)點(diǎn)的法向量方向的一致性去除三維Hough變換提取平面后所含的高植被噪聲點(diǎn)。

經(jīng)過Hough變換提取的平面中會(huì)包含建筑物立面，但在后期提取建筑物輪廓的過程要將建筑物點(diǎn)云投影至XY平面，此時(shí)建筑物的立面點(diǎn)會(huì)被建筑物屋頂面點(diǎn)云所覆蓋。因此，為了提高后期提取建筑物輪廓算法的運(yùn)算效率，本文通過點(diǎn)的法向量方向信息將建筑物立面點(diǎn)去除，提取的建筑物屋頂面點(diǎn)云數(shù)據(jù)如圖6所示。

將提取的建筑物平面點(diǎn)云數(shù)據(jù)投影至XY平面，再利用α-shape算法提取建筑物輪廓點(diǎn)云。設(shè)置滾動(dòng)圓半徑α為0.5 m得到的建筑物輪廓點(diǎn)云和輪廓線如圖7所示。

3 精度評(píng)定

為了對(duì)建筑物輪廓提取精度進(jìn)行定量評(píng)價(jià)，參考文獻(xiàn)[15]并結(jié)合實(shí)際情況,本文提出匹配度、形狀相似度及位置精度共3個(gè)精度評(píng)價(jià)指標(biāo)來評(píng)定建筑物輪廓的提取精度。

3.1 匹配度

本文采用圖像分類精度評(píng)估參數(shù)中的質(zhì)量因子Q作為匹配度的評(píng)價(jià)指標(biāo)，即

(3)

式中，A為本文算法提取的位于原始建筑物輪廓點(diǎn)集中的點(diǎn)個(gè)數(shù)；B為屬于本文算法提取出的輪廓點(diǎn)但不屬于原始的建筑物輪廓點(diǎn)的個(gè)數(shù)；C為屬于原始的建筑物輪廓點(diǎn)但不屬于本文算法提取出的輪廓點(diǎn)的個(gè)數(shù)。

3.2 形狀相似度

通過定義面積差A(yù)d和周長(zhǎng)差Pd描述提取的建筑物形狀相似度,即

(4)

式中，Se為本文算法提取的建筑物面積；Sr為原始的建筑物面積；Ce為本文算法提取的建筑物周長(zhǎng)；Cr為原始的建筑物周長(zhǎng)。

3.3 位置精度

選用建筑物輪廓拐點(diǎn)平面坐標(biāo)的平均距離差Mdd對(duì)位置精度進(jìn)行評(píng)價(jià)，即

(5)

式中，n為每個(gè)建筑物拐點(diǎn)的數(shù)量；(xei，yei)為提取的建筑物拐點(diǎn)坐標(biāo)；(xri，yri)為參考建筑物拐點(diǎn)的坐標(biāo)。

對(duì)整個(gè)測(cè)區(qū)實(shí)例(如圖8(a)所示)所有建筑物采用上述方法進(jìn)行建筑物輪廓提取，提取結(jié)果如圖8(b)所示，同時(shí)計(jì)算上述3個(gè)建筑物輪廓提取評(píng)價(jià)指標(biāo)，結(jié)果見表1。

表1 建筑物輪廓提取精度評(píng)價(jià)

分析表1可知，匹配度質(zhì)量因子Q的平均值達(dá)到0.848，說明本文研究的建筑物輪廓提取方法具有較高的提取準(zhǔn)確度；面積差A(yù)d和周長(zhǎng)差Pd的平均值分別為0.209和0.134，表明本文方法可以較為完整地提取建筑物輪廓形狀；分析位置精度的計(jì)算結(jié)果可知，建筑物輪廓拐點(diǎn)平面坐標(biāo)的平均距離差Mdd優(yōu)于0.7 m，表明提取的建筑物和建筑物的實(shí)際位置在空間上具有很小的差異性，能夠?yàn)楹笃诮ㄖ锶S重建提供準(zhǔn)確的基礎(chǔ)依據(jù)；此外，匹配度質(zhì)量因子、面積差和周長(zhǎng)差的標(biāo)準(zhǔn)差均優(yōu)于0.2，而Mdd的標(biāo)準(zhǔn)差大于0.5，究其原因是本文研究實(shí)例中多棟建筑物連接在一起，致使提取建筑物拐點(diǎn)具有一定的難度。因此，表1中的標(biāo)準(zhǔn)差計(jì)算結(jié)果表明本文方法在提取不同建筑物輪廓時(shí)具有穩(wěn)定性和普遍適用性。

4 結(jié) 語

建筑物輪廓的準(zhǔn)確提取對(duì)建筑物三維重建至關(guān)重要，本文在深入分析已有算法優(yōu)缺點(diǎn)基礎(chǔ)上，提出了一種綜合改進(jìn)的區(qū)域生長(zhǎng)算法、三維Hough變換算法和α-shape算法的建筑物輪廓提取方法，研究了對(duì)機(jī)載LiDAR點(diǎn)云進(jìn)行去噪、濾波、平面提取、輪廓提取及精度評(píng)定等一系列環(huán)節(jié)，形成了一套完整的從機(jī)載LiDAR掃描數(shù)據(jù)中提取建筑物輪廓信息的數(shù)據(jù)處理模式與方法。隨著建筑物設(shè)計(jì)概念的多元化，現(xiàn)有規(guī)則建筑物輪廓提取方法已無法滿足異形建筑物輪廓提取的需要，對(duì)于不規(guī)則建筑物的輪廓提取筆者將另文研究。