固体助推飞行器大气层内闭环制导方法

2023-07-24宋满金唐胜景

宇航学报 2023年6期

郭杰，宋满金，柳青，唐胜景

(1. 北京理工大学宇航学院，北京 100081；2. 北京空天技术研究所，北京 100074)

0 引言

随着助推-吸气式高超声速飞行器朝着多任务、多约束的方向发展[1-2],面对复杂终端约束的大气层内上升段精确制导是目前研究的重点。为了保障助推器与载荷分离时的交班点约束,飞行器需要在上升过程中下压弹道飞行,并在弹道终端实现转平,该低弹道转平特性会导致较大的轴向过载。此外,为了满足载荷后续飞行条件,上升段终端攻角和动压需保持在一定约束范围内。而大气层内上升段存在高度较低、气动扰动较大、动压变化剧烈等问题,极易导致轨迹偏移、抖动甚至发散,给制导带来巨大挑战。

针对运载器上升段制导问题,学者们已经进行了相关研究。摄动制导[3]是一类被广泛应用的标称轨迹制导方法,通过将实际弹道在标准弹道附近泰勒展开,控制目标量与理论值的偏差为零,实现约束满足,但该方法依赖标称轨迹,适应能力和鲁棒性较差[4]。为了提高自主性和制导精度,迭代制导、凸优化制导和能量管理等自适应制导方法陆续被提出。吕新广等[5]、韩雪颖等[6]提出了能够约束终端入轨姿态的迭代制导方法;王洪波等[7]基于改进的序列凸优化设计了固体火箭入轨制导方法,有较高的求解效率;Patha等[8]提出一种交变姿态控制的能量管理方法,陈喆等[9]、陈思远等[10]对能量管理方法进行完善,实现了终端姿态约束;张迁等[11-12]在传统速度管控模型的基础上引入定点制导算法,有效提高了入轨精度。上述自适应算法能够根据飞行器当前状态在线计算制导指令,提高运载器的自主飞行能力。然而,上述方法在进行建模时或不考虑气动力影响,或在推导过程中对气动力进行简化,大多适用于高海拔大气稀薄的飞行阶段,在低海拔气动不确定情况下会有较大偏差。

对于大气层内上升段制导方法,除已经实际应用的标称轨迹跟踪制导外,国内外学者主要集中在基于最优控制理论的运载器大气层内轨迹在线规划与闭环制导的研究上,其求解思路分为间接法与直接法。间接法将最优轨迹上升问题转化为两点边值问题进行求解[13-15],但存在推导繁琐、收敛半径小和对初值敏感等不足,且应用在运载器指数攻角转弯后的飞行阶段。直接法将最优轨迹规划问题转化为有限维参数优化问题,离散状态量或控制变量,利用数值计算方法进行轨迹规划,避免了大量的过程推导。崔乃刚等[16]采用两步梯度投影算法(MGPA)非线性规划方法进行在线轨迹规划,并应用基于线性二次型调节器(LQR)的闭环制导方法实现轨迹跟踪;贺前伟等[17]提出一种基于模型预测静态规划算法的自适应制导方法,实现不确定情况下的末段状态高精度制导。随着凸优化理论日渐完善,其也被应用于大气内上升段制导,郝泽明等[18]针对固体火箭上升段实时轨迹优化问题,提出一种基于邻近-牛顿-康托维奇凸规划的轨迹优化方法,用内点法求解二阶锥规划问题,保障轨迹优化实时性。在智能方法方面,何润林[19]用粒子群优化(PSO)算法对射角、全程攻角进行优化,并设计了自适应动态逆制导策略;周宏宇等[20]提出一种全新的飞行剖面,使用改进PSO求解剖面参数,并引入强化学习机制提高求解效率。上述基于数值求解的优化算法虽能够求得最优或次优轨迹,但一般求解参数较多,且存在求解精度与计算效率的矛盾,在线应用困难。

综上所述,对于上升段制导问题已有较多研究,但存在以下不足:部分方法可以对终端姿态进行约束或能够在线应用,但是不适用于稠密大气层内制导;或者能够通过优化的方法求解大气层内上升段最优轨迹或制导指令,但在线应用困难。本文针对固体助推器大气层内上升段制导问题,考虑推力不可调和耗尽关机特性,重点研究了制导指令设计与终端姿态约束问题,提出一种闭环制导方法。设计了解析形式的俯仰角剖面,能够更新剖面约束、精确控制终端攻角,且只需通过数值方法校正两个剖面参数即可更新整个剖面并直接获得控制量,有效提高在线计算效率;采用在线更新策略,根据飞行器当前状态与目标位置实时更新下一制导时刻俯仰角指令,能够在气动参数不确定情况下有较好的鲁棒性。与传统方法相比,该闭环制导方法在大气层内具有较高的制导精度,且具有工程应用潜力。

