徐家輝 胡 敏 王許煜 李玖陽(yáng) 云朝明

航天工程大學(xué)，北京 101416

0 引言

隨著1957年第1顆衛(wèi)星發(fā)射升空，各國(guó)都開(kāi)始制定衛(wèi)星的發(fā)射計(jì)劃，導(dǎo)致衛(wèi)星和火箭數(shù)目大幅增加，同時(shí)空間目標(biāo)的碰撞加劇了空間碎片[1]增長(zhǎng)的速度。截至2020年2月，尺寸大于10cm的可編目的空間碎片數(shù)目為45122個(gè)。空間碎片的分布不是均勻的，碎片的擁擠程度由可用碎片的空間密度[2]描述。最早由?pik、Wetherill等提出，用來(lái)描述小行星帶內(nèi)天體的分布狀態(tài)[3]，之后，凱斯勒等根據(jù)碎片的軌道特性給出了物體空間密度的計(jì)算方法。圖1給出了碎片的空間密度隨軌道高度的分布情況，中軌道區(qū)域碎片空間密度存在峰值，這正是全球衛(wèi)星導(dǎo)航系統(tǒng)所處的區(qū)域。停留在導(dǎo)航衛(wèi)星附近的廢棄衛(wèi)星及上面級(jí)將給導(dǎo)航衛(wèi)星的安全運(yùn)行帶來(lái)潛在風(fēng)險(xiǎn)。

圖1 碎片空間密度隨軌道高度變化圖

中軌道區(qū)域是導(dǎo)航衛(wèi)星的主要運(yùn)行區(qū)域，其空間密度隨著導(dǎo)航衛(wèi)星的發(fā)射而增加。目前，中軌道區(qū)域還沒(méi)有明確的到壽處置原則，各星座的到壽處置方式也不相同。為此，給導(dǎo)航星座衛(wèi)星的在軌運(yùn)行帶來(lái)安全隱患。另一方面，在地球非球形、日月三體引力和太陽(yáng)光壓等攝動(dòng)力的影響下，廢棄衛(wèi)星的軌道根數(shù)會(huì)發(fā)生長(zhǎng)期變化，有可能導(dǎo)致廢棄衛(wèi)星穿越到導(dǎo)航衛(wèi)星在軌運(yùn)行軌道，給在軌導(dǎo)航衛(wèi)星帶來(lái)碰撞[6]風(fēng)險(xiǎn)。歐空局碎片環(huán)境模型MASTER-2005預(yù)測(cè)將會(huì)有60000個(gè)尺寸大于1cm的空間物體穿過(guò)導(dǎo)航星座區(qū)域，導(dǎo)航衛(wèi)星的在軌安全運(yùn)行問(wèn)題研究顯得較為迫切。

本文以中軌道導(dǎo)航衛(wèi)星為研究對(duì)象，采用半分析方法分離出短周期項(xiàng)，使得積分步長(zhǎng)設(shè)置為1天，兼顧軌道預(yù)報(bào)精度和計(jì)算效率，以此建立了中軌道衛(wèi)星長(zhǎng)期預(yù)報(bào)模型?；陂L(zhǎng)期預(yù)報(bào)模型，對(duì)導(dǎo)航星座衛(wèi)星及上面級(jí)進(jìn)行長(zhǎng)期演化數(shù)值分析，以分析Galileo、BDS、GPS和GLONASS在當(dāng)前處置情況下其廢棄衛(wèi)星和上面級(jí)是否穿越到星座運(yùn)行高度及其所耗時(shí)間。

1 中軌道導(dǎo)航衛(wèi)星長(zhǎng)期預(yù)報(bào)模型

