張博宇,齊 濱1,2,,王晉晉1,2,,梁國龍1,2,

(1.哈爾濱工程大學(xué)水聲技術(shù)全國重點實驗室,哈爾濱 150001;2.海洋信息獲取與安全工信部重點實驗室(哈爾濱工程大學(xué)),工業(yè)和信息化部,哈爾濱 150001;3.哈爾濱工程大學(xué)水聲工程學(xué)院,哈爾濱 150001)

0 引言

水下多目標(biāo)跟蹤廣泛應(yīng)用于軍事和民用領(lǐng)域[1-2]。主動聲吶由于可以同時獲得目標(biāo)的方位和距離信息,因此在水下探測和跟蹤領(lǐng)域備受關(guān)注[3-5]。對于水下弱目標(biāo)檢測來說,通常需要設(shè)置較低的檢測門限,由此使得未知強干擾和噪聲能夠輕易越過門限,導(dǎo)致探測結(jié)果中存在大量的虛警,致使以探測結(jié)果為輸入的跟蹤算法處在密集雜波環(huán)境下。此時,跟蹤算法在初始化新生目標(biāo)時容易將雜波誤認(rèn)為新生目標(biāo),致使跟蹤結(jié)果中存在虛假目標(biāo),從而導(dǎo)致跟蹤算法出現(xiàn)性能退化。

常用的多目標(biāo)跟蹤方法可以分為兩類。第一類是數(shù)據(jù)關(guān)聯(lián)方法,包括多假設(shè)跟蹤[6-7](multiple hypothesis tracker, MHT)和聯(lián)合概率數(shù)據(jù)關(guān)聯(lián)[8-9](joi-nt probability data association, JPDA)等。由于這類方法在跟蹤時需要將測量值與目標(biāo)進行關(guān)聯(lián),因此其計算量隨著測量值和目標(biāo)數(shù)的增加而增大[10]。第二類是基于隨機有限集(random finite set, RFS)理論的多目標(biāo)跟蹤方法,包括概率假設(shè)密度[11](probability hypothesis density, PHD)濾波器,廣義標(biāo)簽多伯努利[12](generalized labeled multi-Bernoulli, GLMB)濾波器和泊松多伯努利混合[13](Poisson multi-Bernoulli mixture, PMBM)濾波器等。這類方法在跟蹤時避免了測量值與目標(biāo)之間的關(guān)聯(lián),與數(shù)據(jù)關(guān)聯(lián)類方法相比,運算效率顯著提高。在基于RFS理論的跟蹤方法中,PHD濾波器具有最低的計算復(fù)雜度[10],因此更適用于解決密集雜波背景下的多目標(biāo)跟蹤問題。PHD濾波器可以通過高斯混合(Gaussian mixture PHD, GM-PHD)[14-16]或序貫蒙特卡羅(sequential Monte Carlo PHD, SMC-PHD)[17- 18]的方式來實現(xiàn)?？紤]到SMC-PHD濾波器根據(jù)聚類算法提取目標(biāo)狀態(tài)的方式會導(dǎo)致提取結(jié)果的不穩(wěn)定[10,14],本文利用GM-PHD濾波器實現(xiàn)水下多目標(biāo)跟蹤。

常規(guī)GM-PHD濾波器在跟蹤時需要先驗已知新生目標(biāo)強度,否則其性能會出現(xiàn)嚴(yán)重退化。由于測量值集合中包含新生目標(biāo)的信息,因此現(xiàn)有的研究方法大多利用測量值集合來估計新生目標(biāo)強度。根據(jù)利用測量值集合的不同,現(xiàn)有的新生目標(biāo)強度估計方法可以分為兩類。第一類是利用當(dāng)前時刻測量值集合估計新生目標(biāo)強度。Fu等[19]利用當(dāng)前時刻測量值集合初始化新生目標(biāo)強度,然后對估計的新生目標(biāo)強度進行融合處理。由于利用當(dāng)前時刻測量值集合只能估計出目標(biāo)的位置信息,為了提高估計精度,第二類方法利用連續(xù)兩個時刻的測量值集合同時估計新生目標(biāo)的位置和速度。Zhang 等[20]根據(jù)連續(xù)兩個時刻的測量值集合,采用兩步初始化方法估計新生目標(biāo)強度。為了充分利用測量值集合中包含的目標(biāo)信息,Zhu等[21]將兩步初始化方法與一步初始化方法相結(jié)合來估計新生目標(biāo)強度。這類方法雖然可以提高新生目標(biāo)強度的估計精度,但是在雜波密集的跟蹤場景中,測量值集合中包含大量的虛警,利用連續(xù)兩個時刻的測量值空間分布來估計新生目標(biāo)強度會出現(xiàn)將雜波誤認(rèn)為新生目標(biāo)的情況,從而導(dǎo)致跟蹤結(jié)果中存在虛假目標(biāo)。

