劉靜云

（上海工程技術(shù)大學(xué) 城市軌道交通學(xué)院，上海 201620）

0 引 言

在林業(yè)資源調(diào)查中，樹木的胸徑（Diameter at Breast Height，DBH）是一項關(guān)鍵基礎(chǔ)數(shù)據(jù)，其對樹木生長趨勢、生物量估算及森林的碳循環(huán)研究有著重要意義。樹木胸徑通常指樹高1.3 m處的直徑數(shù)值。傳統(tǒng)的胸徑量測方式常采用卡尺直接量測或使用皮尺量測圓周大小，并以此計算樹木胸徑數(shù)值。然而這種量測方式時間長，且容易產(chǎn)生由人工目視帶來的量測誤差。隨著測量技術(shù)的發(fā)展，全站儀、經(jīng)緯儀等量測設(shè)備不斷應(yīng)用于林業(yè)調(diào)查中。但是，由于林下環(huán)境復(fù)雜，相關(guān)森林參數(shù)難以被精確獲取，這些設(shè)備的應(yīng)用仍存在一定的局限性。

地面激光掃描（Terrestrial Laser Scanning，TLS）技術(shù)近些年發(fā)展迅速，具有掃描范圍大、采集速度快、數(shù)據(jù)質(zhì)量優(yōu)、自動化程度高等優(yōu)點，非常適合于獲取精細(xì)的林下結(jié)構(gòu)數(shù)據(jù)，現(xiàn)已成為森林調(diào)查、林業(yè)管理和生態(tài)系統(tǒng)研究的重要方法。

近年來，許多學(xué)者針對三維激光點云展開了樹木胸徑量測的相關(guān)研究。大多數(shù)方法的量測思路是：先對點云數(shù)據(jù)高程歸一化處理，再取胸徑高度數(shù)據(jù)將其投影到二維平面中，借助圓或橢圓擬合方法，進(jìn)一步量測出圓的直徑、橢圓的長短半軸數(shù)值作為胸徑值。常見擬合方法有Hough變換、最小二乘擬合、基于隨機(jī)采樣一致性（Random Sample Consensus，RANSAC）的圓擬合等，但這種二維的方式存在傾斜樹木投影時很難被檢測與擬合的缺陷。部分方法還采用圓柱擬合的方式完成樹干的檢測，得到圓柱模型的軸向與半徑參數(shù)，以圓柱的直徑作為樹木胸徑值。此種方式相較于二維的方法能夠避免投影所帶來的問題，魯棒地擬合出樹木模型，但其擬合的精度略有降低。此外，多數(shù)研究關(guān)注于樹木胸徑的量測精度，卻很少兼顧方法的計算效率，而實際應(yīng)用中不僅需要考慮算法結(jié)果的精度，其運行效率也應(yīng)當(dāng)被納入評價體系中。

針對上述問題，本文提出一種自動化的樹木胸徑檢測與量測方法。該方法先通過布料模擬濾波（Cloth Simulation Filter，CSF）算法，得到點云數(shù)據(jù)的地面模型，并對原始數(shù)據(jù)進(jìn)行高程歸一化；再通過點云切片與聚類自動獲取潛在的樹干點云簇，最后基于RANSAC的圓柱擬合方法，檢測點云中的樹干并計算樹干胸徑。

1 方法與原理

本文提出的方法基于地面激光點云實現(xiàn)對樹木胸徑的自動化量測，是一套端到端的方法，整個量測流程可概括為地面模型生成、高程歸一化、點云切片、點云聚類、圓柱擬合等5個步驟完成。

1.1 點云數(shù)據(jù)預(yù)處理

樹木胸徑參數(shù)量測計算時，與地面點云的數(shù)字高程模型（Digital Elevation Model，DEM）密切相關(guān)，地面的起伏不可避免地會對胸徑量測高度的起點產(chǎn)生影響。為了統(tǒng)一點云數(shù)據(jù)中樹木的量測起點，本文對原始點云進(jìn)行高程歸一化，以消除地形帶來的影響。

1.1.1 地面模型生成

點云數(shù)據(jù)的高程歸一化關(guān)鍵在于獲取原始點云地面數(shù)據(jù)的DEM。本文選取CSF算法提取樣地點云地面DEM。該算法是一種提取點云地面模型的算法，其原理是先將點云數(shù)據(jù)反轉(zhuǎn)，再模擬出一塊布料，將布料置于點云上方，待其自然落下，布料的形狀即為網(wǎng)格模型。由于該算法不直接處理非地面點的數(shù)據(jù)，因此運算速率較高，適合用來生成點云數(shù)據(jù)的DEM；CSF的運行時間受輸入數(shù)據(jù)的點數(shù)影響較小，可直接將原始數(shù)據(jù)作為輸入。此外，該算法參數(shù)設(shè)置簡單，只需要通過調(diào)整布料柔軟度、運算迭代次數(shù)、生成柵格大小等參數(shù)即可生成不同的地面模型。為了提取更加準(zhǔn)確的DEM，本文設(shè)置CSF的參數(shù)為：布料柔軟程度為2，進(jìn)行陡坡處理，迭代次數(shù)為500次，柵格大小為0.05 m。

1.1.2 高程歸一化

通過CSF算法獲取點云的DEM后，先將原始點云投影到DEM上，獲取每個點投影后的高度，再對原始點云中所有點減去其投影后的高度，重新生成點云數(shù)據(jù)，即高程歸一化后的點云數(shù)據(jù)。高程歸一化前后的點云效果如圖1所示。

圖1 高程歸一化前后的點云Fig.1 Point cloud before and after elevation normalization

1.2 自動化胸徑量測

1.2.1 點云切片

由于樹木胸徑是1.3 m處的樹干直徑，為了獲取有效的點云數(shù)據(jù)，本文通過切片的方式提取歸一化后點云數(shù)據(jù)中地面高度1.2～1.4 m處的所有點。通過點云切片，可大大減少點云數(shù)據(jù)的處理量，切片區(qū)域的點云數(shù)據(jù)如圖2所示。

1.2.2 點云聚類

為了提高點云聚類時的運算效率，首先對切片后的點云數(shù)據(jù)使用分辨率為0.5 cm的體素網(wǎng)格進(jìn)行體素化處理；再對點云進(jìn)行歐式距離聚類，獲得若干個點云簇；歐式聚類的距離閾值設(shè)置為0.1 m。在歐式距離聚類的過程中，構(gòu)建八叉樹來加速聚類。通過八叉樹的結(jié)構(gòu)，為每個點搜索臨近點，可提升計算效率。由于聚類后仍會存在少數(shù)雜點，可通過設(shè)置點數(shù)閾值，剔除點數(shù)小于100的點云簇，保留點數(shù)多的點云簇作為潛在的樹木對象。

圖2 點云切片F(xiàn)ig.2 Point cloud slicing

1.2.3 圓柱擬合

水平圓與圓柱擬合是常見的樹干模型擬合方式，但由于樹木并不總是垂直于地面，所以投影得到的結(jié)果無法很好地擬合樹干。如果將圓投影到與樹干軸線垂直的方向來計算樹木胸徑，則會增加冗余計算，而用圓柱體估計模型，則可以同時獲取樹干軸線方向與半徑。因此，本文采用圓柱擬合的方式以估計樹木胸徑，并對每個點云簇使用基于RANSAC的圓柱擬合方法來獲取圓柱的參數(shù)。點云簇擬合圓柱模型如圖3所示。

圖3 點云簇擬合圓柱模型Fig.3 Point cloud cluster cylinder model fitting

基于RANSAC的圓柱擬合首先計算輸入點云中每個點的法線，再選用任意兩點及其法線計算擬合圓柱模型，最后將點云中到該圓柱模型的距離閾值以內(nèi)的數(shù)據(jù)點設(shè)定為內(nèi)點，記錄內(nèi)點數(shù)量，完成第一次迭代。通過不斷迭代，將保留最大內(nèi)點數(shù)的模型作為最佳圓柱模型。由于基于RANSAC的模型估計方法僅能估計得到一種模型結(jié)果，當(dāng)數(shù)據(jù)存在兩個模型時則無法迭代收斂，因此本文采用先聚類得到點云簇，再將點云簇分別進(jìn)行RANSAC圓柱擬合的方式來避免這一缺陷。此外，如果按順序依次對點云簇進(jìn)行圓柱擬合，并不高效；而分別對點云簇執(zhí)行圓柱擬合過程，相互之間并不會產(chǎn)生影響，可以采用多線程的并行計算加速這一過程，提高圓柱擬合的效率。本文通過這種方式，快速檢測出每個點云簇中的類圓柱體，將圓柱模型半徑的2倍視為樹木胸徑估測值。

