張 娣,楊 碩

(沈陽化工大學 計算機科學與技術學院,沈陽 110142)

隨著無人機技術的迅速發(fā)展,無人機的應用領域也越來越廣.目前,高效的飛控算法,結合陀螺儀、超聲定位、GPS 定位和計算機視覺算法,使無人機能夠應對較復雜的自然環(huán)境.但是某些特殊情況下,如室內環(huán)境,GPS 失效,無人機精確定位只能依靠陀螺儀和機載攝像頭,陀螺儀僅能保持無人機自身姿態(tài),無法解決無人機位置和朝向漂移的問題.因此,這種情況下,無人機的精確定位就只能依靠計算機視覺算法.但是,計算機視覺定位嚴重依賴于圖像質量,在弱紋理情況下,匹配誤差和置信度較差,使無人機的姿態(tài)參數較難估計,對無人機的控制方面帶來了一定的難度.因此,室內弱紋理環(huán)境下無人機參數的估計非常有必要,準確的參數估計可以更好的控制無人機,使應用場景更多.圖1是無人機參數估計流程圖.

圖1 無人機參數估計流程圖

為了更好的提升無人機的穩(wěn)定性,國內外多名學者進行了研究.王亭亭等[1]采用基于原始ORB 特征的LK 光流法,利用隨機抽樣一致濾波(RANSAC)來提高光流場的精度,再融合里程計、慣導等信息,實現了無人機室內導航,擁有較高的實時性,但其光流速度估計處理較為簡單,采用傳統(tǒng)的前后單追蹤策略,濾除誤匹配效果不理想;齊詠生等[2]針對無人機在停機坪上空飛行時,采集地面視頻,視頻中包含停機坪上的幾何特征,預處理后,采用五步標志提取方法與距離三點法的角點檢測算法相結合,對視覺標志進行提取,提取過后進行匹配,得到單應性矩陣,分解則可得到飛行的參數,該方法針對場景是含有H標志的停機坪,應用場合具有較明顯的特征,在室內面對弱紋理區(qū)域及重復紋理區(qū)域時,沒有明顯的幾何標志,匹配過程中該方法錯誤率較高.周克旻等[3]使用目前流行的基于增強現實(augmented reality,AR)技術的ArUco 標記算法和顏色空間域標記算法,設計了一種多標記的無人機跟蹤系統(tǒng),對無人機在室內飛行過程中進行參數估計.

針對以上部分文獻存在的問題,本文使用的算法在面對室內弱紋理或者重復紋理區(qū)域時,可以準確的估計出無人機飛行時的平移參數及旋轉參數.本方法在面對室內弱紋理區(qū)域時,使用特征提取較快的ORB光流法,配合前后雙向雙追蹤策略和RANSAC 算法,提高參數估計的準確性,降低了誤匹配概率;在旋轉參數估計時,提取單應性矩陣,分解進行估計,得到準確的參數.

1 本文算法

本文主要針對的是室內弱紋理及重復紋理區(qū)域采集的視頻,所以算法主要分成以下3 部分:視頻預處理、平移參數估計和旋轉參數估計,流程圖如圖2所示.算法實現步驟如下.

圖2 本文算法流程圖

步驟1.預處理階段:在處理視頻幀之前,首先將彩色視頻幀轉換成灰度圖,進行后續(xù)處理,節(jié)省時間.拍攝過程中難免存在一定的噪聲,為了不對后續(xù)處理造成影響,先對采集的視頻幀進行去噪處理.

步驟2.平移參數估計階段:根據室內弱紋理區(qū)域可能存在其他物體的情況,對采集視頻圖像進行判斷,根據不同環(huán)境采用不同的算法進行估計.

步驟3.旋轉參數估計階段:通過ORB 算法,提取出4 個特征點,如地板的角點等,結合OpenCV 圖像庫可以很好的計算出一個擬合所有點的單應性矩陣,分解得到旋轉分量

2 預處理

無人機室內飛行過程中,面對地板磚等弱紋理區(qū)域,會采集到該區(qū)域的彩色視頻,首先轉成灰度圖像,消除算法對彩色圖像的依賴,減少計算量,降低程序運行時間.

無人機飛行過程中,由于飛行會使得采集的視頻模糊,具有部分噪聲,本算法采用雙邊濾波的方法去除噪聲.均值濾波[4]直接對鄰域內像素值取平均,對圖像進行模糊,簡單快速,飛行過程本身會產生一部分的模糊,若使用均值濾波會使得視頻幀更加模糊;高斯濾波[4]是一種適用性較強的去噪算法,相比均值濾波,高斯濾波考慮到鄰域內像素點和距離的關系,即距離中心像素越近證明和該像素點越相似,權值對應高一點,這樣效果較好.雙邊濾波[5],結果實際上由他鄰域內的像素點決定.類似于均值濾波后的結果是根據相鄰像素值的均值取得的,只不過各像素值權重為1.雙邊濾波的權重則是根據相鄰像素的距離和相似度得到的.所以,雙邊濾波相比高斯濾波和均值濾波,既能對像素平滑也能保持一定的細節(jié)信息.室內采集視頻時,主要面對地板,水泥地面等弱紋理區(qū)域,會包含部分縫隙,這種區(qū)域在我們后期參數估計時是非常好的特征,雙邊濾波可以在去噪的基礎上保存下邊緣信息,但是由于雙邊濾波計算時間較長,不能保證實時性,因此選用快速雙邊濾波[5],對視頻幀進行去噪,經過多次實驗,x方向參數取10,y方向參數取30,得到的濾波效果是最好的,即保留了邊緣,又去除了噪聲.

3 參數估計

無人機飛行過程中,參數估計主要包含兩方面:平移參數估計和旋轉參數估計.針對平移參數估計,首先判斷視頻類型,是均為弱紋理的室內地面,還是存在非弱紋理區(qū)域,根據區(qū)域不同的特征選擇不同的算法進行估計.

小型無人機室內飛行時,攝像頭采集地面的視頻,面對室內的環(huán)境時,存在部分地板和大理石地面,這些區(qū)域屬于重復紋理區(qū)域,本文將該區(qū)域劃分為弱紋理區(qū)域.室內地面大都是紋理均勻的地板,但是偶爾會存在其他區(qū)域,如圖3所示.針對不同種類的圖像,運動估計的方法也有所差別,所以在估計之前,應先判斷視頻種類,根據視頻選擇對應的方法,節(jié)約時間,加快參數估計速度.因為含有非弱紋理區(qū)域的圖像與均為弱紋理區(qū)域的圖像,所以需要對所屬區(qū)域進行判斷,本方法使用灰度共生矩陣(GLCM)[6]中的特征值進行判斷.由于灰度共生矩陣計算量較大,所以應在不影響圖像紋理特征的前提下,對灰度進行壓縮,壓縮成16 個灰度級,以1為步距、8×8的滑動窗口,采集0°,45°,90°,135°四個方向的灰度共生矩陣,計算每個矩陣的熵和反差分矩陣特征值,對所屬區(qū)域進行判斷,再選擇后續(xù)算法進行參數估計.

圖3 弱紋理區(qū)域及存在非弱紋理區(qū)域場景

3.1 平移參數估計

通過以上對視頻種類的判斷,存在兩種類型視頻,針對這兩種類型的視頻,我們采用兩種方法進行運動矢量估計.當視頻均為弱紋理區(qū)域時,采用ORB 光流特征算法;當圖像存在其他特征時,采用改進Snake 模型,將對應區(qū)域圈出來,計算重心,進而估計出圖像的運動矢量.

3.1.1 非弱紋理區(qū)域參數估計

非弱紋理區(qū)域通俗來說就是弱紋理區(qū)域包含一些其他東西,如室內地面上會存在書本、凳子等物體,將弱紋理的地板進行了遮掩,此時選用Snake 算法,將該區(qū)域分割出來,進而計算重心偏移得出此時的無人機飛行參數.

(1)Snake 模型

實際上,圖像分割有大量算法,例如閾值法、區(qū)域分割法、邊緣分割法等[7].最簡單的閾值法將圖像根據閾值分割為二值圖像,忽略了圖像中的紋理信息,只關注圖像本身像素的值,因此在很多時候效果差強人意;區(qū)域分割方法需要給定一個種子點,進行區(qū)域分割,但是會造成區(qū)域不連續(xù)的效果;邊緣分割法面對的圖像需要是特征豐富的,邊緣明顯的,這樣分割出的效果才是理想的.上述幾種圖像分割方法對噪聲敏感,無人機采集視頻會存在部分噪聲,因此幾種分割方法都會產生不理想的分割效果.

Snake 模型屬于可變模型,使用閉合的參數和曲線對目標區(qū)域進行描繪[7].簡單來講,Snake 模型就是一條可變形的參數曲線及相應的能量函數,以最小化能量目標函數為目標,控制參數曲線變形,具有最小能量的閉合曲線就是目標輪廓[7].Snake 模型由一組節(jié)點組成,如式(1)所示,x(s),y(s)代表圖像中作為控制點的坐標,s是描述邊界的自變量:

在Snakes的控制點上定義能量函數如式(2)所示(反映能量與輪廓之間的關系):

其中,α,β參數理和口操縱著模型的物理行為和局部連續(xù)性.根據Snake 模型的函數,可以知道,想要得到目標的輪廓,需要通過能量函數動態(tài)化,進行多次迭代,逐漸逼近,最終得到理想的輪廓.Snake 模型的優(yōu)點就是不論畫面質量如何,都可以得到封閉、光滑的目標輪廓,相對于其他圖像分割算法,更適合處理差異性較大、結構復雜的圖像.

但是相對較費時,本文使用改進的Snake 算法——分段迭代B-Snake 算法[6].改進算法的原理是通過一次迭代找到目標的曲線單元,并進行分離,在下一次迭代中,只對逼近曲線單元進行迭代;當所有曲線迭代結果接近目標時,對各部分曲線組合生成目標區(qū)域的輪廓.本次采用文獻[8]的改進分段迭代B-Snake 算法,進行5 次分段迭代,可以得到理想的輪廓,如圖4所示,由于現實生活中會存在許多復雜輪廓的物體,因此本次選用440×350的輪廓較不規(guī)則的海星圖像進行實驗,提取出圖像中海星的輪廓特征.改進算法快速收斂到了目標輪廓,節(jié)約時間,減少了迭代次數,提升速度.

圖4 星形圖案及Snake 算法效果圖

(2)重心計算

根據上述Snake 模型算法,計算兩幀圖像的區(qū)域重心,通過重心的比較可以得出此時無人機的平移分量.如圖5所示,為前后兩幀圖像中重心所在區(qū)域截圖,數字代表該位置坐標,陰影部分即為經過計算的重心點,假設上一幀重心點坐標為(34,36),下一幀為(34,37),采用式(3)對坐標相減,即可得到此時無人機位移量,即1 cm.

圖5 圖像重心舉例圖

其中,H1(x1,y1)和H2(x2,y2)為前后兩幀視頻圖像Snake算法圈出的輪廓區(qū)域的重心,H為兩幀重心的偏移量.

3.1.2 弱紋理區(qū)域參數估計

室內地板多為重復紋理和弱紋理區(qū)域,特征不明顯,使用起初的特征匹配產生過多的誤匹配,導致匹配后的參數估計產生大量錯誤,因此不可使用多數文獻中使用的特征匹配估計參數的方法[9].本文在ORB 光流算法的基礎上,結合無人機特定時刻的飛行速度,使用一種適用于無人機面對室內弱紋理地板區(qū)域時的細分匹配規(guī)則算法—改進ORB-LK 算法,對平移參數進行估計.

(1)ORB 特征提取

ORB 特征描述具有較高的實時性,是由描述子和關鍵點兩部分組成.改進的FAST 角點作為關鍵點,BRIEF[9]作為描述子.首先將采集到的視頻幀,建立圖像金字塔,采集每層金字塔中的FAST 角點特征,由于FAST 特征不具備旋轉不變性,本文對FAST 特征添加方向因子;然后,利用BRIEF 算子對特征點進行描述.

(2)LK 光流算法

LK (lucas-kanade)光流算法[10]是利用兩幀之差來追蹤特征點瞬時速度的估計方法,通過迭代計算可以得到一種稀疏光流場.但是LK 光流法必須滿足3 條前提假設[10]:(1)相鄰幀之間的亮度、顏色不發(fā)生變化;(2)相鄰幀之間的像素點運動是足夠小且連續(xù)的;(3)相鄰視頻幀之間,像素的運動一致,擁有相同的光流.因此,可得亮度守恒約束方程[10]如式(4):

其中,vx和vy分別表示圖像在x和y方向上對應的光流值;Ix和Iy表示像素點在x和y方向上對應的偏導數;It表示圖像像素點亮度關于時間的導數.

(3)改進ORB-LK 算法

無人機飛行過程中,光流傳感器開始啟用[11],首先進行特征提取,在接下來計算光流值.當面向室內弱紋理區(qū)域以及重復紋理區(qū)域時,存在部分區(qū)域含有相似的紋理或者特征較少,在參數估計過程中在匹配過程中若使用暴力匹配法這種匹配算法,會造成大量的錯誤匹配,使估計的參數不準確,因此本次選用FLANN匹配算法.

在原始ORB 光流算法[12]基礎上,結合無人機飛行的速度vu,可以計算出光流的閾值如式(5)所示:

其中,h為無人機的飛行高度,根據具體情況而定,假設為一個定值,f表示攝像頭焦距,根據經驗,光流閾值為5.當光流閾值大于 Δ 時,使用ORB 提取特征,使用雙向單追蹤策略對特征點進行篩選,稱為初步匹配;當光流閾值小于 Δ 時,使用FLANN 完成特征匹配,匹配過程中可能會存在誤匹配情況,使用RANSAC 進行剔除,進行二次匹配.初步匹配的目的在于找到視頻幀中盡可能多的特征點,為二次匹配做準備,提高二次匹配的速度;二次匹配目的是優(yōu)化特征匹配,可以得到準確的參數估計,從而對無人機進行精確的定位.(1)初步匹配:首先對圖像使用ORB 算法進行特征提取,得到特征點集M,再使用金字塔LK 算法找到下一幀的圖像特征點集N(前向跟蹤),然后利用金字塔將特征點集N找到前一幀對應的點集O(后向跟蹤),由于噪聲原因,使得M和O的點集不能完全重合,存在一定的誤差距離,遍歷M和O點集,計算距離的誤差p,選取距離小于σ的特征點集,再結合點集N,得出篩選后的點集MiNi+1.(2)二次匹配:當初步匹配過程中光流閾值小于 Δ 時,進行二次匹配.根據初步匹配點集MiNi+1可以得到第i幀的特征點集Ai,ORB 特征繼續(xù)提取第i+1 幀對應的特征點集Ai.然后使用FLANN 進行匹配,匹配會存在誤匹配,此時使用RANSAC 算法剔除誤匹配點.匹配效果如圖6所示.

圖6 FLANN 匹配效果圖

根據上述的結果,再采用前后雙向追蹤策略,如圖7所示.使用第i幀的特征點集Ai進行前向追蹤,得點集Fi+1,再接下來后向追蹤,獲得點集Bi,同樣由于噪聲的存在,導致Ai和Bi不能完全重合,存在一定的誤差,距離誤差記為di;同理,可以得到Ai+1的前向跟蹤點集Fi與后向跟蹤點集Bi+1,同樣可以得到距離誤差di+1.光流法存在3 個假設,所以可以得第i幀點集和Fi點集距離誤差為dj;第i+1 幀,Ai+1點集和Fi+1點集距離誤差為dj+1.經過實驗,當誤差閾值σ取10 時,效果最好,進而得出匹配過濾函數W,如式(6)所示:

圖7 前后雙向追蹤策略

3.2 平移參數估計算法描述

算法1.平移參數估計輸入:拍攝視頻S,前后相鄰幀像素,輸出:平移參數H H1(x1,y1) H2(x2,y2)1 for each frame in S 2 {計算每幀GLCM的熵和反差特征值m,n 3 if(m>3,n<15)4 {由式(1)和式(2)通過B-Snake,得到的非紋理區(qū)域輪廓5 由式(3)計算前后兩幀輪廓重心偏移量H 6 else 7 {提取ORB 特征值8 由式(5)計算光流值Δ 9 if(Δ>5)10 {通過光流與慣導系統(tǒng)得到參數}11 else 12 {FLANN 特征匹配13 RANSAC 剔除錯誤匹配點14 由式(5)計算光流值Δ 15 通過光流與慣導系統(tǒng)得到參數}}}}

3.3 旋轉參數估計

無人機在飛行過程中,不僅會平移產生平移參數,同時也會存在一定的旋轉運動產生旋轉參數,所以參數估計過程中,對旋轉參數也要進行估計[13].單目攝像機拍攝視頻對前后兩幀特征點進行匹配,到處旋轉矩陣,即可得到旋轉參數.在飛行過程中,主要旋轉角度分為繞Z 軸的航向角(ψ),繞X 軸的俯仰角(Ф)和繞Y 軸的翻滾角(θ)[14],如圖8所示.

圖8 無人機旋轉角度示意圖

根據特征匹配,可以得到前后兩幀圖像間的單應性矩陣H,對矩陣H進行分解即可得到旋轉參數,如式(7)所示:

其中,Mn+1和Mn分別為第n+1 幀與第n幀圖像的像素值,矩陣H中包含旋轉矩陣R(rx,ry,rz).根據實驗數據處理結果,對旋轉矩陣進行解算,由于數據誤差,旋轉矩陣不能精準滿足R=RRT=RTR=I.因此對旋轉矩陣進行SVD 分解[14],令R=UDVT,D=I,使旋轉矩陣R強制轉為滿足條件的旋轉矩陣R1=UIVT.特征點在機體坐標系和世界坐標系下對應關系如式(8):

3.4 旋轉參數估計算法描述

算法2.旋轉參數估計輸入:拍攝視頻S輸出:旋轉參數Ф,θ,ψ 1 for each frame in S 2 {提取ORB 特征值3 由式(7)計算單應性矩陣H 4 對單應性矩陣進行分解,得到旋轉矩陣R 5 對R 進行SVD 分解6 由式(8)進行世界坐標系與相機坐標系的轉換

?

4 實驗結果

根據上述算法,本實驗使用分辨率為440×350,焦距為50 mm的攝像機作為圖像采集設備,使用電腦(CPU為Intel I7 3520M,主頻達到2.9 GHz,內存為4 GB)作為圖像處理計算機,結合OpenCV 3.3.1和Ubuntu 16.04環(huán)境下進行仿真.

針對室內弱紋理條件下,選擇了多個場景進行實驗.為了驗證算法的準確性,在飛行器平穩(wěn)飛行過程中,對無人機制造一定的移動,通過算法進行估計,在平移和旋轉參數已知情況下,進行參數估計,將估計量與偏移量進行比較.針對不同的偏移量進行估計,得到如表1–表3的結果.

由表1–表3可以證明,本方法在無人機飛行過程中,平均處理時間達到了25 ms,準確率為95%.文獻[3]提出一種多跟蹤的無人機姿態(tài)參數估計系統(tǒng),在參數估計方向屬于較新的一種方法,本文相較于該方法,運行時間平均縮短了10 ms,準確率有大幅度提升;文獻[14]采用改進的ORB 光流算法進行參數估計,本文在該基礎上進行進一步改進,增加了基本的去噪方法,使得準確率有大幅度的提升.

表1 室內不同地面平移量估計情況比較

表2 室內不同地面旋轉參數估計比較

表3 各方法實驗結果對比

5 結語

通過圖像的預處理,對采集的視頻圖像進行去噪,結合圖像的灰度共生矩陣對視頻種類進行判斷,根據視頻紋理重復類型選擇使用改進ORB 光流算法進行平移參數估計,或者是使用Snake 模型進行估計,聯(lián)合重心偏移計算得出運動平移分量.本文提出的方法結合灰度共生矩陣進行判斷,兩種方法相結合,避免了使用單一方法,全面匹配,估計平移分量的問題.結合單應性矩陣,分解得到旋轉參數.本方法對多種室內弱紋理或重復紋理區(qū)域進行參數估計,擴大了處理范圍,降低了處理時間.

本文算法在精確估計出無人機移動參數情況下,可結合飛控系統(tǒng),將估計出的參數傳入系統(tǒng),在室內GPS 信號較弱情況下,對無人機進行穩(wěn)定控制,保證室內采集到穩(wěn)定清晰的圖像材料.