采用壓縮感知的衛星高頻姿態確定方法

汪 璞，胡祥語，安 瑋，盛衛東，林再平，曾瑤源

(1. 國防科技大學 電子科學學院， 湖南 長沙 410073； 2. 中國人民解放軍93147部隊， 四川 綿陽 621000)

高頻姿態抖動是衛星姿態確定中的常見誤差源。這類高頻振動一般會超過常規陀螺和星敏的頻率測量范圍。所以，當高頻姿態抖動存在時，傳統組合定姿系統難以精確測定衛星姿態[1-5]。

目前主流的高頻姿態抖動測定方法是利用高頻姿態傳感器對姿態進行測量。這類研究常見于高精度對地觀測衛星遙感圖像的應用中[6]，例如，日本的先進對地觀測衛星(advanced land observing satellite， ALOS)基于高采樣率的姿態傳感器實現了對高頻姿態抖動的測量。這類方法關鍵在于具備高性能的姿態傳感器，但高性能姿態傳感器通常難以獲得。在國內，衛星上常用姿態的測量頻率通常低于10 Hz[7]。有一些型號的陀螺能夠在較高頻率進行測量，但其通常較重且價格昂貴。所以，很有必要研究利用低頻陀螺實現高頻姿態確定的方法，實現高性價比的姿態測量。

姿態確定方法大致可分為確定性方法和動態濾波估計算法兩類。動態濾波估計算法精度普遍更優，是目前工程應用中的主流算法，后文主要針對該類算法進行介紹。動態濾波估計算法的關鍵在于選擇合適的狀態變量，根據衛星姿態運動學建立起狀態方程，并根據姿態傳感器的測量特性，建立起描述狀態變量與測量值之間關系的測量方程，通過濾波算法降低誤差，得到高精度的定姿結果。目前常用的動態濾波估計算法有：擴展卡爾曼濾波、無跡卡爾曼濾波、非線性預測濾波和粒子濾波等[8]。這些算法各有優勢，能在各自的應用領域取得較好的效果，但其都只能在姿態傳感器的測量頻率范圍內進行姿態確定，一旦出現超過姿態傳感器測量頻率范圍的高頻姿態，相應算法的姿態確定精度就會下降。

為了能在無高頻姿態傳感器的條件下實現高頻姿態抖動的測量，降低抖動對姿態確定的影響，本文提出了一種使用普通姿態傳感器對高頻段姿態確定的方法。該方法主要基于壓縮感知技術進行實現。衛星的高頻姿態抖動由一系列的頻率分量組成，而且它們通常在頻域稀疏[2]。基于這種稀疏特性，利用壓縮感知中的重構算法，可以從低于奈奎斯特采樣頻率的姿態測量數據中恢復出高頻姿態量。

本文所提出的方法包含三個部分：①獲取多組低頻姿態估計值；②融合估計值；③從融合的估計值中恢復出高頻段姿態數據。為了能獲取原始的姿態估計值，所提方法使用了一組星敏以及三組具有不同采樣頻率的陀螺組合測量系統。三個使用星敏和陀螺數據的卡爾曼濾波器被用于估計飛行器的姿態。使用不同的陀螺數據以及相同的星敏數據，每個濾波器均可以產生一組姿態估計值。每組姿態估計值的時間間隔均不相同。根據壓縮感知理論，重構算法要求有時間上非均勻分布的測量值作為輸入。因此，上述三組非均勻分布的測量值，被融合到一起，用以產生一組非均勻分布測量值。基于融合得到的數據，高頻段的姿態數據將被感知并恢復出來。從本文的仿真結果可見，該高頻姿態確定算法有效而且穩定。

1 基于間接卡爾曼姿態濾波器的低頻姿態確定方法

由陀螺和星敏感器組成的姿態確定系統被廣泛應用于具有嚴苛要求的衛星姿態控制中。陀螺和星敏的測量值通常會配合間接卡爾曼濾波器進行使用[9]。本節介紹了采用基于星敏和陀螺的卡爾曼濾波器進行姿態估計的方法。卡爾曼濾波器的整體結構如圖1所示。

圖1 組合定姿卡爾曼濾波器結構Fig.1 Frame of Kalman filter in attitude determination

1.1 四元數與傳感器模型

