張國亭 王宏 黃勇 鄭世貴 王振興
(1 國防科技大學 空天科學學院,長沙 410073)(2 北京跟蹤與通信技術研究所, 北京 100094) (3 中國科學院上海天文臺,上海 200030)(4 北京空間飛行器總體設計部, 北京 100094)
我國的航天事業發展迅猛,對衛星的定軌精度要求也越來越高。在衛星動力學定軌過程中,有兩種難以精確模制的非引力攝動:大氣阻力攝動和太陽光壓攝動。對于中高軌道衛星,由于沒有大氣阻力攝動且對地球非球形引力攝動也不敏感,因此太陽光壓攝動成為繼地球引力,日月引力之后量級最大的攝動,特別是導航衛星動力學定軌建模誤差最大的攝動力,也是研究的熱點方向[1]。以我國北斗導航衛星為例,相比于中地球軌道(MEO)衛星,地球靜止軌道(GEO)和傾斜地球同步軌道(IGSO)衛星軌道高度較高,受到的太陽輻射壓影響更大[2]。
太陽輻射壓對衛星產生的攝動加速度與太陽輻射強度、衛星相對于太陽的姿態和位置、衛星的面質比、幾何形狀及其表面材質的物理特性等因素有關。由于太陽輻射強度的變化、衛星姿態控制的偏差以及衛星表面材料的老化等因素,太陽輻射壓攝動較難進行模制,它已經成為高軌衛星動力學定軌最主要的誤差源。
太陽輻射壓模型根據建模原理方法的不同可以分為物理分析模型和經驗模型。物理分析模型是指基于衛星的結構和表面的光學特性所建立的模型。經驗模型是指擬合衛星的在軌數據得到的模型。此外,還有一類模型,建模時既使用了衛星的基本信息,又使用了衛星的在軌數據,稱半經驗模型。這幾種模型都曾或仍然被廣泛使用,各有優缺點,對于衛星太陽輻射壓建模具有重要的借鑒意義[3-6]。
目前,導航衛星精密定軌廣泛采用擴展經驗光壓模型(ECOM)等經驗模型[7-8],不僅用于全球定位系統(GPS)的定軌計算,GPS/GLONASS的聯合定軌和北斗導航衛星定軌也常會采用該模型來模制太陽輻射壓攝動[9],并得到較好的結果。ECOM模型使用9個參數分別對三軸上太陽輻射壓攝動的常數項和周期項進行擬合吸收,不僅吸收太陽輻射壓建模的缺陷,而且還能吸收一些其他未模制的力的影響,提高定軌的精度。另外,ECOM模型不使用任何先驗模型,直接估計三軸上的太陽輻射壓攝動,也能取得較好的定軌結果。然而,想要取得較高的定軌精度,必須要有足夠的在軌觀測數據,因此ECOM等經驗模型不適用在軌運行初期的衛星的定軌工作。ECOM參數模型完全是一個經驗模型,需要大量的在軌數據來擬合,待估參數的增多會影響法方程性能。另外,多少時間擬合一組參數,也是一個值得研究的問題。
經驗模型雖然可以取得較高的定軌精度,但是其估計的參數沒有實際的物理意義,無法揭示太陽輻射壓攝動對衛星軌道影響的物理機制。尤其對于在軌運行初期的衛星,沒有足夠的在軌觀測數據,很難通過經驗模型獲取較好的定軌結果。因此,建立高精度的分析型太陽輻射壓模型是非常有必要的。
本文以我國北斗導航衛星作為研究對象,研究了分析型光壓模型在北斗GEO和IGSO等高軌衛星定軌中的應用,首先介紹了分析型光壓的建模原理,然后選取了一顆北斗GEO衛星和一顆IGSO衛星進行了分析,比較了分析型光壓結果和精密定軌結果的差異,并利用分析型光壓模型進行精密定軌,精度約為分米級。


考慮一個“無限小”面元ds:
(1)只考慮光子的吸收和鏡反射,作用于面元上的光壓力為(如圖1所示)
dF=dF1+dF2
(1)
式中:dF1和dF2分別為面元ds受到的吸收作用力和鏡反射作用力,有
(2)
dF2=-ρs|dF1|dr
(3)
式中:n為面元ds的法向單位矢量;d為入射方向的單位矢量;dr為反射方向的單位矢量;ρs為面元ds的鏡反射系數。
F2的方向符合鏡反射規律,但其大小與面元的鏡反射系數有關。于是光壓力又可表達為
(4)
式中:θ為入射方向單位矢量d與面元ds的法向單位矢量n的夾角;dscosθ為面元ds在垂直入射方向的投影。
將dr用已知方向d和n來表達,有
(5)

