基于視覺里程計的森林樣地調(diào)查系統(tǒng)研究

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

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

0 引言

森林占全球陸地生態(tài)系統(tǒng)的30%，為全世界50%以上的動植物提供了棲息地，是生物遺傳多樣性最為豐富的生態(tài)系統(tǒng)[1-5]，在氣候調(diào)節(jié)、碳循環(huán)和水土保持等方面發(fā)揮著重要作用[6-7]。森林清查是自然資源監(jiān)測和管理的重要組成部分，是了解森林資源現(xiàn)狀及其動態(tài)變化的主要方式[8]。森林清查依不同層次的森林管理者及政府決策者的需求不同而不同，為森林蓄積量、生長量、生物量及碳儲量的分析和評估提供系統(tǒng)的森林資源信息[9]。

森林樣地調(diào)查是森林資源清查的重要方式之一,不僅可以直接用于評估森林資源信息，也是空基等數(shù)據(jù)構(gòu)建模型的參考數(shù)據(jù)源。森林樣地通常為森林中具有代表性的小面積區(qū)域，通過對樣地中的樹木進行每木檢尺獲取屬性信息并加以總結(jié)，便可獲取森林總體的評價屬性。森林樣地調(diào)查的單木屬性主要包括樹種、位置、胸徑、樹高等，而生物量、蓄積等屬性可通過所構(gòu)建模型進行估計[10-11]。顯然，樣地調(diào)查的精度及效率直接影響森林清查的可靠性及效率[12]。在傳統(tǒng)森林樣地調(diào)查中，單木屬性信息的獲取主要依賴于簡單的測量工具，如胸徑尺、測高器、測樹槍等。這些方法耗時、費力，且存在較大主觀因素[13-15]。

隨著遙感技術(shù)的發(fā)展，利用傳感器獲取樣地點云，然后從該點云中提取立木參數(shù)的方式被許多研究者關(guān)注[16]。其中，主要使用的傳感器包括激光雷達(Light detection and ranging， LiDAR)和攝影測量(Photogrammetry)相機。激光雷達通過主動測量目標點與激光雷達距離的方式獲取深度信息[17]。由于遮擋等原因?qū)е略摲椒y得的樹木點云密度過小，甚至無法進行擬合以獲取樹木胸徑及位姿[18]。攝影測量是一種被動對周圍紋理進行觀測的工具，通過匹配、位姿估計、稠密化點云等步驟獲取森林三維點云。但由于森林結(jié)構(gòu)復雜，不一定能夠成功獲取同名點，而獲取稠密的樹木點云，當然也無法獲取較高精度的單木屬性[19-22]。

利用攝影測量很難保證所獲稠密點云合格并可以提取所有立木的單木屬性[23-28]。視覺里程計是一種僅使用相機便可以基于攝影測量理論實時估計相機位姿的一種技術(shù)。顯然，這種不需要全球?qū)Ш叫l(wèi)星系統(tǒng)(Global navigation satellite system，GNSS)信號覆蓋條件的技術(shù)適合在林下工作。本文設計一種以視覺里程計估計森林中獲取的連續(xù)圖像序列的位姿、利用非點云的方式獲取立木胸徑及位置的方法。利用該方法構(gòu)建圖像序列處理系統(tǒng)，并將其用于樣地中獲取圖像序列處理，以獲取樣地中立木胸徑及位置。

1 視覺里程計理論

視覺里程計(Visual odometry，VO)為以搭載在運動平臺上的視覺傳感器作為輸入，利用傳感器獲取的時間序列幀數(shù)據(jù)，進行實時估計搭載平臺位姿的技術(shù)[29]。VO系統(tǒng)所搭載的視覺傳感器主要包括單目相機(Monocular camera)、雙目相機(Stereo camera)和深度相機(RGB-D camera)3類[30]。視覺里程計工作所需的傳感器較為簡單，而且在不需要GNSS覆蓋的條件下便可進行相對位姿估計，可以用于機器人室內(nèi)相對導航及宇宙探索中，也適合于林下無GNSS信號條件下估計觀測設備相對位姿的情況[31]。

VO算法主要需要解決的問題是利用不同視覺傳感器獲取的具有時間相關(guān)性的圖像序列進行實時估計搭載平臺的位姿。若設視覺傳感器在N個不同時刻獲取的幀序列為I0:N={I0,I1，…,In}(圖1)，視覺里程計通過相鄰幀對同一環(huán)境進行重疊觀測的特性從而估計兩相鄰幀k-1與k之間位姿的相對變換，該剛性變換可以表示為

(1)

式中Rk-1,k——從k-1時刻到k時刻的旋轉(zhuǎn)矩陣

tk-1,k——從k-1時刻到k時刻的平移向量

所有連續(xù)幀的相對變換集合為T1:n={T0,1,T1,2,…,Tn-1,n}，若設不同時刻的位姿序列為C0:n={C0,C1,…,Cn}，則k時刻的位姿為

Ck=Ck-1Tk-1,k

(2)

且通常設k=0時刻的位姿C0作旋轉(zhuǎn)矩陣為單位陣、平移向量為零向量的假設。基于此便可以估計不同時刻的位姿。

圖1 視覺里程計的本質(zhì)Fig.1 Essence of visual odometer

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

2 森林樣地數(shù)據(jù)處理系統(tǒng)設計

森林樣地數(shù)據(jù)處理系統(tǒng)主要包括視覺里程計及森林調(diào)查系統(tǒng)兩部分(圖2)。視覺里程計部分以外業(yè)所采集的圖像序列為輸入估計各幀圖像的位姿；森林調(diào)查系統(tǒng)需要手動標記樹木位置及胸徑，然后基于同名點自動匹配構(gòu)建樣地坐標系并估計立木位置和胸徑。

2.1 視覺里程計

通常使用的視覺里程計為基于特征的里程計，其利用不同幀中出現(xiàn)多個相同特征構(gòu)建不同幀之間的位姿轉(zhuǎn)換關(guān)系。如圖3所示，為本文所構(gòu)建的基于特征點的視覺里程計估計位姿的基本流程圖，其利用有序的圖像序列，通過特征提取、特征匹配、運動估計和局部光束法平差4個步驟為獲取當前幀的位姿：① 特征點提取部分用于提取圖像中對應于地面的點，好的特征點應具有位置精確(位置及尺度)、可重復性、計算高效及強魯棒性等特征，本文選擇提取速度較快的FAST(Feature from accelerated segment test)方法，提取角點類型的特征點。②特征匹配部分需要建立不同幀中特征點之間的對應關(guān)系，并使用BRIF(Binary robust independent elementary feature)算子描述該特征點的特點且該算子之間的距離可以描述兩個特征點的相似性，便可依次獲取不同圖像、相同場景特征點的對應關(guān)系。③當獲取到匹配點時，即可利用P2P或EPNP估計當前幀的位姿。④利用局部光束法平差從統(tǒng)計優(yōu)化的角度精確化所估計的位姿。

圖3 視覺里程計工作原理Fig.3 Operating principle of visual odometer

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

圖4 基于全局地圖的視覺里程計Fig.4 Visual odometer based on global map

為了減少位姿估計過程中的漂移，本文視覺里程計采用全局稀疏特征點地圖的方式進行研究(圖4)。特征點地圖中特征包含描述子、位置及其標準差信息。當獲取到第1幀圖像的特征點時，假設該幀所在位置為世界坐標原點且無旋轉(zhuǎn)，故直接將所有特征點添加到特征點地圖中并賦予適當?shù)奈恢眉皹藴什睢．敨@取第2幀數(shù)據(jù)時與地圖中的特征點進行匹配，基于匹配的點利用P2P算法、局部光束法平差估計該幀的位姿，獲取到位姿后基于三角測量重新估計特征點地圖中特征點的位置及其標準差；除此之外應將新獲取的特征點添加到特征點地圖。此后獲取的幀數(shù)據(jù)則使用EPNP及局部光束法平差估計該幀圖像的位姿，在更新特征點部分與第2幀相同。一個典型的位姿估計結(jié)果如圖5所示。

