彈道助推段誤差傳播高精度求解方法

王 磊, 鄭 偉, 周 祥

(國防科技大學(xué)航天科學(xué)與工程學(xué)院, 湖南 長沙 410073)

0 引 言

高命中精度始終是戰(zhàn)略彈道導(dǎo)彈發(fā)展的主要目標(biāo)之一。當(dāng)前,普遍將影響導(dǎo)彈命中精度的因素劃分為制導(dǎo)工具誤差和制導(dǎo)方法誤差兩類[1]。制導(dǎo)工具誤差主要包括加速度計和陀螺儀等慣性測量器件的誤差,也是影響導(dǎo)彈落點精度最主要的誤差源。近年來,隨著新型陀螺儀和加速度計的出現(xiàn),如環(huán)形激光陀螺、光纖陀螺、激光加速度計等,慣性測量器件的精度和穩(wěn)定性逐步提升[2-5](陀螺陀螺漂移量可達(dá)10-6°/h,加速度計零偏穩(wěn)定性可達(dá)1×10-5g);同時,針對不同慣性器件的誤差補償模型也愈加精確[6-7],從而使得制導(dǎo)工具誤差對彈道導(dǎo)彈落點精度的影響逐漸減小。

相對地,制導(dǎo)方法誤差的影響逐漸突顯。該類誤差主要包括制導(dǎo)算法誤差、發(fā)動機后效沖量誤差、再入段誤差、地球引力場模型誤差(通常也稱為地球擾動引力場)等。其中,制導(dǎo)算法誤差和再入段誤差本身的影響量級較小,且可采用一定的補償模型進(jìn)行修正;后效沖量誤差可通過改善發(fā)動機性能或增加末修級來減小或者消除;地球擾動引力場對遠(yuǎn)程彈道導(dǎo)彈落點偏差的影響最大可達(dá)1 km量級[8],是制導(dǎo)方法誤差中最重要的組成部分,其對彈道影響的誤差傳播模型及其相應(yīng)的補償方法也成為當(dāng)前亟待解決的課題。

地球擾動引力場對彈道導(dǎo)彈命中精度的影響主要體現(xiàn)在兩方面,即導(dǎo)彈發(fā)射點垂線偏差引起的慣導(dǎo)系統(tǒng)定向誤差與導(dǎo)彈飛行過程中的擾動引力引起的受力誤差。將EGM2008模型[9-10](earth gravitational model 2008)等高階引力場模型直接嵌入到彈載慣導(dǎo)系統(tǒng)中是實現(xiàn)對擾動引力場影響進(jìn)行補償最理想的方法,但對彈載計算機的執(zhí)行效率及容量有極高的要求。為解決這個問題,研究人員提出了一系列地球引力場重構(gòu)模型來快速逼近真實引力場[11-16],并成功應(yīng)用于潛艇、艦船、飛機等低速載體的慣導(dǎo)系統(tǒng)中[17-23]。但對彈道導(dǎo)彈等高速飛行器而言,采用地球引力場重構(gòu)模型對其慣導(dǎo)系統(tǒng)進(jìn)行直接補償?shù)膶嶋H應(yīng)用仍未見報道。

基于地球擾動引力場對彈道的影響量,反饋求解彈道特征量(如關(guān)機點偏導(dǎo)數(shù))的修正值是實現(xiàn)對地球擾動引力影響進(jìn)行補償?shù)牧硪环N常用思路[8]。彈道攝動法是分析系統(tǒng)性擾動對彈道影響最常用的方法。攝動理論最初是由天文學(xué)家在分析月球、木星及土星等天體在受到太陽光壓及其他引力體攝動后的運動特性的研究中發(fā)展起來的。歐拉、拉格朗日、高斯、拉普拉斯等人也對攝動理論的發(fā)展也作出了突出的貢獻(xiàn)[24-25]。我國學(xué)者自20世紀(jì)80年代以來也對擾動引力場作用下的彈道攝動問題進(jìn)行了系統(tǒng)而深入的研究。文獻(xiàn)[26]在不考慮視速度偏差的假設(shè)下研究了導(dǎo)彈飛行過程中擾動引力對遠(yuǎn)程彈道導(dǎo)彈慣性制導(dǎo)系統(tǒng)及命中精度的影響特性;文獻(xiàn)[27]則對彈道導(dǎo)彈被動段誤差傳播攝動方程進(jìn)行了推導(dǎo),并提出了解析的計算方法;文獻(xiàn)[28]推導(dǎo)了垂線偏差、發(fā)射點大地經(jīng)緯度偏差及發(fā)射方位角偏差等定位定向誤差對彈道助推段關(guān)機點狀態(tài)影響的攝動方程。這些方法在擾動引力場對導(dǎo)彈命中精度影響特性分析的應(yīng)用中均取得了較好的效果,但由于沒有考慮狀態(tài)偏差與視加速度之間的耦合影響,其精度在一些初始條件下會有不足;另一方面,現(xiàn)有方法中考慮垂線偏差和飛行過程中擾動引力影響的攝動方程各自獨立,無法對其影響進(jìn)行統(tǒng)一分析。

基于此,本文推導(dǎo)了一種可同時考慮發(fā)射點垂線偏差和導(dǎo)彈飛行過程擾動引力影響的彈道助推段誤差傳播求解方法,該方法也考慮了狀態(tài)偏差引起的視加速度偏差攝動項,求解精度更高。本文正文結(jié)構(gòu)可分為3部分:第一,基于狀態(tài)空間攝動法推導(dǎo)了彈道導(dǎo)彈在受到擾動引力場作用時的攝動方程,這是推導(dǎo)彈道誤差傳播模型的理論基礎(chǔ);第二,通過對攝動方程進(jìn)行分析,深入探索了攝動源的構(gòu)成及其影響機理;最后,通過數(shù)值仿真分析了擾動引力場對彈道助推段的影響量級,并對比驗證了本文所提誤差傳播方法的計算效率和精度。

1 彈道助推段攝動方程

1.1 坐標(biāo)系的定義

彈道助推段誤差傳播模型主要在發(fā)射慣性坐標(biāo)系中描述。如圖1所示,在考慮發(fā)射點垂線偏差的情況下,發(fā)射慣性坐標(biāo)系可以有如下兩種定義方式:

(1)以參考橢球法線為定向基準(zhǔn)的標(biāo)準(zhǔn)發(fā)射慣性系On-xeyeze；

(2)以當(dāng)?shù)刂卮咕€為定向基準(zhǔn)的發(fā)射慣性系On-xayaza。

圖1 發(fā)射慣性系定義示意圖Fig.1 Definition of the launch inertial coordinate system

顯然,坐標(biāo)系On-xeyeze與On-xayaza之間的差別由垂線偏差決定,這二者之間的轉(zhuǎn)換關(guān)系為

(1)

式中,ρe,ρa,ve,va分別為對應(yīng)坐標(biāo)系中的位置和速度矢量;EG為發(fā)射系到地心系的方向余弦陣,且

(2)

λ0,B0,A0為大地經(jīng)緯度和大地方位角;λT,BT,AT為天文經(jīng)緯度和天文方位角。根據(jù)大地測量理論可知它們與垂線偏差子午分量ξ和卯酉分量η之間的關(guān)系為

(3)

1.2 攝動方程

在橢球地球假設(shè)下,導(dǎo)彈助推段彈道運動方程為

(4)

考慮地球擾動引力場影響時,實際彈道運動方程為

(5)

令DG表示On-xayaza到On-xeyeze的坐標(biāo)轉(zhuǎn)換矩陣,根據(jù)式(1)有

(6)

將式(5)在On-xeyeze中表達(dá),可得

(7)

式中,ge=DGga為在On-xeyeze中表達(dá)的真實引力矢量。

求解式(7)與式(4)的等時變分可得

(8)

(9)

式(8)即為擾動引力場影響下彈道助推段運動攝動微分方程。

2 誤差傳播模型推導(dǎo)

為方便分析,將式(8)改寫為矩陣的形式

(10)

其中

(12)

對攝動方程式(10)進(jìn)行分析可知,擾動引力場對彈道的影響主要包括4部分:

(2)由引力模型誤差引起的受力項,即δge;

(3)由狀態(tài)偏差引起的引力加速度耦合項,即A·[δve,δρe]T;

(4)由狀態(tài)偏差引起的視加速度耦合項,即B·[δve,δρe]T。

2.1 幾何項模型

(13)

數(shù)值仿真表明,推力T是產(chǎn)生視加速度的主要來源,忽略氣動力及控制力等小量的影響,幾何項模型可簡化為

(14)

顯然,幾何項的量級與發(fā)射點的垂線偏差及推力的大小直接相關(guān),垂線偏差越大、推力越大,幾何項就越大,反之則越小。

2.2 受力項模型

地球的真實引力無法精確獲得,通常只能通過模型不斷逼近。求解地球擾動引力最常用的方法有點質(zhì)量法和球諧函數(shù)法。前者只能針對固定區(qū)域進(jìn)行求解,后者則可以實現(xiàn)全球范圍內(nèi)任意點的求解。本文仿真中主要采用后者進(jìn)行計算,且采用EGM2008的模型系數(shù)。

擾動引力位在地球外部是一調(diào)和函數(shù),滿足拉普拉斯方程,基于分離變量法,可將其展開成球諧函數(shù)[1]

(15)

擾動引力位對位置的梯度即為擾動引力,即

(16)

由式(16)可知,擾動引力的計算是在地心系內(nèi)完成,彈道計算時需要將其轉(zhuǎn)換到發(fā)射慣性系中,即

(17)

2.3 引力加速度耦合項模型

引力加速度耦合項即是指由狀態(tài)偏差引起的引力矢量偏差。該項重點在于求解發(fā)射慣性系中引力矢量關(guān)于速度、位置矢量的雅克比矩陣,根據(jù)文獻(xiàn)[28]的推導(dǎo),可知矩陣A的表達(dá)式為

(18)

其中

2.4 視加速度耦合項模型

視加速度耦合項即是指由狀態(tài)偏差引起的視加速度矢量偏差。該項重點在于求解發(fā)射慣性系中視加速度矢量關(guān)于速度、位置矢量的雅克比矩陣。忽略控制力的影響,視加速度偏差可表示為

(19)

式中,R和T的表達(dá)式為

(20)

式中,Cx、Cy、Cz分別為阻力系數(shù)、升力系數(shù)和側(cè)力系數(shù);ρ為大氣密度;v為飛行器相對大氣的速度;Sm為彈體最大橫截面積;m為彈體質(zhì)量;ue為排氣速度;Se為噴口截面積;pe為排氣端面壓力;pH為當(dāng)?shù)卮髿鈮毫Α?/p>

記Mv為R對速度矢量的偏導(dǎo)數(shù),Mr為R對位置矢量的偏導(dǎo)數(shù),具體表達(dá)式為

(21)

其中

M1=-Cxρvx;M2=-Cxρvy;M3=-Cxρvz

且有

(22)

記Nr為T對位置矢量的偏導(dǎo)數(shù),具體為

(23)

可得視加速度相對導(dǎo)彈狀態(tài)矢量的雅克比矩陣為

(24)

即

(25)

3 擾動引力場對彈道影響特性分析

3.1 擾動引力場對助推段彈道的影響特性分析

以民兵-3導(dǎo)彈總體參數(shù)、氣動參數(shù)和我國標(biāo)準(zhǔn)大氣參數(shù)為基礎(chǔ)構(gòu)建了導(dǎo)彈誤差傳播仿真系統(tǒng)?；谠摲抡嫦到y(tǒng)首先分析了攝動模型不同攝動項的影響量級,進(jìn)而采用攝動法與彈道求差法對彈道誤差傳播特性進(jìn)行仿真分析,并對比驗證本文所提攝動法的計算效率及精度。

發(fā)射點參數(shù)及垂線偏差信息如表1所示,導(dǎo)彈飛行過程中的擾動引力采用360階的球諧函數(shù)進(jìn)行計算。

表1 導(dǎo)彈發(fā)射點參數(shù)

圖2為4類攝動項的量級及其隨時間的變化曲線?？梢钥闯?在表1所示的初始條件下,幾何項在所有攝動項中所占比例最大,最高可達(dá)500 mgal (1 mgal=10-5m/s2), 且其變化曲線與發(fā)動機推力曲線基本一致,符合前面的分析結(jié)論;受力項的量級僅次于幾何項,且隨著時間的增加而減小;視加速度耦合項隨時間先增大后減小,當(dāng)導(dǎo)彈飛出大氣層外則減小為零;在整個助推段時間周期內(nèi),引力加速度耦合項的量級均在2 mgal以下,因而本文后續(xù)仿真中均忽略了該項。

圖2 不同攝動項隨時間變化曲線圖Fig.2 Curves of different perturbation accelerations over time

下面對擾動引力場影響下的彈道誤差傳播特性進(jìn)行仿真分析。彈道求差法的仿真流程可描述為:在由λ0,B0,A0確定的發(fā)射慣性系內(nèi)和由λT,BT,AT確定的發(fā)射慣性系內(nèi)采用同一時間步長分別進(jìn)行彈道積分,且前者在積分時引力模型只考慮到J2項,后者則考慮到360階擾動引力項。將兩條彈道狀態(tài)量(速度和位置)在同一坐標(biāo)系內(nèi)等時求差即可獲得彈道真實誤差傳播結(jié)果。該方法的誤差源僅為計算機本身累積誤差,可以作為本文所提攝動模型精度評估的參照。為保證積分精度,彈道積分步長應(yīng)足夠小,通常使其與制導(dǎo)步長一致,一般取0.02 s。

與彈道求差法不同,攝動法的積分變量為導(dǎo)彈狀態(tài)的偏差,其量級小且變化緩慢,因而可以采用較大的時間步長進(jìn)行積分,其步長大小基于如下流程確定:

步驟1設(shè)置其步長初值Δt0=0.02 s,并按照當(dāng)前步長采用攝動法求解助推段關(guān)機點狀態(tài)偏差ΔX0;

步驟2按照0.01 s的增量,逐步增加步長值,同時求解相應(yīng)步長下的關(guān)機點狀態(tài)偏差ΔXk;

步驟3計算ΔXk相對ΔX0的誤差百分比ε,并判斷ε是否大于設(shè)定的相對誤差容許值,若成立則終止循環(huán),最終步長取為循環(huán)中上一步的步長,若不成立,則返回到步驟2。

本文仿真中,將相對誤差容許值設(shè)置為0.001,并依據(jù)上述流程求得攝動法可采用的積分步長為1 s。同時,根據(jù)是否忽略視加速度耦合項將攝動法分為兩類,即攝動法-1和攝動法-2。

圖3為彈道求差法、攝動法-1和攝動法-2這3種方法求解彈道助推段誤差傳播的結(jié)果對比曲線,其中圖3(a)～圖3(c)依次為位置矢量偏差隨時間的變化曲線,圖3(d)～圖3(f)為速度矢量偏差隨時間的變化曲線?？梢钥闯?本文提出的誤差傳播模型能較好反應(yīng)地球擾動引力場影響下的彈道助推段誤差傳播特性;同時可知,若忽略攝動方程中的視加速度耦合項,會使誤差傳播模型的精度有明顯的下降。

圖3 沿助推段彈道的導(dǎo)彈狀態(tài)偏差計算結(jié)果Fig.3 Calculation results of the state error in the boost phase

根據(jù)第2.4節(jié)的分析可知,視加速度耦合項反應(yīng)了導(dǎo)彈狀態(tài)偏差引起的氣動力及發(fā)動機推力偏差,且主要引起氣動力偏差。圖2表明該攝動項主要在10～70 s的飛行時間內(nèi)起作用,這是因為導(dǎo)彈在該飛行時間段內(nèi)主要處于稠密大氣層內(nèi),氣動力對導(dǎo)彈狀態(tài)的變化比較敏感;當(dāng)導(dǎo)彈飛出大氣層后,該攝動項則幾乎不再起作用。圖3表明,為了獲得較高精度的誤差傳播模型,視加速度耦合項不應(yīng)該被忽略。

3.2 擾動引力場對導(dǎo)彈射程的影響特性分析

為充分驗證本文所提彈道助推段誤差傳播模型的效率、精度及其對不同初始條件的適應(yīng)性,首先采用遍歷法仿真分析不同導(dǎo)彈發(fā)射方位角條件下兩種攝動模型相對彈道求差法的精度,其次采用蒙特卡羅打靶法分析不同垂線偏差條件下兩種攝動模型相對彈道求差法的精度及其誤差統(tǒng)計特性。

結(jié)果分析時將助推段關(guān)機點狀態(tài)偏差導(dǎo)致的導(dǎo)彈射程偏差ΔL設(shè)為待評價指標(biāo),且ΔL的計算公式為

(26)

式中,等式右邊行向量為導(dǎo)彈射程關(guān)于導(dǎo)彈關(guān)機點狀態(tài)的偏導(dǎo)數(shù),通常稱為關(guān)機點射程偏導(dǎo)數(shù),列向量為導(dǎo)彈關(guān)機點狀態(tài)偏差。本文仿真中關(guān)機點射程偏導(dǎo)數(shù)取文獻(xiàn)[27]中給定的值,即有

3.2.1 發(fā)射方位角遍歷

仿真中,方位角由-180°～175°按照5°等間隔遍歷,共計算72條彈道,其余發(fā)射點參數(shù)與表1所示一致。仿真結(jié)果如圖4和圖5所示。圖4顯示了彈道求差法和攝動法-2計算耗時對比曲線。顯然,完成一次給定初始條件下彈道助推段誤差傳播仿真計算,求差法平均耗時1 025.3 s,而攝動法-2只需要23.4 s。其耗時之比為50.2,與這兩種方法積分步長之比的倒數(shù)基本一致。實際上,采用球諧函數(shù)計算導(dǎo)彈飛行過程中的擾動引力是整個積分計算中最耗時的部分,幾乎占據(jù)了總耗時的99%,且計算時間隨著球諧階數(shù)級數(shù)的增加呈冪級數(shù)增加。這也成為制約彈道求差法計算效率的最主要原因。

圖4 計算時間對比Fig.4 Comparison of time-consuming between the two methods

圖5顯示了擾動引力場引起的射程偏差及攝動法-2的計算精度隨著方位角的變化情況。

圖5 不同方位角條件下的射程偏差求解精度對比Fig.5 Accuracy comparison of range error calculation under different azimuth conditions

3.2.2 垂線偏差隨機打靶

設(shè)置發(fā)射點處垂線偏差的子午分量及卯酉分量都為滿足正態(tài)分布的隨機數(shù),且其均值和方差均取5″。仿真中,在發(fā)射方位角為-160°、20°和140° 3種條件下分別進(jìn)行500次的蒙特卡羅打靶,并統(tǒng)計攝動法-2計算所得射程偏差的相對誤差百分比。

仿真結(jié)果如圖6所示,圖6(a),圖6(c)和圖6(e)分別為3種發(fā)射方位角條件下的攝動法-2相對誤差百分比的散點圖,圖6(b),圖6(d)和圖6(f)則分別為其對應(yīng)的直方圖統(tǒng)計圖。結(jié)果表明,3種方位角條件下攝動法-2的相對誤差百分比都小于10%,且均值分別為1.834 0%,6.663 1%和2.525 4%。

圖6 垂線偏差隨機條件下的射程偏差求解精度對比Fig.6 Accuracy comparison of range error calculation under different deflection of the vertical conditions

4 結(jié) 論

地球擾動引力場對遠(yuǎn)程彈道導(dǎo)彈落點偏差的影響在百米到1 km量級。導(dǎo)彈發(fā)射前快速分析擾動引力場對擬飛行彈道的影響是實現(xiàn)擾動引力場影響精確補償?shù)挠行侄沃弧１疚闹饕槍@一需求,基于狀態(tài)空間攝動法提出一種可同時考慮發(fā)射點垂線偏差和導(dǎo)彈飛行過程擾動引力影響的彈道助推段誤差傳播高精度求解方法。仿真結(jié)果表明,考慮視加速度耦合項的誤差傳播模型在計算精度方面較不考慮視加速度耦合項的情況有明顯提升;同時,相對彈道求差法而言,考慮視加速度耦合項的誤差傳播模型計算效率提升了約50倍,計算相對誤差均小于10%?？梢詾閷崿F(xiàn)擾動引力場影響快速補償提供理論依據(jù)和方法支撐。