范永祥 馮仲科 申朝永 閆 飛 蘇玨穎 王 蔚

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

0 引言

森林資源及其監(jiān)測所獲取的可靠信息可以為不同層級、不同目標(biāo)的林業(yè)部門制定可持續(xù)發(fā)展政策、投資決策等提供支持，從而為森林培育、森林采伐、氣候變化影響評估、火災(zāi)模型和碳儲量估算等奠定基礎(chǔ)[1]。森林資源清查是獲取森林資源信息的重要途徑。由于不同層級的使用者對森林資源信息的使用目的不同，所關(guān)心的森林資源屬性也不盡相同[2]。為滿足清查目的，需要根據(jù)森林資源及所要求的精度制定不同的調(diào)查策略。樣地調(diào)查是一種重要的森林資源清查方法，其通過在研究區(qū)域合理布設(shè)樣地，然后對樣地的森林屬性進(jìn)行調(diào)查、統(tǒng)計，從而反演出整個區(qū)域總體或平均的森林屬性[3]。通常樣地直接調(diào)查的屬性包含樹種、胸徑、樹高及立木位置等，生物量、碳儲量等則可通過以上屬性構(gòu)建的反演模型進(jìn)行估計[4]。傳統(tǒng)樣地調(diào)查中，使用一些簡單的儀器進(jìn)行森林調(diào)查(如使用胸徑尺或卡尺完成胸徑測量，使用羅盤及皮尺等完成立木位置測量)，這些方法不僅耗時耗力，且無法克服觀測者主觀因素的干擾。隨著測距、測角等相關(guān)技術(shù)的發(fā)展，Laser-relascope[5]、電子測樹槍[6]等設(shè)備提供了同時測量樹木胸徑、樹高、位置等的一體化解決方案，但仍然無法克服安裝復(fù)雜、觀測難、受主觀因素干擾等問題。

隨著傳感器技術(shù)的發(fā)展以及計算機運算能力的增強，通過點云提取樣地單木及林分信息的方法得到快速發(fā)展。點云通常通過數(shù)字?jǐn)z影測量技術(shù)、ToF(Time of flight)技術(shù)等獲取。數(shù)字?jǐn)z影測量系統(tǒng)通常以相機作為觀測傳感器，利用相機記錄的紋理信息及角度重建被觀測目標(biāo)的三維采樣及三維模型，但由于森林中植被遮擋嚴(yán)重，很難單純利用攝影測量技術(shù)高效完成樣地調(diào)查；此外，由于林下紋理過于復(fù)雜，很難形成樣地完整三維點云[7-9]。ToF相機[10-13]由于其分辨率低、測距精度低且易受自然光線等影響，也不足以高效完成樣地調(diào)查[14-15]。

相比于普通相機及ToF相機，激光雷達(dá)(Light detection and ranging，LiDAR)具有高效、不受自然光線影響等特性，是近年來比較流行的樣地清查解決方案[16]。地基激光雷達(dá)(Terrestrial laser scanning，TLS)單次掃描由于僅有一個掃描視角，立木遮擋將無法獲取到較大樣地的完整三維數(shù)據(jù)，這導(dǎo)致所提取樣地屬性精度降低或無法獲取到完整的樣地屬性[17-18]。多次掃描雖然避免了數(shù)據(jù)遺漏問題，但合并各站點掃描數(shù)據(jù)具有一定難度，且在大樣地中布設(shè)太多站點耗時費力，TLS不適合在較大的樣地中使用[19-20]。移動激光雷達(dá)系統(tǒng)(Mobile laser scanning，MLS)通常指以GNSS(Global navigation satellite system)、IMU(Inertial measurement unit)、激光雷達(dá)為主要傳感器的運動平臺，其通過在樣地中動態(tài)采集樣地數(shù)據(jù)，后處理時利用GNSS+IMU組合獲取的姿軌曲線拼接激光雷達(dá)數(shù)據(jù)[21-22]；但通常森林郁閉度較高，受樹冠遮擋影響，GNSS接收機很難在林下正常工作。即時定位及建圖(Simultaneous localization and mapping，SLAM)技術(shù)的出現(xiàn)使得在林下無GNSS信號的區(qū)域定位成為可能[23-26]。隨著SLAM算法的改進(jìn)，單臺多線激光雷達(dá)采集數(shù)據(jù)即可完成定位及建圖，單激光雷達(dá)相比于TLS、MLS等具有設(shè)備簡單、成本低等優(yōu)勢，其中一種典型的LiDAR SLAM算法為LOAM(LiDAR odometry and mapping in real-time)[27-30]，該算法以單幀數(shù)據(jù)中面特征、線特征為單幀去畸變、配準(zhǔn)等過程的特征，具有計算速度快等優(yōu)勢；由于森林中線、面特征較少，單幀去畸變過程可能引入較大誤差，導(dǎo)致配準(zhǔn)精度低、構(gòu)建點云噪聲大等，故很難將該算法直接用于森林調(diào)查中。

針對以上問題，本文以LOAM為基礎(chǔ)構(gòu)建森林樣地調(diào)查系統(tǒng)。在SLAM系統(tǒng)工作流中引入二次去畸變、二次配準(zhǔn)等模塊以提高去畸變、配準(zhǔn)的魯棒性及精度；使用32線激光雷達(dá)掃描4塊32 m×32 m的森林樣地，利用改進(jìn)的LOAM系統(tǒng)完成樣地建圖；從點云提取立木位置及胸徑，并通過參考數(shù)據(jù)及LOAM估計結(jié)果對比，完成LiDAR SLAM系統(tǒng)在森林樣地中建圖精度的間接評估。

1 LOAM算法及其變體

在處理連續(xù)單幀LiDAR數(shù)據(jù)時(以第k幀Pk為例)具體步驟為：

(1)基于曲率提取單幀點云數(shù)據(jù)中的線、面特征。

(2)假設(shè)Pk掃描期間雷達(dá)做勻速運動，即若tk時刻至tk+1時刻雷達(dá)位姿變換量為

(1)

則tk時刻至ti時刻雷達(dá)位姿變換量可線性內(nèi)插為

(2)

圖1 LOAM及其變體工作流程圖Fig.1 Workflow of LOAM and its alternatives

從LOAM的結(jié)構(gòu)中可以看出其主要結(jié)構(gòu)為雷達(dá)里程計，并不包含傳統(tǒng)SLAM系統(tǒng)的后端(回環(huán)檢測及圖優(yōu)化)，故并不適用于較大的建圖場景。

LeGO-LOAM(Lightweight and ground-optimized LiDAR odometry and mapping on variable terrain)作為一種重要的LOAM變更體[28]，主要進(jìn)行了如下改進(jìn)：①將單幀點云中地面及周圍點云進(jìn)行了分類。②在特征提取時，提取地面點云中的面特征、周圍點云中的線特征。③在去畸變時，利用面特征完成豎直方向的線元素及翻滾角、俯仰角的優(yōu)化，利用線特征完成其他元素的優(yōu)化。④以ICP(Iterative closest point)算法為基礎(chǔ)構(gòu)建了回環(huán)檢測功能，并添加了圖優(yōu)化功能。相比于LOAM，該系統(tǒng)在室內(nèi)或有建筑物的室外等場景建圖精度更佳，且適用于更大規(guī)模的建圖場景。

2 LiDAR.SLAM算法改進(jìn)

相比于室內(nèi)或有建筑物的室外場景，森林中具有較少的線、面特征，故使用特征進(jìn)行去畸變及配準(zhǔn)精度可能性相對較低。本文以LOAM為基礎(chǔ)構(gòu)建了適合于森林條件的LiDAR SLAM算法，其工作流程如圖2所示。

圖2 森林LOAM工作流程圖Fig.2 Workflow of LOAM for forest inventory

2.1 SLAM系統(tǒng)工作流程設(shè)計

針對森林中具有較少線、面特征的特點，對SLAM系統(tǒng)工作流程進(jìn)行優(yōu)化，以提高SLAM系統(tǒng)位姿估計魯棒性并提高構(gòu)建三維點云精度。相比于LOAM，本文構(gòu)建系統(tǒng)進(jìn)行了優(yōu)化：

(1)在單幀點云特征提取時，將點云分為3類，即用于配準(zhǔn)的下采樣線特征、下采樣面特征及其它點，在保存關(guān)鍵幀數(shù)據(jù)時，將用于配準(zhǔn)的線、面特征存在內(nèi)存中，其它點(占單幀點云總數(shù)約80%)則緩存至本地，相比于傳統(tǒng)LOAM及LeGO-LOAM，該方案可構(gòu)建密度更大、范圍更大的三維點云。

(2)所使用線特征僅保留連續(xù)點形成的曲率較大點，遮擋形成線點被剔除(即僅保留兩側(cè)連續(xù)點與線特征點的深度差小于30 cm的特征)，防止立木切線被作為線特征參與運算。

(4)采用基于里程計的回環(huán)檢測及基于表面的回環(huán)檢測(Scan Context++[29])完成回環(huán)幀的判別，然后使用基于線、面特征的配準(zhǔn)獲取回環(huán)幀與當(dāng)前幀約束；相比于ICP算法，基于線、面特征的特征關(guān)聯(lián)方法更高效且可靠性更強。

(5)不僅為其構(gòu)建了UI界面，通過進(jìn)一步改造，使其擺脫ROS環(huán)境，可在Windows、Ubuntu等操作系統(tǒng)下編譯并運行(圖3)。

圖3 森林LiDAR SLAM系統(tǒng)界面及其操作Fig.3 Forest LiDAR SLAM system interface and its operations

2.2 帶有先驗條件的SLAM優(yōu)化算法設(shè)計

LOAM及其變體進(jìn)行去畸變及配準(zhǔn)優(yōu)化中，將當(dāng)前幀線/面特征與參考幀(或局部地圖)特征進(jìn)行關(guān)聯(lián)后基于點到線、面的距離最小化完成位姿估計，點到線、面的距離可表達(dá)為

(3)

(4)

以上方程均可線性化為廣義平差誤差方程

Av=BX-l

(5)

其中

式中A——單位矩陣

v——點到線或面的偶然誤差

B——線性化矩陣(即dL或dS對位姿線元素及角元素的導(dǎo)數(shù))

l——給定初始位姿X0條件下點到線、面的距離負(fù)值

X——待求初始位姿的改正數(shù)

利用以上廣義平差誤差方程式，便可利用牛頓迭代或買夸特算法等估計待求改正數(shù)，從而逼近位姿真值，實現(xiàn)去畸變或位姿估計。

將激光雷達(dá)的先驗條件等用于“去畸變”及“配準(zhǔn)”算法，以提高位姿估計及建圖精度。在“去畸變”優(yōu)化中，線特征平差誤差公式中參數(shù)可表達(dá)為

(6)

(7)

l=[-dL]

(8)

(9)

式中vdL——線特征到線的偶然誤差

面特征平差誤差公式中參數(shù)可表達(dá)為

(10)

(11)

l=[-dS]

(12)

(13)

式中vdS——面特征到面的偶然誤差

(14)

轉(zhuǎn)換為直角坐標(biāo)系值。極坐標(biāo)觀測值的精度在激光雷達(dá)說明書中已經(jīng)給出，或通過使用前檢?？色@取。若令

(15)

ΣX=MΣpMT

(16)

本文SLAM系統(tǒng)在“配準(zhǔn)”時，將去畸變后線、面特征精度(協(xié)方差陣)代入到配準(zhǔn)優(yōu)化過程中，其中線特征配準(zhǔn)平差誤差公式參數(shù)可表達(dá)為

(17)

(18)

l=[-dL]

(19)

(20)

面特征配準(zhǔn)平差誤差公式參數(shù)可表達(dá)為

(21)

(22)

l=[-dS]

(23)

(24)

3 三維點云構(gòu)建與立木檢尺參數(shù)提取

3.1 研究區(qū)域概況

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

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

3.2 數(shù)據(jù)采集與處理

選擇Velodyne VLP-32C型激光雷達(dá)為數(shù)據(jù)采集設(shè)備，該設(shè)備測距范圍為200 m、典型場景下測量精度為±3 cm、垂直視場角為40°、水平視場角為360°、水平角分辨率最大可達(dá)0.1°(本文選擇最高分辨率下進(jìn)行樣地掃描)。為方便手持并減少遮擋，為其加裝了手持柄；為判斷激光雷達(dá)豎直狀態(tài)，在激光雷達(dá)頂端安裝了萬向水準(zhǔn)器。改裝后設(shè)備如圖4所示。

圖4 用于樣地掃描的激光雷達(dá)系統(tǒng)Fig.4 LiDAR system for scanning forest plots1.水準(zhǔn)儀 2.激光雷達(dá) 3.數(shù)據(jù)采集系統(tǒng) 4.手持柄 5.移動電源

