基于協(xié)方差理論的非關(guān)聯(lián)軌道動態(tài)關(guān)聯(lián)算法

，，，

1. 中國科學(xué)院 國家空間科學(xué)中心，北京 100190 2. 中國科學(xué)院大學(xué)，北京 100190

在國內(nèi)外編目實踐中，一些突發(fā)的情況會對正常的編目工作產(chǎn)生不利的影響，例如太陽爆發(fā)干擾了大氣，使得空間目標(biāo)的預(yù)報精度大大降低，也降低了對近地目標(biāo)的編目維持能力，由于極端空間環(huán)境條件下大氣模型精度有限，影響了大量的軌道預(yù)報精度，使得目標(biāo)丟失。一些突發(fā)的衛(wèi)星事件解體會在很短時間內(nèi)增加數(shù)以千計的可跟蹤目標(biāo)[1-3]。在這種情況下，設(shè)備監(jiān)測到的許多目標(biāo)可能沒有編目，非關(guān)聯(lián)軌道 (Uncorrelated Tracks，UCTs)數(shù)據(jù)和沒有編號的目標(biāo)大大增加，目標(biāo)編目數(shù)據(jù)處理中心的負(fù)擔(dān)立即增大。此外在目標(biāo)解體初期，目標(biāo)的空間密度還可能引起目標(biāo)混淆。這兩種情況必須在處理中心增大處理能力，以及增加額外的關(guān)聯(lián)過程來克服。大型衛(wèi)星星座的部署也會增加許多空間目標(biāo)。有些空間碎片可能已處于設(shè)備探測能力的邊緣，在觀測條件較好時，設(shè)備偶然能得到觀測數(shù)據(jù)，但是，觀測數(shù)據(jù)不足以維持該目標(biāo)的編目，因此目標(biāo)可能丟失。對于光學(xué)觀測設(shè)備，還有可見期的問題[4-10]。

解決這些問題除了增加探測網(wǎng)的探測能力及提高軌道預(yù)報精度外，還需要改進(jìn)關(guān)聯(lián)算法。采用基于協(xié)方差理論關(guān)聯(lián)算法可以提高短弧段雷達(dá)觀測數(shù)據(jù)和稀疏可見期數(shù)據(jù)的關(guān)聯(lián)水平。該計算方法旨在解決UCTs自動關(guān)聯(lián)問題，建立一個更完整的空間對象編目庫。基于協(xié)方差的軌道關(guān)聯(lián)動態(tài)算法是基于軌道不確定性協(xié)方差的動態(tài)算法，利用軌道數(shù)據(jù)的協(xié)方差信息生成完整的軌道狀態(tài)和誤差分布函數(shù)，并用于統(tǒng)計關(guān)聯(lián)軌道。如果協(xié)方差準(zhǔn)確表示實際軌道誤差分布，那么這個算法是最優(yōu)的、最大化的概率關(guān)聯(lián)。協(xié)方差矩陣可以更有效統(tǒng)計和自動方式關(guān)聯(lián)UCTs?；趨f(xié)方差的關(guān)聯(lián)已經(jīng)在精密軌道攝動的模擬中測試，可完成編目10萬個空間目標(biāo)[11]。仿真結(jié)果表明，相對于之前固定門限方法的軌道關(guān)聯(lián)，基于協(xié)方差的軌道關(guān)聯(lián)可降低錯誤關(guān)聯(lián)值達(dá)兩個數(shù)量級?；趨f(xié)方差的數(shù)據(jù)關(guān)聯(lián)方法將提高空間監(jiān)測網(wǎng)低軌道雷達(dá)監(jiān)測數(shù)據(jù)的編目分辨能力。

本文建立了線性和無跡卡爾曼濾波協(xié)方差生成和演化算法，分析了兩種協(xié)方差演化的精度，無跡卡爾曼濾波的協(xié)方差演化方法相對線性協(xié)方差演化方法有較高的演化精度。把協(xié)方差演化方法應(yīng)用到印度一箭104星中的20顆衛(wèi)星軌道關(guān)聯(lián)中，仿真結(jié)果表明采用基于協(xié)方差軌道關(guān)聯(lián)具有很好的關(guān)聯(lián)效果。

