張軍華，陳建林，孫沖，袁建平

1. 西北工業(yè)大學 航天學院航天飛行動力學技術(shù)重點實驗室，西安 710072

2. 西北工業(yè)大學 民航學院，西安 710072

1 引言

隨著商業(yè)航天的發(fā)展，低軌小衛(wèi)星平臺受到越來越廣泛地關注。小衛(wèi)星相對于大衛(wèi)星來說，由于體積小、質(zhì)量輕、成本低，具有更靈活的應用潛力，因此近年來小衛(wèi)星廣泛應用于近地航天任務，已經(jīng)成為航天工程的研究熱點之一[1-2]。

通常，小衛(wèi)星受限于質(zhì)量約束，無法像大衛(wèi)星平臺一樣攜帶大量的燃料，因此高比沖微小推力推進器小衛(wèi)星的常用推進方式，主要包括電推進、太陽帆推進、電帆推進、冷氣推進等[3-10]。隨著航天技術(shù)的發(fā)展，當前的微小衛(wèi)星任務對軌道預測和估計的精度要求不斷提高，這不僅對衛(wèi)星軌道動力學建模精度提出了較高的要求，還對軌道預測和估計方法提出了精度要求。用于軌道的預測和估計的數(shù)值積分方法，如Dormand-Prince方法、Gauss-Jackson方法、Adams-Bashforth-Moulton方法等往往在誤差控制、穩(wěn)定性等方面存在問題[11-15]。而采用簡化的動力學模型往往導致軌道預測的精度無法滿足實際需求[16-19]。

進一步，由于微小推力的量級小，難以嚴格在地面完成高精度的推力標定，需要衛(wèi)星在軌對其推力進行標定。常見的推力標定包括參數(shù)標定、飛輪標定、矢量標定等[20-23]，這些標定方法采用衛(wèi)星載荷對標定所需參數(shù)進行估計。但小衛(wèi)星搭載載荷有限，上述標定方法無法適用于小衛(wèi)星的推力在軌標定。另一方面，微小推進器存在推力損失等不確定性因素，如何保證推力在軌標定的精度，是個亟待解決的問題。

本文針對小衛(wèi)星的高精度軌道預測和微小推力在軌標定問題，對主要的攝動力進行了高精度建模。利用變步長龍格庫塔積分RK78對衛(wèi)星軌道進行了外推，通過與STK仿真結(jié)果對比驗證了模型的精確性；使用無跡卡爾曼濾波器(UKF)同時估計小衛(wèi)星軌道并標定衛(wèi)星推力，證明了在軌估計算法的可行性。本文對軌道遞推和推力標定算法的研究，可以應用于小衛(wèi)星的實時在軌運算中，實時掌握小衛(wèi)星自身的軌道狀態(tài)和推力大小，同時算法的速度較快，有利于小衛(wèi)星在線運算任務的實現(xiàn)。

2 微小推力下小衛(wèi)星動力學模型

2.1 問題描述

小衛(wèi)星在近地軌道上運動時，除了受到地球中心引力場的引力外，還受到地球非球形引力、太陽光壓和第三體引力等攝動影響，小衛(wèi)星采用微小推力進行軌道控制。對微小推力下小衛(wèi)星的在軌運行情況進行預測，不僅需要精確的軌道運動模型，而且對軌道預測和估計的精度也有較高的要求。此外，由于微小推力量級小，且小衛(wèi)星任務過程中存在不確定性，需要將微小推力進行在軌標定。

2.2 動力學模型

首先建立三個坐標系[24]，如圖1所示。

圖1 坐標系示意圖Fig.1 Illustration of coordinate frames

地球慣性坐標系EXYZ(JDate)，坐標系原點在地心E，X軸指向當前時刻春分點，Z軸垂直于地球赤道平面向上，Y軸滿足右手法則。

地球慣性坐標系J2 000，坐標系原點在地心E，X'軸指向J2 000時刻的平春分點，Z軸垂直于地球赤道平面向上，Y'軸滿足右手法則。

衛(wèi)星軌道坐標系Sxyz，坐標系原點在衛(wèi)星質(zhì)心S，x軸沿衛(wèi)星矢徑方向，從地心指向衛(wèi)星，y軸垂直衛(wèi)星矢徑并指向速度方向，z軸垂直于軌道平面，遵循右手法則。

除此3個坐標系外，建立地球固聯(lián)坐標系(ECEF)，其定義為：坐標系原點在地心，X''軸從地心指向赤道與國際參考子午線的交點，Z''軸沿著地球自轉(zhuǎn)軸從地心指向北極點，根據(jù)右手準則Y''軸從地心指向赤道與90°東經(jīng)子午線的交點。

另外，假設衛(wèi)星本身采用三軸對地穩(wěn)定，因此衛(wèi)星的本體坐標系與軌道坐標系是相重合的。

衛(wèi)星受到微小推力加速度aF的作用，且該微小作用力在本體坐標系中的方向不變，如圖2所示。因此，aF在軌道坐標系下可以由θ和φ兩個方位角表示，aF=[aFcosφcosθ,aFcosφsinθ,aFsinφ]T。

圖2 微小推力在軌道坐標系下的示意Fig.2 Illustration of micro-thrust in frame Sxyz

衛(wèi)星受到微小推力、地球扁率(J2、J22等)、太陽光壓以及第三體(日月)引力等攝動，則衛(wèi)星在地心慣性坐標系下的動力學方程可以描述為：

(1)

(2)

式中：ax,nm、ay,nm和az,nm分別表示地球形狀攝動加速度在地球固聯(lián)坐標系中三個方向上的分量：

(3)

式中：R⊕表示地球半徑；地勢系數(shù)Cnm,Snm由EGM96S引力模型給出；Vnm和Wnm分別為：

(4)

(5)

在地球慣性系下的太陽光壓和太陽、月球等第三體引力引起的攝動力為：

(6)

(7)

式中：r⊕,rSun,rMoon,r分別表示地球慣性系下衛(wèi)星到太陽的位置矢量；太陽的位置矢量；月球的位置矢量和衛(wèi)星的位置矢量；G表示萬有引力常數(shù)；AU表示地球到太陽的距離；P表示單位AU距離下的太陽光壓力；ν表示地球的陰影函數(shù)；M表示太陽或者月球的質(zhì)量；Cr表示輻射壓力系數(shù)；m表示衛(wèi)星的質(zhì)量；A表示衛(wèi)星的受照橫截面積。

因此，其他攝動總加速度可以表示為：

(8)

式中：aFx、aFy、aFz、aJx、aJy、aJz、adx、ady、adz分別表示微小推力加速度，地球形狀攝動加速度，其他攝動加速度在慣性坐標系下的三個坐標軸的分量。

3 衛(wèi)星軌道遞推算法

衛(wèi)星常見的星上定軌方式是基于GNSS信號，但是由于GNSS輸出的軌道信息一般在位置和速度上存在誤差，直接應用GNSS信號進行軌道確定，會造成衛(wèi)星軌道定位不精確。為了解決這一問題，提出一種變步長龍格庫塔算法RK78進行軌道外推，具體的算法模型如下：

Runge-Kutta方法是一種間接采用泰勒級數(shù)展開而求解常微分方程初值問題的數(shù)值方法。該方法用某些點上函數(shù)值的不同組合來近似初值問題的解析解，從而避免函數(shù)高階偏導數(shù)的計算，并獲得高階的方法[26]。將Runge-Kutta公式與初值問題解的泰勒展開式進行比較，使之有某些項相同，便可以確定組合系數(shù)，從而建立所需的方法。

由m個函數(shù)值線性組合構(gòu)成的顯式Runge-Kutta公式為：

(9)

式中：K1和Ki分別表示時間段開始時的斜率和每一時間段中點的斜率；ci,αi,βij均為待定系數(shù)；h為步長?？梢钥闯觯匠?9)含有m個函數(shù)值，所以稱為m階Runge-Kutta公式。本文使用的RK78算法表示采用七階-八階變步長Runge-Kutta算法，它用七階方法提供候選解，八階方法控制誤差，是一種自適應步長的常微分方程數(shù)值解法，其整體截斷誤差為O(h8) 。其變步長的原理如下：

對于一步h1對應其誤差估計Δ1，而h0對應其誤差估計Δ0，兩者關系為：

(10)

針對兩個誤差估計值，如果Δ1>Δ0，可以根據(jù)方程(10)確定所要減小步長的大??；如果Δ1<Δ0，可以根據(jù)方程(10)確定能夠在安全范圍內(nèi)增加步長的大小。

針對微小推力軌道外推問題，根據(jù)上述Runge-Kutta公式，使用動力學模型(8)帶入到七階Runge-Kutta公式(9)中。在此基礎上，設置合適的變步長，利用八階Runge-Kutta公式設置截斷誤差，通過積分和迭代，對衛(wèi)星軌道進行外推。

4 基于UKF的推力在軌標定算法

微小推力由于在軌保持時間較長，為了保證微小推力在軌工作的有效性，需要對推力進行在軌標定。本文采用UKF對微小推力進行高精度推力在軌標定。

4.1 UKF的基本原理

無跡卡爾曼濾波(Unscented Kalman Filter，UKF)，是無跡變換(Unscented Transform，UT)與標準卡爾曼濾波體系的結(jié)合，通過無跡變換使非線性系統(tǒng)方程適用于線性假設下的標準卡爾曼體系，該算法的核心思想是采用UT變換，利用一組Sigma采樣點來描述隨機變量的高斯分布，然后通過非線性函數(shù)的傳遞，再利用加權(quán)統(tǒng)計線性回歸技術(shù)來近似非線性函數(shù)的后驗均值和方差[27]。

UT變換的特點是對非線性函數(shù)的概率密度分布進行近似，而不是對非線性函數(shù)進行近似，即使系統(tǒng)模型復雜，也不增加算法實現(xiàn)的難度；而且所得到的非線性函數(shù)的統(tǒng)計量的準確性可以達到三階；除此之外，它不需要計算雅可比矩陣，可以處理不可導非線性函數(shù)。

算法模型如下：

設非線性系統(tǒng)：

(11)

式中：Xk為nx維的系統(tǒng)狀態(tài)向量；Zk為nz維的系統(tǒng)觀測向量；vk為系統(tǒng)噪聲；協(xié)方差矩陣為Pv；nk為觀測噪聲；協(xié)方差矩陣為Pn；式中：vk、nk都是高斯白噪聲，且互不相關。算法具體計算步驟如下：

1)選定濾波初值。

(12)

2)計算k-1時刻的2n+1個σ(Sigma)樣本點，構(gòu)造(2n+1)維向量χi為σ向量：

(13)

λ=α2(n+k)-n

(14)

3)狀態(tài)更新，計算k時刻的一步預測模型值。

(15)

4)計算k時刻的一步預測增廣樣本點。

(16)

5)觀測更新，計算k時刻的觀測變量的一步預測均值、方差和協(xié)方差。

(17)

6)計算增益矩陣。

(18)

7)計算濾波值。

(19)

4.2 在軌標定算法

利用UKF對微小推力的大小及方向進行在軌標定的過程如下：

(1)動力學模型改進

將笛卡爾動力學模型進行狀態(tài)擴維，即包含帶估計的參數(shù)。

(20)

(2)基本方程

利用上一部分的軌道遞推算法，對一個測量周期內(nèi)的航天器狀態(tài)進行預測，寫出預測方程；然后根據(jù)UKF基本原理中的公式分別寫出狀態(tài)估值方程、濾波增益方程、一步預測均方誤差方程和估計均方誤差方程等。

(3)更新

1)設置位置誤差和速度誤差，設計GNSS數(shù)值生成模擬器。

2)融合GNSS的測量值，在UKF濾波器內(nèi)對狀態(tài)預測進行更新。實時估計狀態(tài)x、m、aF、φ、θ。

5 仿真校驗

5.1 軌道遞推仿真結(jié)果

針對微小推力軌道外推問題，使用變步長龍格庫塔RK78算法進行衛(wèi)星軌道預測。初始參數(shù)設置如下：

1)設置初始軌道為圓軌道，參數(shù)軌道半徑為755 3 km，軌道傾角86.5°，偏心率0，升交點赤經(jīng)40°，近地點幅角60°，真近點角30°，地球引力常數(shù)3.986×1014m3/s2。通過軌道要素變換，可以求解衛(wèi)星的初始位置和速度，如表1所示。

表1 初始條件

2)設置衛(wèi)星質(zhì)量為194 kg，微小推力12 mN，平均分配到三軸推力[6.928 2,6.928 2,6.928 2]TmN，推力加速度為[3.571 24,3.571 24,3.571 24]T×10-5m/s2。此外，考慮到微小推力較小，衛(wèi)星推進器的質(zhì)量流量忽略不計。

3)設置太陽光壓攝動參數(shù)：輻射壓力系數(shù)為1.21，衛(wèi)星表面積為3.88 m2，地球到太陽的距離AU=149 597 870 000 m，單位AU距離下的太陽光壓力P=4.56×10-6N/m2。

4)設置太陽引力攝動太陽引力常數(shù)為：1.327 124 38×1020m3/s2，月球引力常數(shù)為：4.902 793 455×1012m3/s2。

利用C++軟件編寫整個推演過程的程序，和STK軟件模擬衛(wèi)星軌道推演24 h，比較不同情況下所編寫的C++程序和STK的結(jié)果對比。根據(jù)對比結(jié)果，有如下結(jié)論：

1)由表2～表6可以看出，在完整考慮4種主要軌道攝動和微小推力作用下，本文所設計的C++程序和STK模擬結(jié)果具有較高的一致性，24 h軌道遞推的位置誤差均在22 m以下，速度誤差小于0.012 m/s，幾乎可以忽略不計。

2)由表2和表3對比可以看出，在僅考慮微小推力情況時(不考慮其他軌道攝動)，本文所設計的C++程序和STK模擬的結(jié)果具有高度一致性，兩種方法得到的24 h軌道遞推的位置差異小于0.3 m，速度差異小于0.001 m/s。而由表4和表2對比可以看出，當考慮4種主要軌道攝動時，兩種方法得到的24 h軌道遞推的位置差異增加至30 m量級，速度差異增加至0.01 m/s量級。

表2 無攝動+無推力結(jié)果對比

表3 無攝動+12 mN推力結(jié)果對比

表4 考慮4種主要軌道攝動+無推力結(jié)果對比

3)由表3與表6對比可以看出，地球形狀攝動力是引起本文所設計的C++程序和STK模擬結(jié)果差異的最大因素，位置差異為22 m左右，主要原因可能是采用的地球引力模型之間存在輕微差異，其他三種主要攝動力并不會引起另種方法之間的較大位置差異。

表5 考慮4種主要軌道攝動+12 mN推力結(jié)果對比

表6 僅考慮地球扁率攝動+12 mN推力結(jié)果對比

總之，盡管本文所設計的C++程序?qū)Φ厍蛐螤顢z動力的處理與STK內(nèi)置的地球形狀攝動力有所差異，但本文所設計的C++程序仍可實現(xiàn)在微小推力作用下的高精度軌道預測。

為了驗證其他影響因素對所設計的C++程序的影響，本文分別針對軌道高度、推力大小和推演時間對程序的結(jié)果進行分析。

圖3和圖4分別給出了不同軌道高度下(200 km～1 500 km)C++程序和STK模擬之間的位置誤差和速度誤差，其中除軌道高度外，其余初始參數(shù)均與本節(jié)第一部分一致，程序考慮12 mN推力且無攝動下的結(jié)果。通過結(jié)果可以看出，不同軌道高度下，程序的位置誤差均在0.33 m以下，速度誤差均在4×10-4m/s以下，這說明本文的程序在不同軌道高度下均可以保證高精度軌道遞推。

圖3 不同軌道高度下位置誤差Fig.3 Position errors at different orbital altitudes

圖4 不同軌道高度下速度誤差Fig.4 Velocity errors at different orbital altitudes

