基于SGP4模型的空間站跟蹤算法

王 霄,李偉超,楊旭海,3,欽偉瑾,王 偉,李志剛

(1. 中國科學院國家授時中心,陜西 西安 710600; 2. 中國科學院精密導航定位與定時技術重點實驗室,陜西 西安 710600; 3. 中國科學院大學天文與空間科學學院,北京 100049)

空間站又稱為太空站或軌道站,是運行在外太空的軌道艙,屬于航天器的一種;與其他人造衛星、空間探測器等航天器的區別在于,空間站更適合人類長時間居住,可以作為航天員在太空停留和開展工作的場所。由于空間站能夠提供地面試驗設施所不能提供的微重力、宇宙空間環境等條件,因此主要被用于開展各種科學研究。研究領域涵蓋物理學、天文學、地理學、生物醫學等多個領域[1]。早在20世紀90年代,我國載人航天工程就確立了“三步走” 發展戰略。20余年來,多次將多名航天員送入太空,全面突破和掌握了載人天地往返、空間出艙及空間交會對接等技術,由此中國載人航天工程已正式步入空間站時代[2]。

2021年4月29日,中國空間站天和核心艙的發射成功拉開了我國空間站建設的帷幕[3],為后續關鍵技術驗證和空間站組裝建設順利實施奠定了基礎。空間站在軌運行期間,實現地面站對空間站的跟蹤觀測是獲取科學數據的前提。2022年2月12日,中國載人航天工程辦公室宣布,為便于各方開展空間目標避碰等軌道分析工作,即日起每天在其官網發布中國空間站兩行軌道根數(two line orbital elements,TLE),供各航天實體參考使用。

中國空間站運行軌道高度為400～450 km,軌道傾角約為42°,運行周期約為1.5 h,運行速度約為7.68 km/s。由于空間站視向速度變化快、過境時間短[4],利用地面站對其進行觀測時,不僅要求地面跟蹤設備的轉動速度跟上空間站的運行速度,還需要地面觀測設備能快速捕獲目標,進行準確跟蹤,盡可能獲得較多的觀測數據。而根據空間站的TLE準確計算地面觀測站的跟蹤角度是實現這些要求的關鍵。

本文從空間目標運行特性的一般規律出發,基于現有的簡化常規攝動(simplified general perturbations,SGP4)模型,研究地面站地平式設備觀測空間站的跟蹤角度的計算算法,以期為地面觀測設備的轉動速度設計和空間站跟蹤觀測提供支持,為后續觀測類似空間站等低軌目標的地面站布局提供參考。

1 航天器軌道參數

如果把航天器和地球都看作質點,那么地心應為航天器運行軌跡的一個焦點,其運動軌跡可用開普勒軌道參數進行描述,即軌道平面半長軸a、軌道偏心率e、近地點角距ω、軌道平面傾角i、升交點赤經Ω及該歷元的真近點角f,這6個參數可以確定它在空間的運動位置,如圖1所示。美國的Celestrak發明創立的TLE,就是以上述6個軌道參數之間的數學關系描述太空航天器運動的位置和速度等信息的[5-7]。衛星、飛行器等空間目標一旦進入太空,即被列入北美防空聯合司令部(North American Aerospace Defense Command,NORAD)衛星星歷編號目錄并被終生跟蹤。同時,NORAD也會定期生成這些空間飛行器的TLE并公開發布,其軌道精度在百米量級。被公開發布的TLE具有較強的通用性,因此常被作為航天器地面跟蹤系統、軌道計算軟件的輸入條件。

圖1 空間飛行器軌道參數

理想情況下,上述6個軌道參數可以完全確定空間站的位置及運動情況。而實際上它在運動過程中受擾動力的影響,如其他天體的引力、地球形狀攝動、大氣阻力、潮汐攝動、太陽輻射壓力等,使得空間站在其軌道面上的運動十分復雜[8]。因此TLE必須與配套模型一起使用才能達到較好的預報精度。當前發布的TLE使用SGP4模型或簡化深空攝動(simplified deep space perturbations,SDP4)模型,考慮了大氣阻力、地球非球形攝動等影響項[9]。其中,SGP4模型適用于運行周期小于225 min的近地空間目標;SDP4模型適用于運行周期大于225 min的空間目標。由于空間站運行軌道較低,因此選用SGP4模型對其位置進行預報和跟蹤。

2 空間站跟蹤算法

地面站在對空間站目標進行觀測時,需提前根據空間站位置計算地面設備觀測的方位角和高度角,以便地面觀測設備不斷調整角度來跟蹤目標。地面站跟蹤觀測空間站的角度由地面站與空間站之間的相對位置關系決定,空間站在以地面跟蹤設備的旋轉中心為原點的站心地平坐標系中的位置可以直觀地反映這一相對關系,即計算地面站觀測空間站的角度,關鍵在于獲得空間站在站心地平坐標系中的位置坐標。

2.1 地心天球坐標系下空間站位置計算

TLE是一段時間內的平均軌道參數,它采用特定方法去除了周期擾動項,本文參考Spacetrack Report No.3[10]中的SGP4模型,由于時間久遠,對該模型作了部分修正和完善。Spacetrack Report No.3中,為簡化計算,對地球平均半徑作了歸一化處理。本文為方便理解,不再作歸一化處理,部分參數的取值與Spacetrack Report No.3中的取值存在差異,具體見表1。原本的參數在取值時,距離以km為單位,時間以min為單位;現距離仍以km為單位,時間則以s為單位。

表1 常數項含義及取值

采用SGP4模型計算空間站在地心天球坐標系下的位置,主要步驟為:①計算相關參數,包括空間站初始平均速度、軌道半長軸、大氣密度參數等;②由空間站TLE恢復出軌道平根數;③由空間站軌道平根數依次計算長期項、長周期項、短周期項的值;④計算空間站目標的位置。本文僅列出與Spacetrack Report No.3中SGP4模型計算中存在差異的公式,而相同的部分不再贅述。

根據空間站TLE中的空間站每天繞地球的圈數ne計算空間站運動的角速度,公式為

(1)

大氣密度參數s和(q0-s)4的取值與觀測目標的近地點高度有關。近地點高度的計算公式為

PERIGEE=a″0(1-e0)-ae

(2)

式中,e0為軌道偏心率。對于空間站而言,近地點高度大于156 km。s和(q0-s)4的值見表1。

由于未對地球平均半徑作歸一化處理,參數C1—C5及D2—D4的計算公式分別調整為

(3)

C1=B*×C2

(4)

(5)

(6)

(7)

(8)

(9)

(10)

式(3)—式(10)中,k3的含義及取值見表1;B*為TLE中的Bstar值;ω0為一段時間內的平均近地點角距,由空間站TLE給出。

在計算軌道參數短周期變化量時,為避免重復計算造成的精度損失,不再直接求u的值,而是利用sinu和cosu對短周期項進行求解,即

(11)

(12)

(13)

(14)

在最終計算所用的參數值中,由于前面未對u直接求解,因此無法獲得uk的值,但可以間接獲得sinuk和cosuk的值,用于最終的位置計算,即

sinuk=sin(u+Δu)=sinucosΔu+cosusinΔu

(15)

cosuk=cos(u+Δu)=cosucosΔu-sinusinΔu

(16)

最終獲得空間站在地心天球坐標系中的位置坐標為

(17)

2.2 地面站觀測空間站角度計算

計算地面站觀測空間站的角度,關鍵在于獲得空間站在以地面觀測站為坐標原點的地平坐標系中的位置坐標,使用SGP4模型已經獲得空間站在地心天球坐標系下的位置坐標,對該坐標系進行旋轉、平移變換,最終可獲得空間站在站心地平坐標系中的位置坐標。

為實現地心天球坐標系與站心地平坐標系之間的轉換,引入一個中間坐標系——地心地球坐標系,它與地球固連在一起,隨地球的運動而運動,因而又稱為地固坐標系。地固坐標系的坐標軸指向隨地球的運動發生變化。由于地球運動十分復雜,具體表現為地球自轉的方向和速度都會隨時改變[11]。地球自轉軸的方向既在空間上變化,又在地球本體內變化,空間變化表現為歲差和章動,在地球本體內的變化表現為極移;地球自轉速度變化導致世界時的不均勻性。考慮到歲差和章動影響,將地心天球坐標系分為瞬時真天球坐標系和瞬時平天球坐標系,瞬時平天球坐標系即協議天球坐標系;考慮到極移影響,將地球坐標系分為瞬時地球坐標系和協議地球坐標系[12]。TLE使用的天球坐標系是以真赤道面作為基本面、X軸由地心指向平春分點的坐標系。該坐標系轉換到協議地球坐標系,需根據當前時刻地球自轉情況及地球自轉軸的極移情況計算旋轉矩陣。由于極移使得X軸方向變化不超過0.4″,Y軸方向變化不超過0.6″,均是一個小量[13],因此為簡化計算,本文暫不考慮極移對坐標系轉換的影響,僅考慮地球自轉的影響。地球自轉角[14](Earth rotation angle,ERA)與世界時(UT1)存在非線性關系,即

ERA=280.460 618 37+360.985 647 366 29×

(JD-2 451 545)+3.879 33×10-4×

(18)

式中,JD表示當前UT1時刻對應的儒略日;T0=(JD-2 451 545)/36 525,表示從2000年1月1日12時(UT1)(這一時刻對應的JD=2451 545.0)起算的儒略世紀數。任意UTC時刻對應的UT1時刻計算方式為

UT1=UTC+(UT1-UTC)

(19)

由于UT1與UTC時刻的差不超過0.9″,因此也可直接用UTC時刻進行計算。

空間站在地心地球坐標系中的坐標為

(20)

如需考慮極移及UT1與UTC時刻的差對坐標系轉換的影響,可在IERS官網下載地球定向參數(Earth orientation parameters,EOP),用于極移和世界時修正。

