基于ET-PHD濾波器和變分貝葉斯近似的擴展目標跟蹤算法

何祥宇，李 靜，楊數強，夏玉杰

（1.洛陽師范學院物理與電子信息學院，河南洛陽 471934；2.洛陽師范學院信息技術學院，河南洛陽 471934）

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

0 引言

目標跟蹤技術是一個備受國內外學者關注的研究課題，應用范圍十分廣泛［1-3］。多目標跟蹤（Multi-Target Tracking，MTT）技術研究是目標跟蹤技術領域內的一個重要研究分支，多目標跟蹤的目的是利用受雜波和噪聲影響的傳感器測量值估計目標的個數及狀態。隨機有限集（Random Finite Set，RFS）理論被認為是處理多目標跟蹤問題的一種簡潔高效的方法，通過把多目標狀態集和測量集定義為隨機有限集，該理論建立了嚴謹的最優貝葉斯多目標濾波框架。

在基于RFS 理論的多目標跟蹤方法中，概率假設密度（Probability Hypothesis Density，PHD）濾波器［4］是一種很有潛力的多目標貝葉斯濾波算法，它不用復雜的數據關聯過程就可以處理復雜環境中的多目標跟蹤問題。標準的PHD 濾波器是基于點目標假設的多目標跟蹤算法，點目標假設即：在每一個時間步，每一個目標假設只能產生至多一個量測值。與點目標不同的是，一個擴展目標在每一個時間步上可產生多個量測值［5］。因此，應用標準的PHD 濾波器跟蹤多個擴展目標會產生不可靠的跟蹤結果。

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

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

為處理未知測量噪聲協方差情況下的多擴展目標跟蹤問題，本文提出了一種基于ET-PHD 濾波器和VB 近似的擴展目標跟蹤算法。本文根據VB 近似理論及標準的ET-PHD 濾波器，提出了ET-PHD 濾波器的一種擴展版本，并利用相對熵最小化方法給出了其在線性高斯假設條件下的解析實現。本文算法是一種遞歸的多擴展目標跟蹤算法，利用預測和更新過程實現后驗強度函數和測量噪聲協方差的計算與估計，進而獲取相應的目標跟蹤結果。

1 擴展的ET-PHD濾波器

定義fk|k-1(xk|xk-1)和gk(zk|xk)分別為目標狀態轉移函數和測量似然，在線性高斯假設條件下，ET-PHD 濾波器的系統模型可分別表示為下面形式：

其中：N(x；m，P)表示均值為m、協方差為P的高斯函數；xk和zk分別表示目標狀態和測量值；Fk-1和Hk分別表示目標狀態轉移矩陣和測量矩陣；Qk-1和Rk分別表示過程噪聲協方差矩陣和測量噪聲協方差矩陣。

根據ET-PHD 濾波器的系統模型，為描述基于VB 近似的ET-PHD 濾波器的公式體系，限定目標狀態xk和測量噪聲協方差Rk是互相獨立的。若定義(xk，Rk)為k時間步擴展目標的增廣狀態變量，則增廣狀態變量(xk，Rk)的聯合轉移函數

其中，fk|k-1(xk|xk-1)和fk|k-1(Rk|Rk-1)分別表示單目標狀態和測量噪聲協方差的轉移函數。

此外，增廣狀態的新生目標RFS限定為泊松隨機集，定義Db，k(xk，Rk)為新生目標RFS的強度函數，擴展的ET-PHD 濾波器的遞歸公式描述如下。

假設Dk-1(xk-1，Rk-1)為k-1 時間步的目標聯合后驗強度函數，則k時間步預測的強度函數Dk|k-1(xk，Rk)可寫為：

若在k時間步的測量集合為Zk，則聯合后驗強度函數

其中：pD，k和κk(zk)分別表示傳感器的檢測概率和泊松雜波RFS 的強度；γ(xk)表示擴展目標xk產生的測量個數期望值；g(zk|xk，Rk)為測量zk的似然函數；p∠Zk表示劃分；|W|表示子集合W中元素個數；ωp和dW分別為劃分p和子集W的非負系數。ωp和dW計算式分別為：

其中：當|W|=1時δ|W|，1=1，否則δ|W|，1=0。

在更新式（5）～（8）中，似然g(zk|xk，Rk)是目標狀態xk和測量噪聲協方差變量Rk的耦合函數，聯合強度函數Dk(xk，Rk)是一種非解析形式。

2 基于VB近似的解析實現

擴展的ET-PHD 濾波器解析實現的一次循環過程包含預測和更新兩個過程，解析實現的基礎除線性高斯模型之外，還需限定目標的后驗強度函數和新目標強度函數為高斯與逆伽馬（Inverse-Gamma，IG）混合分布的形式。由此，定義k時間步的新目標強度Db，k(xk，Rk)為下面的函數形式：

根據上述的假設與定義，擴展的ET-PHD 濾波器解析實現的預測過程和更新過程詳述如下。

2.1 預測過程

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

則k時間步預測的強度函數Dk|k-1(xk，Rk)可寫作：

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

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

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

利用式（18）給出的預測的聯合強度函數Dk|k-1(xk，Rk)的表達式，預測的測量噪聲協方差的計算公式為：

