范永祥 馮仲科 閆 飛 申朝永 關天蒴 蘇玨穎

(1.季華實驗室, 佛山 528000; 2.北京林業(yè)大學精準林業(yè)北京市重點實驗室, 北京 100083;3.貴州省第三測繪院, 貴陽 550004)

0 引言

可靠的森林資源信息是評估森林發(fā)展狀況及制定森林發(fā)展計劃的重要依據(jù)[1]。樣地調查作為一種重要的森林資源信息清查方式,利用抽樣技術在森林中選取適當數(shù)量、尺寸且具有代表性的樣地,通過對樣地每木檢尺、歸納、總結便可反演森林發(fā)展現(xiàn)狀及消長變化[2],是全國森林資源清查(簡稱一類調查)、森林資源規(guī)劃調查(簡稱二類調查)及作業(yè)設計調查(簡稱三類調查)的重要手段[3]。近年來,隨著遙感技術及計算機數(shù)據(jù)處理能力的發(fā)展,星載、機載等遙感手段被用于完成不同尺度森林調查及監(jiān)測,但海量森林遙感數(shù)據(jù)的處理仍依賴于可靠、詳細、完整的森林樣地調查數(shù)據(jù),以便對星載數(shù)據(jù)、機載數(shù)據(jù)及反演模型進行校準[4]。高效、精確地完成森林樣地調查對掌握森林資源現(xiàn)狀與發(fā)展動態(tài)至關重要。

森林樣地調查通過每木檢尺對樣地中所有立木屬性如樹種、胸徑、樹高、第一枝下高、冠幅及立木位置等測定,經(jīng)統(tǒng)計便可獲取平均及總體的樣地林分屬性;而其他的屬性(如葉面積指數(shù)、莖曲線、生物量及碳儲量等)可利用樹種、胸徑、枝下高、冠幅及樹高等建立反演模型進行估計[5-6]。每木檢尺的精度直接決定了樣地調查的精度及各反演模型的精度。傳統(tǒng)樣地調查中,使用一些簡易的工具進行每木檢尺等工作。其中,胸徑采用胸徑尺、輪尺等測量,樹高采用布魯萊斯測高器、超聲波測高器等測定[6]。近年來,一些簡單集成測距及測角傳感器的新型設備被研發(fā)并用于森林樣地清查中[7-11],這些設備雖然解決了每木檢尺需要多人、多種設備配合測量的問題,但仍然是一類以人為主要觀測者的設備,故仍然是一類費時、費力、主觀性強的測量設備。但由于人工對測量結果的監(jiān)督作用,胸徑等測量值仍被認為是可靠的,經(jīng)常作為評估其他調查方法的參考數(shù)據(jù)[5]。

隨著遙感技術及計算機運算能力的發(fā)展,以面陣相機、ToF(Time of flight)相機及激光雷達(Light detection and ranging, LiDAR)等為主要觀測傳感器的近景遙感手段在森林樣地調查中得到廣泛應用[4,12-13]。該類型設備按采集方式可分為靜態(tài)采集設備及動態(tài)采集設備。典型的靜態(tài)采集設備如地面激光雷達(Terrestrial laser scanning, TLS)等,該類型設備需安置于樣地中進行靜態(tài)作業(yè),根據(jù)安置站數(shù)可分為單次掃描及多次掃描兩種作業(yè)模式,其中單次掃描模式容易受遮擋影響而無法獲取樣地完整的點云信息[14-15],多掃描模式則面臨多站數(shù)據(jù)拼接難、作業(yè)效率低、不適合在較大樣地中作業(yè)等問題[16-17]。動態(tài)采集設備將視覺傳感器數(shù)據(jù)與IMU(Inertial measurement unit)及GNSS(Global navigation satellite system)等傳感器數(shù)據(jù)融合,使設備可以在移動中進行數(shù)據(jù)采集并獲取樣地三維數(shù)據(jù)[18-20];特別是,SLAM(Simultaneous localization and mapping)技術的出現(xiàn),使移動測量系統(tǒng)無需依賴于GNSS信號實現(xiàn)了高郁閉度林下定位及建圖[21-23]。相比于傳統(tǒng)設備,該類設備通過在樣地進行數(shù)據(jù)采集后利用計算機對數(shù)據(jù)進行處理,實現(xiàn)每木檢尺及樣地屬性提取,減少了調查人員的工作量且避免了主觀觀測誤差[4]。然而,專業(yè)遙感設備仍面臨價格昂貴、質量較大及無人工監(jiān)督導致魯棒性差等問題。

