基于PHD滤波的临近空间高超声速目标跟踪方法

摘要

本发明公开了一种基于PHD滤波的临近空间高超声速目标跟踪方法，包括步骤如下：使用初始时刻的量测进行对滤波器初始化，将前四个时刻的雷达量测值从雷达站球坐标系转到雷达站ENU直角坐标系，采用两点差分法，得到初始目标强度函数；将初始目标输入初始化后的滤波器中，按照GM‑PHD滤波器的过程进行计算，得到预测目标集；根据当前时刻量测与由上一时刻预测得到的目标位置之间的马氏距离，将量测划分为存活目标量测与杂波量测；利用当前时刻的量测值得到滤波器当前的更新值，完成对存活目标的更新；对步骤S4的更新公式中的高斯项进行剪枝和合并；计算修剪合并后的强度函数、提取目标状态，完成对目标状态估计。本发明实现了强杂波下对未知数目的临近空间飞行器进行精确跟踪。

说明书

技术领域

本发明涉及信号处理技术领域，更具体的，涉及一种基于PHD滤波的临近空间高超声速目标跟踪方法。

背景技术

近年来，随着航天技术的发展，基于助推-滑翔弹道概念的高超声速飞行器逐渐成为国内外研究的热点。这类高超声速飞行器可在大气中滑翔飞行，进行大范围快速机动，具有射程远，速度快，机动性强等优点。对于单目标高超声速飞行器跟踪，以扩展卡尔曼滤波(Extended Kalman Filter,EKF)为代表的非线性滤波算法发挥了重要的作用对高机动目标取得了良好的跟踪效果，实现了对弹道式再入目标的高精度跟踪，如文献[Jawahar A,Koteswara Rao S.Modified Polar Extended Kalman Filter(MP-EKF)for Bearings-Only Target Tracking[J].Indian Journal of ence and Technology,2016,9(26):1-5]、[Belik B V,Belov S G.Using of Extended Kalman Filter for Mobile TargetTracking in the Passive Air Based Radar System☆[J].Procedia ComputerScience,2017,103:280-286]。

在强杂波环境中，杂波量测与目标量测难以区分，初始目标数目未知，同时，针对高超声速飞行器携带多个滑翔弹头，利用诱饵弹增强其突防能力的问题，有必要讨论高超声速多目标跟踪。传统的多目标跟踪以数据关联为核心，代表算法有多级假设检验方法(Multiple Hypothesis Tracking,MHT)和联合概率数据关联(Joint Probabilistic DataAssociation,JPDA)，二者在目标数目和杂波数增加时，算法复杂度呈指数增长，即出现“组合爆炸”。基于随机有限集(random finite sets,RFS)理论的多目标跟踪算法避免了数据关联，其中最重要的工作是Mahler等人提出的概率假设密度(Probability HypothesisDensity,PHD)滤波器[Mahler R.P.S..Multitarget Bayes filtering via first-ordermultitarget moments[J].IEEE Transactions on Aerospace and Electronic Systems,2003,39(4):1152-1178]。其主要实现方法有高斯混合PHD(Gaussian mixture-PHD,GM-PHD)滤波器和序贯蒙特卡洛PHD滤波器(Sequential Monte Carlo-PHD,SMC-PHD)。

目前，使用PHD滤波器进行高超声速多目标跟踪存在以下两个问题：

(1)传统的PHD滤波器通常假定已知初始目标的起始位置信息和新生目标强度函数，在实际应用中，不能直接得到上述信息。

(2)目标量测附近的杂波干扰往往会导致滤波器性能下降，出现错跟，目标数目估计过大的情况。

发明内容

本发明为了解决传统GMPHD算法需要知道初始目标的起始位置信息和新生目标强度函数，且杂波干扰导致滤波器性能下降的问题，提供了一种基于PHD滤波的临近空间高超声速目标跟踪方法，其实现了强杂波下对未知数目的临近空间飞行器进行精确跟踪。