其中diag(?)表示構建對角矩陣的運算。

2.2 更新過程

在更新過程，k時間步的測量集利用距離劃分方法［8］進行劃分，進而實現后驗強度函數Dk(xk，Rk)的計算，將式（6）代入式（5）可得：

若k時間步的預測的聯合強度函數Dk|k-1(xk，Rk)的函數形式與式（18）相同，將式（18）代入式（22）可得：

其中：為方便描述強度函數(xk，Rk)的函數表達式，定義如下符號：

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

利用文獻［13］給出的相對熵最小化結論，聯合強度函數，Rk)可表示為：

為獲取最優的計算結果，聯合強度函數(xk，Rk)表達式中的各參數采用定點迭代算法進行計算，定點迭代算法描述如下：

3）利用下面公式進行若干次的迭代計算：

其中blkdiag(?)表示構建分塊對角矩陣。

定點迭代算法的迭代次數可根據實際需要進行設置，此外式（32）中的其他參數計算如下：

計算得到的強度函數Dk(xk，Rk)用于獲取擴展目標的狀態和個數，并通過修剪和合并限制Dk(xk，Rk)中高斯與逆伽馬混合分量的個數，進而提高運算效率，實現未知噪聲環境中多個擴展目標的實時有效的跟蹤。

3 仿真與結果分析

利用兩個仿真實驗驗證本文算法的有效性，實驗場景設置為包含雜波且目標個數隨機變化的一個二維平面。實驗環境為：英特爾八代酷睿i5-8250 處理器、8 GB 運行內存、Windows 10操作系統和R2012a版本的Matlab仿真軟件。

k時間步的目標狀態向量xk由目標位置和速度構成，測量值zk定義為受噪聲影響的目標位置信息。仿真實驗中采用式（1）和式（2）中描述的目標跟蹤模型，模型中的各個參數定義如下：

抽樣間隔定義為Δ=1 s，目標的存活和檢測概率分別設置為pS，k=0.99 和pD，k=0.99，每個目標產生的測量個數的泊松率設置為γ=10。標準差σq已知，值設置為σq=1.95 m；跟蹤過程中標準差σr的值是未知的。仿真實驗中用到的新目標強度函數定義為：

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

不同算法的目標位置估計精度利用最優子模式分配（Optimal Sub-Pattern Assignment，OSPA）距離［18］進行評估，計算公式為：

其中：X和Y分別表示目標的真實狀態集和狀態估計集；p和c分別設置為p=2和c=3 000 m。

仿真場景中設定有三個目標運動，圖1 為三個目標的真實運動軌跡。圖1所示的仿真場景中，目標1的運動時間設置為50個持續時間步，目標2的運動時間在5到45時間步之間，目標3在時間步10出現。

圖1 真實目標軌跡Fig.1 True target tracks

3.1 固定參數實驗

本實驗中雜波設置為速率為r=10泊松隨機有限集，ETGM-PHD 濾波器中的標準差設置為σr=0.65。圖2 給出的是一次實驗運行中50 個時間步的真實目標軌跡和包含雜波量測的傳感器測量值。

圖2 真實目標軌跡與測量值Fig.2 True target tracks and measurements

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

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

圖4給出的是一次實驗運行中不同算法在50個時間步上的目標位置估計。從圖4 中可以看出，在17、22 和37 時間步上，ET-GM-PHD濾波器未能獲取相應的目標位置，原因是ETGM-PHD濾波器中目標PHD權值的大幅下降造成了這種目標位置估計丟失現象的產生，但本文算法在跟蹤過程中沒有出現目標丟失的現象。

圖4 真實目標軌跡與位置估計Fig.4 True target tracks and position estimations

在實驗中，不同算法的平均性能從平均的目標個數估計和平均的OSPA 距離這兩個方面進行分析，相應的仿真結果分別由圖5 和圖6 給出。從圖5 中可以看出，在11～50 時間步上，ET-GM-PHD 濾波器平均的目標個數估計與目標個數的真實值有著明顯的偏差，該偏差產生的原因是ET-GM-PHD濾波器中的不合適的目標PHD 權值估計。通過對比發現，本文算法平均的目標個數估計非常接近真實的目標個數，偏差很小。圖6 顯示，本文算法的平均OSPA 距離遠小于ET-GM-PHD 濾波器的平均OSPA 距離，這種情況在11～50 時間步之間尤為明顯。

圖5 不同算法平均的目標個數估計Fig.5 Average target number estimations of different algorithms

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

3.2 變化參數實驗

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

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

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

圖8 不同目標量測泊松率下的時間平均的OSPA距離Fig.8 Time averaged OSPA distances with different Poisson rates of target measurement

4 結語

為處理未知測量噪聲協方差情況下的多擴展目標跟蹤問題，利用VB近似理論和高斯與逆伽馬混合分布表示的聯合強度函數，本文提出了一種擴展的ET-PHD 濾波器及其線性高斯解析實現。仿真結果表明：本文算法能夠以可靠的性能跟蹤未知測量噪聲協方差環境中個數未知且隨時間變化的多個擴展目標。在本文研究的基礎之上，下一步將研究ET-SMCPHD 濾波器的算法優化方法，以提高ET-SMC-PHD 濾波器的算法運行速度。