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

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


背景技术：

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


技术实现要素：

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