一種分步式建筑物屋頂面點云高精度提取方法

周欽坤 岳建平 李樂樂 楊 恒

1 河海大學地球科學與工程學院，南京市佛城西路8號，211100 2 蘇州工業園區測繪地理信息有限公司，江蘇省蘇州市蘇虹中路101號，215000

主流的基于機載激光雷達數據的建筑物三維模型重建方法包括基于模型驅動的三維重建與基于數據驅動的三維重建，其中建筑物屋頂面的提取是數據驅動方式中最為關鍵的一個步驟，提取質量的好壞直接關系到后續建筑物拓撲重建的成敗[1]。由于受機載激光雷達點云分布不連續、不規則及建筑物多層次、多面片的復雜結構等因素的影響，從機載激光雷達點云數據中提取高精度建筑物屋頂面仍面臨著較大的挑戰。

學者們圍繞屋頂面的高精度提取開展了眾多研究，提出了數據聚類[2-3]、區域增長[4-5]、模型擬合[6-9]、能量函數最小化[10]等方法。本文針對現有方法存在建筑物屋頂面提取精度較低、適應性較差等問題，提出了一種分步式建筑物屋頂面點云高精度提取方法。

1 分步式提取方法

首先通過預處理選取可靠平面點，利用K-means算法實現可靠點在法向量空間上的聚類；再通過逐步平面估計，提取初始屋頂面片；在此基礎上，進行面片的合并與未標記點的歸屬判斷。

1.1 預處理

點云法向量是建筑物屋頂面提取的重要依據，其在屋脊線等屋頂面交互處通常難以被準確計算[11-12]，在法向量空間域中表現為錯亂分布的點，不利于后續的聚類操作。為保證面片提取的準確性，首先對屋頂點云進行預處理，以得到法向量下可靠的平面點，主要包括以下步驟：

(1)

(2)

2)可靠平面點的選取。每個點均被賦予一個可靠性指標，理論上該值越小，該點越可靠，可選取可靠性指標較好的點組成可靠平面點集seed：

seed={p|ασ<αn,n=L/k,

k=max(Kg,Kw)}

(3)

式中，αn為可靠性指標降序排列時第n個點對應的ασ值，L為建筑物屋頂面點云總個數，Kw為屋頂面片聚類數，數值上等于后續K-means算法中的參數K，參數K隨后續的迭代計算而動態確定。Kg為固定值，可根據屋頂面實際情況在5～7中取值，以防止因去除過多的點而使屋頂面提取困難。

圖1(a)～1(c)分別為建筑物屋頂面影像、屋頂點云可靠性分布及可靠平面點選取結果，圖1(c)中黑色點云為選取的可靠平面點，經過預處理將法向量分布離散的點云去除，對應圖1(d)中紅色點云。典型建筑物屋頂由多個平面構成，每個平面點云的法向量會指向特定方向。圖1(d)為屋頂點云的法向量空間分布，即屋頂點云法向量在高斯球上的分布情況，紅色為根據可靠性指標去除的點云，其他不同顏色代表不同屋頂面的可靠平面點，可見同一屋頂面上可靠平面點在法向量空間上是較為集中的，不同屋頂面上則相互分離，因此可通過點云法向量空間分布提取建筑物屋頂面片。

圖1 可靠性指標分布及可靠平面點選取Fig.1 Distribution of reliability index and selection of reliable plane points

1.2 屋頂面片初步提取

K-means算法是無監督的聚類算法，對于屋頂點云法向量這樣的類球形分布，其聚類效果較好，且實現簡單，調試的參數僅為聚類個數K。本文通過設置聚類指標實現參數K的動態確定，基本原理如下：K-means算法要求簇內各點盡量緊密連接在一起，而簇間距離盡可能大，當K值恰好等于實際所需的聚類數時，屋頂點云在法向量空間上正確聚類；隨著K值的繼續增大，某一屋頂面點云將會被分成數簇，即會存在法向量指向相近的數個簇，此后最大的簇間單位法向量點積會趨近于1且保持穩定。所以可取最大簇間單位法向量點積作為聚類指標，K值即為使聚類指標趨近于1且保持穩定的前一個聚類個數，如圖2中紅色標記點橫坐標為對應的聚類個數。

圖2 K-means算法聚類指標與聚類個數折線圖Fig.2 K-means clustering index and clustering number

通過聚類操作，將法向量指向近似的點云歸為同類，但若建筑物存在平行或近似平行的屋頂，即法向量指向相近但空間上相離的屋頂，則需對其進行進一步區分。本文采用逐步平面估計法進行此類屋頂面的提取，具體步驟如下：

