李天任，雷建長(zhǎng)，王宇航，周 華，黃 佩

0 引 言

近年來(lái)，國(guó)內(nèi)外學(xué)者對(duì)高超聲速飛行器再入軌跡優(yōu)化與制導(dǎo)開(kāi)展了大量研究：一類是傳統(tǒng)的返回艙式低升阻比再入飛行器，常采用傾側(cè)角反轉(zhuǎn)的軌跡優(yōu)化和制導(dǎo)方法；另一類是近年來(lái)發(fā)展的高升阻比滑翔飛行器，一般采用攻角、傾側(cè)角作為控制變量，通過(guò)設(shè)計(jì)兩者在再入過(guò)程中隨時(shí)間、速度、高度的變化規(guī)律來(lái)滿足過(guò)程、終端狀態(tài)約束，同時(shí)達(dá)到某一性能指標(biāo)最優(yōu)。

目前，離線進(jìn)行優(yōu)化的偽譜法是較為成熟的再入體軌跡設(shè)計(jì)方法。偽譜法通過(guò)在整條軌跡的配點(diǎn)處離散，將再入動(dòng)力學(xué)微分方程轉(zhuǎn)化為代數(shù)方程，最后采用序列二次規(guī)劃（Sequential Quadratic Programming，SQP）等方法求解非線性規(guī)劃問(wèn)題。其在單機(jī)上進(jìn)行一次軌跡設(shè)計(jì)的計(jì)算時(shí)間較長(zhǎng)，一般用于離線生成標(biāo)稱軌跡。

為了發(fā)揮高升阻比飛行器的在線機(jī)動(dòng)能力，水尊師[1]和肖紅[2]等人發(fā)展了基于航路點(diǎn)的局部預(yù)測(cè)校正制導(dǎo)方法，但其離線設(shè)計(jì)的航路點(diǎn)是固定的，削弱了飛行器軌跡的變更、再生成能力?；贖arpold[3]的航天飛機(jī)再入軌跡設(shè)計(jì)思路，周浩[4]、李強(qiáng)[5]、劉運(yùn)鵬[6]等人采用常值、分段線性等攻角設(shè)計(jì)方案，用于軌跡的初步設(shè)計(jì)和在線預(yù)測(cè)。另一方面，胡錦川[7]的平穩(wěn)滑翔軌跡以及張冉[8]均要求攻角曲線連續(xù)且可導(dǎo)，在使軌跡平穩(wěn)的同時(shí)降低飛行器姿態(tài)控制器的設(shè)計(jì)要求。

下降段是再入飛行器進(jìn)入大氣層后的初始飛行階段，也被稱為引入段。為了在下降段軌跡的設(shè)計(jì)中簡(jiǎn)化控制量曲線的設(shè)計(jì)，同時(shí)增強(qiáng)控制量剖面的光滑度，易于與平衡滑翔段軌跡光滑對(duì)接，除了現(xiàn)有的常值、分段線性以及復(fù)合指數(shù)函數(shù)的剖面設(shè)計(jì)方案外，有必要引入新的函數(shù)形式。本文采用改進(jìn)的復(fù)合三角函數(shù)，其具有幅值可調(diào)、波峰的位置和峰形可變、導(dǎo)數(shù)連續(xù)且光滑等優(yōu)點(diǎn)，適合在控制量剖面設(shè)計(jì)中采用。

1 再入動(dòng)力學(xué)模型

針對(duì)無(wú)動(dòng)力再入飛行器，在位置系下建立動(dòng)力學(xué)模型，O-XYZ為地心固連坐標(biāo)系，將其繞Z軸正向旋轉(zhuǎn)經(jīng)度θ，再繞新的 Y軸旋轉(zhuǎn)緯度φ得到位置坐標(biāo)系O-xyz，如圖1所示。

圖1 三自由度運(yùn)動(dòng)參數(shù)示意Fig.1 3D Motion Parameters Diagram

采用圓球無(wú)旋地球假設(shè)，飛行器動(dòng)力學(xué)方程[9]為

式中 r，θ，φ，V，γ，ψ分別為地心距、經(jīng)度、緯度、速度、彈道傾角、航向角；g為重力加速度；m為飛行器質(zhì)量；σ為傾側(cè)角；L，D分別為升力和阻力，且，

式中 ρ為大氣密度；S為參考面積； CL， CD分別為升力系數(shù)、阻力系數(shù)。

m，S，CL, CD均參考國(guó)外滑翔飛行器CAV-H，m=907 kg，S=0.4839 m2。對(duì)升力系數(shù) CL和阻力系數(shù) CD進(jìn)行擬合，得到兩者對(duì)攻角α的關(guān)系式[10]：

在再入飛行器初始下降段，常采用零傾側(cè)角、大攻角拉起的方案，在飛行器狀態(tài)參數(shù)滿足一定條件時(shí)，切換到平衡滑翔階段，隨后采用平衡滑翔段的攻角、傾側(cè)角剖面。本文采用的切換條件為

式中1K為彈道傾角的權(quán)系數(shù)；ε為一小量。

采用下降段零傾側(cè)角、單波峰攻角剖面，平衡滑翔段的攻角剖面由擬平衡滑翔條件（Qusi Equilibrium Glide Condition，QEGC）推導(dǎo)得到。根據(jù)擬平衡滑翔條件：

根據(jù)當(dāng)前狀態(tài)已知的 V ,r,g，聯(lián)立式（2）、式（3）、式（6），得到擬平衡滑翔攻角 αQEGC為

式中 A，B為式（3）中攻角的一階項(xiàng)系數(shù)和常數(shù)項(xiàng)。

考慮到氣動(dòng)系數(shù)誤差，在 αQEGC上反饋彈道傾角的誤差項(xiàng) Δ α= K2(γ -γ*)，其中，γ*為平衡滑翔段跟蹤的彈道傾角， K2為反饋增益系數(shù)。由于本文研究的重點(diǎn)不是平衡滑翔段，可取γ*=-0.05°，從而有：

仿真結(jié)果表明，該方法能夠在考慮各項(xiàng)偏差的情況下生成擬平衡滑翔軌跡。

2 下降段控制量剖面設(shè)計(jì)方法

2.1 下降段控制量設(shè)計(jì)思路

下降段攻角剖面常采用最大攻角、分段線性等攻角方案[3～5]，在與滑翔段的對(duì)接處難以連續(xù)和可導(dǎo)，使軌跡平穩(wěn)性降低、同時(shí)對(duì)姿態(tài)控制器的設(shè)計(jì)提出了很高要求。此外，大量的彈道仿真表明不連續(xù)或不可導(dǎo)的控制量剖面常導(dǎo)致對(duì)接點(diǎn)過(guò)載和熱流峰值較大，難以滿足過(guò)程約束。

考慮在再入下降段采用改進(jìn)的復(fù)合三角函數(shù)，以增強(qiáng)對(duì)接點(diǎn)光滑性、降低最大過(guò)載和熱流峰值。

2.2 單波峰攻角剖面生成方法

再入飛行器下降到一定高度，采用如下攻角時(shí)序：

式中mα為飛行器最大可用攻角。根據(jù)CAV-H的氣動(dòng)系數(shù)表[10]，mα取20°，()Ft采用時(shí)間t的二次多項(xiàng)式函數(shù)：

通過(guò)調(diào)節(jié)多項(xiàng)式系數(shù) ai(i = 0 ,1,2)來(lái)調(diào)整波長(zhǎng)、波形和波峰位置，生成廣義正弦函數(shù)簇，見(jiàn)圖2。

圖2 不同波長(zhǎng)、波峰位置的三角函數(shù)簇Fig.2 Trigonometric Functions with Diffierent Wavelength and Peak Position

2.3 搜索過(guò)渡點(diǎn)光滑的攻角剖面

為了增強(qiáng)下降段與平衡滑翔段過(guò)渡點(diǎn)處的攻角剖面光滑性，需要優(yōu)化控制下降段攻角剖面形狀的ai(i = 0 ,1,2)3個(gè)參數(shù)。本文采用改進(jìn)的差分進(jìn)化算法（Differential Evolution，DE）[11]進(jìn)行參數(shù)尋優(yōu)。目標(biāo)函數(shù)見(jiàn)式（11）。

式中 j表示第j個(gè)積分區(qū)間，積分區(qū)間個(gè)數(shù) N = tf/h ，tf，h為積分結(jié)束時(shí)間和積分步長(zhǎng)。本文在Storn[11]的交叉操作之后增加拉回步驟，拉回不在解空間的新個(gè)體，以增強(qiáng)算法的魯棒性。

3 數(shù)值仿真及分析

3.1 攻角剖面光滑性驗(yàn)證

CAV-H 再入點(diǎn)參數(shù)取 H=70 km，V=7.5 km/s，γ=-2°。下降段采用復(fù)合三角函數(shù)的攻角剖面，平衡滑翔段采用擬平衡滑翔攻角，以控制量光滑性為目標(biāo)函數(shù)進(jìn)行參數(shù)優(yōu)化，得到如圖3所示的攻角曲線。

圖3 分段的攻角曲線Fig.3 Segment Curve of Angle of Attack

檢查攻角剖面光滑性，見(jiàn)圖4、圖5。結(jié)果表明，攻角剖面在下降段前期與傳統(tǒng)最大正攻角方案一致，而在過(guò)渡點(diǎn)附近與平衡滑翔段光滑銜接，攻角、攻角一階導(dǎo)數(shù)均連續(xù)。攻角剖面一階導(dǎo)數(shù)大小在過(guò)渡點(diǎn)附近不超過(guò)1 (°)/s，二階導(dǎo)數(shù)大小在過(guò)渡點(diǎn)附近不超過(guò)2 (°)/s2，基本滿足張冉[8]對(duì)攻角光滑性的要求。

圖4 攻角一階導(dǎo)數(shù)Fig.4 First Derivative of Angle of Attack

圖5 攻角二階導(dǎo)數(shù)Fig.5 Second Derivative of Angle of Attack

采用優(yōu)化后的攻角剖面對(duì)CAV-H飛行器的再入下降段和平衡滑翔段進(jìn)行數(shù)值仿真，見(jiàn)圖6、圖7。

圖6 高度-時(shí)間曲線Fig.6 Hight-Time Curve

圖7 彈道傾角及其導(dǎo)數(shù)曲線Fig.7 Curves ofγand its First Derivative

結(jié)果表明，再入飛行器能夠平滑地進(jìn)入擬平衡滑翔狀態(tài)，彈道傾角由初始值過(guò)渡到零附近的過(guò)程連續(xù)且二階可導(dǎo)，一階導(dǎo)數(shù)絕對(duì)值不超過(guò)0.05 (°)/s。

3.2 不同攻角剖面方案的光滑性對(duì)比

采用式（11）中的目標(biāo)函數(shù)f作為控制量曲線的光滑性能指標(biāo)，對(duì)單波峰、分段線性、分段常值3種攻角剖面進(jìn)行光滑性對(duì)比，見(jiàn)圖8和表1。

圖8 攻角曲線對(duì)比Fig.8 Comparision of Diffierent AOA Curves

表1 不同攻角剖面方案的光滑性對(duì)比Tab.1 Comparison of Smoothness with Diffierent AOA Profiles

結(jié)果表明，復(fù)合三角函數(shù)提高了攻角剖面的光滑性，使得曲線在分段過(guò)渡點(diǎn)處保持了較好的連續(xù)性。即使不對(duì)參數(shù)進(jìn)行優(yōu)化，單波峰攻角剖面相比另外兩者有先天性的光滑性優(yōu)勢(shì)。此外，以上 3種剖面參數(shù)化方法的參數(shù)數(shù)量是相同的（均為3個(gè)），因此在軌跡優(yōu)化時(shí)的復(fù)雜度基本一致。

3.3 不同攻角剖面方案的軌跡特性對(duì)比

為了檢驗(yàn)單波峰攻角方案的優(yōu)越性，將仿真結(jié)果與 Harpold[3]分段線性攻角方案、周浩[4]分段常值攻角方案得到的軌跡進(jìn)行對(duì)比。采用胡錦川[7]提出的性能指標(biāo)來(lái)評(píng)價(jià)軌跡平滑性：

式中 aε為縱向加速度，aε=V γ˙；tf為積分時(shí)間，tf=160 s。檢查最大過(guò)載 nmax、峰值熱流 Q˙max、過(guò)渡點(diǎn)動(dòng)壓q等軌跡特性，見(jiàn)表2。

表2 3種攻角剖面的軌跡特性對(duì)比Tab.2 Comparison of Trajectory Characteristics with Three Diffierent AOA Profiles

結(jié)果表明，經(jīng)過(guò)優(yōu)化的單波峰攻角剖面得到的軌跡最平滑，其他軌跡性能指標(biāo)與Harpold[3]攻角方案差別不大，而分段常值攻角剖面的軌跡特性最差。分析可知，分段線性攻角方案在α-t和α-V剖面均連續(xù)但不可導(dǎo)，導(dǎo)致軌跡的平滑性在過(guò)渡點(diǎn)處下降，而分段常值攻角方案難以保證軌跡平滑，無(wú)法滿足再入軌跡精細(xì)化設(shè)計(jì)的任務(wù)需求。

4 結(jié)束語(yǔ)

本文通過(guò)引入新的基函數(shù)，對(duì)原有控制量剖面進(jìn)行了改進(jìn)。新方法在保證參數(shù)數(shù)量不增加的情況下改善了控制量剖面的光滑性，基本滿足胡錦川[7]、張冉[8]對(duì)控制量曲線的光滑性要求。通過(guò)差分進(jìn)化算法對(duì)基函數(shù)進(jìn)行參數(shù)優(yōu)化，使得過(guò)渡點(diǎn)處的攻角曲線光滑性較好，提高了軌跡平滑性，下降段的熱流、過(guò)載等環(huán)境也有一定程度改善。

復(fù)合三角函數(shù)具有在前、后端平緩變化的特點(diǎn)，適合在高速、大動(dòng)壓情況下控制飛行器的姿態(tài)角，可以在其他場(chǎng)景中加以進(jìn)一步的應(yīng)用。

[1] 水尊師, 周軍, 葛致磊. 基于高斯偽譜方法的再入飛行器預(yù)測(cè)校正制導(dǎo)方法研究[J]. 宇航學(xué)報(bào), 2011, 32(6)： 1249-1255.

Shui Zunshi, Zhou Jun, Ge Zhi-lei. On-line predictor-corrector reentry guidance law based on gauss pseudospectral method[J]. Journal of Astronautics, 2011, 32(6)： 1249-1255.

[2] 肖紅, 王勇, 國(guó)海峰, 王征, 臧旭. 基于航路點(diǎn)的高超聲速飛行器預(yù)測(cè)校正制導(dǎo)律設(shè)計(jì)[J]. 電光與控制, 2013, 20(10)： 33-37.

Xiao Hong, Wang Yong, Guo Haifeng, Wang Zheng, Zang Xu. Design of waypoint based predictor-corrector guidance law for hypersonic vehicles[J]. Electronics Optics＆Control, 2013, 20(10)： 33-37.

[3] Harpold J C, Graves C A. Shuttle entry guidance[J]. The Journal of the Astronautical Sciences, 1979, 37(3)： 239-268.

[4] 周浩, 陳萬(wàn)春, 殷興良. 高超聲速飛行器多段彈道優(yōu)化設(shè)計(jì)[J]. 宇航學(xué)報(bào), 2008, 29(4)： 1205-1209.

Zhou Hao, Chen Wanchun, Yin Xingliang. Multi-phase trajectory optimization for a hypersonic vehicle[J]. Journal of Astronautics, 2008,29(4)： 1205-1209.

[5] 李強(qiáng). 高超聲速滑翔飛行器再入制導(dǎo)控制技術(shù)研究[D]. 北京： 北京理工大學(xué), 2015.

Li Qiang. Study on reentry guidance and control method for hypersonic glide vehicle[D]. Beijing： Beijing Institute of Technology, 2015.

[6] 劉運(yùn)鵬, 李伶. 高超聲速導(dǎo)彈高空再入減速段彈道優(yōu)化設(shè)計(jì)[J]. 航天控制, 2010.11, 28(6)： 3-6.

Liu Yunpeng, Li Ling. The optimization design for high-altitude reentry deceleration trajectory of hypersonic missile[J]. Aerospace Control, 2010,28(6)： 3-6.

[7] 胡錦川, 陳萬(wàn)春. 高超聲速飛行器平穩(wěn)滑翔彈道設(shè)計(jì)方法[J]. 北京航空航天大學(xué)學(xué)報(bào), 2015, 41(8)： 1464-1475.

Hu Jinchuan, Chen Wanchun. Steady glide trajectory planning method for hypersonic reentry vehicle[J]. Journal of Beijing University of Aeronautics and Astronautics, 2015, 41(8)： 1464-1475.

[8] 張冉, 李惠峰. 基于平衡滑翔的升力再入攻角設(shè)計(jì)分析及應(yīng)用[J].宇航學(xué)報(bào), 2012, 33(11)： 1557-1563.

Zhang Ran, Li Huifeng. Angle of attack operation zone analysis and application based on equilibrium glide condition for lifting reentry[J].Journal of Astronautics, 2012, 33(11)： 1557-1563.

[9] Vinh N X, Busemann A, et al. Hypersonic and planetary entry flight mechanics[M]. Ann Arbor MI： Univ of Michigan Press, 1980.

[10] Zhang K, Chen W. Reentry vehicle constrained trajectory optimization[R]. AIAA-2011-2231, 2011.

[11] Storn R, Price K. Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces[J]. Journal of Global Optimization, 1997, 11 (4)： 341-359.