基于時(shí)變TVAR模型和CKF濾波的助推器落點(diǎn)預(yù)測(cè)

朱紫陌，陳 龍，魏昌全，李 黎

（1.92941部隊(duì)；2.92419部隊(duì)，遼寧葫蘆島125000）

助推器的落點(diǎn)預(yù)測(cè)對(duì)其快速準(zhǔn)確的回收具有重要意義。目前，計(jì)算質(zhì)點(diǎn)落點(diǎn)的主要方法有數(shù)值積分法、濾波外推法和線性化法等。文獻(xiàn)[1]用彈道方程解算法和系數(shù)解算法，因彈道方程是在一定假設(shè)條件下的簡(jiǎn)化處理模型。質(zhì)點(diǎn)在空氣中受到各種環(huán)境因素的影響，彈道方程只是反映系統(tǒng)內(nèi)部變化的規(guī)律，與實(shí)際模型存在較大差距。文獻(xiàn)[2]提出采用正交多項(xiàng)式擬合彈道和誤差補(bǔ)償預(yù)報(bào)彈道落點(diǎn)的方法，但對(duì)超出離散點(diǎn)范圍的點(diǎn)進(jìn)行預(yù)報(bào)時(shí)，泛化能力有限。濾波外推法能降低噪聲和隨機(jī)誤差對(duì)落點(diǎn)預(yù)測(cè)的影響，但易于發(fā)散，達(dá)不到準(zhǔn)確落點(diǎn)預(yù)測(cè)。

本文從目標(biāo)運(yùn)動(dòng)軌跡的幾何本質(zhì)出發(fā)，利用時(shí)間序列的分析方法建立彈道運(yùn)動(dòng)軌跡。而助推器脫落后運(yùn)動(dòng)軌跡都屬于非平穩(wěn)時(shí)間序列信號(hào)，時(shí)變自回歸（TVAR）模型因引入時(shí)變參數(shù)使得對(duì)非平穩(wěn)序列信號(hào)合理建模，省略了傳統(tǒng)的非平穩(wěn)時(shí)間序列建模中常用的平穩(wěn)化過程，提高模型置信度。本文建立彈道測(cè)量空間位置信息的TVAR 模型，將時(shí)變系數(shù)看作互相關(guān)的多維時(shí)間序列，預(yù)測(cè)助推器未來位置。

1 時(shí)變自回歸（TVAR）模型

自回歸（AR）模型在描述動(dòng)態(tài)數(shù)據(jù)方面有其獨(dú)特的優(yōu)勢(shì)，在所需預(yù)測(cè)周期較短的應(yīng)用情景中，相比自回歸滑動(dòng)平均（ARMA）模型，其參數(shù)估算比較簡(jiǎn)單、精度更高。在隨機(jī)序列中，特定時(shí)間窗口內(nèi)隨機(jī)變量及白噪聲平均值可以決定下一時(shí)間點(diǎn)的隨機(jī)變量。對(duì)于非平穩(wěn)隨機(jī)時(shí)序則利用觀測(cè)數(shù)據(jù)形成的觀測(cè)向量建立TVAR模型。用TVAR模型建模需對(duì)時(shí)變系數(shù)和模型階數(shù)進(jìn)行估計(jì)。

設(shè)0均值、p階TVAR模型如下：

式（1）中：t=1,2,…,N；?t為該模型時(shí)變系數(shù)；p為該模型階數(shù)；εt為均值為0 的白噪聲；Xt為隨機(jī)序列變量。

Xt只與Xt-1，Xt-2，…，Xt-p有線性關(guān)系，如果εt服從均值為0且方差為Q的高斯分布，則Xt是服從正態(tài)分布的隨機(jī)變量。

1.1 階數(shù)的估計(jì)

TVAR模型階數(shù)的選擇關(guān)系到計(jì)算量、實(shí)時(shí)性、精度等問題。因此，須要選取最適合應(yīng)用場(chǎng)景的模型階數(shù)。確定模型階數(shù)的方法有最佳準(zhǔn)則函數(shù)法、殘差方差圖定階法、F檢驗(yàn)定階法等，常使用的是最佳準(zhǔn)則函數(shù)法[4-6]。最佳準(zhǔn)則函數(shù)法包括AIC 準(zhǔn)則、BIC 準(zhǔn)則等。

AIC 準(zhǔn)則的一般形式為：

AIC=-2-ln(模型極大似然度)+2(參數(shù)個(gè)數(shù))。

1.2 時(shí)變系數(shù)的估計(jì)

2 容積卡爾曼濾波彈道模型

本文采用質(zhì)點(diǎn)彈道模型作為外推彈道的狀態(tài)方程[8-10]：

3 助推器落點(diǎn)預(yù)測(cè)

助推器在脫落之后，已經(jīng)無法獲取彈載衛(wèi)星系統(tǒng)轉(zhuǎn)發(fā)器轉(zhuǎn)發(fā)的測(cè)量信息。為了解決未知的需要預(yù)測(cè)的助推器落地之前短暫時(shí)間內(nèi)的位置狀態(tài)集合，本文首次提出一種利用TVAR模型對(duì)測(cè)量值組成的時(shí)間序列進(jìn)行自回歸模型建立的方法，如式（1）所示，經(jīng)過n步預(yù)測(cè)，重復(fù)式（3）、（4），以助推器高度落到以落點(diǎn)為原點(diǎn)的相對(duì)坐標(biāo)系中高度為0 時(shí)作為結(jié)束時(shí)刻，得出一系列預(yù)測(cè)的外彈道偽測(cè)量值z(mì)k～zk+n，從而將濾波預(yù)測(cè)問題轉(zhuǎn)化為估計(jì)問題。

基于CKF 的落點(diǎn)( )x,y,z預(yù)測(cè)算法的具體步驟如下。

步驟1：處理外彈道數(shù)據(jù)，按照節(jié)1 中式（2）～（10）定階較優(yōu)的TVAR模型；

步驟2：進(jìn)行AIC、BIC準(zhǔn)則模型適用性檢驗(yàn)；

步驟3：將檢驗(yàn)后的TVAR模型進(jìn)行下一步預(yù)測(cè)，當(dāng)預(yù)測(cè)值接近設(shè)定的閾值時(shí)，測(cè)量值預(yù)測(cè)結(jié)束；

步驟4：以彈載衛(wèi)星系統(tǒng)最后一刻的轉(zhuǎn)發(fā)測(cè)量值為初始量，進(jìn)行離散化的彈道質(zhì)點(diǎn)模型的狀態(tài)時(shí)間更新，使用TVAR 模型偽測(cè)量值進(jìn)行量測(cè)更新到閾值時(shí)刻的狀態(tài)量(x,y,z)便是落點(diǎn)預(yù)測(cè)值。

4 應(yīng)用驗(yàn)證

為了驗(yàn)證上述算法的有效性，使用一次外場(chǎng)實(shí)測(cè)數(shù)據(jù)進(jìn)行算法驗(yàn)證。選取一組彈載衛(wèi)星系統(tǒng)外測(cè)彈道實(shí)驗(yàn)數(shù)據(jù)，即通過衛(wèi)星系統(tǒng)信標(biāo)差分計(jì)算出時(shí)間序列下的地理坐標(biāo)測(cè)量值。通過坐標(biāo)轉(zhuǎn)換，把它們轉(zhuǎn)換為相對(duì)理論落點(diǎn)為平面的相對(duì)直角坐標(biāo)系(xi,yi,zi)。如圖1 所示，取最后298 個(gè)衛(wèi)星系統(tǒng)彈道測(cè)量點(diǎn)的位置信息進(jìn)行TVAR模型建立。

圖1 外彈道軌跡圖Fig.1 Outer ballistic trajectory

利用Matlab和Eviews軟件進(jìn)行分析，對(duì)外彈道數(shù)據(jù)進(jìn)行平穩(wěn)性檢驗(yàn)。時(shí)間序列X(t)的樣本自相關(guān)系數(shù)如圖2所示。圖3所示X(t)序列的單位根檢驗(yàn)證明了外彈道數(shù)據(jù)是非平穩(wěn)的。

圖2 X(t)序列樣本的自相關(guān)系數(shù)和樣本偏相關(guān)系數(shù)Fig.2 X(t)series autocorrelation coefficient and partial correlation coefficient

圖3 X(t)序列樣本的ADF單位根檢驗(yàn)Fig.3 Unit Root test statistic of ADF for X(t)

