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

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

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

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

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

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

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

本文提出了一種帶傳感器噪聲模型的H-PMHT方法：通過選擇恰當(dāng)?shù)牧繙y模型丟失數(shù)據(jù),實現(xiàn)了引入噪聲模型條件下尋優(yōu)過程的解耦。又通過提出一種類泊松概率質(zhì)量函數(shù)多重卷積簡化計算方法,實現(xiàn)了尋優(yōu)過程的高效計算。最終,在引入噪聲模型條件下實現(xiàn)了計算量隨著目標(biāo)數(shù)線性增長。且同時,這一方法還給出了一種理論化的量化單位尋優(yōu)辦法。

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

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

(1)

(2)

(3)

為第k時刻強度參數(shù)的集合。

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

(4)

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

(5)

(6)

與之對應(yīng)的傳感器輸出為

(7)

(8)

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

(9)

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

(10)

(11)

(12)

(13)

令量化輸出

(14)

(15)

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

(16)

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

令參數(shù)集

Θk={Xk,Λk,η}

(17)

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

(18)

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

2.1 EM方法

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

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

(19)

(20)

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

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

由式(8)可知,計數(shù)值

(21)

其中,

(22)

式中：

(23)

(24)

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

而后,分別定義集合

(25)

i=1,2,…,I}

(26)

這里丟失數(shù)據(jù)為

(27)

2.2 期望步驟

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

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

(28)

式中：等號右側(cè)第1項為Ζk和Ok的對數(shù)聯(lián)合條件概率密度函數(shù)；第2項為Θk的對數(shù)先驗分布。

(29)

(30)

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

(31)

(32)

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

(33)

(34)

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

(35)

式中：C為與Θk無關(guān)的量。

假設(shè)各目標(biāo)狀態(tài)間相互獨立,Xk的先驗密度函數(shù)為

(36)

式中：pX為目標(biāo)狀態(tài)的先驗密度函數(shù)。

假設(shè)Λk的先驗密度函數(shù)為

(37)

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

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

(38)

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

(39)

進(jìn)行整理,輔助函數(shù)Q為

(40)

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

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

(41)

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

2.3 極大化步驟

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

(42)

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

(43)

(44)

(45)

假設(shè)目標(biāo)狀態(tài)的演化規(guī)律可以用一階馬爾科夫模型來描述,則由C-K方程可得

(46)

(47)

(48)

式中：

(49)

(50)

(51)

(52)

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

(53)

(54)

式中：ξ為遺忘因子；Δ為不同采樣時刻時間間隔?？梢缘玫狡涓路匠虨?/p>

(55)

(56)

2.3.2 求Qη的極值

假設(shè)傳感器噪聲為高斯噪聲,即

(57)

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

對Qη表達(dá)式(45)中的對數(shù)項進(jìn)行展開,得

(58)

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

(59)

對式(59)求導(dǎo)并令其為0,進(jìn)而得到η的估計值為

(60)

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

引理定義函數(shù)族

(61)

(62)

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

(63)

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

(64)

(65)

(66)

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

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

(67)

式中：系數(shù)

(68)

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

(69)

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

(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)

式中：系數(shù)c1為

(78)

它由式(75)計算得到。

當(dāng)m1=m2=m時,其條件期望為

(79)

當(dāng)m1≠m2時,其條件期望值為

(80)

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

(81)

式中：當(dāng)m1=m2=m時,系數(shù)cq依次為

(82)

當(dāng)m1≠m2時,為

(83)

進(jìn)一步,利用以上結(jié)果還易得

(84)

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

(85)

總結(jié)以上過程,對Θk進(jìn)行估計,其步驟如下：

步驟1初始設(shè)置

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

設(shè)定迭代次數(shù)計數(shù)值d=0。

步驟2迭代更新

更新迭代次數(shù)計數(shù)值d=d+1。

步驟3回到步驟2,直到估計結(jié)果收斂。

3 仿真驗證

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

(86)

式中：F為狀態(tài)轉(zhuǎn)移矩陣

(87)

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

(88)

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

假設(shè)第m號目標(biāo)的能量分布為高斯密度函數(shù)

(89)

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

(90)

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

(91)

進(jìn)而,定義目標(biāo)的信噪比為

i,j=1,2,…,I

(92)

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

圖1 仿真數(shù)據(jù)示意

3.1 仿真場景

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

3.2 跟蹤實現(xiàn)及性能評價

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

1) 軌跡起始：設(shè)定檢測概率為0.9,在給定的信噪比條件下確定檢測門限。分割第k幀數(shù)據(jù)以得到該時刻的目標(biāo),并與前一幀目標(biāo)關(guān)聯(lián)構(gòu)建新的軌跡。其條件為相應(yīng)的目標(biāo)間的距離小于1.5Δl。

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

(93)

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

圖2 ROC曲線

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

圖3 RMSE圖示

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

該方法在2個目標(biāo)交叉運動時的跟蹤能力。

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

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

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

圖4 場景2中的典型跟蹤結(jié)果

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

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

圖6 場景3中的典型跟蹤結(jié)果

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

4 結(jié) 論

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

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

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