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

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

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

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

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

1.1 改進的區(qū)域生長算法濾波地面點

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

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

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

(2) 計算種子點的平均高程P0，并通過3個種子點生成種子平面。

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

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

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

(6) 為了繼續(xù)生長，將初始種子點中任意點歸為地面點(此處設置為P1)，儲存到地面點集合中，從候選種子點集合中選出最小值ρi所對應的點作為新的種子點P4。

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

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

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

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

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

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

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

(1)

(2)

式中，(Ox，Oy)為滾動點的圓心坐標；(x1，y1)與(x2，y2)分別為待判斷點P1與P2的平面坐標。

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

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

(2) 計算點集S1內其他點到點P0的距離，當所有點的距離都大于α時，將P1、P2定義為輪廓點，并轉至步驟(4)；當出現(xiàn)點的距離小于α的情況時，則轉至步驟(3)。

(3) 對點集S1內下一點重復步驟(1)—(2)，直至S1內所有點均有所判斷方可結束。

(4) 對點集S內未被判斷的點重復步驟(1)—(3)，直至S內所有點均有所判斷方可結束。

2 實例計算

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

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

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

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

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

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

3 精度評定

為了對建筑物輪廓提取精度進行定量評價，參考文獻[15]并結合實際情況,本文提出匹配度、形狀相似度及位置精度共3個精度評價指標來評定建筑物輪廓的提取精度。

3.1 匹配度

本文采用圖像分類精度評估參數(shù)中的質量因子Q作為匹配度的評價指標，即

(3)

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

3.2 形狀相似度

通過定義面積差Ad和周長差Pd描述提取的建筑物形狀相似度,即

(4)

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

3.3 位置精度

選用建筑物輪廓拐點平面坐標的平均距離差Mdd對位置精度進行評價，即

(5)

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

對整個測區(qū)實例(如圖8(a)所示)所有建筑物采用上述方法進行建筑物輪廓提取，提取結果如圖8(b)所示，同時計算上述3個建筑物輪廓提取評價指標，結果見表1。

表1 建筑物輪廓提取精度評價

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

4 結 語

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