何祥宇，李 靜，楊數(shù)強(qiáng)，夏玉杰

（1.洛陽(yáng)師范學(xué)院物理與電子信息學(xué)院，河南洛陽(yáng) 471934；2.洛陽(yáng)師范學(xué)院信息技術(shù)學(xué)院，河南洛陽(yáng) 471934）

（?通信作者電子郵箱xyuhe@foxmail.com）

0 引言

目標(biāo)跟蹤技術(shù)是一個(gè)備受國(guó)內(nèi)外學(xué)者關(guān)注的研究課題，應(yīng)用范圍十分廣泛［1-3］。多目標(biāo)跟蹤（Multi-Target Tracking，MTT）技術(shù)研究是目標(biāo)跟蹤技術(shù)領(lǐng)域內(nèi)的一個(gè)重要研究分支，多目標(biāo)跟蹤的目的是利用受雜波和噪聲影響的傳感器測(cè)量值估計(jì)目標(biāo)的個(gè)數(shù)及狀態(tài)。隨機(jī)有限集（Random Finite Set，RFS）理論被認(rèn)為是處理多目標(biāo)跟蹤問(wèn)題的一種簡(jiǎn)潔高效的方法，通過(guò)把多目標(biāo)狀態(tài)集和測(cè)量集定義為隨機(jī)有限集，該理論建立了嚴(yán)謹(jǐn)?shù)淖顑?yōu)貝葉斯多目標(biāo)濾波框架。

在基于RFS 理論的多目標(biāo)跟蹤方法中，概率假設(shè)密度（Probability Hypothesis Density，PHD）濾波器［4］是一種很有潛力的多目標(biāo)貝葉斯濾波算法，它不用復(fù)雜的數(shù)據(jù)關(guān)聯(lián)過(guò)程就可以處理復(fù)雜環(huán)境中的多目標(biāo)跟蹤問(wèn)題。標(biāo)準(zhǔn)的PHD 濾波器是基于點(diǎn)目標(biāo)假設(shè)的多目標(biāo)跟蹤算法，點(diǎn)目標(biāo)假設(shè)即：在每一個(gè)時(shí)間步，每一個(gè)目標(biāo)假設(shè)只能產(chǎn)生至多一個(gè)量測(cè)值。與點(diǎn)目標(biāo)不同的是，一個(gè)擴(kuò)展目標(biāo)在每一個(gè)時(shí)間步上可產(chǎn)生多個(gè)量測(cè)值［5］。因此，應(yīng)用標(biāo)準(zhǔn)的PHD 濾波器跟蹤多個(gè)擴(kuò)展目標(biāo)會(huì)產(chǎn)生不可靠的跟蹤結(jié)果。

為在RFS理論框架內(nèi)處理多擴(kuò)展目標(biāo)跟蹤問(wèn)題，文獻(xiàn)［6］中提出了標(biāo)準(zhǔn)PHD 濾波器的一個(gè)擴(kuò)展版本，稱之為擴(kuò)展目標(biāo)PHD（Extended Target PHD，ET-PHD）濾波器。目前，ET-PHD濾波器實(shí)現(xiàn)方法主要有兩種，分別是擴(kuò)展目標(biāo)序貫蒙特卡洛PHD（Extended Target Sequential Monte Carlo PHD，ET-SMCPHD）濾波器［7］和擴(kuò)展目標(biāo)高斯混合PHD（Extended Target Gaussian Mixture PHD，ET-GM-PHD）濾波器［8］。ET-SMC-PHD濾波器使用隨機(jī)樣本粒子的加權(quán)和形式近似后驗(yàn)強(qiáng)度函數(shù)，利用聚類及其優(yōu)化算法實(shí)現(xiàn)目標(biāo)狀態(tài)的抽取與估計(jì)，是ETPHD 濾波器一種通用實(shí)現(xiàn)。ET-GM-PHD 濾波器是ET-PHD濾波器在線性高斯假設(shè)條件下的封閉解，利用高斯混合函數(shù)描述目標(biāo)的后驗(yàn)強(qiáng)度函數(shù)。此外，近期提出的擴(kuò)展目標(biāo)箱粒子PHD（Extended Target Box-Particle PHD，ET-Box-PHD）濾波器［9］和橢圓擴(kuò)展目標(biāo)箱粒子EET-BP-PHD（Ellipse Extended Target Box Particle PHD，）濾波器［10］是兩種新穎的擴(kuò)展目標(biāo)跟蹤算法，克服了ET-SMC-PHD 濾波過(guò)程中使用大量樣本粒子近似后驗(yàn)強(qiáng)度函數(shù)時(shí)存在計(jì)算量大的不足，提高了算法的運(yùn)行效率。

標(biāo)準(zhǔn)的ET-PHD 濾波器需要測(cè)量噪聲的先驗(yàn)信息，但這種測(cè)量噪聲的先驗(yàn)信息在許多實(shí)際應(yīng)用環(huán)境中都是未知的。測(cè)量噪聲協(xié)方差的顯著不匹配會(huì)導(dǎo)致ET-PHD 濾波器的目標(biāo)跟蹤性能退化，從而限制了ET-PHD 濾波器的實(shí)際應(yīng)用范圍。處理未知測(cè)量噪聲問(wèn)題的方法主要有交互多模型（Interactive Multiple Models，IMM）方法［11］、粒子方法［12］和變分貝葉斯（Variational Bayesian，VB）近似方法［13］等。由于具有較低的計(jì)算復(fù)雜度，VB 近似方法已經(jīng)成功地被用于基于RFS理論的濾波器中處理未知測(cè)量噪聲協(xié)方差環(huán)境中的多點(diǎn)目標(biāo)跟蹤問(wèn)題［14-17］，并且都取得了令人滿意的跟蹤效果。

為處理未知測(cè)量噪聲協(xié)方差情況下的多擴(kuò)展目標(biāo)跟蹤問(wèn)題，本文提出了一種基于ET-PHD 濾波器和VB 近似的擴(kuò)展目標(biāo)跟蹤算法。本文根據(jù)VB 近似理論及標(biāo)準(zhǔn)的ET-PHD 濾波器，提出了ET-PHD 濾波器的一種擴(kuò)展版本，并利用相對(duì)熵最小化方法給出了其在線性高斯假設(shè)條件下的解析實(shí)現(xiàn)。本文算法是一種遞歸的多擴(kuò)展目標(biāo)跟蹤算法，利用預(yù)測(cè)和更新過(guò)程實(shí)現(xiàn)后驗(yàn)強(qiáng)度函數(shù)和測(cè)量噪聲協(xié)方差的計(jì)算與估計(jì)，進(jìn)而獲取相應(yīng)的目標(biāo)跟蹤結(jié)果。

1 擴(kuò)展的ET-PHD濾波器

定義fk|k-1(xk|xk-1)和gk(zk|xk)分別為目標(biāo)狀態(tài)轉(zhuǎn)移函數(shù)和測(cè)量似然，在線性高斯假設(shè)條件下，ET-PHD 濾波器的系統(tǒng)模型可分別表示為下面形式：

其中：N(x；m，P)表示均值為m、協(xié)方差為P的高斯函數(shù)；xk和zk分別表示目標(biāo)狀態(tài)和測(cè)量值；Fk-1和Hk分別表示目標(biāo)狀態(tài)轉(zhuǎn)移矩陣和測(cè)量矩陣；Qk-1和Rk分別表示過(guò)程噪聲協(xié)方差矩陣和測(cè)量噪聲協(xié)方差矩陣。

根據(jù)ET-PHD 濾波器的系統(tǒng)模型，為描述基于VB 近似的ET-PHD 濾波器的公式體系，限定目標(biāo)狀態(tài)xk和測(cè)量噪聲協(xié)方差Rk是互相獨(dú)立的。若定義(xk，Rk)為k時(shí)間步擴(kuò)展目標(biāo)的增廣狀態(tài)變量，則增廣狀態(tài)變量(xk，Rk)的聯(lián)合轉(zhuǎn)移函數(shù)

其中，fk|k-1(xk|xk-1)和fk|k-1(Rk|Rk-1)分別表示單目標(biāo)狀態(tài)和測(cè)量噪聲協(xié)方差的轉(zhuǎn)移函數(shù)。

此外，增廣狀態(tài)的新生目標(biāo)RFS限定為泊松隨機(jī)集，定義Db，k(xk，Rk)為新生目標(biāo)RFS的強(qiáng)度函數(shù)，擴(kuò)展的ET-PHD 濾波器的遞歸公式描述如下。

假設(shè)Dk-1(xk-1，Rk-1)為k-1 時(shí)間步的目標(biāo)聯(lián)合后驗(yàn)強(qiáng)度函數(shù)，則k時(shí)間步預(yù)測(cè)的強(qiáng)度函數(shù)Dk|k-1(xk，Rk)可寫(xiě)為：

若在k時(shí)間步的測(cè)量集合為Zk，則聯(lián)合后驗(yàn)強(qiáng)度函數(shù)

其中：pD，k和κk(zk)分別表示傳感器的檢測(cè)概率和泊松雜波RFS 的強(qiáng)度；γ(xk)表示擴(kuò)展目標(biāo)xk產(chǎn)生的測(cè)量個(gè)數(shù)期望值；g(zk|xk，Rk)為測(cè)量zk的似然函數(shù)；p∠Zk表示劃分；|W|表示子集合W中元素個(gè)數(shù)；ωp和dW分別為劃分p和子集W的非負(fù)系數(shù)。ωp和dW計(jì)算式分別為：

其中：當(dāng)|W|=1時(shí)δ|W|，1=1，否則δ|W|，1=0。

在更新式（5）～（8）中，似然g(zk|xk，Rk)是目標(biāo)狀態(tài)xk和測(cè)量噪聲協(xié)方差變量Rk的耦合函數(shù)，聯(lián)合強(qiáng)度函數(shù)Dk(xk，Rk)是一種非解析形式。

2 基于VB近似的解析實(shí)現(xiàn)

擴(kuò)展的ET-PHD 濾波器解析實(shí)現(xiàn)的一次循環(huán)過(guò)程包含預(yù)測(cè)和更新兩個(gè)過(guò)程，解析實(shí)現(xiàn)的基礎(chǔ)除線性高斯模型之外，還需限定目標(biāo)的后驗(yàn)強(qiáng)度函數(shù)和新目標(biāo)強(qiáng)度函數(shù)為高斯與逆伽馬（Inverse-Gamma，IG）混合分布的形式。由此，定義k時(shí)間步的新目標(biāo)強(qiáng)度Db，k(xk，Rk)為下面的函數(shù)形式：

根據(jù)上述的假設(shè)與定義，擴(kuò)展的ET-PHD 濾波器解析實(shí)現(xiàn)的預(yù)測(cè)過(guò)程和更新過(guò)程詳述如下。

2.1 預(yù)測(cè)過(guò)程

若k-1 時(shí)間步的強(qiáng)度函數(shù)Dk-1(xk-1，Rk-1)可表示為下面的函數(shù)形式：

則k時(shí)間步預(yù)測(cè)的強(qiáng)度函數(shù)Dk|k-1(xk，Rk)可寫(xiě)作：

其中Db，k(xk，Rk)的表達(dá)式與式（9）相同，且有：

式中，ρ(i，j)表示取值范圍為(0，1]的衰減因子。

根據(jù)式（11）～（17），式（11）描述的預(yù)測(cè)強(qiáng)度函數(shù)Dk|k-1(xk，Rk)可以表示為下面的函數(shù)形式：

利用式（18）給出的預(yù)測(cè)的聯(lián)合強(qiáng)度函數(shù)Dk|k-1(xk，Rk)的表達(dá)式，預(yù)測(cè)的測(cè)量噪聲協(xié)方差的計(jì)算公式為：

