一種改進的直方圖概率多假設多目標跟蹤方法

張奕群,尹立凡,王碩,孫承鋼

北京電子工程總體研究所,北京 100854

未來強對抗的戰場環境對目標跟蹤系統的多目標處理能力提出了更高的要求。而傳統的檢測前跟蹤(Track Before Detect, TBD)方法[1-3]在多目標情況下會遇到“維數災難”問題[4]。為此人們致力于發展高效的多目標處理方法,直方圖概率多假設跟蹤[5](Histogram Probabilistic Multi-Hypothesis Tracking, H-PMHT)就是其中一種。該方法由Streit提出,其計算量與目標數呈線性增長,因而在多目標情況下仍可以保持較高的計算效率[5]。

H-PMHT方法有這一優勢的關鍵在其多項式分布(Multinomial Distribution, MD)量測模型[6-7]。Streit等使用期望極大化(Expectation Maximization,EM)方法[8]在經驗貝葉斯意義下對該模型未知分布參數進行了估計,其結果是包括目標狀態和強度參數在內的未知參數尋優過程可以被解耦為多個易于計算的子尋優過程,由此為H-PMHT方法計算量與目標數呈線性關系奠定了基礎。

在H-PMHT方法基礎上,Davey等以多維泊松分布(Poisson Distribution, PD)替換了多項式分布[9],從而得到了不同目標強度參數互相解耦的改進的量測模型。因而,以此得到的泊松H-PMHT(Poisson H-PMHT,P-HPMHT)方法可以在估計中利用強度參數先驗信息并同時保持計算量與目標數間的線性關系。因為P-HPMHT方法在目標數目較多時仍可以有效的利用先驗信息,其表現比H-PMHT方法更為穩定[10]。

H-PMHT方法的一個主要問題是在低信噪比條件下跟蹤能力較差。Davey等在文獻[1,11]中通過仿真對包括H-PMHT方法在內的幾種TBD方法的跟蹤檢測能力進行了比對研究,研究表明H-PMHT方法的低信噪比跟蹤檢測能力較弱。究其原因,H-PMHT方法的量測方程中并未考慮傳感器噪聲[12]。因此,其跟蹤檢測能力在低信噪比條件下會明顯降低。自2000年提出至今H-PMHT方法對這一問題仍未能解決。其難點在于如果在多項式分布量測模型中引入噪聲模型,則量測量的分布難以處理,進而無法得到有效的跟蹤方法。P-HPMHT方法也存在相同的問題。目前,這一問題仍然沒有解決[12-13]。

在H-PMHT方法中使用MD模型為量測模型,其中隱含著一個前提假設,即存在一個量化單位使得傳感器數據量化結果滿足多項式分布。然而,在實際中該量化單位無法確切得到,或并不存在。因此,若量化單位取值偏差較大,將導致估計結果誤差增加。為解決這一問題,Streit提出了一種“修正先驗信息”(Modified Prior)方法[5],以補償量化單位取值偏差對估計精度的影響。隨后,Davey等也使用了類似的方法解決了P-HPMHT方法存在的相同問題[9]。然而,這種解決方法是一種工程中的變通方法,其缺乏理論依據[14]。Willett等在文獻[15]中以目標的狀態估值為條件對量化單位進行了尋優。但該方法目標的狀態估計受到量化單位估值誤差影響,導致其估計仍存在偏差。

本文提出了一種帶傳感器噪聲模型的H-PMHT方法：通過選擇恰當的量測模型丟失數據,實現了引入噪聲模型條件下尋優過程的解耦。又通過提出一種類泊松概率質量函數多重卷積簡化計算方法,實現了尋優過程的高效計算。最終,在引入噪聲模型條件下實現了計算量隨著目標數線性增長。且同時,這一方法還給出了一種理論化的量化單位尋優辦法。

全文內容安排如下：第1節,給出了使用到的符號與定義,并引出了要解決的問題；第2節,提出了帶噪聲模型的H-PMHT方法；第3節,在單目標低信噪比、雙目標交叉及平行運動的場景下做了仿真實驗,對比了H-PMHT方法、P-HPMHT方法、本文提出的跟蹤方法和經典的聯合概率數據關聯(Joint Probabilistic Data Association,JPDA)方法[16]；第4節,總結了主要結論,并指出了后續還需要關注的問題。

1 符號、定義與問題的引出

(1)

(2)

(3)

為第k時刻強度參數的集合。

假設傳感器空間B分為I個不相交的傳感器分辨單元。記其中第i號分辨單元子空間為Bi?B,第k時刻其上分布的來自目標及雜波的能量一般表示為[1]

(4)

式中：i=1,2,…,I。

(5)

(6)

與之對應的傳感器輸出為

(7)

(8)

H-PMHT方法和P-HPMHT方法的量測模型方程均為

(9)

2 帶傳感器噪聲模型的H-PMHT方法

(10)

(11)

(12)

(13)

令量化輸出

(14)

(15)

式中：*為卷積運算,f(x)|z為函數f在z處取值。

(16)

將式(16)代入式(11),即可以得到傳感器輸出Zk的改進概率模型p(Zk;Xk,Λk,η)。

令參數集

Θk={Xk,Λk,η}

(17)

以下,要解決的主要問題是通過求解式

(18)

來得到Θk的極大后驗估計。式中：p(Zk|Θk)為式(11)；pΘ(Θk)為Θk的先驗分布。

2.1 EM方法

在H-PMHT方法、P-HPMHT方法等一類估計問題中,都通過使用EM方法[7-8]對其模型參數進行尋優。本文也采用相同的方法求解式(18)。

根據EM方法,式(18)的優化問題可以變換為式(19)定義的一系列尋優問題

(19)

(20)

式中：Ok為丟失數據(Missing Data)的集合。

式(20)中,丟失數據Ok的選擇非常重要。只有選擇恰當的丟失數據Ok,才可以使得Q能夠解耦,從而易于計算。

由式(8)可知,計數值

(21)

其中,

(22)

式中：

(23)

(24)

記來自第m號密度分量并落在Bi上的“樣本點”為

而后,分別定義集合

(25)

i=1,2,…,I}

(26)

這里丟失數據為

(27)

2.2 期望步驟

期望步驟的目標是計算輔助函數Q,其定義見式(20)。根據貝葉斯公式,得到該函數表達式中的對數項為

lgp(Ζk,Ok,Θk)=lgp(Ζk,Ok|Θk)+lgp(Θk)

(28)

式中：等號右側第1項為Ζk和Ok的對數聯合條件概率密度函數；第2項為Θk的對數先驗分布。

(29)

(30)

利用以上結果及概率乘法公式,即可以得到丟失數據Ok的條件概率密度函數為

(31)

(32)

將式(32)與式(31)相乘,利用條件概率乘法公式得到Ζk和Ok的聯合條件概率密度函數為

(33)

(34)

對式(34)等號兩邊求對數,可以得到

(35)

式中：C為與Θk無關的量。

假設各目標狀態間相互獨立,Xk的先驗密度函數為

(36)

式中：pX為目標狀態的先驗密度函數。

假設Λk的先驗密度函數為

(37)

式中：pΛ為強度參數的先驗密度函數。

又假設Xk和Λk相互獨立,得到Θk的對數先驗密度函數為

(38)

將式(35)和式(38)一并代入式(28),對其等號兩側求條件期望E[·|·],而后利用約束條件

(39)

進行整理,輔助函數Q為

(40)

在式(40)中,等號右側第1項僅與η相,第2～4項僅與Λk相關,第5～6項僅與Xk相關。易見輔助函數Q的極值求解問題被拆分成多個獨立的極值求解問題。

根據貝葉斯公式知,丟失數據Ok的條件密度函數為

(41)

