CN111488646B - 一种旋转地球下高超声速平稳滑翔弹道的解析求解方法 - Google Patents

一种旋转地球下高超声速平稳滑翔弹道的解析求解方法 Download PDF

Info

Publication number
CN111488646B
CN111488646B CN202010135508.8A CN202010135508A CN111488646B CN 111488646 B CN111488646 B CN 111488646B CN 202010135508 A CN202010135508 A CN 202010135508A CN 111488646 B CN111488646 B CN 111488646B
Authority
CN
China
Prior art keywords
coordinate system
generalized
formula
gned
aircraft
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
Application number
CN202010135508.8A
Other languages
English (en)
Other versions
CN111488646A (zh
Inventor
陈万春
张晚晴
余文斌
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN202010135508.8A priority Critical patent/CN111488646B/zh
Publication of CN111488646A publication Critical patent/CN111488646A/zh
Application granted granted Critical
Publication of CN111488646B publication Critical patent/CN111488646B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Automation & Control Theory (AREA)
  • Navigation (AREA)

Abstract

本发明一种旋转地球下高超声速平稳滑翔弹道的解析求解方法，其特征在于：它包括以下几个步骤：步骤一：建立辅助地心旋转参照系；步骤二：建立基于广义经纬坐标的再入动力学模型；步骤三：再入解析解推导方法。本发明优点是迎合了高超声速平稳滑翔弹道纵程远大于横程的特点，更利于线化处理。获得了以能量为自变量，关于纵程、横程和广义航向角的简化动力学模型，有利于进一步解耦处理。具有精度高、计算量小、普适性强等特点，有利于应用于弹道规划、预测制导、控制设计等各个方面。

Description

一种旋转地球下高超声速平稳滑翔弹道的解析求解方法
技术领域
本发明提供了一种旋转地球下高超声速平稳滑翔弹道的解析求解方法，属于航天技术、武器技术领域
背景技术
高超声速滑翔飞行器具有射程远、突防能力强和毁伤威力大等优势，已成为各国研究的热点。高超声速滑翔飞行器具有十分独特的弹道，大部分时间在20km—100km高度的临近空间内无动力飞行，其飞行速度与高度介于常规固定翼飞行器与航天器之间。
再入弹道解析解研究对于探索再入飞行机理、再入制导设计、飞行器性能评估等具有重要理论意义和应用价值，但是，受再入动力学模型强非线性、强耦合性特性的制约，推导高精度的再入弹道解析解十分困难。随着再入飞行器横向机动能力的提高和机动突防任务的提出，滑翔段三维弹道的快速解算成为了研究热点。但多数解析解由于忽略地球自转，以及地球曲率对航向角的影响，在旋转地球背景下精度不高，存在跟踪误差，难以满足精确制导、分析以及修正的要求，具有一定的局限性。
发明内容
针对高超声速飞行器平稳滑翔解析解上述问题，在旋转地球背景下，基于射面所在的大圆构造一个广义赤道，并以此建立一个辅助地心旋转参照系以及相应的基于广义经纬坐标的新再入动力学模型。在此基础上，对模型进行了适当的换元、简化、以及线化处理，从而获得了以能量为自变量，关于纵程、横程和广义航向角的简化动力学模型。进而，利用正则摄动方法简化横程和航向角之间复杂的耦合关系，然后用基于谱分解的LTV解析求法进行求解，推导出精度更高的滑翔解析解。
本发明一种旋转地球下高超声速平稳滑翔弹道的解析求解方法，其特征在于：它包括以下几个步骤：
步骤一：建立辅助地心旋转参照系；
建立地心赤道旋转坐标系(GER)：原点在地心E，ze轴垂直于地球赤道表面，指向北极；xe轴和ye轴在赤道平面内，且相互垂直。该坐标系随地球一起转动，其绕ze轴的旋转角速度为地球自转角速度ωe。
当地北东下坐标系(NED)：定义原点o在飞行器质心M向地面的垂直投影点；x轴指向当地北方，y轴指向当地东方，z轴铅垂向下指向地心。
为了方便推导滑翔解析解，定义一个随地球转动的辅助地心旋转参照系(AGR)：原点在地心E，
Figure BDA0002397154840000021
轴指向飞行器初始位置，
Figure BDA0002397154840000022
轴位于过飞行器与目标点的大圆平面内，且垂直于
Figure BDA0002397154840000023
轴，
Figure BDA0002397154840000024
轴可根据右手定则确定。
同时，定义当地广义北东下坐标系(GNED)坐标系：原点o在飞行器质心M向地面的垂直投影点；
Figure BDA0002397154840000025
轴铅垂向下指向地心，
Figure BDA0002397154840000026
轴指向AGR坐标系中
Figure BDA0002397154840000027
指向的“北”方，
Figure BDA0002397154840000028
轴由右手定则确定。
以AGR坐标系为参照系建立飞行器运动方程。记飞行器位置矢量、速度矢量和加速度矢量为
Figure BDA0002397154840000029
和
Figure BDA00023971548400000210
则有如下运动方程：
Figure BDA00023971548400000211
Figure BDA00023971548400000212
高超声速飞行器在再入过程中受到重力G和气动力Fair作用。另外，由于地球匀速自转，飞行器还受到惯性力的作用如下：
Figure BDA00023971548400000213
Figure BDA00023971548400000214
其中，
Figure BDA00023971548400000215
是由于地球自转引起的牵连加速度，
Figure BDA00023971548400000216
是由地球自转引起的科氏加速度。
Figure BDA00023971548400000217
是AGR坐标系下的地球自转角加速度矢量。故
Figure BDA00023971548400000218
m是飞行器质量，然而上述方程无法直观展示飞行器几乎贴近地球表面的运动模式，因此，下面在AGR坐标系下建立关于广义经纬坐标的质心运动方程组。
步骤二：建立基于广义经纬坐标的再入动力学模型；
1.广义经度、纬度和高度
在AGR坐标系中建立一套描述飞行器运动状态的广义经纬度坐标系统。广义赤道为xe轴与
Figure BDA00023971548400000219
轴构成的平面与地球表面的交线。广义经线为两个端点均在
Figure BDA00023971548400000220
轴上的半个大圆，其所构成的平面垂直于广义赤道平面。令过飞行器初始位置向地面的垂直投影点的广义经线为广义本初子母线。进而，可以定义广义经度
Figure BDA0002397154840000031
广义纬度
Figure BDA0002397154840000032
Figure BDA0002397154840000033
为海拔高度，
Figure BDA0002397154840000034
是飞行器相对地球的速率，广义弹道倾角
Figure BDA0002397154840000035
和广义航向角
Figure BDA0002397154840000036
来描述物体的位置和运动，以GNED坐标系中
Figure BDA0002397154840000037
轴指向方向为基准。
基于此广义经纬度坐标系统，有：
Figure BDA0002397154840000038
Figure BDA0002397154840000039
Figure BDA00023971548400000310
由于广义速度
Figure BDA00023971548400000311
Re是地球平均半径，大小为6378.137km；广义弹道倾角
Figure BDA00023971548400000312
和广义航向角
Figure BDA00023971548400000313
的导数与加速度矢量有关，较为复杂，下文中进行分节介绍。
2.广义速度
设速度矢量为V(粗体)，在AGR坐标系下的坐标
Figure BDA00023971548400000314
其中，
Figure BDA00023971548400000315
和
Figure BDA00023971548400000316
分别为矢量
Figure BDA00023971548400000317
在AGR坐标系下xe，ye和ze轴的分量。
设
Figure BDA00023971548400000318
在GNED坐标系下的坐标
Figure BDA00023971548400000319
其中，
Figure BDA00023971548400000320
和
Figure BDA00023971548400000321
分别为矢量
Figure BDA00023971548400000322
在GNED坐标系下x，y和z轴的分量。
因为
Figure BDA00023971548400000323
其中，
Figure BDA00023971548400000324
是从AGR坐标系到GNED坐标系的坐标转换矩阵。从AGR坐标系到GNED坐标系需坐标系先绕
Figure BDA00023971548400000325
轴转过
Figure BDA00023971548400000326
角度，再绕
Figure BDA00023971548400000327
轴转过
Figure BDA00023971548400000328
角度。所以
Figure BDA0002397154840000041
将式(12)代入式(11)中，可得
Figure BDA0002397154840000042
Figure BDA0002397154840000043
Figure BDA0002397154840000044
又因为
Figure BDA0002397154840000045
上式对时间求导数可得
Figure BDA0002397154840000046
上式右边表示包括惯性力在被所有合外力产生的沿速度方向的加速度分量。将上式展开获得基于广义经纬坐标的表达式，整理得
Figure BDA0002397154840000047
其中，aAGR和aGNED分别表示加速度矢量a在AGR和GNED坐标系下的矢量表示形式。
根据式(5)，aGNED计算公式如下：
Figure BDA0002397154840000048
其中，
Figure BDA0002397154840000049
为气动力矢量在GNED坐标系下的矢量形式；GGNED为重力矢量在GNED坐标系下的矢量形式；
Figure BDA00023971548400000410
为牵连加速度矢量在GNED坐标系下的矢量形式；
Figure BDA00023971548400000411
为科氏加速度矢量在GNED坐标系下的矢量形式。有如下表达形式：
Figure BDA00023971548400000412
Figure BDA0002397154840000051
其中，D＝ρV2SCD/2为阻力；L＝ρV2SCL/2为升力；γ是弹道倾角；ψ是飞行器航向角，以当地北向为基准；σ是倾侧角。
Figure BDA0002397154840000052
Figure BDA0002397154840000053
其中，
Figure BDA0002397154840000054
为地球自转角速度矢量在GNED坐标系下的表达形式，XGNED为位置矢量在GNED坐标系下的表达形式。有：
Figure BDA0002397154840000055
Figure BDA0002397154840000056
Figure BDA0002397154840000057
其中，ωex，ωey，ωez为地球自转角速度矢量在AGR坐标系下的各轴分量；
Figure BDA0002397154840000058
为地球自转角速度矢量在GNED坐标系下的各轴分量。
经过一系列化简，可得
Figure BDA0002397154840000059
其中，g是重力加速度大小。
3.广义弹道倾角
由广义弹道倾角定义得
Figure BDA00023971548400000510
对上式求导可得
Figure BDA0002397154840000061
将式(27)代入上式，为更清晰的表达
Figure BDA0002397154840000062
将
Figure BDA0002397154840000063
拆分为
Figure BDA0002397154840000064
和
Figure BDA0002397154840000065
两项分别进行计算，表示如下
Figure BDA0002397154840000066
其中，
Figure BDA0002397154840000067
Figure BDA0002397154840000068
式(31)中，
Figure BDA0002397154840000069
将式(24)代入上式中，经过整理后代入式(30)中可得
Figure BDA00023971548400000610
假定飞行器采用BTT飞行方式，最终可得
Figure BDA00023971548400000611
4.广义航向角
根据广义航向角定义，可得
Figure BDA00023971548400000612
对上式求导可得：
Figure BDA00023971548400000613
进行一系列化简后，为更清晰的表达
Figure BDA0002397154840000071
将
Figure BDA0002397154840000072
拆分为
Figure BDA0002397154840000073
和
Figure BDA0002397154840000074
两项分别进行计算，表示如下
Figure BDA0002397154840000075
其中，
Figure BDA0002397154840000076
Figure BDA0002397154840000077
将式(19)代入式(39)，进行整理可得
Figure BDA0002397154840000078
5.AGR坐标系下动力学模型
将式(6)-(8)，式(27)，式(35)和式(41)结合起来，即得到以AGR坐标系为参照系建立飞行器运动方程，整理如下。
Figure BDA0002397154840000079
Figure BDA00023971548400000710
Figure BDA00023971548400000711
Figure BDA00023971548400000712
Figure BDA00023971548400000713
Figure BDA00023971548400000714
步骤三：再入解析解推导方法；
1.运动方程线化
以能量为自变量可以忽略速度方程，达到降维的目的。此外，对于无动力滑翔飞行器，在再入过程中能量单调递减，因此可以作为自变量。绝对能量的表达式为
Figure BDA0002397154840000081
忽略地球自转影响，能量对时间的导数为
Figure BDA0002397154840000082
定义纵向射程
Figure BDA0002397154840000083
横向射程
Figure BDA0002397154840000084
以及航向角误差
Figure BDA0002397154840000085
为了推导解析解，将公式(42)-(43)和(47)除以公式(49)，进而得到以
Figure BDA0002397154840000086
为自变量的运动方程组：
Figure BDA0002397154840000087
Figure BDA0002397154840000088
Figure BDA0002397154840000089
由于平稳滑翔过程中弹道倾角变化率很小，令L1＝Lcosσ，L2＝Lsinσ，假设
Figure BDA00023971548400000810
由式(46)易得
Figure BDA00023971548400000811
将上式代入式(50)-(52)中，但令式(52)中
Figure BDA00023971548400000812
分母中的
Figure BDA00023971548400000813
由于再入制导导引飞行器近似沿着广义赤道飞向目标，有
Figure BDA00023971548400000814
和
Figure BDA00023971548400000815
故假设
Figure BDA00023971548400000816
则(50)-(52)可被简化如下
Figure BDA00023971548400000817
Figure BDA0002397154840000091
Figure BDA0002397154840000092
式(54)-(56)中的阻力D可通过D＝L1/(L1/D)来获得。令平稳滑翔弹道倾角变化率为0，可得L1表达式
Figure BDA0002397154840000093
将上式代入阻力D＝L1/(L1/D)，进而代入式(54)-(56)，可得纵程、横程、航向角变化率关于升阻比的表达式，化简如下：
Figure BDA0002397154840000094
Figure BDA0002397154840000095
Figure BDA0002397154840000096
其中，由于
Figure BDA0002397154840000097
故假设
Figure BDA0002397154840000098
根据式(26)，式(58)-(60)中出现的
Figure BDA0002397154840000099
和
Figure BDA00023971548400000910
在此假设下变为：
Figure BDA00023971548400000911
Figure BDA00023971548400000912
其中，φ0为初始纬度，ψLOS0为初始航向角。
上面两式形式复杂，故化简如下。
定义常数α1和β1为
Figure BDA0002397154840000101
所以有
ωex＝ωecosα1sinβ1 (64)
ωey＝ωecosα1cosβ1 (65)
ωez＝ωesinα1 (66)
代入式(61)-(62)可得：
Figure BDA0002397154840000102
Figure BDA0002397154840000103
得到更为简洁的形式。
此外，令
Figure BDA0002397154840000104
Figure BDA0002397154840000105
Figure BDA0002397154840000106
Figure BDA0002397154840000107
其中，
Figure BDA0002397154840000108
为初始能量。
将式(67)-(68)，式(69)-(72)代入式(58)-(60)，可得：
Figure BDA0002397154840000109
Figure BDA00023971548400001010
Figure BDA00023971548400001011
其中，广义经度可由式
Figure BDA0002397154840000111
计算。
2.正则摄动模型
由式(74)-(75)可以看出，横程与航向角之间存在复杂的耦合关系。本节采用正则摄动方法处理这种耦合关系。根据正则摄动方法，定义参量ε为一种标记，且等于常量k。将纵程xD、横程xC和航向角
Figure BDA0002397154840000112
集合成一个矢量x，有：
Figure BDA0002397154840000113
其中，
Figure BDA0002397154840000114
t表示时间。
令式(73)-(75)为
Figure BDA0002397154840000115
将状态量摄动展开
Figure BDA0002397154840000116
其中，
Figure BDA0002397154840000117
用来表示状态量的0阶，1阶和2阶状态。
将式(77)展开到1阶，如下所示。
0阶动力学方程如下
Figure BDA0002397154840000118
Figure BDA0002397154840000119
Figure BDA00023971548400001110
1阶动力学方程如下，
Figure BDA00023971548400001111
Figure BDA00023971548400001112
0阶和1阶状态初始值为：
Figure BDA0002397154840000121
Figure BDA0002397154840000122
最终解析解的形式为：
Figure BDA0002397154840000123
Figure BDA0002397154840000124
Figure BDA0002397154840000125
3.解析解模型求解
由于
Figure BDA0002397154840000126
故在求解解析解时可忽略
Figure BDA0002397154840000127
变化的影响，用一个平均值代替即
Figure BDA0002397154840000128
记R*＝Re+H*。
考虑实际飞行中，升阻比剖面一般是能量的分段低次多项式函数，故规划参考曲线L1/D和L2/D为能量的二阶多项式如下：
Figure BDA0002397154840000129
Figure BDA00023971548400001210
其中，a2、a1、a0、b2、b1、b0为常值参数。
可通过调节攻角和倾侧角跟踪上述参考曲线。
(1)纵程解析解
从式(77)看出，纵向射程xD是独立的，可以单独求解。将R*＝Re+H*代入式(69)中，可得
Figure BDA00023971548400001211
上式分母中的
Figure BDA00023971548400001212
可以用能量的线性函数拟合：
Figure BDA00023971548400001213
其中，参数h11、h10为常值参数，表达式如下：
Figure BDA00023971548400001214
Figure BDA00023971548400001215
其中，
Figure BDA00023971548400001216
表示终端经度；
Figure BDA00023971548400001217
表示终端能量。
此外，当飞行器严格按照大圆飞行时，式(91)分母中出现的
Figure BDA0002397154840000131
可近似为常值。证明如下：
因为
Figure BDA0002397154840000132
与速度无关，所以可忽略速度影响，视飞行器仅受重力作用做圆周运动。由于重力的力矩为0，令飞行器所受总力矩矢量为
Figure BDA0002397154840000133
故可得
Figure BDA0002397154840000134
在GNED坐标系下，由式(24)-(25)，可得
Figure BDA0002397154840000135
其中，
Figure BDA0002397154840000136
为力矩矢量在GNED坐标系下的表达式；XGNED为位置矢量在GNED坐标系下的表达式；VGNED为速度矢量在GNED坐标系下的表达式。
转到AGR坐标系中可得
Figure BDA0002397154840000137
所以
Figure BDA0002397154840000138
为常值。其中，
Figure BDA0002397154840000139
表示从GNED坐标系到AGR坐标系的坐标转换矩阵；
Figure BDA00023971548400001310
为力矩矢量在AGR坐标系下的表达式。
而结合式(26)，可得
Figure BDA00023971548400001311
表达式：
Figure BDA00023971548400001312
结合式(97)，易知
Figure BDA00023971548400001313
为常值。
但在实际飞行过程中，由于飞行器飞行并不严格沿着大圆飞行，
所以导致
Figure BDA00023971548400001314
并不严格为常值。针对这种情况，可以采用能量的线性函数拟合
Figure BDA00023971548400001315
如下：
Figure BDA00023971548400001316
其中，h21、h20为常值参数，表达式如下：
Figure BDA00023971548400001317
Figure BDA0002397154840000141
式(100)中，
Figure BDA0002397154840000142
以及求解
Figure BDA0002397154840000143
所需的
Figure BDA0002397154840000144
和
Figure BDA0002397154840000145
均可由粗略估计得到。
将式(89)，(92)，(99)代入式(79)，积分可得：
Figure BDA0002397154840000146
其中，
Figure BDA0002397154840000147
Figure BDA0002397154840000148
Figure BDA0002397154840000149
c2＝(a1-c1d0)/d1 (106)
c3＝a0-c2d0 (107)
即为纵向射程解析解。
(2)横程与航向角解析解
由式(77)易知，
Figure BDA00023971548400001410
与
Figure BDA00023971548400001411
相互耦合，构成一个复杂的线性时变系统。
对式(80)-(81)联立，得
Figure BDA00023971548400001412
采用基于谱分解的LTV解析求解法，可得横程与航向角的0阶解析解：
Figure BDA00023971548400001413
Figure BDA00023971548400001414
其中，尽管上述0阶解析解中含有积分项，但由于被积曲线平缓，可采用20点的Gauss-Legendre积分公式计算。
将式(82)-(83)联立，为简化式子，定义函数
Figure BDA0002397154840000151
为
Figure BDA0002397154840000152
可得
Figure BDA0002397154840000153
同样采用基于谱分解的LTV解析求解法，结合式(85)，可得横程与航向角的1阶解析解：
Figure BDA0002397154840000154
Figure BDA0002397154840000155
由于上述两式中的积分项有些复杂，故采用N阶Lagrange插值多项式拟合被积函数。
令
Figure BDA0002397154840000156
其中，i为多项式的项数；pi为多项式中第i项的系数；
Figure BDA0002397154840000157
为多项式自变量；
Figure BDA0002397154840000158
为多项式自变量的i次方；Li为第i项的Lagrange插值多项式拟合函数，定义见式(117)。
Figure BDA0002397154840000159
其中，函数
Figure BDA00023971548400001510
定义见式(111)。
Figure BDA00023971548400001511
定义为
Figure BDA00023971548400001512
其中，参数j表示从0到N且不等于i的正整数。
Figure BDA00023971548400001513
取值为
Figure BDA0002397154840000161
取N＝4时，则
Figure BDA0002397154840000162
的拟合多项式为
Figure BDA0002397154840000163
其中，p14，p13，p12，p11，p10为拟合多项式系数。
取N＝4时，误差项为
Figure BDA0002397154840000164
其中，ζ位于
Figure BDA0002397154840000165
区间，且满足积分中值定理。
同理，令
Figure BDA0002397154840000166
其中，
Figure BDA0002397154840000167
同样用4阶Lagrange插值多项式拟合得：
Figure BDA0002397154840000168
其中，p24，p23，p22，p21，p20为拟合多项式系数。
此时误差项为
Figure BDA0002397154840000169
其中，
Figure BDA00023971548400001610
将式(119)和式(123)分别代入式(113)和(114)中，可得
Figure BDA00023971548400001611
Figure BDA00023971548400001612
当将解析解应用到再入制导律中时，用这样的多项式拟合计算积分项将大大减小计算负担，提高再入制导效率。
将式(125)-(126)，式(109)-(110)代入式(87)-(88)中，可得横程、航向角解析解：
Figure BDA0002397154840000171
Figure BDA0002397154840000172
由上式可以看出，最终解析解由0阶解和1阶解组成。其中，1阶解中的多项式系数需要由0阶解确定。
本发明优点及有益效果在于：
(1)现有的三维弹道解析解多未考虑地球曲率和地球自转，精度较差。为推导更精确的三维弹道解析解，基于射面所在的大圆构造了一个广义赤道，并以此建立一个辅助地心旋转参照系以及相应的基于广义经纬坐标的新再入动力学模型。此动力学模型迎合了高超声速平稳滑翔弹道纵程远大于横程的特点，更利于线化处理。
(2)以解析求解为目标，证明了一些与地球自转相关的复杂公式在沿地球大圆飞行条件保持常值。并以此为基础，对动力学模型进行了适当的换元、简化以及线化处理，从而获得了以能量为自变量，关于纵程、横程和广义航向角的简化动力学模型，有利于进一步解耦处理。
(3)利用正则摄动方法简化横程和航向角之间复杂的耦合关系，并用基于谱分解的LTV解析求法进行解耦处理。再利用拉格朗日插值多项式拟合被积项，得到更为简洁的纵程、横程和航向角的三维解析解。此三维解析解具有精度高、计算量小、普适性强等特点，有利于应用于弹道规划、预测制导、控制设计等各个方面。
附图说明
图1是地心赤道旋转坐标系(GER)和当地北东下坐标系(NED)示意图。
图2是辅助地心旋转参照系(AGR)与辅助当地北东下坐标系(GNED)示意图。
图3是函数
Figure BDA0002397154840000173
与拟合插值多项式对比。
图4是函数
Figure BDA0002397154840000174
与拟合插值多项式对比。
图5是GER坐标系下七种目标状态的经纬度解析解仿真对比曲线。
图6是NED坐标系下七种目标状态的航向角解析解仿真对比曲线。
上述图中，涉及到的符号、代号说明如下：
图1中，
Figure BDA0002397154840000181
表示地心赤道旋转坐标系(GER)，oxyz表示当地北东下坐标系(NED)，North表示北向，East表示东向。h是海拔高度，R是射程，λ是经度，φ是纬度，V是飞行器相对于地球的速率，γ是弹道倾角，ψ是飞行器航向角，以当地北向为基准，σ是倾侧角。图2中，
Figure BDA0002397154840000182
表示辅助地心旋转参照系(AGR)，
Figure BDA0002397154840000183
表示当地北东下坐标系(NED)，Axis表示地球极轴。M是飞行器当前位置点，o是飞行器当前位置点在地球表面的投影点。T是目标点，oT是目标点在地球表面的投影点。ωex表示地球自转角速度在
Figure BDA0002397154840000184
轴的分量，ωey表示地球自转角速度在
Figure BDA0002397154840000185
轴的分量，ωzx表示地球自转角速度在
Figure BDA0002397154840000186
轴的分量。图3中，Function表示函数
Figure BDA0002397154840000187
仿真结果，Lagrange Polynomial表示基于Lagrange插值多项式的拟合结果，Energy表示能量。图4中，Function表示函数
Figure BDA0002397154840000188
仿真结果，Lagrange Polynomial表示基于Lagrange插值多项式的拟合结果，Energy表示能量。图5中，Trajectory Simulation表示弹道仿真结果，Analytical Solution表示解析解计算结果。T1-T7分别表示从目标点1到目标点7。Longitude表示经度，Latitude表示纬度。图6中，Trajectory Simulation表示弹道仿真结果，Analytical Solution表示解析解计算结果。T1-T7分别表示从目标点1到目标点7。Heading Angle表示航向角，Energy表示能量。
具体实施方式
下面将结合附图和实施案例对本发明作进一步的详细说明。
针对高超声速飞行器平稳滑翔解析解上述问题，在旋转地球背景下，基于射面所在的大圆构造一个广义赤道，并以此建立一个辅助地心旋转参照系以及相应的基于广义经纬坐标的新再入动力学模型。在此基础上，对模型进行了适当的换元、简化、以及线化处理，从而获得了以能量为自变量，关于纵程、横程和广义航向角的简化动力学模型。进而，利用正则摄动方法简化横程和航向角之间复杂的耦合关系，然后用基于谱分解的LTV解析求法进行求解，推导出精度更高的滑翔解析解。
本发明一种旋转地球下高超声速平稳滑翔弹道的解析求解方法，其特征在于：它包括以下几个步骤：
步骤一：建立辅助地心旋转参照系；
如图1所示，建立地心赤道旋转坐标系(Geocentric Equatorial Rotatingframe,GER)：原点在地心E，ze轴垂直于地球赤道表面，指向北极；xe轴和ye轴在赤道平面内，且相互垂直。该坐标系随地球一起转动，其绕ze轴的旋转角速度为地球自转角速度ωe。
当地北东下坐标系(local North-East-Down frame,local NED frame)：定义原点o在飞行器质心M向地面的垂直投影点；，x轴指向当地北方，y轴指向当地东方，z轴铅垂向下指向地心。
为了方便推导滑翔解析解，如图2所示，定义一个随地球转动的辅助地心旋转参照系(Auxiliary Geocentric Rotating,AGR frame)：原点在地心E，
Figure BDA0002397154840000191
轴指向飞行器初始位置，
Figure BDA0002397154840000192
轴位于过飞行器与目标点的大圆平面内，且垂直于
Figure BDA0002397154840000193
轴，
Figure BDA0002397154840000194
轴可根据右手定则确定。
同时，定义当地广义北东下坐标系(local GeneralizedNorth-East-Down frame,GNEDframe)坐标系：原点o在飞行器质心M向地面的垂直投影点；
Figure BDA0002397154840000195
轴铅垂向下指向地心，
Figure BDA0002397154840000196
轴指向AGR坐标系中
Figure BDA0002397154840000197
指向的“北”方，
Figure BDA0002397154840000198
轴由右手定则确定。
以AGR坐标系为参照系建立飞行器运动方程。记飞行器位置矢量、速度矢量和加速度矢量为
Figure BDA0002397154840000199
和
Figure BDA00023971548400001910
则有如下运动方程：
Figure BDA00023971548400001911
Figure BDA00023971548400001912
高超声速飞行器在再入过程中受到重力G和气动力Fair作用。另外，由于地球匀速自转，飞行器还受到惯性力的作用如下：
Figure BDA00023971548400001913
Figure BDA00023971548400001914
其中，
Figure BDA00023971548400001915
是由于地球自转引起的牵连加速度，
Figure BDA00023971548400001916
是由地球自转引起的科氏加速度。
Figure BDA00023971548400001917
是AGR坐标系下的地球自转角加速度矢量。故
Figure BDA00023971548400001918
然而，上述方程无法直观展示飞行器几乎贴近地球表面的运动模式，因此，下面在AGR坐标系下建立关于广义经纬坐标的质心运动方程组。
步骤二：建立基于广义经纬坐标的再入动力学模型；
1.广义经度、纬度和高度
在AGR坐标系中建立一套描述飞行器运动状态的广义经纬度坐标系统。广义赤道为
Figure BDA0002397154840000201
轴与
Figure BDA0002397154840000202
轴构成的平面与地球表面的交线。广义经线为两个端点均在
Figure BDA0002397154840000203
轴上的半个大圆，其所构成的平面垂直于广义赤道平面。令过飞行器初始位置向地面的垂直投影点的广义经线为广义本初子母线。进而，可以定义广义经度
Figure BDA0002397154840000204
广义纬度
Figure BDA0002397154840000205
广义高度
Figure BDA0002397154840000206
广义速度
Figure BDA0002397154840000207
广义弹道倾角
Figure BDA0002397154840000208
和广义航向角
Figure BDA0002397154840000209
来描述物体的位置和运动。
基于此广义经纬度坐标系统，我们易知，AGR坐标系中描述位置矢量和速度矢量的方式与GER坐标系一致。但由于AGR坐标系中的广义赤道平面与地球赤道平面不重合，不同于GER坐标系下只有z轴有旋转角速度，地球自转加速度矢量在AGR坐标系下的三个轴上均有分量，故加速度矢量的形式对二者而言并不相同。由于广义经度
Figure BDA00023971548400002010
广义纬度
Figure BDA00023971548400002011
和广义高度
Figure BDA00023971548400002012
对时间的导数仅与位置矢量和速度矢量有关，故其与GER坐标系下形式一致，为：
Figure BDA00023971548400002013
Figure BDA00023971548400002014
Figure BDA00023971548400002015
由于广义速度
Figure BDA00023971548400002016
广义弹道倾角
Figure BDA00023971548400002017
和广义航向角
Figure BDA00023971548400002018
的导数与加速度矢量有关，较为复杂，下文中进行分节介绍。
2.广义速度
设广义速度矢量
Figure BDA00023971548400002019
在AGR坐标系下的坐标
Figure BDA00023971548400002020
设广义速度矢量
Figure BDA00023971548400002021
在GNED坐标系下的坐标
Figure BDA00023971548400002022
因为
Figure BDA0002397154840000211
其中，
Figure BDA0002397154840000212
是从AGR坐标系到GNED坐标系的坐标转换矩阵。从AGR坐标系到GNED坐标系需坐标系先绕
Figure BDA0002397154840000213
轴转过
Figure BDA0002397154840000214
角度，再绕
Figure BDA0002397154840000215
轴转过
Figure BDA0002397154840000216
角度。所以
Figure BDA0002397154840000217
将式(12)代入式(11)中，可得
Figure BDA0002397154840000218
Figure BDA0002397154840000219
Figure BDA00023971548400002110
又因为
Figure BDA00023971548400002111
上式对时间求导数可得
Figure BDA00023971548400002112
上式右边表示包括惯性力在被所有合外力产生的沿速度方向的加速度分量。下面将上式展开获得基于广义经纬坐标的表达式。由公式(13)-(15)可得
Figure BDA00023971548400002113
Figure BDA00023971548400002114
Figure BDA00023971548400002115
将式(18)-(20)代入式(17)，整理得
Figure BDA00023971548400002116
根据式(5)，aGNED计算公式如下：
Figure BDA0002397154840000221
其中，
Figure BDA0002397154840000222
Figure BDA0002397154840000223
Figure BDA0002397154840000224
Figure BDA0002397154840000225
其中
Figure BDA0002397154840000226
Figure BDA0002397154840000227
下面给出
Figure BDA0002397154840000228
求解过程。为求
Figure BDA0002397154840000229
需求
Figure BDA00023971548400002210
由图2易知，GER坐标系下的地球自转角速度
Figure BDA00023971548400002211
表示形式比较简单，为：
Figure BDA00023971548400002212
因此，可以通过
Figure BDA00023971548400002213
利用坐标转换关系得到
Figure BDA00023971548400002214
GER坐标系与AGR坐标系有如下转换关系：
Figure BDA00023971548400002215
其中，
Figure BDA00023971548400002216
是从GER坐标系到AGR坐标系的坐标转换矩阵。设初始时刻飞行器在GER坐标系下的经度为λ0，纬度为φ0，视线航向角为ψLOS0，则从GER坐标系到AGR坐标系需坐标系先绕z轴转过λ0角度，再绕y轴转过-φ0角度，最后绕x轴旋转90-ψLOS0度。所以可得
Figure BDA00023971548400002217
将上式代入式(30)可得
Figure BDA0002397154840000231
其中，ωex，ωey，ωez均为常值。进一步可以获得
Figure BDA0002397154840000232
即为ωe在GNED下的分量，如下：
Figure BDA0002397154840000233
将式(24)-(25)，和式(26)代入式(22)-(23)可得
Figure BDA0002397154840000234
Figure BDA0002397154840000235
将式(34)-(35)代入式(18)中，经过一系列化简，可得
Figure BDA0002397154840000236
上式右边第一项是气动力的影响，第二项是重力的影响，剩下的项是地球自转引起的牵连加速度的影响。
3.广义弹道倾角
由广义弹道倾角定义得
Figure BDA0002397154840000237
对上式求导可得
Figure BDA0002397154840000238
将式(20)与式(27)代入上式，可得
Figure BDA0002397154840000239
其中，
Figure BDA0002397154840000241
Figure BDA0002397154840000242
式(31)中，
Figure BDA0002397154840000243
将式(24)，式(34)-(35)代入上式中，经过整理后代入式(30)中可得
Figure BDA0002397154840000244
假定飞行器采用BTT飞行方式，进行化简，最终可得
Figure BDA0002397154840000245
4.广义航向角
根据广义航向角定义，可得
Figure BDA0002397154840000246
对上式求导可得：
Figure BDA0002397154840000247
将式(18)-(19)代入上式，进行一系列化简可整理为
Figure BDA0002397154840000248
其中，
Figure BDA0002397154840000249
Figure BDA00023971548400002410
将式(34)-(35)，及式(19)代入式(39)，进行整理可得
Figure BDA0002397154840000251
5.模型整理
将上面得到的方程整理如下。
Figure BDA0002397154840000252
Figure BDA0002397154840000253
Figure BDA0002397154840000254
Figure BDA0002397154840000255
Figure BDA0002397154840000256
Figure BDA0002397154840000257
其中，
Figure BDA0002397154840000258
为广义经度，
Figure BDA0002397154840000259
为广义纬度，
Figure BDA00023971548400002510
是海拔高度，
Figure BDA00023971548400002511
是飞行器相对地球的速率，
Figure BDA00023971548400002512
是广义弹道倾角，
Figure BDA00023971548400002513
是广义航向角，以GNED坐标系中
Figure BDA00023971548400002514
轴指向方向为基准。m是飞行器质量，g是重力加速度大小，Re是地球平均半径，大小为6378.137km。
Figure BDA00023971548400002515
为地球自转角速度矢量在GNED坐标系下的各轴分量。L＝ρV2SCL/2为升力，D＝ρV2SCD/2为阻力。
步骤三：再入解析解推导方法；
1.运动方程线化
以能量为自变量可以忽略速度方程，达到降维的目的。此外，对于无动力滑翔飞行器，在再入过程中能量单调递减，因此可以作为自变量。绝对能量的表达式为
Figure BDA0002397154840000261
忽略地球自转影响，能量对时间的导数为
Figure BDA0002397154840000262
定义纵向射程
Figure BDA0002397154840000263
横向射程
Figure BDA0002397154840000264
以及航向角误差
Figure BDA0002397154840000265
为了推导解析解，将公式(42)-(43)和(47)除以公式(49)，进而得到以
Figure BDA0002397154840000266
为自变量的运动方程组：
Figure BDA0002397154840000267
Figure BDA0002397154840000268
Figure BDA0002397154840000269
由于平稳滑翔过程中弹道倾角变化率很小，假设
Figure BDA00023971548400002610
由式(46)易得
Figure BDA00023971548400002611
将上式代入式(50)-(52)中，但令式(52)中
Figure BDA00023971548400002612
分母中的
Figure BDA00023971548400002613
由于再入制导导引飞行器近似沿着广义赤道飞向目标，有
Figure BDA00023971548400002614
和
Figure BDA00023971548400002615
故假设
Figure BDA00023971548400002616
则(50)-(52)可被简化如下
Figure BDA00023971548400002617
Figure BDA00023971548400002618
Figure BDA00023971548400002619
式(54)-(56)中的阻力D可通过D＝L1/(L1/D)来获得。令平稳滑翔弹道倾角变化率为0，可得L1表达式
Figure BDA0002397154840000271
将上式代入阻力D＝L1/(L1/D)，进而代入式(54)-(56)，可得纵程、横程、航向角变化率关于升阻比的表达式，化简如下：
Figure BDA0002397154840000272
Figure BDA0002397154840000273
Figure BDA0002397154840000274
其中，由于
Figure BDA0002397154840000275
故假设
Figure BDA0002397154840000276
根据式(26)，式(58)-(60)中出现的
Figure BDA0002397154840000277
和
Figure BDA0002397154840000278
在此假设下变为：
Figure BDA0002397154840000279
Figure BDA00023971548400002710
上面两式形式复杂，故化简如下。
定义常数α1和β1使得
ωex＝ωecosα1sinβ1 (72)
ωey＝ωecosα1cosβ1 (73)
ωez＝ωesinα1 (74)
代入式(61)-(62)可得：
Figure BDA0002397154840000281
Figure BDA0002397154840000282
得到更为简洁的形式。
此外，令
Figure BDA0002397154840000283
Figure BDA0002397154840000284
Figure BDA0002397154840000285
Figure BDA0002397154840000286
将式(67)-(68)，式(69)-(72)代入式(58)-(60)，可得：
Figure BDA0002397154840000287
Figure BDA0002397154840000288
Figure BDA0002397154840000289
其中，广义经度可由式
Figure BDA00023971548400002810
计算。
2.正则摄动模型
由式(74)-(75)可以看出，横程与航向角之间存在复杂的耦合关系。本节采用正则摄动方法处理这种耦合关系。根据正则摄动方法，定义参量ε为一种标记，且等于常量k，则状态量变为：
Figure BDA00023971548400002811
其中，
Figure BDA00023971548400002812
令式(73)-(75)为
Figure BDA0002397154840000291
将状态量摄动展开
Figure BDA0002397154840000292
其中，
Figure BDA0002397154840000293
用来表示状态量的0阶，1阶和2阶状态。将式(77)展开到1阶，如下所示。
0阶动力学方程如下
Figure BDA0002397154840000294
Figure BDA0002397154840000295
Figure BDA0002397154840000296
1阶动力学方程如下，
Figure BDA0002397154840000297
Figure BDA0002397154840000298
0阶和1阶状态初始值为：
Figure BDA0002397154840000299
Figure BDA00023971548400002910
最终解析解的形式为：
Figure BDA00023971548400002911
Figure BDA00023971548400002912
Figure BDA00023971548400002913
3.模型求解
由于
Figure BDA0002397154840000301
故在求解解析解时可忽略
Figure BDA0002397154840000302
变化的影响，用一个平均值代替即
Figure BDA0002397154840000303
记R*＝Re+H*。
考虑实际飞行中，升阻比剖面一般是能量的分段低次多项式函数，故规划参考曲线L1/D和L2/D为能量的二阶多项式如下：
Figure BDA0002397154840000304
Figure BDA0002397154840000305
可通过调节攻角和倾侧角跟踪上述参考曲线。
(1)纵程解析解
易从式(77)看出，纵程xD是独立的，可以单独求解。将R*＝Re+H*代入式(69)中，可得
Figure BDA0002397154840000306
上式分母中的
Figure BDA0002397154840000307
可以用能量的线性函数拟合：
Figure BDA0002397154840000308
其中，
Figure BDA0002397154840000309
Figure BDA00023971548400003010
此外，当飞行器严格按照大圆飞行时，式(91)分母中出现的
Figure BDA00023971548400003011
可近似为常值。证明如下：
因为
Figure BDA00023971548400003012
与速度无关，所以可忽略速度影响，视飞行器仅受重力作用做圆周运动。由于重力的力矩为0，故可得
Figure BDA00023971548400003013
在GNED坐标系下，由式(24)-(25)，可得
Figure BDA00023971548400003014
转到AGR坐标系中可得
Figure BDA0002397154840000311
所以
Figure BDA0002397154840000312
为常值。
而结合式(26)，可得
Figure BDA0002397154840000313
表达式：
Figure BDA0002397154840000314
结合式(97)，(32)，易知
Figure BDA0002397154840000315
为常值。
但在实际飞行过程中，由于飞行器飞行并不严格沿着大圆飞行，
所以导致
Figure BDA0002397154840000316
并不严格为常值。针对这种情况，可以采用能量的线性函数拟合
Figure BDA0002397154840000317
如下：
Figure BDA0002397154840000318
其中，
Figure BDA0002397154840000319
Figure BDA00023971548400003110
式(100)中，
Figure BDA00023971548400003111
以及求解
Figure BDA00023971548400003112
所需的
Figure BDA00023971548400003113
和
Figure BDA00023971548400003114
均可由粗略估计得到。
将式(89)，(92)，(99)代入式(79)，积分可得：
Figure BDA00023971548400003115
其中，
Figure BDA00023971548400003116
Figure BDA00023971548400003117
Figure BDA0002397154840000321
c2＝(a1-c1d0)/d1 (114)
c3＝a0-c2d0 (115)
即为纵向射程解析解。
(2)横程与航向角解析解
由式(77)易知，
Figure BDA0002397154840000322
与
Figure BDA0002397154840000323
相互耦合，构成一个复杂的线性时变系统。
对式(80)-(81)联立，得
Figure BDA0002397154840000324
采用基于谱分解的LTV解析求解法，可得横程与航向角的0阶解析解：
Figure BDA0002397154840000325
Figure BDA0002397154840000326
其中，尽管上述0阶解析解中含有积分项，但由于被积曲线平缓，可采用20点的Gauss-Legendre积分公式计算。
将式(82)-(83)联立，且令
Figure BDA0002397154840000327
可得
Figure BDA0002397154840000328
同样采用基于谱分解的LTV解析求解法，结合式(85)，可得横程与航向角的1阶解析解：
Figure BDA0002397154840000329
Figure BDA0002397154840000331
由于上述两式中的积分项有些复杂，故采用N阶Lagrange插值多项式拟合被积函数。
令
Figure BDA0002397154840000332
其中，
Figure BDA0002397154840000333
Figure BDA0002397154840000334
定义为
Figure BDA0002397154840000335
Figure BDA0002397154840000336
取值为
Figure BDA0002397154840000337
取N＝4时，则
Figure BDA0002397154840000338
的拟合多项式为
Figure BDA0002397154840000339
取N＝4时，误差项为
Figure BDA00023971548400003310
其中，
Figure BDA00023971548400003311
由于f5 (5)(ζ)形式复杂，故这里采用仿真对比原函数与拟合后的多项式函数来进行误差分析，仿真见图3。
同理，令
Figure BDA00023971548400003312
其中，
Figure BDA00023971548400003313
同样用4阶Lagrange插值多项式拟合得：
Figure BDA00023971548400003314
此时误差项为
Figure BDA0002397154840000341
其中，
Figure BDA0002397154840000342
仿真对比曲线见图4。
将式(119)和式(123)分别代入式(113)和(114)中，可得
Figure BDA0002397154840000343
Figure BDA0002397154840000344
当将解析解应用到再入制导律中时，用这样的多项式拟合计算积分项将大大减小计算负担，提高再入制导效率。
将式(125)-(126)，式(109)-(110)代入式(87)-(88)中，可得横程、航向角解析解：
Figure BDA0002397154840000345
Figure BDA0002397154840000346
由上式可以看出，最终解析解由0阶解和1阶解组成。其中，1阶解中的多项式系数需要由0阶解确定。
实施案例：
为了验证算法的精度，以CAV-H作为再入飞行器模型在七种不同目标状态下进行仿真。在七种不同目标状态下进行仿真。再入初始状态均为h0＝60km，V0＝7200m/s，λ0＝0°，θ0＝45°，γ0＝0°。ψ0依据目标方向取值。采用弹道阻尼控制技术抑制弹道振荡，其中弹道阻尼反馈项系数Kγ＝5。终端约束为ETAEM＝-5.5×107J/kg。仿真结果如图5-图6。
表1终端状态量
Figure BDA0002397154840000347
Figure BDA0002397154840000351
Figure BDA0002397154840000361
图5中，往东飞的弹道由于受到科氏力的影响而产生弯曲。从图5和图6仿真结果可以看出，在升阻比规划为能量多项式的一般情况下，新解析解精度也很高。同时仿真结果表明，新解析解的高精度不受目标状态的影响。当总航程为1万多公里时，解析弹道和数值仿真弹道最终位置误差在2％以内，航向角误差除了T5在6.07％左右，其余均在0.9％以内。
综上所述，通过上述步骤，推导出了本发明方法，即一种旋转地球下高超声速平稳滑翔弹道的解析求解方法，案例仿真结果表明本发明方法能够精准推导得到高精度的解析解，该方法精度高、计算量小、普适性强，十分适用于应用到再入制导、控制设计中。

Claims (1)

1.一种旋转地球下高超声速平稳滑翔弹道的解析求解方法，其特征在于：它包括以下几个步骤：
步骤一：建立辅助地心旋转参照系；
建立地心赤道旋转坐标系GER：原点在地心E，ze轴垂直于地球赤道表面，指向北极；xe轴和ye轴在赤道平面内，且相互垂直；该坐标系随地球一起转动，其绕ze轴的旋转角速度为地球自转角速度ωe；
当地北东下坐标系NED：定义原点o在飞行器质心M向地面的垂直投影点；x轴指向当地北方，y轴指向当地东方，z轴铅垂向下指向地心；
为了方便推导滑翔解析解，定义一个随地球转动的辅助地心旋转参照系AGR：原点在地心E，
Figure FDA0002397154830000011
轴指向飞行器初始位置，
Figure FDA0002397154830000012
轴位于过飞行器与目标点的大圆平面内，且垂直于
Figure FDA0002397154830000013
轴，
Figure FDA0002397154830000014
轴根据右手定则确定；
同时，定义当地广义北东下坐标系GNED坐标系：原点o在飞行器质心M向地面的垂直投影点；
Figure FDA0002397154830000015
轴铅垂向下指向地心，
Figure FDA0002397154830000016
轴指向AGR坐标系中
Figure FDA0002397154830000017
指向的“北”方，
Figure FDA0002397154830000018
轴由右手定则确定；
以AGR坐标系为参照系建立飞行器运动方程；记飞行器位置矢量、速度矢量和加速度矢量为
Figure FDA0002397154830000019
和
Figure FDA00023971548300000110
则有如下运动方程：
Figure FDA00023971548300000111
Figure FDA00023971548300000112
高超声速飞行器在再入过程中受到重力G和气动力Fair作用；另外，由于地球匀速自转，飞行器还受到惯性力的作用如下：
Figure FDA00023971548300000113
Figure FDA00023971548300000114
其中，
Figure FDA00023971548300000115
是由于地球自转引起的牵连加速度，
Figure FDA00023971548300000116
是由地球自转引起的科氏加速度；
Figure FDA00023971548300000117
是AGR坐标系下的地球自转角加速度矢量；故
Figure FDA00023971548300000118
m是飞行器质量，然而上述方程无法直观展示飞行器几乎贴近地球表面的运动模式，因此，下面在AGR坐标系下建立关于广义经纬坐标的质心运动方程组；
步骤二：建立基于广义经纬坐标的再入动力学模型；
2.1广义经度、纬度和高度
在AGR坐标系中建立一套描述飞行器运动状态的广义经纬度坐标系统；广义赤道为
Figure FDA0002397154830000021
轴与
Figure FDA0002397154830000022
轴构成的平面与地球表面的交线；广义经线为两个端点均在
Figure FDA0002397154830000023
轴上的半个大圆，其所构成的平面垂直于广义赤道平面；令过飞行器初始位置向地面的垂直投影点的广义经线为广义本初子母线；进而，定义广义经度
Figure FDA0002397154830000024
广义纬度
Figure FDA0002397154830000025
Figure FDA0002397154830000026
为海拔高度，
Figure FDA0002397154830000027
是飞行器相对地球的速率，广义弹道倾角
Figure FDA0002397154830000028
和广义航向角
Figure FDA0002397154830000029
来描述物体的位置和运动，以GNED坐标系中
Figure FDA00023971548300000210
轴指向方向为基准；
基于此广义经纬度坐标系统，有：
Figure FDA00023971548300000211
Figure FDA00023971548300000212
Figure FDA00023971548300000213
由于Re是地球平均半径，大小为6378.137km；广义弹道倾角
Figure FDA00023971548300000214
和广义航向角
Figure FDA00023971548300000215
的导数与加速度矢量有关，较为复杂，下文中进行分节介绍；
2.2广义速度
设速度矢量为
Figure FDA00023971548300000216
在AGR坐标系下的坐标
Figure FDA00023971548300000217
其中，
Figure FDA00023971548300000218
和
Figure FDA00023971548300000219
分别为矢量
Figure FDA00023971548300000220
在AGR坐标系下xe，ye和ze轴的分量；
设
Figure FDA00023971548300000221
在GNED坐标系下的坐标
Figure FDA00023971548300000222
其中，
Figure FDA00023971548300000223
和
Figure FDA00023971548300000224
分别为矢量
Figure FDA00023971548300000225
在GNED坐标系下x，y和z轴的分量；
因为
Figure FDA0002397154830000031
其中，
Figure FDA0002397154830000032
是从AGR坐标系到GNED坐标系的坐标转换矩阵；从AGR坐标系到GNED坐标系需坐标系先绕
Figure FDA0002397154830000033
轴转过
Figure FDA0002397154830000034
角度，再绕
Figure FDA0002397154830000035
轴转过
Figure FDA0002397154830000036
角度；所以
Figure FDA0002397154830000037
将式(12)代入式(11)中，可得
Figure FDA0002397154830000038
Figure FDA0002397154830000039
Figure FDA00023971548300000310
又因为
Figure FDA00023971548300000311
上式对时间求导数可得
Figure FDA00023971548300000312
上式右边表示包括惯性力在被所有合外力产生的沿速度方向的加速度分量；将上式展开获得基于广义经纬坐标的表达式，整理得
Figure FDA00023971548300000313
其中，aAGR和aGNED分别表示加速度矢量a在AGR和GNED坐标系下的矢量表示形式；
根据式(5)，aGNED计算公式如下：
Figure FDA00023971548300000314
其中，
Figure FDA00023971548300000315
为气动力矢量在GNED坐标系下的矢量形式；GGNED为重力矢量在GNED坐标系下的矢量形式；
Figure FDA00023971548300000316
为牵连加速度矢量在GNED坐标系下的矢量形式；
Figure FDA00023971548300000317
为科氏加速度矢量在GNED坐标系下的矢量形式；有如下表达形式：
Figure FDA0002397154830000041
Figure FDA0002397154830000042
其中，D＝ρV2SCD/2为阻力；L＝ρV2SCL/2为升力；γ是弹道倾角；ψ是飞行器航向角，以当地北向为基准；σ是倾侧角；
Figure FDA0002397154830000043
Figure FDA0002397154830000044
其中，
Figure FDA0002397154830000045
为地球自转角速度矢量在GNED坐标系下的表达形式，XGNED为位置矢量在GNED坐标系下的表达形式；有：
Figure FDA0002397154830000046
Figure FDA0002397154830000047
Figure FDA0002397154830000048
其中，ωex，ωey，ωez为地球自转角速度矢量在AGR坐标系下的各轴分量；
Figure FDA0002397154830000049
为地球自转角速度矢量在GNED坐标系下的各轴分量；
经过化简，可得
Figure FDA00023971548300000410
其中，g是重力加速度大小；
2.3广义弹道倾角
由广义弹道倾角定义得
Figure FDA0002397154830000051
对上式求导可得
Figure FDA0002397154830000052
将式(27)代入上式，为更清晰的表达
Figure FDA0002397154830000053
将
Figure FDA0002397154830000054
拆分为
Figure FDA0002397154830000055
和
Figure FDA0002397154830000056
两项分别进行计算，表示如下
Figure FDA0002397154830000057
其中，
Figure FDA0002397154830000058
Figure FDA0002397154830000059
式(31)中，
Figure FDA00023971548300000510
将式(24)代入上式中，经过整理后代入式(30)中可得
Figure FDA00023971548300000511
设飞行器采用BTT飞行方式，最终可得
Figure FDA00023971548300000512
2.4广义航向角
根据广义航向角定义，可得
Figure FDA0002397154830000061
对上式求导可得：
Figure FDA0002397154830000062
进行一系列化简后，为更清晰的表达
Figure FDA0002397154830000063
将
Figure FDA0002397154830000064
拆分为
Figure FDA0002397154830000065
和
Figure FDA0002397154830000066
两项分别进行计算，表示如下
Figure FDA0002397154830000067
其中，
Figure FDA0002397154830000068
Figure FDA0002397154830000069
将式(19)代入式(39)，进行整理可得
Figure FDA00023971548300000610
2.5 AGR坐标系下动力学模型
将式(6)-(8)，式(27)，式(35)和式(41)结合起来，即得到以AGR坐标系为参照系建立飞行器运动方程，整理如下；
Figure FDA00023971548300000611
Figure FDA00023971548300000612
Figure FDA00023971548300000613
Figure FDA00023971548300000614
Figure FDA0002397154830000071
Figure FDA0002397154830000072
步骤三：再入解析解推导方法；
3.1运动方程线化
以能量为自变量忽略速度方程，达到降维的目的；此外，对于无动力滑翔飞行器，在再入过程中能量单调递减，因此作为自变量；绝对能量的表达式为
Figure FDA0002397154830000073
忽略地球自转影响，能量对时间的导数为
Figure FDA0002397154830000074
定义纵向射程
Figure FDA0002397154830000075
横向射程
Figure FDA0002397154830000076
以及航向角误差
Figure FDA0002397154830000077
为了推导解析解，将公式(42)-(43)和(47)除以公式(49)，进而得到以
Figure FDA0002397154830000078
为自变量的运动方程组：
Figure FDA0002397154830000079
Figure FDA00023971548300000710
Figure FDA00023971548300000711
由于平稳滑翔过程中弹道倾角变化率很小，令L1＝Lcosσ，L2＝Lsinσ，设
Figure FDA00023971548300000712
由式(46)易得
Figure FDA0002397154830000081
将上式代入式(50)-(52)中，但令式(52)中
Figure FDA0002397154830000082
分母中的
Figure FDA0002397154830000083
由于再入制导导引飞行器近似沿着广义赤道飞向目标，有
Figure FDA0002397154830000084
和
Figure FDA0002397154830000085
故设
Figure FDA0002397154830000086
则(50)-(52)被简化如下
Figure FDA0002397154830000087
Figure FDA0002397154830000088
Figure FDA0002397154830000089
式(54)-(56)中的阻力D通过D＝L1/(L1/D)来获得；令平稳滑翔弹道倾角变化率为0，可得L1表达式
Figure FDA00023971548300000810
将上式代入阻力D＝L1/(L1/D)，进而代入式(54)-(56)，可得纵程、横程、航向角变化率关于升阻比的表达式，化简如下：
Figure FDA00023971548300000811
Figure FDA00023971548300000812
Figure FDA00023971548300000813
其中，由于
Figure FDA0002397154830000091
故设
Figure FDA0002397154830000092
根据式(26)，式(58)-(60)中出现的
Figure FDA0002397154830000093
和
Figure FDA0002397154830000094
在此变为：
Figure FDA0002397154830000095
Figure FDA0002397154830000096
其中，φ0为初始纬度，ψLOS0为初始航向角；
上面两式形式复杂，故化简如下；
定义常数α1和β1为
Figure FDA0002397154830000097
所以有
ωex＝ωecosα1sinβ1 (64)
ωey＝ωecosα1cosβ1 (65)
ωez＝ωesinα1 (66)
代入式(61)-(62)可得：
Figure FDA0002397154830000098
Figure FDA0002397154830000099
得到更为简洁的形式；
此外，令
Figure FDA00023971548300000910
Figure FDA00023971548300000911
Figure FDA0002397154830000101
Figure FDA0002397154830000102
其中，
Figure FDA0002397154830000103
为初始能量；
将式(67)-(68)，式(69)-(72)代入式(58)-(60)，可得：
Figure FDA0002397154830000104
Figure FDA0002397154830000105
Figure FDA0002397154830000106
其中，广义经度由式
Figure FDA0002397154830000107
计算；
3.2正则摄动模型
由式(74)-(75)看出，横程与航向角之间存在复杂的耦合关系；采用正则摄动方法处理这种耦合关系；根据正则摄动方法，定义参量ε为一种标记，且等于常量k；将纵程xD、横程xC和航向角
Figure FDA0002397154830000108
集合成一个矢量x，有：
Figure FDA0002397154830000109
其中，
Figure FDA00023971548300001010
t表示时间；
令式(73)-(75)为
Figure FDA00023971548300001011
将状态量摄动展开
Figure FDA00023971548300001012
其中，
Figure FDA00023971548300001013
用来表示状态量的0阶，1阶和2阶状态；
将式(77)展开到1阶，如下所示；
0阶动力学方程如下
Figure FDA0002397154830000111
Figure FDA0002397154830000112
Figure FDA0002397154830000113
1阶动力学方程如下，
Figure FDA0002397154830000114
Figure FDA0002397154830000115
0阶和1阶状态初始值为：
Figure FDA0002397154830000116
Figure FDA0002397154830000117
最终解析解的形式为：
Figure FDA0002397154830000118
Figure FDA0002397154830000119
Figure FDA00023971548300001110
3.3解析解模型求解
由于
Figure FDA00023971548300001111
故在求解解析解时忽略
Figure FDA00023971548300001112
变化的影响，用一个平均值代替即
Figure FDA00023971548300001113
记R*＝Re+H*；
实际飞行中，升阻比剖面是能量的分段低次多项式函数，故规划参考曲线L1/D和L2/D为能量的二阶多项式如下：
Figure FDA00023971548300001114
Figure FDA00023971548300001115
其中，a2、a1、a0、b2、b1、b0为常值参数；
通过调节攻角和倾侧角跟踪上述参考曲线；
(1)纵程解析解
从式(77)看出，纵向射程xD是独立的，单独求解；将R*＝Re+H*代入式(69)中，可得
Figure FDA0002397154830000121
上式分母中的
Figure FDA0002397154830000122
用能量的线性函数拟合：
Figure FDA0002397154830000123
其中，参数h11、h10为常值参数，表达式如下：
Figure FDA0002397154830000124
Figure FDA0002397154830000125
其中，
Figure FDA0002397154830000126
表示终端经度；
Figure FDA0002397154830000127
表示终端能量；
此外，当飞行器严格按照大圆飞行时，式(91)分母中出现的
Figure FDA0002397154830000128
近似为常值；证明如下：
因为
Figure FDA0002397154830000129
与速度无关，所以忽略速度影响，视飞行器仅受重力作用做圆周运动；由于重力的力矩为0，令飞行器所受总力矩矢量为
Figure FDA00023971548300001210
故可得
Figure FDA00023971548300001211
在GNED坐标系下，由式(24)-(25)，可得
Figure FDA00023971548300001212
其中，
Figure FDA00023971548300001213
为力矩矢量在GNED坐标系下的表达式；XGNED为位置矢量在GNED坐标系下的表达式；VGNED为速度矢量在GNED坐标系下的表达式；
转到AGR坐标系中可得
Figure FDA00023971548300001214
所以
Figure FDA00023971548300001215
为常值；其中，
Figure FDA00023971548300001216
表示从GNED坐标系到AGR坐标系的坐标转换矩阵；
Figure FDA00023971548300001217
为力矩矢量在AGR坐标系下的表达式；
而结合式(26)，可得
Figure FDA0002397154830000131
表达式：
Figure FDA0002397154830000132
结合式(97)，易知
Figure FDA0002397154830000133
为常值；
但在实际飞行过程中，由于飞行器飞行并不严格沿着大圆飞行，
所以导致
Figure FDA0002397154830000134
并不严格为常值；针对这种情况，采用能量的线性函数拟合
Figure FDA0002397154830000135
如下：
Figure FDA0002397154830000136
其中，h21、h20为常值参数，表达式如下：
Figure FDA0002397154830000137
Figure FDA0002397154830000138
式(100)中，
Figure FDA0002397154830000139
以及求解
Figure FDA00023971548300001310
所需的
Figure FDA00023971548300001311
和
Figure FDA00023971548300001312
均由粗略估计得到；
将式(89)，(92)，(99)代入式(79)，积分可得：
Figure FDA00023971548300001313
其中，
Figure FDA00023971548300001315
Figure FDA00023971548300001316
Figure FDA00023971548300001314
c2＝(a1-c1d0)/d1 (106)
c3＝a0-c2d0 (107)
即为纵向射程解析解；
(2)横程与航向角解析解
由式(77)易知，
Figure FDA0002397154830000141
与
Figure FDA0002397154830000142
相互耦合，构成一个复杂的线性时变系统；
对式(80)-(81)联立，得
Figure FDA0002397154830000143
采用基于谱分解的LTV解析求解法，可得横程与航向角的0阶解析解：
Figure FDA0002397154830000144
Figure FDA0002397154830000145
其中，尽管上述0阶解析解中含有积分项，但由于被积曲线平缓，采用20点的Gauss-Legendre积分公式计算；
将式(82)-(83)联立，为简化式子，定义函数
Figure FDA0002397154830000146
为
Figure FDA0002397154830000147
可得
Figure FDA0002397154830000148
同样采用基于谱分解的LTV解析求解法，结合式(85)，可得横程与航向角的1阶解析解：
Figure FDA0002397154830000149
Figure FDA00023971548300001410
由于上述两式中的积分项有些复杂，故采用N阶Lagrange插值多项式拟合被积函数；
令
Figure FDA0002397154830000151
其中，i为多项式的项数；pi为多项式中第i项的系数；
Figure FDA0002397154830000152
为多项式自变量；
Figure FDA0002397154830000153
为多项式自变量的i次方；Li为第i项的Lagrange插值多项式拟合函数，定义见式(117)；
Figure FDA0002397154830000154
其中，函数
Figure FDA0002397154830000155
定义见式(111)；
Figure FDA0002397154830000156
定义为
Figure FDA0002397154830000157
其中，参数j表示从0到N且不等于i的正整数；
Figure FDA0002397154830000158
取值为
Figure FDA0002397154830000159
取N＝4时，则
Figure FDA00023971548300001510
的拟合多项式为
Figure FDA00023971548300001511
其中，p14，p13，p12，p11，p10为拟合多项式系数；
取N＝4时，误差项为
Figure FDA00023971548300001512
其中，ζ位于
Figure FDA00023971548300001513
区间，且满足积分中值定理；
同理，令
Figure FDA00023971548300001514
其中，
Figure FDA00023971548300001515
同样用4阶Lagrange插值多项式拟合得：
Figure FDA00023971548300001516
其中，p24，p23，p22，p21，p20为拟合多项式系数；
此时误差项为
Figure FDA0002397154830000161
其中，
Figure FDA0002397154830000162
将式(119)和式(123)分别代入式(113)和(114)中，可得
Figure FDA0002397154830000163
Figure FDA0002397154830000164
当将解析解应用到再入制导律中时，用这样的多项式拟合计算积分项将大大减小计算负担，提高再入制导效率；
将式(125)-(126)，式(109)-(110)代入式(87)-(88)中，可得横程、航向角解析解：
Figure FDA0002397154830000165
Figure FDA0002397154830000166
由上式看出，最终解析解由0阶解和1阶解组成；其中，1阶解中的多项式系数需要由0阶解确定。
CN202010135508.8A 2020-03-02 2020-03-02 一种旋转地球下高超声速平稳滑翔弹道的解析求解方法 Active CN111488646B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010135508.8A CN111488646B (zh) 2020-03-02 2020-03-02 一种旋转地球下高超声速平稳滑翔弹道的解析求解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010135508.8A CN111488646B (zh) 2020-03-02 2020-03-02 一种旋转地球下高超声速平稳滑翔弹道的解析求解方法

Publications (2)

Publication Number Publication Date
CN111488646A CN111488646A (zh) 2020-08-04
CN111488646B true CN111488646B (zh) 2022-04-01

Family

ID=71812455

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010135508.8A Active CN111488646B (zh) 2020-03-02 2020-03-02 一种旋转地球下高超声速平稳滑翔弹道的解析求解方法

Country Status (1)

Country Link
CN (1) CN111488646B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112861249B (zh) * 2020-12-21 2023-03-31 中国人民解放军96901部队23分队 基于攻角和阻力的滑翔弹道随速度变化降阶解计算方法
CN112861251B (zh) * 2020-12-21 2023-03-28 中国人民解放军96901部队23分队 考虑地球自转速率的滑翔弹道随速度变化降阶解计算方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9639085B1 (en) * 2015-08-05 2017-05-02 The United States Of America As Represented By The Secretary Of The Air Force Phugoid peaks trajectory for hypersonic glide vehicles
CN107941087A (zh) * 2017-10-18 2018-04-20 北京航空航天大学 一种基于阻力剖面的高升阻比高超平稳滑翔再入制导方法
CN109446582A (zh) * 2018-09-29 2019-03-08 北京航空航天大学 一种考虑地球自转的高精度降阶平稳滑翔动力学建模方法
CN109508030A (zh) * 2018-11-27 2019-03-22 北京航空航天大学 一种考虑多禁飞区约束的协同解析再入制导方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9639085B1 (en) * 2015-08-05 2017-05-02 The United States Of America As Represented By The Secretary Of The Air Force Phugoid peaks trajectory for hypersonic glide vehicles
CN107941087A (zh) * 2017-10-18 2018-04-20 北京航空航天大学 一种基于阻力剖面的高升阻比高超平稳滑翔再入制导方法
CN109446582A (zh) * 2018-09-29 2019-03-08 北京航空航天大学 一种考虑地球自转的高精度降阶平稳滑翔动力学建模方法
CN109508030A (zh) * 2018-11-27 2019-03-22 北京航空航天大学 一种考虑多禁飞区约束的协同解析再入制导方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Omnidirectionalautonomousentryguidancebasedon3-Danalytical glide formulas;WenbinYu 等;《ISATransactions》;20160919;487–503 *
高超声速飞行器平稳滑翔弹道扰动运动伴随分析;赫泰龙 等;《北京航空航天大学学报》;20190131;109-122 *

Also Published As

Publication number Publication date
CN111488646A (zh) 2020-08-04

Similar Documents

Publication Publication Date Title
CN107966156B (zh) 一种适用于运载火箭垂直回收段的制导律设计方法
Xie et al. Highly constrained entry trajectory generation
Durham Aircraft flight dynamics and control
CN111399531B (zh) 高超声速飞行器滑翔段制导与姿态控制一体化设计方法
CN107861517B (zh) 基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法
Xie et al. Rapid generation of entry trajectories with waypoint and no-fly zone constraints
Zhu et al. Impact time and angle control guidance independent of time-to-go prediction
CN105573337B (zh) 一种满足再入角和航程约束的离轨制动闭路制导方法
CN112198885B (zh) 一种满足机动平台自主降落需求的无人机控制方法
CN111488646B (zh) 一种旋转地球下高超声速平稳滑翔弹道的解析求解方法
CN110015446B (zh) 一种半解析的火星进入制导方法
CN108549785B (zh) 一种基于三维飞行剖面的高超声速飞行器精准弹道快速预测方法
CN109446582B (zh) 一种考虑地球自转的高精度降阶平稳滑翔动力学建模方法
CN109358646B (zh) 带有乘性噪声的导弹自主编队队形随机控制系统建模方法
CN109703768B (zh) 一种基于姿态/轨迹复合控制的软式空中加油对接方法
CN105865455A (zh) 一种利用gps与加速度计计算飞行器姿态角的方法
CN109703769B (zh) 一种基于预瞄策略的空中加油对接控制方法
WO2019089966A1 (en) Improved atmospheric thermal location estimation
Zhang et al. Analytical solutions to three-dimensional hypersonic gliding trajectory over rotating Earth
Desai et al. Six-degree-of-freedom trajectory optimization using a two-timescale collocation architecture
Han et al. Development of unmanned aerial vehicle (UAV) system with waypoint tracking and vision-based reconnaissance
Moutinho Modeling and nonlinear control for airship autonomous flight
CN109445283B (zh) 一种用于欠驱动浮空器在平面上定点跟踪的控制方法
Luo et al. Carrier-based aircraft precision landing using direct lift control based on incremental nonlinear dynamic inversion
CN112507467B (zh) 基于阻力和升阻比的滑翔弹道随速度变化降阶解计算方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant