張國亭 王宏 黃勇 鄭世貴 王振興

(1 國防科技大學(xué) 空天科學(xué)學(xué)院,長沙 410073)(2 北京跟蹤與通信技術(shù)研究所, 北京 100094) (3 中國科學(xué)院上海天文臺,上海 200030)(4 北京空間飛行器總體設(shè)計(jì)部, 北京 100094)

我國的航天事業(yè)發(fā)展迅猛,對衛(wèi)星的定軌精度要求也越來越高。在衛(wèi)星動力學(xué)定軌過程中,有兩種難以精確模制的非引力攝動:大氣阻力攝動和太陽光壓攝動。對于中高軌道衛(wèi)星,由于沒有大氣阻力攝動且對地球非球形引力攝動也不敏感,因此太陽光壓攝動成為繼地球引力,日月引力之后量級最大的攝動,特別是導(dǎo)航衛(wèi)星動力學(xué)定軌建模誤差最大的攝動力,也是研究的熱點(diǎn)方向[1]。以我國北斗導(dǎo)航衛(wèi)星為例,相比于中地球軌道(MEO)衛(wèi)星,地球靜止軌道(GEO)和傾斜地球同步軌道(IGSO)衛(wèi)星軌道高度較高,受到的太陽輻射壓影響更大[2]。

太陽輻射壓對衛(wèi)星產(chǎn)生的攝動加速度與太陽輻射強(qiáng)度、衛(wèi)星相對于太陽的姿態(tài)和位置、衛(wèi)星的面質(zhì)比、幾何形狀及其表面材質(zhì)的物理特性等因素有關(guān)。由于太陽輻射強(qiáng)度的變化、衛(wèi)星姿態(tài)控制的偏差以及衛(wèi)星表面材料的老化等因素,太陽輻射壓攝動較難進(jìn)行模制,它已經(jīng)成為高軌衛(wèi)星動力學(xué)定軌最主要的誤差源。

太陽輻射壓模型根據(jù)建模原理方法的不同可以分為物理分析模型和經(jīng)驗(yàn)?zāi)Ｐ汀Ｎ锢矸治瞿Ｐ褪侵富谛l(wèi)星的結(jié)構(gòu)和表面的光學(xué)特性所建立的模型。經(jīng)驗(yàn)?zāi)Ｐ褪侵笖M合衛(wèi)星的在軌數(shù)據(jù)得到的模型。此外,還有一類模型,建模時(shí)既使用了衛(wèi)星的基本信息,又使用了衛(wèi)星的在軌數(shù)據(jù),稱半經(jīng)驗(yàn)?zāi)Ｐ?。這幾種模型都曾或仍然被廣泛使用,各有優(yōu)缺點(diǎn),對于衛(wèi)星太陽輻射壓建模具有重要的借鑒意義[3-6]。

目前,導(dǎo)航衛(wèi)星精密定軌廣泛采用擴(kuò)展經(jīng)驗(yàn)光壓模型(ECOM)等經(jīng)驗(yàn)?zāi)Ｐ蚚7-8],不僅用于全球定位系統(tǒng)(GPS)的定軌計(jì)算,GPS/GLONASS的聯(lián)合定軌和北斗導(dǎo)航衛(wèi)星定軌也常會采用該模型來模制太陽輻射壓攝動[9],并得到較好的結(jié)果。ECOM模型使用9個(gè)參數(shù)分別對三軸上太陽輻射壓攝動的常數(shù)項(xiàng)和周期項(xiàng)進(jìn)行擬合吸收,不僅吸收太陽輻射壓建模的缺陷,而且還能吸收一些其他未模制的力的影響,提高定軌的精度。另外,ECOM模型不使用任何先驗(yàn)?zāi)Ｐ?直接估計(jì)三軸上的太陽輻射壓攝動,也能取得較好的定軌結(jié)果。然而,想要取得較高的定軌精度,必須要有足夠的在軌觀測數(shù)據(jù),因此ECOM等經(jīng)驗(yàn)?zāi)Ｐ筒贿m用在軌運(yùn)行初期的衛(wèi)星的定軌工作。ECOM參數(shù)模型完全是一個(gè)經(jīng)驗(yàn)?zāi)Ｐ?需要大量的在軌數(shù)據(jù)來擬合,待估參數(shù)的增多會影響法方程性能。另外,多少時(shí)間擬合一組參數(shù),也是一個(gè)值得研究的問題。

經(jīng)驗(yàn)?zāi)Ｐ碗m然可以取得較高的定軌精度,但是其估計(jì)的參數(shù)沒有實(shí)際的物理意義,無法揭示太陽輻射壓攝動對衛(wèi)星軌道影響的物理機(jī)制。尤其對于在軌運(yùn)行初期的衛(wèi)星,沒有足夠的在軌觀測數(shù)據(jù),很難通過經(jīng)驗(yàn)?zāi)Ｐ瞳@取較好的定軌結(jié)果。因此,建立高精度的分析型太陽輻射壓模型是非常有必要的。

本文以我國北斗導(dǎo)航衛(wèi)星作為研究對象,研究了分析型光壓模型在北斗GEO和IGSO等高軌衛(wèi)星定軌中的應(yīng)用,首先介紹了分析型光壓的建模原理,然后選取了一顆北斗GEO衛(wèi)星和一顆IGSO衛(wèi)星進(jìn)行了分析,比較了分析型光壓結(jié)果和精密定軌結(jié)果的差異,并利用分析型光壓模型進(jìn)行精密定軌,精度約為分米級。

1 分析型光壓建模原理

考慮一個(gè)“無限小”面元ds:

(1)只考慮光子的吸收和鏡反射,作用于面元上的光壓力為(如圖1所示)

dF=dF1+dF2

(1)

式中:dF1和dF2分別為面元ds受到的吸收作用力和鏡反射作用力,有

(2)

dF2=-ρs|dF1|dr

(3)

式中:n為面元ds的法向單位矢量;d為入射方向的單位矢量;dr為反射方向的單位矢量;ρs為面元ds的鏡反射系數(shù)。

F2的方向符合鏡反射規(guī)律,但其大小與面元的鏡反射系數(shù)有關(guān)。于是光壓力又可表達(dá)為

(4)

式中:θ為入射方向單位矢量d與面元ds的法向單位矢量n的夾角;dscosθ為面元ds在垂直入射方向的投影。

將dr用已知方向d和n來表達(dá),有

(5)

圖1 面元上的光壓力示意Fig.1 SRP on the face

(2)只考慮光子漫反射,漫反射遵從cos法則,該法則假設(shè)反射方向分布正比于反射方向與面法向夾角的cos值。漫反射關(guān)于面法向?qū)ΨQ,平均反射速度方向即為面法向方向,作用于面元上的光壓力為

(6)

式中:ρe為漫反射系數(shù)。

綜合光子的吸收、鏡反射和漫反射,作用于面元ds上的光壓力為

(7)

該模型包含了光子在材料表面的各種行為,但要準(zhǔn)確獲得表面的光學(xué)性質(zhì)比較困難,而且在空間環(huán)境的長期作用下,材料表面的光學(xué)性質(zhì)會發(fā)生變化。

本文采用的是基于蒙特卡洛射線追蹤算法計(jì)算光壓,蒙特卡洛法又稱隨機(jī)模擬法或統(tǒng)計(jì)試驗(yàn)法,是一種依據(jù)統(tǒng)計(jì)抽樣理論,利用計(jì)算機(jī)研究隨機(jī)變量的數(shù)值計(jì)算方法?；驹頌樵跔顟B(tài)變量的概率分布已知且相互獨(dú)立的條件下,隨機(jī)產(chǎn)生符合狀態(tài)變量概率分布的隨機(jī)數(shù),然后以這些隨機(jī)數(shù)為基礎(chǔ)生成服從某一分布的狀態(tài)函數(shù),進(jìn)行后續(xù)分析。

