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

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

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

1 數據獲取和分析

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

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

2 點云數據的等值線提取

2.1 提取等值線折線

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

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

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

2.2 內部離散點的歸類

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

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

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

圖1 凸包算法Figure 1 Convex hull algorithm

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

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

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

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

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

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

第一次調用凸包算法得到的是一個凸多邊形，對此凸多邊形繼續(xù)使用凸包算法得到的局部區(qū)域也是凸多邊形，但是合并一系列的局部凸多邊形之后得到的等值線卻是凹多邊形。對凹多邊形再次進行剩余點的歸類和調用凸包算法，得到的仍舊是凹多邊形。除第一次調用凸包算法的到的是凸多邊形，以后每次迭代后的合并結果都是凹多邊形。調用凸包算法完成后，剩下的每個待定點和當前折線的拓撲關系是一致的，在折線的內部或者外部?？梢灾?，奇數次迭代之后的待定點都在折線的外部；偶數次迭代之后的待定點都在折線的內部。

3 研究實例

3.1 點云數據的空間分布

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

圖4 算法流程圖Figure 4 Algorithm flowchart

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

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

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

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

3.2 等值線模型的提取

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

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

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

3.3 等值線精度驗證

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

圖7 直徑矯正圖Figure 7 Diameter correction

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

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

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

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

表1 掃描數據和測量數據對比Table 1 Comparison of scan data and actual data

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

4 結語

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

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

應該指出的是，使用此方法提取樹木枝干的等值線對原始數據的要求很高，點云采樣密度不能過大，而且必須嚴格去噪。如何在此等值線模型的基礎上構建樹木枝干的Delauney網格結構將是今后進一步研究的重要內容。

［1］趙陽，余新曉，信忠保，等.地面三維激光掃描技術在林業(yè)中的應用與展望［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］宋宏.地面三維激光掃描測量技術及其應用分析［J］.測繪技術裝備，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］.微計算機信息，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三角網的等值線繪制算法［J］.計算機應用研究，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］吳杭彬，劉春.激光掃描數據的等值線分層提取和多細節(jié)表達［J］.同濟大學學報，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］.測繪學報，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］龔有亮，何玉華，付予傲，等.一種實用的等高線內插算法［J］.測繪學院學報，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.