CN112784506A - 一种基于变结构多模型的再入机动弹道目标跟踪算法 - Google Patents

一种基于变结构多模型的再入机动弹道目标跟踪算法 Download PDF

Info

Publication number
CN112784506A
CN112784506A CN202110130306.9A CN202110130306A CN112784506A CN 112784506 A CN112784506 A CN 112784506A CN 202110130306 A CN202110130306 A CN 202110130306A CN 112784506 A CN112784506 A CN 112784506A
Authority
CN
China
Prior art keywords
model
mode
state
covariance
acceleration
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
CN202110130306.9A
Other languages
English (en)
Other versions
CN112784506B (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.)
Air Force Engineering University of PLA
Original Assignee
Air Force Engineering University of PLA
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 Air Force Engineering University of PLA filed Critical Air Force Engineering University of PLA
Priority to CN202110130306.9A priority Critical patent/CN112784506B/zh
Publication of CN112784506A publication Critical patent/CN112784506A/zh
Application granted granted Critical
Publication of CN112784506B publication Critical patent/CN112784506B/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
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
    • G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
    • G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/08—Probabilistic or stochastic CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • 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)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Fluid Mechanics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Evolutionary Biology (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种基于变结构多模型的再入机动弹道目标跟踪算法，包括步骤S1.构建基于变结构多模型的再入机动弹道目标跟踪算法模型；S2.对状态融合估计模型集中的基本模型进行转换模式互补合并，同时对增广模型进行互补交互；S3.对基本模型模式和增广模型进行第一次状态滤波，对模式概率进行第一次更新；S4.对第一次滤波和更新后的状态
Figure DDA0002924867610000011
与协方差
Figure DDA0002924867610000012
进行第一次融合输出；S5.对第一次融合输出结果进行增广模型匹配度检测；S6.更新模型转换概率矩阵；经验证，本算法在进行再入机动导弹目标跟踪的过程中，具有收敛速度和稳态精度高的特点。

Description

一种基于变结构多模型的再入机动弹道目标跟踪算法
技术领域
本发明涉及信号与信息处理技术领域，具体涉及一种基于变结构多模型的再入机动弹道目标跟踪算法。
背景技术
随着弹道导弹技术的不断发展，再入机动已成为弹道导弹突防中的重要手段，为了有效跟踪再入机动弹道目标，李晓榕等在《Survey of maneuvering targettracking.Part II:Motion models of ballistic and space targets》(IEEETransactions on Aerospace and Electronic Systems，2010，46(1):96-119)中对再入机动弹道目标进行了详细的受力分析并给出了带空气动力参数的加速度动力学模型，描述了目标加速度与位置、速度和空气动力参数的关系，在此基础上给出了可用于滤波估计的分段匀加速度(PCA.Piecewise constant acceleration)模型；根据PCA模型，当目标的空气动力参数已知时，再入目标跟踪问题就转化为非线性状态估计问题；当目标动力学参数未知时，再入目标跟踪问题就转化为非线性状态和未知时变参数的联合估计问题；
目前，再入目标跟踪算法研究主要是建立在PCA模型基础上的，利用PCA模型，Ristic B.等在《Performance bounds and comparison of nonlinear filters fortracking aballistic object on reentry》(IEE Proceedings,Radar,Sonar,andNavigation,2003,150(2):65-70)针对空气动力参数未知的弹道再入目标，对比了几种基于PCA模型的非线性滤波算法的跟踪性能；刘也等在《基于双重酉滤波的再入目标实时跟踪》(国防科技大学学报,2011,33(1):81-86)针对机动再入目标，基于PCA模型研究了双重酉滤波算法；Chen Y.等在《Nonlinear filtering for tracking maneuverableballistic missile targets on reentry》(Proceedings of International RadarConference，2009)针对弹道再入目标跟踪问题，研究了PCA模型过程噪声的调整方法，采用归一化的新息平方构成检测统计量并设定检测门限，根据检测统计量超过门限的程度对过程噪声进行调整，但是这种方法只是一种粗略的过程噪声调整方法，调整的实时性和合理性难以保证，过程噪声频繁跃变也会对滤波算法的稳定性造成不利影响；
从状态建模方法来看，PCA模型将目标加速度作为确定性的控制输入量，因此建立在PCA模型基础上的非线性滤波算法缺乏对目标加速度估计误差的直接的修正能力，目标加速度估计精度完全取决于位置、速度和空气动力参数的估计精度，当目标通过自身动力或空气动力参数变化进行机动时，加速度动力学模型给出的加速度估计误差将会增大，进而引起位置、速度和空气动力参数估计误差的增大，预设固定的过程噪声使滤波算法不能快速合理修正位置、速度和空气动力参数估计误差，因而又会引起加速度误差的进一步增大，快速增大的状态-参数估计误差与滤波算法有限修正能力之间的矛盾会导致跟踪算法状态估计误差迅速增大甚至发散，存在再入跟踪算法抗目标机动能力差，难以有效估计目标时变空气动力参数的缺点。
发明内容
针对上述存在的问题，本发明的目的在于提供一种基于变结构多模型的再入机动弹道目标跟踪算法，本算法通过基于模型增广方法和模型驻留-转换结构，对再入机动弹道目标的轨迹进行计算，能够解决现有再入跟踪算法抗目标机动能力差，难以有效估计目标时变空气动力参数的问题，并且经验证，本算法在进行再入机动导弹目标跟踪的过程中，具有收敛速度和稳态精度高的特点。
为了实现上述目的，本发明所采用的技术方案如下：
一种基于变结构多模型的再入机动弹道目标跟踪算法，包括步骤：
S1.在再入机动弹道目标的跟踪过程中，利用PCA算法和PCJ算法构建基于变结构多模型的再入机动弹道目标跟踪算法模型集；
其中目标跟踪算法模型集包括基本模型集和候选模型集，基本模型集与增广模型融合得到状态融合估计模型集；
S2.对状态融合估计模型集中的基本模型进行转换模式互补合并，同时对增广模型进行互补交互；
S3.对基本模型模式和增广模型进行第一次状态滤波，对模式概率进行第一次更新；
S4.对第一次滤波和更新后的状态
Figure BDA0002924867590000021
与协方差
Figure BDA0002924867590000022
进行第一次融合输出；
S5.对第一次融合输出结果进行增广模型匹配度检测；
S6.更新模型转换概率矩阵：根据所有模型模式的概率和似然值，求得模型的概率和似然值，更新模型转换概率矩阵；通过多个周期由S1至S6循环进行处理，实现对导弹再入机动弹道目标进行跟踪。
优选的，在步骤S1所述的基于变结构多模型的再入机动弹道目标跟踪算法模型集中：
(1)若模型在k+1时刻发生转换，则称k时刻的模型处于“转换模式”；
(2)若模型在k+1时刻不发生转换，则称k时刻模型处于“驻留模式”；
(3)按照k时刻模型已经驻留的采样周期数，将其称为“一级驻留模式”、“二级驻留模式”，达到模型最大驻留时间的称为最高级驻留模式；
其中：在计算过程中，设基本模型集中PCA扩展模型编号为1，PCJ扩展模型编号为2，增广PCJ扩展模型编号为3；
模型1的递推关系为：
Figure BDA0002924867590000023
模型2的递推关系为：
Figure BDA0002924867590000024
在式(1)和式(2)中：pk为空气动力参数向量；
Figure BDA0002924867590000025
和
Figure BDA0002924867590000026
分别是协方差矩阵为
Figure BDA0002924867590000027
和
Figure BDA0002924867590000028
的高斯噪声，
Figure BDA0002924867590000029
Figure BDA00029248675900000210
和
Figure BDA00029248675900000211
为空气动力参数模型的过程噪声，协方差分别为
Figure BDA00029248675900000212
和
Figure BDA00029248675900000213
qCV和qCA为经验参数；
模型3从初始时刻起由基本模型集经模型转换递推求得。
优选的，步骤S2所述的对状态融合估计模型集中的基本模型进行转换模式互补合并，同时对增广模型进行互补交互的具体过程为：
S201.模型1转换模式只向模型2输出空气动力参数，其他状态用模型2的一级驻留模式相应状态填充；
S202.向模型3只输出位置和速度，模型1转换模式状态缺少的加速度采用模型3的加速度状态填充；
S203.模型2转换模式只向模型1输出位置和速度，缺少的空气动力参数采用模型2的一级驻留模式的空气动力参数估计值填充；
S204.向模型3输出位置、速度和加速度；
S205.模型3向模型1转换模式输出位置、速度，缺少的空气动力参数采用模型1的一级驻留模式的空气动力参数填充；向模型2转换模式只输出空气动力参数，缺少的状态采用模型2一级驻留模式的相应状态填充；
S206.即得到k-1时刻的联合状态
Figure BDA0002924867590000031
和联合状态协方差
Figure BDA0002924867590000032
为：
Figure BDA0002924867590000033
根据式(3)可得到模型3向模型1输出的联合状态
Figure BDA0002924867590000034
和联合状态协方差
Figure BDA0002924867590000035
为：
Figure BDA0002924867590000036
其中：在式(4)中，
Figure BDA0002924867590000037
和
Figure BDA0002924867590000038
分别是模型i转换模式的合并概率和模型i最高级驻留模式的合并概率，
Figure BDA0002924867590000039
是模型i转换模式的交互输入概率，(·)T应为前一个变量的转置；对于模型i，用Ti表示其转换模式，
Figure BDA00029248675900000310
表示其第l级驻留模式，即模型i是由模型j转换而来且已经驻留了l个采样周期，l＝1,2,…,βij；
根据式(3)可得到基本模型一级驻留模式的输入为：
Figure BDA0002924867590000041
Figure BDA0002924867590000042
其中：y＝[x；p]是由状态、参数组成的联合状态；{x,P}为模型状态向量及协方差；{p,Pp}为模型参数向量及协方差，C为联合状态协方差；
Figure BDA0002924867590000043
是模型1转向模型2的联合状态向量，由模型2的一级驻留模式状态
Figure BDA0002924867590000044
和模型1的转换模式参数向量
Figure BDA0002924867590000045
组合而成，即
Figure BDA0002924867590000046
Figure BDA0002924867590000047
模型1向模型3输出的状态量为
Figure BDA0002924867590000048
Figure BDA0002924867590000049
是模型2向模型1输出的联合状态向量，由模型2的转换模式状态
Figure BDA00029248675900000410
和模型1的一级驻留模式参数向量
Figure BDA00029248675900000411
组合而成，即
Figure BDA00029248675900000412
Figure BDA00029248675900000413
模型2向模型3输出的状态向量为
Figure BDA00029248675900000414
Figure BDA00029248675900000415
是模型3向模型1输出的联合状态，由模型3的状态和模型1的一级驻留模式参数向量
Figure BDA00029248675900000416
构成，即
Figure BDA00029248675900000417
Figure BDA00029248675900000418
是模型3向模型2输出的联合状态，由模型2一级驻留模式状态和模型3的参数向量
Figure BDA00029248675900000419
构成，即
Figure BDA00029248675900000420
优选的，步骤S3所述的对基本模型模式和增广模型进行第一次状态滤波，对模式概率进行第一次更新的过程包括：
S301.各基本模型模式和增广模型的状态滤波和概率更新的目的是获得模式信息，其形式为{y,C,w}，其中联合状态y及联合状态协方差C可经由式(4)可以得出；
S302.引入基本模型模式和增广模型概率第一次更新的计算公式为：
Figure BDA00029248675900000421
上式中：
Figure BDA00029248675900000422
和
Figure BDA00029248675900000423
分别为基本模型转换模式、驻留模式和增广模型的似然值，
Figure BDA0002924867590000051
Figure BDA0002924867590000052
S303.根据贝叶斯全概率公式，由模型i所有似然值得到模型i的似然值的计算公式为：
Figure BDA0002924867590000053
式(6)中：
Figure BDA0002924867590000054
为转换模式概率，
Figure BDA0002924867590000055
为驻留模式概率，
Figure BDA0002924867590000056
为模型转换概率估计值，Zk-1为k时刻以前的量测集合，
Figure BDA0002924867590000057
经状态滤波和概率更新后，获得基本模型转换模式信息
Figure BDA0002924867590000058
基本模型驻留模式信息
Figure BDA0002924867590000059
和增广模型信息
Figure BDA00029248675900000510
优选的，步骤S4所述的对第一次滤波和更新后的状态
Figure BDA00029248675900000511
与协方差
Figure BDA00029248675900000512
进行第一次融合输出的具体过程包括：
S401.由k+1时刻的所有模式概率和模式状态估计值加权计算可获得状态融合输出：
Figure BDA00029248675900000513
由k+1时刻模型1所有模式概率和模式参数估计值可获得参数融合输出：
Figure BDA00029248675900000514
优选的，步骤S5所述的对第一次融合输出结果进行增广模型匹配度检测的具体过程包括：
S501.将k-1时刻的系统状态
Figure BDA0002924867590000061
和k时刻的量测输入到k-1时刻的增广模型中，得到量测预测和量测预测协方差，采用模型失配检测函数：
Figure BDA0002924867590000062
式中：vk为量测残差，Pzz为滤波器输出的量测预测协方差；当模型预设机动频率与目标实际机动匹配时，Dk服从自由度为m的χ2分布，m为量测维数；
S502.确定k-1时刻增广模型是否匹配k时刻的目标运动模式，即目标真实运动状态的变化规律；
根据得到的Dk值进行判断，如果Dk＜3，认为k时刻目标运动模式与k-1时刻的增广模型相同，状态、协方差的第一次融合输出
Figure BDA0002924867590000063
即为k时刻状态和协方差最终融合输出，同时直接转向步骤S6输出；
S503.如果Dk>3，认为k时刻目标运动模式与k-1时刻的增广模型不同，则进行协同方差判断过程：(1)生成模式共同变量和协方差；(2)选择增广模型；(3)对选择出的增广模型进行第二次状态滤波，同时将基本模型模式和增广模型的概率进行第二次更新；(4)状态和协方差第二次融合输出。
优选的，步骤S503(1)所述的生成模式共同变量和协方差协的具体过程包括：
a.选择加速度作为目标模式sk和动力学模型
Figure BDA0002924867590000064
的共同变量，将第一次融合后的状态和协方差
Figure BDA0002924867590000065
中的加速度状态和协方差分量作为目标模式的共同变量，同时将第一次融合后的状态和协方差
Figure BDA0002924867590000066
中的位置、速度分量和候选空气动力参数集合代入机动再入目标的加速度动力学模型公式：
Figure BDA0002924867590000067
其中，(x,y,z)为位置各分量，
Figure BDA0002924867590000068
为速度各分量，
Figure BDA0002924867590000069
为加速度各分量，p为空气动力参数；
b.通过上式可以得到一组对应不同的空气动力参数的加速度值，相应加速度协方差的计算方法为：根据加速度动力学模型给出的加速度与位置、速度的关系构建变量
Figure BDA0002924867590000071
可简写加速度动力学模型为：
Figure BDA0002924867590000072
式中：
Figure BDA0002924867590000073
c.得到相应的状态协方差矩阵Pn为：
Figure BDA0002924867590000074
式中：Pn中的元素由
Figure BDA0002924867590000075
中的相应元素组成；
d.得到加速度协方差
Figure BDA0002924867590000076
计算方法为：
Figure BDA0002924867590000077
式中：l为差分步长，nxn为xn的维数，Sxn＝chol(Pn)为协方差Pn的乔列斯基分解式，Sxn,j为Sxn的第j列。
优选的，步骤S503(2)所述的选择增广模型的具体过程包括：
设已知量测序列Zk、k时刻以前的状态估计所用的模型集序列Mk-1，Mk-1＝{M1,M2,…,Mk-1}，则k时刻状态估计模型集Mk的自适应方法可表示为：
Figure BDA0002924867590000078
式中：
Figure BDA0002924867590000079
为k时刻基本模型集的递推模型集，
Figure BDA00029248675900000710
为k时刻最优增广模型集，
Figure BDA00029248675900000711
为k时刻候选模型集；
a.根据获得的目标模式加速度、加速度协方差和不同空气动力参数对应的加速度、加速度协方差集合，引入KL信息来度量不同空气动力参数对应的加速度与目标模式加速度之间的差异；
b.假定y是加速度共同变量，p[y|sk,Mk-1,Z]和
Figure BDA00029248675900000712
分别是模式加速度和不同空气动力参数对应的加速度的条件概率密度函数，两者KL意义上的距离为：
Figure BDA0002924867590000081
c.对于模型变量服从高斯分布的模型集，y为高斯分布，所以sk和
Figure BDA0002924867590000082
的共同变量y的条件概率密度函数可用y的均值和协方差来表示：
Figure BDA0002924867590000083
Figure BDA0002924867590000084
式中：
Figure BDA0002924867590000085
Figure BDA0002924867590000086
可变为：
Figure BDA0002924867590000087
式中：n是y的维数，tr[·]是[·]的迹，由该公式可以确定KL意义上候选模型集中与目标运动模式最近似的模型为：
Figure BDA0002924867590000088
d.选择与目标模式加速度KL距离最近的加速度，对应的空气动力参数即为增广模型的空气动力参数取值：
Figure BDA0002924867590000089
Figure BDA00029248675900000810
选择出KL意义上的最优增广模型，得到相应的加速度协方差，提取
Figure BDA00029248675900000811
中的加速度滤波估计值和协方差，利用给出的最优模型选取方法即可得到最优的空气动力参数向量p3，将p3代入模型3即形成新的增广模型。
优选的，步骤S503(3)所述的对选择出的增广模型进行第二次状态滤波，同时将基本模型模式和增广模型的概率也进行第二次更新的具体过程包括：
将初始状态、协方差输入增广模型3中进行滤波运算并计算似然值，根据式(3)第二次更新基本模型模式和增广模型的概率，获得：
Figure BDA00029248675900000812
和
Figure BDA00029248675900000813
优选的，步骤S503(4)所述的状态和协方差第二次融合输出的具体过程包括：
a.由k+1时刻所有模式概率和模式状态估计值加权计算可获得状态融合输出为：
Figure BDA0002924867590000091
b.同理由k+1时刻模型1所有模式概率和模式参数估计值可获得参数融合输出为：
Figure BDA0002924867590000092
本发明的有益效果是：本发明公开了一种基于变结构多模型的再入机动弹道目标跟踪算法，与现有技术相比，本发明的改进之处在于：
针对现有技术中存在的再入跟踪算法抗目标机动能力差，难以有效估计目标时变空气动力参数的问题，本发明提出了一种基于变结构多模型的再入机动弹道目标跟踪算法，本算法通过基于模型增广方法和模型驻留-转换结构，对再入机动弹道目标的轨迹进行计算，能够解决现有再入跟踪算法抗目标机动能力差，难以有效估计目标时变空气动力参数的问题；并且通过验证证明，本算法在进行再入机动导弹目标跟踪计算的过程中，具有收敛速度和稳态精度高的优点。
附图说明
图1为本发明基于变结构多模型的跟踪算法流程图。
图2为本发明互补STC-VSAIMM结构图。
图3为本发明实施例1再入机动弹道目标运动状态变化规律图。
图4为本发明实施例1位置估计均方根误差图。
图5为本发明实施例1速度估计均方根误差图。
图6为本发明实施例1空气动力参数估计均值图。
图7为本发明实施例1模型概率随时间变化的情况图。
图8为本发明实施例1不同算法的模型概率随时间变化的情况。
图9为本发明实施例1量测噪声增大后不同算法的状态估计误差均值图。
图10为本发明实施例1量测噪声增大后不同算法的状态估计均方根误差图。
图11为本发明实施例1量测噪声增大不同算法模型概率随时间变化的情况图。
图12为本发明实施例1量测噪声增大后不同算法的量测误差压缩比图。
图13为本发明实施例1量测噪声增大后不同算法模型概率随时间变化的情况图。
其中：在图3中，图(a)代表目标机动再入轨迹图，图(b)代表目标速度曲线图，图(c)代表目标加速度曲线图；在图4中，图(a)代表位置误差均值，图(b)代表速度误差均值图；在图5中，图(a)代表位置均方根误差图，图(b)代表速度均方根误差图；在图6中，图(a)代表阻力参数估计值图，图(b)代表升力参数估计值图；在图8中，图(a)代表互补STC-VSAIMM算法图，图(b)代表互补STC-AMIMM算法图；在图9中，图(a)代表量测噪声增大后的位置误差均值图，图(b)代表量测噪声增大后的速度误差均值图；在图10中，图(a)代表量测噪声增大后的位置均方根误差图，图(b)代表量测噪声增大后的速度均方根误差图；在图11中，图(a)代表量测噪声增大后的阻力参数估计值图，图(b)代表量测噪声增大后的升力参数估计值图；在图12中，图(a)代表量测噪声增大后的阻力参数估计值图，图(b)代表量测噪声增大后的升力参数估计值图；在图13中，图(a)代表量测噪声增大后的互补STC-VSAIMM算法图，图(b)代表量测噪声增大后的互补STC-AMIMM算法图。
具体实施方式
为了使本领域的普通技术人员能更好的理解本发明的技术方案，下面结合附图和实施例对本发明的技术方案做进一步的描述。
参照附图1-13所示的一种基于变结构多模型的再入机动弹道目标跟踪算法，包括步骤：
S1.构建基于变结构多模型的再入机动弹道目标跟踪算法模型集；
在再入机动弹道目标的跟踪过程中，针对该类目标的气动特性分析，根据PCA算法和PCJ算法的各自特点构建模型集；
在相同的状态估计精度条件下，PCA算法可获得更好的参数估计性能，但是由于PCA算法的抗机动性能较差，当目标状态变化剧烈时难以获得较高精度的状态估计；PCJ算法的抗机动性能较好，但是由于加速度误差被滤波算法实时修正，使得参数与量测的误差相关性下降，降低了参数的收敛速度和稳态估计精度，最终也就降低了算法的稳态跟踪精度；因此，本发明中采用PCA扩展模型和PCJ扩展模型构成模型集，同时利用PCA扩展模型较好的参数估计性能和PCJ模型较好的抗机动性能，则能够实现模型优势互补，兼顾算法的抗机动性能和稳态跟踪精度；根据这种思想，目标跟踪算法模型集包括基本模型集和候选模型集，基本模型集由过程噪声较小的分段PCA扩展模型和过程噪声较大的分段匀Jerk(PCJ)扩展模型构成；增广PCJ扩展模型(下称增广模型)为空气动力参数时变的PCJ扩展模型；将基本模型集与增广模型进行状态融合估计便可得到所需要的状态融合估计模型集；
其中，在本基于变结构多模型的再入机动弹道目标跟踪算法模型集中：假设基本模型集中PCA扩展模型编号为1，PCJ扩展模型编号为2，增广PCJ扩展模型编号为3，模型结构如图2所示；
(1)若模型在k+1时刻发生转换，则称k时刻的模型处于“转换模式”；
(2)若模型在k+1时刻不发生转换，则称k时刻模型处于“驻留模式”；
(3)按照k时刻模型已经驻留的采样周期数，将其称为“一级驻留模式”、“二级驻留模式”，达到模型最大驻留时间的称为最高级驻留模式；
模型1的递推关系为：
Figure BDA0002924867590000101
模型2的递推关系为：
Figure BDA0002924867590000111
在式(1)和式(2)中：pk为空气动力参数向量；
Figure BDA0002924867590000112
和
Figure BDA0002924867590000113
分别是协方差矩阵为
Figure BDA0002924867590000114
和
Figure BDA0002924867590000115
的高斯噪声，
Figure BDA0002924867590000116
Figure BDA0002924867590000117
和
Figure BDA0002924867590000118
为空气动力参数模型的过程噪声，协方差分别为
Figure BDA0002924867590000119
和
Figure BDA00029248675900001110
qCV和qCA为经验参数；
模型3从初始时刻起由基本模型集经模型转换递推求得，其递推方法将在S503(2)的详细说明中给出；
S2.对状态融合估计模型集中的基本模型进行转换模式互补合并，同时对增广模型进行互补交互，其具体过程包括：
S201.模型1转换模式只向模型2输出空气动力参数，其他状态用模型2的一级驻留模式相应状态填充，确保不干扰模型2状态滤波；
S202.模型1转换模式向模型3只输出位置和速度，模型1转换模式状态缺少的加速度采用模型3的加速度状态填充；
S203.模型2转换模式只向模型1输出位置和速度，缺少的空气动力参数采用模型2的一级驻留模式的空气动力参数估计值填充；
S204.模型2转换模式向模型3输出位置、速度和加速度；
S205.模型3向模型1转换模式输出位置、速度，缺少的空气动力参数采用模型1的一级驻留模式的空气动力参数填充；向模型2转换模式只输出空气动力参数，缺少的状态采用模型2一级驻留模式的相应状态填充；
S206.即得到k-1时刻的联合状态
Figure BDA00029248675900001111
和联合状态协方差
Figure BDA00029248675900001112
为：
Figure BDA00029248675900001113
根据式(3)可得到模型3向模型1输出的联合状态
Figure BDA00029248675900001114
和联合状态协方差
Figure BDA00029248675900001115
为：
Figure BDA0002924867590000121
其中：在式(4)中，
Figure BDA0002924867590000122
和
Figure BDA0002924867590000123
分别是模型i转换模式的合并概率和模型i最高级驻留模式的合并概率，
Figure BDA0002924867590000124
是模型i转换模式的交互输入概率，(·)T应为前一个变量的转置；对于模型i，用Ti表示其转换模式，
Figure BDA0002924867590000125
表示其第l级驻留模式，即模型i是由模型j转换而来且已经驻留了l个采样周期，l＝1,2,…,βij；
根据式(3)可得到基本模型一级驻留模式的输入为：
Figure BDA0002924867590000126
Figure BDA0002924867590000127
其中：y＝[x；p]是由状态、参数组成的联合状态；{x,P}为模型状态向量及协方差；{p,Pp}为模型参数向量及协方差，C为联合状态协方差；
Figure BDA0002924867590000128
是模型1转向模型2的联合状态向量，由模型2的一级驻留模式状态
Figure BDA0002924867590000129
和模型1的转换模式参数向量
Figure BDA00029248675900001210
组合而成，即
Figure BDA00029248675900001211
Figure BDA00029248675900001212
模型1向模型3输出的状态量为
Figure BDA00029248675900001213
Figure BDA00029248675900001214
是模型2向模型1输出的联合状态向量，由模型2的转换模式状态
Figure BDA00029248675900001215
和模型1的一级驻留模式参数向量
Figure BDA00029248675900001216
组合而成，即
Figure BDA00029248675900001217
Figure BDA00029248675900001218
模型2向模型3输出的状态向量为
Figure BDA00029248675900001219
Figure BDA00029248675900001220
是模型3向模型1输出的联合状态，由模型3的状态和模型1的一级驻留模式参数向量
Figure BDA00029248675900001221
构成，即
Figure BDA00029248675900001222
Figure BDA00029248675900001223
是模型3向模型2输出的联合状态，由模型2一级驻留模式状态和模型3的参数向量
Figure BDA0002924867590000131
构成，即
Figure BDA0002924867590000132
互补交互输入的最终目的是为了确保三种模型(PCA扩展模型、PCJ扩展模型和增广PCJ扩展模型)互为补充而不相互竞争，其中PCA扩展模型主要作用是对空气动力参数进行小范围调整，为系统提供稳态精度保证，PCJ扩展模型的主要作用是为系统提供一定精度的状态估计并确保目标机动情况下状态估计的快速收敛，增广模型则是通过空气动力参数硬切换实现空气动力参数快速追踪，主要为PCJ扩展模型提供较为准确的空气动力参数促进其状态估计的快速收敛和精度的提升；
S3.对基本模型模式和增广模型进行第一次状态滤波，对模式概率进行第一次更新，具体过程包括：
S301.各基本模型模式和增广模型的状态滤波和概率更新的目的是获得模式信息，其形式为{y,C,w}，其中联合状态y及联合状态协方差C可经由式(4)可以得出；
S302.引入基本模型模式和增广模型概率第一次更新的计算公式为：
Figure BDA0002924867590000133
上式中：
Figure BDA0002924867590000134
和
Figure BDA0002924867590000135
分别为基本模型转换模式、驻留模式和增广模型的似然值，
Figure BDA0002924867590000136
Figure BDA0002924867590000137
S303.根据贝叶斯全概率公式，由模型i所有似然值得到模型i的似然值的计算公式为：
Figure BDA0002924867590000138
式(6)中：
Figure BDA0002924867590000139
为转换模式概率，
Figure BDA00029248675900001310
为驻留模式概率，
Figure BDA00029248675900001311
为模型转换概率估计值，Zk-1为k时刻以前的量测集合，
Figure BDA00029248675900001312
经状态滤波和概率更新后，获得基本模型转换模式信息
Figure BDA00029248675900001313
基本模型驻留模式信息
Figure BDA0002924867590000141
和增广模型信息
Figure BDA0002924867590000142
S4.对第一次滤波和更新后的状态
Figure BDA0002924867590000143
与协方差
Figure BDA0002924867590000144
进行第一次融合输出，具体过程包括：
S401.模式状态、参数互补融合输出为：由k+1时刻的所有模式概率和模式状态估计值加权计算可获得状态融合输出：
Figure BDA0002924867590000145
S402.同理由k+1时刻模型1所有模式概率和模式参数估计值可获得参数融合输出：
Figure BDA0002924867590000146
S5.对第一次融合输出结果进行增广模型匹配度检测，具体过程包括：
S501.将k-1时刻的系统状态
Figure BDA0002924867590000147
和k时刻的量测输入到k-1时刻的增广模型中，得到量测预测和量测预测协方差(跟踪算法是从k-1时刻的滤波值开始的，再加上k时刻的量测值，通过滤波算法就，可以实现k时刻目标滤波了，连续时刻进行，就可以实现目标连续跟踪了)，采用模型失配检测函数：
Figure BDA0002924867590000148
式中：vk为量测残差，Pzz为滤波器输出的量测预测协方差；当模型预设机动频率与目标实际机动匹配时，Dk服从自由度为m的χ2分布，m为量测维数；
S502.确定k-1时刻增广模型是否匹配k时刻的目标运动模式，即目标真实运动状态的变化规律；
根据得到的Dk值进行判断，如果Dk＜3，认为k时刻目标运动模式与k-1时刻的增广模型相同，状态、协方差的第一次融合输出
Figure BDA0002924867590000149
即为k时刻状态和协方差最终融合输出，同时直接转向步骤S6输出；
S503.如果Dk>3，认为k时刻目标运动模式与k-1时刻的增广模型不同，则进行协同方差判断过程：
(1)生成模式共同变量和协方差
a.选择加速度作为目标模式sk和动力学模型
Figure BDA0002924867590000151
的共同变量，将第一次融合后的状态和协方差
Figure BDA0002924867590000152
中的加速度状态和协方差分量作为目标模式的共同变量，同时将第一次融合后的状态和协方差
Figure BDA0002924867590000153
中的位置、速度分量和候选空气动力参数集合代入机动再入目标的加速度动力学模型公式：
Figure BDA0002924867590000154
其中，(x,y,z)为位置各分量，
Figure BDA0002924867590000155
为速度各分量，
Figure BDA0002924867590000156
为加速度各分量，p为空气动力参数；
b.通过上式可以得到一组对应不同的空气动力参数的加速度值，相应加速度协方差的计算方法为：根据加速度动力学模型给出的加速度与位置、速度的关系构建变量
Figure BDA0002924867590000157
可简写加速度动力学模型为：
Figure BDA0002924867590000158
式中：
Figure BDA0002924867590000159
c.得到相应的状态协方差矩阵Pn为：
Figure BDA00029248675900001510
式中：Pn中的元素由
Figure BDA00029248675900001511
中的相应元素组成；
d.得到加速度协方差
Figure BDA00029248675900001512
计算方法为：
Figure BDA00029248675900001513
式中：l为差分步长，nxn为xn的维数，Sxn＝chol(Pn)为协方差Pn的乔列斯基分解式，Sxn,j为Sxn的第j列；
(2)选择增广模型
设已知量测序列Zk、k时刻以前的状态估计所用的模型集序列Mk-1，Mk-1＝{M1,M2,…,Mk-1}，则k时刻状态估计模型集Mk的自适应方法可表示为：
Figure BDA0002924867590000161
式中：
Figure BDA0002924867590000162
为k时刻基本模型集的递推模型集，
Figure BDA0002924867590000163
为k时刻最优增广模型集，
Figure BDA0002924867590000164
为k时刻候选模型集；为了讨论的简便，算法采用单一增广模型，即
Figure BDA0002924867590000165
唯一；
a.根据获得的目标模式加速度、加速度协方差和不同空气动力参数对应的加速度、加速度协方差集合，引入KL(Kullback-Leiber)信息来度量不同空气动力参数对应的加速度与目标模式加速度之间的差异；
b.假定y是加速度共同变量，p[y|sk,Mk-1,Z]和
Figure BDA0002924867590000166
分别是模式加速度和不同空气动力参数对应的加速度的条件概率密度函数，两者KL意义上的距离为：
Figure BDA0002924867590000167
c.对于模型变量服从高斯分布的模型集，y为高斯分布，所以sk和
Figure BDA0002924867590000168
的共同变量y的条件概率密度函数可用y的均值和协方差来表示：
Figure BDA0002924867590000169
Figure BDA00029248675900001610
式中：
Figure BDA00029248675900001611
Figure BDA00029248675900001612
可变为：
Figure BDA00029248675900001613
式中：n是y的维数，tr[·]是[·]的迹，由该公式可以确定KL意义上候选模型集中与目标运动模式最近似的模型为：
Figure BDA0002924867590000171
d.选择与目标模式加速度KL距离最近的加速度，对应的空气动力参数即为增广模型的空气动力参数取值：
Figure BDA0002924867590000172
Figure BDA0002924867590000173
选择出KL意义上的最优增广模型，得到相应的加速度协方差，提取
Figure BDA0002924867590000174
中的加速度滤波估计值和协方差，利用给出的最优模型选取方法即可得到最优的空气动力参数向量p3，将p3代入模型3即形成新的增广模型；
(3)对选择出的增广模型进行第二次状态滤波，同时将基本模型模式和增广模型的概率也进行第二次更新
将初始状态、协方差输入增广模型3中进行滤波运算并计算似然值，根据式(3)第二次更新基本模型模式和增广模型的概率，获得：
Figure BDA0002924867590000175
和
Figure BDA0002924867590000176
(4)状态和协方差第二次融合输出
a.模式状态、参数互补融合输出为：由k+1时刻所有模式概率和模式状态估计值加权计算可获得状态融合输出：
Figure BDA0002924867590000177
b.同理由k+1时刻模型1所有模式概率和模式参数估计值可获得参数融合输出：
Figure BDA0002924867590000178
S6.更新模型转换概率矩阵：根据所有模型模式的概率和似然值，求得模型的概率和似然值，更新模型转换概率矩阵；由此，下一周期由S1至S6循环进行处理，即可实现对导弹再入机动弹道目标的连续跟踪，进而求得再入机动弹道目标的运动轨迹，具体计算流程如图1所示。
实施例1
S7实例验证：机动再入弹道目标跟踪
S701.参数设计：算法参数及状态、参数初始值设置：互补STC-AMIMM算法参数设置：模型驻留时间设定为β12＝3，β21＝1；转换概率初始值为π11＝0.95，π22＝0.95，π12＝0.05，π21＝0.05，i≠j；PCA模型过程噪声方差为0.01，PCJ模型的初始过程噪声方差为50；互补STC-VSAIMM算法的模型驻留时间与互补STC-AMIMM算法相同：β12＝3，β21＝1；初始模型转换概率为π11＝0.95，π12＝0.01，π13＝0.04，π22＝0.9，π21＝0.075，π23＝0.025，π33＝0.8，π31＝0.009，π32＝0.001，其中π11设置较大主要是考虑到PCA扩展模型参数估计过程应尽可能避免受到其他模型干扰，π3i设置较小主要是考虑到存在的参数选取误差；候选参数向量集为
Figure BDA0002924867590000181
其中
Figure BDA0002924867590000182
Figure BDA0002924867590000183
C＝3.734×10-5；
状态初始化方法：初始阻力参数、方差分别设为βd0＝70000kg·m-1·s-2和
Figure BDA0002924867590000184
初始升力参数、方差分别设为βc0＝70000kg·m-1·s-2和
Figure BDA0002924867590000185
如图3所示；
S702.设计目标运动场景：目标初始状态为x0＝432km，z0＝88km，V0＝2290m/s，目标初始速度方向与X轴正向夹角为η0＝190°，；
S703.模拟计算(σr＝50m,σθ＝0.01°)
(1)互补STC-AMIMM算法与互补STC-VSAIMM算法的状态估计误差均值和均方根误差如图4和图5所示，从图中可以看出，在1～60s再入目标的阻力参数和升力参数不变的阶段，互补STC-VSAIMM算法状态估计精度高于互补STC-AMIMM算法；在60～130s，再入目标升力参数、阻力参数多次跳变的阶段，互补STC-VSAIMM算法状态估计的收敛速度和稳态精度均高于互补STC-AMIMM算法，跳变瞬间状态估计峰值误差与互补STC-AMIMM算法相当；在130～225s再入目标的阻力参数和升力参数不变的阶段，互补STC-VSAIMM算法的收敛速度和稳态精度高于互补STC-AMIMM算法；
(2)互补STC-AMIMM算法与互补STC-VSAIMM算法的空气动力参数估计值如图6所示，从图中可以明显看出，互补STC-VSAIMM算法的参数估计性能远远好于互补STC-AMIMM算法，这也是互补STC-VSAIMM算法的状态估计性能好于互补STC-AMIMM算法的根本原因。在目标空气动力参数的每个跳变瞬间，由于PCJ扩展模型给出的加速度估计不准确，互补STC-VSAIMM算法增广模型也无法获得准确的空气动力参数和状态估计，此时状态和空气动力参数的估计完全依赖基本模型集，因此空气动力参数和状态估计误差与互补STC-AMIMM算法相似；在参数跳变后暂时稳定的阶段，由于PCJ扩展模型的加速度能够快速收敛到一定的精度的估计上，因此参数选取误差迅速减小，增广模型的状态估计精度迅速提高，也就提升了算法的状态估计精度；
(3)互补STC-AMIMM算法与互补STC-VSAIMM算法的量测误差压缩比如图7所示，从图中可以看出，在滤波器进入稳定跟踪后，互补STC-VSAIMM算法没有出现跟踪失效，而互补STC-AMIMM算法在第90s附近则出现了跟踪失效；在量测误差压缩能力方面，互补STC-VSAIMM算法在目标跟踪全程均强于互补STC-AMIMM算法；
(4)互补STC-AMIMM算法与互补STC-VSAIMM算法在观测时间内状态估计平均误差、跟踪有效度和计算时间如表1所示：
表1：算法性能对比(σr＝50m,σθ＝0.01°)
Figure BDA0002924867590000191
从表中可以看出，互补STC-VSAIMM算法的状态估计精度和跟踪有效度均高于互补STC-AMIMM算法；相比互补STC-AMIMM算法，互补STC-VSAIMM算法的计算复杂度也更高，这是由于算法需要增加增广模型匹配滤波器，在目标空气动力参数变化时还需运行增广模型选取算法，选取算法的运算量则跟候选模型的数量成正比；
(5)互补STC-AMIMM算法和互补STC-VSAIMM算法内部各模型概率在跟踪过程中的变化情况如图8所示，从图中可以看出，两种算法的PCA扩展模型概率变化情况相近，对比图3(c)可以看出，两种算法的PCA扩展模型在目标加速度较小阶段均获得了最大的概率加权；两种算法内部模型概率变化的区别体现在目标加速度变化较为剧烈的阶段(60～80s,90～110s)，其中互补STC-VSAIMM算法的增广模型概率最大，而互补STC-AMIMM算法的PCJ扩展模型概率最大，即此时互补STC-VSAIMM算法主要依靠准确的状态递推方程提高算法精度，而互补STC-AMIMM算法主要依靠PCJ扩展模型匹配滤波器的修正能力跟踪目标状态变化，由于滤波修正会引入一定的量测噪声，并且滤波增益不能完全合理的补偿状态递推误差，因此互补STC-AMIMM算法的状态跟踪精度低于互补STC-VSAIMM算法；
S704.模拟计算(σr＝200m,σθ＝0.04°)
为了分析量测噪声增大对互补STC-VSAIMM算法性能的影响，将距离和角度量测噪声标准差增大为原来的4倍，即距离噪声标准差增大为200m，角度噪声标准差增大为0.04°，进行以下比较计算；
(1)互补STC-AMIMM算法和互补STC-VSAIMM算法量测噪声增大后的状态估计均值和均方根误差如图9和图10所示，对比图4和图5可以看出，随着量测噪声的增大，互补STC-VSAIMM算法在初始状态突变时刻(第60s附近)的位置误差峰值误差高于互补STC-AMIMM算法；在其他阶段，互补STC-VSAIMM算法的状态估计总体精度均高于互补STC-AMIMM算法；
(2)互补STC-AMIMM算法和互补STC-VSAIMM算法量测噪声增大后的空气动力参数估计值如图11所示，对比图6可以看出，量测噪声的增大对两种算法的参数估计均造成较大影响，其中互补STC-AMIMM算法基本失去了对目标空气动力参数的追踪能力，互补STC-VSAIMM算法的参数估计精度虽然也大幅降低，但是仍然能够追踪目标空气动力参数变化的趋势，并且对变化幅度较大的升力参数估计性能较好；
(3)互补STC-AMIMM算法和互补STC-VSAIMM算法量测噪声增大后的量测误差压缩比如图12所示，对比图7可以看出，在滤波器进入稳定跟踪后，两种算法均实现了全程有效跟踪且量测误差压缩比降低，其中互补STC-VSAIMM算法的量测误差压缩比相对下降更快，表明量测噪声对互补STC-VSAIMM算法造成的影响相对较小；
(4)互补STC-AMIMM算法和互补STC-VSAIMM算法量测噪声增大后的两种算法在观测时间内状态估计平均误差如表2所示：
表2：算法性能对比(σr＝200m,σθ＝0.04°)
Figure BDA0002924867590000201
对比表1可以看出，量测噪声增大引起的互补STC-VSAIMM算法状态估计平均误差增幅小于互补STC-AMIMM，表明了互补STC-VSAIMM算法具有较强的量测噪声抑制能力。
(5)互补STC-AMIMM算法和互补STC-VSAIMM算法量测噪声增大后的两种算法内部各模型概率在跟踪过程中的变化情况如图13所示，对比图8可以看出，随着量测噪声的增大，模型失配引起的状态估计误差在总的状态估计误差中的比重降低，模型似然值差别变小，因此目标空气动力参数跳变阶段模型概率更加接近。
综上各实验结果，均证明了本发明所述基于变结构多模型的再入机动弹道目标跟踪算法在进行再入机动弹道目标跟踪过程中具有好的收敛性和鲁棒性。以上显示和描述了本发明的基本原理、主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (10)

1.一种基于变结构多模型的再入机动弹道目标跟踪算法，其特征在于：包括步骤：
S1.在再入机动弹道目标的跟踪过程中，利用PCA算法和PCJ算法构建基于变结构多模型的再入机动弹道目标跟踪算法模型集；
其中：目标跟踪算法模型集包括基本模型集和候选模型集，基本模型集与增广模型融合得到状态融合估计模型集；
S2.对状态融合估计模型集中的基本模型进行转换模式互补合并，同时对增广模型进行互补交互；
S3.对基本模型模式和增广模型进行第一次状态滤波，对模式概率进行第一次更新；
S4.对第一次滤波和更新后的状态
Figure FDA0002924867580000011
与协方差
Figure FDA0002924867580000012
进行第一次融合输出；
S5.对第一次融合输出结果进行增广模型匹配度检测；
S6.更新模型转换概率矩阵：根据所有模型模式的概率和似然值，求得模型的概率和似然值，更新模型转换概率矩阵；通过多个周期由S1至S6循环进行处理，实现对导弹再入机动弹道目标进行跟踪。
2.根据权利要求1所述的一种基于变结构多模型的再入机动弹道目标跟踪算法，其特征在于：在步骤S1所述的基于变结构多模型的再入机动弹道目标跟踪算法模型集中：
(1)若模型在k+1时刻发生转换，则称k时刻的模型处于“转换模式”；
(2)若模型在k+1时刻不发生转换，则称k时刻模型处于“驻留模式”；
(3)按照k时刻模型已经驻留的采样周期数，将其称为“一级驻留模式”、“二级驻留模式”，达到模型最大驻留时间的称为最高级驻留模式；
其中：在计算过程中，设基本模型集中PCA扩展模型编号为1，PCJ扩展模型编号为2，增广PCJ扩展模型编号为3；
模型1的递推关系为：
Figure FDA0002924867580000021
模型2的递推关系为：
Figure FDA0002924867580000022
在式(1)和式(2)中：pk为空气动力参数向量；
Figure FDA0002924867580000023
和
Figure FDA0002924867580000024
分别是协方差矩阵为
Figure FDA0002924867580000025
和
Figure FDA0002924867580000026
的高斯噪声，
Figure FDA0002924867580000027
Figure FDA0002924867580000028
Figure FDA0002924867580000029
和
Figure FDA00029248675800000210
为空气动力参数模型的过程噪声，协方差分别为
Figure FDA00029248675800000211
和
Figure FDA00029248675800000212
qCV和qCA为经验参数；
模型3从初始时刻起由基本模型集经模型转换递推求得。
3.根据权利要求2所述的一种基于变结构多模型的再入机动弹道目标跟踪算法，其特征在于：步骤S2所述的对状态融合估计模型集中的基本模型进行转换模式互补合并，同时对增广模型进行互补交互的具体过程为：
S201.模型1转换模式只向模型2输出空气动力参数，其他状态用模型2的一级驻留模式相应状态填充；
S202.向模型3只输出位置和速度，模型1转换模式状态缺少的加速度采用模型3的加速度状态填充；
S203.模型2转换模式只向模型1输出位置和速度，缺少的空气动力参数采用模型2的一级驻留模式的空气动力参数估计值填充；
S204.向模型3输出位置、速度和加速度；
S205.模型3向模型1转换模式输出位置、速度，缺少的空气动力参数采用模型1的一级驻留模式的空气动力参数填充；向模型2转换模式只输出空气动力参数，缺少的状态采用模型2一级驻留模式的相应状态填充；
S206.即得到k-1时刻的联合状态
Figure FDA0002924867580000031
和联合状态协方差
Figure FDA0002924867580000032
为：
Figure FDA0002924867580000033
根据式(3)可得到模型3向模型1输出的联合状态
Figure FDA0002924867580000034
和联合状态协方差
Figure FDA0002924867580000035
为：
Figure FDA0002924867580000036
其中：在式(4)中，
Figure FDA0002924867580000037
和
Figure FDA0002924867580000038
分别是模型i转换模式的合并概率和模型i最高级驻留模式的合并概率，
Figure FDA0002924867580000039
是模型i转换模式的交互输入概率，(·)T应为前一个变量的转置；对于模型i，用Ti表示其转换模式，
Figure FDA00029248675800000310
表示其第l级驻留模式，即模型i是由模型j转换而来且已经驻留了l个采样周期，l＝1,2,…,βij；
根据式(3)可得到基本模型一级驻留模式的输入为：
Figure FDA0002924867580000041
Figure FDA0002924867580000042
其中：y＝[x；p]是由状态、参数组成的联合状态；{x,P}为模型状态向量及协方差；{p,Pp}为模型参数向量及协方差，C为联合状态协方差；
Figure FDA0002924867580000043
是模型1转向模型2的联合状态向量，由模型2的一级驻留模式状态
Figure FDA0002924867580000044
和模型1的转换模式参数向量
Figure FDA0002924867580000045
组合而成，即
Figure FDA0002924867580000046
Figure FDA0002924867580000047
模型1向模型3输出的状态量为
Figure FDA0002924867580000048
Figure FDA0002924867580000049
是模型2向模型1输出的联合状态向量，由模型2的转换模式状态
Figure FDA00029248675800000410
和模型1的一级驻留模式参数向量
Figure FDA00029248675800000411
组合而成，即
Figure FDA00029248675800000412
Figure FDA00029248675800000413
模型2向模型3输出的状态向量为
Figure FDA00029248675800000414
Figure FDA00029248675800000415
是模型3向模型1输出的联合状态，由模型3的状态和模型1的一级驻留模式参数向量
Figure FDA00029248675800000416
构成，即
Figure FDA00029248675800000417
Figure FDA00029248675800000418
是模型3向模型2输出的联合状态，由模型2一级驻留模式状态和模型3的参数向量
Figure FDA00029248675800000419
构成，即
Figure FDA00029248675800000420
4.根据权利要求2所述的一种基于变结构多模型的再入机动弹道目标跟踪算法，其特征在于：步骤S3所述的对基本模型模式和增广模型进行第一次状态滤波，对模式概率进行第一次更新的过程包括：
S301.各基本模型模式和增广模型的状态滤波和概率更新的目的是获得模式信息，其形式为{y,C,w}，其中联合状态y及联合状态协方差C可经由式(4)可以得出；
S302.引入基本模型模式和增广模型概率第一次更新的计算公式为：
Figure FDA0002924867580000051
上式中：
Figure FDA0002924867580000052
和
Figure FDA0002924867580000053
分别为基本模型转换模式、驻留模式和增广模型的似然值，
Figure FDA0002924867580000054
Figure FDA0002924867580000055
S303.根据贝叶斯全概率公式，由模型i所有似然值得到模型i的似然值的计算公式为：
Figure FDA0002924867580000056
式(6)中：
Figure FDA0002924867580000057
为转换模式概率，
Figure FDA0002924867580000058
为驻留模式概率，
Figure FDA0002924867580000059
为模型转换概率估计值，Zk-1为k时刻以前的量测集合，
Figure FDA00029248675800000510
经状态滤波和概率更新后，获得基本模型转换模式信息
Figure FDA00029248675800000511
基本模型驻留模式信息
Figure FDA00029248675800000512
和增广模型信息
Figure FDA00029248675800000513
5.根据权利要求2所述的一种基于变结构多模型的再入机动弹道目标跟踪算法，其特征在于：步骤S4所述的对第一次滤波和更新后的状态
Figure FDA00029248675800000514
与协方差
Figure FDA00029248675800000515
进行第一次融合输出的具体过程包括：
S401.由k+1时刻的所有模式概率和模式状态估计值加权计算可获得状态融合输出：
Figure FDA00029248675800000516
S402.同理由k+1时刻模型1所有模式概率和模式参数估计值可获得参数融合输出：
Figure FDA0002924867580000061
6.根据权利要求2所述的一种基于变结构多模型的再入机动弹道目标跟踪算法，其特征在于：步骤S5所述的对第一次融合输出结果进行增广模型匹配度检测的具体过程包括：
S501.将k-1时刻的系统状态
Figure FDA0002924867580000062
和k时刻的量测输入到k-1时刻的增广模型中，得到量测预测和量测预测协方差，采用模型失配检测函数：
Figure FDA0002924867580000063
式中：vk为量测残差，Pzz为滤波器输出的量测预测协方差；当模型预设机动频率与目标实际机动匹配时，Dk服从自由度为m的χ2分布，m为量测维数；
S502.确定k-1时刻增广模型是否匹配k时刻的目标运动模式，即目标真实运动状态的变化规律；
根据得到的Dk值进行判断，如果Dk＜3，认为k时刻目标运动模式与k-1时刻的增广模型相同，状态、协方差的第一次融合输出
Figure FDA0002924867580000064
即为k时刻状态和协方差最终融合输出，同时直接转向步骤S6输出；
S503.如果Dk>3，认为k时刻目标运动模式与k-1时刻的增广模型不同，则进行协同方差判断过程：(1)生成模式共同变量和协方差；(2)选择增广模型；(3)对选择出的增广模型进行第二次状态滤波，同时将基本模型模式和增广模型的概率进行第二次更新；(4)状态和协方差第二次融合输出。
7.根据权利要求6所述的一种基于变结构多模型的再入机动弹道目标跟踪算法，其特征在于：步骤S503(1)所述的生成模式共同变量和协方差协的具体过程包括：
a.选择加速度作为目标模式sk和动力学模型
Figure FDA0002924867580000071
的共同变量，将第一次融合后的状态和协方差
Figure FDA0002924867580000072
中的加速度状态和协方差分量作为目标模式的共同变量，同时将第一次融合后的状态和协方差
Figure FDA0002924867580000073
中的位置、速度分量和候选空气动力参数集合代入机动再入目标的加速度动力学模型公式：
Figure FDA0002924867580000074
其中，(x,y,z)为位置各分量，
Figure FDA0002924867580000075
为速度各分量，
Figure FDA0002924867580000076
为加速度各分量，p为空气动力参数；
b.通过上式可以得到一组对应不同的空气动力参数的加速度值，相应加速度协方差的计算方法为：根据加速度动力学模型给出的加速度与位置、速度的关系构建变量
Figure FDA0002924867580000077
可简写加速度动力学模型为：
Figure FDA0002924867580000078
式中：
Figure FDA0002924867580000079
c.得到相应的状态协方差矩阵Pn为：
Figure FDA0002924867580000081
式中：Pn中的元素由
Figure FDA0002924867580000082
中的相应元素组成；
d.得到加速度协方差
Figure FDA0002924867580000083
计算方法为：
Figure FDA0002924867580000084
式中：l为差分步长，nxn为xn的维数，Sxn＝chol(Pn)为协方差Pn的乔列斯基分解式，Sxn,j为Sxn的第j列。
8.根据权利要求6所述的一种基于变结构多模型的再入机动弹道目标跟踪算法，其特征在于：步骤S503(2)所述的选择增广模型的具体过程包括：
设已知量测序列Zk、k时刻以前的状态估计所用的模型集序列Mk-1，Mk-1＝{M1,M2,…,Mk-1}，则k时刻状态估计模型集Mk的自适应方法可表示为：
Figure FDA0002924867580000085
式中：
Figure FDA0002924867580000086
为k时刻基本模型集的递推模型集，
Figure FDA0002924867580000087
为k时刻最优增广模型集，
Figure FDA0002924867580000088
Figure FDA0002924867580000089
为k时刻候选模型集；
a.根据获得的目标模式加速度、加速度协方差和不同空气动力参数对应的加速度、加速度协方差集合，引入KL信息来度量不同空气动力参数对应的加速度与目标模式加速度之间的差异；
b.假定y是加速度共同变量，p[y|sk,Mk-1,Z]和
Figure FDA00029248675800000810
分别是模式加速度和不同空气动力参数对应的加速度的条件概率密度函数，两者KL意义上的距离为：
Figure FDA0002924867580000091
c.对于模型变量服从高斯分布的模型集，y为高斯分布，所以sk和mk j的共同变量y的条件概率密度函数可用y的均值和协方差来表示：
Figure FDA0002924867580000092
Figure FDA0002924867580000093
式中：
Figure FDA0002924867580000094
Figure FDA0002924867580000095
可变为：
Figure FDA0002924867580000096
式中：n是y的维数，tr[·]是[·]的迹，由该公式可以确定KL意义上候选模型集中与目标运动模式最近似的模型为：
Figure FDA0002924867580000097
d.选择与目标模式加速度KL距离最近的加速度，对应的空气动力参数即为增广模型的空气动力参数取值：
Figure FDA0002924867580000098
Figure FDA0002924867580000099
选择出KL意义上的最优增广模型，得到相应的加速度协方差，提取
Figure FDA00029248675800000910
中的加速度滤波估计值和协方差，利用给出的最优模型选取方法即可得到最优的空气动力参数向量p3，将p3代入模型3即形成新的增广模型。
9.根据权利要求6所述的一种基于变结构多模型的再入机动弹道目标跟踪算法，其特征在于：步骤S503(3)所述的对选择出的增广模型进行第二次状态滤波，同时将基本模型模式和增广模型的概率也进行第二次更新的具体过程包括：
将初始状态、协方差输入增广模型3中进行滤波运算并计算似然值，根据式(3)第二次更新基本模型模式和增广模型的概率，获得：
Figure FDA0002924867580000101
和
Figure FDA0002924867580000102
10.根据权利要求6所述的一种基于变结构多模型的再入机动弹道目标跟踪算法，其特征在于：步骤S503(4)所述的状态和协方差第二次融合输出的具体过程包括：
a.由k+1时刻所有模式概率和模式状态估计值加权计算可获得状态融合输出为：
Figure FDA0002924867580000103
b.同理由k+1时刻模型1所有模式概率和模式参数估计值可获得参数融合输出为：
Figure FDA0002924867580000104
CN202110130306.9A 2021-01-29 2021-01-29 一种基于变结构多模型的再入机动弹道目标跟踪算法 Active CN112784506B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110130306.9A CN112784506B (zh) 2021-01-29 2021-01-29 一种基于变结构多模型的再入机动弹道目标跟踪算法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110130306.9A CN112784506B (zh) 2021-01-29 2021-01-29 一种基于变结构多模型的再入机动弹道目标跟踪算法

Publications (2)

Publication Number Publication Date
CN112784506A true CN112784506A (zh) 2021-05-11
CN112784506B CN112784506B (zh) 2023-04-07

Family

ID=75760031

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110130306.9A Active CN112784506B (zh) 2021-01-29 2021-01-29 一种基于变结构多模型的再入机动弹道目标跟踪算法

Country Status (1)

Country Link
CN (1) CN112784506B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115390560A (zh) * 2022-08-18 2022-11-25 哈尔滨工业大学 一种基于混合网格多模型的地面目标轨迹跟踪方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7026980B1 (en) * 2005-03-04 2006-04-11 Lockheed Martin Corporation Missile identification and tracking system and method
CN106933106A (zh) * 2016-05-26 2017-07-07 哈尔滨工程大学 一种基于模糊控制多模型算法的目标跟踪方法
WO2018010099A1 (zh) * 2016-07-12 2018-01-18 深圳大学 一种用于跟踪转弯机动目标的方法及其系统
CN107704432A (zh) * 2017-07-28 2018-02-16 西安理工大学 一种转移概率自适应的交互多模型目标跟踪方法
CN108762078A (zh) * 2018-06-01 2018-11-06 福州大学 一种曲线轨迹跟踪控制器的设计方法
CN109633590A (zh) * 2019-01-08 2019-04-16 杭州电子科技大学 基于gp-vsmm-jpda的扩展目标跟踪方法
CN111797478A (zh) * 2020-07-27 2020-10-20 北京电子工程总体研究所 一种基于变结构多模型的强机动目标跟踪方法
CN112257259A (zh) * 2020-10-21 2021-01-22 中国人民解放军战略支援部队信息工程大学 基于改进自治多模型的弹道导弹全程弹道估计方法及系统

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7026980B1 (en) * 2005-03-04 2006-04-11 Lockheed Martin Corporation Missile identification and tracking system and method
CN106933106A (zh) * 2016-05-26 2017-07-07 哈尔滨工程大学 一种基于模糊控制多模型算法的目标跟踪方法
WO2018010099A1 (zh) * 2016-07-12 2018-01-18 深圳大学 一种用于跟踪转弯机动目标的方法及其系统
CN107704432A (zh) * 2017-07-28 2018-02-16 西安理工大学 一种转移概率自适应的交互多模型目标跟踪方法
CN108762078A (zh) * 2018-06-01 2018-11-06 福州大学 一种曲线轨迹跟踪控制器的设计方法
CN109633590A (zh) * 2019-01-08 2019-04-16 杭州电子科技大学 基于gp-vsmm-jpda的扩展目标跟踪方法
CN111797478A (zh) * 2020-07-27 2020-10-20 北京电子工程总体研究所 一种基于变结构多模型的强机动目标跟踪方法
CN112257259A (zh) * 2020-10-21 2021-01-22 中国人民解放军战略支援部队信息工程大学 基于改进自治多模型的弹道导弹全程弹道估计方法及系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
FIRAT KUMRU: "Ballistic target tracking using multiple model Kalman filter with a priori ballistic information", 《IEEE XPLORE》 *
许红: "基于自适应的增广状态-交互式多模型的机动目标跟踪算法", 《电子与信息学报》 *
陈映等: "弹道式再入目标跟踪方法对比分析", 《系统工程与电子技术》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115390560A (zh) * 2022-08-18 2022-11-25 哈尔滨工业大学 一种基于混合网格多模型的地面目标轨迹跟踪方法
CN115390560B (zh) * 2022-08-18 2023-09-15 哈尔滨工业大学 一种基于混合网格多模型的地面目标轨迹跟踪方法

Also Published As

Publication number Publication date
CN112784506B (zh) 2023-04-07

Similar Documents

Publication Publication Date Title
CN111985093B (zh) 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法
CN111291471B (zh) 一种基于l1正则无迹变换的约束多模型滤波方法
CN111460636B (zh) 不完全量测数据驱动下的机动扩展目标混合交互式强跟踪滤波方法
CN111693984B (zh) 一种改进的ekf-ukf动目标跟踪方法
CN111913175A (zh) 一种传感器短暂失效下带补偿机制的水面目标跟踪方法
CN111969979B (zh) 一种最小误差熵cdkf滤波器方法
CN113253050A (zh) 一种基于鲸鱼优化卡尔曼滤波算法的行波故障测距方法
CN114819068A (zh) 一种混合型目标航迹预测方法及系统
CN112784506B (zh) 一种基于变结构多模型的再入机动弹道目标跟踪算法
CN112034445B (zh) 基于毫米波雷达的车辆运动轨迹跟踪方法和系统
CN111722213B (zh) 一种机动目标运动参数的纯距离提取方法
CN111562737B (zh) 一种模糊逻辑控制的改进agimm跟踪方法
CN113030940A (zh) 一种转弯机动下的多星凸型扩展目标跟踪方法
CN116047495B (zh) 一种用于三坐标雷达的状态变换融合滤波跟踪方法
CN110989341B (zh) 一种约束辅助粒子滤波方法及目标跟踪方法
CN114445456B (zh) 基于部分模型的数据驱动智能机动目标跟踪方法及装置
CN113190960B (zh) 一种基于非等维状态混合估计的并行imm机动目标跟踪方法
CN115544425A (zh) 一种基于目标信噪比特征估计的鲁棒多目标跟踪方法
CN117784114B (zh) 异常噪声下基于混合熵的不规则扩展目标跟踪方法
Han et al. Improved adaptive unscented Kalman filter algorithm for target tracking
Crouse Strong tracking filters: derivation and improved heuristic
CN111504326B (zh) 一种基于t分布的鲁棒glmb多目标跟踪方法
CN113532422B (zh) 一种基于距离图和数据清洗的多传感器航迹融合方法
Jianxing et al. A Radar Target Tracking Method based on Coordinate Transformation KF
Zhang et al. Multiple Sensor Track Fusion Algorithm Based on LSTM Network

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