大橢圓軌道軌控期間碰撞預(yù)警規(guī)避閾值設(shè)置方法研究

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

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

1 引言

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

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

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

本文針對大橢圓軌道在控制期間使用Box 判據(jù)計算碰撞風(fēng)險漏警和虛警率較高的問題，對影響B(tài)ox 規(guī)避閾值設(shè)置的各種偏差進(jìn)行分析，提出了一種動態(tài)Box 預(yù)警規(guī)避閾值設(shè)置方法，并對該方法進(jìn)行了仿真驗證。

2 Box 預(yù)警規(guī)避閾值設(shè)置影響因素分析

大橢圓軌道在軌控期間，影響B(tài)ox 預(yù)警規(guī)避閾值設(shè)置的主要因素有軌道預(yù)報偏差和控制偏差。

2.1 軌道預(yù)報偏差傳播特性

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

2.1.1 初值偏差影響分析

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

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

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

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

2.1.2 動力學(xué)模型影響分析

動力學(xué)模型偏差因攝動力性質(zhì)的不同而影響不同，對于近地軌道，通?？紤]表2 所示動力學(xué)模型和參數(shù)。 在動力學(xué)模型影響分析過程中，假定沒有初值誤差，以表2 所列動力學(xué)模型作為基準(zhǔn)預(yù)報模型，外推24 h，作為基準(zhǔn)軌道。 分別對表2所列的動力學(xué)模型進(jìn)行了分析，并根據(jù)分析結(jié)果結(jié)合Box 預(yù)警規(guī)避閾值確定大橢圓軌道的預(yù)報動力學(xué)模型。

表2 軌道動力學(xué)模型及參數(shù)Table 2 Orbital dynamics model and parameters

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

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

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

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

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

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

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

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

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

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

2.1.3 構(gòu)建模型

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

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

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

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

4)對軌道預(yù)報偏差進(jìn)行統(tǒng)計，計算U、N、W3個方向的標(biāo)準(zhǔn)差；

5)根據(jù)標(biāo)準(zhǔn)差分布規(guī)律，進(jìn)行建模。

以5000 km 高度大橢圓軌道預(yù)報偏差為例，U、N、W3 個方向的標(biāo)準(zhǔn)差仿真結(jié)果如圖7所示。

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

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

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

2.2 軌控偏差傳播特性

為了節(jié)省燃料，軌道面內(nèi)控制通常選擇在近、遠(yuǎn)地點進(jìn)行，均采用調(diào)整沿跡速度增量實現(xiàn)。 本文在軌控偏差分析中主要研究沿跡速度增量的控制偏差對軌道預(yù)報偏差的影響。

2.2.1 控制軌道偏差影響分析

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

圖8 不同高度軌道預(yù)報偏差圖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 可以看出，由控制偏差引起的軌道預(yù)報偏差規(guī)律性較強(qiáng)，U方向位置偏差隨軌控偏差和預(yù)報時長的增加逐漸增大，N方向位置偏差存在明顯的周期性，最大值隨軌控偏差的增大而增大，W方向位置偏差較小，相對W方向碰撞預(yù)警規(guī)避閾值可忽略不計。 軌控偏差對U方向和N方向預(yù)報偏差影響較大，均大于Box 預(yù)警規(guī)避閾值。 因此，為了更準(zhǔn)確地確定空間目標(biāo)碰撞預(yù)警規(guī)避閾值，需進(jìn)一步分析軌控偏差與預(yù)報偏差的關(guān)系。

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

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

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

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

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

3 Box 動態(tài)閾值設(shè)置研究

3.1 Box 動態(tài)閾值設(shè)置方法提出

由第2 節(jié)的分析可知，軌道偏差和控制偏差對軌道預(yù)報的影響存在一定的規(guī)律性，可以通過多項式擬合的方式實現(xiàn)實時計算。 本文在設(shè)置空間目標(biāo)碰撞預(yù)警規(guī)避閾值時，提出將軌道偏差和控制偏差產(chǎn)生的軌道預(yù)報偏差疊加，實時生成碰撞預(yù)警規(guī)避閾值的方法。

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

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

3.2 仿真驗證

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

圖10 動態(tài)Box 預(yù)警規(guī)避閾值分析圖Fig.10 Analysis of dynamic Box collision risk threshold

4 結(jié)論

1)對Box 預(yù)警規(guī)避閾值設(shè)置影響因素進(jìn)行分析，分析結(jié)果表明軌道初值偏差、動力學(xué)模型偏差和控制偏差導(dǎo)致軌道預(yù)報偏差在U、N、W方向具有隨時間變化的傳播特性；

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