CN104392047B - 一种基于平稳滑翔弹道解析解的快速弹道规划方法 - Google Patents
一种基于平稳滑翔弹道解析解的快速弹道规划方法 Download PDFInfo
- Publication number
- CN104392047B CN104392047B CN201410691412.4A CN201410691412A CN104392047B CN 104392047 B CN104392047 B CN 104392047B CN 201410691412 A CN201410691412 A CN 201410691412A CN 104392047 B CN104392047 B CN 104392047B
- Authority
- CN
- China
- Prior art keywords
- formula
- trajectory
- gamma
- beta
- epsiv
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
- Y02T10/00—Road transport of goods or passengers
- Y02T10/80—Technologies aiming to reduce greenhouse gasses emissions common to all road transportation technologies
- Y02T10/82—Elements for improving aerodynamics
Landscapes
- Navigation (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种基于平稳滑行弹道解析解的快速弹道规划方法,包括以下几个步骤:步骤1:滑翔段弹道规划问题建模;步骤2:滑翔段弹道规划变量设计;步骤3:滑翔段弹道解析解求解;步骤4:滑翔段弹道终端速度控制方案;步骤5:滑翔段弹道再入走廊调整方案;步骤6:滑翔段弹道规划初值生成;步骤7:滑翔段弹道规划流程设计。本发明提出以纵向机动加速度比例系数和横向机动加速度比例系数为滑翔段弹道规划变量,使得运动方程中弹道倾角、弹道偏角、高度、经度与纬度的微分方程中不含速度项;本发明获得了固定纵向机动加速度比例系数和横向机动加速度比例系数对应的滑翔段弹道解析解。
Description
技术领域
本发明涉及一种基于平稳滑行弹道解析解的快速弹道规划方法,属于航天技术、武器技术领域。
背景技术
随着高超声速技术的快速发展,再入弹道规划成为了研究热点。滑翔段是再入弹道的主要构成部分,决定了再入飞行的航程和机动范围。滑翔段弹道快速规划不仅可以分析高超声速飞行器的性能,还可用于在线弹道规划和预测制导,具有较高的研究价值。
由于滑翔段弹道具有飞行时间长,对规划变量高度敏感,且再入过程非线性约束强,使得滑翔弹道的可行域较窄,采用传统优化算法求解时往往需要很长的计算时间。尽管伪谱法、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为飞行器的高度,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为海平面大气密度;β为指数大气模型常数,通常取β=1/7200m-1。
在滑翔段的弹道规划中还需考虑热流密度、动压、过载等过程约束,如下所示
式中,为热流密度,为最大热流密度,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ε<0,则滑翔段的γ单调递减,因此可将运动方程中自变量替换为γ。将式(1)除以式(8)中的弹道倾角微分方程,得
式中,和分别为高度、飞行距离、弹道偏角、经度、纬度和速度对γ的微分。其中,第一行为滑翔段纵向运动微分方程,第二行为滑翔段横向运动微分方程,第三行为速度微分方程。
由式(9)可以看出,若aε和aβ均与V2成正比,则可在纵向运动微分方程和横向运动微分方程中均不显含速度项;同时由于θ和φ的微分方程分别与sinψ和cosψ相关,若要进行解析积分,则ψ的表达式越简单越好,因此aε和aβ之间也成正比。另外,由式(7)可以看出,为了使规划弹道的倾侧角变化幅度较小,aβ还需与ρ成正比。总的来说,可设滑翔段的aε与aβ形式如下,
式中,kε和kβ分别为纵向机动加速度比例系数和横向机动加速度比例系数。将式(10)带入式(9),并取sinγ=γ、cosγ=1(在滑翔段γ为小量),则有
由式(11)至式(12)可知,只要给定kε和kβ,就可对上述微分方程进行解析求解。因此以kε和kβ作为规划变量,则可以将滑翔段弹道规划问题转化为求解非线性方程的问题。
步骤3:滑翔段弹道解析解求解;
设滑翔段的初始弹道倾角为γ0,初始高度为h0,初始弹道偏角为ψ0、初始经度和纬度为分别为θ0和φ0;滑翔段的终点高度为hf,对于给定的kε和kβ,则可对式(11)至(15)式进行解析求解。
第一、高度解析解
将式(4)带入式(11)可得,
对式(16)积分可得,
式中,γ0为初始弹道倾角;ρ0为滑翔段初始大气密度。对式(17)进一步整理可得,
式中,C1为一常数,取值为
第二、飞行距离解析解
将式(18)带入式(12)可得,
设初始飞行距离s0=0,则式(19)积分可得,
s=fs(γ)-fs(γ0) (20)
式中,fs(γ)为由式(19)不定积分获得的函数,其表达式与C1的正负相关,如下所示,
式中,x为自变量。
第三、道偏角解析解
进一步,由式(13)积分可得,
ψ=ψ0+(kβ/kε)(γ-γ0) (21)
式中,ψ0为初始弹道偏角。当滑翔段弹道位于赤道附近时,ψ的取值在π/2附近;若弹道不在赤道附近,则可建立滑翔段起点和终点决定的大圆坐标系,则在大圆坐标系下的弹道偏角ψ还在π/2附近。(如图2所示,坐标原点位于地心,xe1轴由地心指向滑翔起点,ye1位于由地心、滑翔起点、滑翔终点决定的大圆内,且与初始速度方向的夹角小于90deg,ze1由右手定则决定。)
第四、纬度解析解
将cosψ在π/2附近5阶泰勒展开,可得,
取y=π/2-ψ,其中y为积分中间变量,并将式(18)、式(21)和式(22)带入式(15)可得,
式中,C2为一常量,满足C2=γ0+kε(π/2-ψ0)/kβ。取,
式中,Ca0、Ca1和Ca2均为常量。将式(24)带入式(23)可得,
进一步取,
式中,Cb1、Cb2、Cb3和Cb4均为常量。利用式(26)可将式(25)化为,
对式(27)积分可得,
φ=φ0+fφ(y)-fφ(y0) (28)
式中,φ0为初始纬度;y0与初始弹道偏角相关,取值为y0=π/2-ψ0;fφ(y)为与y相关的函数,如下所示
式中,gφ(y)为与y相关的函数,如下所示
第五、经度解析解
式(18)、式(21)和式(29)分别给出了滑翔段弹道的高度、弹道偏角和纬度随弹道倾角变化规律的解析解,将它们分别代入式(14)可得,
可以看出,式(30)的右边仅与γ相关,但复杂度较高,难以解析求解,因此可利用Guass-Legendre求积公式直接获得结果,具体表达式如下,
式中,γi为[γ,γ0]中的第i个高斯节点;Ai为求积系数;n为高斯节点总数。
第六、解析解精度分析
在上面的解析解中,弹道偏角解析解(式21)无任何假设,属于准确解;高度解析解(式18)的误差来源于sinγ一阶泰勒展开,通常γ<2deg,因此高度解析解的相对误差大约为10-4;飞行距离解析解(式20)并未引入新的误差项,因此相对误差与高度相同;纬度解析解(式29)中引入了cosψ的五阶泰勒展开,若ψ的变化范围为[0,π],则展开误差最大为0.0047,因此纬度解析解的相对误差介于0.1%至1%;而求解经度解析解时使用了纬度解析解,并且还存在求积误差,因此其相对误差最大,在1%左右。
步骤4:滑翔段弹道终端速度控制方案;
在步骤3中,给出了给定kε和kβ条件下的滑翔段弹道解析解,但没有给出速度的解析解。为了获得终端速度的大小,需要利用上述解析解对速度进行数值积分。
设负比能量e=μ/r-V2/2,则由式(1)可得,
式中,为负比能量的导数。由式(8)、式(10)和式(32)可得,
由式(33)可以看出,CD是速度求解的关键。通过情况下,CD可写成CL和e的函数,如下所示,
CD=fCD(CL,e) (34)
式中,fCD为升力系数与阻力系数的关系函数。求解式(33)积分的关键为求解CL,由式(6)、(7)式和式(10)可得,
另外,由式(4)和式(18)可得,
由(35)、式(36)和式(37)可求得CL,将CL带入式(34)可得CD的大小。
最终,将r和CD的表达式带入式(33),并采用龙格-库塔法进行积分,可获得终端负比能量,从而解得终端的速度大小。需要指出的是,由式(33)易知CD越大则终端速度越小,而CL越大则CD越大,因此调整CL可以控制终端速度。由于kε通常为一小量,还决定了射程的大小,因此不能用于速度调节;因而只能采用kβ来调整终端速度。图3给出了与Vf的关系,可以看出两者之间的线性化程度非常高,有利于弹道规划时的数值求解。
当规划的弹道包含倾侧反转点时,则反转点前后的kβ需要满足如下连接条件,
式中,kβ-和kβ+分别为反转点前后的横向机动加速度比例系数;ρc、ψc和φc分别为γc处的大气密度、弹道偏角和纬度。
步骤5:滑翔段弹道再入走廊调整方案;
采用初始弹道倾角调整再入弹道在走廊中的位置。由图4可知,初始弹道倾角γ0越大,则再入弹道的高度越高,距离由最大热流密度约束、最大动压约束和最大过载约束确定的再入走廊下边界越远。但考虑到滑翔段攻角的平滑性,初始弹道倾角取值如下所示,
γ0=γ*+Δγ (39)
式中,Δγ为调整弹道再入走廊位置的初始弹道倾角分量;γ*为满足攻角平滑性要求的初始弹道倾角分量,满足,
γ*=-2g/(KNβV0 2) (40)
式中KN=Lcosσ/D,为纵向升阻比。需要指出的是,由于本发明规划的滑翔段弹道满足平稳滑翔条件,因而热流密度、动压和过载变化较为平缓,峰值也较小,因此Δγ通常为零。
步骤6:滑翔段弹道规划初值生成;
在进行滑翔段弹道规划前,需要对sf、γ*、kε和kβ的初值进行估算。
第一、sf估算
球面坐标系下,(θ0,φ0)和(θf,φf)的距离为(如图5所示),
式中,ss为(θ0,φ0)至(θf,φf)的球面距离;和分别满足,
考虑到初始弹道偏角偏离射面的情况,终端飞行距离为,
式中,其中为滑翔段起点子午线与射面所在大圆的夹角。
第二、γ*估算
在滑翔段,纵向升阻比与终端飞行距离存在如下关系,
将式(42)带入式(40),可以求得γ*的估值为,
第三、kε估算
设Δγ=0,则γ0=γ*。由式(17)和式(20)可知,γf、kε和sf满足如下关系,
式中,γf和ρf分别为终端弹道倾角和终端大气密度。由式(44)可知,给定sf,则存在唯一的kε与之对应(如图6所示)。将式(41)得到的sf带入式(44)求解即可获得kε的估值。
第四、kβ估算
由图2可知,在kβ较小时,与Vf的线性关系较好。因此在估算kβ初值时的原则为宁小勿大,这样有利于之后进行迭代修正。由式(21)知,在已知(γ-γ0)/kε的条件下,kβ仅与ψf-ψ0相关。假设ψf-ψ0仅用于修正ψ10的影响,则根据平面圆弧弹道假设有,
ψf-ψ0=2ψ10
从而可得kβ的初始估值为,
kβ=2ψ10CN1/(γ-γ0) (45)
步骤7:滑翔段弹道规划流程设计;
滑翔段弹道规划任务为寻找合适的kε、kβ、Δγ和倾侧角反转点γc,使得再入飞行器从滑翔起点(θ0,φ0,r0,V0,ψ0)至滑翔终点(θf,φf,rf,Vf)。上述四个规划变量分别对应着不同的约束要求,彼此之间几乎不耦合,其中:kε用于调整滑翔段的飞行距离;kβ用于调整滑翔段终点的速度;γc用于调整滑翔段的终点位置;Δγ将用于调整再入弹道高度,以满足过程约束要求。利用滑翔段弹道形状解析解,可将滑翔段弹道规划分拆成弹道形状规划和终端速度校正两部分。其中,弹道形状规划不需要进行弹道积分,解决终端位置约束问题;终端端速度校正则需要对速度微分方程进行积分,解决终端速度约束和过程约束问题。完整的滑翔段弹道流程如图7所示,具体说明如下:
第一、根据步骤6,获取sf、γ*、kε和kβ的规划初值,并假设Δγ=0和γc=γ0。
第二、根据给定的Δγ、γ*、γc、kε和kβ,利用步骤3计算θ(γf)和φ(γf)。其中,γc为反转点弹道倾角;θ(γf)和φ(γf)为规划弹道终端经度和纬度。
第三、修正γc,使得终端纬度满足要求,修正的方法为牛顿迭代法,迭代的终止条件为|φ(γf)-φf|<φlimit。式中φlimit为纬度规划精度,考虑步骤3中的解析弹道计算精度,可取φlimit=0.01deg。若|φ(γf)-φf|≥φlimit则转向步骤7中的第二步。
第四、修正sf(或者kε)使得终端经度满足要求,修正方法为牛顿迭代法,迭代终止条件为|θ(γf)-θf|<θlimit。式中θlimit为经度规划精度,考虑到步骤3中的解析弹道计算精度,可取θlimit=0.1deg。若|θ(γf)-θf|≥θlimit则转向步骤7中的第二步。
第五、根据Δγ、γ*、γc、kε和kβ,利用步骤4积分得到滑翔段弹道终端速度V(γf),并采用牛顿迭代法修正kβ,使得|V(γf)-Vf|<Vlimit。式中Vlimit为速度规划精度,考虑到步骤3解析弹道计算精度,可取Vlimit=10m/s。若|V(γf)-Vf|≥Vlimit则转向步骤7中的第二步。
第六,根据上一步的积分弹道判断是否满足过程约束,并对Δγ进行如下修正,
式中,h为实际飞行高度;hlimit为由最大热流密度约束、最大动压约束和最大过载约束确定的再入走廊高度下边界;min(h-hlimit)<0表示滑翔段弹道超出再入走廊,需要增大Δγ,并转向步骤7中的第二步;min(h-hlimit)≥0表示滑翔段弹道满足再入走廊要求;γlim为min(h-hlimit)取极值时的弹道倾角;Δγ(old)和Δγ(new)分别为修正前和修正后的初始弹道倾角修正量;h(γlim)、s(γlim)和hlimit(γlim)分别为弹道倾角取γlim时的高度、飞行距离和高度下边界。
实施案例:
为了检验解析求解算法的精度,选用CAV作为计算模型,进行数值仿真效验。仿真采用无转动圆球模型和指数大气模型,参数设置如表格1。仿真计算机CPU为Core i3-2120,内存4GB,仿真环境为Matlab。
表格1仿真参数设置
第一,解析解精度效验
在进行弹道规划前,首先对弹道形状解析解精度进行验证,并与Bell解析解进行对比。图8和图9分别给出了纵向弹道解析解与横向弹道解析解与弹道积分结果的对比,可以看出,弹道形状解析解与弹道积分结果几乎完全重合,具有极高的精度。具体来说,采用解析解直接预测滑翔段终点位置,其高度误差小于1m;弹道偏角误差小于0.01deg;大圆坐标系下的纬度误差小于0.1deg,经度误差小于0.3deg;验证了文中的误差分析。图10给出了Bell给出的横向弹道解析解,可以看出Bell解析解只在较小的范围内具有较高的精度,当飞行距离较远时,横程偏差较大。这是由于在Bell解析解中忽略了地球曲率对横向弹道的影响,而采用平面坐标系进行求解,带来了较大的求解误差。
第二,弹道规划可达区域分析
在验证了弹道形状解析解的精度之后,还需分析弹道规划的可达区域。在求解可达区域时,假定弹道偏角为固定值,因而主要通过倾侧反转来调节终端位置,因此无反转的弹道即为弹道规划可达区域的边界。图11至图18给出了边界弹道曲线族。由图11和图12可知,规划弹道的攻角曲线和倾侧角曲线非常光滑且变化幅度较小;由图13可知滑翔段弹道规划的可达区域纵程覆盖范围约为11000千米,横程覆盖范围约为6700千米;由图14和图15可知规划弹道的终端高度和终端速度满足要求;由图16至图18可知规划弹道满足过程约束要求。对于目标在可达区域内的情况,由于仅仅增加了倾侧角反转点来调节终端经纬度,因此其攻角曲线、倾侧曲线、纵向弹道曲线和过程约束曲线的趋势将与图11至图18中给出的结果类似。
在图13的基础上,图19进一步获得了滑翔段规划弹道的可达区域。与平稳滑翔最大可达区域相比,本发明获得的可达区域的要小一些,最大射程损失约为81.7km,最小射程损失约为302.8km,最大横程损失约为923.9km。这是由于本发明选用固定的kβ作为横向机动规划变量,使得在大机动情况下倾侧角存在较大的变化(如图11所示),这会增加横程机动的能量损耗,从而降低横向机动能力。当射程较小或者横程较大时,需要较大的横向机动,因而出现图18所示的可达区域损失。
第三,滑翔段弹道快速规划
弹道规划的初值和终值如前所示,采用图5给出的滑翔段弹道规划流程,对图19给出的可达区域内的目标点进行弹道规划,得到的结果如表格1和图20至图22所示。由表格2可以看出,基于反馈线性化解析解的滑翔段弹道规划方法能够在0.03s左右的时间内规划出一条满足约束要求的再入弹道,并且弹道终端具有较高的位置精度。不足在于规划弹道的终端速度误差较大,这是由于终端速度是根据几何弹道、采用动态逆方法积分获得的,使得弹道形状上的微小偏差通过求解速度的积分得以放大。图20给出了规划弹道的形状;图21给出出了规划弹道的攻角曲线,可以看出攻角变化较为光滑,变化幅度小于4deg;图22给出了倾侧角曲线,同样也较为光滑。
表格2滑翔段弹道规划结果
Claims (1)
1.一种基于平稳滑行弹道解析解的快速弹道规划方法,包括以下几个步骤:
步骤1:滑翔段弹道规划问题建模
设地球为圆球体,不考虑地球自转,则半速度坐标系下的三自由度质点运动方程如下所示,
式中,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为海平面大气密度;β为指数大气模型常数;
过程约束如下所示
式中,为热流密度,为最大热流密度,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ε<0,则滑翔段的γ单调递减,将运动方程中自变量替换为γ;将式(1)除以式(8)中的弹道倾角微分方程,得
式中,和分别为高度、飞行距离、弹道偏角、经度、纬度和速度对γ的微分;其中,第一行为滑翔段纵向运动微分方程,第二行为滑翔段横向运动微分方程,第三行为速度微分方程;
设滑翔段的aε与aβ形式如下,
式中,kε和kβ分别为纵向机动加速度比例系数和横向机动加速度比例系数;将式(10)带入式(9),并取sinγ=γ、cosγ=1,则有
以kε和kβ作为规划变量,将滑翔段弹道规划问题转化为求解非线性方程的问题;
步骤3:滑翔段弹道解析解求解;
设滑翔段的初始弹道倾角为γ0,初始高度为h0,初始弹道偏角为ψ0、初始经度和纬度为分别为θ0和φ0;滑翔段的终点高度为hf,对于给定的kε和kβ,对式(11)至(15)式进行解析求解;
第一、高度解析解
将式(4)带入式(11)可得,
对式(16)积分可得,
式中,γ0为初始弹道倾角;ρ0为滑翔段初始大气密度;对式(17)进一步整理可得,
式中,C1为一常数,取值为
第二、飞行距离解析解
将式(18)带入式(12)可得,
设初始飞行距离s0=0,则式(19)积分可得,
s=fs(γ)-fs(γ0) (20)
式中,fs(γ)为由式(19)不定积分获得的函数,其表达式与C1的正负相关,如下所示,
式中,x为自变量;
第三、道偏角解析解
进一步,由式(13)积分可得,
ψ=ψ0+(kβ/kε)(γ-γ0) (21)
式中,ψ0为初始弹道偏角;
第四、纬度解析解
将cosψ在π/2附近5阶泰勒展开,可得,
取y=π/2-ψ,其中y为积分中间变量,并将式(18)、式(21)和式(22)带入式(15)可得,
式中,C2为一常量,满足C2=γ0+kε(π/2-ψ0)/kβ;取,
式中,Ca0、Ca1和Ca2均为常量;将式(24)带入式(23)可得,
进一步取,
式中,Cb1、Cb2、Cb3和Cb4均为常量;利用式(26)可将式(25)化为,
对式(27)积分可得,
φ=φ0+fφ(y)-fφ(y0) (28)
式中,φ0为初始纬度;y0与初始弹道偏角相关,取值为y0=π/2-ψ0;fφ(y)为与y相关的函数,如下所示
式中,gφ(y)为与y相关的函数,如下所示
第五、经度解析解
式(18)、式(21)和式(29)分别给出了滑翔段弹道的高度、弹道偏角和纬度随弹道倾角变化规律的解析解,将它们分别代入式(14)可得,
利用Guass-Legendre求积公式直接获得结果,具体表达式如下,
式中,γi为[γ,γ0]中的第i个高斯节点;Ai为求积系数;n为高斯节点总数;
步骤4:滑翔段弹道终端速度控制方案;
设负比能量e=μ/r-V2/2,则由式(1)可得,
式中,为负比能量的导数;由式(8)、式(10)和式(32)可得,
由式(33)看出,CD是速度求解的关键;通过情况下,CD可写成CL和e的函数,如下所示,
CD=fCD(CL,e) (34)
式中,fCD为升力系数与阻力系数的关系函数;求解式(33)积分的关键为求解CL,由式(6)、(7)式和式(10)可得,
另外,由式(4)和式(18)可得,
由(35)、式(36)和式(37)可求得CL,将CL带入式(34)可得CD的大小;
最终,将r和CD的表达式带入式(33),并采用龙格-库塔法进行积分,获得终端负比能量,从而解得终端的速度大小;
当规划的弹道包含倾侧反转点时,则反转点前后的kβ需要满足如下连接条件,
式中,kβ-和kβ+分别为反转点前后的横向机动加速度比例系数;ρc、ψc和φc分别为γc处的大气密度、弹道偏角和纬度;
步骤5:滑翔段弹道再入走廊调整方案;
初始弹道倾角取值如下所示,
γ0=γ*+Δγ (39)
式中,Δγ为调整弹道再入走廊位置的初始弹道倾角分量;γ*为满足攻角平滑性要求的初始弹道倾角分量,满足,
式中KN=L cosσ/D,为纵向升阻比;
步骤6:滑翔段弹道规划初值生成;
在进行滑翔段弹道规划前,需要对sf、γ*、kε和kβ的初值进行估算;
第一、sf估算
球面坐标系下,(θ0,φ0)和(θf,φf)的距离为,
式中,ss为(θ0,φ0)至(θf,φf)的球面距离;和分别满足,
则终端飞行距离为,
式中,其中为滑翔段起点子午线与射面所在大圆的夹角;
第二、γ*估算
在滑翔段,纵向升阻比与终端飞行距离存在如下关系,
将式(42)带入式(40),求得γ*的估值为,
第三、kε估算
设Δγ=0,则γ0=γ*;由式(17)和式(20)可知,γf、kε和sf满足如下关系,
式中,γf和ρf分别为终端弹道倾角和终端大气密度;由式(44)可知,给定sf,则存在唯一的kε与之对应;将式(41)得到的sf带入式(44)求解即可获得kε的估值;
第四、kβ估算
假设ψf-ψ0仅用于修正ψ10的影响,则根据平面圆弧弹道假设有,
ψf-ψ0=2ψ10
kβ的初始估值为,
kβ=2ψ10CN1/(γ-γ0) (45)
步骤7:滑翔段弹道规划流程设计;
具体为:
第一、根据步骤6,获取sf、γ*、kε和kβ的规划初值,并假设Δγ=0和γc=γ0;
第二、根据给定的Δγ、γ*、γc、kε和kβ,利用步骤3计算θ(γf)和φ(γf);其中,γc为反转点弹道倾角;θ(γf)和φ(γf)为规划弹道终端经度和纬度;
第三、修正γc,使得终端纬度满足要求,修正的方法为牛顿迭代法,迭代的终止条件为|φ(γf)-φf|<φlimit;式中φlimit为纬度规划精度;若|φ(γf)-φf|≥φlimit则转向步骤7中的第二步;
第四、修正sf或者kε使得终端经度满足要求,修正方法为牛顿迭代法,迭代终止条件为|θ(γf)-θf|<θlimit;式中θlimit为经度规划精度;若|θ(γf)-θf|≥θlimit则转向步骤7中的第二步;
第五、根据Δγ、γ*、γc、kε和kβ,利用步骤4积分得到滑翔段弹道终端速度V(γf),并采用牛顿迭代法修正kβ,使得式中为速度规划精度;若|V(γf)-Vf|≥Vlimit则转向步骤7中的第二步;
第六,根据上一步的积分弹道判断是否满足过程约束,并对Δγ进行如下修正,
式中,h为实际飞行高度;hlimit为由最大热流密度约束、最大动压约束和最大过载约束确定的再入走廊高度下边界;min(h-hlimit)<0表示滑翔段弹道超出再入走廊,需要增大Δγ,并转向步骤7中的第二步;min(h-hlimit)≥0表示滑翔段弹道满足再入走廊要求;γlim为min(h-hlimit)取极值时的弹道倾角;Δγ(old)和Δγ(new)分别为修正前和修正后的初始弹道倾角修正量;h(γlim)、s(γlim)和hlimit(γlim)分别为弹道倾角取γlim时的高度、飞行距离和高度下边界。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410691412.4A CN104392047B (zh) | 2014-11-25 | 2014-11-25 | 一种基于平稳滑翔弹道解析解的快速弹道规划方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410691412.4A CN104392047B (zh) | 2014-11-25 | 2014-11-25 | 一种基于平稳滑翔弹道解析解的快速弹道规划方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104392047A CN104392047A (zh) | 2015-03-04 |
CN104392047B true CN104392047B (zh) | 2017-05-10 |
Family
ID=52609950
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201410691412.4A Active CN104392047B (zh) | 2014-11-25 | 2014-11-25 | 一种基于平稳滑翔弹道解析解的快速弹道规划方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104392047B (zh) |
Families Citing this family (19)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104809271B (zh) * | 2015-03-23 | 2016-04-20 | 北京航天自动控制研究所 | 一种升力式飞行器的再入轨迹的计算方法 |
CN104732106B (zh) * | 2015-04-08 | 2017-11-28 | 中国人民解放军国防科学技术大学 | 考虑不确定因素影响的飞行走廊计算方法 |
CN105022858B (zh) * | 2015-05-08 | 2016-06-01 | 北京航天自动控制研究所 | 一种确定滑翔飞行器阻力加速度走廊边界的方法 |
CN105589068B (zh) * | 2015-12-08 | 2017-09-22 | 河海大学 | 基于三步数值积分的弹道外推方法 |
CN105674804B (zh) * | 2015-12-25 | 2017-06-06 | 北京航空航天大学 | 一种包含法向加速度导数的空射巡航弹下滑段多约束制导方法 |
CN105676638B (zh) * | 2016-01-11 | 2018-06-29 | 北京航空航天大学 | 平稳滑翔/准自然频率跳跃滑翔组合机动突防弹道规划方法 |
CN105718660B (zh) * | 2016-01-21 | 2019-03-01 | 中国工程物理研究院总体工程研究所 | 临近空间大范围机动弹道三维包络计算方法 |
CN106227972A (zh) * | 2016-08-04 | 2016-12-14 | 北京航空航天大学 | 一种高超声速飞行器平稳滑翔弹道的优化方法 |
CN106446466B (zh) * | 2016-11-09 | 2019-07-16 | 沈阳航空航天大学 | 基于可编辑配置参数接口的四旋翼快速建模设计方法 |
CN106643341B (zh) * | 2017-02-24 | 2018-06-01 | 北京临近空间飞行器***工程研究所 | 基于准平衡滑翔原理的力热控制耦合设计方法 |
CN107621198B (zh) * | 2017-08-28 | 2019-04-12 | 北京航空航天大学 | 一种绕多禁飞区的多步决策弹道规划方法 |
CN109190248B (zh) * | 2018-09-03 | 2023-07-18 | 中国运载火箭技术研究院 | 一种用于滑翔飞行器的滑翔射程解析方法及解析*** |
CN110276131B (zh) * | 2019-06-24 | 2022-07-26 | 西北工业大学 | 基于多项式响应面模型的翼身融合水下滑翔机外形优化方法 |
CN110309590B (zh) * | 2019-06-28 | 2021-01-19 | 北京理工大学 | 一种再入飞行器速度-高度再入走廊预测方法 |
CN111348223B (zh) * | 2020-05-25 | 2020-08-21 | 北京星际荣耀空间科技有限公司 | 一种控制弹道顶点高度的闭路制导方法、装置及设备 |
CN111859526B (zh) * | 2020-06-04 | 2024-06-04 | 中国人民解放军国防科技大学 | 助推滑翔导弹总体参数的快速确定方法 |
CN111859527B (zh) * | 2020-06-04 | 2022-08-23 | 中国人民解放军国防科技大学 | 一种助推滑翔导弹全程弹道在线规划方法 |
CN116880527B (zh) * | 2023-07-20 | 2024-02-23 | 中国空气动力研究与发展中心空天技术研究所 | 高超声速飞行器跳跃滑翔飞行射程最大的控制方法及*** |
CN118278224A (zh) * | 2024-06-03 | 2024-07-02 | 西安现代控制技术研究所 | 一种基于能量降阶预测模型的火箭弹滑翔弹道重规划方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103838914A (zh) * | 2013-12-30 | 2014-06-04 | 北京航空航天大学 | 一种高超声速飞行器滑翔段弹道解析求解方法 |
CN103983143A (zh) * | 2014-04-04 | 2014-08-13 | 北京航空航天大学 | 包含速度过程约束和多终端约束的空地导弹投放下滑段制导方法 |
CN104035335A (zh) * | 2014-05-27 | 2014-09-10 | 北京航空航天大学 | 基于高精度纵、横程解析预测方法的平稳滑翔再入制导律 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6877691B2 (en) * | 2002-03-12 | 2005-04-12 | Bae Systems Information And Electronic Systems Integration Inc. | High altitude stripping for threat discrimination |
-
2014
- 2014-11-25 CN CN201410691412.4A patent/CN104392047B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103838914A (zh) * | 2013-12-30 | 2014-06-04 | 北京航空航天大学 | 一种高超声速飞行器滑翔段弹道解析求解方法 |
CN103983143A (zh) * | 2014-04-04 | 2014-08-13 | 北京航空航天大学 | 包含速度过程约束和多终端约束的空地导弹投放下滑段制导方法 |
CN104035335A (zh) * | 2014-05-27 | 2014-09-10 | 北京航空航天大学 | 基于高精度纵、横程解析预测方法的平稳滑翔再入制导律 |
Non-Patent Citations (1)
Title |
---|
基于Radau伪谱法的制导炸弹最优滑翔弹道研究;袁宴波等;《兵工学报》;20140815;第35卷(第8期);第1179-1186页 * |
Also Published As
Publication number | Publication date |
---|---|
CN104392047A (zh) | 2015-03-04 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104392047B (zh) | 一种基于平稳滑翔弹道解析解的快速弹道规划方法 | |
US11286065B2 (en) | Method for designing reentry trajectory based on flight path angle planning | |
US11079239B2 (en) | Method for directly planning reentry trajectory in height-velocity profile | |
Lugo-Cárdenas et al. | Dubins path generation for a fixed wing UAV | |
Karimi et al. | Optimal maneuver-based motion planning over terrain and threats using a dynamic hybrid PSO algorithm | |
Mease et al. | Reduced-order entry trajectory planning for acceleration guidance | |
Lu et al. | Rapid generation of accurate entry landing footprints | |
CN105222780B (zh) | 一种基于Stirling插值多项式逼近的椭球集员滤波方法 | |
Zhang et al. | Entry trajectory planning based on three-dimensional acceleration profile guidance | |
CN109062241B (zh) | 基于线性伪谱模型预测控制的自主全射向再入制导方法 | |
CN104035335A (zh) | 基于高精度纵、横程解析预测方法的平稳滑翔再入制导律 | |
CN106586033A (zh) | 自适应分段的多段线性伪谱广义标控脱靶量再入制导方法 | |
CN105929842A (zh) | 一种基于动态速度调节的欠驱动uuv平面轨迹跟踪控制方法 | |
Yang et al. | Autonomous entry guidance using linear pseudospectral model predictive control | |
CN110750850B (zh) | 强约束复杂任务条件下的三维剖面优化设计方法、***及介质 | |
CN108267953A (zh) | 一种基于领航者-跟踪者水下机器人位置跟踪方法 | |
CN105843224A (zh) | 基于神经动态模型反步法的auv水平面路径跟踪控制方法 | |
CN104597911A (zh) | 空中加油受油机自适应最优对接轨迹跟踪飞行控制方法 | |
Hu et al. | Analytical solution for nonlinear three-dimensional guidance with impact angle and field-of-view constraints | |
Leavitt et al. | Performance of evolved acceleration guidance logic for entry (EAGLE) | |
Gao et al. | Dubins path‐based dynamic soaring trajectory planning and tracking control in a gradient wind field | |
CN105300383A (zh) | 一种基于回溯搜索的无人机空中加油位姿估计方法 | |
CN105718660A (zh) | 临近空间大范围机动弹道三维包络计算方法 | |
CN104809271B (zh) | 一种升力式飞行器的再入轨迹的计算方法 | |
CN104614993B (zh) | 微陀螺仪自适应滑模预设性能控制方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |