大橢圓軌道軌控期間碰撞預警規避閾值設置方法研究

李翠蘭， 劉成軍， 歐陽琦， 孔 靜， 陳 明， 劉 靜

(1.中國科學院國家天文臺， 北京 100012； 2.中國科學院大學， 北京 100049；3.北京航天飛行控制中心， 北京 100094)

1 引言

隨著人類探索太空活動的開展，在近地空間留下了數以千萬計的空間碎片。 截至到2019 年10 月4 日，美國空間監測網監測并編目的在軌空間目標總數已經達到19 779 個，其中衛星5181個，碎片14 598 片[1]，且大部分軌道高度分布在400～2000 km 之間。 空間碎片對在軌航天器的安全構成了嚴重威脅，特別是對長期穿越空間碎片密集區的航天器極大地增加了碰撞風險。

隨著載人航天工程的發展，新一代載人飛船將面向中國近地空間站運營、載人登月等任務需求，向著多功能化、多任務目標適應性的方向發展[2]。 為驗證關鍵技術，新一代載人飛船首飛采用相對較高的橢圓軌道，與2014 年NASA 獵戶座EFT－1 飛船的測試飛行類似[3]，軌道遠地點高度可達8000 km，入軌后經過多次變軌將遠地點從300 km 抬高到8000 km，將穿越空間碎片密集區，對空間目標碰撞風險評估提出了新的挑戰。

目前在空間目標碰撞風險評估中廣泛采用碰撞概率作為評價指標[4]。 文獻[5－8]對區域方法與碰撞概率方法進行了比較分析，結果反映當軌道預報協方差精度較高時，碰撞概率方法比區域方法可靠性高。 李甲龍等[8]分析了不同誤差情況下Box 方法和碰撞概率方法之間的關系，當軌道精度較高時，碰撞概率方法可以減少碰撞預警中的虛警問題；當軌道精度較差時，碰撞概率方法可能會失效，需要結合Box 方法進行綜合評估，減少碰撞預警中的漏警問題。Chan[9]提出當航天器相對距離很小(小于100 m)時，建議結合最小距離進行安全風險評估。 在實際工程應用中，一般情況下先采用Box進行危險目標篩選，然后采用碰撞概率進行最終判定。 但在軌控頻繁階段，如每隔2～3 圈，執行一次軌道機動的航天器，帶控的軌道預報誤差不確定性較大，難以保證碰撞概率計算結果可靠性，通常采用Box 判據作為碰撞風險評估的最終判據。 Box 判據的可靠性與軌道預報精度密切相關。 軌道預報精度受軌道高度、偏心率等的影響，對于大橢圓軌道，近地點附近和遠地點附近軌道預報偏差量級差別較大，且軌道控制通常選擇在近、遠地點控制，導致近、遠地點的軌道偏差不確定性增加。

本文針對大橢圓軌道在控制期間使用Box 判據計算碰撞風險漏警和虛警率較高的問題，對影響Box 規避閾值設置的各種偏差進行分析，提出了一種動態Box 預警規避閾值設置方法，并對該方法進行了仿真驗證。

2 Box 預警規避閾值設置影響因素分析

大橢圓軌道在軌控期間，影響Box 預警規避閾值設置的主要因素有軌道預報偏差和控制偏差。

2.1 軌道預報偏差傳播特性

由于測量精度、軌道模型精度的限制以及空間環境的變化，航天器的初始時刻狀態、動力學模型等都存在不確定性。 而航天器軌道預報則是根據航天器初始時刻狀態進行動力學積分的過程。因此軌道預報偏差特性主要取決于初值偏差和動力學模型偏差。

2.1.1 初值偏差影響分析

在初始偏差影響分析過程中，假定軌道動力學模型是準確的，沒有預報偏差，從而得到初始偏差的傳播特性。

本文假定航天器軌道近地高度不變，約為300 km，遠地點高度分別為2000 km、3500 km、5000 km、6500 km、8000 km 的橢圓軌道，軌道預報時長為24 h。 假定橢圓軌道初值偏差如表1 所示，初值偏差預報偏差如圖1 所示。 從圖1 可以看出，初值偏差主要將導致軌道跡向(U)方向偏差的外包絡呈線性增長，垂跡方向(N)和法向(W)偏差呈周期性變化。 預報24 h 跡向偏差最大約為1 km、垂跡方向偏差最大偏差約為20 m、法向偏差最大偏差約為1.5 m。

圖1 軌道初值偏差預報偏差Fig.1 Orbital propagation error affected by initial orbital error

表1 橢圓軌道初值偏差Table 1 Initial orbital error of the large elliptic orbit

2.1.2 動力學模型影響分析

動力學模型偏差因攝動力性質的不同而影響不同，對于近地軌道，通常考慮表2 所示動力學模型和參數。 在動力學模型影響分析過程中，假定沒有初值誤差，以表2 所列動力學模型作為基準預報模型，外推24 h，作為基準軌道。 分別對表2所列的動力學模型進行了分析，并根據分析結果結合Box 預警規避閾值確定大橢圓軌道的預報動力學模型。

表2 軌道動力學模型及參數Table 2 Orbital dynamics model and parameters

1)地球非球形攝動力。 本文地球非球形引力分別采用32×32 階次和64×64 階次，外推24 h，與基準軌道比較，預報位置偏差如圖2 所示。 從圖2 可以看出，地球非球形引力對預報偏差的影響主要是沿軌道跡向的線性增長，且具有明顯周期性。 使用32×32 階次地球非球形引力，預報24 h跡向偏差最大約為500 m、垂跡方向最大偏差約為15 m、法向最大偏差約為10 m；使用64×64 階次地球非球形引力，預報24 h 跡向偏差最大約為30 m、垂跡方向最大偏差約為0.6 m、法向最大偏差約為0.4 m。 假定橢圓軌道在軌控期間的Box 預警規避閾值以國際空間站的48 h 預警閾值跡向(U)4 km、垂跡方向(N) 0.5 km、法向(W)4 km[10]作為紅色閾值。 相對Box 預警閾值跡向公里量級的標準，綜合考慮預報精度和計算效率的影響，在橢圓軌道預報過程中，地球非球形引力取64×64 階次模型效果最佳。

圖2 非球形引力模型偏差引起的軌道偏差Fig.2 Orbital propagation error affected by non-spherical gravitational perturbation model error

2)日月引力。 在不考慮日月引力的情況下，外推24 h，與基準軌道比較，預報位置偏差如圖3所示。 從圖3 可以看出，日月引力對預報偏差的影響主要是軌道跡向和垂跡方向偏差迅速增大，周期性明顯，外包絡呈線性增長。 預報24 h 跡向偏差最大約為450 m、垂跡方向最大偏差約為200 m、法向最大偏差約為100 m。 相對Box 預警規避閾值跡向和法向公里量級，垂跡方向0.5 km的標準，綜合考慮預報精度和計算效率的影響，在大橢圓軌道預報過程中，需考慮日月引力的影響。

圖3 日月引力模型偏差分析圖Fig.3 Orbital propagation error affected by Solar and Lunar gravitational perturbation model error

3)其他大行星引力。 在不考慮其他大行星引力的情況下，外推24 h，與基準軌道比較，預報位置偏差如圖4 所示。 從圖4 可以看出，預報24 h，其他大行星引力引起的預報偏差約10－4量級。 遠遠小于Box 預警規避閾值的標準，在大橢圓軌道預報過程中，可以忽略不計。

圖4 其他天體引力模型偏差分析圖Fig.4 Orbital propagation error affected by other big planets gravitational perturbation model error

4)太陽光壓。 在太陽光壓考慮10%偏差的情況下，外推24 h，與基準軌道比較，預報位置偏差如圖5 所示。 從圖5 可以看出，預報24 h，太陽光壓考慮10%偏差的情況下，引起的軌道偏差約0.5 m，遠遠小于Box 預警規避閾值的標準。 因此，在大橢圓軌道預報過程中，光壓模型的誤差可以忽略不計。

圖5 太陽光壓偏差分析圖Fig.5 Orbital propagation error affected by sunlight pressure perturbation model error

5)大氣阻力。 在大氣阻力考慮10%偏差的情況下，其他模型采用表2 所列模型，外推24 h，與基準軌道比較，預報位置偏差如圖6 所示。 從圖6 可以看出，預報24 h，大氣阻力考慮10%偏差的情況下，引起的軌道偏差最大約為30 m，遠遠小于Box 預警規避閾值的標準。 因此，在大橢圓軌道預報過程中，大氣阻力模型誤差可以忽略不計。

圖6 大氣阻力偏差分析圖Fig.6 Orbital propagation error affected by atmosphere drag model error

2.1.3 構建模型

對不同高度大橢圓軌道預報偏差分析可知，大橢圓軌道通過合理選取動力學模型，可降低動力學模型偏差，因此影響軌道預報偏差主要來源是初值偏差。 為了模擬大橢圓軌道預報誤差傳播規律，采用蒙特卡洛仿真方法進行10 000 次仿真打靶。 仿真過程如下:

1)假定初值誤差服從正態分布，以表1 所示的初值誤差作為輸入；

2)動力學模型考慮64×64 階次地球非球形引力，日月引力，光壓和大氣阻，光壓和大氣阻力均采用10%誤差；

3)計算打靶軌道與基準軌道在U、N、W3 個方向的軌道預報偏差；

4)對軌道預報偏差進行統計，計算U、N、W3個方向的標準差；

5)根據標準差分布規律，進行建模。

以5000 km 高度大橢圓軌道預報偏差為例，U、N、W3 個方向的標準差仿真結果如圖7所示。

由圖7 可知，軌道預報誤差傳播規律明顯。跡向預報偏差發散較快，外包絡隨時間的平方發散，預報時間越長周期性越明顯；垂跡和法向預報誤差較小，呈明顯的周期性變化，垂跡方向位置偏差小于15 m，法向位置偏差小于1 m。 考慮到跡向周期性變化遠小于其線性變化，垂跡方向和法向值較小，可以忽略周期性變化的影響。 因此，在采用統計擬合方式時，時間作為自變量，預報偏差的外包絡的每圈最大值作為因變量，進行多項式擬合，可建立U、N、W3 個方向的軌道預報誤差解析模型，如式(1)所示:

式中， ΔUr、ΔNr、ΔWr為航天器在U、N、W3方向偏差，a為相關系數，aU0、aU1aU2分別為U方向的零次項系數、一次項系數和二次項系數，aN0、aN1分別為N方向的零次項系數、一次項系數，aW0、aW1分別為W方向的零次項系數、一次項系數，均可通過多項式擬合求解，t為預報時間(單位為h)。

以5000 km 高度大橢圓軌道為例，可建立式(2)所示軌道預報誤差解析模型。

2.2 軌控偏差傳播特性

為了節省燃料，軌道面內控制通常選擇在近、遠地點進行，均采用調整沿跡速度增量實現。 本文在軌控偏差分析中主要研究沿跡速度增量的控制偏差對軌道預報偏差的影響。

2.2.1 控制軌道偏差影響分析

分別對5 種高度(2000 km、3500 km、5000 km、6500 km、8000 km)的橢圓軌道進行了控制預報偏差分析。 假定軌道控制均在近地點，且調整沿跡速度增量執行，以理論控后軌道為基準軌道，在理論控制速度增量上疊加控制偏差后的軌道作為比較軌道，以0.1 m/s 為間隔，0.5 m/s為軌控偏差的上限，分別計算5 種高度的橢圓軌道在不同速度增量下U、N、W方向的位置偏差。 圖8 給出了在0.1 m/s 控制偏差情況下，5種高度的橢圓軌道預報偏差。 圖9 為5000 km 高度下不同控制偏差情況下的軌道預報偏差。

圖8 不同高度軌道預報偏差圖Fig.8 Orbital propagation error affected by control error for different altitude of objects

圖9 不同控制偏差情況下軌道偏差圖Fig.9 Orbital propagation error affected by different control error

圖8 和圖9 可以看出，由控制偏差引起的軌道預報偏差規律性較強，U方向位置偏差隨軌控偏差和預報時長的增加逐漸增大，N方向位置偏差存在明顯的周期性，最大值隨軌控偏差的增大而增大，W方向位置偏差較小，相對W方向碰撞預警規避閾值可忽略不計。 軌控偏差對U方向和N方向預報偏差影響較大，均大于Box 預警規避閾值。 因此，為了更準確地確定空間目標碰撞預警規避閾值，需進一步分析軌控偏差與預報偏差的關系。

從目前的樣本來看，U方向和N方向控制偏差引起軌道預報偏差的傳播具有一定的規律，將軌控偏差和時間作為自變量，最大位置偏差作為因變量，進行線性擬合，可以建立如公式(3)所示的軌控偏差預報模型。

式中，ΔUv、ΔNv為控制誤差引起的U、N、W3方向的預報誤差，δV為控制參數中的速度增量偏差，a為相關系數，aUv0、aUv1分別為U方向的零次項系數、一次項系數，aNv0、aNv1分別為N方向的零次項系數、一次項系數，可通過多項式擬合求解，t為預報時間，T為軌道周期。