2 實驗與分析

本文基于PCL開源庫使用C++語言編程實現(xiàn)整個算法，并且基于OpenMP方法實現(xiàn)并行加速運算過程。實驗環(huán)境為一臺Intel Core i7-10750H CPU＠2.60 GHz、32 GB運行內(nèi)存的計算機(jī)。

2.1 實驗數(shù)據(jù)

本文實驗數(shù)據(jù)為FGI開源基準(zhǔn)數(shù)據(jù)集中的兩塊樣地點云數(shù)據(jù)（圖4），其點數(shù)均超過1億點。該數(shù)據(jù)集收集于芬蘭埃沃附近的森林，植被以蘇格蘭松、挪威云杉為主。每塊樣地面積為32×32 m，其中包含2014年5～7月所有DBH大于5cm的樹木量測數(shù)據(jù)。

圖4 實驗數(shù)據(jù)Fig.4 Experiment data

該數(shù)據(jù)集使用Leica HDS6100地面激光掃描儀進(jìn)行掃描，其量測視野為360°×310°，距離掃描儀25 m處的距離測量精度為±2 mm。該數(shù)據(jù)集從樣地中的5個位置進(jìn)行掃描：中心掃描一次，中心掃描的東北、東南、西南、西北方向掃描4次。邊緣掃描與中心掃描之間的距離約12 m。通過標(biāo)志物，將5站掃描數(shù)據(jù)進(jìn)行配準(zhǔn)，使其中心區(qū)域的樹木點云更加完整，極大地減少了樹木相互遮擋造成的數(shù)據(jù)缺失現(xiàn)象。

2.2 實驗設(shè)計與精度評定

本文將實驗數(shù)據(jù)中的兩塊樣地作為數(shù)據(jù)輸入，經(jīng)過算法處理后得到樹木胸徑預(yù)估值后，與實測值進(jìn)行比較分析。此外，將本文算法與Computree軟件的測算結(jié)果進(jìn)行了對比分析。Computree是一款面向林業(yè)點云數(shù)據(jù)處理的軟件，其中包含類似的樹木胸徑量測功能。因此，本文使用Computree軟件在兩塊樣地中進(jìn)行樹木胸徑值的估測，并與實測值進(jìn)行比較分析。由于該軟件無法直接處理點數(shù)較大的原始數(shù)據(jù)，因此在使用過程中先對樣地數(shù)據(jù)執(zhí)行下采樣過程，將數(shù)據(jù)下采樣至空間點間距為1cm，輸入點數(shù)約2000萬，軟件運算時均使用其默認(rèn)參數(shù)值。

本文采用均方根誤差（）與回歸方程決定系數(shù)（），評定兩種方法對樹木胸徑量測結(jié)果的精度。計算公式如下：

用來衡量實測值與估測值之間的偏差，該數(shù)值越小，表明預(yù)估值越準(zhǔn)確，精度越高。表示估測值與實測值的接近程度，該數(shù)值越接近1，則估測值越接近實測值。

2.3 實驗結(jié)果分析

2.3.1 胸徑量測數(shù)量結(jié)果

兩種方法量測與漏檢的樹木胸徑數(shù)量，以及樣地實測樹木胸徑的數(shù)量以及樹木胸徑數(shù)量的檢測率見表1。從表1可以看出，Computree軟件與本文方法量測胸徑有效值均超過80%，都做到了樹木胸徑的有效量測。其中，本文方法在樣地2中胸徑檢測率最高，Computree方法量測率均低于本文方法。經(jīng)過檢查，實測數(shù)據(jù)樣地中存在部分林下灌木，掃描儀無法掃描到樹干部位，僅能掃描其樹冠處點云，從而造成缺少樹干胸徑點云無法量測。通過分析還發(fā)現(xiàn)雖然兩塊樣地均是5站數(shù)據(jù)配準(zhǔn)后的結(jié)果，但是邊緣處仍存在部分樹木胸徑處點云呈半月狀殘缺，導(dǎo)致漏檢情況出現(xiàn)。此外，Computree軟件計算結(jié)果中樣地1中一棵樹木、樣地2中兩棵樹木與實測值相差過大，將其視為粗差歸為漏檢樹木胸徑值中。