为实现上述本发明目的，采用的技术方案如下：一种基于PHD滤波的临近空间高超声速目标跟踪方法，所述的方法包括步骤如下：

S1：使用初始时刻的量测进行对滤波器初始化，将前四个时刻的雷达量测值从雷达站球坐标系转到雷达站ENU直角坐标系，采用两点差分法，得到初始目标强度函数；

S2：将初始目标输入初始化后的滤波器中，预测步骤按照GM-PHD滤波器的过程进行计算，得到预测目标集；

S3：根据当前时刻量测与由上一时刻预测得到的目标位置之间的马氏距离，将量测划分为存活目标量测与杂波量测；

S4：利用当前时刻的量测值得到滤波器当前的更新值，完成对存活目标的更新；

S5：设置修剪门限、合并门限对步骤S4的更新公式中的高斯项进行剪枝和合并；

S6：计算修剪合并后的强度函数、提取目标状态，完成对目标状态估计。

优选地，步骤S1，初始目标强度函数的表达式如下：

其中，N

式中，x

设所有初始高斯项权值

进一步地，将初始目标输入初始化后的滤波器进行预测，其中初始目标的状态转移的方程表示如下：

X(k+1)＝F(k)X(k)+W(k)

其中，F(k)为状态转移矩阵，有

其中，F

式中：

