樹木枝干點云數(shù)據(jù)的等值線提取

高士增，張懷清，劉 閩，何清平，羅立平

（1.中國林業(yè)科學(xué)研究院 資源信息研究所，北京 100091；2.湖南省攸縣黃豐橋林場 廣黃分場，湖南攸縣 412306）

隨著森林可視化和林業(yè)測繪技術(shù)的不斷發(fā)展和深入，中國數(shù)字化森林的建設(shè)進(jìn)程不斷向前推進(jìn)［1］。地面三維激光掃描作為一種新型的可以自動、連續(xù)、快速地獲取目標(biāo)物的三維點云數(shù)據(jù)的測繪技術(shù)，已在林業(yè)資源調(diào)查［2］、 林分結(jié)構(gòu)研究［3］以及單木三維重建［4］等方面得到了初步應(yīng)用［5］。與傳統(tǒng)森林資源調(diào)查方法作業(yè)周期長、工作效率低、勞動強(qiáng)度大、數(shù)據(jù)單一相比，應(yīng)用地面三維激光掃描技術(shù)獲取測樹因子速度更快、精度更高，而且不會對林木造成損害，有利于保護(hù)生態(tài)環(huán)境［6］。利用三維激光掃描的點云數(shù)據(jù)繪制樹木枝干的等值線圖，不僅可以對點云數(shù)據(jù)進(jìn)行精簡和壓縮，得到的樹木枝干等值線模型還可用于單木的三維重建和樹木的形態(tài)結(jié)構(gòu)研究，對于樹木的可視化模擬具有非常重要的意義。傳統(tǒng)的離散矢量點云數(shù)據(jù)的等值線提取方法可以分為規(guī)則數(shù)據(jù)場的等值線提取和不規(guī)則數(shù)據(jù)場的等值線提取［7－11］。對于樹木枝干點云，不適合采用這些方法，因為對樹木枝干點云數(shù)據(jù)進(jìn)行表面建模，構(gòu)建真三維的Delauney三角網(wǎng)的技術(shù)還不成熟。當(dāng)點云數(shù)量過多時，構(gòu)建的樹木三角網(wǎng)格模型數(shù)據(jù)太大；而當(dāng)點云數(shù)量過少時，構(gòu)建的模型在細(xì)節(jié)上又不能令人滿意。本研究首先根據(jù)樹木的分枝特征將點云數(shù)據(jù)分割成不同的部分，并分析了點云的空間分布情況，然后將點云數(shù)據(jù)按照y坐標(biāo)分層，對于樹木各部分不同層次的點云利用迭代的凸包算法分別提取等值線。本方法解決的關(guān)鍵問題在于，針對由激光掃描得到的點云數(shù)據(jù)多、無拓?fù)湫畔?、低采樣間隔、高采樣細(xì)節(jié)的特點，在未構(gòu)建樹木三維模型的前提下，直接提取樹木枝干的等值線模型。

1 數(shù)據(jù)獲取和分析

采用FARO Laser Scanner Photon地面三維激光掃描儀，在湖南省株洲市黃豐橋林場進(jìn)行數(shù)據(jù)采集實驗。掃描樣木為落葉后的鵝掌楸Liriodendron chinensis。根據(jù)樣木周邊環(huán)境和樣木自身特點，選取3個角度，設(shè)立3站對樣木進(jìn)行掃描。由于要建立樣木精密的三維模型，每站首先采用了較低密度的采樣精度掃描全局，然后使用較高密度的采樣精度對樣木所在的區(qū)域進(jìn)行局部掃描。利用掃描軟件配準(zhǔn)不同站點的掃描數(shù)據(jù)。手動刪除樹木周圍物體的掃描點，去除噪聲點，分離出樹木的點云數(shù)據(jù)。

在對每一層點云進(jìn)行凸包計算時，需要對樹木不同部位的枝干分別計算，因此，需要把樹木枝干的點云數(shù)據(jù)分割成不同的部分。對此，在利用掃描軟件導(dǎo)出點云數(shù)據(jù)時，首先通過坐標(biāo)變換將點云數(shù)據(jù)的坐標(biāo)原點設(shè)置為樹木的根部，且樹木的主干大致平行于y軸，然后依據(jù)樹木的分枝特征，將樹木不同部位的點云數(shù)據(jù)分別導(dǎo)出，編號后合并。

2 點云數(shù)據(jù)的等值線提取

2.1 提取等值線折線

