李 笛，孫騰達(dá)

(衛(wèi)星導(dǎo)航系統(tǒng)與裝備技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，河北 石家莊 050081)

0 引言

衛(wèi)星導(dǎo)航模擬器[1]可以通過(guò)數(shù)據(jù)仿真產(chǎn)生在接收機(jī)各種運(yùn)動(dòng)狀態(tài)下的導(dǎo)航信號(hào)，可以直觀地呈現(xiàn)衛(wèi)星導(dǎo)航系統(tǒng)衛(wèi)星的空間分布與可見(jiàn)性信息，用戶的位置、速度等軌跡信息和接收機(jī)接收觀測(cè)量信息等。因此，衛(wèi)星導(dǎo)航模擬器廣泛用于各種應(yīng)用場(chǎng)景中接收機(jī)的調(diào)試與測(cè)試[2]。低軌衛(wèi)星廣泛用于探測(cè)、通信和導(dǎo)航增強(qiáng)等[3]，發(fā)揮著越來(lái)越重要的作用，在低軌衛(wèi)星任務(wù)的地面仿真中，衛(wèi)星導(dǎo)航模擬器發(fā)揮著不可替代的作用。其中，低軌用戶軌跡長(zhǎng)時(shí)間、高精度仿真是模擬器正確仿真觀測(cè)量的基礎(chǔ)，其仿真精度直接影響低軌用戶數(shù)據(jù)仿真的逼真度。由于低軌衛(wèi)星受力情況非常復(fù)雜，實(shí)現(xiàn)其高精度軌跡仿真特別困難，同時(shí)在模擬器中運(yùn)行還需考慮工程可實(shí)現(xiàn)性，針對(duì)上述問(wèn)題，開(kāi)展低軌用戶高精度軌跡仿真方法研究，并兼顧精度與計(jì)算量，提升工程的可實(shí)現(xiàn)性。

1 低軌用戶動(dòng)力學(xué)模型

低軌用戶軌跡分布取決于空間復(fù)雜的力學(xué)環(huán)境[4]，故建立相對(duì)精確的動(dòng)力學(xué)模型[5]對(duì)低軌用戶軌跡仿真至關(guān)重要。低軌用戶在軌運(yùn)行期間，除了受地球中心引力[6]作用外，還將受到地球非球形引力[7]，太陽(yáng)、月球及其他天體引力的影響，以及太陽(yáng)光壓、大氣阻力和地球潮汐力等因素的影響。

1.1 重力場(chǎng)模型

在初步分析航天器運(yùn)動(dòng)規(guī)律時(shí)，把地球看作理想球體，其對(duì)航天器產(chǎn)生指向地心的、大小與航天器質(zhì)量m成正比而與地心與航天器距離r的平方成反比的引力F[8]：

(1)

由此得到地球引力常數(shù)：

μ=GM=3.986 000 441 5×1014m3/s2。

(2)

于是引力加速度為g=μ/r2，此情況下地球引力場(chǎng)稱為中心引力場(chǎng)。引力加速度為引力勢(shì)Uc的梯度，即Uc=μ/r。

地球并不是完美球體，具有扁平度、梨形和赤道變形等，因此，引力加速度是地心經(jīng)度φ、緯度λ和r的函數(shù)。則地球引力加速度矢量g為地球引力勢(shì)函數(shù)U的梯度：

g=gradU(r，φ，λ)。

(3)

目前引力勢(shì)函數(shù)的普遍表達(dá)[8]為：

([Cn，mcos(mλ)+Sn，msin(mλ)]}，

(4)

式中，RE=6 378.136 3 km為地球赤道半徑；Pn，m(x)為n階m級(jí)的Legendre函數(shù)。為克服各階數(shù)系數(shù)相差太大而引起的計(jì)算誤差，采用的范化式[6]為：

(5)

(6)

式中，

(7)

地球引力勢(shì)模型系數(shù)Jn(即-Cn，0)，Cn，m，Sn，m由測(cè)量得到。許多國(guó)際組織通過(guò)不同手段測(cè)量重力場(chǎng)系數(shù)并創(chuàng)建對(duì)應(yīng)的重力場(chǎng)模型，比較著名的重力場(chǎng)模型包括WGS系列、GEM系列和GRIM系列模型。WGS和GEM模型已經(jīng)合并為EGM96模型，IER2003推薦采用EGM96模型。

1.2 大氣阻力模型

大氣對(duì)低軌用戶的作用包括三部分[6]，與用戶對(duì)大氣流相對(duì)速度方向相反的大氣阻力、升力和副法線作用力，其中大氣阻力占支配地位，其余2種可以忽略。對(duì)低軌用戶而言，大氣阻力為其受到的最大的非引力性攝動(dòng)力。

低軌用戶的大氣阻力加速度[6]為：

(8)

式中，r為低軌用戶位置矢量；CD為大氣阻力系數(shù)；A為低軌用戶迎風(fēng)面；ρ為低軌用戶處的大氣密度；vr為用戶對(duì)大氣流的相對(duì)速度的大小；ev為相對(duì)速度的單位方向向量。

在大氣密度[9]計(jì)算方面，可采用目前最常用的3種經(jīng)驗(yàn)大氣密度模型，即Jacchia系列模型、DTM系列模型和MSIS系列模型。由于大氣變化機(jī)制復(fù)雜，任何經(jīng)驗(yàn)大氣密度模型都很難精確地刻畫(huà)大氣密度的實(shí)際變化。因此，為進(jìn)一步提高大氣模型精度，可利用實(shí)際數(shù)據(jù)對(duì)經(jīng)驗(yàn)?zāi)Ｐ瓦M(jìn)行進(jìn)一步校準(zhǔn)。同時(shí)，由于不同軌道高度大氣密度不同，所以針對(duì)不同軌道類型的低軌衛(wèi)星，還需研究不同經(jīng)驗(yàn)大氣密度模型下的軌道預(yù)報(bào)精度，探索選取最優(yōu)的大氣密度計(jì)算方法。

在大氣阻力系數(shù)[10]計(jì)算方面，一般而言，在低軌衛(wèi)星精密定軌中，常常將大氣阻力系數(shù)作為未知量與衛(wèi)星的運(yùn)動(dòng)狀態(tài)矢量一起解算，解算得到的大氣阻力系數(shù)會(huì)吸收模型誤差，在軌道預(yù)報(bào)中利用解算得到的大氣阻力系數(shù)來(lái)計(jì)算大氣阻力加速度，但大氣阻力系數(shù)的較難準(zhǔn)確地確定是制約解算大氣阻力加速度精度的一大原因。基于目前已有研究，可采用線性回歸分析建立大氣阻力系數(shù)的補(bǔ)償算法，還可基于空間環(huán)境數(shù)據(jù)和神經(jīng)網(wǎng)絡(luò)模型對(duì)空間目標(biāo)大氣阻力系數(shù)進(jìn)行修正。同樣，針對(duì)不同軌道類型的低軌衛(wèi)星還需研究不同大氣阻力系數(shù)算法下的軌道預(yù)報(bào)精度，探索選取最優(yōu)的大氣阻力系數(shù)計(jì)算方法。

在低軌衛(wèi)星迎風(fēng)面積[11]計(jì)算方面，需對(duì)低軌衛(wèi)星的幾何形狀進(jìn)行描述，可采用Cannon-Ball模型和Macro模型。Canno-Ball模型是假設(shè)低軌衛(wèi)星為一個(gè)球體，具有固定的形狀及單一的表面性質(zhì)，對(duì)幾何形狀較為簡(jiǎn)單的低軌衛(wèi)星十分適用。Macro模型也被稱為Box-Wing模型，它會(huì)對(duì)低軌衛(wèi)星每個(gè)面的幾何形狀都進(jìn)行描述，包括面積、單位外法向量以及該面的光學(xué)特性(對(duì)可見(jiàn)光和紅外光的反射和吸收)，特別當(dāng)?shù)蛙壭l(wèi)星在飛行過(guò)程中太陽(yáng)帆板發(fā)生翻轉(zhuǎn)時(shí)，Macro模型亦能做出幾何形狀的描述并反映到阻力計(jì)算中。

1.3 太陽(yáng)光壓模型

照射到衛(wèi)星上的太陽(yáng)光對(duì)衛(wèi)星產(chǎn)生入射作用力，衛(wèi)星吸收一部分太陽(yáng)光轉(zhuǎn)變成熱能和電能，另一部分被反射回去[12]。入射力和反射力統(tǒng)稱為太陽(yáng)輻射壓力，也稱太陽(yáng)光壓力。太陽(yáng)輻射壓力產(chǎn)生的加速度與太陽(yáng)光強(qiáng)度、垂直于入射方向的有效面積、表面反射率、與到太陽(yáng)的距離和光速有關(guān)。由于衛(wèi)星表面材料的老化、太陽(yáng)能量隨太陽(yáng)活動(dòng)的變化以及衛(wèi)星姿態(tài)控制的誤差等因素的影響，使得太陽(yáng)輻射壓攝動(dòng)也是難以精確模型的攝動(dòng)力。

