孫　瑜,吳　楠,孟凡坤,鄧啟正

(1.中國人民解放軍96669部隊,北京 102208;2.信息工程大學(xué) 導(dǎo)航與空天目標工程學(xué)院,鄭州 450001)

?

考慮J2攝動的彈道導(dǎo)彈高精度彈道預(yù)報和誤差傳播分析

孫瑜1,2,吳楠2,孟凡坤2,鄧啟正1

(1.中國人民解放軍96669部隊,北京 102208;2.信息工程大學(xué) 導(dǎo)航與空天目標工程學(xué)院,鄭州 450001)

摘要:為進行高精度的彈道導(dǎo)彈彈道預(yù)報和誤差傳播分析,分析了平根數(shù)法解攝動方程的解形式,提出了一種考慮J2攝動并包含短周期項的解析法,并將其用于彈道預(yù)報,結(jié)合無跡變換構(gòu)造了彈道預(yù)報的誤差傳播分析方法,解決了解析法中開普勒根數(shù)求解過程不可微、偏導(dǎo)數(shù)矩陣難以計算的問題。仿真結(jié)果表明,將解析法用于彈道預(yù)報時,J2攝動短周期項導(dǎo)致的誤差達3個數(shù)量級,不可忽略;該算法預(yù)報的相對誤差只有不到5%,精度與高精度數(shù)值法相當,遠高于傳統(tǒng)解析法,但運算速度是后者的7倍,具有更好的效費比。

關(guān)鍵詞:彈道預(yù)報;無跡變換;J2攝動;解析法;短周期項

彈道預(yù)報與誤差傳播分析是戰(zhàn)略預(yù)警系統(tǒng)的重要研究內(nèi)容。傳統(tǒng)的彈道預(yù)報方法通常有數(shù)值法和解析法兩類。數(shù)值法是在充分考慮導(dǎo)彈飛行過程中各種攝動因素的基礎(chǔ)上,建立導(dǎo)彈動力學(xué)方程,并采用一定的積分準則,利用數(shù)值逼近的方法計算彈道。文獻[1～2]基于導(dǎo)彈動力學(xué)方程采用不同的數(shù)值積分方法計算彈道,驗證了數(shù)值法的有效性。然而,這種方法預(yù)報的精確程度及計算速度與積分步長的選取有很大關(guān)系。小步長積分可以保證較高的預(yù)報精度,但是計算速度較慢;大步長積分計算速度快,但是預(yù)報精度不高。解析法是將導(dǎo)彈自由段飛行軌跡近似為橢圓彈道,通過計算開普勒根數(shù)進行彈道預(yù)報[3]。文獻[4～7]基于橢圓彈道理論,采用解析法實現(xiàn)彈道預(yù)報,并分別從時間和預(yù)報精度等方面提出優(yōu)化算法。盡管采用了優(yōu)化算法,但解析法的精度仍低于數(shù)值法,這是因為解析法將地球視作質(zhì)點,沒有考慮地球是形狀不規(guī)則的扁球體。文獻[8～10]用考慮J2攝動的解析法計算人造衛(wèi)星的軌道,因為人造衛(wèi)星的運行周期很長,所以在解攝動方程的時候為了計算簡便,通常忽略短周期項。然而導(dǎo)彈與人造衛(wèi)星截然不同的是導(dǎo)彈的飛行時間很短,它的短周期項能否忽略值得討論。

本文提出了一種考慮J2攝動并包含短周期項的解析法,并將其用于彈道預(yù)報,與數(shù)值法彈道預(yù)報進行仿真對比預(yù)報的精度及耗時。同時,在誤差傳播分析過程中,結(jié)合無跡變換(Unscented Transformation,UT)構(gòu)造了彈道預(yù)報的誤差傳播分析方法,解決了解析法中開普勒根數(shù)求解過程不可微、偏導(dǎo)數(shù)矩陣難以計算的問題。

1考慮J2攝動并包含短周期項的彈道預(yù)報解析法