采用凸包算法提取每層點云數(shù)據(jù)的等值線折線。每層數(shù)據(jù)中，假設(shè)同一編號枝干的點集P={P0，…，Pn-1}，n為點的數(shù)量，其凸包是一個最小的凸多邊形Q，并且滿足P中的所有點或者在Q上，或者在Q的內(nèi)部。使用不同編號枝干上的點云可以構(gòu)建不同的凸包，凸包的數(shù)量等于點云層中枝條的數(shù)量。凸包算法在計算機(jī)圖形學(xué)、模式識別、圖像處理等領(lǐng)域有著廣泛的應(yīng)用。常用的構(gòu)建凸包的算法是Graham掃描法和Jarvis步進(jìn)法。

當(dāng)點集P有3個或3個以上的點時，凸包算法的過程如下：①計算點集最左邊的點為凸包的頂點的起點，如圖中的P0。②遍歷P中其他所有點，計算有向向量P0Pi。③如果P中的其他點全部在向量P0Pi的同一側(cè)，則Pi為凸包的下一個頂點，如圖1。

此過程執(zhí)行后，按照任意2點的順序把點按極角自動順時針或逆時針排序，而左側(cè)或右側(cè)的判斷則可以用矢量點積性質(zhì)來實現(xiàn)。

2.2 內(nèi)部離散點的歸類

利用凸包算法得到的凸多邊形邊界上的點是點集P的一個子集，可以認(rèn)為P代表任一枝條在某一高度范圍內(nèi)的等值線的一個近似。在此高度范圍內(nèi)，樹枝的表面仍有離散點沒有包括在該折線上，利用該折線代表其等值線過于簡單。因此，為了更加精確地表達(dá)等值線，需要進(jìn)一步將這些離散點歸類到當(dāng)前的折線上。

如圖2所示：QiQi+1是當(dāng)前折線上連接凸包中相鄰2點的一條線段，Pj為凸包內(nèi)部尚未加入到折線上的一個離散點。對于點Pj，應(yīng)該歸類于其最近的折線線段，并且該點在線段的投影Pj0應(yīng)該在線段上。連接 QiPj和 Qi+1Pj這2條線段，組成三角形QiQi+1Pj。

三角形QiQi+1Pj的2個內(nèi)角α和β，依據(jù)平面幾何的知識，可以得到以下的判別條件：假如α和β中有一個角度大于90°，那么投影點Pj0在線段QiQi+1外。假如α和β中有一個角度等于90°，那么投影點Pj0位于Qi或 Qi+1點上。假如α和β全部小于90°，那么投影點Pj0在QiQi+1內(nèi)部。

圖1 凸包算法Figure 1 Convex hull algorithm

圖2 離散點的歸類Figure 2 Classify discrete points

對每個尚未歸類的點Pj，計算該點到當(dāng)前等值線各線段的投影點，利用上述方法判斷投影點是否在其對應(yīng)的線段上，然后利用該點到各線段的距離，對Pj歸類。每個待定點只能歸類于α和β均小于90°并且距離最近的線段。

采用上述方法對當(dāng)前的等值線和待定點進(jìn)行歸類后，對于凸包折線的每條線段，重新計算其和歸類到該線段的待定點的凸包，由此得到局部區(qū)域的凸包。然后將所有的局部凸包在原凸包折線的線段處斷開，得到不閉合的折線段。按順序?qū)⑶昂笳劬€段進(jìn)行連接，即得到下一級折線。該過程如圖3所示。

圖3 迭代的凸包算法Figure 3 Iterative convex hull algorithm

整個算法的流程如圖4所示。

運用此方法可以保證每次迭代完成后的相鄰層之間的等值線不會產(chǎn)生邊緣交叉的問題，可以滿足等值線的基本特點。

第一次調(diào)用凸包算法得到的是一個凸多邊形，對此凸多邊形繼續(xù)使用凸包算法得到的局部區(qū)域也是凸多邊形，但是合并一系列的局部凸多邊形之后得到的等值線卻是凹多邊形。對凹多邊形再次進(jìn)行剩余點的歸類和調(diào)用凸包算法，得到的仍舊是凹多邊形。除第一次調(diào)用凸包算法的到的是凸多邊形，以后每次迭代后的合并結(jié)果都是凹多邊形。調(diào)用凸包算法完成后，剩下的每個待定點和當(dāng)前折線的拓?fù)潢P(guān)系是一致的，在折線的內(nèi)部或者外部。可以知道，奇數(shù)次迭代之后的待定點都在折線的外部；偶數(shù)次迭代之后的待定點都在折線的內(nèi)部。

3 研究實例

3.1 點云數(shù)據(jù)的空間分布

選取1株3年生的鵝掌楸，樹高為 2.26 m。設(shè)全部點云為［Pi（xi，yi，zi），i＝0，…，n－1，n為點的總數(shù)］。坐標(biāo)范圍為（xmin，ymin，zmin）～（xmax，ymax，zmax）。按照坐標(biāo) y方向，且按照一定的高度范圍，對數(shù)據(jù)進(jìn)行分層，得出任意一點 Pi（xi，yi，zi）的層序號為：

圖4 算法流程圖Figure 4 Algorithm flowchart

式（1）中：dz表示層的厚度，ki表示Pi點的層號。

分析合并后樣木點云數(shù)據(jù)，其坐標(biāo)范圍結(jié)果如下：xmin=－0.364 9，xmax=0.322 7，ymin=0，ymax=2.259 4，zmin=－0.266 5，zmax=0.462 8。該批點云的總共點數(shù)為104 808個，根據(jù)掃描精度，取dz=2.5 mm，則計算可知總共分為904層，平均每層點云數(shù)量為116個，所有層的點云數(shù)量分布以及和枝條數(shù)量的關(guān)系如圖5A所示。取y坐標(biāo)為1.000 m（即取1.008 5≤y＜1.001 0）的層數(shù)據(jù)，該層總共點數(shù)為117個，其空間分布如圖5B所示。

圖5 點云的空間分布Figure 5 Spatial distribution of points cloud

由圖5可以看出：使用地面激光掃描儀掃描的點云數(shù)據(jù)，沿樹高方向每層的數(shù)據(jù)量都很大。每層點云的數(shù)量隨枝條個數(shù)的增加有增大的趨勢。并且從點的分布圖中基本可以看出該層點云的等值線形狀。

3.2 等值線模型的提取

使用集成開發(fā)工具Visual Studio.Net，運用C#語言結(jié)合多媒體編程接口DirectX，編程實現(xiàn)樹木枝干等值線的提取。

點云數(shù)據(jù)分層后，每層數(shù)據(jù)中都至少包括一條枝干。當(dāng)點云層的厚度較小時，可以認(rèn)為此層中的點云具有同一高度值。圖6A為提取的整棵樹木的等值線圖。圖6B為當(dāng)高度等于0.80 m（0.798 5≤y＜0.801 0）時相應(yīng)的點云層，此層數(shù)據(jù)中共包括5個樹木枝干，編號分別為①③④⑤⑥。對不同的枝干調(diào)用上述迭代的凸包算法，分別提取出其凸包折線，圖6C為提取出的樹木枝條編號為③的凸包折線。由圖6可以看出在沒有先驗等值線知識和建立點云對象模型的條件下，利用迭代的凸包算法可以有效的對離散點數(shù)據(jù)進(jìn)行連接，形成點云的等值線模型。

圖6 等值線模型實例Figure 6 Contour model instance

3.3 等值線精度驗證

利用提取出的凸包折線的長度計算不同高度樹木枝干的直徑，與實際測量樹木枝干的直徑相對比驗證等值線的精度。其中，因樹干方向大致和凸包折線所在的平面垂直，凸包折線的長度可以當(dāng)作樹干的圓周長，以此計算樹干直徑。但是，枝條的方向和凸包折線所在的平面有一定的夾角，矯正后才可用來計算枝條的直徑。矯正過程如圖7所示。

圖7 直徑矯正圖Figure 7 Diameter correction

圖7中：P為枝條最低端的凸包折線，Q為待矯正凸包折線，以凸包折線中距離最遠(yuǎn)的2點的中點為凸包折線的中心，連接PQ的中心與樹干方向的夾角為θ，d0為矯正前枝條的半徑，d1為矯正后枝條的半徑，由式（2）可求出 d1。

利用上述方法可以提取出樹木不同高度枝干的直徑，部分?jǐn)?shù)據(jù)和實際測量數(shù)據(jù)對比如表1所示。其中，通過等值線模型計算出的數(shù)據(jù)記為掃描數(shù)據(jù)，編號1～15為不同高度樹干直徑，編號16～30為矯正后的樣木不同部位的枝條直徑。

由表1可以看出：掃描數(shù)據(jù)比實測數(shù)據(jù)整體偏大，主要原因是絕大多數(shù)凸包折線的長度大于樹木枝干實際的圓周長，因此，計算得出的直徑會偏大。以掃描數(shù)據(jù)為自變量x，實際測量的數(shù)據(jù)為因變量y，做出樣木枝干直徑的散點圖，并進(jìn)行線性回歸擬合，得到回歸方程。由圖8中的回歸方程可以看出：R2達(dá)到0.996 4，置信度為95%，說明計算數(shù)據(jù)和實測數(shù)據(jù)緊密相關(guān)。

圖8 直徑散點圖Figure 8 Scatter diagram of diameter

表1 掃描數(shù)據(jù)和測量數(shù)據(jù)對比Table 1 Comparison of scan data and actual data

