-
高超声速滑翔飞行器再入运动方程如下:
$$ \left\{\begin{array}{l}\dot{r}=V\mathrm{sin}\gamma \\ \dot{\theta }=\dfrac{V\mathrm{cos}\gamma \mathrm{sin}\psi }{r\mathrm{cos}\varphi }\\ \dot{\varphi }=\dfrac{V\mathrm{cos}\gamma \mathrm{cos}\psi }{r}\\ \dot{V}=-D/m-g\mathrm{sin}\gamma \\ \dot{\gamma }=\dfrac{L{\rm{cos}}\sigma }{mV}-\dfrac{gr-{V}^{2}}{Vr}{\rm{cos}}\gamma \\ \dot{\psi }=\dfrac{L\mathrm{sin}\sigma }{mV\mathrm{cos}\gamma }+\dfrac{V}{r}\mathrm{cos}\gamma \mathrm{sin}\psi \mathrm{tan}\varphi \end{array} \right.$$ (1) 式中:
$ r $ 表示飞行器到地心相对距离;$ \theta $ 、$ \varphi$ 分别为当前经度和纬度;$ V $ 为飞行器相对地球速度;$ \gamma $ 、$ \psi $ 分别为飞行路径角和航向角;$ L $ 和$ D $ 分别为阻力和升力;$ \sigma $ 为倾侧角;$ m $ 为飞行器质量,$ g $ 为地球引力加速度。阻力和升力计算方式如下:
$$ \left\{ \begin{gathered} L = 0.5\rho {V^2}{C_L}{S_r} \hfill \\ D = 0.5\rho {V^2}{C_D}{S_r} \hfill \\ \end{gathered} \right. $$ (2) 式中:
$\; \rho $ 为大气密度;$ h = r - R $ ,$ R $ 为地球平均半径;$ {S_r} $ 为飞行器参考面积;$ {C_L} $ 和$ {C_D} $ 分别为飞行器升力系数和阻力系数。 -
高超声速飞行器再入机动飞行约束条件主要包括终端约束、控制量约束和过程约束。典型的过程约束包含热流率密度约束、动压约束和过载约束。
控制量包括攻角
$ \alpha $ 和倾侧角$ \sigma $ ,控制量约束记为:$$ \left\{ \begin{gathered} {\text{|}}\alpha {\text{|}} \leqslant {\alpha _{{\text{max}}}} \hfill \\ {\text{|}}\sigma {\text{|}} \leqslant {\sigma _{{\text{max}}}} \hfill \\ \end{gathered} \right. $$ (3) 式中:下标max表示控制量的最大值。
高超声速再入轨迹优化目的是控制飞行器到达指定区域或者指定目标点,文中研究过程中终端约束记为:
$$ \left\{ \begin{gathered} |r({t_f}) - {r_f}| \leqslant \Delta r \hfill \\ |S({t_f}) - {S_f}| \leqslant \Delta S \hfill \\ |V({t_f}) - {V_f}| \leqslant \Delta V \hfill \\ \end{gathered} \right. $$ (4) 式中:
$ {t_f} $ 表示终端时刻,下标“$ f $ ”表示终端时刻期望的状态值;$ \Delta r $ 、$ \Delta S $ 和$ \Delta V $ 分别指终端时刻高度、航程和速度的可接受误差。过程约束可记为:
$$ \left\{ \begin{gathered} {N_l} = \sqrt {{L^2} + {D^2}} /m \leqslant {N_{{\text{max}}}} \hfill \\ q = 0.5\rho {V^2} \leqslant {q_{{\text{max}}}} \hfill \\ \dot Q = {K_Q}{\rho ^{0.5}}{V^{3.15}} \leqslant {{\dot Q}_{{\text{max}}}} \hfill \\ \end{gathered} \right. $$ (5) 式中:
$ {N_l} $ 、$ q $ 、$ \dot Q $ 分别表示过载、动压、热流率密度;下标max表示对应约束的最大限值;${K_Q} = 7.968\;6 \times $ $ {10^{ - 5}}\;{\text{J}}{{\text{s}}^{\text{2}}}{\text{/(}}{{\text{m}}^{{\text{3}}{\text{.5}}}}\;{\text{k}}{{\text{g}}^{{\text{0}}{\text{.5}}}}{\text{)}}$ 为热流率常数。此外,高超声速滑翔飞行器在滑翔段飞行时,还需要满足准平衡滑翔条件(Quasi-Equilibrium Glide Condition, QEGC),具体形式如下[19]:
$$ \frac{L}{m}\mathrm{cos}\sigma +\left(\frac{{V}^{2}}{r}-g\right)=0 $$ (6) -
高超声速飞行器飞行速度快的特点,虽然有利于对目标的快速打击,但也会导致气动加热现象严重,从而使的飞行器在高超声速飞行过程中会产生较强的红外辐射。文中拟从轨迹优化的角度,设计一条全过程中红外辐射特性最小的轨迹。
高超声速再入滑翔飞行器滑翔飞行过程中,其红外辐射源主要来自于自身蒙皮气动加热和对太阳光的反射辐射,但由于高速飞行产生的蒙皮气动热辐射远远超过对太阳的反射辐射,因此,分析高超声速滑翔飞行器时的红外特性主要考虑蒙皮加热产生的辐射。理想状态下,蒙皮的气动热带来的红外辐射可由其驻点温度来表征,驻点温度可由下式计算[20]:
$$ {T_s}(h) = {T_0}(h)\left[ {1 + {a_0}{{(\rho /{\rho _1})}^{0.5}}M{a^2}(h)} \right] $$ (7) 式中:
$ {T_s}(h) $ 为当前飞行器表面温度;$ {T_0}(h) $ 为高度$ h $ 处的大气温度;$ {a_0} = 0.164 $ ;$ \;{\rho _1} $ 为海拔$ 19 \;{\text{km}} $ 时对应的大气密度;$ Ma(h) $ 为高度$ h $ 处飞行器飞行马赫数,其计算公式为:$$ Ma(h) = \frac{V}{{20.046\;8 \cdot {{({T_0})}^{0.5}}}} $$ (8) 将飞行器看成一个黑体,由于普朗克公式能够适用于整个电磁波段的黑体辐射计算问题,因此,文中采用普朗克公式计算高超声速再入飞行器驻点的红外辐射强度,可近似表示为:
$$ {W_b} = \frac{{{c_1}}}{{{\lambda ^5}}} \cdot \frac{1}{{{e^{{c_2}/\lambda {T_s}}} - 1}} $$ (9) 式中:
$ \lambda $ 为波长;${c_1} = 3.741\;8 \times {10^8}\;{\rm{W}} \cdot {{\rm{m}}^{ - 2}} \cdot \text{μ}{\rm{m}} {^4}$ ;${c_2} = $ $ 14\;388\; \text{μ}{\rm{m}} \cdot {\rm{K}}$ 。为了减少高超声速滑翔飞行器红外探测暴露的风险,从轨迹优化角度,设计文中目标函数为全程总红外辐射最小,记为:
$$ \min\quad J(\alpha ,\sigma ) = \int_{{t_0}}^{{t_f}} {{W_b}{\rm{d}}t} $$ (10) 综上,高超声速再入轨迹优化问题是寻找最优的控制量变化剖面,以使得公式(10)为目标函数最小,同时满足公式(1)~(6)的约束。
-
WOA是由Mirjalili等在2016年提出的,参考文献[17]指出,WOA相较于其他常见优化算法在收敛性和全局性方面具有一定的优势。但与其他优化算法类似,面对复杂问题时,WOA也会出现计算效率不高和可能陷入局部最优的问题。文中针对该类问题,对标准WOA进行改进,以增强面对复杂问题时的全局搜索能力。
-
WOA是通过模拟自然界鲸鱼群的捕食方式构造出的一种群智能优化方法。捕食过程可简单分为三步:包围猎物、气泡网攻击和寻找猎物。
三个环节分别按照各自的规则进行更新[16]。
包围猎物环节中,其他鲸鱼根据最佳的鲸鱼位置进行更新,进一步包围猎物,更新规则如下:
$$ M(t + 1) = {M^*}(t) - b \cdot D $$ (11) $$ D = \left| {c \cdot {M^*}(t) - M(t)} \right| $$ (12) 式中:
$ M(t) $ 表示第$ t $ 次迭代时鲸鱼的位置;$ {M^*}(t) $ 是当前全局最佳位置;$ b $ ,$ c $ 为步长系数,能够控制向最佳位置移动的步长,分别表示为:$$ b = 2a \cdot {\xi _1} - a $$ (13) $$ c = 2{\xi _1} $$ (14) 式中:
$ {\xi _1} $ 为$ [0,1] $ 中的随机向量;$ a = 2(1 - t/T) $ ,$ T $ 为最大迭代次数。气泡网攻击环节中,鲸鱼吐出气泡以对数螺旋线的方式靠近猎物,更新规则如下:
$$ M(t + 1) = D \cdot {{\rm{e}}^{d{\xi _2}}} \cdot \cos (2\pi {\xi _2}) + {M^*}(t) $$ (15) 式中:
$ d $ 为常数;$ {\xi _2} $ 为$ [0,1] $ 中的随机向量。搜寻猎物阶段,鲸鱼随机选择其他鲸鱼位置进行觅食,进行全局搜索,更新规则如下:
$$ M(t + 1) = {M_r}(t) - b \cdot D $$ (16) $$ D = \left| {c \cdot {M_r}(t) - M(t)} \right| $$ (17) 式中:
$ {M_r}(t) $ 表示当前迭代次数下的鲸鱼种群中随机选择某一个体位置作为更新靠近的目标位置。 -
针对WOA存在的易陷入局部最优的问题,文中提出了一种改进的WOA,在保证收敛速度的同时,增强WOA的全局搜索能力。在参考文献[16-18]工作的启发下,文中的改进思路主要分三个方面:(1) 种群初始化阶段,通过Tent混沌序列和反向学习策略(Opposition-based Learning, OBL)进行种群初始化;(2) 包围猎物环节通过控制参数
$ a $ 的余弦变化,控制搜索步长;(3) 搜寻猎物阶段,选取前两个历史最优解,使得鲸鱼位置更新快速指向潜在最优解。首先,种群初始化阶段,生成Tent混沌序列计算方式如下:
$$ {p_{k + 1}} = \left\{ \begin{gathered} 2{p_k},\quad \quad \quad {p_k} \lt 0.5 \hfill \\ 2(1 - {p_k}),\quad {p_k} \geqslant 0.5 \hfill \\ \end{gathered} \right. $$ (18) 式中:“
$ k $ ”指的是变量的第$ k $ 维,$ k = 1 $ 时;$ {p_1} $ 是$ [0,1] $ 内的随机量。将混沌序列映射到解空间得到混沌种群,计算如下:$$ M_k^p = {B_{lk}} + {p_k} \cdot ({B_{uk}} - {B_{lk}}) $$ (19) 式中:上标“
$ p $ ”指由混沌序列$ p $ 映射得到的值;$ M_k^p $ 指鲸鱼个体第$ k $ 维位置上的取值;$ {B_{uk}} $ 和$ {B_l}_k $ 分别指第$ k $ 维上解空间的上限和下限。将混沌序列映射得到的种群记为$ {U^p} $ ,则可利用反向策略生成种群$ {U^o} $ ,$ {U^o} $ 种个体位置计算如下:$$ M_k^o = {B_l}_k + ({B_u}_k - M_k^p) $$ (20) 式中:
$ M_k^o $ 表示反向种群第$ k $ 个维位置上的取值。设初始化种群大小为$ N $ ,将$ {U^p} $ 和$ {U^o} $ 种群合并,计算适应度并排序。应用精英策略取出适应度最好的$ N $ 个鲸鱼位置组成算法的初始化种群,记为$ {U_0} $ 。通过分析公式(11)、(13)与(16), 可知,控制参数
$ a $ 的变化能够影响算法全局搜索与局部搜索的能力,在种群迭代初期,需要$ a $ 尽量大保持较好的全局搜索能力,在迭代后期,需要$ a $ 尽快减小进行局部搜索,使得算法快速收敛至最优解。原算法中$ a $ 线性减小,对于高维解空间搜索时,易陷入局部最优,因此参考文献[18],对控制参数$ a $ 进行余弦变化的改进,计算形式如下:$$ a = 2 \cdot \cos \left(\frac{\pi }{2} \cdot \frac{t}{T}\right) $$ (21) 分析可知,在迭代前期,
$ a $ 缓慢减小全局搜索充分,算法后期能够快速收敛,专注局部搜索。搜寻猎物阶段,WOA随机选取其他鲸鱼位置最为更新目标,但是这种随机搜索方式在算法后期不利于快速寻优,因此考虑联系当前状态最优解与其他次优解,按照适应度大小排序,选取适应度最优的前两个鲸鱼位置,用于搜寻猎物阶段的位置更新,由公式(16)、(17)得到鲸鱼位置更新:
$$ M(t + 1) = {M_c}(t) - b \cdot \left| {c \cdot {M_c}(t) - M(t)} \right| $$ (22) $$ {M_c}(t) = \frac{{f({M^F}(t))}}{{G(t)}} \cdot {M^F}(t) + \frac{{f({M^S}(t))}}{{G(t)}} \cdot {M^S}(t) $$ (23) $$ G(t) = f({M^F}(t)) + f({M^S}(t)) $$ (24) 式中:
$ {M^F}(t) $ 和$ {M^S}(t) $ 分别为鲸鱼最优位置及次优位置;$ f $ 指在当前位置上的适应度。图1所示为IWOA流程图。 -
文中从轨迹优化的角度,提出一种使得全程总红外辐射最小的轨迹优化算法,进而达到削弱气动热辐射效应对红外探测窗口性能的影响。高超声速飞行器再入轨迹优化问题目的是在找到一条满足各种飞行约束且使得目标函数最优的飞行轨迹,整个问题解决范式是求解攻角剖面和倾侧角剖面。文中首先将其攻角剖面和倾侧角剖面求解转换成一个参数优化问题,进而采用IWOA进行优化求解。
-
文中进行研究时将攻角剖面设计为速度相关的线性分段函数,具体形式如下:
$$ \alpha = \left\{ \begin{gathered} {\alpha _{\max }},\quad \quad {V_1} < V \leqslant {V_0} \hfill \\ \frac{{({\alpha _r} - {\alpha _{\max }})}}{{({V_2} - {V_1})}} \cdot (V - {V_1}) + {\alpha _{\max }},{V_2} < V \leqslant {V_1} \hfill \\ {\alpha _r},\qquad {V_f} < V \leqslant {V_2} \hfill \\ \end{gathered} \right. $$ (25) 式中:
$ {V_0} $ 为再入时刻的速度;$ {\alpha _{\max }} $ 为攻角的最大值;$ {\alpha _r} $ 为最大升阻比对应的攻角值;$ {V_1} $ 和$ {V_2} $ 为算法优化设计攻角剖面相关的参数。再入飞行初期,高超声速滑翔飞行器高度较高,大气稀薄,受到气动力较小,为保证稳定飞行,需尽快进入滑翔飞行段,因此初始阶段,攻角要尽可能大,故选择攻角最大值作为初始值。高升阻比的特性给高超声速飞行器带来了较好的机动能力,故最终阶段选用最大升阻比对应的攻角值。 -
高超声速滑翔飞行器再入飞行过程可分为初始再入段和准平衡滑翔段。准平衡滑翔段飞行以热流率、动压、过载等路径约束和准平衡滑翔条件为边界构成再入飞行走廊,文中研究沿用作者前期工作,采用阻力加速度-速度再入走廊,阻力加速度记为
$ {a_D} $ ,走廊边界具体形式如下:$$ \left\{ \begin{gathered} {a_D} \leqslant \min \left\{ {\frac{{{{\dot Q}_{\max }}{S_r}{C_D}{\rho ^{0.5}}}}{{2{K_Q}m{V^{1.15}}}},\frac{{{q_{\max }}{S_r}{C_D}}}{m},\left. {\frac{{{N_{\max }}{g_0}{D^2}}}{{\sqrt {({L^2} + {D^2})} }}} \right\}} \right. \hfill \\ {a_D} \geqslant \frac{D}{L} \cdot (g - \frac{{{V^2}}}{r}) \hfill \\ \end{gathered} \right. $$ (26) 为快速求解倾侧角剖面,文中受参考文献[14]启发,设计出一种倾侧角一次翻转剖面,具体形式如下:
$$ \sigma = \left\{ \begin{gathered} {\sigma _p},t \leqslant {t_{QEGC}} \hfill \\ {\sigma _p} \cdot \frac{\varsigma }{{{\varsigma _{QEGC}}}},t > {t_{QEGC}} \hfill \\ \end{gathered} \right. $$ (27) 式中:
$ {\sigma _p} $ 为初始再入段飞行时的倾侧角;$ \varsigma $ 为当前时刻航向角偏差;$ {t_{QEGC}} $ 和$ {\varsigma _{QEGC}} $ 分别为由初始再入段到准平衡滑翔段转换点时刻和航向角偏差。初始再入段过程大气稀薄,气动控制力作用有限,因此,此段飞行过程倾侧角大小可设为是一个常值
$ {\sigma _0} $ ,直至到达准平衡滑翔段的转换点。则初始再入段飞行的倾侧角$ {\sigma _p} $ 计算形式如下:$$ {\sigma _p} = {\rm{sign}}({\sigma _0}) \cdot {\sigma _0} $$ (28) 式中:
${\rm{sign}}({\sigma _0})$ 为再入初始时刻的倾侧角符号,可由下式确定:$$ {\rm{sign}}({\sigma _0}) = - {\rm{sign}}({\varsigma _0}) $$ (29) 式中:
$ {\varsigma _0} $ 为初始再入时刻,航向角与到目标点视线方向的偏差,$ {\varsigma _0} = {\psi _0} - {\psi _{los}} $ ,$ {\psi _0} $ 为再入时刻航向角,$ {\psi _{los}} $ 为目标点视线方位角,计算形式如下:$$ \tan {\psi _{los}} = \frac{{\sin ({\theta _f} - \theta )}}{{\cos \varphi {\text{tan}}{\varphi _f} - \sin \varphi \cos ({\theta _f} - \theta )}} $$ (30) 式中:
$ {\theta _f} $ 和$ {\varphi _f} $ 分别为目标点的经度、纬度。因此,初始再入段倾侧角取值结合攻角剖面和初始时刻状态量,不断积分计算集合得到初始再入段轨迹,直至进入阻力加速度-速度再入走廊,到达准平衡滑翔段飞行,即可得到此时转换点时刻$ {t_{QEGC}} $ 和航向角偏差$ {\varsigma _{QEGC}} $ 。初始再入段到准平衡滑翔段的转换点,转换条件为:$$ \left| {\frac{{{\rm{d}}r}}{{{\rm{d}}V}} - {{\left(\frac{{{\rm{d}}r}}{{{\rm{d}}V}}\right)}_{{{QEGC}}}}} \right| < {\delta _1} $$ (31) 式中:
$ {\delta _1} $ 为小量常值。由公式(1)和(6)可得:$$ \left\{ \begin{gathered} \frac{{{\rm{d}}r}}{{{\rm{d}}V}} = \frac{{V\sin \gamma }}{{ - D/m - g\sin \gamma }} \hfill \\ {\left( {\frac{{{\rm{d}}r}}{{{\rm{d}}V}}} \right)_{{\text{QEGC}}}} \doteq \frac{{{C_L}{S_r}\cos {\sigma _0}\rho V + 2mV/r}}{{m{V^2}/{r^2} + \rho {V^2}{C_L}{S_r}\cos {\sigma _c}/2{h_s}}} \hfill \\ \end{gathered} \right. $$ (32) 综上,在倾侧角剖面中设计了一个参数
$ {\sigma _0} $ ,通过不断的积分迭代即可求得每一时刻的状态量和控制量。 -
文中采用罚函数方法处理终端约束,根据终端约束允许误差进行罚函数的设置,则目标函数公式(7)可转换为:
$$ \min \widetilde J(\alpha ,\sigma ) = \int_{{t_0}}^{{t_f}} {{W_b}dt} + {w_1}{b_r} + {w_2}{b_S} + {w_3}{b_V} $$ (33) 式中:
$ {w}_{1} $ ,$ {w_2} $ ,$ {w_3} $ 为相关约束的惩罚因子;$ {b_r} $ ,$ {b_S} $ 和$ {b_V} $ 由公式(4)转化可得:$$ \begin{gathered} {b_r}{\text{ = }}\left\{ \begin{gathered} 0,\quad |r({t_f}) - {r_f}| \leqslant \Delta r \hfill \\ \frac{{|r({t_f}) - {r_f}|}}{{1\;000}},|r({t_f}) - {r_f}| > \Delta r\quad \hfill \\ \end{gathered} \right. \hfill \\ {b_S}{\text{ = }}\left\{ \begin{gathered} 0,\quad |S({t_f}) - {S_f}| \leqslant \Delta S \hfill \\ \frac{{|S({t_f}) - {S_f}|}}{{1\;000}},|S({t_f}) - {S_f}| > \Delta S\quad \hfill \\ \end{gathered} \right. \hfill \\ {b_V}{\text{ = }}\left\{ \begin{gathered} 0,\quad |V({t_f}) - {V_f}| \leqslant \Delta V \hfill \\ \frac{{|V({t_f}) - {V_f}|}}{{1\;000}},|V({t_f}) - {V_f}| > \Delta S\quad \hfill \\ \end{gathered} \right. \hfill \\ \end{gathered} $$ (34) 由此可以看出,文中在考虑终端约束时设置了一定的容许误差,即对终端约束进行一定的松弛,当终端误差小于可接受误差时,约束误差项记为零;当大于可接受误差时,则将其变成一个惩罚项加入目标函数。
综上,高超声速再入轨迹优化问题转化成一各参数优化问题,可以采用文中出的IWOA进行寻优求解。
-
文中求解高超声速再入轨迹优化问题时,共计设计
$ {V_1} $ 、$ {V_2} $ 、$ {\sigma _0} $ 三个参数,采用IWOA对进行优化,得到满足任务目标的攻角剖面和倾侧角剖面。文中进行三个仿真算例验证文中提出算法的有效性:(1)使用IWOA进行再入轨迹优化任务,验证文中算法解决再入轨迹优化问题有效性。(2)使用IWOA进行再入轨迹优化任务,并与原始SSA、PSO和WOA进行对比。(3)进行300次蒙特卡洛打靶仿真实验,验证文中算法的鲁棒性。
-
文中以全程总红外辐射最小为目标函数进行再入轨迹优化仿真实验,采用CAV-H模型进行算法验证,飞行器气动系数模型采用考虑攻角及速度幂指数多项式形式[21],飞行器参考质量907.2 kg,参考横截面积为0.484
$ {{\text{m}}^2} $ ,终端约束允许误差量分别为:$ \Delta r = $ $ 1\;000\;{\text{m}} $ ,$ \Delta S = 10{\text{km}} $ ,$ \Delta V=100\;\text{m/s} $ 。本节实验设置约束条件和设计的三个参数约束范围如表1所示。IWOA种群数量设为30,最大迭代次数为20。表 1 仿真场景参数
Table 1. Simulation scene parameters
State Initial condition Terminal constraint Other parameters Value $ r $/km 60 25 $ {\dot{Q}}_{\text{max}} $/$ {\text{kW/}}{{\text{m}}^{-{\text{2}}}} $ 1200 $ \theta $/(°) 120 160 $ {q}_{\text{max}} $/${\text{kPa} }$ 400 $ \varphi $/(°) 40 65 $ {n}_{\text{max}} $ 6 $ V $/m·s−1 5000 1200 $ {V}_{1} $/m·s−1 $ [4\;000,4\;800] $ $ \gamma $/(°) −1 - $ {V}_{2} $/m·s−1 $ [1\;800,4\;000] $ $ \psi $/(°) 40 - $ |{\sigma _0}| $/(°) $ [0,80] $ 图2~图9分别展示了文中算法再入轨迹优化的结果。图2~图4分别展示了高度、经纬度和速度的变化情况,可以看出文中算法能够制导飞行器精确到达目标点,高度、速度、航程误差分别为1.13 m、1.24 m/s和73.26 m。均满足允许终端高度误差。图5所示为所得轨迹能够顺利进入阻力加速度再入走廊,表明了考虑了飞行器结构性能的前提下,飞行轨迹能够满足路径约束,确保飞行安全。图6~图7分别展示了实验场景下控制量的变化,攻角速度剖面的变化符合所设计的分段线性函数的变化;倾侧角翻转次数较少,且只有一次翻转到之后迅速趋于稳定值,证明文中所设计的倾侧角一次翻转策略的有效性,对控制系统及执行机构具有参考意义。
图8所示为再入飞行过程中飞行器表面驻点温度的变化,飞行过程中最大温度为1921.39 K。图9所示为飞行过程中三种典型红外波段辐射曲线,分析可知,飞行器跳跃滑翔过程中,弹道波谷时,温度急剧升高,辐射强度迅速增强,即辐射强度与气动热产生的温度成正比,辐射强度波峰与弹道波谷重合。初始再入时,飞行速度大,短波的辐射强度远高于中长波,飞行后期时,短波的辐射强度迅速下降,中长波辐射更为剧烈。
表2所示为该算法得到的参数优化结果,切换速度和倾侧角取值均满足算法设置的上下限范围;以短波为例,再入飞行全过程的总红外辐射强度为 29560.64
${\text{kW}}/{\text{sr}}$ 。表 2 四种算法优化结果
Table 2. Optimization results of four algorithms
Variable $ {V_1}$/m·s−1 $ {V}_{2} $/m·s−1 $|{\sigma _0}|$/(°) $ J(u)$/kW·sr−1 Value 4000.01 1800.21 41.48 29560.64 综上,文中提出的基于IWOA高超声速飞行器再入轨迹优化方法能够得到一条精确且安全的满足飞行任务的轨迹,经分析知,飞行过程中,红外辐射强度峰值与弹道波谷重合,辐射强度与气动热产生的温度成正比,最大温度、最大辐射强度及总辐射强度对高超声速再入飞行器的工业设计具有参考意义。
-
以4.1节中实验参数设置为例,分别进行基于IWOA、WOA、SSA、PSO轨迹优化仿真实验,采用再入走廊方法处理路程约束,罚函数方法处理三个终端约束,将四种算法得到结果进行对比。
表3所示为四种优化算法进行仿真得到的设计参数值和目标函数值。可以发现,四种算法得到设计参数值都能够满足上下限约束,文中算法得到的目标函数值为最优,说明了文中算法的全局寻优能力更强。
表 3 四种算法优化结果
Table 3. Simulation results of four algorithms
Algorithms $ {V_1} $/m·s−1 $ {V_2} $/m·s−1 $|{\sigma _0}|$/(°) $J(u)$/kW·sr−1 IWOA 4000.01 1800.21 41.48 29560.64 WOA 4063.99 1951.63 45.47 34563.88 SSA 4000.52 1800.89 43.54 31268.92 PSO 4135.48 1895.52 47.94 38287.97 图10~图12分别展示了四种算法得到高度、经纬度和速度的变化,四种优化算法的高度、速度、航程均能满足终端约束,但由图可知,IWOA算法的到的高度、速度和经纬度更接近设定的目标值,说明文中算法得到的轨迹相较于对比算法更加精确。图13展现四种算法所得轨迹均能进入再入走廊,说明四种方法得到的轨迹均能够满足路径约束,保证了再入飞行过程的安全顺利。图14~图15分别展现了攻角和倾侧角随时间的变化,控制量变化趋势均符合文中提出的轨迹优化算法框架的设计要求。
图16所示为四种算法得到的飞行器驻点温度变化情况,图17以短波为例,展现了四种算法得到再入飞行辐射曲线,可知IWOA得到的飞行轨迹温度峰值最小,对应的全程总辐射强度也最小,相较于PSO、WOA、SSA分别减小23.4%、14.5%和5.5%,符合表3中再入飞行全过程IWOA算法轨迹总红外辐射最小的结果。表4所示为四种算法的终端误差、温度峰值和辐射强度峰值的比较,可以看出文中算法结果相较于对比算法终端误差更小,温度峰值和辐射峰值更低。
表 4 优化结果比较
Table 4. Comparison of optimization results
Terminal errors IWOA WOA SSA PSO Altitude/km 0.077 0.84 1.12 0.74 Velocity/m·s−1 0.0086 0.24 97.56 10.50 Longitude/(°) 0.30 0.72 0.86 1.70 Latitude/(°) 0.07 0.19 0.21 0.42 Peak temperature/K 1921.37 2061.98 1977.75 2173.41 Peak radiation/${\text{kW} } \cdot {\text{s} }{ {\text{r} }^{ - 1} }$ 283.32 368.04 316.12 442.79 综上,图10~图17展现了文中设计的轨迹优化算法框架,即参数化设计分段攻角剖面和倾侧角一次翻转剖面,对常见的群智能优化算法具有一定的适应性,且文中给出的IWOA具有更强的全局寻优能力,得到的再入轨迹更加精确。
-
考虑在初始条件存在扰动情况下进行再入轨迹优化实验仿真,验证文中提出轨迹优化算法的鲁棒性。参考文献[22]设置初始条件扰动数值如表5所示。实验约束设置同4.1节中实验设置,进行300次蒙特卡洛打靶实验。
表 5 蒙特卡洛仿真扰动因素设置
Table 5. Disturbances in the Monte Carlo
Disturbance $ 3\sigma $values Altitude/m 200 Longitude/(°) 1 latitude/(°) 1 Velocity/m·s−1 100 Flight path angle/(°) 0.2 Head angle/(°) 0.2 图18~图23所示为300次的蒙特卡洛仿真实验得到的结果,由图18~图19可知,在初始六个状态量存在一定扰动的情况下,300次打靶实验得到的高度、速度能够满足终端约束误差精度,图20表明文中算法能够克服初始状态扰动带来的影响,准确到达目标点,结合表6中所示的终端速度、高度、经度、纬度误差分析,可表明文中算法能够克服初始条件扰动圆满完成打击任务。图21所示为倾侧角随速度变化情况,可知300次拉偏实验倾侧角变化能够保持一致,全部只进行了一次翻转,验证了文中提出的倾侧角一次翻转策略的有效性和稳定性。图22和图23分别展示了打靶实验得到的驻点温度和红外辐射的变化情况,结合表6所示的温度及辐射数据分析,300次蒙特卡洛实验中,温度峰值最大为2213.6 K,以短波为例的辐射峰值最大为471.73
${\text{kW}}/{\text{sr}}$ ,平均总辐射为31592.07${\text{kW}}/{\text{sr}}$ ,这些数据可为高超声速再入飞行器热防护设计、红外探测窗口结构设计甚至反高超声速飞行器探测系统的设计提供一定的参考价值。表 6 蒙特卡洛仿真结果
Table 6. Simulation results of the Monte Carlo
Terminal error Maximum value Average value Standard deviation Altitude/km 0.99 0.71 0.21 Velocity/m·s−1 99.99 94.88 9.68 Longitude/(°) 3.71 1.12 0.42 Latitude/(°) 1.00 0.23 0.20 Peak temerature/K 2213.60 1962.54 178.02 Peak radiation/${\text{kW} } \cdot {\text{s} }{ {\text{r} }^{ - 1} }$ 471.73 314.34 52.20 Total radiation/${\text{kW} } \cdot {\text{s} }{ {\text{r} }^{ - 1} }$ 40278.91 31592.07 3516.83 综上,文中提出的基于IWOA的高超声速飞行器全程红总外辐射最小再入轨迹优化算法,在初始状态扰动条件下能够准确完成再入任务,且具有一定的鲁棒性。
Trajectory optimization of hypersonic glide vehicle with minimum total infrared radiation (Invited)
-
摘要: 为降低高超声速飞行器再入过程中,产生的气动热辐射对红外探测窗口性能的影响,从轨迹优化的角度,以再入飞行全程驻点总红外辐射为目标函数,提出了一种基于改进鲸鱼优化算法(Whale Optimization Algorithm,WOA)的高超声速飞行器轨迹优化算法。首先,通过Tent混沌映射和控制因子余弦变化改进WOA,改进算法位置更新时的位置指向性,增强算法全局搜索能力;同时,将再入轨迹优化问题转化为控制量剖面参数优化问题,采用倾侧角一次翻转策略,利用普朗克公式计算驻点红外辐射,并设计目标函数,利用阻力加速度再入走廊处理路径约束,采用罚函数法将终端约束同目标函数相结合;最后,利用改进的WOA对设计的控制量剖面进行参数寻优,获得使目标函数最优的解。仿真实验表明: 文中改进的WOA能够有效完成全程总红外辐射最小的再入轨迹优化任务,全局搜索能力强,且具有较好的鲁棒性。Abstract: In order to reduce the influence of aerodynamic thermal radiation generated during the reentry process of hypersonic vehicle on the performance of infrared detection window, from the perspective of trajectory optimization, a trajectory optimization algorithm of hypersonic vehicle based on improved whale optimization algorithm was proposed with the total infrared radiation of reentry flight stagnation point as the objective function. Firstly, the Whale optimization algorithm was improved by Tent chaotic map and control factor chord change. The position directivity of the algorithm was improved when the position was updated, and the global search ability of the algorithm was enhanced. At the same time, the reentry trajectory optimization problem was transformed into the parameter optimization problem of the control profile, and a one-time reversal strategy of the inclination angle was proposed. The Planck formula was used to calculate the infrared radiation of the stagnation point, and the objective function was designed. The resistance acceleration was used to enter the corridor to deal with the path constraint. The penalty function method was used to combine the terminal constraint with the objective function. Finally, the improved Whale optimization algorithm was used to optimize the parameters of the designed control profile to obtain the optimal solution of the objective function. The simulation results showed that the improved whale algorithm can effectively complete the reentry trajectory optimization task with the minimum total infrared radiation, and has strong global search ability and good robustness.
-
表 1 仿真场景参数
Table 1. Simulation scene parameters
State Initial condition Terminal constraint Other parameters Value $ r $ /km60 25 $ {\dot{Q}}_{\text{max}} $ /$ {\text{kW/}}{{\text{m}}^{-{\text{2}}}} $ 1200 $ \theta $ /(°)120 160 $ {q}_{\text{max}} $ /${\text{kPa} }$ 400 $ \varphi $ /(°)40 65 $ {n}_{\text{max}} $ 6 $ V $ /m·s−15000 1200 $ {V}_{1} $ /m·s−1$ [4\;000,4\;800] $ $ \gamma $ /(°)−1 - $ {V}_{2} $ /m·s−1$ [1\;800,4\;000] $ $ \psi $ /(°)40 - $ |{\sigma _0}| $ /(°)$ [0,80] $ 表 2 四种算法优化结果
Table 2. Optimization results of four algorithms
Variable $ {V_1}$ /m·s−1$ {V}_{2} $ /m·s−1$|{\sigma _0}|$ /(°)$ J(u)$ /kW·sr−1Value 4000.01 1800.21 41.48 29560.64 表 3 四种算法优化结果
Table 3. Simulation results of four algorithms
Algorithms $ {V_1} $ /m·s−1$ {V_2} $ /m·s−1$|{\sigma _0}|$ /(°)$J(u)$ /kW·sr−1IWOA 4000.01 1800.21 41.48 29560.64 WOA 4063.99 1951.63 45.47 34563.88 SSA 4000.52 1800.89 43.54 31268.92 PSO 4135.48 1895.52 47.94 38287.97 表 4 优化结果比较
Table 4. Comparison of optimization results
Terminal errors IWOA WOA SSA PSO Altitude/km 0.077 0.84 1.12 0.74 Velocity/m·s−1 0.0086 0.24 97.56 10.50 Longitude/(°) 0.30 0.72 0.86 1.70 Latitude/(°) 0.07 0.19 0.21 0.42 Peak temperature/K 1921.37 2061.98 1977.75 2173.41 Peak radiation/ ${\text{kW} } \cdot {\text{s} }{ {\text{r} }^{ - 1} }$ 283.32 368.04 316.12 442.79 表 5 蒙特卡洛仿真扰动因素设置
Table 5. Disturbances in the Monte Carlo
Disturbance $ 3\sigma $ valuesAltitude/m 200 Longitude/(°) 1 latitude/(°) 1 Velocity/m·s−1 100 Flight path angle/(°) 0.2 Head angle/(°) 0.2 表 6 蒙特卡洛仿真结果
Table 6. Simulation results of the Monte Carlo
Terminal error Maximum value Average value Standard deviation Altitude/km 0.99 0.71 0.21 Velocity/m·s−1 99.99 94.88 9.68 Longitude/(°) 3.71 1.12 0.42 Latitude/(°) 1.00 0.23 0.20 Peak temerature/K 2213.60 1962.54 178.02 Peak radiation/ ${\text{kW} } \cdot {\text{s} }{ {\text{r} }^{ - 1} }$ 471.73 314.34 52.20 Total radiation/ ${\text{kW} } \cdot {\text{s} }{ {\text{r} }^{ - 1} }$ 40278.91 31592.07 3516.83 -
[1] Tian B L, Li Z Y, Wu S Y, et al. Reentry trajectory optimization, guidance and control methods for resuable launch vehicles review [J]. Acta Aeronautica et Astronautica Sinica, 2020, 41(11): 624972. (in Chinese) [2] Zhang Y L, Xie Y. Review of trajectory planning and guidance methods for gliding vehicles [J]. Acta Aeronautica et Astronautica Sinica, 2020, 41(1): 45. (in Chinese) [3] Wang Y H, Wang Q, Gao L, et al. Aero-thermo-radiation of a hypersonic vehicle [J]. Infrared and Laser Engineering, 2013, 42(6): 1399-1403. (in Chinese) [4] Li Y H, Hu H Y, Wang Q. Radiative transmission property of infrared window in hypersonic vehicle [J]. Infrared and Laser Engineering, 2020, 49(4): 0404002. (in Chinese) doi: 10.3788/IRLA202049.0404002 [5] Wang Y H, Wang Q, Zhang B C, et al. Experiment of the thermo-radiation characteristic of infrared window of hypersonic vehicles [J]. Infrared and Laser Engineering, 2015, 44(6): 1716-1720. (in Chinese) [6] Zhang K, Chen Z G, Zhao Y Y. Dome protecting technologies for overseas high-velocity guided missiles [J]. Infrared and Laser Engineering, 2013, 42(1): 154-158. (in Chinese) [7] Tan H P, Xia X L, Liu L H. Numberical Calculation on Infrared Radiative Properties and Transfer Calculating Thermal Radiation [M]. Harbin: Harbin Institute of Technology Press, 2006. (in Chinese) [8] Li J, Jiang Z Y. Online trajectory planning algorithm for hypersonic glide re-entry problem [J]. Journal of Beijing University of Aeronautics and Astronautics, 2020, 46(3): 579. (in Chinese) [9] Zong Q, Li Z Y, Ye L Q, et al. Variable trust region se-quential convex programming for RLV online reentry trajectory reconstruction [J]. Journal od Harbin Institute of Technology, 2020, 52(3): 147. (in Chinese) [10] Gao Y, Cai G B, Zhang S X, et al. Reentry maneuver guidance for hypersonic glide vehicles under multiple no-fly zones [J]. Journal of Ordnance Equipment Engineering, 2019, 40(8): 32. (in Chinese) [11] Gong Z F, Liu G, Song R, et al. Reentry tracking control of hypersonic vehicle with fly zone constraints [J]. Journal of National University of Defense Technology, 2020, 42(1): 66. (in Chinese) [12] Li X Q, Ma R, Zhang S, et al. Improved design of ant colony algorithm and its application in path plan-ning [J]. Acta Aeronautica et Astronautica Sinica, 2020, 41(S2): 724381. (in Chinese) [13] Xu H, Cai G B, Zhang S X, et al. Hypersonic reentry optimization by using improved sparrow search algorithm and control parametrization method [J]. Advance in Space Research, 2022, 69(6): 2512-2514. doi: 10.1016/j.asr.2021.12.030 [14] Zhou H Y, Wang X G, Bai B, et al. Reentry guidance with constrained impact for hypersonic weapon by novel particle swarm optimization [J]. Aerospace Science and Technology, 2018, 78(6): 205. doi: 10.1016/j.ast.2018.04.024 [15] Li Z H, Hu C, Ding C, et al. Stochastic gradient particle swarm optimization based entry trajectory rapid planning for hypersonic glide vehicles [J]. Aerospace Science and Technology, 2018, 76(5): 176. [16] Cheng Z, Zhang L Q, Li X L, et al. Chaotic hybrid particle swarm optimization algorithm based on Tent map [J]. Systems Engineering and Electronics, 2007, 29(1): 103-106. (in Chinese) [17] Wu K, Tan X C. Path planning of UAV based on im-proved whale optimization algorithm [J]. Acta Aeronautica et Astronautica Sinica, 2020, 41(S2): 724286. (in Chinese) [18] Huang Q B, Li J X, Song C N, et al. Whale optimization algorithm based on cosine control factor and polynomial mutation [J]. Control and Decision, 2020, 35(3): 559-568. (in Chinese) [19] Gao Y, Cai G B, Xu H, et al. Reentry maneuver guidance of hypersonic glide vehicle under virtual multi-tentacle detection [J]. Acta Aeronautica et Astronautica Sinica, 2021, 41(11): 623703. (in Chinese) [20] Liu T Y, Liu H, Zou J. Modified calculation of dynamic infrared radiation for hypersonic reentry target [J]. Laser and Infrared, 2021, 51(3): 328-332. (in Chinese) [21] Xu H, Cai G B, Zhang S X. Modified aerodynamic coefficient fitting models of hypersonic gliding vehicle in reentry phase [J]. Journal of Astronautics, 2021, 42(9): 1139. (in Chinese) [22] Zhou H Y, Wang X G, Zhao Y L, et al. Online guidance for aerospace vehicle in return-gliding phase [J]. Journal of Astronautics, 2021, 42(2): 175. (in Chinese)