俞建國(guó) 劉 梅 王 駿

①(哈爾濱工業(yè)大學(xué)電子與信息工程學(xué)院 哈爾濱 150001)

②(中國(guó)電子科技集團(tuán)第29研究所 成都 610036)

1 引言

彈道導(dǎo)彈(Ballistic Missile, BM)在現(xiàn)代戰(zhàn)爭(zhēng)中具有穿透力強(qiáng)，威懾力大等特點(diǎn)，是各國(guó)空間預(yù)警和監(jiān)視系統(tǒng)重點(diǎn)關(guān)注的目標(biāo)。BM 為了提高突防能力，在再入過(guò)程中往往伴隨著多彈頭分導(dǎo)技術(shù)，并且在多批次飽和攻擊模式下，密集程度較高，與傳統(tǒng)多目標(biāo)數(shù)據(jù)關(guān)聯(lián)[1,2]相比，對(duì)分導(dǎo)彈頭數(shù)據(jù)關(guān)聯(lián)難點(diǎn)在于：(1)分導(dǎo)彈頭數(shù)目以及分導(dǎo)時(shí)間未知，彈頭之間距離較小，速度非常相近，易于形成團(tuán)狀；(2)彈頭再入速度大，距離打擊目標(biāo)近，攔截時(shí)間短。因此，必須尋求新的數(shù)據(jù)關(guān)聯(lián)方法用于解決彈道導(dǎo)彈分導(dǎo)特定應(yīng)用背景問(wèn)題。

目前運(yùn)用較為成熟的多目標(biāo)數(shù)據(jù)關(guān)聯(lián)方法為貝葉斯類(lèi)數(shù)據(jù)互聯(lián)。貝葉斯典型方法包括單步貝葉斯和多步貝葉斯方法。單步貝葉斯方法只對(duì)最新的量測(cè)集合進(jìn)行關(guān)聯(lián)，是一種次優(yōu)貝葉斯方法，比較有代表性的算法有最近鄰[3](Nearest-Neighbor, NN)和聯(lián)合概率數(shù)據(jù)互聯(lián)算法[4](Joint Probability Data Association, JPDA)。NN方法在目標(biāo)密集、數(shù)量變化以及有雜波和漏警的環(huán)境下效果難以令人滿(mǎn)意，而JPDA計(jì)算關(guān)聯(lián)概率權(quán)值是個(gè)NP-hard問(wèn)題。這兩種方法都不能處理目標(biāo)建航和航跡撤銷(xiāo)問(wèn)題。多步貝葉斯方法對(duì)當(dāng)前時(shí)刻以前的量測(cè)值進(jìn)行研究，給出每個(gè)量測(cè)序列的概率，以多假設(shè)跟蹤[5](Multiple Hypothesis Tracking, MHT)算法最具代表性。它通過(guò)對(duì)當(dāng)前時(shí)刻之前的量測(cè)序列進(jìn)行假設(shè)關(guān)聯(lián)，返回最大后驗(yàn)概率的假設(shè)關(guān)聯(lián)作為最優(yōu)關(guān)聯(lián)，是目前多目標(biāo)關(guān)聯(lián)應(yīng)用較為成熟的算法。它能夠解決不定數(shù)目的目標(biāo)建航和終結(jié)問(wèn)題，不足之處是通過(guò)枚舉的方法獲得所有關(guān)聯(lián)假設(shè)，計(jì)算量較大。

在此應(yīng)用背景下本文提出了一種改進(jìn)的實(shí)時(shí)滑窗 MCMC方法用于解決彈道導(dǎo)彈分導(dǎo)數(shù)據(jù)關(guān)聯(lián)問(wèn)題。該算法將彈頭分導(dǎo)視為平穩(wěn)的馬爾可夫過(guò)程，通過(guò)蒙特卡洛采樣的方法對(duì)監(jiān)控區(qū)域的測(cè)量集合進(jìn)行組合優(yōu)化，獲得最大的后驗(yàn)概率密度進(jìn)而逼近馬氏鏈的平穩(wěn)分布，結(jié)合彈頭分導(dǎo)的實(shí)際應(yīng)用背景進(jìn)行繼承性?xún)?yōu)化以及關(guān)聯(lián)假設(shè)概率權(quán)值重新分配，在保持較高的關(guān)聯(lián)概率下極大地提高了實(shí)時(shí)性。

2 濾波狀態(tài)方程與測(cè)量方程

2.1 狀態(tài)方程

彈道目標(biāo)的作戰(zhàn)窗口一般在被動(dòng)段，因此本文所用的濾波模型主要針對(duì)彈道目標(biāo)的被動(dòng)段，用f表示。BM 在被動(dòng)段主要受地球引力、大氣阻力以及慣性力作用，這里慣性力包括科里奧利力和向心力?？紤]地球的橢球性，保留球引力二階以下帶諧項(xiàng)，結(jié)合BM在該階段所受的4個(gè)力可得到地心固連參考坐標(biāo)系(Earth-Centered Earth-Fixed, ECEF)下的6個(gè)常微分方程用于描述被動(dòng)段模型[6,7]，方程組如下：

式(1)中m為導(dǎo)彈質(zhì)量，x,y,z和vx,vy,vz分別代表導(dǎo)彈在ECEF坐標(biāo)系下的3個(gè)方向的位置和速度，x,y,z分別為3個(gè)方向位置單位矢量，F(xiàn)g為重力，具體表達(dá)式可參考文獻(xiàn)[8]。||˙||為速度幅值，ω為非慣性坐標(biāo)的自轉(zhuǎn)角速度大小為 7 .292115 × 10-5rad/s。β為質(zhì)阻比，表征彈頭在再入過(guò)程中速度衰減快慢的物理量，一般與彈道目標(biāo)的質(zhì)量和幾何形狀有關(guān)，在再入過(guò)程中是變化的，因此將其作為其中的一個(gè)狀態(tài)變量進(jìn)行估計(jì)。ρ(h)為空氣密度，與當(dāng)?shù)氐暮０胃叨萮有關(guān)，本文采用美國(guó)標(biāo)準(zhǔn)大氣層模型1976(USSA1976)。

假設(shè)第k時(shí)刻系統(tǒng)狀態(tài)向量為X(k)=[x(k),vx(k),y(k),vy(k),z(k),vz(k),β(k)]T，對(duì)其離散化可得狀態(tài)方程如式(2)所示。

其中F為狀態(tài)轉(zhuǎn)移矩陣，G為輸入控制矩陣，u(k)為加速度控制項(xiàng)，它由式(1)確定，V(k)為各狀態(tài)變量對(duì)應(yīng)的過(guò)程噪聲，這里認(rèn)為它服從零均值的高斯分布，協(xié)方差為 0 .01·[I7×7]。F,G表達(dá)式可參考文獻(xiàn)[9]。

2.2 測(cè)量方程

假設(shè)目標(biāo)被地面雷達(dá)站所捕獲，則獲得的數(shù)據(jù)應(yīng)為天東北坐標(biāo)系(Up-East-North, UEN)的觀測(cè)值Z=[r,a,e]T，其中r為目標(biāo)與雷達(dá)站之間的距離，a和e分別為站心坐標(biāo)下方位角和俯仰角。為了后續(xù)統(tǒng)一處理的方便，將UEN下的測(cè)量值轉(zhuǎn)換到ECEF坐標(biāo)下[6,10]，獲得坐標(biāo)變換后測(cè)量值ZTrans=[x,y,z]T,這樣便能獲得線(xiàn)性的觀測(cè)方程如下：

H為測(cè)量矩陣，W(k+1)為測(cè)量噪聲，不同坐標(biāo)系測(cè)量噪聲的轉(zhuǎn)換過(guò)程可參考文獻(xiàn)[10]，至此完成了狀態(tài)方程和測(cè)量方程的建立。

3 滑窗MCMC問(wèn)題描述及算法實(shí)現(xiàn)

3.1 后驗(yàn)概率模型

以面向測(cè)量數(shù)據(jù)為指導(dǎo)思想，多目標(biāo)數(shù)據(jù)關(guān)聯(lián)問(wèn)題就是將測(cè)量集合進(jìn)行劃分，進(jìn)而獲得后驗(yàn)概率密度最大的組合。

令Z(t)={zi(t):i=1,…,n(t)}為t時(shí)刻雷達(dá)所獲得的目標(biāo)測(cè)量值，則t0時(shí)刻之前所有的測(cè)量集合為= {Z(t) :1 ≤t≤t0}。令Ω為所有可能的關(guān)聯(lián)組合空間，對(duì)于其中任意一個(gè)關(guān)聯(lián)假設(shè)w∈Ω，對(duì)它進(jìn)行如下約束：

(1)w= {τ0,τ1, … ,τk}，其中τk為第k個(gè)子集劃分，即由相應(yīng)的測(cè)量集合構(gòu)成了第k條航跡。因航跡是由測(cè)量值集合經(jīng)過(guò)濾波獲得，若不加特別說(shuō)明本文中航跡與相應(yīng)的測(cè)量集合等價(jià)；

(3)τ0代表一組虛警；

(4)|τi∩z(t) |≤1,i∈ [ 1,k],t∈ [1,t0]；表明每個(gè)測(cè)量點(diǎn)跡只能屬于其中的一條航跡或者屬于虛警；

(5)|τi|≥ 2 ,i∈[1,k]，表示|τi|航跡的長(zhǎng)度至少為2，否則無(wú)法與測(cè)量點(diǎn)跡區(qū)分。

對(duì)于一個(gè)給定的假設(shè)關(guān)聯(lián)w，虛警τ0和真實(shí)航跡組合{τ1,… ,τk}被完全確定，后續(xù)便可對(duì)每條航跡進(jìn)行狀態(tài)估計(jì)。通過(guò)貝葉斯準(zhǔn)則，對(duì)于一個(gè)給定的假設(shè)關(guān)聯(lián)組合w在已知測(cè)量集合Zk的前提下，它的后驗(yàn)概率P(w|Zk)的表達(dá)式如下：

P(Zk)表示測(cè)量集合Zk出現(xiàn)的概率，表達(dá)式與場(chǎng)景具體參數(shù)相關(guān)，如檢測(cè)概率(PD)和監(jiān)視空域內(nèi)虛警目標(biāo)個(gè)數(shù)(假設(shè)它滿(mǎn)足均值為λV的泊松分布)等，具體可參考文獻(xiàn)[11]。P(Zk|w)為關(guān)聯(lián)組合w對(duì)測(cè)量集合Zk的似然函數(shù)，可以通過(guò)動(dòng)力學(xué)方程和測(cè)量模型表達(dá)如下：

式(5)表示關(guān)聯(lián)組合中所有真實(shí)航跡的似然函數(shù)之和，其中C為歸一化常值，|τ|為航跡長(zhǎng)度，并且認(rèn)為航跡是滿(mǎn)足以其濾波值為均值，誤差協(xié)方差為方差的高斯分布。通過(guò)式(4)和式(5)便可得到后驗(yàn)概率P(w|Zk)的解。在給定測(cè)量集Zk的前提下，尋找最大后驗(yàn)概率的關(guān)聯(lián)假設(shè)w等價(jià)如下：

求解式(6)所得的關(guān)聯(lián)組合w即為多目標(biāo)跟蹤中點(diǎn)跡航跡關(guān)聯(lián)的解。

3.2 MCMC原理

Bergman等人[12]于2000年首次將MCMC 算法應(yīng)用于多目標(biāo)跟蹤，Oh Songhwai等人[11,13]進(jìn)一步提出了詳細(xì)的關(guān)聯(lián)假設(shè)步驟，后續(xù)的研究大多集中于靜態(tài)的全局關(guān)聯(lián)實(shí)際應(yīng)用[14,15]，并未對(duì)MCMC關(guān)聯(lián)算法的動(dòng)態(tài)繼承性以及實(shí)時(shí)性進(jìn)一步研究和闡述。

MCMC是一種從分布為π的狀態(tài)空間Ω采樣構(gòu)造一個(gè)馬爾可夫鏈M(其狀態(tài)w，服從平穩(wěn)分布π(w) )的方法。如果當(dāng)前的關(guān)聯(lián)狀態(tài)為w∈Ω，則提出的假設(shè)關(guān)聯(lián)w'∈Ω服從提議分布q(w,w')，其接受概率為A(w,w')由Metropolis Hastings方法獲得。

如果新的關(guān)聯(lián)假設(shè)被否定，則保持上一次關(guān)聯(lián)狀態(tài)w。若M具有不可約性和非周期性，由各態(tài)歷經(jīng)定理可知M收斂于其平穩(wěn)分布π(w)。因此，對(duì)于任意給定的一個(gè)有界函數(shù)f，采樣狀態(tài)經(jīng)過(guò)f作用后的統(tǒng)計(jì)均值可通過(guò)E(f(w) )來(lái)逼近。

MCMC算法的流程如下：

(1)輸入測(cè)量集合Zk，蒙特卡羅仿真次數(shù)mc，任意給出初始關(guān)聯(lián)w=winit;

(2)由提議分布q(w,w')得到一種嘗試關(guān)聯(lián)組合w'；

(3)生成U(0,1)隨機(jī)數(shù)u，若u≤A(w,w')，令w=w'，，否則回到步驟(2)；

(4)直到mc次循環(huán)完成，最終得到的w為最優(yōu)關(guān)聯(lián)組合。

通過(guò)以上算法流程便可得到最優(yōu)的關(guān)聯(lián)組合。從MCMC算法流程看，它累積了足夠的測(cè)量值Zk因此比單步貝葉斯關(guān)聯(lián)算法如 NN關(guān)聯(lián)正確率更高，同時(shí)每次關(guān)聯(lián)只需保存最優(yōu)的關(guān)聯(lián)組合，與多步貝葉斯關(guān)聯(lián)算法MHT相比，能節(jié)約大量存儲(chǔ)空間和運(yùn)行時(shí)間。

提議分布q(w,w')表示兩個(gè)關(guān)聯(lián)組合w和w'通過(guò)單步變化所構(gòu)成的操作組合[11,13]，主要包括新航跡的產(chǎn)生，老航跡的終結(jié)以及航跡更新等等。

3.3 滑窗MCMC算法實(shí)現(xiàn)

傳統(tǒng)的 MCMC算法對(duì)于處理靜態(tài)的測(cè)量集合具有很好的效果，但對(duì)于實(shí)時(shí)的動(dòng)態(tài)變化的測(cè)量集還存在著許多未解決的問(wèn)題。首先是實(shí)時(shí)性問(wèn)題，隨著測(cè)量數(shù)據(jù)的增加，如果采用全局測(cè)量數(shù)據(jù)進(jìn)行組合會(huì)出現(xiàn)“組合爆炸”的情況；其次是關(guān)聯(lián)組合的繼承性和關(guān)聯(lián)假設(shè)概率分配問(wèn)題，一方面要考慮關(guān)聯(lián)的正確率，另一方面要減少不必要的關(guān)聯(lián)假設(shè)的次數(shù)。根據(jù)“分導(dǎo)”這一物理背景對(duì)傳統(tǒng)的MCMC算法進(jìn)行加窗改進(jìn)，包括關(guān)聯(lián)組合繼承性以及關(guān)聯(lián)假設(shè)的概率分配，繼承性改進(jìn)示意如圖1所示。從圖中可見(jiàn)，k時(shí)刻在觀測(cè)掃描窗口[Tk,Tk+Win]內(nèi)測(cè)量集合最優(yōu)組合為w', Win為時(shí)間窗口的大小，下一時(shí)刻時(shí)間窗口滑至[Tk+1,Tk+1+Win]，選擇w2=w'-w1作為下一時(shí)刻馬爾可夫鏈的初始狀態(tài)，這里w1表示在時(shí)間窗[Tk,Tk+1]內(nèi)w'的子集，在剔除w1的同時(shí)將新的測(cè)量值Znew加入滑窗集合，即

圖1 滑窗MCMC算法繼承性示意圖

通過(guò)上述處理，將測(cè)量集合按窗口進(jìn)行劃分，類(lèi)似于降維處理，減小關(guān)聯(lián)組合數(shù)量，同時(shí)將上一窗口的最優(yōu)關(guān)聯(lián)組合遺傳給下一個(gè)窗口作為新的馬爾可夫鏈的初始狀態(tài)，能夠使其快速逼近穩(wěn)態(tài)分布。

另一方面，傳統(tǒng)的 MCMC算法對(duì)于提議分布中的關(guān)聯(lián)假設(shè)是按等概率均勻選取的，而針對(duì)本文的彈頭分導(dǎo)問(wèn)題，很明顯航跡的產(chǎn)生的概率是最大的，而消亡則相對(duì)較小，相應(yīng)地將航跡產(chǎn)生概率調(diào)整為其它假設(shè)的2倍，使它更接近真實(shí)物理場(chǎng)景，快速達(dá)到平穩(wěn)分布。

4 仿真結(jié)果與分析

為驗(yàn)證本文所提出的滑窗 MCMC算法的有效性，設(shè)計(jì)了相應(yīng)的仿真場(chǎng)景，分析了分導(dǎo)前后關(guān)聯(lián)正確概率以及運(yùn)算時(shí)間，并在同等條件下比較了MHT方法效果，驗(yàn)證了本文所提出的算法的優(yōu)越性。

本文所采用的濾波器為 UKF[16](Unscented Kalman Filter)，濾波初始狀態(tài)Xinit可用前兩個(gè)時(shí)刻的測(cè)量值進(jìn)行確定，即

其中x,y,z為初始位置測(cè)量值，Δt為采樣間隔，β0為先驗(yàn)的質(zhì)阻比值，當(dāng)目標(biāo)為近中遠(yuǎn)程彈道導(dǎo)彈它的量級(jí)分別為2000 kg/s2, 5000 kg/s2和8000 kg/s2。初始協(xié)方差求解可參考文獻(xiàn)[10]，滑窗長(zhǎng)度選擇要在關(guān)聯(lián)正確率和關(guān)聯(lián)時(shí)間折中，實(shí)驗(yàn)中發(fā)現(xiàn)窗口長(zhǎng)度大于5個(gè)采樣周期代價(jià)函數(shù)可靠性較好，因此本文滑窗長(zhǎng)度取 5，即統(tǒng)計(jì)結(jié)果從第 5個(gè)采樣點(diǎn)開(kāi)始計(jì)算。

本文所有仿真實(shí)驗(yàn)均在Matlab7.0軟件平臺(tái)，以及CPU為AMD Athlon 3000+，內(nèi)存為DDR400 1 G的硬件平臺(tái)實(shí)現(xiàn)。

4.1 仿真場(chǎng)景設(shè)置

本文模擬某中程彈道導(dǎo)彈軌跡建模，導(dǎo)彈再入高度為70 km，在再入10 s后進(jìn)行彈頭分導(dǎo)，持續(xù)跟蹤 50個(gè)采樣周期。假設(shè)雷達(dá)測(cè)距噪聲標(biāo)準(zhǔn)差為300 m，方位角和俯仰角測(cè)量噪聲標(biāo)準(zhǔn)差均為 0.1 mrad，服從均值為零的正態(tài)分布，采樣間隔Δt=1 s。仿真參數(shù)選取參考美國(guó)太平洋靶場(chǎng)實(shí)驗(yàn)的地基雷達(dá)在跟蹤戰(zhàn)術(shù)彈道導(dǎo)彈的工作參數(shù)[17]。為更真實(shí)地模擬彈頭分導(dǎo)的實(shí)際情況，設(shè)置了3種不同的仿真場(chǎng)景分別為稀疏、密集以及異常密集場(chǎng)景，每個(gè)母彈頭隨機(jī)分出0～2個(gè)子彈頭。多批次彈道航跡的產(chǎn)生過(guò)程如下：以標(biāo)準(zhǔn)彈道再入點(diǎn)的狀態(tài)為中心，分別在3個(gè)方向的位置和速度加上2000 m和100 m/s高斯分布的噪聲作為其余彈道的初始狀態(tài)，并通過(guò)式(1)外推獲得真實(shí)軌跡。分導(dǎo)彈頭位置繼承母彈頭的位置，在速度上以母彈頭速度為基礎(chǔ)，分別受到150 m/s, 75 m/s和75 m/s為標(biāo)準(zhǔn)差的高斯分布零均值的隨機(jī)速度沖量，這樣分導(dǎo)彈頭既能夠保證位置的連續(xù)性，同時(shí)也能夠與母彈頭分開(kāi)，在一定程度上反應(yīng)了分導(dǎo)的物理過(guò)程。

4.2 結(jié)果與分析

因彈頭分導(dǎo)數(shù)目未知，很難對(duì)跟蹤航跡的位置和速度進(jìn)行均方根誤差統(tǒng)計(jì)。本文主要以關(guān)聯(lián)正確率以及計(jì)算時(shí)間來(lái)評(píng)價(jià)算法的優(yōu)劣。關(guān)聯(lián)正確率定義如下：

平均關(guān)聯(lián)正確率則為跟蹤時(shí)間內(nèi)關(guān)聯(lián)正確率的均值，仿真中涉及的參數(shù)有檢測(cè)概率PD，虛警數(shù)目λV以及蒙特卡洛仿真次數(shù)mc。

實(shí)驗(yàn)1(稀疏場(chǎng)景PD=1，λV=0， mc=500)：再入彈道目標(biāo)為10個(gè)，再入10 s后每個(gè)母彈頭隨機(jī)產(chǎn)生0～2的分導(dǎo)彈頭，最終得到17個(gè)彈道目標(biāo)，仿真場(chǎng)景以及全局關(guān)聯(lián)概率如圖2，圖3所示。

從圖3可看出，在該場(chǎng)景下，兩種算法均能得到較高的關(guān)聯(lián)正確概率，大約在92%以上。整體上MHT算法略?xún)?yōu)于滑窗MCMC，因?yàn)樵谠搱?chǎng)景中分導(dǎo)彈頭數(shù)量較少，由分導(dǎo)產(chǎn)生的交叉平行等系統(tǒng)干擾較小，次優(yōu)的滑窗 MCMC只對(duì)掃描窗內(nèi)的數(shù)據(jù)進(jìn)行關(guān)聯(lián)，而 MHT則是進(jìn)行全局關(guān)聯(lián)，充分利用了所有的測(cè)量信息。但 MHT算法波動(dòng)較大，穩(wěn)定性不如滑窗MCMC。MHT和滑窗MCMC計(jì)算時(shí)間分別為35 s與22 s，處于同一量級(jí)，MHT計(jì)算量稍大。由于該場(chǎng)景目標(biāo)較少，分導(dǎo)后對(duì)兩種算法的關(guān)聯(lián)正確率影響并不明顯，仿真結(jié)果也驗(yàn)證了在系統(tǒng)干擾較小的時(shí)候全局最優(yōu)算法具有較好性能。

實(shí)驗(yàn)2(密集場(chǎng)景PD=1,λV=0, mc= 1000)：再入彈道目標(biāo)為30個(gè)，分導(dǎo)參數(shù)設(shè)置不變，最終得到57個(gè)彈道目標(biāo)，仿真場(chǎng)景以及全局關(guān)聯(lián)概率如圖4，圖5所示。

圖2 稀疏仿真場(chǎng)景

圖3 稀疏場(chǎng)景全局關(guān)聯(lián)正確率

圖4 密集仿真場(chǎng)景

從圖4可見(jiàn)，該場(chǎng)景在分導(dǎo)后有較多的交叉現(xiàn)象，場(chǎng)景比較復(fù)雜。從圖5看出該場(chǎng)景下滑窗MCMC算法顯示出明顯的優(yōu)勢(shì)，關(guān)正確率保持在0.9左右，而MHT則在0.85左右，并且前者具有更好的穩(wěn)定性。分導(dǎo)后，兩種算法的關(guān)聯(lián)正確率都有不同程度的下降，MHT更為明顯，正確率下降至0.7左右。分導(dǎo)結(jié)束后，隨著時(shí)間增加，目標(biāo)數(shù)目趨于穩(wěn)定，相互之間相關(guān)性變?nèi)?，兩種算法的關(guān)聯(lián)正確率逐步恢復(fù)至正常水平。在該場(chǎng)景中，MHT和滑窗MCMC運(yùn)行時(shí)間分別為492 s與117 s，滑窗MCMC實(shí)時(shí)性明顯優(yōu)于MHT。這說(shuō)明了最優(yōu)估計(jì)算法雖然能最大程度地利用了測(cè)量信息，但存在以下兩方面不足：第一是計(jì)算量大；第二是對(duì)未知的系統(tǒng)干擾非常敏感。該結(jié)果驗(yàn)證當(dāng)未知系統(tǒng)干擾過(guò)大，最優(yōu)估計(jì)算法性能難于保證。

