CN107544067B - 一种基于高斯混合近似的高超声速再入飞行器跟踪方法 - Google Patents

Abstract

本发明公开了一种基于高斯混合近似的高超声速再入飞行器跟踪方法，在东北上惯性参考坐标系下建立高超声速机动再入飞行器的运动跟踪模型，并分别建立目标跟踪系统的雷达观测模型和目标先验模型集；根据先验模型集和雷达观测模型，并结合雷达回波数据,获得目标状态平滑估计值和协方差矩阵平滑估计值；本发明引入高斯混合近似思想，设计多模型下的高斯混合平滑滤波器，采用高斯分量自适应合并策略，同时合理选择平滑窗长，可避免粒子滤波存在的计算大和粒子退化问题，计算过程简单，易于实现，对于突然机动有着很好的估计精度，提高了高超声速再入飞行器的快速高精度跟踪能力。

Description

一种基于高斯混合近似的高超声速再入飞行器跟踪方法

【技术领域】

本发明属于雷达目标跟踪技术领域，尤其涉及一种基于高斯混合近似的高超声速再入飞行器跟踪方法。

【背景技术】

近年来，对于高超声速飞行器再入机动段的主动跟踪逐渐成为军事领域的研究热点。高超声速飞行器采用著名的“钱学森弹道”或“Sanger弹道”，其运动过程一般分为主动段、自由段、再入段。再入段是指目标重新返回进入稠密大气层的阶段，通常认为飞行器离地面高度80-100km即进入再入段。再入段飞行器相比飞机和普通导弹等飞行器有着特殊之处，它飞行快、速度高，再入过程受到的空气动力影响更为突出，加上对机动性的要求，再入过程中实施机动变轨飞行，机动方式种类繁多、复杂多变，机动加速度可高达5-15g。目标运动呈现出很强的非线性特性，气动模型及其复杂，微小的模型误差便会影响跟踪器的性能和准确性，这些对雷达精确跟踪和定位提出了严峻的挑战。

目前，对于再入机动目标跟踪问题主要集中在运动建模和非线性滤波两个方面。在运动建模方面，常见模型有协同转弯、匀加速，Singer、当前统计、Jerk等模型。上述模型将加速度建模为随机过程，未能利用飞行器的气动参数和空气动力学，不能准确描述目标真实运动特性，而基于阻力系数，升阻比等气动参数的动力学模型可精细地描述目标运动，提高目标跟踪精度。在滤波方面，常见有扩展卡尔曼滤波(EKF)，不敏卡尔曼滤波(UKF),粒子滤波(PF)，高斯混和近似(GMA)等非线性滤波算法。滤波算法中EKF需要线性化且误差较大；UKF不适用于解决状态高维问题；PF算法依赖大量粒子数和重采样方法，计算量大；高斯混合近似方法可逼近任意非线性关系，适用于再入目标的状态强非线性问题，需要解决高斯个数的自适应问题。同时，对于机动目标跟踪问题，多模型是目前最有效方法之一。

【发明内容】

本发明的目的是提供一种基于高斯混合近似的高超声速再入飞行器跟踪方法，以提高机动目标跟踪精度并保持合适计算量，为高超声速机动再入目标跟踪应用提供参考价值。

本发明采用以下技术方案：一种基于高斯混合近似的高超声速再入飞行器跟踪方法，其特征在于，具体包括以下步骤：

步骤一、根据地球圆球模型和拟合大气模型，结合高超声速再入过程中目标的受力情况，在东北上惯性参考坐标系下建立高超声速机动再入飞行器的运动跟踪模型；

步骤二、在东北上惯性参考坐标系下，根据步骤一中的运动跟踪模型，分别建立目标跟踪系统的雷达观测模型和目标先验模型集；

步骤三、根据步骤二中的先验模型集和雷达观测模型，并结合雷达回波数据,获得目标状态平滑估计值和协方差矩阵平滑估计值。

进一步地，步骤一具体方法为：

步骤1.1、根据地球圆球模型和拟合大气模型，得出目标的加速度模型为：

其中，a_x,a_y,a_z分别为目标在东北上惯性参考坐标系中x轴、y轴、z轴的加速度；h为目标在东北上惯性参考坐标系中的高度，ρ(h)为拟合大气模型，v为目标在东北上惯性参考坐标系中的速度，

为弹体坐标系到东北上惯性参考坐标系的转移矩阵，α_v,α_t,α_c分别为气动阻力参数、气动转弯力参数、气动爬升力参数，μ为地球引力常量，x、y、z分别为目标在东北上惯性参考坐标系中的x轴、y轴、z轴的坐标值，r为地心与目标的径向距离，r₀为地球半径，H_R为雷达高度，ω_e为地球自转角速度，B为雷达位置的纬度；

步骤1.2、根据步骤1.1中的目标加速度模型，结合气动力参数得出目标的运动跟踪模型：

X_k+1＝F(θ_k+1)X_k+Ga_k(X_k)+w_k，

其中，F(θ_k+1)为k+1时刻的状态转移矩阵，θ_k+1为k+1时刻的机动时间常数，X_k为k时刻的目标状态，G为过程噪声分布矩阵，a_k为k时刻的目标加速度，w_k为k时刻的零均值过程白噪声。

进一步地，步骤二具体过程为：

雷达观测模型具体为：

其中，x_k、y_k、z_k分别为目标在东北上惯性参考坐标系中k时刻的x轴、y轴、z轴的坐标值，Z_k为直到k时刻的所有量测值，v_k为k时刻的零均值白噪声；

目标先验模型集为：

其中，m为先验模型集的个数，m为大于等于2的整数，i、j分别表示第i、个和第j个先验模型。

进一步地，步骤三的具体方法为：

步骤3.1、将混合概率密度函数，一步预测概率密度函数和条件后验概率密度函数分别用N个高斯分量表示：

混合概率密度函数：

一步预测概率密度函数：

条件后验概率密度函数

其中，p表示概率密度函数，θ_k为k时刻的模式，i和j分别代表第i和第j个先验模型，Z^k为

和

为每个高斯分量在k时刻的权值，α和β分别表示第α和第β个分量，且满足以下公式：

和

和

分别为上述三个公式的一二阶矩；

步骤3.2、根据步骤3.1可得高斯混合模型存在下的输入交互公式：

其中，

μ_k,i|j为目标在k时刻从模型i转移到j的条件概率，μ_k,i为目标在k时刻处于模型i的概率，将所有高斯分量按权值大小排序，按照权值大小，取前m-1个高斯分量，且合并后面的N_t-(m-1)个高斯分量得到一个新的高斯分量；

合并后的第m个高斯分量的状态和协方差表达式分别为：

其中，

为高斯合并后归一化的权值；

步骤3.3、根据步骤3.2及卡尔曼滤波原理，得到状态的一步预测估计为：

计算上述公式的一二阶矩，得到状态预测和协方差预测分别为：

状态预测值：

协方差预测值：

根据k+1时刻的量测，更新k+1时刻的状态估计和协方差估计；

则其后验概率密度函数为下式：

其中，第j个模型的第β个分量的量测预测值为：

第j个模型的第β个分量的量测预测协方差为：

第j个模型的第β个分量的状态与量测间的协方差为：

权值：

第j个模型的第β个分量的状态估计值：

第j个模型的第β个分量的协方差估计值：

第j个模型的第β个分量的增益为：

同时，可得到k+1时刻第j个模型的似然值如下：

步骤3.4、利用区间[k,k+M]的估计结果求取k时刻的状态平滑估计值，具体为：P_GMA(x_k|θ_k＝j,Z^M)，

其中，M为时间窗长；

则从k+1时刻推导出k时刻的状态平滑结果：

其中：

第j个模型的第β个分量的权值平滑结果为：

第j个模型的第β个分量的平滑增益为：

第j个模型的第β个分量的状态平滑结果为：

第j个模型的第β个分量的协方差平滑结果为：

根据步骤3.3中获得的k+1时刻的状态估计值为：

其中，权值

根据上述公式可到区间长度为M中的状态前向平滑值P_GMA(x_k|θ_k＝j,Z^M)：

结合步骤3.2可得到合并后的估计值：

其中，权值为：

计算上式的一二阶矩得到，第j个模型的状态平滑估计值和协方差平滑估计值：

第j个模型的状态平滑估计值：

第j个模型的协方差平滑估计值：

步骤3.5、根据交互式多模型的输出交互方法，得出最终目标状态平滑估计值和协方差平滑估计值：

最终状态平滑估计值为：

最终协方差平滑估计值为：

本发明的有益效果是：在飞行器气动模型的基础上，采用机动时间常数表征目标的机动特性，对高超声速飞行器运动进行建模，形成目标模型库，相对平坦地球模型建立的气动模型(MaRV)和常规运动模型如加速模型，辛格模型，当前统计模型等，可更好的描述再入飞行器的物理特性，提高建模精度；引入高斯混合近似思想，设计多模型下的高斯混合平滑滤波器，采用高斯分量自适应合并策略，同时合理选择平滑窗长，可避免粒子滤波存在的计算大和粒子退化问题，计算过程简单，易于实现，对于突然机动有着很好的估计精度，提高了高超声速再入飞行器的快速高精度跟踪能力。

【附图说明】

图1为本发明一种基于高斯混合近似的高超声速再入飞行器跟踪方法的流程图；

图2为本发明中高超声速再入飞行器CAV-H在再入过程中的受力分析示意图；

图3为本发明中三维空间雷达传感器测量系统示意图，X,Y,Z分别为惯性参考坐标系中x,y,z轴的坐标，α为俯仰角，β为方位角，r为目标距离雷达的距离，T为采样周期，V为目标飞行速度，O为惯性参考坐标系的中心；

图4为本发明中东北上惯性参考坐标系下飞行器的运动轨迹图；

图5为本发明中东北上惯性参考坐标系下的X轴方向飞行器的真实速度曲线图；

图6为本发明中东北上惯性参考坐标系下的Y轴方向飞行器的真实速度曲线图；

图7为本发明中东北上惯性参考坐标系下的Z轴方向飞行器的真实速度曲线图；

图8为本发明中东北上惯性参考坐标系下的X轴方向飞行器的真实加速度曲线图；

图9为本发明中东北上惯性参考坐标系下的Y轴方向飞行器的真实加速度曲线图；

图10为本发明中东北上惯性参考坐标系下的Z轴方向飞行器的真实加速度曲线图；

图11为本发明中飞行器在X轴方向本发明与MaRV相比跟踪结果的均方根误差示意图；

图12为本发明中飞行器在Y轴方向本发明与MaRV相比跟踪结果的均方根误差示意图；

图13为本发明中飞行器在Z轴方向本发明与MaRV相比跟踪结果的均方根误差示意图；

图14为本发明中飞行器在X轴方向本发明与粒子滤波和非平滑算法跟踪结果的均方根误差示意图；

图15为本发明中飞行器在Y轴方向本发明与粒子滤波和非平滑算法跟踪结果的均方根误差示意图；

图16为本发明中飞行器在Z轴方向本发明与粒子滤波和非平滑算法跟踪结果的均方根误差示意图。

【具体实施方式】

下面结合附图和具体实施方式对本发明进行详细说明。

本发明公开了一种基于高斯混合近似的高超声速再入飞行器跟踪方法，如图1所示，以高超声速再入飞行器的气动力学运动模型为基础，分析不同地球模型和力对目标运动建模的影响程度，建立了利用机动时间常数描述气动参数的运动跟踪模型，并采取高斯混合近似平滑滤波器对目标(即飞行器)的状态进行估计。

具体包括以下步骤：

步骤一、根据地球圆球模型和拟合大气模型，结合高超声速再入过程中目标的受力情况，受力包括气动阻力、转弯力、爬升力、地球引力及哥氏力，结合地球圆球模型和拟合大气模型在东北上(ENU)惯性参考坐标系下建立高超声速机动再入飞行器的运动跟踪模型(HGRV)。

其中，其中东北上惯性参考坐标系原点为雷达中心，X轴和Y轴是地球参考椭球的切线，分别指向东和北，Z轴垂直于当地水平面向上。

由于地球模型从椭球到圆球的简化误差很小，从圆球到平坦模型的简化误差显著，且在高空(60km)时地球自转引起的哥氏力与气动力大小相当，低空(30km)时指数大气模型会带来一定气动力偏差。本发明考虑目标受空气阻力和地球引力作用，采用考虑哥氏力的圆球模型和拟合大气模型用于高超声飞行器再入段跟踪问题。由于机动时间常数可表示目标机动的能力，故气动参数演化过程表示为气动时间常数的函数，过程噪声也与机动时间常数有关，具有一定的物理意义；同时避免使用椭球模型带来的计算复杂度，既保证了建模精度又便于快速计算。

步骤1.1、通过分析，如图2所示，地球模型从椭球到圆球的简化误差很小，从圆球到平坦模型的简化误差显著，且在高空(60km)时地球自转引起的哥氏力与气动力大小相当，低空(30km)时指数大气模型会带来一定气动力偏差。

因此，本发明考虑目标受空气阻力和地球引力作用，采用考虑哥氏力的圆球模型和拟合大气模型用于高超声飞行器再入段跟踪问题。

拟合大气模型为：

ρ(h)＝ρ₀e^-ch，

其中，ρ₀为海平面大气密度，h为目标高度，c为大气密度计算公式的拟合系数。

东北上惯性参考坐标系下的(目标)飞行器速度为：

其中，v_x,v_y,v_z分别为目标在东北上惯性参考坐标系中三个轴的速度分量，v为目标的速度。

建立弹体坐标系，弹体坐标系的原点0为导弹质心，OX轴为弹体外壳对称轴，指向导弹头部；OY轴位于导弹纵对称面内，OY轴垂直于OX轴；OZ轴垂直于导弹纵对称面，在发射瞬时顺着发射方向看过去，OZ指向右方。OZ、OX、OY组成右手直角坐标系，且OZ、OX、OY与弹体的相对方向关系在导弹飞行过程中保持不变，这主要是因为现代飞行器很可能进行翻转运动，此时OZ可能指向左方，但三轴方向与弹体的相对关系在发射之初就已确定，并不再改变。在弹体坐标系上，建立飞行器的气动参数模型，气动阻力参数α_v沿OX轴反方向，气动转弯力α_t参数沿着OZ轴方向，气动爬升力参数α_c沿着OY轴方向。

弹体坐标系(BF)到东北上坐标系(ENU)的转移矩阵为：

其中，

则可得到加速度模型为：

其中，a_x,a_y,a_z分别为目标在东北上惯性参考坐标系中x轴、y轴、z轴的加速度；h为目标在东北上惯性参考坐标系中的高度，ρ(h)为拟合大气模型，v为目标在东北上惯性参考坐标系中的速度，

为弹体坐标系到东北上惯性参考坐标系的转移矩阵，α_v,α_t,α_c分别为气动阻力参数、气动转弯力参数、气动爬升力参数，μ为地球引力常量，x、y、z分别为目标在东北上惯性参考坐标系中的x轴、y轴、z轴的坐标值，r为地心与目标的径向距离，r₀为地球半径，H_R为雷达高度，ω_e为地球自转角速度，B为雷达位置的纬度。

步骤1.2、气动力参数假设为机动时间常数θ的函数，其演化过程满足指数衰减形式的一阶马尔科夫过程，T为雷达采样周期，气动力参数函数离散化的公式为：

假定气动力参数方差为σ_α，则气动力参数噪声

的方差为：

Q_α＝σ_α ²(1-exp(-2T/θ))≈2σ_α ²(T/θ)，

则目标的待估计状态为：

X＝[x y z v_x v_y v_z α_v α_t α_c]，

根据步骤1.1中的目标加速度模型，结合气动力参数得出目标的运动跟踪模型：

X_k+1＝F(θ_k+1)X_k+Ga_k(X_k)+w_k，

其中，F(θ_k+1)为k+1时刻的状态转移矩阵，θ_k+1为k+1时刻的机动时间常数，X_k为k时刻的目标状态，G为过程噪声分布矩阵，a_k为k时刻的目标加速度，w_k为k时刻的零均值过程白噪声。

状态转移矩阵F(θ_k+1)为：

过程噪声分布矩阵为：

其中，I_3×3,0_3×3分别为3*3单位矩阵和全零矩阵。

步骤二、在东北上惯性参考坐标系下，根据步骤一中的运动跟踪模型，分别建立目标跟踪系统的雷达观测模型和目标先验模型集。

雷达传感器为地基X波段雷达，可用于高超声速飞行器等再入目标的末端拦截，如图3所示，探测系统可获得目标的径向距、方位角、俯仰角信息。

建立目标的先验模型集，主要由不同的机动时间常数表示，机动时间常数越大，目标飞行越平稳，反之则运动越波动，也就是机动性能强。

雷达传感器可得到目标的径向距、方位角、俯仰角，假定量测噪声为零均值过程白噪声，量测模型为：

其中，x_k、y_k、z_k分别为目标在东北上惯性参考坐标系中k时刻的x轴、y轴、z轴的坐标值，v_k为k时刻的零均值过程白噪声。

根据气动力参数离散化的公式和适用于跟踪的运动模型，假定机动时间常数取值共m个，即θ₁,θ₂,Lθ_m，每个模型的初始概率相同均为1/m，m个模型之间的转移概率满足主对角线占优原则，即目标先验模型集为：

其中，先验模型集由不同的机动时间常数表示，m为机动时间常数的个数，m大于等于2的整数，i、j分别表示第i个和第j个模型。

步骤三、根据步骤二中的先验模型集和雷达观测模型，并结合雷达回波数据,采用高斯混合近似滤波器进行计算，并设定平滑窗长，在多模型框架下获得目标位置和速度及气动参数的状态平滑估计值和协方差矩阵平滑估计值。

根据高超声速再入飞行器的运动模型和雷达传感器的量测模型，并采用前两拍数据进行差分得到状态初始值，假定每一时刻的状态可由m个模型按照不同权值表征，同时每个模型下由N个高斯分量表征状态。通过计算每个模型下的所有高斯分量的状态估计结果，并与多个模型相乘，得到N*m个状态估计值，通过权值排序进行高斯项合并，得到目标的状态融合估计值和协方差矩阵融合结果。

经过一定时间窗长M后，得到k+M时刻的估计值，与其他多模型和高斯混合模型不同的是，本发明需要进行状态和协方差的平滑计算。设定固定时间窗长M，由k+M时刻的结果回溯至k时刻，得到高斯混合滤波器的状态平滑结果。

在高超声速再入飞行器的跟踪过程中，为保证估计的精度和实时性，需要进行多次仿真，得到最佳的高斯分量个数N和平滑时间窗长M。

其具体方法为：

步骤3.1、考虑到多模型框架下存在模型切换和模型非线性的耦合问题，将混合概率密度函数，一步预测概率密度函数和条件后验概率密度函数分别用N个高斯分量表示：

混合概率密度函数：

一步预测概率密度函数：

条件后验概率密度函数

其中，p表示概率密度函数，θ_k为k时刻的模式，i和j分别代表第i和第j个模型，Z^k为

和

为每个高斯分量在k时刻的的权值，α和β分别表示第α和第β个分量，且满足以下公式：

和

和

分别为上述三个公式的一二阶矩；

步骤3.2、与交互式多模型的输入交互类似，本发明推导出高斯混合模型存在下的输入交互公式，具体如下：

其中，

μ_k,i|j为目标在k时刻从模型i转移到j的条件概率，μ_k,i为目标在k时刻处于模型i的概率。从公式可以看出，由于模型交互作用，上述公式中的高斯分量个数变成N_t＝N×m，为了减少计算量，将所有高斯分量按权值大小排序，按照权值大小，取前m-1个高斯分量，且合并后面的N_t-(m-1)个高斯分量得到一个新的高斯分量，最终仍保持m个高斯分量，实现了高斯分解与合并自适应；

合并后的第m个高斯分量的状态和协方差表达式分别为：

其中，

为高斯合并后归一化的权值；

步骤3.3、经过输入交互后，根据步骤3.2及卡尔曼滤波原理，得到状态的一步预测估计为：

计算上述公式的一二阶矩，得到状态预测和协方差预测分别为：

状态预测值：

协方差预测值：

根据k+1时刻的量测，更新k+1时刻的状态估计和协方差估计；

则其后验概率密度函数为下式：

其中，第j个模型的第β个分量的量测预测值为：

第j个模型的第β个分量的量测预测协方差为：

第j个模型的第β个分量的状态与量测间的协方差为：

权值：

第j个模型的第β个分量的状态估计值：

第j个模型的第β个分量的协方差估计值：

第j个模型的第β个分量的增益为：

同时，可得到k+1时刻第j个模型的似然值如下：

步骤3.4、利用区间[k,k+M]的估计结果求取k时刻的状态平滑估计值，具体为：P_GMA(x_k|θ_k＝j,Z^M)，

其中，M为时间窗长；

则从k+1时刻推导出k时刻的状态平滑结果：

其中：

第j个模型的第β个分量的权值平滑结果为：

第j个模型的第β个分量的平滑增益为：

第j个模型的第β个分量的状态平滑结果为：

第j个模型的第β个分量的协方差平滑结果为：

根据步骤3.3中获得的k+1时刻的状态估计值为：

其中，权值满足和为1的条件，即

根据上述公式可到区间长度为M中的状态前向平滑值P_GMA(x_k|θ_k＝j,Z^M)：

结合步骤3.2中的高斯分解与合并自适应策略，可得到合并后的估计值：

其中，权值为：

计算上式的一二阶矩得到，第j个模型的状态平滑估计值和协方差平滑估计值：

第j个模型的状态平滑估计值：

第j个模型的协方差平滑估计值：

步骤3.5、根据交互式多模型的输出交互方法，得出最终目标状态平滑估计值和协方差平滑估计值：

最终状态平滑估计值为：

最终协方差平滑估计值为：

仿真场景如图4所示，雷达采样周期T为1秒，雷达坐在经度为17度，纬度为139度，雷达测量的标准误差分别为：径向距100米，方位角0.01弧度，俯仰角0.01弧度；海平面大气密度ρ₀为1.2258千克/立方米，拟合系数c取值为-1.3785*10^-4，地球半径为6.34*10⁶公里，地球引力常量μ取3.986004418*10^14，地球自转角速度ω_e取7.292*10^-5弧度/秒。

本发明参考马丁公司设计的高超声速飞行器CAV-H的模型参数，飞行器质量为907.2千克，初始位置经纬高设定为135度,35度，80公里，初始攻角设为10度，初始速度5马赫，航向角135度，航迹角2度；

先验模型集个数共取m＝3个，分别为：θ₁＝5s,θ₂＝5*10³,θ₃＝10⁴s，分别表示目标处于剧烈运动、缓和运动、平稳运动，各模型的初始概率为1/3，飞行过程中，目标可从任一状态跳转到其他状态。

设仿真全过程共320秒，为了说明再入飞行器的机动性能，分别给出了速度和加速度曲线图；

如图5、图6、图7所示，给出了飞行器在三个轴方向上的真实速度曲线图，坐标原点设为(0,0,0)。

如图8、图9、图10所示，给出了飞行器在三个轴方向上的真实加速度曲线图，可以看出目标在100-230秒之间机动较大，其他段为非机动段。

为了对比分析说明本发明中的建模精确性和算法的有效性，进行如下两个实验：

仿真试验中，对比算法除了算法核心参数不同，其他参数均为相同，并进行100次蒙特卡洛仿真。

任一参数的均方根误差计算公式为：

其中，n为蒙特卡洛仿真次数，

为k时刻的真实参数信息，

为k时刻的滤波器估计的参数值；

实验一：将本发明中的机动时间常数表示的气动模型(HGRV)与传统机动再入模型(MaRV)进行对比，其中MaRV模型采用平坦地球模型，指数大气模型且不考虑地球自转的影响；

