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

（1.中國星網(wǎng)網(wǎng)絡(luò)創(chuàng)新研究院，北京 100029；2.中國科學(xué)院國家天文臺，北京 100012；3.中國科學(xué)院國家空間科學(xué)中心，北京 100190）

1 引言

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

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

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

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

2 星歷數(shù)據(jù)處理

星鏈?zhǔn)敲绹仗剿骷夹g(shù)公司的一個項目，太空探索技術(shù)公司計劃在2019年至2024年間在太空搭建由約1.2萬顆衛(wèi)星組成的“星鏈”網(wǎng)絡(luò)提供互聯(lián)網(wǎng)服務(wù)，其中1584 顆將部署在地球上空550千米處的近地軌道，并從2020年開始工作。據(jù)有關(guān)文件顯示，該公司還準(zhǔn)備再增加3 萬顆，使衛(wèi)星總量達到約4.2萬顆[11-13]。

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

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

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

3 衛(wèi)星機動的識別

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

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

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

式中，Δa為半長軸的變化量，Δath為半長軸變化的均值，σ為半長軸變化的標(biāo)準(zhǔn)差。

4 機動加速度的計算方法

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

4.1 二分迭代法

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

式中，f為衛(wèi)星的速度矢量，v為衛(wèi)星的速度標(biāo)量。

由該機動力產(chǎn)生的機動加速度計算公式為：

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

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

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

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

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

（3）根據(jù)機動時刻后的星歷數(shù)據(jù)進行軌道確定，獲得機動結(jié)束時刻的軌道根數(shù)；

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

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

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

4.2 能量轉(zhuǎn)換法

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

即，

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

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

s由數(shù)值法軌道預(yù)報積分得出，當(dāng)衛(wèi)星軌道為近圓軌道時，s可以近似為圓弧。這種情況下，s可以由以下公式簡單計算，

式中，r為機動時段衛(wèi)星掃過的角度。

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

5 實例分析

5.1 機動加速度計算

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

星歷文件名稱為：

MEME_44713_STARLINk-1007_3102130_Operational_1352064660_UNCLASSIFIED.txt

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

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

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

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

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

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

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

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

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

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

表2 機動信息表Table 2 Maneuver information

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

5.2 計算結(jié)果驗證

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

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

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

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

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

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

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

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

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

6 結(jié)論

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