1 基于軌道協(xié)方差UCTs動態(tài)關(guān)聯(lián)算法模擬分析方法

基于協(xié)方差的UCTs動態(tài)關(guān)聯(lián)算法是把軌道協(xié)方差作為關(guān)聯(lián)時的影響因素，關(guān)聯(lián)量是復(fù)合等概率橢球。當(dāng)探測器的測量誤差是高斯分布并且軌道動力學(xué)為線性時，協(xié)方差關(guān)聯(lián)是最佳關(guān)聯(lián)量。采用此關(guān)聯(lián)方法，相同坐標(biāo)系下兩衛(wèi)星的狀態(tài)量和協(xié)方差須預(yù)報到同一時刻。例如，可以從先前t1時刻的軌道和協(xié)方差預(yù)報到t2時刻軌道，兩軌道狀態(tài)差為:

X2,1=X2(t2)-X1(t2)

(1)

式中：X1(t2)和X2(t2)為空間目標(biāo)在t2時刻的軌道狀態(tài)，包括位置和速度。關(guān)聯(lián)過程中通過計算統(tǒng)計量k來衡量兩衛(wèi)星的關(guān)聯(lián)程度。k屬于χ2統(tǒng)計，也被稱為馬氏距離：

k2=[X2,1]T[P1(t2)+P2(t2)]-1[X2,1]

(2)

式中：P1和P2分別為兩衛(wèi)星的軌道協(xié)方差矩陣。

為了實際模擬基于協(xié)方差的軌道關(guān)聯(lián)，需要模擬空間目標(biāo)監(jiān)測網(wǎng)的監(jiān)測模型和編目過程。決定從真實雷達(dá)測站出發(fā)，不同雷達(dá)測量數(shù)據(jù)隨機(jī)加入測量噪聲，模擬雷達(dá)測距及測角數(shù)據(jù)。定軌模塊對這些測量數(shù)據(jù)進(jìn)行初始軌道確定及最小二乘處理生成軌道及協(xié)方差數(shù)據(jù)。關(guān)聯(lián)模塊根據(jù)定軌目標(biāo)的協(xié)方差數(shù)據(jù)進(jìn)行演化實現(xiàn)與編目庫和非關(guān)聯(lián)目標(biāo)的關(guān)聯(lián)。

1.1 真實軌道

真實軌道是由數(shù)值軌道傳播器生成的。動力學(xué)方程模型包括太陽輻射光壓，空氣阻力密度模型，第三體的影響作用力和21×21階的重力場[12]。50個物體的軌道狀態(tài)取自軌道碎片編目。這些物體是從空間目標(biāo)目錄里半長軸為6860～6890之間抽樣的。目標(biāo)的橫截面積設(shè)為常值，面積和質(zhì)量與碎片編目的值相同。設(shè)置每一個物體的阻力系數(shù)和光壓反射系數(shù)為默認(rèn)值?？臻g目標(biāo)軌道攝動動力學(xué)模型采用的是EGM2008的21階地球重力場模型，NRLMSISE2000大氣阻力模型，球形光壓模型。依據(jù)軌道攝動方程演化計算理論上空間目標(biāo)的真實軌道[13]。

1.2 雷達(dá)探測器模擬

依據(jù)現(xiàn)有的雷達(dá)跟蹤測量站，模擬空間目標(biāo)的跟蹤測量數(shù)據(jù)。用國外公布的雷達(dá)參數(shù)對雷達(dá)跟蹤站進(jìn)行建模。表1所示這些參數(shù)來自美國巴利亞多(Vallado)[14]。對雷達(dá)的觀測數(shù)據(jù)進(jìn)行模擬，包括測距范圍、仰角、方位角、雷達(dá)高度、測距標(biāo)準(zhǔn)差、仰角標(biāo)準(zhǔn)差以及方位角標(biāo)準(zhǔn)差。在觀測數(shù)據(jù)中加入了高斯噪聲[15]。每一步都用標(biāo)準(zhǔn)差生成隨機(jī)偏差模擬高斯噪聲。采用這種標(biāo)準(zhǔn)差是因為根據(jù)一段時間對雷達(dá)探測器的分析顯示日常偏差標(biāo)準(zhǔn)差與測量標(biāo)準(zhǔn)差大致相等。如果觀測數(shù)據(jù)超出了方位角和仰角的限制就要去掉[16]。假定某個空間目標(biāo)截面積為1 m2時觀測站最大測量距離為ρmax。雷達(dá)的返回值與1/ρ4成正比，其中ρ是測距，所以觀測量應(yīng)滿足：