圖1 面元上的光壓力示意Fig.1 SRP on the face
(2)只考慮光子漫反射,漫反射遵從cos法則,該法則假設反射方向分布正比于反射方向與面法向夾角的cos值。漫反射關于面法向對稱,平均反射速度方向即為面法向方向,作用于面元上的光壓力為
(6)
式中:ρe為漫反射系數。
綜合光子的吸收、鏡反射和漫反射,作用于面元ds上的光壓力為
(7)
該模型包含了光子在材料表面的各種行為,但要準確獲得表面的光學性質比較困難,而且在空間環境的長期作用下,材料表面的光學性質會發生變化。
本文采用的是基于蒙特卡洛射線追蹤算法計算光壓,蒙特卡洛法又稱隨機模擬法或統計試驗法,是一種依據統計抽樣理論,利用計算機研究隨機變量的數值計算方法。基本原理為在狀態變量的概率分布已知且相互獨立的條件下,隨機產生符合狀態變量概率分布的隨機數,然后以這些隨機數為基礎生成服從某一分布的狀態函數,進行后續分析。
射線追蹤算法是追蹤射線路徑上與射線相交的所有面,計算出射線與面元的所有交點,然后只考慮距離射線出發點最近的交點。
射線追蹤一般步驟如下:
(1)確定一個平行六面體封閉盒,該盒包含全部幾何模型;
(2)把封閉盒分成許多等尺寸的方塊;
(3)對于每一個方塊,列出落在方塊中的所有面,這些面至少有一個點落在方塊中;
(4)依照射線前進方向,依次計算射線路徑上的方塊;
(5)檢查射線與當前方塊中的所有面是否有交點,如果存在交點且交點位于當前方塊中,執行(6),否則返回(4)計算下一個方塊;
(6)找出距離射線出發點最近的交點完成一次光壓計算。
離散空間中的射線追蹤順序示意見圖2。
可把太陽入射光看成一定數量光線的集合,每條光線相互獨立。由于航天器距離太陽非常遠,可認為所有光線為平行光線。如果光線數量足夠就能夠模擬太陽入射光產生的壓力效應。計算光線的一次撞擊,同時考慮光子與表面的鏡反射、漫反射與吸收效應。采用射線追蹤算法的計算誤差約為射線數開平方分之一。如10000條射線,則光壓力的誤差為1%?;谏渚€追蹤算法的光壓分析程序流程如圖3所示。

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

圖3 基于射線追蹤算法的光壓分析程序流程Fig.3 Steps of the analytical SRP model program based on ray tracing approach
影響光壓分析精度的因素有衛星幾何模型、表面材料光學參數、太陽輻射強度、光壓模型、數值計算方法(本文方法為射線數)等。
隨著伽利略(Galileo)衛星導航系統和北斗衛星導航系統的逐步建設和完善,以及GPS系統和GLONASS系統逐步現代化推出更多新的民用信號,國際GNSS服務組織(IGS)于2012年著手在全球范圍內構建多模GNSS系統試驗網絡(Multi-GNSS EXperiment,MGEX),實現對GPS、GLONASS、Galileo、北斗和其他衛星導航系統的跟蹤和監視,為新興衛星導航系統的高精度數據處理提供了數據條件。基于MGEX全球監測站的北斗IGSO和MEO衛星定軌結果徑向(R)誤差優于10cm,三維軌道誤差約為幾十厘米,可以用于精密單點定位服務。但是全球跟蹤網的增加,GEO衛星的軌道精度不會得到明顯提升,約在米級。MGEX的定軌結果基于ECOM等經驗光壓模型[10]。
北斗GEO-6于2012年發射,星下點經度為80.3°E,北斗IGSO-5于2011年發射,星下點經度(赤道)為96°E。從MGEX網站下載了2016年3月3日至4月2日約1個月GEO-6和IGSO-5的精密星歷(武漢大學分析中心結果),sp3c格式的衛星軌道為衛星地固系(WGS-84)下的X、Y、Z軸坐標值。將精密星歷的X-Y-Z作為觀測量對軌道進行逐天定軌,定軌參數包括初始軌道根數和經驗型的ECOM5光壓模型參數,統計定軌精度,并輸出太陽輻射壓攝動加速度。
圖4為GEO-6和IGSO-5的定軌精度,GEO-6的定軌精度大部分時間優于0.2m,IGSO-5的定軌精度優于0.1m(時間起點為2016年3月3日)。
按照第1節分析型光壓建模的原理,結合北斗衛星的幾何尺寸和表面材料特征,以及探測器和太陽的相對位置等信息,可以計算得到北斗衛星的太陽光壓攝動力。圖5給出了利用分析型光壓模型得到的GEO-6和IGSO-5衛星所受到的太陽輻射壓攝動加速度,坐標系為RTN(軌道的徑向R、橫向T和法向N)。前述精密定軌輸出的太陽輻射壓攝動加速度時間序列與之類似,在RTN三個方向分別約為10-7、10-7和10-8m/s2量級,N方向受力比其他兩個方向小一個數量級。頻譜分析結果表明:GEO-6和IGSO-5衛星太陽輻射壓攝動的主要周期為24h,和軌道周期一致。
將該值和前述精密定軌輸出的太陽輻射壓加速度進行比較,差異結果如圖6和圖7所示。

