CN111931287A - 一种临近空间高超声速目标轨迹预测方法 - Google Patents

一种临近空间高超声速目标轨迹预测方法 Download PDF

Info

Publication number
CN111931287A
CN111931287A CN202010640882.3A CN202010640882A CN111931287A CN 111931287 A CN111931287 A CN 111931287A CN 202010640882 A CN202010640882 A CN 202010640882A CN 111931287 A CN111931287 A CN 111931287A
Authority
CN
China
Prior art keywords
model
signal
modeling
order
prediction
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
CN202010640882.3A
Other languages
English (en)
Other versions
CN111931287B (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.)
Harbin Institute of Technology
Beijing Institute of Electronic System Engineering
Original Assignee
Harbin Institute of Technology
Beijing Institute of Electronic System Engineering
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 Harbin Institute of Technology, Beijing Institute of Electronic System Engineering filed Critical Harbin Institute of Technology
Priority to CN202010640882.3A priority Critical patent/CN111931287B/zh
Publication of CN111931287A publication Critical patent/CN111931287A/zh
Application granted granted Critical
Publication of CN111931287B publication Critical patent/CN111931287B/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/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Automation & Control Theory (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

一种临近空间高超声速目标轨迹预测方法，属于临近空间高超声速目标轨迹预测领域，本发明为了解决现有技术中通过制导规律在线辨识、拟合外推或模板匹配等方法实现轨迹预测，在临近高超声速目标防御过程中，目标的动力学模型等是未知的，且高超声速目标制导律复杂多变，在线估计困难且误差较大的问题，本发明通过将所得历史弹道数据分解为弹道趋势信号和弹道周期跳跃信号，并分别对趋势信号进行建模，对周期跳跃信号进行建模，最后对步骤二建立的趋势信号模型和步骤三建立的周期跳跃信号模型进行叠合得到弹道完整参数化模型，然后基于完整模型外推实现轨迹预测，本发明主要适用于高超声速目标横向机动轨迹预测、目标速度预测等方面。

Description

一种临近空间高超声速目标轨迹预测方法
技术领域
本发明属于临近空间高超声速目标轨迹预测领域，具体涉及一种临近空间高超声速目标轨迹预测方法。
背景技术
临近空间高超声速飞行器是指飞行于20～100km高度空域，飞行马赫数大于5，具有执行快速打击任务能力的飞行器。相对于传统弹道导弹目标，该类飞行器具有速度快、航程远、机动能力强、打击精度高的特点。借助于优越的气动性能，该类飞行器具有强大的横向机动、纵向跳跃机动能力，对现代防御体系构成了巨大威胁。轨迹精确预测是临近空间高超声速目标有效拦截的基础。因此，开展轨迹预测方法研究对临近空间高超声速目标防御具有十分重要的战略意义。
由于气动力等非保守力作用，传统弹道导弹轨迹预测方法难以满足高超声速目标轨迹预测的要求。现在对临近空间高超声速目标轨迹预测的研究还很少，且大多基于目标动力学模型等先验信息基础上，通过制导规律在线辨识、拟合外推或模板匹配等方法实现轨迹预测。然而，在临近高超声速目标防御过程中，目标的动力学模型等是未知的，且高超声速目标制导律复杂多变，在线估计困难且误差较大。因此，本发明把可观测的目标弹道数据作为研究对象，分析其内在规律并给出参数化描述，实现对目标运动规律建模，从而开发一种不依赖目标动力学模型等先验信息和制导律估计的轨迹预测方法。
针对数据内在规律的挖掘，信号处理领域有着丰富的研究成果。高超声速目标弹道数据表现出明显的趋势下降和周期跳跃特点，非平稳信号处理中的集成经验模态分解法能把非平稳数据分解为不同时间尺度的局部特征信号，适合于分析非线性、非平稳信号序列。而自回归(AR)建模方法能有效处理平稳信号并给出其参数化模型。在非平稳信号处理领域，时变自回归(TVAR)建模方法能给出动态系统的参数化模型，在非平稳随机信号分析和建模方面具有优势被广泛应用。因此，本发明使用上述信号处理方法对历史弹道数据进行分析和分解，并对不同模态的弹道数据分别建模以实现目标轨迹预报精度。
发明内容
本发明为了解决现有技术中通过制导规律在线辨识、拟合外推或模板匹配等方法实现轨迹预测，在临近高超声速目标防御过程中，目标的动力学模型等是未知的，且高超声速目标制导律复杂多变，在线估计困难且误差较大的问题，进而提供一种临近空间高超声速目标轨迹预测方法：
一种临近空间高超声速目标轨迹预测方法，所述步骤是通过如下步骤实现的:
步骤一：将所得历史弹道数据分解为弹道趋势信号和弹道周期跳跃信号；
步骤二：对趋势信号进行建模；
步骤三：对周期跳跃信号进行建模；
步骤四：对步骤二建立的趋势信号模型和步骤三建立的周期跳跃信号模型进行叠合得到弹道完整参数化模型，然后基于完整模型外推实现轨迹预测；
进一步地，所述步骤一的具体内容为：采用集成经验模态分解方法将所得历史弹道数据分解为弹道趋势信号和弹道周期跳跃信号，其中集成经验模态分解方法如下：
步骤一一：在弹道序列X(n)中分别加入K组高斯白噪声序列Gk，得到加噪信号Xk(n)，即：
Xk(n)＝X(n)+Gk (1)
式中，k＝1,…,K
步骤一二：对每组加噪信号Xk(n)进行EMD分解，得到I个固有模态函数IMF分量和残余项：
Figure RE-GDA0002712173150000021
步骤一三：把K次EMD分解的输出均值作为集成经验模态分解的分解结果：
Figure RE-GDA0002712173150000022
步骤一四：原始弹道数据可表示为：
X(n)＝T(n)+R(n) (4)
Figure RE-GDA0002712173150000023
式中，残余项R(n)即为弹道数据的趋势信号，T(n)即为弹道数据的非平稳周期跳跃信号；
进一步地，所述步骤二的具体内容为：利用自回归(AR)建模方法对趋势信号进行建模，模型阶数和模型系数决定了建模精度，其中自回归(AR)建模方法如下：
步骤二一：自回归模型系数估计；
步骤二二：自回归模型阶数选择。
进一步地，所述步骤二一的具体内容如下：
步骤A：令向前、向后预测误差功率之和为最小：
Figure RE-GDA0002712173150000031
式中，上标f表示前向预测(Forward prediction)，上标b表示后向预测(Backwardprediction)，上式中
Figure RE-GDA0002712173150000032
式中：上标f表示前向预测(Forward prediction)，上标b表示后向预测(Backwardprediction)，
Figure RE-GDA0002712173150000033
为前向预测误差，
Figure RE-GDA0002712173150000034
为后向预测误差，ρf为前向预测误差功率，ρb为后向预测误差功率；
则向前、向后观测误差分别为：
Figure RE-GDA0002712173150000035
上式中为了保证不超出给定数据范围，求和范围为p≤n≤N-1；
反射系数kp的递推公式形式如下：
Figure RE-GDA0002712173150000036
式中：rxx(p)为自相关函数，
Figure RE-GDA0002712173150000037
为预测误差，
使预测误差功率之和为最小时的反射系数kp，令：
Figure RE-GDA0002712173150000038
将式(8)带人式(12)中，可得：
Figure RE-GDA0002712173150000041
步骤B：根据步骤A的推导过程，利用Burg递推算法求AR模型参数：
步骤b1：给定初始条件：
Figure RE-GDA0002712173150000042
步骤b2：利用式(13)计算kp：
步骤b3：更新预测误差平均功率：
ρp＝(1-|kp|2)ρp-1 (13)
步骤b4：更新模型系数：
Figure RE-GDA0002712173150000043
步骤b5：更新向前、向后观测误差：
Figure RE-GDA0002712173150000044
步骤b6：令p＝p+1，返回步骤b2，直到所需模型阶数为止。
进一步地，步骤二二中的模型系数直接关系着建模和预测精度，本步骤基于Burg算法确定自回归模型系数，具体过程为：
步骤C：确实模型形式，趋势信号RD(n)是线性、慢变、低频信号，因此自回归模型的形式如下：
Figure RE-GDA0002712173150000045
式中，x(n)是均值为零的信号序列，n＝1,…,N，p为模型阶数，ak为模型系数，u(n)为零均值、方差为
Figure RE-GDA0002712173150000046
的白噪声序列，该模型的含义是，模型当前输出由当前输入和过去p 个时刻的输出决定；
步骤D：确定模型阶数，模型阶数是建模精度的重要决定因素，本步骤以建模绝对误差均值最小为指标搜索确定模型阶数：
Figure RE-GDA0002712173150000047
模型阶数搜索从1阶开始，若po+1阶模型的绝对误差均值MAE(po+1)大于MAE(p)，则po即为最优模型阶数。
进一步地，所述步骤三的具体内容为：利用时变自回归(TVAR)建模方法对周期跳跃信号进行建模，目标弹道周期跳跃信号是周期和幅值都随时间变化的非平稳信号，时变自回归模型用一组基函数的线性加权和来近似模型的时变参数，从而实现非平稳信号的精确建模，其中时变自回归(TVAR)建模方法如下：
步骤三一：时变自回归模型结构和基函数选择；
步骤三二：利用最小二乘法估计基函数展开系数；
步骤三三：确定模型阶数和基函数维数
进一步地，所述步骤三一的具体内容如下：
步骤E：自回归模型结构为：
Figure RE-GDA0002712173150000051
式中，y(n)是均值为零的非平稳信号序列，n＝1,…,N，p为模型阶数，cm(n)为模型的时变参数，e(n)为零均值、方差为
Figure RE-GDA0002712173150000052
的白噪声序列；
用一组基函数的线性加权和近似模型的时变参数bm(n)，即
Figure RE-GDA0002712173150000053
式中，L为基函数维数，cm,l为基函数的展开系数；
步骤F：把式(19)代入式(18)可得y(n)；
Figure RE-GDA0002712173150000054
进一步地，所述步骤三二是利用最小二乘法估计基函数展开系数，其具体内容如下：
定义矢量
Figure RE-GDA0002712173150000055
Figure RE-GDA0002712173150000056
利用
Figure RE-GDA0002712173150000057
可以把式(18)写成以下形式：
Figure RE-GDA0002712173150000058
式中C为基函数展开系数矢量：
C＝[c10,c11,…,c1L|…|cp0,cp1,…,cpL]T (23)
利用已有n时刻前的先验信号对y(n)进行预测：
Figure RE-GDA0002712173150000061
则预测误差e(n)和误差方差ζ可表示为：
Figure RE-GDA0002712173150000062
Figure RE-GDA0002712173150000063
根据最小二乘原理，最优基函数展开系数矢量C应该使预测误差e(n)的方差ζ最小，对式(26)求导并令其等于零，即可求得展开系数矢量C的估计值
Figure RE-GDA0002712173150000064
Figure RE-GDA0002712173150000065
进一步地，所述步骤三三，确定模型阶数和基函数维数，其具体内容如下：
步骤G：以建模均值误差最小为指标搜索模型阶数p和基函数维数L：
Figure RE-GDA0002712173150000066
式(28)中变量为正整数；模型阶数和基函数维数搜索流程如下：
步骤g1：给定基函数维数L*，利用信息论准则确定模型阶数的近似最优值p*，适用于TVAR模型的信息论准则有以下形式：
Figure RE-GDA0002712173150000067
式中，LR(p)是TVAR模型阶数的似然比，vp根据选择的信息准则有如下形式：
Figure RE-GDA0002712173150000068
式中：AIC为Akaike信息准则(Akaike information criterion)，MDL为最短描述长度准则(Minimum description length criterion)，MAP为最大后验概率准则 (Maximuma posteriori criterion)；
步骤g2：在区间[p*-k，p*+k]内搜索最优值popt，如果MAE(p+1,L*)不小于 MAE(p,L*)则认为p为最优模型阶数popt；
步骤g3：在区间[L*-s，L*+s]内搜索最优值Lopt，如果MAE(popt,L+1)不小于 MAE(popt,L+1)则认为L为最优基函数维数Lopt
进一步地，所述步骤四的具体内容为对步骤二建立的趋势信号模型和步骤三建立的周期跳跃信号模型进行叠合得到弹道完整参数化模型，然后基于完整模型外推实现轨迹预测，其叠合内容如下：
目标弹道完整模型通过趋势信号RD(n)的AR模型和跳跃信号T(n)的TVAR模型叠加得到：
Figure RE-GDA0002712173150000071
式中：n＝1,…,N；
利用目标弹道完整模型(31)实现j＝N+1,…,Np时刻的弹道预测，得到预测轨迹为：
Figure RE-GDA0002712173150000072
本发明与现有技术相比具有以下有益效果：
1、本发明提供的一种临近空间高超声速目标轨迹预测方法把可观测弹道数据作为建模对象，不依赖目标动力学模型等先验信息和制导律估计，适用于各种临近空间高超声速目标轨迹预测且具有很高的预测精度。
2、本发明提供的一种临近空间高超声速目标轨迹预测方法根据目标弹道特点，采用集成经验模态分解方法把弹道分解为趋势信号和周期跳跃信号分别建模，提高了建模和预测精度。
3、本发明提供的一种临近空间高超声速目标轨迹预测方法考虑了周期跳跃信号的非平稳特性，使用了时变自回归方法进行周期跳跃信号建模，有效解决了建模中存在的周期漂移问题，提高了周期跳跃信号的建模和预测精度。
4、本发明提供的一种临近空间高超声速目标轨迹预测方法把建模绝对误差均值作为指标，设计了一种时变自回归模型阶数和基函数维数搜索方法，提高了建模和预测精度。
附图说明
图1Burg算法估计AR模型参数流程图；
图2TVAR模型阶数和基函数维数搜索流程图；
图3地心距量测信号时域图；
图4地心距趋势信号时域图；
图5地心距周期跳跃信号时域图；
图6地心距趋势信号建模结果；
图7地心距周期跳跃信号建模结果；
图8地心距预测结果；
图9地心距预测误差。
具体实施方式
具体实施方式一：，本实施方式提供了一种临近空间高超声速目标轨迹预测方法，所述步骤是通过如下步骤实现的：
步骤一：将所得历史弹道数据分解为弹道趋势信号和弹道周期跳跃信号；
步骤二：对趋势信号进行建模；
步骤三：对周期跳跃信号进行建模；
步骤四：对步骤二建立的趋势信号模型和步骤三建立的周期跳跃信号模型进行叠合得到弹道完整参数化模型，然后基于完整模型外推实现轨迹预测。
本实施方式中，步骤一中，采用集成经验模态分解方法将所得历史弹道数据分解为弹道趋势信号和弹道周期跳跃信号。高超声速目标弹道数据是非线性、非平稳信号，具有明显的趋势和周期跳跃的特点，为了提高参数化建模精度，本发明采用经验模态分解方法，把历史弹道数据分解为弹道趋势信号和弹道周期跳跃信号以便于对不同模态信号进行精确建模；
步骤二中，利用自回归(AR)建模方法对趋势信号进行建模。弹道趋势信号为线性、慢变、低频的平稳信号，采用自回归(AR)建模方法进行平稳信号参数化建模时，模型阶数和模型系数决定了建模精度。因此，步骤2主要包括：以建模绝对误差均值最小为指标搜索自回归模型结构和采用Burg算法估计自回归模型系数两部分；
步骤三中，利用时变自回归(TVAR)建模方法对周期跳跃信号进行建模。目标弹道周期跳跃信号是周期和幅值都随时间变化的非平稳信号，时变自回归模型用一组基函数的线性加权和来近似模型的时变参数，从而实现非平稳信号的精确建模。时变自回归 (TVAR)建模方法建模过程主要包括三个方面的内容：1)基函数选择；2)利用最小二乘法估计基函数展开系数；3)采用建模绝对误差均值最小为指标，基于信息论准改进模型阶数和基函数维数搜索方法，确定最优模型阶数和基函数维数，以减小建模误差；
步骤四中，对步骤二建立的趋势信号模型和步骤三建立的周期跳跃信号模型进行叠合得到弹道完整参数化模型，然后基于完整模型外推实现轨迹预测。
具体实施方式二：本实施方式是对具体实施方式一所述的步骤一作进一步限定，本实施方式中，所述步骤一的具体内容为：采用集成经验模态分解方法将所得历史弹道数据分解为弹道趋势信号和弹道周期跳跃信号，其中集成经验模态分解方法如下：
步骤一一：在弹道序列X(n)中分别加入K组高斯白噪声序列Gk，得到加噪信号Xk(n)，即：
Xk(n)＝X(n)+Gk (1)
式中，k＝1,…,K
步骤一二：对每组加噪信号Xk(n)进行EMD分解，得到I个固有模态函数IMF分量和残余项：
Figure RE-GDA0002712173150000091
步骤一三：把K次EMD分解的输出均值作为集成经验模态分解的分解结果：
Figure RE-GDA0002712173150000092
步骤一四：原始弹道数据可表示为：
X(n)＝T(n)+R(n) (4)
Figure RE-GDA0002712173150000093
式中，残余项R(n)即为弹道数据的趋势信号，T(n)即为弹道数据的非平稳周期跳跃信号。其它组成及连接方式与具体实施方式一相同。
本实施方式中，根据目标弹道特点，采用集成经验模态分解方法把弹道分解为趋势信号和周期跳跃信号，并利用两种信号分别建模，提高了建模精度和预测精度。
具体实施方式三：本实施方式是对具体实施方式一中所述的步骤二作进一步限定，本实施方式中，所述步骤二的具体内容为：利用自回归(AR)建模方法对趋势信号进行建模，模型阶数和模型系数决定了建模精度，其中自回归(AR)建模方法如下：
步骤二一：自回归模型系数估计；
步骤二二：自回归模型阶数选择。
其它组成及连接方式与具体实施方式二相同。
本实施方式中，趋势信号RD(n)是线性、慢变、低频信号，周期和幅值都随时间变化的平稳信号，因此采用自回归(AR)建模方法足以准确的保证建模的稳定和准确性。
具体实施方式四：参照图1说明本实施方式，本实施方式是对具体实施方式三所述的步骤二一作进一步限定，本实施方式中，所述步骤二一的具体内容如下：
步骤A：令向前、向后预测误差功率之和为最小：
Figure RE-GDA0002712173150000101
式中，上标f表示前向预测(Forward prediction)，上标b表示后向预测(Backwardprediction)，上式中
Figure RE-GDA0002712173150000102
式中：上标f表示前向预测(Forward prediction)，上标b表示后向预测(Backwardprediction)，
Figure RE-GDA0002712173150000103
为前向预测误差，
Figure RE-GDA0002712173150000104
为后向预测误差，ρf为前向预测误差功率，ρb为后向预测误差功率；
则向前、向后观测误差分别为：
Figure RE-GDA0002712173150000105
上式中为了保证不超出给定数据范围，求和范围为p≤n≤N-1；
反射系数kp的递推公式形式如下：
Figure RE-GDA0002712173150000111
式中：rxx(p)为自相关函数，
Figure RE-GDA0002712173150000112
为预测误差，
使预测误差功率之和为最小时的反射系数kp，令：
Figure RE-GDA0002712173150000113
将式(8)带人式(12)中，可得：
Figure RE-GDA0002712173150000114
步骤B：根据步骤A的推导过程，利用Burg递推算法求AR模型参数：
步骤b1：给定初始条件：
Figure RE-GDA0002712173150000115
步骤b2：利用式(13)计算kp：
步骤b3：更新预测误差平均功率：
ρp＝(1-|kp|2)ρp-1 (13)
步骤b4：更新模型系数：
Figure RE-GDA0002712173150000116
步骤b5：更新向前、向后观测误差：
Figure RE-GDA0002712173150000117
步骤b6：令p＝p+1，返回步骤b2，直到所需模型阶数为止。其它组成及连接方式与具体实施方式三相同。
具体实施方式五：参照图1说明本实施方式，本实施方式是对具体实施方式三所述的步骤二二作进一步限定，本实施方式中，步骤C：确实模型形式，趋势信号RD(n)是线性、慢变、低频信号，因此自回归模型的形式如下：
Figure RE-GDA0002712173150000121
式中，x(n)是均值为零的信号序列，n＝1,…,N，p为模型阶数，ak为模型系数，u(n)为零均值、方差为
Figure RE-GDA0002712173150000122
的白噪声序列，该模型的含义是，模型当前输出由当前输入和过去p 个时刻的输出决定；
步骤D：确定模型阶数，模型阶数是建模精度的重要决定因素，本步骤以建模绝对误差均值最小为指标搜索确定模型阶数：
Figure RE-GDA0002712173150000123
模型阶数搜索从1阶开始，若po+1阶模型的绝对误差均值MAE(po+1)大于MAE(p)，则po即为最优模型阶数。
其它组成及连接方式与具体实施方式四相同。
具体实施方式六：参照图2说明本实施方式，本实施方式是对具体实施方式一所述的步骤三作进一步限定，本实施方式中，所述步骤三的具体内容为：利用时变自回归(TVAR)建模方法对周期跳跃信号进行建模，目标弹道周期跳跃信号是周期和幅值都随时间变化的非平稳信号，时变自回归模型用一组基函数的线性加权和来近似模型的时变参数，从而实现非平稳信号的精确建模，其中时变自回归(TVAR)建模方法如下：
步骤三一：时变自回归基函数选择；
步骤三二：利用最小二乘法估计基函数展开系数；
步骤三三：确定模型阶数和基函数维数。
其它组成及连接方式与具体实施方式五相同。
本实施方式中，目标弹道数据周期跳跃信号T(n)是周期和幅值都随时间变化的非平稳信号，TVAR模型是非平稳信号建模的有效方法之一。TVAR模型与AR模型相比最大特征是其模型参数具有时变性，常用于TVAR模型的基函数有Legendre基函数、沃尔什基基函数、傅里叶基函数、DCT基函数、小波基函数等。根据信号特点选择不同的基函数， Legendre基函数适用于缓变的时间序列建模，傅立叶基函数和DCT基函数适用于准周期变化信号。
具体实施方式七：参照图2说明本实施方式，本实施方式是对具体实施方式六所述的步骤三一作进一步限定，本实施方式中所述步骤三一的具体内容如下：
步骤E：自回归模型结构为：
Figure RE-GDA0002712173150000131
式中，y(n)是均值为零的非平稳信号序列，n＝1,…,N，p为模型阶数，cm(n)为模型的时变参数，e(n)为零均值、方差为
Figure RE-GDA0002712173150000132
的白噪声序列；
用一组基函数的线性加权和近似模型的时变参数bm(n)，即
Figure RE-GDA0002712173150000133
式中，L为基函数维数，cm,l为基函数的展开系数；
步骤F：把式(19)代入式(18)可得y(n)；
Figure RE-GDA0002712173150000134
其它组成及连接方式与具体实施方式六相同。
具体实施方式八：参照图2说明本实施方式，本实施方式是对具体实施方式六所述的步骤三二作进一步限定，本实施方式中，所述步骤三二是利用最小二乘法估计基函数展开系数，其具体内容如下：
定义矢量
Figure RE-GDA0002712173150000135
Figure RE-GDA0002712173150000136
利用
Figure RE-GDA0002712173150000137
可以把式(18)写成以下形式：
Figure RE-GDA0002712173150000138
式中C为基函数展开系数矢量：
C＝[c10,c11,…,c1L|…|cp0,cp1,…,cpL]T (23)
利用已有n时刻前的先验信号对y(n)进行预测：
Figure RE-GDA0002712173150000139
则预测误差e(n)和误差方差ζ可表示为：
Figure RE-GDA00027121731500001310
Figure RE-GDA0002712173150000141
根据最小二乘原理，最优基函数展开系数矢量C应该使预测误差e(n)的方差ζ最小，对式(26)求导并令其等于零，即可求得展开系数矢量C的估计值
Figure RE-GDA0002712173150000142
Figure RE-GDA0002712173150000143
其它组成及连接方式与具体实施方式七相同。
具体实施方式九：参照图2说明本实施方式，本实施方式是对具体实施方式六所述的步骤三三作进一步限定，本实施方式中，所述步骤三三，确定模型阶数和基函数维数，其具体内容如下：
步骤G：以建模均值误差最小为指标搜索模型阶数p和基函数维数L：
Figure RE-GDA0002712173150000144
式(28)中变量为正整数；模型阶数和基函数维数搜索流程如下：
步骤g1：给定基函数维数L*，利用信息论准则确定模型阶数的近似最优值p*，适用于TVAR模型的信息论准则有以下形式：
Figure RE-GDA0002712173150000145
式中，LR(p)是TVAR模型阶数的似然比，vp根据选择的信息准则有如下形式：
Figure RE-GDA0002712173150000146
式中：AIC为Akaike信息准则(Akaike information criterion)，MDL为最短描述长度准则(Minimum description length criterion)，MAP为最大后验概率准则
(Maximum a posteriori criterion)；
步骤g2：在区间[p*-k，p*+k]内搜索最优值popt，如果MAE(p+1,L*)不小于 MAE(p,L*)则认为p为最优模型阶数popt；
步骤g3：在区间[L*-s，L*+s]内搜索最优值Lopt，如果MAE(popt,L+1)不小于MAE(popt,L+1)则认为L为最优基函数维数Lopt。
其它组成及连接方式与具体实施方式八相同。
具体实施方式十：本实施方式是对具体实施方式一所述的步骤四作进一步限定，本实施方式中，所述步骤四的具体内容为对步骤二建立的趋势信号模型和步骤三建立的周期跳跃信号模型进行叠合得到弹道完整参数化模型，然后基于完整模型外推实现轨迹预测，其叠合内容如下：
目标弹道完整模型通过趋势信号RD(n)的AR模型和跳跃信号T(n)的TVAR模型叠加得到：
Figure RE-GDA0002712173150000151
式中：n＝1,…,N；
利用目标弹道完整模型(31)实现j＝N+1,…,Np时刻的弹道预测，得到预测轨迹为：
Figure RE-GDA0002712173150000152
实施例
本发明实施例仿真采用的计算机配置为：CPU为i3-2100，主频3.1GHz,内存3GB。仿真基于远程高超声速滑翔式再入飞行器为模型，飞行器参数采用美国通用航空飞行器(Common Aero Vehicle,CAV)的设计参数，飞行器质量907.2kg，特征面积0.4837m2，初始高度60km，速度6km/s，当地速度倾角0°，飞行时长700s，采样周期0.1s，以最大升阻比飞行弹道的地心距为例(图3)进行仿真验证。
以690s作为预测起点，对弹道数据进行加噪声处理，模拟生成目标地心距历史量测数据。首先进行步骤一，采用集成经验模态分解方法，提取地心距趋势信号(见图4)和地心距周期信号(见图5)。分解得到的趋势信号和周期跳跃信号表现出不同的形态和规律，使用不同建模方法构建信号参数化模型，以提高建模精度。
执行步骤二，由图4可见，地心距趋势信号是线性、平稳、低频信号。利用自回归(AR)建模方法对趋势信号进行建模。模型系数搜索以建模绝对误差均值最小为指标，在每步搜索过程中采用Burg方法估计模型系数。二阶模型建模绝对误差均值小于三阶模型建模绝对误差均值，因此最优模型阶数为二阶，得到地心距趋势信号参数化模型为：
x(n)＝-1.99x(n-1)+x(n-2) (33)
执行步骤三。由图5可见，地心距周期跳跃信号是幅值和周期随时间变化的非平稳、高频信号。利用时变自回归(TVAR)建模方法对地心距周期跳跃信号进行建模，建模结果见图7。首先根据信号特征选择傅里叶基函数拟合模型自回归系数：
Figure RE-GDA0002712173150000161
式中：k＝0,1,2,…。
然后，执行图2所示的最优模型阶数和基函数维数搜索流程。在每步搜索过程中，使用3.2介绍的基于最小二乘法基函数展开系数估计方法进行基函数展开系数估计。首先，令基函数维数为L*等于12，采用AIC准则确定的近似模型阶数p*为三阶。然后，在区间[1，5]内搜索最优模型阶数。二阶模型建模绝对误差均值小于三阶模型建模绝对误差均值，因此最优模型阶数为二阶。最后，在区间[6，18]搜索最优基函数维数，确定最优基函数维数为9维。因此，得到地心距周期跳跃信号参数化模型：
Figure RE-GDA0002712173150000162
式中：展开系数估计值
Figure RE-GDA0002712173150000163
Figure RE-GDA0002712173150000164
由见图6和图7建模结果可知，采用自回归方法对趋势信号建模精度很高，而对周期跳跃信号建模存在模型周期漂移问题，导致建模误差很大。而采用时变自回归方法对周期跳跃信号建模能有效的解决上述问题，且建模精度极高。此外，相对于周期跳跃信号建模误差，趋势信号建模误差极小。因此周期跳跃信号建模精度直接决定了弹道数据建模和轨迹预测精度。
执行步骤四。地心距趋势信号参数化模型时(33)和地心距周期跳跃信号参数化模型(35)，得到地心距信号完整参数化模型：
Figure RE-GDA0002712173150000165
通过对模型(36)外推实现地心距轨迹预测，预测结果与预测误差分布见图8和图9。仿真中把不同模态数据均使用自回归(AR)建模方法进行建模预测作为对比项。
如图9所示，本发明的地心距预测误差在1.56km以内，预测绝对误差均值为0.5882km，而自回归模型预测绝对误差均值为1.3959km。仿真结果表明本发明算法能精确的提炼和拟合目标机动规律，实现了对高超声速目标机动轨迹的精确预报，验证了本发明的临近空间高超声速目标轨迹预测方法的优越性。
由上可知，本发明实施例实现了不基于先验信息的目标纵向地心距机动轨迹精确预测。但以上所述仅为本发明的较佳实施例而已，并不用于限制本发明，本发明同样适用于高超声速目标横向机动轨迹预测、目标速度预测等方面。凡在本发明的原则和精神之内所作的任何修改、等同替换和改进等，均就包含在本发明的保护范围之内。

