CN104965418B - 一种基于弹道阻尼控制和热流解析预测的引入段制导方法 - Google Patents

Abstract

本发明一种基于弹道阻尼控制和热流解析预测的引入段制导方法，包括以下步骤：(1)引入段制导问题建模，包括动力学方程建模、过程约束和终端约束建模；(2)热流密度约束转换关系求解；通过对给定纵向升力系数下的引入段弹道进行积分，获得最大热流密度约束与最小纵向升力系数之间的显式关系；(3)定阻尼微分反馈控制方案，保证引入段的终端高度、弹道倾角的收敛；(4)升力系数和倾侧角分配，保证引入段终端攻角和倾侧角与滑翔段的期望攻角和倾侧角无缝衔接。本发明提出的引入段制导方法能够闭环解析的考虑包括热流密度过程约束和多个终端约束，具有计算简单，鲁棒性好的特点。

Description

一种基于弹道阻尼控制和热流解析预测的引入段制导方法

技术领域

本发明提供一种基于弹道阻尼控制和热流解析预测的引入段制导方法，它是涉及平稳滑翔弹道阻尼控制方法和最大热流密度解析预测方法，属于航天技术、武器技术领域。

背景技术

高超声速飞行器引入段是指从助推段结束到进入平稳滑翔状态之前的过渡弹道，它面临着复杂的环境变化和严酷的热流考验，式再入任务中的难点之一，因此有必要进行深入研究。

当前的引入段制导多采用简单的程序制导，如等攻角、等弹道倾角飞行，然后上方法并不能使得再入飞行器准确进入平稳滑翔状态。为了设计出符合引入段任务要求的弹道，间接法、遗传算法和伪谱法等优化算法被用于引入段弹道优化，但离实现在线制导还有一定距离。

弹道阻尼是最新发展的一种再入弹道控制方法，具有很强的鲁棒性。该方法最早被用于炮射火箭弹的滑翔增程，之后被引入再入滑翔弹道控制。利用弹道阻尼技术，可将引入段的制导问题简化为平稳滑翔初始状态偏差的收敛问题。而热流密度约束则可转化为高度速度走廊的下边界，从而通过预测引入段弹道的最低点来估算最大热流密度。

发明内容

本发明的目的是提供一种基于弹道阻尼控制和热流解析预测的引入段制导方法，它通过预测给定升力系数下引入段弹道的最低点，得到满足最大热流密度约束的纵向升力系数边界，并结合定阻尼微分反馈控制方法，最终获得了满足引入段制导任务需求的解析制导方法。

本发明一种基于弹道阻尼控制和热流解析预测的引入段制导方法，它包括以下几个步骤：

步骤1：引入段制导问题建模，包括动力学方程建模、过程约束和终端约束建模：

(1)动力学方程如下：

式中，h、V、γ、s分别为再入飞行器的高度、速度、弹道倾角和飞行距离； 分别为再入飞行器的高度、速度、弹道倾角和飞行距离对时间的导数；r为从地心至飞行器的径向距离，与高度的关系为h＝r-R₀，其中R₀为地球半径；g为重力加速度；L₁和D分别为升力加速度纵向分量和阻力加速度；

(2)过程约束如下：

α_min≤α≤α_max |σ|≤σ_max

式中，为热流密度；为最大热流密度；k为常数，取值为k＝5.188×10^-8；α和σ分别为攻角和倾侧角，是引入段制导问题的控制变量；α_min和α_max分别为最小攻角和最大攻角；σ_max为最大倾侧角；

(3)终点约束如下：

式中，V_f为终端速度；和分别为期望纵向升力系数和纵向升阻比；k_c为组合常数；β_r为指数大气模型常数；f_a和f_V分别为速度相关的函数，如下所示，

上式中，为对V的偏导数，f_a(V_f)和f_V(V_f)是终端速度V_f对应的f_a和f_V；

步骤2：热流密度约束转换关系求解；通过对给定纵向升力系数下的引入段弹道进行积分，获得最大热流密度约束与最小纵向升力系数之间的显式关系，包含如下内容：

(1)高度与弹道倾角关系