衛(wèi)星軌道預(yù)報(bào)[7-8]根據(jù)時(shí)間尺度可分為長(zhǎng)期和短期預(yù)報(bào)，基本思想都是基于當(dāng)前時(shí)刻的軌道狀態(tài)及攝動(dòng)力對(duì)軌道狀態(tài)積分外推。劉林提出了一種衛(wèi)星軌道預(yù)報(bào)的分析方法，用平均根數(shù)表示t時(shí)刻衛(wèi)星的位置和速度，并直接用衛(wèi)星直角坐標(biāo)的位置和速度分量表示地球非球形攝動(dòng)的周期項(xiàng)，可以有效避免由于小偏心率和軌道傾角而出現(xiàn)的奇點(diǎn)問(wèn)題。王大為對(duì)不同攝動(dòng)力進(jìn)行半分析法處理，并將完整攝動(dòng)力模型的半分析方法應(yīng)用于不同類型軌道，以研究在不同初始軌道根數(shù)半長(zhǎng)軸和軌道傾角的組合下，半分析法在100年內(nèi)軌道預(yù)報(bào)的精度。張斌斌以碎片所在的軌道高度和面質(zhì)比對(duì)碎片進(jìn)行分類，通過(guò)平均法去除短周期項(xiàng)，獲得空間碎片平均分點(diǎn)根數(shù)變化率的一階變化率方程，在保證軌道預(yù)報(bào)進(jìn)度的情況下顯著提高了計(jì)算效率。

1.1 攝動(dòng)力模型

中軌道區(qū)域的軌道動(dòng)力學(xué)與低軌道區(qū)域的軌道動(dòng)力學(xué)有很大不同，主要是由于不需要考慮大氣阻力[12]。在近地空間的這一區(qū)域，其它擾動(dòng)(例如地球非球形、三體引力和太陽(yáng)光壓)主導(dǎo)并確定了軌道的演化。導(dǎo)航衛(wèi)星在軌過(guò)程還受到其它微小攝動(dòng)力的作用，但其加速度的量級(jí)太小，在中軌道導(dǎo)航衛(wèi)星軌道長(zhǎng)期預(yù)報(bào)時(shí)可以忽略。因此，中軌道區(qū)域的攝動(dòng)力模型如式(1)：

asd=a0+ans+aS+aM+asr

(1)

式中：a0表示地球中心引力項(xiàng)；ans表示地球非球形攝動(dòng)項(xiàng)；asr表示太陽(yáng)光壓攝動(dòng)項(xiàng)；aS和aM分別表示太陽(yáng)、月球引力攝動(dòng)項(xiàng)。

1.2 長(zhǎng)期預(yù)報(bào)模型

對(duì)中軌道區(qū)域航天器軌道進(jìn)行長(zhǎng)期預(yù)報(bào)，需要按一定的積分步長(zhǎng)對(duì)航天器的軌道根數(shù)進(jìn)行外推。航天器在攝動(dòng)力作用下，其軌道狀態(tài)的變化較為復(fù)雜，包含短周期、長(zhǎng)周期和長(zhǎng)期變化3部分。在對(duì)航天器進(jìn)行軌道長(zhǎng)期預(yù)報(bào)時(shí)，一般只考慮軌道的長(zhǎng)周期和長(zhǎng)期變化，即平均軌道根數(shù)的變化，需要對(duì)軌道狀態(tài)平均化以去除短周期變化。

采用分點(diǎn)根數(shù)(Equinoctial Elements)描述航天器的軌道狀態(tài)可以有效避免運(yùn)動(dòng)方程出現(xiàn)奇異的情況。本文研究對(duì)象為中軌道導(dǎo)航衛(wèi)星，因此只考慮順行軌道，分點(diǎn)根數(shù)記為(a，h，k，p，q，λ)，與開(kāi)普勒軌道根數(shù)的轉(zhuǎn)換關(guān)系如式(2)：

(2)

在慣性坐標(biāo)系下航天器的運(yùn)動(dòng)方程為：

(3)

式中：μ為地球引力常數(shù)，R是一種類似勢(shì)能的函數(shù)，是由保守?cái)z動(dòng)力得到的攝動(dòng)函數(shù)，fnon為航天器受到的非保守?cái)z動(dòng)力加速度。