智能手機作為一種廉價、便攜及高普及性的設備,是一種理想的森林樣地調查設備。近年來,手機芯片、傳感器等技術得到長足進步,使消費級面陣相機、ToF相機、激光雷達及IMU等傳感器均可被嵌入智能手機設備上,Tango、AR Core、AR Kite等SLAM算法可在手機上運行并支持增強現(xiàn)實的實現(xiàn),這些均為智能手機在森林樣地調查中提供了新的機會[24-26]。文獻[27]利用內嵌有ToF相機的智能手機掃描并獲取了森林樣地三維點云,獲取立木胸徑及位置均方根誤差分別為1.15～1.91 cm及0.200～1.135 m。文獻[28]利用內嵌有LiDAR的Apple iPad掃描并獲取了樣地三維點云,并提取到立木胸徑估計值均方根誤差為3.10～3.40 cm及0.109～0.218 m。文獻[29]利用內嵌有LiDAR的iPhone及iPad開發(fā)了ForestScanner,用于實時估計立木胸徑及位置,胸徑及立木間相對距離估計值均方根誤差為2.3 cm及3.1～4.4 m。課題組利用內嵌有ToF相機的智能手機研發(fā)了增強現(xiàn)實森林樣地調查系統(tǒng),該系統(tǒng)利用RGB-D SLAM系統(tǒng)實時數(shù)據(jù)實現(xiàn)胸徑、立木位置及樹高測量,利用基于樹干回環(huán)檢測提高了立木位置精度,利用增強現(xiàn)實場景實時人工監(jiān)督測量結果并交互修正誤差,胸徑及立木位置的估計值均方根誤差為1.00 cm及0.078～0.086 m[30-32]。然而,目前主流智能手機并非內嵌有ToF相機及LiDAR等傳感器,這將影響智能手機在森林樣地調查中的推廣,故亟需研發(fā)基于單目SLAM系統(tǒng)(僅使用面陣相機及IMU作為數(shù)據(jù)源)的森林樣地調查系統(tǒng);此外,該系統(tǒng)在胸徑估計中采用單幀深度圖進行圓擬合估計,不僅魯棒性差且不適合對傾斜立木進行精確胸徑與立木位置估計。本文利用帶有單目SLAM系統(tǒng)的手機構建森林樣地調查系統(tǒng);構建基于單目SLAM數(shù)據(jù)過濾胸高點云數(shù)據(jù)的方法;基于點到圓柱體表面距離及圓柱體切線到圓柱體表面距離構建胸徑與立木位置精確估計算法;利用增強現(xiàn)實技術進行結果實時顯示與操作交互。然后,該系統(tǒng)在5塊32 m×32 m的方形樣地中進行測試,并對立木位置及胸徑估計結果進行評估。

1 單目SLAM技術及深度圖估計方法

1.1 單目SLAM技術

本文單目SLAM系統(tǒng)以面陣相機為觀測傳感器、IMU為運動傳感器,其中:面陣相機的優(yōu)勢在于長時間、長距離位姿估計誤差較小,IMU的優(yōu)勢在于短時間、短距離內的位姿估計較準(即使在運動快速變化時),兩種傳感器數(shù)據(jù)融合可有效彌補各自的不足,是目前較為主流的SLAM系統(tǒng)之一。若設0～K時刻面陣相機對路標點{m}的觀測序列為

