考慮有色噪聲影響的脈沖星導(dǎo)航2級強(qiáng)跟蹤差分濾波器

許強(qiáng)，崔洪亮，丁邦平，趙愛罡，趙陽

青州高新技術(shù)研究所 測試控制系，濰坊 262500

X射線脈沖星是一種能夠周期性發(fā)射穩(wěn)定X射線脈沖信號的中子星［1］。利用這一脈沖信號實現(xiàn)定位、授時功能的導(dǎo)航方法稱為X射線脈沖星導(dǎo)航。它相比于傳統(tǒng)的慣性導(dǎo)航、無線電導(dǎo)航和星光導(dǎo)航等導(dǎo)航方式具有累積誤差較小、潛在應(yīng)用情景廣泛、信號不易受地面干擾源和天體運動影響、兼具授時功能等優(yōu)點。目前已經(jīng)有大量的國內(nèi)外學(xué)者研究了其在近地衛(wèi)星導(dǎo)航任務(wù)和地月軌道轉(zhuǎn)移、火星探測等行星際探測任務(wù)中的應(yīng)用問題［2-4］。

根據(jù)相關(guān)原理，X射線脈沖星導(dǎo)航需要使用到包括中心天體在內(nèi)的眾多太陽系天體的星歷數(shù)據(jù)［5］。而以近地衛(wèi)星為例，目前已發(fā)表的相關(guān)文獻(xiàn)證明1 km的地球軌道誤差就會對導(dǎo)航結(jié)果帶來將近千米級的定位誤差［6-7］。所以X射線脈沖星導(dǎo)航技術(shù)的應(yīng)用對太陽系內(nèi)天體的星歷數(shù)據(jù)精度提出了較高要求。此外，目前絕大多數(shù)的研究文獻(xiàn)均將導(dǎo)航濾波器模型中的過程噪聲與觀測噪聲認(rèn)為是高斯白噪聲［8-11］。過程噪聲主要為除中心天體外其他天體對航天器的攝動影響，而觀測噪聲為探測器對X射線光子的探測噪聲。實際上，這2部分噪聲在天體運動等各方面因素的影響下不僅很難保證為白噪聲，而且兩者間還可能存在相關(guān)的情況。如果采用基于白噪聲假設(shè)的數(shù)據(jù)模型進(jìn)行建模解算，那么同樣會影響計算精度。

目前，太陽系內(nèi)行星星歷的計算主要依賴于甚長基線干涉測量技術(shù)（Very Long Baseline Interferometry， VLBI）技術(shù)和不同空間任務(wù)返回的觀測數(shù)據(jù)，測定精度著實有限。以目前應(yīng)用較為廣泛的由美國噴氣推進(jìn)實驗室發(fā)布的DE（De?velopment Ephemerides）星歷為例，地球的位置精度只能精確到千米以內(nèi)，木星更是僅能到幾十千米［12-13］。因此，如何降低太陽系內(nèi)星歷誤差對X射線脈沖星導(dǎo)航的影響已經(jīng)成為一項亟待解決的實際問題。有學(xué)者研究了增廣估計算法以排除星歷誤差的干擾［7］。但是該方法會導(dǎo)致系統(tǒng)的狀態(tài)量維數(shù)增加，帶來較重的運算負(fù)擔(dān)和數(shù)值病態(tài)問題。而且，增廣估計算法將太陽系內(nèi)星歷誤差認(rèn)為是不變的?？紤]到實際星歷誤差會存在緩慢的變化，增廣估計算法長時間運行時的跟蹤性能無法保證。其他學(xué)者還研究了不同形式的差分導(dǎo)航算法［14-16］。但是這些導(dǎo)航算法大多沒有考慮差分計算過程中會導(dǎo)致過程噪聲與觀測噪聲相關(guān)的問題。雖然文獻(xiàn)［14］考慮了該問題，但它與其他研究文獻(xiàn)一樣，對原始模型中的過程噪聲和觀測噪聲都非常理想地認(rèn)為是高斯白噪聲。而文獻(xiàn)［16］采用的方法是多航天器間觀測數(shù)據(jù)的差分計算。該方法雖然具有一定效果，但并不適合單個航天器的導(dǎo)航系統(tǒng)。

