CN114357807A - 一种用于大气层外拦截的最优中制导方法及装置 - Google Patents

一种用于大气层外拦截的最优中制导方法及装置 Download PDF

Info

Publication number
CN114357807A
CN114357807A CN202210234975.5A CN202210234975A CN114357807A CN 114357807 A CN114357807 A CN 114357807A CN 202210234975 A CN202210234975 A CN 202210234975A CN 114357807 A CN114357807 A CN 114357807A
Authority
CN
China
Prior art keywords
interception
vector
terminal
equation
missile
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
CN202210234975.5A
Other languages
English (en)
Other versions
CN114357807B (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 CN202210234975.5A priority Critical patent/CN114357807B/zh
Publication of CN114357807A publication Critical patent/CN114357807A/zh
Application granted granted Critical
Publication of CN114357807B publication Critical patent/CN114357807B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • 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
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Geometry (AREA)
  • Operations Research (AREA)
  • Evolutionary Computation (AREA)
  • Algebra (AREA)
  • Computer Hardware Design (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Aiming, Guidance, Guns With A Light Source, Armor, Camouflage, And Targets (AREA)

Abstract

本发明提供了一种用于大气层外拦截的最优中制导方法及装置，其中，该用于大气层外拦截的最优中制导方法包括：基于拦截弹和目标导弹的轨道要素，建立终端时间可变的兰伯特问题的数学模型，基于建立的数学模型推导拦截弹期望关机速度方向矢量方程；基于拦截弹角动量、拦截弹性能参数以及预设的假设条件，预测助推段终端结束时刻的位置矢量和速度矢量；采用牛顿迭代法，基于兰伯特问题数学模型以及拦截弹期望关机速度方向矢量方程，求解拦截弹期望关机速度方向矢量；基于助推段终端位置矢量和速度矢量以及拦截弹期望关机速度方向矢量，求解拦截弹最优制导指令。可以提高拦截精度。

Description

一种用于大气层外拦截的最优中制导方法及装置
技术领域
本发明涉及制导控制技术领域，具体而言，涉及一种用于大气层外拦截的最优中制导方法及装置。
背景技术
随着武器技术的不断发展，导弹已成为现代战争中最为重要的攻击武器之一。为了打赢现代战争，导弹防御技术被越来越多的国家所重视。目前导弹为了追求更大的射程，绝大多数均采用弹道式飞行，导弹的射程越远，弹道的最高点会越高，导弹在大气层外被动飞行的时间也越长。导弹在大气层外的被动飞行段中，由于主发动机已经关闭，且不受气动力，因而，沿着开普勒轨道飞行，此时导弹的过载较小，飞行轨迹容易预测，是反导拦截的重要阶段。
目前，中段反导拦截的大部分过程都是无动力惯性滑行段，因此，拦截弹的关机点参数将直接影响拦截精度。现有用于大气层外拦截的解析制导律，是基于弹目重力为剩余飞行时间的二次函数的假设推导得到的，在滑行段距离较远时，该假设的精度将越来越低，直接影响拦截弹助推段关机点参数的选择，进而影响拦截精度。而目前几乎所有的反导拦截弹均使用固体火箭发动机，其受限于结构质量和复杂度，无法实现随时关机，因此其关机点速度大小一般无法控制，只能控制速度方向。综上所述，为提高使用固体火箭发动机的导弹拦截能力，需要设计一种适合固体火箭发动机的大气层外拦截最优中制导律。
发明内容
有鉴于此，本发明的目的在于提供用于大气层外拦截的最优中制导方法及装置，应用于装载有固体火箭发动机的拦截弹，以提高拦截精度。
第一方面，本发明实施例提供了用于大气层外拦截的最优中制导方法，包括：
基于拦截弹和目标导弹的轨道要素，建立终端时间可变的兰伯特问题的数学模型，基于建立的数学模型推导拦截弹期望关机速度方向矢量方程；
基于拦截弹角动量、拦截弹性能参数以及预设的假设条件，预测助推段终端结束时刻的位置矢量和速度矢量；
采用牛顿迭代法，基于兰伯特问题数学模型以及拦截弹期望关机速度方向矢量方程，求解拦截弹期望关机速度方向矢量；
基于助推段终端位置矢量和速度矢量以及拦截弹期望关机速度方向矢量，求解拦截弹最优制导指令。
结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方式，其中，所述基于拦截弹和目标导弹的轨道要素，建立终端时间可变的兰伯特问题的数学模型，推导拦截弹期望关机速度方向矢量方程，包括：
基于当前制导时刻的位置矢量及速度矢量，获取目标的轨道要素；
基于轨道要素中的终端真近点角、近地点幅角、升交点赤经以及轨道倾角，获取预测拦截点的位置矢量以及地心距的表达式，以及，基于当前时刻以及轨道要素中的偏近点角，获取预测拦截时刻；
基于关机点位置矢量和速度矢量，获取拦截弹的终端位置矢量和速度矢量表达式，基于拦截弹的终端位置矢量和速度矢量表达式，构建拦截弹兰伯特问题的数学模型。
结合第一方面的第一种可能的实施方式，本发明实施例提供了第一方面的第二种可能的实施方式，其中，所述基于拦截弹的终端位置矢量和速度矢量表达式，构建拦截弹兰伯特问题的数学模型，包括：
依据轨道力学知识，对拦截弹的终端位置矢量和速度矢量表达式中的拉格朗日系数表达式进行变换，依据变换前后的拉格朗日系数表达式相等的原则，得到拉格朗日系数方程；
基于轨道力学知识、当前时刻、预测拦截时刻以及拉格朗日系数方程，构建终端时间可变的兰伯特问题数学模型，该模型为包含拦截弹关机点弹道倾角、滑行段偏近点角增量和终端真近点角的代数方程组。
结合第一方面的第二种可能的实施方式，本发明实施例提供了第一方面的第三种可能的实施方式，其中，所述基于建立的数学模型推导拦截弹期望关机速度方向矢量方程，包括：
基于终端时间可变的兰伯特问题数学模型，获取拦截弹关机点弹道倾角、滑行段偏近点角增量和终端真近点角；
基于获取的拦截弹关机点弹道倾角、滑行段偏近点角增量、终端真近点角以及拦截弹的终端位置矢量、关机点位置矢量和速度表达式，获取拦截器期望关机速度方向矢量方程。
结合第一方面、第一方面的第一种可能的实施方式至第三种可能的实施方式中的任一种可能的实施方式，本发明实施例提供了第一方面的第四种可能的实施方式，其中，所述基于拦截弹角动量、拦截弹性能参数以及预设的假设条件，预测助推段终端结束时刻的位置矢量和速度矢量，包括：
基于拦截弹在地心惯性坐标系下的角动量，获取以轨道倾角和升交点赤经定义的轨道平面；
在轨道平面内，基于拦截弹性能参数，构建助推段动力学模型；
基于预设的假设条件，对助推段动力学模型中拦截弹速度导数进行积分，得到导弹助推段终端速度方程；
构建表征助推段弹道倾角的中间变量，对其微分方程进行积分，得到助推段结束时刻中间变量及助推速度积分项方程；
基于中间变量，获取助推段终端弹道倾角，基于助推段终端弹道倾角以及速度，获取与助推段速度矢量相关的积分项；
使用N个节点的高斯-勒让德积分公式，对与助推段速度矢量相关的积分项进行求解，得到与助推段速度矢量相关的积分项值；
以无量纲时间为自变量，对拦截弹弹道倾角进行拟合，得到拦截弹弹道倾角拟合方程；
基于拦截弹弹道倾角拟合方程以及助推段动力学模型，获取导弹助推段终点的位置和导弹位置矢量与地心到升交点连线的夹角；
基于导弹助推段终点的位置、导弹位置矢量与地心到升交点连线的夹角、轨道平面、助推段终端弹道倾角以及导弹助推段终端速度方程，确定拦截弹助推段结束时刻的位置矢量和速度矢量。
结合第一方面、第一方面的第一种可能的实施方式至第三种可能的实施方式中的任一种可能的实施方式，本发明实施例提供了第一方面的第五种可能的实施方式，其中，所述采用牛顿迭代法，基于兰伯特问题数学模型以及拦截弹期望关机速度方向矢量方程，求解拦截弹期望关机速度方向矢量，包括：
基于兰伯特问题数学模型的因变量数组分量，获取对拦截弹关机点弹道倾角、滑行段偏近点角的增量和终端真近点角的偏导数方程；
依据因变量数组各分量对拦截弹关机点弹道倾角、滑行段偏近点角的增量和终端真近点角的偏导数方程，获取兰伯特问题数学模型的雅克比矩阵；
基于雅克比矩阵，使用牛顿迭代法，获取兰伯特问题数学模型的迭代公式；
基于迭代终止条件以及迭代变量初始值，对兰伯特问题数学模型进行迭代运算，得到拦截弹期望关机速度方向矢量。
结合第一方面的第五种可能的实施方式，本发明实施例提供了第一方面的第六种可能的实施方式，其中，获取所述迭代变量初始值，包括：
采用预设的零弹目重力差模型获取预测拦截时刻表达式；
基于目标当前时刻的位置矢量和速度矢量，获取拦截弹关机时刻目标位置矢量和速度矢量方程；
基于拦截弹关机时刻目标位置矢量和速度矢量方程以及预测拦截时刻表达式，获取预测拦截时刻的初始猜想值；
基于预测拦截时刻的初始猜想值、目标由当前位置到终端位置的转移时间，获取终端真近点角的初始猜想值。
第二方面，本发明实施例还提供了一种用于大气层外拦截的最优中制导装置，包括：
期望矢量方程构建模块，用于基于拦截弹和目标导弹的轨道要素，建立终端时间可变的兰伯特问题的数学模型，基于建立的数学模型推导拦截弹期望关机速度方向矢量方程；
矢量预测模块，用于基于拦截弹角动量、拦截弹性能参数以及预设的假设条件，预测助推段终端结束时刻的位置矢量和速度矢量；
期望矢量求解模块，用于采用牛顿迭代法，基于兰伯特问题数学模型以及拦截弹期望关机速度方向矢量方程，求解拦截弹期望关机速度方向矢量；
最优指令求解模块，用于基于助推段终端位置矢量和速度矢量以及拦截弹期望关机速度方向矢量，求解拦截弹最优制导指令。
第三方面，本申请实施例提供了一种计算机设备，包括存储器、处理器及存储在所述存储器上并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述方法的步骤。
第四方面，本申请实施例提供了一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序，所述计算机程序被处理器运行时执行上述的方法的步骤。
本发明实施例提供的用于大气层外拦截的最优中制导方法及装置，通过基于拦截弹和目标导弹的轨道要素，建立终端时间可变的兰伯特问题的数学模型，基于建立的数学模型推导拦截弹期望关机速度方向矢量方程；基于拦截弹角动量、拦截弹性能参数以及预设的假设条件，预测助推段终端结束时刻的位置矢量和速度矢量；采用牛顿迭代法，基于兰伯特问题数学模型以及拦截弹期望关机速度方向矢量方程，求解拦截弹期望关机速度方向矢量；基于助推段终端位置矢量和速度矢量以及拦截弹期望关机速度方向矢量，求解拦截弹最优制导指令。这样，通过推导拦截弹期望关机速度方向矢量方程，基于牛顿迭代法求解拦截弹最优制导指令，能够使得拦截弹在助推段内推力矢量的偏角均保持在较小量级，从而提升拦截精度。
为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。
图1示出了本发明实施例所提供的用于大气层外拦截的最优中制导方法流程示意图；
图2示出了本发明实施例所提供的终端时间可变的兰伯特问题的数学模型示意图；
图3示出了本发明实施例所提供的轨道要素示意图；
图4示出了本发明实施例所提供的推力矢量方向示意图；
图5示出了本发明实施例所提供的用于大气层外拦截的最优中制导装置结构示意图；
图6为本申请实施例提供的一种计算机设备600的结构示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明实施例提供了一种用于大气层外拦截的最优中制导方法及装置，下面通过实施例进行描述。
图1示出了本发明实施例所提供的用于大气层外拦截的最优中制导方法流程示意图。如图1所示，该方法包括：
步骤101，基于拦截弹和目标导弹的轨道要素，建立终端时间可变的兰伯特问题的数学模型，基于建立的数学模型推导拦截弹期望关机速度方向矢量方程；
本发明实施例中，使用拦截弹和目标导弹的轨道要素，将拦截弹的终端位置矢量和速度矢量表示为关机点位置矢量和速度矢量的线性组合，从而推导建立目标开普勒轨道及拦截弹兰伯特问题的数学模型，得到拦截弹期望关机速度方向矢量方程。
本发明实施例中，作为一可选实施例，基于拦截弹和目标导弹的轨道要素，建立终端时间可变的兰伯特问题的数学模型，推导拦截弹期望关机速度方向矢量方程，包括：
A11，基于当前制导时刻的位置矢量及速度矢量，获取目标的轨道要素；
图2示出了本发明实施例所提供的终端时间可变的兰伯特问题的数学模型示意图。如图2所示，根据轨道力学的相关知识，沿开普勒轨道运行的目标（包括拦截弹以及目标导弹）的当前位置矢量信息，可由如下轨道要素完全定义：
Figure F_220228184538052_052705001
（1）
式中，
Figure F_220228184538130_130837002
为升交点赤经，
Figure F_220228184538266_266589003
为轨道倾角，
Figure F_220228184538461_461912004
是轨道半长轴,
Figure F_220228184538620_620099005
是轨道偏心率，
Figure F_220228184538682_682619006
为近地点幅角，
Figure F_220228184538776_776350007
为目标由近地点到当前点所经过的时间。
图3示出了本发明实施例所提供的轨道要素示意图。如图3所示，作为一可选实施例，各轨道要素可通过目标当前制导时刻的位置矢量及速度矢量计算得到。
A12，基于轨道要素中的终端真近点角、近地点幅角、升交点赤经以及轨道倾角，获取预测拦截点的位置矢量以及地心距的表达式，以及，基于当前时刻以及轨道要素中的偏近点角，获取预测拦截时刻；
本发明实施例中，目标的真近点角定义为沿运动方向测得的由轨道近地点到目标点的角度。根据开普勒轨道方程，预测拦截点即目标的位置矢量以及其地心距的表达式为：
Figure F_220228184538855_855931008
（2）
式中，
Figure F_220228184538934_934047009
为目标的位置矢量，
Figure F_220228184539031_031228010
为地心距，
Figure F_220228184539125_125035011
为终端真近点角；
Figure F_220228184539205_205552012
为终端位置矢量，
Figure F_220228184539283_283687013
、
Figure F_220228184539361_361814014
、
Figure F_220228184539426_426727015
分别为目标在地心惯性坐标系下的三个分量。
本发明实施例中，基于终端真近点角、近地点幅角、升交点赤经以及轨道倾角，获取拦截点的终端位置矢量：
Figure F_220228184539520_520480016
（3）
本发明实施例中，根据目标的开普勒轨道方程，目标由当前位置到终端位置的转移时间可以表示为：
Figure F_220228184539616_616709017
（4）
式中：
Figure F_220228184539710_710433018
为当前时刻，预测拦截时刻
Figure F_220228184539807_807644019
，
Figure F_220228184539933_933083020
为转移时间，
Figure F_220228184540013_013206021
为地心引力常数。
其中，偏近点角与终端真近点角的几何关系为：
Figure F_220228184540091_091294022
（5）
通过式（2）和（4）可知，预测拦截点
Figure F_220228184540169_169410023
由预测拦截时刻
Figure F_220228184540251_251944024
唯一确定。然而，式（4）所示的预测拦截时刻
Figure F_220228184540314_314445025
与偏近点角
Figure F_220228184540411_411106026
之间的关系式为超越方程。因此，无法由
Figure F_220228184540489_489229027
解析求解
Figure F_220228184540586_586380028
。
A13，基于关机点位置矢量和速度矢量，获取拦截弹的终端位置矢量和速度矢量表达式，基于拦截弹的终端位置矢量和速度矢量表达式，构建拦截弹兰伯特问题的数学模型。
本发明实施例中，由于开普勒轨道被限制在一个平面内，则拦截弹的终端位置矢量
Figure F_220228184540665_665043029
和速度矢量
Figure F_220228184540727_727529030
表达式可以表示为关机点位置矢量
Figure F_220228184540808_808085031
和速度矢量
Figure F_220228184540886_886251032
的线性组合：
Figure F_220228184540964_964336033
（6）
Figure F_220228184541045_045379034
和
Figure F_220228184541139_139152035
分别为
Figure F_220228184541226_226634036
和
Figure F_220228184541351_351552037
相对时间
Figure F_220228184541445_445318038
的导数。
Figure F_220228184541523_523457039
、
Figure F_220228184541604_604467040
、
Figure F_220228184541682_682618041
和
Figure F_220228184541760_760744042
为拉格朗日系数，具体表达式为：
Figure F_220228184541904_904775043
（7）
式中，
Figure F_220228184542078_078603044
为拦截弹滑行段偏近点角的增量，
Figure F_220228184542238_238772045
为拦截弹发动机关机时刻，
Figure F_220228184542332_332533046
为拦截弹轨道半长轴，其可根据关机点的状态量计算得到，表达式如下所示：
Figure F_220228184542413_413562047
（8）
本发明实施例中，作为一可选实施例，基于拦截弹的终端位置矢量和速度矢量表达式，构建拦截弹兰伯特问题的数学模型，包括：
A131，依据轨道力学知识，对拦截弹的终端位置矢量和速度矢量表达式中的拉格朗日系数表达式进行变换，依据变换前后的拉格朗日系数表达式相等的原则，得到拉格朗日系数方程；
本发明实施例中，由于关机速度大小
Figure F_220228184542507_507316048
及关机点地心距
Figure F_220228184542587_587359049
均是不可调的，因而，在求解过程中，可以将
Figure F_220228184542665_665551050
视为常值。这样，拉格朗日系数可以表示为：
Figure F_220228184542744_744139051
（9）
式中，
Figure F_220228184542844_844742052
为拦截弹轨道半通径，
Figure F_220228184542938_938469053
为拦截弹关机点位置矢量
Figure F_220228184543023_023892054
和终端位置矢量
Figure F_220228184543133_133284055
之间的夹角。
根据轨道力学知识有：
Figure F_220228184543198_198228056
（10）
Figure F_220228184543291_291973057
（11）
式中，
Figure F_220228184543390_390077058
为拦截弹关机点弹道倾角，
Figure F_220228184543483_483899059
、
Figure F_220228184543546_546896060
和
Figure F_220228184543671_671861061
分别为拦截弹关机点位置单位矢量在地心坐标系下的投影分量。
由式（7）及式（9）相等，可得如下拉格朗日系数方程：
Figure F_220228184543749_749971062
（12）
A132，基于轨道力学知识、当前时刻、预测拦截时刻以及拉格朗日系数方程，构建终端时间可变的兰伯特问题数学模型，该模型为包含拦截弹关机点弹道倾角、滑行段偏近点角增量和终端真近点角的代数方程组。
本发明实施例中，将式（10）及（4）代入上式并经过一系列代数运算，可以得到如下关于截弹关机点弹道倾角、滑行段偏近点角增量和终端真近点角的代数方程组：
Figure F_220228184543892_892585063
（13）
式中，
Figure F_220228184544005_005356064
为拦截弹发动机剩余工作时间。
为了后续的推导、求解方便，可以将上述代数方程组记为：
Figure F_220228184544083_083494065
（14）
其中：
Figure F_220228184544161_161614066
为因变量数组，
Figure F_220228184544245_245109067
、
Figure F_220228184544323_323241068
、
Figure F_220228184544404_404768069
分别为其分量，
Figure F_220228184544482_482890070
为自变量数组。
本发明实施例中，作为一可选实施例，基于建立的数学模型推导拦截弹期望关机速度方向矢量方程，包括：
B11，基于终端时间可变的兰伯特问题数学模型，获取拦截弹关机点弹道倾角、滑行段偏近点角增量和终端真近点角；
B12，基于获取的拦截弹关机点弹道倾角、滑行段偏近点角增量、终端真近点角以及拦截弹的终端位置矢量、关机点位置矢量和速度表达式，获取拦截器期望关机速度方向矢量方程。
上述方程组给出了目标开普勒轨道及拦截弹兰伯特问题的数学模型。通过求解该方程组，可以确定拦截弹滑行阶段的期望轨道。
将
Figure F_220228184544545_545434071
,
Figure F_220228184544627_627429072
和
Figure F_220228184544705_705536073
的值代入式（6），由拉格朗日系数可得拦截器期望关机速度方向矢量方程为：
Figure F_220228184544785_785117074
（15）
式中，
Figure F_220228184544926_926259075
为拦截器期望关机速度方向矢量。
步骤102，基于拦截弹角动量、拦截弹性能参数以及预设的假设条件，预测助推段终端结束时刻的位置矢量和速度矢量；
本发明实施例中，作为一可选实施例，基于拦截弹角动量、拦截弹性能参数以及预设的假设条件，预测助推段终端结束时刻的位置矢量和速度矢量，包括：
A21，基于拦截弹在地心惯性坐标系下的角动量，获取以轨道倾角和升交点赤经定义的轨道平面；
本发明实施例中，假设拦截弹处于标控状态，即推力方向与当前速度矢量方向相同，此时导弹的轨道平面不会改变。可以用轨道倾角和升交点赤经定义轨道平面，其计算方法为：
Figure F_220228184545009_009257076
（16）
Figure F_220228184545087_087377077
（17）
式中，
Figure F_220228184545165_165505078
、
Figure F_220228184545245_245586079
和
Figure F_220228184545354_354973080
为地心惯性坐标系下的角动量分量，
Figure F_220228184545536_536130081
为角动量大小。
A22，在轨道平面内，基于拦截弹性能参数，构建助推段动力学模型；
本发明实施例中，在轨道平面内可以建立如下助推段动力学模型：
Figure F_220228184545728_728990082
（18）
Figure F_220228184545810_810054083
（19）
Figure F_220228184545888_888191084
（20）
Figure F_220228184545966_966321085
（21）
式中，
Figure F_220228184546047_047370086
为导弹推力大小，
Figure F_220228184546125_125490087
为导弹初始质量，
Figure F_220228184546203_203631088
为导弹质量流量，
Figure F_220228184546281_281726089
为导弹发动机工作时间，
Figure F_220228184546359_359859090
为导弹所受重力加速度大小，
Figure F_220228184546440_440440091
为导弹弹道倾角，
Figure F_220228184546518_518564092
为导弹位置，
Figure F_220228184546582_582482093
为导弹速度，
Figure F_220228184546661_661148094
为导弹位置矢量与地心到升交点连线的夹角。
A23，基于预设的假设条件，对助推段动力学模型中拦截弹速度导数进行积分，得到导弹助推段终端速度方程；
本发明实施例中，因为助推段时间较短，所以重力的轴向分量的变化也很小，可以忽略。在该假设条件下，助推段动力学模型中拦截弹速度的导数可以表示为：
Figure F_220228184546754_754881095
（22）
式中，
Figure F_220228184546882_882335096
和
Figure F_220228184546960_960462097
分别表示当前时刻的重力加速度和弹道倾角。
对上式积分，可得导弹助推段终端速度方程：
Figure F_220228184547350_350605098
（23）
式中，
Figure F_220228184547526_526375099
为导弹助推段终端速度，
Figure F_220228184547671_671865100
为导弹当前时刻的速度，
Figure F_220228184547950_950703101
是导弹主发动机剩余工作时间，
Figure F_220228184548078_078629102
和
Figure F_220228184548536_536065103
分别为速度无量纲因子和时间无量纲因子，
Figure F_220228184548817_817876104
为初始轴向重推比，其表达式分别为：
Figure F_220228184549187_187467105
（24）
A24，构建表征助推段弹道倾角的中间变量，对其微分方程进行积分，得到助推段结束时刻中间变量及助推速度积分项方程；
本发明实施例中，为简化推导，定义中间变量
Figure F_220228184549297_297355106
为：
Figure F_220228184549411_411122107
（25）
将上式代入式（19），并忽略助推段地心距变化的影响，可得：
Figure F_220228184549520_520510108
（26）
对上式积分，可得助推段结束时刻
Figure F_220228184549654_654304109
的值的方程为：
Figure F_220228184549831_831545110
（27）
式中，
Figure F_220228184549972_972169111
为
Figure F_220228184550288_288081112
的初值，
Figure F_220228184550566_566849113
和
Figure F_220228184551108_108408114
是与
Figure F_220228184551243_243645115
有关的积分项，其表达式分别为：
Figure F_220228184551403_403336116
（28）
A25，基于中间变量，获取助推段终端弹道倾角，基于助推段终端弹道倾角以及速度，获取与助推段速度矢量相关的积分项；
通过式（25）可知助推段终端弹道倾角方程可表示为：
Figure F_220228184551559_559580117
（29）
将式（23）代入式（28）可得与助推段速度矢量相关的积分项：
Figure F_220228184551733_733406118
（30）
Figure F_220228184551956_956065119
（31）
式中，
Figure F_220228184552131_131343120
为无量纲化的发动机剩余工作时间，
Figure F_220228184552268_268057121
为无量纲化的导弹当前速度。
A26，使用N个节点的高斯-勒让德积分公式，对与助推段速度矢量相关的积分项进行求解，得到与助推段速度矢量相关的积分项值；
本发明实施例中，虽然式（31）中的积分项较为简单，但也无法直接求得解析解，使用N个节点的高斯-勒让德积分公式对与助推段速度矢量相关的积分项进行求解，可得积分项值：
Figure F_220228184552538_538089122
（32）
式中，
Figure F_220228184552761_761731123
为第i个积分权重，
Figure F_220228184552889_889165124
是第i个高斯节点处的无量纲化时间。一般情况下，当
Figure F_220228184553069_069820125
时高斯-勒让德积分的精度就已经足够。
A27，以无量纲时间为自变量，对拦截弹弹道倾角进行拟合，得到拦截弹弹道倾角拟合方程；
本发明实施例中，因为助推段动力学模型中的拦截弹位置公式（式（20））和拦截弹位置矢量与地心到升交点连线的夹角公式（式（21））中包含拦截弹弹道倾角的三角函数，所以很难直接求得拦截弹位置和拦截弹位置矢量与地心到升交点连线的夹角的解析解。但是因为助推段的时间较短，在助推段内拦截弹弹道倾角的变化也较小，可以分别使用以无量纲时间为自变量的二次函数的形式对拦截弹弹道倾角的正弦和余弦函数进行拟合，拦截弹弹道倾角拟合方程如下：
Figure F_220228184553288_288578126
（33）
式中，
Figure F_220228184553383_383266127
和
Figure F_220228184553446_446299128
分别为拟合系数。
通过端点条件可以解得：
Figure F_220228184553540_540036129
（34）
A28，基于拦截弹弹道倾角拟合方程以及助推段动力学模型，获取导弹助推段终点的位置和导弹位置矢量与地心到升交点连线的夹角；
本发明实施例中，将式（34）代入式（20）和式（21）并积分，可得导弹助推段终点的位置
Figure F_220228184553620_620126130
和导弹位置矢量与地心到升交点连线的夹角
Figure F_220228184553698_698248131
为：
Figure F_220228184553776_776376132
（35）
Figure F_220228184553901_901368133
（36）
式中，
Figure F_220228184553998_998041134
为初始时刻拦截弹位置矢量大小，
Figure F_220228184554076_076174135
为助推段平均地心距。
A29，基于导弹助推段终点的位置、导弹位置矢量与地心到升交点连线的夹角、轨道平面、助推段终端弹道倾角以及导弹助推段终端速度方程，确定拦截弹助推段结束时刻的位置矢量和速度矢量。
本发明实施例中，根据式（16）、式（17）、式（23）、式（29）、式（35）和式（36），可以确定拦截弹助推段结束时刻的位置矢量和速度矢量。
步骤103，采用牛顿迭代法，基于兰伯特问题数学模型以及拦截弹期望关机速度方向矢量方程，求解拦截弹期望关机速度方向矢量；
本发明实施例中，为提高求解速度，可通过推导代数方程组（兰伯特问题数学模型）的雅可比矩阵，并预测求解变量的初始猜想值，最终通过牛顿迭代法求解该代数方程组，得到拦截弹期望关机速度方向矢量。
本发明实施例中，作为一可选实施例，采用牛顿迭代法，基于兰伯特问题数学模型以及拦截弹期望关机速度方向矢量方程，求解拦截弹期望关机速度方向矢量，包括：
A31，基于兰伯特问题数学模型的因变量数组分量，获取对拦截弹关机点弹道倾角、滑行段偏近点角的增量和终端真近点角的偏导数方程；
本发明实施例中，对于式（13），迭代变量
Figure F_220228184554138_138675136
及
Figure F_220228184554219_219237137
是显含于函数表达式中。因此，易得
Figure F_220228184554297_297374138
对
Figure F_220228184554393_393033139
及
Figure F_220228184554471_471700140
的偏导数方程分别为：
Figure F_220228184554549_549812141
（37）
Figure F_220228184554646_646965142
（38）
本发明实施例中，由式（2）及式（11）可知，
Figure F_220228184554771_771982143
及
Figure F_220228184554853_853031144
均由
Figure F_220228184554931_931160145
决定，则
Figure F_220228184554996_996024146
对
Figure F_220228184555071_071786147
的偏导数可表示为：
Figure F_220228184555134_134263148
（39）
分别对式（2）及式（11）的等式两边取对终端真近点角的偏导数，可得：
Figure F_220228184555215_215821149
（40）
Figure F_220228184555309_309581150
（41）
式中，
Figure F_220228184555390_390585151
、
Figure F_220228184555484_484879152
、
Figure F_220228184555565_565430153
分别为关机点位置矢量
Figure F_220228184555645_645513154
在地心坐标系下的投影分量。与式（37）-（39）类似，
Figure F_220228184555723_723636155
对迭代变量
Figure F_220228184555804_804689156
,
Figure F_220228184555882_882812157
及
Figure F_220228184555960_960941158
的偏导数可表示为：
Figure F_220228184556136_136725159
（42）
由式（5），得
Figure F_220228184556233_233401160
对
Figure F_220228184556311_311534161
的偏导数为：
Figure F_220228184556392_392537162
（43）
函数
Figure F_220228184556471_471195163
对
Figure F_220228184556549_549332164
,
Figure F_220228184556714_714851165
及
Figure F_220228184556862_862802166
求偏导，易得
Figure F_220228184557011_011738167
对这三个变量的偏导数为：
Figure F_220228184557105_105948168
（44）
A32，依据因变量数组各分量对拦截弹关机点弹道倾角、滑行段偏近点角的增量和终端真近点角的偏导数方程，获取兰伯特问题数学模型的雅克比矩阵；
由式（37）-（39）、式（42）、式（44）可知，代数方程组
Figure F_220228184557216_216318169
的雅克比矩阵为：
Figure F_220228184557294_294421170
（45）
A33，基于雅克比矩阵，使用牛顿迭代法，获取兰伯特问题数学模型的迭代公式；
使用牛顿迭代法求解方程组（14）的迭代公式为：
Figure F_220228184557406_406238171
（46）
其中，上标
Figure F_220228184557500_500008172
代表第
Figure F_220228184557562_562491173
次迭代的变量。
A34，基于迭代终止条件以及迭代变量初始值，对兰伯特问题数学模型进行迭代运算，得到拦截弹期望关机速度方向矢量。
本发明实施例中，如果相邻两次迭代变量的差小于给定的收敛阀值
Figure F_220228184557642_642579174
则结束迭代，即迭代终止条件为：
Figure F_220228184557720_720710175
（47）
本发明实施例中，作为一可选实施例，获取迭代变量初始值，包括：
A341，采用预设的零弹目重力差模型获取预测拦截时刻表达式；
由式（46）可知，为了启动牛顿迭代，需要迭代变量的一个初始猜想值。由式（4）及式（5）易知，目标的真近点角与预测拦截时间有关。因此，有必要估计预测拦截时刻的值。采用零弹目重力差模型估计预测拦截时刻，其表达式为：
Figure F_220228184557802_802242176
（48）
式中，
Figure F_220228184557880_880380177
及
Figure F_220228184557958_958500178
分别为拦截弹关机时刻弹目相对位置及相对速度矢量，右上角标T表示矩阵的转置。
A342，基于目标当前时刻的位置矢量和速度矢量，获取拦截弹关机时刻目标位置矢量和速度矢量方程；
本发明实施例中，对于拦截弹关机时刻自身位置矢量和速度矢量已通过步骤102预测得到，对于拦截弹关机时刻目标位置矢量和速度矢量，对应的拦截弹关机时刻目标位置矢量和速度矢量方程可近似为：
Figure F_220228184558040_040533179
（49）
式中，
Figure F_220228184558134_134285180
、
Figure F_220228184558202_202145181
分别为目标当前时刻的位置矢量和速度矢量。
A343，基于拦截弹关机时刻目标位置矢量和速度矢量方程以及预测拦截时刻表达式，获取预测拦截时刻的初始猜想值；
A344，基于预测拦截时刻的初始猜想值、目标由当前位置到终端位置的转移时间，获取终端真近点角的初始猜想值。
本发明实施例中，在获得预测拦截时刻的初始猜想值后，便可根据式（4）及式（5）计算终端真近点角的初始猜想值。
本发明实施例中，作为一可选实施例，基于预测拦截时刻的初始猜想值、目标由当前位置到终端位置的转移时间，获取终端真近点角的初始猜想值，包括：
A3441，基于目标由当前位置到终端位置的转移时间，获取平近点角；
本发明实施例中，由于式（4）所示的方程为关于偏近点角
Figure F_220228184558280_280286182
的超越方程，因此不存在解析解，只能通过数值迭代计算偏近点角的值。
本发明实施例中，定义平近点角
Figure F_220228184558391_391072183
，由式（4）易得，其值可由
Figure F_220228184558469_469733184
计算得到：
Figure F_220228184558532_532229185
（50）
A3442，基于平近点角，利用迭代方法，对目标终端偏近点角进行迭代运算，得到目标终端偏近点角；
本发明实施例中，可由以下迭代公式计算目标终端偏近点角的值：
Figure F_220228184558627_627449186
（51）
与式（46）同理，上标
Figure F_220228184558705_705550187
代表第
Figure F_220228184558786_786579188
次迭代的变量。
A3443，基于目标终端偏近点角，依据目标终端偏近点角与终端真近点角的几何关系，得到终端真近点角的初始猜想值。
本发明实施例中，获得
Figure F_220228184558865_865247189
的初始猜想值后, 相应的
Figure F_220228184558943_943372190
的初始猜想值可由式（5）计算得到：
Figure F_220228184559024_024898191
（52）
本发明实施例中，作为一可选实施例，该方法还包括：
A3444，依据助推段终端弹道倾角方程，获取助推段终端弹道倾角的初始猜想值；
A3445，依据助推段终端弹道倾角的初始猜想值、终端真近点角的初始猜想值以及因变量数组分量，获取
Figure F_220228184559103_103036192
的初始猜想值。
本发明实施例中，
Figure F_220228184559182_182582193
的初值猜想值可选为通过式（29）计算得到的名义值。至于
Figure F_220228184559261_261234194
的初始猜想值，可由
Figure F_220228184559339_339338195
和
Figure F_220228184559420_420405196
的初始值求解
Figure F_220228184559514_514164197
获得，其表达式为：
Figure F_220228184559592_592297198
（53）
式中，
Figure F_220228184559723_723642199
和
Figure F_220228184559803_803729200
可根据
Figure F_220228184559882_882320201
的初始猜想值利用式（2）及（11）计算得到。
需要注意的是，本发明实施例给出的迭代变量的初始猜想值只在第一个制导周期使用，之后的制导周期可采用上一周期计算得到的结果作为迭代初值。由于上一制导周期的最优解与此制导周期的最优解很接近，这样可以有效的减少迭代次数，提高计算速度。
步骤104，基于助推段终端位置矢量和速度矢量以及拦截弹期望关机速度方向矢量，求解拦截弹最优制导指令。
本发明实施例中，建立弹道倾角及航向角的摄动方程，选择性能泛函，根据一阶必要条件可解得协态变量和控制变量的函数，再通过高斯-勒让德积分公式求制导系数的近似解。
本发明实施例中，作为一可选实施例，基于助推段终端位置矢量和速度矢量以及拦截弹期望关机速度方向矢量，求解拦截弹最优制导指令，包括：
A41，基于关机点弹道倾角和关机点航向角，获取拦截弹关机点速度方向表达式；
本发明实施例中，拦截弹关机点速度方向可由关机点弹道倾角（助推段终端弹道倾角）
Figure F_220228184559960_960449202
和关机点航向角
Figure F_220228184600071_071791203
表示，其拦截弹关机点速度方向表达式分别为：
Figure F_220228184600231_231936204
（54）
式中，
Figure F_220228184600407_407724205
、
Figure F_220228184600550_550781206
和
Figure F_220228184600618_618666207
分别为拦截弹关机点速度单位矢量在地心坐标系下的投影分量。
A42，基于拦截弹关机点速度方向表达式以及拦截弹动力学方程，获取关机点弹道倾角和关机点航向角的摄动方程；
本发明实施例中，根据拦截弹动力学方程可得其弹道倾角和航向角的摄动方程为：
Figure F_220228184600696_696785208
（55）
式中：
Figure F_220228184600774_774914209
和
Figure F_220228184600856_856463210
分别为拦截弹推力方向单位矢量
Figure F_220228184600934_934580211
在垂直面和水平面内的投影分量。
当已知
Figure F_220228184601015_015143212
和
Figure F_220228184601093_093275213
后，则
Figure F_220228184601171_171398214
在当前速度方向的投影分量
Figure F_220228184601267_267589215
为：
Figure F_220228184601345_345722216
（56）
图4示出了本发明实施例所提供的推力矢量方向示意图。如图4所示，推力方向矢量
Figure F_220228184601444_444350217
三个分量对应方向的单位矢量
Figure F_220228184601522_522471218
、
Figure F_220228184601603_603506219
、
Figure F_220228184601681_681652220
可以分别表示为：
Figure F_220228184601775_775401221
（57）
式中，
Figure F_220228184601863_863254222
和
Figure F_220228184601972_972649223
分别为拦截弹当前位置矢量和速度矢量。
A43，依据推力矢量分量，获取推力方向矢量；
本发明实施例中，如果得到推力矢量分量
Figure F_220228184602053_053731224
和
Figure F_220228184602131_131833225
，则推力方向矢量
Figure F_220228184602196_196283226
可表示为：
Figure F_220228184602274_274426227
（58）
A44，获取关机点弹道倾角和关机点航向角的摄动方程需满足的终端约束，所述终端约束为基于拦截弹关机点弹道倾角和航向角的预测值的约束；
本发明实施例中，式（55）需要满足的终端约束为：
Figure F_220228184602368_368170228
（59）
式中，
Figure F_220228184602449_449226229
和
Figure F_220228184602527_527341230
分别为拦截弹关机点弹道倾角和航向角的期望值。
A45，基于最优推力矢量方向，获取性能指标函数；
本发明实施例中，为了获得稳定的控制输入，基于拦截弹最优推力方向矢量在垂直面和水平面内投影分量的平方和设计性能指标函数为：
Figure F_220228184602593_593712231
（60）
式中，
Figure F_220228184602719_719234232
为权重系数，权重项
Figure F_220228184602801_801755233
是用来对
Figure F_220228184602879_879889234
及
Figure F_220228184602942_942391235
曲线进行整形的。如果选择
Figure F_220228184603025_025395236
，那么拦截弹关机时对
Figure F_220228184603103_103515237
及
Figure F_220228184603219_219240238
的终端值将收敛为零。
A46，依据求解拦截弹最优制导指令的哈密顿函数以及关机点弹道倾角和关机点航向角的摄动方程，获取拦截弹关机时刻弹道倾角偏差及航向角偏差的表达式；
本发明实施例中，最优控制问题的哈密顿函数为：
Figure F_220228184603297_297369239
（61）
其中，
Figure F_220228184603393_393520240
和
Figure F_220228184603472_472176241
为协态变量。根据一阶必要条件，可得：
Figure F_220228184603550_550299242
（62）
Figure F_220228184603708_708497243
（63）
将上式代入式（55）并进行积分，可得拦截弹关机时刻弹道倾角偏差及航向角偏差的表达式为：
Figure F_220228184603909_909668244
（64）
A47，基于拦截弹关机时刻弹道倾角偏差及航向角偏差的表达式以及终端约束，获取协态量值表达式；
本发明实施例中，结合式（59）所示的终端约束，可求得协态量值表达式为：
Figure F_220228184604073_073251245
（65）
式中，
Figure F_220228184604166_166995246
和
Figure F_220228184604247_247566247
为中间变量，其表达式为：
Figure F_220228184604325_325694248
（66）
式中，
Figure F_220228184604422_422860249
,
Figure F_220228184604516_516614250
和
Figure F_220228184604597_597160251
分别为拦截弹当前时刻质量，速度大小及弹道倾角。其中拦截弹关机点速度
Figure F_220228184604675_675302252
可通过式（23）计算得到。
A48，使用以无量纲时间为自变量的二次函数，对助推段终端弹道倾角进行拟合，得到倾角拟合公式；
本发明实施例中，由于拦截弹助推段的时间较短，在助推段内弹道倾角的变化也较小，可以使用以无量纲时间
Figure F_220228184604769_769048253
为自变量的二次函数的形式对弹道倾角的余弦函数进行拟合，倾角拟合公式如下：
Figure F_220228184604833_833990254
（67）
式中，
Figure F_220228184604927_927739255
、
Figure F_220228184604994_994114256
、
Figure F_220228184605072_072750257
分别为拟合系数。
A49，基于导弹助推段终端速度方程、倾角拟合公式以及协态量值表达式，得到协态量中间参数表达式；
本发明实施例中，将式（23）和式（67）代入式（66），可得协态量中间参数表达式：
Figure F_220228184605150_150880258
（68）
A50，使用N个节点的高斯-勒让德积分公式，对协态量中间参数表达式进行求解，得到协态量中间参数值；
本发明实施例中，式（68）无法求得解析解，因而，使用N个节点的高斯-勒让德积分公式，对其进行求解，可得协态量中间参数值：
Figure F_220228184605269_269046259
（69）
其中，
Figure F_220228184605403_403773260
为第i个高斯积分节点
Figure F_220228184605517_517056261
对应的时刻，
Figure F_220228184605603_603511262
是第i个高斯节点处的无量纲化时间，
Figure F_220228184605681_681648263
为第i个高斯积分权重系数。
A51，基于协态量值表达式以及哈密顿函数，获取包括最优推力矢量方向的拦截弹最优制导指令。
将式（65）代入式（63）可得
Figure F_220228184605775_775398264
及
Figure F_220228184605855_855950265
的最优值为：
Figure F_220228184605934_934090266
（70）
式中，
Figure F_220228184606030_030286267
和
Figure F_220228184606124_124033268
分别为
Figure F_220228184606186_186522269
和
Figure F_220228184606272_272956270
的最优值；
Figure F_220228184606351_351076271
和
Figure F_220228184606436_436530272
分别为制导系数。其具体表达式为：
Figure F_220228184606530_530275273
（71）
至此，已经推导得出大气层外拦截弹最优制导指令。
本发明实施例的一种适合固体火箭发动机的大气层外拦截最优中制导方法，以大气层外远距离拦截为背景，利用开普勒轨道的性质和兰伯特问题的特点，将拦截弹终端位置矢量和速度矢量表示为关机点位置矢量和速度矢量的线性组合，基于固体火箭发动机无法随时关机的特点，预测关机点参数，通过推导代数方程组的雅克比矩阵简化牛顿迭代法的计算量，提升求解速度，最后依据最优控制问题的一阶必要条件推导得出拦截弹最优制导指令，使得拦截弹在助推段内推力矢量的偏角均保持在较小量级，进而降低对拦截弹控制系统的要求，提升拦截精度。
图5示出了本发明实施例所提供的用于大气层外拦截的最优中制导装置结构示意图。如图5所示，该装置包括：
期望矢量方程构建模块501，用于基于拦截弹和目标导弹的轨道要素，建立终端时间可变的兰伯特问题的数学模型，基于建立的数学模型推导拦截弹期望关机速度方向矢量方程；
本发明实施例中，期望矢量方程构建模块501包括：
要素获取单元（图中未示出），用于基于当前制导时刻的位置矢量及速度矢量，获取目标的轨道要素；
拦截点定位单元，用于基于轨道要素中的终端真近点角、近地点幅角、升交点赤经以及轨道倾角，获取预测拦截点的位置矢量以及地心距的表达式，以及，基于当前时刻以及轨道要素中的偏近点角，获取预测拦截时刻；
模型构建单元，用于基于关机点位置矢量和速度矢量，获取拦截弹的终端位置矢量和速度矢量表达式，基于拦截弹的终端位置矢量和速度矢量表达式，构建拦截弹兰伯特问题的数学模型。
本发明实施例中，作为一可选实施例，基于拦截弹的终端位置矢量和速度矢量表达式，构建拦截弹兰伯特问题的数学模型，包括：
依据轨道力学知识，对拦截弹的终端位置矢量和速度矢量表达式中的拉格朗日系数表达式进行变换，依据变换前后的拉格朗日系数表达式相等的原则，得到拉格朗日系数方程；
基于轨道力学知识、当前时刻、预测拦截时刻以及拉格朗日系数方程，构建终端时间可变的兰伯特问题数学模型，该模型为包含拦截弹关机点弹道倾角、滑行段偏近点角增量和终端真近点角的代数方程组。
本发明实施例中，作为一可选实施例，基于建立的数学模型推导拦截弹期望关机速度方向矢量方程，包括：
基于终端时间可变的兰伯特问题数学模型，获取拦截弹关机点弹道倾角、滑行段偏近点角增量和终端真近点角；
基于获取的拦截弹关机点弹道倾角、滑行段偏近点角增量、终端真近点角以及拦截弹的终端位置矢量、关机点位置矢量和速度表达式，获取拦截器期望关机速度方向矢量方程。
矢量预测模块502，用于基于拦截弹角动量、拦截弹性能参数以及预设的假设条件，预测助推段终端结束时刻的位置矢量和速度矢量；
本发明实施例中，作为一可选实施例，矢量预测模块502包括：
平面构建单元（图中未示出），用于基于拦截弹在地心惯性坐标系下的角动量，获取以轨道倾角和升交点赤经定义的轨道平面；
模型生成单元，用于在轨道平面内，基于拦截弹性能参数，构建助推段动力学模型；
模型求解单元，用于基于预设的假设条件，对助推段动力学模型中拦截弹速度导数进行积分，得到导弹助推段终端速度方程；
变换单元，用于构建表征助推段弹道倾角的中间变量，对其微分方程进行积分，得到助推段结束时刻中间变量方程；
积分项构建单元，用于基于中间变量，获取助推段终端弹道倾角，基于助推段终端弹道倾角以及速度，获取与助推段速度矢量相关的积分项；
积分项值获取单元，使用N个节点的高斯-勒让德积分公式，对与助推段速度矢量相关的积分项进行求解，得到与助推段速度矢量相关的积分项值；
拟合单元，用于以无量纲时间为自变量，对拦截弹弹道倾角进行拟合，得到拦截弹弹道倾角拟合方程；
夹角获取单元，用于基于拦截弹弹道倾角拟合方程以及助推段动力学模型，获取导弹助推段终点的位置和导弹位置矢量与地心到升交点连线的夹角；
位速矢量获取单元，用于基于导弹助推段终点的位置、导弹位置矢量与地心到升交点连线的夹角、轨道平面、助推段终端弹道倾角以及导弹助推段终端速度方程，确定拦截弹助推段结束时刻的位置矢量和速度矢量。
期望矢量求解模块503，用于采用牛顿迭代法，基于兰伯特问题数学模型以及拦截弹期望关机速度方向矢量方程，求解拦截弹期望关机速度方向矢量；
本发明实施例中，作为一可选实施例，期望矢量求解模块503包括：
偏导单元（图中未示出），用于基于兰伯特问题数学模型的因变量数组分量，获取对对拦截弹关机点弹道倾角、滑行段偏近点角的增量和终端真近点角的偏导数方程；
雅克比矩阵获取单元，用于依据因变量数组各分量对拦截弹关机点弹道倾角、滑行段偏近点角的增量和终端真近点角的偏导数方程，获取兰伯特问题数学模型的雅克比矩阵；
迭代单元，用于基于雅克比矩阵，使用牛顿迭代法，获取兰伯特问题数学模型的迭代公式；
期望矢量求解单元，用于基于迭代终止条件以及迭代变量初始值，对兰伯特问题数学模型进行迭代运算，得到拦截弹期望关机速度方向矢量。
本发明实施例中，作为一可选实施例，获取所述迭代变量初始值，包括：
采用预设的零弹目重力差模型获取预测拦截时刻表达式；
基于目标当前时刻的位置矢量和速度矢量，获取拦截弹关机时刻目标位置矢量和速度矢量方程；
基于拦截弹关机时刻目标位置矢量和速度矢量方程以及预测拦截时刻表达式，获取预测拦截时刻的初始猜想值；
基于预测拦截时刻的初始猜想值、目标由当前位置到终端位置的转移时间，获取终端真近点角的初始猜想值。
最优指令求解模块504，用于基于助推段终端位置矢量和速度矢量以及拦截弹期望关机速度方向矢量，求解拦截弹最优制导指令。
本发明实施例中，作为一可选实施例，最优指令求解模块504具体用于：
基于关机点弹道倾角和关机点航向角，获取拦截弹关机点速度方向表达式；
基于拦截弹关机点速度方向表达式以及拦截弹动力学方程，获取关机点弹道倾角和关机点航向角的摄动方程；
依据关机点弹道倾角和关机点航向角的摄动方程，获取最优推力矢量方向；
获取关机点弹道倾角和关机点航向角的摄动方程需满足的终端约束，所述终端约束为基于拦截弹关机点弹道倾角和航向角的预测值的约束；
基于最优推力矢量方向，获取性能指标函数；
依据求解拦截弹最优制导指令的哈密顿函数以及关机点弹道倾角和关机点航向角的摄动方程，获取拦截弹关机时刻弹道倾角偏差及航向角偏差的表达式；
基于拦截弹关机时刻弹道倾角偏差及航向角偏差的表达式以及终端约束，获取协态量值表达式；
使用以无量纲时间为自变量的二次函数，对助推段终端弹道倾角进行拟合，得到倾角拟合公式；
基于导弹助推段终端速度方程、倾角拟合公式以及协态量值表达式，得到协态量中间参数表达式；
使用N个节点的高斯-勒让德积分公式，对协态量中间参数表达式进行求解，得到协态量中间参数值；
基于协态量值表达式以及哈密顿函数，获取包括最优推力矢量方向的拦截弹最优制导指令。
如图6所示，本申请一实施例提供了一种计算机设备600，用于执行图1中的用于大气层外拦截的最优中制导方法，该设备包括存储器601、与存储器601通过总线相连的处理器602及存储在该存储器601上并可在该处理器602上运行的计算机程序，其中，上述处理器602执行上述计算机程序时实现上述用于大气层外拦截的最优中制导方法的步骤。
具体地，上述存储器601和处理器602能够为通用的存储器和处理器，这里不做具体限定，当处理器602运行存储器601存储的计算机程序时，能够执行上述用于大气层外拦截的最优中制导方法。
对应于图1中的用于大气层外拦截的最优中制导方法，本申请实施例还提供了一种计算机可读存储介质，该计算机可读存储介质上存储有计算机程序，该计算机程序被处理器运行时执行上述用于大气层外拦截的最优中制导方法的步骤。
具体地，该存储介质能够为通用的存储介质，如移动磁盘、硬盘等，该存储介质上的计算机程序被运行时，能够执行上述用于大气层外拦截的最优中制导方法。
在本申请所提供的实施例中，应该理解到，所揭露系统和方法，可以通过其它的方式实现。以上所描述的系统实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，又例如，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些通信接口，系统或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。
所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。
另外，在本申请提供的实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。
所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本申请的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备（可以是个人计算机，服务器，或者网络设备等）执行本申请各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器（Read-Only Memory， ROM）、随机存取存储器（Random Access Memory ，RAM）、磁碟或者光盘等各种可以存储程序代码的介质。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释，此外，术语“第一”、“第二”、“第三”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
最后应说明的是：以上所述实施例，仅为本申请的具体实施方式，用以说明本申请的技术方案，而非对其限制，本申请的保护范围并不局限于此，尽管参照前述实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本申请实施例技术方案的精神和范围。都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应所述以权利要求的保护范围为准。

Claims (10)

1.一种用于大气层外拦截最优中制导方法，其特征在于，包括：
基于拦截弹和目标导弹的轨道要素，建立终端时间可变的兰伯特问题的数学模型，基于建立的数学模型推导拦截弹期望关机速度方向矢量方程；
基于拦截弹角动量、拦截弹性能参数以及预设的假设条件，预测助推段终端结束时刻的位置矢量和速度矢量；
采用牛顿迭代法，基于兰伯特问题数学模型以及拦截弹期望关机速度方向矢量方程，求解拦截弹期望关机速度方向矢量；
基于助推段终端位置矢量和速度矢量以及拦截弹期望关机速度方向矢量，求解拦截弹最优制导指令。
2.根据权利要求1所述的方法，其特征在于，所述基于拦截弹和目标导弹的轨道要素，建立终端时间可变的兰伯特问题的数学模型，基于建立的数学模型推导拦截弹期望关机速度方向矢量方程，包括：
基于当前制导时刻的位置矢量及速度矢量，获取目标的轨道要素；
基于轨道要素中的终端真近点角、近地点幅角、升交点赤经以及轨道倾角，获取预测拦截点的位置矢量以及地心距的表达式，以及，基于当前时刻以及轨道要素中的偏近点角，获取预测拦截时刻；
基于关机点位置矢量和速度矢量，获取拦截弹的终端位置矢量和速度矢量表达式，基于拦截弹的终端位置矢量和速度矢量表达式，构建拦截弹兰伯特问题的数学模型。
3.根据权利要求2所述的方法，其特征在于，所述基于拦截弹的终端位置矢量和速度矢量表达式，构建拦截弹兰伯特问题的数学模型，包括：
依据轨道力学知识，对拦截弹的终端位置矢量和速度矢量表达式中的拉格朗日系数表达式进行变换，依据变换前后的拉格朗日系数表达式相等的原则，得到拉格朗日系数方程；
基于轨道力学知识、当前时刻、预测拦截时刻以及拉格朗日系数方程，构建终端时间可变的兰伯特问题数学模型，该模型为包含拦截弹关机点弹道倾角、滑行段偏近点角增量和终端真近点角的代数方程组。
4.根据权利要求3所述的方法，其特征在于，所述基于建立的数学模型推导拦截弹期望关机速度方向矢量方程，包括：
基于终端时间可变的兰伯特问题数学模型，获取拦截弹关机点弹道倾角、滑行段偏近点角增量和终端真近点角；
基于获取的拦截弹关机点弹道倾角、滑行段偏近点角增量、终端真近点角以及拦截弹的终端位置矢量、关机点位置矢量和速度表达式，获取拦截器期望关机速度方向矢量方程。
5.根据权利要求1至4任一项所述的方法，其特征在于，所述基于拦截弹角动量、拦截弹性能参数以及预设的假设条件，预测助推段终端结束时刻的位置矢量和速度矢量，包括：
基于拦截弹在地心惯性坐标系下的角动量，获取以轨道倾角和升交点赤经定义的轨道平面；
在轨道平面内，基于拦截弹性能参数，构建助推段动力学模型；
基于预设的假设条件，对助推段动力学模型中拦截弹速度导数进行积分，得到导弹助推段终端速度方程；
构建表征助推段弹道倾角的中间变量，对其微分方程进行积分，得到助推段结束时刻中间变量及助推速度积分项方程；
基于中间变量，获取助推段终端弹道倾角，基于助推段终端弹道倾角以及速度，获取与助推段速度矢量相关的积分项；
使用N个节点的高斯-勒让德积分公式，对与助推段速度矢量相关的积分项进行求解，得到与助推段速度矢量相关的积分项值；
以无量纲时间为自变量，对拦截弹弹道倾角进行拟合，得到拦截弹弹道倾角拟合方程；
基于拦截弹弹道倾角拟合方程以及助推段动力学模型，获取导弹助推段终点的位置和导弹位置矢量与地心到升交点连线的夹角；
基于导弹助推段终点的位置、导弹位置矢量与地心到升交点连线的夹角、轨道平面、助推段终端弹道倾角以及导弹助推段终端速度方程，确定拦截弹助推段结束时刻的位置矢量和速度矢量。
6.根据权利要求1至4任一项所述的方法，其特征在于，所述采用牛顿迭代法，基于兰伯特问题数学模型以及拦截弹期望关机速度方向矢量方程，求解拦截弹期望关机速度方向矢量，包括：
基于兰伯特问题数学模型的因变量数组分量，获取对拦截弹关机点弹道倾角、滑行段偏近点角的增量和终端真近点角的偏导数方程；
依据因变量数组各分量对拦截弹关机点弹道倾角、滑行段偏近点角的增量和终端真近点角的偏导数方程，获取兰伯特问题数学模型的雅克比矩阵；
基于雅克比矩阵，使用牛顿迭代法，获取兰伯特问题数学模型的迭代公式；
基于迭代终止条件以及迭代变量初始值，对兰伯特问题数学模型进行迭代运算，得到拦截弹期望关机速度方向矢量。
7.根据权利要求6所述的方法，其特征在于，获取所述迭代变量初始值，包括：
采用预设的零弹目重力差模型获取预测拦截时刻表达式；
基于目标当前时刻的位置矢量和速度矢量，获取拦截弹关机时刻目标位置矢量和速度矢量方程；
基于拦截弹关机时刻目标位置矢量和速度矢量方程以及预测拦截时刻表达式，获取预测拦截时刻的初始猜想值；
基于预测拦截时刻的初始猜想值、目标由当前位置到终端位置的转移时间，获取终端真近点角的初始猜想值。
8.一种用于大气层外拦截的最优中制导装置，其特征在于，包括：
期望矢量方程构建模块，用于基于拦截弹和目标导弹的轨道要素，建立终端时间可变的兰伯特问题的数学模型，基于建立的数学模型推导拦截弹期望关机速度方向矢量方程；
矢量预测模块，用于基于拦截弹角动量、拦截弹性能参数以及预设的假设条件，预测助推段终端结束时刻的位置矢量和速度矢量；
期望矢量求解模块，用于采用牛顿迭代法，基于兰伯特问题数学模型以及拦截弹期望关机速度方向矢量方程，求解拦截弹期望关机速度方向矢量；
最优指令求解模块，用于基于助推段终端位置矢量和速度矢量以及拦截弹期望关机速度方向矢量，求解拦截弹最优制导指令。
9.一种计算机设备，其特征在于，包括：处理器、存储器和总线，所述存储器存储有所述处理器可执行的机器可读指令，当计算机设备运行时，所述处理器与所述存储器之间通过总线通信，所述机器可读指令被所述处理器执行时执行如权利要求1至7任一所述的用于大气层外拦截的最优中制导方法的步骤。
10.一种计算机可读存储介质，其特征在于，该计算机可读存储介质上存储有计算机程序，该计算机程序被处理器运行时执行如权利要求1至7任一所述的用于大气层外拦截的最优中制导方法的步骤。
CN202210234975.5A 2022-03-11 2022-03-11 一种用于大气层外拦截的最优中制导方法及装置 Active CN114357807B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210234975.5A CN114357807B (zh) 2022-03-11 2022-03-11 一种用于大气层外拦截的最优中制导方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210234975.5A CN114357807B (zh) 2022-03-11 2022-03-11 一种用于大气层外拦截的最优中制导方法及装置

Publications (2)

Publication Number Publication Date
CN114357807A true CN114357807A (zh) 2022-04-15
CN114357807B CN114357807B (zh) 2022-05-27

Family

ID=81094690

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210234975.5A Active CN114357807B (zh) 2022-03-11 2022-03-11 一种用于大气层外拦截的最优中制导方法及装置

Country Status (1)

Country Link
CN (1) CN114357807B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114771877A (zh) * 2022-05-26 2022-07-22 哈尔滨工业大学 一种考虑导航误差的最优拦截制导方法
CN115092421A (zh) * 2022-06-22 2022-09-23 哈尔滨工业大学 一种基于轨道预报和Lambert变轨方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5804812A (en) * 1996-10-29 1998-09-08 Mcdonnell Douglas Corporation Multiple node lambert guidance system
CN105716470A (zh) * 2016-03-22 2016-06-29 北京航空航天大学 一种微分对策反拦截机动突防/精确打击导引方法
CN106643298A (zh) * 2016-11-29 2017-05-10 北京宇航系统工程研究所 一种基于预置命中点的大气层内反导拦截器中制导方法
CN108363299A (zh) * 2018-01-26 2018-08-03 北京航空航天大学 一种外大气层拦截最优末制导方法
CN109506517A (zh) * 2018-11-21 2019-03-22 中国人民解放军空军工程大学 一种带约束的中制导弹道优化方法
CN109857145A (zh) * 2018-11-27 2019-06-07 北京航空航天大学 一种基于迭代预测命中点的增程型拦截弹预测制导方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5804812A (en) * 1996-10-29 1998-09-08 Mcdonnell Douglas Corporation Multiple node lambert guidance system
CN105716470A (zh) * 2016-03-22 2016-06-29 北京航空航天大学 一种微分对策反拦截机动突防/精确打击导引方法
CN106643298A (zh) * 2016-11-29 2017-05-10 北京宇航系统工程研究所 一种基于预置命中点的大气层内反导拦截器中制导方法
CN108363299A (zh) * 2018-01-26 2018-08-03 北京航空航天大学 一种外大气层拦截最优末制导方法
CN109506517A (zh) * 2018-11-21 2019-03-22 中国人民解放军空军工程大学 一种带约束的中制导弹道优化方法
CN109857145A (zh) * 2018-11-27 2019-06-07 北京航空航天大学 一种基于迭代预测命中点的增程型拦截弹预测制导方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
SHILEI ZHAO 等: "Optimal Guidance Law with Impact-Angle Constraint and Acceleration Limit for Exo-Atmospheric Interception", 《AEROSPACE》, vol. 8, no. 12, 23 November 2021 (2021-11-23), pages 1 - 28 *
呼卫军 等: "J2 项摄动下的远程拦截耗尽关机中制导律设计", 《宇航学报》, vol. 38, no. 7, 30 July 2017 (2017-07-30), pages 694 - 703 *
杜润乐 等: "大气层外拦截器最优中制导律设计", 《航天控制》, vol. 29, no. 4, 15 August 2011 (2011-08-15), pages 27 - 31 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114771877A (zh) * 2022-05-26 2022-07-22 哈尔滨工业大学 一种考虑导航误差的最优拦截制导方法
CN114771877B (zh) * 2022-05-26 2022-11-18 哈尔滨工业大学 一种考虑导航误差的最优拦截制导方法
CN115092421A (zh) * 2022-06-22 2022-09-23 哈尔滨工业大学 一种基于轨道预报和Lambert变轨方法
CN115092421B (zh) * 2022-06-22 2023-04-11 哈尔滨工业大学 一种基于轨道预报和Lambert变轨方法

Also Published As

Publication number Publication date
CN114357807B (zh) 2022-05-27

Similar Documents

Publication Publication Date Title
CN114357807B (zh) 一种用于大气层外拦截的最优中制导方法及装置
Theodoulis et al. Guidance and control design for a class of spin-stabilized fin-controlled projectiles
Taub et al. Intercept angle missile guidance under time varying acceleration bounds
Liu et al. Solving the maximum-crossrange problem via successive second-order cone programming with a line search
Cho et al. Optimal impact angle control guidance law based on linearization about collision triangle
Li et al. A dimension-reduction solution of free-time differential games for spacecraft pursuit-evasion
de Celis et al. Guidance and control for high dynamic rotating artillery rockets
Fresconi et al. Practical assessment of real-time impact point estimators for smart weapons
Zhu et al. Impact time and angle control guidance independent of time-to-go prediction
Xiong et al. Hyperbolic tangent function weighted optimal intercept angle guidance law
CN109446582B (zh) 一种考虑地球自转的高精度降阶平稳滑翔动力学建模方法
Reisner et al. Optimal guidance-to-collision law for an accelerating exoatmospheric interceptor missile
Fresconi et al. High maneuverability projectile flight using low cost components
CN108073742B (zh) 基于改进粒子滤波算法的拦截导弹末段飞行状态估计方法
Hu et al. Analytical solution for nonlinear three-dimensional guidance with impact angle and field-of-view constraints
Zhang et al. Optimal guidance law design for impact with terminal angle of attack constraint
Ollerenshaw et al. Model predictive control of a direct fire projectile equipped with canards
Celis et al. GNSS/IMU laser quadrant detector hybridization techniques for artillery rocket guidance
Zhou et al. Generalized quasi-spectral model predictive static programming method using Gaussian quadrature collocation
Lukacs et al. Trajectory-shape-varying missile guidance for interception of ballistic missiles during the boost phase
Gaudet et al. Integrated and Adaptive Guidance and Control for Endoatmospheric Missiles via Reinforcement Meta-Learning
Fresconi et al. Theory, guidance, and flight control for high maneuverability projectiles
Głębocki et al. Sensitivity analysis and flight tests results for a vertical cold launch missile system
Fresconi et al. Flight performance of a small diameter munition with a rotating wing actuator
Doroshin Evolution of the precessional motion of unbalanced gyrostats of variable structure

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