自适应分段的多段线性伪谱广义标控脱靶量再入制导方法
阅读：1037发布：2020-07-17
专利汇可以提供自适应分段的多段线性伪谱广义标控脱靶量再入制导方法专利检索，专利查询，专利分析的服务。并且本 发明 公开了自适应分段的多段线性伪谱广义标控脱靶量再入制导方法，该制导方法不仅适用于高 升阻比 再入 飞行器 ，且能满足多种终端约束和过程约束。该方法结合了非线性模型预测控制，多段线性伪谱法以及变分法的基本思想，与Yu方法最大的不同在于，处理的对象为考虑耦合效果的降阶非线性再入动 力 学，并且除了考虑终端 位置 要求和 倾侧 角 要求外，还考虑了终端高度和弹道倾角约束。，下面是自适应分段的多段线性伪谱广义标控脱靶量再入制导方法专利的具体信息内容。
1.自适应分段的多段线性伪谱广义标控脱靶量再入制导方法，其特征在于，具体步骤如下：
(1)初始化：设置初始、终端的计算仿真参数，通过离线弹道优化等获得合理的初始控制参数和分段时间初值；
(2)下落段制导：以最大攻角和零倾侧角下落，直到高度的变化率为零，进入滑翔段阶段；
(3)滑翔段制导倾侧反转判断：根据当前能量状态与分段时间的关系，判断是否进行倾侧反转，当当前能量大于分段时刻的能量时，不进行反转，进入步骤(4)，当当前能量小于分段时刻的能量时，且不为最后一个反转点时，重新建立非线性参数控制问题，进入步骤(4)，当为最后一个反转点时，进入步骤(9)；
(4)滑翔段制导预测弹道积分：使用当前时刻的状态量作为积分初值，所述的初始控制参数和分段时间初值作为控制输入，其中当前时刻所对应的控制参数记为U0，进行预测弹道积分，获得广义标控脱靶量即终端状态偏差dψ以及全局多段弹道信息，Xk，Uk；
(5)精度判定：判断所述广义标控脱靶量的精度范围，如果所述dψ满足所述步骤(1)设置的终端计算仿真参数的精度要求，进入步骤(6)，如果不满足精度要求，进入步骤(8)；
(6)过程边界控制：判断所述的U0是否满足多种过程约束界限，当所述的U0超过边界控制Umax时，所述的U0等于Umax，否则，所述的U0等于U0，进入步骤(7)；
(7)执行控制指令：将所述的U0作用到弹道阻尼控制方程中，获得平稳滑翔的控制指令Uc，并将其作用到实际的系统中去，返回步骤(3)，并且将当前的分段时间作为倾侧反转的判断基准，当前的U0和分段时间作为下一步弹道积分的控制输入；
(8)滑翔段制导控制参数和分段时间更新：分别对预测的每段积分弹道进行线性化处理，结合多段伪谱法和变分法原理，获得满足终端标控脱靶量修正的控制参数和分段时间解析修正关系，进而获得控制参数和分段时间的更新，并进入步骤(6)；
(9)末制导段显示制导：分别根据当前状态和终端约束条件在纵、横平面采用比例导引和多项式制导，获得能够导引飞行器安全地飞向目标且满足多种终端约束的制导指令。
2.如权利要求1所述的再入制导方法，其特征在于，所述步骤(8)具体包括：
首先，构建降阶的再入动力学方程
其中，所述e是相对能量，所述E是绝对能量；所述 是经度、纬度和航向角；所述F是总升阻比函数，是以相对能量为自变量的函数且在攻角确定情况下能够通过平稳滑翔条件离线获得；所述u为控制变量表示纵向升阻比，是总升阻比在垂直方向上的一个分量；所述kld是由惯组输出的升阻比偏差估计，并且在预测积分过程中是常值；所述μ是引力常数，所述r为地心距；
然后，将控制指令参数化；将纵向升阻比曲线设计成分为三段的参数函数，其具体表达形式如公式(2)所示：
其中，所述E0和Ef分别代表初始和终端绝对能量，所述ef代表终端相对能量，所述Ere1和Ere2分别表示第一次倾侧反转和第二次倾侧反转发生时的绝对能量，所述k1为常数，所述c1和c2是二次曲线的参数，其取值通过分段函数光滑，导数连续两个条件确定；
其次，将多段非线性方程线性化，将具有终端约束的降阶再入运动学方程如公式(3)所示
其中，所述 当初始的k1，Ere1和Ere2给定，通过欧拉法或龙格库塔数值积分方法预测全局的状态量，获得预测的终端状态与实际所需的终端状态之间的终端误差表示为δxf＝x(Ef)-xf，将所述公式(3)在预测弹道周围进行高阶泰勒展开，获得一组以状态偏差为自变量的误差传播动力学方程如公式(4)所示；
其中，所述δx是状态量的偏差，所述δu是控制量的偏差，为标量；所述A1、A2、A3是系数矩阵，为3×3的矩阵，所述B1、B2、B3是系数矩阵，为3×1的向量；
最后，使用多段线性伪谱法求解控制量的偏差δu，修正控制序列：
先推导一般的线性伪谱校正公式，将线性动力学方程在LG节点上离散为如公式(5)所示：
其中，所述D是微分逼近矩阵，下标k、l分别代表第k个与第l个插值点；
以所述公式(5)中的第一段代数约束关系式为例，第一段的终端状态偏差δxre1可以通过Gauss型积分公式而表示为初始状态偏差和在LG节点上的状态偏差的显示函数关系如公式(6)所示，
其中，ω是Gauss型积分的权函数，所述δxre1为关于初始状态偏差δx0和控制偏差δu的显式解析表达式，其表达形式如公式(7)所示，
其中，所述 是系数矩阵，为3×3的矩阵，所述 是3×1的向量，同理可以获得系数矩阵 和 和 根据链式法则建立起关于初始状态偏差δx0、控制偏差δu以及终端状态偏差δxf之间的函数关系，其具体表达式如公式(8)所示，
直接获得关于所述控制偏差δu与所述终端状态偏差δxf之间的修正关系如公式(9)所示，
用来修正最近的倾侧反转点从而消除积累的纬度误差的校正方法如下：
终端经度偏差 与终端纬度偏差 的关系如公式(10)所示，
其中，所述 与 分别表示终端经度约束与终端纬度约束；所述δEk表示第k次指令修正过程中的第一次反转点的绝对能量偏差； 与 表示终端状态量的函数
关系
于是，由反转点改变δEk所引起的原反转点状态的等时变分 表示为如公式(11)所示，
进一步的，由反转点改变δEk所引起的终端状态偏差表示为如公式(12)所示，
所以，作为控制量的纵向升阻比的修正δu和倾侧反转能量时刻的修正δEk表达为如公式(13)所示，
3.如权利要求1或2所述的再入制导方法，其特征在于，所述步骤(9)具体包括：
设计纵、横向制导平面的指令分配算法，通过调整攻角和倾侧角实现两个制导平面内的制导指令，其中所述横向制导平面采用比例导引如公式(14)所示，
其中，所述n1是横向制导指令，所述Npc是比例导引常数，所述 是视线角速度，所述v是投影在所述横向平面的速度；
所述纵向制导平面采用多项式制导方法，制导指令u0为
其中，所述sf，hf和γf分别为需要在TAEM交班点满足的剩余飞行器距离，高度和弹道倾角，并且所述s0，h0和γ0分别为飞行器当前的相应状态量，所述m和n为导引增益；
当保留地球曲率并且忽略地球自转的再入动力学方程，纵向加速度nvertical和横向加速度 分别的表示为：
其中，所述γ、ψ、φ分别表示当前弹道倾角、航向角和纬度；
由于侧向加速度和纵向加速度都是总加速度的分量，获得指令倾侧角的表达式如公式(17)所示：
之后，采用牛顿迭代确定所需的指令攻角如公式(18)所示：
其中，Cl为升力系数，Clα为升力系数对攻角的导数；k表示第k次迭代过程。
自适应分段的多段线性伪谱广义标控脱靶量再入制导方法
技术领域
[0001] 本发明涉及航空飞行领域，具体涉及一种自适应分段的多段线性伪谱广义标控脱靶量再入制导方法。
背景技术
[0002] 众所周知，再入飞行器在再入飞行阶段需要经受严酷的再入环境，而再入制导系统即——提供导引飞行器安全准确地飞向目标的制导指令，对于再入飞行至关重要。随着阿波罗飞船和航天飞机的成功飞行，再入制导算法的研究取得了巨大的进步。
[0003] 尽管现在存在很多类型的制导算法，最著名的制导算法当属航天飞机制导律，这也是所有跟踪制导算法的基础。航天飞机再入飞机制导律主要分为两部分，一部分为在安全再入走廊之内，规划一个满足射程要求的阻力加速度剖面；另一部分为基于简化的阻力动力学线化模型，设计剖面跟踪制导算法。调整飞行器的倾侧角幅值用于消除阻力加速度偏差，而通过预设的门限走廊调整倾侧角的符号用以控制航向角误差。为了增强制导算法的制导性能，Mease和Kremer在反馈线性化的框架下，重新研究了航天飞机再入制导算法，并指出如果制导模型建立足够精确，动力学中的非线性项可以得到补偿。Bharadwaj和Rao采用逼近反馈线性化理论，设计再入跟踪制导律，其中，纵、横程同时作为跟踪输出参数。Mease同样也研究了基于加速度规划的降阶再入弹道规划方法，用于提供考虑纵、横向运动的阻力加速度规划。Dukeman采用线性二次调节理设计纵向剖面的跟踪方法，仿真结果表明，LQR算法对于不同的任务具有很强的鲁棒性。
[0004] 尽管参考弹道跟踪算法取得了巨大的成功，但是该类算法太依赖于离线生成的参考弹道，为了克服这个弊端，增大飞行器应对更大扰动偏差的能力，无数的学者将研究重点转向预测校正再入制导律上，该算法是通过在线迭代更新控制参数，用以消除预测的模型偏差。Zimmerman提出了一种智能规划方法，用于在线计算满足热流需求的可行再入弹道。Shen和Lu基于拟平衡滑翔假设，发展了快速在线规划算法用于生成满足路径约束的可行弹道，但是，现有的预测校正算法都是采用传统的航向误差角走廊，进行横向弹道控制，对于中、高升阻比的再入飞行器，该方法将导致较大的横向误差，因此，Shen研究了采用终端横向距离，作为走廊的横向制导控制策略，提高了中等升阻比飞行器的横向制导精度。
[0005] 需要格外指出的是，这些预测校正算法的时间消耗，将随控制参数个数的增加而呈指数倍增加，这是由于用于更新控制参数的梯度信息是通过采用数值积分的割线法获得，因此，为了减少算法的复杂性，并增强算法的鲁棒性，这些预测校正算法通常只考虑了一二个控制参数，这样将只能保证较少的终端约束。此外，这些算法都是采用门限走廊限制策略来较少航向角误差，对于高升阻比飞行器，这样不仅将极大地限制飞行器的横向机动能力，还会带来较大的横向误差和较多的倾侧反转次数。同时，为了更好地完成中、末交班，需要考虑更多的终端约束条件。
发明内容
[0006] 本发明的目的是针对高升阻比滑翔飞行器开发一种能够充分发挥飞行器横向机动能力且能保证多种终端约束的再入制导方法，即提供一种自适应分段的多段线性伪谱广义标控脱靶量再入制导方法，该方法解决的关键技术问题是快速确定倾侧角反转时刻的多段在线再入弹道规划算法。
[0007] 为实现上述目的，本发明采用以下技术方案：
[0008] 自适应分段的多段线性伪谱广义标控脱靶量再入制导方法，具体步骤如下：
[0009] (1)初始化：设置初始、终端的计算仿真参数，通过离线弹道优化等获得合理的初始控制参数和分段时间初值；
[0010] (2)下落段制导：以最大攻角和零倾侧角下落，直到高度的变化率为零，进入滑翔段阶段；
[0011] (3)滑翔段制导倾侧反转判断：根据当前能量状态与分段时间的关系，判断是否进行倾侧反转，当当前能量大于分段时刻的能量时，不进行反转，进入步骤(4)，当当前能量小于分段时刻的能量时，且不为最后一个反转点时，重新建立非线性参数控制问题，进入步骤(4)，当为最后一个反转点时，进入步骤(9)；
[0012] (4)滑翔段制导预测弹道积分：使用当前时刻的状态量作为积分初值，所述的初始控制参数和分段时间初值作为控制输入，其中当前时刻所对应的控制参数记为U0，进行预测弹道积分，获得广义标控脱靶量即终端状态偏差dψ以及全局多段弹道信息，Xk，Uk；
[0013] (5)精度判定：判断所述广义标控脱靶量的精度范围，如果所述dψ满足所述步骤(1)设置的终端计算仿真参数的精度要求，进入步骤(6)，如果不满足精度要求，进入步骤(8)；
[0014] (6)过程边界控制：判断所述的U0是否满足多种过程约束界限，当所述的U0超过边界控制Umax时，所述的U0等于Umax，否则，所述的U0等于U0，进入步骤(7)；
[0015] (7)执行控制指令：将所述的U0作用到弹道阻尼控制方程中，获得平稳滑翔的控制指令Uc，并将其作用到实际的系统中去，返回步骤(3)，并且将当前的分段时间作为倾侧反转的判断基准，当前的U0和分段时间作为下一步弹道积分的控制输入；
[0016] (8)滑翔段制导控制参数和分段时间更新：分别对预测的每段积分弹道进行线性化处理，结合多段伪谱法和变分法原理，获得满足终端标控脱靶量修正的控制参数和分段时间解析修正关系，进而获得控制参数和分段时间的更新，并进入步骤(6)；
[0017] (9)末制导段显示制导：分别根据当前状态和终端约束条件在纵、横平面采用比例导引和多项式制导，获得能够导引飞行器安全地飞向目标且满足多种终端约束的制导指令。
[0018] 如上所述的再入制导方法，优选地，所述步骤(8)具体包括：
[0019] 首先，构建降阶的再入动力学方程
[0020]
[0021] 其中，所述e是相对能量，所述E是绝对能量；所述 是经度、纬度和航向角；所述F是总升阻比函数，是以相对能量为自变量的函数且在攻角确定情况下能够通过平稳滑翔条件离线获得；所述u为控制变量表示纵向升阻比，是总升阻比在垂直方向上的一个分量；所述kld是由惯组输出的升阻比偏差估计，并且在预测积分过程中是常值；所述μ是引力常数，所述r为地心距；
[0022] 然后，将控制指令参数化；将纵向升阻比曲线设计成分为三段的参数函数，其具体表达形式如公式(2)所示：
[0023]
[0024] 其中，所述E0和Ef分别代表初始和终端绝对能量，所述ef代表终端相对能量，所述Ere1和Ere2分别表示第一次倾侧反转和第二次倾侧反转发生时的绝对能量，所述k1为常数，所述c1和c2是二次曲线的参数，其取值通过分段函数光滑，导数连续两个条件确定；
[0025] 其次，将多段非线性方程线性化，将具有终端约束的降阶再入运动学方程如公式(3)所示
[0026]
[0027] 其中，所述 当初始的k1，Ere1和Ere2给定，通过欧拉法或龙格库塔数值积分方法预测全局的状态量，获得预测的终端状态与实际所需的终端状态之间的终端误差表示为δxf＝x(Ef)-xf，将所述公式(3)在预测弹道周围进行高阶泰勒展开，获得一组以状态偏差为自变量的误差传播动力学方程如公式(4)所示；
[0028]
[0029] 其中，所述δx是状态量的偏差，所述δu是控制量的偏差，为标量；所述A1、A2、A3是系数矩阵，为3×3的矩阵，所述B1、B2、B3是系数矩阵，为3×1的向量；
[0030] 最后，使用多段线性伪谱法求解控制量的偏差δu，修正控制序列：
[0031] 先推导一般的线性伪谱校正公式，将线性动力学方程在LG节点上离散为如公式(5)所示：
[0032]
[0033] 其中，所述D是微分逼近矩阵，下标k、l分别代表第k个与第l个插值点；
[0034] 以所述公式(5)中的第一段代数约束关系式为例，第一段的终端状态偏差δxre1可以通过Gauss型积分公式而表示为初始状态偏差和在LG节点上的状态偏差的显示函数关系如公式(6)所示，
[0035]
[0036] 其中，ω是Gauss型积分的权函数，所述δxre1为关于初始状态偏差δx0和控制偏差δu的显式解析表达式，其表达形式如公式(7)所示，
[0037]
[0038] 其中，所述 是系数矩阵，为3×3的矩阵，所述 是3×1的向量，同理可以获得系数矩阵 和 和 根据链式法则建立起关于初始状态偏差δx0、控制偏差δu以及终端状态偏差δxf之间的函数关系，其具体表达式如公式(8)所示，
[0039]
[0040] 直接获得关于所述控制偏差δu与所述终端状态偏差δxf之间的修正关系如公式(9)所示，
[0041]
[0042] 用来修正最近的倾侧反转点从而消除积累的纬度误差的校正方法如下：
[0043] 终端经度偏差 与终端纬度偏差 的关系如公式(10)所示，
[0044]
[0045] 其中，所述 与 分别表示终端经度约束与终端纬度约束；所述δEk表示第k次指令修正过程中的第一次反转点的绝对能量偏差； 与 表示终端状态量的函数关系
[0046] 于是，由反转点改变δEk所引起的原反转点状态的等时变分 表示为如公式(11)所示，
[0047]
[0048] 进一步的，由反转点改变δEk所引起的终端状态偏差表示为如公式(12)所示，[0049]
[0050] 所以，作为控制量的纵向升阻比的修正δu和倾侧反转能量时刻的修正δEk表达为如公式(13)所示，
[0051]
[0052] 如上所述的再入制导方法，优选地，所述步骤(9)具体包括：
[0053] 设计纵、横向制导平面的指令分配算法，通过调整攻角和倾侧角实现两个制导平面内的制导指令，其中所述横向制导平面采用比例导引如公式(14)所示，
[0054]
[0055] 其中，所述n1是横向制导指令，所述Npc是比例导引常数，所述 是视线角速度，所述v是投影在所述横向平面的速度；
[0056] 所述纵向制导平面采用多项式制导方法，制导指令u0为
[0057]
[0058] 其中，所述sf，hf和γf分别为需要在TAEM交班点满足的剩余飞行器距离，高度和弹道倾角，并且所述s0，h0和γ0分别为飞行器当前的相应状态量，所述m和n为导引增益；
[0059] 当保留地球曲率并且忽略地球自转的再入动力学方程，纵向加速度nvertical和横向加速度nlatral分别的表示为：
[0060]
[0061] 其中，所述γ、ψ、φ分别表示当前弹道倾角、航向角和纬度；
[0062] 由于侧向加速度和纵向加速度都是总加速度的分量，获得指令倾侧角的表达式如公式(17)所示：
[0063]
[0064] 之后，采用牛顿迭代确定所需的指令攻角如公式(18)所示：
[0065]
[0066] 其中，Cl为升力系数，Clα为升力系数对攻角的导数；k表示第k次迭代过程。
[0067] 本发明提供的自适应分段的多段线性伪谱广义标控脱靶量再入制导方法，该制导方法不仅适用于高升阻比再入飞行器，且能满足多种终端约束和过程约束。该方法结合了非线性模型预测控制，多段线性伪谱法以及变分法的基本思想，与Yu方法最大的不同在于，处理的对象为考虑耦合效果的降阶非线性再入动力学，并且除了考虑终端位置要求和倾侧角要求外，还考虑了终端高度和弹道倾角约束，能够处理多段非线性控制问题的线性伪谱在线规划，并且该方法通过规划控制变量将其表示为几个参数的形式，以至于不需要考虑协态变量的影响，此外，该方法还能处理状态和控制不连续的问题。
[0068] 本发明的方法在线性伪谱广义标控制导控制方法和变分原理的基础上，提出了一种新型的滑翔段预测-校正制导算法，该算法首先通过弹道积分获得全弹道信息，并在预测弹道周围进行多段非线性系统的线性化，之后通过变分原理将反转时刻能量的影响，转化为对反转时刻状态的影响，并结合广义标控脱靶量就能确定纵向升阻比和倾侧反转能量时刻的偏差，对终端状态偏差的多段误差传播方程，结合多段伪谱法离散该方程，就能获得关于终端广义标控脱靶量修正的解析表达式，因而，仅仅通过一次弹道积分预测就能同时提供，用以消除纵向和横向误差的纵向升阻比和倾侧反转时刻的更新，大大提高了计算效率。
[0069] 本发明方法分别设计了作用于球面地球假设下的比例导引和多项式导引，它们在最后一次反转后，分别在纵横向方向保证各自平面内的多种终端约束，此外，该方法的适用性，有效性和制导性能将通过高敏的数值算例进行验证，该方法的鲁棒性也通过考虑较大不确定性和随机干扰的蒙特卡洛打靶试验进行验证，仿真结果表明该算法不仅具有很快的计算效率，而且具有很高的计算精度，即使在恶劣的再入环境下也能提供稳定和鲁棒的再入制导性能。
[0070] 本发明方法还具有以下特点：
[0071] (1)能够基于广义标控脱靶量获得解析的标控修正关系，因而计算速度快，计算精度高；
[0072] (2)能在线确定反转时刻，从而充分发挥高升阻比再入飞行器的巨大横向机动能力；
[0074] 图1为本发明滑翔段伪谱预测校正制导算法和末段调整段的制导方法的实施流程图。
[0075] 图2为状态变化示意图。
[0076] 图3为蒙特卡洛打靶地面轨迹曲线。
[0077] 图4为蒙特卡洛打靶高度-速度曲线。
[0078] 图5为蒙特卡洛打靶攻角曲线。
[0079] 图6为蒙特卡洛打靶倾侧角曲线。
[0080] 图7为蒙特卡洛打靶航向角误差曲线。
[0081] 图8为蒙特卡洛打靶速度曲线。
[0082] 图9为蒙特卡洛打靶热流曲线。
[0083] 图10为蒙特卡洛打靶动压曲线。
[0084] 图11为蒙特卡洛打靶过载曲线。
[0085] 图12为蒙特卡洛打靶弹道倾角曲线。
[0086] 图13为蒙特卡洛打靶升阻比变化曲线。
[0087] 图14为蒙特卡洛打靶倾侧角-高度误差散布。
[0088] 图15为蒙特卡洛打靶航向角-速度误差散布。
[0089] 图16为蒙特卡洛打靶弹道倾角散布。
具体实施方式
[0090] 本发明的自适应分段的多段线性伪谱广义标控脱靶量再入制导方法，其中，关于滑翔段伪谱预测校正制导算法和末段调整段的制导方法的实施流程图，如图1所示。其中，滑翔段模块表示滑翔段的多段线性伪谱预测-校正再入制导算法的实施流程，末段模块表示末段调整段的实施流程，需要指出的是，在第一次倾侧反转前，预测的积分弹道和线性动力学方程由三部分组成，并且多段线性伪谱预测-校正再入制导算法根据预测误差修正第一个倾侧反转点和纵向升阻比的大小。当第一次倾侧反转之后，预测的弹道积分和线性动力学方程将缩减到只由二部分组成，且多段线性伪谱在线规划再入制导方法，根据预测误差修正第二个倾侧反转点和纵向升阻比的大小。当第二次倾侧反转发生时，飞行器将进入到末段调整段，通过比例导引和多项式制导律分别获得纵向和横向制导指令，并通过指令分配环节，将制导指令转化为攻角指令和倾侧角指令，并且最终导引飞行器飞向TAEM交班点，且满足多个终端约束。
[0091] 下面将结合附图和实例对本发明作进一步的详细说明。
[0092] 实施例1
[0093] 自适应分段的多段线性伪谱广义标控脱靶量再入制导方法，具体步骤如下：
[0094] (1)初始化：设置初始、终端的计算仿真参数，通过离线弹道优化等获得合理的初始控制参数和分段时间初值。
[0095] 即初始化建模，具体如下：
[0096] 1.)再入动力学方程
[0097] 在圆球、自转地球假设条件下，三自由度质点再入动力学方程可以表示为[0098]
[0099] 其中，所有的动力学方程都是对时间的导数，r表示为飞行器质心到地球圆心的地心距离，单位为m；θ和φ分别表示飞行器所在位置的经度和纬度，其单位为°；ν表示飞行器相对于地球的速度，其单位为m/s；γ表示为飞行器相对于地球速度矢量与当地水平面的夹角，称之为弹道倾角，ψ表示为飞行器相对于地球速度矢量在当地水平面的投影与正北方向的夹角，并以顺时针旋转为正，称之为航向角，它们的单位都为°；m表示飞行器的质量，其单位为kg；g＝μ/r2表示为飞行器所受的重力加速度，其中，μ是地球重力常数，重力加速度的单位为m/s2；σ表示飞行器沿速度方向旋转的角度，称之为倾侧角，其单位为°；ω表示地球-5自转角速度，其值等于7.2921×10 rad/s；L和D分别表示飞行器所承受的升力和阻力，其表达式表示为：
[0100]
[0101] 其中，ρ＝ρ0exp(-h/H)表示大气密度，ρ0表示为海平面标准大气压，其值等于1.225kg/m3，h表示飞行器质心到地球表面的高度，且地球的半径为6378245m；H表示大气密度常数，其值等于7200；Sref表示飞行器的特征面积，Cl和Cd分别表示飞行器的升力系数和阻力系数，它们只与马赫数(Ma)和攻角(AOA)有关。
[0102] 2.)再入弹道约束
[0103] 对于再入飞行器，典型的过程约束包括驻点热流密度、动压以及法向过载，它们通常可以表示为如下的形式：
[0104]
[0105] 其中，第一个方程表示驻点热流密度约束，它只与高度和速度有关；第二个方程表示为动压约束，同样，它只与高度和速度有关；第三个方程为过载约束，它与高度、速度以及攻角有关； qmax和nmax分别为飞行器所能承受的最大热流、动压和过载，它们与飞行器的结构配置以及任务需求有关，这三类约束在飞行过程中是必须满足的约束，因此，也称之为“硬约束”。
[0106] 另一类需要特别说明的约束为终端约束，终端约束条件的设置与传统水平着陆再入飞行器(如X-33，航天飞机)的终端约束条件一样，该条件为飞行器在离目标点位置特定的距离结束再入飞行段，并且将飞行器顺利地交班到终端区域能量管理(TAEM)制导系统中，因此，这需要飞行器在TAEM交班时刻满足正确的终端状态约束，本发明中所考虑的终端状态约束为高度、速度、弹道倾角和倾侧角，除此之外，还需保证飞行器在TAEM交班时刻的速度矢量方向要与目标视线方向的夹角尽量小，具体的约束条件如下：
[0107]
[0108] 其中，ef表示在TAEM交班时刻的相对能量，它的表达式为e＝-μ/r+v2/2：sf表示在TAEM交班时刻与目标之间的距离；Δψ表示航向误差角，它是航向角与目标视线角之间的夹角。
[0109] 3.)飞行器表述及其模型假设
[0110] 通用航空飞行器(CAV)是迄今为止最具代表性的高升阻比再入飞行器，仅仅依靠气动力控制，就能在无动力条件下滑翔穿越大气层。从Phillips的研究报告中获知，当前存在高升力CAV和低升力CAV两类通用航空飞行器，高升力CAV——CAV-H具有更大的升力系数和更高的升阻比，因此，本发明将以它作为研究模型，用以测试再入制导算法。CAV-H的质量为907Kg，特征面积为0.4839m2，最大升阻比在10度攻角附近，且接近3.5。为了使本发明的研究更为直观、简单，本发明将CAV-H的升力系数和阻力系数采用高阶的多项式函数进行拟合，获得以马赫数和攻角为输入的升力系数函数和阻力系数函数。除此之外，因为最大升阻比所在的攻角为10度附近，本发明将攻角的变化范围拓展到至5～20度。CAV-H在飞行过程中所能承受的最大热流、动压和过载分别为500W/cm2、100Kpa和2g。在再入飞行过程中，将飞行器的攻角变化设计为以相对能量为自变量的分段函数。
[0111]
[0112] 其中，emid定义为飞行器最后一次倾侧反转时刻的能量估计，aα1和aα2分别为多项式系数项，它们的取值通过以下两个条件确定，在emid时刻，攻角为10度，它对能量的导数为零。值得指出的是当emid不是制导输出时，emid是预设的固定值，在制导仿真过程中，其速度设为3500m/s，高度为35000m。
[0113] 4.)辅助惯性坐标系(AGI)
[0114] 后面将在建立的辅助惯性坐标系(AGI)中讨论飞行器的运动方程，因此，现主要介绍辅助惯性坐标系。该坐标系将固连在惯性空间，并且在每个制导周期根据飞行器的状态与估计的虚拟目标点的位置进行更新。虚拟目标点的更新也是根据当前飞行器的状态、终端速度和地球自转角速度进行。如下所示
[0115] λp＝λT+ω·tgo,φp＝φT,hp＝hT,tgo＝Sgo/(v+vf) (1.6)
[0116] 其中，λp、φp和hp分别表示虚拟目标点的经、纬度和高度；tgo表示剩余飞行时间，Sgo表示当前飞行器所在的位置与虚拟目标点位置所形成的大圆弧夹角，于是，以飞行器位置与虚拟目标位置在地球表面连线作为广义赤道就能确定该辅助惯性坐标系，实际上，AGI坐标系与飞行器所在的当地东北天坐标系(NED)只相差一个绕地心方向旋转的目标视线角。在AGI坐标系下的飞行器运动方程可以表示如下：
[0117]
[0118] 其中，和 分别为飞行器在AGI坐标中的经纬度；表示为绝对速度，是相对速度v矢量与地球自转角速度在飞行器所在位置所形成的牵连速度矢量的和；和 分别表示绝对速度矢量在AGI坐标系中的弹道倾角和航向角，它们的定义与东北天坐标系中的类似；绝对能量表示为如下形式：
[0119]
[0120] 根据AGI坐标系下的动力学关系，它的导数可以简单的表示为：
[0121]
[0122] (2)下落段制导：以最大攻角和零倾侧角下落，直到高度的变化率为零，进入滑翔段阶段。
[0123] 具体地，初始下降段是指飞行器从初始高度下落到高度变化率为零的一段高度范围，该段结束时，飞行器所受到的气动升力已经充分建立起来并能够支撑飞行器进行滑翔飞行。在这一阶段，由于飞行器速度很高并且大气密度急剧增加，攻角和倾侧角的调整将对飞行器受到的局部最大驻点热流和之后的滑翔弹道产生巨大的影响。一般情况下，局部最大热流将发生在该段结束时刻，因此，为了保证飞行器的局部驻点热流最小且飞行器具有最大的能量调整裕度，飞行器在该段将按最大攻角和零度倾侧角下落，如公式1.10所示。
[0124]
[0125] 其中αcmd是指令攻角，σcmd是指令倾侧角，αmax是最大攻角。
[0126] 初始下降段结束时刻的条件如下
[0127]
[0128] 其中，v和γ是飞行器当前的速度和弹道倾角，可以通过惯组(IMU)输出。
[0129] (3)滑翔段制导倾侧反转判断：根据当前能量状态与分段时间的关系，判断是否进行倾侧反转，当当前能量大于分段时刻的能量时，不进行反转，进入步骤(4)，当当前能量小于分段时刻的能量时，且不为最后一个反转点时，重新建立非线性参数控制问题，进入步骤(4)，当为最后一个反转点时，进入步骤(9)。
[0130] (4)滑翔段制导预测弹道积分：使用当前时刻的状态量作为积分初值，所述的初始控制参数和分段时间初值作为控制输入，其中当前时刻所对应的控制参数记为U0，进行预测弹道积分，获得广义标控脱靶量即终端状态偏差dψ以及全局多段弹道信息，Xk，Uk。
[0131] (5)精度判定：判断所述广义标控脱靶量的精度范围，如果所述dψ满足所述步骤(1)设置的终端计算仿真参数的精度要求，进入步骤(6)，如果不满足精度要求，进入步骤(8)。
[0132] (6)过程边界控制：判断所述的U0是否满足多种过程约束界限，当所述的U0超过边界控制Umax时，所述的U0等于Umax，否则，所述的U0等于U0，进入步骤(7)。
[0133] (7)执行控制指令：将所述的U0作用到弹道阻尼控制方程中，获得平稳滑翔的控制指令Uc，并将其作用到实际的系统中去，返回步骤(3)，并且将当前的分段时间作为倾侧反转的判断基准，当前的U0和分段时间作为下一步弹道积分的控制输入。
[0134] (8)滑翔段制导控制参数和分段时间更新：分别对预测的每段积分弹道进行线性化处理，结合多段伪谱法和变分法原理，获得满足终端标控脱靶量修正的控制参数和分段时间解析修正关系，进而获得控制参数和分段时间的更新，并进入步骤(6)。
[0135] 具体地，此步骤中采用一种全新的滑翔段预测-校正再入制导算法，该方法分别对预测的每段积分弹道进行线性化处理，结合多段伪谱法和变分法原理，获得满足终端标控脱靶量修正的控制参数和分段时间解析修正关系，进而获得控制参数和分段时间的更新。该方法不仅具有很快的计算效率——不需要求解非线性规划问题，而且只需很少的离散节点就能获得很高的计算精度。此外，在一次弹道积分条件下它能够同时提供关于纵向升阻比和倾侧角反转时刻的指令更新。
[0136] 1).降阶的再入动力学方程
[0137] 根据Yu和Mease的研究工作，将地球自转所产生的科式力和牵连加速度的影响考虑到AGI中的运动方程之中，并且假定飞行器在固定的球面上运动且弹道倾角对时间的导数足够小以至于能够忽略，飞行器的运动方程就能够转化为一组以绝对能量为自变量的降阶的运动方程，即构建降阶的再入动力学方程具体如下：
[0138]
[0139] 其中，e是相对能量，E是绝对能量； 是经度、纬度和航向角；F是总升阻比函数，它是以相对能量为自变量的函数，且在攻角确定情况下能够通过平稳滑翔条件离线获得；控制变量u，表示纵向升阻比，它是总升阻比在垂直方向上的一个分量；kld是由惯组输出的升阻比偏差估计，并且在预测积分过程中认为是常值；μ是引力常数，r为地心距，在计算过程中认为是常值，并且选择当前地心距与终端地心距的平均值作为该值。
[0140] 2).控制指令的参数化
[0141] 控制指令是根据平稳滑翔弹道的特点直接进行规划处理，因为攻角是通过离线规划获得，控制指令的规划也就是指将纵向升阻比曲线离散为几个参数的函数。实际上，因为纵向升阻比和倾侧角之间具有简单的一一映射函数关系，离散纵向升阻比曲线也就是离散倾侧角曲线。在滑翔过程中，将安排两次倾侧反转以保证飞行器具有足够的能力消除积累起来的横向误差，第一次倾侧反转安排发生在飞行中段，第二次倾侧反转安排发生在离TAEM段不远处。因此，纵向升阻比曲线就可以设计成如下一个分为三段的参数函数，其具体表达形式如下所示：
[0142]
[0143] 其中，E0和Ef分别代表初始和终端绝对能量，ef代表终端相对能量，Ere1和Ere2分别表示第一次倾侧反转和第二次倾侧反转发生时的绝对能量，k1为一个常数，表示两次倾侧反转之前，纵向升阻比都按常值设计，选择常值的主要原因就是保证飞行器能够以平稳的状态滑翔，并且留有足够大的调整能力从而克服纵向存在的误差。在第二次倾侧反转之后，纵向升阻比的曲线设计为一个典型的二次曲线，并且它的终端值等于总升阻比曲线的终端值，c1和c2是二次曲线的参数，它们的取值通过分段函数光滑，导数连续两个条件确定，具体如下：
[0144]
[0145] 于是，获得关于定c1和c2的解析表达式：
[0146]
[0147] 因而，当k1，Ere1和Ere2给定，纵向升阻比曲线也就确定，当初始的状态已知时，就能够通过对降阶的非线性运动方程进行积分而确定全部再入滑翔弹道。
[0148] 3).多段非线性方程的线性化
[0149] 将具有终端约束的降阶再入运动学方程，它能够表示为一般的非线性动力学方程的形式：
[0150]
[0151] 其中，所述 如果初始的k1，Ere1和Ere2给定，就能通过欧拉法或龙格库塔数值积分方法预测全局的状态量，那么预测的终端状态与实际所需的终端状态之间的终端误差就可以表示为δxf＝x(Ef)-xf，将方程(1.16)在预测弹道周围进行高阶泰勒展开，于是一组以状态偏差为自变量的误差传播动力学方程就能通过忽略二阶以上高阶项获得，具体如下：
[0152]
[0153] 其中，x＝xref-δx,u＝uref-δu，矩阵A为3×3的矩阵，即A1、A2、A3是系数矩阵，为3×3的矩阵，矩阵B为3×1的向量，即B1、B2、B3是系数矩阵，为3×1的向量；δu为标量。通过推导，获得矩阵A和矩阵B中各个元素的具体表达式如下：
[0154]
[0155] 和
[0156]
[0157] 由于存在倾侧反转，矩阵B1中第三个元素的第一项符号与矩阵B3中相应的符号一致，与矩阵B2中相应的符号相反。此外，δx(E0)是与积分弹道之间的初始状态偏差，并且，由于弹道积分预测在每一个制导周期中都是以当前状态为初值进行积分，于是，初始状态偏差为零。
[0158] 4).多段线性伪谱校正算法
[0159] 对于线性伪谱校正，首先，将非线性动力学方程在预测的积分弹道附近进行拟线性处理，从而建立一个以偏差为自变量的线性微分方程，之后对该线性微分方程的状态变量和控制变量都通过一组拉格朗日插值多项式进行逼近，微分动力学约束也将通过正交配点转化为一组多变量的代数约束，最终将求解非线性方程的控制问题转化为一个连续求解线性代数方程的问题，该方法相对于直接求解非线性规划问题只需很短的时间就能完成求解。在本发明中，线性伪谱校正将拓展到求解具有控制变量参数化的多段多点边值问题。正如之前介绍的一样，当前伪谱法的研究主要是基于三类正交节点，即Legendre-Gauss节点(LG)、Legendre-Gauss-Radau节点(LGR)以及Legendre-Gauss-兰伯特节点(LGL)，并且，这些节点都是Legendre多项式或者Legendre多项式与它的导数所形成的新的多项式的根，当前的研究表明，由于LGL的微分逼近矩阵具有奇异性，LGL伪谱离散将会具有一个不完备的求解空间，而LG伪谱离线和LGR伪谱离散由于微分逼近矩阵是满秩而不存在这样的缺陷。因此，选择LG正交离散(Gauss伪谱法)作为离散的基本方法开展研究，下面的内容将具体地介绍离散过程。
[0160] 因为拉格朗日插值多项式的节点都是正交节点，且都是分布在[-1，1]之间，离散的第一步就是将各段的时域从[t0，tf]转换至[-1，1]上，转换函数为简单的线性插值函数，如下所示：
[0161]
[0162] 于是，线性动力学方程(1.17)就能够转化为下面的形式
[0163]
[0164] 为了直观理解该方程和方便推导，其中的 表示为如下的形式：
[0165]
[0166] 本发明将LN(τ)定义为N阶的拉格朗日插值多项式，τi表示N阶Legendre多项式的根，也就是之前提到的LG节点。需要注意的是，这里没有解析求解N阶Legendre多项式根的方法，必需通过数值算法才能获得LG节点。于是，就可以将线性动力学方程(1.21)中的状态量表示为一组由LG节点为支撑点所形成的拉格朗日插值多项式基来线性组合如下：
[0167]
[0168] 根据拉格朗日插值多项式的基本理论，LN(τ)满足以下的基本性质
[0169]
[0170] 并且
[0171] δxN(τl)＝δx(τl) (1.25)
[0172] 之后，LG节点上的状态导数值就能通过微分逼近矩阵转化为一组代数方程，并且建立起与各个LG节点上状态量之间的显示关系。其中，微分逼近矩阵D是一个N×N+1的矩阵，通过对拉格朗日插值多项式的各个元素分别求导获得，D矩阵各元素的具体表达形式以及状态量与其导数的代数关系如下所示：
[0173]
[0174] 其中，
[0175]
[0176] 假定每一段的状态量在LG节点上表示为如下形式
[0177]
[0178] 将方程(1.26)带入到方程(1.21)中，那么，线性动力学方程(1.21)不仅能够转化为一组代数约束，并且能够表示为LG节点上各个状态的函数，具体如下：
[0179]
[0180] 其中，k＝1，2，……，N。以方程(1.29)中的第一段代数约束关系式为列，重新组合微分逼近矩阵D，一个以状态矢量形式表示且更为简单的关系式如下所示：
[0181]
[0182] 其中，D1和D2:n的具体表达形式如下：
[0183]
[0184] 其中，s表示为状态变量的个数，矩阵A1和B1分别表示如下：
[0185]
[0186] 再次对方程(1.30)进行重新组合，在LG节点上的各个状态量可以解析地表示为：
[0187]
[0188] 需要指出的是在Gauss伪谱离散过程中，动力学约束关系式只约束在LG节点上而不包含终端节点，因此，δx也只包含了LG节点上的状态偏差值而不包含终端状态偏差。实际上，终端状态偏差可以通过Gauss型积分公式而表示为初始状态偏差和在LG节点上的状态偏差的显示函数关系如下所示：
[0189]
[0190] 其中，ω是Gauss型积分的权函数，它是通过数值算法计算获得，将方程(1.33)带入到方程(1.34)中并重新整理，第一段的终端状态偏差就能表示为如下的形式：
[0191]
[0192] 其中，1表示的为单位列向量，W1是Gauss型积分公式的权函数矩阵，其结构如下所示：
[0193]
[0194] 最终，将方程(1.33)带入到方程(1.35)中，第一段的终端状态偏差就能表示为关于初始状态偏差δx0和控制偏差δu的显式解析表达式，其表达形式如下：
[0195]
[0196] 其中，系数矩阵 是一个3×3的方阵，而系数矩阵 是一个3×1的向量，它们中的元素能够表示为：
[0197]
[0198] 同样，类似于第一段的推导，获得第二段终端状态偏差的显示表达式，并且将第一段的终端状态偏差作为第二段的初始状态偏差，其中系数矩阵 和 的表达形式如下：
[0199]
[0200] 同样，第三段终端状态偏差中的系数矩阵 和 也能表示为：
[0201]
[0202] 于是，当方程(1.38)、(1.39)和(1.40)确定时，就能根据链式法则容易地建立起关于初始状态偏差、控制偏差以及终端状态偏差之间的函数关系，其具体表达式如下：
[0203]
[0204] 重新整理可以表示如下：
[0205]
[0206] 预测的弹道积分都是从当前状态开始的，因而，初始的状态偏差在每个制导周期都为零，于是，就能直接通过系数矩阵获得关于控制偏差与终端状态偏差之间的修正关系，通过计算获得控制偏差就能在线性化的模型框架下修正所存在的终端状态偏差如下：
[0207]
[0208] 一般来说，该修正方法与多段动力学方程的分段个数无关，因而，该方法能够推广到具有强终端状态约束的任意段多段多边值修正问题中。
[0209] 观察纵向升阻比曲线的方程(1.13)，控制变量u被设计为仅仅只有一个参数，k1的参数函数，因此，调整该控制变量只能保证一个终端状态约束的修正，实际上，纵向升阻比对于经度的灵敏度远高于其他的状态量，自然而然地选择控制变量对经度进行修正，其修正关系如下：
[0210]
[0211] 其实，纵向的距离要求通过调节纵向升阻比就能满足，为了保证飞行器能够准确的飞向目标，还需要保证飞行器的横向距离要求，可以通过调节最近一个倾侧反转点的发生时机来保证。结合之前推导的线性伪谱校正公式以及变分法，一个全新的用来修正最近的倾侧反转点从而消除积累的纬度误差的校正方法如下：
[0212]
[0213] 假设在小量δu和δE的作用下，终端状态量可以表示为：
[0214]
[0215] 当忽略掉高阶项，并且假定在校正更新的作用下，终端的状态为所需的经纬度时，可以获得如下的表达形式：
[0216]
[0217] 需要注意的是，其中 和 分别为不考虑反转点变化时，纵向升阻比对终端经纬度的偏导数，可以通过之前的线性伪谱校正公式直接给出，具体如下：
[0218]
[0219] 而 和 分别为不考虑纵向升阻比变化时，倾侧反转点对终端经度和纬度的偏导数，其表达式的推导过程如下所示。
[0220] 假设在反转发生时刻，反转点能量的变化将引起反转点初始状态的改变，如图2所示。
[0221] 在线性动力学方程假设情况下，倾侧反转前后，由反转点改变所引起的状态量偏差分别为 和 那么，由反转点改变所引起的原反转点状态的等时变分可以表示为
[0222]
[0223] 于是，终端状态偏差可以表示为
[0224]
[0225] 进而，
[0226]
[0227] 最终，结合公式(1.47)、(1.48)和(1.51)就能获得关于纵向升阻比和倾侧反转能量时刻的修正就可以解析的表达为：
[0228]
[0229] 迭代更新产生的纵向升阻比和最近的倾侧反转点就可以表示为：
[0230]
[0231] 当控制修正量确定之后，就能根据当前的升阻比确定倾侧角的幅值大小：
[0232]
[0233] 其中，(L/D)IMU为通过惯组单元测量获得的当前升阻比大小。
[0234] 重复上面的预测积分和控制修正步骤就能获得满足终端要求的控制量和最近的倾侧反转点，该方法实施过程简单，且将多段伪谱法与多段非线性动力学系统的线性化、模型预测控制和变分原理结合起来，因而，不仅保留了正交多项式离散的高精确性，也具有不需要求解非线性规划问题的高效性，此外，该方法能够在一个弹道预测积分内同时完成关于纵向升阻比和最近的倾侧反转点的修正。
[0235] 5).过程约束控制策略
[0236] 本发明在执行过程中，将驻点热流、动压或过载等过程约束将转化为倾侧角的幅值约束，通过限制实际倾侧角的幅值来达到对过程约束的严格控制。以最大驻点热流约束为列，最大允许的倾侧角幅值通过以下公式计算获得
[0237]
[0238] 其中，rmin为当前速度给定情况下，通过驻点热流约束函数所确定的最小高度。其中， L为加速度计测量获得的当前升力大小，为通过惯组单元估计的当前升力大小，dr/dv和(dr/dv)|Q分别表示飞行器质点动力学方程与驻点热流密度函数所确定的高度对速度的导数，当热流密度约束起作用时，kcst为一个正值，其它时刻为零。实际上，其它约束函数也都为高度与速度的函数，相应的最大倾侧角幅值边界也都能通过同样的方法确定，因此，实际的倾侧角幅值将限制到过程约束所确定的倾侧角幅值边界之内。
[0239] σ(v)≤σmax(v) (1.56)
[0240] 6).弹道阻尼抑制策略
[0241] 再入飞行器(尤其是高升阻比的再入飞行器)的再入弹道具有弱阻尼的振荡特性，直接将获得攻角指令和倾侧角指令作用到弹体将会产生弹道的周期振荡，于是，之前用于计算制导指令的降阶动力学方程假设将不再使用，由Yu和Chen提出的弹道阻尼控制技术能够很高的解决抑制弹道跳跃的问题，该方法简单易行，也是EGAS的主要组成部分。在这一章内，该方法将用于提供能够保持滑翔飞行的指令攻角和指令倾侧角，考虑忽略地球自转影响的弹道倾角二阶动力学方程。
[0242]
[0243] 在考虑飞行器在再入飞行的任何阶段，当攻角指令和倾侧角指令给定时，一个能够维持弹道倾角二阶导数为零的特殊的弹道倾角可以通过公式(1.58)确定。
[0244]
[0245] 其中，上标为“*”的项表示由攻角指令所确定的升力与阻力或升力系数与阻力系数，之后，弹道跳跃的抑制就通过在纵向升力方向加入该特殊弹道倾角的负反馈，并保持横向升力不变而实现，实际作用的倾侧角和升力系数就通过公式(1.59)确定。
[0246]
[0247] 其中，下标为“2”的项表示实际需要作用的倾侧角和以实际作用的攻角所确定的升力系数，Kγ表示为确定弹道抑制效果的负反馈增益，并且通过几次简单的仿真就能确定。之后，再次整理公式(1.59)就能获得升力系数与倾侧角的解耦显示表达式：
[0248]
[0249] 之后，不像Yu在文献中采用泰勒级数展开的方式计算实际作用攻角，牛顿法将用于直接求解所需的实际作用攻角，其中，升力系数已经拟合为一个关于攻角和马赫数的多项式函数：
[0250]
[0251] 当使用上一步获得攻角作为初始值代入计算，只需消耗很少的时间(几毫秒)，仅仅几步计算就能获得满足精度要求的解，值得指出的是在实际的再入过程中，升力系数还需乘以一个由加速度计和惯组单元所估计的升力偏差系数。
[0252] (9)末制导段显示制导：分别根据当前状态和终端约束条件在纵、横平面采用比例导引和多项式制导，获得能够导引飞行器安全地飞向目标且满足多种终端约束的制导指令。具体地，
[0253] 1).横向比例制导律
[0254] 比例导引作为最著名的且最具鲁棒性的制导律，将在这一章节内被推广到末端调整段的横向导引中，从而导引飞行器在横向平面飞向目标，并且保证终端的横向加速度收敛到零，正如比例导引的显示表达式所示，比例导引的制导加速度指令是通过视线角变化率、投影在横向平面的速度以及比例导引常数所确定。
[0255]
[0256] 2).纵向多项式制导律
[0257] 在这一节中，由Kim提出的多项式制导策略将用于导引飞行器在纵向平面飞向TAEM交班点，并且满足高度、弹道倾角等纵向状态约束。其制导指令表示如下：
[0258]
[0259] 其中，sf为TAEM交班点与目标之间的距离，并且制导增益m和n分别为满足n>m>0的常值，这样，能够导引飞行器在纵向平面飞向TAEM交班点，并且满足终端高度和弹道倾角约束的纵向指令就已经获得，值得注意的是，由于希望弹道的形状尽量简单、光滑从而使飞行器的飞行状态平稳，导引增益的选择也将是使指令多项式的阶数尽可能的少，因此，在之后的所有仿真过程中，我们将导引增益设置为m＝2，n＝1。
[0260] 3).纵、横向制导平面的指令分配
[0261] 当采用比例导引和多项式制导方法分别在纵向和横向两个平面内确定所需的制导指令后，剩下的问题就是如何调整攻角和倾侧角以实现两个制导平面内的制导指令，让我们考虑保留地球曲率并且忽略地球自转的再入动力学方程，纵向加速度和横向加速度就能够分别的表示为
[0262]
[0263] 由于侧向加速度和纵向加速度都是总加速度的分量，于是，获得指令倾侧角的表达式：
[0264]
[0265] 之后，采用牛顿迭代就可以很快的确定所需的指令攻角如下所示：
[0266]
[0267] 其中，Cl为升力系数，Clα为升力系数对攻角的导数；k表示第k次迭代过程。
[0268] 实施例2
[0269] 在本实施例中，考虑模型不确定性和随机拉偏的Monte Carlo仿真打靶试验将用于验证本发明实施例1中提出的自适应分段的多段线性伪谱广义标控脱靶量在线规划再入制导方法的性能。仿真模型为实施例1中，初始化步骤第3)小节中介绍的CAV-H。考虑的随机偏差包括初始的状态偏差(高度、经纬度、速度、弹道倾角和航向角)、气动拉偏以及大气环境拉偏，初始状态和终端状态如表1所示，此外，还必须保证终端的倾侧角和航向误差角收敛到零，随机干扰的3-sigma值如表2所示。制导指令的更新周期是1Hz，所有程序都在3.3Ghz的个人电脑中安装的MATLAB2008b中运行，采用更加高效的编译环境能够提高计算效率。
[0270] 表1初始状态和终端状态
[0271]
[0272] 表2蒙特卡洛打靶的随机干扰设置
[0273]
[0274] 需要指出的是，虽然不确定性和随机干扰模型比较简单，但足以校验所设计制导律的制制导性能和鲁棒性，最大的升阻比拉偏将达到26.09％(-15％的升力拉偏和15％的阻力拉偏)，于是，1000次Monte Carlo打靶的结果如图3-16所示。其中，图3为1000次蒙特卡洛打靶的地面轨迹，显而易见，所有的弹道都飞向TAEM交班位置，并且在各种拉偏影响下，横向运动散布很大，最大的横向机动距离到达15度(1670km)，最小的横向机动距离仅仅为5度(557km)。图4为1000次蒙特卡洛打靶的高度随速度变化曲线，可见，在综合拉偏影响下，飞行器会按照合适的滑翔高度飞行，因此，高度散布也比较大，达到近5km的范围，但最终都收敛到预设的终端高度上。图5为1000次蒙特卡洛打靶的攻角曲线，可见，滑翔段的攻角在规划的10度附近变化，并且能够清晰的看出倾侧角反转的影响，由于末段调整段的初始高度散布较大，这段时期的攻角也会具有较大的散布，从而保证终端高度。图6为1000次蒙特卡洛打靶的倾侧角曲线，可见，为了适应不同的综合拉偏，倾侧角的幅值在30度到60度之间变化，并且两次倾侧发转的发生时刻也能在这图上清晰的看到，终端的倾侧角都收敛到零。图7为1000次蒙特卡洛打靶的航向误差角随能量的变化，由于倾侧反转将导致航向误差角的导数符号发生变化，因此，在这图上能够容易的看出两次倾侧反转发生的能量高度，并且在综合拉偏影响下，航向误差角的变化散布很大，最大的航向误差角达到了45度左右，而最小的航向误差角才在10度附近，且航向误差角在终端都收敛到零。图8为1000次蒙特卡洛打靶的速度随时间变化曲线，可见，初始速度具有较大的散布，但是，速度最终都在2500m/s附近。图9、图10和图11分别为1000次蒙特卡洛打靶的驻点热流密度、动压和过载随时间的变
2
化曲线，三种过程约束都处于安全范围之内，最大热流密度已经贴着500W/cm的热流极限，最大的动压不超过80Kpa，最大的过载也不超过1.6g。图12为1000次蒙特卡洛打靶的弹道倾角随时间的变化曲线，终端结果都收敛到零。
2
化曲线，三种过程约束都处于安全范围之内，最大热流密度已经贴着500W/cm的热流极限，最大的动压不超过80Kpa，最大的过载也不超过1.6g。图12为1000次蒙特卡洛打靶的弹道倾角随时间的变化曲线，终端结果都收敛到零。
[0275] 图13为1000次蒙特卡洛打靶的升阻比随能量的变化曲线，可见，在综合拉偏下，总升阻比(可用升阻比)在2.5到4.5之间变化，而纵向升阻比也在2到2.5之间变化，可见，所设置的拉偏还是对再入飞行具有相当大的影响。图14、图15和图16分别为1000次蒙特卡洛打靶的终端倾侧角和终端高度的散布统计，终端航向误差角和终端速度偏差以及终端弹道倾角的散布统计。很明显，终端的倾侧角偏差不会超过1deg，所有的高度偏差也不会超过3m，由于采用了比例导引，终端的航向角误差都小于10-2deg，终端速度偏差也都处于马赫数0.1之内。值得指出的是，高度、倾侧角和弹道倾角的终端偏差会出现偏置是由于选择TAEM交班点之后50m处作为制导解算的虚拟目标点所致，表3为蒙特卡洛打靶终端状态的统计值，其终端误差都非常小。
[0276] 表3蒙特卡洛打靶终端状态统计
[0277] 均值 方差/均方差 50％几率落入值
δhf(m) 1.2663 0.0885/0.2975 1.2584
δvf(m/s) -9.5902 76.7625/8.7614 7.0148
δγf(Deg) -0.0158 1.3352E-5/0.0037 0.0158
δψf(Deg) 0.0033 5.6597E-7/7.5231E-4 0.0032
δσf(Deg) 0.2013 0.0158/0.1258 0.1864
δhf(m) 1.2663 0.0885/0.2975 1.2584
δvf(m/s) -9.5902 76.7625/8.7614 7.0148
δγf(Deg) -0.0158 1.3352E-5/0.0037 0.0158
δψf(Deg) 0.0033 5.6597E-7/7.5231E-4 0.0032
δσf(Deg) 0.2013 0.0158/0.1258 0.1864
[0278] 通过上述算例表明，本发明不仅具有很快的计算速度和求解精度，而且能够保证多种终端误差，同时，还能在线确定倾侧角反转时刻，从而充分发挥飞行器的巨大横向机动能力，相对传统方法将大大较少倾侧反转次数。
| 标题 | 发布/更新时间 | 阅读量 |
|---|---|---|
| 火星大气进入段轨迹最优跟踪制导方法 | 2020-05-08 | 310 |
| 一种无人机 | 2020-05-13 | 879 |
| 一种多约束火星大气进入预测制导方法 | 2020-05-15 | 251 |
| 一种高超声速飞行器及其弹道设计 | 2020-05-13 | 699 |
| 强约束复杂任务条件下的三维剖面优化设计方法、系统及介质 | 2020-05-14 | 743 |
| 一种垂直起降水陆两栖地效飞行器 | 2020-05-11 | 566 |
| 可回收探空仪、系统及回收方法 | 2020-05-13 | 63 |
| 一种基于马格努斯效应的分离圆柱式减摇装置 | 2020-05-15 | 928 |
| 一种高升阻比的低速无人机翼型 | 2020-05-08 | 316 |
| 一种风力助航转筒 | 2020-05-08 | 669 |
高效检索全球专利专利汇是专利免费检索,专利查询,专利分析-国家发明专利查询检索分析平台，是提供专利分析，专利查询，专利检索等数据服务功能的知识产权数据服务商。
我们的产品包含105个国家的1.26亿组数据，免费查、免费专利分析。
分析报告
专利汇分析报告产品可以对行业情报数据进行梳理分析，涉及维度包括行业专利基本状况分析、地域分析、技术分析、发明人分析、申请人分析、专利权人分析、失效分析、核心专利分析、法律分析、研发重点分析、企业专利处境分析、技术处境分析、专利寿命分析、企业定位分析、引证分析等超过60个分析角度，系统通过AI智能系统对图表进行解读，只需1分钟，一键生成行业专利分析报告。
QQ群二维码
意见反馈
意见反馈