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

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

1 引 言

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

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

2 單株樹骨架提取算法

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

圖1 樹木骨架提取流程

2.1 模糊c-均值聚類

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

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

對數(shù)據(jù)集X進(jìn)行模糊分區(qū)，主要算法如下：

(1)

(2)

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

(3)

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

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

2.2 B樣條曲線擬合

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

(4)

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

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

2.3 法矢變化

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

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

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

圖4 法矢變化及節(jié)點識別

3 試驗及其結(jié)果分析

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

3.1.1 數(shù)據(jù)獲取及處理

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

圖5 地面LiDAR樹木點云圖

圖6 抽稀后的樹木點云圖

圖7 聚類后的樹木點云圖

3.1.2 骨架提取

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

圖8 聚類團(tuán)質(zhì)心

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

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

圖10 完整的樹木骨架

3.2 試驗2：機(jī)載點云樹木骨架提取

3.2.1 數(shù)據(jù)獲取及處理

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

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

3.2.2 骨架提取

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

圖11 樹木位置、樹干識別

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

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

圖13 機(jī)載條件下提取的樹木骨架

3.3 定性分析

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

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

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

4 結(jié)束語

骨架提取在圖形可視化及建模的應(yīng)用中是一個比較基本的問題，它能夠反映物體的結(jié)構(gòu)特征及拓?fù)湫畔ⅰＯ∈椟c云數(shù)據(jù)的精度較低，存在樹木數(shù)據(jù)結(jié)構(gòu)缺失的情況，選擇一種方法最大限度地提取出樹木骨架信息至關(guān)重要。

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

參考文獻(xiàn)：

[1] 程章林，陳寶權(quán).大規(guī)模樹木三維建模研究[J].先進(jìn)技術(shù)研究通報，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維點云的拓?fù)浣Y(jié)構(gòu)表征與計算技術(shù)[J].中國圖象圖形學(xué)報，2008，13(8)：1576-1587.

[7] 黃文偉.基于拉普拉斯算子的點云骨架提取[D].大連：大連理工大學(xué)，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] 劉春，陳華云，吳杭彬.激光三維遙感的數(shù)據(jù)處理與特征提取[M].北京：科學(xué)出版社，2010：8-14.

[12] 徐祖艦，王滋政，陽鋒.機(jī)載激光雷達(dá)測量技術(shù)及工程應(yīng)用實踐[M].武漢：武漢大學(xué)出版社，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] 劉雪梅，莊晉林，張樹生，等.利用自適應(yīng)模糊橢球聚類實現(xiàn)點云分區(qū)[J].計算機(jī)工程與應(yīng)用，2007，43(15)：33-35.

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

[20] 成思源，謝韶旺.Geomagic Studio 逆向工程技術(shù)及應(yīng)用[M].北京：清華大學(xué)出版社，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] 石銀濤，程效軍，張鴻飛.基于參數(shù)L-系統(tǒng)的三維樹木仿真[J].同濟(jì)大學(xué)學(xué)報(自然科學(xué)版)，2011，39(12)：1871-1876.