基于三維激光點云的樹冠體積計算方法研究

李宏星，劉玉春

(1.湖北生態工程職業技術學院園林與建工學院，湖北 武漢 430200； 2.武漢市測繪研究院，湖北 武漢 430010)

1 引 言

三維激光掃描技術作業速度快，測量精度高，能夠獲取豐富的樹木點云數據。高精度的樹木點云數據信息量大，對林業、生態環境等研究具有重要的衍生價值，在城市森林調查、園林綠化、樹木管理等工作中得到了許多成功應用，例如，綠化數據的采集[1，2]、城市綠地建庫[3]、園林綠化工程竣工測量[4]、城市樹木特征提取[5]、三維綠量估計[6]等。

樹冠體積是評價樹木生長狀況、生態效益的一個重要因子。樹冠是不規則的非實心體，且內部充滿空隙，導致樹冠體積難以準確定義和計算。常用做法是將包裹整個樹冠的最小包絡體體積作為樹冠的體積。傳統方法將樹冠近似看作規則幾何體來估算體積，精度較差[7]。隨著三維激光掃描技術的引入，基于三維點云計算樹冠體積成為一種效率高、結果可靠的方法，得到了廣泛應用[8～11]。為此，國內外學者提出了許多種樹冠體積計算方法。除了前述的簡單幾何體近似法，還有體元法和計算幾何法。體元法將三維空間劃分成一個個的小立方體，統計出內部含有數據點的小立方體數量，再乘上小立方體的體積得到樹冠體積[8]。這種方法的優點是算法簡單，可以排除樹冠內部的空洞，缺點是計算量較大，無法顧及因遮擋產生的數據空缺現象。此外，計算的體積與體元尺度相關，體元越大，得到的樹冠體積越大。計算計算幾何法直接根據三維點云構建不規則形體，并將其體積作為樹冠體積，例如凸包法[12]、α-shape法[11]。為了提高這類方法的計算精度，通常將樹冠進行分層，每層看作是一個臺體，各層體積總和即為樹冠體積。分層方法的優點是可較好地避免體積高估，但是層數的劃分對估計結果有較大影響。分層法的一個關鍵是計算每層切面的面積，常用的凸包多邊形算法仍會導致高估[13，14]。為克服凸包算法估計結果偏大的問題，筆者提出了一種基于扇形分割的算法[15，16]。由于樹冠體積真值未知，無法評價各種體積算法的絕對精度。實際應用中只能根據主觀因素(如研究目的和期望的準確性)和客觀因素(如樹冠特征、時間和成本因素)來選擇適當的樹冠體積計算方法。這需要對各種算法的特性有清楚的認識，為此，本文分析幾種常用的樹冠體積計算方法的優缺點，并通過實例對各種方法的計算結果進行比較。

2 樹冠體積計算方法

2.1 簡單幾何體法

簡單幾何體法將樹冠看作一個規則的三維幾何實體，如圖1(b)，其體積可以幾何體的體積公式計算。常用的幾何體有圓錐體、半球體、橢球體、圓柱體和拋物體等，它們的體積公式如下：

(1)

式中，V為樹冠體積，D和H分別為冠徑和冠高。

這種方法只需要測量冠幅和冠高作為輸入，計算簡單，是傳統測樹學最常用的樹冠體積估計方法。但是，由于樹冠形狀不規則，近似幾何體的選擇取決于樹種、樹木的生物特征，以及研究人員的主觀判斷。不同樹種的樹冠形態各異，沒有哪一種幾何體適合逼近所有樹種的樹冠。因此，簡單幾何體法通常比較適用于計算形狀比較規則的樹冠體積。對于形狀復雜的樹冠，該方法的估計結果一般比實際體積嚴重偏大。

2.2 體元法

體元法將樹冠所在的空間范圍劃分成邊長為d的小立方體，遍歷這些小立方體，判斷每個立方體是否包含有點云數據，有則將其視為有效體元并保留，否則視為無效體元進行刪除，如圖1(c)所示。把最后剩下的有效體元體積之和作為樹冠體積，即：

V=N·d3

(2)

其中N是有效體元的數量。

