儲雪峰，吳 楠，王 鋒，皇甫列鋒

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

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

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

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

1 簡化CS-CKF的彈道估計方法

1.1 CS模型

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

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

過程噪聲矩陣

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

式中amax為機動加速度的最大值。

1.2 濾波算法初始化

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

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

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

1.3 CS-CKF算法

時間更新階段:

1)構(gòu)造容積點

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

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

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

測量更新階段:

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

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

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

7)計算濾波增益

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

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

1.4 簡化CS-CKF算法

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

算法的一個周期分為兩個階段。

時間更新階段:

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

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

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

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

2 仿真實驗

2.1 仿真流程

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

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

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

3)利用前3個時刻等效探測位置，對濾波算法進行初始化。

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

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

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

圖2 仿真流程

2.2 仿真場景

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

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

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

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

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

1)估計誤差

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

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

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

2)一步預(yù)測標準差

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

圖5 x方向一步預(yù)測標準差的差值

3)運算時間

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

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

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

3 結(jié)束語

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