衛星導航模擬器星座軌道外推方法研究

葉紅軍，潘 峰,李 笛

(1.衛星導航系統與裝備技術國家重點實驗室，河北 石家莊050081; 2.北京衛星導航中心,北京100094)

0 引言

衛星導航模擬器可以通過仿真的方式生成各種動態條件接收機天線端的導航信號，可以更加直觀地認識軌道的空間分布、運行周期以及衛星系統在地面上的DOP值分布等，結合信號處理算法可以十分逼真地對北斗、GPS等導航系統的觀測值信號進行仿真，生成模擬的載波和偽距觀測值，實現對接收機的調試及測試[1-4]。逼真的模擬導航星座運行及空間信號的各種特征，是模擬器開展接收機測試的前提[5-7]，其中長時間高精度的實現運行軌道的模擬是模擬器產生導航電文的基礎，北斗導航系統由于我國國土的特點在星座設計上是一種多軌道面混合組網的方式，對其高精度的導航模擬軌道外推與GPS星座相比更加困難，同時在模擬器運行中還需要綜合考慮工程的可實現性，針對此問題，開展針對北斗混合星座的長時間精密軌道遞推研究，并實現精度與計算量的統籌，提高工程可實現性。

1 導航衛星軌道誤差因素

導航衛星在空間軌道運行中將受到各種攝動力的影響，使得其運行軌道變得很復雜，以致初始時刻的位置和速度(或6個軌道參數)不再是理想的二體問題情況下的常數，而是時間的復雜函數。衛星受力一般分為兩大類:中心力與非中心力[8-10]。保守攝動力包括地球非球體引力位攝動、N體攝動、因日月引力引起的地球固體潮攝動及海潮攝動、大氣潮攝動、地球自轉形變攝動以及因相對論效應引起的攝動等。非保守攝動力包括大氣阻力攝動、太陽直射輻射壓攝動和地球反照輻射壓攝動等[11]。

1.1 重力場模型分析

衛星受地球非球形引力是否精確取決于地球重力場模型的精確程度。確定引力場模型的方法主要有3種：通過對衛星軌道的跟蹤觀測、利用陸地重力測定法和利用高度計數據確定引力場模型。

在近40多年的時間里，經過研究發現對于高度計數據誤差影響最大的就是衛星軌道誤差，對于軌道誤差影響最大的就是地球重力場模型的精確程度，因此在重力場模型上投入了廣泛的、深入的研究，從1959年第一次精密確定J2項系數以來，到目前先后有Kozai、Itzsak、Kaula、Rapp、Gaposchkin、GEM-9、GEM-10B、GRIM3-L1、GRIM4-2S2、GEM-T3、JGM-3、EGM96S。JGM-3對于精密定軌已經是一個非常精確的重力場模型，是70階70次的模型[12]。盡管JGM3相對于精密定軌已經足夠，人們又發展了新的重力場模型，利用了40顆衛星的跟蹤數據，采用GPS和TDRSS衛星對衛星跟蹤系統，使引力場模型又有了進一步的提高，獲得了EGM96S(70×70)和EGM96(360×360)模型。

因此，可以通過選擇JGM-3、EGM96S和EGM96模型來計算衛星軌道，能達到比較理想的軌道外推效果。

1.2 大氣阻尼模型分析

大氣阻尼模型目前難以建立十分精確的模型。迄今為止雖然已經建立了很多大氣模型，但各模型之間的差異比較大，在200 km的高度上，各模型典型密度差為20%，在較高的高度上相差更高。目前大氣模型有：Harris-Priester、Jacchia-71、Jacchia-Roberts、Jacchia-Lineberry、Jacchia-Gill、Jacchia-77、Jacchia-Lafontaine、MSIS-77、MSIS-86、TD-88和DTM[13]。

J71模型提出了一個合理的大氣密度模型的描述，因此廣泛的應用于精密定軌和軌道外推研究中。J71模型考慮了周日和27天太陽效應，根據大氣密度是角度ψd(衛星和周日峰之間的角度)的函數，通過對衛星觀測推導而來。在高度h大氣密度被表示為：

(1)

式中，

logeρ′(h)=-16.021-1.985×10-3h+6.383e-0.026h；

F10.7為10.7cm太陽輻射流量；

(2)

除了J71模型的廣泛使用外，各種其他密度模型也依然有其應用價值，其中也不乏簡單的、容易執行的和精巧的理論的模型，為了比較這些模型，表1是以J71模型為基準計算的各模型平均大氣密度、最大大氣密度之差以及在CPU計算時間上的差異。

表1 各種大氣密度模型與J71的比較CPU

2 高穩定軌道積分方法

在衛星的精密軌道外推中，積分器的設計是一項基礎工作。它的作用體現在2個方面：如果建立了衛星的力模型，則可以獲得衛星的運動方程。通過對衛星的運動方程進行積分就可以獲得衛星的運動狀態(位置和速度向量)；同時積分器要完成變分方程的積分，獲得衛星狀態對軌道參數和力模型參數的偏導數用于軌道外推。程序運行的大量時間是在積分過程，因此除了動力學模型需要很精確外，外推計算中也需要高精度的積分器。

目前積分器所采用的比較成熟的算法有以下幾種：① Runge-Kutta算法；② Adams算法；③Cowell算法。

2.1 Runge-Kutta算法

Runge-Kutta方法可以用來解如下的初值問題：

(3)

式中，X0為變量X的初值；F為變量t和X的函數。

如果步長為h，則Runge-Kutta算法可以用來計算變量X在t0+h的值X(t0+h)。重復這種步驟就可以獲得一系列的解X(t0+h)，X(t0+2h)，X(t0+3h)，……，X(t0+nh)，其中，n是一個整數。則利用Taylor級數可以將X(tn+nh)在點tn處展開如下：

(4)

本質上，Runge-Kutta算法是對同階的Taylor級數展開的逼近。將此算法向前推進一步就需要處多次計算函數F的值，在多數實際應用中，這是非常耗時間的。因此，稱Runge-Kutta算法為單步算法，該方法一般用作多步法的起步算法。

積分的誤差將取決于步長的大小和函數F的性質。為了保證軌道積分具有合理的精度，在計算中對積分步長進行自適應控制是一種可取的辦法。由于軌道的周期性，因此，在軌道積分中僅僅需要對一些特殊的周期采用步長的控制。Press于1992年提出了一種新的控制方法，每次積分進行2次計算，一次是采用全步長進行，然后以半步長進行2次；將這2種方法獲得的結果進行比較就可以判斷是否需要對步長進行改變。

2.2 Adams算法

對初值問題，其解為：

(5)

Adams算法使用Newton向后差分計算來擬合函數F，表達式如下：

(6)

式中，Fn為函數在tn處的取值；h為步長；▽kF為函數F的n階向后差分。

在多數情況下Adams-Moulton方法要比Adams-Bashforth方法的精度更高；但是，在計算函數值Fn+1前必須要有Xn+1，因此，在使用Adams-Moulton方法時要使用遞歸的方法。一種比較簡單的方法是利用Adams-Bashforth方法計算Xn+1的近似值，然后利用Xn+1的近似值計算Fn+1的近似值。經驗表明，上述方法一般可以達到足夠高的精度。

2.3 Cowell算法

對初值問題

(7)

其解可以寫為：

(8)

式中，X為衛星的位置向量。或者說，函數F只是問題的函數而與速度無關。

與Adams-Bashforth方法類似，函數F由式(9)擬合：

(9)

由于使用了Fn+1來擬合函數，因此，Cowell方法可以達到比Stormer方法更高的精度。同樣，計算Fn+1需要Xn+1上的值，因此，解決此問題需要進行迭代。可以使用Stormer方法計算Xn+1的近似值，然后利用該近似值計算Fn+1從而開始迭代過程。

3 算法選擇與流程設計

為了解決定軌問題中的初值問題，上面分析了3種積分方法。Runge-Kutta方法是一種單步算法，不同階數的Runge-Kutta方法具有不同的計算公式。即使在同階的Runge-Kutta方法中，其系數的確定也不是唯一的，在每次向前推進積分時，需要多次計算函數值，這在計算效率上是非常不利的。Runge-Kutta方法一個很重要的特點在于它是自起步的，因此，一般Runge-Kutta方法用作多步算法的起步算法。

Adams算法是一種多步算法。由于其系數間存在簡單的對應關系，因此，如果需要提高該算法的階數則很方便。但由于該算法是不能自起步的，在起步后，每次積分都只需要計算一次函數值，因此，特別對那些具有復雜關系的初值問題具有很高的計算效率。如果想達到更高的精度，可以通過在迭代過程中綜合使用Adams-Bashforth和Adams-Moulton方法獲得[14-16]。

Cowell方法也是一種多步算法，其階數也可以很方便地提高，但該方法也需要其他方法提供起步值。研究表明，Cowell方法具有比同階Adams方法更高的精度，但是，根據本文推導過程可以發現,該方法只適用于右函數僅僅是時間和衛星位置向量函數的特殊情況。但是，像空氣阻力就與衛星的速度有關。因此，Cowell方法只能在部分問題中使用，其他綜合了Cowell算法的方法也具有這種特點。

通過上面的討論，很明顯需要將衛星運動過程中受到的力進行分類:一類力只與時間和衛星的位置有關，其他的力為另一類。第1類可以使用Cowell方法進行積分，而第2類就可以采用Adams方法。2種方法的起步初始值均由Runge-Kutta方法提供。

至于如何選擇積分步長和積分器的階數則取決于軌道遞推所需要的精度。一般，在軟件設計中，這2個參數會設置成可調節的，這樣，在幾次迭代運行后就可以確定一個比較合適的值。在實用中，一般采用8階Runge-Kutta方法和12階Adams方法就足夠了。但并不是步長越小越好，也不是階數越高越好。

表2對幾種常用在衛星軌道遞推中的數值積分器特征進行了對比。根據積分器計算下一步點用到前面一個還是多個步點的信息，可把積分器分為單步法和多步法；根據積分過程中積分步長是固定不變還是自動調控，可把積分器分為固定步長法和變步長法；根據多步法積分過程中是否用到差分算子，可把多步法積分器分為求和型和非求和型(也叫一般型)；根據積分器是對一階還是二階微分方程積分，可把積分器分為一次型和二次型。一般來說，多步法優于單步法，因為隨著數值方法的階數增加，后者計算右函數的次數將增加，而前者基本保持不變，對于軌道遞推這樣的動力學問題，方程右函數的計算是十分復雜，數值積分耗時主要表現在右函數的計算上。變步長法優于固定步長法，因為變步長法對橢圓型軌道更有效；求和型優于非求和型，因為求和型更能減少計算的省略誤差。二次型優于一次型，因為在相同的精度要求下，對二階微分方程求解所允許的步長比對一階微分方程求解時使用的步長大的多。

每種積分器既有其優點，也有其缺點。在實際的應用實踐過程中，要根據具體問題來選擇適當的積分方法，既要考慮積分的精度、穩定性和收斂性，也要考慮計算的速度，同時還要考慮程序設計的復雜性。

表2 數值積分方法特征

數值積分方法單步法/多步法變步長/固定步長求和型/非求和型一次型/二次型RK(Runge-Kutta)方法單步法固定步長NA一次型RKF(Runge-Kutta-Fehlberg)方法單步法變步長NA一次型Adams方法多步法固定步長非求和型一次型Cowell方法多步法固定步長非求和型二次型

由于北斗衛星星座為由GEO、IGSO和MEO衛星組成的混合星座，因此在實際導航模擬器軌道外推算法上采用RKF和Adams相結合的方法來進行實現。軌道積分器主要完成對GPS/BDS/GLONASS/GALILEO衛星的軌道積分外推工作。數據仿真技術流程圖如圖1所示。

圖1 衛星導航模擬器軌道積分器流程

該模塊輸入的是衛星的初始歷元以及力模型的基本參數，結果為積分得到的參考軌道以及相應的對各動力參數的偏導數。該積分器起步算法采用的是Runge-Kutta-fehlberg方法(RKF)，它是一種嵌套的Runge-Kutta方法，容易在計算機上實現，而且變步長方便，能保證所需要的精度，穩定度也較好。起步后，采用基于Adams顯式Adams-Bashfort公式和隱式Adams-Moulton公式，以及Cowell公式的預報校正多步線性積分算法進行數值積分。以上算法非常成熟，而且基本上都可以滿足精度要求。但是，由于低軌衛星受力狀況比較復雜，某些動力模型參數需要分段估計，因此需要積分器能處理力模型參數突變的能力。

4 仿真結果及分析

衛星導航模擬器軌道仿真模塊與STK10.0的高精度軌道預報模塊(HPOP)積分設置如表3所示，力學模型設置如表4所示，采用相同的地球定向參數文件。考慮的力學模型有地球引力、地球固體潮汐影響、日月引力、太陽輻射壓力與相對論效應影響，其余影響較小的力學模型未考慮在內。

導航星座軌道仿真位置差異如表5所示。北斗衛星星座軌道仿真三維位置差異均方根如圖2所示。

表3 導航星座軌道仿真積分方法設置

積分設置項現有軌道仿真STK高精度軌道預報積分開始時刻UTC時間2010年1月1日0時0分0秒積分初始值地心固定坐標系積分坐標系地心慣性坐標系積分方法RKF6(7)啟動Adams11階預估校正RKF7(8)積分步長固定60s相對誤差控制最小步長1s最大步長86400s輸出步長60s積分時長24h,86400s輸出坐標系地心固定坐標系

表4 導航星座軌道仿真力學模型設置

參量模型設置地球引力WGS84_EGM96引力位,21階21次地球固體潮汐影響考慮日月影響,截斷到地球引力位的階數和次數,包含時間依賴項[17-19]三體引力考慮太陽、月亮,采用JPL的DE405星歷太陽輻射壓力球形模型,太陽輻射壓力系數1.0,衛星面積質量比值0.02m2/kg,考慮太陽光線從太陽到地球的傳播時間,采用雙錐地影模型[20-22]相對論效應影響考慮