(3)

式中：A為空間目標(biāo)截面積，兩次測量之間為10 s間隔時，觀測數(shù)據(jù)設(shè)置限制在2′弧段內(nèi)。每個雷達(dá)探測器被分配最大數(shù)量的可跟蹤目標(biāo)。

表1 雷達(dá)探測器參數(shù)

1.3 軌道確定協(xié)方差

為了測試k值是否可以用來關(guān)聯(lián)編目軌道，選取物體的一個子集部分構(gòu)成了模擬目錄。選用目錄里每個物體的第一軌道，使用測距、方位角、仰角三要素方法進(jìn)行初始軌道確定。軌道狀態(tài)估計作為最小二乘批處理程序的先驗狀態(tài)。當(dāng)執(zhí)行最小二乘批差分校正時，每一個目錄物體所有前三天的觀測數(shù)據(jù)合并在一起。每個物體的面積和質(zhì)量都設(shè)為默認(rèn)值(1m2和100 kg)，但是要估計阻力系數(shù)及光壓反射參數(shù)的修正值[17-18]。在軌道計算中用數(shù)值計算軌道的傳播。在未來大于三天的時間里，物體估計狀態(tài)隨時間傳播且儲存為星歷表、用來做為一條空間物體目錄。50個空間物體只有30個創(chuàng)建了編目星歷，其他物體的星歷未知。

在先驗?zāi)夸泟?chuàng)建后，所有的真實軌道在未來24 h傳播。用真實軌道可以形成觀測數(shù)據(jù)，且已知物體的每個軌道都被設(shè)置了編號。最小二乘差分校正用來獲得每個軌道的狀態(tài)和協(xié)方差估計。由于這些軌道相當(dāng)?shù)亩?120 s)，這些單個軌道的估計有很大的不確定性，這在協(xié)方差上得到了反映。

1.4 數(shù)值積分模型的軌道協(xié)方差演化

1.4.1 誤差演化線性化方法

(4)

(5)

對式(5)微分：

(6)

(7)

e(t)=Φ(t,ti)ei+wd(t,ti)

(8)

其中，狀態(tài)轉(zhuǎn)換矩陣Φ(t,ti)滿足：

(9)

線性協(xié)方差轉(zhuǎn)換矩陣Φ(t,ti)計算過程中忽略了高階項的影響，只考慮了一階項。線性協(xié)方差演化是一種分析方法，其優(yōu)點是計算簡化高效。但是對于非線性程度高、協(xié)方差演化周期長或初始協(xié)方差較大情況其計算精度較低[19]，不適用于非高斯協(xié)方差演化。

1.4.2 UKF協(xié)方差演化方法

UKF是基于確定性樣點計算的濾波算法，以Unscented變換為基礎(chǔ)，采用卡爾曼線性濾波框架和確定性采樣方式進(jìn)行濾波。UKF采用一組有效的采樣點來將狀態(tài)近似為高斯隨機(jī)變量，捕捉高斯隨機(jī)變量的均值和協(xié)方差。當(dāng)狀態(tài)量通過實際的非線性系統(tǒng)之后，后驗均值和協(xié)方差可以精確到三階(泰勒展開)。由于不需要推導(dǎo)和計算Jacobian矩陣，基于UT的UKF算法更便于實現(xiàn)。

(1)構(gòu)造Sigma點

(10)

(2)利用Sigma點進(jìn)行非線性傳播

根據(jù)計算所得Sigma點直接按照離散化非線性模型進(jìn)行狀態(tài)和觀測預(yù)測，對應(yīng)產(chǎn)生的Sigma樣本點為：

