強(qiáng)雜波背景下基于變分貝葉斯推理的機(jī)載雷達(dá)目標(biāo)跟蹤算法

李淑慧，鄧志紅，馮肖雪，潘 峰，2

（1.北京理工大學(xué)自動(dòng)化學(xué)院，北京 100081；2.昆明北理工產(chǎn)業(yè)技術(shù)研究院有限公司，云南昆明 650000）

1 引言

機(jī)載雷達(dá)作為一種對(duì)地觀測(cè)的重要傳感器，具有全天候、全天時(shí)、高分辨的特點(diǎn)，能實(shí)現(xiàn)目標(biāo)的探測(cè)、識(shí)別和估計(jì)等功能，因此受到軍事領(lǐng)域的廣泛關(guān)注.但是，機(jī)載雷達(dá)常常工作于復(fù)雜多變的惡劣環(huán)境，使得目標(biāo)跟蹤面臨的不確定性問題異常突出：（1）機(jī)載雷達(dá)一般處于平視或下視工作狀態(tài)，必然受到比地基雷達(dá)更為嚴(yán)重的地（海）雜波［1］.雷達(dá)雜波的脈沖干擾會(huì)使機(jī)載雷達(dá)系統(tǒng)的量測(cè)噪聲呈現(xiàn)復(fù)雜、未知的非高斯長(zhǎng)拖尾特性［2，3］.（2）實(shí)際作戰(zhàn)環(huán)境中，機(jī)載雷達(dá)跟蹤的敵方目標(biāo)（如坦克、戰(zhàn)車、無人機(jī)等）往往具有較強(qiáng)的機(jī)動(dòng)性.敵方目標(biāo)的強(qiáng)機(jī)動(dòng)將會(huì)使目標(biāo)跟蹤模型中的狀態(tài)量發(fā)生突變，進(jìn)而誘導(dǎo)未知的非高斯長(zhǎng)拖尾系統(tǒng)噪聲［4］.（3）載機(jī)的運(yùn)動(dòng)導(dǎo)致其雜波分布范圍廣、強(qiáng)度大，尤其在城市和山區(qū)地帶，副瓣雜波通常會(huì)淹沒落入雜波區(qū)的目標(biāo)，使機(jī)載雷達(dá)無法檢測(cè)到目標(biāo)，出現(xiàn)隨機(jī)的目標(biāo)航跡丟失（即量測(cè)丟失）現(xiàn)象［5］.總之，由于復(fù)雜的作戰(zhàn)環(huán)境以及目標(biāo)機(jī)動(dòng)特性帶來的不確定性因素給機(jī)載雷達(dá)目標(biāo)跟蹤帶來了巨大的挑戰(zhàn).因此，研究含量測(cè)丟失和未知長(zhǎng)拖尾非高斯噪聲的機(jī)載雷達(dá)目標(biāo)跟蹤問題具有十分重要的意義.

目前，已有大量的文獻(xiàn)研究含未知噪聲的狀態(tài)估計(jì)問題［6］.例如，Ardeshiri 等［6］針對(duì)過程噪聲協(xié)方差和量測(cè)噪聲協(xié)方差均未知的狀態(tài)估計(jì)問題，采用變分貝葉斯方法設(shè)計(jì)了近似貝葉斯平滑器.然而，該算法只適用于高斯噪聲場(chǎng)景.為此，一些學(xué)者也針對(duì)含未知長(zhǎng)拖尾非高斯噪聲的非線性系統(tǒng)狀態(tài)估計(jì)問題開展了研究［7～10］.趙俊波等［8］針對(duì)含未知過程噪聲協(xié)方差和量測(cè)噪聲協(xié)方差的非線系統(tǒng)，設(shè)計(jì)了魯棒的無跡卡爾曼估計(jì)器.黃玉龍等［10］針對(duì)含長(zhǎng)拖尾過程和測(cè)量噪聲的非線性系統(tǒng)，提出了魯棒的狀態(tài)平滑器.該算法憑借變分推理逼近真實(shí)后驗(yàn)PDF（Probability Density Function）的優(yōu)勢(shì)，聯(lián)合估計(jì)出長(zhǎng)拖尾非高斯噪聲的自由度、尺度矩陣和系統(tǒng)狀態(tài).上述算法均不適用量測(cè)丟失的場(chǎng)景.

