基于改進α-shape算法的三維點云樹冠體積計算方法

程 鋼 王敬宇 楊 杰 趙宗澤 王 磊

(1.河南理工大學測繪與國土信息工程學院，焦作 454000；2.河南理工大學教務處，焦作 454000)

0 引言

森林是重要的生態與經濟資源，包含了地球上約80%的植被生物量，對全球碳循環平衡具有非常重要的意義[1]。隨著三維綠量[2]、冠積指數[3]等植被生物量統計參數概念的提出，三維樹木因子的測定也變得愈發重要。樹冠作為樹木的重要組成部分，其體積參數的確定對森林生物量和三維綠量的統計分析具有重要作用。

傳統的樹冠體積量測方法[2-5]通過構建樹冠體積預估方程進行量測，這種方法要考慮樹種等因素的影響，其作業周期長、擬合形狀差，測量精度很難保證，易造成計算結果不穩定、不準確等方面的問題。近年來，三維激光掃描技術已逐漸被應用到逆向工程[6]、文物保護[7]以及傳統測繪[8]等領域，該技術可以在不接觸物體的情況下對物體進行遠距離測量，這為無損、快速、準確地量測樹冠體積提供了契機。基于激光點云數據進行樹冠體積計算的方法主要分為基于內部結構的樹冠體積計算和基于外部輪廓的樹冠體積計算兩類。第1類計算方法，如文獻[9-10]提出以固定的體元來模擬不規則樹冠形狀的體元模擬法，該方法對于濃密程度較高的樹木由于遮擋原因會導致其利用三維點云數據構建的結構內部產生“空洞”現象，使測算結果相較于真值偏小。第2類計算方法，如文獻[11]通過改進Delaunay算法機制，采用一種基于空間分割分塊優先機制的三角網表面重建算法來提取樹冠三維信息，具有較好的穩定性;文獻[12]通過提取樹冠冠體邊緣點云，應用不規則三角網TIN計算冠體體積，由于每次提取的輪廓點僅有4個，故樹冠體積計算結果不夠準確;文獻[13]將樹冠分割為多個不規則的臺體，但是將每個點云切片的橫截面看作圓形或橢圓形，導致計算結果偏差較大;文獻[14]對文獻[13]中的臺體法進行改進，通過計算每個樹冠點云切片的凸包面積來求取臺體體積，但計算結果仍有較大偏差;文獻[15]運用三維凸包算法計算樹冠凸包體積，由于未考慮樹冠邊緣間的空隙，故計算結果較真值偏差仍較大。

總的來看，體元法基于樹冠內部結構計算樹冠體積，有利于消除樹冠內部無效體積，對于內部結構遮擋較少的稀疏枝葉樹冠體積計算有一定優勢，但對于枝葉較為濃密的樹冠，由于缺少樹冠內部點云數據，其計算結果往往不準確。因此，由于遮擋的緣故，體元法通常存在樹冠體積低估的問題。臺體法將三維樹冠體積計算轉換為多個二維切片面積計算，使其可以在充分利用樹冠外部幾何特征的基礎上來求解三維樹冠信息，是一種簡單、快速、有效的樹冠體積量測方法，其準確度除與樹冠內部結構有關外，最關鍵的步驟在于臺體外輪廓的精確計算。為了提高臺體法的準確度，文獻[16]提出了一種基于過濾三角網的樹冠邊界精確提取算法，對最優分層間距與過濾閾值進行了討論，在考慮樹冠外部空隙的情況下得到更為準確的樹冠體積；文獻[17]將邊長閾值作為圓直徑，在切片凸包的基礎上搜尋“疑似邊界點集”中與直徑端點組成最大角的點作為新邊界點，但未對邊長閾值做進一步討論。上述2種方法均通過設置固定邊長閾值來消除樹冠邊緣存在的空隙，閾值過小，容易使整個切片碎片化，閾值過大，則會高估切片面積。考慮到動態閾值對于點云外輪廓形狀描述的重要性，本文在α-shape算法[18-22]的基礎上，借鑒三角網生長算法的思想[23-26]，提出一種考慮邊界點云密度的動態閾值改進α-shape算法，以期提取更為準確的點云切片邊界，并利用臺體法的基本思想獲取更為準確的樹冠體積。