射線追蹤算法是追蹤射線路徑上與射線相交的所有面,計(jì)算出射線與面元的所有交點(diǎn),然后只考慮距離射線出發(fā)點(diǎn)最近的交點(diǎn)。

射線追蹤一般步驟如下:

(1)確定一個(gè)平行六面體封閉盒,該盒包含全部幾何模型;

(2)把封閉盒分成許多等尺寸的方塊;

(3)對于每一個(gè)方塊,列出落在方塊中的所有面,這些面至少有一個(gè)點(diǎn)落在方塊中;

(4)依照射線前進(jìn)方向,依次計(jì)算射線路徑上的方塊;

(5)檢查射線與當(dāng)前方塊中的所有面是否有交點(diǎn),如果存在交點(diǎn)且交點(diǎn)位于當(dāng)前方塊中,執(zhí)行(6),否則返回(4)計(jì)算下一個(gè)方塊;

(6)找出距離射線出發(fā)點(diǎn)最近的交點(diǎn)完成一次光壓計(jì)算。

離散空間中的射線追蹤順序示意見圖2。

可把太陽入射光看成一定數(shù)量光線的集合,每條光線相互獨(dú)立。由于航天器距離太陽非常遠(yuǎn),可認(rèn)為所有光線為平行光線。如果光線數(shù)量足夠就能夠模擬太陽入射光產(chǎn)生的壓力效應(yīng)。計(jì)算光線的一次撞擊,同時(shí)考慮光子與表面的鏡反射、漫反射與吸收效應(yīng)。采用射線追蹤算法的計(jì)算誤差約為射線數(shù)開平方分之一。如10000條射線,則光壓力的誤差為1%。基于射線追蹤算法的光壓分析程序流程如圖3所示。

圖2 射線追蹤順序示意Fig.2 Ray tracing sequence

圖3 基于射線追蹤算法的光壓分析程序流程Fig.3 Steps of the analytical SRP model program based on ray tracing approach

影響光壓分析精度的因素有衛(wèi)星幾何模型、表面材料光學(xué)參數(shù)、太陽輻射強(qiáng)度、光壓模型、數(shù)值計(jì)算方法(本文方法為射線數(shù))等。

2 北斗衛(wèi)星分析型光壓模型精度評估

隨著伽利略(Galileo)衛(wèi)星導(dǎo)航系統(tǒng)和北斗衛(wèi)星導(dǎo)航系統(tǒng)的逐步建設(shè)和完善,以及GPS系統(tǒng)和GLONASS系統(tǒng)逐步現(xiàn)代化推出更多新的民用信號,國際GNSS服務(wù)組織(IGS)于2012年著手在全球范圍內(nèi)構(gòu)建多模GNSS系統(tǒng)試驗(yàn)網(wǎng)絡(luò)(Multi-GNSS EXperiment,MGEX),實(shí)現(xiàn)對GPS、GLONASS、Galileo、北斗和其他衛(wèi)星導(dǎo)航系統(tǒng)的跟蹤和監(jiān)視,為新興衛(wèi)星導(dǎo)航系統(tǒng)的高精度數(shù)據(jù)處理提供了數(shù)據(jù)條件。基于MGEX全球監(jiān)測站的北斗IGSO和MEO衛(wèi)星定軌結(jié)果徑向(R)誤差優(yōu)于10cm,三維軌道誤差約為幾十厘米,可以用于精密單點(diǎn)定位服務(wù)。但是全球跟蹤網(wǎng)的增加,GEO衛(wèi)星的軌道精度不會得到明顯提升,約在米級。MGEX的定軌結(jié)果基于ECOM等經(jīng)驗(yàn)光壓模型[10]。

北斗GEO-6于2012年發(fā)射,星下點(diǎn)經(jīng)度為80.3°E,北斗IGSO-5于2011年發(fā)射,星下點(diǎn)經(jīng)度(赤道)為96°E。從MGEX網(wǎng)站下載了2016年3月3日至4月2日約1個(gè)月GEO-6和IGSO-5的精密星歷(武漢大學(xué)分析中心結(jié)果),sp3c格式的衛(wèi)星軌道為衛(wèi)星地固系(WGS-84)下的X、Y、Z軸坐標(biāo)值。將精密星歷的X-Y-Z作為觀測量對軌道進(jìn)行逐天定軌,定軌參數(shù)包括初始軌道根數(shù)和經(jīng)驗(yàn)型的ECOM5光壓模型參數(shù),統(tǒng)計(jì)定軌精度,并輸出太陽輻射壓攝動加速度。

圖4為GEO-6和IGSO-5的定軌精度,GEO-6的定軌精度大部分時(shí)間優(yōu)于0.2m,IGSO-5的定軌精度優(yōu)于0.1m(時(shí)間起點(diǎn)為2016年3月3日)。

按照第1節(jié)分析型光壓建模的原理,結(jié)合北斗衛(wèi)星的幾何尺寸和表面材料特征,以及探測器和太陽的相對位置等信息,可以計(jì)算得到北斗衛(wèi)星的太陽光壓攝動力。圖5給出了利用分析型光壓模型得到的GEO-6和IGSO-5衛(wèi)星所受到的太陽輻射壓攝動加速度,坐標(biāo)系為RTN(軌道的徑向R、橫向T和法向N)。前述精密定軌輸出的太陽輻射壓攝動加速度時(shí)間序列與之類似,在RTN三個(gè)方向分別約為10-7、10-7和10-8m/s2量級,N方向受力比其他兩個(gè)方向小一個(gè)數(shù)量級。頻譜分析結(jié)果表明:GEO-6和IGSO-5衛(wèi)星太陽輻射壓攝動的主要周期為24h,和軌道周期一致。

將該值和前述精密定軌輸出的太陽輻射壓加速度進(jìn)行比較,差異結(jié)果如圖6和圖7所示。

圖4 基于經(jīng)驗(yàn)光壓模型的GEO-6和IGSO-5衛(wèi)星的定軌精度Fig.4 Orbital accuracy of GEO-6 and IGSO-5 satellites based on ECOM

圖5 分析型模型得到的太陽輻射壓攝動加速度Fig.5 Analytical model of SRP perturbation acceleration

圖6 GEO-6衛(wèi)星的分析型光壓結(jié)果和精密定軌結(jié)果差異Fig.6 Differences between analytical SRP results and precision orbit determination results of GEO-6 satellite

圖7 IGSO-5衛(wèi)星的分析型光壓結(jié)果和精密定軌結(jié)果差異Fig.7 Differences between analytical SRP results and precision orbit determination results of IGSO-5 satellite

從上述分析結(jié)果可以看出:精密定軌輸出的太陽輻射壓加速度和分析型光壓模型的結(jié)果比較,對GEO-6衛(wèi)星,R和T方向差異在10-8m/s2量級,N方向在10-9m/s2量級,相對誤差在10%左右。對IGSO-5衛(wèi)星,差異要小1/2左右,約為5%?？紤]到分析型光壓模型的建模過程和精密定軌無關(guān),而精密定軌輸出的太陽輻射壓加速度和精密定軌密切相關(guān),IGSO衛(wèi)星相比GEO衛(wèi)星,分析型光壓模型和精密定軌所采用的經(jīng)驗(yàn)?zāi)Ｐ头铣潭雀?這可能是由于IGSO衛(wèi)星的定軌精度更高引起的。

3 分析型光壓模型在導(dǎo)航衛(wèi)星精密定軌中的應(yīng)用

進(jìn)一步分析分析型光壓模型在北斗高軌衛(wèi)星精密定軌的應(yīng)用,一方面可以直接應(yīng)用分析模型輸出的太陽輻射壓加速度時(shí)間序列來代替定軌軟件中的光壓模型;另一方面考慮到分析模型可能還存在一定的誤差;可以對其增加若干調(diào)節(jié)參數(shù),類似于對星載加速度計(jì)數(shù)據(jù)的處理方式,如下

