空間機動目標跟蹤的改進自適應IMM算法

尹聚祺, 楊 震, 羅亞中, 周劍勇

(國防科技大學空天科學學院, 湖南 長沙 410073)

0 引 言

隨著航天技術的發展,在軌空間目標數量急劇增長,使得太空環境越來越擁擠[1-3]。未來,具有自主機動能力的空間非合作目標對己方航天器的自主異常接近將對己方航天器安全產生重大威脅,天基在軌自主空間態勢感知能力對己方航天器在軌安全防御顯得至關重要[4-6]。通過軌道機動,空間目標可按需要逼近其他目標。一般而言,空間目標的機動方式可建模為兩種[7]:① 脈沖機動;② 有限推力機動。本研究重點關注前者,即空間非合作目標在脈沖機動控制下的狀態估計問題。

近年來,空間非合作機動目標跟蹤問題受到廣泛關注,常用研究方法包括單模型方法和多模型方法。由于非合作性,目標機動發生時刻和脈沖推力大小均為未知量,傳統單模型濾波方法往往在機動發生時由于模型失配造成濾波發散或精度降低。Fitzgerald[8]分析了當被跟蹤目標具有機動屬性時,由于標稱動力學與實際機動不匹配,傳統卡爾曼及擴展卡爾曼濾波(extended Kalman filter,EKF)算法的跟蹤性能下降甚至發散。Zhou等[9]基于殘差序列正交化的思想提出強跟蹤濾波(residual-normalized strong tracking filter,RNSTF)算法,提高了單模型濾波算法對機動信息的魯棒適應能力,但對脈沖機動敏感性不足而易導致濾波精度下降。Jiang等[10]針對這一問題進一步提出殘差歸一化的RNSTF方法,能更及時地檢測機動并提高機動后的跟蹤精度,但整體跟蹤精度仍有待提到。由于單模型方法不能有效適應目標運動過程中的各種運動狀態,在非合作機動目標跟蹤問題上應用受限,因此發展出了多模型方法。

最著名的多模型方法是Blom等[11]提出的交互多模型(interacting multiple model,IMM)算法。IMM算法通過選擇恰當的模型集來匹配目標運動過程中的狀態變化,通過模型轉移概率將各模型輸出結果進行交互,在機動目標跟蹤方面具有適應能力強的優勢,被廣泛應用到機動目標跟蹤問題中[12-18]。在空間機動目標跟蹤領域,Lee等[19]基于IMM算法給出一種空間目標機動檢測和機動后軌道表征算法。Goff等基于IMM算法和協方差膨脹思想來解決脈沖機動目標狀態跟蹤問題,提出一種非合作機動目標實時跟蹤算法[20],但存在漏檢情況。針對連續推力機動目標跟蹤問題,Goff等將IMM算法與可變狀態維度濾波器相結合,提出一種自適應變維濾波估計算法[21]。Xiong等[22]將魯棒卡爾曼濾波與IMM算法相結合,提出一種魯棒多模自適應算法。Huang等[23]基于IMM估計算法,結合Singer模型及無跡卡爾曼濾波方法,給出了超音速滑翔目標跟蹤方法。以上文獻進行了富有價值的研究,但基于傳統IMM算法設計的估計方法為提高算法適應能力存在模型集設計復雜,且估計精度不高的問題。

在傳統IMM方法中,模型轉移概率參數一般設計為先驗給定的固定值。由于目標運動的不確定性,該固定先驗假設的模型概率轉換方式是模型切換與未切換情況下的折衷[24]。模型轉移概率的先驗不確定性會造成濾波精度的損失,因此對其進行改進的在線自適應性估計方法一直是國內外的研究熱點。Lee等[25]通過將未知機動信息建模為特定條件下狀態變化問題,從而使得模型轉移概率自適應變化,然而這仍需要一定的先驗信息。國內學者針對該問題進行了詳細研究[26-30]。郭志[29]和許登榮[30]分別從模型概率和模型似然函數的角度出發,分別將平方根容積卡爾曼濾波器(square-root cubature Kalman filter,SRCKF)、強跟蹤修正輸入估計(strong tracking modified input estimation,STMIE)及勻速運動(constant velocity,CV)模型引入IMM算法中,提出了IMM-SRCKF及IMM-STMIECV算法對模型轉移概率進行自適應性改造,改善了模型切換速度和濾波精度。模型似然函數值從量的角度衡量了該模型的匹配情況,但簡單基于模型似然函數比來對模型轉移概率進行自適應修改在目標機動時刻存在奇異問題。

