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

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

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

1 引言

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

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

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

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

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

2.1 問題描述

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

2.2 動力學(xué)模型

首先建立三個坐標(biāo)系[24]，如圖1所示。

圖1 坐標(biāo)系示意圖Fig.1 Illustration of coordinate frames

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

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

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

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

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

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

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

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

(1)

(2)

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

(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分別表示微小推力加速度，地球形狀攝動加速度，其他攝動加速度在慣性坐標(biāo)系下的三個坐標(biāo)軸的分量。

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

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

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

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

(9)

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

對于一步h1對應(yīng)其誤差估計(jì)Δ1，而h0對應(yīng)其誤差估計(jì)Δ0，兩者關(guān)系為：

(10)

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

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

4 基于UKF的推力在軌標(biāo)定算法

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

4.1 UKF的基本原理

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

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

算法模型如下：

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

(11)

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

1)選定濾波初值。

(12)

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

(13)

λ=α2(n+k)-n

(14)

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

(15)

4)計(jì)算k時刻的一步預(yù)測增廣樣本點(diǎn)。

(16)

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

(17)

6)計(jì)算增益矩陣。

(18)

7)計(jì)算濾波值。

(19)

4.2 在軌標(biāo)定算法

利用UKF對微小推力的大小及方向進(jìn)行在軌標(biāo)定的過程如下：

(1)動力學(xué)模型改進(jìn)

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

(20)

(2)基本方程

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

(3)更新

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

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

5 仿真校驗(yàn)

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

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

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

表1 初始條件

2)設(shè)置衛(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)星推進(jìn)器的質(zhì)量流量忽略不計(jì)。

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

4)設(shè)置太陽引力攝動太陽引力常數(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種主要軌道攝動和微小推力作用下，本文所設(shè)計(jì)的C++程序和STK模擬結(jié)果具有較高的一致性，24 h軌道遞推的位置誤差均在22 m以下，速度誤差小于0.012 m/s，幾乎可以忽略不計(jì)。

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

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

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

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

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

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

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

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

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

圖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 在軌標(biāo)定仿真結(jié)果

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

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

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

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

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

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

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

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

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

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

6 結(jié)論

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