，，，

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

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

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

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

(1)

式中：X1(t2)和X2(t2)為空間目標在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)，需要模擬空間目標監(jiān)測網(wǎng)的監(jiān)測模型和編目過程。決定從真實雷達測站出發(fā)，不同雷達測量數(shù)據(jù)隨機加入測量噪聲，模擬雷達測距及測角數(shù)據(jù)。定軌模塊對這些測量數(shù)據(jù)進行初始軌道確定及最小二乘處理生成軌道及協(xié)方差數(shù)據(jù)。關(guān)聯(lián)模塊根據(jù)定軌目標的協(xié)方差數(shù)據(jù)進行演化實現(xiàn)與編目庫和非關(guān)聯(lián)目標的關(guān)聯(lián)。

1.1 真實軌道

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

1.2 雷達探測器模擬

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

(3)

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

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

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

為了測試k值是否可以用來關(guān)聯(lián)編目軌道，選取物體的一個子集部分構(gòu)成了模擬目錄。選用目錄里每個物體的第一軌道，使用測距、方位角、仰角三要素方法進行初始軌道確定。軌道狀態(tài)估計作為最小二乘批處理程序的先驗狀態(tài)。當執(zhí)行最小二乘批差分校正時，每一個目錄物體所有前三天的觀測數(shù)據(jù)合并在一起。每個物體的面積和質(zhì)量都設(shè)為默認值(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é)方差估計。由于這些軌道相當?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ǔ)，采用卡爾曼線性濾波框架和確定性采樣方式進行濾波。UKF采用一組有效的采樣點來將狀態(tài)近似為高斯隨機變量，捕捉高斯隨機變量的均值和協(xié)方差。當狀態(tài)量通過實際的非線性系統(tǒng)之后，后驗均值和協(xié)方差可以精確到三階(泰勒展開)。由于不需要推導(dǎo)和計算Jacobian矩陣，基于UT的UKF算法更便于實現(xiàn)。

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

(10)

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

根據(jù)計算所得Sigma點直接按照離散化非線性模型進行狀態(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值小于限定閾值，則認為該軌道匹配且軌道的標識或者標簽保持不變。如果k值大于限定閾值，則更改標簽以表明它是不相關(guān)軌道即UCTs。然后就是把UCTs與所有編目里的物體進行比較。為了減少比較數(shù)量，采用預(yù)測濾波的方法只選擇軌道根數(shù)大致相當?shù)男菤v。如果比較值k小于限定閾值，則該軌道被標記為屬于物體編目。接著，所有被分配到編目里的新軌道與先前確定的軌道結(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é)方差演化分析

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

假設(shè)在J2000坐標系下空間站位置三方向初值誤差均為10 m，速度三方向誤差均為0.01 m/s?？臻g站協(xié)方差演化時長為一天，選取演化時長中的100個點，比較線性方法和UKF方法演化的100個點與HPOP協(xié)方差演化的誤差。圖1中前3個圖采用是兩種方法與HPOP對R、S、W三個方向的演化誤差進行比較，最后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誤差精度高于線性方法兩個量級。對于低軌道空間目標，軌道高度越高，UKF協(xié)方差演化誤差精度越高。

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

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

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

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

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

圖4是20個衛(wèi)星中，相同衛(wèi)星不同歷元軌道根數(shù)進行協(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)當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)在雷達觀測數(shù)據(jù)模擬過程中根據(jù)不同雷達測站的噪聲加入了測量誤差。

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

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

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

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)，當軌道關(guān)聯(lián)量閾值選取為200時，20顆衛(wèi)星軌道可以全部關(guān)聯(lián)。

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