基于RGB-D SLAM手機的森林樣地調(diào)查系統(tǒng)研究

范永祥 馮仲科 陳盼盼 高 祥 申朝永

(1.北京林業(yè)大學(xué)精準林業(yè)北京市重點實驗室， 北京 100083; 2.安徽農(nóng)業(yè)大學(xué)理學(xué)院， 合肥 230036)

0 引言

森林具有多種環(huán)境功能和社會經(jīng)濟功能，與人類生存、發(fā)展息息相關(guān)。在環(huán)境方面，森林具有固碳、水土保持、調(diào)節(jié)氣候、保持生物多樣性等多方面的作用；在社會經(jīng)濟方面，森林為人類提供了食物、木材、能源等[1-2]。

森林資源信息能夠為不同層面的決策者提供管理森林、計劃造林等依據(jù)。森林清查的主要目的是收集不同決策者所需要的森林屬性，通常調(diào)查屬性包括生物量、蓄積量、樹種等[3]。傳統(tǒng)的地面清查方法是在調(diào)查區(qū)域合理布設(shè)樣地，并對樣地的森林屬性進行清查，通過這些點估計來反演整個研究區(qū)域平均或總體的林業(yè)屬性[4]。通常樣地調(diào)查屬性包括樹種、胸徑、樹高，生物量、蓄積量等屬性一般通過構(gòu)建以樹種、胸徑、樹高為自變量的模型進行估計[5]。除這些屬性外，樣地樹木位置是解算林分空間結(jié)構(gòu)參數(shù)如角尺度、混交度、大小比數(shù)等因子所必需的屬性[6]。傳統(tǒng)樣地調(diào)查中，使用胸徑尺或卡尺測量胸徑，使用測高器測量樹高，使用羅盤、皮尺或全站儀測量樹木位置。這些方法不僅耗時耗力，而且無法克服觀測者主觀因素的干擾。Laser-relascope[7]、電子測樹槍[8]等設(shè)備提供了同時測量樹木胸徑、樹高、位置等一體化解決方案，但仍然無法克服安裝復(fù)雜、觀測難、主觀因素干擾等問題。

激光雷達(Light detection and ranging, LiDAR)技術(shù)的出現(xiàn)及計算機運算能力的提高為樣地清查提供了新的解決方案[9-11]。利用激光雷達對樣地進行掃描并獲取樣地三維采樣點云，然后從點云中客觀地提取樣地屬性。地基激光掃描(Terrestrial laser scanning，TLS)作為一種地基激光雷達，已經(jīng)被許多學(xué)者用于樣地清查，并通過算法對森林屬性進行提取和評估[5,12]。TLS通常以兩種方式被使用：單次掃描和合并多次掃描。單次掃描即在樣地中心建立站點并對樣地進行掃描[13]。該方式數(shù)據(jù)處理簡單，但在較密的森林中由于遮擋的原因經(jīng)常發(fā)生漏測。多掃描即在樣地中建立多個掃描站點，然后利用算法將多次的單掃描數(shù)據(jù)進行合并[14]。該方式雖然避免了遺漏問題，但合并單次掃描數(shù)據(jù)的難度較大。顯然，TLS不適合在較大的樣地中使用。移動激光掃描(Mobile laser scanning，MLS)的出現(xiàn)解決了TLS多次掃描數(shù)據(jù)難以合并的問題，從而可以在較大的樣地中進行森林屬性清查[15-17]。移動激光雷達通常為載有全球衛(wèi)星導(dǎo)航系統(tǒng)(Global navigation satellite system，GNSS)、慣性測量單元(Inertial measurement unit，IMU)、激光雷達組合的運動平臺。其中，GNSS+IMU的組合能夠在移動激光雷達系統(tǒng)運動過程中對LiDAR的位姿進行測量，利用所測得的位姿可以對各個時刻的掃描數(shù)據(jù)進行合并。但在林下沒有被GNSS覆蓋的區(qū)域，MLS系統(tǒng)很難構(gòu)建全局一致的點云。即時定位及構(gòu)圖(Simultaneous localization and mapping，SLAM) 技術(shù)的出現(xiàn)使得在林下無GNSS信號的區(qū)域定位成為可能。SLAM系統(tǒng)在運動過程中，利用單目相機、雙目相機或RGB-D相機對周圍環(huán)境進行觀測，從而獲取觀測序列，然后使用該視覺觀測序列對周圍環(huán)境建圖，并估計SLAM系統(tǒng)的位姿[18-19]。目前的研究中，將激光SLAM和雙目SLAM作為MLS系統(tǒng)的輔助相對定位系統(tǒng)，從而使MLS系統(tǒng)在無GNSS信號樣地中掃描樣地時仍能獲取全局一致的點云[20-22]，但MLS系統(tǒng)構(gòu)造復(fù)雜、質(zhì)量較大、價格昂貴，尤其對于輪式系統(tǒng)還需要較好的林下工作條件。

