一种实时解析构造的空天飞行器再入制导方法与流程

本发明涉及一种实时解析构造的空天飞行器再入制导方法，属于制导控制技术领域，特别适用于空天飞行器、可重复使用运载器、高超声速飞行器的再入段的制导。

背景技术：

空天飞行器(aerospacevehicle,asv)是既能在航空领域又能在航天领域工作的新型飞行器，它结合了航空技术与航天技术。空天飞行器能像普通飞机一样水平起飞，以高超声速在大气层内飞行，并可直接加速进入地球轨道，成为航天飞行器，返回大气层后，像飞机一样在机场着陆。

再入段是空天飞行器返回着陆过程的一个重要阶段，是空天飞行器安全返回与顺利着陆的重要保障。再入段位于离轨段之后，末端能量管理段之前，其主要目的是逐步耗散飞行器的动能和势能量，控制飞行器以合适的高度、速度和航向到达末端能量管理窗口。典型的初期段飞行起始于速度约7511m/s、海拔高度约120000m，结束于末端能量管理窗口，终端标称速度约765m/s、终端标称高度约26700m。

空天飞行器再入段飞行高度和速度变化范围大，飞行状态变化剧烈。同时，再入动压、过载和热流密度等限制严重，决定了飞行走廊窄，还存高空大气环境和气动力系数的不确定性大的问题，过程约束条件强以及参数不确定性大的特点，制导控制难度大。

再入制导方法通常分为标准轨迹制导和预测校正制导两大类。标准轨迹制导方法通过设计导引律以实时跟踪事先确定的标准再入轨迹，对gnc控制器的性能要求比较低，但是较大的初始轨迹偏差以及再入过程环境扰动等因素的影响较大，可能会导致制导性能严重退化。预测校正制导方法则是通过预测再入飞行器的落点位置偏差来对制导指令进行实时修正，能够显著降低初始散布误差和飞行过程中各种偏差对制导性能的影响，提高制导精度。因此，预测校正制导成为空天飞行器再入制导方法的发展趋势。然而，目前大量的预测校正制导方法均建立在数值积分弹道方程的基础上，该类方法的计算量大，给工程实施造成了巨大的难度。

在工程实施中，航天飞机采用基于阻力加速度剖面的标准弹道再入制导方法计算量小，对于初始误差和过程干扰较小的再入任务可以实现较好的制导效果，并进一步应用在了hermes飞行器、x-33飞行器。航天飞机采用的基于阻力加速度剖面的标准弹道再入制导方法难以摆脱鲁棒性和适应性差的固有缺陷。为满足新一代空天飞行器对自主性、精确性、安全性和可靠性的苛刻要求，世界各国开展了大量的新型升力式再入制导方法的研究和验证，比较典型的为改进的加速度再入制导算法，但该方法选用三段线性曲线组成的阻力加速度剖面限制了飞行器的再入飞行能力，同时由于三段曲线不连续，不利于飞行器跟踪参考阻力剖面，此外，用到的搜索算法比较消耗机时。

技术实现要素：

本发明的技术解决问题是：克服现有技术的不足之处，提供一种实时解析构造的空天飞行器再入制导方法，该方法能够获得高的制导精度，提高了空天飞行器的纵横航程在线调整能力，同时计算简单，工程实现容易。

本发明的技术方案是：

一种实时解析构造的空天飞行器再入制导方法，包括步骤如下：

1)判断飞行器弹道倾角是否为零，若飞行器弹道倾角不为零，则进入步骤2)，若飞行器弹道倾角为零，则进入步骤3)；

2)将保持飞行器弹道倾角为零作为制导信息，将所述制导信息传递给飞行器控制系统，直至飞行器弹道倾角零后，进入步骤3)；

3)确定飞行器能量-阻力加速度走廊上界和走廊下界，根据所述飞行器能量-阻力加速度走廊上界和走廊下界，生成飞行器能量阻力加速度特征剖面；

4)修正飞行器能量阻力加速度特征剖面,获得飞行器能量阻力加速度特征剖面的修正量；

5)根据飞行器能量阻力加速度特征剖面的修正量，解析构造再入制导模型；

6)根据所述解析构造的再入制导模型，确定飞行器倾侧角指令；

7)将所述确定的倾侧角指令作为制导信息发送给飞行器控制系统，完成飞行器再入制导工作。

所述步骤3)确定飞行器能量-阻力加速度走廊上界根据下列等式确定，具体为：

所述步骤3)飞行器能量-阻力加速度走廊下界根据下列等式确定，具体为：

其中，qmax为最大动压约束，nmax为最大过载约束，为最大热流密度约束，dqmax为动压约束下的飞行器能量-阻力加速度走廊上界，dnmax为过载约束下的飞行器能量-阻力加速度走廊上界，为热流密度约束约束下的飞行器能量-阻力加速度走廊上界，deqmin为平衡滑翔约束下的飞行器能量-阻力加速度走廊下界，cl为气动升力系数，kq为与飞行器外形相关的常值参数，r为飞行器地心距，sref为飞行器参考面积，m为飞行器质量，cd为气动阻力系数，kq为与飞行器外形相关的常值参数，v为飞行器速度，g为引力加速度。

所述步骤3)生成的飞行器能量阻力加速度特征剖面，具体为：

根据所述飞行器能量-阻力加速度走廊上界和走廊下界，以飞行器能量为横坐标，以飞行器阻力加速度为纵坐标，生成飞行器能量阻力加速度特征剖面；

a)当e≥ec时，所述飞行器能量阻力加速度特征剖面包括n段二次曲线，n为正整数；其中，第i段二次曲线对应的飞行器能量阻力加速度特征剖面，具体为：

di0＝ci1e²+ci2e+ci3，eif<e≤eib；

其中，i＝1,2,3,...,n，eif、eib分别为第i条阻力加速度二次曲线的终端能量和起始能量，ci1、ci2和ci3为第i条阻力加速度二次曲线的系数，ec为二次曲线和一次曲线分界点能量，e为飞行器能量；

b)当e<ec时，所述飞行器能量阻力加速度特征剖面为一次曲线，具体为；

其中，dc为ec点对应的飞行器阻力加速度，df为ef点对应的飞行器阻力加速度，ef为再入终端的标称能量，vf为再入终端标称速度，hf为再入终端标称高度，gf为再入终端标称引力加速度，ρf为再入终端标称大气密度。

所述步骤4)修正飞行器能量阻力加速度特征剖面的方法，具体为：

41)根据所述飞行器能量阻力加速度特征剖面，确定飞行器能量阻力加速度特征剖面对应的飞行航程sd0；

42)确定飞行器初始待飞航程st0；

43)修正所述步骤3)生成的飞行器能量阻力加速度特征剖面，使得所述步骤41)确定的sd0与所述步骤42)确定的st0一致，获得修正后的所述飞行器能量阻力加速度特征剖面。

所述步骤41)飞行器能量阻力加速度特征剖面对应的飞行航程sd0，具体为:

所述步骤42)飞行器初始待飞航程st0，具体为:

st0＝r0arccos(sinφ0sinφf+cosφ0cosφfcos(λf-λ0))，

其中，φf为再入终端标称纬度，λf为再入终端标称经度，φ0为弹道倾角变化为零时飞行器所处的纬度，λ0为弹道倾角变化为零时飞行器所处的经度，r0为弹道倾角变化为零时飞行器的地心距。

所述步骤5)解析构造再入制导模型，具体为:

sdt＝sd-st，

st＝rarccos(sinφsinφf+cosφcosφfcos(λf-λ))，

其中，θ为飞行器弹道倾角，u＝cosσ，d0为当前能量对应的初始阻力加速度剖面下的阻力加速度，d为飞行器阻力加速度，δd为飞行器能量阻力加速度特征剖面的修正量，l为飞行器升力加速度，hs为大气密度特征参数，dn＝dl0+δd，r为当前飞行器地心距；φ为当前飞行器纬度；λ为当前飞行器经度；

en确定方法如下

a)当e＞ec时：

en＝ec，

b)当ef<en≤ec时：

en＝en，

所述步骤6)，确定飞行器倾侧角指令具体为：

σ＝σabssign(σ)，

σabs＝arccos(u)，