本文針對以上問題提出了適合繞某一中心天體軌道運行的2級強(qiáng)跟蹤差分卡爾曼濾波器（Two-stage Strong Tracking Differential Kal?man Filter，TSTDKF）。首先將變化較為緩慢的中心天體星歷誤差作為獨立偏差濾波器的狀態(tài)量予以實時估計?？紤]到星歷的周年變化，構(gòu)建了多重自適應(yīng)調(diào)節(jié)因子以保證算法的跟蹤精度。為解決原始模型中可能存在的有色噪聲問題，將第1級濾波器設(shè)計為一種改進(jìn)的差分卡爾曼濾波器，提高其對有色噪聲的魯棒性能，拓寬X射線脈沖星的應(yīng)用范圍。

1 X射線脈沖星導(dǎo)航原理

以繞地球運行的航天器為例，X射線脈沖星導(dǎo)航的原理可描述為如圖 1所示。其中，c為光速；rsat為衛(wèi)星在國際天球參考坐標(biāo)系（Interna?tional Celestial Reference System， ICRS）［17］中的位置矢量；rE和νE分別為地球在ICRS坐標(biāo)系中的位置和速度矢量；r和ν分別為航天器在地球質(zhì)心坐標(biāo)系（Earth Centered Inertial， ECI）中的位置矢量和速度矢量，且滿足rsat=rE+r；tSSB和tSC分別為太陽系質(zhì)心SSB（Solar System Bary?center，SSB）和航天器處的脈沖信號到達(dá)時間TOA（Time-of-Arrival，TOA）；n為 脈 沖 星 在ICRS坐標(biāo)系中的單位方向矢量。

圖1 X射線脈沖星導(dǎo)航原理Fig. 1 X-ray pulsar navigation principle

安裝在航天器上的X射線探測器配合相應(yīng)的信號恢復(fù)算法，可以獲得tSC。而tSSB可以根據(jù)對應(yīng)脈沖星的相位時間模型進(jìn)行預(yù)測［18］。tSC與tSSB之間的差值Δt便是脈沖信號在航天器與SSB之間沿脈沖星方向傳播的時間，也就等價于兩點間距離在該脈沖星方向的投影。如果選用3顆不同方向的脈沖星，則通過坐標(biāo)變換即可得到航天器在ICRS坐標(biāo)系中的絕對位置。配合卡爾曼濾波等迭代算法即可獲得穩(wěn)定持續(xù)的導(dǎo)航信息。

根據(jù)幾何關(guān)系，tSSB和tSC之間的時間差Δt與rsat滿足：

考慮到脈沖信號在傳播過程中會受到引力延遲效應(yīng)和周年視差效應(yīng)等因素的影響，其高階模型實際上可近似表示為［5］

式中：D為脈沖星到SSB的距離；b是SSB在ICRS坐標(biāo)系中相對于太陽質(zhì)心的位置矢量；b和rsat分別為b和rsat的長度模量；G是萬有引力常量；Msun是太陽質(zhì)量。

式（2）便是目前TOA計算應(yīng)用最為廣泛的一種高階模型，也是脈沖星導(dǎo)航卡爾曼濾波器中的觀測模型。而狀態(tài)模型則可用航天器繞中心天體飛行所滿足的動力學(xué)方程。對近地航天器而言，其在ECI坐標(biāo)系中滿足的動力學(xué)方程為

式 中：r=[rx，ry，rz]T；v=[vx，vy，vz]T；μ=GMEarth是地球引力常數(shù)；MEarth為地球質(zhì)量；J2為地球的二階帶諧項；REarth為地球半徑。ΔFx、ΔFy、ΔFz為其他天體對航天器的引力攝動，是過程噪聲的主要組成部分，f(?)為式（3）中的函數(shù)關(guān)系。若將系統(tǒng)的狀態(tài)量x選為航天器在ECI坐標(biāo)系中的位置r與速度v，觀測量確定為tSSB和tSC之間的時間差Δt，則系統(tǒng)的狀態(tài)空間方程為

式中：w為系統(tǒng)的過程噪聲，除其他天體的引力擾動外，還包括太陽光壓，地球潮汐等攝動影響；x=[rkvk]T為狀態(tài)量；z為觀測量，即不同脈沖星觀測的時間差Δt；h(?)為式（2）等號右側(cè)的函數(shù)關(guān)系式；V為觀測噪聲。

為滿足線性濾波器的應(yīng)用，常將式（4）利用雅克比矩陣進(jìn)行離散并線性化，將式（5）采取近似省略和泰勒展開方式進(jìn)行線性化。離散線性化后的狀態(tài)空間方程可表示為［19-20］

式中：Φk|k?1、Gk、Xk、Wk、Zk、Vk分別為離散化后的狀態(tài)轉(zhuǎn)移矩陣、噪聲驅(qū)動矩陣、狀態(tài)量、過程噪聲、觀測量及觀測噪聲，且Xk=[rkvk]T=[rxk，ryk，rzk，vxk，vyk，vzk]T；Δti為脈沖星i在航天器與SSB之間的時間差；rsatk?1為離散化后ICRS坐標(biāo)系中航天器在k-1時刻的位置矢量，rsatk?1為其標(biāo)量形式；Ek=[rEkvEk]T為離散化后ICRS坐標(biāo)系中地球的位置及速度矢量；rk、vk為離散化后ECI坐標(biāo)系中航天器的位置及速度矢量，且滿足rsatk=rEk+rk；ni與Di分別為不同脈沖星的單位方向矢量和到SSB的距離；Hk為觀測矩陣。

2 誤差分析

2.1 星歷誤差

在國際地球自轉(zhuǎn)服務(wù)組織 （International Earth Rotation and Reference Systems Service，IERS）2003、2010規(guī) 范 中，分 別 將DE405和DE421歷表作為ICRS坐標(biāo)系中的動力實現(xiàn)［21-22］。受不同時期觀測技術(shù)的限制，DE系列星歷的精度也是在逐漸提升的。2008年發(fā)布的DE421相對于1995年的DE405而言，其定位精度就得到了極大地改善［23］。雖然無法獲得天體真實的位置數(shù)據(jù)，但是通過不同星歷數(shù)據(jù)間的比較作差也能較為直觀地展示出精度上的提升。以地球為例，分別利用DE405和DE421計算其從2015-10-1704：00：00到2016-10-1604：00：00在ICRS坐標(biāo)系中的位置變化情況，采樣點間隔設(shè)置為5 h。2種星歷計算結(jié)果間的差值如圖 2所示。

通過圖 2，基本上可認(rèn)為DE421相對于DE405軌道誤差降低了數(shù)十千米。但是DE421仍然是不完美的。其中金星、地球和火星的軌道精度僅為亞千米級，而水星的軌道誤差約為數(shù)千米，木星和土星更是達(dá)到數(shù)十千米［12］。

圖2 DE405與DE421間地球星歷位置誤差Fig. 2 Earth position error between DE405 and DE421

根據(jù)導(dǎo)航原理和文獻(xiàn)［7］的數(shù)據(jù)分析可以判斷，能夠?qū)?dǎo)航結(jié)果產(chǎn)生顯著影響的主要為中心天體的星歷位置誤差。因此，本文的星歷誤差分析均以中心天體的星歷位置誤差為對象，所涉及的仿真實驗均為應(yīng)用場景最多的近地衛(wèi)星，星歷計算和誤差設(shè)置也均參考DE421。選擇地球以外的其他中心天體或參考星歷對本文分析結(jié)論并無實質(zhì)影響。

假設(shè)地球真實的位置矢量rE與根據(jù)星歷計算得到的位置矢量E滿足：

式中：δrE為地球星歷誤差。

若導(dǎo)航系統(tǒng)使用帶誤差的位置矢量E進(jìn)行導(dǎo)航解算，便會得到帶有誤差的觀測量Δ：

式中：Δ為帶有模型誤差的時間差。

由于sat滿足：

則此時真實的時間差Δt應(yīng)當(dāng)滿足［14］：

式中：B為地球星歷位置誤差導(dǎo)致的模型誤差；Γk為地球星歷誤差的觀測矩陣。

目前眾多已發(fā)表的文獻(xiàn)證明，在X射線脈沖星導(dǎo)航的數(shù)百乃至上千秒的單個觀測周期內(nèi)可認(rèn)為天體的星歷誤差δrE和對應(yīng)的模型誤差B是不變的［7，14-16］。但是長時間來看，這部分誤差卻因天體運動而存在緩慢且大幅的變化，存在導(dǎo)致算法發(fā)散的隱患。

為驗證以上分析，利用STK軟件產(chǎn)生一條衛(wèi)星仿真軌道，參數(shù)如表 1所示。假設(shè)該軌道上運行的某一航天器裝有X射線脈沖星導(dǎo)航系統(tǒng)，且其觀測的脈沖星為B0531+21、B1821-24、B1937+21，具體參數(shù)如表2所示［24-25］。

表1 衛(wèi)星軌道參數(shù)Table 1 Parameters of satellite orbit

表2 脈沖星參數(shù)［24-25］Table 2 Parameters of pulsars［24-25］

由于DE421中地球軌道精度為亞千米級，且大小變化受天體周期運動影響大多呈現(xiàn)類似的三角函數(shù)關(guān)系。因此，本文的數(shù)據(jù)仿真中假設(shè)ICRS坐標(biāo)系中地球的星歷位置誤差變化滿足：

式中：torbit為地球距離2015-10-1704：00：00的運行時間；Tearth為地球公轉(zhuǎn)一周的時間。

導(dǎo)航使用擴(kuò)展卡爾曼濾波器（Extended Kal?man Filter，EKF），觀測周期tobs為300 s，探測器面積等其他導(dǎo)航參數(shù)見表 3。初始誤差δX0取為［1 km， 1 km， 1 km， 0.02 m/s， 0.02 m/s， 0.02 m/s］。

表3 其他參數(shù)Table 3 Other parameters

脈沖星的觀測噪聲方差σV可確定［5］為

式中：Sdec為有效觀測面積；d為Pwidth與Ptime之比。

仿真過程中導(dǎo)航系統(tǒng)的均方根誤差（Root Mean Squared Error，RMSE）、地球星歷誤差、模型誤差c·B變化情況如圖 3～圖 6所示。

圖3 位置均方根誤差Fig. 3 Position RMSE of navigation system

圖4 導(dǎo)航系統(tǒng)速度均方根誤差Fig. 4 Velocity RMSE of navigation system

表3、圖 3～圖 6驗證了理論分析結(jié)果，證明對中心天體星歷誤差進(jìn)行實時估計時，有必要考慮算法的跟蹤性能，以解決星歷誤差緩慢變化的影響。

圖5 地球星歷誤差Fig. 5 Earth ephemeris error

圖6 星歷誤差導(dǎo)致的模型誤差Fig. 6 Model errors caused by ephemeris error

2.2 有色噪聲

式（4）所示的狀態(tài)方程中，將地球以外的第三天體引力ΔFx、ΔFy、ΔFz作為白噪聲計入過程噪聲，實際環(huán)境中這部分噪聲是難以保證為白噪聲的。同樣以表1中的衛(wèi)星軌道為例，利用DE421星歷計算衛(wèi)星受到的第三天體引力。第三天體引力模型為［5］

式中：μ3rd?body為不同天體的引力常量；rsat3rd?body為航天器相對第三天體的位置矢量；r3rd?bodyEarth為第三天體相對地球的位置矢量；rsat3rd?body和r3rd?bodyEarth分別為對應(yīng)矢量的長度模量。

所考慮的第三天體包括：太陽、水星、金星、火星、木星、土星、天王星、海王星、冥王星及月球。ICRS坐標(biāo)系中不同軸向的引力攝動及頻譜分析結(jié)果如圖 7和圖 8所示。顯然，該仿真軌道上的第三天體引力攝動都不符合高斯白噪聲的特點。

圖7 不同軸向引力攝動Fig. 7 Gravitational perturbation in different axes

圖8 引力攝動頻譜分析Fig. 8 Spectrum analysis of gravitational perturbation

對于觀測噪聲而言，X射線探測器在工作時受自身精度及電流擾動等因素影響，也會產(chǎn)生設(shè)備噪聲，該部分噪聲往往是有色噪聲。除此之外，式（2）的高階模型實際上也是省略了部分高階項的近似非線性觀測模型。省略的高階項為太陽系內(nèi)除太陽與地球外其他天體運動對脈沖星觀測的影響。結(jié)合以上觀測噪聲的分析，這部分噪聲極有可能也為有色噪聲。且由于都與第三天體運動有關(guān)，還可能出現(xiàn)過程噪聲與觀測噪聲相關(guān)的情況。因此，為保證算法的可靠性，有必要針對普通卡爾曼濾波器進(jìn)行改進(jìn)設(shè)計。

3 TSTDKF算法設(shè)計

2級卡爾曼濾波器（Two-stage Kalman fil?ter，TKF）的設(shè)計初衷便是代替增廣算法解決線性系統(tǒng)中的定常偏差問題［26］。它具有收斂速度快、運算量低等優(yōu)點，也適用于對普通卡爾曼濾波器的優(yōu)化改進(jìn)［27］。TKF由2部分組成，第1級為無偏濾波器，實際上是不考慮偏差影響的普通卡爾曼濾波器，第2級為獨立偏差濾波器，單獨用于實現(xiàn)偏差量的估計及結(jié)果補償。

根據(jù)式（11）～式（16）的分析，因Γk作為觀測矩陣可根據(jù)天體參數(shù)實時求解，則可采用2級濾波方法將地球星歷誤差δrE作為獨立偏差濾波器的狀態(tài)量予以估計。但由于δrE存在緩慢變化的性質(zhì)，因此在TKF基礎(chǔ)上本文設(shè)計了TSTDKF濾波器。該濾波器將無偏濾波器改進(jìn)為一種考慮相關(guān)噪聲的差分無偏濾波器，并在獨立偏差濾波器中添加了多重自適應(yīng)調(diào)節(jié)因子以增強(qiáng)其對緩慢變化量的跟蹤性能，解決傳統(tǒng)TKF只對定常偏差有效的局限性。

3.1 差分無偏濾波器

假設(shè)式（6）和式（7）中的過程噪聲Wk與觀測噪聲Vk為滿足以下一階自回歸模型的有色噪聲：

式中：ξk、ek為互不相關(guān)的高斯白噪聲，方差分別為σξ、σe；L、J、F均為對角相關(guān)系數(shù)矩陣，表征有色噪聲的相關(guān)性強(qiáng)弱；假設(shè)此相關(guān)系數(shù)矩陣均已通過參數(shù)辨識等手段獲得；當(dāng)相關(guān)系數(shù)矩陣均為零矩陣時說明過程噪聲與觀測噪聲為白噪聲。

根據(jù)以上關(guān)系可構(gòu)造新的觀測量Z′k如下所示：

顯然此時新構(gòu)建的觀測方程中除HkXk、ξk外均為已知量，而觀測噪聲也變?yōu)榱税自肼暒蝛。但是結(jié)合式（22）可見，此時觀測噪聲與過程噪聲仍然存在相關(guān)關(guān)系。因此，根據(jù)式（6）、式（22）及式（24）構(gòu)建如下方程：

若要求過程噪聲與觀測噪聲不相關(guān)，則應(yīng)當(dāng)要求變量Ωk滿足：

即

式中：I為單位矩陣；HkGkJ+I需滿足為非奇異矩陣。

因此，新構(gòu)建的方程式（25）即可認(rèn)為是差分無偏濾波器的狀態(tài)方程。結(jié)合式（24），可將新的差分無偏濾波器狀態(tài)空間方程表述為

式中：Φ′k|k?1為新的狀態(tài)轉(zhuǎn)移矩陣；W′k為過程噪聲，且與觀測噪聲ξk無關(guān)；Ψk為已知輸入量。

顯然，經(jīng)過重構(gòu)得到的狀態(tài)空間方程滿足觀測噪聲ξk為高斯白噪聲，過程噪聲W′k與觀測噪聲ξk無關(guān)的條件。但此時的過程噪聲W′k仍為有色噪聲。由于僅當(dāng)過程噪聲為有色噪聲時，對卡爾曼濾波的影響僅體現(xiàn)在狀態(tài)量一步預(yù)測方差PXk|k?1的計算中［28-29］，因此下面結(jié)合以上關(guān)系推導(dǎo)PXk|k?1新的迭代方程。

定義：

則此時狀態(tài)量預(yù)測方差PXk|k?1應(yīng)滿足：

式中：E（·）為求期望運算。

而根據(jù)卡爾曼濾波的迭代關(guān)系，重構(gòu)后的最優(yōu)估計值應(yīng)滿足如下關(guān)系：

即

進(jìn)而求得協(xié)方差PXk?1，W′k為

結(jié)合式（21）、式（26）及式（32），可將式（38）整理為

式中：

其中：PVk?2，Vk?1、PQk?1，Qk為對應(yīng)變量協(xié)方差。

同理，結(jié)合式（32）可求得PW′k滿足：

式中：PVk?1、PQk為對應(yīng)變量方差。

而根據(jù)一階自回歸模型的相關(guān)性質(zhì)可得

將式（38）～式（46）結(jié)論代入式（35）即可求得PXk|k?1。

結(jié)合以上推導(dǎo)及普通卡爾曼濾波的相關(guān)知識，可將差分無偏濾波器的迭代方程總結(jié)為

