基于林地激光點云的樹木胸徑自動提取方法

吳杭彬，王旭飛，劉 春

（同濟(jì)大學(xué)測繪與地理信息學(xué)院，上海 200092）

森林生態(tài)系統(tǒng)是陸地生態(tài)的一項關(guān)鍵組成成分，定期進(jìn)行森林資源調(diào)查對其發(fā)展與保護(hù)有重要意義。林業(yè)資源調(diào)查中，胸徑通常指樹木距離地面1.3m 處的直徑大小，它是樹木的一項基礎(chǔ)參數(shù)，不僅與樹木冠幅、材積有極大的相關(guān)性，還可以用于推測樹木蓄積量、生物量等其他關(guān)鍵林分參數(shù)［1］。傳統(tǒng)的樹木胸徑參數(shù)常采取人工方式進(jìn)行調(diào)查量測，借助圍尺或卡尺等工具［2］，對調(diào)查森林樣地中的每一棵樹木依次量測與記錄，但這種調(diào)查方式耗費較多的人力、物力。近些年，一些新興的量測設(shè)備運用于胸徑量測中，例如電子經(jīng)緯儀［3］、全站儀［4］、CCD超站儀［5］等，但這些設(shè)備的量測方式與傳統(tǒng)方式類似，仍需要一棵棵樹木進(jìn)行量測，無法實現(xiàn)樹木批量自動化的胸徑提取。

三維激光掃描技術(shù)近年來發(fā)展迅速，激光雷達(dá)設(shè)備能以非接觸的方式獲取林木結(jié)構(gòu)參數(shù)，在林業(yè)與生態(tài)領(lǐng)域得到了廣泛應(yīng)用［6-7］。其中，機(jī)載激光雷達(dá)具備獲取大范圍森林資源結(jié)構(gòu)的能力［8］，但更多的應(yīng)用于林上結(jié)構(gòu)參數(shù)反演；而地面激光掃描技術(shù)，其精準(zhǔn)性、無損性為林下結(jié)構(gòu)的測量提供了新的技術(shù)手段［9］。近年來，許多學(xué)者展開了基于TLS 數(shù)據(jù)的樹木胸徑參數(shù)提取研究。其中，許多研究基于Hough 變換的算法提取了樹木胸徑參數(shù)，其提取步驟一般是：先對林區(qū)點云進(jìn)行高程歸一化，再對點云進(jìn)行分層并柵格化，使用Hough 變換算法對每層點云進(jìn)行圓的識別，擬合二維圓形并增加判斷條件檢測點云中的樹木位置及胸徑值［10］。但是該類方法將點云數(shù)據(jù)轉(zhuǎn)換為柵格數(shù)據(jù)時不可避免地會產(chǎn)生精度的損失，從而影響最終的提取結(jié)果。Calders 等［11］與Tansey等［12］則采用基于最小二乘圓擬合的方式對提取的單樹點云樹木胸徑進(jìn)行了計算量測。但他們的方法需要對林地數(shù)據(jù)實現(xiàn)單木分割的預(yù)處理步驟，這增加了額外的計算處理時間。并且，他們以單木分割結(jié)果的最低點為胸徑高度的量測起點，這降低了胸徑高度的計算精度。此外，Olofsson 等［13］改進(jìn)了基于RANSAC模型對二維圓的擬合算法，并以此提取樹木胸徑值。然而當(dāng)樹木存在傾斜角度時，以上所述二維圓擬合的方法則會存在進(jìn)行二維投影時數(shù)據(jù)冗余、模型誤差增大的缺陷。

針對上述缺陷，本文提出了一種基于隨機(jī)RANSAC 模型的樹木胸徑自動提取方法。面向大量高精度的地面激光點云數(shù)據(jù)，先通過CSF 算法提取數(shù)字地面模型，并根據(jù)該模型提取樹木胸徑高度的點云層，然后使用歐式距離對點云層聚類，最后對聚類后的點云簇分別采用隨機(jī)RANSAC 模型進(jìn)行圓柱擬合，以圓柱模型的直徑作為該樹木胸徑。

1 樹木胸徑自動提取算法

1.1 樹木胸徑自動提取算法流程

為實現(xiàn)自動化的地面激光點云樹木胸徑參數(shù)提取，給出了如圖1 的總體算法流程。該算法主要步驟為：基于CSF算法的點云濾波與林地DEM構(gòu)建、基于平行DEM的樹木胸徑點云提取與聚類、基于隨機(jī)RANSAC模型的圓柱擬合與胸徑計算。

圖1 樹木胸徑量測算法流程Fig.1 Flow chart of tree DBH measurement algorithm

1.2 基于CSF算法的點云濾波與林地DEM構(gòu)建

通過地面激光掃描儀采集的激光點云數(shù)據(jù)常包含地形起伏、地表崎嶇不平等情況，樹木胸徑的高度與地形密切相關(guān)，這為準(zhǔn)確計算每棵樹木的胸徑高度帶了一定困難。為獲取樹木胸徑附近的點云即樹高1.3m 處附近的點云，則需依靠準(zhǔn)確的DEM。CSF 算法［14］是一種高效的點云地面濾波算法，它的原理是：先將處理的點云數(shù)據(jù)反轉(zhuǎn)即點云中x，y坐標(biāo)不變，z坐標(biāo)變?yōu)橄喾磾?shù)，所有點的坐標(biāo)由（x，y，z）轉(zhuǎn)變?yōu)椋▁，y，-z）。然后模擬出一塊布料模型，根據(jù)點云場景的地形來設(shè)置布料的柔軟程度，起伏越大柔軟程度設(shè)置越大，反之越平坦的地形，柔軟程度設(shè)置越小。此外，還需要設(shè)置柵格大小決定生成模型布料粒子的數(shù)量。最后，將布料置于反轉(zhuǎn)點云上方，模擬自然下落的過程，落到點云則停止下落。再對布料粒子計算其受重力影響產(chǎn)生的位移和粒子間引力、排斥力產(chǎn)生的位移，遍歷所有布料粒子，迭代終止時布料形狀即為DEM。

CSF 算法參數(shù)設(shè)置簡單，生成模型準(zhǔn)確、魯棒［15］，滿足本文生成DEM 的需求。如圖2 所示，林地激光掃描數(shù)據(jù)經(jīng)過CSF 算法處理后，得到了樹木點云、地面點云數(shù)據(jù)與該林地的DEM 結(jié)果，其中樹木點云進(jìn)一步用于后續(xù)的處理。

圖2 基于CSF算法的點云濾波結(jié)果Fig.2 Point cloud filtering results based on CSF algorithm

1.3 基于平行DEM的樹木胸徑點云提取與聚類

在獲取林地的DEM后，考慮到胸徑計算需要在距離地面1.3m 處完成，因此本文提出了基于平行DEM 的樹木胸徑點云提取和聚類方法。這一方法避免了傳統(tǒng)胸徑計算過程中的單株樹木提取問題，也有利于減少林地中灌木的影響。先基于原始DEM構(gòu)建了平行的曲面，兩者的距離為1.3m，將曲面（設(shè)曲面厚度為d）與樹木點云進(jìn)行疊加運算，其交集可視為樹木胸徑的截面。如圖3 所示，左側(cè)為單株樹木點云，下方為構(gòu)建的林地DEM 網(wǎng)格，中間為平行DEM的曲面模型，右側(cè)為樹木胸徑點云的提取結(jié)果。

圖3 單棵樹木點云及胸徑處點云Fig.3 One tree point cloud and point cloud at DBH

提取林地中樹木胸徑處的點云后，本文采用歐式聚類算法［16］對點云聚類，將林地中不同樹木的樹干胸徑點云按單株進(jìn)行區(qū)分，便于進(jìn)行后續(xù)的樹木胸徑計算。如圖4 所示，上方為林地點云的胸徑處數(shù)據(jù)，下方為聚類后的各株樹干胸徑點云，不同灰度代表不同類別。

圖4 樹木胸徑處點云歐式聚類Fig.4 Point cloud at DBH of trees Euclidean clustering

1.4 基于隨機(jī)RANSAC圓柱擬合與胸徑計算

經(jīng)過1.2 和1.3 兩步處理獲得了每棵樹木的胸徑點云，考慮到樹干的生態(tài)生理結(jié)構(gòu)常為圓柱形且胸徑處上下直徑變化較小，因此本文采用圓柱模型對樹干進(jìn)行描述，其表達(dá)式為

式中：圓柱模型包含7 個參數(shù)，(x0，y0，z0)為圓柱軸線上一點；(α，β，γ)為圓柱軸線的單位方向向量；r為圓柱半徑。

本文對胸徑點云使用隨機(jī)RANSAC 模型［17］進(jìn)行圓柱擬合，隨機(jī)RANSAC 模型是基于RANSAC模型的一種改進(jìn)模型，不同之處在于內(nèi)點選取的過程中，沒有像RANSAC模型用全部點計算是否為內(nèi)點，而是分為兩步：先隨機(jī)選取部分點觀察其是否為內(nèi)點，再根據(jù)部分點的結(jié)果對剩余點進(jìn)行計算檢查。其圓柱模型擬合的計算步驟如下：

