關震宇，李杰，楊歡，徐蓓蓓，劉暢

(北京理工大學 機電學院，北京100081)

0 引言

微小型無人飛行器(MAV)并不是常規(guī)飛行器的簡單縮小，由于受到體積和載荷重量等方面的限制，其在結構設計、氣動設計、動力系統(tǒng)和導航控制方法等方面與傳統(tǒng)飛行器存在較多不同。尺寸的微小化對飛行器姿態(tài)測量系統(tǒng)的設計提出了挑戰(zhàn)，傳統(tǒng)的高精度、低漂移的陀螺由于體積和重量的原因并不適合微小型飛行器使用，而基于微機電技術開發(fā)的微陀螺和微加速度計等器件受限于現(xiàn)有的精度、穩(wěn)定性、溫漂等因素制約，在長期使用時會產生較大的積累誤差［1］。由于上述傳統(tǒng)器件的局限性，需要探索一種新型的微小型飛行器姿態(tài)信息估計方法對傳統(tǒng)方法進行替代和補充。

早期飛機并沒有復雜的傳感系統(tǒng)，飛行員對于飛機姿態(tài)的獲取直觀地來源于其視野中地平線位置的變化;由于任務使命的關系，MAV 通常攜帶有攝像設備進行偵察和航拍作業(yè)，在偵察相機的視場中，地平線又是自然界最明顯的標志物，因此借助地平線位置來測量飛行器姿態(tài)就成為了一種可行的思路。目前關于地平線的各種提取方法已經相當成熟，可以分為基于邊緣特征的地平線提?。?-3］、基于區(qū)域統(tǒng)計特征的地平線提?。?-4］和基于機器學習的地平線提取［5］;在使用地平線信息進行姿態(tài)角計算方面，同樣開展了大量的研究工作。參考文獻［6 -7］通過地平線與圖像水平線夾角進行飛行器滾轉角估計，參考文獻［8 -9］使用分割后天地各自的質心位置來估算滾轉角，并使用天地分割線兩側面積的變化來估計飛行器的俯仰角，參考文獻［3］也采用了類似的方法來估計俯仰角，但天地線兩側面積比例會由于滾轉角的變化而發(fā)生較大的變動，因此會對俯仰角的估計帶來不小的估計誤差，參考文獻［10］具體指出了這一點問題，并提出上述方法甚至可能會造成嚴重的飛行事故。參考文獻［11］使用了攝像機模型對俯仰角和滾轉角的計算進行了理論推導，提出了一種基于直線模型的俯仰角和滾轉角的計算方法，但不足之處在于該方法需要預先對所使用攝像機進行標定獲取其內參數矩陣，而且從其仿真結果來看該方法測量精度不高。除此之外，上述所有方法都只能計算飛行器的俯仰角和滾轉角，對于飛行器俯仰角速度、滾轉角速度和偏航角速度等重要姿態(tài)參數或無法直接測量，或需要通過差分運算進行估計，在實際使用中存在一定的局限性。參考文獻［12］針對上述局限提出了一種基于光流的滾轉角速度和俯仰角速度提取方法，但該方法需要計算圖像全局光流場，由于全局光流場的計算量較大，因此該方法在實時應用方面存在問題。參考文獻［13］總結了基于視覺的無人飛行器姿態(tài)估計方法，指出基于光流的估計方法可以實現(xiàn)俯仰、滾轉和偏航3 個通道姿態(tài)信息的估計。

為了克服上述不足，本文提出了一種圖像中的直線稀疏光流場計算方法，并在此基礎上提出了基于直線稀疏光流場的飛行器姿態(tài)信息估計方法。該方法可以用于無法安裝微機電陀螺的微小型飛行器進行姿態(tài)信息獲取，也可以通過與傳統(tǒng)導航手段進行融合，輔助飛行器進行導航控制。該方法不依賴計算圖像的全局光流，僅需要計算地平線附近的稀疏光流，從而降低了算法計算量。經過理論推導與仿真分析，該方法可以有效地對飛行器俯仰角速度、滾轉角速度和偏航角速度等姿態(tài)信息進行估計，并在飛行器角速度較低的情況下具有較高的精度。

1 直線稀疏光流場的定義與提取算法