圖4 基于經驗光壓模型的GEO-6和IGSO-5衛星的定軌精度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衛星的分析型光壓結果和精密定軌結果差異Fig.6 Differences between analytical SRP results and precision orbit determination results of GEO-6 satellite

圖7 IGSO-5衛星的分析型光壓結果和精密定軌結果差異Fig.7 Differences between analytical SRP results and precision orbit determination results of IGSO-5 satellite
從上述分析結果可以看出:精密定軌輸出的太陽輻射壓加速度和分析型光壓模型的結果比較,對GEO-6衛星,R和T方向差異在10-8m/s2量級,N方向在10-9m/s2量級,相對誤差在10%左右。對IGSO-5衛星,差異要小1/2左右,約為5%??紤]到分析型光壓模型的建模過程和精密定軌無關,而精密定軌輸出的太陽輻射壓加速度和精密定軌密切相關,IGSO衛星相比GEO衛星,分析型光壓模型和精密定軌所采用的經驗模型符合程度更好,這可能是由于IGSO衛星的定軌精度更高引起的。
進一步分析分析型光壓模型在北斗高軌衛星精密定軌的應用,一方面可以直接應用分析模型輸出的太陽輻射壓加速度時間序列來代替定軌軟件中的光壓模型;另一方面考慮到分析模型可能還存在一定的誤差;可以對其增加若干調節參數,類似于對星載加速度計數據的處理方式,如下
aci=ai+ki×ami
(8)
式中:i=1,2,3(對應于RTN軸),aci為經過標校后的i方向的加速度值;ai為加速度常數偏差參數;ki為i方向的尺度因子參數;ami為觀測值。標校參數在定軌時予以解算。
對GEO-6和IGSO-5衛星,如果直接使用分析型光壓模型,以MGEX軌道作為參考軌道,定軌精度如圖8所示(時間起點為2016年3月3日),GEO-6位置誤差約在2.5m左右,IGSO-5位置誤差約1.5m左右。

圖8 直接使用分析型光壓模型定軌精度Fig.8 Orbit accuracy is determined directly using analytical SRP model
增加標校參數(偏差和尺度因子),GEO-6和IGSO-5定軌精度如圖9所示(時間起點為2016年3月3日),GEO-6位置誤差小于1m,平均誤差約為0.5m。解算出來的尺度因子接近于1,偏差參數約在10-9~10-8m/s2量級。相應IGSO-5的定軌精度比GEO-6的定軌精度均要高,增加標校參數后位置誤差小于0.1m。
上述分析表明:分析型光壓模型可能存在一定的系統性偏差,可以通過建模的方式對系統誤差參數予以解算以提高定軌精度。從定軌結果看,基于分析型光壓模型的IGSO衛星的定軌精度和常規經驗模型定軌精度相當,而GEO衛星的定軌精度略低于常規經驗模型,初步驗證了分析性光壓模型在北斗衛星精密定軌中的應用。

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