羅楠 許錄平 張華

（西安電子科技大學(xué)電子工程學(xué)院，西安710071）

1 引言

脈沖星是一種高速、穩(wěn)定自轉(zhuǎn)的中子星，其輻射周期信號的長期穩(wěn)定度可與銫原子鐘媲美，被譽為自然界最穩(wěn)定的頻率基準[1]。美國科學(xué)家于1974年首次提出了利用脈沖星進行航天器自主軌道確定，開創(chuàng)了脈沖星導(dǎo)航先河[2]。將脈沖星作為航天器的導(dǎo)航信源有其獨特的優(yōu)點：脈沖星信號是自然存在的，任何時候都可以免費獲取；脈沖星分布于整個天空，不同的航天器可以同時進行時間差測量；每顆脈沖星輻射輪廓都是唯一并可預(yù)測的；能為近地及星際空間的多種航天器提供導(dǎo)航信息。基于脈沖星輻射定時模型的導(dǎo)航是一種頗具潛力的航天器自主導(dǎo)航方法，它不僅可以脫離地基測控的支持，為近地遙感及深空探測的各種衛(wèi)星提供長期穩(wěn)定可靠的自主導(dǎo)航服務(wù)，而且能為需要多衛(wèi)星數(shù)據(jù)融合的衛(wèi)星系統(tǒng)提供星座間位置與時間的聯(lián)合標定。但由于脈沖星信號弱，在目前技術(shù)條件下導(dǎo)航精度較低，因此組合不同類型的導(dǎo)航手段，用信息融合方法來完成精確、穩(wěn)定的導(dǎo)航任務(wù)是必要的。

本文在分析星敏感器和脈沖星導(dǎo)航兩種不同導(dǎo)航機制的基礎(chǔ)上，提出了一種基于無跡卡爾曼濾波（Unscented Kalman Filter，UKF）的信息融合天文自主導(dǎo)航方法。結(jié)合兩者在導(dǎo)航信息連續(xù)性和精確性上的互補優(yōu)勢，一方面用脈沖星提供連續(xù)不斷的導(dǎo)航信息；另一方面利用星敏感器敏感折射星光，并通過大氣對星光折射的數(shù)學(xué)模型和誤差補償技術(shù)精確敏感地平，提高這種組合導(dǎo)航方法的導(dǎo)航精度。

2 敏感地平測量模型

星光仰角是指航天器觀測的恒星星光矢量方向與地心矢量方向之間的夾角，它是航天器在地心坐標系中位置矢量的函數(shù)。如圖1所示，當(dāng)衛(wèi)星到地心與地球邊緣切線方向的夾角已知，星敏感器的量測模型[2]可表示為

式中α為航天器到地心方向與到恒星方向的夾角；rE為地心坐標系中衛(wèi)星位置矢量；為地心慣性坐標系中衛(wèi)星到導(dǎo)航恒星視線方向單位矢量；vα為量測噪聲，包括安裝誤差、地平敏感誤差等。

公式 （1）表示成量測方程形式為

式中Z為量測向量；h為非線性轉(zhuǎn)換方程；X（）t為衛(wèi)星的狀態(tài)向量；δZ為噪聲向量，主要包括安裝誤差和地平敏感器量測誤差。

圖1 星敏感器／脈沖星聯(lián)合導(dǎo)航示意圖Fig.1 Schematic diagram of integrated navigation utilizing star sensor and pulsar

3 脈沖星到達時間測量模型

脈沖星導(dǎo)航技術(shù)是利用探測器同時精確測量多顆特性已知的脈沖星輻射脈沖到達太陽系質(zhì)心（Solar System Barycentre，SSB）和航天器的時間差，經(jīng)過必要的修正測得航天器在脈沖星視線方向距SSB的距離，結(jié)合軌道動力學(xué)模型采用相應(yīng)的信號處理和導(dǎo)航算法實現(xiàn)對航天器位置、速度、時間的計算和預(yù)報。采用SSB慣性坐標系作為參考架，脈沖星脈沖到達SSB相位模型[3]可表示為

式中φ0為參考歷元t0時刻的脈沖相位；ν（t）為脈沖星參考歷元t時的自轉(zhuǎn)頻率。該模型中的參數(shù)通過對脈沖星的長期觀測得到。將脈沖星輻射短期觀測和累積得到的脈沖星脈沖輪廓模型，在星載原子鐘保持時間（固有時）下進行比對，可以得到脈沖星到達用戶星和SSB的時間差[4]。如圖1所示，由于時間差是航天器位置矢量的函數(shù)，SSB坐標系下其一階表達式可表示為

式中rssb＝rEssb＋rE；為SSB到脈沖星視線方向的矢量；tU，i為航天器測得的第i顆脈沖星輻射脈沖到達時間（Time of Arrival，TOA）；tb，i為同一脈沖在SSB處的TOA，i表示第i顆脈沖星。公式（4）的一階方程忽略了相對論等效應(yīng)對TOA測量的影響，要在工程上實現(xiàn)精確的脈沖星定時模型，必須將實測TOA的固有時通過各項修正轉(zhuǎn)化到質(zhì)心力學(xué)時（Barycentric Dynamical Time，TDB），求得TDB下的時間差：

式中ΔRΘ為Roemer延遲，包括一階多普勒延遲和視差效應(yīng)；ΔEΘ為相對論效應(yīng)修正，太陽系中ΔEΘ＜1ns可忽略；ΔSΘ為太陽系Shapiro延遲；Δf，i和δt分別為第i顆脈沖星色散效應(yīng)誤差和衛(wèi)星時鐘偏差。時間差測量方程可簡化為

式中D0為脈沖星到SSB的距離；V為脈沖星固有運動速度；ΔtN＝tN-t0表示從脈沖星初始時刻開始到輻射第N個脈沖的時間間隔，由于D0?VΔtN，仿真計算時VΔtN可以忽略；b表示SSB相對太陽質(zhì)心的位置矢量；μsun為太陽引力常數(shù)。

式（6）用測量方程簡化表示為

4 軌道動力學(xué)模型

定軌問題中軌道狀態(tài)常用開普勒軌道要素或三維慣性坐標系下的位置、速度來描述。開普勒軌道六參數(shù)有5個是常量，定軌精度高，適用于特定軌道；對于具有變軌任務(wù)的衛(wèi)星，使用位置、速度描述更為方便。地心坐標系下設(shè)向量X＝｛rT，vT｝T表示衛(wèi)星狀態(tài)，衛(wèi)星運動用二體方程描述：

式中GE為地球引力常數(shù)；a為軌道攝動力，a＝anonspher＋a3body＋adrag＋aSR＋aother，其中anonspher為地球非球形引力攝動，a3body為第三體攝動，adrag為大氣阻力攝動，aSR為太陽直接輻射壓攝動，aother為其他攝動力。地球質(zhì)心引力與各種攝動力共同決定了衛(wèi)星的在軌運動。

5 基于UKF的信息融合方法

定軌過程首先是對衛(wèi)星受力建模，并對運動方程積分，得到預(yù)測軌道狀態(tài)；然后，在量測量支持下，依賴軌道的估計采用擬合或濾波的方法進行差分改正或軌道狀態(tài)改進。從以上各節(jié)分析可知，基于式（2）、式（7）的量測方程和式（9）的軌道動力學(xué)方程可構(gòu)建迭代濾波對軌道狀態(tài)進行改進。在系統(tǒng)噪聲滿足高斯分布，統(tǒng)計特性的先驗知識足夠的前提下，Kalman濾波是基于H2范數(shù)導(dǎo)出的最優(yōu)線性迭代濾波方法。由于軌道改進過程的量測方程和狀態(tài)方程均是非線性的，通常使用次優(yōu)的擴展卡爾曼濾波（EKF），它使用泰勒級數(shù)的一階項對觀測方程和狀態(tài)方程線性化。EKF有兩個缺陷：對強非線性方程線性化會引入較大的誤差，導(dǎo)致估計誤差增大；線性化過程會導(dǎo)致假設(shè)為高斯分布的先驗和后驗概率的均值和方差估計錯誤。而無跡卡爾曼濾波（UKF）被認為是性能優(yōu)于EKF的非線性濾波方法[6]，它使用一組確定的離散Sigma采樣點來參數(shù)化高斯隨機狀態(tài)變量的均值和方差，并直接使用非線性狀態(tài)方程來傳播后驗概率均值和方差，而不需要線性化，但UKF方法穩(wěn)定性較EKF稍差。

文獻[7-8]分析了X射線脈沖星導(dǎo)航技術(shù)在航天器行星際空間飛行中的應(yīng)用，解決了相位導(dǎo)航存在的整周期模糊度問題。在軌道動力學(xué)模型估計航天器位置的基礎(chǔ)上，將其真實位置和估計位置的脈沖星脈沖相位差作為反饋，使用UKF濾波來得到高精度的導(dǎo)航信息。本文通過分析星敏感器和脈沖星導(dǎo)航的工作機制，將來自這兩種不同導(dǎo)航方式的信息源數(shù)據(jù)進行相關(guān)、組合，充分利用多個傳感器提供的導(dǎo)航資源，合理使用它們在空間和時間上的冗余互補信息，設(shè)計了UKF聯(lián)合濾波器來完成信息融合，以提高衛(wèi)星自主導(dǎo)航的精確度。

設(shè)UKF算法的離散狀態(tài)方程和測量方程為

式中xk為n維向量狀態(tài)變量；f（xk-1）為狀態(tài)傳播函數(shù)；h（xk）為測量函數(shù)；wk和vk為非相關(guān)高斯白噪聲，且分布滿足wk～N（0，Qk）和vk～N（0，Rk），Qk和Rk分別為wk和vk的方差。

設(shè)UKF的算法實現(xiàn)過程為

4）時間更新：

5）量測更新：

聯(lián)合濾波系統(tǒng)中，子濾波器1基于狀態(tài)方程式（9）與觀測方程式（2）；子濾波器2基于狀態(tài)方程式（9）與觀測方程式（7）。融合濾波器將兩個子濾波器的輸出結(jié)果進行信息融合，完成時間更新，并將融合結(jié)果反饋到各子濾波器，作為下一級濾波的初值。鑒于UKF的穩(wěn)定性較差，將一組EKF作為備份濾波器加入組合定姿系統(tǒng)中，用于在UKF異常時維持系統(tǒng)正常工作，導(dǎo)航信息融合方法如圖2所示。

圖2 信息融合方法Fig.2 Methodology of information fusion

由于UKF的穩(wěn)定性較EKF差，子濾波器中只反饋狀態(tài)變量。衛(wèi)星在軌運動由于受地球等天體遮擋影響，星敏感器觀測到的導(dǎo)航恒星和X射線探測器探測到的脈沖星個數(shù)是不斷變化的，因此兩個子濾波器在運算過程中觀測方程和觀測噪聲協(xié)方差陣也隨之變化。

6 仿真分析

6.1 仿真系統(tǒng)

為了驗證X射線脈沖星導(dǎo)航原理和UKF信息融合算法的可行性和有效性，利用激光光子輻射模擬脈沖星X射線光子輻射，設(shè)計了一種基于激光光量子探測的X射線脈沖星導(dǎo)航地面仿真系統(tǒng)，其硬件結(jié)構(gòu)如圖3所示。信號模擬計算機提取各參數(shù)庫中信號模型生成脈沖星信號模擬數(shù)據(jù)和星敏感器觀測數(shù)據(jù)；導(dǎo)航算法部分模擬星載計算機，用于X射線脈沖星信號數(shù)據(jù)和星敏感器數(shù)據(jù)處理并用于信息融合導(dǎo)航。時間保持單元通過現(xiàn)場可編程門陣列 （Field-programmable Gate Array，F(xiàn)PGA）綜合GPS定時接收機提供的秒脈沖信號和溫補晶振的時鐘信號，用于標定光子發(fā)射時間；光子發(fā)射單元由激光調(diào)制與驅(qū)動、670nm半導(dǎo)體激光器和光學(xué)發(fā)射天線構(gòu)成，用于將模擬數(shù)據(jù)轉(zhuǎn)換為激光光子；光子探測單元包括接收光學(xué)天線、衰減器、670nm帶通濾光器和光子探測與計數(shù)器；光學(xué)接收天線模擬X射線準直器，用于定向和噪聲初步抑制，其接收的光子經(jīng)過衰減器衰減后由670nm帶通濾光器濾除日光中其他光譜成分，然后由時間保持單元2進行時間標定后，送至導(dǎo)航算法單元。

該仿真系統(tǒng)激光發(fā)射功率與探測到的光子數(shù)之間的關(guān)系[9]為

圖3 仿真系統(tǒng)結(jié)構(gòu)框圖Fig.3 Frame diagram of the simulation system

式中L為光學(xué)接收天線直徑；d為激光器和探測器間距離；θ為激光光束發(fā)散角；ρ為可控衰減系數(shù)，用可調(diào)衰減器實現(xiàn)；ω為其他衰減；η為探測器效率；T為脈沖持續(xù)時間；h為普朗克常量；c為光速；λ為光波長；SI為閃爍因子，典型值為0.4～1.0。若設(shè)觀測到的光子數(shù)為λS，背景輻射噪聲為λn，激光發(fā)射功率的概率模型為[9]

式中 [·]為取整操作；Pt為激光發(fā)射功率。仿真時系統(tǒng)根據(jù)模擬的軌道位置和時間確定信號的發(fā)射時間，再根據(jù)不同脈沖星空間觀測信號強度和背景輻射噪聲用式（16）計算激光發(fā)射功率，并對激光進行模擬調(diào)制。

6.2 仿真條件

衛(wèi)星軌道動力學(xué)模型考慮地球非球形引力 （攝動函數(shù)取J2～J6項），日、月引力及大氣阻力等3種攝動力，積分過程使用四階龍格-庫塔數(shù)值積分方法，積分步長設(shè)定為10s。用戶星軌道初值從北美防空聯(lián)合司令部提供的兩行軌道根數(shù)集中選取，并使用簡化常規(guī)／深空擾動的近似解析解模型[10]程序讀取。用戶星初始軌道誤差為1.5km，并假設(shè)用戶星到SSB的脈沖星輻射脈沖周期模糊數(shù)已知，衛(wèi)星嚴格對地定向并且三軸姿態(tài)穩(wěn)定。

星敏感器與紅外地平儀測量1σ精度分別設(shè)為6″和0.05°，視場角5°×5°，并設(shè)觀測噪聲服從零均值高斯分布。恒星星歷取自第谷星表數(shù)據(jù)庫，仿真中濾波器從備擇恒星中選取較亮的10顆作為導(dǎo)航星。X射線脈沖星選用B1937＋21、B1821-24、B1509-58和B0531＋21，并設(shè)脈沖星自身運動速度為零。4顆脈沖星的信號積分時間與距離測量精度關(guān)系見表1，其特征屬性參數(shù)見表2，其中pf為脈沖輻射流量與平均輻射流量比，本文仿真中積分時間選500s。脈沖星TOA量測值仿真利用式（5）的簡化形式加上以測距精度為方差的零均值高斯噪聲完成[11]，TOA測量使用泰勒快速傅里葉變換（Taylor FFT）算法。

表1 脈沖星觀測時間對定位精度的影響Tab.1 Pulsar range measurement accuracy over 500s，1000sand 5000sobservation time

表2 用于定軌的X射線脈沖星參數(shù)Tab.2 Parameters of pulsars for orbit determination

仿真系統(tǒng)中半導(dǎo)體激光器功率5mW，中心波長670nm；光學(xué)發(fā)射天線透光率約80%，發(fā)射角（π／6）rad；光學(xué)接收天線接收面直徑5cm，透光率約70%；衰減器使用兩片衰減深度為10-3的可調(diào)衰減器構(gòu)成，帶通濾光片中心波長（670±10）nm，峰值透過率55%；光子探測與計數(shù)器有效光敏面直徑20μm，暗計數(shù)50個／s，670nm波長光子探測效率為10%，光學(xué)發(fā)射天線與光學(xué)接收天線相距3m，采樣間隔（1／1024）s。以脈沖星B1937＋21為例做信號累積試驗，累積500s得到累積脈沖輪廓與標準輪廓的比較如圖4所示。計算二者相關(guān)系數(shù)為0.996，可以說明該系統(tǒng)可以較好地模擬出脈沖星輻射信號。

6.3 仿真結(jié)果

選取LEO衛(wèi)星MEGSAT-1（衛(wèi)星號26546U）做定軌試驗。圖5是一段分別使用星敏感器、脈沖星和信息融合方法導(dǎo)航，濾波周期為500s的仿真結(jié)果。從圖5的仿真結(jié)果可以看出，濾波收斂后，融合方法收斂速度更快，導(dǎo)航精度整體優(yōu)于單獨使用脈沖星或星敏感器方法的精度。

圖4 脈沖星累積輪廓與標準輪廓Fig.4 Standard and accumulated profiles of pulsar

圖5 星敏感器、脈沖星及融合方法定軌位置和速度誤差Fig.5 Position and velocity errors of star-sensor，pulsars and information fusion method

在相同仿真條件下，表3給出了使用EKF和UKF方法進行信息融合時，3種導(dǎo)航方法位置、速度的估計誤差結(jié)果。由表3可知，基于UKF的信息融合方法位置估計精度比單獨使用脈沖星和星敏感器的導(dǎo)航方法分別提高了52.7%和43.6%，速度估計精度分別提高了82.2%和70.5%。

表3 脈沖星、星敏感器、融合方法進行位置、速度估計的誤差對照Tab.3 Estimation errors of position and velocity using pulsars，star-sensor and information fusion method respectively

從圖5和表3給出的仿真結(jié)果可得出結(jié)論：1）脈沖星導(dǎo)航由于脈沖星輻射信號弱的限制，估計精度要略低于星敏感器；2）結(jié)合脈沖星和星敏感器的融合導(dǎo)航方法估計精度顯然高于單一導(dǎo)航方法，尤其是使用UKF時，導(dǎo)航誤差更小，大幅度地提高了導(dǎo)航精度。

圖6給出了相同的仿真條件下用EKF和UKF進行信息融合時，位置和速度估計誤差的比較。從圖中可以看出UKF融合方法能得到比EKF融合方法更好的性能。結(jié)合表3給出的估計誤差結(jié)果，從圖6（a）和（b）中的曲線對比可以發(fā)現(xiàn)：1）UKF方法相比于EKF來說，收斂得更快，收斂后估計精度更高；2）基于EKF的信息融合方法位置和速度估計精度雖然比基于UKF的信息融合方法差，但仍然優(yōu)于僅使用星敏感器或脈沖星的導(dǎo)航方法。

圖6 UKF、EKF信息融合方法誤差Fig.6 Error of UKF and EKF information fusion

仿真試驗結(jié)果表明，無論使用那種濾波器，融合方法均提高了導(dǎo)航精度，是完成精確導(dǎo)航任務(wù)的更佳選擇。但由于UKF必須保證迭代過程中Pk-1的正定性，其濾波穩(wěn)定性卻不如EKF，導(dǎo)航系統(tǒng)可使用UKF作為主濾波器，EKF作為副濾波器，以提高信息融合系統(tǒng)的可靠性。

7 結(jié)束語

基于UKF的信息融合方法，有效地組合了基于星敏感器和基于脈沖星兩種不同的導(dǎo)航技術(shù)，性能明顯優(yōu)于僅使用星敏感器或脈沖星的導(dǎo)航方法，且優(yōu)于EKF融合方法。設(shè)計了利用激光光子模擬脈沖星X射線光子的半實物仿真系統(tǒng)，利用該系統(tǒng)驗證了融合方法的有效性。信息融合導(dǎo)航方法具有完全自主性，可作為星載導(dǎo)航系統(tǒng)的備份，提高系統(tǒng)可靠性并可集成于星載處理器，不僅能用于地球軌道衛(wèi)星導(dǎo)航，還能為深空探測和太陽系空間航天器提供導(dǎo)航定位服務(wù)。

[1]JOSEPH H，TAYLOR J H.Millisecond pulsars：nature′s most stable clocks [J].Proceedings of the IEEE，1991，79（7）：1054-1062.

[2]李建勛，柯熙政.基于脈沖星定時模型的自主導(dǎo)航定位方法 [J].中國科學(xué)（G輯：物理學(xué) 力學(xué) 天文學(xué)），2009，39 （2）：311-317.LI JIANXUN，KE XIZHENG.Study on autonomous navigation based on pulsar timing model[J].Science in China（Series G：Physics，Mechanics and Astronomy），2009，39（2）：311-317.

[3]DOWNS G S.Interplanetary navigation using pulsating radio sources[R].NASA，N74-34150，1974：1-2.

[4]帥平，陳忠貴.關(guān)于X射線脈沖星導(dǎo)航的軌道力學(xué)問題 [J].中國科學(xué)（E輯：技術(shù)科學(xué)），2009，29（3）：556-561.SHUAI PING，CHEN ZHONGGUI.Several problems of orbital mechanics about X-ray pulsar navigation[J].Science in China （Series E：Technological Sciences），2009，29（3）：556-561.

[5]SHEIKH S I，PINES D J.Recursive estimation of spacecraft position and velocity using X-ray pulsar time of arrival measurement[J].Navigation，Journal of the Institute of Navigation，2006，53（3）：149-166.

[6]WAN E A，VAN DER MERWE R.The unscented Kalman filter for nonlinear estimation [C].Adaptive Systems for Signal Processing，Communications，and Control Symposium 2000，AS-SPCC，the IEEE 2000，2000：153-158.

[7]楊博，郭星燦，楊勇.X射線脈沖星導(dǎo)航在行星際軌道上的應(yīng)用 [J].北京航空航天大學(xué)學(xué)報，2009，35（11）：1384-1387.YANG BO，GUO XINGCAN，YANG YONG.Use of X-ray pulsar-based navigation method on interplanetary trajectory [J].Journal of Beijing University of Aeronautics and Astronautics，2009，35（11）：1384-1387.

[8]YANG BO，GUO XINGCAN，YANG YONG.The use of X-ray pulsar-based navigation method for interplanetary flight[C].International Symposium on Photoelectronic Detection and Imaging 2009：Terahertz and high energy radiation detection technologies and applications，Beijing，2009：17-19.

[9]張華，許錄平.基于光量子探測的XPNAV半物理仿真及建模 [J].光電子·激光，2011，22（6）：905-910.ZHANG HUA，XU LUPING.Hardware in-the-loop simulation and mode ling for XPNAV based on light quantum detection [J].Journal of Optoelectronics·Laser，2011，22 （6）：905-910.

[10]HOOTS F R，ROEHRICH R L.Spacetrack report No.3models for propagation of NORAD element sets[R].Department of Defense，Defense Documentation Center，1980.

[11]SHEIKH S I，PINES D J.Recursive estimation of spacecraft position using X-ray pulsar time of arrival measurements[C].ION 61st Annual Meeting，Cambridge，MA，2005：27-29.