考慮阻力系數(shù)時(shí)變的下壓段半解析時(shí)間預(yù)測(cè)方法

崔正達(dá), 魏明英, 李運(yùn)遷

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

0 引 言

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

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

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

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

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

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

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

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

高超聲速飛行器在航跡末段快速俯沖,其飛行時(shí)間較短、飛行高度和射程較低、縱側(cè)向耦合不深,因此可獨(dú)立設(shè)計(jì)兩個(gè)維度的軌跡并假設(shè)地球?yàn)椴恍D(zhuǎn)的勻質(zhì)圓球,得到攔截器下壓段縱平面內(nèi)動(dòng)力學(xué)方程:

(1)

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

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

ρ=ρ0e-β h

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

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

2 高度射程剖面設(shè)計(jì)

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

(2)

其中,

a1RL+a0

可將下壓段非線性特性最強(qiáng)的空氣密度簡(jiǎn)化為航程的多項(xiàng)式形式便于解析降階,有:

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

(3)

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

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

(4)

(5)

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

(6)

(7)

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

(8)

式中:

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

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

(9)

用攻角跟蹤彈道所需升力Y時(shí)有如下關(guān)系:

(10)

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

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

(11)

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

(12)

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

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

(13)

根據(jù)動(dòng)力學(xué)方程可得

(14)

有:

(15)

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

(16)

(17)

整理得

(18)

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

(19)

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

式(19)可以寫(xiě)為

(20)

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

4 仿真分析

4.1 3種仿真算例

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

在本節(jié)中,改變彈道參數(shù)h1、h2、θH可得到不同的彈道形狀,模擬3種不同的彈道形式:第1種是滑翔段飛行彈道,第2種是高度大范圍變化的末段俯沖下壓彈道,第3種是先上揚(yáng)再快速下壓的末段俯沖突防彈道。對(duì)這3種典型彈道分別比較文獻(xiàn)[17]所應(yīng)用的常值阻力假設(shè)、本文提出的簡(jiǎn)化微分方程和彈道數(shù)值積分得到的預(yù)報(bào)結(jié)果,得到的仿真曲線如圖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ù)假設(shè)得到的速度變化與本文提出的預(yù)報(bào)方法以及彈道仿真曲線十分接近,因?yàn)樵诨瓒伪疚奶岢龅姆椒ㄍ卣沽拷詾樾×?方法退化為采用常值阻力假設(shè)的解析方法;在俯沖下壓段由于產(chǎn)生了負(fù)攻角以及較大的彈道傾角,在重力和空氣動(dòng)力的共同作用下速度的變化規(guī)律呈現(xiàn)出了很強(qiáng)的非線性,本文提出的方法在設(shè)定阻力系數(shù)模型的時(shí)候考慮的因素較多,能夠比常值阻力假設(shè)更好地揭示出真實(shí)變化規(guī)律,故本文提出的拓展方法所預(yù)報(bào)的速度與仿真結(jié)果更為相近,預(yù)報(bào)誤差較小。另外,僅在縱平面內(nèi)借助重力勢(shì)能和動(dòng)能之間的轉(zhuǎn)換規(guī)劃彈道形狀可以減小協(xié)同過(guò)程中的能量消耗。彈道仿真也說(shuō)明彈道形狀、起始-終止的高度落差對(duì)導(dǎo)彈的抵達(dá)時(shí)間以及終端速度產(chǎn)生較大影響。

4.2 比較提出方法的計(jì)算時(shí)間效率

準(zhǔn)確預(yù)報(bào)速度變化可以為抵達(dá)時(shí)間準(zhǔn)確估計(jì)打下堅(jiān)實(shí)基礎(chǔ),對(duì)速度-剩余射程進(jìn)行一階積分,以3種飛行器典型彈道為例對(duì)抵達(dá)時(shí)間進(jìn)行預(yù)測(cè),并記錄系統(tǒng)計(jì)算耗時(shí),得到飛行時(shí)間估計(jì)結(jié)果和計(jì)算耗時(shí)如表1～表3所示。

表1 近似平衡滑翔段計(jì)算效率比較

表2 俯沖下壓段計(jì)算效率比較

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

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

4.3 蒙特卡羅仿真結(jié)果

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

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

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

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

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

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

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

仿真結(jié)果表明,采用本文提出的方法可以有效地在各種彈道形狀下估計(jì)遭遇時(shí)間,并且估計(jì)散布小于常值阻力假設(shè)。常值假設(shè)標(biāo)準(zhǔn)差為0.421 92,本文提出的方法估計(jì)標(biāo)準(zhǔn)差為0.222 49。

圖13 動(dòng)壓預(yù)測(cè)曲線對(duì)比Fig.13 Comparison of dynamic pressure estimation curve

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

5 結(jié) 論

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