一種利用衛星星歷計算機動加速度的方法

崔雙星，楊志濤，劉衛，王榮蘭，向開恒，賀泉，吳中灝

（1.中國星網網絡創新研究院，北京 100029；2.中國科學院國家天文臺，北京 100012；3.中國科學院國家空間科學中心，北京 100190）

1 引言

軌道機動是地球衛星進行軌道控制的主要方式，尤其對于近地軌道衛星，由于受大氣影響較大，軌道衰減明顯，使得軌道機動比較頻繁。相比傳統脈沖式發動機，電推進具有高比沖、高效率和小推力的特點。典型的電推進系統比沖1500～3000s，推力多在毫牛至數百毫牛間，工質多為氙、氪或碘[1-3]。采用高效的小推力發動機實現航天器機動，可以有效降低任務過程中的燃料消耗。隨著商業航天的發展，電推進技術廣泛用于航天器姿態控制、位置保持和軌道機動等方面[4,5]。

關于小推力電推進多集中在控制策略、優化方法等方面的研究[6-9]，鮮有其逆向問題，即由軌道反演電推力的相關研究工作。由軌道反演推力也多集中在一段時間內的累積推力或平均推力的計算或定性分析上。羅冰顯等人針對2022 年2 月美國太空探索技術公司（SpaceX）的星鏈衛星受地磁暴影響墜毀大氣層事件，展開大氣阻力與推進策略相關研究工作，得出軌道推力無法補償大氣阻力可能是導致衛星再入的主要原因[10]。相比傳統軌道反演累積推力或平均推力方法，本文提出一種識別軌道機動起止時刻的精確方法，從而可進一步較精確地計算出機動時間內的機動推力。這對于分析小推力發動機的性能、衛星軌道維持與規避策略設計等面都具有重要的意義。

目前，隨著世界各國對太空探索任務的逐步深入，航天任務日益頻繁，空間目標種類和數量不斷增加，太空環境變得愈來愈復雜，為保障空間目標運行安全，軌道機動的快速檢測尤為重要。同時，為降低地磁暴等大氣層事件對低軌衛星的再次影響，對歷史軌道機動識別與推力的計算，也能為未來衛星的設計與運行提供必要支持。

本文通過分析星鏈衛星的星歷數據，對衛星的軌道機動進行分析，研究機動時刻確定及機動加速度計算的具體方法。

2 星歷數據處理

星鏈是美國太空探索技術公司的一個項目，太空探索技術公司計劃在2019年至2024年間在太空搭建由約1.2萬顆衛星組成的“星鏈”網絡提供互聯網服務，其中1584 顆將部署在地球上空550千米處的近地軌道，并從2020年開始工作。據有關文件顯示，該公司還準備再增加3 萬顆，使衛星總量達到約4.2萬顆[11-13]。

美國Space-Track網站，每天定時發布星鏈衛星的星歷數據，發布時間間隔為8h，每個星歷文件時間跨度為3 天[17]。星歷文件由衛星的位置速度和協方差矩陣組成。由星歷文件中的協方差矩陣，可以得出衛星星歷的沿跡誤差，由星歷文件的位置速度，可以推得衛星的瞬時軌道根數。圖1是2022年11月7日發布的Starlink 1007衛星星歷沿跡誤差和傾角演化圖。

圖1 Starlink 1007星歷沿跡誤差和傾角演化圖Fig. 1 Starlink 1007 along-track error and inclination changes with time

從誤差演化圖可以看出，前半段數據誤差隨時間增長顯著變大，這符合低軌道物體的誤差演化規律，后半段數據的誤差在最大值的基礎上保持不變，應該是誤差取了極大值。由傾角的變化圖也可以明顯看出前半段數據與后半段數據的差異，前半段數據可以明顯看出短周期項的變化，后半段數據只反應了長周期的變化，推斷前半段數據生成采用了高精度的軌道預報模型，后半段數據則是采用了簡化的軌道預報模型。所以前半段數據更為準確，本文僅采用前16小時的數據進行分析。

3 衛星機動的識別

