CN105549387A - 一种助推段广义标控脱靶量解析制导方法 - Google Patents

Abstract

一种助推段广义标控脱靶量解析制导方法，它有五大步骤：步骤1：助推段制导问题建模，包括标控弹道运动建模和制导修正一阶摄动建模和终端约束建模；步骤2：助推段标控弹道正则摄动求解建模；步骤3：基于正则摄动的助推段标控弹道零阶项解析求解；步骤4：基于正则摄动的助推段标控弹道一阶项解析求解；步骤5：基于标控脱靶量的助推段最优制导指令求解。本发明以零攻角为助推段弹道的标准控制，并采用正则摄动方法进行解析求解，从而获得了解析的标控脱靶量；最后采用最优控制方法求解了修正标控脱靶量的制导指令，从而获得了一种满足助推段任务需求的非线性解析最优制导方法。

Description

一种助推段广义标控脱靶量解析制导方法

技术领域

本发明提供一种助推段广义标控脱靶量解析制导方法，属于航天技术、武器技术领域。

背景技术

助推段弹道通常可分为大气层内的程序飞行段和大气层外的制导飞行段，其中程序飞行段采用依据时间或者速度设定的开环控制指令，通常不进行反馈调整；制导飞行段则采用摄动制导或者迭代制导进行闭环修正。但对于助推—滑翔高超声速飞行器来说，助推段的中的高度越高则再入拉起段的最大热流密度越大。这就要求助推段降低终端高度，从而使得助推段弹道的大部分或者全部位于大气层内。当前大气层内的闭环制导方法通常采用最优控制理论进行建模，并采用数值优化算法进行求解，如修正打靶法、配点法、最小值算法等。但上述方法均存在着计算量大、难于在线应用等问题。

零控脱靶量是研究制导问题的一个重要概念，狭义的零控脱靶量指的是在导弹和目标均不进行机动的条件下，两者之间可达到的最小距离；广义的零控脱靶量则还包含终端角度误差、终端速度误差等。在经典的制导方法中，制导指令通常为零控脱靶量的负反馈，如比例导引、广义显式制导方法、分段最优制导方法等。尽管零控脱靶量的概念在解决含终端位置约束和角度约束的制导问题上取得了很大的成功，当随着制导问题的非线性程度增加和约束的增多，零控脱靶量已经不能完全满足制导方法设计的需求。

广义标控脱靶量是在零控脱靶量的基础上发展出来的一个制导相关概念，指的是飞行器按照给定的标准控制指令飞行所得轨迹的终端状态与约束之间的偏差，用于解决非线性系统的制导问题。由于广义标控脱靶量通常需要通过求解非线性微分方程获得，目前多采用数值积分的方法进行求解。在获得广义标控脱靶量之后，既可以通过构造相关的负反馈对当前的标准控制指令进行修正，又可以采用优化算法对标准控制指令进行更新。

发明内容

本发明的目的是提供一种助推段广义标控脱靶量解析制导方法。该方法以零攻角为助推段弹道的标准控制；并采用正则摄动方法进行解析求解，从而获得了解析的标控脱靶量；最后采用最优控制方法求解了修正标控脱靶量的制导指令，从而获得了一种满足助推段任务需求的非线性解析最优制导方法。

本发明为一种助推段广义标控脱靶量解析制导方法，包括以下几个步骤：

步骤1：助推段制导问题建模，包括标控弹道运动建模和制导修正一阶摄动建模和终端约束建模。忽略地球自转，航迹坐标系下的助推器纵向运动方程如下所示：

$\left\{\begin{array}{c}\stackrel{\·}{V}=(P\cos \mathrm{\α}-D)/m-g\sin \mathrm{\γ}\\ \stackrel{\·}{\mathrm{\γ}}=(P\sin \mathrm{\α}+L)/\left(mV\right)+\[V/(R_0+h)-g/V\]\cos \mathrm{\γ}\\ \stackrel{\·}{h}=V\sin \mathrm{\γ}\\ \stackrel{\·}{m}=q_m\end{array}\right.$

上式中，V、γ、h和m分别为助推器的速度、弹道倾角、高度和质量；和分别为助推器速度、弹道倾角、高度和质量对时间的导数；P为助推器的推力，与发动机的质量秒流量和当地大气压强相关；α为助推器的攻角；L、D分别为助推器的升力和阻力，与飞行器的攻角、马赫数及动压相关；q_m为助推器发动机的质量秒流量；g为重力加速度，与助推器高度相关；R₀为地球半径。

助推段通常采用零攻角作为标准控制，此时有L＝0；若进一步假定质量秒流量q_m为常数，则可得标控弹道的运动模型为：

$\left\{\begin{array}{c}{\stackrel{\·}{V}}_b=(P_b-D_b)/(m_0-q_{m}t)-g_b(e^{x_b}-1)/(e^{x_b}+1)\\ {\stackrel{\·}{x}}_b=2\[V_b/(R_0+h_b)-g_b/V_b\]\\ {\stackrel{\·}{h}}_b=V_b(e^{x_b}-1)/(e^{x_b}+1)\end{array}\right.$

上式中，m₀和t分别为助推器的初始质量和飞行时间；下标‘b’表征标控弹道，V_b、h_b、g_b、P_b和D_b分别为标控弹道的速度、高度、重力加速度、推力和阻力；x_b为与标控弹道的弹道倾角相关的变量，表达式如下：

x_b＝ln[(1+sinγ_b)/(1-sinγ_b)]

上式中，γ_b为标控标控弹道的弹道倾角。

将运动方程在标控弹道附近求解一阶摄动，并忽略速度摄动项和弹道倾角的自相关项，可得制导修正一阶摄动方程为：

$\left\{\begin{array}{c}\mathrm{\Δ}\stackrel{\·}{\mathrm{\γ}}=f_{ba}{a}_n\\ \mathrm{\Δ}\stackrel{\·}{h}=f_{bh}\mathrm{\Δ}\mathrm{\γ}\end{array}\right.---\left(1\right)$

上式中，Δγ和Δh分别为助推段的弹道倾角增量和高度增量，和分别为它们对时间的导数；f_ba和f_bh均为由标控弹道确定的时变参数；a_n为助推器的法向加速度，由推力法向分量和升力共同决定。

此外，助推段终端高度和弹道倾角还需满足约束值，即终端摄动修正量需满足如下关系：

$\left\{\begin{array}{c}{\mathrm{\Δ\γ}}_f={\mathrm{\γ}}_{cf}-{\mathrm{\γ}}_{bf}\\ {\mathrm{\Δh}}_f=h_{cf}-h_{bf}\end{array}\right.---\left(2\right)$

上式中，γ_cf和h_cf分别为由制导任务确定的助推段终端弹道倾角和高度；γ_bf和h_bf分别为标控弹道的终端弹道倾角和高度；Δγ_f和Δh_f为期望终端弹道倾角增量和高度增量，它们共同构成了助推段的广义标控脱靶量。

步骤2：助推段标控弹道正则摄动求解建模。包含助推段标控弹道速度、弹道倾角及高度的零阶微分方程和一阶微分方程，分别如下：

$\left\{\begin{array}{c}{\stackrel{\·}{V}}_b^{\left(0\right)}=\stackrel{\‾}{P}/(m_0-q_{m}t)\\ {\stackrel{\·}{x}}_b^{\left(0\right)}=2\[V_b^{\left(0\right)}/(R_0+h_0)-g_i/V_b^{\left(0\right)}\]\\ {\stackrel{\·}{h}}_b^{\left(0\right)}=V_b^{\left(0\right)}\[{\mathrm{sin\γ}}_0+c_2(x_b^{\left(0\right)}-x_0)\]\end{array}\right.$

$\left\{\begin{array}{c}\mathrm{\ϵ}{\stackrel{\·}{V}}_b^{\left(1\right)}=a_{\mathrm{\ϵ}}^{\left(0\right)}\\ \mathrm{\ϵ}{\stackrel{\·}{x}}_b^{\left(1\right)}=k_{\mathrm{\γ}1}^{\left(0\right)}{\mathrm{\ϵV}}_b^{\left(1\right)}+k_{\mathrm{\γ}2}^{\left(0\right)}\\ \mathrm{\ϵ}{\stackrel{\·}{h}}_b^{\left(1\right)}=k_{r1}^{\left(0\right)}\left({\mathrm{\ϵV}}_b^{\left(1\right)}\right)+k_{r2}^{\left(0\right)}\left({\mathrm{\ϵx}}_b^{\left(1\right)}\right)+k_{r3}^{\left(0\right)}\end{array}\right.$

上式中，上标‘(0)’和‘(1)’分别表征正则摄动的零阶项和一阶项；和分别为助推段标控弹道速度、弹道倾角相关变量和高度的零阶项，和分别为它们对时间的一阶导数；和分别为助推段标控弹道速度、弹道倾角相关变量和高度的一阶项，和分别为它们对时间的一阶导数。ε为正则摄动符号常数；为与初始推力、阻力和弹道倾角相关的常数；h₀、γ₀、x₀和g_i分别为助推段标控弹道的初始高度、弹道倾角、弹道倾角相关变量和重力加速度，c₂为与它们相关的常数。和分别为与助推段标控弹道零阶项相关的系数。

步骤3：基于正则摄动的助推段标控弹道零阶项解析求解。对上一步得到的助推段正则摄动零阶微分方程解析积分，可得速度、弹道倾角相关变量和高度的零阶项的解析解分别如下：

$\left\{\begin{array}{c}{V}_b^{\left(0\right)}=V_0-V_e\ln (1-t/\mathrm{\τ})\\ x_b^{\left(0\right)}=x_0+\frac{2c_1}{V_e(R_0+h_0)}{f}_{\mathrm{\γ}1}-\frac{2c_1{g}_i}{V_e}{f}_{\mathrm{\γ}2}\\ h_b^{\left(0\right)}=h_0+c_1\[e^{-V_0/V_e}(V_0+V_e){\mathrm{sin\γ}}_0-e^{-V_b/V_e}(V_b+V_e)({\mathrm{sin\γ}}_0-c_2{c}_0+c_2{x}_b^{\left(0\right)})\]+\\ \frac{c_1^2{c}_2}{V_e}(\frac{1}{R_0+h_0}{f}_{h1}+\frac{V_e}{R_0+h_0}{f}_{h2}+\frac{g_i{V}_e}{2}(e^{-2V_b/V_e}-e^{-2V_0/V_e})-g_i{V}_e{f}_{h3})\end{array}\right.$

上式中，V₀、x₀和h₀分别为助推段标控弹道的初始速度、弹道倾角相关变量和高度；V_e为与平均推力和质量秒流量确定的常数；τ为与助推器初始质量和质量秒流量相的常数；c₁为与V_e、τ和初始速度相关的常数；f_γ1、f_γ2、f_h1、f_h2和f_h3均为与助推段正则摄动速度零阶项相关的变量。

步骤4：基于正则摄动的助推段标控弹道一阶项解析求解。利用上一步中得到助推段标控弹道零阶项的解析解，对助推段标控弹道一阶摄动微分方程中的系数进行多项式拟合，最后进行积分可得速度、弹道倾角相关变量和高度的一阶项的解析解分别如下：

${\mathrm{\ϵV}}_b^{\left(1\right)}=\frac{c_{v3}}{4}{t}^4+\frac{c_{v2}}{3}{t}^3+\frac{c_{v1}}{2}{t}^2+c_{v0}t$

$\begin{array}{c}{\mathrm{\ϵx}}_b^{\left(1\right)}=\left(\frac{c_{v3}{c}_{\mathrm{\γ}12}}{28}\right)t^7+(\frac{c_{v3}{c}_{\mathrm{\γ}11}}{24}+\frac{c_{v2}{c}_{\mathrm{\γ}12}}{18})t^6+(\frac{c_{v3}{c}_{\mathrm{\γ}10}}{20}+\frac{c_{v2}{c}_{\mathrm{\γ}11}}{15}+\frac{c_{v1}{c}_{\mathrm{\γ}12}}{10})t^5+\\ (\frac{c_{v2}{c}_{\mathrm{\γ}10}}{12}+\frac{c_{v1}{c}_{\mathrm{\γ}11}}{8}+\frac{c_{v0}{c}_{\mathrm{\γ}12}}{4})t^4+(\frac{c_{v1}{c}_{\mathrm{\γ}10}}{6}+\frac{c_{v0}{c}_{\mathrm{\γ}11}+c_{\mathrm{\γ}22}}{3})t^3+\left(\frac{c_{v0}{c}_{\mathrm{\γ}10}+c_{\mathrm{\γ}21}}{2}\right)t^2+c_{\mathrm{\γ}20}t\end{array}$

${\mathrm{\ϵh}}_b^{\left(1\right)}=c_{10h}{t}^{10}+c_{9h}{t}^9+c_{8h}{t}^8+c_{7h}{t}^7+c_{6h}{t}^6+c_{5h}{t}^5+c_{4h}{t}^4+c_{3h}{t}^3+c_{2h}{t}^2+c_{1h}t$

上式中，c_v3、c_v2、c_v1和c_v0为速度相关参数的三次多项式拟合系数；c_γ10、c_γ11、c_γ12、c_γ20、c_γ21、c_γ22为弹道倾角相关参数的二次多项式拟合系数；c_1h、c_2h、c_3h、c_4h、c_5h、c_6h、c_7h、c_8h、c_9h、c_10h为高度相关参数的拟合系数。

步骤5：基于标控脱靶量的助推段最优制导指令求解。将步骤3和步骤4获得的助推段标控弹道解析解带入式(2)可以获得广义标控脱靶量，同时对式(1)中的f_ba和f_bh进行多项式拟合，并以下式为求解助推段最优制导指令的目标函数，

$J={\∫}_{t_r}^{t_f}\frac{a_n^2}{2(t_f-t)}dt$

最终可得助推段的最优制导指令为，

$a_n=\frac{f_{r2f}{\mathrm{\Δh}}_f-f_{r4f}{\mathrm{\Δ\γ}}_f}{f_{5f}{f}_{r2f}-f_{r3f}{f}_{r4f}}{V}_r{t}_{go}$

上式中，a_n为助推段法向加速度，为制导指令；f_r2f、f_r3f、f_r4f和f_r5f分别为t＝t_f时f_r2、f_r3、f_r4和f_r5的取值；t_f为助推段的终端时刻；V_r为助推段的当前速度；t_go为助推段的剩余飞行时间，t_go＝t_f-t。

本发明的优点在于：

(1)提出了以零攻角作为助推器制导修正段的标准控制，不但简化了标控弹道的运动方程，而且使得制导修正所需的攻角位于零附近，降低了助推器制导修正带来的速度损失。

(2)采用正则摄动方法解析求解了助推段纵向弹道，获得了速度、弹道倾角和高度的高精度解析解。

(3)给出了一种基于广义标控脱靶量的解析最优制导方法，能够同时满足助推段终端高度和弹道倾角约束，同时使得终端攻角收敛到零。

附图说明

图1是本发明流程图。

图2a是系数的二次多项式拟合结果示意图。

图2b是系数的二次多项式拟合结果示意图。

图2c是系数的二次多项式拟合结果示意图。

图2d是系数的二次多项式拟合结果示意图。

图2e是系数的二次多项式拟合结果示意图。

图2f是系数的二次多项式拟合结果示意图。

图3是标控弹道速度解析解精度校验示意图。

图4是标控弹道的弹道倾角解析解精度校验示意图。

图5是标控弹道的高度解析解精度校验示意图。

图6a是标控弹道速度解析解绝对误差示意图。

图6b是标控弹道速度解析解相对误差示意图。

图7a是标控弹道的弹道倾角解析解绝对误差示意图。

图7b是标控弹道的弹道倾角解析解相对误差示意图。

图8a是标控弹道高度解析解绝对误差示意图。

图8b是标控弹道高度解析解相对误差示意图。

图9a是助推段纵向弹道打靶曲线示意图。

图9b是助推段弹道倾角打靶曲线示意图。

图9c是助推段速度打靶曲线示意图。

图9d是助推段攻角打靶曲线示意图。

图9e是助推段法向过载打靶曲线示意图。

图9f是助推段动压打靶曲线示意图。

图10a助推段打靶终端高度和弹道倾角偏差散布示意图。

图10b助推段打靶终端速度和攻角散布示意图。

图中符号说明如下：

n_z代表由升力产生的法向过载；q代表动压，单位为kPa；

为修正轴向加速度；和分别为与助推段标控弹道零阶项相关的系数。

具体实施方式

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

针对高超声速飞行器助推段的制导问题，本发明以零攻角零攻角作为助推器制导修正段的标准控制，并对标控弹道运动方程、制导修正一阶摄动方程和终端约束进行了建模；然后采用正则摄动方法对标控弹道进行解析求解，获得了速度、弹道倾角和高度的零阶项和一阶项的解析解；最后采用最优控制方法获得了基于标控脱靶量的最优制导指令。整个发明流程如图1所示。

本发明是一种助推段广义标控脱靶量解析制导方法，包括以下几个步骤：

步骤1：助推段制导问题建模，包括标控弹道运动建模和制导修正一阶摄动建模和终端约束建模。

忽略地球自转，航迹坐标系下的助推器纵向运动方程如下所示：

$\left\{\begin{array}{c}\stackrel{\·}{V}=(P\cos \mathrm{\α}-D)/m-g\sin \mathrm{\γ}\\ \stackrel{\·}{\mathrm{\γ}}=(P\sin \mathrm{\α}+L)/\left(mV\right)+\[V/(R_0+h)-g/V\]\cos \mathrm{\γ}\\ \stackrel{\·}{h}=V\sin \mathrm{\γ}\\ \stackrel{\·}{m}=q_m\end{array}\right.$

上式中，V、γ、h和m分别为助推器的速度、弹道倾角、高度和质量；和分别为助推器速度、弹道倾角、高度和质量对时间的导数；α为助推器的攻角；R₀为地球半径；g为重力加速度，与助推器高度相关；q_m为助推器发动机的质量秒流量；P为助推器的推力，与发动机的质量秒流量和当地大气压强相关，表达式如下：

P＝I_spq_mg₀-P_aS_N

上式中，I_sp为助推器发动机比冲；S_N为发动机尾喷管面积；g₀为地面重力加速度；P_a为当地大气压强，与飞行高度相关。L、D分别为助推器的升力和阻力，表达式如下：

L＝ρV²SC_L/2D＝ρV²SC_D/2

上式中，C_L和C_D分别为升力系数和阻力系数；ρ为当地大气密度，与飞行高度相关；S为助推器气动参考面积。

由于助推器通常为轴对称飞行器，其零升攻角通常为0deg，若助推段采用零攻角作为标准控制，则有Psinα_b+L_b＝0,α_b＝0deg，为标控攻角；L_b为标准控制对应的升力，从而可得标控弹道的运动模型为：

$\left\{\begin{array}{c}{\stackrel{\·}{V}}_b=(P_b-D_b)/(m_0-q_{m}t)-g_b(e^{x_b}-1)/(e^{x_b}+1)\\ {\stackrel{\·}{x}}_b=2\[V_b/(R_0+h_b)-g_b/V_b\]\\ {\stackrel{\·}{h}}_b=V_b(e^{x_b}-1)/(e^{x_b}+1)\end{array}\right.---\left(3\right)$

上式中，m₀和t分别为助推器的初始质量和飞行时间；下标‘b’表征标控弹道，V_b、h_b、g_b、P_b和D_b分别为标控弹道的速度、高度、重力加速度、推力和阻力；x_b为与标控弹道的弹道倾角相关的变量，表达式如下：

x_b＝ln[(1+sinγ_b)/(1-sinγ_b)]

上式中，γ_b为标控标控弹道的弹道倾角。

将运动方程在标控弹道附近求解一阶摄动，并忽略速度摄动项和弹道倾角的自相关项，可得制导修正一阶摄动方程为：

$\left\{\begin{array}{c}\mathrm{\Δ}\stackrel{\·}{\mathrm{\γ}}=f_{ba}{a}_n\\ \mathrm{\Δ}\stackrel{\·}{h}=f_{bh}\mathrm{\Δ}\mathrm{\γ}\end{array}\right.---\left(4\right)$

上式中，Δγ和Δh分别为助推段的弹道倾角增量和高度增量，和分别为它们对时间的导数；f_ba和f_bh均为由标控弹道确定的时变参数，表达式分别如下：

f_ba＝1/V_bf_bh＝V_bcosγ_b

a_n为助推器的法向加速度，由推力法向分量和升力共同决定，表达式如下：

a_n＝(Psinα+L)/m

此外，助推段终端高度和弹道倾角还需满足约束值，即终端摄动修正量需满足如下关系：

$\left\{\begin{array}{c}{\mathrm{\Δ\γ}}_f={\mathrm{\γ}}_{cf}-{\mathrm{\γ}}_{bf}\\ {\mathrm{\Δh}}_f=h_{cf}-h_{bf}\end{array}\right.---\left(5\right)$

上式中，γ_cf和h_cf分别为由制导任务确定的助推段终端弹道倾角和高度；γ_bf和h_bf分别为标控弹道的终端弹道倾角和高度；Δγ_f和Δh_f为期望终端弹道倾角增量和高度增量，它们共同构成了助推段的广义标控脱靶量。

步骤2：助推段标控弹道正则摄动求解建模。

将标控弹道的速度微分方程等号右边分为两个部分，如下所示：

${\stackrel{\·}{V}}_b=a_{\stackrel{\‾}{P}}+a_{\mathrm{\ϵ}}---\left(6\right)$

上式中，为由平均推力产生的轴线加速度，a_ε为其余轴线力产生的摄动修正轴向加速度，它们的表达式分别如下：

$a_{\stackrel{\‾}{P}}=\stackrel{\‾}{P}/(m_0-q_{m}t)---\left(7\right)$

$a_{\mathrm{\ϵ}}=\frac{S_N(P_{a0}-P_a)+D_0-D_b+m_0{g}_i{\mathrm{sin\γ}}_0}{m_0-q_{m}t}-g_b\left(\frac{e^{x_b}-1}{e^{x_b}+1}\right)---\left(8\right)$

上式中，P_a0、D₀和γ₀分别为助推段标控弹道的初始大气压、零升阻力和弹道倾角，g_i为初始重力加速度；为与初始推力、阻力和弹道倾角相关的常数，表达式如下：

$\stackrel{\‾}{P}=g_0{I}_{sp}{q}_m-S_N{P}_{a0}-D_0-m_0{g}_i{\mathrm{sin\γ}}_0$

由于助推段的推力通常远大于重力的轴向分量和阻力的变化，即有，因此可采用正则摄动的方法对标控弹道进行解析求解。

设V_b、x_b和h_b的表达式分别如下，

$V_b=V_b^{\left(0\right)}+{\mathrm{\ϵV}}_b^{\left(1\right)}---\left(9\right)$

$x_b=x_b^{\left(0\right)}+{\mathrm{\ϵx}}_b^{\left(1\right)}---\left(10\right)$

$h_b=h_b^{\left(0\right)}+{\mathrm{\ϵh}}_b^{\left(1\right)}---\left(11\right)$

上式中，上标‘(0)’和‘(1)’分别表征正则摄动的零阶项和一阶项；和分别为助推段标控弹道速度、弹道倾角相关变量和高度的零阶项；和分别为助推段标控弹道速度、弹道倾角相关变量和高度的一阶项。将式(9)、(10)和(11)带入式(3)，可得，

${\stackrel{\·}{V}}_b^{\left(0\right)}+\mathrm{\ϵ}{\stackrel{\·}{V}}_b^{\left(1\right)}=a_{\stackrel{\‾}{P}}-\mathrm{\ϵ}\frac{a_{\mathrm{\ϵ}}}{A}---\left(12\right)$

${\stackrel{\·}{x}}_b^{\left(0\right)}+\mathrm{\ϵ}{\stackrel{\·}{x}}_b^{\left(1\right)}=\[\frac{V_b^{\left(0\right)}+{\mathrm{\ϵV}}_b^{\left(1\right)}}{R_0+h_b^{\left(0\right)}+{\mathrm{\ϵh}}_b^{\left(1\right)}}-\frac{g_b}{(V_b^{\left(0\right)}+{\mathrm{\ϵV}}_b^{\left(1\right)})}\]---\left(13\right)$

${\stackrel{\·}{h}}_b^{\left(0\right)}+\mathrm{\ϵ}{\stackrel{\·}{h}}_b^{\left(1\right)}=(V_b^{\left(0\right)}+{\mathrm{\ϵV}}_b^{\left(1\right)})\left(\frac{e^{x_b^{\left(0\right)}+{\mathrm{\ϵx}}_b^{\left(1\right)}}-1}{e^{x_b^{\left(0\right)}+{\mathrm{\ϵx}}_b^{\left(1\right)}}+1}\right)---\left(14\right)$

上式中，A＝1/ε，为正则摄动求解符号常数；和分别为助推段标控弹道速度、弹道倾角相关变量和高度的零阶项对时间的一阶导数；和分别为助推段标控弹道速度、弹道倾角相关变量和高度的一阶项对时间的一阶导数。

由式(14)可得速度的正则摄动零阶项和一阶修正项分别满足如下微分方程：

${\stackrel{\·}{V}}_b^{\left(0\right)}=\stackrel{\‾}{P}/(m_0-qt)---\left(15\right)$

$\mathrm{\ϵ}{\stackrel{\·}{V}}_b^{\left(1\right)}=a_{\mathrm{\ϵ}}^{\left(0\right)}---\left(16\right)$

上式中，为将正则摄动零阶状态带入式(8)计算得到的摄动修正轴向加速度。

将式(13)的等号右边对ε进行一阶泰勒展开，并且忽略高度变化的影响，可得弹道倾角相关变量的正则摄动零阶项和一阶项的微分方程分别如下：

${\stackrel{\·}{x}}_b^{\left(0\right)}=2\[\frac{V_b^{\left(0\right)}}{R_0+h_0}-\frac{g_i}{V_b^{\left(0\right)}}\]---\left(17\right)$

$\mathrm{\ϵ}{\stackrel{\·}{x}}_b^{\left(1\right)}=k_{\mathrm{\γ}1}^{\left(0\right)}{\mathrm{\ϵV}}_b^{\left(1\right)}+k_{\mathrm{\γ}2}^{\left(0\right)}---\left(18\right)$

上式中，h₀和g_i分别为助推段标控弹道的初始高度和重力加速度；和均为与助推段标控弹道零阶项相关的系数，表达式如下：

$k_{\mathrm{\γ}1}^{\left(0\right)}=\frac{2}{R_0+h_b^{\left(0\right)}}+\frac{2g_b^{\left(0\right)}}{{\left(V_b^{\left(0\right)}\right)}^2}$

$k_{\mathrm{\γ}2}^{\left(0\right)}=2(\frac{V_b^{\left(0\right)}}{R_0+h_b^{\left(0\right)}}-\frac{V_b^{\left(0\right)}}{R_0+h_0}-\frac{g_b^{\left(0\right)}-g_i}{V_b^{\left(0\right)}})$

由于助推器制导段的弹道倾角通常较小，因此可对式(14)中的弹道倾角相关项作如下拟合，

$\frac{e^{x_b^{\left(0\right)}}-1}{e^{x_b^{\left(0\right)}}+1}\≈{\mathrm{sin\γ}}_0+c_2(x_b^{\left(0\right)}-x_0)---\left(19\right)$

上式中，x₀为助推段标控弹道的初始弹道倾角相关变量；c₂为拟合常数，表达式如下：

$c_2=\frac{e^{x_{bf}^{\left(0\right)}-1-(e^{x_{bf}^{\left(0\right)}}+1){\mathrm{sin\γ}}_0}}{(x_{bf}^{\left(0\right)}-x_0^{\left(0\right)})(e^{x_{bf}^{\left(0\right)}}+1)}$

上式中，为助推段标控弹道的弹道倾角相关变量零阶项终值。将式(19)带入式(14)可得，

${\stackrel{\·}{h}}_b^{\left(0\right)}=V_b^{\left(0\right)}\[{\mathrm{sin\γ}}_0+c_2(x_b^{\left(0\right)}-x_0)\]---\left(20\right)$

$\mathrm{\ϵ}{\stackrel{\·}{h}}_b^{\left(1\right)}=k_{r1}^{\left(0\right)}\left({\mathrm{\ϵV}}_b^{\left(1\right)}\right)+k_{r2}^{\left(0\right)}\left({\mathrm{\ϵx}}_b^{\left(1\right)}\right)+k_{r3}^{\left(0\right)}---\left(21\right)$

上式中，和均为与助推段标控弹道零阶项相关的系数，表达式如下：

$k_{h1}^{\left(0\right)}=(e^{x_b^{\left(0\right)}}-1)/(e^{x_b^{\left(0\right)}}+1)$

$k_{h2}^{\left(0\right)}=2e^{x_b^{\left(0\right)}}{V}_b^{\left(0\right)}/{(e^{x_b^{\left(0\right)}}+1)}^2$

$k_{h3}^{\left(0\right)}=V_b^{\left(0\right)}\{\frac{e^{x_b^{\left(0\right)}}-1}{e^{x_b^{\left(0\right)}}+1}-\[{\mathrm{sin\γ}}_0+c_2(x_b^{\left(0\right)}-x_0)\]\}$

步骤3：基于正则摄动的助推段标控弹道零阶项解析求解。

利用步骤2得到的助推段正则摄动零阶项的微分方程，分别进行解析积分，即可获得速度、弹道倾角相关变量和高度的零阶项解析解。首先，对式(15)积分可得，

$V_b^{\left(0\right)}=V_0-V_{e}ln(1-t/\mathrm{\τ})---\left(22\right)$

上式中，V₀为助推段标控弹道的初始速度；τ和V_e均为助推器发动机相关的常数，表达式如下：

$\begin{array}{cc}\mathrm{\τ}=m_0/q_m& V_e=\stackrel{\‾}{P}/q_m\end{array}$

将式(22)带入式(15)可得，

${\stackrel{\·}{V}}_b^{\left(0\right)}=\frac{V_e}{c_1{e}^{-V_b^{\left(0\right)}/V_e}}---\left(23\right)$

上式中，c₁为常数，c₁＝τexp(V₀/V_e)。将式(23)带入式(17)可得，

$\frac{{\mathrm{dx}}_b^{\left(0\right)}}{{\mathrm{dV}}_b^{\left(0\right)}}=\frac{2\[V_b^{\left(0\right)}/(R_0+h_0)-g_i/V_b^{\left(0\right)}\]c_1{e}^{-V_b^{\left(0\right)}/V_e}}{V_e}$

对上式进行积分可得，

$x_b^{\left(0\right)}=x_0+\frac{2c_1}{V_e(R_0+h_0)}{f}_{\mathrm{\γ}1}-\frac{2c_1{g}_i}{V_e}{f}_{\mathrm{\γ}2}---\left(24\right)$

上式中，f_γ1和f_γ2均为与助推段正则摄动速度零阶项相关的变量，表达式如下：

$f_{\mathrm{\γ}1}={\∫}_{V_0}^{V_b}{\mathrm{Ve}}^{-V/V_e}dV=V_e^2\[e^{-V_0/V_P}(V_0/V_e+1)-e^{-V_b/V_e}(V_b/V_e+1)\]$

$f_{\mathrm{\γ}2}={\∫}_{V_0}^{V_b}\frac{e^{-V/V_e}}{V}dV=ln\frac{V_b}{V_0}+\underset{k=1}{\overset{\mathrm{\∞}}{\Σ}}\frac{(V_b^k-V_0^k)/{(-V_e)}^k}{k\·k!}$

将式(23)带入式(20)可得，

$\frac{{\mathrm{dh}}_b^{\left(0\right)}}{{\mathrm{dV}}_b^{\left(0\right)}}=\frac{c_1}{V_e}{V}_b^{\left(0\right)}{e}^{-V_b^{\left(0\right)}/V_e}\[{\mathrm{sin\γ}}_0+c_2(x_b^{\left(0\right)}-x_0)\]$

对上式分步积分可得，

$\begin{array}{c}{h}_b^{\left(0\right)}=h_0+c_1\[e^{-V_0/V_e}(V_0+V_e){\mathrm{sin\γ}}_0-e^{-V_b/V_e}(V_b+V_e)({\mathrm{sin\γ}}_0-c_2{x}_0+c_2{x}_b^{\left(0\right)})\]+\\ \frac{c_1^2{c}_2}{V_e}(\frac{1}{R_0+h_0}{f}_{h1}+\frac{V_e}{R_0+h_0}{f}_{h2}+\frac{g_i{V}_e}{2}(e^{-2V_b/V_e}-e^{-2V_0/V_e})-g_i{V}_e{f}_{h3}\end{array}---\left(25\right)$

上式中，f_h1、f_h2和f_h3分别为与助推段正则摄动速度零阶项相关的变量，表达式如下：

$f_{h1}={\∫}_{V_0}^{V_b}{V}^2{e}^{-2V/V_e}dV=\frac{V_e^3}{4}\[(\frac{2V_0^2}{V_e^2}+\frac{2V_0}{V_e}+1)e^{-2V_0/V_e}-(\frac{2V_b^2}{V_e^2}+\frac{2V_b}{V_e}+1)e^{-2V_b/V_e}\]$

$f_{h2}={\∫}_{V_0}^{V_b}{\mathrm{Ve}}^{-2V/V_e}dV=\frac{V_e^2}{4}\[(\frac{2V_0}{V_e}+1)e^{-2V_0/V_e}-(\frac{2V_b}{V_e}+1)e^{-2V_b/V_e}\]$

$f_{h3}={\∫}_{V_0}^{V_b}\frac{e^{-2V/V_e}}{V}dV=\ln \frac{V_b}{V_0}+\underset{k=1}{\overset{\mathrm{\∞}}{\Σ}}\frac{2^k(V_b^k-V_0^k)/{(-V_e)}^k}{k\·k!}$

步骤4：基于正则摄动的助推段标控弹道一阶项解析求解。

将助推段标控弹道的零阶项解析解带入式(16)、(18)和(21)，由于表达式过于复杂，不能直接进行积分获得一阶修正项析解，需要首先对其中的系数进行多项式拟合。由图2a-图2f可知，和近似为时间的二次多项式；而则近似为时间的三次多项式。令t分别为0、t_m和t_f，分别带入式(22)、(24)和(25)，可以求得速度、弹道倾角相关变量和高度的三个取值，将其带入上述系数的表达式，即可得到它们在t＝0,t_m,t_f时间节点上的取值，从而可以拟合得到和的二次多项式。而则还需它在t＝0时的导数信息来进行三次多项式拟合，具体如下，

${\stackrel{\·}{a}}_{\mathrm{\ϵ}0}^{\left(0\right)}=-\frac{S_N}{m_0}\frac{\∂P_{a0}}{\∂h}{\stackrel{\·}{h}}_0-g_i\left({\mathrm{cos\γ}}_0\right){\stackrel{\·}{\mathrm{\γ}}}_0-\frac{2D_0}{m_0}(\frac{\∂{\mathrm{\ρ}}_i}{\∂h}\frac{{\stackrel{\·}{h}}_0}{{\mathrm{\ρ}}_i}+\frac{{\stackrel{\·}{V}}_0}{V_0}+\frac{\∂C_{D0}}{\∂Ma}\frac{{\stackrel{\·}{V}}_0}{C_{D0}c})+\frac{2\mathrm{\μ}}{{(R_0+h_0)}^3}{\stackrel{\·}{h}}_0{\mathrm{sin\γ}}_0$

上式中，ρ_i、C_D0和D₀分别为助推段标控弹道的初始大气密度、阻力系数和阻力；和分别为标控弹道初始大气压强和密度对高度的偏导数；为标控弹道初始阻力系数对马赫数的偏导数；c为声速；Ma为马赫数；和分别为标控弹道初始速度、弹道倾角和高度对时间的导数，可由式(3)求解获得。最终，可以获得如下拟合结果，

$\left\{\begin{array}{c}{a}_{\mathrm{\ϵ}}^{\left(0\right)}=c_{v3}{t}^3+c_{v2}{t}^2+c_{v1}t+c_{v0}\\ k_{\mathrm{\γ}1}^{\left(0\right)}=c_{\mathrm{\γ}12}{t}^2+c_{\mathrm{\γ}11}t+c_{\mathrm{\γ}10}\\ k_{\mathrm{\γ}2}^{\left(0\right)}=c_{\mathrm{\γ}22}{t}^2+c_{\mathrm{\γ}21}t+c_{\mathrm{\γ}20}\\ k_{h1}^{\left(0\right)}=c_{h12}{t}^2+c_{h11}t+c_{h10}\\ k_{h2}^{\left(0\right)}=c_{h22}{t}^2+c_{h21}t+c_{h20}\\ k_{h3}^{\left(0\right)}=c_{h32}{t}^2+c_{h31}t+c_{h30}\end{array}\right.---\left(26\right)$

上式中，c_v3、c_v2、c_v1和c_v0为速度相关参数的三次多项式拟合系数；c_γ10、c_γ11、c_γ12、c_γ20、c_γ21、c_γ22为弹道倾角相关参数的二次多项式拟合系数；c_h10、c_h11、c_h12、c_h20、c_h21、c_h22、c_h30、c_h31和c_h32为高度相关参数的二次多项式拟合系数。图2a-图2f给出了原始数据与拟合结果的对比，可以看出，上述拟合具有较高的精度。

最后，将式(26)给出的多项式拟合结果分别带入式(16)、(18)和(21)，可得积分获得助推段标控弹道一阶项的解析解如下，

${\mathrm{\ϵV}}_b^{\left(1\right)}=\frac{c_{v3}}{4}{t}^4+\frac{c_{v2}}{3}{t}^3+\frac{c_{v1}}{2}{t}^2+c_{v0}t---\left(27\right)$

$\begin{array}{c}{\mathrm{\ϵx}}_b^{\left(1\right)}=\left(\frac{c_{v3}{c}_{\mathrm{\γ}12}}{28}\right)t^7+(\frac{c_{v3}{c}_{\mathrm{\γ}11}}{24}+\frac{c_{v2}{c}_{\mathrm{\γ}12}}{18})t^6+(\frac{c_{v3}{c}_{\mathrm{\γ}10}}{20}+\frac{c_{v2}{c}_{\mathrm{\γ}11}}{15}+\frac{c_{v1}{c}_{\mathrm{\γ}12}}{10})t^5+\\ (\frac{c_{v2}{c}_{\mathrm{\γ}10}}{12}+\frac{c_{v1}{c}_{\mathrm{\γ}11}}{8}+\frac{c_{v0}{c}_{\mathrm{\γ}12}}{4})t^4+(\frac{c_{v1}{c}_{\mathrm{\γ}10}}{6}+\frac{c_{v0}{c}_{\mathrm{\γ}11}+c_{\mathrm{\γ}22}}{3})t^3+\left(\frac{c_{v0}{c}_{\mathrm{\γ}10}+c_{\mathrm{\γ}21}}{2}\right)t^2+c_{\mathrm{\γ}20}t\end{array}---\left(28\right)$

${\mathrm{\ϵh}}_b^{\left(1\right)}=c_{10h}{t}^{10}+c_{9h}{t}^9+c_{8h}{t}^8+c_{7h}{t}^7+c_{6h}{t}^6+c_{5h}{t}^5+c_{4h}{t}^4+c_{3h}{t}^3+c_{2h}{t}^2+c_{1h}t---\left(29\right)$

上式中，c_1h、c_2h、c_3h、c_4h、c_5h、c_6h、c_7h、c_8h、c_9h、c_10h为高度相关参数的拟合系数，它们的表达式分别如下：

$\begin{array}{cc}{c}_{10h}=\frac{c_{v3}{c}_{\mathrm{\γ}12}}{280}{c}_{h22}& c_{9r}=\frac{c_{v3}{c}_{\mathrm{\γ}12}}{252}{c}_{h21}+(\frac{c_{v3}{c}_{\mathrm{\γ}11}}{216}+\frac{c_{v2}{c}_{\mathrm{\γ}12}}{162})c_{h22}\end{array}$

$c_{8r}=\frac{c_{v3}{c}_{y12}}{224}{c}_{h20}+(\frac{c_{v3}{c}_{y11}}{192}+\frac{c_{v2}{c}_{y12}}{144})c_{h21}+(\frac{c_{v3}{c}_{y10}}{160}+\frac{c_{v2}{c}_{y11}}{120}+\frac{c_{v1}{c}_{y12}}{80})c_{h22}$

$\begin{array}{c}{c}_{7r}=(\frac{c_{v3}{c}_{\mathrm{\γ}11}}{168}+\frac{c_{v2}{c}_{\mathrm{\γ}12}}{126})c_{h20}+(\frac{c_{v3}{c}_{\mathrm{\γ}10}}{140}+\frac{c_{v2}{c}_{\mathrm{\γ}11}}{105}+\frac{c_{v1}{c}_{\mathrm{\γ}12}}{70})c_{h21}+\\ (\frac{c_{v2}{c}_{\mathrm{\γ}10}}{84}+\frac{c_{v1}{c}_{\mathrm{\γ}11}}{56}+\frac{c_{v0}{c}_{\mathrm{\γ}12}}{28})c_{h22}+\frac{c_{v3}{c}_{h12}}{28}\end{array}$

$\begin{array}{c}{c}_{6r}=(\frac{c_{v3}{c}_{\mathrm{\γ}10}}{120}+\frac{c_{v2}{c}_{\mathrm{\γ}11}}{90}+\frac{c_{v1}{c}_{\mathrm{\γ}12}}{60})c_{h20}+(\frac{c_{v2}{c}_{\mathrm{\γ}10}}{72}+\frac{c_{v1}{c}_{\mathrm{\γ}11}}{48}+\frac{c_{v0}{c}_{\mathrm{\γ}12}}{24})c_{h21}\\ +(\frac{c_{v1}{c}_{\mathrm{\γ}10}}{36}+\frac{c_{v0}{c}_{\mathrm{\γ}11}+c_{\mathrm{\γ}22}}{18})c_{h22}+\frac{c_{v3}{c}_{h11}}{24}+\frac{c_{v2}{c}_{h12}}{18}\end{array}$

$\begin{array}{c}{c}_{5r}=(\frac{c_{v2}{c}_{\mathrm{\γ}10}}{60}+\frac{c_{v1}{c}_{\mathrm{\γ}11}}{40}+\frac{c_{v0}{c}_{\mathrm{\γ}12}}{20})c_{h20}+(\frac{c_{v1}{c}_{\mathrm{\γ}10}}{30}+\frac{c_{v0}{c}_{\mathrm{\γ}11}+c_{\mathrm{\γ}22}}{15})c_{h21}+\\ \left(\frac{c_{v0}{c}_{\mathrm{\γ}10}+c_{\mathrm{\γ}21}}{10}\right)c_{h22}+\frac{c_{v3}{c}_{h10}}{20}+\frac{c_{v2}{c}_{h11}}{15}+\frac{c_{v1}{c}_{h12}}{10}\end{array}$

$c_{4r}=(\frac{c_{v1}{c}_{\mathrm{\γ}10}}{24}+\frac{c_{v0}{c}_{\mathrm{\γ}11}+c_{\mathrm{\γ}22}}{12})c_{h20}+\frac{(c_{v0}{c}_{\mathrm{\γ}10}+c_{\mathrm{\γ}21})c_{h21}}{8}+\frac{c_{\mathrm{\γ}20}{c}_{h22}}{4}+\frac{c_{v0}{c}_{h12}}{4}+\frac{c_{v2}{c}_{h10}}{12}+\frac{c_{v1}{c}_{h11}}{8}$

$c_{3r}=\frac{(c_{v0}{c}_{\mathrm{\γ}10}+c_{\mathrm{\γ}21})c_{h20}}{6}+\frac{c_{\mathrm{\γ}20}{c}_{h21}}{3}+\frac{c_{v1}{c}_{h10}}{6}+\frac{c_{v0}{c}_{h11}+c_{h32}}{3}$

为了校验上述标控弹道解析解的精度，设助推器的初始高度为35.86km，初始弹道倾角为19.98deg，初始速度为876.9m/s，分别采用数值积分和解析解的方法对标控弹道进行求解，得到的结果如图3至图5所示。由图可以看出，仅采用正则摄动零阶解时，速度和高度的解析解存在着较大的误差；在考虑正则摄动一阶修正项的影响后，速度、高度和弹道倾角的解析解均具有非常高的精度。图6a-图6b至图8a-图8b进一步给出了解析解的误差分析。由图可以看出，制导段速度解析解的相对误差小于0.8％；弹道倾角解析解的相对误差小于0.3％；高度解析解的相对误差小于1％。

步骤5：基于标控脱靶量的助推段最优制导指令求解。

将式(24)、(25)、(28)和(29)获得的助推段标控弹道的弹道倾角和高度解析解带入式(5)，即可获得助推段的广义标控脱靶量。若Δγ_f≠0或者Δh_f≠0，则需要用标控弹道附近的一阶摄动方程求解制导指令。定义制导修正目标函数如下：

$J={\∫}_{t_0}^{t_f}\frac{a_n^2}{2(t_f-t)}dt---\left(30\right)$

上式中，t₀和t_f分别为制导段的初始时刻和终端时刻。由式(5)和(30)可得上述优化问题的哈密尔顿函数为：

$H={\mathrm{\λ}}_1{f}_{ba}{a}_n+{\mathrm{\λ}}_2{f}_{bh}\mathrm{\Δ}\mathrm{\γ}+0.5a_n^2/(t_f-t)---\left(31\right)$

上式中，λ₁和λ₂为协态变量。将上式分别对Δγ、Δh和a_n求偏导，可得协态方程和控制方程分别如下：

$\begin{array}{cc}{\stackrel{\·}{\mathrm{\λ}}}_1=-{\mathrm{\λ}}_2{f}_{bh}& {\stackrel{\·}{\mathrm{\λ}}}_2=0\end{array}---\left(32\right)$

a_n＝-λ₁f_ba(t_f-t)(33)

对式(32)积分，并带入式(33)可得，

a_n＝(λ₂₀f_r1-λ₁₀)f_ba(t_f-t)(34)

上式中，λ₁₀和λ₂₀均为协态变量初值；f_r1为由f_bh积分获得的时间相关函数，如下所示：

$f_{r1}={\∫}_{t_0}^t{f}_{bh}dt$

将式(34)带入(5)中的弹道倾角增量微分方程，积分可得，

Δγ＝λ₂₀f_r2+λ₁₀f_r3(35)

上式中，f_r2和f_r3均为时间相关函数，表达式如下：

$\begin{array}{cc}{f}_{r2}={\∫}_{t_0}^t{f}_{ba}^2{f}_{r1}(t_f-t)dt& f_{r3}=-{\∫}_{t_0}^t{f}_{ba}^2(t_f-t)dt\end{array}$

将式(35)带入式(5)中的高度增量微分方程，积分可得，

Δh＝λ₂₀f_r4+λ₁₀f_r5(36)

上式中，f_r4和f_r5均为时间相关函数，表达式如下：

$\begin{array}{cc}{f}_{r4}={\∫}_{t_r}^t{f}_{bh}{f}_{r2}dt& f_{r5}={\∫}_{t_r}^t{f}_{bh}{f}_{r3}dt\end{array}$

最终，由式(34)、(35)和(36)可得制导所需的控制量为：

$a_n=\frac{f_{r2f}{\mathrm{\Δh}}_f-f_{r4f}{\mathrm{\Δ\γ}}_f}{f_{5f}{f}_{r2f}-f_{r3f}{f}_{r4f}}{V}_0{t}_{go}---\left(37\right)$

上式中，f_r2f、f_r3f、f_r4f和f_r5f分别为t＝t_f时f_r2、f_r3、f_r4和f_r5的取值；t_go为助推段的剩余飞行时间，t_go＝t_f-t₀。

实施案例

为了校验基于广义标控脱靶量的助推段制导方法，采用米诺陶5火箭的第一级和第二级为助推器模型，其中第一级采用程序俯仰角制导，第二级采用基于广义标控脱靶量的助推段制导方法，并且两级之间增加无动力滑翔段。仿真参数设置如表1所示，拉偏参数设置如表2所示。

表1助推段制导仿真参数设置

表2仿真参数拉偏设置

参数	拉偏值(3σ)	参数	拉偏值(3σ)
				升力系数	10％	大气密度	10％
阻力系数	10％	质量秒流量	5％

图9a-图9f给出了助推段打靶结果，由图可以看出，纵向弹道和弹道倾角曲线散布均较小，仅在第二级初始阶段需要较大的攻角，但随着第二级剩余飞行时间的减少，攻角会逐渐收敛到零附近；整个助推段的最大过载小于0.5，最大动压在90kPa左右，且第二级的动压近似为常数。图10(a)、(b)给出了助推段打靶的终端状态分布，表3进一步给出了它们的统计结果。可以看出，终端高度偏差小于5m；终端弹道倾角偏差小于0.05deg；终端攻角小于0.3deg；终端速度散布的标准差约为20m/s。

表3助推段打靶结果统计

Claims (1)

1.一种助推段广义标控脱靶量解析制导方法，其特征在于：它包括以下步骤：

步骤1：助推段制导问题建模，包括标控弹道运动建模和制导修正一阶摄动建模和终端约束建模；忽略地球自转，航迹坐标系下的助推器纵向运动方程如下所示：

$\left\{\begin{array}{c}\stackrel{\·}{V}=(Pcos\mathrm{\α}-D)/m-gsin\mathrm{\γ}\\ \stackrel{\·}{\mathrm{\γ}}=(Psin\mathrm{\α}+L)/\left(mV\right)+\[V/(R_0+h)-g/V\]cos\mathrm{\γ}\\ \stackrel{\·}{h}=V\sin \mathrm{\γ}\\ \stackrel{\·}{m}=q_m\end{array}\right.$

上式中，V、γ、h和m分别为助推器的速度、弹道倾角、高度和质量； 和分别为助推器速度、弹道倾角、高度和质量对时间的导数；P为助推器的推力，与发动机的质量秒流量和当地大气压强相关；α为助推器的攻角；L、D分别为助推器的升力和阻力，与飞行器的攻角、马赫数及动压相关；q_m为助推器发动机的质量秒流量；g为重力加速度，与助推器高度相关；R₀为地球半径；

助推段通常采用零攻角作为标准控制，此时有L＝0；若进一步假定质量秒流量q_m为常数，则得到标控弹道的运动模型为：

$\left\{\begin{array}{c}{\stackrel{\·}{V}}_b=(P_b-D_b)/(m_0-q_{m}t)-g_b(e^{x_b}-1)/(e^{x_b}+1)\\ {\stackrel{\·}{x}}_b=2\[V_b/(R_0+h_b)-g_b/V_b\]\\ {\stackrel{\·}{h}}_b=V_b(e^{x_b}-1)/(e^{x_b}+1)\end{array}\right.$

上式中，m₀和t分别为助推器的初始质量和飞行时间；下标‘b’表征标控弹道，V_b、h_b、g_b、P_b和D_b分别为标控弹道的速度、高度、重力加速度、推力和阻力；x_b为与标控弹道的弹道倾角相关的变量，表达式如下：

x_b＝ln[(1+sinγ_b)/(1-sinγ_b)]

上式中，γ_b为标控标控弹道的弹道倾角；

将运动方程在标控弹道附近求解一阶摄动，并忽略速度摄动项和弹道倾角的自相关项，得制导修正一阶摄动方程为：

$\left\{\begin{array}{c}\mathrm{\Δ}\stackrel{\·}{\mathrm{\γ}}=f_{ba}{a}_n\\ \mathrm{\Δ}\stackrel{\·}{h}=f_{bh}\mathrm{\Δ}\mathrm{\γ}\end{array}\right.---\left(1\right)$

上式中，Δγ和Δh分别为助推段的弹道倾角增量和高度增量，和分别为它们对时间的导数；f_ba和f_bh均为由标控弹道确定的时变参数；a_n为助推器的法向加速度，由推力法向分量和升力共同决定；

此外，助推段终端高度和弹道倾角还需满足约束值，即终端摄动修正量需满足如下关系：

$\left\{\begin{array}{c}{\mathrm{\Δ\γ}}_f={\mathrm{\γ}}_{cf}-{\mathrm{\γ}}_{bf}\\ {\mathrm{\Δh}}_f=h_{cf}-h_{bf}\end{array}\right.---\left(2\right)$

上式中，γ_cf和h_cf分别为由制导任务确定的助推段终端弹道倾角和高度；γ_bf和h_bf分别为标控弹道的终端弹道倾角和高度；Δγ_f和Δh_f为期望终端弹道倾角增量和高度增量，它们共同构成了助推段的广义标控脱靶量；

步骤2：助推段标控弹道正则摄动求解建模；包含助推段标控弹道速度、弹道倾角及高度的零阶微分方程和一阶微分方程，分别如下：

$\left\{\begin{array}{c}{\stackrel{\·}{V}}_b^{\left(0\right)}=\stackrel{\‾}{P}/(m_0-q_{m}t)\\ {\stackrel{\·}{x}}_b^{\left(0\right)}=2\[V_b^{\left(0\right)}/(R_0+h_0)-g_i/V_b^{\left(0\right)}\]\\ {\stackrel{\·}{h}}_b^{\left(0\right)}=V_b^{\left(0\right)}\[sin{\mathrm{\γ}}_0+c_2(x_b^{\left(0\right)}-x_0)\]\end{array}\right.$

$\left\{\begin{array}{c}\mathrm{\ϵ}{\stackrel{\·}{V}}_b^{\left(1\right)}=a_{\mathrm{\ϵ}}^{\left(0\right)}\\ \mathrm{\ϵ}{\stackrel{\·}{x}}_b^{\left(1\right)}=k_{\mathrm{\γ}1}^{\left(0\right)}{\mathrm{\ϵV}}_b^{\left(1\right)}+k_{\mathrm{\γ}2}^{\left(0\right)}\\ \mathrm{\ϵ}{\stackrel{\·}{h}}_b^{\left(1\right)}=k_{r1}^{\left(0\right)}\left({\mathrm{\ϵV}}_b^{\left(1\right)}\right)+k_{r2}^{\left(0\right)}\left({\mathrm{\ϵx}}_b^{\left(1\right)}\right)+k_{r3}^{\left(0\right)}\end{array}\right.$

上式中，上标‘(0)’和‘(1)’分别表征正则摄动的零阶项和一阶项；和分别为助推段标控弹道速度、弹道倾角相关变量和高度的零阶项，和分别为它们对时间的一阶导数；和分别为助推段标控弹道速度、弹道倾角相关变量和高度的一阶项，和分别为它们对时间的一阶导数；ε为正则摄动符号常数；为与初始推力、阻力和弹道倾角相关的常数；h₀、γ₀、x₀和g_i分别为助推段标控弹道的初始高度、弹道倾角、弹道倾角相关变量和重力加速度，c₂为与它们相关的常数；和分别为与助推段标控弹道零阶项相关的系数；

步骤3：基于正则摄动的助推段标控弹道零阶项解析求解；对上一步得到的助推段正则摄动零阶微分方程解析积分，得速度、弹道倾角相关变量和高度的零阶项的解析解分别如下：

$\left\{\begin{array}{c}{V}_b^{\left(0\right)}=V_0-V_e\ln (1-t/\mathrm{\τ})\\ x_b^{\left(0\right)}=x_0+\frac{2c_1}{V_e(R_0+h_0)}{f}_{\mathrm{\γ}1}-\frac{2c_1{g}_i}{V_e}{f}_{\mathrm{\γ}2}\\ h_b^{\left(0\right)}=h_0+c_1\[e^{-V_0/V_e}(V_0+V_e){\mathrm{sin\γ}}_0-e^{-V_b/V_e}(V_b+V_e)({\mathrm{sin\γ}}_0-c_2{x}_0+c_2{x}_b^{\left(0\right)})\]+\\ \frac{c_1^2{c}_2}{V_e}(\frac{1}{R_0+h_0}{f}_{h1}+\frac{V_e}{R_0+h_0}{f}_{h2}+\frac{g_i{V}_e}{2}(e^{-2V_b/V_e}-e^{-2V_0/V_e})-g_i{V}_e{f}_{h3})\end{array}\right.$

上式中，V₀、x₀和h₀分别为助推段标控弹道的初始速度、弹道倾角相关变量和高度；V_e为与平均推力和质量秒流量确定的常数；τ为与助推器初始质量和质量秒流量相的常数；c₁为与V_e、τ和初始速度相关的常数；f_γ1、f_γ2、f_h1、f_h2和f_h3均为与助推段正则摄动速度零阶项相关的变量；

步骤4：基于正则摄动的助推段标控弹道一阶项解析求解；利用上一步中得到助推段标控弹道零阶项的解析解，对助推段标控弹道一阶摄动微分方程中的系数进行多项式拟合，最后进行积分得速度、弹道倾角相关变量和高度的一阶项的解析解分别如下：

${\mathrm{\ϵV}}_b^{\left(1\right)}=\frac{c_{v3}}{4}{t}^4+\frac{c_{v2}}{3}{t}^3+\frac{c_{v1}}{2}{t}^2+c_{v0}t$

$\begin{array}{c}{\mathrm{\ϵx}}_b^{\left(1\right)}=\left(\frac{c_{v3}{c}_{\mathrm{\γ}12}}{28}\right)t^7+(\frac{c_{v3}{c}_{\mathrm{\γ}11}}{24}+\frac{c_{v2}{c}_{\mathrm{\γ}12}}{18})t^6+(\frac{c_{v3}{c}_{\mathrm{\γ}10}}{20}+\frac{c_{v2}{c}_{\mathrm{\γ}11}}{15}+\frac{c_{v1}{c}_{\mathrm{\γ}12}}{10})t^5+\\ (\frac{c_{v2}{c}_{\mathrm{\γ}10}}{12}+\frac{c_{v1}{c}_{\mathrm{\γ}11}}{8}+\frac{c_{v0}{c}_{\mathrm{\γ}12}}{4})t^4+(\frac{c_{v1}{c}_{\mathrm{\γ}10}}{6}+\frac{c_{v0}{c}_{\mathrm{\γ}11}+c_{\mathrm{\γ}22}}{3})t^3+\left(\frac{c_{v0}{c}_{\mathrm{\γ}10}+c_{\mathrm{\γ}21}}{2}\right)t^2+c_{\mathrm{\γ}20}t\end{array}$

${\mathrm{\ϵh}}_b^{\left(1\right)}=c_{10h}{t}^{10}+c_{9h}{t}^9+c_{8h}{t}^8+c_{7h}{t}^7+c_{6h}{t}^6+c_{5h}{t}^5+c_{4h}{t}^4+c_{3h}{t}^3+c_{2h}{t}^2+c_{1h}t$

上式中，c_v3、c_v2、c_v1和c_v0为速度相关参数的三次多项式拟合系数；c_γ10、c_γ11、c_γ12、c_γ20、c_γ21、c_γ22为弹道倾角相关参数的二次多项式拟合系数；c_1h、c_2h、c_3h、c_4h、c_5h、c_6h、c_7h、c_8h、c_9h、c_10h为高度相关参数的拟合系数；

步骤5：基于标控脱靶量的助推段最优制导指令求解；将步骤3和步骤4获得的助推段标控弹道解析解带入式 $\left\{\begin{array}{c}{\mathrm{\Δ\γ}}_f={\mathrm{\γ}}_{cf}-{\mathrm{\γ}}_{bf}\\ {\mathrm{\Δh}}_f=h_{cf}-h_{bf}\end{array}\right.---\left(2\right)$ 获得广义标控脱靶量，同时对式 $\left\{\begin{array}{c}\mathrm{\Δ}\stackrel{\·}{\mathrm{\γ}}=f_{ba}{a}_n\\ \mathrm{\Δ}\stackrel{\·}{h}=f_{bh}\mathrm{\Δ}\mathrm{\γ}\end{array}\right.---\left(1\right)$ 中的f_ba和f_bh进行多项式拟合，并以下式为求解助推段最优制导指令的目标函数，

$J={\∫}_{t_r}^{t_f}\frac{a_n^2}{2(t_f-t)}dt$

最终得助推段的最优制导指令为，

$a_n=\frac{f_{r2f}{\mathrm{\Δh}}_f-f_{r4f}{\mathrm{\Δ\γ}}_f}{f_{5f}{f}_{r2f}-f_{r3f}{f}_{r4f}}{V}_r{t}_{go}$

上式中，a_n为助推段法向加速度，为制导指令；f_r2f、f_r3f、f_r4f和f_r5f分别为t＝t_f时f_r2、f_r3、f_r4和f_r5的取值；t_f为助推段的终端时刻；V_r为助推段的当前速度；t_go为助推段的剩余飞行时间，t_go＝t_f-t。