一种基于平稳滑翔弹道解析解的快速弹道规划方法技术

本发明专利技术公开了一种基于平稳滑行弹道解析解的快速弹道规划方法，包括以下几个步骤：步骤1：滑翔段弹道规划问题建模；步骤2：滑翔段弹道规划变量设计；步骤3：滑翔段弹道解析解求解；步骤4：滑翔段弹道终端速度控制方案；步骤5：滑翔段弹道再入走廊调整方案；步骤6：滑翔段弹道规划初值生成；步骤7：滑翔段弹道规划流程设计。本发明专利技术提出以纵向机动加速度比例系数和横向机动加速度比例系数为滑翔段弹道规划变量，使得运动方程中弹道倾角、弹道偏角、高度、经度与纬度的微分方程中不含速度项；本发明专利技术获得了固定纵向机动加速度比例系数和横向机动加速度比例系数对应的滑翔段弹道解析解。

【技术实现步骤摘要】

本专利技术涉及一种基于平稳滑行弹道解析解的快速弹道规划方法，属于航天技术、武器

技术介绍
随着高超声速技术的快速发展，再入弹道规划成为了研究热点。滑翔段是再入弹道的主要构成部分，决定了再入飞行的航程和机动范围。滑翔段弹道快速规划不仅可以分析高超声速飞行器的性能，还可用于在线弹道规划和预测制导，具有较高的研究价值。由于滑翔段弹道具有飞行时间长，对规划变量高度敏感，且再入过程非线性约束强，使得滑翔弹道的可行域较窄，采用传统优化算法求解时往往需要很长的计算时间。尽管伪谱法、SQP等优化算法在计算效率上有较大的提高，其计算耗时依然在秒级以上，并且当过程约束碰边界时，其计算效率会大大降低，甚至影响到结果的收敛性。为了提高滑翔段弹道规划的速度，一些特殊的约束被引入滑翔段弹道规划，如阻力曲线、等热流曲线、等动压曲线、平衡滑翔条件等。这些约束的共同特点是将滑翔段弹道限定在某种特殊的平稳滑翔状态，从而降低了弹道规划问题的敏度。但由于在弹道规划时还需要进行数值积分，限制了规划速度的提高。为了摆脱对弹道的积分的依赖，动态逆被引入滑翔段弹道规划。该方法直接规划几何弹道形状，并通过调整弹道形状来满足约束要求。但由于滑翔段前后的环境变化剧烈，使得滑翔段全程弹道难以用低阶曲线拟合。
技术实现思路
本专利技术的目的是通过求解平稳滑翔弹道高精度解析解，获得不依赖弹道积分的滑翔段弹道快速规划方法，为分析高超声速飞行器性能和预测制导提供技术支持。一种基于平稳滑行弹道解析解的快速弹道规划方法，包括以下几个步骤：步骤1：滑翔段弹道规划问题建模；步骤2：滑翔段弹道规划变量设计；步骤3：滑翔段弹道解析解求解；步骤4：滑翔段弹道终端速度控制方案；步骤5：滑翔段弹道再入走廊调整方案；步骤6：滑翔段弹道规划初值生成；步骤7：滑翔段弹道规划流程设计。本专利技术的优点在于：(1)提出以纵向机动加速度比例系数和横向机动加速度比例系数为滑翔段弹道规划变量，使得运动方程中弹道倾角、弹道偏角、高度、经度与纬度的微分方程中不含速度项；(2)获得了固定纵向机动加速度比例系数和横向机动加速度比例系数对应的滑翔段弹道解析解，包括高度解析解、飞行距离解析解、弹道偏角解析解、经度解析解和纬度解析解；(3)给出了纵向机动加速度比例系数和横向机动加速度比例系数的初值生成方法；(4)提出以纵向机动加速度比例系数、横向机动加速度比例系数、横向反转弹道倾角和初始弹道倾角增量为滑翔段弹道规划变量，其中：纵向机动加速度比例系数对应滑翔段飞行距离；横向机动加速度比例系数对应终端速度；横向反转弹道倾角对应终端经纬度；初始弹道倾角增量对应过程约束，上述四者之间相互独立。(5)本专利技术采用滑翔弹道解析解规确定弹道，而积分仅仅用于终端速度校正，使得弹道规划具有非常快的速度。附图说明图1是平稳滑翔弹道规划建模流程图；图2是大圆坐标系；图3是与Vf的关系；图4是Δγ对过程约束的影响；图5是球面飞行距离估计；图6是kε与sf的关系；图7是滑翔段弹道快速规划方法流程；图8是纵向弹道解析解精度验证；图9是横向弹道解析解精度验证；图10是Bell解析解；图11是可达区域边界弹道攻角曲线族；图12是可达区域边界弹道倾侧角曲线族；图13是可达区域边界横向弹道曲线族；图14是可达区域边界纵向弹道曲线族；图15是可达区域边界弹道速度曲线族；图16是可达区域边界弹道热流密度曲线族；图17是可达区域边界弹道动压曲线族；图18是可达区域边界弹道过载曲线族；图19是滑翔段弹道规划可达区域；图20是平面规划弹道；图21是规划弹道攻角曲线；图22是规划弹道倾侧角曲线。具体实施方式下面将结合附图和实施例对本专利技术作进一步的详细说明。本专利技术在分析运动方程的耦合条件的基础上，提出了一种以固定纵向机动加速度比例系数和横向机动加速度比例系为控制变量的弹道模式。该弹道模式下，速度与其他状态变量之间解耦，从而获得了解耦的弹道形状微分方程，并利用该微分方程求解得到了高精度的三维弹道形状解析解。此基础上，构造了滑翔段弹道的快速规划方法：首先利用弹道形状解析解规划出满足位置约束的弹道；然后通过调整横向机动大小来满足终端速度约束；最后通过调整初始弹道倾角来满足过程约束要求。本专利技术是一种基于平稳滑行弹道解析解的快速弹道规划方法，规划方法建模流程如图1所示，规划流程如图7所示，包括以下几个步骤：步骤1：滑翔段弹道规划问题建模为了便于理论分析，假设地球为圆球体，不考虑地球自转的影响，则半速度坐标系下的三自由度质点运动方程如下所示， h · = V sin γ θ · = V cos γ sin ψ / ( r cos φ ) ]]> s · = V cos γ R 0 / r φ · = V cos γ cos ψ / r ]]> 本文档来自技高网...

