無再入觀測彈道導彈氣動參數和落點聯合預報

吳 楠,王 鋒,孟凡坤

(信息工程大學 數據與目標工程學院,河南 鄭州 450001)

彈道導彈落點預報是導彈預警系統中的一個重要環節,落點預報精度是預警雷達的重要技術指標,提高落點預報和誤差分析精度具有很強的工程意義[1-2]。提高落點預報精度的關鍵是建立準確的預報模型[3],如孫瑜等[4]建立了考慮地球J2攝動的彈道預報模型,潘曉剛等[5]建立了考慮J4攝動的彈道預報模型,但這2種方法均沒有考慮導彈再入時所受氣動阻力的影響。由于再入氣動阻力是影響落點預報精度僅次于地球自轉的因素[6],較多的研究方法是在再入過程中通過引入觀測數據快速地辨識氣動參數,然后再用于落點預報,如陳映等[7]利用交互多模型方法估計彈道系數,史金光等[8]建立了氣動復合系數模型并借助擴展卡爾曼濾波進行辨識,趙捍東等[9]借助觀測數據用徑向基函數神經網絡逼近彈道方程,張金鳳等[10]將再入狀態估計和彈道系數辨識進行聯合優化,均進一步提高了落點預報精度。然而在導彈的早期預警階段,由于導彈并未再入,無法獲得再入觀測數據,且氣動參數信息也無法預先準確獲知,此時如何考慮再入氣動的影響,并建立高精度的落點預報和分析模型是研究的難點。

本文在無再入觀測數據條件下,基于彈道導彈打擊目標時對終端速度變化范圍的限制,反推出氣動阻力參數的變化范圍,然后將氣動阻力參數作為隨機過程建模,并增廣為目標的狀態,建立考慮氣動阻力、地球自轉和J2攝動的高精度落點預報方程,采用數值積分和無跡變換的方法對落點經緯度和落地時刻進行預報和誤差分析。

1 考慮氣動阻力的導彈落點計算方程

(1)

(2)

(3)

(4)

2 氣動阻力參數變化范圍預估

(5)

(6)

(7)

(8)

式中:δβmax(k)為相對于βmax(k)的小量,可取為δβmax(k)=10-8βmax(k)。

3 基于無跡變換的落點均值預報和誤差傳播計算

氣動阻力參數是導彈的固有屬性,具有確定性,但防御方無法獲取該參數準確的變化規律,因此在彈道和落點預報時,只能將氣動力參數作為隨機量處理。假設用Wiener過程描述氣動力參數[11],相應的氣動阻力參數模型可表示為

(9)

式中:w(t)為白噪聲過程。將式(9)離散化,得:

β(k+1)=β(k)+wk

(10)

式中:wk為離散化驅動白噪聲序列,方差取決于隨機序列{βk}的分布模型。由于沒有任何氣動阻力參數的先驗信息,且為了表征氣動力參數的劇烈變化情況,假設氣動阻力參數在變化范圍內服從均勻分布,則氣動阻力參數的均值為

(11)

方差為

(12)

(13)

當系統中含有隨機量時,對其統計特性的分析方法有Monte Carlo抽樣統計法和協方差分析描述函數法等,其中前者具有推導計算簡便、結果直觀可靠等優點,但缺點是耗時過長且結果不是精確解,后者則具有省時、高精度且便于理論分析等優點,但缺點是統計線性化過程過于繁瑣,且要求系統模型全域一階可微,當系統模型較為復雜且不滿足可微性條件時,協方差分析描述函數法將不再適用。無跡變換(unscented transform,UT)是用于計算經過非線性變換的隨機變量統計特性的新方法[12],UT不需要對非線性方程線性化,而是對狀態矢量的概率密度函數進行近似,因此將UT應用于非線性隨機系統的狀態預報,與協方差分析描述函數法相比具有2個優點:①不需要求解Jacobian矩陣,不要求系統模型具有一階可微性,算法更易實現,應用范圍進一步擴大;②預報獲得的后驗均值和方差都能夠精確到二階,具有更高的預報精度。本節采用UT方法,基于再入段增廣狀態方程式(13),對考慮再入氣動阻力的落點進行預報和誤差傳播分析。定義落點信息矢量z=(tfLfBf)T,利用式(13)和式(2)～式(4)可以獲得落點信息矢量和初始增廣狀態矢量間的函數關系H(·):

(14)

(15)

式中:O為零矩陣。利用UT進行均值和誤差傳播計算,首先基于式(15)構造2n+1個Sigma點和相應權重,對于本文情況,y0的維數n=7。

(16)

z(i)=H(χ(i))i=0,1,…,2n

(17)

最后利用統計的方法計算落點信息矢量的均值和協方差矩陣:

(18)

4 落點誤差橢圓計算

(19)

然后將計算得出的誤差橢圓旋轉平移至大地經緯度平面。

(20)

式中:φ∈[0,2π),l為標準差的倍數,通常取l=3(落入該誤差橢圓的概率為97.3%)。

5 仿真分析

5.1 仿真條件設置

基于某雷達針對某型彈道導彈的探測數據,進行濾波和彈道估計,以濾波穩定時刻t0=540 s時的彈道估計值作為預報初值,此時時間為2011年3月7日8時53分19秒(協調世界時),狀態矢量(地固系位置和速度矢量,單位分別為m和m/s)估計均值為

對應的估計誤差協方差矩陣Px0為

表1 導彈的真實落點信息

5.2 未知氣動阻力參數預估

圖1 氣動阻力參數預測均值與真實值

從結果可以看出,再入過程中真實氣動阻力參數是時變的,變化范圍為6.2×10-5～6.8×10-5m2/kg,但變化幅度不是很大,且由于再入過程持續時間較短,所以將其作為常值進行近似處理是較為合理的,而本文基于先驗信息預測獲得的氣動阻力參數均值與真實值相比,處于同一數量級,且相對誤差約為10%,說明本文算法通過先驗信息預測的氣動阻力參數是較為準確而有效的。

5.3 落點預報誤差分析

表2 3種方法落點預報結果比較

從計算結果可以看出:①解析法雖然耗時最短,約為其他2種算法的1/10,但是無論是落點還是時間預報精度都較差;②積分法落點預報精度較高,但由于其時間預報誤差較大,造成時空不匹配,大大增加了落點預報誤差;③本文算法耗時與積分法相當,由于考慮了氣動阻力的影響而且對阻力參數進行了有效的預測,落地時間精度顯著提高,所以時空匹配較好,具有較高的落點和時間聯合預報精度;④本文方法預報時間的標準差要大于另2種方法,這說明當考慮了氣動阻力參數后,再入段速度變化較為劇烈,導致導彈的預報落地時間具有較大的散布,較另2種方法與實際情況更為符合;⑤3種算法計算的誤差橢圓大小和形狀類似,橢圓的長軸(表征縱向標準差)近似沿射向方向,縱向標準差(長軸)顯著大于橫向標準差(短軸),前者約為后者的5倍,這是因為導彈飛行時速度縱向分量顯著大于橫向分量。與前2種方法相比,本文算法計算的誤差橢圓縱向標準差(長軸)有所增加,這是因為考慮氣動阻力的影響后增加了導彈速度縱向分量變化的不確定性,從而增大了落點縱向散布。

以Monte Carlo打靶仿真2 000次獲得的落點散布為基準,分別與3種算法計算的誤差橢圓進行對比,以驗證算法的準確性和有效性,比較結果如圖2所示。

圖2 打靶(2 000次)落點散布與計算誤差橢圓的比較

從結果可以看出:①解析幾何法計算的誤差橢圓精度最差,縱向和橫向都存在偏差,縱向偏差是由于沒有考慮再入氣動阻力導致,因此整體向前偏移(實際落點由于受到阻力作用射程會更近一些),橫向偏差則主要是沒有考慮地球引力J2攝動影響導致;②數值積分法的誤差橢圓精度優于解析幾何法,由于考慮了地球引力J2攝動影響,橫向偏差基本消除,但是由于沒有考慮再入氣動阻力的影響,仍然存在縱向偏差;③本文提出的考慮再入氣動阻力模型法預報精度最高,與實際的打靶結果最為符合,這是因為該方法既考慮了地球引力J2攝動影響,又考慮了再入氣動阻力影響,有效消除了縱向偏差和橫向偏差。

同樣,以Monte Carlo仿真打靶2 000次獲得的落地時間散布為基準,分別與3種算法計算的落地時間均值和±3σ范圍進行比較,如圖3所示,圖中,N為打靶次數。

圖3 打靶(2 000次)落地時間與預測均值方差的比較

結合表2和圖3可以看出,解析法和積分法的結果幾乎相同,由于其預報時間均值顯著小于真值,且標準差結果偏小(與實際不符),導致打靶結果幾乎都超出其預報均值±3σ范圍之外,這說明誤差分析結果與實際統計特性不匹配,這是因為2種方法均沒有考慮氣動阻力的減速作用,導致預報的導彈落地時間均有所提前,從而帶來更大的落點預報誤差;而本文算法由于考慮了再入氣動阻力減速的作用,預報的落地時間與打靶結果符合較好,預報均值與打靶均值較為接近(相差1 s),且打靶結果均落在其預報的均值±3σ范圍以內,說明本文算法的誤差分析結果與實際統計特性更為匹配,可靠度和準確度更高。

6 結論

在彈道導彈早期預警階段,氣動阻力參數對落點預報精度具有顯著影響,但通常無法獲知且無再入觀測數據對其進行辨識。本文基于導彈末端速度變化范圍的限制,有效地預估了氣動阻力參數的變化范圍,然后建立了考慮氣動阻力、地球自轉和J2攝動項的落點預報和誤差分析模型,并計算了落點誤差散布橢圓。仿真結果表明,與傳統的解析法和積分法相比,本文方法由于有效地預估出氣動阻力參數范圍,落點預報精度尤其是預報時間精度顯著提高,計算的誤差散布橢圓與打靶結果最為吻合。對于導彈早期預警階段的彈道和落點預報工作具有較好的應用價值。