1.1.1 四元數

姿態用四元數進行表示：

(1)

其中，qe=esin(θ/2)，q4=cos(θ/2),單位向量e表示衛星旋轉軸在某指定姿態參數坐標系(如東南地)中的指向，θ是旋轉角度。

基于四元數的衛星轉動方程可以表示成如下形式：

(2)

(3)

1.1.2 傳感器模型

陀螺輸出yg和星敏輸出ys定義如下：

(4)

1.2 姿態預測與更新

濾波器的目的是從yg和ys中估計姿態四元數q。在本文中，基于間接卡爾曼濾波器進行估計。首先，利用式(5)進行姿態四元數預測：

(5)

利用式(6)更新q：

(6)

δq并不依賴于ωg(角速度)但依賴于陀螺漂移b以及隨機噪聲ηg。所以，盡管旋轉角度較大，可以假設δq較小而且其可以被近似表示為：

(7)

間接卡爾曼濾波器并不直接估計姿態四元數q，而是對δq進行估計，然后利用δq的估計值，結合姿態傳感器的測量值計算得到實際姿態值。

通過對式(5)和式(6)的線性離散化，可以得到：

(8)

1.3 狀態方程

對一個間接卡爾曼濾波器而言，其衛星姿態誤差的狀態向量可以定義為：

(9)

此狀態向量包含六個變量。對卡爾曼濾波器而言，其誤差狀態方程可以定義為：

(10)

其中：

ηb表示陀螺漂移誤差。

為了應用離散時間的卡爾曼濾波器公式，需要將上面的誤差預測公式離散化：

Xk=Mk-1Xk-1+Wk

(11)

其中：Mk-1表示狀態轉移矩陣從時間tk-1向前預測到時間tk，

Mk-1≈I6×6+FΔt

(12)

Δt等于陀螺測量值的采樣間隔；Wk是高斯白噪聲過程，其均值和方程函數定義為

(13)

1.4 測量方程

(14)

(15)

由于姿態角度差異極小(通常小于60角秒，即?3×10-4rad)，誤差矩陣中的高階項均為可忽略小量，它們可以被近似表示為：

(16)

基于式(16)，這些姿態差異角能夠通過式(17)與狀態方程聯系起來：

(17)

通過對式(17)的線性離散化，可以得到：

Zk=HkXk+Vk

(18)

其中，Vk是測量誤差且具有如下特性：

(19)

測量矩陣Hk可定義[9]如下：

(20)

1.5 卡爾曼濾波器

標準卡爾曼濾波器的主要步驟如算法1所示[13]。

算法1 卡爾曼濾波器處理步驟

2 壓縮感知簡介

壓縮感知是21世紀初興起的一種信號復原技術。基于該技術，可以利用遠低于奈奎斯特采樣率的采樣值，精確重構出原始信號[14]。為便于闡述，以如下問題為例來介紹壓縮感知中的基本概念。假設需要從一個不完整測量值中，恢復一個N×1維的信號x，該向量的測量值為：

y=Φx

(21)

其中，Φ是一個M×N的測量矩陣，且M?N；y表示M×1的測量值向量。由于y的維數遠遠低于x的維數，方程(21)是欠定的，難以實現對信號x的重構。然而，如果原始信號x是K稀疏的(即x中僅有K個元素不為0)，則x有可能被精確求解出，其如圖2所示。

圖2 對稀疏信號的線性測量Fig.2 Linear measurement of sparse signal

在圖2中，白色方塊表示0值，彩色方塊表示非0值。在該種情況下，只需確定K個非0元素在x中的位置并針對這K個元素進行求解，而不需要求解x中所有元素的值。這極大地降低了對采樣點數量的要求，即便y中元素數量遠少于x，也可能準確恢復出x。根據壓縮感知理論，當Φ滿足約束等距性準則[15]時，則信號x可以由測量值y通過求解最優l0范數問題精確重構：

(22)

但在實際應用中，大部分自然信號都不是稀疏的，在對這類信號進行重構之前，首先需要通過某種變化Ψ進行稀疏表示，通過對其稀疏系數的重構，間接重構出原始信號。考慮一個非稀疏的自然信號f的重構問題，其不完整測量值y=Φf，f可以被稀疏表示，其稀疏表示系數為x，即：

f=Ψx

(23)

具體如圖3所示。

圖3 信號的稀疏表示Fig.3 Sparse representation of signal

關于信號f的測量方程可以被重寫為：

(24)

圖4 壓縮感知的線性測量Fig.4 Linear measurement in compressive sensing

(25)

在確定壓縮感知的測量方程后，剩下的核心問題就是如何從式(24)中求解出x，該求解過程被稱為信號重構。在信號重構領域，相關研究人員已提出了一系列的信號重構算法，主要有最小l1范數法、匹配追蹤系列算法、迭代閾值法和最小全變分法等。各類方法均有各自特點，需要根據應用環境的不同，選取相應的算法。其中，匹配追蹤類算法由于其計算復雜度相對較低，更適宜于二維信號的恢復重構，本文中即采用該類算法[18-19]。

3 高頻段姿態確定

高頻段姿態確定包含三個主要步驟。首先，利用三個間接卡爾曼濾波器進行姿態估計。由于輸入的陀螺數據有不同的采樣頻率，所以三個濾波器的估計值也具有不同時間間隔。然后，所有濾波器的姿態估計值被融合到一起，用以模擬對初始姿態的非均勻采樣。最后，利用壓縮感知從融合數據中恢復出高頻段姿態數據。受限于計算復雜度，該方法每次僅能恢復出一部分高頻姿態數據，如果要完整恢復出全部的姿態，需要反復使用所提方法。本文所提方法的框架如圖5所示。

圖5 基于壓縮感知的濾波方法Fig.5 Scheme of filter based on compressive sensing

3.1 數據融合

為了模擬對平臺姿態數據的非均勻性采樣，三組數據按照其中各個數據的時間戳點進行了融合。

(26)

其中，t0,t1,…,tk表示姿態數據的時間戳，并且滿足t0

(27)

其中，As表示經過融合后的數據，這個融合后的數據被展示在圖6中。

圖6 融合數據Fig.6 Merged data

As中數據在時間域中非均勻分布。所以，融合數據能被視為一種有效的獲取非均勻分布的姿態序列的方法。

3.2 高頻姿態數據重構

高頻段姿態Ah的重構可以被看作是壓縮感知中一個典型的信號重構問題。它可以被表示為如下形式：

Ash=ΦAh+z

(28)

其中：Φ是一個k×N的測量矩陣；z是誤差項；Ash是Ah的采樣集合。

由于Ah在時域不稀疏，所以其并不能直接通過式(28)進行重構，需要被重寫為如下形式：

Ash=ΦΨxh+z=Θxh+z

(29)

(30)

目前有大量的稀疏重構算法可供選擇。基追蹤算法是其中一種常用的信號恢復算法，它能達到較好的精度而且比較穩定。所以，本文使用基追蹤算法進行重構。由于稀疏基Ψ和測量矩陣Φ對重構比較重要，下面先介紹一下它們的定義。

1)稀疏基。基于前期的研究[12]，航天器的姿態角可以用式(31)進行描述：

(31)

其中，A(t)代表了一個軸上的姿態角，Ai表示正弦函數的幅度，φi是正弦函數的相位角，fi表示該函數的頻率。

由式(31)可知，A(t)經過傅里葉變換后在頻域稀疏。所以，傅里葉基被選為稀疏基，其定義如下：

(32)

2)測量矩陣。令Ts表示As中元素所對應的時間，該測量矩陣直接與時間序列Ts相關聯。假設Ts=(t0,t1,…,tk)，且Δt是三個陀螺采樣時間間隔(Δt1,Δt2,Δt3)的最大公約數，通過將Ts除以Δt能夠獲得一個位置向量Sk。其定義如下所示：

(33)

基于Sk，可以獲得K個行向量b1,b2,…,bk：

(34)

在向量bi中，只有第pi個元素為1，其他所有元素均為0。

然后，該測量矩陣可以被表示為：

(35)

4 仿真實驗測試

為了驗證本文所提姿態確定方法的有效性，利用仿真數據對所提算法進行了測試。將其與間接卡爾曼濾波器進行比較測試。為了進一步闡述所提方法的有效性，在不同噪聲條件下對該方法進行了測試，并對不同范圍下的頻率進行了測試。

4.1 仿真條件

