一种弹道规划方法与流程

本发明涉及一种弹道规划
技术领域：
，尤其涉及一种通过曲线外形计算弹道参数的草图交互虚拟域动态逆弹道规划方法。
背景技术：
：随着导弹防御系统全面部署的临近，各军事强国目前竞相研究具有一定机动能力的导弹以提高导弹突防能力从而保持导弹打击能力。然而这类导弹多在大气层内飞行，在满足发射约束与任务需求的同时，必须考虑能量、过载、热流、动压等约束，从而大大增加了导弹飞行过程复杂度。如何对飞行弹道进行规划，在满足复杂约束条件下提高作战效能，是此类导弹作战应用面临的重大挑战。常规弹道规划方法应用于此类导弹弹道规划时存在如下困难：(1)以控制量为弹道规划设计变量，通过调整控制量改变导弹受力状态，从而实现弹道修改，无法根据任务需求直接对弹道位置特征进行调整，导致弹道规划过程过于复杂；(2)在导弹受力状态改变之后，一般需要通过弹道积分获得导弹的飞行特性与约束条件的满足情况，由于弹道规划中积分运算需要大量重复迭代，导致弹道规划过程效率较低；(3)通过改变控制量间接改变弹道的方式，在实际应用中难以获得与目标弹道完全一致的规划结果，导致弹道规划过程无法全面反映设计者意愿。LuPing[1]将动态逆思想应用于制导领域，使制导算法更加简洁，提高了航天飞机上升段制导效率。动态逆方法由于能够从路径得到控制量，因而在机器人路径规划、战斗机和无人机航迹规划及巡航导弹弹道规划等领域受到关注，并取得初步的阶段成果[2]-[4]。由于动态逆思想需要将轨迹方程不断微分以得到速度、加速度等必要信息，所以要求系统方程各阶导数在时域中有特定含义，而这样做又难以满足“轨迹可调节”的意愿。对此，Yakimenko[5],[6]通过引入虚拟路径，将轨迹优化问题从时域转化到虚拟域，再用动态逆的方法优化轨迹，提出了虚拟域动态逆(InverseDynamicsintheVirtualDomain，IDVD)方法。该方法通过虚拟速度将虚拟域路径和时域路径对应起来。虚拟域动态逆方法首先将虚拟域参数微分，得到虚拟域速度、加速度和过载等信息，但这些信息没有时域的时间标签。接下来利用虚拟速度将虚拟域加速度和过载等信息对应到时域上，并加上时域的时间标签。这样做可以得到解析的参数求解方程，因此计算速度非常快。由于虚拟域轨迹的设计是通过多项式拟合的，虚拟轨迹并不能保证全局最优，进而也无法保证时域轨迹的全局最优。JohnA.Lukacs[7],IanD.Cowling[8],GeorgeBoyarko等[9]分别与Yakimenko合作，研究了虚拟域动态逆方法在四旋翼无人机轨迹规划和制导、拦截弹在线制导、飞行器交会等的应用。NAkhtar等[10]研究了虚拟域动态逆方法在滑翔机实时轨迹规划仿真中的应用，结果表明该方法计算速度快，轨迹规划方案有利于节省燃料。闫梁等[3]使用变节点算法对虚拟域动态逆方法进行了改进，在节点总数不变的情况下提高了运算效率，但算法在大部分情况下，可能需经多次迭代积分确定节点位置，运算复杂度高。ZhangYu等[2]将虚拟域动态逆方法应用在无人战斗机的轨迹规划仿真中，获得了较高的计算速度和较好的收敛特性。MarcoCiarcià等[11]基于虚拟域动态逆方法提出了一种近忧制导策略，并应用于航天器交会对接地面试验中，通过对比实验证明了方法的有效性。JacopoVentura等[12]研究了快速最优姿态轨迹生成问题，评估了不同情况下虚拟域动态逆方法的计算性能。技术实现要素：传统弹道规划方法应用于具有变向能力的导弹，例如助推滑翔导弹、有机动能力的弹道导弹和巡航导弹等存在的不适应性。如现有虚拟域动态逆方法在最优控制理论框架下提出，其目的是为了解决制导问题，采用了处理两点边值问题的方法解决轨迹规划问题，一般能够对起点和终点的位置参数进行调节，但对过程参数如飞行轨迹不能进行人机交互调整，因此不具备对弹道的交互调控能力。针对现有技术存在的以上缺陷，本发明提出了一种弹道规划方法。该方法引入人的智能因素，基于解析曲线直接设计弹道特征，逆向计算控制量与约束条件，并通过草图交互技术将设计人员意愿与计算机规划过程有机结合，显著改善弹道规划过程的直观性，提高复杂约束下助推滑翔导弹弹道规划的效率。本发明采用的技术方案是：一种弹道规划方法，包括以下步骤：S1.建立助推滑翔导弹运动模型；设地平面为一平面，导弹飞行侧滑角为零，得到助推滑翔导弹三自由度动力学计算方程：x·1=Vcosγ·cosψx·2=Vcosγ·sinψx·3=-V·sinψV·=g(nx-sinγ)γ·=gV-1(ny-cosγ)ψ·=gnzV-1·cos-1γ---(1)]]>其中，x1、x2、x3分别表示助推滑翔导弹在北东地坐标系下沿三个坐标轴的位置，V为导弹飞行速度，γ为弹道倾角，ψ为弹道偏角，g为重力加速度；nx、ny、nz分别为弹道坐标系下沿坐标轴三个方向的过载：nx=Tcosα-Dmgny=Tsinα+Lmgnz=Zmg---(2)]]>其中，T为导弹发动机推力，α为导弹飞行攻角，D为导弹飞行阻力，L为导弹飞行升力，Z为导弹飞行侧向力，m为助推滑翔导弹的质量；T、α和m均为时间的函数，通过插值或拟合后的公式求得：T＝T(t)(3)α＝α(t)(4)m＝m(t)(5)速度系下导弹的阻力、升力、侧向力为D=12ρV2SCD---(6)]]>L=12ρV2SCL---(7)]]>Z=12ρV2SCZ---(8)]]>其中ρ为大气密度，S为特征面积，CD、CL、CZ分别为阻力系数、升力系数和侧向力系数，通过插值或拟合后的公式得到：CD＝CD(α,V,-x3)(9)CL＝CL(α,V,-x3)(10)CZ＝CZ(α,V,-x3)(11)S2.约束条件建模(1)驻点热流密度驻点热流密度qws采用Kemp-Riddell公式计算：qws=131884.2RN(ρ∞ρo)1/2(v∞vc)3.25(1-hw/hs)≤qwsmax---(12)]]>式中ρο＝1.225(kg/m3)；νc＝7900(m/s)；RN－驻点曲率半径(m)；qws－驻点热流密度(kW/m2)；hw－壁面焓值(kJ/kg)；hs－驻点焓值(kJ/kg)；qwsmax－最大驻点热流密度(kW/m2)；ρ∞－来流密度；ν∞－来流速度；(2)过载nnx≤nxmaxny≤nymaxnz≤nzmax---(13)]]>其中nxmax、nymax、nzmax为助推滑翔导弹在弹道坐标系下三个方向可承受的最大过载；(3)动压qq=12ρV2≤qmax---(14)]]>qmax为助推滑翔导弹可承受的最大动压；(4)控制量当仅考虑助推滑翔导弹纵平面运动时，取导弹飞行攻角α为控制量，需要对其幅值进行约束，如下：αmin≤α≤αmax(15)αmin为最小攻角，αmax为最大攻角；(5)起始条件和终端约束助推滑翔导弹的初始位置x1=x10,x2=x20,x3=x30---(16)]]>和为给定的发射点坐标；S3.虚拟域草图快速生成方法；S3.1助推滑翔导弹垂直上升段草图生成方法取发射点坐标为原点，即O(0,0,0)，那么可通过控制点V0(0,0,V0z)来定义垂直上升段轨迹，其中V0z为控制点V0的z轴坐标，O和V0即为垂直上升段弹道的控制点；不失一般性，设τ∈[a,a+1)，其中a＝0,1,2,…，τ为虚拟域的自变量即虚拟域时间，那么在北东地坐标系Oxyz下，虚拟域轨迹为：x(τ)=0y(τ)=0z(τ)=V0z(τ-a)2---(17)]]>其中x(τ)、y(τ)、z(τ)分别为北东地坐标系Oxyz中x轴、y轴和z轴的坐标；然后得到x(τ)、y(τ)、z(τ)在虚拟域的一阶导数和二阶导数：x′(τ)=0y′(τ)=0z′(τ)=2V0z(τ-a)x′′(τ)=0x′′(τ)=0z′′(τ)=2V0z---(18)]]>S3.2助推滑翔导弹转弯及滑翔段草图生成方法助推滑翔导弹转弯及滑翔段的草图可通过C1连续曲线、C2连续曲线或者C3连续曲线生成；S4.草图交互虚拟域动态逆弹道参数求解方法；通过S3获得一条完整曲线后，可将曲线表示成为只含有一个自变量的函数r(τ)，以虚拟域时刻τ为自变量，τ的取值范围是[0,n]，其中n为曲线段数；假设有函数f(·)，τ和t分别为虚拟域时刻和时域时刻，τ和t有数学关系：df(τ)dt=df(τ)dτdτdt---(34)]]>在此定义λ(τ)=dτdt---(35)]]>λ(τ)为虚拟速度，并约定表示f(·)在时域内求导，即f′(·)表示f(·)在虚拟域求导，即根据(34)、(35)两式有：f·(τ)=f′(τ)λ(τ)---(36)]]>进而有：f··(τ)=d(f′(τ)λ(τ))dτdτdt=f′′(τ)λ2(τ)+f′(τ)λ(τ)λ′(τ)---(37)]]>确定虚拟域曲线r(τ)后，在τ的取值范围取N个节点将虚拟域曲线r(τ)离散；若τ的取值范围为[0,τf]，则可通过如下方法获得节点之间在虚拟域的时间间隔Δτ=τfN-1---(38)]]>每一个节点对应的虚拟域时刻可表示为τj＝τj-1+Δτ,j＝2,…,N(39)其中N为偶数；助推滑翔导弹采用固体发动机作为动力，时间-质量关系m(t)、时间-推力关系T(t)视为已知；首先求出当前节点处虚拟域基本状态xi、x′i、x″i(i＝1,2,3)，由r(τ)以及其一阶、二阶导数得到虚拟域基本状态，整理后结果形式如下：x1,j-1x1,j-1′x1,j-1′′x2,j-1x2,j-1′x2,j-1′′x3,j-1x3,j-1′x3,j-1′′=r(τj-1)r′(τj-1)r′′(τj-1)---(40)]]>当τj-1∈[0,1)时r(τj-1)＝r1(τj-1)，当τj-1∈[1,2]时r(τj-1)＝r2(τj-1)，其中xi,j-1表示虚拟域基本状态xi在节点(j-1)处的值，即xi(τj-1)；同样的x′i,(j-1)表示虚拟域基本状态x′i在节点(j-1)处的值，x″i,(j-1)表示虚拟域基本状态x″i在节点(j-1)处的值；其中各节点处的状态量和控制量的计算方法如下：(1)弹道倾角γj-1=arctan-x3,j-1′x1,j-1′2+x2,j-1′2---(41)]]>(2)弹道偏角ψj-1=arctanx2,j-1′x1,j-1′---(42)]]>(3)弹道倾角虚拟域导数γ′j-1和弹道偏角虚拟域导数ψ′j-1γj-1′=x3,j-1′(x1,j-1′x1,j-1′′+x2,j-1′x2,j-1′′)-x3,j-1′′(x1,j-1′2+x2,j-1′2)(x1,j-1′2+x2,j-1′2)32cos2γj-1---(43)]]>ψj-1′=x2,j-1′′x1,j-1′-x2,j-1′x1,j-1′′x1,j-1′2cos2ψj-1---(44)]]>(4)y方向过载ny,j-1=Vj-1g-1λj-1x3,j-1′(x1,j-1′x1,j-1′′+x2,j-1′x2,j-1′′)-x3,j-1′′(x1,j-1′2+x2,j-1′2)(x1,j-1′2+x2,j-1′2)32·cos2γj-1+cosγj-1---(45)]]>(5)z方向过载nz,j-1=Vj-1g-1λj-1x2,j-1′′x1,j-1′-x2,j-1′x1,j-1′′x1,j-1′2cos2ψj-1cosγj-1---(46)]]>(6)大气密度hj-1＝-x3,j-1(47)ρj-1＝fair(hj-1)(48)其中fair(hj-1)为大气密度函数；(7)升力采用攻角的一次函数表示升力系数即：CL,j-1＝l1αj-1+l2(49)那么升力Lj-1=12ρj-1Vj-12S(l1αj-1+l2)---(50)]]>其中S为导弹特征面积，假设这一数值在飞行过程中不变；其中l1和l2为拟合系数；(8)攻角αj-1=ny,j-1mj-1g-12ρj-1Vj-12Sl2Tj-1+12ρj-1Vj-12Sl1---(51)]]>(9)阻力nx,j-1=Tj-1cosαj-1-12ρj-1Vj-12S(d1Vj-12+d2Vj-1+Vj-1)mj-1g---(52)]]>其中CD,j-1=d1Vj-12+d2Vj-1+Vj-1---(53)]]>Dj-1=12ρj-1Vj-12S(d1Vj-12+d2Vj-1+Vj-1)---(54)]]>d1和d2为拟合系数；(8)速度的虚拟域导数V′(τ)=V·(τ)/λ(τ)---(55)]]>(9)下一节点处速度Vj＝Vj-1+g(nx,j-1-sinγj-1)/λj-1Δτ(56)(10)节点(j-1)到节点j的飞行时间Δtj=Σi=13(xi,j-xi,j-1)2(Vj+Vj-1)/2---(57)]]>(11)下一节点时域时刻tj＝tj-1+Δtj(58)(12)Δt2的估计值Δt~2=TboostN-1---(59)]]>(13)λ1的估计值λ~1=ΔτΔt2≈ΔτΔt~2=τf(N-1)(N-1)Tboost=τfTboost.---(60)]]>本发明S3.2中，助推滑翔导弹转弯及滑翔段的草图通过C1连续曲线生成的方法如下：首先定义下列3个关于t的多项式为带形状参数λb和μb的三次调配函数：X0(t)=1-λbt+(2λb-3)t2+(2-λb)t3X1(t)=λbt+(μb-2λb)t2+(λb-μb)t3X2(t)=(3-μb)t2+(μb-2)t3,t∈[0,1]---(19)]]>其中0≤λb,μb≤3；然后给定二维或三维空间中3个控制顶点Vi(i＝0,1,2)，称曲线r(t)为可调控的三次参数曲线：r(t)＝X0(t)V0+X1(t)V1+X2(t)V2,t∈[0,1](20)其中Xi(t)(i＝0,1,2)为按式(19)定义的调配函数；或者，将上式(20)写成矩阵形式为：r(t)=1tt2t3100-λbλb02λb-3μb-2λb3-λb2-λbλb-μbμb-2V0V1V2,t∈[0,1].---(21)]]>本发明S3.2中，助推滑翔导弹转弯及滑翔段的草图通过C2连续曲线生成的方法如下：给定一组控制点<V0,V1,…,Vn＞，以其构成控制切线多边形，取重构控制点：b2i=Vib2i+1=λciVi+(1-λci)Vi+1,i=0,1,...,n-1---(22)]]>式(22)中λci为切点控制参数，且0＜λci＜1；那么由重构控制点b0,b1,…,b2n构成的曲线为：ri(t)=Σj=03Bj(t)bi+j,i=0,1,...,2n-2,0≤t≤αc---(23)]]>其中αc为调节参数，且0≤αc≤π；B0(t),B1(t),B2(t),B3(t)为C-B样条基函数，分别定义如下：B0(t)=(αc-t)-sin(αc-t)2αc(1-C)B3(t)=t-sin(t)2αc(1-C)B1(t)=B3(t)-2B0(t)+2(αc-t)(1-C)2αc(1-C)B2(t)=B0(t)-2B3(t)+2t(1-C)2αc(1-C)---(24)]]>其中C＝cosαc，S＝sinαc；或者，将C2连续的C-B样条曲线写成矩阵形式为：ri(t)=Σj=03Bj(t)bi+j=12αc(1-C)sintcostt1]]>C-(1+2C)2+C-1-S2S-S0-11+2C-(1+2C)1αc-2αcCαc0bibi+1bi+2bi+3---(25)]]>其中i＝0,1,…,2n。本发明S3.2中，助推滑翔导弹转弯及滑翔段的草图通过C3连续曲线生成的方法如下：给定一组控制点<V0,V1,…,Vn>，以其构成控制多边形，构造边矢量：ai＝Vi-Vi-1,i＝1,2,…,n(27)则可计算控制多边形每个顶点处的切矢量T0=t0(-a2)+(1-t0)a1Ti=tiai+(1-ti)ai+1Tn=tnan+(1-tn)(-an-1),i=1,2,...,n-1---(28)]]>令Vi*=Vi+|ai+1×Ti+1||Ti×Ti+1|Ti,i=1,2,...,n-1---(29)]]>及l0=|V0*-V0|li=min{|Vi*-Vi|,|Vi-Vi-1*|}ln=|Vn-Vn-1*|,i=1,2,...,n-1---(30)]]>取重构控制点d4i=Vi-λiliTi|Ti|d4i+1=Vi-12λiliTi|Ti|d4i+2=Vi+12λiliTi|Ti|d4i+3=Vi+λiliTi|Ti|,i=0,1,...,n-1---(31)]]>其中λi为控制调节参数，且0＜λi＜1；那么由重构控制点可构成四次B样条曲线ri(t)=Σj=04Ej(t)di+j,i=0,1,...,4n-1,0≤t≤1---(32)]]>其中E0(t)=(t4-4t3+6t2-4t+1)/24E1(t)=(-4t4+12t3-6t2-12t+11)/24E2(t)=(6t4-12t3-6t2+12t+11)/24E3(t)=(-4t4+4t3+6t2+4t+1)/24E4(t)=t4/24.]]>本发明的有益效果在于：(1)本发明脱离最优控制理论框架，基于完整的解析曲线直接设计弹道特征，逆向计算控制量与约束条件，不仅能够对起点和终点的位置参数进行调节，更能对过程参数如飞行轨迹进行调整。与现有最好技术相比，本发明不局限于解决导弹的制导问题，除此之外还可解决导弹的弹道规划、弹道设计、实现弹道能力评估等问题。(2)本发明将人类智慧通过草图弹道的形式引入算法中，与现有技术相比，实现了人与计算机的高效交互沟通，极大方便了弹道规划者利用计算机实现本人意愿，降低了因计算机理解偏差导致复算的可能性，提高了计算效率。(3)现有技术是通过调整导弹飞行的控制量实现对弹道的调整，调整过程复杂、直观性差，难以进行形式复杂弹道的调整或规划。本发明实现了规划者对弹道的直接调整，简化了调整过程、直观性强，易于进行形式复杂弹道的调整或规划。(4)本发明使用解析的曲线不断微分得到弹道参数，与目前常用的积分法相比省去了积分环节，通过仿真实验验证，每条弹道的平均计算时长为0.013s，较积分法提高三个数量级。附图说明图1为本发明的实现过程示意图；图2为本发明的流程图；具体实施方式下面结合附图和具体实施方式对本发明作进一步说明。一种弹道规划方法，包括以下步骤：S1.建立助推滑翔导弹运动模型；设地平面为一平面，导弹飞行侧滑角为零，可得到助推滑翔导弹三自由度动力学计算方程：x·1=Vcosγ·cosψx·2=Vcosγ·sinψx·3=-V·sinψV·=g(nx-sinγ)γ·=gV-1(ny-cosγ)ψ·=gnzV-1·cos-1γ---(1)]]>其中，x1、x2、x3分别表示助推滑翔导弹在北东地坐标系下沿三个坐标轴的位置，V为导弹飞行速度，γ为弹道倾角，ψ为弹道偏角，g为重力加速度，文中统一取为9.8m/s2；nx、ny、nz分别为弹道坐标系下沿坐标轴三个方向的过载：nx=Tcosα-Dmgny=Tsinα+Lmgnz=Zmg---(2)]]>其中，T为导弹发动机推力，α为导弹飞行攻角，D为导弹飞行阻力，L为导弹飞行升力，Z为导弹飞行侧向力，m为助推滑翔导弹的质量。T、α和m均为时间的函数，一般通过插值或拟合后的公式求得：T＝T(t)(3)α＝α(t)(4)m＝m(t)(5)速度系(定义见文献[15])下导弹的阻力、升力、侧向力为D=12ρV2SCD---(6)]]>L=12ρV2SCL---(7)]]>Z=12ρV2SCZ---(8)]]>其中ρ为大气密度，S为特征面积，CD、CL、CZ分别为阻力系数、升力系数和侧向力系数，它们与攻角α、飞行速度V和飞行高度h(h＝-x3)有关，一般在进行实验后制作成插值表，通过插值或拟合后的公式得到：CD＝CD(α,V,-x3)(9)CL＝CL(α,V,-x3)(10)CZ＝CZ(α,V,-x3)(11)S2.约束条件建模(1)驻点热流密度因为驻点是助推滑翔导弹加热较为严重的部分，因此通常作为弹道规划的重要约束条件。本发明中采用Kemp-Riddell公式计算：qws=131884.2RN(ρ∞ρo)1/2(v∞vc)3.25(1-hw/hs)≤qwsmax---(12)]]>式中ρο＝1.225(kg/m3)；νc＝7900(m/s)；RN－驻点曲率半径(m)；qws－驻点热流密度(kW/m2)；hw－壁面焓值(kJ/kg)；hs－驻点焓值(kJ/kg)；qwsmax－最大驻点热流密度(kW/m2)；ρ∞－来流密度；ν∞－来流速度；(2)过载n为使导弹结构在飞行过程中不发生破坏，一般需要对过载进行约束。nx≤nxmaxny≤nymaxnz≤nzmax---(13)]]>其中nxmax、nymax、nzmax为助推滑翔导弹在弹道坐标系下三个方向可承受的最大过载。(3)动压q动压一定程度上反映助推滑翔导弹飞行环境的优劣，一般也作为重要的约束条件给出。q=12ρV2≤qmax---(14)]]>qmax为助推滑翔导弹可承受的最大动压。(4)控制量当仅考虑助推滑翔导弹纵平面运动时，一般取导弹飞行攻角α为控制量，这里需要对其幅值进行约束。αmin≤α≤αmax(15)αmin为最小攻角，αmax为最大攻角。(5)起始条件和终端约束根据任务要求，可以对位置、速度、速度倾角、攻角等进行约束。大多数情况下，弹道规划问题中都含有助推滑翔导弹的初始位置这一起始条件。x1=x10,x2=x20,x3=x30---(16)]]>和为给定的发射点坐标。上述(1)至(5)约束条件中，不等式两端的值，如qwsmax、αmin等均为导弹参数，由导弹的飞行性能决定，在导弹定型后，这些参数也随之确定，可视为已知。S3.虚拟域草图快速生成方法；(1)助推滑翔导弹垂直上升段草图生成方法取发射点坐标为原点，即O(0,0,0)，那么可通过控制点V0(0,0,V0z)来定义垂直上升段轨迹，其中V0z为控制点V0的z轴坐标，O和V0即为垂直上升段弹道的控制点。控制点为控制一条曲线外形的关键节点，可以通过改变它们的位置调整曲线形状，从而构建所需的曲线。不失一般性，设τ∈[a,a+1)(a＝0,1,2,…)为虚拟域的自变量，即虚拟域时间，那么在北东地坐标系Oxyz下，虚拟域轨迹为：x(τ)=0y(τ)=0z(τ)=V0z(τ-a)2---(17)]]>其中x(τ)、y(τ)、z(τ)分别为北东地坐标系Oxyz中x轴、y轴和z轴的坐标。然后得到x(τ)、y(τ)、z(τ)在虚拟域的一阶导数和二阶导数：x′(τ)=0y′(τ)=0z′(τ)=2V0z(τ-a)x′′(τ)=0x′′(τ)=0z′′(τ)=2V0z---(18)]]>(2)助推滑翔导弹转弯及滑翔段草图生成方法首先对曲线连续性作必要说明。设有m段多项式曲线构成一条完整曲线C(τ)，其中第i条多项式曲线为Ci(τ),i＝1,2,…,m.若为Ci(τ)的j阶导矢，τi为Ci(τ)的末端点自变量，如果0≤j≤k对每段曲线都成立，则称曲线C(τ)在τi处是Ck连续的[23],[24]。转弯及滑翔段的草图可通过C1连续曲线、C2连续曲线或者C3连续曲线生成，下面分别介绍它们的生成方法。(a)C1连续曲线草图生成方法贝塞尔曲线[25]具有操作直观，连接平滑等特点，适于解决本发明中的弹道表示问题。这里介绍其中一类可调控的三次多项式曲线。定义1下列3个关于t的多项式为带形状参数λb和μb的三次调配函数：X0(t)=1-λbt+(2λb-3)t2+(2-λb)t3X1(t)=λbt+(μb-2λb)t2+(λb-μb)t3X2(t)=(3-μb)t2+(μb-2)t3,t∈[0,1]---(19)]]>其中0≤λb,μb≤3。定义2给定二维或三维空间中3个控制顶点Vi(i＝0,1,2)，称曲线r(t)为可调控的三次参数曲线：r(t)＝X0(t)V0+X1(t)V1+X2(t)V2,t∈[0,1](20)其中Xi(t)(i＝0,1,2)为按式(19)定义的调配函数。上式(20)改写成矩阵形式为：r(t)=1tt2t3100-λbλb02λb-3μb-2λb3-λb2-λbλb-μbμb-2V0V1V2,t∈[0,1]---(21)]]>(b)C2连续曲线草图生成方法给定一组控制点<V0,V1,…,Vn＞，以其构成控制切线多边形，取重构控制点[31][32]b2i=Vib2i+1=λciVi+(1-λci)Vi+1,i=0,1,...,n-1---(22)]]>式(22)中λci为切点控制参数[33]，且0＜λci＜1。那么由重构控制点b0,b1,…,b2n构成的曲线为：ri(t)=Σj=03Bj(t)bi+j,i=0,1,...,2n-2,0≤t≤αc---(23)]]>其中αc为调节参数，且0≤αc≤π。B0(t),B1(t),B2(t),B3(t)为C-B样条基函数，分别定义如下：B0(t)=(αc-t)-sin(αc-t)2αc(1-C)B3(t)=t-sin(t)2αc(1-C)B1(t)=B3(t)-2B0(t)+2(αc-t)(1-C)2αc(1-C)B2(t)=B0(t)-2B3(t)+2t(1-C)2αc(1-C)---(24)]]>这里C＝cosαc，S＝sinαc。也可将C2连续的C-B样条曲线写成矩阵形式ri(t)=Σj=03Bj(t)bi+j=12αc(1-C)sintcostt1]]>C-(1+2C)2+C-1-S2S-S0-11+2C-(1+2C)1αc-2αcCαc0bibi+1bi+2bi+3---(25)]]>其中i＝0,1,…,2n。若要使曲线通过端点处的控制点可以令初始两个控制点和结尾末尾两个控制点重合[34]。V0=V1Vn-1=Vn,n≥4---(26)]]>(c)C3连续草图生成方法给定一组控制点<V0,V1,…,Vn>，以其构成控制多边形[35],[36]。构造边矢量ai＝Vi-Vi-1,i＝1,2,…,n(27)则可以计算控制多边形每个顶点处的切矢量T0=t0(-a2)+(1-t0)a1Ti=tiai+(1-ti)ai+1Tn=tnan+(1-tn)(-an-1),i=1,2,...,n-1---(28)]]>令Vi*=Vi+|ai+1×Ti+1||Ti×Ti+1|Ti,i=1,2,...,n-1---(29)]]>及l0=|V0*-V0|li=min{|Vi*-Vi|,|Vi-Vi-1*|}ln=|Vn-Vn-1*|,i=1,2,...,n-1---(30)]]>取重构控制点d4i=Vi-λiliTi|Ti|d4i+1=Vi-12λiliTi|Ti|d4i+2=Vi+12λiliTi|Ti|d4i+3=Vi+λiliTi|Ti|,i=0,1,...,n-1---(31)]]>其中λi为控制调节参数，且0＜λi＜1。那么由重构控制点可构成四次B样条曲线ri(t)=Σj=04Ej(t)di+j,i=0,1,...,4n-1,0≤t≤1---(32)]]>其中E0(t)=(t4-4t3+6t2-4t+1)/24E1(t)=(-4t4+12t3-6t2-12t+11)/24E2(t)=(6t4-12t3-6t2+12t+11)/24E3(t)=(-4t4+4t3+6t2+4t+1)/24E4(t)=t4/24---(33)]]>S4.草图交互虚拟域动态逆弹道参数求解方法；通过S3获得一条完整曲线后，可将曲线表示成为只含有一个自变量的函数r(τ)，以虚拟域时刻τ为自变量，τ的取值范围是[0,n]，(n为曲线段数)，它不是时域内的时刻。r(τ)虽然能够反映导弹飞行轨迹的形状，但不能描述轨迹上每一点状态与时域内时刻的对应关系。下面根据之前介绍的方法先行建立轨迹上每一点状态与时域内时刻的对应关系。假设有函数f(·)，τ和t分别为虚拟域时刻和时域时刻，τ和t有数学关系：df(τ)dt=df(τ)dτdτdt---(34)]]>这里定义λ(τ)=dτdt---(35)]]>λ(τ)为虚拟速度，由于虚拟域时刻和时域时刻是关于计算节点(下文简称节点)的单调递增函数，故λ(τ)＞0，即λ(τ)具有正定性。并约定表示f(·)在时域内求导，即f′(·)表示f(·)在虚拟域求导，即下文中其它函数任意阶数的导数以此类推。根据(34)、(35)两式有：f·(τ)=f′(τ)λ(τ)---(36)]]>更进一步有：f··(τ)=d(f′(τ)λ(τ))dτdτdt=f′′(τ)λ2(τ)+f′(τ)λ(τ)λ′(τ)---(37)]]>确定虚拟域曲线r(τ)后，为方便数值计算，一般的做法是对连续函数进行离散，这里是在τ的取值范围取N个节点将曲线离散。节点数目可自行定义，但需注意节点数目增大时计算量也随之增大。本发明中为方便助推段分段取N为偶数。若τ的取值范围为[0,τf]，则可通过如下方法获得节点之间在虚拟域的时间间隔Δτ=τfN-1---(38)]]>每一个节点对应的虚拟域时刻可表示为τj＝τj-1+Δτ,j＝2,…,N(39)其中N为偶数。助推滑翔导弹采用固体发动机作为动力，时间-质量关系m(t)、时间-推力关系T(t)视为已知。首先求出当前节点处虚拟域基本状态xi、x′i、x″i(i＝1,2,3)。由r(τ)以及其一阶、二阶导数得到虚拟域基本状态，整理后结果形式如下：x1,j-1x1,j-1′x1,j-1′′x2,j-1x2,j-1′x2,j-1′′x3,j-1x3,j-1′x3,j-1′′=r(τj-1)r′(τj-1)r′′(τj-1)---(40)]]>当τj-1∈[0,1)时r(τj-1)＝r1(τj-1)，当τj-1∈[1,2]时r(τj-1)＝r2(τj-1)。其中x1,j-1表示x1在节点(j-1)处的值，即x1(τj-1)，其余类似下标写法意义相同，本文在不致混淆的前提下均在用此记法。另外需要说明，当前节点(j-1)对应时域时刻tj-1，在计算上一个节点(j-2)各状态量和控制量时(本发明将这些量统称节点变量)给出，下面介绍具体方法。(1)弹道倾角γj-1=arctan-x3,j-1′x1,j-1′2+x2,j-1′2---(41)]]>(2)弹道偏角ψj-1=arctanx2,j-1′x1,j-1′---(42)]]>(3)弹道倾角虚拟域导数γ′j-1和弹道偏角虚拟域导数ψ′j-1由于此步较为繁琐，为使算式更加清晰，省略γ、xi的下标(j-1)。γj-1′=x3,j-1′(x1,j-1′x1,j-1′′+x2,j-1′x2,j-1′′)-x3,j-1′′(x1,j-1′2+x2,j-1′2)(x1,j-1′2+x2,j-1′2)32cos2γj-1---(43)]]>ψj-1′=x2,j-1′′x1,j-1′-x2,j-1′x1,j-1′′x1,j-1′2cos2ψj-1---(44)]]>(4)y方向过载ny,j-1=Vj-1g-1λj-1x3,j-1′(x1,j-1′x1,j-1′′+x2,j-1′x2,j-1′′)-x3,j-1′′(x1,j-1′2+x2,j-1′2)(x1,j-1′2+x2,j-1′2)32·cos2γj-1+cosγj-1---(45)]]>(5)z方向过载nz,j-1=Vj-1g-1λj-1x2,j-1′′x1,j-1′-x2,j-1′x1,j-1′′x1,j-1′2cos2ψj-1cosγj-1---(46)]]>(6)大气密度使用杨炳尉的大气参数模型[38]得到大气密度ρj-1。hj-1＝-x3,j-1(47)ρj-1＝fair(hj-1)(48)其中fair(hj-1)为大气密度函数，已知高度可直接求取大气密度。(7)升力本文采用攻角的一次函数表示升力系数即：CL,j-1＝l1αj-1+l2(49)那么升力Lj-1=12ρj-1Vj-12S(l1αj-1+l2)---(50)]]>其中S为导弹特征面积，假设这一数值在飞行过程中不变。其中l1和l2为拟合系数。(8)攻角αj-1=ny,j-1mj-1g-12ρj-1Vj-12Sl2Tj-1+12ρj-1Vj-12Sl1---(51)]]>(9)阻力nx,j-1=Tj-1cosαj-1-12ρj-1Vj-12S(d1Vj-12+d2Vj-1+Vj-1)mj-1g---(52)]]>其中CD,j-1=d1Vj-12+d2Vj-1+Vj-1---(53)]]>Dj-1=12ρj-1Vj-12S(d1Vj-12+d2Vj-1+Vj-1)---(54)]]>d1和d2为拟合系数。(8)速度的虚拟域导数V′(τ)=V·(τ)/λ(τ)---(55)]]>(9)下一节点处速度Vj＝Vj-1+g(nx,j-1-sinγj-1)/λj-1Δτ(56)其中Δτ已通过(38)式给出。(10)节点(j-1)到节点j的飞行时间Δtj=Σi=13(xi,j-xi,j-1)2(Vj+Vj-1)/2---(57)]]>(11)下一节点时域时刻tj＝tj-1+Δtj(58)(12)Δt2的估计值Δt~2=TboostN-1---(59)]]>(13)λ1的估计值λ~1=ΔτΔt2≈ΔτΔt~2=τf(N-1)(N-1)Tboost=τfTboost.---(60)]]>实际计算时，λ1可在附近取值。在此附上本发明涉及的参考文献：[1]LuP.InverseDynamicsApproachtoTrajectoryOptimizationforanAerospacePlane[J].JournalofGuidance,Control,andDynamics,1993,16(4):726-732.[2]ZhangYu,ChenJing,ShenLincheng.Real-timetrajectoryplanningforUCAVair-to-surfaceattackusinginversedynamicsoptimizationmethodandrecedinghorizoncontrol[J].ChineseJournalofAeronautics,2013,26(4):1038–1056.[3]闫梁,李辕,赵继广,杜小平.基于变节点虚拟域动态逆的轨迹实时优化[J].航空学报,2013,34(12):2794-2803.[4]陈法龙.高超声速滑翔飞行器弹道快速规划研究[D].长沙:国防科学技术大学研究生院,2012.[5]YakimenkoO.DirectMethodforRapidPrototypingofNear-OptimalAircraftTrajectories[J].JournalofGuidanceControlandDynamics,2002,3(5):865-875.[6]YakimenkoO.A.,XuY.J.,GarethB.,ComputingShort-TimeAircraftManeuversUsingDirectMethodsAIAAGuidance,NavigationandControlConferenceandExhibit,Honolulu,Hawaii,18-21August2008,AIAA2008-6632.[7]LukacsJ.A.,YakimenkoO.A.Trajectory-Shape-VaryingMissileGuidanceforInterceptionofBallisticMissilesduringtheBoostPhase[C].AIAAGuidance,NavigationandControlConferenceandExhibit,HiltonHead,SouthCarolina,20-23August2007,AIAA2007-6538.[8]CowlingI.D.,YakimenkoO.A.,WhidborneJ.F.,CookeA.K.DirectMethodBasedControlSystemforanAutonomousQuadrotor[J].IntellRobotSyst,2010,60:285–316.[9]GeorgeB.,YakimenkoO.A.,MarcelloR.Real-Time6DoFGuidanceForofSpacecraftProximityManeuveringandCloseApproachwithaTumblingObject[C].AIAA/AASAstrodynamicsSpecialistConferenceToronto,OntarioCanada,2-5August2010,AIAA2010-7666.[10]AkhtarN.,WhidborneJ.F.,CookeA.K.Real-timeoptimaltechniquesforunmannedairvehiclesfuelsaving[J].PartG:AerospaceEngineering,2011,226:1315-1328.[11]MarcoC.,AlessioG.,MarcelloR.Anear-optimalguidanceforcooperativedockingmaneuvers[J].ActaAstronautica,2014,102:367–377.[12]JacopoV.,MarcelloR.,UlrichW.Performanceevaluationoftheinversedynamicsmethodforoptimalspacecraftreorientation[J].ActaAstronautica,2015,110:266-278.[13]YakimenkoO.DirectMethodforRapidPrototypingofNear-OptimalAircraftTrajectories[J].JournalofGuidanceControlandDynamics,2002,3(5):865-875.[14]钱杏芳,林瑞雄,赵亚男.导弹飞行力学[M].北京:北京理工大学出版社,2008:28-77.[15]贾沛然,陈克俊,何力.远程火箭弹道学[M].长沙:国防科技大学出版社,1993:57-73.[16]谢愈.复杂约束下高超声速滑翔飞行器弹道规划方法研究[D].长沙:国防科学技术大学研究生院,2012.[17]YLi,TMHosedales,YZSong,SGong.Free-handsketchrecognitionbymulti-kernelfeaturelearning[J].ComputerVision&ImageUnderstanding,2015,137:1-11.[18]吴玲达,邓维,张友根等.在线草图识别研究综述[J].计算机应用研究,2015,32(6):1601-1607.[19]HammondT,PaulsonB.Recognizingsketchedmultistrokeprimitives[J].ACMTransonInteractiveIntelligentSystems,2011,4(1):1-34.[20]DelayeA,AnquetilE.Fuzzyrelativepositioningtemplatesforsymbolrecognition[C].Procofthe11thInternationalConferenceonDocumentAnalysisandRecognition.WashingtonDC:IEEEComputerSociety,2011:1220-122.[21]张友根,吴玲达,宋汉辰等.手绘非规则军标图形的结构化识别方法[J].国防科技大学学报,2013,35(3):72-77.[22]张友根,吴玲达,邓维等.融合时空上下文的手绘笔画图文分类[J].电子与信息学报,2013,35(1):113-118.[23]莫蓉,常智勇.计算机辅助几何造型技术(第二版)[M].北京:科学出版社,2009:27-91.[24]LesPiegl,WayneTiller.非均匀有理B样条(第二版)[M].北京:清华大学出版社,2010:32-33.[25]吴晓勤,韩旭里.三次Bézier曲线的扩展[J].工程图学学报,2005,26(6):98-102.[26]吴晓勤.带形状参数的Bézier曲线[J].中国图象图形学报,2006,11(2):269-274.[27]李军成.一类可调控的三次多项式曲线[J].计算机工程与科学,2010,32(4):52-61.[28]DelayeA,LiuChenglin.Contextualtext/non-textstrokeclassificationinonlinehandwrittennoteswithconditionalrandomfields[J].PatternRecognition,2014,47(3):959-968.[29]PetersonEJ,StahovichTF,DoiE,etal.Groupingstrokesintoshapesinhand-drawndiagrams[C].Procofthe24thAAAIConferenceonArtificialIntelligence.2010:974-979.[30]胡伟峰,赵江洪.用户期望意象驱动的汽车造型基因进化[J].机械工程学报,2011,47(16):176-181.[31]FangKui,CaiFang,TanJian-Rong.C2andC3ContinuousBéziersplinecurveWithgiventangentpolygon.JournalofComputer-AidedDesign&ComputerGraphics,2000,12(5):330-332.[32]王成伟.带有给定切线多边形的C2连续的C-B样条曲线[J].计算机辅助设计与图形学学报,2001,13(12):1-4.[33]赵江洪,谭浩,谭征宇等.汽车造型设计:理论、方法与应用[M].北京理工大学出版社,2010.[34]韩旭里,刘圣军.二次Bézier曲线的扩展[J].中南工业大学学报(自然科学版),2003,34(2):214-17.[35]FangKui.Closed(G2-Continuous)BéziercurveWithgiventangentpolygon.ChineseJournalofNumericalMathematics,1991,13(2):34-37.[36]王成伟,姚云.C3连续的保凸四次B样条插值曲线[J].北京服装学院学报(自然科学版),2002,22(1):66-70.[37]LHering.Closed(C2-andC3-Continuous)BézierandB-splinecurveWithgiventangentpolygons.CAD,1983,15(1):3-6.[38]杨炳尉.标准大气参数的公式表示[J].宇航学报,1983,(1).当前第1页1 2 3