（1）假定一株樹木胸徑點云A（數(shù)量為n）作為圓柱模型擬合的輸入點云，先計算點云A 中所有點的法線向量并隨機(jī)選取A內(nèi)兩點，通過坐標(biāo)、法線值建立初始圓柱模型［18］。

（2）對A內(nèi)剩余點隨機(jī)選取部分點集B（數(shù)量為m，m＜＜n-2）計算其是否符合初始圓柱模型的內(nèi)點要求，如果B中所有點符合內(nèi)點要求（到擬合模型的距離滿足閾值要求），則對點云A中剩余的n-2-m個點檢查是否符合內(nèi)點要求，記錄圓柱模型參數(shù)的7 個參數(shù)與內(nèi)點數(shù)量，并隨機(jī)取兩點展開一次新的循環(huán)。如果B 中有任意一點不符合內(nèi)點要求，則跳出計算并開始下一次循環(huán)。當(dāng)循環(huán)結(jié)束后，內(nèi)點數(shù)量最多的圓柱模型即為最佳擬合圓柱模型。

隨機(jī)RANSAC模型加快了模型擬合速度，且能夠保證與標(biāo)準(zhǔn)RANSAC 算法相同的擬合精度。因此，本文采用隨機(jī)RANSAC模型分別對每一簇點云進(jìn)行圓柱擬合，見圖5，將圓柱模型的直徑作為該樹木的胸徑值。

2 實驗與分析

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

實驗區(qū)域位于上海市青浦區(qū)某區(qū)域，共計2 塊樣地，大小分別為70m×30m、60m×50m，樣地中主要樹種分別為鵝掌楸、無患子（圖6）。

圖6 實驗數(shù)據(jù)采集環(huán)境Fig.6 Experiment data environment

兩塊樣地的點云數(shù)據(jù)由Z+F 5 010C 地面激光掃描儀所獲取，該掃描儀量測視野為360°×310°，量測距離為0.3～80m，量測精度在50 m 處可達(dá)±3mm。點云數(shù)據(jù)采集時先對樣地中心進(jìn)行了單站掃描，再對中心站周圍東北、東南、西南、西北四個方向各進(jìn)行一次數(shù)據(jù)掃描［19］，中心站與四周各站平均距離約15m，現(xiàn)場共布置12 個靶球作為配準(zhǔn)標(biāo)志物。掃描結(jié)束后，在Z+F control軟件中通過標(biāo)志物將四周的掃描站配準(zhǔn)到中心站，并融合得到了完整的樣地掃描數(shù)據(jù)。多站融合的點云數(shù)據(jù)有利于減少相互遮擋帶來的影響，獲得更為完整的樹木林下結(jié)構(gòu)，兩塊樣地五站融合后的點云數(shù)量均接近兩億點。

2.2 實驗方法與參數(shù)設(shè)置

基于PCL開源庫［20］使用C++編程語言實現(xiàn)了本文算法，并基于OpenMP 第三方庫［21］在CSF 算法生成DEM 的過程與多簇點云擬合模型的過程實現(xiàn)了并行加速運算。實驗運行的硬件環(huán)境為計算機(jī)內(nèi)存32GB，CPU 型號Intel（R）Core（TM）i7-9 750H CPU@2.60GHZ。

本文按圖1 所示流程進(jìn)行了實驗，將兩塊樣地激光點云分別作為輸入，通過所提算法獲得了兩塊樣地中樹木胸徑的預(yù)估值。為驗證和對比本文算法計算結(jié)果，采用兩種對比方法。第一種方法的基本流程與本文相同，僅將圓柱擬合的步驟由隨機(jī)RANSAC的擬合方法替換為最小二乘的擬合方法。第二種方法為基于Hough 變換的算法［12］，其投影的點云高度為1.28m～1.32m。

在本文算法運行過程中，CSF算法的參數(shù)為：布料柔軟程度設(shè)置為2，進(jìn)行陡坡處理，柵格大小設(shè)置為0.5m，算法迭代次數(shù)設(shè)置為500 次。平行DEM的厚度d設(shè)置為0.1m。歐式距離聚類的搜索距離閾值設(shè)置為0.1m，約束最小點云簇的點數(shù)設(shè)置為50。

在采集數(shù)據(jù)時，現(xiàn)場使用卷尺量測并記錄了所有大于5cm的樹木胸徑數(shù)據(jù)，每棵樹木量測兩次，取量取胸徑的平均值作為該樹木的胸徑值，將其視為真值，用于各算法的比較驗證。

采用檢測率、胸徑偏差及運行時間三項指標(biāo)對三種算法進(jìn)行了評價。檢測率為算法檢測到樣地樹木占實際量測樹木數(shù)量的比率，其計算公式為

式中：n代表算法檢測樹木數(shù)量值，n?代表實際量測樹木數(shù)量值。檢測率的大小可以反應(yīng)算法的魯棒性與有效性。

胸徑偏差代表真實胸徑與不同算法計算胸徑之間差異的絕對值，其計算式為

式中：x表示算法提取樹木胸徑數(shù)值，x?代表實際量測樹木胸徑數(shù)值。偏差值的大小能夠反映算法量測精度的高低，與真實值的偏差越小，則精度越高，反之則越低。

2.3 DEM構(gòu)建與胸徑點云提取

對兩樣地點云數(shù)據(jù)執(zhí)行基于CSF 算法的點云濾波后，生成了對應(yīng)的DEM 模型，其中樣地1 的林地點云與DEM 模型如圖7 所示。圖中還給出了樣地的地面局部點云，可發(fā)現(xiàn)林地點云數(shù)據(jù)中地面地形并不平坦，而是存在起伏，且樹與樹之間存在溝壑現(xiàn)象，DEM的提取精度會直接影響胸徑自動計算的精度。

圖7 實驗點云數(shù)據(jù)、DEM及局部溝壑Fig.7 Experiment point cloud data,DEM and local ravines

得到準(zhǔn)確的DEM 模型后，構(gòu)建一個平行于DEM模型的曲面模型，兩者相距1.3m，通過曲面模型提取了1.25～1.35m處的胸徑點云層。針對該點云層，通過歐式距離聚類得到了許多不同的點云簇，聚類結(jié)果如圖8 所示，不同的灰度代表了不同的點云簇。

圖8 林地胸徑點云歐式聚類結(jié)果Fig8 Point cloud at DBH of trees Euclidean clustering results

2.4 樹木胸徑檢測數(shù)量結(jié)果

將本文算法、基于Hough變換的算法以及基于最小二乘的算法提取的胸徑數(shù)量按式（2）計算檢測率進(jìn)行評估，其結(jié)果如表1所示。結(jié)果表明，三種算法在兩塊實驗樣地中的檢測率均超過了95%，達(dá)到了有效檢測。本文提出的算法在兩塊實驗樣地中的檢測率均為100%，即檢測到了所有樹木?；贖ough變換的算法與基于最小二乘的算法在樣地1中的檢測率相同，均漏檢了2棵樹木。在樣地2中，其余兩種方法略低于本文所提算法，基于最小二乘的算法漏檢了6棵樹木。在圓柱模型擬合時部分點云簇中存在小部分噪聲點，基于最小二乘的算法可能無法將其忽略，擬合模型失敗而導(dǎo)致漏檢樹木。

表1 不同方法的樹木檢測數(shù)量結(jié)果Tab.1 The detection results of trees by different methods

2.5 樹木胸徑量測結(jié)果

對三種不同方法在兩樣地中的胸徑計算結(jié)果按式（3）分別計算偏差，并統(tǒng)計偏差的均值、最大值和最小值。

表2結(jié)果表明，三種方法的偏差最大值差異較大，其中基于最小二乘的算法在樣地2 中的最大偏差為42.07cm，為三種方法在兩塊樣地中的最大偏差，本文算法的最大偏差在兩塊樣地中均為最小。

三種方法的偏差最小值在兩塊樣地中的表現(xiàn)差異較小，均未超過0.05cm。三種方法偏差均值均小于4cm，其中基于最小二乘的算法與本文算法偏差均值均小于1cm，取得了良好的提取效果。整體來看，本文算法的量測結(jié)果優(yōu)于其余兩種算法。

進(jìn)一步分析三種方法在兩樣地中所有樹木胸徑與實測結(jié)果的偏差值，并將其從大到小排序展示，如圖9和10所示。

通過偏差結(jié)果可以發(fā)現(xiàn)，基于Hough 變換的算法在兩塊樣地中的量測誤差較其余兩種方法更大，本文圖1所示算法框架即其余兩種方法所遵循框架的量測精度更高。如圖9b、9c和圖10b、10c所示，因為最小二乘的算法與本文方法采用了相同框架，量測偏差分布較為類似，但隨機(jī)RANSAC模型擬合圓柱模型的樹木胸徑量測值結(jié)果與實際量測的偏差綜合來看更小，因此，可以證明本文方法優(yōu)于最小二乘的方法。通過比較三種方法的偏差，驗證了本文所提出的基于隨機(jī)RANSAC 模型樹木胸徑提取算法 優(yōu)于其余兩種算法。

圖9 樣地1的胸徑偏差分布Fig.9 Distribution of bias in tree DBH measurement in Plot 1

圖10 樣地2的胸徑偏差分布Fig.10 Distribution of bias in tree DBH measurement in Plot 2

此外，考慮到本文研究結(jié)果的實際應(yīng)用性，對三種算法的量測精度進(jìn)行了分析比較。按照國家標(biāo)準(zhǔn)《森林資源規(guī)劃設(shè)計調(diào)查技術(shù)規(guī)程》（GB/T 26424-2010）［22］中森林資源調(diào)查精度的胸徑允許誤差要求，本文對三種方法每棵樹木胸徑量測的誤差進(jìn)行了統(tǒng)計與評估，其結(jié)果如表3 所示。胸徑量測誤差是由上文中胸徑偏差與真實胸徑量測值的比值計算所得。表3 中等級A 代表了胸徑量測誤差小于5%的樹木占總體樹木數(shù)量的比例，等級B 代表了胸徑量測誤差為5%～10%的樹木占總體樹木數(shù)量的比例，等級C代表了胸徑量測誤差不超過10%～15%的樹木占總體樹木數(shù)量的比例，等級D 代表了胸徑量測誤差大于15%的樹木占總體樹木數(shù)量的比例，其中量測誤差小于15%的樹木符合國家標(biāo)準(zhǔn)精度等級的要求。表3 結(jié)果表明，基于Hough 變換的算法在兩樣地中樹木量測精度明顯低于其余兩種方法，量測誤差占比分布均勻，誤差偏大的樹木較多，在樣地2 中超過一半的樹木量測誤差超過了15%，整體表現(xiàn)不佳?；谧钚《说乃惴繙y精度整體較高，在兩樣地中超過3/4誤差均達(dá)到了5%的要求，但是仍有部分樹木的量測誤差超過了15%。本文算法較其余兩種算法量測精度表現(xiàn)更加良好，在兩樣地中量測誤差小于5%的樹木數(shù)量均達(dá)到最高，并且在樣地1 中所有樹木的量測誤差均小于15%，符合國家標(biāo)準(zhǔn)中對量測精度的要求；在樣地2中，超過95%的樹木量測誤差符合國家標(biāo)準(zhǔn)的最低要求，整體上本文算法保證了較高的量測精度，有較高的應(yīng)用可行性。

表3 不同方法的胸徑量測誤差等級占比Tab.3 Proportion of error grade of DBH measurement by different methods

2.6 效率比較

許多樹木胸徑的研究聚焦于胸徑的量測精度，卻很少關(guān)注方法的計算效率，而在實際應(yīng)用中算法的運行效率也是一項重要的評價指標(biāo)，因此表4 給出了三種方法在兩樣地的運行時間。時間結(jié)果表明，基于Hough 變換的算法與基于最小二乘的算法在兩塊樣地中運行時間均超過了1 分鐘，而本文所提算法均未超過30s，面對同樣的數(shù)據(jù)運行時間更短，其運算效率更高。

表4 不同方法的運行時間Tab.4 The running time of different methods

3 結(jié)語

提出了一種基于隨機(jī)RANSAC 模型的樹木胸徑提取算法，該算法包括基于CSF 算法的點云濾波與林地DEM 構(gòu)建、基于平行DEM 的樹木胸徑點云提取與聚類、基于隨機(jī)RANSAC模型的圓柱擬合與胸徑計算三個步驟。使用上海市青浦區(qū)的兩塊林區(qū)樣地點云數(shù)據(jù)對該算法進(jìn)行了驗證，與實際人工測量樹木胸徑的平均偏差分別為0.79cm 和0.52cm，95%及以上的樹木量測結(jié)果符合國家標(biāo)準(zhǔn)的量測精度要求，這一結(jié)果可支持該算法的可行性。從檢測率、實測數(shù)據(jù)的偏差結(jié)果與運行時間性能角度分析，本文算法均優(yōu)于基于Hough變換的算法與基于最小二乘的算法。本文算法不僅精度達(dá)到了要求，且運算效率高，具有良好的應(yīng)用前景。

作者貢獻(xiàn)聲明：

吳杭彬：研究方向確定、算法設(shè)計、論文修改。

王旭飛：算法設(shè)計、數(shù)據(jù)采集、數(shù)據(jù)分析、論文撰寫與修改。

劉春：數(shù)據(jù)分析、論文修改。