基于CKF的多雷達(dá)分布式再入彈道目標(biāo)實(shí)時(shí)跟蹤算法

李春月，廖育榮，倪淑燕，陳 帥

（航天工程大學(xué) 光電裝備系，北京 101416）

0 引 言

現(xiàn)階段，單雷達(dá)系統(tǒng)對(duì)再入彈道目標(biāo)的實(shí)時(shí)跟蹤精度已不能滿足導(dǎo)彈攔截需求，人們開始研究采用多個(gè)雷達(dá)對(duì)彈道目標(biāo)進(jìn)行協(xié)同跟蹤。多個(gè)雷達(dá)對(duì)彈道導(dǎo)彈的實(shí)時(shí)跟蹤問題本質(zhì)上即為傳感器網(wǎng)絡(luò)[1]中的最優(yōu)狀態(tài)估計(jì)問題。由于傳統(tǒng)的集中式估計(jì)方法帶有一個(gè)數(shù)據(jù)融合中心，多個(gè)測站的量測數(shù)據(jù)均發(fā)往融合中心進(jìn)行計(jì)算與狀態(tài)估計(jì)，大大增加了融合中心的通信量與計(jì)算量。近年來分布式估計(jì)方法受到越來越多的關(guān)注，該方法中每個(gè)雷達(dá)僅與其相鄰雷達(dá)進(jìn)行信息交互，無融合中心，可擴(kuò)展性強(qiáng)，特殊環(huán)境下生存能力強(qiáng)[2-3]。

在解決分布式估計(jì)問題的眾多方法中，文獻(xiàn)[4-6]提出的基于一致性分布式卡爾曼濾波（Distributed Kalman Filter，DKF）算法引起了學(xué)者們的廣泛關(guān)注。該算法利用信息濾波器（Information Kalman Filter，IKF）與平均一致性濾波器級(jí)聯(lián)的方法，實(shí)現(xiàn)了對(duì)離散線性高斯系統(tǒng)中目標(biāo)狀態(tài)的分布式估計(jì)。引入擴(kuò)展卡爾曼濾波（Extended Kalman Filter，EKF）[7]算法可以將 DKF 推廣到高斯非線性系統(tǒng)，但由于EKF算法對(duì)非線性狀態(tài)方程進(jìn)行的是一階泰勒級(jí)數(shù)展開，非線性程度不高，在強(qiáng)非線性系統(tǒng)中濾波精度較低。無跡卡爾曼濾波（Unscented Kalman Filter，UKF）[8]是基于UT變換的采樣策略，采用Sigma點(diǎn)近似表示非線性函數(shù)的分布，有效提高了目標(biāo)狀態(tài)估計(jì)精度，由此基于統(tǒng)計(jì)線性誤差傳播理論的分布式UKF算法也被提出，但是該算法仍未解決UKF算法對(duì)于高維非線性系統(tǒng)（維數(shù)n≥4）濾波不穩(wěn)定的問題。文獻(xiàn)[9-10]通過采用Spherical-Radial規(guī)則近似非線性函數(shù)傳遞的后驗(yàn)均值和協(xié)方差的方法，依據(jù)高斯濾波框架提出容積卡爾曼濾波（Cubature Kalman Filter，CKF）算法。相比于EKF與UKF算法，CKF算法有嚴(yán)格的數(shù)學(xué)推導(dǎo)過程，在沒有顯著增加計(jì)算量的前提下，有效提升了算法的估計(jì)精度與計(jì)算穩(wěn)定性，在各個(gè)領(lǐng)域得到廣泛應(yīng)用[11-12]，因此CKF算法更加適合作為分布式非線性濾波的基礎(chǔ)算法。

本文提出一種基于CKF的多雷達(dá)分布式彈道目標(biāo)跟蹤算法，通過將CKF的信息濾波形式嵌入DKF濾波框架，得到分布式容積信息濾波器。該算法無融合中心，信息交互僅存在于具有通信鏈路的相鄰雷達(dá)之間，每個(gè)雷達(dá)單獨(dú)執(zhí)行濾波計(jì)算，獲得與集中式方法相近的跟蹤精度，仿真結(jié)果驗(yàn)證了算法的有效性。

