大規模星座近圓無奇異同步快速軌道預報方法

聶 濤，潘政雨，李金峰，徐 瑞，秦 凡，崔平遠

(1. 北京理工大學宇航學院，北京 100081；2. 深空自主導航與控制工信部重點實驗室，北京 100081)

0 引 言

20世紀90年代“銥星”(Iridium)等提供通信和網絡服務的低軌星座工程問世以來,低軌衛星星座便廣受關注,這主要是因為軌道高度低導致傳輸時延短、路徑損耗小,但是大部分早期項目都以失敗而告終,尤其是銥星系統的破產為低軌星座的發展蒙上了陰影。近些年,隨著集成化和自動化技術的發展,制造衛星成本以及發射成本逐漸降低,市場需求量也不斷擴大,以美國SpaceX公司的42 000顆“星鏈”(Starlink)巨型衛星互聯網星座為代表,各國紛紛出臺數百顆、數千顆星的低軌巨型星座計劃,低軌巨型星座再次掀起前所未有的熱潮[1]。

低軌巨型星座是指在高度400～2 000 km近地軌道上分布,由幾百顆及以上(可高達幾萬顆)衛星構成,具有廣覆蓋、高通量、低延時等特征的大型衛星星座,根據星座功能不同,可以分為巨型互聯網星座、巨型導航星座、巨型遙感星座等。目前,國外最具代表性的低軌巨型星座主要包括美國SpaceX公司的Starlink星座和英國OneWeb公司的OneWeb星座[2]。OneWeb星座于2014年提出,初期部署648顆均勻分布在18個軌道高度為1 200 km、軌道傾角為87.9°的近極圓軌道上,旨在建立覆蓋全球的高速電信網絡,截止到2022年12月,OneWeb已經發射部署502顆衛星。與OneWeb星座相比,Starlink星座規模更加龐大,發展更為迅速,目前已經發射部署3 000多顆衛星,主要分布在多個軌道高度為550 km左右的傾斜軌道上。

中國在2020年將衛星互聯網計劃納入“新基建”范疇,標志著衛星互聯網項目成為國家戰略性工程。同年9月,中國向國際電信聯盟ITU提交了由12 992顆衛星組成、以“GW”為代號的低軌巨型星座申請網絡資料,由此,開啟了中國的低軌巨型星座建設新的篇章。

對照歷史數據并參照巨型星座發展計劃,預計未來10年內發射的小衛星總數將超過3萬顆,并且絕大多數將采用300 km至2 000 km的低軌軌道,巨型互聯網星座計劃和發射數量呈爆炸式增長,逼近近地空間的承載極限,極易與空間物體尤其是空間碎片發生碰撞[3-4]。據ESA估計,目前在軌超過10 cm的碎片約有34 000個,大小在1 cm至10 cm之間的碎片個數約為90萬,而在1 mm至1 cm之間的碎片數量則超過1.2億個。

碰撞風險評估的前提是對星座衛星軌道進行軌道預報,目前國內外對于單顆星的軌道預報方法比較多,而如何對多顆星進行同步預報的研究比較少。不考慮攝動影響時,航天器將沿著不變的圓錐曲線運動,而在攝動影響下,衛星的軌跡變得復雜,軌道要素會隨著時間發生變化,常見的攝動力主要包括中心引力場的非球形引力攝動、月太三體攝動、太陽光壓攝動、大氣阻力攝動等,與中心引力相比,攝動力屬于小量,故而描述軌道運動的攝動方程對應于小參數方程,對于攝動方程的求解主要包括數值積分法、解析法和半解析法三類,數值方法可以對攝動力進行精確建模,依賴于數值積分方法對軌道動力學微分方程進行求解,為了獲得高精度結果,需要采用比較小的積分步長,而解析方法通過泰勒級數以及平均方法等近似處理可以獲得軌道動力學方程的解析解,計算量比較小并且可以獲得軌道運動的一般規律,而半解析方法綜合了解析法與數值法的特點,在計算精度與計算效率之間取得折中。

對于低軌衛星來講,J2攝動占主導,Kozai[5]在1959年提出平均軌道根數方法對J2攝動進行求解,對于J2攝動來講,它主要會引起升交點赤經漂移以及近地點旋轉,而軌道半長軸、偏心率以及軌道傾角長期保持不變,Kuiack等[6]注意到這一點,提出了一種考慮J2攝動的非線性相對運動計算方法,國內也有一些學者針對攝動下的軌道動力學問題進行了研究,徐光延等[7]提出了一種任意階帶諧非球形引力攝動的相對運動建模方法,蒼中亞等[8]提出根據衛星觀測的極光能量注入數據改進地磁擾動期間引起的軌道預報誤差方法,聶濤等[9]建立了考慮前4階非球形引力攝動的平均軌道動力學模型,同時針對三體攝動給出了半解析計算方法[10-11],孔繁澤等[12]利用圖形處理器并行計算方法對多星軌道進行加速計算,克服傳統軌道預報速度的瓶頸,文援蘭等[13]研究了星座衛星的軌道確定問題,此外,也有一些學者針對星座設計問題[14-17]、相對運動控制問題[18]以及碰撞風險評估[19]等問題進行了研究。

本文對大規模星座軌道預報問題開展研究,考慮大規模星座衛星數量多、軌道攝動復雜等問題,提出了一種基于平均速率矩陣的多星同步軌道快速預報方法。首先,建立了J2攝動下的平均軌道動力學模型,并針對小偏心率存在奇異性問題,提出了一種近圓無奇異解析同步快速軌道預報模型,與文獻[6]方法同時含有真近點角和平近點角不同,本文給出了僅含平近點角的平瞬轉換關系,簡化了軌道預報過程。在此基礎上,建立了基于平均速率矩陣的多星軌道同步預報算法,能夠實現對多顆衛星軌道進行同步預報。最后,通過仿真驗證了算法的有效性。

1 J2平均動力學模型

對于低軌衛星來講,它受到的主要攝動力為J2攝動,本節將介紹J2攝動的平均軌道動力學理論,包括平均軌道動力學方程以及平均軌道根數與瞬時軌道根數之間的轉換。考慮J2攝動的哈密頓函數可以表示為[5]

(1)

式中:μ是地球引力常數;Re是地球赤道半徑,r是衛星軌道半徑;J2是二階帶諧項變量;f是真近點角;g為近地點幅角;H,G,L為Delaunay軌道要素中的共軛動量變量。

根據Brouwer理論,先后消除哈密頓函數中的平近點角和近地點幅角,可以確定出平均哈密頓函數為[5]

(2)

(3)

平均軌道根數與瞬時軌道根數之間的轉換關系為

(4)

式中[6]:

3sin2i(ecosf+1)3cos(2f+2ω)}

4(1-3cos2i)(η3-(ecosf+1)3)-

12sin2icos(2f+2ω)(η2-(ecosf+1)3)-

2ω))+3cos(2f+2ω)]-

e[3sin(f+2ω)+sin(3f+2ω)-6sinf]+

ecosf-η2+1)-sin(3f+2ω)(3e2cos2f+9ecosf+

η2+6)]+2sinf(3cos2i-1)(e2cos2f+3ecosf+

η2+2)+3e(5cos2i+3)(esinf+f-M)+

e(3-5cos2i)[e(3sin(f+2ω)+sin(3f+2ω))+

2)cos6i-5(31e2+38)cos4i+(43e2+

[sin2isin(3f+2ω)(3e2cos2f+9ecosf+η2+6)+

2sinf(3cos2i-1)(e2cos2f+3ecosf+η2+2)-

3sin2isin(f+2ω)((ecosf+1)2+

(5)

這些代表瞬時軌道根數與平均軌道根數之間轉換的周期修正項。

2 近圓無奇異解析軌道預報方法

根據式(5)可知,當e=0以及i=63.43°或116.57°,平瞬軌道轉換會出現奇異。考慮到星座衛星軌道大多采用小偏心率的近圓軌道,本節將重點針對小偏心率存在奇異的問題,提出一種小偏心率無奇異的解析軌道預報方法。此外,現有的平瞬軌道轉換算法同時含有真近點角和平近點角,本節基于Fourier-Bessel級數理論,推導出一種僅含平近點角的平瞬軌道轉換算法。

(6)

而其他軌道參數的平均速率為常數,為此,

(7)

為了獲得瞬時軌道根數,需要進行平瞬轉換,而根據平瞬轉換關系式(4)和式(5)可知,平瞬轉換中含有真近點角,而平均軌道動力學計算得到的是平近點角,基于Fourier-Bessel級數理論[20]

(8)

式中:Jk(ke)為第k階Bessel函數

(9)

將式(8)代入式(5)中并只保留到e的一階項,可以確定出

e(3cos2i-1)cosM}

4[-3e×cos(2M+2ω)+3cos(M+2ω)+

2ω)]+3cos(2M+2ω)}

18sinM]+3sin(2M+2ω)}

5cos2i-3)+sin(3M+2ω)(sin2i(3cos2M+

13)-5cos2i+3)+sinM(3cos2i(cos2M+20)+

3sin(M+2ω))+(9cos2i+3)sinM]+

3(3-5cos2i)sin(2ω+2M)+9sin2isin2M×

sin(M-2ω)-40e(3cos2i+2)sin2ω+(5cos2i+

3)[3e2sin(5M+2ω)-3(11e2+4)sin(M+

2ω)+(41e2+28)sin(3M+2ω)+18esin

8e2+6)+e2cos2M+8e2+36ecos2icosM+6]+

(10)

觀察平瞬轉換關系式(5)和式(10)可以注意到表達式中偏心率出現在分母中,對于小偏心率近圓軌道來講,它會出現奇異性,為了消除小偏心率奇異性,引入一組新的變量:

a,i,Ω,u,p,q

(11)

(12)

e(9-17cos2i)sin(3M+2ω)+6e(33cos2i+

19)sinM+(6-30cos2i)sin(2M+2ω)]

(13)

而對于變量p,一階平瞬軌道轉換近似為

(14)

將式(5)代入式(14)可以確定出式(15),同理,可以確定出式(16):

41sin(3M+2ω)+3sin(5M+2ω)+3sin(M-

2ω))+2(3cos2i+1)sinM(cos2M+8)]+

cos2i+3)cos2M+15cos2i+13)+9(3cos2i+

1)(5cos2i+3)]-8sin2isinMsin2ω×

[48(5cos2i+3)cos2M+255cos2i+157]]+

(15)

41sin(3M+2ω)+3sin(5M+2ω)+3sin(M-

2ω)]+2(3cos2i+1)sinM(cos(2M)+8)}-

cosM-6(5cos2i+3)cos(3M))-6(5cos2i+

3)sin3Mcos2ω)-8(3cos2i+1)sinM(-5×

sin2icos2ω+15cos2i+9)]-16sin2isin

(16)

從式(13)、式(15)和式(16)可以看到,表達式不存在奇異問題,偏心率不會存在分母中,根據式(15)和式(16)可以確定瞬時軌道偏心率為

(17)

而瞬時平近點角為

(18)

進而可以確定出瞬時近地點幅角為

ω=u-M

(19)

在這種情況下,已知初值軌道根數,利用式(7)可以計算出任意時刻的平均軌道根數狀態,再利用式(4)和式(10)計算出瞬時軌道半長軸、瞬時軌道傾角以及瞬時升交點赤經,最后利用式(13)及式(15)～(19)算出瞬時軌道偏心率、瞬時平近點角以及瞬時近地點幅角。

3 基于平均速率矩陣的多星軌道同步預報方法

前文給出了一種解析的計算單顆衛星軌道狀態的方法,在此基礎上,下文將給出一種多星軌道同步預報方法。

(20)

對于n顆衛星,它們的平均軌道根數為

Xt=X0+Vt

(21)

式中:初始平均軌道根數狀態矩陣為

(22)

速率矩陣為

(23)

在確定所有衛星的平均軌道根數之后,再利用平瞬軌道根數轉換關系即可計算出瞬時軌道根數。

對于多星同步軌道預報算法的計算流程為:

1) 給定所有衛星的初始狀態,利用式(4)確定出衛星的初始平均軌道根數;

2) 基于式(22)和式(23)對衛星初始狀態矩陣X0以及V0進行初始化;

3) 給定任意時刻t,利用式(21)計算衛星的平均狀態矩陣;

4) 利用式(4)和式(10)計算出瞬時軌道半長軸、瞬時軌道傾角以及瞬時升交點赤經;

5) 利用式(13)及式(15)～(19)算出瞬時軌道偏心率、瞬時平近點角以及瞬時近地點幅角。

4 仿真校驗

針對本文提出的大規模星座軌道同步高效預報方法,以美國SpaceX公司的Starlink星座為仿真對象,首先根據公開的雙行根數(Two-line Element,TLE)對Starlink在軌數據進行分析,確定Starlink星座運行軌道數據,之后運用本文提出的軌道預報方法對星座軌道進行預報,將結果與高精度預報模型進行對比,進而驗證所提算法的有效性。

4.1 Starlink星座在軌數據分析

截至2022年11月,SpaceX公司已累計發射65批次共計3 558顆衛星,進入2022年下半年,衛星部署進程明顯加快。根據衛星數據網站CelesTrak上Starlink星座的TLE數據,可以確定出衛星的軌道根數分布,其中軌道高度以及軌道傾角分布圖像如圖1和圖2所示,Starlink在軌衛星數量為3 251顆,這3 251顆衛星中可以分為正在離軌、正在部署、在運行軌道以及偏離軌道等幾種情況。由軌道傾角分布圖像可以看出軌道傾角分布比較集中,主要分布在4個區間[ 52.99°, 53.10°],[ 69.99°, 70.01°],[ 97.51°, 97.70°],[ 53.20°, 53.23°],在這4個不同區間的衛星數量分別為1 491顆、51顆、187顆以及1 552顆。

