董文豪 達凱 宋志勇 付強

（國防科技大學(xué)電子科學(xué)學(xué)院ATR國防重點實驗室，湖南長沙 410073）

1 引言

多普勒雷達通過多普勒效應(yīng)來探測目標(biāo)的位置和相對速度，廣泛應(yīng)用在機載預(yù)警、導(dǎo)彈制導(dǎo)、戰(zhàn)場偵察和天氣預(yù)警等領(lǐng)域中，在強雜波背景和一些干擾的情況下也可實現(xiàn)對目標(biāo)的探測［1］。隨著隱身技術(shù)的不斷發(fā)展，目標(biāo)雷達散射截面積（Radar Cross Section，RCS）呈數(shù)量級遞減趨勢，低信噪比（Signal-to-Noise Ratio，SNR）下的多目標(biāo)檢測與狀態(tài)估計成為了多普勒雷達發(fā)展中的基本任務(wù)和重大挑戰(zhàn)［2］。

傳統(tǒng)的多目標(biāo)檢測估計方法為跟蹤前檢測（Detect Before Track，DBT），在一定的虛警率下，基于Newman-Pearson 準(zhǔn)則利用目標(biāo)與背景間的能量差異來得到目標(biāo)檢測結(jié)果?；跀?shù)據(jù)的傳統(tǒng)多目標(biāo)跟蹤方法主要有聯(lián)合概率數(shù)據(jù)關(guān)聯(lián)（Joint Proba?bilistic Data Association，JPDA）［3］和多假設(shè)跟蹤（Multi-Hypothesis Tracking，MHT）［4］，低SNR 情況下需要降低判決門限來保留目標(biāo)信息，確保目標(biāo)通過門限的同時也帶來了大量虛警，這類方法由于數(shù)據(jù)關(guān)聯(lián)導(dǎo)致計算復(fù)雜度過高而難以有效檢測和估計目標(biāo)。Mahler 開創(chuàng)的有限集統(tǒng)計學(xué)理論（Finite Set Statistics，F(xiàn)ISST）可有效避免數(shù)據(jù)關(guān)聯(lián)，利用隨機有限集（Random Finite Set，RFS）描述多目標(biāo)狀態(tài)和量測，與貝葉斯框架完美融合，發(fā)展出的多目標(biāo)跟蹤方法有概率假設(shè)密度（Probability Hypothesis Den?sity，PHD）［5］、集勢概率假設(shè)密度（Cardinalized PHD，CPHD）［6］和多伯努利濾波器（Multi-Bernoulli，MBR）［7］。PHD 濾波器假定分布為多目標(biāo)泊松過程，通過一階矩近似可完全描述多目標(biāo)分布，CPHD濾波器假定分布為獨立同分布群過程（Independent and Identically Distributed Cluster，IIDC），通過二階矩近似可同時傳播勢分布和PHD，與PHD濾波器相比準(zhǔn)確估計目標(biāo)數(shù)目的同時也增加了計算復(fù)雜度。PHD 和CPHD 濾波器都使用單個密度函數(shù)來表示多個單目標(biāo)狀態(tài)，MBR 濾波器使用多個伯努利RFS對多目標(biāo)狀態(tài)建模，每一個伯努利項利用存在概率和狀態(tài)密度函數(shù)來表示。針對MBR 濾波器的勢偏置問題，Vo［8］提出了勢均衡多伯努利濾波器（Cardi?nality Balanced Multi-target Multi-Bernoulli，CBMeM?Ber）。為了形成目標(biāo)軌跡，Vo［9］提出了廣義標(biāo)簽多伯努利（Generalized Labeled Multi-Bernoulli，GLMB）的可計算近似解，簡化了計算量。FISST 在多傳感器多目標(biāo)跟蹤領(lǐng)域也有新的發(fā)展，文獻［10］通過得到所有目標(biāo)軌跡的概率分布來解決GLMB多傳感器多目標(biāo)場景中的軌跡斷裂問題，有關(guān)多傳感器的最新研究進展可參考文獻［11］。

以上的RFS 濾波器文獻都是基于標(biāo)準(zhǔn)傳感器，標(biāo)準(zhǔn)傳感器是指一個目標(biāo)只能產(chǎn)生一個量測結(jié)果或者不產(chǎn)生量測，且每個量測只能由單個目標(biāo)或者背景產(chǎn)生。疊加式傳感器是指每個量測都以疊加的方式受到多個目標(biāo)影響，即多個目標(biāo)對測量的貢獻等于單獨每個目標(biāo)單獨存在時對量測的貢獻之和，且每個目標(biāo)可以影響任意數(shù)目的量測，因此量測之間不獨立。多普勒雷達回波信號為各個目標(biāo)產(chǎn)生的多普勒回波疊加形成的時域信號，通過對時域回波的分析處理來判斷是否存在目標(biāo)并估計出其多普勒頻率、幅相特征，屬于典型的疊加式傳感器應(yīng)用。除此之外，疊加式傳感器的應(yīng)用場景有陣列天線到達角（Direction of Arrival，DOA）估計［12］、無線通信網(wǎng)絡(luò)多用戶檢測［13］、射頻層分析跟蹤［14］和基于熱源定位的熱電堆陣列［15］等?；诏B加式傳感器的多目標(biāo)RFS 濾波器在預(yù)測步驟與基于疊加式傳感器的經(jīng)典多目標(biāo)RFS 濾波器相同，都是把多目標(biāo)狀態(tài)建模為RFS，而兩種濾波器的更新步驟不同，前者把量測建模為隨機矢量，而后者把量測建模為RFS。馬勒推導(dǎo)了一般疊加式CPHD 的精確濾波方程［16］，但該方程由于計算復(fù)雜而難以實現(xiàn)。針對該問題，Hauschildt 提出了一種近似閉式解［15］，但也只適用于少量密集目標(biāo)情形。Thouin、Nannuru和Coates于2011年提出了加性似然矩（additive like?lihood Moment，ALM）濾波器用于疊加式PHD 濾波器遞推，該近似方法有效但其中的觀測更新方程只適用于離散狀態(tài)空間，馬勒在此基礎(chǔ)上推廣得到了近似疊加式CPHD濾波器［17］，隨后與Nannuru 與Coates 匯總了該濾波器及其特例（即PHD 濾波器），并給出了可計算處理的近似疊加式PHD 和疊加式CPHD 的輔助粒子實現(xiàn)形式［18］。疊加式RFS傳感器在MBR 濾波器也得到了發(fā)展，與疊加式CPHD 濾波器的推導(dǎo)方法類似，Nannuru 和Coates 通過為每一個伯努利RFS 定義一個條件PHD，然后在用與疊加式CPHD 濾波器類似的推導(dǎo)方法來推導(dǎo)疊加式MBR濾波器［19］，并給出了其粒子實現(xiàn)形式［20］。考慮到CPHD 濾波器中的勢分布可準(zhǔn)確估計目標(biāo)數(shù)目，Nannuru 和Coates 推導(dǎo)了混合疊加MBR 和CPHD 濾波器的更新方程及輔助粒子實現(xiàn)，驗證了該濾波器在模擬聲學(xué)傳感器和射頻層分析中的有效性［21］。國內(nèi)基于疊加式傳感器的RFS 多目標(biāo)濾波器研究較少，范紅旗的譯著［22］第19 章介紹了部分相關(guān)進展，但至今還沒有將其應(yīng)用到多普勒雷達目標(biāo)檢測估計領(lǐng)域的相關(guān)文獻。

