吳 楠,王 鋒,孟凡坤
(信息工程大學 數(shù)據(jù)與目標工程學院,河南 鄭州 450001)
彈道導彈落點預報是導彈預警系統(tǒng)中的一個重要環(huán)節(jié),落點預報精度是預警雷達的重要技術指標,提高落點預報和誤差分析精度具有很強的工程意義[1-2]。提高落點預報精度的關鍵是建立準確的預報模型[3],如孫瑜等[4]建立了考慮地球J2攝動的彈道預報模型,潘曉剛等[5]建立了考慮J4攝動的彈道預報模型,但這2種方法均沒有考慮導彈再入時所受氣動阻力的影響。由于再入氣動阻力是影響落點預報精度僅次于地球自轉的因素[6],較多的研究方法是在再入過程中通過引入觀測數(shù)據(jù)快速地辨識氣動參數(shù),然后再用于落點預報,如陳映等[7]利用交互多模型方法估計彈道系數(shù),史金光等[8]建立了氣動復合系數(shù)模型并借助擴展卡爾曼濾波進行辨識,趙捍東等[9]借助觀測數(shù)據(jù)用徑向基函數(shù)神經(jīng)網(wǎng)絡逼近彈道方程,張金鳳等[10]將再入狀態(tài)估計和彈道系數(shù)辨識進行聯(lián)合優(yōu)化,均進一步提高了落點預報精度。然而在導彈的早期預警階段,由于導彈并未再入,無法獲得再入觀測數(shù)據(jù),且氣動參數(shù)信息也無法預先準確獲知,此時如何考慮再入氣動的影響,并建立高精度的落點預報和分析模型是研究的難點。
本文在無再入觀測數(shù)據(jù)條件下,基于彈道導彈打擊目標時對終端速度變化范圍的限制,反推出氣動阻力參數(shù)的變化范圍,然后將氣動阻力參數(shù)作為隨機過程建模,并增廣為目標的狀態(tài),建立考慮氣動阻力、地球自轉和J2攝動的高精度落點預報方程,采用數(shù)值積分和無跡變換的方法對落點經(jīng)緯度和落地時刻進行預報和誤差分析。
(1)
(2)
(3)
(4)
(5)
(6)
(7)
(8)
式中:δβmax(k)為相對于βmax(k)的小量,可取為δβmax(k)=10-8βmax(k)。
氣動阻力參數(shù)是導彈的固有屬性,具有確定性,但防御方無法獲取該參數(shù)準確的變化規(guī)律,因此在彈道和落點預報時,只能將氣動力參數(shù)作為隨機量處理。假設用Wiener過程描述氣動力參數(shù)[11],相應的氣動阻力參數(shù)模型可表示為
(9)
式中:w(t)為白噪聲過程。將式(9)離散化,得:
β(k+1)=β(k)+wk
(10)
式中:wk為離散化驅動白噪聲序列,方差取決于隨機序列{βk}的分布模型。由于沒有任何氣動阻力參數(shù)的先驗信息,且為了表征氣動力參數(shù)的劇烈變化情況,假設氣動阻力參數(shù)在變化范圍內服從均勻分布,則氣動阻力參數(shù)的均值為
(11)
方差為
(12)
(13)
當系統(tǒng)中含有隨機量時,對其統(tǒng)計特性的分析方法有Monte Carlo抽樣統(tǒng)計法和協(xié)方差分析描述函數(shù)法等,其中前者具有推導計算簡便、結果直觀可靠等優(yōu)點,但缺點是耗時過長且結果不是精確解,后者則具有省時、高精度且便于理論分析等優(yōu)點,但缺點是統(tǒng)計線性化過程過于繁瑣,且要求系統(tǒng)模型全域一階可微,當系統(tǒng)模型較為復雜且不滿足可微性條件時,協(xié)方差分析描述函數(shù)法將不再適用。無跡變換(unscented transform,UT)是用于計算經(jīng)過非線性變換的隨機變量統(tǒng)計特性的新方法[12],UT不需要對非線性方程線性化,而是對狀態(tài)矢量的概率密度函數(shù)進行近似,因此將UT應用于非線性隨機系統(tǒng)的狀態(tài)預報,與協(xié)方差分析描述函數(shù)法相比具有2個優(yōu)點:①不需要求解Jacobian矩陣,不要求系統(tǒng)模型具有一階可微性,算法更易實現(xiàn),應用范圍進一步擴大;②預報獲得的后驗均值和方差都能夠精確到二階,具有更高的預報精度。本節(jié)采用UT方法,基于再入段增廣狀態(tài)方程式(13),對考慮再入氣動阻力的落點進行預報和誤差傳播分析。定義落點信息矢量z=(tfLfBf)T,利用式(13)和式(2)~式(4)可以獲得落點信息矢量和初始增廣狀態(tài)矢量間的函數(shù)關系H(·):
(14)
(15)
式中:O為零矩陣。利用UT進行均值和誤差傳播計算,首先基于式(15)構造2n+1個Sigma點和相應權重,對于本文情況,y0的維數(shù)n=7。
(16)
z(i)=H(χ(i))i=0,1,…,2n
(17)
最后利用統(tǒng)計的方法計算落點信息矢量的均值和協(xié)方差矩陣:
(18)
(19)
然后將計算得出的誤差橢圓旋轉平移至大地經(jīng)緯度平面。
(20)
式中:φ∈[0,2π),l為標準差的倍數(shù),通常取l=3(落入該誤差橢圓的概率為97.3%)。
基于某雷達針對某型彈道導彈的探測數(shù)據(jù),進行濾波和彈道估計,以濾波穩(wěn)定時刻t0=540 s時的彈道估計值作為預報初值,此時時間為2011年3月7日8時53分19秒(協(xié)調世界時),狀態(tài)矢量(地固系位置和速度矢量,單位分別為m和m/s)估計均值為
對應的估計誤差協(xié)方差矩陣Px0為
表1 導彈的真實落點信息
圖1 氣動阻力參數(shù)預測均值與真實值
從結果可以看出,再入過程中真實氣動阻力參數(shù)是時變的,變化范圍為6.2×10-5~6.8×10-5m2/kg,但變化幅度不是很大,且由于再入過程持續(xù)時間較短,所以將其作為常值進行近似處理是較為合理的,而本文基于先驗信息預測獲得的氣動阻力參數(shù)均值與真實值相比,處于同一數(shù)量級,且相對誤差約為10%,說明本文算法通過先驗信息預測的氣動阻力參數(shù)是較為準確而有效的。
表2 3種方法落點預報結果比較
從計算結果可以看出:①解析法雖然耗時最短,約為其他2種算法的1/10,但是無論是落點還是時間預報精度都較差;②積分法落點預報精度較高,但由于其時間預報誤差較大,造成時空不匹配,大大增加了落點預報誤差;③本文算法耗時與積分法相當,由于考慮了氣動阻力的影響而且對阻力參數(shù)進行了有效的預測,落地時間精度顯著提高,所以時空匹配較好,具有較高的落點和時間聯(lián)合預報精度;④本文方法預報時間的標準差要大于另2種方法,這說明當考慮了氣動阻力參數(shù)后,再入段速度變化較為劇烈,導致導彈的預報落地時間具有較大的散布,較另2種方法與實際情況更為符合;⑤3種算法計算的誤差橢圓大小和形狀類似,橢圓的長軸(表征縱向標準差)近似沿射向方向,縱向標準差(長軸)顯著大于橫向標準差(短軸),前者約為后者的5倍,這是因為導彈飛行時速度縱向分量顯著大于橫向分量。與前2種方法相比,本文算法計算的誤差橢圓縱向標準差(長軸)有所增加,這是因為考慮氣動阻力的影響后增加了導彈速度縱向分量變化的不確定性,從而增大了落點縱向散布。
以Monte Carlo打靶仿真2 000次獲得的落點散布為基準,分別與3種算法計算的誤差橢圓進行對比,以驗證算法的準確性和有效性,比較結果如圖2所示。
圖2 打靶(2 000次)落點散布與計算誤差橢圓的比較
從結果可以看出:①解析幾何法計算的誤差橢圓精度最差,縱向和橫向都存在偏差,縱向偏差是由于沒有考慮再入氣動阻力導致,因此整體向前偏移(實際落點由于受到阻力作用射程會更近一些),橫向偏差則主要是沒有考慮地球引力J2攝動影響導致;②數(shù)值積分法的誤差橢圓精度優(yōu)于解析幾何法,由于考慮了地球引力J2攝動影響,橫向偏差基本消除,但是由于沒有考慮再入氣動阻力的影響,仍然存在縱向偏差;③本文提出的考慮再入氣動阻力模型法預報精度最高,與實際的打靶結果最為符合,這是因為該方法既考慮了地球引力J2攝動影響,又考慮了再入氣動阻力影響,有效消除了縱向偏差和橫向偏差。
同樣,以Monte Carlo仿真打靶2 000次獲得的落地時間散布為基準,分別與3種算法計算的落地時間均值和±3σ范圍進行比較,如圖3所示,圖中,N為打靶次數(shù)。
圖3 打靶(2 000次)落地時間與預測均值方差的比較
結合表2和圖3可以看出,解析法和積分法的結果幾乎相同,由于其預報時間均值顯著小于真值,且標準差結果偏小(與實際不符),導致打靶結果幾乎都超出其預報均值±3σ范圍之外,這說明誤差分析結果與實際統(tǒng)計特性不匹配,這是因為2種方法均沒有考慮氣動阻力的減速作用,導致預報的導彈落地時間均有所提前,從而帶來更大的落點預報誤差;而本文算法由于考慮了再入氣動阻力減速的作用,預報的落地時間與打靶結果符合較好,預報均值與打靶均值較為接近(相差1 s),且打靶結果均落在其預報的均值±3σ范圍以內,說明本文算法的誤差分析結果與實際統(tǒng)計特性更為匹配,可靠度和準確度更高。
在彈道導彈早期預警階段,氣動阻力參數(shù)對落點預報精度具有顯著影響,但通常無法獲知且無再入觀測數(shù)據(jù)對其進行辨識。本文基于導彈末端速度變化范圍的限制,有效地預估了氣動阻力參數(shù)的變化范圍,然后建立了考慮氣動阻力、地球自轉和J2攝動項的落點預報和誤差分析模型,并計算了落點誤差散布橢圓。仿真結果表明,與傳統(tǒng)的解析法和積分法相比,本文方法由于有效地預估出氣動阻力參數(shù)范圍,落點預報精度尤其是預報時間精度顯著提高,計算的誤差散布橢圓與打靶結果最為吻合。對于導彈早期預警階段的彈道和落點預報工作具有較好的應用價值。