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

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

1 引言

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

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

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

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

2 問題描述

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

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

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

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

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

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

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

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

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

3 基于變分貝葉斯推理的機載雷達機動目標跟蹤算法

3.1 近分層概率圖模型

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

由于式（12）中似然PDFp(yk|xk)為加和的形式，很難直接利用變分貝葉斯方法估計出未知的狀態(tài).于是，將伯努利隨機變量?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 基于變分貝葉斯方法的機載雷達機動目標跟蹤算法設(shè)計

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

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

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

3.3 近似后驗PDF的計算

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

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

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

（a）前向傳播

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

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

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

其中，

（5）當Θ1：K中的元素為ν時，利用式（18）得

其中，

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

3.4 基于變分貝葉斯推理的機動目標跟蹤算法

本文所提用于量測丟失的機動目標跟蹤算法總結(jié)如算法1所示.

4 仿真實驗驗證

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

4.1 算法估計性能驗證

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

其中，K=200表示總的仿真步長.為了驗證所提RVBSD 算法的有效性，本文分別針對上述5 種情況，采用如下的算法進行對比：標準的容積卡爾曼平滑器（Cubature Kalman Smoother，CKS）［21］，基于量測丟失的容積卡爾曼平滑器（Cubature Kalman Smoother for measurement Drop，CKSD）［11］，基于未知過程和量測噪聲協(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)隨機選擇初始的狀態(tài)估計蒙特卡洛迭代次數(shù)和變分迭代次數(shù)N分別為500和10.為了評價上述算法的性能，選用均方根誤差（Root Mean Square Errors，RMSEs）、平均均方根誤差（Average Root Mean Square Errors，ARMSEs）、絕對誤差（Absolute Value Biases，AVBs）以及平均絕對誤差（Averaged Absolute Value Biases，AAVBs）作為性能評價指標［22］.上述所有的狀態(tài)估計算法都在Matlab2020a 仿真平臺上運行，選用的計算機型號為：Intel Core i7-9750H CPU@2.60 GHz.

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

圖2 雷達目標跟蹤軌跡

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

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

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

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

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

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

4.2 不同檢測概率下算法估計性能驗證

為了驗證不同檢測概率下本文算法的估計性能，除情形1外，還采用如下幾種形式的目標檢測概率進行對比分析.

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

表3 不同目標檢測情形下狀態(tài)估計的RMSEs對比

表4 不同目標檢測情形下狀態(tài)估計的AAVBs對比

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

5 結(jié)論

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