郁 凱,王驤予涵,盧鴻謙
(1. 哈爾濱工業大學航天學院,哈爾濱 150000;2. 哈爾濱工程大學航天與建筑工程學院,哈爾濱 150001)
隨著技術的不斷更新,現代飛行器的速度顯著提高。因此,以往將飛行器的運動狀態簡化為勻速運動的做法已經不再具有代表性。如果地面設備用于觀測飛行器,考慮到駕駛員的操控行為和控制指令,目標可能隨時進行轉彎、閃避或采取其他特殊的攻擊姿態等機動動作[1]。通過使用單個或多個主動或被動站點,探測目標(散射體或輻射體)并獲取時間信息TOA(time of arrival)以及角度信息AOA(angle of arrival)等數據[2],然后結合適當的數據處理方法(例如三角定位、時差定位、頻差定位等),可以確定機動目標在三維空間中的位置[3]。這是探測系統用于對機動目標進行定位和跟蹤的基本過程。
最初的卡爾曼濾波適用于線性系統,利用之前時刻的狀態對當前時刻的狀態進行估計,并用觀測量對估計狀態進行修正。這種算法在高斯條件下不僅滿足極大似然估計,還是最小方差估計的最優結果。目前,對于機動目標定位算法研究里面最普遍的方法仍是基于上述方法及其優化方法[4]。優化方法主要包括擴展卡爾曼濾波器(extended Kalman filter,EKF)、粒子濾波器等[5]。羅笑冰等[6]利用雷達提供的量測信息求取噪聲參數,進而利用噪聲信息優化Singer模型加速度先驗信息,提高了該模型的適用性范圍。在此基礎上為了更好地融合多站提供的量測信息,又發展出了多模型算法[7],這種算法在每次迭代濾波過程中充分考慮其他各個濾波器模型估計結果,引入基于貝葉斯估計的交互理念,實現更高精度的目標跟蹤。典型的有交互式多模型算法 (interacting multiple models,IMM)。該方法可以根據存儲器靈活選取馬爾可夫階數,完成更高精度的運動建模與濾波估計融合[8]。臧勤等[9]為了解決對地目標定位過程中出現的不連續收斂現象,將時間信息引入交互式多模型,但是這種方法對時間噪聲較為敏感,適用于遮擋物較少的空間跟蹤。Phani和Arye[10]從IMM貝葉斯原理出發,使用順序貝葉斯濾波完成跟蹤過程,該算法極大降低了原有濾波方法的計算復雜度和觀測空間復雜度,但是也舍棄了一定精度。為了削弱單站定位的故障性影響,提高系統可靠性和穩定性,葉瑾等[11]首先對傳感器的量測信息進行精度方面的加權融合,然后使用對非線性量測關系一階泰勒展開,在交互式多模型的框架下完成了初步的信息融合跟蹤,但是該算法子濾波器運動模型單一,在目標做逃避飛行等運動時,精度降低。
上述方法對于機動目標的跟蹤大多基于單個運動模型或者機動性不強的運動模型組來計算,無法降低模型失配帶來的跟蹤損失[12]。為此,本文首先將勻速直線運動、勻加速直線運動、慢機動Singer模型和快機動Singer模型引入交互式多模型來預測目標機動過程。其中勻速直線運動和勻加速直線運動適用于目標加速度穩定變化的目標運動,快Singer和慢Singer模型則從時間的相關領域模擬目標加速度機動程度。然后利用TOA-AOA幾何關系求解的目標位置作為濾波初始值狀態,結合擴展卡爾曼子濾波器構建交互式多模型算法。由于單純利用TOA-AOA幾何關系就可以求解得到目標位置,所以相較于其他濾波隨機生成的高斯初值模型,本文濾波算法從一開始就可以收斂,且收斂穩定。最后,本文算法在目標機動情況下,通過交互式多模型框架切換模型子濾波器應對目標非線性狀態轉移,進而有效跟蹤目標。
假設主動站O的坐標為(0,0,0),發出電磁波的時間為T0,(xA,yA,zA)為被動站A的坐標,目標M的坐標為(x,y,z) ,A對M的觀測量為(φ,ψ,TA)′,TA為被動站A接收波的時間,φ為被動站A接收波方位角,ψ為被動站A接收波俯仰角,量測模型如下
(1)
(2)
(TA-T0)·c=
(3)
其中c為光速。

Xk+1=φXk+W
(4)
式中,k代表時刻,φ代表狀態轉移矩陣,W為系統噪聲。
觀測方程為

(5)
(6)

觀測方程為非線性方程,根據泰勒展開式保留一階項的方法求得其擴展卡爾曼濾波量測矩陣為
(7)

本章算法以TOA-AOA幾何關系求解得到的目標位置作為勻速直線運動擴展卡爾曼子濾波器、勻加速直線運動擴展卡爾曼子濾波器、快Singer模型擴展卡爾曼子濾波器和慢Singer模型擴展卡爾曼子濾波器的初始值,直接保證了算法收斂性;接著將各個模型子濾波器的狀態估計作為交互式多模型的輸入,根據高斯似然比分配各個子濾波器權值比重,最后加權得到下一時刻的目標位置估計。通過模型切換子濾波器組合的非線性關系,從加速度層面實時模擬目標可能出現的任何機動運動。
在目標機動運動時,使用單一加速度建模方式無法避免模型失配帶來的損失。動態多模型算法可以緩解該矛盾,其中突出的算法就是IMM算法。IMM法來源于廣義貝葉斯估計,在各類濾波器模型的基礎上,根據高斯置信度模型實時選擇最優機動濾波器模型[13]。
若量測噪聲滿足高斯假設,則根據卡爾曼濾波的量測模型,可確定兩隨機變量量測參數Z、狀態參數X之間的關系,可由如下概率密度函數表示
(8)
其中n代表量測參數Z的維度。

1)狀態估計的交互處理:



(9)
其中
Poj(k-1|k-1)=
(10)
2)模型修正:

(11)
Poj(k|k-1)=φPoj(k-1|k-1)φT+
Γ(k-1)Q(k-1)Γ(k-1)T
(12)
Kj(k)=Poj(k|k-1)Hj(k)T·
(Hj(k)Poj(k|k-1)Hj(k)T)-1
(13)
Hj(k)Xoj(k|k-1))
(14)
Pj(k|k)=(I-Kj(k)Hj(k))Poj(k|k-1)
(15)
3)模型可能性計算:

(16)
(17)
(18)
4)模型概率更新:
(19)
(20)

(21)
(22)

2.2.1 離散勻速直線運動模型
目標運動過程一般是平穩運動和機動運動的結合,平穩運動狀態可以用離散勻速(constant velocity, CV)模型來表述,下文叫CV模型。
其狀態方程可以寫為
XT,k+1=φXT,k+Γvk
(23)

其中δT為采樣步長,本文統一為1 s,v1k,v2k,v3k為相互獨立的高斯白噪聲。
因此狀態噪聲的方差陣為
(24)
記

2.2.2 離散勻加速直線運動模型
對于大部分目標的不劇烈機動運動狀態建模,為了計算方便,常以用離散勻加速(constant acceleration, CA)模型來表述,也就是CA模型。
它的狀態方程為
XT,k+1=φXT,k+Γvk
(25)
式中

vk[v1k,v2k,v3k]T各元素均服從高斯噪聲,它們代表目標三軸運動的加速度變化量。于是狀態噪聲協方差陣為
(26)

2.2.3 Singer運動模型
Singer運動模型最早來源于全局統計模型[14],該模型根據大量實測機動數據和領域加速度技術歸納而來,基本滿足各類運動模式。但是這種算法對于經驗參數要求比較高,目標機動程度不一樣,Singer參數也需要變動。為此,設計快機動Singer模型和慢機動Singer模型。
設定目標機動加速度為α(t), 其相關函數為
(27)

(28)
式中,P0為非機動概率;典型的經驗取值范圍[15]:大氣擾動κ1=1,慢速轉彎κ2=1/60,逃避機動κ3=1/20。大氣擾動參數通常指的是用于描述大氣環境對目標運動的不確定性和擾動的參數。大氣擾動參數可以包括風速、風向、氣壓變化等與大氣環境相關的因素。這些因素會對目標的運動軌跡產生影響,導致目標的實際運動與理想模型存在偏差。
由于a(t)是色噪聲,因此可以用白噪聲w(t)通過一個成形濾波器[15]來產生,產生的關系為
(29)
式中w(t)的相關函數為
(30)
目標連續運動狀態模型為
(31)
將其離散化,該模型可以變為
XδT,k+1=φ(δT,α)XδT,k+Uk
(32)
φ(T,α)=L-1[(sI-N)-1]
=Fj
(33)

因此,噪聲Uk的協方差矩陣為
(34)

Qk的精確表達式(Qk為對稱陣)為
(35)
交互式多模型濾波和擴展卡爾曼濾波算法的收斂速度和收斂穩定性與目標位置初始有關系。可以利用TOA-AOA幾何關系,根據量測信息求取目標位置解析解,以此作為濾波算法狀態迭代的初值。目標M與主動站O、被動站A的空間位置如圖1所示。

圖1 目標解算模型Fig.1 Target solution model

(36)
如果存在β=π-u的幾何關系,那么
f=LOA·cos(β)
(37)
b=LOA·sin(β)
(38)
(39)

M的坐標,即各個子濾波器濾波算法位置初值為
(40)
速度初值和加速度初值統一設置為0。
根據主-被動結合定位方法的基本原理,針對不同機動特性目標進行運動估計。在對機動目標進行跟蹤時,IMM算法可以在各類運動模型的擴展卡爾曼子濾波器中實時轉換,從而實現更高精度的目標定位跟蹤。但是,在一般IMM算法中,各個模型匹配的都是CV和CA模型的擴展卡爾曼濾波器。當目標發生強機動時,跟蹤精度迅速降低,嚴重影響跟蹤效果。本章針對目標飛行完整過程(包含起飛與降落)將融合勻速直線運動、勻加速直線運動、慢機動Singer1模型和快機動Singer2模型的IMM算法,簡稱IMM2,與單個EKF狀態運動模型進行比對仿真。
設置主動站O位于(0,0,0),被動站A位于(-10,0,0) km,被動站B位于(0,10,0) km,被動站C位于(10,0,0) km,被動站D位于(0,-10,0) km,仿真時間長2 600 s,采樣時間間隔2 s,蒙特卡羅仿真300次。目標最大水平速度60 m/s,最大爬升率7 m/s,最大下降率8 m/s,最大使用高度5 km,最大俯仰角±50°,初始位置(-562,-153,-421) m,表1是目標部分機動時間表。

表1 目標部分機動時間表
設定IMM子濾波器初始參數:Singer1模型的xyz三軸機動頻率分別為1/60,1/60,1,單位:Hz。三軸機動加速度的極大值分別為8 m/s2,5 m/s2,3 m/s2,三軸最大機動概率為0.1,0.1,0.1,三軸非機動概率為0.5,0.5,0.5;Singer2模型的xyz三軸機動頻率分別為1/20,1/20,1/60,單位:Hz,三軸機動加速度的極大值分別為8 m/s2,8 m/s2,3 m/s2,三軸最大機動概率為0.2,0.3,0.1,三軸非機動概率為0.2,0.1,0.1;CA模型的三軸狀態標準差分別為6 m/s2,4 m/s2,2 m/s2;CV模型的三軸狀態標準差分別為5 m/s2,3 m/s2,2 m/s2。
設置測角誤差為1°,時間誤差為20 ns,模型互相關轉移矩陣、模型置信度分別如下
表2所示為各類算法從0時刻開始到達該時間段內的位置定位誤差均值。對比仿真圖2、圖3和表2可以發現,IMM2模型在機動時刻的定位跟蹤效果更好,但是存在收斂速度較慢的問題,這是由于前半段時間,目標上升運動導致俯仰角變化太大。此外,不同濾波器參數的設置也會影響濾波效果。可以看出來,針對該軌跡Singer1模型的機動參數設置比Singer2模型的設置更加合理。接下來換一個上升運動輕緩的運動軌跡,濾波器參數設置不變,如圖4,5,6所示。

表2 IMM與各類算法濾波誤差均值對比表

圖2 IMM與各類算法濾波效果圖Fig.2 Filtering effect of IMM and various algorithms

圖3 目標劇烈俯仰運動的IMM與各類算法X軸定位誤差圖Fig.3 IMM of vigorous pitch motion of target and X-axis positioning error plot of various algorithms

圖4 目標輕微俯仰運動的IMM與各類算法X軸定位誤差圖Fig.4 IMM of slight pitch motion of target and X-axis positioning error plot of various algorithms

圖5 IMM與各類算法濾波效果圖Fig.5 Filtering effect of IMM and various algorithms

圖6 IMM三軸跟蹤誤差圖Fig.6 IMM triaxial tracking error
可以發現,在飛行目標針對基站俯仰角變化不劇烈時,IMM2模型收斂速度明顯提升且定位誤差在50 m左右。
接下來使用第一個目標軌跡研究IMM濾波參數和噪聲變化對定位影響,如表3所示。

表3 IMM不同參數下的目標跟蹤誤差
其中
μ1=[0.2,0.4,0.3,0.1],μ2=[0.4,0.2,0.1,0.3]

由表3可知,模型狀態方程的初始狀態轉移矩陣對于跟蹤效果有直接影響,但是隨著迭代更新,算法會根據各個子濾波器殘差和似然檢驗比實時選擇濾波器,所以算法依然可以有效跟蹤機動目標,且誤差波動不大;模型誤差相較于量測信息誤差,對于目標的定位精度影響更大。
本文在高斯白噪聲條件下,根據主被動站提供的與目標之間的距離和方位角,俯仰角量測信息進行對機動目標的定位跟蹤研究。鑒于使用單一運動模型的擴展卡爾曼濾波在目標機動時刻跟蹤誤差變大的問題,以TOA-AOA解析解為初始值,設計了一種融合勻速直線運動、勻加速直線運動、快機動Singer2運動和慢機動Singer1運動的IMM算法。通過仿真表明:
1)融合了CV、CA、慢機動Singer1和快機動Singer2的IMM濾波算法相較于單一運動模型的EKF算法,在目標與基站幾十千米范圍、時間量測誤差20 ns、角度誤差2°的條件下,對機動目標的跟蹤收斂誤差提升30%。
2)Singer模型對于機動運動的狀態轉移預測性要好于勻加速直線運動模型20%。
3)IMM濾波算法具有很強的對象特性。如果知道目標先驗加速度和機動參數范圍可以有效跟蹤機動目標;如果目標先驗參數信息不足,Singer模型的針對性會降低,需要根據量測信息估計目標加速度等信息的范圍。
4)本實驗針對高斯白噪聲,缺乏非高斯噪聲下的優化設計工作,擬根據粒子濾波原理,在非高斯白噪聲條件下,設計適用于任何噪聲的交互式多模型算法。