為保證樣地掃描點云完整性，并最大限度利用SLAM系統(tǒng)回環(huán)檢測減少位姿漂移，以提高建圖精度，設(shè)計了固定的樣地掃描路徑(圖5)，即掃描起點為樣地中心，沿西北方向開始掃描；到達(dá)西北角點后開始類航空攝影測量航線式掃描，掃描平行軌跡間距為6 m；完成類航線掃描后到達(dá)東南角，最后還需回到樣地中心。修正掃描路徑：①相鄰平行軌跡間利用回環(huán)檢測實現(xiàn)局部位姿漂移修正。②軌跡起始段及終止段軌跡與類航線相交，利用回環(huán)檢測可實現(xiàn)局部位姿漂移修正。③軌跡終止點與起始點重合，利用回環(huán)檢測實現(xiàn)全局掃描路徑位姿漂移修正。

圖5 樣地掃描路徑Fig.5 Plot scanning path

為保證激光雷達(dá)坐標(biāo)系與其他參考坐標(biāo)系的轉(zhuǎn)換，具體樣地數(shù)據(jù)采集過程為：①在樣地中心及正北方向邊界處放置激光反射標(biāo)志(圖6)。②手持激光雷達(dá)至樣地中心，保持水平后開啟數(shù)據(jù)采集。③沿規(guī)劃路徑完成樣地掃描。

圖6 激光反射標(biāo)志Fig.6 Laser reflection mark

在完成森林樣地掃描后，將數(shù)據(jù)導(dǎo)入本文構(gòu)建LOAM系統(tǒng)完成樣地三維點云構(gòu)建；然后，使用LiDAR 360軟件通過坐標(biāo)變換、地面提取、DEM生成與高度歸一化、胸高提取及胸高圓柱體擬合等，完成樣地立木胸徑及位置提取(圖7)。

圖7 森林樣地點云后處理Fig.7 Post-processing of forest plot point cloud

3.3 研究方法

利用LiDAR SLAM系統(tǒng)及LOAM系統(tǒng)構(gòu)建了森林樣地三維點云。然后，將點云中提取立木胸徑及立木位置與參考值對比，實現(xiàn)對SLAM系統(tǒng)生產(chǎn)點云精度的間接評估；通過對比SLAM系統(tǒng)與LOAM系統(tǒng)后處理數(shù)據(jù)統(tǒng)計結(jié)果，評價本文所構(gòu)建SLAM系統(tǒng)。本文使用胸徑尺測量胸徑作為胸徑參考值，全站儀(Leica TS60型)測量立木位置數(shù)據(jù)與胸徑參考值聯(lián)合計算立木位置為立木位置參考值。在使用全站儀對立木位置進(jìn)行測量時，將全站儀架設(shè)于樣地中心(即與布設(shè)于樣地中心的激光反射標(biāo)志十字中心對齊)，通過瞄準(zhǔn)布設(shè)于樣地正北方向的反射標(biāo)志完成北向初始化。若由于遮擋，部分立木無法被觀測，可通過“合并多站”方式完成所有立木位置觀測。本文采用由立木位置誤差統(tǒng)計獲取的均值向量μ及二維協(xié)方差陣Σ對立木位置估計值精度進(jìn)行評價，其定義式為

(25)

(26)

式中 (xi，yi)——立木位置估計值

(xir，yir)——參考值

n——立木總數(shù)

此外，協(xié)方差陣Σ的最大特征值σmax可表示為

(27)

式中 eigenvalues()——Σ所有特征值

該值為立木位置最大變異性方向的標(biāo)準(zhǔn)差，本文用其評價立木位置精度。

使用偏差(BIAS)、均方根誤差(Root mean squared error, RMSE)、相對偏差(relBIAS)及相對均方根誤差(relRMSE)對胸徑估計值進(jìn)行評估(其中RMSE也用于立木位置精度評估)。

4 試驗結(jié)果分析

4.1 SLAM后處理數(shù)據(jù)定性評估

利用SLAM系統(tǒng)及LOAM系統(tǒng)分別對樣地原始幀數(shù)據(jù)進(jìn)行后處理，構(gòu)建森林三維密集點云數(shù)據(jù)。從定性角度看：①2種SLAM解算獲取的姿軌曲線并非完全重合，故2種解算方法獲取森林點云具有一定區(qū)別(圖8紅色軌跡為森林SLAM解算姿軌曲線，藍(lán)色軌跡為LOAM解算姿軌曲線；灰色點為樣地點云)。②SLAM剔除了遮擋線特征，避免立木與視點切線被作為線特征參與運算，采用回環(huán)檢測方法修正位姿圖，建圖漂移更小，不同關(guān)鍵幀拼接而成的胸高點云分層較LOAM系統(tǒng)小，故點云厚度較薄；而LOAM中將立木與視點切線特征作為線特征將產(chǎn)生誤匹配，胸高圓柱體中心往往被填充(圖9)。

