一种适应推力可调模式的飞行器精确软着陆闭环制导方法与流程

本发明涉及运载火箭具有推力可调能力且采用垂直回收模式的末制导技术，该发明同时也涉及火星及月球软着陆末制导技术。

背景技术：

火箭助推级垂直定点回收并重复使用是一种降低发射成本有效的途径。由于火箭动力着陆段飞行时间短，要求落点精度高，因此火箭垂直定点着陆过程对制导控制方法提出了很高的要求。由于火箭垂直返回定点着陆的制导任务十分复杂，寻找能够满足过程约束以及终端约束并且在收敛特性上具有快速型以及鲁棒性的制导算法，一直是众多学者探索并且研究的方向。近年来，一些火箭垂直回收与定点着陆制导算法相继被提出，在多约束制导算法方面取得了很大的进步，然而依然存在有一些问题亟待解决：

目前只有美国spacex公司通过猎鹰-9火箭对其一级助推器真正实现了进行的多次垂直回收，其末制导阶段采用了一种基于凸优化的轨迹规划算法实现在线闭环制导。既根据火箭动力学模型将直接法轨迹优化(非凸)问题转化为由凸函数描述的凸优化问题，从而使原始的非线性规划问题的求解，具有一般数值优化算法所具有的收敛域大的优点，并且还具有更快的收敛速度。在实际计算过程中，凸优化算法具有多项式时间复杂度，从而保证在给定精度条件下具有求解时间可控的收敛效率，能够在预先确定的迭代次数上限内收敛且迭代次数几乎与问题的维数无关。另外，求解凸优化问题时采用的原始-对偶内点法，对于一些特定的凸优化问题该求解器不依赖于任何初始猜想条件并且能够可靠求解。

技术实现要素：

要解决的技术问题

本发明主要解决目前采用凸优化技术存在的不足，既当回收火箭的动力学特性具有高度非线性且强干扰时，考虑过程约束是凸化难度大。此外该方法求解采用数值内点法进行求解，计算量大，当不确定干扰时会导致收敛时间变慢，甚至无法进行快速求解。

为了克服以上凸优化方法确定运载火箭垂直回收段制导律所存在的实时性、收敛性及求解方式复杂等问题，本发明结合实际运载火箭可调推力的特性，从最优控制的角度提出了一种实时的闭环解析制导方案。

技术方案

一种适应推力可调模式的飞行器精确软着陆闭环制导方法，其特征在于步骤如下：

步骤1：计算视线坐标内的三个速度分量vsx,vsy,vsz：根据探测系统测量得到在地面雷达坐标系内箭体的三个速度分量vtx,vty,vtz以及三个位置分量xt,yt,zt，并通过等式计算：

式中弹目连线的高低角和方位角

步骤2：根据地面雷达坐标系内箭体的三个速度分量vtx,vty,vtz计算弹道坐标系内的速度倾角θ以及抗奇异倾角θ′：

步骤3：根据地面雷达坐标系内箭体的三个速度分量vtx,vty,vtz计算弹道坐标系内的速度偏角σ以及抗奇异偏角σ'：

步骤4：根据地面雷达坐标系内箭体的三个速度分量vtx,vty,vtz及高低角λd计算视线坐标系内高低角速率和方位角速率

式中

步骤5：计算剩余时间的近似值tgo＝rmt/(-vsx)：

步骤6：根据高低角λd和方位角λt、高低角速率和方位角速率剩余时间tgo计算视线角速率和

式中垂直着陆终端视线角约束为γdf＝0°；γtf＝90°；

步骤7：根据视线角速率和方位角λt计算弹道倾角速率和弹道偏角速率

步骤8：根据速度倾角θ以及抗奇异倾角θ′、速度偏角σ以及抗奇异偏角σ'计算速度滚转角γx：

步骤9：计算弹道坐标系的纵向过载指令nyc和侧向过载指令nzc：根据惯性导航器件测量得到飞行器当前速度大小v以及步骤2和步骤3计算得到的θ,θ'和σ,σ'，通过等式计算：

式中g为重力加速度；

步骤10：计算飞行器的理论推力fthe与剩余着陆时间s：根据惯性导航器件测量得到飞行器当前高度h以及飞行速度v，并求解下式二维非线性方程组得到理论推力与着陆时间：

式中为飞行器的额定秒流量；

步骤11：计算节流阀调节系数κ：根据飞行器动力系统的比冲vex和额定秒流量以及步骤10的理论推力值fthe，并通过等式计算：

步骤12：输出步骤9计算的纵向过载指令nyc和侧向过载指令nzc，输出步骤10的理论推力fthe，并按照步骤11计算的κ调节发动机节流阀。

