一种考虑禁飞区约束的天基再入飞行器可达域快速计算方法

1.本发明涉及一种考虑禁飞区约束的天基再入飞行器可达域快速计算方法，属于天基再入飞行器技术 领域。


背景技术：

2.天基再入滑翔飞行器的可达区域指的是飞行器能够到达的着陆点的分布范围。可达区域在飞行器前 期的设计中能够对飞行器的飞行能力和机动范围进行分析，为飞行器的任务规划提供关键信息，另一方 面当飞行器在飞行过程中出现突发故障而不能维持标称状态飞行时能够为飞行器备选目标点提供重要 依据。可达区域的计算一方面在飞行器前期的设计中能够对飞行器的飞行能力和机动范围进行分析以及 为飞行器的任务规划提供关键信息，另一方面当飞行器在飞行过程中出现突发故障而不能维持标称状态 飞行时能够为飞行器备选着陆点提供重要依据。航天器可达区域的分析对完成各项任务具有重要参考价 值，同时也能为各种轨道优化设计提供相应的理论基础。
3.对一组初始条件，飞行器的可达区域是可行的目标点的集合。目标点是可行的是指飞行器有足够的 能力到达该目标点，并且在飞行过程中不违背任何约束条件。可达区域能为飞行任务规划提供重要的参 考信息，因此需要对可达区域进行计算，有助于在任务规划阶段得到更加合理的任务策略。传统的飞行 器可达域计算，通常忽略离轨段机动能力而仅考虑大气层内的飞行可达域。少有研究考虑地球形状为椭 球，且飞行器在不同的初始射向下可达域形状的变化规律，而是将不同初始射向下的可达域范围近似为 某一固定射向下的边界。考虑禁飞区对可达域的影响，有助于得到更加可行的任务规划结果，传统针对 助推-滑翔飞行器的可达域计算方法已经难以适应新任务的需求。因此，为了在任务规划阶段考虑禁飞 区的影响，提出一种考虑禁飞区约束的天基再入飞行器可达域快速计算方法，具有较强的工程实际意义。


技术实现要素：

4.本发明提出了一种考虑禁飞区约束的天基再入飞行器可达域快速计算方法，基于凸优化方法求解离 轨制动段、再入滑翔段可达域，基于椭圆近似法，提出再入滑翔段可达域快速求解方法，并基于解析预 测航程法，给出禁飞区影响下的可达域修正方法，以解决现有技术中存在的问题。
5.一种考虑禁飞区约束的天基再入飞行器可达域快速计算方法，所述考虑禁飞区约束的天基再入飞行 器可达域快速计算方法包括以下步骤：
6.s100、首先，提出基于凸优化的离轨段轨迹规划方法，将离轨制动问题转化为凸优化问题，使用二 阶锥规划方法进行求解，得到单条离轨段最优轨迹；
7.s200、其次，在单条离轨段最优轨迹的基础上，采用遍历边界的方式进行可达域计算，将可达域边 界单侧等分为n份，设计三维离轨轨迹规划及可达域计算流程，计算出可达域的最近点、最远点、左边 界和右边界；
8.s300、然后，提出基于凸优化的再入滑翔段轨迹规划方法，将再入段轨迹优化问题转化为凸优化问题，得到单条再入滑翔段最优轨迹；
9.s400、进一步地，改变性能指标，通过遍历边界得到再入段可达域；
10.s500、最后，合并离轨段与再入滑翔段可达域，结合解析预测法计算禁飞区影响的再入滑翔段可达能力，最终得到考虑禁飞区的天基再入滑翔飞行器可达域边界。
11.进一步的，在s100中，具体包括以下步骤：
12.s110、状态方程归一化：
13.定义为推力—加速度向量，其中t为推力与发动机初始重量之比，m为飞行器当前质量与初始质量之比，定义z＝lnm，将r，v分别对r0和进行归一化，并将时间t对进行归一化，得到更改后的状态方程，
[0014][0015]
其中，σ为一松弛变量，代表推力-加速度，状态量为x＝(r,v,z)
t
，控制量为u＝(τ,σ)
t
，将状态方程写为矩阵形式，
[0016][0017]
其中，v
ex
＝i
sp
g0，r＝‖r‖，发动机推力加速度上界及定义松弛关系约束分别为，
[0018]
0≤σ≤t
max
e-z
(3)
[0019]
||τ||≤σ(4)
[0020]
当使用迭代方法求解最优解时，式(4)看作线性方程；
[0021]
s120、求解初始轨迹：
[0022]
首先需要一条初始轨迹，在此基础进行优化，在此选取到达固定目标点，且指定离轨过程时间问题的燃料最优解，
[0023]
tf＝t0(5)
[0024]
将离轨过程分为n段，每一段时间间隔为求解n 1个状态变量及控制变量，飞行器初始位置、速度与终端时刻位置、速度约束可表示为：
[0025][0026][0027]
离轨过程约束为式(2)-式(7)，由于地心距r为归一化后的变量，在离轨过程中近似为1，在对结果 精确度要求不高的情况下，式(2)中的r近似为1，变为线性方程，燃料指标函数为：
[0028][0029]
由于所有约束及目标函数均为凸的，此问题为一标准的凸优化问题；
[0030]
s130、迭代求解时间最优离轨制动问题：
[0031]
引入时间变量，此时式(2)与式(7)均为非线性等式约束，同时，需要在初始解的基础上进行多次迭 代求解，
[0032]
当凸优化问题中含有非线性等式约束hi(y)＝0,i＝1,...,q时，假设当前迭代结果为y
[k]
，其中y＝(x,u,tf)，求解y
[k 1]
的过程如下：
[0033]
在y
[k]
处对函数进行泰勒展开，hi(y)＝0可以写为
[0034][0035]
其中，为二阶修正，需要对其评估，
[0036]
同时，又有展开式，
[0037][0038]
其中，为高阶修正，只需近似求解便近似为
[0039]
凸优化问题的约束为，
[0040][0041]
||y-y
[k]
||≤ρ(12)
[0042]
优化结果为y
p
，将在y
[k]
处一阶展开
[0043][0044]
等式两侧同乘(y
p-y
[k]
)
t
/2，得到，
[0045][0046]
将式(13)左侧作为对的近似，新的约束为，
[0047][0048]
||y-y
[k]
||≤ρ(15)
[0049]
将非凸约束hi(y)＝0转换为凸约束，
[0050]
时间最优的离轨制动问题中，除了考虑式(2)、式(3)、式(4)之外，其余约束为，
[0051][0052]
z(n 1)＞zf(17)
[0053]
指标函数为：
[0054]
j＝t(18)
[0055]
式(2)的离散形式以及式(16)为非线性约束，将其按前述可得到凸函数的约束条件，并附加约束：
[0056][0057]
即求得y
[k 1]
以及u
[k 1]
，多次迭代直至收敛，得到时间最优离轨制动结果。
[0058]
进一步的，在s200中，具体的，采用遍历边界的方式进行可达域计算，假设将可达域边界单侧等分为n份，设计三维离轨轨迹规划及可达域计算流程，包括以下步骤：
[0059]
s210、计算可达域最近点：以轨道面内射程最小为优化指标，求解射程最近的轨迹，记录最小射程r
min
；
[0060]
s220、计算可达域最远点：以轨道面内射程最大为优化指标，求解射程最远的轨迹，记录最大射程r
max
；
[0061]
s230、计算左边界：首先基于有限推力制导算法进行仿真，得到飞行器剩余燃料质量mr，将mr等分为n份，设燃料耗尽时飞行器质量为m0，对于第i(i＝1,2...n-1)个边界点，最近边界点编号为i＝1，近似认为剩余燃料中有i
·
mr/n的燃料用于改变飞行方向，产生横向运动，设置飞行器首先进行横向制动角分量β＝90
°
，符号恒为正的运动，当消耗完该部分燃料后，再进行纵向的离轨制动，纵向离轨制动算法为上节算法；
[0062]
s240、计算右边界：采用s230的方法，令横向制动角分量β符号恒为负，即得到右边界轨迹。
[0063]
进一步的，在s300中，具体的，
[0064]
考虑地球自转的影响，建立如下再入飞行器运动方程，
[0065][0066]
飞行器再入过程中，表示为，
[0067][0068]
其中，q和n分别为热流密度、动压及过载，kq为热流密度系数，q
max
和n
max
分别为 峰值热流密度、峰值动压及过载峰值，
[0069]
至此，得到该最优控制问题如下，
[0070][0071]
对于非线性系统，
[0072][0073]
存在状态变量和控制变量的边界约束，
[0074][0075]
边界约束条件，
[0076][0077]
其中，ψ
ol
、ψ
ou
分别为初始时刻ψ的上下界，ψ
fl
，ψ
fu
分别为终端时刻ψ的上下界，
[0078]
路径约束条件：
[0079]cl
≤c[x(τ),u(τ),τ]≤cu(26)
[0080]
其中，c
l
与cu分别为路径约束的上界与下界，
[0081]
采用基于信赖域约束的序列凸优化方法，通过求解一系列的凸优化子问题，来逐渐逼近原有问题的解，
[0082]
线性化后得到，
[0083][0084]
将状态方程在xk处一阶泰勒展开，线性化得到的动力学方程为
[0085][0086]
式(26)中的基本约束也可通过上述凸化方法处理，序列凸优化的迭代寻优过程是在给定的参考轨迹附件线性化，为保证子问题向原问题逼近，需增加信赖域约束||x-xk||≤δ，
[0087]
在式(28)中增加一个虚拟控制量，得，
[0088][0089]
在迭代过程中，在目标函数中控制虚拟控制量的范围，使其在迭代后期为小量，从而保证问题的准确性，即，
[0090]j′
＝j δ||cw||2(30)
[0091]
求解无限维最优控制问题的数值解，对其进行离散，将时间n等分，得到离散后的动力学方程为，
[0092][0093]
其中，上角标k表示第k次迭代，下角标i表示ti时刻，将式(31)重新排列组合得，
[0094][0095]
其中
[0096][0097][0098]
经上述方法将再入段轨迹优化问题转化为凸优化问题。
[0099]
进一步的，在s400中，具体包括以下步骤：
[0100]
s410、进行轴线射程最近点计算：以初始射向上纵程最近为性能指标，计算得出该弹道，该弹道记为l1；
[0101]
s420、进行轴线射程最远点计算：以射程最远为性能指标，计算得出该弹道，该弹道记为l2；
[0102]
s430、进行左右边界计算：以l1和l2的终端射程为两个端点，平均n等分，从小到大均匀取n-1个点为终端射程约束，分别以横程最大、最小为性能指标进行优化，得到左右边界各n-1条，记左侧边界弹道为l3～l
n 1
，右侧边界弹道为l
n 2
～l2×n；
[0103]
s440、轴线计算：以l1和l2的终端射程为两个端点，平均n等分，从小到大均匀取n-1个点为终端射程约束，终端横程约束为0，以飞行时间为性能指标进行优化，可以得到轴线上n-1条弹道，分别为l2×
n 1
～l3×
n-1；
[0104]
s450、近边界计算：以初始射向上左侧射程最近为性能指标，计算得到弹道l3×n，以初始射向上右侧射程最近为性能指标，计算得到弹道l3×
n 1
，分别以l3×n和l1的终端射程、l3×
n 1
和l1的终端射程为两个端点，平均m等分，从小到大均匀取m-1个点为终端射程约束，分别计算初始射向上单侧固定纵程下横程的最大最小值，得到弹道l3×
n 2
～l3×
n (m-1)
×4；
[0105]
采用椭圆近似拟合再入飞行器可达域的边界，其中椭圆的四个特征点计算过程包括：
[0106]
s451、以初始射向上纵程最小为性能指标，计算得出该弹道，该弹道记为l1，
[0107]
j＝min(sr)(35)
[0108]
s452、以初始射向上纵程最大为性能指标，计算得出该弹道，该弹道记为l2，
[0109]
j＝min(-sr)(36)
[0110]
通过最大最小射程点来计算固定射程s计算公式如下所示：
[0111]
s＝(sm (s
f-sm)/2)
·
c(37)
[0112]
式中c为经验系数，其取值范围为0≤c≤1；
[0113]
s453、固定射程s，计算射程s下以横程最大为性能指标的弹道，该弹道记为l3，
[0114]
j＝min(-cr)(38)
[0115]-0.001≤d
r-r≤0.001(39)
[0116]
s454、固定射程s，计算射程s下以横程最小为性能指标的弹道，该弹道记为l4，
[0117]
j＝min(cr)(40)
[0118]-0.001≤d
r-r≤0.001(41)按照椭圆计算公式如式(42)所示，
[0119][0120]
得到椭圆的中心e，长半轴a，短半轴b1、b2分别为
[0121][0122]
进一步的，在s500中，具体的：
[0123]
为了衡量飞行器再入滑翔段受禁飞区影响下的可达能力，基于解析法进行禁飞区影响下可达域的定量计算，
[0124]
采用dubins曲线近进行禁飞区规避轨迹的计算，假设飞行器在转弯过程为无侧滑平衡飞行，速度倾角近似为零：
[0125][0126]
升力的分量产生向心力，假设转弯半径为rz，可得：
[0127][0128]
联立式(44)和(45)，在水平面内忽略离心力，得到飞行器转弯半径的近似解析解：
[0129][0130]
则维持转弯半径rz平衡飞行的倾侧角为：
[0131][0132]
假设转弯前飞行高度h1与大气密度ρ1，转弯后的数值为h2、ρ2，由大气密度拟合公式估算，
[0133][0134]
转弯后为维持平衡飞行，大气密度需要满足，
[0135][0136]
得到高度变化量与倾侧角的关系，
[0137][0138]
得到基于能量的航程预测为，
[0139][0140]
由式(51)解析计算出由于转弯而引起的射程减小量。
[0141]
本发明的有益效果：本发明提出了一种考虑禁飞区约束的天基再入飞行器可达域快速计算方法，基 于凸优化方法求解了三维离轨制动段、再入滑翔段可达域，相比于现有的可达域计算方法，考虑了离轨 段的横向机动能力，拓宽了可达域范围。基于椭圆近似法提出了再入滑翔段可达域快速求解方法，并通 过解析法预测航程，给出了禁飞区影响下的可达域修正方法，最终实现了考虑禁飞区约束的多阶段可达 域快速计算需求。因此，本发明相比于现有基于轨迹规划方法得到了更大范围的可达域，同时相比于现 有基于轨迹优化方法求解可达域求解速度更快，具有广阔的工程应用前景。
附图说明
[0142]
图1为三维可达域示意图；
[0143]
图2为椭圆近似可达域示意图；
[0144]
图3为离轨制动段可达域计算结果，其中，图3(a)为高度变化曲线，图3(b)为速度变化曲线， 图3(c)为速度倾角变化曲线，图3(d)为速度偏角变化曲线，图3(e)为三维轨迹，图3(f)质量 变化曲线，图3(g)为纵横程变化曲线，图3(h)为离轨段可达域；
[0145]
图4为再入滑翔段可达域计算结果，其中，图4(a)为三维轨迹，图4(b)为速度变化曲线，图4 (c)为经纬度变化曲线，图4(d)为地面轨迹变化曲线，图4(e)为经纬度可达域边界，图4(f)纵 横程可达域边界；
[0146]
图5为可达域求解结果对比图；
[0147]
图6为无禁飞区可达域计算；
[0148]
图7为转弯距离与损失射程关系；
[0149]
图8为禁飞区影响下可达域边界。
具体实施方式
[0150]
下面将结合本发明实施例中的附图对本发明实施例中的技术方案进行清楚、完整
地描述，显然，所 描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技 术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0151]
参照图1-图8所示，本发明提出了一种考虑禁飞区约束的天基再入飞行器可达域快速计算方法，基 于凸优化方法求解离轨制动段、再入滑翔段可达域，基于椭圆近似法，提出再入滑翔段可达域快速求解 方法，并基于解析预测航程法，给出禁飞区影响下的可达域修正方法，以解决现有技术中存在的问题。
[0152]
一种考虑禁飞区约束的天基再入飞行器可达域快速计算方法，所述考虑禁飞区约束的天基再入飞行 器可达域快速计算方法包括以下步骤：
[0153]
s100、首先，提出基于凸优化的离轨段轨迹规划方法，将离轨制动问题转化为凸优化问题，使用二 阶锥规划方法进行求解，得到单条离轨段最优轨迹；
[0154]
s200、其次，在单条离轨段最优轨迹的基础上，采用遍历边界的方式进行可达域计算，将可达域边 界单侧等分为n份，设计三维离轨轨迹规划及可达域计算流程，计算出可达域的最近点、最远点、左边 界和右边界；
[0155]
s300、然后，提出基于凸优化的再入滑翔段轨迹规划方法，将再入段轨迹优化问题转化为凸优化问 题，得到单条再入滑翔段最优轨迹；
[0156]
s400、进一步地，改变性能指标，通过遍历边界得到再入段可达域；
[0157]
s500、最后，合并离轨段与再入滑翔段可达域，结合解析预测法计算禁飞区影响的再入滑翔段可达 能力，最终得到考虑禁飞区的天基再入滑翔飞行器可达域边界。
[0158]
进一步的，在s100中，在时间最优的离轨制动问题中，给定起始时刻的高度、速度，要求飞行器 以最短时间到达某一高度，并满足终端的速度大小及倾角要求。使用凸优化对问题进行求解的步骤如下：
[0159]
s110、状态方程归一化：
[0160]
定义为推力—加速度向量，其中t为推力与发动机初始重量之比，m为飞行器当前质量与初 始质量之比，定义z＝lnm，将r，v分别对r0和进行归一化，并将时间t对进行归一 化，得到更改后的状态方程，
[0161][0162]
其中，σ为一松弛变量，代表推力-加速度，状态量为x＝(r,v,z)
t
，控制量为u＝(τ,σ)
t
，将状 态方程写为矩阵形式，
[0163][0164]
其中，v
ex
＝i
sp
g0，r＝‖r‖，发动机推力加速度上界及定义松弛关系约束分别为，
[0165]
0≤σ≤t
max
e-z
(3)
[0166]
||τ||≤σ(4)
[0167]
当使用迭代方法求解最优解时，式(2)看作线性方程；
[0168]
s120、求解初始轨迹：
[0169]
使用迭代方法求解时间最优离轨制动问题时，首先需要一条初始轨迹，在此基础进行优化，在此选取到达固定目标点，且指定离轨过程时间问题的燃料最优解，
[0170]
tf＝t0(5)
[0171]
求解燃料最优离轨问题，需要将离轨过程分为n段，每一段时间间隔为求解n 1个状态变量及控制变量，飞行器初始位置、速度与终端时刻位置、速度约束可表示为：
[0172][0173][0174]
离轨过程约束为式(2)-式(7)，由于地心距r为归一化后的变量，在离轨过程中近似为1，在对结果精确度要求不高的情况下，式(2)中的r近似为1，变为线性方程，燃料指标函数为：
[0175][0176]
由于所有约束及目标函数均为凸的，此问题为一标准的凸优化问题；
[0177]
s130、迭代求解时间最优离轨制动问题：
[0178]
求解时间最优离轨制动问题不指定固定的再入点以及离轨时间，只要求到达固定
的再入高度时满足指定的速度大小及倾角，因此，需要引入时间变量。此时式(2)与式(7)均为非线性等式约束。同时，为保证结果的精确性，需要在初始解的基础上进行多次迭代求解。
[0179]
当凸优化问题中含有非线性等式约束hi(y)＝0,i＝1,...,q时，假设当前迭代结果为y
[k]
，其中y＝(x,u,tf)，求解y
[k 1]
的过程如下：
[0180]
在y
[k]
处对函数进行泰勒展开，hi(y)＝0可以写为
[0181][0182]
其中，为二阶修正，需要对其评估，
[0183]
同时，又有展开式，
[0184][0185]
其中，为高阶修正，只需近似求解便近似为凸优化问题的约束为，
[0186][0187]
||y-y
[k]
||≤ρ(12)
[0188]
优化结果为y
p
，将在y
[k]
处一阶展开
[0189][0190]
等式两侧同乘(y
p-y
[k]
)
t
/2，得到，
[0191][0192]
将式(13)左侧作为对的近似，新的约束为，
[0193][0194]
||y-y
[k]
||≤ρ(15)
[0195]
因此便将非凸约束hi(y)＝0转换为凸约束。
[0196]
时间最优的离轨制动问题中，除了考虑式(2)、式(3)、式(4)之外，其余约束为，
[0197][0198]
z(n 1)＞zf(17)
[0199]
指标函数为：
[0200]
j＝t(18)
[0201]
式(2)的离散形式以及式(16)为非线性约束，将其按前述可得到凸函数的约束条件，并附加约束：
[0202][0203]
即求得y
[k 1]
以及u
[k 1]
，多次迭代直至收敛，得到时间最优离轨制动结果。
[0204]
进一步的，在s200中，具体的，将剩余燃料使飞行器前期先加速，后期再减速，从而满足快速性 的需求。通过纵向平面中的运动，可以耗散飞行器剩余质量，类似地，可以在横侧向进行机动，从而拓 宽飞行器在离轨制动阶段的可达范围。三维空间中，制动发动机不再只产生纵向制动力，可通过横向的 摆角产生侧向制动力分量。
[0205]
采用遍历边界的方式进行可达域计算，假设将可达域边界单侧等分为n份，设计三维离轨轨迹规划 及可达域计算流程如下：
[0206]
s210、计算可达域最近点：以轨道面内射程最小为优化指标，求解射程最近的轨迹，记录最小射程 r
min
；
[0207]
s220、计算可达域最远点：以轨道面内射程最大为优化指标，求解射程最远的轨迹，记录最大射程 r
max
；
[0208]
s230、计算左边界：首先基于有限推力制导算法进行仿真，得到飞行器剩余燃料质量mr，将mr等 分为n份，设燃料耗尽时飞行器质量为m0，对于第i(i＝1,2...n-1)个边界点，最近边界点编号为i＝1， 近似认为剩余燃料中有i
·
mr/n的燃料用于改变飞行方向，产生横向运动，设置飞行器首先进行横向制 动角分量β＝90
°
，符号恒为正的运动，当消耗完该部分燃料后，再进行纵向的离轨制动，纵向离轨制 动算法为上节算法；
[0209]
s240、计算右边界：采用s230的方法，令横向制动角分量β符号恒为负，即得到右边界轨迹。
[0210]
进一步的，在s300中，具体的，
[0211]
考虑地球自转的影响，建立如下再入飞行器运动方程，
[0212][0213]
飞行器再入过程中，为保证飞行任务成功，需满足热流、动压、过载等路径约束，可表示为，
[0214][0215]
其中，q和n分别为热流密度、动压及过载，kq为热流密度系数，q
max
和n
max
分别为峰值热流密度、峰值动压及过载峰值，
[0216]
至此，得到该最优控制问题如下，
[0217][0218]
对于非线性系统，
[0219][0220]
存在状态变量和控制变量的边界约束，
[0221][0222]
边界约束条件，
[0223][0224]
其中，ψol、ψ
ou
分别为初始时刻ψ的上下界，ψ
fl
，ψ
fu
分别为终端时刻ψ的上下界，
[0225]
路径约束条件：
[0226]cl
≤c[x(τ),u(τ),τ]≤cu(26)
[0227]
其中，c
l
与cu分别为路径约束的上界与下界，
[0228]
针对这种复杂的多约束最优控制问题基本上均是通过线性化处理来求解。在此采用基于信赖域约束的序列凸优化方法。通过求解一系列的凸优化子问题，来逐渐逼近原有问题的解。
[0229]
线性化后得到，
[0230][0231]
将状态方程在xk处一阶泰勒展开，线性化得到的动力学方程为
[0232][0233]
式(26)中的基本约束也可通过上述凸化方法处理，序列凸优化的迭代寻优过程是
在给定的参考轨迹附件线性化，为保证子问题向原问题逼近，需增加信赖域约束||x-xk||≤δ，
[0234]
为避免迭代初期“伪不可行”的问题，降低对初值的依赖程度，同时弥补线性化带来的误差，在式(28)中增加一个虚拟控制量，可得，
[0235][0236]
在迭代过程中，在目标函数中控制虚拟控制量的范围，使其在迭代后期为小量，从而保证问题的准确性，即，
[0237]j′
＝j δ||cw||2(30)
[0238]
为求解无限维最优控制问题的数值解，必须对其进行离散。为了简单起见，状态变量像通常一样离散成等距的数值点，各个数值点出的控制量则通过b样条控制点的表达式表示，控制点个数的选取可在实际仿真中调整，优化过程中的积分项均采用梯形积分表示。仿真结果显示这种方法可以十分有效地抑制数值求解中控制量的锯齿化现象，更加有利于迭代收敛，并且减少了需要求解的变量个数。将时间n等分，可得离散后的动力学方程为，
[0239][0240]
其中，上角标k表示第k次迭代，下角标i表示ti时刻，将式(31)重新排列组合得，
[0241][0242]
其中
[0243]
[0244][0245]
将再入段轨迹优化问题转化为凸优化问题后，可离线求解单条再入滑翔轨迹。
[0246]
进一步的，在s400中，计算可达域需要得到其左边界、右边界、最远点、近边界的闭合边界，即可得到可达区域。再入段可达域求解步骤为：
[0247]
s410、进行轴线射程最近点计算：以初始射向上纵程最近为性能指标，计算得出该弹道，该弹道记为l1；
[0248]
s420、进行轴线射程最远点计算：以射程最远为性能指标，计算得出该弹道，该弹道记为l2；
[0249]
s430、进行左右边界计算：以l1和l2的终端射程为两个端点，平均n等分，从小到大均匀取n-1个点为终端射程约束，分别以横程(包含符号，左侧为正，右侧为负)最大、最小为性能指标进行优化，得到左右边界各n-1条，记左侧边界弹道为l3～l
n 1
，右侧边界弹道为l
n 2
～l2×n；
[0250]
s440、轴线计算：以l1和l2的终端射程为两个端点，平均n等分，从小到大均匀取n-1个点为终端射程约束，终端横程约束为0，以飞行时间为性能指标进行优化，可以得到轴线上n-1条弹道，分别为l2×
n 1
～l3×
n-1
；
[0251]
s450、近边界计算：以初始射向上左侧射程最近为性能指标，计算得到弹道l3×n，以初始射向上右侧射程最近为性能指标，计算得到弹道l3×
n 1
，分别以l3×n和l1的终端射程、l3×
n 1
和l1的终端射程为两个端点，平均m等分，从小到大均匀取m-1个点为终端射程约束，分别计算初始射向上单侧固定纵程下横程的最大最小值，得到弹道l3×
n 2
～l3×
n (m-1)
×4；
[0252]
为了快速求解再入段可达域，可采用椭圆近似拟合再入飞行器可达域的边界。可达域快速计算方法需要计算四个特征点，来拟合两个椭圆边界，将其近似为再入飞行器可达域边界，其中椭圆的四个特征点计算过程如下所示：
[0253]
s451、以初始射向上纵程最小为性能指标，计算得出该弹道，该弹道记为l1，
[0254]
j＝min(sr)(35)
[0255]
s452、以初始射向上纵程最大为性能指标，计算得出该弹道，该弹道记为l2(且记最小射程为sf，记最大纵程为r
max
)，
[0256]
j＝min(-sr)(36)
[0257]
通过最大最小射程点来计算固定射程s计算公式如下所示：
[0258]
s＝(sm (s
f-sm)/2)
·
c(37)
[0259]
式中c为经验系数，其取值范围为0≤c≤1；
[0260]
s453、固定射程s，计算射程s下以横程最大为性能指标的弹道，该弹道记为l3(记固定纵程为r，记最大横程为z
max
)，
[0261]
j＝min(-cr)(38)
[0262]-0.001≤d
r-r≤0.001(39)
[0263]
s454、固定射程s，计算射程s下以横程最小为性能指标的弹道，该弹道记为l4(记固定纵程为r，记最小横程为z
min
)，
[0264]
j＝min(cr)(40)
[0265]-0.001≤d
r-r≤0.001(41)
[0266]
按照椭圆计算公式如式(42)所示，
[0267][0268]
得到椭圆的中心e，长半轴a，短半轴b1、b2分别为
[0269][0270]
进一步的，在s500中，具体的：
[0271]
为了衡量飞行器再入滑翔段受禁飞区影响下的可达能力，基于解析法进行禁飞区影响下可达域的定 量计算。
[0272]
已知飞行器最大飞行能力，得到其射程-能量关系，存在禁飞区时，为了规避禁飞区，需要通过横 向机动的方式，势必增加地面轨迹的长度，当禁飞区导致的地面轨迹长度大于最大飞行能力时，必然失 败，目标点不可达；当禁飞区导致的地面轨迹长度大于最大飞行能力与横向运动导致的纵向能力损失值 差值时，同样失败，目标点不可达；当禁飞区导致的地面轨迹长度小于该差值时，任务成功，目标点可 达。因此，问题的关键是，如何定量分析横向运动导致的纵向能力损失值，通过大量仿真分析或者通过 解析计算，得到横向某个变量对能量/射程的影响，进行计算与补偿。
[0273]
采用dubins曲线近进行禁飞区规避轨迹的计算，假设飞行器在转弯过程为无侧滑平衡飞行，速度 倾角近似为零：
[0274][0275]
升力的分量产生向心力，假设转弯半径为rz，可得：
[0276][0277]
联立式(44)和(45)，在水平面内忽略离心力，得到飞行器转弯半径的近似解析解：
[0278][0279]
则维持转弯半径rz平衡飞行的倾侧角为：
[0280][0281]
假设转弯前飞行高度h1与大气密度ρ1，转弯后的数值为h2、ρ2，由大气密度拟合公式估算，
[0282][0283]
转弯后为维持平衡飞行，大气密度需要满足，
[0284][0285]
得到高度变化量与倾侧角的关系，
[0286][0287]
得到基于能量的航程预测为，
[0288][0289]
由式(51)解析计算出由于转弯而引起的射程减小量。
[0290]
以某天基再入飞行器的对地可达域计算为例，说明本发明的有效性。
[0291]
(1)离轨制动段可达域仿真
[0292]
根据离轨制动段可达域计算流程，设置离轨制动段最大仿真时长为1500s，期望再入角为θf＝-3.5
°
， 其余仿真条件不变，仿真结果如下图3所示。
[0293]
该组仿真条件下，离轨制动段横向机动能力约为400km，飞行器所携带燃料均耗尽，横向机动产生 位移的同时，能够改变飞行器原有射向，根据横向机动范围的不同，射向可改变1
°
～3
°
，从而可拓宽再 入滑翔段的横向飞行能力。通常情况下，飞行器在离轨制动阶段可达域的地面投影近似为椭圆形状，本 算例中，因为离轨制动段存在最大离轨时间约束，因此可达域远边界中并不呈现椭圆形状，若去掉或者 设置较大离轨时间限制，椭圆形状将更加明显。
[0294]
(2)再入滑翔段无禁飞区可达域仿真
[0295]
对以再入段轨迹优化问题选择不同性能指标，可以计算得到特定飞行条件下的可达域。飞行过程头 部热流密度约束为总加热量不大于5000mj/m2，动压约束末段气 动过载不大于15g。考虑控制量幅值约束为气动数据表攻角、侧滑角最大幅值(采 用内插值的方式)，即α∈[-16
°
,16
°
]，β∈[-16
°
,16
°
]，γc∈[-90
°
,90
°
]。实际系统中，控制量不 仅受到幅值约束，还要受到速率约束，控制量速率约束尚未要求，先自行设定为：先自行设定为：
[0296]
可达域为多条轨迹曲线终点空间位置的集合，连接每条曲线的终端位置便得到了飞行器在该组初始 条件下的可达域。
[0297]
为了进一步说明椭圆拟合可达域在准确性和快速性方面的优势，分别使用轨迹规划方法、凸优化方 法、椭圆近似快速计算方法进行仿真，仿真结果如图5和表1所示。
[0298][0299]
表1计算时间对比
[0300]
从中可以看出，基于椭圆近似的可达域快速计算方法具有较高的计算精度，其与基于凸优化的可达 域计算方法相比边界几乎一致，说明采用椭圆拟合精度较高，另外与基于轨迹规划的可达域计算方法相 比，边界明显较大；从图5中也可以看出，本节所提的可达域快速计算方法在计算时间上明显缩短，虽 然仍比基于轨迹规划方法计算时间长，但已具备在线计算的能力。
[0301]
(3)再入滑翔段禁飞区约束可达域仿真
[0302]
采用解析预测剩余航程进行仿真计算，假设飞行器经过最近点后，按照dubins曲线进行转弯，有/ 无禁飞区可达域计算结果如图6-图8所示。
[0303]
由图7可知，随着转弯距离的增加，射程损失增加，两者乘正相关关系，且随着转弯距离的增加， 射程损失效果明显，这是由于转弯导致飞行高度降低，阻力增大，飞行器能量损失过快引起的。当禁飞 区位于可达域上侧时，对下侧无影响，禁飞区仅对禁飞区后方区域有影响。在禁飞区后方随机设置目标 点，引导飞行器向目标点飞行，由图8可以看出，由于禁飞区的影响，可达域边界在禁飞区后“内陷”， 基于解析预测航程的方法能够快速得到禁飞区影响下的可达域边界及目标点是否可达。
[0304]
由以上的实例验证可知，本发明给出的一种考虑禁飞区约束的天基再入飞行器可达域快速计算方法 是合理且有效的。