基于預(yù)報(bào)偏差的LEO航天器軌道異常檢測

張濤濤 白顯宗 郝嘉 陳磊

（國防科學(xué)技術(shù)大學(xué)航天與材料工程學(xué)院，長沙410073）

1 引言

當(dāng)航天器受到軌道機(jī)動(dòng)、空間天氣突變、解體或者極端情況下的碰撞等因素影響時(shí)，航天器的軌道會(huì)發(fā)生突變，此時(shí)，軌道異常事件檢測信息可以為航天器所有者和操作者采取正確的措施提供依據(jù)。國內(nèi)外學(xué)者針對航天器軌道異常問題，提出了多種檢測方法。Patera提出了一種基于移動(dòng)窗口曲線擬合 （Moving Window Curve Fit，MWCF）技術(shù)的空間事件檢測方法，該方法在一定長度的時(shí)間窗口內(nèi)對數(shù)據(jù)進(jìn)行擬合，利用編目值和擬合值相減得到偏差，通過窗口的移動(dòng)得到偏差數(shù)據(jù)序列，最后通過統(tǒng)計(jì)偏差數(shù)據(jù)的標(biāo)準(zhǔn)差設(shè)定檢測門限進(jìn)行異常檢測［1］。Kelecy、Swartz等人基于Patera的曲線擬合技術(shù)對軌道異常檢測方法進(jìn)行了改進(jìn)［2-3］。基于曲線擬合的檢測方法存在移動(dòng)窗口的長度和擬合曲線的階次設(shè)定比較困難的不足。此外，董云峰等人提出了一種基于小波分析的軌道機(jī)動(dòng)檢測方法［4］。楊旭等人提出了單個(gè)衛(wèi)星相鄰根數(shù)時(shí)間差控制加綜合判據(jù)的衛(wèi)星軌道異常分析方法［5］。上述檢測方法基本上從數(shù)據(jù)處理入手，利用數(shù)學(xué)方法直接分析軌道數(shù)據(jù)，不利于分析航天器軌道變化的物理規(guī)律。本文將航天器軌道預(yù)報(bào)與數(shù)學(xué)檢測方法相結(jié)合進(jìn)行航天器軌道異常的分析。

北美防空司令部（NORAD）發(fā)布的兩行軌道要素（Two Line Elements，TLE）是目前最完整的空間目標(biāo)編目數(shù)據(jù)，其對應(yīng)的SGP4（Simplified General Perturbation 4）軌道預(yù)報(bào)模型是解析模型。作為例子，本文對航天器歷史TLE數(shù)據(jù)進(jìn)行分析，使用SGP4模型的長期項(xiàng)部分預(yù)報(bào)軌道根數(shù)的長期變化并基于馬氏距離對航天器軌道異常進(jìn)行判定。此外，本文所提出的方法同樣適用于其他類型的軌道數(shù)據(jù)以及相應(yīng)預(yù)報(bào)模型。

2 特征軌道要素長期變化預(yù)報(bào)模型

（1）特征軌道要素的選擇

特征軌道要素是一個(gè)或多個(gè)可以表征航天器軌道異常的軌道要素或?qū)С鰠?shù)，選取合適的特征軌道要素是進(jìn)行軌道異常檢測的基礎(chǔ)。航天器軌道異常主要可以分為兩類：一類是軌道平面內(nèi)軌道形狀和大小發(fā)生突變，主要表現(xiàn)為半長軸、偏心率等參數(shù)的變化；另一類是軌道面的空間位置發(fā)生突變，主要表現(xiàn)在傾角、升交點(diǎn)赤經(jīng)等參數(shù)的變化。為能夠充分地表征這兩類異常，本文選擇軌道半長軸a和軌道傾角i作為特征軌道要素。

（2）基于SGP4模型的軌道要素長期變化預(yù)報(bào)

SGP4是一種考慮了地球非球形攝動(dòng)、大氣阻力攝動(dòng)、第三體引力攝動(dòng)等因素綜合影響的軌道傳播計(jì)算模型，其可以配合NORAD發(fā)布的TLE平根數(shù)進(jìn)行航天器軌道預(yù)報(bào)［6-9］。