光流的概念最早由Gibson 于1950年提出，是指圖像中模式運動的速度。光流場是一種二維(2D)的瞬時速度場，其中2D 速度矢量是景物中可見的三維速度矢量在成像表面的投影［14］。光流計算是由Horn 等［14］首先提出，之后出現(xiàn)了一系列的有關全局光流計算方法［15-18］。全局光流理論上可以估計出圖像上每一個像素點的光流矢量值，其不足之處在于需要大量的迭代計算來獲取，算法時間效率不高，且易受到圖像局部噪聲的影響。為了解決上述問題，稀疏光流場的概念應運而生。目前文獻中常涉及的稀疏光流場主要指角點的稀疏光流場，其獲取方法首先使用基于特征的方法對相鄰兩幅圖像進行角點匹配，隨后使用匹配角點的運動矢量來近似替代其光流矢量。這種計算方法相較于全局光流計算具有運算速度快、對噪聲不敏感等優(yōu)點，其缺點是僅能獲得圖像中局部(角點處)的光流信息，但因為角點作為圖像中重要的特征，它的局部光流也含有重要信息，因此稀疏光流場在圖像識別、運動檢測、目標提取等方面具有廣泛應用。

所謂直線稀疏光流場，就是在角點稀疏光流場基礎上提出的一種圖像稀疏光流場，它是指相鄰兩幀圖像中共線點的光流矢量的集合，如圖1所示。與角點稀疏光流場不同的是，直線稀疏光流場是通過提取兩幀圖像中相匹配的直線，并通過計算直線上每一個點的運動矢量而獲得的。因為圖像中的邊緣信息常用直線描述，因此直線稀疏光流場就特別適合描述圖像中邊緣的變化。

圖1 直線稀疏光流場定義Fig.1 Sparse line-optical flow field

直線稀疏光流場的計算方法如下:

設攝像機幾何模型為透視模型，建立OXYZ 空間坐標系，oxy 為圖像坐標系，o 為攝像機光心，光軸與Z 軸重合，圖像坐標系原點設為在OXY 坐標系(0，0，f)處，其中f 為焦距，x 軸和y 軸分別與X 軸和Y 軸平行，如圖2所示。

在成像平面上，Li+1和Li是空間一條直線在成像平面的連續(xù)兩幀投影，p 是空間一點P 在成像平面上的投影。對于直線Li上一點pi，其在直線Li+1對應點計為pi+1.

圖2 攝像機透視投影模型Fig.2 Camera perspective projection model

不失一般性，設f=1，對于該攝像機模型有

對空間一點P=(X，Y，Z)T，其在成像平面投影p=(x，y)T，點P 在空間的運動滿足如下運動方程:

式中:矩陣R 和矩陣T 分別描述點P 的旋轉運動和平移運動;Δt 為兩幀圖像間的時間間隔，當Δt 很小時，可對(2)式作如下近似:

式中:v = (v1，v2，v3)T，ωz分別表征P 點的平移速度矢量和繞光軸逆時針方向的旋轉角速度矩陣，將其代入(2)式，得

經整理，即可得

展開(5)式最后一行，得

則

當點P 沿光軸的運動相比于點P 到光心的距離非常小時，即滿足

時，有

此時，對(5)式兩邊同除Z(t)，

經整理，即可得

則有

式中:w(x，y，t)表示圖像上p 點在t 時刻的二維運動速度，即時變圖像的光流，這樣就建立了圖像上點位移場與光流場之間的關系。

在圖像坐標系下，設直線Li和Li+1的直線方程分別為

式中:ki、bi、ki+1、bi+1分別為直線Li、Li+1的斜截式直線方程參數。點pi和點pi+1的坐標分別為(xi，yi)，(xi+1，yi+1).

點pi處圖像灰度滿足光流恒等式為

式中:Ix、Iy、It分別是圖像上(xi，yi)點處的灰度值在x、y、t 3 個方向上的差分值。

對于連續(xù)的第i 幀和第i +1 幀圖像，點pi的稀疏光流場為

聯(lián)立(12)式、(13)式、(14)式，得

式中:除(xi+1，yi+1)坐標外，其余均為已知。求解(15)式并代入(14)式即可求取點(xi，yi)處的光流矢量:

遍歷直線Li上的點，即可計算出對應的直線稀疏光流場。

2 基于直線稀疏光流場的飛行器姿態(tài)測量

在無人機航拍圖像中，地平線往往是最明顯的標志物，在圖像經過地平線提取算法處理后，地平線更是可以成為視野中唯一的直線，因此特別適合使用直線稀疏光流場對于其運動進行分析。本文將使用地平線附近的稀疏光流場作為例子說明基于直線稀疏光流場的飛行器姿態(tài)測量方法。

在圖2模型的基礎上建立地平線的投影關系模型，如圖3所示，建立大地坐標系OwXwYwZw，機體坐標系OXYZ，圖像坐標系oxy，大地坐標系是以地平線上一點Ow為原點建立的右手系，取豎直向上方向為Yw軸，地平線方向為Xw軸;攝像機與機體頭部固定連接，以攝像機光心O 為原點建立右手系OXYZ，其中Z 軸與光軸重合，取豎直向下方向為Y軸。在焦平面建立圖像坐標系oxy，以圖像中心為原點o，x 軸和y 軸分別與OXYZ 坐標系的X 軸和Y軸平行。

圖3 地平線投影模型Fig.3 Horizon projection model

以坐標系OXYZ 為慣性參考系，在地平線上任取一點P，設其在OXYZ 坐標系下坐標為P(X，Y，Z)，其在焦平面上的對應點為p'(x，y，z)，則點P(X，Y，Z)的運動方程可描述為

式中:ηx、ηy、ηz是點P 的平移運動參數;ωx、ωy、ωz是飛行器滾轉角速度。分別表征俯仰角速度、偏航角速度和滾轉角速度。

將(1)式代入(17)式，得

由于地平線距離很遠，因此在本模型中，可認為(8)式中條件可以得到滿足，這樣由(11)式可得

將(19)式代入(18)式，得

對(20)式整理得

消去Z'得

即可建立光流場與圖像運動之間的關系。

因為點P 在地平線上，而地平線距離觀察點的距離可近似為無窮遠，因此可以近似看作Z→∞，因此消去(22)式中含Z 項，即可得

(23)式是一個關于ωx、ωy、ωz的線性方程組，只要獲取不少于地平線上不少于3 個點的光流矢量值，即可解算出飛行器的3 個角速度。

特別說明的是，本文雖然以地平線作為研究對象說明基于直線稀疏光流場的飛行器姿態(tài)測量方法，但是本方法并不局限于地平線存在的情況，只要無人機攝像機視場中存在典型直線特征(如機場跑道、道路、樓宇等)，本方法同樣適用。

3 仿真與實驗

3.1 直線稀疏光流場計算仿真

使用連續(xù)兩幀測試圖像對于直線稀疏光流場的計算進行仿真，測試圖像如圖4所示，大小為416×240 像素，黑白圖像。兩幀圖像幀間為0.04 s，首先使用Sobel 算子對兩幅圖像進行邊緣提取，如圖5所示。

圖4 連續(xù)兩幀測試圖像Fig.4 Two consecutive images

對第17 幀和第18 幀圖像使用第2 節(jié)中直線稀疏光流場提取算法計算直線稀疏光流場，(16)式中Ix、Iy、It的計算采用如下差分格式:

圖5 測試圖像邊緣提取的結果Fig.5 Edge extraction of test images

式中:δx、δy、δt 為差分步長。在配備Intel i7 CPU 的計算機上，使用MatLab2012 平臺進行數值仿真，計算出的直線光流場如圖6所示，耗時0.223 6 s. 由于沒有該圖像的標準光流場可以對比，使用Horn 方法在同樣平臺上計算第17 和18 幀全局光流場，耗時3.704 4 s，如圖7所示。

圖6 直線稀疏光流場提取算法獲得的稀疏光流場Fig.6 Sparse line-optical flow field calculated by the proposed method

由上述仿真結果可以看出，直線稀疏光流場與Horn 算法一樣，都可以精確描述圖像中邊緣附近的運動情況。由于沒有該序列圖像的標準光流場，使用參考文獻［15］中的光流場誤差測量算法對該算法相對于Horn 算法的進行誤差測量，從而評估該算法。使用本算法計算出的直線稀疏光流場與Horn算法計算出的全局光流場相比，取其邊緣附近對應點的光流值，并按照(25)式計算各個點的角誤差，

圖7 采用Horn 算法計算的全局光流場Fig.7 Global optical flow field calculated by Horn algorithm

式中:wH為Horn 算法計算的光流矢量;wSLOF為本算法計算的光流矢量。

計算本算法得到光流場相對于Horn 算法得到的光流場的相對平均角誤差(RAAE)和相對標準差(RSTD)，見表1.

表1 本算法相對于Horn 算法獲取的光流場誤差Tab.1 RAAE and RSTD errors of globe optical flow field calculated by the proposed method

從上述誤差分析結果來看，使用直線稀疏光流場提取算法獲取的光流場在圖像邊緣附近的精度與Horn 算法獲得的較為接近，平均誤差為4.733 9° ±2.875 2°，因此可以替代后者對圖像中直線特征的邊緣區(qū)域進行光流分析;同時，由于使用本方法的時間開銷要遠遠小于Horn 算法，耗時僅相當于后者的6%，因此本算法更適合諸如無人飛行器視覺導航等對于算法實時性要求較高的場合應用。

3.2 基于直線稀疏光流場的姿態(tài)提取實驗

為了驗證直線稀疏光流場的姿態(tài)提取算法，使用無人機搭載攝像機進行航拍實驗，同時利用機載導航設備記錄姿態(tài)信息。機載導航設備使用ADIS16405 陀螺測量飛行器實時的角速度信息，其動態(tài)響應范圍為±300°/s，采樣率為100 Hz，三軸陀螺儀的角速度測量誤差為±0.05°，加速度計量程±18 g，測量誤差±9 mg. 航拍攝像機采用1/3 in 的CCD 攝像機，視場角是30°，攝像機有效像素1 024×768.

對航拍圖像序列離線計算其直線稀疏光流場，估計飛行器的姿態(tài)信息，并與機載導航設備獲取的姿態(tài)信息進行比較，驗證算法的可行性。

首先對航拍圖像進行地平線提取，提取算法采用Sobel 邊緣提取復合旋轉投影算法提取地平線［2］，在此不再贅述。部分幀地平線提取結果如圖8所示。

圖8 第36 幀和第37 幀地平線提取結果Fig.8 Horizons extracted in 36th frame and 37th frame

依據直線稀疏光流場提取算法，可得兩幀之間的直線光流場如圖8所示，用時0.328 7 s.

圖9為使用本文方法獲取的直線稀疏光流場。為了對比，同時使用Horn 算法計算兩幀圖像間的全局光流場，耗時4.713 s，光流場如圖10所示。由上述仿真結果可以看出，使用直線稀疏光流場提取算法獲取的圖像地平線光流相較于使用Horn 方法計算而得到的全局光流，更不易受到諸如地面水體反射、空中云霧等噪聲的影響，因此可以更準確地估計出圖像中直線部分的光流信息，同時相比Horn 方法具有更快的計算速度。

圖9 使用直線稀疏光流提取得到的稀疏光流場Fig.9 Sparse line-optical flow field calculated by the proposed method

獲取圖像的直線稀疏光流場后，利用(23)式即可完成飛行器俯仰角速度、滾轉角速度和偏航角速度的解算。對100 幀航拍圖像使用上述方法進行計算，并對比機載導航系統(tǒng)角速度陀螺的輸出，如圖11～圖13所示。

圖10 使用Horn 方法計算的圖像全局光流場Fig.10 Global optical flow field calculated by Horn algorithm

圖11 使用稀疏光流場進行的滾轉角速度估計及其估計誤差曲線Fig.11 Roll rate and its errors estimated with sparse optical flow field

結合圖11～圖13可看出，使用稀疏光流場可以用來估計俯仰、偏航、滾轉三通道的飛行器角速度信息，其最大估計誤差除滾轉通道小于±10°/s 外，其余兩通道最大估計誤差均小于±5°/s. 還可發(fā)現(xiàn)，稀疏光流場估計姿態(tài)角的精度與該通道的角速度正相關，在角速度比較小的情況下，使用稀疏光流場還是可以獲得較高的估計精度，而當飛行器角速度增大時，估計誤差隨即增大。

造成這種估計誤差較大的一個因素是圖像信息的采樣率要低于角速度陀螺的采樣率，本飛行實驗采用的陀螺采樣率達到100 Hz，而使用航拍圖像的圖像幀率只有25 幀/s，由于光流估計的數據采樣率不足，從而帶來了一定的估計誤差;另一個因素是本文中只使用了地平線一條直線計算其附近的稀疏光流場，帶來了一定的估計誤差，理論上在進行合理的直線匹配之后，本算法可以獲取視場中多條典型直線特征附近的稀疏光流場，利用多條直線的稀疏光流場進行姿態(tài)角估計將提高本算法的估計精度。

圖12 使用稀疏光流場進行的俯仰角速度估計及其估計誤差曲線Fig.12 Pitch rate and errors estimated with sparse optical flow field

4 總結

本文提出了一種基于直線的稀疏光流場計算方法，并進行了理論分析與仿真，確定了該方法在計算圖像局部光流問題上的有效性;并提出了一種基于直線稀疏光流場的飛行器姿態(tài)角信息提取方法，是現(xiàn)有無人飛行器導航手段的有效替代和補充。通過該方法，可以實時提取飛行器的俯仰角速度、滾轉角速度和偏航角速度信息。對實際飛行實驗獲得航拍數據進行仿真并對比同時采集的導航數據，結果顯示，該算法在滾轉通道角速度誤差在±10°/s，在俯仰和偏航通道角速度估計誤差均小于±5°/s.

