李 濤,黃 昊,陳 磊

(1. 國防科技大學 空天科學學院, 湖南 長沙 410073; 2. 軍事科學院 國防科技創(chuàng)新研究院, 北京 100071)

在當前空間目標數(shù)量日益增長且多為非合作目標的背景下，及時發(fā)現(xiàn)航天器軌道機動、碰撞解體等空間事件具有重要意義。一方面，及時獲知空間事件信息對及時采取合適的應對措施至關(guān)重要；另一方面，基于歷史機動信息可以分析預測航天器的活動規(guī)律和機動行為，有效支撐空間態(tài)勢感知。

軌道機動會引起空間目標的軌道異常，檢測特定軌道參數(shù)的異常變化是從歷史軌道數(shù)據(jù)中獲取機動信息的有效途徑。美國戰(zhàn)略司令部發(fā)布的兩行軌道要素(Two Line Elements, TLE)數(shù)據(jù)是目前為止公開發(fā)布的最完備的空間目標軌道信息編目數(shù)據(jù)。國內(nèi)外學者已提出一系列從航天器歷史TLE數(shù)據(jù)中檢測軌道機動的方法。概括而言，這些方法分為兩類。一類直接對編目軌道要素或其導出量的歷史數(shù)據(jù)進行統(tǒng)計分析，以檢出統(tǒng)計意義上的異常編目數(shù)據(jù)[1-4]；另一類通過比較軌道要素編目值與相應預報值判斷編目值是否異常[5-7]，其中預報值通常通過與TLE數(shù)據(jù)配套使用的SGP4/SDP4軌道預報模型[8]得到，編目值與預報值的差，即預報誤差則被作為異常值檢測的指標。由于TLE本身包含誤差且軌道預報模型也會引入一部分誤差，因此預報誤差值具有不確定性且其分布特性與預報時間有關(guān)?，F(xiàn)有的機動檢測方法并未考慮其概率分布特性，僅利用預報誤差的統(tǒng)計均值和方差信息確定異常值檢測的門限。方法中往往暗含預報誤差服從高斯分布的假設(shè)，這一未經(jīng)驗證的預設(shè)條件可能會對預報誤差異常值檢測的造成干擾，進而影響機動檢測的準確率。

本文提出一種基于預報誤差分布擬合的機動檢測方法。通過分析預報誤差檢測軌道參數(shù)的異常編目值，但不預設(shè)預報誤差的概率分布形式，而是利用歷史數(shù)據(jù)進行擬合。對目標的歷史TLE數(shù)據(jù)兩兩預報求差獲得預報誤差的樣本數(shù)據(jù)后，利用期望最大化(Expectation Maximization, EM)算法對樣本數(shù)據(jù)擬合得到以高斯混合模型表達的預報誤差的概率分布模型。隨后，以此為基礎(chǔ)設(shè)計軌道參數(shù)異常值檢測方法，并利用軌道機動與軌道異常之間的對應關(guān)系實現(xiàn)機動檢測。最后，選取典型活動衛(wèi)星為測試目標進行機動檢測，分析驗證所提方法的可行性和機動檢測性能。需要指出的是，在所設(shè)計的軌道機動檢測方法中，用于分析的軌道參數(shù)的選擇并不是唯一的，本文僅以衛(wèi)星的軌道半長軸為例對機動檢測方法進行說明，實際上也可以選擇軌道傾角、偏心率、軌道能量等其他參數(shù)。

1 基于歷史數(shù)據(jù)的預報誤差分布擬合

1.1 樣本數(shù)據(jù)的生成

預報誤差樣本數(shù)據(jù)的生成方法如圖1所示。利用SGP4/SDP4模型，ti時刻的TLE編目軌道參數(shù)ai將被順次外推至隨后m個編目值所在的歷元時刻ti+1，ti+2，…，ti+m，以得到相應的預報值ai+1,i，ai+2,i，…，ai+m,i。將這些預報值與相應的編目值相減，就得到m個預報誤差數(shù)據(jù)Δai+1,i，Δai+2,i，…，Δai+m,i，以及對應的預報時間Δti+1,i=ti+1-ti，Δti+2,i=ti+2-ti，…，Δti+m,i=ti+m-ti。

圖1 生成預報誤差樣本數(shù)據(jù)的示意圖Fig.1 Sample data generation of the prediction error

對所有容許的軌道參數(shù)編目值順次重復上述過程，可以得到一系列與不同預報時間相對應的預報誤差數(shù)據(jù)。它們構(gòu)成了預報誤差的樣本數(shù)據(jù)，記為集合S。在生成樣本數(shù)據(jù)時，可以根據(jù)實際情況確定參數(shù)m的值以獲得合適的樣本數(shù)據(jù)量。

由于預報誤差的分布特性與預報時間密切相關(guān)，樣本數(shù)據(jù)被用于分布擬合之前需要按其預報時間進行分組?？紤]到同一空間目標的絕大部分TLE之間的發(fā)布時間間隔一般為軌道周期的整數(shù)倍[9]，集合S中各預報誤差對應的預報時間一般也近似等于軌道周期的整數(shù)倍，故各預報時間均可以與其軌道周期T的整數(shù)比值表示，再將具有相同預報時間的誤差數(shù)據(jù)進行分組。假設(shè)樣本數(shù)據(jù)最終被分為N個數(shù)據(jù)組，即

S={Sr1,…,Srk,…,SrN}

(1)

則Srk包含所有預報時間為rk倍軌道周期的預報誤差數(shù)據(jù)。

1.2 預報誤差的分布擬合

對預報誤差的分布擬合包括兩步：第一步是建立以高斯混合模型描述的預報誤差的概率分布數(shù)學模型；第二步是利用EM算法從樣本數(shù)據(jù)中學習概率分布模型的模型參數(shù)。

高斯混合模型常被用作描述隨機變量概率分布的參數(shù)模型，它由有限個高斯分布的加權(quán)和構(gòu)成。理論上來說，只要高斯分布的個數(shù)足夠多，利用高斯混合模型可以對任意形式的概率分布進行數(shù)學描述?；诟咚够旌夏Ｐ停A報時間為rk倍軌道周期時，預報誤差Δa的概率密度函數(shù)可表示如下

(2)

(3)

在式(2)表達的概率密度函數(shù)中，待確定的參數(shù)包括高斯分布的個數(shù)J以及各高斯分布對應的參數(shù)αj、μj和σj。文獻[10]利用高斯混合模型對TLE位置預報誤差的概率分布進行了擬合，結(jié)果表明高斯分布的個數(shù)取3足以得到滿意的擬合效果。以此為參照，這里直接取J=3，其他參數(shù)則利用EM算法從預報誤差樣本數(shù)據(jù)中學習得到。

EM算法是一種迭代算法，1977年由Dempster等總結(jié)提出[12]，用于含有隱變量的概率模型參數(shù)的極大似然估計，或極大后驗概率估計。EM算法的每次迭代由兩步組成：E步，求期望；M步，求極大。利用EM算法從預報誤差數(shù)據(jù)組Srk中估計frk(Δa)中參數(shù)的迭代計算過程[12]如下。

輸入：數(shù)據(jù)組Srk={Δa1,Δa2,…,ΔaK}及高斯混合模型frk(Δa)。

輸出：高斯混合模型的描述參數(shù)αj、μj和σj。

1)取參數(shù)的初始值開始迭代。

2)E步：依據(jù)當前模型參數(shù)，計算第j個高斯成分對觀測數(shù)據(jù)Δak的響應度

(4)

3)M步：計算新一輪迭代的模型參數(shù)

(5)

(6)

(7)

4)重復第2和第3步，直到收斂。

