范永祥 馮仲科 申朝永 閆 飛 蘇玨穎 王 蔚
(1.季華實(shí)驗(yàn)室, 佛山 528000; 2.北京林業(yè)大學(xué)精準(zhǔn)林業(yè)北京市重點(diǎn)實(shí)驗(yàn)室, 北京 100083;3.貴州省第三測(cè)繪院, 貴陽(yáng) 550004)
森林資源及其監(jiān)測(cè)所獲取的可靠信息可以為不同層級(jí)、不同目標(biāo)的林業(yè)部門(mén)制定可持續(xù)發(fā)展政策、投資決策等提供支持,從而為森林培育、森林采伐、氣候變化影響評(píng)估、火災(zāi)模型和碳儲(chǔ)量估算等奠定基礎(chǔ)[1]。森林資源清查是獲取森林資源信息的重要途徑。由于不同層級(jí)的使用者對(duì)森林資源信息的使用目的不同,所關(guān)心的森林資源屬性也不盡相同[2]。為滿(mǎn)足清查目的,需要根據(jù)森林資源及所要求的精度制定不同的調(diào)查策略。樣地調(diào)查是一種重要的森林資源清查方法,其通過(guò)在研究區(qū)域合理布設(shè)樣地,然后對(duì)樣地的森林屬性進(jìn)行調(diào)查、統(tǒng)計(jì),從而反演出整個(gè)區(qū)域總體或平均的森林屬性[3]。通常樣地直接調(diào)查的屬性包含樹(shù)種、胸徑、樹(shù)高及立木位置等,生物量、碳儲(chǔ)量等則可通過(guò)以上屬性構(gòu)建的反演模型進(jìn)行估計(jì)[4]。傳統(tǒng)樣地調(diào)查中,使用一些簡(jiǎn)單的儀器進(jìn)行森林調(diào)查(如使用胸徑尺或卡尺完成胸徑測(cè)量,使用羅盤(pán)及皮尺等完成立木位置測(cè)量),這些方法不僅耗時(shí)耗力,且無(wú)法克服觀測(cè)者主觀因素的干擾。隨著測(cè)距、測(cè)角等相關(guān)技術(shù)的發(fā)展,Laser-relascope[5]、電子測(cè)樹(shù)槍[6]等設(shè)備提供了同時(shí)測(cè)量樹(shù)木胸徑、樹(shù)高、位置等的一體化解決方案,但仍然無(wú)法克服安裝復(fù)雜、觀測(cè)難、受主觀因素干擾等問(wèn)題。
隨著傳感器技術(shù)的發(fā)展以及計(jì)算機(jī)運(yùn)算能力的增強(qiáng),通過(guò)點(diǎn)云提取樣地單木及林分信息的方法得到快速發(fā)展。點(diǎn)云通常通過(guò)數(shù)字?jǐn)z影測(cè)量技術(shù)、ToF(Time of flight)技術(shù)等獲取。數(shù)字?jǐn)z影測(cè)量系統(tǒng)通常以相機(jī)作為觀測(cè)傳感器,利用相機(jī)記錄的紋理信息及角度重建被觀測(cè)目標(biāo)的三維采樣及三維模型,但由于森林中植被遮擋嚴(yán)重,很難單純利用攝影測(cè)量技術(shù)高效完成樣地調(diào)查;此外,由于林下紋理過(guò)于復(fù)雜,很難形成樣地完整三維點(diǎn)云[7-9]。ToF相機(jī)[10-13]由于其分辨率低、測(cè)距精度低且易受自然光線等影響,也不足以高效完成樣地調(diào)查[14-15]。
相比于普通相機(jī)及ToF相機(jī),激光雷達(dá)(Light detection and ranging,LiDAR)具有高效、不受自然光線影響等特性,是近年來(lái)比較流行的樣地清查解決方案[16]。地基激光雷達(dá)(Terrestrial laser scanning,TLS)單次掃描由于僅有一個(gè)掃描視角,立木遮擋將無(wú)法獲取到較大樣地的完整三維數(shù)據(jù),這導(dǎo)致所提取樣地屬性精度降低或無(wú)法獲取到完整的樣地屬性[17-18]。多次掃描雖然避免了數(shù)據(jù)遺漏問(wèn)題,但合并各站點(diǎn)掃描數(shù)據(jù)具有一定難度,且在大樣地中布設(shè)太多站點(diǎn)耗時(shí)費(fèi)力,TLS不適合在較大的樣地中使用[19-20]。移動(dòng)激光雷達(dá)系統(tǒng)(Mobile laser scanning,MLS)通常指以GNSS(Global navigation satellite system)、IMU(Inertial measurement unit)、激光雷達(dá)為主要傳感器的運(yùn)動(dòng)平臺(tái),其通過(guò)在樣地中動(dòng)態(tài)采集樣地?cái)?shù)據(jù),后處理時(shí)利用GNSS+IMU組合獲取的姿軌曲線拼接激光雷達(dá)數(shù)據(jù)[21-22];但通常森林郁閉度較高,受樹(shù)冠遮擋影響,GNSS接收機(jī)很難在林下正常工作。即時(shí)定位及建圖(Simultaneous localization and mapping,SLAM)技術(shù)的出現(xiàn)使得在林下無(wú)GNSS信號(hào)的區(qū)域定位成為可能[23-26]。隨著SLAM算法的改進(jìn),單臺(tái)多線激光雷達(dá)采集數(shù)據(jù)即可完成定位及建圖,單激光雷達(dá)相比于TLS、MLS等具有設(shè)備簡(jiǎn)單、成本低等優(yōu)勢(shì),其中一種典型的LiDAR SLAM算法為L(zhǎng)OAM(LiDAR odometry and mapping in real-time)[27-30],該算法以單幀數(shù)據(jù)中面特征、線特征為單幀去畸變、配準(zhǔn)等過(guò)程的特征,具有計(jì)算速度快等優(yōu)勢(shì);由于森林中線、面特征較少,單幀去畸變過(guò)程可能引入較大誤差,導(dǎo)致配準(zhǔn)精度低、構(gòu)建點(diǎn)云噪聲大等,故很難將該算法直接用于森林調(diào)查中。
針對(duì)以上問(wèn)題,本文以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)完成樣地建圖;從點(diǎn)云提取立木位置及胸徑,并通過(guò)參考數(shù)據(jù)及LOAM估計(jì)結(jié)果對(duì)比,完成LiDAR SLAM系統(tǒng)在森林樣地中建圖精度的間接評(píng)估。
在處理連續(xù)單幀LiDAR數(shù)據(jù)時(shí)(以第k幀Pk為例)具體步驟為:
(1)基于曲率提取單幀點(diǎn)云數(shù)據(jù)中的線、面特征。
(2)假設(shè)Pk掃描期間雷達(dá)做勻速運(yùn)動(dòng),即若tk時(shí)刻至tk+1時(shí)刻雷達(dá)位姿變換量為
(1)
則tk時(shí)刻至ti時(shí)刻雷達(dá)位姿變換量可線性?xún)?nèi)插為
(2)
圖1 LOAM及其變體工作流程圖Fig.1 Workflow of LOAM and its alternatives
從LOAM的結(jié)構(gòu)中可以看出其主要結(jié)構(gòu)為雷達(dá)里程計(jì),并不包含傳統(tǒng)SLAM系統(tǒng)的后端(回環(huán)檢測(cè)及圖優(yōu)化),故并不適用于較大的建圖場(chǎng)景。
LeGO-LOAM(Lightweight and ground-optimized LiDAR odometry and mapping on variable terrain)作為一種重要的LOAM變更體[28],主要進(jìn)行了如下改進(jìn):①將單幀點(diǎn)云中地面及周?chē)c(diǎn)云進(jìn)行了分類(lèi)。②在特征提取時(shí),提取地面點(diǎn)云中的面特征、周?chē)c(diǎn)云中的線特征。③在去畸變時(shí),利用面特征完成豎直方向的線元素及翻滾角、俯仰角的優(yōu)化,利用線特征完成其他元素的優(yōu)化。④以ICP(Iterative closest point)算法為基礎(chǔ)構(gòu)建了回環(huán)檢測(cè)功能,并添加了圖優(yōu)化功能。相比于LOAM,該系統(tǒng)在室內(nèi)或有建筑物的室外等場(chǎng)景建圖精度更佳,且適用于更大規(guī)模的建圖場(chǎng)景。
相比于室內(nèi)或有建筑物的室外場(chǎng)景,森林中具有較少的線、面特征,故使用特征進(jìn)行去畸變及配準(zhǔn)精度可能性相對(duì)較低。本文以LOAM為基礎(chǔ)構(gòu)建了適合于森林條件的LiDAR SLAM算法,其工作流程如圖2所示。
圖2 森林LOAM工作流程圖Fig.2 Workflow of LOAM for forest inventory
針對(duì)森林中具有較少線、面特征的特點(diǎn),對(duì)SLAM系統(tǒng)工作流程進(jìn)行優(yōu)化,以提高SLAM系統(tǒng)位姿估計(jì)魯棒性并提高構(gòu)建三維點(diǎn)云精度。相比于LOAM,本文構(gòu)建系統(tǒng)進(jìn)行了優(yōu)化:
(1)在單幀點(diǎn)云特征提取時(shí),將點(diǎn)云分為3類(lèi),即用于配準(zhǔn)的下采樣線特征、下采樣面特征及其它點(diǎn),在保存關(guān)鍵幀數(shù)據(jù)時(shí),將用于配準(zhǔn)的線、面特征存在內(nèi)存中,其它點(diǎn)(占單幀點(diǎn)云總數(shù)約80%)則緩存至本地,相比于傳統(tǒng)LOAM及LeGO-LOAM,該方案可構(gòu)建密度更大、范圍更大的三維點(diǎn)云。
(2)所使用線特征僅保留連續(xù)點(diǎn)形成的曲率較大點(diǎn),遮擋形成線點(diǎn)被剔除(即僅保留兩側(cè)連續(xù)點(diǎn)與線特征點(diǎn)的深度差小于30 cm的特征),防止立木切線被作為線特征參與運(yùn)算。
(4)采用基于里程計(jì)的回環(huán)檢測(cè)及基于表面的回環(huán)檢測(cè)(Scan Context++[29])完成回環(huán)幀的判別,然后使用基于線、面特征的配準(zhǔn)獲取回環(huán)幀與當(dāng)前幀約束;相比于ICP算法,基于線、面特征的特征關(guān)聯(lián)方法更高效且可靠性更強(qiáng)。
(5)不僅為其構(gòu)建了UI界面,通過(guò)進(jìn)一步改造,使其擺脫ROS環(huán)境,可在Windows、Ubuntu等操作系統(tǒng)下編譯并運(yùn)行(圖3)。
圖3 森林LiDAR SLAM系統(tǒng)界面及其操作Fig.3 Forest LiDAR SLAM system interface and its operations
LOAM及其變體進(jìn)行去畸變及配準(zhǔn)優(yōu)化中,將當(dāng)前幀線/面特征與參考幀(或局部地圖)特征進(jìn)行關(guān)聯(lián)后基于點(diǎn)到線、面的距離最小化完成位姿估計(jì),點(diǎn)到線、面的距離可表達(dá)為
(3)
(4)
以上方程均可線性化為廣義平差誤差方程
Av=BX-l
(5)
其中
式中A——單位矩陣
v——點(diǎn)到線或面的偶然誤差
B——線性化矩陣(即dL或dS對(duì)位姿線元素及角元素的導(dǎo)數(shù))
l——給定初始位姿X0條件下點(diǎn)到線、面的距離負(fù)值
X——待求初始位姿的改正數(shù)
利用以上廣義平差誤差方程式,便可利用牛頓迭代或買(mǎi)夸特算法等估計(jì)待求改正數(shù),從而逼近位姿真值,實(shí)現(xiàn)去畸變或位姿估計(jì)。
將激光雷達(dá)的先驗(yàn)條件等用于“去畸變”及“配準(zhǔn)”算法,以提高位姿估計(jì)及建圖精度。在“去畸變”優(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)觀測(cè)值的精度在激光雷達(dá)說(shuō)明書(shū)中已經(jīng)給出,或通過(guò)使用前檢校可獲取。若令
(15)
ΣX=MΣpMT
(16)
本文SLAM系統(tǒng)在“配準(zhǔn)”時(shí),將去畸變后線、面特征精度(協(xié)方差陣)代入到配準(zhǔn)優(yōu)化過(guò)程中,其中線特征配準(zhǔn)平差誤差公式參數(shù)可表達(dá)為
(17)
(18)
l=[-dL]
(19)
(20)
面特征配準(zhǔn)平差誤差公式參數(shù)可表達(dá)為
(21)
(22)
l=[-dS]
(23)
(24)
選擇位于北京市近郊西山試驗(yàn)林場(chǎng)(39°58′N(xiāo),116°11′E)為研究區(qū)域,該林場(chǎng)主要以人工林為主體。選擇其中4塊32 m×32 m的方形樣地為研究對(duì)象。所選樣地地面具有較少灌木且容易到達(dá),不同徑階組比較均衡。表1總結(jié)了不同樣地的基本屬性。
表1 樣地屬性的描述統(tǒng)計(jì)Tab.1 Summary statistics of plot attributes
選擇Velodyne VLP-32C型激光雷達(dá)為數(shù)據(jù)采集設(shè)備,該設(shè)備測(cè)距范圍為200 m、典型場(chǎng)景下測(cè)量精度為±3 cm、垂直視場(chǎng)角為40°、水平視場(chǎng)角為360°、水平角分辨率最大可達(dá)0.1°(本文選擇最高分辨率下進(jìn)行樣地掃描)。為方便手持并減少遮擋,為其加裝了手持柄;為判斷激光雷達(dá)豎直狀態(tài),在激光雷達(dá)頂端安裝了萬(wàn)向水準(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.移動(dòng)電源
為保證樣地掃描點(diǎn)云完整性,并最大限度利用SLAM系統(tǒng)回環(huán)檢測(cè)減少位姿漂移,以提高建圖精度,設(shè)計(jì)了固定的樣地掃描路徑(圖5),即掃描起點(diǎn)為樣地中心,沿西北方向開(kāi)始掃描;到達(dá)西北角點(diǎn)后開(kāi)始類(lèi)航空攝影測(cè)量航線式掃描,掃描平行軌跡間距為6 m;完成類(lèi)航線掃描后到達(dá)東南角,最后還需回到樣地中心。修正掃描路徑:①相鄰平行軌跡間利用回環(huán)檢測(cè)實(shí)現(xiàn)局部位姿漂移修正。②軌跡起始段及終止段軌跡與類(lèi)航線相交,利用回環(huán)檢測(cè)可實(shí)現(xiàn)局部位姿漂移修正。③軌跡終止點(diǎn)與起始點(diǎn)重合,利用回環(huán)檢測(cè)實(shí)現(xiàn)全局掃描路徑位姿漂移修正。
圖5 樣地掃描路徑Fig.5 Plot scanning path
為保證激光雷達(dá)坐標(biāo)系與其他參考坐標(biāo)系的轉(zhuǎn)換,具體樣地?cái)?shù)據(jù)采集過(guò)程為:①在樣地中心及正北方向邊界處放置激光反射標(biāo)志(圖6)。②手持激光雷達(dá)至樣地中心,保持水平后開(kāi)啟數(shù)據(jù)采集。③沿規(guī)劃路徑完成樣地掃描。
圖6 激光反射標(biāo)志Fig.6 Laser reflection mark
在完成森林樣地掃描后,將數(shù)據(jù)導(dǎo)入本文構(gòu)建LOAM系統(tǒng)完成樣地三維點(diǎn)云構(gòu)建;然后,使用LiDAR 360軟件通過(guò)坐標(biāo)變換、地面提取、DEM生成與高度歸一化、胸高提取及胸高圓柱體擬合等,完成樣地立木胸徑及位置提取(圖7)。
圖7 森林樣地點(diǎn)云后處理Fig.7 Post-processing of forest plot point cloud
利用LiDAR SLAM系統(tǒng)及LOAM系統(tǒng)構(gòu)建了森林樣地三維點(diǎn)云。然后,將點(diǎn)云中提取立木胸徑及立木位置與參考值對(duì)比,實(shí)現(xiàn)對(duì)SLAM系統(tǒng)生產(chǎn)點(diǎn)云精度的間接評(píng)估;通過(guò)對(duì)比SLAM系統(tǒng)與LOAM系統(tǒng)后處理數(shù)據(jù)統(tǒng)計(jì)結(jié)果,評(píng)價(jià)本文所構(gòu)建SLAM系統(tǒng)。本文使用胸徑尺測(cè)量胸徑作為胸徑參考值,全站儀(Leica TS60型)測(cè)量立木位置數(shù)據(jù)與胸徑參考值聯(lián)合計(jì)算立木位置為立木位置參考值。在使用全站儀對(duì)立木位置進(jìn)行測(cè)量時(shí),將全站儀架設(shè)于樣地中心(即與布設(shè)于樣地中心的激光反射標(biāo)志十字中心對(duì)齊),通過(guò)瞄準(zhǔn)布設(shè)于樣地正北方向的反射標(biāo)志完成北向初始化。若由于遮擋,部分立木無(wú)法被觀測(cè),可通過(guò)“合并多站”方式完成所有立木位置觀測(cè)。本文采用由立木位置誤差統(tǒng)計(jì)獲取的均值向量μ及二維協(xié)方差陣Σ對(duì)立木位置估計(jì)值精度進(jìn)行評(píng)價(jià),其定義式為
(25)
(26)
式中 (xi,yi)——立木位置估計(jì)值
(xir,yir)——參考值
n——立木總數(shù)
此外,協(xié)方差陣Σ的最大特征值σmax可表示為
(27)
式中 eigenvalues()——Σ所有特征值
該值為立木位置最大變異性方向的標(biāo)準(zhǔn)差,本文用其評(píng)價(jià)立木位置精度。
使用偏差(BIAS)、均方根誤差(Root mean squared error, RMSE)、相對(duì)偏差(relBIAS)及相對(duì)均方根誤差(relRMSE)對(duì)胸徑估計(jì)值進(jìn)行評(píng)估(其中RMSE也用于立木位置精度評(píng)估)。
利用SLAM系統(tǒng)及LOAM系統(tǒng)分別對(duì)樣地原始幀數(shù)據(jù)進(jìn)行后處理,構(gòu)建森林三維密集點(diǎn)云數(shù)據(jù)。從定性角度看:①2種SLAM解算獲取的姿軌曲線并非完全重合,故2種解算方法獲取森林點(diǎn)云具有一定區(qū)別(圖8紅色軌跡為森林SLAM解算姿軌曲線,藍(lán)色軌跡為L(zhǎng)OAM解算姿軌曲線;灰色點(diǎn)為樣地點(diǎn)云)。②SLAM剔除了遮擋線特征,避免立木與視點(diǎn)切線被作為線特征參與運(yùn)算,采用回環(huán)檢測(cè)方法修正位姿圖,建圖漂移更小,不同關(guān)鍵幀拼接而成的胸高點(diǎn)云分層較LOAM系統(tǒng)小,故點(diǎn)云厚度較?。欢鳯OAM中將立木與視點(diǎn)切線特征作為線特征將產(chǎn)生誤匹配,胸高圓柱體中心往往被填充(圖9)。
圖8 SLAM與LOAM后處理姿軌曲線定性對(duì)比Fig.8 Comparison of post-processing trajectory curves between forest SLAM and LOAM
圖9 SLAM與LOAM后處理胸高點(diǎn)云數(shù)據(jù)定性對(duì)比Fig.9 Comparison of post-processing DBH point clouds between forest SLAM and LOAM
由LiDAR SLAM系統(tǒng)及LOAM系統(tǒng)獲取立木位置估計(jì)值,如圖10所示,顯然LiDAR SLAM估計(jì)值較LOAM偏差小。表2統(tǒng)計(jì)結(jié)果表明:①SLAM系統(tǒng)解算結(jié)果中,各樣地x、y兩軸估計(jì)值的平均誤差為-0.029~0.027 m,即無(wú)明顯偏差;而LOAM解算結(jié)果雖然總體接近無(wú)偏差,但各樣地平均誤差較大(-0.101~0.092 m)。②SLAM系統(tǒng)解算結(jié)果中,x、y兩軸協(xié)方差值接近于0,說(shuō)明兩軸間估計(jì)值無(wú)明顯相關(guān)性,LOAM解算結(jié)果有類(lèi)似結(jié)論。③SLAM系統(tǒng)解算結(jié)果中,最大誤差方向的協(xié)方差σmax(0.066~0.098 m)及RMSE(0.052~0.104 m)均明顯小于LOAM計(jì)算結(jié)果(σmax為0.112~0.148 m,RMSE為0.114~0.144 m),顯然森林SLAM估計(jì)值變異性較小(圖11)。
圖10 立木位置圖Fig.10 Estimated and reference stem positions
表2 立木位置估計(jì)值的精度Tab.2 Accuracies of stem position estimations
圖11 所有樣地立木位置誤差Fig.11 Position errors of all trees in plots
由LiDAR SLAM系統(tǒng)及LOAM解算立木胸徑估計(jì)值如圖12所示,SLAM估計(jì)值較LOAM估計(jì)值更接近于兩軸中線,顯然SLAM胸徑估計(jì)值可靠性更強(qiáng)。表3統(tǒng)計(jì)結(jié)果顯示,SLAM解算各樣地胸徑估計(jì)值BIAS為-0.08~0.65 cm(relBIAS為-0.59%~2.42%),與LOAM估計(jì)結(jié)果(BIAS為1.89~2.62 cm,relBIAS為10.09%~14.42%)相比具有較小偏差;RMSE為0.087~1.23 cm(relRMSE為3.61%~4.94%),顯然較LOAM估計(jì)值(RMSE為3.37~3.86 cm,relRMSE為14.48%~25.08%)精度更高。從圖13可以看出,不同徑階SLAM立木胸徑估計(jì)值均具有較小誤差且變異性比較一致;而不同徑階LOAM胸徑估計(jì)值雖變異性比較一致,但有約2 cm的估計(jì)偏差。
圖12 胸徑估計(jì)值散點(diǎn)圖Fig.12 Scatter plot of estimated DBH
圖13 不同徑階組胸徑誤差箱形圖Fig.13 Errors of DBH observations under different DBH
表3 胸徑估計(jì)值Tab.3 Accuracies of DBH estimations
(1)針對(duì)森林環(huán)境中線、面特征少,無(wú)法精確建圖的缺點(diǎn),構(gòu)建了LiDAR SLAM森林樣地調(diào)查系統(tǒng),以便僅利用單臺(tái)多線激光雷達(dá)可精確完成森林樣地調(diào)查。該SLAM算法利用二次去畸變、二次配準(zhǔn)改善了位姿估計(jì)及點(diǎn)云地圖的魯棒性;該系統(tǒng)將激光雷達(dá)測(cè)量精度、位姿估計(jì)精度等先驗(yàn)信息融入去畸變及配準(zhǔn)算法中,提高了位姿估計(jì)及點(diǎn)云地圖精度。
(2)該系統(tǒng)在4塊32 m×32 m方形樣地中進(jìn)行了測(cè)試,通過(guò)從該系統(tǒng)產(chǎn)生點(diǎn)云中提取的立木位置和胸徑對(duì)系統(tǒng)精度間接評(píng)估。經(jīng)與參考值進(jìn)行對(duì)比表明,該系統(tǒng)所獲取點(diǎn)云在立木位置及胸徑估計(jì)時(shí)均能獲取比LOAM算法精度更高的結(jié)果。顯然,改進(jìn)的LiDAR SLAM算法使單臺(tái)多線激光雷達(dá)高精度完成森林樣地調(diào)查成為可能。