表1 兩種方法樹木胸徑估測數(shù)量分析Tab.1 Analysis of the estimation number of trees DBH by two methods

2.3.2 胸徑量測結(jié)果

將Computree軟件與本文算法計算的樹木胸徑估測值分別按上文評價標(biāo)準(zhǔn)進(jìn)行精度評定，兩種方法在不同樣地中的回歸分析結(jié)果如圖5所示。從圖5中可以發(fā)現(xiàn)，兩種方法在兩塊樣地中的計算結(jié)果均接近于1，均小于3 cm，有較好的預(yù)估效果。本文算法在樣地2中的估測結(jié)果更好，Computree軟件在樣地1上估測結(jié)果最差。此外本文算法的結(jié)果相較于Computree軟件計算結(jié)果的更加接近于1，數(shù)值更小。

2.3.3 胸徑量測誤差分析

本文將兩種算法估測值與實測值相減再取絕對值的結(jié)果視為單棵樹木的量測誤差，圖6展示了兩塊樣地胸徑量測誤差的箱型圖。從圖6中可以看出，兩種方法的估測結(jié)果誤差主體部分均位于0～2 cm之間，少數(shù)異常值分布于3～8 cm。Computree軟件在兩塊樣地中的最大誤差均大于本文方法。其中，樣地1中的誤差均值與本文相近，樣地2的誤差均值大于本文方法。

2.3.4 時間性能分析

除精度外，算法的運行時間也是一個關(guān)鍵的衡量指標(biāo)。由于本文算法采用了多線程加速計算，在樣地1中進(jìn)行了單線程與多線程的時間測試，統(tǒng)計了1～6線程算法運行5次所消耗時間的平均值，如圖7所示。從圖7中可以看出，使用單線程計算時，算法執(zhí)行時間最慢（17.36 s）；而多線程運算均產(chǎn)生了時間效率的提升，其中6線程計算耗時最短（10.64 s），比單線程縮短了接近7 s的時間，由此證明并行運算大幅提升了時間效率。

圖5 估測DBH與實測DBH的散點圖及其關(guān)系Fig.5 Scatter plot and relationship between the estimated DBH and true DBH

圖6 DBH誤差箱型圖Fig.6 Box plot of the DBH error

圖7 運行時間與線程數(shù)的關(guān)系Fig.7 The relationship between execution time and the number of threads

此外，本文對兩種方法的運行時間進(jìn)行了統(tǒng)計，兩塊樣地中消耗時間如圖8所示。即便Computree軟件的輸入數(shù)據(jù)是下采樣后的點云數(shù)據(jù)，但在兩塊樣地中仍消耗時間超過2500 s。本文方法采用6線程并行計算，其量測樹木胸徑平均耗時10 s左右，在樹木胸徑量測的時間效率方面展現(xiàn)了一定的潛力。

圖8 兩種方法的運行時間Fig.8 Execution time of two methods

3 結(jié)束語

目前，三維激光掃描技術(shù)在森林調(diào)查工作中的應(yīng)用越來越多，基于地面激光點云展開對樹木胸徑的量測工作具有較強(qiáng)的實際價值與研究意義。本文提出了一套端到端的自動化樹木胸徑提取算法，能夠以多站配準(zhǔn)點云數(shù)據(jù)作為輸入，快速精確地獲取樹木胸徑值。通過樹木胸徑量測對比實驗，表明本文算法不僅在精度上優(yōu)于Computree軟件中的算法，又具有明顯的時間性能優(yōu)勢。針對多站配準(zhǔn)的、數(shù)據(jù)量較大的林業(yè)數(shù)據(jù)，該方法很好地解決了樹木胸徑量測的時間成本問題。然而，本文算法也存在許多需要改進(jìn)的方面，目前該算法面向的林業(yè)場景僅適用于樹干明顯的數(shù)據(jù)，對于樹干嚴(yán)重遮擋的森林還需要進(jìn)一步探索。