基于高斯伪谱法的天基再入飞行器滑翔轨迹优化设计研究

周 亮,王昊宇,尚海滨,曹 莉,兰国峰

(1. 上海机电工程研究所，上海 201109； 2. 上海航天技术研究院北京研发中心，北京 100081； 3. 北京理工大学 宇航学院，北京 100081)

0 引 言

天基再入飞行器离开天基平台后以极高的速度再入地球大气层，利用高升阻比气动特性可以在临近空间进行无动力跳跃滑翔飞行。天基再入飞行器具有响应速度快、机动性能好、突防能力强、区域覆盖广等诸多优点，可实现远程快速精确打击和力量投送，被认为是一种具有广阔应用前景的再入飞行器[1-2]。

天基平台运行于近地轨道，由于地球自转的影响，目标点大部分时间不在轨道平面内，此时飞行器可到达的最大横向距离成为评估覆盖范围的重要指标[3]。天基再入飞行器的横向距离可从两个方面进行控制：一是在再入飞行器离开天基平台到进入大气层之前的离轨段，通过施加速度切向方向的推力来改变相对轨道面的横向距离；二是在进入大气层之后的无动力滑翔阶段，依靠气动力控制飞行速度方向，从而改变横向距离。本文主要关注无动力滑翔段的再入飞行器最远横向滑翔轨迹优化问题。

飞行器再入滑翔阶段由于超高马赫数导致气动力变化较为剧烈，同时飞行器自身动力学模型是一组时变非线性微分方程组，无法得到解析解，给再入轨迹的设计和优化带来很大的难度。目前轨迹优化问题根据求解方式不同分为间接法和直接法。间接法根据庞特里亚金极大值原理，将轨迹优化问题转化为两点边值问题进行求解，LuPing[4]、李惠峰[5]等使用间接法对高马赫数轨迹优化问题进行了分析，但间接法受限于初值猜测等问题，不能很好地应用于工程之中。直接法是将轨迹优化问题进行离散化，转换为一系列的参数优化问题进行求解，李瑜[6]通过直接打靶法进行了弹道优化，涂良辉[7]使用直接配点法研究了再入轨迹优化问题，Rao[8]使用Legendre伪谱法研究了再入飞行器轨迹控制问题，雍恩米[9]讨论了高斯伪谱法在再入轨迹优化中的应用。这些研究表明直接法在高马赫数再入轨迹优化中有很大优势，其中伪谱法作为配点法的一种，相比于其他方法具有参数规模小、收敛性能好等优点，被广泛用于高马赫数轨迹优化问题。尚海滨[10]提出求解解析雅可比矩阵可有效提高高斯伪谱法收敛速度和精度。

目前再入飞行器多采用升力体形式进行轨迹优化，本文针对轴对称再入飞行器给出了气动力分解方式，并采用高斯伪谱法对再入轨迹进行优化，同时给出了计算过程中的解析雅可比矩阵的求解方式，提高了算法的收敛精度和收敛性，最终得到了满足约束的横向距离最大再入轨迹。

1 轨迹优化问题的描述

高马赫数飞行器再入大气层后以无动力滑翔方式飞行，主要受重力和空气动力的影响。为方便描述飞行状态，在半速度坐标系中进行再入轨迹分析。该坐标系原点在再入飞行器的质心o，ox轴沿再入飞行器的速度方向，oy轴在速度矢量v与地心构成的平面内，垂直于ox轴，向上为正，oz轴垂直于xoy平面。

假设地球为匀速自转的标准球体，则轴对称再入飞行器在半速度坐标系中的动力学模型可以表示为[11]

(1)

式中：r为再入飞行器质心和地心的距离；γ和φ分别为再入飞行器的地心经度和纬度；v为再入飞行器速度；θ为弹道倾角(再入飞行器速度方向和当地水平面夹角)；σ为弹道偏角(再入飞行器速度方向与正北夹角，顺时针为正)；we为地球自转速度；Fx、Fy和Fz分别为气动阻力、升力和侧向力，控制量攻角α、侧滑角β隐含在其中，计算方式为

(2)

式中：S为轴对称飞行器特征面积；ρ为飞行器飞行高度处的大气密度；Cx、Cy和Cz分别为气动阻力系数、升力系数和侧向力系数。

气动力系数是马赫数M、攻角α、侧滑角β和高度h的函数，下面给出其推导方式。

定义飞行器本体纵轴指向头部方向与飞行器速度方向之间的夹角η为总攻角，则阻力系数和总升力系数可以分别表示为

(3)

式中：CN和CA分别是飞行器本体坐标系下的法向力气动系数和轴向力气动系数。

根据球面三角形公式可推出总攻角η与攻角α和侧滑角β的关系为

cosη=cosαcosβ

(4)

由于气动力在半速度坐标系中的分解为

(5)

式中：D为阻力；L为升力。

因此可以得到半速度坐标系中三轴的气动系数分别为

(6)

CL和CD是马赫数M和攻角α的函数，根据气动参数表进行二元拟合得到，具体拟合式为

(7)

实际任务中，再入飞行器的目标点一般不在初始轨道面内，此时需要飞行器利用气动力做出横向机动，横向机动的能力代表飞行器的有效覆盖范围，性能指标为[11]

minJ=-|sinδ|

(8)

式中，δ表示飞行器再入轨迹的横程对应的地心角，δ和横程Rcr可分别由式(9)、式(10)求出。

sinδ=-sinσ1sinφ2cosφ1cosΔλ-cosσ1cosφ1sinΔλ+

sinσ1cosφ2sinφ1

(9)

Rcr=RE|δ|

(10)

式中：σ1、φ1是再入点的飞行器弹道偏角和纬度；φ2是终点的纬度；Δλ是终点经度与起始点经度之差；RE为地球平均半径。

本文讨论的是横程最大问题，因此|δ|越大，则对应的横程越大。由于再入飞行器的机动能力是有限的，有|δ|<90°。在该条件下，|δ|越大，则|sinδ|越大，因此可以取公式(8)作为性能指标。

(11)

(12)

(13)

(14)

式中：R为球头半径；ρs为自由流密度；vs为自由流速度；ρ0为海平面大气密度1.225 kg/m3；v0为第一宇宙速度7 900 m/s。

对于终端点不确定的横向飞行距离最大问题，终端约束为终端高度、速度和弹道倾角等，可表示为

(15)

根据上述分析，轴对称再入飞行器横向飞行距离最大的优化问题P0可表示为

minJ=-|sinδ| s.t. (1),(11)～(15)成立

2 高斯伪谱法原理分析

采用高斯伪谱法对上述优化模型进行分析求解。高斯伪谱法作为直接法中的一种，是基于配点法发展而来的[11]，其核心思想是采用插值多项式基函数的实根作为配点，同时对控制变量和状态变量进行离散，利用离散点上的轨迹状态值构造拉格朗日多项式，对状态和控制进行全局的拟合，将连续空间的最优控制问题转化为非线性规划问题，从而方便使用数值方法求解。

2.1 时域变换

采用勒让德-高斯-洛巴托(Legendre-Gauss-Lobatto)点，简称LGL点，对连续的优化问题进行离散化，离散区间为[-1,1]。原动力学中状态变量和控制变量的时域为[t0,tf]，需要进行时域变换，将原优化问题的时间区间从t∈[t0,tf]转换到τ∈[-1,1]，建立仿射关系

(17)

可将动力学方程(1)表示为简化形式

则原动力学方程通过时域变换后的形式为

(18)

2.2 插值计算

定义集合Γ为N+1个LGL点，且集合中的 LGL点严格单调递增，以LGL点为插值节点，分别构建N+1个拉格朗日插值多项式对状态变量和控制变量进行插值拟合，如式(19)～(20)所示。

(19)

(20)

式中：x(τi)和u(τi)分别为τi点对应的真实状态变量和控制变量；Li(τ)为N次拉格朗日插值基函数，表示为

(21)

式中：PN为N次勒让德多项式。

显然Li(τk)具有如下性质

(22)

代入式(19)～(20)中，可得到在LGL插值点上估计值是精确等于真实值的，即

X(τi)=x(τi)

(23)

U(τi)=u(τi)

(24)

拉格朗日插值拟合的状态方程(19)的微分形式为

(25)

(26)

式中，k=0,1,2,…，N。

根据式(18)和式(25)构建微分偏差方程

(27)

原优化问题中的动力学约束可转为ξk=0的约束，即只与LGL节点处状态变量和控制变量有关的约束。同理原优化问题中的过程约束(11)～(14)也可按照上述方法进行插值并转化为与LGL节点处状态变量和控制变量有关的约束[12]，定义为

gk(x(τk),u(τk),τk)≤0

(28)

根据式(6)～(7)可看出，性能指标只与初始状态变量和末端状态变量有关，已为离散状态。

至此原优化问题P0转换为非线性规划问题P1，优化变量为每个离散点的状态变量x(τi)、控制变量u(τi)，可以表示为