運(yùn)動(dòng)方程(3)是由航天器的位置和速度描述的，將式(3)轉(zhuǎn)化為軌道根數(shù)變化率隨軌道根數(shù)變化的參數(shù)化運(yùn)動(dòng)方程(Variation of Parameters Equations of Motion)。Danielson等[14]用(a1,a2,…,a6)分別表示分點(diǎn)根數(shù)(a，h，k，p，q，λ)，參數(shù)化運(yùn)動(dòng)方程用分點(diǎn)根數(shù)描述的形式如下：

(4)

式中：i=1,2,3,4,5,6；n為衛(wèi)星運(yùn)行的平均角速度；δi6為克羅內(nèi)克函數(shù)，其值為(0,0,0,0,0,1)。

參數(shù)化運(yùn)動(dòng)方程(4)軌道根數(shù)變化率由3部分組成，右式第1項(xiàng)為二體引力部分；第2項(xiàng)為非保守?cái)z動(dòng)部分，也叫項(xiàng)；第3部分為保守?cái)z動(dòng)部分，也叫拉格朗日項(xiàng)。對(duì)不同攝動(dòng)力求平均攝動(dòng)函數(shù)時(shí)，可以得到保守力攝動(dòng)作用下解析的平均攝動(dòng)函數(shù)；而非保守?cái)z動(dòng)力平均攝動(dòng)函數(shù)只能通過(guò)在一個(gè)周期內(nèi)求平均的方法獲得，從而得到不同攝動(dòng)作用下平均分點(diǎn)根數(shù)的一階變化率。平均攝動(dòng)函數(shù)的具體求解還需要考慮很多因素，具體結(jié)果可參考Danielson對(duì)半分析模型的總結(jié)報(bào)告。最后，將不同攝動(dòng)作用下的分點(diǎn)根數(shù)一階變化率相加，再利用數(shù)值積分方法對(duì)航天器軌道狀態(tài)外推更新，可以計(jì)算出航天器預(yù)報(bào)時(shí)刻的平均軌道根數(shù)。

2 中軌道衛(wèi)星長(zhǎng)期預(yù)報(bào)模型精度驗(yàn)證

利用半分析方法軌道長(zhǎng)期預(yù)報(bào)模型，對(duì)中軌道區(qū)域航天器的軌道根數(shù)進(jìn)行長(zhǎng)期外推更新。在軌道長(zhǎng)期預(yù)報(bào)模型的精度驗(yàn)證上，與STK中的LOP模型進(jìn)行對(duì)比，以評(píng)估其精度性能。選取中軌道區(qū)域GLONASS衛(wèi)星的典型軌道進(jìn)行100年長(zhǎng)期預(yù)報(bào)，選取的軌道初始參數(shù)在表1中給出。

表1 選取軌道長(zhǎng)期預(yù)報(bào)的初始軌道根數(shù)

根據(jù)導(dǎo)航衛(wèi)星的特點(diǎn)，衛(wèi)星的質(zhì)量設(shè)為1000kg，截面積為20m2。半分析法軌道長(zhǎng)期預(yù)報(bào)模型的攝動(dòng)力包括地球非球形攝動(dòng)力、日月三體引力和太陽(yáng)光壓。重力場(chǎng)勢(shì)函數(shù)的階數(shù)為4×4，太陽(yáng)光壓系數(shù)取為1.3。在利用半分析法軌道長(zhǎng)期預(yù)報(bào)模型對(duì)軌道根數(shù)進(jìn)行外推更新時(shí)，采用4階龍格-庫(kù)塔積分器求解太陽(yáng)光壓攝動(dòng)力作用下的平均根數(shù)的變化率，積分步長(zhǎng)設(shè)置為0.1rad。軌道根數(shù)積分外推則采用4階Adams-Cowell預(yù)估校正方法進(jìn)行求解，積分步長(zhǎng)設(shè)置為1天。STK中LOP模型的攝動(dòng)力和參數(shù)設(shè)置與半分析法軌道長(zhǎng)期預(yù)報(bào)模型一樣。

