丁海龍，趙溫波

(解放軍陸軍軍官學(xué)院， 合肥 230031)

?

·信號/數(shù)據(jù)處理·

雷達(dá)組網(wǎng)GMPHDF關(guān)鍵參數(shù)研究

丁海龍，趙溫波

(解放軍陸軍軍官學(xué)院， 合肥 230031)

高斯混合概率假設(shè)密度濾波具有嚴(yán)密的數(shù)學(xué)基礎(chǔ)，適合跟蹤弱信噪比多目標(biāo)，但其目標(biāo)分布協(xié)方差P和高斯元素裁剪門限T至今未有合理取值規(guī)則，影響了跟蹤效果，且殘差協(xié)方差S參與增益計算時需要對其進(jìn)行逆計算，如果S為非正定，會導(dǎo)致計算發(fā)散。針對上述問題，通過概率統(tǒng)計方法推導(dǎo)了參數(shù)P和T的取值規(guī)則，通過Cholesky和QR分解，確定了參數(shù)S的計算規(guī)則。仿真比較分析表明：文中提出的目標(biāo)分布協(xié)方差P、裁剪門限T和殘差協(xié)方差S的計算規(guī)則用于雷達(dá)組網(wǎng)高斯混合概率假設(shè)密度濾波跟蹤弱信噪比多目標(biāo)時，能較高精度地跟蹤到所有目標(biāo)，且沒有帶來多余計算負(fù)擔(dān)。

雷達(dá)組網(wǎng)；高斯混合概率假設(shè)密度濾波；參數(shù)

0 引 言

現(xiàn)代戰(zhàn)爭條件下，反輻射導(dǎo)彈、隱身目標(biāo)、超低空突防和綜合性電子干擾等裝備廣泛使用，使得雷達(dá)系統(tǒng)獲得的目標(biāo)回波能量大為減少，對目標(biāo)的檢測跟蹤效果也大為降低，目標(biāo)信息急劇減少，雜波信息卻迅速增加，這種情況我們稱之為雷達(dá)弱信噪比(Weak Signal to Noise Ratio, WSNR)目標(biāo)檢測跟蹤[1]。雷達(dá)組網(wǎng)系統(tǒng)[2-3]為多源數(shù)據(jù)融合系統(tǒng)，一定程度上可以提高空域目標(biāo)探測的空間/時間分辨率，但WSNR條件下，由于大量測量噪聲和雜波信息的存在，弱信息目標(biāo)很容易淹沒在雜波虛警當(dāng)中，當(dāng)檢測跟蹤目標(biāo)數(shù)未知的多個目標(biāo)時，雷達(dá)檢測到的數(shù)據(jù)并非就是目標(biāo)信息，為了準(zhǔn)確跟蹤到目標(biāo)，需要對目標(biāo)數(shù)目和狀態(tài)同時進(jìn)行估計。傳統(tǒng)的方法是將雷達(dá)測量信息與估計目標(biāo)進(jìn)行數(shù)據(jù)關(guān)聯(lián)計算，數(shù)據(jù)關(guān)聯(lián)會造成大量的計算量，給工程應(yīng)用帶來不便。概率假設(shè)密度濾波(Probability Hypothesis Density Filter, PHDF)[4-5]以隨機(jī)有限集和貝葉斯濾波為其理論基礎(chǔ)，借鑒貝葉斯單目標(biāo)跟蹤辦法，將多目標(biāo)狀態(tài)和測量信息分別作為有限集來處理，避免了復(fù)雜的數(shù)據(jù)關(guān)聯(lián)。高斯混合概率假設(shè)密度濾波(Gaussian Mixture PHDF, GMPHDF)[6-8]是PHDF的一種解析解實(shí)現(xiàn)形式，實(shí)現(xiàn)了PHDF的工程實(shí)際應(yīng)用。然而，對于GMPHDF當(dāng)中多目標(biāo)方差P和裁剪門限T的選取，至今未提出合理的選取規(guī)則，僅憑經(jīng)驗(yàn)進(jìn)行選取，很容易增大跟蹤誤差，甚至造成跟蹤失敗，并且在計算濾波增益K時，需要求殘差協(xié)方差的逆矩陣 ，如果S是非正定矩陣，會造成計算發(fā)散。本文通過數(shù)學(xué)分析，推導(dǎo)了P、T合理的選取規(guī)則，通過Cholesky和〗QR分解，解決了計算增益時的發(fā)散問題，并通過成功的雷達(dá)組網(wǎng)弱信噪比環(huán)境跟蹤多目標(biāo)仿真實(shí)驗(yàn)，驗(yàn)證了P、T選取的合理性和S計算增益的有效性，得到了很好的效果。

1 高斯混合概率假設(shè)密度濾波

GMPHDF的理論基礎(chǔ)是隨機(jī)有限集和貝葉斯濾波理論，它是在高斯條件下通過傳播多目標(biāo)隨機(jī)有限集的一階矩(概率假設(shè)密度)來近似后驗(yàn)分布，不考慮一階矩以上的統(tǒng)計特性。其對應(yīng)的迭代等式為

Dk+1|k(X)= ∫ps,k+1(ζ)fk+1|k(X|ζ)Dk(ζ)dζ+

∫βk+1|k(X|ζ)Dk(ζ)dζ+γk+1(X)

(1)

Dk+1(X)= [1-pD,k+1(X)]Dk+1|k(X)+

(2)

式中：X為目標(biāo)狀態(tài)的隨機(jī)有限集；Z為測量值的隨機(jī)有限集；ps,k+1(ζ)為k時刻狀態(tài)為ζ的目標(biāo)在當(dāng)前k+1時刻存活的概率；fk+1|k(X|ζ)為狀態(tài)轉(zhuǎn)移概率密度；βk+1|k(X|ζ)為k時刻狀態(tài)為ζ的目標(biāo)在k+1時刻衍生的目標(biāo)概率假設(shè)密度；γk+1(X)為k+1時刻新生目標(biāo)概率假設(shè)密度；pD,k+1(X)為當(dāng)前k+1時刻檢測到目標(biāo)狀態(tài)X的概率；Lk+1(Z|X)為似然函數(shù)；κk+1(Z)為雜波隨機(jī)有限集密度。

根據(jù)參考文獻(xiàn)[6-8]，GMPHDF濾波流程為

(3)

(4)

(5)

(6)

(7)

Dk+1|k(x)=Dγ,k+1|k(x)+Dβ,k+1|k(x)+

DS,k+1|k(x)

(8)

Dk+1(x)= (1-pD)Dk+1|k(x)+

(9)

(10)

(11)

(12)

(13)

(14)

式(3)～式(5)分別為計算當(dāng)前K+1時刻新生、衍生、存在目標(biāo)的預(yù)測概率假設(shè)密度(Probability Hypothesis Density, PHD)，式(6)、式(7)分別為計算衍生、存在目標(biāo)的預(yù)測協(xié)方差矩陣，式(8)是新生、衍生、存在目標(biāo)預(yù)測PHD的高斯混合，式(9)為計算更新后所有目標(biāo)高斯元素的PHD，式(10)是被檢測到的更新后高斯元素，式(11)～式(14)分別計算高斯元素PHD更新時用到的權(quán)值、狀態(tài)值、協(xié)方差、增益。

濾波更新以后高斯元素會大量增加，導(dǎo)致過重的計算負(fù)擔(dān)，需要對更新后的高斯元素進(jìn)行裁剪與合并

(15)

(16)

式中：T為PHD裁剪門限，U為合并門限。將不滿足式(15)的高斯元素裁剪掉，不參與下一步計算，再將滿足式(15)的高斯元素根據(jù)式(16)進(jìn)行合并。

最后，提取權(quán)值大于0.5的元素為估計的目標(biāo)數(shù)與狀態(tài)。

2 高斯混合概率假設(shè)密度濾波參數(shù)研究