圖8 SLAM與LOAM后處理姿軌曲線定性對比Fig.8 Comparison of post-processing trajectory curves between forest SLAM and LOAM

圖9 SLAM與LOAM后處理胸高點云數(shù)據(jù)定性對比Fig.9 Comparison of post-processing DBH point clouds between forest SLAM and LOAM

4.2 立木位置精度評估

由LiDAR SLAM系統(tǒng)及LOAM系統(tǒng)獲取立木位置估計值，如圖10所示，顯然LiDAR SLAM估計值較LOAM偏差小。表2統(tǒng)計結(jié)果表明：①SLAM系統(tǒng)解算結(jié)果中，各樣地x、y兩軸估計值的平均誤差為-0.029～0.027 m，即無明顯偏差；而LOAM解算結(jié)果雖然總體接近無偏差，但各樣地平均誤差較大(-0.101～0.092 m)。②SLAM系統(tǒng)解算結(jié)果中，x、y兩軸協(xié)方差值接近于0，說明兩軸間估計值無明顯相關(guān)性，LOAM解算結(jié)果有類似結(jié)論。③SLAM系統(tǒng)解算結(jié)果中，最大誤差方向的協(xié)方差σmax(0.066～0.098 m)及RMSE(0.052～0.104 m)均明顯小于LOAM計算結(jié)果(σmax為0.112～0.148 m，RMSE為0.114～0.144 m)，顯然森林SLAM估計值變異性較小(圖11)。

圖10 立木位置圖Fig.10 Estimated and reference stem positions

表2 立木位置估計值的精度Tab.2 Accuracies of stem position estimations

圖11 所有樣地立木位置誤差Fig.11 Position errors of all trees in plots

4.3 胸徑精度評估

由LiDAR SLAM系統(tǒng)及LOAM解算立木胸徑估計值如圖12所示，SLAM估計值較LOAM估計值更接近于兩軸中線，顯然SLAM胸徑估計值可靠性更強。表3統(tǒng)計結(jié)果顯示，SLAM解算各樣地胸徑估計值BIAS為-0.08～0.65 cm(relBIAS為-0.59%～2.42%)，與LOAM估計結(jié)果(BIAS為1.89～2.62 cm，relBIAS為10.09%～14.42%)相比具有較小偏差；RMSE為0.087～1.23 cm(relRMSE為3.61%～4.94%)，顯然較LOAM估計值(RMSE為3.37～3.86 cm，relRMSE為14.48%～25.08%)精度更高。從圖13可以看出，不同徑階SLAM立木胸徑估計值均具有較小誤差且變異性比較一致；而不同徑階LOAM胸徑估計值雖變異性比較一致，但有約2 cm的估計偏差。

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

圖13 不同徑階組胸徑誤差箱形圖Fig.13 Errors of DBH observations under different DBH

表3 胸徑估計值Tab.3 Accuracies of DBH estimations

5 結(jié)論

(1)針對森林環(huán)境中線、面特征少，無法精確建圖的缺點，構(gòu)建了LiDAR SLAM森林樣地調(diào)查系統(tǒng)，以便僅利用單臺多線激光雷達(dá)可精確完成森林樣地調(diào)查。該SLAM算法利用二次去畸變、二次配準(zhǔn)改善了位姿估計及點云地圖的魯棒性；該系統(tǒng)將激光雷達(dá)測量精度、位姿估計精度等先驗信息融入去畸變及配準(zhǔn)算法中，提高了位姿估計及點云地圖精度。

(2)該系統(tǒng)在4塊32 m×32 m方形樣地中進(jìn)行了測試，通過從該系統(tǒng)產(chǎn)生點云中提取的立木位置和胸徑對系統(tǒng)精度間接評估。經(jīng)與參考值進(jìn)行對比表明，該系統(tǒng)所獲取點云在立木位置及胸徑估計時均能獲取比LOAM算法精度更高的結(jié)果。顯然，改進(jìn)的LiDAR SLAM算法使單臺多線激光雷達(dá)高精度完成森林樣地調(diào)查成為可能。