【技术保护点】
一种基于平稳滑行弹道解析解的快速弹道规划方法，包括以下几个步骤：步骤1：滑翔段弹道规划问题建模设地球为圆球体，不考虑地球自转，则半速度坐标系下的三自由度质点运动方程如下所示，h.=Vsinγθ.=Vcosγsinψ/(r cosφ)]]>s.=VcosγR0/rφ.=Vcosγcosψ/r]]>                                                   (1)γ.=[Lcosσ+(V2/r-g)cosγ]/Vψ.=[Lsinσ/cosγ+(V2/r)cosγsinψtanφ]/V]]>V.=-D-gsinγ]]>式中，h为飞行器的高度，r为从地心到飞行器质心的距离，h＝r‑R0，其中R0为地球半径；θ、φ分别为经度和纬度；s为飞行距离；V为相对地球的速度；γ为当地弹道倾角；ψ为弹道偏角；和分别为高度、飞行距离、经度、纬度、弹道倾角、弹道偏角和速度对时间的导数；σ为倾侧角；g＝μ/r2为当地重力加速度，μ为地球重力加速度相关的常数；L和D分别为升力加速度和阻力加速度，其表达式为，L＝ρV2SrefCL/(2m)   (2)D＝ρV2SrefCD/(2m)   (3)式中，m为飞行器质量；Sref为气动参考面积；CL和CD分别为升力系数和阻力系数，与马赫数Ma和攻角α相关；ρ为大气密度，如下所示，ρ＝ρseae‑βh   (4)式中，ρsea为海平面大气密度；β为指数大气模型常数；过程约束如下所示Q.=kρV3.15≤Q.maxa=12ρV2≤qmaxn=L2+D2/g0≤nmax---(5)]]>式中，为热流密度，为最大热流密度，k为常系数；q为动压，qmax为最大动压；n为总过载，nmax为最大总过载，g0＝9.81m/s2；步骤2：滑翔段弹道规划变量设计引入纵向机动加速度aε和横向机动加速度aβ，表达式如下，aε＝Lcosσ+(V2/r‑g)cosγ   (6)aβ＝Lsinσ/cosγ+(V2/r)cosγsinψtanφ   (7)式中，aε和aβ分别为纵向机动加速度和横向机动加速度；将式(6)和式(7)分别带入式(1)，则弹道倾角微分方程和弹道偏角微分方程为，γ.=aϵ/Vψ.=aβ/V---(8)]]>假设aε＜0，则滑翔段的γ单调递减，将运动方程中自变量替换为γ；将式(1)除以式(8)中的弹道倾角微分方程，得dhdγ=V2sinγaϵdsdγ=R0V2cosγraϵ]]>dψdγ=aβaϵdθdγ=V2cosγsinψraϵcosφdφdγ=V2cosγcosψraϵ---(9)]]>dVdγ=-(D+gsinγ)Vaϵ]]>式中，和分别为高度、飞行距离、弹道偏角、经度、纬度和速度的微分方程；其中，第一行为滑翔段纵向运动微分方程，第二行为滑翔段横向运动微分方程，第三行为速度微分方程；设滑翔段的aε与aβ形式如下，aϵ=ρV2Srefkϵ2maβ=ρV2Srefkβ2m---(10)]]>式中，kε和kβ分别为纵向机动加速度比例系数和横向机动加速度比例系数；将式(10)带入式(9)，并取sinγ＝γ、cosγ＝1，则有dhdγ=2mρSrefkϵγ---(11)]]>dsdγ=2mρSrefkϵR0r---(12)]]>dψdγ=kβkϵ---(13)]]>dθdγ=2mρSrefkϵsinψrcosφ---(14)]]>dφdγ=2mρSrefkϵcosψr---(15)]]>以kε和kβ作为规划变量，将滑翔段弹道规划问题转化为求解非线性方程的问题；步骤3：滑翔段弹道解析解求解；设滑翔段的初始弹道倾角为γ0，初始高度为h0，初始弹道偏角为ψ0、初始经度和纬度为分别为θ0和φ0；滑翔段的终点高度为hf，对于给定的kε和kβ，对式(11)至(15)式进行解析求解；第...

【技术特征摘要】
1.一种基于平稳滑行弹道解析解的快速弹道规划方法，包括以下几个步骤：
步骤1：滑翔段弹道规划问题建模
设地球为圆球体，不考虑地球自转，则半速度坐标系下的三自由度质点运动方程如下所
示，
h . = V sin γ θ . = V cos γ sin ψ / ( r cos φ ) ]]> s . = V cos γ R 0 / r φ . = V cos γ cos ψ / r ]]>                                                   (1)
γ . = [ L cos σ + ( V 2 / r - g ) cos γ ] / V ψ . = [ L sin σ / cos γ + ( V 2 / r ) cos γ sin ψ tan φ ] / V ]]> V . = - D - g sin γ ]]>式中，h为飞行器的高度，r为从地心到飞行器质心的距离，h＝r-R0，其中R0为地球
半径；θ、φ分别为经度和纬度；s为飞行距离；V为相对地球的速度；γ为当地弹道倾角；
ψ为弹道偏角；和分别为高度、飞行距离、经度、纬度、弹道倾角、
弹道偏角和速度对时间的导数；σ为倾侧角；g＝μ/r2为当地重力加速度，μ为地球重力加
速度相关的常数；L和D分别为升力加速度和阻力加速度，其表达式为，
L＝ρV2SrefCL/(2m)   (2)
D＝ρV2SrefCD/(2m)   (3)
式中，m为飞行器质量；Sref为气动参考面积；CL和CD分别为升力系数和阻力系数，与
马赫数Ma和攻角α相关；ρ为大气密度，如下所示，
ρ＝ρseae-βh   (4)
式中，ρsea为海平面大气密度；β为指数大气模型常数；
过程约束如下所示
Q . = k ρ V 3.15 ≤ Q . max a = 1 2 ρ V 2 ≤ q max n = L 2 + D 2 / g 0 ≤ n max - - - ( 5 ) ]]>式中，为热流密度，为最大热流密度，k为常系数；q为动压，qmax为最大动压；
n为总过载，nmax为最大总过载，g0＝9.81m/s2；
步骤2：滑翔段弹道规划变量设计
引入纵向机动加速度aε和横向机动加速度aβ，表达式如下，
aε＝Lcosσ+(V2/r-g)cosγ   (6)
aβ＝Lsinσ/cosγ+(V2/r)cosγsinψtanφ   (7)
式中，aε和aβ分别为纵向机动加速度和横向机动加速度；将式(6)和式(7)分别带入
式(1)，则弹道倾角微分方程和弹道偏角微分方程为，
γ . = a ϵ / V ψ . = a β / V - - - ( 8 ) ]]>假设aε＜0，则滑翔段的γ单调递减，将运动方程中自变量替换为γ；将式(1)除以式
(8)中的弹道倾角微分方程，得
dh dγ = V 2 sin γ a ϵ ds dγ = R 0 V 2 cos γ ra ϵ ]]> dψ dγ = a β a ϵ dθ dγ = V 2 cos γ sin ψ ra ϵ cos φ dφ dγ = V 2 cos γ cos ψ ra ϵ - - - ( 9 ) ]]> dV dγ = - ( D + g sin γ ) V a ϵ ]]>式中，和分别为高度、飞行距离、弹道偏角、经度、纬
度和速度的微分方程；其中，第一行为滑翔段纵向运动微分方程，第二行为滑翔段横向运动
微分方程，第三行为速度微分方程；
设滑翔段的aε与aβ形式如下，
a ϵ = ρ V 2 S ref k ϵ 2 m a β = ρ V 2 S ref k β 2 m - - - ( 10 ) ]]>式中，kε和kβ分别为纵向机动加速度比例系数和横向机动加速度比例系数；将式(10)
带入式(9)，并取sinγ＝γ、cosγ＝1，则有
dh dγ = 2 m ρ S ref k ϵ γ - - - ( 11 ) ]]> ds dγ = 2 m ρ S ref k ϵ R 0 r - - - ( 12 ) ]]> dψ dγ = k β k ϵ - - - ( 13 ) ]]> dθ dγ = 2 m ρ S ref k ϵ sin ψ r cos φ - - - ( 14 ) ]]> dφ dγ = 2 m ρ S ref k ϵ cos ψ r - - - ( 15 ) ]]>以kε和kβ作为规划变量，将滑翔段弹道规划问题转化为求解非线性方程的问题；
步骤3：滑翔段弹道解析解求解；
设滑翔段的初始弹道倾角为γ0，初始高度为h0，初始弹道偏角为ψ0、初始经度和纬度为
分别为θ0和φ0；滑翔段的终点高度为hf，对于给定的kε和kβ，对式(11)至(15)式进行
解析求解；
第一、高度解析解
将式(4)带入式(11)可得，
e - βh dh = 2 m ρ sea S ref k ϵ γdγ - - - ( 16 ) ]]>对式(16)积分可得，
γ 2 - γ 0 2 = S ref k ϵ mβ ( ρ 0 - ρ ) - - - ( 17 ) ]]>式中，γ0为初始弹道倾角；ρ0为滑翔段初始大气密度；对式(17)进一步整理可得，
m ρ S ref k ϵ = 1 β ( C 1 - γ 2 ) - - - ( 18 ) ]]>式中，C1为一常数，取值为 C 1 = S ref k ϵ ρ 0 / ( mβ ) + γ 0 2 ; ]]>第二、飞行距离解析解
将式(18)带入式(12)可得，
ds dγ = - 2 R 0 βr ( γ 2 - C 1 ) - - - ( 19 ) ]]>设初始飞行距离s0＝0，则式(19)积分可得，
s＝fs(γ)-fs(γ0)   (20)
式中，fs(γ)为由式(19)不定积分获得的函数，其表达式与C1的正负相关，如下所示，
f s ( x ) = - 2 R 0 βr - C 1 arctan ( x - C 1 ) C 1 < 0 2 R 0 / ( βrx ) C 0 = 0 - R 0 βr C 1 ln ( x - C 1 x + C 1 ) C 1 > 0 ]]>式中，x为自变量；
第三、道偏角解析解
进一步，由式(13)积分可得，
ψ＝ψ0+(kβ/kε)(γ-γ0)   (21)
式中，ψ0为初始弹道偏角；
第四、纬度解析解
将cosψ在π/2附近5阶泰勒展开，可得，
cos ψ = ( π 2 - ψ ) - 1 6 ( π 2 - ψ ) 3 + 1 120 ( π 2 - ψ ) 5 - - - ( 22 ) ]]>取y＝π/2-ψ，其中y为积分中间变量，并将式(18)、式(21)和式(22)带入式(15)
可得，
dφ dy = ( k ϵ rβ k β ) 2 y - y 3 / 3 + y 5 / 60 [ C 2 - ( k ϵ / k β ) y ] 2 - C 1 - - - ( 23 ) ]]>式中，C2为一常量，满足C2＝γ0+kε(π/2-ψ0)/kβ；取，
Ca0＝kβ/(60rβkε)Ca1＝2C2kβ/kε                                          (24)
C a 2 = ( C 2 2 - C 1 ) ( k β / k ϵ ) 2 ]]>式中，Ca0、Ca1和Ca2均为常量；将式(24)带入式(23)可得，
dφ dx = C a 0 ( y 5 - 50 y 3 + 120 y ) y 2 - C a 1 y + C a 2 - - - ( 25 ) ]]>进一步取，
C b 1 = C a 1 2 - C a 2 - 20 ]]> C b 2 = C a 1 3 - 2 C a 1 C a 2 - 20 C a 1 ]]>                                           (26)
C b 3 = C a 1 4 - 3 C a 1 2 C a 2 - 20 C a 1 2 + C a 2 2 + 20 C a 2 + 120 ]]> C b 4 = - C a 1 3 C a 2 + 2 C a 1 C a 2 2 + 20 C a 1 C a 2 ]]>式中，Cb1、Cb2、Cb3和Cb4均为常量；利用式(26)可将式(25)化为，
dφ dy = C a 0 ( y 3 + C a 1 y 2 + C b 1 y + C b 2 ) + ( C a 0 C b 3 / 2 ) ( 2 y - C a 1 ) y 2 - C a 1 y + C a 2 + C a 0 ( C b 3 C a 1 / 2 + C ...

【专利技术属性】
技术研发人员：陈万春，胡锦川，
申请(专利权)人：北京航空航天大学，
类型：发明
国别省市：北京;11