基于弹道解析解的运载火箭助推段制导方法及装置

1.本发明涉及航天技术、武器技术领域，尤其是涉及一种基于弹道解析解 的运载火箭助推段制导方法及装置。


背景技术：

2.目前运载火箭助推段弹道通常可分为大气层内的程序飞行段和大气层 外的制导飞行段，其中程序飞行段采用依据时间或者速度设定的开环控制 指令，通常不进行反馈调整；制导飞行段则采用摄动制导或者迭代制导进行 闭环修正。当前大气层内的闭环制导方法通常采用最优控制理论进行建模， 并采用数值优化算法进行求解，如修正打靶法、配点法、最小值算法等。但 上述方法在实现过程中计算量较大、导致效率较低，且实用性较低。


技术实现要素：

3.有鉴于此，本发明的目的在于提供一种基于弹道解析解的运载火箭助 推段制导方法及装置，以提高制导效率及实用性。
4.第一方面，本发明实施例提供了一种基于弹道解析解的运载火箭助推 段制导方法，包括：建立运载火箭的助推段弹道分步求解模型；基于助推段 弹道分步求解模型，利用正则摄动理论求取运载火箭的攻角为零的情况下， 助推段的弹道解析解；基于助推段弹道分步求解模型及运载火箭的攻角为 零的情况下助推段弹道的解析解，利用正则摄动方法求取运载火箭的攻角 不为零的情况下，助推段的弹道解析解；根据弹道解析解以及预先确定的制 导策略，求取运载火箭在助推段的最优制导指令。
5.结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方 式，其中，建立运载火箭的助推段弹道分步求解模型的步骤，包括：基于预 先确定的运载火箭的助推段动力学模型，建立运载火箭的攻角为零的情况 下，助推段弹道的动力学模型；基于运载火箭的助推段动力学模型，建立运 载火箭的攻角不为零的情况下，助推段弹道的增量动力学模型；将助推段弹 道的动力学模型及增量动力学模型确定为助推段弹道分步求解模型。
6.结合第一方面的第一种可能的实施方式，本发明实施例提供了第一方 面的第二种可能的实施方式，其中，基于助推段弹道分步求解模型，利用正 则摄动理论求取运载火箭的攻角为零的情况下，助推段的弹道解析解的步 骤，包括：基于运载火箭的攻角为零的情况下，助推段弹道的动力学模型， 利用正则摄动理论建立第一正则摄动模型；基于第一正则摄动模型，生成零 攻角弹道的零阶项微分方程及零攻角弹道的一阶项微分方程；基于零攻角 弹道的零阶项微分方程，得到零攻角弹道的零阶项解析解；基于零攻角弹 道的一阶项微分方程，得到零攻角弹道的一阶项解析解；基于零攻角弹道的 零阶项解析解及零攻角弹道的一阶项解析解，生成运载火箭的攻角为零的 情况下，助推段的近似弹道解析解；将近似弹道解析解确定为运载火箭的攻 角为零的情况下，助推段的弹道解析解。
7.结合第一方面的第二种可能的实施方式，本发明实施例提供了第一方 面的第三种可能的实施方式，其中，基于助推段弹道分步求解模型及运载火 箭的攻角为零的情况下
助推段弹道的解析解，利用正则摄动方法求取运载 火箭的攻角不为零的情况下，助推段的弹道解析解的步骤，包括：基于运 载火箭的攻角不为零的情况下，助推段弹道的增量动力学模型，利用正则摄 动理论建立第二正则摄动模型；基于第二正则摄动模型，生成弹道增量的零 阶项微分方程及弹道增量的一阶项微分方程；基于弹道增量的零阶项微分 方程，得到弹道增量的零阶解析解；基于弹道增量的一阶项微分方程，得到 弹道增量的一阶解析解；基于弹道增量的零阶项解析解及弹道增量的一阶 项解析解，生成弹道增量的近似解析解；将弹道增量的近似解析解与运载火 箭的攻角为零的情况下，助推段的弹道解析解之和确定为运载火箭的攻角 不为零的情况下，助推段的弹道解析解。
8.结合第一方面，本发明实施例提供了第一方面的第四种可能的实施方 式，其中，根据弹道解析解以及预先确定的制导策略，求取运载火箭在助推 段的最优制导指令的步骤，包括：基于弹道解析解，得到终端状态与攻角曲 线之间的关系表达式；基于关系表达式，建立终端状态偏差的修正模型；通 过切比雪夫插值多项式对预先确定的运载火箭在助推段的性能指标进行近 似，得到离散的性能指标；基于终端状态偏差的修正模型以及离散的性能指 标，建立最优制导修正的解析求解模型；基于最优制导修正的解析求解模型 以及预先确定的制导策略，求解最优制导指令。
9.第二方面，本发明实施例还提供一种基于弹道解析解的助推段制导装 置，包括：模型建立模块，用于建立运载火箭的助推段弹道分步求解模型； 第一求解模块，用于基于助推段弹道分步求解模型，利用正则摄动理论求取 运载火箭的攻角为零的情况下，助推段的弹道解析解；第二求解模块，用于 基于助推段弹道分步求解模型及运载火箭的攻角为零的情况下助推段弹道 的解析解，利用正则摄动方法求取运载火箭的攻角不为零的情况下，助推段 的弹道解析解；最优制导模块，用于根据弹道解析解以及预先确定的制导策 略，求取运载火箭在助推段的最优制导指令。
10.结合第二方面，本发明实施例提供了第二方面的第一种可能的实施方 式，其中，上述模型建立模块还用于：基于预先确定的运载火箭的助推段动 力学模型，建立运载火箭的攻角为零的情况下，助推段弹道的动力学模型； 基于运载火箭的助推段动力学模型，建立运载火箭的攻角不为零的情况下， 助推段弹道的增量动力学模型；将助推段弹道的动力学模型及增量动力学 模型确定为助推段弹道分步求解模型。
11.结合第二方面，本发明实施例提供了第二方面的第二种可能的实施方 式，其中，上述第一求解模块还用于：基于运载火箭的攻角为零的情况下， 助推段弹道的动力学模型，利用正则摄动理论建立第一正则摄动模型；基于 第一正则摄动模型，生成零攻角弹道的零阶项微分方程及零攻角弹道的一 阶项微分方程；基于零攻角弹道的零阶项微分方程，得到零攻角弹道的零阶 项解析解；基于零攻角弹道的一阶项微分方程，得到零攻角弹道的一阶项解 析解；基于零攻角弹道的零阶项解析解及零攻角弹道的一阶项解析解，生成 运载火箭的攻角为零的情况下，助推段的近似弹道解析解；将近似弹道解析 解确定为运载火箭的攻角为零的情况下，助推段的弹道解析解。
12.第三方面，本发明实施例还提供一种电子设备，包括处理器和存储器， 存储器存储有能够被处理器执行的机器可执行指令，处理器执行机器可执 行指令以实现上述方法。
13.第四方面，本发明实施例还提供一种机器可读存储介质，机器可读存储 介质存储有机器可执行指令，机器可执行指令在被处理器调用和执行时，机 器可执行指令促使处理
器实现上述方法。
14.本发明实施例带来了以下有益效果：
15.本发明实施例提供了一种基于弹道解析解的运载火箭助推段制导方法 及装置，首先建立运载火箭的助推段弹道分步求解模型，然后基于助推段弹 道分步求解模型，利用正则摄动理论求取运载火箭的攻角为零的情况下，助 推段的弹道解析解，进一步基于助推段弹道分步求解模型及运载火箭的攻 角为零的情况下助推段弹道的解析解，利用正则摄动方法求取运载火箭的 攻角不为零的情况下，助推段的弹道解析解；最后根据弹道解析解以及预先 确定的制导策略，求取运载火箭在助推段的最优制导指令。该方式通过求取 运载火箭在助推段的弹道解析解，并基于此解析解建立了最优制导指令的 解析求解模型，避免了进行数值积分预测，提高了制导效率，且实用性较高。
16.本发明的其他特征和优点将在随后的说明书中阐述，并且，部分地从说 明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优 点在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。
17.为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施 例，并配合所附附图，作详细说明如下。
附图说明
18.为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下 面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍， 显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通 技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其 他的附图。
19.图1本发明实施例提供的一种基于弹道解析解的运载火箭助推段制导 方法的流程图；
20.图2为本发明实施例提供的函数ξv(v)与拟合插值多项式对比图；
21.图3为本发明实施例提供的函数f
ρ
(v)与拟合插值多项式对比图；
22.图4为本发明实施例提供的四种不同终端约束的高度vs弹道倾角曲 线图；
23.图5为本发明实施例提供的四种不同终端约束下每一个制导周期的计 算时间示意图；
24.图6为本发明实施例提供的制导周期的计算时间的示意图；
25.图7为本发明实施例提供的一种基于弹道解析解的运载火箭助推段制 导方法装置的结构示意图；
26.图8为本发明实施例提供的一种电子设备的结构示意图。
具体实施方式
27.为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发 明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述， 显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常 在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来 布置和设计。
28.因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限 制要求保
护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发 明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得 的所有其他实施例，都属于本发明保护的范围。
29.目前助推段弹道通常可分为大气层内的程序飞行段和大气层外的制导 飞行段，其中程序飞行段采用依据时间或者速度设定的开环控制指令，通常 不进行反馈调整；制导飞行段则采用摄动制导或者迭代制导进行闭环修正。 但对于助推—滑翔高超声速飞行器来说，助推段的中的高度越高则再入拉 起段的最大热流密度越大。这就要求助推段降低终端高度，从而使得助推段 弹道的大部分或者全部位于大气层内。
30.当前大气层内的闭环制导方法通常采用最优控制理论进行建模，并采 用数值优化算法进行求解，如修正打靶法、配点法、最小值算法等。但上述 方法均存在着计算量大、难于在线应用等问题。
31.基于此，本发明实施例提供的一种基于弹道解析解的助推段制导方法、 装置以及电子设备，可以于各种飞行器的助推段制导场景中。
32.本发明实施例提供了一种基于弹道解析解的运载火箭助推段制导方法， 如图1所示，该方法包括以下步骤：
33.步骤s100，建立运载火箭的助推段弹道分步求解模型。
34.通常情况下，首先基于预先确定的运载火箭的助推段动力学模型，建立 运载火箭的攻角为零的情况下，助推段弹道的动力学模型；然后基于运载火 箭的助推段动力学模型，建立运载火箭的攻角不为零的情况下，助推段弹道 的增量动力学模型；最后将助推段弹道的动力学模型及增量动力学模型确 定为助推段弹道分步求解模型。
35.步骤s102，基于助推段弹道分步求解模型，利用正则摄动理论求取运 载火箭的攻角为零的情况下，助推段的弹道解析解。
36.在具体实现过程中，可以基于运载火箭的攻角为零的情况下，助推段弹 道的动力学模型，利用正则摄动理论建立第一正则摄动模型；然后基于第一 正则摄动模型，生成零攻角弹道的零阶项微分方程及零攻角弹道的一阶项 微分方程；再基于零攻角弹道的零阶项微分方程，得到零攻角弹道的零阶项 解析解；并基于零攻角弹道的一阶项微分方程，得到零攻角弹道的一阶项解 析解；进而基于零攻角弹道的零阶项解析解及零攻角弹道的一阶项解析解， 生成运载火箭的攻角为零的情况下，助推段的近似弹道解析解；最后将近似 弹道解析解确定为运载火箭的攻角为零的情况下，助推段的弹道解析解。
37.步骤s104，基于助推段弹道分步求解模型及运载火箭的攻角为零的情 况下助推段弹道的解析解，利用正则摄动方法求取运载火箭的攻角不为零 的情况下，助推段的弹道解析解。
38.在具体实现时，首先基于运载火箭的攻角不为零的情况下，助推段弹道 的增量动力学模型，利用正则摄动理论建立第二正则摄动模型；然后基于第 二正则摄动模型，生成弹道增量的零阶项微分方程及弹道增量的一阶项微 分方程；再基于弹道增量的零阶项微分方程，得到弹道增量的零阶解析解； 并基于弹道增量的一阶项微分方程，得到弹道增量的一阶解析解；进而基于 弹道增量的零阶项解析解及弹道增量的一阶项解析解，生成弹道增量的近 似解析解；最后将弹道增量的近似解析解与运载火箭的攻角为零的情况下， 助推段的弹道解析解之和确定为运载火箭的攻角不为零的情况下，助推段 的弹道解析解。
39.步骤s106，根据弹道解析解以及预先确定的制导策略，求取运载火箭 在助推段的最优制导指令。
40.在具体实现过程中，首先基于弹道解析解，得到终端状态与攻角曲线之 间的关系表达式；然后基于关系表达式，建立终端状态偏差的修正模型；再 通过切比雪夫插值多项式对预先确定的运载火箭在助推段的性能指标进行 近似，得到离散的性能指标；进而基于终端状态偏差的修正模型以及离散的 性能指标，建立最优制导修正的解析求解模型；最后基于最优制导修正的解 析求解模型以及预先确定的制导策略，求解最优制导指令。
41.本发明实施例提供了一种基于弹道解析解的运载火箭助推段制导方法， 首先建立运载火箭的助推段弹道分步求解模型，然后基于助推段弹道分步 求解模型，利用正则摄动理论求取运载火箭的攻角为零的情况下，助推段的 弹道解析解，进一步基于助推段弹道分步求解模型及运载火箭的攻角为零 的情况下助推段弹道的解析解，利用正则摄动方法求取运载火箭的攻角不 为零的情况下，助推段的弹道解析解；最后根据弹道解析解以及预先确定的 制导策略，求取运载火箭在助推段的最优制导指令。该方式通过求取运载火 箭在助推段的弹道解析解，并基于此解析解建立了最优制导指令的解析求 解模型，避免了进行数值积分预测，提高了制导效率，且实用性较高。
42.针对高超声速飞行器助推段制导的问题，本发明提出了另一种基于弹 道解析解的助推段制导方法，该方法在图1所示的方法基础上实现。该方 法中，首先采用正则摄动方法推导了助推段弹道的解析解，此解析解由零攻 角弹道解析解及非零攻角引起的弹道增量的解析解两部分组成，其表示为 关于插值节点处攻角值的函数，可以以较高的精度预测助推段的终端速度、 弹道倾角及高度。随后，基于此解析解设计了满足强终端约束的最优制导方 法。此制导方法利用解析解将最优控制问题离散为最优化问题，将所求变量 由攻角曲线转化为插值节点处的攻角值，随后利用线性近似迭代求得优化 变量，最后利用切比雪夫插值多项式拟合出攻角曲线。
43.上述方法包括以下几个步骤：
44.步骤1：建立助推段动力学模型：
45.忽略地球自转，航迹坐标系下的助推器纵向运动方程如下述公式(1) 所示：
[0046][0047]
其中，v、γ、h和m分别为助推器的速度、弹道倾角、高度和质量； 为助推器的速度对时间的导数；为助推器的弹道倾角对时间的导数；为助推器的高度对时间的导数；为助推器的质量对时间的导数；α为助推 器的攻角；r为助推器质心与地心的距离；qm为助推器的发动机的质量秒流 量；p为助推器的推力，p与发动机的质量秒流量和当地大气压强相关；qm为助推器的质量秒流量，也就是助推器质量每秒钟减少的数量，这里的qm是正数，因此质量对时间的导数就是l和d分别为助推器的升力 和阻力，与飞行器的攻角、马赫数及动压相关；g是空中实时的重力加速度。
[0048]
其中，p通过下述公式表示：p＝i
spqmg0-pase；
[0049]
式中，i
sp
为助推器发动机比冲；se为发动机尾喷管面积；g0为地面重力 加速度；pa为当地大气压强，与飞行高度相关。
[0050]
l、d分别为助推器的升力和阻力，表达式如下：
[0051][0052]
上式中，ρ为当地大气密度，与飞行高度相关；s为助推器气动参考面 积；c
l
和cd分别为升力系数和阻力系数，两者可以分别拟合为攻角的一次 函数和二次函数，即：
[0053][0054]
其中，为升力线系数，c
d0
为零升阻力系数，为诱导阻力系数，这 些系数均可拟合为关于马赫数ma的函数。
[0055]
步骤2：助推段零攻角弹道解析求解，包括正则摄动模型的建立以及零 阶项和一阶项的解析求解：
[0056]
2.1建立正则摄动模型
[0057]
当助推段零攻角飞行时，显然有利用 了这一关系，ρ为当地大气密度，与飞行高度相关，s为助推器气 动参考面积，c
l
为升力系数。从而可得零攻角弹道的运动学方程(也称为运 动模型)为：
[0058][0059]
其中，m0和t分别为助推器的初始质量和飞行时间；下标
‘
b’表征零 攻角弹道，vb、hb、gb和ρb为零攻角弹道的速度、高度、重力加速度和大气 密度；θb为与零攻角弹道的弹道倾角γb的倾角相关变量， θb＝ln[(1 sinγb)/(1-sinγb)]；c
d0
为零升阻力系数；rb为零攻角弹道的质心与 地心的距离，γb为标控标控弹道的弹道倾角，为推力的平均值，中的 e＝2.7，是自然数对数。
[0060]
由公式(3)中可以看出vb、θb和hb之间存在着复杂的耦合关系，因此 无法直接解析求解。这里采用正则摄动方法来解决这一难题，引入一个小参 数ε并设为等于一个常数k。根据正则摄动方法，公式(3)的状态方程 需要改写为：
[0061][0062]
其中，xb＝[v
b θ
b hb]
t
，xb为零攻角弹道的状态量，为xb对时间的导 数。
[0063]
公式(3)中，零攻角弹道的速度微分方程等号右边可分为两个部分， 如下所示：
[0064][0065]
上式中，f
ave
为由平均轴向力，f
ε
为摄动轴向力，它们的表达式分别如 下：
[0066][0067][0068]
上式中，m
ave
、γ
ave
、v
ave
、ρ
ave
和g
ave
分别是平均质量、平均弹道攻角、 平均速度、平均密度以及平均重力加速度，这些值可以看作是初值和终端值 的平均值。在助推段，推力通常远远大于气动阻力和重力轴向分量的变化。 因此，摄动轴向力足够小，可以视为修正项。此外，零攻角时弹道倾角的变 化较小。因此，可以将sinγb近似为如下所示关于δθb(δθb＝θ
b-θ0)的线性 函数：
[0069]
sinγb≈sinγ0 c1δθbꢀꢀꢀ
(8)；
[0070]
式中，c1＝(sinγ
bf-sinγ0)/(θ
bf-θ0)
ꢀꢀꢀ
(9)；
[0071]
其中,γ
bf
和θ
bf
分别是γb和θb的终端值，θ0为θb的初始 值。
[0072]
那么，零攻角弹道的动力学方程改写为(即公式3可以改写为)：
[0073][0074]
根据正则摄动理论，vb、θb和hb可以表示为如下所示的关于ε的多项式：
[0075][0076]
上式中，上标(i)表征正则摄动的第i阶项。将式(11)代入式(10)并 进行泰勒级数展开，可以得到各阶动力学方程，其中零攻角弹道的零阶动力 学方程如公式(12)所示：
[0077][0078]
需要注意的是公式(12)中，为了将θb和hb的微分方程解耦，由于地心 距r的变化量比其初始值r0小得多，与地心距r相关的项在初始值处r0进行 泰勒展开。
[0079]
零攻角弹道的一阶动力学方程如公式(13)所示：
[0080]
[0081]
其中，上标
‘
(0)’和
‘
(1)’分别表征正则摄动的零阶项和一阶项；v
b(0)
、和分别为助推段标控弹道速度、弹道倾角相关变量和高度的零阶项，和分别为它们对时间的一阶导数；v
b(1)
、和分别为助推段 标控弹道速度、弹道倾角相关变量和高度的一阶项，和分别为 它们对时间的一阶导数。f
ave
为与平均推力、阻力和弹道倾角相关的常数， 即f
ave
是平均轴向力，为常数；r0、γ0和gi分别为助推段标控弹道的初始地心 距、弹道倾角和重力加速度，c1为与它们相关的常数。是弹道倾角的零阶 项，可以由根据下式计算而得：
[0082][0083]
需要注意的是，由于需要注意的是，由于的微分方程中忽略了与有关 的项。
[0084]
此外，零阶项和一阶项的初始值分别为：
[0085][0086][0087]
其中，v0、θ0和h0分别为速度、弹道倾角相关项以及高度的初始值。
[0088]
2.2零阶项解析解
[0089]
利用得到的助推段正则摄动零阶项的微分方程(12)，分别进行解析积 分，即可获得速度、弹道倾角相关变量和高度的零阶项解析解。首先，对式 (13)的第一个方程积分可得：
[0090][0091]
式中，t和ve均为助推器发动机相关的常数，表达式如下：
[0092][0093]
显然，其他变量的微分方程都与v
b(0)
有关，因此这些变量的解析解将表 示为v
b(0)
的函数。为了简化，引入一个无量纲的变量，(ve为与助推 器发动机相关的常数)，并将在之后的推导过程中将其视为独立变量。由式 (16)可得，的表达式为：
[0094][0095]
其中，是变量的初值。将式代入式可得，关于时间t的导 数可表示为：
[0096][0097]
其中，为关于时间的导数。
[0098]
零阶动力学方程的公式(12)的第二和第三个微分方程分别除以式(19) 可得：
[0099][0100][0101]
式中，k1和k2为常数。
[0102]
对式(20)进行积分可得的解析解为：
[0103][0104]
式中，函数的表达式为：
[0105][0106]
其中，n、m代表函数参数。
[0107]
是如下所示的积分项：
[0108][0109]
其中，是积分中自变量的符号，代表其n次幂。
[0110]
虽然积分项很简单，但其不存在解析解。这里可以用切比雪夫插值多项 式来近似被积函数。令：
[0111][0112]
其中，li(x)(i＝0,...,n)定义为：
[0113][0114]
式中，插值节点xi(i＝0,...,n)∈(-1,1)是切比雪夫多项式的根，其表达式 为：
[0115]
xi＝cos[(2i 1)π/(2n 2)],i＝0,...,n
ꢀꢀꢀ
(27)；
[0116]
根据高斯-切比雪夫积分公式可得的表达式为：
[0117][0118]
其中，是积分变量在插值节点xi处的值；是积分权 值系数，其具体表达式可查阅相关数学工具书。
[0119]
图2中，将对式(24)数值积分得到的结果与式(28)所示利用6阶切 比雪夫插值多项式得到的结果进行了对比，可以看出利用切比雪夫插值多 项式得到的结果具有较高的精度。将式(28)代入式(22)可得的解析 解为：
[0120][0121]
其中，是函数，其表达式已在式(23)中给出，是的 解析解。
[0122]
为了推导的解析解，引入一个如下所示的积分项：
[0123][0124]
利用分步积分公式，可得的解析解为：
[0125][0126]
(31)；其中，i
θ
是式(30)所示的积分项的符号；
[0127]
其中，为中间变量，其表达式为：
[0128][0129]
那么，对式(21)进行积分可得的解析解为：
[0130][0131]
式中，
[0132][0133]
2.3一阶项解析解
[0134]
一阶微分方程的公式(13)除以公式(19)可得，v
b(1)
、和对的导数为：
[0135][0136]
式中，和均为无量纲的变量；均为无量纲的变量；
和为常值系数。
[0137]
对式(35)的第一项进行积分，可得v
b(1)
的解析解为：
[0138][0139]
其中，是助推器零升阻力系数的平均值。
[0140][0141]
由于不是关于的线性函数，因此积分项无法 进行解析求解。不过，其被积函数只有有关，因此，我们可以利用切比雪 夫插值多项式对其近似。令：
[0142][0143]
式中，如图3所示，与 原函数(即)相比，近似多项式具有足够的精度。 那么，的解析解可表示为：
[0144][0145]
对式(35)的第二个方程进行积分，可得的解析解为：
[0146][0147]
其中，和均为关于的函数，其表达式分别为：
[0148][0149]
[0150][0151]
式中，
[0152][0153][0154][0155][0156]
式(40)中同样存在两个积分项，其具体的表达式无法直接解析求得。 为获得其表达式，我们同样利用切比雪夫插值多项式对被积函数进行近似。 令：
[0157][0158]
式中，系数的表达式为：
[0159][0160]
那么，式(47)中的积分项的解可表示为
[0161][0162][0163]
式中，函数及为如下积分项：
[0164][0165][0166]
利用分步积分公式，可得及的解析解为：
[0167]
[0168][0169]
式中，和的表达式为：
[0170][0171][0172]
对于的解析解，其被积函数中的所有变量如以 及的解析解均已得到并表示为关于的函数，因此，我们可以利用切比 雪夫插值多项式对进行近似。令：
[0173][0174]
那么，的解析解可表示为：
[0175][0176]
步骤3：非零攻角弹道增量解析求解，包括弹道增量的动力学方程的推 导、建立正则摄动模型、零阶项解析求解以及一阶项解析求解：
[0177]
3.1弹道增量的动力学方程
[0178]
定义助推器以非零攻角飞行时，速度、弹道倾角以及高度相对于零攻角 解的增量分别为δv、δγ和δh，其动力学方程可为助推器实际动力学方程和 零攻角动力学方程之差得到。为了简化推导，在求解过程中采用了如下泰勒 展开近似。
[0179][0180]
得到如公式(59)所示的弹道增量的动力学方程：
[0181][0182]
其中，cd为阻力系数。在该公式中由于r＞＞δh，上式中忽略了高度的 增量对1/r的影响。公式(59)除以公式(19)，可得δv、δγ和δh对的导 数分别为：
[0183][0184]
其中，f
v1
～f
v6
、f
γ1
～f
γ5
和f
h1
～f
h4
均为的函数，其表达式分别为：
[0185][0186][0187][0188][0189][0190][0191][0192][0193]
对式(58)进行积分便可得弹道增量的解析解，但攻角曲线一般表示为 与飞行时间t的关系，这里需要将其转换为关于的函数关系。由式(18) 可得，t和之间的关系可表示为：
[0194][0195]
那么，利用式(62)可以得到攻角关于的函数关系。
[0196]
3.2建立正则摄动模型
[0197]
由式(8)可知，δv、δγ以及δh是互相耦合的，因此无法直接求得其 解析解。这里，同样采用正则摄动方法去解决耦合的问题。根据正则摄动理 论，弹道增量的动力学方程需要写成如下形式：
[0198][0199]
式中，δx＝[δv δγ δh]
t
。通过对式(60)中三个微分方程的数量级进 行分析，可将式(60)所示的动力学方程改写为：
[0200][0201]
根据正则摄动理论，将δv、δθ和δh表示为关于参数ε的多项式：
[0202][0203]
上式中，上标(i)表征正则摄动的第i阶项。将式(65)代入式(64)并 进行泰勒级数展开，可以得到各阶动力学方程，其中零阶项微分方程为：
[0204][0205]
类似地，一阶修正项的微分方程可表示为：
[0206][0207]
零阶项和一阶修正项的初值均为零。
[0208]
3.3零阶项解析求解
[0209]
显然，式(66)所示的零阶项微分方程是不耦合的。因此，按顺序对这 些微分方程
进行积分，可以得到δγ
(0)
、δv
(0)
和δh
(0)
的解析解分别为：
[0210][0211][0212][0213]
虽然上述各式中的积分项比较复杂，但被积函数只与自变量v有关。因 此，我们同样也可以利用切比雪夫插值多项式来近似被积函数以求得积分 项的解析解。那么，δγ
(0)
、δv
(0)
和δh
(0)
解析解可表示为：
[0214][0215][0216][0217]
为了便于后续的应用，这里将上述三式表示为矩阵的形式。其中，δγ
(0)
的解析解可表示为：
[0218][0219]
式中，向量α为插值节点处助推器攻角值组成的向量。
[0220][0221]
是系数矩阵，其表达式为：
[0222][0223]
式中，向量为积分权函数组成的向量，其表达式如下：
[0224][0225]fγ1
是一个对角矩阵，其第i个对角线的值为函数f
γ1
在插值节点x
vi
处的 值,即：
[0226][0227]
将式(70)代入式(69)的第二个表达式，δv
(0)
的解析解可表示为：
[0228][0229]
式中，向量α
squ
为插值节点处助推器攻角的平方组成的向量，其表达式 如下：
[0230][0231]
和是系数矩阵，其表达式为：
[0232][0233][0234]
式中，矩阵a
cb
为积分权值矩阵，其表达式如下：
[0235][0236]fv1
和f
v3
均为对角矩阵，其第i个对角线的值分别为函数f
v1
和f
v3
在插 值节点处的值，即：
[0237][0238]
将式(70)和(75)代入式(69)的第三个表达式，δh
(0)
的解析解可表 示为：
[0239][0240]
式中，系数矩阵和的表达式分别为：
[0241][0242][0243]
式(81)和(82)中的f
h1
和f
h2
均为对角矩阵，其第i个对角线的值分 别为函数f
h1
和f
h2
在插值节点处的值，即：
[0244][0245]
3.4一阶修正项解析求解
[0246]
对于一阶项，可以通过对式(67)所示的三个微分方程进行积分得到其 解析解。与零阶项求解过程类似，这里同样利用切比雪夫插值多项式近似被 积函数。那么，δγ
(1)
的解析解可表示为：
[0247][0248]
与式(80)类似，式(84)中的第一项可表示为
[0249][0250]
式中，
[0251][0252][0253]
对于式(84)中的第二项，将式(72)和(80)代入其中，可得其表达 式为：
[0254][0255]
式中，和分别为如下所示的系数矩阵：
[0256][0257][0258]
其中，i(i)是n 1维单位矩阵的第i行。那么，结合式(85)和(87), 可得：
[0259][0260]
利用切比雪夫插值多项式逼近被积函数，可得δv
(1)
可表示为：
[0261][0262]
与式(80)类似，式(90)中的第一项可表示为：
[0263][0264]
式中，
[0265][0266][0267]
对于式(90)中的第二项，将式(72)和(80)代入其中，可得其表达 式为：
[0268][0269]
式中，
[0270][0271][0272][0273]
那么，结合式(91)和(93),可得δv
(1)
的解析解为：
[0274][0275]
与式(84)和(90)类似，δh
(1)
的解析解可表示为：
[0276][0277]
将式(95)代入式(96)的第一项，其解析表达式可写为：
[0278][0279]
式中，
[0280][0281][0282][0283][0284][0285]
与式(97)类似，式(96)的第二项可表示为：
[0286][0287]
式中，
[0288][0289][0290][0291][0292]
对于式(96)的第三项，将式(72)和(75)代入其中，其可表示为：
[0293][0294]
式中，
[0295][0296][0297]
那么，综合式(97)、(99)和(101)，δh
(1)
的解析解可表示为：
[0298][0299]
式中，
[0300][0301][0302][0303][0304]
综合零阶项和一阶项的解析解，弹道增量的近似解析解可表示为
[0305][0306]
式中，a和a
squ
为如下所示的分块对角矩阵：
[0307]
a＝diag{α,α,α},a
squ
＝diag{α
squ
,α
squ
,α
squ
}
ꢀꢀꢀ
(106)；
[0308]
系数矩阵为：
[0309][0310]
其中，每一项系数值均可由零阶项和一阶项之和确定。
[0311]
步骤4：基于解析解的助推段最优制导指令求解：
[0312]
助推段制导的目的是为后续飞行创造良好的初始条件。因此，助推段的 终端状态应满足一些约束条件。基于模型预测控制思想，以助推器当前状态 为初始状态，可以将制导问题转化为求解一系列具有终端约束的最优控制 问题。
[0313]
定义的最优控制问题的性能指标为：
[0314][0315]
其中，j为自定义的性能指标，t0和tf分别为制导段的初始时刻和终端 时刻；系数n是一个可选变量，通过改变n的取值可以对攻角曲线进行整 形。
[0316]
助推器飞行过程中需要满足公式(1)所示的动力学约束以及如公式 (109)所示的终端约束：ψ(x(tf))＝0
ꢀꢀ
(109)。
[0317]
公式(1)所示的动力学约束可以利用步骤2和步骤3推导的解析解转 化为如公式
(110)所示的代数约束：
[0318][0319]
同样，利用切比雪夫插值多项式近似被积函数，式(108)所示的积分 型性能泛函可表示为：
[0320][0321]
式中，q∈r
(n 1)
×
(n 1)
是一个正定对角矩阵，其对角线上的元素为：
[0322][0323]
因此，非线性最优控制问题被转化为非线性约束优化问题，其优化变量 是向量α，目标函数为公式(111)，约束条件为公式(109)和(110)。基于 线性逼近的思想，非线性约束优化问题可以进一步转化为一系列二次规划 问题，通过迭代求解这一系列二次规划问题，便可以得到了优化变量的最优 值。
[0324]
将公式(109)和(110)在参考弹道处进行泰勒展开并忽略高阶项，可 得如式(113)所示的以实际攻角α与参考攻角α
p
的偏差为变量的线性终端 约束函数。
[0325][0326]
式中，α
p
是参考攻角向量，其各元素为参考攻角α
p
在插值节点处x
vi
的 值，即：
[0327][0328]
δα＝α-α
p
实际攻角α与参考攻角α
p
在插值节点处的差；是终端 约束函数ψ(
·
)对终端状态量x(tf)的偏导数，其表达式与ψ(
·
)的具体表达式有 关；是终端状态量对插值节点处的攻角值的偏导数，其表达式为：
[0329][0330]
式中，式中，x
pf
是参 考攻角作用下终端状态量的值，其表达式为：
[0331][0332]
对于式(111)所示的目标函数，将δα作为优化变量，其同样可改写为 如下二次型的表达式：
[0333][0334]
利用拉格朗日乘子法可以求得此二次规划问题的解，其解满足如下线 性方程组：sz＝k
ꢀꢀꢀ
(118)；
[0335]
式中，矩阵s和k分别为线性方程组的系数矩阵，表达式如下：
[0336][0337]
k＝-[qα
p
；ψ(x
pf
)]
ꢀꢀꢀ
(119)；
[0338]
z＝[δα
t
,υ
t
]
t
为待求的变量，其中υ为拉格朗日乘子向量。求解上述线 性方程组便可得到δα，它是实际攻角α与参考攻角α
p
在插值节点处的差，则 攻角在插值节点处的值为：α＝α
p
δα
ꢀꢀꢀ
(120)；
[0339]
之后，将α
p
更新为α并不断迭代计算δα直到||ψ(x(tf))||＜δ(其中，δ是 期望的误差范围)。需要注意的是，终端状态x(tf)可以表示为插值节点上的 攻角值和一些与攻角无关的系数的函数，因此在一个制导周期的多次迭代 中，只需进行一个解析积分。
[0340]
最终，t时刻的攻角值可通过切比雪夫插值多项式近似得到，其表达式 为：
[0341][0342]
式中，是t时刻变量的值，其可以通过式(18)计算得到。
[0343]
为了校验基于解析解的助推段制导方法，采用某二级运载火箭为助推 器模型，其中第一级采用程序俯仰角制导，第二级采用基于解析解的助推段 制导方法，因此仿真时只给出了第二级的结果，第二级的参数如表1所示。 仿真时初始值设置为v0＝875m/s、γ0＝20deg和h0＝32km，终端约束如表2所 示。仿真过程中，参数n的取值为0.5。
[0344]
表1 助推器模型的参数
[0345]
[0346][0347]
表2 助推段导仿真终端约束
[0348]
终端约束case 1case 2case 3case 4弹道倾角(deg)0055高度(km)40454540
[0349]
表3给出了四种情况的终端偏差值，图4及图5给出了高度-弹道倾角 曲线以及攻角曲线。从表3和图5中可以看出本文提出的制导律可以导引 助推器到达到期望的终端弹道倾角和终端高度。图5比较了本文制导律的 攻角曲线和利用高斯伪谱法得到的开环最优控制得到的攻角曲线。可以看 出，尽管模型中存在一定的偏差，但两者几乎重合，验证了所提制导律的最 优性。
[0350]
表3 助推段制导终端偏差
[0351][0352]
图6进一步给出了每一个制导周期的计算时间。从图6可以看出，每 一个制导周期的计算时间都在0.01s内，大部分分布在0.001s-0.003s之间， 远远小于制导周期，在此，制导周期的时间长度为0.1s。因此，该方法可以 很好地应用于在线制导。
[0353]
上述基于解析解的助推段制导方法具有以下优点：
[0354]
(1)提出了一种利用插值多项式表征攻角曲线与终端状态之间关系的 方法，基于此方法可将终端状态表示为插值节点处攻角值的函数。
[0355]
(2)采用正则摄动方法解析求解了助推段纵向非零攻角弹道，获得了速 度、弹道倾角和高度的高精度解析解。
[0356]
(3)给出了一种基于弹道解析解的助推段最优制导方法，能够同时满足 助推段终端高度和弹道倾角约束，同时使得攻角曲线几乎为最优值。
[0357]
对应于上述方法实施例，本发明实施例还提供一种基于弹道解析解的 助推段制导装置，如图7所示，该装置包括：
[0358]
模型建立模块700，用于建立运载火箭的助推段弹道分步求解模型；
[0359]
第一求解模块702，用于基于助推段弹道分步求解模型，利用正则摄动 理论求取运载火箭的攻角为零的情况下，助推段的弹道解析解；
[0360]
第二求解模块704，用于基于助推段弹道分步求解模型及运载火箭的攻 角为零的情况下助推段弹道的解析解，利用正则摄动方法求取运载火箭的 攻角不为零的情况下，助推段的弹道解析解；
[0361]
最优制导模块706，用于根据弹道解析解以及预先确定的制导策略，求 取运载火箭在助推段的最优制导指令。
[0362]
进一步地，上述模型建立模块还用于：基于预先确定的运载火箭的助推 段动力学模型，建立运载火箭的攻角为零的情况下，助推段弹道的动力学模 型；基于运载火箭的助推段动力学模型，建立运载火箭的攻角不为零的情况 下，助推段弹道的增量动力学模型；将助推段弹道的动力学模型及增量动力 学模型确定为助推段弹道分步求解模型。
[0363]
本发明实施例提供的基于弹道解析解的助推段制导装置，与上述实施 例提供的基于弹道解析解的运载火箭助推段制导方法具有相同的技术特征， 所以也能解决相同的技术问题，达到相同的技术效果。
[0364]
本发明实施例还提供了一种电子设备，参见图8所示，该电子设备包 括处理器130和存储器131，该存储器131存储有能够被处理器130执行的 机器可执行指令，该处理器130执行机器可执行指令以实现上述基于弹道 解析解的运载火箭助推段制导方法。
[0365]
进一步地，图8所示的电子设备还包括总线132和通信接口133，处理 器130、通信接口133和存储器131通过总线132连接。
[0366]
其中，存储器131可能包含高速随机存取存储器(ram，random accessmemory)，也可能还包括非不稳定的存储器(non-volatile memory)，例如至 少一个磁盘存储器。通过至少一个通信接口133(可以是有线或者无线)实 现该系统网元与至少一个其他网元之间的通信连接，可以使用互联网，广域 网，本地网，城域网等。总线132可以是isa总线、pci总线或eisa总线 等。总线可以分为地址总线、数据总线、控制总线等。为便于表示，图8中 仅用一个双向箭头表示，但并不表示仅有一根总线或一种类型的总线。
[0367]
处理器130可能是一种集成电路芯片，具有信号的处理能力。在实现 过程中，上述方法的各步骤可以通过处理器130中的硬件的集成逻辑电路 或者软件形式的指令完成。上述的处理器130可以是通用处理器，包括中 央处理器(central processing unit，简称cpu)、网络处理器(network processor， 简称np)等；还可以是数字信号处理器(digital signal processing，简称dsp)、 专用集成电路(application specific integrated circuit，简称asic)、现成可编 程门阵列(field-programmable gate array，简称fpga)或者其他可编程逻辑 器件、分立门或者晶体管逻辑器件、分立硬件组件。可以实现或者执行本发 明实施例中的公开的各方法、步骤及逻辑框图。通用处理器可以是微处理器 或者该处理器也可以是任何常规的处理器等。结合本发明实施例所公开的 方法的步骤可以直接体现为硬件译码处理器执行完成，或者用译码处理器 中的硬件及软件模块组合执行完成。软件模块可以位于随机存储器，闪存、 只读存储器，可编程只读存储器或者电可擦写可编程存储器、寄存器等本领 域成熟的存储介质中。该存储介质位于存储器131，处理器130读取存储器 131中的信息，结合其硬件完成前述实施例的方法的步骤。
[0368]
本发明实施例还提供了一种机器可读存储介质，该机器可读存储介质 存储有机器可执行指令，该机器可执行指令在被处理器调用和执行时，该机 器可执行指令促使处理器实现上述基于弹道解析解的运载火箭助推段制导 方法，具体实现可参见方法实施例，在此不再赘述。
[0369]
本发明实施例所提供的基于弹道解析解的运载火箭助推段制导方法、 装置和电子设备的计算机程序产品，包括存储了程序代码的计算机可读存 储介质，程序代码包括的指令可用于执行前面方法实施例中的方法，具体实 现可参见方法实施例，在此不再赘述。
[0370]
功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时， 可以存
储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术 方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以 以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包 括若干指令用以使得一台计算机设备(可以是个人计算机，网关电子设备， 或者网络设备等)执行本发明各个实施例方法的全部或部分步骤。而前述的 存储介质包括：u盘、移动硬盘、只读存储器(rom，read-only memory)、 随机存取存储器(ram，random access memory)、磁碟或者光盘等各种可 以存储程序代码的介质。
[0371]
最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对 其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通 技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修 改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换， 并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。