一种考虑多禁飞区约束的协同解析再入制导方法与流程

本发明涉及一种考虑多禁飞区约束的协同解析再入制导方法，属于航天技术、武器技术、制导控制领域。

背景技术：

协同编队飞行技术对于提升体系作战效能具有重要意义，拥有广阔的应用前景，目前已经在无人机领域得到了广泛而深入的研究。但是，对于高超声速滑翔飞行器，编队飞行十分困难。因为此类飞行器初速很高，但无动力，需要通过适当的横向机动精确管理飞行能量。但对于具有不同起滑状态、受到不同程度干扰的飞行器，其速度衰减曲线和横向机动幅度均不相同，因而难以组成相对位置比较稳定的飞行编队。

在协同编队飞行技术方面，常见的编队形式有楔队、梯队、横队、纵队和v形等，可以实现协同侦察、防御和进攻等复杂任务。此外，通过利用领机产生旋涡中的上升气流，编队飞行还可以有效节省燃料，增加航程。但是，受速度和横向机动幅度的时变特性的制约，目前还没有可以有效应对高超声速滑翔飞行器协同编队飞行问题的理论和方法。

yuw.，chenw.，analyticalentryguidanceforno-fly-zoneavoidance(针对禁飞区规避的解析再入制导律)[j]aerospacescienceandtechnology，2018，72：426-442.提出了具有高精度的地球旋转补偿模型，获得了性能更好的三维再入弹道解析解，然后基于解析解设计了渐增式的倾侧反转序列规划方案，实现了以较少的倾侧反转代价应对大量的禁飞区约束的目标。在本发明的以下说明中简称为文献一。

gaxleyc.，atmosphericentry(大气层再入)[m].therandcorporation，1960.针对固定地球模型，给出了常值纵向升阻比条件下的飞行时间解析解。在本发明的以下说明中简称为文献二。

技术实现要素：

本发明的目的是为了解决目前还没有可以有效应对高超声速滑翔飞行器协同编队飞行的问题，提出一种考虑多禁飞区约束的协同解析再入制导方法，解决在多禁飞区环境中导引多个高超声速滑翔飞行器同时到达目标的协同飞行问题。

本发明的主要内容分为三个方面：(1)针对旋转地球模型，推导了高精度的飞行时间解析解；(2)针对单个飞行器，设计了考虑多禁飞区和抵达时间约束的解析再入制导方法；(3)针对多个飞行器，设计了多禁飞区环境中导引多个高超声速滑翔飞行器同时到达目标的协同飞行方案。其中，在解析再入制导方法设计时，控制终端速度和抵达时间均由纵向升阻比控制，因此在控制上是一个欠驱动问题，本文利用纵程解析解和时间解析解设计了一种能够兼顾能量管理要求和抵达时间要求的纵向参考剖面规划方法，可以有效解决此问题。

本发明一种考虑多禁飞区约束的协同解析再入制导方法，包括以下几个步骤：

步骤1：再入制导问题描述

这里采用圆球形旋转地球下的三维再入质点动力学模型，其中利用经度λ、纬度φ和高度h描述飞行器位置信息，而用速率v、弹道倾角γ和航向角ψ描述速度信息。

极高的飞行速度导致力热环境极其恶劣。因此，为了保证飞行器的各个子系统正常工作，再入飞行轨迹需要满足热流密度来流动压q、过载n约束。考虑到飞控系统能力有限，需要对攻角和倾侧角范围及变化率进行限制。除了常规约束，这里进一步考虑多个圆形禁飞区约束。

再入段飞行在飞行器到目标的水平距离为staem时终止。此时期望的终端高度为htaem，终端速度为vtaem，终端航向误差|δψtaem|和终端倾侧角|σtaem|满足约束条件。此外，这里还指定期望的终端时刻为ttaem。

步骤2：基于旋转地球模型的飞行时间解析解

为了应对抵达时间约束，这里推导飞行时间解析解。能量定义如下

其中，v是速度，h是高度，re是地球平均半径，μ是引力常数。能量对时间t的导数为

式中g为重力加速度。在旋转地球背景下，有如下复杂的非线性方程

其中，γ是弹道倾角，d是阻力，m是质量，ωe是地球自转角速度，φ是纬度，ψ是航向角。公式(4)中右边第3项第4项是由地球自转引起的惯性力分量。本文将沿速度方向的惯性力分量视为为附加阻力，如下

定义等效阻力为将公式(3)～(5)代入公式(2)，并取其倒数，可得

定义纵向升阻比为其中l1＝lcos(σ)是升力在纵向平面内的分量，规划为

其中a0，a1，a2为二次多项式系数。

现在利用剖面估计平稳滑翔状态下的等效阻力(记为)。弹道倾角时间导数如下

假设可估算平稳滑翔纵向升力，如下

其中，δl1是惯性力沿垂向的分量，被视为附加纵向升力，公式如下

与此同时，采用上述假设，δd也可以被简化为

则可以由下式估算

将公式(12)代入公式(6)，并利用e替换v可得

由于h＜＜re，故可令h取滑翔高度的中间值h^*，进而令r^*＝re+h^*。尽管惯性力是小量，但是当速度接近第一宇宙速度时，惯性力有可能会造成上式分母为零，从而发生奇异。为了避免发生上述奇异情况，这里将惯性力看作小量，并进行一阶taylor展开，如下

其中，简化符号hz1、hm的定义如下

hm＝2e+μ/r^*(16)

根据δl1和δd的表达式，可将hz1分为两个部分hpaf1和hpaf2，如下

hz1＝hpaf1+hpaf2(17)

其中，

hpaf1＝-2r^*vωecos(φ)sin(ψ)(18)

由于hpaf2的影响要远小于hpaf1，因此，可以粗略地采用线性函数对hpaf2进行近似，如下

