CN108036676A - 一种基于三维再入弹道解析解的全射向自主再入制导方法 - Google Patents

一种基于三维再入弹道解析解的全射向自主再入制导方法 Download PDF

Info

Publication number
CN108036676A
CN108036676A CN201711260656.7A CN201711260656A CN108036676A CN 108036676 A CN108036676 A CN 108036676A CN 201711260656 A CN201711260656 A CN 201711260656A CN 108036676 A CN108036676 A CN 108036676A
Authority
CN
China
Prior art keywords
angle
formula
attack
resistance ratio
lift
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201711260656.7A
Other languages
English (en)
Other versions
CN108036676B (zh
Inventor
余文斌
陈万春
赵鹏雷
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN201711260656.7A priority Critical patent/CN108036676B/zh
Publication of CN108036676A publication Critical patent/CN108036676A/zh
Application granted granted Critical
Publication of CN108036676B publication Critical patent/CN108036676B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F41—WEAPONS
    • F41G—WEAPON SIGHTS; AIMING
    • F41G3/00—Aiming or laying means
    • F41G3/22—Aiming or laying means for vehicle-borne armament, e.g. on aircraft
    • F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
    • F41—WEAPONS
    • F41G—WEAPON SIGHTS; AIMING
    • F41G3/00—Aiming or laying means
    • F41G3/08—Aiming or laying means with means for compensating for speed, direction, temperature, pressure, or humidity of the atmosphere
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Automation & Control Theory (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Traffic Control Systems (AREA)

Abstract

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

Description

一种基于三维再入弹道解析解的全射向自主再入制导方法
技术领域
本发明涉及一种基于三维再入弹道解析解的全射向自主再入制导方法，属于航天技术、 武器技术、制导控制领域。
背景技术
CAV是一种高升阻比的高超声速滑翔飞行器，其升阻比高达3。相比之下，常规再入飞 行器的升阻比大约仅为1，例如航天飞机和X-33RLV等。高升阻比使得CAV在再入阶段拥有更远的滑行距离和更广泛的横向机动范围。在再入阶段内，飞行器的滑翔高度在 20km-100km之间，马赫数由20逐渐衰减至6附近。不同于寻的制导方法，用于此阶段的再 入制导方法不但要将飞行器导引向目标，还需要通过执行适当的横向机动来合理耗散飞行能 量。另外，为了满足热防护系统和飞行器结构强度的限制，再入轨迹需要满足热流密度、来 流动压和可用过载等约束。
航天飞机再入制导方法是最早得到成功应用的适用于滑翔式飞行器的再入制导方法。此 制导方法成功的关键在于采用倾侧角而非攻角作为主要调节手段来管理飞行能量，因为气动 阻力与倾侧角之间存在明确、单调的关系。此制导方法将制导流程划分为四部分：1、确定满 足所有过程约束的再入走廊；2、在走廊内规划一条满足终端条件的参考轨迹；3、调节倾侧 角跟踪参考轨迹；4、设计倾侧反转门限，以便消除航向误差。这四个部分构成了所有再入制 导方法的基本流程。然而，虽然此传统再入制导方法能够可靠导引常规再入飞行器，例如航 天飞机和X-33RLV等(升阻比约为1)，但是无法有效导引CAV。这是因为后者的升阻比远 高于前者，而高升阻比会引起弱阻尼长周期弹道振荡。这不仅导致热流密度的峰值过大，也 不利于对参考轨迹的跟踪，特别是存在大干扰的情况下。传统再入制导方法并不能有效消除 弹道振荡。
Lu P.,Entry guidance:a unified method[J].Journal of Guidance Control&Dynamics,2014, 37(3):713-728将弹道阻尼控制技术与预测校正制导律结合，开发了一种能够可靠导引CAV的 再入制导方法，简记为PCGG，但制导精度较低，终端航向误差仍然会轻易超出要求，有时 甚至高达20deg。在本发明的以下说明中简称为文献一。
Yu W.,Chen W.,Entry guidance with real-time planning of referencebased on analytical solutions[J].Advances in Space Research,2015,55(9):2325–2345开发了三维再入弹道解析解， 然后利用此解析解设计了一种能够自行规划参考弹道的再入制导方法，记为EGAS，不过， 经过深入研究发现，由于三维再入弹道解析解没有考虑地球自转的影响，而弹道规划又没有 对此影响进行充分补偿，EGAS并不能应对所有飞行方向的任务，特别是沿南北方向的任务。 在本发明的以下说明中简称为文献二。
发明内容
本发明的目的是为了解决上述问题，提出一种基于三维再入弹道解析解的全射向自主再 入制导方法。其中，纵程解析解被用于在线自主规划纵向参考剖面，而横程解析解被用于规 划反转点。这些解析解考虑了地球曲率的影响，但是忽略了地球自转的影响。为了补偿地球 自转的影响，这里将由地球自转引起的惯性力与实际气动力组合为等效气动力，并利用解析 解规划等效气动剖面。为了避免攻角和倾侧角出现大幅剧烈调整，通过合理分析强非线性的 惯性力剖面变化规律，提取了其中近似为常值的项，并提出了基于反比例函数的等效气动剖 面拟合公式。进一步，开发了基于复杂等效气动剖面形式的再入弹道解析解，并给出剖面参 数解算方法。实施例验证了此制导方法具有高精度、全方位适应性和强鲁棒性。
步骤1：建立运动学与动力学方程
在研究再入制导问题时，常用经度λ、纬度φ和高度H来描述飞行器位置信息，而用速 率V、弹道倾角γ和航向角ψ来描述速度信息。考虑地球曲率和自转的影响，基于上述状态 变量的再入运动方程如下
其中，t为飞行时间，m是飞行器质量，为常值，kg。σ是倾侧角，rad。Re是地球平均半径，大小为6356.766km。L是升力，D是阻力，g为重力加速度，ωe为地球自转角速度。
一般，高超声速飞行器在飞行过程中需要满足热流密度来流动压q和机动过载约束n， 如下所示
其中，ρ是大气密度，kQ是与热流密度相关的常数，qmax是最大来流动压，是最大热流 密度，nmax是最大机动过载，g0是海平面高度的重力加速度。
另外，考虑到飞控系统能力的限制，对攻角指令和倾侧角指令变化率进行约束。攻角指 令约束：和倾侧角指令约束为：|σ|≤80deg、和
当飞行器抵达一个以目标为中心50km为半径的圆上时，再入段终止。此时需要满足的 终端条件是：期望的终点剩余飞行距离STAEM＝50km，期望的终端高度HTAEM＝25km，期望的终端速度VTAEM＝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)开始重复上述过程，直到平稳滑翔阶段结束。
由于合理补偿地球自转的影响，上述策略可以保证攻角和倾侧角在此阶段近似为常值， 这有三个好处：(1)合理耗散飞行能量，避免制导指令在飞行后期饱和；(2)有利于飞行控 制系统精确跟踪指令；(3)由于气动拉偏系数对攻角敏感，常值攻角有利于准确评估气动拉 偏的影响，从而提高制导精度和鲁棒性。下面详细介绍制导策略。
假设CE是地心，M代表飞行器，T代表目的地。在由此三点构成的平面内构造一个以CE为圆心，地球平均半径Re为半径的大圆GC1。飞行器的速度矢量记为V。V与大圆GC1 的夹角记为Δψ，以大圆切线在当地水平面内沿逆时针方向转向V为正方向。因此，Δψ此 时取负值。假设P1是轨迹上的任意一点。过P1点作一个新的大圆GC2，其垂直于大圆GC1。 这两个大圆的交点记为P2。那么定义P1点的纵程为大圆弧的长度，而P1点的横程为大 圆弧的长度。
定义L1和L2分别是升力在纵平面和当地水平面内的分量，即L1＝Lcos(σ)和 L2＝Lsin(σ)。定义纵向升阻比为L1与阻力D的比值，即而横向升阻比 为L2与D的比值，即L2/D＝Lsin(σ)/D。定义能量E为
其中，μ是地球引力常数，为3.986005×1014m3/s2。
文献二以纵向升阻比和横向升阻比为控制量，开发了如下三维再入弹道解析 解，其可以预测纵程xD、横程xC和航向角误差Δψ。这些解析解考虑了地球曲率的影响，但 忽略了地球自转的影响。
其中，E0代表初始能量，E代表当前能量，xE为积分变量，定义域为(E,E0)。xC0为初始横 程，Δψ0为初始航向角误差。R*为常数，且有R*＝Re+H*，其中H*一般取飞行过程中高度 的中间值。f1(xE)为简记函数表达式，如下
从解析解中可以看出：在忽略地球自转的情况下，纵程与纵向升阻比几乎成正比，而横 程与横向升阻比几乎成正比。
(1)基准攻角和基准升阻比设计
以能量E为自变量，设计剖面函数。基准攻角αbsl的形式如下
其中，α1＝10deg近似对应最大升阻比攻角，以获得最远滑行距离，α2＝6deg为基准攻 角参数，Eα＝-5.55×107J/kg为临近平稳滑翔阶段和高度调整阶段的交班点，ETAEM为再入过 程中期望的终端能量。在有些扰动情况下，有可能出现E＜ETAEM的情况。此时，为了防止攻 角出现负值，设定限制：αbsl≥5deg。则可画出基准攻角对应的升阻比曲线，称为基准升阻 比曲线
(2)纵、横向升阻比设计
为了合理耗散能量，并减轻飞控系统负担，在平稳滑翔阶段倾侧角要尽可能保持常值。 按此要求，设计如下纵向升阻比曲线
其中，和是两个需要设计的待定升阻比参数，即气动剖面参数。由基准升阻 比曲线可以看出，在E≥Eα时令纵向升阻比保持常值可满足尽可能保持倾侧角大 小为常值的要求。在后续步骤中会介绍如何根据剩余飞行距离确定的大小；当E＜Eα 时，纵向升阻比要平滑地减小到但期望的终端倾侧角应近似为0，因此，
令其中为补偿气拉偏影响的参 数。是利用气动辨识技术测得的当前时刻实际升阻比，而是当前状态对 应的理想升阻比。则可得参数化的横向升阻比模值曲线为
(3)设计等效纵、横向升阻比并求解气动剖面参数
将惯性力等效为三项附加气动力ΔL1、ΔL2和ΔD。根据公式(4)(5)(6)有。
由于γ≈0，假设sin(γ)＝0和cos(γ)＝1，则式(18－20)可被简化为
由于V一般是ωe(Re+H)cos(φ)的数倍、甚至数十倍，即2V>>ωe(Re+H)cos(φ)，则附加 气动力ΔL1和ΔL2的公式可以进一步被简化为
ΔL1＝2mVωecos(φ)sin(ψ) (24)
ΔL2＝2mVωesin(φ) (25)
定义等效纵向升阻比横向升阻比为
上式的一阶Taylor展开为
从上式中可以看到，如果高度越高，大气密度越小，D越小，则会对和产生很 大的影响。由于在大部分时间里飞行器平稳滑翔，所以D取平稳滑翔高度对应的阻力值，记 为DSG。则有
忽略地球自转的影响，并且在平稳滑翔状态下假设γ≈0和进而由公式(5)和(10)得
其中，R*＝Re+H。
利用纵向升阻比可以计算出平稳滑翔高度对应的阻力值DSG，如下
其中，L1(SG)为平稳滑翔高度对应的升力值。下面分析等效纵向升阻比和等效横向升 阻比的变化趋势，并进行适当的拟合。将公式(23–25)和(33)代入到公式(30–31) 中，并整理可得
其中
其中
其中，h1,h2,h3,hz1,hz2,hz3,hm为拟合参数。
为了应用解析解，需要将h1–h3拟合成关于能量的表达式。由于分母已经是E的函数，这 里只需要拟合上述公式的分子。
通过观察仿真结果，发现cos(φ)sin(ψ)近似保持常值。在公式(23－25)中，V和sin(φ)分 别采用过两端点的线性函数来拟合，这表明Vsin(φ)适用利用二次多项式来拟合。
由于2V>>ωe(Re+H)cos(φ)，附加气动力ΔD(公式(23))要远远小于附加气动力ΔL1和 ΔL2。因此，尽管附加气动力ΔD变化无常，但是其影响较小，本发明粗略地利用线性公式来 近似附加气动力ΔD。因此，本发明利用线性公式拟合ΔL1和ΔD，而利用二次曲线拟合ΔL2。 根据上述分析，本发明采用线性公式拟合hz1和hz2，而利用二次多项式拟合hz3，如下
其中，xE为积分变量，物理意义为飞行器能量，上划线“-”表示此变量为拟合值，kh1(0)、 kh1(1)、kh2(0)、kh2(1)、kh3(0)、kh3(1)、kh3(2)为相应拟合参数系数。和的参数满足这样的条 件：拟合线段经过原函数的两个端点，第一个端点由当前状态确定，第二个端点由期望的终 端状态确定。下面说明拟合参数系数的计算方法。
其中，
其中，是上一制导周期获得的当前纵向升阻比，是期望的终端 升阻比，φTAEM是目标的纬度，是飞行器的航向角ψ在某种意义下的加权平均值，是 飞行器的航向角ψ在终点处的加权平均值。
根据公式(25)和(39)，通过分别拟合V和sin(φ)，可以得到的表达式
展开上式右得拟合公式参数
接下来，利用纵程解析解确定的大小，以满足纵程和终端能量要求。根据公式(16) 和(34)，规划的等效纵向升阻比曲线可记为如下形式
将式(11)中的纵向升阻比用上式的等效升阻比替换，并积分可得基于反比例 函数气动剖面的纵程解析表达式。由于是分段函数，根据积分前后能量E1，E2与交班点 能量Eα的大小关系，下面需要考虑如下三种情况
1)当Eα≤E2≤E1时，将公式(52)中的第一个子式代入公式(11)，并积分可得纵程为
2)当E2≤E1≤Eα时，将公式(52)中的第二个子式代入公式(11)，并积分可得纵程为
其中，a0,a1,a2为简记符号，其值如下
3)当E2≤Eα≤E1时，则纵程则可以由上面两个公式组合计算，如下
xD(E2,E1)＝xD(Eα,E1)+xD(E2,Eα) (58)
为了计算待定升阻比参数需要首先计算剩余飞行距离sgo。建立地心赤道旋转坐 标系CE-xeyeze(GER)，其原点在地心CE，xe和ye轴在赤道平面内，并且xe轴指向本初子午线，ze轴指向北极。M点代表飞行器，N点为M点向赤道平面的垂直投影点，T代表目 标。剩余飞行距离就是过M和T的大圆弧长度。定义为从CE到M的向量，定义为 从CE到T的向量。剩余飞行距离由下式计算
sgo＝Reη (59)
其中η是和之间的夹角。下面求解η。
记xEM是的单位向量，xET是的单位向量。由于与赤道平面的夹角刚好就 是纬度φ，xEM在赤道平面内的分量大小为cos(φ)，而沿ze轴分量大小为sin(φ)。进而，利用 线段CEN与xe轴的夹角为经度λ，可以得到单位向量xEM沿xe轴的分量为cos(λ)cos(φ)，沿ye轴的分量为sin(λ)cos(φ)。因此，xEM在地心赤道旋转坐标系GER中的坐标为
其中，上标“GER”代表以GER系为坐标系，而上标“T”代表向量转置。同理可以得到xET在GER坐标系中的坐标
其中，λT和φT分别是目标点处的经度和纬度。从几何意义上，两个单位向量的点积刚好 等于它们之间夹角的余弦。因此，η可以由下式计算
接下来，利用剩余飞行距离sgo可以确定待定升阻比参数由于等效纵向升阻比曲 线是分段函数，下面需要分两种情况考虑：
1)当E＞Eα时，有如下等式关系
xD(ETAEM,Eα)+xD(Eα,E)＝sgo-STAEM (63)
其中，xD(ETAEM,Eα)由公式(54)计算，而xD(Eα,E)由公式(53)计算。通过求解公式(63)，可以得到
其中，kxD1，kxD2为简记符号，其值如下
2)当E≤Eα时
此时飞行主要处于高度调整阶段，将采用其它制导策略。因此，无需再更新
(4)倾侧反转点规划
在飞行中，飞行器需要通过调节倾侧角来跟踪纵向升阻比曲线，但是这会导致大幅横向 机动。为了消除由此引起的横程误差，飞行器需要在适当的时候进行倾侧反转。本发明仅安 排两次反转，以减轻飞控系统负担。记倾侧反转对应的能量节点分别为EBR1和EBR2。在第一 次反转发生之前，令EBR2＝Eα，然后利用横程解析解来规划EBR1。在第一次反转发生之后， 采用一种基于在线弹道仿真的迭代修正算法微调EBR2。
在确定升阻比曲线之后，考虑到倾侧反转，利用公式(35)，等效横向升阻比曲 线可以被表述为
其中，sgn是一个取值±1的符号函数，用于确定初始倾侧方向。这里sgn的取值使得飞 行器在初始时刻向航向误差减小的方向倾侧。kBR代表已经发生反转的次数，由公式 (17)计算。
记飞行器当前能量为E，此时飞行器的横程为xC，航向角误差为Δψ，将公式(67)代入到横程解析解中(公式(12))中，并令飞行器能量从E变化到ETEAM，则可预测飞行器的 终端横程xCf为
其中，xE为积分变量。xD(ETAEM,E)为飞行器能量从E积分到ETEAM时的纵程，xD(ETAEM,xE)为 飞行器能量从xE积分到ETEAM时的纵程，这两项可由公式(11)计算得到。F(EBR1,E)、 F(EBR2,EBR1)、F(ETAEM,EBR2)为函数F(xE2,xE1)分别在定义域(EBR1,E)、(EBR2,EBR1)、 (ETAEM,EBR2)上的积分值。
函数F(xE2,xE1)的定义为
其中简记函数表达式f2(xE)的形式如下
此外，由于规划的等效纵向升阻比曲线是分段函数，则纵程xD(ETAEM,xE)同样是分 段函数，如下所示
其中，xD(ETAEM,Eα),xD(ETAEM,xE)可由公式(54)计算，xD(Eα,xE)可由公式(53)计算。
由于期望终端横程xCf＝0，第一次倾侧反转对应的能量节点EBR1的数值可通过牛顿法迭 代求解方程xCf(EBR1)＝0得到，如下。
其中，上标“(k)”、“(k+1)”为能量节点EBR1的迭代次数，xCf对EBR1的导数x′Cf为
由于xCf随EBR1单调变化，一般只需要迭代数次就可以满足终端横程为 了减轻弹载计算机计算负担，本发明并不实时修正EBR1，而是仅在第一次反转发生之前的几 个适当时间节点进行修正。
(5)基准倾侧角设计
飞行器需要通过调节基准倾侧角来跟踪纵向升阻比曲线，并实现倾侧反转。基准倾 侧角σbsl如下
其中，是当前时刻利用惯导气动数据辨识的实际升阻比。ΔE是考虑飞行器的滚转 速率限幅引入的补偿量，如下
其中，是自动驾驶仪限制的最大倾侧角变化率，aD是当前的阻力加速度，可以利用 含加速度计的惯导组件测量。
(6)生成攻角指令和倾侧角指令
若直接采用公式(15)的基准攻角αbsl和公式(74)的基准倾侧角作为指令攻角αcmd和σcmd， 则滑翔弹道会存在弱阻尼、长周期振荡。因此，本发明采用弹道阻尼技术来抑制弹道振荡。 攻角指令αcmd和倾侧角指令σcmd为
αcmd＝αbsl+cos(σbsl)kγ(γSG-γ) (76)
其中，kγ是增益系数，α1是基准攻角剖面的一个参数，γSG为平稳滑翔弹道对应的弹道 倾角。γSG的高精度近似公式如下
其中，Dbsl＝CD(bsl)qSref是αbsl对应的阻力，CD(bsl)是αbsl对应的阻力系数，ρ是大气密度， q是动压，Sref是气动参考面积。系数d1和d2的计算公式如下
其中，CL(bsl)是基准攻角αbsl对应的升力系数。dCL(bsl)/dE可提前计算好，d[cos(σbsl)]/dE 可利用差分估算。
为了保证热流密度等过程约束，本发明将过程约束转化成倾侧角指令范围限制。对于平 稳滑翔弹道，大的倾侧角会导致高度降低加剧，从而大气密度增大，引起热流密度、来流动 压及过载的增大。最大的可用倾侧角发生在滑翔高度下限Hmin处，此时大气密度达到最大， 升力也达到最大，升力在垂直方向上的分量近似满足平稳滑翔条件。因此平稳滑翔状态下的 最大可用倾侧角可由公式(81)右边第一项计算。当有弹道振荡时，如果飞行器快速掉高， 需要减小σmax，以防止飞行器掉到Hmin以下。此时，添加反馈补偿，即公式(81)右边第二 项。
其中，L1是升力在纵平面内的分量，Lmax为可达到的最大升力，参数kσ＝-50，且有
其中，ΔL1由公式(21)计算，是利用气动辨识技术测量得到的当前实际升力系数。 进而采取如下措施
步骤5：高度调整段制导策略
在最后一次反转发生之后，飞行器进入高度调整阶段。在此阶段，为了达到期望的终端 高度，飞行器的基准攻角逐渐被调整至α2(在公式(15)中定义)。在此阶段，由于弹道倾 角变化率较大，再入滑翔解析解精度不能满足终端要求。为此，这里本发明采用基于比例导 引律的制导策略，以消除终端航向误差。但是，不再跟踪纵向升阻比曲线会造成终端速度误 差。经理论分析可知，终点速度随最后一次反转点单调变化。为了满足终端速度要求，本发 明采取这样一个方案：在最后一次反转发生之前，利用基于弹道仿真的修正算法对第二次反 转点进行修正，以消除反转之前积累的速度误差；在最后一次反转之后，通过调节攻角来跟 踪剩余飞行距离—能量参考曲线，以克服干扰。由于这里仅对剩余一小段弹道进行数次仿真， 而非是实时仿真，所以并不会造成较大的计算负担。另外，修正算法还会对基准攻角曲线进 行略微调整，以提高终端高度精度。
(1)基准攻角和最后一次反转点修正
在进入高度调整阶段之前，需要利用基于弹道仿真的迭代修正算法，对基准攻角参数α2和第二次反转点的能量EBR2进行微调，以满足终端速度和高度约束。为了减轻计算负担，本 发明在第二次反转发生之前仅进行两轮修正：第一轮发生在第一次反转之后不久；第二轮大 约发生在第二次反转发生之前的60s。
首先介绍在线弹道仿真。在第n次弹道仿真中，首先根据已测的气动拉偏信息对气动模 型粗略修正，然后根据上一步的迭代结果更新第二次反转点能量和基准攻角参数，即和 最后以本发明的再入制导方法得到的控制律，对运动方程公式(1–6)进行数值积分， 直到为止。这里及下文，上标“(n)”表示此变量为第n次弹道仿真的值，n＝1,2,3,…。
迭代修正算法包含两个子算法：基准攻角修正算法和反转点修正算法。首先，调用基准 攻角修正算法修正基准攻角曲线，以便提高终端高度精度。此任务仅需要执行一次弹道仿真。 接着，本发明利用反转点修正算法精确微调最后一次反转点EBR2，以便满足终端速度约束。 这两种算法的详细介绍如下。
1)基准攻角修正算法
此算法通过执行一次弹道仿真，预测包括终端高度在内的终端状态，然后，通过比 较和HTAEM来修正基准攻角参数α2，以提高终端高度精度。注意区分下标“f”和“TAEM”， 这里及下文中的“f”表示此变量为仿真所获得的终端状态，“TAEM”表示此变量为期望的终 端状态。具体做法如下
在执行完一次弹道仿真后，由于和近似于0，假设和并 且假设然后利用公式(5)可得
其中，根据气动辨识技术测量到的当前气动拉偏信息所修正的终点升力 系数。
飞行器期望的终端状态应当满足
其中，表示第二次弹道仿真时应当使用的基准攻角，即第一次弹道仿真的基准攻角 修正之后的值。由于滑翔高度远远小于地球半径，本发明假设另外， 还假设然后将公式(85)、(86)两式相减，并经过代数变换，可得
线化可得
其中，且就是第一次弹道仿真后获得的 应当对基准攻角参数进行的修正量。是在考虑气动拉偏修正之后的升力线斜率。通 过求解公式(87–88)，可得
从而，修正之后的基准攻角参数为此值将作为下一次弹道仿真中的基准攻 角参数。
2)反转点的修正算法
这里利用反转点修正算法微调EBR2，以便满足终端速度要求。
首先定性分析EBR2与Vf的关系：若EBR2减小，最后一次反转推迟，进而在最后一次反转 时刻的航向角误差增大。因此，在高度调整阶段，比例导引律会产生更大的倾侧角来消除航 向角误差。这使得纵向升阻比减小，从而导致终点速度Vf减小。
基于上述单调关系，首先利用弹道仿真预测终点速度然后利用secant法修正EBR2， 如下
重复上述过程，直到
(2)高度调整阶段的基准倾侧角设计
此阶段，本发明采用比例导引进行基准倾侧角的设计。视线方位角变化率为
比例导引律产生的需用横向机动加速度aL2为
其中，ΔL2是由公式(22)计算。上式右边第二项是为了补偿地球自转的影响。为了防 止初始倾侧角饱和，这里令kPN从2逐渐变化到4，如下
另一方面，在平稳滑翔情况下，纵向平面内升力加速度aL1与重力、离心力和惯性力平衡， 即
其中ΔL1是由公式(21)计算。则高度调整阶段的基准倾侧角为
(3)生成高度调整阶段的攻角和倾侧角指令
高度调整阶段的攻角和倾侧角指令如下
其中，是剩余飞行距离—能量参考曲线，由步骤5(1)中的在线弹道仿真获得。 kα＝5π/(1.8×107)，其意味着当实际剩余飞行器距离偏离参考曲线1km时，对攻角修正0.05 度。在此阶段，各种干扰仍然会引起速度误差。由于调节攻角也可调节升阻比，为了保证终 端速度约束，这里通过调节攻角来跟踪由于此阶段的航程较短，干扰引起的速度误 差较小，上述跟踪策略并不会引起攻角的大幅调整。公式(97)右边第二项用来抑制弹道振荡。 为了保证热流密度等过程约束，这里亦采取如下修正措施
其中，由公式(81)计算确定。
本发明的优点在于：
(1)相较于传统的制导方法，该制导方法可以有效导引大升阻比飞行器的再入制导， 能够有效抑制高升阻比引起的弱阻尼长周期弹道振荡，降低热流密度的峰值，有利于对参考 轨迹的跟踪；
(2)由于解析规划参考剖面效率高，指令生成时间短，本制导方法具备快速响应紧急 任务的能力；
(3)该制导方法考虑地球自转的影响，将由地球自转引起的惯性力与实际气动力组合 为等效气动力，并利用解析解规划等效气动剖面，使飞行器具有全向再入任务处理能力；
(4)该制导方法提出了基于反比例函数的等效气动剖面拟合公式，并进一步推导了基 于复杂等效气动剖面形式的再入弹道解析解，具有高精度、全方位适应性和强鲁棒性。
附图说明
图1是再入轨迹示意图。
图2是纵程与横程的定义示意图。
图3是基准升阻比曲线。
图4是将牵连加速度和科氏加速度等效为附加气动力的示意图。
图5是剩余飞行距离示意图。
图6是传统制导方法的轨迹地面投影曲线
图7是传统制导方法的高度－速度曲线
图8是传统制导方法的倾侧角曲线
图9是本发明制导方法的轨迹地面投影曲线。
图10是本发明制导方法的高度－速度曲线。
图11是本发明制导方法的攻角随时间的变化曲线。
图12是本发明制导方法的倾侧角曲线。
图13是本发明制导方法的热流密度曲线。
图14是本发明制导方法的来流动压曲线。
图15是本发明制导方法的机动过载曲线。
具体实施方式
下面将结合附图和实施例对本发明做进一步的详细说明。
本发明是一种基于三维再入弹道解析解的全射向自主再入制导方法，此制导方法利用文 献二中的三维再入弹道解析解在线规划参考轨迹。此解析解仅考虑了地球曲率的影响，没有 考虑地球自转的影响。为了补偿地球自转的影响，这里以地面为参照物，将由地球自转引起 的惯性力和实际气动力组合为等效气动力，并将等效气动曲线作为参考剖面。为了合理耗散 飞行能量，防止指令出现剧烈变化，甚至饱和，本发明通过合理分析强非线性的惯性力剖面 变化规律，提取了其中近似为常值的项，并提出了基于反比例函数的等效气动剖面拟合公式。 然后，基于反比例函数形式气动剖面，本发明开发了更加复杂的新三维再入弹道解析解。进 而可以利用纵程解析解规划满足射程与终端能量要求的纵向参考剖面，并利用横程解析解规 划反转点，以消除横程误差。整个过程包括以下几个步骤：
步骤1：建立运动学与动力学方程
在研究再入制导问题时，常用经度λ、纬度φ和高度H来描述飞行器位置信息，而用速 率V、弹道倾角γ和航向角ψ来描述速度信息。考虑地球曲率和自转的影响，基于上述状态 变量的再入运动方程如下
其中，t为飞行时间，m是飞行器质量，为常值，kg。σ是倾侧角，rad。Re是地球平均半径，大小为6356.766km。L是升力，D是阻力，g为重力加速度，ωe为地球自转角速度。
一般，高超声速飞行器在飞行过程中需要满足热流密度来流动压q和机动过载约束n， 如下所示
其中，ρ是大气密度，kQ是与热流密度相关的常数，qmax是最大来流动压，是最大热流 密度，nmax是最大机动过载，g0是海平面高度的重力加速度。
另外，考虑到飞控系统能力的限制，对攻角指令和倾侧角指令变化率进行约束。攻角指 令约束：和倾侧角指令约束为：|σ|≤80deg、和
当飞行器抵达一个以目标为中心50km为半径的圆上时，再入段终止。此时需要满足的 终端条件是：期望的终点剩余飞行距离STAEM＝50km，期望的终端高度HTAEM＝25km，期望的终端速度VTAEM＝2000m/s，期望的终端航向误差|ΔψTAEM|≤5deg和期望的终端倾侧角 |σTAEM|≤30deg。这里及全文中下标“TAEM”代表期望的终端状态。
步骤2：本发明再入制导方法概述
图1为一条在所设计制导方法导引下的再入弹道轨迹，其中上子图是高度-纵程曲线，下 图是轨迹的地面投影曲线。根据弹道特点，再入过程一般被划分为三个阶段：下降段、平稳 滑翔阶段和高度调整阶段。在下降段，初始高度较高，大气密度稀薄，飞行器快速掉高。在 平稳滑翔阶段，由于大气密度可以提供足够的升力来平衡重力及离心力，滑翔高度平缓下降。 在此阶段，本制导方法利用三维再入弹道解析解实时解析规划等效气动参考剖面，并在数个 合适的时间点迭代调整倾侧反转点，以消除横程误差。在最后一次倾侧反转之后，此时飞行 器距离目标足够近，再入过程进入高度调整阶段。此时飞行器通过调节攻角曲线来达到期望 的飞行高度。下面步骤3-5详细介绍各个阶段所采用的制导策略。
步骤3：下降段制导策略
由于初始高度过高，升力过小，飞行器快速掉高。随着飞行器扎进稠密大气层，热流密 度会迅速增大，并在谷底达到局部峰值。为了保证热流密度约束，这里采用最保守的制导方 案，其采用最大可用攻角、零倾侧角下滑，从而尽可能地提升滑翔高度。当升力足以支撑飞 行器平稳滑翔时，攻角由最大值平滑地过渡到基准攻角，随后进入平稳滑翔阶段。
步骤4：平稳滑翔段制导策略
平稳滑翔阶段的制导策略流程如下：
1)设计一条基准攻角剖面，并确定与之对应的基准升阻比剖面。
2)设计纵向升阻比参数化剖面，其外形与基准升阻比剖面相似。这样有利于近似保持倾 侧角大小为常值。与此同时可以得到横向升阻比参数化剖面。气动剖面参数需要在步骤3) 与步骤4)中进行求解。
3)为了补偿地球自转的影响，将惯性力与实际气动力组合为等效气动力。通过合理评估 惯性力的影响，在上述常规气动剖面的基础上设计反比例函数形式的等效纵、横向升阻比曲 线。
4)推导基于反比例函数复杂形式等效气动剖面的新再入弹道解析解，然后根据纵程和终 端能量要求，利用纵程解析解求解出步骤2)中的气动剖面参数。
5)为了消除横程误差，利用横程解析确定倾侧反转点。为了减轻计算负担，本发明仅在 几个特殊时间节点修正反转点，而非实时修正。
6)调节基准倾侧角跟踪步骤2)中纵向升阻比曲线，并实现倾侧反转。
7)为了抑制弹道振荡，将弹道阻尼反馈引入到基准攻角和倾侧角中，从而得到制导指令。
8)将过程约束转化为倾侧角约束，进而限制倾侧角指令大小，以便保证过程约束。
9)从步骤2)开始重复上述过程，直到平稳滑翔阶段结束。
由于合理补偿地球自转的影响，上述策略可以保证攻角和倾侧角在此阶段近似为常值， 这有三个好处：(1)合理耗散飞行能量，避免制导指令在飞行后期饱和；(2)有利于飞行控 制系统精确跟踪指令；(3)由于气动拉偏系数对攻角敏感，常值攻角有利于准确评估气动拉 偏的影响，从而提高制导精度和鲁棒性。下面详细介绍制导策略。
图2中，CE是地心，M代表飞行器，T代表目的地。在由此三点构成的平面内构造一个以CE为圆心，地球平均半径Re为半径的大圆GC1。飞行器的速度矢量记为V。V与大圆GC1 的夹角记为Δψ，以大圆切线在当地水平面内沿逆时针方向转向V为正方向。因此，图中所 示的Δψ此时取负值。图中的曲线是一条再入飞行轨迹示意曲线，P1是轨迹上的任意一点。 过P1点作一个新的大圆GC2，其垂直于大圆GC1。这两个大圆的交点记为P2。那么定义P1点的纵程为大圆弧的长度，而P1点的横程为大圆弧的长度。
定义L1和L2分别是升力在纵平面和当地水平面内的分量，即L1＝Lcos(σ)和 L2＝Lsin(σ)。定义纵向升阻比为L1与阻力D的比值，即而横向升阻比 为L2与D的比值，即L2/D＝Lsin(σ)/D。定义能量E为
其中，μ是地球引力常数，为3.986005×1014m3/s2。
文献二以纵向升阻比和横向升阻比为控制量，开发了如下三维再入弹道解析 解，其可以预测纵程xD、横程xC和航向角误差Δψ。这些解析解考虑了地球曲率的影响，但 忽略了地球自转的影响。
其中，E0代表初始能量，E代表当前能量，xE为积分变量，定义域为(E,E0)。xC0为初始横 程，Δψ0为初始航向角误差。R*为常数，且有R*＝Re+H*，其中H*一般取飞行过程中高度 的中间值。f1(xE)为简记函数表达式，如下
从解析解中可以看出：在忽略地球自转的情况下，纵程与纵向升阻比几乎成正比，而横 程与横向升阻比几乎成正比。
(1)基准攻角和基准升阻比设计
以能量E为自变量，设计剖面函数。基准攻角αbsl的形式如下
其中，α1＝10deg近似对应最大升阻比攻角，以获得最远滑行距离，α2＝6deg为基准攻 角参数，Eα＝-5.55×107J/kg为临近平稳滑翔阶段和高度调整阶段的交班点，ETAEM为再入过 程中期望的终端能量。在有些扰动情况下，有可能出现E＜ETAEM的情况。此时，为了防止攻 角出现负值，设定限制：αbsl≥5deg。
图3为基准攻角对应的升阻比曲线，称为基准升阻比曲线，其中为基准升阻比。
(2)纵、横向升阻比设计
为了合理耗散能量，并减轻飞控系统负担，在平稳滑翔阶段倾侧角要尽可能保持常值。 按此要求，设计如下纵向升阻比曲线
其中，和是两个需要设计的待定升阻比参数，即气动剖面参数。由基准升阻 比曲线可以看出，在E≥Eα时令纵向升阻比保持常值可满足尽可能保持倾侧角大 小为常值的要求。在后续步骤中会介绍如何根据剩余飞行距离确定的大小；当E＜Eα 时，纵向升阻比要平滑地减小到但期望的终端倾侧角应近似为0，因此，
令其中
为补偿气拉偏影响的参数。是利用气动辨识技术测得的当前时刻实际升阻比，而 是当前状态对应的理想升阻比。则可得参数化的横向升阻比模值曲线为
(3)设计等效纵、横向升阻比并求解气动剖面参数
图4中，将惯性力等效为三项附加气动力ΔL1、ΔL2和ΔD。根据公式(102)、(103)、(104) 有。
由于γ≈0，假设sin(γ)＝0和cos(γ)＝1，则式(116－118)可被简化为
由于V一般是ωe(Re+H)cos(φ)的数倍、甚至数十倍，即2V>>ωe(Re+H)cos(φ)，则附加 气动力ΔL1和ΔL2的公式可以进一步被简化为
ΔL1＝2mVωecos(φ)sin(ψ) (122)
ΔL2＝2mVωesin(φ) (123)
定义等效纵向升阻比横向升阻比为
上式的一阶Taylor展开为
从上式中可以看到，如果高度越高，大气密度越小，D越小，则会对和产生很 大的影响。由于在大部分时间里飞行器平稳滑翔，所以D取平稳滑翔高度对应的阻力值，记 为DSG。则有
忽略地球自转的影响，并且在平稳滑翔状态下假设γ≈0和进而由公式(103)和 (108)得
其中，R*＝Re+H。
利用纵向升阻比可以计算出平稳滑翔高度对应的阻力值DSG，如下
其中，L1(SG)为平稳滑翔高度对应的升力值。下面分析等效纵向升阻比和等效横向升 阻比的变化趋势，并进行适当的拟合。将公式(121–123)和(131)代入到公式(128–129) 中，并整理可得
其中
其中
其中，h1,h2,h3,hz1,hz2,hz3,hm为拟合参数。
为了应用解析解，需要将h1–h3拟合成关于能量的表达式。由于分母已经是E的函数，这 里只需要拟合上述公式的分子。
通过观察仿真结果，发现cos(φ)sin(ψ)近似保持常值。在公式(121－123)中，V和sin(φ) 分别采用过两端点的线性函数来拟合，这表明Vsin(φ)适用利用二次多项式来拟合。
由于2V>>ωe(Re+H)cos(φ)，附加气动力ΔD(公式(121))要远远小于附加气动力ΔL1和 ΔL2。因此，尽管附加气动力ΔD变化无常，但是其影响较小，本发明粗略地利用线性公式来 近似附加气动力ΔD。因此，本发明利用线性公式拟合ΔL1和ΔD，而利用二次曲线拟合ΔL2。 根据上述分析，本发明采用线性公式拟合hz1和hz2，而利用二次多项式拟合hz3，如下
其中，xE为积分变量，物理意义为飞行器能量，上划线“-”表示此变量为拟合值，kh1(0)、kh1(1)、kh2(0)、kh2(1)、kh3(0)、kh3(1)、kh3(2)为相应拟合参数系数。和的参数满足这样的条 件：拟合线段经过原函数的两个端点，第一个端点由当前状态确定，第二个端点由期望的终 端状态确定。下面说明拟合参数系数的计算方法。
其中，
其中，是上一制导周期获得的当前纵向升阻比，是期望的终端 升阻比，φTAEM是目标的纬度，是飞行器的航向角ψ在某种意义下的加权平均值，是 飞行器的航向角ψ在终点处的加权平均值。
根据公式(123)和(137)，通过分别拟合V和sin(φ)，可以得到的表达式
展开上式右得拟合公式参数
接下来，利用纵程解析解确定的大小，以满足纵程和终端能量要求。根据公式(114) 和(34)，规划的等效纵向升阻比曲线可记为如下形式
将式(109)中的纵向升阻比用上式的等效升阻比替换，并积分可得基于反比 例函数气动剖面的纵程解析表达式。由于是分段函数，根据积分前后能量E1，E2与交班 点能量Eα的大小关系，下面需要考虑如下三种情况
1)当Eα≤E2≤E1时，将公式(150)中的第一个子式代入公式(109)，并积分可得纵程为
2)当E2≤E1≤Eα时，将公式(150)中的第二个子式代入公式(109)，并积分可得纵程为
其中，a0,a1,a2为简记符号，其值如下
3)当E2≤Eα≤E1时，则纵程则可以由上面两个公式组合计算，如下
xD(E2,E1)＝xD(Eα,E1)+xD(E2,Eα) (156)
为了计算待定升阻比参数需要首先计算剩余飞行距离sgo。图5为地心赤道旋转 坐标系CE-xeyeze(GER)，其原点在地心CE，xe和ye轴在赤道平面内，并且xe轴指向本初子午线，ze轴指向北极。M点代表飞行器，N点为M点向赤道平面的垂直投影点，T代表 目标。剩余飞行距离就是过M和T的大圆弧长度。定义为从CE到M的向量，定义为从CE到T的向量。剩余飞行距离由下式计算
sgo＝Reη (157)
其中η是和之间的夹角。下面求解η。
记xEM是的单位向量，xET是的单位向量。由于与赤道平面的夹角刚好就 是纬度φ，xEM在赤道平面内的分量大小为cos(φ)，而沿ze轴分量大小为sin(φ)。进而，利用 线段CEN与xe轴的夹角为经度λ，可以得到单位向量xEM沿xe轴的分量为cos(λ)cos(φ)，沿ye轴的分量为sin(λ)cos(φ)。因此，xEM在地心赤道旋转坐标系GER中的坐标为
其中，上标“GER”代表以GER系为坐标系，而上标“T”代表向量转置。同理可以得到xET在GER坐标系中的坐标
其中，λT和φT分别是目标点处的经度和纬度。从几何意义上，两个单位向量的点积刚好 等于它们之间夹角的余弦。因此，η可以由下式计算
接下来，利用剩余飞行距离sgo可以确定待定升阻比参数由于等效纵向升阻比曲 线是分段函数，下面需要分两种情况考虑：
1)当E＞Eα时，有如下等式关系
xD(ETAEM,Eα)+xD(Eα,E)＝sgo-STAEM (161)
其中，xD(ETAEM,Eα)由公式(152)计算，而xD(Eα,E)由公式(151)计算。通过求解公 式(161)，可以得到
其中，kxD1，kxD2为简记符号，其值如下
2)当E≤Eα时
此时飞行主要处于高度调整阶段，将采用其它制导策略。因此，无需再更新
(4)倾侧反转点规划
在飞行中，飞行器需要通过调节倾侧角来跟踪纵向升阻比曲线，但是这会导致大幅横向 机动。为了消除由此引起的横程误差，飞行器需要在适当的时候进行倾侧反转。本发明仅安 排两次反转，以减轻飞控系统负担。记倾侧反转对应的能量节点分别为EBR1和EBR2。在第一 次反转发生之前，令EBR2＝Eα，然后利用横程解析解来规划EBR1。在第一次反转发生之后， 采用一种基于在线弹道仿真的迭代修正算法微调EBR2。
在确定升阻比曲线之后，考虑到倾侧反转，利用公式(133)，等效横向升阻比曲线可以被表述为
其中，sgn是一个取值±1的符号函数，用于确定初始倾侧方向。这里sgn的取值使得飞 行器在初始时刻向航向误差减小的方向倾侧。kBR代表已经发生反转的次数，由公式 (17)计算。
记飞行器当前能量为E，此时飞行器的横程为xC，航向角误差为Δψ，将公式(165)代 入到横程解析解中(公式(110))中，并令飞行器能量从E变化到ETEAM，则可预测飞行器的终端横程xCf为
其中，xE为积分变量。xD(ETAEM,E)为飞行器能量从E积分到ETEAM时的纵程，xD(ETAEM,xE)为 飞行器能量从xE积分到ETEAM时的纵程，这两项可由公式(109)计算得到。F(EBR1,E)、 F(EBR2,EBR1)、F(ETAEM,EBR2)为函数F(xE2,xE1)分别在定义域(EBR1,E)、(EBR2,EBR1)、 (ETAEM,EBR2)上的积分值。
函数F(xE2,xE1)的定义为
其中简记函数表达式f2(xE)的形式如下
此外，由于规划的等效纵向升阻比曲线是分段函数，则纵程xD(ETAEM,xE)同样是分 段函数，如下所示
其中，xD(ETAEM,Eα),xD(ETAEM,xE)可由公式(152)计算，xD(Eα,xE)可由公式(151)计算。
由于期望终端横程xCf＝0，第一次倾侧反转对应的能量节点EBR1的数值可通过牛顿法迭 代求解方程xCf(EBR1)＝0得到，如下。
其中，上标“(k)”、“(k+1)”为能量节点EBR1的迭代次数，xCf对EBR1的导数x′Cf为
由于xCf随EBR1单调变化，一般只需要迭代数次就可以满足终端横程为 了减轻弹载计算机计算负担，本发明并不实时修正EBR1，而是仅在第一次反转发生之前的几 个适当时间节点进行修正。
(5)基准倾侧角设计
飞行器需要通过调节基准倾侧角来跟踪纵向升阻比曲线，并实现倾侧反转。基准倾 侧角σbsl如下
其中，是当前时刻利用惯导气动数据辨识的实际升阻比。ΔE是考虑飞行器的滚转 速率限幅引入的补偿量，如下
其中，是自动驾驶仪限制的最大倾侧角变化率，aD是当前的阻力加速度，可以利用 含加速度计的惯导组件测量。
(6)生成攻角指令和倾侧角指令
若直接采用公式(113)的基准攻角αbsl和公式(172)的基准倾侧角作为指令攻角αcmd和 σcmd，则滑翔弹道会存在弱阻尼、长周期振荡。因此，本发明采用弹道阻尼技术来抑制弹道 振荡。攻角指令αcmd和倾侧角指令σcmd为
αcmd＝αbsl+cos(σbsl)kγ(γSG-γ) (174)
其中，kγ是增益系数，α1是基准攻角剖面的一个参数，γSG为平稳滑翔弹道对应的弹道 倾角。γSG的高精度近似公式如下
其中，Dbsl＝CD(bsl)qSref是αbsl对应的阻力，CD(bsl)是αbsl对应的阻力系数，ρ是大气密度， q是动压，Sref是气动参考面积。系数d1和d2的计算公式如下
其中，CL(bsl)是基准攻角αbsl对应的升力系数。dCL(bsl)dE可提前计算好，d[cos(σbsl)]dE 可利用差分估算。
为了保证热流密度等过程约束，本发明将过程约束转化成倾侧角指令范围限制。对于平 稳滑翔弹道，大的倾侧角会导致高度降低加剧，从而大气密度增大，引起热流密度、来流动 压及过载的增大。最大的可用倾侧角发生在滑翔高度下限Hmin处，此时大气密度达到最大， 升力也达到最大，升力在垂直方向上的分量近似满足平稳滑翔条件。因此平稳滑翔状态下的 最大可用倾侧角可由公式(179)右边第一项计算。当有弹道振荡时，如果飞行器快速掉高， 需要减小σmax，以防止飞行器掉到Hmin以下。此时，添加反馈补偿，即公式(179)右边第二 项。
其中，L1是升力在纵平面内的分量，Lmax为可达到的最大升力，参数kσ＝-50，且有
其中，ΔL1由公式(119)计算，是利用气动辨识技术测量得到的当前实际升力系数。 进而采取如下措施
步骤5：高度调整段制导策略
在最后一次反转发生之后，飞行器进入高度调整阶段。在此阶段，为了达到期望的终端 高度，飞行器的基准攻角逐渐被调整至α2(在公式(113)中定义)。在此阶段，由于弹道倾 角变化率较大，再入滑翔解析解精度不能满足终端要求。为此，这里本发明采用基于比例导 引律的制导策略，以消除终端航向误差。但是，不再跟踪纵向升阻比曲线会造成终端速度误 差。经理论分析可知，终点速度随最后一次反转点单调变化。为了满足终端速度要求，本发 明采取这样一个方案：在最后一次反转发生之前，利用基于弹道仿真的修正算法对第二次反 转点进行修正，以消除反转之前积累的速度误差；在最后一次反转之后，通过调节攻角来跟 踪剩余飞行距离—能量参考曲线，以克服干扰。由于这里仅对剩余一小段弹道进行数次仿真， 而非是实时仿真，所以并不会造成较大的计算负担。另外，修正算法还会对基准攻角曲线进 行略微调整，以提高终端高度精度。
(1)基准攻角和最后一次反转点修正
在进入高度调整阶段之前，需要利用基于弹道仿真的迭代修正算法，对基准攻角参数α2和第二次反转点的能量EBR2进行微调，以满足终端速度和高度约束。为了减轻计算负担，本 发明在第二次反转发生之前仅进行两轮修正：第一轮发生在第一次反转之后不久；第二轮大 约发生在第二次反转发生之前的60s。
首先介绍在线弹道仿真。在第n次弹道仿真中，首先根据已测的气动拉偏信息对气动模 型粗略修正，然后根据上一步的迭代结果更新第二次反转点能量和基准攻角参数，即和 最后以本发明再入制导方法得到的控制律，对运动方程公式(99–104)进行数值积分， 直到为止。这里及下文，上标“(n)”表示此变量为第n次弹道仿真的值，n＝1,2,3,…。
迭代修正算法包含两个子算法：基准攻角修正算法和反转点修正算法。首先，调用基准 攻角修正算法修正基准攻角曲线，以便提高终端高度精度。此任务仅需要执行一次弹道仿真。 接着，本发明利用反转点修正算法精确微调最后一次反转点EBR2，以便满足终端速度约束。 这两种算法的详细介绍如下。
1)基准攻角修正算法
此算法通过执行一次弹道仿真，预测包括终端高度在内的终端状态，然后，通过比 较和HTAEM来修正基准攻角参数α2，以提高终端高度精度。注意区分下标“f”和“TAEM”， 这里及下文中的“f”表示此变量为仿真所获得的终端状态，“TAEM”表示此变量为期望的终 端状态。具体做法如下
在执行完一次弹道仿真后，由于和近似于0，假设和并 且假设然后利用公式(103)可得
其中，根据气动辨识技术测量到的当前气动拉偏信息所修正的终点升力 系数。
飞行器期望的终端状态应当满足
其中，表示第二次弹道仿真时应当使用的基准攻角，即第一次弹道仿真的基准攻角 修正之后的值。由于滑翔高度远远小于地球半径，本发明假设另外， 还假设然后将公式(183)、(184)两式相减，并经过代数变换，可得
线化可得
其中，且就是第一次弹道仿真后获得的 应当对基准攻角参数进行的修正量。是在考虑气动拉偏修正之后的升力线斜率。通 过求解公式(185–186)，可得
从而，修正之后的基准攻角参数为此值将作为下一次弹道仿真中的基准攻 角参数。
2)反转点的修正算法
这里利用反转点修正算法微调EBR2，以便满足终端速度要求。
首先定性分析EBR2与Vf的关系：若EBR2减小，最后一次反转推迟，进而在最后一次反转 时刻的航向角误差增大。因此，在高度调整阶段，比例导引律会产生更大的倾侧角来消除航 向角误差。这使得纵向升阻比减小，从而导致终点速度Vf减小。
基于上述单调关系，首先利用弹道仿真预测终点速度然后利用secant法修正EBR2， 如下
重复上述过程，直到
(2)高度调整阶段的基准倾侧角设计
此阶段，本发明采用比例导引进行基准倾侧角的设计。视线方位角变化率为
比例导引律产生的需用横向机动加速度aL2为
其中，ΔL2是由公式(120)计算。上式右边第二项是为了补偿地球自转的影响。为了防 止初始倾侧角饱和，这里令kPN从2逐渐变化到4，如下
另一方面，在平稳滑翔情况下，纵向平面内升力加速度aL1与重力、离心力和惯性力平衡， 即
其中ΔL1是由公式(119)计算。则高度调整阶段的基准倾侧角为
(3)生成高度调整阶段的攻角和倾侧角指令
高度调整阶段的攻角和倾侧角指令如下
其中，是剩余飞行距离—能量参考曲线，由步骤5(1)中的在线弹道仿真获得。 kα＝5π/(1.8×107)，其意味着当实际剩余飞行器距离偏离参考曲线1km时，对攻角修正0.05 度。在此阶段，各种干扰仍然会引起速度误差。由于调节攻角也可调节升阻比，为了保证终 端速度约束，这里通过调节攻角来跟踪由于此阶段的航程较短，干扰引起的速度误 差较小，上述跟踪策略并不会引起攻角的大幅调整。公式(195)右边第二项用来抑制弹道振荡。 为了保证热流密度等过程约束，这里亦采取如下修正措施
其中，由公式(179)计算确定。
实施例
传统的再入制导方法没有充分考虑地球自转的影响，不具备全射向再入制导能力，尤其 是在执行南北方向飞行的任务时，存在显著的跟踪误差，且终端约束误差也比本发明制导律 高出一个数量级。在本实施例中，提供七不同射向的例子，验证本发明全射向自主再入制导 方法的高精度、全方位适应性和强鲁棒性，同时给出文献二中再入制导方法的仿真结果作为 对比。设定初始条件为H0＝80km,λ0＝00,φ0＝500,V0＝7000m/s,γ0＝00。根据目标位置的差 异，表1列出了不同射向下的初始航向角ψ0及终端经度、纬度要求。当飞行器抵达一个以目 标为中心50km为半径的圆上时，再入段终止。此时需要满足的终端条件是：STAEM＝50km， HTAEM＝25km，VTAEM＝2000m/s，|ΔψTAEM|≤5deg和|σTAEM|≤30deg。
表1
从仿真结果表2可以看出，在传统再入制导方法导引下，飞行器并没有成功处理T3和 T4的任务，终端速度、高度以及航向角远远偏离期望值，造成很大的终端误差。这主要是由 于未充分考虑地球自转的影响，因此在执行T3和T4这两个接近南北方向飞行的任务时，跟 踪误差较大，同时由于在飞行中未能修正基准攻角参数α2，造成终端误差较大。图6是传统 再入制导方法的再入轨迹在地面投影曲线，图7是对应的高度－速度曲线，可以看出，T4和 T5情况下的热流密度在一段时间内超出了过程约束限制。图8是倾侧角曲线，可以看出，由 于T3和T4情况任务失败，传统再入制导方法的倾侧角曲线在最后的高度调整阶段也出现异 常现象。
表2
从表3可以看出，由于考虑了地球自转的影响，将由地球自转引起的惯性力与实际气动 力组合为等效气动力，并利用解析解规划等效气动剖面，并进一步推导了基于复杂等效气动 剖面形式的再入弹道解析解，本发明制导方法成功应对了沿不同方向的再入任务，并获得了 很高的制导精度，特别是由于在高度调整阶段采用了比例导引律，终端航向误差非常小。从 表3中可以看出，为了提高终端高度精度，基于弹道仿真的迭代修正算法对基准攻角曲线参 数α2进行了微调。另外，EBR2分布比较集中，近似等于Eα。可见迭代修正算法仅对EBR2进行 了微调，这间接验证了解析解的高精度和参考轨迹规划方法的合理性。图9展示了所有轨迹 的地面投影曲线，并给出了T1情况在终点处的局部放大图。从放大图可以看出，所设计的制 导方法成功将飞行器导引到距离目标50km的圆上。其它情况的终止条件与T1情况相同。图10是高度-速度曲线，从中可见，在平稳滑翔阶段和高度调整阶段高度下降很平稳，弱阻尼长 周期振荡得到有效抑制。图11是攻角曲线，从中可以清晰地分辨出飞行的三个阶段。图12 是倾侧角曲线。从中可以看出，倾侧角在平稳滑翔阶段近似为常值，而在高度调整阶段逐渐 收敛到零附近，相对传统制导方法也较平滑。图13是热流密度曲线，从中可见，情况T4和 T5的热流峰值也同样有一段时间超出了限制要求。这两种情况的共同点是飞行器均向西飞 行，而且纬度较低。因此，公式(103)中惯性力的主导项2Vωecos(φ)sin(ψ)为负值，此时惯 性力方向向下，从而压低飞行高度，造成向西飞行将面临更加严峻的热流考验。图14–15分 别是来流动压和机动过载曲线。从中可见，由于动压直接影响过载，动压和过载曲线比较相 似。另外，在平稳滑翔阶段动压和过载近似线性单调增长，这主要是由于速度减小导致地球 曲率引起的离心力减小，从而需要增大升力以平衡重力。
表3。

Claims (1)

1.一种基于三维再入弹道解析解的全射向自主再入制导方法，其特征在于：
步骤1：建立运动学与动力学方程
在研究再入制导问题时，用经度λ、纬度φ和高度H来描述飞行器位置信息，而用速率V、弹道倾角γ和航向角ψ来描述速度信息；考虑地球曲率和自转的影响，基于上述状态变量的再入运动方程如下：
其中，t为飞行时间，m是飞行器质量，为常值，kg；σ是倾侧角，rad；Re是地球平均半径，大小为6356.766km；L是升力，D是阻力，g为重力加速度，ωe为地球自转角速度；
高超声速飞行器在飞行过程中需要满足热流密度来流动压q和机动过载约束n，如下所示：
其中，ρ是大气密度，kQ是与热流密度相关的常数，qmax是最大来流动压，是最大热流密度，nmax是最大机动过载，g0是海平面高度的重力加速度；
另外，考虑到飞控系统能力的限制，对攻角指令和倾侧角指令变化率进行约束；攻角指令约束：和倾侧角指令约束为：|σ|≤80deg、和
当飞行器抵达一个以目标为中心50km为半径的圆上时，再入段终止；此时需要满足的终端条件是：期望的终点剩余飞行距离STAEM＝50km，期望的终端高度HTAEM＝25km，期望的终端速度VTAEM＝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)开始重复上述过程，直到平稳滑翔阶段结束；
假设CE是地心，M代表飞行器，T代表目的地；在由此三点构成的平面内构造一个以CE为圆心，地球平均半径Re为半径的大圆GC1；飞行器的速度矢量记为V；V与大圆GC1的夹角记为Δψ，以大圆切线在当地水平面内沿逆时针方向转向V为正方向；因此，Δψ此时取负值；假设P1是轨迹上的任意一点；过P1点作一个新的大圆GC2，其垂直于大圆GC1；这两个大圆的交点记为P2；那么定义P1点的纵程为大圆弧的长度，而P1点的横程为大圆弧的长度；
定义L1和L2分别是升力在纵平面和当地水平面内的分量，即L1＝Lcos(σ)和L2＝Lsin(σ)；定义纵向升阻比为L1与阻力D的比值，即而横向升阻比为L2与D的比值，即L2/D＝Lsin(σ)/D；定义能量E为：
其中，μ是地球引力常数，为3.986005×1014m3/s2；
背景技术中文献二以纵向升阻比和横向升阻比为控制量，开发了如下三维再入弹道解析解，其预测纵程xD、横程xC和航向角误差Δψ；这些解析解考虑了地球曲率的影响，但忽略了地球自转的影响；
其中，E0代表初始能量，E代表当前能量，xE为积分变量，定义域为(E,E0)；xC0为初始横程，Δψ0为初始航向角误差；R*为常数，且有R*＝Re+H*，其中H*取飞行过程中高度的中间值；f1(xE)为简记函数表达式，如下：
从解析解中看出：在忽略地球自转的情况下，纵程与纵向升阻比几乎成正比，而横程与横向升阻比几乎成正比；
(1)基准攻角和基准升阻比设计
以能量E为自变量，设计剖面函数；基准攻角αbsl的形式如下：
其中，α1＝10deg近似对应最大升阻比攻角，以获得最远滑行距离，α2＝6deg为基准攻角参数，Eα＝-5.55×107J/kg为临近平稳滑翔阶段和高度调整阶段的交班点，ETAEM为再入过程中期望的终端能量；在有些扰动情况下，出现E＜ETAEM的情况；此时，为了防止攻角出现负值，设定限制：αbsl≥5deg；则画出基准攻角对应的升阻比曲线，称为基准升阻比曲线；
(2)纵、横向升阻比设计
为了合理耗散能量，并减轻飞控系统负担，在平稳滑翔阶段倾侧角要保持常值；按此要求，设计如下纵向升阻比曲线：
其中，和是两个需要设计的待定升阻比参数，即气动剖面参数；由基准升阻比曲线看出，在E≥Eα时令纵向升阻比保持常值满足保持倾侧角大小为常值的要求；在后续步骤中会介绍如何根据剩余飞行距离确定的大小；当E＜Eα时，纵向升阻比要平滑地减小到但期望的终端倾侧角应近似为0，因此，令其中为补偿气拉偏影响的参数；是利用气动辨识技术测得的当前时刻实际升阻比，而是当前状态对应的理想升阻比；则得到参数化的横向升阻比模值曲线为
(3)设计等效纵、横向升阻比并求解气动剖面参数
将惯性力等效为三项附加气动力ΔL1、ΔL2和ΔD；根据公式(4)、(5)、(6)有：
由于γ≈0，假设sin(γ)＝0和cos(γ)＝1，则式(18－20)被简化为：
由于V是ωe(Re+H)cos(φ)的数倍、甚至数十倍，即2V＞＞ωe(Re+H)cos(φ)，则附加气动力ΔL1和ΔL2的公式进一步被简化为：
ΔL1＝2mVωecos(φ)sin(ψ) (24)
ΔL2＝2mVωesin(φ) (25)
定义等效纵向升阻比横向升阻比为
上式的一阶Taylor展开为：
从上式中看到，如果高度越高，大气密度越小，D越小，则会对和产生很大的影响；由于在大部分时间里飞行器平稳滑翔，所以D取平稳滑翔高度对应的阻力值，记为DSG；则有
忽略地球自转的影响，并且在平稳滑翔状态下假设γ≈0和进而由公式(5)和(10)得：
其中，R*＝Re+H；
利用纵向升阻比计算出平稳滑翔高度对应的阻力值DSG，如下：
其中，L1(SG)为平稳滑翔高度对应的升力值；下面分析等效纵向升阻比和等效横向升阻比的变化趋势，并进行适当的拟合；将公式(23–25)和(33)代入到公式(30–31)中，并整理得到：
其中，
其中，
其中，h1,h2,h3,hz1,hz2,hz3,hm为拟合参数；
为了应用解析解，需要将h1–h3拟合成关于能量的表达式；由于分母已经是E的函数，这里只需要拟合上述公式的分子；
通过观察仿真结果，发现cos(φ)sin(ψ)近似保持常值；在公式(23－25)中，V和sin(φ)分别采用过两端点的线性函数来拟合，这表明Vsin(φ)适用利用二次多项式来拟合；由于2V＞＞ωe(Re+H)cos(φ)，附加气动力ΔD要远远小于附加气动力ΔL1和ΔL2；因此，尽管附加气动力ΔD变化无常，但是其影响较小，利用线性公式来近似附加气动力ΔD；因此，利用线性公式拟合ΔL1和ΔD，而利用二次曲线拟合ΔL2；根据上述分析，采用线性公式拟合hz1和hz2，而利用二次多项式拟合hz3，如下：
其中，xE为积分变量，物理意义为飞行器能量，上划线“-”表示此变量为拟合值，kh1(0)、kh1(1)、kh2(0)、kh2(1)、kh3(0)、kh3(1)、kh3(2)为相应拟合参数系数；和的参数满足这样的条件：拟合线段经过原函数的两个端点，第一个端点由当前状态确定，第二个端点由期望的终端状态确定；下面说明拟合参数系数的计算方法；
其中，
其中，是上一制导周期获得的当前纵向升阻比，是期望的终端升阻比，φTAEM是目标的纬度，是飞行器的航向角ψ在某种意义下的加权平均值，是飞行器的航向角ψ在终点处的加权平均值；
根据公式(25)和(39)，通过分别拟合V和sin(φ)，得到的表达式：
展开上式右得拟合公式参数：
接下来，利用纵程解析解确定的大小，以满足纵程和终端能量要求；根据公式(16)和(34)，规划的等效纵向升阻比曲线记为如下形式：
将式(11)中的纵向升阻比用上式的等效升阻比替换，并积分得基于反比例函数气动剖面的纵程解析表达式；由于是分段函数，根据积分前后能量E1，E2与交班点能量Eα的大小关系，下面需要考虑如下三种情况：
1)当Eα≤E2≤E1时，将公式(52)中的第一个子式代入公式(11)，并积分得到纵程为：
2)当E2≤E1≤Eα时，将公式(52)中的第二个子式代入公式(11)，并积分得到纵程为：
其中，a0,a1,a2为简记符号，其值如下：
3)当E2≤Eα≤E1时，则纵程则由上面两个公式组合计算，如下：
xD(E2,E1)＝xD(Eα,E1)+xD(E2,Eα) (58)
为了计算待定升阻比参数需要首先计算剩余飞行距离sgo；建立地心赤道旋转坐标系GER CE-xeyeze，其原点在地心CE，xe和ye轴在赤道平面内，并且xe轴指向本初子午线，ze轴指向北极；M点代表飞行器，N点为M点向赤道平面的垂直投影点，T代表目标；剩余飞行距离就是过M和T的大圆弧长度；定义为从CE到M的向量，定义为从CE到T的向量；剩余飞行距离由下式计算：
sgo＝Reη (59)
其中η是和之间的夹角；下面求解η；
记xEM是的单位向量，xET是的单位向量；由于与赤道平面的夹角刚好就是纬度φ，xEM在赤道平面内的分量大小为cos(φ)，而沿ze轴分量大小为sin(φ)；进而，利用线段CEN与xe轴的夹角为经度λ，得到单位向量xEM沿xe轴的分量为cos(λ)cos(φ)，沿ye轴的分量为sin(λ)cos(φ)；因此，xEM在地心赤道旋转坐标系GER中的坐标为
其中，上标“GER”代表以GER系为坐标系，而上标“T”代表向量转置；同理得到xET在GER坐标系中的坐标
其中，λT和φT分别是目标点处的经度和纬度；两个单位向量的点积刚好等于它们之间夹角的余弦；因此，η由下式计算：
接下来，利用剩余飞行距离sgo确定待定升阻比参数由于等效纵向升阻比曲线是分段函数，下面需要分两种情况考虑：
1)当E＞Eα时，有如下等式关系：
xD(ETAEM,Eα)+xD(Eα,E)＝sgo-STAEM (63)
其中，xD(ETAEM,Eα)由公式(54)计算，而xD(Eα,E)由公式(53)计算；通过求解公式(63)，得到
其中，kxD1，kxD2为简记符号，其值如下：
2)当E≤Eα时
此时飞行处于高度调整阶段，将采用其它制导策略；因此，无需再更新
(4)倾侧反转点规划
在飞行中，飞行器需要通过调节倾侧角来跟踪纵向升阻比曲线，但是这会导致大幅横向机动；为了消除由此引起的横程误差，飞行器需要进行倾侧反转；该制导方法仅安排两次反转，以减轻飞控系统负担；记倾侧反转对应的能量节点分别为EBR1和EBR2；在第一次反转发生之前，令EBR2＝Eα，然后利用横程解析解来规划EBR1；在第一次反转发生之后，采用一种基于在线弹道仿真的迭代修正算法微调EBR2；
在确定升阻比曲线之后，考虑到倾侧反转，利用公式(35)，等效横向升阻比曲线被表述为：
其中，sgn是一个取值±1的符号函数，用于确定初始倾侧方向；这里sgn的取值使得飞行器在初始时刻向航向误差减小的方向倾侧；kBR代表已经发生反转的次数，由公式(17)计算；
记飞行器当前能量为E，此时飞行器的横程为xC，航向角误差为Δψ，将公式(67)代入到横程解析解中即公式(12)，并令飞行器能量从E变化到ETEAM，则预测飞行器的终端横程xCf为：
其中，xE为积分变量；xD(ETAEM,E)为飞行器能量从E积分到ETEAM时的纵程，xD(ETAEM,xE)为飞行器能量从xE积分到ETEAM时的纵程，这两项由公式(11)计算得到；F(EBR1,E)、F(EBR2,EBR1)、F(ETAEM,EBR2)为函数F(xE2,xE1)分别在定义域(EBR1,E)、(EBR2,EBR1)、(ETAEM,EBR2)上的积分值；
函数F(xE2,xE1)的定义为：
其中简记函数表达式f2(xE)的形式如下：
此外，由于规划的等效纵向升阻比曲线是分段函数，则纵程xD(ETAEM,xE)同样是分段函数，如下所示：
其中，xD(ETAEM,Eα),xD(ETAEM,xE)由公式(54)计算，xD(Eα,xE)由公式(53)计算；
由于期望终端横程xCf＝0，第一次倾侧反转对应的能量节点EBR1的数值通过牛顿法迭代求解方程xCf(EBR1)＝0得到，如下；
其中，上标“(k)”、“(k+1)”为能量节点EBR1的迭代次数，xCf对EBR1的导数x′Cf为
由于xCf随EBR1单调变化，只需要迭代数次就满足终端横程为了减轻弹载计算机计算负担，并不实时修正EBR1，而是仅在第一次反转发生之前的几个适当时间节点进行修正；
(5)基准倾侧角设计
飞行器需要通过调节基准倾侧角来跟踪纵向升阻比曲线，并实现倾侧反转；基准倾侧角σbsl如下：
其中，是当前时刻利用惯导气动数据辨识的实际升阻比；ΔE是考虑飞行器的滚转速率限幅引入的补偿量，如下：
其中，是自动驾驶仪限制的最大倾侧角变化率，aD是当前的阻力加速度，利用含加速度计的惯导组件测量；
(6)生成攻角指令和倾侧角指令
若直接采用公式(15)的基准攻角αbsl和公式(74)的基准倾侧角作为指令攻角αcmd和σcmd，则滑翔弹道会存在弱阻尼、长周期振荡；因此，采用弹道阻尼技术来抑制弹道振荡；攻角指令αcmd和倾侧角指令σcmd为：
αcmd＝αbsl+cos(σbsl)kγ(γSG-γ) (76)
其中，kγ是增益系数，α1是基准攻角剖面的一个参数，γSG为平稳滑翔弹道对应的弹道倾角；γSG的高精度近似公式如下：
其中，Dbsl＝CD(bsl)qSref是αbsl对应的阻力，CD(bsl)是αbsl对应的阻力系数，ρ是大气密度，q是动压，Sref是气动参考面积；系数d1和d2的计算公式如下：
其中，CL(bsl)是基准攻角αbsl对应的升力系数；dCL(bsl)/dE提前计算好，d[cos(σbsl)]/dE利用差分估算；
为了保证热流密度等过程约束，将过程约束转化成倾侧角指令范围限制；对于平稳滑翔弹道，大的倾侧角会导致高度降低加剧，从而大气密度增大，引起热流密度、来流动压及过载的增大；最大的用倾侧角发生在滑翔高度下限Hmin处，此时大气密度达到最大，升力也达到最大，升力在垂直方向上的分量近似满足平稳滑翔条件；因此平稳滑翔状态下的最大用倾侧角由公式(81)右边第一项计算；当有弹道振荡时，如果飞行器快速掉高，需要减小σmax，以防止飞行器掉到Hmin以下；此时，添加反馈补偿，即公式(81)右边第二项；
其中，L1是升力在纵平面内的分量，Lmax为达到的最大升力，参数kσ＝-50，且有
其中，ΔL1由公式(21)计算，是利用气动辨识技术测量得到的当前实际升力系数；公式如下：
步骤5：高度调整段制导策略
(1)基准攻角和最后一次反转点修正
在进入高度调整阶段之前，需要利用基于弹道仿真的迭代修正算法，对基准攻角参数α2和第二次反转点的能量EBR2进行微调，以满足终端速度和高度约束；为了减轻计算负担，在第二次反转发生之前仅进行两轮修正：第一轮发生在第一次反转之后不久；第二轮大约发生在第二次反转发生之前的60s；
首先介绍在线弹道仿真；在第n次弹道仿真中，首先根据已测的气动拉偏信息对气动模型粗略修正，然后根据上一步的迭代结果更新第二次反转点能量和基准攻角参数，即和最后以再入制导方法得到的控制律，对运动方程公式(1–6)进行数值积分，直到为止；上标“(n)”表示此变量为第n次弹道仿真的值，n＝1,2,3,…；
迭代修正算法包含两个子算法：基准攻角修正算法和反转点修正算法；首先，调用基准攻角修正算法修正基准攻角曲线，以便提高终端高度精度；此任务仅需要执行一次弹道仿真；接着，利用反转点修正算法精确微调最后一次反转点EBR2，以便满足终端速度约束；这两种算法的详细介绍如下；
1)基准攻角修正算法
通过执行一次弹道仿真，预测包括终端高度在内的终端状态，然后，通过比较和HTAEM来修正基准攻角参数α2，以提高终端高度精度；注意区分下标“f”和“TAEM”，“f”表示此变量为仿真所获得的终端状态，“TAEM”表示此变量为期望的终端状态；具体做法如下：
在执行完一次弹道仿真后，由于和近似于0，假设和并且假设然后利用公式(5)得到
其中，根据气动辨识技术测量到的当前气动拉偏信息所修正的终点升力系数；
飞行器期望的终端状态应当满足
其中，表示第二次弹道仿真时应当使用的基准攻角，即第一次弹道仿真的基准攻角修正之后的值；由于滑翔高度远远小于地球半径，假设另外，还假设然后将公式(85)、(86)两式相减，并经过代数变换，得到：
线化得到：
其中，且就是第一次弹道仿真后获得的应当对基准攻角参数进行的修正量；是在考虑气动拉偏修正之后的升力线斜率；通过求解公式(87–88)，得到：
从而，修正之后的基准攻角参数为此值将作为下一次弹道仿真中的基准攻角参数；
2)反转点的修正算法
利用反转点修正算法微调EBR2，以便满足终端速度要求；
首先定性分析EBR2与Vf的关系：若EBR2减小，最后一次反转推迟，进而在最后一次反转时刻的航向角误差增大；因此，在高度调整阶段，比例导引律会产生更大的倾侧角来消除航向角误差；这使得纵向升阻比减小，从而导致终点速度Vf减小；
基于上述单调关系，首先利用弹道仿真预测终点速度然后利用secant法修正EBR2，如下：
重复上述过程，直到
(2)高度调整阶段的基准倾侧角设计
此阶段，采用比例导引进行基准倾侧角的设计；视线方位角变化率为：
比例导引律产生的需用横向机动加速度aL2为：
其中，ΔL2是由公式(22)计算；上式右边第二项是为了补偿地球自转的影响；为了防止初始倾侧角饱和，这里令kPN从2逐渐变化到4，如下
另一方面，在平稳滑翔情况下，纵向平面内升力加速度aL1与重力、离心力和惯性力平衡，即：
其中ΔL1是由公式(21)计算；则高度调整阶段的基准倾侧角为：
(3)生成高度调整阶段的攻角和倾侧角指令
高度调整阶段的攻角和倾侧角指令如下：
其中，是剩余飞行距离—能量参考曲线，由步骤5(1)中的在线弹道仿真获得；kα＝5π/(1.8×107)，其意味着当实际剩余飞行器距离偏离参考曲线1km时，对攻角修正0.05度；在此阶段，各种干扰仍然会引起速度误差；由于调节攻角也能调节升阻比，为了保证终端速度约束，通过调节攻角来跟踪由于此阶段的航程短，干扰引起的速度误差小，上述跟踪策略并不会引起攻角的大幅调整；公式(97)右边第二项用来抑制弹道振荡；为了保证热流密度等过程约束，这里亦采取如下修正措施：
其中，由公式(81)计算确定。
CN201711260656.7A 2017-12-04 2017-12-04 一种基于三维再入弹道解析解的全射向自主再入制导方法 Active CN108036676B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711260656.7A CN108036676B (zh) 2017-12-04 2017-12-04 一种基于三维再入弹道解析解的全射向自主再入制导方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711260656.7A CN108036676B (zh) 2017-12-04 2017-12-04 一种基于三维再入弹道解析解的全射向自主再入制导方法

Publications (2)

Publication Number Publication Date
CN108036676A true CN108036676A (zh) 2018-05-15
CN108036676B CN108036676B (zh) 2019-08-23

Family

ID=62095187

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711260656.7A Active CN108036676B (zh) 2017-12-04 2017-12-04 一种基于三维再入弹道解析解的全射向自主再入制导方法

Country Status (1)

Country Link
CN (1) CN108036676B (zh)

Cited By (23)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108710383A (zh) * 2018-05-25 2018-10-26 哈尔滨工业大学 一种基于航点规划与跟踪的四旋翼飞行器自主降落控制方法
CN109002576A (zh) * 2018-06-11 2018-12-14 北京航空航天大学 一种线性高阶比例制导系统脱靶量的幂级数解
CN109062241A (zh) * 2018-09-29 2018-12-21 北京航空航天大学 基于线性伪谱模型预测控制的自主全射向再入制导方法
CN109190248A (zh) * 2018-09-03 2019-01-11 中国运载火箭技术研究院 一种用于滑翔飞行器的滑翔射程解析方法及解析系统
CN109270960A (zh) * 2018-12-05 2019-01-25 中南大学 基于Radau伪谱法的在线最优反馈再入制导方法
CN109446582A (zh) * 2018-09-29 2019-03-08 北京航空航天大学 一种考虑地球自转的高精度降阶平稳滑翔动力学建模方法
CN109708525A (zh) * 2018-12-12 2019-05-03 中国人民解放军陆军工程大学 一种导弹飞行弹道的解算方法、系统及终端设备
CN109740198A (zh) * 2018-12-14 2019-05-10 中国人民解放军国防科技大学 一种基于解析预测的滑翔飞行器三维再入制导方法
CN110147521A (zh) * 2019-04-25 2019-08-20 北京航空航天大学 一种高超声速飞行器跳跃滑翔弹道解析求解方法
CN110733670A (zh) * 2019-11-05 2020-01-31 中国人民解放军国防科技大学 一种短航程低过载的再入轨迹设计方法
CN111306989A (zh) * 2020-03-12 2020-06-19 北京航空航天大学 一种基于平稳滑翔弹道解析解的高超声速再入制导方法
CN111351401A (zh) * 2018-12-21 2020-06-30 北京理工大学 应用于捷联导引头制导飞行器的防侧偏制导方法
CN111366044A (zh) * 2019-12-29 2020-07-03 湖北航天飞行器研究所 一种平飞过渡段制导控制方法
CN111397441A (zh) * 2019-01-03 2020-07-10 北京理工大学 带有捷联导引头的远程制导飞行器的全射程覆盖制导系统
CN109508030B (zh) * 2018-11-27 2020-08-04 北京航空航天大学 一种考虑多禁飞区约束的协同解析再入制导方法
CN111580555A (zh) * 2020-05-13 2020-08-25 北京控制工程研究所 一种高超声速飞行器上升段分段自适应预测校正制导方法
CN112198894A (zh) * 2020-07-31 2021-01-08 北京理工大学 旋翼无人机自主移动降落制导方法及系统
CN112507466A (zh) * 2020-12-21 2021-03-16 中国人民解放军96901部队23分队 考虑地球自转速率的滑翔弹道随能量变化降阶解
CN112541305A (zh) * 2020-12-10 2021-03-23 中国航空工业集团公司沈阳飞机设计研究所 一种基于全局变量求导的飞机气动力特性分析方法
CN113008080A (zh) * 2021-01-26 2021-06-22 河北汉光重工有限责任公司 一种基于刚性原理对近海目标进行火控解算方法
US20210216688A1 (en) * 2020-01-13 2021-07-15 Microsoft Technology Licensing, Llc Configuring aerodynamic simulation of a virtual object
CN114740894A (zh) * 2022-05-13 2022-07-12 北京航空航天大学 基于注意力机制与门控循环单元的飞行器制导方法和系统
CN115542746A (zh) * 2022-12-05 2022-12-30 北京航空航天大学 高超声速飞行器的能量管控再入制导方法及装置

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20120045258A (ko) * 2010-10-29 2012-05-09 삼성중공업 주식회사 선박용 프로펠러 설계방법, 이에 의하여 설계된 선박용 프로펠러, 및 이에 의하여 설계된 선박용 프로펠러를 갖는 선박
CN102880187A (zh) * 2012-09-21 2013-01-16 北京控制工程研究所 一种跳跃式再入飞行器初次再入段横向制导方法
CN103863579A (zh) * 2014-03-31 2014-06-18 北京控制工程研究所 一种深空探测返回过程的预测校正制导方法
CN104035335A (zh) * 2014-05-27 2014-09-10 北京航空航天大学 基于高精度纵、横程解析预测方法的平稳滑翔再入制导律

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
KR20120045258A (ko) * 2010-10-29 2012-05-09 삼성중공업 주식회사 선박용 프로펠러 설계방법, 이에 의하여 설계된 선박용 프로펠러, 및 이에 의하여 설계된 선박용 프로펠러를 갖는 선박
CN102880187A (zh) * 2012-09-21 2013-01-16 北京控制工程研究所 一种跳跃式再入飞行器初次再入段横向制导方法
CN103863579A (zh) * 2014-03-31 2014-06-18 北京控制工程研究所 一种深空探测返回过程的预测校正制导方法
CN104035335A (zh) * 2014-05-27 2014-09-10 北京航空航天大学 基于高精度纵、横程解析预测方法的平稳滑翔再入制导律

Cited By (34)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108710383A (zh) * 2018-05-25 2018-10-26 哈尔滨工业大学 一种基于航点规划与跟踪的四旋翼飞行器自主降落控制方法
CN109002576B (zh) * 2018-06-11 2021-11-02 北京航空航天大学 一种线性高阶比例制导系统脱靶量的幂级数解
CN109002576A (zh) * 2018-06-11 2018-12-14 北京航空航天大学 一种线性高阶比例制导系统脱靶量的幂级数解
CN109190248B (zh) * 2018-09-03 2023-07-18 中国运载火箭技术研究院 一种用于滑翔飞行器的滑翔射程解析方法及解析系统
CN109190248A (zh) * 2018-09-03 2019-01-11 中国运载火箭技术研究院 一种用于滑翔飞行器的滑翔射程解析方法及解析系统
CN109446582A (zh) * 2018-09-29 2019-03-08 北京航空航天大学 一种考虑地球自转的高精度降阶平稳滑翔动力学建模方法
CN109446582B (zh) * 2018-09-29 2023-04-07 北京航空航天大学 一种考虑地球自转的高精度降阶平稳滑翔动力学建模方法
CN109062241A (zh) * 2018-09-29 2018-12-21 北京航空航天大学 基于线性伪谱模型预测控制的自主全射向再入制导方法
CN109508030B (zh) * 2018-11-27 2020-08-04 北京航空航天大学 一种考虑多禁飞区约束的协同解析再入制导方法
CN109270960A (zh) * 2018-12-05 2019-01-25 中南大学 基于Radau伪谱法的在线最优反馈再入制导方法
CN109708525A (zh) * 2018-12-12 2019-05-03 中国人民解放军陆军工程大学 一种导弹飞行弹道的解算方法、系统及终端设备
CN109740198A (zh) * 2018-12-14 2019-05-10 中国人民解放军国防科技大学 一种基于解析预测的滑翔飞行器三维再入制导方法
CN109740198B (zh) * 2018-12-14 2022-12-23 中国人民解放军国防科技大学 一种基于解析预测的滑翔飞行器三维再入制导方法
CN111351401A (zh) * 2018-12-21 2020-06-30 北京理工大学 应用于捷联导引头制导飞行器的防侧偏制导方法
CN111351401B (zh) * 2018-12-21 2022-12-23 北京理工大学 应用于捷联导引头制导飞行器的防侧偏制导方法
CN111397441A (zh) * 2019-01-03 2020-07-10 北京理工大学 带有捷联导引头的远程制导飞行器的全射程覆盖制导系统
CN110147521A (zh) * 2019-04-25 2019-08-20 北京航空航天大学 一种高超声速飞行器跳跃滑翔弹道解析求解方法
CN110733670A (zh) * 2019-11-05 2020-01-31 中国人民解放军国防科技大学 一种短航程低过载的再入轨迹设计方法
CN111366044A (zh) * 2019-12-29 2020-07-03 湖北航天飞行器研究所 一种平飞过渡段制导控制方法
US20210216688A1 (en) * 2020-01-13 2021-07-15 Microsoft Technology Licensing, Llc Configuring aerodynamic simulation of a virtual object
CN111306989A (zh) * 2020-03-12 2020-06-19 北京航空航天大学 一种基于平稳滑翔弹道解析解的高超声速再入制导方法
CN111306989B (zh) * 2020-03-12 2021-12-28 北京航空航天大学 一种基于平稳滑翔弹道解析解的高超声速再入制导方法
CN111580555A (zh) * 2020-05-13 2020-08-25 北京控制工程研究所 一种高超声速飞行器上升段分段自适应预测校正制导方法
CN112198894A (zh) * 2020-07-31 2021-01-08 北京理工大学 旋翼无人机自主移动降落制导方法及系统
CN112541305A (zh) * 2020-12-10 2021-03-23 中国航空工业集团公司沈阳飞机设计研究所 一种基于全局变量求导的飞机气动力特性分析方法
CN112541305B (zh) * 2020-12-10 2024-01-02 中国航空工业集团公司沈阳飞机设计研究所 一种基于全局变量求导的飞机气动力特性分析方法
CN112507466B (zh) * 2020-12-21 2023-05-09 中国人民解放军96901部队23分队 考虑地球自转速率的滑翔弹道随能量变化降阶解计算方法
CN112507466A (zh) * 2020-12-21 2021-03-16 中国人民解放军96901部队23分队 考虑地球自转速率的滑翔弹道随能量变化降阶解
CN113008080A (zh) * 2021-01-26 2021-06-22 河北汉光重工有限责任公司 一种基于刚性原理对近海目标进行火控解算方法
CN113008080B (zh) * 2021-01-26 2023-01-13 河北汉光重工有限责任公司 一种基于刚性原理对近海目标进行火控解算方法
CN114740894B (zh) * 2022-05-13 2022-08-26 北京航空航天大学 基于注意力机制与门控循环单元的飞行器制导方法和系统
CN114740894A (zh) * 2022-05-13 2022-07-12 北京航空航天大学 基于注意力机制与门控循环单元的飞行器制导方法和系统
CN115542746A (zh) * 2022-12-05 2022-12-30 北京航空航天大学 高超声速飞行器的能量管控再入制导方法及装置
CN115542746B (zh) * 2022-12-05 2023-02-24 北京航空航天大学 高超声速飞行器的能量管控再入制导方法及装置

Also Published As

Publication number Publication date
CN108036676B (zh) 2019-08-23

Similar Documents

Publication Publication Date Title
CN108036676B (zh) 一种基于三维再入弹道解析解的全射向自主再入制导方法
CN107966156B (zh) 一种适用于运载火箭垂直回收段的制导律设计方法
CN109508030A (zh) 一种考虑多禁飞区约束的协同解析再入制导方法
CN104035335B (zh) 基于高精度纵、横程解析预测方法的平稳滑翔再入制导方法
CN104392047B (zh) 一种基于平稳滑翔弹道解析解的快速弹道规划方法
CN110908396B (zh) 可重复使用运载器全阶段再入返回制导方法
Slegers et al. Model predictive control of a parafoil and payload system
CN103838914B (zh) 一种高超声速飞行器滑翔段弹道解析求解方法
CN106586033A (zh) 自适应分段的多段线性伪谱广义标控脱靶量再入制导方法
CN108387140A (zh) 一种考虑多个禁飞区约束的解析再入制导方法
CN113479347B (zh) 一种火箭垂直回收着陆段轨迹控制方法
CN109062241B (zh) 基于线性伪谱模型预测控制的自主全射向再入制导方法
CN107121929B (zh) 基于线性协方差模型预测控制的鲁棒再入制导方法
CN107861517A (zh) 基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法
CN110015446B (zh) 一种半解析的火星进入制导方法
CN104648695A (zh) 一种基于倾侧角可用性的再入走廊最优规划方法
CN105352502B (zh) 一种微惯性航姿参考系统的姿态获取方法
Jwo et al. Development of a strapdown inertial navigation system simulation platform
CN106643341A (zh) 基于准平衡滑翔原理的力热控制耦合设计方法
CN107270896A (zh) 一种行人定位与轨迹跟踪方法和系统
CN110186478A (zh) 用于捷联式惯导系统的惯性传感器选型方法及系统
CN115265532A (zh) 一种用于船用组合导航中的辅助滤波方法
CN105973237B (zh) 基于实际飞行数据插值的仿真动态轨迹解析生成方法
Brunner et al. Comparison of numerical predictor-corrector and Apollo skip entry guidance algorithms
CN111488646A (zh) 一种旋转地球下高超声速平稳滑翔弹道的解析求解方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant