張 煒，崔 文，田 鑫，雷 園，劉 興

(西安衛(wèi)星測控中心，西安 710043)

1 引言

天宮一號是我國第一個目標飛行器和空間實驗室，于2011年9月29日發(fā)射入軌，先后與神舟八號、神舟九號、神舟十號成功對接完成相關任務。2016年3月16日天宮一號目標飛行器正式終止數(shù)據(jù)服務，無法繼續(xù)維持正常運行軌道，在大氣阻力的影響下軌道不斷衰減，最終于2018年4月2日8時15分左右再入大氣層。軌道壽命末期天宮一號重約8 t，是我國非受控再入的第一個超大型航天器。大質(zhì)量航天器在再入過程中并不能被完全燒毀，仍有10%～40%的殘骸返回地球表面[1]，對地表的生命群體、建筑設施、生態(tài)環(huán)境等具有很大威脅，因此有必要進行精確的再入時間和落點預報[2]。

天宮一號的再入過程有以下幾個特點：

1) 再入過程中飛行器姿態(tài)經(jīng)歷了穩(wěn)定到不穩(wěn)定、緩慢翻滾到快速翻滾的變化；

2) 再入過程中太陽活動基本處于極低水平；

3) 再入的最后階段地磁環(huán)境處于持續(xù)擾動狀態(tài)。

因此，天宮一號的再入案例對于開展軌道動力學和再入預報研究具有特殊意義。

大氣阻力是即將再入航天器所受的最主要非保守攝動力，精確的面質(zhì)比及合理的大氣阻力特性建模是準確計算大氣阻力加速度、預報再入時間的關鍵[3]。對于絕大多數(shù)空間目標，其準確外形、質(zhì)量、姿態(tài)和表面材料等都是未知的，分別確定大氣阻力系數(shù)、迎風面積和質(zhì)量難度很大，因此通常引入彈道系數(shù)B進行統(tǒng)一處理，彈道系數(shù)B的定義為式(1)[4]：

(1)

其中，CD為大氣阻力系數(shù)，A為迎風面積，m為質(zhì)量。

大部分非受控再入航天器都處于姿態(tài)失控狀態(tài)，研究天宮一號在低軌甚至超低軌道情況下的軌道動力學特性，有必要先判別飛行器姿態(tài)，針對不同研究目的選擇不同姿態(tài)階段的數(shù)據(jù)。傳統(tǒng)的姿態(tài)判別方法是基于目標的雷達特征(包括窄帶雷達散射截面積、一維距離像、二維/三維雷達圖像等)[5]或光學特征[6]進行判別，這些判別方法對觀測條件具有較高的要求，且一個地基觀測弧段往往只有幾分鐘，在航天器翻滾周期大于觀測弧段時長的情況下往往難以準確判斷飛行器姿態(tài)變化，因此需要其他輔助判別方法。

本文從姿態(tài)變化引起迎風面積變化這一特點出發(fā)，提出一種基于彈道系數(shù)的天宮一號目標飛行器姿態(tài)判別方法：使用美國戰(zhàn)略司令部公開發(fā)布的雙行根數(shù)(Two Line Element，TLE)作為數(shù)據(jù)源，首先對TLE數(shù)據(jù)進行篩選，剔除精度較差的異常TLE，然后基于多組TLE數(shù)據(jù)使用最小二乘法計算彈道系數(shù)，對天宮一號的飛行姿態(tài)進行判別。

2 TLE數(shù)據(jù)預處理

本文選擇美國戰(zhàn)略司令部發(fā)布的TLE[7]作為數(shù)據(jù)源。為了比較數(shù)據(jù)服務終止前后彈道系數(shù)的變化，選擇2015年1月1日作為數(shù)據(jù)起始時間。由于數(shù)據(jù)服務終止前天宮一號可能發(fā)生軌道機動，發(fā)布的TLE可能包括精度差甚至錯誤結(jié)果，因此TLE使用前需先進行預處理，檢測軌道機動，并剔除異常TLE。根據(jù)大型空間目標TLE根數(shù)精度的特點，預處理按以下幾步進行[8]：

1) 剔除發(fā)布后立即被修正的TLE；

2) 剔除B*為負的TLE(B*為TLE中的阻力系數(shù)項)；

3) 軌道機動檢測；

4) 基于平運動變化剔除異常TLE。

本節(jié)主要介紹軌道機動檢測和基于平運動變化剔除異常TLE的方法。

低軌航天器的軌道機動一般分為軌道面內(nèi)機動和軌道面外機動；軌道面內(nèi)機動會引起平運動、半長軸、能量和(或)偏心率的突然變化，而軌道面外機動改變的則是衛(wèi)星軌道的傾角或升交點赤經(jīng)[9]。比較相鄰兩段數(shù)據(jù)的軌道特征量(半長軸、能量或傾角等)，若特征量出現(xiàn)較大變化，則認為兩段數(shù)據(jù)之間發(fā)生軌道機動[10]。圖1展示了基于TLE檢測軌道機動方法的基本原理，本文將連續(xù)6條TLE數(shù)據(jù)組成數(shù)據(jù)段，使用最小二乘法擬合對應一階多項式，依據(jù)連續(xù)兩段數(shù)據(jù)的多項式變化進行軌道機動檢測。檢測到2015年1月1日以后天宮一號共進行兩次軌道機動，一次發(fā)生在2015年4月10日，另一次發(fā)生在2015年12月16日，均為軌道面內(nèi)軌道機動。

