基于機載激光雷達點云的交互式樹木分割與建模方法研究

張玥焜，余文杰，趙習之，呂艷星，馮偉桓，李崢嶸，胡少軍

(1.西北農林科技大學信息工程學院，陜西 楊凌 712100；2.北京星閃世圖科技有限公司，北京 100085)

樹木建模在計算機圖形學、虛擬現實、農林業(yè)等領域具有廣泛地研究與應用價值。一方面，通過樹木建模獲取其生物量、葉面積指數和產量等表型參數，有助于掌握森林資源信息，并進行高效管理與決策；另一方面，計算機生成的真實感樹木模型可增加虛擬現實場景的沉浸感，可直觀向非專業(yè)人員展示農林業(yè)復雜模型，并為機器人噴藥等提供重要的空間信息。

樹木建模經過近50 年的研究，目前已形成基于點云[1-6]、圖像[7-14]、規(guī)則或過程[15-19]及基于草圖的方法[20-23]，其中最新的文獻為楊垠暉和王銳[24]的“樹木的真實感建模與繪制綜述”。最近，國內學者在該領域的研究非常活躍，WANG 等[25]采用變分優(yōu)化方法生成了真實感強的樹模型；WANG 等[26]提出的基于插值混合模型可從2 棵給定模型中構建新的樹模型；WANG 等[27]利用樹形空間的統(tǒng)計建模方法，降低了樹木建模復雜度，隨后WANG 等[28]又提出了一種更加完善和精確的樹木建模方法，即通過樹形空間對樹木模型進行統(tǒng)計分析，實現了從少量樹木生成多棵樹木的目的；張曉鵬團隊研究了基于樹的可見點云合成非可見點云數據并生成三維模型的新方法[6]，進一步提出在多視圖像生成稠密樹點云后采用基于過程的方法構建真實感植物模型的新方法[29]。

早期樹木的建模主要采用基于規(guī)則的方法，如文獻[15]開發(fā)的L 系統(tǒng)，隨著激光雷達掃描技術的蓬勃發(fā)展，最近基于點云的建模方法成為了新的研究熱點。文獻[17]采用空間殖民算法從目標樹冠體中隨機分布的點中構建了真實感強的樹結構模型。文獻[19]和文獻[23]拓展了文獻[17]的方法，提出了更加高效的方法用以生成三維樹模型。由于空間殖民算法提取的骨架點不是從原始點云中進行選擇，難以保留樹的分枝細節(jié)。文獻[12]根據圖像序列重建樹點云，并通過復制可見分枝來生成被遮擋部分的分枝。文獻[1]基于地面激光掃描點云生成樹的分枝結構，首先利用最短路徑方法構建子圖，通過連接相鄰簇質心的方式提取樹的骨架結構，最后對點云缺失部分的分枝進行骨架合成。文獻[2]提出了針對多棵樹點云的自動全局優(yōu)化算法，與文獻[1]的算法相比，重建時間由分鐘縮短至秒。隨后，文獻[3]又提出了一種新的基于葉簇的樹模型生成方法。

上述基于點云的建模方法要求點云中樹枝骨架清晰，因此大部分研究采用地面激光雷達獲取的高密度樹點云，其方法采集到的樹點云密度低，難以有效實現重建。文獻[4]提出了一種基于雙約束的貪婪算法實現了稀疏激光雷達樹點云的自動幾何重建，但該方法建模魯棒性較差，難以解決殘缺樹點云的建模問題。

在樹木點云分割方面，可采用基于區(qū)域增長、基于分水嶺、基于最小生成樹和基于規(guī)范的分割方法實現[30]。文獻[2]將三維點云投影到地面上，并將高密度點識別為根節(jié)點，然后利用Dijkstra 最短路徑算法求解生成樹，最后用零權值去除根邊實現了單棵樹的分割，但基于最小生成樹的方法不適于樹間距小且點云稀疏情況下的單木分割。文獻[4]針對稀疏機載激光雷達(airborne lidar scanning，ALS)樹點云，提出了一種基于規(guī)范分割的方法自動實現了大規(guī)模稀疏樹點云的分割。但上述2 種自動分割方法難以避免過分割現象，無法滿足小規(guī)模樹點云局部分割和深度學習標注過程中對單棵樹點云分割精度的要求。

本文主要針對ALS 采集的低密度樹點云，面向單棵樹點云的較高精度分割與魯棒性建模應用領域，提出一種交互式分割和建模新方法，避免自動化分割和建模可能造成的過分割及缺失點云建模失敗等問題。首先，基于高度圖的交互式分割方法，分層從一片稀疏點云中提取出單棵樹點云，然后基于改進的空間殖民算法，通過交互式調節(jié)相關建模參數實現單棵樹的建模，使得原始點云細節(jié)特征得以保留，最終獲得忠實于原始點云的樹模型。

1 點云數據獲取

ALS 樹點云采用無人機Swiss Drones’ Dragon (載重35 kg，最大飛行高度2 000 m，滿載續(xù)航時間1 h)搭載激光雷達掃描設備Riegl Vux-Sys (激光脈沖頻率最高550 kHz，掃描速度10～200 line/s，位置精度0.05～0.30 m)獲取，采集的樹點云密度為40～50 點/m2。為便于后期交互式分割處理，首先采用lastools (https://lastools.github.io)進行地面點云和樹木點云的提取。如圖1(a)所示，ALS 系統(tǒng)具有一次性掃描大規(guī)模點云的優(yōu)勢，但采集的樹點云內部出現枝干稀疏或缺失等現象(圖1(b))，從而影響后期分割與重建。

圖1 機載激光雷達掃描系統(tǒng)獲取的樹點云((a)掃描系統(tǒng)獲取的山脈樹點云片斷；(b)分割出的單棵樹點云) Fig.1 Tree point clouds acquired by ALS system ((a) A slice of mountain point clouds with thousands of trees scanned by ALS system; (b) Segmented single tree point cloud)

2 交互式樹點云層次分割

由于ALS 樹點云具有密度低、枝葉交叉錯綜復雜等問題，采用自動分割方法會造成過分割，難以滿足未來基于深度學習樹點云分割研究中對人工標注方面的精度需求。因此，本文提出一種交互式分割方法，可利用高度圖和點云切片技術，通過交互式調整每個切片生成精度較高的單棵樹點云。主要包括生成高度圖、交互式粗選單棵樹區(qū)域和交互式精細調整切片凸包3 個步驟。

2.1 生成高度圖

ALS 系統(tǒng)采集的樹點云中樹冠部分較多，而一般樹冠形狀呈現樹干中心位置高、周圍低且從中心到四周高度呈遞減趨勢，基于這一特點，可將點云的最高位置投影到地平面上，顏色用對應高度值表示形成高度圖，從而便于分辨出單棵樹點云。

圖2 選取一片ALS 樹林點云，點云中任意一點的位置和地面高度分別記為Pi=(Pix,Piy,Piz)和Gi(i=1,2,…,n)，Pmin和Pmax分別為點云包圍盒最小和最大位置點，Gmin和Gmax分別為所有樹木點云相對地面的高度最小值和最大值，Δs為分辨率，將點云離散化投影的二維高度圖長寬分別為W=(Pmax…x-Pmin…x)/Δs和H=(Pmax…z-Pmin…z)/Δs，則對應高度圖生成算法為：

圖2 基于機載雷達點云樹點云生成的點云高度圖 Fig.2 Height map generated from ALS tree point clouds

根據算法1 生成的高度圖(圖2)，可清晰發(fā)現每棵樹的大致中心位置和覆蓋區(qū)域，為下一步進行單棵樹的粗選提供了依據。

2.2 交互式粗選單棵樹區(qū)域

在生成的高度圖中可直觀發(fā)現單棵樹的大致區(qū)域，基于本文開發(fā)的交互式分割系統(tǒng)，鼠標單擊高度圖區(qū)域選擇樹冠中心(a,b)，然后在樹的邊緣位置再次單擊鼠標可確定樹冠的大致區(qū)域半徑r，對應到三維空間系統(tǒng)將聯動生成一個圓柱形區(qū)域，其中圓柱上下底中心位置為(aΔs+Pmax…x,Pmax…y,bΔs+Pmax…z)和(aΔs+Pmin…x,Pmin…y,bΔs+Pmin…z)，圓柱半徑為rΔs。計算落入圓柱體內的點并記為分別為圓柱體內點云包圍盒最小和最大坐標值點。設定Δh為高度方向切片分辨率，可將落入圓柱體內的點分為個切片，落入第i個切片的點Qi(i=1,2,…,N)對應高度范圍為針對切片i內的點Qi，采用單調鏈凸包算法求得對應凸包，步驟如下：

(1) 將給定Qi按x坐標升序排序，若x坐標值相同則按y坐標升序排序；

(2) 創(chuàng)建凸包上部。將排序后的點按照x坐標從小到大加入凸包S，若新加入點使得S不再是凸多邊形，則逆序刪除之前加入S的點，直到S重新成為凸多邊形為止；

(3) 創(chuàng)建凸包下部。將排序后的點按照x坐標從大到小加入凸包V，若新加入點使得V不再是凸多邊形，則逆序刪除之前加入V的點，直到V重新成為凸多邊形為止；

(4) 合并集合S=S∪V形成最終Qi對應凸包。

圖3 右上角藍色線段形成區(qū)域即為切片i中的點形成的凸包，觀察可見，粗選點云分割精度低，出現了不相關2 簇點云被分到同一切片的問題。

圖3 根據高度圖交互式粗選單棵樹點云柱形區(qū)域生成對應點云切片與凸包 Fig.3 Corresponding 2D convex envelope and slice of a point cloud generated from a rough single tree point cloud in a cylinder space chosen from a height map

2.3 交互式精細調整切片凸包

為解決圖3 中點云切片誤分割問題，針對凸包，系統(tǒng)增加了鼠標交互式調整或刪除凸包多邊形頂點的操作，并動態(tài)判斷哪些點落入了新多邊形區(qū)域。設定多邊形頂點為Si(i=1,2,…,N)，判斷點p是否在多邊形內部，見算法2。

在調整凸包頂點的過程中，需動態(tài)判斷Qi中所有點是否在新的多邊形內，如圖4 所示，若在多邊形內部，則全部選定為當前單棵樹切片中的點云，至上而下依次調整所有切片如圖5(a)所示，并依次連接相鄰切片形成如圖5(b)所示的三維包絡，導出其點云即形成分割后的單棵樹點云。

圖4 交互式調整切片凸包確定單棵樹點云區(qū)域 Fig.4 Refine single tree point cloud in a selected slice

圖5 連接各層調整后的切片點云形成單棵樹點云包絡 ((a)調整好的所有點云切片多邊形區(qū)域覆蓋的點即為目標點；(b)對應切片生成的三維包絡) Fig.5 Generate 3D hull of a single tree point cloud by connecting adjusted slices ((a) Slices of a tree point cloud; (b) 3D hull generated from the slices)

3 交互式樹木建模

由于ALS 樹點云密度低、枝干結構不清晰且常在樹干部分存在缺失信息，傳統(tǒng)的用于高密度樹點云的重建算法不能直接用于稀疏ALS 樹點云的魯棒重建，空間殖民算法具有算法簡單、魯棒性好的優(yōu)點[17]，后面經過改進已成功在基于草圖的建模中得到有效應用[19,23]，因此本文采用該算法用于生成樹枝骨架結構。

3.1 改進的空間殖民算法

空間殖民算法基于樹木生長空間競爭的原理，能較好地提取被冠層遮擋的枝干部分骨架信息。假設p為初始骨架點，對應可生長空白空間點集合為S(p)。空間殖民算法首先需計算骨架點到其受影響的每一空白空間點q之間的方向向量，然后對方向向量求和并標準化為向量n，向量n用于確定新樹枝的生長方向，設定骨架點間距為d，則新的骨架點位置可表示為p+dn，然后設定刪除閾值dk，若空間待搜索點與新骨架點的距離小于刪除閾值，則將待搜索點刪除。此外，默認的空白空間點搜索區(qū)域是一個球域，容易生成角度過大的樹枝分杈且偏離點云位置，因此需改進該算法[5,19]，將搜索空間約束為半頂角角度θp、影響半徑為rp的錐形區(qū)域，即

其中，vp?p和vpq分別為上一級樹枝骨架和當前骨架的生長方向。改進的空間殖民算法3 流程如下：

3.2 樹枝粗度計算與幾何表示

空間殖民算法類只能提取生長空間點或樹點云對應的三維骨架，為獲取骨架對應的樹木幾何模型，本文基于異速生長模型計算樹枝粗度[1]，然后采用廣義圓柱體進行樹枝幾何重建，從而實現樹枝的幾何表示，樹葉根據葉序規(guī)則添加[4]，最終得到樹木的三維模型。

假設樹干最大粗度為droot，當前樹枝長度為Lp，其上第i個子樹枝長度為Lp,i，則基于異速生長模型計算的當前樹枝粗度dp為

其中，Subtree(p)為樹枝p上的所有子樹枝集合，λ為值在1～2 之間的粗度調節(jié)因子，由樹種決定。設置樹干最大粗度后，采用深度優(yōu)先搜索算法可計算所有樹枝粗度。

3.3 交互式建模

自動化建模魯棒性較差，默認的參數組合難以解決殘缺樹點云和自然樹木的真實感建模問題。因此本文設計了如圖6 所示的用戶界面，通過交互式調節(jié)約束角度θp、刪除閾值dk和影響半徑rp等3 個參數，用戶可快速生成不同參數組合下形狀各異的重建樹幾何模型。此外用戶還可在圖6 對話框中交互式設置樹點云重采樣密度、樹干最大粗度droot、樹葉類型和葉序規(guī)則[7]，通過少量參數設置即可實時生成不同類型的樹木模型，為建模人員提供了靈活選擇。

圖6 交互式樹木建模用戶界面(框區(qū)域內為可交互式調節(jié)的約束角度、刪除閾值和影響半徑參數) Fig.6 User interface of interactive tree modeling (The interactively adjustable constraint angle,kill distance and influence radius are shown in the red rectangle)

4 實驗結果分析與討論

本試驗平臺為Windows 7 操作系統(tǒng)，計算機配置為CPU i5-7200U 2.50 GHz，內存8.0 G，顯卡 NVIDIA GeForce 940MX，所有實驗樹點云對象均采用第2 節(jié)介紹的ALS 系統(tǒng)采集。交互式樹點云分割和建模軟件演示視頻參見鏈接:https://cie.nwsuaf.edu.cn/docs/2020-07/08679ab69a1d49a3aebf95a6a7f 79214.zip.

4.1 交互式樹點云分割結果

基于本文算法，測試了如圖7 所示難以自動分割、從單一角度觀察難以識別的小片ALS 樹點云，實驗結果表明，交互式分割能夠有效實現單棵樹點云的分割。

圖7 交互式分割測試結果((a)調整好的點云切片；(b)形成的三維包絡；(c)提取的分割后單棵樹點云) Fig.7 Result of interactive segmentation ((a) Adjusted slices of a tree point cloud;(b) 3D hull generated from the slices; (c) Extracted single tree point cloud from the 3D hull)

此外，以2 棵相鄰激光雷達樹點云為測試對象，與文獻[2]提出的最小生成樹分割(圖8(a))和文獻[4]提出的規(guī)范割分割(圖8(b))相比，本文交互式分割方法(圖8(c))較好地避免了誤分割。

圖8 不同方法分割結果比較((a)最小生成樹分割；(b)規(guī)范化分割；(c)交互式分割) Fig.8 Comparison of various segmentation results from different methods ((a) Minimum spanning tree segmentation; (b) Normalized-cut;(c) Interactive segmentation)

4.2 交互式樹木建模結果

圖9(a)為一棵ALS 樹點云，實驗中固定交互式調節(jié)約束角度、刪除閾值和影響半徑中2 個參數，僅改變其中一個參數，重建結果表明：①在刪除閾值和影響半徑不變的情況下，約束角度從整體上抑制了樹枝的走向，在無角度約束時，部分枝條會出現向下生長和枝條開角不合理的情況，約束角度θp越大，枝條開角隨之增大，生成的枝條數量越多(圖9(b,c))，不合理角度的枝條被修剪，但約束角度不能無限增大，當約束角度大于其實際所有枝條中的最大角度，約束便不再有意義，故找調整樹木對應的最佳生長約束角度對樹木的幾何重建非常重要，經大量實驗證明，本系統(tǒng)測試的樹點云最佳生長約束角度在80°～100°之間；②在約束角度和影響半徑不變的情況下，刪除閾值dk主要影響細枝的生成數量，如圖9(c,d)所示。

圖9 改進的空間殖民算法中約束角度θp、刪除閾值dk 和影響半徑rp 對建模結果的影響((a)輸入點云；(b) θp=85°，dk=2d，rp=20d；(c) θp=95°，dk=2d，rp=20d；(d) θp=95°，dk=3d，rp=20d；(e) θp=85°，dk=2d，rp=10d)Fig.9 The impact of the different combinations of constraint angle,kill distance and influence radius in the improved Space Colonization Algorithm (SCA) ((a) Input point cloud;(b) θp=85°,dk=2d,rp=20d;(c) θp=95°,dk=2d,rp=20d; (d) θp=95°,dk=3d,rp=20d;(e) θp=85°,dk=2d,rp=10d)

為檢驗采用交互式操作是否能找到較好保留原始樹枝點云特征的參數組合，以一棵裸眼能觀察到部分樹枝結構的樹點云為例(圖10(a)～(c))，從藍色區(qū)域重建結果可知，通過設置較小的刪除閾值和影響半徑，可重建出與原始樹枝點云更加吻合的樹枝幾何模型。

圖10 不同參數組合下重建結果與原始樹枝點云對比((a)原始點云；(b) θp=90°，dk=2d，rp=20d；(c) θp=90°，dk=1.8d，rp=12d)Fig.10 Reconstructed results of branch details with different combinations of parameters in SCA ((a) Input point cloud; (b) θp=90°,dk=2d,rp=20d;(c) θp=90°,dk=1.8d,rp=12d)

為驗證改進算法的魯棒性，測試了如圖11 所示的2 組殘缺樹點云，重建結果表明，雖然缺失大量點云信息，改進的空間殖民算法仍然能跨過缺失點云區(qū)域，穩(wěn)定生成較自然的樹枝形態(tài)。進一步研究發(fā)現，即使對于明顯無樹形特征的“Stanford Bunny”及“Happy Buddha”點云，通過交互式調整合適參數，亦能重建出兼具點云特征和樹形特征的樹模型(圖12)。

圖11 殘缺點云重建結果 Fig.11 Reconstructed results of two incomplete point clouds with different parameters ((a) θp=90°,dk=2d,rp=17d;(b) θp=90°,dk=2d,rp=20d)

圖12 “Stanford Bunny”及“Happy Buddha”重建結果 Fig.12 Reconstructed results of“Stanford Bunny”and“Happy Buddha”point clouds with different parameters ((a) θp=90°,dk=1.2d,rp=18d;(b) θp=90°,dk=1.2d,rp=7d)

圖13(a)和(b)為航拍真實樹木照片，圖13(c)～ (f)為對應ALS 樹木點云，采用本文交互式重建算法，可生成與對應點云和真實照片接近的樹木幾何模型(圖13(e)和(g))。但是，樹木幾何模型重建質量與樹點云的質量緊密相關，因為采集的樹點云密度低，所以重建結果與真實照片相比在細節(jié)上依然有較大差異。

圖13 重建結果與真實航拍樹木圖片對比((a)航拍區(qū)域及對應(b)單棵樹木照片；(c)與區(qū)域(a)對應的機載激光雷達掃描點云；(d)單棵樹點云；(e)采用本文算法對樹點云(d)重建的樹枝結構俯視圖；(f)樹點云(d)對應的側視圖；(g)重建結果) Fig.13 Comparison of reconstructed tree model and aerial photograph and point cloud ((a) The area of photograph with (b) A single tree;(c) The corresponding point clouds of the area (a);(d) The top view of the single tree point cloud;(e) Top view of the reconstructed tree model;(f) Side view of the single tree point cloud (d);(g) The reconstructed tree model)

圖14(上)展示了6 棵形狀不同、點云數量不同(點云數量對應為6 534，4 809，8 407，10 234，2 163和4 462)、密度不同的樹點云，設置相同的重建參數(θp=90°，dk=2d，rp=20d)，采用改進的空間殖民算法重建樹枝模型并添加樹葉后的效果如圖14(下)所示，結果表明，改進的空間殖民算法能夠較好地重建出樹的自然形態(tài)。

圖14 6 種不同形態(tài)(上)樹點云對應的(下)重建結果展示 Fig.14 (Top) Six tree point clouds with various shapes and their corresponding (Bottom) reconstructed tree models

在重建效率方面，對圖11 和圖14 中對應的8棵樹點云數量及建模時間進行統(tǒng)計(圖15)，可以發(fā)現，當點云數量小于5 000 時，重建時間約在1.2 s之內，當點云數量大于等于5 000 時，重建時間明顯增長，因此后期需進一步對算法優(yōu)化，降低時間復雜度。

圖15 點云數量與重建時間關系 Fig.15 Relationship between the number of points and the reconstructed time

5 結束語

本文針對ALS 樹點云密度低、點云缺失導致的自動分割與建模困難等問題，提出了一種基于高度圖和切片技術相結合的交互式分割方法，與最小生成樹分割及規(guī)范自動分割方法對比發(fā)現，手動交互式分割有效減少了樹木重疊部分對分割結果的影響，避免過分割問題，且對于單一角度難以識別的復雜點云，可借助多角度調整實現精細分割。

在完成單棵樹點云的分割之后，基于改進的空間殖民算法，通過交互式調節(jié)約束角度、刪除閾值和影響半徑等3 個參數的不同組合，可實現樹枝茂密程度、光滑度以及和原始點云吻合度的調整，從而最終決定樹木的生長形態(tài)，此外對于殘缺的稀疏樹點云，該算法亦能較好地重建出與點云吻合度高的自然樹模型。

由于內在生物的多樣性和外在生長環(huán)境的差異性，自然界樹的形態(tài)千差萬別，很難找到一種通用的方法精確描述自然界中所有的樹種，本文通過實驗發(fā)現，幾乎沒有任何一組參數，可以讓所有樹木點云都能得到完美重建。但對于同一樹種，其外在的差異性相對較小，今后，一方面將基于學習的方法研究不同樹種的參數取值空間，例如根據樹種的不同自動選取相應的參數組合，提高重建自動化效果；另一方面，將結合kd樹和并行計算優(yōu)化空間殖民算法，提高重建效率。最后，交互式分割存在費時、費力等問題，今后擬結合基于區(qū)域增長或規(guī)范的分割方法，在二維切片上先進行自動預分割，然后在預分割基礎上進行交互式調整，實現半自動分割，提高操作效率。