基于三維激光點(diǎn)云的樹冠體積計(jì)算方法研究

李宏星，劉玉春

(1.湖北生態(tài)工程職業(yè)技術(shù)學(xué)院園林與建工學(xué)院，湖北 武漢 430200； 2.武漢市測(cè)繪研究院，湖北 武漢 430010)

1 引 言

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

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

2 樹冠體積計(jì)算方法

2.1 簡(jiǎn)單幾何體法

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

(1)

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

這種方法只需要測(cè)量冠幅和冠高作為輸入，計(jì)算簡(jiǎn)單，是傳統(tǒng)測(cè)樹學(xué)最常用的樹冠體積估計(jì)方法。但是，由于樹冠形狀不規(guī)則，近似幾何體的選擇取決于樹種、樹木的生物特征，以及研究人員的主觀判斷。不同樹種的樹冠形態(tài)各異，沒(méi)有哪一種幾何體適合逼近所有樹種的樹冠。因此，簡(jiǎn)單幾何體法通常比較適用于計(jì)算形狀比較規(guī)則的樹冠體積。對(duì)于形狀復(fù)雜的樹冠，該方法的估計(jì)結(jié)果一般比實(shí)際體積嚴(yán)重偏大。

2.2 體元法

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

V=N·d3

(2)

其中N是有效體元的數(shù)量。

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

2.3 凸包法

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

圖1 樹冠體積算法示意圖

(3)

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

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

2.4 扇形分割法

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

(4)

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

(5)

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

3 實(shí)驗(yàn)結(jié)果與分析

3.1 實(shí)驗(yàn)數(shù)據(jù)

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

表1 樣樹的基本信息

圖2 20棵樣樹的點(diǎn)云

3.2 結(jié)果分析

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

圖3 不同分層間距對(duì)應(yīng)的樹冠體積相對(duì)變化量(步長(zhǎng)0.01 m)

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

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

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

表2 不同方法樹冠體積計(jì)算結(jié)果的回歸方程、決定系數(shù)(R2)和相關(guān)系數(shù)

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

圖5 以體元法為參考的線性回歸關(guān)系

4 結(jié) 論

結(jié)果分析表明，雖然不同樹冠體積計(jì)算方法得到的數(shù)值結(jié)果有較大差異，但是在統(tǒng)計(jì)上有顯著的相關(guān)性。簡(jiǎn)單幾何體法和凸包法計(jì)算的樹冠體積相對(duì)偏大，而扇形分割法和體元法的結(jié)果相對(duì)偏小。若以體元法結(jié)果作為參考，扇形分割法的準(zhǔn)確性最高。相比于體元法的時(shí)間復(fù)雜度，扇形分割法更具優(yōu)勢(shì)，是一種能夠兼顧精度和計(jì)算效率的樹冠體積計(jì)算方法。