顯然，根據(jù)Φ′k|k?1、Ψk及W′k的定義，當(dāng)觀測噪聲與過程噪聲為不相關(guān)的白噪聲，即L、J、F陣均為零陣時，以上過程退化為EKF。

3.2 強(qiáng)跟蹤獨立偏差濾波器

結(jié)合式（15）、式（28）和式（29），若普通獨立偏差濾波器的狀態(tài)量bk確定為地球星歷位置誤差δrE，則迭代方程可表述為［30］

式中：Λk為獨立偏差濾波器對無偏濾波器狀態(tài)量的糾正矩陣；Sk和Uk為該級濾波器迭代的中間量；Mk相當(dāng)于該級濾波器狀態(tài)量的協(xié)方差；為該級濾波器增益矩陣。

通過分析式（57）可以發(fā)現(xiàn)，獨立偏差濾波器的狀態(tài)更新仍需用到無偏濾波器的殘差信息。對于TKF而言，無偏濾波器的殘差變化一定程度體現(xiàn)了獨立偏差濾波器狀態(tài)量bk的變化趨勢。因為除噪聲因素外，無偏濾波器的殘差主要由bk，即系統(tǒng)偏差決定。因此，本文采用無偏濾波器的殘差方差作為強(qiáng)化獨立偏差濾波器跟蹤性能的依據(jù)，在增益矩陣的計算中構(gòu)建多重自適應(yīng)調(diào)節(jié)因子λk如式（59）～式（61）所示：

式中：ρ=[ρ1，ρ2，ρ3]為尺度調(diào)節(jié)參數(shù)，可由先驗知識確定，無先驗知識情況下ρ1、ρ2、ρ3可設(shè)置為1；ΔZ′k為根據(jù)狀態(tài)量預(yù)測值所確定的觀測殘差；τ為記憶參數(shù)，用于確定λk計算的漸消記憶窗口長度；(?)i，i表示矩陣的第i個對角線元素。

實際上，尺度調(diào)節(jié)參數(shù)ρ是調(diào)整在線跟蹤性能的指標(biāo)參數(shù)。尺度調(diào)節(jié)參數(shù)越大，跟蹤性越強(qiáng)。在已經(jīng)出現(xiàn)大范圍誤差的情況下可以更快速的完成誤差糾正，但是在已經(jīng)接近實際結(jié)果的穩(wěn)態(tài)下，該參數(shù)越大，越會增加穩(wěn)態(tài)下估計結(jié)果的波動性。因此，該參數(shù)可以在有先驗知識情況下進(jìn)行調(diào)整，例如可定期對使用脈沖星導(dǎo)航的衛(wèi)星軌道進(jìn)行測定，以確定不同階段尺度調(diào)節(jié)參數(shù)的大小，進(jìn)而遠(yuǎn)程進(jìn)行參數(shù)注入。若偏離誤差較大，可適當(dāng)增加尺度調(diào)節(jié)參數(shù)大小；若已達(dá)到誤差允許范圍內(nèi)的穩(wěn)態(tài)，可適當(dāng)減小該值。取值范圍一般在（0，50］區(qū)間即可。

4 仿真分析

為驗證算法的有效性，使用表 1中參數(shù)所確定的衛(wèi)星軌道進(jìn)行仿真分析，橫向比較有色噪聲條件下，存在地球星歷誤差時EKF、TKF及TSTDKF的導(dǎo)航性能。3種算法的觀測周期為300 s，位置及速度誤差采用均方根誤差計算，其他參數(shù)仍同表2和表3。對于TSTDFK而言，ρ=[8， 15， 9]，τ=100。

假設(shè)觀測噪聲與過程噪聲滿足式（21）～式（23）所示的一階自回歸模型，具體參數(shù)取值如下所示：

σξ由式（19）確定，σe在2.2節(jié)數(shù)據(jù)分析基礎(chǔ)上，綜合考慮其他參數(shù)的影響設(shè)置為

式中：e1=(10?8)2；e2=(10?10)2。

其他相關(guān)變量的初始值設(shè)置為

當(dāng)星歷誤差變化滿足式（17）關(guān)系時，3種算法的仿真結(jié)果對比如圖 9所示。通過對圖 9的分析可以發(fā)現(xiàn)，在噪聲相關(guān)且存在變化星歷誤差條件下，TSTDKF的位置及速度均方根誤差均低于TKF和EKF。在星歷誤差估計方面，TSTDKF也比TKF收斂更快，效果更好。