實(shí)驗(yàn)3(異常密集場(chǎng)景PD=1,λV=0, mc=1000)：再入彈道目標(biāo)為100個(gè)，分導(dǎo)參數(shù)設(shè)置與前面實(shí)驗(yàn)相同，最終得到193個(gè)彈道目標(biāo)，仿真場(chǎng)景以及全局關(guān)聯(lián)概率如圖6，圖7所示。

從圖6可見(jiàn)，該場(chǎng)景目標(biāo)密集程度非常高，特別是在分導(dǎo)后，目標(biāo)數(shù)量幾乎增加了一倍，并且交叉現(xiàn)象非常嚴(yán)重，屬于復(fù)雜場(chǎng)景。分析圖7可得到，在發(fā)生分導(dǎo)后，兩種算法的關(guān)聯(lián)概率都有明顯下降，滑窗MCMC降至0.8左右，而MHT則達(dá)到0.6。隨后滑窗MCMC較快恢復(fù)至正常水平，而MHT正確率則不斷振蕩，經(jīng)過(guò)12 s延遲后才達(dá)到正常水平，穩(wěn)定性比前者相差較大?；?MCMC關(guān)聯(lián)正確率大多保持在 85%以上，而 MHT明顯低于滑窗MCMC。而在該仿真場(chǎng)景下，MHT運(yùn)行時(shí)間到達(dá)8914 s，而滑窗 MCMC為478 s，幾乎相差20倍，實(shí)時(shí)性明顯優(yōu)于MHT。導(dǎo)致兩種算法性能差別的原因與實(shí)驗(yàn)2類(lèi)似，不過(guò)多闡述。

針對(duì)某一確定的場(chǎng)景，選擇不同的檢測(cè)概率和虛警數(shù)目參數(shù)對(duì)算法性能影響不同，采用某一固定的參數(shù)得出的結(jié)論并不全面，因此前3個(gè)實(shí)驗(yàn)均認(rèn)為檢測(cè)概率為1，虛警數(shù)目為0，專(zhuān)門(mén)增加了實(shí)驗(yàn)4，選擇密集場(chǎng)景，針對(duì)不同的檢測(cè)概率和虛警數(shù)目進(jìn)行平均關(guān)聯(lián)正確率的統(tǒng)計(jì)，探討不同檢測(cè)概率和虛警數(shù)目下兩種算法的性能。

實(shí)驗(yàn)4(密集場(chǎng)景)：場(chǎng)景設(shè)置與實(shí)驗(yàn)2完全一致，得到在不同檢測(cè)概率PD和虛警數(shù)目λV下的關(guān)聯(lián)結(jié)果如圖8，圖9所示。

從圖 8可看出，不同檢測(cè)概率下滑窗 MCMC獲得的平均正確關(guān)聯(lián)概率均比 MHT要高，同時(shí)兩種算法的關(guān)聯(lián)正確率都隨著檢測(cè)概率的增加而提高；從圖9可見(jiàn)，隨著虛警數(shù)目的增加，滑窗MCMC與MHT關(guān)聯(lián)效果均下降，MHT下降速度遠(yuǎn)快于滑窗MCMC，當(dāng)虛警數(shù)目大于50, MHT關(guān)聯(lián)正確率小于0.1，而滑窗MCMC還保存在0.4以上，這說(shuō)明了在惡劣環(huán)境下滑窗MCMC效果優(yōu)于MHT。

圖5 密集場(chǎng)景全局關(guān)聯(lián)正確率

圖6 異常密集仿真場(chǎng)景

圖7 異常密集場(chǎng)景全局關(guān)聯(lián)正確率

圖8 不同檢測(cè)概率下的關(guān)聯(lián)結(jié)果

圖9 不同虛警數(shù)目下的關(guān)聯(lián)結(jié)果

5 結(jié)論

本文首先介紹了彈道被動(dòng)段的動(dòng)力學(xué)特性，對(duì)傳統(tǒng)的 MCMC算法的繼承性及關(guān)聯(lián)假設(shè)概率權(quán)值進(jìn)行改進(jìn)，提出一種改進(jìn)的實(shí)時(shí)滑窗 MCMC算法用于解決多批次彈道導(dǎo)彈分導(dǎo)團(tuán)目標(biāo)數(shù)據(jù)關(guān)聯(lián)問(wèn)題。仿真結(jié)果證實(shí)，在惡劣環(huán)境下，本文所提算法不僅關(guān)聯(lián)概率要高于 MHT算法，并且具有較好的穩(wěn)定性，更重要的是本文算法的實(shí)時(shí)性遠(yuǎn)優(yōu)于MHT，適用于密集環(huán)境下的變數(shù)目的實(shí)時(shí)數(shù)據(jù)關(guān)聯(lián)。

[1]Hamouda Y E M and Phillips C. Adaptive sampling for energy-efficient collaborative multi-target tracking in wireless sensor networks [J].IET Wireless Sensor Systems, 2011,1(1): 15-25.

[2]Mansouri M, Snoussi H, and Richard C. Channel estimation and multiple target tracking in wireless sensor networks based on quantised proximity sensors [J].IET Wireless Sensor Systems, 2011, 1(1): 7-14.

[3]Li XR and Bar-Shalom Y. Tracking in clutter with nearest neighbor filters: analysis and performance [J].IEEE Transactions on Aerospace and Electronic Systems, 1996,32(3): 995-1010.

[4]Chang Kuo-chu and Bar-Shalom Y. Joint probabilistic data association for multitarget tracking with possibly unresolved measurements and maneuvers [J].IEEE Transactions on Automatic Control, 1984, 29(7): 585-594.

[5]Evangeline P, Benjamin P, and Michele R. Hybrid algorithms for multitarget tracking using MHT and GM-CPHD [J].IEEE Transactions on Aerospace and Electronic Systems,2011, 47(2): 832-847.

[6]Benavoli A, Chisci L, and Farina A. Tracking of a ballistic missile with a-prior information[J].IEEE Transactions on Aerospace and Electronic Systems, 2007, 43(3): 1000-1016.

[7]Minvielle P. Decades of improvements in re-entry ballistic vehicle tracking[J].IEEE Aerospace and Electronic Systems Magazine, 2005, 20(8): Part 1, CF-1-CF-14.

[8]Li X R and Jilkov V P. Survey of maneuvering target tracking. Part II: motion models of ballistic and space targets[J].IEEE Transactions on Aerospace and Electronic Systems,2010, 46(1): 96-119.

[9]權(quán)太范. 目標(biāo)跟蹤新理論與技術(shù)[M]. 北京: 國(guó)防工業(yè)出版社,2009: 206-213.

Quan Tai-fan. Target Tracking: Advanced Theory and Techniques[M]. Beijing: National Defence Industry Press,2009: 206-213.

[10]何友, 修建娟, 張晶煒, 等. 雷達(dá)數(shù)據(jù)處理及應(yīng)用[M]. 第2版.北京: 電子工業(yè)出版社, 2009: 35-37, 66-74.

He You, Xiu Jian-juan, Zhang Jing-wei,et al.. Radar Data Process and Application[M]. 2nd Edition, Beijing: Electronic Industry Press, 2009: 35-37, 66-74.

[11]Oh Songhwai, Russell S, and Sastry C. Markov chain Monte Carlo data association for general multiple-target tracking problems [C]. 43rd IEEE conference on Decision and Control,Atlantis, Paradise Island, Bahamas, 2004, Vol.1: 735-742.

[12]Bergman N and Doucet A. Markov chain Monte Carlo data association for target tracking [C]. 2000 IEEE International Conference on Acoustics, Speech, and Signal Processing,Istanbul, Turkey, 2000, Vol.2: 705-708.

[13]Oh Songhwai, Russell S, and Sastry S. Markov chain Monte Carlo data association for multi-target tracking [J].IEEE Transactions on Automatic Control, 2009, 54(3): 481-497.

[14]Lian F, Han C Z, Liu W F,et al.. Sequential Monte Carlo implementation and state extraction of the group probability hypothesis density filter for partly unresolvable group targets-tracking problem[J].IET Radar,Sonar&Navigation,2010, 4(5): 685-702.

[15]齊照輝, 劉雪梅, 梁偉. 基于MCMC的導(dǎo)彈主動(dòng)段突防仿真及靈敏度分析 [J]. 系統(tǒng)工程與電子技術(shù), 2010, 32(4):803-806.

Qi Zhao-hui, Liu Xue-mei, and Liang Wei. Survival simulation and sensitivity analysis of ballistic missile in boost phase based on MCMC[J].Journal of Systems Engineering and Electronics, 2010, 32(4): 803-806.

[16]Julier S J and Uhlmann J K. A new extension of the Kalman filter to nonlinear systems [J].Proceedings of SPIE, 1997,3068(182):182-193.

[17]Farrell W J. Interacting multiple model filter for tactical ballistic missile tracking [J].IEEE Transactions on Aerospace and Electronic Systems, 2008, 44(2): 418-426.