IIDC RFS 僅使用一個狀態(tài)密度函數(shù)對所有目標(biāo)進行建模，而MBR RFS可單獨對每個目標(biāo)建模并更新其狀態(tài)信息，MBR RFS 不僅結(jié)果簡單，而且更適合處理極端非線性情況，但其缺點是目標(biāo)數(shù)目估計精度不如IIDC RFS 準(zhǔn)確。此外，由于目標(biāo)數(shù)量隨著時間的推移而變化，MBR 濾波器由于需要在每個時間步添加多個伯努利分量來對新生目標(biāo)建模而導(dǎo)致計算效率較低，使用IIDC RFS 對新生目標(biāo)建模，只有在新生目標(biāo)出現(xiàn)時才會導(dǎo)致IIDC RFS 執(zhí)行聚類操作，其他情況則不需要執(zhí)行復(fù)雜度較高的聚類操作。綜合以上因素，本文將多普勒雷達回波信號建模為疊加式傳感器模型，考慮使用混合MBRCPHD 濾波器，給出了貝葉斯框架下該濾波器的預(yù)測和更新方程，并基于輔助粒子濾波方式實現(xiàn)了多普勒雷達多目標(biāo)的聯(lián)合檢測與狀態(tài)估計。具體內(nèi)容安排如下：首先描述了目標(biāo)狀態(tài)RFS 和觀測隨機矢量，緊接著分別介紹了MBR RFS 和IIDC RFS，以及混合MBR-CPHD 濾波器的預(yù)測方程和更新方程。然后給出了該濾波器的輔助粒子濾波實現(xiàn)步驟，并對多普勒雷達回波信號仿真，驗證濾波器的有效性。最后對全文進行總結(jié)。

2 問題描述

2.1 狀態(tài)RFS

多普勒雷達回波信號的一般形式如下：

其中，a為信號幅度，f為由雷達平臺和目標(biāo)徑向運動引起的多普勒頻率，φ為信號初始相位，單目標(biāo)狀態(tài)矢量可用向量x=[φ，f，a]T表示，k時刻目標(biāo)狀態(tài)RFS 可表示為，Nk≥0 表示k時刻的目標(biāo)數(shù)目，單目標(biāo)相位和頻率關(guān)系與距離和速度關(guān)系類似，可由CV 模型來描述，因此相位和頻率關(guān)系可表示為：

其中，λ=c/f0為載頻波長，c為光速，f0為載波頻率，噪聲ωc為有目標(biāo)徑向加速度引起的零均值高斯白噪聲，ωc～Ν(0，)。

假定目標(biāo)幅度不變，即目標(biāo)為非起伏目標(biāo)，采用一階AR模型來描述非起伏目標(biāo)：

其中ωa是零均值高斯白噪聲，分布為ωa～Ν(0，)，則狀態(tài)轉(zhuǎn)移模型表示為：

2.2 觀測隨機矢量

k時刻的觀測矢量可表示為Zk=[zk，1，zk，2，…，]，對于多普勒回波信號，其觀測模型為：

其中，觀測噪聲vk為獨立的復(fù)高斯矢量，其概率密度函數(shù)如下：

其中fs為信號采樣頻率，則似然函數(shù)可表示為：

3 隨機有限集

3.1 MBR RFS

單目標(biāo)伯努利的PHD為：

MBR RFS是N個統(tǒng)計獨立的單目標(biāo)伯努利RFS 的并集，即Ξ=Ξ1∪Ξ2∪...∪ΞN，則多目標(biāo)伯努利RFS的PHD為［7］：

3.2 IIDC RFS

IIDC RFS 用勢分布πc(n)和密度函數(shù)qc(x)表示，IIDC的PHD表示為：

其中μc是勢分布均值，對應(yīng)的二階矩分別為［17］：

3.3 混合MBR-CPHD濾波器

考慮到MBR 濾波器結(jié)構(gòu)簡單，且適用于式（7）這種極端非線性觀測模型，而CPHD 濾波器的勢分布可以用來準(zhǔn)確估計目標(biāo)數(shù)目，因此采用混合MBR-CPHD 濾波器基于貝葉斯準(zhǔn)則實現(xiàn)多普勒雷達多目標(biāo)聯(lián)合檢測與估計，具體步驟如下。

預(yù)測

將k時刻的后驗分布建模為MBR RFS，參數(shù)表示為：

更新

更新步驟基于貝葉斯規(guī)則傳播三個參量，分別為：MBR 分量的PHD、IIDC 分量的PHD 以及IIDC 分量的勢分布，具體表示為：

式（22）至式（25）的推導(dǎo)可參考文獻［17］和文獻［18］。

IIDC分量的PHD更新為：

式（29）和式（30）的推導(dǎo)可參考文獻［17］和文獻［18］。

IIDC分量的勢更新為：

式（33）和式（34）的推導(dǎo)參考文獻［17］和文獻［18］?？梢哉J(rèn)為MBR 濾波器是混合MBR-CPHD 濾波器的特例，通過將混合MBR-CPHD 濾波器中的IIDC 分量置為空集即可得到MBR濾波器。

4 MBR-CPHD輔助粒子濾波實現(xiàn)

由于式（7）的高度非線性，混合MBR-CPHD 濾波器需采用粒子實現(xiàn)形式。粒子濾波器首先需要設(shè)計重要性分布（也叫提議分布）進行粒子采樣，粒子濾波器的經(jīng)典代表自舉（Boot-Strap）濾波器就是把狀態(tài)轉(zhuǎn)移函數(shù)作為重要性分布，由于沒有考慮觀測值，可能導(dǎo)致該濾波器的提議分布與似然分布差別較大，只有極個別采樣粒子位于高似然區(qū)，進而引起粒子退化問題。為此，我們采用輔助粒子濾波器（Auxiliary Particle Filter，APF）［23］，利用觀測值構(gòu)建輔助變量，在此基礎(chǔ)上設(shè)計提議分布可以增加高似然區(qū)域內(nèi)的采樣粒子數(shù)，具體實現(xiàn)步驟如下。

步驟1 提議分布設(shè)置：

設(shè)輔助分布為：

步驟2 預(yù)測：

對于存活目標(biāo)，對應(yīng)RFS 為MBR 分量，基于狀態(tài)轉(zhuǎn)移密度對粒子狀態(tài)進行預(yù)測，即

MBR分量的存在概率預(yù)測為：

對于新生目標(biāo)，對應(yīng)RFS 為IIDC 分量，粒子在狀態(tài)空間中均勻采樣，勢分布為離散的均勻分布。

步驟3 更新：

對粒子權(quán)重進行更新：

步驟4 IIDC RFS近似及重采樣

由于使用伯努利隨機有限集更好地表示單個目標(biāo)，因此使用多個伯努利分量近似更新的IIDC分量，最終的后驗分布使用MBR RFS建模。先利用最大似然估計器估計出IIDC 分量中最大勢，然后利用K 均值算法對IIDC 分量進行聚類，將其劃分為Nck+1個存在概率為1 的伯努利分量，然后對所有的MBR分量重采樣，即

步驟5 航跡管理

設(shè)置存在概率門限τ，當(dāng)MBR 分量中的第i個存在概率若滿足

則刪除第i個伯努利項，反之則保留作為存活目標(biāo)RFS。

5 仿真

多普勒雷達信號的基本參數(shù)設(shè)置如表1所示。

表1 仿真參數(shù)Tab.1 Simulation parameters

仿真時長為50 幀，每一幀的觀測矢量維度為L=5ms×50kHz=250，觀測噪聲協(xié)方差矩陣∑z=2 ?IL×L，共設(shè)置三個目標(biāo)，目標(biāo)1始終存在，目標(biāo)2在第9幀出現(xiàn)，第35幀消失，目標(biāo)3在第17幀出現(xiàn)，第24 幀消失?；旌螹BR-CPHD 濾波器的最大目標(biāo)數(shù)設(shè)為15，每個目標(biāo)粒子數(shù)為1000，粒子存活概率和新生概率分別設(shè)置為0.9、0.1，且均與粒子狀態(tài)無關(guān)，過程噪聲方差設(shè)置為σc=0.5，σa=0.5，存在概率門限設(shè)為τ=0.1。

當(dāng)信噪比為-4 dB 時，MBR 濾波器和混合MBR-CPHD 濾波器的單次運行結(jié)果如圖1 至圖4所示。