彈道導(dǎo)彈關(guān)機后進入自由段依靠慣性飛行,而自由段的高度一般在200 km以上,此時空氣阻力的影響可以忽略不計,導(dǎo)彈主要外力源為地球中心引力,可以將導(dǎo)彈繞地球的運動視為二體運動,從理論上可以認為其彈道為橢圓彈道,落點為橢圓彈道與地球的交點。解析法通過計算橢圓彈道的6個開普勒根數(shù)來預(yù)報導(dǎo)彈任意時刻的狀態(tài)。然而地球作為一個形狀不規(guī)則的扁球體,其重力場對自由段的影響很大,二階帶諧系數(shù)J2對彈道導(dǎo)彈自由段落點偏差的影響最大可達十幾km,四階帶諧系數(shù)J4對自由段落點偏差的影響可達數(shù)十m[11-12]。文獻[13]采用數(shù)值積分法，分析了不同積分步長、不同緯度和發(fā)射方位角下J2和J4項以及不同階數(shù)的擾動重力對導(dǎo)彈自由段落點的影響。J2項導(dǎo)致的攝動加速度約為10-3量級,而J4項攝動與大氣阻力、第三體引力、太陽光壓等其他攝動均為10-6量級及以下。若考慮J4項攝動,則大氣阻力、三體引力、太陽光壓等攝動均不能忽略,此時將大幅增加計算耗時。因此,綜合考慮預(yù)報的精度和耗時,本文只考慮J2項攝動的影響。

1.1J2攝動運動方程

地球是一個形狀不規(guī)則的扁球體,要精確地積分計算出其勢函數(shù),必須知道地球的表面形狀及內(nèi)部的密度分布。目前,通常采用球函數(shù)展開式推導(dǎo)出地球勢函數(shù)的標準表達式[14]:

(Cnmcosmλ+Snmsinmλ)Pnm(sinφ)

(1)

式中:f為萬有引力常數(shù);me為地球質(zhì)量;r為地球外地心距;ae為地球赤道平均半徑;Jn為帶諧系數(shù),Jn=-Cn0;n≠m時,Cnm,Snm稱為田諧系數(shù);n=m時,Cnm,Snm稱為扇諧系數(shù);Pn(sinφ)稱為勒讓德函數(shù);Pnm(sinφ)稱為締合勒讓德函數(shù);φ,λ為地球的地心緯度和經(jīng)度。

物體繞地球運動的微分方程可寫成:

(2)

(Cnmcosmλ+Snmsinmλ)Pnm(sinφ)

(3)

當n=2時,即為考慮J2項攝動,取J2=1.082 63×10-3。

以6個開普勒根數(shù)為基本變量的拉格朗日型攝動運動方程[15]為

(4)

式中:a為軌道半長軸;e為偏心率;i為軌道傾角;Ω為升交點赤經(jīng);A為近地點幅角;M為平近點角;ω為平均角速度。設(shè)σ=(aeiΩAM)為開普勒根數(shù)構(gòu)成的向量。

角析法中,平均根數(shù)法解攝動運動方程的特征是將攝動變化中不同性質(zhì)的項區(qū)分開,相應(yīng)的冪級數(shù)解的形式為

(5)

式中:

1.2開普勒根數(shù)轉(zhuǎn)移

對于衛(wèi)星而言,由于其在軌運行時間長,在計算其攝動運動方程時,為了計算簡便通常忽略短周期項。對于彈道導(dǎo)彈,其飛行時間很短,所以在計算彈道導(dǎo)彈的攝動運動方程的時候不能忽略短周期項。而對于具體取到多少階的攝動解,則要取決于預(yù)報的精度需求。文獻[13]分析了不同階數(shù)對彈道導(dǎo)彈自由段落點的影響,中低階部分的影響從幾十m到數(shù)百m不等,對于自由段射程約6 000km的彈道導(dǎo)彈,考慮到30階左右可以達到m級的計算精度。本文重點在于綜合考慮預(yù)報的精度和預(yù)警時間,因此本文只考慮一階攝動解:包括所有的一階攝動項(長期項、長周期項和短周期項)和二階長期項。

考慮J2項攝動的平均根數(shù)法一階解完整結(jié)果如下:

(6)

2基于無跡變換的誤差傳播分析

計算協(xié)方差矩陣時,純粹的解析法需要計算笛卡爾坐標對開普勒根數(shù)的偏導(dǎo)數(shù)矩陣以及開普勒根數(shù)的轉(zhuǎn)移矩陣,然而在計算開普勒根數(shù)轉(zhuǎn)移矩陣時會出現(xiàn)一個問題:如果考慮短周期項,則開普勒根數(shù)求解過程不可微,偏導(dǎo)數(shù)矩陣難以計算;如果不考慮短周期項,計算轉(zhuǎn)移矩陣所用的模型就與狀態(tài)方程中所用的模型不一致,得到的解是完全錯誤的,這個結(jié)論在文獻[16]中得到證實。而UT用采樣點來近似狀態(tài)向量的后驗概率密度,結(jié)合UT的解析法在對彈道導(dǎo)彈進行誤差傳播分析時,就不需要計算笛卡爾坐標對開普勒根數(shù)的偏導(dǎo)數(shù)矩陣以及開普勒根數(shù)的轉(zhuǎn)移矩陣。

(7)

(8)

通過橢圓彈道理論計算每個χ點對應(yīng)的開普勒根數(shù)σ,然后考慮J2攝動采用平根數(shù)法解拉格朗日攝動運動方程,計算出t時刻對應(yīng)的開普勒根數(shù)σt。通過這些非線性變換得到樣點的輸出:

Ψj=f(χj)j=0,1,…,12

(9)

最后,計算t時刻位置速度的均值和協(xié)方差矩陣:

(10)

基于UT變換來處理協(xié)方差矩陣的非線性傳遞,可以避開計算笛卡爾坐標對開普勒根數(shù)的偏導(dǎo)數(shù)矩陣以及開普勒根數(shù)的轉(zhuǎn)移矩陣,并且對任意t時刻的協(xié)方差矩陣可以直接計算,不需要計算t時刻之前的協(xié)方差矩陣,能夠有效提高彈道預(yù)報效率。

圖1　系統(tǒng)流程框圖

3仿真分析

仿真將通過普通解析法、基于UT的解析法和數(shù)值法對彈道導(dǎo)彈自由段進行預(yù)報,忽略空氣阻力等其他攝動力,只考慮J2攝動的影響。普通解析法是直接解算攝動運動方程,通過計算開普勒根數(shù)變化率,解算出預(yù)報終點的位置速度;基于UT的解析法即上文所提的方法;數(shù)值法是通過對目標的運動方程進行積分計算彈道的方法,數(shù)值積分采用四階龍格-庫塔算法,仿真步長設(shè)為1s。首先分析解析法中短周期項對彈道預(yù)報的影響,然后分析基于UT的解析法計算的協(xié)方差矩陣和普通解析法進行蒙特卡洛仿真得到的均方根誤差的一致性,最后對比分析解析法和數(shù)值法的預(yù)報精度和運算速度。設(shè)在地心慣性坐標系下,探測到導(dǎo)彈在自由段500s時刻的位置為x=2 163 451m,y=-6 342 006m,z=2 191 877m,速度vx=594.31m/s,vy=711.49m/s,vz=5 453.41m/s,協(xié)方差為

預(yù)報終點t=1 287s。

3.1短周期項對彈道預(yù)報的影響

有短周期項的解析法與數(shù)值法的均值差變化情況及無短周期項的解析法與數(shù)值法的均值變化情況如圖2所示。

圖2　均值之差的變化情況

從圖2中可以看到,無短周期項的解析法與數(shù)值法計算的均值差在位置的x軸、y軸方向及速度的z軸方向均成發(fā)散趨勢,速度的x軸、y軸方向及位置的z軸方向有成周期變化的趨勢。位置的x軸方向終點處偏差約12km,y軸方向終點處偏差達到25km,z軸方向偏差最大達到15km,終點處約14km;速度的x軸方向偏差最大有3m/s,終點處為0.12m/s,y軸方向偏差最大達到14m/s,終點處為12.38m/s,z軸方向終點處偏差為38.08m/s。而與無短周期項相比,考慮了短周期項的解析法與數(shù)值法之差則小得多,接近于0。t=1 287s時刻,均值差的具體結(jié)果如表1所示。由此可見,由于彈道導(dǎo)彈飛行時間短,短周期項對彈道預(yù)報的影響很大,位置偏差能夠達到3個數(shù)量級,所以在進行高精度彈道預(yù)報時不能忽略短周期項。

表1　均值之差仿真結(jié)果

3.2誤差一致性分析

采用UT變換解析法計算的標準差和普通解析法1 000次蒙特卡洛仿真得到的均方根誤差的變化情況如圖3所示。

圖3　普通解析法均方根誤差與基于UT變換解析法的標準差變化情況

從圖3中可以看到,二者的趨勢基本相同,且預(yù)報的協(xié)方差位置和速度均略小于均方根誤差。在預(yù)報的終點處,位置x軸的相對誤差為1%,位置y軸的相對誤差為5.7%,位置z軸的相對誤差為3.03%,速度x軸的相對誤差為1.7%,速度y軸的相對誤差為7.1%,速度z軸的相對誤差為2.6%?？梢?二者的相對誤差非常小,而UT變換計算的協(xié)方差矩陣比蒙特卡洛仿真計算的均方根誤差略小,是因為協(xié)方差矩陣是截斷誤差,忽略了高階項,屬于正常誤差,可以說二者是一致的。

3.3彈道預(yù)報和誤差傳播分析

圖4給出了考慮J2攝動解析法的標準差和數(shù)值法1 000次蒙特卡洛仿真得到的均方根誤差的變化情況。從圖中可以看到,二者的變化趨勢基本一致,解析法的協(xié)方差矩陣在x軸、y軸方向的位置和速度均略大于數(shù)值法的均方根誤差,在z軸方向的位置和速度均略小于數(shù)值法的均方根誤差,說明解析法的精度在x軸和y軸略低于數(shù)值法,而z軸略高于數(shù)值法。在預(yù)報的終點處,位置x軸的相對誤差為4.6%,位置y軸的相對誤差為4.2%,位置z軸的相對誤差為2.5%,速度x軸的相對誤差為4.8%,速度y軸的相對誤差為3.5%,速度z軸的相對誤差為1.1%。無論位置還是速度,兩種方法的相對誤差都不超過5%,說明解析法和數(shù)值法的預(yù)報精度相當,基本在同一數(shù)量級。

考慮J2攝動的解析法的總耗時為128.075ms,其中每個χ點的總耗時為9.851ms。每個χ點在傳播的過程中,500s初始時刻的位置速度轉(zhuǎn)換成開普勒根數(shù)的耗時最大,為4.478ms,占其傳播總耗時的45%左右;計算500s時刻平均根數(shù)的耗時為2.001ms,占20%左右;計算長期項耗時為0.883ms,占9%左右;計算1 287s時刻的平均根數(shù)的耗時為0.616ms,占6%左右;計算短周期項和長周期項的耗時分別為0.755ms和0.698ms,分別占8%和7%左右。可見,在解析法中位置速度與開普勒根數(shù)之間的轉(zhuǎn)換耗時最大。數(shù)值法積分步長為1s,每步長的平均耗時約為1.048ms,預(yù)報自由段時長為788s,數(shù)值法總耗時為826.602ms,是解析法的7倍左右。

預(yù)報飛行時長800s左右的彈道得到相同級別的精度,數(shù)值法計算耗時可以達到考慮J2攝動解析法的7倍左右。原因是數(shù)值法在計算任意時刻導(dǎo)彈狀態(tài)的時候,需要按照一定的步長逐步積分,耗時與所選取的步長成反比,保證一定的精度就需要選取較小的步長,而以耗時的增加為代價;而解析法中,在觀測到導(dǎo)彈初始狀態(tài)后,就能夠快速計算出彈道6個開普勒根數(shù)的變化情況,只要給出任意時刻點,就可以直接計算導(dǎo)彈在該時刻的狀態(tài),所以計算速度快。

圖4　數(shù)值法均方根誤差與考慮J2攝動解析法標準差的變化情況

綜上所述,總結(jié)得到如下結(jié)論:

①由于彈道導(dǎo)彈飛行時間短,短周期項對彈道預(yù)報影響很大,位置偏差能夠達到3個數(shù)量級,在進行高精度彈道預(yù)報的時候不能忽略;

②結(jié)合UT的解析法進行誤差傳播分析,能夠有效避開計算笛卡爾坐標對開普勒根數(shù)的偏導(dǎo)數(shù)矩陣以及開普勒根數(shù)的轉(zhuǎn)移矩陣;

③考慮J2攝動的解析法與積分步長為1s的數(shù)值法相比,彈道預(yù)報的精度在同一數(shù)量級,并且在預(yù)報800s左右的彈道時,耗時只有數(shù)值法的1/7左右。

4結(jié)束語

本文通過仿真對解析法和數(shù)值法在彈道預(yù)報中的性能進行了對比分析。解析法在進行誤差傳播分析時結(jié)合UT,能夠有效避開計算笛卡爾坐標對開普勒根數(shù)的偏導(dǎo)數(shù)矩陣以及開普勒根數(shù)的轉(zhuǎn)移矩陣。傳統(tǒng)的解析法只是簡單地近似成二體問題,沒有考慮其他攝動因素,所以計算精度不高;而考慮了J2攝動之后的解析法,計算精度和誤差傳播發(fā)散程度都與積分步長為1s的數(shù)值法相當,在預(yù)報800s左右的彈道時運算速度是數(shù)值法的7倍左右。但是考慮J2攝動的解析法在彈道預(yù)報中,不能夠像計算衛(wèi)星軌道一樣簡單地忽略短周期項,短周期項對彈道預(yù)報影響很大,在進行高精度的預(yù)報彈道時不能忽略。

參考文獻

[1]張龍杰,謝曉方.通用變步長彈道仿真方法研究[J].彈道學(xué)報,2011,23(2):42-46.

ZHANG Long-jie,XIE Xiao-fang.Research on general variable step-size ballistic simulation method[J].Journal of Ballistics,2011,23(2):42-46.(in Chinese)

[2]程誠,張小兵.某制導(dǎo)炮彈二位兩相流內(nèi)彈道性能分析與數(shù)值模擬研究[J].兵工學(xué)報,2015,36(1):58-63.

CHENG Cheng,ZHANG Xiao-bing.Two-dimensional numerical simulation on two-phase flow interior ballistic performance of a guided projectile[J].Acta Armamentarii,2015,36(1):58-63.(in Chinese)

[3]張毅,楊慧耀,李俊莉.彈道導(dǎo)彈彈道學(xué)[M].長沙:國防科技大學(xué)出版社,1999.

ZHANG Yi,YANG Hui-yao,LI Jun-li.Ballistic missile ballistics[M].Changsha:National University of Defense Technology Press,1999.(in Chinese)

[4]李森,劉紅軍,程仲.提高彈道導(dǎo)彈落點預(yù)報精度方法[J].電子科技,2012,25(9):85-87.

LI Sen,LIU Hong-jun,CHENG Zhong.Improvement of TBM’s impact-point prediction precision[J].Electronic Sci & Tech,2012,25(9):85-87.(in Chinese)

[5]沈慧娜,徐振來.提高彈道導(dǎo)彈落點預(yù)報精度的動弧平滑平均法[J].現(xiàn)代雷達,2009,31(7):55-58.

SHEN Hui-na,XU Zhen-lai.Moving arc smoothing averaging algorithm for improving the prediction precision of ballistic missile impact point[J].Modern Radar,2009,31(7):55-58.(in Chinese)

[6]劉仁,王愛華,郭桂治.基于關(guān)機點狀態(tài)的戰(zhàn)術(shù)彈道導(dǎo)彈落點估計[J].空軍工程大學(xué)學(xué)報,2010,11(1):27-30.

LIU Ren,WANG Ai-hua,GUO Gui-zhi.Impact point estimation of tactical ballistic missile based on the state of burnout point[J].Journal of Air Force Engineering University,2010,11(1):27-30.(in Chinese)

[7]劉彥君,喬士東,黃金才,等.一種高精度彈道導(dǎo)彈落點預(yù)測方法[J].彈道學(xué)報,2012,24(1):23-26.

LIU Yan-jun,QIAO Shi-dong,HUANG Jin-cai,et al.A method of impact point prediction of ballistic missile[J].Journal of Ballistics,2012,24(1):23-26.(in Chinese)

[8]劉林.一種人造地球衛(wèi)星的攝動計算方法[J].天文學(xué)報,1975,16(1):5-80.

LIU Lin.A calculating method of perturbation of artificial satellite[J].Acta Astronomica Sinica,1975,16(1):5-80.(in Chinese)