TLE數(shù)據(jù)包括空間目標(biāo)的軌道要素以及其他相關(guān)信息，如目標(biāo)編號(hào)、國際編號(hào)、大氣阻力系數(shù)等。其中用于SGP4模型預(yù)報(bào)的參數(shù)包括：

式中t0為歷元時(shí)刻；n0為歷元時(shí)刻的平均運(yùn)動(dòng)角速度 （轉(zhuǎn)／天）；e0為歷元時(shí)刻的平均軌道偏心率；i0為歷元時(shí)刻的平均軌道傾角；ω0為歷元時(shí)刻的平均近地點(diǎn)角距；Ω0為歷元時(shí)刻的平均升交點(diǎn)赤經(jīng)；M0為歷元時(shí)刻的平均平近點(diǎn)角；B＊為大氣阻力系數(shù)；這些軌道要素中的下標(biāo) “0”表示歷元時(shí)刻的值。由TLE數(shù)據(jù)直接得到的n0是Kozai格式，需要轉(zhuǎn)化為Brouwer格式平均角速度才可以使用SGP4模型進(jìn)行軌道預(yù)報(bào)［9］：

對于LEO航天器，地球非球形攝動(dòng)J2項(xiàng)對于半長軸的長期影響為零，同時(shí)由于其軌道高度較低，第三體引力、太陽光壓等對半長軸的長期攝動(dòng)可以忽略。影響LEO航天器半長軸長期變化的主要因素是大氣阻力。盡管在這樣的高度大氣已經(jīng)極為稀薄，但是由于航天器的飛行速度大、持續(xù)時(shí)間長，并且大氣阻力不是保守力而是耗散力，因此航天器半長軸將緩慢減小。在大氣阻力攝動(dòng)影響下，半長軸長期變化可由下式描述［9］：

式中t-t0為預(yù)報(bào)時(shí)間間隔；C1，D2，D3，D4為相關(guān)的系數(shù)，詳見參考文獻(xiàn)［9］。

對于LEO航天器所處軌道中的非太陽同步軌道，其傾角受地球非球形攝動(dòng)、大氣阻力攝動(dòng)、第三體引力攝動(dòng)影響產(chǎn)生的長期變化量很小，可以近似為零。對于太陽同步軌道上的航天器，由于其軌道升交點(diǎn)赤經(jīng)的進(jìn)動(dòng)角速度同地球公轉(zhuǎn)角速度相同，太陽與軌道面的相對位置不變，因此太陽引力會(huì)導(dǎo)致軌道傾角的長期變化，可以表示為［10］

式中ns為地球繞太陽運(yùn)動(dòng)的平均角速度；n為航天器運(yùn)動(dòng)角速度；is為黃赤夾角；us為太陽在赤道上的平經(jīng)度；Ω為航天器升交點(diǎn)赤經(jīng)。

3 偏差數(shù)據(jù)生成

（1）軌道預(yù)報(bào)

軌道預(yù)報(bào)的基本方法如下：利用tk-1時(shí)刻的TLE編目軌道要素σk-1，通過SGP4模型長期項(xiàng)預(yù)報(bào)得到tk時(shí)刻軌道要素長期變化的預(yù)期值σk／（k-1），同時(shí)利用tk時(shí)刻的編目根數(shù)σk得到ak和ik，如圖1。其中，下標(biāo)“k／（k-1）”表示由tk-1歷元時(shí)刻預(yù)報(bào)得到的tk歷元時(shí)刻的預(yù)報(bào)值，下標(biāo)“k”表示tk歷元時(shí)刻的編目值。

圖1 基于SGP4長期項(xiàng)模型的軌道預(yù)報(bào)示意Fig.1 Orbit prediction based on SGP4model secular term

（2）數(shù)據(jù)平滑及偏差生成

TLE作為一種有效的空間目標(biāo)編目數(shù)據(jù)，被廣泛用于空間目標(biāo)監(jiān)視領(lǐng)域，例如航天器機(jī)動(dòng)檢測、碰撞交會(huì)分析、校核大氣密度模型等［3］。由于測量、定軌以及模型精度 （例如地球模型、大氣模型、預(yù)報(bào)模型）等因素影響，在得到的編目和預(yù)報(bào)數(shù)據(jù)中包含了一定的誤差，上述誤差會(huì)造成編目值和預(yù)報(bào)值在真值周圍抖動(dòng)，如圖2。

圖2 特征參數(shù)的抖動(dòng)Fig.2 Dithering of characteristic parameters

這種抖動(dòng)可能導(dǎo)致正常情況下偏差數(shù)據(jù)增大從而造成虛警，為了消除該小抖動(dòng)并盡可能的保持原特征參數(shù)的變化，使用三次樣條曲線對數(shù)據(jù)進(jìn)行平滑處理。建立擬合目標(biāo)函數(shù)如式（2），其中第一項(xiàng)代表擬合偏差，第二項(xiàng)代表擬合的平滑程度。

式中X為待平滑處理的數(shù)據(jù)序列；f為平滑得到的函數(shù)；L為所有歷元時(shí)刻間隔中最小的時(shí)間段。

ω為決定數(shù)據(jù)擬合效果的權(quán)系數(shù)，取值范圍為［0，1］。ω越大，擬合誤差越小，但是擬合平滑效果越差；反之，ω越小，擬合得到的曲線越光滑，但是擬合偏差越大。為了在保持較小擬合誤差的前提下獲取較好的擬合效果，可以根據(jù)軌道情況的不同變化選取合適的ω，一般情況下，對于瞬時(shí)變化比較劇烈的情況，應(yīng)該首要考慮控制擬合偏差，即選擇較大的ω對于軌道參數(shù)變化較平緩的情況，可以適當(dāng)?shù)倪x擇較小的ω，以便獲得較好的擬合效果。本文中對二者進(jìn)行了綜合考慮，利用以下經(jīng)驗(yàn)公式確定ω：

式中tj為第j個(gè)時(shí)間間隔；為所有歷元時(shí)刻間隔的平均值。

依據(jù)目標(biāo)函數(shù)，在Matlab中使用三次樣條曲線對軌道半長軸和傾角的編目值和預(yù)報(bào)值進(jìn)行平滑處理可得偏差數(shù)據(jù)序列：

式中

4 異常檢測

偏差數(shù)據(jù)中包含半長軸和傾角兩個(gè)參數(shù)，此時(shí)的檢測可視為多維數(shù)據(jù)異常檢測。在多維異常檢測中，當(dāng)各維分量代表不同質(zhì)的量或者各分量的差異很大時(shí) （如本方法中兩個(gè)特征參數(shù)偏差量級相差較大），利用普通的歐氏距離計(jì)算會(huì)出現(xiàn) “大數(shù)吃小數(shù)”的情況。馬氏距離是對各維數(shù)據(jù)經(jīng)過加權(quán)處理后描述兩個(gè)數(shù)據(jù)點(diǎn)之間接近程度的量，可以有效避免該問題的發(fā)生。設(shè)總體均值為μ，協(xié)方差陣為Σ，ri為一個(gè)樣本矢量，則樣本ri到總體的馬氏距離為［11］

本文進(jìn)行偏差數(shù)據(jù)異常分析時(shí)，總體的均值μ與協(xié)方差矩陣Σ均是未知的，可用它們的無偏估計(jì)即樣本均值和樣本協(xié)方差陣C代替，此時(shí)數(shù)據(jù)點(diǎn)ri到總體的馬氏距離為

馬氏距離實(shí)質(zhì)上描述了考慮概率分布的數(shù)據(jù)點(diǎn)與分布中心偏離程度，樣本值到總體的馬氏距離越大表示產(chǎn)生此樣本的概率越小，反之，馬氏距離越小表示產(chǎn)生此樣本的概率越大。

異常檢測的偏差門限值dT可以依據(jù)檢測精度選取，dT越大則誤警率越低，但是漏警率越高；dT越小則漏警率越低，但是誤警率越高。一般情況下可取為3～5。如果某樣本值到總體的馬氏距離大于給定的門限值dT，則認(rèn)為該樣本值為異常值。綜上可得具體的算法流程如圖3所示。

圖3 基于預(yù)報(bào)偏差的軌道異常檢測流程Fig.3 Flow chart of the detection method based on prediction dispersion

5 算例

Terra衛(wèi)星是一顆美國的太陽同步軌道衛(wèi)星，其裝載的5種傳感器能同時(shí)采集地球大氣、陸地、海洋和太陽能量平衡的信息從而實(shí)現(xiàn)對地球的觀測。下面以其為例，用本文所提出的方法進(jìn)行軌道異常檢測分析。表1為已知的Terra衛(wèi)星2010年的軌道機(jī)動(dòng)歷史數(shù)據(jù)［12］。

表1 Terra衛(wèi)星2010年軌道機(jī)動(dòng)事件Tab.1 Maneuver history of Terra in 2010

對Terra衛(wèi)星，利用NORAD發(fā)布的2010年的TLE數(shù)據(jù)［13］，進(jìn)行分析可得到其2010年半長軸和傾角的變化情況如圖4，圖中豎線表示Terra衛(wèi)星2010年的機(jī)動(dòng)歷史。

圖4 Terra衛(wèi)星2010年半長軸、傾角變化以及機(jī)動(dòng)歷史Fig.4 Semi-major axis，inclination history and recorded maneuvers of Terra in 2010

通過本文第3節(jié)所建立的軌道預(yù)報(bào)模型對Terra衛(wèi)星的2010年軌道數(shù)據(jù)進(jìn)行處理，可以得到其特征軌道要素的預(yù)報(bào)值，使用數(shù)據(jù)平滑處理方法對編目值和預(yù)報(bào)值進(jìn)行處理可以得到平滑值。為便于說明問題，下面僅就Terra在2010年的第222天至228天之間的半長軸數(shù)據(jù)處理情況進(jìn)行詳細(xì)說明（見圖5），軌道傾角的處理與此類似。

對Terra衛(wèi)星的半長軸和軌道傾角進(jìn)行預(yù)報(bào)、平滑處理后根據(jù)公式（3）可以得到其偏差數(shù)據(jù)，利用公式（4）可以得到每個(gè)歷元時(shí)刻的馬氏距離，并選取檢測門限值為5對馬氏距離數(shù)據(jù)進(jìn)行異常分析，檢測結(jié)果如圖6和圖7，圖6中的①～⑥是檢測到的6次軌道異常。

由計(jì)算結(jié)果可以看出，Terra衛(wèi)星在2010年的異常值集中在6個(gè)時(shí)間點(diǎn)附近。分析原因可知：

1）軌道維持機(jī)動(dòng)往往不是一次調(diào)整完成，而是在一段時(shí)間內(nèi)不斷的調(diào)整和測量達(dá)到設(shè)定值。

2）TLE定軌數(shù)據(jù)受到歷史數(shù)據(jù)的影響，當(dāng)發(fā)生軌道異常時(shí)不可能一次完成精確的定軌，而是利用多組新的軌道數(shù)據(jù)才能進(jìn)行穩(wěn)定的定軌。

綜上可以判斷Terra衛(wèi)星在2010發(fā)生了6次軌道異常。檢測結(jié)果的次數(shù)和時(shí)間均與已知的軌道機(jī)動(dòng)歷史相符。此外通過分析可知第1、3、4次為軌道高度維持，第2、5、6次為軌道傾角機(jī)動(dòng)，而且其在傾角機(jī)動(dòng)的同時(shí)對軌道高度進(jìn)行了提升。

圖6 Terra衛(wèi)星2010年機(jī)動(dòng)檢測結(jié)果Fig.6 Maneuver detection results of Terra in 2010

圖7 Terra衛(wèi)星2010年半長軸、傾角偏差分布圖Fig.7 Semi-major axis dispersion and inclination dispersion distributing of Terra in 2010

6 結(jié)束語

本文研究軌道異常檢測問題，通過TLE數(shù)據(jù)分析了特征軌道要素的變化情況，在此基礎(chǔ)上，利用軌道預(yù)報(bào)模型得到特征軌道要素的預(yù)報(bào)值，通過三次樣條曲線擬合方法對編目值和預(yù)報(bào)值中的小擾動(dòng)進(jìn)行平滑處理后得到偏差數(shù)據(jù)，基于馬氏距離對多維偏差數(shù)據(jù)進(jìn)行異常檢測，最后以Terra衛(wèi)星2010年的機(jī)動(dòng)事件檢測為例進(jìn)行了分析。

