李巍岳,劉春,吳杭彬,孫偉偉
(1.同濟(jì)大學(xué) 測繪與地理信息學(xué)院,上海 200092;2.現(xiàn)代工程測量國家測繪地理信息局 重點實驗室,上海 200092;3.礦山空間信息技術(shù)國家測繪地理信息局重點實驗室,焦作 454000)
樹木三維建模一直是計算機(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)化方法。
樹木信息一般較為復(fù)雜,采用激光掃描技術(shù)所獲得的點云受到漫反射和其他信息的干擾,點云比較稀疏且精度有限。雖然所獲取的點云數(shù)據(jù)分布稀疏離散,但點云之間存在密度關(guān)聯(lián)與較好的層次性。鑒于此,本文主要的思想是:大多數(shù)的激光點反映的是樹木不同位置的葉子信息,葉子點云內(nèi)具有團(tuán)聚的特點,將每一個聚類中心取代周圍的葉子點云,作為樹木骨架的末端;還有部分點云信息反映了樹木的主干及枝條,采用自由型曲線進(jìn)行模擬,可以最佳的還原樹木自身的結(jié)構(gòu),同時,控制自由型曲線走向的法向量得到了樹木枝條的節(jié)點;最終得到單株樹的骨架(圖1)。
圖1 樹木骨架提取流程
葉子點云數(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-均值點云聚類
樹木主干及枝干結(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)勢,選用其作為主干及枝干重建曲線。
點云根據(jù)B樣條擬合得到的任一點p(u,w)處,構(gòu)造一個垂直與樣條曲面的矢量,該矢量稱為法矢。單位法矢可以根據(jù)該點處切矢pu和pw的矢量積得到。
單位法矢在幾何造型中是不可缺少的,在本研究中主要用其方向衡量B樣條曲線的彎曲程度,點云法矢方向變化率(Δn)大,曲線波動大。主干及枝干在分支或者節(jié)點處,Δn較大;根據(jù)給定的初始閾值,可識別出法矢方向變化率較大的位置(圖4)。此外,規(guī)定法線方向在樣條曲線走向的順時針方向上。
圖3 離散點B樣條曲線擬合
圖4 法矢變化及節(jié)點識別
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.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ī)載條件下提取的樹木骨架
通過地面及機(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 不同樹木骨架提取方法比較
骨架提取在圖形可視化及建模的應(yīng)用中是一個比較基本的問題,它能夠反映物體的結(jié)構(gòu)特征及拓?fù)湫畔ⅰO∈椟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.