馬 娟，范玉珠，李 樂

(1. 中國電子科技集團(tuán)公司第三十八研究所，安徽合肥 230088； 2. 中國人民解放軍63629部隊(duì)，北京 100162)

0 引言

PD雷達(dá)常用于機(jī)載預(yù)警和火控雷達(dá)領(lǐng)域，利用目標(biāo)與雷達(dá)之間的相對運(yùn)動(dòng)產(chǎn)生的多普勒效應(yīng)進(jìn)行目標(biāo)信息提取和處理。機(jī)載PD雷達(dá)具有良好的下視和上視探測能力、精確的多普勒速度量測能力和同時(shí)探測數(shù)十批甚至上千批目標(biāo)的能力。雷達(dá)數(shù)據(jù)處理系統(tǒng)對量測數(shù)據(jù)進(jìn)行關(guān)聯(lián)和濾波，形成可靠的雷達(dá)情報(bào)[1-2]。航跡跟蹤中濾波算法的性能直接影響雷達(dá)情報(bào)質(zhì)量。

機(jī)載雷達(dá)在以載機(jī)為中心的運(yùn)動(dòng)平臺上獲得量測數(shù)據(jù)，量測數(shù)據(jù)直接描述的是雷達(dá)和目標(biāo)的相對位置變化[3]，一般情況下雷達(dá)和目標(biāo)的相對運(yùn)動(dòng)狀態(tài)方程是非線性方程，很難直接進(jìn)行濾波和預(yù)測。機(jī)載PD雷達(dá)濾波時(shí)一般將量測數(shù)據(jù)轉(zhuǎn)換到以地面測量站為中心的慣性坐標(biāo)系(簡稱地面慣性坐標(biāo)系)下，再采用與地面雷達(dá)相同的處理方法進(jìn)行濾波。該方法存在以下缺點(diǎn)：目標(biāo)在某一時(shí)刻只能有一種假設(shè)的運(yùn)動(dòng)模型，如果目標(biāo)運(yùn)動(dòng)狀態(tài)與運(yùn)動(dòng)模型不匹配，會(huì)導(dǎo)致濾波器嚴(yán)重發(fā)散；量測數(shù)據(jù)轉(zhuǎn)換到地面慣性坐標(biāo)系下進(jìn)行濾波，用經(jīng)過非線性變換后的數(shù)據(jù)濾波，數(shù)據(jù)的噪聲特性失真，導(dǎo)致濾波器精度下降；作為PD雷達(dá)特色的目標(biāo)多普勒速度在濾波過程中沒有得到利用，損失了一維重要的信息。

目前機(jī)動(dòng)目標(biāo)跟蹤的主流算法是交互多模型(Interactive Multi-Model，IMM)濾波算法[4-6]。本文方法采用IMM濾波器，模型集由勻速直線運(yùn)動(dòng)(CV)模型、當(dāng)前統(tǒng)計(jì)(CS)模型和估計(jì)角速度的勻速轉(zhuǎn)彎運(yùn)動(dòng)(CT)模型組成。其中，CV模型和CS模型采用線性卡爾曼濾波實(shí)現(xiàn)，CT模型采用容積卡爾曼濾波實(shí)現(xiàn)[7]。濾波器量測誤差協(xié)方差矩陣由天線坐標(biāo)系下的量測誤差和天線坐標(biāo)系下的量測值進(jìn)行無偏轉(zhuǎn)換近似計(jì)算。利用機(jī)載PD雷達(dá)多普勒速度量測精度高的特性，將其應(yīng)用到模型概率初始估計(jì)和模型參數(shù)初始估計(jì)中，提高模型初始參數(shù)的精度。

1 基于CV模型、CS模型和CT模型的交互多模型濾波器

基于CV模型、CS模型和CT模型的交互多模型濾波器周期運(yùn)行，由數(shù)據(jù)驅(qū)動(dòng)，每一次關(guān)聯(lián)成功后，濾波器的工作步驟如下：

1) 量測數(shù)據(jù)預(yù)處理

接收到每一個(gè)周期的探測數(shù)據(jù)后進(jìn)行量測數(shù)據(jù)的預(yù)處理，包括量測誤差協(xié)方差矩陣計(jì)算和濾波器輸入量測數(shù)據(jù)計(jì)算。

2) 模型初始化

濾波器第一次調(diào)用時(shí)，初始化模型概率、模型參數(shù)和Markova轉(zhuǎn)移矩陣。

3) 狀態(tài)向量和協(xié)方差初始化

濾波器第一次調(diào)用時(shí)，根據(jù)先驗(yàn)信息和已獲得的量測數(shù)據(jù)初始化各模型的狀態(tài)向量和協(xié)方差 (j=1,2,3分別對應(yīng)CS、CV、CT模型，下同)。

4) 相互作用

計(jì)算與每一個(gè)模型匹配的濾波器的混合初始條件。模型j在k時(shí)刻的初始狀態(tài)向量為

(1)

相應(yīng)的協(xié)方差為

(2)

5) 卡爾曼濾波

6) 計(jì)算模型的可能性

(3)

7) 模型概率更新

(4)

(5)

8) 組合估計(jì)

狀態(tài)向量和協(xié)方差的最終組合結(jié)果為

(6)

(7)

2 量測數(shù)據(jù)預(yù)處理

濾波器啟動(dòng)前需要先進(jìn)行量測數(shù)據(jù)的預(yù)處理，包括量測誤差協(xié)方差矩陣計(jì)算和濾波器輸入量測數(shù)據(jù)計(jì)算。計(jì)算量測誤差協(xié)方差矩陣時(shí)，先由直接測量得到的目標(biāo)在天線坐標(biāo)系下的距離、方位和載機(jī)姿態(tài)角信息獲取對應(yīng)的距離量測誤差、方位量測誤差和載機(jī)慣性坐標(biāo)系下的極坐標(biāo)量測值，再采用無偏轉(zhuǎn)換方法計(jì)算量測誤差協(xié)方差矩陣。載機(jī)慣性極坐標(biāo)系下的量測值和對應(yīng)時(shí)刻的載機(jī)大地坐標(biāo)位置進(jìn)行坐標(biāo)變換獲得地面慣性坐標(biāo)系下的位置作為濾波器輸入量測數(shù)據(jù)[8]。

由式(8)計(jì)算目標(biāo)在天線坐標(biāo)系下的直角坐標(biāo)位置P1=[x1,y1,z1]T。其中r1，a1，e1為根據(jù)第k個(gè)周期雷達(dá)探測的目標(biāo)原始回波解算出的位置坐標(biāo)，分別表示目標(biāo)距離(單位：m)、目標(biāo)方位角(以機(jī)頭為中心，單位：rad)和目標(biāo)俯仰角(單位：rad)。

(8)

根據(jù)探測點(diǎn)跡距離r1和方位角a1，由雷達(dá)特性和統(tǒng)計(jì)分布等先驗(yàn)信息，獲得距離量測誤差dr和方位量測誤差da。計(jì)算量測誤差協(xié)方差矩陣為

(9)

其中，

(cosh(2·da2)-cosh(da2))+

(sin(a1))2·(sinh(2·da2)-sinh(da2)))+

dr2·exp(-2·da2)·((cos(a1))2·

(2·cosh(2·da2)-cosh(da2))+

(sin(a1))2·(2·sinh(2·da2)-sinh(da2)))

(cosh(2·da2)-cosh(da2))+

(cos(a1))2·(sinh(2·da2)-sinh(da2)))+

dr2·exp(-2·da2)·((sin(a1))2·

(2·cosh(2·da2)-cosh(da2))+

(cos(a1))2·(2·sinh(2·da2)-sinh(da2)))

r12=sin(a1)·cos(a1)·exp(-4·da2)·

r21=r12

由式(10)計(jì)算旋轉(zhuǎn)變換公式T1，其中，α、β、γ分別為載機(jī)的偏航角、俯仰角和橫滾角(單位：rad)。

T1=

(10)

由式(11)計(jì)算目標(biāo)在載機(jī)慣性坐標(biāo)系下的直角坐標(biāo)值P2=[x2,y2,z2]T。

P2=T1*P1

(11)

再由式(12)計(jì)算目標(biāo)在載機(jī)慣性坐標(biāo)系下的極坐標(biāo)值。

(12)

已知P2、載機(jī)位置和地面測量站的位置，利用大地坐標(biāo)變換公式計(jì)算目標(biāo)在測站坐標(biāo)系下的直角坐標(biāo)位置Zk=[x,y,z]T，Zk即是第k個(gè)周期濾波器的輸入量測值。

3 濾波器實(shí)現(xiàn)

3.1 模型的初始概率和參數(shù)初始化

交互多模型濾波器中各模型的初始概率和模型參數(shù)根據(jù)速度、加速度和徑向速度估計(jì)，加速度和徑向速度的變化對應(yīng)目標(biāo)的機(jī)動(dòng)情況，對一批目標(biāo)第一次濾波時(shí)，根據(jù)徑向速度、加速度和目標(biāo)機(jī)動(dòng)的先驗(yàn)分布，獲得模型初始概率和初始機(jī)動(dòng)參數(shù)。

1) 模型概率初始化

(13)

2) 初始化模型預(yù)置參數(shù)

由dλ和CS模型機(jī)動(dòng)因子的先驗(yàn)分布獲得機(jī)動(dòng)因子α的最大后驗(yàn)估計(jì)。CT模型的轉(zhuǎn)彎角速度w初始化為0。

α=argmaxp(α|dλ)

(14)

3) 初始化Markova轉(zhuǎn)移矩陣

濾波器第一次調(diào)用時(shí)，初始化模型組的Mar-kova轉(zhuǎn)移矩陣P，其中pij(i=1,2,3;j=1,2,3)表示模型i到模型j的轉(zhuǎn)移概率。轉(zhuǎn)移矩陣中CV模型與CS模型可相互轉(zhuǎn)換，CV模型與CT模型可相互轉(zhuǎn)換，CS模型與CT模型之間不能相互轉(zhuǎn)換。

(15)

3.2 容積卡爾曼濾波器

時(shí)間更新的計(jì)算步驟如下：

(16)

(17)

(18)

4) 通過經(jīng)CT模型狀態(tài)方程傳遞后的2n個(gè)cubature點(diǎn)的狀態(tài)來估計(jì)k時(shí)刻的狀態(tài)

(19)

5) 計(jì)算k時(shí)刻2n個(gè)cubature點(diǎn)的狀態(tài)對應(yīng)的協(xié)方差的一步預(yù)測

(20)

式中，Qk-1為過程噪聲協(xié)方差矩陣。

觀測更新的計(jì)算步驟如下：

1) 對協(xié)方差的一步預(yù)測Pk,k-1進(jìn)行喬列斯基分解：

Pk,k-1=Lk,k-1*(Lk,k-1)T

(21)

(22)

(23)

(24)

3) 求k時(shí)刻的新息協(xié)方差：

(25)

求k時(shí)刻狀態(tài)預(yù)測誤差與量測預(yù)測誤差的互協(xié)方差矩陣：

(26)

計(jì)算容積卡爾曼濾波的濾波器增益:

(27)

計(jì)算k時(shí)刻的新息：

(28)

進(jìn)行k時(shí)刻的狀態(tài)更新：

(29)

進(jìn)行k時(shí)刻的協(xié)方差更新：

(30)

4 仿真實(shí)驗(yàn)

為了驗(yàn)證算法的有效性，結(jié)合目標(biāo)仿真系統(tǒng)進(jìn)行了仿真實(shí)驗(yàn)，實(shí)驗(yàn)結(jié)果顯示本文方法明顯提高了濾波精度，對機(jī)動(dòng)和非機(jī)動(dòng)狀態(tài)的目標(biāo)均有穩(wěn)定的濾波效果。

4.1 實(shí)驗(yàn)一

實(shí)驗(yàn)一采用MATLAB仿真，仿真參數(shù)如下：

載機(jī)運(yùn)動(dòng)軌跡：載機(jī)與目標(biāo)初始距離100 km，載機(jī)以100 m/s的速度向正北方向進(jìn)行勻速直線運(yùn)動(dòng)。

目標(biāo)運(yùn)動(dòng)軌跡：目標(biāo)以200 m/s的速度進(jìn)行40 s的勻速直線運(yùn)動(dòng)，之后以0.1 rad/s的角速度進(jìn)行60 s的勻速轉(zhuǎn)彎運(yùn)動(dòng)，之后勻速運(yùn)動(dòng)30 s，再以-0.1 rad/s的角速度進(jìn)行60 s勻速轉(zhuǎn)彎運(yùn)動(dòng)，最后勻速運(yùn)動(dòng)30 s。

濾波器參數(shù)設(shè)置：距離量測誤差150 m；方位角量測誤差0.1°；幀周期10 s；蒙特卡洛仿真次數(shù)100。

圖1所示為目標(biāo)真實(shí)運(yùn)動(dòng)軌跡的采樣值，圖2為觀測值和濾波值的絕對距離誤差的均方根。由圖可見，在目標(biāo)機(jī)動(dòng)或非機(jī)動(dòng)狀態(tài)下，濾波結(jié)果都優(yōu)于觀測結(jié)果，濾波器具有良好的自適應(yīng)性能。

圖1 實(shí)驗(yàn)一目標(biāo)軌跡

圖2 實(shí)驗(yàn)一目標(biāo)濾波前后誤差

4.2 實(shí)驗(yàn)二

實(shí)驗(yàn)數(shù)據(jù)來源于某機(jī)載雷達(dá)仿真系統(tǒng)，目標(biāo)進(jìn)行隨機(jī)機(jī)動(dòng)，目標(biāo)軌跡如圖3所示，濾波前后目標(biāo)位置的絕對距離均方根誤差如圖4所示。濾波后精度明顯優(yōu)于測量精度。

圖3 實(shí)驗(yàn)二目標(biāo)軌跡

圖4 實(shí)驗(yàn)二目標(biāo)濾波前后誤差

5 結(jié)束語

本方法避免了模型與濾波器不匹配導(dǎo)致的濾波器精度嚴(yán)重下降，減少了平臺運(yùn)動(dòng)和坐標(biāo)變換對量測誤差的影響。將多普勒速度應(yīng)用到模型概率初始化和模型參數(shù)初始化中，進(jìn)一步提高了濾波器的初始精度和可靠性。線性濾波采用線性卡爾曼濾波，非線性濾波采用容積卡爾曼濾波，既保證了濾波器的精度和數(shù)值穩(wěn)定性，又使運(yùn)算量最小。仿真實(shí)驗(yàn)表明該算法提高了機(jī)載PD雷達(dá)的航跡情報(bào)質(zhì)量，尤其是提高了航跡濾波的精度和穩(wěn)定性。