其中分母和分子項依次見式(11)及式(34)。利用該結果將式(40)做化簡,其等號右側第3行中的條件期望可以整理為式(33)。

2.3 極大化步驟

接下來,計算式(40)的極值。首先,將式(40)整理為

(42)

式中：各分量為子輔助函數,其定義分別為

(43)

(44)

(45)

假設目標狀態的演化規律可以用一階馬爾科夫模型來描述,則由C-K方程可得

(46)

(47)

(48)

式中：

(49)

(50)

(51)

(52)

因而,借鑒Granstr?m等提出的一種泊松強度參數濾波方法[19]對式(52)尋優。假設

(53)

(54)

式中：ξ為遺忘因子；Δ為不同采樣時刻時間間隔。可以得到其更新方程為

(55)

(56)

2.3.2 求Qη的極值

假設傳感器噪聲為高斯噪聲,即

(57)

式中：μ為噪聲均值;σ2為噪聲方差。

對Qη表達式(45)中的對數項進行展開,得

(58)

計算式(58)展開結果中各項條件期望E[·|·],得到Qη展開為

(59)

對式(59)求導并令其為0,進而得到η的估計值為

(60)

式(60)計算過程中涉及到的2個期望值的計算問題,將在2.4節中討論。另外,以上為高斯假設下得到的估計方法,對于非高斯情況,Qη的極大化問題有待進一步研究。

引理定義函數族

(61)

(62)

式中：cq為系數,其具體值由λa、λb、qa及qb計算得到。

(63)

因為2個函數卷積的特征函數等于其各自特征函數的乘積,所以證明等式(62)成立等價于證明存在系數cq使得

(64)

(65)

(66)

不失一般性,假設0≤qa≤qb。利用以上特征函數表達式,通過合并同類項可得

Pqa(τ;λa)×Pqb(τ;λb)=e(λa+λb)(ejτ-1)·

(67)

式中：系數

(68)

由式(65)～式(68)結果知,式(64)成立,又等價于存在系數cq使得

(69)

當m>n時,S(n,m)=0,將式(66)代入式(69) 等號左側,即可得到

(70)

對照式(70)與式(67)知,若cq存在,即向量

(71)

存在,矩陣方程

(72)

可解,式中：向量

(73)

矩陣U為(qa+qb+1)×(qa+qb+1)的上三角陣,且其中第i行、第j列元素為

(74)

因為,S(j-1,i-1)=1、λa和λb都大于0,在矩陣U對角線上各項非0,所以其行列式非0。因此,存在逆矩陣U-1。也因此,向量c存在且

(75)

因此,引理證明成立。

(76)

式(76)中第2個等號成立是因為式(13)成立。

(77)

式中：系數c1為

(78)

它由式(75)計算得到。

當m1=m2=m時,其條件期望為

(79)

當m1≠m2時,其條件期望值為

(80)

則使用前述引理可將上述2種情況都化簡為

(81)

式中：當m1=m2=m時,系數cq依次為

(82)

當m1≠m2時,為

(83)

進一步,利用以上結果還易得

(84)

因而,還可以得到更簡化的計算式

(85)

總結以上過程,對Θk進行估計,其步驟如下：

步驟1初始設置

設定量化單位初始值η(0)為上一時刻估計值。

設定迭代次數計數值d=0。

步驟2迭代更新

更新迭代次數計數值d=d+1。

步驟3回到步驟2,直到估計結果收斂。

3 仿真驗證

在仿真中,假設傳感器空間B為xOy二維坐標平面上的正方形區域,其右上頂點為坐標原點O且兩直角邊分別在x軸及y軸上。設定傳感器共包含20×20個分辨單元。以(i,j)表示沿x軸向第i行,沿y軸向第j列的分辨單元,設定其空間為正方形區域[(i-1)Δl,(j-1)Δl]×[iΔl,jΔl],其中Δl為每個分辨單元的邊長。設定Δl=0.5。

(86)

式中：F為狀態轉移矩陣

(87)

ωk-1為高斯過程噪聲,其均值為0,協方差陣為

(88)

式中：qω表征噪聲強度,設定qω=0.001。

假設第m號目標的能量分布為高斯密度函數

(89)

記式(84)在第(i,j)號分辨單元上的積分為

(90)

仿真中,模擬第(i,j)號分辨單元的輸出為

(91)

進而,定義目標的信噪比為

i,j=1,2,…,I

(92)

根據以上給出的參數,在目標信噪比為10 dB 和6 dB場合使用到的2幀傳感器數據的示意圖見圖1。其中,目標位于左下角圓圈中心,目標強度擴散3×3個分辨單元。在10 dB時,目標還可能被分辨出來,但6 dB時,已經不能分辨出其位置。

圖1 仿真數據示意

3.1 仿真場景

場景1單個目標做勻速直線運動。目標的起始狀態為x0=Δl[16.5,-1,5.5,1]T,即目標初始位置在分辨單元(17,6)的中心位置且其運動速度沿x軸及y軸分別為-1和1分辨單元/幀。目標的真實軌跡序列為X={xk:k=1,2,…,K},其中K設定為11幀。設定目標信噪比SNR分別為10、6和3 dB,根據式(91)模擬生成1 000組×11幀的傳感器數據。再以無目標情況模擬生成1 000幀僅有噪聲的傳感器數據。

3.2 跟蹤實現及性能評價

在場景1中,以H-PMHT方法、P-HPMHT方法及本文提出的H-PMHT/SN方法跟蹤目標。其中,第k時刻本文方法的工作分為軌跡起始和軌跡維持2個部分。

1) 軌跡起始：設定檢測概率為0.9,在給定的信噪比條件下確定檢測門限。分割第k幀數據以得到該時刻的目標,并與前一幀目標關聯構建新的軌跡。其條件為相應的目標間的距離小于1.5Δl。

利用M/N準則[6]管理軌跡,其中M為2,N為3,即連續3幀軌跡質量的測度大于門限,則維持這一軌跡,否則將其刪除。其中,軌跡質量,即信噪比估計值按下式計算[6]

(93)

根據仿真結果,從跟蹤檢測能力和估計精度2個方面對跟蹤方法進行評價。其中,以特征曲線(Receiver Operating Characteristic,ROC)評價跟蹤檢測能力；以位置估值的均方根誤差(Root Mean Square Error,RMSE)評價各方法的估計精度。ROC定義為檢測比率相對誤跟蹤比率的函數曲線。首先使用某一跟蹤方法處理1 000幀僅有噪聲的傳感器數據,則檢測到“目標”的一幀稱為“誤跟蹤”幀。以“誤跟蹤”幀數比總的幀數為該跟蹤方法的誤跟蹤比率。其中“誤跟蹤”幀的判據是：某一幀和與之相連的前2幀中估計得到的3幀軌跡,其中有任2幀的目標位置估值與目標真實位置間的距離誤差大于2σsp。而后使用給定的跟蹤方法處理不同給定信噪比條件下生成的1 000 組×11幀數據,其中正確檢測到目標的總的幀數比總的幀數,即為檢測比率。由此,得到的ROC曲線為圖2。

圖2 ROC曲線

從結果看在低信噪比條件下,本文跟蹤方法較H-PMHT方法和P-HPMHT方法表現更好,尤其是在信噪比為6 dB及3 dB時。以誤跟蹤比率為10-3時結果進行分析,在信噪比為6 dB時,本文方法檢測比率提升了20%,信噪比為3 dB時,提升了10%。另外,各方法RMSE見圖3。可以發現當跟蹤上目標時,各方法的精度相近。

圖3 RMSE圖示

設置場景1的目的是對H-PMHT/SN方法的低信噪比跟蹤能力進行評價。場景2則是檢驗

該方法在2個目標交叉運動時的跟蹤能力。