上式中，γ和h为当前的弹道倾角和高度，而γ_ter和h_ter为预测的引入段弹道倾角和高度；K₁和K₂均为与当前状态相关的常数，它们的表达式如下所示，

K₁＝ρ_seaSC_L1/2m K₂＝(g/V²-1/r)cosγ

上式中，ρ_sea为海平面大气密度；C_L1为纵向升力系数；

(2)纵向升力系数与最大热流密度关系

上式中，h_min为当前速度下最大热流密度对应的最小高度；C_L1min为满足过程约束的最小纵向升力系数；C_Lmin为最小攻角对应的升力系数；

步骤3：定阻尼微分反馈控制方案；引入段满足终端约束制导所需的纵向升力系数为，

上式中，C_L1g为制导所需的纵向升力系数；为滑翔段的期望纵向升力系数；σ^*为滑翔段的期望倾侧角；ζ_c为期望阻尼，通常取ζ_c＝0.7；为期望纵向升阻比；Δγ为当前弹道倾角与平稳滑翔弹道倾角的差，如下所示，

Δγ＝γ-γ_sg

上式中，γ_sg为平稳滑翔弹道倾角，如下所示，

上式中，f_a(V)和f_V(V)为当前速度V对应的f_a和f_V的值；

步骤4：升力系数和倾侧角分配策略；引入段制导所需的升力系数如下，

上式中，C_L1为实际制导的纵向升力系数，C_L1max最大纵向升力系数，由最大攻角决定；C_L1min为最小纵向升力系数；

升力系数和倾侧角分配如下，

上式中，C_Lg和σ_g分别为实际制导所需的升力系数和倾侧角；为滑翔段期望攻角对应的升力系数；

通过上述4个步骤，即通过预测给定升力系数下引入段弹道的最低点，得到满足最大热流密度约束的纵向升力系数边界，并结合定阻尼微分反馈控制方法，最终获得了满足引入段制导任务需求的解析制导方法。

本发明的优点在于：

(1)通过求解固定纵向升力系数下拉起段弹道高度的解析解，获得了该条件下的弹道最小高度，将其与最大热流密度约束对应的H-V走廊高度下边界进行对比，获得了满足最大热流密度约束的最小纵向升力系数。

(2)采用定阻尼微分反馈的方法作为引入段的控制策略，使得引入段弹道快速收敛到给定攻角和倾侧角的平稳滑翔状态。结合上面的控制边界，从而能够获得了同时满足终端高度、弹道倾角、攻角和倾侧角约束及最大热流密度约束的闭环解析制导方法。

(3)本发明简单实用，并且具有较强的抵抗拉偏能力，与伪谱法的最优弹道差别也较小。

附图说明

图1是本发明所述方法流程图。

图2是引入段解析预测弹道示意图。

图3是标称下制导方法获得的攻角曲线。

图4是标称下制导方法获得的倾侧角曲线。

图5是标称下制导方法获得的弹道曲线。

图6是标称下制导方法获得的热流密度曲线。

图7(a)是升力系数拉偏获得的攻角曲线。

图7(b)是升力系数拉偏获得的纵向弹道曲线。

图7(c)是升力系数拉偏获得的最大热流密度曲线。

图7(d)是升力系数拉偏获得的速度曲线。

图8(a)是阻力系数拉偏获得的攻角曲线。

图8(b)是阻力系数拉偏获得的纵向弹道曲线。

图8(c)是阻力系数拉偏获得的最大热流密度曲线。

图8(d)是阻力系数拉偏获得的速度曲线。

图9是与伪谱法的攻角曲线对比。

图10是与伪谱法的弹道曲线对比。

图中符号、代号说明如下：

图2中h表示高度；t表示时间；s代表时间单位秒。图3中，α代表攻角；deg代表攻角的单位度。图4中，σ代表倾侧角。图5中，s代表飞行距离，单位为千米。图6中，dQ代表热流密度，W/cm²(瓦每平方厘米)为其单位。图7(a)-(d)中，V代表速度，单位为m/s；图中图例+15％表示升力系数增大15％；-15％表示升力系数减小15％；0％表示升力系数不变。图8(a)-(d)中，图中图例+15％表示阻力系数增大15％；-15％表示阻力系数减小15％；0％表示阻力系数不变。图10中，s代表飞行距离。

具体实施方式

下面将结合附图和实施案例对本发明作进一步的详细说明。

针对引入段多终端约束和强过程约束的制导问题，本发明首先通过预测给定纵向升力系数的引入段弹道最低点，将热流密度约束转化成为了纵向升力系数下边界约束；在此基础上，采用定阻尼微分反馈控制策略使得高度和弹道倾角快速收敛到平稳滑翔状态，并设计了相应的升力系数分配策略，从而获得了满足多约束要求的引入段制导方法，具体详见图1。

本发明一种基于平稳滑翔弹道阻尼控制和最大热流密度解析预测的引入段制导方法，包括以下几个步骤：

步骤1：引入段制导问题建模

(1)动力学方程

由于再入拉前段横向弹道影响较小，并且地球自转的影响较小，因此可忽略其影响，则动力学方程简化为，

式中，h、V、γ、s分别为再入飞行器的高度、速度、弹道倾角和飞行距离； 分别为再入飞行器的高度、速度、弹道倾角和飞行距离对时间的导数；r为从地心至飞行器的径向距离，与高度的关系为h＝r-R₀，其中R₀为地球半径；g为重力加速度；L₁和D分别为升力加速度纵向分量和阻力加速度，它们的表达式为，

式中，C_L1为纵向升力系数分量，C_L1＝C_Lcosσ；C_L和C_D分别为升力系数和阻力系数；σ为倾侧角；m为飞行器质量；S为气动参考面积；ρ为大气密度，通常采用指数大气模型，如下所示，

式中，ρ_sea为海平面大气密度；β_r为指数大气模型常数，通常取β_r＝1/7200.

(2)过程约束和终端约束

在引入段，飞行器的动压和过载均较小，而最大热流密度约束则影响较大，其表达式如下，

式中，为最大热流密度；k为常数，取值为k＝5.188×10^-8；除此之外，攻角和倾侧角需要满足如下约束，

α_min≤α≤α_max|σ|≤σ_max (8)

式中，α_min和α_max分别为最小攻角和最大攻角；σ_max为最大倾侧角；

对于给定的期望攻角α^*和倾侧角σ^*，引入段终点约束如下，

式中，和分别为由α^*和σ^*求解获得的期望纵向升力系数和纵向升阻比；k_c为组合常数，k_c＝ρ_seaS/(2m)；f_a和f_V分别为速度相关的系数，表达式如下：

f_V＝2g/(gV-V³/r) (12)

式中，为对V的偏导数；f_a(V_f)和f_V(V_f)是终端速度V_f对应的f_a和f_V；

步骤2：热流密度约束转换关系求解；

(1)高度与弹道倾角关系

由式(1)和式(3)可得高度和弹道倾角关系如下，

将式(5)和式(6)带入式(13)可得，

在引入段，γ通常较小，因此可假定cosγ≈1；除此之外，在引入段的速度变化通常较小(高度下降过程中的重力加速可部分补偿阻力的减速的影响)，因此1/r-g/V²的变化较为缓慢，可认为其为常量；若进一步假定C_L1为常量，则可得，

上式中，K₁和K₂均为常数，它们的表达式如下所示，

K₁＝ρ_seaSC_L1/2m K₂＝(g/V²-1/r)cosγ

对式(15)进行积分可得，

上式中，γ和h为当前的弹道倾角和高度，而γ_ter和h_ter为预测的引入段弹道倾角和高度。

图2给出了解析预测弹道与实际弹道的对比，图中初始状态为h₀＝85km、γ₀＝0deg、V₀＝6500m/s，攻角为20deg，倾侧角为0deg。由图可以看出，式(16)预测的再入弹道最小高度略大于实际弹道的最小高度，但误差非常小。

(2)纵向升力系数与最大热流密度关系

将式(6)带入式(7)可得满足最大热流密度约束的最低高度为，

式中，h_min为满足热流密度约束的最小高度。因此，假定引入段的当前速度为V，若该段弹道的最小高度大于h_min，则符合最大热流密度约束的要求。在式(16)中，当γ_ter＝0时，则所获得的h_ter为给定C_L1下引入段弹道的最小高度，如下所示，

