劉 鵬， 倪鄭鴻遠， 趙 巍*， 王志強

0 引 言

從第一顆人造衛(wèi)星的升空到成功踏上月球的土地，人類對外太空的探索取得了可觀的成績?nèi)缃?，人類已將探索的目光從大型天體轉(zhuǎn)向小行星、彗星等小天體.近幾十年來，美國、日本、歐空局都相繼開展了小天體探測任務[1-4].

由于小天體體積小、距離地球遙遠，無法在地球上得到其物理參數(shù)、形狀以及表面細節(jié)的可靠估計，并且由于小天體形狀不規(guī)則，觀測難度大，探測器很難獲取其大量的高分辨率的圖像，所以利用圖像中小天體的輪廓信息來恢復三維模型的方法具有重要的研究意義.對小天體旋轉(zhuǎn)參數(shù)的估計包括旋轉(zhuǎn)軸方向估計和自轉(zhuǎn)角速度估計，典型方法包括光變曲線法和視覺運動參數(shù)估計法.用光變曲線估計小天體的旋轉(zhuǎn)參數(shù)假設小天體為橢球體，通過光變曲線的變化估計小天體的旋轉(zhuǎn)軸指向和自轉(zhuǎn)速度.物體運動參數(shù)估計也是計算機視覺中的一項重要研究內(nèi)容，在探測器導航著陸、目標跟蹤等方面都會涉及到這一問題.

本文研究一種基于小天體探測接近過程中使用圖像上特征點跟蹤和擴展卡爾曼濾波器的小天體旋轉(zhuǎn)參數(shù)估計方法.該方法首先建立小天體旋轉(zhuǎn)關系模型，表示小天體在相機坐標系中的姿態(tài)變化；然后定義小天體旋轉(zhuǎn)的狀態(tài)方程，通過對觀測圖像序列上的特征點跟蹤，利用擴展卡爾曼濾波器方法得到小天體旋轉(zhuǎn)軸指向和自轉(zhuǎn)角速度估計.在實驗中，分析了相機坐標系小天體坐標系之間的四元數(shù)初值，圖像上特征點跟蹤精度，相機的觀測指向等因素對小天體旋轉(zhuǎn)參數(shù)估計精度的影響.實驗結(jié)果表明，本文提出的基于圖像特征點跟蹤和擴展卡爾曼濾波器的小天體旋轉(zhuǎn)參數(shù)估計方法能夠得到具有較高精度的估計結(jié)果.

1 相關工作

1.1 基于光變曲線的小天體旋轉(zhuǎn)參數(shù)估計

當小天體在視野中成像為小面積擴點目標狀態(tài)時，探測器相機無法觀察到小天體的輪廓形狀和具體細節(jié)，此時光變曲線的信息可以為旋轉(zhuǎn)軸的估計提供參考.假設小天體的形狀用橢球體表示，記橢球體的三個半軸長度分別為a,b,c且a≥b≥c，假設c為旋轉(zhuǎn)軸可以使橢球體的動能最小，即處于穩(wěn)定狀態(tài).在入射光為平行光且能量相同的情況下，小天體的反射亮度會與其視截面的面積大小成正比，而其視截面的面積大小會隨著小天體的旋轉(zhuǎn)呈周期性變化.同時，當相機的視線方向與小天體旋轉(zhuǎn)軸的夾角發(fā)生變化時，視截面的面積也會發(fā)生相應的變化，進而導致光變曲線的結(jié)果不同.通過光變曲線的振幅反演可以得到小天體的旋轉(zhuǎn)軸指向和旋轉(zhuǎn)周期[5].

1.2 基于視覺的旋轉(zhuǎn)參數(shù)估計方法

運動參數(shù)估計是計算機視覺領域的一項重要的研究內(nèi)容，在相機標定、探測器導航著陸、目標跟蹤等方面都會涉及到這一問題.按照估計的原理可以分為直接法、光流法和基于特征的方法三類[6].

直接法是最簡單的運動參數(shù)估計方法，它既不需要提前進行特征提取和匹配也不需要計算光流，直接利用圖像中每個像素點的亮度來估計物體的運動.文獻[7]提出了一個直接利用亮度梯度來估計物體運動信息的方法，該方法對兩幀圖像選擇一個任意的固定點，這個點包含一個面片(patch)和對應的速度.通過解算兩個線性方程組，可以得到像素平移參數(shù)，將這些平移參數(shù)應用于第二幀圖像中的每一個像素點就得到了固定的第二幀圖像.最后求解3×3矩陣的特征值和特征向量，利用最小特征值對應的特征向量，第一幀圖像以及固定的第二幀圖像可以求出平移方向、旋轉(zhuǎn)速度和相對深度.Heel提出了基于卡爾曼濾波器的直接估計場景結(jié)構(gòu)和相機相對于場景的位姿的方法.隨著時間的推移，卡爾曼濾波器可以逐漸地提升預測的效果.在更新階段，卡爾曼濾波器會結(jié)合當前的深度估計和測量值對估計做出校正.

Gibson提出了光流的概念：物體在空間中運動，由相機成像后，圖像上像素運動的瞬時速度稱為光流.對場景中的運動物體進行連續(xù)拍攝得到圖像序列，相鄰圖像中物體運動具有一定的相關性，光流法利用這種相關性，通過分析像素隨時間產(chǎn)生的變化，可以計算出物體的運動信息.根據(jù)引入的約束條件的不同，分析方法也不同，例如基于梯度的Lucas-Kanade算法[8]、區(qū)域匹配法等.由于光流法假設物體亮度不變，且計算量較大無法進行實時預測，所以較難應用于小天體的運動參數(shù)估計中.

基于特征的方法是目前位姿參數(shù)估計領域的研究熱門，也是成功應用于小天體運動參數(shù)估計的方法，首先要從圖像中提取的特征，然后進行特征匹配，根據(jù)跟蹤的特征點信息來估計物體的運動參數(shù).根據(jù)使用的特征的不同，方法的適用場景也不相同.常用的特征提取算法可以分為尺度不變特征點算法和角點特征點算法.SIFT特征[9]是尺度不變特征檢測算法的代表，這一類特征的優(yōu)點是對圖像的尺度、旋轉(zhuǎn)以及亮度變化不敏感，但是執(zhí)行效率較低；而以Harris角點算法計算簡單，但是由于不具有尺度不變性.基于特征的方法可以分為兩類：批處理算法和遞歸算法.批處理算法不滿足實時性，而遞歸算法只需要上一次的狀態(tài)估計和當前的測量值，就能得到當前的狀態(tài)估計，占用的存儲資源少，執(zhí)行效率高.BROIDA提出了基于循環(huán)擴展卡爾曼濾波器(IEKF)的遞歸算法[10].IU等人提出的方法將旋轉(zhuǎn)平移表示為任意階多項式的形式，將問題轉(zhuǎn)化為多項式系數(shù)的估計問題[11].

2 小天體旋轉(zhuǎn)參數(shù)估計

2.1 姿態(tài)表示和觀測方程

小天體旋轉(zhuǎn)參數(shù)估計主要涉及三個坐標系：世界坐標系Ow-XYZ、小行星坐標系OS-xyz、相機坐標系Oc-xyz，如圖1所示.

圖1 坐標系的定義Fig.1 Definition of coordinate systems

用四元數(shù)表示相機坐標系與小天體坐標系之間的姿態(tài)關系.對于一個四元數(shù)q，其數(shù)學代數(shù)形式為：q=q1+iq2+jq3+kq4，去滿足正交化約束

(1)

設t時刻小天體坐標系原點Os在相機坐標系Oc-xyz中的坐標定義為

S0(t)=[x0(t)y0(t)z0(t)]T

(2)

在小天體坐標系OS-xyz中的小天體表面一點Si0,i=1,2,…,M,M為特征點的數(shù)量.在相機坐標系Oc-xyz中的坐標可以表示為：

Si(t)=S0(t)+R(t)Si0(t)=[xiyizi]T

(3)

t時刻Si(t)對應的圖像坐標為：

(4)

2.2 狀態(tài)變量和狀態(tài)方程

M個特征點在小天體坐標系中的坐標表示為：

[x1y1z1]T,[x2y2z2]T,…,[xMyMzM]T

定義狀態(tài)變量s(t)：

(5)

狀態(tài)變量的描述將表1.

表1 狀態(tài)變量的描述Tab.1 Description of state variables

狀態(tài)變量中包含了小天體和相機之間的旋轉(zhuǎn)關系，在相機坐標系下的小天體坐標系原點位置、角速度和線速度，以及選取的特征點在小天體坐標系下的位置.利用小天體坐標系原點在相機坐標系下的z軸位置z0(t)作為狀態(tài)變量中所有坐標的歸一化因子.定義了狀態(tài)變量s(t)后，小天體位移-旋轉(zhuǎn)狀態(tài)轉(zhuǎn)移方程可以用式(6)表示：

st+1=f(st)+wt

(6)

其中wt為均值為0，協(xié)方差矩陣為Qt的高斯隨機噪聲，f(·)為狀態(tài)轉(zhuǎn)移函數(shù).狀態(tài)變量的微分近似為：

(7)

其中：

=s3(t)-s1(t)s5(t)

=s4(t)-s2(t)s5(t)

=-s3(t)s5(t)

=-s4(t)s5(t)

=0.5[s12(t)s7(t)-s11(t)s8(t)+s10(t)s9(t)]

=0.5[-s12(t)s6(t)-s10(t)s8(t)+s11(t)s9(t)]

=0.5[s11(t)s6(t)-s10(t)s7(t)+s12(t)s9(t)]

=0.5[-s10(t)s6(t)-s11(t)s8(t)-s12(t)s9(t)]

…

=-s3M+1(t)s5(t)

=-s3M+2(t)s5(t)

=-s3M+3(t)s5(t)

定義測量值為M個特征點對應的圖像坐標為：

p(t)=[X1(t)Y1(t) …XM(t)YM(t)]T

(8)

上式中的元素表示為：

(9)

將

x0(t)/z0(t)=s1

y0(t)/z0(t)=s2

[xiyizi]T/z0(t)=[s3i+10s3i+11s3i+12]T

[r11r12r13]=

[r31r32r33]=

代入式(8)，得觀測方程：

p(t)=h(st)+vt

(10)

其中vt為均值為0，協(xié)方差矩陣為Rt的高斯噪聲矢量.

2.3 擴展卡爾曼濾波器估計狀態(tài)

將狀態(tài)轉(zhuǎn)移方程(6)和觀測方程(10)寫成標準的擴展卡爾曼濾波器形式：

(11)

擴展卡爾曼濾波器的計算過程如圖2所示：

圖2 擴展卡爾曼濾波器計算過程Fig.2 Calculation process of extended Kalman filter

3 實驗結(jié)果

采用真實的67彗星三維模型分別在探測器距目標100 km和20 km兩種的情況下進行實驗，探討了四元數(shù)初值，圖像特征點精度和相機觀測方向等對于小天體旋轉(zhuǎn)參數(shù)估計精度的影響.

3.1 實驗條件

目標：使用67p三維模型，位于世界坐標系原點，繞z軸每幀轉(zhuǎn)1°.

相機：光軸指向世界坐標系原點，渲染圖像的分辨率為1024×1024，像元尺寸為8.6 μm，焦距為630.243 mm，共180幀.

特征點選?。涸谌S模型上選擇4個特征點，如圖3所示，這些特征點在小天體坐標系下的坐標：

p1(-2354 cm, 14124 cm, -4268 cm)，

p2(7581 cm, 7663 cm, -15264 cm)，

p3(-7745 cm, 2840 cm, -12237 cm)，

p4(5660 cm, 1347 cm, -19400 cm).

圖3 在67p彗星上選取的特征點Fig.3 Selected feature points on comet 67p

觀測值：已知特征點的三維坐標，利用三維重建后得到的相機外參數(shù)和相機本身的內(nèi)參數(shù)，可以計算出特征點在每一幅圖像上的坐標.

旋轉(zhuǎn)軸指向誤差：擴展卡爾曼濾波器預測得到的旋轉(zhuǎn)軸方向與模型真實旋轉(zhuǎn)軸方向的夾角(相機坐標系下).

角速度誤差：模型設置的角速度與濾波器預測角速度的差值(在相機坐標系下3個坐標軸方向).每一次實驗均進行100次，取結(jié)果均值.

3.2 四元數(shù)初值對旋轉(zhuǎn)參數(shù)估計精度的影響

在設置67彗星坐標系到相機坐標系的三軸旋轉(zhuǎn)關系為(-45°,0°,0°)，對四元數(shù)逐步加入5°、10°、20°的偏差，對67p旋轉(zhuǎn)軸估計誤差如圖4所示.實驗結(jié)果表明，在距離目標100 km時，用本文方法估計67p旋轉(zhuǎn)軸指向?qū)τ谒脑獢?shù)初值誤差不敏感.

圖4 相機與67p彗星相距100 km時，四元數(shù)初值對旋轉(zhuǎn)參數(shù)估計精度的影響Fig.4 When the distance between the camera and comet 67p is 100 km, the effect of initial quaternion value on the estimation accuracy of rotation parameters is as follows

3.3 特征點跟蹤精度對旋轉(zhuǎn)參數(shù)估計精度的影響

對觀測的特征點圖像坐標加入均值為0，標準差分別為0.04 mm，0.08 mm，0.12 mm的高斯噪聲(像素尺寸為8.6 μm)，67彗星旋轉(zhuǎn)軸指向估計誤差如圖5所示.實驗結(jié)果表明，在距離目標100 km時，特征點跟蹤誤差的標準差不應超過0.04 mm.

圖5 相機與67p彗星相距100 km時，圖像特征點跟蹤精度對旋轉(zhuǎn)參數(shù)估計精度的影響Fig.5 When the distance between the camera and comet 67p is 100 km, the influence of the tracking accuracy of image feature points on the estimation accuracy of rotation parameters

3.4 相機光軸方向?qū)πD(zhuǎn)參數(shù)估計精度的影響

假設在100 km位置已經(jīng)得到關于小天體旋轉(zhuǎn)軸指向，分別設置探測器在小天體坐標系以下3個位置觀測：(0,0,-20 km)、(0,14 km,-14 km)、(0,20 km,0)(即與旋轉(zhuǎn)軸夾角為0°、45°、90°).實驗結(jié)果如圖6所示.實驗結(jié)果表明，在探測器距離小行星20 km時，在兩極(觀測方向與小天體旋轉(zhuǎn)軸夾角為0°)，觀測可以獲得最好的估計結(jié)果，在赤道處(觀測方向與小天體旋轉(zhuǎn)軸夾角為90°)觀測雖然也可以獲得不錯的結(jié)果，但是由于特征點的消失和加入導致預測過程出現(xiàn)波動，需要較長的迭代次數(shù)才能收斂.

4 結(jié) 論

本文研究的基于圖像特征點跟蹤的擴展卡爾曼濾波器小天體旋轉(zhuǎn)參數(shù)估計方法，定義了用于小天體旋轉(zhuǎn)參數(shù)估計的狀態(tài)變量，推到了狀態(tài)轉(zhuǎn)移方程和測量方程，詳細給出了計算過程.用67P彗星三維模型進行了實驗.實驗結(jié)果表明本文方法能夠在相機與目標小天體相聚100 km距離時，估計出小天體的旋轉(zhuǎn)參數(shù)；當相機與目標小天體相聚20 km距離時，將相機主光軸方向與小天體旋轉(zhuǎn)軸方向的夾角控制在45°以內(nèi)，本文方法能夠魯邦地估計小天體旋轉(zhuǎn)軸指向和自轉(zhuǎn)速度.