基于單站三維激光掃描數(shù)據(jù)的林木胸徑提取方法

柯 敏，解 靜，劉 昔

（天津水運(yùn)工程勘察設(shè)計(jì)院 天津市水運(yùn)工程測(cè)繪技術(shù)重點(diǎn)實(shí)驗(yàn)室，天津 300456）

胸徑又稱干徑，指喬木主干離地表面1.3 m處的直徑，斷面畸形時(shí)，測(cè)取最大值和最小值的平均值[1]。在綠化工程投資審計(jì)過(guò)程中，林木胸徑大小作為審核重點(diǎn)之一，其精度直接影響工程造價(jià)[2]。傳統(tǒng)的胸徑測(cè)量方法主要是用標(biāo)準(zhǔn)測(cè)量工具，如測(cè)樹(shù)胸徑尺、游標(biāo)卡尺等，在地面向上1.3 m處測(cè)量樹(shù)木的直徑。該方法效率低下，數(shù)據(jù)獲取周期長(zhǎng)，時(shí)效性差，且勞動(dòng)強(qiáng)度較大，不能滿足實(shí)際調(diào)查需求[3-4]。

為快速獲取林木胸徑數(shù)據(jù)，近年來(lái)國(guó)內(nèi)技術(shù)人員嘗試?yán)冒ㄈS激光掃描技術(shù)在內(nèi)的許多技術(shù)進(jìn)行試驗(yàn)研究。在技術(shù)應(yīng)用方面，關(guān)曉平等[5]利用電子全站儀，應(yīng)用三角高程法及截面圓法分別計(jì)算樹(shù)高及胸徑，成果精度可靠。曹先革等[6]通過(guò)分析影響地面三維激光掃描精度的因素及誤差，提出外業(yè)時(shí)掃描儀距離被掃描物體之間的距離盡量保持在50 m以內(nèi)，點(diǎn)云數(shù)據(jù)拼接應(yīng)優(yōu)先選擇標(biāo)靶拼接。在算法研究方面，樊麗等[7]利用二維激光掃描儀獲取單木胸徑處點(diǎn)云數(shù)據(jù)，分別采用最小二乘圓擬合算法、Nlinfit算法、Taubin算法、Pratt算法對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行擬合處理，得出基于最小方差迭代圓擬合算法結(jié)果精度最高，算法穩(wěn)定。孫浩等[8]利用二維激光掃描5種不同管徑的PVC管，應(yīng)用幾何法與擬合法，結(jié)合4種角度補(bǔ)償算法計(jì)算管徑并分析精度，表明擬合法在所有計(jì)算方法中精度與穩(wěn)定性最高。劉偉樂(lè)等[9]利用三維激光掃描技術(shù)對(duì)少量楊屬Populus植物進(jìn)行多站掃描、配準(zhǔn)后分析，提出林木胸徑點(diǎn)云最理想的厚度為1～2 cm，采用快速凸包算法提取胸徑效果最好。王祺等[10]對(duì)采用三維激光掃描技術(shù)獲得的樹(shù)木胸徑點(diǎn)云數(shù)據(jù)，利用蟻群算法和B樣條曲線擬合技術(shù)計(jì)算樹(shù)木胸徑，相對(duì)傳統(tǒng)算法精度取得較大提高。王寧寧等[11]利用三維掃描技術(shù)獲取高密度人工林中楊樹(shù)的冠形特征數(shù)據(jù)，通過(guò)提取胸徑處截面面積，利用圓面積公式計(jì)算胸徑值并進(jìn)行回歸分析，結(jié)果表明測(cè)量精度可靠。李丹等[12]利用多站TSL數(shù)據(jù)通過(guò)配準(zhǔn)、過(guò)濾、抽稀等處理，生成灰度影像，結(jié)合Hough變換和橢圓擬合方法提取單木胸徑等參數(shù)，并與外業(yè)實(shí)際測(cè)量數(shù)據(jù)對(duì)比，結(jié)果一致性較高。

上述研究中基于點(diǎn)云數(shù)據(jù)進(jìn)行林木胸徑提取主要包括以下特點(diǎn)：（1）多為利用標(biāo)靶裝置，通過(guò)對(duì)多站激光掃描數(shù)據(jù)進(jìn)行拼接后，提取胸徑處較為完整的截面數(shù)據(jù)進(jìn)行胸徑提取。（2）研究林木樣本數(shù)量較少，數(shù)據(jù)量較小，因而可利用相關(guān)算法快速進(jìn)行最優(yōu)解研究，保證林木胸徑提取精度。對(duì)于城市綠化林木，主要位于道路兩側(cè)，與森林、公園等林木呈片狀分布不同，主要呈條帶狀密集分布。因而，城市綠化林木胸徑調(diào)查工作具有林木數(shù)量多、密度大、通視條件差、調(diào)查精度要求高等特點(diǎn)，林木胸徑提取難度較大。如采用上述方法進(jìn)行胸徑提取，容易造成標(biāo)靶數(shù)量設(shè)置過(guò)多、利用率低，從而影響調(diào)查工作進(jìn)度，同時(shí)對(duì)于獲取的海量數(shù)據(jù)進(jìn)行提取及處理，也比較復(fù)雜及困難。

針對(duì)綠化林木的特點(diǎn)，通過(guò)現(xiàn)場(chǎng)試驗(yàn)，采用單站獲取的三維激光掃描數(shù)據(jù)，利用最小二乘圓擬合算法進(jìn)行胸徑快速提取，并對(duì)提取精度進(jìn)行分析統(tǒng)計(jì)，擬解決傳統(tǒng)作業(yè)方法效率低、時(shí)效差、周期長(zhǎng)的問(wèn)題，優(yōu)化多站數(shù)據(jù)拼接作業(yè)方法數(shù)據(jù)量龐大、處理難度高的問(wèn)題，實(shí)現(xiàn)快速提取綠化林木胸徑的目的，以期為城市綠化林木調(diào)查統(tǒng)計(jì)分析奠定基礎(chǔ)。

1 試驗(yàn)方法與數(shù)據(jù)采集

1.1 采集方式選取

目前，搭載三維激光掃描設(shè)備的平臺(tái)有多種，常見(jiàn)的包括機(jī)載測(cè)量平臺(tái)、車載測(cè)量平臺(tái)、地面設(shè)站測(cè)量平臺(tái)及背包測(cè)量平臺(tái)。根據(jù)三維激光掃描技術(shù)定位方式及掃描距離、精度等技術(shù)指標(biāo)的不同，目前可應(yīng)用于綠化林木測(cè)量的平臺(tái)主要包括機(jī)載測(cè)量平臺(tái)、背包測(cè)量平臺(tái)及設(shè)站測(cè)量平臺(tái)，其定位方式及技術(shù)指標(biāo)見(jiàn)下表1。其中表1所列精度為激光掃描儀相對(duì)于某一特定距離條件下精度，對(duì)于脈沖式掃描儀一般標(biāo)定距離為100 m。

表1 激光掃描測(cè)量平臺(tái)參數(shù)Table1 Parameter comparison table of laser scanning measurement platform

三維激光掃描測(cè)量誤差主要包括：定位誤差、測(cè)距誤差、測(cè)角誤差、集成誤差及其他誤差[13]。機(jī)載測(cè)量平臺(tái)及背包測(cè)量平臺(tái)屬于移動(dòng)測(cè)量平臺(tái)，其位置及姿態(tài)處于不斷運(yùn)動(dòng)變化過(guò)程中，平臺(tái)內(nèi)部集成了慣導(dǎo)系統(tǒng)，受慣導(dǎo)精度、定位精度及系統(tǒng)集成精度等因素綜合影響，測(cè)量精度均低于設(shè)站測(cè)量平臺(tái)。為滿足綠化林木胸徑測(cè)量精度要求，本試驗(yàn)選用點(diǎn)位精度更高的地面設(shè)站式三維激光測(cè)量平臺(tái)[13-15]。

1.2 現(xiàn)場(chǎng)數(shù)據(jù)采集

試驗(yàn)區(qū)為天津市市政綠化區(qū)域，位于濱海新區(qū)海濱大道臨港工業(yè)區(qū)段西側(cè)，整個(gè)區(qū)域呈帶狀分布，南北長(zhǎng)約1 km，東西寬約20 m。區(qū)域內(nèi)林木眾多、分布密集，其中西側(cè)為喬木，東側(cè)為灌木，本次僅對(duì)喬木進(jìn)行數(shù)據(jù)采集。喬木間距約為2.5 m×2.5 m，樹(shù)種為‘107楊’Populus×euramericana‘74/76’及刺槐Robinia pseudoacacia，種植時(shí)間均為2013年4月，所用苗木樹(shù)齡約為3 a。喬木自西向東呈排狀分布，其中西側(cè)3排為‘107楊’，其余5排為刺槐[16]。

2018年4月，用澳大利亞Maptek公司的Maptek I-Site 8820型三維激光掃描儀、美國(guó)Trimble公司的Trimble R8 GPS雙頻接收機(jī)采集數(shù)據(jù)。根據(jù)測(cè)區(qū)林木密集程度，從測(cè)區(qū)北側(cè)開(kāi)始，設(shè)置每測(cè)站間距為50 m，通過(guò)逐次設(shè)站掃描獲取數(shù)據(jù)，外業(yè)三維激光掃描作業(yè)步驟如下：

（1）將三維激光掃描儀架設(shè)在設(shè)計(jì)站位，對(duì)儀器進(jìn)行整平，將GPS RTK天線安裝在掃描儀指定位置。

（2）利用BHCORS提供的網(wǎng)絡(luò)RTK服務(wù)，分別獲取掃描設(shè)站點(diǎn)及后視點(diǎn)的三維坐標(biāo)信息。

（3）將三維激光掃描儀照準(zhǔn)后視點(diǎn)，并將掃描設(shè)站點(diǎn)及后視點(diǎn)坐標(biāo)信息輸入三維激光掃描儀控制手簿中，開(kāi)始掃描測(cè)量。掃描區(qū)域角度設(shè)置為360°，掃描質(zhì)量設(shè)置為中檔。

（4）本站掃描完成后，檢查數(shù)據(jù)完整性，如有缺失則補(bǔ)充掃描，如沒(méi)有則搬站至下一站位掃描，直至完成整個(gè)測(cè)區(qū)掃描作業(yè)。

1.3 數(shù)據(jù)預(yù)處理

激光點(diǎn)云數(shù)據(jù)預(yù)處理包括兩部分，一是對(duì)各測(cè)站掃描數(shù)據(jù)進(jìn)行裁切、拼接，二是對(duì)點(diǎn)云數(shù)據(jù)進(jìn)行過(guò)濾處理。

由于采用單站掃描點(diǎn)云數(shù)據(jù)進(jìn)行胸徑測(cè)量，為保證掃描點(diǎn)云密度及精度，以外業(yè)掃描相鄰兩站點(diǎn)連線的中垂線為界線，與綠化林木東西兩側(cè)邊界構(gòu)成的封閉區(qū)域?yàn)閱握军c(diǎn)云數(shù)據(jù)裁切范圍，利用處理軟件剔除范圍以外點(diǎn)云數(shù)據(jù)，得到尺寸為50 m×20 m的矩形單站點(diǎn)云數(shù)據(jù)。然后將各測(cè)站經(jīng)裁切后的點(diǎn)云數(shù)據(jù)，利用設(shè)站點(diǎn)坐標(biāo)進(jìn)行拼接，得到整個(gè)實(shí)驗(yàn)區(qū)域點(diǎn)云數(shù)據(jù)。

對(duì)試驗(yàn)區(qū)域點(diǎn)云數(shù)據(jù)進(jìn)行地面點(diǎn)分類，剔除地表以上點(diǎn)云數(shù)據(jù)，獲得地面點(diǎn)云數(shù)據(jù)。以地面點(diǎn)云數(shù)據(jù)構(gòu)建數(shù)字地面模型（DTM，Digital Terrain Model），以DTM依據(jù)，截取其垂直方向上胸徑處指定厚度的點(diǎn)云數(shù)據(jù)，過(guò)濾其他點(diǎn)云數(shù)據(jù)。胸徑處點(diǎn)云數(shù)據(jù)的截取厚度分為0 cm，2 cm，4 cm及6 cm共4個(gè)等級(jí)，提取數(shù)量不同的點(diǎn)云數(shù)據(jù)，分析胸徑測(cè)量的精度，結(jié)果見(jiàn)圖1。由圖1可知，截取厚度為0 cm且距離超過(guò)30 m時(shí)，因點(diǎn)云數(shù)量不夠無(wú)法計(jì)算胸徑值，其余截取厚度情況下均可計(jì)算胸徑值。根據(jù)掃描距離、點(diǎn)云數(shù)量及胸徑計(jì)算精度，胸徑點(diǎn)云最適宜的厚度為4 cm。

圖1 胸徑點(diǎn)云數(shù)據(jù)不同截取厚度Figure1 Different thickness values for intercepted DBH point cloud data of DBH point cloud data

2 胸徑計(jì)算與精度統(tǒng)計(jì)與分析

2.1 最小二乘圓擬合算法

假設(shè)（x0，y0）為圓心坐標(biāo)，R為圓半徑，則圓方程展開(kāi)后可表示為：

令a=-2x0，b=-2y0，c=x02+y02-R2，則式（1）可以變換為：

求解出參數(shù)a，b，c，便可得到圓心坐標(biāo)及半徑值。假設(shè)圓弧上有限的采集點(diǎn)坐標(biāo)為（xi，yi），采集點(diǎn)至圓心的距離為di，差值為ΔI=di2-R2（i=1，…，n）則有：

令F（a，b，c）為△i2之和，則有：

最小二乘圓擬合可表示為求F（a，b，c）最小。對(duì)a，b，c分別求偏導(dǎo)，令偏導(dǎo)等于0，分析比較極值點(diǎn)的函數(shù)值可求得最小值，公式如下：

求解式（5）可計(jì)算出參數(shù)a，b，c值，進(jìn)而得到圓心坐標(biāo)（x0，y0）及圓半徑。

2.2 林木胸徑計(jì)算

將提取的原胸徑點(diǎn)云數(shù)據(jù)導(dǎo)入Auto CAD軟件中，剔除其他干擾點(diǎn)后導(dǎo)出保存為DWG格式文件。利用LISP語(yǔ)言編制基于最小二乘圓擬合算法的程序，對(duì)各株林木胸徑點(diǎn)云數(shù)據(jù)進(jìn)行計(jì)算、統(tǒng)計(jì)后，得到各株林木胸徑處圓心坐標(biāo)及胸徑值。胸徑提取見(jiàn)圖2。

圖2 胸徑提取示意圖Figure2 The schematic diagram of DBH calculation

2.3 精度統(tǒng)計(jì)與分析

為驗(yàn)證本方法計(jì)算的林木胸徑的準(zhǔn)確性，對(duì)實(shí)驗(yàn)區(qū)域林木胸徑現(xiàn)場(chǎng)抽樣進(jìn)行人工測(cè)量，將三維激光掃描點(diǎn)云提取的與現(xiàn)場(chǎng)測(cè)量的林木胸徑值進(jìn)行比對(duì)，并記錄相關(guān)檢測(cè)結(jié)果。依據(jù)規(guī)范要求[17]，共隨機(jī)抽測(cè)104株，其中‘107楊’39株，刺槐65株，基本覆蓋整個(gè)實(shí)驗(yàn)區(qū)域，胸徑計(jì)算值與現(xiàn)場(chǎng)量測(cè)值差值統(tǒng)計(jì)見(jiàn)表2。

在胸徑計(jì)算值與量測(cè)值對(duì)比的104組數(shù)據(jù)中，差值大于2 cm的有5株，＞1～2 cm的17株，≤1 cm的82株，粗差率不足5%，誤差分布基本合理?，F(xiàn)場(chǎng)對(duì)林木進(jìn)行查看，發(fā)現(xiàn)該5株粗差林木均在西側(cè)樹(shù)皮嚴(yán)重缺失，因點(diǎn)云數(shù)據(jù)不完整造成差值較大。其余林木差值均小于2 cm，主要由儀器測(cè)量誤差及胸徑算法精度造成。不考慮上述5株林木因樹(shù)皮缺失造成數(shù)據(jù)差值較大，對(duì)剩余99株林木胸徑計(jì)算值按樹(shù)種分別進(jìn)行線性回歸分析，結(jié)果見(jiàn)圖3。

表2 胸徑計(jì)算值與量測(cè)值差值統(tǒng)計(jì)Table2 DBHby calculationandmeasurement

圖3 胸徑計(jì)算精度回歸分析Figure3 Regression analysis on accuracy of DBH calculation

根據(jù)測(cè)量誤差理論，分別對(duì)上述隨機(jī)檢測(cè)的‘107楊’及刺槐胸徑測(cè)量差值進(jìn)行精度統(tǒng)計(jì)，三維激光掃描儀計(jì)算的林木胸徑值中誤差σ為：

經(jīng)計(jì)算，‘107楊’計(jì)算中誤差σ1為0.8 cm，刺槐計(jì)算中誤差σ2為0.7 cm，整體計(jì)算中誤差σ為0.7 cm。從林木胸徑線性回歸分析結(jié)果及中誤差計(jì)算結(jié)果可知，本次刺槐與‘107楊’胸徑計(jì)算精度基本一致，整體計(jì)算中誤差優(yōu)于1 cm，精度較為理想。

3 結(jié)論與討論

通過(guò)利用計(jì)算機(jī)程序?qū)Φ孛嫒S激光掃描儀獲取的單站點(diǎn)云數(shù)據(jù)進(jìn)行綠化林木胸徑提取與精度分析，得到主要結(jié)論如下：

（1）利用設(shè)站式三維激光掃描儀對(duì)綠化林木進(jìn)行現(xiàn)場(chǎng)掃描，采用基于最小二乘圓擬合算法對(duì)單站掃描點(diǎn)云數(shù)據(jù)快速進(jìn)行胸徑計(jì)算，方法可行、精度可靠，可推廣應(yīng)用于實(shí)際綠化林木審計(jì)、調(diào)查項(xiàng)目中。

（2）傳統(tǒng)林木胸徑調(diào)查多采用人工量測(cè)方式，效率較低、數(shù)據(jù)查驗(yàn)繁瑣；現(xiàn)階段采用三維激光掃描調(diào)查多為利用標(biāo)靶球等對(duì)多站掃描數(shù)據(jù)進(jìn)行拼接、提取、計(jì)算，數(shù)據(jù)量龐大，處理難度大；本文采用單站掃描數(shù)據(jù)進(jìn)行胸徑提取，相較前兩種方式，數(shù)據(jù)量少、效率高、可隨時(shí)查驗(yàn)。

（3）分析本項(xiàng)目存在粗差的5株樹(shù)木情況，現(xiàn)場(chǎng)調(diào)查發(fā)現(xiàn)對(duì)應(yīng)樹(shù)木樹(shù)皮在某一方向存在較大缺失現(xiàn)象，且現(xiàn)場(chǎng)受林木密集、物體遮擋等因素影響，造成該林木胸徑處獲取的點(diǎn)云數(shù)據(jù)不足，無(wú)法辨識(shí)胸徑截面在不同方向上存在較大的差異，最終影響計(jì)算精度。后續(xù)將進(jìn)一步研究相關(guān)解決方案。