高超声速飞行器多约束再入滑翔机动弹道优化设计

刘晓慧，聂万胜，杨新垒

（装备学院，北京，101416）

高超声速飞行器多约束再入滑翔机动弹道优化设计

刘晓慧，聂万胜，杨新垒

（装备学院，北京，101416）

建立了升力体再入滑翔飞行器的气动模型和多约束模型。多约束模型除了包括热流密度、气动过载、动压和终端约束等典型约束外，还建立了更符合实际任务的路径点和禁飞区约束模型，并利用路径点、禁飞区和终端约束划分弹道，在各段分别使用高斯伪谱法进行弹道求解，将多段多约束的最优控制问题转换为非线性规划问题。改进的准平衡滑翔条件保证了弹道平缓。最后通过Matlab仿真计算验证了所用分段高斯伪谱法规划弹道比传统的高斯伪谱法具有更精确的优化结果和更高的优化效率。

高超声速飞行器；机动弹道；再入滑翔；分段高斯伪谱法

0 引 言

实际任务往往要求高超声速飞行器能够平缓再入大气层和大范围滑翔机动飞行。再入段弹道平缓可便于控制系统操控[1]，而滑翔机动飞行则可以规避拦截和进行防御，因此有必要根据实际需求进行高超声速飞行器再入滑翔机动弹道的优化设计。高超声速飞行器的弹道优化涉及众多非线性运动方程和约束条件，是一种复杂的非线性多约束最优控制问题[2]，需要通过数值方法进行计算。目前数值方法主要有两类：

a）间接优化目标函数，通过变分法或最小值原理将最优控制问题转换为汉密尔顿边界值问题（Hamilton Boundary Value Problem，HBVP）进行求解。HBVP计算精确度高，但对未知边界条件极其敏感而且求解复杂困难[3]。

b）直接优化目标函数，如直接打靶法和伪谱法，通过在特定时间节点对微分方程进行离散，将最优控制问题转换为非线性规划（Nonlinear Programming，NLP）问题进行数值求解。NLP操作简易，对复杂问题有更高的优化效率[4]，但因未计算共态变量，求解精确度相对会下降。

很多学者就此类问题进行了研究。Shen Zuojun使用准平衡滑翔条件，将包括热流密度、气动过载和动压在内的路径约束转换为倾侧角约束[5]。Reddien首次使用高斯伪谱法（Gauss Pseudospectral Method，GPM）求解最优控制问题，使用全局正交多项式在特定节点同时离散状态量和控制量[4]。Benson详细推导了积分和微分形式的GPM[6]。Huntington改进了GPM，提出了计算边界控制量的新方法[7]。Jorris和 Cobb使用GPM 首次研究了路径点和禁飞区约束下高超声速飞行器的三维弹道优化[8]。

在实际任务中，飞行器往往需要经过特定任务路径点并顺利避开敌方拦截或勿入区域。因此，除了典型的热流密度、气动过载、动压、终端条件等约束外，还要考虑路径点和禁飞区约束，这给传统的GPM求解带来一定难度：

a）到达路径点的时刻未知，不利于状态量和控制量离散节点的选取；

b）离散节点的分布在整个积分的两端密集，中间稀疏，拟合结果会降低禁飞区附近的精确度；

c）多阶次的插值多项式也会导致状态量和控制量出现较大的拟合误差。

本文提出的分段GPM，利用指定路径点和禁飞区拟接触点来划分弹道，能够很好地解决上述难点，有效优化高超声速飞行器再入滑翔弹道并满足包括典型约束和路径点、禁飞区约束在内的多种复杂约束。

1 数学模型

1.1 飞行器气动模型

飞行器采用升力体构型，气动外形如图1所示。

对于高超声速飞行器，其入口条件取来流状态参数（速度、温度、压强），出口条件采用外推方式获得，飞行器表面满足无穿透、无滑移的壁面边界条件[9]。

查标准大气属性表得到 35 km处的大气压强为558.87 Pa，温度为237 K，密度为8.2×10-3kg/m3，选用标准k-ε湍流模型分别计算得到Ma为6、8、10和12时飞行器升力系数CL、阻力系数CD和升阻比L/D随攻角α的变化曲线，如图2所示。

由图 2可知，最大升阻比出现在α为 6°左右且Ma在10以上，CL、CD和L/D随Ma变化极小，可近似认为只与α有关，与高马赫数无关理论一致[10]。由于弹道再入滑翔段的飞行速度是Ma为10左右，因此可将CL和CD分别拟合为α的一次函数和二次函数：

同时采用美国1976标准大气模型来计算随高度变化的地球重力加速度g、大气密度ρ、温度T和声速c。

1.2 飞行器滑翔动力学模型

由于高升阻比飞行器其气动力和地球重力分别约为科氏惯性力和地球自转引起的惯性离心力的100倍和1 000倍，因此可以假设地球为非旋转椭球体，坐标系原点O位于飞行器质心，OX轴沿飞行速度方向，OY轴位于包含速度向量的纵向对称平面内并垂直于OX轴指向上，OZ轴与OXY平面垂直构成右手准则，飞行器的三自由度质点动力学模型如下[11]：

式中V为飞行器速度；m为飞行器质量，m=907 kg；D为飞行器气动阻力，D= 0.5ρV2SmCD；L为飞行器气动升力，L= 0.5ρV2SmCL；Sm为飞行器参考面积，Sm= 0.48 m2；θ为飞行路径角；ϕ为飞行航向角（定义从北向南顺时针方向为正）；β为飞行倾侧角；λ和ψ分别为地球经度和纬度；H为飞行高度；Re为标准地球平均半径，Re= 6 478 km；r为飞行器到地心的径向距离，r=H+Re；LV为飞行航程。对于在临近空间飞行的高超声速飞行器，其飞行高度H远小于地球平均半径Re，因此航程的导数可进一步简化为[12]

1.3 约束模型

1.3.1 控制约束

以飞行攻角α和倾侧角β为控制变量，为使飞行器平稳飞行，有如下约束方程：

1.3.2 终端约束

为保证在一定的落角（终端路径角）和速度下将飞行器引导至指定终端目标位置，约束方程为

式中LVf表示指定航程；Hf为终端高度；λf和ψf分别为终端目标的经度和纬度；θf为指定终端路径角；Vf为终端速度。

为满足末制导阶段初始条件，需要对滑翔飞行的终端高度和速度进行约束，本文假设高度约束为35 km以上，速度约束为3 km/s。

1.3.3 路径约束

a）物理量约束。

热流密度气动过载n和动压q的约束方程如下：

式中k为常数，k= 5.21×10-8；为最大热流密度，根据文献[13]，取= 300 W/cm2；nmax为最大气动过载，nmax= 4g；qmax为最大动压，qmax= 40 kPa。

b）路径点约束。

设飞行器在时刻ti到达第i个路径点时的坐标为（x(ti),y(ti)），指定路径点的坐标为（xi,yi），则约束方程为[11]

c）禁飞区约束。

假设禁飞区为圆柱型区域，第j个禁飞区中心位置的坐标为（xcj,ycj），其半径为Rj，飞行器到第j个禁飞区中心的距离为 ∆xj=x(t) −xcj， ∆yj=y(t) −ycj。则规避禁飞区的约束方程为

2 最优弹道求解方法

2.1 准平衡滑翔条件的改进

由于弹道再入滑翔段飞行路径角较小且变化相对较缓，因此可假设cosθ= 1，，将式（3）的第2个方程化简为

式（10）称为准平衡滑翔条件（Quasi-Equilibrium Glide Condition，QEGC）[14]。理论上，只要气动升力足够，式（10）能够使弹道平缓沿直线前行。但实际上，飞行器再入后的速度会随着时间变小，为了保持准平衡滑翔，就需要更大的攻角，这样气动升力才能够抵消掉重力的影响。然而，攻角达到其上限后将不能再增大，致使可用的升力不足以保持准平衡滑翔，使得，θ迅速减小，这将导致严重不稳定[15]。因此对式（10）QEGC进行了一定改进后得到：

式中ε为一极小负预选值。这样，可以使θ在整个滑翔段稳定平缓地减小，即使攻角α达到上限，也能够保证θ减小得更少。

由于CL是α的函数，QEGC建立了α和β两个控制变量间的关系。已知β，则α可用β表示，反之亦然。很明显，QEGC减少了优化变量，提高了优化效率。

2.2 高斯伪谱法及其改进

2.2.1 传统高斯伪谱法

弹道优化属于最优控制问题，可通过数值模拟的方法确定近似解，将连续性最优控制问题转换为NLP问题。其中一种常用方法就是GPM，传统的GPM求解优化问题的基本思路如下[4]。

首先定义状态量为x(t) ∈Rn，控制量为u(t) ∈Rm，初始时间为t0，最终时间为tf，则在广义区间τ∈ [τ0,τf] =[− 1,1]上，有如下转换关系：

最优控制问题要求通过控制变量使目标函数最小。选取N个勒让德-高斯（LG）点τk（k= 1,…,N），加上2个边界点τ0和τf，总共N+2个离散点。因此通过N+1个拉格朗日插值多项式Ls来逼近状态量，同时通过N个拉格朗日插值多项式Lc来逼近控制量：

将式（13）微分求导，得到：

每个拉格朗日多项式在 LG点的导数都可表示为微分形式的近似矩阵D∈RN(N+1)，该矩阵的各项为

通过该矩阵将动力学约束改写为代数约束：

终端状态Xf可通过高斯积分用Xk≡X(τk) ∈Rn，Uk≡U(τk)∈Rm和初始状态X0表示：

目标函数J的积分部分也用高斯积分近似：

式中ωk为高斯权重。

边界约束为

LG点的路径约束为

至此，目标函数和代数约束转换为NLP问题。

2.2.2 分段高斯伪谱法

分段高斯伪谱法在传统GPM的基础上加以改进，以适用于路径点和禁飞区约束的弹道优化，基本思路为：首先根据指定路径点和禁飞区拟接触点将整个弹道规划分为多个阶段；在每个阶段使用传统GPM将状态量和控制量分别在 LG点离散，转换为代数约束；终端状态由初始状态和高斯积分表示；相邻两个阶段之间满足连续性连接条件和路径点及禁飞区约束条件，保证状态量和控制量的连续性；将所有阶段组合共同确定目标函数达到全局最优。

3 仿真计算与分析

设飞行器初始及滑翔段终端条件、对应的路径点和禁飞区约束条件分别如表1、表2和表3所示。

表1 飞行器初始端和滑翔段终端条件

表2 路径点约束条件

表3 禁飞区约束条件

3.1 无路径点和禁飞区约束的弹道

作为对比，先忽略路径点和禁飞区约束，目标函数为J= max(tf)，采用30个LG点来离散最优控制问题，使用传统GPM求解得到的三维弹道如图3所示。

计算结果显示，使用传统GPM能够得到满足终端约束条件的再入滑翔段弹道：滑翔终端位置为[56.331 2°W，29.680 3°N]，高度为35.012 km，速度为 3 007.9 m/s；落点位置为[25.095 2°W，25.001 3°N]，速度为800.4 m/s，飞行路径角65.22°，航向角-31.26°，总飞行时间为2 369 s。从图3的θ曲线可以看出，在改进的QEGC下，再入和滑翔段θ变化很小，弹道基本保持平缓，便于控制系统操控，利于减少热烧蚀。

3.2 满足路径点约束的弹道

在第3.1节的算例中，飞行器的初始端和滑翔终端条件不变，在其间增加2个路径点约束，进而将整个弹道分为初始端至路径点1、路径点1至路径点2、路径 点2至滑翔终端、滑翔终端至目标4段，相邻两段中，前一阶段的最终状态值即为后一阶段的初始状态值，并通过路径点约束将相邻2个阶段的状态连接起来。

先使用传统GPM求解弹道，采用30个LG点进行离散；再使用分段GPM，将弹道分为4段，每段采用7个LG点进行离散，得到的弹道如图4和图5所示。

