CN104035335B - 基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法 - Google Patents
基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法 Download PDFInfo
- Publication number
- CN104035335B CN104035335B CN201410228163.5A CN201410228163A CN104035335B CN 104035335 B CN104035335 B CN 104035335B CN 201410228163 A CN201410228163 A CN 201410228163A CN 104035335 B CN104035335 B CN 104035335B
- Authority
- CN
- China
- Prior art keywords
- cos
- phi
- gamma
- psi
- sin
- 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
- 238000000034 method Methods 0.000 title claims abstract description 159
- 230000008569 process Effects 0.000 claims abstract description 27
- 230000007423 decrease Effects 0.000 claims abstract description 16
- 238000004088 simulation Methods 0.000 claims description 84
- 239000013598 vector Substances 0.000 claims description 36
- 238000012937 correction Methods 0.000 claims description 23
- 230000006870 function Effects 0.000 claims description 22
- 230000014509 gene expression Effects 0.000 claims description 21
- 239000011159 matrix material Substances 0.000 claims description 20
- 238000004458 analytical method Methods 0.000 claims description 17
- 230000001133 acceleration Effects 0.000 claims description 16
- 125000004432 carbon atom Chemical group C* 0.000 claims description 15
- 230000000694 effects Effects 0.000 claims description 14
- 238000004364 calculation method Methods 0.000 claims description 13
- 230000010355 oscillation Effects 0.000 claims description 12
- 230000008859 change Effects 0.000 claims description 11
- 230000008901 benefit Effects 0.000 claims description 8
- 238000009795 derivation Methods 0.000 claims description 7
- 230000009466 transformation Effects 0.000 claims description 7
- 230000005484 gravity Effects 0.000 claims description 6
- 230000010354 integration Effects 0.000 claims description 6
- 238000006243 chemical reaction Methods 0.000 claims description 5
- 238000013461 design Methods 0.000 claims description 5
- 230000003595 spectral effect Effects 0.000 claims description 5
- 241000764238 Isis Species 0.000 claims description 4
- 238000013459 approach Methods 0.000 claims description 4
- 230000000630 rising effect Effects 0.000 claims description 4
- 230000007704 transition Effects 0.000 claims description 4
- 238000000354 decomposition reaction Methods 0.000 claims description 3
- 230000003111 delayed effect Effects 0.000 claims description 3
- NCGICGYLBXGBGN-UHFFFAOYSA-N 3-morpholin-4-yl-1-oxa-3-azonia-2-azanidacyclopent-3-en-5-imine;hydrochloride Chemical compound Cl.[N-]1OC(=N)C=[N+]1N1CCOCC1 NCGICGYLBXGBGN-UHFFFAOYSA-N 0.000 claims description 2
- 150000001875 compounds Chemical class 0.000 claims description 2
- 230000000881 depressing effect Effects 0.000 claims description 2
- 230000000994 depressogenic effect Effects 0.000 claims description 2
- 230000004069 differentiation Effects 0.000 claims description 2
- 238000009499 grossing Methods 0.000 claims description 2
- 238000007726 management method Methods 0.000 claims description 2
- 230000004048 modification Effects 0.000 claims description 2
- 238000012986 modification Methods 0.000 claims description 2
- 239000002245 particle Substances 0.000 claims description 2
- 229920006395 saturated elastomer Polymers 0.000 claims description 2
- 238000001228 spectrum Methods 0.000 claims description 2
- 230000003247 decreasing effect Effects 0.000 claims 2
- RCZJVHXVCSKDKB-OYKKKHCWSA-N tert-butyl 2-[(z)-[1-(2-amino-1,3-thiazol-4-yl)-2-(1,3-benzothiazol-2-ylsulfanyl)-2-oxoethylidene]amino]oxy-2-methylpropanoate Chemical compound N=1C2=CC=CC=C2SC=1SC(=O)\C(=N/OC(C)(C)C(=O)OC(C)(C)C)C1=CSC(N)=N1 RCZJVHXVCSKDKB-OYKKKHCWSA-N 0.000 claims 2
- 230000003068 static effect Effects 0.000 claims 1
- 238000010586 diagram Methods 0.000 description 10
- RZVHIXYEVGDQDX-UHFFFAOYSA-N 9,10-anthraquinone Chemical compound C1=CC=C2C(=O)C3=CC=CC=C3C(=O)C2=C1 RZVHIXYEVGDQDX-UHFFFAOYSA-N 0.000 description 2
- 230000002829 reductive effect Effects 0.000 description 2
- 241000255925 Diptera Species 0.000 description 1
- 238000000342 Monte Carlo simulation Methods 0.000 description 1
- 241000392831 Pharyngostrongylus lambda Species 0.000 description 1
- OAICVXFJPJFONN-UHFFFAOYSA-N Phosphorus Chemical compound [P] OAICVXFJPJFONN-UHFFFAOYSA-N 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 238000013016 damping Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000002401 inhibitory effect Effects 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000002452 interceptive effect Effects 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 230000000452 restraining effect Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
- 230000007306 turnover Effects 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Landscapes
- Navigation (AREA)
Abstract
基于高精度纵、横程解析预测方法的平稳滑翔再入制导律,该方法有三大步骤:步骤1:规划快速下降段的制导方案;步骤2:规划平稳滑翔段的制导方案;步骤3:规划弹道下压段的制导方案;基于能量的纵向射程、横向射程和航向角的解析解,求解解析解的具体步骤如下:步骤1:建立飞行器再入问题的运动方程;步骤2:建立再入问题的过程约束和终端约束;步骤3:引入广义赤道和广义子午线;步骤4:运动方程线性化;步骤5:求解解析解;步骤6:验证上述解析解可行性。此制导方法能够导引飞行器飞行到距离目标一定距离的指定高度,并能够满足相应的过程约束,末速度约束和方位角误差约束,具有很好的鲁棒性。
Description
技术领域
本发明属于航天技术、高超声速飞行器再入制导技术领域,涉及一种基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法,尤其涉及一种新的基于能量-纵程及能量-横程高精度解析预测方法的通用飞行器(CAV)的平稳滑翔再入制导方法。具体给出能满足过程约束,终端的速度、高度、航向角误差和射程约束要求的再入制导方法。
背景技术
通用飞行器(CAV)是一种高超声速飞行器,它通过运载火箭或者弹道导弹助推到亚轨道,之后再入到大气层中无动力飞行。CAV的飞行过程可以大致分成两个阶段:再入滑翔段和末端制导段。CAV与运载器分离之后,开始再入滑翔段,CAV在再入过程中需要在满足热流率约束、动压约束和过载约束等过程约束的前提下,进行横向机动,一直飞行到距离目标一定距离和高度时,进入末制导段。
飞行器再入制导有两种形式:跟踪参考轨迹和预测校正制导。航天飞机再入制导是典型的参考轨迹跟踪制导,航天飞机再入飞行前先设计一个阻力-速度参考剖面来满足射程要求,其在实际飞行中通过调节倾侧角来跟踪阻力-速度参考剖面。当航向角误差超过预先设计的误差门限曲线时,通过倾侧角翻转来实现横向射程控制。另一种形式的再入制导方法是预测校正制导,是以消除实际轨道的预报落点和预定落点位置之间的偏差为目的的制导方法。关于高超声速飞行器滑翔再入问题,还有很多研究给出解析形式的结果,这些结果大多假定升阻比和弹道倾角为常数,或为某些系数的函数形式,并忽略地球曲率的影响。
发明内容
1、目的:本发明给出了一种新的基于能量-纵程及能量-横程高精度解析预测方法的通用飞行器(CAV)的平稳滑翔再入制导方法,是一种高超声速飞行器的再入制导方法。在推导过程中,通过将升阻比L/D和倾侧角设计成与能量有关的函数,并利用线性化方法得到一种特殊类型的变系数线性系统。为求解此系统,提出一种新的基于谱分解的求解方法,然后根据此求解方法得到飞行器再入飞行射程、横向射程和航向角的解析解。随后基于这些解析解,给出新的再入制导方法,此制导方法能够导引飞行器飞行到距离目标一定距离的指定高度,并能够满足相应的过程约束,末速度约束和方位角误差约束,具有很好的鲁棒性。
2、技术方案:在给出本发明的制导方法前,对于高超声速飞行器再入飞行问题,本发明给出一种新的基于能量的纵向射程、横向射程和航向角的解析解,求解解析解的具体步骤如下:
步骤1:建立飞行器再入问题的运动方程;
如附图1所示,在推导飞行器运动方程之前,首先定义两个坐标系。
地心赤道旋转(GER)坐标系:坐标系原点在地心,ze轴指向北极;xe和ye轴是赤道平面内相互垂直的两条轴线。GER坐标系以ze为中心,以与地球自转角速度相同的角速度旋转,地球自转角速度ωe为7.292116×10-5rad/s。
当地北东下(NED)坐标系:坐标系原点o为飞行器质心M和地心连线与地球表面的交点,即飞行器在地面的垂直投影点;x轴指向当地的北向,y轴指向当地的东向,z轴铅垂向下。
由于大气可以认为是相对地球静止的,选择在GER坐标系中处理再入问题较合适。将再入飞行器看作质点,在球形地球模型下,飞行器运动方程为
在上述运动方程中,如图1所示,t是时间,λ是经度,φ是纬度,H是高度,V是对旋转地球的相对速度,γ是飞行的弹道倾角,ψ是航向角,m是飞行器质量,σ是倾侧角,Re是地球的平均半径,其值为6356.766km,L和D分别为升力和阻力,g是当地的重力加速度。GER坐标系到NED坐标系的坐标转换矩阵为
步骤2:建立再入问题的过程约束和终端约束;
通常,飞行器再入飞行过程中需要满足热流率动压q和过载n的约束,约束条件为
对于本文所要研究的CAV模型,kQ=1.5×10-8,qmax=150kpa,nmax=2,其中g0为海平面重力加速度,ρ为大气密度。
当CAV飞行到距离目标水平距离STAEM=50km时,再入滑翔段结束,转入末制导段。设计的滑翔段末端条件为:航向角误差|ΔψTAEM|≤5deg,对旋转地球的相对速度为VTAEM=2000m/s,高度为HTAEM=25km,倾侧角|σTAEM|≤30deg。其中STAEM,ΔψTAEM和VTAEM是判断制导方法优劣的主要要素。下标“TAEM”代表可重复使用飞行器(RLV)中的末端能量管理段(TAEM),在此,为了方便,仍如此沿用。
步骤3:引入广义赤道和广义子午线;
为了求解出解析解,先建立辅助地心惯性(AGI)坐标系,并定义广义赤道和子午线。
在惯性空间中,目标T的位置是随着旋转地球一起旋转的,所以,如附图2所示,飞行器M飞行到预测命中点P,并不是飞行器再入飞行开始阶段的目标位置。预测命中点P可以大致预测为
在此,λP、φP和HP分别为碰撞点P的经纬度和高度。λT、φT和HT分别为目标的经纬度和高度。tgo是剩余飞行时间,可以由下式近似计算
式中,sgo是剩余射程,假设飞行器亚轨道再入点和地心的连线与地球表面交点为M,交点M与目标位置T之间的地球表面大圆弧长度即为剩余射程。由于单位向量的点乘即是两向量之间夹角的余弦值,可得
上式中,和分别为和的单位向量,为地心E指向飞行器位置M的向量,为地心E指向目标位置T的向量。上标“GER”代表其在GER坐标系中的分量。如图1所示,在ze轴的分量是sin(φ),在赤道面的分量是cos(φ),可得在xe和ye轴的分量为cos(λ)cos(φ)和sin(λ)cos(φ)。故为
上式中,上标“T”代表矩阵或者向量的转置。类似的,可以得出为
为地心E指向P的向量,的单位向量为
如图2所示,定义AGI坐标系如下:坐标系原点在地心;轴指向飞行器位置M;轴在由M、E和P组成的平面内,并垂直于轴;轴与上述两轴组成右手坐标系。AGI坐标系中轴的单位向量为
由于轴与轴垂直,运用格莱姆-施密特正交化方法,可以得到AGI坐标系的轴单位向量为
最后可以得到轴的单位向量为
在此,地心E、飞行器当前位置M和预测命中点P组成的平面为MEP平面,定义MEP平面与地球表面的交线为广义赤道平面,与广义赤道正交的大圆弧的一半定义为广义子午线,可知广义子午线端点在轴上。通过此定义,再入飞行器(CAV)的位置可以用广义经度广义纬度和广义高度表示。为分析方便,假设广义的0°子午线通过CAV初始再入点,故CAV初始位置表示为
尽管由于AGI坐标系随着CAV的运动而旋转,在此为方便分析,可认为是静止的。CAV相对于AGI坐标系的速度,用表示,是相对地球速度和牵连速度的矢量和,初始速度为
在此,上标“NED”代表向量在NED坐标系中的分量形式。定义的大小为广义速度定义和当地水平面夹角为广义弹道倾角定义的水平分量于当地广义子午线之间夹角为广义航向角由式(21)可得到,初始的和为
由于亚轨道再入初始点M点沿广义子午线上的切线与平行,因此与当地水平分量之间夹角即为初始航向角为
上式中的和是由公式(18)和(19)得到的,由如下公式得到
上式中,是由公式(7)给出。由图2可知,P点的广义经纬度和高度为
前面说明了对旋转地球模型下的再入制导的末端条件,这些条件需要转移到AGI坐标系中,定义为
通过运用格莱姆-施密特正交化方法,得到与正交的向量为
定义x3为与过P点与水平面、广义赤道交线都平行的单位向量。x3在NEDP坐标系中可表示为
上式中,缩写词“NEDP”代表在位置P处的NED坐标系,将λ=λP和φ=φP带入到公式(7)中,即得到从GER坐标系到NEDP坐标系的坐标转换矩阵
定义相对于AGI坐标系的最终速度向量为在AGI坐标系中,再入制导方法导引飞行器近似沿着广义赤道到达P,并且使得广义弹道倾角接近0,在此为方便分析,如图3所示,假设与x3平行。尽管这样假设不准确,会造成一定的误差,但是制导方法也能够通过最后一次倾侧角翻转来消除其引起的误差。定义相对于旋转地球的飞行器最终速度为VTAEM。除此之外,为得到地球旋转引起的牵连速度假设滑翔段结束距离P的距离STAEM远小于地球半径。指向当地东向,大小为
速度向量的大小为在以上前提情况下,可得
由图3可得
将上式展开得
上式中,是在y轴上的分量。进而可得到的估计值为
定义相对于惯性空间AGI坐标系的单位质量机械能为可表示为
式中,μ是重力常数,则最终的末端约束条件转化成相对于惯性空间的能量为
步骤4:运动方程线性化;
在分析时,AGI坐标系可以认为是静止的,得到运动方程为
其中能量对时间的导数为
定义射程xD和横向射程xC为和将公式(39,40,44)与公式(45)联立得到
令L1=Lcos(σ),L2=Lsin(σ)。假设飞行器以平稳滑翔条件(Steady GlideCondition,SGC)滑翔飞行,飞行器的广义弹道倾角变化率为0,即由公式(43)可得
经过转换可得
由式(37)可得
在AGI坐标系中,由于制导方法导引飞行器近似沿着广义赤道飞向P,可知 因此可以假设和tan(xC/Re)=xC/Re,将式(50)(51)带入到使(46-48),运用上述假设,可得
为减小微分计算的复杂性,式(48)中的分母可假设为1。
步骤5:求解解析解;
由于 对式(52-54)的解影响很小。故在分析时,可以设为一常量令R*=Re+H*。定义纵向升阻比为L1/D=Lcos(σ)/D,将L1/D写成如下形式
定义水平升阻比为L2/D=Lsin(σ)/D,将L2/D写成如下形式
在此,设计L1/D和L2/D曲线,可以通过同时调节攻角和倾侧角来跟踪参考升阻比曲线。
1.射程表达式
将式(55)带入式(52),积分可得射程的表达式为
2.横向射程和方位角表达式
将式(53)和(54)结合,得到
令
式(58)可写成如下形式
下面我们引入谱分解方法来求解上述特殊类型的变系数线性系统。定义为
在式(64)两边乘以得到
上式(66)可以写成
为求取其微分,反向利用点乘法则,可得
对式(68)两边积分,得
注意到其中02×2是2×2阶零矩阵,I2×2为2×2单位阵。的逆为
将式(70)左乘以式(69)得
上式中,为状态转移矩阵。通过对矩阵A谱分解,由矩阵论知识可得
上式中,λ1=i和λ2=-i是矩阵A的特征值,其中G1和G2是矩阵A的谱投影,为
进而可得
和
由式(60)和(62),有
这样就可以求得,横向射程和航向角的解析公式为
在此,需要提到的是在求解再入制导方法时,飞行器飞行的攻角为方案攻角,通过调节倾侧角来跟踪规划的L1/D曲线,这意味着L2/D曲线不能任意规划,也不能像式(56)一样用有限次数多项式表示。然而由于L2/D的表达式只包含并没有在上述微分中体现,因此上述两个解析解仍然适用。
步骤6:验证上述解析解可行性;
下面给出一个解析解和弹道仿真的对比实例。仿真中使用的是CAV模型,初始条件为终端条件为仿真中考虑两种情形,在第一种情形中,CAV再入滑翔段距离较长,设计的参考曲线为
L1/D=2.5 (81)
在第二种情形中,CAV再入滑翔段较短,设计的参考曲线为
L1/D=1.5 (83)
文章“A Closed-Form Solution to Lifting Reentry”中,Bell在假设L/D和倾侧角σ为常数的前提下,求解出射程,横向射程和航向角的解析解形式,其方法忽略了地球曲率对航向角产生的影响,即式(43)中二次项。在此为作对比,Bell方法的解析解分别表示为和为了便于区别,将本文中给出的解析解用Yu-Chen法(YCGF)表示,Bell的解析解用Bell法(BGF)表示。根据上述参考剖面,可将BGF修改为如下形式
其中,
上式中,通过弹道仿真可算得在时,大约为5802.29m/s。弹道仿真最终的末速度为在第一种情形中,ksgn=1,在第二种情形中,ksgn=-1。在计算解析解积分时采用20个点的高斯-勒让得积分法。在弹道仿真中,采用了消除长周期振动的弹道阻尼技术。除此之外,为了能够更好的跟踪参考曲线,仿真中忽略了对倾侧角变化速率的限制。仿真程序用C语言编写,在装有因特尔酷睿-2双核处理器的计算机上运行得到下面的结果。
仿真结果如附图4-6所示,由下表1中可以看出,YCGF具有很好的精度,而其运算速度比弹道仿真小两个数量级。
首先,我们分析射程方向结果。由于在做微分运算时假设这就意味着忽略了航向角误差,由图中看出,YGCF法的射程结果比弹道仿真的结果要大。而在实际运用时,再入制导方法可以通过多次倾侧角翻转来减小航向角误差,提高YCGF法射程的准确性。在这里也给出了BGF法的仿真结果作为对比。由于BGF法在推导射程解析式时考虑了航向角误差,BGF法在射程上的精度要比YCGF法要高,在第二种情形中比较明显。然而,相对于YCGF法,在在线实时规划L1/D剖面时BGF法需要更大的运算载荷,因为在线规划曲线是需要考虑射程和航向角的交互影响,而且需要数值方法计算积分。
接下来分析横向射程和航向角的结果。从图表的结果中可以看出来BGF法在横向射程和航向角上的精度要比YCGF方法低,尤其是在射程较大时显得更加明显。这是由于BGF法在做微分计算时,省略了式(44)括弧中二次项。除此之外,在YCGF法中射程表达式中含有一个积分项,由于YCGF法使得倾侧角翻转只需做几次更新,而不是实时规划,所以运算速度更快。
表1 弹道仿真与解析解的结果对比
本发明在给出上述新的高超声速飞行器再入飞行问题射程、横向射程和航向角的解析解后,根据这些解析解本发明给出一种新的再入制导方法。首先,根据飞行器再入飞行的弹道特点,可以将整个再入飞行过程分成三个阶段:快速下降段(Descent Phase,DP),平稳滑翔段(Steady Glide Phase,SGP),弹道下压段(Altitude Adjustment Phase,AAP)。制导方法分别给出了每段的相应的制导方法。具体阐述如下:
步骤1:规划快速下降段的制导方案;
CAV在与运载器分离之后进入无动力滑翔状态,开始进入快速下降段,直到开始满足平稳滑翔条件(SGC)时结束。在这一段,由于大气密度ρ非常小,飞行器快速下降,高度快速降低。随着高度下降,大气密度升高,飞行器的热流率急剧升高,而最大热流率大致出现在这段的结束。为了能满足热流率约束,设计此段飞行器飞行时保持最大攻角(αmax)并保持倾侧角为0,这样可以使得在下降段飞行的更高,从而降低热流率。当时,说明此时开始升力已经足够大,可以满足平稳滑翔条件。为保证飞行器能从DP平滑转移到SGP,下式(88)给出攻角方案。指令攻角和倾侧角为
σcmd=σplan=0deg (89)
Δγ=γplan-γ (90)
上式中,αplan是SGP的规划的攻角,在式(91)中给出。γplan是平稳滑翔的弹道倾角,将在后面给出说明。γ为当前飞行器的弹道倾角,Δγ1是当时的Δγ,kγ是值为5的常数。由式(88)可知,Δγ从Δγ1逐渐趋近于0的同时,指令攻角αcmd从αmax平滑趋向于αplan。当Δγ=0,时,下降段DP结束,因为之后飞行器将向上飞。在此处飞行器相对于惯性空间的单位质量接机械能表示为
步骤2:规划平稳滑翔段的制导方案;
在平稳滑翔段,飞行器沿着规划的攻角曲线飞行。制导方法利用射程方向的解析表达式实时更新L1/D曲线,并通过调节倾侧角跟踪纵向升阻比曲线。制导方法利用横向射程的表达式来确定倾侧角的第一次翻转节点。下面给出详细结果。
1.规划的攻角(αplan)和升阻比(L/D)plan
(L/D)plan是方案攻角(αplan)对应的升阻比曲线。L/D通常随着攻角和马赫数(Ma)变化。Ma受大气运动影响,但是在此可以认为大气是相对旋转地球静止的。定义相对于旋转地球的单位质量机械能E,如式(91)所示,在此,将方案攻角αplan和(L/D)plan表示为能量E的函数形式。这样做的优势在于不同飞行任务的αplan和(L/D)plan不需要再进行调整修改。
为了能够飞出更远的射程,飞行器在平稳滑翔段大部分时间内以α1=10deg飞行,这样可以近似保持最大的升阻比(L/D)max。当飞行器飞行到轨迹末段时,平滑的将αplan减小到α2=6deg,可以将高度H降低到HTAEM进入末制导段,这样可以保证飞行器有足够的动压满足末端机动飞行。在此设计αplan为
上式中E0是初始的能量,平稳滑翔飞行结束时的E2可设为-5.6×107J/kg,Ef是由终端约束所决定的能量,可以通过将VTAEM和HTAEM带入式(91)中得到。
当αplan曲线确定后,就确定了高度走廊,如图7所示。高度走廊的下边界(Hmin)由前面所述的过程约束决定。除此之外忽略地球旋转的影响,假设倾侧角σ=0,这样可以通过求解式(5)得到高度走廊的上边界(Hmax)。
下面给出相对能量E的升阻比(L/D)plan曲线。由于αplan曲线已经确定,(L/D)plan是马赫数Ma的函数,而Ma是速度V和高度H的函数。所以,如果确定了高度曲线,就可以得到相应的(L/D)plan曲线。实际上,当飞行器高速滑翔时,升阻比L/D很轻微的随Ma变化,但是由于高度曲线被限定在很小的走廊中,H对Ma影响很小,所以高度对(L/D)plan曲线影响很小。在此,可以假设飞行器沿着高度走廊的中间飞行,高度为
当给定E和H时,可以得到V和Ma,然后利用αplan和Ma,可以得到相对于E的(L/D)plan曲线,如图8所示。从图中可以看出,最大的升阻比可达3。
在利用解析解快速生成参考弹道的过程中,为了确定第一次翻转节点需要知道在惯性空间中的基准升阻比曲线这就需要建立E与之间的转换关系。由式(22,37,91)中可知,E与之间的关系为
由于在实际飞行中,飞行任务剩下的弹道曲线并不明确,故上式中的关系并不实用。通常E0=-3.5×104kJ/kg,Ef=-6×104kJ/kg,V0ωe(Re+H0)≈3.3×103kJ/kg,Vfωe(Re+Hf)≈0.93×103kJ/kg和ωe 2(Re+H0)2≈0.22×103kJ/kg。从中可以看出式(94)中的非线性项远小于线性项,即线性项占主要作用,因此在实际仿真结果中与E看起来近似直线。这里与E之间近似采用如下线性转化关系
上式中,E与为当前相对旋转地球和相对于惯性空间的能量值,Ef与(式(37))为相对旋转地球和相对于惯性空间的终端能量值。xE和任意时刻相对旋转地球和相对于惯性空间的能量值。由于(L/D)plan曲线一大部分可认为是常值,由式(95)确立的误差较小,并且误差可以通过第二次倾侧角翻转来修正。在此,令与相关的(L/D)plan曲线为可以由下式得到
在下式(97)中,当指令攻角αplan曲线和高度曲线确立之后,可以将其转化为最大倾侧角σmax约束。当飞行器以平稳滑翔条件飞行时,增大σ会使得高度H降低。当将H降低到Hmin时,即可得到最大的σmax。如式(98)所示,在平稳滑翔时,纵向升力L1与重力和向心力近似平衡,由式(97)中右边第一项,可以求得σmax。若飞行轨迹有振荡,式(97)右边第二项可以消除振荡。因此,若高度H快速下降,到达(H)E-(Hmin)E>0,式(97)右边第二项会降低σmax,使得飞行器能够拉起,保持在Hmin之上。在此,为书写方便,函数f对E求导,写成(f)E。
Lmax=CLplan[0.5ρ(Hmin)V2]Sref(99)
上式中,kσ可设为-50。CLplan为αplan的升力系数,由式(3,4,91)可知,E的变化率为
忽略地球旋转影响,由式(3)和式(100),可得
2.规划纵向升阻比L1/D曲线
再入制导方法利用射程上的解析解实时规划纵向升阻比曲线,然后通过调节倾侧角跟踪规划的升阻比曲线,以满足终端速度约束。如图8所示,当E>E2时,(L/D)plan可以认为是常量。故时,设计纵向升阻比L1/D曲线为水平段。当时,(L/D)plan随着E下降而下降,导致L1/D不能认为是常数。但是,在此可认为(L1/D)1到(L1/D)2的转变是线性的,L1/D曲线设计为
上式中,是将xE=E2带入到式(95)中得到的。(L1/D)1和(L1/D)2是需要实时修正的参量。由于倾侧角σ的余弦为纵向升阻比与总升阻比之比,而且最终的倾侧角约束要求接近0。(L1/D)2则由下式实时计算
其中,(L/D)est(E)是由惯性导航系统(Inertial Navigation System,INS)得到的,(L/D)est(E)是在当前能量下的规划升阻比,(L/D)plan(Ef)是在终点能量下的规划升阻比。在AGI坐标系中,剩余射程为
上式中,由式(26)求取,下面计算(L1/D)1:
a.当
将式(102)带入式(52),然后积分,得
其中当时有
当时有
其中,将带入式(107)可得通过求解式(105),得
其中
b.当
在这种情形下,由于大部分时间飞行器处于AAP阶段,不需要跟踪纵向升阻比L1/D曲线,因此其不需要再更新。
3.确定第一次倾侧角翻转节点
再入制导方法利用横向射程的解析解表达式(式79),来决定第一次倾侧角翻转的能量节点倾侧角翻转是规划的,规划的横向升阻比L2/D曲线为
上式中,是第二次倾侧角翻转的能量节点。在此假设剩余飞行中的L/D微分由INS给出,预测的为
sgn的符号由再入初始条件决定
从上面的结果可以看出,这里安排两次倾侧角翻转,第一次在第二次在并且第二次倾侧角翻转离终点状态较近,这样选取接近末端的有下列好处:的调节可以修正之前积累的误差,并且使得后面剩余较短射程内的误差较小,近似可以忽略。为方便利用横向射程解,构造如下积分函数
上述积分函数与横向射程解析解积分形式类似。由于纵向升阻比L1/D曲线是分段函数,也可以写成如下分段函数形式
其中,和由式(107)求得,可由式(106)求得。由于假设最终的横向射程解析式仅是的函数。将初始条件式(20-23)和L2/D曲线式(111-113)带入横向射程解析式(79)得到最终形式为
其导数为
在此,运用牛顿法求解的解
为了提升运算效率,做三次运算。第一次运算是在平稳滑翔段的开始,得到当T1=200s,时,第二次计算,得到当T2=30s,时,第三次计算得到在这里,是因为倾侧角变化率的限制所加的延迟,在式(121)中给出。aD是由惯导系统所测得的当前的阻力加速度。
4.确定基准倾侧角
由于L1=Lcos(σ),L1/D=L/Dcos(σ),为跟踪纵向升阻比L1/D曲线,基准的倾侧角为
其中,是考虑到飞行器的滚转速率的限制,提前时间Δt翻转基准倾侧角,其值为
上式中,Δt完成倾侧角翻转所需时间的一半,可以用下式估计
其中,σ是当前倾侧角,是最大滚转速率。
5.攻角指令αcmd和倾侧角指令σcmd
若将基准攻角αplan和倾侧角σplan直接用作攻角和倾侧角指令,再入飞行轨迹将会有微弱的长周期振动(Phugoid Oscillation)。在此,下面给出一种消除长周期振动的方法。如图9所示,升力系数的垂直分量CL1和水平分量CL2,分别由基准攻角αplan和倾侧角σplan决定。为消除振动,在垂直方向上附加一个与垂直速度分量相反的升力,即在CL1上加上ΔCL1,设计ΔCL1为
其中,是升力线斜率。γplan为平稳滑翔所需要的近似弹道倾角,在后面式(132)中给出。kγ为增益系数,值为5。若飞行器上升较快,γplan-γ<0,ΔCL1<0,造成升力L1降低,从而抑制飞行器上升过快,同样的可以抑制下降过快。
a.指令攻角αcmd推导
对于给定的马赫数Ma,攻角是升力系数CL的函数,α=fα(CL),其一阶Taylor展开为
注意到有
将式(123,125)代入到式(124)中,得
αcmd=αplan+cos(σplan)kγ(γplan-γ) (126)
b.指令倾侧角σcmd推导
由图9可知
其一阶Taylor近似为
其中,对于CAV,令αplan=10deg。为了保证飞行中满足各个过程约束,倾侧角有最大界限为
其中,σmax为最大倾侧角约束,由式(97)得到。
c.基准弹道倾角γplan推导
γplan为平稳滑翔的近似弹道倾角。由式(5),假设忽略地球自转影响,可得
将上式对时间求导得
其中,由式(100)可得
为简化式(131),假设γplan和g的导数为0。除此之外,γplan≈0,令cos(γplan)=1,sin(γplan)=γplan。将式(3-4)带入式(131)中,忽略地球自转影响,可以得到
上式中,Dplan=CDplanqSref,CDplan是攻角为αplan时的阻力系数,ρ是大气密度,Sref为参考面积,d1和d2分别表示如下
前面提到,将函数f对E的导数定为(f)E。因为CLplan(E)已经规划,(CLplan)E是已知的。在平稳滑翔段可用式(136)求得[cos(σplan)]E。然而,在弹道下压段,由于[cos(σplan)]E的计算公式较复杂,可用相邻两时刻的cos(σplan)的差分计算。
其中,(L/D)plan(E)是规划曲线,((L/D)plan)E是已知的,(L1/D)E可由式(95)和(102)得到。
步骤3:规划弹道下压段的制导方案;
飞行器经过第二次倾侧角翻转之后,进入弹道下压段(AAP)。通过降低攻角,飞行器快速进入稠密大气层,这样能够使得飞行器在末制导段有较大的动压进行机动。这就使得在这段中,不能近似接近0,所以使得上述解析解难以满足终端约束条件。在此用一种与前面不同的制导方案:在最后一次倾侧角翻转之前,为得到需求的最终速度,在线的弹道仿真修正能量节点在最后一次倾侧角翻转之后,采用比例导引法(ProportionalNavigation,PN)导引,下面首先介绍弹道下压段的制导方案,再给出修正的方法。
1.基准倾侧角
在弹道下压段,αplan曲线与平稳滑翔段相同,而倾侧角由比例导引律确定。当飞行器接近目标时,可以忽略地球曲率的影响。如图2所示,MP视线的方位角的变化率为
其中,由比例导引律得到的需用横向机动加速度为
其中,为了防止初始倾侧角饱和,令kPN从2逐渐变化到4,为
其中,xDf是第二次倾侧角翻转的能量节点。另外,升力的垂直分量需要与重力和向心力平衡,可得
进而可以得到基准倾侧角为
其中,由式(121)求得。
2.指令攻角和指令倾侧角
在飞行器进入弹道下压段之前,需要获得射程-能量参考曲线并且需要修正在弹道下压段,通过调节攻角跟踪曲线,如式(142)所示,可以消除由于干扰引起的最终速度误差。注意到在弹道下压段中,大部分时间升阻比L/D对攻角导数大于0,意味着增加攻角会增加升阻比,造成最终速度升高。为了消除长周期振动,σcmd与式(128)一样。
其中,kα值为(5π/18)×10-6,也就是说0.5deg的攻角修正,会引起10km的射程散布。在弹道下压段中,在计算式(133)中的γplan时,由于[cos(σplan)]E表达式较复杂,其对E的导数由差分求得。同时,为了满足过程约束,仍然需要最大倾侧角约束
3.第二次倾侧角翻转
在弹道下压段,不需要跟踪L1/D曲线,需要调整翻转能量节点来满足最终速度要求。下面分析与最终速度Vf之间的关系:若降低则翻转时间推迟,时刻的航向角误差增加,为消除航向角误差,需要更大的横向机动加速度,因而倾侧角的大小会增大,造成纵向升阻比L1/D减小,导致最终速度降低。所以,降低造成最终速度Vf降低。
在第一次倾侧角翻转之后,在此采用预测校正法修正通过弹道仿真预测最终速度,之后运用secant法修正当前的状态作为仿真的初始状态,由于INS能够实时测量当前气动加速度,通过对比实际和理想气动加速度,可以得到当前气动系数的偏差。弹道预测仿真假设当前估计的偏差即为剩余飞行段的气动系数偏差。除此之外,仿真中忽略式(142)中右边的第二项。当Sgo=STAEM时,仿真结束,可以得到预测的最终速度修正的secant法式为
其中,为第k次仿真的为减小运算次数,的修正有两次:第一次预测校正过程在第一次倾侧角翻转之后,得到T1=60s时,此时进行第二次预测校正过程,得到之后可令
当CAV飞行到Sgo=STAEM时,仿真结束。
3、优点及功效:本发明的优点在于:
(1)本发明给出新的关于能量的射程、横向射程和航向角的解析解,求取解析解时考虑了地球曲率的影响,精度更高。
(2)相较于传统的跟踪参考轨迹和预测校正制导,本发明通过新的解析解规划再入制导方法,能很好的避免对整个飞行过程中的重复积分,降低运算载荷,也便于工程实现。
(3)本发明将再入飞行过程分为三段,分别规划每段的再入制导方法。为了减小飞行控制系统(FCS)的负荷,再入飞行设计成只有两侧倾侧角翻转,相比与传统方法的多次翻转,有很大的提升。
(4)由于典型的再入飞行轨迹有一定的长周期振动,本发明给出三维的抑制长周期振动的方法,使得飞行轨迹更加平稳。
(5)由于参考剖面和倾侧角翻转是在线规划的,本发明给出的制导方法能够很好的在一系列气动和大气参数散布情况下,保持很好的鲁棒性。
附图说明
图1是GER坐标系:E-xeyeze和NED坐标系o-xyz的示意图。
图2是辅助地心惯性(AGI)坐标系的示意图。
图3是在AF坐标系下的终点速度近似计算示意图。
图4是解析解和弹道仿真结果中的射程对比结果图。
图5是解析解和弹道仿真结果中的横向射程对比结果图。
图6是解析解和弹道仿真结果中的航向角对比结果图。
图7是规划的高度走廊图。
图8是基准升阻比曲线图。
图9是附加升力系数抑制弹道振动的示意图。
图10是制导方案流程框图。
图11是三种理想情形仿真结果的高度-速度曲线图。
图12是三种理想情形仿真结果的地面投影曲线图。
图13是三种理想情形仿真结果的弹道倾角曲线图。
图14是三种理想情形仿真结果的攻角曲线图。
图15是三种理想情形仿真结果的倾侧角曲线图。
图16是气动拉偏情形下的升阻比拉偏百分比曲线图。
图17是气动拉偏情形下仿真结果的高度-速度曲线图。
图18是气动拉偏情形下仿真结果的地面投影曲线图。
图19是气动拉偏情形下仿真结果的攻角曲线图。
图20是气动拉偏情形下仿真结果的倾侧角曲线图。
图21是YCEG法Monto Carlo仿真的高度-速度曲线图。
图22是LGG法Monto Carlo仿真的高度-速度曲线图。
图23是YCEG法Monto Carlo仿真的地面投影曲线图。
图24是LGG法Monto Carlo仿真的地面投影曲线图。
图25是YCEG法Monto Carlo仿真的攻角曲线图。
图26是LGG法Monto Carlo仿真的攻角曲线图。
图27是YCEG法Monto Carlo仿真的倾侧角曲线图。
图28是LGG法Monto Carlo仿真的倾侧角曲线图。
图29是YCEG法Monto Carlo仿真的热流率曲线图。
图30是LGG法Monto Carlo仿真的热流率曲线图。
图31是YCEG法Monto Carlo仿真的最终状态图。
图32是LGG法Monto Carlo仿真的最终状态图。
图中符号说明如下:
图1中:xe、ye和ze轴分别为地心赤道旋转(GER)坐标系的三个坐标轴,x、y和z轴分别为当地北东下(NED)坐标系的三个坐标轴,E为地心,M是飞行器质心,o为飞行器在地面投影,λ是经度,φ是纬度,H是高度,V是对旋转地球的相对速度,VH为V的水平分量,γ是飞行的弹道倾角,ψ是航向角,ωe为地球自转角速度。
图2中:E为地心,M是飞行器质心,T为飞行器再入时目标位置,P为预测命中点。和轴分别为辅助地心惯性(AGI)坐标系三个坐标轴。
图3中:P为预测命中点,xP和yP轴分别为P点NED坐标系的坐标轴,为最终速度向量,VTAEM为相对于旋转地球模型的飞行器最终速度,为地球旋转引起的牵连速度,θP为与yP之间夹角,x3为与过P点与水平面和广义赤道交线平行的单位向量。
图4-6中:横坐标的为相对惯性空间的单位质量机械能。
图7-8中:横坐标的E为相对旋转地球的单位质量机械能,图7中的Hmin为高度走廊的下边界,Hmax为高度走廊的上边界,图8中的L/Dplan为参考的升阻比。
图9中:M为飞行器质心,CL1和CL2分别为升力系数的垂直分量和水平分量。ΔCL1为抑制长周期振动附加的升力系数垂直分量。
图13中:γ为弹道倾角,γplan为平稳滑翔所需要的近似弹道倾角。
图11、17、21和22中的Hmin为高度走廊的下边界。
具体实施方式
下面具体结合附图以及仿真实例对本发明做进一步详细的说明。
图1是GER坐标系:E-xeyeze和NED坐标系o-xyz的示意图,形象显示飞行器再入模型的相关位置和角度关系。图2是辅助地心惯性(AGI)坐标系的示意图,显示AGI系的定义。图3是在AF坐标系下的终点速度近似计算示意图,帮助求解旋转地球模型下的相关速度关系。图4-6分别是解析解和弹道仿真结果中的纵向射程对比、横向射程以及航向角对比结果图,可以很清楚显示本文解析解的可行性。图7是规划的高度走廊图,是由相关的过程约束所确定的。图8是基准升阻比曲线图,对于不同的飞行任务,基准升阻比曲线近似不变。图9是附加升力系数抑制弹道振动的示意图。为形象说明本发明的流程,图10给出了本发明的制导方案流程框图。
下面具体结合几个仿真案例来展示该制导方法的可行性及优势。
本发明用于再入制导仿真的模型是CAV,其最大升阻比可达3。由于飞行控制系统的限制,约束攻角的变化率为:和倾侧角约束为:|σ|≤80deg、和飞行器的初始条件为:H0=80km,λ0=0rad,φ0=0rad,V0=6500m/s,γ0=0rad,
案例1:无干扰理想情形
首先观察三种理想的情形,不考虑干扰的影响,目标之间距离较远。情形1,目标位置为λT1Re=7500km,φT1Re=1000km;情形2,λT2Re=10000km,φT2Re=0;情形3,λT3Re=12000km,φT3Re=2000km。仿真结果如图11-15所示,较小射程情形下,为跟踪纵向升阻比L1/D曲线,会降低倾侧角σ的大小,从而降低高度H。图11是三种理想情形仿真结果的高度-速度曲线图,可以看出高度速度曲线始终保持在下边界之上,能满足相关的过程约束。图12是三种理想情形仿真结果的地面投影曲线图,显示飞行器能够准确飞行到距离目标一定距离上空,并近似指向目标,航向角误差很小。图13是三种理想情形仿真结果的弹道倾角曲线图,显示γplan是对γ的很好的近似。图14是三种理想情形仿真结果的攻角曲线图,为抑制长周期振动,攻角并没有保持规划的值,而是有很小的反馈修正。图15是三种理想情形仿真结果的倾侧角曲线图。仿真结果还显示,消除长周期振动的方法能够很好的消除由于倾侧角翻转造成的轨迹振荡。
表2 三种理想情形的仿真结果
案例2:气动拉偏情形
在实际飞行中,存在气动参数散布。一般情况下,气动拉偏系数是攻角和马赫数的函数。为简化考虑,下面给出线性的气动力拉偏系数模型
上式中,δCL是升力系数CL拉偏系数,δCD是阻力系数CD拉偏系数。由于L/D对射程有重要影响,下面给出升阻比L/D较大偏差两种情形。这里仅考虑阻力系数拉偏的情况,即δCL=0。情形1,令δCD0=-20%,和情形2,令δCD0=25%,和目标位置为λT1Re=9000km,φT1Re=1000km。气动拉偏结果如附图16所示,仿真结果如附图17-20所示。图17是气动拉偏情形下仿真结果的高度-速度曲线图,高度速度曲线始终保持在下边界之上。图18是气动拉偏情形下仿真结果的地面投影曲线图,显示飞行器能够准确飞行到距离目标一定距离上空,并近似指向目标,航向角误差很小。图19是气动拉偏情形下仿真结果的攻角曲线图。图20是气动拉偏情形下仿真结果的倾侧角曲线图。仿真结果显示若CD增加,导致L/D降低,倾侧角σ大小降低,从而降低横向机动能力,为满足横向射程要求,第一次翻转时间推迟。经过第二次倾侧角翻转之后,为跟踪剩余射程-能量曲线,需要轻微的调整α。
表3 气动拉偏结果
案例3:Monte Carle仿真情形
表4 正态分布的随机变量数值特性
Monte Carlo仿真可以检验再入制导方法的鲁棒性。如表4所示,给出了再入初始条件,气动模型和大气模型的散布。式(146-147)给出了气动拉偏模型。δρ代表大气密度的百分比拉偏,代表东西向的风速,代表南北向的风速。在每次仿真中,δρ、和为常数。尽管散布模型较为简单,并不与实际情形相符合,但是这些情形比实际情况更为严格,能够很好的验证制导方法的鲁棒性。作为对比,同时给出文献“ConstrainedPredictor-Corrector Entry Guidance”中的制导方法Monte Carlo仿真结果。为区分两种方法,本文中的制导方法称为Yu-Chen’s Entry Guidance(YCEG),文献中的制导方法称为Lu’s Gliding Guidance(LGG)。为满足热流率约束,对LGG做轻微修改,在飞行器快速下降段,攻角最大,倾侧角为0。除此之外,LGG的攻角曲线也是通过式(92)给出的,由于LGG最终的倾侧角趋近于30deg,为达到HTAEM,将α2设为7deg。多次仿真绘制仿真结果如图21-32所示。图21和22分别是YCEG法和LGG法Monto Carlo仿真的高度-速度曲线图。图23和24分别是YCEG法和LGG法Monto Carlo仿真的地面投影曲线图,可以看出两种方法都能使得飞行器准确飞行到距离目标一定距离上空,并近似指向目标,航向角误差较小。图25和26分别是YCEG法和LGG法Monto Carlo仿真的攻角曲线图,YCEG法中攻角的轻微的反馈修正能够很好的为抑制飞行轨迹的长周期振动。图27和28分别是YCEG法和LGG法Monto Carlo仿真的倾侧角曲线图。从仿真结果可以看出,YCEG能够很有效的消除由倾侧角翻转造成的振荡。在LGG中,由于消除振荡的增益随时间下降,在飞行末端,振荡不能够很好的抑制。在YCEG中,整个飞行过程只有两次倾侧角翻转,而且位置比较集中。但是在LGG中,需要由3-5次倾侧角翻转,而且由于L/D的散布,次数和位置不确定,更大的L/D会增加倾侧角翻转的次数。图29和30分别为YCEG法和LGG法Monto Carlo仿真的热流率曲线图,热流率都没有超过约束值。图31和32分别是YCEG法和LGG法Monto Carlo仿真的最终状态图。从图31-32中,可以看出来,YCEG的最终航向角误差远小于LGG。尽管有几个仿真结果,可能是由于大气密度ρ的下降,造成出现H稍小于Hmin的情况,但是过程约束仍然能满足。
Claims (1)
1.基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法,实施该制导方法前需要进行求解解析,具体步骤如下:
步骤1:建立飞行器再入问题的运动方程;
首先定义两个坐标系,
地心赤道旋转即GER坐标系:坐标系原点在地心,ze轴指向北极;xe和ye轴是赤道平面内相互垂直的两条轴线;GER坐标系以ze为中心,以与地球自转角速度相同的角速度旋转,地球自转角速度ωe为7.292116×10-5rad/s;
当地北东下即NED坐标系:坐标系原点o为飞行器质心M和地心连线与地球表面的交点,即飞行器在地面的垂直投影点;x轴指向当地的北向,y轴指向当地的东向,z轴铅垂向下;
由于大气认为是相对地球静止的,选择在GER坐标系中处理再入问题较合适;将再入飞行器看作质点,在球形地球模型下,飞行器运动方程为
在上述运动方程中,t是时间,λ是经度,φ是纬度,H是高度,V是对旋转地球的相对速度,γ是飞行的弹道倾角,Ψ是航向角,m是飞行器质量,σ是倾侧角,Re是地球的平均半径,其值为6356.766km,L和D分别为升力和阻力,g是当地的重力加速度,GER坐标系到NED坐标系的坐标转换矩阵为
步骤2:建立再入问题的过程约束和终端约束;
通常,飞行器再入飞行过程中需要满足热流率动压q和过载n的约束,约束条件为
对于CAV模型,kQ=1.5×10-8,qmax=150kpa,nmax=2,其中g0为海平面重力加速度,ρ为大气密度;
当CAV飞行到距离目标水平距离STAEM=50km时,再入滑翔段结束,转入末制导段;设计的滑翔段末端条件为:航向角误差|ΔψTAEM|≤5deg,对旋转地球的相对速度为VTAEM=2000m/s,高度为HTAEM=25km,倾侧角|σTAEM|≤30deg;其中,STAEM,ΔΨTAEM和VTAEM是判断制导方法优劣的主要要素;下标“TAEM”代表可重复使用飞行器RLV中的末端能量管理段TAEM;
步骤3:引入广义赤道和广义子午线;
为了求解出解析解,先建立辅助地心惯性即AGI坐标系,并定义广义赤道和子午线;
在惯性空间中,目标T的位置是随着旋转地球一起旋转的,飞行器M飞行到预测命中点P,并不是飞行器再入飞行开始阶段的目标位置,预测命中点P大致预测为
在此,λP、φP和HP分别为碰撞点P的经纬度和高度;λT、φT和HT分别为目标的经纬度和高度,tgo是剩余飞行时间,由下式近似计算
式中,sgo是剩余射程,假设飞行器亚轨道再入点和地心的连线与地球表面交点为M,交点M与目标位置T之间的地球表面大圆弧长度即为剩余射程;由于单位向量的点乘即是两向量之间夹角的余弦值,得
上式中,和分别为和的单位向量,为地心E指向飞行器位置M的向量,为地心E指向目标位置T的向量;上标“GER”代表其在GER坐标系中的分量,在ze轴的分量是sin(φ),在赤道面的分量是cos(φ),可得在xe和ye轴的分量为cos(λ)cos(φ)和sin(λ)cos(φ);故为
上式中,上标“T”代表矩阵或者向量的转置,类似的,可以得出为
为地心E指向P的向量,的单位向量为
定义AGI坐标系如下:坐标系原点在地心;轴指向飞行器位置M;轴在由M、E和P组成的平面内,并垂直于轴;轴与上述两轴组成右手坐标;AGI坐标系中轴的单位向量为
由于轴与轴垂直,运用格莱姆-施密特正交化方法,得到AGI坐标系的轴单位向量为
最后可以得到轴的单位向量为
在此,地心E、飞行器当前位置M和预测命中点P组成的平面为MEP平面,定义MEP平面与地球表面的交线为广义赤道平面,与广义赤道正交的大圆弧的一半定义为广义子午线,可知广义子午线端点在轴上;通过此定义,再入飞行器CAV的位置用广义经度广义纬度和广义高度表示;为分析方便,假设广义的0°子午线通过CAV初始再入点,故CAV初始位置表示为
尽管由于AGI坐标系随着CAV的运动而旋转,在此为方便分析,可认为是静止的;CAV相对于AGI坐标系的速度,用表示,是相对地球速度和牵连速度的矢量和,初始速度为
在此,上标“NED”代表向量在NED坐标系中的分量形式;定义的大小为广义速度定义和当地水平面夹角为广义弹道倾角定义的水平分量于当地广义子午线之间夹角为广义航向角由式(21)得到,初始的和为
由于亚轨道再入初始点M点沿广义子午线上的切线与平行,因此与当地水平分量之间夹角即为初始航向角为
上式中的和是由公式(18)和(19)得到的,由如下公式得到
上式中,是由公式(7)给出P点的广义经纬度和高度为
前面说明了对旋转地球模型下的再入制导的末端条件,这些条件需要转移到AGI坐标系中,定义为
通过运用格莱姆-施密特正交化方法,得到与正交的向量为
定义x3为与过P点与水平面、广义赤道交线都平行的单位向量,x3在NEDP坐标系中表示为
上式中,缩写词“NEDP”代表在位置P处的NED坐标系,将λ=λP和φ=φP带入到公式(7)中,即得到从GER坐标系到NEDP坐标系的坐标转换矩阵
定义相对于AGI坐标系的最终速度向量为在AGI坐标系中,再入制导方法导引飞行器近似沿着广义赤道到达P,并且使得广义弹道倾角接近0,在此为方便分析,假设与x3平行,制导律也能够通过最后一次倾侧角翻转来消除其引起的误差;定义相对于旋转地球的飞行器最终速度为VTAEM,除此之外,为得到地球旋转引起的牵连速度假设滑翔段结束距离P的距离STAEM远小于地球半径,指向当地东向,大小为
速度向量的大小为在以上前提情况下,可得
将上式展开得
上式中,是在y轴上的分量,进而得到的估计值为
定义相对于惯性空间AGI坐标系的单位质量机械能为表示为
式中,μ是重力常数,则最终的末端约束条件转化成相对于惯性空间的能量为
步骤4:运动方程线性化;
在分析时,AGI坐标系认为是静止的,得到运动方程为
其中能量对时间的导数为
定义射程xD和横向射程xC为和将公式(39)(40)(44)与公式(45)联立得到
令L1=Lcos(σ),L2=Lsin(σ);假设飞行器以平稳滑翔条件即Steady GlideCondition,SGC滑翔飞行,飞行器的广义弹道倾角变化率为0,即由公式(43)得
经过转换可得
由式(37)可得
在AGI坐标系中,由于制导律导引飞行器近似沿着广义赤道飞向P,可知 因此假设和tan(xC/Re)=xC/Re,将式(50)、(51)带入到式(46)、(47)和(48),运用上述假设,得
为减小微分计算的复杂性,式(48)中的分母假设为1;
步骤5:求解解析解;
由于 对式(52)(53)和(54)的解影响很小,故在分析时,设为一常量令R*=Re+H*;定义纵向升阻比为L1/D=Lcos(σ)/D,将L1/D写成如下形式
定义水平升阻比为L2/D=Lsin(σ)/D,将L2/D写成如下形式
在此,设计L1/D和L2/D曲线,通过同时调节攻角和倾侧角来跟踪参考升阻比曲线;
①射程表达式
将式(55)带入式(52),积分可得射程的表达式为
②横向射程和方位角表达式
将式(53)和(54)结合,得到
令
式(58)写成如下形式
引入谱分解方法来求解上述特殊类型的变系数线性系统,定义为
在式(64)两边乘以得到
上式(66)写成
为求取其微分,反向利用点乘法则,可得
对式(68)两边积分,得
注意到其中02×2是2×2阶零矩阵,I2×2为2×2单位阵,的逆为
将式(70)左乘以式(69)得
上式中,为状态转移矩阵,通过对矩阵A谱分解,由矩阵论知识得
上式中,λ1=i和λ2=-i是矩阵A的特征值,其中G1和G2是矩阵A的谱投影,为
进而可得
和
由式(60)和(62),有
这样就求得,横向射程和航向角的解析公式为
在此,需要提到的是在求解再入制导方法时,飞行器飞行的攻角为方案攻角,通过调节倾侧角来跟踪规划的L1/D曲线,这意味着L2/D曲线不能任意规划,也不能像式(56)一样用有限次数多项式表示;然而由于L2/D的表达式只包含并没有在上述微分中体现,因此上述两个解析解仍然适用;
步骤6:验证上述解析解可行性;
该基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法,具体步骤如下:
步骤A:规划快速下降段的制导方案;
CAV在与运载器分离之后进入无动力滑翔状态,开始进入快速下降段,直到开始满足平稳滑翔条件即SGC时结束;在这一段,由于大气密度ρ非常小,飞行器快速下降,高度快速降低;随着高度下降,大气密度升高,飞行器的热流率急剧升高,而最大热流率大致出现在这段的结束;为了能满足热流率约束,设计此段飞行器飞行时保持最大攻角(αmax)并保持倾侧角为0,这样使得在下降段飞行的更高,从而降低热流率;当时,说明此时开始升力已经足够大,满足平稳滑翔条件;为保证飞行器能从DP平滑转移到SGP,下式(81)给出攻角方案,指令攻角和倾侧角为
σcmd=σplan=0deg (82)
Δγ=γplan-γ (83)
上式中,αplan是SGP的规划的攻角,在式(84)中给出;γplan是平稳滑翔的弹道倾角,γ为当前飞行器的弹道倾角,Δγ1是当时的Δγ,kγ是值为5的常数;由式(81)可知,Δγ从Δγ1逐渐趋近于0的同时,指令攻角αcmd从αmax平滑趋向于αplan;当Δγ=0,时,下降段DP结束,因为之后飞行器将向上飞,在此处飞行器相对于惯性空间的单位质量接机械能表示为
步骤B:规划平稳滑翔段的制导方案;
在平稳滑翔段,飞行器沿着规划的攻角曲线飞行;制导律利用射程方向的解析表达式实时更新L1/D曲线,并通过调节倾侧角跟踪纵向升阻比曲线,制导律利用横向射程的表达式来确定倾侧角的第一次翻转节点,下面给出详细结果;
(1)规划的攻角(αplan)和升阻比(L/D)plan
(L/D)plan是方案攻角(αplan)对应的升阻比曲线,(L/D)通常随着攻角和马赫数(Ma)变化,Ma受大气运动影响,在此认为大气是相对旋转地球静止的;定义相对于旋转地球的单位质量机械能E,如式(84)所示,在此,将方案攻角αplan和(L/D)plan表示为能量E的函数形式;这样做的优势在于不同飞行任务的αplan和(L/D)plan不需要再进行调整修改;
为了能够飞出更远的射程,飞行器在平稳滑翔段大部分时间内以α1=10deg飞行,这样近似保持最大的升阻比(L/D)max;当飞行器飞行到轨迹末段时,平滑的将αplan减小到α2=6deg,将高度H降低到HTAEM进入末制导段,这样保证飞行器有足够的动压满足末端机动飞行,在此设计αplan为
上式中E0是初始的能量,平稳滑翔飞行结束时的E2可设为-5.6×107J/kg,Ef是由终端约束所决定的能量,通过将VTAEM和HTAEM带入式(84)中得到;
当αplan曲线确定后,就确定了高度走廊,高度走廊的下边界(Hmin)由前面所述的过程约束决定;除此之外忽略地球旋转的影响,假设倾侧角σ=0,这样通过求解式(5)得到高度走廊的上边界(Hmax);
下面给出相对能量E的升阻比(L/D)plan曲线,由于αplan曲线已经确定,(L/D)plan是马赫数Ma的函数,而Ma是速度V和高度H的函数;所以,确定了高度曲线,就得到相应的(L/D)plan曲线;实际上,当飞行器高速滑翔时,升阻比(L/D)很轻微的随Ma变化,但是由于高度曲线被限定在很小的走廊中,H对Ma影响很小,所以高度对(L/D)plan曲线影响很小;在此,假设飞行器沿着高度走廊的中间飞行,高度为
当给定E和H时,得到V和Ma,然后利用αplan和Ma,得到相对于E的(L/D)plan曲线,最大的升阻比可达3;
在利用解析解快速生成参考弹道的过程中,为了确定第一次翻转节点需要知道在惯性空间中的基准升阻比曲线这就需要建立E与之间的转换关系,由式22、37和91可知,E与之间的关系为
由于在实际飞行中,飞行任务剩下的弹道曲线并不明确,故上式中的关系并不实用,通常
从中看出式(87)中的非线性项远小于线性项,即线性项占主要作用,因此在实际仿真结果中与E看起来近似直线;这里与E之间近似采用如下线性转化关系
上式中,E与为当前相对旋转地球和相对于惯性空间的能量值,Ef与即式(37)为相对旋转地球和相对于惯性空间的终端能量值;xE和任意时刻相对旋转地球和相对于惯性空间的能量值,由于(L/D)plan曲线一大部分是常值,由式(89)确立的误差较小,并且误差通过第二次倾侧角翻转来修正,在此,令与相关的(L/D)plan曲线为由下式得到
在下式(91)中,当指令攻角αplan曲线和高度曲线确立之后,将其转化为最大倾侧角σmax约束;当飞行器以平稳滑翔条件飞行时,增大σ会使得高度H降低,当将H降低到Hmin时,即得到最大的σmax;如式(92)所示,在平稳滑翔时,纵向升力L1与重力和向心力近似平衡,由式(91)中右边第一项,求得σmax;若飞行轨迹有振荡,式(91)右边第二项可以消除振荡,因此,若高度H快速下降,到达(H)E-(Hmin)E>0,式(91)右边第二项会降低σmax,使得飞行器能够拉起,保持在Hmin之上;在此,函数f对E求导,写成fE;
Lmax=CLplan[0.5ρ(Hmin)V2]Sref (93)
上式中,kσ设为-50,CLplan为αplan的升力系数,由式(3)(4)和(84)可知,E的变化率为
忽略地球旋转影响,由式(3)和式(94),得
(2)规划纵向升阻比L1/D曲线
再入制导方法利用射程上的解析解实时规划纵向升阻比曲线,然后通过调节倾侧角跟踪规划的升阻比曲线,以满足终端速度约束;当E>E2时,(L/D)plan认为是常量,故时,设计纵向升阻比L1/D曲线为水平段;当时,(L/D)plan随着E下降而下降,导致L1/D不能认为是常数;但是,在此可认为(L1/D)1到(L1/D)2的转变是线性的,L1/D曲线设计为
上式中,是将xE=E2带入到式(89)中得到的,(L1/D)1和(L1/D)2是需要实时修正的参量,由于倾侧角σ的余弦为纵向升阻比与总升阻比之比,而且最终的倾侧角约束要求接近0,(L1/D)2则由下式实时计算
其中,(L/D)est(E)是由惯性导航系统即Inertial Navigation System,INS得到的,(L/D)est(E)是在当前能量下的规划升阻比,(L/D)plan(Ef)是在终点能量下的规划升阻比;在AGI坐标系中,剩余射程为
上式中,由式(26)求取,下面计算(L1/D)1:
a.当
将式(96)带入式(52),然后积分,得
其中当时有
当时有
其中,将带入式(101)得通过求解式(99),得
其中
b.当
在这种情形下,由于大部分时间飞行器处于AAP阶段,不需要跟踪纵向升阻比L1/D曲线,因此其不需要再更新;
(3)确定第一次倾侧角翻转节点
再入制导方法利用横向射程的解析解表达式(79),来决定第一次倾侧角翻转的能量节点倾侧角翻转是规划的,规划的横向升阻比L2/D曲线为
上式中,是第二次倾侧角翻转的能量节点,在此假设剩余飞行中的L/D微分由INS给出,预测的为
sgn的符号由再入初始条件决定
从上面的结果看出,这里安排两次倾侧角翻转,第一次在第二次在并且第二次倾侧角翻转离终点状态较近,这样选取接近末端的有下列好处:的调节修正之前积累的误差,并且使得后面剩余较短射程内的误差较小,近似可以忽略;为方便利用横向射程解,构造如下积分函数
上述积分函数与横向射程解析解积分形式类似,由于纵向升阻比L1/D曲线是分段函数,写成如下分段函数形式
其中,和由式(101)求得,由式(100)求得;由于假设最终的横向射程解析式仅是的函数,将初始条件式(20)、(21)、(22)、(23)和L2/D曲线式(105)、(106)、(107)带入横向射程解析式(79)得到最终形式为
其导数为
在此,运用牛顿法求解的解
为了提升运算效率,做三次运算;第一次运算是在平稳滑翔段的开始,得到当T1=200s,时,第二次计算,得到当T2=30s,时,第三次计算得到在这里,是因为倾侧角变化率的限制所加的延迟,在式(115)中给出,aD是由惯导系统所测得的当前的阻力加速度;
(4)确定基准倾侧角
由于L1=Lcosσ,L1/D=L/Dcosσ,为跟踪纵向升阻比L1/D曲线,基准的倾侧角为
其中,是考虑到飞行器的滚转速率的限制,提前时间Δt翻转基准倾侧角,其值为
上式中,Δt完成倾侧角翻转所需时间的一半,用下式估计
其中,σ是当前倾侧角,是最大滚转速率;
(5)攻角指令αcmd和倾侧角指令σcmd
若将基准攻角αplan和倾侧角σplan直接用作攻角和倾侧角指令,再入飞行轨迹将会有微弱的长周期振动,在此,下面给出一种消除长周期振动的方法;升力系数的垂直分量CL1和水平分量CL2,分别由基准攻角αplan和倾侧角σplan决定;为消除振动,在垂直方向上附加一个与垂直速度分量相反的升力,即在CL1上加上ΔCL1,设计ΔCL1为
其中,是升力线斜率,γplan为平稳滑翔所需要的近似弹道倾角,在后面式(126)中给出;kγ为增益系数,值为5,若飞行器上升较快,γplan-γ<0,ΔCL1<0,造成升力L1降低,从而抑制飞行器上升过快,同样的可以抑制下降过快;
c.指令攻角αcmd推导
对于给定的马赫数Ma,攻角是升力系数CL的函数,α=fα(CL),其一阶Taylor展开为
注意到有
将式(117)和(119)代入到式(118)中,得
αcmd=αplan+cos(σplan)kγ(γplan-γ) (120)
d.指令倾侧角σcmd推导
其一阶Taylor近似为
其中,对于CAV,令αplan=10deg;为了保证飞行中满足各个过程约束,倾侧角有最大界限为
其中,σmax为最大倾侧角约束,由式(91)得到;
e.基准弹道倾角γplan推导
γplan为平稳滑翔的近似弹道倾角,由式(5),假设忽略地球自转影响,得
将上式对时间求导得
其中,由式(94)得
为简化式(125),假设γplan和g的导数为0,除此之外,γplan≈0,令cos(γplan)=1,sin(γplan)=γplan,将式(3)和(4)带入式(125)中,忽略地球自转影响,得到
上式中,Dplan=CDplanqSref,CDplan是攻角为αplan时的阻力系数,ρ是大气密度,Sref为参考面积,d1和d2分别表示如下
前面提到,将函数f对E的导数定为fE;因为CLplan(E)已经规划,(CLplan)E是已知的;在平稳滑翔段用式(130)求得[cos(σplan)]E;然而,在弹道下压段,由于[cos(σplan)]E的计算公式较复杂,用相邻两时刻的cos(σplan)的差分计算;
其中,(L/D)plan(E)是规划曲线,[(L/D)plan]E是已知的,(L1/D)E可由式(89)和(96)得到;
步骤C:规划弹道下压段的制导方案;
飞行器经过第二次倾侧角翻转之后,进入弹道下压段即AAP;通过降低攻角,飞行器快速进入稠密大气层,这样能够使得飞行器在末制导段有较大的动压进行机动;这就使得在这段中,不能近似接近0,所以使得上述解析解难以满足终端约束条件;在此用一种与前面不同的制导方案:在最后一次倾侧角翻转之前,为得到需求的最终速度,在线的弹道仿真修正能量节点在最后一次倾侧角翻转之后,采用比例导引法即ProportionalNavigation,PN导引,下面首先介绍弹道下压段的制导方案,再给出修正的方法;
①基准倾侧角
在弹道下压段,αplan曲线与平稳滑翔段相同,而倾侧角由比例导引律确定;当飞行器接近目标时,忽略地球曲率的影响,MP视线的方位角的变化率为
其中,由比例导引律得到的需用横向机动加速度为
其中,为了防止初始倾侧角饱和,令kPN从2逐渐变化到4,为
其中,xDf是第二次倾侧角翻转的能量节点;另外,升力的垂直分量需要与重力和向心力平衡,可得
进而可以得到基准倾侧角为
其中,由式(115)求得;
②指令攻角和指令倾侧角
在飞行器进入弹道下压段之前,需要获得射程-能量参考曲线并且需要修正在弹道下压段,通过调节攻角跟踪曲线,如式(136)所示,消除由于干扰引起的最终速度误差;注意到在弹道下压段中,大部分时间升阻比L/D对攻角导数大于0,意味着增加攻角会增加升阻比,造成最终速度升高,为了消除长周期振动,σcmd与式(122)一样;
其中,kα值为(5π/18)×10-6,也就是说0.5deg的攻角修正,会引起10km的射程散布;在弹道下压段中,在计算式(127)中的γplan时,由于[cos(σplan)]E表达式较复杂,其对E的导数由差分求得;同时,为了满足过程约束,仍然需要最大倾侧角约束
③第二次倾侧角翻转
在弹道下压段,不需要跟踪L1/D曲线,需要调整翻转能量节点来满足最终速度要求;下面分析与最终速度Vf之间的关系:若降低则翻转时间推迟,时刻的航向角误差增加,为消除航向角误差,需要更大的横向机动加速度,因而倾侧角的大小会增大,造成纵向升阻比L1/D减小,导致最终速度降低;所以,降低造成最终速度Vf降低;
在第一次倾侧角翻转之后,在此采用预测校正法修正通过弹道仿真预测最终速度,之后运用secant法修正当前的状态作为仿真的初始状态,由于INS能够实时测量当前气动加速度,通过对比实际和理想气动加速度,得到当前气动系数的偏差;弹道预测仿真假设当前估计的偏差即为剩余飞行段的气动系数偏差,除此之外,仿真中忽略式(136)中右边的第二项,当Sgo=STAEM时,仿真结束,得到预测的最终速度修正的secant法式为
其中,为第k次仿真的为减小运算次数,的修正有两次:第一次预测校正过程在第一次倾侧角翻转之后,得到T1=60s时,此时进行第二次预测校正过程,得到之后可令
当CAV飞行到Sgo=STAEM时,仿真结束。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410228163.5A CN104035335B (zh) | 2014-05-27 | 2014-05-27 | 基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410228163.5A CN104035335B (zh) | 2014-05-27 | 2014-05-27 | 基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104035335A CN104035335A (zh) | 2014-09-10 |
CN104035335B true CN104035335B (zh) | 2017-01-04 |
Family
ID=51466147
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201410228163.5A Active CN104035335B (zh) | 2014-05-27 | 2014-05-27 | 基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104035335B (zh) |
Families Citing this family (41)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104392047B (zh) * | 2014-11-25 | 2017-05-10 | 北京航空航天大学 | 一种基于平稳滑翔弹道解析解的快速弹道规划方法 |
CN104376225A (zh) * | 2014-11-27 | 2015-02-25 | 江西洪都航空工业集团有限责任公司 | 一种风标式攻角传感器的攻角修正计算方法 |
CN104634183B (zh) * | 2014-12-18 | 2016-01-13 | 北京控制工程研究所 | 一种基于升阻比实时估计的自适应制导方法 |
CN104714553B (zh) * | 2015-01-14 | 2017-03-29 | 西北工业大学 | 基于几何规划的滑翔飞行器末端能量管理轨迹规划方法 |
CN104656450B (zh) * | 2015-01-20 | 2017-04-19 | 北京航空航天大学 | 一种高超声速飞行器平稳滑翔再入弹道设计方法 |
CN104787361B (zh) * | 2015-04-02 | 2015-11-04 | 北京航天自动控制研究所 | 升力式飞行器再入制导的阻力加速度变化率的确定方法 |
CN105022403B (zh) * | 2015-05-11 | 2016-03-23 | 北京航天自动控制研究所 | 滑翔飞行器的纵向轨迹控制增益的确定方法 |
CN104965418B (zh) * | 2015-06-01 | 2017-05-10 | 北京航空航天大学 | 一种基于弹道阻尼控制和热流解析预测的引入段制导方法 |
CN106021628B (zh) * | 2015-07-03 | 2019-06-18 | 中国运载火箭技术研究院 | 一种运载火箭垂直返回弹道设计方法 |
CN105115512B (zh) * | 2015-09-23 | 2017-10-10 | 北京理工大学 | 一种火星大气进入段侧向预测校正制导方法 |
CN105718660B (zh) * | 2016-01-21 | 2019-03-01 | 中国工程物理研究院总体工程研究所 | 临近空间大范围机动弹道三维包络计算方法 |
CN105740506B (zh) * | 2016-01-21 | 2018-12-11 | 中国工程物理研究院总体工程研究所 | 沿临近空间大范围机动弹道空间包络的扰动引力逼近方法 |
CN105760573B (zh) * | 2016-01-21 | 2019-07-02 | 中国工程物理研究院总体工程研究所 | 沿临近空间大范围机动弹道的扰动引力延拓逼近方法 |
CN105923172B (zh) * | 2016-04-18 | 2017-03-22 | 北京航天自动控制研究所 | 一种升力式飞行器的倾斜转弯翻转制导方法 |
CN105759830B (zh) * | 2016-04-29 | 2017-03-22 | 北京航天自动控制研究所 | 一种升力式飞行器高动态下压段制导方法 |
CN106020216B (zh) * | 2016-05-13 | 2017-03-22 | 北京航天自动控制研究所 | 一种攻角约束下的平衡滑翔制导力分配方法 |
CN106292701B (zh) * | 2016-08-16 | 2018-12-21 | 北京控制工程研究所 | 一种基于扰动补偿思想的rlv进场着陆段制导律获取方法 |
CN106371312B (zh) * | 2016-09-12 | 2019-09-20 | 中国人民解放军国防科学技术大学 | 基于模糊控制器的升力式再入预测-校正制导方法 |
CN106643341B (zh) * | 2017-02-24 | 2018-06-01 | 北京临近空间飞行器系统工程研究所 | 基于准平衡滑翔原理的力热控制耦合设计方法 |
CN107065933B (zh) * | 2017-04-19 | 2020-04-21 | 中国人民解放军海军航空大学 | 基于气动模型的临近空间高超声速目标跟踪方法 |
CN107132765B (zh) * | 2017-06-01 | 2020-04-28 | 烟台南山学院 | 一种基于轨迹规划的攻击角度与攻击时间控制方法 |
CN107121015B (zh) * | 2017-06-16 | 2018-10-16 | 湖北航天技术研究院总体设计所 | 一种快速弹上弹道在线规划方法 |
CN107168374B (zh) * | 2017-07-06 | 2020-07-21 | 中国人民解放军军械工程学院 | 横向平面的自适应比例微分导引方法 |
CN107844128B (zh) * | 2017-10-13 | 2018-11-16 | 北京航空航天大学 | 一种基于复合比例导引的高超声速飞行器巡航段制导方法 |
CN108036676B (zh) * | 2017-12-04 | 2019-08-23 | 北京航空航天大学 | 一种基于三维再入弹道解析解的全射向自主再入制导方法 |
CN108387140B (zh) * | 2018-01-19 | 2020-05-05 | 北京航空航天大学 | 一种考虑多个禁飞区约束的解析再入制导方法 |
CN108548541B (zh) * | 2018-03-13 | 2020-09-18 | 北京控制工程研究所 | 一种以开伞高度为控制目标的大气进入制导方法 |
CN108549785B (zh) * | 2018-05-03 | 2021-09-24 | 中国人民解放军国防科技大学 | 一种基于三维飞行剖面的高超声速飞行器精准弹道快速预测方法 |
CN109508030B (zh) * | 2018-11-27 | 2020-08-04 | 北京航空航天大学 | 一种考虑多禁飞区约束的协同解析再入制导方法 |
CN111290427B (zh) * | 2018-12-06 | 2021-07-09 | 北京理工大学 | 抗高过载的飞行器侧偏修正系统 |
CN109941460B (zh) * | 2019-04-09 | 2020-08-07 | 北京空间技术研制试验中心 | 航天器亚轨道返回再入过载降低设计方法 |
CN110147521B (zh) * | 2019-04-25 | 2021-02-02 | 北京航空航天大学 | 一种高超声速飞行器跳跃滑翔弹道解析求解方法 |
CN110425943B (zh) * | 2019-08-06 | 2021-05-07 | 西北工业大学 | 面向变质心飞行器的工程化再入制导方法 |
CN110908407B (zh) * | 2019-11-15 | 2021-06-22 | 南京航空航天大学 | 一种rlv再入热流率跟踪的改进预测制导方法 |
CN110989644B (zh) * | 2019-11-29 | 2021-04-23 | 上海宇航系统工程研究所 | 一种考虑目标点多终端约束的飞行器轨迹规划方法 |
CN111306989B (zh) * | 2020-03-12 | 2021-12-28 | 北京航空航天大学 | 一种基于平稳滑翔弹道解析解的高超声速再入制导方法 |
CN112069605B (zh) * | 2020-11-10 | 2021-01-29 | 中国人民解放军国防科技大学 | 一种带有攻击时间约束的比例导引律设计方法 |
CN114167888B (zh) * | 2021-11-19 | 2023-06-20 | 湖北航天技术研究院总体设计所 | 一种滑翔高超声速飞行器末端位置和速度控制方法 |
CN116361926B (zh) * | 2023-06-01 | 2023-08-18 | 西安现代控制技术研究所 | 一种弹道导弹滑翔增程段初始机械能闭环调整方法 |
CN117647933B (zh) * | 2024-01-26 | 2024-03-29 | 中国人民解放军国防科技大学 | 一种面向精度提升的轨迹规划方法 |
CN118313054A (zh) * | 2024-03-19 | 2024-07-09 | 上海华模科技有限公司 | 气动模型的确定方法及装置、飞行模拟机及存储介质 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
GB1435203A (en) * | 1972-07-14 | 1976-05-12 | Sperry Rand Corp | Aircraft automatic flight control system |
EP0530924A2 (en) * | 1991-09-03 | 1993-03-10 | The Boeing Company | Aircraft flare control system utilizing an envelope limiter |
EP0708393A1 (en) * | 1994-10-19 | 1996-04-24 | Honeywell Inc. | Self-adaptive limiter for automatic control of approach and landing |
CN103728976A (zh) * | 2013-12-30 | 2014-04-16 | 北京航空航天大学 | 一种基于广义标控脱靶量概念的多过程约束和多终端约束末制导律 |
-
2014
- 2014-05-27 CN CN201410228163.5A patent/CN104035335B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
GB1435203A (en) * | 1972-07-14 | 1976-05-12 | Sperry Rand Corp | Aircraft automatic flight control system |
EP0530924A2 (en) * | 1991-09-03 | 1993-03-10 | The Boeing Company | Aircraft flare control system utilizing an envelope limiter |
EP0708393A1 (en) * | 1994-10-19 | 1996-04-24 | Honeywell Inc. | Self-adaptive limiter for automatic control of approach and landing |
CN103728976A (zh) * | 2013-12-30 | 2014-04-16 | 北京航空航天大学 | 一种基于广义标控脱靶量概念的多过程约束和多终端约束末制导律 |
Non-Patent Citations (6)
Title |
---|
Bézier approximation based inverse dynamic guidance for entry glide trajectory;Tawfiqur Rahman,et al.;《Control Conference (ASCC)》;20130131;第1-6页 * |
Entry Guidance with Real Time Planning of Reference based on Analytical Solutions;Wenbin Yu,et al.;《Advances in Space Research》;20150228;第1-21页 * |
Quasi-equilibrium glide trajectory design of waverider-based hypersonic vehicle;Ping Li,et al.;《2010 International Conference on Mechanical and Electrical Technology (ICMET 2010)》;20100930;第408-412页 * |
高超声速飞行器再入多段导引方法研究;刘冠南等;《飞行力学》;20120831;第30卷(第4期);第337-344页 * |
高超声速飞行器多约束多种机动突防模式弹道规划;张科南等;《弹道学报》;20120930;第24卷(第3期);第85-90页 * |
高超声速飞行器滑行航迹优化;周浩等;《北京航空航天大学学报》;20060531;第32卷(第5期);第513-517页 * |
Also Published As
Publication number | Publication date |
---|---|
CN104035335A (zh) | 2014-09-10 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104035335B (zh) | 基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法 | |
CN108036676B (zh) | 一种基于三维再入弹道解析解的全射向自主再入制导方法 | |
CN111306989B (zh) | 一种基于平稳滑翔弹道解析解的高超声速再入制导方法 | |
Li et al. | Review and prospect of guidance and control for Mars atmospheric entry | |
Brunner et al. | Skip entry trajectory planning and guidance | |
Xie et al. | Highly constrained entry trajectory generation | |
Saraf et al. | Design and evaluation of an acceleration guidance algorithm for entry | |
Li et al. | Improved artificial potential field based lateral entry guidance for waypoints passage and no-fly zones avoidance | |
CN107941087A (zh) | 一种基于阻力剖面的高升阻比高超平稳滑翔再入制导方法 | |
CN110015446B (zh) | 一种半解析的火星进入制导方法 | |
CN106383994A (zh) | 一种基于脉冲和气动辅助结合的低轨轨道面转移方法 | |
CN105573337A (zh) | 一种满足再入角和航程约束的离轨制动闭路制导方法 | |
CN106054604A (zh) | 基于模型预测控制理论的再入飞行器鲁棒最优制导方法 | |
Zhang et al. | Entry guidance for high-L/D hypersonic vehicle based on drag-vs-energy profile | |
Liang et al. | Waypoint constrained guidance for entry vehicles | |
Furfaro et al. | Mars atmospheric entry guidance via multiple sliding surface guidance for reference trajectory tracking | |
Figueroa-González et al. | Landing guidance strategy of space rider | |
Cihan et al. | Analytical Earth-aerocapture guidance with near-optimal performance | |
CN107102547A (zh) | 一种基于滑模控制理论的rlv着陆段制导律获取方法 | |
Bojun et al. | High-precision adaptive predictive entry guidance for vertical rocket landing | |
CN109977543A (zh) | 基于侧向优先的三维剖面覆盖区域计算方法、系统及介质 | |
CN102508492A (zh) | 一种飞行器在等高航路点间的定高度大圆飞行实现方法 | |
CN107796401B (zh) | 跳跃式再入飞行器线性伪谱参数修正横向制导方法 | |
Zhou et al. | A simple reentry trajectory generation and tracking scheme for common aero vehicle | |
Yomchinda et al. | Autonomous control and path planning for autorotation of unmanned helicopters |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant |