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

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

0 引 言

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

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

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

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

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

1 相對狀態(tài)預測模型

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

(1)

式中:X(k)∈Rn為n維為狀態(tài)向量;y(k)∈Rp為p維觀測向量;f(X(k)):Rn→Rn為一步狀態(tài)預測方程;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記號函數(shù),當t=k+1時,δt,k+1=1,t≠k+1時,δt,k+1=0;U為追蹤星控制輸入量,δt,k+1GU表示追蹤星大小及時間均未知的機動信息。

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

(2)

式中:

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

2 相對狀態(tài)觀測模型

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

(3)

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

對應觀測敏感性矩陣為

(4)

H各分量為

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

(5)

(6)

式中:H各非0分量為

3 改進自適應IMM估計算法

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

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)和模型轉(zhuǎn)移概率γij表示:

(9)

模型概率mi初值及模型轉(zhuǎn)移概率γij為經(jīng)驗給定值。

步驟 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模型概率更新

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

(18)

(19)

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

(20)

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

(21)

步驟 4綜合輸出

(22)

(23)

3.2 模型轉(zhuǎn)移概率自適應算法

(24)

(25)

(26)

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

3.3 改進自適應IMM算法

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

為解決此問題,本文在第3.2節(jié)的基礎(chǔ)上提出非線性變換函數(shù)對模型似然函數(shù)進行變換。變換前,傳統(tǒng)模型似然函數(shù)有以下不足:

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

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

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

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

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

(27)

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

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

Λj=exp(Λj1)

(28)

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

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

exp((Λscale)η)≤Valuemax

(29)

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

(30)

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

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

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

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

4 仿真算例分析

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

4.1 仿真場景及參數(shù)設(shè)置

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

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

設(shè)激光測距儀測距誤差標準差為1 m,光學相機角度測量誤差標準差為0.001 rad,測量頻率為1 Hz。目標星初始軌道根數(shù)如表1所示。在以目標星質(zhì)心為原點的VVLH系下,追蹤星相對目標星X(0)=[-60 000 m,5 000 m,20 000 m, 30 m/s,-5 m/s,-5 m/s]T,且在仿真過程中,假設(shè)追蹤星在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 目標星初始軌道根數(shù)

多模型方法模型轉(zhuǎn)移概率矩陣設(shè)為

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

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)

模型轉(zhuǎn)移概率修正參數(shù)分別設(shè)置為γ=1/10,Th=0.85,η=1/10。

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

4.2 仿真結(jié)果分析

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

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

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

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

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

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

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

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

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

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

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

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

圖8給出了考慮速率測量的多模型濾波仿真結(jié)果。引入速率測量信息,濾波開始后的RMS位置及速度估計誤差明顯減小,且改善了最大峰值誤差,機動后的收斂速度也有改善。當目標速度發(fā)生變化時,由于距離及角度測量量非瞬態(tài)響應量,目標的機動信息不能在機動時刻體現(xiàn),故存在檢測延遲,進而導致較大的峰值誤差。對于多模型濾波,加入速率量測數(shù)據(jù)后,機動時刻,濾波器殘差敏感到速率的明顯變化后,EKF1模型概率由1變?yōu)?,此刻EKF2的濾波結(jié)果主導綜合輸出結(jié)果,由于模型不匹配造成的RMS速度誤差的峰值誤差較大的現(xiàn)象得到抑制,這種變化在發(fā)生較大機動量時更能體現(xiàn)。表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 結(jié) 論

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