本文基于IMM算法,并結合文獻[30]模型轉移概率自適應算法,研究了追蹤星(空間非合作目標)采用脈沖機動方式接近目標星(己方航天器),目標星采用測角加測距的量測方式跟蹤追蹤星,在雙方信息不透明情況下追蹤星的相對狀態估計問題。本文的創新點在于提出了一種改進自適應IMM-EKF估計方法,具體為:① 將白噪聲模型與EKF算法相結合,設計了精簡的模型集以適應目標機動和非機動時的狀態變化;② 航天器相對運動采用直接離散化處理的非線性相對軌道動力學方程,相比于近距離、圓參考軌道及線性化假設下的Clohessy-Wiltshire(CW)方程提高了預報精度且增加了適用范圍;③ 提出了一種變換函數進而克服了簡單使用模型似然函數比對模型轉移概率修正時存在奇異的問題。最后對比分析了加入速率測量信息對濾波結果的改善作用。仿真結果表明，本文提出的改進自適應IMM-EKF方法相比傳統IMM濾波方法、IMM-SRCKF、RNSTF濾波及容積卡爾曼濾波(cubature Kalman filter, CKF)方法跟蹤效果更好,引入速率量測信息后,IMM-EKF方法最大峰值誤差及估計精度得到了改善。

1 相對狀態預測模型

考慮離散時間非線性動力學及觀測方程:

(1)

式中:X(k)∈Rn為n維為狀態向量;y(k)∈Rp為p維觀測向量;f(X(k)):Rn→Rn為一步狀態預測方程;h(X(k)):Rn→Rp為觀測方程;w(k)及v(k)分別為n維、p維零均值高斯白噪聲,滿足E[w(k)w(k)T]=Q(k),E[v(k)v(k)T]=R(k),且B∈Rn×l,G∈Rn×m均為常矩陣;δt,k+1為Kronecker記號函數,當t=k+1時,δt,k+1=1,t≠k+1時,δt,k+1=0;U為追蹤星控制輸入量,δt,k+1GU表示追蹤星大小及時間均未知的機動信息。

以目標星質心為原點建立VVLH(vehicle velocity, local horizontal)坐標系,其中x軸指向目標速度方向,z軸指向地心,y軸滿足右手系定則。在VVLH系下建立相對運動方程來描述追蹤星和目標星的相對運動關系,則脈沖機動下的相對狀態預報方程可表示為

(2)

式中:

ΔT為采樣時間間隔;α為采樣時刻目標星軌道角加速度;ω為采樣時刻目標星軌道角速度;n為目標星平均軌道角速度;e為目標星軌道偏心率;θ為采樣時刻目標星真近點角;Rc為采樣時刻目標星軌道半徑;μ為地球引力常數。注意到,式(2)為對二體相對非線性動力學方程的直接離散化處理,相比于近距離、圓參考軌道及線性化假設下的CW方程提高了預報精度且增加了適用范圍,即目標星軌道滿足橢圓情況。

2 相對狀態觀測模型

以運行在橢圓軌道上的目標星質心為原點建立VVLH坐標系,目標星利用自身攜帶設備對追蹤星進行測量跟蹤,測量量為相對距離、俯仰角和方位角,建立天基觀測模型。在VVLH系下,如圖1所示。俯仰角α定義為追蹤星相對目標星的視線方向與xy平面的夾角,取值范圍為(-π/2,π/2);方位角β定義為視線方向在xy平面上的投影與x軸方向的夾角,以x軸為起點逆時針轉動為正,取值范圍為(-π,π)，觀測方程為

(3)

圖1 觀測坐標系Fig.1 Observation coordinate system

對應觀測敏感性矩陣為

(4)

H各分量為

引入相對速率量測信息后,相應的觀測方程及觀測敏感矩陣分別為

(5)

(6)

式中:H各非0分量為

3 改進自適應IMM估計算法

當目標由非機動狀態轉入機動狀態時,傳統單模型算法一般會由于模型失配而導致濾波發散。IMM估計算法可包含多個模型以適應不同工況,通過交互融合多種模型濾波結果,在保持濾波穩定性的同時達到精度要求。

3.1 標準IMM算法

