[发明专利]一种基于三维再入弹道解析解的全射向自主再入制导方法

摘要

本发明一种基于三维再入弹道解析解的全射向自主再入制导方法，包括如下步骤：1、建立运动学与动力学方程；2、再入制导方法概述；3、下降段制导策略；4、平稳滑翔段制导策略；5、高度调整段制导策略。优点在于：该方法可以有效导引大升阻比飞行器的再入制导，能够有效抑制高升阻比引起的弱阻尼长周期弹道振荡，降低热流密度的峰值，有利于对参考轨迹的跟踪；该方法将由地球自转引起的惯性力与实际气动力组合为等效气动力，利用解析解规划等效气动剖面，使飞行器具有全向再入任务处理能力；该方法提出了基于反比例函数的等效气动剖面拟合公式，并进一步推导了基于复杂等效气动剖面形式的再入弹道解析解，具有高精度、全方位适应性和强鲁棒性。

主权项

1.一种基于三维再入弹道解析解的全射向自主再入制导方法，其特征在于：步骤1：建立运动学与动力学方程在研究再入制导问题时，用经度λ、纬度φ和高度H来描述飞行器位置信息，而用速率V、弹道倾角γ和航向角ψ来描述速度信息；考虑地球曲率和自转的影响，再入运动方程如下：其中，t为飞行时间，m是飞行器质量，为常值，kg；σ是倾侧角，rad；R_e是地球平均半径，大小为6356.766km；L是升力，D是阻力，g为重力加速度，ω_e为地球自转角速度；高超声速飞行器在飞行过程中需要满足热流密度来流动压q和机动过载约束n，如下所示：其中，ρ是大气密度，k_Q是与热流密度相关的常数，q_max是最大来流动压，是最大热流密度，n_max是最大机动过载，g₀是海平面高度的重力加速度；另外，考虑到飞控系统能力的限制，对攻角指令和倾侧角指令变化率进行约束；攻角指令约束：和倾侧角指令约束为：|σ|≤80deg、和当飞行器抵达一个以目标为中心50km为半径的圆上时，再入段终止；此时需要满足的终端条件是：期望的终点剩余飞行距离S_TAEM＝50km，期望的终端高度H_TAEM＝25km，期望的终端速度V_TAEM＝2000m/s，期望的终端航向误差|Δψ_TAEM|≤5deg和期望的终端倾侧角|σ_TAEM|≤30deg；这里及全文中下标“TAEM”代表期望的终端状态；步骤2：再入制导方法概述根据弹道特点，再入过程被划分为三个阶段：下降段、平稳滑翔阶段和高度调整阶段；在下降段，初始高度高，大气密度稀薄，飞行器快速掉高；在平稳滑翔阶段，由于大气密度提供足够的升力来平衡重力及离心力，滑翔高度平缓下降；在此阶段，利用三维再入弹道解析解实时解析规划等效气动参考剖面，并在数个时间点迭代调整倾侧反转点，以消除横程误差；在最后一次倾侧反转之后，此时飞行器距离目标足够近，再入过程进入高度调整阶段；此时飞行器通过调节攻角曲线来达到期望的飞行高度；下面步骤3‑5详细介绍各个阶段所采用的制导策略；步骤3：下降段制导策略由于初始高度过高，升力过小，飞行器快速掉高；随着飞行器扎进稠密大气层，热流密度会迅速增大，并在谷底达到局部峰值；为了保证热流密度约束，采用最保守的制导方案，其采用最大可用攻角、零倾侧角下滑，从而提升滑翔高度；当升力足以支撑飞行器平稳滑翔时，攻角由最大值平滑地过渡到基准攻角，随后进入平稳滑翔阶段；步骤4：平稳滑翔段制导策略平稳滑翔阶段的制导策略流程如下：1)设计一条基准攻角剖面，并确定与之对应的基准升阻比剖面；2)设计纵向升阻比参数化剖面，其外形与基准升阻比剖面相似；这样有利于近似保持倾侧角大小为常值；与此同时得到横向升阻比参数化剖面；气动剖面参数和需要在步骤3)与步骤4)中进行求解；3)为了补偿地球自转的影响，将惯性力与实际气动力组合为等效气动力；通过合理评估惯性力的影响，在气动剖面的基础上设计反比例函数形式的等效纵、横向升阻比曲线；4)推导基于反比例函数复杂形式等效气动剖面的新再入弹道解析解，然后根据纵程和终端能量要求，利用纵程解析解求解出步骤2)中的气动剖面参数和5)为了消除横程误差，利用横程解析确定倾侧反转点；为了减轻计算负担，仅在几个特殊时间节点修正反转点，而非实时修正；6)调节基准倾侧角跟踪步骤2)中纵向升阻比曲线，并实现倾侧反转；7)为了抑制弹道振荡，将弹道阻尼反馈引入到基准攻角和倾侧角中，从而得到制导指令；8)将过程约束转化为倾侧角约束，进而限制倾侧角指令大小，以便保证过程约束；9)从步骤2)开始重复上述过程，直到平稳滑翔阶段结束；假设C_E是地心，M代表飞行器，T代表目的地；在由此三点构成的平面内构造一个以C_E为圆心，地球平均半径R_e为半径的大圆GC1；飞行器的速度矢量记为与大圆GC1的夹角记为Δψ，以大圆切线在当地水平面内沿逆时针方向转向为正方向；因此，Δψ此时取负值；假设P₁是轨迹上的任意一点；过P₁点作一个新的大圆GC2，其垂直于大圆GC1；这两个大圆的交点记为P₂；那么定义P₁点的纵程为大圆弧MP₂的长度，而P₁点的横程为大圆弧P₁P₂的长度；定义L₁和L₂分别是升力在纵平面和当地水平面内的分量，即L₁＝Lcos(σ)和L₂＝Lsin(σ)；定义纵向升阻比为L₁与阻力D的比值，即而横向升阻比为L₂与D的比值，即L₂/D＝Lsin(σ)/D；定义能量E为：其中，μ是地球引力常数，为3.986005×10¹⁴m³/s²；以纵向升阻比和横向升阻比为控制量，开发了如下三维再入弹道解析解，其预测纵程x_D、横程x_C和航向角误差Δψ；这些解析解考虑了地球曲率的影响，但忽略了地球自转的影响；其中，E₀代表初始能量，E代表当前能量，x_E为积分变量，定义域为(E,E₀)；x_C0为初始横程，Δψ₀为初始航向角误差；R^*为常数，且有R^*＝R_e+H^*，其中H^*取飞行过程中高度的中间值；f₁(x_E)为简记函数表达式，如下：从解析解中看出：在忽略地球自转的情况下，纵程与纵向升阻比几乎成正比，而横程与横向升阻比几乎成正比；(1)基准攻角和基准升阻比设计以能量E为自变量，设计剖面函数；基准攻角α_bsl的形式如下：其中，α₁＝10deg近似对应最大升阻比攻角，以获得最远滑行距离，α₂＝6deg为基准攻角参数，E_α＝‑5.55×10⁷J/kg为临近平稳滑翔阶段和高度调整阶段的交班点，E_TAEM为再入过程中期望的终端能量；在有些扰动情况下，出现E＜E_TAEM的情况；此时，为了防止攻角出现负值，设定限制：α_bsl≥5deg；则画出基准攻角对应的升阻比曲线，称为基准升阻比曲线；(2)纵、横向升阻比设计为了合理耗散能量，并减轻飞控系统负担，在平稳滑翔阶段倾侧角要保持常值；按此要求，设计如下纵向升阻比曲线：其中，和是两个需要设计的待定升阻比参数，即气动剖面参数；由基准升阻比曲线看出，在E≥E_α时令纵向升阻比保持常值满足保持倾侧角大小为常值的要求；在后续步骤中会介绍如何根据剩余飞行距离确定的大小；当E＜E_α时，纵向升阻比要平滑地减小到但期望的终端倾侧角应近似为0，因此，令其中为补偿气拉偏影响的参数；是利用气动辨识技术测得的当前时刻实际升阻比，而是当前状态对应的理想升阻比；则得到参数化的横向升阻比模值曲线为(3)设计等效纵、横向升阻比并求解气动剖面参数将惯性力等效为三项附加气动力ΔL₁、ΔL₂和ΔD；根据公式(4)、(5)、(6)有：由于γ≈0，假设sin(γ)＝0和cos(γ)＝1，则式(18－20)被简化为：由于V是ω_e(R_e+H)cos(φ)的数倍、甚至数十倍，即2V＞＞ω_e(R_e+H)cos(φ)，则附加气动力ΔL₁和ΔL₂的公式进一步被简化为：ΔL₁＝2mVω_ecos(φ)sin(ψ)             (24)ΔL₂＝2mVω_esin(φ)              (25)定义等效纵向升阻比横向升阻比为上式的一阶Taylor展开为：从上式中看到，如果高度越高，大气密度越小，D越小，则会对和产生很大的影响；由于在大部分时间里飞行器平稳滑翔，所以D取平稳滑翔高度对应的阻力值，记为D_SG；则有忽略地球自转的影响，并且在平稳滑翔状态下假设γ≈0和进而由公式(5)和(10)得：其中，R^*＝R_e+H；利用纵向升阻比计算出平稳滑翔高度对应的阻力值D_SG，如下：其中，L_1(SG)为平稳滑翔高度对应的升力值；下面分析等效纵向升阻比和等效横向升阻比的变化趋势，并进行适当的拟合；将公式(23–25)和(33)代入到公式(30–31)中，并整理得到：其中，其中，其中，h₁,h₂,h₃,h_z1,h_z2,h_z3,h_m为拟合参数；为了应用解析解，需要将h₁–h₃拟合成关于能量的表达式；由于分母已经是E的函数，这里只需要拟合上述公式的分子；通过观察仿真结果，发现cos(φ)sin(ψ)近似保持常值；在公式(23－25)中，V和sin(φ)分别采用过两端点的线性函数来拟合，这表明Vsin(φ)适用利用二次多项式来拟合；由于2V＞＞ω_e(R_e+H)cos(φ)，附加气动力ΔD要远远小于附加气动力ΔL₁和ΔL₂；因此，尽管附加气动力ΔD变化无常，但是其影响较小，利用线性公式来近似附加气动力ΔD；因此，利用线性公式拟合ΔL₁和ΔD，而利用二次曲线拟合ΔL₂；根据上述分析，采用线性公式拟合h_z1和h_z2，而利用二次多项式拟合h_z3，如下：其中，x_E为积分变量，物理意义为飞行器能量，上划线“‑”表示此变量为拟合值，k_h1(0)、k_h1(1)、k_h2(0)、k_h2(1)、k_h3(0)、k_h3(1)、k_h3(2)为相应拟合参数系数；和的参数满足这样的条件：拟合线段经过原函数的两个端点，第一个端点由当前状态确定，第二个端点由期望的终端状态确定；下面说明拟合参数系数的计算方法；其中，其中，是上一制导周期获得的当前纵向升阻比，是期望的终端升阻比，φ_TAEM是目标的纬度，是飞行器的航向角ψ在某种意义下的加权平均值，是飞行器的航向角ψ在终点处的加权平均值；根据公式(25)和(39)，通过分别拟合V和sin(φ)，得到的表达式：展开上式右得拟合公式参数：接下来，利用纵程解析解确定的大小，以满足纵程和终端能量要求；根据公式(16)和(34)，规划的等效纵向升阻比曲线记为如下形式：将式(11)中的纵向升阻比用上式的等效升阻比替换，并积分得基于反比例函数气动剖面的纵程解析表达式；由于是分段函数，根据积分前后能量E₁，E₂与交班点能量E_α的大小关系，下面需要考虑如下三种情况：1)当E_α≤E₂≤E₁时，将公式(52)中的第一个子式代入公式(11)，并积分得到纵程为：2)当E₂≤E₁≤E_α时，将公式(52)中的第二个子式代入公式(11)，并积分得到纵程为：其中，a₀,a₁,a₂为简记符号，其值如下：3)当E₂≤E_α≤E₁时，则纵程则由上面两个公式组合计算，如下：x_D(E₂,E₁)＝x_D(E_α,E₁)+x_D(E₂,E_α)         (58)为了计算待定升阻比参数需要首先计算剩余飞行距离s_go；建立地心赤道旋转坐标系GER C_E‑x_ey_ez_e，其原点在地心C_E，x_e和y_e轴在赤道平面内，并且x_e轴指向本初子午线，z_e轴指向北极；M点代表飞行器，N点为M点向赤道平面的垂直投影点，T代表目标；剩余飞行距离就是过M和T的大圆弧长度；定义为从C_E到M的向量，定义为从C_E到T的向量；剩余飞行距离由下式计算：s_go＝R_eη             (59)其中η是和之间的夹角；下面求解η；记x_EM是的单位向量，x_ET是的单位向量；由于与赤道平面的夹角刚好就是纬度φ，x_EM在赤道平面内的分量大小为cos(φ)，而沿z_e轴分量大小为sin(φ)；进而，利用线段C_EN与x_e轴的夹角为经度λ，得到单位向量x_EM沿x_e轴的分量为cos(λ)cos(φ)，沿y_e轴的分量为sin(λ)cos(φ)；因此，x_EM在地心赤道旋转坐标系GER中的坐标为其中，上标“GER”代表以GER系为坐标系，而上标“T”代表向量转置；同理得到x_ET在GER坐标系中的坐标其中，λ_T和φ_T分别是目标点处的经度和纬度；两个单位向量的点积刚好等于它们之间夹角的余弦；因此，η由下式计算：接下来，利用剩余飞行距离s_go确定待定升阻比参数由于等效纵向升阻比曲线是分段函数，下面需要分两种情况考虑：1)当E＞E_α时，有如下等式关系：x_D(E_TAEM,E_α)+x_D(E_α,E)＝s_go‑S_TAEM           (63)其中，x_D(E_TAEM,E_α)由公式(54)计算，而x_D(E_α,E)由公式(53)计算；通过求解公式(63)，得到其中，k_xD1，k_xD2为简记符号，其值如下：2)当E≤E_α时此时飞行处于高度调整阶段，将采用其它制导策略；因此，无需再更新L₁/D₁；(4)倾侧反转点规划在飞行中，飞行器需要通过调节倾侧角来跟踪纵向升阻比曲线，但是这会导致大幅横向机动；为了消除由此引起的横程误差，飞行器需要进行倾侧反转；该制导方法仅安排两次反转，以减轻飞控系统负担；记倾侧反转对应的能量节点分别为E_BR1和E_BR2；在第一次反转发生之前，令E_BR2＝E_α，然后利用横程解析解来规划E_BR1；在第一次反转发生之后，采用一种基于在线弹道仿真的迭代修正算法微调E_BR2；在确定升阻比曲线之后，考虑到倾侧反转，利用公式(35)，等效横向升阻比曲线被表述为：其中，sgn是一个取值±1的符号函数，用于确定初始倾侧方向；这里sgn的取值使得飞行器在初始时刻向航向误差减小的方向倾侧；k_BR代表已经发生反转的次数，由公式(17)计算；记飞行器当前能量为E，此时飞行器的横程为x_C，航向角误差为Δψ，将公式(67)代入到横程解析解中即公式(12)，并令飞行器能量从E变化到E_TEAM，则预测飞行器的终端横程x_Cf为：其中，x_E为积分变量；x_D(E_TAEM,E)为飞行器能量从E积分到E_TEAM时的纵程，x_D(E_TAEM,x_E)为飞行器能量从x_E积分到E_TEAM时的纵程，这两项由公式(11)计算得到；F(E_BR1,E)、F(E_BR2,E_BR1)、F(E_TAEM,E_BR2)为函数F(x_E2,x_E1)分别在定义域(E_BR1,E)、(E_BR2,E_BR1)、(E_TAEM,E_BR2)上的积分值；函数F(x_E2,x_E1)的定义为：其中简记函数表达式f₂(x_E)的形式如下：此外，由于规划的等效纵向升阻比曲线是分段函数，则纵程x_D(E_TAEM,x_E)同样是分段函数，如下所示：其中，x_D(E_TAEM,E_α),x_D(E_TAEM,x_E)由公式(54)计算，x_D(E_α,x_E)由公式(53)计算；由于期望终端横程x_Cf＝0，第一次倾侧反转对应的能量节点E_BR1的数值通过牛顿法迭代求解方程x_Cf(E_BR1)＝0得到，如下；其中，上标“(k)”、“(k+1)”为能量节点E_BR1的迭代次数，x_Cf对E_BR1的导数x′_Cf为由于x_Cf随E_BR1单调变化，只需要迭代数次就满足终端横程为了减轻弹载计算机计算负担，并不实时修正E_BR1，而是仅在第一次反转发生之前的几个适当时间节点进行修正；(5)基准倾侧角设计飞行器需要通过调节基准倾侧角来跟踪纵向升阻比曲线，并实现倾侧反转；基准倾侧角σ_bsl如下：其中，是当前时刻利用惯导气动数据辨识的实际升阻比；ΔE是考虑飞行器的滚转速率限幅引入的补偿量，如下：其中，是自动驾驶仪限制的最大倾侧角变化率，a_D是当前的阻力加速度，利用含加速度计的惯导组件测量；(6)生成攻角指令和倾侧角指令若直接采用公式(15)的基准攻角α_bsl和公式(74)的基准倾侧角作为指令攻角α_cmd和σ_cmd，则滑翔弹道会存在弱阻尼、长周期振荡；因此，采用弹道阻尼技术来抑制弹道振荡；攻角指令α_cmd和倾侧角指令σ_cmd为：α_cmd＝α_bsl+cos(σ_bsl)k_γ(γ_SG‑γ)               (76)其中，k_γ是增益系数，α₁是基准攻角剖面的一个参数，γ_SG为平稳滑翔弹道对应的弹道倾角；γ_SG的高精度近似公式如下：其中，D_bsl＝C_D(bsl)qS_ref是α_bsl对应的阻力，C_D(bsl)是α_bsl对应的阻力系数，ρ是大气密度，q是动压，S_ref是气动参考面积；系数d₁和d₂的计算公式如下：其中，C_L(bsl)是基准攻角α_bsl对应的升力系数；dC_L(bsl)/dE提前计算好，d[cos(σ_bsl)]/dE利用差分估算；为了保证热流密度等过程约束，将过程约束转化成倾侧角指令范围限制；对于平稳滑翔弹道，大的倾侧角会导致高度降低加剧，从而大气密度增大，引起热流密度、来流动压及过载的增大；最大的用倾侧角发生在滑翔高度下限H_min处，此时大气密度达到最大，升力也达到最大，升力在垂直方向上的分量近似满足平稳滑翔条件；因此平稳滑翔状态下的最大用倾侧角由公式(81)右边第一项计算；当有弹道振荡时，如果飞行器快速掉高，需要减小σ_max，以防止飞行器掉到H_min以下；此时，添加反馈补偿，即公式(81)右边第二项；其中，L₁是升力在纵平面内的分量，L_max为达到的最大升力，参数k_σ＝‑50，且有其中，ΔL₁由公式(21)计算，是利用气动辨识技术测量得到的当前实际升力系数；公式如下：步骤5：高度调整段制导策略(1)基准攻角和最后一次反转点修正在进入高度调整阶段之前，需要利用基于弹道仿真的迭代修正算法，对基准攻角参数α₂和第二次反转点的能量E_BR2进行微调，以满足终端速度和高度约束；为了减轻计算负担，在第二次反转发生之前仅进行两轮修正：第一轮发生在第一次反转之后不久；第二轮大约发生在第二次反转发生之前的60s；首先介绍在线弹道仿真；在第n次弹道仿真中，首先根据已测的气动拉偏信息对气动模型粗略修正，然后根据上一步的迭代结果更新第二次反转点能量和基准攻角参数，即和最后以再入制导方法得到的控制律，对运动方程公式(1–6)进行数值积分，直到为止；上标“(n)”表示此变量为第n次弹道仿真的值，n＝1,2,3,…；迭代修正算法包含两个子算法：基准攻角修正算法和反转点修正算法；首先，调用基准攻角修正算法修正基准攻角曲线，以便提高终端高度精度；此任务仅需要执行一次弹道仿真；接着，利用反转点修正算法精确微调最后一次反转点E_BR2，以便满足终端速度约束；这两种算法的详细介绍如下；1)基准攻角修正算法通过执行一次弹道仿真，预测包括终端高度在内的终端状态，然后，通过比较和H_TAEM来修正基准攻角参数α₂，以提高终端高度精度；注意区分下标“f”和“TAEM”，“f”表示此变量为仿真所获得的终端状态，“TAEM”表示此变量为期望的终端状态；具体做法如下：在执行完一次弹道仿真后，由于和近似于0，假设和并且假设然后利用公式(5)得到其中，根据气动辨识技术测量到的当前气动拉偏信息所修正的终点升力系数；飞行器期望的终端状态应当满足其中，表示第二次弹道仿真时应当使用的基准攻角，即第一次弹道仿真的基准攻角修正之后的值；由于滑翔高度远远小于地球半径，假设另外，还假设然后将公式(85)、(86)两式相减，并经过代数变换，得到：线化得到：其中，且就是第一次弹道仿真后获得的应当对基准攻角参数进行的修正量；是在考虑气动拉偏修正之后的升力线斜率；通过求解公式(87–88)，得到：从而，修正之后的基准攻角参数为此值将作为下一次弹道仿真中的基准攻角参数；2)反转点的修正算法利用反转点修正算法微调E_BR2，以便满足终端速度要求；首先定性分析E_BR2与V_f的关系：若E_BR2减小，最后一次反转推迟，进而在最后一次反转时刻的航向角误差增大；因此，在高度调整阶段，比例导引律会产生更大的倾侧角来消除航向角误差；这使得纵向升阻比减小，从而导致终点速度V_f减小；首先利用弹道仿真预测终点速度然后利用secant法修正E_BR2，如下：重复上述过程，直到(2)高度调整阶段的基准倾侧角设计此阶段，采用比例导引进行基准倾侧角的设计；视线方位角变化率为：比例导引律产生的需用横向机动加速度a_L2为：其中，ΔL₂是由公式(22)计算；上式右边第二项是为了补偿地球自转的影响；为了防止初始倾侧角饱和，这里令k_PN从2逐渐变化到4，如下另一方面，在平稳滑翔情况下，纵向平面内升力加速度a_L1与重力、离心力和惯性力平衡，即：其中ΔL₁是由公式(21)计算；则高度调整阶段的基准倾侧角为：(3)生成高度调整阶段的攻角和倾侧角指令高度调整阶段的攻角和倾侧角指令如下：其中，是剩余飞行距离—能量参考曲线，由步骤5(1)中的在线弹道仿真获得；k_α＝5π/(1.8×10⁷)，其意味着当实际剩余飞行器距离偏离参考曲线1km时，对攻角修正0.05度；在此阶段，各种干扰仍然会引起速度误差；由于调节攻角也能调节升阻比，为了保证终端速度约束，通过调节攻角来跟踪由于此阶段的航程短，干扰引起的速度误差小，跟踪并不会引起攻角的大幅调整；公式(97)右边第二项用来抑制弹道振荡；为了保证热流密度等过程约束，这里亦采取如下修正措施：其中，由公式(81)计算确定。