梁 軒,王耀力,常 青,翟銀枝
(1.太原理工大學信息與計算機學院,山西晉中 030600;2.山西省晉中市榆次區國有烏金山林場,山西 晉中 030619)
傳統林木形態測量技術依靠人力測量,費時費力。通過無人機或地面車輛搭載相機[1-4],短時間內收集大量圖像,為實現自動化林木監測提供數據基礎。同時定位與建圖(SLAM)和運動推斷結構(SFM)方法[5-7]適用于靜態場景的重建,但無法解決動態變化的林木重建問題。構建時空地圖[8-9]可以提高時空數據關聯的準確性,這通常為自動駕駛應用設計,解決不同視角與不同時間的數據關聯方法。董靖[10]實現了農作物的時空重建,但前端實時提取SIFT 特征將對數據收集前端造成較大計算負擔。文獻[11-12]通過計算農作物重建結果的幾何特征,得出農作物生長趨勢數據,但分析模型根據農作物設計,難以直接應用于樹木重建結果分析中。
為此,提出了基于時空重建的森林林木形態監測算法,該算法由三部分組成:第一部分是搭建基于同步定位與建圖的視覺數據收集前端,用于估計相機位姿;第二部分是后端建立圖像間魯棒數據關聯,構建完整的時空重建;第三部分是設計通過稠密點云預處理,進行樹木形態變化監測。
V-SLAM 數據收集前端的傳感器是RGB-D 深度相機,RGB-D 深度相機基于針孔相機模型。該模型包括四個坐標系:世界坐標系、相機坐標系、成像平面坐標系與像素平面坐標系。假設現實世界空間點的三維坐標為(Xw,Yw,Zw)T,記為Pw。對應相機平面坐標系坐標為(Xc,Yc,Zc)T,記為Pc。對應成像平面坐標系的坐標為(Xi,Yi)T,記為Pi。對應像素平面坐標系的坐標為(u,v)T,記為P。設成像平面到相機光心O的距離(焦距)為f。
世界坐標系與相機坐標系轉換如圖1 所示。

圖1 世界坐標系與相機坐標系轉換
由模型推導出,世界坐標系的空間點轉換到相機平面坐標系下的公式為:
R為3×3 矩陣表征相機旋轉量,t表為3×1 矩陣表征相機平移量,都是系統待計算的數據,假設相機初始位置為空間原點時,旋轉與平移矩陣實際上代表相機位姿矩陣T。
相機平面坐標系下空間點坐標轉換到像素平面坐標系下齊次坐標的公式為:
K為相機內參矩陣,根據實際情況進行標定得到。圖2、3 是相機平面到像素平面的轉換過程示意圖。dx與dy是成像平面到像素平面橫軸與縱軸的放縮倍數,u0與v0是像素平面坐標原點在橫軸與縱軸平移量。

圖2 相機平面到成像平面轉換

圖3 成像平面到像素平面轉換
ORB 特征[13]在相機的旋轉與圖像的放縮情況下都有魯棒性,且在相機自動增益、自動曝光與明暗變化的情況具有良好不變性。
提取ORB 特征比提取SIFT 特征[14]耗費更少時間,具體見后面的實驗。數據收集前端運行過程中,系統實時性與魯棒性制約相機移動速度,常用方法是借助GPU 加速并行處理,但在數據收集前端中使用GPU 過于奢侈。搭載相機的無人機或者手持嵌入式設備方式下使用GPU 也不便捷。提取ORB 特征可以適應普通CPU 的系統需求,降低前端設備成本,具有更好的實時性。
ORB 特征提取后,采用PnP 匹配模式估計相機位姿變化,通過隨機抽樣一致算法去除大部分匹配錯誤的特征。采用的前端相機可以直接輸出深度圖,由此求出三維空間特征點坐標,克服了單目相機采用PnP 匹配模式時需要初始化的缺點。采用非線性方法來計算PnP 匹配下的相機位姿估計,具體細節見1.3 節。
在相機位姿估計中,采用PnP 匹配模式求解,把相機位姿估計看成一個非線性最小二乘問題。如圖4 所示,求解最小二乘問題最優解,實際上就是最小化空間點重投影誤差的過程。重投影誤差是估計值與測量值的差值,估計值是重投影的像素坐標,測量值是相機測得的像素坐標P,由以下公式表示:

圖4 重投影誤差原理
將多空間點重投影誤差最小二乘項相加求最小值,整體誤差取得最小值時,求出的相機位姿就是最優解。經過RANSAC 排除錯誤匹配對,且優化前相機位姿與真實值比較接近,采用最小二乘法最小化重投影誤差不容易陷入局部最優解中。
數據收集前端包括跟蹤模塊與局部建圖模塊兩個核心模塊。前端系統結構如圖5 所示。

圖5 前端系統結構
跟蹤模塊在相機當前幀中發現ORB 特征,完成PnP 匹配估計相機位姿。相機當前幀輸入前需要預處理,因為采用的是基于機器人操作平臺(ROS)的RGB-D 深度相機,RGB 模塊與深度模塊工作頻率通常不同,采用ROS 的多傳感器時間軟同步機制實現兩模塊輸入在相近時間內的非嚴格對準。跟蹤模塊對每幀圖像做關鍵幀準入判定,關鍵幀主要保存了ORB 特征與相機估計位姿。
局部建圖模塊接收跟蹤模塊中的關鍵幀,建立局部稀疏地圖。稀疏地圖由能在多關鍵幀都能觀測到的ORB 特征產生。通過維護全局地圖容器,該容器包含了最基本的地圖信息、關鍵幀信息及其連接信息,局部建圖模塊可以與之交互,將局部稀疏地圖通過必要的整合操作構成全局稀疏地圖,并對全局地圖中的關鍵幀做準出判定,減少需要維護的關鍵幀。
SFM[5-7]與SLAM 類似,前者更關注于特征點估計以完整恢復重建結構,后者更關注相機位姿估計以獲取更準確的定位信息。前端準確估計了相機位姿,后端將結合當前幀RGB-D 圖像,完成稠密點云重建。
SFM 是全局式重建,但前端收集數據卻是增量式的,這導致了累積誤差。采用基于詞袋的回環檢測[15]來降低累積誤差。基于關鍵幀的詞袋由每幀ORB 特征生成的描述向量經過k-means 聚類構成,通過樹狀數據結構查詢匹配,最后由TF-IDF進行相似度計算確認回環。確認回環后,使用2.3節的逆向投影圖形變化有界搜索方法,測出相機位姿作為先驗,重新完成時空數據關聯,解決從不同時間與不同視角觀察同一場景時圖像數據關聯的關鍵問題。
為建立三維重建之間的時空關聯,采取逆向投影有界搜索算法在不同時間不同視角測得的關鍵幀圖像之間建立數據關聯。其原理公式如下,設a、b為兩個時刻:
式(4)-(5)表示,對于空間同一點,已知深度因子、相機內參矩陣以及兩時刻相機位姿矩陣時,兩時刻對應像素坐標是確定的。為減少誤差影響,在該確定像素坐標(ub,vb)上設定一個有界范圍δ重新搜索該匹配像素點(u,v),以建立魯棒數據關聯。
后端會把輸入的相機位姿與稀疏地圖采用基于圖優化理論[15]的G2o 庫進行一次全局的光束法平差,其原理遵循1.3 節,最小化重投影誤差僅用于局部稀疏地圖建立,此處是全局稀疏地圖優化。如圖6 所示,虛線代表相機的運動軌跡,三角形代表單次位姿估計,表征相機運動模型。點虛線代表相機對ORB特征點對應空間點的觀測(如Pw2可以在b、c、d時刻的相機位姿中觀測到),觀測模型對應1.3 節提到的估計值。

圖6 圖優化模型
如圖7 所示,后端將前端輸入數據進行優化后開始稠密重建,通過時空重建模塊完成稠密點云時空數據關聯,最后輸出帶有時間信息的稠密點云到點云后期處理模塊。

圖7 后端系統結構與后期處理模塊
時空重建后端輸出的時空稠密點云帶有許多重建噪聲,通過基于半徑的異常點濾波[17]去除大部分重建噪聲。指定內點搜索半徑與相鄰最小點數,遍歷稠密點云中所有點,點云某點搜索半徑內點數少于相鄰最小點數,則認為該點為噪聲點予以排除。
經過異常點濾波之后的點云仍存在部分重建噪聲,來源是位姿估計誤差下重建點的偏移,以遠點代替了真正需要重建的點,遠點往往重建出來帶有特殊的RGB 信息,據此將可疑點濾除。
實驗測試計算機采用Intel CORE i5 10thCPU,Ubuntu 18.04系統。相機為Intel RealsensetmD455 深度相機。實驗流程是使用ROS 發布訂閱話題以獲取相機輸入,采用基于圖優化的G2o 庫,以及PCL 庫做點云處理。采集室外實際植物數據以及TUM 數據集構建測試數據集。
TUM 數據集中挑取任意同一幀圖像分別提取SIFT 特征與ORB 特征。提取出3 431 個SIFT 特征,耗時0.57 s,平均每個特征耗時為166.13 μs。提取出500 個ORB 特征,用時9 978 μs,平均每個特征耗時為19.956 μs。提取速度提高接近一個數量級,這僅是提取單幀圖像得到的提升,當圖像規模高達幾百幀時,將節省更多時間。ORB 特征在同一幀圖像中數量上比SIFT 特征要少,這在特征匹配步驟中能節省時間。
采用自建的室外樹木數據集進行三維重建后進行基于半徑的異常點濾波與基于RGB 的重建噪聲濾波,如圖8 所示,上圖為濾波前,下圖為濾波后,可以看到大部分環境中的噪聲點都被濾除,保留下的點云更加緊湊與清晰。

圖8 基于半徑的異常點濾波與基于RGB的重建噪聲濾波
采用自建室外數據集進行時空重建結果變化檢測,結果如圖9 所示。低矮樹植高度變化明顯,三次不同時間的重建能在同一場景建立起數據關聯,得到高度結果如表1 所示。當收集長期樹木數據集時,也可以進行對應的形態監測,結果如圖10 所示。

表1 低矮植株重建測量

圖9 低矮樹植數據關聯重建結果(上:3月30日;中:4月13日;下:4月27日)

圖10 樹木數據關聯重建結果(上:4月13日;下:4月27日)
該研究具有三方面優勢:改進了數據收集前端,跟蹤匹配ORB 特征而不是SIFT 特征,大大減少跟蹤時的計算量;在后端實現了重建結果的準確時空數據關聯;在重建基礎上盡量還原樹木信息,測量結果誤差在5%以內。后續工作可以通過樹木形態監測結果估算樹木凋落可燃物量。