基于縱向濾波的星敏感器低頻誤差在線估計

葉立軍，劉付成，尹海寧，徐櫻，寶音賀西

1.清華大學 航天航空學院，北京 100084

2.上海航天控制技術研究所，上海 201109

3.上海市空間智能控制技術重點實驗室，上海 201109

人們對衛星姿態確定精度的要求越來越高，目前，衛星姿態測量敏感器精度最高的星敏感器[1-5]，是眾多航天任務的首選[6]，衛星一般會配置兩臺或多臺星敏感器相互備份，為了進一步提高姿態測量精度，一種可行的方式是將兩臺或多臺星敏感器進行數據融合，以獲得高于每臺星敏感器的姿態確定精度。而多臺星敏感器數據融合的前提是消除星敏感器之間的低頻誤差。

星敏感器低頻誤差按頻率主要可分為3類，第1類表現為星敏感器視場周期[4,7]，主要來源是鏡頭畸變殘差，此類誤差很難從系統層面補償，主要通過單機設計減小，其典型誤差幅值一般為3″～8″[7]；第2類表現為軌道周期[4-12]，主要來源是在軌熱彈性變形，其誤差幅值可達60″[12]，某些衛星幅值甚至高達500″[13]，一般可通過地面后處理補償；第3類表現為年周期，主要來源是光行差[7]及地球繞太陽周年變化引起的熱環境變化等，光行差可在線校正后可忽略，年周期熱環境變化引起的低頻誤差，短期內一般認為不變，當作常值誤差處理。

目前，研究軌道周期低頻誤差估計的文獻很多：龐博等[5]闡述了星敏感器低頻誤差特點及產生原因，介紹了當前幾種主要的星敏感器低頻誤差的抑制與補償方法，指出了當前大多數星敏感器低頻誤差是通過地面處理方式，無法實現誤差實時校正，研究適合在軌實時校正方法是具有實際工程意義的。Winkler[14]、Calhoun[15]、Qiu[16]等提出將星敏感器低頻誤差建模為一階馬爾可夫過程，通過協方差調整對星敏感器的低頻誤差進行濾除，該方法可用于在軌處理不規律的星敏感器低頻誤差，但算法精度與地面基于傅立葉級數處理算法相比較差。此外，當星敏感器出現故障后再接入時，由于算法低頻誤差需要重新估計并收斂，會出現姿態跳變現象，不利于在軌應用。Calhoun[15]、Wang[17]、趙琳[18]、雷琦[19]、熊凱[20]、徐櫻[21]、鐘金鳳[22]等將低頻誤差建模為傅立葉級數，以離線方式，估計出低頻誤差傅立葉系數，實現星敏感器低頻誤差補償。對于多倍頻軌道周期低頻誤差，傅立葉系數較多，算法計算量大。此外，由于年周期的低頻誤差存在，為了保證低頻誤差補償精度，需要定期對傅立葉系數重新估計并更新，增加了地面工作量。Lei和Yang[23-24]參考了Sarkka[25]的部分成果研究了利用無極卡爾曼濾波(UKF)融合Rauch-Tung-Striebel算法，形成URTS固定區間平滑算法，采用事后處理的方式提高了衛星姿態確定精度，但該算法不適用于在線應用。Iwatat[26]、龐博[27]等通過設置地標觀測點，以數據處理的方式實現了星敏感器低頻誤差的估計和補償，提高了星敏感器的姿態確定精度，這種方法需要設置地標觀測點，只適合地面事后處理。

針對現有公開文獻提出算法的不足，本文提出了基于縱向濾波的星敏感器軌道周期慢變系統誤差在線估計與補償算法，可解決星敏感器間軌道周期低頻誤差在線的在線估計問題，利于實現多星敏感器數據在線高精度融合。

1 星敏感器誤差特性

目前針對星敏感器測量誤差分析的學者很多，文獻[15,17-21]采用多階傅立葉級數來模擬軌道周期低頻誤差并作為仿真輸入，但均未考慮星敏感器視場周期和年周期低頻誤差特性因素，不利于客觀評價算法的有效性。本文將采用在軌數據作為仿真輸入，分析星敏感器低頻誤差特性。

經曝光時差校正后星敏感器測量四元數保留了星敏感器間低頻誤差原始信息，根據文獻[28]可計算得兩星敏感器光軸夾角Z(光軸平轉角，體現了兩星敏感器非光軸方向測量誤差)以及兩星敏感器光軸相對旋轉角(光軸旋轉角，體現了兩星敏感器光軸方向測量誤差)：

(1)

(2)

(3)

式中：?表示四元數乘法；q-1為四元數的逆[29]；qs1s2為星敏感器s1到星敏感器s2的旋轉四元數；qis1為星敏感器s1測量四元數，代表慣性系到星敏感器s1本體系旋轉四元數：qis2為星敏感器s2測量四元數，代表慣性系到星敏感器s2本體系旋轉四元數。

根據某衛星在軌遙測下傳的兩星敏感器測量四元數，經式(1)～式(3)計算出兩星敏感器光軸平轉角和光軸旋轉角，與兩星敏感器相應的理論角作差后，得兩星敏感器1和2(s1和s2)沿光軸旋轉角的總誤差如圖1所示。圖中綠色曲線為采用文獻[4]中滑動窗口法計算得低頻誤差(LFE)項。通過圖1(a)和圖1(b)可以看出，兩星敏感器之間低頻誤差主要體現為軌道周期特性，低頻誤差峰峰值約為30″。

進一步地，將總誤差減去低頻誤差，可以得到兩星敏感器測量噪聲，如圖2所示。

值得說明的是，圖1和圖2中的星敏感器誤差(含總誤差、低頻誤差、測量噪聲)為星敏感器2(待修正星敏感器)相對星敏感器s1(基準星敏感器)的誤差。

可以看出，兩臺星敏感器非光軸方向測量噪聲約為7.4″(3σ)，光軸方向測量噪聲約為50.6″(3σ)。根據噪聲不相關原理，且認為兩臺星敏感器測量精度一致，則可得到單臺星敏感器非光軸方向測量噪聲約為5.2″(3σ)，光軸方向測量噪聲約為35.7″(3σ)，與單機指標相符，且與文獻[30-31]中的結論一致：星敏感器光軸方向測量噪聲為非光軸方向測量噪聲的6～10倍。

圖1 兩星敏感器(s1和s2)相對誤差

圖2 兩星敏感器測量噪聲

2 算法設計

2.1 估計和補償策略

星敏感器低頻誤差估計示意如圖3所示，forbit代表軌道圈次。將星敏感器低頻誤差估計，將軌道幅角等間隔均勻離散為N段，將每段中心點作為特征點，這些點的集合稱為星敏感器間低頻誤差估計表，簡稱“估計表”。每個估計表的元素均可看作一個估計子模塊，每個模塊僅估計不同軌道序列相同軌道相位的星敏感器低頻誤差，并分別估計結果進行存儲，最終可將低頻誤差估計問題變為近常值誤差估計問題，采用工程上常用的低通濾波即可實現低頻誤差無偏估計，得到該點星敏感器間低頻誤差值。

圖3 星敏感器低頻誤差估計示意圖

得到估計表后，利用線性插值算法，通過估計表反演出任意時刻當前系統誤差并補償。

2.2 估計和補償算法

星敏感器間低頻誤差在線辨識流程如圖4所示。圖4(a)為星敏感器低頻誤差估計的主流程，圖4(b)為核心算法，是受圖4(a)調用的子流程。

算法流程大致分為以下5個步驟：

步驟1對測量點進行采樣

在相位點R附近對稱取一段數據求平均，提高觀測精度，以更進一步提高估計表的估計精度和收斂速度。

相位點R對應的軌道相位計算方法為

圖4 星敏感器低頻誤差估計流程

(4)

式中：R=0,1,…,N-1。

步驟2低頻誤差估計算法

設基準星敏感器為星敏感器1，待補償星敏感器為星敏感器2，星敏感器2的誤差四元數為

(5)

式中：qs′2s′1為從星敏感器2到星敏感器1的理論轉換四元數，為已知常值；qs′1s1為單位四元數；化簡后可得

(6)

說明：式(6)需要歸一化，并進行標量為正的處理。

步驟3基于縱向濾波估計得星敏感器低頻誤差

第R補償表更新不再是基于傳統的時域濾波，而是基于歷史軌道的同一點進行縱向濾波，其表達式為

(7)

步驟4采用線性插值法計算當前時刻星敏感器2的偏差四元數為

(8)

步驟5對星敏感器2測量四元數進行補償：

(9)

值得說明的是，由于采用線性插值算法，通過線性插值法遍歷軌道上所有點，需要0,1,…,N共N+1個點，物理上第N點與第0點為同一點，因此，當更新第0點時，需要同時更新第N點。

2.3 插值誤差分析

由于軌道上任何位置均處于“估計表”中某兩個相鄰點之間，可采用常規的插值算法反演得到任意軌道點的星敏感器低頻誤差估計值。真實誤差曲線為圓弧形，算法采用線性插值算法計算任意點誤差，計算誤差曲線為折線形，這之間存在的誤差稱之為“建模誤差”。

圖5 星敏感器低頻誤差插值誤差示意圖

(10)

式中：α為點P所處正弦函數的相位，范圍為0～2π。不難看出，α=π/2時，L2取最大值為

(11)

可見，當組成星敏感器低頻誤差的軌道周期基頻階數g越高，補償表個數N越少，則插值誤差越大，當星敏感器低頻誤差特性一定時(g不變)，隨著補償表個數N的增加，插值誤差將以N平方反比關系減小。例如，當g=3時，模型誤差要求小于測量誤差0.01倍，則根據L2<0.01，可得N>66.6。

注：此處用正弦曲線模擬星敏感器在軌軌道周期低頻誤差曲線，為近似做法，僅用于評估插值算法誤差的影響因素和量級，可指導算法參數設計，無法精確計算插值算法誤差，工程設計時需要留好足夠裕量。

2.4 補償點附近采樣區選取

補償點附近采樣區以該補償點為中心左右均勻鋪開，弧段損失建模如圖6所示。

如圖6所示，擴大補償點附近的采樣區，意味著可采更多的采樣點，將采樣點求平均后可提高該補償點的估計精度；但反過來，采樣區越大，相應的弧段損失也越大，導致該補償點估計精度降低，因此需要合理選擇采樣點個數，使補償點精度最高。

(12)

式(12)化簡可得

(13)

經泰勒展開，Δγ較小，忽略高階小量，式(13)可寫為

(14)

可看出，所用弧段越小，則弧段損失越小。

(15)

假設衛星軌道周期為Torbit，則衛星軌道角速度為

(16)

弧段內采樣點個數為

(17)

式中：f為系統采樣頻率。

代入式(15)，有

(18)

采樣點測量誤差為

(19)

對采樣點測量誤差S求二階導，可得

(20)

S″>0，說明S有最小值，令

(21)

得最優弧長表達式為

(22)

表1～表3給出一些工程典型工況最優采樣弧長對比。

表1 測量誤差對最優采樣弧長的影響

通過表1對比分析，測量誤差越大，同等條件下需要更多采樣點以提高精度，所需最優弧長越長。

表2 采樣頻率對最優采樣弧長的影響

通過表2對比分析，采樣頻率越高，同等條件時采樣點越多，可降低對采樣弧長需求，即最優弧長越短。

表3 軌道周期對最優采樣弧長的影響

通過表3對比分析，軌道周期越長，同等條件下采樣點越多，可降低對采樣弧長需求，即最優弧長越短。

以Δγ=1°為例說明弧段損失量級：

(23)

可看出，當采樣弧段為(-1°,1°)，由于弧段損失導致的精度損失比低頻誤差小約6個量級，可忽略不計。

由表1可以看出，經過優化后，采樣值加權平均后，其精度比單機本身測量精度高8～360倍，利于該特征點低頻誤差精確估計。

2.5 補償表個數選取

根據2.3節分析，補償表個數越多，插值誤差越小，并可根據精度需要，確定補償表個數的最低要求。

除了需要滿足插值誤差足夠小要求以外，還需要考慮另外一個問題：一般要求當圈估計的補償表值需要到一個軌道周期后才能使用，這樣通過插值計算的低頻誤差會是平滑的。否則，當算法首次使用或補償表值未完全補償到位時，由于同時使用當圈補償值與上圈補償進行插值計算，會引起計算的低頻誤差出現鋸齒波跳變，不利于工程應用。

上述問題可通過設計采樣弧長與補償點個數N來解決，如圖7所示。

圖7 采樣弧長與補償點個數N的關系

圖7中，兩補償點間任意軌道相位的軌道周期低頻誤差需要通過相鄰兩點補償表插值得到，兩相鄰點補償點間軌道弧長為

(24)

顯然，要避免當軌道圈次估計的補償元素(圖7中的第2點)不與上軌道圈次補償元素(圖7中的第3點)混合使用，需滿足：

Δλ<Δγ

(25)

此外，為了保證各估計表元素對應的數據采集子模塊運行時序，要求兩相鄰數據采集區不能相互重疊，采用單線采集顯然不能滿足上述要求，需采用多線采集方式。

根據圖7不難看出，最少需要3線(即圖7中的Ⅰ，Ⅱ，Ⅲ)采集方式，可以滿足上述約束。即：對于Ⅰ線，負責更新第3M點補償表；對于Ⅱ線，負責更新第3M+1點補償表；對于Ⅲ線，負責更新第3M+2點補償表；其中，M=0,1,…，N/3。

采樣區邊界Δγ取值范圍為

(26)

式中：Δλ′=360°/max(M)。

式(26)化簡可得

(27)

綜合式(24)、式(25)和式(27)，可得

(28)

式中：N為3的整數倍的正整數。

根據表1工程典型工況，取Δγ=1.254 5°具有較好的普適性，建議工程選取N=360。

2.6 算法及其穩定性分析

對于某個補償點來說，星敏感器低頻誤差估計為離散系統，離散周期為軌道周期Torbit。令

(29)

則式(10)可表示為

(30)

式(30)在頻域的表達式為

(31)

可看出，第R補償點估計值逼近理論值的過程本質上是一階慣性環節，一階慣性環節無諧振波，具有絕對收斂性的特性。

此外根據式(29)，補償系數K越大，估計系統時間常數Γ越小。其中，Г為系統修正的時間常數，即估計值達到理論值63.21%所需的時間[32]。

2.7 誤差分析及補償系數選取

根據上述分析，可給出低頻誤差估計殘差δ的表達式，分為兩大部分：δ1為一階系統單位階躍的逼近過程中的偏差；δ2為星敏感器測量輸入的隨機誤差濾波后的誤差。

δ=δ1+δ2

(32)

(33)

(34)

式中：t為算法運行時間；BC為某相位上的低頻誤差真值；σ為星敏感器測量誤差。

可看出，修正的時間常數Г越大，第1部分系統誤差δ1越大，同時第2部分隨機誤差δ2越小。

此外，實際衛星在軌應用時，隨著衛星軌道面進動與地球繞太陽公轉的共同影響，衛星光照條件還有年周期的特性。本算法還需及時響應低頻誤差年周期變化特性，例如：要求Tt時間內，修正后低頻誤差殘差與修正前低頻誤差之比達到ξ，則根據式(33)有

(35)

將式(35)代入式(29)，可得

(36)

實際工程使用時，當低頻誤差中的長期項變化較慢(如年周期)，可適當減小K的取值，以進一步提高濾波精度；反之，當低頻誤差中的長期項變化較快(如天周期)，可適當增加K的取值(如K=1)，以進一步提高對變化的跟蹤響應速度。

3 仿真與分析

3.1 模型誤差仿真驗證

為了驗證算法模型誤差，在仿真時僅保留低頻誤差，不加入基準星敏感器和待修星敏感器測量白噪聲，取補償點個數N=90，補償表更新系數K=1，修正前后星敏感器殘誤大小如圖8所示。

圖8 插值誤差仿真圖

圖8中，第1軌道為修正前星敏感器測量誤差，第2軌道為修正后星敏感器測量誤差，顯然，第2軌道星敏感器測量誤差是插值誤差，可看出，仿真得到的插值誤差與原始誤差之比為

(37)

理論分析得到的算法誤差與原始誤差之比為

(38)

可看出，仿真與理論計算的算法誤差約為0.26%，證明了理論分析的正確性。

3.2 算法仿真驗證

為了驗證算法有效性，采用某在軌衛星兩星敏感器的原始測量四元數，結合兩星敏感器理論安裝四元數，根據本文算法，得到星敏感器2低頻誤差補償過程，仿真時，為了更好體現中間修正過程，補償表更新系數K=0.5，根據式(29)，系統修正的時間常數Г=Torbit。

經過3個軌道周期(即3Г)辨識與修正，根據式(33)，軌道周期低頻誤差殘差為初始值的12.5%，結合第1節星敏感器軌道周期低頻誤差峰值，本案例仿真數據軌道周期低頻誤差峰值將從30″減小至約3.75″，與星敏感器視場周期低頻誤差(3″～8″)相當，基本達到理論上的最佳估計狀態，具體仿真情況如圖9所示。

圖9 星敏感器低頻誤差的辨識與補償過程

圖9(a)為兩星敏感器光軸夾角，圖9(b)為兩星敏感器沿光軸旋轉角，可以看出，經補償，星敏感器低頻誤差中，常值誤差和軌道周期誤差均明顯消除，僅剩下星敏感器視場周期慢變誤差。

修正前后，星敏感器測量誤差峰值統計值對比如表4所示。

表4 修正前后兩星敏感器相對低頻誤差幅值對比

由表4可看出，基于在軌兩星敏感器原始測量四元數，經過本文算法補償后，星敏感器軌道周期低頻誤差基本消除。假設兩星敏低頻誤差特性相同且不相關，則修正后單臺星敏非光軸方向低頻誤差為5″，光軸方向低頻誤差為5.6″，總低頻誤差為7.5″，修正后誤差與星敏感器視場周期低頻誤差幅值相當，修正效果達到預期。

4 結 論

1) 采取縱向的濾波方法，將軌道平均離散為足夠多點，將軌道周期交變的星敏感器低頻誤差估計問題轉換為常值誤差估計問題，并采用一階慣性環節為模型的估計算法，算法收斂，且工程實現簡單可靠。

2) 基于本文算法，當某星敏感器長時間故障，僅會使對應補償表本次數據更新暫停，但該點歷史估計數據(理論上為常值)仍有效，任何時刻星敏感器恢復有效并接入本文算法，不會引起算法重新收斂，待估計值無波動，星敏感器之間基準切換均不會引起姿態跳變，利于工程應用。

3) 工程上可以實現軌道周期和年周期(含遠大于軌道周期)低頻誤差的在線無差估計并實時補償；不適用于星敏感器視場周期(小于軌道周期)的低頻誤差估計。

4) 基于本文算法，星敏感器軌道周期低頻誤差會被消除，殘余的星敏感器低頻誤差即為視場周期低頻誤差，因此本文算法也可用于評估星敏感器視場周期低頻誤差大小。

5) 在本文算法基礎上，配合多星敏感器數據融合算法，可進一步減小系統測量噪聲。