Claims (10)

1.一种临近空间高超声速目标轨迹预测方法，其特征在于：所述步骤是通过如下步骤实现的：
步骤一：将所得历史弹道数据分解为弹道趋势信号和弹道周期跳跃信号；
步骤二：对趋势信号进行建模；
步骤三：对周期跳跃信号进行建模；
步骤四：对步骤二建立的趋势信号模型和步骤三建立的周期跳跃信号模型进行叠合得到弹道完整参数化模型，然后基于完整模型外推实现轨迹预测。
2.根据权利要求1中所述的一种临近空间高超声速目标轨迹预测方法，其特征在于：所述步骤一的具体内容为：采用集成经验模态分解方法将所得历史弹道数据分解为弹道趋势信号和弹道周期跳跃信号，其中集成经验模态分解方法如下：
步骤一一：在弹道序列X(n)中分别加入K组高斯白噪声序列Gk，得到加噪信号Xk(n)，即：
Xk(n)＝X(n)+Gk (1)
式中，k＝1,…,K
步骤一二：对每组加噪信号Xk(n)进行EMD分解，得到I个固有模态函数IMF分量和残余项：
Figure RE-FDA0002712173140000011
步骤一三：把K次EMD分解的输出均值作为集成经验模态分解的分解结果：
Figure RE-FDA0002712173140000012
步骤一四：原始弹道数据可表示为：
X(n)＝T(n)+R(n) (4)
Figure RE-FDA0002712173140000013
式中，残余项R(n)即为弹道数据的趋势信号，T(n)即为弹道数据的非平稳周期跳跃信号。
3.根据权利要求2中所述的一种临近空间高超声速目标轨迹预测方法，其特征在于：所述步骤二的具体内容为：利用自回归AR建模方法对趋势信号进行建模，模型阶数和模型系数决定了建模精度，其中自回归AR建模方法如下：
步骤二一：自回归模型系数估计；
步骤二二：自回归模型阶数选择。
4.根据权利要求3中所述的一种临近空间高超声速目标轨迹预测方法，其特征在于：所述步骤二一的具体内容如下：
步骤A：令向前、向后预测误差功率之和为最小；
Figure RE-FDA0002712173140000021
式中，上标f表示前向预测，上标b表示后向预测，上式中
Figure RE-FDA0002712173140000022
式中：上标f表示前向预测，上标b表示后向预测，
Figure RE-FDA0002712173140000023
为前向预测误差，
Figure RE-FDA0002712173140000024
为后向预测误差，ρf为前向预测误差功率，ρb为后向预测误差功率；
则向前、向后观测误差分别为：
Figure RE-FDA0002712173140000025
上式中为了保证不超出给定数据范围，求和范围为p≤n≤N-1；
反射系数kp的递推公式形式如下：
Figure RE-FDA0002712173140000026
式中：rxx(p)为自相关函数，
Figure RE-FDA0002712173140000027
为预测误差，
确定预测误差功率之和为最小时的反射系数kp，令：
Figure RE-FDA0002712173140000031
将式(8)带人式(12)中，可得：
Figure RE-FDA0002712173140000032
步骤B：根据步骤A的推导过程，利用Burg递推算法求AR模型参数：
步骤b1：给定初始条件：
Figure RE-FDA0002712173140000033
步骤b2：利用式(10)计算kp：
步骤b3：更新预测误差平均功率：
ρp＝(1-|kp|2)ρp-1 (13)
步骤b4：更新模型系数：
Figure RE-FDA0002712173140000034
步骤b5：更新向前、向后观测误差：
Figure RE-FDA0002712173140000035
步骤b6：令p＝p+1，返回步骤b2，直到所需模型阶数为止。
5.根据权利要求4中所述的一种临近空间高超声速目标轨迹预测方法，其特征在于：所述步骤二二中的模型系数直接关系着建模和预测精度，本步骤基于Burg算法确定自回归模型系数，具体过程为：
步骤C：确实模型形式：趋势信号RD(n)是线性、慢变、低频信号，因此自回归模型的形式如下：
Figure RE-FDA0002712173140000036
式中，x(n)是均值为零的信号序列，n＝1,…,N，p为模型阶数，ak为模型系数，u(n)为零均值、方差为
Figure RE-FDA0002712173140000037
的白噪声序列，该模型的含义是，模型当前输出由当前输入和过去p个时刻的输出决定；
步骤D：确定模型阶数：模型阶数是建模精度的重要决定因素，本步骤以建模绝对误差均值最小为指标搜索确定模型阶数：
Figure RE-FDA0002712173140000041
模型阶数搜索从1阶开始，若po+1阶模型的绝对误差均值MAE(po+1)大于MAE(p)，则po即为最优模型阶数。
6.根据权利要求5中所述的一种临近空间高超声速目标轨迹预测方法，其特征在于：所述步骤三的具体内容为：利用时变自回归TVAR建模方法对周期跳跃信号进行建模，目标弹道周期跳跃信号是周期和幅值都随时间变化的非平稳信号，时变自回归模型用一组基函数的线性加权和来近似模型的时变参数，从而实现非平稳信号的精确建模，其中时变自回归TVAR建模方法如下：
步骤三一：时变自回归基函数选择；
步骤三二：利用最小二乘法估计基函数展开系数；
步骤三三：确定模型阶数和基函数维数。
7.根据权利要求6中所述的一种临近空间高超声速目标轨迹预测方法，其特征在于：所述步骤三一的具体内容如下：
步骤E：自回归模型结构为：
Figure RE-FDA0002712173140000042
式中，y(n)是均值为零的非平稳信号序列，n＝1,…,N，p为模型阶数，cm(n)为模型的时变参数，e(n)为零均值、方差为
Figure RE-FDA0002712173140000043
的白噪声序列；
用一组基函数的线性加权和近似模型的时变参数bm(n)，即
Figure RE-FDA0002712173140000044
式中，L为基函数维数，cm,l为基函数的展开系数；
步骤F：把式(19)代入式(18)可得y(n)；
Figure RE-FDA0002712173140000045
8.根据权利要求7中所述的一种临近空间高超声速目标轨迹预测方法，其特征在于：所述步骤三二是利用最小二乘法估计基函数展开系数，其具体内容如下：
定义矢量
Figure RE-FDA0002712173140000051
Figure RE-FDA0002712173140000052
利用
Figure RE-FDA0002712173140000053
可以把式(18)写成以下形式：
Figure RE-FDA0002712173140000054
式中C为基函数展开系数矢量：
C＝[c10,c11,…,c1L|…|cp0,cp1,…,cpL]T (23)
利用已有n时刻前的先验信号对y(n)进行预测：
Figure RE-FDA0002712173140000055
则预测误差e(n)和误差方差ζ可表示为：
Figure RE-FDA0002712173140000056
Figure RE-FDA0002712173140000057
根据最小二乘原理，最优基函数展开系数矢量C应该使预测误差e(n)的方差ζ最小，对式(26)求导并令其等于零，即可求得展开系数矢量C的估计值
Figure RE-FDA0002712173140000058
Figure RE-FDA0002712173140000059
9.根据权利要求8中所述的一种临近空间高超声速目标轨迹预测方法，其特征在于：所述步骤三三，确定模型阶数和基函数维数，其具体内容如下：
步骤G：以建模均值误差最小为指标搜索模型阶数p和基函数维数L：
Figure RE-FDA00027121731400000510
式(28)中变量为正整数；模型阶数和基函数维数搜索流程如下：
步骤g1：给定基函数维数L*，利用信息论准则确定模型阶数的近似最优值p*，适用于TVAR模型的信息论准则有以下形式：
Figure RE-FDA00027121731400000511
式中，LR(p)是TVAR模型阶数的似然比，vp根据选择的信息准则有如下形式：
Figure RE-FDA0002712173140000061
式中：AIC为Akaike信息准则(Akaike information criterion)，MDL为最短描述长度准则(Minimum description length criterion)，MAP为最大后验概率准则(Maximum aposteriori criterion)；
步骤g2：在区间[p*-k，p*+k]内搜索最优值popt，如果MAE(p+1,L*)不小于MAE(p,L*)则认为p为最优模型阶数popt；
步骤g3：在区间[L*-s，L*+s]内搜索最优值Lopt，如果MAE(popt,L+1)不小于MAE(popt,L+1)则认为L为最优基函数维数Lopt。
10.根据权利要求9中所述的一种临近空间高超声速目标轨迹预测方法，其特征在于：所述步骤四的具体内容为对步骤二建立的趋势信号模型和步骤三建立的周期跳跃信号模型进行叠合得到弹道完整参数化模型，然后基于完整模型外推实现轨迹预测，其叠合内容如下：
目标弹道完整模型通过趋势信号RD(n)的AR模型和跳跃信号T(n)的TVAR模型叠加得到：
Figure RE-FDA0002712173140000062
式中：n＝1,…,N；
利用目标弹道完整模型(31)实现j＝N+1,…,Np时刻的弹道预测，得到预测轨迹为：
Figure RE-FDA0002712173140000063
CN202010640882.3A 2020-07-06 2020-07-06 一种临近空间高超声速目标轨迹预测方法 Active CN111931287B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010640882.3A CN111931287B (zh) 2020-07-06 2020-07-06 一种临近空间高超声速目标轨迹预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010640882.3A CN111931287B (zh) 2020-07-06 2020-07-06 一种临近空间高超声速目标轨迹预测方法

