楊佳琪, 范大昭, 楊佳賓,, 楊幸彬, 紀 松
(1. 信息工程大學地理空間信息學院,河南 鄭州 450001; 2. 北京字跳網絡技術有限公司,北京 100086)
無人機遙感是近年來興起的一種新型遙感系統,具有快速、靈活、低成本、高分辨率等特點[1]。傳統雷達(light detection and ranging, LiDAR)點云由于其成本較高,難以得到廣泛應用。隨著航攝影像質量的不斷提高,利用無人機影像處理軟件[2]獲取三維信息已成為一種新型測繪手段。然而由于無人機影像的數據量大、軟件處理算法復雜,數據處理常依賴于離線工作流[3],難以在線獲取所需三維信息,因此在自然災害的快速應急響應中難以發揮良好的作用。
近年來,利用同步定位與地圖構建(simultaneous localization and mapping,SLAM)技術在線獲取無人機影像位姿普及[4]。相較于LiDAR點云,影像密集匹配點云具有豐富的紋理特征信息,更易于高效處理[5],且成本相對較低。為了實現高質量的影像密集匹配,文獻[6]提出一種改進代價的半全局匹配算法,以獲取視差圖,通過限制視差范圍減少了計算資源。文獻[7]采用雙目視覺方法實時獲取三維點云,實現了對遠距離目標三維信息的實時獲取。文獻[8]提出一種基于紋理優化的動態約束視差半全局匹配(semi-global matching,SGM)算法,通過對半全局立體匹配代價積聚過程的優化,大大減少了深度估計過程的耗時。文獻[9]提出一種基于多視影像間的重疊關系與方向確定高分辨率航空多視影像密集匹配算法,完成了密集點云的生成與融合。物體邊緣往往具有較大的視差變化,文獻[10]提出了一種顧及影像紋理信息的視差圖膨脹腐蝕算法,結合金字塔影像進一步約束視差搜索范圍,減少了算法所需的內存開銷和計算時間。上述幾種方法盡管在一定程度上優化了算法耗時和效果,然而僅能考慮兩視圖信息,在應急響應中難以滿足需求。文獻[11]針對無人機短時間應對自然災害設計了一個在線三維重建框架,使用開源的塊匹配算法逐步增量生成密集的地理參考三維點云和數字表面模型(digital surface model, DSM),但在影像分辨率過高時無法實時運行。文獻[12]設計了一種實時測圖框架,基于視差平面的快速立體匹配算法實現了影像密集匹配過程,然而獲取的深度和點云信息精度有限。為了最大限度減少計算時間、提高計算精度,文獻[13]提出了一種基于物方的半全局匹配方法直接獲取DSM,但在計算效率和精度方面優化較少。文獻[14]提出了一種基于半全局鉛垂線軌跡匹配的直接DSM生成方法,可以考慮多視信息在物方空間中直接生成DSM,但無法實現在線計算。
為了解決上述問題,本文提出一種無人機在線密集點云和DSM生成算法,結合增量式的多視圖半全局在線密集匹配算法獲取實時點云與DSM。首先,利用SLAM和實時差分定位(real time kinematic,RTK)松耦合方式在線估計影像位姿,在此基礎上選擇當前幀的多視圖參考幀。然后,在深度空間計算離散采樣深度,將采樣深度投影到多視圖參考幀上計算聯合匹配代價,并考慮多視圖遮擋關系對代價進行優化。最后,基于增量式半全局積聚策略得到積聚代價體,利用“勝者為王”算法結合拋物線擬合得到最優深度圖,利用深度圖融合得到密集點云和DSM。利用3組典型地區的無人機影像對本文方法進行在線驗證,基于定性和定量評價驗證本文方法生成密集點云和DSM的有效性。
基于無人機影像在線生成密集點云和DSM的流程主要包括:影像在線位姿估計、影像在線密集匹配獲取深度圖、深度圖融合和DSM生成4部分。為了在線估計影像位姿,在ORB SLAM3[15]系統框架基礎上進行了改進。無人機系統上搭載的RTK可實時獲取無人機位置和軌跡,通過對視覺SLAM和RTK獲取的位置信息進行松耦合綁定,從而在線獲取地理坐標框架下的無人機影像位姿信息。
針對僅利用兩視圖信息可能會因視角變化存在遮擋變形和直接將場景深度作為深度空間搜索范圍進行求解時存在較大困難的問題,本文提出一種增量式的在線影像密集匹配方法,利用多視圖信息直接在深度空間計算聯合代價。另外,為節省在線深度估計計算時間,增量式地將前一幀深度映射至待求解幀,增量式約束代價體計算和積聚過程。具體算法如下:

(2)多視圖聯合代價體計算。為了計算當前幀深度圖,直接在深度空間進行解算,基本思路為利用圖像位姿將不同深度投影至多幀圖像計算匹配代價。首先根據無人機飛行時相對航高H和相機大致傾角φ,利用z=Harcsinφ確定深度空間搜索范圍[zmin,zmax];然后利用式(1)將最小和最大深度投影至多個參考幀圖像,計算兩個投影點之間像素距離的最大值。
(1)
式中,zl為給定深度值;p為當前圖像像素點齊次坐標,p′l為像點p在參考幀圖像的投影點坐標;K、R、T分別為相機內參和相機到世界坐標系的旋轉和平移矩陣。為在深度空間求解,本文在深度空間內進行非等間隔采樣,離散采樣深度值zl(l=0, 1,…,L)計算公式如下
(2)
按照式(2)將像點p按照不同zl投影至參考幀圖像得到極線上像點坐標p′l,后續過程中在D軸采樣方向上將每個采樣深度位置對應一個匹配代價值,計算過程如圖1所示。

圖1 多視圖聯合代價體計算過程
為了計算像點p在不同深度下的匹配代價,本文對文獻[16]的census代價進行改進,采用中心對稱思想避免中心像素點噪聲對代價的影響,同時節省了一半的計算耗時。考慮視圖間像素遮擋概率,聯合多視圖信息計算代價權重,通過權重聯合匹配代價計算代價矩陣C。中心對稱census比特鏈的計算函數式為
P(u+i,v+j)]
(3)

(4)
式中,WCTfj(u′,v′)為P(u,v)所對應投影點的中心像素灰度P′(u′,v′)的比特值。
(3) 增量式代價積聚與深度計算(考慮深度圖映射)。考慮到Census局部代價本身的約束,為了獲得更高的深度圖精度與完整度,結合半全局匹配方法中掃描線優化策略,對代價矩陣C進行全局能量的聚集,使代價矩陣更加精確可靠。鑒于本文方法在深度空間上直接計算多視圖聯合代價,因此將傳統全局能量函數E(l)[17]修改為
(5)
式中,lp為像點p的深度采樣位置;lq為像點p臨近像點q的深度采樣位置;c(p,lp)為像點p在lp處的多視圖聯合匹配代價;T()為判斷函數,其目的是對像點p臨近視差差值等于1和大于1的像點q代價分別增加P1、P2懲罰項。

(6)
在代價積聚過程中,對fi幀圖像上有先驗深度的區域按照一定的容差范圍對采樣位置進行代價積聚;對于無先驗深度的區域,按照正常積聚策略完成代價積聚,具體計算過程為
(7)

以上計算僅考慮了一個方向r的代價聚合值Lr(p,lp),聚合后的信息量較少,因此考慮從8個方向對代價矩陣進行代價積聚(如圖2所示),最終生成積聚后的代價矩陣。

(8)

圖3顯示了利用兩視圖半全局立體匹配方法[17](簡稱為兩視圖法)和本文方法生成的深度圖和點云。可以看出,本文方法相較于兩視圖法更好地保護了深度圖的邊緣與棱角信息(如圖3(a)所示);此外,在弱紋理處(如圖3(b)所示)能夠較大程度上改善代價計算的準確度,從而使深度圖更加平滑精確。
為了在線獲取測區密集點云,首先將前一幀fi-1深度投影至當前幀fi得到先驗深度;然后將fi估計的深度與先驗深度作比較,剔除相對誤差大于5%閾值范圍的噪聲,對剔除后的深度進行等權融合得到新的深度圖;最后按照式(8)將新的深度圖增量式地投影至三維空間,即可得到稠密點云結果。
(9)

試驗選取鄉村、山地、城區3組典型無人機影像測試本文算法生成稠密點云和DSM的效果。鄉村區域包含弱紋理建筑物和邊緣輪廓清晰的伐木場,用于評定本文方法在弱紋理和物體邊緣的匹配效果;山地區域地形起伏較大,不同視角拍攝的圖像存在一定變形和遮擋區域,用于評價本文方法在高差大的山地區域點云生成的穩健性;城區數據用于評定本文方法在城區的效果。3個測區拍攝的影像位置(SLAM在線位姿估計得到)和本文方法生成的稠密點云如圖4所示。

