CN109933847A - 一种改进的主动段弹道估计算法 - Google Patents
一种改进的主动段弹道估计算法 Download PDFInfo
- Publication number
- CN109933847A CN109933847A CN201910091116.3A CN201910091116A CN109933847A CN 109933847 A CN109933847 A CN 109933847A CN 201910091116 A CN201910091116 A CN 201910091116A CN 109933847 A CN109933847 A CN 109933847A
- Authority
- CN
- China
- Prior art keywords
- thrust
- vector
- missile target
- angle
- 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
Links
- 238000004422 calculation algorithm Methods 0.000 title claims abstract description 34
- 239000013598 vector Substances 0.000 claims abstract description 132
- 238000000034 method Methods 0.000 claims abstract description 79
- 239000011159 matrix material Substances 0.000 claims abstract description 72
- 238000001514 detection method Methods 0.000 claims abstract description 45
- 230000001133 acceleration Effects 0.000 claims abstract description 43
- 230000008569 process Effects 0.000 claims abstract description 40
- 238000001914 filtration Methods 0.000 claims abstract description 14
- 230000003044 adaptive effect Effects 0.000 claims abstract description 8
- 239000000446 fuel Substances 0.000 claims description 30
- 238000005070 sampling Methods 0.000 claims description 13
- 230000008859 change Effects 0.000 claims description 10
- 238000005259 measurement Methods 0.000 claims description 9
- 230000009466 transformation Effects 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 6
- 230000003190 augmentative effect Effects 0.000 claims description 5
- 239000000523 sample Substances 0.000 claims description 4
- 230000003416 augmentation Effects 0.000 abstract description 6
- 238000006243 chemical reaction Methods 0.000 abstract 1
- 230000005648 markovian process Effects 0.000 abstract 1
- 238000013179 statistical model Methods 0.000 description 7
- 230000005653 Brownian motion process Effects 0.000 description 6
- 238000004088 simulation Methods 0.000 description 6
- 230000009286 beneficial effect Effects 0.000 description 4
- 238000004364 calculation method Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 238000009795 derivation Methods 0.000 description 2
- 239000002904 solvent Substances 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 238000012795 verification Methods 0.000 description 2
- 238000005311 autocorrelation function Methods 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 238000005315 distribution function Methods 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000008447 perception Effects 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Landscapes
- Aiming, Guidance, Guns With A Light Source, Armor, Camouflage, And Targets (AREA)
Abstract
本发明提供一种改进的主动段弹道估计算法。该方法包括：根据导弹目标的初始质量、各级发动机的喷气速度和关机时刻等计算推力加速度模板；将推力方向角增广为状态矢量分量，利用增广后的状态矢量和推力加速度模板建立精细的参数化动力学模型；根据导弹目标和两颗预警卫星在地球固联坐标系中的位置矢量及地球固联坐标系到每颗预警卫星轨道坐标系的坐标转换矩阵，建立预警卫星探测模型；根据前3组等效探测位置，估算导弹目标初始状态矢量和初始状态协方差矩阵；利用一阶Markov过程描述推力方向角，构造自适应推力方向角变化的过程噪声矩阵；利用UKF滤波算法进行滤波处理，估计导弹目标的运动状态。本发明可提高预警卫星探测弹道导弹弹道估计的精度。
Description
技术领域
本发明涉及空间目标信号感知与处理技术领域，尤其涉及一种改进的主动段弹道估计方法。
背景技术
面对日益增大的弹道导弹威胁，各国大力发展天基红外弹道导弹预警系统，其典型代表是美国的天基红外系统(Space Based Infrared System，SBIRS)，它通过探测弹道导弹主动段发动机辐射红外信号来估计弹道参数，包括：导弹的发射点、射向、关机时间、关机点位置速度和落点位置等众多信息，实质上属于机动目标跟踪问题。天基红外探测属于被动探测，估计目标运动状态的前提是对目标运动进行建模。由于弹道导弹主动段受到推力、气动力和地心引力等多种力的影响，受力复杂，导弹质量实时变化，且有可能进行机动，所以，运动建模比较困难是天基红外系统探测弹道导弹主动段弹道估计面临的主要问题之一。
目前针对主动段弹道导弹的运动建模主要有两类方法：基于模板建模和基于模型建模。基于模型建模不使用任何先验信息，鲁棒性好，但估计精度较低。基于模板建模通过利用导弹的先验信息，提高估计精度，目前基于模板建模大多依赖标准弹道模板，估计精度由标准弹道模板库的准确性和完备程度决定，难以实际应用。针对基于标准弹道模板建模实用性不高的问题，有的利用推力加速度模板建模，但假设速度矢量与推力矢量方向相同符合单级近程弹道导弹主动段的运动，但不适用于多级远程弹道导弹在真空中飞行的情况。
发明内容
针对如何通过引入有限先验信息来提高弹道导弹主动段弹道估计精度的问题，本发明提供一种改进的主动段弹道估计算法，通过引入导弹的部分特征参数，可提高预警卫星探测弹道估计的精度。
本发明提供一种改进的主动段弹道估计算法，包括以下步骤：
步骤1、根据导弹目标的初始总质量、载荷质量、各级发动机的喷气速度、舱体质量、燃料质量、燃料秒耗量和关机时刻计算导弹目标的推力加速度模板；
步骤2、将导弹目标的推力方向角增广为导弹目标的状态矢量分量，利用增广后的导弹目标的状态矢量和所述推力加速度模板建立导弹目标的主动段参数化动力学模型；
步骤3、根据导弹目标在地球固联坐标系中的位置矢量、两颗预警卫星在地球固联坐标系中的位置矢量和地球固联坐标系分别到两颗预警卫星轨道坐标系的坐标转换矩阵，建立天基红外系统的预警卫星探测模型；
步骤4、根据所述两颗预警卫星探测到的导弹目标的前3组等效探测位置，计算导弹目标的初始状态矢量和初始状态协方差矩阵；
步骤5、利用一阶Markov过程描述推力方向角参数，并根据所述推力方向角参数构造自适应推力方向角变化的过程噪声矩阵；
步骤6、根据所述初始状态矢量和初始状态协方差矩阵，利用设定的UKF 滤波算法基于所述主动段运动模型和所述预警卫星探测模型进行滤波递推处理，估计导弹目标的运动状态。
进一步地，所述步骤1具体包括：
步骤1.1、根据第i级发动机燃料的质量mfi、第i级发动机舱体的质量mci和载荷的质量mload计算导弹目标的初始总质量M为发动机级数；
步骤1.2、根据第i级发动机燃料的质量mfi、第i-1级发动机和第i级发动机的关机时刻和计算第i级发动机燃料的秒耗量
步骤1.3、根据第i级发动机的喷气速度uei和第i级发动机燃料的秒耗量计算第i级发动机推力Fi；
步骤1.4、根据导弹目标的初始总质量和每一级发动机燃料的质量、发动机舱体的质量、发动机的关机时刻、发动机燃料的秒耗量及发动机推力计算导弹目标的推力加速度。
进一步地，所述步骤2具体为：
步骤2.1、根据式(5)计算导弹目标的总加速度模板：
其中，和γ分别为导弹目标的推力在地球固联坐标系中的俯仰角、偏航角和滚转角，ae为地球的长半轴，ωe为地球自转角速度，J2为地球非球形摄动，μ为地球引力常数，ap为推力加速度ap的模，[x y z vx vy vz]T为导弹目标的六维状态矢量，[x y z]T为导弹目标在地球固联坐标系中的位置矢量，[vx vy vz]T为导弹目标在地球固联坐标系中的速度矢量；
步骤2.2、将所述俯仰角和偏航角ψ增广为所述六维状态矢量的分量，建立状态矢量微分方程(6)：
步骤2.3、定义导弹目标在k时刻的状态矢量对所述状态矢量微分方程(6)进行离散化处理，得到导弹目标的主动段参数化动力学模型(7)：
Xk+1＝ΦkXk+Ukak (7)
其中，k表示离散时间点，T为采样周期。
进一步地，所述步骤3具体为：
步骤3.1、按照式(9)计算导弹目标在两颗预警卫星轨道坐标系中的位置矢量：
其中，为第l个预警卫星在地球固联坐标系中的位置矢量，[x y z]T为导弹目标在地球固联坐标系中的位置矢量，为地球固联坐标系到第l个预警卫星轨道坐标系的坐标转换矩阵，l＝1,2；
步骤3.2、根据所述导弹目标在两颗预警卫星轨道坐标系中的位置矢量按照式(10)计算两组预警卫星的方位角Al和俯仰角El：
其中，方位角Al为位置矢量在平面内的投影与轴的夹角，俯仰角El为位置矢量与平面的夹角；
步骤3.3、将得到的两组所述预警卫星的方位角Al和俯仰角El转换为导弹目标的等效探测位置。
进一步地，所述步骤4具体为：
步骤4.1、定义探测方程(11)：
YS＝hS·θ2+νS (11)
其中，为k＝0,1,2时刻的等效探测噪声矢量，YS为k＝0,1,2时刻的等效探测位置矢量组成的9维矢量，为待估计的k＝2时刻的状态矢量，r2为位置矢量，速度矢量，为加速度矢量，为测量矩阵， T为采样周期，0i×j为i行j列的零矩阵，I为单位矩阵；
步骤4.2、按照式(12)计算等效探测噪声矢量vS的等效噪声协方差矩阵：
其中，[Pr0,Pr1,Pr2]为k＝0,1,2时刻的位置矢量的协方差矩阵；
步骤4.3、利用加权最小二乘法按照式(13)计算θ2的估计
步骤4.4、按照式(14)计算的协方差矩阵
步骤4.5、将k＝2时刻的位置矢量[x2 y2 z2]T作为初始位置矢量并按照式 (15)计算导弹目标的推力的初始俯仰角和初始偏航角ψ2：
步骤4.6、定义导弹目标的初始状态矢量估计值并按照式(16)计算初始状态协方差矩阵
其中，为矩阵前6行、前6列组成的子矩阵，为导弹目标的推力的初始俯仰角的方差，为导弹目标的推力的初始偏航角的方差，0i×j为 i行j列的零矩阵。
进一步地，所述步骤5具体为：
步骤5.1、定义导弹目标的推力的俯仰角和偏航角ψ的参数模型(17)：
步骤5.2、对所述俯仰角和偏航角ψ的参数模型(17)按照式(18)进行离散化处理，得到：
步骤5.3、按照式(19)计算推力的俯仰角过程噪声和推力的偏航角过程噪声Qψ：
其中，和τψ分别为俯仰角和偏航角ψ的机动时间常数，和wψk分别为用于辨识俯仰角和偏航角ψ的过程白噪声序列，T为采样周期；
步骤5.4、根据所述推力的俯仰角过程噪声和推力的偏航角过程噪声Qψ按照式(20)计算自适应推力方向角变化的过程噪声矩阵Qk：
其中，0i×j为i行j列的零矩阵。
进一步地，所述步骤6具体为：
步骤6.1、构造2n+1个Sigma点χk(i)和权重Wi，i＝0,1,2,…,2n；
步骤6.2、χk(i)基于状态方程的非线性传播，得到
步骤6.3、计算一步状态预测的均值自适应推力方向角变化的过程噪声矩阵Qk和协方差Pk+1|k；
步骤6.4、基于所述均值重新构造Sigma点χk+1|k(i)和权重Wi'；
步骤6.5、基于探测方程的非线性传播求出
步骤6.6、计算探测值的预估值自协方差矩阵和互协方差矩阵
步骤6.7、计算基于探测量的状态更新与协方差更新Pk+1|k+1。
本发明的有益效果：
本发明实施例提供的一种改进的主动段弹道估计算法，具有以下有益效果：
(1)具有较强实用性。建模采取利用不完备的加速度模板，与传统的基于标准弹道模板建模相比，有更强的实用性。
(2)弹道估计精度高。与传统的基于模型建模相比，利用了一些导弹的先验信息，提高了弹道估计的精度。
(3)状态增广实现联合估计。将导弹推力方向角增广为目标状态分量，在估计导弹状态的同时，可以获得推力方向角较为准确的估计。
(4)对推力姿态角初始化有理论依据。利用导弹20秒转弯不大的特性，用估计的初始位置地心矢径，对推力方向角进行初始化，比较合理，有效减小了初始误差。
(5)过程噪声可以随着状态的变化进行自适应调整。采用一阶Markov过程描述推力方向角，较Wiener过程更符合实际，获得的过程噪声可以随着目标状态变化进行自适应的调整。
(6)可信度强。通过正确的目标运动模型建模、严谨的数学推导、合理的工程近似手段和充分的仿真验证试验，提出基于不完备加速度模板的弹道估计方法，具有较强的可信度。
附图说明
图1为本发明实施例提供的一种改进的主动段弹道估计算法的流程示意图。
图2为本发明实施例提供的导弹与卫星的几何关系示意图；
图3为本发明实施例提供的UKF滤波处理流程示意图；
图4为本发明实施例提供的推力加速度变化曲线图；
图5为本发明实施例提供的主动段弹道三维曲线图；
图6为本发明实施例提供的x方向位置估计的均方根误差-时间曲线图；
图7为本发明实施例提供的x方向速度估计的均方根误差-时间曲线图；
图8为本发明实施例提供的推力俯仰角和偏航角估计值与真实值对比示意图。
具体实施方式
为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
图1为本发明实施例提供的一种改进的主动段弹道估计算法的流程示意图。该算法包括以下步骤：
S101、根据导弹目标的初始总质量、载荷质量、各级发动机的喷气速度、舱体质量、燃料质量、燃料秒耗量和关机时刻计算导弹目标的推力加速度模板；
S102、将导弹目标的推力方向角增广为导弹目标的状态矢量分量，利用增广后的导弹目标的状态矢量和所述推力加速度模板建立导弹目标的主动段参数化动力学模型；
S103、根据导弹目标在地球固联坐标系中的位置矢量、两颗预警卫星在地球固联坐标系中的位置矢量和地球固联坐标系分别到两颗预警卫星轨道坐标系的坐标转换矩阵，建立天基红外系统的预警卫星探测模型；
S104、根据两颗所述预警卫星探测到的导弹目标的前3组等效探测位置，计算导弹目标的初始状态矢量和初始状态协方差矩阵；
具体地，所述初始状态矢量包括初始位置矢量、初始速度矢量和初始推力方向角；
S105、利用一阶Markov过程描述推力方向角参数，并根据所述推力方向角参数构造自适应推力方向角变化的过程噪声矩阵；
S106、根据所述初始状态矢量和初始状态协方差矩阵，利用设定的UKF 滤波算法基于所述主动段运动模型和所述预警卫星探测模型进行滤波递推处理，估计导弹目标的运动状态。
本发明实施例提供的一种改进的主动段弹道估算方法，具有以下有益效果：
(1)本发明实施例通过引入导弹目标的一些先验信息，求出不完备推力加速度模板，然后将推力方向角增广为目标状态分量，并基于推力加速度模板信息建立精细的参数化动力学模型，相比于现有技术中的基于标准弹道模板建模，本发明实施例具有较强实用性且有效提高了弹道估计的精度和稳定性。同时将导弹推力方向角增广为目标状态分量，在估计导弹状态的同时，可以获得推力方向角较为准确的估计。
(2)传统状态增广法用Wiener过程来描述增广的参数分量，而推力方向角具有非随机性和分段相关性，用Wiener过程描述推力方向角与实际差别较大，针对该问题，本发明实施例采用一阶Markov过程描述推力方向角，实现了过程噪声自适应调整。
在上述实施例的基础上，本发明提供有一种实施例，其流程包括以下步骤：
S201、计算导弹目标的推力加速度模板。对同一种类型的弹道导弹，推力加速度的大小有固定的变化规律，可以利用弹道导弹发动机的喷气速度、燃料质量、秒耗量等特征参数，求出只有大小没有方向的不完备推力加速度模板，具体地，包括以下子步骤：
S2011、弹道导弹通常分为发动机和有效载荷两个部分，发动机由燃料和装燃料的舱体组成。因此，可根据第i级发动机燃料的质量mfi、第i级发动机舱体的质量mci和载荷的质量mload按照式(1)计算导弹目标的初始总质量 M为发动机级数：
例如，3级弹道导弹的初始总质量
S2012、根据第i级发动机燃料的质量mfi、第i-1级发动机和第i级发动机的关机时刻和计算第i级发动机燃料的秒耗量
具体地，在实际应用中，可假设各级发动机燃料的秒耗率均为常数(即根据式(2)可计算得到各级发动机燃料的耗率，i＝1,2,…,M，M为发动机级数：
例如，以3级弹道导弹为例，各级发动机的关机时刻依次为和秒，则：
S2013、根据第i级发动机的喷气速度uei和第i级发动机燃料的秒耗量计算第i级发动机推力Fi；
具体地，天基红外系统探测弹道导弹的大部分时间内，导弹在稀薄大气或真空中飞行，可按照式(3)计算得到各级发动机推力，i＝1,2,…,M，M为发动机级数：
其中，g为重力加速度常数。
例如，以三级弹道导弹为例，各级发动机的喷气速度依次为ue1、ue2、ue3，则各级发动机推力可以近似为：
S2014、根据导弹目标的初始总质量和每一级发动机燃料的质量、发动机舱体的质量、发动机的关机时刻、发动机燃料的秒耗量及发动机推力计算导弹目标的推力加速度。
具体地，作为一种可实施方式，可按照式(4)计算计算推力加速度ap：
例如，以三级弹道导弹为例，t时刻推力加速度ap(t)
S202、将推力方向角增广为状态分量，建立主动段参数化动力学模型。具体包括以下子步骤：
S2021、根据式(5)计算导弹目标的总加速度模板；
具体地，由于天外红外系统探测弹道导弹的大部分时间，气动力比较小可以忽略。因此在地球固联坐标系中进行描述，主动段弹道导弹加速度a可以简化为：
a＝ap+ag+ae+ac
其中，ap为推力加速度、ag为地心引力加速度、ae为离心力加速度，ac为科氏力加速度。因此，导弹目标的总加速度：
其中，和γ分别为导弹目标的推力在地球固联坐标系中的俯仰角、偏航角和滚转角，ae为地球的长半轴，ωe为地球自转角速度，J2为地球非球形摄动，μ为地球引力常数，ap为推力加速度ap的模，[x y z vx vy vz]T为导弹目标的六维状态矢量，[x y z]T为导弹目标在地球固联坐标系中的位置矢量，[vx vy vz]T为导弹目标在地球固联坐标系中的速度矢量。从(5)式可以看出，总加速度仅含有俯仰角偏航角ψ两个未知量，与滚转角γ无关。
S2022、将所述俯仰角和偏航角ψ增广为所述六维状态矢量的分量，建立状态矢量微分方程(6)：
S2023、定义导弹目标在k时刻的状态矢量对所述状态矢量微分方程(6)进行离散化处理，得到导弹目标的主动段参数化动力学模型(7)：
Xk+1＝ΦkXk+Ukak (7)
其中，k表示离散时间点，T为采样周期。
S203、描述天基红外系统的探测模型。若天基红外系统中的L个预警卫星同时探测到弹道导弹，探测方程为
Zk+1(A1,E1,A2,E2)＝h(Xk+1) (8)
其中，(A1,E1,A2,E2)为k+1时刻两颗卫星的观测值，也称角探测量。Xk+1为k+1时刻的状态矢量；该探测方程为描述k+1时刻两颗卫星的观测值与k+1时刻的状态矢量关系的隐式方程。
具体包括以下子步骤：
S2031、按照式(9)计算导弹目标在两颗预警卫星轨道坐标系中的位置矢量：
其中，为第l个预警卫星在地球固联坐标系中的位置矢量，[x y z]T为导弹目标在地球固联坐标系中的位置矢量，为地球固联坐标系到第l个预警卫星轨道坐标系的坐标转换矩阵，l＝1,2；
S2032、根据所述导弹目标在两颗预警卫星轨道坐标系中的位置矢量按照式 (10)计算两组预警卫星的方位角Al和俯仰角El：
其中，方位角Al为位置矢量在平面内的投影与轴的夹角，俯仰角El为位置矢量与平面的夹角。
例如，作为一种可实施方式，若天基红外系统中的两颗预警卫星同时探测到弹道导弹，则探测方程为Zk+1(A1,E1,A2,E2)＝h(Xk+1)。为了简化表述，以下的公式均省略下标k+1。
假设导弹目标在地球固联坐标系中的位置矢量为[x y z]T，两颗卫星在地球固联坐标系中的位置矢量分别为和地球固联坐标系到每个卫星轨道坐标系的坐标转换矩阵分别为和则导弹目标在两个卫星轨道坐标系中的位置矢量分别为
由导弹目标的位置矢量可以计算出预警卫星红外探测器的两个角探测量：
S2033、将得到的两组所述预警卫星的方位角Al和俯仰角El转换为导弹目标的等效探测位置；
具体地，若有两颗卫星可以同时对目标进行观测，则利用双星观测数据可以实现目标的立体定位。如图2所示，理论情况下两颗卫星对导弹的视线矢量应该在由卫星和导弹所决定的平面内，但由于观测数据存在误差，视线矢量与该平面存在略微角度，不过在粗定位情况下，仍可认为视线矢量与该平面在同一平面内，由三角形ΔS1S2P内的三角几何关系，根据双星观测数据可以计算目标的位置。即获得X＝h-1(Z)的表达式。
已知卫星S1和S2在地球固联坐标系中的位置矢量分别为和在三角形ΔS1S2P中，定义三个角分别为S1、S2和P，所对的边分别为n、m和矢量c在卫星S1轨道坐标系中可表示为
而c的长度为
c＝||b-a||
在视线上任取一点P10(不妨取)，则矢量e在卫星S1轨道坐标系中描述为：
则角S1可由e和c计算：
类似的可得：
在三角形ΔS1S2P中有：
式中m、n和c分别为矢量m、n和c的模，整理得：
利用一元二次方程的求根公式可以计算矢量m的模m。矢量m在卫星S1轨道坐标系中可表示为
将其转换到地球固联坐标系中，可得：
交点P在地球固联坐标系中的坐标矢量dECF为：
dECF＝mECF+aECF
上式计算的位置坐标dECF即为根据双星观测数据计算的等效探测位置，其协方差矩阵为：
其中，偏导数矩阵可采用数值方法求解。
S204、滤波算法初始化，求初始状态矢量和初始状态协方差矩阵。具体地，若采用匀加速模型求导弹目标的初始状态，需获得k＝0～2的三组等效位置探测数据和对应的协方差矩阵。具体包括以下子步骤：
S2041、设待估矢量θ2为k＝2时的状态矢量，定义探测方程(11)：
YS＝hS·θ2+νS (11)
其中，为k＝0,1,2时刻的等效探测噪声矢量，YS为k＝0,1,2时刻的等效探测位置矢量组成的9维矢量，为待估计的k＝2时刻的状态矢量，r2为位置矢量，速度矢量，为加速度矢量，为测量矩阵， T为采样周期，0i×j为i行j列的零矩阵，I为单位矩阵，
具体地，已知等效探测位置矢量YS，采样周期T，测量矩阵hS以及等效探测噪声矢量vS，根据可得待估计的k＝2时刻的状态矢量θ2。
S2042、按照式(12)计算等效探测噪声矢量vS的等效噪声协方差矩阵RS：
其中，[Pr0,Pr1,Pr2]为k＝0,1,2时刻的位置矢量的协方差矩阵。
S2043、利用加权最小二乘法按照式(13)计算θ2的估计
S2044、按照式(14)计算的协方差矩阵
S2045、由于天基红外系统首次探测到弹道导弹，约在导弹发射20s左右，导弹转弯的角度不大，可以用初始估计位置地心矢径的方向近似推力的方向。即将k＝2时刻的位置矢量[x2 y2 z2]T作为初始位置矢量并按照式(15)计算导弹目标的推力的初始俯仰角和初始偏航角ψ2：
S2046、定义导弹目标的初始状态矢量估计值由于俯仰角、偏航角与位置、速度分量没有关系，且俯仰角与偏航角也不相关。因此，按照式(16)计算初始状态协方差矩阵其中，为矩阵前6行、前6列组成的子矩阵，为导弹目标的推力的初始俯仰角的方差，为导弹目标的推力的初始偏航角的方差，0i×j为 i行j列的零矩阵。
S205、构造自适应的过程噪声。传统的状态增广法一般用Wiener过程来描述增广的参数分量，Wiener过程表示的参数序列互不相关，但导弹推力方向角具有非随机性和分段相关性，因此用自相关函数为指数衰减形式的一阶Markov 过程描述推力方向角参数更符合实际。具体包括以下子步骤：
S2051、定义导弹目标的推力的俯仰角和偏航角ψ的参数模型(17)：
S2052、对所述俯仰角和偏航角ψ的参数模型(17)按照式(18)进行离散化处理，得到：
S2053、按照式(19)计算推力的俯仰角过程噪声和推力的偏航角过程噪声Qψ：
其中，和τψ分别为俯仰角和偏航角ψ的机动时间常数，和wψk分别为用于辨识俯仰角和偏航角ψ的过程白噪声序列，T为采样周期；
S2054、根据所述推力的俯仰角过程噪声和推力的偏航角过程噪声Qψ按照式(20)计算自适应推力方向角变化的过程噪声矩阵Qk：
其中，0i×j为i行j列的零矩阵。和分别反映了俯仰角和偏航角的变化范围，二者均代表了目标的机动能力；和τψ则反映了目标的实时机动策略，不同的和τψ对应于不同的实时运动模式。
从简化后的式(19)可以看出，推力的俯仰角过程噪声推力的偏航角过程噪声Qψ与俯仰角方差和采样周期T成正比，与机动时间常数成反比，物理意义较为明确，且能够根据推力俯仰角的变化情况进行自适应调整。
S206、设计用于弹道估计的滤波算法，进行滤波处理。
目前，用于弹道估计的滤波算法主要有扩展卡尔曼滤波(Extended KalmanFilter，EKF)和无迹卡尔曼滤波(Unscented Kalman Filter，UKF)。天基红外系统探测模型和弹道导弹参数化动力学模型都是非线性的，尤其是将两个推力方向角增广为状态分量后，导弹运动模型的非线性进一步增强。EKF是用Taylor展开的一次项近似非线性系统，忽略高阶分量，用于弹道估计可能会导致较大的截断误差。UKF利用无迹变换(UnscentedTransform，UT)构建一组Sigma取样点，用取样点来近似后验概率分布函数；而后利用非线性状态方程和探测方程，获得状态矢量经非线性传播后的统计特性，可以获得精确到二阶的后验均值和协方差。具体包括以下子步骤：
在获得k时刻的状态估计值以及相应的协方差矩阵Pk|k后，一个完整的运算周期图3所示，主要分为以下七步：
S2061、构造2n+1个Sigma点χk(i)和权重Wi；
具体地，可按照式(21)进行构造。
其中，为对协方差矩阵Pk进行Cholesky分解获得的上三角阵的第i行， n为状态矢量维数，n+κ＝3，若n>3则标量参数κ取相应的负值。
S2062、χk(i)基于状态方程(即公式(7))的非线性传播，得到
具体地，
S2063、计算一步状态预测的均值自适应推力方向角变化的过程噪声矩阵Qk和协方差Pk+1|k；
具体地，和协方差Pk+1|k可按照式(23)进行计算：
过程噪声矩阵Qk按照式(20)进行计算。
S2064、基于所述均值重新构造Sigma点χk+1|k(i)和权重Wi'；
具体地，可按照式(24)进行构造。
S2065、基于探测方程(即公式(8))的非线性传播求出
S2066、计算探测值的预估值自协方差矩阵和互协方差矩阵
具体地，可按照式(26)进行计算：
S2067、计算基于角探测量的状态更新与协方差更新Pk+1|k+1。
本发明实施例提供的一种改进的主动段弹道估计算法，具有以下有益效果：
(1)具有较强实用性。建模采取利用不完备的加速度模板，与传统的基于标准弹道模板建模相比，有更强的实用性。
(2)弹道估计精度高。与传统的基于模型建模相比，利用了一些导弹的先验信息，提高了弹道估计的精度。
(3)状态增广实现联合估计。将导弹推力方向角增广为目标状态分量，在估计导弹状态的同时，可以获得推力方向角较为准确的估计。
(4)对推力姿态角初始化有理论依据。利用导弹20秒转弯不大的特性，用估计的初始位置地心矢径，对推力方向角进行初始化，比较合理，有效减小了初始误差。
(5)过程噪声可以随着状态的变化进行自适应调整。采用一阶Markov过程描述推力方向角，较Wiener过程更符合实际，获得的过程噪声可以随着目标状态变化进行自适应的调整。
(6)可信度强。通过正确的目标运动模型建模、严谨的数学推导、合理的工程近似手段和充分的仿真验证试验，提出基于不完备加速度模板的弹道估计方法，具有较强的可信度。
下面对本发明提供的一种改进的主动段弹道估计算法进行仿真实验，进一步验证本发明算法的有效性。
以某三级导弹为例，假设推力加速度变化曲线如图4所示。导弹主动段飞行曲线如图5所示。分别部署在东经208°和东经110°的两颗SBIRS-GEO发现导弹目标，首次探测到导弹的时间为20s，探测间隔为1s，最后探测到导弹的时间为 196s，视线测量误差取50μrad。对本发明算法和基于当前统计模型CS的传统算法蒙特卡洛仿真500次。仿真参数设置如下：
基于当前统计模型CS的传统算法：机动频率为0.001Hz，x轴、y轴、z轴方向加速度最大值均为200m/s2。
本发明算法：推力俯仰角机动时间常数为200s，推力偏航角机动时间常数τψ为200s。由于导弹在主动段缓慢转弯，发射20秒后导弹纵轴与地心矢量有一定偏差，根据经验，导弹发射20s时推力俯仰角标准差设为50°，偏航角标准差σψ0设为50°。
以x轴方向状态估计误差为例，其它方向基本相同，仿真结果如图6和图 7所示。实线描述的是本发明算法估计结果，虚线描述的是基于当前统计模型 CS的弹道估计结果。从图6和图7可以看出，本发明算法估计精度更高。位置估计的均方根误差比基于当前统计模型CS约小25％；速度估计的均方根误差约比基于当前统计模型CS小50％。另外，从图7可以看出，本发明算法收敛速度更快，同样收敛到50m/s，本发明算法仅需要17s左右，基于当前统计模型CS 的传统算法则需要24s。
其次，当前统计模型CS在68s和145s后速度估计的均方根误差突然增大，出现尖峰，是因为导弹发动机在两个时刻进行分离。而本发明算法在导弹发动机分级时，估计误差没有明显变化，原因有两点：因为本发明算法已知推力加速度在这两个时刻的变化情况，建立的运动模型较为精确；另外，本发明算法用一阶 Markov过程描述推力方向角过程噪声，使过程噪声能够根据推力方向角变化情况进行自适应调整，有效提高了滤波的稳定性。
最后，本发明算法可以获得推力方向角较为准确的估计。图8是本发明算法单次仿真的结果，可实时得到推力俯仰角和偏航角较为准确的估计值。使用初始估计位置的地心矢径对推力方向角进行初始化，俯仰角和偏航角与真实值分别偏差35°和12°，偏差比较小，且通过滤波，俯仰角和偏航角均快速收敛于真实值。
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。
Claims (7)
1.一种改进的主动段弹道估计算法，其特征在于，包括：
步骤1、根据导弹目标的初始总质量、载荷质量、各级发动机的喷气速度、舱体质量、燃料质量、燃料秒耗量和关机时刻计算导弹目标的推力加速度模板；
步骤2、将导弹目标的推力方向角增广为导弹目标的状态矢量分量，利用增广后的导弹目标的状态矢量和所述推力加速度模板建立导弹目标的主动段参数化动力学模型；
步骤3、根据导弹目标在地球固联坐标系中的位置矢量、两颗预警卫星在地球固联坐标系中的位置矢量和地球固联坐标系分别到两颗预警卫星轨道坐标系的坐标转换矩阵，建立天基红外系统的预警卫星探测模型；
步骤4、根据两颗所述预警卫星探测到的导弹目标的前3组等效探测位置，计算导弹目标的初始状态矢量和初始状态协方差矩阵；
步骤5、利用一阶Markov过程描述推力方向角参数，并根据所述推力方向角参数构造自适应推力方向角变化的过程噪声矩阵；
步骤6、根据所述初始状态矢量和初始状态协方差矩阵，利用设定的UKF滤波算法基于所述主动段运动模型和所述预警卫星探测模型进行滤波递推处理，估计导弹目标的运动状态。
2.根据权利要求1所述的方法，其特征在于，所述步骤1具体包括：
步骤1.1、根据第i级发动机燃料的质量mfi、第i级发动机舱体的质量mci和载荷的质量mload计算导弹目标的初始总质量M为发动机级数；
步骤1.2、根据第i级发动机燃料的质量mfi、第i-1级发动机和第i级发动机的关机时刻和计算第i级发动机燃料的秒耗量
步骤1.3、根据第i级发动机的喷气速度uei和第i级发动机燃料的秒耗量计算第i级发动机推力Fi；
步骤1.4、根据导弹目标的初始总质量和每一级发动机燃料的质量、发动机舱体的质量、发动机的关机时刻、发动机燃料的秒耗量及发动机推力计算导弹目标的推力加速度。
3.根据权利要求1所述的方法，其特征在于，所述步骤2具体为：
步骤2.1、根据式(5)计算导弹目标的总加速度模板：
其中，ψ和γ分别为导弹目标的推力在地球固联坐标系中的俯仰角、偏航角和滚转角，ae为地球的长半轴，ωe为地球自转角速度，J2为地球非球形摄动，μ为地球引力常数，ap为推力加速度ap的模，[x y z vx vy vz]T为导弹目标的六维状态矢量，[x y z]T为导弹目标在地球固联坐标系中的位置矢量，[vx vy vz]T为导弹目标在地球固联坐标系中的速度矢量；
步骤2.2、将所述俯仰角和偏航角ψ增广为所述六维状态矢量的分量，建立状态矢量微分方程(6)：
步骤2.3、定义导弹目标在k时刻的状态矢量对所述状态矢量微分方程(6)进行离散化处理，得到导弹目标的主动段参数化动力学模型(7)：
Xk+1＝ΦkXk+Ukak (7)
其中，k表示离散时间点，T为采样周期。
4.根据权利要求1所述的方法，其特征在于，所述步骤3具体为：
步骤3.1、按照式(9)计算导弹目标在两颗预警卫星轨道坐标系中的位置矢量：
其中，为第l个预警卫星在地球固联坐标系中的位置矢量，[x y z]T为导弹目标在地球固联坐标系中的位置矢量，为地球固联坐标系到第l个预警卫星轨道坐标系的坐标转换矩阵，l＝1,2；
步骤3.2、根据所述导弹目标在两颗预警卫星轨道坐标系中的位置矢量按照式(10)计算两组预警卫星的方位角Al和俯仰角El：
其中，方位角Al为位置矢量在平面内的投影与轴的夹角，俯仰角El为位置矢量与平面的夹角；
步骤3.3、将得到的两组所述预警卫星的方位角Al和俯仰角El转换为导弹目标的等效探测位置。
5.根据权利要求3所述的方法，其特征在于，所述步骤4具体为：
步骤4.1、定义探测方程(11)：
YS＝hS·θ2+νS (11)
其中，为k＝0,1,2时刻的等效探测噪声矢量，YS为k＝0,1,2时刻的等效探测位置矢量组成的9维矢量，为待估计的k＝2时刻的状态矢量，r2为位置矢量，速度矢量，为加速度矢量，为测量矩阵，T为采样周期，0i×j为i行j列的零矩阵，I为单位矩阵；
步骤4.2、按照式(12)计算等效探测噪声矢量vS的等效噪声协方差矩阵RS：
其中，[Pr0,Pr1,Pr2]为k＝0,1,2时刻的位置矢量的协方差矩阵；
步骤4.3、利用加权最小二乘法按照式(13)计算θ2的估计
步骤4.4、按照式(14)计算的协方差矩阵
步骤4.5、将k＝2时刻的位置矢量[x2 y2 z2]T作为初始位置矢量并按照式(15)计算导弹目标的推力的初始俯仰角和初始偏航角ψ2：
步骤4.6、定义导弹目标的初始状态矢量估计值并按照式(16)计算初始状态协方差矩阵
其中，为矩阵前6行、前6列组成的子矩阵，为导弹目标的推力的初始俯仰角的方差，为导弹目标的推力的初始偏航角的方差，0i×j为i行j列的零矩阵。
6.根据权利要求1所述的方法，其特征在于，所述步骤5具体为：
步骤5.1、定义导弹目标的推力的俯仰角和偏航角ψ的参数模型(17)：
步骤5.2、对所述俯仰角和偏航角ψ的参数模型(17)按照式(18)进行离散化处理，得到：
步骤5.3、按照式(19)计算推力的俯仰角过程噪声和推力的偏航角过程噪声Qψ：
其中，和τψ分别为俯仰角和偏航角ψ的机动时间常数，和wψk分别为用于辨识俯仰角和偏航角ψ的过程白噪声序列，T为采样周期；
步骤5.4、根据所述推力的俯仰角过程噪声和推力的偏航角过程噪声Qψ按照式(20)计算自适应推力方向角变化的过程噪声矩阵Qk：
其中，0i×j为i行j列的零矩阵。
7.根据权利要求1所述的方法，其特征在于，所述步骤6具体为：
步骤6.1、构造2n+1个Sigma点χk(i)和权重Wi，i＝0,1,2,…,2n；
步骤6.2、χk(i)基于状态方程的非线性传播，得到
步骤6.3、计算一步状态预测的均值自适应推力方向角变化的过程噪声矩阵Qk和协方差Pk+1|k；
步骤6.4、基于所述均值重新构造Sigma点χk+1|k(i)和权重Wi'；
步骤6.5、基于探测方程的非线性传播求出
步骤6.6、计算探测值的预估值自协方差矩阵和互协方差矩阵
步骤6.7、计算基于探测量的状态更新与协方差更新Pk+1|k+1。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910091116.3A CN109933847B (zh) | 2019-01-30 | 2019-01-30 | 一种改进的主动段弹道估计算法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910091116.3A CN109933847B (zh) | 2019-01-30 | 2019-01-30 | 一种改进的主动段弹道估计算法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109933847A true CN109933847A (zh) | 2019-06-25 |
| CN109933847B CN109933847B (zh) | 2022-09-16 |
Family
ID=66985357
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910091116.3A Expired - Fee Related CN109933847B (zh) | 2019-01-30 | 2019-01-30 | 一种改进的主动段弹道估计算法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109933847B (zh) |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110728026A (zh) * | 2019-09-16 | 2020-01-24 | 南京理工大学 | 一种基于角速度量测的末端弹道目标被动跟踪方法 |
| CN110807270A (zh) * | 2019-11-13 | 2020-02-18 | 北京环境特性研究所 | 一种基于尾焰辐射线型反演发动机参数及预估弹道的方法 |
| CN111462182A (zh) * | 2020-03-31 | 2020-07-28 | 南京航空航天大学 | 一种基于红外预警图像的弹道导弹三维轨迹估计方法 |
| CN111475767A (zh) * | 2020-03-18 | 2020-07-31 | 中国科学院紫金山天文台 | 考虑地球自转影响的最小能量弹道严格构造方法 |
| CN112269174A (zh) * | 2020-10-21 | 2021-01-26 | 中国人民解放军战略支援部队信息工程大学 | 基于交互式多模型信息融合的助推滑翔飞行器目标估计方法及系统 |
| CN112595319A (zh) * | 2020-11-18 | 2021-04-02 | 中国西安卫星测控中心 | 一种模型自适应补偿的返回弹道估计算法 |
| CN113962057A (zh) * | 2021-06-29 | 2022-01-21 | 南京航空航天大学 | 基于时序交会的远程导弹主动段运动参数修正方法 |
| CN116992553A (zh) * | 2023-05-25 | 2023-11-03 | 中国人民解放军32804部队 | 助推滑翔飞行器全程弹道估计方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20120316819A1 (en) * | 2008-04-22 | 2012-12-13 | United States Government, As Represented By The Secretary Of The Navy | Process for estimation of ballistic missile boost state |
| CN107421550A (zh) * | 2017-07-25 | 2017-12-01 | 北京航空航天大学 | 一种基于星间测距的地球‑Lagrange联合星座自主定轨方法 |
| CN107862279A (zh) * | 2017-11-03 | 2018-03-30 | 中国电子科技集团公司第三研究所 | 一种脉冲声信号识别分类方法 |
-
2019
- 2019-01-30 CN CN201910091116.3A patent/CN109933847B/zh not_active Expired - Fee Related
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20120316819A1 (en) * | 2008-04-22 | 2012-12-13 | United States Government, As Represented By The Secretary Of The Navy | Process for estimation of ballistic missile boost state |
| CN107421550A (zh) * | 2017-07-25 | 2017-12-01 | 北京航空航天大学 | 一种基于星间测距的地球‑Lagrange联合星座自主定轨方法 |
| CN107862279A (zh) * | 2017-11-03 | 2018-03-30 | 中国电子科技集团公司第三研究所 | 一种脉冲声信号识别分类方法 |
Non-Patent Citations (1)
| Title |
|---|
| 张峰等: "预警卫星对弹道导弹主动段状态估计算法", 《红外与激光工程》 * |
Cited By (14)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110728026B (zh) * | 2019-09-16 | 2022-10-21 | 南京理工大学 | 一种基于角速度量测的末端弹道目标被动跟踪方法 |
| CN110728026A (zh) * | 2019-09-16 | 2020-01-24 | 南京理工大学 | 一种基于角速度量测的末端弹道目标被动跟踪方法 |
| CN110807270A (zh) * | 2019-11-13 | 2020-02-18 | 北京环境特性研究所 | 一种基于尾焰辐射线型反演发动机参数及预估弹道的方法 |
| CN110807270B (zh) * | 2019-11-13 | 2023-09-29 | 北京环境特性研究所 | 一种基于尾焰辐射线型反演发动机参数及预估弹道的方法 |
| CN111475767A (zh) * | 2020-03-18 | 2020-07-31 | 中国科学院紫金山天文台 | 考虑地球自转影响的最小能量弹道严格构造方法 |
| CN111475767B (zh) * | 2020-03-18 | 2021-03-16 | 中国科学院紫金山天文台 | 考虑地球自转影响的最小能量弹道严格构造方法 |
| CN111462182A (zh) * | 2020-03-31 | 2020-07-28 | 南京航空航天大学 | 一种基于红外预警图像的弹道导弹三维轨迹估计方法 |
| CN112269174A (zh) * | 2020-10-21 | 2021-01-26 | 中国人民解放军战略支援部队信息工程大学 | 基于交互式多模型信息融合的助推滑翔飞行器目标估计方法及系统 |
| CN112595319A (zh) * | 2020-11-18 | 2021-04-02 | 中国西安卫星测控中心 | 一种模型自适应补偿的返回弹道估计算法 |
| CN112595319B (zh) * | 2020-11-18 | 2023-09-22 | 中国西安卫星测控中心 | 一种模型自适应补偿的返回弹道估计方法 |
| CN113962057A (zh) * | 2021-06-29 | 2022-01-21 | 南京航空航天大学 | 基于时序交会的远程导弹主动段运动参数修正方法 |
| CN113962057B (zh) * | 2021-06-29 | 2022-06-24 | 南京航空航天大学 | 基于时序交会的远程导弹主动段运动参数修正方法 |
| CN116992553A (zh) * | 2023-05-25 | 2023-11-03 | 中国人民解放军32804部队 | 助推滑翔飞行器全程弹道估计方法 |
| CN116992553B (zh) * | 2023-05-25 | 2024-02-06 | 中国人民解放军32804部队 | 助推滑翔飞行器全程弹道估计方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109933847B (zh) | 2022-09-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109933847B (zh) | 一种改进的主动段弹道估计算法 | |
| CN107544067B (zh) | 一种基于高斯混合近似的高超声速再入飞行器跟踪方法 | |
| CN103528587B (zh) | 自主组合导航系统 | |
| CN106885570A (zh) | 一种基于鲁棒sckf滤波的紧组合导航方法 | |
| CN106885569A (zh) | 一种强机动条件下的弹载深组合arckf滤波方法 | |
| CN103940433B (zh) | 一种基于改进的自适应平方根ukf算法的卫星姿态确定方法 | |
| CN111474960B (zh) | 基于控制量特征辅助的临空高超声速目标跟踪方法 | |
| CN107255924A (zh) | 基于扩维模型的容积卡尔曼滤波提取捷联导引头制导信息的方法 | |
| CN108073742B (zh) | 基于改进粒子滤波算法的拦截导弹末段飞行状态估计方法 | |
| CN104296753A (zh) | 一种基于多模型滤波的空间目标定位方法 | |
| Zhao et al. | A filtering approach based on MMAE for a SINS/CNS integrated navigation system | |
| CN112710298A (zh) | 基于动力学模型辅助的旋转弹用地磁卫星组合导航方法 | |
| CN109781374A (zh) | 一种实时在线快速估计飞行器推力的方法 | |
| CN110231619B (zh) | 基于恩克法的雷达交接时刻预报方法及装置 | |
| CN116611160A (zh) | 基于一段实测弹道参数的无控飞行器在线实时特征参数辨识与弹道预报方法 | |
| de Celis et al. | Neural network-based controller for terminal guidance applied in short-range rockets | |
| CN115774928A (zh) | 基于改进Laplace模型的空间碎片短弧仅测角初定轨优化方法 | |
| Zhao et al. | Multiple model adaptive estimation algorithm for SINS/CNS integrated navigation system | |
| CN112949150A (zh) | 基于变结构自适应多模型箱粒子滤波弹道目标跟踪方法 | |
| CN118484652B (zh) | 基于elm与ukf的飞行器气动参数在线辨识方法 | |
| Machala et al. | Quasi-LPV modelling of a projectile’s behaviour in flight | |
| CN116909303B (zh) | 一种用于临近空间目标跟踪的过程噪声自适应调节方法 | |
| CN116992553B (zh) | 助推滑翔飞行器全程弹道估计方法 | |
| Pei et al. | Mahalanobis distance based adaptive unscented Kalman filter and its application in GPS/MEMS-IMU integration | |
| Manarvi et al. | Application of Kalman Filters in Orbit Determination: A Literature Survey |
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 | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20220916 |