[9]曹文華,李傳榮,李子揚,等.基于攝動模型的衛(wèi)星軌道預(yù)測算法精度評估[J].遙感信息,2012,27(6):21-27.

CAO Wen-hua,LI Chuan-rong,LI Zi-yang,et al.Accuracy assessment of satellite orbit prediction based onJ2perturbation model[J].Remote Sensing Information,2012,27(6):21-27.(in Chinese)

[10]張俊.基于開普勒二體運動修正地球扁率J2攝動項算法[J].航天控制,2014,32(6):22-25.

ZHANG Jun.Correction algorithm ofJ2perturbations of the earth oblateness based on Kepler two-body motion[J].Aerospace Control,2014,32(6):22-25.(in Chinese)

[11]鄭偉,湯國建.上升段引力常數(shù)變化及其對彈道導(dǎo)彈運動的影響[J].導(dǎo)彈與航天運載技術(shù),2008,295(3):12-14.

ZHENG Wei,TANG Guo-jian.Gravitation constant variety in ascend phase and its effect on ballistic missile[J].Missile and Space Vehcile,2008,295(3):12-14.(in Chinese)

[12]鄭偉,湯國建.彈道導(dǎo)彈自由段解算的等高約束解析解[J].宇航學(xué)報,2007,28(2):269-272.

ZHENG Wei,TANG Guo-jian.Contour restricted analytical solution for free flight trajectory of ballistic missile[J].Journal of Astronautics,2007,28(2):269-272.(in Chinese)

[13]馮偉,劉根友,郝曉光.重力場對彈道導(dǎo)彈自由段落點影響的仿真分析[J].測繪科學(xué),2009,34(5):118-120.

FENG Wei,LIU Gen-you,HAO Xiao-guang.Simulated analysis of the effect of gravity field to the fall point in free flight trajectory of ballistic missile[J],Science of Surveying and Mapping,2009,34(5):118-120.(in Chinese)

[14]賈沛然,陳克俊,何力.遠程火箭彈道學(xué)[M].長沙:國防科技大學(xué)出版社,1993.

JIA Pei-ran,CHEN Ke-jun,HE Li.Long-range rocket ballistics[M].Changsha:National University of Defense Technology Press,1993.(in Chinese)

[15]李濟生,葉杰,王家松,等.航天器軌道確定[M].北京:國防工業(yè)出版社,2003.

LI Ji-sheng,YE Jie,WANG Jia-song,et al.Determination of spacecraft[M].Beijing:National Defense Industry Press,2003.(in Chinese)

[16]MAY J A.A study of the effects of state transition matrix approximations,NASA CP 21123[R].1980.

Analysis on High Precision Trajectory Prediction and Error Propagation of Ballistic Missile ConsideringJ2Perturbation

SUN Yu1,2,WU Nan2,MENG Fan-kun2,DENG Qi-zheng1

(1.Unit 96669 of PLA,Beijing 102208,China;2.College of Navigation and Aerospace Target,Information Engineering University,Zhengzhou 450001,China)

Abstract:To carry out high-precision trajectory prediction and error propagation analysis of ballistic missile,the solution form of mean element algorithms for perturbation equation was analyzed,and an analytical algorithm consideringJ2perturbation and containing the short period term was proposed for trajectory prediction.By constructing the error propagation analysis method of trajectory prediction based on unscented transform,the problem that the Kepler elements is not differentiable and the partial derivative matrix is difficult to calculate was solved.The simulation results show that the error caused by theJ2perturbation short period term,which can not be ignored,reaches 3 orders of magnitude when analytical algorithm is applied in the trajectory prediction.The relative error predicted by the proposed algorithm is only less than 5%,and the precision is far higher than the precision of traditional analytical algorithm,which is comparable with the precision of high-precision numerical algorithm,but the operation speed is 7 times of the latter,and the algorithm has a better effectiveness-cost ratio.

Key words:trajectory prediction;unscented transformation;J2perturbation;analytical algorithm;short-period term

收稿日期:2015-12-30

作者簡介:孫瑜(1986- ),男,助理工程師,碩士,研究方向為空天目標預(yù)警信息處理。E-mail:woshizhongsy@163.com。

中圖分類號:TJ760.1

文獻標識碼:A

文章編號:1004-499X(2016)02-0018-07