于彥君，岳程斐，李化義，陳雪芹，劉培

（1.哈爾濱工業(yè)大學(xué)衛(wèi)星技術(shù)研究所，哈爾濱 150001；2.哈爾濱工業(yè)大學(xué)（深圳）空間科學(xué)與應(yīng)用技術(shù)研究院，深圳 518055；3.上海衛(wèi)星工程研究所，上海 710119）

0 引言

軌道不確定性傳播是指在給定軌道初始狀態(tài)及其相關(guān)不確定性的基礎(chǔ)上，預(yù)測(cè)未來(lái)的軌道狀態(tài)及其統(tǒng)計(jì)特性。軌道不確定性傳播在空間態(tài)勢(shì)感知與決策等相關(guān)任務(wù)中起著重要作用，例如軌跡跟蹤、碰撞概率分析、傳感器資源管理和異常檢測(cè)等［1］。對(duì)于航天器而言，其軌道狀態(tài)矢量的真值是無(wú)法確切知道的，當(dāng)存在測(cè)量信息時(shí)，軌道狀態(tài)誤差可以通過(guò)濾波算法［2］進(jìn)行估計(jì)，該過(guò)程通常是收斂的；當(dāng)不存在測(cè)量信息時(shí)，可利用濾波算法的預(yù)測(cè)步驟進(jìn)行不確定性估計(jì)，即進(jìn)行軌道不確定性傳播，該過(guò)程通常是發(fā)散的。在軌道不確定性傳播問(wèn)題中，初始軌道不確定性在任何坐標(biāo)系中都不是絕對(duì)高斯的，但目前通常采用高斯假設(shè)。高斯誤差超橢圓體隨著傳播時(shí)間的增加而放大和旋轉(zhuǎn)，并且由于動(dòng)力學(xué)的非線性性質(zhì)而逐漸變?yōu)榉歉咚沟?。在笛卡爾坐?biāo)系中，非高斯軌道不確定性通常位于“香蕉型”弧線上，該弧線沿航天器的標(biāo)稱軌跡彎曲［3］。軌道不確定性的典型傳播過(guò)程如圖1所示。

圖1 軌道不確定性傳播示意圖［1］Fig.1 Illustration of uncertainty propagation process［1］

軌道不確定性傳播方法包含局部線性化方法和蒙特卡洛法（Monte Carlo method，MC）等［4-5］。然而，由于軌道動(dòng)力學(xué)的非線性特質(zhì)，局部線性化方法過(guò)低的精度和蒙特卡洛方法的高計(jì)算量極大地阻礙了這2種方法的應(yīng)用。為了在降低計(jì)算量的同時(shí)保證合適的精度，針對(duì)不同場(chǎng)景，學(xué)者們提出了多種適用于非線性軌道不確定性傳播的方法，例如高斯和模型（Gaussian mixture model，GMM）［6-7］、無(wú)跡變換法（Unscented transformation，UT）［8］、狀態(tài)轉(zhuǎn)移張量法［9］、微分代數(shù)法［10］、多項(xiàng)式展開(kāi)法［11］、坐標(biāo)轉(zhuǎn)換法［12-13］以及將上述方法相互結(jié)合的混合方法等［14-17］。在這些方法中，UT法［18］解決了線性化方法和蒙特卡洛法不能兼顧計(jì)算效率和精度這一難題。其通過(guò)非線性積分從初始分布中確定性地選擇樣本來(lái)近似未來(lái)時(shí)間的概率分布以取代復(fù)雜的高階分析。這一過(guò)程中，通過(guò)近似概率分布代替近似非線性變換，提高了計(jì)算精度；通過(guò)較少的樣本計(jì)算代替MC 法中的大樣本計(jì)算，提高了計(jì)算效率。坐標(biāo)轉(zhuǎn)換法則通過(guò)選擇在非線性較弱的坐標(biāo)系中建立動(dòng)力學(xué)模型，在相同計(jì)算量條件下，達(dá)到大幅延長(zhǎng)不確定性真實(shí)度（Uncertainty realism，UR）［19］保持時(shí)間的效果。坐標(biāo)轉(zhuǎn)換法與UT 法相結(jié)合可用于軌道不確定性的傳播，如文獻(xiàn)［20］和［21］。

隨著低軌巨型衛(wèi)星星座的發(fā)展，衛(wèi)星間碰撞概率顯著增加，尤其當(dāng)單顆衛(wèi)星同時(shí)面臨與多個(gè)衛(wèi)星碰撞的風(fēng)險(xiǎn)時(shí)，快速估計(jì)分析軌道不確定性傳播結(jié)果并評(píng)估碰撞風(fēng)險(xiǎn)成為亟待突破的問(wèn)題。特別是考慮星上有限計(jì)算能力條件下，在確保不確定性估計(jì)精度的基礎(chǔ)上，快速獲得某一特定時(shí)刻或時(shí)段內(nèi)的軌道狀態(tài)及其統(tǒng)計(jì)特性具有重要的應(yīng)用價(jià)值。

針對(duì)這一問(wèn)題，本文提出了一種基于半解析法的軌道不確定性估計(jì)方法。該方法是一種混合方法，具體特點(diǎn)如下：首先選擇在軌道要素坐標(biāo)系中進(jìn)行不確定性傳播，延長(zhǎng)UR 的保持時(shí)間，跟蹤笛卡爾坐標(biāo)系中非高斯形式的軌道不確定性。其次，采用解析模型結(jié)合迭代法實(shí)現(xiàn)對(duì)任意時(shí)刻的軌道要素的快速計(jì)算，并將該方法應(yīng)用于軌道不確定性估計(jì)，大幅提升計(jì)算效率。此外，采用球形單邊采樣策略而非目前常用的對(duì)稱采樣策略，進(jìn)一步降低由采樣點(diǎn)數(shù)量帶來(lái)的計(jì)算量，并結(jié)合球形單邊采樣策略特性和高斯攝動(dòng)方程對(duì)采樣點(diǎn)進(jìn)行修正，再估計(jì)傳播后的軌道不確定性，以彌補(bǔ)采用解析模型傳播樣本點(diǎn)導(dǎo)致的精度降低。最后，將提出的算法與高斯和模型相結(jié)合，進(jìn)一步改善軌道不確定性估計(jì)精度。

1 基于半解析法的軌道要素快速計(jì)算方法

不確定性傳播算法的有效性可以通過(guò)使用非線性程度較低的坐標(biāo)系得到顯著提高。例如，相較于笛卡爾坐標(biāo)系中的衛(wèi)星軌道動(dòng)力學(xué)，采用軌道要素作為動(dòng)力學(xué)模型能夠吸收非線性動(dòng)力學(xué)中最主要的項(xiàng)（即引力中的r-2項(xiàng)），使得軌道要素空間的高斯分布能夠更準(zhǔn)確地捕捉到在笛卡爾空間中經(jīng)常觀察到的非高斯行為，并能夠更長(zhǎng)時(shí)間地保持UR［22］。本文采用經(jīng)典軌道要素進(jìn)行軌道不確定性傳播，針對(duì)低軌衛(wèi)星給出了考慮J2 攝動(dòng)與大氣阻力攝動(dòng)的軌道要素模型和解析解。

1.1 軌道要素長(zhǎng)期項(xiàng)與周期項(xiàng)建模

瞬時(shí)軌道要素通?？煞譃殚L(zhǎng)期項(xiàng)、長(zhǎng)周期項(xiàng)和短周期項(xiàng)。而上述長(zhǎng)期項(xiàng)和周期項(xiàng)可以基于高斯攝動(dòng)方程［23］（Gauss variational equations，GVE）進(jìn)行推導(dǎo)。本文采用NTW坐標(biāo)系下的描述形式。NTW坐標(biāo)系以航天器為中心，T軸與軌道相切，以速度矢量所指方向?yàn)檎较?；W軸沿角動(dòng)量矢量方向；N軸垂直于T軸和W軸，方向由右手定則確定。當(dāng)0 ＜e＜1，i≠0，i≠180°時(shí)，對(duì)應(yīng)GVE為：

令參數(shù)向量α?[a，e，i，Ω，ω，M]T，則在考慮J2攝動(dòng)和指數(shù)模型大氣阻力的影響下，長(zhǎng)期項(xiàng)的時(shí)間導(dǎo)數(shù)［24］為：

