張敬艷 董 凱，2 孫 順

（1.海軍航空大學(xué)信息融合研究所 煙臺 264000）（2.中國電子科學(xué)研究院 北京 100041）

1 引言

艦載防空火控多普勒雷達(dá)系統(tǒng)的任務(wù)除了跟蹤飛機(jī)、艦船等傳統(tǒng)威脅之外，還要應(yīng)對高機(jī)動小目標(biāo)（巡航導(dǎo)彈等）和低慢小目標(biāo)（無人機(jī)等）。對空中來襲目標(biāo)進(jìn)行跟蹤時，一般在直角坐標(biāo)系下對目標(biāo)的狀態(tài)信息進(jìn)行描述，在極坐標(biāo)系下對雷達(dá)的觀測信息進(jìn)行描述，因此，目標(biāo)的狀態(tài)信息與雷達(dá)的觀測信息之間呈現(xiàn)非線性函數(shù)關(guān)系［1］，需要使用非線性濾波器。目前常見的非線性濾波估計(jì)算法有擴(kuò)展卡爾曼濾波（Extended Kalman Filter，EKF）、不敏卡爾曼濾波（Unscented Kalman Filter，UKF）和容積卡爾曼濾波（Cubature Kalman Filter，CKF）。文獻(xiàn)［2］中提到利用EKF算法對非線性系統(tǒng)存在的濾波問題進(jìn)行處理，該算法是利用卡爾曼濾波算法框架，通過泰勒級數(shù)展開將非線性函數(shù)線性化，但是該算法存在一些問題，首先會產(chǎn)生二次項(xiàng)以上的截?cái)嗾`差，影響系統(tǒng)的定位精度，其次算法的迭代更新需要計(jì)算雅可比矩陣，一階矩計(jì)算量小，可用于弱非線性系統(tǒng)，而二階及以上矩陣計(jì)算量較大，影響跟蹤的實(shí)時性，不適用于強(qiáng)非線性系統(tǒng)，且算法的初始狀態(tài)不易確定，初始狀態(tài)選擇不當(dāng)會導(dǎo)致濾波器發(fā)散，系統(tǒng)魯棒性差［3］。文獻(xiàn)［4］針對EKF算法存在的問題，提出利用UKF算法解決非線性系統(tǒng)的目標(biāo)跟蹤問題；UKF算法是根據(jù)目標(biāo)狀態(tài)的均值和協(xié)方差，利用UT變換確定采樣點(diǎn)，通過非線性函數(shù)傳遞對狀態(tài)來近似后驗(yàn)概率密度分布，并得到目標(biāo)狀態(tài)估計(jì)的均值和協(xié)方差，算法復(fù)雜度降低的同時，目標(biāo)狀態(tài)估計(jì)更為準(zhǔn)確。此外，由于不需要計(jì)算雅克比矩陣，故可對不可導(dǎo)的非線性函數(shù)進(jìn)行處理。但是當(dāng)階數(shù)大于3時，會損失部分采樣點(diǎn)對非線性函數(shù)后驗(yàn)分布特性的統(tǒng)計(jì)，導(dǎo)致系統(tǒng)精度降低。加拿大學(xué)者Arasaratnam和Haykin在文獻(xiàn)［5］中首次提出CKF算法，該算法基于三階球面徑向容積規(guī)則，利用2n個容積點(diǎn)對概率密度分布函數(shù)進(jìn)行統(tǒng)計(jì)逼近［6］，并且CKF算法的容積點(diǎn)和權(quán)值僅由狀態(tài)的維數(shù)唯一確定，可以提前計(jì)算和存儲，從而降低運(yùn)算量［7］，解決了UKF算法在強(qiáng)非線性和高階濾波時估計(jì)精度和濾波穩(wěn)定性低的問題。

針對艦載防空火控多普勒雷達(dá)目標(biāo)跟蹤過程中，量測信息的強(qiáng)非線性對目標(biāo)跟蹤精度帶來的影響。本文在容積卡爾曼濾波算法中引入多普勒雷達(dá)中的量測信息，在量測方程中聯(lián)合利用距離，角度和徑向速度，通過非線性容積點(diǎn)采樣，降低非線性量測對目標(biāo)跟蹤精度的影響。仿真結(jié)果表明，所提方法對目標(biāo)的位置與速度具有更高的估計(jì)精度。在此基礎(chǔ)上進(jìn)一步分析了多普勒量測誤差以及距離量測與多普勒量測的相關(guān)系數(shù)對目標(biāo)跟蹤性能的影響，為多普勒信息的引入提出參考條件。

2 問題描述

在兩維笛卡爾坐標(biāo)系中，對運(yùn)動目標(biāo)的狀態(tài)進(jìn)行描述，目標(biāo)的運(yùn)動模型一般可表示為［8］

式中，目標(biāo)的狀態(tài)向量X(k)=[x(k)，y(k)，(k)，(k)]Τ，x(k)和y(k)分別為目標(biāo)在x和y兩個方向上的位置分量，(k)和(k)分別為目標(biāo)在x和y兩個方上的速度分量，F(xiàn)k∈Rm×n為狀態(tài)轉(zhuǎn)移矩陣，G(k)為適當(dāng)維數(shù)的系數(shù)矩陣，u(k)為確定性輸入向量，V(k)是均值為零且方差為Q(k)的Gauss白噪聲序列。

假設(shè)一部在極坐標(biāo)系下位于坐標(biāo)原點(diǎn)的兩坐標(biāo)雷達(dá)，其包括目標(biāo)距離、方位角和徑向速度信息的量測方程可以表示為

式中，ρ(k)、θ(k)和(k)分別為距離、方位角和徑向速度的量測值，x(k)、y(k)、(k)和(k)為目標(biāo)真實(shí)位置和速度，為距離、方位角和徑向速度的真實(shí)值，為相應(yīng)的量測噪聲，假定其均為零均值高斯白噪聲序列，方差分別為σ2ρ、σ2θ和σ2，且不相關(guān)，不相關(guān)，但是是相關(guān)的，且其相關(guān)系數(shù)設(shè)為η。即有

則k時刻量測噪聲協(xié)方差在極坐標(biāo)系下為

3 DCKF濾波算法

3.1 容積點(diǎn)選取準(zhǔn)則

三階球面-徑向容積規(guī)則是根據(jù)目標(biāo)狀態(tài)的先驗(yàn)均值和協(xié)方差，通過容積規(guī)則選取容積點(diǎn)，利用非線性函數(shù)對容積點(diǎn)進(jìn)行傳遞，最后對傳遞后的容積點(diǎn)加權(quán)處理，得到目標(biāo)狀態(tài)的近似后驗(yàn)均值和協(xié)方差。

1）容積準(zhǔn)則

（1）笛卡爾坐標(biāo)中，積分形式如下［9］：

令x=ry，其中r為球體半徑r∈[0，+∞)，y為方向向量且yΤy=1，構(gòu)成Un={y∈Rn|yΤy=1} 的單位球體表面。易知xΤx=yΤrry=r2，滿足下式：

式中σ(y)表示球面方向矢量y所對應(yīng)的積分區(qū)域。

（2）球面-徑向坐標(biāo)積分

將積分分解為球面積分S(r)和徑向積分R(r)，得到

根據(jù)上式可以看出，容積準(zhǔn)則［10］是將原始積分轉(zhuǎn)化為某個n維幾何體的容積的積分式。

2）球面-徑向容積準(zhǔn)則

通過數(shù)值積分將球面積分與徑向積分近似為

式中，{yi，ωs，i}為計(jì)算面積積分的積分點(diǎn)集合與相對應(yīng)權(quán)重，Ms為相應(yīng)積分點(diǎn)數(shù)，{rj，ωr，j} 為計(jì)算徑向積分的積分點(diǎn)集合與相對應(yīng)權(quán)重，Mr為相應(yīng)積分點(diǎn)數(shù)。

根據(jù)上式，在球面-徑向容積準(zhǔn)則［11］下原積分的近似計(jì)算為

對于三階球面-徑向容積準(zhǔn)則，Mr=1，Ms=2n，容積點(diǎn)數(shù)是狀態(tài)向量的維數(shù)的兩倍。用三階球面-徑向容積準(zhǔn)則表示的標(biāo)準(zhǔn)高斯加權(quán)積分為

3.2 基于多普勒量測的CKF濾波實(shí)現(xiàn)

算法實(shí)現(xiàn)步驟如下。

1）時間更新

已知k時刻的狀態(tài)X(k)誤差協(xié)方差為P(k)，對P(k)進(jìn)行Cholesky分解

選擇狀態(tài)量容積點(diǎn)為

狀態(tài)容積點(diǎn)一步預(yù)測：

狀態(tài)一步預(yù)測：

狀態(tài)估計(jì)協(xié)方差一步預(yù)測：

2）量測更新

對P(k+1|k)進(jìn)行Cholesky分解：

計(jì)算狀態(tài)預(yù)測容積點(diǎn)：

量測容積點(diǎn)的一步預(yù)測：

量測的一步預(yù)測：

新息協(xié)方差：

3）濾波增益

相關(guān)協(xié)方差一步預(yù)測：

濾波增益：

4）方程更新

狀態(tài)更新：

協(xié)方差更新：

4 仿真分析

考慮目標(biāo)為在二維空間做勻速直線運(yùn)動物體，假定目標(biāo)的初始位置為（40，40）km，目標(biāo)的初始速度為（-30，-30）m/s，雷達(dá)的掃描周期為T=5s，雷達(dá)的測距誤差σρ=100m，測角誤差為σθ=5mrad，步長為200。將引入多普勒量測的容積卡爾曼濾波（DCKF）與引入多普勒量測的擴(kuò)展卡爾曼濾波（DEKF）、引入多普勒量測的不敏卡爾曼濾波（DUKF）、未引入多普勒量測的擴(kuò)展卡爾曼濾波（EKF）、未引入多普勒量測的不敏卡爾曼濾波（UKF）和未引入多普勒量測的容積卡爾曼濾波（CKF）對比，進(jìn)行1000次蒙特卡洛實(shí)驗(yàn)，考察濾波器對目標(biāo)位置和速度的均方根誤差（Root Meaning Square Error，RMSE）［12］：

4.1 兩種場景下的算法仿真分析

表1 參數(shù)設(shè)置

圖1 場景1條件下不同算法誤差

圖2 場景2條件下不同算法誤差

進(jìn)一步綜合圖1和圖2的仿真結(jié)果可知，當(dāng)η相同時，算法的濾波效果會受到的影響。初步仿真結(jié)果表明，較小時，所提算法濾波效果改善明顯，反之，目標(biāo)跟蹤效果改善不明顯。

4.2 量測誤差與相關(guān)系數(shù)對算法性能影響分析

為進(jìn)一步分析η相同、不同時的濾波性能，令η=0.1，取 0.1m/s～12m/s，其他仿真條件不變，得到仿真結(jié)果如圖3（a）所示；令η=0.9，取0.1m/s～12m/s，其它仿真條件不變，得到仿真結(jié)果如圖3（b）所示，其中算法性能使用目標(biāo)跟蹤過程中最后100個時刻的位置RMSE的均值來進(jìn)行評估。

圖3 不同條件下的算法性能對比

由圖3（a）可知，當(dāng)η=0.1，為0.1m/s～0.3m/s，由文獻(xiàn)［7］可知，由于狀態(tài)的初始誤差協(xié)方差較大，量測噪聲較小，各算法的狀態(tài)估計(jì)偏差較大，性能較差，無法對系統(tǒng)進(jìn)行有效的狀態(tài)濾波；隨著的增大，未引入多普勒信息的濾波算法性能穩(wěn)定不變，引入多普勒信息的濾波方法的目標(biāo)位置誤差逐漸增大，并趨近于未引入多普勒信息的濾波算法，其中本文所提出的DCKF算法對位置的估計(jì)優(yōu)于其他濾波算法。

由圖3（b）可知，當(dāng)η=0.9時，隨著的增大，未引入多普勒信息的濾波算法性能穩(wěn)定不變，引入多普勒信息的濾波算法對位置的估計(jì)性能均有所下降，且在，濾波效果要差于未引入多普勒信息的濾波算法，但本文所提出的DCKF算法對位置的估計(jì)仍優(yōu)于引入多普勒信息的其他濾波算法。

圖4 不同η條件下的算法性能對比

綜上所述，未引入多普勒信息的濾波算法與η和無關(guān)，而引入多普勒信息的濾波算法受η與的雙重影響，在本文給定的仿真條件下，可以得到如下結(jié)論。

1）當(dāng)η較小，隨著增大，引入多普勒信息的濾波算法趨近于未引入多普勒信息的濾波算法，其中所提的DCKF算法對位置的估計(jì)精度明顯優(yōu)于其他濾波算法；

2）當(dāng)η較大，較小時，所提出的DCKF算法對位置的估計(jì)優(yōu)于其他濾波算法；

3）當(dāng)η較大，較大時，所提出的DCKF算法對位置的估計(jì)優(yōu)于引入多普勒信息的其他濾波算法。此時，未引入多普勒信息的濾波算法的位置估計(jì)效果相對較好，說明在這種情況下濾波算法中不宜引入多普勒信息。

5 結(jié)語

本文針對基于多普勒雷達(dá)徑向速度信息的目標(biāo)跟蹤問題，提出了一種基于容積卡爾曼濾波（Cu?bature Kalman Filter，CKF）的多普勒雷達(dá)目標(biāo)跟蹤算法。在量測方程中聯(lián)合利用距離、角度和徑向速度信息，并采用CKF算法降低非線性量測方程對目標(biāo)跟蹤精度的影響。通過對不同算法的仿真結(jié)果對比分析驗(yàn)證了所提方法的有效性和優(yōu)越性。其中，當(dāng)徑向速度誤差精度較高時，所提DCKF方法的性能對相關(guān)系數(shù)不敏感；反之，當(dāng)徑向速度誤差精度較低時，選取較小的相關(guān)系數(shù)，可使所提算法性能最佳；當(dāng)徑向速度誤差較大，選取較大的相關(guān)系數(shù)時，不宜引入多普勒信息。