1 研究方法

1.1 臺體法計算樹冠體積

臺體法計算樹冠體積是把單棵樹木整個樹冠從下到上看作是由許多個臺體和一個錐體所組成的幾何體，如圖1所示，通過對樹冠點云數據基于高程進行切片并利用離散化格林公式求解每個切片面積S(圖2)，進而計算出每個臺體體積，然后將這些臺體和錐體體積進行累加計算，即可算出整個樹冠的體積。計算公式為

(1)

(2)

式中m——端點數量

(xk，yk)、(xk+1，yk+1)——第k、k+1個端點的平面坐標

Si、Si+1——第i、i+1層點云切片面積

n——切片數量

hi——相鄰切片之間的間隔

hn+1——圓錐體高

V——樹冠體積

計算過程中，點云切片面積的準確性至關重要，直接影響到后續樹冠體積的計算，合理、準確地求取樹冠點云切片面積成為臺體法計算樹冠體積的關鍵問題。

1.2 改進α-shape算法

α-shape算法作為經典的邊緣探測算法之一，最早是由EDELSBRUNNER等[18]提出并被用于獲取散亂點集的輪廓，本質上使用該算法所獲取的圖形為歐幾里得空間中點集三角剖分的子集。文獻[19]通過引入加權思想來控制點集中每一個點對周圍點的影響力進而實現捕捉不同層次的局部空間特征，使得加權α-shape算法被廣泛應用到局部特征提取領域，文獻[20-21]均采用這種加權算法來實現對局部特征圖形結構的探測。文獻[27]通過計算檢測點的k個鄰近點平均距離和增設調節因子，自適應設置滾動圓半徑α，但未考慮實測點云中激光掃描噪聲和散亂點云點間距離散程度對α的影響。在提取樹冠散亂點云輪廓方面，文獻[28]通過使用α-shape算法獲取了樹冠體積并對結果進行了驗證，但其獲取輪廓的α是固定的，理論上仍存在提升樹冠體積精度的空間。

二維空間中，α-shape算法可以想象成一個半徑為α的圓在一堆無序點集P中進行滾動，當α足夠大時，這個圓不會滾動至點集內部，此時的滾動痕跡即為該點集的邊界[27]，如圖3所示。本文將α-shape算法運用到點云切片外輪廓的獲取，根據不同切片的點云結構，動態設定α。

根據文獻[18]可知，當半徑α趨于無窮小時，點集P中的所有點都有可能是邊界點；當半徑α趨于無窮大時，會形成一個包含點集P中所有點的凸包。理論上，存在一個閾值αb，當α≥αb時，其所形成的邊界多邊形可以包含點集P中的所有點。

本文通過迭代方式來找尋每個點云切片的閾值αb并獲取該閾值下切片外輪廓多邊形的面積。迭代公式為

αi+1=αi+Δα(i=1,2,…,n)

(3)

式中αi、αi+1——第i、i+1次迭代時圓的半徑

Δα——迭代步長i——迭代次數

設樹冠點云位于空間直角坐標系OXYZ中，OZ軸為高程軸。對樹冠點云基于高程進行等間距切片處理，設其中一個點云切片的點集為D，則改進α-shape算法構建切片外輪廓多邊形及求取面積的流程圖如圖4所示，具體如下：

(1)輸入點云切片并設置初始α。

(2)構建鄰域點集。將點集D投影至XOY平面，獲取投影點集P中每一個點Pi的半徑為2α的圓形鄰域內的點集，其中i=1,2,…,n。為了提高其查詢速度，本文借助了kd-tree結構來獲取點集P中每一個點的鄰域點集。