標準IMM算法包含模型輸入交互、模型估計、模型概率更新和綜合輸出4個遞推模塊,下面給出從k-1時刻到k時刻的遞推步驟。本研究中模型集均采用白噪聲模型,濾波器為兩個EKF器分別記為EKF1和EKF2。其中EKF1過程噪聲Q1較小,為非機動模型;EKF2將未知的機動信息建模為較大的過程噪聲Q2。

步驟 1模型輸入交互

(7)

(8)

式中:mi|j(k-1)為模型i到模型j的混合概率,可由模型概率mi(k-1)和模型轉移概率γij表示:

(9)

模型概率mi初值及模型轉移概率γij為經驗給定值。

步驟 2模型估計

(10)

Pj(k|k-1)=F(k-1)P0j(k-1)F(k-1)T+
BQj(k-1)BT

(11)

(12)

(13)

(14)

Pj(k)=[I-K(k)H(k)]Pj(k|k-1)

(15)

K(k)=Pj(k|k-1)H(k)T·
[H(k)Pj(k|k-1)H(k)T+R(k)]-1

(16)

(17)

步驟 3模型概率更新

假設模型j的殘差vj(k)服從高斯分布,其協方差為Sj(k),則模型似然函數Λj(k)和模型概率mj(k)可更新為

(18)

(19)

模型殘差vj(k)及其協方差Sj(k)可表示為

(20)

Sj(k)=H(k)Pj(k|k-1)H(k)T+R(k)

(21)

步驟 4綜合輸出

(22)

(23)

3.2 模型轉移概率自適應算法

(24)

(25)

(26)

否則保持其所在行的元素不變。

3.3 改進自適應IMM算法

將第3.2節方法應用到脈沖機動目標的相對狀態跟蹤問題中,當采用非機動模型EKF1濾波器跟蹤目標的運動時,非機動時刻具有較高精度的跟蹤效果;而當目標機動時,非機動模型不能匹配目標的狀態變化。這時非機動模型似然函數值接近0,采用似然函數比放大或縮小上一時刻的模型轉移概率會出現奇異現象,文獻[30]所提IMM模型狀態轉移概率自適應修正方法失效。

為解決此問題,本文在第3.2節的基礎上提出非線性變換函數對模型似然函數進行變換。變換前,傳統模型似然函數有以下不足:

(1)非機動模型似然函數值在機動時刻數值趨于0,導致模型似然函數比奇異;

(2)模型似然函數值在整個濾波過程中數值較大(本文算例為104～105量級),在設計可避免奇異的變換函數時,需考慮其影響。

為解決因非機動模型似然函數值在機動時刻數值趨于0而造成模型似然函數比在機動時刻奇異問題,本文采用在[0,+∞)區間上函數值不過0的自然指數函數(e)對其進行變換。然而直接采用指數函數進行變換會造成變換后的模型似然函數值無窮大的情況,進而引起數值計算錯誤,所以對變換函數進行如下設計。

首先，對模型似然函數進行冪函數縮放變換:

Λj1=(Λj0)η,0<η<ηmax

(27)

式中:η為縮放因子,通過選擇η的大小來縮放初始模型似然函數值至一合理范圍內,避免進行e指數函數變換后Λj值過大的情況。

其次，對縮放變換后的模型似然函數進行指數函數變換。

Λj=exp(Λj1)

(28)

綜合式(27)和式(28),ηmax可按以下步驟確定。

步驟 1令η=0,則IMM-EKF方法退化為傳統IMM方法,可確定變換前模型似然函數值量級范圍Λscale。進而由計算機浮點數運算上界Valuemax確定關系方程。

exp((Λscale)η)≤Valuemax

(29)

步驟 2由式(29)求解η范圍:

(30)

本文算例中Λscale取保守值106,64位計算機浮點數運算上界Valuemax約為1.797 6×10308,則可得ηmax≈0.475 2。不同η的取值將影響濾波的估計精度及機動后的收斂速度。

由于以e為底的指數函數和冪指數大于0的冪函數都是增函數,故其復合函數也是增函數,所以可以通過選擇參數縮放因子η的值,在不改變似然函數變化趨勢的情況下使得變換后的模型似然函數值Λj≥1,避免采用第3.2節方法對模型轉移概率進行自適應修正時出現奇異現象。

圖2給出了研究的模型轉移概率自適應IMM算法框圖。實線表示k-1時刻算法計算流程,虛線表示由k-1時刻更新至k時刻,進行下一步計算。

圖2 模型轉移概率自適應修正算法框圖Fig.2 Diagram of adaptive correction algorithm for modeltransition probability

4 仿真算例分析

為驗證改進方法的有效性,首先仿真分析了所提變換函數的有效性;其次,將與傳統IMM算法、IMM-SRCKF算法、RNSTF濾波及CKF濾波方法對比分析,說明改進后的IMM算法在跟蹤效果上的優勢;最后,仿真對比分析引入速率測量信息對濾波延遲及最大峰值誤差的改善效果。

4.1 仿真場景及參數設置

假設目標星運行在地球同步軌道上,追蹤星在脈沖機動控制下接近目標星,目標星采用光學相機和激光測距儀對追蹤星進行觀測測量,實時獲得相對距離、俯仰角和偏航角參數,估計兩者間相對狀態,如圖3所示。

圖3 仿真場景示意Fig.3 Simulation scenario

設激光測距儀測距誤差標準差為1 m,光學相機角度測量誤差標準差為0.001 rad,測量頻率為1 Hz。目標星初始軌道根數如表1所示。在以目標星質心為原點的VVLH系下,追蹤星相對目標星X(0)=[-60 000 m,5 000 m,20 000 m, 30 m/s,-5 m/s,-5 m/s]T,且在仿真過程中,假設追蹤星在t=300 s和t=800 s時作兩次機動,脈沖為[10 m/s,0 m/s,5 m/s]T、[0 m/s,2.8 m/s,-9.8 m/s]T。

表1 目標星初始軌道根數

多模型方法模型轉移概率矩陣設為

各模型j(j=1,2) 初始模型概率、濾波估計初值、初始協方差矩陣及過程噪聲協方差矩陣設置為

mj(0)=0.5

Pj(0)=diag(10, 10, 10, 1, 1, 1)

Q1(k)=diag(10-4, 10-4, 10-4, 10-8, 10-8, 10-8)

Q2(k)=diag(10-2, 10-2, 10-2, 102, 102, 102)

模型轉移概率修正參數分別設置為γ=1/10,Th=0.85,η=1/10。

文獻[30]所提RNSTF濾波方法仿真場景參數設置同IMM濾波方法一致,漸消因子設為ρ=0.95,CKF仿真參數同上。

4.2 仿真結果分析

對以上方法進行Monte Carlo仿真分析,結果如圖4和圖5所示。

圖4 變換前模型似然函數變化情況Fig.4 Change of model likelihood function before transforming

圖5 變換后模型似然函數變化情況Fig.5 Change of model likelihood function after transforming

圖4和圖5分別給出了對IMM算法模型集中各模型似然函數進行變換前后的仿真結果。變換前,各模型似然函數值波動劇烈,波動范圍在104～105量級,且模型1似然函數值在追蹤星機動后趨于0,這使利用模型似然函數比對模型轉移概率自適應修正時出現奇異的現象,導致濾波發散。采用本文所提方法對模型似然函數值進行變換后,各模型似然函數值波動范圍大大減小,且在機動后最小值趨于1,避免了造成奇異現象,提高了算法穩定性。

表2給出了變換函數中參數縮放因子η取不同值時對所提IMM-EKF算法跟蹤性能的影響統計情況,并對每類工況進行100次Monte Carlo仿真分析，統計均方根(root mean square, RMS)估計誤差均值。表2中“—” 表示濾波不成功。例如,當η=1時,即變換函數為簡單e指數變換,這使得變換后的模型似然函數無窮大,進而導致濾波失敗。當η的值不斷取小時,RMS位置及速度誤差均值不斷增大,但機動后的收斂速度不斷加快。例如,當η值趨于0時,變換后的模型似然函數值趨于1,使得似然函數比接近1,基于似然函數比對模型轉移概率矩陣不再具有修正作用,IMM-EKF濾波效果退化接近標準IMM算法。綜合考慮濾波估計誤差和機動后的收斂速度,本文將η的值取為0.1。

表2 不同η取值對算法跟蹤性能影響的統計情況

圖6給出了未考慮速率測量的多模型濾波RMS誤差仿真結果。由圖6分析可得,隨著仿真時間推進,即兩航天器間的相對距離減小,濾波估計結果的RMS位置誤差呈減小趨勢。相比于多模型方法,單模型RNSTF方法在全濾波時段位置和速度估計精度較低,且在追蹤星機動后產生較大的峰值誤差,但具有較快的收斂速度。CKF算法在全濾波階段位置及速度估計精度較低,且對機動信息不敏感。IMM-SRCKF算法在本算例工況下無論是濾波精度還是機動發生后的收斂速度都較標準IMM算法所得結果相差不大。較于IMM、IMM-SRCKF、IMM-EKF,本文改進的模型轉移概率自適應IMM-EKF算法在估計精度上具有更好效果。

圖6 未考慮速率測量的RMS誤差對比分析Fig.6 Comparative analysis of RMS error without considering rate measurement

表3給出了比較的4種方法在全觀測時段內RMS位置及速度平均誤差,最大峰值誤差及濾波結束時刻估計誤差統計數據,從量化的角度說明了以上分析結論。就對比的3類指標而言,單模型RNSTF方法及CKF算法估計均具有較大的估計誤差,IMM-EKF方法估計誤差最小,表明了提出的IMM-EKF方法良好的估計效果。

表3 未考慮速率測量的算法跟蹤性能比較

圖7給出了IMM-EKF方法模型概率值變化情況。當追蹤星未機動時,非機動模型EKF1占據主要地位,對應模型概率接近數值1,非機動模型EKF2模型概率接近數值0;當追蹤星機動時,模型集中非機動模型EKF1預測值和測量值產生較大偏差,模型似然函數值減小,對應模型概率發生突變減小,模型概率接近數值0,機動模型EKF2模型概率接近數值1,模型切換。因此,IMM-EKF方法中模型概率值在機動發生后突變這一規律可用于目標機動檢測的判斷依據,實現跟蹤過程中實時的機動檢測,這對天基態勢感知十分必要。

圖7 IMM-EKF方法模型概率變化情況Fig.7 Change of model probability of IMM-EKF

圖8給出了考慮速率測量的多模型濾波仿真結果。引入速率測量信息,濾波開始后的RMS位置及速度估計誤差明顯減小,且改善了最大峰值誤差,機動后的收斂速度也有改善。當目標速度發生變化時,由于距離及角度測量量非瞬態響應量,目標的機動信息不能在機動時刻體現,故存在檢測延遲,進而導致較大的峰值誤差。對于多模型濾波,加入速率量測數據后,機動時刻,濾波器殘差敏感到速率的明顯變化后,EKF1模型概率由1變為0,此刻EKF2的濾波結果主導綜合輸出結果,由于模型不匹配造成的RMS速度誤差的峰值誤差較大的現象得到抑制,這種變化在發生較大機動量時更能體現。表4定量給出了速率測量信息對所提IMM-EKF方法跟蹤性能的影響情況。引入速率測量信息后,IMM-EKF方法RMS位置及速度估計平均誤差分別減少27.67%、60.07%,最大峰值誤差分別減少7.27%、11.71%,跟蹤能力明顯提升。

圖8 考慮速率測量的多模型RMS誤差對比分析Fig.8 Comparative analysis of multi-model RMS error considering rate measurement

表4 速率測量對IMM-EKF方法跟蹤性能的影響

5 結 論

空間博弈對抗條件下,對機動目標實時精準的相對狀態估計是制定博弈策略的重要先決條件。本文在空間近距離追逃博弈背景下,基于IMM方法研究了目標星對脈沖機動控制下的追蹤星進行相對狀態估計的問題,提出了模型轉移概率改進自適應IMM算法,克服了模型似然函數比在機動時刻奇異的問題。仿真結果表明:① 改進的模型轉移概率自適應IMM算法相比于傳統IMM、IMM-SRCKF、單模型RNSTF及CKF濾波方法,在位置和速度估計精度上明顯提高;② 當追蹤星機動時,模型集中非機動模型預測值和測量值產生較大偏差,模型似然函數值減小,對應模型概率發生突變減小,這一規律可作為目標機動檢測依據,可在空間目標態勢感知任務中用于實時機動檢測與告警;③ 在本文提出的模型轉移概率自適應IMM算法框架下,進一步引入速率量測信息,有效抑制濾波延遲問題,減小了多模型濾波估計方法的最大峰值誤差,同時提高了跟蹤精度,可見增加速率測量器件對提高機動目標跟蹤精度具有重要意義。