利用TLE數(shù)據(jù)判別天宮一號(hào)目標(biāo)飛行器姿態(tài)

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

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

1 引言

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

天宮一號(hào)的再入過程有以下幾個(gè)特點(diǎn)：

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

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

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

因此，天宮一號(hào)的再入案例對(duì)于開展軌道動(dòng)力學(xué)和再入預(yù)報(bào)研究具有特殊意義。

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

(1)

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

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

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

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

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

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

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

3) 軌道機(jī)動(dòng)檢測(cè)；

4) 基于平運(yùn)動(dòng)變化剔除異常TLE。

本節(jié)主要介紹軌道機(jī)動(dòng)檢測(cè)和基于平運(yùn)動(dòng)變化剔除異常TLE的方法。

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

圖1 軌道機(jī)動(dòng)檢測(cè)的基本原理Fig.1 The basic principle of maneuver detection

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

(2)

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

圖2 基于平運(yùn)動(dòng)變化的異常TLE剔除方法Fig.2 The outlier filtering method based on mean motion

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

圖3 預(yù)處理前后天宮一號(hào)的B*序列Fig.3 The pre and post preprocess of B* of TIANGONG-1

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

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

(3)

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

(4)

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

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

(5)

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

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

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

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

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

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

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

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

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

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

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

3) 2017年10月12日后

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

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

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

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

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

5 結(jié)論

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

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