表5 導航星座軌道仿真位置差異均方根

導航星座位置差異均方根/m限差/mCOMPASS0.013GPS0.018GLONASS0.026Galileo0.0170.1

圖2 北斗衛星星座軌道仿真三維位置差異均方根

通過圖2可以看出在相同參數設置下，通過該方法仿真的結果與STK軟件軌道數據偏差遠小于0.1m，可以滿足衛星導航模擬器數學仿真軟件的要求并基于此產生相應的導航星歷和電文。

5 結束語

通過調研分析國外已有成熟的精密定軌及軌道外推技術，結合衛星導航模擬器工程化實現需求，深入研究了衛星導航星座軌道外推仿真的數學模型和算法，完成了各種主要誤差因素對軌道仿真精度的影響分析，確定了衛星導航星座軌道外推的方法和軟件方案，完成了導航系統星軌道仿真軟件流程設計與實現，仿真成果與STK軟件軌道數據偏差小于0.1m，達到了較高的精度，輸出結果包括衛星軌道坐標(根數)、衛星運動速度，可實現48h的軌道外推，解決了以往軌道外推隨著時間逐漸偏差加大的問題，對生成連續的衛星導航星歷支持衛星導航接收機開展長時間測試驗證具有直接的支持作用。

[1] 冉承其.“北斗”衛星導航系統建設與發展[J].國際太空，2013(10)：11-15．

[2] 譚述森.衛星導航定位工程(第2版)[M].北京：國防工業出版社，2010．

[3] 李雋.衛星導航信號模擬器體系結構分析[J].無線電工程，2006，36(8)：30-39.

[4] 葉紅軍.多模式衛星導航模擬器設計與實現[J].無線電工程，2014，44(7)：43-46．

[5] 黃勇，胡小工，王小亞，等.中高軌衛星廣播星歷精度分析[J].天文學進展，2006，24(1)：81-87．

[6] 高玉東，郗曉寧，王威.導航衛星廣播星歷擬合改進算法設計[J].國防科技大學學報，2007，29(5)：18-21．

[7] 樓益棟，劉萬科，張小紅.GPS衛星星歷精度分析[J].測繪信息與工程，2003，28(6)：4-6.

[8] 趙娜，董崢.衛星星座運行管理方法研究[J].無線電工程，2010，40(6)：62-64.

[9] 劉季.北斗GEO衛星軌道算法研究[J].測繪地理信息，2012，37(5)：33-36.

[10] 吳靜，常青，吳今培，等.高動態GPS信號模擬器衛星星歷產生方法研究[J].無線電工程，2004，34(5)：42-60.

[11] 劉子令，姚志成.衛星/慣性組合導航信號仿真器設計

[J].無線電工程，2014，44(7)：39-50．

[12] 郁聰沖，邊少鋒.現階段北斗衛星導航系統可用性分析[J].海洋測繪.2012，32(5)：74-76．

[13] 高為廣，蘇牡丹，李軍正，等.北斗衛星導航系統試運行服務性能評估[J].武漢大學學報信息科學版，2012，37(11)：352-355．

[14] 雷浩，廉保旺，何偉，等.STK北斗二代衛星導航系統在亞太地區DOP值仿真分析[J].火力與指揮控制，2014，39(6)：52-55．

[15] 周兵.北斗衛星導航系統發展現狀與建設構想[J].無線電工程，2016，46(4)：1-4．

[16]WANGMengli，SUNGuangfu，WANGFeixue，etal.WeightedGeometricDilutionofPrecision’sAnalysisforMixedConstellationNavigationSystem[J].ChineseSpaceScienceandTechnology，2007，10(5)：50-56．

[17]CAIHongliang，LIXing.WeaknessAnalysisofNavigationConstellationBasedonServiceVailability[J].CSNC2015，2015162(3)：150-156．

[18] 高孝杰，張晶晶，高微，等.基于網絡模式的北斗高精度定位數據播發[J].計算機工程，2017，43(10)：1-4.

[19] 李飛琦，鮑泓，潘峰，等.智能車導航中的路口軌跡生成策略[J].計算機工程，2017，43(8)：1-7.

[20] 黃建生，王曉玲，王敬艷，等.GPS導航定位設備測試技術研究[J].電子技術與軟件工程，2013(6)：36-37.

[21] 羅大成，劉巖，劉延飛，等.星間鏈路技術的研究現狀與發展趨勢[J].電訊技術，2014，54(7)：1016-1024.

[22] 胡方強，呂濤，包亞萍.改進的自適應Kalman濾波在SINS/GPS組合導航中的應用[J].計算機工程與應用，2017，53(4)：1-6.