一種稀疏點云環境下的單株樹骨架提取算法

李巍岳，劉春，吳杭彬，孫偉偉

(1.同濟大學 測繪與地理信息學院，上海 200092；2.現代工程測量國家測繪地理信息局 重點實驗室，上海 200092；3.礦山空間信息技術國家測繪地理信息局重點實驗室，焦作 454000)

1 引 言

樹木三維建模一直是計算機圖形圖像學及仿真學研究的熱點之一，其模型重構技術已在虛擬現實、林業調查、生態環境建設、古樹名木保護以及農林業等方面有著廣泛的應用。多數學者主要用生物形態法及幾何構造法兩種方法來構建模型。前者主要是從植物學的角度模擬樹木形成、生長及衰敗的過程[1]，但該方法沒有考慮樹木的形態及結構信息，所以很難構建出真實樹木的三維模型；后者通過計算機圖形學的手段，參照樹木的拓撲結構信息來生成模型[2]，彌補了前者在建模上的不足，但該方法沒有嚴格遵循植物學的規律。樹木重構首先要提取樹木的骨架信息[3]，樹木的骨架信息對描述形狀具有重要的意義，不僅能夠保持樹木拓撲結構、方向及局部特征，而且可以以直觀、緊湊的形式有效地表示3D模型，是樹木建模的關鍵問題之一[4-6]。骨架提取目前沒有嚴格的概念定義，無法確定提取算法的優劣，這是因為骨架提取針對不同的研究領域，質量的評價目前還沒有一個統一的標準[7]。在樹木骨架提取的研究中，應滿足骨架的拓撲一致性、居中性及細性等特點。已有的研究[8-10]多數基于圖像及手工繪制草圖的交互式建模，提取樹木的骨架信息。這部分研究中，除了基于圖像的建模外，大部分的研究都是建立虛擬三維樹木，沒有實現真實樹木的數字化。由于許多樹木形體高大，枝條眾多，樹葉與枝干間存在遮擋，單純的利用圖像很難實現樹木準確的骨架提取。

激光掃描測距技術(Light Detection and Ranging，LiDAR)是一種快速、直接地獲取地形表面模型的技術；與傳統的光學及微波遙感不同，LiDAR能夠快速、精確地獲取地面特征在水平和垂直方向上的位置，克服了傳統觀測手段中樹木部分結構枝葉間遮擋的問題，是目前能測定森林覆蓋地區地面高程的可行技術之一[11-12]。近些年來，利用點云信息提取樹木骨架的研究逐漸增多，Pyysalo等采用矢量多邊形的方法，在點云中識別出云杉、松樹、樺樹三類樹種，分別計算點云的算數中心及各點與中心點的連線、夾角，按照逐漸增大的角度順次連接各點形成單木樹的骨架[13]；Gorte等基于點云與柵格影像，分別選擇濾波、數學形態學、最短路徑等算法，在構建的體素模型中提取3D樹木骨架[14]；Delagrange等在垂直方向上將點云按照1cm間隔進行水平切片，找出每一層數據點構成的中心，然后從下向上通過最小生成樹將這些中心點連接得到樹木骨架[15]；Bucksch等利用八叉樹數據結構來組織點云數據，通過體素模型標注物體的延伸方向，再結合圖論方法及約束規則來提取骨架[16]；Raumonen等采用圓柱體模型，在收集到的地面LiDAR數據中，快速自動地構建出單株樹的樹枝及樹干[17]。目前針對激光掃描數據提取樹木骨架的方法中，多數建立在三維網絡模型或體素化模型的基礎上。由于點云模型獲取困難，噪聲多，缺少模型的連接信息，在提取過程中，需要進行點云去噪，體素剖分等處理，過程費時，結果較難預料[18]；特別是在稀疏點云的條件下，影響的因素更多，嚴重地影響到樹木的三維重建，因此有必要尋找一種方法能夠解決上面提到的問題。本文根據樹木的生長規律，分別闡述了單株樹的主干、節點及枝條信息的獲取方法，總結出在稀疏點云(地面及機載)條件下提取樹木骨架的最優化方法。

2 單株樹骨架提取算法

樹木信息一般較為復雜，采用激光掃描技術所獲得的點云受到漫反射和其他信息的干擾，點云比較稀疏且精度有限。雖然所獲取的點云數據分布稀疏離散，但點云之間存在密度關聯與較好的層次性。鑒于此，本文主要的思想是：大多數的激光點反映的是樹木不同位置的葉子信息，葉子點云內具有團聚的特點，將每一個聚類中心取代周圍的葉子點云，作為樹木骨架的末端；還有部分點云信息反映了樹木的主干及枝條，采用自由型曲線進行模擬，可以最佳的還原樹木自身的結構，同時，控制自由型曲線走向的法向量得到了樹木枝條的節點；最終得到單株樹的骨架(圖1)。

圖1 樹木骨架提取流程

2.1 模糊c-均值聚類

葉子點云數據較為離散，不同的樣本類屬存在不確定性，模糊c-均值分類方法使用目標函數將數據分成互相分離的類別，其主要算法為：

設數據集X={x1,x2,…,xn}?Rr，每個數據樣本xi由m個特征定義，m選擇根據諸克軍等提出的遺傳-迭代自組織分析技術(GA-ISODATA)來進行優選[19]。

對數據集X進行模糊分區，主要算法如下：

(1)

(2)

vi由m個特征來表述，對于點云中每個類的中心坐標按照下面公式計算：

(3)

j是特征空間中的一個變量，即j=1，2，…，m。通過經典模糊聚類方法，可以識別出三維離散點云中密度相對集中的區域(圖2)，將離散的點組織為若干個組別。在樹木識別中，通過距離的選擇，可以初步判斷骨架點云與葉子點云的分布。葉子點云中每個類的中心能夠反映枝條的末端。

圖2 基于模糊c-均值點云聚類

2.2 B樣條曲線擬合

樹木主干及枝干結構較為復雜，單純采用直線段來模擬，不美觀也不符合自然規律，常采用帶有弧段變化的曲線來表現。同時，選取的曲線需要滿足控制樹木形態及易于實現的特點。B樣條方法能夠表示與設計物體自由型的曲線與曲面(圖3)，是逆向工程中廣泛應用的數學描述方法[20]，其曲線方程可寫為：

(4)

其中，pi(i=0，1，…，n)為控制頂點，順序連接成的折線成為B樣條控制多邊形；Ni，k(u)(i=0，1，…，n)稱為k次規范B樣條基函數[21]。

B樣條基是多項式樣條空間具有最小支承的一組基，可通過改變控制點個數設計一條曲線，不需要改變多項式的次數。鑒于上述B樣條的優勢，選用其作為主干及枝干重建曲線。

2.3 法矢變化

點云根據B樣條擬合得到的任一點p(u，w)處，構造一個垂直與樣條曲面的矢量，該矢量稱為法矢。單位法矢可以根據該點處切矢pu和pw的矢量積得到。

單位法矢在幾何造型中是不可缺少的，在本研究中主要用其方向衡量B樣條曲線的彎曲程度，點云法矢方向變化率(Δn)大，曲線波動大。主干及枝干在分支或者節點處，Δn較大；根據給定的初始閾值，可識別出法矢方向變化率較大的位置(圖4)。此外，規定法線方向在樣條曲線走向的順時針方向上。

圖3 離散點B樣條曲線擬合

圖4 法矢變化及節點識別

3 試驗及其結果分析

3.1 試驗1：地面稀疏點云樹木骨架提取

3.1.1 數據獲取及處理

試驗采用FARO公司的Focus 3D高精度三維激光掃描儀，該儀器以每秒最大976000點的速率可掃描最長為153.49m，具有三維虛擬重現、速度控制、高精確度及視野范圍大等特點[22]。首先，設計各掃描站及標靶的位置，滿足每測站之間至少有3個標靶重合，使點云數據能夠統一到相同的儀器坐標系下；然后，對各標靶進行精確掃描，統一每個標靶的標識，滿足相同的標靶在不同測站中的標識必須相同。通過這幾個原則，對實地樹木進行采集，得到點云圖(圖5)。數據中包括大量冗余，采用系統抽稀的方法，間隔50個點作為抽樣間隔，其原理是在50個樣本點中隨機選擇一個樣本數據點，圖6為壓縮后的點云數據。通過這種方法來隨機模擬稀疏點云的環境。

圖5 地面LiDAR樹木點云圖

圖6 抽稀后的樹木點云圖

圖7 聚類后的樹木點云圖

3.1.2 骨架提取

利用模糊c-均值聚類方法對抽稀后的點云數據進行分割，根據式(3)，當距離取為2cm，GA-ISODATA算法葉子點云被分為50類時，得到較為理想的分割密度(圖7，不同的顏色表示不同的聚類團)。其中，骨架點云與葉子點云通過不同的顏色初步得以區分；另外，不同葉子點云表現出一定的團聚特點。不同的聚類團由相應的枝條相串聯，根據樹木生理特點，其枝條的末端應分布在葉子點云聚類團的質心處(圖8，質心用紅色點表示)，提取每個聚類團的質心點為枝條節點，最外圍的質心數據作為枝條末端。

圖8 聚類團質心

將提取出的部分骨架點云，利用1.2節及1.3節的辦法，進行法矢方向變化率判斷及B樣條曲線擬合。由于地面LiDAR主要以近水平方式掃描樹木，樹木枝干點云存在大量的冗余，骨架起始節點選擇樹木點云最下端形成的流形曲面中心(圖9中A點)。調整曲線法向量變化的閾值，在搜索中遵循點與點之間的最短距離，選擇進行B樣條曲線擬合的點，完成該部分點云樹木骨架的提取(圖9)，在節點處(圖9中B、C、D、E、F、G、H、I點)，法矢方向變化率較大。按照同樣的方法，將部分點云樹木骨架的末端與聚類團的質心進行法矢方向變化率判斷及B樣條曲線擬合，形成完整的樹木骨架(圖10)。該方法能夠解決點云數據冗余及遮擋的問題，提取出的骨架能夠真實地反映樹木的拓撲結構。

圖9 部分點云樹木骨架提取

圖10 完整的樹木骨架

3.2 試驗2：機載點云樹木骨架提取

3.2.1 數據獲取及處理

機載LiDAR數據采用LiteMapper 5600系統獲取，其中激光掃描儀為Riegl LMS-Q560，波長為1550nm，激光脈沖的長度為3.5ns，激光脈沖發散角小于等于0.5mrad。獲取的離散點云數據采用WGS84坐標系，UTM投影北半球6度帶第50分帶。其中，點云之間的平均點間隔為0.15m，平均點密度為6.5個/m2。機載點云數據反映了物體垂直范圍的情況，在樹木信息提取中，主要可以獲取單木樹高、冠幅、林區的平均數高、體積、密度等。

通過點云數據的濾波處理，分離出地面點與非地面點；根據植被與人工地物的反射特性及空間幾何關系，利用數學形態學方法將植被及人工地物區分開；同時利用分水嶺算法，分割了樹木冠層高度模型，提取出表述單株樹的點云數據。

3.2.2 骨架提取

機載LiDAR通過較遠的距離及垂直掃描方式獲取樹木信息，其點云數據普遍稀少，且由于葉子的遮擋，缺少對枝干特征的點云數據描述。試驗中利用相關文獻[23-24]里提到的樹木位置、樹干的提取方法：通過低通濾波處理的樹木冠層高度模型中觀察局部最大點，從而得到樹木的位置；每棵樹的樹冠由平行的水平多邊形包圍，樹干由地面到樹頂的垂直矢量表示(圖11)。利用模糊c-均值對提取的樹木點云進行聚類，當距離取50cm葉子點云被分為30類得到較為理想的分割結果(圖12)。

圖11 樹木位置、樹干識別

圖12 聚類后的樹木點云分布圖

選擇離樹干最近的4個點A、B、C、D作為B樣條曲線模擬的起始點，參照1.3節的方法追蹤法矢方向變化率，在搜索中遵循點與點之間的最短距離，與聚類團的質心擬合得到樹木一級枝條；然后從一級枝條的下一個節點再進行法矢方向變化率及最短距離搜索，與聚類團的質心擬合次一級枝條，最終遍歷整個點云數據，得到了機載條件下的樹木骨架(圖13)。

圖13 機載條件下提取的樹木骨架

3.3 定性分析

通過地面及機載方式獲取的樹木信息方式不同，前者以近水平的方式掃描樹木，激光波長較小，樹木的骨架信息比較完整，但存在數據冗余及部分點云信息遮擋的問題；后者以近垂直的方式掃描樹木，激光波長較大，主要可以反映樹木的高度、冠幅信息，骨架信息嚴重缺失，影響到后期的樹木建模。

樹木骨架提取目前較為成熟的方法是參數化L-系統[25]，傳統的參數化L-系統是根據樹木的類型選定合理的公理和產生式，設計不同的初始狀態與樹木描述規則，通過不斷地迭代產生一系列字符串，然后逐一讀取字符，按照執行規則的不同逐一對讀取的字符串進行解釋。針對機載LiDAR樹木建模，Pyysalo[13]等在2002年提出了一種構建樹木矢量模型的方法，依次按照樹冠之間的1m間隔將整個點云分成若干層，在每個高度層內的點按遞增的方位角排列，從而得到樹木的水平多邊形，該方法已經投入到了大范圍林木林地模型的生產中。將本文提出的方法與這兩種方法進行比較，結果如表1所示。本文的方法在算法上簡單、易于實現；支持地面及機載LiDAR數據的樹木骨架提取，骨架具有一定的彎曲度，接近真實樹木的情況；適用于壓縮后的數據，可以為單木樹及大范圍樹木建模提供一定的理論基礎。

表1 不同樹木骨架提取方法比較

4 結束語

骨架提取在圖形可視化及建模的應用中是一個比較基本的問題，它能夠反映物體的結構特征及拓撲信息。稀疏點云數據的精度較低，存在樹木數據結構缺失的情況，選擇一種方法最大限度地提取出樹木骨架信息至關重要。

本文首先基于模糊c-均值聚類的方法，對點云數據進行分割，根據樹木的枝葉分布特性，選取聚類團的質心擬合枝干形狀的自由型結構，并通過曲線的法矢變化率判斷樹木的節點，最終實現了樹木骨架的提取，并分別通過地面及機載LiDAR數據的試驗，驗證了該方法適用于稀疏點云條件。

參考文獻：

[1] 程章林，陳寶權.大規模樹木三維建模研究[J].先進技術研究通報，2009，3(6)：40-44.

[2] SU Z X，ZHAO Y D，ZHAO C J，et al.Skeleton extraction for tree models[J].Mathematical and Computer Modelling，2011，54(3-4)，1115-1120.

[3] LIVNY Y，YAN F，OLSON M，et al.Automatic reconstruction of tree skeletal structures from point clouds[J].ACM Transactions on Graphics.2010，29(6)：1-8.

[4] AU O K C，TAI C L，CHU H K，et al.Skeleton extraction by mesh contraction[J].ACM Transactions on Graphics，2008，27(3)：1-10.

[5] HENNING J G，RADTKE P J.Detailed stem measurements of standing trees from ground-based scanning LiDAR[J].Forest Science，2006，52(1)：67-80.

[6] 張義寬，張曉鵬，查紅彬，等.3維點云的拓撲結構表征與計算技術[J].中國圖象圖形學報，2008，13(8)：1576-1587.

[7] 黃文偉.基于拉普拉斯算子的點云骨架提取[D].大連：大連理工大學，2009：2-4.

[8] TENG C H，CHEN Y S，HSU W H.Constructing a 3D trunk model from two images[J].Graphical Models，2007，69(1)：33-56.

[9] CHENG Z L，ZHANG X P，THIERRY F.Tree skeleton extraction from a single range image[C].The 2ndInternational Symposium on Plant Growth Modeling，Simulation，Visualization and Application，2006，274-281.

[10] OKABE M，OWADA S，IGARASHI T.Interactive design of botanical trees using freehand sketches and example-based editing[J].Computer Graphics Forum，Proceedings of Eurographics，2005，24(3)：487-496.

[11] 劉春，陳華云，吳杭彬.激光三維遙感的數據處理與特征提取[M].北京：科學出版社，2010：8-14.

[12] 徐祖艦，王滋政，陽鋒.機載激光雷達測量技術及工程應用實踐[M].武漢：武漢大學出版社，2009：2-6.

[13] PYYSALO U，HYYPPA，H.Reconstruction tree crowns from laser scanner data for feature extraction[R].International Society for Photogrammetry and Remote Sensing-ISPRS Commission III Symposium，Graz，Austria，2002.

[14] GORTE B，PFEIFER N.3D image processing to reconstruct trees from laser scans[R].Proceedings of the 10thAnnual Conference of the Advanced School for Computing and Imaging(ASCI)，Ouddorp，the Netherlands，2004.

[15] DELAGRANGE S，ROCHON P.Reconstruction and analysis of a deciduous sapling using digital photographs or terrestrial-LiDAR technology[J].Annual Botany，2011，108(6)：991-1000.

[16] BUCKSCH A，LINDENBERGH R，MENENTI M.Skeltre：Robust skeleton extraction from imperfect point clouds[J].The Visual Computer，2010，26(10)：1283-1300.

[17] RAUMONEN P，KAASALAINEN M，AKERBLOM M，et al.Fast automatic precision tree models from terrestrial laser scanner data[J].Remote Sensing，2013，(5)：491-520.

[18] 劉雪梅，莊晉林，張樹生，等.利用自適應模糊橢球聚類實現點云分區[J].計算機工程與應用，2007，43(15)：33-35.

[19] 諸克軍，蘇順華，黎金玲.模糊C均值中的最優聚類與最佳聚類數[J].系統工程理論與實踐，2005，(3)：52-61.

[20] 成思源，謝韶旺.Geomagic Studio 逆向工程技術及應用[M].北京：清華大學出版社，2010：19-21.

[21] WU Z K，ZHOU M Q，WANG X C.Interactive modeling of 3D tree with ball B-spline curves[J].The International Journal of Virtual Reality，2009，8(2)：101-107.

[22] AZMY S N，SAH S A，SHAFIE N J，et al.Counting in the dark：Non-intrusive laser scanning for population counting and identifying roosting bats[J].Scientific Reports，2010，(2)：524.

[23] YU X W，HYYPA J，HARRI K，et al.Automatic detection of harvested trees and determination of forest growth using airborne laser scanning[J].Remote Sensing of Environment，2004，90(4)：451-462.

[24] HYYPPA J，INKINEN M.Detecting and estimating attributes for single trees using laser scanner[J].The Photogrammetric Journal of Finland，1999，16：27-42.

[25] 石銀濤，程效軍，張鴻飛.基于參數L-系統的三維樹木仿真[J].同濟大學學報(自然科學版)，2011，39(12)：1871-1876.