基于地基激光雷達(dá)的活立木材積提取算法

熊妮娜，王佳

(1.北京市城市管理研究院，北京 100028；2.北京林業(yè)大學(xué)精準(zhǔn)林業(yè)北京市重點(diǎn)實(shí)驗(yàn)室，北京 100083；3.北京林業(yè)大學(xué)測繪與3S技術(shù)中心，北京 100083)

傳統(tǒng)的立木材積主要依靠胸徑尺、測高器等工具對標(biāo)準(zhǔn)樣地進(jìn)行調(diào)查，以標(biāo)準(zhǔn)木胸徑、樹高查詢材積表進(jìn)行估測，該方法效率低、人為干擾因素多，受材積表精度影響較大[1]。目前全站儀作為立木材積精準(zhǔn)量測的儀器已經(jīng)開展了相關(guān)實(shí)驗(yàn)和應(yīng)用研究，該方法技術(shù)成熟，但同樣存在操作復(fù)雜、不能海量獲取坐標(biāo)點(diǎn)等問題。而地基激光雷達(dá)技術(shù)作為現(xiàn)代測量的新興技術(shù)手段之一，近年來被廣泛地應(yīng)用在單木計(jì)測和森林調(diào)查各個方面，激光雷達(dá)掃描可以數(shù)秒內(nèi)一次性采集到上萬點(diǎn)，在獲取目標(biāo)三維坐標(biāo)位置方面具有先天優(yōu)勢。目前在林業(yè)中應(yīng)用于森林生態(tài)系統(tǒng)、森林資源調(diào)查的各領(lǐng)域，包括樹種、樹高、胸徑、生物量等單木屬性信息的提取，葉面積指數(shù)的估算，冠層空隙度，冠層輻射，冠層結(jié)構(gòu)，葉面積分布，樹干曲線，競爭指數(shù)等[2]。隨著地基激光雷達(dá)體積逐漸減小，更加便于攜帶，同時生產(chǎn)成本的降低帶來價(jià)格的降低，地基激光雷達(dá)將會更加普遍地應(yīng)用于林業(yè)調(diào)查中。

目前國內(nèi)外學(xué)者開展了利用地基激光雷達(dá)提取單木參數(shù)的一系列研究。Thies等[3]、Astrup等[4]提出了一種基于立木的地面激光掃描數(shù)據(jù)重建樹干三維表面的方法和算法。Tansey等[5]根據(jù)地面激光掃描儀數(shù)據(jù)估算了科西嘉島松林地的樹木材積。鄧向瑞等[6]利用地基激光雷達(dá)對活立木進(jìn)行了精確掃描測量，利用軟件提取樹木的直徑、樹高、冠幅等測樹因子，利用材積模型計(jì)算立木材積，并與伐倒木的材積作對比。王劍等[7]通過激光點(diǎn)云數(shù)據(jù)構(gòu)建一系列圓柱模型片段疊加而成樹干的三維模型，該方法對樹干的體積計(jì)算也提供了一種便利。蔣佳文等[8]采用地面激光雷達(dá)進(jìn)行多站點(diǎn)掃描獲取立木的點(diǎn)云信息，提取有關(guān)點(diǎn)云分布的特征參數(shù)，建立基于特征參數(shù)的材積模型。參考國內(nèi)外學(xué)者的研究，本研究基于地基激光雷達(dá)點(diǎn)云數(shù)據(jù)，改進(jìn)Hough變換方法，通過提取不同高度處的樹干直徑，運(yùn)用斷面積區(qū)分求積法計(jì)算樹干材積，以期為無損精準(zhǔn)測木提供新的方法。

1 數(shù)據(jù)處理與獲取

1.1 外業(yè)數(shù)據(jù)采集

本研究采用FARO 120型地基激光雷達(dá)，掃描儀穩(wěn)定可靠，測量距離為0.6～120 m，分辨率為0.1 mm，數(shù)據(jù)獲取率≤50 800像素/s，25 m內(nèi)誤差≤2 mm，水平和垂直視野范圍為360°和320°。

在樹干胸高處朝南或朝北方向貼標(biāo)靶紙，選取樣木周圍無遮擋的地方均勻放置3個參考球，選擇能同時見到3個參考球的均勻分布的3個位置設(shè)站，理想的測站位置間隔角度為120°(以待掃描樹為參考)，安裝地基激光雷達(dá)并連接好相關(guān)設(shè)備，對掃描范圍、掃描分辨率等參數(shù)進(jìn)行設(shè)置后，進(jìn)行樹木掃描。本研究設(shè)置掃描區(qū)域?yàn)樗椒较?60°，垂直方向155°，掃描分辨率為10 000點(diǎn)/周，每株樣木需從不同角度掃描至少3次，耗時約10 min/株。

1.2 內(nèi)業(yè)數(shù)據(jù)處理

本研究點(diǎn)云數(shù)據(jù)處理參考張浩[9]、張冬等[10]采用的枝葉分離方法，主要包括4個步驟，結(jié)果如圖1所示：

1)點(diǎn)云數(shù)據(jù)預(yù)處理。將地基激光雷達(dá)掃描數(shù)據(jù)導(dǎo)入計(jì)算機(jī)中，通過其配套軟件Faro Scene，對點(diǎn)云數(shù)據(jù)進(jìn)行加載、配準(zhǔn)、濾波、剔除、導(dǎo)出等內(nèi)業(yè)處理。打開Faro Scene軟件，直接把測站數(shù)據(jù)加載到軟件中，第一次打開軟件時，掃描數(shù)據(jù)的Mio掃描點(diǎn)數(shù)值默認(rèn)為62，這個數(shù)值是根據(jù)計(jì)算機(jī)的內(nèi)存自動顯示的，本研究涉及的點(diǎn)云提取使用的計(jì)算機(jī)內(nèi)存為2 G，經(jīng)過在不同內(nèi)存的計(jì)算機(jī)中進(jìn)行點(diǎn)云數(shù)據(jù)提取實(shí)驗(yàn)，結(jié)果顯示，默認(rèn)的Mio掃描點(diǎn)的數(shù)值均能滿足精度要求。

2)樹木點(diǎn)云數(shù)據(jù)主方向計(jì)算。樹木點(diǎn)云數(shù)據(jù)主方向計(jì)算需要借助近鄰點(diǎn)坐標(biāo)信息及近鄰點(diǎn)相應(yīng)的法向量信息，原始的樹木點(diǎn)云只包含點(diǎn)云坐標(biāo)信息，本研究借助于主成分分析進(jìn)行點(diǎn)云法向量計(jì)算；法截線擬合法計(jì)算主方向?qū)υ肼朁c(diǎn)具有一定的免疫力，是目前計(jì)算點(diǎn)云主方向較好的一種方法，計(jì)算出的枝干部分的點(diǎn)云主方向與枝干中軸線平行，冠層部分的點(diǎn)云主方向雜亂無章。

圖1 內(nèi)業(yè)處理過程Fig. 1 The interior work process

3)樹木點(diǎn)云數(shù)據(jù)枝葉分割。首先根據(jù)點(diǎn)云數(shù)據(jù)各點(diǎn)及相應(yīng)的主方向構(gòu)建圓柱，圓柱半徑定為2倍采樣間距，圓柱高度定為10倍采樣間距。然后將圓柱包圍盒內(nèi)的點(diǎn)投影到圓柱軸線上，根據(jù)圓柱內(nèi)的點(diǎn)與圓柱軸線的夾角、投影點(diǎn)在軸線上的分布密度進(jìn)行判斷，區(qū)分冠層葉子點(diǎn)和枝干點(diǎn)。研究中分割時閾值?1設(shè)為0.7，閾值?2設(shè)為30°，閾值?3設(shè)為0.8，其結(jié)果較好。

4)枝干骨架特征點(diǎn)提取。枝干部分依據(jù)點(diǎn)云測地圖劃分水平集，基于最小二乘法對水平集進(jìn)行二維圓擬合，同時通過目視判斷對擬合結(jié)果進(jìn)行修正。

2 立木材積提取算法

2.1 任意處直徑提取算法

本研究利用自動檢測任意處直徑的方法，首先截取任意處直徑的點(diǎn)云數(shù)據(jù)水平切片，并進(jìn)行柵格化，然后再使用圓擬合結(jié)合Hough變換，這是一個特征提取技術(shù)在特定形狀不完善的情況下由一個決策程序找到對象的規(guī)則形狀的方法，最終將圓檢測出來，進(jìn)而計(jì)算直徑。該方法分為兩個步驟：

圖2 胸高部分水平切片截取算法和胸徑量測算法流程Fig. 2 The flow chart of DBH of individual trees retrieval algorithm

1)任意處直徑部分水平切片截取方法及步驟

任意處直徑的點(diǎn)云數(shù)據(jù)水平切片是從構(gòu)成單木的點(diǎn)云數(shù)據(jù)中選取高于任意處直徑0.1 m至低于直徑處0.1 m的點(diǎn)云數(shù)據(jù)。以Matlab實(shí)現(xiàn)算法，具體過程如圖2所示。單木樹干任意直徑處上下0.1 m的小切片被截取下來，將此部分點(diǎn)云數(shù)據(jù)投影到xy平面，以便進(jìn)行下一步邊緣檢測，并確定檢測圓的半徑。

2)任意直徑處量測算法及步驟

在林業(yè)測量中，常遇到直徑切面為不規(guī)則圓的情況，故本研究在不規(guī)則的直徑切面點(diǎn)云投影數(shù)據(jù)上檢測最能滿足截面輪廓的標(biāo)準(zhǔn)圓，該檢測圓的直徑即視為樹干任意處直徑。圓的檢測和半徑的計(jì)算通過Hough變換方法，由于Hough變換檢測圓程序需要讀入圖片，所以首先對水平切片投影后的數(shù)據(jù)進(jìn)行柵格化處理，結(jié)合格網(wǎng)大小，換算出真實(shí)圓的直徑[11-12]。變換圓的方程為x2+y2+2ax+2by+c=0，用兩個累加器A(a，b)和B((a2+b2-c)1/2)來進(jìn)行統(tǒng)計(jì)可能的圓心和半徑。

本研究采用了一種改進(jìn)的Hough變換算法，通過降低隨機(jī)采樣的點(diǎn)數(shù)，以及利用采樣點(diǎn)的特征信息來判斷是否進(jìn)行累積，來避免大量的無目標(biāo)的隨機(jī)采樣與無用累積[13-14]，使算法性能大大提高。具體算法如圖2所示。經(jīng)測試，將最小圓周上邊緣點(diǎn)的個數(shù)定義為2 000較為合適。由于原始的點(diǎn)云數(shù)據(jù)是5個站點(diǎn)配準(zhǔn)疊加的數(shù)據(jù)，數(shù)據(jù)量較大，如果參數(shù)設(shè)置太小，會導(dǎo)致檢測圓的效果不佳，檢測出多個小圓。第一次隨機(jī)選取3個邊緣點(diǎn)，將這3個邊緣點(diǎn)代入變換圓公式x2+y2+2ax+2by+c=0，用兩個累加器A(a，b)和B[(a2+b2-c)1/2]來進(jìn)行統(tǒng)計(jì)可能的圓心和半徑。然后清除屬于該檢測圓的所有邊緣點(diǎn)。由于這里視樹干直徑為一個標(biāo)準(zhǔn)圓，故所有小圓的圓心位置視為樹干直徑標(biāo)準(zhǔn)大圓的邊緣點(diǎn)，進(jìn)行邊緣檢測，當(dāng)某次循環(huán)結(jié)束但產(chǎn)生邊緣點(diǎn)的統(tǒng)計(jì)個數(shù)達(dá)到或超過閾值的構(gòu)造圓時，則圓檢測結(jié)束。檢測結(jié)果不是單個圓，而是以多個小圓圓心為邊緣點(diǎn)組成多個半徑一致、圓心位置相近的圓。所得結(jié)果對于獲得樹干直徑無影響，若想獲得圓心位置，將得到的多個圓圓心位置均值可看作所求圓心半徑，即為單木的位置。最后結(jié)合柵格圖像格網(wǎng)大小，得到單木的真實(shí)任意樹干處直徑。

以實(shí)驗(yàn)提取的一棵山楊為例，提取的樹干不同高度的點(diǎn)云切片見圖3a，選取胸徑處的點(diǎn)云切片數(shù)據(jù)為例進(jìn)行后面的直徑提取見圖3b，選取位置水平切片投影后的單木直徑點(diǎn)云數(shù)據(jù)分布見圖3c，柵格化后的水平切片見圖3d，利用Hough變換檢測出的圓效果見圖3e，圓心坐標(biāo)、半徑長度和真實(shí)單木直徑在Matlab軟件中可自動計(jì)算出來。對本研究的單木，其圓半徑為56，再結(jié)合預(yù)先設(shè)置的柵格格網(wǎng)大小0.1 cm，樹干直徑即為56×0.1×2=11.2 cm。

圖3 利用Hough變換檢測出圓的過程Fig. 3 Detect circle by using Hough transform

2.2 單木材積計(jì)算模型

本研究采用中央斷面積區(qū)分求積法計(jì)算單木材積，將樹干按一定長度(本研究采用2 m)分段，量出每段中央直徑和最后不足一個區(qū)分段梢頭底端直徑，利用中央斷面近似求積式，求算各分段的材積(V1，V2，V3，…，Vn)，并合計(jì)(V)。

V=V1+V2+V3+…+Vn+V′=g1l+g2l+g3l+

(1)

2.3 實(shí)驗(yàn)精度驗(yàn)證方法

對于研究立木材積精度驗(yàn)證最可靠的方法為：將點(diǎn)云數(shù)據(jù)提取的立木材積和樹木伐倒后計(jì)算的材積相對比，但由于受限于政策原因，無法對樹木進(jìn)行此類破壞性試驗(yàn)。因此，本研究利用全站儀活立木材積測量方法對實(shí)驗(yàn)樹木提取立木材積作為材積近似真值，全站儀無損測量立木的原理如下，首先用胸徑尺測量立木樹干的地徑D0和胸徑D1處(距離地面1.3 m處)的胸高h(yuǎn)1，然后用全站儀測量儀器中心到胸徑D1處的斜距S，再用全站儀將胸徑以上的樹干分成n段并測量相應(yīng)位置的水平角αi和天頂距γi(下角標(biāo)i為分段號)，后運(yùn)用三角高程原理計(jì)算樹高H，并按照圓臺累加法計(jì)算得到立木樹干的總材積V，具體方法如圖4所示。目前已有研究表明，全站儀立木材積提取方法獲取的材積與伐倒木真實(shí)材積誤差在5%以內(nèi)[15-16]，可以近似作為真值檢驗(yàn)其他方法。

注：H為樹高，m；h1, …, hi, …, hn為各分段高度，m；V為總材積，m3；V1, …, Vi, …, Vn為各分段材積，m3；D0、D1分別為地徑和胸徑，cm；D2, …, Di, …, Dn為各分段直徑，cm；α1,…, αi, …, αn為各分段處的水平角，(°)；γ1, …, γi, …, γn為各分段處的天頂距，(°)；S為儀器中心到胸徑處的斜距，m。 圖4 全站儀測量立木材積原理Fig. 4 Schematic diagram for measuring standing tree volume using total station

e=V激光V全站

(2)

(3)

(4)

式中：V全站為全站儀提取材積數(shù)據(jù)，m3；V激光為地基激光雷達(dá)提取材積數(shù)據(jù)，m3。

3 結(jié)果與分析

本研究實(shí)測共56棵山楊樹，采用現(xiàn)場直接量測的樹木胸徑做為參考值，將地基激光雷達(dá)掃描提取的胸徑與其進(jìn)行對比，其結(jié)果見表1。其中絕對誤差最大的為0.023 3 m，最小值為0.008 0 m，沒有絕對誤差的有9棵，占總數(shù)的16.07%；絕對誤差為正值的共有35棵，占總數(shù)的62.50%；為負(fù)值的為12棵，占總數(shù)的21.43%。從相對誤差來看，最大的為8.82%，最小的為0.00%，平均相對誤差為1.81%，可以看出地基激光雷達(dá)掃描提取胸徑與現(xiàn)場直接量測的樹木胸徑十分接近，精度較高。