以下采用軌跡丟失率(Track Loss Ratio,TLR)對跟蹤結果進行評價。TLR定義為“誤跟蹤”的軌跡數比上總的蒙特卡洛仿真數。其中,“誤跟蹤”與仿真1中給出的定義相同,仿真數為1 000次。

在10 dB和6 dB信噪比條件下,比對跟蹤方法的軌跡丟失率,見圖4。當信噪比為10 dB時,H-PMHT、P-HPMHT和H-PMHT/SN這3種方法總體上保持了較低的軌跡丟失率,且并沒有受到軌跡交叉的較大影響。但JPDA方法在軌跡交叉時,其軌跡丟失率明顯升高。出現這一結果的原因是H-PMHT、P-HPMHT和H-PMHT/SN這3種方法是TBD算法,它們更充分地考慮了目標接近時的相互影響,也較少受到軌跡交叉影響。

當信噪比為6 dB時,H-PMHT和P-HPMHT方法的軌跡丟失率都出現了明顯的上升。這是因為H-PMHT一類方法的收斂結果與迭代初值密切相關,而初值由前一時刻的估值所確定。因為噪聲影響,6 dB時初值誤差隨時間積累加快,會導致這類方法軌跡丟失率增大。從結果看到,H-PMHT/SN方法跟蹤效果也受初值誤差累積的影響,但它利用了噪聲信息,因而仍保持了較低的軌跡丟失率。該結果還可以由圖4體現,H-PMHT和P-HPMHT方法的軌跡在末尾處偏差較大。另外,可以由圖5發現JPDA方法在6 dB時的表現最差,其軌跡丟失率在軌跡交叉之后并未降低。這是因為在6 dB時該方法很難再糾正軌跡交叉時被“拉偏”的軌跡。

圖4 場景2中的典型跟蹤結果

圖5 場景2中的軌跡丟失率

場景3與2類似。所不同的是,2個目標的運動軌跡保持平行。在場景3中,一個典型跟蹤結果見圖6。圖中,目標1和目標2從左下角向右上角運動,其間距離保持為2個分辨單元。同場景2一樣,分別在信噪比10 dB和6 dB下,比對4種跟蹤方法的軌跡丟失率,見圖7。當信噪比為10 dB時,前3種方法都保持了較低的軌跡丟失率。但JPDA方法卻表現較差。究其原因,目標在場景3中始終保持較小的間距,因而通過JPDA方法跟蹤目標,軌跡有更大的概率被另一目標點跡“拉走”。這一現象也可以在圖6中體現,其中不同目標軌跡發生交叉。當信噪比為6 dB 時,噪聲對跟蹤結果的影響更為突出。H-PMHT和P-HPMHT這2種方法的軌跡丟失率都出現了十分明顯的上升,H-PMHT/SN方法卻仍有能力保持較低的軌跡丟失率。這時該方法利用噪聲統計特性的能力得到了體現。同樣的,還可以發現這種程度的信噪比條件下JPDA方法的表現最差。

圖6 場景3中的典型跟蹤結果

圖7 場景3中的軌跡丟失率

4 結 論

1) 本文提出了一種帶噪聲模型的H-PMHT方法,即H-PMHT/SN方法。該方法較傳統的H-PMHT方法和P-HPMHT方法在低信噪比條件下有著更強跟蹤檢測能力。同時,在多目標情況下,其計算量與目標數仍能保持線性關系。這為其實際應用提供了廣闊空間。

2) 該方法對目標狀態、強度參數及量化單位同時進行了估計,其中量化單位的估計準則為極大似然。在理論層面,它解決了量化單位真實值未知而導致的估計結果偏差問題。同時,還為工程中避免再使用“修正先驗信息”這一變通方法提供了一種可行的解決辦法。

3) 量化單位的估計值可能被用作構建懲罰函數來辨識目標總數M[14,15]。另外,本文并沒有給出Θk估計誤差的協方差陣,這將影響其濾波效果。這些問題有待進一步研究。