ToF (Time of flight) 相機[23-24]作為一種測量深度的消費級產(chǎn)品，具有體積小、價格便宜、功率低等優(yōu)勢，是激光雷達的一種替代品，且經(jīng)常作為RGB-D SLAM的輸入。隨著智能手機運算能力的增強及SLAM算法的改進，RGB-D SLAM系統(tǒng)可以在手機上運行，從而使手機具有實時定位并構(gòu)圖的功能[25-26]。TOMATK等[27]使用這種手機掃描樣地，并在線下從獲取的點云中提取了樣地樹木位置、胸徑。FAN等[28]使用這種手機實現(xiàn)了實時獲取樣地的單木屬性——位置、胸徑和樹高，但并未評估樣地屬性。

本文利用帶有RGB-D SLAM系統(tǒng)的手機構(gòu)建森林樣地調(diào)查系統(tǒng)，主要實現(xiàn)樣地構(gòu)建、每木檢尺及林分/樣地屬性計算，并使用增強現(xiàn)實技術(shù)進行結(jié)果實時顯示與交互。應(yīng)用該系統(tǒng)在18塊半徑為7.5 m的圓形樣地中進行測試，對估計結(jié)果進行評估。

1 理論與技術(shù)

1.1 SLAM理論

SLAM為以控制及觀測為輸入，在未知環(huán)境中未知位置構(gòu)建周圍三維地圖的同時實時估計當(dāng)前運動平臺位姿的過程。即在SLAM過程中，運動平臺位姿及所有路標點的位置在沒有任何先驗信息的條件下被實時估計。從概率分布的角度講[29-30]，控制輸入可以描述為運動模型P(xk|xk-1,uk)，其中，xk-1、xk為k-1、k時刻的位姿，uk為控制輸入。

觀測輸入可以被描述為觀測模型P(zk|xk,m)。其中，zk為k時刻的觀測，m為所有路標點。

SLAM求解過程通常被視為具有馬爾可夫性，即未來的位姿僅與前一時刻的位姿有關(guān)系而與過去的位姿沒有關(guān)系，從而運動模型與觀測模型可以通過貝葉斯過濾(Bayes filter)來求解當(dāng)前位姿的后驗分布，先驗位姿的估計被描述為預(yù)測更新

(1)

式中Z0:k-1——0到k-1時刻的觀測

U0:k——0到k時刻的控制輸入

x0——初始時刻的位姿

利用先驗分布及觀測模型求解后驗位姿估計被描述為觀測更新

(2)

貝葉斯過濾僅從概率角度解決了SLAM問題，在實際使用時需要給出運動模型及觀測模型的具體形式，目前主流的解決方案包括擴展卡爾曼濾波(Extended Kalman filter，EKF)及其變體[31]、非線性優(yōu)化(Sparse non-linear optimization method)[30]及粒子濾波(Particle filter)[31]。

1.2 RGB-D相機技術(shù)

RGB相機是一種用于獲取周圍環(huán)境紋理信息的工具，常用于SLAM系統(tǒng)的觀測輸入，但由于其僅具有測角而不具有深度測量的功能，所以無法為SLAM系統(tǒng)所估計的位姿及地圖提供尺度。深度相機因具有測距功能，彌補了RGB相機不能獲取深度的問題。RGB-D相機作為RGB相機與深度相機硬連接組成體，可以同時獲取周圍環(huán)境的紋理信息及深度信息。

目前,深度相機進行深度估計主要基于飛行時間(Time of flight)原理，這種相機被稱為ToF相機。飛行時間一般通過使用脈沖信號或調(diào)制連續(xù)波來獲取，但目前主要使用基于調(diào)制連續(xù)波的技術(shù)。這種技術(shù)通過使用發(fā)光二極管(Light-emitting diodes, LED)向目標發(fā)射近紅外光(Near-infrared light, NIR)，當(dāng)光被反射回來時利用傳感器接收并獲取當(dāng)前相位差及信號強度，用于計算測量目標深度及其可靠性。ToF相機測距原理基本與激光雷達一致，但激光雷達采用逐點掃描的方式對周圍目標深度進行采樣，即在測量深度時需要記錄角度才能獲取所采樣點的坐標；而ToF相機使用普通相機中心投影、通過像素坐標的方式測量角度。

將ToF相機與RGB相機硬連接后的組合體標定其內(nèi)方位元素及相對外方位元素后便可以成為真正意義的RGB-D相機，然后便可將深度相機獲取的像素深度映射到RGB相機像素對應(yīng)的深度

(3)

其中

ud=[udvd1]T

(4)

(5)

(6)

式中dd——深度相機像素點坐標ud對應(yīng)的深度

Kc、Kd——RGB相機與深度相機的內(nèi)方位元素矩陣

fx、fy——x、y方向以像素為單位的焦距

cx、cy——投影中心投影在像素平面坐標系的坐標

Rc_d、tc_d——深度相機坐標系相對于RGB相機坐標系的旋轉(zhuǎn)矩陣及平移向量

dc——RGB相機像素點坐標uc對應(yīng)的深度

ud、uc——深度相機像素點坐標、RGB相機像素點坐標

1.3 Google Tango技術(shù)

雖然近年來SLAM算法得到改進，芯片運算能力得到提高，要使移動設(shè)備能夠?qū)崟r估計位姿仍然具有難度。Google Tango通過使用特殊硬件使智能手機能夠運行SLAM系統(tǒng)。如圖1所示，Google Tango以運動跟蹤相機(RGB相機)及ToF相機(深度相機)作為觀測輸入，9軸的加速度傳感器、重力傳感器、羅盤傳感器估計控制輸入，并使用專用的低功率計算機視覺處理器加快數(shù)據(jù)處理過程。Google Tango手機使用基于外觀的實時構(gòu)圖(Real-time appearance-based mapping，RTAB-Map)[32-33]SLAM系統(tǒng)，如圖2所示，RTAB-Map系統(tǒng)分為前端和后端兩部分。前端為視覺里程計，該部分使用RGB-D及IMU獲取的數(shù)據(jù)逐幀估計手機的相對位姿。隨著估計過程的進行，估計誤差累積使位姿漂移逐漸加重，后端的目的正是為了減少這種誤差，后端由局部光束法平差(Local bundle adjustment)模塊與回環(huán)檢測模塊組成， 局部平差模塊用于對最近的關(guān)鍵幀位姿進行平差修正；回環(huán)檢測通過檢測訪問過的區(qū)域使用位姿圖優(yōu)化(Pose graph optimization)對各時刻的位姿進行糾正。

圖1 Google Tango手機(Lenovo Phab 2 Pro)的結(jié)構(gòu)Fig.1 Structure of Google Tango phone (Lenovo Phab 2 Pro)