GMPHDF實(shí)現(xiàn)了PHDF的工程實(shí)踐應(yīng)用，但至今未給出GMPHDF裁剪處理當(dāng)中門限T和新生、衍生目標(biāo)分布協(xié)方差Pλ、Pβ具體取值規(guī)則，多是由經(jīng)驗(yàn)進(jìn)行取值。裁剪門限T取值過大會丟失有用信息，容易導(dǎo)致檢測跟蹤失敗，取值過小會引入過量高斯元素(包括大量雜波高斯元素)參與下一步高斯元素合并計算，增大計算量，同時影響合并計算精度。GMPHDF各高斯元素服從高斯分布N(x;m,P)，裁剪門限T用于與高斯元素PHD進(jìn)行比較，因此，T具體取值決定于高斯元素的PHD，各元素PHD受分布協(xié)方差P影響。從這個意義上說，同樣需要先確定各高斯元素分布協(xié)方差的取值。

2.1 確定目標(biāo)分布方差取值規(guī)則

新生、衍生、存在目標(biāo)的分布協(xié)方差Pλ、Pβ、PS描述高斯元素在均值附近的發(fā)散情況，以更準(zhǔn)確地表示真實(shí)目標(biāo)狀態(tài)分布，由于存在目標(biāo)分布協(xié)方差PS與過程噪聲有關(guān)，取決于狀態(tài)運(yùn)動方程，下面主要是確定新生、衍生目標(biāo)分布協(xié)方差Pλ、Pβ的選取規(guī)則。

式(14)中，參與計算狀態(tài)更新增益的殘差協(xié)方差為

(17)

2.2 確定裁剪門限T取值規(guī)則

新生、衍生目標(biāo)方差Pλ、Pβ數(shù)量級與測量方差R相關(guān)，令其對應(yīng)標(biāo)準(zhǔn)差為σ，則高斯元素m服從高斯分布N(x;m,σ2)，其概率密度為

(18)

在高斯元素裁剪處理過程中，當(dāng)高斯元素概率假設(shè)密度滿足式(19)時，此元素就被裁剪，不參與下一步計算

D=ωN(x;m,σ2)

(19)

因此，門限T的取值要保證絕大多數(shù)(99.7%)有效高斯元素能參與計算。但門限T不能太小，否則會引入大量權(quán)值小的雜波元素，雜波元素參與下一步元素合并計算將影響其精度，并且增大計算量；T也不能太大，太大則造成只有少數(shù)高斯元素參與計算，丟失有用信息。式(18)概率密度f對應(yīng)的概率函數(shù)可轉(zhuǎn)換為標(biāo)準(zhǔn)正態(tài)分布的概率函數(shù)[10]

(20)

T取值要保證99.7%高斯元素參與計算，因此

(21)

查標(biāo)準(zhǔn)正態(tài)分布函數(shù)值表可得高斯元素狀態(tài)x取值范圍為x∈(m-3σ,m+3σ)，如圖1所示。

圖1 高斯元素取值范圍

將狀態(tài)x取值范圍代入式(19)，因此可確定裁剪門限T取值規(guī)則

(22)

2.3 殘差協(xié)方差分解計算

在式(14)、(17)中的殘差協(xié)方差S如果是非正定矩陣，其參與增益K的計算時會導(dǎo)致計算發(fā)散，造成濾波失敗，為解決這個問題，可對其進(jìn)行Cholesky和QR分解，令

(23)

再對其進(jìn)行QR分解

QrCholS=qr(CholS)

(24)

則殘差協(xié)方差可轉(zhuǎn)換為

(25)

增益就可通過下式進(jìn)行計算

(26)

由于QrCholS是經(jīng)過QR分解得到的正定上三角矩陣，不會造成計算發(fā)散，因此，解決了計算增益K時的發(fā)散問題。

3 仿真實(shí)驗(yàn)與分析

