喬清青，陳萬春

(北京航空航天大學(xué)宇航學(xué)院，北京100191)

0 引言

高超聲速導(dǎo)彈由于具有飛行高度高、速度快、機(jī)動(dòng)性好等特點(diǎn)，能有效縮短對(duì)目標(biāo)的反應(yīng)時(shí)間，具有較高的突防成功率［1］。高超聲速導(dǎo)彈飛行距離較遠(yuǎn)，需要一種中制導(dǎo)律將導(dǎo)彈導(dǎo)引至目標(biāo)附近，為末端作戰(zhàn)保存能量，滿足各種約束條件，并完成中末制導(dǎo)的平滑交接班。

導(dǎo)彈的中制導(dǎo)問題可歸結(jié)為一個(gè)非線性兩點(diǎn)邊值問題。由于導(dǎo)彈動(dòng)力學(xué)系統(tǒng)的高度非線性，兩點(diǎn)邊值問題無法求出解析解。奇異攝動(dòng)是解決復(fù)雜非線性最優(yōu)控制問題的一種經(jīng)典方法，通過抑制動(dòng)力學(xué)方程中的小參數(shù)使系統(tǒng)階數(shù)降低，可容易地獲得相應(yīng)的降階解，再計(jì)算小參數(shù)獲得邊界層解，以修正降階解。使用奇異攝動(dòng)可得到性能較理想的近似解，并可根據(jù)需要獲得開環(huán)或閉環(huán)解。奇異攝動(dòng)方法在航空航天領(lǐng)域得到了廣泛的應(yīng)用。由于飛行器動(dòng)力學(xué)方程中小參數(shù)不明顯，因此通常使用“強(qiáng)迫奇異攝動(dòng)”方法，即先根據(jù)估算的各變量的時(shí)間常數(shù)以及應(yīng)用經(jīng)驗(yàn)區(qū)分變化快慢不同的變量，再在快變量前插入小參數(shù)ε(0＜ε＜1)構(gòu)成奇異攝動(dòng)系統(tǒng)。令小參數(shù)為零，即假設(shè)較快變量能在瞬間達(dá)到平衡，從而使系統(tǒng)階數(shù)降低，并求解出降階解，然后通過邊界層修正計(jì)算考慮被忽略的快變量及其邊界條件。用奇異攝動(dòng)方法求解飛行器制導(dǎo)非線性兩點(diǎn)邊值問題，可獲得算法較簡單，并可在彈上實(shí)時(shí)解算的近最優(yōu)制導(dǎo)律［2-3］。

文獻(xiàn)［4］針對(duì)超聲速飛行器的最大航程滑行問題，將動(dòng)力學(xué)方程分成三個(gè)層次，用奇異攝動(dòng)方法求解了一種最大航程閉環(huán)導(dǎo)引律。文獻(xiàn)［5］考慮了球形地球上的飛行器質(zhì)點(diǎn)動(dòng)力學(xué)模型，推導(dǎo)了奇異攝動(dòng)最大航程制導(dǎo)律，給出了該算法與用序列二次規(guī)劃(SQP)算法解算出的最優(yōu)彈道和熱流率曲線。無熱流約束的彈道均為跳躍-滑翔式彈道，飛行中產(chǎn)生的熱流率也較大。

針對(duì)上述問題，本文利用文獻(xiàn)［6］中的高超聲速飛行器的氣動(dòng)模型，用奇異攝動(dòng)和動(dòng)態(tài)逆方法為高超聲速導(dǎo)彈推導(dǎo)了一種使熱流率較低的最大末速中制導(dǎo)律。

1 導(dǎo)彈質(zhì)點(diǎn)動(dòng)力學(xué)模型

式中，r，h，E，γ，m，n 和 R0分別為航程、高度、比能、彈道傾角、質(zhì)量、過載和地球平均半徑。問題的初始條件為:

式中，CD0和CDt分別為零升阻力系數(shù)和升致阻力系數(shù)。動(dòng)壓定義為q=ρV2/2。各氣動(dòng)力系數(shù)見文獻(xiàn)［4］。導(dǎo)彈迎角為:

式中，m，g，S分別為導(dǎo)彈的質(zhì)量、當(dāng)?shù)刂亓铀俣群蛥⒖济娣e，“～”表示變量的實(shí)際值。

令導(dǎo)彈飛行時(shí)產(chǎn)生的最大熱流率為:

式中，熱流率單位為W/cm2。

文獻(xiàn)［6］用奇異攝動(dòng)方法為高超聲速飛行器推導(dǎo)了一種最遠(yuǎn)航程制導(dǎo)律，并用其與配點(diǎn)法和SQP方法解算出的最優(yōu)彈道進(jìn)行比較，證明了高超聲速導(dǎo)彈的無約束彈道為跳躍-滑翔彈道。當(dāng)導(dǎo)彈飛行至較高高度時(shí)，指令過載很小，甚至趨近于0。跳躍式彈道可很好地為導(dǎo)彈保存能量，但產(chǎn)生的熱流率較大，不適于工程實(shí)際應(yīng)用。

根據(jù)文獻(xiàn)［7］中的結(jié)論，要減小熱流率，就應(yīng)使飛行彈道盡量平滑。用一些尋優(yōu)方法，經(jīng)過反復(fù)迭代運(yùn)算可以得到較理想的結(jié)果，但其所需運(yùn)算時(shí)間長，無法在飛行器上實(shí)時(shí)進(jìn)行。另外，目標(biāo)的運(yùn)動(dòng)，以及導(dǎo)彈在實(shí)際飛行中受到的各種干擾，均會(huì)導(dǎo)致飛行狀態(tài)發(fā)生改變，從而無法滿足跟蹤最優(yōu)值所需的條件。

2 奇異攝動(dòng)中制導(dǎo)律

由于高超聲速導(dǎo)彈飛行時(shí)間較長，需要在中制導(dǎo)段盡量保存能量，以保證末端攔截有足夠的能量。問題變成使微分方程式(1)～式(4)所描述的系統(tǒng)從初始狀態(tài)轉(zhuǎn)移到預(yù)測攔截點(diǎn)在笛卡爾坐標(biāo)系中的末端位置(rf，hf)，使得末端比能最大。假設(shè)目標(biāo)以當(dāng)前速度作勻速直線運(yùn)動(dòng)，則經(jīng)過剩余飛行時(shí)間tgo后，其所處位置作為預(yù)測攔截點(diǎn)。剩余飛行時(shí)間估算如下:

式中，Rtm為導(dǎo)彈與目標(biāo)的距離;Vc為接近速度;a為導(dǎo)彈沿視線方向的平均加速度，認(rèn)為與導(dǎo)彈軸向加速度相等。

用文獻(xiàn)［8］中介紹的方法估算各狀態(tài)變量的時(shí)間常數(shù)。根據(jù)其大小可將變量分成:慢變量(r，E)、中間變量h和快變量γ。在中間變量和快變量前分別插入小參數(shù) ε1(0＜ε1＜1)和 ε2(0＜ε2＜1)，可得到導(dǎo)彈的質(zhì)心動(dòng)力學(xué)奇異攝動(dòng)模型為:

為使末端速度最大，可令性能指標(biāo)為:J=－Ef。

2.1 降階系統(tǒng)

最慢變量構(gòu)成了降階系統(tǒng)，由降階系統(tǒng)求解出的解稱為外解。令ε1=ε2=0，得到γ=0，n=1－V2/(R0+h)。帶入式(7)～式(10)，得到降階系統(tǒng)為:，其中:D=D0+K［1－V2/(R0+h)］2/q;高度h為零階系統(tǒng)的控制量。Hamilton函數(shù)寫成:H=λrVR0/(R0+h)－λEDV。由伴隨方程，可知 λr為常數(shù)。由橫截條件 λE(tf)= －1和H(tf)=0，得:

系統(tǒng)中各變量與時(shí)間無關(guān)，因此H(t)=0，可求得:

由控制方程?H/?h=0，帶入式(11)和式(12)可得:

式(13)為h和E的隱函數(shù)。為加快彈上運(yùn)算速度，可先將不同比能下的?D/?h與高度h的關(guān)系列成表格，這樣通過插值即可求出與某指定能量下對(duì)應(yīng)的最優(yōu)高度h*。根據(jù)第1節(jié)中導(dǎo)彈氣動(dòng)力系數(shù)，計(jì)算出的比能與降階最優(yōu)高度關(guān)系如圖1所示。

圖1 降階系統(tǒng)最優(yōu)高度與比能的關(guān)系Fig.1 Relationship between reduced order optimal altitude and specific energy

本節(jié)中的V，D，h均為降階值，下文中用上標(biāo)“o”表示外解中的變量。降階系統(tǒng)中忽略了高度變量的邊界值。為使導(dǎo)彈能夠達(dá)到預(yù)定位置，對(duì)降階最優(yōu)高度作如下修正:

式中，kh為待選系數(shù);rtm為兩者之間的水平距離。

2.2 第一邊界層

在飛行器質(zhì)量不變和一定速度條件下，飛行器受到的零升阻力隨海拔高度的增加而減小，由相同過載產(chǎn)生的阻力隨海拔高度的增加而增加。因此，當(dāng)指令過載不為0時(shí)，存在一個(gè)使阻力達(dá)到極小值的高度。

從物理上看，式(13)的意義在于:在某能量狀態(tài)下，在降階最優(yōu)高度上作水平飛行時(shí)，飛行器受到的阻力最小。為得到一較平滑的彈道，邊界層修正得到的系統(tǒng)控制量應(yīng)使導(dǎo)彈能在降階最優(yōu)高度上飛行。

令 ε1=1，ε2=0，則 n=cos γ［1 － V2/(R0+h)］，得到的第一邊界層動(dòng)態(tài)由高度動(dòng)力學(xué)方程決定。該層的控制量為彈道傾角γ。用動(dòng)態(tài)逆方法對(duì)該層進(jìn)行求解，第一邊界層的系統(tǒng)方程和輸出方程為:

該系統(tǒng)并非仿射非線性系統(tǒng)，因此將控制量改設(shè)為:

式中，ωh為待定系數(shù)。系統(tǒng)的控制為:

將式(14)帶入，整理后可得到期望彈道傾角的表達(dá)式為:

對(duì)導(dǎo)彈系統(tǒng)的時(shí)間尺度分析可知，速度和高度變量在相同的時(shí)間尺度上，因此該層中的速度變量可看作當(dāng)前導(dǎo)彈速度V?？疾焓?15)所示的一階系統(tǒng)，其性能主要取決于h的時(shí)間常數(shù)。

下文中用上標(biāo)“i1”表示第一邊界層中的變量。該層計(jì)算中忽略了控制量彈道傾角的邊界條件。為了使導(dǎo)彈在中制導(dǎo)末端，導(dǎo)引頭能夠順利截獲目標(biāo)，需對(duì)最優(yōu)進(jìn)行修正:

其中:

2.3 第二邊界層

令ε1=ε2=1，可得到第二邊界層動(dòng)力學(xué)方程，為系統(tǒng)方程完整形式。其控制變量為指令過載n。第二邊界層系統(tǒng)動(dòng)態(tài)主要由γ的動(dòng)力學(xué)方程決定，控制量為n。其系統(tǒng)方程和輸出方程為:

第二邊界層系統(tǒng)為原系統(tǒng)中最快的一層，認(rèn)為其中的變量均可取為當(dāng)前值?？煽闯觯鲜剿鞠到y(tǒng)是一個(gè)仿射非線性系統(tǒng)，可直接構(gòu)造逆系統(tǒng)并求解控制律。求得控制量為:

式中，ωγ為待定系數(shù)。

2.4 中制導(dǎo)算法

通過式(13)得到降階最優(yōu)后，用式(16)和式(17)即可計(jì)算得到導(dǎo)彈的指令過載。

根據(jù)動(dòng)態(tài)逆理論，式(16)和式(17)中的待定系數(shù)ωh和ωγ應(yīng)分別為高度變量和彈道傾角變量的時(shí)間常數(shù)的倒數(shù)。用文獻(xiàn)［4］介紹的方法可估算變量的時(shí)間常數(shù)，但該方法依賴人工選取的典型值，因此估算出的時(shí)間常數(shù)在某區(qū)間內(nèi)均可取值。文獻(xiàn)［9］中介紹了如何根據(jù)二階系統(tǒng)階躍響應(yīng)分析方法選取合適的時(shí)間常數(shù)值的方法，然而由于高超聲速導(dǎo)彈飛行速度快，可達(dá)到很高的高度，導(dǎo)致其受到的阻力對(duì)指令過載的變化相當(dāng)敏感。因此，當(dāng)導(dǎo)彈處于很大的高度時(shí)，其指令過載趨近于0，這樣才不會(huì)使飛行器能量迅速損失。

對(duì)于本文中的例子，導(dǎo)彈在初始一段時(shí)間內(nèi)從較大的高度下落，需對(duì)初始下降段的降階最優(yōu)高度和第一邊界層最優(yōu)彈道傾角進(jìn)行修正:當(dāng)h=92 008 m時(shí)，將導(dǎo)彈以當(dāng)前狀態(tài)作自由下落的彈道高度作為降階最優(yōu)高度;將導(dǎo)彈以當(dāng)前狀態(tài)作自由下落的彈道傾角作為第一邊界層最優(yōu)彈道傾角，以保證導(dǎo)彈在最初下落段指令過載趨近于0。

3 仿真算例

本文用Matlab軟件對(duì)制導(dǎo)律進(jìn)行了仿真驗(yàn)證。導(dǎo)彈的初始狀態(tài)為:r0=0 m，h0=112 080 m，V0=5164.9 m/s。設(shè)目標(biāo)初始狀態(tài) ht(t0)=10 km，rt(t0)=5 500 km。目標(biāo)以Ma=2的速度向?qū)椨鎰蛩俚雀唢w行。仿真過程中，令迎角α≤25°，導(dǎo)彈和目標(biāo)的距離小于等于10 km時(shí)進(jìn)入末制導(dǎo)，仿真結(jié)束。仿真結(jié)果見圖2～圖8。

圖2為導(dǎo)彈和目標(biāo)彈道，可見該制導(dǎo)律可將導(dǎo)彈導(dǎo)引至目標(biāo)附近。

圖2 導(dǎo)彈和目標(biāo)彈道曲線Fig.2 Trajectories of missile and target curve

圖3 和圖4分別為高度和彈道傾角曲線。可以看出，基于動(dòng)態(tài)逆的奇異攝動(dòng)制導(dǎo)律的本質(zhì)是使導(dǎo)彈跟隨降階的最優(yōu)值。

圖5為比能曲線。可以看出，在下落段，導(dǎo)彈能量幾乎沒有損失;在平穩(wěn)飛行段，能量亦平穩(wěn)減小，能夠較好地為末端攔截保存能量。

圖3 高度曲線Fig.3 Altitude curve

圖4 彈道傾角曲線Fig.4 Flight path angle curve

圖5 比能曲線Fig.5 Specific energy curve

圖6 為指令過載曲線。可以看出，在導(dǎo)彈下落的初始段，指令過載趨近于0，使該段能量損失減少到最小。

圖6 指令過載曲線Fig.6 Command overload curve