圖4 3個測區拍攝時相機位置和稠密點云結果
本文算法在Windows 10 64位系統下基于VC++2015開發,其中影像位姿估計過程結合了ORB-SLAM3在線獲取,影像密集匹配、增量式點云融合與DSM生成過程為本文核心算法,密集匹配和點云融合過程使用OpenCL進行了并行加速。試驗使用的硬件平臺為華碩筆記本,處理器為Intel? CoreTMi7-8700 CPU、最高主頻為3.2 GHz,CPU內存為15.6 GB,圖形顯卡為NVIDIA GTX1060。
為驗證本文方法生成的稠密點云效果,在3組數據的基礎上選取局部典型地物稠密點云進行放大顯示,如圖5—圖7所示。由圖5可知,無論是伐木場還是弱紋理房屋屋頂的三維點云都具有較好的精度和完整度;由圖6可知,即使是在地形突變的區域(圖6中②③)點云形狀和走勢也符合實際地表形態,說明本文方法能夠在一定程度上抗拒視角和遮擋影響,穩健性較好;由圖7可知,點云完整度較高,由于影像是從視頻中抽幀得到,因此點云精度相比于鄉村和山地略低,但從整體上看本文方法生成的點云在建筑物和樹木上具有較好的完整度和精度,符合實際地物形態。

圖5 鄉村區域稠密點云結果

圖6 山地區域稠密點云結果

圖7 城區稠密點云結果
圖8為本文增量式DSM生成的效果圖,可以看出無論是在弱紋理房屋的鄉村區域還是地面起伏較大的山地區域,本文方法生成的DSM均能保持完整地表形態,且平滑性較好;從建筑物的邊緣細節和地形斷面處可以看出,本文方法生成的DSM具有清晰的輪廓。

圖8 3個測區DSM
為進一步定量評價本文方法生成的稠密點云絕對精度,導入3組數據的位姿數據到SURE軟件,離線生成稠密點云且與本文生成結果進行疊加,通過絕對值誤差評定本文方法精度。圖9為3組數據的誤差分布。

圖9 3組數據誤差分布
由圖9可以看出,本文方法生成的稠密點云在大部分區域內誤差較小,鄉村地區的誤差最小,相比于SURE軟件生成的稠密點云,誤差小于3 cm的像素數占比超過90%。
對圖9誤差進行統計分析,剔除不符合正態分布3σ的誤差后,計算一定范圍內最大絕對誤差、絕對值均值誤差和均方根誤差來定量評定本文方法生成的點云精度,結果見表1。由均方根誤差可知,本文方法生成的點云在鄉村地區較小,精度最高,在山區和密集建筑物的城區誤差相對較大;由最大絕對誤差和絕對值均值誤差可知,受視頻影像質量影響城區最大絕對誤差和絕對值均值誤差最大,山地次之。

表1 3個測區的點云精度評定結果 m
為定量評價本文算法的性能,將本文方法計算的位姿結果導入商業軟件SURE和開源軟件COLMAP后離線生成稠密點云和DSM,然后分別統計兩種方法與本文方法的計算耗時,結果見表2,需要說明的是,3種方法均在相同的硬件環境下測試,使用COLMAP[19]時關閉了幾何一致性深度計算。

表2 算法的計算性能對比結果
由表2統計結果可知,本文方法平均一幀計算耗時僅約為0.4 s,具有較高的計算效率。由于本文方法和SURE軟件(密集匹配采用類似半全局立體匹配的思想)的輸入輸出數據相同,因此具有較好的可比性。由3組數據的對比結果可知,本文方法在精度沒有損失太多的情況下計算效率要遠高于SURE軟件(計算效率是SURE軟件的約25倍),平均計算效率基本可以達到2~3幀/s,足以滿足無人機在線生成稠密點云和DSM的要求。
傳統無人機影像獲取密集點云和DSM的方法耗時過多、點云質量和處理速度難以有較好的平衡,無法實現對自然災害的應急響應。本文針對此問題提出了一種無人機大場景在線密集點云和DSM生成方法。通過直接在深度空間內引入多視圖中心對稱census代價,并考慮視圖遮擋關系計算聯合代價體,減少一半計算耗時的同時提升了代價準確性,在此基礎上提出增量式代價積聚策略,聯系當前幀與歷史幀深度進一步優化耗時且提高深度圖精度,將相鄰幀深度進行檢查并剔除噪聲后,快速在線生成密集點云與DSM。選擇典型的鄉村、山地、城區3組測試數據進行試驗。結果表明,本文方法能夠以2~3幀/s的平均處理速度在線生成測區稠密點云和DSM信息,同時生成的點云精度與SURE軟件相比,可控制在為5 cm誤差范圍內,DSM上地物邊緣清晰、平滑度較好,可以為自然災害的快速應急響應提供精確實時的三維信息。