aci=ai+ki×ami

(8)

式中:i=1,2,3(對應(yīng)于RTN軸),aci為經(jīng)過標(biāo)校后的i方向的加速度值;ai為加速度常數(shù)偏差參數(shù);ki為i方向的尺度因子參數(shù);ami為觀測值。標(biāo)校參數(shù)在定軌時(shí)予以解算。

對GEO-6和IGSO-5衛(wèi)星,如果直接使用分析型光壓模型,以MGEX軌道作為參考軌道,定軌精度如圖8所示(時(shí)間起點(diǎn)為2016年3月3日),GEO-6位置誤差約在2.5m左右,IGSO-5位置誤差約1.5m左右。

圖8 直接使用分析型光壓模型定軌精度Fig.8 Orbit accuracy is determined directly using analytical SRP model

增加標(biāo)校參數(shù)(偏差和尺度因子),GEO-6和IGSO-5定軌精度如圖9所示(時(shí)間起點(diǎn)為2016年3月3日),GEO-6位置誤差小于1m,平均誤差約為0.5m。解算出來的尺度因子接近于1,偏差參數(shù)約在10-9～10-8m/s2量級。相應(yīng)IGSO-5的定軌精度比GEO-6的定軌精度均要高,增加標(biāo)校參數(shù)后位置誤差小于0.1m。

上述分析表明:分析型光壓模型可能存在一定的系統(tǒng)性偏差,可以通過建模的方式對系統(tǒng)誤差參數(shù)予以解算以提高定軌精度。從定軌結(jié)果看,基于分析型光壓模型的IGSO衛(wèi)星的定軌精度和常規(guī)經(jīng)驗(yàn)?zāi)Ｐ投ㄜ壘认喈?dāng),而GEO衛(wèi)星的定軌精度略低于常規(guī)經(jīng)驗(yàn)?zāi)Ｐ?初步驗(yàn)證了分析性光壓模型在北斗衛(wèi)星精密定軌中的應(yīng)用。

圖9 分析型光壓模型+標(biāo)校參數(shù)定軌精度Fig.9 Analytical SRP model and calibration parameters orbit accuracy

4 結(jié)束語

本文分析了分析型光壓模型在北斗高軌衛(wèi)星精密定軌中的應(yīng)用,對GEO-6和IGSO-5衛(wèi)星建立了分析型光壓模型,選取了2016年3月3日至4月2日約1個(gè)月的數(shù)據(jù),將分析型光壓模型結(jié)果和精密定軌輸出的太陽輻射壓加速度進(jìn)行了比較,對GEO-6衛(wèi)星,R和T方向差異在10-8m/s2量級,N方向在10-9m/s2量級,相對誤差在10%左右。對IGSO-5衛(wèi)星,差異要小1/2左右,約為5%。進(jìn)一步分析分析型光壓模型在北斗高軌衛(wèi)星精密定軌的應(yīng)用,直接使用分析型光壓模型,以MGEX軌道作為參考軌道,GEO-6位置誤差約在2.5m左右,IGSO-5位置誤差約1.5m左右;增加偏差和尺度因子標(biāo)校參數(shù)后,位置精度分別提高到0.5m和0.1m左右。本文分析結(jié)果表明:分析型光壓模型可以應(yīng)用于導(dǎo)航衛(wèi)星精密定軌,以北斗GEO和IGSO衛(wèi)星為例,可以實(shí)現(xiàn)分米級的定軌精度。該結(jié)果也可以應(yīng)用于其他有精密定軌需求的高軌衛(wèi)星,如果衛(wèi)星剛進(jìn)入軌道或者衛(wèi)星進(jìn)行軌道控制后,難以積累大量實(shí)測數(shù)據(jù)建立高精度的經(jīng)驗(yàn)光壓模型,在衛(wèi)星發(fā)射前建立分析型光壓模型,也是提高定軌精度的一種解決手段。