仿真算例中，用迎角作為導(dǎo)彈的控制量。由圖7所示的迎角曲線可看出，在短暫的自由下落結(jié)束后，迎角立即增大到最大值，以產(chǎn)生足夠的過載使導(dǎo)彈高度調(diào)整至降階最優(yōu)高度位置。

圖7 迎角曲線Fig.7 Angle of attack curve

導(dǎo)彈沿降階最優(yōu)高度飛行，形成了平滑的彈道，使大部分飛行時(shí)間內(nèi)由圖8所示的熱流率曲線沒有發(fā)生大的變化。在彈道末端，由于調(diào)整導(dǎo)彈飛行至目標(biāo)附近，產(chǎn)生了較大的指令過載，熱流率也有所上升。整個(gè)中制導(dǎo)段，最大熱流率為551.02 W/cm2，大大低于跳躍-滑翔彈道的最大熱流率值。通過調(diào)整動(dòng)態(tài)逆控制器中的系數(shù)ωh和ωγ，可獲得更小的最大熱流率，以達(dá)到工程應(yīng)用的要求。

圖8 最大熱流率曲線Fig.8 Maximum heating rate

4 結(jié)束語

本文針對(duì)高超聲速導(dǎo)彈的無約束中制導(dǎo)彈道產(chǎn)生熱流率過大的問題，用奇異攝動(dòng)和動(dòng)態(tài)逆方法推導(dǎo)了一種最大末速中制導(dǎo)律。用動(dòng)態(tài)逆方法跟蹤由降階系統(tǒng)計(jì)算得出的最優(yōu)高度和由第一邊界層系統(tǒng)計(jì)算得出的最優(yōu)彈道傾角，得到了一個(gè)閉環(huán)近最優(yōu)中制導(dǎo)律。仿真結(jié)果可見，本文推導(dǎo)的中制導(dǎo)律除可將導(dǎo)彈導(dǎo)引至目標(biāo)附近，并為末端攔截保存能量外，產(chǎn)生了占中制導(dǎo)絕大部分時(shí)間的平滑彈道，使熱流率曲線在該段時(shí)間內(nèi)沒有突變，并且與無熱流約束的最優(yōu)滑行彈道相比，雖然保存能量能力稍弱，但大大減小了最大熱流率，符合工程應(yīng)用的要求。

［1］ 張麗靜，劉東升，于存貴，等.高超聲速飛行器［J］.航空兵器，2010，(2):13-16.

［2］ Naidu D S，Calise A.Singular perturbations and time scales in guidance，navigation and control of aerospace system［R］.AIAA-95-3321，1995.

［3］ Cheng V H L，Gupta N K.Advance midcourse guidance for air-to-air missile ［J］.Journal of Guidance，1986，9(2):135-142.

［4］ Sheu D L，Chen Y M，Chang Y J.Optimal glide for maximum range［R］.AIAA-98-4462，1998.

［5］ 喬清青，陳萬春，李佳峰.高超聲速飛行器最大航程滑行奇異攝動(dòng)導(dǎo)引律［C］//臨近空間飛行器及防御技術(shù)發(fā)展論壇文集.上海:中國航天科技集團(tuán)公司第八研究院第八設(shè)計(jì)部，2010:119-127.

［6］ Morimoto H，Chuang C H.Minimum-fuel trajectory along entire flight profile for a hypersonic vehicle with constraint［R］.AIAA-98-4122，1998.

［7］ Yu Wenbin，Chen Wanchun.Guidance scheme for glide range maximization of a hypersonic vehicle［R］.AIAA-2011-6714，2011.

［8］ Ardema M D，Rajan N.Slow and fast state variables for three-dimensional flight dynamics［J］.Journal of Guidance，1985，8(4):532-535.

［9］ 喬清青，陳萬春.基于動(dòng)態(tài)逆的空空導(dǎo)彈奇異攝動(dòng)中制導(dǎo)律［J］.北京航空航天大學(xué)學(xué)報(bào)，2011，37(11):1365-1371.