圖5 位姿估計結(jié)果示例Fig.5 Example of pose estimation results

2.2 立木位置與胸徑估計系統(tǒng)構(gòu)建

利用樣地中拍攝的有序圖像數(shù)據(jù)估計樣地中立木的位置及胸徑，通常所需要的立木位置使用以x軸水平向東、y軸正方向指向正北、z軸垂直向上的右手坐標系(樣地坐標系)進行描述。然而，視覺里程計的全局坐標系是以第1幀圖像的位姿構(gòu)建的初始化坐標系。另外，視覺里程計假設第2幀圖像與初始化坐標系原點的距離為單位距離。顯然視覺里程計初始化坐標系描述下的位姿尺度不正確，且與目標坐標系——樣地坐標系不一致。顯然，在估計樣地中立木位置前需要糾正其尺度，并構(gòu)建從初始化坐標系到樣地坐標系的轉(zhuǎn)換，從而使所有幀數(shù)據(jù)能夠在樣地坐標系中描述，從而提取基于樣地坐標系下的立木位置。然后，不同于傳統(tǒng)基于MVS(Muti-view stereo)構(gòu)建稠密點云并利用點云提取立木信息的方法，本文試圖利用在圖像上直接標注立木位置及胸徑的方法估計待估單木信息。

2.2.1樣地坐標系構(gòu)建

為構(gòu)建樣地坐標系，在掃描樣地前首先在樣地中放置如圖6所示的特殊棋盤格，并使棋盤格上對應于樣地中心的點與樣地中心重合，箭頭方向指示正北方向。在內(nèi)業(yè)處理時，按照軟件提示選擇相關(guān)幀棋盤格對應的樣地中心點及正北一點(圖7a)，然后使用極線搜索和塊匹配技術(shù)[32]搜索其他相關(guān)幀對應于這兩點的像素坐標，并通過三角測量的方法估計這兩點在初始化坐標系的坐標值，然后便可依據(jù)棋盤格對應這兩點的實際距離修正所有幀的實際位姿。尺度計算公式為

(3)

式中L——棋盤格上兩點的實際距離

xnorth、xcenter——棋盤格上正北一點及中心點在初始化坐標系中的坐標值

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

在修正尺度后便可獲取正確尺度的樣地中心及正北一點，為確定樣地坐標系與初始化坐標系之間的變換關(guān)系還需要獲取豎直方向的一個向量，本文通過假設樣地中立木樹干豎直向上估計該向量。如圖7b所示，通過構(gòu)建的森林樣地調(diào)查系統(tǒng)在圖像中選擇立木地徑及樹干上一點，并使用極線搜索、塊匹配技術(shù)搜索自動估計相應的在初始化坐標系中的坐標，則豎直方向的單位向量為

(4)

式中xbottom_i、xtop_i——第i棵被標記立木的地徑處一點及樹干上一點的坐標

N——被標記立木總數(shù)

由此，便可以估計樣地坐標系x軸單位向量在初始化坐標系中為

(5)

y軸的單位向量為

Iy=IzIx

(6)

由此可以獲取由樣地坐標系到初始化坐標系的旋轉(zhuǎn)矩陣為

Rinit_plot=[IxIyIz]

(7)

及平移向量

tinit_plot=xcenter

(8)

可以將所有圖像位姿以樣地坐標系進行表示。

2.2.2立木位置及胸徑估計

圖7 森林調(diào)查系統(tǒng)操作界面Fig.7 Forest survey system operation interfaces

(9)

式中xplot——胸高處一點在樣地坐標系中的坐標值，其可由地徑處點坐標z軸值加1.3 m獲取

Rcamplot、tcamplot——圖像由樣地坐標系變換到相機坐標系的旋轉(zhuǎn)矩陣及平移向量，即圖像在樣地坐標系中的位姿的逆變換