(3)設置邊界初始點。由凸包性質可知，點集中任意一個凸包極值點均為該點集的邊界點，為了便于獲取可以將初始點設置為OX或OY方向上其中一個極值點即可，本文以OY方向上的最小值點PyMin為例，將其設置為初始邊界點并加入到邊界點集L中。

(4)獲取整個切片的邊界點集L。包括以下幾個步驟：

①尋找下一個邊界點。以當前邊界點Li為基礎，從其鄰域點集Ni中選取任意一點n1，求出過點Li和n1且半徑為α的兩個圓o和o′；遍歷點集Ni，依次判斷鄰域中其他點是否落入圓o和o′內部，如果存在一個不包含其他鄰域點的圓(“空圓”準則，見圖5)，則可以判斷出點n1是點集P的邊界點；對Ni內下一點重復上述過程，直至對Ni內所有點均使用“空圓”準則判別之后，統計其邊界點的個數。

其中，初始點PyMin(L1)較為特殊，該點在尋找下一個邊界點的過程中不需要統計其鄰域內的邊界點個數，只需要在NyMin中找到一個符合“空圓”準則的邊界點即可執行步驟①，反之，如果一個邊界點也找不到，則跳轉至步驟(5)。

②判斷邊界點的唯一性。如果新找到的邊界點有且僅有1個，則執行步驟③。如果新找到的邊界點個數為0或多于1個，則跳轉至步驟(5)，如圖6a中，邊界點7的鄰域點集為空，新邊界點的個數為0，邊界無法回到起始點。圖6b中，邊界點3存在2個新邊界點，邊界出現2個可能的走向。圖6這2種情況均說明當前閾值無效，需要更新α。

③判斷唯一邊界點是否為初始點。如果該邊界點為初始點PyMin(圖7)，則執行步驟④；反之，則將該邊界點添加至點集L中，并執行步驟①。

④判斷點集L是否包含點集P中的極值點。獲取點集P中多個方向上的凸包極值點[29]，如圖8a所示，根據文獻[18]可知，相同點集下，該點集的凸包極值點亦是其邊界凹包的極值點，故這些極值點一定為點集P的邊界點。如果點集L未包含這些極值點，則該凹包僅僅為點集P的一個局部凹包，如圖8b所示，則跳轉至步驟(5)；反之，輸出相應α值下的凹包并執行步驟(6)，如圖8c所示。

(5)更新α并判斷α是否達到終止條件。利用式(3)迭代更新α，之后將α與迭代閾值αmax相比較，如果大于αmax則停止迭代過程，輸出一個包含點集P的凸包面積；反之，執行步驟(1)重新設置α。設置終止條件的原因在于：一方面，當α過大時，提取的邊界形狀趨于凸包，另一方面，可以避免程序陷入死循環。

(6)求解切片外輪廓多邊形面積。利用離散化格林公式求解出切片外輪廓多邊形面積的絕對值并輸出，則整個邊界提取過程結束。

最后，將整個樹冠的所有點云切片均重復步驟(1)～(6)以獲取其相應的切片外輪廓面積，求取對應的臺體體積，進而通過累加求出更為準確的樹冠體積。

2 分層間距與迭代步長的確定

2.1 樹冠點云獲取

本次實驗研究區域為某校園內，選擇晴朗天氣采集樹木點云數據，采集平臺為Leica公司MS60型全站儀，掃描分辨率設置為0.01 m×0.01 m(與儀器激光發射中心某個距離的平面上，點云點間距在水平和垂直方向上均為0.01 m)，每秒測量3 600個點。選取校內14棵樹進行掃描，獲取樹木點云數據后，對數據進行去噪和濾波處理，僅留下純粹的樹冠點云數據，如圖9所示。算法實驗在一臺CPU為Intel Core i5-9300H、頻率為2.4 GHz的便攜式計算機上進行，開發平臺為Matlab R2018a，α的初始值設置為0.01 m。

2.2 分層間距及迭代步長

使用改進α-shape算法計算樹冠體積時，不同的迭代步長與樹冠分層間距會直接影響樹冠體積的計算效率與準確度，并且考慮到方法的實用性等因素，需要為改進α-shape算法確定最優的迭代步長以及分層間距。由微積分原理可知，樹冠分層間距越小，使用該方法計算的樹冠體積越準確。因此，將迭代步長為0.01 m、分層間距為0.1 m時的計算結果作為體積參照。為了消除不同樹木之間計算時間與體積的差異，更為合理表征方法的計算效率與準確度，本文對不同樹木樣本的樹冠體積與計算時間進行歸一化處理，計算公式為

(4)

(5)

式中nt、nV——計算時間、樹冠體積歸一化值，nt越小表示計算效率越高，nV越大表示體積損失越大

t、Vα——使用改進α-shape算法計算單木樹冠體積所用時長與體積

t0、V0——迭代步長為0.01 m、分層間距為0.1 m時所用時長和體積

采用不同迭代步長和分層間距對樣本數據進行實驗和統計，結果如圖10所示。分析迭代步長對運算時間的影響及其造成的體積損失(圖10a、10b)可知，當迭代步長Δα約為0.1 m時，計算效率的提升效果開始趨于緩和。考慮到隨著迭代步長的增大，樹冠體積損失也在增加，因此迭代步長Δα=0.05 m是一個較為理想的選擇。在Δα=0.05 m的基礎上，分析分層間距對時間和體積的影響(圖10c、10d)可知，當分層間距Δh=0.2 m時，樹木樣本體積損失情況與Δh<0.2 m時的體積損失情況相當，并且其計算效率亦有顯著的提升，故選擇分層間距Δh=0.2 m。

對樣本樹木樹冠體積計算的時長和體積進行統計和分析，見表1。分析表1可知，對女貞樹2來說，采用選定的Δα與Δh時，體積V1比V0僅損失了1.360 7%，而時間僅為原來的11.851 3%，運算效率大大提升。對白楊樹3來說，采用選定的Δα與Δh時，體積V1比V0僅損失了2.512 6%，而時間為原來的29.080 2%，運算效率明顯提升。體積損失最高的是銀杏樹3，V1比V0損失了9.167 3%，時間為原來的21.845 9%，運算效率亦有明顯提升。可見，在V1與V0的相對誤差不超過9.167 3%的情況下，改進α-shape算法的用時僅為參照數據用時的11.460 6%～31.581 9%，顯著提高了方法的運算效率和實用性，可以滿足樹冠體積測算的實際需求。

表1 樹冠樣本計算結果

3 實驗結果與分析

3.1 樹冠體積計算

分別使用體元累加法、Graham算法[14]和改進α-shape算法對所有樹冠樣本進行體積計算。參數設置方面，將改進α-shape算法的初值、迭代步長、分層間距、最大迭代閾值αmax分別設置為0.01、0.05、0.2 m、2；Graham算法的分層間距亦設置為0.2 m；參照文獻[9]，本文對不同體元邊長與體積增幅G之間的關系進行比較分析，從圖11可見，當體元邊長為10 cm左右時，體積增幅已趨于穩定，故將體元累加法中的體元邊長設置為10 cm(數據獲取方法同3.2節)。

(6)

式中Bi、Bj+1——體元邊長為j、j+1時所計算的樹冠體積

Gj——體元邊長為j+1與j的樹冠體積比值

將體元累加法、Graham算法和改進α-shape算法所計算出的樹冠體積分別記為v0、v1和v2，結果如表2所示。

通過對不同樹種不同方法獲得的單木樹冠體積(表2)進行分析，可以得出以下規律:

