考虑阻力系数时变的下压段半解析时间预测方法

崔正达, 魏明英, 李运迁

(1. 北京电子工程总体研究所, 北京 100854; 2. 北京仿真中心, 北京 100083)

0 引 言

高超声速滑翔飞行器通过升力增程的方式在高度20～100 km的临近空间飞行,该空域空气稀薄,飞行器受到的阻力较小,但仍具备机动能力,且受制于地球曲率,雷达难以远距离探测和预警[1]。诸多优势使得高超声速飞行器一经问世就是学者们研究的热点,被誉为航空工业皇冠上的明珠。高超声速滑翔飞行器运动方程形式复杂、非线性特性强、再入过程对控制量高度敏感、过程约束多,受过载、热流密度等约束影响较大、对终端状态要求高等限制因素使得其轨迹设计成为该领域的难点和热点之一。

为了在诸多约束条件下规划航迹,精准达成作战目标,目前主要有3种技术途径:标称轨迹法、预测校正法和前两者的混合制导方法[2]。标称轨迹法是将上述问题转化为多约束非线性规划问题,经过大量的离线运算或者在线优化进行求解,得到参考飞行轨迹。标称轨迹法在20世纪70年代首先被发展并应用在航天飞机的再入返回过程之中[3-4],1977年Page和Rogers针对大气层外飞行器再入问题提出了基于三次曲线的弹道参考制导律[5],该制导率可以满足终端角度约束,类似的还有基于圆弧制导轨迹的方法[6]。但大量的仿真实验表明该类制导律精度欠佳[7]。1979年Harpold和Graves提出了经典的基于速度多项式的阻力剖面再入制导方法[8],研究者在此基础上基于不同的关注方向,发展出了表征飞行器运动特性的很多种广义飞行剖面,在飞行任务中只需跟踪迭代生成的参考轨迹即可达到预定效果。较常见的有升阻比-能量剖面[9-10]、阻力-能量剖面[11-12]、阻力-倾侧角指令剖面[13]、阻力-速度剖面[8]、攻角-速度剖面[14-15]、动压-高度剖面[16]、高度-射程剖面[17]等。

预测校正法的核心思想是利用弹载计算机的运算能力在线预测飞行轨迹以及终端状态,并求解出其与期望值的偏差,然后校正制导指令不断消除预测飞行轨迹与理想值的偏差。该方法相较于标称轨迹法具有更高的鲁棒性,而根据轨迹预测的方式不同,预测校正法又可以分为数值解[18-22]和解析解法[23-26]。解析预测制导需要常系数假设[24]或者将轨迹调制到特定形式获得近似解析解预测[17];数值预测的基本思路是先由预先指定的剖面求得控制指令序列,再代入弹道微分方程求解数值积分对终端状态进行预测,该方法在线计算量很大,对于弹载计算机的计算能力有较大的依赖。

上述两类典型制导方法各存在优点和不足,有学者尝试将二者互相取长补短形成混合制导方法兼顾标称轨迹法的快速性以及预测校正法的鲁棒性优点:如Hu等[27]就将滑翔段再分成平衡段和线性段,分段采用不同的方法以获得二者的优势。王青等[28]、Wang等[29]也是采用分段的思想将两种制导方法混合,兼顾速度及精度。

多飞行器协同制导概念出现以后,时间控制需求迫切,针对滑翔段的时间预测算法也相继出现[30-32]。但是滑翔段的优化算法由于使用平衡滑翔或定阻力系数假设,不再适用于复杂变化的下压段。

快速下压弹道的主要任务是在复杂约束下将速度、高度和弹道倾角等飞行器运动参数控制到指定值,其间要经历非常强的气动非线性影响。国内外在下压段的研究以弹道整形制导律[33-36]和带落角约束的制导律[37-39]为主。这些经典理论提出的时间都比较早,时间不是主要的受控目标。上述制导律可以满足落角和落点约束,但无法对遭遇时间进行准确的估计和控制,不能满足时间协同的新兴需求。

基于上述情况,本文提出一种考虑阻力系数时变的下压段时间预测方法,设计随攻角和马赫数变化的阻力系数模型,并简化了飞行器运动方程以精确快速预报导弹运动参数和时间。本文按如下顺序编排:首先建立飞行器下压段数学模型;然后基于解析理论给出以剩余射程为自变量、满足终端约束的高度-射程剖面表达式,将空气密度随高度变化规律线化为剩余射程的多项式函数,并且结合气动拟合结果,在阻力系数表达式中引入攻角和马赫数的影响项,基于上述解析理论推导降阶弹道微分方程,得到了以剩余飞行距离为自变量的一阶微分方程,该方程可以通过数值积分快速求解;最后以典型弹道为例给出对比传统方法的数值仿真结果以及对结果的讨论。

1 滑翔飞行器数学模型

高超声速飞行器在航迹末段快速俯冲,其飞行时间较短、飞行高度和射程较低、纵侧向耦合不深,因此可独立设计两个维度的轨迹并假设地球为不旋转的匀质圆球,得到拦截器下压段纵平面内动力学方程:

(1)

式中:V为飞行器相对于地面的速度标量;r为飞行器质心到地心的距离;θ为当地弹道倾角,定义为速度方向与当地水平面的夹角,向上为正;m为飞行器质量;g为当地引力加速度;Y和D为飞行器所受的升力及阻力;CD为飞行器阻力系数;CL为飞行器升力系数,该系数可认为是攻角和马赫数的表达式;Sref为飞行器参考面积,由于只考虑纵平面内的运动所以控制量选为攻角α。

考虑气动力时,采用如下的指数大气密度模型:

ρ=ρ0e-β h

其中，β=1.406 4×10-4m-1;ρ0=1.225 kg·m-2。

定义从当前位置到终端位置间剩余距离在地面上的投影长度为射程,记为RL。

2 高度射程剖面设计

快速下压段需要在约束范围内将高度和弹道倾角规划到指定区域,所以针对下降段距离近、弹道形状复杂、约束条件严苛的应用背景下拓展了文献[17]提出的高度-射程剖面表达式,在弹道震荡平缓预设下定义飞行轨迹为剩余射程的六次函数:

(2)

其中,

a1RL+a0

可将下压段非线性特性最强的空气密度简化为航程的多项式形式便于解析降阶,有:

ρ=ρ0e-β h=ρ0F(RL)

(3)

式(2)所确定的弹道形状需要确定a0～a6共7个弹道参数,选取初始高度h0、初始飞行路径角θ0、1/3处高度h1、2/3处高度h2、终端高度hf、终端路径角θf以及弹道h1点处路径角θH。线性解算矩阵同文献[17]所述,选取次高点倾角θH的原因是该值可以有效控制过载分布,充分利用飞行器机动能力。若不如此拓展,根据高度射程表达式的数学特性,会出现高空可用过载能力小但需用过载能力大的情况,不适合实际工程应用。

下压段中当地弹道倾角θ、倾侧角σ均为小量,则飞行路径角的解析解为

(4)

(5)

由飞行器运动方程可知弹道倾角及其变化率可以与飞行器飞行高度变化率及升力联系起来:

(6)

(7)

联立式(1)、式(2)、式(5)和式(7)可得到升力Y的简化表达式:

(8)

式中:

3 考虑攻角及马赫数变化的阻力系数模型

现有弹道参数解析预测方法在推导弹道解析解的时候,为简化推导过程习惯假定阻力系数为常值,这样的假设对于高度和速度变化相对稳定的滑翔段是合适的。但是,当阻力系数的两个主要影响因素即马赫数和攻角发生改变时,该假设与真实的偏差会对结果产生较大影响,影响预报精度。设气动力系数CL和CD可表示为如下形式的攻角、马赫数函数:

(9)

用攻角跟踪弹道所需升力Y时有如下关系:

(10)

当攻角α=0时不产生升力,即Y≈0,由式(9)得CL0≈0。

然后将攻角表达式(10)代入式(9),考虑攻角α对阻力系数CD的影响得到:

(11)

联立升力值表达式(8):

(12)

分析上述表达式,式(12)中第1部分即常值部分,是传统预测方法所考虑的部分;第2部分是马赫数变化带来的阻力系数的变化量;第3部分是改变弹道形状和平衡重力所需的升力,进而诱导产生的额外阻力。

根据下压弹道常见高度范围,得到以剩余射程为自变量的阻力系数拓展表达式:

(13)

根据动力学方程可得

(14)

有:

(15)

由式(15)联立式(13)可得

(16)

(17)

整理得

(18)

联立下压段高度表达式(3):

(19)

式(19)就是解析理论简化得到的一阶微分方程,数值积分可以得到准确的速度随射程变化规律。方程中第1项兼顾了传统常值假设下所考虑的情形,在其基础上第2项是由马赫数变化主导的影响项,第3项则是由攻角主导的影响项。第4项为重力势能和动能的转换项。为了简化表述定义如下常数:

式(19)可以写为

(20)

观察以RL为自变量,V为因变量的微分方程,该速度-剩余射程的方程在弹道解析理论下变为一阶微分方程,可以用数值解法来预报其运动状态,从而以剩余飞行距离RL为自变量对飞行器速度和时间进行估计。若省略后3项,则该方程退化为原始文献形式。

4 仿真分析

4.1 3种仿真算例

根据文献[30],设置仿真初始条件为:经纬度λ0=φ0=0°,剩余距离RL为210 km,初始高度为h0=35 km,速度V0=2 000 m/s,速度倾角θ0=0°。仿真计算机所用的CPU型号为Intel(R) Core i5-4590,主频3.29 GHz。软件为Matlab 2012b。

在本节中,改变弹道参数h1、h2、θH可得到不同的弹道形状,模拟3种不同的弹道形式:第1种是滑翔段飞行弹道,第2种是高度大范围变化的末段俯冲下压弹道,第3种是先上扬再快速下压的末段俯冲突防弹道。对这3种典型弹道分别比较文献[17]所应用的常值阻力假设、本文提出的简化微分方程和弹道数值积分得到的预报结果,得到的仿真曲线如图1～图6所示。

图1 近似平衡滑翔弹道Fig.1 Trajectory in quasi-equilibrium glide phase

图2 近似平衡滑翔速度变化曲线Fig.2 Velocity change curve in quasi-equilibrium glide phase

图3 俯冲下压弹道Fig.3 Trajectory in dive phase

图4 俯冲下压速度变化曲线Fig.4 Velocity change curve in dive phase

图5 俯冲突防弹道Fig.5 Trajectory in penetration phase

图6 俯冲突防速度变化曲线Fig.6 Velocity in penetration phase

在滑翔段速度变化线性下降,由常值阻力系数假设得到的速度变化与本文提出的预报方法以及弹道仿真曲线十分接近,因为在滑翔段本文提出的方法拓展量皆为小量,方法退化为采用常值阻力假设的解析方法;在俯冲下压段由于产生了负攻角以及较大的弹道倾角,在重力和空气动力的共同作用下速度的变化规律呈现出了很强的非线性,本文提出的方法在设定阻力系数模型的时候考虑的因素较多,能够比常值阻力假设更好地揭示出真实变化规律,故本文提出的拓展方法所预报的速度与仿真结果更为相近,预报误差较小。另外,仅在纵平面内借助重力势能和动能之间的转换规划弹道形状可以减小协同过程中的能量消耗。弹道仿真也说明弹道形状、起始-终止的高度落差对导弹的抵达时间以及终端速度产生较大影响。

4.2 比较提出方法的计算时间效率

准确预报速度变化可以为抵达时间准确估计打下坚实基础,对速度-剩余射程进行一阶积分,以3种飞行器典型弹道为例对抵达时间进行预测,并记录系统计算耗时,得到飞行时间估计结果和计算耗时如表1～表3所示。

表1 近似平衡滑翔段计算效率比较

表2 俯冲下压段计算效率比较

表3 俯冲突防段计算效率比较Table 3 Computational efficiency comparison in penetration phase s

对比上述结果可知,文献[17]预报的结果在下压段和俯冲段与弹道仿真的区别比较大,常值假设无法适应,会出现5～10 s的预报误差(相对误差4.5%～9%),而简化解析方程仅有不到1 s的预报误差(相对误差1%～0.9%)。通过与另外两种预报方法进行综合比对,本文所提方法可以在不显著提升计算时间复杂度的同时显著提升预报精度。

4.3 蒙特卡罗仿真结果

为验证方法有效性,对各情况中气动参数加入10%的拉偏量,各仿真1 000次得到的预测偏差如图7～图12所示。

图7 平衡滑翔-常值阻力系数预测误差分布Fig.7 Estimated error distribution in quasi-equilibrium glide phase with constant drag coefficient

图8 平衡滑翔-提出方法预测误差分布Fig.8 Estimated error distribution in quasi-equilibrium glide phase with proposed method

图9 俯冲下压-常值阻力系数预测误差分布Fig.9 Estimated error distribution in dive phase with constant drag coefficient

图10 俯冲下压-提出方法预测误差分布Fig.10 Estimated error distribution in dive phase with proposed method

图11 俯冲突防-常值阻力系数预测误差分布Fig.11 Estimated error distribution in penetration phase with constant drag coefficient

图12 俯冲突防-提出方法预测误差分布Fig.12 Estimated error distribution in penetration phase with proposed method

仿真结果表明,采用本文提出的方法可以有效地在各种弹道形状下估计遭遇时间,并且估计散布小于常值阻力假设。常值假设标准差为0.421 92,本文提出的方法估计标准差为0.222 49。

图13 动压预测曲线对比Fig.13 Comparison of dynamic pressure estimation curve

本文采取的剖面形式可准确解析飞行器所在位置的大气密度,因此对速度的精准预报提升了动压的解析预测精度。动压的精准预报可以从可用过载以及最大动压两个方面为优化算法判断约束条件提供轻便准确的模型,加快优化速度及精度,有利于在弹道规划过程中更大程度地发挥飞行器过载能力。

5 结 论

本文深入分析高超声速滑翔飞行器下压段速度变化规律,对精度影响最大的两个因素,即攻角、马赫数,分别在阻力系数表达式中进行拓展,给出了新的微分方程形式。基于解析解理论对复杂的弹道微分方程进行简化,得到仅有一阶的简化微分方程,该方程可以通过数值积分快速求解。通过典型弹道仿真表明,该方法可将俯冲下压段的时间预报精度从10 s提高到1 s左右,同时不显著提升计算复杂度,实现滑翔飞行器俯冲段时间的精准、快速预报,提高全弹道协同规划能力。