1 再入彈道目標(biāo)跟蹤數(shù)學(xué)模型

1.1 再入彈道動(dòng)力學(xué)模型

考慮再入彈道目標(biāo)的一般運(yùn)動(dòng)狀態(tài)，對(duì)再入彈道目標(biāo)進(jìn)行受力分析，其主要受到的空氣動(dòng)力表現(xiàn)為大氣阻力，大氣阻力加速度方向與再入速度方向相反。定義以測量雷達(dá)為原點(diǎn)的北天東（North-Up-East）坐標(biāo)系O-xyz，Ox指向正北方向，Oy垂直于地球表面指向正上方，Oz指向正東方向。由于再入彈道導(dǎo)彈速度快，再入時(shí)間短，忽略地球自轉(zhuǎn)角速度的影響。對(duì)再入目標(biāo)建模如圖1所示，假設(shè)地球?yàn)闃?biāo)準(zhǔn)球體，Oe為地球質(zhì)心，ag和af分別為引力加速度和大氣阻力加速度，彈道目標(biāo)的位置、速度分量和彈道系數(shù)構(gòu)成7維狀態(tài)向量，給出再入彈道目標(biāo)的系統(tǒng)狀態(tài)方程：

式中：CD為氣動(dòng)阻力系數(shù)；m為彈道目標(biāo)質(zhì)量；A為彈道目標(biāo)有效截面積。為了保證彈道系數(shù)的非負(fù)性，防止濾波器發(fā)散，采用式（2）的指數(shù)模型，γ(t)的變化率可以用零均值高斯白噪聲表示，即：

圖1 雷達(dá)觀測再入彈道目標(biāo)示意圖Fig.1 Schematic diagram of reentry ballistic target observed by radar

由此可以得到再入彈道目標(biāo)的7維狀態(tài)方程、狀態(tài)量。為了保證濾波器的精度，采用四階龍格-庫塔方法對(duì)狀態(tài)方程進(jìn)行離散，離散的基本模型如下：

離散后的結(jié)果及相應(yīng)參數(shù)為：

最終得到離散非線性系統(tǒng)狀態(tài)方程為：

1.2 多雷達(dá)量測模型

多雷達(dá)量測系統(tǒng)中，假設(shè)共有5個(gè)雷達(dá)測站，分別采用集中式模式與分布式模式。雷達(dá)測站間的通信拓?fù)浣Y(jié)構(gòu)如圖2所示。集中式模式中，雷達(dá)測站僅擔(dān)任量測任務(wù)，量測數(shù)據(jù)統(tǒng)一發(fā)往數(shù)據(jù)融合中心進(jìn)行跟蹤計(jì)算；分布式模式中，各雷達(dá)測站間同時(shí)擔(dān)任量測與目標(biāo)跟蹤任務(wù)，通信鏈路僅存在于相鄰雷達(dá)測站間。

圖2 雷達(dá)間的通信拓?fù)浣Y(jié)構(gòu)示意圖Fig.2 Schematic diagram of communication topology structure among radars

在北天東坐標(biāo)系下建立雷達(dá)觀測再入彈道目標(biāo)的量測模型。如圖1所示，雷達(dá)對(duì)彈道目標(biāo)的距離為r，俯仰角為η，方位角為ε，有z=[r,ε,η]Τ。由觀測值相對(duì)彈道目標(biāo)位置的幾何關(guān)系建立非線性量測方程：

2 多雷達(dá)分布式彈道目標(biāo)跟蹤算法

多雷達(dá)分布式彈道目標(biāo)跟蹤算法本質(zhì)為傳感器網(wǎng)絡(luò)中的目標(biāo)估計(jì)問題，需要通過圖論實(shí)現(xiàn)各雷達(dá)測站間通信拓?fù)涞臄?shù)學(xué)語言描述。結(jié)合平均一致性算法，通過相鄰雷達(dá)間進(jìn)行雙向信息交互，使濾波過程中的某個(gè)中間量實(shí)現(xiàn)全局的一致，最終完成彈道目標(biāo)的分布式估計(jì)。該模式不存在數(shù)據(jù)融合中心，濾波跟蹤算法在各雷達(dá)測站下分別進(jìn)行。首先考慮常規(guī)離散非線性系統(tǒng)：

式中：噪聲ωk，vk相互獨(dú)立，統(tǒng)計(jì)特性分別滿足ωk～N(ωk;0,Qk)，vk～N(vk;0,Rk)。

2.1 基礎(chǔ)理論

用無向圖G=(V,E)對(duì)由n個(gè)雷達(dá)測站的通信拓?fù)浣Y(jié)構(gòu)進(jìn)行數(shù)學(xué)建模，其中，V=[1,2,…,n]表示雷達(dá)測站集合，無向圖G中的通信鏈路由表示，此時(shí)雷達(dá)s與雷達(dá)j連通。雷達(dá)s的連通雷達(dá)集合記為所以有Ns?V。無向圖G中，集合Ns中元素的個(gè)數(shù)記為雷達(dá)s的度，有無向圖G中，如果任意兩雷達(dá)之間存在一條通信路徑（可以為多跳），則稱該無向圖是連通的。

傳感器網(wǎng)絡(luò)中的一致性策略規(guī)定了每個(gè)節(jié)點(diǎn)在僅與鄰居節(jié)點(diǎn)通信條件下更新自身狀態(tài)的過程，其中最直接的方法即為將節(jié)點(diǎn)自身的狀態(tài)更新為自身和鄰居節(jié)點(diǎn)狀態(tài)的加權(quán)線性組合，記為[14]：

式中：ζs(t)為節(jié)點(diǎn)狀態(tài)信息矩陣；μsj(t)為Metropolis加權(quán)系數(shù)[15]，可由式（12）計(jì)算得到。

由式（12）可以看出，加權(quán)系數(shù)僅和節(jié)點(diǎn)自身與其鄰居節(jié)點(diǎn)的度有關(guān)，無需全局信息，基于此可進(jìn)行分布式濾波算法設(shè)計(jì)。

引理 1[16]：假設(shè)ζs(0)為節(jié)點(diǎn)s的初始狀態(tài)，若對(duì)于所有節(jié)點(diǎn)，無向圖G都連通，則根據(jù)式（11）所示的一致性算法，所有節(jié)點(diǎn)的狀態(tài)將收斂于初始值的線性組和，即：

2.2 擴(kuò)展信息濾波器

由于卡爾曼濾波的信息濾波形式可以用信息狀態(tài)向量與信息矩陣分別累加的形式實(shí)現(xiàn)多元信息融合，尤其適合作為分布式算法中的信息融合方法，所以這里首先介紹非線性卡爾曼濾波的信息濾波遞推公式。擴(kuò)展信息濾波器（Extended Information Filter，EIF）實(shí)際上為EKF的等價(jià)信息表示，具體方法為：將EKF中的狀態(tài)向量與誤差協(xié)方差矩陣求逆，得到信息矩陣與信息狀態(tài)向量EIF 算法的主要執(zhí)行步驟為：

1）時(shí)間更新

2）量測更新

3）狀態(tài)更新

上述即為EIF算法的遞推更新過程。

2.3 容積信息濾波器

這里首先介紹CKF算法。CKF算法通過球面-相徑容積準(zhǔn)則選取一組等權(quán)值的容積點(diǎn)來近似狀態(tài)函數(shù)的后驗(yàn)概率分布，對(duì)于狀態(tài)維數(shù)為n的非線性高斯系統(tǒng)，容積點(diǎn)的選取規(guī)則為：

然后進(jìn)行時(shí)間更新過程：

最后進(jìn)行量測更新過程：

計(jì)算k時(shí)刻濾波增益：

完成狀態(tài)向量與協(xié)方差矩陣的更新：

以EIF為框架進(jìn)行容積信息濾波器（Cubature Information Filter，CIF）設(shè)計(jì)。首先要明確EIF中時(shí)間更新與量測更新以逆協(xié)方差矩陣為基礎(chǔ)，并涉及線性觀測矩陣，而CKF算法中狀態(tài)方程與量測方程均為非線性系統(tǒng)，不存在線性觀測矩陣，這時(shí)需要利用線性誤差傳播理論，利用協(xié)方差與互協(xié)方差等效表示EIF中的線性觀測矩陣，即：

將式（33）代入式（17），式（18）得到CIF算法的量測更新：

2.4 DCKF算法實(shí)現(xiàn)

考慮由ns個(gè)雷達(dá)測站組成的再入彈道目標(biāo)實(shí)時(shí)跟蹤系統(tǒng)，第s個(gè)測站的量測方程為：

若采用有數(shù)據(jù)融合中心的集中式濾波方法，所有量測數(shù)據(jù)均由融合中心處理，則式（34），式（35）可以改寫為：

此時(shí)重新定義平均觀測向量與平均逆協(xié)方差矩陣，有：

則式（37），式（38）可化簡為：

結(jié)合一致性策略與引理1觀察式（41），式（42），可以發(fā)現(xiàn)若各個(gè)雷達(dá)測站有相同的初始信息狀態(tài)，則可以通過具有通信鏈路的相鄰雷達(dá)進(jìn)行信息交互，然后結(jié)合平均一致性算法使平均觀測向量與平均逆協(xié)方差矩陣達(dá)到一致，從而以分布式的方式實(shí)現(xiàn)多雷達(dá)集中式彈道目標(biāo)跟蹤算法的近似估計(jì)。進(jìn)而總結(jié)出基于CKF的多雷達(dá)分布式再入彈道目標(biāo)實(shí)時(shí)跟蹤算法執(zhí)行步驟如下（雷達(dá)測站s）：

1）初始化

2）時(shí)間更新

式（21）～式（25）為時(shí)間更新過程，得到一步狀態(tài)估計(jì)與先驗(yàn)誤差協(xié)方差矩陣。

求時(shí)間更新的信息狀態(tài)向量與信息矩陣：

3）量測更新

設(shè)定平均觀測向量與平均逆協(xié)方差矩陣初值：

循環(huán)t=0,1,2,…,τ，τ∈N為達(dá)到一致時(shí)的迭代次數(shù)：

輸出：

4）狀態(tài)更新

通過式（19），式（20）求得更新后的后驗(yàn)狀態(tài)向量與誤差協(xié)方差矩陣。

3 仿真分析

在Matlab R2010b環(huán)境下對(duì)多雷達(dá)分布式觀測再入彈道目標(biāo)系統(tǒng)進(jìn)行建模仿真，對(duì)提出的基于CKF的多雷達(dá)分布式再入彈道目標(biāo)實(shí)時(shí)跟蹤算法進(jìn)行驗(yàn)證。仿真過程中，首先基于狀態(tài)模型，在雷達(dá)站坐標(biāo)系下生成帶有狀態(tài)噪聲的再入彈道目標(biāo)的真實(shí)軌跡，給定再入初值為：

然后用四階龍格-庫塔方法對(duì)狀態(tài)方程進(jìn)行離散時(shí)，離散步長h=1 s。根據(jù)量測模型生成帶有量測噪聲的測量信息用于濾波。狀態(tài)噪聲協(xié)方差矩陣、量測噪聲協(xié)方差矩陣分別為：

生成的再入彈道目標(biāo)真實(shí)軌跡如圖3所示。

圖3 再入彈道目標(biāo)的運(yùn)動(dòng)軌跡Fig.3 Motion trail of reentry ballistic target

設(shè)定蒙特卡洛打靶次數(shù)為500次，分別以單站CKF濾波算法、基于CKF的多站集中式跟蹤算法與基于CKF的多站分布式實(shí)時(shí)跟蹤算法對(duì)上述彈道軌跡進(jìn)行實(shí)時(shí)跟蹤仿真。分布式算法的通信拓?fù)浣Y(jié)構(gòu)如圖2b）所示，采用的是無中心的環(huán)形結(jié)構(gòu)，每個(gè)測站與相鄰的兩個(gè)測站間具有雙向通信鏈路；集中式CKF算法采用的是圖2a）所示的通信拓?fù)浣Y(jié)構(gòu)，數(shù)據(jù)融合中心僅接收量測信息，不參與量測過程。

給出3種算法的濾波初始狀態(tài)向量與誤差協(xié)方差矩陣：

通過統(tǒng)計(jì)各個(gè)時(shí)刻的位置均方根誤差與速度均方根誤差來驗(yàn)證所提算法的有效性，其中位置均方根誤差計(jì)算公式為：

仿真結(jié)果如圖4，圖5所示，給定相同的估計(jì)初值與噪聲矩陣情況下，其中圖4為3種算法的速度均方根誤差，圖5為統(tǒng)計(jì)的3種算法的位置均方根誤差。通過觀察誤差曲線，發(fā)現(xiàn)多雷達(dá)系統(tǒng)跟蹤載入彈道目標(biāo)的收斂速度明顯高于單雷達(dá)系統(tǒng)。其中，集中式算法的實(shí)時(shí)跟蹤精度最高，分布式算法相較于集中式跟蹤精度有稍微下降，但是遠(yuǎn)高于單雷達(dá)測站的測量精度。

圖4 彈道目標(biāo)速度均方根誤差Fig.4 RMSE of ballistic target velocity

圖5 彈道目標(biāo)位置均方根誤差Fig.5 RMSE of ballistic target position

表1統(tǒng)計(jì)了3種算法最后100 s仿真結(jié)果的平均值。從表1中能直觀地看出，基于CKF的分布式再入彈道目標(biāo)實(shí)時(shí)跟蹤算法在降低通信量與計(jì)算量的前提下保持了較高的跟蹤精度，相比于單站CKF算法，速度RMSE降低了0.1 m·s-1，位置RMSE降低了14 m；與多站集中式CKF算法精度相近。

表1 150～250 s速度與位置均方根誤差均值Table 1 Mean of RMSE of position and velocity with in 150～250 s

表2為相同仿真條件下，蒙特卡洛打靶次數(shù)為300次時(shí)，2種算法進(jìn)行250 s再入目標(biāo)跟蹤時(shí)的耗時(shí)統(tǒng)計(jì)，側(cè)面反映出了2種算法的計(jì)算量大小。通過對(duì)比發(fā)現(xiàn)，DCKF算法的計(jì)算用時(shí)較集中式CKF算法有明顯減少，表明對(duì)同一目標(biāo)進(jìn)行相同時(shí)長跟蹤時(shí)，DCKF算法具有顯著降低運(yùn)算量的優(yōu)勢。

表2 跟蹤250 s時(shí)2種算法計(jì)算耗時(shí)比較Table 2 Comparison of time consumption of two algorithms after reentry target tracking for 250 s

4 結(jié) 論

本文提出一種基于CKF的多雷達(dá)分布式彈道目標(biāo)跟蹤算法，通過統(tǒng)計(jì)線性誤差傳播的方法將CKF嵌入擴(kuò)展信息濾波器得到容積信息濾波器，然后利用一致性算法將多量測值集中式濾波進(jìn)行分布式等價(jià)表示，以更新狀態(tài)估計(jì)與誤差協(xié)方差矩陣得到DCKF算法。該算法中，信息交互僅存在于具有通信鏈路的相鄰雷達(dá)之間，每個(gè)雷達(dá)單獨(dú)執(zhí)行濾波計(jì)算，獲得與集中式方法相近的跟蹤精度，仿真結(jié)果驗(yàn)證了本文算法的有效性。