圖1 軌道機動檢測的基本原理Fig.1 The basic principle of maneuver detection

基于平運動變化剔除異常TLE是軌道機動檢測方法的簡化和延伸。使用平運動n作為特征量，尾部數(shù)據(jù)段長度仍為6，而頭部數(shù)據(jù)段長度為1，僅擬合尾部數(shù)據(jù)段的多項式，計算平運動的相對變化量TR如式(2)[8]：

(2)

其中，Δnc為頭部數(shù)據(jù)n與外推值之間的差，Δnp為多項式外推至頭部數(shù)據(jù)歷元時n的變化量。方法如圖2所示。

圖2 基于平運動變化的異常TLE剔除方法Fig.2 The outlier filtering method based on mean motion

本文TR的閾值取為0.6，即認為TR大于0.6的TLE是異常數(shù)據(jù)。比較預處理前后天宮一號TLE序列的B*，結(jié)果如圖3所示?？梢钥闯觯A處理后B*序列的平滑性有了很大改善，認為預處理后天宮一號的TLE整體軌道精度有了很大改善。

圖3 預處理前后天宮一號的B*序列Fig.3 The pre and post preprocess of B* of TIANGONG-1

3 基于TLE的彈道系數(shù)計算方法

TLE數(shù)據(jù)并不能直接提供天宮一號的彈道系數(shù)信息，較為相關的是B*項，B*與彈道系數(shù)B之間的關系可以近似表示為[11]：B=12.741×B*，但是此結(jié)果與真實彈道系數(shù)仍存在較大誤差?；赥LE計算航天器的彈道系數(shù)，最簡單的辦法是根據(jù)兩條TLE的半長軸的衰減進行計算，即通過調(diào)整彈道系數(shù)，使較老TLE預報至較新TLE歷元時預報的半長軸與較新TLE的實際半長軸一致。彈道系數(shù)的修正可以表示為式(3)[12]：

(3)

其中，ΔaTLE為兩條TLE的半長軸實際衰減量，ΔaB為使用初始彈道系數(shù)B0進行軌道預報的半長軸衰減量，B′為修正后的彈道系數(shù)。這種方法的問題是TLE精度對彈道系數(shù)精度影響很大。本文對此方法進行改進，使用最小二乘方法對N(N≥3)條TLE進行彈道系數(shù)的擬合，使各TLE的實際半長軸與計算半長軸之間的誤差最小。基本過程如下：

(4)

其中，ρ為大氣密度，本文使用MSIS-90大氣模型計算，V為航天器相對于大氣的運動速度，e為偏心率，f為真近點角，ωE為地球自轉(zhuǎn)角速度，i為軌道傾角，μ為引力常數(shù)，M為平近點角。

4) 最后一條TLE計算結(jié)束后，使用公式(5)計算彈道系數(shù)修正值；

(5)

大氣阻力僅有長期作用效果，且采用開普勒平根數(shù)作為根數(shù)系統(tǒng)，因此進行軌道積分時的積分步長可取為軌道周期的整數(shù)倍，以提高計算效率。

4 姿態(tài)判別結(jié)果

彈道系數(shù)是空間目標的固有屬性，不同空間目標的彈道系數(shù)可能會有較大差異，同一空間目標在不同姿態(tài)下彈道系數(shù)也可能不同。

本文的彈道系數(shù)是基于軌道數(shù)據(jù)的解算結(jié)果，彈道系數(shù)吸收了軌道誤差、大氣環(huán)境建模誤差及航天器姿態(tài)變化等，因此解算的彈道系數(shù)并不是常值。對于天宮一號而言，預處理后的TLE用于解算彈道系數(shù)時精度已經(jīng)足夠，對彈道系數(shù)精度影響很??；大氣環(huán)境建模的誤差在大氣環(huán)境平靜時影響不超過10%，擾動期間影響可能超過20%[13]，是姿態(tài)穩(wěn)定時造成彈道系數(shù)波動的最主要因素；空間目標姿態(tài)穩(wěn)定時面質(zhì)比不發(fā)生變化，快速翻滾的航天器在翻滾方向不發(fā)生改變的情況下也可以認為姿態(tài)對彈道系數(shù)無影響。天宮一號外形結(jié)構(gòu)比較復雜，姿態(tài)變化引起的彈道系數(shù)變化較為明顯。使用前文所述方法對2015年1月1日起美國戰(zhàn)略司令部公開發(fā)布的TLE進行預處理及彈道系數(shù)計算，結(jié)果如圖4所示。

