趙 蒙, 端軍紅, 王明宇, 殷雙斌, 魯 娜

(1. 空軍工程大學防空反導學院, 陜西 西安 710051; 2. 93792部隊, 河北 廊坊 065000)

0 引 言

導彈制導律的設(shè)計研究一直是精確制導控制領(lǐng)域內(nèi)關(guān)注的熱點課題之一[1-6],導彈制導過程一般分為初制導、中制導和末制導,中制導主要作用是控制攔截器飛行至目標附近空域,為末制導提供良好的交班環(huán)境,確保末制導的順利實施[7-9]。對于彈道導彈目標的攔截,傳統(tǒng)破片式殺傷難以有效摧毀彈頭目標,必須使攔截彈直接通過動能碰撞方式實現(xiàn)對彈頭目標的有效攔截,這對攔截彈制導律的制導精度提出了更高要求[10-11]。中制導與末制導有著本質(zhì)的不同,中制導終端約束條件不是攔截器與目標的相對距離為零,而是使攔截彈在中段發(fā)動機關(guān)機后處于“零控攔截狀態(tài)”,大部分學者在中制導律研究方面,主要是使用最優(yōu)化理論[12-15]等方法以消除零控脫靶量為目標[16-19]。文獻[20]將目標和攔截器的引力加速度差簡化后,提出了一種基于零控脫靶量的中制導律,但該制導律只適用于近距離目標攔截情況,在遠距離目標攔截中無法保證制導精度。文獻[21]設(shè)計了一種彈道導彈助推段攔截最優(yōu)制導律,由于無法精確獲得彈道導彈目標推力加速度,該制導律并不實用。文獻[22]基于零控脫靶量設(shè)計了制導律,文獻[23]提出了一種攔截主動防御目標的微分對策制導律,但是在計算過程中都涉及到大量的微積分運算,對彈載計算機的數(shù)據(jù)處理能力要求非常高。文獻[24]設(shè)計了一種基于零控脫靶量的大氣層外超遠程攔截制導方法,在制導過程中由于該制導律的推力方向變化較為劇烈,對攔截器姿態(tài)控制系統(tǒng)要求非常高,不利于工程設(shè)計。文獻[25]提出了一種基于線性重力差模型的攔截彈中制導方法,論文中采用延遲點火的策略去修正預(yù)測模型與實際情況的偏差,但是延時點火的具體時間難以精確計算。

本文針對大氣層外的中段反導作戰(zhàn),基于零控脫靶量采用最優(yōu)控制理論,以脫靶量和攔截彈發(fā)動機的推力加速度為性能指標,設(shè)計了一種攔截彈發(fā)動機固定推進時間的中制導策略,仿真結(jié)果表明,在脫靶時間和脫靶量計算足夠精確條件下,本文設(shè)計的反導攔截彈中制導律脫靶量可控制在百米量級內(nèi),為攔截彈以動能碰撞方式攔截目標創(chuàng)造了良好的末制導條件。

1 相對運動模型建立

在地心慣性坐標系下,設(shè)目標位置矢量為RT、速度矢量為VT,攔截彈位置矢量為RI、速度矢量為VI,則目標相對于攔截彈的位置矢量和速度矢量R=RT-RI、V=VT-VI,反導作戰(zhàn)過程中攔截彈和目標的幾何關(guān)系如圖1所示。

圖1 攔截過程幾何示意圖Fig.1 Geometric diagram of interception process

大氣層外攔截彈和目標的運動方程為

(1)

(2)

式中:rT，rI分別為RT和RI的大小;μ為地球引力常數(shù);aI為攔截彈發(fā)動機推力加速度矢量;f為軌道攝動加速度矢量[26-27]。在地心慣性坐標系下f在3個軸的加速度分量分別為

(3)

式中:x、y、z為空間目標的位置分量;J2=1.082 64×10-3;re為地球赤道半徑。

由式(1)和式(2)可得目標相對于攔截彈運動方程為

(4)

由于反導作戰(zhàn)過程中攔截彈與彈道目標飛行時間相對較短,相對重力加速度而言攝動力產(chǎn)生的加速度非常小,因此目標與攔截彈的攝動加速度之差則可以忽略不計,也就是fT-fI≈0。另外,假設(shè)在目標與攔截彈的相對運動過程中,將攔截彈和目標所受到的地球引力加速度之差簡化為線性模型[20],即

(5)

式中:tf為攔截時刻;G(t0)=GT(t0)-GI(t0)為目標與攔截彈的重力差。

設(shè)tbo(tf≥tbo≥t0)為攔截彈發(fā)動機關(guān)機時刻,則攔截彈在固體發(fā)動機燃料耗盡之后,攔截彈推力加速度aI=0,聯(lián)立式(4)和式(5)可解得攔截彈和目標相對運動方程:

(6)

(7)

攔截彈中制導發(fā)動機關(guān)機后,后續(xù)時間直至末制導攔截彈將無控自由地向目標飛行,因此中制導終端約束不是使攔截彈與目標之間的相對距離為零,而是在中制導結(jié)束后使攔截彈在無控狀態(tài)下飛行至目標,即零控脫靶量為零。

根據(jù)脫靶量定義以及彈目相對運動方程,攔截彈在無控飛行狀態(tài)下,脫靶時刻可根據(jù)下式求解:

R(tf)·V(tf)=0

(8)

將式(6)、式(7)代入式(8)得

(9)

其中,脫靶時刻tf可通過牛頓迭代法求解式(9)獲得,最終可得零控脫靶量預(yù)測值為

(10)

2 攔截彈制導律設(shè)計

將目標相對于攔截彈運動方程改寫成狀態(tài)空間表示形式:

(11)

根據(jù)攔截彈和目標的相對運動方程,令t0=tbo可知

(12)

所以

X(tf)=P1X(tbo)+P2

(13)

式中：

在反彈道導彈作戰(zhàn)中,要實現(xiàn)對彈道導彈的有效攔截,則攔截彈的脫靶量必須足夠小,以保證攔截器與目標的直接碰撞,另外由于攔截彈的有效載荷有限,所以其攜帶的固體燃料也是非常有限的,在中制導過程中還要控制能量消耗。因此,以攔截彈與目標相對距離和攔截彈發(fā)動機的推力加速度作為優(yōu)化指標,設(shè)性能指標函數(shù)[28-29]為

(14)

根據(jù)狀態(tài)空間方程和性能指標函數(shù),可得哈密頓函數(shù)[30-31]為

(15)

正則方程為

(16)

式(16)為一階齊次微分方程組,可解得

λ(t)=e-ATtλ(0)=

(17)

終端指標函數(shù)可以表示成攔截彈發(fā)動機關(guān)機時刻預(yù)測出脫靶時刻攔截彈與目標相對位置的平方形式,因此終端條件可表示為

(18)

所以橫截條件

(19)

所以

(20)

令

根據(jù)最優(yōu)控制理論,極小值條件可通過求解控制方程獲得,由控制方程?H/?U(t)=0可得

U(t)=-BTλ(t)=-BTω(t,tbo)λ(tbo)

(21)

聯(lián)立方程式(19)和式(21)可知

(22)

將式(20)、式(21)代入到狀態(tài)方程式(11)中得

(23)

解方程式(23)可得

(24)

令t=tbo可得

(25)

將式(25)代入式(22),化簡后得到制導加速度:

(26)

式中:tbo≥t≥0;tg=tf-t為剩余飛行時間;tb=tbo-t為發(fā)動機剩余工作時間。

在上述計算過程中,將攔截彈和目標的引力加速度差利用線性化重力模型替代,若目標與攔截彈的相對距離較近,可將目標和攔截彈引力加速度差視為零,即GI(t0)≈GT(t0),則攔截彈發(fā)動機關(guān)機后目標相對于攔截彈的運動方程為

V(t)=V(t0)

(27)

R(t)=V(t0)(t-t0)+R(t0)

(28)

則預(yù)測的脫靶量和脫靶時刻為

(29)

同理,令式(11)中U=-aI(t)、式(13)中X(tf)=P1X(tbo),再次使用最優(yōu)化理論可解得

(30)

將式(10)、式(29)代入式(26)、式(30),可得到通用形式的制導律:

(31)

式中：Rmiss(t)表示為t(tbo≥t≥0)時刻攔截彈與目標之間的零控脫靶量。當攔截彈的推力發(fā)動機工作時間tbo設(shè)定好之后,目標速度和位置狀態(tài)可根據(jù)Kepler理論進行實時估計,攔截彈的位置和速度可由彈載測量設(shè)備實時獲取,攔截彈發(fā)動機點火后可根據(jù)式(31)近似估算的加速度指令制導控制攔截彈飛行。

3 脫靶量的計算

根據(jù)第2節(jié)內(nèi)容可知,制導加速度最終可以表達為脫靶時刻tf、攔截彈發(fā)動機工作時間tbo以及脫靶量Rmiss的函數(shù),在發(fā)動機工作時間以及脫靶時刻確定以后,接下來就是如何精確預(yù)測脫靶量Rmiss。制導律設(shè)計中,終端指標函數(shù)表示為攔截彈脫靶時刻攔截彈與目標相對位置的平方,由于論文中通過近似等效方式獲得攔截彈與目標之間的相對運動方程,與真實的相對運動方程必然存在著差異。

當設(shè)定攔截時間tf、發(fā)動機工作時長tbo后,關(guān)于脫靶量Rmiss的預(yù)測一般有兩種方法。一種是利用等效近似的解析方法來求解脫靶量的大小,這種方法優(yōu)勢是計算速度相對較快,但是計算精度低,例如式(10)、式(29)就是利用近似等效方法求解。另一種是利用Kepler橢圓軌道理論分別迭代計算攔截彈和目標的位置矢量,然后做差來求解脫靶量,這種方法計算精度很高,但是在制導過程中需要持續(xù)迭代求解目標和攔截彈的位置。

近似等效的解析方法論文中已經(jīng)有相關(guān)計算方法,下面重點就如何利用Kepler理論迭代計算脫靶量大小的方法進行介紹,計算流程圖如圖2所示,具體計算步驟如下。

步驟 1初始化,根據(jù)t0時刻攔截彈和雷達獲取的目標運動狀態(tài),在攔截彈推力發(fā)動機工作之后開始計算。

步驟 2根據(jù)初值,計算起始時刻的零控脫靶量。

步驟 3根據(jù)計算出的零控脫靶量計算攔截彈的推力加速度。

步驟 4以Δt為步長,迭代更新計算攔截彈和目標位置矢量和速度矢量,其中Δt越小,精度越高,但是迭代次數(shù)也相應(yīng)增加。

步驟 5判斷該時刻t與tbo的大小,若t≤tbo,則返回步驟 2,令t=t0+Δt繼續(xù)迭代計算;否則結(jié)束計算。

圖2 脫靶量計算流程圖Fig.2 Miss distance calculation flow chart

4 仿真驗證

4.1 仿真條件

制導律仿真驗證時,一般都是以攔截時刻的脫靶量為指標進行衡量,下面就本文提出制導律的有效性進行驗證,參照文獻[10]中的數(shù)據(jù),攔截彈和目標初始運動參數(shù)如表1所示,彈目相距1 560 km,攔截彈中制導段推力發(fā)動機工作時間為15 s。

表1 攔截彈和目標的狀態(tài)信息

4.2 脫靶量仿真

制導加速度最終表達為脫靶時刻tf、攔截彈發(fā)動機工作時間tbo以及脫靶量Rmiss的函數(shù),下面以表1中目標和攔截彈的運動數(shù)據(jù),對脫靶量進行仿真計算。仿真時間從中制導發(fā)動機工作開始至發(fā)動機關(guān)機結(jié)束,共計工作15 s,中制導統(tǒng)一使用式(31)進行制導,將不同的脫靶量計算方法應(yīng)用到攔截彈中制導過程,制導過程中脫靶量的變化結(jié)果如圖3～圖5所示。其中預(yù)測圖3使用式(10)、圖4使用式(29)、圖5使用Kepler迭代法計算脫靶量,真實脫靶量是通過軌道積分計算而得到的。中制導結(jié)束后,不同方法估計的脫靶量和真實脫靶量如表2所示,由仿真結(jié)果可見,Kepler迭代法預(yù)測出的脫靶量精度最高,式(29)預(yù)測的精度最差,式(10)介于兩者之間。

圖3 式(10)預(yù)測出的脫靶量變化曲線Fig.3 Variation curve of miss distance predicted by equation (10)

圖4 式(29)預(yù)測出的脫靶量變化曲線Fig.4 Variation curve of miss distance predicted by equation (29)

圖5 Kepler迭代法預(yù)測出的脫靶量變化曲線Fig.5 Variation curve of miss distance predicted by Kepler iterative method

4.3 仿真數(shù)據(jù)

從圖6可以看出，自慣性坐標系下，加入中制導后攔截彈逐漸向目標運動。下面著重對不同模型下制導律的性能進行仿真比較,仿真計算中,脫靶時刻分別由式(9)、式(29)迭代計算,如圖7所示,其中t1利用式(29)計算,t2利用式(9)計算。

表2 脫靶量計算結(jié)果

圖6 慣性坐標系下攔截彈和目標的運動軌跡Fig.6 Trajectory of interceptor and target in inertial coordinate system

圖7 脫靶時刻的迭代計算結(jié)果Fig.7 Iterative calculation results of miss time

