密集雜波背景下的水下多目標跟蹤方法

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

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

0 引言

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

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

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

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

1 系統模型

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

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

(1)

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

(2)

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

(3)

zXY,k=Hxk+wXY,k

(4)

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

(5)

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

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

2.1 預測過程

預測過程包括預測存活目標和新生目標,接下來分別介紹存活目標和新生目標的預測過程。

2.1.1 預測存活目標

假設k-1時刻目標的強度函數為[14, 25]

(6)

(7)

其中,ψk表示目標存活概率。

2.1.2 預測新生目標

(8)

(9)

(10)

(11)

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

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

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

(12)

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

φ=1,…,φ-1

(13)

(14)

(15)

表1 新生目標強度估計方法

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

2.2 更新過程

預測過程得到的目標強度函數為

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

(16)

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

(17)

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

2.3 提取目標狀態向量

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

(18)

(19)

3 仿真實驗

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

3.1 算法有效性驗證

圖2所示為目標機動場景下的運動態勢圖,其中聲吶靜止在原點。

圖2 目標與聲吶的運動態勢圖Fig.2 Real trajectories of targets and sonar

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

(a) 方位估計結果

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

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

(a) OSI-GMPHD

(b) TSI-GMPHD

(c) TSI-OSI-GMPHD

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

3.2 算法統計性能分析

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

(a) OSPA距離誤差

(b) 目標數目估計結果圖5 蒙特卡羅實驗統計結果Fig.5 Statistical results of monte carlo experiments

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

表3 平均OSPA誤差統計結果

4 結論

本文通過研究,得到如下結論:

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

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

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