需要指出的是，該方法不僅能夠獲得地平線附近的直線稀疏光流場，通過合理的直線匹配，該方法也可以計算出視場中存在的多條典型直線特征附近的稀疏光流場。使用上述多條直線的稀疏光流場來估計飛行器的姿態(tài)信息，將可以提高本算法的估計精度。

圖13 使用稀疏光流場進行的偏航角速度估計及其估計誤差曲線Fig.13 Yaw rate and its errors estimated with sparse optical flow field

References)

［1］ 蔡瑜，葉雄英，朱榮，等. 用于微小飛行器姿態(tài)測量的紅外地平儀研制［J］. 儀表技術與傳感器，2009 (7):33 -35.CAI Yu，YE Xiong-ying，ZHU Rong，et al. Infrared horizon detector making for measuring micro air vehicle attitude［J］. Instrument Technique and Sensor，2009(7):33 -35. (in Chinese)

［2］ Bao G Q，Xiong S S，Zhou Z Y. Vision-based horizon extraction for micro air vehicle flight control［J］. IEEE Transactions on Instrumentation and Measurement，2005，54(3):1067 -1072.

［3］ Pereira G A S，Iscold P，Torres L A B. Airplane attitude estimation using computer vision:simple method and actual experiments［J］.Electronics Letters，2008，44(22):1303 -1305.

［4］ Ettinger S M，Nechyba M C，Ifju P G，et al. Vision guided flight stability and control for micro air vehicles［J］. Advanced Robotics，2003，17(7):617 -640.

［5］ Todorovic S，Nechyba M C，Ifju P G. Sky/ground modeling for autonomous MAV flight［C］∥IEEE International Conference on Robotics and Automation. Taipei:IEEE，2003:1422 -1427.

［6］ Fefilatyev S，Smarodzinava V，Hall L O，et al. Horizon detection using machine learning techniques［C］∥Proceedings of 5th International Conference on Machine Learning and Applications. Orlando，US:IEEE，2006:17 -21.

［7］ Chiu C C，Lo C T. Vision-only automatic flight control for small uavs［J］. IEEE Transactions on Vehicular Technology，2011，60(6):2425 -2437.

［8］ Cornall T D，Egan G K，Price A. Aircraft attitude estimation from horizon video［J］. Electronics Letters，2006，42(13):744 -745.

［9］ Thurrowgood S，Moore R J D，Bland D，et al. UAV attitude control using the visual horizon［C］∥Proceedings of the 2010 Australasian Conference on Robotics and Automation. Brisbane，Australia:IRAA，2010.

［10］ Grzywna J W，Jain A，Plew J，et al. Rapid development of visionbased control for MAVs through a virtual flight testbed［C］∥IEEE International Conference on Robotics and Automation. Barcelona，Spain:IEEE，2005:3696 -3702.

［11］ 程序，郝群，宋勇，等. 基于直線模型的衛(wèi)星飛行器姿態(tài)角計算［J］. 北京理工大學學報，2010，30(7):798 -802.CHENG Xu，HAO Qun，SONG Yong，et al.Method for calculating MAV attitude angles based on straight line model［J］. Transactions of Beijing Institute of Technology，2010，30(7):798 -802. (in Chinese)

［12］ Dusha D，Mejias L，Walker R. Fixed-wing attitude estimation using temporal tracking of the horizon and optical flow［J］. Journal of Field Robotics，2011，28(3):355 -372.

［13］ Shabayek A E R，Demonceaux C，Morel O，et al. Vision based UAV attitude estimation:progress and insights［J］. Journal of Intelligent ＆ Robotic Systems，2012，65(1/2/3/4):295 -308.

［14］ Horn B K，Schunck B G. Determining optical flow［J］.Artificial Intelligence，1981，17:185 -203.

［15］ Lucas B D，Kanade T. An iterative image registration technique with an application to stereo vision［C］∥Proceedings of 7th International Joint Conference on Artificial Intelligence. Vancouver. British Columbia，Canada:San Francico，CA:Morgan Kaufmann Pubishers Inc，1981:674 -679.

［16］ Terzopoulos D. Regularization of inverse visual problems involving discontinuities［J］. IEEE Transactions on Pattern Analysis and Machine Intelligence，1986，8(4):413 -424.

［17］ Hildreth E C. The measurement of visual motion［M］. Cambridge，Massachusetts:MIT Press，1984.

［18］ Barron J L，F(xiàn)leet D J，Beauchemin S S. Performance of optical flow techniques［J］. International Journal of Computer Vision，1994，12(1):43 -77.