χk(i)=f[χk-1(i),k]

(11)

式中：χk(i)(i=0，…,2n)為預(yù)測的對稱分布Sigma點。

(3)根據(jù)樣本點計算狀態(tài)均值和方差

(12)

(13)

從上面的3個步驟中可以看出，Unscented變換不需要計算非線性函數(shù)的雅可比矩陣，變換中僅涉及函數(shù)計算，即使系統(tǒng)的模型比較復(fù)雜的情況下，算法實現(xiàn)的難度也并不會增加。

1.5 UCTs關(guān)聯(lián)

處理這些軌道的第一步是比較星歷看軌道是否屬于已經(jīng)存在的編目物體，然后更新編目星歷。對于與已經(jīng)編目在案的軌道相符的軌道，采用軌道狀態(tài)值和協(xié)方差計算關(guān)聯(lián)量k值。如果k值小于限定閾值，則認(rèn)為該軌道匹配且軌道的標(biāo)識或者標(biāo)簽保持不變。如果k值大于限定閾值，則更改標(biāo)簽以表明它是不相關(guān)軌道即UCTs。然后就是把UCTs與所有編目里的物體進(jìn)行比較。為了減少比較數(shù)量，采用預(yù)測濾波的方法只選擇軌道根數(shù)大致相當(dāng)?shù)男菤v。如果比較值k小于限定閾值，則該軌道被標(biāo)記為屬于物體編目。接著，所有被分配到編目里的新軌道與先前確定的軌道結(jié)合形成新的星歷。這些綜合的觀測數(shù)據(jù)用來估計一個更新的編目軌道并創(chuàng)立一個新的、為期三天的預(yù)報編目星歷。下一步就是UCTs與UCTs之間的關(guān)聯(lián)，所有的UCTs形成一個候選軌道數(shù)據(jù)庫，不能與編目庫關(guān)聯(lián)的軌道再與UCTs候選軌道關(guān)聯(lián)。

2 模擬分析

2.1 協(xié)方差演化分析

低軌道目標(biāo)分別選取400 km軌道高度的近圓軌道目標(biāo)國際空間站和1 500 km軌道高度的近圓軌道目標(biāo)Ajisai衛(wèi)星。對國際空間站和Ajisai衛(wèi)星設(shè)定初值協(xié)方差，采用線性協(xié)方差演化方法和UKF協(xié)方差演化方法計算演化時間為一天時，空間站和衛(wèi)星的軌道協(xié)方差在RSW坐標(biāo)系下的分量。RSW坐標(biāo)系的R軸(徑向)是從地球指向航天器方向的矢量，S軸(橫向)在軌道平面內(nèi)垂直于R軸，W軸(垂跡)垂直于軌道面。把演化結(jié)果同STK軟件中高精度軌道(HPOP)協(xié)方差演化結(jié)果比較，分析兩種方法協(xié)方差演化精度。軌道攝動力主要考慮地球引力、大氣阻力、光壓攝動力及第三體引力。地球引力采用的是21階的EGM2008模型，大氣阻力采用的是NRLMSISE2000模型，大氣阻力系數(shù)Cd=2，光壓系數(shù)Cr=1.8。

假設(shè)在J2000坐標(biāo)系下空間站位置三方向初值誤差均為10 m，速度三方向誤差均為0.01 m/s?？臻g站協(xié)方差演化時長為一天，選取演化時長中的100個點，比較線性方法和UKF方法演化的100個點與HPOP協(xié)方差演化的誤差。圖1中前3個圖采用是兩種方法與HPOP對R、S、W三個方向的演化誤差進(jìn)行比較，最后1個圖是采用STK中的HPOP對空間站軌道協(xié)方差的演化結(jié)果。

相比較發(fā)現(xiàn)，線性方法在S方向一天的協(xié)方差演化誤差是米量級，UKF方法演化誤差是分米量級。UKF協(xié)方差演化誤差精度高于線性方法1個量級。