(29)

2.3 解析雅可比矩阵

求解上述非线性规划问题时，需要多次对性能指标和约束雅可比矩阵进行求解，因此提高雅可比计算精度可以有效地提升非线性规划问题的收敛性[13]，在此给出解析雅可比矩阵的推导过程。

由式(9)可以看出，性能指标只与始末状态的经纬度和弹道偏角有关，因此性能指标对优化变量的雅可比矩阵可表示为

(30)

根据式(27)，ξk与每个LGL点的优化变量都有关，以第k点的约束为例，雅可比矩阵可表示为

(31)

式中：Fk为第k点的状态和控制变量。

对于式(28)中的过程约束，由原约束式(11)～(14)可以看出，各个约束量离散后都只与当前离散点的优化变量有关，约束式(15)只与最后一个离散点有关，同理可以求出解析雅可比矩阵，此处不再列举具体形式。

3 仿真与分析

表1 初始状态Tab. 1 Initial state

表2 末端状态约束Tab. 2 Terminal state restriction

图1 高度随时间变化曲线Fig.1 Time history of height

图2 速度随时间变化曲线Fig.2 Time history of velocity

图3 弹道倾角随时间变化曲线Fig.3 Time history of trajectory oblique angle

图1～3分别为飞行器再入过程中高度r、速度v以及弹道倾角θ随时间变化的仿真曲线。由图1可以看出轴对称飞行器再入过程呈现出跳跃式轨迹；结合图2可以发现再入初始阶段飞行器速度较快，飞行器高度变化较为剧烈；随着速度的降低，根据图3可以看出飞行器进入相对平稳滑翔阶段，弹道倾角保持在[-1°,1°]之间波动，飞行的最后阶段由于要满足末端点的状态约束，飞行器在滑行飞行末段进行了姿态调整，弹道倾角波动范围增大，末端点处为-2°，满足所给出的约束。

3种不同初始状态下的仿真结果对比如表3所示，可以看出再入时间与再入速度呈正相关，即再入速度越大，满足约束条件的再入时间越长，同时优化性能指标横向距离也越大。

表3 仿真结果Tab.3 Results of simulation

图4～5为控制变量随时间变化曲线，初始阶段弹道倾角为负数，因此采用大攻角再入，快速调整弹道倾角使热流密度峰值满足约束要求，之后滑翔阶段攻角缓慢增加，符合约束要求并且利于工程实现。初始阶段以较大侧滑角再入可以快速改变弹道偏角，使飞行器速度转向最大横向距离方向，之后侧滑角逐渐减小到0°，利用气动力进行速度方向调整，使得飞行器再入的横向距离增大。图6～8分别给出了热流密度、过载和动压3个主要路径约束随时间的变化曲线。

图4 攻角随时间变化曲线Fig.4 Time history of attack angle

图5 侧滑角随时间变化曲线Fig.5 Time history of sideslip angle

图6 热流密度随时间变化Fig.6 Time history of heat flux

图7 过载随时间变化Fig.7 Time history of overload

图8 动压随时间变化Fig.8 Time history of kinetic pressure

可以看出，3种再入状态下的路径约束都满足了设定的约束条件，说明高斯伪谱法均能得到满足约束并且使横向距离最大的最优控制曲线。

3个路径约束只有热流密度达到了约束边界，说明影响飞行过程的最主要因素为热流密度峰值约束。由热流密度约束式(11)可以看出热流密度主要受大气密度和飞行器速度影响，初始阶段飞行器高度快速下降导致大气密度增加，而飞行器速度下降并不明显，导致热流密度在此处达到峰值，因此飞行器初始阶段飞行条件较为恶劣。此阶段对飞行器控制系统、隔热系统等有较高要求，这也符合工程应用实际，从而进一步说明此方法得到的最优解接近实际情况。

综上分析，对于给定的3种初始状态，高斯伪谱法均能得到满足约束并且使横向距离最大的最优控制曲线，说明了算法的有效性和鲁棒性。

4 结束语

本文以天基再入飞行器轨迹优化设计为背景，对轴对称高马赫数飞行器再入问题进行建模分析，提出用高斯伪谱法解决此类问题并对算法进行详细说明。对不同的再入状态进行仿真计算，得到了较为平滑的控制曲线，并且满足各种约束条件，实现了轴对称高马赫数飞行器再入轨迹优化问题的求解。因此，本方法对于高马赫数飞行器再入轨迹的设计和优化有一定的工程意义。