hpaf2(e)≈kpaf2(1)v+kpaf2(0)(20)

其中，系数kpaf2(1)和kpaf2(0)由两个端点确定，如下

其中，hpaf2(e)由当前状态代入公式(19)计算得到，而hpaf2(etaem)是将期望的终端状态代入公式(19)计算得到。

实际上，高超声速滑翔飞行器围绕记为广义赤道的某一个大圆飞向目标，飞行器围绕广义赤道左右横向机动，广义赤道的方位角可以看作是飞行器航向角的某种加权平均值，因此，在后续航向角曲线尚且未知的情况下，本发明利用广义赤道的方位角替代hpaf2(e)和hpaf2(etaem)公式中的航向角。至此，在公式(14)中，除了自变量e之外，其它参量均为常值，进而可以对公式(14)进行积分。经大量推导、整理，可得飞行时间解析解的表达式如下

其中，系数kt(1)、kt(2)、kt(3)、kt(4)、kt(5)、kt(6)、kt(7)的表达式如下

步骤3：考虑多禁飞区和抵达时间约束的解析再入制导方法

此步针对单个飞行器，设计考虑多禁飞区和抵达时间约束的再入制导方法。根据弹道特点，再入过程一股分为三个阶段：下降段、平稳滑翔阶段和高度调整阶段。

s31：下降段

为了避免掉进稠密大气层而引起热流密度过大，下降段飞行器以最大可用攻角、零倾侧角下滑。当升力足以支撑飞行器平稳滑翔时，攻角平滑过渡到基准攻角，随即进入平稳滑翔阶段。

s32：平稳滑翔阶段

平稳滑翔阶段是最长、最重要、也是最复杂的再入飞行阶段，此阶段的制导方案利用三维再入弹道解析解和飞行时间解析解在线快速规划满足禁飞区和抵达时间约束的参考弹道。平衡滑翔阶段的制导流程设计如下：

s321、设计一条基准攻角剖面，并确定相应的基准升阻比剖面；

s322、考虑能量管理要求和抵达时间约束，合理设计参数化的纵向升阻比剖面，同时可以确定参数化的横向升阻比剖面；

s323、为了补偿地球自转的影响，在和剖面的基础上，建立等效纵、横向升阻比(和)剖面；

s324、按照能量管理和抵达时间要求，利用文献一中的纵程解析解和所述步骤2得到的飞行时间解析解求解纵向剖面参数；

s325、每隔60s，利用纵程解析解及文献一中的横程解析解根据禁飞区和终端位置约束更新一次倾侧反转序列。

s326、调节基准倾侧角跟踪等效纵向升阻比剖面，并按计划倾侧反转；

s327、为了抑制弹道振荡，将弹道阻尼反馈引入到基准攻角和倾侧角中，从而得到制导指令。

s328、按照过程约束限制倾侧角指令；

s329、从步骤s322开始重复上述过程，直到平稳滑翔阶段结束。

s33：高度调整阶段

在最后一次反转点(记此节点对应的飞行器能量为ebr(ntr))之后，飞行器进入高度调整阶段，但是在最后一次反转点之前的某个节点(记为对应的飞行器能量为est)，飞行器就开始了高度调整阶段的准备工作。这里本发明设计了一种基于在线弹道仿真的多目标数值迭代规划方案。在此方案中，利用方向导数改进了拟牛顿法迭代算法，从而减少了弹道仿真次数，可大幅提高迭代收敛速度。所述的多目标数值迭代规划方案的流程为：

s331、利用弹道仿真预测终端时间和终端状态；

s332、判断终端速度或终端时间是否满足要求；若是则跳转步骤s333，否则跳转步骤s334；

s333、根据终端速度和终端时间要求，利用拟牛顿迭代法调整纵向升阻比剖面和最后一次反转点，跳转步骤s331；

s334、判断终端高度是否满足要求；若是，则跳转步骤s336，否则跳转步骤s335；

s335、根据终端高度误差修正基准攻角参数，跳转步骤s331；

s336、完成迭代规划。

当e＝est时，本制导方法利用此多目标数值迭代规划方案规划纵向升阻比、剩余飞行距离、以及飞行时间相对能量的参考剖面，并确微调后续轨迹，以满足终端时间、速度、高度要求。由于仅需要对后续一小段弹道进行数次仿真，所以并不会造成较大的计算负担。当est＞e＞ebr(ntr)时，本制导方法跟踪在e＝est时刻利用多目标数值迭代方案规划获得的纵向升阻比、剩余飞行距离、以及飞行时间相对能量的参考剖面，而当ebr(ntr)＞e时，则利用比例导引律确定基准倾侧角，以消除航向误差，并通过微调攻角跟踪飞行时间和剩余飞行距离参考剖面，以保证满足终端时间和速度要求。

步骤4：多飞行器抵达时间协同方案

一股情况下，飞行器抵达目的地的时间受两个重要因素影响：发射时间以及再入飞行过程中对时间的控制，其中前者用于对抵达时间进行大范围调整，而后者用于克服干扰，对时间进行精确微调。但是，由于本发明并不考虑助推段和下压段的飞行过程，在这里利用再入段的起始时间替代发射时间，利用再入段的终端时间替代抵达时间。在实际情况中，则需要综合考虑这三段的飞行过程。尽管助推过程中的不确定因素会引起再入段起始时刻和起始状态的偏差，但是利用上文介绍的制导方法可以在再入飞行中有效克服这些偏差因素。

常见的时间协同问题有：同步到达同一目标，异步到达同一目标，同步到达不同目标，异步到达不同目标，以及它们之间的组合问题。由于这些问题的应对策略类似，下面以同步到达同一目标的问题为例介绍协同方案。

假设有nhgv个飞行器，记为hgvi，i＝1，2，…，nhgv。在给定发射点和目标点之后，根据助推火箭性能和所选主动段制导方案，可以确定飞行器的再入段起始状态。之后，仅根据能量管理要求，利用纵程解析解确定相应的纵向升阻比剖面参数。进而，利用离线弹道仿真可以预测所有飞行器的飞行时间tef(hgvi)，i＝1，2，…，nhgv，从中选出时间最长的一个飞行时间tef(max)，并预留发射准备时间tpre和助推飞行时间tboost，可以确定期望的抵达时间ttaem＝tpre+tboost+tef(max)，进而利用ttaem-tef(hgvi)可以确定各个飞行器的再入段起始时间。之后，根据所选的助推段制导方法，可以确定相应的发射时间。

本发明的优点在于：

(1)发明了基于旋转地球模型的高超声速滑翔弹道飞行时间的精确解析解，预测误差保持在3％以内，可适用于非常值纵向升阻比剖面的情况。而传统的预测方法由于未考虑地球旋转的影响，且只适用于常值纵向升阻比的情况，误差较大；

(2)可基于三维再入弹道解析解和飞行时间解析解，在线快速规划满足禁飞区、抵达时间、终端能量等多约束的参考弹道，解决了在多禁飞区环境中导引多个高超声速滑翔飞行器同时到达目标的协同飞行问题。

(3)设计了考虑终端时间、速度和高度要求的基于在线弹道仿真的多目标数值迭代规划方案。同时，为了减轻计算负担，利用方向导数改进了拟牛顿法，减少了弹道仿真次数。

附图说明

图1是本发明工作流程示意图；

图2是基准升阻比剖面曲线；

图3是高度调整阶段的多目标数值迭代规划方案流程图；

图4是飞行器轨迹的地面投影；

图5是飞行器的高度-速度曲线；

图6是飞行器的攻角-时间曲线；

图7是飞行器的倾侧角-时间曲线；

图8是飞行器的剩余飞行距离-能量曲线；

图9是飞行器的飞行时间-能量曲线。

图10a、b是montecarlo仿真的飞行器再入轨迹；

图11a、b是montecarlo仿真的飞行器高度-速度曲线；

图12a、b是montecarlo仿真的飞行器攻角曲线；

图13a、b是montecarlo仿真的飞行器倾侧角曲线；

图14a、b是montecarlo仿真的飞行器终端速度和终端航向误差分布情况；

图15a、b是montecarlo仿真的飞行器到达时间分布情况。

具体实施方式

下面将结合附图和实施例对本发明做进一步的详细说明。

本发明的主要内容分为三个方面：(1)针对旋转地球模型，推导了高精度的飞行时间解析解；(2)针对单个飞行器，设计了考虑多禁飞区和抵达时间约束的解析再入制导方法；(3)针对多个飞行器，设计了多禁飞区环境中导引多个高超声速滑翔飞行器同时到达目标的协同飞行方案。本发明工作流程示意图如图1所示，其中，在解析再入制导方法设计时，控制终端速度和抵达时间均由纵向升阻比控制，因此在控制上是一个欠驱动问题，本文利用纵程解析解和时间解析解设计了一种能够兼顾能量管理要求和抵达时间要求的纵向参考剖面规划方法，可以有效解决此问题。

整个过程包括以下几个步骤：

步骤1：再入制导问题描述

这里采用圆球形旋转地球下的三维再入质点动力学模型，其中利用经度λ、纬度φ和高度h描述飞行器位置信息，而用速率v、弹道倾角γ和航向角ψ描述速度信息。

极高的飞行速度导致力热环境极其恶劣。因此，为了保证飞行器的各个子系统正常工作，再入飞行轨迹需要满足热流密度来流动压q、过载n约束。考虑到飞控系统能力有限，需要对攻角和倾侧角范围及变化率进行限制。除了常规约束，这里进一步考虑多个圆形禁飞区约束。

再入段飞行在飞行器到目标的水平距离为staem时终止。此时期望的终端高度为htaem，终端速度为vtaem，终端航向误差|δψtaem|和终端倾侧角|σtaem|满足约束条件。此外，这里还指定期望的终端时刻为ttaem。

步骤2：基于旋转地球模型的飞行时间解析解

为了应对抵达时间约束，这里推导飞行时间解析解。能量定义如下

其中，v是速度，h是高度，re是地球平均半径，μ是引力常数。能量对时间t的导数为

式中g为重力加速度。在旋转地球背景下，有如下复杂的非线性方程

其中，γ是弹道倾角，d是阻力，m是质量，ωe是地球自转角速度，φ是纬度，ψ是航向角。公式(4)中右边第3项第4项是由地球自转引起的惯性力分量。本文将沿速度方向的惯性力分量视为为附加阻力，如下

定义等效阻力为将公式(33)～(35)代入公式(32)，并取其倒数，可得

定义纵向升阻比为其中l1＝lcos(σ)是升力在纵向平面内的分量，并规划纵向升阻比如下，

其中a0，a1，a2为二次多项式系数。

现在利用剖面估计平稳滑翔状态下的等效阻力(记为)。弹道倾角时间导数如下

假设γ≈0和可估算平稳滑翔纵向升力l1(sg)，如下

其中，δl1是惯性力沿垂向的分量，被视为附加纵向升力，公式如下

与此同时，采用上述假设，δd也可以被简化为

则平稳滑翔状态下的等效阻力可以由下式估算

将公式(42)代入公式(36)，并利用e替换v可得

由于h＜＜re，故可令h取滑翔高度的中间值h^*，进而令平均地心距r^*＝re+h^*。尽管惯性力是小量，但是当速度接近第一宇宙速度时，惯性力有可能会造成上式分母为零，从而发生奇异。为了避免发生上述奇异情况，这里将惯性力看作小量，并进行一阶taylor展开，如下