為了建立TVAR(p)模型，p擬定選取30，即從1～30建立30個(gè)參數(shù)。利用AIC和BIC準(zhǔn)則，由圖4～6可看出，對(duì)X、Z、Y通道數(shù)據(jù)，p值較好的是10、10、9，但是考慮時(shí)變參數(shù)m次方的選擇，最后p值確定為6。經(jīng)式（7）基函數(shù)正交最小二乘方法計(jì)算，建立X軸、Y軸和Z軸的TVAR模型。

地球曲率造成一定軌跡高度以下的外彈道數(shù)據(jù)無法通過外測(cè)手段進(jìn)行獲得，所以無法確定相對(duì)于預(yù)示落點(diǎn)的落點(diǎn)坐標(biāo)。根據(jù)時(shí)間序列的自回歸擬合模型，可以通過設(shè)置高度為0 來計(jì)算終止閾值進(jìn)行軌跡預(yù)測(cè)，得到X、Y、Z3個(gè)方向的外彈道偽測(cè)量值。取彈載衛(wèi)星系統(tǒng)彈道均勻時(shí)間步長(zhǎng)測(cè)量的298個(gè)位置測(cè)量信息，與其TVAR 模型相比，由圖7～9 可以看出，兩者的擬合程度非常接近，但是由于外彈道數(shù)據(jù)自身存在噪聲，導(dǎo)致預(yù)測(cè)得到的外彈道偽測(cè)量值也存在較大的噪聲波動(dòng)?；讷@取的外彈道數(shù)據(jù)偽測(cè)量值作為CKF濾波的量測(cè)信息，分別利用質(zhì)點(diǎn)理論模型建立狀態(tài)方程、對(duì)偽測(cè)量值建立量測(cè)方程并估計(jì)其測(cè)量誤差矩陣，通過不斷遞推狀態(tài)方程和量測(cè)方程更新獲取對(duì)X、Y、Z3 個(gè)方向的濾波值。圖7～9 表明，濾波后軌跡噪聲顯著減小。

圖4 X 坐標(biāo)的TVAR模型的AIC、BIC階數(shù)Fig.4 AIC and BIC order of X coordinate TVAR model

圖5 Z 坐標(biāo)的TVAR模型的AIC、BIC階數(shù)Fig.5 AIC and BIC order of Z coordinate TVAR model

圖6 Y 坐標(biāo)的TVAR模型的AIC、BIC階數(shù)Fig.6 AIC and BIC order of Y coordinate TVAR model

圖7 X 坐標(biāo)的TVAR模型和觀測(cè)值曲線Fig.7 X coordinate TVAR model and observed value curve

圖8 Z 坐標(biāo)的TVAR模型和觀測(cè)值曲線Fig.8 Z coordinate TVAR model and observed value curve

圖9 Y 坐標(biāo)的TVAR模型和觀測(cè)值曲線Fig.9 Y coordinate TVAR model and observed value curve

依照上述方法進(jìn)行了20 次的CKF 濾波試驗(yàn)，試驗(yàn)統(tǒng)計(jì)落點(diǎn)的坐標(biāo)偏差分布。統(tǒng)計(jì)表明，TVAR 和CKF 聯(lián)合預(yù)測(cè)的落點(diǎn)預(yù)測(cè)精度略高，平均標(biāo)準(zhǔn)差為TVAR預(yù)測(cè)的74.54%，如圖10所示。

圖10 TVAR預(yù)測(cè)與TVAR和CKF聯(lián)合預(yù)測(cè)落點(diǎn)分布的比較Fig.10 Comparison of impact point between TVAR and TVAR-CKF

5 結(jié)論

時(shí)變TVAR 能較好地用于非平穩(wěn)時(shí)序的建模，CKF對(duì)非線性系統(tǒng)濾波性能好、精度高。該方法結(jié)合時(shí)變TVAR 和CKF 兩者的優(yōu)點(diǎn)對(duì)助推器落點(diǎn)進(jìn)行預(yù)測(cè)。實(shí)驗(yàn)結(jié)果表明，雖然TVAR 具有較強(qiáng)的擬合數(shù)據(jù)能力，也可以短時(shí)預(yù)測(cè)最后的落點(diǎn)，但單純使用TVAR預(yù)測(cè)的彈落點(diǎn)沒有和融合CKF 濾波估計(jì)的助推器落點(diǎn)更接近實(shí)際落點(diǎn)，這為助推器落點(diǎn)預(yù)測(cè)提供了一種可行的方法。