此方法的優點在于適用于任何不規則的樹冠，算法簡單易于實現，僅需要判斷體元的有效性，有效體元為樹冠體積部分，無效體元為樹冠之外或樹冠內部空隙部分。通過刪除無效體元可以剔除樹冠外圍的凹陷部分和內部空隙，避免過大估計樹冠體積。但是，體元法的缺點也同樣明顯。首先是通過點云數據無法區分內部空隙是否是真的不存在植被元素，還是因為被遮蔽導致沒有激光反射點。例如枝干內部空間就不可能形成反射點，如果體元設置太小，一些大的枝干內部空間將不被計入樹冠體積。另一個問題是算法雖然簡單但是比較耗時，需要遍歷數量巨大的體元，而且對每個體元都要遍歷一次點云數據以判斷是否為有效體元。更重要的是體元法計算的體積與采用的體元大小高度相關。大的體元難以刻畫樹冠結構細節，但可以補償遮蔽區。小體元能精細描繪樹冠結構，但是無法補償遮蔽區體積。如何確定合適的體元大小是關鍵，也是一個難點。韋雪花等[8]利用圓柏進行實驗指出體元邊長的最佳取值為D/10，因其采用的樣本平均冠徑為 2.39 m，故采用的體元邊長為 0.2 m。

2.3 凸包法

凸包法可分為三維凸包法和二維凸包法。三維凸包法構建包含整個樹冠三維點集的最小凸包體，以這個三維凸包體的體積作為樹冠體積，如圖1(d)所示。也可以首先將樹冠進行分層，對每層的點云數據構建三維凸包，并計算它們的體積之和。二維凸包法是先將樹冠等間隔分層切片，如圖1(e)采用二維凸包算法對每層切片構建最小凸多邊形并計算其面積，樹冠體積采用下列公式計算[9，14]：

圖1 樹冠體積算法示意圖

(3)

式中，h是分層間距，n是切片層數，Si為第i層切片面積。

常用的凸包算法有增量算法、包裹法、Graham掃描法等。這些傳統凸包算法的不足是會將樹冠邊緣的凹陷空間計入樹冠體積，導致體積估值偏大，三維凸包算法尤為嚴重。為了改進二維凸包法，董亞涵等提出了迭代漸進的凸包算法[14]，但改進算法提取得到的外輪廓實際不是凸多邊形。對于分層算法，不同的分層間距值將對計算結果產生影響。理論上，這些間距取值越小，體積計算越準確。但是，由于樹冠是離散采樣，劃分過細不僅降低計算效率，還可能得出不合理的結果。林松探討了不同分層間距對樹冠體積計算結果的影響，不同形狀的樹冠影響方式不同[9]。

2.4 扇形分割法

如圖1(f)，扇形分割法基于分層切片的思想，先將樹冠分層，再把每層切面分為若干小扇區。取切面中心為極點，扇區內最大極距為扇區半徑，對應的數據點作為輪廓點得到切面的外輪廓，切面面積就是這些扇形面積的總和，可表示為[15]：

(4)

式中，m是每層劃分的扇區數量，△θ是扇區的圓心角，Si，j為第i層切面第j個扇區的面積，ri，j是對應的最大極距。將式(3)代入式(4)可得體積計算式為：

(5)

除了分層間距h外，扇形分割法中還有一個關鍵參數是扇區圓心角△θ。當△θ=2π時切面面積最大，為每層數據點的外接圓面積。隨著△θ取值減小，扇區劃分越多，得到的外輪廓趨近于樹冠外緣，得到的樹冠體積越精確。然而，與分層間距類似，并非扇區劃分得越小就越好，需要根據點云數據具體特征選擇合適的圓心角大小。

3 實驗結果與分析

3.1 實驗數據