近年來，研究人員也針對(duì)含量測(cè)丟失的非線性系統(tǒng)狀態(tài)估計(jì)問題進(jìn)行研究［11～14］.馬柯茂等［11］針對(duì)含高斯噪聲的非線性系統(tǒng)，采用量測(cè)丟失補(bǔ)償?shù)姆椒ㄔO(shè)計(jì)了用于量測(cè)丟失的狀態(tài)估計(jì)方法.然而，該算法難以應(yīng)用于長(zhǎng)拖尾非高斯噪聲場(chǎng)景.江順等［12］針對(duì)含延時(shí)和隨機(jī)量測(cè)丟失的網(wǎng)絡(luò)化非線性系統(tǒng)，設(shè)計(jì)了故障估計(jì)器.雖然該算法適用于非高斯噪聲，但它要求噪聲協(xié)方差有界且已知，量測(cè)丟失概率恒定已知，且狀態(tài)估計(jì)的精度較低.上述算法只適用于噪聲已知的場(chǎng)景，尚不能處理量測(cè)丟失概率和噪聲均未知的情況.

由于機(jī)動(dòng)目標(biāo)跟蹤場(chǎng)景復(fù)雜，未知的長(zhǎng)拖尾非高斯噪聲和量測(cè)丟失不可避免地同時(shí)存在于真實(shí)系統(tǒng)中.現(xiàn)有的算法均無法解決同時(shí)含量測(cè)丟失和未知長(zhǎng)拖尾非高斯噪聲的狀態(tài)估計(jì)問題.為此，本文研究強(qiáng)雜波背景下含量測(cè)丟失的機(jī)動(dòng)目標(biāo)跟蹤問題.

2 問題描述

考慮到復(fù)雜的作戰(zhàn)環(huán)境、強(qiáng)雜波以及目標(biāo)強(qiáng)機(jī)動(dòng)等因素的影響，本文研究了如下的機(jī)載雷達(dá)目標(biāo)跟蹤系統(tǒng)［4，15，16］：

其中，xk∈Rn是k時(shí)刻目標(biāo)的狀態(tài)，yk∈Rm是笛卡爾坐標(biāo)系下k時(shí)刻雷達(dá)接收到的測(cè)量，fk-1(·)和hk(·)分別為已知的函數(shù).表示直接利用全球定位系統(tǒng)（Global Positioning System，GPS）得到機(jī)載平臺(tái)的位置.隨機(jī)變量?k服從概率未知的伯努利隨機(jī)分布，表示量測(cè)數(shù)據(jù)的丟失情況.假設(shè)初始狀態(tài)向量x0服從均值為協(xié)方差為P0|0的高斯分布，且與系統(tǒng)噪聲和量測(cè)噪聲互不相關(guān).

由于閃爍噪聲服從長(zhǎng)拖尾的非高斯分布［15］，同時(shí)考慮到目標(biāo)的閃爍、雜波、強(qiáng)電磁干擾、以及目標(biāo)強(qiáng)機(jī)動(dòng)等因素的影響，本文采用閃爍噪聲來模擬未知的長(zhǎng)拖尾過程噪聲∈Rn和量測(cè)噪聲∈Rm.此外，由于學(xué)生t分布可以視為無窮多個(gè)高斯分布與長(zhǎng)拖尾Gamma分布的混合，具有長(zhǎng)拖尾非高斯特性，常將閃爍噪聲建模為學(xué)生t分布［17］.于是分別寫成如下等式：

由于長(zhǎng)拖尾的非高斯噪聲未知，常常無法獲得準(zhǔn)確的尺度矩陣和自由度參數(shù).假設(shè)尺度矩陣Q，R的先驗(yàn)分布分別服從逆Wishart分布［10］，即

其中，t0，u0分別表示自由度參數(shù)，T0，U0分別為逆尺度參數(shù).未知的自由度參數(shù)ω，ν分別服從Gamma先驗(yàn)分布，即

此外，x0，Q，R，ω，v的聯(lián)合先驗(yàn)概率分布服從如下等式：

此外，由于地（海）雜波以及載機(jī)運(yùn)動(dòng)的影響，副瓣雜波通常會(huì)淹沒落入雜波區(qū)的目標(biāo)，使機(jī)載雷達(dá)無法檢測(cè)的目標(biāo)，出現(xiàn)隨機(jī)的量測(cè)丟失現(xiàn)象.當(dāng)?k=1 時(shí)，表示k時(shí)刻雷達(dá)檢測(cè)到目標(biāo).當(dāng)?k=0 時(shí)，表示k時(shí)刻雷達(dá)丟失目標(biāo).隨機(jī)變量?k滿足如下PDF：

其中，e0∈[0，1]和e1∈[0，1]為Beta 分布的先驗(yàn)形狀參數(shù).

本文的主要目的在于針對(duì)含量測(cè)丟失和未知長(zhǎng)拖尾非高斯噪聲的機(jī)載雷達(dá)目標(biāo)跟蹤系統(tǒng)式（1），利用變分貝葉斯方法設(shè)計(jì)出強(qiáng)雜波背景下含量測(cè)丟失的機(jī)載雷達(dá)機(jī)動(dòng)目標(biāo)跟蹤算法，進(jìn)而估計(jì)出目標(biāo)的狀態(tài).

3 基于變分貝葉斯推理的機(jī)載雷達(dá)機(jī)動(dòng)目標(biāo)跟蹤算法

3.1 近分層概率圖模型

本節(jié)的目的是設(shè)計(jì)分層概率圖模型，為后續(xù)變分貝葉斯推理做準(zhǔn)備.當(dāng)量測(cè)丟失概率已知時(shí)，似然概率密度p(yk|xk)為

由于式（12）中似然PDFp(yk|xk)為加和的形式，很難直接利用變分貝葉斯方法估計(jì)出未知的狀態(tài).于是，將伯努利隨機(jī)變量?k的概率轉(zhuǎn)化為乘積形式的概率質(zhì)量函數(shù)，

于是，似然PDFp(yk|xk)可寫為

根據(jù)式（1）和式（2）得

根據(jù)式（14）和式（15）得

于是，式（4）～式（7）、式（11）、式（13）、式（16）和式（17）共同構(gòu)成了如圖1所示的分層概率圖模型.

圖1 分層概率圖模型

3.2 基于變分貝葉斯方法的機(jī)載雷達(dá)機(jī)動(dòng)目標(biāo)跟蹤算法設(shè)計(jì)

因?yàn)槭剑?8）表達(dá)形式復(fù)雜，所以聯(lián)合PDFp(Θ1：K，y1：K)不存在解析解.本文利用平均場(chǎng)理論［18］和變分貝葉斯方法，通過尋找一個(gè)解析、易分解的近似PDF來逼近真實(shí)的平滑后驗(yàn)PDF，即

令θk為Θ1：K中的任意元素，即θk∈Θ1：K.式（19）中的q(θk)是θk的近似后驗(yàn)PDF.通過最小化近似平滑后驗(yàn)PDFq(Θ1：K) 與真實(shí)后驗(yàn) PDFp(Θ1：K|y1：K) 之間的Kullback-Leibler散度就可以獲得q(θk)的表達(dá)形式［10］，即

其中，E(·) 表示期望，cθ是與θk無關(guān)的常數(shù).由于q(λ1：K)的變分參數(shù)相互耦合，不能直接求解式（20），而需要利用固定點(diǎn)迭代法求解式（20）.固定點(diǎn)迭代法在迭代更新中的任意PDF 時(shí)，需要使用其他PDF 上一次的迭代值［18］.

3.3 近似后驗(yàn)PDF的計(jì)算

為了計(jì)算近似的平滑后驗(yàn)PDF，對(duì)式（18）取對(duì)數(shù)，即log(p(Θ1：K，y1：K)).利用式（20）可依次計(jì)算各個(gè)近似后驗(yàn)PDFq(θk)的變分參數(shù).

（1）當(dāng)Θ1：K中的元素為x0：K時(shí)，利用式（18）和式（20），q(i+1)(x0：K)計(jì)算如下：

根據(jù)式（22）可知，q(i+1)(x0：K)等同于具有修正似然PDFN(yk；hk(xk)，)的非線性后驗(yàn)平滑PDF.于是，將含有量測(cè)丟失的非線性系統(tǒng)就轉(zhuǎn)化為不含量測(cè)丟失等價(jià)的非線性高斯系統(tǒng).基于高斯平滑器設(shè)計(jì)原理［19］，狀態(tài)的更新過程主要包括前向傳播和后向傳播.因此，可得如下等式：

（a）前向傳播

其中，tr(·)表示跡操作.

令θ=?1：K，利用式（18）得

根據(jù)式（38）和式（40），q(i+1)(?1：K)可以更新為伯努利分布.隨機(jī)變量?k取0和1時(shí)的概率分別為

其中，

（5）當(dāng)Θ1：K中的元素為ν時(shí)，利用式（18）得

其中，

接下來，計(jì)算求解近似后驗(yàn)PDF 時(shí)所需的期望.根據(jù)Gamma分布、伯努利分布和Beta分布的性質(zhì)可得

3.4 基于變分貝葉斯推理的機(jī)動(dòng)目標(biāo)跟蹤算法

本文所提用于量測(cè)丟失的機(jī)動(dòng)目標(biāo)跟蹤算法總結(jié)如算法1所示.

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

以機(jī)載雷達(dá)跟蹤空中目標(biāo)為例，驗(yàn)證本文所提RVBSD 算法的有效性.參考文獻(xiàn)［15，16，20］設(shè)置了本文的仿真場(chǎng)景.假設(shè)空中目標(biāo)做恒速轉(zhuǎn)彎模型運(yùn)動(dòng)，即

4.1 算法估計(jì)性能驗(yàn)證

本小節(jié)設(shè)置如下形式的目標(biāo)檢測(cè)概率：

其中，K=200表示總的仿真步長(zhǎng).為了驗(yàn)證所提RVBSD 算法的有效性，本文分別針對(duì)上述5 種情況，采用如下的算法進(jìn)行對(duì)比：標(biāo)準(zhǔn)的容積卡爾曼平滑器（Cubature Kalman Smoother，CKS）［21］，基于量測(cè)丟失的容積卡爾曼平滑器（Cubature Kalman Smoother for measurement Drop，CKSD）［11］，基于未知過程和量測(cè)噪聲協(xié)方差的變分貝葉斯容積卡爾曼平滑器（Variational Bayesian based Cubature Kalman Smoother with unknown parameters QR，VBCKS_QR）［6］，基于未知自由度參數(shù)的變分貝葉斯容積卡爾曼平滑器（Variational Bayesian based Cubature Kalman Smoother with unknown parameters wv，VBCKS_wv）［10］，基于噪聲參數(shù)未知的變分貝葉斯平滑器（Variational Bayesian based Cubature Kalman Smoother with unknown parameters QRwv，VBCKS_QRwv）［10］.CKS 和CKSD 方法的噪聲協(xié)方差分別設(shè)為VBCKS_QR 方法的初始噪聲參數(shù)為t0=6，T0=u0=4，U0=VBCKS_QRwv 方法的初始參數(shù)為s0=r0=g0=3 且h0=3.本文所提RVBSD 算法初始參數(shù)值分別為r0=h0=g0=s0=3，T0=U0=t0=6，u0=4，e1=e0=0.5.根據(jù)高斯分布N(x0，P0|0)隨機(jī)選擇初始的狀態(tài)估計(jì)蒙特卡洛迭代次數(shù)和變分迭代次數(shù)N分別為500和10.為了評(píng)價(jià)上述算法的性能，選用均方根誤差（Root Mean Square Errors，RMSEs）、平均均方根誤差（Average Root Mean Square Errors，ARMSEs）、絕對(duì)誤差（Absolute Value Biases，AVBs）以及平均絕對(duì)誤差（Averaged Absolute Value Biases，AAVBs）作為性能評(píng)價(jià)指標(biāo)［22］.上述所有的狀態(tài)估計(jì)算法都在Matlab2020a 仿真平臺(tái)上運(yùn)行，選用的計(jì)算機(jī)型號(hào)為：Intel Core i7-9750H CPU@2.60 GHz.