直接的太陽(yáng)光壓力[6]可以表示為：

(9)

式中，P⊙為太陽(yáng)常數(shù)，等于4.560 4×10-6N/m2，其物理意義是完全吸收的物體在距太陽(yáng)一個(gè)天文單位(AU)處，單位面積上所受的輻射壓力；CR為太陽(yáng)光壓系數(shù)；A為垂直于衛(wèi)星與太陽(yáng)方向的截面積；ν為蝕因子(地影區(qū)ν=0，衛(wèi)星光照區(qū)ν=1，陰影區(qū)0<ν<1)；r⊙為太陽(yáng)的位置矢量。

1.4 N體攝動(dòng)模型

太陽(yáng)、月球和行星的引力攝動(dòng)可近似用點(diǎn)質(zhì)量來(lái)描述，在地心慣性系中，天體對(duì)衛(wèi)星的N體引力加速度為：

(10)

式中，GMi為質(zhì)點(diǎn)的引力常數(shù)；ri為質(zhì)點(diǎn)在地心慣性系中的位置矢量。

2 坐標(biāo)系及其轉(zhuǎn)換

為了減少牽連運(yùn)動(dòng)所引起的附加速度，使衛(wèi)星的運(yùn)動(dòng)方程相對(duì)簡(jiǎn)化，一般在慣性坐標(biāo)系[5]中建立和求解衛(wèi)星的運(yùn)動(dòng)方程。而地球?qū)πl(wèi)星的引力是在地固坐標(biāo)系[13]中定義的。此外，衛(wèi)星攝動(dòng)力分析計(jì)算要涉及到衛(wèi)星軌道坐標(biāo)系，把加速度轉(zhuǎn)換到慣性坐標(biāo)系或地固坐標(biāo)系要涉及衛(wèi)星固定坐標(biāo)系。因此要建立衛(wèi)星運(yùn)動(dòng)方程和參數(shù)估計(jì)的觀測(cè)方程，必須給出這些坐標(biāo)系統(tǒng)的定義及其相互轉(zhuǎn)換關(guān)系。

2.1 坐標(biāo)系統(tǒng)定義

2.1.1 協(xié)議慣性坐標(biāo)系

具有絕對(duì)意義的慣性坐標(biāo)系統(tǒng)是無(wú)法獲得的，人類只能根據(jù)科學(xué)任務(wù)的需要和一些約定建立相應(yīng)精度的慣性坐標(biāo)系統(tǒng)。通常使用的慣性坐標(biāo)系統(tǒng)是協(xié)議慣性系[14](CIS)，其具體實(shí)現(xiàn)就是國(guó)際協(xié)議天球參考框架(ICRF)。在CIS下，物體的運(yùn)動(dòng)可以以很高的精度滿足Newton力學(xué)定理，這對(duì)于動(dòng)力學(xué)低軌軌跡仿真的相關(guān)研究十分重要。

2.1.2 協(xié)議地球坐標(biāo)系

協(xié)議地球坐標(biāo)系地固坐標(biāo)系統(tǒng)[14](CTS)是為了描述地面觀測(cè)站位置所建立的與地球體聯(lián)結(jié)的坐標(biāo)系統(tǒng)，依表達(dá)方式又可分為直角坐標(biāo)系及大地坐標(biāo)系。CTS是使用起來(lái)最直觀的坐標(biāo)系統(tǒng)，它以特定的協(xié)議與地球固聯(lián)，并由一組全球參考站維持。

2.1.3 衛(wèi)星軌道坐標(biāo)系

為了便于對(duì)特定攝動(dòng)力進(jìn)行分析，將作用于低軌用戶的力分解到衛(wèi)星軌道坐標(biāo)系統(tǒng)[5](RTN)下進(jìn)行分析。此外，在RTN系統(tǒng)下，還便于對(duì)各項(xiàng)誤差進(jìn)行分析，如利用HL-SST模式研究恢復(fù)地球重力場(chǎng)時(shí)，徑向的軌道誤差占非常重要地位。

2.1.4 衛(wèi)星固定坐標(biāo)系

衛(wèi)星固定坐標(biāo)系[4](SBF)的主要作用是定義衛(wèi)星在慣性空間中的姿態(tài)，同時(shí)建立各相關(guān)傳感器坐標(biāo)系與慣性系的關(guān)系。如加速度計(jì)是固聯(lián)在衛(wèi)星質(zhì)心處，其軸系相應(yīng)地平行于衛(wèi)星固聯(lián)坐標(biāo)系，在轉(zhuǎn)換加速度計(jì)觀測(cè)量到CIS或CTS時(shí)都需要使用SBF 作為中間過(guò)渡坐標(biāo)系統(tǒng)。

2.2 坐標(biāo)系統(tǒng)轉(zhuǎn)換方法

2.2.1 協(xié)議地球坐標(biāo)系到協(xié)議慣性坐標(biāo)系的轉(zhuǎn)換

如果以RCIS表示某點(diǎn)在歷元J2000.0對(duì)應(yīng)的CIS中的坐標(biāo)，以RCTS表示CTS中的坐標(biāo)，則有[8]：

RCIS=P(t)N(t)S(t)Pm(t)RCTS，

(11)

式中，P(t)，N(t)，S(t)，Pm(t)分別為歲差矩陣、章動(dòng)矩陣、周日自轉(zhuǎn)矩陣和極移矩陣。極移矩陣將歷元平地球坐標(biāo)系，即CTS轉(zhuǎn)換到瞬時(shí)極地球坐標(biāo)系。地球自轉(zhuǎn)矩陣將瞬時(shí)極地球坐標(biāo)系轉(zhuǎn)換到真天球坐標(biāo)系。章動(dòng)矩陣將真平天球坐標(biāo)系轉(zhuǎn)換到瞬平天球坐標(biāo)系。歲差矩陣將瞬平天球坐標(biāo)系轉(zhuǎn)換到CIS[J2000.0]。根據(jù)IERS規(guī)范(McCarthy，1996)，有[8]：

(12)

(13)

(14)

(15)

2.2.2 協(xié)議慣性坐標(biāo)系到衛(wèi)星軌道坐標(biāo)系的轉(zhuǎn)換

RTN的3個(gè)坐標(biāo)軸在CIS下的單位向量[14]為：

(16)

即：

RCIS=MRRTN

(17)

式中，RRTN為質(zhì)點(diǎn)RTN坐標(biāo)系中的三維位置；M為從RTN到CIS的轉(zhuǎn)換矩陣：

M=

(18)

2.2.3 衛(wèi)星固定坐標(biāo)系到協(xié)議慣性坐標(biāo)系的轉(zhuǎn)換

從SBF到CIS的轉(zhuǎn)換可以通過(guò)太陽(yáng)和衛(wèi)星在協(xié)議慣性坐標(biāo)系中的位置矢量來(lái)實(shí)現(xiàn)。更高精度轉(zhuǎn)換一般需要用到衛(wèi)星的姿態(tài)數(shù)據(jù)。姿態(tài)數(shù)據(jù)給出的是由星體固定坐標(biāo)系到協(xié)議慣性系轉(zhuǎn)換的四元數(shù)(或稱歐拉參數(shù))q1，q2，q3，q4。由四元數(shù)表示的坐標(biāo)轉(zhuǎn)換矩陣[14]為：

Q=

(19)

所以，將SBF下的向量RSBF轉(zhuǎn)換到CIS下的向量RCIS的公式為：

RCIS=QRSBF。

(20)

3 算法流程設(shè)計(jì)

在利用低軌衛(wèi)星力學(xué)模型和坐標(biāo)系統(tǒng)轉(zhuǎn)換建立低軌衛(wèi)星軌跡信息(位置+速度)狀態(tài)方程的基礎(chǔ)上，兼顧軌跡準(zhǔn)確度與工程計(jì)算量的要求，利用Runge-Kutta算法完成數(shù)值積分，得到各歷元時(shí)刻的低軌衛(wèi)星軌跡信息。

3.1 Runge-Kutta數(shù)值積分算法

Runge-Kutta算法[6]是一種工程上廣泛應(yīng)用的高精度單步算法。已知方程導(dǎo)數(shù)與初值信息：

(21)

式中，y0為矢量y的初值；f(t，y)為矢量y關(guān)于變量t的一階偏導(dǎo)數(shù)。則根據(jù)Runge-Kutta算法：

y(tn+1)≈y(tn)+h·Φ=η(tn+1) ,

(22)

