CN111783358B - 一种基于贝叶斯估计的高超速飞行器长期轨迹预报方法 - Google Patents

Abstract

本发明涉及轨迹预报技术，是一种基于贝叶斯估计的高超速飞行器长期轨迹预报方法，为解决高超声速飞行器轨迹预报的主要技术难点：构建合适的目标意图代价函数，准确反映目标飞行意图；开发目标飞行的可行区快速结算方法；基于贝叶斯估计设计合理的多步递推预测算法，预测目标攻击意图和轨迹预测，实时建立预测轨迹簇；本发明主要包括四部分内容：多指标意图函数建立、高超声速飞行器可行区分析、贝叶斯估计下意图推测算法和基于蒙特卡洛打靶的高超声速飞行器轨迹预报方法一方面推测目标飞行意图，降低目标机动致使的轨迹不确定性；另一方面设计实时未来轨迹多步递推算法，充分挖掘高超声速飞行器运动的潜在规律，实现目标轨迹的长时间预报。

Description

一种基于贝叶斯估计的高超速飞行器长期轨迹预报方法

技术领域

本发明涉及轨迹预报技术，是一种基于贝叶斯估计的高超速飞行器长期轨迹预报方法。

背景技术

高超声速飞行器是一类飞行马赫数大于5，依靠气动升力控制在临近空间内长时间远距离飞行的飞行器。从防御方角度而言，针对高速大机动目标通常采用预测命中点法进行制导拦截。在拦截器接近目标时，以相反的速度矢量正面迎击目标，使得拦截器能够以小的速度、过载能力完成拦截任务。这种方式对目标预测精度要求高，目前成熟的轨迹预测方法多是针对弹道导弹和空间再入目标，采用的预测方法以其轨迹的相对固定性为理论依据。高超声速目标利用其高升阻比外形机动能形成一个较大的打击区域，与传统的惯性目标相比，具有强机动能力、无固定飞行轨迹，一般性方法难以对其弹道预测。此外，考虑目标飞行器可能在探测到拦截目标后发生机动、亦或是未来某时刻改变进攻目标，为了增加防御成功率，防御方需要合理分析当前空中目标态势并进行威胁评估，推测目标飞行意图，以降低目标未来机动不确定性的影响。据此，一种长期的高超声速飞行器轨迹预测方法将是防御方研究的重点。

贝叶斯估计理论是使用输入测量和数学过程模型随时间递归估计未知概率密度函数的一般概率方法。在目标预测中，从量测数据提取对意图推断有用信息，结合已知先验，在贝叶斯框架下完成对空中目标态势与威胁评估进而预测轨迹。这种方法实质是根据飞行意图推断未来时刻目标在特定空域概率，已在机器人行为预测、民航轨迹监测和飞航导弹拦截等领域上广泛应用，具有可行性。

基于意图推断的轨迹预测算法从原理上可以分为两大部分：意图估计与轨迹预测。意图估计部分完成的是对目标行为模式的判断，因此设计一个合理的意图代价函数尤为重要，决定了预测准确性。轨迹预测部分完成的是对未来时刻空间概率的估计，需要多步递推未来状态，因而确定好未来可行空间是预测的先决条件。同时，基于贝叶斯递推估计未来状态时，会存在随预测步数增加，计算复杂度呈指数增大的情况，不利于实时预测。

所以，高超声速飞行器轨迹预报存在以下主要技术难点：

1、构建合适的意图代价函数，准确反映目标飞行意图；

2、快速、正确解算目标飞行的可行区；

3、设计合理的目标意图推测算法和目标轨迹预测算法。多步递推预测算法，基于贝叶斯估计能实时建立预测轨迹簇。

发明内容

发明目的：为了克服高超声速飞行器轨迹预报存在以下主要技术难点，本发明提供一种基于贝叶斯估计的高超速飞行器长期轨迹预报方法。

本发明是通过以下方案实施的：一种基于贝叶斯估计的高超速飞行器长期轨迹预报方法，它包括以下步骤：

步骤一：建立意图代价函数；

在水平面内分析轨迹预测问题，在飞行区域中，存在若干攻击目标，假设高超飞行器作战意图主要体现在航向误差、弹目距离以及战场环境上，并结合航向误差、弹目距离变化以及战场环境定义飞行器意图代价函数；

步骤二：进行可行域分析；

(1)以“区域”的概念将飞行环境离散，将对位置参数的概率估计转化为对区域的概率估计；

(2)采用简化动力学模型快速预测所有可行区间；

步骤三：贝叶斯估计框架下，进行目标意图推测；

贝叶斯估计框架下，基于已有的观测信息确定攻击目标的概率，进行意图推测；

步骤四：基于贝叶斯估计和蒙特卡洛打靶法进行轨迹预测；

所述步骤一具体为：

在飞行区域ψ中，存在若干攻击目标θ，定义飞行器意图代价函数为：

归一化系数K计算方式如下：

式中θ_η∈Θ表示某个已知的攻击目标，φ_k∈Φ表示战场区域内某一禁飞区，χ⁺表示可行集”；

公式(0.1)中右侧第一项a·(δ(x_i,X_i+1,θ_η)-δ(x_i,θ_η))表示弹目距离变化，用δ(·)表示括号内元素最短距离，有以下关系：

δ(a,b,c)＝δ(a,b)+δ(b,c) (0.3)

式中，a表示弹目距离的权重系数、b表示航向角的权重系数和c表示禁飞区的权重系数；

该项值恒为非负数，当飞行器沿最短路径接近目标时，为极限值0；权重系数a表征了飞行器以最小路径前往目标的潜在意图，当权重系数a→+∞，朝着最短路径移动会获得最大条件概率；当权重系数a→0，朝着各个方向移动条件概率相当；

公式(0.1)中右侧第二项b·|ν(x_i,X_i+1)-ν(x_i,θ_η)|/δ(x_i,θ_η)表示航向角偏差，ν(·)表示括号内元素连线与正北方向夹角；用Δν(·)表示该项分子，其值为0时，代表飞行器速度方向指向了该目标，这一移动会获得最大的条件概率；同时，考虑到当飞行器离目标较远时，运动方向更自由，在该项分母引入了弹目距离δ(x_i,θ_η)，使其获得自适应的权重系数；越接近目标，与意图函数中其他项相比，航向角偏差项权重更大，

公式(0.1)中右侧第三项c·∑I_φ(x_i,X_i+1,φ_k)表示禁飞区约束影响，定义为：

为第k个禁飞区半径,

、

为当前速度方向与禁飞区切线夹角，当速度矢量在切线角内时，夹角值定义为正，若飞行器朝向禁飞区飞行，

越大，代价也越大，会获得较小条件概率，

为安全裕度，当飞行器离禁飞区较远时，运动自由，安全裕度较大，任意朝向均不受禁飞区约束影响；靠近禁飞区时，安全裕度减小，只有避开禁飞区方向才能获得较大概率，下面给出

值计算方式；

若瞬时速度方向指向禁飞区中心，则

；当飞行器恰好可以极限转弯半径r_p从禁飞区边缘飞过，定义此时速度方向与切线夹角为

根据三角函数关系可知：

高超飞行器侧向机动能力主要由升力分量提供，因此转弯半径可由下式估算：

该式在应用时，状态参数取跟踪时刻末值，攻角α、倾侧角σ取为极限值，以估计极限转弯半径r_p，m为飞行器质量，v为飞行器速度，L为升力，C_L为飞行器升力系数，ρ为当前大气密度，S为飞行器特征面积，σ为飞行器倾侧角，

攻角α与飞行器升力系数C_L存在关系C_L＝f(α，Ma)，其中f为函数名，Ma为飞行器当前马赫数；

除以上因素外，目标重要程度也应对概率推断产生影响，作为非状态相关项，通过设置非均一分布的初始概率来体现目标重要程度，定义条件概率初值：

式中I_θ(·)表示为打击目标的重要程度。

再进一步地，所述步骤二中(1)以“区域”的概念将飞行环境离散，将对位置参数的概率估计转化为对区域的概率估计的具体过程为：

对于实际长度为L×W大小的战场ψ，将其重新划分为M×N单位长度的区域，区域最小单元为边长为e的正方形；

飞行器每一时刻的实际位置即可用当前所在区域表示，采用<ν_i,ν_j>来表示区域i到区域j之间的转移过程。

进一步地，所述步骤二中(2)采用简化动力学模型快速预测所有可行区间的具体过程为：

在二维平面内，飞行器简化动力学方程形式：

式中θ、φ、ψ为状态变量，通过积分获得未来取值；

其他参数中，速度

基于当前观测数据，以最小二乘拟合的形式得到预测时间处的值；路径角γ₀、地心距r₀取最后观测数据相应值，在积分中保持为常数；

为了生成所有可行区间，在预测中保持攻角为极值α_max，倾侧角分别取[-σ_max,-0.5σ_max,0,0.5σ_max,σ_max]，令预测时间步长为dT_p，通过线性插值上述5条轨迹在预测时间n·dT_p处的位置参数，获得每一预测时间点的可行区间，得到对应的区域集合，再通过禁飞区约束，去除部分实际不可达区域后，记第n步预测可行区间集合为Uⁿ，其中元素满足

对于一步状态转移过程

认为集合Uⁿ中所有元素都能映射到Uⁿ⁺¹中，转移概率由意图代价函数确定；

如果某一目标点出现在了插值区域上，则在积分过程中获得相应的仿真时长；如果目标点不在可行区域内，则加以判断。

再进一步地，所述步骤二中如果目标点不在可行区域内，则加以判断的方法为：按照数学的原理确定点是否位于多边形内，用“光线投射算法”进行演算，从点出发生成一条射线，当射线与多边形边相交次数为奇数时，点位于多边形内，偶数时则位于多边形外。

进一步地，所述步骤三的具体过程为：

基于已有的观测信息x_1:i确定目标θ的概率，进行意图推测；

在贝叶斯框架下，Pr(θ|X_1:i＝x_1:i)表述为目标关于观测信息的后验概率，递推关系：

假设系统的状态转移服从一阶马尔可夫模型：当前时刻状态x_i只与上一时刻状态x_i-1有关，进而下列关系成立：

Pr(X_i|X_1:i-1＝x_1:i-1,θ)＝Pr(X_i|X_i-1＝x_i-1,θ) (0.10)

上式(0.9)中，分母项Pr(X_i|X_1:i-1＝x_1:i-1)完全由量测噪声模型确定，为一常数，计算中可舍去，因此后验更新公式：

Pr(θ|X_1:i＝x_1:i)∝Pr(X_i|X_i-1＝x_i-1,θ)·Pr(θ|X_1:i-1＝x_1:i-1) (0.11)

式(0.11)右侧第一项Pr(X_i|X_i-1＝x_i-1,θ)为利用步骤一定义结合观测数据计算得到似然概率，在初始概率分布Pr(θ|X₁＝x₁)下，每一观测时刻都可通过迭代计算出当前所有目标后验概率分布。

所述步骤三意图推测的算法流程如下：

(1)更新当前时刻t_i获得的量测数据x_i，现有观测数据集为X_1:i＝x_1:i；

(2)利用t_i-1时刻量测数据集X_1:i-1，构建速度拟合曲线，确定一步积分初值，生成一步观测时间内可行区域Uⁱ，并判断可达目标；

(3)分别针对所有可达目标，利用前一时刻观测数据x_i-1与可行域Uⁱ计算归一化系数K；利用当前观测x_i与归一系数K，计算在不同目标θ下实际运动对应的似然概率Pr(X_i|X_1:i-1＝x_1:i-1,θ)；

(4)结合前一步对不同目标θ的后验估计Pr(θ|X_1:i-1＝x_1:i-1)，利用式(0.11)完成概率更新，选择后验概率最大者为可能打击目标。

进一步地，所述步骤四的具体过程为：

采用步骤二过程离散了飞行环境后，使轨迹预测变成了区域预测，依据似然概率一步外推，得到下一时刻状态的概率分布：

上式中，第一行由全概率公式得到，第二、三行由联合概率与条件概率转换关系得到，从第三行到第四行应用一阶马尔可夫假设得到；

右侧第一项Pr(X_i+1|X_i＝x_i,θ＝θ_η)根据意图代价函数式(0.1)在当前时刻一步递推计算；计算当前步到未来时刻的似然概率；

右侧第二项Pr(θ＝θ_η|X_1:i＝x_1:i)为目标θ_η对应的后验概率，在贝叶斯框架下，基于意图估计得到轨迹预测；

多步预测公式在此基础上展开：

求和项中χ-表示对于某个特定位置X_i+j+1，所有一步转移可行的x_i+j所在区域；该式推导方式同上；右侧第二项Pr(X_i+j＝x_i+j|X_1:i＝x_1:i)为左侧项Pr(X_i+j+1|X_1:i＝x_1:i)的上一步值，因此该公式具有迭代结构，当得到一步预测概率分布式(0.12)时，得到所有时刻概率分布；

右侧第一项Pr(X_i+j+1|X_i+j＝x_i+j)为状态转移概率，理论上由系统状态方程完全确定，在计算时可由全概率公式获得：

对于式中第一项Pr(X_i+j+1|X_i+j＝x_i+j,θ＝θ_η)表示未来时刻的似然概率，在运算中只分析出了各步可行域，利用意图代价函数得出；

对于后验概率项Pr(θ＝θ_η|X_i+j＝x_i+j)，进行如下假设：在预测未来状态转移时，目标概率保持不变，即：

Pr(θ＝θ_η|X_i+j＝x_i+j)≈Pr(θ＝θ_η|X_1:i＝x_1:i) (0.15)

基于客观的观测数据，推测出了当前的目标后验概率，将此应用于当前时刻对未来状态预测；

依据多步预测公式逻辑，提出一种随机采样打靶算法进行后续预测。

再进一步地，所述一种随机采样打靶算法，算法流程如下：

(1)更新当前时刻t_i获得的量测数据x_i，现有观测数据集为X_1:i＝x_1:i；

(2)采用步骤三设计的“意图推测算法”，完成对不同目标θ的后验估计更新Pr(θ|X_1:i＝x_1:i)；

(3)利用当前t_i时刻量测数据集X_1:i，构建速度拟合曲线，并确定多步积分初值；

以dT_p为一步预测间隔，生成从观测位置至各目标点θ的一步转移可行区域

并预测各目标到达时间

(4)初始化M×N×T维的计数矩阵，其中M×N维对应轨迹预测时进行离散后所有飞行区域，T对应抵达所有目标时最大间隔步数；

(5)定义蒙特卡洛打靶次数为N_all，针对不同目标θ基于后验估计分配打靶次数N_j＝Pr(θ|X_1:i＝x_1:i)×N_all；

对于这N_j次打靶，轨迹终点保持为目标θ，轨迹起点为观测点所对应区域；

一次打靶中，每步状态转移＜ν_m,ν_n＞依据似然概率Pr(X_n|X_m＝x_m,θ＝θ_η)随机进行；

完成一步转移后，为当前时间步对应的M×N维计数阵中ν_n处元素加1；

(6)在所有打靶完成后，未来区域概率分布即为计数阵中每层元素值与总打靶次数N_all比值。

发明原理：多指标意图函数建立、高超声速飞行器可行域分析、贝叶斯估计下意图推测算法和基于蒙特卡洛打靶的高超声速飞行器轨迹预报方法。其中多指标意图函数建立与可行域分析两部分，是意图推测与轨迹预测内容的基础，从如何表述运动概率与如何确定运动区间两方面完整描述了飞行器在意图驱使下的运动过程。

多指标意图函数以似然概率的形式定义了目标运动意图，综合考虑了航向误差、弹目距离、目标重要度、避开危险区域四种常见作战因素。航向误差考虑的是，飞行器飞向目标点过程中，随弹目距离接近，弹目连线与运动方向夹角会逐渐减小；弹目距离项则基于以下事实，为降低速度损失，飞行器往往会采取最短路径接近目标；避开危险区域项，考虑飞行器机动能力有限，在朝目标飞行前会先避开禁飞区飞行；目标重要度，体现在目标点初始概率分布不同上。基于此，在同时存在多个目标点时，依据上述准则可判断出下一步运动的似然概率分布。

高超声速飞行器可行域分析解决如何量化未来可行区域的问题。在贝叶斯框架下，对未来状态递推估计时，需要确定所有可行的位置。从数学意义而言，对于非平凡问题，当且仅当时间离散并且每一时间步只存在有限次可能的状态转移时，理论方程才是可计算的。因此，该部分首先提出了以区域代替位置坐标的离散方法，将一般轨迹预测中未知状态预测转化为对区域概率的估计。其次，根据跟踪数据以极限倾侧角、攻角，理论计算了未来可行空域包络，并以此初步确定可行目标点。

贝叶斯估计下意图推测算法完成的是基于已有的观测信息，确定目标概率分布的过程。具体而言，利用上一时刻状态计算下一步似然概率分布，结合当前实际跟踪状态更新后验估计。在这一部分能获得随跟踪时间变化的目标意图曲线，将作为轨迹预测的基准。

基于蒙特卡洛打靶的高超声速飞行器轨迹预报以概率的形式给出了未来可能空域的分布。基于贝叶斯多步递推预测方法，计算复杂度较高，因而此处依据意图的后验概率设计了一种新的打靶预测方法。在每一次打靶中，未来一步状态转移概率依似然函数确定，转移过程随机生成，仿真点分布情况以轨迹簇形式体现。点越集中，说明该轨迹对应概率较高，推断为目标可能的飞行路径。

有益效果：

提供一种贝叶斯框架下的高超声速飞行器轨迹预报方法，可完成对目标飞行意图的评估和长时间的轨迹预测。本发明旨在从战略层面上，为空天防御体系设计者提供一种全新模式、符合工程应用的轨迹预报方法，一方面揣测目标飞行意图，降低未来机动致使的轨迹不确定性；另一方面设计实时未来轨迹多步递推算法，充分挖掘高超声速飞行器运动的潜在规律，实现目标轨迹的长时间预报。

附图说明

图1为轨迹预测示意图；

图2为禁飞区示意图；

图3为飞行可行域图示；

图4为目标点判断图示；

图5为飞行仿真环境；

图6(a)为100秒时刻意图推测结果的意图推测曲线，图6(b)为100秒时刻意图推测结果的识别目标曲线；

图7(a)为100秒时刻预测空间概率分布，图7(b)为100秒时刻预测最大概率分布，图7(c)为100秒时刻预测10步概率分布，图7(d)为100秒时刻不同目标路径最高概率变化；

图8(a)为300秒时刻预测空间概率分布，图8(b)为300秒时刻预测最大概率分布，图8(c)为300秒时刻预测10步概率分布，图8(d)为300秒时刻不同目标路径最高概率变化；

图9(a)为400秒时刻预测空间概率分布，图9(b)为400秒时刻预测最大概率分布，图9(c)为400秒时刻预测10步概率分布，图9(d)为400秒时刻不同目标路径最高概率变化。

具体实施方式

具体实施方式一：一种基于贝叶斯估计的高超速飞行器长期轨迹预报方法，其特征在于：它包括以下步骤：

步骤一：建立意图代价函数；

在水平面内分析轨迹预测问题，在飞行区域中，存在若干攻击目标，假设高超飞行器作战意图主要体现在航向误差、弹目距离以及战场环境上，并结合航向误差、弹目距离变化以及战场环境定义飞行器意图代价函数；

步骤二：进行可行域分析；

(1)以“区域”的概念将飞行环境离散，将对位置参数的概率估计转化为对区域的概率估计；

(2)采用简化动力学模型快速预测所有可行区间；

步骤三：贝叶斯估计框架下，进行目标意图推测；

贝叶斯估计框架下，基于已有的观测信息确定攻击目标的概率，进行意图推测；

步骤四：基于贝叶斯估计和蒙特卡洛打靶法进行轨迹预测；；

具体实施方式二：所述步骤一具体为：

在飞行区域ψ中，存在若干攻击目标θ，定义飞行器意图代价函数为：

归一化系数K计算方式如下：

式中θ_η∈Θ表示某个已知的攻击目标，φ_k∈Φ表示战场区域内某一禁飞区，χ⁺表示可行集”；

公式(0.1)中右侧第一项a·(δ(x_i,X_i+1,θ_η)-δ(x_i,θ_η))表示弹目距离变化，用δ(·)表示括号内元素最短距离，有以下关系：

δ(a,b,c)＝δ(a,b)+δ(b,c) (0.3)

式中，a表示弹目距离的权重系数、b表示航向角的权重系数和c表示禁飞区的权重系数；

该项值恒为非负数，当飞行器沿最短路径接近目标时，为极限值0；权重系数a表征了飞行器以最小路径前往目标的潜在意图，当权重系数a→+∞，朝着最短路径移动会获得最大条件概率；当权重系数a→0，朝着各个方向移动条件概率相当；

公式(0.1)中右侧第二项b·|ν(x_i,X_i+1)-ν(x_i,θ_η)|/δ(x_i,θ_η)表示航向角偏差，ν(·)表示括号内元素连线与正北方向夹角；用Δν(·)表示该项分子，其值为0时，代表飞行器速度方向指向了该目标，这一移动会获得最大的条件概率；同时，考虑到当飞行器离目标较远时，运动方向更自由，在该项分母引入了弹目距离δ(x_i,θ_η)，使其获得自适应的权重系数；越接近目标，与意图函数中其他项相比，航向角偏差项权重更大，

公式(0.1)中右侧第三项c·∑I_φ(x_i,X_i+1,φ_k)表示禁飞区约束影响，定义为：

为第k个禁飞区半径,

、

为当前速度方向与禁飞区切线夹角，当速度矢量在切线角内时，夹角值定义为正，若飞行器朝向禁飞区飞行，

越大，代价也越大，会获得较小条件概率，

为安全裕度，当飞行器离禁飞区较远时，运动自由，安全裕度较大，任意朝向均不受禁飞区约束影响；靠近禁飞区时，安全裕度减小，只有避开禁飞区方向才能获得较大概率，下面给出

值计算方式；

若瞬时速度方向指向禁飞区中心，则

；当飞行器恰好可以极限转弯半径r_p从禁飞区边缘飞过，定义此时速度方向与切线夹角为

根据三角函数关系可知：

高超飞行器侧向机动能力主要由升力分量提供，因此转弯半径可由下式估算：

该式在应用时，状态参数取跟踪时刻末值，攻角α、倾侧角σ取为极限值，以估计极限转弯半径r_p，m为飞行器质量，v为飞行器速度，L为升力，C_L为飞行器升力系数，ρ为当前大气密度，S为飞行器特征面积，σ为飞行器倾侧角，

攻角α与飞行器升力系数C_L存在关系C_L＝f(α，Ma)，其中f为函数名，Ma为飞行器当前马赫数；

除以上因素外，目标重要程度也应对概率推断产生影响，作为非状态相关项，通过设置非均一分布的初始概率来体现目标重要程度，定义条件概率初值：

式中I_θ(·)表示为打击目标的重要程度。

其他具体实施方式与具体实施方式一相同。

具体实施方式三：所述步骤二中(1)以“区域”的概念将飞行环境离散，将对位置参数的概率估计转化为对区域的概率估计的具体过程为：

对于实际长度为L×W大小的战场ψ，将其重新划分为M×N单位长度的区域，区域最小单元为边长为e的正方形；

飞行器每一时刻的实际位置即可用当前所在区域表示，采用<ν_i,ν_j>来表示区域i到区域j之间的转移过程。

其他具体实施方式与具体实施方式一相同。

具体实施方式四：所述步骤二中(2)采用简化动力学模型快速预测所有可行区间的具体过程为：

在二维平面内，飞行器简化动力学方程形式：

式中θ、φ、ψ为状态变量，通过积分获得未来取值；

其他参数中，速度v基于当前观测数据，以最小二乘拟合的形式得到预测时间处的值；路径角γ₀、地心距r₀取最后观测数据相应值，在积分中保持为常数；

为了生成所有可行区间，在预测中保持攻角为极值α_max，倾侧角分别取[-σ_max,-0.5σ_max,0,0.5σ_max,σ_max]，令预测时间步长为dT_p，通过线性插值上述5条轨迹在预测时间n·dT_p处的位置参数，获得每一预测时间点的可行区间，得到对应的区域集合，再通过禁飞区约束，去除部分实际不可达区域后，记第n步预测可行区间集合为Uⁿ，其中元素满足

对于一步状态转移过程

认为集合Uⁿ中所有元素都能映射到Uⁿ⁺¹中，转移概率由意图代价函数确定；

如果某一目标点出现在了插值区域上，则在积分过程中获得相应的仿真时长；如果目标点不在可行区域内，则加以判断。

其他具体实施方式与具体实施方式一相同。

具体实施方式五：所述步骤二中如果目标点不在可行区域内，则加以判断的方法为：按照数学的原理确定点是否位于多边形内，用“光线投射算法”进行演算，从点出发生成一条射线，当射线与多边形边相交次数为奇数时，点位于多边形内，偶数时则位于多边形外。

其他具体实施方式与具体实施方式一相同。

具体实施方式六：所述步骤三的具体过程为：

基于已有的观测信息x_1:i确定目标θ的概率，进行意图推测；

在贝叶斯框架下，Pr(θ|X_1:i＝x_1:i)表述为目标关于观测信息的后验概率，递推关系：

假设系统的状态转移服从一阶马尔可夫模型：当前时刻状态x_i只与上一时刻状态x_i-1有关，进而下列关系成立：

Pr(X_i|X_1:i-1＝x_1:i-1,θ)＝Pr(X_i|X_i-1＝x_i-1,θ) (0.10)

上式(0.9)中，分母项Pr(X_i|X_1:i-1＝x_1:i-1)完全由量测噪声模型确定，为一常数，计算中可舍去，因此后验更新公式：

Pr(θ|X_1:i＝x_1:i)∝Pr(X_i|X_i-1＝x_i-1,θ)·Pr(θ|X_1:i-1＝x_1:i-1) (0.11)

式(0.11)右侧第一项Pr(X_i|X_i-1＝x_i-1,θ)为利用步骤一定义结合观测数据计算得到似然概率，在初始概率分布Pr(θ|X₁＝x₁)下，每一观测时刻都可通过迭代计算出当前所有目标后验概率分布。

其他具体实施方式与具体实施方式一相同。

具体实施方式七：所述步骤三意图推测的算法流程如下：

(1)更新当前时刻t_i获得的量测数据x_i，现有观测数据集为X_1:i＝x_1:i；

(2)利用t_i-1时刻量测数据集X_1:i-1，构建速度拟合曲线，确定一步积分初值，生成一步观测时间内可行区域Uⁱ，并判断可达目标；

(3)分别针对所有可达目标，利用前一时刻观测数据x_i-1与可行域Uⁱ计算归一化系数K；利用当前观测x_i与归一系数K，计算在不同目标θ下实际运动对应的似然概率Pr(X_i|X_1:i-1＝x_1:i-1,θ)；

(4)结合前一步对不同目标θ的后验估计Pr(θ|X_1:i-1＝x_1:i-1)，利用式(0.11)完成概率更新，选择后验概率最大者为可能打击目标。

其他具体实施方式与具体实施方式一相同。

具体实施方式八：所述步骤四的具体过程为：

采用步骤二过程离散了飞行环境后，使轨迹预测变成了区域预测，依据似然概率一步外推，得到下一时刻状态的概率分布：

上式中，第一行由全概率公式得到，第二、三行由联合概率与条件概率转换关系得到，从第三行到第四行应用一阶马尔可夫假设得到；

右侧第一项Pr(X_i+1|X_i＝x_i,θ＝θ_η)根据意图代价函数式(0.1)在当前时刻一步递推计算；计算当前步到未来时刻的似然概率；

右侧第二项Pr(θ＝θ_η|X_1:i＝x_1:i)为目标θ_η对应的后验概率，在贝叶斯框架下，基于意图估计得到轨迹预测；

多步预测公式在此基础上展开：

求和项中χ^-表示对于某个特定位置X_i+j+1，所有一步转移可行的x_i+j所在区域；该式推导方式同上；右侧第二项Pr(X_i+j＝x_i+j|X_1:i＝x_1:i)为左侧项Pr(X_i+j+1|X_1:i＝x_1:i)的上一步值，因此该公式具有迭代结构，当得到一步预测概率分布式(0.12)时，得到所有时刻概率分布；

右侧第一项Pr(X_i+j+1|X_i+j＝x_i+j)为状态转移概率，理论上由系统状态方程完全确定，在计算时可由全概率公式获得：

对于式中第一项Pr(X_i+j+1|X_i+j＝x_i+j,θ＝θ_η)表示未来时刻的似然概率，在运算中只分析出了各步可行域，利用意图代价函数得出；

对于后验概率项Pr(θ＝θ_η|X_i+j＝x_i+j)，进行如下假设：在预测未来状态转移时，目标概率保持不变，即：

Pr(θ＝θ_η|X_i+j＝x_i+j)≈Pr(θ＝θ_η|X_1:i＝x_1:i) (0.15)

基于客观的观测数据，推测出了当前的目标后验概率，将此应用于当前时刻对未来状态预测；

依据多步预测公式逻辑，提出一种随机采样打靶算法进行后续预测。

其他具体实施方式与具体实施方式一相同。

具体实施方式九：所述一种随机采样打靶算法，算法流程如下：

(1)更新当前时刻t_i获得的量测数据x_i，现有观测数据集为X_1:i＝x_1:i；

(2)采用步骤三设计的“意图推测算法”，完成对不同目标θ的后验估计更新Pr(θ|X_1:i＝x_1:i)；

(3)利用当前t_i时刻量测数据集X_1:i，构建速度拟合曲线，并确定多步积分初值；

以dT_p为一步预测间隔，生成从观测位置至各目标点θ的一步转移可行区域

并预测各目标到达时间

(4)初始化M×N×T维的计数矩阵，其中M×N维对应轨迹预测时进行离散后所有飞行区域，T对应抵达所有目标时最大间隔步数；

(5)定义蒙特卡洛打靶次数为N_all，针对不同目标θ基于后验估计分配打靶次数N_j＝Pr(θ|X_1:i＝x_1:i)×N_all；

对于这N_j次打靶，轨迹终点保持为目标θ，轨迹起点为观测点所对应区域；

一次打靶中，每步状态转移<ν_m,ν_n>依据似然概率Pr(X_n|X_m＝x_m,θ＝θ_η)随机进行；

完成一步转移后，为当前时间步对应的M×N维计数阵中ν_n处元素加1；

(6)在所有打靶完成后，未来区域概率分布即为计数阵中每层元素值与总打靶次数N_all比值。

其他具体实施方式与具体实施方式八相同。

具体实施例：

1、意图代价函数

在飞行器接近目标途中，当前的飞行状态一定程度上隐含着飞行意图。在民航监测领域，当既定航线出现偏移时，地面站从航向偏移、位置偏移、TTG偏移、速度偏移等信息预测飞机意图。建立基于状态信息的意图函数，使得预测飞行轨迹成了可能。本发明中，轨迹预测问题只在水平面内分析，原因在于高超飞行器在飞行中常采用跳跃弹道模式，纵向难以推测出有效信息。

如附图1所示，在飞行区域ψ中，存在若干攻击目标θ，假设高超飞行器作战意图主要体现在航向误差、弹目距离以及战场环境上。综合以上因素，定义飞行器意图代价函数为：

归一化系数K计算方式如下：

式中θ_η∈Θ表示某个已知的攻击目标，φ_k∈Φ表示战场区域内某一禁飞区。右侧第一项表示弹目距离变化，用δ(·)表示括号内元素最短距离，有以下关系：

δ(a,b,c)＝δ(a,b)+δ(b,c) (0.3)

显然，该项值恒为非负数，当飞行器沿最短路径接近目标时，为极限值0。权重系数a表征了飞行器以最小路径前往目标的潜在意图，当权重系数α→+∞，朝着最短路径移动会获得最大条件概率；当权重系数α→0，朝着各个方向移动条件概率相当。

第二项表示航向角偏差，ν(·)表示括号内元素连线与正北方向夹角。用Δν(·)表示该项分子，其值为0时，代表飞行器速度方向指向了该目标，这一移动会获得最大的条件概率。同时，考虑到当飞行器离目标较远时，运动方向更自由，因此在该项分母引入了弹目距离，使其获得自适应的权重系数。越接近目标，与意图函数中其他项相比，航向角偏差项权重更大。

第三项表示禁飞区约束影响，定义为：

如附图2所示，

、

为当前速度方向与禁飞区切线夹角，当速度矢量在切线角内时，夹角值定义为正。若飞行器朝向禁飞区飞行，

越大，代价也越大，会获得较小条件概率。

为安全裕度，当飞行器离禁飞区较远时，运动自由，安全裕度较大，任意朝向均不受禁飞区约束影响；靠近禁飞区时，安全裕度减小，只有避开禁飞区方向才能获得较大概率。下面给出

值计算方式。

附图2中，若瞬时速度方向指向禁飞区中心，则

。当飞行器恰好可以极限转弯半径r_p从禁飞区边缘飞过，定义此时速度方向与切线夹角为

根据三角函数关系可知：

高超飞行器侧向机动能力主要由升力分量提供，因此转弯半径可由下式估算：

该式在应用时，状态参数取跟踪时刻末值，攻角、倾侧角取为极限值，以估计最小转弯半径。

除以上因素外，目标重要程度也应对概率推断产生影响。作为非状态相关项，不在意图函数中引入重要度，而通过设置非均一分布的初始概率来体现这一因素，定义条件概率初值：

式中I_θ(·)表示为打击目标的重要程度。

2、可行域分析

意图函数式(0.1)在应用过程会存在以下困难：1)运动过程是连续的，如何量化估计状态参数；2)在计算K值时，如何确定一步时间内的可行集χ⁺。从数学角度而言，“对于非平凡问题，当且仅当时间离散并且每一时间步只存在有限次可能的状态转移时，理论方程才是可计算的。”因而，本发明首先以“区域”的概念将飞行环境ψ离散，将对位置参数的概率估计转化为对区域的概率估计。

对于实际长度为L×W大小的战场ψ，将其重新划分为M×N单位长度的区域，区域最小单元为边长为e的正方形。飞行器每一时刻的实际位置即可用当前所在区域表示，采用<v_i,ν_j>来表示区域i到区域j之间的转移过程。这种处理后这一过程将是有限、可分析的。

对于问题2，一般性的处理方式是在每一时间点利用状态估计积分预测下一时间点新的状态。但在可行区域数量较大情况下，这种计算过于冗长。本发明采用简化动力学模型快速预测所有可行区间。二维平面内，飞行器简化动力学方程形式如下：

式(0.8)中θ、φ、ψ为状态变量，通过积分获得未来取值。其他参数中，速度

基于当前观测数据，以最小二乘拟合的形式得到预测时间处的值；路径角γ₀、地心距r₀取最后观测数据相应值，在积分中保持为常数。为了生成所有可行区间，在预测中保持攻角为极值α_max，倾侧角分别取[-σ_max,-0.5σ_max,0,0.5σ_max,σ_max]，令预测时间步长为dT_p，通过线性插值上述5条轨迹在预测时间n·dT_p处的位置参数，就获得了每一预测时间点的可行区间，得到对应的区域集合，如附图3所示。再通过禁飞区约束，去除部分实际不可达区域后，记第n步预测可行区间集合为Uⁿ，其中元素满足

对于一步状态转移过程

认为集合Uⁿ中所有元素都能映射到Uⁿ⁺¹中，转移概率由意图代价函数确定。应当声明以此方式描述转移过程是不够严谨的，并未考虑相邻步间元素在转移时存在的动力学约束。但此关系在一次积分方法下难以确定，因而使用概率来限制元素间转移，使得不满足动力学约束的转移获得较低概率。如意图函数中第一项与第二项，当转移发生在距离较远、方向差异较大的两个区域间时，距离变化函数与方向角误差函数确保了这一转移模式不受青睐。

本节最后讨论预测到达时间与可能目标的问题。前文设计了通过5条标准轨迹来生成所有可行区间的方法，如果某一目标点出现在了插值区域上，相应的仿真时长就很容易在积分过程中获得。但相反，如果目标点不在可行区域内，需加以判断。从数学上，此问题实质是确定点是否位于多边形内。“光线投射算法”能有效解决这一问题，其描述如下：从点出发生成一条射线，当射线与多边形边相交次数为奇数时，点位于多边形内，偶数时则位于多边形外。如附图4所示。

3、贝叶斯估计下意图推测算法

基于已有的观测信息x_1:i确定攻击目标θ的概率，这一过程称之为意图推测。在贝叶斯框架下，Pr(θ|X_1:i＝x_1:i)表述为目标关于观测信息的后验概率，有如下递推关系：

本发明中，假设系统的状态转移服从一阶马尔可夫模型：当前时刻状态x_i只与上一时刻状态x_i-1有关。因此，下式成立：

Pr(X_i|X_1:i-1＝x_1:i-1,θ)＝Pr(X_i|X_i-1＝x_i-1,θ) (0.10)

式(0.9)中，分母项Pr(X_i|X_1:i-1＝x_1:i-1)完全由量测噪声模型确定，为一常数，计算中可舍去。因此后验更新公式：

Pr(θ|X_1:i＝x_1:i)∝Pr(X_i|X_i-1＝x_i-1,θ)·Pr(θ|X_1:i-1＝x_1:i-1) (0.11)

式(0.11)中右侧第一项为似然概率，利用第一节定义结合观测数据能够计算得到。在初始概率分布Pr(θ|X₁＝x₁)下，每一观测时刻都可通过迭代计算出当前所有目标后验概率分布。意图推测算法流程如下：

(1)更新当前时刻t_i获得的量测数据x_i，现有观测数据集为X_1:i＝x_1:i；

(2)利用t_i-1时刻量测数据集X_1:i-1，构建速度拟合曲线，确定一步积分初值，生成一步观测时间内可行区域Uⁱ，并判断可达目标；

(3)分别针对所有可达目标，利用前一时刻观测数据x_i-1与可行域Uⁱ计算归一化系数K；利用当前观测x_i与归一系数K，计算在不同目标θ下实际运动对应的似然概率Pr(X_i|X_1:i-1＝x_1:i-1,θ)；