K——相機內(nèi)方位元素矩陣

zcam——xplot在相機坐標系下的z軸坐標值

在獲取到胸高處一點，基于極線搜索和塊匹配精確化胸徑處點，確保該點落在樹干上。然后，需要在森林調(diào)查系統(tǒng)軟件中點擊樹木胸高處左、右側(cè)邊緣以完成樹木位置及胸徑的估計(圖7d)。在點擊立木左右邊緣后，一個以該圖像的位姿原點為原點、z軸水平指向胸徑處點、y軸豎直向上的右手輔助相機坐標系被構(gòu)建以便于參數(shù)估計，從而求解問題可以被描述到oAC-xAC-zAC平面內(nèi)(圖8)。則立木的平面位置及胸徑計算式為

(10)

(11)

其中xC=[xCzC]TxL=[xLzL]T

xR=[xRzR]T

圖8 樹木位置及胸徑估計問題Fig.8 Tree position and breast diameter estimation problem

式中xC、xL、xR——樹木位置在oAC-xAC-zAC平面的坐標

d——胸徑

xL、xR——樹干邊緣處點的平面坐標

xcorr_bh——xL和xR的中心點坐標

其由胸徑處點xp=[xpzp]T及左右側(cè)邊緣的像素獲取得

(12)

式中Rac_c——從相機坐標系到輔助坐標系的旋轉(zhuǎn)矩陣

在獲取到立木位置的平面位置xC后設置其高度為胸徑處點的高，便可換算到樣地坐標系從而獲取到樣地坐標系描述下的立木位置。

3 實驗材料與方法

3.1 研究區(qū)域概況

以北京市西山實驗林場(39°58′N，116°11′E)為研究區(qū)域。該林場海拔為300～400 m，平均降水量為630 mm。研究區(qū)主要包括人工純林和針闊混交林。選擇了12塊半徑為7.5 m的圓形樣地為研究對象，被選擇樣地地面灌草較少，以盡量減少環(huán)境對相機視野的阻擋。樣地的基本情況如表1所示。

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

3.2 實驗方法

采用大疆靈眸Osmo型口袋云臺相機，其成像傳感器為1/2.3”CMOS，該相機具有質(zhì)量小、結(jié)構(gòu)緊湊、智能、便攜等優(yōu)勢，另配置有較小的三軸機械增穩(wěn)云臺，抖動抑制精度高達±0.005°。通過Matlab 2014a軟件，利用張正友標定法對相機進行標定，所獲取相機內(nèi)參及畸變系數(shù)如表2所示。

為了驗證森林樣地內(nèi)業(yè)處理系統(tǒng)所估計單木參數(shù)的精度，使用胸徑尺測量胸徑作為胸徑參考值，立木位置使用全站儀測量值聯(lián)合胸徑計算獲取的立木位置作為樹木位置參考值。本文使用偏差(BIAS)、均方根誤差(Root mean squared error, RMSE)、相對偏差(relBIAS)及相對均方根誤差(relRMSE)對各觀測值進行評估，其計算公式為

表2 相機內(nèi)參標定數(shù)據(jù)Tab.2 Experimental result for camera calibration of intrinsic parameters

(13)

(14)

(15)

(16)

式中xi——測量值

xir——相對于測量值xi的參考值

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

BIAS——估計值偏差

Brel——估計值相對偏差

RMSE——估計值均方根誤差

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

4 結(jié)果與分析

4.1 立木位置估計

立木位置估計值在x、y軸方向分別為0.04、-0.03 m，且最大偏差不超過0.10 m(表3)。雖然樣地6在x軸方向的RMSE達到0.40 m，但總體而言估計值的RMSE(0.21、0.17 m)較小。圖9為所有樣地中立木位置的誤差，由圖9可知，雖然單個樣地中似乎立木位置存在系統(tǒng)誤差，但總體而言其誤差趨于0 m，且最大誤差小于0.7 m。

4.2 胸徑估計

經(jīng)與參考值對比，胸徑估計值具有較小偏差(0.09 cm, 0.51%)及較小的RMSE(0.88 cm，5.03%)(表4)。其中，最大偏差及RMSE分別為-0.95 cm及1.16 cm。由圖10可見，胸徑估計值與參考值相接近，且無偏差過大點；且不同徑階估計值也無明顯偏差且具有相近的變異性(圖11)。

5 討論

首先利用視覺里程計基本原理恢復了在森林樣地中獲取圖像序列的位姿，然后利用構(gòu)建的森林樣地調(diào)查系統(tǒng)通過構(gòu)建樣地坐標系、三角測量等過程估計了樣地中立木胸徑及位置。測試結(jié)果表明：立木位置估計值在x軸和y軸方向的BIAS分別為0.04、-0.03 m，RMSE分別為0.21、0.17 m；立木胸徑估計值的BIAS及RMSE分別為0.09 cm(0.51%)和0.88 cm(5.03%)。

表3 立木位置估計值的精度Tab.3 Accuracies of tree position m

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

樣地號RMSE/cmrelRMSE/%BIAS/cmrelBIAS/%10.914.900.442.3320.673.130.080.3530.593.880.130.8840.472.500.160.8551.166.820.633.6761.086.680.875.3570.854.01-0.25-1.1880.895.190.362.0791.047.020.392.65100.844.890.241.42111.046.42-0.95-5.85120.844.92-0.71-4.17平均0.885.030.090.51

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

圖11 胸徑估計值箱線圖Fig.11 Errors of DBH observations for different DBH values

目前對立木位置的研究相對較多，但較少明確給出精度。TANG等[32]利用移動激光雷達對不同的樣地進行了掃描，其運動軌跡的RMSE為0.32、0.29 m，但未評估樹木的精度。FAN等[10]利用RGB-D SLAM手機對立木位置進行了實時估計，其x軸及y軸方向均方根誤差(RMSE)均為0.12 m。LIU等[33]設計了一種實時動態(tài)多功能CCD連續(xù)攝影的地面觀測儀器，通過點云提取立木位置，其偏差最大為0.19 m，最小為0.137 m，均方根誤差(RMSE)最大為0.201 m，最小為0.162 m。顯然，本文所構(gòu)建調(diào)查系統(tǒng)所獲取立木位置精度與以上方法所獲取結(jié)果相似。

目前對立木胸徑獲取方法的研究較多。LIANG等[34]通過地面激光雷達獲得點云數(shù)據(jù)并利用圓柱體擬合的方法來提取胸徑，其胸徑估計值的BIAS為-0.18～0.76 cm，RMSE為0.74～2.41 cm；范永祥等[35]基于RGB-D SLAM手機所構(gòu)建的森林樣地調(diào)查系統(tǒng)對立木胸徑進行了評估，結(jié)果表明：胸徑估計值的BIAS及RMSE分別為0.36 cm(1.93%)及0.69 cm(3.89%)；QIU等[36]利用連續(xù)攝影測量原理及圖像識別技術(shù)研發(fā)了PDA(Personal digital assistant)攝影測樹儀器，該設備估計胸徑最大偏差為2.1 cm，絕大多數(shù)樣地測量偏差小于1.5 cm。本文所構(gòu)建系統(tǒng)所獲取胸徑精度與這些設備也相似。

6 結(jié)論

(1)設計了以視覺里程計估計、從森林樣地中獲取的有序圖像序列的位姿；并構(gòu)建了森林樣地數(shù)據(jù)處理系統(tǒng)，以實現(xiàn)尺度修正、建立樣地坐標系、立木位置及胸徑標記、獲取樣地中立木位置及胸徑。

(2)選取12塊半徑為7.5 m的圓形樣地為研究對象，利用相機獲取樣地中連續(xù)圖像序列。經(jīng)所設計系統(tǒng)處理，所獲取胸徑及立木位置均具有較高精度。