式中，η(tn+1)為y(tn+1)的近似解；Φ為增量函數(shù)，h=tn+1-tn。根據(jù)經(jīng)典RK4 Runge-Kutta算法[6]，增量方程為：

(23)

式中，

(24)

根據(jù)Oliver Montenbruck和Eberhard Gill的研究，RK4解的局部截?cái)嗾`差[6]有以下特征：

eRK4=|y(tn+1)-η(tn+1))≤const·h5，

(25)

且RK4解與四階泰勒展開(kāi)多項(xiàng)式精度相當(dāng)：

(26)

式中，y(tn)(i)，(i=2，3，4)為矢量y關(guān)于變量t的i階偏導(dǎo)數(shù)。

3.2 低軌用戶軌跡仿真算法設(shè)計(jì)

低軌用戶軌跡狀態(tài)量為其位置與速度[15]，即

(27)

式中，R為低軌用戶三維位置矢量；V為其三維速度矢量。對(duì)于任意已知R，V，可低軌衛(wèi)星受力模型與坐標(biāo)轉(zhuǎn)化方法得到其對(duì)應(yīng)的加速度α，則軌跡狀態(tài)量關(guān)于時(shí)間t的一階偏導(dǎo)數(shù)為：

(28)

在此基礎(chǔ)上，在給定初值R0，V0的前提下，利用RK4 Runge-Kutta算法即可計(jì)算得到各個(gè)歷元時(shí)刻的軌跡信息。

低軌軌跡仿真算法流程如圖1所示。

圖1 低軌軌跡仿真算法流程Fig.1 Flow chart of LEO track simulation algorithm

4 算法仿真及驗(yàn)證

為了測(cè)試低軌用戶軌跡仿真算法的逼真性與可靠性[16-17]，以STK10.0軟件仿真的低軌衛(wèi)星軌道數(shù)據(jù)作為理論參考值進(jìn)行測(cè)試驗(yàn)證。

4.1 軌跡參數(shù)設(shè)定

軌跡仿真開(kāi)始于UTC時(shí)間2010年1月1日0時(shí)0分0秒，持續(xù)86 400 s，仿真步長(zhǎng)為10 s，輸出地心地固坐標(biāo)系的位置速度。低軌用戶的六參數(shù)設(shè)置如表1所示。

表1 低軌用戶軌跡初值設(shè)置Tab.1 Initial value setting of LEO user track

地球重力場(chǎng)選用70階EGM96模型，不考慮固體潮影響；大氣密度計(jì)算選用NRLMSISE2000模型；太陽(yáng)光壓選用圓柱地影模型；考慮太陽(yáng)和月球的三體引力影響，選用JPL的DE405星歷計(jì)算日月位置；考慮相對(duì)論效應(yīng)。

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

本次仿真的軌跡與STK10.0軟件仿真軌跡比較結(jié)果顯示，經(jīng)過(guò)統(tǒng)計(jì)24 h后其三維位置誤差保持在10.0 m之內(nèi)，速度誤差保持在0.012 m/s之內(nèi)。

三維位置誤差分布如圖2所示，三維速度誤差分布如圖3所示。

圖2 三維位置誤差分布Fig.2 Error distribution graph of three-dimensional position

圖3 三維速度誤差分布Fig.3 Error distribution graph of three-dimensional velocity

通過(guò)上述測(cè)試結(jié)果分析，低軌用戶軌跡仿真精度能到達(dá)較高水平，典型攝動(dòng)力模型下逼近國(guó)外STK10.0軟件軌道仿真精度水平。因此，該仿真算法在逼真性和可靠性上完全可以滿足衛(wèi)星導(dǎo)航模擬器數(shù)學(xué)仿真的需要。

5 結(jié)束語(yǔ)

在調(diào)研國(guó)內(nèi)外已有軌跡仿真技術(shù)方案的基礎(chǔ)上，結(jié)合衛(wèi)星導(dǎo)航模擬器應(yīng)用實(shí)際需求，深入研究了低軌用戶軌跡仿真的數(shù)學(xué)模型，確定了低軌用戶軌跡仿真的算法，完成了低軌用戶軌跡仿真的流程設(shè)計(jì)與軟件實(shí)現(xiàn)。經(jīng)過(guò)仿真與分析，低軌用戶軌跡仿真軟件與STK10.0軟件的仿真結(jié)果偏差小于10 m，達(dá)到了很高的精度，提升了衛(wèi)星導(dǎo)航模擬器的性能，有力支持了低軌衛(wèi)星地面驗(yàn)證工作。