CN105740506A - 沿临近空间大范围机动弹道空间包络的扰动引力逼近方法 - Google Patents
沿临近空间大范围机动弹道空间包络的扰动引力逼近方法 Download PDFInfo
- Publication number
- CN105740506A CN105740506A CN201610042316.6A CN201610042316A CN105740506A CN 105740506 A CN105740506 A CN 105740506A CN 201610042316 A CN201610042316 A CN 201610042316A CN 105740506 A CN105740506 A CN 105740506A
- Authority
- CN
- China
- Prior art keywords
- phi
- cos
- lambda
- sin
- trajectory
- 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.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Navigation (AREA)
Abstract
本发明公开了一种沿临近空间大范围机动弹道空间包络的扰动引力逼近方法,包括以下步骤:建立弹道包络侧向宽度与纵程的数学模型;建立换极坐标系;建立换极坐标系中飞行器动力学模型;计算换极坐标系弹道三维包络;换极坐标系全域空域剖分;建立一般坐标系空间包络扰动引力重构模型;当前弹道位置所在网格判断;建立计算网格内局域坐标系;网格内部扰动引力网函数逼近计算。本发明涉及的弹道空间包络计算方法可覆盖任意情况下的临近空间大范围机动弹道,根据本发明方法进行沿弹道的扰动引力计算,具有计算规模小、计算速度快、计算精度高的特征,能满足弹道规划、制导计算等对计算量和计算速度的要求。
Description
技术领域
本发明涉及一种飞行器动力学建模领域的沿临近空间大范围机动弹道空间包络的扰动引力快速计算方法,尤其涉及一种沿临近空间大范围机动弹道空间包络的扰动引力逼近方法。
背景技术
临近空间大范围机动弹道特指高超声速滑翔飞行器包含初始下降段弹道和滑翔段弹道在内的弹道,其起点为再入点,终点为滑翔段弹道终点。该类弹道长时间处于临近空间,具有侧向大范围机动的特性,受地球外部空间扰动引力影响显著。扰动引力通过影响导航解算而影响飞行器当前状态量的确定,进而对制导计算和弹道精度产生影响。有计算表明,由扰动引力引起的临近空间大范围机动弹道终端位置偏差可达几十公里,飞行器可能因无法进入末制导有效半径而导致飞行任务失败,因此有必要在导航解算和制导计算中考虑扰动引力的影响。
对于采用纯惯性导航方案的飞行器,弹载设备无法测量扰动引力,引力项精度只能通过射前构建扰动引力模型并进行飞行过程中实时计算来保证。射前模型的重构速度受飞行器发射快速性指标的约束,实时计算方法的计算量、计算速度和存储量等受弹载计算机数据存储能力和计算能力的约束。现有地球物理和大地测量领域的传统计算方法,通常具有数据存储量大、计算规模大、计算速度慢的特征,难以满足射前引力模型快速构建及飞行过程中快速计算的要求。与沿弹道式弹道的扰动引力快速计算问题相比,沿临近空间大范围机动弹道的扰动引力快速计算具有更加严格的要求。首先,弹道侧向机动范围可达几千公里,弹道形态因飞行任务不同而迥异,因此建立的方法须能适应射前飞行任务临时变更的情况。其次,由于弹道具有大范围侧向机动的特性,需合理确定逼近空域范围以防止实际弹道超出逼近区间。最后,由于弹道处于扰动引力变化复杂的临近空间,对逼近方法的精度具有更高要求。
目前尚没有一种成熟的方法,能够达到沿任意临近空间大范围机动弹道空间包络的扰动引力快速计算的要求。
发明内容
本发明的目的就在于为了解决上述问题而提供一种逼近精度高、计算速度快、数据存储量小的沿临近空间大范围机动弹道空间包络的扰动引力逼近方法。
本发明通过以下技术方案来实现上述目的:
一种沿临近空间大范围机动弹道空间包络的扰动引力逼近方法,包括以下步骤:
(1)建立弹道包络侧向宽度与纵程的数学模型;
(2)建立换极坐标系;
(3)建立换极坐标系中飞行器动力学模型;
(4)计算换极坐标系弹道三维包络;
(5)换极坐标系全域空域剖分;
(6)建立一般坐标系空间包络扰动引力重构模型;
(7)当前弹道位置所在网格判断;
(8)建立计算网格内局域坐标系;
(9)网格内部扰动引力网函数逼近计算。
作为优选,所述步骤(1)的方法为:
令弹道再入点为I,其经度为λI、地心纬度为φI,弹道落点为T,其经度为λT、地心纬度为φT;以由再入点和落点确定的再入大圆弧平面为对称面,称沿射向方向对称面以左弹道为左侧机动弹道,沿射向方向对称面以右弹道为右侧机动弹道;记最大左侧机动弹道距对称面的最大距离为左边界Bl,记最大右侧机动弹道距对称面的最大距离为右边界Br;
根据给定的飞行器再入点飞行状态参数、弹道终端飞行状态参数、飞行过程约束条件及终端约束条件,针对再入点为(λI,φI)=(0°,0°)、落点为(λki,0°)(i=1,2,3,4,…)的低空大范围机动弹道进行弹道计算,记各落点对应的弹道纵程为Li,由下式(1)计算Li:
Li=Rearccos(sinφIsinφT+cosφIcosφTcos(λT-λI))(1)
其中,Re为地球半径;
通过改变再入点初始速度方位角获得最大左侧机动弹道和最大右侧向机动弹道,记各落点对应的弹道侧向包络宽度为Wio,由下式(2)计算Wio:
Wio=Bl-Br(2)
取Cw倍Wio为实际侧向包络宽度Wi,对Li(i=1,2,3,4,…)和Wi(i=1,2,3,4,…)进行拟合,可建立如下式(3)所示的弹道包络侧向宽度与纵程的数学模型,确定模型系数为ai:
所述步骤(2)的方法为:
引入一个换极坐标系:用表示换极坐标系中各物理量,用X表示一般坐标系中各物理量,按如下方式建立换极坐标系:
首先,定义一个再入大圆弧平面作为换极赤道平面:对目标点确定的情况,将再入点和目标点地心矢径构成的再入大圆弧平面作为换极赤道平面;对于目标点未确定的情况,根据再入点位置及速度方位角确定的再入大圆弧平面作为换极赤道面;
然后,基于换极赤道平面定义换极坐标系OE为地心,轴沿再入点地心矢径方向,轴在换极赤道面内垂直于轴指向目标点方向,轴与轴、轴构成右手系;
所述步骤(3)的方法为:
在换极坐标系中建立以时间为自变量的滑翔飞行器动力学方程,其飞行状态量为换极后的经度地心纬度航迹偏航角速度速度倾角和地心距
其中,Cσ、Cθ为哥氏加速度项,和为牵连加速度项,
其中,
其中,ωe为地球旋转加速度矢量,λp和φp为换极后极点P的经度和地心纬度,AP为P的方位角;
根据换极坐标系定义,一般坐标系与换极坐标系中地心距、当地速度倾角及速度的定义一致,
定义
其中,ψf为点F的方位角,
由一般坐标系中λ和φ确定换极坐标系中和的表达式为:
由换极坐标系中λ和φ确定一般坐标系中λ和φ的表达式为:
由一般坐标系中σ确定换极坐标系中的表达式为:
其中,
所述步骤(4)的方法为:
根据步骤(3)所述坐标转换关系,根据弹道再入点经度λI、再入点地心纬度φI、落点经度λT、落点地心纬度φT计算得到换极弹道参数,由计算得 记换极弹道落点经度为
由式(15)计算换极弹道纵程
其中,Re为地球半径;
将代入式(16),可计算得到换极弹道侧向包络宽度
其中,模型系数ai由步骤(1)确定;
将换极弹道侧向包络描述为长度为宽度为的矩形,其边界由式(17)确定:
其中,为换极系侧向包络东向下边界,为换极系侧向包络东向上边界;为换极系侧向包络北向下边界,为换极系侧向包络北向上边界;
侧向包络对应的经度范围Δλ及地心纬度范围Δφ由式(18)确定:
以弹道高度范围的Cr倍作为弹道包络的垂向范围
所述步骤(5)的方法为:
记换极系空间包络确定的空域为Ω,按式(19)计算Ω边界:
其中,为换极系Ω天向下边界,为换极系Ω天向上边界,为换极系Ω东向下边界,为换极系Ω东向上边界,为换极系Ω北向下边界,为换极系Ω北向上边界,dr、dλ和dφ分别为空域剖分间隔;
按照天向东向北向的间隔将由式(19)确定的空间包络空域Ω均匀剖分为q个互不重叠的子域Ωe(e=1,2,…,q),记网格节点坐标为
所述步骤(6)的方法为:
根据换极系网格节点坐标由步骤(3)所述坐标转换关系计算一般坐标系网格节点坐标
由球谐函数方法计算一般坐标系网格节点扰动引力位TG,
其中,
与TG相应的引力即为扰动引力加速度δg,即
δg=gradTG(23)
则扰动引力加速度δg在天东北坐标系OE-REN中的三个分量δgR、δgE、δgN为:
存储一般坐标系节点位置三分量和节点扰动引力三分量,完成空间包络扰动引力重构模型构建;
所述步骤(7)的方法为:
令当前弹道位置坐标为A(r,λ,φ),当满足如式(25)所示的条件时,
判断A(r,λ,φ)位于由节点N(rgi,λgj,φgk)、N(rgi,λgj+1,φgk)、N(rgi,λgj,φgk+1)、N(rgi,λgj+1,φgk+1)、N(rgi+1,λgj,φgk)、N(rgi+1,λgj+1,φgk)、N(rgi+1,λgj,φgk+1)、N(rgi+1,λgj+1,φgk+1)所确定的网格中;
所述步骤(8)的方法为:
局部坐标系由半径rl=rgi+dr/2的球面、经度λl=λgj+dλ/2的子午面、纬度φl=φgk+dφ/2的纬圈的交线组成,原点l(rl,λl,φl)为三交线的交点,原点局部坐标为l(0,0,0),ξ、η、ζ分别沿原点l的天向、北向和东向,则单元内任意点A(r,λ,φ)的局部坐标A(ξ,η,ζ)为:
单元顶点Ai(ri,λi,φi)的局部坐标Ai(ξi,ηi,ζi)为:
所述步骤(9)的方法为:
记某计算单元8个顶点对应的扰动引力值分别为gi,1、gi,2、gi,3、gi,4、gi+1,1、gi+1,2、gi+1,3、gi+1,4,称12条棱为计算单元上的1-网,定义在其上的函数为1-网函数,记为fi(ξ,η,ζ)(i=0,1,…,11);令L(ξ)、L(η)、L(ζ)分别为关于ξ、η、ζ的一次Lagrange插值算符,其插值基函数为:
令l(3)为3维1-网函数插值算符,则:
l(3)=L(ξ)L(η)+L(η)L(ζ)+L(ζ)L(ξ)-2L(ξ)L(η)L(ζ)
将l(3)作用于1-网函数,即可求得单元内任意一点A(ξ,η,ζ)的值δg(ξ,η,ζ),
至此,经过上述九步可最终建立一种沿临近空间大范围机动弹道空间包络的扰动引力网函数逼近方法,方法可适应沿空间包络内任意弹道的扰动引力快速计算要求。
本发明的有益效果在于:
1、本发明首次针对临近空间大范围机动弹道的扰动引力快速计算问题,提出一种基于弹道空间包络和网函数逼近理论的扰动引力计算方法,方法的计算规模、计算速度明显优于传统扰动引力计算方法,可适应沿弹道空间包络内任意弹道的扰动引力快速计算,计算精度满足临近空间大范围机动弹道的精度要求;
2、本方法基于弹道空间包络计算构建扰动引力重构模型,可有效避免实际弹道因侧向大范围机动而超出逼近区间的问题,能覆盖考虑禁飞区及最大侧向机动等极端情况,对各种情况下的弹道均具有良好的适应性;
3、本方法采用的弹道空间包络计算方法能够根据再入点和目标点快速确定弹道可达域,计算得到的包络可覆盖具各种情况下的机动弹道;
4、本方法引入极点变换的思想以简化弹道空间包络边界的确定以及空域剖分建模过程,引入局部坐标系以简化逼近方法的建立;
5、本发明首次将网函数逼近理论应用到针对临近空间大范围机动弹道的扰动引力快速计算中,网函数的插值函数类具有Coons型结构以及明显的统计特征,便于计算机实现,且具有较高的逼近精度。
附图说明
图1为不同射向弹道侧向包络示意图;
图2为东向不同纵程弹道侧向包络示意图;
图3为滑翔弹道三维包络示意图;
图4为换极坐标系及飞行状态量定义示意图;
图5为弹道再入点I和换极后极点P的位置关系示意图;
图6为换极坐标系航迹偏航角的位置关系示意图;
图7-1为换极坐标系空间包络空域剖分节点示意图之俯视图;
图7-2为换极坐标系空间包络空域剖分节点示意图之侧视图;
图8-1为一般坐标系空间包络空域剖分节点示意图之俯视图;
图8-2为一般坐标系空间包络空域剖分节点示意图之侧视图;
图9为空间包络节点扰动引力天向分量示意图;
图10为空间包络节点扰动引力北向分量示意图;
图11为空间包络节点扰动引力东向分量示意图;
图12为网函数计算单元示意图;
图13为沿不同区域弹道扰动引力逼近结果对比示意图;
图14为沿不同区域弹道扰动引力逼近误差对比示意图;
图15为沿不同射向弹道扰动引力逼近结果对比示意图;
图16为沿不同射向弹道扰动引力逼近误差对比示意图;
图17为沿不同射程弹道扰动引力逼近结果对比示意图;
图18为沿不同射程弹道扰动引力逼近误差对比示意图;
图19为扰动引力引起的弹道终端状态偏差(初始方位角0°~360°)示意图;
图20为逼近误差引起的弹道终端状态偏差(初始方位角0°~360°)示意图。
具体实施方式
下面结合实施例和附图对本发明作进一步说明:
实施例:
基于CAV-H飞行器模型进行仿真,基本仿真条件设置为:①质量m=908kg,参考面积Sm=0.48375m2;②再入点状态参数:速度VI=6500m/s,速度倾角θI=0°,高度HI=80km;③滑翔段终端状态参数:速度Vk=2500m/s,速度倾角θk=0°,高度Hk=30km,终端航迹偏航角偏差不大于±5°;④飞行约束条件:最大热流密度最大动压qmax=100kPa,最大过载nmax=3g;⑤终端结束时距目标点的待飞航程:Stogo=100km。
仿真计算机配置为:Intel(R)Core(TM)i5-3470CPU3.20GHz,内存为3.46GB。软件环境为WindowXP操作系统,计算程序基于VC++6.0开发。
本发明所述沿临近空间大范围机动弹道空间包络的扰动引力逼近方法,包括以下步骤:
(1)建立弹道包络侧向宽度与纵程的数学模型,具体方法为:
令弹道再入点为I,其经度为λI、地心纬度为φI,弹道落点为T,其经度为λT、地心纬度为φT;以由再入点和落点确定的再入大圆弧平面为对称面,称沿射向方向对称面以左弹道为左侧机动弹道,沿射向方向对称面以右弹道为右侧机动弹道;记最大左侧机动弹道距对称面的最大距离为左边界Bl,记最大右侧机动弹道距对称面的最大距离为右边界Br;
针对不同射向的临近空间大范围机动弹道进行仿真,其再入点经度为λI=0°、再入点地心纬度为φI=0°,落点经度λT和落点地心纬度φT见表1:
表1不同射向弹道落点坐标
射向 | 正北 | 正东 | 正南 | 正西 |
经度(°) | 0 | 60 | 0 | -60 |
地心纬度(°) | 60 | 0 | -60 | 0 |
计算得到不同射向弹道侧向包络示意图见图1,不同射向弹道侧向包络边界见表2。由图1和表2可见,由于地球自转和非球形因素的影响,东向弹道包络最大、南北向弹道包络次之,西向弹道包络最小;南北向弹道包络向东有小范围偏移。
表2不同射向弹道侧向包络边界
射向 | 正东 | 正南 | 正西 | 正北 |
左边界Bl(km) | 1597 | 1072 | 1072 | 1072 |
右边界Br(km) | 1597 | 1513 | 1072 | 1528 |
针对射向为正东、纵程不同的临近空间大范围机动弹道进行仿真。弹道再入点经度为λI=0°,再入点地心纬度为φI=0°;落点地心纬度为0°,落点经度分别为50°、60°、70°和80°。
根据给定的飞行器再入点飞行状态参数、弹道终端飞行状态参数、飞行过程约束条件及终端约束条件,针对再入点为(λI,φI)=(0°,0°)、落点为(λki,0°)(i=1,2,3,4,…)的低空大范围机动弹道进行弹道计算,记各落点对应的弹道纵程为Li,由下式(1)计算Li:
Li=Rearccos(sinφIsinφT+cosφIcosφTcos(λT-λI))(1)
其中,Re为地球半径;
计算得到不同落点对应的纵程Li见表3。
表3不同落点对应的纵程
落点经度(°) | 50 | 60 | 70 | 80 |
纵程Li(km) | 5560 | 6671 | 7784 | 8896 |
通过改变再入点初始速度方位角获得最大左侧机动弹道和最大右侧向机动弹道,不同纵程弹道侧向包络示意图见图2。记各落点对应的弹道侧向包络宽度为Wio,由下式(2)计算Wio:
Wio=Bl-Br(2)
考虑建模误差,取Cw倍Wio为实际侧向包络宽度Wi,这里取1.5倍Wio为实际侧向包络宽度Wi。计算得到不同纵程弹道侧向包络宽度见表4。
表4不同纵程弹道侧向包络宽度
对Li(i=1,2,3,4,…)和Wi(i=1,2,3,4,…)进行拟合,可建立如下式(3)所示的弹道包络侧向宽度与纵程的数学模型,确定模型系数为确定模型系数为ai:
具体为:W=-1.064813L+0.000879L2-1.225047e-7L3+4.607571e-12L4,临近空间大范围机动弹道三维包络示意图见图3。
(2)建立换极坐标系,具体方法为:
引入一个换极坐标系:用表示换极坐标系中各物理量,用X表示一般坐标系中各物理量,按如下方式建立换极坐标系:
首先,定义一个再入大圆弧平面作为换极赤道平面:对目标点确定的情况,将再入点和目标点地心矢径构成的再入大圆弧平面作为换极赤道平面;对于目标点未确定的情况,根据再入点位置及速度方位角确定的再入大圆弧平面作为换极赤道面;
然后,基于换极赤道平面定义换极坐标系OE为地心,轴沿再入点地心矢径方向,轴在换极赤道面内垂直于轴指向目标点方向,轴与轴、轴构成右手系,见图4。
(3)建立换极坐标系中飞行器动力学模型,具体方法为:
在换极坐标系中建立以时间为自变量的滑翔飞行器动力学方程,其飞行状态量为换极后的经度地心纬度航迹偏航角速度速度倾角和地心距
其中,Cσ、Cθ为哥氏加速度项,和为牵连加速度项,
其中,
其中,ωe为地球旋转加速度矢量,λp和φp为换极后极点P的经度和地心纬度,AP为P的方位角,见图5。
根据换极坐标系定义,一般坐标系与换极坐标系中地心距、当地速度倾角及速度的定义一致,
定义
其中,ψf为点F的方位角,
由一般坐标系中λ和φ确定换极坐标系中和的表达式为:
由换极坐标系中λ和φ确定一般坐标系中λ和φ的表达式为:
由一般坐标系中σ确定换极坐标系中的表达式为:
其中,
式(13)的位置关系见图6。
(4)计算换极坐标系弹道三维包络,具体方法为:
针对再入点经度λI=100°、再入点地心纬度φI=20°、落点经度λT=165°、落点地心纬度φT=50°的临近空间大范围机动弹道进行计算。根据步骤(3)所述坐标转换关系,计算得到换极弹道参数,由计算得
由式(15)计算换极弹道纵程
其中,Re为地球半径;
由计算得将代入式(16),可计算得到换极弹道侧向包络宽度
其中,模型系数ai由步骤(1)确定;
由计算得将换极弹道侧向包络描述为长度为宽度为的矩形,其边界由式(17)确定:
其中,为换极系侧向包络东向下边界,为换极系侧向包络东向上边界;为换极系侧向包络北向下边界,为换极系侧向包络北向上边界
侧向包络对应的经度范围Δλ及地心纬度范围Δφ由式(18)确定:
以弹道高度范围的Cr倍作为弹道包络的垂向范围
由计算可得,考虑机动和偏差,取
(5)换极坐标系全域空域剖分,具体方法为:
记换极系空间包络确定的空域为Ω,按式(19)计算Ω边界:
其中,为换极系Ω天向下边界,为换极系Ω天向上边界,为换极系Ω东向下边界,为换极系Ω东向上边界,为换极系Ω北向下边界,为换极系Ω北向上边界,dr、dλ和dφ分别为空域剖分间隔。
取计算得到的换极系空间包络空域范围见表5。
表5换极系空间包络空域范围
按照天向东向北向的间隔将由式(19)确定的空间包络空域Ω均匀剖分为q个互不重叠的子域Ωe(e=1,2,…,q),记网格节点坐标为
换极坐标系空间包络空域剖分节点示意图见图7。
(6)建立一般坐标系空间包络扰动引力重构模型,具体方法为:
根据换极系网格节点坐标由步骤(3)所述坐标转换关系计算一般坐标系网格节点坐标N(rgi,λgj,φgk)(gi,gj,gk=0,1,2,…)。一般坐标系空间包络空域剖分节点示意图见图8。
由球谐函数方法计算一般坐标系网格节点扰动引力位TG,
其中,
与TG相应的引力即为扰动引力加速度δg,即
δg=gradTG(23)
则扰动引力加速度δg在天东北坐标系OE-REN中的三个分量δgR、δgE、δgN为:
存储一般坐标系节点位置三分量和节点扰动引力三分量,完成空间包络扰动引力重构模型构建。空间包络节点扰动引力示意图见图9、图10和图11。
(7)当前弹道位置所在网格判断,具体方法为:
令当前弹道位置坐标为A(r,λ,φ),当满足如式(25)所示的条件时,
判断A(r,λ,φ)位于由节点N(rgi,λgj,φgk)、N(rgi,λgj+1,φgk)、N(rgi,λgj,φgk+1)、N(rgi,λgj+1,φgk+1)、N(rgi+1,λgj,φgk)、N(rgi+1,λgj+1,φgk)、N(rgi+1,λgj,φgk+1)、N(rgi+1,λgj+1,φgk+1)所确定的网格中。
(8)建立计算网格内局域坐标系,具体方法为:
局部坐标系由半径rl=rgi+dr/2的球面、经度λl=λgj+dλ/2的子午面、纬度φl=φgk+dφ/2的纬圈的交线组成,原点l(rl,λl,φl)为三交线的交点,原点局部坐标为l(0,0,0),ξ、η、ζ分别沿原点l的天向、北向和东向,则单元内任意点A(r,λ,φ)的局部坐标A(ξ,η,ζ)为:
单元顶点Ai(ri,λi,φi)的局部坐标Ai(ξi,ηi,ζi)为:
(9)网格内部扰动引力网函数逼近计算,具体方法为:
基于网函数逼近理论进行单元内部的逼近计算。记某计算单元8个顶点对应的扰动引力值分别为gi,1、gi,2、gi,3、gi,4、gi+1,1、gi+1,2、gi+1,3、gi+1,4,称12条棱为计算单元上的1-网,定义在其上的函数为1-网函数,记为fi(ξ,η,ζ)(i=0,1,…,11)。网函数计算单元示意图见图12。
令L(ξ)、L(η)、L(ζ)分别为关于ξ、η、ζ的一次Lagrange插值算符,其插值基函数为:
令l(3)为3维1-网函数插值算符,则:
l(3)=L(ξ)L(η)+L(η)L(ζ)+L(ζ)L(ξ)-2L(ξ)L(η)L(ζ)
将l(3)作用于1-网函数,即可求得单元内任意一点A(ξ,η,ζ)的值δg(ξ,η,ζ),
取节点扰动引力采用1080阶球谐函数进行赋值。以1080阶球谐函数计算得到的扰动引力值作为近似真值,对比分析逼近方法在不同情况下的赋值情况。
针对纵程为6700km,初始方位角为90°的滑翔弹道进行计算。弹道再入点和落点分别为:①平原地区λI_A1=115°,φI_A1=35°,λT_A1=180°,φT_A1=16.5°;②丘陵地区λI_A2=110°,φI_A2=27°,λT_A2=173°,φT_A2=13°;③特大山区λI_A3=80°,φI_A3=42°,λT_A3=147°,φT_A3=19.3°。沿不同区域弹道扰动引力逼近结果示意图见图13,逼近误差见图14和表6。
表6不同区域弹道扰动引力逼近精度统计
针对纵程为7000km,弹道起点位于特大山区(λI=80°,φI=42°),初始方位角分别为0°、180°和270°的滑翔弹道进行计算。滑翔弹道终点分别为:①正北λT_σ1=80°,φT_σ1=102.3°;②正南λT_σ2=80°,φT_σ2=-18.3°;③正西λT_σ3=12.5°,φT_σ3=20°。沿不同射向弹道扰动引力逼近结果示意图见图15,逼近误差见图16和表7。
表7不同射向弹道扰动引力逼近精度统计
针对弹道起点位于平原地区(λI=115°,φI=35°),初始方位角为90°,射程分别为6000km、7000km和8000km的滑翔弹道进行计算。滑翔弹道终点分别为:①射程6000km弹道λT_r1=174.3°,φT_r1=19.8°;②射程7000km弹道λT_r2=182.2°,φT_r2=15°;③射程8000km弹道λT_r3=189.9°,φT_r3=10°。沿不同射程弹道扰动引力逼近结果示意图见图17,逼近误差见图18和表8。
表8不同射程弹道扰动引力逼近精度统计
针对弹道起点为λI=0°,φI=0°,初始方位角为i×10°(i=0,1,2,…,35)的7000km射程滑翔弹道进行计算。由扰动引力引起的弹道终端状态偏差见图19。分别采用逼近方法和1080阶球谐函数方法求解弹道计算中的扰动引力,视两者落点偏差之差为逼近误差引起的落点偏差。计算得到逼近误差引起的弹道终端状态偏差见图20。
由图13~图20及表6~表8可见,本发明提出的方法能够适用于不同滑翔弹道的计算。在的网格划分下,沿滑翔弹道的平均逼近误差可控制在0.1mgal~1mgal量级。由逼近误差导致的滑翔弹道终端位置偏差在10m量级,使弹道精度提高约1e-3个量级。
表9采用不同扰动引力计算方法所需的弹道计算时间
计算方法 | 计算时间(s) |
网函数逼近方法 | 61.257 |
1080阶球谐函数方法 | 125970.109 |
Stokes方法 | 104653.258 |
采用不同扰动引力计算方法所需的弹道计算时间见表9。由表9可见,采用本发明方法进行弹道计算中的扰动引力计算,单条弹道的计算时间远小于采用传统扰动引力计算方法所需的弹道计算时间。
综合上述仿真结果可获得以下结论:
1、由本发明确定的方法建立沿临近空间大范围机动弹道空间包络的扰动引力重构模型,对位于不同区域,具有不同射程和不同射向的弹道均有较好的适应性,弹道积分无奇点。
2、通过本发明计算得到的扰动引力具有较高的逼近精度。在10km/1°/1°的网格划分下,可使沿弹道的平均扰动引力逼近误差控制在0.1mgal~1mgal量级,由逼近误差导致的滑翔弹道终端位置偏差控制在10m量级,弹道精度提高约1e-3个量级。
3、本方法具有计算速度快的特征。采用本发明方法进行弹道计算中的扰动引力计算,单条弹道的计算时间远小于采用传统扰动引力计算方法所需的弹道计算时间。
上述实施例只是本发明的较佳实施例,并不是对本发明技术方案的限制,只要是不经过创造性劳动即可在上述实施例的基础上实现的技术方案,均应视为落入本发明专利的权利保护范围内。
Claims (2)
1.一种沿临近空间大范围机动弹道空间包络的扰动引力逼近方法,其特征在于:包括以下步骤:
(1)建立弹道包络侧向宽度与纵程的数学模型;
(2)建立换极坐标系;
(3)建立换极坐标系中飞行器动力学模型;
(4)计算换极坐标系弹道三维包络;
(5)换极坐标系全域空域剖分;
(6)建立一般坐标系空间包络扰动引力重构模型;
(7)当前弹道位置所在网格判断;
(8)建立计算网格内局域坐标系;
(9)网格内部扰动引力网函数逼近计算。
2.根据权利要求1所述的沿临近空间大范围机动弹道空间包络的扰动引力逼近方法,其特征在于:
所述步骤(1)的方法为:
令弹道再入点为I,其经度为λI、地心纬度为φI,弹道落点为T,其经度为λT、地心纬度为φT;以由再入点和落点确定的再入大圆弧平面为对称面,称沿射向方向对称面以左弹道为左侧机动弹道,沿射向方向对称面以右弹道为右侧机动弹道;记最大左侧机动弹道距对称面的最大距离为左边界Bl,记最大右侧机动弹道距对称面的最大距离为右边界Br;
根据给定的飞行器再入点飞行状态参数、弹道终端飞行状态参数、飞行过程约束条件及终端约束条件,针对再入点为(λI,φI)=(0°,0°)、落点为(λki,0°)(i=1,2,3,4,…)的低空大范围机动弹道进行弹道计算,记各落点对应的弹道纵程为Li,由下式(1)计算Li:
Li=Rearccos(sinφIsinφT+cosφIcosφTcos(λT-λI))(1)
其中,Re为地球半径;
通过改变再入点初始速度方位角获得最大左侧机动弹道和最大右侧向机动弹道,记各落点对应的弹道侧向包络宽度为Wio,由下式(2)计算Wio:
Wio=Bl-Br(2)
取Cw倍Wio为实际侧向包络宽度Wi,对Li(i=1,2,3,4,…)和Wi(i=1,2,3,4,…)进行拟合,可建立如下式(3)所示的弹道包络侧向宽度与纵程的数学模型,确定模型系数为ai:
所述步骤(2)的方法为:
引入一个换极坐标系:用表示换极坐标系中各物理量,用X表示一般坐标系中各物理量,按如下方式建立换极坐标系:
首先,定义一个再入大圆弧平面作为换极赤道平面:对目标点确定的情况,将再入点和目标点地心矢径构成的再入大圆弧平面作为换极赤道平面;对于目标点未确定的情况,根据再入点位置及速度方位角确定的再入大圆弧平面作为换极赤道面;
然后,基于换极赤道平面定义换极坐标系OE为地心,轴沿再入点地心矢径方向,轴在换极赤道面内垂直于轴指向目标点方向,轴与轴、轴构成右手系;
所述步骤(3)的方法为:
在换极坐标系中建立以时间为自变量的滑翔飞行器动力学方程,其飞行状态量为换极后的经度地心纬度航迹偏航角速度速度倾角和地心距
其中,Cσ、Cθ为哥氏加速度项,和为牵连加速度项,
其中,
其中,ωe为地球旋转加速度矢量,λp和φp为换极后极点P的经度和地心纬度,AP为P的方位角;
根据换极坐标系定义,一般坐标系与换极坐标系中地心距、当地速度倾角及速度的定义一致,
定义
其中,ψf为点F的方位角,
由一般坐标系中λ和φ确定换极坐标系中和的表达式为:
由换极坐标系中λ和φ确定一般坐标系中λ和φ的表达式为:
由一般坐标系中σ确定换极坐标系中的表达式为:
其中,
所述步骤(4)的方法为:
根据步骤(3)所述坐标转换关系,根据弹道再入点经度λI、再入点地心纬度φI、落点经度λT、落点地心纬度φT计算得到换极弹道参数,由计算得 记换极弹道落点经度为
由式(15)计算换极弹道纵程
其中,Re为地球半径;
将代入式(16),可计算得到换极弹道侧向包络宽度
其中,模型系数ai由步骤(1)确定;
将换极弹道侧向包络描述为长度为宽度为的矩形,其边界由式(17)确定:
其中,为换极系侧向包络东向下边界,为换极系侧向包络东向上边界;为换极系侧向包络北向下边界,为换极系侧向包络北向上边界;
侧向包络对应的经度范围Δλ及地心纬度范围Δφ由式(18)确定:
以弹道高度范围的Cr倍作为弹道包络的垂向范围
所述步骤(5)的方法为:
记换极系空间包络确定的空域为Ω,按式(19)计算Ω边界:
其中,为换极系Ω天向下边界,为换极系Ω天向上边界,为换极系Ω东向下边界,为换极系Ω东向上边界,为换极系Ω北向下边界,为换极系Ω北向上边界,dr、dλ和dφ分别为空域剖分间隔;
按照天向东向北向的间隔将由式(19)确定的空间包络空域Ω均匀剖分为q个互不重叠的子域Ωe(e=1,2,…,q),记网格节点坐标为
所述步骤(6)的方法为:
根据换极系网格节点坐标由步骤(3)所述坐标转换关系计算一般坐标系网格节点坐标
由球谐函数方法计算一般坐标系网格节点扰动引力位TG,
其中,
与TG相应的引力即为扰动引力加速度δg,即
δg=gradTG(23)
则扰动引力加速度δg在天东北坐标系OE-REN中的三个分量δgR、δgE、δgN为:
存储一般坐标系节点位置三分量和节点扰动引力三分量,完成空间包络扰动引力重构模型构建;
所述步骤(7)的方法为:
令当前弹道位置坐标为A(r,λ,φ),当满足如式(25)所示的条件时,
判断A(r,λ,φ)位于由节点N(rgi,λgj,φgk)、N(rgi,λgj+1,φgk)、N(rgi,λgj,φgk+1)、N(rgi,λgj+1,φgk+1)、N(rgi+1,λgj,φgk)、N(rgi+1,λgj+1,φgk)、N(rgi+1,λgj,φgk+1)、N(rgi+1,λgj+1,φgk+1)所确定的网格中;
所述步骤(8)的方法为:
局部坐标系由半径rl=rgi+dr/2的球面、经度λl=λgj+dλ/2的子午面、纬度φl=φgk+dφ/2的纬圈的交线组成,原点l(rl,λl,φl)为三交线的交点,原点局部坐标为l(0,0,0),ξ、η、ζ分别沿原点l的天向、北向和东向,则单元内任意点A(r,λ,φ)的局部坐标A(ξ,η,ζ)为:
单元顶点Ai(ri,λi,φi)的局部坐标Ai(ξi,ηi,ζi)为:
所述步骤(9)的方法为:
记某计算单元8个顶点对应的扰动引力值分别为gi,1、gi,2、gi,3、gi,4、gi+1,1、gi+1,2、gi+1,3、gi+1,4,称12条棱为计算单元上的1-网,定义在其上的函数为1-网函数,记为fi(ξ,η,ζ)(i=0,1,…,11);令L(ξ)、L(η)、L(ζ)分别为关于ξ、η、ζ的一次Lagrange插值算符,其插值基函数为:
令l(3)为3维1-网函数插值算符,则:
l(3)=L(ξ)L(η)+L(η)L(ζ)+L(ζ)L(ξ)-2L(ξ)L(η)L(ζ)
将l(3)作用于1-网函数,即可求得单元内任意一点A(ξ,η,ζ)的值δg(ξ,η,ζ),
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610042316.6A CN105740506B (zh) | 2016-01-21 | 2016-01-21 | 沿临近空间大范围机动弹道空间包络的扰动引力逼近方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201610042316.6A CN105740506B (zh) | 2016-01-21 | 2016-01-21 | 沿临近空间大范围机动弹道空间包络的扰动引力逼近方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105740506A true CN105740506A (zh) | 2016-07-06 |
CN105740506B CN105740506B (zh) | 2018-12-11 |
Family
ID=56247536
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201610042316.6A Active CN105740506B (zh) | 2016-01-21 | 2016-01-21 | 沿临近空间大范围机动弹道空间包络的扰动引力逼近方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105740506B (zh) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106547991A (zh) * | 2016-11-25 | 2017-03-29 | 中国工程物理研究院总体工程研究所 | 沿滑翔弹道的扰动引力重构模型优化方法 |
CN110044210A (zh) * | 2019-04-22 | 2019-07-23 | 中国人民解放军国防科技大学 | 考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法 |
CN110046439A (zh) * | 2019-04-22 | 2019-07-23 | 中国人民解放军国防科技大学 | 考虑高阶扰动引力影响的弹道偏差解析预报算法 |
CN111348223A (zh) * | 2020-05-25 | 2020-06-30 | 北京星际荣耀空间科技有限公司 | 一种控制弹道顶点高度的闭路制导方法、装置及设备 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103838914A (zh) * | 2013-12-30 | 2014-06-04 | 北京航空航天大学 | 一种高超声速飞行器滑翔段弹道解析求解方法 |
CN104035335A (zh) * | 2014-05-27 | 2014-09-10 | 北京航空航天大学 | 基于高精度纵、横程解析预测方法的平稳滑翔再入制导律 |
CN104751012A (zh) * | 2015-04-23 | 2015-07-01 | 中国人民解放军国防科学技术大学 | 沿飞行弹道的扰动引力快速逼近方法 |
CN105138808A (zh) * | 2015-10-19 | 2015-12-09 | 中国人民解放军国防科学技术大学 | 基于摄动理论的滑翔弹道误差传播分析方法 |
-
2016
- 2016-01-21 CN CN201610042316.6A patent/CN105740506B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103838914A (zh) * | 2013-12-30 | 2014-06-04 | 北京航空航天大学 | 一种高超声速飞行器滑翔段弹道解析求解方法 |
CN104035335A (zh) * | 2014-05-27 | 2014-09-10 | 北京航空航天大学 | 基于高精度纵、横程解析预测方法的平稳滑翔再入制导律 |
CN104751012A (zh) * | 2015-04-23 | 2015-07-01 | 中国人民解放军国防科学技术大学 | 沿飞行弹道的扰动引力快速逼近方法 |
CN105138808A (zh) * | 2015-10-19 | 2015-12-09 | 中国人民解放军国防科学技术大学 | 基于摄动理论的滑翔弹道误差传播分析方法 |
Non-Patent Citations (5)
Title |
---|
HUAN,XUEYING AN, ETC.: "Fast local representation of gravity anomaly along flight trajectory", 《PROCEEDINGS OF THE INSTITUTION OF MECHANICAL ENGINEERS PART G JOURNAL OF AEROSPACE ENGINEERING》 * |
YU XIE, LUHUA LIU,ETC.: "Rapid generation of entry trajectories with waypoint and no-fly zone constraints", 《ACTA ASTRONAUTICA》 * |
郑伟, 钱山, 汤国建: "弹道导弹制导计算中扰动引力的快速赋值", 《飞行力学》 * |
郑伟,谢愈,汤国建: "自由段弹道扰动引力计算的球谐函数极点变换", 《宇航学报》 * |
钱山,郑伟,张士峰,蔡洪: "一种弹道导弹再入弹道解析方法", 《飞行力学》 * |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106547991A (zh) * | 2016-11-25 | 2017-03-29 | 中国工程物理研究院总体工程研究所 | 沿滑翔弹道的扰动引力重构模型优化方法 |
CN106547991B (zh) * | 2016-11-25 | 2020-11-03 | 中国工程物理研究院总体工程研究所 | 沿滑翔弹道的扰动引力重构模型优化方法 |
CN110044210A (zh) * | 2019-04-22 | 2019-07-23 | 中国人民解放军国防科技大学 | 考虑任意阶地球非球型引力摄动的闭路制导在线补偿方法 |
CN110046439A (zh) * | 2019-04-22 | 2019-07-23 | 中国人民解放军国防科技大学 | 考虑高阶扰动引力影响的弹道偏差解析预报算法 |
CN111348223A (zh) * | 2020-05-25 | 2020-06-30 | 北京星际荣耀空间科技有限公司 | 一种控制弹道顶点高度的闭路制导方法、装置及设备 |
CN111348223B (zh) * | 2020-05-25 | 2020-08-21 | 北京星际荣耀空间科技有限公司 | 一种控制弹道顶点高度的闭路制导方法、装置及设备 |
Also Published As
Publication number | Publication date |
---|---|
CN105740506B (zh) | 2018-12-11 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN103838914B (zh) | 一种高超声速飞行器滑翔段弹道解析求解方法 | |
CN111306989B (zh) | 一种基于平稳滑翔弹道解析解的高超声速再入制导方法 | |
CN104035335B (zh) | 基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法 | |
CN106547991B (zh) | 沿滑翔弹道的扰动引力重构模型优化方法 | |
CN107544067A (zh) | 一种基于高斯混合近似的高超声速再入飞行器跟踪方法 | |
CN107679655A (zh) | 一种航天发射火箭落点预测系统 | |
CN109344449B (zh) | 航天器月地转移轨道逆向设计方法 | |
CN105740506A (zh) | 沿临近空间大范围机动弹道空间包络的扰动引力逼近方法 | |
CN105718660B (zh) | 临近空间大范围机动弹道三维包络计算方法 | |
CN107861517A (zh) | 基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法 | |
CN104751012A (zh) | 沿飞行弹道的扰动引力快速逼近方法 | |
CN105893659A (zh) | 一种卫星访问预报快速计算方法 | |
CN104215242A (zh) | 一种基于横向游移坐标系的极区惯性导航方法 | |
CN106371312A (zh) | 基于模糊控制器的升力式再入预测‑校正制导方法 | |
CN103453907B (zh) | 基于分层大气模型的行星进入段导航滤波方法 | |
CN105300387B (zh) | 一种火星大气进入段非线性非高斯秩滤波方法 | |
CN105138808A (zh) | 基于摄动理论的滑翔弹道误差传播分析方法 | |
Huifeng et al. | Footprint problem with angle of attack optimization for high lifting reentry vehicle | |
CN105760573B (zh) | 沿临近空间大范围机动弹道的扰动引力延拓逼近方法 | |
Dalle et al. | Flight envelope calculation of a hypersonic vehicle using a first principles-derived model | |
de La Cámara et al. | Polar night vortex breakdown and large-scale stirring in the southern stratosphere | |
CN109870906A (zh) | 一种基于bbo优化人工势场的高速旋翼飞行器路径规划方法 | |
CN103921957B (zh) | 一种探月飞船跳跃式再入的跃起点能量管理方法 | |
CN102508492A (zh) | 一种飞行器在等高航路点间的定高度大圆飞行实现方法 | |
CN103344245A (zh) | 火星进入段imu和甚高频无线电组合导航的ud-skf方法 |
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 |