從圖1 至圖3 可以看出，三個目標(biāo)的相位在中間時刻相互交叉，即使目標(biāo)相位在仿真時間段內(nèi)不斷跳變，MBR 濾波器和混合MBR-CPHD 濾波器幾乎都能準(zhǔn)確估計出不同時刻不同目標(biāo)的真實相位；由于假設(shè)狀態(tài)模型為CV 模型，所以不同目標(biāo)頻率保持不變，三個目標(biāo)的信號幅度相同，MBR 濾波器和混合MBR-CPHD 濾波器也都可以較為準(zhǔn)確地估計出目標(biāo)的真實頻率和真實幅度，但MBR濾波器存在一定的延遲效應(yīng)，前兩次目標(biāo)新生時MBR濾波器會延遲1～3 幀。這是由于CPHD 可以通過勢分布準(zhǔn)確估計出目標(biāo)數(shù)目，所以混合MBR-CPHD 濾波器可以準(zhǔn)確估計目標(biāo)新生和消失，并準(zhǔn)確估計出對應(yīng)的目標(biāo)狀態(tài)，正如圖4所示。

為了有效評估RFS 多目標(biāo)濾波器的狀態(tài)估計精度與目標(biāo)數(shù)目估計精度，采用最優(yōu)子模式分配度量（Optimal Sub Pattern Assignment，OSPA）［7］作為性能評價指標(biāo)，其中的勢代價因子c決定目標(biāo)數(shù)目估計精度相對于相對狀態(tài)估計精度的重要性，靈敏度因子p決定濾波器對虛警或者漏檢的敏感度。由于目標(biāo)狀態(tài)變量量綱不同，單獨統(tǒng)計不同信噪比下不同目標(biāo)狀態(tài)的OSPA，分別做100 次蒙特卡洛實驗，結(jié)果分別如圖5 至圖7 所示。隨著信噪比的提高，MBR 濾波器和混合MBR-CPHD 的性能都有所提升，但混合MBR-CPHD 濾波器對每一個目標(biāo)狀態(tài)變量的估計性能幾乎都優(yōu)于MBR濾波器，僅當(dāng)信噪比為1 dB 時MBR-CPHD 濾波器的頻率OSPA 統(tǒng)計值略高于MBR 濾波器，可能原因是本身MBR 濾波器相比CPHD 濾波器更適用于本文高度非線性應(yīng)用，隨著信噪比的提高，導(dǎo)致混合MBR-CPHD 的估計性能改善略低于MBR濾波器。從整體對比發(fā)現(xiàn)，MBR濾波器受信噪比影響更大，當(dāng)信噪比為-4 dB 時，其各個狀態(tài)變量的OSPA 距離均遠大于混合MBRCPHD濾波器。

為了比較MBR 濾波器和混合MBR-CPHD 濾波器的計算復(fù)雜度，進行100次蒙特卡洛實驗，統(tǒng)計每次蒙特卡洛試驗的執(zhí)行時間直方圖對比如圖8 所示。MBR 濾波器中的新生步驟由于需要多個伯努利項實現(xiàn)，導(dǎo)致MBR 分量計算復(fù)雜度提高，CPHD濾波器中的粒子聚類步驟同樣也需要大量的計算，但混合MBR-CPHD 濾波器中只有在目標(biāo)新生時才會有聚類步驟，當(dāng)目標(biāo)數(shù)目不變或者消亡時不需要執(zhí)行此步驟，而MBR 濾波器時刻需要執(zhí)行新生步驟，因此混合MBR-CPHD 濾波器的執(zhí)行時間小于MBR濾波器。

6 結(jié)論

為了避免低信噪比下的目標(biāo)信息丟失，本文直接基于原始多普勒雷達回波信號并將其建模為疊加式傳感器模型，利用混合MBR-CPHD 濾波器估計多普勒雷達多目標(biāo)的個數(shù)及其狀態(tài)，并給出了該濾波器的APF 實現(xiàn)流程。仿真實驗結(jié)果表明，與基于疊加式傳感器的MBR 濾波器相比，混合MBR-CPHD 濾波器在低信噪比下的目標(biāo)數(shù)目和狀態(tài)估計精度更高，同時降低了計算復(fù)雜度。