圖2 展示了采用上述6 種算法所得的目標(biāo)軌跡跟蹤示意圖.圖3 分別描述了上述6 種算法的位置RMSEs、速度RMSEs 以及加速度RMSEs 結(jié)果圖.圖4 是與之相對(duì)應(yīng)的位置AVBs、速度AVBs 以及加速度AVBs 結(jié)果圖.表1 和表2 分別是基于500 次蒙特卡洛仿真得出的位置ARMSEs、速度ARMSEs、加速度ARMSEs、位置AAVBs、速度AAVBs 以及加速度AAVBs.此外，表1 還給出了上述6種算法的計(jì)算代價(jià).

圖2 雷達(dá)目標(biāo)跟蹤軌跡

從圖2可以看出，本文所提RVBSD算法估計(jì)的軌跡更貼近真實(shí)的目標(biāo)運(yùn)動(dòng)軌跡.此外，從圖3、圖4可以看出，本文所提的RVBSD算法具有較小的位置RMSEs、速度RMSEs、加速度RMSEs、位置AVBs、速度AVBs以及加速度AVBs.需要特別指出的是，當(dāng)采用標(biāo)準(zhǔn)的CKS算法估計(jì)含量測(cè)丟失和長(zhǎng)拖尾非高斯噪聲的系統(tǒng)狀態(tài)時(shí)，不可避免地會(huì)發(fā)生濾波發(fā)散.濾波發(fā)散主要的原因在于標(biāo)準(zhǔn)的CKS 算法只適用于高斯系統(tǒng)且對(duì)丟失的量測(cè)未做任何的處理.雖然基于量測(cè)丟失補(bǔ)償?shù)腃KSD算法采用了量測(cè)預(yù)測(cè)補(bǔ)償?shù)牟呗裕怯捎趯⑾到y(tǒng)噪聲等價(jià)為高斯噪聲，進(jìn)行狀態(tài)估計(jì)時(shí)并不能很好地刻畫非高斯噪聲的長(zhǎng)拖尾特性.與CKS 算法和CKSD 算法相比，剩余的三種對(duì)比算法均是將系統(tǒng)噪聲建模為學(xué)生t分布.在某種程度上，它們能刻畫非高斯噪聲的長(zhǎng)拖尾特性.但是VBCKS_QR 算法和VBCKS_wv算法分別針對(duì)未知學(xué)生t分布的尺度參數(shù)和自由度參數(shù)進(jìn)行估計(jì)，不適用于二者均未知的場(chǎng)景.雖然VBCKS_QRwv 算法可以聯(lián)合估計(jì)未知學(xué)生t 分布的尺度和自由度參數(shù)，但與VBCKS_QR算法和VBCKS_wv 算法類似，未考慮量測(cè)丟失的情況.而本文所提的RVBSD 算法不僅能聯(lián)合估計(jì)未知學(xué)生t分布中的自由度參數(shù)和尺度參數(shù)，還考慮了量測(cè)丟失對(duì)狀態(tài)估計(jì)的影響.綜上所述，本文所提的RVBSD算法能更好地跟蹤目標(biāo)的軌跡.

圖3 狀態(tài)估計(jì)RMSEs對(duì)比

圖4 狀態(tài)估計(jì)AVBs對(duì)比

此外，如表1 和表2 所示，本文所提RVBSD 算法有較小的位置ARMSEs、速度ARMSEs、加速度ARMSEs、位置AAVBs、速度AAVBs 以及加速度AAVBs.這也驗(yàn)證了本文所提算法的有效性.雖然本文所提的RVBSD算法與其他5種算法相比耗時(shí)相對(duì)較長(zhǎng)，但是為了獲得較高的估計(jì)精度，犧牲一定的時(shí)間成本也是允許的.總體而言，本文所提RVBSD算法是有效的.