為了驗(yàn)證本文推導(dǎo)的正確性，采用Matlab軟件仿真的方法來測試驗(yàn)證。仿真硬件環(huán)境為：Pentium(R) Dual-Core CPU，HP E5200，主頻2.5 GHz，2.00 GB的內(nèi)存。仿真場景如下：組網(wǎng)雷達(dá)兩部，分別為雷達(dá)1和雷達(dá)2，雷達(dá)1的距離精度為130 m，方位角精度為0.3°，俯仰角精度為0.2°，配置位置為[118° 29° 120 m]T。雷達(dá)2的距離精度為90 m，方位角精度為0.4°，俯仰角精度為0.1°，配置位置為[117° 31° 20 m]T。融合中心的配置位置為[117° 30° 170 m]T。為了驗(yàn)證計算的簡便性，特別假定兩組網(wǎng)雷達(dá)對目標(biāo)進(jìn)行等間隔(1 s)探測，融合中心每秒能接收到一次探測數(shù)據(jù)，全程探測時間為100 s，所探測的目標(biāo)處于強(qiáng)雜波環(huán)境，并且目標(biāo)數(shù)量未知，探測目標(biāo)全程飛行高度6 000 m不變，探測區(qū)域?yàn)槎S區(qū)域[-105, 105]×[-105105]，目標(biāo)存活概率PS=0.99，存活目標(biāo)服從線性高斯分布，狀態(tài)轉(zhuǎn)移矩陣和過程噪聲方差分別為

(27)

(28)

式中：時間間隔Δt=1s；過程噪聲標(biāo)準(zhǔn)差σ=5m/s2。新生、衍生目標(biāo)都服從泊松分布，其高斯混合形式分別為

(29)

βk+1|k(x|ζ)=0.05N(x;ζ,Qβ)

(30)

Kk(z)=λcVu(z)

(31)

在26s時目標(biāo)1衍生出目標(biāo)3，目標(biāo)3做勻速直線運(yùn)動，速度為V3=[-100 -20]T，濾波裁剪合并高斯元素過程中合并門限U=4，允許參與計算的最多高斯元素J=100。下面，分別對5種取值情況進(jìn)行仿真比較，每種情況殘差協(xié)方差S均應(yīng)用本文的分解處理。跟蹤開始時，將目標(biāo)1、2作為新生目標(biāo)處理，跟蹤以后，跟蹤點(diǎn)跡與目標(biāo)真實(shí)點(diǎn)跡誤差用最優(yōu)子分配(Optimal Subpattern Assignment, OSPA)距離來度量，OSPA距離計算參考文獻(xiàn)[11-12]，實(shí)驗(yàn)中距離上限設(shè)置為200 m。為了克服隨機(jī)因素影響，進(jìn)行100次仿真統(tǒng)計計算。

(1)根據(jù)本文取值規(guī)則，新生、衍生目標(biāo)分布協(xié)方差取值與融合中心測量噪聲虛擬方差R一致，裁剪門限T根據(jù)式(22)計算

Pγ=diag([105, 105, 25, 25]T)

(32)

Pβ=diag([105, 105, 400, 400]T)

(33)

T=10-7

(34)

跟蹤效果圖如圖2、圖3所示。

圖2 GMPHDF跟蹤、目標(biāo)真實(shí)點(diǎn)跡與當(dāng)前時刻測量值

圖2是跟蹤點(diǎn)跡與真實(shí)點(diǎn)跡的OSPA距離圖，OSPA距離均方根為25m，全程計算耗時89.18s。從圖2可以看出，在第26s從目標(biāo)1衍生出目標(biāo)3，在強(qiáng)雜波環(huán)境下，選用本文參數(shù)取值規(guī)則的GMPHDF跟蹤多目標(biāo)，經(jīng)過一定點(diǎn)跡信息積累以后，能跟蹤到繼續(xù)存在目標(biāo)和新衍生目標(biāo)，從圖3可以看出，本文GMPHDF跟蹤誤差較小，且在一定時間內(nèi)，跟蹤誤差方差(OSPA均方根)隨時間增長越來越小，適合于工程實(shí)際應(yīng)用。

圖3 X、Y方向GMPHDF跟蹤點(diǎn)跡與測量值

(2)與條件(1)相比，只增大裁剪門限T，其他條件相同：T=10-6。

與跟蹤條件(1)相比，當(dāng)裁剪門限T取值大于本文取值規(guī)則確定的值時，從圖4可以看出，GMPHDF不能跟蹤到繼續(xù)存在目標(biāo)和新衍生目標(biāo)，跟蹤失敗。

圖4 GMPHDF跟蹤、目標(biāo)真實(shí)點(diǎn)跡與當(dāng)前時刻測量值

(3)只減小裁剪門限T，其他條件與(1)相同：T=10-8。

從圖5可以看出，與條件(1)相比，當(dāng)裁剪門限T減小時，同樣能跟蹤到所有目標(biāo)，但由圖6可知，由于裁剪掉的高斯元素變少，進(jìn)入裁剪以后濾波計算的高速元素增加，其中，包括大量雜波元素，增大了計算負(fù)擔(dān)，計算耗時與條件(1)的89.18s相比，增加到185.61s，并且影響了計算精度，圖6中OSPA距離均方根增加到30m，跟蹤誤差方差并未像條件(1)隨時間增加而減小。其結(jié)果與本文分析一致。

圖5 GMPHDF跟蹤、目標(biāo)真實(shí)點(diǎn)跡與當(dāng)前時刻測量值

圖6 跟蹤點(diǎn)跡與真實(shí)目標(biāo)點(diǎn)跡OSPA距離

(4)只增大新生、衍生目標(biāo)分布協(xié)方差Pγ、Pβ，其他條件與(1)相同

Pγ=diag([106, 106, 25,25]T)

(35)

Pβ=diag([106, 106, 400, 400]T)

(36)

此種條件下，同樣導(dǎo)致跟蹤失敗，其跟蹤效果圖如圖4、圖5所示。

(5)只減小新生、衍生目標(biāo)分布協(xié)方差 ，其他條件與(1)相同

Pγ=diag([104, 104, 25,25]T)

(37)

Pβ=diag([104, 104, 400, 400]T)

(38)

從圖7可以看出，目標(biāo)分布方差減小時，經(jīng)過一定點(diǎn)跡信息積累以后，只能跟蹤到繼續(xù)存在目標(biāo)，不能檢測跟蹤到新衍生目標(biāo)。

圖7 GMPHDF跟蹤、目標(biāo)真實(shí)點(diǎn)跡與當(dāng)前時刻測量值

4 結(jié)束語

針對GMPHDF弱信噪比環(huán)境下跟蹤多目標(biāo)時，目標(biāo)分布方差和裁剪門限按經(jīng)驗(yàn)選取，影響跟蹤效果，以及殘差協(xié)方差參與增益計算容易導(dǎo)致計算發(fā)散問題，從概率統(tǒng)計角度確定了分布方差和裁剪門限選值規(guī)則，通過Cholesky和QR分解保證了殘差協(xié)方差正定性，避免了取逆計算時發(fā)散問題。5種情況仿真對比實(shí)驗(yàn)表明，本文確定的參數(shù)方差P、裁剪門限T、殘差協(xié)方差S計算規(guī)則用于GMPHDF跟蹤弱信噪比多目標(biāo)時，能成功跟蹤到新生、衍生和存在目標(biāo)，跟蹤精度較高，OSPA距離均方根誤差25 m，且均方根誤差一定時間內(nèi)隨濾波步數(shù)增長，越來越小。因此，本文確定的GMPHDF參數(shù)取值規(guī)則科學(xué)合理，適用于跟蹤弱信噪比多目標(biāo)。

[1] 戰(zhàn)立曉，湯子躍，朱振波. 雷達(dá)微弱目標(biāo)檢測前跟蹤算法綜述[J]. 現(xiàn)代雷達(dá), 2013, 35(4): 45-52，57. Zhan Lixiao, Tang Ziyue, Zhu Zhenbo. An overview on track-before-detect algorithms for radar weak targets[J]. Modern Radar,2013,35(4): 45-52,57.

[2] 王文星，印士波，張 幟,等. 雷達(dá)組網(wǎng)應(yīng)用的體系結(jié)構(gòu)研究[J]. 現(xiàn)代雷達(dá), 2014,36(5): 26-30. Wang Wenxing, Yin Shibo, Zhang Zhi,et al. A study on radar network application architecture[J]. Modern Radar, 2014, 36(5): 26-30.

[3] 楊躍輪，李洪梅. 分布式干擾下雷達(dá)組網(wǎng)探測效能分析[J]. 指揮控制與仿真, 2012,34(3): 41-43，66. Yang Yuelun, Li Hongmei. Analysis of radar nets' operational effectiveness on the condition of distributed jammers[J]. Command Control & Simulation,2012,34(3): 41-43,66.

[4] Mahler R P S. Multitarget bayes filtering via first-order multitarget moments[J]. IEEE Transactions on Aerospace and Electronic Systems, 2003, 39(4): 1152-1178.

[5] 吳 偉, 黃 冶. 目標(biāo)起始和維持規(guī)則輔助的SMC-PHD濾波方法[J]. 現(xiàn)代雷達(dá), 2013,35(9): 24-29,35. Wu Wei, Huang Ye. An improved SMC-PHD algorithm assisted by target initiation and continuity rules[J]. Modern Radar, 2013, 35(9): 24-29,35.

[6] Vo B N, Ma W K. The Gaussian mixture probability hypothesis density filter[J]. IEEE Transactions on Signal Processing, 2006, 54(11): 4091-4104.

[7] Li W L, Jia Y M. Gaussian mixture PHD filter for jump Markov models based on best-fitting Gaussian approximation[J]. Signal Processing, 2011, 91(4): 1036-1042.

[8] 劉宗香，謝維信，王 品，等. 一種具有信息保持能力的GM-PHD濾波器[J]. 電子學(xué)報, 2013, 41(8): 1603-1608. Liu Zongxiang, Xie Weixin,Wang Pin, et al. A gaussian mixture PHD filter with the capability of information hold[J]. Acta Electronica Sinica, 2013, 41(8): 1603-1608.

[9] 趙溫波，都基焱. 組網(wǎng)雷達(dá)噪聲慣性坐標(biāo)系誤差統(tǒng)計特性研究[J]. 炮兵學(xué)院學(xué)報, 2010(5): 91-95. Zhao Wenbo, Du Jiyan. Study on statistical properties of radar network noise in inertial coordinate system[J]. Journal of Artillery Academy, 2010(5): 91-95.

[10] 龍永紅. 概率論與數(shù)理統(tǒng)計[M]. 北京：高等教育出版社, 2000. Long Yonghong. Probability theory and mathematical statistics[M]. Beijing: Higher Education Press, 2000.

[11] Hoffman J R, Mahler R P S. Multitarget miss distance via optimal assignment[J]. IEEE Transactions on Systems, Man, and Cybernetics, Part A: Systems and Humans, 2004, 34(3): 327-336.

[12] Schuhmacher D, Vo B T, Vo B N. A consistent metric for performance evaluation of multi-object filters[J]. IEEE Transactions on Signal Processing, 2008,56(8): 3447-3457.

丁海龍 男，1987年生，碩士研究生。研究方向?yàn)槔走_(dá)數(shù)據(jù)處理。

趙溫波 男，1972年生，副教授，碩士生導(dǎo)師。研究方向?yàn)榻y(tǒng)計信號處理。

A Study on Some Key Parameters of GMPHDF in Radar Network

DING Hailong，ZHAO Wenbo

(Amy Officer Academy of PLA, Hefei 230031, China)

Gaussian mixture probability hypothesis density filter(PHDF), which is suitable for tracking weak signal to noise ratio(WSNR) multi-target, has rigorous mathematical foundation. But there is no reasonable calculation rules about distribution covariancePand truncation thresholdTof Gaussian elements in PHDF as yet, which bring bad influence to PHDF. Because the residual covarianceSwill be inverse calculated when it is involved in the gain calculation. IfSis non-positive definite, it would lead to divergence calculation. Probability statistics derivation is used to determine thePandTcalculation rules. Cholesky and QR decomposition is used to solve theScalculation problem. The simulation comparison demonstrates that PHDF in radar network, using the proposed calculation rules ofP,TandS, can precisely track WSNR multi-target, containing exist, birth and spawn targets, and bring no extra calculation burden.

radar network; Gaussian mixture probability hypothesis density filter; parameter

10. 16592/ j. cnki. 1004-7859. 2015. 09. 011

國家自然科學(xué)基金資助項(xiàng)目(61273001)；安徽省自然科學(xué)基金資助項(xiàng)目(11040606M130)

丁海龍 Email：656797226@qq.com

2015-04-26

2015-05-27

TN953

A

1004-7859(2015)08-0044-06