修正鐘差的脈沖星方位誤差估計(jì)算法

肖永強(qiáng),王宏力,馮 磊,由四海,何貽洋,許 強(qiáng)

(1.火箭軍工程大學(xué)導(dǎo)彈工程學(xué)院,陜西 西安 710025;2.青州高新技術(shù)研究所測(cè)試控制系,山東 青州 262500)

0 引 言

X射線(xiàn)脈沖星導(dǎo)航是一種新興的自主導(dǎo)航方式,其利用脈沖星輻射的信號(hào)對(duì)航天器進(jìn)行定位、授時(shí)等服務(wù)。其自主性強(qiáng)、可靠性高,有望擺脫航天器對(duì)全球定位系統(tǒng)(global positioning system,GPS)等衛(wèi)星導(dǎo)航的依賴(lài),在民用和軍事領(lǐng)域有著巨大的發(fā)展前景[1-10]。但是目前要將其真正應(yīng)用到實(shí)際工程中仍存在很大的差距,其中一個(gè)很重要的原因就是脈沖星方位誤差的影響[11-16]。

在脈沖星導(dǎo)航過(guò)程中,0.001″的脈沖星方位誤差會(huì)造成幾百米的定位誤差[17]。因此,為了提高脈沖星方位誤差估計(jì)精度,國(guó)內(nèi)學(xué)者采用了基于信標(biāo)衛(wèi)星的估計(jì)[18]、魯棒濾波估計(jì)[19-21]、組合導(dǎo)航[22-24]等方法對(duì)脈沖星方位誤差進(jìn)行估計(jì)。此外,許強(qiáng)在利用信標(biāo)衛(wèi)星進(jìn)行估計(jì)時(shí)考慮了實(shí)際存在的衛(wèi)星位置誤差的影響[25]。但是上述研究都沒(méi)有考慮航天器的時(shí)鐘鐘差影響。雖然孫守明等研究了脈沖星導(dǎo)航以及脈沖星與慣性、多普勒等組合導(dǎo)航中的鐘差修正[26-28]問(wèn)題,但是在利用信標(biāo)衛(wèi)星進(jìn)行脈沖星方位誤差估計(jì)時(shí),依然會(huì)受到時(shí)鐘鐘差的影響。然而從目前公開(kāi)的文獻(xiàn)來(lái)看,利用信標(biāo)衛(wèi)星進(jìn)行脈沖星方位誤差估計(jì)中的鐘差修正問(wèn)題尚未有學(xué)者研究。當(dāng)利用信標(biāo)衛(wèi)星進(jìn)行脈沖星方位誤差估計(jì)時(shí),由于航天器長(zhǎng)時(shí)間運(yùn)行,時(shí)鐘會(huì)發(fā)生漂移,進(jìn)而會(huì)造成系統(tǒng)偏差,從而使脈沖星方位誤差估計(jì)精度降低。

為了解決脈沖星方位誤差估計(jì)時(shí)的鐘差影響,本文通過(guò)分析得出時(shí)鐘鐘差對(duì)脈沖星方位誤差估計(jì)的影響為緩變的過(guò)程,可看作系統(tǒng)偏差,常用的處理方式是利用增廣狀態(tài)法。但是若利用增廣狀態(tài)法處理,當(dāng)系統(tǒng)狀態(tài)與系統(tǒng)偏差維數(shù)相近時(shí),會(huì)使濾波計(jì)算量增加,且易出現(xiàn)數(shù)值不穩(wěn)定的問(wèn)題[29]。因此,為了解決上述問(wèn)題,本文提出了一種考慮鐘差修正的兩步卡爾曼濾波(two-stage Kalman filter,TSKF)算法。

1 脈沖星方位誤差估計(jì)原理

圖1 脈沖星方位誤差估計(jì)原理Fig.1 Principle of pulsar position error estimation

轉(zhuǎn)換過(guò)程[1]為

(1)

(2)

式中,tsat為脈沖到達(dá)航天器的真實(shí)時(shí)間;n為真實(shí)的脈沖星單位方向矢量。設(shè)脈沖星的赤經(jīng)為α,赤緯為β,則滿(mǎn)足:

(3)

設(shè)(Δα,Δβ)為脈沖星方位誤差,則真實(shí)的脈沖星方位與帶誤差的脈沖星方位滿(mǎn)足

(4)

將式(4)代入式(3),進(jìn)行泰勒展開(kāi)并忽略二階及以上小項(xiàng)可得

(5)

式中,帶誤差的脈沖星方向矢量為

(6)

記脈沖星單位方向矢量誤差Δn為

(7)

則可以將式(5)表示為

n-Δn=n′

(8)

(9)

取狀態(tài)變量為X=[ΔαΔβ]T,則傳統(tǒng)脈沖星方位誤差估計(jì)算法的狀態(tài)方程和觀(guān)測(cè)方程為

(10)

(11)

式中,狀態(tài)轉(zhuǎn)移矩陣為

觀(guān)測(cè)矩陣為

式中，Wk和ηk為系統(tǒng)噪聲和量測(cè)噪聲。

2 時(shí)鐘鐘差影響分析

由于航天器時(shí)鐘頻率和相位的漂移,脈沖到達(dá)航天器的真實(shí)時(shí)間和航天器時(shí)鐘測(cè)得的脈沖到達(dá)航天器的時(shí)間之間存在時(shí)鐘鐘差,設(shè)航天器的時(shí)鐘鐘差為δt,則滿(mǎn)足

(12)

時(shí)鐘鐘差可通過(guò)估計(jì)相對(duì)于標(biāo)準(zhǔn)時(shí)間的鐘差、鐘差漂移率和鐘差漂移率的變化率獲得,因此可將衛(wèi)星時(shí)鐘鐘差模型[30]表示為

(13)

式中,x1、x2和x3分別為鐘差、鐘差漂移率和鐘差漂移率變化率;τ為時(shí)間間隔。取狀態(tài)變量為x=[x1，x2，x3]T,則鐘差離散過(guò)程模型可表示為

(14)

式中,狀態(tài)轉(zhuǎn)移矩陣為

(15)

式中,q1、q2和q3為噪聲的功率譜密度。

仿真取時(shí)鐘誤差漂移率為3.637 979×10-12,時(shí)鐘漂移率的變化率為6.66×10-18/s;根據(jù)銣原子鐘模型[26],取時(shí)鐘的噪聲譜密度分別為q1=1.11×10-22s,q2=2.22×10-32/s和q3=6.66×10-45/s3。給定時(shí)鐘初始時(shí)刻的鐘差為0,取時(shí)間間隔為1 s,則可得到鐘差隨時(shí)間的變化如圖2所示。

圖2 鐘差隨時(shí)間的變化Fig.2 Change of clock error with time

由圖2可得,雖然鐘差漂移較小,但隨著航天器的長(zhǎng)時(shí)間運(yùn)行,鐘差接近5×10-6s,該值與光速相乘理論上會(huì)引起1 500 s左右的定位誤差,將會(huì)對(duì)量測(cè)方程式(11)造成較大的影響,進(jìn)而會(huì)對(duì)脈沖星方位誤差估計(jì)造成嚴(yán)重的影響。

使用傳統(tǒng)脈沖星方位誤差估計(jì)算法進(jìn)行仿真,分別驗(yàn)證其在有無(wú)鐘差影響下的估計(jì)效果。以脈沖星B0531+21作為觀(guān)測(cè)脈沖星,其參數(shù)如表1所示。

表1 脈沖星B0531+21參數(shù)Table 1 Parameters of pulsar B0531+21

其中,P為脈沖周期,W為脈沖寬度,Fx為脈沖輻射光子流量,pf為一個(gè)脈沖周期內(nèi)脈沖輻射流量與平均輻射流量之比。

脈沖星的觀(guān)測(cè)噪聲方差[2]為

(16)

式中,A為探測(cè)器有效面積,本仿真中設(shè)為1 m2;Bx=0.005 ph·cm-2·s-1,為宇宙背景噪聲;d為脈沖寬度W與脈沖周期P之比;tobs為觀(guān)測(cè)時(shí)間,設(shè)為1 000 s;則可計(jì)算[25]得到σ=(77.69 m)2。脈沖星方位誤差設(shè)為(2 mas,2 mas),初始狀態(tài)設(shè)為0。

使用同一顆衛(wèi)星,分別在有無(wú)鐘差影響的情況下進(jìn)行仿真,衛(wèi)星軌道參數(shù)如表2所示,仿真結(jié)果如圖3所示。

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

圖3 有無(wú)鐘差的脈沖星方位誤差估計(jì)結(jié)果Fig.3 Estimation results of pulsar position error without correction of clock error

對(duì)比分析圖3可得,當(dāng)無(wú)鐘差影響時(shí),傳統(tǒng)脈沖星方位誤差估計(jì)算法能較為精確地估計(jì)方位誤差。但當(dāng)存在鐘差時(shí),估計(jì)結(jié)果誤差較大,尤其是赤緯誤差,接近80 mas;而且赤經(jīng)和赤緯誤差估計(jì)結(jié)果都沒(méi)有收斂到一個(gè)定值。可見(jiàn),時(shí)鐘鐘差會(huì)對(duì)脈沖星方位誤差估計(jì)精度產(chǎn)生較大的影響,因此非常有必要對(duì)航天器時(shí)鐘鐘差進(jìn)行修正。

3 修正鐘差的脈沖星方位誤差估計(jì)算法

當(dāng)考慮鐘差時(shí),將式(1)和式(2)相減可得新的觀(guān)測(cè)模型為

(17)

由第1.2節(jié)分析可知,鐘差變化是一個(gè)緩慢過(guò)程,常用的處理方法是增廣狀態(tài)法,即將系統(tǒng)狀態(tài)和鐘差組合為新的狀態(tài)變量進(jìn)行濾波解算。但是這樣就會(huì)使系統(tǒng)狀態(tài)由2維變?yōu)?維,不僅會(huì)增加運(yùn)算負(fù)擔(dān),還容易造成數(shù)值不穩(wěn)定的問(wèn)題[29]。因此,為了解決上述問(wèn)題,本文采用TSKF算法。

TSKF算法是由Friedland提出用于處理系統(tǒng)常值偏差問(wèn)題[31],Ignagni 在這基礎(chǔ)上將其應(yīng)用到處理緩變偏差上[32]。該算法解耦狀態(tài)估計(jì)與系統(tǒng)偏差估計(jì),可有效減小濾波過(guò)程中的計(jì)算量,提高數(shù)值計(jì)算穩(wěn)定性。王奕迪等在前期研究中也將其應(yīng)用到了脈沖星導(dǎo)航中[33-34]。

根據(jù)TSKF算法原理取第一步濾波的狀態(tài)量為X=[Δα,Δβ]T,第二步濾波狀態(tài)量b為鐘差δt,可將方位誤差估計(jì)的TSKF算法方程寫(xiě)為

一步濾波方程:

Xk/k-1=AkXk-1

(18)

(19)

(20)

Xk=Xk/k-1+Kk(Zk-HkXk/k-1)

(21)

Pk=(I-KkHk)Pk/k-1

(22)

二步濾波方程:

Uk=Ak-1Vk-1+Bk-1

(23)

Sk=HkUk+Ck

(24)

(25)

Vk=Uk-KkSk

(26)

(27)

(28)

最終估計(jì)結(jié)果為

(29)

式中,Qk為系統(tǒng)噪聲Wk的方差;Pk為狀態(tài)量的協(xié)方差;Rk為觀(guān)測(cè)噪聲ηk的方差;Bk為δt在式(11)中的驅(qū)動(dòng)方程;Ck為δt在觀(guān)測(cè)方程(12)中的驅(qū)動(dòng)方程。由上述分析可得,Bk=0,Ck=c·Φ。

仿真條件同第2節(jié),估計(jì)結(jié)果如圖4所示。

圖4 TSKF算法估計(jì)結(jié)果Fig.4 Estimation results of TSKF algorithm

為了進(jìn)一步驗(yàn)證所提算法的有效性,選取不同的鐘差對(duì)比傳統(tǒng)算法與本文所提算法的估計(jì)結(jié)果,如表3所示。

表3 不同鐘差條件下TSKF算法估計(jì)偏差Table 3 TSKF estimation bias under different clock errors

由圖4和表3可得,當(dāng)不考慮鐘差的傳統(tǒng)算法估計(jì)偏差較大時(shí),本文所提出的考慮鐘差的TSKF算法估計(jì)精度較高,赤經(jīng)和赤緯誤差估計(jì)精度均能保持在0.1 mas以?xún)?nèi),且都能收斂到2 mas左右?？梢?jiàn)本文所提算法能有效隔離時(shí)鐘鐘差的影響,使估計(jì)精度保持在無(wú)鐘差影響下的水平。

為比較TSKF與增廣狀態(tài)濾波算法(augmented state Kalman filter algorithm,ASKF)的精度,仿真條件同第2節(jié),得到仿真結(jié)果如圖5所示。由圖5可得,在10～50天內(nèi),TSKF赤經(jīng)估計(jì)結(jié)果振幅略高于ASKF,二者都能很快收斂到2 mas左右,最終估計(jì)結(jié)果精度相當(dāng),均在0.01 mas以?xún)?nèi)。0～50天,TSKF赤緯估計(jì)結(jié)果偏差大于ASKF,但處于可控范圍內(nèi),且TSKF收斂速度略快于ASKF,最終估計(jì)結(jié)果TSKF精度略低于ASKF,但偏差在0.06 mas以?xún)?nèi),仍在可控范圍內(nèi)。其他條件不變,對(duì)比兩種算法在不同鐘差條件下的估計(jì)結(jié)果,如表4所示。表4也驗(yàn)證了TSKF算法估計(jì)結(jié)果接近ASKF,使方位誤差估計(jì)精度保持在較高的水平。為比較兩者的計(jì)算量,統(tǒng)計(jì)上述仿真實(shí)驗(yàn)中兩種算法Matlab程序中的浮點(diǎn)計(jì)算次數(shù),得到TSKF算法133 480 860次,ASKF算法316 957 544次,可得TSKF運(yùn)算量遠(yuǎn)小于ASKF,運(yùn)算量減小了57.89%,有效地減小了計(jì)算負(fù)擔(dān)。

圖5 兩種算法估計(jì)結(jié)果對(duì)比Fig.5 Estimated results comparison of two algorithms

表4 不同鐘差條件下兩種算法估計(jì)偏差Table 4 Estimation bias of two algorithms under different

由上述分析可知,本文提出的TSKF算法估計(jì)精度較高,雖略低于ASKF算法估計(jì)精度,但處于可控范圍內(nèi),而TSKF算法計(jì)算量顯著低于ASKF算法,具有更高的運(yùn)算效率。

為驗(yàn)證TSKF算法的魯棒性,分析當(dāng)仿真用的噪聲方差與數(shù)據(jù)真實(shí)方差不一致時(shí)的估計(jì)結(jié)果。仿真時(shí)將真實(shí)值的觀(guān)測(cè)噪聲增大3倍,其他條件不變,即仿真中觀(guān)測(cè)噪聲方差R=(0.077 69 km)2,而真實(shí)值的觀(guān)測(cè)噪聲為0.233 07 km,估計(jì)結(jié)果如圖6所示。

分析圖6可知,當(dāng)仿真用的噪聲方差與數(shù)據(jù)真實(shí)方差不一致時(shí),在0～50天赤經(jīng)估計(jì)誤差略大于圖4估計(jì)結(jié)果,但估計(jì)結(jié)果也能在第80天左右收斂到2 mas左右。在0～100天赤緯估計(jì)結(jié)果偏差大于圖4赤緯估計(jì)結(jié)果,但在第100天后也能收斂到2 mas左右?？梢?jiàn),TSKF算法在仿真用的噪聲方差與數(shù)據(jù)真實(shí)方差不一致時(shí),仍能使估計(jì)結(jié)果收斂到無(wú)鐘差影響的狀態(tài)下。

圖6 噪聲不一致時(shí)的估計(jì)結(jié)果Fig.6 Estimated results with inconsistent noise

由于時(shí)鐘鐘差主要由初始鐘差、鐘差漂移率和鐘差漂移率變化率確定,因此為進(jìn)一步驗(yàn)證本文算法的魯棒性,在原鐘差模型的基礎(chǔ)上添加周期項(xiàng):

(30)

式中,A0和B0為幅值;Ts為衛(wèi)星軌道周期。仿真時(shí),鐘差真實(shí)值采用上述帶周期項(xiàng)的模型,幅值A(chǔ)0和B0都取為1×10-6s,其他條件不變,得到估計(jì)結(jié)果如圖7所示。

圖7 增加周期項(xiàng)的估計(jì)結(jié)果Fig.7 Estimated results of increasing period item

其他條件不變,選取不同幅值A(chǔ)0和B0的大小進(jìn)行仿真實(shí)驗(yàn),每次仿真實(shí)驗(yàn)進(jìn)行20次,并選取最后6天的估計(jì)結(jié)果取平均值,估計(jì)結(jié)果如表5所示。

表5 不同幅值條件下的估計(jì)偏差Table 5 Estimation bias under different amplitude conditions

分析上述圖表可得,當(dāng)鐘差中加入周期變化值而與模型值不一致時(shí),雖然仿真前期偏差較大,但隨著時(shí)間的推移,估計(jì)結(jié)果仍能收斂到2 mas左右。雖然與模型一致時(shí)會(huì)在2 mas附近上下抖動(dòng),但這種抖動(dòng)幅度是微小的,可以通過(guò)取平均值的方法來(lái)提高精度,也驗(yàn)證了TSKF算法具有較強(qiáng)的魯棒性。

4 結(jié) 論

(1) 利用航天器(信標(biāo)衛(wèi)星)進(jìn)行脈沖星方位誤差估計(jì)時(shí),航天器的時(shí)鐘鐘差會(huì)造成系統(tǒng)偏差,對(duì)估計(jì)結(jié)果產(chǎn)生影響,且這種偏差和影響是不容忽略的。

(2) 考慮鐘差修正的TSKF算法能有效隔離時(shí)鐘鐘差對(duì)脈沖星方位誤差估計(jì)的影響,在保證濾波解算穩(wěn)定性的同時(shí),能使估計(jì)精度保持在無(wú)鐘差影響下的水平。

(3) 與ASKF算法相比,TSKF算法具有相近的估計(jì)精度和更高的運(yùn)算效率,且TSKF算法具有較強(qiáng)的魯棒性。