其中，简化符号hz1、hm的定义如下

hm＝2e+μ/r^*(46)

根据δl1和δd的表达式，可将hz1分为两个部分hpaf1和hpaf2，如下

hz1＝hpaf1+hpaf2(47)

其中，

hpaf1＝-2r^*vωecos(φ)sin(ψ)(48)

由于hpaf2的影响要远小于hpaf1，因此，可以粗略地采用线性函数对hpaf2进行近似，如下

hpaf2(e)≈kpaf2(1)v+kpaf2(0)(50)

其中，系数kpaf2(1)和kpaf2(0)由两个端点确定，如下

其中，hpaf2(e)由当前状态代入公式(49)计算得到，而hpaf2(etaem)是将期望的终端状态代入公式(49)计算得到。

实际上，高超声速滑翔飞行器围绕记为广义赤道的某一个大圆飞向目标，飞行器围绕广义赤道左右横向机动，广义赤道的方位角可以看作是飞行器航向角的某种加权平均值，因此，在后续航向角曲线尚且未知的情况下，本发明利用广义赤道的方位角替代hpaf2(e)和hpaf2(etaem)公式中的航向角。至此，在公式(44)中，除了自变量e之外，其它参量均为常值，进而可以对公式(44)进行积分。经推导可得飞行时间解析解的表达式如下

其中，系数kt(1)、kt(2)、kt(3)、kt(4)、kt(5)、kt(6)、kt(7)的表达式如下

本发明的飞行时间解析解与传统方法及数值弹道仿真的结果对比分析如实施例一所示。

步骤3：考虑多禁飞区和抵达时间约束的解析再入制导方法

此步针对单个飞行器，设计考虑多禁飞区和抵达时间约束的再入制导方法。根据弹道特点，再入过程一股分为三个阶段：下降段、平稳滑翔阶段和高度调整阶段。下降段的初始高度较高，大气稀薄，飞行器快速掉高。在平稳滑翔阶段，大气密度适中，升力足以平衡重力，滑翔高度平缓下降。得益于高升力和高速度，飞行器在此阶段的滑翔距离可达上万公里。当飞行器距离目标足够近时，进入高度调整阶段。此时飞行器通过适当调节攻角来获得期望飞行高度。相比于之前的制导方法(文献一)，这里主要有两个方面的改进：(1)根据能量管理要求和抵达时间要求，利用飞行时间解析解和纵程解析解改进了参考剖面的规划方法；(2)为了对高度调整阶段的轨迹进行精确微调，根据终端时间、速度、高度约束，设计了数值迭代规划方案以及相应的跟踪策略。下面步骤4～6详细介绍这三个阶段的制导方案。

步骤4：下降段制导方法设计

为了避免掉进稠密大气层而引起热流密度过大，下降段飞行器以最大可用攻角、零倾侧角下滑。当升力足以支撑飞行器平稳滑翔时，攻角平滑过渡到基准攻角，随即进入平稳滑翔阶段。

步骤5：平稳滑翔阶段制导方法设计

平稳滑翔阶段是最长、最重要、也是最复杂的再入飞行阶段，此阶段的制导方案利用三维再入弹道解析解和飞行时间解析解在线快速规划满足禁飞区和抵达时间约束的参考弹道。平衡滑翔阶段的制导流程如下：

s51设计一条基准攻角剖面，并确定相应的基准升阻比剖面；

s52考虑能量管理要求和抵达时间约束，合理设计参数化的纵向升阻比剖面，同时可以确定参数化的横向升阻比剖面；

s53为了补偿地球自转的影响，在和剖面的基础上，建立等效纵、横向升阻比(和)剖面；

s54按照能量管理和抵达时间要求，利用文献一中的纵程解析解和所述步骤2得到的飞行时间解析解求解纵向剖面参数；

s55每隔60s，利用纵程解析解以及文献一中的横程解析解根据禁飞区和终端位置约束更新一次倾侧反转序列。

s56调节基准倾侧角跟踪等效纵向升阻比剖面，并按计划倾侧反转；

s57为了抑制弹道振荡，将弹道阻尼反馈引入到基准攻角和倾侧角中，从而得到制导指令。

s58按照过程约束限制倾侧角指令；

s59从步骤s52开始重复上述过程，直到平稳滑翔阶段结束。

(1)基准攻角和基准升阻比

为了发挥飞行器最大能力，采取如下基准攻角αbsl剖面

其中，eα＝-5.55×10⁷j/kg位于平稳滑翔和高度调整阶段的交班点附近。etaem是再入段期望的终端能量。α1＝10°为升阻比最大时的攻角值。在高度调整阶段，为了获得期望的终端高度，将基准攻角逐渐减小至攻角α2。这里，α2被初步设置为6°，并在进入高度调整阶段前会对其进行精确微调。为了在大干扰情况下维持滑翔，人为限制αbsl≥5°。与基准攻角αbsl对应的基准升阻比剖面如图2所示。

(2)改进参考剖面

由于飞行能量管理和抵达时间控制相互耦合、相互制约，仅利用纵向升阻比剖面控制终端能量和时间并不容易。经大量探索，本发明发现如下形式的纵向升阻比剖面可以在保证满足终端能量要求的前提下，实现对抵达时间的有效调节。

其中，e代表当前时刻的能量，xe代表任意时刻的能量。函数fl1/d(1)(xe)、fl1/d(2)(xe)、fl1/d(3)(xe)的定义如下

和是三个能量节点对应的纵向升阻比，也是剖面的三个参数，其中，前两个参数用于管理能量并控制飞行时间，在下文中会利用解析解进行求解。为了使得倾侧角最终收敛到零，令第三个参数其中，是期望终端状态对应的理想升阻比，系数用于粗略补偿气动拉偏的影响，是利用气动辨识技术获得的当前实际升阻比，而是当前状态对应的理想升阻比。在下文中，当不特指某一能量节点对应的函数值时，为方便起见省去自变量的书写，如简写为

在确定基准升阻比剖面和纵向升阻比剖面之后，可以由几何关系确定参数化的横向升阻比剖面如下

其中，nr是已经发生的反转次数，sgn是符号函数，用于确定飞行器的初始倾侧方向，是的模值，由下式估算

为了补偿地球自转的影响，将惯性力与气动力组合为等效气动力。通过合理分析惯性力的变化规律，在上述和剖面的基础上，可以建立如下反比例函数组合形式的等效纵向升阻比剖面横向升阻比剖面

这里，

其中，hm如公式(46)所示，和是一些复杂的惯性力组合项的近似公式，如下

近似公式的系数kh1(0)、kh1(1)、kh2(0)、kh2(1)、kh3(0)、kh3(1)、kh3(2)、kh4(0)、kh4(1)由当前状态和期望的终端状态确定，计算方法如下。

1)计算系数kh1(0)和kh1(1)

其中

是上一制导周期获得的当前纵向升阻比，是期望的终端升阻比，φtaem是目标的纬度，是飞行器的广义赤道方位角航向角，是飞行器的广义赤道方位角航向角在终点处的值。

2)计算系数kh2(0)和kh2(1)

其中，

3)计算系数kh3(0)、kh3(1)和kh3(2)

4)计算系数kh4(0)和kh4(1)

其中

(3)纵向剖面参数解算

这里，本发明根据能量管理和抵达时间要求，利用纵程解析解和飞行时间解析解来确定公式(62)中的参数和

记a为一股形式的纵向升阻比多项式，即公式(37)的系数数组，如下

则利用文献一中的纵程解析解公式，可得在此纵向升阻比控制下，飞行器从能量节点e0到能量节点e的纵程为

其中，函数fxd(1)(e，e0，a)，fxd(2)(e，e0，a)和fxd(3)(e，e0，a)的定义如下

记a1、a2、a3和a4分别为与公式(62)-(63)相关的多项式系数数组，如下

则针对公式(62)所示的纵向升阻比剖面，可预测终端纵程为

同时，针对公式(62)所示的纵向升阻比剖面，利用公式(53)可预测终端时间如下

则公式(62)中的参数和可以通过求解如下二元一次方程组获得

其中，sgo是剩余飞行距离，ttaem期望的抵达时间。

但一方面受限于能量管理要求的制约，另一方面未能充分考虑横向机动的影响，方程组(90)的解会引起倾侧角大幅振荡，妨碍禁飞区规避，不利于实际应用。因此，本发明提出通过大幅调整纵向升阻比剖面来克服由横向机动引起的时间偏差。经探索，对纵向升阻比剖面参数进行修正，令

其中，和是理想情况下标称轨迹对应的纵向剖面参数，和则是根据实际情况需要进行修正的参数。将公式(91)代入(88)和(89)，并整理可得

其中，xd(bsl)(etaem，e)和t(etaem，e)是标称轨迹对应的部分，而δxd(etaem，e)和δt(etaem，e)则是对实际情况的修正，公式如下

标称情况下，与应使飞行器的纵程满足终端约束，即

xd(bsl)(etaem，e)＝sgo-staem(97)

则令将公式(93)代入公式(97)可解得

注意，上式分母中的表达式利用了如下的积分性质

另一方面，和用于消除终端时间误差δtf，但不能引起附加航程。利用公式(94)和(96)可得

解得

其中，简记符号d0的表达如下

为了预测终端时间误差δtf，需要事先离线规划标称弹道，并获得两个参考剖面，即射程剖面sgo(ref)(e)和时间剖面t(ref)(e)，其中，标称弹道的纵向升阻比剖面仅根据能量管理要求确定，而终端时间则由发射时间控制。但实际情况的剖面可能是sgo(e)和t(e)。因此预测终端时间误差δtf的第一个组成部分是δtf(1)＝t(e)-t(ref)(e)，第二个组成部分则是由航程误差δsgo＝sgo(e)-sgo(ref)(e)引起的δtf(2)。现在推导在基准控制下由δsgo引起的时间误差δtf(2)。将公式(98)代入公式(95)可得

利用差分公式可得由δsgo引起的时间误差δtf(2)

至此，完成了纵向剖面参数解算工作。

(4)倾侧反转序列规划

由于高超声速滑翔飞行器主要通过调整倾侧角模值跟踪纵向参考剖面，所以在确定纵向升阻比剖面之后，倾侧角模值剖面基本确定。因此，飞行器仅能通过适时改变倾侧方向控制横向机动，以便规避禁飞区域并最终抵达目标。本发明采用解析迭代方案规划倾侧反转序列。由于优先调整已有反转点应对禁飞区约束，此策略所需反转次数较少，可充分发挥飞行器的横向机动能力，应对大量禁飞区约束。下面介绍策略流程。

1)为了抗干扰，人为扩大所有禁飞区半径，从而使得飞行器与真实禁飞区边界之间保留一个安全距离。考虑到弹道误差的累积效应，这里将半径增量设置为关于飞行器到禁飞区距离的线性函数，并且会随着飞行器接近目标而逐渐减小；

2)仅根据终端位置约束，利用横程解析解规划两个反转点；

3)从最近禁飞区开始，利用解析解搜寻距离禁飞区最近的航迹点，以判断该禁飞区约束是否得到满足。如果是，则继续判断下一个禁飞区，否则，执行4)；

4)判断该禁飞区是位于倒数第二个反转点之前、最后两个反转点之间、还是最后一次反转点之后。如果是位于倒数第二个反转点之前，则执行5)；如果是位于最后两个反转点之间，则执行6)；如果是位于最后一次反转点之后，则执行7)；

5)如果此禁飞区位于倒数第二次反转点之前，则表明此禁飞区距离终点较远，因此可以不考虑终端位置约束。这里先利用禁飞区规避策略调整反转点，应对当前的禁飞区约束，然后利用终端横程控制策略1调整后续反转点，以保证满足终端位置要求。跳转执行8)。下文中会对禁飞区规避策略、终端横程控制策略1以及6)中的终端横程控制策略2进行说明；

6)如果禁飞区位于最后两次反转点之间，此时规划倾侧反转序列需要兼顾禁飞区约束和终端位置要求。这里同样首先利用禁飞区规避策略应对禁飞区约束，但是由于此时距离目标较近，之后会采用终端横程控制策略2规划后续反转点，以满足终端横程要求。跳转执行8)；

7)如果禁飞区位于最后一次反转点之后，由于禁飞区距离终点太近，这里忽略此禁飞区约束，并保留原有反转点序列。跳转执行8)；

8)确定下一个禁飞区，回到3)继续判断禁飞区约束是否得到满足，直到检测完所有禁飞区约束。

下面说明禁飞区规避策略、终端横程控制策略1和终端横程控制策略2的原理。

1)禁飞区规避策略

此策略分两种情况应对禁飞区约束：当前禁飞区位于倒数第二次反转点之前；当前禁飞区约束位于最后两次反转点之间。

在第一种情况中，由于距离目标较远，在处理当前禁飞区约束的过程中无需考虑终端位置约束。本发明利用解析解优先调整已有反转点，以使得飞行器从靠近禁飞区的一侧绕过禁飞区。如果不成功，则进一步尝试从禁飞区的另一侧绕过去，或者新增一个反转点应对禁飞区约束。在成功规避当前禁飞区之后，由于反转序列调整会改变之前的轨迹，所以需要检查之前的禁飞区约束是否仍然满足。如果不满足，则进一步判断哪个威胁等级更高，以决定是否采纳新反转点序列。

在第二种情况中，由于距离目标较近，此时反转点序列规划需要兼顾禁飞区规避和终端位置约束。本发明优先调整已有反转点应对禁飞区约束，但不同的是将尽量使得飞行器沿禁飞区靠近目标的一侧绕过禁飞区。后续过程与上面相同。

2)终端横程控制策略1

此策略主要是针对当前禁飞区位于倒数第二次反转点之前的情况。在利用禁飞区规避策略成果应对当前禁飞区约束之后，由于轨迹调整，后续反转点序列并不满足终端位置约束。因此，终端横程控制策略1将放弃后续的反转点，并根据终端位置约束，重新规划最后两次反转。

3)终端横程控制策略2

此策略主要是针对当前禁飞区约束位于最后两次反转点之间的情况。同样，在成功调用禁飞区规避策略之后，放弃后续的反转点。由于此时距离目标较近，这里首先调整最后一次反转点应对终端位置约束，然后判断该反转点是否满足平稳滑翔阶段和高度调整阶段的交班条件。如果不满足，则通过分析横向机动规律适当增加1-2次反转。

(5)平稳滑翔阶段的基准倾侧角

利用纵向升阻比与总升阻比之间的几何关系，可以确定基准倾侧角。记第nr次倾侧反转对应的能量为ebr(nr)。当ebr(nr)+δe≥e≥ebr(nr+1)+δe时，基准倾侧角σbsl为

其中，补偿值δe使得倾侧反转适当提前，从而补偿由滚转速率限制引起的响应滞后。

(6)平稳滑翔阶段的指令攻角和倾侧角

为了抑制弹道振荡，引入三维弹道阻尼控制技术，如下

αcmd＝αbsl+cos(σbsl)kγ(γsg-γ)(106)

其中，αcmd和σcmd是指令攻角和倾侧角，kγ是反馈增益系数，γsg是平稳滑翔弹道倾角，可由当前状态确定。

为了满足热流密度、来流动压和过载等约束，可将过程约束转化为倾侧角约束，进而限制倾侧角指令大小。

步骤6：高度调整阶段制导方法设计

在最后一次反转点(记此节点对应的飞行器能量为ebr(ntr))之后，飞行器进入高度调整阶段，但是在最后一次反转点之前的某个节点(记为对应的飞行器能量为est)，飞行器就开始了高度调整阶段的准备工作。这里本发明设计了一种基于在线弹道仿真的多目标数值迭代规划方案。当e＝est时，利用多目标数值迭代规划方案获得纵向升阻比、剩余飞行距离、以及飞行时间相对能量的参考剖面，并精确微调后续轨迹，以满足终端时间、速度、高度要求。由于仅需要对后续一小段弹道进行数次仿真，所以并不会造成较大的计算负担。当est＞e＞ebr(ntr)时，本制导方法跟踪在e＝est时由多目标数值迭代规划方案获得的纵向升阻比、剩余飞行距离、以及飞行时间相对能量的参考剖面，而当ebr(ntr)＞e时，则利用比例导引律确定基准倾侧角，以消除航向误差，并通过微调攻角跟踪飞行时间和剩余飞行距离参考剖面，以保证满足终端时间和速度要求。高度调整阶段的具体制导方案如下：

(1)多目标数值迭代规划方案

多目标数值迭代规划方案考虑了终端时间、速度和高度约束，其中终端时间和速度由参数和ebr(ntr)控制，终端高度由基准攻角参数α2控制。如图3所示，多目标数值迭代规划方案流程如下：

s611、利用弹道仿真预测终端时间和终端状态；

s612、判断终端速度或终端时间是否满足要求；若是则跳转步骤s613，否则跳转步骤s614；

s613、根据终端速度和终端时间要求，利用拟牛顿迭代法调整纵向升阻比剖面和最后一次反转点，跳转步骤s611；

s614、判断终端高度是否满足要求；若是，则跳转步骤s616，否则跳转步骤s615；

s615、根据终端高度误差修正基准攻角参数，跳转步骤s611；

s616、完成迭代规划。

在步骤s613中，针对不同情况，本发明的应对策略如下：

1)如果经弹道仿真预测发现终端速度和终端时间大于期望值，则可以适当增大减小和ebr(ntr)；

2)如果终端速度大于期望值而终端时间小于期望值，则可以适当减小和ebr(ntr)；

3)如果终端速度小于期望值而终端时间大于期望值，则可以适当增大和ebr(ntr)；

4)如果终端速度和终端时间均小于期望值，则可以适当减小增大和ebr(ntr)。

由于与ebr(ntr)同步变化，可将作为一个辅助调节手段，并设计如下关系式

其中，和δebr(ntr)分别表示与ebr(ntr)的变化量，kf＝2×10^-7为调节系数。上式表示若ebr(ntr)改变0.1×10⁷kj/kg，则相应改变0.2。这样减小了一个搜索维度，有助于降低算法复杂度。

下面利用牛顿迭代法搜索满足终端时间、速度和高度约束的和ebr(ntr)。定义向量x和向量函数g

其中，和分别为参数和ebr(ntr)对应的在线弹道仿真得到的终端速度和时间。求解方程g(x)＝0的牛顿迭代公式为

其中，上标“(k)”代表第k次迭代，矩阵jg是向量函数g的雅克比矩阵，如公式(81)所示，而是jg的逆矩阵。

由于和需要通过积分获得，jg无法解析求解。如果利用常规的差分方法求解，则在每一步需要执行3次弹道积分，这会造成较大的计算负担。

为了减轻计算负担，本发明提出一种基于方向导数的jg求解方案。定义向量p^(k)＝x^(k)-x^(k-1)，定义向量q^(k)是垂直于p^(k)、模长为0.01的一个小向量。则有如下方向导数的近似计算公式

其中，||p^(k)||和||q^(k)||分别表示向量p^(k)和q^(k)的模长。利用方向导数与偏导数之间的关系可得

其中，θx(p^(k))和θy(p^(k))是向量p^(k)分别与x轴、y轴的夹角，而θx(q^(k))和θy(q^(k))则是向量q^(k)分别与x轴、y轴的夹角。利用公式(113)可解得

在当前迭代中，由于只需要计算g(x^(k)+q^(k))和g(x^(k))，这里仅需执行2次弹道仿真，从而大幅降低计算量。由于每一步迭代时，的变化量的数量级大约为0.1～1，而δebr(ntr)的数量级为10⁶，上述方程通常是病态的，但可以利用换元法解决此问题。

在步骤s615中则利用以下公式(85)微调基准攻角剖面参数α2。由于公式精度较高，这里仅需要修正一次。α2需要的改变量δα2为

其中，clf(est)是通过气动辨识技术获得的终端升力系数，是估计的升力线斜率，sref是飞行器参考面积，ef、qf、vf、hf、φf、ψf分别是在线弹道仿真后得到的终端能量、动压、速度、高度、纬度、航向角的值，qtaem、vtaem分别是期望的终端动压与速度值。

(2)高度调整阶段的指令攻角和倾侧角

当est＞e＞ebr(ntr)时，飞行器处于高度调整阶段之前的准备阶段。此时制导方案跟踪由多目标数值迭代规划方案获得的纵向升阻比剖面、剩余飞行距离剖面、以及飞行时间剖面，分别记为sgo(ref)(e)、tref(e)。基准倾侧角σbsl的计算公式如下

其中，这里纵向升阻比主要由组成，但是为了在有干扰情况下保证能量管理和终端时间要求，设计并跟踪考虑时间误差修正的剩余飞行距离剖面，如下

sgo(ref)(2)(e)＝sgo(ref)(e)-ksgov(t-tref(e))(117)

其中，ksgo＝0.2sgo/sgo(st)。sgo(st)是执行多目标数值迭代规划方案时刻对应的剩余飞行距离。

此后，将公式(116)代入公式(107)，并结合公式(106)，就可指令攻角和倾侧角。

当e＜ebr(ntr)时，飞行器进入高度调整阶段，基准攻角依然如公式(61)所示，而基准倾侧角则有比例导引律确定。视线方位角变化率为

其中，δψ是航向误差。比例导引律产生的需用横向机动加速度为

其中，δl2是惯性力沿l2方向的分量。上式右边第二项-δl2/m是为了补偿地球自转的影响。为了防止初始倾侧角饱和，这里令有效导引比kpn随剩余飞行距离从2逐渐变化到4。另一方面，在平稳滑翔情况下，纵向平面内升力加速度al1与重力、离心力和惯性力平衡，即

则基准倾侧角为

进而，设计攻角和倾侧角指令如下

αcmd＝αbsl+kα[sgo-sgo(ref)(2)(e)](122)

步骤7：多飞行器抵达时间协同方案设计

步骤3-6完成了针对单个飞行器的、考虑多禁飞区和终端时间约束的解析制导方法设计，这里进一步针对多个飞行器研究抵达时间协同方案。一股情况下，飞行器抵达目的地的时间受两个重要因素影响：发射时间以及再入飞行过程中对时间的控制，其中前者用于对抵达时间进行大范围调整，而后者用于克服干扰，对时间进行精确微调。但是，由于本发明并不考虑助推段和下压段的飞行过程，本发明中利用再入段的起始时间替代发射时间，利用再入段的终端时间替代抵达时间。注意：在实际情况中需要综合考虑这三段的飞行过程。尽管助推过程中的不确定因素会引起再入段起始时刻和起始状态的偏差，但是利用前面介绍的制导方法可以在再入飞行中有效克服这些偏差因素。

常见的时间协同问题有：同步到达同一目标，异步到达同一目标，同步到达不同目标，异步到达不同目标，以及它们之间的组合问题。由于这些问题的应对策略类似，下面以同步到达同一目标的问题为例介绍协同方案。

假设有nhgv个飞行器，分别记为hgvi，i＝1，2，…，nhgv。在给定发射点和目标点之后，根据助推火箭性能和所选主动段制导方案，可以确定飞行器的再入段起始状态。之后，同公式(98)，令纵向升阻比剖面参数仅根据能量管理要求，利用纵程解析解确定相应的纵向升阻比剖面参数，如公式(94)所示。这里，下标“hgv(i)”代指第i个飞行器。

进而，利用离线弹道仿真可以预测每个飞行器的飞行时间tef(hgvi)，i＝1，2，…，nhgv，从中选出时间最长的一个，并预留发射准备和助推飞行时间，可以确定期望的抵达时间ttaem，进而利用ttaem-tef(hgvi)可以确定各个飞行器的再入段起始时间。之后，根据所选的助推段制导方法，可以确定相应的发射时间。

实施例：

实施例一

本实施例进行本发明的飞行时间解析解与传统方法及数值弹道仿真结果对比，验证本发明的飞行时间解析解的预测精度。

文献二中常值纵向升阻比条件下的飞行时间公式为

设置飞行器的初始经度λ0＝0deg、初始纬度φ0＝50deg，初始能量e0＝-3.8602×10⁴kj/kg和终端能量ef＝-5.5×10⁴kj/kg。考虑5个不同的飞行方向：ψ0＝100deg、180deg、-100deg、20deg、-20deg。设置文献二中的纵向升阻比为横向升阻比为

仿真结果如表1所示。从仿真结果中可以看出，受地球自转的影响，飞行器沿不同方向的飞行时间并不相同，但是文献二的解析解未能考虑由地球自转引起的时间差异。通过与弹道仿真结果进行对比可以看出：由于本发明解析解对地球自转的影响进行了补偿，其精度较高，预测误差保持在3％以内。在计算效率上，时间解析解的计算耗时至少比弹道仿真小5个数量级。

表1

实施例二

本实施例为理想无干扰的情况，解决在多禁飞区环境中导引多个高超声速滑翔飞行器(v1，v2，v3)从不同地点发射但同时到达同一目标的协同飞行问题。在飞行区域内布置了64个半径为200km的圆形禁飞区，飞行器的初始条件如表1所示，目标点位置为经度λt＝130deg、纬度φt＝-20deg。利用步骤7的多飞行器抵达时间协同方案，可选取ttaem＝2900s，并确定各个飞行器的再入起始时刻，分别为109.29s、358.69s和751.62s。

表2

再入段飞行在飞行器到目标的水平距离为staem＝50km时终止。此时期望的终端高度为htaem＝25km，终端速度为vtaem＝2000m/s，终端航向误差满足|δψtaem|≤5deg和终端倾侧角满足|σtaem|≤30deg。

图4展示了飞行器在本发明制导方法控制下成功规避了所有禁飞区，并抵达目标。图5是高度-速度曲线，其中hmin曲线是由过程约束确定的高度下界。从图中可以看出，由于在高速阶段力热环境比较恶劣，飞行器容易贴近高度下界飞行。此外还可以看出终端速度和高度满足终端要求。图6展示了飞行器的攻角规律，从中可以看出三个飞行器具有不同的再入起始时刻。图7展示了倾侧角剖面。图8展示了三种情况对应的参考剩余飞行距离剖面。图9展示了三种情况对应的参考飞行时间剖面。

实施例三

本实施例通过montecarlo仿真验证本发明制导方法在飞行器模型有拉偏情况下的鲁棒性。其中，气动系数采用如下线性拉偏模型

其中，δcl和δcd分别是升力系数和阻力系数拉偏百分比，随攻角α和马赫数ma变化。δcl0、是相关的零均值正态分布随机参数。风力和大气密度摄动模型如下

δρ＝kρδρ(max)(130)

其中，是沿东西方向的风速，是沿南北方向的风速，δρ是大气密度拉偏百分比。vwind(max)是最大可能风速，随高度变化，在高空可达170m/s。δρ(max)是大气密度最大可能拉偏百分比，在高空可达50％。和kρ是相应的零均值正态分布随机参数。上述随机参数及初始条件误差的统计特性如表3所示。

表3

这里考虑从两个不同地点发射两个飞行器(v1，v2)进行协同飞行的情景，标称条件如表4所示，目标位置为经度λt＝120deg和纬度φt＝0deg。禁飞区分布如图10所示。设定期望抵达时间为2900s，则利用多飞行器抵达时间协同方案可已确定标称情况下再入段起始时间为426.68s和404.29s。注意，由于是事先离线计算，标称轨迹以及制导律所需要的剩余飞行距离、飞行时间参考剖面是基于理想气动模型和理想大气模型获得的，并未考虑、也无法考虑任何在线测量到的拉偏数据。

表4

图10展示了飞行器的再入轨迹，从中可见，在存在不确定干扰环境中，本制导方法能够成功导引飞行器规避禁飞区，并顺利抵达目标。图11展示了高度-速度曲线，其中hmin是由热流密度等过程约束确定的高度下界。由图可见，本制导方法可保证飞行安全。图12-13展示了攻角曲线和倾侧角曲线。图14展示了终端速度和终端航向误差的分布情况，均在要求范围之内。图15统计了到达时间分布情况，飞行器v1的抵达时间误差在±5s以内的占比为79.01％；飞行器v2的抵达时间误差在±5s以内的占比为74.36％。