實驗樣本采集自圭亞那和印度尼西亞的熱帶森林，由荷蘭瓦赫寧根大學土地利用和碳排放數據網站(http：//lucid.wur.nl/datasets/terrestrial-lidar-of-tropical-forests)提供[17]。樹木點云數據采集使用的儀器是RIEGL VZ-400 3D地面激光掃描儀，該儀器使用的激光波長為 1 550 nm，角分辨率為0.06°，可實現360°全方位掃描，最大天頂角可達100°。表1列出了樣樹的基本信息，GUY表示圭亞那，IND表示印度尼西亞，共20棵樹木，包括11個樹種。為了獲得高質量的點云數據，每棵數設置了13個掃描測站，每站掃描儀沿垂直軸和水平軸旋轉進行了2次掃描。獲得的高密度點云如圖2所示，其中最小的點云也超過24萬個采樣點。由圖2可以看到這些樹木的形態差異較大，樹高最大值、最小值分別為 38.3 m、21.2 m，平均值為 29.75 m。利用卷尺實測的冠徑最大值為 21.82 m，最小值為 5.5 m，均值為 14.8 m。

表1 樣樹的基本信息

圖2 20棵樣樹的點云

3.2 結果分析

分別采用圓錐體法、半球體法、體元法、三維凸包法、二維凸包法和扇形分割法計算了20棵樹木的樹冠體積。由于點云數據量大，顧及計算精度和計算效率，體元法的體元邊長取為 0.4 m。對于二維凸包法和扇形分割法，為了確定合適的分層間距，以 0.01 m的步長改變分層間距，得到樹冠體積估值的相對變化量，取相對變化量小于給定閾值時對應的分層間距作為最佳值。圖3是20棵樹木樹冠體積相對變化量均值隨分層間距的變化曲線，取相對變化量的閾值為0.2%，對應的分層間距為 0.5 m。利用類似的方法我們確定了扇形分割法中△θ的最優取值為2°。各種方法估算的樹冠體積如圖4所示，表2給出了不同方法結果的線性回歸方程、決定系數(R2)和相關系數。

圖3 不同分層間距對應的樹冠體積相對變化量(步長0.01 m)

圖4表明在選用的6種方法中三維凸包法估算的樹冠體積嚴重偏大，對于大型樹冠尤為嚴重，例如第1、7～9、16和20號樹木。三維凸包算法將樹冠外圍凹陷空間計入樹冠體積，樹冠越大外圍的凹陷空間可能就越多，因此，三維凸包法對于大型樹冠或外圍存在較多凹陷的樹冠容易高估其體積。兩種簡單幾何體法(圓錐體和半球體)的估值也偏大，但是兩種方法的結果還算比較接近，相關系數達到0.907。第7～9號樹木的樹冠較為扁平，冠徑比冠高大得多，導致兩種幾何體方法的估值相差較大。對兩種凸包法，二維凸包法得到的結果明顯小于三維算法的結果，但它們之間的相關系數高達0.97，顯著高于其他方法與三維凸包法的相關性。體元法和扇形分割法估算的體積較為接近，在全部方法它們的結果相對偏小。

圖4 不同方法估算的樹冠體積

由表2可以看出，不同方法計算得到的樹冠體積在統計上具有較強相關性，兩兩之間的相關系數都在0.8以上。表2中的回歸模型也展示了不同方法結果之間良好的線性關系，線性擬合的R2值大多都超過0.8。其中，只有體元法和簡單幾何體法之間的R2值相對較小(<0.8)。原因主要是樣本樹冠形態較為不規則，所用的圓錐體或半球體不能很好地刻畫樹冠形狀，而其他幾種方法基本能反映樹冠的輪廓形狀。扇形分割法與二維凸包法之間的相關性最強，相關系數為0.989，這可能是由于兩種方法均基于分層切片的緣故。采用扇形分割法計算的樹冠體積要小于二維凸包法的計算結果，說明扇形分割法能夠有效克服傳統凸包算法將樹冠邊緣的凹陷空間計入樹冠體積的缺陷。

表2 不同方法樹冠體積計算結果的回歸方程、決定系數(R2)和相關系數

對于高密度的完整樹冠點云，采用較小體元的體元法樹冠體積估值可以近似看作真值[9]。從樣樹點云的采集方式和表1給出的點云數量可知樣樹的點云密度極高，因此，可以將體元法得到的樹冠體積為真值作為參考，來評價其他5種方法的效果。對比5種方法，扇形分割法與體元法的相關系數最大為0.958，說明扇形分割法結果與真值的相關性最強。從圖5給出的各方法體積估計值與體元法結果的線性回歸關系可以看出，扇形分割法與體元法的結果最為接近。而且它們的回歸模型R2值為0.918，也是5種方法中最大的，說明線性擬合度最優。由此可見，扇形分割法是5種方法中最為精確的樹冠體積計算方法。

圖5 以體元法為參考的線性回歸關系

4 結 論

結果分析表明，雖然不同樹冠體積計算方法得到的數值結果有較大差異，但是在統計上有顯著的相關性。簡單幾何體法和凸包法計算的樹冠體積相對偏大，而扇形分割法和體元法的結果相對偏小。若以體元法結果作為參考，扇形分割法的準確性最高。相比于體元法的時間復雜度，扇形分割法更具優勢，是一種能夠兼顧精度和計算效率的樹冠體積計算方法。