一种高超声速飞行器平稳滑翔再入弹道设计方法

摘要

一种高超声速飞行器平稳滑翔再入弹道设计方法，它包括以下步骤：1、动力学模型建立；2、平稳滑翔弹道提出；3、平稳滑翔弹道解析解求解；4、平稳滑翔弹道动态特性分析；5、平稳滑翔反馈控制设计。本发明提出了满足纵向加速度变化率最小的平稳滑翔弹道概念，建立了平稳滑翔高度动态微分方程，并给出了较为简洁的弹道振荡自然频率和阻尼表达式，证明了纯微分反馈是抑制再入滑翔弹道振荡的最佳方案，将定阻尼微分反馈用于引入段弹道的生成，简化了再入滑翔弹道的纵平面制导。

说明书

技术领域

本发明涉及一种高超声速飞行器平稳滑翔再入弹道设计方法，属于航空航天及武器技术 领域。

背景技术

高超声速飞行器是指飞行马赫数大于或者等于5的飞行器，它具有飞行速度快、突防能 力强、毁伤威力大、全球到达等明显优势，已经成为当今世界各国武器研制的热点和焦点。

平衡滑翔指的是纵向受力近似平衡的一种飞行状态。它是再入的一种重要飞行模式，具 有高度变化平缓、热流密度和动压峰值小、攻角和倾侧角曲线光滑等优点，被广泛的应用于 再入制导。但是由于再入弹道对控制变量高度敏感，在弹道优化时往往难以获得平衡滑翔弹 道；同时再入过程中的各种扰动也容易使得弹道偏离平衡滑翔状态，出现弹道振荡。因此， 为了使高超声速飞行器在再入过程中保持平衡滑翔，不仅需要获得平衡滑翔参考弹道，还需 在飞行过程中抑制扰动带来的弹道振荡。

发明内容

针对上述高超声速飞行器平衡滑翔再入问题，提出了满足纵向加速度变化率最小的平稳 滑翔弹道概念，并给出平稳滑翔弹道的高度、弹道倾角、纵向加速的解析解和弹道动态特性 分析，从而设计了一种控制平滑、弹道振荡收敛速度快、鲁棒性好的反馈控制方法。

本发明一种高超声速飞行器平稳滑翔再入弹道设计方法，包括以下几个步骤：

步骤1：动力学模型建立

假设地球为非旋转均质圆球，建立再入纵向动力学方程，得到了再入飞行器的高度、速 度、弹道倾角和射程对时间的一阶导数；

步骤2：平稳滑翔弹道提出

以攻角α和倾侧角σ为规划变量进行弹道规划，选取为目标函数，所获得的 再入弹道为平稳滑翔弹道，其中a_ε为纵向加速度，为纵向加速度变化率；T为 时间常数；

步骤3：平稳滑翔弹道解析解求解

通过对高度一阶导数再求导建立再入高度的二阶动态微分方程，考虑纵向加速 度的影响，采用正则摄动的方法逐步逼近真实的平稳滑翔高度，然后将平稳滑翔高度解析解 对时间求导，可获得平稳滑翔弹道倾角和纵向加速度的解析解；

步骤4：平稳滑翔弹道动态特性分析

定义并分析平稳滑翔高度增量，将其带入再入高度的二阶动态微分方程，可获得平稳滑 翔弹道的动态特性，包括稳定性分析，自然频率和阻尼分析；

步骤5：平稳滑翔反馈控制设计

为了使振荡的弹道快速收敛到平稳滑翔状态，对平稳滑翔高度增量动态系统采用微分反 馈进行极点配置，抑制弹道震荡，同时给出了定阻尼下的微分反馈系数。

其中，在步骤1中所述的“假设地球为非旋转均质圆球，建立再入纵向动力学方程，得 到了再入飞行器的高度、速度、弹道倾角和射程对时间的一阶导数”，其具体作法如下：

假设地球为非旋转均质圆球，建立再入纵向动力学方程如下所示，

$\stackrel{·}{V}=-D-g\sin \gamma ---\left(2\right)$

$\stackrel{·}{\gamma }=\frac{1}{V}[L_1+(V^2/r-g)\cos \gamma ]---\left(3\right)$

$\stackrel{·}{s}=\mathrm{rV}\cos \gamma /R_0---\left(4\right)$

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

其中，在步骤2中所述的“以攻角α和倾侧角σ为规划变量进行弹道规划，选取为目标函数，所获得的再入弹道为平稳滑翔弹道，”，其具体作法如下：

以攻角α和倾侧角σ为规划变量进行弹道规划，选取为目标函数，所获得的 再入弹道为平稳滑翔弹道，其中a_ε为纵向加速度，为纵向加速度变化率；T为 时间常数。它具有如下特点：

(1)当且仅当再入弹道的攻角曲线和倾侧角曲线均连续可导时，才能获得对应的平稳滑 翔弹道，不连续或者不可导的攻角曲线(或者倾侧角曲线)均只能得到振荡的再入弹道。

(2)平稳滑翔弹道是给定攻角曲线和倾侧角曲线所对应的再入弹道族的中心，其它弹道 均围绕平稳滑翔弹道振荡。

其中，在步骤3中所述的“通过对高度一阶导数再求导建立再入高度的二阶动态微分方 程，考虑纵向加速度的影响，采用正则摄动的方法逐步逼近真实的平稳滑翔高度。然后将平 稳滑翔高度解析解对时间求导，可获得平稳滑翔弹道倾角和纵向加速度的解析解”，其具体作 法如下：

