一種深度學習網(wǎng)絡的樹木點云骨架重建方法

徐風， 張博， 袁星月， 安鋒， 陳幫乾， 周立軍， 云挺*

(1.南京林業(yè)大學信息科學技術學院， 南京 210037； 2.南京林業(yè)大學機械電子工程學院， 南京 210037； 3.中國熱帶農業(yè)科學院橡膠研究所， 海口 570000)

樹木骨架空間模型的精準重建在研究林木資源信息化和表型結構特征反演中發(fā)揮著至關重要的作用[1-2]。近些年，針對樹木模型重建的技術大致可以分為兩類，分別是基于圖像與激光點云的樹木重建。這兩種技術在數(shù)據(jù)獲取上存在著差異性，但提取樹木空間特征的方式上又具有一定的相似性。

基于圖像的樹木枝干重建方法通常采用雙目視覺[3]、多目視覺[4]及深度相機的方法捕獲樹木的紋理與景深信息[5]，結合曲率約束實現(xiàn)樹木枝干重建，或利用多目視覺的多張圖像通過交互式編輯重建樹木枝干[6]。基于圖像的樹干重建方法雖然簡單易行，但提取的樹木枝干停留在二維層面，很難直觀反映樹木的三維幾何形態(tài)，且提取結果會被多樣的環(huán)境背景、復雜的樹干空間拓撲結構、立體匹配、誤差、視角遮擋等多方面因素干擾[7]。

激光掃描具有高精度、稠密度等特性[8]，同時，可結合計算機圖形學、機器視覺等理論算法開展林木三維點云的處理分析，例如，利用圖論將采集到的點云數(shù)據(jù)存儲為八叉樹結構[9]，采用結合點法線的PROSAC算法建模提取林分樹干[10]；或使用空間填補的三角面片構造枝干表面來建立樹木冠層模型[11]；或通過建立方向場引導并結合圖連通算法從激光點云重建樹木骨架[12]；其他還有利用提取枝干特征點構建模型的方式，主要有：采用水平切片[13]、超體素分割[14]、基于矢量的枝干軸線提取模型[15]和基于拉普拉斯的骨架收縮方法[16]，從稠密的掃描點云中找尋局部的方向性[17]、結構性[18]、與不同的空間特征[16]也是樹木枝干重建的要點。與基于二維圖像的樹木枝干重建方式相比[7]，基于激光點云數(shù)據(jù)的樹木建模得到的樹木幾何精度高，普遍能夠較好地呈現(xiàn)出樹木的空間真實狀態(tài)。但也存在一些問題，例如，樹木具有多個一級二級分枝，冠層內部拓撲結構復雜，不能精準分類葉子和樹木枝干點云分類精度不高等，這些因素對樹木骨架后期重建會產(chǎn)生較多干擾。

近年來，面向點云處理分析的深度學習網(wǎng)絡逐漸被學者們關注，如PointNet[19]、VoxelNetVOX[20]、PointCNN[21]、分層K-d樹ContextNet[22]和SPLATNet[23]。前幾種網(wǎng)絡目前主要應用于工業(yè)中產(chǎn)品和規(guī)則幾何物體的檢測，以及車輛自動駕駛道路上目標識別，而應用于復雜拓撲結構的林木點云中的器官識別并不多見。同時，目前深度學習算法處理復雜點云還有很多未解決的問題，主要集中在3個方面。

(1)深度學習既需要考慮識別物體的全局特征，又需要融入物體的局部特征，這涉及局部鄰域基團的特征計算與對象整體特征表達在深度學習網(wǎng)絡中融合的問題。

(2)空間點云由排列散亂且無序的點構成，難以進行規(guī)則的卷積運算。如何面向空間點云開展有效的卷積運算并且提取特征，是待研究的問題。

(3)多尺度分解可以整合目標對象的整體與細節(jié)特征，如何有效地將多尺度分解應用到點云的處理分析中是研究的方向。

鑒于以上種種因素，現(xiàn)運用晶格投影策略的深度學習網(wǎng)絡與計算機圖形學算法，處理地面激光樹木點云數(shù)據(jù)。該深度學習網(wǎng)絡構造點云鄰域基團，通過旋轉不變模塊與多尺度晶格投影得到該點在不同尺度下的二維特征，以便于卷積運算，并結合深度學習的池化與多層感知機操作，自動實現(xiàn)復雜林木點云的枝葉分類操作。以海南多類樹木為研究對象，開展精準的枝葉分離并重建樹木空間枝干的三維模型。同時，對不同級別枝干和相應的葉團簇進行參數(shù)提取，并與真實測量數(shù)據(jù)比對以驗證算法的有效性。

1 激光點云數(shù)據(jù)獲取與預處理

1.1 研究區(qū)與樣本概況

研究區(qū)位于海南大學校園(20°3′N，110°20′E)與儋州市橡膠種植園(19°32′N，109°28′E)，試驗地地理位置在低緯度熱帶北緣，屬于熱帶海洋氣候，春季溫暖少雨多旱，夏季高溫多雨，秋季多臺風暴雨，全年日照時間長，輻射能量大，年平均氣溫23.8 ℃，最高平均氣溫28 ℃左右，最低平均氣溫18 ℃左右，年平均降水量1 664 mm。選取研究區(qū)內多株不同樹種作為實驗對象，如橡膠樹(品種為PR107和CATAS 7-20-59)、香樟樹、紅楓樹、櫻花樹等。

1.2 數(shù)據(jù)獲取與預處理

點云數(shù)據(jù)的獲取使用Leica Scanstation C10 型號三維激光掃描儀，其掃描角度為360°× 270°，掃描速率為5萬點/s，各類樹木的數(shù)據(jù)采集時間為2019年8月15日。考慮到樹木的冠層內部結構復雜、樹葉相互遮擋[24-25]，為了獲取完整的樹木三維點云數(shù)據(jù)，以每棵樹冠為中心，分別對其進行對稱的2站掃描，每站掃描儀距掃描樹冠的中心3 m遠，掃描精度設置為中等精度，具體掃描參數(shù)為：角度精度(水平/垂直)：中等0.057°/0.057°、視場角(水平/垂直)：360°/270°、掃描點間距：距離3 m時最小<3 mm、標靶掃描精度：3～5 mm。接著采用 Cyclone 軟件結合標靶數(shù)據(jù)對原始兩站點云進行了手動拼接，從而得到完整的目標樹葉子與枝干點云數(shù)據(jù)，掃描的點云密度平均為12 000 點/m2。由于獲取到的樹木三維點云數(shù)據(jù)中含有噪點，采用文獻[26]的方法對采集的點云進行了去噪處理，接著把標記好枝葉信息的點云數(shù)據(jù)帶入晶格投影的深度學習網(wǎng)絡中訓練。

2 實驗方法

2.1 晶格投影的深度學習網(wǎng)絡實現(xiàn)樹木枝葉分離

圖1 枝葉分類的晶格投影的深度學習網(wǎng)絡Fig.1 Deep learning network for the classification of wood and leaf point clouds

2.1.1 旋轉不變性模塊

構造的旋轉不變性模塊的具體如下：首先將點云數(shù)據(jù)看作一個整體，分步驟旋轉，運用三個旋轉矩陣：

(1)

(2)

(3)

式中：RX、RY、RZ分別為繞X、Y、Z軸的旋轉矩陣；α、β、γ分別為繞各自軸的旋轉角度，輸入以pi為中心點，直徑為0.1 m呈立方體領域的點云基團，該點云的坐標為n×3的矩陣pj(xj，yj，zj)，j=1，2，…，n，矩陣pj經(jīng)旋轉變換后得到新坐標矩陣：

p′j=pjRXRYRZ

(4)

旋轉變換結果如圖1(b)所示，一個不規(guī)則放置的局部枝干和葉子經(jīng)旋轉變換后，枝干的軸向幾乎與Z軸平行(即直立)，且葉子平行于XOZ平面。

2.1.2 晶格投影與重心插值

將旋轉后的樹木局部點云p′j作為輸入置于一個立方體之中。該立方體的XOY平面、XOZ平面和YOZ平面上都排布有晶格頂點，單個平面上的頂點有共σ2個，呈σ×σ二維矩陣狀規(guī)則排列。這里σ為晶格尺度，它是一個常數(shù)，且三層晶格投影模塊的σ不相同，是逐步減小的。定義第一層的σ為10，那么第二層的σ為「10/2?，第三層的σ為「10/4?，如圖1 (a)所示。「?指向上取整。晶格的大小是通過晶格尺度σ來控制的，σ的變化會改變晶格頂點之間的距離，可以類比為對晶格進行等比例放大的操作。σ越大，晶格尺寸越大，點云細節(jié)表現(xiàn)更多。

(5)

式(5)中：mod為求余操作；{}內左邊為被除數(shù)；右邊為除數(shù)；「?指向上取整；?」指向下取整。

晶格投影完成后，使用重心插值將對投影生成的點進行濺射操作，將一個點的數(shù)據(jù)散開至最近的三個晶格頂點處。

(6)

式(6)中：只有bj，0、bj，1和bj，2三個未知數(shù)，解得插值系數(shù)為

(7)

以一個局部枝干和一片局部樹葉為例，上述操作結束后能夠得到如圖2和圖3所示結果。圖2(d)與圖3(d)是使用到的樹葉局部點云圖與枝干局部點云圖，圖2(a)～圖2(c)是樹葉點云在不同晶格尺度與不同平面的晶格投影與重心插值結果，圖3(a)～圖3(c)是枝干點云在不同晶格尺度與不同平面的晶格投影與重心插值結果。各圖像對應的晶格尺度σ已標注在分圖題中，對應的投影平面已標注在各圖像左上方。三條顏色條分別對應10、5和3這三種晶格尺度下，晶格投影與重心插值結束后，晶格點所累加的插值系數(shù)大小。

將所有點在同一個平面上生成的重心插值系數(shù)存儲在同一個矩陣內，因為投影到三個平面，且平面數(shù)量與存儲插值系數(shù)的矩陣數(shù)量相同。所以存儲插值系數(shù)的矩陣編號k=0、1、2，共有三個存儲插值系數(shù)的矩陣，每個矩陣記為

圖2 樹葉點云的晶格投影圖Fig.2 Schematic diagram showing the lattice projection for leaf point clouds

圖3 枝干點云的晶格投影圖Fig.3 Schematic diagram showing the lattice projection for branch point clouds

(8)

2.1.3 激活函數(shù)與數(shù)據(jù)集

生成的輸出數(shù)據(jù)根據(jù)使用到的晶格尺度的不同，對應使用不同的多層感知器進行卷積。具體來說：晶格尺度為10時，得到的3×n×3的輸出數(shù)據(jù)全部輸入卷積模板為MLP(3， 16)的多層感知機中，卷積后得到3×n×16的數(shù)據(jù)。晶格尺度為5時，得到的3×n×3的輸出數(shù)據(jù)全部輸入卷積模板為MLP(3， 32)的多層感知機中，卷積后得到3×n×32的數(shù)據(jù)。晶格尺度為3時，得到的3×n×3的輸出數(shù)據(jù)全部輸入卷積模板為MLP(3， 64)的多層感知機中，卷積后得到3×n×64的數(shù)據(jù)。接著鏈接三組輸出數(shù)據(jù)，得到一個3×n×112的特征，進行累加，得到n×336個特征。再使用全局池化操作，形成一個336的向量，最后使用多層感知機MLP(64，32，2)與激活函數(shù)，得到2×1的枝和葉置信度值。

枝葉分類問題是一個二分類問題，為了處理該問題，深度學習網(wǎng)絡在最后使用sigmoid作為激活函數(shù)，計算公式為

(9)

sigmoid函數(shù)將網(wǎng)絡輸出的2×1的枝和葉置信度值轉化至0～1，并取較大值代表對應輸入點云的類別。

數(shù)據(jù)集由4種樹木組成，分別是橡膠樹(CATAS 7-20-59)、橡膠樹(PR 107)、香樟樹和紅楓樹。橡膠樹是熱帶樹種，具有很大的經(jīng)濟價值， CATAS 7-20-59與PR 107相比，較為低矮。香樟樹是亞熱帶樹種，樹形高大，多用作綠化與美化景觀。紅楓是亞熱帶樹種，樹形優(yōu)美，多用作觀賞植物。數(shù)據(jù)集樹木參數(shù)如表1所示。

使用機器學習[27]，輔助人工標注對數(shù)據(jù)集進行了枝葉分離，給每個點標好標簽以便深度學習網(wǎng)絡訓練時使用。部分枝葉分離后的訓練集如圖4所示。

表1 深度學習數(shù)據(jù)集樹木參數(shù)Table 1 The tree parameters of the deep learning dataset

棕色為分割后的枝干點云; 綠色為分割后的葉子點云圖4 部分深度學習訓練集 Fig.4 Partial training sets for our deep learning network

2.2 樹木骨架重建

2.2.1 樹干高度分層與中心點求取

圖5 部分橡膠樹干點云及每層聚類中心點提取Fig.5 Part of the scanned rubber tree branch point cloud and the corresponding extracted central points of each layer

2.2.2 枝干鏈表確立及圓柱體擬合

圖6 樹木骨架邊緣點與分枝節(jié)點提取Fig.6 Extraction of the bifurcation and edge points from each tree skeleton

直到該鏈表遇到邊緣結點時，從根結點開始的鏈的中心點被分類為主枝干結點。同時，在主枝干鏈上延伸的分枝結點稱為一級分枝結點，進而定位整個樹干的其他一級分枝。

(10)

(11)

式(11)中：q1、q2分別為擬合直線L1上的任意2個點，‖‖2為2范數(shù)。

重復上述過程，即從當前枝干段上隨機采樣兩個點，并不斷優(yōu)化得到d的最小值，直至d不再變化。

圖7 生長角度變化最小準則確立主枝干鏈表Fig.7 Tree trunk chain determination using the rule of minimal variation in the growth angle

圖8 自適應半徑計算的圓柱體分段枝干點云擬合Fig.8 Adaptive radius calculation of cylindrical subsection branch point cloud fitting

2.2.3 一級枝干對應葉團簇分割

最后，以分離的主枝干和各一級分枝掃描點云為聚類中心，根據(jù)空間分水嶺聚類算法[28]，實現(xiàn)不同枝干所對應的葉子點云分類，并運用Alphashape算法[30]計算單株樹不同分類點云對應的空間葉團簇的體積。

3 結果與分析

3.1 枝葉分離與骨架重建

深度學習網(wǎng)絡運算在英特爾i7-7700 CPU@2.80 GHz處理器和16 GB內存的Windows 10 64位PC上執(zhí)行，在PyCharm和MATLAB軟件平臺上執(zhí)行了晶格投影策略的深度學習網(wǎng)絡和樹木重建程序，同時使用NVIDIA RTX 2080Ti GPU代替CPU來減少訓練時間。在晶格投影策略的深度學習網(wǎng)絡模型中，學習率為0.000 1，批量大小設置為16，迭代次數(shù)為300。設置的總訓練時間約為80 h。

具體訓練過程中的準確率和損失函數(shù)值如圖9所示，隨著學習迭代次數(shù)的不斷增加，訓練樣本的分類準確率呈上升趨勢，而損失函數(shù)呈下降趨勢，這表明所使用的深度學習模型參數(shù)滿足全局優(yōu)化且收斂。同時，訓練期間神經(jīng)元網(wǎng)絡在每個批次中遇到了一些復雜的樣本，如點云存在局部遮擋和數(shù)據(jù)缺失、一些葉片掃描不全呈現(xiàn)枝干形態(tài)而被當前的深度網(wǎng)絡權值進行誤判等，進而導致回歸損失函數(shù)值的局部波動。

在經(jīng)歷了100個Epoch之后，準確性和損失函數(shù)值分別到達0.90和0.05，表明模型參數(shù)較優(yōu)。通過晶格投影的深度學習網(wǎng)絡的分割結果如圖10所示。

綠色代表葉子點云；棕色代表枝干點云圖10 不同樹木枝葉分割結果圖Fig.10 Wood-leaf classification results for scanned points of different tree species

深度學習方法與傳統(tǒng)機器視覺算法枝葉分類性能對比如表2所示，分別在分類精度、IOU與分類時間三個參量上闡述算法的優(yōu)劣。其中IOU定義為

(12)

面向不同樹種枝干數(shù)據(jù)，根據(jù)每層高度的中心點與枝干掃描點云進行圓柱體擬合以構造樹木骨架模型，部分結果如圖11所示。

表2 機器學習與深度學習的點云枝葉分類性能表Table 2 Performance evaluation of the wood-leaf point classification using a deep learning approach versus the machine learning method

3.2 不同分枝與葉團簇提取

面向不同樹種，運用本文方法提取的不同分枝及對應的葉團簇如圖12和13所示。提取的不同一級枝干顯示為不同顏色，而用不同顏色的最小凸包來展示相應分割的葉團簇。

3.3 樹木參數(shù)反演結果比較

表3羅列了單株紅楓樹和香樟樹以及不同品種的多棵橡膠樹的生長參數(shù)實際測量值和本文方法的計算值，其中包括樹高、冠積、葉團簇體積、胸徑、一級枝干、二級枝干直徑和主枝干與一級枝干之間夾角，同時，把實地測量值與本算法計算值進行對比驗證，定量化表達本文方法的有效性。

圖14給出了運用本文方法與實地驗證的橡膠樹比較結果。圖14(a)顯示了4棵PR107基因型和4棵CATAS 7-20-59基因型橡膠樹計算的一級枝干直徑與實測值的比較結果，具體為PR107 (決定系數(shù)R2=0.93， 均方根誤差RMSE=0.77 cm， 相對均方根誤差rRMSE=5.75%)，CATAS 7-20-59 (R2=0.92， RMSE=0.51 cm， rRMSE=7.21%)。 圖14(b)顯示PR107基因型和CATAS 7-20-59基因型橡膠樹主枝干與一級枝干夾角的計算與測量比對結果，可以看出PR107型橡膠樹具有較大的分枝角，與實測值比對為(R2=0.92， RMSE=4.86°， rRMSE=7.56%)；而CATAS 7-20-59型橡膠樹的分枝角相對較小，反演結果為(R2=0.91， RMSE=1.77°， rRMSE=8.06%)。實驗整體結果表明算法在樹木枝干直徑與分枝角度估算上具有較高的精度。圖14(c)顯示了兩種橡膠樹一級枝干直徑和對應的葉團簇體積呈現(xiàn)正相關性，其中PR107型橡膠樹具有更多的一級枝干和對應的葉團簇，整體樹冠結構為發(fā)散形，而CATAS 7-20-59具有較少的一級枝干，樹冠結構為倒花瓶形。同時，圖14(c)表明較粗的枝干可以支撐起較大的葉團簇，這與更多葉花果等器官需要更粗的枝干進行營養(yǎng)傳輸和重力支撐的原理相吻合。

不同顏色代表不同枝干分段圓柱體擬合結果圖11 圓柱體分段擬合的樹木骨架重建結果圖Fig.11 Skeleton reconstruction results of tree species using many cylinder models to fit each segment of tree branch points

不同顏色對應不同一級分枝圖12 不同樹木一級枝干提取結果圖Fig.12 Every first-order branch extraction results using our approach for different tree species

棕色代表枝干點云其他顏色代表分割后的葉團簇圖13 一級分枝對應的葉團簇分割結果圖Fig.13 Foliage clump borne on the first-order branch of the same tree was segmented

表3 本文方法獲得的林分參數(shù)與實地測量值對比Table 3 Tree property retrieval using our method versus field measurements

圖14 算法結果與實測數(shù)據(jù)的比對分析Fig.14 The retrieved growth properties by our method versus the field measurements

4 結論

樹木骨架的精準重建對于分析樹木的表型結構、樹體自身的物理特性、地理環(huán)境因子的影響都起到了信息支撐的作用。本文研究面向地面激光雷達數(shù)據(jù)開展多類闊葉樹骨架形態(tài)的建模，通過面向枝葉點云分類的深度學習模型設計、空間鏈表結構的主枝干與一級枝干查找、參數(shù)自適應調整的圓柱體擬合、空間分水嶺分類的葉團簇提取，精準且智能的重建了林木的枝干骨架結構與參數(shù)獲取。

本文研究設計了晶格投影策略的深度學習網(wǎng)絡，并融入了點云局部基團特征、旋轉不變模塊與多尺度晶格投影等模塊，解決了局部特征提取、空間點云卷積和多尺度特征融合等問題。

雖然深度學習網(wǎng)絡枝葉分類準確率達到91.31%，但導致準確率降低的因素主要有：實際掃描中枝葉互相遮擋，造成數(shù)據(jù)缺失；葉子和枝干點云點云緊密貼合，特征混淆而難以分離；點云訓練集與測試集中存在噪聲，影響了深度學習網(wǎng)絡最終分割性能，這些都是下一步研究主要提升的方向。