圖2 RTAB-Map系統(tǒng)位姿估計、地圖構(gòu)建過程Fig.2 Process of estimating pose and mapping in RTAB-Map system

2 樣地調(diào)查系統(tǒng)設(shè)計與原理

該樣地調(diào)查系統(tǒng)主要依托于SLAM系統(tǒng)。如圖3所示，SLAM系統(tǒng)以RGB-D相機及IMU為輸入并實時估計獲取的實時RGB圖像、點云數(shù)據(jù)、位姿等信息；樣地調(diào)查系統(tǒng)使用位姿及深度信息來估計單木參數(shù)，RGB圖像及位姿則通過OpenGL構(gòu)建三維虛擬場景，通過將SLAM坐標系統(tǒng)與OpenGL坐標系統(tǒng)對齊從而通過手機屏幕為使用者提供增強現(xiàn)實(Augmented reality, AR)場景，并使用戶能通過屏幕與AR場景進行交互。

圖3 樣地調(diào)查系統(tǒng)基本結(jié)構(gòu)Fig.3 Structure of forest inventory system

樣地調(diào)查系統(tǒng)主要包含定義樣地坐標系、構(gòu)建全局一致的稀疏地圖、每木檢尺及參數(shù)計算4部分，如圖4所示。定義的樣地坐標系用于描述樣地中每一棵樹木的位置；構(gòu)建全局一致性稀疏地圖使每木檢尺時所獲取的手機位姿通過回環(huán)檢測減少漂移，進而約束樹木位置的估計誤差；每木檢尺過程觀測樣地中的所有樹木；參數(shù)計算過程計算樣地所代表區(qū)域的林分參數(shù)等。

圖4 樣地調(diào)查系統(tǒng)工作流程圖Fig.4 Workflow of forest inventory system

SLAM系統(tǒng)會定義一個以系統(tǒng)啟動時設(shè)備位置為原點、x軸正向水平向右、z軸正向垂直向上、y軸正向指向深度方向的右手坐標系，即初始化坐標系。SLAM系統(tǒng)的位姿將以該坐標系為參考坐標系。樣地坐標系是以樣地中心為原點、x軸水平向東、y軸正向指向正北、z軸垂直向上的右手坐標系,可以更好地描述樣地中樹木位置。因此，需要獲取兩個坐標系的相對關(guān)系以實現(xiàn)將以初始化坐標系為參考系的樹木位置轉(zhuǎn)換為以樣地坐標系為參考系。兩坐標系之間的相對關(guān)系由平移和旋轉(zhuǎn)共6個自由度來確定，僅需獲取2個點在不同坐標系中的坐標值便可以獲取兩坐標系的變換矩陣。本文通過采集樣地中心點xi_center以及正北方向一點xi_north來獲取變換矩陣,其旋轉(zhuǎn)矩陣Ri_p及平移矩陣ti_p分別為

(7)

ti_p=xi_center

(8)

(9)

為了自動獲取樣地中心及正北方向一點，本文設(shè)計了棋盤如圖5所示，用其中3×3圓形格網(wǎng)指示樣地中心、4×4棋盤格中心指示樣地正北方向點，然后通過使用OpenCV庫棋盤格檢測的方法檢測兩棋盤格并獲取圓形格網(wǎng)及棋盤格中心坐標，從而獲取xi_center及xi_north，并經(jīng)過轉(zhuǎn)換獲取樣地坐標系(圖6a)。

圖5 確定樣地坐標系所使用的棋盤格Fig.5 Chessboard used for building plot coordinate system

SLAM系統(tǒng)通過回環(huán)檢測及位姿圖優(yōu)化可以減少前端視覺里程計位姿估計產(chǎn)生的漂移。若在觀測前獲取了樣地全局一致的地圖，則在每木檢尺時手機位姿漂移將通過后端得到優(yōu)化，所以本文在每木檢尺之前進行了樣地掃描從而構(gòu)建樣地稀疏地圖的步驟。全局一致性地圖的準確性與掃描軌跡有著密切聯(lián)系，掃描過程適當(dāng)?shù)膽?yīng)用回環(huán)檢測則可獲取準確的地圖，本文采用螺旋形路徑進行掃描以盡可能多的構(gòu)建回環(huán)檢測(圖7)，圖中綠色三角形表示樣地中心，也是掃描起點；紅色方形為掃描終點，該點接近于樣地邊界。為了確保能夠獲取完整的樣地地圖且不至于浪費時間掃描樣地外的地面， 或者在每木檢尺時避免采集樣地外樹木的信息，確定樣地邊界極為重要。本文使用增強現(xiàn)實技術(shù)在OpenGL坐標系內(nèi)構(gòu)建了樣地邊界，當(dāng)觀測者接近樣地邊界時，手機屏幕將會顯示樣地邊界的位置(圖6b)。

圖6 樣地觀測過程中樣地調(diào)查系統(tǒng)狀態(tài)Fig.6 Different statuses of plot survey system during observation

每木檢尺過程需要對樣地中所有立木的胸徑、樹高和樹種等進行觀測。為確定胸徑所在位置，首先需要確定地徑處一點以確定胸高位置(圖6c、6d)，然后利用胸高處的點云過濾、擬合可獲取立木胸徑及位置信息[28]，經(jīng)過映射到OpenGL坐標系便可以將結(jié)果展示在屏幕上(圖6e)；在完成胸徑、位置估計后，若需要觀測樹高則移動到能夠觀測到樹梢的位置，在屏幕上點擊樹梢，便可計算樹高 (圖6f)。

圖7 樣地掃描路徑Fig.7 Scanning trajectory for building plot map

每木檢尺過程獲取了立木的胸徑、樹高、立木位置、樹種信息(圖6g、6h)，所有樣地、樣地所代表的林分等信息便可由此計算(圖6i)。林分平均胸徑為

(10)

式中di——樣地中第i棵樹的胸徑，cm

N——樣地中立木總數(shù)

林分平均樹高為

(11)

式中hi——樣地中第i棵樹的樹高,m

第i棵立木材積為

(12)

式中fε——實驗形數(shù)

每公頃蓄積量(單位：m3/hm2)為

(13)

式中S——樣地面積，m2

每公頃橫斷面積(單位：m2/hm2)為

(14)

林分株樹密度(單位：株/hm2)為

(15)

若樣地中立木的位置平均值μ及協(xié)方差Σ分別為

(16)

(17)

其中

xi=(xi,yi,zi)T

式中xi——樣地中第i棵樹在樣地坐標系中的位置坐標

則樣地近似平面的法向量為協(xié)方差矩陣Σ最小特征值所對應(yīng)的特征向量nmin_eigenvalue，則樣地的坡度及坡向分別為

(18)

(19)

其中

(20)

式中Iz——樣地坐標系中z軸方向的單位向量

Iy——樣地坐標系中y軸方向的單位向量

E——將nmin_eigenvalue投影到Oxy平面的投影矩陣

3 實驗區(qū)域與方法

3.1 實驗區(qū)域概況

實驗區(qū)域為北京市西山實驗林場，該林場位于北京市近郊小西山(39°58′N，116°11′E)。選擇18塊半徑為7.5 m的圓形樣地為研究對象。為了使實驗具有代表性，所選的樣地具有不同的坡向、坡度，林分具有不同密度、樹種。所選樣地地面具有較少灌木且容易到達。表1為不同樣地立木的基本屬性。

表1 樣地屬性Tab.1 Summary statistics of plot attributes

3.2 實驗方法