p1＝(2-2αT+α

q1＝(αT-1+e

r1＝(1-e

s1＝e

W(k)表示均值为0，协方差为Q(k)的高斯白噪声序列；

预测步骤具体计算公式如下：

v

其中

其中，y

再进一步地，步骤S3，划分量测具体如下：

从时刻2开始，对k时刻的量测进行划分，记：

其中，s

其中，R

则落在门限值内的量测集合作为存活目标量测，

落在门限值外的量测集合作为杂波量测，

其中T表示门限值。

再进一步地，所述的门限值由下式确定，若P

T＝-2ln(1-P

再进一步地，步骤S4，使用量测集合

其中，

其中，H

设V为整个观测区域面积，λ为杂波平均数，V

则杂波强度

再进一步地，步骤S5，所述的剪枝和合并，具体如下：

令l＝0，

当集合I非空，重复以下过程：

l＝l+1，

I＝I\L，直到

其中，U_merg表示合并门限。

再进一步地，步骤S5，修剪合并后的强度函数为：

再进一步地，步骤S5，根据下式提取目标状态：

本发明的有益效果如下：

本发明能够实现对多个高机动目标的跟踪，尤其在目标飞行轨迹高度非线性的情况表现良好，同时还可以克服传统GMPHD算法需要初始目标信息的缺陷，并且能够在强杂波环境下对未知数目临近空间飞行器进行跟踪。

附图说明

图1是本实施例所述的临近空间高超声速目标跟踪方法的步骤流程图。

图2是本实施例所述的助推-滑翔弹道的示意图。

图3是通过本实施例所述的方法的跟踪结果的示意图。

图4是通过本实施例所述的方法得到的前100秒的目标数目估计的效果图，图中，A表示传统EK-GMPHD方法，B表示本实施例所述的目标跟踪方法。

图5是采用Wasserstein距离评价目标状态估计的效果图。

具体实施方式

下面结合附图和具体实施方式对本发明做详细描述。

实施例1

如图1所示，一种基于PHD滤波的临近空间高超声速目标跟踪方法，所述的方法包括步骤如下：

步骤S1：使用初始时刻的量测进行对滤波器初始化，将前四个时刻的雷达量测值从雷达站球坐标系转到雷达站ENU直角坐标系，采用两点差分法，得到初始目标强度函数；

所述的初始目标强度函数的表达式如下：

其中，N

式中，x

设所有初始高斯项权值

步骤S2：将初始目标输入初始化后的滤波器中，预测步骤按照GM-PHD滤波器的过程进行计算，得到预测目标集；

其中初始目标的状态转移的方程表示如下：

X(k+1)＝F(k)X(k)+W(k)

其中，F(k)为状态转移矩阵，有

其中，F

式中：

p1＝(2-2αT+α

q1＝(αT-1+e

r1＝(1-e

s1＝e

W(k)表示均值为0，协方差为Q(k)的高斯白噪声序列。

预测步骤具体计算公式如下：

v

其中

其中，y

步骤S3：根据当前时刻量测与由上一时刻预测得到的目标位置之间的马氏距离，将量测划分为存活目标量测与杂波量测；

从时刻2开始，对k时刻的量测进行划分，记：

其中，s

其中，R

则落在门限值内的量测集合作为存活目标量测，

落在门限值外的量测集合作为杂波量测，

其中T表示门限值。

所述的门限值由下式确定，若P

T＝-2ln(1-P

S4：利用当前时刻的量测值得到滤波器当前的更新值，完成对存活目标的更新；

使用量测集合

其中，

其中，H

上述方程组中杂波强度K

设V为整个观测区域面积，λ为杂波平均数，V

则杂波强度

步骤S5：设置修剪门限、合并门限对步骤S4的更新公式中的高斯项进行剪枝和合并；

令l＝0，

当集合I非空，重复以下过程：

l＝l+1，

I＝I\L，直到

其中，U_merg表示合并门限。

S6：计算修剪合并后的强度函数、提取目标状态，完成对目标状态估计。

修剪合并后的强度函数为：

根据下式提取目标状态：

为了进一步对本实施例所述的临近空间高超声速目标跟踪方法的技术效果，具体举例说明如下：

在雷达站ENU三维直角坐标系下跟踪两个高机动目标，扫描周期T＝1s，目标初始状态如表1：

表1目标初始状态

临近空间超高声速飞行器的空气动力模型是以下一组非线性常微分方程组，

其中，u为地球重力常数，u＝3.986 005×10

ρ(h)＝ρ

其中ρ

使用四阶龙哥库塔法求解微分方程组的数值解，得到仿真的弹道，微分方程组初值如表1，得到助推-滑翔弹道如图2所示。

设置机动频率α＝0.5，目标在x，y，z方向上的加加速度方差

(1)按照上述的仿真参数设置对变量进行初始化，以及按照本实施例所述的方法的步骤S1对滤波器初始化；

(2)按照本实施例所述的方法的步骤S2得到预测目标集；

(3)按照本实施例所述的方法的步骤S3对得到的量测进行划分；

(4)按照本实施例所述的方法的步骤S4利用当前时刻的量测值得到滤波器当前的更新值，其中量测方程h(·)为

其中，Z

其中，

(5)按照本实施例所述的方法的步骤S5对高斯项进行剪枝和合并；

(6)按照本实施例所述的方法的步骤S6提取目标状态；

(7)循环步骤S2-S6直到结束。

最终的跟踪结果见附图3，最终得到的前100秒的目标数目估计如附图4所示，其中A代表传统EK-GMPHD方法，B代表本实施例所述的目标跟踪方法。

可以看出，传统EK-GMPHD方法会造成目标数目过估，采用目标估计绝对误差的时间平均作为目标数目估计的评价指标，即：

其中，

采用Wasserstein距离评价目标状态估计精度，Wasserstein距离是度量两组点集相似程度的指标，其数值越小，目标估计精度越高，实验结果见附图5，其中A代表传统EK-GMPHD方法，B代表本实施例所述的目标跟踪方法。

可以看出，使用本实施例所述的目标跟踪方法目标状态估计精度要优于传统方法。传统EK-GMPHD方法的Wasserstein距离的时间平均值为377 04，本实施例所述的目标跟踪方法的Wasserstein距离的时间平均值为193 40，仅为传统EK-GMPHD方法的51.29％。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。