A.高度动态微分方程建立

式(1)给出了高度对时间的一阶倒数，将其对时间进一步求导，带入式(2)和式(3)，并考 虑γ为小量，cosγ≈1且g(sinγ)²≈0。取可得，

$\stackrel{··}{h}+(D/V)\stackrel{·}{h}-L_1+g-V^2/r=0---\left(5\right)$

式(5)中，V为长周期项，将其视为时变系数，则获得了再入高度的二阶动态微分方程。

B.基于正则摄动的平稳滑翔高度解析解求解

考虑a_ε的影响，相对于h_sg来说，a_ε带来的高度偏差Δh_eq为一小量，因此可采用正则摄 动的方法逐步逼近真实的平稳滑翔高度。引入小参数ε，则式(5)可化为如下形式，

$k_c{C}_{L1}^{*}{V}^2{e}^{-{\beta }_{r}h}-g+\frac{V^2}{r}=ϵ(A_1\stackrel{··}{h}+A_2\stackrel{·}{h})---\left(6\right)$

式中，ε为一常数小量，A₁＝1/ε，为纵向升 阻比的倒数，为升力系数纵向分量，β_r为指数大气模型常数，通常取β_r＝1/7200。

设平稳滑翔高度为h_sg，根据正则摄动法则，可将其写为如下所示，

h_sg(t,ε)＝h_e0+εh_e1+ε³h_e2+...            (7)

式中，h_e0、h_e1和h_e2分别为平稳滑翔高度的零阶项、一阶项和二阶项。将式(7)带入式 (6)可得，

$\left(\begin{array}{c}{k}_c{C}_{L1}^{*}{V}^2{e}^{-{\beta }_r{h}_{e0}}\{1-{\beta }_r[ϵh_{e1}+{ϵ}^2{h}_{e2}+...]+{\beta }_r^2{[ϵh_{e1}+{ϵ}^2{h}_{e2}+...]}^2/2+...\}-g+V^2/r\\ =A_1ϵ[{\stackrel{··}{h}}_{e0}+ϵ{\stackrel{··}{h}}_{e1}+{ϵ}^2{\stackrel{··}{h}}_{e2}+...]+A_2ϵ[{\stackrel{·}{h}}_{e0}+ϵ{\stackrel{·}{h}}_{e1}+{ϵ}^2{\stackrel{·}{h}}_{e2}+...]\}\end{array}\right)---\left(8\right)$

式中，和分别为平稳滑翔高度的零阶项、一阶项和二阶项的一阶导数，和分别为平稳滑翔高度的零阶项、一阶项和二阶项的二阶导数。

由上式可得，平稳滑翔高度的零阶项即为平衡滑翔高度，即h_e0＝h_eq。

高度的一阶项为，

$h_{e1}=-\frac{A_1{\stackrel{··}{h}}_{e0}+A_2{\stackrel{·}{h}}_{e0}}{(g-V^2/r){\beta }_r}---\left(9\right)$

高度的二阶项为，

$h_{e2}=\frac{{\beta }_r{h}_{e1}^2}{2}-\frac{A_1{\stackrel{··}{h}}_{e1}+A_2{\stackrel{·}{h}}_{e1}}{{\beta }_r(g-V^2/r)}---\left(10\right)$

C.平稳滑翔弹道倾角和纵向加速度解析解求解

通过对平稳滑翔高度解析解求导，并认为和f_K保持不变，忽略g sinγ和a_εsg对的 影响，可获得平稳滑翔弹道倾角和纵向加速度的解析解，如下所示，

${\gamma }_{\mathrm{sg}}=\frac{-2g{f}_K}{V^2{\beta }_r}---\left(11\right)$

$a_{\mathrm{ϵsg}}=\frac{4g{f}_K^2(g-V^2/r)}{V^2{\beta }_r}---\left(12\right)$

式中，γ_sg和a_εsg分别为平稳滑翔弹道倾角和纵向加速度。

其中，在步骤4中所述的“定义并分析平稳滑翔高度增量，将其带入再入高度的二阶动 态微分方程，可获得平稳滑翔弹道的动态特性，包括稳定性分析，自然频率和阻尼分析”，其 具体作法如下：

A.稳定性分析

定义平稳滑翔高度增量如下，

Δh＝h-h_sg                 (13)

通过分析Δh，可获得平稳滑翔弹道的动态特性。将式(13)带入式(5)可得，

$({\stackrel{··}{h}}_{\mathrm{sg}}+\Delta \stackrel{··}{h})+\frac{f_K{L}_1^{*}{e}^{-{\beta }_r\mathrm{\Delta h}}}{V}({\stackrel{·}{h}}_{\mathrm{sg}}+\Delta \stackrel{·}{h})-L_1^{*}{e}^{-{\beta }_r\mathrm{\Delta h}}+g-V^2/r=0---\left(14\right)$

式中，为平稳滑翔高度增量的二阶导数，为平稳滑翔纵向升力加速度。将阶展开，并忽略高度二阶项，可得，

$\Delta \stackrel{··}{h}+\frac{f_K{L}_1^{*}}{V}\Delta \stackrel{·}{h}+L_1^{*}{\beta }_r(1-f_K{\gamma }_{\mathrm{sg}})\mathrm{\Delta h}=0---\left(15\right)$

式(15)为Δh的二阶线性微分方程，通过参数V的根轨迹分析Δh的自然特性，可知系统 稳定；随着速度的增大，稳定性降低，即速度越大弹道越容易振荡。

B.自然频率和阻尼分析

利用式(15)，并忽略a_εsg和γ_sg等小量的影响，可得平稳滑翔高度增量振荡的自然频率和 阻尼分别如下所示，

${\omega }_n=\sqrt{{\beta }_r(g-V^2/r)}---\left(16\right)$

$\zeta =\frac{f_K}{2}\sqrt{\frac{g/V^2-1/r}{{\beta }_r}}---\left(17\right)$

可知，平稳滑翔高度增量振荡的自然频率仅与当前的飞行速度相关，而阻尼则与速度及 纵向升阻比相关。

其中，在步骤5中所述的“为了使振荡的弹道快速收敛到平稳滑翔状态，对平稳滑翔高 度增量动态系统采用微分反馈进行极点配置，抑制弹道震荡，同时给出了定阻尼下的微分反 馈系数”，其具体作法如下：

A.反馈系数选取

对平稳滑翔高度增量动态系统的极点进行配置。在传统的再入跟踪制导律中，通常采用 比例反馈(Δh或者ΔD)和微分反馈(或者Δγ)来配置系统的极点，将它们引入平稳滑 翔高度增量动态系统后，动态方程变为，

$\Delta \stackrel{··}{h}+\frac{f_K{L}_1^{*}}{V}\Delta \stackrel{·}{h}+L_1^{*}{\beta }_r(\frac{\cos \sigma }{\cos {\sigma }_c}-f_K{\gamma }_{\mathrm{sg}})\mathrm{\Delta h}=-k_1\mathrm{\Delta h}-k_2\Delta \stackrel{·}{h}---\left(18\right)$

式中，k₁和k₂为反馈系数；σ_c为期望倾侧角；σ为实际倾侧角。

比例反馈并不能加快振荡弹道的收敛，却需要较大的反馈加速度。因此在再入弹道跟踪 时可仅采用微分反馈来抑制弹道振荡，即k₁＝0。

B.定阻尼微分反馈系数

需要根据当前速度实时调整反馈系数，设阻尼取定值，则反馈系数k₂如下，

$k_2=-(f_K{L}_1^{*}/V+2{\zeta }_c^2{\beta }_r\mathrm{V\Delta \gamma })+\sqrt{4{\zeta }_c^2{\beta }_r[{\zeta }_c^2{\beta }_r{V}^2\Delta {\gamma }^2+L_1^{*}+L_1^{*}{f}_K(\mathrm{\Delta \gamma }-{\gamma }_{\mathrm{sg}})]}---\left(19\right)$

式中，ζ_c为指定阻尼，通常取ζ_c＝0.707。

本发明的优点在于：

(1)提出了满足纵向加速度变化率最小的平稳滑翔弹道概念，且平稳滑翔弹道是所有相 同攻角和倾侧角再入弹道族的中心，其余弹道均围绕它振荡。

(2)建立了平稳滑翔高度动态微分方程，并获得了平稳滑翔弹道高度、弹道倾角和纵向 加速度的解析解，利用解析解可快速获得期望升力系数和阻力系数对应的平稳滑翔弹道。

(3)再入过程中，平稳滑翔高度增量动态系统是稳定的，并给出了较为简洁的弹道振荡 自然频率和阻尼表达式，其中自然频率仅与当前的飞行速度相关，而阻尼则还与纵向升阻比 成反比，均与实际弹道仿真结果相吻合。

(4)证明了纯微分反馈是抑制再入滑翔弹道振荡的最佳方案，具有控制平滑、弹道振荡 收敛速度快、鲁棒性好等特点。

(5)将定阻尼微分反馈用于引入段弹道的生成，简化了再入滑翔弹道的纵平面制导。

附图说明

图1是攻角连续性对平稳滑翔弹道的影响，

图1(a)是再入攻角曲线，

图1(b)是再入弹道曲线，

图1(c)是再入弹道倾角曲线，

图1(d)是再入纵向加速度曲线；

图2是不同初始高度下的滑翔段弹道的高度偏差和弹道倾角偏差，

图2(a)是高度偏差，

图2(b)是弹道倾角偏差；

图3是平衡滑翔高度与平稳滑翔高度之差；

图4是平稳滑翔高度解析解误差分析；

图5是平稳滑翔弹道倾角数值解与解析解，

图5(a)是算例1弹道倾角数值解与解析解，

图5(b)是算例2弹道倾角数值解与解析解，

图5(c)是算例3弹道倾角数值解与解析解,

图5(d)是平稳滑翔弹道倾角解析解相对误差；

图6是平稳滑翔纵向加速度数值解与解析解，

图6(a)是算例1纵向加速度数值解与解析解，

图6(b)是算例2纵向加速度数值解与解析解，

图6(c)是算例3纵向加速度数值解与解析解,

图6(d)是平稳滑翔纵向加速度解析解相对误差；

图7是平稳滑翔高度增量动态系统结构；

图8是参数V的根轨迹；

图9是初始高度对自然频率的影响(定攻角)，

图9(a)是不同初始高度下弹道曲线，

图9(b)是不同初始高度下自然频率曲线；

图10是初始弹道倾角对自然频率的影响(定攻角)，

图10(a)是不同初始弹道倾角下弹道曲线，

图10(b)是不同初始弹道倾角下自然频率曲线；

图11是攻角振荡对自然频率的影响，

图11(a)是震荡攻角曲线，

图11(b)是跳跃滑翔和平稳滑翔弹道曲线，

图11(c)是跳跃滑翔和平稳滑翔弹道高度差，

图11(d)是自然频率对比；

图12是不同方法获得的自然频率；

图13是倾侧角对弹道阻尼的影响，

图13(a)是不同倾侧角下的弹道震荡曲线，

图13(b)是不同倾侧角下的弹道阻尼曲线；

图14是k1和k2对参数根轨迹的影响；

图15是不同初始高度下k₁Δh和所对应的反馈加速度比较；

图16是不同反馈系数下的弹道(小扰动)；

图17是不同反馈系数下的高度差(小扰动)；

图18是不同反馈系数下的倾侧角(小扰动)；

图19是k1的影响分析(大扰动)；

图20是存在初始状态偏差的滑翔弹道，

图20(a)是初始高度偏差下的滑翔弹道，

图20(b)是初始弹道倾角偏差下的滑翔弹道；

图21是引入段的攻角曲线；

图22是引入段的弹道曲线；

图23是引入段的热流密度曲线；

图24是最大热流密度比较；

图25是基于定阻尼反馈的仿真结果，

图25(a)是倾侧角曲线，

图25(b)是攻角曲线；

图26是基于定阻尼微分反馈的引入段和滑翔段弹道，

图26(a)是再入弹道曲线，

图26(b)是热流密度曲线；

图27本发明所述方法流程图。

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

α、γ、σ、h、s、t、V、a_ε分别是弹道的攻角、弹道倾角、倾侧角、高度、射程、飞行时 间、速度、纵向加速度；Δh、Δγ分别是实际弹道与平稳滑翔弹道飞行高度差、弹道倾角差， 分别是飞行高度差的一阶、二阶导数；h₀、γ₀分别是初始高度、初始弹道倾角；ω_n、ξ 分别是自然频率和弹道阻尼；L₁^*为平稳滑翔纵向升力加速度,γ_sg平稳滑翔弹道倾角,β_r为指 数大气模型常数,f_K为纵向升阻比的倒数,Δa是反馈加速度，k₁、k₂是反馈系数，σ_c是期望 倾侧角，dQ是热流密度，dQ_max是最大热流密度。

具体实施方式

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

本发明通过分析给定攻角曲线和倾侧角曲线再入弹道族的特点，提出以纵向加速度导数 的平方的积分来衡量弹道的平滑程度，并以此为基础定义了平稳滑翔弹道。然后对平稳滑翔 弹道的解析解和弹道动态特性进行分析，从而获得了快速生成平稳滑翔弹道的方法。

本发明是一种高超声速飞行器平稳滑翔再入弹道设计方法，见图27所示，它包括以下几 个步骤：

步骤1：动力学模型建立

假设地球为非旋转均质圆球，建立再入纵向动力学方程如下所示，

$\stackrel{·}{V}=-D-g\sin \gamma ---\left(2\right)$

$\stackrel{·}{\gamma }=\frac{1}{V}[L_1+(V^2/r-g)\cos \gamma ]---\left(3\right)$

$\stackrel{·}{s}=\mathrm{rV}\cos \gamma /R_0---\left(4\right)$

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

$\left(\begin{array}{cc}{L}_1=\frac{\rho V^{2}S{C}_{L1}}{2m}& D=\frac{\rho V^{2}S{C}_D}{2m}\end{array},---\left(5\right)\right)$

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

$\rho ={\rho }_0{e}^{-{\beta }_{r}h}---\left(6\right)$

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

步骤2：定义平稳滑翔弹道

以攻角α和倾侧角σ为规划变量进行弹道规划，选取为目标函数，所获得的 再入弹道为平稳滑翔弹道，其中a_ε为纵向加速度，为纵向加速度变化率；T为 时间常数。它具有如下特点：

(1)当且仅当再入弹道的攻角曲线和倾侧角曲线均连续可导时，才能获得对应的平稳滑 翔弹道，不连续或者不可导的攻角曲线(或者倾侧角曲线)均只能得到振荡的再入弹道(如 图1(a)、(b)、(c)、(d)所示，图中σ＝0deg)。

(2)平稳滑翔弹道是给定攻角曲线和倾侧角曲线所对应的再入弹道族的中心，其它弹道 均围绕平稳滑翔弹道振荡。图2(a)、(b)给出了不同初始高度下的滑翔段弹道的高度偏差和 弹道倾角偏差。图中攻角为15deg，倾侧角为0deg；Δh＝h-h_sg，Δγ＝γ-γ_sg，h和γ分别为 弹道实际高度和弹道倾角，γ_sg为平稳滑翔弹道倾角。由图2(a)、(b)可以看出，Δh和Δγ均 在零附近振荡，并且振幅逐渐减小，证明了上述结论。

步骤3：平稳滑翔弹道解析解求解

A.高度动态微分方程建立

式(1)给出了高度对时间的一阶倒数，将其对时间进一步求导，并带入式(2)和式(3)，可 得，

再入过程中，γ为小量，因此上式中cosγ≈1且g(sinγ)²≈0。取带入式 (7)，可得，

$\stackrel{··}{h}+(D/V)\stackrel{·}{h}-L_1+g-V^2/r=0---\left(8\right)$

式(8)中，V为长周期项，将其视为时变系数，则获得了再入高度的二阶动态微分方程。 对比式(8)和式(3)还可得，

$a_{ϵ}=\stackrel{··}{h}+(D/V)\stackrel{·}{h}---\left(9\right)$

设为纵向升阻比的倒数，则阻力可写为，

D＝f_K(g+a_ε-V²/r)       (10)

将式(5)和式(10)带入式(8)可得高度动态微分方程

$k_c{C}_{L1}^{*}{V}^2{e}^{-{\beta }_{r}h}-g+\frac{V^2}{r}=\stackrel{··}{h}+\frac{f_K(g+a_{ϵ}-V^2/r)}{V}\stackrel{·}{h}---\left(11\right)$

式中，k_c＝ρ₀S/(2m)。

B.基于正则摄动的平稳滑翔高度解析解求解

为了实现对平稳滑翔弹道的快速规划，需要对平稳滑翔弹道的高度、弹道倾角和纵向加 速度进行解析求解。平稳滑翔弹道的纵向加速度a_ε为小量，采用平衡滑翔条件来估算平稳滑 翔弹道的高度。平衡滑翔条件假定a_ε＝0，其表达式如下，

$k_c{C}_{L1}^{*}{V}^2{e}^{-{\beta }_r{h}_{\mathrm{eq}}}-g+V^2/r=0---\left(12\right)$

式中，h_eq为平衡滑翔高度，由式(14)可得其表达式为，

$h_{\mathrm{eq}}=-\ln [(g/V^2-1/r)/\left(k_c{C}_{L1}^{*}\right)]/{\beta }_r---\left(13\right)$

然而，由于平衡滑翔条件忽略了a_ε的影响，使得平衡滑翔高度与平稳滑翔高度之间存在 一定的偏差，如图3所示。图中Δh_eq＝h_eq-h_sg，h_sg为平稳滑翔高度。可以看出当速度较大 时，两者之差较小，采用式(13)估算平稳滑翔高度具有较高的精度；当速度较小时，两者之 差较大，采用式(13)估算平稳滑翔高度精度较差。除此之外，高超声速飞行器的纵向升阻比 也影响式(13)估算精度，纵向升阻比越小则估算精度越低。

为了获得精度更高的平稳滑翔高度解析解，需要进一步考虑a_ε的影响。由上面的分析知， 相对于h_sg来说，a_ε带来的高度偏差Δh_eq为一小量，因此可采用正则摄动的方法逐步逼近真实 的平稳滑翔高度。引入小参数ε，则式(11)可化为如下形式，

$k_c{C}_{L1}^{*}{V}^2{e}^{-{\beta }_{r}h}-g+\frac{V^2}{r}=ϵ(A_1\stackrel{··}{h}+A_2\stackrel{·}{h})---\left(14\right)$

式中，ε为一常数小量，A₁＝1/ε，A₂＝A₁f_K(g+a_ε-V²/r)/V。

设平稳滑翔高度为h_sg，根据正则摄动法则，可将其写为如下所示，

h_sg(t,ε)＝h_e0+εh_e1+ε³h_e2+...            (15)

式中，h_e0、h_e1和h_e2分别为平稳滑翔高度的零阶项、一阶项和二阶项。将式(15)带入式 (14)可得，

$\left(\begin{array}{c}{k}_c{C}_{L1}^{*}{V}^2{e}^{-{\beta }_r{h}_{e0}}\{1-{\beta }_r[ϵh_{e1}+{ϵ}^2{h}_{e2}+...]+{\beta }_r^2{[ϵh_{e1}+{ϵ}^2{h}_{e2}+...]}^2/2+...\}-g+V^2/r\\ =A_1ϵ[{\stackrel{··}{h}}_{e0}+ϵ{\stackrel{··}{h}}_{e1}+{ϵ}^2{\stackrel{··}{h}}_{e2}+...]+A_2ϵ[{\stackrel{·}{h}}_{e0}+ϵ{\stackrel{·}{h}}_{e1}+{ϵ}^2{\stackrel{·}{h}}_{e2}+...]\}\end{array}\right)---\left(16\right)$

式中，和分别为平稳滑翔高度的零阶项、一阶项和二阶项的一阶导数，和分别为平稳滑翔高度的零阶项、一阶项和二阶项的二阶导数。

由上式可得，平稳滑翔高度的零阶项即为平衡滑翔高度，即h_e0＝h_eq，一阶项和二阶项则 分别满足如下关系，

$-k_c{C}_{L1}^{*}{V}^2{e}^{-{\beta }_r{h}_{e0}}{\beta }_r{h}_{e1}=A_1{\stackrel{··}{h}}_{e0}+A_2{\stackrel{·}{h}}_{e0}---\left(17\right)$

$k_c{C}_{L1}^{*}{V}^2{e}^{-{\beta }_r{h}_{e0}}({\beta }_r^2{h}_{e1}^2/2-{\beta }_r{h}_{e2})=A_1{\stackrel{··}{h}}_{e1}+A_2{\stackrel{·}{h}}_{e1}---\left(18\right)$

将h_e0对时间求一阶导数可得，

${\stackrel{·}{h}}_{e0}=\frac{1}{{\beta }_r}[f_a+f_V]\stackrel{·}{V}---\left(19\right)$

式中，f_a和f_V分别为速度相关系数，$f_a=(\partial C_{L1}^{*}/\partial V)/C_{L1}^{*},f_V=2g/(\mathrm{gV}-V^3/r).$

进一步将式(19)对时间求导可得，

${\stackrel{··}{h}}_{e0}=\frac{1}{{\beta }_r}[f_a+f_V]\stackrel{··}{V}+\frac{1}{{\beta }_r}[\frac{\partial f_a}{\partial V}+\frac{\partial f_V}{\partial V}]{\stackrel{·}{V}}^2---\left(20\right)$

式中，和分别为f_a和f_V对速度的偏导数。式(19)和式(20)中还包含着速 度的一阶导数和二阶导数它们分别如下所示，

$\stackrel{·}{V}=-\frac{g+a_{ϵ}-V^2/r}{K_1^{*}}-\frac{g}{V}{\stackrel{·}{h}}_{e0}---\left(21\right)$

$\stackrel{··}{V}=\frac{2V\stackrel{·}{V}}{r{K}_1^{*}}+\frac{(g+a_{ϵ1}-V^2/r)}{{\left(K_1^{*}\right)}^2}\frac{\partial K_1^{*}}{\partial V}\stackrel{·}{V}+(g/V^2){\stackrel{·}{h}}_{e0}\stackrel{·}{V}-(g/V){\stackrel{··}{h}}_{e0}---\left(22\right)$

式中，为给定纵向升阻比，为给定纵向升阻比对速度的偏导数。

将式(12)、式(19)和式(20)带入式(17)，可得高度的一阶项为，

$h_{e1}=-\frac{A_1{\stackrel{··}{h}}_{e0}+A_2{\stackrel{·}{h}}_{e0}}{(g-V^2/r){\beta }_r}---\left(23\right)$

将h_e1分别对时间求一阶导数和二阶导数，并带入式(18)，可得高度的二阶项如下所示，

$h_{e2}=\frac{{\beta }_r{h}_{e1}^2}{2}-\frac{A_1{\stackrel{··}{h}}_{e1}+A_2{\stackrel{·}{h}}_{e1}}{{\beta }_r(g-V^2/r)}---\left(24\right)$

图4给出了不同阶数的平稳滑翔高度解析解与数值方法获得的平稳滑翔高度之间的误差 (仿真攻角为15deg，倾侧角为0deg)。由图可知，仅考虑高度零阶项时，解析解与数值解之 间的最大误差可达40m；在考虑一阶项后，最大误差降低到1m左右；但二阶项仅使最大误差 减小0.2m；因此，在实际应用时可只考虑到高度一阶项。

C.平稳滑翔弹道倾角和纵向加速度解析解求解

通过对平稳滑翔高度解析解求导，可获得平稳滑翔弹道倾角和纵向加速度的解析解，如 下所示，

${\gamma }_{\mathrm{sg}}={\stackrel{·}{h}}_{\mathrm{sg}}/V---\left(25\right)$

$a_{\mathrm{ϵsg}}=\frac{{\stackrel{··}{h}}_{\mathrm{sg}}+(g-V^2/r)f_K{\gamma }_{\mathrm{sg}}}{1-f_K{\gamma }_{\mathrm{sg}}}---\left(26\right)$

式中，γ_sg和a_εsg分别为平稳滑翔弹道倾角和纵向加速度；和分别为平稳滑翔高度 的一阶导数和二阶导数，如下所示，

$\left(\begin{array}{cc}{\stackrel{·}{h}}_{\mathrm{sg}}={\stackrel{·}{h}}_{e0}+ϵ{\stackrel{·}{h}}_{e1}& {\stackrel{··}{h}}_{\mathrm{sg}}={\stackrel{··}{h}}_{e0}\end{array},+ϵ{\stackrel{··}{h}}_{e1}---\left(27\right)\right)$

在式(25)和式(26)中，若和f_K保持不变，并忽略g sinγ和a_εsg对的影响，则γ_sg和 a_εsg可简化为如下形式，

${\gamma }_{\mathrm{sg}}=\frac{-2g{f}_K}{V^2{\beta }_r}---\left(28\right)$

$a_{\mathrm{ϵsg}}=\frac{4g{f}_K^2(g-V^2/r)}{V^2{\beta }_r}---\left(29\right)$

图5(a)、(b)、(c)、(d)和图6(a)、(b)、(c)、(d)分别给出了平稳滑翔弹道倾角和纵 向加速度的解析解，‘X’表示解析解，线条表示数值解。图中仿真算例的攻角曲线和倾侧角 曲线如表1所示，

表1

编号 攻角(deg) 倾侧角(deg) 算例1 α＝17-0.001V σ＝0 算例2 α＝15 σ＝0 算例3 α＝8+0.001V σ＝0 算例4 α＝15+5sin[0.0008π(V-2000)] σ＝0

由图5(a)、(b)、(c)、(d)和图6(a)、(b)、(c)、(d)可以看出，两者的解析解均与数 值解非常接近，均随速度的增大而增大，其中平稳滑翔弹道倾角解析解的相对误差小于1％， 而平稳滑翔纵向加速度解析解的相对误差则小于2％。

步骤4：平稳滑翔弹道动态特性分析

A.稳定性分析

为了获得平稳滑翔弹道，飞行高度和弹道倾角必须分别满足式(15)和式(25)，然而实际 再入过程中滑翔段的起始高度和弹道倾角很可能不满足上述关系，并且飞行过程中扰动也会 使弹道偏离平稳滑翔状态，从而使得弹道振荡。定义平稳滑翔高度增量如下，

Δh＝h-h_sg                    (30)

通过分析Δh，可获得平稳滑翔弹道的动态特性。将式(30)带入式(8)可得，

$({\stackrel{··}{h}}_{\mathrm{sg}}+\Delta \stackrel{··}{h})+\frac{f_K{L}_1^{*}{e}^{-{\beta }_r\mathrm{\Delta h}}}{V}({\stackrel{·}{h}}_{\mathrm{sg}}+\Delta \stackrel{·}{h})-L_1^{*}{e}^{-{\beta }_r\mathrm{\Delta h}}+g-V^2/r=0---\left(31\right)$

式中，为平稳滑翔高度增量的二阶导数，为平稳滑翔纵向升力加速度，如下所示，

$L_1^{*}=g+a_{\mathrm{ϵsg}}-V^2/r---\left(32\right)$

在式(31)中，将一阶展开，并忽略高度二阶项，可得，

$\Delta \stackrel{··}{h}+\frac{f_K{L}_1^{*}}{V}\Delta \stackrel{·}{h}+L_1^{*}{\beta }_r(1-f_K{\gamma }_{\mathrm{sg}})\mathrm{\Delta h}=0---\left(33\right)$

式(33)为Δh的二阶线性微分方程，其系统结构如图7所示。由于和 均仅与速度相关，因此可通过参数V的根轨迹分析Δh的自然特性(如图8所示)。由 图可以看出，整个再入过程中(2000m/s<V<7000m/s)根轨迹均位于虚轴左半边，说明系统稳 定；随着速度的增大，根轨迹逐渐向虚轴靠拢，稳定性降低，即速度越大弹道越容易振荡。

B.自然频率和阻尼分析

利用式(33)，可得平稳滑翔高度增量振荡的自然频率和阻尼分别如下所示，

${\omega }_n=\sqrt{{L_1^{*}\beta }_r(1-f_K{\gamma }_{\mathrm{sg}})}---\left(34\right)$

$\zeta =\frac{f_K{L}_1^{*}}{2V\sqrt{L_1^{*}{\beta }_r(1-f_K{\gamma }_{\mathrm{sg}})}}---\left(35\right)$

将式(32)带入式(34)和式(35)，并忽略a_εsg和γ_sg等小量的影响，可得，

${\omega }_n=\sqrt{{\beta }_r(g-V^2/r)}---\left(36\right)$

$\zeta =\frac{f_K}{2}\sqrt{\frac{g/V^2-1/r}{{\beta }_r}}---\left(37\right)$

由式(36)和式(37)可知，平稳滑翔高度增量振荡的自然频率仅与当前的飞行速度相关， 而阻尼则与速度及纵向升阻比相关。

图9(a)、(b)和图10(a)、(b)分别给出了不同初始高度和初始弹道倾角(算例2) 下的平稳滑翔高度增量的自然频率。可以看出，两者对自然频率的影响均非常小，实际弹道 仿真获得的自然频率与式(36)给出的解析结果几乎完全重合。图11(a)、(b)、(c)、(d)则 分析了攻角(算例4)的变化对自然频率的影响，可以看出它对自然频率的影响很小。图12 进一步将本方法结果与Harpold在航天飞机再入制导律中获得的再入弹道振荡自然频率解析 解(式38)进行了对比，可以看出两者几乎完全重合，均与实际弹道高度增量的振荡频率相 吻合，但本方法给出的结论更为简洁。

${\omega }_{n·\mathrm{Harpold}}^2=3{\stackrel{·}{D}}_{\mathrm{ref}}(\frac{{\stackrel{·}{D}}_{\mathrm{ref}}}{D_{\mathrm{ref}}^2}-\frac{1}{V})+\frac{4D_{\mathrm{ref}}^2}{V^2}-\beta (\frac{V^2}{r}-g)-\frac{{\stackrel{··}{D}}_{\mathrm{ref}}}{D_{\mathrm{ref}}}-\frac{{\stackrel{·}{C}}_{\mathrm{Dref}}^2}{C_{\mathrm{Dref}}^2}+\frac{{\stackrel{··}{C}}_{\mathrm{Dref}}}{C_{\mathrm{Dref}}}---\left(38\right)$

式中，D_ref为参考阻力，分别为参考阻力的一阶导数、二阶导数。ω_n·Harpold为 航天飞行再入制导律所获得的再入弹道振荡自然频率。

图13(a)、(b)给出了阻尼曲线，由于无法从弹道获得阻尼大小，因此只能进行定性分 析。由图可以看出，相同攻角下，倾侧角越大，则对应弹道的高度振荡收敛越快，说明阻尼 越大，与式(39)给出的结论相符。

步骤5：平稳滑翔反馈控制设计

A.反馈系数选取

由于再入弹道的阻尼非常小，使得弹道在偏离平稳滑翔的弹道振荡不易收敛。为了使振 荡的弹道快速收敛到平稳滑翔状态，需要对平稳滑翔高度增量动态系统的极点进行配置。在 传统的再入跟踪制导律中，通常采用比例反馈(Δh或者ΔD)和微分反馈(或者Δγ)来 配置系统的极点，将它们引入平稳滑翔高度增量动态系统后，动态方程变为，

$\Delta \stackrel{··}{h}+\frac{f_K{L}_1^{*}}{V}\Delta \stackrel{·}{h}+L_1^{*}{\beta }_r(\frac{\cos \sigma }{\cos {\sigma }_c}-f_K{\gamma }_{\mathrm{sg}})\mathrm{\Delta h}=-k_1\mathrm{\Delta h}-k_2\Delta \stackrel{·}{h}---\left(39\right)$

式中，k₁和k₂为反馈系数；σ_c为期望倾侧角；σ为实际倾侧角，满足如下关系，

$\cos \sigma =[1-(k_1\mathrm{\Delta h}+k_2\Delta \stackrel{·}{h})/L_1^{*}]\cos {\sigma }_c---\left(40\right)$

图14(a)、(b)、(c)、(d)给出了k₁和k₂对式(39)给出的平稳滑翔高度增量动态系统参 数根轨迹的影响。可以看出，改变k₁会使根轨迹上下移动，主要影响系统的自然频率；改变k₂则使得根轨迹左右移动，主要影响系统的阻尼。由图(a)可知，在k₁＝0时，改变k₂也能使 根轨迹明显远离虚轴，能够明显改变系统的阻尼；由图(b)可知，增大k₁会降低k₂的调节 效率；由图(c)、(d)可以看出，改变k₁不会改变根轨迹与虚轴的相对位置，反而会减小系 统的阻尼。

图15(a)、(b)比较了k₁Δh和所对应的反馈加速度，可以看出要远远小于k₁Δh。 图16至图18给出了不同反馈系数下的弹道仿真结果，由图16可以看出，在倾侧反转(小扰 动)下，弹道会在平稳滑翔高度附近振荡。由图17和图18可知，当k₂＝0时，高度差几乎 不收敛，并且倾侧角振荡厉害；而当k₂＝0.03时，高度差快速收敛到零，倾侧角也快速收敛 到期望值。图19(a)、(b)则进一步分析了k₁的影响，由图可以看出，k₁取0和0.001时， 弹道的收敛速度几乎一致，但k₁＝0对应的倾侧角的变化更为平缓。

总之，比例反馈并不能加快振荡弹道的收敛，却需要较大的反馈加速度。因此在再入弹 道跟踪时可仅采用微分反馈来抑制弹道振荡，即k₁＝0。

B.定阻尼微分反馈系数

当再入高度增量动态系统的阻尼为0.4至0.8时，系统的调节时间较短、超调量较小。 但对于给定的k₂，系统的阻尼随速度变化较大，因此需要根据当前速度实时调整反馈系数。 设阻尼取定值，则反馈系数k₂如下，

$k_2=-(f_K{L}_1^{*}/V+2{\zeta }_c^2{\beta }_r\mathrm{V\Delta \gamma })+\sqrt{4{\zeta }_c^2{\beta }_r[{\zeta }_c^2{\beta }_r{V}^2\Delta {\gamma }^2+L_1^{*}+L_1^{*}{f}_K(\mathrm{\Delta \gamma }-{\gamma }_{\mathrm{sg}})]}---\left(41\right)$

式中，ζ_c为指定阻尼，通常取ζ_c＝0.707。将上式带入式(40)，即可获得抑制弹道振荡 所需的倾侧角。若由式(40)获得的倾侧角超出其可行范围，则需要通过调整攻角来获得所需 要的反馈加速度。

图20(a)、(b)给出了存在初始状态偏差情况下采用定阻尼微分反馈获得的滑翔段弹道， 图中s_togo为剩余射程。可以看出，不论初始状态如何，在定阻尼微分反馈的作用下，弹道均 快速收敛到平稳滑翔状态。

实施例：

将定阻尼微分反馈应用于引入段制导，引入段的起点和终点的速度、高度和弹道倾角均 给定，其中伪谱法以引入段射程最大为目标函数。由图21可知，采用定阻尼微分反馈方法获 得攻角曲线较为光滑，并且在引入段的高度快速下降阶段保持攻角最大，这样能最大限度的 降低引入段的最大热流密度(如图23所示)。由图22可知，采用定阻尼微分反馈方法获得的 引入段弹道更为光滑，但射程略微小于伪谱法结果。图24比较了采用定阻尼微分反馈方法获 得的引入段的最大热流密度与飞行器最大攻角飞行(飞行器极限能力)获得的最大热流密度， 可以看出两者在最大热流密度较大时高度重合，证明定阻尼微分反馈方法能够满足引入段的 最大热流密度约束要求，从而实现了引入段和滑翔段的制导控制方法的统一。

图25(a)、(b)和图26(a)、(b)给出了基于定阻尼微分反馈的引入段和滑翔段弹道仿 真结果。由图25(a)、(b)可以看出，在经过引入段的调整后，攻角和倾侧角均快速收敛到 期望值。由图26(a)、(b)可以看出，再入弹道只在引入段存在着由于飞行器气动力不足带 来的“波谷”，之后的滑翔段弹道均非常光滑，并且热流密度在“波谷”达到峰值之后一直递 减，满足约束要求。