劉 代, 趙永波, 高 劍, 李 偉
(1. 西安電子科技大學雷達信號處理國家重點實驗室, 陜西 西安 710071;2. 西安電子工程研究所, 陜西 西安 710100)
在雷達數據處理中,將概率數據關聯法(probabilistic data association, PDA)算法和Kalman濾波算法相結合使用越來越廣泛[1-5]。PDA算法能在跟蹤波門內對多個過門限的量測數據進行數據關聯,確定各個量測來源于目標的概率,并利用其對新息加權獲得目標的狀態估計。在雷達目標跟蹤中,Kalman濾波算法是最小方差意義下的最優算法[6]。將PDA算法與Kalman算法結合使用很好地得到目標的狀態估計并使目標的最小方差最優。文獻[1-5]將PDA算法和Kalman濾波算法相結合應用到MUSIC信號處理、碼分多址(code division multiple access, CDMA)通信和多輸入多輸出(multiple-input multiple-output, MIMO)信道,得到很好的濾波效果。文獻[7]利用徑向速度觀測誤差與距離觀測誤差相關條件下,提出序貫處理的濾波算法。文獻[7-8]只是在Kalman濾波中引入了徑向速度,卻沒有利用更加準確的動目標檢測(moving targets detection, MTD)測出的徑向速度,未與PDA算法結合使用。目前解MTD速度模糊方法較多[9-13],本文基于MTD測出的不模糊徑向速度,提出了一種利用徑向速度的Kalman-PDA算法(定義為KalmanV-PDA),利用此徑向速度建立速度波門,推導了速度波門表達式,在量測方程中引入徑向速度維,利用MTD測出的徑向速度實時更新目標觀測值中的徑向速度,而不是利用目標三軸速度和位置的數學表達式來計算量測徑向速度。利用徑向速度的Kalman-PDA算法相比標準的Kalman-PDA濾波算法目標跟蹤性能得到很大提升。
假設目標在XYZ平面內做勻速直線運動。其狀態方程為
xk=Fk|k-1xk-1+wk-1
(1)
目標量測模型
zk=h(xk)+vk
(2)
式中,zk為k時刻的觀測值,zk=[xTk,yTk,zTk,vrk]T,xTk=rTkcosθTkcosφTk,yTk=rTksinθTkcosφTk,zTk=rTksinφTk,rTk為k時刻雷達測得的目標徑向距離,θTk為k時刻雷達測得的目標方位角,φTk為k時刻雷達測得的目標俯仰角,在目標觀測模型中增加徑向速度維,構成了四維觀測向量,vrk代表目標徑向速度;h(·)表示觀測函數;vk為k時刻的觀測噪聲,其協方差矩陣為R。
步驟1預測k時刻目標的狀態
(3)
步驟2計算第i個測量值對應的新息
(4)
計算真實新息協方差矩陣:
(5)
式中,Hk=(ΔxkhT(xk))T為雅克比矩陣;Rk為k時刻的測量誤差協方差矩陣。
計算狀態預測協方差矩陣:
pk|k-1=FPk-1|k-1FT+Qk-1
(6)
步驟3計算第i個量測值源于目標的條件概率βi[15]
(7)
步驟4根據計算的條件概率βi,將各個觀測的新息進行組合
(8)
即可獲取目標狀態的濾波結果:
(9)
濾波器的增益矩陣:
(10)
步驟5計算濾波后的狀態協方差矩陣
Pk|k=[I-KkH]Pk|k-1
(11)
上述算法是將標準Kalman濾波算法與PDA算法結合起來使用,PDA算法做數據關聯, Kalman濾波做均方誤差最小收斂。
脈沖多普勒(pulse Doppler, PD)雷達是一種常用雷達體制,利用多普勒效應檢測目標信息,利用相參處理的方式提取目標徑向速度。采用窄帶多普勒濾波器組覆蓋整個多普勒頻率范圍,如圖1所示。窄帶多普勒濾波器組可以通過快速傅里葉變換完成。
圖1 多普勒濾波器組特性Fig.1 Characteristics of Doppler filter bank
假設多普勒通道個數N,目標的回波位于M號濾波器,脈沖重復周期TR,相參處理周期為CPI,脈沖重復頻率fr, 脈沖重復周期為PRT,雷達載頻fc,雷達波長λ,電磁波傳播的速度c,可以得出:
目標的多普勒頻移為
(12)
目標的徑向速度為
(13)
MTD測速最大不模糊速度為
(14)
MTD測速精度為
(15)
通過上面公式可以看出只要選擇合適的雷達載頻和和脈沖相參處理周期就能保證MTD測速具有很好的測速精度,通常雷達設計時MTD測速精度都較高。
目標k時刻徑向速度與k-1時刻徑向速度差值為
(16)
假設k-1時刻目標速度為vk-1,則
(17)
代入式(16)可得
(18)
則徑向速度差值波門為
(19)
式中,vmax為目標最大運動速度。
Zk|k-1=H(k)Xk|k-1+Vk
(20)
其中
在KalmanV-PDA算法中用MTD測出的徑向速度和目標位置信息組成目標觀測值zk=[xTk,yTk,zTk,vrk]T,其中(xTk,yTk,zTk)為k時刻目標在直角坐標系中的坐標,Vrk為k時刻MTD測出的徑向速度,zk是由1,2,…,lk個觀測值組成。
利用KalmanV-PDA算法步驟如下:
步驟1預測k時刻目標的狀態
(21)
計算第i個測量值對應的新息:
(22)
計算真實新息協方差矩陣:
(23)
計算狀態預測協方差矩陣:
pk|k-1=FPk-1|k-1FT+Qk-1
(24)
步驟3計算第i個量測值源于目標的條件概率βi[15]
(25)
步驟4根據條件概率βi,將各個觀測的新息進行組合
(26)
目標狀態的濾波結果:
(27)
濾波器的增益矩陣
(28)
步驟5計算濾波后的狀態協方差矩陣
Pk|k=[I-KkH]Pk|k-1
(29)
為了驗證文中所提KalmanV-PDA算法的有效性,本文進行了如下場景仿真,假設目標遠離雷達飛行,從初始位置為(12 000 m,12 000 m,1 000 m)飛出,X軸方向速度200 m/s,Y軸方向速度100 m/s作勻速直線運動,高度保持1 000 m不變。雷達相參處理周期(coherent process interval recurrence, CPI)為80 ms,脈沖重復頻率(pulse recurrence frequency,PRF)為1 kHz,雷達載頻fc為2 000 MHz,電磁波傳播的速度c為3×108m/s,每1 s觀測一次數據,觀測時間為100 s。則MTD測速精度為0.937 m/s,仿真中取測量距離、方位角誤差分別為50 m和8 mrad,徑向速度誤差分別取1 m/s和5 m/s。作200次蒙特卡羅仿真。目標狀態轉移矩陣為
將本文提出的KalmanV-PDA與標準Kalman-PDA濾波算法進行了對比,分別從目標位置的均方根誤差(root mean square error, RMSE)和速度RMSE以及收斂速度三方面進行了對比,從圖3目標距離RMSE曲線可以看出KalmanV-PDA算法性能最好,具體體現在收斂速度最快,位置RMSE最小。從圖4目標速度RMSE曲線可以看出KalmanV-PDA也是性能最好的,具體體現在目標速度的RMSE最小,收斂速度最快。
本文分析了MTD測速誤差對算法的影響,從圖3和圖4可以看出在測速誤差取為測速精度的5倍時,KalmanV-PDA算法比標準Kalman-PDA濾波算法性能優越,MTD測速誤差對濾波算法性能影響較小,可以忽略不計。分析KalmanV-PDA算法性能好的原因就是由于采用MTD檢測出的目標精確徑向速度,更新了觀測方程中的目標觀測值,提高了目標跟蹤性能。
圖2 目標運動軌跡Fig.2 Motion trajectory of target
圖3 目標位置RMSE曲線Fig.3 RMSE curve of target
圖4 目標速度RMSE曲線Fig.4 RMSE curve of target velocity
利用MTD測速的KalmanV-PDA算法,相比標準的Kalman-PDA算法,收斂速度快,目標位置RMSE和速度RMSE小,提高了目標預測精度,改善了目標跟蹤性能,而且通過分析 MTD測速誤差對跟蹤性能影響較小。利用徑向速度波門可以應用到航跡相關中,這可以作為后續研究內容。