圖5和圖6分別給出了不同推力大小下(0 mN-20 mN)C++程序和STK模擬之間的位置誤差和速度誤差，其中除推力大小外，其余初始參數(shù)均與本節(jié)第一部分一致，程序考慮存在推力且無攝動下的結(jié)果。通過結(jié)果可以看出，不同推力大小下，程序的位置誤差均在0.6 m以下，速度誤差均在6×10-4m/s以下，這說明本文的程序在不同推力大小下均可以保證高精度軌道遞推。

圖5 不同推力大小下位置誤差Fig.5 Position errors with different thrusts

圖6 不同推力大小下速度誤差Fig.6 Velocity errors with different thrusts

圖7和圖8分別給出了不同推演時間下(1～14 d)C++程序和STK模擬之間的位置誤差和速度誤差，其中除推演時間外，其余初始參數(shù)均與本節(jié)第一部分一致，程序考慮存在4種攝動和12 mN推力下的結(jié)果。通過結(jié)果可以看出，隨著推演時間的增加位置誤差和速度誤差都會變大，14 d內(nèi)的位置誤差在140 m以下，位置誤差在0.35 m/s以下，這說明本文的程序隨著推演時間的增大，誤差會累積增大，但還能保證正常的精度要求。

圖7 不同推演時間下位置誤差Fig.7 Position errors with different propagation time

圖8 不同推演時間下速度誤差Fig.8 Velocity errors with different propagation time

5.2 在軌標定仿真結(jié)果

根據(jù)上文的算法，利用C/C++對微小推力大小和其本體坐標系下的方向在3 h內(nèi)進行在軌標定的仿真驗證。

針對本文中的微小推力在軌標定問題，使用無跡卡爾曼濾波器(UKF)進行衛(wèi)星軌道估計和推力標定，通過預測和更新，完成對微小推力的標定。衛(wèi)星軌道遞推的初始參數(shù)與5.1節(jié)第一部分初始參數(shù)一致。

利用所編寫的C++程序在3 h時間間隔內(nèi)對衛(wèi)星微小推力進行推力在軌標定，初始位置各分量誤差為1 000 m，速度各分量誤差為0.05 m/s，微小推力各方向初始誤差為0.000 01 m/s2，測量精度為1 m，濾波器更新時間間隔為10 s。表7展示了一次軌道標定結(jié)果，可以看出位置誤差、速度誤差和推力誤差均收斂，推力精度占比收斂至0.429 4%。

表7 推力在軌標定3 h結(jié)果

如圖9～圖12以圖示的形式給出了UKF推力在軌標定的結(jié)果，圖9展示了3 h內(nèi)的位置估計誤差，圖10為速度誤差，圖11為推力誤差，圖12為推力精度占比，由上圖可以看出，各誤差和占比均收斂。

圖9 推力在軌標定的位置誤差Fig.9 Position errors of on-orbit thrust calibration

圖10 推力在軌標定的速度誤差Fig.10 Velocity errors of on-orbit thrust calibration

圖11 推力在軌標定的推力誤差Fig.11 Thrust errors of on-orbit thrust calibration

圖12 推力在軌標定的推力精度占比Fig.12 Thrust accuracy ratio of on-orbit thrust calibration

此外，通過改變C++程序中的位置誤差、速度誤差和位置誤差，該程序仍能保證各誤差和占比均收斂。

6 結(jié)論

對于微小推力下的小衛(wèi)星軌道遞推和推力在軌標定問題，本文提出了基于變步長龍格庫塔的軌道遞推算法，利用高精度的軌道動力學模型，能夠?qū)π⌒l(wèi)星軌道進行高精度的推演。在此基礎上，設計基于UKF推力標定算法，對小衛(wèi)星施加的微小推力進行在軌標定。仿真實例中對軌道遞推算法與STK進行比較，位置誤差均在22 m以下，速度誤差均在0.02 m/s以下；同時推力在軌標定各誤差均收斂，且精度得到保證，仿真數(shù)據(jù)說明了本文所設計算法的可行性。本文設計的算法基于C++程序，相較于其他程序算法具有更快的計算速率，有利于實際任務中實現(xiàn)星上的在線計算。在下一步的研究中，致力于消除地球形狀攝動力對算法帶來的誤差。