1)取某一聚類全部點云單位法向量加權得到的聚類法向量n，并選此聚類中ασ最小的點pmin(ασ)(可認為此點是平面性最可靠的點),pmin(ασ)∈R3,構建平面方程P1：

n·(x,y,z)T-n·(pmin(ασ))T=0

(4)

2)計算聚類中各點到平面P1的距離，若小于給定閾值，則判斷此點為平面上一點，進行聯通性分析，取最大聯通子集為當前平面的內點；

3)對內點進行最小二乘擬合，確定平面P2；

4)計算聚類中各點到平面P2的距離，重復步驟2)的操作，記第i次估計平面P2時P2內點個數為num2i；

5)若相鄰2次平面P2檢測的內點數量趨近恒定，即|num2i-num2i-1|

6)從聚類中移除內點，并重復步驟1)～5)進行迭代，由于屋頂平行面片個數一般有限，迭代次數一般不大于10。

本文提出的逐步平面估計是一個尋找最佳屋頂平面的過程，由于法向量的計算誤差與聚類中近似平行面片的存在，從整個聚類中得到的某一屋頂平面較為粗糙。通過步驟1)法向量加權與最可靠點的選取，使得P1大體逼近屋頂面片，P1的內點即為某一屋頂面上的點云，由此部分點可擬合出更優的屋頂平面P2，得到更多的屋頂面內點，通過數次內點擬合即可逐步修正平面P2，使其成為最佳平面。上述操作進行了2次連通性分析，第1次是為了剔除由于平面P1的誤差納入的其他屋頂點，進而依據剩余內點擬合得到更優平面；第2次是為了去除在數學意義上共面而在空間上分離的點，得到純凈的單一屋頂面片。逐步平面估計屋頂面的提取結果見圖3(c)，可以看出，該算法能較好地分離平行屋頂面。與RANSAC算法相比，該方法不用設置最小面片點數且不用隨機選點以構成平面模型，避免了由于選點不在同一平面而造成的大量不必要的計算。

圖3 聚類結果及逐步平面估計結果Fig.3 The results of clustering and stepwise plane estimation

1.3 面片后處理

面片后處理主要包括面片合并及未標記點的歸類。由于法向量計算不準確及K-means聚類誤差等因素，部分屋頂面經過粗提取后可能形成多個面片，如圖4(a)所示。由于屋頂面A、B近似平行且存在上述誤差，使得屋頂面B中的一部分與屋頂面A聚為一類，在逐步平面估計后，屋頂面B將會形成2個破碎面片，如圖4(b)所示。為得到完整的屋頂面，需合并屬于同一屋頂面的多個面片。令已找到的屋頂面片為M={m1，m2，…，mn}，尋找面片mi的相鄰面片，首先定義面片鄰域Ner，即由面片中每個點的鄰域點面片號構成的無重復標號集合，搜索面片mi的Ner，其中每個標號都代表面片mi的一個相鄰面片，計算面片mi與各個相鄰面片的最小二乘擬合平面殘差平方和，若殘差平方和小于閾值，則將面片mi與相應相鄰面片合并。重復上述過程2～3次，即可實現屋頂面片的完整提取，面片合并結果見圖4(c)。

圖4 屋頂面片合并Fig.4 Merging of roof patch

圖5 未標記點歸類Fig.5 Classification of unmarked points

2 實驗與分析

2.1 實驗數據與方法

為驗證本文算法的可行性與穩定性，選取Vaihingen機載LiDAR數據集中部分建筑物的點云數據進行實驗。由于本文研究主體為建筑物屋頂，因此地面點、植被點及建筑物立面點等數據均已被剔除。

2.2 實驗結果及定性分析

各算法的提取結果如圖6所示，圖中不同顏色代表提取的不同屋頂面，其中紅色點為未能提取出的屋頂面或屋頂附屬物點。通過與航空影像進行對比可以看出，RANSAC算法及區域增長算法的提取效果與本文算法相比差異較大。

圖6 單個建筑物屋頂面提取結果Fig.6 Extraction results of roof surface for single building

區域增長算法僅依靠種子點與待判斷點之間的法向量夾角及粗糙度等局部特征進行生長，所以分割結果與局部特征的計算精度息息相關。由于局部特征計算的不穩定性，導致建筑物屋頂將存在不同程度的邊界模糊問題，如建筑物4(矩形框區域，下同)。同時，對于較小的屋頂面如建筑物3，計算出的法向量指向雜亂，無法實現有效的聚類，而在過渡平緩的相鄰屋頂面處，由于局部特征變化較小，會形成屋頂面的欠分割，如建筑物1和2。此外，在鄰接關系復雜的區域容易形成過分割，如建筑物6。