利用掃描數(shù)據(jù)，通過回歸方程計算直徑的理論值，根據(jù)公式（理論值－實測值）/實測值，得到誤差為3.4%，誤差數(shù)值滿足林業(yè)測樹的要求，說明利用此種方法提取的等值線模型可以用于樹木枝干任何高度處的直徑的測量。

4 結(jié)語

針對利用地面三維激光掃描技術(shù)得到的樹木枝干點云數(shù)據(jù)建模困難的問題，提出了采用分層提取等值線的方法建立樹木枝干的等值線模型，所建立的模型包含樹木枝干的原始信息，可以用于樹木基本測樹因子的測量。

使用凸包算法提取每層點云中樹木不同枝干的等值線，然后將剩余的離散點歸類到當(dāng)前的等值線線段，迭代掉用凸包算法，得到該層的等值線模型。此方法不僅算法簡單，而且提取的等值線不存在邊緣交叉問題，提取的等值線的精度可以滿足林業(yè)測樹的要求。

應(yīng)該指出的是，使用此方法提取樹木枝干的等值線對原始數(shù)據(jù)的要求很高，點云采樣密度不能過大，而且必須嚴(yán)格去噪。如何在此等值線模型的基礎(chǔ)上構(gòu)建樹木枝干的Delauney網(wǎng)格結(jié)構(gòu)將是今后進(jìn)一步研究的重要內(nèi)容。

［1］趙陽，余新曉，信忠保，等.地面三維激光掃描技術(shù)在林業(yè)中的應(yīng)用與展望［J］.世界林業(yè)研究，2010，23（4）：41－45.ZHAO Yang，YU Xinxiao，XIN Zhongbao，et al.Application and outlook of terrestrial 3D laser scanning technology in forestry［J］.World For Res，2010，23（4）∶41－45.

［2］JAKOB W.Application and statistical analysis of terrestrial laser scanning and forest growth simulations to determine selected characteristics of Douglas-fir stands［J］.Folia For Polon Ser A，2009，51（2）∶123－137.

［3］GABOR B，GEZA K.Algorithms for stem mapping by means of terrestrial laser scanning［J］.Acta Silv Lign Hung，2009，5（1）∶119－130.

［4］TANSEY K,SELMES N.Estimating tree and stand variables in a corsican pine woodland from terrestrial laser scanner data［J］.Int J Rem Sens，2009，30（19）∶5195－5209.

［5］宋宏.地面三維激光掃描測量技術(shù)及其應(yīng)用分析［J］.測繪技術(shù)裝備，2008，10（2）：40－43.SONG Hong.Analysis and utilization of terrain 3D laser scanning survey technique［J］.Geom Technol Equ，2008，10（2）∶40－43.

［6］王劍，周國民.基于激光掃描儀的樹干三維重建方法研究［J］.微計算機(jī)信息，2009，25（8－3）：228－230 WANG Jian，ZHOU Guomin.3D Reconstruction of tree stems based on laser scanning of standing［J］.Microcomput Info，2009，25（8－3）：228－230.

［7］蔣瑜，杜斌，盧軍，等.基于Delaunay三角網(wǎng)的等值線繪制算法［J］.計算機(jī)應(yīng)用研究，2010，27（1）：101－103.JIANG Yu，DU Bin，LU Jun，et al.Algorithm of drawing isoline based on Delaunay triangle net［J］.Appl Res Comput，2010，27（1）∶101－103.

［8］LORENSEN W，CLINE H.Marching cubes∶A high resolution 3D surface construction algorithm［J］.Computer Graph，1987，21（4）∶163－169.

［9］吳杭彬，劉春.激光掃描數(shù)據(jù)的等值線分層提取和多細(xì)節(jié)表達(dá)［J］.同濟(jì)大學(xué)學(xué)報，2009，37（2）：267－271.WU Hangbin，LIU Chun.Point cloud-based stratified contour extraction and its multi-lod expression with ground laser range scanning［J］.J Tongji Univ，2009，37（2）∶267－271.

［10］趙夫來，郝向陽.根據(jù)地性線繪制等高線的研究［J］.測繪學(xué)報，1999，28（4）：345－349.ZHAO Fulai，HAO Xiangyang.A research on contour plotting based on bone lines［J］.Acta Geod Cartograph Sin，1999，28（4）∶345－349.

［11］龔有亮，何玉華，付予傲，等.一種實用的等高線內(nèi)插算法［J］.測繪學(xué)院學(xué)報，2002（3）：36－39.GONG Youliang，HE Yuhua，F(xiàn)U Yu’ao，et al.A practical contour interpolation algorithm［J］.J Inst Survey Map，2002（3）∶36－39.