(1)對于枝葉密度稀疏的樹冠(如銀杏樹1～3)，v1、v2相較于v0均偏差較大，存在高估現象，v2相比于v1更接近v0。分析原因在于，對于枝葉密度稀疏的樹冠其內部的空隙可以被體元法更好地識別和去除，體積明顯小于使用外輪廓進行體積估算的方法；由于Graham算法沒有顧及到樹冠邊緣的空隙問題(圖12)，本文算法相較于Graham算法所計算的結果v1會更為準確。

(2)枝葉茂密的情況下(如女貞樹1～3、白楊樹1～8)，3種方法得到的樹冠體積更為接近。分析表2可知，v2恰好處于v0與v1之間，表明Graham算法的高估問題仍然存在，改進α-shape算法雖然略高于體元累加法，但相較于Graham算法更接近于體元累加法。分析原因在于，對于枝葉密度較高的樹冠，體元法對于內部空隙的去除作用減弱。

表2 不同方法計算樹冠體積

(3)通過表2可知，計算不同樣本樹冠體積時，求得了不同的α，說明本文方法對不同樹種樹冠的適用性較好。

綜合分析來看，樹冠體積計算準確度與樹冠自身的枝葉結構關系緊密。本文算法可以適用于不同樹種的樹冠，無論是低密度還是高密度樹冠，其結果相對于Graham算法可以去除凸包邊界帶來的高估問題，相對于體元累計法則更有利于樹冠總體占用空間的計算。

3.2 改進α-shape算法穩定性測試

計算結果的穩定性是衡量一個算法好壞的重要指標，從樣本數據中任意選取一個樹冠樣本進行算法穩定性測試。對原始點云數據經去噪濾波處理后，進行均勻抽稀。抽稀過程以樹冠點云的立方體包圍盒為基礎，通過構建不同大小的正方體體元，保留體元內距離體元幾何中心的最近點，剔除體元內其他點。使用這種方式可以將原始點云分別抽稀為分辨率為1～10 cm的樹冠點云(分辨率k表示每個邊長為k的體元內僅保留不多于1個激光點)，用以測試不同算法計算樹冠體積的效果。對實驗結果進行歸一化處理，即

(7)

(8)

式中Kv、Ks——抽稀后樹冠樣本的體積相對誤差和抽稀后點云占比。Kv越大則體積損失越大；Ks越小，則表示抽稀后剩余點云數據占原始點云數據比例越少

Fk、F1——點云體元邊長為k、1 cm時樹冠體積

Ck、C1——點云體元邊長為k、1 cm時點數量

測試結果見表3。表3第1列，點云體元表示1個點所占的立方體空間，數值越大，點云分辨率越低。分析表3可知，在設置相同參數情況下，隨著點云分辨率降低，不同算法計算的樹冠體積均呈下降趨勢。其中，體元法受點云分辨率的影響波動較大，當點云體元邊長為10 cm時，損失了原始樹冠體積的51.016 6%；Graham算法與本文算法體積下降的幅度相較于去除點云的比例相對較小，在點云體元邊長為10 cm，且抽稀后點云數量僅為原始點云樣本的3.478 0%的情況下，損失了原始樹冠體積的13.746 1%和11.804 6%。這說明本文算法在相同的迭代步長與分層間距下，計算結果受點云分辨率影響相對較小，高密度點云使用的相關參數對于低密度數據亦具有一定適用性，在一定程度上說明了本文算法的穩定性。

表3 抽稀與計算結果

4 結論

(1)樹冠體積計算的準確性與樹冠內部枝葉結構和點云密度均有關系。

(2)本文算法對不同樹種樹冠體積計算具有較好的適用性，無論對于高密度還是低密度樹冠，采用改進α-shape算法的樹冠體積計算結果具有良好的穩定性，在過濾樣本為原始點云樣本3.478 0%的情況下，僅損失了原始樹冠體積的11.804 6%，而且相較于已有其他方法更為準確，既避免了Graham凸包算法的高估問題，相較于體元累加法也更利于樹冠總體占用空間的計算。