表1中航天器的軌道狀態(tài)與中軌道BDS衛(wèi)星相符，對(duì)軌道進(jìn)行100年向后推演，并對(duì)比LOP長(zhǎng)期預(yù)報(bào)模型，各軌道根數(shù)演化結(jié)果如圖2～5所示。圖中虛線表示LOP長(zhǎng)期預(yù)報(bào)模型計(jì)算結(jié)果，實(shí)線表示半分析法軌道長(zhǎng)期預(yù)報(bào)模型的計(jì)算結(jié)果，點(diǎn)劃線為兩模型偏差的從絕對(duì)值。圖2為半長(zhǎng)軸100年演化結(jié)果，可以看出半分析法軌道長(zhǎng)期預(yù)報(bào)模型與LOP模型的半長(zhǎng)軸的長(zhǎng)期項(xiàng)變化為0，LOP模型的短周期項(xiàng)振幅較大，半長(zhǎng)軸最大偏差控制在4km以內(nèi)。圖3為軌道偏心率長(zhǎng)期演化結(jié)果，最大偏差控制在0.0012以內(nèi)。圖4為軌道傾角的長(zhǎng)期演化結(jié)果，呈現(xiàn)長(zhǎng)期增長(zhǎng)短期震蕩的變化趨勢(shì)，軌道傾角最大偏差小于1°。

圖2 半長(zhǎng)軸在100年內(nèi)的長(zhǎng)期演化結(jié)果

圖3 偏心率在100年內(nèi)的長(zhǎng)期演化結(jié)果

圖4 軌道傾角在100年內(nèi)的長(zhǎng)期演化結(jié)果

升交點(diǎn)赤經(jīng)100年演化結(jié)果如圖5所示，半分析法軌道長(zhǎng)期預(yù)報(bào)模型的演化結(jié)果與LOP模型的變化趨勢(shì)是一致的，升交點(diǎn)赤經(jīng)的偏差隨著演化時(shí)間的增長(zhǎng)而出現(xiàn)積累，100年內(nèi)的最偏差控制在20°以內(nèi)。

圖5 升交點(diǎn)赤經(jīng)在100年內(nèi)的長(zhǎng)期演化結(jié)果

3 長(zhǎng)期演化安全性分析

中軌道區(qū)域是三3廣的區(qū)域，雖然其碎片平均空間密度還較小，但在導(dǎo)航星座軌道高度附近有空間密度峰值，給導(dǎo)航衛(wèi)星的長(zhǎng)期在軌運(yùn)行帶來(lái)安全隱患，因此需要分析導(dǎo)航星座長(zhǎng)期演化的安全性。中軌道導(dǎo)航星座長(zhǎng)期演化的安全性主要分析廢棄衛(wèi)星和上面級(jí)100年演化時(shí)間內(nèi)是否會(huì)穿越到運(yùn)行衛(wèi)星高度。如果會(huì)，則穿越到運(yùn)行衛(wèi)星高度所耗時(shí)越長(zhǎng)，對(duì)在軌運(yùn)行衛(wèi)星的潛在安全威脅越小。

四大全球衛(wèi)星導(dǎo)航系統(tǒng)星座構(gòu)型均是Walker星座，軌道面有3個(gè)或6個(gè)，不同軌道面上的衛(wèi)星在攝動(dòng)力作用下其軌道根數(shù)的變化趨勢(shì)不盡相同。在軌運(yùn)行衛(wèi)星仍然保持穩(wěn)定的星座構(gòu)型，可以按軌道面分別進(jìn)行長(zhǎng)期演化分析。上面級(jí)及廢棄衛(wèi)星因已經(jīng)偏離運(yùn)行軌道面且分布較為分散，不需要按軌道面分類。

