熊妮娜,王佳
(1.北京市城市管理研究院,北京 100028;2.北京林業(yè)大學精準林業(yè)北京市重點實驗室,北京 100083;3.北京林業(yè)大學測繪與3S技術中心,北京 100083)
傳統(tǒng)的立木材積主要依靠胸徑尺、測高器等工具對標準樣地進行調查,以標準木胸徑、樹高查詢材積表進行估測,該方法效率低、人為干擾因素多,受材積表精度影響較大[1]。目前全站儀作為立木材積精準量測的儀器已經開展了相關實驗和應用研究,該方法技術成熟,但同樣存在操作復雜、不能海量獲取坐標點等問題。而地基激光雷達技術作為現代測量的新興技術手段之一,近年來被廣泛地應用在單木計測和森林調查各個方面,激光雷達掃描可以數秒內一次性采集到上萬點,在獲取目標三維坐標位置方面具有先天優(yōu)勢。目前在林業(yè)中應用于森林生態(tài)系統(tǒng)、森林資源調查的各領域,包括樹種、樹高、胸徑、生物量等單木屬性信息的提取,葉面積指數的估算,冠層空隙度,冠層輻射,冠層結構,葉面積分布,樹干曲線,競爭指數等[2]。隨著地基激光雷達體積逐漸減小,更加便于攜帶,同時生產成本的降低帶來價格的降低,地基激光雷達將會更加普遍地應用于林業(yè)調查中。
目前國內外學者開展了利用地基激光雷達提取單木參數的一系列研究。Thies等[3]、Astrup等[4]提出了一種基于立木的地面激光掃描數據重建樹干三維表面的方法和算法。Tansey等[5]根據地面激光掃描儀數據估算了科西嘉島松林地的樹木材積。鄧向瑞等[6]利用地基激光雷達對活立木進行了精確掃描測量,利用軟件提取樹木的直徑、樹高、冠幅等測樹因子,利用材積模型計算立木材積,并與伐倒木的材積作對比。王劍等[7]通過激光點云數據構建一系列圓柱模型片段疊加而成樹干的三維模型,該方法對樹干的體積計算也提供了一種便利。蔣佳文等[8]采用地面激光雷達進行多站點掃描獲取立木的點云信息,提取有關點云分布的特征參數,建立基于特征參數的材積模型。參考國內外學者的研究,本研究基于地基激光雷達點云數據,改進Hough變換方法,通過提取不同高度處的樹干直徑,運用斷面積區(qū)分求積法計算樹干材積,以期為無損精準測木提供新的方法。
本研究采用FARO 120型地基激光雷達,掃描儀穩(wěn)定可靠,測量距離為0.6~120 m,分辨率為0.1 mm,數據獲取率≤50 800像素/s,25 m內誤差≤2 mm,水平和垂直視野范圍為360°和320°。
在樹干胸高處朝南或朝北方向貼標靶紙,選取樣木周圍無遮擋的地方均勻放置3個參考球,選擇能同時見到3個參考球的均勻分布的3個位置設站,理想的測站位置間隔角度為120°(以待掃描樹為參考),安裝地基激光雷達并連接好相關設備,對掃描范圍、掃描分辨率等參數進行設置后,進行樹木掃描。本研究設置掃描區(qū)域為水平方向360°,垂直方向155°,掃描分辨率為10 000點/周,每株樣木需從不同角度掃描至少3次,耗時約10 min/株。
本研究點云數據處理參考張浩[9]、張冬等[10]采用的枝葉分離方法,主要包括4個步驟,結果如圖1所示:
1)點云數據預處理。將地基激光雷達掃描數據導入計算機中,通過其配套軟件Faro Scene,對點云數據進行加載、配準、濾波、剔除、導出等內業(yè)處理。打開Faro Scene軟件,直接把測站數據加載到軟件中,第一次打開軟件時,掃描數據的Mio掃描點數值默認為62,這個數值是根據計算機的內存自動顯示的,本研究涉及的點云提取使用的計算機內存為2 G,經過在不同內存的計算機中進行點云數據提取實驗,結果顯示,默認的Mio掃描點的數值均能滿足精度要求。
2)樹木點云數據主方向計算。樹木點云數據主方向計算需要借助近鄰點坐標信息及近鄰點相應的法向量信息,原始的樹木點云只包含點云坐標信息,本研究借助于主成分分析進行點云法向量計算;法截線擬合法計算主方向對噪聲點具有一定的免疫力,是目前計算點云主方向較好的一種方法,計算出的枝干部分的點云主方向與枝干中軸線平行,冠層部分的點云主方向雜亂無章。
圖1 內業(yè)處理過程Fig. 1 The interior work process
3)樹木點云數據枝葉分割。首先根據點云數據各點及相應的主方向構建圓柱,圓柱半徑定為2倍采樣間距,圓柱高度定為10倍采樣間距。然后將圓柱包圍盒內的點投影到圓柱軸線上,根據圓柱內的點與圓柱軸線的夾角、投影點在軸線上的分布密度進行判斷,區(qū)分冠層葉子點和枝干點。研究中分割時閾值?1設為0.7,閾值?2設為30°,閾值?3設為0.8,其結果較好。
4)枝干骨架特征點提取。枝干部分依據點云測地圖劃分水平集,基于最小二乘法對水平集進行二維圓擬合,同時通過目視判斷對擬合結果進行修正。
本研究利用自動檢測任意處直徑的方法,首先截取任意處直徑的點云數據水平切片,并進行柵格化,然后再使用圓擬合結合Hough變換,這是一個特征提取技術在特定形狀不完善的情況下由一個決策程序找到對象的規(guī)則形狀的方法,最終將圓檢測出來,進而計算直徑。該方法分為兩個步驟:
圖2 胸高部分水平切片截取算法和胸徑量測算法流程Fig. 2 The flow chart of DBH of individual trees retrieval algorithm
1)任意處直徑部分水平切片截取方法及步驟
任意處直徑的點云數據水平切片是從構成單木的點云數據中選取高于任意處直徑0.1 m至低于直徑處0.1 m的點云數據。以Matlab實現算法,具體過程如圖2所示。單木樹干任意直徑處上下0.1 m的小切片被截取下來,將此部分點云數據投影到xy平面,以便進行下一步邊緣檢測,并確定檢測圓的半徑。
2)任意直徑處量測算法及步驟
在林業(yè)測量中,常遇到直徑切面為不規(guī)則圓的情況,故本研究在不規(guī)則的直徑切面點云投影數據上檢測最能滿足截面輪廓的標準圓,該檢測圓的直徑即視為樹干任意處直徑。圓的檢測和半徑的計算通過Hough變換方法,由于Hough變換檢測圓程序需要讀入圖片,所以首先對水平切片投影后的數據進行柵格化處理,結合格網大小,換算出真實圓的直徑[11-12]。變換圓的方程為x2+y2+2ax+2by+c=0,用兩個累加器A(a,b)和B((a2+b2-c)1/2)來進行統(tǒng)計可能的圓心和半徑。
本研究采用了一種改進的Hough變換算法,通過降低隨機采樣的點數,以及利用采樣點的特征信息來判斷是否進行累積,來避免大量的無目標的隨機采樣與無用累積[13-14],使算法性能大大提高。具體算法如圖2所示。經測試,將最小圓周上邊緣點的個數定義為2 000較為合適。由于原始的點云數據是5個站點配準疊加的數據,數據量較大,如果參數設置太小,會導致檢測圓的效果不佳,檢測出多個小圓。第一次隨機選取3個邊緣點,將這3個邊緣點代入變換圓公式x2+y2+2ax+2by+c=0,用兩個累加器A(a,b)和B[(a2+b2-c)1/2]來進行統(tǒng)計可能的圓心和半徑。然后清除屬于該檢測圓的所有邊緣點。由于這里視樹干直徑為一個標準圓,故所有小圓的圓心位置視為樹干直徑標準大圓的邊緣點,進行邊緣檢測,當某次循環(huán)結束但產生邊緣點的統(tǒng)計個數達到或超過閾值的構造圓時,則圓檢測結束。檢測結果不是單個圓,而是以多個小圓圓心為邊緣點組成多個半徑一致、圓心位置相近的圓。所得結果對于獲得樹干直徑無影響,若想獲得圓心位置,將得到的多個圓圓心位置均值可看作所求圓心半徑,即為單木的位置。最后結合柵格圖像格網大小,得到單木的真實任意樹干處直徑。
以實驗提取的一棵山楊為例,提取的樹干不同高度的點云切片見圖3a,選取胸徑處的點云切片數據為例進行后面的直徑提取見圖3b,選取位置水平切片投影后的單木直徑點云數據分布見圖3c,柵格化后的水平切片見圖3d,利用Hough變換檢測出的圓效果見圖3e,圓心坐標、半徑長度和真實單木直徑在Matlab軟件中可自動計算出來。對本研究的單木,其圓半徑為56,再結合預先設置的柵格格網大小0.1 cm,樹干直徑即為56×0.1×2=11.2 cm。
圖3 利用Hough變換檢測出圓的過程Fig. 3 Detect circle by using Hough transform
本研究采用中央斷面積區(qū)分求積法計算單木材積,將樹干按一定長度(本研究采用2 m)分段,量出每段中央直徑和最后不足一個區(qū)分段梢頭底端直徑,利用中央斷面近似求積式,求算各分段的材積(V1,V2,V3,…,Vn),并合計(V)。
V=V1+V2+V3+…+Vn+V′=g1l+g2l+g3l+
(1)
對于研究立木材積精度驗證最可靠的方法為:將點云數據提取的立木材積和樹木伐倒后計算的材積相對比,但由于受限于政策原因,無法對樹木進行此類破壞性試驗。因此,本研究利用全站儀活立木材積測量方法對實驗樹木提取立木材積作為材積近似真值,全站儀無損測量立木的原理如下,首先用胸徑尺測量立木樹干的地徑D0和胸徑D1處(距離地面1.3 m處)的胸高h1,然后用全站儀測量儀器中心到胸徑D1處的斜距S,再用全站儀將胸徑以上的樹干分成n段并測量相應位置的水平角αi和天頂距γi(下角標i為分段號),后運用三角高程原理計算樹高H,并按照圓臺累加法計算得到立木樹干的總材積V,具體方法如圖4所示。目前已有研究表明,全站儀立木材積提取方法獲取的材積與伐倒木真實材積誤差在5%以內[15-16],可以近似作為真值檢驗其他方法。
注: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全站為全站儀提取材積數據,m3;V激光為地基激光雷達提取材積數據,m3。
本研究實測共56棵山楊樹,采用現場直接量測的樹木胸徑做為參考值,將地基激光雷達掃描提取的胸徑與其進行對比,其結果見表1。其中絕對誤差最大的為0.023 3 m,最小值為0.008 0 m,沒有絕對誤差的有9棵,占總數的16.07%;絕對誤差為正值的共有35棵,占總數的62.50%;為負值的為12棵,占總數的21.43%。從相對誤差來看,最大的為8.82%,最小的為0.00%,平均相對誤差為1.81%,可以看出地基激光雷達掃描提取胸徑與現場直接量測的樹木胸徑十分接近,精度較高。
針對實測的56棵山楊樹,運用本研究地基激光雷達掃描提取的材積與全站儀提取材積進行對比,其結果見表1。其中,絕對誤差最大的為 0.070 6 m3,最小值為-0.110 2 m3,絕對誤差為正值的共有37棵,占總數的66.1%,為負值的為19棵,占總數的33.9%。從相對誤差來看,最大的為15.07%,最小的為0.06%,平均相對誤差為5.86%,可以看出地基激光雷達掃描提取材積與全站儀提取材積十分接近,誤差較小。同時分析實驗結果發(fā)現,地基激光雷達掃描提取的材積常常會比全站儀提取的稍大,特別是個別樹木出現了10%以上的相對誤差。通過對誤差較大樹木的點云數據觀察和分析發(fā)現,產生該問題的原因在于利用全站儀觀測時,由于人為操作,如某一樹高處,遮擋嚴重觀測不到樹干特征點,則可變換位置進行觀測,能夠比較好地提取干形,而地面激光雷達則只是固定的三站,掃描的點云數據無法完整包括不同高度處樹干的特征信息,特別是樹葉較為密集處提取直徑時,可能將部分樹葉和枝干的點云數據參與計算,進而造成了點云數據提取的材積普遍大于全站儀提取的材積。
表1 地基激光雷達掃描提取的胸徑和材積Table 1 Data analysis of DBH and volume based on ground-based laser scanner
將樣木實測數據導入SPSS軟件中,建立兩種方法的一元線性模型,公式如下:
V全站=1.02V激光-0.012
(5)
計算該線性模型的相關系數R為0.998,決定系數為0.996,說明兩種方法具有顯著的相關性。該模型回歸系數的標準誤差為0.022,說明模型回歸系數較好,進一步驗證了兩種方法相關性極高。
由此,可以得出地基激光雷達掃描方法可以作為活立木材積提取方法。
1)本研究采用改進的Hough變換方法提取樹干任意處直徑,避免了以前隨機Hough 變換的無效累積,提高了算法性能。實驗采用直接測量的胸徑值作為參考,通過地基激光雷達量測的胸徑值平均相對誤差為1.81%;實驗通過采用基于地基激光雷達提取材積與全站儀提取材積進行比較,材積十分接近,平均相對誤差為5.86%,并且存在顯著相關性。有效驗證了基于地基激光雷達提取活立木材積的可行性。
2)本研究以地基激光雷達這一快速、準確獲取被測物體三維空間結構信息的技術為基礎,實現了對樣木的三維重構。研究選取了56棵山楊樹進行地基激光雷達掃描,樣木之間遮擋不嚴重,枝葉密度不高,樹干點云容易提取,掃描效果較好。未來的研究中應選擇多個樹種開展實驗,驗證該方法的可行性。
3)與全站儀獲取立木材積方法比較,地基激光雷達雖然同樣需要多次設站掃描樹木,后期還需對點云數據各種處理、研建模型等過程,精度相似,過程更麻煩。但應用地基激光雷達技術的真正意義在于一次掃描可以提取包括活立木材積在內的諸多森林調查因子,例如葉面積指數、冠層空隙度、枝葉密度等全站儀難以獲取的參數,而樹木健康狀況、樹種識別等需要紋理信息的參數全站儀更是無法做到,在這些方面地基激光雷達的應用已經得到了很好的實驗驗證。同時,地基激光雷達比全站儀操作簡單,測量時間短,全站儀通常測量1棵立木需要30 min以上,地基激光雷達單次掃描3~5 min,1棵樹10 min即可完成,效率有很大的提升。而內業(yè)數據處理部分,雖然地基激光雷達相對麻煩,但是通過編程開發(fā)計算程序以后可以重復使用。因此,利用本研究提出的方法提取活立木材積,在森林資源監(jiān)測中地基激光雷達有廣闊的應用前景。