如图11、图12、图13所示，给出了飞行器在三个轴方向本发明提出的模型与MaRV模型相比跟踪结果的均方根误差示意图。可以看到，在X轴方向上，本发明中的HGRV模型与MaRV模型估计精度相当，这是由于图5所示，目标在X轴方向上机动过载变化不大，两种模型均可以较好的描述；在Y轴和Z轴方向上，本发明的模型相比MaRV模型精度有了明显的提升，同时也可以从图6-7看出，这两个方向上的加速度变化剧烈，机动性能更强；所以可以得出结论：本发明中的模型更适用于机动较强的飞行器，对于机动较弱的情况与传统MaRV模型精度相当。

实验二：为了说明本发明提出的高斯混合平滑滤波器(IMM-GMAS)的效果，与粒子滤波算法(PF)和不带平滑的高斯混合滤波器(IMM-GMA)进行对比,其中粒子滤波的粒子数设为P＝1000；

如图14、图15、图16所示，给出了飞行器在三个轴方向本发明与粒子滤波和非平滑算法跟踪结果的均方根误差示意图；可以看出，首先将不带平滑的IMM-GMA算法与粒子滤波算法进行对比，本发明中的高斯混合(GMA)滤波器比粒子滤波估计精度更高，且降低了峰值误差，表明本发明算法应对机动情况比粒子滤波要好，进一步说明本发明中基于多模型的高斯混合滤波器在处理目标机动时具有很好的性能；

其次，将带平滑的高斯混合滤波器(IMM-GMAS)与不带平滑的高斯混合滤波器(IMM-GMA)比较，设置平滑窗长为2，从图14-16可以看出，IMM-GMAS的估计误差整体上小于IMM-GMA，说明本发明的IMM-GMAS的平滑策略可以进一步提高估计精度，这在本发明首次提出多模型框架下的高斯混合平滑滤波器，具有很好的估计性能。

在计算复杂度上，各种算法的100次蒙特卡罗仿真计算时间归一化后如表1所示：

表1

算法	计算时间
		IMM-GMAS(HGRV)	2.13
IMM-GMAS(MaRV)	2.07
		IMMGMA(HGRV)	1
PF(P＝1000,HGRV)	2.69

参照表1，可以看出本发明的IMM-GMAS算法计算量适中，比粒子滤波要低，由于平滑策略的存在计量为不带平滑的IMM-GMAS两倍左右，即使如此，本发明算法的每拍计算时间均低于采样周期T＝1，可保证实时估计，符合实际情况。使用本发明提出的一种基于高斯混合近似的高超声速再入飞行器跟踪方法，在跟踪CAV-H飞行器再入段时计算量适中，而且估计精度较高，对于高超声速再入飞行器的实时稳定跟踪是个很好的选择，具有工程实用价值。

Claims (3)

1.一种基于高斯混合近似的高超声速再入飞行器跟踪方法，其特征在于，具体包括以下步骤：

步骤一、根据地球圆球模型和拟合大气模型，结合高超声速再入过程中目标的受力情况，在东北上惯性参考坐标系下建立高超声速机动再入飞行器的运动跟踪模型；

步骤二、在东北上惯性参考坐标系下，根据步骤一中的所述运动跟踪模型，分别建立目标跟踪系统的雷达观测模型和目标先验模型集；

步骤三、根据步骤二中的所述先验模型集和雷达观测模型，并结合雷达回波数据,获得目标状态平滑估计值和协方差矩阵平滑估计值；

步骤三的具体方法为：

步骤3.1、将混合概率密度函数，一步预测概率密度函数和条件后验概率密度函数分别用N个高斯分量表示：

混合概率密度函数：

一步预测概率密度函数：

条件后验概率密度函数

其中，p表示概率密度函数，θ_k为k时刻的模式，i和j分别代表第i和第j个先验模型，Z^k为

和

为每个高斯分量在k时刻的权值，α和β分别表示第α和第β个分量，且满足以下公式：

和

和

分别为上述三个公式的一二阶矩；

步骤3.2、根据步骤3.1可得高斯混合模型存在下的输入交互公式：

其中，

μ_k,i|j为目标在k时刻从模型i转移到j的条件概率，μ_k,i为目标在k时刻处于模型i的概率，将所有高斯分量按权值大小排序，按照权值大小，取前m-1个高斯分量，且合并后面的N_t-(m-1)个高斯分量得到一个新的高斯分量；