(4)结合前一步对不同目标θ的后验估计Pr(θ|X_1:i-1＝x_1:i-1)，利用式(0.11)完成概率更新，选择后验概率最大者为可能打击目标。

仿真环境设计如下，如附图5所示，在经纬度范围分别为[100,135]与[30,35]的区域内，存在3个目标点与3个禁飞区。高超声速飞行器模型参考CAV-H的基本参数，图中黑色虚线代表飞行器真实的飞行轨迹，通过轨迹规划方法生成，其实际意图是攻击目标2，飞行时间约为675.3秒。由于本发明不涉及具体的数据跟踪算法，因此观测数据由实际值混合白噪声得到，误差极值设计为1km。观测段数据采样周期为1秒，每隔10秒进行一次意图更新。各目标重要度初值相同，飞行器前100秒航行为初始观测段，不进行意图推测，目标概率保持不变。自100秒后，每隔10秒意图迭代更新一次，直到飞行器抵达目标点。相关参数如下：

表1仿真参数

目标概率随时间变化如附图6(a)所示，识别出的可能打击目标如附图6(b)所示。从意图推测曲线可以看出目标1的概率全程恒定为0，这是因为通过算法第2步判断出了该目标点不在可达区内。同理，500秒以后的跟踪末期，目标3因为飞行器已经从旁边经过，不在可达区内，推测概率恒为0。当两目标推测概率差值大于一定阈值时，认为能够显著识别进攻目标，否则二者均为可能目标。本发明中，阈值取ΔPr＝20％。因此在附图6(a)中，100秒至250秒段，认为目标2与目标3均是可能进攻的目标；250秒以后，目标3推测概率显著大于目标2，认为其是唯一可能进攻目标。此算例中，飞行器真实意图为攻击目标2，以正确推测时间与全程推测时间的比值来衡量算法的准确性，可以得到该算法效率为：

因此，利用该算法推测目标进攻意图具有一定实用价值。

4、基于蒙特卡洛打靶的轨迹预测算法

采用第2节方法离散了飞行环境后，轨迹预测实际上变成了区域预测。依据似然概率一步外推，可以得到下一时刻状态的概率分布：

上式中第一行由全概率公式得到，第二、三行由联合概率与条件概率转换关系得到，从第三行到第四行应用一阶马尔可夫假设可得。右侧第一项根据意图代价函数式(0.1)在当前时刻可一步递推计算。与前一节在意图估计中计算所不同之处在于，前者在当前时刻计算的是依据上一步到当前步的似然概率，后者计算的则是当前步到未来时刻的似然概率。第二项为目标θ_η对应的后验概率，上节算法在每一观测时刻均对该值进行了更新。因此，在贝叶斯框架下，轨迹预测基于意图估计得到。多步预测公式在此基础上展开：

求和项中χ^-表示对于某个特定位置X_i+j+1，所有一步转移可行的x_i+j所在区域。该式推导方式同上。右侧第二项为左侧项在上一步值，因此该公式具有迭代结构，当得到一步预测概率分布式(0.12)时，即可得到所有时刻概率分布。右侧第一项为状态转移概率，理论上由系统状态方程完全确定。在计算时可由全概率公式获得：

对于式中第一项Pr(X_i+j+1|X_i+j＝x_i+j,θ＝θ_η)表示未来时刻的似然概率，在运算中只要分析出了各步可行域，利用意图代价函数直接可得。但对于后验概率项Pr(θ＝θ_η|X_i+j＝x_i+j)，因为无法获得未来时刻的观测数据，故意图推测无法迭代更新。因此，为了计算方便，此时应进行如下假设：在预测未来状态转移时，目标概率保持不变。即：

Pr(θ＝θ_η|X_i+j＝x_i+j)≈Pr(θ＝θ_η|X_1:i＝x_1:i) (0.15)

应当声明，这一假设是十分合理且必要的。基于客观的观测数据，推测出了当前的目标后验概率，将此应用于当前时刻对未来状态预测。对于无法预知的未来，这一方式显然优于基于主观揣测未来状态。

至此，采用概率形式通过当前观测数据，预测未来概率分布已经可操作。但在实际应用中，受限于计算效率，多步预测式(0.13)存在以下问题：随着预测步增加，计算一步区域转移概率时，需要考虑的可行区间χ^-是呈指数增大的。产生这种现象原因在于，前文定义了“一步状态转移过程

时，认为集合Uⁿ中所有元素都能映射到Uⁿ⁺¹中”。假设在第n步共生成了N个可行区间，那么第n+1步可行区间数量应与前一步相关，设为β·N个，为完成一步计算需要的计算次数将为β·N²。在多步预测后期可行区间数量庞大情况下，这种计算使得轨迹预测丧失了实时性的特点。为此，本发明依据多步预测公式逻辑，提出了一种随机采样打靶算法，能在兼顾计算效率与求解精度上逼近真实概率分布情况。算法流程如下：

(1)更新当前时刻t_i获得的量测数据x_i，现有观测数据集为X_1:i＝x_1:i；

(2)采用步骤三设计的“意图推测算法”，完成对不同目标θ的后验估计更新Pr(θ|X_1:i＝x_1:i)；

(3)利用当前t_i时刻量测数据集X_1:i，构建速度拟合曲线，并确定多步积分初值；

以dT_p为一步预测间隔，生成从观测位置至各目标点θ的一步转移可行区域

并预测各目标到达时间

(4)初始化M×N×T维的计数矩阵，其中M×N维对应轨迹预测时进行离散后所有飞行区域，T对应抵达所有目标时最大间隔步数；

(5)定义蒙特卡洛打靶次数为N_all，针对不同目标θ基于后验估计分配打靶次数N_j＝Pr(θ|X_1:i＝x_1:i)×N_all；

对于这N_j次打靶，轨迹终点保持为目标θ，轨迹起点为观测点所对应区域；

一次打靶中，每步状态转移<ν_m,ν_n>依据似然概率Pr(X_n|X_m＝x_m,θ＝θ_η)随机进行；

完成一步转移后，为当前时间步对应的M×N维计数阵中ν_n处元素加1；

(6)在所有打靶完成后，未来区域概率分布即为计数阵中每层元素值与总打靶次数N_all比值。

算法说明如下：1)总体而言，整个轨迹预测算法需要对目标的后验估计信息，因此在第二步中，联合了前一节提及的意图推测算法；2)预测算法中第三步再次进行了可行区域确定，这与意图推测算法中所不同之处在于，此处利用的是当前观测信息来多步预测至目标点区域，并且考虑到预测精度要求，在轨迹预测算法中会重新以不同边长e划分飞行区域；3)基于多步预测公式(0.13)，迭代过程主要依靠转移概率函数Pr(X_i+j+1|X_i+j＝x_i+j)完成，因此在算法第五步中，采取了打靶次数分配与似然概率计算两种方法来近似这一过程。4)一次打靶过程，完全依据概率生成一条从观测点至目标点的轨迹，因此在第六步中通过统计这些轨迹位置分布情况，来代表真实的概率分布。

飞行环境仍如附图5所示，自观测段某一时刻起，开始预测到目标飞行轨迹。具体形式以空间位置概率分布的方式呈现，概率高的区域代表飞行器最有可能以此途径抵达目标。预测频率为每50秒更新一次概率分布，预测时步长设置为20秒，预测阶段最小区域单元边长设置为0.2°。轨迹预测环节相关参数如下：

表2仿真参数

轨迹预测结果如附图7(a)-9(d)所示。以图7(a)-7(d)为例，表现的是在目标飞行第100秒时预报至目标点的结果。其中图7(a)预测的是未来位置的概率分布，为了更好区别高、低概率空间，此处在体现可行区域时，按高概率Pr＞0.5、较高概率0.2＜Pr＜0.5、一般概率0.01＜Pr＜0.2、低概率Pr＜0.01四种程度划分空间，并对应深蓝、蓝、浅蓝、浅色四种不同程度颜色。图中扇状分布区域最广，呈浅色，表明这些位置是预测可行点，但运动经过概率较低。对于蓝色、深蓝色区域，是未来较大概率会经过的空域，体现了飞行器自当前位置运动到目标点的一种趋势性，如图中存在两种抵达目标的方式。

将各预测步中概率最高点连线，预测步中间区域概率做插值处理，就能得到更为直观的预测轨迹，如图7(b)所示。与一般预测方法推导出的轨迹含义不同，通过概率方式得到的曲线，仅代表在未来每一时刻，该轨迹上的点拥有最大概率，未来可能的飞行轨迹并不是唯一的。

该曲线每步具体对应的概率值绘制在了7(d)中。可以看出预测轨迹上的概率值并非一成不变。在预测的前几步时，由于可行域较少，并且到两目标的轨迹前部分重合，因此预测曲线上概率值较大Pr＞0.5。随预测时间增加，轨迹差异性增大，俩种曲线上概率变化各不相同，但最终抵达目标时会回归到该种情况下的概率极限值Pr(θ|X_1:i＝x_1:i)。

因此，预测路径上最高概率值的变化趋势一般为由高至低再至高，这与可行区域数量变化是相一致的。在预测刚开始阶段，由于预测时间短，未来可行区域较小，概率分布较为集中，轨迹预测容易；随预测时间增长，可行区域数量逐渐增大，概率均分在了多个区域上，因此概率极值下降，这也代表了这一阶段想准确估计出轨迹难度较高；在预测末期，虽然可行区间数量更大，但受到目标点约束，概率分布再次集中，这阶段准确的轨迹预测也较为容易。

图7(c)体现的是具体某一预测步时，在各个可行区域上概率分布情况，其横坐标为可行区域从左至右的编号排序。可以看出这些区域中大部分是低概率区域，然后是一般概率区域，极少高概率区域。整体呈现出单峰的特点，这与当前时刻是预测初期，预测轨迹差异性较小相关。

附图8(a)-(d)对应在目标飞行第300秒时预报至目标点的结果。其中空间概率分布曲线整体趋势与100秒预测时相同，存在两条通向目标点的可行路径。但路径上最高概率存在明显差异：300秒预测时，针对目标3路径上的概率值明显整体低于目标2。这一现象是因为轨迹预测算法是基于意图推断得来的，通过前一节算法分析，在第250秒已经判断得出目标2的意图概率显著高于目标3，因此预测结果也体现了这一判断。

至于具体到某一步的概率分布情况，附图8(c)又与图7(c)不同。可以看到，300秒时刻，同样是第10步预测，呈现出的是双峰的特点。因为此时预测已接近中期，针对不同目标，轨迹的差异性已经展现。

附图9(a)-(d)对应在目标飞行第400秒时预报至目标点的结果。其空间概率分布，现仅存一条指向目标的轨迹，这是因为目标3已不包括在飞行器的可达区内。此时预测退化为单目标的轨迹推测。从最高概率变化曲线可以看出，仍是体现出预测中段分布概率低，两头分布概率高的特点。

综上所述，采用本发明设计的基于意图推测的轨迹预测算法，能很好的解决多目标点预测问题，给出了意图推测结果曲线、轨迹概率分布曲线，具有较高准确性与实时性，满足工程应用需求。

以上所述仅为本发明的较佳实施例而已，并不用于限制本发明，对于机器人领域多目标轨迹预测、车辆行为推测等本发明同样适用。凡在本发明原则和精神之内所作的任何修改、等同替换和改进等，均就包含在本发明的保护范围之内。

Claims (6)

1.一种基于贝叶斯估计的高超速飞行器长期轨迹预报方法，其特征在于：它包括以下步骤：

步骤一：建立意图代价函数；

在水平面内分析轨迹预测问题，在飞行区域中，存在若干攻击目标，假设高超飞行器作战意图主要体现在航向误差、弹目距离以及战场环境上，并结合航向误差、弹目距离变化以及战场环境定义飞行器意图代价函数；

步骤二：进行可行域分析；

(1)以“区域”的概念将飞行环境离散，将对位置参数的概率估计转化为对区域的概率估计；

(2)采用简化动力学模型快速预测所有可行区间；

步骤三：贝叶斯估计框架下，进行目标意图推测；

贝叶斯估计框架下，基于已有的观测信息确定攻击目标的概率，进行意图推测；

步骤四：基于贝叶斯估计和蒙特卡洛打靶法进行轨迹预测；

所述步骤一具体为：

在飞行区域ψ中，存在若干攻击目标θ，定义飞行器意图代价函数为：

归一化系数K计算方式如下：

式中θ_η∈Θ表示某个已知的攻击目标，φ_k∈Φ表示战场区域内某一禁飞区，χ⁺表示可行集”；

公式(0.1)中右侧第一项a·(δ(x_i,X_i+1,θ_η)-δ(x_i,θ_η))表示弹目距离变化，用δ(·)表示括号内元素最短距离，有以下关系：

δ(a,b,c)＝δ(a,b)+δ(b,c) (0.3)

式中，a表示弹目距离的权重系数、b表示航向角的权重系数和c表示禁飞区的权重系数；

该项值恒为非负数，当飞行器沿最短路径接近目标时，为极限值0；权重系数a表征了飞行器以最小路径前往目标的潜在意图，当权重系数a→+∞，朝着最短路径移动会获得最大条件概率；当权重系数a→0，朝着各个方向移动条件概率相当；

公式(0.1)中右侧第二项b·|ν(x_i,X_i+1)-ν(x_i,θ_η)|/δ(x_i,θ_η)表示航向角偏差，ν(·)表示括号内元素连线与正北方向夹角；用Δν(·)表示该项分子，其值为0时，代表飞行器速度方向指向了该目标，这一移动会获得最大的条件概率；同时，考虑到当飞行器离目标较远时，运动方向更自由，在该项分母引入了弹目距离δ(x_i,θ_η)，使其获得自适应的权重系数；越接近目标，与意图函数中其他项相比，航向角偏差项权重更大，

公式(0.1)中右侧第三项c·∑I_φ(x_i,X_i+1,φ_k)表示禁飞区约束影响，定义为：

为第k个禁飞区半径,

为当前速度方向与禁飞区切线夹角，当速度矢量在切线角内时，夹角值定义为正，若飞行器朝向禁飞区飞行，

越大，代价也越大，会获得较小条件概率，

为安全裕度，当飞行器离禁飞区较远时，运动自由，安全裕度较大，任意朝向均不受禁飞区约束影响；靠近禁飞区时，安全裕度减小，只有避开禁飞区方向才能获得较大概率，下面给出

值计算方式；

若瞬时速度方向指向禁飞区中心，则

当飞行器恰好可以极限转弯半径r_p从禁飞区边缘飞过，定义此时速度方向与切线夹角为

根据三角函数关系可知：

高超飞行器侧向机动能力主要由升力分量提供，因此转弯半径可由下式估算：

该式在应用时，状态参数取跟踪时刻末值，攻角α、倾侧角σ取为极限值，以估计极限转弯半径r_p，m为飞行器质量，v为飞行器速度，L为升力，C_L为飞行器升力系数，ρ为当前大气密度，S为飞行器特征面积，σ为飞行器倾侧角，

攻角α与飞行器升力系数C_L存在关系C_L＝f(α，Ma)，其中f为函数名，Ma为飞行器当前马赫数；

除以上因素外，目标重要程度也应对概率推断产生影响，作为非状态相关项，通过设置非均一分布的初始概率来体现目标重要程度，定义条件概率初值：

式中I_θ(·)表示为打击目标的重要程度；

所述步骤二中(1)以“区域”的概念将飞行环境离散，将对位置参数的概率估计转化为对区域的概率估计的具体过程为：