有益效果

本发明制导方案从总体步骤来看只需要根据弹目连线计算纵侧向过载；再根据实时高度、速度计算所需推力，计算过程是简单并且易于理解的。从对制导方案的实施条件来看，其具有闭环解析的形式，这使得改进过程得以层层推进并具有相应物理意义，且不存在凸优化当中所存在的制导律收敛性问题。其次，由于推力指令的作用方向始终位于速度的反方向上，这使得基于当前攻角、侧滑角，能够进行火箭姿态的实时解算，换句话说就是不存在制导控制一体化的问题。另外，其计算量在制导周期以内较优化方式设计的制导律要小得多，具有实时性好的优点。最后，通过对可用中间推力的轨迹跟踪，能够使得制导方案能够满足推力边界要求，同时具有较强的鲁棒性。

附图说明

图1控落角导引律示意图

图2控落角导引律分平面示意图

图3制导计算流程示意图

图4组合偏差状态仿真动力曲线簇

图5组合偏差状态仿真气动角曲线簇

图6组合偏差状态仿真纵向角度曲线簇

图7组合偏差状态仿真速度高度曲线簇

图8组合偏差状态仿真动压与纵向过载曲线簇

具体实施方式

现结合实施例、附图对本发明作进一步描述：

为满足定点与末端轨迹垂直水平面的要求，制导方案在运载火箭的纵侧向过载指令计算上采用变换了制导坐标系的控落角制导律；为满足运载火箭末端速度、高度为零的终端要求，基于原一维软着陆问题的“bang-bang”问题，实施考虑气动力及推力边界的实时推力计算方法，使速度方向的制导指令的提供能在满足末端指标要求的同时令总体制导方案具有较好的鲁棒性。

1、控落角制导律在应用中的改进形式

从图1和图2可得弹道倾角及偏航角速率由以下角度关系给出，

2、更换控落角制导指令生成坐标系以避免奇异问题

前一节当中弹目连线是通过先偏航λt，再俯仰λd的方式得到的。现为了避免弹目连线垂直于水平面时奇异问题的发生(λt此时无定义)，采用先俯仰λt，再偏航λd的方式。相应的，弹道坐标系也需要进行不同选择顺序的改变。在这样的旋转顺序下则原制导指令的计算方式变为以下形式，

这里的下标2代表更换旋转顺序后与新弹道坐标系定义的相关量，其中新系下的弹道倾角和方向角如下，

新弹道倾角与方向角角速率的计算如下，

新弹道坐标系到旧弹道坐标系的转换由一个滚转角γx进行，其计算方式为，

则旧弹道坐标系过载指令为，

更换旋转顺序后，为使运载垂直落地，可令，

3、任意高度速度下的推力计算问题

原最短下落时间制导律描述的是一维的直线下落问题，其解为在开机高度之前零推力，而在开机高度时及其以后按最大推力给出，低于开机高度则没有定义。这样的制导方式虽然能够满足燃料最省指标，但是若考虑气动力及其他系统偏差造成影响时，其鲁棒性表现较差，且做不到在任意高度速度下给出一个合理推力的效果。

本发明对以上方程组进行数值求解，即固定当前飞行高度、速度，采用标称情况下的加速度、比冲等系统参数，则得到一个关于剩余时间s以及推力a的二元非线性方程组。通过阻尼牛顿法或n-sor等收敛域大且收敛快速的迭代方法得出符合任意高度、速度下相应的理论最快下落时间制导律所对应的推力大小fthe。需要注意的是，这里的飞行高度在应用中使用的是当前箭体到目标着陆点的相对距离。

4、气动力补偿问题

针对原高度方向制导律模型对气动力影响的忽略，本发明采用气动力系数与当前状态算出当前气动阻力(速度方向)的估计大小，并将其与前一节得到的理论计算推力相减，实现推力对气动力的补偿。即：

5、两种制导律应用上的互补整合

5.1纵侧向过载补偿

纵侧向过载的补偿采用线性补偿的方式如下，

其中，nyz、nyzc分别代表实际与指令纵侧向总过载，n1代表f1作用下最大攻角带来的过载大小。

5.2考虑推力上下界的轨迹修正

通过固定非线性方程组的推力为可用推力中值，剩余时间为根据当前高度速度进行推力计算时得到的估计剩余时间，结合标称条件下的重力加速度、比冲等参数可以计算出推力取值更接近于推力中值的高度与速度大小，通过构建与其相关的误差项达到轨迹修正的效果，并保证所需推力远离边界，位于中值附近。比例项系数的取值大小可以理解为每分母的误差大小给予分子大小的推力。