圖9 不同算法的仿真結(jié)果Fig. 9 Simulation results of different algorithms

為進(jìn)一步驗證算法對地球星歷誤差的跟蹤性能，在相同參數(shù)條件下將觀測周期延長為1 a，并將星歷誤差變化趨勢在式（17）基礎(chǔ)上，增加式（69）、式（70）2種趨勢，觀察TKF與TSTDKF的跟蹤情況，如圖 10所示。

圖10 不同算法對星歷誤差的跟蹤情況Fig. 10 Tracking accuracy of ephemeris errors by different algorithms

結(jié)合圖10，雖然長期來看TSTDKF也會在某些拐點存在一定的失真，但均能夠有效糾正，總體精度較為穩(wěn)定，說明本文所設(shè)計的TSTDKF在應(yīng)對變化的星歷誤差方面是有效的。

為說明TSTDKF在應(yīng)對不同噪聲情況下的魯棒性，根據(jù)表 4中的實驗設(shè)置分別對3種算法各進(jìn)行50次蒙特卡洛仿真，時間30 d。星歷誤差仍根據(jù)式（17）設(shè)置，最終不同導(dǎo)航算法的均方根誤差統(tǒng)計如表 5所示。

表4 實驗條件設(shè)置情況Table 4 Setting of experimental conditions

通過對表 5數(shù)據(jù)的分析可以直觀發(fā)現(xiàn)：不同噪聲條件下，TSTDKF的導(dǎo)航精度均優(yōu)于其他算法，且性能表現(xiàn)相對穩(wěn)定；位置誤差對噪聲條件的敏感度比速度誤差高；在實驗1條件下，TSTDKF優(yōu)勢最為突出，位置和速度均方根誤差分別比EKF提升56.49%和27.66%，比TKF提升35.18%和17.07%；相比于過程噪聲，觀測噪聲為有色噪聲時對算法影響更大；在實驗5條件下，TKF結(jié)果與TSTDKF結(jié)果相當(dāng)。

表5 實驗結(jié)果統(tǒng)計Table 5 Statistics of experimental results

之所以在實驗5條件下TKF和TSTDKF結(jié)果相當(dāng)，是因為各噪聲為白噪聲時，TSTDKF的無偏濾波器退化為與TKF完全一致。此時兩者的差別僅體現(xiàn)在獨立偏差濾波器的設(shè)計上。而根據(jù)圖9和圖10可見，TKF算法在開始的30 d內(nèi)雖然收斂較慢，但仍可緩慢地接近標(biāo)準(zhǔn)值。但是長時間運行時TKF往往難以保證較好的跟蹤精度，因此隨著時間增加，導(dǎo)航誤差將會逐漸增大。

這是因為TKF的設(shè)計初衷便是應(yīng)對定常偏差問題。雖然在其他文獻(xiàn)中也指出針對緩慢變化的誤差具有一定的有效性，但是這僅僅是短期及小范圍內(nèi)。如圖9所示，在初期TKF確實可以在一定程度上實現(xiàn)對星歷誤差的追蹤。這是因為地球星歷誤差變化緩慢且小幅，基本可認(rèn)為常值。但是隨著時間延遲、變化幅度增大，TKF缺乏有效跟蹤性能的問題就會暴露，如圖10所示。

5 結(jié)論

1）通過誤差分析，證明有色噪聲和中心天體的星歷誤差是制約X射線脈沖星導(dǎo)航精度的影響因素之一，有必要針對該誤差因素設(shè)計相應(yīng)的算法予以補償。

2）在TKF濾波算法基礎(chǔ)上，提出了TSTDKF算法。該算法將無偏濾波器設(shè)計為一種改進(jìn)的差分濾波器以提高其應(yīng)對有色噪聲時的魯棒性，并在獨立偏差濾波器中結(jié)合觀測殘差構(gòu)建多重自適應(yīng)調(diào)節(jié)因子以提高算法的跟蹤性能。

3）通過仿真實驗，證明存在有色噪聲和變化的地球星歷誤差時，TSTDKF算法的導(dǎo)航性能優(yōu)于EKF和TKF，位置和速度均方根誤差最大可比EKF提升56.49%和27.66%，比TKF提升35.18%和17.07%，對星歷誤差的估計精度也好于TKF。