表1 狀態(tài)估計(jì)的RMSEs對(duì)比

表2 狀態(tài)估計(jì)的AAVBs對(duì)比

4.2 不同檢測(cè)概率下算法估計(jì)性能驗(yàn)證

為了驗(yàn)證不同檢測(cè)概率下本文算法的估計(jì)性能，除情形1外，還采用如下幾種形式的目標(biāo)檢測(cè)概率進(jìn)行對(duì)比分析.

鑒于4.1 節(jié)已經(jīng)驗(yàn)證本文算法的優(yōu)勢(shì)，為節(jié)省篇幅，本小節(jié)省略了其他算法的結(jié)果分析.除了目標(biāo)檢測(cè)概率，其他參數(shù)的設(shè)置均同4.1節(jié).表3和表4分別給出了不同目標(biāo)檢測(cè)情形下位置ARMSEs、速度ARMSEs、加速度ARMSEs、位置AAVBs、速度AAVBs 以及加速度AAVBs的統(tǒng)計(jì)結(jié)果.

表3 不同目標(biāo)檢測(cè)情形下狀態(tài)估計(jì)的RMSEs對(duì)比

表4 不同目標(biāo)檢測(cè)情形下狀態(tài)估計(jì)的AAVBs對(duì)比

首先，值得肯定的是，本文所提的算法在情形1～3下取得了較好的跟蹤結(jié)果.這進(jìn)一步驗(yàn)證了本文所提算法的有效性.此外，從目標(biāo)檢測(cè)概率πk在K/3 ≤k＜2K/3區(qū)間取不同概率值來看，隨著目標(biāo)檢測(cè)概率πk的減少，狀態(tài)估計(jì)的性能大致呈下降趨勢(shì).從相同的目標(biāo)檢測(cè)概率πk在不同區(qū)間的取值來看，越長(zhǎng)的采樣區(qū)間遭受強(qiáng)雜波（即較低的目標(biāo)檢測(cè)概率），狀態(tài)估計(jì)的性能總體呈下降的趨勢(shì).需要特別注意的是，當(dāng)目標(biāo)檢測(cè)概率分別取情形4和情形6時(shí)，都出現(xiàn)了濾波發(fā)散問題.這主要是因?yàn)殡S著目標(biāo)檢測(cè)概率的降低，收集的量測(cè)中包含有用的狀態(tài)估計(jì)信息缺失嚴(yán)重.同樣，當(dāng)目標(biāo)檢測(cè)概率受強(qiáng)雜波影響的區(qū)間越多，強(qiáng)雜波淹沒目標(biāo)的有用信息越多，進(jìn)而使雷達(dá)傳感器難以接收到有效的量測(cè)對(duì)狀態(tài)進(jìn)行估計(jì)，出現(xiàn)了濾波發(fā)散問題.因此，為了有效的進(jìn)行目標(biāo)狀態(tài)的估計(jì)，目標(biāo)的檢測(cè)概率不能太低，遭受強(qiáng)雜波的時(shí)間（即區(qū)間）也不能太長(zhǎng)，否則無法收集到目標(biāo)的有用信息對(duì)狀態(tài)進(jìn)行估計(jì).總體上來說，本文提出的算法能魯棒地、有效地應(yīng)對(duì)強(qiáng)雜波的干擾，對(duì)目標(biāo)的狀態(tài)進(jìn)行估計(jì).

5 結(jié)論

本文設(shè)計(jì)了強(qiáng)雜波背景下含量測(cè)丟失的目標(biāo)跟蹤算法.該算法采用學(xué)生t分布來模擬非高斯噪聲的長(zhǎng)拖尾特性.通過引入伯努利隨機(jī)變量，將求和形式的后驗(yàn)概率密度函數(shù)轉(zhuǎn)換成乘積形式的概率質(zhì)量函數(shù)，并構(gòu)建了分層狀態(tài)空間模型.在此基礎(chǔ)上，設(shè)計(jì)了用于量測(cè)丟失的魯棒變分貝葉斯平滑器.以機(jī)載雷達(dá)跟蹤空中目標(biāo)為例驗(yàn)證了本文算法的有效性.