為獲得空間站在站心地平坐標系中的位置,需提前已知觀測站在地心地球坐標系中的位置。一般利用三維坐標或經緯度、高程表示測站位置,測站三維坐標與經緯度、高程之間可相互轉換。測站經緯度、高程轉換為測站三維坐標的公式分別為

f=1/298.257 224

(21)

(22)

S=(1-2f+f2)×C

(23)

(24)

空間站在以該測站為坐標原點的地平坐標系中的位置坐標為

(25)

根據空間站在站心地平坐標系中的坐標,可以得到地面站觀測空間站的方位角和高度角分別如下

(26)

(27)

式(26)計算的方位角范圍為[-90°,270°]。如果觀測設備方位軸轉動范圍為0°～360°,只需對式(26)計算出的方位角小于0°的角度加上360°即可。對于高度角而言,地面設備在觀測時的高度角大于0°時,目標才有可能被地面設備觀測到,即空間站可觀測條件是空間站在地平坐標系中Z方向的坐標值大于0。

TLE提供的空間站軌道精度約在百米量級,使用SGP4模型對空間站位置進行預報時,隨著預報時長的增加,空間站位置誤差不斷變大。假設空間站在站心地平坐標系中各個坐標方向的位置誤差分別為Δx、Δy、Δz,根據式(26)、式(27),空間站位置誤差引起的觀測角度誤差為

(28)

(29)

可以看出,地面站觀測空間站的角度誤差不僅與空間站位置誤差(即軌道誤差)有關,還與其所處位置有關。地面站觀測高度角升高時,地面觀測站與空間站距離逐漸變小,空間站在Z方向軌道誤差對高度角誤差的影響減小。方位角接近90°或270°時,方位角計算誤差主要受空間站Y方向軌道誤差影響,高度角計算誤差受空間站X方向軌道誤差影響;方位角接近0°或180°時,方位角計算誤差則主要受空間站X方向軌道誤差影響,高度角計算誤差受空間站Y方向軌道誤差影響。

3 空間站跟蹤角度計算結果及分析

為驗證本文算法的正確性,選取空間站在2021年11月8日的TLE,分別采用本文算法和STK軟件計算西安觀測空間站一周的方位角和高度角,并對兩種方法得到的結果作差。對于方位角差值,兩種方法得到的方位角直接相減得到的僅是在水平面上方位角的差值。而在實際觀測中,不同的高度角會導致方位角所對應的弧長不同。因此對于不同的觀測高度角,方位角誤差導致的指向角誤差也不相同。在計算方位角誤差時,不僅要考慮上述兩者的差值,還要考慮該觀測歷元的高度角。圖2為考慮高度角后的方位角差值,圖3為兩種方法計算的高度角差值。可以看出,方位角和高度角差值均在0.001°以內,證明了該算法的正確性。

圖2 兩種方法計算的方位角差值

圖3 兩種方法計算的高度角差值

此外,為驗證空間站TLE的有效預報時長,利用上述算法,根據2021年11月9—14日(年積日313—318)的空間站TLE,計算西安觀測空間站的方位角和高度角,并與上述根據2021年11月8日的空間站TLE外推計算的一周的結果進行對比。如圖4、圖5所示,方位角預報誤差和高度角預報誤差隨預報時長的增加逐漸增大,在年積日313—315內(即11月9—11日),方位角和高度角的預報效果較好,誤差在2°以內,11日后預報效果較差,方位角和高度角預報差值最大超過8°。因此在地面站觀測空間站時,地面站所使用的TLE的更新周期最長不宜超過3 d。

圖4 方位角預報差值

圖5 高度角預報差值

地面站在觀測空間站時,隨著觀測高度角的升高,空間站越來越接近近地點,地面觀測站與空間站之間的距離越來越小,根據開普勒第二定律,空間站的運動速度越來越大,在通過近地點時速度達到最大。圖6和圖7為地面站觀測空間站時,在部分觀測弧段內的方位角速度和高度角速度隨觀測高度角的變化。在觀測高度角超過80°時,方位角速度急劇增加,最大超過16°/s。由于空間站在高度角超過80°時運行速度過大,使得地面觀測設備的轉動速度無法滿足空間站的運動速度要求,因此一般情況下僅觀測高度角低于80°的弧段。

圖6 方位角速度隨高度角變化

圖7 高度角速度隨高度角變化

4 結 語

根據空間站軌道參數研究地面站的觀測角度是實現空間站跟蹤觀測的前提。本文以空間站TLE為輸入條件,根據低軌衛星運行特點,利用修正后的SGP4模型計算得到空間站在地心天球坐標系下的位置,經過坐標系旋轉等變換后,得到空間站在站心地平坐標系中的位置,并最終計算得到地面站觀測空間站目標的跟蹤角度;將以上算法的計算結果與STK軟件的計算結果進行對比發現,兩者差值小于0.001°,證明了本文算法的正確性,分析了利用TLE對空間站進行軌道預報時其有效時長,TLE更新時長不宜超過3 d。此外,本文提出的跟蹤算法不僅適用于空間站目標,還適用于其他運行周期小于225 min的低軌空間目標。