喬清青,陳萬春
(北京航空航天大學宇航學院,北京100191)
高超聲速導彈由于具有飛行高度高、速度快、機動性好等特點,能有效縮短對目標的反應時間,具有較高的突防成功率[1]。高超聲速導彈飛行距離較遠,需要一種中制導律將導彈導引至目標附近,為末端作戰保存能量,滿足各種約束條件,并完成中末制導的平滑交接班。
導彈的中制導問題可歸結為一個非線性兩點邊值問題。由于導彈動力學系統的高度非線性,兩點邊值問題無法求出解析解。奇異攝動是解決復雜非線性最優控制問題的一種經典方法,通過抑制動力學方程中的小參數使系統階數降低,可容易地獲得相應的降階解,再計算小參數獲得邊界層解,以修正降階解。使用奇異攝動可得到性能較理想的近似解,并可根據需要獲得開環或閉環解。奇異攝動方法在航空航天領域得到了廣泛的應用。由于飛行器動力學方程中小參數不明顯,因此通常使用“強迫奇異攝動”方法,即先根據估算的各變量的時間常數以及應用經驗區分變化快慢不同的變量,再在快變量前插入小參數ε(0<ε<1)構成奇異攝動系統。令小參數為零,即假設較快變量能在瞬間達到平衡,從而使系統階數降低,并求解出降階解,然后通過邊界層修正計算考慮被忽略的快變量及其邊界條件。用奇異攝動方法求解飛行器制導非線性兩點邊值問題,可獲得算法較簡單,并可在彈上實時解算的近最優制導律[2-3]。
文獻[4]針對超聲速飛行器的最大航程滑行問題,將動力學方程分成三個層次,用奇異攝動方法求解了一種最大航程閉環導引律。文獻[5]考慮了球形地球上的飛行器質點動力學模型,推導了奇異攝動最大航程制導律,給出了該算法與用序列二次規劃(SQP)算法解算出的最優彈道和熱流率曲線。無熱流約束的彈道均為跳躍-滑翔式彈道,飛行中產生的熱流率也較大。
針對上述問題,本文利用文獻[6]中的高超聲速飛行器的氣動模型,用奇異攝動和動態逆方法為高超聲速導彈推導了一種使熱流率較低的最大末速中制導律。

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


式中,CD0和CDt分別為零升阻力系數和升致阻力系數。動壓定義為q=ρV2/2。各氣動力系數見文獻[4]。導彈迎角為:

式中,m,g,S分別為導彈的質量、當地重力加速度和參考面積,“~”表示變量的實際值。
令導彈飛行時產生的最大熱流率為:

式中,熱流率單位為W/cm2。
文獻[6]用奇異攝動方法為高超聲速飛行器推導了一種最遠航程制導律,并用其與配點法和SQP方法解算出的最優彈道進行比較,證明了高超聲速導彈的無約束彈道為跳躍-滑翔彈道。當導彈飛行至較高高度時,指令過載很小,甚至趨近于0。跳躍式彈道可很好地為導彈保存能量,但產生的熱流率較大,不適于工程實際應用。
根據文獻[7]中的結論,要減小熱流率,就應使飛行彈道盡量平滑。用一些尋優方法,經過反復迭代運算可以得到較理想的結果,但其所需運算時間長,無法在飛行器上實時進行。另外,目標的運動,以及導彈在實際飛行中受到的各種干擾,均會導致飛行狀態發生改變,從而無法滿足跟蹤最優值所需的條件。
由于高超聲速導彈飛行時間較長,需要在中制導段盡量保存能量,以保證末端攔截有足夠的能量。問題變成使微分方程式(1)~式(4)所描述的系統從初始狀態轉移到預測攔截點在笛卡爾坐標系中的末端位置(rf,hf),使得末端比能最大。假設目標以當前速度作勻速直線運動,則經過剩余飛行時間tgo后,其所處位置作為預測攔截點。剩余飛行時間估算如下:

式中,Rtm為導彈與目標的距離;Vc為接近速度;a為導彈沿視線方向的平均加速度,認為與導彈軸向加速度相等。
用文獻[8]中介紹的方法估算各狀態變量的時間常數。根據其大小可將變量分成:慢變量(r,E)、中間變量h和快變量γ。在中間變量和快變量前分別插入小參數 ε1(0<ε1<1)和 ε2(0<ε2<1),可得到導彈的質心動力學奇異攝動模型為:

為使末端速度最大,可令性能指標為:J=-Ef。
最慢變量構成了降階系統,由降階系統求解出的解稱為外解。令ε1=ε2=0,得到γ=0,n=1-V2/(R0+h)。帶入式(7)~式(10),得到降階系統為:,其中:D=D0+K[1-V2/(R0+h)]2/q;高度h為零階系統的控制量。Hamilton函數寫成:H=λrVR0/(R0+h)-λEDV。由伴隨方程,可知 λr為常數。由橫截條件 λE(tf)= -1和H(tf)=0,得:

系統中各變量與時間無關,因此H(t)=0,可求得:

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

式(13)為h和E的隱函數。為加快彈上運算速度,可先將不同比能下的?D/?h與高度h的關系列成表格,這樣通過插值即可求出與某指定能量下對應的最優高度h*。根據第1節中導彈氣動力系數,計算出的比能與降階最優高度關系如圖1所示。

圖1 降階系統最優高度與比能的關系Fig.1 Relationship between reduced order optimal altitude and specific energy
本節中的V,D,h均為降階值,下文中用上標“o”表示外解中的變量。降階系統中忽略了高度變量的邊界值。為使導彈能夠達到預定位置,對降階最優高度作如下修正:

式中,kh為待選系數;rtm為兩者之間的水平距離。
在飛行器質量不變和一定速度條件下,飛行器受到的零升阻力隨海拔高度的增加而減小,由相同過載產生的阻力隨海拔高度的增加而增加。因此,當指令過載不為0時,存在一個使阻力達到極小值的高度。
從物理上看,式(13)的意義在于:在某能量狀態下,在降階最優高度上作水平飛行時,飛行器受到的阻力最小。為得到一較平滑的彈道,邊界層修正得到的系統控制量應使導彈能在降階最優高度上飛行。
令 ε1=1,ε2=0,則 n=cos γ[1 - V2/(R0+h)],得到的第一邊界層動態由高度動力學方程決定。該層的控制量為彈道傾角γ。用動態逆方法對該層進行求解,第一邊界層的系統方程和輸出方程為:

該系統并非仿射非線性系統,因此將控制量改設為:


式中,ωh為待定系數。系統的控制為:

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

對導彈系統的時間尺度分析可知,速度和高度變量在相同的時間尺度上,因此該層中的速度變量可看作當前導彈速度V。考察式(15)所示的一階系統,其性能主要取決于h的時間常數。
下文中用上標“i1”表示第一邊界層中的變量。該層計算中忽略了控制量彈道傾角的邊界條件。為了使導彈在中制導末端,導引頭能夠順利截獲目標,需對最優進行修正:

其中:

令ε1=ε2=1,可得到第二邊界層動力學方程,為系統方程完整形式。其控制變量為指令過載n。第二邊界層系統動態主要由γ的動力學方程決定,控制量為n。其系統方程和輸出方程為:

第二邊界層系統為原系統中最快的一層,認為其中的變量均可取為當前值。可看出,上式所示系統是一個仿射非線性系統,可直接構造逆系統并求解控制律。求得控制量為:

式中,ωγ為待定系數。
通過式(13)得到降階最優后,用式(16)和式(17)即可計算得到導彈的指令過載。
根據動態逆理論,式(16)和式(17)中的待定系數ωh和ωγ應分別為高度變量和彈道傾角變量的時間常數的倒數。用文獻[4]介紹的方法可估算變量的時間常數,但該方法依賴人工選取的典型值,因此估算出的時間常數在某區間內均可取值。文獻[9]中介紹了如何根據二階系統階躍響應分析方法選取合適的時間常數值的方法,然而由于高超聲速導彈飛行速度快,可達到很高的高度,導致其受到的阻力對指令過載的變化相當敏感。因此,當導彈處于很大的高度時,其指令過載趨近于0,這樣才不會使飛行器能量迅速損失。
對于本文中的例子,導彈在初始一段時間內從較大的高度下落,需對初始下降段的降階最優高度和第一邊界層最優彈道傾角進行修正:當h=92 008 m時,將導彈以當前狀態作自由下落的彈道高度作為降階最優高度;將導彈以當前狀態作自由下落的彈道傾角作為第一邊界層最優彈道傾角,以保證導彈在最初下落段指令過載趨近于0。
本文用Matlab軟件對制導律進行了仿真驗證。導彈的初始狀態為:r0=0 m,h0=112 080 m,V0=5164.9 m/s。設目標初始狀態 ht(t0)=10 km,rt(t0)=5 500 km。目標以Ma=2的速度向導彈迎面勻速等高飛行。仿真過程中,令迎角α≤25°,導彈和目標的距離小于等于10 km時進入末制導,仿真結束。仿真結果見圖2~圖8。
圖2為導彈和目標彈道,可見該制導律可將導彈導引至目標附近。

圖2 導彈和目標彈道曲線Fig.2 Trajectories of missile and target curve
圖3 和圖4分別為高度和彈道傾角曲線。可以看出,基于動態逆的奇異攝動制導律的本質是使導彈跟隨降階的最優值。
圖5為比能曲線。可以看出,在下落段,導彈能量幾乎沒有損失;在平穩飛行段,能量亦平穩減小,能夠較好地為末端攔截保存能量。

圖3 高度曲線Fig.3 Altitude curve

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

圖5 比能曲線Fig.5 Specific energy curve
圖6 為指令過載曲線。可以看出,在導彈下落的初始段,指令過載趨近于0,使該段能量損失減少到最小。

圖6 指令過載曲線Fig.6 Command overload curve
仿真算例中,用迎角作為導彈的控制量。由圖7所示的迎角曲線可看出,在短暫的自由下落結束后,迎角立即增大到最大值,以產生足夠的過載使導彈高度調整至降階最優高度位置。

圖7 迎角曲線Fig.7 Angle of attack curve
導彈沿降階最優高度飛行,形成了平滑的彈道,使大部分飛行時間內由圖8所示的熱流率曲線沒有發生大的變化。在彈道末端,由于調整導彈飛行至目標附近,產生了較大的指令過載,熱流率也有所上升。整個中制導段,最大熱流率為551.02 W/cm2,大大低于跳躍-滑翔彈道的最大熱流率值。通過調整動態逆控制器中的系數ωh和ωγ,可獲得更小的最大熱流率,以達到工程應用的要求。

圖8 最大熱流率曲線Fig.8 Maximum heating rate
本文針對高超聲速導彈的無約束中制導彈道產生熱流率過大的問題,用奇異攝動和動態逆方法推導了一種最大末速中制導律。用動態逆方法跟蹤由降階系統計算得出的最優高度和由第一邊界層系統計算得出的最優彈道傾角,得到了一個閉環近最優中制導律。仿真結果可見,本文推導的中制導律除可將導彈導引至目標附近,并為末端攔截保存能量外,產生了占中制導絕大部分時間的平滑彈道,使熱流率曲線在該段時間內沒有突變,并且與無熱流約束的最優滑行彈道相比,雖然保存能量能力稍弱,但大大減小了最大熱流率,符合工程應用的要求。
[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] 喬清青,陳萬春,李佳峰.高超聲速飛行器最大航程滑行奇異攝動導引律[C]//臨近空間飛行器及防御技術發展論壇文集.上海:中國航天科技集團公司第八研究院第八設計部,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] 喬清青,陳萬春.基于動態逆的空空導彈奇異攝動中制導律[J].北京航空航天大學學報,2011,37(11):1365-1371.