Publications (2)

Publication Number Publication Date
CN111931287A true CN111931287A (zh) 2020-11-13
CN111931287B CN111931287B (zh) 2023-02-24

Family

ID=73314051

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010640882.3A Active CN111931287B (zh) 2020-07-06 2020-07-06 一种临近空间高超声速目标轨迹预测方法

Country Status (1)

Country Link
CN (1) CN111931287B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112861262A (zh) * 2021-02-08 2021-05-28 北京理工大学 用于高动态载体的弹道预测的自组织数据驱动建模方法
CN113269363A (zh) * 2021-05-31 2021-08-17 西安交通大学 一种高超声速飞行器的轨迹预测方法、系统、设备及介质

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110309909A (zh) * 2019-06-26 2019-10-08 北京控制工程研究所 一种高速大范围机动目标轨迹的智能实时预测方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110309909A (zh) * 2019-06-26 2019-10-08 北京控制工程研究所 一种高速大范围机动目标轨迹的智能实时预测方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
HAN C等: "Method of trajectory prediction for unpowered gliding hypersonic vehicle in gliding phase", 《2016 IEEE ADVANCED INFORMATION MANAGEMENT, COMMUNICATES, ELECTRONIC AND AUTOMATION CONTROL CONFERENCE (IMCEC)》 *
KAI Z等: "Multi-Level Recursive Trajectory Prediction for Hypersonic Gliding Reentry Vehicle", 《MODERN DEFENSE TECHNOLOGY》 *
张凯等: "基于意图推断的高超声速滑翔目标贝叶斯轨迹预测", 《宇航学报》 *
张召等: "谱估计理论在弹道数据参数化建模中的应用", 《宇航学报》 *
李庚泽等: "基于轨迹预测的高超声速飞行器拦截中/末制导研究", 《上海航天》 *
韩春耀等: "高超声速滑翔飞行器轨迹预测分析", 《火力与指挥控制》 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112861262A (zh) * 2021-02-08 2021-05-28 北京理工大学 用于高动态载体的弹道预测的自组织数据驱动建模方法
CN113269363A (zh) * 2021-05-31 2021-08-17 西安交通大学 一种高超声速飞行器的轨迹预测方法、系统、设备及介质
CN113269363B (zh) * 2021-05-31 2023-06-30 西安交通大学 一种高超声速飞行器的轨迹预测方法、系统、设备及介质

Also Published As

Publication number Publication date
CN111931287B (zh) 2023-02-24

Similar Documents

Publication Publication Date Title
Julier et al. Unscented filtering and nonlinear estimation
Gordon et al. Bayesian state estimation for tracking and guidance using the bootstrap filter
CN111722214B (zh) 雷达多目标跟踪phd实现方法
CN111983927B (zh) 一种最大协熵mcc准则的椭球集员滤波方法
CN111931287B (zh) 一种临近空间高超声速目标轨迹预测方法
CN111783307A (zh) 一种高超声速飞行器状态估计方法
CN111474960A (zh) 基于控制量特征辅助的临空高超声速目标跟踪方法
Liu et al. Modified unscented Kalman filter using modified filter gain and variance scale factor for highly maneuvering target tracking
Yang et al. UGHF for acoustic tracking with state-dependent propagation delay
Li et al. Auxiliary truncated particle filtering with least-square method for bearings-only maneuvering target tracking
Zuo et al. Simplified unscented particle filter for nonlinear/non-Gaussian Bayesian estimation
RU2570111C1 (ru) Устройство радиолокационного распознавания воздушно-космических объектов
CN116047495B (zh) 一种用于三坐标雷达的状态变换融合滤波跟踪方法
Chen et al. IMM tracking of a 3D maneuvering target with passive TDOA system
Havangi An adaptive particle filter based on PSO and fuzzy inference system for nonlinear state systems
Han et al. Maneuvering target tracking with unknown acceleration using retrospective-cost-based adaptive input and state estimation
Singh et al. Tracking of ballistic target on re-entry using ensemble Kalman filter
Li et al. Motion modeling and state estimation in range-squared coordinate
Peng et al. An IMM-VB algorithm for hypersonic vehicle tracking with heavy tailed measurement noise
Yadaiah et al. Fuzzy Kalman Filter based trajectory estmation
Luo et al. Joint estimation of target location and relative altitude from angle measurements
Nie et al. Adaptive tracking algorithm based on 3D variable turn model
Nanda et al. Performance analysis of Cubature rule based Kalman filter for target tracking
CN112819303B (zh) 基于pce代理模型的飞行器追踪效能评估方法及系统
Agate et al. Particle filtering algorithm for tracking multiple road-constrained targets

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