葉紅軍，潘 峰,李 笛

(1.衛(wèi)星導(dǎo)航系統(tǒng)與裝備技術(shù)國家重點實驗室，河北 石家莊050081; 2.北京衛(wèi)星導(dǎo)航中心,北京100094)

0 引言

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

1 導(dǎo)航衛(wèi)星軌道誤差因素

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

1.1 重力場模型分析

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

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

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

1.2 大氣阻尼模型分析

大氣阻尼模型目前難以建立十分精確的模型。迄今為止雖然已經(jīng)建立了很多大氣模型，但各模型之間的差異比較大，在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模型提出了一個合理的大氣密度模型的描述，因此廣泛的應(yīng)用于精密定軌和軌道外推研究中。J71模型考慮了周日和27天太陽效應(yīng)，根據(jù)大氣密度是角度ψd(衛(wèi)星和周日峰之間的角度)的函數(shù)，通過對衛(wèi)星觀測推導(dǎo)而來。在高度h大氣密度被表示為：

(1)

式中，

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

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

(2)

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

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

2 高穩(wěn)定軌道積分方法

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

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

2.1 Runge-Kutta算法

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

(3)

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

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

(4)

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

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

2.2 Adams算法

對初值問題，其解為：

(5)

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

(6)

式中，F(xiàn)n為函數(shù)在tn處的取值；h為步長；▽kF為函數(shù)F的n階向后差分。

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

2.3 Cowell算法

對初值問題

(7)

其解可以寫為：

(8)

式中，X為衛(wèi)星的位置向量?；蛘哒f，函數(shù)F只是問題的函數(shù)而與速度無關(guān)。

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

(9)

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

3 算法選擇與流程設(shè)計

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

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

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

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

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

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

每種積分器既有其優(yōu)點，也有其缺點。在實際的應(yīng)用實踐過程中，要根據(jù)具體問題來選擇適當?shù)姆e分方法，既要考慮積分的精度、穩(wěn)定性和收斂性，也要考慮計算的速度，同時還要考慮程序設(shè)計的復(fù)雜性。

表2 數(shù)值積分方法特征

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

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

圖1 衛(wèi)星導(dǎo)航模擬器軌道積分器流程

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

4 仿真結(jié)果及分析

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

導(dǎo)航星座軌道仿真位置差異如表5所示。北斗衛(wèi)星星座軌道仿真三維位置差異均方根如圖2所示。

表3 導(dǎo)航星座軌道仿真積分方法設(shè)置

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

表4 導(dǎo)航星座軌道仿真力學(xué)模型設(shè)置

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

表5 導(dǎo)航星座軌道仿真位置差異均方根

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

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

通過圖2可以看出在相同參數(shù)設(shè)置下，通過該方法仿真的結(jié)果與STK軟件軌道數(shù)據(jù)偏差遠小于0.1m，可以滿足衛(wèi)星導(dǎo)航模擬器數(shù)學(xué)仿真軟件的要求并基于此產(chǎn)生相應(yīng)的導(dǎo)航星歷和電文。

5 結(jié)束語

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

[1] 冉承其.“北斗”衛(wèi)星導(dǎo)航系統(tǒng)建設(shè)與發(fā)展[J].國際太空，2013(10)：11-15．

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

[3] 李雋.衛(wèi)星導(dǎo)航信號模擬器體系結(jié)構(gòu)分析[J].無線電工程，2006，36(8)：30-39.

[4] 葉紅軍.多模式衛(wèi)星導(dǎo)航模擬器設(shè)計與實現(xiàn)[J].無線電工程，2014，44(7)：43-46．

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

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

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

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

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

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

[11] 劉子令，姚志成.衛(wèi)星/慣性組合導(dǎo)航信號仿真器設(shè)計

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

[12] 郁聰沖，邊少鋒.現(xiàn)階段北斗衛(wèi)星導(dǎo)航系統(tǒng)可用性分析[J].海洋測繪.2012，32(5)：74-76．

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

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

[15] 周兵.北斗衛(wèi)星導(dǎo)航系統(tǒng)發(fā)展現(xiàn)狀與建設(shè)構(gòu)想[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] 高孝杰，張晶晶，高微，等.基于網(wǎng)絡(luò)模式的北斗高精度定位數(shù)據(jù)播發(fā)[J].計算機工程，2017，43(10)：1-4.

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

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

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

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