儲(chǔ)雪峰，吳 楠，王 鋒，皇甫列鋒

(1.中國人民解放軍戰(zhàn)略支援部隊(duì)信息工程大學(xué)，河南 鄭州 450001;2.中國人民解放軍32682部隊(duì)，山東 濟(jì)南 250000)

天基紅外導(dǎo)彈預(yù)警系統(tǒng)探測彈道導(dǎo)彈，主要是通過探測發(fā)動(dòng)機(jī)尾焰輻射信號(hào)的角度信息來估計(jì)目標(biāo)運(yùn)動(dòng)狀態(tài)，其典型代表是美國的天基紅外系統(tǒng)(Space Based Infrared System，SBIRS)［1］。由于探測器存在系統(tǒng)誤差，且長距離探測受到大氣層較強(qiáng)的干擾，探測數(shù)據(jù)中同時(shí)含有隨機(jī)誤差。為了實(shí)時(shí)估計(jì)目標(biāo)的運(yùn)動(dòng)狀態(tài)，需要對(duì)導(dǎo)彈進(jìn)行運(yùn)動(dòng)建模，并設(shè)計(jì)濾波算法對(duì)導(dǎo)彈運(yùn)動(dòng)進(jìn)行濾波。助推段彈道導(dǎo)彈機(jī)動(dòng)性強(qiáng)，當(dāng)前統(tǒng)計(jì)(Current Statistic，CS)模型［2］將導(dǎo)彈機(jī)動(dòng)加速度看作一階時(shí)間相關(guān)過程，能較好地匹配導(dǎo)彈在助推段的運(yùn)動(dòng)情況。目前，常用的非線性濾波算法主要有擴(kuò)展卡爾曼濾波(Extended Kalman Filtering，EKF)［3］、無跡卡爾曼濾波(Uncented Kalman Filtering，UKF)［4］和容積卡爾 曼 濾 波(Cubature Kalman Filtering，CKF)［5］等。SBIRS探測模型的非線性較強(qiáng)，采用EKF進(jìn)行估計(jì)可能會(huì)產(chǎn)生較大的截?cái)嗾`差［6］。導(dǎo)彈狀態(tài)矢量維數(shù)較高，采用UKF進(jìn)行估計(jì)可能出現(xiàn)狀態(tài)協(xié)方差矩陣不正定的情況，導(dǎo)致濾波中斷［7］。CKF算法是利用容積點(diǎn)非線性傳播的后驗(yàn)統(tǒng)計(jì)特性來近似高斯積分，可以克服上述問題，但時(shí)間更新和測量更新兩個(gè)階段均需構(gòu)造2n個(gè)容積點(diǎn)(n為機(jī)動(dòng)目標(biāo)狀態(tài)矢量的維數(shù))，當(dāng)n較大時(shí)運(yùn)算量較大。由于導(dǎo)彈在助推段運(yùn)動(dòng)時(shí)間短，彈道估計(jì)對(duì)濾波算法時(shí)效性要求比較高，有必要從運(yùn)算效率上對(duì)CKF算法進(jìn)行改進(jìn)。

本文提出一種簡化CS-CKF的助推段彈道估計(jì)方法，測量更新采用構(gòu)造容積點(diǎn)非線性傳播的方式實(shí)現(xiàn)，而狀態(tài)和協(xié)方差一步預(yù)測按照線性狀態(tài)方程計(jì)算。

下面依次介紹基于CS-CKF的彈道估計(jì)方法和基于簡化CS-CKF的彈道估計(jì)方法，通過仿真實(shí)驗(yàn)來驗(yàn)證簡化CS-CKF彈道估計(jì)方法的優(yōu)越性。

1 簡化CS-CKF的彈道估計(jì)方法

1.1 CS模型

用CS模型描述導(dǎo)彈在助推段的運(yùn)動(dòng)狀態(tài)，離散形式的狀態(tài)方程為

當(dāng)探測時(shí)間間隔較小時(shí)，一個(gè)時(shí)間間隔內(nèi)導(dǎo)彈的加速度可認(rèn)為是常數(shù)，即輸入控制矩陣Gk為0。因此，CS模型可以簡化為分段常加速度模型為

過程噪聲矩陣

機(jī)動(dòng)頻率α 為機(jī)動(dòng)時(shí)間常數(shù)的倒數(shù)，假設(shè)導(dǎo)彈在整個(gè)助推段都在機(jī)動(dòng)，則機(jī)動(dòng)時(shí)間常數(shù)就等于助推段飛行時(shí)間。CS模型采用修正的瑞利分布來描述機(jī)動(dòng)加速度均值，機(jī)動(dòng)加速度方差的表達(dá)式為

式中amax為機(jī)動(dòng)加速度的最大值。

1.2 濾波算法初始化

按照文獻(xiàn)［8］提出的方法，求出k=0～2時(shí)刻導(dǎo)彈等效探測位置r0、r1、r2，使用r0、r1、r2對(duì)濾波器進(jìn)行初始化。具體方法如下:

假設(shè)導(dǎo)彈前三個(gè)時(shí)刻做常加速運(yùn)動(dòng)，X2為k=2時(shí)刻的九維狀態(tài)矢量，則k=0～2時(shí)刻的探測方程可寫成

式中，νS為k=0～2時(shí)刻等效位置探測噪聲組成的矢量，相應(yīng)的協(xié)方差矩陣)。令權(quán)重矩陣，利用加權(quán)最小二乘算法計(jì)算k=2時(shí)刻X2的估計(jì)值，則相應(yīng)的狀態(tài)協(xié)方差矩陣，式中的下標(biāo)2均表示k=2時(shí)刻。

1.3 CS-CKF算法

時(shí)間更新階段:

1)構(gòu)造容積點(diǎn)

對(duì)Pk矩陣作Cholesky分解，即，構(gòu)造容積點(diǎn)

2)容積點(diǎn)基于狀態(tài)方程的非線性傳播

其中，Qk為CS模型構(gòu)造的過程噪聲矩陣。

測量更新階段:

4)基于一步狀態(tài)預(yù)測重新構(gòu)造容積點(diǎn)

5)容積點(diǎn)基于探測方程的非線性傳播

6)計(jì)算探測量一步預(yù)測均值

7)計(jì)算濾波增益

其中Rk為探測噪聲矩陣。而后，計(jì)算k+1時(shí)刻濾波器的濾波增益

8)計(jì)算k+1時(shí)刻狀態(tài)估計(jì)值和狀態(tài)誤差協(xié)方差矩陣

1.4 簡化CS-CKF算法

圖1 簡化CS-CKF的遞推處理流程

算法的一個(gè)周期分為兩個(gè)階段。

時(shí)間更新階段:

1)狀態(tài)一步預(yù)測

2)狀態(tài)協(xié)方差一步預(yù)測

測量更新階段與CS-CKF算法相同，不再贅述。

從以上步驟可以看出，單個(gè)周期內(nèi)，簡化CS-CKF算法的時(shí)間更新階段僅需計(jì)算1次，而CS-CKF算法時(shí)間更新階段需計(jì)算2n次。

2 仿真實(shí)驗(yàn)

2.1 仿真流程

仿真流程如圖2所示，主要分為6個(gè)步驟。

1)仿真一條助推段彈道，設(shè)置雙星位置、探測時(shí)間間隔和視線測量誤差，利用雙星探測模型構(gòu)造一組含噪聲的探測數(shù)據(jù)。

2)利用雙星位置，將雙星探測數(shù)據(jù)轉(zhuǎn)化為導(dǎo)彈等效探測位置。

3)利用前3個(gè)時(shí)刻等效探測位置，對(duì)濾波算法進(jìn)行初始化。

4)設(shè)置CS模型相關(guān)參數(shù)，構(gòu)造過程噪聲協(xié)方差矩陣。

5)利用簡化CS-CKF算法進(jìn)行濾波遞推處理，計(jì)算導(dǎo)彈各個(gè)時(shí)刻的估計(jì)值，進(jìn)而得到估計(jì)的彈道數(shù)據(jù)。

6)將估計(jì)的彈道數(shù)據(jù)與仿真的彈道數(shù)據(jù)進(jìn)行對(duì)比分析，驗(yàn)證算法的有效性。

圖2 仿真流程

2.2 仿真場景

以某三級(jí)彈道導(dǎo)彈為例，假設(shè)導(dǎo)彈發(fā)射點(diǎn)的地理坐標(biāo)為160°W、30°N，海拔為40 m，發(fā)射方位角為302.28°，導(dǎo)彈最大負(fù)攻角為6.57°。通過計(jì)算，得到導(dǎo)彈助推段三維彈道曲線，如圖3所示。

分別定點(diǎn)于162°W和110°E的兩顆GEO發(fā)現(xiàn)導(dǎo)彈目標(biāo)，首次探測到導(dǎo)彈的時(shí)間為20 s，探測間隔為1 s，最后探測到導(dǎo)彈的時(shí)間為196 s，視線測量誤差為50 μ rad。

2.3 參數(shù)設(shè)置

CS模型中，導(dǎo)彈機(jī)動(dòng)頻率為0.005 Hz，在地球固聯(lián)坐標(biāo)系中描述，導(dǎo)彈x軸、y軸、z軸方向加速度最大值均為95 m/s2，衛(wèi)星探測時(shí)間間隔為1 s。

2.4 仿真結(jié)果及分析

1)估計(jì)誤差

圖3 導(dǎo)彈助推段三維彈道曲線

對(duì)CS-CKF和簡化CS-CKF算法蒙特卡洛仿真1000次，以x方向?yàn)槔?，狀態(tài)估計(jì)誤差曲線如圖4所示。用實(shí)線描述CS-CKF算法的估計(jì)均方根誤差，用點(diǎn)虛線描述簡化CS-CKF算法的估計(jì)均方根誤差。從圖4可以看出，兩種算法在x方向狀態(tài)估計(jì)均方根誤差差別非常小。y方向、z方向與x方向情況相同，說明兩種算法的狀態(tài)估計(jì)誤差相當(dāng)。

圖4 x方向狀態(tài)估計(jì)誤差曲線

2)一步預(yù)測標(biāo)準(zhǔn)差

對(duì)CS-CKF和簡化CS-CKF算法進(jìn)行單次仿真，對(duì)每個(gè)時(shí)刻一步狀態(tài)協(xié)方差預(yù)測矩陣Pk+1/k的元素Pk+1/k(1，1)、Pk+1/k(4，4)、Pk+1/k(7，7)進(jìn)行開平方，分別得到一步預(yù)測x方向位置、速度、加速度的標(biāo)準(zhǔn)差。對(duì)各個(gè)時(shí)刻兩種算法的一步預(yù)測標(biāo)準(zhǔn)差進(jìn)行相減，得到兩種算法x方向一步預(yù)測標(biāo)準(zhǔn)差的差值情況，如圖5所示，與狀態(tài)估計(jì)值相比，兩種算法的一步預(yù)測標(biāo)準(zhǔn)差差值非常小，可以忽略。其他方向一步預(yù)測標(biāo)準(zhǔn)差的差別情況與x方向相同，由此可知，兩種算法一步預(yù)測的一致性比較強(qiáng)，本文算法是合理的。

圖5 x方向一步預(yù)測標(biāo)準(zhǔn)差的差值

3)運(yùn)算時(shí)間

實(shí)驗(yàn)平臺(tái)計(jì)算機(jī)CPU為Core(TM)i7-4790(3.60 GHz)，內(nèi)存為8 GB，顯卡為AMD Radeon R7350;軟件編程環(huán)境為Matlab2017b。對(duì)CS-CKF算法和簡化CSCKF算法進(jìn)行1000次蒙特卡洛仿真，簡化CS-CKF濾波函數(shù)運(yùn)算時(shí)間為34.76 s，CS-CKF濾波函數(shù)運(yùn)算時(shí)間為48.09 s，簡化CS-CKF濾波函數(shù)運(yùn)算時(shí)間比CS-CKF濾波函數(shù)運(yùn)算時(shí)間縮短38%。

進(jìn)一步對(duì)兩種濾波函數(shù)的代碼進(jìn)行比較，耗費(fèi)時(shí)間相差較大的步驟如表1所示。由表1可以看出，兩種濾波函數(shù)運(yùn)算時(shí)間的主要差別在時(shí)間更新階段。CS-CKF濾波函數(shù)時(shí)間更新階段耗費(fèi)的總時(shí)間為13.19 s，其中構(gòu)造容積點(diǎn)、基于容積點(diǎn)進(jìn)行狀態(tài)一步預(yù)測和協(xié)方差一步預(yù)測耗費(fèi)的時(shí)間分別3.06 s、3.97 s和6.16 s。簡化CS-CKF濾波算法的時(shí)間更新是直接利用線性狀態(tài)方程實(shí)現(xiàn)的，因此運(yùn)算時(shí)間較短。

表1 兩種濾波函數(shù)耗費(fèi)時(shí)間差別較大的步驟運(yùn)算時(shí)間表

3 結(jié)束語

針對(duì)CS-CKF的彈道估計(jì)運(yùn)算量較大的問題，提出一種簡化CS-CKF的助推段彈道估計(jì)方法。依次介紹運(yùn)動(dòng)模型、估計(jì)方法、濾波初始化方法和遞推處理步驟等，并進(jìn)行仿真實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果表明，與傳統(tǒng)CS-CKF的彈道估計(jì)算法相比，本文算法合理有效，在估計(jì)誤差相當(dāng)?shù)臈l件下，運(yùn)算時(shí)間大約縮短38%。對(duì)于實(shí)時(shí)性要求很高的助推段彈道估計(jì)來說，簡化CS-CKF的算法優(yōu)勢比較明顯。