分別利用式(26)、式(30)和式(31)制導控制攔截彈的飛行狀態(tài),使用式(31)制導時,脫靶量的大小利用Kepler理論迭代計算獲得,具體計算流程如圖2所示。中制導過程中攔截彈與目標的脫靶量(彈目相對距離)變化曲線如圖8、圖9所示,其中R1、R2、R3、R4分別表示使用式(30)制導、使用式(26)制導、使用式(31)制導以及攔截彈在無制導時脫靶量變化曲線,攔截效果如表3所示。由仿真結(jié)果可知,使用Kepler迭代法求解出的脫靶量最小,僅為4.48 m,為攔截彈以直接碰撞方式實現(xiàn)對目標的攔截創(chuàng)造了優(yōu)良的末制導條件;使用線性重力差模型得到的脫靶量為330.8 m,還需末制導繼續(xù)修正;使用零重力差模型得到的脫靶量為7 724.8 m,末制導可能無修正能力,可見該模型對于遠距離攔截效果差。仿真結(jié)果表明,中制導的制導精度與脫靶量估計精度密切相關(guān),若脫靶量計算精準,本文所提出的攔截彈中制導律能夠使制導攔截彈以直接碰撞方式實現(xiàn)對目標的攔截。與文獻[10]所提方法相比,本文所設(shè)計的制導律制導精度更高,在迭代計算下制導精度可達米級。

圖8 脫靶量仿真結(jié)果Fig.8 Miss distance simulation results

圖9 脫靶量仿真結(jié)果局部放大圖Fig.9 Partial enlarged view of miss distance simulation results

表3 攔截效果

圖10是使用不同制導公式時攔截彈制導加速度變化曲線,其中a1、a2、a3分別是利用式(30)、式(26)、式(31)計算出的攔截彈發(fā)動機推力加速度大小。由圖可見,a2與a3大小相差不大,但是都與a1相差較大,這是由于式(30)在計算過程中將目標與攔截彈的重力差簡化為零的結(jié)果,這也就是圖9中在脫靶時刻導致脫靶量大的原因,因此為了精確估計零控脫靶量,攔截彈與目標間的重力差必須要合理考慮。

在上述仿真中脫靶時刻是根據(jù)攔截彈和目標運動狀態(tài)預(yù)測出的,一旦攔截彈和目標的運動參數(shù)確定,脫靶時刻也就基本確定了。然而實際反導作戰(zhàn)中,脫靶時刻并不固定,而是根據(jù)攔截彈的發(fā)射時刻、運動參數(shù)等設(shè)定的,脫靶時刻直接影響著攔截彈的運動特征,下面著重討論脫靶時刻對攔截彈加速度指令的影響。分別設(shè)置脫靶時刻為130 s、150 s、170 s、190 s、210 s,利用式(31)進行中制導,仿真結(jié)果如圖11和圖12所示。

圖10 不同制導律下攔截彈加速度變化曲線Fig.10 Acceleration variation curve of interceptor under different guidance laws

圖11 不同脫靶時刻下攔截彈的運動軌跡Fig.11 Trajectory of interceptor at different miss times

圖12 不同脫靶時刻下攔截彈的加速度變化曲線Fig.12 Acceleration variation curve of interceptor at different miss times

其中,圖11為不同攔截時刻攔截彈與目標的運動軌跡,對應(yīng)的脫靶量分別為8.6 m、7.9 m、6.8 m、5.2 m、4.2 m,如表4所示,這進一步驗證了制導算法的有效性。圖12為不同攔截時刻攔截彈的加速度曲線,可見脫靶時刻在130～210 s之間時,攔截彈的機動過載在3g～37g范圍變化,實際在機動過載過大時攔截彈可能無法承受。假設(shè)發(fā)動機能夠提供最大的推力加速度為amax,對于采用固定推進時間的中制導發(fā)動機,為了保證制導過程中發(fā)動機推力加速度不超過amax,需要合理設(shè)計脫靶時刻來實現(xiàn)這一目標。

表4 不同脫靶時刻攔截效果

5 結(jié)束語

本文基于零控脫靶量采用最優(yōu)控制理論,對大氣層外攔截彈和彈道導彈目標的交戰(zhàn)空間進行了理論推導,設(shè)計了一種攔截彈中制導指令,通過仿真計算發(fā)現(xiàn)攔截彈的制導精度與零控脫靶量的預(yù)測精度密切相關(guān),針對固定推進時間的中制導過程討論了調(diào)節(jié)攔截彈機動過載大小的實現(xiàn)途徑。仿真結(jié)果表明,在脫靶量預(yù)測精確條件下,本文設(shè)計的反導攔截彈中制導律脫靶量可控制在百米量級內(nèi),制導精度比較高,在中段制導結(jié)束后,為攔截彈順利實施末制導提供了優(yōu)良的條件。下一步將主要圍繞零控脫靶量的精準估算開展研究,期望相關(guān)研究結(jié)論對反導攔截彈制導律設(shè)計及相應(yīng)的仿真技術(shù)研究提供一定的借鑒和參考。