假定h_ter＝h_min，则可得满足最大热流密度约束的最小纵向升力系数为，

上式中，C_L1min满足最大热流密度约束的最小纵向升力系数。需要指出的是，当≥0时，弹道的高度增加同时速度减小，使得热流密度减小，因而不需考虑热流密度约束的影响，此时的最小纵向升力系数由最大倾侧角σ_max决定。综上所述，满足过程约束的最小纵向升力系数分量为，

上式中，C_Lmin为最小攻角对应的升力系数。

步骤3：定阻尼微分反馈控制方案；

定阻尼微分反馈的滑翔弹道制导策略能够适用于状态偏差很大的情况。设弹道动态系统的期望阻尼ζ_c＝0.7，则微分反馈系数为，

式中，k₂为微分反馈系数。在获得微分反馈系数后，反馈控制所需的纵向升力系数为，

式中，C_L1g为定阻尼微分反馈获得的纵向升力系数；ρ为当前的大气密度；Δγ为当前弹道倾角与平稳滑翔弹道倾角的差，

Δγ＝γ-γ_sg

式中，γ_sg为平稳滑翔弹道倾角，如下所示，

上式中，f_a(V)和f_V(V)为当前速度V对应的f_a和f_V的值。

步骤4：升力系数和倾侧角分配策略；

在获得了制导所需的纵向升力系数之后，还需考虑可用控制边界的影响。式(20)给出了C_L1的下界；而C_L1则由最大攻角决定，攻角取α_max所对应的升力系数C_Lmax即为C_L1的上界。最终可得制导所需的纵向升力系数为，

式中，C_L1为实际制导的纵向升力系数，C_Lmax最大升力系数。

在获得制导所需的纵向升力系数之后，还需对升力系数和倾侧角进行分配。若引入段不进行横向机动(σ＝0)，则C_L＝C_L1。若需要进行横向机动，则升力系数和倾侧角分配如下，

式中，C_Lg和σ_g分别为实际制导所需的升力系数和倾侧角；为滑翔段期望攻角对应的升力系数。

实施案例

为了检验解析求解算法的精度，选用CAV作为计算模型，进行数值仿真效验。仿真过程中的参数设置如表1所示，任务设置如表2所示。

表1仿真参数设置

表2期望攻角和倾侧角

参数	Case1	Case2	Case3
				期望攻角α*	15deg	15deg	15deg
期望倾侧角σ*	20deg	40deg	60deg

首先在标称情况下对引入段闭环解析制导方法进行仿真校验，得到的结果如图3至图4所示。由图3和图4可以看出，实际的攻角和倾侧角均逐渐逼近期望值，并且两者随时间变化非常光滑。图5给出了弹道曲线，可以看出引入段弹道的高度在经过一个“波谷”后，快速的收敛到平稳滑翔高度，而出现“波谷”的原因则是攻角饱和。图6则给出了热流密度曲线，可以看出，在闭环考虑热流密度的影响后，该制导方法获得的引入段弹道最大热流密度满足要求。

图7(a)-(d)和图8(a)-(d)分别给出了升力系数拉偏和阻力系数拉偏下的制导方法仿真结果，期望攻角和倾侧角与表2中的Case2相同。可以看出，本文给出的制导方法在拉偏的情况下也能够圆满完成引入段的制导任务。其中减小升力系数会增大引入段的最大热流密度(但不超过约束值)，而改变阻力系数则对引入段的最大热流密度几乎无影响。

最后，将本发明获得的引入段弹道与伪谱法优化获得的结果进行对比。其中伪谱法的目标函数如下所示，

式中，J表示目标函数；V_f为引入段的终端速度；k_α为弹道整形系数，作用是使攻角曲线光滑；t₀和t_f分别为引入段的起始时刻和终端时刻；为引入段攻角对时间的导数。仿真中的参考攻角为10deg，并且整个仿真过程倾侧角均为零，得到的结果见表3、图9和图10。

表3仿真结果对比

