李明磊：雙層優化的激光雷達點雲場景分割方法

《測繪學報》

構建與學術的橋樑 拉近與權威的距離

雙層優化的激光雷達點雲場景分割方法

李明磊1, 劉少創2, 楊歡3, 亓晨3

1. 南京航空航天大學電子信息工程學院, 江蘇 南京 211106;

2. 中國科學院遙感與數字地球研究所, 北京 100101;

3. 武漢大學測繪學院, 湖北 武漢 430079

收稿日期：2017-09-01；修回日期：2017-11-29

基金項目：江蘇省自然科學基金（BK20170781）

摘要：對激光雷達掃描的非結構化點雲進行分割處理，是進行數據組織、重構和信息提取的重要步驟。本文根據點雲表面的局部可微的性質，提出了一種遞進形式的雙層優化分割演算法。首先在黎曼幾何框架下計算點的拓撲關係和距離度量特性，以k均值聚類的方法獲得過分割體素，作為底層分割結果。然後，將點雲的體素模式化為節點，構建最小生成樹，提取節點的高級特徵信息，利用圖優化得到對點雲細節自適應的區域分割效果。通過真實數據進行驗證，並與現有方法比較，證明所提演算法的可行性和先進性。

關鍵詞：點雲分割黎曼幾何超體素最小生成樹特徵提取

Bilevel Optimization for Scene Segmentation of LiDAR Point Cloud

LI Minglei1, LIU Shaochuang2, YANG Huan3, QI Chen3

Abstract: The segmentation of point clouds obtained by light detection and ranging (LiDAR) systems is a critical step for many tasks, such as data organization, reconstruction and information extraction.In this paper, we propose a bilevel progressive optimization algorithm based on the local differentiability.First, we define the topological relation and distance metric of points in the framework of Riemannian geometry, and in the point-based level using k-means method generates over-segmentation results, e.g.super voxels.Then these voxels are formulated as nodes which consist a minimal spanning tree.High level features are extracted from voxel structures, and a graph-based optimization method is designed to yield the final adaptive segmentation results.The implementation experiments on real data demonstrate that our method is efficient and superior to state-of-the-art methods.

Key words:point cloud segmentationRiemannian geometrysuper voxelminimal spanning treefeature extraction

三維點雲分割的目標是依據點的屬性獲得空間同質區域的聚類，分割處理是移動設備導航定位、目標提取識別和場景幾何建模等應用的重要底層技術[1-2]。然而，由於數據的採樣密度非均勻、點雲非結構化分布、雜訊和外點等影響，完全地自動化分割在計算機視覺、計算機圖形學、攝影測量與遙感等研究領域仍然是一項艱巨的任務[3-5]。本文根據原理不同將點雲分割演算法歸為以下4類：

(1) 基於區域增長法(region growing)的分割技術。該方法從種子點出髮根據相似性度量和距離度量對鄰域點判斷是否擴展合併為同一區域[6-7]。區域增長演算法相對易於實現，但它對雜訊敏感，演算法中的增長準則對不同局部區域的細節差異的自適應能力弱，而且種子點選取的不確定性將會影響到分割結果的可靠性。

(2) 基於圖(graph-based)的分割演算法。這類演算法通過將點雲模式化為由節點和反映節點關係的邊組成的圖模型。假設存在一種條件隨機場模型，依據極大後驗估計準則使用圖割演算法將圖模型的部分連接斷開形成獨立的區域[8-10]。該方法對前景和背景進行分割，適用於指定的目標提取，或以監督分類方式實現多目標提取。

(3) 基於模型擬合(model-fitting)的分割演算法。例如以隨機抽樣一致性檢驗(random sample consensus，RANSAC)為依據的參數化幾何元素(平面、橢球和圓柱等)擬合演算法[11-12]，被廣泛應用於建築物立面提取[13-14]。與此類似還有基於霍夫變換(Hough transform)的元素提取技術[15]，該類演算法提取的幾何結構精簡緊湊，可視化效果好。但在參數求解時通常使用統一的閾值從而缺乏尺度適應能力，此外一次分割只能提取出一類基本元素，當場景具有複雜多樣的結構形態時並不適用。

(4) 以機器學習(machine learning)為基礎的分割演算法。文獻[16]使用支持向量機將幾何特徵歸類，文獻[17—18]提出一種以神經網路學習的方法將分割與分類同步結合進行目標識別，文獻[19]將點雲轉化為體素(volume)數據來提取特徵再使用AdaBoost進行訓練和目標提取。這類演算法在訓練過程中需要大量的已有標記數據做訓練樣本，而且分割對象限定為已經訓練過的目標物，因此使該演算法難以靈活地推廣。

前兩類演算法可以歸為自下而上的分割方法，而後兩類則是基於模式引導的自上而下的方法，它們的適用性都存在一定的局限。為了提高點雲分割演算法的抗噪性和對多種類型場景的適用性，解決三維點雲特徵描述困難的問題，本文針對激光雷達(light detection and ranging, LiDAR)數據的特性，在黎曼幾何框架下考慮將體素過分割方法和基於圖的方法相結合，設計分層優化的自動分割演算法，該演算法邏輯清晰易於實現，而且運算速度快，對點雲的局部細節具有良好的自適應性。

1 演算法設計

本文雙層優化策略的點雲分割演算法包括以下3個基本步驟：①通過定義局部坐標系來建立特徵提取的參考基準；②以八叉樹節點為種子點，結合離散點的特徵信息，通過迭代優化聚類實現以超體素為輸出的底層分割；③以超體素為節點建立最小生成樹圖(minimal spanning tree，MST)，以MST為對象進行頂層優化聚類實現場景目標分割。

1.1 局部坐標系建立

LiDAR點雲能反映被觀測目標的表面形態特性，與影像重建的點雲，如多視立體重建(multi view stereo)相比，LiDAR點雲採樣均勻且密度高。儘管屬於非連續空間採樣，但緻密的表面點滿足建立局部歐氏坐標系的條件，可以用黎曼度量來研究局部可微性質，此時點的鄰域關係不再局限於全局坐標系的組織。對點雲中的任意一點，以其切平面為基礎構建局部歐氏坐標框架，統計鄰域內的點集的協方差，由此自適應調整距離度量並提取點特徵向量，使特徵參數可以考慮到數據局部結構的細節差異，避免了全局統一參數約束，有利於提高演算法適用性。

在處理非結構化的LiDAR點雲數據時，描述離散點的特徵的方式具有多種多樣[20-22]，各種方法首先都需要建立起空間索引和鄰域搜索機制。本文中設p是一個查詢點，以Np：{qi|i=1, 2, …,k}表示查詢輸出的一組鄰域點集，其中每個點qi與搜索點p的距離度量滿足‖qi,p‖n≤dmax，為此採用近似最鄰近查詢方法(approximate nearest neighbor，ANN)[23]作為檢索結構實現快速的鄰域搜索。

點雲表面某點處的切平面Π方程可以由點p和對應法向量np表示，其中p表示鄰域點集合Np的重心點，即

。具體而言，本文採用主成分分析(principal component analysis，PCA)演算法求解一個切平面的最小二乘估計，鄰域點集合Np的協因數矩陣C的特徵值和特徵向量的公式如下

(1)

式中，ζi表示點qi的權重，等權觀測情況下通常設為1;C是一個對稱半正定矩陣，也是一個帶權的協方差矩陣，它的特徵值為實數λj,j=0, 1, 2；而它的3個特徵向量vj對應鄰域點集合Np的主分量，構成一個空間正交的框架。通過奇異值分解獲得矩陣C的特徵向量，那麼最小特徵值λ對應的特徵向量v在幾何意義上就與切平面法向量共線，對其進行歸一化即得到單位法向量的估計。

圖 1給出一組點雲表面局部坐標系構建的示意圖，設平面Π為過一點p的切平面，當已知點p的表面法向量估計值時，即，則局部坐標系的Z軸指向與共線，X-Y指向基準定義為

(2)

式中，。

1.2 底層優化超體素(super voxel)分割

本文將點雲分割過程設計分為上下兩層遞進的組織形式，本節介紹底層優化方法，1.3節內容將介紹頂層優化過程。底層結構的分割效果是利用改化的k均值聚類演算法獲得過分割的(over-segmented)點雲聚類，可以表示為點雲體素(voxel)的形式，具體包含以下3步。

(1) 定義相似性度量參數。點雲聚類演算法通常以法向量、曲率和離散度等屬性信息構建特徵向量來進行相似性度量，當具有影像數據時，顏色和紋理信息也被包含其中。本文中僅考慮原始LiDAR點雲形式的數據，此時底層特徵採用法向量和曲率組成的特徵向量表示，具體的一點p的特徵向量表示形式為fp∈R5:[nx,ny,nz,k1,k2]T，其中k1,k2表示表面位置的兩個主曲率值。

對於距離度量而言，本演算法不採用全局三維坐標定義的歐氏距離dp,qi，因為全局歐氏距離不顧及點的局部空間走勢，不利於細節信息的區分。相反，本文利用局部坐標系統，根據點的曲率改進距離度量，記為dproj。如圖 2所示，點p的鄰域點qi到p的距離是利用投影點q′i到查詢點p的歐氏距離進行計算的，其中ρ為曲率半徑。如此，如果某點處的切平面和局部點雲擬合越好，則ρ→∞，qi越接近q′i，改進距離dproj越接近於歐氏距離dp,qi；反之，如果曲率越大，切平面與局部點雲擬合效果差，則鄰域點的距離度量值也將增大。

(2) 聚類中心初始化。對於底層體素分割的另一個關鍵問題就是分割數目的選取和聚類中心的初始化。本文採用自適應八叉樹索引結構來創立聚類中心點集合S：{si|i=1, 2, …,m}，與點p類似每一個種子點si也有一組特徵fsi∈R5。對於三維數據而言，八叉樹可以快速實現鄰域點查詢，而葉子節點本身即是一類以空間距離做度量的廣義上的聚類。演算法集成過程中利用了泊松表面重建演算法[24]中的自適應八叉樹結構做引導，根據經驗值以第6、7級節點(與點密度和場景範圍相關)作為局部聚類中心效果理想，而節點的寬度dr即作為下一步搜索半徑參數。

(3) 聚類搜索閾設計。不同於傳統的k均值聚類演算法將每個點與所有聚類中心做判斷，本演算法中原始點只與到自身距離在3倍搜索半徑(3×dr)內的聚類中心點做相似性判斷。之後運算與k均值聚類相似，利用相似性度量和距離度量的加權距離判別查詢點p所屬的聚類中心sp。所有點判別結束後，根據聚類結果更新聚類中心sp的位置並進行下一輪迭代聚類，直到聚類中心位置收斂。聚類判別公式為

(3)

式中，

和組成聯合度量函數；λ為兩類度量的歸一化參數。取得滿足最小度量值的中心si即為點p所屬的聚類中心sp。通過改進距離度量和聚類中心判別策略，演算法大大提高了聚類結果的可靠性，而且基於八叉樹選擇種子點引導分割的體素片段在局部表現出很好的尺度一致性。

1.3 頂層優化圖聚類

底層優化得到的是過分割的點雲體素集合，這類集合表達了點雲小尺度的相似性關係，為目標級別的分割提供支撐。頂層優化的數據操作對象不再是獨立的點而是體素，具體過程分為以下3步。

(1) 挖掘體素的高級特徵，而體素的特徵提取目前仍是一項難題。此外，體素的鄰接關係的定義不再單純以距離為依據，因而鄰接關係的確定也相對複雜。在底層分割時，每一個體素對應產生了一個聚類中心點，該中心點被用做體素鄰域判斷的初始條件。所有聚類中心集合形成了一個新的小型點雲集合，因此同樣利用ANN演算法，可以快速的獲得聚類中心的鄰域中心集合，而這些中心鎖定了一組初始的體素的對應關係。此時的ANN演算法僅僅利用聚類中心的距離做搜索依據，而體素的固有屬性遠不止距離信息，因此在初始鄰域搜索時利用較大的搜索半徑獲得一個數量冗餘的鄰域集合{V′neighbor}。下一步通過高維特徵的相似性縮減鄰域集合的數目，建立緊湊可靠的體素鄰域關係。

一個體素的首要特徵是利用其中所有點依據最小二乘法擬合得到的平面法向量Nvoxel。第2類特徵是體素所有點集的協方差矩陣的3個特徵值[λ，λ1，λ2]T，它反映了體素的空間分布特性。兩類特徵組合構成了一個6維的特徵向量，據此每個體素在初始鄰域{V′neighbor}中進行相似性判別，僅保留屬性最接近的k個體素做最終的鄰域集合{Vneighbor}。

(2) 當鄰域關係確立後，進行拓撲圖構建。本文利用Kruskal坐標系[25]建立以體素為節點的MST，實現對體素結構的拓撲關係索引。MST的建立步驟是首先設每一個體素為一個節點e，將其所有鄰域節點{Vneighbor}按邊的距離權值從小到大排序，循環從權值最小的邊開始遍歷每條邊直至圖中所有的節點都在同一個連通分量中，該連通分量所經過的所有路徑和節點即構成了一個MST。

(3) 最後從MST中的起點出發，根據邊的權值按區域增長法進行聚類判別，當遇到增長終止條件時，自動開始新的區域聚類，從而實現頂層非監督分割聚類。

2 試驗與分析2.1 數據處理與演算法實現

試驗環節採用了室內和室外兩組點雲數據分別進行處理分析，運算平台為win10 64位系統，處理器為i7-6500U搭配8 GB內存。演算法功能採用C++編程實現，利用QT搭建軟體界面，OpenGL作為三維繪圖引擎，並開發實現八叉樹索引和RANSAC幾何元素檢測等相關功能。

圖 3通過室內數據試驗展示了處理過程的中間結果。該數據是一組利用Leica ScanStation C10掃描儀獲取的卧室點雲，共包含45.1萬個空間點，掃描距離在1.5 m到3.5 m，採樣密度約為9000點/m2。圖 3(a)給出原始點雲的一個快照；圖 3(b)是初始化的過分割聚類中心點，可以發現聚類中心點概括了原始點雲的布局；圖 3(c)是經過迭代優化後的底層過分割聚類結果，其中不同的體素以隨機生成的顏色進行區分標識；圖 3(d)展示了以體素為基本節點的MST圖構建情況，該圖是區域增長法的路徑尋優基礎。通過建立以超分割體素為節點的MST圖，一方面保持了原有數據的幾何信息和拓撲信息，另一方面從點雲到體素節點極大地減少了數據量，為提高數據處理效率提供了重要的支撐。

2.2 試驗比較分析

為驗證本文演算法的先進性，試驗中將基於RANSAC的參數化幾何元素擬合分割演算法[11]與本文雙層優化演算法進行了比較，其中RANSAC點擬合距離閾值設為0.1 m。兩種演算法的分割結果如圖 4所示，其中圖 4(a)、(b)分別表示RANSAC分割結果和本文演算法的結果，圖 4(c)、(e)是一個沙發的局部放大效果，圖 4(d)展示了本文演算法的超體素分割的局部效果。由於沙發表面並非規則的幾何形狀，從圖 4(c)、(e)中可以比較發現，RANSAC擬合參數化表面的分割結果不理想，儘管通過設置距離和法向量的擬合殘差閾值，可以實現對場景的局部分段分割，但單一模式的平面表達會使分割得到的對象失去原有的細節特性。相反，本文演算法利用點雲過分割和體素聚類分割兩層優化步驟，將掃描點雲的局部可微性質保留在體素內部，並在適當的過渡區域斷開MST圖，實現點雲的非監督分類，具有較好的細節保真度。

兩種演算法的定量比較統計結果見表 1，包含演算法分割耗時、分割數目和擬合殘差指標。從表 1中可以得出本文演算法比較RANSAC擬合分割演算法有優異表現。在效率方面，RANSAC演算法需要多次隨機採樣確定擬合參數，而且每提取一個分割面都需要單獨全局擬合計算，計算耗時相對較多；而本文利用超體素作為圖節點減少數據量，極大地降低了運算量。另外，RANSAC演算法分割的目標局限為參數化的基本幾何形狀，因此分割數目有限，無法對複雜的場景形狀進行有效表達。兩種方法的擬合殘差接近，這是由於演算法本身進行了擬合距離閾值設置，限定了最大的擬合殘差，因此具有統一的擬合殘差量級。

表 1兩種分割演算法比較數據統計Tab. 1Statistics of two different algorithms

表選項

此外，為驗證本文演算法的廣泛適用性，利用車載LiDAR點雲數據進行了室外場景分割測試，原始點雲和分割結果如圖 5所示。演算法通過自適應八叉樹可以在突變區域提取可靠種子點，根據種子點的鄰域點集的局部可微性進行聚類，從而生成有效的超體素，避免丟失場景中的突變地形和地物信息。如圖 5(b)所示，路燈和電線杆被正確地分割提取，充分證明了本文演算法對多種場景的適用性。本文演算法的局限性在於缺乏語義層次的判別，僅從幾何層面分析可能將一類目標物分割為幾個局部連貫的片段，如圖 5(b)圖框選區域中將一輛卡車的車頭和車身分割成兩個區域。

3 結論

本文設計了一種面向LiDAR點雲的自動分割演算法，首先對底層點雲進行迭代均值聚類獲得過分割，形成超體素集合。然後，在頂層以體素為節點建立拓撲圖模型，利用MST結構極大地簡化了圖的結構，提高區域增長演算法的效率。本文演算法不需要限定分割結構是平面或柱面等幾何類型，因此適用於多類型目標分割，能夠同時兼顧點雲的局部細節和全局整體分布情況，適合多類應用場景的推廣。下一步的研究工作將提取分割結果的高級特徵進行語義層次的判別，實現點雲目標識別。

【引文格式】李明磊，劉少創，楊歡，等。雙層優化的激光雷達點雲場景分割方法[J]. 測繪學報，2018，47(2)：269-274. DOI: 10.11947/j.AGCS.2018.20170493

喜歡這篇文章嗎？立刻分享出去讓更多人知道吧！