两种方法均能得到经过路径点且满足终端约束条件的再入滑翔段弹道。使用传统 GPM的弹道的路径点1实际位置为[99.987 0°W，33.485 9°N]，路径点2实际位置为[70.008 3°W，33.989 0°N]；滑翔终端位置为[58.618 5°W，34.263 7°N]，高度为34.959 km，速度为 2 993.1 m/s；落点位置为[25.078 2°W，25.021 5°N]，速度为 797.0 m/s，飞行路径角22.41°，航向角-19.11°，总飞行时间为2518 s。使用分段GPM的弹道的路径点 1实际位置为[99.999 8°W，33.500 1°N]，路径点 2实际位置为[69.999 8°W，34.100 2°N]；滑翔终端位置为[56.250 0°W，33.976 7°N]，高度为35.006 km，速度为3 092.7 m/s；落点位置为[25.340 2°W，25.011 8°N]，速度为809.9 m/s，飞行路径角63.28°，航向角-50.36°，总飞行时间为2 476 s。

相比于传统GPM求解的弹道，分段GPM选择了更少的 LG点，且只在初始位置、终端位置和路径点位置分布密集，使得弹道能够更加严格的满足约束条件，且滑翔终端速度更大，落点误差更小，飞行时间更短，因此有效的提高了优化效率。另外，从图5中θ与ϕ曲线的对比可以看出，使用分段GPM求解得到的飞行路径角变化更小，弹道更为平缓。

3.3 满足路径点约束和禁飞区约束的弹道

在第 3.2节的算例中，分别在初始端至路径点 1和路径点2至滑翔终端2段中加入2个禁飞区。同样，先使用传统GPM求解弹道，采用30个LG点进行离散；再使用分段GPM，将弹道分为6段，每段采用5个LG点进行离散，得到结果如图6所示。

计算结果显示，传统GPM未能避开禁飞区，无法 求解到满足禁飞区约束条件的弹道；分段GPM则能够顺利求解到符合各种约束条件的弹道，飞行经过2个路径点，顺利避开了2个禁飞区并紧贴其边缘以减少额外能耗。禁飞区1的接触点位置为[107.100 3°W，33.010 9°N]，禁飞区2的接触点位置为[77.155 5°W，33.823 7°N]；路径点1实际位置为[100.001 5°W，33.499 6°N]，路径点 2实际位置为[70.000 9°W，34.099 5°N]；滑翔终端位置为[55.328 1°W，33.910 9°N]，高度为35.015 km，速度为3176.1 m/s；落点位置为[西经25.035 6°W，北纬24.994 2°N]，速度为 814.4 m/s，飞行路径角 53.31°，航向角为-58.03°，飞行总时间为2 597s。

图7为使用分段GPM求解得到的路径点和禁飞区约束弹道的飞行路径角θ、航向角ϕ、飞行速度V、热流密度气动过载n及动压q的变化情况。

图7中热流密度在飞行初始段和飞行末段均控制在较低水平，这是因为，初始高度处大气极其稀薄，ρ很小，而飞行末段V大大降低，在滑翔段出现最大值，接近，即达到最大值；相应地，由图2可知，在飞行初始段先增大到最大值16°后迅速减小并保持在最大升阻比攻角6°附近，就是为了保证弹道平缓，减少热烧蚀，控制热流密度。气动过载n一直控制在较低水平，直到飞行末段接近最大值，这是由末制导段飞行高度下降快且大气密度急剧增大导致的，动压q亦同理为逐渐增大并在飞行末段达到最大值；分析对应攻角变化，在飞行末段，α先迅速增大以满足n和q的约束，之后迅速减小至0.76°，以压低弹道，保证攻击目标时具有一定的落角。通过分析可以发现，热流密度对飞行器滑翔段的约束更强，气动过载和动压则是对飞行末段的约束更强。

4 结 论

本文在准平衡滑翔假设下，提出了一种分段高斯伪谱法来优化高超声速弹道，实现平缓再入和机动滑翔，并且满足包括典型约束、路径点约束和禁飞区约束在内的多种约束条件。该方法优点如下：

a）在改进的准平衡滑翔条件假设下，再入段弹道平缓，避免了高度振荡，便于控制系统操控。

b）使用分段 GPM，根据指定路径点将弹道划分为多段，无需考虑到达每个路径点的时刻，而是将路径点视作每段的终端约束条件，再使用GPM求解。相比传统GPM，分段GPM可以更少的离散节点完成弹道求解，降低了插值多项式的阶次，进而减小了拟合误差，优化结果更精确，优化效率更高。

c）使用分段 GPM，根据禁飞区的拟接触点来划分弹道，在接触点附近的离散节点更密集，弹道平滑，可顺利避开禁飞区并通过路径点，且各种典型约束都能控制在较低水平，而传统GPM则无法得到可行解。

[1] 和争春, 何开锋, 朱国林. 过载限制条件下升力体再入大气层弹道优化设计[J]. 弹道学报, 2007,19(2): 8-12.

[2] 张烨琛, 包为民, 闵昌万, 张宁宁. 高超声速无动力滑翔飞行器最优弹道研究[J]. 导弹与航天运载技术, 2014(4): 5-8.

[3] 李佳峰, 周浩, 陈万春. 高超声速飞行器滑行段最优弹道的间接算法[J].飞行力学, 2009,27(3): 41-49.

[4] Reddien G W. Collocation at Gauss points as a discre-tization in optimal control[J]. Control Optim, 1979, 17(2): 298-306.

[5] Shen Zuojun. Onboard generation of three-dimensional constrained entry trajectories[J]. Guidance, Control and Dynamics, 2003,26(1): 111-121.

[6] Benson D A. Direct trajectory optimization and costate estimation via an orthogonal collocation method[J]. Guid Control Dyn, 2006,29(6): 1435-1440.

[7] Huntington G T. Advancement and analysis of a Gauss pseudospectral transcription for optimal control prob-lems dissertation[D]. Massachusetts Institute of Technology, 2007.

[8] Jorris T R, Cobb G C. Three-dimensional trajectory op-timization satisfying waypoint and no-fly zone con-straints[J]. Guidance，Control and Dynamics, 2009,32(2): 543-553.

[9] 丁洪波, 蔡洪, 张士峰, 李安梁. 高起声速滑翔式再入飞行器最大航程飞行轨迹分析[J]. 国防科技大学学报, 2009,31(6): 67-72.

[10] Philips T H. A common aero vehicle (CAV) model, de-scription, and employment guide[R]. Schafer Corporation for AFRL and AFSPC, 2003.

[11] 黄育秋, 何麟书. 升力式再入飞行器多约束多阶段弹道优化设计[J].哈尔滨工业大学学报, 2011,43(7): 144-148.

[12] 方群, 李新三. 临近空间高超声速无动力滑翔飞行器最优轨迹设计及制导研究[J]. 宇航学报, 2008,29(5): 1485-1491.

[13] 刘欣, 李建林, 葛健全, 杨涛. 滑翔式飞行器再人弹道设计[J]. 弹箭与制导学报, 2011,31(6): 161-164.

[14] 张科南, 周浩, 陈万春. 高超声速飞行器多约束多种机动突防模式弹道规划[J]. 弹道学报, 2012,23(3): 85-90.

[15] 陈小庆, 侯中喜, 刘建震. 高超声速滑翔飞行器弹道特性分析[J]. 导弹与航天运载技术, 2011(2): 5-9.

Optimization and Design of Reentry Gliding Maneuvering Trajectory with Multi-Constraints for Hypersonic Vehicle

Liu Xiao-hui, Nie Wan-sheng, Yang Xin-lei
(Equipment Academy, Beijing, 101416)

The aerodynamics model of the reentry gliding vehicle and the multi-constraints model including heating rate, overload, dynamic pressure, way point and no-fly zone constraints are built. The study is based on the assumption of quasi-equilibrium glide condition to keep the trajectory smooth. The trajectory is divided into multi-phases according to the waypoints, no-fly zones and terminal constraints and the optimal trajectory is solved by employing multi-phases gauss pseudospectral method, which transfers the multi-constraints and multi-phases trajectory optimization problem into a nonlinear programming problem. The simulations can well demonstrate that the improved method in this paper is more precise and efficient for the optimization than the traditional method.

Hypersonic vehicle; Maneuvering trajectory; Reentry gliding; Multi-phases Gauss pseudospectral method

V421

A

1004-7182(2017)02-0006-06

10.7654/j.issn.1004-7182.20170202

2016-07-08；

2017-01-10

刘晓慧（1991-），女，硕士研究生，主要研究方向为高超声速飞行器总体技术