生雪莉，陳洋，郭龍祥，郝豪言，周媛媛，殷敬偉

(1.哈爾濱工程大學(xué) 水聲技術(shù)重點(diǎn)實(shí)驗(yàn)室，黑龍江 哈爾濱 150001； 2.海洋信息獲取與安全工信部重點(diǎn)實(shí)驗(yàn)室(哈爾濱工程大學(xué))，工業(yè)和信息化部，黑龍江 哈爾濱 150001； 3.哈爾濱工程大學(xué) 水聲工程學(xué)院，黑龍江 哈爾濱 150001)

水聲目標(biāo)跟蹤一直都是聲吶領(lǐng)域中重要研究的課題。近幾十年來(lái)，隨著海洋船舶技術(shù)和艦船降噪技術(shù)的快速發(fā)展，以及人類海洋活動(dòng)的頻繁，目標(biāo)輸入信號(hào)的信噪比不斷下降，導(dǎo)致水聲目標(biāo)跟蹤技術(shù)朝著分布式融合方向進(jìn)行了發(fā)展[1]。目前，傳統(tǒng)的目標(biāo)跟蹤算法，如多假設(shè)跟蹤(multiple hypotheses tracking, MHT)算法、聯(lián)合概率數(shù)據(jù)關(guān)聯(lián)(joint probabilistic data association, JPDA)等算法通常需要進(jìn)行復(fù)雜的關(guān)聯(lián)計(jì)算，因此當(dāng)量測(cè)個(gè)數(shù)較多時(shí)，算法的計(jì)算復(fù)雜度大幅增加，嚴(yán)重時(shí)可能無(wú)法進(jìn)行計(jì)算[2-3]。然而，基于RFS理論的概率假設(shè)密度(probability hypothesis density, PHD)濾波器由于不需要進(jìn)行關(guān)聯(lián)計(jì)算，因此計(jì)算效率大幅提高。例如高斯混合概率假設(shè)密度(gaussian mixture probability hypothesis density, GMPHD)濾波器[4]，其計(jì)算復(fù)雜度與量測(cè)個(gè)數(shù)和次數(shù)均呈線性關(guān)系[5]。因此，這類方法在假目標(biāo)干擾嚴(yán)重的場(chǎng)景中解決水聲目標(biāo)探測(cè)問(wèn)題具有重要的應(yīng)用價(jià)值[6]。例如，洛馬公司的反潛多目標(biāo)跟蹤識(shí)別系統(tǒng)和海軍研究辦公室的被動(dòng)聲吶反潛信息融合和目標(biāo)識(shí)別系統(tǒng)[7]。但是，這種方法也存在一個(gè)明顯的缺點(diǎn)，那就是只能生成點(diǎn)跡估計(jì)，而無(wú)法提供目標(biāo)的航跡信息，這給后續(xù)的航跡融合、航跡管理等工作帶來(lái)了困難。

在解決PHD算法航跡問(wèn)題的方法中，主流的做法有2個(gè)方向[8]。一類是通過(guò)關(guān)聯(lián)算法來(lái)解決PHD的目標(biāo)航跡生成問(wèn)題。文獻(xiàn)[9-10]利用數(shù)據(jù)關(guān)聯(lián)算法對(duì)序貫蒙特卡洛概率假設(shè)密度(sequential monte carlo probability hypothesis density, SMCPHD)濾波器進(jìn)行改進(jìn)，解決了SMCPHD濾波器的目標(biāo)航跡生成問(wèn)題。文獻(xiàn)[11]則對(duì)比了SMCPHD和MHT算法性能，并發(fā)現(xiàn)SMCPHD生成的航跡相比傳統(tǒng)MHT算法生成的航跡效果更好。Clark還對(duì)SMCPHD和GMPHD跟蹤濾波器的航跡生成問(wèn)題給出了多種解決方案[12-13]，并成功解決了基于聲吶圖像數(shù)據(jù)的多目標(biāo)跟蹤問(wèn)題[14]。

另一類是給PHD濾波器的目標(biāo)狀態(tài)增加一個(gè)標(biāo)記，通過(guò)標(biāo)記信息確定目標(biāo)的航跡。例如文獻(xiàn)[16]通過(guò)增加一個(gè)額外的標(biāo)記來(lái)區(qū)分不同的航跡，解決了航跡生成問(wèn)題。文獻(xiàn)[17]不僅利用標(biāo)記信息解決了航跡問(wèn)題，還提高了PHD濾波器剪枝算法的效率。文獻(xiàn)[18]則利用高斯項(xiàng)的標(biāo)記集進(jìn)行標(biāo)記的假設(shè)關(guān)聯(lián)，并采用匈牙利算法尋求最優(yōu)的假設(shè)關(guān)系，為GMPHD濾波器提供一個(gè)更加穩(wěn)定的目標(biāo)航跡。文獻(xiàn)[19]則對(duì)GMPHD濾波器中的高斯項(xiàng)進(jìn)行標(biāo)記，并利用樹狀結(jié)構(gòu)傳遞這些標(biāo)記，完成目標(biāo)航跡標(biāo)記的更新，解決航跡生成問(wèn)題。此外，Karl[20]和Han[21]等則在擴(kuò)展目標(biāo)的航跡問(wèn)題進(jìn)行深入研究，并提出了相應(yīng)的解決方法。

借鑒上述方法經(jīng)驗(yàn)，本文以GMPHD濾波器為基礎(chǔ)，通過(guò)將數(shù)據(jù)關(guān)聯(lián)方法與標(biāo)記法相結(jié)合，提出了一種適用于多聲吶多目標(biāo)航跡的生成融合方法，旨在解決單傳感器受時(shí)空非均勻傳播問(wèn)題，降低環(huán)境噪聲的影響。此外，本文還利用標(biāo)記信息構(gòu)建了一種簡(jiǎn)易的延時(shí)濾波器，解決了PHD濾波器估計(jì)遺漏的假目標(biāo)問(wèn)題。

1 PHD理論

在目標(biāo)跟蹤問(wèn)題中，貝葉斯濾波方法是當(dāng)今最為常用的方法。令fk|k-1(xk|xk-1)表示狀態(tài)x從k-1時(shí)刻到k時(shí)刻的馬爾可夫轉(zhuǎn)移密度，gk(·|·)表示似然函數(shù)，pk表示k時(shí)刻的后驗(yàn)概率密度，z表示量測(cè)，那么根據(jù)貝葉斯濾波理論有:

(1)

(2)

將其推廣到在多目標(biāo)的探測(cè)環(huán)境下，首先根據(jù)RFS理論，可以將k時(shí)刻的目標(biāo)狀態(tài)和傳感器的量測(cè)分別建立為一個(gè)隨機(jī)有限集，即：

Xk={xk,1,xk,2,…,xk,N(k)}∈F(χ)

(3)

Zk={zk,2,zk,2,…,zk,M(k)}∈F(Z)

(4)

式中：x表示目標(biāo)狀態(tài)；z表示傳感器的量測(cè)；F(Z)和F(χ)為單目標(biāo)量測(cè)z和狀態(tài)x的全部有限子集的集合。

考慮到k時(shí)刻的目標(biāo)狀態(tài)可能有3種：新生、衍生以及滅亡，那么目標(biāo)狀態(tài)集Xk還可以表示為：

(5)

式中：Sk|k-1、Bk|k-1和Γk分別表示為從k-1時(shí)刻存活、衍生到k時(shí)刻的目標(biāo)和k時(shí)刻新生目標(biāo)。

同時(shí)，在量測(cè)過(guò)程中還可能伴隨漏檢、錯(cuò)誤量測(cè)以及干擾的影響，則k時(shí)刻的量測(cè)集可表示為：

(6)

式中：Kk和Θk分別表示為第k時(shí)刻的錯(cuò)誤量測(cè)和干擾。

因此，可將單目標(biāo)的貝葉斯濾波推廣到多目標(biāo)隨機(jī)有限集的貝葉斯濾波，即:

(7)

(8)

式中μs為F(X)的一個(gè)參考指標(biāo)。

最后，通過(guò)MMSE(minimum mean squared error)或MAP(maximum a posteriori)準(zhǔn)則便可估計(jì)出目標(biāo)k時(shí)刻的狀態(tài)?？梢?，PHD濾波器避免了傳統(tǒng)跟蹤算法中的關(guān)聯(lián)過(guò)程，提高了計(jì)算性能，但其遞推公式卻存在解析解難求的問(wèn)題。為此，Ba-Ngu Vo提出了采用粒子近似的SMCPHD濾波器和采用高斯混合近似的GMPHD濾波器。本文假設(shè)系統(tǒng)服從高斯模型，因此選用計(jì)算效率更高的GMPHD濾波器作為本文的核心濾波器，其具體實(shí)現(xiàn)過(guò)程可參見文獻(xiàn)[4]。

2 標(biāo)記關(guān)聯(lián)的多聲吶航跡融合

2.1 分布式多聲吶數(shù)據(jù)融合結(jié)構(gòu)

如圖1所示，在分布式多傳感器融合方法框架下，本文首先利用GMPHD濾波器為每個(gè)聲吶的量測(cè)進(jìn)行獨(dú)立濾波，得到局部估計(jì)(點(diǎn)跡)。其次，為局部估計(jì)增加一個(gè)特殊的標(biāo)記，并利用關(guān)聯(lián)算法確定局部估計(jì)間的關(guān)聯(lián)關(guān)系，完成標(biāo)記更新并形成局部航跡。再次，利用關(guān)聯(lián)算法確定不同聲吶局部估計(jì)之間的關(guān)聯(lián)關(guān)系，并利用標(biāo)記關(guān)聯(lián)關(guān)系優(yōu)化關(guān)聯(lián)計(jì)算過(guò)程以及對(duì)錯(cuò)誤標(biāo)記進(jìn)行糾錯(cuò)。最后，根據(jù)關(guān)聯(lián)關(guān)系，通過(guò)多傳感器融合算法完成航跡融合，得到融合后的全局航跡。此外，利用全局航跡中的標(biāo)記信息，本文還給出了一種簡(jiǎn)易的延時(shí)濾波器，可以濾掉存活時(shí)間較短的假目標(biāo)航跡，令全局航跡估計(jì)結(jié)果更加可靠。

圖1 分布式多聲吶數(shù)據(jù)融合結(jié)構(gòu)Fig.1 Distributed multi-sonar data fusion structure diagram

2.2 基于標(biāo)記的航跡關(guān)聯(lián)方法

在本文中，主要選用一種基于馬氏距離的最近領(lǐng)域法(nearest neighbor，NN)來(lái)解決數(shù)據(jù)關(guān)聯(lián)問(wèn)題，即確定點(diǎn)跡與點(diǎn)跡、點(diǎn)跡與航跡以及航跡之間的關(guān)聯(lián)關(guān)系[22]。若用x1和x2表示待關(guān)聯(lián)的2個(gè)狀態(tài)，p1和p2表示其對(duì)應(yīng)的誤差協(xié)方差矩陣，則可以計(jì)算出2個(gè)狀態(tài)之間的馬氏距離α：

α=(x1-x2)T(p1+p2)(x1-x2)

(9)

那么，通過(guò)設(shè)立門限Tα，當(dāng)α小于門限Tα?xí)r，則認(rèn)為2個(gè)狀態(tài)是相互關(guān)聯(lián)的，其中α值越小，2個(gè)狀態(tài)的關(guān)聯(lián)程度越高。

2.2.1 標(biāo)記的局部聲吶航跡關(guān)聯(lián)

(10)

(11)

2.2.2 標(biāo)記的全局多聲吶航跡關(guān)聯(lián)

圖2 分布式多聲吶航跡關(guān)聯(lián)示意Fig.2 Distributed multi-sonar track association diagram

表1 全局關(guān)聯(lián)表AHTable 1 Globally associated table AH

可以看出，利用AH表可以優(yōu)化關(guān)聯(lián)計(jì)算問(wèn)題，提高關(guān)聯(lián)計(jì)算效率，即通過(guò)“行”排列便可找出同一聲吶(包含融合中心)中各個(gè)航跡的標(biāo)記；通過(guò)“列”排列便可找出同一航跡對(duì)應(yīng)不同聲吶中的局部航跡標(biāo)記的關(guān)聯(lián)關(guān)系。

再次，若局部估計(jì)為一個(gè)新目標(biāo)的標(biāo)記(即AH表中未確定該標(biāo)記與全局標(biāo)記間的關(guān)聯(lián)關(guān)系)，則可以利用式(9)～(11)，采用多次外推的方法確定該新生目標(biāo)的局部標(biāo)記與全局標(biāo)記間的關(guān)聯(lián)關(guān)系。若無(wú)法找到對(duì)應(yīng)的關(guān)聯(lián)關(guān)系，則認(rèn)為系統(tǒng)檢測(cè)到一個(gè)新生的目標(biāo)航跡，并對(duì)AH表的記錄進(jìn)行更新;否則，則根據(jù)AH表進(jìn)行標(biāo)記糾正。

最后，由于在全局航跡中，其標(biāo)記的第2項(xiàng)參數(shù)已經(jīng)記錄了多聲吶第1次檢測(cè)到該目標(biāo)航跡的起始時(shí)刻，因此，僅需要為該全局航跡添加一個(gè)結(jié)束時(shí)刻的標(biāo)記，即可實(shí)現(xiàn)對(duì)該航跡存活時(shí)間的記錄，即每當(dāng)存在該航跡時(shí)，則對(duì)該航跡的結(jié)束時(shí)刻進(jìn)行更新。

2.3 基于標(biāo)記的多聲吶數(shù)據(jù)融合和延時(shí)濾波器

在多傳感器數(shù)據(jù)融合方法中，簡(jiǎn)單凸組合融合算法是一種容易工程實(shí)現(xiàn)，且當(dāng)兩傳感器局部估計(jì)誤差不相關(guān)時(shí)，可以得到最優(yōu)估計(jì)的融合方法[23]。但是，該方法并沒有考慮聲吶漏檢等問(wèn)題。因此，本文主要根據(jù)AH表中的關(guān)聯(lián)關(guān)系，將加權(quán)融合策略和簡(jiǎn)單凸組合融合方法相結(jié)合，完成多聲吶的數(shù)據(jù)融合，得到全局航跡估計(jì)。

2.3.1 簡(jiǎn)單凸組合融合技術(shù)

(12)

將式(12)推廣到多聲吶(n>2)的情況時(shí)，則有：

(13)

2.3.2 加權(quán)融合策略

為了實(shí)現(xiàn)整個(gè)探測(cè)系統(tǒng)具有最大的檢測(cè)能力，本文的加權(quán)融合策略為：若多個(gè)聲吶的局部估計(jì)可以相互關(guān)聯(lián)且關(guān)聯(lián)程度最小時(shí)，則利用簡(jiǎn)單凸組合融合方法進(jìn)行數(shù)據(jù)融合，得到全局航跡估計(jì)；否則，由于無(wú)法判斷該目標(biāo)狀態(tài)是否為其他聲吶漏檢的目標(biāo)，因此本文將所有沒有關(guān)聯(lián)上(α>Tα)的局部航跡也保留作為全局航跡估計(jì)。需要說(shuō)明的是，這個(gè)加權(quán)融合策略可能會(huì)導(dǎo)致全局航跡中存在假目標(biāo)航跡，因此需要進(jìn)行后續(xù)的濾波處理。

2.3.3 基于標(biāo)記的延時(shí)濾波器

通常，真實(shí)目標(biāo)的航跡應(yīng)該是相對(duì)連續(xù)的，而干擾的假目標(biāo)的航跡是相對(duì)孤立的。因此，本文借助全局航跡的標(biāo)記信息，通過(guò)計(jì)算每個(gè)全局航跡的存活時(shí)間和被檢測(cè)到的次數(shù)，便可以比較為容易地濾除掉假目標(biāo)。例如，可以統(tǒng)計(jì)第k到第k-s幀內(nèi)，同一目標(biāo)航跡被檢測(cè)到的次數(shù)。若統(tǒng)計(jì)的次數(shù)少于aPd·s，則判定為假目標(biāo)，直接濾除即可，其中a為常數(shù)，Pd為檢測(cè)概率。

3 數(shù)值仿真

本文將借助OSPA(optimal sub-patten assignment)指標(biāo)對(duì)本算法性能進(jìn)行檢驗(yàn)。根據(jù)OSPA的定義，若令X={x1,x2,…,xn}和G={g1,g2,…,gm}分別表示估計(jì)出的目標(biāo)狀態(tài)集和目標(biāo)狀態(tài)真值集。其中，x和g分別代表目標(biāo)狀態(tài)和真值，那么，OSPA指標(biāo)可以被定義為[24]：

(14)

式中:dc(x,x′)=min{c,d(x,x′)}為相關(guān)程度的截止距離度量;c為關(guān)聯(lián)截至半徑;p為一個(gè)無(wú)量綱實(shí)數(shù)。

3.1 常規(guī)觀測(cè)場(chǎng)景

在[-1 000 m,1 000 m]×[-1 000 m,1 000 m]的監(jiān)視場(chǎng)景下，假設(shè)2聲吶均可以對(duì)該范圍所有的出現(xiàn)的目標(biāo)進(jìn)行探測(cè)。其檢測(cè)概率為Pd=0.9。目標(biāo)狀態(tài)用xk=[px,kpy,kvx,kvy,k]T來(lái)表示，其中(px,k,py,k)為目標(biāo)在k時(shí)刻的位置狀態(tài)，(vx,k,vy,k)為對(duì)應(yīng)的矢量速度。在連續(xù)觀測(cè)100幀數(shù)據(jù)的試驗(yàn)條件下，共模擬3個(gè)勻速直線運(yùn)動(dòng)的目標(biāo)。其中，目標(biāo)1在第1幀位置(-900,-900)出現(xiàn)，在第80幀終止在位置(600,600)處；目標(biāo)2在第1幀位置(-900,0)出現(xiàn)，在第100幀終止在位置(500,-900)處；目標(biāo)3在第60幀以速度[-20;0](m/s)在位置(300,0)處進(jìn)行移動(dòng)。算法考察近5幀的歷史數(shù)據(jù)，即s=5。

假設(shè)在監(jiān)視范圍內(nèi)，共模擬2個(gè)聲吶進(jìn)行聯(lián)合探測(cè)。其中，每個(gè)聲吶的探測(cè)過(guò)程都是相互獨(dú)立的，并且聲吶的量測(cè)數(shù)據(jù)中每幀最多包含50個(gè)假目標(biāo)干擾，且假目標(biāo)位置服從均勻分布。馬爾可夫轉(zhuǎn)移狀態(tài)矩陣F、過(guò)程噪聲協(xié)方差矩陣Q和量測(cè)噪聲協(xié)方差矩陣R分別為：

R=diag(100,100,100,100)(m2)。

式中dt=1 s為觀測(cè)周期。

如圖3所示為本次仿真試驗(yàn)結(jié)果。可以看出，圖3(a)和(b)體現(xiàn)了GMPHD具有良好的目標(biāo)跟蹤濾波能力；圖3(c)的仿真結(jié)果不僅證明了本方法能夠?qū)崿F(xiàn)多聲吶的數(shù)據(jù)融合，得到全局估計(jì)結(jié)果，還很好地形成了目標(biāo)航跡；圖3(d)的OSPA分析結(jié)果則表明了全局估計(jì)結(jié)果相比任意局部聲吶的估計(jì)跟蹤結(jié)果更加精確、穩(wěn)定(均值更小、方差更小)。

圖3 常規(guī)觀測(cè)場(chǎng)景下的仿真結(jié)果Fig.3 Simulation results under conventional observation

3.2 間斷觀測(cè)場(chǎng)景