3.1 GPS

GPS星座共發(fā)射72顆衛(wèi)星，其中在軌運(yùn)行衛(wèi)星31顆，廢棄衛(wèi)星41顆。圖6給出了GPS廢棄衛(wèi)星100年內(nèi)遠(yuǎn)、近地點(diǎn)地心距的變化，從圖中可以看到GPS廢棄衛(wèi)星因?yàn)槌跏计穆瘦^大，軌道高度范圍也大，部分衛(wèi)星其遠(yuǎn)、近地點(diǎn)已經(jīng)穿越到BDS和GLONASS運(yùn)行區(qū)域，40年左右會(huì)穿越的Galileo運(yùn)行區(qū)域，100年內(nèi)穩(wěn)定性較差。

圖6 GPS廢棄衛(wèi)星遠(yuǎn)近地點(diǎn)地心距100年長(zhǎng)期演化

GPS上面級(jí)包括DELTA 4 R/B和ATLAS 5 CENTAUR R/B兩種型號(hào)。同時(shí)，上面級(jí)的處置策略采用抬高軌道高度的方式，部分上面級(jí)已經(jīng)穿越到BDS導(dǎo)航星座運(yùn)行區(qū)域，影響B(tài)DS衛(wèi)星的在軌安全運(yùn)行。對(duì)GPS上面級(jí)進(jìn)行100年演化，其遠(yuǎn)、近地點(diǎn)地心距演化如圖7，GPS上面級(jí)穿越到BDS和GPS運(yùn)行區(qū)域較為顯著，45年左右將會(huì)穿越到GLONASS和Galileo運(yùn)行區(qū)域。

圖7 GPS上面級(jí)遠(yuǎn)近地點(diǎn)地心距100年長(zhǎng)期演化

3.2 GLONASS

GLONASS區(qū)域共發(fā)射了134顆衛(wèi)星，其中在軌運(yùn)行衛(wèi)星25顆，廢棄衛(wèi)星109顆。對(duì)109顆GLONASS廢棄衛(wèi)星進(jìn)行長(zhǎng)期演化，圖8給出了其100年時(shí)間內(nèi)遠(yuǎn)、近地點(diǎn)地心距的變化，從中可以看到GLONASS廢棄衛(wèi)星因?yàn)槌跏计穆瘦^小，100年內(nèi)其偏心率增長(zhǎng)也較小。GLONASS廢棄衛(wèi)星約65年左右才能穿越到GPS運(yùn)行區(qū)域，約95年左右才能穿越到BDS運(yùn)行區(qū)域，100年內(nèi)不會(huì)穿越到Galileo運(yùn)行區(qū)域，穩(wěn)定性較好。

圖8 GLONASS廢棄衛(wèi)星遠(yuǎn)近地點(diǎn)地心距100年長(zhǎng)期演化

GLONASS上面級(jí)包括SL-12 R/B(2)和FREGAT R/B兩種型號(hào)，早期GLONASS上面級(jí)沒(méi)有離軌處置，因此還滯留在運(yùn)行軌道附近，2011年后采用了抬高軌道高度的處置策略。目前，GLONASS上面級(jí)只穿越到GLONASS運(yùn)行區(qū)域，尚未穿越到其它導(dǎo)航星座運(yùn)行區(qū)域。對(duì)GLONASS上面級(jí)進(jìn)行100年軌道演化，其遠(yuǎn)、近地點(diǎn)地心距演化結(jié)果如圖9，上面級(jí)穿越到GPS運(yùn)行區(qū)域需要40年，80年左右將會(huì)穿越到BDS運(yùn)行區(qū)域。

圖9 GLONASS上面級(jí)遠(yuǎn)近地點(diǎn)地心距100年長(zhǎng)期演化

3.3 GALILEO

