董文豪 達(dá)凱 宋志勇 付強(qiáng)

（國(guó)防科技大學(xué)電子科學(xué)學(xué)院ATR國(guó)防重點(diǎn)實(shí)驗(yàn)室，湖南長(zhǎng)沙 410073）

1 引言

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

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

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

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

2 問題描述

2.1 狀態(tài)RFS

多普勒雷達(dá)回波信號(hào)的一般形式如下：

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

其中，λ=c/f0為載頻波長(zhǎng)，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 觀測(cè)隨機(jī)矢量

k時(shí)刻的觀測(cè)矢量可表示為Zk=[zk，1，zk，2，…，]，對(duì)于多普勒回波信號(hào)，其觀測(cè)模型為：

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

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

3 隨機(jī)有限集

3.1 MBR RFS

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

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

3.2 IIDC RFS

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

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

3.3 混合MBR-CPHD濾波器

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

預(yù)測(cè)

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

更新

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

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

IIDC分量的PHD更新為：

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

IIDC分量的勢(shì)更新為：

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

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

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

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

設(shè)輔助分布為：

步驟2 預(yù)測(cè)：

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

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

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

步驟3 更新：

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

步驟4 IIDC RFS近似及重采樣

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

步驟5 航跡管理

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

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

5 仿真

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

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

仿真時(shí)長(zhǎng)為50 幀，每一幀的觀測(cè)矢量維度為L(zhǎng)=5ms×50kHz=250，觀測(cè)噪聲協(xié)方差矩陣∑z=2 ?IL×L，共設(shè)置三個(gè)目標(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，每個(gè)目標(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 時(shí)，MBR 濾波器和混合MBR-CPHD 濾波器的單次運(yùn)行結(jié)果如圖1 至圖4所示。

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

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

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

6 結(jié)論

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