崔正達, 魏明英, 李運遷

(1. 北京電子工程總體研究所, 北京 100854; 2. 北京仿真中心, 北京 100083)

0 引 言

高超聲速滑翔飛行器通過升力增程的方式在高度20～100 km的臨近空間飛行,該空域空氣稀薄,飛行器受到的阻力較小,但仍具備機動能力,且受制于地球曲率,雷達難以遠距離探測和預警[1]。諸多優(yōu)勢使得高超聲速飛行器一經(jīng)問世就是學者們研究的熱點,被譽為航空工業(yè)皇冠上的明珠。高超聲速滑翔飛行器運動方程形式復雜、非線性特性強、再入過程對控制量高度敏感、過程約束多,受過載、熱流密度等約束影響較大、對終端狀態(tài)要求高等限制因素使得其軌跡設計成為該領域的難點和熱點之一。

為了在諸多約束條件下規(guī)劃航跡,精準達成作戰(zhàn)目標,目前主要有3種技術途徑:標稱軌跡法、預測校正法和前兩者的混合制導方法[2]。標稱軌跡法是將上述問題轉(zhuǎn)化為多約束非線性規(guī)劃問題,經(jīng)過大量的離線運算或者在線優(yōu)化進行求解,得到參考飛行軌跡。標稱軌跡法在20世紀70年代首先被發(fā)展并應用在航天飛機的再入返回過程之中[3-4],1977年Page和Rogers針對大氣層外飛行器再入問題提出了基于三次曲線的彈道參考制導律[5],該制導率可以滿足終端角度約束,類似的還有基于圓弧制導軌跡的方法[6]。但大量的仿真實驗表明該類制導律精度欠佳[7]。1979年Harpold和Graves提出了經(jīng)典的基于速度多項式的阻力剖面再入制導方法[8],研究者在此基礎上基于不同的關注方向,發(fā)展出了表征飛行器運動特性的很多種廣義飛行剖面,在飛行任務中只需跟蹤迭代生成的參考軌跡即可達到預定效果。較常見的有升阻比-能量剖面[9-10]、阻力-能量剖面[11-12]、阻力-傾側(cè)角指令剖面[13]、阻力-速度剖面[8]、攻角-速度剖面[14-15]、動壓-高度剖面[16]、高度-射程剖面[17]等。

預測校正法的核心思想是利用彈載計算機的運算能力在線預測飛行軌跡以及終端狀態(tài),并求解出其與期望值的偏差,然后校正制導指令不斷消除預測飛行軌跡與理想值的偏差。該方法相較于標稱軌跡法具有更高的魯棒性,而根據(jù)軌跡預測的方式不同,預測校正法又可以分為數(shù)值解[18-22]和解析解法[23-26]。解析預測制導需要常系數(shù)假設[24]或者將軌跡調(diào)制到特定形式獲得近似解析解預測[17];數(shù)值預測的基本思路是先由預先指定的剖面求得控制指令序列,再代入彈道微分方程求解數(shù)值積分對終端狀態(tài)進行預測,該方法在線計算量很大,對于彈載計算機的計算能力有較大的依賴。

上述兩類典型制導方法各存在優(yōu)點和不足,有學者嘗試將二者互相取長補短形成混合制導方法兼顧標稱軌跡法的快速性以及預測校正法的魯棒性優(yōu)點:如Hu等[27]就將滑翔段再分成平衡段和線性段,分段采用不同的方法以獲得二者的優(yōu)勢。王青等[28]、Wang等[29]也是采用分段的思想將兩種制導方法混合,兼顧速度及精度。

多飛行器協(xié)同制導概念出現(xiàn)以后,時間控制需求迫切,針對滑翔段的時間預測算法也相繼出現(xiàn)[30-32]。但是滑翔段的優(yōu)化算法由于使用平衡滑翔或定阻力系數(shù)假設,不再適用于復雜變化的下壓段。

快速下壓彈道的主要任務是在復雜約束下將速度、高度和彈道傾角等飛行器運動參數(shù)控制到指定值,其間要經(jīng)歷非常強的氣動非線性影響。國內(nèi)外在下壓段的研究以彈道整形制導律[33-36]和帶落角約束的制導律[37-39]為主。這些經(jīng)典理論提出的時間都比較早,時間不是主要的受控目標。上述制導律可以滿足落角和落點約束,但無法對遭遇時間進行準確的估計和控制,不能滿足時間協(xié)同的新興需求。

基于上述情況,本文提出一種考慮阻力系數(shù)時變的下壓段時間預測方法,設計隨攻角和馬赫數(shù)變化的阻力系數(shù)模型,并簡化了飛行器運動方程以精確快速預報導彈運動參數(shù)和時間。本文按如下順序編排:首先建立飛行器下壓段數(shù)學模型;然后基于解析理論給出以剩余射程為自變量、滿足終端約束的高度-射程剖面表達式,將空氣密度隨高度變化規(guī)律線化為剩余射程的多項式函數(shù),并且結合氣動擬合結果,在阻力系數(shù)表達式中引入攻角和馬赫數(shù)的影響項,基于上述解析理論推導降階彈道微分方程,得到了以剩余飛行距離為自變量的一階微分方程,該方程可以通過數(shù)值積分快速求解;最后以典型彈道為例給出對比傳統(tǒng)方法的數(shù)值仿真結果以及對結果的討論。

1 滑翔飛行器數(shù)學模型

高超聲速飛行器在航跡末段快速俯沖,其飛行時間較短、飛行高度和射程較低、縱側(cè)向耦合不深,因此可獨立設計兩個維度的軌跡并假設地球為不旋轉(zhuǎn)的勻質(zhì)圓球,得到攔截器下壓段縱平面內(nèi)動力學方程:

(1)

式中:V為飛行器相對于地面的速度標量;r為飛行器質(zhì)心到地心的距離;θ為當?shù)貜椀纼A角,定義為速度方向與當?shù)厮矫娴膴A角,向上為正;m為飛行器質(zhì)量;g為當?shù)匾铀俣?Y和D為飛行器所受的升力及阻力;CD為飛行器阻力系數(shù);CL為飛行器升力系數(shù),該系數(shù)可認為是攻角和馬赫數(shù)的表達式;Sref為飛行器參考面積,由于只考慮縱平面內(nèi)的運動所以控制量選為攻角α。

考慮氣動力時,采用如下的指數(shù)大氣密度模型:

ρ=ρ0e-β h

其中，β=1.406 4×10-4m-1;ρ0=1.225 kg·m-2。

定義從當前位置到終端位置間剩余距離在地面上的投影長度為射程,記為RL。

2 高度射程剖面設計

快速下壓段需要在約束范圍內(nèi)將高度和彈道傾角規(guī)劃到指定區(qū)域,所以針對下降段距離近、彈道形狀復雜、約束條件嚴苛的應用背景下拓展了文獻[17]提出的高度-射程剖面表達式,在彈道震蕩平緩預設下定義飛行軌跡為剩余射程的六次函數(shù):

(2)

其中,

a1RL+a0

可將下壓段非線性特性最強的空氣密度簡化為航程的多項式形式便于解析降階,有:

ρ=ρ0e-β h=ρ0F(RL)

(3)

式(2)所確定的彈道形狀需要確定a0～a6共7個彈道參數(shù),選取初始高度h0、初始飛行路徑角θ0、1/3處高度h1、2/3處高度h2、終端高度hf、終端路徑角θf以及彈道h1點處路徑角θH。線性解算矩陣同文獻[17]所述,選取次高點傾角θH的原因是該值可以有效控制過載分布,充分利用飛行器機動能力。若不如此拓展,根據(jù)高度射程表達式的數(shù)學特性,會出現(xiàn)高空可用過載能力小但需用過載能力大的情況,不適合實際工程應用。

下壓段中當?shù)貜椀纼A角θ、傾側(cè)角σ均為小量,則飛行路徑角的解析解為

(4)

(5)

由飛行器運動方程可知彈道傾角及其變化率可以與飛行器飛行高度變化率及升力聯(lián)系起來:

(6)

(7)

聯(lián)立式(1)、式(2)、式(5)和式(7)可得到升力Y的簡化表達式:

(8)

式中:

3 考慮攻角及馬赫數(shù)變化的阻力系數(shù)模型

現(xiàn)有彈道參數(shù)解析預測方法在推導彈道解析解的時候,為簡化推導過程習慣假定阻力系數(shù)為常值,這樣的假設對于高度和速度變化相對穩(wěn)定的滑翔段是合適的。但是,當阻力系數(shù)的兩個主要影響因素即馬赫數(shù)和攻角發(fā)生改變時,該假設與真實的偏差會對結果產(chǎn)生較大影響,影響預報精度。設氣動力系數(shù)CL和CD可表示為如下形式的攻角、馬赫數(shù)函數(shù):

(9)

用攻角跟蹤彈道所需升力Y時有如下關系:

(10)

當攻角α=0時不產(chǎn)生升力,即Y≈0,由式(9)得CL0≈0。

然后將攻角表達式(10)代入式(9),考慮攻角α對阻力系數(shù)CD的影響得到:

(11)

聯(lián)立升力值表達式(8):

(12)

分析上述表達式,式(12)中第1部分即常值部分,是傳統(tǒng)預測方法所考慮的部分;第2部分是馬赫數(shù)變化帶來的阻力系數(shù)的變化量;第3部分是改變彈道形狀和平衡重力所需的升力,進而誘導產(chǎn)生的額外阻力。

根據(jù)下壓彈道常見高度范圍,得到以剩余射程為自變量的阻力系數(shù)拓展表達式:

(13)

根據(jù)動力學方程可得

(14)

有:

(15)

由式(15)聯(lián)立式(13)可得

(16)

(17)

整理得

(18)

聯(lián)立下壓段高度表達式(3):

(19)

式(19)就是解析理論簡化得到的一階微分方程,數(shù)值積分可以得到準確的速度隨射程變化規(guī)律。方程中第1項兼顧了傳統(tǒng)常值假設下所考慮的情形,在其基礎上第2項是由馬赫數(shù)變化主導的影響項,第3項則是由攻角主導的影響項。第4項為重力勢能和動能的轉(zhuǎn)換項。為了簡化表述定義如下常數(shù):

式(19)可以寫為

(20)

觀察以RL為自變量,V為因變量的微分方程,該速度-剩余射程的方程在彈道解析理論下變?yōu)橐浑A微分方程,可以用數(shù)值解法來預報其運動狀態(tài),從而以剩余飛行距離RL為自變量對飛行器速度和時間進行估計。若省略后3項,則該方程退化為原始文獻形式。

4 仿真分析

4.1 3種仿真算例

根據(jù)文獻[30],設置仿真初始條件為:經(jīng)緯度λ0=φ0=0°,剩余距離RL為210 km,初始高度為h0=35 km,速度V0=2 000 m/s,速度傾角θ0=0°。仿真計算機所用的CPU型號為Intel(R) Core i5-4590,主頻3.29 GHz。軟件為Matlab 2012b。

在本節(jié)中,改變彈道參數(shù)h1、h2、θH可得到不同的彈道形狀,模擬3種不同的彈道形式:第1種是滑翔段飛行彈道,第2種是高度大范圍變化的末段俯沖下壓彈道,第3種是先上揚再快速下壓的末段俯沖突防彈道。對這3種典型彈道分別比較文獻[17]所應用的常值阻力假設、本文提出的簡化微分方程和彈道數(shù)值積分得到的預報結果,得到的仿真曲線如圖1～圖6所示。

圖1 近似平衡滑翔彈道Fig.1 Trajectory in quasi-equilibrium glide phase

圖2 近似平衡滑翔速度變化曲線Fig.2 Velocity change curve in quasi-equilibrium glide phase

圖3 俯沖下壓彈道Fig.3 Trajectory in dive phase

圖4 俯沖下壓速度變化曲線Fig.4 Velocity change curve in dive phase

圖5 俯沖突防彈道Fig.5 Trajectory in penetration phase

圖6 俯沖突防速度變化曲線Fig.6 Velocity in penetration phase

在滑翔段速度變化線性下降,由常值阻力系數(shù)假設得到的速度變化與本文提出的預報方法以及彈道仿真曲線十分接近,因為在滑翔段本文提出的方法拓展量皆為小量,方法退化為采用常值阻力假設的解析方法;在俯沖下壓段由于產(chǎn)生了負攻角以及較大的彈道傾角,在重力和空氣動力的共同作用下速度的變化規(guī)律呈現(xiàn)出了很強的非線性,本文提出的方法在設定阻力系數(shù)模型的時候考慮的因素較多,能夠比常值阻力假設更好地揭示出真實變化規(guī)律,故本文提出的拓展方法所預報的速度與仿真結果更為相近,預報誤差較小。另外,僅在縱平面內(nèi)借助重力勢能和動能之間的轉(zhuǎn)換規(guī)劃彈道形狀可以減小協(xié)同過程中的能量消耗。彈道仿真也說明彈道形狀、起始-終止的高度落差對導彈的抵達時間以及終端速度產(chǎn)生較大影響。

4.2 比較提出方法的計算時間效率

準確預報速度變化可以為抵達時間準確估計打下堅實基礎,對速度-剩余射程進行一階積分,以3種飛行器典型彈道為例對抵達時間進行預測,并記錄系統(tǒng)計算耗時,得到飛行時間估計結果和計算耗時如表1～表3所示。

表1 近似平衡滑翔段計算效率比較

表2 俯沖下壓段計算效率比較

表3 俯沖突防段計算效率比較Table 3 Computational efficiency comparison in penetration phase s

對比上述結果可知,文獻[17]預報的結果在下壓段和俯沖段與彈道仿真的區(qū)別比較大,常值假設無法適應,會出現(xiàn)5～10 s的預報誤差(相對誤差4.5%～9%),而簡化解析方程僅有不到1 s的預報誤差(相對誤差1%～0.9%)。通過與另外兩種預報方法進行綜合比對,本文所提方法可以在不顯著提升計算時間復雜度的同時顯著提升預報精度。

4.3 蒙特卡羅仿真結果

為驗證方法有效性,對各情況中氣動參數(shù)加入10%的拉偏量,各仿真1 000次得到的預測偏差如圖7～圖12所示。

圖7 平衡滑翔-常值阻力系數(shù)預測誤差分布Fig.7 Estimated error distribution in quasi-equilibrium glide phase with constant drag coefficient

圖8 平衡滑翔-提出方法預測誤差分布Fig.8 Estimated error distribution in quasi-equilibrium glide phase with proposed method

圖9 俯沖下壓-常值阻力系數(shù)預測誤差分布Fig.9 Estimated error distribution in dive phase with constant drag coefficient

圖10 俯沖下壓-提出方法預測誤差分布Fig.10 Estimated error distribution in dive phase with proposed method

圖11 俯沖突防-常值阻力系數(shù)預測誤差分布Fig.11 Estimated error distribution in penetration phase with constant drag coefficient

圖12 俯沖突防-提出方法預測誤差分布Fig.12 Estimated error distribution in penetration phase with proposed method

仿真結果表明,采用本文提出的方法可以有效地在各種彈道形狀下估計遭遇時間,并且估計散布小于常值阻力假設。常值假設標準差為0.421 92,本文提出的方法估計標準差為0.222 49。

圖13 動壓預測曲線對比Fig.13 Comparison of dynamic pressure estimation curve

本文采取的剖面形式可準確解析飛行器所在位置的大氣密度,因此對速度的精準預報提升了動壓的解析預測精度。動壓的精準預報可以從可用過載以及最大動壓兩個方面為優(yōu)化算法判斷約束條件提供輕便準確的模型,加快優(yōu)化速度及精度,有利于在彈道規(guī)劃過程中更大程度地發(fā)揮飛行器過載能力。

5 結 論

本文深入分析高超聲速滑翔飛行器下壓段速度變化規(guī)律,對精度影響最大的兩個因素,即攻角、馬赫數(shù),分別在阻力系數(shù)表達式中進行拓展,給出了新的微分方程形式。基于解析解理論對復雜的彈道微分方程進行簡化,得到僅有一階的簡化微分方程,該方程可以通過數(shù)值積分快速求解。通過典型彈道仿真表明,該方法可將俯沖下壓段的時間預報精度從10 s提高到1 s左右,同時不顯著提升計算復雜度,實現(xiàn)滑翔飛行器俯沖段時間的精準、快速預報,提高全彈道協(xié)同規(guī)劃能力。