式中：ρ0為軌道半徑r0處的大氣密度；K2=(CDS2/m)，K1=(CDS1/m)Q；H為標(biāo)高；Q=(1 -rPωe；S1和S2分別是垂直于T軸和N軸方向 的橫截面積；rP和vP分別代表近地點(diǎn)半徑和近地點(diǎn)速度；ωe=7.292 115 855 3 × 10-5rad·s-1為地球角速度，RE為地球半徑。

類比文獻(xiàn)［26］中的推導(dǎo)方式，將文獻(xiàn)［24］中給出的大氣阻力長(zhǎng)周期項(xiàng)遞推方程進(jìn)一步推導(dǎo)，可得大氣阻力攝動(dòng)長(zhǎng)周期項(xiàng)解析解αdl（α）為：

式中：ni和nω分別為軌道傾角和近心點(diǎn)角的平均變化率，近似為：

J2 攝動(dòng)影響下的一階長(zhǎng)周期項(xiàng)αl（α）、短周期項(xiàng)αs（α）和大氣阻力攝動(dòng)影響下的短周期項(xiàng)αds（α）在參考文獻(xiàn)［24-25］中已給出，在此不再贅述。

1.2 任意時(shí)刻軌道要素快速計(jì)算方法

假設(shè)初始時(shí)刻為t0，需要獲取t時(shí)刻的瞬時(shí)軌道要素。一般情況下，可通過(guò)數(shù)值方法，如龍格庫(kù)塔法對(duì)式（1）直接進(jìn)行積分獲取。該方式需要從t0時(shí)刻逐步積分到t時(shí)刻。然而，在實(shí)際情況下，通常只對(duì)未來(lái)某一時(shí)刻或時(shí)段內(nèi)的軌道要素感興趣。例如，在進(jìn)行碰撞概率計(jì)算時(shí)，在整個(gè)軌道周期內(nèi)計(jì)算碰撞概率顯然是不必要的，通常需要首先進(jìn)行碰撞時(shí)間區(qū)間估計(jì)，隨后僅在該時(shí)間區(qū)間內(nèi)進(jìn)行碰撞分析［27］。

基于軌道要素長(zhǎng)期項(xiàng)和周期項(xiàng)模型能夠快速計(jì)算任意時(shí)刻軌道要素。在1.1節(jié)中將瞬時(shí)軌道要素分解為長(zhǎng)期項(xiàng)，長(zhǎng)周期項(xiàng)和短周期項(xiàng)。由于長(zhǎng)期項(xiàng)是線性項(xiàng)，t時(shí)刻的軌道要素長(zhǎng)期項(xiàng)，即t時(shí)刻平均軌道要素約為：

其中，上標(biāo)表示對(duì)應(yīng)時(shí)刻，αˉt0可通過(guò)下式獲得：

對(duì)于軌道要素周期項(xiàng)的解析解，需要采用瞬時(shí)軌道要素作為自變量，而式（8）給出的是t時(shí)刻平均軌道要素。因此為了提高精度，t時(shí)刻瞬時(shí)軌道要素可通過(guò)迭代法求解，其步驟為：

1）依據(jù)t時(shí)刻平均軌道要素計(jì)算瞬時(shí)軌道要素初值：

2 軌道不確定性半解析快速估計(jì)方法

在第1節(jié)中給出了快速計(jì)算任意時(shí)刻軌道要素的方法。因此，在給定t0時(shí)刻瞬時(shí)軌道要素的均值和協(xié)方差后（假設(shè)為高斯噪聲），即可選取特定的西格瑪點(diǎn)，通過(guò)改進(jìn)UT 變換進(jìn)行軌道不確定性傳播。方便起見(jiàn)，將上述算法稱為改進(jìn)半解析無(wú)跡變換（Improved semianalytical unscented transformation，ISUT）軌道不確定性估計(jì)方法；ISUT 算法與GMM 結(jié)合能夠進(jìn)一步提高精度，將該算法稱為GMMISUT算法。

2.1 球形單邊采樣策略

目前采用UT 法進(jìn)行軌道不確定性傳播的文獻(xiàn)通常采用傳統(tǒng)對(duì)稱采樣策略，即對(duì)于一個(gè)ns維空間，需要選擇2ns+1 個(gè)西格瑪點(diǎn)?？紤]到UT 的計(jì)算成本與使用的西格瑪點(diǎn)數(shù)成正比，對(duì)于高維問(wèn)題，最小化西格瑪點(diǎn)的數(shù)量可以最小化計(jì)算成本。文獻(xiàn)［18，28］中給出了球形單邊采樣策略，它僅需ns+2個(gè)西格瑪點(diǎn)即可捕獲與傳統(tǒng)對(duì)稱采樣策略相當(dāng)?shù)姆植夹畔ⅲ_(dá)到與截?cái)喽A濾波器相同的預(yù)測(cè)能力。因此，本文采用球形單邊采樣策略進(jìn)行西格瑪點(diǎn)的選取。球形單邊采樣方法通過(guò)迭代的方式獲得采樣點(diǎn)，其步驟可總結(jié)如下：

1）選擇位于原點(diǎn)（ι=0）和以原點(diǎn)為中心的超球體（ι≠0）上的采樣點(diǎn)權(quán)重：

2）使用比例無(wú)跡變換調(diào)整權(quán)重：

式中：比例因子σ∈（0，1］。

3）定義各個(gè)采樣點(diǎn)權(quán)重：

式中：β是影響第0采樣點(diǎn)協(xié)方差權(quán)重的參數(shù)。

4）當(dāng)維數(shù)大于1 維時(shí)，通過(guò)式（14）迭代獲得采樣點(diǎn)：

其中，各個(gè)向量序列初值為：

2.2 ISUT方法

根據(jù)2.1小節(jié)中的球形單邊采樣策略獲得采樣點(diǎn)后，選擇初始西格瑪點(diǎn)集：

采用1.2節(jié)中的任意時(shí)刻軌道要素快速計(jì)算方法傳播初始西格瑪點(diǎn)集，可獲得t時(shí)刻的西格瑪點(diǎn)集xι（t）。由于該方法計(jì)算西格瑪點(diǎn)集所采用的動(dòng)力學(xué)基于GVE 的解析解，存在一定的截?cái)嗾`差。為了使得計(jì)算結(jié)果更接近于GVE 的結(jié)果，需要對(duì)xι（t）的分布進(jìn)行調(diào)整，即對(duì)采樣點(diǎn)進(jìn)行修正。

由式（14）和式（16）可知，位于超球體原點(diǎn)的第0 西格瑪點(diǎn)x0(t0)=為不確定性分布的均值。而其余ns+1 個(gè)西格瑪點(diǎn)位于以原點(diǎn)為中心的超球體上。因此，可通過(guò)軌道遞推獲得GVE 動(dòng)力學(xué)對(duì)應(yīng)的t時(shí)刻的第0 西格瑪點(diǎn)，計(jì)算其與x0（t）的差值Δx，并依據(jù)軌道要素間的關(guān)系對(duì)解析方式獲得的西格瑪點(diǎn)集xι（t）進(jìn)行修正，從而減小誤差。

本文所采用的動(dòng)力學(xué)模型中平近點(diǎn)角M為快變量，其余5 個(gè)軌道要素為慢變量。在進(jìn)行不確定性傳播時(shí)，M發(fā)散速度較快，其余5 變量發(fā)散速度較慢。因此，對(duì)于各西格瑪點(diǎn)的前5個(gè)軌道要素，采用依據(jù)差值Δx進(jìn)行平移的修正策略，即：

對(duì)于發(fā)散較快的M，設(shè)計(jì)如下修正策略進(jìn)行補(bǔ)償：

基于修正后的西格瑪點(diǎn)集，t時(shí)刻的瞬時(shí)軌道要素均值和誤差協(xié)方差矩陣為：

則ISUT 軌道不確定性估計(jì)算法的流程如圖2所示。

圖2 ISUT軌道不確定性估計(jì)算法框圖Fig.2 Diagram of ISUT uncertainty propagation algorithm

最終輸出的t時(shí)刻瞬時(shí)軌道要素均值和誤差協(xié)方差矩陣將用于構(gòu)建t時(shí)刻軌道不確定性的概率密度函數(shù)（Probability density function，PDF），即：

式中：N（x；μU，PU）表示正態(tài)分布；x，μU=和PU=Pt分別表示輸入變量、均值和誤差協(xié)方差矩陣。

2.3 GMMISUT方法

根據(jù)式（20）可知，ISUT 法傳播的不確定性PDF是高斯形式的，即僅能夠傳播前2階統(tǒng)計(jì)矩（均值與方差）。對(duì)于非線性系統(tǒng)，在進(jìn)行不確定性傳播時(shí)，不確定性通常會(huì)呈現(xiàn)出非線性特性，即包含了高階統(tǒng)計(jì)矩，如峰度和偏度。因此，為了能夠進(jìn)一步提高半解析法的估計(jì)精度，可在初始時(shí)刻將高斯形式的不確定性近似為GMM。GMM 是對(duì)高斯PDF 的直接拓展，其形式為多個(gè)高斯PDF的加權(quán)總和：

式中：p（x）為待分解的PDF為p（x）的GMM 近似；ηG，μG，PG分別表示第G分量的權(quán)重、均值和協(xié)方差。NGM表示GMM 組件的數(shù)量。為了確保PDF 的非負(fù)性，并確保PDF 曲線下的總面積為1，權(quán)重必須是非負(fù)的，且總和為1，即：

采用GMM 方法對(duì)PDF 進(jìn)行描述在保留高斯PDF 特性的同時(shí)擴(kuò)展了高斯PDF 的適用性。GMM可以近似多種多樣的PDF，文獻(xiàn)［29］中證明了采用GMM 對(duì)PDF 進(jìn)行近似時(shí)，當(dāng)GMM 組件數(shù)量趨向無(wú)窮大時(shí)均勻收斂。這是一個(gè)非常直觀的結(jié)果，因?yàn)殡S著組件協(xié)方差趨近于零，GMM 的每個(gè)組件都趨向于脈沖函數(shù)。因此，通過(guò)減小組件的協(xié)方差、增加組件的數(shù)量，并適當(dāng)分布組件的均值，可以輕松地近似大部分PDF，且在進(jìn)行不確定傳播時(shí)能夠同時(shí)傳播低階和高階統(tǒng)計(jì)矩，進(jìn)一步提升不確定性估計(jì)精度。

目前，拆分高斯PDF 有2 種方法，第1 種為在多個(gè)方向進(jìn)行拆分，第2 種為使用單變量拆分庫(kù)在一個(gè)特定的方向進(jìn)行拆分［30］。采用第1種方法通常需要求解L2 最小化優(yōu)化問(wèn)題，會(huì)引入額外的計(jì)算量。因此，本文采用第2種方法，將多變量高斯PDF在其最大特征值方向進(jìn)行拆分，步驟如下：

1）將標(biāo)準(zhǔn)正態(tài)分布p（x；0，1）拆分為NGM個(gè)均質(zhì)高斯分布組件：

表1 3組件拆分庫(kù)Table 1 Three-component splitting library

圖3 3組件拆分庫(kù)組件及其總和與標(biāo)準(zhǔn)高斯分布的比較Fig.3 Components of the 3-component splitting libraries and their sum as compared with the standard Gaussian distribution

2）采用譜分解將協(xié)方差矩陣PU分解為：

式中：Λ是PU的特征值矩陣，λk是第k大的特征值。選擇沿第k個(gè)特征向量的方向分解p（x；μU，PU），則式（21）中的參數(shù)計(jì)算方法為：

式中：υk是PU的第k個(gè)特征向量。對(duì)于每個(gè)GMM 組件，均采用ISUT算法進(jìn)行傳播，即構(gòu)成GMMISUT 方法。在應(yīng)用GMMISUT 實(shí)現(xiàn)不確定性傳播的過(guò)程中，應(yīng)合理選擇NGM從而平衡計(jì)算效率與精度。

3 數(shù)值仿真

3.1 軌道不確定性真實(shí)度的似然一致性度量

似然一致性度量L 描述了2 個(gè)PDF 之間的重疊程度，對(duì)于彼此更一致的PDF來(lái)說(shuō)，L會(huì)更大。在將數(shù)據(jù)點(diǎn)樣本與PDF進(jìn)行比較的情況下，2個(gè)PDF之間的似然一致性度量［32］（Likelihood agreement measure，LAM）定義為：

由于需要分析軌道不確定性PDF 與蒙特卡洛法生成的樣本的一致性，可假設(shè)q(x)由狄拉克模型（Dirac mixture model，DMM）給出：

式中：δ(x-μj)是以μj為中心的Dirac-δ分布，權(quán)重γj=K-1。若p（x）為高斯PDF，將式（20）和式（27）代入式（26），可得：

若p（x）為GMM的PDF，將式（21）和式（27）代入式（26），可得：

因此，通過(guò)DMM 給定1 組樣本點(diǎn)，并與高斯/高斯和模型進(jìn)行比較，DMM 與高斯/高斯和模型具有相同分布的可能性可以分別通過(guò)式（28）和式（29）計(jì)算。LAM的值越高，則給定的高斯/高斯和模型更有可能生成DMM。該方法允許將多個(gè)高斯/高斯和模型與1 組DMM 進(jìn)行比較。其中最準(zhǔn)確的高斯/高斯和模型具有最高的LAM值。

3.2 初始條件和參數(shù)設(shè)置

衛(wèi)星的傳感器通常能夠提供衛(wèi)星的瞬時(shí)位置、速度和時(shí)間信息，且無(wú)需處理原始偽距或載波信息。為了更貼近真實(shí)情況，在地心慣性坐標(biāo)系中給出衛(wèi)星初始測(cè)量誤差協(xié)方差及衛(wèi)星位置速度均值。假設(shè)地心慣性坐標(biāo)系中衛(wèi)星的三軸位置和速度誤差標(biāo)準(zhǔn)差分別為100 m 和1 m/s。位置均值rECI=［2 505.35 -6 439.95 1 857.00］ km，速度均值vECI=［2.81 -0.96 -6.84］ km/s。選取樣本點(diǎn)數(shù)目為MC=10 000，通過(guò)坐標(biāo)轉(zhuǎn)換［22］，可得初始瞬時(shí)軌道要素的均值和誤差協(xié)方差如表2所示。

表2 瞬時(shí)軌道要素的初始均值與協(xié)方差Table 2 Initial mean and covariance of osculating orbital elements

假設(shè)初始大氣密度ρ0=2.34 × 10-13kg/m3，密度標(biāo)高H=68.7 km，阻力系數(shù)CD=2.2。航天器垂直于切線方向截面積S1=0.2 m2，垂直于次法線方向的截面積S2=0.2 m2，航天器質(zhì)量m=10 kg。其他參數(shù)：σ=1，β=2，W0=0.25。

3.3 仿真結(jié)果與分析

為了比較UT、ISUT和GMMISUT法的計(jì)算效率，定義UT 法和ISUT 法的計(jì)算消耗時(shí)間之比εT（t）、UT法和GMMISUT法的計(jì)算消耗時(shí)間之比εTG（t）為：

式中：τUT（t）為采用UT算法獲取t時(shí)刻軌道不確定性所需的時(shí)間，τISUT（t）為采用ISUT 法獲取t時(shí)刻軌道不確定性所需的時(shí)間，τGMMISUT（t）為采用GMMISUT法獲取t時(shí)刻軌道不確定性所需的時(shí)間。顯然，當(dāng)εT(t) ＞1時(shí)，ISUT 法的計(jì)算效率高于UT法，當(dāng)εTG(t) ＞1時(shí)，GMMISUT法的計(jì)算效率高于UT法。

采用ODE45 方法執(zhí)行UT 法的積分步驟，并與ISUT 法和GMMISUT 法進(jìn)行比較，結(jié)果如圖4 所示。從圖4 中可以看出，當(dāng)t＜1d時(shí)，所需獲取的軌道不確定性所在時(shí)刻t越大，ISUT 算法和GMMISUT 算法的計(jì)算效率相較于UT算法越高。當(dāng)t≥1d后，εT（t）和εTG（t）趨于穩(wěn)定，ISUT 法的計(jì)算效率約為UT 法的7.9倍，GMMISUT法的計(jì)算效率約為UT法的2.6倍。本文中對(duì)比分析所采用的UT 法采用了單邊球形采樣策略，采樣點(diǎn)數(shù)目為8個(gè)，均需進(jìn)行積分；對(duì)于ISUT法，采樣點(diǎn)數(shù)目為8個(gè)，需對(duì)其中1 個(gè)采樣點(diǎn)進(jìn)行積分；對(duì)于GMMISUT法，采樣點(diǎn)數(shù)目為24個(gè)，需對(duì)其中3 個(gè)采樣點(diǎn)進(jìn)行積分。因此可知，上述算法的計(jì)算負(fù)擔(dān)大小主要取決于需要進(jìn)行積分的采樣點(diǎn)數(shù)目。采用ISUT法和GMMISUT 法能夠減小需要數(shù)值積分計(jì)算的點(diǎn)的數(shù)目，因此計(jì)算效率更高。

圖4 εT（t）和εTG（t）變化情況Fig.4 The variation of εT（t）and εTG（t）

為了分析采用UT、SUT、ISUT 和GMMISUT 法進(jìn)行軌道不確定性傳播的精度，根據(jù)3.1 小節(jié)中給出的LAM 計(jì)算方法，分別計(jì)算4種算法的LAM 值。在對(duì)比分析UT法，SUT 法和ISUT 法時(shí)，以ISUT 法的LAM 值為標(biāo)準(zhǔn)進(jìn)行歸一化，則3 種算法的歸一化LAM 值隨時(shí)間變化情況如圖5 所示?？梢钥闯?，沒(méi)有進(jìn)行修正的SUT 法的歸一化LAM 值始終較低，對(duì)于軌道不確定性跟蹤效果較差。在第1 天內(nèi)，ISUT法的精度與UT 法相近，甚至整體略優(yōu)于UT 法。而在第2天之后，ISUT法的精度相較于UT法有顯著的降低。ISUT 算法的誤差大于UT 算法主要是由計(jì)算解析解過(guò)程中的截?cái)嘁鸬?。相較于SUT法，ISUT法的采樣點(diǎn)修正過(guò)程彌補(bǔ)了一部分解析計(jì)算造成的精度降低。

圖5 UT、SUT和ISUT法的歸一化LAM隨時(shí)間變化情況Fig.5 Normalized LAM of UT method，SUT method and ISUT method over time

在對(duì)比分析UT、ISUT 和GMMISUT 法時(shí)，以UT法的LAM 值為標(biāo)準(zhǔn)進(jìn)行歸一化，則3 種算法的歸一化LAM 值隨時(shí)間變化情況如圖6 所示。可以看出，在第1 天內(nèi)，ISUT 法與GMMISUT 法效果相近，GMMISUT 法的精度略微高于ISUT 法。第2 天之后，GMMISUT 法的精度相較于ISUT 法有顯著提升，且逐漸優(yōu)于UT 法。這是由于GMM 法的引入使得ISUT 法的精度得以提升。基于GMM 方法的性能與每個(gè)組件的協(xié)方差有關(guān)，協(xié)方差越小，結(jié)果越好。這一結(jié)論具有直觀意義，因?yàn)楫?dāng)GMM 中每個(gè)分量的協(xié)方差趨于零時(shí)，每個(gè)分量就會(huì)變得越來(lái)越像脈沖函數(shù)。從而更輕松地近似各種PDF。結(jié)合圖4 和圖6 可知，GMMISUT 法在計(jì)算效率優(yōu)于UT 法的同時(shí)，在大部分時(shí)段給出了優(yōu)于UT法的結(jié)果。

圖6 UT、ISUT和GMMISUT法的歸一化LAM變化情況Fig.6 Normalized LAM of UT method，ISUT method and GMMISUT method over time

4 結(jié)論

對(duì)于軌道不確定性傳播問(wèn)題，提出基于半解析法和球形單邊采樣無(wú)跡變換的軌道不確定性快速估計(jì)方法。該方法的特色為采用球形單邊采樣無(wú)跡變換結(jié)合高斯攝動(dòng)方程及其解析解進(jìn)行采樣點(diǎn)的傳播和修正，從而實(shí)現(xiàn)軌道不確定性快速估計(jì)。其中解析解考慮了J2 攝動(dòng)和大氣阻力攝動(dòng)的解析形式，被用于實(shí)現(xiàn)樣本點(diǎn)的快速傳播，加快計(jì)算效率；球形單邊采樣在保證二階精度的同時(shí)進(jìn)一步降低計(jì)算量；基于高斯攝動(dòng)方程和采樣點(diǎn)特性的采樣點(diǎn)修正緩解了由采用解析解進(jìn)行樣本點(diǎn)傳播導(dǎo)致的精度降低問(wèn)題。將該方法與GMM 結(jié)合能夠進(jìn)一步提高精度。仿真分析表明，所提算法在保證合理精度的同時(shí)，大幅提升了不確定性估計(jì)的效率，適用于在軌不確定性快速傳播，可為在軌避障分析等提供初始輸入。