在未来军事作战环境中，超音速无人机拥有出色的机动性能及突防能力，可以替代有人战斗机快速执行对地攻击和目标侦查的高风险任务，得到了国内外许多研究者的关注^[1-3]. 随着军事和科学应用研究的深入，将超音速飞机小型化，作为低成本的飞行试验系统变得至关重要^[4-5]. 对于小型的超音速无人机而言，由于机身空间结构的限制并兼顾良好的经济性能，对气动外形设计与推进系统的要求极高. 当接近声速时，激波阻力急剧增加，涡喷发动机推力有限，跨音速加速性能较差，则通过规划在该飞行阶段的合适航迹，将在很大程度上缓解对推进系统严苛要求的挑战^[6-9]. 从经济性能角度来看，选择合适推力且不配备加力燃烧室的涡喷发动机是低成本小型超音速无人机的较优方案. 本文提出配置高推重比的无加力涡喷发动机，利用马赫俯冲机动实现跨音速飞行从而突破音障的方案，可以有效降低发动机成本，提高该型无人机的经济性能. 在既定先进的气动外形构型下，研究无人机可靠且快速突破音障的航迹具有现实工程意义与广泛应用价值.

目前，轨迹规划方法主要有最优控制法和智能优化算法，其中最优控制法包括间接法和直接法^[10-11]. 间接法收敛半径小，难以估计共轭变量的初值，且超音速无人机轨迹优化问题存在多个复杂约束，将会使推导求解十分困难，故利用参数化方法将连续空间的最优控制问题转化为非线性规划问题（nonlinear program，NLP）的直接法被广泛应用于飞行器轨迹优化领域^[12-15]. 另一类研究方法是利用传统的分而治之模型与智能优化算法进行求解. 陈晓等^[16]用能量法建立爬升、平飞和下滑段轨迹优化模型，利用再开始FR（Fletcher-Revees）共轭梯度法求解最优轨迹. 安兵辉等^[3]用能量高度法构建超音速飞机进入俯冲、俯冲直线段和改出俯冲3段俯冲段模型，使用遗传算法进行轨迹规划. 运用传统分而治之模型，仅是在各段分别寻优，且这类智能优化算法在求解复杂约束问题时局部搜索能力不足，易发生早熟收敛^[11]. 近年来，结合伪谱法与hp型有限元法的hp自适应伪谱法受到广泛关注，该算法由Darby^[17]提出，在优化过程中对区间网格密度h（配点数）和多项式阶数p进行自适应调整，相较于其他伪谱法，对状态与控制变量的近似程度更高，适用于状态或控制变量变化较快的优化问题，表现出很高的效率及计算精度^[18].

为了探究小型无加力燃烧室的超音速无人机能够顺利从高亚音速逐步过渡到超音速飞行的可行轨迹，保证该飞行段中的耗油量或飞行时间最小，本文通过建立该型无人机的运动模型，考虑在该飞行段的可操作性限制、初值状态约束、控制约束及路径约束等一系列约束条件，采用hp自适应伪谱法研究跨音速飞行段的最小油耗和最小时间2种轨迹优化问题，通过与传统飞行方案所得的结果进行比较，可以为工程领域的实际应用提供重要参考.

无人机加速从高亚音速到超音速的飞行过程，一般限制在铅垂平面内，可以将无人机作为质点处理，超音速无人机在地面坐标系下的动力学模型^[19]为

(1) $ \left. \begin{array}{*{20}{l}} \dot x = v\cos \; \gamma , \hfill \\ \dot h = v\sin \;\gamma , \hfill \\ \dot v = \dfrac{{T\cos \;(\alpha + {\varphi _{\rm{T}}}) - D}}{m} - g\sin \gamma , \hfill \\ \dot \gamma = \dfrac{{T\sin \;(\alpha + {\varphi _{\rm{T}}}) + L}}{{mv}} - \dfrac{{g\cos \gamma }}{v}, \hfill \\ \dot m = - s_{\rm{fc }} T/36{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 000. \hfill \\ \end{array} \right\} $

式中： $ \left( {x,h} \right) $为超音速无人机的位置坐标； $ v$为飞行速度； $ \alpha $为迎角； $ \gamma $为航迹倾角； $ m $为无人机质量； $ g $为重力加速度； $ L、D $ 和 $ T $分别为升力、阻力和涡喷发动机推力； $ {\varphi _{\rm{T}}} $为发动机安装角，假设发动机推力方向与飞行方向一致，则 ${\varphi _{\rm{T}}} = 0$； $s_{\rm{fc}} $为发动机燃油消耗率，与飞行高度h、飞行马赫数Ma有关.

$ L $和 ${{D}}$可以表示为

(2) $ \left. \begin{array}{*{20}{l}} {{L = }}\;\rho {v^2}S{C_{\rm{L}}}/2, \hfill \\ D = \rho {v^2}S{C_{\rm{D}}}/2. \hfill \\ \end{array} \right\} $

式中： $ \;\rho $为空气密度，可以通过1976年美国标准大气模型获得对应高度的空气密度； $ S $为无人机参考面积； $ {C_{\rm{L}}} $和 $ {C_{\rm{D}}} $分别为升力系数和阻力系数.

$ {C_{\rm{L}}} $、 $ {C_{\rm{D}}} $受翼型、 $ \alpha $和Ma的影响，设 $ \alpha $在较小范围内变动，对于轴对称无人机，升力系数可以看作是 $ \alpha $的线性函数，阻力系数为 $ \alpha $的二次函数，通过下式计算：

(3) $ \left. \begin{array}{*{20}{l}} {C_{\rm{L}}}\left( {\alpha ,Ma} \right) = {C_{{\rm{L}}\alpha }}\left( {Ma} \right)\alpha , \hfill \\ {C_{\rm{D}}}\left( {\alpha ,Ma} \right) = {C_{{\rm{D}}0}}\left( {Ma} \right) + {A_i}C_{\rm{L}}^2\left( {\alpha ,Ma} \right). \hfill \\ \end{array} \right\} $

式中： $ {C_{{\rm{L}}\alpha }} $为升力线斜率； $ {C_{{\rm{D}}0}} $和 $ {A_i} $分别为零升阻力系数和诱导阻力因子，可以通过风洞试验或计算流体力学方法得到.

发动机推力、燃油消耗率与飞行高度、飞行速度、转速等密切相关，模型较复杂^[20]. 在高亚音速到超音速阶段，为了使无人机能够顺利突破音障，该阶段发动机均应处于最大转速状态. 根据发动机参数，最大转速下推力、燃油消耗率与 $ h $、Ma之间的关系如图1所示.

为了使无人机从高亚音速到超音速阶段飞行时间最小或耗油量最小，可以通过最优控制理论进行优化，优化目标主要包括以下2个方面. 1）从初始爬升点能够到达巡航飞行点；2）飞行过程中的耗油量或飞行时间最小. 综合考虑无人机的复杂约束，最优控制轨迹规划问题可以表示为如下数学形式.

1）初始爬升点的约束条件为

(4) $ \left. \begin{array}{*{20}{l}} x\left( {{t_0}} \right) = {x_0},\;h\left( {{t_0}} \right) = {h_0}, \hfill \\ m\left( {{t_0}} \right) = {m_0}, \hfill \\ Ma\left( {{t_0}} \right) = M{a_0},\;\gamma \left( {{t_0}} \right) = {\gamma _0}. \hfill \\ \end{array} \right\} $

式中： $ {t_0} $为初始爬升点时刻， $ \left( {{x_0},{h_0}} \right) $为初始爬升点坐标， $ {m_0} $为初始爬升点质量， $ M{a_0} $为初始爬升点马赫数， $ {\gamma _0} $为初始爬升点航迹倾角.

2）超音速巡航飞行点的约束条件为

(5) $ \left. \begin{array}{*{20}{l}} x\left( {{t_{\rm{f}}}} \right) = {x_{\rm{f}}},\; h\left( {{t_{\rm{f}}}} \right) = {h_{\rm{f}}}, \hfill \\ m\left( {{t_{\rm{f}}}} \right) = {m_{\rm{f}}}, \hfill \\ Ma\left( {{t_{\rm{f}}}} \right) = M{a_{\rm{f}}},\; \gamma \left( {{t_{\rm{f}}}} \right) = {\gamma _{\rm{f}}}. \hfill \\ \end{array} \right\} $

