陳爾康, 高長生, 荊武興

(哈爾濱工業(yè)大學(xué) 航天學(xué)院, 黑龍江 哈爾濱 150001)

彈道導(dǎo)彈在空間中[1]會(huì)釋放出包括彈頭、誘餌、碎片在內(nèi)的多種空間目標(biāo)，這些目標(biāo)形成了目標(biāo)群。在攻防對抗系統(tǒng)仿真中不論是對目標(biāo)群進(jìn)行跟蹤、識別還是攔截，都需要其較為精確的軌跡[2]和姿態(tài)數(shù)據(jù)，因此對空間目標(biāo)的軌跡和姿態(tài)數(shù)據(jù)生成方法進(jìn)行研究是非常必要的。

在大多數(shù)相關(guān)研究中，目標(biāo)的軌跡數(shù)據(jù)一般根據(jù)橢圓彈道原理計(jì)算得到[3-4]。但這類方法忽略了地球扁率和自轉(zhuǎn)的影響，生成的軌跡數(shù)據(jù)不夠精確。與真實(shí)場景類似，目標(biāo)的彈道由發(fā)射諸元確定，因此軌跡數(shù)據(jù)生成問題實(shí)際上是發(fā)射諸元的計(jì)算問題[5-6]。發(fā)射諸元的計(jì)算方法主要有彈道積分、對射表擬合、解析算法[7-8]和軌跡優(yōu)化[9]等，但這些方法在計(jì)算諸元之前需要進(jìn)行大量的計(jì)算作為準(zhǔn)備。除此之外，目標(biāo)群還還包括分導(dǎo)彈頭和誘餌。這類目標(biāo)在中段由母艙釋放并作被動(dòng)飛行[10]，因此只需計(jì)算出釋放后的初始狀態(tài)即可得到軌跡數(shù)據(jù)。誘餌的初始狀態(tài)可根據(jù)其具體類型設(shè)定釋放相對速度矢量后計(jì)算得到。而分導(dǎo)彈頭的初始狀態(tài)則需滿足命中指定目標(biāo)點(diǎn)的要求，需按照一定的分導(dǎo)策略進(jìn)行計(jì)算[11]。文獻(xiàn)[12]對若干種不同分導(dǎo)策略進(jìn)行了研究并得出了最優(yōu)分導(dǎo)策略。文獻(xiàn)[13]則應(yīng)用極小值原理推導(dǎo)了分導(dǎo)速度增量的近似解析解，但無法得到能夠精確命中設(shè)定目標(biāo)的軌跡數(shù)據(jù)。隨著誘餌技術(shù)的不斷發(fā)展，單純的軌跡信息難以識別出真彈頭，必須考慮目標(biāo)運(yùn)動(dòng)細(xì)節(jié)信息[14]。文獻(xiàn)[15-16]將目標(biāo)除質(zhì)心平動(dòng)以外的微小運(yùn)動(dòng)定義為微動(dòng)并分析了其特點(diǎn)，建立了振動(dòng)、旋轉(zhuǎn)、翻滾和錐旋四類典型微動(dòng)的模型，并發(fā)展成為中段目標(biāo)識別的重要手段[17]。但相關(guān)研究只關(guān)注雷達(dá)信號的結(jié)果，缺少對目標(biāo)姿態(tài)數(shù)據(jù)生成方法的研究。本文從目標(biāo)數(shù)據(jù)生成的角度出發(fā)，根據(jù)彈道導(dǎo)彈的動(dòng)力學(xué)特性和飛行程序設(shè)計(jì)了4種類型的基準(zhǔn)彈道并給出了對應(yīng)的諸元計(jì)算方法。并在此基礎(chǔ)之上，提出了一種分導(dǎo)彈頭最優(yōu)速度增量的實(shí)時(shí)計(jì)算方法。最后針對目標(biāo)的微動(dòng)特性，設(shè)計(jì)了姿態(tài)數(shù)據(jù)的快速生成算法。

1 數(shù)學(xué)模型

1.1 質(zhì)心運(yùn)動(dòng)模型

考慮地球扁率和自轉(zhuǎn)，在發(fā)射系內(nèi)建立彈道導(dǎo)彈質(zhì)心的運(yùn)動(dòng)模型[18]如下

(1)

式中,ρ為導(dǎo)彈相對發(fā)射系的位置;v為導(dǎo)彈相對發(fā)射系的速度;P為推力;R為氣動(dòng)力;g為引力加速度;ae為牽連加速度;ak為哥氏加速度。

氣動(dòng)力在速度坐標(biāo)系中分解為如下形式

(2)

引力加速度沿地心矢徑r方向和地球自轉(zhuǎn)角速度ωe矢量方向分解如下

g=grr+gωeωe

(3)

式中，gr和gωe的表達(dá)式分別為

(4)

式中，μ為地球的引力常數(shù);r為地心距;ae為地球赤道平均半徑;J2為帶諧系數(shù);φ為地心緯度。

1.2 主動(dòng)段飛行程序

彈道導(dǎo)彈的軌跡主要由其主動(dòng)段決定,由于彈道導(dǎo)彈的軌跡近似位于射面以內(nèi),因此其飛行程序主要由其俯仰角變化規(guī)律決定。工程上將其主動(dòng)段飛行程序分為3段[18]:

第1段為垂直上升段,導(dǎo)彈俯仰角保持90°不變。第2段為轉(zhuǎn)彎段,該段內(nèi)攻角變化規(guī)律如下:

(5)

式中，αm為攻角絕對值的最大值,稱為限制攻角。

第3段為瞄準(zhǔn)段，此段中攻角保持固定值αx不變,稱為瞄準(zhǔn)角。

2 基準(zhǔn)彈道設(shè)計(jì)

2.1 彈道類型分析

在沒有瞄準(zhǔn)段能量管理,即瞄準(zhǔn)角為零的情況下,導(dǎo)彈的射程隨限制攻角的增大呈現(xiàn)先增大后減小的變化趨勢。因此某個(gè)中間射程會(huì)對應(yīng)2個(gè)限制攻角值,但這兩個(gè)值所對應(yīng)的彈道高度是不一樣的,一高一低,分別對應(yīng)高彈道和低彈道。

而在瞄準(zhǔn)角不為零時(shí),彈道高度下壓,起到了調(diào)整彈道高度的作用,稱為固定高度彈道。此外,同一個(gè)射程對應(yīng)著多個(gè)限制攻角和瞄準(zhǔn)角的組合,其中存在導(dǎo)彈關(guān)機(jī)點(diǎn)機(jī)械能最小的組合,該組合對應(yīng)的彈道即為最小能量彈道。

2.2 發(fā)射諸元計(jì)算方法

考慮到彈道導(dǎo)彈主要在射面內(nèi)飛行且其射程隨限制攻角先增大后減小的特點(diǎn),采用大步長遍歷粗搜索與精確迭代相結(jié)合的方法計(jì)算發(fā)射諸元,如圖1所示。首先通過球面三角算出發(fā)射方位角的粗搜索值,然后進(jìn)行限制攻角和發(fā)射方位角的遍歷,最后通過精確迭代計(jì)算出發(fā)射諸元。

圖1 發(fā)射諸元計(jì)算流程

3 分導(dǎo)彈道設(shè)計(jì)

3.1 理論分析

為便于分析和仿真計(jì)算,本文采用脈沖變軌的方式并忽略再入段。這樣分導(dǎo)彈道的設(shè)計(jì)就簡化為經(jīng)典Lambert問題[19]。

由分導(dǎo)點(diǎn)掃過地心角θ至目標(biāo)點(diǎn)的時(shí)間ΔT可由拉格朗日轉(zhuǎn)移時(shí)間方程表示[20]

(6)

式中,μ為引力常數(shù);a為軌道半長軸;α和β為拉格朗日參數(shù),cosα=1-s/a,cosβ=1-(s-c)/a,s=(rk+rb+c)/2,c為分導(dǎo)點(diǎn)與目標(biāo)點(diǎn)之間的距離。

因此飛行時(shí)間ΔT與半長軸a之間的關(guān)系如圖2所示。

圖2 ΔT與a的關(guān)系

在ΔT<ΔTm時(shí),ΔT隨a增大而單調(diào)遞減;在ΔT>ΔTm時(shí),ΔT隨a增大而單調(diào)遞增。

而軌道偏心率為[19]

(7)

記分導(dǎo)點(diǎn)處基準(zhǔn)彈道的速度為vc1,分導(dǎo)彈道的速度為v1,則機(jī)動(dòng)速度增量為

(8)

式中

因此增速ΔV與半長軸a之間的關(guān)系如圖3所示。

圖3 ΔV與a的關(guān)系

由圖2和圖3可知:在ΔT<ΔTm時(shí),ΔV隨a增大先減小后增大;在ΔT>ΔTm時(shí),ΔV隨a增大而單調(diào)遞增。因此存在飛行時(shí)間ΔTvmin使得增速ΔV最小

ΔV(ΔTvmin)=minΔV(ΔT)

(9)

3.2 增速計(jì)算方法

根據(jù)3.1節(jié)的分析結(jié)果,本文設(shè)計(jì)了一種搜索最小增速的雙層迭代策略。

對某一固定飛行時(shí)間,構(gòu)造如下指標(biāo)函數(shù)

J(X)=(GT(X)TG(X)

(10)

式中,G(X)為殘差,取為

G(ΔVx,ΔVy,ΔVz)=[rbx-rfx,rby-rfy,rbz-rfz]T

(11)

式中，[rbx,rby,rbz]T為b點(diǎn)的地心矢徑,[rfx,rfy,rfz]T為分導(dǎo)彈頭由分導(dǎo)點(diǎn)K飛行ΔT時(shí)間后的地心矢徑。

采用如下迭代格式使(12)式中指標(biāo)最小

Xk+1=Xk-A-1B

(12)

式中

(13)

經(jīng)過迭代后可得增速為

(14)

對于外層迭代,在增速最小時(shí)有

ΔV′(ΔT)=0

(15)

運(yùn)用牛頓迭代法有迭代公式

(16)

如此反復(fù)迭代即可計(jì)算得到最小增速ΔVmin及其對應(yīng)的飛行時(shí)間ΔTvmin。

4 微動(dòng)姿態(tài)數(shù)據(jù)生成

4.1 理論分析

由于中段目標(biāo)一般為自旋穩(wěn)定,因此其微動(dòng)特征主要是其受干擾力矩作用產(chǎn)生的進(jìn)動(dòng)。在釋放過程中,中段目標(biāo)會(huì)受到?jīng)_量Q的作用[21],可將其按橫向和縱向分解為

(17)

式中,η為沖量Q與目標(biāo)對稱軸的夾角。

進(jìn)動(dòng)周期和進(jìn)動(dòng)角為

(18)

式中,It為目標(biāo)橫向轉(zhuǎn)動(dòng)慣量,Is為目標(biāo)縱向轉(zhuǎn)動(dòng)慣量。

由(17)式和(18)式可知,雖然進(jìn)動(dòng)特性涉及參數(shù)很多,但主要特征為進(jìn)動(dòng)周期和進(jìn)動(dòng)角。因此數(shù)據(jù)生成算法應(yīng)根據(jù)這2個(gè)量生成目標(biāo)的姿態(tài)數(shù)據(jù)。

4.2 姿態(tài)數(shù)據(jù)計(jì)算方法

在仿真中目標(biāo)姿態(tài)一般以相對于地球固聯(lián)坐標(biāo)系姿態(tài)角的形式給出,具體如圖4所示。其中o-XEYEZE為地球固聯(lián)坐標(biāo)系,o-x1y1z1為目標(biāo)的彈體坐標(biāo)系,ψe為偏航角,φe為俯仰角。而滾轉(zhuǎn)角γe為oz1軸與ox1軸、oZE軸所確定平面的夾角。由于本文考慮的目標(biāo)均做自旋運(yùn)動(dòng),因此根據(jù)自旋角速度計(jì)算滾轉(zhuǎn)角即可。

圖4 姿態(tài)示意圖

記沿ox1、oy1軸的單位向量為i,j,經(jīng)坐標(biāo)轉(zhuǎn)換后在地固系的坐標(biāo)分別為[ixe,iye,ize]T、[jxe,jye,jze]T。

則俯仰角、偏航角為

(19)

為便于計(jì)算建立角動(dòng)量慣性坐標(biāo)系,以角動(dòng)量方向(即進(jìn)動(dòng)軸方向)為X軸,任意選擇與X軸構(gòu)成右手坐標(biāo)系的Y軸和Z軸。與體系間的姿態(tài)角分別為進(jìn)動(dòng)角θ、章動(dòng)角ψ和自轉(zhuǎn)角φ。體系到角動(dòng)量慣性系的坐標(biāo)變換矩陣為

(20)

角動(dòng)量慣性系到釋放瞬時(shí)彈體系的坐標(biāo)轉(zhuǎn)換矩陣為

(21)

釋放瞬時(shí)彈體系到地固系的轉(zhuǎn)換陣為

(22)

釋放瞬時(shí)地固系到地固系的轉(zhuǎn)換關(guān)系只與自轉(zhuǎn)角λ=ωeΔT有關(guān)

(23)

由以上坐標(biāo)變化關(guān)系可計(jì)算彈體系到地固系的坐標(biāo)轉(zhuǎn)換矩陣

Ce1=CefCfaCapCp1

(24)

至此已求出了i,j在地固系中的坐標(biāo),代入(19)式即可計(jì)算出姿態(tài)角。

5 仿真分析

5.1 基準(zhǔn)彈道設(shè)計(jì)

對基準(zhǔn)彈道設(shè)計(jì)方法進(jìn)行仿真驗(yàn)證,以某導(dǎo)彈為仿真對象,發(fā)射點(diǎn)設(shè)為(155.462°W,19.653°N),目標(biāo)點(diǎn)設(shè)為(116.395°E,39.922°N),初步計(jì)算目標(biāo)射程為8 481.88 km。本文提出方法(方法1)的彈道設(shè)計(jì)結(jié)果分別如表1和圖5所示。為便于比較,文獻(xiàn)[8]中方法(方法2)的彈道設(shè)計(jì)結(jié)果如表2所示。

表1 基準(zhǔn)彈道設(shè)計(jì)結(jié)果(方法1)

表2 基準(zhǔn)彈道設(shè)計(jì)結(jié)果(方法2)

圖5 基準(zhǔn)彈道

由表1和表2可知，2種方法計(jì)算出的諸元非常相近，說明本文提出的方法能夠給出足夠準(zhǔn)確的彈道數(shù)據(jù)。而在計(jì)算時(shí)間上，雖然方法2的迭代時(shí)間小于方法1的總耗時(shí)，但方法2事先需要大量時(shí)間計(jì)算彈道進(jìn)行擬合，在處理不同的多個(gè)發(fā)射點(diǎn)時(shí)計(jì)算耗時(shí)會(huì)進(jìn)一步增大，因此方法1在應(yīng)用靈活性和計(jì)算耗時(shí)上優(yōu)于方法2。由圖5，4種類型的彈道均以較高的精度命中目標(biāo)點(diǎn)，且由于彈道類型不同其彈道高度也不同，其中固定高度彈道和最小能量彈道的彈道高度位于高彈道和低彈道之間。

5.2 分導(dǎo)彈道設(shè)計(jì)

以5.1節(jié)中的高彈道為基準(zhǔn)彈道進(jìn)行分導(dǎo)彈道設(shè)計(jì)的仿真。設(shè)置分導(dǎo)時(shí)刻為發(fā)射后350 s，分導(dǎo)目標(biāo)點(diǎn)為(117.216°E,39.128°N)，設(shè)計(jì)結(jié)果如圖6所示。其中橫程偏差為0.8 m，縱程偏差為8.36 m，計(jì)算耗時(shí)為0.894 s。利用最優(yōu)分導(dǎo)方法計(jì)算分導(dǎo)彈道，其橫程偏差為35.4 m，縱程偏差為247.6 m，計(jì)算耗時(shí)為0.115 s。雖然該方法計(jì)算耗時(shí)更短，但由于該方法未使用精確模型迭代計(jì)算，因此精度不如本文提出的方法，難以用于系統(tǒng)仿真中。

圖6 分導(dǎo)彈道

可以看到分導(dǎo)彈道以較小偏差命中分導(dǎo)目標(biāo)點(diǎn)，且計(jì)算耗時(shí)較少，能夠滿足仿真需求。

5.3 姿態(tài)數(shù)據(jù)生成

同樣以5.1節(jié)中的高彈道為基準(zhǔn)彈道進(jìn)行誘餌姿態(tài)數(shù)據(jù)生成的仿真。設(shè)置進(jìn)動(dòng)角為4°，自旋周期為0.1 s，錐旋周期為0.5 s，釋放時(shí)間為發(fā)射后400 s。分別使用本文方法(方法1)和傳統(tǒng)的姿態(tài)動(dòng)力學(xué)方程積分方法(方法2)計(jì)算姿態(tài)數(shù)據(jù)，仿真結(jié)果如圖7和圖8所示。

圖7 偏航角

圖8 俯仰角

可以看到本文提出的方法與直接積分動(dòng)力學(xué)方程的計(jì)算結(jié)果一致，但本文方法無需進(jìn)行積分計(jì)算，在計(jì)算時(shí)間上有優(yōu)勢，適于系統(tǒng)仿真使用。

6 結(jié) 論

以系統(tǒng)仿真中的目標(biāo)數(shù)據(jù)生成為研究對象，考慮仿真的特點(diǎn)，對各類目標(biāo)的特性進(jìn)行分析并建立了相應(yīng)的模型，給出仿真數(shù)據(jù)的生成方法：

1)根據(jù)導(dǎo)彈飛行程序特點(diǎn)和彈道特性將基準(zhǔn)彈道分為4類，并設(shè)計(jì)了相應(yīng)的彈道諸元計(jì)算方法。該方法將諸元計(jì)算分為粗搜索和精搜索兩部分，無需事先準(zhǔn)備工作即可快速準(zhǔn)確地生成基準(zhǔn)彈道數(shù)據(jù)。仿真結(jié)果表明，該方法能夠生成準(zhǔn)確的軌跡數(shù)據(jù)；雖然總計(jì)算耗時(shí)大于傳統(tǒng)方法，但傳統(tǒng)方法在迭代前的準(zhǔn)備工作需要耗費(fèi)大量時(shí)間，而本文提出的方法并不需要任何準(zhǔn)備工作，因此該方法能夠更好地滿足系統(tǒng)仿真的要求。

2)將分導(dǎo)彈道的設(shè)計(jì)問題轉(zhuǎn)化為Lambert問題，分析了飛行時(shí)間與分導(dǎo)增速的關(guān)系?；诖岁P(guān)系，提出了一種最小增速分導(dǎo)彈道規(guī)劃方法。這一兩層迭代算法能夠計(jì)算處分導(dǎo)彈頭的最小增速。利用這一增速即可計(jì)算出分導(dǎo)彈頭的彈道數(shù)據(jù)。與傳統(tǒng)方法相比，該方法雖然計(jì)算耗時(shí)小幅增加，但精度明顯提高，更適合系統(tǒng)仿真使用。

3)分析了中段目標(biāo)的姿態(tài)微動(dòng)特性，設(shè)計(jì)了一種通過坐標(biāo)轉(zhuǎn)換矩陣計(jì)算姿態(tài)數(shù)據(jù)的方法，能夠快速準(zhǔn)確的生成包含目標(biāo)微動(dòng)特性的姿態(tài)運(yùn)動(dòng)數(shù)據(jù)。仿真結(jié)果表明，該方法計(jì)算得到的姿態(tài)數(shù)據(jù)與動(dòng)力學(xué)模型積分計(jì)算的結(jié)果一致，能夠準(zhǔn)確反映目標(biāo)微動(dòng)特性，且計(jì)算耗時(shí)短，使用方便，符合仿真的要求。