基于雙濾波器的高精度目標(biāo)跟蹤估計(jì)方法

秦健凱，李鵬

北京交通大學(xué) 電子信息工程學(xué)院，北京 100044

1 引言

航天器是當(dāng)代獲取空間信息的主要方式，如何獲取太空中航天器的運(yùn)動狀態(tài)成為當(dāng)前重點(diǎn)研究的問題之一。

空間目標(biāo)在推進(jìn)器作用下改變原軌跡進(jìn)入新軌道，即軌道機(jī)動。空間目標(biāo)軌道機(jī)動分為連續(xù)推力和脈沖機(jī)動兩類，前一類可用于構(gòu)造機(jī)動模型并同時估計(jì)模型參數(shù)與目標(biāo)狀態(tài)，后一類可通過機(jī)動檢測策略來確定目標(biāo)機(jī)動時刻，從而做出應(yīng)對。對于存在脈沖機(jī)動的非合作目標(biāo)，目前多采用天基平臺進(jìn)行觀測，并設(shè)計(jì)相應(yīng)的機(jī)動檢測策略。天基平臺的觀測值具有非線性的特點(diǎn)，因此非線性濾波成為空間目標(biāo)跟蹤估計(jì)的主要方法之一[1]。

現(xiàn)有研究基于擴(kuò)展卡爾曼濾波、雙擴(kuò)展卡爾曼濾波及改進(jìn)的雙重?zé)o跡卡爾曼濾波(DUKF)進(jìn)行軌道或目標(biāo)參數(shù)估計(jì)，獲得了效果提升。但注意到簡化攝動力的模型會導(dǎo)致累積誤差，使用雙擴(kuò)展卡爾曼濾波對復(fù)雜非線性模型失效，或?qū)顟B(tài)和參數(shù)分別采用不同的觀測方程估計(jì)，會導(dǎo)致解算復(fù)雜[2-4]。為此，本文考慮只用一個觀測方程來精簡描述，即選取加速度作為觀測信息，方程既包含狀態(tài)量又有參數(shù)量。在考慮太陽光壓的基礎(chǔ)上，基于DUKF進(jìn)行估計(jì)，從而克服建模不準(zhǔn)確及解算復(fù)雜的缺點(diǎn)，實(shí)現(xiàn)合作目標(biāo)運(yùn)動狀態(tài)以及機(jī)動加速度的高效估計(jì)。

對非合作航天器，無法獲得其加速度觀測值以及機(jī)動信息。現(xiàn)有研究采用相對運(yùn)動方程與濾波器切換策略實(shí)現(xiàn)目標(biāo)跟蹤，并對多種濾波算法的精度和收斂速度進(jìn)行對比研究[5-7]。文獻(xiàn)[8-9]以半長軸變化為變軌依據(jù)，通過聯(lián)合策略實(shí)現(xiàn)機(jī)動下的半長軸檢測，簡化處理偏心率變化的情況。此外，對于地基觀測，地面測控存在時延，且對于一些無法在全球布置測控站的地區(qū)還存在信號遮擋問題[10]。文獻(xiàn)[11-12]綜合考慮兩個軌道根數(shù)，將偏心率作為輔助特征參數(shù)，設(shè)計(jì)目標(biāo)狀態(tài)噪聲自適應(yīng)調(diào)整策略。但要注意部分觀測方法的限制較強(qiáng)，且復(fù)雜的模型方程會對濾波效率帶來負(fù)面影響。為此，本文采用天基觀測平臺獲得非合作目標(biāo)的觀測值，用半正焦弦代替半長軸，當(dāng)檢測到目標(biāo)機(jī)動后進(jìn)行濾波器切換，從而克服求解半長軸時可能無窮大的缺陷及對地基觀測的依賴。進(jìn)而，采用雙EKF實(shí)現(xiàn)基于半正焦弦機(jī)動檢測策略的非合作目標(biāo)的高精度位置速度跟蹤估計(jì)。

論文總體方案流程如圖1所示。采用DUKF算法對存在連續(xù)推力的合作目標(biāo)進(jìn)行跟蹤估計(jì)；采用雙EKF結(jié)合機(jī)動檢測策略對存在未知機(jī)動的非合作目標(biāo)進(jìn)行跟蹤估計(jì)。為降低虛警，設(shè)定當(dāng)連續(xù)檢測到3個點(diǎn)滿足要求時，認(rèn)為目標(biāo)發(fā)生機(jī)動。

圖1 總體方案流程Fig.1 The flow chart of the overall scheme

2 理論基礎(chǔ)

合作航天器可通過加速度計(jì)獲取信息，而非合作航天器在信息獲取方面不配合，可以將一個合作航天器作為天基觀測平臺來獲取非合作目標(biāo)的觀測信息，建立相對運(yùn)動方程，描述其運(yùn)動狀態(tài)。

2.1 航天器運(yùn)動方程

(1)合作航天器運(yùn)動方程

考慮可以產(chǎn)生連續(xù)推力的合作航天器，受到中心天體引力以及太陽光壓攝動的影響，假設(shè)推力加速度方向與參考坐標(biāo)系下三軸方向一致。目標(biāo)的狀態(tài)量X=[rvε]T，其中r為目標(biāo)在參考系下的位置矢量，v為速度矢量，ε為反射因子。假設(shè)太陽光照有效面積不變，且考慮太陽光不能時刻垂直照射太陽翼，故用ε作為變量體現(xiàn)這一變化。目標(biāo)的參數(shù)量θ=[m]其中m為航天器的質(zhì)量。太陽光壓攝動加速度為:

式中：AU為天文單位，1AU=1.496×108km；S為航天器受輻射有效面積；ρ=4.56×10-6N/m2為太陽輻射壓強(qiáng)；CR=1+ε為太陽光反射系數(shù)；ε∈[0,1]。

機(jī)動加速度為：

式中：F=[FxFyFz]T為推力矢量；m0為航天器初始質(zhì)量；Δt為機(jī)動持續(xù)時間。Isp為發(fā)動機(jī)比沖，gT為參考點(diǎn)引力加速度。合作航天器通過其攜帶的加速度計(jì)，可以獲得除引力外的加速度大小，且初始質(zhì)量、推進(jìn)器類型等信息已知。

(2)非合作航天器運(yùn)動方程

天基觀測平臺作為合作目標(biāo)，可通過獲得其運(yùn)動狀態(tài)建立天基平臺下的軌道坐標(biāo)系，通過相對運(yùn)動方程描述非合作航天器的運(yùn)動狀態(tài)。天基平臺通?；谖⒉ɡ走_(dá)、光學(xué)相機(jī)等獲得目標(biāo)的相對距離及角度信息。對機(jī)動未知的非合作航天器，通過半正焦弦的變化來檢測機(jī)動，之后切換濾波器，可在短時間內(nèi)收斂并具有較高的濾波精度。

式中：ρ=[xyz]T為軌道坐標(biāo)系中非合作航天器的相對位置；n為角速度；μ為中心天體引力常數(shù)；d為天基平臺軌道半徑；相對加速度f=fT-fC=[fxfyfz]T，則觀測方程為：

式中：Yk為觀測矢量；A為偏航角；H為高度角；R為相對距離；εk為已知高斯白噪聲；x=[xkykzk]T為航天器軌道坐標(biāo)系下的位置矢量。

2.2 濾波算法

選取適當(dāng)?shù)臑V波算法，利用自主導(dǎo)航處理測量信息高精度估計(jì)軌道信息的核心環(huán)節(jié)。本文基于擴(kuò)展卡爾曼濾波(EKF)和無跡卡爾曼濾波(UKF)[13]，通過融合雙濾波及機(jī)動檢測策略實(shí)現(xiàn)高效跟蹤估計(jì)。

考慮非線性狀態(tài)方程如下[14]：

xk+1=f1(xk,θk)+Vk

zk=h(xk,θk)+nk

式中：nk～N(0,Rk)為觀測噪聲。

UKF算法簡述如下。

(1)預(yù)測過程

構(gòu)造總數(shù)為2L+1的sigma點(diǎn)集：

(i=1,2,…,L)

(i=L+1,L+2,…,2L)

通過狀態(tài)轉(zhuǎn)移函數(shù)f(·)，將sigma點(diǎn)集映射到新sigma點(diǎn)集：

權(quán)重公式為：

λ=α2(L+K)-L

其中參數(shù)取值：α=10-3，K=0，β=2，加權(quán)后的sigma點(diǎn)集用于預(yù)測狀態(tài)的估計(jì)值和協(xié)方差：

(2)更新過程

通過觀測函數(shù)h(·)，將sigma點(diǎn)集映射到新sigma點(diǎn)集：

加權(quán)新的sigma點(diǎn)集后用于預(yù)測觀測的估計(jì)值和協(xié)方差：

狀態(tài)測量協(xié)方差矩陣用于計(jì)算卡爾曼增益：

狀態(tài)更新如下：

圖2 DUKF算法流程Fig.2 The flow chart of dual unscented Kalman filter

UKF1為狀態(tài)濾波器，UKF2為參數(shù)濾波器；x為狀態(tài)變量，θ為參數(shù)變量。DUKF中平行運(yùn)行2個UKF濾波器，使用相同觀測方程，分別用于狀態(tài)估計(jì)和參數(shù)估計(jì)，且在估計(jì)過程中實(shí)時校正參數(shù)，并對運(yùn)動狀態(tài)的估計(jì)進(jìn)行反饋，提升精度。

2.3 基于DUKF的合作目標(biāo)運(yùn)動狀態(tài)估計(jì)

對于雙重?zé)o跡卡爾曼濾波器，其狀態(tài)濾波器為：

即εk+1=εk+q3。其中，q1、q2、q3均為已知的高斯白噪聲。

參數(shù)濾波器如下：

式中：q4為已知的高斯白噪聲。

觀測方程如下：

a=af+asun+q5

式中：q5為已知的高斯白噪聲。

對上述狀態(tài)濾波器，質(zhì)量m為已知量，由參數(shù)濾波器給出；對參數(shù)濾波器，反射因子?為已知量，由狀態(tài)濾波器給出。兩個濾波器之間相互配合來提高濾波精度。

2.4 基于半正焦弦機(jī)動檢測策略的非合作目標(biāo)運(yùn)動狀態(tài)估計(jì)

對于非合作目標(biāo)，在濾波跟蹤估計(jì)的過程中可能會存在主動且未知機(jī)動，導(dǎo)致濾波精度下降甚至?xí)l(fā)散。因此需要通過機(jī)動檢測的手段，來應(yīng)對隨時可能發(fā)生的未知機(jī)動，故有必要進(jìn)行非合作目標(biāo)的機(jī)動檢測策略的設(shè)計(jì)[15-16]。

針對軌道動力學(xué)非線性的特性，可采用EKF來解決隨機(jī)連續(xù)非線性系統(tǒng)問題[17]。本節(jié)在EKF的基礎(chǔ)上使用雙EKF實(shí)現(xiàn)存在一次或多次機(jī)動的非合作目標(biāo)跟蹤。脈沖推力產(chǎn)生的機(jī)動加速度可表示為[18]：

apulse=Δvδ(t-tc)

式中：Δv為脈沖推力產(chǎn)生的速度變換量；t為機(jī)動時刻；δ(t-tc)為tc的沖激函數(shù)。

對脈沖推力動力學(xué)方程兩邊求tc-τ到tc+τ的定積分，當(dāng)τ趨于0時，可得：

rf-ro=0

vf-vo=Δv

式中：下標(biāo)o表示機(jī)動前，f表示機(jī)動后。

由此可見，機(jī)動前后瞬間目標(biāo)位置矢量基本不變，而速度矢量有較大變化Δv，且與推力同向。

當(dāng)非合作目標(biāo)航天器主動機(jī)動后，其軌道根數(shù)中變化最為明顯的是半長軸a，但軌道根數(shù)求解半長軸可能無窮大，因此本文使用半正焦弦p作為特征參數(shù)。

圖3為橢圓示意圖。

圖3 橢圓示意Fig.3 Schematic diagram of ellipse

圖3中O為極坐標(biāo)原點(diǎn)，F(xiàn)為右焦點(diǎn)，|OF| =c，橢圓半長軸為a，短半軸為b，M為橢圓上一點(diǎn)(r,φ)，|MF|=r，|FK|=Q。中心天體位于點(diǎn)F處，偏心率e=c/a，由x2/a2+y2/b2=1，當(dāng)x=c時，y=±b2/a，取正值即為半正焦弦：

p=b2/a=(a2-c2)/a=a(1-e2)

可見，當(dāng)a和e其中一個變化或二者均變化時，p均會改變。e∈[0,1]，a對p的影響更為明顯，對于后續(xù)公式(1)中閾值γ選取時要綜合考慮以上因素。

當(dāng)檢測到機(jī)動時，由于受力改變，仍用原來濾波器進(jìn)行濾波可能會發(fā)散，因此要及時切換到另外一個不同參數(shù)的濾波器，再進(jìn)行后續(xù)濾波?？臻g中的航天器會受到多種因素的干擾，即使無機(jī)動時特征參數(shù)也會在一定范圍內(nèi)發(fā)生小幅度變化，因此設(shè)計(jì)如下的檢測策略。

假設(shè)前n1秒無機(jī)動，根據(jù)這段時間內(nèi)的半正焦弦估計(jì)值統(tǒng)計(jì)其均值及方差，代入式(1)。

(1)

當(dāng)滿足式(1)時，表示檢測到機(jī)動，需要更換濾波器參數(shù)繼續(xù)濾波；不滿足時，繼續(xù)使用原濾波器。為了減少誤報(bào)，可以設(shè)置檢測到連續(xù)若干點(diǎn)滿足上式時表示機(jī)動開始。本文設(shè)置連續(xù)檢測到3個點(diǎn)滿足上式時，視為機(jī)動開始。

此處的雙EKF不同于前文的DUKF，前者指的是兩個獨(dú)立的EKF，由于機(jī)動前后運(yùn)動狀態(tài)發(fā)生改變，所需的最優(yōu)濾波參數(shù)一定不同，因此使用兩個不同參數(shù)的EKF，檢測到機(jī)動后通過切換濾波器實(shí)現(xiàn)短時間收斂。濾波參數(shù)指的是初始誤差協(xié)方差、狀態(tài)和觀測噪聲協(xié)方差等，不同于DUKF中的參數(shù)變量。DUKF的狀態(tài)方程與狀態(tài)量和參數(shù)量均相關(guān)，參數(shù)對于狀態(tài)的估計(jì)有很大的影響，若參數(shù)變化較大，卻把它考慮為一個定值就會使估計(jì)結(jié)果有很大的偏差，因此使用一個濾波器對參數(shù)量進(jìn)行同步估計(jì)、實(shí)時校正，通過兩個濾波器間的相互嵌套和信息利用對狀態(tài)估計(jì)進(jìn)行反饋校正，從而提高濾波精度。

3 仿真分析

3.1 合作目標(biāo)跟蹤估計(jì)仿真結(jié)果

考慮具有連續(xù)常推力軌道機(jī)動，通過UKF和DUKF實(shí)現(xiàn)對于合作目標(biāo)的跟蹤估計(jì)。目標(biāo)航天器初始軌道根數(shù)為：半長軸a0=706 300 km，偏心率e0=0.146，軌道傾角i0=30.394/(180π)°，近地點(diǎn)幅角w0=74.373/(180π)°，升交點(diǎn)赤經(jīng)Ω0=42.790/(180π)°，真近點(diǎn)角η0=45.937/(180π)°，初始反射因子ε0=0.5，航天器初始總質(zhì)量m0=1 000 kg，常推力10 N。濾波周期設(shè)定為T=1 s，仿真時間2 000 s。初始誤差協(xié)方差矩陣為：

P0=diag[1 1 1 0.01 0.01 0.01 1×10-4]

狀態(tài)過程噪聲協(xié)方差矩陣為：

Q1=diag[9×10-18, 9×10-18,9×10-18,

1×10-12,1×10-12,1×10-12, 1×10-4]

參數(shù)過程噪聲協(xié)方差矩陣為：

Q2=diag[1×10-5]

觀測噪聲協(xié)方差矩陣為：

R=diag[3×10-13, 3×10-13, 3×10-13]

模型中存在不確定參數(shù)時可考慮使用擴(kuò)維UKF，即在UKF的基礎(chǔ)上把參數(shù)也作為狀態(tài)變量一起進(jìn)行估計(jì)。擴(kuò)維UKF和DUKF濾波仿真結(jié)果如圖4和圖5所示。

圖4 UKF與DUKF的三軸位置估計(jì)誤差Fig.4 Three-axes position estimation errors using UKF and DUKF

圖5 UKF與DUKF的三軸速度估計(jì)誤差Fig.5 Three-axes velocity estimation errors using UKF and DUKF

圖4和圖5表明，擴(kuò)維UKF和DUKF估計(jì)的位置和速度均收斂，為了比較兩種方法的精度與收斂速度，將x方向的位置誤差和y方向的速度誤差仿真圖放大，如圖6所示。

圖6 UKF與DUKF在x軸放大的位置估計(jì)誤差和y軸放大的速度估計(jì)誤差Fig.6 Magnified position estimation errors along axis x and velocity estimation errors along axis y of UKF and DUKF

擴(kuò)維UKF最大位置誤差約為2.4 km，DUKF約為2.1 km；擴(kuò)維UKF最大速度誤差約為0.02 km/s，收斂時間約為348 s，DUKF約為0.02 km/s和218 s。DUKF相對于擴(kuò)維UKF收斂更快，因?yàn)閰?shù)濾波器分擔(dān)了一部分狀態(tài)濾波器的工作，降低了其維數(shù)，且協(xié)同工作使得最大誤差值變小，精度提高。從圖7可以看出，DUKF對于質(zhì)量(加速度)的估計(jì)收斂更快、精度更高。

圖7 質(zhì)量估計(jì)誤差對比Fig.7 Quality estimation errors of UKF and DUKF

綜上可知：將質(zhì)量等因素?cái)U(kuò)展到原狀態(tài)量中的擴(kuò)維UKF濾波方法在理論上更直觀，但實(shí)際上增加了維數(shù)和計(jì)算量；DUKF雙重濾波算法將狀態(tài)信息與參數(shù)信息進(jìn)行融合[19]，利用聯(lián)合密度函數(shù)邊緣化處理將組合估計(jì)問題分解成狀態(tài)估計(jì)子問題和參數(shù)估計(jì)子問題；且兩個濾波器協(xié)同運(yùn)行可有效提高濾波精度。

3.2 非合作目標(biāo)跟蹤估計(jì)仿真結(jié)果

本節(jié)對非合作目標(biāo)相對導(dǎo)航算法進(jìn)行仿真。設(shè)前200 s無機(jī)動，目標(biāo)航天器分別在501 s、1 001 s、2 001 s處設(shè)置機(jī)動，大小分別為10-6km/s2、3×10-6km/s2、5×10-6km/s2。天基平臺初始軌道根數(shù)為：

ac=7 122.30 km，ec=0,ic=(97.750 π/180)°,wc=(311.559 π/180)°,Ωc=(354.829 π/180)°,θc=(160.18 π/180)°。

在天基平臺軌道坐標(biāo)系下，相對位置和速度初值[x,y,z,vx,vy,vz]為：

[0.025,0.025,0.050,-0.388 9×10-3,

-0.439 2×10-3,0.824 6×10-3]

濾波周期設(shè)定為T=1 s，仿真時間3 000 s，初始誤差協(xié)方差矩陣為:

Pk0=diag[1×10-8， 1×10-8， 1×10-8，

1×10-12， 1×10-12， 1×10-12]

過程噪聲協(xié)方差矩陣為：

Qk1=diag[1×10-14， 1×10-14， 1×10-14，

1×10-16，1×10-16， 1×10-16]

觀測噪聲協(xié)方差矩陣為：

Rk1=diag[1×10-12， 1×10-12， 1×10-8]

切換后的濾波器參數(shù)為：

Qk2=diag[1×10-12， 1×10-12， 1×10-12，

1×10-14，1×10-14， 1×10-14]

Rk2=diag[1×10-12， 1×10-12， 1×10-16]

仿真結(jié)果如圖8、圖9所示。將z方向相對位置和速度誤差放大，如圖10所示。

圖8 有/無機(jī)動檢測的相對位置估計(jì)誤差Fig.8 Relative position estimation errors with and without detection strategy

圖9 有/無機(jī)動檢測的相對速度估計(jì)誤差Fig.9 Relative velocity estimation errors with and without detection strategy

圖10 在z軸放大的相對位置和速度估計(jì)誤差Fig.10 Magnified estimation errors along axis z for relative position and velocity

從圖8～圖10可以看出，若無機(jī)動檢測機(jī)制，機(jī)動發(fā)生后濾波器再次收斂速度較慢，約300 s。這是因?yàn)闄C(jī)動使航天器的運(yùn)動狀態(tài)發(fā)生較大改變，原參數(shù)無法繼續(xù)使濾波器達(dá)到最優(yōu)效果，雖然經(jīng)過較長時間后也能收斂，但在實(shí)際工程中可能造成目標(biāo)丟失。仿真最大相對位置誤差約為6.5 m，最大相對速度誤差約為0.12 m/s。

在設(shè)置機(jī)動檢測的情況下，分別在500 s、1 000 s、2 000 s左右共檢測到3次機(jī)動，只需很短時間即可實(shí)現(xiàn)機(jī)動檢測并切換到另一個濾波器后再次收斂，最大相對位置誤差約為0.1 m，最大相對速度誤差約為0.04 m/s，最大位置和速度誤差相比于無機(jī)動檢測的結(jié)果均減小。在3次機(jī)動情況下，設(shè)置機(jī)動檢測可保證快速收斂，且濾波精度相比于無檢測明顯提高。這主要因?yàn)闄C(jī)動值變大后，運(yùn)動狀態(tài)發(fā)生更大改變，無檢測策略的濾波算法收斂慢，因此產(chǎn)生較大誤差。相比之下，設(shè)置機(jī)動檢測的策略可有效抑制這一點(diǎn)，可以在檢測到連續(xù)3個超出半正焦弦值統(tǒng)計(jì)范圍時切換濾波器，由此降低了誤報(bào)的可能性，保證了濾波跟蹤的精度和速度。

4 結(jié)論

本文研究了合作與非合作航天器的高精度目標(biāo)跟蹤估計(jì)方法。針對帶加速度計(jì)測量的合作航天器，考慮太陽光壓攝動，采用DUKF對其運(yùn)動狀態(tài)及質(zhì)量進(jìn)行雙重估計(jì)，實(shí)現(xiàn)對目標(biāo)的跟蹤估計(jì)；與擴(kuò)維UKF相比，DUKF獲得了更好的濾波精度和收斂速度。針對機(jī)動信息未知的非合作航天器，在天基平臺軌道系下建立相對運(yùn)動方程，采用半正焦弦機(jī)動檢測策略檢測機(jī)動并切換濾波器，結(jié)合雙EKF實(shí)現(xiàn)目標(biāo)跟蹤估計(jì)，其與無檢測策略相比，收斂速度更快，精度更高，且可實(shí)現(xiàn)多次機(jī)動的檢測。

本文對非合作航天器的估計(jì)是在使用DUKF獲得天基平臺(合作目標(biāo))運(yùn)動狀態(tài)的基礎(chǔ)上進(jìn)行的，所采用的基于半正焦弦的檢測策略還難以實(shí)現(xiàn)對于半長軸和偏心率同時具有較大變化的檢測，值得后續(xù)深入系統(tǒng)的研究。