Z0:K={Z0,Z1,…,ZK}

IMU運動估計序列為

u1:K={u1,u2,…,uK}

則SLAM過程將相機觀測數(shù)據(jù)與IMU數(shù)據(jù)進行融合,并對各時刻位姿

x1:K={x1,x2,…,xK}

及路標點{m}進行實時估計(圖1);從概率論角度而言,SLAM過程表達為后驗概率推斷

圖1 SLAM過程貝葉斯網(wǎng)Fig.1 Bayesian network of SLAM process

(1)

SLAM解決方案正是通過對該后驗概率推斷采用不同的簡化假設(如馬爾可夫性假設等)并利用各種濾波及滑窗優(yōu)化技術實現(xiàn)對位姿及路標點的后驗推斷[33-34]。

1.2 深度圖估計方法及立木樹干深度圖特點

深度圖的估計以SLAM過程獲取的位姿及稀疏地圖為輸入,在朗伯面假設、光度一致性假設及可視性約束下,通過構建光度一致性非線性優(yōu)化模型實現(xiàn)[35-36]。該非線性模型可簡單表示為

(2)

式中IR(i,j)、Ik(i,j)——當前Patch及第k幀相應Patch中(i,j)偏移像素點光度

ck——第k幀顏色尺度

E——當前幀附近k幀相應Patch光度殘差之和

通過對該光度殘差平方和極小化,完成對當前Patch中心像素的深度等進行優(yōu)化估計并通過區(qū)域生長獲取深度圖及其置信圖。通過該方法獲取立木樹干深度圖及由深度圖轉換的點云,如圖2所示。由于森林中灌木等干擾、深度圖估計中對樹干圓柱體相切位置點估計不精確等原因,深度圖中存在大量噪點,從而在利用該深度圖進行立木位置及胸徑等估計時造成影響,本文將構建相應的過濾算法及胸徑等估計算法解決該問題。

圖2 利用多視圖影像重建立木樹干深度圖及單幀點云Fig.2 Depth images and single-frame point cloud of tree trunk reconstructed by multi-view stereo technology

2 測樹系統(tǒng)架構與核心算法

2.1 測樹系統(tǒng)工作框架

本文測樹系統(tǒng)基本工作框架及界面參考課題組已構建的RGB-D SLAM測樹系統(tǒng)工作架構進行,該系統(tǒng)主要由單目SLAM系統(tǒng)、測樹核心系統(tǒng)及測樹用戶交互系統(tǒng)3個模塊組成,其中:SLAM系統(tǒng)用于位姿估計與稠密深度圖估計,測樹核心系統(tǒng)用于接收用戶指令并利用SLAM輸出數(shù)據(jù)完成每木檢尺等工作,測樹用戶交互模塊用于為用戶提供增強現(xiàn)實窗口并通過手勢等完成與測樹核心系統(tǒng)交互。本測樹系統(tǒng)增強現(xiàn)實交互界面如圖3所示,每木檢尺操作工作流程如圖4所示。相比于RGB-D SLAM測樹系統(tǒng),新型測樹系統(tǒng)算法利用多次、多視角測量數(shù)據(jù)對立木位置及胸徑進行更精確估計,故在交互界面中增加“確認”按鈕以協(xié)助該功能的完成。

圖3 測樹系統(tǒng)增強現(xiàn)實交互界面Fig.3 Augmented reality interface of tree measurement system

圖4 森林樣地調查中測樹系統(tǒng)手勢操作及計算流程圖Fig.4 Gesture operation and calculation flowchart of tree measuring system in plot survey

2.2 立木位置與胸徑估計算法

在進行立木位置與胸徑估計前,本文以被點擊像素點坐標以及該時刻深度圖為輸入,先對深度圖中噪點進行過濾,然后對深度圖中立木圓柱體切點進行分類與精確估計,具體流程為:根據(jù)測量時設備與立木距離經(jīng)驗值,移除深度圖中深度小于0.2 m或大于1.3 m的估計值;根據(jù)像素周圍點深度的平滑度過濾樹干周圍灌木等遮擋噪點及樹干圓柱體相切位置不合理點,平滑度計算公式為

