基于改進LOAM的森林樣地調查系統設計與試驗

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

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

0 引言

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

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

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

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

1 LOAM算法及其變體

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

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

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

(1)

則tk時刻至ti時刻雷達位姿變換量可線性內插為

(2)

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

從LOAM的結構中可以看出其主要結構為雷達里程計，并不包含傳統SLAM系統的后端(回環檢測及圖優化)，故并不適用于較大的建圖場景。

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

2 LiDAR.SLAM算法改進

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

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

2.1 SLAM系統工作流程設計

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

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

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

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

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

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

2.2 帶有先驗條件的SLAM優化算法設計

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

(3)

(4)

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

Av=BX-l

(5)

其中

式中A——單位矩陣

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

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

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

X——待求初始位姿的改正數

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

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

(6)

(7)

l=[-dL]

(8)

(9)

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

面特征平差誤差公式中參數可表達為

(10)

(11)

l=[-dS]

(12)

(13)

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

(14)

轉換為直角坐標系值。極坐標觀測值的精度在激光雷達說明書中已經給出，或通過使用前檢校可獲取。若令

(15)

ΣX=MΣpMT

(16)

本文SLAM系統在“配準”時，將去畸變后線、面特征精度(協方差陣)代入到配準優化過程中，其中線特征配準平差誤差公式參數可表達為

(17)

(18)

l=[-dL]

(19)

(20)

面特征配準平差誤差公式參數可表達為

(21)

(22)

l=[-dS]

(23)

(24)

3 三維點云構建與立木檢尺參數提取

3.1 研究區域概況

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

表1 樣地屬性的描述統計Tab.1 Summary statistics of plot attributes

3.2 數據采集與處理

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

圖4 用于樣地掃描的激光雷達系統Fig.4 LiDAR system for scanning forest plots1.水準儀 2.激光雷達 3.數據采集系統 4.手持柄 5.移動電源

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

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

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

圖6 激光反射標志Fig.6 Laser reflection mark

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

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

3.3 研究方法

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

(25)

(26)

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

(xir，yir)——參考值

n——立木總數

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

(27)

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

該值為立木位置最大變異性方向的標準差，本文用其評價立木位置精度。

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

4 試驗結果分析

4.1 SLAM后處理數據定性評估

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

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

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

4.2 立木位置精度評估

由LiDAR SLAM系統及LOAM系統獲取立木位置估計值，如圖10所示，顯然LiDAR SLAM估計值較LOAM偏差小。表2統計結果表明：①SLAM系統解算結果中，各樣地x、y兩軸估計值的平均誤差為-0.029～0.027 m，即無明顯偏差；而LOAM解算結果雖然總體接近無偏差，但各樣地平均誤差較大(-0.101～0.092 m)。②SLAM系統解算結果中，x、y兩軸協方差值接近于0，說明兩軸間估計值無明顯相關性，LOAM解算結果有類似結論。③SLAM系統解算結果中，最大誤差方向的協方差σmax(0.066～0.098 m)及RMSE(0.052～0.104 m)均明顯小于LOAM計算結果(σ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系統及LOAM解算立木胸徑估計值如圖12所示，SLAM估計值較LOAM估計值更接近于兩軸中線，顯然SLAM胸徑估計值可靠性更強。表3統計結果顯示，SLAM解算各樣地胸徑估計值BIAS為-0.08～0.65 cm(relBIAS為-0.59%～2.42%)，與LOAM估計結果(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 結論

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

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