應用GNSS定位數據的微小衛星自主導航方案設計

劉燎 吳愛國 陶釗榕 孫華苗 謝成清

(1 哈爾濱工業大學(深圳)，廣東深圳 518055) (2 深圳航天東方紅海特衛星有限公司，廣東深圳 518064)

隨著低成本全球衛星導航系統(GNSS)接收機(包括美國的GPS及中國的“北斗”導航系統)的應用[1]，在微小衛星上利用GNSS接收機進行實時軌道確定進而提高其自主能力，已成為一種發展趨勢[2]。一個完整的自主導航系統應包括衛星實時定軌和軌道預報功能，使衛星能自主完成星上任務規劃。通常，軌道預報的計算方法有數值方法和解析方法2種。數值方法計算精度較高，適用于實時的狀態遞推；缺點是耗費大量機時，且難以完成預報。解析方法利用解析模型快速地預報任意時刻衛星的位置、速度；缺點是精度較低[3]。目前，應用GPS定位數據的衛星實時軌道確定主要采用卡爾曼(Kalman)濾波[4]及基于軌道動力學模型的軌道確定技術。但是，對于采用這種軌道確定技術的導航系統，其結果很難用于長期的在軌星歷預報(例如幾個軌道周期后甚至幾天后的軌道預報)，且往往需要占用大量的星載計算機機時，因而不適用于長期軌道預報[5]。文獻[6]中提出了一種采用簡化的動力學模型和一種嵌套插值算法的積分器，但是也需要高性能的星載計算機，不適合應用于微小衛星的軌道預報。文獻[7]中對基于GPS定位數據的衛星軌道預報算法進行研究，但并未對衛星實時軌道確定方法進行研究。

針對目前常用的近地圓軌道，本文設計一種應用GNSS定位數據的微小衛星自主導航方案，可進行自主軌道確定和預報。自主導航模塊根據GNSS定位數據，采用一種改進的數值積分方法，可進行在軌實時軌道確定，且計算量相對較小，可解決采用解析動力學模型精度較低的問題。同時，GNSS定位數據發送給星歷模型生成模塊，利用基于UD分解的Kalman濾波技術進行星歷參數的在軌估計，在GNSS接收機故障或者定位數據無效的情況下，利用星歷參數進行軌道預報，提高衛星的自主能力。

1 導航方案設計

本文提出的自主導航方案系統框圖見圖1。首先，GNSS接收機生成的定位數據發送給自主導航模塊，采用基于4階龍格-庫塔法和Richardson外插法及擴展Kalman濾波方法進行導航運算[7]，對衛星在慣性坐標系下的位置和速度矢量進行估計，根據狀態估計值生成飛行軌跡插值多項式，再將該插值多項式發送到軌道預報模塊。軌道預報模塊根據自動導航模塊發送的軌跡插值多項式，按制導和控制周期要求生成當前時刻的軌道狀態數據。考慮到備份的要求，GNSS接收機生成的定位數據同時發給星歷模型生成器，后者根據解析軌道模型，采用Bierman UD分解形式的Kalman濾波技術對相關星歷模型參數進行在軌估計[8]。該導航方案可對任意時間間隔的衛星星歷進行預報，不必按步長對軌道進行積分，其計算量相對較小，對星載計算機的性能要求不高，可滿足微小衛星的使用要求。當GNSS接收機出現故障后，估計得到的軌道參數可作為備份[9]，確保衛星在軌飛行時能獲得軌道信息。

圖1 導航方案系統框圖Fig.1 Navigation scheme system block diagram

在上述導航方案中，需要解決自主導航模塊中對于軌道狀態的計算和星歷模型生成器中對于星歷參數的在軌估計問題。

1.1 軌道狀態計算

自主導航模塊根據GNSS定位數據，采用基于4階龍格-庫塔法和Richardson外插法及擴展Kalman濾波進行導航運算，對衛星在慣性坐標系下的位置和速度矢量進行估計。

1)軌道動力學模式選擇

分析表明，在低軌近圓軌道上，太陽、月球引力攝動及太陽輻射壓力的影響可忽略不計，大氣阻力的影響相對較大，但由于測量信息頻度較高，可隨時更正，因此仿真分析中可忽略[8]。地球引力場模型采用4×4階的JGM-3模型已能滿足精度要求，考慮余量，本文中地球引力場模型取6×6階。

2)運動方程及數值積分算法

衛星所受地球引力攝動加速度為

(1)

(2)

式中：地心距r=‖r‖，r為衛星在J2000.0地心赤道坐標系中的位置矢量矩陣形式；λ和φ分別為衛星地心經度和緯度；GM和RE分別為地球引力場系數和地球赤道平均半徑；非球形諧函數Cnm和Snm表示地球非球形引力場；Pnm為n階m次的締合勒讓德多項式。

(3)

式中：衛星的狀態方程f可表示為狀態變量y和時間t的函數。

本文采用一種基于4階龍格-庫塔法和Richardson外插法的數值積分算法(RK4R)。在[t,t+2h](t和h分別為積分時間和積分步長)內進行3次RK4數值積分，積分步長分別為[t,t+h]，[t+h,t+2h]，[t,t+2h]。完成數值積分計算后，根據Richardson外插法，對軌道狀態變量y(t+h)和y(t+2h)進行修正，從而提高計算精度。

3)軌道參數濾波算法

(4)

Kalman增益矩陣表達式為

(5)

式中：R為測量誤差方差矩陣；雅克比矩陣Gi=?zi/?y，zi為衛星狀態的測量值。

狀態測量值的更新及其相應的協方差矩陣為

(6)

式中：g為測量方程。

(7)

1.2 星歷參數在軌估計

本節利用GNSS定位數據對星歷參數進行在軌估計，得到星歷參數可進行軌道計算和預報，作為軌道參數信息的備份和進行軌道預報等相關信息計算。圖2為星歷參數估計流程。

圖2 星歷參數估計流程Fig.2 Ephemeris parameter estimation flow

1.2.1 算法說明

本文算法的計算過程可描述為：根據標志判斷是否需要初始化；若此時需要進行初始化，根據GNSS定位數據進行濾波器初始化；若本次計算過程中初始化已完成，則進行Kalman濾波計算。利用GNSS定位數據對歷元時刻軌道平根數濾波器進行初始化，該初始化過程可描述為：將GNSS測量得到的位置和速度矢量轉換為J2000坐標系下；將J2000坐標系下的位置和速度矢量轉換為瞬時軌道根數；通過迭代過程計算平均根數，并給出對應的星歷模型參數；根據提供的參數計算狀態誤差協方差矩陣，并進行UD分解；輸出GNSS時刻對應的星歷模型參數以及初始U0陣和D0陣。

1.2.2 濾波算法

濾波流程如圖3所示。

圖3 濾波流程Fig.3 Flow of filter

利用基于Bierman UD分解Kalman濾波算法對星歷參數進行濾波，得到更新的星歷參數以及協方差矩陣對應的U陣和D陣，主要計算步驟如下。

(1)根據濾波得到的星歷模型參數X-以及其對應的歷元時刻t0，進行軌道預報，得到當前GNSS定位數據時刻tGNSS對應的衛星位置矢量。

(2)計算殘差。若殘差大于給定的容限，則判斷為野值，退出當前濾波過程；否則，執行步驟(3)。

(8)

式中：RZ為繞J2000坐標系Z軸的旋轉矩陣；S為tGNSS時刻對應的格林威治恒星時角。

(9)

(10)

2 仿真驗證

本文采用內部不帶導航算法的GNSS接收機。這種類型的接收機完成對GPS衛星信號的接收后，僅進行簡單的處理，直接輸出定位數據。其優點是初始化時間短、可靠性較高；缺點是定位解的速度精度較差。用MATLAB軟件對本文方案進行數學仿真驗證，仿真條件如表1所示。

對于仿真條件，不同軌道高度的濾波(定位)精度如圖4所示。對于500 km軌道高度，軌道濾波精度在40 m以內；600 km軌道高度，濾波精度在30 m以內；700 km軌道高度，濾波精度在25 m以內。

在GNSS數據無效的情況下，采用軌道遞推的形式，濾波時間和軌道預報精度如圖5所示。可以看出：濾波時間為2天，預報時間5天以內，軌道預報精度基本在15 km以內，可滿足中低精度的微小衛星任務需求。其中，濾波時間2天情況下，對于3種軌道高度，軌道預報精度時間曲線如圖6所示。

表1 仿真條件

圖4 3種軌道高度的濾波精度Fig.4 Filter accuracy for three orbit altitudes

圖5 3種軌道高度的預報精度Fig.5 Prediction accuracy for three orbit altitudes

圖6 濾波2天的預報精度Fig.6 Prediction accuracy for 2-day filter time

從以上仿真結果分析，700 km濾波精度在25 m以內，可滿足一般微小衛星0.1°指向精度的要求。預報精度基本在15 km以內，也可滿足微小衛星1°指向精度的要求和對進出地影區和地面站時間的預報精度要求。隨著軌道高度的增加，相同濾波時間和預報時間下，預報精度會提高，主要是由于隨著軌道高度的增加，氣動力對軌道的干擾在減小，1000 km以上影響基本可以忽略不計。因此，本文提出的導航方案可以滿足大部分微小衛星的中等精度指向精度和軌道預報任務要求，同時可顯著提高微小衛星的自主性，降低地面測控系統的壓力和成本，具有良好的工程應用前景。

3 結束語

本文針對微小衛星長期自主運行需求，研究了應用GNSS定位數據的微小衛星自主導航方案。該方案主要根據GNSS定位數據提高衛星在軌的自主性：根據GNSS定位數據進行濾波生成軌跡插值多項式(初始化時間短)，可生成當前時刻的軌道狀態數據和進行軌道預報。同時，又考慮到備份的要求，根據解析軌道模型來估計平均軌道參數，當GNSS接收機出現故障后，確保衛星在軌飛行時能獲得軌道信息。在軌道動力學選擇過程中，本文地球引力場模型取6×6階，在選取高階引力場模型的基礎上，在軌實時軌道確定能否達到10 m以內和預報精度能夠達到10 km以內，并且如何減少對計算機計算量的需求，可作為下一步研究的內容。