RANSAC算法的每次隨機采樣估計模型會獲取點數最多的屋頂面，兩屋頂面相交處的點被優先劃分到較大的屋頂面，形成屋頂競爭現象，造成邊界不明確，如建筑物5。由于RANSAC算法需要設置最小面片點數閾值，而全局最優的閾值難以找到，當采用較大的閾值時，較小的屋頂面由于達不到閾值條件而無法提取，如建筑物3和10；當閾值設置較小時，較大的屋頂面可能會因為過早地滿足閾值條件而變得破碎，如建筑物9。此外，受限于算法本身的原理，對于曲率變化較大的屋頂區域，如建筑物7，無法實現有效提取。

相較于以上2種算法，本文算法的提取結果更加優異，拓撲關系更加清晰，且不會受屋頂面大小差異、數量多少及復雜程度等的影響，具有更高的穩定性。但本文算法本質上是一種平面提取算法，對于如建筑物7的屋頂面，與RANSAC算法一樣，提取效果較差。此外，本文算法依賴于可靠平面點，若屋頂面過小且點密度太低，導致計算出的法向量不準確，則屋頂面將不存在可靠平面點，該屋頂面的提取也會失敗，如建筑物8。

2.3 定量分析與評價

建筑物屋頂面在空間上為平面，在機載LiDAR點云數據中表現為點集，故本文從點與面2個維度對提取結果進行定量分析，在面維度上借鑒文獻[13]中的質量評價標準，在點維度上參考文獻[3]與[14]的評價方法。

由表1(單位%)可知，RANSAC算法與區域增長算法在提取完整率上與本文算法精度相仿，其原因為對兩者增加了連通性分析與未歸類點的二次判斷等措施；但在正確率與提取質量上，RANSAC算法及區域增長算法與本文算法依然存在較大差距，尤其是區域增長算法，說明這2種算法存在較為嚴重的過分割與欠分割等問題。結合表2可以發現，本文算法能正確提取絕大多數屋頂面，而RANSAC算法與區域增長算法的提取效果較差。觀察建筑物影像可知，本文算法未能正確提取的主要為面積較小(包含的腳點數較少)的屋頂面或曲率變化較大的非平面屋頂面。此外，RANSAC算法與區域增長算法對過渡平緩及拓撲關系復雜的屋頂面也存在錯誤提取的情況。表1和2中維度的精度統計結果共同驗證了§2.2中的定性分析結論，并進一步證明本文算法對不同復雜程度的建筑物屋頂面的提取具有適應性與準確性。

表1 屋頂面點維度提取精度統計

表2 屋頂面面維度提取精度統計

為比較不同算法的提取效率，統計各算法對不同屋頂面的提取時間，結果如圖7所示。由圖可知，本文算法的提取效率與區域增長算法相近，明顯優于RANSAC算法。雖然本質上同為平面提取算法，但本文算法首先將近似平行屋頂聚類，然后通過逐步平面估計等措施提取屋頂面，避免了RANSAC算法因對建筑物屋頂面整體選點構建平面模型而造成的大量不必要計算，故提取時間較RANSAC算法明顯縮短。

圖7 各算法提取時間統計Fig.7 Statistics for extraction time of each algorithm

通過定性分析與定量評價可知，本文算法明顯優于其他2種算法，對于屋頂結構復雜的各類建筑物也有較好的提取效果，其主要原因為：1)將可靠平面點與非可靠平面點分別進行處理，在將可靠平面點聚類得到平行屋頂面時，可避免非可靠點對聚類的影響；2)采用由粗到精的分步式處理方式，即先將平行屋頂聚類，然后通過逐步平面估計等措施提取初始屋頂面，再進行細化處理得到建筑物完整屋頂面；3)面片合并較好地解決了屋頂面的過分割，保證了屋頂面的完整性，對未標記點進行歸屬判斷，避免了屋頂面的競爭現象，可得到更加清晰的面片邊界。

3 結 語

本文針對復雜建筑物屋頂面提取存在的問題，提出一種分步式建筑物屋頂面點云高精度提取方法。實驗結果表明，相較于傳統方法，本文方法的提取結果更優，且提取效率較高，對不同復雜程度的建筑物屋頂面均能取得較好的提取效果，具有較強的適應性與準確性。由于本文算法本質上也是平面提取算法，故無法有效提取曲面屋頂，另外本文算法依賴于可靠平面點的選取，不能有效提取不存在可靠平面點的小屋頂面。如何準確估計點云法向量以得到小屋頂面的可靠平面點，將是后續研究的重點。