考慮有色噪聲影響的脈沖星導航2級強跟蹤差分濾波器

許強，崔洪亮，丁邦平，趙愛罡，趙陽

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

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

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

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

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

1 X射線脈沖星導航原理

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

圖1 X射線脈沖星導航原理Fig. 1 X-ray pulsar navigation principle

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

根據幾何關系，tSSB和tSC之間的時間差Δt與rsat滿足：

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

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

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

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

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

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

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

2 誤差分析

2.1 星歷誤差

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

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

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

根據導航原理和文獻［7］的數據分析可以判斷，能夠對導航結果產生顯著影響的主要為中心天體的星歷位置誤差。因此，本文的星歷誤差分析均以中心天體的星歷位置誤差為對象，所涉及的仿真實驗均為應用場景最多的近地衛星，星歷計算和誤差設置也均參考DE421。選擇地球以外的其他中心天體或參考星歷對本文分析結論并無實質影響。

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

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

若導航系統使用帶誤差的位置矢量E進行導航解算，便會得到帶有誤差的觀測量Δ：

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

由于sat滿足：

則此時真實的時間差Δt應當滿足［14］：

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

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

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

表1 衛星軌道參數Table 1 Parameters of satellite orbit

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

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

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

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

表3 其他參數Table 3 Other parameters

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

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

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

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

圖4 導航系統速度均方根誤差Fig. 4 Velocity RMSE of navigation system

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

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

圖6 星歷誤差導致的模型誤差Fig. 6 Model errors caused by ephemeris error

2.2 有色噪聲

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

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

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

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

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

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

3 TSTDKF算法設計

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

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

3.1 差分無偏濾波器

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

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

根據以上關系可構造新的觀測量Z′k如下所示：

顯然此時新構建的觀測方程中除HkXk、ξk外均為已知量，而觀測噪聲也變為了白噪聲ξk。但是結合式（22）可見，此時觀測噪聲與過程噪聲仍然存在相關關系。因此，根據式（6）、式（22）及式（24）構建如下方程：

若要求過程噪聲與觀測噪聲不相關，則應當要求變量Ωk滿足：

即

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

因此，新構建的方程式（25）即可認為是差分無偏濾波器的狀態方程。結合式（24），可將新的差分無偏濾波器狀態空間方程表述為

式中：Φ′k|k?1為新的狀態轉移矩陣；W′k為過程噪聲，且與觀測噪聲ξk無關；Ψk為已知輸入量。

顯然，經過重構得到的狀態空間方程滿足觀測噪聲ξk為高斯白噪聲，過程噪聲W′k與觀測噪聲ξk無關的條件。但此時的過程噪聲W′k仍為有色噪聲。由于僅當過程噪聲為有色噪聲時，對卡爾曼濾波的影響僅體現在狀態量一步預測方差PXk|k?1的計算中［28-29］，因此下面結合以上關系推導PXk|k?1新的迭代方程。

定義：

則此時狀態量預測方差PXk|k?1應滿足：

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

而根據卡爾曼濾波的迭代關系，重構后的最優估計值應滿足如下關系：

即

進而求得協方差PXk?1，W′k為

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

式中：

其中：PVk?2，Vk?1、PQk?1，Qk為對應變量協方差。

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

式中：PVk?1、PQk為對應變量方差。

而根據一階自回歸模型的相關性質可得

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

結合以上推導及普通卡爾曼濾波的相關知識，可將差分無偏濾波器的迭代方程總結為

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

3.2 強跟蹤獨立偏差濾波器

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

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

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

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

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

4 仿真分析

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

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

σξ由式（19）確定，σe在2.2節數據分析基礎上，綜合考慮其他參數的影響設置為

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

其他相關變量的初始值設置為

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

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

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

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

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

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

表4 實驗條件設置情況Table 4 Setting of experimental conditions

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

表5 實驗結果統計Table 5 Statistics of experimental results

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

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

5 結論

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

2）在TKF濾波算法基礎上，提出了TSTDKF算法。該算法將無偏濾波器設計為一種改進的差分濾波器以提高其應對有色噪聲時的魯棒性，并在獨立偏差濾波器中結合觀測殘差構建多重自適應調節因子以提高算法的跟蹤性能。

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