其中diag(?)表示構(gòu)建對(duì)角矩陣的運(yùn)算。

2.2 更新過(guò)程

在更新過(guò)程，k時(shí)間步的測(cè)量集利用距離劃分方法［8］進(jìn)行劃分，進(jìn)而實(shí)現(xiàn)后驗(yàn)強(qiáng)度函數(shù)Dk(xk，Rk)的計(jì)算，將式（6）代入式（5）可得：

若k時(shí)間步的預(yù)測(cè)的聯(lián)合強(qiáng)度函數(shù)Dk|k-1(xk，Rk)的函數(shù)形式與式（18）相同，將式（18）代入式（22）可得：

其中：為方便描述強(qiáng)度函數(shù)(xk，Rk)的函數(shù)表達(dá)式，定義如下符號(hào)：

其中⊕表示垂直向量連接。

利用文獻(xiàn)［13］給出的相對(duì)熵最小化結(jié)論，聯(lián)合強(qiáng)度函數(shù)，Rk)可表示為：

為獲取最優(yōu)的計(jì)算結(jié)果，聯(lián)合強(qiáng)度函數(shù)(xk，Rk)表達(dá)式中的各參數(shù)采用定點(diǎn)迭代算法進(jìn)行計(jì)算，定點(diǎn)迭代算法描述如下：

3）利用下面公式進(jìn)行若干次的迭代計(jì)算：

其中blkdiag(?)表示構(gòu)建分塊對(duì)角矩陣。

定點(diǎn)迭代算法的迭代次數(shù)可根據(jù)實(shí)際需要進(jìn)行設(shè)置，此外式（32）中的其他參數(shù)計(jì)算如下：

計(jì)算得到的強(qiáng)度函數(shù)Dk(xk，Rk)用于獲取擴(kuò)展目標(biāo)的狀態(tài)和個(gè)數(shù)，并通過(guò)修剪和合并限制Dk(xk，Rk)中高斯與逆伽馬混合分量的個(gè)數(shù)，進(jìn)而提高運(yùn)算效率，實(shí)現(xiàn)未知噪聲環(huán)境中多個(gè)擴(kuò)展目標(biāo)的實(shí)時(shí)有效的跟蹤。

3 仿真與結(jié)果分析

利用兩個(gè)仿真實(shí)驗(yàn)驗(yàn)證本文算法的有效性，實(shí)驗(yàn)場(chǎng)景設(shè)置為包含雜波且目標(biāo)個(gè)數(shù)隨機(jī)變化的一個(gè)二維平面。實(shí)驗(yàn)環(huán)境為：英特爾八代酷睿i5-8250 處理器、8 GB 運(yùn)行內(nèi)存、Windows 10操作系統(tǒng)和R2012a版本的Matlab仿真軟件。

k時(shí)間步的目標(biāo)狀態(tài)向量xk由目標(biāo)位置和速度構(gòu)成，測(cè)量值z(mì)k定義為受噪聲影響的目標(biāo)位置信息。仿真實(shí)驗(yàn)中采用式（1）和式（2）中描述的目標(biāo)跟蹤模型，模型中的各個(gè)參數(shù)定義如下：

抽樣間隔定義為Δ=1 s，目標(biāo)的存活和檢測(cè)概率分別設(shè)置為pS，k=0.99 和pD，k=0.99，每個(gè)目標(biāo)產(chǎn)生的測(cè)量個(gè)數(shù)的泊松率設(shè)置為γ=10。標(biāo)準(zhǔn)差σq已知，值設(shè)置為σq=1.95 m；跟蹤過(guò)程中標(biāo)準(zhǔn)差σr的值是未知的。仿真實(shí)驗(yàn)中用到的新目標(biāo)強(qiáng)度函數(shù)定義為：

強(qiáng)度函數(shù)中高斯分量的修剪閾值為Tw=10-3，合并閾值為U=4，高斯分量個(gè)數(shù)最大值設(shè)置為Jmax=100。逆伽馬分布的參數(shù)初始值設(shè)置為α0=β0=1，退化因子的值定義為ρ=0.9。

不同算法的目標(biāo)位置估計(jì)精度利用最優(yōu)子模式分配（Optimal Sub-Pattern Assignment，OSPA）距離［18］進(jìn)行評(píng)估，計(jì)算公式為：

其中：X和Y分別表示目標(biāo)的真實(shí)狀態(tài)集和狀態(tài)估計(jì)集；p和c分別設(shè)置為p=2和c=3 000 m。

仿真場(chǎng)景中設(shè)定有三個(gè)目標(biāo)運(yùn)動(dòng)，圖1 為三個(gè)目標(biāo)的真實(shí)運(yùn)動(dòng)軌跡。圖1所示的仿真場(chǎng)景中，目標(biāo)1的運(yùn)動(dòng)時(shí)間設(shè)置為50個(gè)持續(xù)時(shí)間步，目標(biāo)2的運(yùn)動(dòng)時(shí)間在5到45時(shí)間步之間，目標(biāo)3在時(shí)間步10出現(xiàn)。

圖1 真實(shí)目標(biāo)軌跡Fig.1 True target tracks

3.1 固定參數(shù)實(shí)驗(yàn)

本實(shí)驗(yàn)中雜波設(shè)置為速率為r=10泊松隨機(jī)有限集，ETGM-PHD 濾波器中的標(biāo)準(zhǔn)差設(shè)置為σr=0.65。圖2 給出的是一次實(shí)驗(yàn)運(yùn)行中50 個(gè)時(shí)間步的真實(shí)目標(biāo)軌跡和包含雜波量測(cè)的傳感器測(cè)量值。

圖2 真實(shí)目標(biāo)軌跡與測(cè)量值Fig.2 True target tracks and measurements

圖3為不同算法在一次實(shí)驗(yàn)運(yùn)行中目標(biāo)的PHD權(quán)值。從圖3 中可以看出，在17 和22 等時(shí)間步上，ET-GM-PHD 濾波器估計(jì)的不同目標(biāo)的PHD 權(quán)值存在不同程度的減小，但本文算法的PHD 權(quán)值估計(jì)卻十分可靠。原因是ET-GM-PHD 濾波器使用了不匹配的測(cè)量噪聲參數(shù)進(jìn)行濾波計(jì)算，導(dǎo)致了濾波性能的下降；但本文算法在跟蹤過(guò)程能夠獲得合適的測(cè)量噪聲參數(shù)，進(jìn)而取得了可靠的目標(biāo)PHD權(quán)值。

圖3 本文算法和ET-GM-PHD濾波器中不同目標(biāo)的PHD權(quán)值Fig.3 PHD weights of different targets of proposed algorithm and ET-GM-PHD filter

圖4給出的是一次實(shí)驗(yàn)運(yùn)行中不同算法在50個(gè)時(shí)間步上的目標(biāo)位置估計(jì)。從圖4 中可以看出，在17、22 和37 時(shí)間步上，ET-GM-PHD濾波器未能獲取相應(yīng)的目標(biāo)位置，原因是ETGM-PHD濾波器中目標(biāo)PHD權(quán)值的大幅下降造成了這種目標(biāo)位置估計(jì)丟失現(xiàn)象的產(chǎn)生，但本文算法在跟蹤過(guò)程中沒(méi)有出現(xiàn)目標(biāo)丟失的現(xiàn)象。

圖4 真實(shí)目標(biāo)軌跡與位置估計(jì)Fig.4 True target tracks and position estimations

在實(shí)驗(yàn)中，不同算法的平均性能從平均的目標(biāo)個(gè)數(shù)估計(jì)和平均的OSPA 距離這兩個(gè)方面進(jìn)行分析，相應(yīng)的仿真結(jié)果分別由圖5 和圖6 給出。從圖5 中可以看出，在11～50 時(shí)間步上，ET-GM-PHD 濾波器平均的目標(biāo)個(gè)數(shù)估計(jì)與目標(biāo)個(gè)數(shù)的真實(shí)值有著明顯的偏差，該偏差產(chǎn)生的原因是ET-GM-PHD濾波器中的不合適的目標(biāo)PHD 權(quán)值估計(jì)。通過(guò)對(duì)比發(fā)現(xiàn)，本文算法平均的目標(biāo)個(gè)數(shù)估計(jì)非常接近真實(shí)的目標(biāo)個(gè)數(shù)，偏差很小。圖6 顯示，本文算法的平均OSPA 距離遠(yuǎn)小于ET-GM-PHD 濾波器的平均OSPA 距離，這種情況在11～50 時(shí)間步之間尤為明顯。

圖5 不同算法平均的目標(biāo)個(gè)數(shù)估計(jì)Fig.5 Average target number estimations of different algorithms

圖6 不同算法的平均OSPA距離Fig.6 Average OSPA distances of different algorithms

3.2 變化參數(shù)實(shí)驗(yàn)

本仿真實(shí)驗(yàn)中，不同算法的仿真結(jié)果為不同參數(shù)取值情況下進(jìn)行的1 000 次仿真實(shí)驗(yàn)中獲取的平均值，ET-GM-PHD濾波器中標(biāo)準(zhǔn)差σr取值分別為0.65、0.75、0.95 和1.65 四種情況。首先，目標(biāo)量測(cè)的泊松率固定為γ=10，雜波速率在5到50 之間以等變化間隔5 進(jìn)行取值，不同算法在不同的雜波速率取值情況下時(shí)間上平均的OSPA 距離的對(duì)比結(jié)果如圖7所示。然后，雜波速率固定為r=10，目標(biāo)量測(cè)的泊松率在3到12 之間變化，不同目標(biāo)量測(cè)泊松率時(shí)的時(shí)間平均的OSPA距離由圖8給出。

圖7 顯示，不同算法的跟蹤精度隨著雜波速率的增大而下降，但本文算法時(shí)間平均的OSPA 距離在不同的雜波速率上總是一直小于ET-GM-PHD 濾波器相應(yīng)的OSPA 距離。圖8顯示，當(dāng)目標(biāo)量測(cè)的泊松率在10 附近取值時(shí)，所有算法都可以獲取最佳的跟蹤精度；但對(duì)于所有的目標(biāo)量測(cè)泊松率，本文算法都可以獲得精度更好的目標(biāo)位置估計(jì)。根據(jù)圖7 和圖8給出的仿真對(duì)比結(jié)果可以得到，在未知測(cè)量噪聲協(xié)方差環(huán)境中，本文算法能夠獲得可信賴的多擴(kuò)展目標(biāo)跟蹤結(jié)果。

圖7 不同雜波速率下的時(shí)間平均的OSPA距離Fig.7 Time averaged OSPA distances with different clutter rates

圖8 不同目標(biāo)量測(cè)泊松率下的時(shí)間平均的OSPA距離Fig.8 Time averaged OSPA distances with different Poisson rates of target measurement

4 結(jié)語(yǔ)

為處理未知測(cè)量噪聲協(xié)方差情況下的多擴(kuò)展目標(biāo)跟蹤問(wèn)題，利用VB近似理論和高斯與逆伽馬混合分布表示的聯(lián)合強(qiáng)度函數(shù)，本文提出了一種擴(kuò)展的ET-PHD 濾波器及其線性高斯解析實(shí)現(xiàn)。仿真結(jié)果表明：本文算法能夠以可靠的性能跟蹤未知測(cè)量噪聲協(xié)方差環(huán)境中個(gè)數(shù)未知且隨時(shí)間變化的多個(gè)擴(kuò)展目標(biāo)。在本文研究的基礎(chǔ)之上，下一步將研究ET-SMCPHD 濾波器的算法優(yōu)化方法，以提高ET-SMC-PHD 濾波器的算法運(yùn)行速度。