張 鶴,李東升,陳愛軍
(中國計量大學計量測試工程學院,浙江 杭州 310018)
近年來,在國家節能減排政策的大力推動下,我國碳交易市場得到了快速發展,對林業碳匯提出了更高的要求。傳統的砍伐式林木生長量計算方法已經無法滿足當前碳匯精度的需要。為此,研發一種樹木活體材積計量方法正勢在必行。本文應用三維激光掃描技術對樹干和樹枝進行模擬測量,使用的模擬樹干和樹枝為可精確測量體積的標準管體,這樣就給誤差評定帶來很好的方便性[1]。
由于現有的三維激光掃描儀并不是為專門應用于樹木測量設計,而且由于樹枝的相互遮擋,利用三維激光點云數據進行樹枝內部結構體積參數計算尚存在難點,因此如何獲得精確的立木材積量需要應用特殊的算法。目前可用于計算樹木體積的三維模型主要有兩種:1)幾何模型,主要計算樹冠外包體積,樹木的樹枝部分以及樹枝分叉內部的黏連體積無法消除,對于分枝特征表現不足,計算得到的結構參數也更加單一,只能計算樹木冠層形態,無法全面地獲得樹木樹枝的體積數值[2];2)體素模型,建立體素模型,將待測區域分割為體素,包含激光點的有效體素個數乘以體素體積即得總體積[3]。這種模型在樹枝葉茂密時易受內部點云缺失影響,模型易形成內部空洞,也無法全面地指示樹枝體積形態[4]。
基于此,本文提出一種基于Alpha-Shapes的分段改進算法模擬樹枝三維體積測量的新方法。其中Alpha-Shapes算法由Edelsbrunner首先提出[5-6],最初用于點集輪廓的相關構建。目前該算法應用于網格生成、醫學圖像分析和建筑物外輪廓的邊緣提取,但如何將Alpha-Shapes算法應用于提取樹干樹枝點的點云輪廓,進行樹枝三維模型構建并進行體積計算還缺少深入的系統研究[7]。
該方法采用標準塑料直管段、三通管段分叉和四通管段分叉模擬樹枝自然情況,使用三維激光掃描系統獲取三種塑料管段的三維空間點陣數據,建立模擬樹枝的三維模型,利用計算復雜不規則物體體積的Alpha-Shapes算法對三維激光掃描模型進行處理,分段計算出塑料管段體積并以及推導用于樹枝體積,有助于提高樹木體積的計算精度,為林業碳匯量的計算提供了數據參考[8-11]。
基于Alpha-Shapes的分段改進算法主要應用了Alpha-Shapes算法并適當進行了點云數據的分段檢測。Alpha-Shapes算法適用于從一堆無序的點集中提取物體的邊緣進行包絡并可以進一步進行體積計算。包絡示意圖如圖1所示[12]。假設點集中半徑為 α的圓由點集內任意兩點P1(x1,y1)、P2(x2,y2)唯一確定,若圓內無其他點,則P1、P2為邊界點,線段P1P2為邊界線段。如圖2所示,并可以得到過這兩點的圓的圓心P3,即求出與P1點、P2點距離為 α 的點P3(x3、y3)。
圖1 包絡示意圖
圖2 Alpha-Shapes模型示意圖
由已有距離交匯算法得:
式中:
獲得圓心后,通過其他點到圓心的距離是否小于α值,并以此判斷P1、P2是否為邊界點。若無任何點到圓心的距離小于α值,則P1、P2為邊界點,線段P1P2為邊界線段[13]。
三維Alpha-Shapes算法的公式原理與二維原理相同,通過三點作半徑為α的球來判斷邊界點,并在得到的邊界點處建立三角片面,重構出曲面。該算法可適用于多面體內外輪廓線提取和體積計算,通過對α值的設定,來控制輪廓線的精細程度[14-15]。對于凹多邊多面體,如分叉樹枝部分,需要考慮凹邊拐角處和內部情況,若α值設定不合理,則拐角分叉處會被鈍化引起嚴重變形,故α值應根據實際情況進行取值[16]。
使用Alpha-Shapes算法進行邊緣檢測以及體積計算的時候,會產生黏連現象,因此在Z軸上采用分段算法進行體積分割計算。點Pi(Xi,Yi,Zi)(i=1,2,···,n,n為Z軸坐標點的總數)坐標范圍為Zmin~Zmax,按照Z軸方向以及一定高度β值對數據分段,則任意一點P(Xi,Yi,Zi)在Z軸分段的序號為
式中:β——分段高度;
Ci——Pi點以Z軸為分段方向的序號。
在分段過程中在最高層處可能存在數據不足、數據缺失等問題,因此需要對最高層點云進行點云數據判斷并對于最高層點云數量不足的地方加以調整。公式可表示為
式中:δ——Z軸點經過高度為β值的分段后剩余不足為一段的點數;
ε——滿足一個分段需要補足的點數[17]。
在δ個點中進行坐標比較,選取Z坐標值最大值,也就是最高層處進行補點,補充ε個與最高層重合的點,這樣既不會改變掃描目標的建模形狀,也能夠補足點云數量,讓點云數據足以分段。點云數據補充完整后則可以進行體積分段計算,按分段的序號依次進行體積計算,最后將體積累加。
Alpha-Shapes分段改進算法體積計算過程時選擇合適的檢測半徑α值和分段高度β值是一個重要問題[18]。檢測半徑α值和分段高度β值的大小將直接影響點云數據建模精度和體積計算誤差。因此,在實際的應用中需要在點云數據建模精度和體積計算誤差之間尋找一個平衡點,使其達到最好的計算效果。
針對樹枝建模及體積模擬,因樹枝情況復雜,體積分散,不便于驗證算法誤差值。綜合考慮采用標準塑料直管段對樹干進行模擬,采用三通管段和四通管段對樹枝的分叉情況進行模擬。
三維掃描實驗如圖3所示,在進行三維點云掃描實驗時,對上述管段進行豎直放置,并調整三維激光掃描儀和標靶的位置,其中標靶作為點云數據拼接的基準,其放置位置尤為重要,不僅影響站點布設的數量,同時也影響拼接精度。標靶應為3個,布置在公共掃描區域內,呈銳角三角形放置。在標靶布設過程中,需要根據測量精度要求控制標靶與掃描儀之間的距離,本實驗設計標靶與掃描儀之間的距離為30~60 cm。同時避免將3個標靶放置在同一直線上。記錄相應的位置坐標以便于后續實驗數據處理。
圖3 三維掃描實驗圖
本文以瑞士徠卡HDS7000三維激光掃描儀作為數據獲取平臺,采用基于標靶采集的方法進行三維掃描實驗:首先設置掃描密度為高等分辨率得到點云數據,并利用Cyclone軟件完成點云數據的拼接,處理點云數據,刪除噪聲點,在點云中確定掃描目標的起止點高度并記錄,最后將處理好點云坐標數據文件以txt格式導出。
以直徑為50 mm的標準塑料直管段為例,樣本的真實照片和掃描處理后得出的管段點云數據如圖4所示。得到管段點云數據后則導入Matlab中間進行體積計算與驗證。
圖4 標準塑料直管段點云數據圖
本試驗針對50 mm直徑標準塑料直管段、標準三通管段、標準四通管段模擬樹枝分叉進行三維激光掃描實驗,并分別設置Alpha-Shapes算法半徑參數α以及分段高度包含的激光點數β進行體積計算,并得到體積值V,并對其中產生的單層誤差體積值V1′及其累計誤差和體積值誤差 ?V1進行分析。
首先分析在使用標準塑料三通管段和四通管段進行模擬樹枝分叉掃描試驗時出現的黏連問題和鏤空現象,由于標準塑料三通管段和四通管段的黏連問題和鏤空現象原理相同,故此處僅分析標準塑料三通管段。
對標準塑料三通管段和四通管段進行體積計算時,β2值不變,當α2值較大時,在標準塑料三通管段分叉處會產生黏連問題,造成三維模型體積值V2增大,如圖5(a)所示;當α2值較小時,因點云數據只存在于標準塑料三通管段分叉表面,掃描得到的三維模型孔隙隨之加大,產生鏤空的現象,造成V2減小,如圖5(b)所示。因此,當α2取值較大或較小時,最終會導致體積值誤差 ?V2增加,需選取合適的α2值來降低黏連問題或鏤空現象造成的 ?V2,如圖5(c)所示。
圖5 不同β2值與α2值的三通管段建模示意圖
α2值不變,β2值改變時,當β2值增大時,三維建模受到噪聲點影響隨之增加,且在分叉處會產生較大黏連問題,如圖5(d)所示;當β2值減小時,分段計算次數隨之增多,會導致單層誤差體積值V2′值疊加,誤差增大,且建模也會因分段次數過多而不清晰從而導致 ?V2增大,如圖5(e)所示。由此可知,當β2值取值較大或較小時,最終會導致 ?V2增加,因此需選用合適的β2值來減低噪聲點影響或建模不清晰造成的體積值誤差,如圖5(f)所示。
針對標準塑料直管段,進行三維掃描實驗,建立三維模型并得到標準塑料直管段三維點坐標數據,輸入Matlab中設置不同α1值以及β1值進行Alpha-Shapes算法的體積計算,得到數據如表1所示。
表1 標準直管段不同α1值和β1值以及對應體積值V1 10-3 m3
由表1可知,對標準塑料直道進行計算時,β1值不變,α1大于0.25則誤差過大不予討論。當α1值由0.25向0.1變化時,隨著α1值減小則單層體積間隙會更明顯,V1′相對會增加,因此 ?V1隨著 α1的減小而增加;當 α1值由0.1向0.05變化時,?V1變化幅度較小,此時Alpha-Shapes算法計算體積較穩定,受噪聲點影響較小,由此可知該范圍是α1最佳取值范圍;而當α1值由0.05向0.025變化時,?V1在不同β1值中呈現誤差都較大,誤差值在0.935~0.999之間,由于鏤空現象造成體積值V1減少,因此產生較大 ?V1。
α1值不變,β1值改變時,因無分叉影響,所以分段會造成較大誤差,所以 ?V1隨著β1值的增加而增加。
通過Alpha-Shapes算法計算α1值取不同值,β1值取不分段的整體的體積值Vz1,進行單層誤差體積值計算,公式可表示為:
式中:i——計算體積值V1時的分段次數;
V1′——單層誤差體積值。
計算后繪制不同α1值以及β1值情況下的V1′值折線圖如圖6所示。
圖6 標準直段各分段單層誤差體積值V1′折線圖
由圖6可知,β1值不變,α1值改變時,單層誤差體積值V1′隨 α1值的減小而減小。α1值由 0.25 向0.1 變化時,V1值變化趨勢較明顯;α1值由0.075向0.05變化時,V1值變化平穩,并能取得V1′最小值;α1由0.05向0.025變化時,V1′值接近1,此時體積值誤差較大,予以舍去。
α1值不變,β1值改變時,V1′隨 β1值的減小而減小。β1值在100~500間時,V1′值較小,但由于分段次數過多,會導致V1′累計誤差增大,?V1增大;β1值在500~700間時,V1值變化較平穩;β1值在700~1 100間時,V1值變化較大。
綜上所述,對標準塑料直管段模擬樹干進行體積計算時,α1值由0.075向0.05變化,β1值在500~700范圍內變化時,體積值V1計算較為理想,V1′變化較平穩,?V1較小。
標準塑料三通管段與四通管段模擬樹枝分叉實物圖如圖7所示,三通管段進行三維掃描實驗并經過體積計算后可得體積值V2,數據如表2所示。四通管段進行三維掃描實驗并經過體積計算后可得體積值V3,數據如表3所示。
圖7 標準三通管段與四通管段模擬樹枝分叉實物圖
表2 標準三通管段不同α2值和β2值以及對應體積值V2 10-3 m3
表3 標準四通管段模擬樹枝不同α3值和β3值以及對應體積值V3 10-3 m3
由附錄表2可知,在三通管段數據中,當β2值在100~500范圍內變化時,由于黏連問題使得V2增加,V2′的累計使得V2減小,此時受V2′累計誤差的影響較大,因此整體體積偏小。當α2值在0.25~0.1范圍內變化時,?V2較小是因為分叉處黏連問題增加的體積抵消了因V2′累計誤差而減小的體積,但黏連問題和V2′累積問題并沒有減小,因此該α2值范圍不能采用;當α2值由0.1向0.025變化時,由于鏤空現象造成 ?V2較大,不能采用。
當β2值在500~1 100范圍內變化時,此時受V2′累計誤差的影響較小,受黏連問題的影響較大,V2偏大,可以通過改變α2值進行修正。當α2值由0.25向0.075變化時,?V2較大,但隨著 α2值的減小而降低;當α2值由0.075向0.05變化時,三維建模黏連問題較小,?V2較小;當 α2值由0.05向0.025變化時,由于鏤空現象造成 ?V2較大,不能采用。
而由表3可知,在四通管段數據中,α3值與β3值的變化情況與三通管段數據變化情況相同,則不再贅述。
同理,通過Alpha-Shapes算法計算三通管段α2值取不同值,β2值取不分段的整體的體積值Vz2,進行V2′計算并繪制折線圖,如圖8所示。而四通管段的單層誤差體積值V3′如圖9所示。
圖8 標準塑料三通管段各個單層誤差值V2′折線圖
圖9 標準塑料四通管段各個分段單層誤差體積值V3′折線圖
由圖 8(a)可知,β2值不變,當 α2在 10~0.25 范圍內變化時,單層誤差體積值值V2′隨α2的減小而減小;當α2值由0.25向0.025變化時,V2′趨于穩定,因此對該范圍著重分析。如圖8(b),當α2由0.25向 0.075變化時,V2′持續下降,當 α2由 0.075向0.05變化時,V2′趨于穩定。當α2由0.05向0.025變化時,V2′陡然減小。
α2值不變,V2′隨 β2值減小而減小,綜合對比,當 β2值700~900范圍內變化時,V2′值在不同 α2取值值中都比較穩定,?V2較小。
綜上所述,對標準塑料三通管段模擬樹枝分叉進行體積計算時,在分叉處選用合適的β值進行體積計算能有效減小體積值誤差,因此,β2值在700~900范圍內變化,α2值由0.075向0.05變化時,V2計算較為理想,V2′累計誤差影響較小,?V2較小。
同理可得,對標準塑料四通管段模擬樹枝分叉進行體積計算時,在分叉處選用合適的β3值進行體積計算能有效減小體積值誤差 ?V3。因此,β3值在700~900范圍內變化,α3值由0.075向0.05變化時,V3計算較為理想,V3′累計誤差影響較小,?V3較小。
本次實驗對真實樹枝進行掃描,實物圖和掃描圖如圖10所示,通過三維掃描實驗并進行體積計算后可得體積值V4,數據如表4所示。
圖10 真實樹枝實物與三維掃描圖
表4 真實樹枝不同α4值和β4值以及對應體積值V4 10-3 m3
由表4可知,V4的變化趨勢與標準塑料三通管段和四通管段的體積值變化相似。當β4值在100~500范圍內變化時,體積值誤差 ?V4較大,而且綜合考慮鏤空現象以及黏連現象,則該范圍不能采用。而當β4值在500~1 100范圍內變化時,此時受單層誤差體積值V4′累計誤差的影響較小,當α4值由0.1向0.05變化時,三維建模黏連問題較小,體積值誤差 ?V4較小。
同理,通過Alpha-Shapes算法計算α4值取不同值,β4值取不分段的整體的體積值Vz4,進行單層誤差體積值V4′計算并繪制折線圖,如圖11所示。
圖11 真實樹枝各個分段單層誤差體積值折線圖
由圖 11可知,當α4由 0.25向0.075變化時,V4′持續下降,當 α4由0.075向0.05變化時,V4′趨勢較理想。當α4由0.05向0.025變化時,V4′陡然變化并減小。
α4值不變,V4′隨 β4值減小而減小,綜合對比,當 β4值700~900范圍內變化時,V4′值在不同 α3取值中均比較穩定,?V4較小。
綜上所述,對真實樹枝分叉進行體積計算時,體積計算情況、α4值與β4值范圍選取以及誤差分布情況與標準塑料三通管段和四通管段的情況相似,說明該算法在真實樹木掃描體積計算中也能夠參考應用。
針對不同α值和β值情況下,標準塑料直管段、標準塑料三通管段和四通管段以及真實樹枝的體積值誤差如圖12所示。
圖12 各個塑料管段以及真實樹枝誤差體積值?V折線圖
由圖12(a)可知,在標準塑料直管段的體積計算中當 β1取值在 100~1 100 間,隨著 α1從 0.025~0.25間取值的增加而 ?V1增加,當α1值由0.075向0.05變化時,?V1變化幅度較小。而當 β1取值在700~900間時,?V1變化幅度最小,而且 ?V1也較小,因此此時為7.46%最佳值。
由圖12(b) (c)(d)可知標準塑料三通管段和四通管段、真實樹枝的體積計算中,當β值在100~500范圍內,α值由0.25向0.1變化時,分叉處黏連問題增加的體積抵消了因V′累計誤差而減小的體積,使得體積值誤差 ?V較小,但黏連問題和V′累積問題并沒有減小因此該范圍不能采用;當α值由0.1向0.025變化時,由于鏤空現象造成?V較大,不能采用。當β值在500~1 100范圍內變化時,當α2值由0.25向0.075變化時,?V較大,但隨著α值的減小而降低;當α值由0.075向0.05變化時,黏連問題較小,?V較小;當α值由0.05向0.025變化時,由于鏤空現象造成 ?V較大,不能采用。而當β值在700~900范圍內變化,α值由0.075向0.05變化時,?V取得最小值,范圍為4%~5%。
除了黏連現象和鏤空現象等建模算法造成的誤差外,在三維掃描時產生的噪聲點在體積計算時也會造成體積值誤差,而且在標準塑料管道以及真實樹枝的邊緣部分會產生邊緣性誤差造成點云數據不準確,邊緣提前錯誤等問題。在后期工作中可以在誤差的減小方面進行進一步研究。
本文提出的Alpha-Shapes分段改進算法,通過標準塑料直管段模擬樹干實驗分析可知,?V1隨β1值增加而增加,α1值在 0.05~0.075 間,β1值在 700~900間能得到較好的V1;通過標準塑料三通管段和四通管段模擬樹枝分叉實驗分析可知,在分叉處選用合適的α2、β2和α3、β3值進行體積計算能得到較好的體積值V2和V3,同時降低黏連問題和鏤空現象的影響,減少V2′和V3′的累積誤差,有效減小?V2和 ?V3,由實驗數據可知,β2和 β3值在700~900間,α2和α3值在0.05~0.075間時能得到較好的V2和V3,得到的體積誤差值在4%~5%之間。通過對真實樹枝的實驗分析可知真實樹枝分叉部分附和標準塑料三通管段和四通管段選取實驗最佳參數。
利用Alpha-Shapes分段改進算法對標準塑料直管段、三通管段和四通管段以及真實樹枝進行三維重建時,在α值和β值變化過程中,構建的三維模型也在有規律地變化。可以看出,不同α值和β值構建的三維模型能夠反映樹枝樹干的空間伸展形態、曲面體積、分枝特性及分叉處規則性等結構參數,研究結果可作為樹枝樹干結構體積計算的研究基礎。