合并后的第m个高斯分量的状态和协方差表达式分别为：

其中，

为高斯合并后归一化的权值；

步骤3.3、根据步骤3.2及卡尔曼滤波原理，得到状态的一步预测估计为：

计算上述公式的一二阶矩，得到状态预测和协方差预测分别为：

状态预测值：

协方差预测值：

根据k+1时刻的量测，更新k+1时刻的状态估计和协方差估计；

则其后验概率密度函数为下式：

其中，第j个模型的第β个分量的量测预测值为：

第j个模型的第β个分量的量测预测协方差为：

第j个模型的第β个分量的状态与量测间的协方差为：

权值：

第j个模型的第β个分量的状态估计值：

第j个模型的第β个分量的协方差估计值：

第j个模型的第β个分量的增益为：

同时，可得到k+1时刻第j个模型的似然值如下：

步骤3.4、利用区间[k,k+M]的估计结果求取k时刻的状态平滑估计值，具体为：P_GMA(x_k|θ_k＝j,Z^M)，

其中，M为时间窗长；

则从k+1时刻推导出k时刻的状态平滑结果：

其中：

第j个模型的第β个分量的权值平滑结果为：

第j个模型的第β个分量的平滑增益为：

第j个模型的第β个分量的状态平滑结果为：

第j个模型的第β个分量的协方差平滑结果为：

根据步骤3.3中获得的k+1时刻的状态估计值为：

其中，权值

根据上述公式可到区间长度为M中的状态前向平滑值P_GMA(x_k|θ_k＝j,Z^M)：

结合步骤3.2可得到合并后的估计值：

其中，权值为：

计算上式的一二阶矩得到，第j个模型的状态平滑估计值和协方差平滑估计值：

第j个模型的状态平滑估计值：

第j个模型的协方差平滑估计值：

步骤3.5、根据交互式多模型的输出交互方法，得出最终目标状态平滑估计值和协方差平滑估计值：

最终状态平滑估计值为：

最终协方差平滑估计值为：

2.如权利要求1所述的基于高斯混合近似的高超声速再入飞行器跟踪方法，其特征在于，步骤一具体方法为：

步骤1.1、根据地球圆球模型和拟合大气模型，得出目标的加速度模型为：

其中，a_x,a_y,a_z分别为目标在东北上惯性参考坐标系中x轴、y轴、z轴的加速度；h为目标在东北上惯性参考坐标系中的高度，ρ(h)为拟合大气模型，v为目标在东北上惯性参考坐标系中的速度，

为弹体坐标系到东北上惯性参考坐标系的转移矩阵，α_v,α_t,α_c分别为气动阻力参数、气动转弯力参数、气动爬升力参数，μ为地球引力常量，x、y、z分别为目标在东北上惯性参考坐标系中的x轴、y轴、z轴的坐标值，r为地心与目标的径向距离，r₀为地球半径，H_R为雷达高度，ω_e为地球自转角速度，B为雷达位置的纬度；

步骤1.2、根据步骤1.1中的目标加速度模型，结合气动力参数得出目标的运动跟踪模型：

X_k+1＝F(θ_k+1)X_k+Ga_k(X_k)+w_k，

其中，F(θ_k+1)为k+1时刻的状态转移矩阵，θ_k+1为k+1时刻的机动时间常数，X_k为k时刻的目标状态，G为过程噪声分布矩阵，a_k为k时刻的目标加速度，w_k为k时刻的零均值过程白噪声。

3.如权利要求1所述的基于高斯混合近似的高超声速再入飞行器跟踪方法，其特征在于，步骤二具体过程为：

所述雷达观测模型具体为：

其中，x_k、y_k、z_k分别为目标在东北上惯性参考坐标系中k时刻的x轴、y轴、z轴的坐标值，Z_k为直到k时刻的所有量测值，v_k为k时刻的零均值白噪声；

所述目标先验模型集为：

其中，m为先验模型集的个数，m为大于等于2的整数，i、j分别表示第i、个和第j个先验模型。

Images