李 濤,黃 昊,陳 磊

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

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

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

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

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

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

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

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

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

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

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

(1)

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

1.2 預(yù)報(bào)誤差的分布擬合

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

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

(2)

(3)

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

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

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

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

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

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

(4)

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

(5)

(6)

(7)

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

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

2 機(jī)動檢測方法

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

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

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

但是，考慮到預(yù)報(bào)誤差的真實(shí)分布，即擬合得到的概率分布函數(shù)frk(Δa)并不一定滿足高斯分布，因此直接應(yīng)用n-σ法則進(jìn)行異常值檢測并不合適。這里參照n-σ法則提出一種改進(jì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依次對應(yīng)P1=0.682 7、P2=0.954 5、P3=0.997 3。區(qū)間[lbn,ubn]則根據(jù)概率Pn及擬合得到的概率密度函數(shù)frk(Δa)計(jì)算得到。記與frk(Δa)對應(yīng)的累積概率分布函數(shù)為Frk，則lbn和ubn的計(jì)算公式如下

(8)

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

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

2.2 軌道機(jī)動的檢測

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

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

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

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

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

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

3 算例分析

3.1 算例設(shè)計(jì)

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

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

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

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

3.2 結(jié)果分析

以27424號衛(wèi)星為例對機(jī)動檢測的詳細(xì)過程說明如下。利用2003年1月1日至2016年12月31日的歷史TLE數(shù)據(jù)，共得到128925個(gè)軌道參數(shù)預(yù)報(bào)誤差數(shù)據(jù)。這些數(shù)據(jù)按照其對應(yīng)的預(yù)報(bào)時(shí)間被分成196個(gè)數(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)星典型預(yù)報(bào)誤差數(shù)據(jù)組的柱狀圖Fig.3 Histograms of typical prediction error data groups for object 27424

隨后，這些預(yù)報(bào)誤差數(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)星典型預(yù)報(bào)誤差數(shù)據(jù)組進(jìn)行擬合得到的概率密度函數(shù)Fig.4 Fitted probability density functions of typical prediction error data groups for object 27424

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

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

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

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

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

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

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

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

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

4 結(jié)論

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