本文结构如下:首先,设计了一种多项式形式的俯仰角剖面,通过合理地选择剖面终端值保证了终端攻角约束。然后,以数值算法校正俯仰角剖面参数保证高度和弹道倾角约束,同时设计剖面在线更新流程。最后,进行仿真分析验证该方法的制导精度与鲁棒性。

1 上升段数学模型

1.1 运动模型

本文研究固体助推器大气层内上升段制导问题,整个上升过程高度较低、耗时较短,可忽略地球自转,基于“瞬时平衡”假设,只考虑纵向平面运动,得到飞行器在速度坐标系下的质心运动方程为:

(1)

式中:r为位置;G为重力;T为推力;L为升力;D为阻力;m为质量;mc为秒耗量。

运动方程的控制量为攻角α。重力、推力、升力和阻力可以写为飞行时间、位置、速度和攻角α的函数。其中,重力表达式为

(2)

式中:μ为地球引力常数;R为地心到初始发射位置的矢量。

推力表达式为

T=T(t)Γ(α)Iv

(3)

式中:T为推力幅值;Γ(x)为旋转矩阵;Iv为速度的单位矢量,表达式为

(4)

式中:v为速度矢量。Γ(x)表达式如下:

(5)

升力和阻力表达式为

(6)

式中:L,D分别为升力幅值和阻力幅值,表达式如下:

(7)

式中:ρ为大气密度;V为速度幅值;S为飞行器参考面积;Ma为马赫数;CL,CD分别为升力系数和阻力系数,是马赫数和攻角的函数,表达式如下:

(8)

1.2 约束条件

大气层内上升段气动环境复杂,大气密度随高度变化范围大,对流层和同温层以内大气密度变化超过40倍;且助推段马赫数变化剧烈,经历亚声速、跨声速、超声速甚至高超声速等阶段,为了保证飞行器飞行安全与结构完整性,飞行器大气层内上升段需要满足一定的过程约束。此外,飞行器在助推器关机时需要满足交班点条件等任务要求,因此还需满足一定的终端约束。

(1)过程约束

考虑飞行器结构强度和安全等因素,飞行过程中需要考虑过载约束:

(9)

式中:nx,ny分别为弹体坐标系下的轴向过载和法向过载;G为飞行器重量。

(2)终端约束

从飞行器的飞行任务考虑,助推段终端需要满足高度约束和弹道倾角约束。同时为了满足载荷后续工作条件,还需要终端攻角、马赫数-动压保持在一定窗口约束范围内,综上可得终端约束如下:

(10)

式中:h为飞行器高度;θ为弹道倾角;tf为终端时间;q为动压,qfmin和qfmax分别为动压窗口的下界与上界。

(3)约束转化

助推器上升过程中,法向过载大小决定着细长轴体所受弯矩大小,过大会导致轴体断裂,造成严重后果,因此法向过载是重点关注的约束量。由式(7)和式(9)可得:

(11)

定义CLcosα+CDsinα为“过载系数”,式(11)可写为

(12)

上升段为一压低弹道的过程,在大动压段攻角基本为负值,且过载系数是速度和攻角的函数,其变化规律如图1所示,表明在一定范围内,过载系数与速度和攻角呈现单调函数关系。在速度一定时,过载系数绝对值与攻角绝对值正相关,因此可通过限制攻角大小来满足法向过载约束。

图1 过载系数与速度和攻角的关系

在飞行器实际飞行过程中,已知当前高度H和速度V,可通过求解式(13)得到满足约束的当前最小攻角值。

(13)

若当前α不满足幅值约束,则通过式(14)进行限幅处理

(14)

式中:αt为实际攻角值;αmin为当前满足约束的攻角最小值。

2 上升段闭环制导方法

2.1 控制量选择

由固体助推飞行器的动力学方程可知,攻角和俯仰角都可作为飞行控制量,因此需要从数值计算的角度分析攻角和俯仰角对飞行轨迹的影响能力,确定剖面选择[21]。

图2为攻角改变时飞行器体轴以及轨迹改变情况示意图,攻角从α1改变到α2,弹道倾角从θ1改变到θ2。除飞行器体轴改变造成的正增量Δα2外,还有速度方向改变产生的攻角负增量Δα1;图3为俯仰角改变时飞行器体轴以及轨迹改变情况示意图,俯仰角从φ1改变到φ2,弹道倾角同样从θ1改变到θ2,但是过程中只有体轴变化产生的俯仰角正增量Δφ。

图2 攻角改变时轨迹变化

图3 俯仰角改变时轨迹变化

2.2 飞行剖面设计

基于上节分析,对于大气层内上升段制导,考虑大气密度分布特点和飞行器速度规律,本文设计了如下分段多项式形式的俯仰角剖面:

(15)

式中:ki(i=1,2,3,4)分别为3次多项式系数;t1为垂直上升段结束时间;tf为飞行总时间,剖面形状如图4所示。确定式(15)中的3次多项式系数需要4组(t,φ)形式的约束条件,其求解方法如下。

图4 俯仰角剖面

初始上升段中,t1可通过人为选择给定,考虑到飞行器发射高度低,大气密度大,初始飞行阶段处于大气稠密区,气动力和推力之间存在较强耦合,选择合适的t1可以使飞行器在跨声速之后进行转弯,避免受到较大的气动载荷和气动扰动,降低载荷、控制系统设计难度。根据上升段姿态角关系φ=α+θ,可通过终端弹道倾角约束与终端攻角约束得到终端俯仰角大小:

(16)

设θ′f为上一制导周期预测得到的终端弹道倾角值,则θ′f会随着制导进程逐渐接近实际弹道倾角终端值θf,因此可在第一个制导周期内采用式(16)计算剖面终端值,后续采用下式更新剖面终端值,实现Δθf的误差消除:

φ(tf)=θ′f+(αfmin+αfmax)/2

(17)

综上,可得到两组约束条件:(t1,π/2)、(tf,φ(tf))。

另外两组约束条件通过设计剖面参数获得。在上述分析的基础上,增加(t2,φ2)、(t3,φ3)两个控制点。其中t2,t3为人为给定的时间节点,且满足t1≤t2≤t3≤tf;φ2,φ3为2个待设计剖面参数。一旦确定φ2与φ3,则可得到如下方程组:

(18)

求解式(18)即可得到式(15)中的3次多项式系数,进而确定整个俯仰角剖面。该剖面分为两段,其中第一段常值俯仰角用于将飞行器渡过跨声速阶段;第二段多项式形式俯仰角充分利用飞行器性能,进行转弯以及调整飞行状态来满足终端约束。

对于t2,t3的选择,理论上只要满足t2≠t3的条件,就能够求解式(18)得到完整的俯仰角剖面。但若|t2-t3|在数值上过小,求解过程中可能会出现矩阵奇异现象,导致无解;同理,|t1-t2|和|t3-tf|的值也不能过小。因此,本文选择t2,t3的策略如下:

(19)

2.3 剖面参数更新策略

传统基于在线轨迹规划的闭环制导思路为:在飞行过程中根据当前状态和终端约束规划轨迹,将轨迹设计由离线过程变为在线过程,并设计制导律对轨迹进行跟踪[22]。这类方法在实际应用中需要合理选择轨迹规划周期与制导周期,才能兼顾计算效率与制导精度。本文所设计的俯仰角剖面为解析形式,且完全由2个剖面参数(φ2,φ3)确定,并且实际飞行中飞行器可直接根据剖面得到制导指令,即将闭环制导问题转化为剖面参数(φ2,φ3)的搜索问题,大幅提高求解效率。

在每个制导周期内,给定初始剖面参数(φ2,φ3),求解得到完整的俯仰角剖面;以当前状态为初始状态,积分式(1)得到终端状态,从而得到终端高度h(tf)和终端弹道倾角θ(tf)。由于(φ2,φ3)决定了整个俯仰角剖面,也即决定了终端高度与终端弹道倾角,因此可以利用二元非线性方程组求根的牛顿迭代法搜索得到满足终端约束的剖面参数。

在当前时刻与终端时刻确定的情况下,终端高度和终端弹道倾角是两个剖面参数(φ2,φ3)的函数:

(20)

本文在每个制导周期内,以高度、弹道倾角双约束,基于俯仰角剖面同时对两个剖面参数进行校正:

(21)