以5000 km 橢圓軌道為例，表3 給出了U方向外包絡軌道預報偏差隨軌控偏差和時間的變化。 表4 給出了N方向最大軌道預報偏差隨軌控偏差和時間的變化。 通過多項式擬合，U方向和N方向軌道預報偏差與控制偏差之間存在如式(4)所示關系。

表3 U 方向軌道預報偏差隨軌控偏差和時間的變化Table 3 Orbital propagation error in U direction affected by velocity error and time

表4 N 方向軌道預報偏差隨軌控偏差和時間的變化Table 4 Orbital propagation error in N direction affected by velocity error and time

3 Box 動態閾值設置研究

3.1 Box 動態閾值設置方法提出

由第2 節的分析可知，軌道偏差和控制偏差對軌道預報的影響存在一定的規律性，可以通過多項式擬合的方式實現實時計算。 本文在設置空間目標碰撞預警規避閾值時，提出將軌道偏差和控制偏差產生的軌道預報偏差疊加，實時生成碰撞預警規避閾值的方法。

假設ΔUr1、ΔNr1、ΔWr1為主航天器3 方向由初值偏差和力模型偏差引起的預報偏差， ΔUv、ΔNv、ΔWv為由控制偏差引起的主航天器三向預報偏差，ΔUr2、ΔNr2、ΔWr2為由初值偏差和力模型偏差引起的空間碎片三向預報偏差，則動態BOX預警閾值如式(5)所示:

在工程應用中，空間碎片的預報偏差難以逐一求解，且主航天器的預報偏差依據第2 節對初值偏差和動力學模型偏差的分析，預報偏差小于4 km×0.5 km×4 km。 因此，在碰撞預警規避閾值設置時，在公式(5)的基礎上進行了簡化，在碰撞預警規避閾值4 km×0.5 km×4 km 的基礎上，在U和N方向上疊加由控制偏差導致的軌道預報偏差，生成Box 動態碰撞預警規避閾值(4＋ΔUv) km× (0.5＋ΔNv) km×4 km。

3.2 仿真驗證

依據3.1 節提出的簡化空間目標碰撞預警規避閾值設置方法，空間碎片數據使用公開的雙行根數數據，對5000 km 的橢圓軌道在12 h 內的碰撞風險進行了計算，結果見圖10。 圖10 中粗淺藍色線為以4 km 作為BOX 預警規避閾值，黑色實線為0.1 m/s 控制偏差對應的動態BOX 預警規避閾值，綠色虛線為0.2 m/s 控制偏差對應的動態BOX 預警規避閾值，紅色點劃線為0.3 m/s控制偏差對應的動態BOX 預警規避閾值，藍色粗實線為0.4 m/s 控制偏差對應的動態BOX 預警規避閾值，紫色點虛線為0.5 m/s 控制偏差對應的動態BOX 預警規避閾值，黑色空心點為危險事件。 在線下方標明為危險事件，線上方為虛警事件。 從圖中可以看出，傳統的Box 預警規避門限，即粗淺藍色線為預警門限，80%的碰撞事件發生在上方，如果實際軌控偏差為0.5 m/s，則會導致較多漏警事件；同樣如果采用0.5 m/s 控制偏差，即紫色虛線作為預警門限，則U方向預警規避門限將達到100 km，N方向預警規避門限將達到4 km，會有較多目標進入Box 預警門限，而實際軌控偏差小于0.5 m/s，則會導致較多虛警事件。在實際工程任務中，可根據實時標定的軌控偏差，計算動態Box 預警閾值，降低虛警事件和漏警事件。

圖10 動態Box 預警規避閾值分析圖Fig.10 Analysis of dynamic Box collision risk threshold

4 結論

1)對Box 預警規避閾值設置影響因素進行分析，分析結果表明軌道初值偏差、動力學模型偏差和控制偏差導致軌道預報偏差在U、N、W方向具有隨時間變化的傳播特性；

2)本文提出了動態Box 預警閾值設置方法，即將各種偏差導致的軌道預報偏差疊加作為Box預警閾值的方法，進一步對閾值設置方法進行了簡化，即在現有碰撞預警規避閾值4 km×0.5 km×4 km 的基礎上，在U和N方向上疊加由控制偏差導致的軌道預報偏差，生成Box 動態碰撞預警規避閾值(4 ＋ΔUv) km× (0.5＋ΔNv) km×4 km。 動態Box 預警閾值設置方法極大地降低使用Box 預警判據計算的碰撞風險虛警率和漏警率。