抗差自適應(yīng)UKF算法在地基光學(xué)跟蹤空間目標(biāo)中的應(yīng)用

劉光明, 徐帆江

(中國科學(xué)院軟件研究所天基綜合信息系統(tǒng)重點(diǎn)實(shí)驗(yàn)室, 北京 100190)

0 引 言

空間目標(biāo)跟蹤定軌技術(shù)是空間目標(biāo)監(jiān)視系統(tǒng)應(yīng)用的基礎(chǔ),在當(dāng)今航天精確定軌領(lǐng)域,改進(jìn)測控設(shè)備的測控精度并且提高空間目標(biāo)的軌道確定精度是研究熱點(diǎn)[1-2]。抗差自適應(yīng)濾波算法的處理流程為[3-4]:如果觀測值存在異常,就根據(jù)抗差估計(jì)原則對值進(jìn)行處理;如果動(dòng)力學(xué)模型發(fā)生誤差,則對模型的各個(gè)分量使用自適應(yīng)的方法進(jìn)行調(diào)節(jié),或者將模型當(dāng)作整體,使用統(tǒng)一的自適應(yīng)因子的方法改變模型對狀態(tài)參數(shù)的影響[5]。非線性濾波是空間目標(biāo)實(shí)時(shí)跟蹤的有效方法之一[6],在實(shí)時(shí)跟蹤任務(wù)中,測量異常值(野值)的存在可能導(dǎo)致目標(biāo)跟蹤誤差突增甚至丟失,文獻(xiàn)[6]討論了實(shí)時(shí)定軌的抗差估計(jì)方法,針對無跡卡爾曼濾波(unscented Kalman filter,UKF)算法在未知時(shí)變噪聲下出現(xiàn)的濾波精度下降或者發(fā)散的問題。文獻(xiàn)[7]在其基礎(chǔ)上提出了SPRHKF算法,此算法對模型的瞬時(shí)干擾誤差和不確定性具有一定的魯棒性。上述算法雖然在不同程度上改進(jìn)了UKF算法的自適應(yīng)能力,但是仍然存在一些不足[8-9]:算法普遍存在計(jì)算復(fù)雜并且引入其他限制條件而導(dǎo)致難以工程實(shí)現(xiàn)的問題,另外雖然解決了UKF算法的發(fā)散問題,但是其收斂性較差,例如對初始狀態(tài)誤差要求盡可能小[10]。對于傳統(tǒng)UKF算法在跟蹤濾波的過程中由于噪聲統(tǒng)計(jì)特性時(shí)變的影響而出現(xiàn)濾波精度下降乃至發(fā)散的問題,本文提出了新的抗差自適應(yīng)UKF算法,該算法基于極大后驗(yàn)估計(jì),仿真實(shí)例結(jié)果驗(yàn)證了算法在不同條件下的有效性和可行性。

1 目標(biāo)觀測模型

地面站對空間非合作目標(biāo)的觀測幾何如圖1所示。由于僅有一個(gè)跟蹤目標(biāo)的觀測平臺,每個(gè)時(shí)刻對應(yīng)的觀測方程都是非正定的。所以需要使用軌道動(dòng)力學(xué)模型求解目標(biāo)軌道,使用在觀測時(shí)間序列中測量的數(shù)據(jù)進(jìn)行計(jì)算,實(shí)際處理時(shí)多使用基于軌道動(dòng)力學(xué)的濾波算法[11]。

圖1 地面站對空間非合作目標(biāo)的觀測幾何Fig.1 Observation geometry of ground station to space non cooperative target

地基光學(xué)測角跟蹤中,觀測平臺獲得空間目標(biāo)的測角信息(赤經(jīng)β、赤緯ε),赤經(jīng)和赤緯可以由J2000地心慣性坐標(biāo)系下目標(biāo)位置矢量r=[x,y,z]T和地面站位置矢量R=[X,Y,Z]T確定[12]。

(1)

則觀測矢量Y=[β,ε]T可以表示為狀態(tài)矢量X=[rT,RT]T的非線性函數(shù)：

Y=H(X)+v

(2)

式中,v為測量噪聲矢量。

2 目標(biāo)跟蹤系統(tǒng)模型

一段觀測時(shí)間序列下,構(gòu)造目標(biāo)狀態(tài)轉(zhuǎn)移函數(shù)xk=f(xk-1),目標(biāo)的角度測量序列為yk,則可建立地面站對空間非合作目標(biāo)的光學(xué)測角跟蹤模型[13]為

(3)

設(shè)離散非線性隨機(jī)系統(tǒng)(3)中的wk和vk是互不相關(guān)的非零均值高斯白噪聲,且具有時(shí)變統(tǒng)計(jì)特性[14]，即

(4)

式中,Qk是非負(fù)定對稱的過程噪聲的協(xié)方差矩陣;Rk是正定對稱的觀測噪聲的協(xié)方差矩陣。

設(shè)初始狀態(tài)x0與wk、vk互不相關(guān),服從高斯正態(tài)分布,其均值和協(xié)方差矩陣[15]為

(5)

所以要解決的問題就是在量測噪聲和系統(tǒng)噪聲的統(tǒng)計(jì)特性時(shí)變或者未知的情況下,基于角度測量值yk設(shè)計(jì)具有噪聲統(tǒng)計(jì)估計(jì)器的抗差自適應(yīng)UKF算法,計(jì)算目標(biāo)位置速度矢量。

3 帶時(shí)變噪聲統(tǒng)計(jì)特性估計(jì)器的抗差自適應(yīng)UKF算法

定理1根據(jù)極大后驗(yàn)估計(jì)理論[16],在滿足式(4)和式(5)的情況下具有時(shí)變噪聲統(tǒng)計(jì)特性的次優(yōu)極大后驗(yàn)估計(jì)公式為

(6)

J=p[qk-1,Qk-1,rk,Rk,Xk|Yk]=

p[Yk|qk-1,Qk-1,rk,Rk,Xk]·

p[Xk|qk-1,Qk-1,rk,Rk]p[qk-1,Qk-1,rk,Rk]/p[Yk]

(7)

由假設(shè)知時(shí)變噪聲wk、vk相互獨(dú)立,為方便推導(dǎo),不妨?xí)簳r(shí)設(shè)wk、vk具有常值統(tǒng)計(jì)特性:E[wk]=q,Cov[(wk-qk)(wj-qj)T]=Qδkj,E[vk]=r,Cov[(vk-rk)(vj-rj)T]=Rδkj,則根據(jù)條件概率分布密度的乘法定理[17],可得

A1·|P0|-1/2·|Q|-k/2·

(8)

(9)

從而,有

J=A1·A2·|P0|-1/2·|Q|-k/2·|R|-k/2·

{p[qk-1,Qk-1,rk,Rk]/p[Yk]}·

A·|Q|-k/2|R|-k/2·

(10)

式中

A=A1·A2·{p[qk-1,Qk-1,rk,Rk]/p[Yk]}·

為常數(shù)項(xiàng)。

J和lnJ極值相同,式(10)可改寫為

(11)

(12)

(13)

(14)

(15)

(16)

(17)

類似地有

(18)

(19)

證畢

(20)

進(jìn)而有

(21)

(22)

具體計(jì)算步驟如下:

步驟1濾波初始化

步驟2計(jì)算Sigma樣本點(diǎn)

式中,(P)i表示取矩陣P的第i行;

步驟3時(shí)間更新

χi,k|k-1=f(χi,k-1)

步驟4量測更新

yk|k-1=h(χj)

Kk,j=Pxkyk(Pvkvk)-1

步驟5計(jì)算抗差自適應(yīng)因子βk[20]

βk=

步驟6再次進(jìn)行量測更新

4 抗差自適應(yīng)UKF算法流程

圖2 抗差自適應(yīng)UKF算法流程圖Fig.2 Flowchart of robust adaptive UKF algorithm

5 仿真實(shí)驗(yàn)與分析

設(shè)地基觀測站Station的大地坐標(biāo)為:北緯40.1°、東經(jīng)116.2°、高程0.8 km,空間目標(biāo)SatT的軌道根數(shù)如表1所示,星歷時(shí)間為2016年1月1日12:00:00。

表1 空間目標(biāo)的軌道根數(shù)

利用STK模擬生成了作為標(biāo)稱軌道的空間目標(biāo)運(yùn)行軌道。另外,濾波器的軌道預(yù)報(bào)模型和標(biāo)稱軌道計(jì)算模型都只考慮非球形J2、J3、J4項(xiàng)攝動(dòng),以對比檢驗(yàn)本文提出的抗差自適應(yīng)UKF算法的性能。設(shè)仿真時(shí)長為Ta=5 000 s,濾波周期Ts=2 s,N表示濾波點(diǎn)數(shù),蒙特卡羅仿真次數(shù)取M=50次。數(shù)值實(shí)驗(yàn)中,觀測數(shù)據(jù)中加入粗差,噪聲統(tǒng)計(jì)特性q、Q、r、R未知且Q時(shí)變。仿真結(jié)果評價(jià)指標(biāo)采用誤差均方根(root mean square error,RMSE)、平均誤差均方根(average root mean square error,ARMSE)及收斂時(shí)間Tc。

設(shè)初始狀態(tài)誤差為[5 km 5 km 5 km 1 m 1 m 1 m],觀測數(shù)據(jù)誤差取40 μrad。其中,常值系統(tǒng)誤差為10 μrad,隨機(jī)誤差為30 μrad,則有r=[10-510-5],R=diag[9×10-109×10-10];加入的過程噪聲統(tǒng)計(jì)特性為q=[50 50 50 0 0 0],Q=diag[10210210210-610-610-6]。仿真過程中，在t=3 800 s和t=4 000 s時(shí)觀測數(shù)據(jù)中加入400 μrad的粗差,在t=3 200 s后,過程噪聲統(tǒng)計(jì)協(xié)方差陣為Q=diag[1.5×1041.5×1041.5×10410-610-610-6]。

設(shè)置初始觀測噪聲均值和協(xié)方差矩陣分別為:r0=[0 0],R0=diag[10-1010-10],過程噪聲均值和協(xié)方差矩陣分別為:q0=[0 0 0 0 0 0],Q0=diag[1 1 1 10-610-610-6]。

本文設(shè)計(jì)2個(gè)場景。場景1:Station對SatT采用普通UKF算法進(jìn)行跟蹤濾波;場景2:Station對SatT采用抗差自適應(yīng)UKF算法進(jìn)行跟蹤濾波，引入抗差自適應(yīng)因子,初值設(shè)為1,構(gòu)造自適應(yīng)因子采用的常量k0=1.3,k1=6.0,遺忘因子ω=0.95。仿真結(jié)果如圖3～圖6和表2所示。

圖3 位置誤差均方根的比較Fig.3 Comparison of the root mean square of position error

圖4 速度誤差均方根的比較Fig.4 Comparison of the root mean square of velocity error

圖6 抗差自適應(yīng)UKF的自適應(yīng)因子變化曲線(場景2)Fig.6 Adaptive factor change curve of robust adaptive UKF (Case 2)

由以上仿真結(jié)果可知:

(1) 在t=3 200 s前，即過程噪聲的協(xié)方差矩陣還沒有改變前,UKF算法的非線性濾波效果正常,但是在收斂時(shí)RMSE_R和RMSE_V曲線振蕩較大,并且收斂時(shí)間也相對較長，Tc=255 s。當(dāng)過程噪聲的協(xié)方差矩陣改變后,UKF算法逐漸失效,其濾波精度的振蕩降低甚至有發(fā)散的趨勢。特別是在t=3 800 s和t=4 000 s時(shí)觀測數(shù)據(jù)中加入粗差后,振蕩幅值快速增大,RMSE_R最高達(dá)到2 900 m,RMSE_V最高達(dá)到2 m/s。盡管經(jīng)過一段時(shí)間誤差曲線的振蕩幅值逐漸收窄,但是極易受到擾動(dòng)影響而出現(xiàn)發(fā)散。

(2) 抗差自適應(yīng)UKF可以克服過程噪聲協(xié)方差陣時(shí)變和觀測異常的問題,其RMSE_R和RMSE_V的收斂值最終分別保持在940.22 m和0.685 m/s,并且振蕩幅值很小,具有較高的穩(wěn)定性。

(3) 由圖5可知,抗差自適應(yīng)UKF在t=3 200 s后,觀測噪聲的均值的模和協(xié)方差陣的矩陣范數(shù)分別階躍到了0.7×10-4和5×10-8附近振蕩,說明抗差自適應(yīng)UKF將過程噪聲協(xié)方差陣時(shí)變影響通過時(shí)變噪聲統(tǒng)計(jì)估計(jì)器傳遞到了觀測噪聲的統(tǒng)計(jì)特性估計(jì)上,提高了濾波過程中狀態(tài)向量估計(jì)的穩(wěn)定性。

(4) 由圖6可知，抗差自適應(yīng)UKF的自適應(yīng)因子在濾波整個(gè)階段都在頻繁調(diào)整,并且對過程噪聲協(xié)方差陣的時(shí)變性更為敏感。自適應(yīng)因子在t=3 200 s時(shí)降低到0.1,調(diào)整幅值相對最大,這也在一定程度上減輕了時(shí)變噪聲統(tǒng)計(jì)估計(jì)器的工作負(fù)擔(dān)。自適應(yīng)因子在觀測異常時(shí)自身幅值調(diào)整并不明顯,這與觀測粗差不算太“惡劣”有一定關(guān)系,或者說抗差自適應(yīng)UKF的自適應(yīng)因子主要在于整體上提高濾波對觀測異常和狀態(tài)擾動(dòng)的抵制能力。

6 結(jié) 論

傳統(tǒng)UKF算法的濾波精度高于擴(kuò)展卡爾曼濾波算法,但在濾波前需已知噪聲的先驗(yàn)統(tǒng)計(jì)特性[21]。由于在地面站對空間非合作目標(biāo)的光學(xué)測角跟蹤濾波實(shí)際應(yīng)用中,噪聲先驗(yàn)信息一般未知或者不準(zhǔn)確,即使已知但卻有時(shí)變性,會(huì)使傳統(tǒng)UKF算法的濾波精度下降乃至發(fā)散。提出基于極大后驗(yàn)估計(jì)的抗差自適應(yīng)UKF算法,將抗差自適應(yīng)因子和時(shí)變噪聲統(tǒng)計(jì)特性在線估計(jì)相結(jié)合,解決了傳統(tǒng)UKF在噪聲統(tǒng)計(jì)未知時(shí)變下濾波性能不佳的問題, 且有抵制觀測異常和模型擾動(dòng)的能力。仿真結(jié)果表明，該算法可以在未知噪聲先驗(yàn)統(tǒng)計(jì)的情況下,在線估計(jì)未知時(shí)變噪聲的統(tǒng)計(jì)特性,對噪聲變化具有一定的自適應(yīng)能力,使用的抗差自適應(yīng)因子可以控制跟蹤濾波的穩(wěn)定性。

[1] KONG Q, GUO J, SUN Y, et al. Centimeter-level precise orbit determination for the HY-2A satellite using DORIS and SLR tracking data[J]. Acta Geophysica, 2017,65(1):1-12.

[2] JU B, GU D, HERRING T A, et al. Precise orbit and baseline determination for maneuvering low earth orbiters[J]. GPS Solutions, 2017, 21(1):53-64.

[3] ZHANG L, YIN X, NING Z, et al. Robust filtering for a class of networked nonlinear systems with switching communication channels[J].IEEE Trans.on Cybernetics,2017,47(3):671-682.

[4] ABOLHASANI M, RAHMANI M. Robust Kalman filtering for discrete-time systems with stochastic uncertain time-varying parameters[J]. Electronics Letters, 2017, 53(3):146-148.

[5] MAJUMDER R, SADHU S. Robust extended Kalman filter for ballistic object tracking during re-entry[C]∥Proc.of the IEEE India Conference, 2017:1-6.

[6] GAO N, WANG M Y, ZHAO L. A novel robust Kalman filter on AHRS in the magnetic distortion environment[J]. Advances in Space Research, 2017: 201-209.

[7] SONG Q, HAN J D. An adaptive UKF algorithm for the state and parameter estimation of a mobile robot[J]. Acta Automatica Sinica, 2008, 34(1): 72-79.

[8] JULIER S J, UHLMANN J K. Unscented filtering and nonlinear estimation[J]. Proceedings of the IEEE,2004,92(3):401-422.

[9] KHALID S S, REHMAN N U, ABRAR S. Robust stochastic integration filtering for nonlinear systems under multivariate t-distributed uncertainties[J].Signal Processing,2017,140(4):53-59.

[10] BEIDAGHI S, JALALI A A, SEDIGH A K, et al. RobustH∞, filtering for uncertain discrete-time descriptor systems[J]. International Journal of Control Automation and Systems,2017,15(3):995-1002.

[11] HAM B, CHO M, PONCE J. Robust guided image filtering using nonconvex potentials[J].IEEE Trans.on Pattern Analysis and Machine Intelligence,2017:2070-2077.

[12] JIA B, PHAM K D, BLASCH E, et al. Cooperative space object tracking using space-based optical sensors via consensus-based filters[J]. IEEE Trans.on Aerospace and Electronic Systems, 2016, 52(4):1908-1936.

[13] JIA B, BLASCH E, PHAM K D, et al. Space object tracking and maneuver detection via interacting multiple model cubature Kalman filters[C]∥Proc.of the IEEE Aerospace Conference, 2015:1-8.

[14] CHEN C, LI S, QIN H, et al. Real-time and robust object tracking in video via low-rank coherency analysis in feature space[J]. Pattern Recognition, 2015, 48(9):2885-2905.

[15] WEN X, ZHANG X W, ZHU Y P. Design of fault detection observer based on hyper basis function[J]. Tsinghua Science and Technology, 2015, 20(2):200-204.

[16] WAN P, CHEUNG G, FLORENCIO D, et al. Image bit-depth enhancement via maximum-a-posteriori estimation of AC signal[J]. IEEE Trans.on Image Processing, 2016, 25(6):2896-2909.

[17] YU T, LAI Y H. Generalized maximum a posteriori, spectral amplitude estimation for speech enhancement[J]. Speech Communication, 2016, 76(C):112-126.

[18] QUINTANA D, GARCARODRGUEZ S, CINCOTTI S, et al. Combining RMT-based filtering with time-stamped resampling for robust optimization[J]. International Journal of Computational Intelligence Systems, 2015, 8(5):874-885.

[19] 趙琳, 王小旭, 孫明, 等. 基于極大后驗(yàn)估計(jì)和指數(shù)加權(quán)的自適應(yīng)UKF濾波算法[J]. 自動(dòng)化學(xué)報(bào), 2010,36(7):57-61.

ZHAO L, WANG X X, SUN M, et al. Adaptive UKF filtering algorithm based on maximum a posterior estimation and exponential weighting[J]. Automatica Sinica,2010,36(7):57-61.

[20] SNIEKERS S, VAART A V D. Adaptive Bayesian credible sets in regression with a Gaussian process prior[J].Statistics,2015,9(2): 1504-1512.

[21] 劉光明. 基于天基測角信息的空間非合作目標(biāo)跟蹤算法及相關(guān)技術(shù)研究[D]. 長沙:國防科技大學(xué), 2011.

LIU G M. Non-cooperative space target tracking algorithm based on angular measurement information and related technology[D].Changsha: National University of Defense Technology,2011.