5.3节流阀精度的补偿

针对节流阀精度的影响下对标称速度跟踪效果不好而导致末端指标难以满足要求的结果。通过增加标称速度误差的积分项，可以达到良好的修偏效果。

为了防止积分项出现饱和的问题，可以在速度跟踪误差绝对值小于一定阈值时才加入积分项，否则积分项清零。

具体实施步骤：

步骤1：计算视线坐标内的三个速度分量vsx,vsy,vsz：根据探测系统测量得到在地面雷达坐标系内箭体的三个速度分量vtx,vty,vtz以及三个位置分量xt,yt,zt，并通过等式计算：

式中弹目连线的高低角和方位角

步骤2：计算弹道坐标系内的速度倾角θ以及抗奇异倾角θ′：根据探测系统测量得到在地面雷达坐标系内箭体的三个速度分量vtx,vty,vtz，并通过等式计算：

步骤3：计算弹道坐标系内的速度偏角σ以及抗奇异偏角σ'：根据探测系统测量得到在地面雷达坐标系内箭体的三个速度分量vtx,vty,vtz，并通过等式计算：

步骤4：计算视线坐标系内高低角速率和方位角速率根据步骤1计算得到的视线坐标内的三个速度分量vsx,vsy,vsz以及高低角λd，并通过等式计算：

式中

步骤5：计算剩余时间的近似值tgo＝rmt/(-vsx)。

步骤6：计算视线角速率和根据步骤1的高低角λd和方位角λt、步骤2计算得到的高低角速率和方位角速率以及步骤2计算得到的剩余时间tgo，通过下式：

式中垂直着陆终端视线角约束为γdf＝0°；γtf＝90°

步骤7：根据视线角速率和方位角λt计算弹道倾角速率和弹道偏角速率

步骤8：计算速度滚转角γx：根据步骤2和步骤3计算得到的θ,θ'和σ,σ'，并通过等式：

步骤9：计算弹道坐标系的纵向过载指令nyc和侧向过载指令nzc：根据惯性导航器件测量得到飞行器当前速度大小v以及步骤2和步骤3计算得到的θ,θ'和σ,σ'，通过等式计算：

式中g为重力加速度。

步骤10：计算飞行器的理论推力fthe与剩余着陆时间s：根据惯性导航器件测量得到飞行器当前高度h以及飞行速度v，并求解下式二维非线性方程组得到理论推力与着陆时间：

式中为飞行器的额定秒流量。

步骤11：计算节流阀调节系数κ：根据飞行器动力系统的比冲vex和额定秒流量以及步骤10的理论推力值fthe，并通过等式计算：

步骤12：输出步骤9计算的纵向过载指令nyc和侧向过载指令nzc，输出步骤10的理论推力fthe，并按照步骤11计算的κ调节发动机节流阀。

实施实例为：

1)初始条件：回收初始位置(x,y,z)＝(161.521,2204.74,-731.565),初始弹道倾角-102.4407°,初始航向偏角-62.6448°,速度251.6982(m/s),目标位置(x,y,z)＝(0,0,0)，推力值取198kn，节流阀系数可变化范围为[0.7,1]。

2)偏差条件：动力偏差的评估采用初始质量正负5‰、秒流量正负3％、比冲正负19.6m/s和节流阀精度正负2％进行拉偏，另外取推力时滞1s。

采用本发明的制导方法，测试得出的结果图4-8。

图4至图8展示的是一个运载火箭的一级发动机在不同的拉偏条件下进行的近地(高度2km左右)软着陆实例效果。图4可以看到，指令节流阀系数在前半段飞行过程中，纵侧向过载不够时将节流阀系数拉满，而在后半段节流阀系数的变化较平缓，末端接近推力中值所在位置，并且能够针对不同拉偏情况对推力进行相应的调节。这说明根据所设计制导律计算得出的实时推力变化能够满足适应性以及推力边界的要求。从图5可以看出，飞行过程中的攻角与侧滑角变化连续性较好，并且末端值趋近于零，这说明末端的纵侧向过载需求小，反映了据制导律得出的纵侧向过载能满足定点着陆的需求。图6的俯仰角及弹道倾角在末端逐渐收束至-90度，这说明制导律的控落角效果能够实现飞行器的垂直水平面的着陆轨迹。图7中表明飞行器在发动机开机后速度和高度都能在末端趋近于零，这说明制导律能够满足软着陆的需求。图8所展示的动压及纵向过载值说明飞行过程中能够满足飞行器的最大动压及过载限制。