基于伪谱法的再入可达域影响因素分析

李兆亭， 周 祥， 张洪波， 汤国建

(国防科技大学 空天科学学院，长沙 410073)

再入可达域是指可重复使用飞行器在满足约束条件下，到达着陆或者交班状态下所能允许的区域范围.它是飞行器纵向机动能力与横向机动能力的综合反映，体现了飞行器的整体机动性能.可达域求解可为飞行器轨迹规划与制导、应急着陆点选择、制导律性能评估提供依据，因此具有重要意义.

再入可达域的求解方法按照原理可以分为两类：参数化方法和非参数化方法.参数化方法是指将再入轨迹规划问题的状态量或控制量等进行参数化，采用某些优化方法对参数化后的问题进行优化求解；非参数化方法则不包含此类参数化、问题优化过程，常见的基于剖面的可达域求解方法就属于此类.如曾夕娟等[1]在“r-V”空间内引入约束的对数表达式简化再入走廊的表达，将再入轨迹分为3段设计，引入准平衡滑翔条件获取解析制导律以跟踪所规划轨迹，通过平移“边界跟踪段”轨迹使其布满再入走廊，从而得到最终的可达域.常松涛等[2]在基于E剖面的再入走廊下，引入了倾侧角反向策略，改进了阻力加速度曲线插值方式和最大可行阻力加速度曲线的生成方法，使计算得到的可达区域更精确、边界更均匀、弹道更平滑.He等[3]推导得到一个具有不确定性效应的再入走廊，其中攻角剖面可调节，纵向阻力剖面设计为上、下拟合安全边界的插值结果，横向升力-阻力剖面则利用准平衡滑翔条件得到，最后通过比例积分微分(PID)跟踪器来跟踪剖面，得到可达域.Zhang等[4]提出一种基于几何预测轨迹的可达域计算方法，在D-E走廊下，对横向运动状态、转弯半径和角度进行分析求解，得到D-E剖面的最远可达边界，并提出几何预测的轨迹；通过重复计算所有D-E剖面的可达边界，生成可达域.上述基于剖面的可达域求解方法通常需要人工设计相关剖面，而剖面参数决定了其可达域范围，因此得到的可达域通常只是实际可达域的一部分.此外，Liu等[5]提出一种使用多项式来描述可达域边界的方法，使用反向传播神经网络快速确定空对地飞行器的可达域.Zhang等[6]通过快速解析的方法求解3个不同的最优控制问题，得到可达域的内、外边界.章吉力等[7]对考虑禁飞区条件下的空天飞机再入可达区域问题进行研究，并基于极限绕飞轨迹提出一种不限禁飞区位置的可达区域求解方法.吴楠等[8]基于平衡滑翔假设和最优化飞行准则，通过数值积分获得最大纵程和横程弹道，基于椭圆近似法利用3个末端点构建可达域椭圆边界.

可达域的参数化求解方法包含的种类较多，常见方法主要包括基于虚拟目标点[9-11]、基于伪谱法[12-14]和基于凸优化方法[15]等.如赵吉松等[16]提出一种基于稀疏差分法和网格细化技术的快速、高精度轨迹求解方法，并利用可达域的求解验证了方法的有效性.蔺君等[17]对攻角、倾侧角进行参数化，利用带约束的差分进化算法求解满足再入过程约束和终端约束的再入轨迹.梁巨平等[18]选取倾侧角为时间的分段常值函数，采用遗传算法来求解一系列的轨迹优化问题.赵江等[19]设计了一种参数化的倾侧角剖面，利用约束粒子群优化算法求解满足再入过程约束和末端约束的最优滑翔轨迹.相比基于剖面的求解方法，可达域的参数化求解方法不需要设计相关剖面，其可达域范围更大且更接近实际值，但受限于具体算法的性能，其全局最优性仍难以保证.孙勇等[20]基于拟能量概念提出了多段轨迹优化方法，将轨迹按照拟能量进行分段，对各段分别进行优化，减少了每次优化的变量个数，计算效率得到较大提高.

伪谱法具有良好的适用性和强大的求解能力，基于伪谱法求解可达域的方案较多.如王涛等[21]提出一种基于Gauss伪谱法的再入可达域计算方法，采用固定的攻角剖面，仅对倾侧角进行单变量寻优.李柯等[22]基于Radau伪谱法求解得到飞行器的可达域，着重分析了升阻比、终端速度以及终端倾角对可达域的影响.Wang等[23]在准平衡滑行条件下，使用高斯伪谱法计算得到了返回航天器的可达域.解永锋等[24]通过定义加权的横程、纵程组合性能指标函数，将可达域求解问题转化为组合性能指标最优的控制问题，采用勒让德伪谱法快速计算可达域.

本文继承发展了文献[21]的求解思路，将攻角、倾侧角同时作为控制量代入到优化问题中，采用伪谱法进行轨迹规划，得到的可达域范围更广且更接近实际值.此外，描述了基于伪谱法的可达域快速求解方法，并对伪谱法进行了简要阐述；然后基于上述方法对影响可达域的相关因素进行仿真研究.仿真结果表明，飞行器质量、气动参考面积、大气密度等会对短纵程轨迹产生明显影响；同时，升阻比大小与可达域范围成正相关.

1 再入飞行器及其可达域描述

介绍再入飞行器的动力学模型，并在换极坐标系下给出利用伪谱法求解可达域的思路.

1.1 动力学方程及约束

在求解可达域时，通常将飞行器的运动转换到换极坐标系[25]中.在换极坐标系中，飞行器起点的经纬度均为0°，终点的经度和纬度分别描述了再入纵程和再入横程.利用这些特性，可方便地标定伪谱法中的参数范围，简化弹道规划算法.

在换极坐标系下，可以推导得到飞行器再入运动方程：

(1)

(2)

式中：ρ为大气密度；M为飞行器质量；S为参考面积；CL和CD分别为升力系数和阻力系数，均与攻角α有关.同时，再入飞行过程中还需考虑多种过程约束，如热流密度约束、过载约束、动压约束等，具体为

(3)

在本文中，目标飞行器为CAV-H，其相关参数和气动数据可参考文献[26].

1.2 换极坐标系下的可达域求解

为了加快优化速度和提高优化质量，文献[21]提出了一种基于Gauss伪谱法的再入可达域计算方法，采用固定的攻角剖面，将攻角作为状态变量，仅对倾侧角进行单变量寻优.本文延续该求解思路，但与之不同的是，攻角剖面并不是固定的，即攻角、倾侧角将同时作为待优化变量.

再入可达域问题可表述为求解不同纵程条件下的最大横程.在换极坐标系内，落点的横程可用纬度表征:

J=min(cosφpf)

(4)

式中：φpf为换极坐标系中终点的纬度，如图1所示.

图1中，Lset为设定的纵程，λpf为换极坐标系中终点的经度，二者数值相等.不妨将Lset视为终端约束，通过改变Lset值得到不同纵程条件下的最大横程轨迹.同时将换极坐标系下，φpf>0°的所有轨迹称之为北半球轨迹，反之称之为南半球轨迹.在进行伪谱法的参数设置时，考虑到地球自转的影响，即微分方程(1)中的科氏加速度项和牵连加速度项均不为0，不妨将把北半球轨迹的纬度变化范围设置为 [-2°, 90°]，将南半球轨迹的纬度变化范围设置为 [-90°, 2°].

图1 不同纵程下的最大横程Fig.1 Maximum cross-ranges of different longitudinal ranges

求解上述不同纵程下的最大横程问题，即可得到有效的可达域.但实际仿真中，同样需要确定最大、最小纵程(Lmax和Lmin)，以便确定Lset的取值范围.因此，最终的规划流程如图2所示.首先求解最大、最小纵程，标定Lmax和Lmin的具体数值；然后从纵程Lmax-ΔL开始，分别求解两个最大横程问题，依次减小纵程ΔL，直到达到Lmin；最后即可得到飞行器的再入可达域.

图2 数值求解可达域的一般流程Fig.2 General procedure for numerically solving footprint

2 伪谱轨迹规划方法

介绍伪谱法的原理，描述将原问题转化为非线性规划(NLP)问题的过程，给出伪谱法求解问题的一般流程.

2.1 伪谱离散化方法

伪谱法作为直接法中的典型代表，其原理与直接法一致：基于离散化一组节点上的控制变量或状态变量，将最优控制问题转换为NLP问题.不同的是，伪谱法使用全局多项式同时对状态量和控制量进行参数化，并通过正交配点来近似微分代数方程.其一般流程如下所示.

首先引入新的时间变量τ，对原问题做如下转换：

(5)

式中：t0为初始时间；tf为终端时间.转换后得到的状态量、控制量均采用拉格朗日多项式进行逼近.

(6)

式中：τ∈[-1, +1]；x为实际的状态量；X为多项式拟合的状态量；lj(τ),j=1, 2, …,Nk+1为拉格朗日多项式；τ1,τ2, …,τNk为配点，通常可以选择Legendre-Gauss(LG)、Legendre-Gauss-Radau(LGR)、Legendre-Gauss-Lobatto(LGL)等正交配点[27].对于拉格朗日多项式，显然有：

(7)

对拉格朗日多项式进行微分，可以得到微分矩阵算子：

(8)

DX-f(X)=0

(9)

伪谱法中对于动力学微分方程的处理如式(9)所示，其作为等式约束与其他约束一起组成最终的NLP问题.可以看到，由于引入了微分方程的右侧项f(X)，此约束一般而言是非线性约束，所以最终形成的离散问题，也就是NLP问题.关于再入过程中的其他约束，如式(3)所示，可以通过状态量、控制量的拟合式(6)进行显示表达.

最终，原再入轨迹规划问题就被转化为NLP问题.实际上，考虑到轨迹规划问题的复杂性，求解时不仅需要增加多项式的阶数，通常也需要增加多项式的个数，即为常用的ph网格细化方法，将归一化后的总时间[-1, 1]划分为多个区间，每一个区间对应一个拉格朗日多项式.

2.2 NLP问题及其求解流程

利用上述伪谱离散化方法，可以将原问题转化为NLP问题.通常可以将NLP问题描述为

minΦ(Z)

(10)

s.t.Fmin≤F(Z)≤Fmax

Zmin≤Z≤Zmax

式中：Φ为性能指标函数；Z为自变量，通常包括状态量、控制量、时间变量、参数变量等；Zmax和Zmin分别为上述各变量的上、下限；F为等式和不等式约束，通常包含关于动力学微分方程的等式约束、动压约束、过载约束、热流约束等；Fmax和Fmin分别对应约束的上、下限，特别地，等式约束对应的上、下限均为0.

得到如式(10)所示的NLP问题后，可以采用通用的非线性规划求解器对问题进行求解，如稀疏非线性优化器(SNOPT)[28]或内点优化器(IPOPT)[29]等.为使NLP问题能够被求解，需要即时计算一些关键信息，如梯度、雅可比矩阵、海森矩阵等.此类信息的解算速度严重影响了最终的计算时间，为此，文献[30]推导出了NLP问题的梯度、雅可比矩阵和拉格朗日海森矩阵的明确表达式，确定了NLP问题的稀疏结构，可以显著提高NLP问题的计算效率和可靠性.

在NLP问题被求解后，通常来说此时的解很难满足期望的精度要求.因此，需要对网格进行细化，在更细密的网格下离散原问题，并求解得到更大规模的NLP问题.如此，上述网格被不断细化，直到满足迭代截止条件.

综上所述，伪谱法求解优化问题的流程如图3所示.图中，原问题首先在初始网格条件下进行离散化得到NLP问题，然后调用相关求解器求解此类问题，最后不断加密网格并求解离散化得到的NLP问题，直到满足截止条件.需要注意的是，由于稀疏网格条件下NLP问题与原问题不等价，NLP问题的收敛容限可以适当增加以减少计算时间、提高求解效率.

图3 伪谱法求解优化问题的一般流程Fig.3 General flow of pseudospectral method for solving optimization problems

3 仿真分析

将上述伪谱法应用于图3所示的可达域求解方案中，即可快速求解再入飞行器的可达域.

3.1 标准情况下的可达域求解

飞行器在初始时刻的经度、纬度被设定为(0°，0°)，仿真采用的大气环境为指数型大气模型，其他参数如表1所示.仿真计算机硬件条件为AMD Ryzen 7 5800H@3.2 GHz，伪谱法求解器SWIFT由本小组独立开发.

表1中给出了伪谱法中一些参数的设置范围.其中，h为飞行器的滑翔高度，φS为求解换极坐标系下南半球的飞行轨迹时的纬度，φN为求解换极坐标系下北半球的飞行轨迹时的纬度.为保证问题的全局最优性，部分参数范围取得较宽泛，如速度方位角、速度倾角等.

表1 伪谱法求解可达域时的部分参数设置Tab.1 Parameters for solving footprint by pseudospectral method

此外，一些其他参数设置如下：初始网格设置为2×2，即将总时间分为2段，每一段的轨迹采用2阶多项式进行逼近；最大迭代次数设置为5；原问题期望的收敛容限设置为10-4；最大飞行时间设置为 4 500 s；过载不超过3 g，动压不超过100 kPa，热流不超过150 kW/m2.最终得到的可达域仿真结果如图4所示.

