汪 雷,刘 欣,杨 涛,张梦樱
(国防科学技术大学 航天科学与工程学院,长沙410073)
高超声速滑翔式再入飞行器外形通常采用升力体或乘波体,可实现高升阻比再入,利用空气动力控制飞行轨迹,能够完成远距离的非弹道式再入机动飞行。高超声速滑翔式飞行器以其在增大射程、突破导弹防御系统和再入段机动等方面具备的强大优势而备受瞩目,成为近年来的研究热点,美国的HTV-2是其典型代表。高超声速滑翔式飞行器再入后通过气动力控制弹道,可实现大范围的机动飞行。本文将飞行器可通过机动飞行到达的区域称为目标覆盖范围,国外的文献[1]中一般将这个区域称为“footprint”,此区域可初步反映飞行器的射程能力和机动能力。另外,当飞行中飞行器出现故障或飞行任务变更需改变落点时,必须知道当前弹道参数下飞行器的目标覆盖范围,以便选择可行的落点。Lu Ping等人[2]采用伪平衡滑翔条件,研究了再入飞行过程中的目标覆盖范围的快速生成。
本文首先介绍了高超声速滑翔式飞行器的再入动力学模型和再入过程中受到的弹道约束;然后给出了利用弹道优化算法求解目标覆盖范围的基本流程,在合理选择优化指标与约束条件的情况下,选用伪谱法来完成优化计算。鉴于上述方法计算量过大,本文提出了一种基于高度-速度(h-v)剖面跟踪的目标覆盖范围求解方法,通过对不同h-v曲线的跟踪和侧向机动得到完整的射程覆盖范围。以典型高超声速滑翔式再入飞行器为对象,验证了所提出方法的有效性。
由于高超声速滑翔式飞行器外形为升力体,弹体扁平,常采用倾斜(BTT)转弯技术,故可假设飞行中侧滑角保持为零。假设地球为不旋转圆球,建立三自由度运动模型,其中位置参数以地心距r、经度λ和纬度φ3个参数来描述;速度参数由速度大小v、速度倾角θ和速度偏角σ确定。速度倾角θ是速度向量与当地水平面的夹角,速度向量指向水平面上方时θ为正。速度偏角σ是速度向量在当地水平面的投影与正北方向的夹角,从正北方向到速度向量顺时针旋转时σ为正;ν为倾侧角。无动力三自由度再入运动方程为[3]
式中:m为飞行器质量;飞行器高度h=r-R0,R0为平均地球半径;重力加速度g=μ/r2,μ=3.986 005×1014m3/s2为 地 球 引 力 系 数分别为飞行器的气动阻力和升力;Cx,Cy分别为阻力系数、升力系数;ρ为大气密度;S为飞行器参考面积。
飞行器再入过程是一个非常复杂的飞行过程,必须考虑热流、动压、过载及机动能力等因素对飞行器的影响,即这些因素不能超过飞行器的最大承受能力。
驻点热流约束:
式中:ks为取决于飞行器头部形状的热流传递系数。
法向过载约束:
式中:ny为法向过载,α为攻角。
动压约束:
式中:q为动压,qmax为最大动压。
另外,基于再入制导控制能力的考虑,为保证弹道的可控性,必须保证沿弹道飞行时飞行器可获得的最大升力始终能够平衡其它力,此约束被称之为平衡滑翔约束,可取约束条件为
定义再入走廊为再入飞行器安全返回所必须满足的各种约束条件的交集,再入走廊可直观地反映飞行器再入过程约束。取指数大气密度模型:ρ(H)=ρ0e-h/hs,其中,ρ0为海平面大气密度,hs=7.11km。根据驻点热流、动压、法向过载约束以及平衡滑翔约束,可得h-v剖面内与各约束条件对应的再入走廊边界:
式中:Cy,max表示对应速度下的最大升力系数。
代入各约束的具体取值,可计算得再入走廊,如图1所示。
此外,弹道约束还包括控制量约束和终端弹道参数约束等。作为对地打击的再入飞行器,根据作战需要,通常还对落地速度和落地弹道倾角有要求,即
式中:θf为落地的最小弹道倾角,vf为落地最小速度,tf表示弹道终端时刻。
如进行目标覆盖能力的分析,则终端位置约束只考虑r(tf)=rf。
图1 高度-速度剖面再入走廊
目标覆盖范围可通过计算一系列最优再入弹道确定。弹道优化的算法可采用目前飞行器再入轨迹优化领域较为流行的伪谱法。
弹道优化问题可描述为一般的最优控制问题,即在时间区间[t0,tf]中寻找最优控制变量u(t)∈Rl,最小化性能指标:
且使状态变量x(t)∈Rn、初始时间t0、终端时间tf满足运动微分方程约束:
以及边界条件(端点约束):
和过程约束:
弹道优化的算法可采用目前飞行器再入轨迹优化领域较为流行的伪谱法。伪谱法的基本原理是依据Legendre-Gauss(LG)等离散点对时间轴进行离散化处理,并利用Lagrange多项式作为基函数来近似状态变量与控制变量,然后引入导数矩阵将动力学方程转化为等式约束,最终将最优控制问题转化为非线性规划问题求解,其基本原理可见文献[4]。目前已出现一些基于伪谱法的优化软件包,本文采用基于Matlab的Gpops软件包[5]完成优化计算。
目标覆盖范围的优化求解流程如下。
①计算“最偏”的2条弹道。
分别将优化指标取为J=min[±σB(tf)],其中,σB(tf)为再入点与弹道终点的连线与正北方向的夹角,其取值可根据再入点与弹道终点的经纬度计算:
将优化所得的σB(tf)最大值、最小值分别记为σB,fmax和σB,fmin,引入新的弹道终端约束条件:
式中:w∈[0,1],为权重系数。
②计算目标覆盖范围的远边界。
加入式(13)所示的新的弹道终端约束,其余弹道约束不变,取不同的w值分别进行最大射程优化,即优化指标取为J=min(-L(tf))(L表示射程),可得到一组优化弹道,其落点将组成飞行器目标覆盖范围的远边界。
③计算目标覆盖范围的近边界。
在加入式(13)的约束情况下,再取不同的w值分别进行最小射程优化,即优化指标取为J=minL(tf),可得到一组优化弹道,其落点将组成飞行器射程覆盖范围的近边界。
最后,将前面得到的最优弹道的落点连接起来,得到目标覆盖范围。
取飞行器的弹道起点为再入点,参数为:地心距r0=(R0+90)km,经度λ0=0°,纬度φ0=0°,速度v0=4 600m/s,速度倾角θ0=-6.5°,速度偏角σ0=45°,弹道终端高度为0,速度为1 000m/s,倾角为-60°,弹道约束考虑驻点热流、动压、过载以及平衡滑翔约束。计算得到“最偏”的2条弹道的地面轨迹如图2最上端2条曲线所示,最大射程弹道的地面轨迹如图3所示,最小射程弹道的地面轨迹如图2所示。最终得到的目标覆盖范围为椭圆形区域,如图4所示,在部分最小射程弹道中,飞行器实现了大角度的迂回飞行,使得近边界曲线向内凹陷。
图2 组成近边界弹道的地面轨迹
图3 组成远边界弹道的地面轨迹
图4 目标覆盖范围(优化)
由图5可见,热流、动压等再入过程约束对最大射程弹道基本不产生约束作用,最大射程弹道主要受飞行器气动性能和落点约束的影响,落点速度降低或飞行器最大升阻比提升,则可达的射程变大,目标覆盖范围的远边界更远。由图6可知,最小射程弹道主要受热流、动压等过程约束影响,若放松这些约束,则可实现更小的射程,近边界更近。
图5 组成远边界弹道的高度-速度曲线
图6 组成近边界弹道的高度-速度曲线
通过最优弹道的计算能获得较精确的目标覆盖范围,但在整个求解过程中需要完成数十次的优化计算,计算耗时较大,只能用于飞行性能的离线分析,无法应用于飞行任务在线规划等对计算快速性要求较高的场合。
本文提出的基于弹道跟踪的目标覆盖范围求解方法的基本思路是:首先将再入过程约束转化为hv平面内的再入走廊,在再入走廊内产生可行的h-v曲线;然后通过反馈线性法对h-v曲线进行跟踪以获得飞行攻角的取值,取不同的侧倾角,实现侧向机动,其中在零侧倾角情况下对最大升阻比平衡滑翔曲线和再入走廊下边界的跟踪,得到纵向最大射程和最小射程弹道,对应目标覆盖范围的最远端和最近端;最后通过对不同h-v曲线的跟踪和侧向机动得到完整的目标覆盖范围。
再入走廊下边界由动压、热流和过载约束边界组成。再入走廊的上边界通过平衡滑翔约束得到,表达式为
将式中的Cy,max改为其他状态下的升力系数,可得再入走廊内可行的平衡滑翔曲线。再入走廊边界与平衡滑翔曲线如图7所示。
图7 再入走廊边界与平衡滑翔曲线
可将这些曲线拟合成速度的三次函数,用于弹道跟踪:
为实现对h-v曲线的跟踪,本文采用反馈线性化方法来设计跟踪制导律。该方法的基本思想是:首先用代数变换将一个非线性系统的动态特性转换成线性的动态特性;然后用熟知的线性控制理论进行设计。在推导过程中,仅考虑纵平面内运动。动力学方程为
该纵向动力学系统的状态变量为(rvθ)。由于飞行器再入后仅由气动力控制轨迹,选取纵向升力作为动力学系统的输入:
为了实现对h-v曲线的跟踪,选取系统的输出为
式中:h=r-R0。
通过以上输入、输出的选取,h-v曲线的跟踪问题就变为输入为u、输出为h的单输入单输出非线性系统跟踪控制问题。为了产生输出与输入之间的线性微分关系,可重复地对输出函数y进行微分直到结果中出现输入u的显式表达式,然后设计新的输入去抵消非线性。对y=h两边求二阶导数,可得:
式中:θT为需跟踪的弹道速度倾角。
将简化的纵平面动力学方程代入式(19)可得:
可见,经过2次微分后,在式(20)中显式地出现了输出y与输入u的关系。由于返回过程中b(x)不为0,则可以选择新输入u*,使其与输入u具有如下关系:
则式中非线性被消去,输出y与新输入的关系可表示为一个简单的线性双积分器关系:
可将闭环系统的跟踪误差设计为一个二阶振荡环节:
式中:h为当前飞行高度;为参考高度;ζ为阻尼系数,一般取为0.7;ω为自然频率,其值可由飞行器的可用控制能力给出[6]。
根据误差方程可知:而==a(x)+b(x)u,则可得到系统的控制输入:
可以看出,u作为非线性输入,包含r,v,θ3个状态变量,它可通过反馈的方式将动力学系统关于高度h的非线性跟踪控制转化为线性双积分器的控制形式。由于r,v,θ3个状态变量可以在实际飞行过程中测量得到,且标准弹道的二阶导数一般也易求得,因此,控制输入是可以得到的。
进一步可得到跟踪h-v剖面所需的攻角或侧倾角,本文取侧倾角为常值ν0,以攻角为主要控制量,此时跟踪h-v剖面所需的升力系数值为
可通过对升力系数反插值获得攻角。通过以上反馈线性化方法,可实现对h-v剖面的有效跟踪。进行离线弹道计算时,可通过上述跟踪方法得到相应的控制指令,然后积分三自由度运动方程,得到其他弹道参数。
仿真条件与3.3节保持一致。利用基于弹道跟踪的目标覆盖范围求解方法进行计算,计算结果如图8~图11所示。
图8 跟踪弹道的地面轨迹
飞行路程的远近取决于飞行器纵向升阻比的大小和再入弹道约束。从优化计算的结果可以看出,组成目标覆盖范围边界的弹道在h-v曲线上具有一些共同点,远边界弹道的h-v曲线基本沿最大升阻比攻角对应的平衡滑翔曲线波动(图5),而近边界弹道的h-v曲线则基本紧贴着再入走廊下边界(图6)。零侧倾角情况下,通过对最大升阻比平衡滑翔曲线和再入走廊下边界的跟踪,可近似纵向最大射程弹道和最小射程弹道,如图8、图9所示,图中,FL,max为最大升力,FL/FD,max为最大升阻比。踪方法得到的目标覆盖范围比弹道优化方法得到的范围要小,主要原因是由于采用常值侧倾角来进行侧向机动,没有充分发挥飞行器的横向机动能力。
基于弹道跟踪的目标覆盖范围求解方法,在弹道跟踪规律确定的情况下,得到目标覆盖范围只需完成数十次的弹道积分运算,相比弹道优化方法,计算耗时大大减少。
图9 跟踪弹道的高度-速度曲线
图10 侧向机动弹道的地面轨迹
图11 2种方法的目标覆盖范围比较
取常值侧倾角即可实现侧向弹道机动。如果在对最大升阻比平衡滑翔曲线进行跟踪的同时分别取30°和-30°的常值侧倾角,所得弹道的地面轨迹如图10所示。通过对再入走廊内若干h-v曲线的跟踪,采用常值侧倾角产生侧向机动,可得一系列跟踪弹道,将所得弹道的地面轨迹的终点连接起来,可得到近似的目标覆盖范围,如图11所示。采用弹道跟
本文给出的2种计算方法均能完成目标覆盖范围的计算。基于弹道优化的计算方法,在合理设置优化指标和终端弹道约束的情况下,通过求取一系列最优弹道,最终获得目标覆盖范围的边界,其计算精度较高。基于弹道跟踪的计算方法,对再入走廊边界曲线的跟踪,在计算速度方面具有优势,计算简单方便,能够快速地计算出再入落点区域。2种计算方法各具特点,能满足不同场合的计算需求。
[1]BOLLINO K P.High-fidelity real-time trajectory optimization for reusable launch vehicles[D].Monterey,California,USA:Naval Ostgraduate School,2006.
[2]LU Ping,XUE Song-bai.Rapid generation of accurate entry landing footprints[C]//AIAA Guidance,Navigation and Control Conference.Chicago,Illinois:AIAA,2009.
[3]赵汉元.飞行器再入动力学和制导[M].长沙:国防科技大学出版社,1997.ZHAO Han-yuan.Dynamics and guidance of a reentry vehicle[M].Changsha:National University of Defense Technology Press,1997.(in Chinese)
[4]HUNTINGTON G T.Advancement and analysis of a Gauss pseudospectral transcription for optimal control problems[D].Massachusetts,USA:Massachusetts Institute of Technology,2007:51-57,115-143.
[5]RAO A V,BENSON D A,DARBY C L.Algorithm:GPOPS,a Matlab software for solving multiple-phase optimal control problems using the Gauss pseudospectral method[J].ACM Transactions on Mathematical Software,2010,37(2):1-39.
[6]闫晓东,唐硕.基于反馈线性化的 H—V返回轨道跟踪方法[J].宇航学报,2008,29(5):1 546-1 550.YAN Xiao-dong,TANG Shuo.Entry trajectory tracking law for suborbital launch vehicle via feedback linearization[J].Journal of Astronautics,2008,29(5):1 546-1 550(in Chinese)