对于实际长度为L×W大小的战场ψ，将其重新划分为M×N单位长度的区域，区域最小单元为边长为e的正方形；

飞行器每一时刻的实际位置即可用当前所在区域表示，采用<ν_i,ν_j>来表示区域i到区域j之间的转移过程；

所述步骤四的具体过程为：

采用步骤二过程离散了飞行环境后，使轨迹预测变成了区域预测，依据似然概率一步外推，得到下一时刻状态的概率分布：

上式中，第一行由全概率公式得到，第二、三行由联合概率与条件概率转换关系得到，从第三行到第四行应用一阶马尔可夫假设得到；

右侧第一项Pr(X_i+1|X_i＝x_i,θ＝θ_η)根据意图代价函数式(0.1)在当前时刻一步递推计算；计算当前步到未来时刻的似然概率；

右侧第二项Pr(θ＝θ_η|X_1:i＝x_1:i)为目标θ_η对应的后验概率，在贝叶斯框架下，基于意图估计得到轨迹预测；

多步预测公式在此基础上展开：

求和项中χ^-表示对于某个特定位置X_i+j+1，所有一步转移可行的x_i+j所在区域；该式推导方式同上；右侧第二项Pr(X_i+j＝x_i+j|X_1:i＝x_1:i)为左侧项Pr(X_i+j+1|X_1:i＝x_1:i)的上一步值，因此该公式具有迭代结构，当得到一步预测概率分布式(0.12)时，得到所有时刻概率分布；

右侧第一项Pr(X_i+j+1|X_i+j＝x_i+j)为状态转移概率，理论上由系统状态方程完全确定，在计算时可由全概率公式获得：

对于式中第一项Pr(X_i+j+1|X_i+j＝x_i+j,θ＝θ_η)表示未来时刻的似然概率，在运算中只分析出了各步可行域，利用意图代价函数得出；

对于后验概率项Pr(θ＝θ_η|X_i+j＝x_i+j)，进行如下假设：在预测未来状态转移时，目标概率保持不变，即：

Pr(θ＝θ_η|X_i+j＝x_i+j)≈Pr(θ＝θ_η|X_1:i＝x_1:i) (0.15)

基于客观的观测数据，推测出了当前的目标后验概率，将此应用于当前时刻对未来状态预测。

2.根据权利要求1所述的一种基于贝叶斯估计的高超速飞行器长期轨迹预报方法，其特征在于：所述步骤二中(2)采用简化动力学模型快速预测所有可行区间的具体过程为：

在二维平面内，飞行器简化动力学方程形式：

式中θ、φ、ψ为状态变量，通过积分获得未来取值；

其他参数中，速度

基于当前观测数据，以最小二乘拟合的形式得到预测时间处的值；路径角γ₀、地心距r₀取最后观测数据相应值，在积分中保持为常数；

为了生成所有可行区间，在预测中保持攻角为极值α_max，倾侧角分别取[-σ_max,-0.5σ_max,0,0.5σ_max,σ_max]，令预测时间步长为dT_p，通过线性插值上述5条轨迹在预测时间n·dT_p处的位置参数，获得每一预测时间点的可行区间，得到对应的区域集合，再通过禁飞区约束，去除部分实际不可达区域后，记第n步预测可行区间集合为Uⁿ，其中元素满足

对于一步状态转移过程

认为集合Uⁿ中所有元素都能映射到Uⁿ⁺¹中，转移概率由意图代价函数确定；

如果某一目标点出现在了插值区域上，则在积分过程中获得相应的仿真时长；如果目标点不在可行区域内，则加以判断。

3.根据权利要求2所述的一种基于贝叶斯估计的高超速飞行器长期轨迹预报方法，其特征在于：所述步骤二中如果目标点不在可行区域内，则加以判断的方法为：按照数学的原理确定点是否位于多边形内，用“光线投射算法”进行演算，从点出发生成一条射线，当射线与多边形边相交次数为奇数时，点位于多边形内，偶数时则位于多边形外。

4.根据权利要求1所述的一种基于贝叶斯估计的高超速飞行器长期轨迹预报方法，其特征在于：所述步骤三的具体过程为：

基于已有的观测信息x_1:i确定目标θ的概率，进行意图推测；

在贝叶斯框架下，Pr(θ|X_1:i＝x_1:i)表述为目标关于观测信息的后验概率，递推关系：

假设系统的状态转移服从一阶马尔可夫模型：当前时刻状态x_i只与上一时刻状态x_i-1有关，进而下列关系成立：

Pr(X_i|X_1:i-1＝x_1:i-1,θ)＝Pr(X_i|X_i-1＝x_i-1,θ) (0.10)

上式(0.9)中，分母项Pr(X_i|X_1:i-1＝x_1:i-1)完全由量测噪声模型确定，为一常数，计算中可舍去，因此后验更新公式：

Pr(θ|X_1:i＝x_1:i)∝Pr(X_i|X_i-1＝x_i-1,θ)·Pr(θ|X_1:i-1＝x_1:i-1) (0.11)

式(0.11)右侧第一项Pr(X_i|X_i-1＝x_i-1,θ)为利用步骤一定义结合观测数据计算得到似然概率，在初始概率分布Pr(θ|X₁＝x₁)下，每一观测时刻都可通过迭代计算出当前所有目标后验概率分布。

5.根据权利要求1所述的一种基于贝叶斯估计的高超速飞行器长期轨迹预报方法，其特征在于：所述步骤三意图推测的算法流程如下：

(1)更新当前时刻t_i获得的量测数据x_i，现有观测数据集为X_1:i＝x_1:i；

(2)利用t_i-1时刻量测数据集X_1:i-1，构建速度拟合曲线，确定一步积分初值，生成一步观测时间内可行区域Uⁱ，并判断可达目标；

(3)分别针对所有可达目标，利用前一时刻观测数据x_i-1与可行域Uⁱ计算归一化系数K；利用当前观测x_i与归一系数K，计算在不同目标θ下实际运动对应的似然概率Pr(X_i|X_1:i-1＝x_1:i-1,θ)；

(4)结合前一步对不同目标θ的后验估计Pr(θ|X_1:i-1＝x_1:i-1)，利用式(0.11)完成概率更新，选择后验概率最大者为可能打击目标。

6.根据权利要求5所述的一种基于贝叶斯估计的高超速飞行器长期轨迹预报方法，其特征在于：当前时刻对未来状态预测具体为：依据多步预测公式逻辑，提出一种随机采样打靶算法进行后续预测；

所述一种随机采样打靶算法，算法流程如下：

(1)更新当前时刻t_i获得的量测数据x_i，现有观测数据集为X_1:i＝x_1:i；

(2)采用步骤三设计的“意图推测算法”，完成对不同目标θ的后验估计更新Pr(θ|X_1:i＝x_1:i)；

(3)利用当前t_i时刻量测数据集X_1:i，构建速度拟合曲线，并确定多步积分初值；

以dT_p为一步预测间隔，生成从观测位置至各目标点θ的一步转移可行区域

并预测各目标到达时间

(4)初始化M×N×T维的计数矩阵，其中M×N维对应轨迹预测时进行离散后所有飞行区域，T对应抵达所有目标时最大间隔步数；

(5)定义蒙特卡洛打靶次数为N_all，针对不同目标θ基于后验估计分配打靶次数N_j＝Pr(θ|X_1:i＝x_1:i)×N_all；

对于这N_j次打靶，轨迹终点保持为目标θ，轨迹起点为观测点所对应区域；

一次打靶中，每步状态转移<ν_m,ν_n>依据似然概率Pr(X_n|X_m＝x_m,θ＝θ_η)随机进行；

完成一步转移后，为当前时间步对应的M×N维计数阵中ν_n处元素加1；

(6)在所有打靶完成后，未来区域概率分布即为计数阵中每层元素值与总打靶次数N_all比值。

Images