CN103245257B - 基于Bezier曲线的多约束飞行器导引方法 - Google Patents

Abstract

本发明提供了一种基于Bezier曲线的多约束飞行器导引方法，属于飞行器动力学与制导技术领域。本发明方法根据具体的飞行器运动模型，三阶Bezier曲线方程模拟弹道，用到四个控制点，中间两个控制点用参数k₁和k₂描述，确定飞行器的攻角和倾侧角；通过优化计算，获得满足脱靶量、碰撞角和末端攻角要求的参数k₁和k₂的范围；找到末端速度最大和末端速度最小所对应的k₁和k₂值，根据给定的末端速度，确定参数k₁和k₂。本发明方法基于逆动力学，同时考虑了末端碰撞角、攻角和速度约束，能够同时保证碰撞角、末端攻角和末端碰撞速度的要求，可实现对固定目标的精确打击。

Description

基于Bezier曲线的多约束飞行器导引方法

技术领域

本发明属于飞行器动力学与制导技术领域，涉及一种飞行器末段制导方法，具体涉及一种飞行器自主攻击地面固定目标时使用的一种导引方法。

背景技术

导弹制导规律即导引律是空战中实现战机追踪/拦截导引的火控***关键技术之一。导引律的选择对导弹始发欧精确打击目标至关重要。经典的导引律主要考虑的因素是脱靶量。随着技术的发展，无论反舰导弹还是地对地导弹，在保证脱靶量的同时，按照期望的碰撞角交会目标，并且交会时能保持期望的姿态和速度是设计制导律的目标。导引律主要有基于解析形式的显式导引律和基于数值优化的导引律。

参考文献[1]～[10]用不同的表达形式实现了不同约束的显式导引律。Lu设计了一种基于比例导引的自适应闭环导引律，可以按照期望的方向高精度的命中目标(参考文献[1]:P.Lu,D.B.Doman,and J.D.Schierman,“Adaptive Terminal Guidance for Hypervelocity Impact inSpecified Direction,”Journal of Guidance Control and Dynamics,Vol.29,No.2,2006,pp.269-278)；A.Ratnoo提出了一种基于比例导引的制导律，可以实现在平面内按照所有可能的碰撞角命中固定目标(参考文献[2]：A.Ratnoo,and D.Ghose,“Impact Angle Constrained Interception ofStationary Targets,”Journal of Guidance Control and Dynamics,Vol.31,No.6,2008,pp.1816-1821)；C.Ryoo设计了一种满足终端碰撞角的能量最优导引律(参考文献[3]：C.Ryoo,H.Cho,and M.Tahk,“Optimal Guidance Laws with Impact Angle Constraint,”Journal of GuidanceControl and Dynamics,Vol.28,No.4,2005,pp.724-732)；Y.Lee设计了一种最优导引律(参考文献[4]：Y.Lee,C.Ryoo,and E.Kim,“Optimal Guidance with Constraints on Impact Angle andTerminal Acceleration,”AIAA Guidance,Navigation,and Control Conference and Exhibit,2003,Austin,Texas,pp.1-7)，在满足脱靶量和终端碰撞角要求的同时可以调整终端加速度；R.York则设计了一种可同时满足碰撞角和终端攻角要求的最优导引律(参考文献[5]：R.J.York,and H.L.Pastrick,“Optimal Terminal Guidance with Constraints at Final Time,”Journal of Spacecraft,Vol.14,No.6,1977,pp.381-383)；I.R.Manchester设计了一种有碰撞角约束的平面交会导引律，这种导引律是采用了跟踪一个圆弧至目标的思路，不需要剩余航程的信息(参考文献[6]:I.R.Manchester,and A.V.Savkin,“Circular-Navigation-GuidanceLawfor Precision Missile/TargetEngagements,”Journal of Guidance Control and Dynamics,Vol.29,No.2,2006,pp.314-320)。B.S.Kim设计了一种偏置的比例导引律，此导引律是在传统的比例导引的基础上增加一时变的偏置项，可实现以期望的姿态命中目标(参考文献[7]:B.S.Kim,J.G.Lee,and H.S.Han,“Biased PNG Law for impact with Angular Constraint,”IEEE Transactions on Aerospace andElectronics Systems,Vol.34,No.1,1998,pp.277-288)。R.W.Morgan设计了一种法向加速度约束的的能量最优的导引律(参考文献[8]:R.W.Morgan,H.Tharp,and T.L.Vincent,“MinimumEnergy Guidance for Aerodynamically Controlled Missiles,”IEEE Transactions on AutomaticControl,Vol.56,No.9,2011,pp.2026-2037)。D.Sang设计了一种考虑导引头视场角限制的导引律，通过转换逻辑可保证目标始终在导引头锁定视场内(参考文献[9]:D.Sang,C.Ryoo,and M.Tahk,“A Guidance Law with a Switching Logic for Maintaining Seeker’s Lock-on for StationaryTargets,”KSAS International Journal,Vol.9,No.2,2008,pp.87-97)。A.Naghash基于逆动力学方法设计了一种次优显式导引律。通过离线优化描述弹道的三阶Bezier曲线的相关参数，获得末速最大的弹道，然后将此参数用于导引律，但是此方法没有考虑碰撞角约束以及末端攻角等约束(参考文献[10]:A.Naghash,R.Esmaelzadeh,M.Mortazavi,and R.Jamilnia,“Near OptimalGuidance Law for Decent to a Point Using Inverse Problem Approach,”Aerospace Science andTechnology,No.12,2008,pp.241-247)。

参考文献[11]～[13]基于在线数值优化的导引律实现对固定目标以期望的角度精确碰撞。K.P.Bollino设计了一种基于伪谱方法在线优化的导引律，通过在线优化来不断修正扰动和不确定性带来的偏差，此方法对弹载计算机的计算能力有较高的要求(参考文献[11]:K.P.Bollino,I.M.Ross,and D.D.Doman,“Optimal Nonlinear Feedback Guidance for Reentry Vehicles,”AIAAGuidance,Navigation,and Control Conference and Exhibit,2006,Keystone,Colorado,pp.1-20)。H.B.Oza成功应用非线性模型预测静态规划方法于空对地导弹的制导上，此方法在保证脱靶量的同时，可以保证碰撞角的约束，以法向过载最小作为规划的目标(参考文献[12]:H.B.Oza,and R.Padhi,“Impact-Angle-Constrained Suboptimal Model Predictive StaticProgrammingGuidance of Air-to-Ground Missiles,”Journal of Guidance Control and Dynamics,Vol.35,No.1,2012,pp.153-164)。A.Ratnoo将二维平面内碰撞角约束的导引问题描述成非线性规划问题，然后利用SDRE技术进行了求解(参考文献[13]:A.Ratnoo,and D.Ghose,“State-DependentRiccati-Equation-Based Guidance Law for Impact-Angle-Constrained Trajectories,”Journal ofGuidance Control and Dynamics,Vol.32,No.1,2009,pp.320-325)。

但上述现有技术中都还没有记载能够很好地同时保证碰撞角、末端攻角以及调节末端碰撞速度的导引方法。

发明内容

本发明的目的是为了解决上述问题，提出一种基于逆动力学方法的考虑末端碰撞角、攻角和速度约束的导引方法，即基于Bezier曲线的多约束飞行器导引方法，本发明方法能够同时保证碰撞角、末端攻角以及调节末端碰撞速度。

本发明提出的基于Bezier曲线的多约束飞行器导引方法，包括以下几个步骤：

第一步、根据具体的飞行器运动模型和三阶Bezier曲线方程，确定导引控制指令，导引控制指令采用飞行器的攻角α和倾侧角σ作为控制变量；

第二步、通过优化计算，获得满足脱靶量、碰撞角和末端攻角要求的k₁和k₂的范围；

第三步、在第二步得到的范围内，找到末端速度最大和末端速度最小所对应的k₁和k₂值，根据给定的末端速度，确定参数k₁和k₂。

所述的第一步确定飞行器的攻角α和倾侧角σ的方法是：

首先，采用三阶Bezier曲线来拟合铅垂面和水平面的弹道，具体公式如下：

$\left\{\begin{array}{c}x={(1-\mathrm{\τ})}^3{x}_0+3\mathrm{\τ}{(1-\mathrm{\τ})}^2{x}_1+3{\mathrm{\τ}}^2(1-\mathrm{\τ})x_2+{\mathrm{\τ}}^3{x}_f\\ y={(1-\mathrm{\τ})}^3{y}_0+3\mathrm{\τ}{(1-\mathrm{\τ})}^2{y}_1+3{\mathrm{\τ}}^2(1-\mathrm{\τ})y_2+{\mathrm{\τ}}^3{y}_f\\ z={(1-\mathrm{\τ})}^3{z}_0+3\mathrm{\τ}{(1-\mathrm{\τ})}^2{z}_1+3{\mathrm{\τ}}^2(1-\mathrm{\τ})z_2+{\mathrm{\τ}}^3{z}_f\end{array}\right.$

其中，x、y、z是拟合的飞行器的空间弹道上的任一位置坐标；参数τ∈[0,1]，τ＝0代表起点，τ＝1代表终点；(x₀,y₀,z₀)、(x₁,y₁,z₁)、(x₂,y₂,z₂)、(x_f,y_f,z_f)为控制弹道形状的四个控制点，(x₀,y₀,z₀)和(x_f,y_f,z_f)分别为弹道的起点和终点；四个控制点(x₀,y₀,z₀)、(x₁,y₁,z₁)、(x₂,y₂,z₂)和(x_f,y_f,z_f)映射在XOY平面上的点分别为(x₀,y₀)、(x₁,y₁)、(x₂,y₂)和(x_f,y_f)，则点(x₀,y₀)与(x₁,y₁)的连线与点(x₂,y₂)和(x_f,y_f)的连线的交点为(x_m,y_m)；设参数参数则中间两个控制点(x₁,y₁,z₁)和(x₂,y₂,z₂)用参数k₁和k₂来描述并获取：

$\left\{\begin{array}{c}{x}_1=x_0+k_1(x_m-x_0)\\ y_1=y_0+\frac{\tan {\mathrm{\γ}}_0}{\cos {\mathrm{\ψ}}_0}(x_1-x_0)\\ z_1=z_0-\tan {\mathrm{\ψ}}_0(x_1-x_0)\end{array}\right.;$ $\left\{\begin{array}{c}{x}_2=x_m+k_2(x_f-x_m)\\ y_2=y_f-\frac{\tan {\mathrm{\γ}}_f}{\cos {\mathrm{\ψ}}_f}(x_f-x_2)\\ z_2=z_f+\tan {\mathrm{\ψ}}_f(x_f-x_2)\end{array}\right.;k_1\∈\left[\mathrm{0,1}\right],k_2\∈\left[\mathrm{0,1}\right]$

进一步，确定弹道坐标系下的法向加速度的指令a_yB和侧向加速度的指令a_zB：

a_yB＝gcosγ+v²cosγcosψdγ/dx

a_zB＝-v²cos²γcosψdψ/dx

其中，g是重力加速度，γ是弹道倾角，v是飞行器的速度，ψ是弹道偏角。

其次，确定弹道坐标系下的法向加速度a_y和侧向加速度a_z：

$\begin{array}{cc}\left\{\begin{array}{c}{a}_y=a_{\mathrm{yB}}\left(0\right)\·n_{\⊥\max }\·g/\sqrt{a_{\mathrm{yB}}{\left(0\right)}^2+a_{\mathrm{zB}}{\left(0\right)}^2}\\ a_z=a_{\mathrm{zB}}\left(0\right)\·n_{\⊥\max }\·g/\sqrt{a_{\mathrm{yB}}{\left(0\right)}^2+a_{\mathrm{zB}}{\left(0\right)}^2}\end{array}\right.& \mathrm{if}\sqrt{a_{\mathrm{yB}}{\left(0\right)}^2+a_{\mathrm{zB}}{\left(0\right)}^2}/g>n_{\⊥\max }\\ \left\{\begin{array}{c}{a}_y=a_{\mathrm{yB}}\left(0\right)\\ a_z=a_{\mathrm{zB}}\left(0\right)\end{array}\right.& \mathrm{if}\sqrt{a_{\mathrm{yB}}{\left(0\right)}^2+a_{\mathrm{zB}}{\left(0\right)}^2}/g\≤n_{\⊥\max }\end{array}$

其中，a_yB(0)是在初始时刻所对应的法向加速度的指令，a_zB(0)是在初始时刻所对应的侧向加速度的指令，n_⊥max是飞行器的可用法向过载。

然后，获得倾侧角指令σ_c和攻角指令α_c： $\left\{\begin{array}{c}{\mathrm{\σ}}_c=\arctan (a_z/a_y)\\ {\mathrm{\α}}_c={f_L}^{-1}(M_a,{\mathrm{ma}}_y/\left(\mathrm{qS}\cos \mathrm{\σ}\right))\end{array}\right.;$ 其中，M_a为马赫数，M_a＝v/a，a为当地音速，m表示飞行器的质量，动压q＝ρv²/2，ρ是大气密度，S是参考面积，f_L是升力系数C_L的表示函数。

最后，确定飞行器的倾侧角σ：

$\mathrm{\&sigma;}=\left\{\begin{array}{cc}0& \mathrm{if}{a}_y=0\\ \left.\begin{array}{c}\left.\begin{array}{cc}{\mathrm{\&sigma;}}_c& \mathrm{if}{\mathrm{\&sigma;}}_{\min }\&le;{\mathrm{\&sigma;}}_c\&le;{\mathrm{\&sigma;}}_{\max }\\ {\mathrm{\&sigma;}}_{\max }& \mathrm{if}{\mathrm{\&sigma;}}_c>{\mathrm{\&sigma;}}_{\max }\\ {\mathrm{\&sigma;}}_{\min }& \mathrm{if}{\mathrm{\&sigma;}}_c<{\mathrm{\&sigma;}}_{\min }\end{array}\right\}\mathrm{if}\left|\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}\right|\&le;{\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }\\ \left.\begin{array}{cc}{\mathrm{\&sigma;}}_{\mathrm{prev}}+\left({t-t}_{\mathrm{prev}}\right){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }& \mathrm{if}{\mathrm{\&sigma;}}_{\min }\&le;{\mathrm{\&sigma;}}_{\mathrm{prev}}+\left({t-t}_{\mathrm{prev}}\right){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }\&le;{\mathrm{\&sigma;}}_{\max }\\ {\mathrm{\&sigma;}}_{\max }& \mathrm{if}{\mathrm{\&sigma;}}_{\mathrm{prev}}+(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }>{\mathrm{\&sigma;}}_{\max }\\ {\mathrm{\&sigma;}}_{\min }& \mathrm{if}{\mathrm{\&sigma;}}_{\mathrm{prev}}+\left({t-t}_{\mathrm{prev}}\right){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }<{\mathrm{\&sigma;}}_{\min }\end{array}\right\}\\ \left.\begin{array}{cc}{\mathrm{\&sigma;}}_{\mathrm{prev}}-\left({t-}_{\mathrm{prev}}\right){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }& \mathrm{if}{\mathrm{\&sigma;}}_{\min }\&le;{\mathrm{\&sigma;}}_{\mathrm{prev}}-(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }\&le;{\mathrm{\&sigma;}}_{\max }\\ {\mathrm{\&sigma;}}_{\max }& \mathrm{if}{\mathrm{\&sigma;}}_{\mathrm{prev}}-\left({t-t}_{\mathrm{prev}}\right){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }>{\mathrm{\&sigma;}}_{\max }\\ {\mathrm{\&sigma;}}_{\min }& \mathrm{if}{\mathrm{\&sigma;}}_{\mathrm{prev}}-\left({t-t}_{\mathrm{prev}}\right){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }<{\mathrm{\&sigma;}}_{\min }\end{array}\right\}\mathrm{if}\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}<-{\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }\end{array}i\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}>{\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }>\right\}& \mathrm{if}{a}_y\&NotEqual;0\end{array}\right.$

其中，是设定的常数，表示倾侧角对时间求导的最大数值；σ_min和σ_max是设定的常数，分别表示侧倾角需设定的最小值和最大值；σ_prev表示前一制导周期的倾侧角，t_prev表示前一制导周期的时间，t表示当前时间。

根据得到的倾侧角σ，能够确定攻角指令α_c，进一步确定飞行器的攻角α：

$\mathrm{\&alpha;}=\left\{\begin{array}{cc}\left.\begin{array}{cc}{\mathrm{\&alpha;}}_c& \mathrm{if}{\mathrm{\&alpha;}}_{\min }\&le;{\mathrm{\&alpha;}}_c\&le;{\mathrm{\&alpha;}}_{\max }\\ {\mathrm{\&alpha;}}_{\max }& \mathrm{if}{\mathrm{\&alpha;}}_c>{\mathrm{\&alpha;}}_{\max }\\ {\mathrm{\&alpha;}}_{\min }& \mathrm{if}{\mathrm{\&alpha;}}_c<{\mathrm{\&alpha;}}_{\min }\end{array}\right\}& \mathrm{if}\left|\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}\right|\&le;{\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }\\ \left.\begin{array}{cc}{\mathrm{\&alpha;}}_{\mathrm{prev}}+(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }& \mathrm{if}{\mathrm{\&alpha;}}_{\min }\&le;{\mathrm{\&alpha;}}_{\mathrm{prev}}+(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }\&le;{\mathrm{\&alpha;}}_{\max }\\ {\mathrm{\&alpha;}}_{\max }& \mathrm{if}{\mathrm{\&alpha;}}_{\mathrm{prev}}+(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }>{\mathrm{\&alpha;}}_{\max }\\ {\mathrm{\&alpha;}}_{\min }& \mathrm{if}{\mathrm{\&alpha;}}_{\mathrm{prev}}+(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }<{\mathrm{\&alpha;}}_{\min }\end{array}\right\}& \mathrm{if}\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}>{\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }\\ \left.\begin{array}{cc}{\mathrm{\&alpha;}}_{\mathrm{prev}}-(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }& \mathrm{if}{\mathrm{\&alpha;}}_{\min }\&le;{\mathrm{\&alpha;}}_{\mathrm{prev}}-(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }\&le;{\mathrm{\&alpha;}}_{\max }\\ {\mathrm{\&alpha;}}_{\max }& \mathrm{if}{\mathrm{\&alpha;}}_{\mathrm{prev}}-(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }>{\mathrm{\&alpha;}}_{\max }\\ {\mathrm{\&alpha;}}_{\min }& \mathrm{if}{\mathrm{\&alpha;}}_{\mathrm{prev}}-(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }\&le;{\mathrm{\&alpha;}}_{\min }\end{array}\right\}& \mathrm{if}\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}<-{\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }\end{array}\right.$

其中，是设定的常数，表示攻角对时间求导的最大数值；α_min和α_max是设定的常数，分别表示攻角需设定的最小值和最大值；α_prev表示前一制导周期的攻角。

所述的第三步中，在标称条件下，要实现期望的末速，首先固定参数k₁的值，然后采用Secant法根据下式来实时修正k₂的值：

$k_2(i+1)=k_2\left(i\right)+\frac{(v_{\mathrm{fdes}}-v_f\left(i\right))(k_2\left(i\right)-k_2(i-1))}{(v_f\left(i\right)-v_f(i-1))}$

k₂(i)表示第i次在线迭代过程中的k₂值，v_fdes表示期望的末速值，v_f(i)表示第i次仿真预测的末速值；

经过在线迭代，确定参数k₂，并收敛至期望的末速。

本发明的飞行器导引方法的优点和积极效果在于：

(1)可实现对固定目标的精确打击；

(2)满足末端碰撞角的要求；

(3)满足末端小攻角的要求；

(4)满足攻角、法向过载等过程量的边界要求。

附图说明

图1是本发明提供的飞行器导引方法的整体步骤流程示意图；

图2是本发明步骤一中用Bezier曲线估计的弹道的示意图；

图3是不同末端弹道倾角情况下的弹道倾角、速度、攻角变化曲线；

图4是不同末端弹道倾角情况下的铅垂面弹道曲线；

图5是满足脱靶量以及末端弹道倾角要求的k₁和k₂取值范围的示意图；

图6是末端速度随k₁和k₂变化的曲面；

图7是末端速度随k₁和k₂变化的曲线；

图8是末端速度、末端弹道倾角和脱靶量随k₂变化的曲线。

具体实施方式

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

本发明针对无动力再入滑翔飞行器末段制导问题，提出的一种基于Bezier曲线的多约束飞行器导引方法，流程如图1所示，包括以下几个步骤：

第一步、根据飞行器运动模型和三阶Bezier曲线的方程，推导出导引控制指令。

步骤1.1：建立再入滑翔飞行器运动模型，该模型包括运动方程和飞行器运动过程中所受的约束条件。

由于高超声速飞行器在下降段飞行距离较小，可假设地球为一平面大地，飞行器在下降段被看作无动力的质点，其下降的运动方程如下所示：

$\begin{array}{c}\stackrel{\&CenterDot;}{v}=(D-\mathrm{mg}\sin \mathrm{\&gamma;})/m\\ \stackrel{\&CenterDot;}{\mathrm{\&gamma;}}=(L\cos \mathrm{\&sigma;}-\mathrm{mg}\cos \mathrm{\&gamma;})/\left(\mathrm{mv}\right)\\ \stackrel{\&CenterDot;}{\mathrm{\&psi;}}=L\sin \mathrm{\&sigma;}/(-\mathrm{mv}\cos \mathrm{\&gamma;})\\ \stackrel{\&CenterDot;}{x}=v\cos \mathrm{\&gamma;}\cos \mathrm{\&psi;}\\ \stackrel{\&CenterDot;}{y}=v\sin \mathrm{\&gamma;}\\ \stackrel{\&CenterDot;}{z}=-v\cos \mathrm{\&gamma;}\sin \mathrm{\&psi;}\end{array}---\left(1\right)$

其中升力L和阻力D的定义如下：

L＝qSC_L D＝qSC_D   (2)

其中，动压q＝ρv²/2。m表示飞行器的质量，g表示重力加速度，γ是弹道倾角，σ是飞行器的倾侧角，v是飞行器的速度，ψ是弹道偏角，x、y、z是飞行器的位置坐标。字符上加点表示导数。所述的弹道是指飞行器的飞行轨迹。大气密度ρ的计算采用1976年美国标准大气COESA模型。S是参考面积，C_L是升力系数，C_D是阻力系数，C_L和C_D由给定的数据拟合得到，可描述成如下形式：

C_L＝f_L(M_a,α) C_D＝f_D(M_a,α)   (3)

其中，f_L是升力系数C_L的表示函数，f_D是阻力系数C_D的表示函数，α是飞行器的攻角，M_a为马赫数，M_a＝v/a，a为当地音速。

飞行器在下降的运动过程中所受的约束条件有：落角约束、落速约束、控制约束以及法向过载约束。

(1)由于战斗部的启爆要求以及具体任务的需要，末端弹道倾角需要达到指定的角度，则弹道倾角γ的落角约束条件为：

γ_f min≤γ_f≤γ_f max   (4)

其中，γ_f min和γ_f max是一指定的常数，一般设置在-90度至0度之间，根据实际需要选定，设t_f表示导弹飞行的终端时刻，γ_f表示终端时刻的弹道倾角的值。

(2)出于任务要求及制导***的需求，对飞行器的落速也提出了严格的要求，落速约束为：

v_f min≤v_f≤v_f max   (5)

其中v_f min和v_f max是一指定的常数，v_f max大于v_f min，一般都是大于0的值，具体根据实际需要选定，例如可设置两个常数为1000m/s左右。v_f表示终端时刻t_f时飞行器的速度。

(3)本发明中采用攻角α和倾侧角σ作为控制变量，根据飞行器的总体要求，需要攻角和倾侧角保证在一定的范围内：

α_min≤α≤α_max σ_min＜σ＜σ_max   (6)

其中，α_min、α_max、σ_min、σ_max是给定的常数。α_min和α_max分别表示攻角需设定的最小值和最大值。σ_min和σ_max分别表示侧倾角需设定的最小值和最大值。

由于考虑到动态过程和可实现性，对攻角和倾侧角的变化率也有着相应的约束:

$\left|\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}\right|\&le;{\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }\left|\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}\right|\&le;{\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }---\left(7\right)$

其中，是给定的常数，单位为：deg/s，一般设置为0.1到90之间，可根据实际需要设置表示倾侧角对时间求导的最大数值。表示攻角对时间求导的最大数值。

(4)在下降过程中飞行器的法向过载n_⊥的绝对值不能大于飞行器的可用法向过载n_⊥max，此约束描述如下：

|n_⊥|≤n_⊥max   (8)

在下降过程中飞行器的热流、动压也需要符合一定的约束，但是在研究中，发现热流和动压均能符合要求，不是此段的主要矛盾，故不单独考虑。

步骤1.2：获取基于三阶Bezier曲线的多约束导引指令。

将飞行器下降段弹道投影至地面坐标系的Ox_dy_d的铅垂平面和Oy_dz_d的水平平面。弹道在这两个平面按照给定的Bezier曲线规律运动。

n阶的Bezier曲线P(τ)定义为：

$P\left(\mathrm{\&tau;}\right)=\underset{i=0}{\overset{n}{\mathrm{\&Sigma;}}}{B}_i{J}_{n,i}\left(\mathrm{\&tau;}\right)---\left(9\right)$

其中，B_i表示第i个控制点，Bernstein基函数参数τ∈[0,1]，τ＝0代表起点，τ＝1代表终点。

多项式(9)的r阶导数可写成

$\frac{d^r}{d{\mathrm{\&tau;}}^r}P\left(\mathrm{\&tau;}\right)=\frac{n!}{(n-r)!}\underset{i=0}{\overset{n-r}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;}}^r{B}_i{J}_{n-r,i}---\left(10\right)$

其中，Δ⁰B_i＝B_i；Δ^kB_i＝Δ^k-1B_i+1-Δ^k-1B_i，k＝0,...,r。Δ^kB_i表示曲线P(τ)在第k阶求导时的第i个控制点。

本发明方法中采用三阶Bezier曲线来拟合铅垂面和水平面的弹道，具体公式如下：

$\left\{\begin{array}{c}x={(1-\mathrm{\&tau;})}^3{x}_0+3\mathrm{\&tau;}{(1-\mathrm{\&tau;})}^2{x}_1+3{\mathrm{\&tau;}}^2(1-\mathrm{\&tau;})x_2+{\mathrm{\&tau;}}^3{x}_f\\ y={(1-\mathrm{\&tau;})}^3{y}_0+3\mathrm{\&tau;}{(1-\mathrm{\&tau;})}^2{y}_1+3{\mathrm{\&tau;}}^2(1-\mathrm{\&tau;})y_2+{\mathrm{\&tau;}}^3{y}_f\\ z={(1-\mathrm{\&tau;})}^3{z}_0+3\mathrm{\&tau;}{(1-\mathrm{\&tau;})}^2{z}_1+3{\mathrm{\&tau;}}^2(1-\mathrm{\&tau;})z_2+{\mathrm{\&tau;}}^3{z}_f\end{array}\right.---\left(11\right)$

(x,y,z)是拟合的飞行器的空间弹道上的任位置坐标，(x₀,y₀,z₀)、(x₁,y₁,z₁)、(x₂,y₂,z₂)、(x_f,y_f,z_f)为控制弹道形状的四个控制点，其中，(x₀,y₀,z₀)、(x_f,y_f,z_f)分别为弹道的起点和终点，中间两个控制点(x₁,y₁,z₁)和(x₂,y₂,z₂)不在弹道上。中间两个控制点用参数k₁和k₂来描述并获取。如图2所示，四个控制点(x₀,y₀,z₀)、(x₁,y₁,z₁)、(x₂,y₂,z₂)和(x_f,y_f,z_f)映射在XOY平面上的点分别为(x₀,y₀)、(x₁,y₁)、(x₂,y₂)和(x_f,y_f)，则点(x₀,y₀)与(x₁,y₁)的连线与点(x₂,y₂)和(x_f,y_f)的连线的交点为(x_m,y_m)，根据平面几何关系，能够得到起点处的弹道倾角γ₀和终点处的弹道倾角γ_f具有如下关系：

$\tan {\mathrm{\&gamma;}}_0=\frac{y_m-y_0}{x_m-x_0};$ $\tan {\mathrm{\&gamma;}}_f=\frac{y_f-y_m}{x_f-x_m};$

进一步，可确定坐标(x_m,y_m)： $\left\{\begin{array}{c}{x}_m=\frac{x_f\tan {\mathrm{\&gamma;}}_f-x_0\tan {\mathrm{\&gamma;}}_0+y_0-y_f}{\tan {\mathrm{\&gamma;}}_f-\tan {\mathrm{\&gamma;}}_0}\\ y_m=-(x_f-x_m)\tan {\mathrm{\&gamma;}}_f+y_f\end{array}\right.$

设参数然后得到两个控制点(x₁,y₁,z₁)和(x₂,y₂,z₂)如下式：

$\left\{\begin{array}{c}{x}_1=x_0+k_1(x_m-x_0)\\ y_1=y_0+\frac{\tan {\mathrm{\&gamma;}}_0}{\cos {\mathrm{\&psi;}}_0}(x_1-x_0)\\ z_1=z_0-\tan {\mathrm{\&psi;}}_0(x_1-x_0)\end{array}\right.---\left(12\right)$

$\left\{\begin{array}{c}{x}_2=x_m+k_2(x_f-x_m)\\ y_2=y_f-\frac{\tan {\mathrm{\&gamma;}}_f}{\cos {\mathrm{\&psi;}}_f}(x_f-x_2)\\ z_2=z_f+\tan {\mathrm{\&psi;}}_f(x_f-x_2)\end{array}\right.---\left(13\right)$

一般情况下，k₁∈[0,1]，k₂∈[0,1]；只要选择合适的k₁和k₂，根据公式(12)和(13)，就可以获得相应的中间两个控制点的坐标，这样就可以控制弹道的形状了。用Bezier曲线估计的弹道如图2所示。

将运动方程(1)写成对x求导的形式，具体方程如下：

dv/dx＝(D-mgsinγ)/(mvcosγcosψ)

dγ/dx＝(Lcosσ-mgcosγ)/(mv²cosγcosψ)

dψ/dx＝Lsinσ/(-mv²cos²γcosψ)   (14)

dy/dx＝tanγ/cosψ

dz/dx＝-tanψ

弹道坐标系下的法向加速度a_y和侧向加速度a_z可定义为：

a_y＝Lcosσ/m；a_z＝Lsinσ/m   (15)

根据动力学方程(14)可求解出法向加速度的指令a_yB和侧向加速度的指令a_zB为：

a_yB＝gcosγ+v²cosγcosψdγ/dx

                        (16)

a_zB＝-v²cos²γcosψdψ/dx

为了进一步求解式(16)中的dγ/dx和dψ/dx，由动力学方程可得，

tanγ＝cosψdy/dx；tanψ＝-dz/dx   (17)

对上面(17)中的两式求导，可得

dγ/dx＝(d²y/dx²cosψ-dy/dx·dψ/dx·sinψ)cos²γ   (18)

dψ/dx＝-cos²ψ·d²z/dx²

而根据三阶Bezier曲线公式(11)，可得其对τ的一阶导数和二阶导数如式(19)和(20)所示。

$\left\{\begin{array}{c}\mathrm{dx}/\mathrm{d\&tau;}=-3{(1-\mathrm{\&tau;})}^2{x}_0+3x_1(3{\mathrm{\&tau;}}^2-4\mathrm{\&tau;}+1)+3x_2(-3{\mathrm{\&tau;}}^2+2\mathrm{\&tau;}))+3{\mathrm{\&tau;}}^2{x}_f\\ \mathrm{dy}/\mathrm{d\&tau;}=-3{(1-\mathrm{\&tau;})}^2{y}_0+3y_1(3{\mathrm{\&tau;}}^2-4\mathrm{\&tau;}+1)+3y_2(-3{\mathrm{\&tau;}}^2+2\mathrm{\&tau;}))+3{\mathrm{\&tau;}}^2{y}_f\\ \mathrm{dz}/\mathrm{d\&tau;}=-3{(1-\mathrm{\&tau;})}^2{z}_0+3z_1(3{\mathrm{\&tau;}}^2-4\mathrm{\&tau;}+1)+3z_2(-3{\mathrm{\&tau;}}^2+2\mathrm{\&tau;}))+3{\mathrm{\&tau;}}^2{z}_f\end{array}\right.---\left(19\right)$

$\left\{\begin{array}{c}{d}^{2}x/d{\mathrm{\&tau;}}^2=6(1-\mathrm{\&tau;})x_0+6x_1(3\mathrm{\&tau;}-2)+6x_2(-3\mathrm{\&tau;}+1)+6\mathrm{\&tau;}{x}_f\\ d^{2}y/d{\mathrm{\&tau;}}^2=6(1-\mathrm{\&tau;})y_0+6y_1(3\mathrm{\&tau;}-2)+6y_2(-3\mathrm{\&tau;}+1)+6\mathrm{\&tau;}{y}_f\\ d^{2}z/d{\mathrm{\&tau;}}^2=6(1-\mathrm{\&tau;})z_0+6z_1(3\mathrm{\&tau;}-2)+6z_2(-3\mathrm{\&tau;}-1)+6\mathrm{\&tau;}{z}_f\end{array}\right.---\left(20\right)$

故可求出y和z对x的一阶导数和二阶导数，分别如式(21)和(22)所示。

$\frac{\mathrm{dy}}{\mathrm{dx}}=\frac{\mathrm{dy}/\mathrm{d\&tau;}}{\mathrm{dx}/\mathrm{d\&tau;}};\frac{\mathrm{dz}}{\mathrm{dx}}=\frac{\mathrm{dz}/\mathrm{d\&tau;}}{\mathrm{dx}/\mathrm{d\&tau;}}---\left(21\right)$

$\frac{d^{2}y}{d{x}^2}=\frac{\frac{d^{2}y}{d{\mathrm{\&tau;}}^2}-\frac{d^{2}x}{d{t}^2}\frac{\mathrm{dy}/\mathrm{dt}}{\mathrm{dx}/\mathrm{dt}}}{{\left(\frac{\mathrm{dx}}{\mathrm{dt}}\right)}^2};\frac{d^{2}z}{d{x}^2}=\frac{\frac{d^{2}z}{d{\mathrm{\&tau;}}^2}-\frac{d^{2}x}{d{t}^2}\frac{\mathrm{dz}/\mathrm{dt}}{\mathrm{dx}/\mathrm{dt}}}{{\left(\frac{\mathrm{dx}}{\mathrm{dt}}\right)}^2}---\left(22\right)$

故公式(16)的法向加速度指令得以求解，但是此指令是根据Bezier曲线弹道获得的。在制导回路中，需要实时根据当前的状态和给定的末端状态来求得法向加速度a_yB和侧向加速度a_zB，但实际使用的只是每次求得的法向加速度和侧向加速度在初始时刻的值，即τ＝0时所对应的值a_yB(0)和a_zB(0)。每个制导周期都会根据不同的当前状态来更新Bezier曲线，从而获得新的制导指令。

由于法向过载是有界的，故考虑过载约束的法向加速度如下：

$\begin{array}{cc}\left\{\begin{array}{c}{a}_y=a_{\mathrm{yB}}\left(0\right)\&CenterDot;n_{\&perp;\max }\&CenterDot;g/\sqrt{a_{\mathrm{yB}}{\left(0\right)}^2+a_{\mathrm{zB}}{\left(0\right)}^2}\\ a_z=a_{\mathrm{zB}}\left(0\right)\&CenterDot;n_{\&perp;\max }\&CenterDot;g/\sqrt{a_{\mathrm{yB}}{\left(0\right)}^2+a_{\mathrm{zB}}{\left(0\right)}^2}\end{array}\right.& \mathrm{if}\sqrt{a_{\mathrm{yB}}{\left(0\right)}^2+a_{\mathrm{zB}}{\left(0\right)}^2}/g>n_{\&perp;\max }\\ \left\{\begin{array}{c}{a}_y=a_{\mathrm{yB}}\left(0\right)\\ a_z=a_{\mathrm{zB}}\left(0\right)\end{array}\right.& \mathrm{if}\sqrt{a_{\mathrm{yB}}{\left(0\right)}^2+a_{\mathrm{zB}}{\left(0\right)}^2}/g\&le;n_{\&perp;\max }\end{array}---\left(23\right)$

根据升力及升力系数的定义(2)和(3)，可得

σ_c＝arctan(a_z/a_y)

                              (24)

α_c＝f_L ^-1(M_a, ma_y/(qScosσ))

σ_c表示倾侧角指令，α_c表示攻角指令。

但由于攻角和倾侧角的大小以及其变化率均需满足约束条件，故实际使用的倾侧角σ和攻角α由式(25)和(26)确定。

$\mathrm{\&sigma;}=\left\{\begin{array}{cc}0& \mathrm{if}{a}_y=0\\ \left.\begin{array}{c}\left.\begin{array}{cc}{\mathrm{\&sigma;}}_c& \mathrm{if}{\mathrm{\&sigma;}}_{\min }\&le;{\mathrm{\&sigma;}}_c\&le;{\mathrm{\&sigma;}}_{\max }\\ {\mathrm{\&sigma;}}_{\max }& \mathrm{if}{\mathrm{\&sigma;}}_c>{\mathrm{\&sigma;}}_{\max }\\ {\mathrm{\&sigma;}}_{\min }& \mathrm{if}{\mathrm{\&sigma;}}_c<{\mathrm{\&sigma;}}_{\min }\end{array}\right\}\mathrm{if}\left|\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}\right|\&le;{\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }\\ \left.\begin{array}{cc}{\mathrm{\&sigma;}}_{\mathrm{prev}}+\left({t-t}_{\mathrm{prev}}\right){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }& \mathrm{if}{\mathrm{\&sigma;}}_{\min }\&le;{\mathrm{\&sigma;}}_{\mathrm{prev}}+\left({t-t}_{\mathrm{prev}}\right){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }\&le;{\mathrm{\&sigma;}}_{\max }\\ {\mathrm{\&sigma;}}_{\max }& \mathrm{if}{\mathrm{\&sigma;}}_{\mathrm{prev}}+(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }>{\mathrm{\&sigma;}}_{\max }\\ {\mathrm{\&sigma;}}_{\min }& \mathrm{if}{\mathrm{\&sigma;}}_{\mathrm{prev}}+\left({t-t}_{\mathrm{prev}}\right){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }<{\mathrm{\&sigma;}}_{\min }\end{array}\right\}\\ \left.\begin{array}{cc}{\mathrm{\&sigma;}}_{\mathrm{prev}}-\left({t-}_{\mathrm{prev}}\right){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }& \mathrm{if}{\mathrm{\&sigma;}}_{\min }\&le;{\mathrm{\&sigma;}}_{\mathrm{prev}}-(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }\&le;{\mathrm{\&sigma;}}_{\max }\\ {\mathrm{\&sigma;}}_{\max }& \mathrm{if}{\mathrm{\&sigma;}}_{\mathrm{prev}}-\left({t-t}_{\mathrm{prev}}\right){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }>{\mathrm{\&sigma;}}_{\max }\\ {\mathrm{\&sigma;}}_{\min }& \mathrm{if}{\mathrm{\&sigma;}}_{\mathrm{prev}}-\left({t-t}_{\mathrm{prev}}\right){\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }<{\mathrm{\&sigma;}}_{\min }\end{array}\right\}\mathrm{if}\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}<-{\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }\end{array}i\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}>{\stackrel{\&CenterDot;}{\mathrm{\&sigma;}}}_{\max }>\right\}& \mathrm{if}{a}_y\&NotEqual;0\end{array}\right.---\left(25\right)$

$\mathrm{\&alpha;}=\left\{\begin{array}{cc}\left.\begin{array}{cc}{\mathrm{\&alpha;}}_c& \mathrm{if}{\mathrm{\&alpha;}}_{\min }\&le;{\mathrm{\&alpha;}}_c\&le;{\mathrm{\&alpha;}}_{\max }\\ {\mathrm{\&alpha;}}_{\max }& \mathrm{if}{\mathrm{\&alpha;}}_c>{\mathrm{\&alpha;}}_{\max }\\ {\mathrm{\&alpha;}}_{\min }& \mathrm{if}{\mathrm{\&alpha;}}_c<{\mathrm{\&alpha;}}_{\min }\end{array}\right\}& \mathrm{if}\left|\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}\right|\&le;{\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }\\ \left.\begin{array}{cc}{\mathrm{\&alpha;}}_{\mathrm{prev}}+(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }& \mathrm{if}{\mathrm{\&alpha;}}_{\min }\&le;{\mathrm{\&alpha;}}_{\mathrm{prev}}+(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }\&le;{\mathrm{\&alpha;}}_{\max }\\ {\mathrm{\&alpha;}}_{\max }& \mathrm{if}{\mathrm{\&alpha;}}_{\mathrm{prev}}+(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }>{\mathrm{\&alpha;}}_{\max }\\ {\mathrm{\&alpha;}}_{\min }& \mathrm{if}{\mathrm{\&alpha;}}_{\mathrm{prev}}+(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }<{\mathrm{\&alpha;}}_{\min }\end{array}\right\}& \mathrm{if}\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}>{\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }\\ \left.\begin{array}{cc}{\mathrm{\&alpha;}}_{\mathrm{prev}}-(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }& \mathrm{if}{\mathrm{\&alpha;}}_{\min }\&le;{\mathrm{\&alpha;}}_{\mathrm{prev}}-(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }\&le;{\mathrm{\&alpha;}}_{\max }\\ {\mathrm{\&alpha;}}_{\max }& \mathrm{if}{\mathrm{\&alpha;}}_{\mathrm{prev}}-(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }>{\mathrm{\&alpha;}}_{\max }\\ {\mathrm{\&alpha;}}_{\min }& \mathrm{if}{\mathrm{\&alpha;}}_{\mathrm{prev}}-(t-t_{\mathrm{prev}}){\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }\&le;{\mathrm{\&alpha;}}_{\min }\end{array}\right\}& \mathrm{if}\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}<-{\stackrel{\&CenterDot;}{\mathrm{\&alpha;}}}_{\max }\end{array}\right.---\left(26\right)$

其中，σ_prev表示前一制导周期的倾侧角，t_prev表示前一制导周期的时间，α_prev表示前一制导周期的攻角，t表示当前时间。

所以只要选择好控制Bezier曲线形状的参数k₁和k₂，就可根据式(25)和(26)得到相应的控制，这样只需要得到合适的Bezier曲线的控制参数即可获得导引指令。

第二步、通过优化计算，获得满足脱靶量、碰撞角和末端攻角要求的k₁、k₂的范围；