式中： $ {t_{\rm{f}}} $为巡航飞行点时刻， $ \left( {{x_{\rm{f}}},{h_{\rm{f}}}} \right) $为巡航飞行点坐标， $ {m_{\rm{f}}} $为巡航飞行点质量， $ M{a_{\rm{f}}} $为巡航飞行点马赫数， $ {\gamma _{\rm{f}}} $为巡航飞行点航迹倾角.

3）控制约束条件。涡喷发动机在跨音速飞行工况下应均处于最大转速状态，可以不考虑发动机转速n的影响。无人机失速、过载限制的飞行特性与 $ \alpha $密切相关，且工程应用中易于测量，故设定 $ \alpha $为控制变量。考虑到 $ \alpha $在较小范围内变动，设定控制变量的最大值为 $ {\alpha _{\max }} $.

(6) $ \left| \alpha \right| \leqslant {\alpha _{\max }}. $

4）实时路径约束条件。无人机作马赫俯冲机动必须严格满足实时路径约束，否则生成的轨迹不可行，甚至可能会影响机体的结构与操纵安全^[21].

(7) $ \left.\begin{split} & {h_0} \leqslant h(t) \leqslant {h_{\max }},M{a_0} \leqslant Ma(t) \leqslant M{a_{\rm{f}}}, \hfill \\ & {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {m_{\rm{f}}} \leqslant m(t) \leqslant {m_0}, \hfill \\ & {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left| {\gamma (t)} \right| \leqslant {\gamma _{\max }},\; |{n_y}| \leqslant {n_{y\max }}. \end{split}\right\} $

式中：t为飞行时间， $ {h_{\max }} $为爬升最大高度， $ {\gamma _{\max }} $为最大航迹倾角， $ {n_y} $为无人机法向过载， $ {n_{y\max }} $为最大法向过载.

5）目标函数. 根据建立的超音速无人机模型，选取 $ h、V$、 $ \gamma $、m为状态变量， $\alpha $为控制变量，在适当的目标函数和必要的约束条件下，建立跨音速飞行段轨迹优化模型. 选择的目标函数即式（8）、（9），分别给出最小燃料和最小爬升时间轨迹优化问题，则燃料消耗量即为 $ {m_{{\rm{FC}}}} = m\left( {{t_0}} \right) - m\left( {{t_{\rm{f}}}} \right) $. 最省燃料是侦察型超音速无人机追求的一个重要指标，关系到经济性能；最小爬升时间是攻击型无人机的关键要素，对快速打击敌方目标至关重要.

(8) $ {J_{\min ,{\text{fuel}}}} = - m{\kern 1pt} {\kern 1pt} ({t_{\rm{f}}}). $

(9) $ {J_{\min ,{\text{time}}}} = \int_0^{{t_{\rm{f}}}} {{\rm{d}}t} . $

通过使用hp自适应伪谱法，将该过程的轨迹优化问题转化为带有复杂约束的非线性规划问题，利用序列二次规划法进行求解.

采用的hp自适应伪谱法是融合Gauss伪谱法与hp型有限元法，将状态变量与控制变量在Legendre-Gauss（LG）点上离散，以这些离散点为节点构造全局插值多项式，从而近似状态变量与控制变量，对插值多项式进行求导得到微分矩阵. 利用微分矩阵将动力学微分方程约束转化为代数约束，将小型超音速无人机连续时间最优控制问题转换为非线性规划问题. 设e为离散化动力学约束和路径约束的最大容许误差，若某一离散区间的容许误差高于e，则对区间网格密度h和全局插值多项式的阶次p进行自适应调整. 相比于其他伪谱法，如Legendre伪谱法、Radau伪谱法，Gauss伪谱法具有更高的求解精度和收敛速度，在处理含初始和终端约束的问题上更具优势^[22-23].

Gauss伪谱法对超音速无人机马赫俯冲机动轨迹优化问题的离散过程如下. 将最优控制问题的时间区间 $ \left[ {{t_0},{t_{\rm{f}}}} \right] $划分为 $ K $个网格，并转换到区间 $ \left[ { - 1,1} \right] $，进行映射变换，有

(10) $ t = \frac{{{t_k} - {t_{k - 1}}}}{2}\tau + \frac{{{t_k} + {t_{k - 1}}}}{2}; \; {{t_k} > {t_{k - 1}}} .$

式中： $t\in \left[{t}_{k-1},{t}_{k}\right]，\tau \in \left[-1,1\right]，k=1,\cdot \cdot \cdot ,K,\;{t}_{0} < \cdot \cdot \cdot < $ $ {t_K} = {t_{\rm{f}}} .$

在Gauss伪谱法的开区间 $ \left( { - 1,1} \right) $上取N个离散点（LG配点），即 $ \left\{ {{\tau _1},{\tau _2}, \cdot \cdot \cdot ,{\tau _i}, \cdot \cdot \cdot ,{\tau _N}} \right\} $ $ \left( {i = } \right. $ $ \left. {1,2, \cdot \cdot \cdot ,N} \right) $，这些LG配点是N阶Legendre多项式 $ {G_N}\left( \tau \right) $ 的零点，即

(11) $ {G_N}\left( \tau \right) = \frac{1}{{{2^N}N!}} \cdot \frac{{{{\rm{d}}^N}}}{{{\rm{d}}{\tau ^N}}}\left[ {{{\left( {{\tau ^2} - 1} \right)}^N}} \right]. $

由于研究对象关乎到初始约束和终端约束，须引入初始时刻和终端时刻2个节点. 将状态变量在前N+1个LG配点处离散，以N+1阶Lagrange插值多项式 $ L_j^{\left( k \right)}\left( \tau \right) $作为基函数，将状态变量近似为

(12) $ \left.\begin{split} & {x^{\left( k \right)}}\left( \tau \right) \approx {X^{\left( k \right)}}\left( \tau \right) = \displaystyle\sum\limits_{j = 0}^N {X_j^{^{\left( k \right)}}} \left( \tau \right)L_j^{^{\left( k \right)}}\left( \tau \right), \hfill \\ & L_j^{^{\left( k \right)}}\left( \tau \right) = \displaystyle\prod\limits_ {\begin{subarray}{l} l = 0 \\ l \ne j \end{subarray} } ^N {\frac{{\tau - \tau _l^{^{\left( k \right)}}}}{{{\tau _j} - \tau _l^{^{\left( k \right)}}}}} ; \;j = 0,1, \cdot \cdot \cdot ,N. \end{split} \right\} $

类似地，将控制变量在N个LG配点处离散，以N阶Lagrange插值多项式近似得到 $ {u^{\left( k \right)}}\left( \tau \right) $.

(13) $ \left.\begin{split} & {u^{\left( k \right)}}\left( \tau \right) \approx {U^{\left( k \right)}}\left( \tau \right) = \sum\limits_{j = 1}^N {U_j^{^{\left( k \right)}}} \left( \tau \right)\tilde L_j^{^{\left( k \right)}}\left( \tau \right), \hfill \\ & \tilde L_j^{^{\left( k \right)}}\left( \tau \right) = \prod\limits_{\begin{subarray}{l} l = 1 \\ l \ne j \end{subarray} }^N {\frac{{\tau - \tau _l^{^{\left( k \right)}}}}{{{\tau _j} - \tau _l^{^{\left( k \right)}}}}} ; \;j = 1,2, \cdot \cdot \cdot ,N. \end{split} \right\}$

式（12）未定义终端时刻状态，且初始时刻已知 $ X\left( {{\tau _0}} \right) = x\left( {{t_0}} \right) $，可以用积分运算离散化来近似，即

(14) $ \begin{split} x\left( {{\tau _{\rm{f}}}} \right) \approx & X\left( {{\tau _{\rm{f}}}} \right) = X\left( {{\tau _0}} \right) + \frac{{{t_k} - {t_{k - 1}}}}{2}\sum\limits_{k = 1}^K {\sum\limits_{i = 1}^N {\omega _i^{\left( k \right)}} } \times \hfill \\ & f\left[ {{X^{\left( k \right)}}\left( {\tau _i^{}} \right),{U^{\left( k \right)}}\left( {\tau _i^{}} \right),\tau ;\;{t_{k - 1}},{t_k}} \right]. \end{split} $

式中： $ {\tau _i} $ $ \left( {i = 1,2, \cdot \cdot \cdot ,N} \right) $为LG点， $ \omega _i^{\left( k \right)} $为高斯权重因子.

将状态变量的时间导数近似为

(15) $\begin{split} & \sum\limits_{j = 0}^N {\dot L_j^{\left( k \right)}} \left( {{\tau _i}} \right)x_j^{\left( k \right)}\left( {{\tau _i}} \right) = \frac{{{t_k} - {t_{k - 1}}}}{2}{V^{\left( k \right)}}\left( {{\tau _i}} \right)\cos\; {\gamma ^{\left( k \right)}}\left( {{\tau _i}} \right), \hfill \\ & \sum\limits_{j = 0}^N {\dot L_j^{\left( k \right)}} \left( {{\tau _i}} \right)h_j^{\left( k \right)}\left( {{\tau _i}} \right) = \frac{{{t_k} - {t_{k - 1}}}}{2}{V^{\left( k \right)}}\left( {{\tau _i}} \right)\sin\; {\gamma ^{\left( k \right)}}\left( {{\tau _i}} \right), \hfill \\ & \sum\limits_{j = 0}^N {\dot L_j^{\left( k \right)}} \left( {{\tau _i}} \right)V_j^{\left( k \right)}\left( {{\tau _i}} \right) = \frac{{{t_k} - {t_{k - 1}}}}{2}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \times \hfill \\ &\qquad \left( {\frac{{{T^{\left( k \right)}}\left( {{\tau _i}} \right)\cos \;{\alpha ^{\left( k \right)}}\left( {{\tau _i}} \right) - {D^{\left( k \right)}}\left( {{\tau _i}} \right)}}{m} - g\sin \;{\gamma ^{\left( k \right)}}\left( {{\tau _i}} \right)} \right), \hfill \\ & \sum\limits_{j = 0}^N {\dot L_j^{\left( k \right)}} \left( {{\tau _i}} \right)\gamma _j^{\left( k \right)}\left( {{\tau _i}} \right) = \frac{{{t_k} - {t_{k - 1}}}}{2}\left( {\frac{{{T^{\left( k \right)}}\left( {{\tau _i}} \right)\sin {\alpha ^{\left( k \right)}}\left( {{\tau _i}} \right)}}{{mg}}} \right. + \hfill \\ &\qquad {\kern 1pt} \left. {\frac{{{L^{\left( k \right)}}\left( {{\tau _i}} \right)}}{{mg}} - \cos {\gamma ^{\left( k \right)}}\left( {{\tau _i}} \right)} \right)\frac{g}{{{V^{\left( k \right)}}\left( {{\tau _i}} \right)}}, \hfill \\ & \sum\limits_{j = 0}^N {\dot L_j^{\left( k \right)}} \left( {{\tau _i}} \right)m_j^{\left( k \right)}\left( {{\tau _i}} \right) = - \frac{{{t_k} - {t_{k - 1}}}}{2} \times \hfill \\ & \qquad {\kern 1pt} {s_{\rm{fc}}^{\left( k \right)}}\left( {{\tau _i}} \right) {T^{\left( k \right)}}\left( {{\tau _i}} \right)/36{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 000. \end{split} $

将实时路径约束离散化，可得

(16) $ C\left[ {{X^{\left( k \right)}}\left( {\tau _i^{}} \right),{U^{\left( k \right)}}\left( {\tau _i^{}} \right),\tau ;\;{t_{k - 1}},{t_k}} \right] \leqslant 0. $

将目标函数离散化，积分项用数值积分求解，可得如下形式：

(17) $ \begin{split} & J = \varPhi \left( {X\left( {{\tau _0}} \right),{t_0},X\left( {{\tau _{\rm{f}}}} \right),{t_{\rm{f}}}} \right) + \frac{{{t_k} - {t_{k - 1}}}}{2}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \times \hfill \\ & \sum\limits_{k = 1}^K {\sum\limits_{i = 1}^N {\omega _i^{\left( k \right)}} L\left[ {{X^{\left( k \right)}}\left( {\tau _i^{}} \right),{U^{\left( k \right)}}\left( {\tau _i^{}} \right),\tau ;\;{t_{k - 1}},{t_k}} \right]} . \end{split} $

为了确定网格是否应细化或者近似多项式的次数是否应增加，须计算第 $ k $个网格区间内第 $ i $个状态分量近似值 $ {X^{\left( k \right)}}\left( {\tau _i^{}} \right) $的曲率 $ {\kappa ^{\left( k \right)}}\left( \tau \right) $为

(18) $ {\kappa ^{\left( k \right)}}\left( \tau \right) = \frac{{\left| {{{\ddot X}^{\left( k \right)}}\left( {\tau _i^{}} \right)} \right|}}{{\left| {{{\left[ {1 + {{\dot X}^{\left( k \right)}}{{\left( {\tau _i^{}} \right)}^2}} \right]}^{3/2}}} \right|}}. $

将最大曲率 $\kappa _{\max }^{\left( k \right)} $与平均曲率 $ {\bar \kappa ^{\left( k \right)}} $之比定义为

(19) $ {r_k} = {{\kappa _{\max }^{\left( k \right)}}}/{{{{\bar \kappa }^{\left( k \right)}}}}. $

若 ${r_k} > {r_{k\max }}$，则网格间隔相对于其余的网格间隔存在大的曲率，应细化网格，增加网格的配点数目；若 ${r_k} < {r_{k\max }}$，则应增加多项式的阶次，以在网格区间上获得更好地逼近. ${r_{k\max }} > 0$为用户自定义常数.

更新后第 $ k $个网格区间内多项式的阶次 $ {N_k} $由下式决定：

(20) $ {N_k} = {M_{\rm{d}}} + {\rm{ceil}}\;\left[ {\lg\; \left( {e_{\max }^{\left( k \right)}/e} \right)} \right] + {A_{\rm{d}}}. $

式中： ${M_{\rm{d}}}$为网格中多项式的初始阶次， ${A_{\rm{d}}} > 0$为控制网格中配置点数目增长的整数常数， $ e_{\max }^{\left( k \right)} $为最大误差， ${\rm{ceil}}\left( \cdot \right)$为舍入到下一个最高整数的运算符.

更新后整条轨迹的网格间隔新数量 $ {n_k} $由下式决定：

(21) $ {n_k} = {\rm{ceil}}\;\left[ {{c_{\rm{e}}}\lg\;\left( {e_{\max }^{\left( k \right)}/e} \right)} \right].$

式中： ${c_{\rm{e}}}$为控制网格间隔数量增长的整数常数. 使用曲率密度函数的积分来确定新网格点的位置，密度的函数表达式为

(22) $ p\left( \tau \right) = {c_\tau }\kappa {\left( \tau \right)^{1/3}}. $

常值参数 $ {c_\tau } $满足：

(23) $ \int_{ - 1}^1 {p\left( \zeta \right)} {\rm{d}}\zeta = 1. $

则累积分布函数为

(24) $ F\left( \tau \right) = \int_{ - 1}^\tau {p\left( \zeta \right)} {\rm{d}}\zeta . $

对应的下一迭代中 $ {n_k} $个点的位置可由下式决定：

(25) $ F\left( {{\tau _i}} \right) = \frac{{i - 1}}{{{n_k}}}; \;1 \leqslant i \leqslant {n_k} + 1. $

若 $ {n_k} = 1 $，则不创建新的子间隔. $ {n_k} $的最小值至少为2.

hp自适应伪谱法算法的流程如下.

1）形成初始化网格.

2）对每个网格 $ k = 1, \cdot \cdot \cdot ,K $使用固定阶次的近似多项式，解决产生的非线性规划（NLP）问题.

3）若 $ e_{\max }^{\left( k \right)} \leqslant e $，则继续执行下一个网格 $ k $.

4）若 $ {r_k} \geqslant {r_{k\max }} $或 ${N_k} > {M_{\rm{d}}}$，则将第 $ k $个网格间隔细化为 $ {n_k} $个子间隔，设置任一子间隔的多项式阶次为 ${M_{\rm{d}}}$，继续执行下一个网格 $ k $.

5）否则，设置第 $ k $个子间隔的多项式阶次为 $ {N_k} $.

6）在网格参数更新完毕后，返回步骤2），开展下一次迭代.

经过以上步骤，可以将跨音速飞行段的小型超音速无人机轨迹优化问题转化为一系列非线性参数优化问题，采用序列二次规划算法进行求解. 选取超音速无人机的迎角作为优化设计变量，考虑初始、终端与路径约束，结合标准大气和气动力模型，进行轨迹优化仿真. 若性能指标达到最优，则求解结束；否则利用优化方法重新进行设计，直至获得性能指标最优方案. 超音速无人机的最优轨迹设计框架如图2所示.

以某小型超音速无人机为研究对象，结合实际飞行情况，该型无人机的参数设定为： $ {m_0} = $ $420{\kern 1pt} {\kern 1pt} {\kern 1pt} {\rm{kg}}, S = $ $ 1.33{\kern 1pt} {\kern 1pt} {{\rm{m}}^2}$. 将初始约束条件设定为： ${x_0} = 0,{h_0} = 8{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 000{\kern 1pt} {\kern 1pt} {\rm{m}}, $ $ M{a_0} = 0.8,{\gamma _0} = 0$. 将终端约束条件设定为： ${h_{\rm{f}}} = $ $ 11{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 000{\kern 1pt} {\kern 1pt} {\rm{m}},M{a_{\rm{f}}} = 1.6$， $\;{\gamma _{\rm{f}}} = 0,{m_{\rm{f}}} \geqslant 300{\kern 1pt} {\kern 1pt} {\rm{kg}}$. 将控制约束条件设定为： $ \left| \alpha \right| \leqslant {4^ \circ } $. 将实时路径约束条件设定为： $ \left| \gamma \right| \leqslant {40^ \circ }$， $8{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 000{\kern 1pt} {\kern 1pt} {\rm{m}} \leqslant h\left( t \right) \leqslant 14{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 000{\kern 1pt} {\kern 1pt} {\rm{m}}$， $ 0.8Ma \leqslant Ma\left( t \right) \leqslant $ $1.6{\kern 1pt} {\kern 1pt} Ma, \;\left| {{n_y}} \right| \leqslant 2$.

小型超音速无人机在高亚音速模式向超音速模式过渡的情况下，激波阻力出现，零升阻力系数急剧增大，但突破音障后，随着Ma的增加，激波强度减弱，最小燃料及最小爬升时间的轨迹可能非常复杂. 最小燃料轨迹优化是关系无人机经济性能的关键问题，最小爬升时间轨迹优化是关系无人机快速打击的关键问题，因此分别对这两种方案下的轨迹进行优化仿真. 上述方法提及的最大容许误差e将影响计算精度与计算效率，e分别取为10⁻⁶、10⁻⁴、10⁻³和10⁻²，在最小燃料消耗情况下的仿真结果如表1所示。表中，t_c为计算时间.

随着最大容许误差要求的不断缩小， $ {t_{\rm{c}}} $变长，结果的精度没有显著提高. 当最大容许误差取10⁻²时，t有一定的偏差，故将e设置为10⁻³. 图3给出以最小燃料消耗和最小爬升时间为性能指标的最优轨迹.

沿最小燃油轨迹（minimum fuel trajectory，MFT）与沿最小爬升时间轨迹（minimum climb time trajectory，MTT）的飞行趋势相同，即无人机在跨音速段均以马赫俯冲机动达到超音速，但MFT比MTT加速爬升的高度更高，俯冲落差更大. 从图3（b）可知，无人机沿着MTT爬升得更快，需要更大的动压来提供足够的推力. 从图3（c）可知，迎角均处于正迎角状态，且在临近终端状态时，随着无人机再一次作小幅俯冲，迎角均增大. 从图3（d）可知，MFT的油耗为20.13 kg，比MTT减少了2.4 kg，降低了10.65%，但飞行时间增加了9.21%，因此MFT燃油消耗较少，可以充分发挥该型无人机高空高速巡航的性能优势，MTT可以适用于短程快速打击.

超音速无人机在跨音速段实际飞行过程中，传统的飞行方案为无人机先等速爬升至一定高度，而后平飞加速至突破音速，最后以等 $ \gamma $加速下降. 由于该飞行阶段要使无人机燃油消耗量尽可能低，高空气动阻力较小，可以将无人机等速爬升至较高处后转入平飞加速段，采用的分而治之模型因高效性与稳定性而被广泛应用. 但传统飞行方案的制定通常是经验行为，不能保证轨迹的最优性. 为了验证所研究的最小燃料及最小爬升时间的轨迹，与传统飞行方案进行比较，结果如图4所示. 图中，传统飞行方案1为平飞加速至1.4 Ma后以 $ \gamma {{ = - 2}}{{\text{0}}^ \circ } $等 $ \gamma $加速下降，传统飞行方案2为平飞加速至1.5 Ma后以 $ \gamma {{ = - 4}}{{\text{0}}^ \circ } $等 $ \gamma $加速下降.

对比图4可以看出，本文得到的最小燃料消耗及最小爬升时间轨迹更合理，性能指标均优于传统飞行方案轨迹，最小燃料消耗轨迹相较于传统飞行方案1，油耗减少了2.52 kg，降低了 11.13%. 最小爬升时间轨迹相较于传统飞行方案1，飞行时间减少了231.39 s，降低了46.03%. 由此可知，该型无人机在跨音速段通过马赫俯冲机动突破音速，相比于平飞加速到突破音速，再利用俯冲加速达到特定巡航马赫数，耗油量更少，飞行时间更短，更能够发挥该型无人机侦察与作战的优势.

基于燃油消耗最少的原则，探究该型无人机在不同起飞质量、爬升最大高度及终端飞行高度约束条件下对飞行轨迹的影响. 对于 $ {m_0} \in \left\{ {220\; {\rm{kg}},} \right. $ $ \left. {320\; {\rm{kg}}, 420\; {\rm{kg}}, 520\; {\rm{kg}}} \right\} $，各种起飞质量下的MFT轨迹曲线如图5所示.

从图5（a）、（d）可以看出，无人机的起飞质量越小，则加速爬升的高度越高，俯冲落差越大，完成相同指标所耗费的时间越少，耗油量越小. 从图5（b）可以看出，俯冲加速时马赫数变化率相较于超音速加速飞行时更大，印证了无人机在1 Ma附近时飞行状态极其复杂. 从图5（c）可以看出，起飞质量越小，则迎角的变化相对越频繁，须不断改变飞行姿态以提供足够的升力，机动性能更强. 通过讨论表明，该型无人机不同起飞质量对燃油消耗影响较大，随着起飞质量的增加，在最小燃油轨迹（MFT）上油耗显著增加，对初步设计及控制具有重要的指导作用.

为了研究MFT与爬升最大高度及终端飞行高度间的关系，对于 ${h_{\max }} \in \left\{ {11 \;{\kern 1pt}{\rm{ km}}} \right.\left. {,12\; {\kern 1pt}{\rm{ km}} {\kern 1pt} ,13 \;{\kern 1pt}{\rm{ km}} {\kern 1pt} } \right\}$以及 ${h_{\rm{f}}} \in \left\{ {11 \;{\kern 1pt}{\rm{ km}} {\kern 1pt} ,12 \;{\kern 1pt}{\rm{ km}} {\kern 1pt} } \right\}$，给出该型无人机各种爬升最大高度及终端飞行高度的优化结果，如图6所示.

对比图6可以看出，在终端飞行高度不变的情况下，爬升最大高度越高，耗油量越小. 这是因为当爬升最大高度为11 000 m时，无人机达到约1.3 Ma后作平飞加速. 由3.3节可知，平飞加速过程突破音障更加耗油，经济性能较差. 在爬升最大高度不变的情况下，终端飞行高度越高，耗油量越大. 这是因为两条轨迹前期的运动状态是基本一致的，当临近终端飞行高度时，高度越高，须继续作加速爬升；高度越低，则仅须作小幅度俯冲加速即可达到指定终端飞行高度. 由此可知，该型无人机的MFT与爬升最大高度及终端飞行高度之间密切相关，保持在一定的高度下才能使该型无人机作马赫俯冲机动完成跨音速阶段燃料更省.

本文提出小型超音速无人机在无加力状态下，利用马赫俯冲加速突破跨音速波阻后进行超音速巡航飞行的轨迹优化方法. 为了优化该型无人机的上升轨迹，在必要的约束条件下，建立最小燃料消耗和最小爬升时间的轨迹优化问题，并转化为带有非线性约束的最优控制问题，利用hp自适应伪谱法求解. 结果表明，得到的最小燃料消耗及最小爬升时间轨迹均优于传统飞行方案轨迹，最小燃油轨迹可以充分发挥该型无人机高空高速侦查的优势，最小爬升时间轨迹可以适用于短程快速打击. 研究发现，该型无人机最小燃油轨迹（MFT）与起飞质量、爬升最大高度及终端飞行高度之间密切相关，起飞重量越大，油耗显著增加，且爬升最大高度及终端飞行高度应限定在一定的高度下才能使油耗较小. 本文工作为小型超音速无人机的方案设计提供了重要的技术支撑，也为该型无人机在实际工程领域中的航迹规划研究提供了新的思路.