在以上研究的基礎(chǔ)上，得出了如下結(jié)論：

1）本文所提出的方法能夠確定軌道異常的次數(shù)、時(shí)間，選擇軌道半長軸和軌道傾角作為特征軌道要素能夠確定軌道的異常類型 （共面或異面）。

2）通過對TLE編目數(shù)據(jù)和預(yù)報(bào)數(shù)據(jù)進(jìn)行平滑處理，可以有效地抑制了各種小擾動(dòng)造成虛警。但是由于TLE數(shù)據(jù)本身精度和SGP4模型的預(yù)報(bào)精度有限，因此其不適合分析量級變化較小的機(jī)動(dòng)，對于變化量級更小的軌道異常事件，使用更加精密的軌道預(yù)報(bào)模型和相應(yīng)的數(shù)據(jù)進(jìn)行檢測將是下一步工作的重點(diǎn)。

［1］PATERA R P.Space event detection method ［J］.Journal of Spacecraft and Rockets，2008，45（3）：554-559.

［2］KELECY T，HALL D，HAMADA K，et al.Satellite maneuver detection using Two-line Element（TLE）data ［C］.Proceedings of the Advanced Maui Optical and Space Surveillance Technologies Conference，Maui，Hawaii，September，2007.

［3］SWARTZ R L JR，COGGI J，et al.A Swift SIFT for Satellite Event Detection ［C］.AIAA／AAS Astrodynamic Specialist Conference，Toronto，Ontario Canada，August 2-5，2010.

［4］董云峰，蘇建敏.利用小波分析識(shí)別空間目標(biāo)的軌道機(jī)動(dòng) ［J］.宇航學(xué)報(bào)，2004，25（2）：213-218.DONG YUNFENG，SU JIANMIN.Detection of space target orbit maneuver by wavelet analysis ［J］.Journal of Astronautics，2004，25（2）：213-218.

［5］楊旭，劉靜，吳相彬，等.利用TLE數(shù)據(jù)分析LEO衛(wèi)星軌道異常的新方法——綜合判據(jù)法 ［J］.空間科學(xué)學(xué)報(bào)，2011，31（2）：223-228.YANG XU，LIU JING，WU XIANGBIN，et al.New method to analyze the orbital abnormal of LEO satellite using TLE data——compositive criterion［J］.Chinese Journal of Space Science，2011，31（2）：223-228.

［6］HOOTS F R，ROEHRICH R L.Space Track Report No.3：Models for Propagation of NORAD Element sets［R］.Peterson：Aerospace Defense Command，United States Air Force，1980.

［7］VALLADO D A，CRAWFORD P，HUJSAK R，et al.Revisiting Spacetrack Report＃3 ［C］.AIAA／AAS Astrodynamics Specialist Conference and Exhibit，Keystone，Colorado，August 21-24，2006.

［8］VALLADO D A，CRAWFORD P.SGP4Orbit Determination ［C］.AIAA／AAS Astrodynamics Specialist Conference and Exhibit，Honolulu，Hawaii，August 18-21，2008.

［9］HOOTS F R，SCHUMACHER R W，GLOVER R A.History of analytical orbit modeling in the U.S.space surveillance system ［J］.Journal of Guidance，Control，and Dynamics，2004，27（2）：174-185.

［10］楊嘉墀.航天器軌道動(dòng)力學(xué)與控制.［M］.北京：宇航出版社，2001.

［11］吳翊，李永樂，胡慶軍.應(yīng)用數(shù)理統(tǒng)計(jì) ［M］.長沙：國防科技大學(xué)出版社，1995.

［12］UNITED STATES STRATEGIC COMMAND（USSTRATCOM）.Space-Track-Catalog Number Query Results.［EB／OL］.［2012-04-16］.http：∥www.space-track.org／perl／id＿query.pl.

［13］ATMOSPHERIC SCIENCE DATA CENTER.Terra Spacecraft Loss of Accurate Orbit Data Record［EB／OL］.［2012-04-16］.http：∥eosweb.larc.nasa.gov∥PRODOCS／misr／terra＿maneuvers.html.