采用同樣方法對Ajisai衛(wèi)星軌道協(xié)方差演化，結(jié)果如圖2所示，相同的初值協(xié)方差A(yù)jisai衛(wèi)星通過HPOP方法協(xié)方差演化一天的值小于空間站的演化值。這主要是空間站所受的大氣阻力大于Ajisai衛(wèi)星。對于Ajisai衛(wèi)星S方向演化一天的誤差，線性方法是米量級，UKF方法演化誤差是厘米量級。UKF誤差精度高于線性方法兩個量級。對于低軌道空間目標(biāo)，軌道高度越高，UKF協(xié)方差演化誤差精度越高。

從圖1和圖2可以看出，對于空間站和Ajisai衛(wèi)星，通過UKF方法演化的R方向協(xié)方差誤差大于線性方法演化誤差。這是由于軌道協(xié)方差隨時間演化呈現(xiàn)非高斯性引起的。采用蒙特卡洛方法可以驗證當(dāng)軌道協(xié)方差演化時長為一天時，空間目標(biāo)的軌道協(xié)方差不再是高斯分布的橢球，而是呈月牙狀。月牙的彎曲方向就是在徑向(R方向)。圖3是采用蒙特卡洛方法演化的軌道協(xié)方差分布圖。對Ajisai衛(wèi)星初始根數(shù)加上位置標(biāo)準(zhǔn)差10 m、速度標(biāo)準(zhǔn)差0.01 m的隨機(jī)誤差，構(gòu)成10 000個采樣點，含隨機(jī)誤差的采樣點在與HPOP相同模型攝動力作用下向前傳播一天，把演化后的10 000個采樣點與標(biāo)準(zhǔn)精密軌道相差得到蒙特卡洛誤差采樣點分布。從圖3可以看出蒙特卡洛采樣點分布不再是標(biāo)準(zhǔn)橢球形，采樣點分布在R方向呈現(xiàn)彎曲的非高斯性分布。這種非高斯性彎曲分布導(dǎo)致R方向誤差協(xié)方差值變大，而UKF方法可以比較明顯地表示這種協(xié)方差演化特性。

2.2 UCTs關(guān)聯(lián)模擬

選取2017年2月印度發(fā)射的一箭104星中的20個目標(biāo)衛(wèi)星做UCTs關(guān)聯(lián)，衛(wèi)星軌道半長軸6 860～6 890 km，軌道傾角97.5°，衛(wèi)星軌道區(qū)域分布比較集中。各顆衛(wèi)星之間具有相似的軌道類型，在與箭體分離后一段時間各個衛(wèi)星在同一軌道面上，軌道之間只存在較小的相位差異，采用之前半長軸、軌道傾角篩選再進(jìn)行空間目標(biāo)編目存在困難，適用于驗證基于協(xié)方差算法軌道關(guān)聯(lián)。

每顆衛(wèi)星選取2～4個不同時間點的軌道星歷，軌道間隔小于兩天。采用美國巴利亞多(Vallado)提供的雷達(dá)參數(shù)模擬20個目標(biāo)衛(wèi)星的雷達(dá)觀測數(shù)據(jù)。并用HerrickGibbs方法對短弧段各個衛(wèi)星初始軌道進(jìn)行定軌，用最小二乘法對軌道進(jìn)行精密定軌，生成衛(wèi)星在某點的軌道根數(shù)及協(xié)方差信息。應(yīng)用軌道預(yù)報方法對軌道根數(shù)進(jìn)行數(shù)值預(yù)報，通過線性協(xié)方差演化方法和UKF協(xié)方差演化方法對衛(wèi)星初始協(xié)方差進(jìn)行演化。

20個UCTs模擬真實軌道是基于space-track提供的兩行軌道根數(shù)，由于公布的兩行軌道根數(shù)不包含軌道協(xié)方差信息，因此假設(shè)其為真實軌道重新定軌生成協(xié)方差。而兩行軌道根數(shù)本身含有軌道誤差，這可能是關(guān)聯(lián)過程中關(guān)聯(lián)閾值較大的主要原因。同時兩行軌道根數(shù)是以真赤道平春分點(True Equator Mean Equinox,TEME)作為參考坐標(biāo)系，需要轉(zhuǎn)變?yōu)镴2000坐標(biāo)系分析。

圖4是20個衛(wèi)星中，相同衛(wèi)星不同歷元軌道根數(shù)進(jìn)行協(xié)方差關(guān)聯(lián)的關(guān)聯(lián)量k值的統(tǒng)計分布情況?？梢钥闯?0顆衛(wèi)星UCTs關(guān)聯(lián)量分布主要在60～80。分析發(fā)現(xiàn)當(dāng)UCTs關(guān)聯(lián)閾值限定為200時，20顆衛(wèi)星中相同衛(wèi)星不同歷元的軌道根數(shù)都可以成功關(guān)聯(lián)。關(guān)聯(lián)值大的原因作者分析如下：

1)首先文章中選用美國發(fā)布的兩行軌道根數(shù)，這種根數(shù)本身存在較大的軌道誤差。由于模擬中選取20顆衛(wèi)星的兩行軌道根數(shù)作為真實軌道，忽略了兩行軌道根數(shù)數(shù)據(jù)本身所含有的誤差，所以計算的關(guān)聯(lián)閾值較大。

2)在雷達(dá)觀測數(shù)據(jù)模擬過程中根據(jù)不同雷達(dá)測站的噪聲加入了測量誤差。

3)在精密定軌過程中由于弧段較短時不能進(jìn)行精密定軌確定協(xié)方差，因此在一些衛(wèi)星之間的關(guān)聯(lián)模擬過程中選用了較長的弧段進(jìn)行定軌，弧段長定軌協(xié)方差小，而馬氏距離計算時協(xié)方差的逆較大，關(guān)聯(lián)值較大。

假設(shè)選取精密的軌道根數(shù)作為真實軌道，短弧段情況下可以計算出軌道協(xié)方差，則關(guān)聯(lián)閾值會大幅下降。

圖5表示相同衛(wèi)星不同歷元間隔對應(yīng)的關(guān)聯(lián)量?？梢钥闯?，隨著軌道時間間隔增加關(guān)聯(lián)量逐漸下降，相同軌道時間間隔，不同衛(wèi)星軌道所計算的關(guān)聯(lián)量也不相同。關(guān)于關(guān)聯(lián)閾值的選取需要針對不同衛(wèi)星軌道做進(jìn)一步的分析。

3 結(jié)束語

本文研究了基于協(xié)方差理論的UCTs動態(tài)關(guān)聯(lián)算法，對比分析了基于雅克比轉(zhuǎn)換的線性協(xié)方差演化方法和基于UT轉(zhuǎn)換的UKF協(xié)方差演化方法。以400 km和1k500kkm軌道高度為例，初始軌道協(xié)方差相同的情況下，基于UT變換的協(xié)方差演化精度要高于基于雅克比轉(zhuǎn)換的線性協(xié)方差方法演化精度。應(yīng)用UKF方法軌道協(xié)方差演化一天的位置誤差分別是分米量級和厘米量級，應(yīng)用線性方法演化一天的位置誤差是米級。

把協(xié)方差演化方法應(yīng)用到UCTs動態(tài)關(guān)聯(lián)中去，模擬印度一箭104星中的20顆衛(wèi)星的軌道關(guān)聯(lián)，當(dāng)軌道關(guān)聯(lián)量閾值選取為200時，20顆衛(wèi)星軌道可以全部關(guān)聯(lián)。

軌道協(xié)方差演化有很好的應(yīng)用前景，本文在空間目標(biāo)編目方面，為處理各非關(guān)聯(lián)軌道提出了一種動態(tài)算法，其關(guān)聯(lián)量源于每個非關(guān)聯(lián)軌道的軌道誤差協(xié)方差矩陣。另外空間目標(biāo)碰撞評估的一個重要環(huán)節(jié)是計算兩個空間物體的碰撞概率。碰撞概率的計算使用了相對狀態(tài)概率密度函數(shù)的數(shù)值積分算法。每種算法需要為兩個對象輸入狀態(tài)向量的誤差信息。狀態(tài)向量的誤差信息通過協(xié)方差矩陣的方式來表示。本研究可以為將來的空間目標(biāo)編目和空間碎片預(yù)警提供方法參考。