為了驗證設(shè)計的系統(tǒng)在樣地調(diào)查中的精度，使用胸徑尺測量胸徑作為胸徑參考值，使用全站儀(Leica Flexline TS06plus)測量樹高作為樹高參考值，使用全站儀測量立木位置作為立木位置的參考值。樣地所代表區(qū)域林分參數(shù)的參考值均基于胸徑、樹高及立木位置參考值進行計算。為了給出所設(shè)計系統(tǒng)的估計精度，本文使用偏差(BIAS)、均方根誤差(Root mean squared error, RMSE)、相對偏差(relative BIAS)及相對均方根誤差(relative RMSE)對各測量值進行評估，計算公式如下

(21)

(22)

(23)

(24)

式中xi——測量值

xir——測量值xi對應(yīng)的參考值

n——測量值總數(shù)

BIAS——估計值偏差

Brel——估計值相對偏差

RMSE——估計值均方根誤差

Rrel——估計值相對均方根誤差

圖8 胸徑估計值散點圖Fig.8 Scatter plot of estimated DBH values

4 實驗結(jié)果

所設(shè)計系統(tǒng)通過樣地調(diào)查估計林分因子及樣地因子的關(guān)鍵在于測量樣地中立木胸徑、樹高及位置的精度。如圖8所示，胸徑估計值是無偏的且隨著胸徑增加估計值的誤差離散度不變。由表2可以看出，胸徑估計值具有較小的偏差(0.32 cm)及較小的離散度(1.27 cm)。圖9說明了樹高估計值是無偏的(8.0 cm)，而估計值的離散度隨樹高增加而變大，但由表2可知總體仍然具有較小的離散度(95 cm)。圖10描述了樹木在水平面上的位置誤差，總體而言立木位置是無偏的(x方向偏差為3.0 cm；y方向為-1.7 cm)，且離散度較小(x方向為3.8 cm；y軸方向為3.0 cm)，由表2可知，高程方向的偏差(0.4 cm)及離散度(1.1 cm)相較于x、y方向較小。

表2 單木屬性精度Tab.2 Estimations of tree-based properties

圖9 樹高估計值散點圖Fig.9 Scatter plot of estimated tree height values

圖10 立木平面位置估計誤差Fig.10 Position errors of all trees in plots

表3為18塊樣地的林分參數(shù)，所有的林分參數(shù)均與參考值相接近。僅樣地3株樹密度與參考值不一致。如表4所示，總體而言所有林分估計值具有較小偏差及較小的離散度。表5比較了樣地坡度及坡向測量值與估計值，結(jié)果表明坡度的估計值與參考值接近；當(dāng)坡度較小時，由于SLAM系統(tǒng)位姿估計誤差導(dǎo)致坡向估計值與參考值差異較大?？傮w而言，坡度及坡向的測量值仍然可靠。

表3 林分屬性估計結(jié)果Tab.3 Estimation results of forest-based properties

表4 林分/樣地屬性估計值的精度Tab.4 Accuracies of forest-based and sample-based properties

5 結(jié)論

(1)基于RGB-D SLAM手持式設(shè)備構(gòu)建了森林樣地調(diào)查系統(tǒng)，該系統(tǒng)使用SLAM系統(tǒng)估計了設(shè)備的相對位姿，利用ToF相機獲取周圍三維采樣，依此設(shè)計了計算立木胸徑、樹高及位置的算法，從而形成了構(gòu)建樣地、每木檢尺及林分/樣地參數(shù)估計的調(diào)查系統(tǒng)?；贏ndroid平臺顯示系統(tǒng)設(shè)計了增強現(xiàn)實功能，用于實時顯示調(diào)查結(jié)果，為觀測者提供檢查及交互功能。

表5 樣地屬性估計結(jié)果Tab.5 Estimation results of plot-based properties (°)

(2)該系統(tǒng)在18塊半徑為7.5 m的圓形樣地中進行了測試。與參考值對比表明，該系統(tǒng)在林分及樣地屬性估計時均能獲取較精確的結(jié)果。