其中，a0＞0,a1＞0,a2＞0,ε0＞0，其中g^-1为g的逆矩阵，sign()为符号函数，σi-1为上一制导周期的倾侧角指令，δψd为侧向方位边界，δψ为侧向方位

本发明与现有技术相比的优点在于：

1)本发明采用阻力加速度特征剖面自主生成，以及基于此的阻力加速度剖面实时调整的航程误差调整方式，使得阻力加速度剖面曲线连续且易于跟踪，起到了航程偏差在线实时削除的作用，提高了制导精度；

2)本发明采用基于实时构造的再入制导方式，起到了提高了鲁棒性和适应性，避免了传统的标准弹道再入制导方法的固有缺陷，飞行器的纵横航程在线调整能力强和制导精度高；

3)本发明采用了解析计算实时制导参数的方式，减少了飞控机数据运算能力的要求，避免了数值预测校正计算量大，使得该方法易于工程实施。

附图说明

图1为方法流程图；

图2为本发明实施例的飞行器能量-阻力加速度走廊示意图；

图3为本发明实施例的侧向方位边界示意图；

图4为本发明实施例的速度曲线；

图5为本发明实施例的全程高度曲线；

图6为本发明实施例的全程地轨迹曲线；

图7为本发明实施例的全程弹道倾角曲线；

图8为本发明实施例的全程的倾侧角指令；

图9为本发明实施例的全程的阻力加速度曲线。

具体实施方式

从再入制导开始直到弹道倾角大小减小到零，为飞行器航迹拉起段；从飞行器航迹拉起段结束到飞行器能量减少至再入终端能量，为飞行器滑翔飞行段。

如图1所示，一种实时解析构造的空天飞行器再入制导方法，包括步骤如下：

1)判断飞行器弹道倾角是否为零，若飞行器弹道倾角不为零，则飞行器的飞行阶段为航迹拉起段，进入步骤2)，若飞行器弹道倾角为零，则飞行器的飞行阶段为滑翔飞行段，进入步骤3)；

2)将保持飞行器弹道倾角为零作为制导信息，将所述制导信息传递给飞行器控制系统，直至飞行器弹道倾角零后，进入步骤3)；

3)确定飞行器能量-阻力加速度走廊上界和走廊下界，根据所述飞行器能量-阻力加速度走廊上界和走廊下界，生成飞行器能量阻力加速度特征剖面；

4)修正飞行器能量阻力加速度特征剖面,获得飞行器能量阻力加速度特征剖面的修正量；

5)根据飞行器能量阻力加速度特征剖面的修正量，即根据获得的修正后当前飞行器能量阻力加速度特征剖面对航程偏差的影响，解析构造再入制导模型；

6)根据所述解析构造的再入制导模型，确定飞行器倾侧角指令；

7)将所述确定的倾侧角指令作为制导信息发送给飞行器控制系统，完成飞行器再入制导工作。

步骤3)确定飞行器能量-阻力加速度走廊上界根据下列等式确定，具体为：

所述步骤3)飞行器能量-阻力加速度走廊下界根据下列等式确定，具体为：

其中，qmax为最大动压约束，nmax为最大过载约束，为最大热流密度约束，dqmax为动压约束下的飞行器能量-阻力加速度走廊上界，dnmax为过载约束下的飞行器能量-阻力加速度走廊上界，为热流密度约束约束下的飞行器能量-阻力加速度走廊上界，deqmin为平衡滑翔约束下的飞行器能量-阻力加速度走廊下界，cl为气动升力系数，kq为与飞行器外形相关的常值参数，r为飞行器地心距，sref为飞行器参考面积，m为飞行器质量，cd为气动阻力系数，kq为与飞行器外形相关的常值参数，v为飞行器速度，g为引力加速度。

飞行器能量阻力加速度特征剖面是以飞行器能量为横坐标，飞行器阻力加速度为纵坐标的一条连续光滑曲线，具体由多个光滑连续的二次曲线和一次曲线组成，光滑曲线位于步骤3确定飞行器能量-阻力加速度走廊上界和走廊下界范围内部。其中，v为飞行器速度，g为引力加速度，h为飞行器高度，大气密度ρ0为海平面大气密度，hs为大气特征常数，sref为飞行器参考面积，m为飞行器质量，cd为气动阻力系数。

步骤3)生成的飞行器能量阻力加速度特征剖面，具体为：根据所述飞行器能量-阻力加速度走廊上界和走廊下界，以飞行器能量为横坐标，以飞行器阻力加速度为纵坐标，生成飞行器能量阻力加速度特征剖面；

a)当e≥ec时，所述飞行器能量阻力加速度特征剖面包括n段二次曲线，n为正整数；其中，第i段二次曲线对应的飞行器能量阻力加速度特征剖面，具体为：

di0＝ci1e²+ci2e+ci3，eif<e≤eib；

其中，i＝1,2,3,...,n，eif、eib分别为第i条阻力加速度二次曲线的终端能量和起始能量，ci1、ci2和ci3为第i条阻力加速度二次曲线的系数，ec为二次曲线和一次曲线分界点能量，e为飞行器能量；ci1、ci2和ci3由eif对应的阻力加速度dif、eib对应的阻力加速度dib以及eim(eif<eim≤eib)对应的阻力加速度dim，这三个特征点确定。利用该阻力加速度曲线，可计算出对应的飞行航程为：

b)当e<ec时，所述飞行器能量阻力加速度特征剖面为一次曲线，具体为；

其中，dc为ec点对应的飞行器阻力加速度，df为ef点对应的飞行器阻力加速度，ef为再入终端的标称能量，vf为再入终端标称速度，hf为再入终端标称高度，gf为再入终端标称引力加速度，ρf为再入终端标称大气密度。利用该阻力加速度曲线，可计算出对应的飞行航程为：

步骤4)修正飞行器能量阻力加速度特征剖面的方法，具体为：

41)根据所述飞行器能量阻力加速度特征剖面，确定飞行器能量阻力加速度特征剖面对应的飞行航程sd0；

42)确定飞行器初始待飞航程st0；

43)修正所述步骤3)生成的飞行器能量阻力加速度特征剖面，使得所述步骤41)确定的sd0与所述步骤42)确定的st0一致，获得修正后的所述飞行器能量阻力加速度特征剖面。

步骤41)飞行器能量阻力加速度特征剖面对应的飞行航程sd0，具体为:

步骤42)飞行器初始待飞航程st0，具体为:

st0＝r0arccos(sinφ0sinφf+cosφ0cosφfcos(λf-λ0))，

其中，φf为再入终端标称纬度，λf为再入终端标称经度，φ0为弹道倾角变化为零时飞行器所处的纬度，λ0为弹道倾角变化为零时飞行器所处的经度，r0为弹道倾角变化为零时飞行器的地心距。

步骤43)修正所述步骤3)生成的飞行器能量阻力加速度特征剖面，使得所述步骤41)确定的sd0与所述步骤42)确定的st0一致，获得修正后的所述飞行器能量阻力加速度特征剖面，具体为：

431)通过调整各条二次系数中eif对应的阻力加速度dif、eib对应的阻力加速度dib以及eim(eif<eim≤eib)对应的阻力加速度dim，最终使得生成的飞行器能量阻力加速度特征剖面对应的飞行航程sd0与待飞航程st0一致，即sd0＝st0。同时，应保证选取的特征点拟合出的阻力加速度曲线之间具有连续性；

432)得到航程匹配对应的eif对应的阻力加速度dif、eib对应的阻力加速度dib以及eim(eif<eim≤eib)对应的阻力加速度dim，根据上述系数确定,ci1、ci2和ci3为第i条阻力加速度二次曲线的系数。得到特征剖面的最终结果。

步骤5)根据飞行器能量阻力加速度特征剖面的修正量，解析构造再入制导模型，具体为：

51)当前飞行器能量阻力加速度特征剖面和再实时调整dddc后的飞行器能量阻力加速度特征剖面公式为：

511)当e≥ec时，当前阻力加速度剖面中第i条阻力加速度曲线为：

di＝di0+δd,eif<e≤eib；

其中，δd为飞行器能量阻力加速度特征剖面的修正量，当前飞行器能量阻力加速度特征剖面再实时调整dddc后的第i条阻力加速度剖面为：

dic＝di0+δd+dddc，eif<e≤eib；

512)当e<ec时，当前阻力加速度剖面为：

其中，ec为预先设定的二次曲线和一次曲线分界点能量，en为当前能量，当en≤ef时，dl＝df；当ef<en≤ec时，en＝en；en＞ec时，en＝ec。δd为en处飞行器能量阻力加速度特征剖面的修正量。在en处再实时调整dddc，并保持终端能量点处阻力加速度不变，实时调整后的一次阻力加速度曲线为：

513)当前阻力加速度剖面实时修正dddc对航程偏差的影响为：

其中，sdt＝sd-st，sd为当前阻力加速度对应的飞行航程，st为当前待飞航程，为导航系统提供的阻力加速度变化率。

st＝rarccos(sinφsinφf+cosφcosφfcos(λf-λ))，

其中，r为当前飞行器地心距；φ为当前飞行器纬度；λ为当前飞行器经度，当飞行器当前能量en≤eif时，ein＝eif；当eif<en≤eib时，ein＝en；en＞eib时，ein＝eib。

52)根据实时调整阻力加速度剖面对飞行器航程偏差的影响，解析构造再入纵向制导模型，具体方法为：

sdt＝sd-st，

st＝rarccos(sinφsinφf+cosφcosφfcos(λf-λ))，

其中，θ为飞行器弹道倾角，u＝cosσ，d0为当前能量对应的初始阻力加速度剖面下的阻力加速度，d为飞行器阻力加速度，δd为飞行器能量阻力加速度特征剖面的修正量，l为飞行器升力加速度，hs为大气密度特征参数，dn＝dl0+δd，r为当前飞行器地心距；φ为当前飞行器纬度；λ为当前飞行器经度；

en确定方法如下：

a)当e＞ec时：

en＝ec，

b)当ef<en≤ec时：

en＝en，

步骤6)，确定飞行器倾侧角指令具体为：

σ＝σabssign(σ)，

σabs＝arccos(u)，

其中，a0＞0,a1＞0,a2＞0,ε0＞0，其中g^-1为g的逆矩阵，sign()为符号函数，σi-1为上一制导周期的倾侧角指令，δψd为侧向方位边界，δψ为侧向方位。

实施例

初始条件为：初始高度120000m，初始速度7511.4m/s，初始弹道倾角-1.1°，初始航迹方位角46.3°，初始经度15.589°，初始纬度8.370°。

终端条件为：终端高度26740m，终端经度85.668075°，终端纬度42.194275°，终端地速765.8m/s，终端弹道倾角-8.5°，终端航迹方位角0°。

过程约束：最大动压约束12000pa，最大过载约束2.5g，驻点热流密度约束800kw/m²。图5为本实施例全程高度曲线，图6为本实施例全程地轨迹曲线，图7为本实施例全程弹道倾角曲线，图8为本实施例的全程的航迹偏角曲线，图9为本实施例的全程的阻力加速度曲线。

本发明的具体步骤如下：

步骤1)飞行器能量-阻力加速度走廊的确定，具体为：

其中，qmax＝12000pa，nmax＝2.5g，飞行器能量-阻力加速度走廊具体如图2所示。

步骤3)生成飞行器能量阻力加速度特征剖面简称阻力加速度剖面，具体方法为：

31)当e≥ec(ec＝4.4×10⁶)时，为二次曲线段，由多条二次曲线拼接而成，具体为:

利用该阻力加速度曲线，可计算出对应的飞行航程为：

32)当e<ec时，为一次曲线段，由一条一次曲线组成，具体为：

其中，dc＝12.45，df＝7.3081，ef＝5.5306×10⁵。利用该阻力加速度曲线，可计算出对应的飞行航程为：

sf0＝3.6710×10⁵，

步骤41)利用如下公式计算该阻力加速度曲线下的飞行航程：

sd0＝3.47539×10⁶，

步骤42)利用如下公式计算初始待飞航程st0：

st0＝r0arccos(sinφ0sinφf+cosφ0cosφfcos(λf-λ0))，

其中：φf＝42.194275°，λf＝85.668075°，φ0＝24.858637°，λ0＝37.355623°，r0＝6.457264956×10⁶。

步骤43)通过调整各条二次系数中eif对应的阻力加速度dif、eib对应的阻力加速度dib以及eim(eif<eim≤eib)对应的阻力加速度dim，最终使得生成的飞行器能量阻力加速度特征剖面对应的飞行航程sd0与待飞航程st0一致，即sd0＝st0。同时，应保证选取的特征点拟合出的阻力加速度曲线之间具有连续性。

最终选取的再入阻力加速度特征点及飞行器能量阻力加速度特征剖面具体如图2所示。

所述步骤4)、步骤5)根据当前飞行器能量阻力加速度特征剖面实时修正对航程偏差的影响，实时调整飞行器能量阻力加速度特征剖面，解析构造再入纵向制导模型，具体方法为：

51)根据当前飞行器能量阻力加速度特征剖面对应的飞行航程与当前待飞航程的偏差的情况，实时调整飞行器能量阻力加速度特征剖面，并解算出对飞行器航程偏差的影响。当前飞行器能量阻力加速度特征剖面和再实时调整dddc后的飞行器能量阻力加速度特征剖面公式为：

511)当e≥ec时，当前阻力加速度剖面中第i条阻力加速度曲线为：

di＝di0+δd，eif<e≤eib；

其中，δd为当前阻力加速度剖面相对于阻力加速度特征剖面的已修正量。当前飞行器能量阻力加速度特征剖面再实时调整dddc后的第i条阻力加速度剖面为：

dic＝di0+δd+dddc，eif<e≤eib；

512)当e<ec时，当前阻力加速度剖面为：

其中，当en≤ef(en为当前能量)时，dl＝df；当ef<en≤ec时，en＝en；en＞ec时，en＝ec。δd为en处当前阻力加速度相对于阻力加速度特征剖面的已修正量。在en处再实时调整dddc，并保持终端能量点处阻力加速度不变，实时调整后的一次阻力加速度曲线为：

当前阻力加速度剖面实时修正dddc对航程偏差的影响为：

其中，sdt＝sd-st，sd为当前阻力加速度对应的飞行航程，st为当前待飞航程。

st＝rarccos(sinφsinφf+cosφcosφfcos(λf-λ))，

其中，r为当前飞行器地心距，φ为当前飞行器纬度，λ为当前飞行器经度，sin为能量大于等于ec部分中的第i条阻力加速度剖面对应的飞行航程，具体为：

其中，当飞行器当前能量en≤eif时，ein＝eif；当eif<en≤eib时，ein＝en；en＞eib时，ein＝eib。

sfn为能量小于ec部分对应的一次曲线阻力加速度剖面对应的飞行航程，具体为：

根据实时调整阻力加速度剖面对飞行器航程偏差的影响，解析构造再入纵向制导模型，具体方法为：

根据实时调整阻力加速度剖面对飞行器航程偏差的影响，得出解析构造的再入制导模型为：

其中，g＝μdg0，θ为飞行器弹道倾角，u＝cosσ，d0为当前能量对应的初始阻力加速度剖面下的阻力加速度，其中，dn＝dl0+δd，飞行器速度v曲线如图4所示。

而的解析求解公式分别为：

a)能量大于等于ec部分中的第i条阻力加速度剖面曲线：

b)能量小于ec部分：

步骤6)根据解析构造的再入纵向制导模型，确定飞行器纵向制导律，并解算出倾侧角指令大小，具体方法为：

61)利用解析构造的再入制导模型，得到的再入纵向制导律为：

其中，a0＝0.0001，a1＝0.1，a2＝0.01，ε0＝1。

62)根据u，计算倾侧角指令大小计算公式为：

σabs＝arccos(u)，

所述步骤6)使用侧向制导律，根据侧向方位误差，实时确定倾侧角指令的正负值，具体方法为：

其中，sign()为符号函数，σi-1为上一制导周期的倾侧角指令，δψd为侧向方位边界具体如图3所示，δψ为侧向方位。

所述步骤6)最后根据倾侧角指令大小和倾侧角指令的正负，得出倾侧角指令结果如图8所示，具体方法为：

σ＝σabssign(σ)。

本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。