图4(a)～(c)中， 红色、绿色、蓝色曲线分别代表最大纵程轨迹、北半球最小纵程轨迹、南半球最小纵程轨迹，其中最大纵程约为136.4°，南北半球最小纵程约为3.62°；32条黑色曲线代表不同纵程条件下的最值横程轨迹，其中所有轨迹的最值横程约为±54.52°.将所有轨迹的终点连接在一起，即可形成最终的飞行器可达域.可见，标准情况下此类飞行器的可达域接近于椭圆形.同时给出文献[21]中带有攻角剖面的可达域计算结果，如图4(d)中黄色实线所示，可见受攻角剖面的限制，这种情况下得到的可达域相对较小.此时，最大纵程约为112.6°，相比136.4°，损失约17.45%.最值横程约为±41.53°，相比±54.52°，损失23.83%.

整体而言，对于构成可达域的图4(a)～(c)中的35条轨迹，其速度、高度、速度方位角等随时间变化较为规律.对于速度而言，纵程越大，速度减小得越慢.同时从4(b)也可以看到，速度减小变慢的原因在于飞行器的滑翔高度较高；从4(c)可知，整个可达域条件下的速度方位角在-90°～270°之间变化.

图4 各条轨迹的状态量变化Fig.4 State change for each trajectory

给出上述35条规划轨迹的控制量变化，如图5所示.其中，攻角的变化总是以一个较大的值开始，逐步过渡到最大升阻比攻角附近，最后阶段受表1中终端速度倾角范围设置的影响，攻角快速拉升到30°附近，以抬高速度倾角至设定范围.倾侧角的变化则明显分为南、北半球两部分.当规划北半球轨迹时，倾侧角多数情况下小于0°；南半球轨迹则相反.同时可以看到，纵程越大的轨迹，起点处的倾侧角绝对值越小，相对时间的平均斜率的绝对值也越小；而纵程越小的轨迹，起点处的倾侧角绝对值越大，相对时间的平均斜率的绝对值也越大.

图5 各条轨迹的控制量变化Fig.5 Control changes for each trajectory

3.2 参数拉偏条件下的可达域

研究某些参数对可达域的影响，将某些参数进行拉偏，在其他参数不变的条件下求解可达域.一些拉偏情形如表2所示；仿真得到的可达域随拉偏参数的变化如图6所示.

表2 5种参数拉偏情形Tab.2 Five cases of parameter deviation

图6(a)中，在质量拉偏±5%之间，飞行器的可达域基本没有变化；而当质量拉偏+10%时，短纵程轨迹的某些约束达到上限，其最小纵程有了明显变化，但最大纵程基本没有变化.同样的现象出现在图6(b)和6(c)中，气动参考面积和大气密度在±5%之间波动时，可达域不发生变化；当拉偏在 -8% 时，其最小纵程有了明显的损失，而最大纵程不变.

图6 可达域随参数拉偏情况的变化Fig.6 Variation of footprint with parameter deviation

结合式(1)～(3)可知，飞行器的质量、气动参考面积、大气密度仅出现在式(2)和式(3)中，并不显含在动力学微分方程(1)中.并且质量、面积、大气密度的变化仅会导致升力、阻力的同步变化，因此对于飞行器某一时刻的飞行状态而言，通过适当改变攻角、倾侧角，就可以使得式(1)中的左侧微分项完全一致，则在一定范围内，可达域是不变的.而超出此范围后，受到热流、过载、动压等约束的影响，部分轨迹，特别是短纵程轨迹，可能存在达到约束上界的情况，因此必然要通过增加纵程以降低过程约束.

此外，图6(d)为升力和阻力系数拉偏条件下的可达域计算结果.其中，实线代表升力系数拉偏的可达域，虚线代表阻力系数拉偏的可达域.可见，升力系数的增加和阻力系数的减小均会导致可达域范围增加，反之亦然.同时，升力系数拉偏+10%和阻力系数拉偏-10% 的结果一致，升力系数拉偏-10% 和阻力系数拉偏 +10% 的结果也一致，说明升力、阻力系数的变化对可达域的影响是等价的.更准确地，升阻比的大小影响了可达域的大小：升阻比越大，可达域越大；升阻比越小，可达域越小.

4 结语

本文描述了基于伪谱法的可达域快速求解方法，并基于此方法对影响可达域的相关因素进行了研究与仿真.仿真结果表明，在给定参数条件下，飞行器质量、气动参考面积、大气密度等在±5%范围内不会导致可达域的改变；超出一定范围后，如飞行器质量增加10%、气动参考面积减少8%、大气密度减少8%等，均会对短纵程轨迹产生明显影响，且体现在可达域的左半部分，其右半部分不受影响；升阻比对可达域的影响较大，其大小与可达域范围成正相关.