針對上述問題,本文提出一種滑動窗兩步初始化GM-PHD(sliding window two step initialization GM-PHD, SWTSI-GMPHD)濾波器。該濾波器首先利用滑動時間窗構(gòu)造測量值集合,然后根據(jù)滑動窗內(nèi)測量值的空間分布采用兩步初始化方法對測量值進行平滑處理,最后用平滑得到的新生目標(biāo)信息來估計目標(biāo)強度。將滑動窗兩步初始化方法嵌入GM-PHD濾波器,可以提高密集雜波背景下的多目標(biāo)跟蹤精度。本文方法與基于當(dāng)前時刻或連續(xù)兩個時刻測量值估計新生目標(biāo)強度的GM-PHD濾波器相比,可以有效減少虛假目標(biāo)數(shù)目。

1 系統(tǒng)模型

高斯混合概率假設(shè)密度濾波器作為多目標(biāo)貝葉斯濾波器的近似,將復(fù)雜的多重積分轉(zhuǎn)移到單目標(biāo)狀態(tài)空間,在降低計算復(fù)雜度的同時實現(xiàn)了多目標(biāo)跟蹤[14]。因此,首先介紹單目標(biāo)系統(tǒng)模型,即狀態(tài)模型和量測模型。目標(biāo)的狀態(tài)模型采用如下遞歸高斯模型[22]

xk=Fk-1xk-1+nk-1

(1)

其中,目標(biāo)的狀態(tài)矢量為xk=[px,k,vx,k,py,k,vy,k]T;(px,k,py,k)表示目標(biāo)的位置;(vx,k,vy,k)表示目標(biāo)的速度。Fk-1表示狀態(tài)轉(zhuǎn)移矩陣;nk-1表示均值為0,方差為Qk-1的高斯白噪聲。Fk-1和Qk-1的表達式分別為

(2)

其中,T表示采樣間隔;σv表示過程噪聲的標(biāo)準(zhǔn)差。由于主動聲吶可以獲得目標(biāo)的方位和距離信息,因此量測模型可以表示為[23]

(3)

zXY,k=Hxk+wXY,k

(4)

其中,偽線性量測向量zXY,k和偽線性量測矩陣H的表達式分別為

(5)

2 滑動窗兩步初始化高斯混合概率假設(shè)密度濾波器

盡管GM-PHD濾波器具有較低的計算復(fù)雜度,但是在跟蹤時需要先驗已知新生目標(biāo)的強度函數(shù),否則其性能將出現(xiàn)嚴(yán)重退化。在實際跟蹤場景中,新生目標(biāo)的強度通常未知,為此文獻[19-21]提出了自適應(yīng)估計新生目標(biāo)強度的GM-PHD濾波器。然而現(xiàn)有的跟蹤方法均未考慮雜波干擾對新生目標(biāo)強度的影響,因此在雜波密集的水下多目標(biāo)跟蹤場景中,這些跟蹤方法會將雜波誤認(rèn)為新生目標(biāo),從而導(dǎo)致跟蹤結(jié)果中存在虛假目標(biāo)。針對上述問題,提出一種滑動窗兩步初始化方法來預(yù)測新生目標(biāo)強度,將該方法嵌入GM-PHD濾波器,可以提高密集雜波背景下的多目標(biāo)跟蹤精度。

2.1 預(yù)測過程

預(yù)測過程包括預(yù)測存活目標(biāo)和新生目標(biāo),接下來分別介紹存活目標(biāo)和新生目標(biāo)的預(yù)測過程。

2.1.1 預(yù)測存活目標(biāo)

假設(shè)k-1時刻目標(biāo)的強度函數(shù)為[14, 25]

(6)

(7)

其中,ψk表示目標(biāo)存活概率。

2.1.2 預(yù)測新生目標(biāo)

(8)

(9)

(10)

(11)

其中,H表示偽線性量測矩陣;R(zXY,k)表示偽線性量測噪聲的方差;min(·)表示最小值函數(shù)。取集合ZXY,k與ZsXY,k的差集得到由新生目標(biāo)測量值和雜波組成的測量值集合ZreXY,k,即ZreXY,k=ZXY,k-ZsXY,k。由于雜波在量測空間任意分布,而新生目標(biāo)測量值分布在目標(biāo)真實位置附近,因此提出一種滑動窗兩步初始化方法來提取集合ZreXY,k中的新生目標(biāo)測量值。

假設(shè)滑動窗的窗長為φ,則在滑動窗內(nèi)的測量值集合ZXY,win可以表示為

ZXY,win={ZreXY,k,ZXY,k+1,…,ZXY,k+φ-1}

(12)

由于水下目標(biāo)的運動速度較慢,目標(biāo)在滑動窗內(nèi)的運動狀態(tài)可以近似為勻速直線運動,因此利用式(13)提取新生目標(biāo)測量值

φ=1,…,φ-1

(13)

(14)

(15)

表1 新生目標(biāo)強度估計方法

圖1 新生目標(biāo)強度估計方法的流程圖Fig.1 Flowchart of newborn target intensity estimation method

2.2 更新過程

預(yù)測過程得到的目標(biāo)強度函數(shù)為

vk|k-1(xk)=vs,k|k-1(xk)+vγ,k(xk)

(16)

其中,Jk|k-1表示預(yù)測的目標(biāo)數(shù)目。利用k時刻的偽測量值集合ZXY,k對預(yù)測的強度函數(shù)vk|k-1(xk)進行更新,可以得到更新的強度函數(shù)vk(xk)的表達式為[14]

(17)

其中,?k表示探測概率。

2.3 提取目標(biāo)狀態(tài)向量

為了提高濾波器的運算效率,采用文獻[14]中提出的方法對更新得到的強度函數(shù)vk(xk)中的高斯分量進行剪枝與合并,即保留vk(xk)中權(quán)值較大的高斯分量,并對狀態(tài)相似的高斯分量進行合并。假設(shè)經(jīng)過剪枝與合并過程后得到的強度函數(shù)vPM,k(xk)的表達式為

(18)

(19)

3 仿真實驗

為了充分驗證提出方法的有效性,將其跟蹤結(jié)果與一步初始化GM-PHD(one step initialization GM-PHD, OSI-GMPHD)濾波器[19],兩步初始化GM-PHD(two step initialization GM-PHD, TSI-GMPHD)濾波器[20],兩步與一步初始化相結(jié)合的GM-PHD(two step initialization-one step initialization GM-PHD, TSI-OSI-GMPHD)濾波器[21]進行對比。表2所示為仿真實驗參數(shù)設(shè)置,其中距離門限dth和滑動窗長度φ決定了新生目標(biāo)強度估計的準(zhǔn)確度。距離門限dth用來提取測量值集合ZXY,win中的新生目標(biāo)測量值,并由目標(biāo)的運動速度決定。當(dāng)dth的取值較小時,提取的新生目標(biāo)測量值存在漏報的情況,從而導(dǎo)致跟蹤算法的精度下降;當(dāng)dth的取值較大時,提取的新生目標(biāo)測量值存在雜波,致使跟蹤結(jié)果中存在虛假目標(biāo)?？紤]到水下目標(biāo)的運動速度通常較慢,因此在仿真實驗中將dth設(shè)置為30 m?；瑒哟伴L度φ決定測量值集合ZXY,win中包含的測量值數(shù)目。當(dāng)φ的取值較小時,利用測量值的空間分布提取新生目標(biāo)測量值時,會導(dǎo)致雜波包含在新生目標(biāo)測量值集合ZXY,new中,從而導(dǎo)致跟蹤結(jié)果中出現(xiàn)虛假目標(biāo);當(dāng)φ的取值較大時,目標(biāo)的機動導(dǎo)致其在滑動窗內(nèi)的運動狀態(tài)不滿足勻速運動的假設(shè)條件,致使提取的新生目標(biāo)測量值集合ZXY,new存在漏報的情況,從而導(dǎo)致估計的新生目標(biāo)強度函數(shù)與真實值之間存在偏差。綜上,將滑動窗長度φ設(shè)置為5。

3.1 算法有效性驗證

圖2所示為目標(biāo)機動場景下的運動態(tài)勢圖,其中聲吶靜止在原點。

圖2 目標(biāo)與聲吶的運動態(tài)勢圖Fig.2 Real trajectories of targets and sonar

假設(shè)雜波在[0,360]deg×[0,1 000]m的區(qū)間內(nèi)滿足泊松分布。雜波強度為Kk(zXY,k)=mc/V,其中mc表示雜波數(shù)目;V表示雜波分布的空間面積。將測量值由非線性量測空間轉(zhuǎn)換到偽線性量測空間后,其分布在半徑為1 000 m的圓形區(qū)域內(nèi),因此雜波分布的空間面積為V=π×106m2。將雜波數(shù)mc設(shè)置為20,可以得到如圖3所示的方位和距離估計結(jié)果。從圖3可以看出,在密集雜波環(huán)境下,測量結(jié)果中包含較多的雜波干擾,并且來源于目標(biāo)的測量值被淹沒在雜波中。

(a) 方位估計結(jié)果

(b) 距離估計結(jié)果圖3 單次實驗的方位和距離估計結(jié)果Fig.3 Bearing and range estimation results of a single experiment

將圖3所示的方位和距離估計結(jié)果作為跟蹤算法的輸入,可以得到如圖4所示的跟蹤結(jié)果和目標(biāo)數(shù)目估計結(jié)果。OSI-GMPHD濾波器利用當(dāng)前時刻測量值集合來估計新生目標(biāo)強度,受到測量值集合中的雜波干擾,其跟蹤結(jié)果中存在大量的虛假目標(biāo);TSI-GMPHD濾波器利用連續(xù)兩個采樣時刻的測量值集合估計新生目標(biāo)強度,相比于OSI-GMPHD濾波器,其跟蹤結(jié)果中的虛假目標(biāo)數(shù)目較少;TSI-OSI-GMPHD濾波器將兩步初始化與一步初始化方法相結(jié)合,充分利用測量值集合來估計新生目標(biāo)強度,具有較高的跟蹤精度。然而,當(dāng)雜波在兩個采樣時刻均出現(xiàn)的情況下,TSI-OSI-GMPHD濾波器會將雜波誤認(rèn)為新生目標(biāo),因此其跟蹤結(jié)果中仍然存在虛假目標(biāo)。本文提出的SWTSI-GMPHD濾波器,利用滑動窗兩步初始化方法可以提高新生目標(biāo)強度估計的準(zhǔn)確度,減少跟蹤結(jié)果中的虛假目標(biāo)數(shù)目。

(a) OSI-GMPHD

(b) TSI-GMPHD

(c) TSI-OSI-GMPHD

(d) SWTSI-GMPHD圖4 不同算法的跟蹤結(jié)果和目標(biāo)數(shù)目估計結(jié)果Fig.4 Tracking and target number estimation results of different algorithms

3.2 算法統(tǒng)計性能分析

為了進一步驗證提出方法的有效性,將上述仿真場景進行800次蒙特卡羅實驗,利用最優(yōu)子模式分配(optimal sub-pattern assignment, OSPA)準(zhǔn)則[26]和平均目標(biāo)數(shù)估計結(jié)果來衡量不同算法的跟蹤精度。OSPA的相關(guān)參數(shù)設(shè)置為:距離敏感性參數(shù)p=1,關(guān)聯(lián)敏感參數(shù)c=100。圖5所示為不同算法的OSPA誤差統(tǒng)計結(jié)果和目標(biāo)數(shù)目估計結(jié)果。從圖5中可以看出,OSI-GMPHD濾波器的OSPA誤差最大,并且估計的目標(biāo)數(shù)目與真實值之間存在較大偏差,因此OSI-GMPHD濾波器的跟蹤精度最低;相比于OSI-GMPHD濾波器,TSI-GMPHD濾波器的OSPA誤差較低,目標(biāo)數(shù)目估計結(jié)果與真實值之間的偏差較小;TSI-OSI-GMPHD濾波器將兩步初始化與一步初始化方法相結(jié)合,充分利用測量值集合中包含的信息,因此其跟蹤精度高于TSI-GMPHD濾波器;SWTSI-GMPHD濾波器可以有效提高新生目標(biāo)強度估計的準(zhǔn)確度,因此具有最低的OSPA誤差和最接近真實值的目標(biāo)數(shù)目估計結(jié)果。

(a) OSPA距離誤差

(b) 目標(biāo)數(shù)目估計結(jié)果圖5 蒙特卡羅實驗統(tǒng)計結(jié)果Fig.5 Statistical results of monte carlo experiments

表3所示為不同算法的平均OSPA誤差統(tǒng)計結(jié)果,與OSI-GMPHD,TSI-GMPHD和TSI-OSI-GMPHD算法相比,SWTSI-GMPHD將跟蹤精度分別提高69.84%,52.62%和41.05%。

表3 平均OSPA誤差統(tǒng)計結(jié)果

4 結(jié)論

本文通過研究,得到如下結(jié)論:

1)針對水下密集雜波環(huán)境導(dǎo)致多目標(biāo)跟蹤算法輸出的跟蹤結(jié)果中包含較多虛假目標(biāo)的問題,提出一種滑動窗兩步初始化高斯混合概率假設(shè)密度濾波器。

2)該方案利用滑動時間窗構(gòu)造測量值集合,然后根據(jù)滑動窗內(nèi)測量值的空間分布提取新生目標(biāo)信息,進而估計新生目標(biāo)強度。將滑動窗兩步初始化方法嵌入GM-PHD濾波器,可以有效減少跟蹤結(jié)果中的虛假目標(biāo)數(shù)目。

3)仿真實驗結(jié)果表明,相較于其他方法,提出方法將跟蹤精度提高69.84%,52.62%和41.05%。