基于滑窗MCMC的彈道導彈分導團目標數據關聯研究

俞建國 劉 梅 王 駿

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

②(中國電子科技集團第29研究所 成都 610036)

1 引言

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

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

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

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

2.1 狀態(tài)方程

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

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

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

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

2.2 測量方程

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

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

3 滑窗MCMC問題描述及算法實現

3.1 后驗概率模型

以面向測量數據為指導思想，多目標數據關聯問題就是將測量集合進行劃分，進而獲得后驗概率密度最大的組合。

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

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

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

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

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

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

P(Zk)表示測量集合Zk出現的概率，表達式與場景具體參數相關，如檢測概率(PD)和監(jiān)視空域內虛警目標個數(假設它滿足均值為λV的泊松分布)等，具體可參考文獻[11]。P(Zk|w)為關聯組合w對測量集合Zk的似然函數，可以通過動力學方程和測量模型表達如下：

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

求解式(6)所得的關聯組合w即為多目標跟蹤中點跡航跡關聯的解。

3.2 MCMC原理

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

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

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

MCMC算法的流程如下：

(1)輸入測量集合Zk，蒙特卡羅仿真次數mc，任意給出初始關聯w=winit;

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

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

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

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

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

3.3 滑窗MCMC算法實現

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

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

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

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

4 仿真結果與分析

為驗證本文所提出的滑窗 MCMC算法的有效性，設計了相應的仿真場景，分析了分導前后關聯正確概率以及運算時間，并在同等條件下比較了MHT方法效果，驗證了本文所提出的算法的優(yōu)越性。

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

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

本文所有仿真實驗均在Matlab7.0軟件平臺，以及CPU為AMD Athlon 3000+，內存為DDR400 1 G的硬件平臺實現。

4.1 仿真場景設置

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

4.2 結果與分析

因彈頭分導數目未知，很難對跟蹤航跡的位置和速度進行均方根誤差統計。本文主要以關聯正確率以及計算時間來評價算法的優(yōu)劣。關聯正確率定義如下：

平均關聯正確率則為跟蹤時間內關聯正確率的均值，仿真中涉及的參數有檢測概率PD，虛警數目λV以及蒙特卡洛仿真次數mc。

實驗1(稀疏場景PD=1，λV=0， mc=500)：再入彈道目標為10個，再入10 s后每個母彈頭隨機產生0～2的分導彈頭，最終得到17個彈道目標，仿真場景以及全局關聯概率如圖2，圖3所示。

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

實驗2(密集場景PD=1,λV=0, mc= 1000)：再入彈道目標為30個，分導參數設置不變，最終得到57個彈道目標，仿真場景以及全局關聯概率如圖4，圖5所示。

圖2 稀疏仿真場景

圖3 稀疏場景全局關聯正確率

圖4 密集仿真場景

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

實驗3(異常密集場景PD=1,λV=0, mc=1000)：再入彈道目標為100個，分導參數設置與前面實驗相同，最終得到193個彈道目標，仿真場景以及全局關聯概率如圖6，圖7所示。

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

針對某一確定的場景，選擇不同的檢測概率和虛警數目參數對算法性能影響不同，采用某一固定的參數得出的結論并不全面，因此前3個實驗均認為檢測概率為1，虛警數目為0，專門增加了實驗4，選擇密集場景，針對不同的檢測概率和虛警數目進行平均關聯正確率的統計，探討不同檢測概率和虛警數目下兩種算法的性能。

實驗4(密集場景)：場景設置與實驗2完全一致，得到在不同檢測概率PD和虛警數目λV下的關聯結果如圖8，圖9所示。

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

圖5 密集場景全局關聯正確率

圖6 異常密集仿真場景

圖7 異常密集場景全局關聯正確率

圖8 不同檢測概率下的關聯結果

圖9 不同虛警數目下的關聯結果

5 結論

本文首先介紹了彈道被動段的動力學特性，對傳統的 MCMC算法的繼承性及關聯假設概率權值進行改進，提出一種改進的實時滑窗 MCMC算法用于解決多批次彈道導彈分導團目標數據關聯問題。仿真結果證實，在惡劣環(huán)境下，本文所提算法不僅關聯概率要高于 MHT算法，并且具有較好的穩(wěn)定性，更重要的是本文算法的實時性遠優(yōu)于MHT，適用于密集環(huá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]權太范. 目標跟蹤新理論與技術[M]. 北京: 國防工業(yè)出版社,2009: 206-213.

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

[10]何友, 修建娟, 張晶煒, 等. 雷達數據處理及應用[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的導彈主動段突防仿真及靈敏度分析 [J]. 系統工程與電子技術, 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.