(3)

其中,d為當前像素的深度,dij為x、y方向偏移量分別為i、j像素的深度,本文中S>0.008時深度估計被認為無效;查詢距離被點擊像素最近的、深度有效的像素點標記為種子點,查詢該種子點周圍8個像素中深度有效的像素并標記為新的種子點,然后繼續(xù)索引新種子點周圍有效深度像素及種子標記,直至無新種子被添加為止,最后將所有種子作為樹干上的點;將深度圖左右兩側邊緣點作為切點候選點,利用相機位姿及內參將候選點變換到世界坐標系后,將各點按照極坐標角度排序,將中位數(shù)點作為提取到的樹干圓柱體切線點(圖5)。將深度圖中各點及切點轉換到局部坐標系后,該掃描預處理過程獲取數(shù)據(jù)最終可表示為{T,{x},xL,xR,θ},其中:T表征該時刻相機的位置t及姿態(tài)R,{x}為位于樹干圓柱體表面的估計點集合,xL、xR為樹干左右兩側切點,θ為觀測點對樹干圓柱體的觀測角,公式為

圖5 樹干圓柱體深度圖噪聲過濾結果Fig.5 Noise filtering result of tree trunk cylinder depth map

(4)

(5)

其中

式中θk——第k次掃描觀測角

nk——第k次掃描中樹干圓柱體總點數(shù)

x=tk+τnk,i

(6)

其中

式中τ——自變量

x——切線上任意一點

tk——第k次掃描時相機所在位置

nk,i——切線方向向量

由該切線構建的損失函數(shù)為

(7)

其中

式中jDk——立木胸徑

利用所構建損失函數(shù)進行優(yōu)化后,便可獲取到樹干圓柱體參數(shù),然后將估計結果構建到增強現(xiàn)實場景中進行監(jiān)督即可。

3 研究材料與方法

3.1 研究區(qū)域概況

選擇位于北京市近郊小西山的西山試驗林場(39°58′N,116°11′E)為研究區(qū)域,該林場主要以人工林為主體。本文選擇其中4塊32 m×32 m的方形樣地為研究對象。所選樣地地面具有較少灌木且容易到達,不同徑階組比較均衡。表1總結了不同樣地的基本屬性。

表1 樣地屬性描述統(tǒng)計Tab.1 Summary statistics of plot attributes

3.2 數(shù)據(jù)采集

為充分利用回環(huán)檢測修正立木位置,本文數(shù)據(jù)采集中使用如圖6所示測量軌跡(圖中棕色圓圈為立木,數(shù)字表示觀測順序,綠色箭頭表示行徑方向,紅色箭頭表示閉環(huán)觀測),即從樣地中心出發(fā)并測量距離樣地中心最近的立木;行經(jīng)到西北角并對行徑路線上鄰近立木進行每木檢尺;按仿航線法逐行完成每木檢尺;回到樣地中心,并對行徑路線上立木進行二次觀測。為充分研究新型設備胸徑觀測次數(shù)對測量精度的影響,本文試驗中選擇單次觀測、正交觀測、對稱觀測及環(huán)繞觀測4種觀測方式分別完成同一樣地的每木檢尺,4種觀測方法對胸高圓柱體觀測視角如圖7所示(圖中棕色圓圈為胸高圓柱體,紅色方框為設備觀測位置,黑色虛線為設備觀測方向)。

圖6 每木檢尺中行徑軌跡Fig.6 Path track in forest plot survey

圖7 不同觀測方式中觀測視點位置與方向Fig.7 Viewpoint position and direction by different observation methods

3.3 研究方法

為驗證本文中新型設備在森林樣地調查中的精度,使用胸徑尺測量值作為胸徑參考值、使用全站儀(Leica flexline TS06plus)對立木位置進行測量并作為參考值。使用偏差BIAS、均方根誤差(Root mean squared error,RMSE)、相對偏差Brel及相對均方根Rrel對各測量值的偏差及變異性進行評估。