圖1 在軌衛星軌道高度分布圖像

圖2 在軌衛星軌道傾角分布圖像

進一步對第一組軌道傾角在[52.99°, 53.08°]的1 491顆衛星分布進行分析,得到的軌道高度分布圖如圖3所示。由圖像可以分析衛星的軌道高度主要集中在[547.1,547.4] km,在此區間內的衛星數量為1 445顆,區間外的衛星數量為46顆,同時發現區間外的衛星偏心率以及近地點幅角也不集中,因此,我們認為軌道高度在[547.1,547.4] km的1 445顆衛星為正常運行衛星,而其他衛星可能為正在離軌的廢棄衛星或正在部署的衛星。

圖3 軌道傾角處于[52.99°, 53.08°]衛星的軌道高度分布

最后,對1 445顆衛星按升交點赤經大小對衛星進行排序,得到的仿真圖如圖4所示。由圖可以看到這1 445顆衛星的軌道高度、軌道偏心率、軌道傾角近似相等,升交點赤經分布均勻,不同軌道面的赤經差為5度,得到的衛星數量統計分布圖如圖5所示,每個軌道面的衛星數量為20顆,這與Starlink公開的構型數據相一致,結果驗證了統計分析方法的正確性。

圖4 軌道高度處于[547.1, 547.4]km衛星的升交點赤經分布

圖5 各升交點赤經衛星數量統計

采用同樣的數據統計分析方法對其他殼層(具有相同軌道半長軸、偏心率以及軌道傾角的衛星構成的空間)的數據進行分析,得到的結果如表1所示,結合Starlink公開的發射部署計劃,可以確定出各殼層的部署完成情況,最早開始部署的殼層1(軌道傾角為53°、軌道高度為547 km)已完成部署的91.22%。

表1 星鏈在軌數據分析結果

4.2 Starlink星座演化分析

根據Starlink在軌數據,利用本文提出的算法可以對Starlink衛星在軌演化進行分析,得到的衛星在三維空間以及二維地圖上的分布如圖6和圖7所示,從分布圖像可以看到衛星主要分布在南北緯60度范圍內。本文提出的方法是一種基于平均軌道動力學理論的解析計算方法,同時基于速率矩陣可以同步計算任意時刻的星座的所有衛星軌道狀態,這種方法對在軌的3 251顆衛星計算一次所需的時間為0.005 s(采用Intel(R) Core(TM) i7-10700 CPU @ 2.90 GHz的CPU,基于Matlab R2019b進行仿真)。在相同的仿真環境下,進行對比驗證:如果采用相同的軌道動力學模型但通過循環3 251次計算星座衛星軌道,計算所需的時間為0.046 s,這就意味著本文提出的速率矩陣同步計算策略可以將計算速度提高一個數量級;如果考慮前15階非球形引力攝動以及大氣阻力、月/太三體攝動和太陽光壓等攝動影響,其中衛星的阻力系數和質量以及面積分別設置為2.4,6 kg,0.15 m2[9],利用非線性無簡化動力學模型,并采用4階龍格庫塔積分方法進行高精度軌道遞推,高精度數值方法對3 251顆衛星進行軌道遞推一周,所需的時間為幾十小時。在星座中任選一顆衛星,獲得與高精度數值方法的誤差對比結果如圖8所示,可以注意到一周內,解析方法的積分誤差小于2.7 km。

圖6 星鏈三維空間分布

圖7 星鏈二維空間分布

5 結 論

本文針對大規模星座軌道預報方法開展研究,考慮J2攝動影響,建立了基于平均軌道根數的近圓無奇異解析軌道預報模型,并提出了基于平均速率矩陣的多星軌道同步預報方法,最后,以Starlink星座為仿真對象,基于公開的TLE數據對星座軌道進行了分析,結果表明目前Starlink星座主要分布在4個不同軌道傾角的殼層,其中軌道傾角為53°、軌道高度為547 km的殼層已經完成計劃部署的91.22%,并利用本文提出的同步軌道預報方法對Starlink星座衛星軌道進行預報,從計算速度角度看,相比于單次預報方法,本文提出的同步預報方法可將計算速度提高一個數量級,而從計算精度的角度看,解析預報方法對軌道預報7天積分誤差小于2.7 km。