針對實(shí)測的56棵山楊樹，運(yùn)用本研究地基激光雷達(dá)掃描提取的材積與全站儀提取材積進(jìn)行對比，其結(jié)果見表1。其中，絕對誤差最大的為 0.070 6 m3，最小值為-0.110 2 m3，絕對誤差為正值的共有37棵，占總數(shù)的66.1%，為負(fù)值的為19棵，占總數(shù)的33.9%。從相對誤差來看，最大的為15.07%，最小的為0.06%，平均相對誤差為5.86%，可以看出地基激光雷達(dá)掃描提取材積與全站儀提取材積十分接近，誤差較小。同時分析實(shí)驗(yàn)結(jié)果發(fā)現(xiàn)，地基激光雷達(dá)掃描提取的材積常常會比全站儀提取的稍大，特別是個別樹木出現(xiàn)了10%以上的相對誤差。通過對誤差較大樹木的點(diǎn)云數(shù)據(jù)觀察和分析發(fā)現(xiàn)，產(chǎn)生該問題的原因在于利用全站儀觀測時，由于人為操作，如某一樹高處，遮擋嚴(yán)重觀測不到樹干特征點(diǎn)，則可變換位置進(jìn)行觀測，能夠比較好地提取干形，而地面激光雷達(dá)則只是固定的三站，掃描的點(diǎn)云數(shù)據(jù)無法完整包括不同高度處樹干的特征信息，特別是樹葉較為密集處提取直徑時，可能將部分樹葉和枝干的點(diǎn)云數(shù)據(jù)參與計(jì)算，進(jìn)而造成了點(diǎn)云數(shù)據(jù)提取的材積普遍大于全站儀提取的材積。

表1 地基激光雷達(dá)掃描提取的胸徑和材積Table 1 Data analysis of DBH and volume based on ground-based laser scanner

將樣木實(shí)測數(shù)據(jù)導(dǎo)入SPSS軟件中，建立兩種方法的一元線性模型，公式如下：

V全站=1.02V激光-0.012

(5)

計(jì)算該線性模型的相關(guān)系數(shù)R為0.998，決定系數(shù)為0.996，說明兩種方法具有顯著的相關(guān)性。該模型回歸系數(shù)的標(biāo)準(zhǔn)誤差為0.022，說明模型回歸系數(shù)較好，進(jìn)一步驗(yàn)證了兩種方法相關(guān)性極高。

由此，可以得出地基激光雷達(dá)掃描方法可以作為活立木材積提取方法。

4 結(jié)論與討論

1)本研究采用改進(jìn)的Hough變換方法提取樹干任意處直徑，避免了以前隨機(jī)Hough 變換的無效累積，提高了算法性能。實(shí)驗(yàn)采用直接測量的胸徑值作為參考，通過地基激光雷達(dá)量測的胸徑值平均相對誤差為1.81%；實(shí)驗(yàn)通過采用基于地基激光雷達(dá)提取材積與全站儀提取材積進(jìn)行比較，材積十分接近，平均相對誤差為5.86%，并且存在顯著相關(guān)性。有效驗(yàn)證了基于地基激光雷達(dá)提取活立木材積的可行性。

2)本研究以地基激光雷達(dá)這一快速、準(zhǔn)確獲取被測物體三維空間結(jié)構(gòu)信息的技術(shù)為基礎(chǔ)，實(shí)現(xiàn)了對樣木的三維重構(gòu)。研究選取了56棵山楊樹進(jìn)行地基激光雷達(dá)掃描，樣木之間遮擋不嚴(yán)重，枝葉密度不高，樹干點(diǎn)云容易提取，掃描效果較好。未來的研究中應(yīng)選擇多個樹種開展實(shí)驗(yàn)，驗(yàn)證該方法的可行性。

3)與全站儀獲取立木材積方法比較，地基激光雷達(dá)雖然同樣需要多次設(shè)站掃描樹木，后期還需對點(diǎn)云數(shù)據(jù)各種處理、研建模型等過程，精度相似，過程更麻煩。但應(yīng)用地基激光雷達(dá)技術(shù)的真正意義在于一次掃描可以提取包括活立木材積在內(nèi)的諸多森林調(diào)查因子，例如葉面積指數(shù)、冠層空隙度、枝葉密度等全站儀難以獲取的參數(shù)，而樹木健康狀況、樹種識別等需要紋理信息的參數(shù)全站儀更是無法做到，在這些方面地基激光雷達(dá)的應(yīng)用已經(jīng)得到了很好的實(shí)驗(yàn)驗(yàn)證。同時，地基激光雷達(dá)比全站儀操作簡單，測量時間短，全站儀通常測量1棵立木需要30 min以上，地基激光雷達(dá)單次掃描3～5 min，1棵樹10 min即可完成，效率有很大的提升。而內(nèi)業(yè)數(shù)據(jù)處理部分，雖然地基激光雷達(dá)相對麻煩，但是通過編程開發(fā)計(jì)算程序以后可以重復(fù)使用。因此，利用本研究提出的方法提取活立木材積，在森林資源監(jiān)測中地基激光雷達(dá)有廣闊的應(yīng)用前景。