4 試驗結果

4.1 位置估計精度

新型設備在5塊樣地中立木位置估計誤差統(tǒng)計結果如表2所示,結果顯示:X軸及Y軸方向的估計值偏差范圍為-0.014～0.020 m,顯然無偏;兩軸方向估計值均方根誤差小于0.09 m,顯然變異性較小;立木位置最大觀測誤差小于0.25 m,且魯棒性較好(圖8、9);相比于單次觀測,正交觀測、對稱觀測及環(huán)繞觀測方法中位置估計均方根誤差、最大觀測誤差等略有減小,其中環(huán)繞觀測均方根精度最高,但精度提升不顯著。顯然,從效率角度講,單次觀測方法最佳,且精度能滿足日常需求;從精度及效率角度而言,正交觀測及對稱觀測均可取得較高精度;此外,隨著觀測次數(shù)的增加,可一定程度上提高精度,但觀測次數(shù)超過2次(如環(huán)繞觀測)時,精度提升效果不顯著。

表2 立木位置估計值精度

圖8 立木位置誤差分布Fig.8 Position errors of all trunks in plots

圖9 立木位置估計誤差箱狀圖Fig.9 Errors of tree position observations

4.2 胸徑估計精度

使用不同觀測方法獲取5塊樣地中胸徑估計誤差統(tǒng)計結果,如表3所示。結果顯示:不同測量結果與參考值一致且具有較小偏差(-0.85～-0.03 cm,圖10),且均方根誤差較小(1.32～2.51 cm);在不同徑階組下,胸徑估計值的偏差、均方根表現(xiàn)較為一致(圖11);隨著觀測次數(shù)的增加,測量值的偏差逐漸減小,顯然胸高圓柱體數(shù)據(jù)的完善有助于對胸徑進行更高精度的估計(圖10、11);隨著觀測次數(shù)的增加,測量結果魯棒性提高(圖11),這是由于當樹干無法近似為圓柱體時,單次測量值很難獲取樹干胸高圓柱體完整信息,不同側面觀測數(shù)據(jù)加入優(yōu)化算法后有效抑制數(shù)據(jù)出現(xiàn)噪聲。

表3 立木胸徑估計值精度Tab.3 Accuracies of DBH estimations

圖10 胸徑估計值散點圖Fig.10 Scatter plot of DBH estimation

圖11 胸徑估計值誤差箱狀圖Fig.11 Errors of DBH observations under different DBH

5 結論

(1)以內嵌有普通相機及IMU的智能手機為硬件系統(tǒng)、單目SLAM系統(tǒng)為數(shù)據(jù)源,研究立木胸徑及位置的測量算法,以利用內嵌這些普通傳感器的手機完成森林樣地調查。首先基于平滑度設計了從多視圖幾何重構深度圖中過濾胸高圓柱體數(shù)據(jù)的方法,并獲取高魯棒性胸高圓柱體表面點云數(shù)據(jù)及切線;然后,利用多組胸高圓柱體表面點云及切線數(shù)據(jù)設計并構建了胸高圓柱體擬合損失函數(shù),實現(xiàn)對立木胸徑及位置精確擬合的算法;最后,以該算法為基礎在手機平臺上構建了交互式增強現(xiàn)實測樹系統(tǒng),并用于森林樣地調查。

(2)新型測樹系統(tǒng)在5塊32 m×32 m方形樣地中進行了測試,每塊樣地設計并使用了4種不同的觀測方法對立木胸高圓柱體觀測。經(jīng)與參考值進行對比表明,新型測樹系統(tǒng)可高魯棒性測量立木位置及胸徑,且通過兩次以上對胸高圓柱體觀測可提高胸徑估計精度(特別是不可近似為圓柱體的立木樹干)。新型測樹系統(tǒng)可完成森林樣地調查。