由表3可以看出，本文制导方法的终端速度略小于伪谱法，但飞行距离略远于伪谱法的结果。将终端速度差转化成对应的射程差，再加上两者的飞行距离差，可得本文制导方法仅比伪谱法损失20km的射程。此外，两者的最大热流密度也相差不大。图9和图10分别给出了两种方法的攻角曲线和弹道曲线，可以看出两者非常相近。

Claims (1)

1.一种基于弹道阻尼控制和热流解析预测的引入段制导方法，其特征在于：它包括以下几个步骤：

步骤1：引入段制导问题建模，包括动力学方程建模、过程约束和终端约束建模：

(1)动力学方程如下：

$\begin{array}{cc}\stackrel{\·}{h}=Vsin\mathrm{\γ}& \stackrel{\·}{s}=rVcos\mathrm{\γ}/R_0\end{array}$

$\stackrel{\·}{\mathrm{\γ}}=\frac{1}{V}\[L_1+(V^2/r-g)cos\mathrm{\γ}\]$

$\stackrel{\·}{V}=-D-gsin\mathrm{\γ}$

式中，h、V、γ、s分别为再入飞行器的高度、速度、弹道倾角和飞行距离； 分别为再入飞行器的高度、速度、弹道倾角和飞行距离对时间的导数；r为从地心至飞行器的径向距离，与高度的关系为h＝r-R₀，其中R₀为地球半径；g为重力加速度；L₁和D分别为升力加速度纵向分量和阻力加速度；

(2)过程约束如下：

$\stackrel{\·}{Q}=k\sqrt{\mathrm{\ρ}}{V}^{3.}\≤{\stackrel{\·}{Q}}_{max}$

α_min≤α≤α_max；|σ|≤σ_max

式中，为热流密度；为最大热流密度；k为常数，取值为k＝5.188×10^-8；ρ为大气密度；α和σ分别为攻角和倾侧角，是引入段制导问题的控制变量；α_min和α_max分别为最小攻角和最大攻角；σ_max为最大倾侧角；

(3)终点约束如下：

$h_f=-ln\[(g/V_f^2-1/r)/\left(k_c{C}_{L1}^{*}\right)\]/{\mathrm{\β}}_r$

${\mathrm{\γ}}_f=-\frac{\[f_a\left(V_f\right)+f_V\left(V_f\right)\](g-V^2/r)}{{\mathrm{\β}}_r{\mathrm{VK}}_1^{*}+{\mathrm{gK}}_1^{*}\[f_a\left(V_f\right)+f_V\left(V_f\right)\]}$

式中，V_f为终端速度；和分别为期望纵向升力系数和纵向升阻比；k_c为组合常数；β_r为指数大气模型常数；f_a和f_V分别为速度相关的函数，如下所示，

$f_a=(\∂C_{L1}^{*}/\∂V)/C_{L1}^{*};f_V=2g/(gV-V^3/r)$

上式中，为对V的偏导数，f_a(V_f)和f_V(V_f)是终端速度V_f对应的f_a和f_V；

步骤2：热流密度约束转换关系求解；通过对给定纵向升力系数下的引入段弹道进行积分，获得最大热流密度约束与最小纵向升力系数之间的显式关系，包含如下内容：

(1)高度与弹道倾角关系

${\mathrm{cos\γ}}_{ter}-cos\mathrm{\γ}=\frac{K_1}{{\mathrm{\β}}_r}(e^{-{\mathrm{\β}}_r{h}_{ter}}-e^{-{\mathrm{\β}}_{r}h})+K_2(h_{ter}-h)$

上式中，γ和h为当前的弹道倾角和高度，而γt_er和h_ter为预测的引入段弹道倾角和高度；K₁和K₂均为与当前状态相关的常数，它们的表达式如下所示，

K₁＝ρ_seaSC_L1/2m；K₂＝(g/V²-1/r)cosγ

上式中，ρ_sea为海平面大气密度；C_L1为纵向升力系数；S为气动参考面积；m为飞行器质量；

(2)纵向升力系数与最大热流密度关系

$C_{L1\min }=\left\{\begin{array}{cc}\frac{2{\mathrm{m\&beta;}}_r\&lsqb;(1-cos\mathrm{\&gamma;})-K_2(h_{min}-h)\&rsqb;}{{\mathrm{S\&rho;}}_{sea}(e^{-{\mathrm{\&beta;}}_r{h}_{\min }}-e^{-{\mathrm{\&beta;}}_{r}h})}& \stackrel{\&CenterDot;}{h}<0\\ C_{Lmin}{\mathrm{cos\&sigma;}}_{max}& \stackrel{\&CenterDot;}{h}\&GreaterEqual;0\end{array}\right.$

上式中，h_min为当前速度下最大热流密度对应的最小高度；C_L1min为满足过程约束的最小纵向升力系数；C_Lmin为最小攻角对应的升力系数；

步骤3：定阻尼微分反馈控制方案；引入段满足终端约束制导所需的纵向升力系数为，

$C_{L1g}=C_{L1}^{*}-\&lsqb;\frac{4{\mathrm{m\&zeta;}}_c}{\mathrm{\&rho;}VS{\mathrm{cos\&sigma;}}^{*}}\sqrt{(g-V^2/r){\mathrm{\&beta;}}_r}-\frac{2m(g-V^2/r)}{{\mathrm{\&rho;V}}^2{\mathrm{SK}}_1^{*}{\mathrm{cos\&sigma;}}^{*}}\&rsqb;C_{L1}^{*}\mathrm{\&Delta;}\mathrm{\&gamma;}$

上式中，C_L1g为制导所需的纵向升力系数；为滑翔段的期望纵向升力系数；σ^*为滑翔段的期望倾侧角；ζ_c为期望阻尼，取ζ_c＝0.7；为期望纵向升阻比；Δγ为当前弹道倾角与平稳滑翔弹道倾角的差，如下所示，

Δγ＝γ-γ_sg

上式中，γ_sg为平稳滑翔弹道倾角，如下所示，

${\mathrm{\&gamma;}}_{sg}=-\frac{\&lsqb;f_a\left(V\right)+f_V\left(V\right)\&rsqb;(g-V^2/r)}{{\mathrm{\&beta;}}_r{\mathrm{VK}}_1^{*}+{\mathrm{gK}}_1^{*}\&lsqb;f_a\left(V\right)+f_V\left(V\right)\&rsqb;}$

上式中，f_a(V)和f_V(V)为当前速度V对应的f_a和f_V的值；

步骤4：升力系数和倾侧角分配策略；引入段制导所需的升力系数如下，

$C_{L1}=\left\{\begin{array}{cc}{C}_{L1\max }& C_{L1g}>C_{L1\max }\\ C_{L1g}& C_{L1\min }\&le;C_{L1g}\&le;C_{L1\max }\\ C_{L1\min }& C_{L1g}<C_{L1\min }\end{array}\right.$

上式中，C_L1为实际制导的纵向升力系数，C_L1max最大纵向升力系数，由最大攻角决定；C_L1min为最小纵向升力系数；

升力系数和倾侧角分配如下，

$C_{Lg}=\left\{\begin{array}{cc}{C}_{L1}& C_{L1}>C_L^{*}\\ C_L^{*}& C_L^{*}{\mathrm{cos\&sigma;}}_{max}\&le;C_{L1}\&le;C_L^{*}\\ C_{L1}/{\mathrm{cos\&sigma;}}_{max}& C_{L1}<C_L^{*}{\mathrm{cos\&sigma;}}_{max}\end{array}\right.$

${\mathrm{\&sigma;}}_g=\left\{\begin{array}{cc}0& C_{L1}>C_L^{*}\\ \arccos (C_{L1}/C_L^{*})& C_L^{*}{\mathrm{cos\&sigma;}}_{max}\&le;C_{L1}\&le;C_L^{*}\\ {\mathrm{\&sigma;}}_{\max }& C_{L1}<C_L^{*}{\mathrm{cos\&sigma;}}_{max}\end{array}\right.$

上式中，C_Lg和σ_g分别为实际制导所需的升力系数和倾侧角；为滑翔段期望攻角对应的升力系数；

通过上述4个步骤，即通过预测给定升力系数下引入段弹道的最低点，得到满足最大热流密度约束的纵向升力系数边界，并结合定阻尼微分反馈控制方法，最终获得了满足引入段制导任务需求的解析制导方法。