圖4 天宮一號的彈道系數(shù)結(jié)果Fig.4 The ballistic coefficient estimation results of TIANGONG-1

其中，綠色虛線為天宮一號終止數(shù)據(jù)服務時間，即2016年3月16日，紅色虛線前后(大約2017年10月8日)天宮一號的彈道系數(shù)有明顯差異。分段進行天宮一號彈道系數(shù)分析：

1) 天宮一號數(shù)據(jù)服務終止前后

從圖5可以看到，數(shù)據(jù)服務終止前彈道系數(shù)的均值約為0.0105 m2/kg，數(shù)據(jù)服務終止后至紅色虛線所示日期前彈道系數(shù)的均值約為0.0113 m2/kg，前后相差不到8%，差異主要由空間環(huán)境周期性變化引起。數(shù)據(jù)服務終止前后天宮一號的彈道系數(shù)變化特點基本一致，即使空間環(huán)境穩(wěn)定也存在周期性變化，變化周期約為5天，初步判斷天宮一號數(shù)據(jù)服務終止后姿態(tài)未發(fā)生明顯變化。

圖5 天宮一號數(shù)據(jù)服務終止前后彈道系數(shù)特性Fig.5 The characteristic of ballistic coefficients of TIANGONG-1 before and after data service termination

2) 2017年10月8日前后

2017年10月5日至10月13日天宮一號的彈道系數(shù)如圖6所示，期間天宮一號的彈道系數(shù)逐漸減小，10月6日前天宮一號的彈道系數(shù)均值約為0.0113 m2/kg，10月12日后天宮一號彈道系數(shù)的均值約為0.0069 m2/kg。期間太陽活動及地磁環(huán)境均穩(wěn)定，無磁擾、磁暴等異常情況發(fā)生，可以排除彈道系數(shù)變化由大氣環(huán)境建模誤差引起。分析認為天宮一號姿態(tài)變化引起彈道系數(shù)持續(xù)減小，2017年10月6日至10月12日，天宮一號姿態(tài)逐漸變化。

圖6 2017年10月5日至10月13日天宮一號的彈道系數(shù)Fig.6 The ballistic coefficients of TIANGONG-1 between Oct. 5 and Oct. 13, 2017

圖7 2017年10月12日后天宮一號的彈道系數(shù)結(jié)果Fig.7 The ballistic coefficients of TIANGONG-1 after Oct. 12, 2017

3) 2017年10月12日后

圖7為2017年10月12日后天宮一號的彈道系數(shù)結(jié)果及地磁指數(shù)情況?？梢钥闯?，彈道系數(shù)與地磁環(huán)境具有明顯相關性，可以認為地磁變化是造成彈道系數(shù)變化的主要原因，這是因為現(xiàn)有的大氣密度模型不夠精確，地磁環(huán)境擾動情況下的大氣密度計算誤差大于平靜時期的大氣密度計算誤差。期間天宮一號的翻滾速度可能繼續(xù)加快，但由于翻滾方向不變的情況下平均等效迎風面積不變，無法再基于彈道系數(shù)判斷翻滾速度變化情況。

綜上，可以將天宮一號的姿態(tài)變化分為以下幾個階段：

1) 穩(wěn)定階段：2016年3月16日至2017年10月6日。天宮一號的彈道系數(shù)變化特點與終止數(shù)據(jù)服務前基本一致，認為天宮一號處于穩(wěn)定狀態(tài)；

2) 逐漸變化階段：2017年10月6日至2017年10月12日。天宮一號的姿態(tài)從穩(wěn)定轉(zhuǎn)為緩慢翻滾，彈道系數(shù)逐漸減?。?/p>

3) 翻滾階段：2017年10月12日以后。天宮一號處于持續(xù)翻滾狀態(tài)，翻滾速度變化對彈道系數(shù)的影響很小，彈道系數(shù)基本穩(wěn)定。

5 結(jié)論

本文從姿態(tài)變化引起迎風面積變化這一特點出發(fā)，使用美國戰(zhàn)略司令部公開發(fā)布的雙行根數(shù)TLE作為數(shù)據(jù)源，提出基于彈道系數(shù)的天宮一號姿態(tài)判別方法，判別出天宮一號姿態(tài)穩(wěn)定、逐漸變化及翻滾三個階段。實測的雷達特征數(shù)據(jù)發(fā)現(xiàn)，天宮一號于2017年10月中旬開始有明顯翻滾現(xiàn)象，結(jié)合本文的分析結(jié)果，可以將天宮一號姿態(tài)變化的時間精確到10月6日至10月12日之間。本文提出的方法可以作為傳統(tǒng)姿態(tài)判別方法的有效補充。

值得指出的是，本文提出的方法也存在一定局限性。在飛行器結(jié)構(gòu)較為簡單、姿態(tài)失穩(wěn)后彈道系數(shù)變化不明顯的情況下，可能無法進行準確判斷。