假設有三組陀螺和一組星敏被安裝于一個三軸穩定的衛星之上。每組陀螺采樣間隔均不相同，分別為[55 ms(18.2 Hz)85 ms(11.8 Hz)95 ms(10.5 Hz)]，三組陀螺的測量時間間隔為5 ms。為簡便起見，假設一組之中的三個陀螺平行于飛行器的三軸且其性能均相同。星敏的采樣頻率設為1 s(1 Hz)。

仿真時間持續了100 s，其仿真的步長為5 ms。 星敏和陀螺的測量數據通過數字仿真獲得。陀螺所采用的誤差為角度的隨機游走誤差以及角速率的隨機游走誤差，其誤差率分別為5 μrad/s 和0.5 μrad/s，其中μrad表示微弧度，即百萬分之一弧度。三組陀螺有相同的誤差參數。星敏的測量誤差模型僅包括正態分布的隨機誤差(15 μrad)。

衛星姿態可以通過式(31)進行描述。假設姿態包括14個頻率分量，仿真參數被列在表1中。衛星的高頻振動分量通常幅度較小，所以在仿真中將其設為一個小量。

表1 高頻段姿態參數Tab.1 Parameters of high frequency attitude

在后續的魯棒性測試和帶寬測試中，部分條件進行了改動。為了測試所提出方法的魯棒性，調整了陀螺和星敏誤差大小，此外還測試了該方法在姿態突變情況下的性能。在帶寬測試中，將頻率14在20～200 Hz范圍內進行了調整。

4.2 與間接卡爾曼濾波器的比較

在本節中，將卡爾曼濾波器與本文所提方法進行了比較。卡爾曼濾波器使用1 Hz的星敏數據和18.2 Hz的陀螺數據。圖7展示了真實的姿態角、通過間接卡爾曼濾波器獲得的姿態估計值以及本文所提方法獲得的姿態估計值。為了更精確地展示各種方法的效果，僅將部分圖像放大進行了展示。由圖7可知，使用間接卡爾曼濾波器得到的結果中缺乏高頻量，而本文所提姿態估計方法得到的結果幾乎與真實值一致。

圖7 間接卡爾曼濾波器與本文所提方法的比較Fig.7 Comparison between indirect Kalman filter and proposed method

間接卡爾曼濾波器和本文所提方法的三軸姿態估計誤差被展示在圖8中。為了更清楚地展示兩種方法性能，圖中只展示了部分結果(59～60 s)。從圖8(a)中可以看出，間接卡爾曼濾波器的估計誤差主要由高頻姿態抖動組成。這類傳統的姿態確定方法不能提供高頻段的姿態估計信息。間接卡爾曼濾波器估計誤差的均值和均方根分別為[0.59×10-5-1.12×10-50.24×10-5]rad和[1.66×10-31.63×10-31.67×10-3]rad。圖8(b)顯示了基于壓縮感知的姿態確定方法的估計誤差主要是隨機噪聲，其中的高頻誤差量已經被消掉。誤差的均值和均方根值分別為[0.76×10-5-1.13×10-50.34×10-5]rad和[2.57×10-52.54×10-56.10×10-5]rad。兩種方法的誤差均值均為接近0的小量，本文所提方法的誤差均方根有明顯降低。

通過圖7和圖8可以發現，本文所提方法可提高對高頻姿態量的估計精度。

(a) 間接卡爾曼濾波器(a) Indirect Kalman filter

(b) 本文所提方法(b) Proposed method圖8 不同方法的估計誤差Fig.8 Estimation error of different methods

4.3 魯棒性測試

4.3.1 不同星敏陀螺測量參數

為驗證算法的普適性，在4.1節仿真條件的基礎上，通過對星敏陀螺測量頻率的調整(其他條件不變)重新配置了三組仿真場景，以測試在不同星敏陀螺條件下算法的性能。仿真條件如表2所示。

表2 星敏陀螺仿真參數Tab.2 Parameters of star sensor & gyros simulation

根據表2中設定的星敏陀螺條件，模擬測試了高頻定姿算法在這三種條件下的姿態確定精度，得到的結果如表3所示。

表3 不同星敏陀螺測量條件下姿態確定精度Tab.3 Attitude determination precise of different scenes

由以上仿真結果可知，在不同的星敏陀螺測量條件下，高頻定姿算法的姿態確定殘差依然保持在同一量級，達到較高的姿態確定精度，體現了其良好的普適性。

4.3.2 不同噪聲干擾測試

為了展示本文所提方法的可靠性，在不同噪聲條件下對本方法進行了測試。在不同的噪聲測試中，所有其他的仿真條件都和前面的仿真實驗相同，僅對噪聲的幅度進行了調整。相應的噪聲水平和噪聲幅度如表4所示。

表4 噪聲水平定義Tab.4 Definition of noise level

本文所提方法是通過姿態估計誤差平均值以及均方根誤差對方法進行評價的。不同噪聲水平下的平均誤差和均方根誤差如圖9所示。

(a) 俯仰角(a) Pitch angle

(b) 滾動角(b) Roll angle

(c) 偏航角(c) Yaw angle圖9 噪聲測試結果Fig.9 Result of noise test

在圖9中，姿態估計誤差的均值在0附近而均方根誤差隨噪聲水平的增長而線性增加。輸入誤差的微小改變僅僅會導致輸出的微小改變，所以可以認為本文所提的方法是穩定的。

4.3.3 姿態突變測試

在實際應用中，有時會出現姿態突變情況,姿態會在一個很短的時間內出現劇烈變化。由于這種姿態突變在頻域中并不稀疏，其將會影響姿態重構的精度。所以，很有必要測試本文所提算法在出現姿態突變時的效果。

在該測試中，除添加了一次姿態突變以外，所有的仿真條件都和第一次仿真測試相同。姿態突變發生在50 s左右，其幅度被設為0.02 rad。

姿態突變條件下，間接卡爾曼濾波器和本文所提方法的三軸姿態估計值如圖10所示，估計誤差如圖11所示。

圖10 姿態突變條件下估計值Fig.10 Estimation result in attitude jump

圖11展示了在突變存在的情況下，本文所提方法的估計誤差出現了劇烈變化，但是這些誤差很快消退，而且這些估計誤差小于間接卡爾曼濾波器。盡管突變會影響本文所提算法的效能，但是該方法能夠很快恢復。

(a) 間接卡爾曼濾波器(a) Indirect Kalman filter

(b) 本文所提方法(b) Proposed method圖11 姿態突變條件下估計誤差Fig.11 Estimation error in attitude jump

4.4 帶寬測試

為了確定文中所提姿態確定算法的帶寬，該方法被用以測試不同的高頻分量。高頻姿態抖動主要是由衛星通信天線、反射鏡轉動以及動量輪旋轉等引起。這些因素引起的誤差振動通常小于200 Hz，而在本方法中使用的最高陀螺頻率接近20 Hz，所以將測試的高頻范圍限定在20～200 Hz。其他仿真條件與第一次仿真相同，只有頻率分量14在20～200 Hz之間進行了調整，每次的調整步長為1 Hz。 不同高頻分量所對應的姿態估計誤差的均值以及均方誤差被展示在圖12中。

隨著頻率增加，姿態估計誤差的均值和均方誤差值都無明顯變化。實驗中均方誤差通常在10 μrad以下，在最壞的情況中，均方誤差也是在20 μrad附近。圖12展示了所提出方法在20～200 Hz范圍內能夠對超出姿態傳感器測量范圍的姿態進行確定。姿態采樣頻率小于20 Hz時，由于在姿態傳感器的測量范圍內，必然可以被測量到，所以未做相應測試。

(a) 俯仰角(a) Pitch angle

(b) 滾動角(b) Roll angle

(c) 偏航角(c) Yaw angle圖12 在不同頻率的測試誤差Fig.12 Estimation error at different frequency

5 結論

本文提出了一種新的姿態確定方法。該方法僅用普通采樣頻率的姿態傳感器就可實現航天器高頻姿態估計，且能達到較高精度。該方法利用了航天器高頻姿態量在頻域的稀疏性，通過壓縮感知技術，可以從低頻的姿態測量數據中恢復出衛星高頻姿態數據。通過仿真結果可以發現，該方法能夠測量出高頻姿態，并且對誤差具有較好的魯棒性。此外，通過增加姿態傳感器，改進重構算法，可以進一步改善高頻姿態測量的精度。