公式(25)、(26)所描述的导引律是与参数k₁和k₂直接相关的，因为k₁、k₂的范围为[0,1]，所以可以通过枚举的方法循环计算来获得k₁和k₂所要求的范围，只要k₁和k₂位于此范围内，就可满足脱靶量、碰撞角和末端攻角的要求。

第三步、根据给定的末端速度，确定参数k₁和k₂。

从Bezier曲线的选取时就考虑了末端位置和碰撞角，所以只要飞行器具有足够的能力可以实现此Bezier曲线所描述的弹道，那么脱靶量和碰撞角均可满足。一般通过仿真，来获得k₁和k₂的范围，如图5所示的实施例。在特定的初始条件和末端条件下，可以获得保证脱靶量和末端碰撞角约束的k₁和k₂的范围。然后在此范围内，通过对k₁和k₂的优化，找到末端速度最大和末端速度最小所对应的k₁和k₂值。如果期望末端速度最大，只需选择末端速度最大所对应的k₁和k₂即可。

从仿真中可知末端速度对k₂更加敏感，对k₁的影响不大。故可固定k₁的值，靠k₂即可调节末速的大小。如图8中末速v_f与k₂的曲线图所示，末速随k₂单调递增。所以在标称条件下，要实现期望的末速，采用Secant法根据公式(27)来实时修正k₂的值，经过几次在线迭代，即可收敛至期望的值。

$k_2(i+1)=k_2\left(i\right)+\frac{(v_{\mathrm{fdes}}-v_f\left(i\right))(k_2\left(i\right)-k_2(i-1))}{(v_f\left(i\right)-v_f(i-1))}---\left(27\right)$

k₂(i)表示第i次在线迭代过程中的k₂值，v_fdes表示期望的末速值，v_f(i)表示第i次仿真预测的末速值。

实施例：

实施例中，以某高超声速升力体飞行器模型为例，此飞行器重907kg，空气动力系数与马赫数和攻角相关。实验1是验证在考虑法向过载约束、攻角和倾侧角约束、碰撞角约束以及终端攻角的前提下的制导律的适应性问题；实验2是进一步考虑终端速度约束的制导律问题。

实验1：不同碰撞角条件下的末速最大；

仿真初始条件和终端条件如表1所示。

表1 仿真初始条件和期望的末端条件

	x(km)	y(km)	z(km)	γ(°)	ψ(°)	α(°)	σ(°)	v(m/s)
									初始条件	0	35	0	0	0	0	0	2000
末端条件	200	0	0	-60，-75	0	[-2,2]		Max

采用本发明方法提供的制导方法，优点是在保证命中精度的同时能较好的保证末端碰撞角和末端攻角值很小。根据表1给出的条件，采用matlab 2008a中fmincon命令，优化参数k₁和k₂，使得末端速度最大，在不同末端弹道倾角条件下获得的参数值以及末端速度值如下表所示。

表2 两种末端弹道倾角条件下的仿真结果

	γf(°)	miss distance(m)	vf(m/s)	αf(°)	(k1,k2)
						情况1	-59.997	0.004	1129.6	-0.21	(0.85,0.35)
情况2	-74.995	0.02	1124.5	-0.26	(0.85,0.45)

miss distance表示脱靶量，就是终端时刻的位置偏差。根据以上条件得到的弹道如图3和图4所示。图3是两种末端弹道倾角下的弹道倾角γ、速度v、攻角α随时间变化的曲线示意图。图4是两种末端弹道倾角下的铅垂面弹道曲线示意图。图3和图4中，γ_f的值分别约写为-60deg和-75deg。

实验2：末端速度调节及k₁和k₂的选择问题；

在目标静止和选定末端碰撞角的情况下，可以通过调节k₁和k₂的值来调节末端速度的大小。在无扰动情况下，首先要找到保证脱靶量为0以及满足末端的弹道倾角约束的k₁和k₂的区间；然后在此区间的基础上求得对应指定末端速度的值。

仿真初始条件和终端条件如表3所示。

表3 仿真初始条件和期望的末端条件

	x(km)	y(km)	z(km)	γ(°)	ψ(°)	α(°)	σ(°)	v(m/s)
									初始条件	0	30	5	0	0	0	0	1850
末端条件	100	0	0	-89.99	0	[-2,2]	0

通过对k₁和k₂进行组合仿真，获得满足脱靶量约束和末端弹道倾角约束的k₁和k₂的范围，如图5所示。

进一步可在以上范围内，求得合适的k₁和k₂以满足末端速度需求。末速与k₁和k₂的关系如图6所示。

通过优化求得末端速度最小值点：(k₁,k₂)＝(0.80,0.1878)，对应的末速为：699.9m/s；末端速度最大值点：(k₁,k₂)＝(0.8241,0.5416)，对应的末速为：1196.2m/s。所以，要控制末端速度在[699.9,1196.2]之间，只需要找到合适的(k₁,k₂)点即可。从图7可知，当k₁取不同值时，对应的末端速度最大值和最小值变化都不大，也就是说末速主要与k₂相关。

本实例中选择k₁＝0.9；通过优化，可求得对应的末端速度最小值v_min的k₂值k₂_v_min＝0.2388，对应的末端速度最小值v_min为：699.94m/s，对应末端速度最大值v_max的k₂值k₂_v_max＝0.6297，对应的末端速度最大值v_max为：1195.3m/s。从优化结果可看出，获得的最大最小末端速度和同时优化k₁和k₂时获得的结果相近；并且从图7也可看出k₁在图5所示的范围内选择某一固定的值，只变化k₂就可获得接近最大和接近最小的末端速度。因此为了方便起见，在后面的分析中直接选择k₁＝0.9。

如果k₁＝0.9，k₂∈[k₂_v_min,k₂_v_max]，脱靶量和末端弹道倾角均可以很好地满足，如图8所示，并且末速随k₂在此区间单调递增，所以只要在此区间选择不同的k₂值就可得到不同的末端速度；k₂小于k₂_v_min或者大于k₂_v_max时，由于控制能力的限制，无法提供满足期望弹道的法向加速度，故无法保证脱靶量和末端碰撞角的要求；而由于落速v_f在[k₂_v_min,k₂_v_max]区间是单调递增的，所以要实现任意属于[699.94,1195.3]m/s之间的末端速度均可通过插值找到相应的参数k₂。

也就是说，在无扰动情况下，要想控制速度v_f∈[699.94,1195.3]m/s，根据图8找到合适的k₂即可。实际制导问题中，采用Secant法经过在线几次迭代即可找过合适的k₂，从而获得期望的末端速度。

仿真初始条件和终端条件如表3所示，并在此基础上要求末速为1100m/s；选择k₁＝0.9，采用迭代公式(27)，求k₂，选择k₂的初值为：0.5，迭代停止条件为|v_fdes-v_f|＜0.1，经过5次迭代即可收敛至k₂＝0.428；仿真结果如表4所示。

表4 仿真结果

	miss distance(m)	Δvf(m/s)	Δγf(°)	αf(°)	ayf
						仿真误差	0.004	0.05	2e-7	-0.35	0

其中，miss distance表示脱靶量，就是终端时刻的位置偏差，Δv_f表示实际终端速度与期望终端速度的偏差，Δγ_f表示实际弹道倾角与期望弹道倾角的偏差、α_f表示终端时刻的攻角，a_yf表示终端时刻的法向加速度。

Claims (2)

1.一种基于Bezier曲线的多约束飞行器导引方法，其特征在于，包括如下步骤：

第一步、根据飞行器运动模型和三阶Bezier曲线方程，确定导引控制指令，导引控制指令采用飞行器的攻角α和倾侧角σ作为控制变量，具体确定方法是：

首先，采用三阶Bezier曲线来拟合铅垂面和水平面的弹道，具体公式如下：

$\left\{\begin{array}{c}x={(1-\mathrm{\&tau;})}^3{x}_0+3\mathrm{\&tau;}{(1-\mathrm{\&tau;})}^2{x}_1+3{\mathrm{\&tau;}}^2(1-\mathrm{\&tau;})x_2+{\mathrm{\&tau;}}^3{x}_f\\ y={(1-\mathrm{\&tau;})}^3{y}_0+3\mathrm{\&tau;}{(1-\mathrm{\&tau;})}^2{y}_1+3{\mathrm{\&tau;}}^2(1-\mathrm{\&tau;})y_2+{\mathrm{\&tau;}}^3{y}_f\\ z={(1-\mathrm{\&tau;})}^3{z}_0+3\mathrm{\&tau;}{(1-\mathrm{\&tau;})}^2{z}_1+3{\mathrm{\&tau;}}^2(1-\mathrm{\&tau;})z_2+{\mathrm{\&tau;}}^3{z}_f\end{array}\right.$

其中，x、y、z是拟合的飞行器的空间弹道上的任一位置坐标；参数τ∈[0,1]，τ＝0代表起点，τ＝1代表终点；(x₀,y₀,z₀)、(x₁,y₁,z₁)、(x₂,y₂,z₂)、(x_f,y_f,z_f)为控制弹道形状的四个控制点，(x₀,y₀,z₀)和(x_f,y_f,z_f)分别为弹道的起点和终点；四个控制点(x₀,y₀,z₀)、(x₁,y₁,z₁)、(x₂,y₂,z₂)和(x_f,y_f,z_f)映射在XOY平面上的点分别为(x₀,y₀)、(x₁,y₁)、(x₂,y₂)和(x_f,y_f)，则点(x₀,y₀)与(x₁,y₁)的连线与点(x₂,y₂)和(x_f,y_f)的连线的交点为(x_m,y_m)；设参数参数则中间两个控制点(x₁,y₁,z₁)和(x₂,y₂,z₂)用参数k₁和k₂来描述并获取：

$\left\{\begin{array}{c}{x}_1=x_0+k_1(x_m-x_0)\\ y_1=y_0+\frac{\tan {\mathrm{\&gamma;}}_0}{\cos {\mathrm{\&psi;}}_0}(x_1-x_0)\\ z_1=z_0-\tan {\mathrm{\&psi;}}_0(x_1-x_0)\end{array}\right.;$ $\left\{\begin{array}{c}{x}_2=x_m+k_2(x_f-x_m)\\ y_2=y_f-\frac{\tan {\mathrm{\&gamma;}}_f}{\cos {\mathrm{\&psi;}}_f}(x_f-x_2)\\ z_2=z_f+\tan {\mathrm{\&psi;}}_f(x_f-x_2)\end{array}\right.;$ k₁∈[0,1]，k₂∈[0,1]

进一步，确定弹道坐标系下的法向加速度的指令a_yB和侧向加速度的指令a_zB：

a_yB＝gcosγ+v²cosγcosψdγ/dx

a_zB＝-v²cos²γcosψdψ/dx

其中，g是重力加速度，γ是弹道倾角，v是飞行器的速度，ψ是弹道偏角；

其次，确定弹道坐标系下的法向加速度a_y和侧向加速度a_z：

$\left\{\begin{array}{c}{a}_y=a_{\mathrm{yB}}\left(0\right)\&CenterDot;n_{\&perp;\max }\&CenterDot;g/\sqrt{a_{\mathrm{yB}}{\left(0\right)}^2+a_{\mathrm{zB}}{\left(0\right)}^2}\\ a_z=a_{\mathrm{zB}}\left(0\right)\&CenterDot;n_{\&perp;\max }\&CenterDot;g/\sqrt{a_{\mathrm{yB}}{\left(0\right)}^2+a_{\mathrm{zB}}{\left(0\right)}^2}\end{array}\right.\mathrm{if}\sqrt{a_{\mathrm{yB}}{\left(0\right)}^2+a_{\mathrm{zB}}{\left(0\right)}^2}/g>n_{\&perp;\max }$

$\left\{\begin{array}{c}{a}_y=a_{\mathrm{yB}}\left(0\right)\\ a_z=a_{\mathrm{zB}}\left(0\right)\end{array}\right.\mathrm{if}\sqrt{a_{\mathrm{yB}}{\left(0\right)}^2+a_{\mathrm{zB}}{\left(0\right)}^2}/g\&le;n_{\&perp;\max }$

其中，a_yB(0)是在初始时刻所对应的法向加速度的指令，a_zB(0)是在初始时刻所对应的侧向加速度的指令，n_⊥max是飞行器的可用法向过载；

然后，获得倾侧角指令σ_c和攻角指令α_c： $\left\{\begin{array}{c}{\mathrm{\&sigma;}}_c=\arctan (a_z/a_y)\\ {\mathrm{\&alpha;}}_c={f_L}^{-1}(M_a,{\mathrm{ma}}_y/\left(\mathrm{qS}\cos \mathrm{\&sigma;}\right))\end{array}\right.;$ 其中，M_a为马赫数，M_a＝v/a，a为当地音速，m表示飞行器的质量，动压q＝ρv²2，ρ是大气密度，S是参考面积，f_L是升力系数C_L的表示函数；

最后，确定飞行器的倾侧角σ：

其中，是设定的常数，表示倾侧角对时间求导的最大数值；σ_min和σ_max是设定的常数，分别表示侧倾角需设定的最小值和最大值；σ_prev表示前一制导周期的倾侧角，t_prev表示前一制导周期的时间，t表示当前时间；

根据得到的倾侧角σ，能够确定攻角指令α_c，进一步确定飞行器的攻角α：

其中，是设定的常数，表示攻角对时间求导的最大数值；α_min和α_max是设定的常数，分别表示攻角需设定的最小值和最大值；α_prev表示前一制导周期的攻角；

第二步、通过优化计算，获得满足脱靶量、碰撞角和末端攻角要求的k₁和k₂的范围；

第三步、在第二步得到的范围内，找到末端速度最大和末端速度最小所对应的k₁和k₂值，根据给定的末端速度，确定参数k₁和k₂。

2.根据权利要求1所述的飞行器导引方法，其特征在于，所述的第三步中，在标称条件下，要实现期望的末速，首先固定参数k₁的值，然后采用Secant法根据下式来实时修正k₂的值：

$k_2(i+1)=k_2\left(i\right)+\frac{(v_{\mathrm{fdes}}-v_f\left(i\right))(k_2\left(i\right)-k_2(i-1))}{(v_f\left(i\right)-v_f(i-1))}$

k₂(i)表示第i次在线迭代过程中的k₂值，v_fdes表示期望的末速值，v_f(i)表示第i次仿真预测的末速值；

经过在线迭代，确定参数k₂，并收敛至期望的末速。