對數(shù)據(jù)集S中的各數(shù)據(jù)組依次重復上述擬合過程，即可得到預報時間為r1，r2，…，rN時預報誤差的概率密度函數(shù)，分別記為fr1(Δa)，fr2(Δa)，…，frN(Δa)?？紤]到樣本數(shù)據(jù)過小時會影響誤差分布的擬合效果，在上述擬合過程中忽略了數(shù)據(jù)量不足的數(shù)據(jù)組。

2 機動檢測方法

2.1 軌道參數(shù)異常值的檢測

通過分析預報誤差對軌道參數(shù)異常值進行檢測是一種行之有效的方法，檢測問題可概述為：根據(jù)將軌道參數(shù)序列的第i個編目值預報至第j個(i

根據(jù)1.1節(jié)所述樣本數(shù)據(jù)的生成方法，預報誤Δaj,i必定屬于樣本數(shù)據(jù)集S中的某數(shù)據(jù)組。假定第i個和第j個軌道參數(shù)編目值的歷元時間差滿足Δtj,i=rk，即等于rk倍軌道周期，則Δaj,i是數(shù)據(jù)組Srk中的元素。因此，問題變成判斷Δaj,i在數(shù)據(jù)組Srk中是否屬于異常值。

但是，考慮到預報誤差的真實分布，即擬合得到的概率分布函數(shù)frk(Δa)并不一定滿足高斯分布，因此直接應用n-σ法則進行異常值檢測并不合適。這里參照n-σ法則提出一種改進的異常值檢測法則，并命名為n-P法則，如圖2所示。

圖 2 n-P法則示意圖Fig. 2 Diagram of the n-P rule

n-P法則將超出區(qū)間[lbn,ubn]的數(shù)據(jù)視作異常值，類似地，n=1、2、3依次對應P1=0.682 7、P2=0.954 5、P3=0.997 3。區(qū)間[lbn,ubn]則根據(jù)概率Pn及擬合得到的概率密度函數(shù)frk(Δa)計算得到。記與frk(Δa)對應的累積概率分布函數(shù)為Frk，則lbn和ubn的計算公式如下

(8)

如果Δaj,i位于區(qū)間[lbn,ubn]之內(nèi)，則認為第j個軌道參數(shù)編目值是正常的，否則它將被判定為異常值。

考慮到由任意兩個滿足要求(i

2.2 軌道機動的檢測

軌道機動會使軌道參數(shù)發(fā)生異常變化，表現(xiàn)為軌道參數(shù)編目值序列中的異常值。理論上來說，檢出軌道參數(shù)編目值序列中的異常值，即可實現(xiàn)機動檢測。然而，由于TLE本身包含數(shù)據(jù)噪聲(這一噪聲來源于軌道確定過程中用來生成TLE數(shù)據(jù)的軌道狀態(tài)觀測值誤差以及建模過程中的不確定性)，即使某個編目值被檢出為異常值，并不能認定它必然是由軌道機動造成的。

在利用預報誤差進行機動檢測的背景下，為消除噪聲干擾，一個有效的方法是將軌道參數(shù)編目值分別預報至后續(xù)多個編目值所在歷元。然后根據(jù)所得的預報誤差判斷后續(xù)各編目值是否異常，只有隨后的大部分編目值都被判定為異常時才認為這是由軌道機動造成的。

基于該思想，從衛(wèi)星的軌道參數(shù)編目值序列中檢測歷史軌道機動的數(shù)據(jù)處理流程確定如下。

步驟1：依次選取編目值序列中的各元素為預報基準。對于每一基準數(shù)據(jù)，將其順次外推至后續(xù)m個編目值所在歷元(這里的m值與1.1節(jié)中的m值相同)，根據(jù)所得的預報誤差分別判斷這m個編目值是否異常，并記錄異常值個數(shù)。最后可以得到一個序列ID，序列ID中第i個元素記錄了以軌道參數(shù)編目值序列的第i個元素為預報基準時，后續(xù)m個編目值中被檢出為異常值的個數(shù)。

步驟2：對于編目值序列中與噪聲對應的單個異常值，在第1步中分別以在它之前的m-1個編目值為預報基準進行異常檢測時，序列ID中對應的元素值均為1。當以該異常值本身為預報基準時，序列ID中對應的元素值則為m。因此，噪聲引起的單個異常值一般對應于序列ID中長度為m+1、均值約為2的非零元素數(shù)據(jù)段。為消除TLE數(shù)據(jù)噪聲的干擾，從序列ID中找出由非零元素構(gòu)成且元素均值不超過3的數(shù)據(jù)段，并將這些數(shù)據(jù)段中的元素值重置為0，最終得到一個新序列，記為IDm。

步驟3：IDm中的非零元素可以認為必定是由軌道機動造成的。從該序列中提取出所有由非零元素組成的數(shù)據(jù)段，每一個數(shù)據(jù)段被認為代表一個軌道機動，并且數(shù)據(jù)段中最大元素值對應的預報基準所在歷元時刻被視為機動時間。

3 算例分析

3.1 算例設(shè)計

為評估所提機動檢測方法的檢測性能，選取7個典型的LEO軌道衛(wèi)星作為測試目標，如表1所示。

表1 選定衛(wèi)星概況Tab.1 Summary of selected satellites

對每一個衛(wèi)星，首先將軌道參數(shù)編目值序列中所有元素分別外推至其后續(xù)15個編目值所在歷元，即取m=15，以獲得預報誤差的樣本數(shù)據(jù)。然后，利用這些樣本數(shù)據(jù)擬合得到不同預報時間下預報誤差的概率分布模型。然后，建立基于擬合分布模型和2-P法則建立軌道參數(shù)異常值的檢測門限。最后，應用2.2節(jié)所述的檢測方法對軌道參數(shù)編目值序列進行分析，檢出軌道機動事件。

此外，為了驗證軌道參數(shù)預報誤差分布的非高斯特性，并評估該特性對機動檢測性能的影響，還設(shè)計了假設(shè)預報誤差服從高斯分布的對比試驗。與本文所提方法相比，對比試驗的唯一不同之處在于：在得到各衛(wèi)星的軌道參數(shù)預報誤差數(shù)據(jù)并按預報時間完成分組后，將采用高斯模型而不是高斯混合模型對預報誤差的概率分布進行擬合。對比試驗得到的檢測結(jié)果將被用于與本文所提方法得到的檢測結(jié)果之間進行對比，從而得出結(jié)論。

3.2 結(jié)果分析

以27424號衛(wèi)星為例對機動檢測的詳細過程說明如下。利用2003年1月1日至2016年12月31日的歷史TLE數(shù)據(jù)，共得到128925個軌道參數(shù)預報誤差數(shù)據(jù)。這些數(shù)據(jù)按照其對應的預報時間被分成196個數(shù)據(jù)組。圖3給出了典型數(shù)據(jù)組的柱狀圖。

(a) Δti+m,i=88T (b) Δti+m,i=102T

(c) Δti+m,i=116T (d) Δti+m,i=131T圖3 27424號衛(wèi)星典型預報誤差數(shù)據(jù)組的柱狀圖Fig.3 Histograms of typical prediction error data groups for object 27424

隨后，這些預報誤差數(shù)據(jù)組被用于擬合得到分別由高斯混合模型及高斯模型表示的概率分布模型，它們構(gòu)成了軌道參數(shù)異常值檢測的基礎(chǔ)。對前述數(shù)據(jù)組擬合得到的概率密度函數(shù)如圖4所示。

(a) Δti+m,i=88T (b) Δti+m,i=102T

(c) Δti+m,i=116T (d) Δti+m,i=131T圖4 對27424號衛(wèi)星典型預報誤差數(shù)據(jù)組進行擬合得到的概率密度函數(shù)Fig.4 Fitted probability density functions of typical prediction error data groups for object 27424

與擬合得到的預測誤差的兩種概率密度函數(shù)相對應，可以得到兩種軌道參數(shù)異常值的檢測門限。相應地，利用2.2節(jié)所述數(shù)據(jù)處理方法對27424號衛(wèi)星的軌道參數(shù)編目值序列進行處理，得到兩種機動檢測結(jié)果(如圖5和圖6所示，其中檢出的軌道機動以“*”標注)。對其余6顆衛(wèi)星重復上述機動檢測流程，即可得到所有衛(wèi)星的機動檢測結(jié)果(見表2)。

由表2可知，利用人工方法檢測出7個衛(wèi)星共進行了520次機動。而本文所提方法成功檢出了其中的478次，成功率達到了92%左右，在漏檢42次機動的同時發(fā)生了7次誤檢?？偟膩碚f，該方法的機動檢測準確率較好且能有效避免誤檢。

相較而言，對比試驗只檢測出了520次機動中的439次，成功率大致為84%，并且誤檢次數(shù)和漏檢次數(shù)均達到本文所提方法的兩倍左右。在對比試驗中，實際上隱含了預報誤差服從高斯分布的假設(shè)。相應地，軌道參數(shù)異常值的檢測門限基于2-σ法則確定。圖4及圖5表明，預報誤差的真實分布實際上與高斯分布相去甚遠，可知假定預報誤差服從高斯分布并不合適。在高斯分布假設(shè)下，盡管大部分量級較大的機動事件能被正確檢出，但對于大部分小量級機動事件及部分大量級機動事件，這一假設(shè)將嚴重限制算法的檢測性能，導致較大概率的誤檢和漏檢。反之，基于預報誤差分布擬合的機動檢測方法則能取得較好的檢測效果，出現(xiàn)誤檢和漏檢的概率也更低。

(a) 衛(wèi)星的半長軸時間-歷史(a) Semi-major axis time-history

(b) 檢測到的異常值個數(shù) (b) Detected number of anomaly points圖5 基于本文所述機動檢測方法的機動檢測結(jié)果Fig.5 Maneuver detection result of our method

(a) 衛(wèi)星的半長軸時間-歷史(a) Semi-major axis time-history

(b) 檢測到的異常值個數(shù)(b) Detected number of anomaly points圖6 對比試驗得出的機動檢測結(jié)果Fig.6 Maneuver detection result of contrast test

表2所有衛(wèi)星的機動檢測結(jié)果
Tab.2 Maneuver detection results of all satellites

衛(wèi)星編號(SSN)機動次數(shù)(人工檢測)機動次數(shù)(本文方法)機動次數(shù)(對比試驗)誤檢漏檢成功誤檢漏檢成功259946402620135126620452243611342742412301211121710528376104011930208429499362135413529107782126621464291087012681466合計5207424781581439

4 結(jié)論

本文提出了一種基于預報誤差分布擬合的軌道機動檢測方法。該方法首先通過對衛(wèi)星軌道參數(shù)編目值和預報值進行比較并分析預報誤差來識別軌道參數(shù)異常值，進而得到相應的歷史機動信息。在利用預報誤差進行軌道參數(shù)異常值檢測時，并不預設(shè)預報誤差的概率分布形式，而是利用歷史數(shù)據(jù)擬合得到。隨后，基于擬合得到的分布模型推導得到對軌道參數(shù)異常值的檢測方法。最后，根據(jù)軌道機動和軌道參數(shù)異常編目值的對應關(guān)系，提出從軌道參數(shù)編目值序列中獲得歷史機動的具體方法。對典型衛(wèi)星的機動檢測結(jié)果表明，假定軌道參數(shù)預報誤差服從高斯分布并不符合實際情況，而利用高斯混合模型則可以對預報誤差的概率分布進行很好的擬合。以之為基礎(chǔ)構(gòu)建的機動檢測方法能夠在準確檢出絕大部分歷史機動的同時有效避免誤檢事件的發(fā)生。