Galileo系統(tǒng)共發(fā)射了26顆衛(wèi)星，早期發(fā)射的2顆衛(wèi)星已經(jīng)偏離了運(yùn)行軌道面。

在Galileo導(dǎo)航星座附近還有因發(fā)射任務(wù)而產(chǎn)生的上面級(jí)，Galileo上面級(jí)包括FREGAT R/B和ARIANE 5 R/B兩種型號(hào)。同時(shí)，Galileo上面級(jí)的處置情況也不同，早期上面級(jí)采取抬高軌道高度的處置方式，2016年后則將上面級(jí)處置在運(yùn)行軌道高度下方。對(duì)Galileo上面級(jí)進(jìn)行100年長(zhǎng)期演化，其遠(yuǎn)、近地點(diǎn)地心距隨時(shí)間演化結(jié)果如圖10。歸功于Galileo較小的初始偏心率控制，Galileo上面級(jí)的廢棄軌道很穩(wěn)定，很少穿越到Galileo運(yùn)行軌道高度區(qū)域，更不會(huì)穿越到其它導(dǎo)航星座運(yùn)行軌道區(qū)域。

圖10 Galileo上面級(jí)遠(yuǎn)近地點(diǎn)地心距100年長(zhǎng)期演化

3.4 BDS

BDS在中軌道區(qū)域共發(fā)射29顆衛(wèi)星，其中有一顆衛(wèi)星已經(jīng)偏離了運(yùn)行軌道面。BDS上面級(jí)的型號(hào)為YZ-1 R/B，其處置方式為抬高軌道高度，圖11給出了BDS上面級(jí)100年內(nèi)遠(yuǎn)、近地點(diǎn)的演化情況，目前，已經(jīng)有部分BDS上面級(jí)穿越到Galileo和BDS星座的運(yùn)行區(qū)域。對(duì)于尚未穿越到運(yùn)行區(qū)域的BDS上面級(jí)，在100年內(nèi)，隨著偏心率的增大，上面級(jí)逐漸穿越到BDS，Galileo、GPS和GLONASS運(yùn)行區(qū)域。

圖11 BDS上面級(jí)遠(yuǎn)近地點(diǎn)地心距100年長(zhǎng)期演化

從上述分析可知，衛(wèi)星的軌道不是穩(wěn)定不變的，在攝動(dòng)力作用下會(huì)不斷穿越到在軌衛(wèi)星軌道上。因此針對(duì)未來(lái)導(dǎo)航衛(wèi)星及上面級(jí)的到壽處置策略，需要控制廢棄軌道遠(yuǎn)、近地點(diǎn)高度，可以通過(guò)廢棄軌道高度和初始偏心率來(lái)控制，使航天器不能穿越到其它導(dǎo)航星座運(yùn)行區(qū)域。其次，不僅要關(guān)注其它導(dǎo)航星座對(duì)導(dǎo)航星座的干擾，也要防止對(duì)自身運(yùn)行衛(wèi)星的干擾，因此應(yīng)該劃分導(dǎo)航衛(wèi)星及其上面級(jí)廢棄軌道處置高度帶，同時(shí)控制初始偏心率，使廢棄衛(wèi)星長(zhǎng)期穩(wěn)定滯留在廢棄軌道。

4 結(jié)論

中軌道區(qū)域空間目標(biāo)數(shù)量的不斷增加，不僅使中軌道區(qū)域的軌道資源縮減，也使得在軌衛(wèi)星發(fā)生碰撞的可能性增大。本文建立了中軌道衛(wèi)星軌道長(zhǎng)期預(yù)報(bào)模型，分析了中軌道導(dǎo)航星座100年長(zhǎng)期演化的安全性。仿真結(jié)果表明：Galileo穩(wěn)定性最好，GPS對(duì)其它星座的安全影響最大，GLONASS則主要穿越到自身運(yùn)行區(qū)域，BDS穩(wěn)定性也較好。