1.一种考虑多个禁飞区约束的解析再入制导方法，其特征在于：它包括以下几个步骤：
步骤一：高超声速飞行器建模；
在旋转地球背景下，采用标准大气模型和平方反比引力场模型，建立高超声速飞行器的六自由度动力学方程；其中，利用经度λ、纬度φ和高度H来描述飞行器位置信息，而用速率V、弹道倾角γ和航向角ψ来描述速度信息；
飞行器在飞行过程中要受到约束；
1.1过程约束：
在飞行过程中再入弹道需要满足热流密度 来流动压q、过载n的约束，以保证飞行安全；此外，考虑到飞控系统能力有限，需要对攻角和倾侧角范围及变化率进行限制；除了上述约束，还要考虑多个圆形禁飞区约束；记禁飞区的个数为nNFZ，分别编号为
圆心为 半径为 记Cpi的经度为λCpi，纬度为φCpi；
1.2终端条件：
当飞行器到目标距离为STAEM＝50km时，再入段终止；此时期望的终端高度为HTAEM＝
25km，终端速度为VTAEM＝2000m/s，终端航向误差满足|ΔψTAEM|≤5deg和终端倾侧角满足|σTAEM|≤30deg；
步骤二：再入制导方法设计；
根据弹道特点将再入过程划分为三个阶段：下降段、平稳滑翔阶段和高度调整阶段；
下降段：
以最大可用攻角、零倾侧角下滑；当升力足以支撑飞行器平稳滑翔时，攻角平滑过渡到基准攻角，随即进入下一个阶段；
平稳滑翔阶段：
平稳滑翔阶段是最长、最重要、也是最复杂的再入飞行阶段；由于需要利用改进的三维再入弹道解析解在线快速规划满足多个禁飞区约束的参考弹道，下面介绍三维再入弹道解析解；
三维再入弹道解析解：
通过对不考虑地球自转的再入运动模型进行适当的简化和线化，得到了一个线性时变LTV系统，然后通过提出基于矩阵谱分解的LTV系统解析法获得了如下通用形式的、可以预测纵程xD、横程xC和航向误差Δψ的三维再入弹道解析解：
其中，xD代表纵向射程，Re代表地球半径，μ代表海平面处重力加速度与地球半径的平方的乘积，是一个常数，xC代表横向射程，xC0代表初始时刻的横向射程，Δψ代表横程误差，Δψ0代表初始时刻的横程误差；f1(xE)是自定义的函数，H代表飞行器距海平面高度，V是飞行器相对于地球的速率；E代表当前能量，E0代表初始能量，xE代表任意时刻的能量；R*为常数，大小为Re+H*，其中H*取滑翔高度平均值；L1和L2分别是升力在纵平面和当地水平面内的分量，即L1＝Lcos(σ)和L2＝Lsin(σ)；定义纵向升阻比L1/D为L1与阻力D的比值，而横向升阻比L2/D为L2与D的比值；从解析解中看出，L1/D和L2/D分别对纵程和横程有重要影响，并且作为控制量被设计成关于能量的函数；
为了补偿地球自转的影响，将由地球自转引起的惯性力与气动力组合为等效气动力，然后采用等效纵、横向升阻比作为参考剖面；通过合理预测惯性力的趋势，将等效纵、横向升阻比剖面设计成反比例函数；在设计过程中，采用了如下公式计算平稳滑翔阻力DSG：
其中，ΔL1是惯性力沿L1方向的分量；进而采用线化处理方法得到如下改进的气动关系；m表示飞行器质量，g表示重力加速度；
其中，h1，h2，h3，h4为自定义参数，具体表达式如下；
hz1，hz2，hz3，hz4，hm为自定义参数，具体表达式如下；
h4则是与ΔL1相关的修正项；相关的分子和分母表达式如下：
其中，ΔL2是惯性力沿L2方向的分量，ΔD是惯性力沿D方向的分量；ΔL1、ΔL2和ΔD的公式如下：
其中，ψ是飞行器航向角，以当地北向为基准；φ是纬度；ωe为地球自转角速度；
为了便于应用解析解，需要获得h1、h2、h3和h4与能量的关系；由于hm已经是关于能量的函数，不需要对其进行处理；对hz1、hz2、hz3和hz4进行如下近似处理；
利用当前状态和期望的终端状态确定上述多项式系数kh1(0)、kh1(1)、kh2(0)、kh2(1)、kh3(0)、kh3(1)、kh3(2)、kh4(0)和kh4(1)；
将 和 代入公式(1)-(3)，然后积分公式(1)-(3)，得到考虑地球自转补偿的三
维再入弹道解析解；不失一般性，令：
a0，a1，a2，b0，b1，b2为多项式系数；
则改进的解析解如下：
其中，fxD(1)，fxD(2)，fxD(3)，f2(xE)，f3(xE)为自定义函数，具体表达式如下平稳滑翔阶段制导流程：
步骤2.1：设计一条基准攻角剖面，并确定与之对应的基准升阻比剖面；
步骤2.2：参照基准升阻比剖面的外形，设计纵向升阻比剖面；
步骤2.3：在L1/D和L2/D的基础上，进一步获得关于 和 的参数化剖面；
步骤2.4：通过利用新纵程解析解求解方程xD(ETAEM,E)＝sgo-STAEM来计算参数L1/D1；
其中，ETAEM表示终端能量；sgo表示剩余射程；STAEM表示终端射程；
步骤2.5：获得考虑地球自转补偿的横程解析解；进一步考虑多个禁飞区约束，然后利用横程解析解每隔60s规划一轮反转点；
步骤2.6：调节基准倾侧角跟踪等效纵向升阻比剖面，并实现倾侧反转；
步骤2.7：为了抑制弹道振荡，采用弹道阻尼反馈；
步骤2.8：根据过程约束限制倾侧角指令；
步骤2.9：从步骤2.2开始重复上述过程，直到平稳滑翔阶段结束，即当最后一次反转发生时；
高度调整阶段：
在最后一次反转发生之后，飞行器进入高度调整阶段；通过调节基准攻角来调整飞行高度，并利用比例导引律确定基准倾侧角，从而消除终端航向误差；在最后一次反转发生之前利用基于在线弹道仿真的迭代修正方法对最后一次反转点进行微调，以满足终端速度要求；
步骤三：考虑多禁飞区约束的倾侧反转点解析迭代规划方案；
借助广义经纬度坐标系统来描述再入运动；为此，下面首先确定禁飞区的广义经纬坐标；接着，为了判断轨迹是否满足禁飞区约束，需要利用解析解迭代解算距离禁飞区最近的航迹点；最后，设计一种反转点解析规划策略，根据禁飞区规避要求适当调整或增加反转点，直到获得一条满足禁飞区和终端位置约束，且所需反转次数少的参考弹道；
3.1禁飞区的广义经纬坐标
假设有nNFZ个圆形禁飞区，分别编号为 其圆心分别记为
半径为 记禁飞区pi的圆心经纬坐标为(λCpi,φCpi)；记从地心CE指向M、T和Cpi的单位向量为 和 由下式计算：
记 为垂直于广义赤道平面的单位向量，由下式计算：
过Cpi作一个与广义赤道平面相垂直的大圆，其与广义赤道的交点记为Bpi；定义 为从CE指向Bpi的单位向量；由于 与 在广义赤道平面内的分量同向，则 由下式计算：
则Cpi的广义经度 为 与 的夹角，由下式计算：
Cpi的广义纬度 为 与广义赤道平面的夹角，由下式计算：
计算当纵程等于 时的能量ECpi；ECpi在后面将被用作解析迭代求解距离禁飞区最近航迹点Tpi的初值；由于 是分段函数，纵程解析解如下：
其中，Eα表示基准攻角剖面分界点；E代表当前能量，Es代表之后某一时刻的能量；
xD对Es的偏导数为：
其中，h1(ECpi)和h4(ECpi)由公式(9)计算；利用Newton法求解 得到
ECpi，如下：
表示距离禁飞区最近航迹点Tpi的能量的第k次迭代值； 表示距离禁飞区最近航
迹点Tpi的能量的第k+1次迭代值；
3.2距离禁飞区最近航迹点
为了判断再入轨迹是否满足禁飞区约束，需要寻找轨迹上距离禁飞区中心最近的点Tpi，并计算其到禁飞区中心的距离；
为了解析迭代搜寻最近点，需要进一步给出考虑多次倾侧反转的横程解析解以及相关偏导数，如下：
如果E和Es之间没有反转点，则有：
其中，xD(Es,E)由公式(33)计算；f2(xE)的表达式如公式(26)所示；nR是已经反转的次数，也就是E之前的反转次数；sgn是符号函数，用于确定初始倾侧方向；F函数的定义如公式(37)所示；F函数的第三个输入与积分区间无关，而是一个与被积函数相关的参数；
通过比较公式(27)和(38)，可见f4(xE)与f3(xE)之间存在细微的不同；h2(xE)、h3(xE)和h4(xE)表达式如公式(9)所示；相比较于公式(21)，公式(36)不含与xC0相关的项；这是因为广义赤道在每一个制导周期根据当前状态更新一次，从而导致了xC0＝0；此时xC对Es的偏导数为：
如果E和Es之间有一个反转点
xC对Es的偏导数为：
如果E和Es之间有ns≥2个反转点：
表示反转nR+nS次时的能量， 表示反转nR+k+1次时的能量, 表示
反转nR+k次时的能量；
xC对Es的偏导数为：
定义一个与轨迹到禁飞区中心距离相关的量，如下：
RCpi对Es的偏导数为：
记Tpi对应的能量为ETpi；当Es＝ETpi时RCpi对应极小值，则有R′Cpi(ETpi)＝0；利用Secant法求解此方程就得到ETpi，如下：
表示第k-1次迭代的ETpi能量值；
迭代所需的两个初值在ECpi附近选取；
3.3渐增式反转点解析迭代规划方案
总体方案
设计反转次数渐增的解析迭代规划方案；将满足 要求的禁飞区按照由近及
远的顺序重新排列，其中， 是目标的广义经度，并记由此获得的有效禁飞区个数为nNFZ(valid)；在首次规划中仅根据横程要求安排两次反转，然后逐个检测所有有效禁飞区约束满足情况；如果存在禁飞区约束不满足情况，则进一步分三种情况处理此禁飞区约束：
情况一：如果此禁飞区位于倒数第二次反转点之前，即ETpi≥EBR(nTR-1)，则表明此禁飞区距离终点较远；此时，首先利用禁飞区规避策略调整已有反转点或增加新反转点，以规避禁飞区，然后采用终端横程控制策略1规划最后两次反转点，以满足终端横程要求；为了尽量减少反转次数，禁飞区规避策略优先调整已有反转点；
情况二：如果禁飞区位于最后两次反转点之间，即EBR(nTR-1)＞ETpi≥EBR(nTR)，则首先同样利用禁飞区规避策略来满足禁飞区约束；但是由于此时距离目标较近，之后会采用终端横程控制策略2规划最后1至3次反转点，从而满足终端横程要求；
情况三：如果禁飞区位于最后一次反转点之后，由于禁飞区距离终点太近，忽略此禁飞区约束，并保留原有反转点；
随后继续检测下一个禁飞区约束，直到完成所有有效禁飞区约束的检测，即当pi＞nNFZ(valid)时；在每处理完一个禁飞区约束之后就立即考虑终端横程约束，而不是在考虑完所有禁飞区约束之后再考虑终端横程约束；这样有利于防止轨迹为了规避禁飞区而偏离目标太远，从而导致最终无法抵达目标；
为了应对干扰，在参考弹道与禁飞区边界之间设置一个安全距离；记禁飞区pi对应的安全距离为Δrpi；在规划过程中，首先通过将禁飞区半径扩大Δrpi构造虚拟禁飞区，然后根据这些虚拟禁飞区约束进行弹道规划；考虑到误差累积效应，Δrpi被设计成关于飞行器到禁飞区距离的函数，并随着飞行器接近禁飞区而逐渐减小，如下：
其中，spi代表飞行器到禁飞区pi边界的距离，Δr*＝100km和s*＝10000km分别是两个参数，根据经验设定；对于相邻的两个虚拟禁飞区相交而对应的真实禁飞区不相交的情况，需要修正两个禁飞区对应的安全距离；
禁飞区规避策略：
下面详细介绍禁飞区规避策略；此策略采用三种方式规避禁飞区：方式一：新增一次反转点；方式二：调节已有的一个反转点；方式三：调节已有的两个反转点，其中，优先调节已有反转点来规避禁飞区；由于高超声速飞行器的转弯半径高达数千公里，极有可能出现无法规避某些禁飞区的情况；在此情况下，策略将会尽量规避威胁等级较高的禁飞区，但是允许飞行器进入威胁等级较低的禁飞区；
禁飞区规避方式一：
通过在禁飞区pi之前增加一次反转来规避禁飞区；记nR(pi)为原反转点序列中位于禁飞区pi之前的反转点次序，那么新增反转点的次序为nR(pi)+1；下面需要确定飞行器是沿禁飞区上边界绕飞，还是下边界绕飞；换而言之就是确定期望的规避横程 由飞行力学可知：如果倾侧角为正，则横程减小；反之，则横程增大；如果 则新增反转将
使得倾侧角变负，从而增大横程；因此，这时期望飞行器沿第pi个禁飞区上边界绕过去，即反之，当 时，这时期望飞行器沿禁飞区右侧边界绕过去
这里， 和 由下式计算：
代表期望的规避横程， 代表第pi个禁飞区上边界， 代表第pi
个禁飞区下边界；
其中， 是禁飞区pi中心对应的横程，rpi是禁飞区pi的半径；利用Newton法迭
代调节新增的第nR(pi)+1次反转点EBR(nRpi+1)，如下所示：
其中， 利用公式(40)计算，并且
禁飞区规避方式二：
用于E与ETpi之间至少存在一个反转点的情况，通过调整禁飞区之前的一次反转点EBR(nRpi)来规避禁飞区；当禁飞区pi距离目标较远，即ETpi≥EBR(nTR-1)时，如果xC(ETpi)与较近，则首先尝试沿禁飞区上边界绕飞，即令 反之，如果xC
(ETpi)与 较近时，则尝试令 当禁飞区pi距离目标较近，即
EBR(nTR-1)＞ETpi≥EBR(nTR)时，则需要兼顾目标抵达要求；此时，如果 则首先尝试从禁飞区下边界绕过去；如果 则尝试从禁飞区上边界绕过去；如果上述尝试不成功，则进一步尝试从禁飞区另一侧绕过去；
其中， 由公式(40)或(42)计算，并且
禁飞区规避方式三：
在调用禁飞区规避方式二调整失败的情况下，进一步尝试利用禁飞区规避方式一来同时调整禁飞区之前的两个反转点EBR(nRpi-1)和EBR(nRpi)，以规避禁飞区；
EBR(nRpi)和EBR(nRpi-1)表示针对禁飞区pi增加的反转点；
首先将EBR(nRpi)右移至ETpi；此时进一步将EBR(nRpi-1)左移来扩大正倾侧角区间；如果EBR(nRpi)已经被左移到了EBR(nRpi-1)，则无法通过调节EBR(nRpi-1)来扩大负倾侧角区间；因此，禁飞区规避方式三仅适用于需要将EBR(nRpi-1)左移的情况；进而确定
如果 则将EBR(nRpi-1)向左移动会扩大正倾侧角区间，从而使得xC(ETpi)减
小，因此令 反之则令 令EBR(nRpi)＝ETpi，并利用如
下公式调整EBR(nRpi-1)
其中， 由公式(40)计算，并且
终端横程控制策略：
在处理完禁飞区pi约束后，需要随即考虑终端横程要求；这里分两种情况考虑：(1)如果禁飞区pi位于原反转点序列里的倒数第二个反转点之前，则说明禁飞区pi距离终点还很远，此时采用终端横程控制策略1来消除终端横程误差；(2)如果禁飞区pi位于原反转点序列里的最后两个反转点之间，则需要采用终端横程控制策略2来消除终端横程误差；
首先介绍终端横程控制策略1，首先判断在处理禁飞区pi约束的过程中是否对原反转点序列进行了调整；如果没调整，则保留原反转点序列即可；如果调整了，则弃置禁飞区pi之后的反转点，然后重新安排两次反转点，并利用终端横程控制策略1进行规划；如果规划成功，则采用新反转点序列；否则，忽略禁飞区pi，并保留原反转点序列；
为了保留抗干扰的能力，这里安排两次反转；前面介绍过，记总反转次数为nTR，记最后两次反转对应的能量为EBR(nTR-1)和EBR(nTR)；这里根据高度调整阶段的需要，令EBR(nTR)＝Eα，然后利用secant法解析迭代规划EBR(nTR-1)，如下：
其中， 仍然由公式(42)计算，并且有
接下来介绍终端横程控制策略2，此策略适用于禁飞区pi位于原反转点序列里的最后两个反转点之间的情况；这里同样首先判断在处理禁飞区pi约束的过程中是否对原反转点序列进行了调整；如果没调整，则保留原反转点序列即可；如果调整了，则尝试利用终端横程控制策略2重新规划最后一次反转点，以满足终端位置要求；如果不成功，则忽略禁飞区pi，并保留原反转点序列；如果成功，则进一步观察调整后EBR(nTR)与Eα的大小关系；为了满足高度调整阶段制导要求，根据经验，希望|EBR(nTR)-Eα|＜0.3×107J/kg；如果不满足此条件，则进一步需要考虑如下两种情况：
情况一：如果EBR(nTR)-Eα≥0.3×107J/kg，则为了保持EBR(nTR)＝Eα，在ETpi之后增加两次反转点EBR(nTR-2)和EBR(nTR-1)，然后利用终端横程控制策略1调整EBR(nTR-1)；这样同样达到终端横程调节目标；
情况二：如果EBR(nTR)-Eα≤-0.3×107J/kg，则同样为了保持EBR(nTR)＝Eα，在ETpi之后增加一次反转点EBR(nTR-1)，然后利用终端横程控制策略1调整EBR(nTR-1)，以满足终端横程要求；
通过调节最后一次反转来消除终端横程误差，如下：
其中， 由公式(40)或(42)计算，并且有