衛星在運動過程中受地球非球形引力、大氣阻力、光壓、日月引力等多種攝動的影響，引起大量的周期攝動，由于在小推力下的機動推力非常小，基本上與大氣阻力攝動屬于同一量級，在瞬時軌道根數下，小推力淹沒在其它攝動力的作用中，很難識別出小推力的變化，因此，很難直接通過衛星星歷識別出衛星機動的起止時間。為了有效識別衛星的機動，就需要剔除衛星其它攝動力的影響，過濾掉其它攝動力引起的短期變化。

由于軌道平均根數不包含軌道攝動的各種短周期項，因此將瞬時根數轉換為平均根數，也可以有效剔除其它攝動力的影響。在瞬時根數轉平均根數過程中需要充分考慮多種攝動因素的影響[14,15]，本文的方法中采用了8 階GGM03C 引力場模型，考慮了日月引力和光壓攝動模型。

在平均根數條件下軌道機動的判別有多種方法，最簡單的辦法是直接通過半長軸的變化進行判斷[16]。判別標準是：

式中，Δa為半長軸的變化量，Δath為半長軸變化的均值，σ為半長軸變化的標準差。

4 機動加速度的計算方法

本文采用兩種計算機動加速度的方法：二分迭代法和能量轉換法。其中二分法迭代從衛星的運動力學模型出發，采用數值法軌道預報推算機動時段內的加速度，是較精確的方法，但由于存在多次的數值法軌道預報，計算速度較慢，不適用于大規模的應用，但可以作為其它方法準確性的驗證。能量轉換法，從衛星能量變換出發，并進行適當近似，推算機動加速度。由于存在近似計算，該方法精度較二分迭代法低，但運行速度較快。

4.1 二分迭代法

采用軌道預報迭代反演的方法可以推算機動時段的機動加速度，軌道預報采用數值法模型，機動力模型采用沿速度方向恒定推力的方式，通常小推力的方向與衛星運動方向基本一致，因此本文簡單地按只有沿速度方向的分量來處理，不過針對推力方向與速度方向差別較大的情況，若能準確掌握二者的夾角，仍然可通過類似的方法分析。由于小推力機動過程中，衛星本身的質量變化很小，可以近似認為質量不變。在此條件下，假定衛星時刻在速度方向上遭受了恒定的推力，機動時段內衛星的機動推力向量可按以下公式計算：

式中，f為衛星的速度矢量，v為衛星的速度標量。

由該機動力產生的機動加速度計算公式為：

式中，a為機動加速度，為機動開始時刻，為機動結束時刻。

將以上推力模型加入數值法軌道預報算法中，通過二分迭代法反演軌道機動推力，計算流程如圖2。具體如下：

圖2 機動加速度反演流程Fig. 2 Flow chart of the maneuver acceleration inversion process

（1）根據長半軸變化法判斷軌道機動的起止時刻，將軌道按照機動時段分為機動時段前軌道和機動時段后軌道；

（2）根據機動時刻前的星歷數據進行軌道確定，獲得機動開始時刻的軌道根數；

（3）根據機動時刻后的星歷數據進行軌道確定，獲得機動結束時刻的軌道根數；

（4）設置機動加速度的初值和邊界值，邊界值可以根據報道的機動推力和質量估算，即最大可能推力與最低質量的比值（小推力一般在幾mN 到幾百mN 之間，衛星質量一般不小于200kg，假定最大推力為500mN，因此可設定初始加速的邊界值為0.0025m/s2），加速度的初值設為0m/s2；

（5）根據機動加速度得到機動力模型，將模型疊加到軌道預報模型，將機動開始時刻的軌道根數預報至機動結束時刻；

（6）對比預報時刻的軌道半長軸與真實根數的半長軸，與設置的半長軸差的閾值進行比較，滿足，則輸出機動加速度，否則調整機動加速度，返回步驟（5）重新計算。

4.2 能量轉換法

當機動時段比較短時，可以忽略由速度變化而引起的大氣、三體等攝動力引起的能量變化，這時，可以認為機動情況與非機動情況衛星動能和勢能的變化，近似等于機動力對衛星所做的功，因此有：

即，

式中，m為衛星的質量，GM為引力常數，s為機動時段衛星運動軌跡的長度，為機動后衛星的勢能，為機動后衛星的動能，為無機動條件下衛星在機動結束時刻的勢能，為無機動條件下衛星在機動結束時刻的動能，v′和v是分別為機動結束時刻無機動條件與有機動條件下衛星的速度，r和分別為機動結束時刻無機動條件與有機動條件下衛星的地心距。

由于機動時段一般比較短，r的變化很小，可以近似為常數，所以有，

s由數值法軌道預報積分得出，當衛星軌道為近圓軌道時，s可以近似為圓弧。這種情況下，s可以由以下公式簡單計算，

式中，r為機動時段衛星掃過的角度。

若機動時段較長，r可能變化較大，這時候可對機動時段進行數值法軌道預報，在預報過程中通過積分器在同步計算得出r，也可以直接由機動時段內的星歷對上述公式進行離散化，在每個離散時段內再將r近似為常數計算得出。

5 實例分析

5.1 機動加速度計算

這里選取2022 年11 月7 日發布的Starlink 1007 目標的星歷，利用以上方法進行分析。Starlink 1007衛星的目標編號為44713，軌道高度550km左右，軌道傾角53.06°。

星歷文件名稱為：

MEME_44713_STARLINk-1007_3102130_Operational_1352064660_UNCLASSIFIED.txt

提取星歷文件中前16小時的星歷數據，轉換為開普勒瞬時根數，再將瞬時根數轉換為平均根數。瞬時根數和平均根數下半長軸的演變規律如圖3和圖4所示。由圖可知，平均根數條件下，半長軸出現明顯的跳變，圖4 右是半長軸跳變處的局部放大圖。

圖3 瞬時根數半長軸演化圖Fig. 3 The evolution diagram of the semi-major axis of the osculation orbit

圖4 平均根數半長軸變化圖Fig. 4 The evolution diagram of the semi-major axis of the mean orbit

通過半長軸變化判據，得出衛星機動時段為：2022/11/07 05:08:42.000 UTC 至 2022/11/07 05:11:42.000UTC，機動持續時長180s。分別使用機動時段前的數據和機動后的數據進行定軌，可以得到衛星機動前（機動起始時刻）與機動后（機動結束時刻）的軌道根數。這里定軌方法采用的是基于最小二乘的批處理算法，定軌結果如表1。

表1 機動前后的軌道參數Table 1 Orbit elements before and after maneuver

采用考慮30 階GGM03C 引力場、日月引力、光壓、大氣、潮汐、相對論效應等攝動的數值法軌道預報模型，將機動前的軌道根數疊加小推力模型，進行軌道預報，將軌道預報結果與發布結果的長半軸進行比較，進一步采用二分法迭代，迭代31次半長軸誤差可小于1m，計算可得機動加速度為：2.62381×10-4m/s2。

采用能量轉換法，利用機動前的定軌根數進行軌道預報，得到無機動條件下的軌道根數，在此基礎上，獲得機動終止時刻的動能和勢能，與發布的機動終止時刻的動能和勢能進行比較，得到小推力作用下的能量變化，在此基礎上根據公式（5）計算得到的機動加速度為：2.59545×10-4m/s2，與二分迭代法得到的結果非常接近，偏差約為1%。

采用上述方法對其它星歷進行同樣的分析，選取NORAD編號44714、44715、44716、44717、44718 和編號50844、50819、50832、50827 和50836的10個目標，對應星鏈衛星Starlink 1008-1012、Starlink 3326-3330，提取11月7日的星歷數據，仍然選擇前16個小時的數據進行分析。對應的星歷文件為：

MEME_44714_STARLINk-1008_3102159_Operational_1352066400_UNCLASSIFIED.txt

MEME_44715_STARLINk-1009_3102142_Operational_1352065380_UNCLASSIFIED.txt

MEME_44716_STARLINk-1010_3102108_Operational_1352063340_UNCLASSIFIED.txt

MEME_44717_STARLINk-1011_3102124_Operational_1352064300_UNCLASSIFIED.txt

MEME_44718_STARLINk-1012_3102156_Operational_1352066220_UNCLASSIFIED.txt

MEME_50844_STARLINk-3326_3110423_Operational_1352089440_UNCLASSIFIED.txt

MEME_50819_STARLINk-3327_3110436_Operational_1352090220_UNCLASSIFIED.txt

MEME_50832_STARLINk-3328_3110425_Operational_1352089560_UNCLASSIFIED.txt

MEME_50827_STARLINk-3329_3110436_Operational_1352090220_UNCLASSIFIED.txt

MEME_50836_STARLINk-3330_3110419_Operational_1352089200_UNCLASSIFIED.txt

采用長半軸法判定機動時段，在此基礎上利用能量轉換法進行分析，計算得到的機動信息如表2。

表2 機動信息表Table 2 Maneuver information

從計算結果可知，這些目標的機動時長為180～240s左右，當機動時長為180s時，機動加速度在0.00025m/s2左右，當機動時長為240s時，機動加速度為0.00019m/s2左右。每次機動的推力和時長成反比，推力大時，機動時間短。

5.2 計算結果驗證

根據上節計算的機動加速度構建機動模型，將機動模型加入預報模型，使用機動前的軌道進行預報，將預報星歷與發布星歷進行對比，這樣就可以得出機動加速度的計算與星歷的符合情況，對計算結果進行驗證。

圖5 是衛星44713在無機動和有機動的條件下，衛星預報星歷與發布星歷之間的半長軸的差值，圖6是這兩種條件在J2000坐標系下xyz三個方向的差值變化圖。由圖7可知，加入機動模型后，衛星預報與原始星歷之間的誤差大大降低，半長軸誤差由80多米下降到3m左右，在J2000坐標系下，三個軸的分量的誤差，由公里級下降到米級，可見加入機動推力能很好地符合原始的軌道變化。圖5中在機動時刻出現波動，分析原因可能為識別的機動起止時刻與真實機動起止時刻之間的誤差以及插值方法所致，因為星歷文件是間隔60s的離散的時刻點，因此，實際的軌道機動開始時間與終止時間可能與計算值存在0～60s范圍內的誤差。

圖5 有無軌道機動預報與原始星歷半長軸差值變化圖(左：無軌道機動，右：有軌道機動）Fig. 5 The difference of the semi-major axis of the prediction orbit without maneuver and with maneuver

圖6 在有無軌道機動下預報星歷與原始星歷差值和方向的投影Fig. 6 The difference in x,y,z-axis direction

圖7 有無軌道機動預報與原始星歷半長軸差值變化圖Fig. 7 The difference of the semi-major axis of the prediction orbit without maneuver and with maneuver

下面給出在無機動和有機動的條件下，部分衛星預報星歷與發布星歷之間的半長軸的差值變化圖和在xyz 方向的差值變化圖。圖7、圖8、圖9、圖10分別為NORAD編號44714、44715和編號50819、50827衛星的半長軸差值變化圖和在xyz方向的差值變化圖。由圖可知通過能量轉換法計算的推力可以與原始星歷很好地吻合。

圖8 在有無軌道機動下預報星歷與原始星歷差值在x方向的投影Fig. 8 The difference in x-axis direction

圖9 在有無軌道機動下預報星歷與原始星歷差值在y方向的投影Fig. 9 The difference in y-axis direction

圖10 在有無軌道機動下預報星歷與原始星歷差值在z方向的投影Fig. 10 The difference in z-axis direction

6 結論

本文利用公開發布的星鏈星歷數據，對基于星歷進行衛星軌道機動識別及機動加速度計算方法進行分析，通過對軌道半長軸變化，確定軌道機動的時段的識別方法，并采用數值積分反演的方式和衛星能量轉化的方法，確定軌道機動時段的推力。通過該方法仿真分析了多顆星鏈衛星在2022年11月7日的機動，并計算了機動推力。仿真結果表明，通過本方法解算得到的機動加速度疊加到軌道預報模型后預報產生的星歷與發布的星歷基本一致，證明了軌道機動時段確定和機動加速度解算的正確性。這些結果可為我們認識星鏈星座的推進策略提供幫助。