式中:hf为期望终端高度,θf为期望终端弹道倾角,二者可由终端约束得到。上述两个变量在每个制导周期内均为已知量,则方程组式(21)可化为关于俯仰角参数(φ2,φ3)的二元非线性方程组:

(22)

式中:F(φ2,φ3)=h(tf)-hf,G(φ2,φ3)=θ(tf)-θf。利用牛顿迭代法可快速求解上式,且能得到收敛解[23-24]：

(23)

式中:F,G分别为式(22)中F(φ2,φ3),G(φ2,φ3)函数;φ2,i和φ3,i为第i次迭代时的φ2,φ3值;F′φ2,F′φ3,G′φ2,G′φ3为F,G函数分别对φ2,φ3的偏导数。

考虑到非线性函数的复杂性,无法直接求得准确的偏导数表达式,本文采用有限差分法近似得到偏导数值,且精度满足如下要求。

(24)

式中:Δ为一大于0的小量。实际计算过程中,通过式(23)得到下一次迭代的φ2,φ3值,重复迭代直至F<ε1且G<ε2,即可得到满足精度要求的φ2,φ3解。ε1,ε2为可接受的终端高度与终端弹道倾角误差。若在飞行过程中因参数扰动过大,或因控制量限幅导致后续制导周期得到不满足精度的解,此时取目标函数:

(25)

利用参数优化的牛顿迭代法可快速求解使目标函数式(25)最小的剖面参数[25]：

(26)

(27)

式中:a1,a2,a3,b1和b2为滤波器待定参数。式(27)的时域表达式如下:

(28)

式中:x(t)为输入信号;y(t)为输出信号。将上式转换成对应的离散域差分方程,并将时域离散信号转变为序列信号可得:

y(n)=λ1x(n)+λ2x(n-1)+λ3x(n-2)+

λ4y(n-1)+λ5y(n-2)

(29)

式中:λi(i=1,2,3,4,5)为相应的系数,表达式如下:

(30)

式中:Ts为采样周期。使用式(29)对所求得的俯仰角进行滤波,最终可以得到平滑的制导指令。整个基于俯仰角剖面更新的闭环制导算法总体流程如图5所示。

图5 闭环制导方法流程

在制导过程中,为了提高预测环节计算效率,减少对箭上计算能力的占用,前期可采用定步数积分,积分步长Δh取(tf-t)/N,t为当前时刻,N为积分总步数,可取200;当后期Δh减小到与制导周期一致时,保持不变。另外,为减少剖面更新迭代次数,可取上一制导周期剖面参数作为迭代初值,提高搜索效率;还可以在上升段前期多个制导周期更新一次俯仰角剖面。随着飞行时间增加,剖面更新耗时减少,可减少为每个制导周期内更新一次剖面。本文采取剖面更新频率如下:

(31)

上式表示κ个制导周期更新一次剖面参数。

上述基于俯仰角剖面参数更新的闭环制导核心思路是将牛顿迭代与剖面设计相结合,将每个制导周期内制导指令求解转化为终端高度、终端弹道倾角双约束的双剖面参数搜索问题。这种思路还可以适用于离线轨迹规划,既可基于攻角-时间剖面,也可基于俯仰角-高度剖面,甚至基于俯仰角-速度剖面进行设计。剖面参数也可以选择任意两个同时影响整个剖面形状的参数。

3 仿真校验

以固体火箭助推飞行器为研究对象,地球半径取6 371 km;地球引力常数μ取3.986 005×1014m3/s2;g为重力加速度;声速模型采用美国标准大气1976,在50 km以下具有较高的拟合精度。

飞行器总飞行时间为90 s;垂直上升段时间t1取20 s;初始质量m0为4 150 kg;秒耗量mc为30 kg/s;推力幅值T为95 kN;参考面积S为0.8 m2;制导周期取0.1 s;初始剖面参数的迭代初值由线性插值得到,φ2取60°,φ3取30°;滤波器参数a1,a2,a3,b1和b2分别取0.1, 0.2, 1, 1, 1.4,采样周期Ts取0.1 s。

以飞行器出发射筒时状态作为初始状态,如表1所示;发动机关机时飞行器需处于转平状态,终端状态约束如表2所示。过载约束为:轴向过载最大值nxmax为3;法向过载最大值nymax为2。

表1 初始状态

表2 终端状态

设置3组仿真条件:第1组为标准条件;第2组为大气密度、升力系数和阻力系数各拉偏-10%,推力拉偏-5%,质量秒耗量拉偏-2%;第3组为大气密度、升力系数和阻力系数各拉偏+10%,推力拉偏+5%,质量秒耗量拉偏+2%。

本文数值计算在Intel i5 3.10 GHz台式机电脑下进行,程序均在MATLAB 2020b环境下编译运行。

3.1 滤波效果分析

为验证本文所提出的闭环制导方法中的滤波效果,分别进行有滤波器和无滤波器的对比仿真。第3组条件下仿真曲线如图6所示,所有组仿真结果见表3。

表3 滤波前后仿真结果

图6 滤波前后对比仿真曲线

由表3可知,滤波器使用前后,飞行器终端状态接近且均具有较高精度,而由图6中(a)、(b)可知,滤波器使用前,俯仰角和攻角出现震荡现象,这是因为在干扰条件下,为保证终端精度,相邻两次剖面更新所得的俯仰角剖面存在差异,导致时序制导指令产生波动。加入滤波器后,俯仰角和攻角较为平滑,符合实际飞行需求。

3.2 对比仿真分析

为验证本文所提方法的有效性,选取基于导引系数快速求取的摄动制导[3]、基于PD的最优轨迹跟踪[26]和本文闭环制导3种制导方法进行仿真对比。其中最优轨迹以终端速度最大作为目标函数,通过Gauss伪谱+序列二次规划(SQP)方法在标准条件下离线规划。第3组条件下仿真曲线如图7所示,所有组仿真结果见表4。

表4 不同仿真条件下的仿真结果

图7 第3组拉偏情况下仿真曲线

由表4中3组仿真结果的高度误差可知,3种制导方法在标准状况下均有较高精度,但在干扰较大的情况下,基于导引系数快速求取的摄动制导高度误差接近80 m,说明简化气动力的该方法不适用于大气层内上升段制导。基于PD的最优轨迹跟踪和闭环制导方法均能较好地约束终端高度和终端弹道倾角,其中闭环制导方法的终端高度误差在10 m以内,终端弹道倾角误差不超过0.1°,相比于PD跟踪均提高了一个数量级。这是因为本文提出的闭环制导方法,能够在每一个制导周期内,通过搜索双剖面参数实现对终端高度和弹道倾角的精确控制,鲁棒性较强。由马赫数和动压值仿真结果可知,飞行全程动压值较小,3种制导方法均满足终端窗口约束;3组终端攻角也均满足窗口约束,且在标准条件下攻角波动均不超过0.1°。在有干扰的情况下,闭环制导的终端攻角波动比PD跟踪更小,说明本文提出的更新剖面终端值方法能更好地约束终端攻角。

由图7(a)(d)可知,标称最优轨迹与3种制导方法所得轨迹平滑、变化平稳;图7(b)(e)(f)表明所有轨迹均满足过程约束;图7(c)表明4种方法得到的攻角曲线连续且没有抖振出现,其中基于PD的最优轨迹跟踪攻角曲线变化趋势与标称最优攻角接近,且与闭环制导攻角曲线差异较大。这是因为PD跟踪依赖标称轨迹,在干扰存在时调整能力有限,而闭环制导自主性较强,能够根据状态实时调整控制量,充分发挥飞行器飞行能力;俯仰角剖面平均更新时间为28 ms,证明该方法具有在线应用价值。

3.3 制导精度与鲁棒性分析

针对气动系数偏差、大气密度偏差、推力偏差和质量秒耗量偏差进行200次蒙特卡洛仿真。假设各项参数扰动符合正态分布,扰动参数如表5所示。仿真结果如图8所示,终端约束平均值见表6。

表5 扰动参数

表6 终端约束蒙特卡洛仿真平均值

由图8(a)(d)可知,高度和弹道倾角曲线平滑、轨迹无震荡;图8(b)(e)(f)表明在气动不确定情况下,所有仿真均满足过程约束,且动压、过载变化平缓,终端动压均在窗口范围内;图8(c)表明所有仿真攻角曲线无抖振出现,终端攻角均满足窗口约束;图8(g)(h)表明200组仿真终端高度误差均不超过10 m、终端弹道倾角误差均不超过0.15°,具有较高精度;图8(i)为终端攻角分布,各种随机扰动下,所有终端攻角的最大波动值不超过0.3°,证明本文所提方法能够较好地约束终端攻角,且具有较强的鲁棒性。