基于雙濾波器的高精度目標跟蹤估計方法

秦健凱，李鵬

北京交通大學 電子信息工程學院，北京 100044

1 引言

航天器是當代獲取空間信息的主要方式，如何獲取太空中航天器的運動狀態成為當前重點研究的問題之一。

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

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

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

論文總體方案流程如圖1所示。采用DUKF算法對存在連續推力的合作目標進行跟蹤估計；采用雙EKF結合機動檢測策略對存在未知機動的非合作目標進行跟蹤估計。為降低虛警，設定當連續檢測到3個點滿足要求時，認為目標發生機動。

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

2 理論基礎

合作航天器可通過加速度計獲取信息，而非合作航天器在信息獲取方面不配合，可以將一個合作航天器作為天基觀測平臺來獲取非合作目標的觀測信息，建立相對運動方程，描述其運動狀態。

2.1 航天器運動方程

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

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

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

機動加速度為：

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

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

天基觀測平臺作為合作目標，可通過獲得其運動狀態建立天基平臺下的軌道坐標系，通過相對運動方程描述非合作航天器的運動狀態。天基平臺通常基于微波雷達、光學相機等獲得目標的相對距離及角度信息。對機動未知的非合作航天器，通過半正焦弦的變化來檢測機動，之后切換濾波器，可在短時間內收斂并具有較高的濾波精度。

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

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

2.2 濾波算法

選取適當的濾波算法，利用自主導航處理測量信息高精度估計軌道信息的核心環節。本文基于擴展卡爾曼濾波(EKF)和無跡卡爾曼濾波(UKF)[13]，通過融合雙濾波及機動檢測策略實現高效跟蹤估計。

考慮非線性狀態方程如下[14]：

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

zk=h(xk,θk)+nk

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

UKF算法簡述如下。

(1)預測過程

構造總數為2L+1的sigma點集：

(i=1,2,…,L)

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

通過狀態轉移函數f(·)，將sigma點集映射到新sigma點集：

權重公式為：

λ=α2(L+K)-L

其中參數取值：α=10-3，K=0，β=2，加權后的sigma點集用于預測狀態的估計值和協方差：

(2)更新過程

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

加權新的sigma點集后用于預測觀測的估計值和協方差：

狀態測量協方差矩陣用于計算卡爾曼增益：

狀態更新如下：

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

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

2.3 基于DUKF的合作目標運動狀態估計

對于雙重無跡卡爾曼濾波器，其狀態濾波器為：

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

參數濾波器如下：

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

觀測方程如下：

a=af+asun+q5

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

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

2.4 基于半正焦弦機動檢測策略的非合作目標運動狀態估計

對于非合作目標，在濾波跟蹤估計的過程中可能會存在主動且未知機動，導致濾波精度下降甚至會發散。因此需要通過機動檢測的手段，來應對隨時可能發生的未知機動，故有必要進行非合作目標的機動檢測策略的設計[15-16]。

針對軌道動力學非線性的特性，可采用EKF來解決隨機連續非線性系統問題[17]。本節在EKF的基礎上使用雙EKF實現存在一次或多次機動的非合作目標跟蹤。脈沖推力產生的機動加速度可表示為[18]：

apulse=Δvδ(t-tc)

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

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

rf-ro=0

vf-vo=Δv

式中：下標o表示機動前，f表示機動后。

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

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

圖3為橢圓示意圖。

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

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

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

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

當檢測到機動時，由于受力改變，仍用原來濾波器進行濾波可能會發散，因此要及時切換到另外一個不同參數的濾波器，再進行后續濾波。空間中的航天器會受到多種因素的干擾，即使無機動時特征參數也會在一定范圍內發生小幅度變化，因此設計如下的檢測策略。

假設前n1秒無機動，根據這段時間內的半正焦弦估計值統計其均值及方差，代入式(1)。

(1)

當滿足式(1)時，表示檢測到機動，需要更換濾波器參數繼續濾波；不滿足時，繼續使用原濾波器。為了減少誤報，可以設置檢測到連續若干點滿足上式時表示機動開始。本文設置連續檢測到3個點滿足上式時，視為機動開始。

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

3 仿真分析

3.1 合作目標跟蹤估計仿真結果

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

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

狀態過程噪聲協方差矩陣為：

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

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

參數過程噪聲協方差矩陣為：

Q2=diag[1×10-5]

觀測噪聲協方差矩陣為：

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

模型中存在不確定參數時可考慮使用擴維UKF，即在UKF的基礎上把參數也作為狀態變量一起進行估計。擴維UKF和DUKF濾波仿真結果如圖4和圖5所示。

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

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

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

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

擴維UKF最大位置誤差約為2.4 km，DUKF約為2.1 km；擴維UKF最大速度誤差約為0.02 km/s，收斂時間約為348 s，DUKF約為0.02 km/s和218 s。DUKF相對于擴維UKF收斂更快，因為參數濾波器分擔了一部分狀態濾波器的工作，降低了其維數，且協同工作使得最大誤差值變小，精度提高。從圖7可以看出，DUKF對于質量(加速度)的估計收斂更快、精度更高。

圖7 質量估計誤差對比Fig.7 Quality estimation errors of UKF and DUKF

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

3.2 非合作目標跟蹤估計仿真結果

本節對非合作目標相對導航算法進行仿真。設前200 s無機動，目標航天器分別在501 s、1 001 s、2 001 s處設置機動，大小分別為10-6km/s2、3×10-6km/s2、5×10-6km/s2。天基平臺初始軌道根數為：

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

在天基平臺軌道坐標系下，相對位置和速度初值[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]

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

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

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

過程噪聲協方差矩陣為：

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

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

觀測噪聲協方差矩陣為：

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

切換后的濾波器參數為：

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]

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

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

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

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

從圖8～圖10可以看出，若無機動檢測機制，機動發生后濾波器再次收斂速度較慢，約300 s。這是因為機動使航天器的運動狀態發生較大改變，原參數無法繼續使濾波器達到最優效果，雖然經過較長時間后也能收斂，但在實際工程中可能造成目標丟失。仿真最大相對位置誤差約為6.5 m，最大相對速度誤差約為0.12 m/s。

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

4 結論

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

本文對非合作航天器的估計是在使用DUKF獲得天基平臺(合作目標)運動狀態的基礎上進行的，所采用的基于半正焦弦的檢測策略還難以實現對于半長軸和偏心率同時具有較大變化的檢測，值得后續深入系統的研究。