在3.1節(jié)的仿真場(chǎng)景下，本文還考察了一種間斷的觀測(cè)場(chǎng)景。其中，假設(shè)聲吶1從第20幀～第70幀受到海洋信道等因素影響，導(dǎo)致該聲吶無(wú)法探測(cè)到目標(biāo)1，形成對(duì)目標(biāo)1的間斷觀測(cè)，其仿真結(jié)果如圖4所示。

對(duì)比圖4(a)和(c)可以看出，雖然聲吶1在一段時(shí)間內(nèi)都無(wú)法探測(cè)到目標(biāo)1，但是當(dāng)該聲吶再次探測(cè)到目標(biāo)1時(shí)，理論上會(huì)判定為一個(gè)新目標(biāo)的航跡(長(zhǎng)時(shí)間未檢測(cè)到該目標(biāo))。但是，由于本算法在航跡生成融合過(guò)程中，借助AH表對(duì)航跡的標(biāo)記進(jìn)行了糾錯(cuò)，因此可以得到一個(gè)正確的局部航跡，如圖中4(c)航跡1所示。從圖4(b)、(e)和(f)中可以發(fā)現(xiàn)，當(dāng)2個(gè)聲吶至少有一個(gè)聲吶可以探測(cè)到目標(biāo)時(shí)，本方法便可以得到一個(gè)相對(duì)穩(wěn)定、連續(xù)的目標(biāo)航跡。

3.3 延時(shí)濾波器

由于本文的加權(quán)融合策略可能會(huì)導(dǎo)致全局航跡中存在比較明顯的假目標(biāo)航跡。為此，在本節(jié)中主要通過(guò)仿真驗(yàn)證延時(shí)濾波器的作用。在上一節(jié)仿真條件下，延時(shí)濾波器的檢測(cè)常數(shù)a=0.5。其仿真結(jié)果如圖5所示。

在圖5(a)中，聲吶1在(-450,-150)處附近產(chǎn)生了一個(gè)明顯的錯(cuò)誤估計(jì)。由于本文的加權(quán)融合策略，導(dǎo)致了全局航跡估計(jì)中產(chǎn)生了對(duì)應(yīng)的錯(cuò)誤航跡(如圖5(c)所示)。但是，利用本文的延時(shí)濾波器，可以濾掉將這種存活時(shí)間較短的假目標(biāo)航跡，降低了加權(quán)融合策略所導(dǎo)致的全局估計(jì)中存在假目標(biāo)的可能性，進(jìn)一步保證了算法估計(jì)的正確性。當(dāng)然，也可以利用MHT或JPDA方法進(jìn)行二次濾波來(lái)濾掉這種存活時(shí)間比較短的假目標(biāo)，但是，這種做法顯然沒有統(tǒng)計(jì)航跡標(biāo)記個(gè)數(shù)簡(jiǎn)單、易實(shí)現(xiàn)。

圖4 間斷觀測(cè)場(chǎng)景下的仿真結(jié)果Fig.4 Simulation results under intermittent observation scene

圖5 延時(shí)濾波仿真結(jié)果Fig.5 Simulation results of delay filtering

4 結(jié)論

1)即使某個(gè)聲吶無(wú)法對(duì)目標(biāo)進(jìn)行連續(xù)觀測(cè)，也可以正確地確定該目標(biāo)航跡信息。

2)基于標(biāo)記的延時(shí)濾波器不僅可以進(jìn)一步濾除局部聲吶中的假目標(biāo)干擾，還具有工程易實(shí)現(xiàn)等優(yōu)點(diǎn)。在OSPA評(píng)估算法的幫助下，仿真試驗(yàn)從常規(guī)觀測(cè)和間斷觀測(cè)2個(gè)場(chǎng)景檢驗(yàn)了該方法的性能，證明了該方法具有形成準(zhǔn)確、連續(xù)目標(biāo)航跡的能力。此外，延時(shí)濾波器的仿真也說(shuō)明了這種濾波器能夠有效地濾除存活時(shí)間較短的假目標(biāo)干擾。